36
36
void local_init(void)
38
int i, j, k, ij, stat, a, b, r, l, I, L;
39
int nmo, nso, nocc, nocc_all, nvir, noei, nirreps, nfzc;
40
double **C; /* AO -> localized MO transformation matrix */
41
double **Ci; /* localized MO -> AO transformation matrix */
42
double **D; /* 1/2 SCF closed-shell density matrix (AO) */
43
double **Rt, **Rt_full; /* Projected, redundant virtual transform (R-tilde) */
44
double **S; /* AO overlap */
45
double **St; /* Projected virtual overlap */
46
double **Xt; /* Projected, non-redundant virtual transform (X-tilde) */
47
double ***V; /* MO -> projected, redundant virtual transform */
48
double **Fmo;/* MO basis Fock matrix */
49
double **F; /* AO basis Fock matrix */
50
double **Ft; /* Projected, redundant virtual Fock matrix */
51
double **Fbar; /* Projected, non-redundant virtual Fock matrix */
52
double ***W; /* Transformation matrix from tilde -> bar for each ij pair*/
53
double *eps_occ; /* occupied orbital energies for local denominators */
54
double **eps_vir; /* virtual orbital energies for local denominators */
56
double *evals, **evecs;
57
double *eps_all; /* All MO energies */
60
int natom, atom, am, offset, nshell, shell_length, next_atom;
61
int row, col, max, m, errcod, cnt;
62
int *rank, *boolean, *ipiv;
63
int *l_length, *aostart, *aostop, *ao2atom;
65
int **domain, *domain_len, **pairdomain, *pairdom_len, *pairdom_nrlen;
67
double *fR, cutoff, *charge, *SR, *Z, tmp, *ss;
69
int print_test, num_entries, entry_len, orbital;
70
int t1_length, t2_length, puream, weak;
74
double **U, **Ui, **R, **WW, **RS;
75
double **evecs_u, **evecs_v, *work, *sigma;
79
38
chkpt_init(PSIO_OPEN_OLD);
81
natom = chkpt_rd_natom();
82
nshell = chkpt_rd_nshell();
83
eps_all = chkpt_rd_evals();
84
stype = chkpt_rd_stype();
85
snuc = chkpt_rd_snuc();
39
local.natom = chkpt_rd_natom();
90
/* C1 symmetry only */
91
nirreps = moinfo.nirreps;
93
fprintf(outfile, "\nError: localization must use C1 symmetry.\n");
94
exit(PSI_RETURN_FAILURE);
98
nmo = moinfo.nmo; /* should be the same as nso */
100
fprintf(outfile, "\nError: NMO != NSO! %d != %d\n", nmo, nso);
101
exit(PSI_RETURN_FAILURE);
104
nocc = moinfo.occpi[0]; /* active doubly occupied orbitals */
105
nfzc = moinfo.frdocc[0]; /* frozen doubly occupied orbitals */
106
nocc_all = nocc + nfzc; /* all doubly occupied orbitals */
107
nvir = moinfo.virtpi[0]; /* active virtual orbitals */
114
/* A couple of scratch arrays */
115
X = block_matrix(nso, nso);
116
Y = block_matrix(nso, nso);
119
Ci = block_matrix(nso, nso);
120
for(i=0; i < nso; i++)
121
for(j=0; j < nso; j++)
124
invert_matrix(C, Ci, nso, outfile);
126
for(i=0; i < nso; i++)
127
for(j=0; j < nso; j++)
131
fprintf(outfile, "\n\tC inverse (Ci):\n");
132
print_mat(Ci, nso, nso, outfile);
135
/* Compute the AO-basis overlap integrals */
137
S = block_matrix(nso,nso);
138
C_DGEMM('t','n',nso,nso,nso,1.0,&(Ci[0][0]),nso,&(Ci[0][0]),nso,
142
/* Get the overlap integrals -- these should be identical to AO S */
143
noei = nso*(nso+1)/2;
144
ss = init_array(noei);
145
stat = iwl_rdone(PSIF_OEI,PSIF_SO_S,ss,noei,0,0,outfile);
146
S = block_matrix(nso,nso);
147
for(i=0,ij=0; i < nso; i++)
148
for(j=0; j <= i; j++,ij++) {
149
S[i][j] = S[j][i] = ss[ij];
154
fprintf(outfile, "\n\tAO Overlap (S)\n");
155
print_mat(S, nso, nso, outfile);
158
/* Build the SCF closed-shell density matrix/2 */
159
D = block_matrix(nso,nso);
160
for(i=0; i < nso; i++)
161
for(j=0; j < nso; j++)
162
for(k=0; k < nocc_all; k++)
163
D[i][j] += C[i][k] * C[j][k];
166
fprintf(outfile, "\n\tAO-basis SCF Density (D):\n");
167
print_mat(D, nso, nso, outfile);
171
/* Compute the length of each AM block */
172
ip_boolean("PUREAM", &(puream), 0);
173
l_length = init_int_array(LIBINT_MAX_AM);
175
for(l=1; l < LIBINT_MAX_AM; l++) {
176
if(puream) l_length[l] = 2 * l + 1;
177
else l_length[l] = l_length[l-1] + l + 1;
180
/* Set up the atom->AO and AO->atom lookups */
181
aostart = init_int_array(natom);
182
aostop = init_int_array(natom);
183
for(i=0,atom=-1,offset=0; i < nshell; i++) {
185
shell_length = l_length[am];
187
if(atom != snuc[i] - 1) {
188
if(atom != -1) aostop[atom] = offset-1;
190
aostart[atom] = offset;
192
offset += shell_length;
194
aostop[atom] = offset-1;
196
ao2atom = init_int_array(nso);
197
for(i=0; i < natom; i++)
198
for(j=aostart[i]; j < aostop[i]; j++) ao2atom[j] = i;
203
/************* Build the orbital domains ************/
205
domain = init_int_matrix(nocc, natom);
206
domain_len = init_int_array(nocc);
207
charge = init_array(natom);
208
rank = init_int_array(natom);
209
boolean = init_int_array(natom);
210
SR = init_array(nso);
212
ipiv = init_int_array(nso);
213
fR = init_array(nocc);
215
for(i=nfzc; i < nocc_all; i++) {
217
/* Compute the contribution of each atom to this orbital's charge/population */
218
for(j=0; j < natom; j++) {
220
for(k=aostart[j]; k <= aostop[j]; k++) {
222
for(l=0; l < nso; l++) tmp += S[k][l] * C[l][i];
228
/* Rank the atomic contributions to the orbital's charge */
229
for(j=0; j < natom; j++) { rank[j] = 0; boolean[j] = 0; }
230
for(j=0,max=0; j < natom; j++) /* find the overall maximum */
231
if(fabs(charge[j]) >= fabs(charge[max])) max = j;
232
rank[0] = max; boolean[max] = 1;
233
for(j=1; j < natom; j++) {
235
while(boolean[max]) max++; /* find an unused max */
236
for(k=0; k < natom; k++)
237
if((fabs(charge[k]) >= fabs(charge[max])) && !boolean[k]) max = k;
238
rank[j] = max; boolean[max] = 1;
241
/* Build the orbital's domain starting in order of decreasing charge contribution */
242
for(j=0; j < nso; j++) {
244
for(k=0; k < nso; k++)
245
SR[j] += S[j][k] * C[k][i];
248
domain[i-nfzc][rank[0]] = 1; /* at least one atom must be in the domain */
249
domain_len[i-nfzc] = 1;
253
while(fabs(fR[i-nfzc]) > local.cutoff) {
255
/* Completeness check */
256
for(j=0,row=0; j < natom; j++) {
257
if(domain[i-nfzc][j]) {
258
for(k=aostart[j]; k <= aostop[j]; k++,row++) {
262
for(l=0,col=0; l < natom; l++) {
263
if(domain[i-nfzc][l]) {
265
for(m=aostart[l]; m <= aostop[l]; m++,col++)
266
X[row][col] = S[k][m];
275
errcod = C_DGESV(row, 1, &(X[0][0]), nso, &(ipiv[0]), &(Z[0]), nso);
277
fprintf(outfile, "\nError in DGESV return in orbital domain construction.\n");
278
exit(PSI_RETURN_FAILURE);
282
for(j=0,row=0; j < natom; j++) {
283
if(domain[i-nfzc][j]) {
284
for(k=aostart[j]; k <= aostop[j]; k++,row++)
285
for(l=0; l < nso; l++) fR[i-nfzc] -= Z[row] * S[k][l] * C[l][i];
289
/* Augment the domain if necessary */
290
if(fabs(fR[i-nfzc]) > local.cutoff) {
291
domain[i-nfzc][rank[next_atom++]] = 1;
292
domain_len[i-nfzc]++;
299
/* Allow user input of selected domains */
300
if(ip_exist("DOMAINS",0)) {
301
ip_count("DOMAINS", &num_entries,0);
302
for(i=0; i < num_entries; i++) {
303
ip_count("DOMAINS", &entry_len, 1, i);
304
ip_data("DOMAINS", "%d", &orbital, 2, i, 0);
306
/* Clear out the current domain for this orbital */
307
for(j=0; j < natom; j++) domain[orbital][j] = 0;
308
domain_len[orbital] = 0;
310
for(j=1; j < entry_len; j++) {
311
errcod = ip_data("DOMAINS","%d", &atom,2,i,j);
312
domain[orbital][atom] = 1;
313
domain_len[orbital]++;
318
/* Print the orbital domains */
320
for(i=0; i < nocc; i++)
321
if(domain_len[i] > max) max = domain_len[i];
323
fprintf(outfile, "\n ****** Occupied Orbital Domains ******\n");
324
fprintf(outfile, " Orbital Domain");
325
for(i=0; i < max-2; i++) fprintf(outfile, " "); /* formatting junk */
326
fprintf(outfile, " Completeness\n");
327
fprintf(outfile, " ------- ------");
328
for(i=0; i < max-2; i++) fprintf(outfile, "---"); /* more formatting junk */
329
fprintf(outfile, " ------------\n");
330
for(i=0; i < nocc; i++) {
331
fprintf(outfile, " %2d ",i);
332
for(j=0,cnt=0; j < natom; j++) if(domain[i][j]) { fprintf(outfile, " %2d", j); cnt++; }
333
if(cnt < max) for(; cnt < max; cnt++) fprintf(outfile, " ");
334
fprintf(outfile, " %7.5f\n", fR[i]);
338
/* Build the pair domains */
339
pairdomain = init_int_matrix(nocc*nocc,natom);
340
pairdom_len = init_int_array(nocc*nocc);
341
for(i=0,ij=0; i < nocc; i++)
342
for(j=0; j < nocc; j++,ij++)
343
for(k=0; k < natom; k++) {
344
if(domain[i][k] || domain[j][k]) {
345
pairdomain[ij][k] = 1;
346
pairdom_len[ij] += aostop[k] - aostart[k] + 1;
350
/* Identify and/or remove weak pairs */
351
weak_pairs = init_int_array(nocc*nocc);
352
fprintf(outfile, "\n");
353
for(i=0,ij=0; i < nocc; i++)
354
for(j=0; j < nocc; j++,ij++) {
356
for(k=0; k < natom; k++)
357
if(domain[i][k] && domain[j][k]) weak = 0;
359
if(weak && strcmp(local.weakp,"NONE")) {
362
if(!strcmp(local.weakp,"NEGLECT")) {
363
fprintf(outfile, "\tPair %d %d = [%d] is weak and will be deleted.\n", i, j, ij);
366
else weak_pairs[ij] = 0;
369
/* Compute the total number of singles and doubles */
370
/* replacement code from TDC on 11-5-02 */
371
t1_length = t2_length = 0;
372
for(i=0,ij=0; i < nocc; i++) {
373
for(k=0; k < natom; k++) {
375
for(a=aostart[k]; a <= aostop[k]; a++) t1_length++;
377
for(j=0; j < nocc; j++,ij++) {
378
for(k=0; k < natom; k++) {
379
for(l=0; l < natom; l++) {
380
if(pairdomain[ij][k] && pairdomain[ij][l] && !weak_pairs[ij]) {
381
for(a=aostart[k]; a <= aostop[k]; a++)
382
for(b=aostart[l]; b <= aostop[l]; b++)
390
/* Print excitation space reduction info */
391
fprintf(outfile, "\n\tL1 Length = %d (local), %d (canonical)\n",
392
t1_length, nocc*nvir);
393
fprintf(outfile, "\tL2 Length = %d (local), %d (canonical)\n\n",
394
t2_length, nocc*nocc*nvir*nvir);
397
local.domain = domain;
398
local.pairdomain = pairdomain;
399
local.pairdom_len = pairdom_len;
400
local.weak_pairs = weak_pairs;
401
local.aostart = aostart;
402
local.aostop = aostop;
415
ip_boolean("DOMAIN_PRINT",&(print_test),0);
417
fprintf(outfile, "Printing of orbital domains requested...exiting.\n\n");
418
exit(PSI_RETURN_FAILURE);
421
/************* Orbital Domains Complete ***************/
423
/* Compute the complete virtual space projector */
424
Rt_full = block_matrix(nso,nso);
425
for(i=0; i < nso; i++) Rt_full[i][i] = 1.0;
427
C_DGEMM('n','n',nso,nso,nso,-1.0,&(D[0][0]),nso,&(S[0][0]),nso,
428
1.0,&(Rt_full[0][0]),nso);
431
fprintf(outfile, "\n\tVirtual-Space Projector (R-tilde):\n");
432
print_mat(Rt_full, nso, nso, outfile);
435
/* Compute the norm of each PAO */
436
for(i=0; i < nso; i++) {
438
for(j=0; j < nso; j++) {
439
norm += Rt_full[j][i] * Rt_full[j][i];
443
fprintf(outfile, "\tNorm of orbital %4d = %20.12f...deleteing\n", i, norm);
444
for(j=0; j < nso; j++) Rt_full[j][i] = 0.0;
447
fprintf(outfile, "\n");
450
/* Grab the MO-basis Fock matrix */
451
Fmo = block_matrix(nso, nso);
452
for(i=0; i < nfzc; i++) Fmo[i][i] = eps_all[i];
453
dpd_file2_init(&fock, CC_OEI, 0, 0, 0, "fIJ");
454
dpd_file2_mat_init(&fock);
455
dpd_file2_mat_rd(&fock);
456
for(i=0; i < nocc; i++)
457
for(j=0; j < nocc; j++)
458
Fmo[i+nfzc][j+nfzc] = fock.matrix[0][i][j];
459
dpd_file2_mat_close(&fock);
460
dpd_file2_close(&fock);
462
dpd_file2_init(&fock, CC_OEI, 0, 1, 1, "fAB");
463
dpd_file2_mat_init(&fock);
464
dpd_file2_mat_rd(&fock);
465
for(i=0; i < nvir; i++)
466
for(j=0; j < nvir; j++)
467
Fmo[i+nfzc+nocc][j+nfzc+nocc] = fock.matrix[0][i][j];
468
dpd_file2_mat_close(&fock);
469
dpd_file2_close(&fock);
472
fprintf(outfile, "\n\tMO Basis Fock matrix:\n");
473
print_mat(Fmo, nso, nso, outfile);
476
/* Build the AO-basis Fock matrix */
477
F = block_matrix(nso,nso);
478
C_DGEMM('t','n',nso,nso,nso,1.0,&(Ci[0][0]),nso,&(Fmo[0][0]),nso,
480
C_DGEMM('n','n',nso,nso,nso,1.0,&(X[0][0]),nso,&(Ci[0][0]),nso,
483
/* Build the occupied orbital energy list */
484
eps_occ = init_array(nocc);
485
for(i=0;i < nocc; i++) eps_occ[i] = Fmo[i+nfzc][i+nfzc];
488
fprintf(outfile, "\n\tAO-Basis Fock Matrix:\n");
489
print_mat(F, nso, nso, outfile);
492
/* Compute R^+ S for virtual orbitals */
493
RS = block_matrix(nvir,nso);
494
for(a=0; a < nvir; a++)
495
for(i=0; i < nso; i++)
496
X[i][a] = C[i][a+nocc_all];
497
C_DGEMM('t','n',nvir,nso,nso,1.0,&(X[0][0]),nso,&(S[0][0]),nso,0.0,&(RS[0][0]),nso);
499
/* Build the virtual metric and W transforms for each pair domain */
500
Rt = block_matrix(nso, nso);
501
W = (double ***) malloc(nocc * nocc * sizeof(double **));
502
V = (double ***) malloc(nocc * nocc * sizeof(double **));
503
eps_vir = (double **) malloc(nocc * nocc * sizeof(double *));
504
pairdom_nrlen = init_int_array(nocc * nocc); /* dimension of non-redundant basis */
507
for(ij=0; ij < nocc * nocc; ij++) {
509
if(pairdom_len[ij]) {
511
zero_mat(Rt, nso, nso);
513
/* Build the virtual space projector for this pair */
514
for(k=0,L=0; k < natom; k++) {
515
if(pairdomain[ij][k]) {
516
for(l=aostart[k]; l <= aostop[k]; l++,L++) {
517
for(m=0; m < nso; m++) {
518
Rt[m][L] = Rt_full[m][l];
524
/* Compute the MO -> projected virtual transformation matrix */
525
V[ij] = block_matrix(nvir,pairdom_len[ij]);
526
C_DGEMM('n','n',nvir,pairdom_len[ij],nso,1.0,&(RS[0][0]),nso,&(Rt[0][0]),nso,0.0,
527
&(V[ij][0][0]),pairdom_len[ij]);
530
fprintf(outfile, "\nV[%d]:\n", ij);
531
fprintf(outfile, "======\n");
532
print_mat(V[ij], nvir, pairdom_len[ij], outfile);
535
/* Virtual space metric */
536
St = block_matrix(pairdom_len[ij],pairdom_len[ij]);
537
C_DGEMM('n','n',nso,pairdom_len[ij],nso,1.0,&(S[0][0]),nso,&(Rt[0][0]),nso,
539
C_DGEMM('t','n',pairdom_len[ij],pairdom_len[ij],nso,1.0,&(Rt[0][0]),nso,&(X[0][0]),nso,
540
0.0,&(St[0][0]),pairdom_len[ij]);
543
fprintf(outfile, "\n\tVirtual-Space Metric (S-tilde) for ij = %d:\n", ij);
544
print_mat(St, pairdom_len[ij], pairdom_len[ij], outfile);
547
/* Diagonalize metric */
548
evals = init_array(pairdom_len[ij]);
549
evecs = block_matrix(pairdom_len[ij],pairdom_len[ij]);
550
sq_rsp(pairdom_len[ij],pairdom_len[ij],St,evals,1,evecs,1e-12);
552
/* Count the number of zero eigenvalues */
553
for(i=0,cnt=0; i < pairdom_len[ij]; i++) if(evals[i] <= 1e-6) cnt++;
555
pairdom_nrlen[ij] = pairdom_len[ij]-cnt;
558
fprintf(outfile, "\n\tS-tilde eigenvalues for ij = %d:\n", ij);
559
for(i=0; i < pairdom_len[ij]; i++) fprintf(outfile, "\t%d %20.12f\n", i, evals[i]);
561
fprintf(outfile, "\n\tS-tilde eigenvectors for ij = %d:\n", ij);
562
print_mat(evecs,pairdom_len[ij],pairdom_len[ij],outfile);
565
/* Build the projected, non-redundant transform (X-tilde) */
566
Xt = block_matrix(pairdom_len[ij],pairdom_nrlen[ij]);
567
for(i=0,I=0; i < pairdom_len[ij]; i++) {
568
if(evals[i] > 1e-6) {
569
for(j=0; j < pairdom_len[ij]; j++)
570
Xt[j][I] = evecs[j][i]/sqrt(evals[i]);
578
fprintf(outfile, "\n\tTransform to non-redundant, projected virtuals (X-tilde) for ij = %d:\n", ij);
579
print_mat(Xt, pairdom_len[ij], pairdom_nrlen[ij], outfile);
585
/* Build the projected (redundant) virtual Fock matrix */
586
Ft = block_matrix(pairdom_len[ij], pairdom_len[ij]);
587
C_DGEMM('t','n',pairdom_len[ij],nso,nso,1.0,&(Rt[0][0]),nso,&(F[0][0]),nso,
589
C_DGEMM('n','n',pairdom_len[ij],pairdom_len[ij],nso,1.0,&(X[0][0]),nso,&(Rt[0][0]),nso,
590
0.0,&(Ft[0][0]),pairdom_len[ij]);
592
/* Project the Fock matrix into the non-redundant virtual space */
593
Fbar = block_matrix(pairdom_nrlen[ij],pairdom_nrlen[ij]);
594
C_DGEMM('t','n',pairdom_nrlen[ij],pairdom_len[ij],pairdom_len[ij],1.0,
595
&(Xt[0][0]),pairdom_nrlen[ij],&(Ft[0][0]),pairdom_len[ij],0.0,&(X[0][0]),nso);
596
C_DGEMM('n','n',pairdom_nrlen[ij],pairdom_nrlen[ij],pairdom_len[ij],1.0,
597
&(X[0][0]),nso,&(Xt[0][0]),pairdom_nrlen[ij],0.0,&(Fbar[0][0]),pairdom_nrlen[ij]);
600
fprintf(outfile, "\n\tFbar matrix for ij = %d:\n", ij);
601
print_mat(Fbar,pairdom_nrlen[ij],pairdom_nrlen[ij],outfile);
604
/* Diagonalize Fbar */
605
evals = init_array(pairdom_nrlen[ij]);
606
evecs = block_matrix(pairdom_nrlen[ij],pairdom_nrlen[ij]);
607
sq_rsp(pairdom_nrlen[ij],pairdom_nrlen[ij],Fbar,evals,1,evecs,1e-12);
610
fprintf(outfile, "\n\tFbar eigenvectors for ij = %d:\n", ij);
611
print_mat(evecs,pairdom_nrlen[ij],pairdom_nrlen[ij],outfile);
614
/* Finally, build the W matrix */
615
W[ij] = block_matrix(pairdom_len[ij],pairdom_nrlen[ij]);
616
C_DGEMM('n','n',pairdom_len[ij],pairdom_nrlen[ij],pairdom_nrlen[ij],1.0,
617
&(Xt[0][0]),pairdom_nrlen[ij],&(evecs[0][0]),pairdom_nrlen[ij],
618
0.0,&(W[ij][0][0]),pairdom_nrlen[ij]);
622
fprintf(outfile, "\n\tW Transformation Matrix for ij = %d:\n", ij);
623
print_mat(W,pairdom_len[ij],pairdom_nrlen[ij],outfile);
626
/* build the orbital energy list */
627
eps_vir[ij] = init_array(pairdom_nrlen[ij]);
628
for(i=0; i < pairdom_nrlen[ij]; i++)
629
eps_vir[ij][i] = evals[i]; /* virtual orbital energies */
632
fprintf(outfile, "\n\tVirtual orbital Energies for ij = %d:\n", ij);
633
for(i=0; i < pairdom_nrlen[ij]; i++)
634
fprintf(outfile, "%d %20.12f\n", i, eps_vir[ij][i]);
645
} /* if(pairdom_len[ij]) */
663
local.eps_occ = eps_occ;
664
local.eps_vir = eps_vir;
665
local.pairdom_nrlen = pairdom_nrlen;
42
local.nso = moinfo.nso;
43
local.nocc = moinfo.occpi[0]; /* active doubly occupied orbitals */
44
local.nvir = moinfo.virtpi[0]; /* active virtual orbitals */
667
46
fprintf(outfile, "\tLocalization parameters ready.\n\n");
677
50
void local_done(void)
682
for(i=0; i < local.nocc*local.nocc; i++) {
683
if(local.pairdom_len[i]) {
684
free_block(local.W[i]);
685
free_block(local.V[i]);
686
free(local.eps_vir[i]);
696
free_int_matrix(local.pairdomain, local.nocc*local.nocc);
697
free_int_matrix(local.domain, local.nocc);
698
free(local.pairdom_len);
699
free(local.pairdom_nrlen);
700
free(local.weak_pairs);
702
52
fprintf(outfile, "\tLocal parameters free.\n");
705
55
void local_filter_T1(dpdfile2 *T1)
709
int *pairdom_len, *pairdom_nrlen;
710
double ***V, ***W, *eps_occ, **eps_vir;
711
59
double *T1tilde, *T1bar;
713
62
nocc = local.nocc;
714
63
nvir = local.nvir;
717
eps_occ = local.eps_occ;
718
eps_vir = local.eps_vir;
719
pairdom_len = local.pairdom_len;
720
pairdom_nrlen = local.pairdom_nrlen;
65
local.pairdom_len = init_int_array(nocc*nocc);
66
local.pairdom_nrlen = init_int_array(nocc*nocc);
67
local.eps_occ = init_array(nocc);
68
psio_read_entry(CC_INFO, "Local Pair Domain Length", (char *) local.pairdom_len,
69
nocc*nocc*sizeof(int));
70
psio_read_entry(CC_INFO, "Local Pair Domain NR Length", (char *) local.pairdom_nrlen,
71
nocc*nocc*sizeof(int));
72
psio_read_entry(CC_INFO, "Local Occupied Orbital Energies", (char *) local.eps_occ,
75
local.W = (double ***) malloc(nocc * nocc * sizeof(double **));
76
local.V = (double ***) malloc(nocc * nocc * sizeof(double **));
77
local.eps_vir = (double **) malloc(nocc * nocc * sizeof(double *));
79
for(ij=0; ij < nocc*nocc; ij++) {
80
local.eps_vir[ij] = init_array(local.pairdom_nrlen[ij]);
81
psio_read(CC_INFO, "Local Virtual Orbital Energies", (char *) local.eps_vir[ij],
82
local.pairdom_nrlen[ij]*sizeof(double), next, &next);
85
for(ij=0; ij < nocc*nocc; ij++) {
86
local.V[ij] = block_matrix(nvir,local.pairdom_len[ij]);
87
psio_read(CC_INFO, "Local Residual Vector (V)", (char *) local.V[ij][0],
88
nvir*local.pairdom_len[ij]*sizeof(double), next, &next);
91
for(ij=0; ij < nocc*nocc; ij++) {
92
local.W[ij] = block_matrix(local.pairdom_len[ij],local.pairdom_nrlen[ij]);
93
psio_read(CC_INFO, "Local Transformation Matrix (W)", (char *) local.W[ij][0],
94
local.pairdom_len[ij]*local.pairdom_nrlen[ij]*sizeof(double), next, &next);
722
97
dpd_file2_mat_init(T1);
723
98
dpd_file2_mat_rd(T1);
790
204
X2 = block_matrix(nvir,nso);
791
205
T2tilde = block_matrix(nso,nso);
792
206
T2bar = block_matrix(nvir, nvir);
793
208
for(i=0,ij=0; i < nocc; i++) {
794
209
for(j=0; j < nocc; j++,ij++) {
796
if(!weak_pairs[ij]) {
211
if(!local.weak_pairs[ij]) {
798
213
/* Transform the virtuals to the redundant projected virtual basis */
799
C_DGEMM('t', 'n', pairdom_len[ij], nvir, nvir, 1.0, &(V[ij][0][0]), pairdom_len[ij],
214
C_DGEMM('t', 'n', local.pairdom_len[ij], nvir, nvir, 1.0, &(local.V[ij][0][0]), local.pairdom_len[ij],
800
215
&(T2->matrix[0][ij][0]), nvir, 0.0, &(X1[0][0]), nvir);
801
C_DGEMM('n', 'n', pairdom_len[ij], pairdom_len[ij], nvir, 1.0, &(X1[0][0]), nvir,
802
&(V[ij][0][0]), pairdom_len[ij], 0.0, &(T2tilde[0][0]), nso);
216
C_DGEMM('n', 'n', local.pairdom_len[ij], local.pairdom_len[ij], nvir, 1.0, &(X1[0][0]), nvir,
217
&(local.V[ij][0][0]), local.pairdom_len[ij], 0.0, &(T2tilde[0][0]), nso);
804
219
/* Transform the virtuals to the non-redundant virtual basis */
805
C_DGEMM('t', 'n', pairdom_nrlen[ij], pairdom_len[ij], pairdom_len[ij], 1.0,
806
&(W[ij][0][0]), pairdom_nrlen[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
807
C_DGEMM('n', 'n', pairdom_nrlen[ij], pairdom_nrlen[ij], pairdom_len[ij], 1.0,
808
&(X2[0][0]), nso, &(W[ij][0][0]), pairdom_nrlen[ij], 0.0, &(T2bar[0][0]), nvir);
220
C_DGEMM('t', 'n', local.pairdom_nrlen[ij], local.pairdom_len[ij], local.pairdom_len[ij], 1.0,
221
&(local.W[ij][0][0]), local.pairdom_nrlen[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
222
C_DGEMM('n', 'n', local.pairdom_nrlen[ij], local.pairdom_nrlen[ij], local.pairdom_len[ij], 1.0,
223
&(X2[0][0]), nso, &(local.W[ij][0][0]), local.pairdom_nrlen[ij], 0.0, &(T2bar[0][0]), nvir);
810
225
/* Divide the new amplitudes by the denominators */
811
for(a=0; a < pairdom_nrlen[ij]; a++) {
812
for(b=0; b < pairdom_nrlen[ij]; b++) {
813
T2bar[a][b] /= (eps_occ[i] + eps_occ[j] - eps_vir[ij][a] - eps_vir[ij][b]);
226
for(a=0; a < local.pairdom_nrlen[ij]; a++) {
227
for(b=0; b < local.pairdom_nrlen[ij]; b++) {
228
T2bar[a][b] /= (local.eps_occ[i] + local.eps_occ[j]
229
- local.eps_vir[ij][a] - local.eps_vir[ij][b]);
817
233
/* Transform the new T2's to the redundant virtual basis */
818
C_DGEMM('n', 'n', pairdom_len[ij], pairdom_nrlen[ij], pairdom_nrlen[ij], 1.0,
819
&(W[ij][0][0]), pairdom_nrlen[ij], &(T2bar[0][0]), nvir, 0.0, &(X1[0][0]), nvir);
820
C_DGEMM('n','t', pairdom_len[ij], pairdom_len[ij], pairdom_nrlen[ij], 1.0,
821
&(X1[0][0]), nvir, &(W[ij][0][0]), pairdom_nrlen[ij], 0.0, &(T2tilde[0][0]), nso);
234
C_DGEMM('n', 'n', local.pairdom_len[ij], local.pairdom_nrlen[ij], local.pairdom_nrlen[ij], 1.0,
235
&(local.W[ij][0][0]), local.pairdom_nrlen[ij], &(T2bar[0][0]), nvir, 0.0, &(X1[0][0]), nvir);
236
C_DGEMM('n','t', local.pairdom_len[ij], local.pairdom_len[ij], local.pairdom_nrlen[ij], 1.0,
237
&(X1[0][0]), nvir, &(local.W[ij][0][0]), local.pairdom_nrlen[ij], 0.0, &(T2tilde[0][0]), nso);
823
239
/* Transform the new T2's to the MO basis */
824
C_DGEMM('n', 'n', nvir, pairdom_len[ij], pairdom_len[ij], 1.0,
825
&(V[ij][0][0]), pairdom_len[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
826
C_DGEMM('n', 't', nvir, nvir, pairdom_len[ij], 1.0, &(X2[0][0]), nso,
827
&(V[ij][0][0]), pairdom_len[ij], 0.0, &(T2->matrix[0][ij][0]), nvir);
240
C_DGEMM('n', 'n', nvir, local.pairdom_len[ij], local.pairdom_len[ij], 1.0,
241
&(local.V[ij][0][0]), local.pairdom_len[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
242
C_DGEMM('n', 't', nvir, nvir, local.pairdom_len[ij], 1.0, &(X2[0][0]), nso,
243
&(local.V[ij][0][0]), local.pairdom_len[ij], 0.0, &(T2->matrix[0][ij][0]), nvir);
829
245
else /* This must be a neglected weak pair; force it to zero */
830
246
memset((void *) T2->matrix[0][ij], 0, nvir*nvir*sizeof(double));
added
removed
src/bin/cclambda/local.c
