3
\brief Enter brief description of file here
10
#include <libipv1/ip_lib.h>
11
#include <libciomr/libciomr.h>
12
#include <libpsio/psio.h>
13
#include <libiwl/iwl.h>
14
#include <libint/libint.h>
15
#include <libchkpt/chkpt.h>
17
#include <libdpd/dpd.h>
25
namespace psi { namespace cis {
28
** local_init(): Set up parameters of local excitation domains.
30
** The orbital domains constructed here are based on those described
31
** in Broughton and Pulay, J. Comp. Chem. 14, 736-740 (1993). The
32
** localization of the occupied orbitals is done elsewhere (see the
33
** program "local"). Pair domains are defined as the union of pairs
34
** of single occupied orbital domains. "Weak pairs", which are
35
** defined as pair domains whose individual occupied orbital domains
36
** have no atoms in common, are identified (cf. int *weak_pairs).
43
int i, j, k, ij, stat, a, b, r, l, I, L;
44
int nmo, nso, nocc, nocc_all, nvir, noei, nirreps, nfzc;
45
double **C; /* AO -> localized MO transformation matrix */
46
double **Ci; /* localized MO -> AO transformation matrix */
47
double **D; /* 1/2 SCF closed-shell density matrix (AO) */
48
double **Rt, **Rt_full; /* Projected, redundant virtual transform (R-tilde) */
49
double **S; /* AO overlap */
50
double **St; /* Projected virtual overlap */
51
double **Xt; /* Projected, non-redundant virtual transform (X-tilde) */
52
double ***V; /* MO -> projected, redundant virtual transform */
53
double **Fmo;/* MO basis Fock matrix */
54
double **F; /* AO basis Fock matrix */
55
double **Ft; /* Projected, redundant virtual Fock matrix */
56
double **Fbar; /* Projected, non-redundant virtual Fock matrix */
57
double ***W; /* Transformation matrix from tilde -> bar for each ij pair*/
58
double *eps_occ; /* occupied orbital energies for local denominators */
59
double **eps_vir; /* virtual orbital energies for local denominators */
61
double *evals, **evecs;
62
double *eps_all; /* All MO energies */
63
dpdfile2 fock, FMI, FAE;
65
int natom, atom, am, offset, nshell, shell_length, next_atom;
66
int row, col, max, m, errcod, cnt;
67
int *rank, *boolean, *ipiv;
68
int *l_length, *aostart, *aostop, *ao2atom;
70
int **domain, *domain_len, **pairdomain, *pairdom_len, *pairdom_nrlen;
72
double *fR, cutoff, *charge, *SR, *Z, tmp, *ss;
74
int print_test, num_entries, entry_len, orbital;
75
int t1_length, t2_length, puream, weak;
79
double **U, **R, **WW, **RS;
81
chkpt_init(PSIO_OPEN_OLD);
83
natom = chkpt_rd_natom();
84
nshell = chkpt_rd_nshell();
85
puream = chkpt_rd_puream();
86
eps_all = chkpt_rd_evals();
87
stype = chkpt_rd_stype();
88
snuc = chkpt_rd_snuc();
93
/* C1 symmetry only */
94
nirreps = moinfo.nirreps;
96
fprintf(outfile, "\nError: localization must use C1 symmetry.\n");
97
exit(PSI_RETURN_FAILURE);
101
nmo = moinfo.nmo; /* should be the same as nso */
103
fprintf(outfile, "\nError: NMO != NSO! %d != %d\n", nmo, nso);
104
exit(PSI_RETURN_FAILURE);
107
nocc = moinfo.occpi[0]; /* active doubly occupied orbitals */
108
nfzc = moinfo.frdocc[0]; /* frozen doubly occupied orbitals */
109
nocc_all = nocc + nfzc; /* all doubly occupied orbitals */
110
nvir = moinfo.virtpi[0]; /* active virtual orbitals */
117
/* A couple of scratch arrays */
118
X = block_matrix(nso, nso);
119
Y = block_matrix(nso, nso);
122
Ci = block_matrix(nso, nso);
123
for(i=0; i < nso; i++)
124
for(j=0; j < nso; j++)
127
invert_matrix(C, Ci, nso, outfile);
129
for(i=0; i < nso; i++)
130
for(j=0; j < nso; j++)
134
fprintf(outfile, "\n\tC inverse (Ci):\n");
135
print_mat(Ci, nso, nso, outfile);
138
/* Compute the AO-basis overlap integrals */
140
S = block_matrix(nso,nso);
141
C_DGEMM('t','n',nso,nso,nso,1.0,&(Ci[0][0]),nso,&(Ci[0][0]),nso,
145
/* Get the overlap integrals -- these should be identical to AO S */
146
noei = nso*(nso+1)/2;
147
ss = init_array(noei);
148
stat = iwl_rdone(PSIF_OEI,PSIF_SO_S,ss,noei,0,0,outfile);
149
S = block_matrix(nso,nso);
150
for(i=0,ij=0; i < nso; i++)
151
for(j=0; j <= i; j++,ij++) {
152
S[i][j] = S[j][i] = ss[ij];
157
fprintf(outfile, "\n\tAO Overlap (S)\n");
158
print_mat(S, nso, nso, outfile);
161
/* Build the SCF closed-shell density matrix/2 */
162
D = block_matrix(nso,nso);
163
for(i=0; i < nso; i++)
164
for(j=0; j < nso; j++)
165
for(k=0; k < nocc_all; k++)
166
D[i][j] += C[i][k] * C[j][k];
169
fprintf(outfile, "\n\tAO-basis SCF Density (D):\n");
170
print_mat(D, nso, nso, outfile);
174
/* Compute the length of each AM block */
175
l_length = init_int_array(LIBINT_MAX_AM);
177
for(l=1; l < LIBINT_MAX_AM; l++) {
178
if(puream) l_length[l] = 2 * l + 1;
179
else l_length[l] = l_length[l-1] + l + 1;
182
/* Set up the atom->AO and AO->atom lookups */
183
aostart = init_int_array(natom);
184
aostop = init_int_array(natom);
185
for(i=0,atom=-1,offset=0; i < nshell; i++) {
187
shell_length = l_length[am];
189
if(atom != snuc[i] - 1) {
190
if(atom != -1) aostop[atom] = offset-1;
192
aostart[atom] = offset;
194
offset += shell_length;
196
aostop[atom] = offset-1;
198
ao2atom = init_int_array(nso);
199
for(i=0; i < natom; i++)
200
for(j=aostart[i]; j <= aostop[i]; j++) ao2atom[j] = i;
205
/************* Build the orbital domains ************/
207
domain = init_int_matrix(nocc, natom);
208
domain_len = init_int_array(nocc);
209
charge = init_array(natom);
210
rank = init_int_array(natom);
211
boolean = init_int_array(natom);
212
SR = init_array(nso);
214
ipiv = init_int_array(nso);
215
fR = init_array(nocc);
217
for(i=nfzc; i < nocc_all; i++) {
219
/* Compute the contribution of each atom to this orbital's charge/population */
220
for(j=0; j < natom; j++) {
222
for(k=aostart[j]; k <= aostop[j]; k++) {
224
for(l=0; l < nso; l++) tmp += S[k][l] * C[l][i];
230
/* Rank the atomic contributions to the orbital's charge */
231
for(j=0; j < natom; j++) { rank[j] = 0; boolean[j] = 0; }
232
for(j=0,max=0; j < natom; j++) /* find the overall maximum */
233
if(fabs(charge[j]) >= fabs(charge[max])) max = j;
234
rank[0] = max; boolean[max] = 1;
235
for(j=1; j < natom; j++) {
237
while(boolean[max]) max++; /* find an unused max */
238
for(k=0; k < natom; k++)
239
if((fabs(charge[k]) >= fabs(charge[max])) && !boolean[k]) max = k;
240
rank[j] = max; boolean[max] = 1;
243
/* Build the orbital's domain starting in order of decreasing charge contribution */
244
for(j=0; j < nso; j++) {
246
for(k=0; k < nso; k++)
247
SR[j] += S[j][k] * C[k][i];
250
domain[i-nfzc][rank[0]] = 1; /* at least one atom must be in the domain */
251
domain_len[i-nfzc] = 1;
255
while(fabs(fR[i-nfzc]) > local.cutoff) {
257
/* Completeness check */
258
for(j=0,row=0; j < natom; j++) {
259
if(domain[i-nfzc][j]) {
260
for(k=aostart[j]; k <= aostop[j]; k++,row++) {
264
for(l=0,col=0; l < natom; l++) {
265
if(domain[i-nfzc][l]) {
267
for(m=aostart[l]; m <= aostop[l]; m++,col++)
268
X[row][col] = S[k][m];
277
errcod = C_DGESV(row, 1, &(X[0][0]), nso, &(ipiv[0]), &(Z[0]), nso);
279
fprintf(outfile, "\nError in DGESV return in orbital domain construction.\n");
280
exit(PSI_RETURN_FAILURE);
284
for(j=0,row=0; j < natom; j++) {
285
if(domain[i-nfzc][j]) {
286
for(k=aostart[j]; k <= aostop[j]; k++,row++)
287
for(l=0; l < nso; l++) fR[i-nfzc] -= Z[row] * S[k][l] * C[l][i];
291
/* Augment the domain if necessary */
292
if(fabs(fR[i-nfzc]) > local.cutoff) {
293
domain[i-nfzc][rank[next_atom++]] = 1;
294
domain_len[i-nfzc]++;
301
/* Allow user adjustment of selected domains */
302
if(ip_exist("DOMAINS",0)) {
303
ip_count("DOMAINS", &num_entries,0);
304
for(i=0; i < num_entries; i++) {
305
ip_count("DOMAINS", &entry_len, 1, i);
306
ip_data("DOMAINS", "%d", &orbital, 2, i, 0);
308
/* Clear out the current domain for this orbital */
310
for(j=0; j < natom; j++) domain[orbital][j] = 0;
311
domain_len[orbital] = 0;
314
for(j=1; j < entry_len; j++) {
315
errcod = ip_data("DOMAINS","%d", &atom,2,i,j);
316
domain[orbital][atom] = 1;
317
domain_len[orbital]++;
322
/* Print the orbital domains */
324
for(i=0; i < nocc; i++)
325
if(domain_len[i] > max) max = domain_len[i];
327
fprintf(outfile, "\n ****** Occupied Orbital Domains ******\n");
328
fprintf(outfile, " Orbital Domain");
329
for(i=0; i < max-2; i++) fprintf(outfile, " "); /* formatting junk */
330
fprintf(outfile, " Completeness\n");
331
fprintf(outfile, " ------- ------");
332
for(i=0; i < max-2; i++) fprintf(outfile, "---"); /* more formatting junk */
333
fprintf(outfile, " ------------\n");
334
for(i=0; i < nocc; i++) {
335
fprintf(outfile, " %2d ",i);
336
for(j=0,cnt=0; j < natom; j++) if(domain[i][j]) { fprintf(outfile, " %2d", j); cnt++; }
337
if(cnt < max) for(; cnt < max; cnt++) fprintf(outfile, " ");
338
fprintf(outfile, " %7.5f\n", fR[i]);
342
/* Build the pair domains */
343
pairdomain = init_int_matrix(nocc*nocc,natom);
344
pairdom_len = init_int_array(nocc*nocc);
345
for(i=0,ij=0; i < nocc; i++)
346
for(j=0; j < nocc; j++,ij++)
347
for(k=0; k < natom; k++) {
348
if(domain[i][k] || domain[j][k]) {
349
pairdomain[ij][k] = 1;
350
pairdom_len[ij] += aostop[k] - aostart[k] + 1;
354
/* Identify and/or remove weak pairs */
355
weak_pairs = init_int_array(nocc*nocc);
356
fprintf(outfile, "\n");
357
for(i=0,ij=0; i < nocc; i++)
358
for(j=0; j < nocc; j++,ij++) {
360
for(k=0; k < natom; k++)
361
if(domain[i][k] && domain[j][k] && local.ghost != k) weak = 0;
363
if(weak && strcmp(local.weakp,"NONE")) {
366
if(!strcmp(local.weakp,"MP2"))
367
fprintf(outfile, "\tPair %d %d [%d] is weak and will be treated with MP2.\n", i, j, ij);
368
else if(!strcmp(local.weakp,"NEGLECT")) {
369
fprintf(outfile, "\tPair %d %d = [%d] is weak and will be deleted.\n", i, j, ij);
372
else weak_pairs[ij] = 0;
375
/* Compute the total number of singles and doubles */
376
t1_length = t2_length = 0;
377
for(i=0,ij=0; i < nocc; i++) {
378
t1_length += domain_len[i];
379
for(j=0; j < nocc; j++,ij++) {
380
for(k=0; k < natom; k++) {
381
for(l=0; l < natom; l++) {
382
if(pairdomain[ij][k] && pairdomain[ij][l] && !weak_pairs[ij]) {
383
for(a=aostart[k]; a <= aostop[k]; a++)
384
for(b=aostart[l]; b <= aostop[l]; b++)
392
/* Print excitation space reduction info */
393
fprintf(outfile, "\n\tT1 Length = %d (local), %d (canonical)\n",
394
t1_length, nocc*nvir);
395
fprintf(outfile, "\tT2 Length = %d (local), %d (canonical)\n\n",
396
t2_length, nocc*nocc*nvir*nvir);
399
local.domain = domain;
400
local.pairdomain = pairdomain;
401
local.pairdom_len = pairdom_len;
402
local.weak_pairs = weak_pairs;
403
local.aostart = aostart;
404
local.aostop = aostop;
417
ip_boolean("DOMAIN_PRINT",&(print_test),0);
419
fprintf(outfile, "Printing of orbital domains requested...exiting.\n\n");
420
exit(PSI_RETURN_FAILURE);
423
/************* Orbital Domains Complete ***************/
425
/* Compute the complete virtual space projector */
426
Rt_full = block_matrix(nso,nso);
427
for(i=0; i < nso; i++) Rt_full[i][i] = 1.0;
429
C_DGEMM('n','n',nso,nso,nso,-1.0,&(D[0][0]),nso,&(S[0][0]),nso,
430
1.0,&(Rt_full[0][0]),nso);
433
fprintf(outfile, "\n\tVirtual-Space Projector (R-tilde):\n");
434
print_mat(Rt_full, nso, nso, outfile);
437
/* Compute the norm of each PAO */
438
for(i=0; i < nso; i++) {
440
for(j=0; j < nso; j++) {
441
norm += Rt_full[j][i] * Rt_full[j][i];
445
fprintf(outfile, "\tNorm of orbital %4d = %20.12f...deleteing\n", i, norm);
446
for(j=0; j < nso; j++) Rt_full[j][i] = 0.0;
449
fprintf(outfile, "\n");
452
/* Grab the MO-basis Fock matrix */
453
Fmo = block_matrix(nso, nso);
454
for(i=0; i < nfzc; i++) Fmo[i][i] = eps_all[i];
455
dpd_file2_init(&fock, CC_OEI, 0, 0, 0, "fIJ");
456
dpd_file2_mat_init(&fock);
457
dpd_file2_mat_rd(&fock);
458
for(i=0; i < nocc; i++)
459
for(j=0; j < nocc; j++)
460
Fmo[i+nfzc][j+nfzc] = fock.matrix[0][i][j];
461
dpd_file2_mat_close(&fock);
462
dpd_file2_close(&fock);
464
dpd_file2_init(&fock, CC_OEI, 0, 1, 1, "fAB");
465
dpd_file2_mat_init(&fock);
466
dpd_file2_mat_rd(&fock);
467
for(i=0; i < nvir; i++)
468
for(j=0; j < nvir; j++)
469
Fmo[i+nfzc+nocc][j+nfzc+nocc] = fock.matrix[0][i][j];
470
dpd_file2_mat_close(&fock);
471
dpd_file2_close(&fock);
474
fprintf(outfile, "\n\tMO Basis Fock matrix:\n");
475
print_mat(Fmo, nso, nso, outfile);
478
/* Build the AO-basis Fock matrix */
479
F = block_matrix(nso,nso);
480
C_DGEMM('t','n',nso,nso,nso,1.0,&(Ci[0][0]),nso,&(Fmo[0][0]),nso,
482
C_DGEMM('n','n',nso,nso,nso,1.0,&(X[0][0]),nso,&(Ci[0][0]),nso,
485
/* Build the occupied orbital energy list */
486
eps_occ = init_array(nocc);
487
for(i=0;i < nocc; i++) eps_occ[i] = Fmo[i+nfzc][i+nfzc];
490
fprintf(outfile, "\n\tAO-Basis Fock Matrix:\n");
491
print_mat(F, nso, nso, outfile);
494
/* Compute R^+ S for virtual orbitals */
495
RS = block_matrix(nvir,nso);
496
for(a=0; a < nvir; a++)
497
for(i=0; i < nso; i++)
498
X[i][a] = C[i][a+nocc_all];
499
C_DGEMM('t','n',nvir,nso,nso,1.0,&(X[0][0]),nso,&(S[0][0]),nso,0.0,&(RS[0][0]),nso);
501
/* Build the virtual metric and W transforms for each pair domain */
502
Rt = block_matrix(nso, nso);
503
W = (double ***) malloc(nocc * nocc * sizeof(double **));
504
V = (double ***) malloc(nocc * nocc * sizeof(double **));
505
eps_vir = (double **) malloc(nocc * nocc * sizeof(double *));
506
pairdom_nrlen = init_int_array(nocc * nocc); /* dimension of non-redundant basis */
509
for(ij=0; ij < nocc * nocc; ij++) {
511
if(pairdom_len[ij]) {
513
zero_mat(Rt, nso, nso);
515
/* Build the virtual space projector for this pair */
516
for(k=0,L=0; k < natom; k++) {
517
if(pairdomain[ij][k]) {
518
for(l=aostart[k]; l <= aostop[k]; l++,L++) {
519
for(m=0; m < nso; m++) {
520
Rt[m][L] = Rt_full[m][l];
526
/* Compute the MO -> projected virtual transformation matrix */
527
V[ij] = block_matrix(nvir,pairdom_len[ij]);
528
C_DGEMM('n','n',nvir,pairdom_len[ij],nso,1.0,&(RS[0][0]),nso,&(Rt[0][0]),nso,0.0,
529
&(V[ij][0][0]),pairdom_len[ij]);
531
/* Virtual space metric */
532
St = block_matrix(pairdom_len[ij],pairdom_len[ij]);
533
C_DGEMM('n','n',nso,pairdom_len[ij],nso,1.0,&(S[0][0]),nso,&(Rt[0][0]),nso,
535
C_DGEMM('t','n',pairdom_len[ij],pairdom_len[ij],nso,1.0,&(Rt[0][0]),nso,&(X[0][0]),nso,
536
0.0,&(St[0][0]),pairdom_len[ij]);
539
fprintf(outfile, "\n\tVirtual-Space Metric (S-tilde) for ij = %d:\n", ij);
540
print_mat(St, pairdom_len[ij], pairdom_len[ij], outfile);
543
/* Diagonalize metric */
544
evals = init_array(pairdom_len[ij]);
545
evecs = block_matrix(pairdom_len[ij],pairdom_len[ij]);
546
sq_rsp(pairdom_len[ij],pairdom_len[ij],St,evals,1,evecs,1e-12);
548
/* Count the number of zero eigenvalues */
549
for(i=0,cnt=0; i < pairdom_len[ij]; i++) if(evals[i] <= 1e-6) cnt++;
551
pairdom_nrlen[ij] = pairdom_len[ij]-cnt;
554
fprintf(outfile, "\n\tS-tilde eigenvalues for ij = %d:\n", ij);
555
for(i=0; i < pairdom_len[ij]; i++) fprintf(outfile, "\t%d %20.12f\n", i, evals[i]);
557
fprintf(outfile, "\n\tS-tilde eigenvectors for ij = %d:\n", ij);
558
print_mat(evecs,pairdom_len[ij],pairdom_len[ij],outfile);
561
/* Build the projected, non-redundant transform (X-tilde) */
562
Xt = block_matrix(pairdom_len[ij],pairdom_nrlen[ij]);
563
for(i=0,I=0; i < pairdom_len[ij]; i++) {
564
if(evals[i] > 1e-6) {
565
for(j=0; j < pairdom_len[ij]; j++)
566
Xt[j][I] = evecs[j][i]/sqrt(evals[i]);
574
fprintf(outfile, "\n\tTransform to non-redundant, projected virtuals (X-tilde) for ij = %d:\n", ij);
575
print_mat(Xt, pairdom_len[ij], pairdom_nrlen[ij], outfile);
581
/* Build the projected (redundant) virtual Fock matrix */
582
Ft = block_matrix(pairdom_len[ij], pairdom_len[ij]);
583
C_DGEMM('t','n',pairdom_len[ij],nso,nso,1.0,&(Rt[0][0]),nso,&(F[0][0]),nso,
585
C_DGEMM('n','n',pairdom_len[ij],pairdom_len[ij],nso,1.0,&(X[0][0]),nso,&(Rt[0][0]),nso,
586
0.0,&(Ft[0][0]),pairdom_len[ij]);
588
/* Project the Fock matrix into the non-redundant virtual space */
589
Fbar = block_matrix(pairdom_nrlen[ij],pairdom_nrlen[ij]);
590
C_DGEMM('t','n',pairdom_nrlen[ij],pairdom_len[ij],pairdom_len[ij],1.0,
591
&(Xt[0][0]),pairdom_nrlen[ij],&(Ft[0][0]),pairdom_len[ij],0.0,&(X[0][0]),nso);
592
C_DGEMM('n','n',pairdom_nrlen[ij],pairdom_nrlen[ij],pairdom_len[ij],1.0,
593
&(X[0][0]),nso,&(Xt[0][0]),pairdom_nrlen[ij],0.0,&(Fbar[0][0]),pairdom_nrlen[ij]);
596
fprintf(outfile, "\n\tFbar matrix for ij = %d:\n", ij);
597
print_mat(Fbar,pairdom_nrlen[ij],pairdom_nrlen[ij],outfile);
600
/* Diagonalize Fbar */
601
evals = init_array(pairdom_nrlen[ij]);
602
evecs = block_matrix(pairdom_nrlen[ij],pairdom_nrlen[ij]);
603
sq_rsp(pairdom_nrlen[ij],pairdom_nrlen[ij],Fbar,evals,1,evecs,1e-12);
606
fprintf(outfile, "\n\tFbar eigenvectors for ij = %d:\n", ij);
607
print_mat(evecs,pairdom_nrlen[ij],pairdom_nrlen[ij],outfile);
610
/* Finally, build the W matrix */
611
W[ij] = block_matrix(pairdom_len[ij],pairdom_nrlen[ij]);
612
C_DGEMM('n','n',pairdom_len[ij],pairdom_nrlen[ij],pairdom_nrlen[ij],1.0,
613
&(Xt[0][0]),pairdom_nrlen[ij],&(evecs[0][0]),pairdom_nrlen[ij],
614
0.0,&(W[ij][0][0]),pairdom_nrlen[ij]);
618
fprintf(outfile, "\n\tW Transformation Matrix for ij = %d:\n", ij);
619
print_mat(W,pairdom_len[ij],pairdom_nrlen[ij],outfile);
622
/* build the orbital energy list */
623
eps_vir[ij] = init_array(pairdom_nrlen[ij]);
624
for(i=0; i < pairdom_nrlen[ij]; i++)
625
eps_vir[ij][i] = evals[i]; /* virtual orbital energies */
628
fprintf(outfile, "\n\tVirtual orbital Energies for ij = %d:\n", ij);
629
for(i=0; i < pairdom_nrlen[ij]; i++)
630
fprintf(outfile, "%d %20.12f\n", i, eps_vir[ij][i]);
641
} /* if(pairdom_len[ij]) */
659
local.eps_occ = eps_occ;
660
local.eps_vir = eps_vir;
661
local.pairdom_nrlen = pairdom_nrlen;
663
local.weak_pair_energy = 0.0;
665
fprintf(outfile, "\tLocalization parameters ready.\n\n");
675
void local_done(void)
680
for(i=0; i < local.nocc*local.nocc; i++) {
681
if(local.pairdom_len[i]) {
682
free_block(local.W[i]);
683
free_block(local.V[i]);
684
free(local.eps_vir[i]);
694
free_int_matrix(local.pairdomain);
695
free_int_matrix(local.domain);
696
free(local.pairdom_len);
697
free(local.pairdom_nrlen);
698
free(local.weak_pairs);
700
fprintf(outfile, "\tLocal parameters free.\n");
703
void local_filter_T2(dpdbuf4 *T2)
705
int ij, i, j, a, b, ab;
707
int *pairdom_len, *pairdom_nrlen, *weak_pairs;
708
double ***V, ***W, *eps_occ, **eps_vir;
709
double **X1, **X2, **T2tilde, **T2bar;
716
eps_occ = local.eps_occ;
717
eps_vir = local.eps_vir;
718
pairdom_len = local.pairdom_len;
719
pairdom_nrlen = local.pairdom_nrlen;
720
weak_pairs = local.weak_pairs;
722
/* Grab the MO-basis T2's */
723
dpd_buf4_mat_irrep_init(T2, 0);
724
dpd_buf4_mat_irrep_rd(T2, 0);
726
X1 = block_matrix(nso,nvir);
727
X2 = block_matrix(nvir,nso);
728
T2tilde = block_matrix(nso,nso);
729
T2bar = block_matrix(nvir, nvir);
730
for(i=0,ij=0; i < nocc; i++) {
731
for(j=0; j < nocc; j++,ij++) {
733
/* Transform the virtuals to the redundant projected virtual basis */
734
C_DGEMM('t', 'n', pairdom_len[ij], nvir, nvir, 1.0, &(V[ij][0][0]), pairdom_len[ij],
735
&(T2->matrix[0][ij][0]), nvir, 0.0, &(X1[0][0]), nvir);
736
C_DGEMM('n', 'n', pairdom_len[ij], pairdom_len[ij], nvir, 1.0, &(X1[0][0]), nvir,
737
&(V[ij][0][0]), pairdom_len[ij], 0.0, &(T2tilde[0][0]), nso);
739
/* Transform the virtuals to the non-redundant virtual basis */
740
C_DGEMM('t', 'n', pairdom_nrlen[ij], pairdom_len[ij], pairdom_len[ij], 1.0,
741
&(W[ij][0][0]), pairdom_nrlen[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
742
C_DGEMM('n', 'n', pairdom_nrlen[ij], pairdom_nrlen[ij], pairdom_len[ij], 1.0,
743
&(X2[0][0]), nso, &(W[ij][0][0]), pairdom_nrlen[ij], 0.0, &(T2bar[0][0]), nvir);
745
/* Divide the new amplitudes by the denominators */
746
for(a=0; a < pairdom_nrlen[ij]; a++) {
747
for(b=0; b < pairdom_nrlen[ij]; b++) {
748
T2bar[a][b] /= (eps_occ[i] + eps_occ[j] - eps_vir[ij][a] - eps_vir[ij][b]);
752
/* Transform the new T2's to the redundant virtual basis */
753
C_DGEMM('n', 'n', pairdom_len[ij], pairdom_nrlen[ij], pairdom_nrlen[ij], 1.0,
754
&(W[ij][0][0]), pairdom_nrlen[ij], &(T2bar[0][0]), nvir, 0.0, &(X1[0][0]), nvir);
755
C_DGEMM('n','t', pairdom_len[ij], pairdom_len[ij], pairdom_nrlen[ij], 1.0,
756
&(X1[0][0]), nvir, &(W[ij][0][0]), pairdom_nrlen[ij], 0.0, &(T2tilde[0][0]), nso);
758
/* Transform the new T2's to the MO basis */
759
C_DGEMM('n', 'n', nvir, pairdom_len[ij], pairdom_len[ij], 1.0,
760
&(V[ij][0][0]), pairdom_len[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
761
C_DGEMM('n', 't', nvir, nvir, pairdom_len[ij], 1.0, &(X2[0][0]), nso,
762
&(V[ij][0][0]), pairdom_len[ij], 0.0, &(T2->matrix[0][ij][0]), nvir);
772
/* Write the updated MO-basis T2's to disk */
773
dpd_buf4_mat_irrep_wrt(T2, 0);
774
dpd_buf4_mat_irrep_close(T2, 0);
777
void local_filter_U2(dpdbuf4 *T2, double lambda)
779
int ij, i, j, a, b, ab;
781
int *pairdom_len, *pairdom_nrlen, *weak_pairs;
782
double ***V, ***W, *eps_occ, **eps_vir;
783
double **X1, **X2, **T2tilde, **T2bar;
790
eps_occ = local.eps_occ;
791
eps_vir = local.eps_vir;
792
pairdom_len = local.pairdom_len;
793
pairdom_nrlen = local.pairdom_nrlen;
794
weak_pairs = local.weak_pairs;
796
/* Grab the MO-basis T2's */
797
dpd_buf4_mat_irrep_init(T2, 0);
798
dpd_buf4_mat_irrep_rd(T2, 0);
800
X1 = block_matrix(nso,nvir);
801
X2 = block_matrix(nvir,nso);
802
T2tilde = block_matrix(nso,nso);
803
T2bar = block_matrix(nvir, nvir);
804
for(i=0,ij=0; i < nocc; i++) {
805
for(j=0; j < nocc; j++,ij++) {
807
/* Transform the virtuals to the redundant projected virtual basis */
808
C_DGEMM('t', 'n', pairdom_len[ij], nvir, nvir, 1.0, &(V[ij][0][0]), pairdom_len[ij],
809
&(T2->matrix[0][ij][0]), nvir, 0.0, &(X1[0][0]), nvir);
810
C_DGEMM('n', 'n', pairdom_len[ij], pairdom_len[ij], nvir, 1.0, &(X1[0][0]), nvir,
811
&(V[ij][0][0]), pairdom_len[ij], 0.0, &(T2tilde[0][0]), nso);
813
/* Transform the virtuals to the non-redundant virtual basis */
814
C_DGEMM('t', 'n', pairdom_nrlen[ij], pairdom_len[ij], pairdom_len[ij], 1.0,
815
&(W[ij][0][0]), pairdom_nrlen[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
816
C_DGEMM('n', 'n', pairdom_nrlen[ij], pairdom_nrlen[ij], pairdom_len[ij], 1.0,
817
&(X2[0][0]), nso, &(W[ij][0][0]), pairdom_nrlen[ij], 0.0, &(T2bar[0][0]), nvir);
819
/* Divide the new amplitudes by the denominators */
820
for(a=0; a < pairdom_nrlen[ij]; a++) {
821
for(b=0; b < pairdom_nrlen[ij]; b++) {
822
T2bar[a][b] /= (eps_occ[i] + eps_occ[j] - eps_vir[ij][a] - eps_vir[ij][b] + lambda);
826
/* Transform the new T2's to the redundant virtual basis */
827
C_DGEMM('n', 'n', pairdom_len[ij], pairdom_nrlen[ij], pairdom_nrlen[ij], 1.0,
828
&(W[ij][0][0]), pairdom_nrlen[ij], &(T2bar[0][0]), nvir, 0.0, &(X1[0][0]), nvir);
829
C_DGEMM('n','t', pairdom_len[ij], pairdom_len[ij], pairdom_nrlen[ij], 1.0,
830
&(X1[0][0]), nvir, &(W[ij][0][0]), pairdom_nrlen[ij], 0.0, &(T2tilde[0][0]), nso);
832
/* Transform the new T2's to the MO basis */
833
C_DGEMM('n', 'n', nvir, pairdom_len[ij], pairdom_len[ij], 1.0,
834
&(V[ij][0][0]), pairdom_len[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
835
C_DGEMM('n', 't', nvir, nvir, pairdom_len[ij], 1.0, &(X2[0][0]), nso,
836
&(V[ij][0][0]), pairdom_len[ij], 0.0, &(T2->matrix[0][ij][0]), nvir);
846
/* Write the updated MO-basis T2's to disk */
847
dpd_buf4_mat_irrep_wrt(T2, 0);
848
dpd_buf4_mat_irrep_close(T2, 0);
852
}} // namespace psi::cis
added
removed
src/bin/cis/local.cc
