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>
26
** local_init(): Set up parameters of local excitation domains.
28
** The orbital domains constructed here are based on those described
29
** in Broughton and Pulay, J. Comp. Chem. 14, 736-740 (1993). The
30
** localization of the occupied orbitals is done elsewhere (see the
31
** program "local"). Pair domains are defined as the union of pairs
32
** of single occupied orbital domains. "Weak pairs", which are
33
** defined as pair domains whose individual occupied orbital domains
34
** have no atoms in common, are identified (cf. int *weak_pairs).
39
void domain_print(int, int, int *, int **, double *);
40
void transpert(char *);
41
void sort_pert(char *, double **, double **, double **, int, int, int);
42
void build_F_RHF(double);
43
void build_B_RHF(double);
46
void local_polar(char*, int **, int *, int, int *, int *);
47
void local_magnetic(char*, int **, int *, int, int *, int *);
51
int i, j, k, ij, stat, a, b, l, I, L;
52
int nmo, nso, nocc, nocc_all, nvir, noei, nirreps, nfzc;
53
double domain_tot, domain_ave;
54
double **C; /* AO -> localized MO transformation matrix */
55
double **Ci; /* localized MO -> AO transformation matrix */
56
double **D; /* 1/2 SCF closed-shell density matrix (AO) */
57
double **Rt, **Rt_full; /* Projected, redundant virtual transform (R-tilde) */
58
double **S; /* AO overlap */
59
double **St; /* Projected virtual overlap */
60
double **Xt; /* Projected, non-redundant virtual transform (X-tilde) */
61
double ***V; /* MO -> projected, redundant virtual transform */
62
double **Fmo;/* MO basis Fock matrix */
63
double **F; /* AO basis Fock matrix */
64
double **Ft; /* Projected, redundant virtual Fock matrix */
65
double **Fbar; /* Projected, non-redundant virtual Fock matrix */
66
double ***W; /* Transformation matrix from tilde -> bar for each ij pair*/
67
double *eps_occ; /* occupied orbital energies for local denominators */
68
double **eps_vir; /* virtual orbital energies for local denominators */
70
double *evals, **evecs;
71
double *eps_all; /* All MO energies */
74
int natom, atom, am, offset, nshell, shell_length, next_atom;
75
int row, col, max, m, errcod, cnt;
76
int *rank, *boolean, *ipiv;
77
int *l_length, *aostart, *aostop, *ao2atom;
79
int **domain, *domain_len, **pairdomain, *pairdom_len, *pairdom_nrlen;
80
int **domain_bp, *domain_len_bp;
82
double *fR, *charge, *SR, *Z, tmp, *ss;
84
int print_test, num_entries, entry_len, orbital;
85
int t1_length, t2_length, puream, weak;
91
chkpt_init(PSIO_OPEN_OLD);
93
natom = chkpt_rd_natom();
94
nshell = chkpt_rd_nshell();
95
eps_all = chkpt_rd_evals();
96
stype = chkpt_rd_stype();
97
snuc = chkpt_rd_snuc();
102
/* C1 symmetry only */
103
nirreps = moinfo.nirreps;
105
fprintf(outfile, "\nError: localization must use C1 symmetry.\n");
106
exit(PSI_RETURN_FAILURE);
110
nmo = moinfo.nmo; /* should be the same as nso */
112
fprintf(outfile, "\nError: NMO != NSO! %d != %d\n", nmo, nso);
113
exit(PSI_RETURN_FAILURE);
116
nocc = moinfo.occpi[0]; /* active doubly occupied orbitals */
117
nfzc = moinfo.frdocc[0]; /* frozen doubly occupied orbitals */
118
nocc_all = nocc + nfzc; /* all doubly occupied orbitals */
119
nvir = moinfo.virtpi[0]; /* active virtual orbitals */
126
/* A couple of scratch arrays */
127
X = block_matrix(nso, nso);
128
Y = block_matrix(nso, nso);
131
Ci = block_matrix(nso, nso);
132
for(i=0; i < nso; i++)
133
for(j=0; j < nso; j++)
136
invert_matrix(C, Ci, nso, outfile);
138
for(i=0; i < nso; i++)
139
for(j=0; j < nso; j++)
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
domain_bp = init_int_matrix(nocc, natom);
208
domain_len_bp = 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];
292
/* Augment the domain if necessary */
293
if(fabs(fR[i-nfzc]) > local.cutoff) {
294
domain[i-nfzc][rank[next_atom++]] = 1;
295
domain_len[i-nfzc]++;
300
for(i=0; i<nocc; i++) {
301
domain_len_bp[i] = domain_len[i];
302
for(k=0; k<natom; k++)
303
domain_bp[i][k] = domain[i][k];
305
/* Print the orbital domains */
306
fprintf(outfile, "\n ****** Boughton-Pulay Occupied Orbital Domains ******\n");
307
domain_print(nocc, natom, domain_len_bp, domain_bp, fR);
309
/* Identify and/or remove weak pairs -- using Bougton-Pulay domains */
310
weak_pairs = init_int_array(nocc*nocc);
311
if(!strcmp(local.pairdef,"BP")) {
312
fprintf(outfile, "\n");
313
for(i=0,ij=0; i < nocc; i++)
314
for(j=0; j < nocc; j++,ij++) {
316
for(k=0; k < natom; k++)
317
if(domain[i][k] && domain[j][k]) weak = 0;
319
if(weak && strcmp(local.weakp,"NONE")) {
322
if(!strcmp(local.weakp,"MP2"))
323
fprintf(outfile, "\tPair %d %d [%d] is weak and will be treated with MP2.\n", i, j, ij);
324
else if(!strcmp(local.weakp,"NEGLECT")) {
325
fprintf(outfile, "\tPair %d %d = [%d] is weak and will be deleted.\n", i, j, ij);
328
else weak_pairs[ij] = 0;
332
/* If this is a response calculation, augment domains using polarized orbitals */
333
fprintf(outfile, "\n");
334
if(local.domain_polar) {
335
fprintf(outfile, "\tGenerating electric-field CPHF solutions for local-CC.\n");
338
sort_pert("Mu", moinfo.MUX, moinfo.MUY, moinfo.MUZ,
339
moinfo.irrep_x, moinfo.irrep_y, moinfo.irrep_z);
341
if(local.domain_sep) {
345
local_polar("X", domain, domain_len, natom, aostart, aostop);
346
fprintf(outfile, "\nMu_X (%lf)\n", 0);
347
domain_print(nocc, natom, domain_len, domain, fR);
348
for(i=0; i<nocc; i++) {
349
domain_len[i] = domain_len_bp[i];
350
for(k=0; k<natom; k++)
351
domain[i][k] = domain_bp[i][k];
355
local_polar("Y", domain, domain_len, natom, aostart, aostop);
356
fprintf(outfile, "\nMu_Y (%lf)\n", 0);
357
domain_print(nocc, natom, domain_len, domain, fR);
358
for(i=0; i<nocc; i++) {
359
domain_len[i] = domain_len_bp[i];
360
for(k=0; k<natom; k++)
361
domain[i][k] = domain_bp[i][k];
365
local_polar("Z", domain, domain_len, natom, aostart, aostop);
366
fprintf(outfile, "\nMu_Z (%lf)\n", 0);
367
domain_print(nocc, natom, domain_len, domain, fR);
368
for(i=0; i<nocc; i++) {
369
domain_len[i] = domain_len_bp[i];
370
for(k=0; k<natom; k++)
371
domain[i][k] = domain_bp[i][k];
375
build_F_RHF(params.omega[0]);
377
local_polar("X", domain, domain_len, natom, aostart, aostop);
378
fprintf(outfile, "\nMu_X (%lf)\n", params.omega[0]);
379
domain_print(nocc, natom, domain_len, domain, fR);
380
for(i=0; i<nocc; i++) {
381
domain_len[i] = domain_len_bp[i];
382
for(k=0; k<natom; k++)
383
domain[i][k] = domain_bp[i][k];
387
local_polar("Y", domain, domain_len, natom,
389
fprintf(outfile, "\nMu_Y (%lf)\n", params.omega[0]);
390
domain_print(nocc, natom, domain_len, domain, fR);
391
for(i=0; i<nocc; i++) {
392
domain_len[i] = domain_len_bp[i];
393
for(k=0; k<natom; k++)
394
domain[i][k] = domain_bp[i][k];
398
local_polar("Z", domain, domain_len, natom, aostart, aostop);
399
fprintf(outfile, "\nMu_Z (%lf)\n", params.omega[0]);
400
domain_print(nocc, natom, domain_len, domain, fR);
401
for(i=0; i<nocc; i++) {
402
domain_len[i] = domain_len_bp[i];
403
for(k=0; k<natom; k++)
404
domain[i][k] = domain_bp[i][k];
408
build_F_RHF(-params.omega[0]);
410
local_polar("X", domain, domain_len, natom, aostart, aostop);
411
fprintf(outfile, "\nMu_X (%lf)\n", -params.omega[0]);
412
domain_print(nocc, natom, domain_len, domain, fR);
413
for(i=0; i<nocc; i++) {
414
domain_len[i] = domain_len_bp[i];
415
for(k=0; k<natom; k++)
416
domain[i][k] = domain_bp[i][k];
420
local_polar("Y", domain, domain_len, natom,
422
fprintf(outfile, "\nMu_Y (%lf)\n", -params.omega[0]);
423
domain_print(nocc, natom, domain_len, domain, fR);
424
for(i=0; i<nocc; i++) {
425
domain_len[i] = domain_len_bp[i];
426
for(k=0; k<natom; k++)
427
domain[i][k] = domain_bp[i][k];
431
local_polar("Z", domain, domain_len, natom, aostart, aostop);
432
fprintf(outfile, "\nMu_Z (%lf)\n", -params.omega[0]);
433
domain_print(nocc, natom, domain_len, domain, fR);
434
for(i=0; i<nocc; i++) {
435
domain_len[i] = domain_len_bp[i];
436
for(k=0; k<natom; k++)
437
domain[i][k] = domain_bp[i][k];
443
local_polar("X", domain, domain_len, natom, aostart, aostop);
445
local_polar("Y", domain, domain_len, natom, aostart, aostop);
447
local_polar("Z", domain, domain_len, natom, aostart, aostop);
451
if(local.domain_mag) {
452
fprintf(outfile, "\tGenerating magnetic-field CPHF solutions for local-CC.\n");
455
sort_pert("L", moinfo.LX, moinfo.LY, moinfo.LZ,
456
moinfo.irrep_x, moinfo.irrep_y, moinfo.irrep_z);
458
if(local.domain_sep) {
462
local_magnetic("X", domain, domain_len, natom, aostart, aostop);
463
fprintf(outfile, "\nL_X (%lf)\n", 0);
464
domain_print(nocc, natom, domain_len, domain, fR);
465
for(i=0; i<nocc; i++) {
466
domain_len[i] = domain_len_bp[i];
467
for(k=0; k<natom; k++)
468
domain[i][k] = domain_bp[i][k];
472
local_magnetic("Y", domain, domain_len, natom,
474
fprintf(outfile, "\nL_Y (%lf)\n", 0);
475
domain_print(nocc, natom, domain_len, domain, fR);
476
for(i=0; i<nocc; i++) {
477
domain_len[i] = domain_len_bp[i];
478
for(k=0; k<natom; k++)
479
domain[i][k] = domain_bp[i][k];
483
local_magnetic("Z", domain, domain_len, natom, aostart, aostop);
484
fprintf(outfile, "\nL_Z (%lf)\n", 0);
485
domain_print(nocc, natom, domain_len, domain, fR);
486
for(i=0; i<nocc; i++) {
487
domain_len[i] = domain_len_bp[i];
488
for(k=0; k<natom; k++)
489
domain[i][k] = domain_bp[i][k];
493
build_B_RHF(params.omega[0]);
495
local_magnetic("X", domain, domain_len, natom, aostart, aostop);
496
fprintf(outfile, "\nL_X (%lf)\n", params.omega[0]);
497
domain_print(nocc, natom, domain_len, domain, fR);
498
for(i=0; i<nocc; i++) {
499
domain_len[i] = domain_len_bp[i];
500
for(k=0; k<natom; k++)
501
domain[i][k] = domain_bp[i][k];
505
local_magnetic("Y", domain, domain_len, natom, aostart, aostop);
506
fprintf(outfile, "\nL_Y (%lf)\n", params.omega[0]);
507
domain_print(nocc, natom, domain_len, domain, fR);
508
for(i=0; i<nocc; i++) {
509
domain_len[i] = domain_len_bp[i];
510
for(k=0; k<natom; k++)
511
domain[i][k] = domain_bp[i][k];
515
local_magnetic("Z", domain, domain_len, natom, aostart, aostop);
516
fprintf(outfile, "\nL_Z (%lf)\n", params.omega[0]);
517
domain_print(nocc, natom, domain_len, domain, fR);
518
for(i=0; i<nocc; i++) {
519
domain_len[i] = domain_len_bp[i];
520
for(k=0; k<natom; k++)
521
domain[i][k] = domain_bp[i][k];
525
build_B_RHF(-params.omega[0]);
527
local_magnetic("X", domain, domain_len, natom, aostart, aostop);
528
fprintf(outfile, "\nL_X (%lf)\n", -params.omega[0]);
529
domain_print(nocc, natom, domain_len, domain, fR);
530
for(i=0; i<nocc; i++) {
531
domain_len[i] = domain_len_bp[i];
532
for(k=0; k<natom; k++)
533
domain[i][k] = domain_bp[i][k];
537
local_magnetic("Y", domain, domain_len, natom, aostart, aostop);
538
fprintf(outfile, "\nL_Y (%lf)\n", -params.omega[0]);
539
domain_print(nocc, natom, domain_len, domain, fR);
540
for(i=0; i<nocc; i++) {
541
domain_len[i] = domain_len_bp[i];
542
for(k=0; k<natom; k++)
543
domain[i][k] = domain_bp[i][k];
547
local_magnetic("Z", domain, domain_len, natom, aostart, aostop);
548
fprintf(outfile, "\nL_Z (%lf)\n", -params.omega[0]);
549
domain_print(nocc, natom, domain_len, domain, fR);
550
for(i=0; i<nocc; i++) {
551
domain_len[i] = domain_len_bp[i];
552
for(k=0; k<natom; k++)
553
domain[i][k] = domain_bp[i][k];
559
local_magnetic("X", domain, domain_len, natom, aostart, aostop);
561
local_magnetic("Y", domain, domain_len, natom, aostart, aostop);
563
local_magnetic("Z", domain, domain_len, natom, aostart, aostop);
567
/* Allow user input of selected domains */
568
if(ip_exist("DOMAINS",0)) {
569
ip_count("DOMAINS", &num_entries,0);
570
for(i=0; i < num_entries; i++) {
571
ip_count("DOMAINS", &entry_len, 1, i);
572
ip_data("DOMAINS", "%d", &orbital, 2, i, 0);
574
/* Clear out the current domain for this orbital */
575
for(j=0; j < natom; j++) domain[orbital][j] = 0;
576
domain_len[orbital] = 0;
578
for(j=1; j < entry_len; j++) {
579
errcod = ip_data("DOMAINS","%d", &atom,2,i,j);
580
domain[orbital][atom] = 1;
581
domain_len[orbital]++;
586
/* Recheck Completeness */
587
for(i=nfzc; i < nocc_all; i++) {
589
/* Build the orbital's domain starting in order of decreasing charge contribution */
590
for(j=0; j < nso; j++) {
592
for(k=0; k < nso; k++)
593
SR[j] += S[j][k] * C[k][i];
596
for(j=0,row=0; j < natom; j++) {
597
if(domain[i-nfzc][j]) {
598
for(k=aostart[j]; k <= aostop[j]; k++,row++) {
602
for(l=0,col=0; l < natom; l++) {
603
if(domain[i-nfzc][l]) {
605
for(m=aostart[l]; m <= aostop[l]; m++,col++)
606
X[row][col] = S[k][m];
615
errcod = C_DGESV(row, 1, &(X[0][0]), nso, &(ipiv[0]), &(Z[0]), nso);
617
fprintf(outfile, "\nError in DGESV return in orbital domain construction.\n");
618
exit(PSI_RETURN_FAILURE);
622
for(j=0,row=0; j < natom; j++) {
623
if(domain[i-nfzc][j]) {
624
for(k=aostart[j]; k <= aostop[j]; k++,row++) {
625
for(l=0; l < nso; l++) fR[i-nfzc] -= Z[row] * S[k][l] * C[l][i];
633
/* Print the orbital domains */
634
fprintf(outfile, "\n ****** Final Occupied Orbital Domains ******\n");
635
if(!local.domain_sep)
636
domain_print(nocc, natom, domain_len, domain, fR);
638
/* Build the pair domains */
639
pairdomain = init_int_matrix(nocc*nocc,natom);
640
pairdom_len = init_int_array(nocc*nocc);
641
for(i=0,ij=0; i < nocc; i++)
642
for(j=0; j < nocc; j++,ij++)
643
for(k=0; k < natom; k++) {
644
if(domain[i][k] || domain[j][k]) {
645
pairdomain[ij][k] = 1;
646
pairdom_len[ij] += aostop[k] - aostart[k] + 1;
650
/* Identify and/or remove weak pairs -- for CPHF "response" domains */
651
if(local.domain_polar || local.domain_mag || ip_exist("DOMAINS",0)) {
652
fprintf(outfile, "\n");
653
for(i=0,ij=0; i < nocc; i++)
654
for(j=0; j < nocc; j++,ij++) {
656
for(k=0; k < natom; k++)
657
if(domain[i][k] && domain[j][k]) weak = 0;
659
if(weak && strcmp(local.weakp,"NONE")) {
662
if(!strcmp(local.weakp,"MP2"))
663
fprintf(outfile, "\tPair %d %d [%d] is weak and will be treated with MP2.\n", i, j, ij);
664
else if(!strcmp(local.weakp,"NEGLECT")) {
665
fprintf(outfile, "\tPair %d %d = [%d] is weak and will be deleted.\n", i, j, ij);
668
else weak_pairs[ij] = 0;
672
/* Compute the total number of singles and doubles */
673
/* replacement code from TDC on 11-5-02 */
674
t1_length = t2_length = 0;
675
for(i=0,ij=0; i < nocc; i++) {
676
for(k=0; k < natom; k++) {
678
for(a=aostart[k]; a <= aostop[k]; a++) t1_length++;
680
for(j=0; j < nocc; j++,ij++) {
681
for(k=0; k < natom; k++) {
682
for(l=0; l < natom; l++) {
683
if(pairdomain[ij][k] && pairdomain[ij][l] && !weak_pairs[ij]) {
684
for(a=aostart[k]; a <= aostop[k]; a++)
685
for(b=aostart[l]; b <= aostop[l]; b++)
693
/* Print excitation space reduction info */
694
fprintf(outfile, "\n\tT1 Length = %d (local), %d (canonical)\n",
695
t1_length, nocc*nvir);
696
fprintf(outfile, "\tT2 Length = %d (local), %d (canonical)\n\n",
697
t2_length, nocc*nocc*nvir*nvir);
700
local.domain = domain;
701
local.domain_len = domain_len;
702
local.pairdomain = pairdomain;
703
local.pairdom_len = pairdom_len;
704
local.weak_pairs = weak_pairs;
705
local.aostart = aostart;
706
local.aostop = aostop;
708
free_int_matrix(domain_bp, nocc);
722
ip_boolean("DOMAIN_PRINT",&(print_test),0);
724
fprintf(outfile, "Printing of orbital domains requested...exiting.\n\n");
725
exit(PSI_RETURN_FAILURE);
728
/************* Orbital Domains Complete ***************/
730
/* Compute the complete virtual space projector */
731
Rt_full = block_matrix(nso,nso);
732
for(i=0; i < nso; i++) Rt_full[i][i] = 1.0;
734
C_DGEMM('n','n',nso,nso,nso,-1.0,&(D[0][0]),nso,&(S[0][0]),nso,
735
1.0,&(Rt_full[0][0]),nso);
738
fprintf(outfile, "\n\tVirtual-Space Projector (R-tilde):\n");
739
print_mat(Rt_full, nso, nso, stdout);
742
/* Compute the norm of each PAO */
743
for(i=0; i < nso; i++) {
745
for(j=0; j < nso; j++) {
746
norm += Rt_full[j][i] * Rt_full[j][i];
749
if(norm < 0.2 && strcmp(local.freeze_core,"FALSE")) {
750
fprintf(outfile, "\tNorm of orbital %4d = %20.12f...deleteing\n", i, norm);
751
for(j=0; j < nso; j++) Rt_full[j][i] = 0.0;
754
fprintf(outfile, "\n");
757
/* Grab the MO-basis Fock matrix */
758
Fmo = block_matrix(nso, nso);
759
for(i=0; i < nfzc; i++) Fmo[i][i] = eps_all[i];
760
dpd_file2_init(&fock, CC_OEI, 0, 0, 0, "fIJ");
761
dpd_file2_mat_init(&fock);
762
dpd_file2_mat_rd(&fock);
763
for(i=0; i < nocc; i++)
764
for(j=0; j < nocc; j++)
765
Fmo[i+nfzc][j+nfzc] = fock.matrix[0][i][j];
766
dpd_file2_mat_close(&fock);
767
dpd_file2_close(&fock);
769
dpd_file2_init(&fock, CC_OEI, 0, 1, 1, "fAB");
770
dpd_file2_mat_init(&fock);
771
dpd_file2_mat_rd(&fock);
772
for(i=0; i < nvir; i++)
773
for(j=0; j < nvir; j++)
774
Fmo[i+nfzc+nocc][j+nfzc+nocc] = fock.matrix[0][i][j];
775
dpd_file2_mat_close(&fock);
776
dpd_file2_close(&fock);
779
fprintf(outfile, "\n\tMO Basis Fock matrix:\n");
780
print_mat(Fmo, nso, nso, outfile);
783
/* Build the AO-basis Fock matrix */
784
F = block_matrix(nso,nso);
785
C_DGEMM('t','n',nso,nso,nso,1.0,&(Ci[0][0]),nso,&(Fmo[0][0]),nso,
787
C_DGEMM('n','n',nso,nso,nso,1.0,&(X[0][0]),nso,&(Ci[0][0]),nso,
790
/* Build the occupied orbital energy list */
791
eps_occ = init_array(nocc);
792
for(i=0;i < nocc; i++) eps_occ[i] = Fmo[i+nfzc][i+nfzc];
795
fprintf(outfile, "\n\tAO-Basis Fock Matrix:\n");
796
print_mat(Fmo, nso, nso, outfile);
799
/* Compute R^+ S for virtual orbitals */
800
RS = block_matrix(nvir,nso);
801
for(a=0; a < nvir; a++)
802
for(i=0; i < nso; i++)
803
X[i][a] = C[i][a+nocc_all];
804
C_DGEMM('t','n',nvir,nso,nso,1.0,&(X[0][0]),nso,&(S[0][0]),nso,0.0,&(RS[0][0]),nso);
806
/* Build the virtual metric and W transforms for each pair domain */
807
Rt = block_matrix(nso, nso);
808
W = (double ***) malloc(nocc * nocc * sizeof(double **));
809
V = (double ***) malloc(nocc * nocc * sizeof(double **));
810
eps_vir = (double **) malloc(nocc * nocc * sizeof(double *));
811
pairdom_nrlen = init_int_array(nocc * nocc); /* dimension of non-redundant basis */
814
for(ij=0; ij < nocc * nocc; ij++) {
816
zero_mat(Rt, nso, nso);
818
/* Build the virtual space projector for this pair */
819
for(k=0,L=0; k < natom; k++) {
820
if(pairdomain[ij][k]) {
821
for(l=aostart[k]; l <= aostop[k]; l++,L++) {
822
for(m=0; m < nso; m++) {
823
Rt[m][L] = Rt_full[m][l];
829
/* Compute the MO -> projected virtual transformation matrix */
830
V[ij] = block_matrix(nvir,pairdom_len[ij]);
831
C_DGEMM('n','n',nvir,pairdom_len[ij],nso,1.0,&(RS[0][0]),nso,&(Rt[0][0]),nso,0.0,
832
&(V[ij][0][0]),pairdom_len[ij]);
835
fprintf(outfile, "\nV[%d]:\n", ij);
836
fprintf(outfile, "======\n");
837
print_mat(V[ij], nvir, pairdom_len[ij], outfile);
840
/* Virtual space metric */
841
St = block_matrix(pairdom_len[ij],pairdom_len[ij]);
842
C_DGEMM('n','n',nso,pairdom_len[ij],nso,1.0,&(S[0][0]),nso,&(Rt[0][0]),nso,
844
C_DGEMM('t','n',pairdom_len[ij],pairdom_len[ij],nso,1.0,&(Rt[0][0]),nso,&(X[0][0]),nso,
845
0.0,&(St[0][0]),pairdom_len[ij]);
848
fprintf(outfile, "\n\tVirtual-Space Metric (S-tilde) for ij = %d:\n", ij);
849
print_mat(St, pairdom_len[ij], pairdom_len[ij], outfile);
852
/* Diagonalize metric */
853
evals = init_array(pairdom_len[ij]);
854
evecs = block_matrix(pairdom_len[ij],pairdom_len[ij]);
855
sq_rsp(pairdom_len[ij],pairdom_len[ij],St,evals,1,evecs,1e-12);
857
/* Count the number of zero eigenvalues */
858
for(i=0,cnt=0; i < pairdom_len[ij]; i++) if(evals[i] <= 1e-6) cnt++;
860
pairdom_nrlen[ij] = pairdom_len[ij]-cnt;
863
fprintf(outfile, "\n\tS-tilde eigenvalues for ij = %d:\n", ij);
864
for(i=0; i < pairdom_len[ij]; i++) fprintf(outfile, "\t%d %20.12f\n", i, evals[i]);
866
fprintf(outfile, "\n\tS-tilde eigenvectors for ij = %d:\n", ij);
867
print_mat(evecs,pairdom_len[ij],pairdom_len[ij],outfile);
870
/* Build the projected, non-redundant transform (X-tilde) */
871
Xt = block_matrix(pairdom_len[ij],pairdom_nrlen[ij]);
872
for(i=0,I=0; i < pairdom_len[ij]; i++) {
873
if(evals[i] > 1e-6) {
874
for(j=0; j < pairdom_len[ij]; j++)
875
Xt[j][I] = evecs[j][i]/sqrt(evals[i]);
883
fprintf(outfile, "\n\tTransform to non-redundant, projected virtuals (X-tilde) for ij = %d:\n", ij);
884
print_mat(Xt, pairdom_len[ij], pairdom_nrlen[ij], outfile);
890
/* Build the projected (redundant) virtual Fock matrix */
891
Ft = block_matrix(pairdom_len[ij], pairdom_len[ij]);
892
C_DGEMM('t','n',pairdom_len[ij],nso,nso,1.0,&(Rt[0][0]),nso,&(F[0][0]),nso,
894
C_DGEMM('n','n',pairdom_len[ij],pairdom_len[ij],nso,1.0,&(X[0][0]),nso,&(Rt[0][0]),nso,
895
0.0,&(Ft[0][0]),pairdom_len[ij]);
897
/* Project the Fock matrix into the non-redundant virtual space */
898
Fbar = block_matrix(pairdom_nrlen[ij],pairdom_nrlen[ij]);
899
C_DGEMM('t','n',pairdom_nrlen[ij],pairdom_len[ij],pairdom_len[ij],1.0,
900
&(Xt[0][0]),pairdom_nrlen[ij],&(Ft[0][0]),pairdom_len[ij],0.0,&(X[0][0]),nso);
901
C_DGEMM('n','n',pairdom_nrlen[ij],pairdom_nrlen[ij],pairdom_len[ij],1.0,
902
&(X[0][0]),nso,&(Xt[0][0]),pairdom_nrlen[ij],0.0,&(Fbar[0][0]),pairdom_nrlen[ij]);
905
fprintf(outfile, "\n\tFbar matrix for ij = %d:\n", ij);
906
print_mat(Fbar,pairdom_nrlen[ij],pairdom_nrlen[ij],outfile);
909
/* Diagonalize Fbar */
910
evals = init_array(pairdom_nrlen[ij]);
911
evecs = block_matrix(pairdom_nrlen[ij],pairdom_nrlen[ij]);
912
sq_rsp(pairdom_nrlen[ij],pairdom_nrlen[ij],Fbar,evals,1,evecs,1e-12);
915
fprintf(stdout, "\n\tFbar eigenvectors for ij = %d:\n", ij);
916
print_mat(evecs,pairdom_nrlen[ij],pairdom_nrlen[ij],outfile);
919
/* Finally, build the W matrix */
920
W[ij] = block_matrix(pairdom_len[ij],pairdom_nrlen[ij]);
921
C_DGEMM('n','n',pairdom_len[ij],pairdom_nrlen[ij],pairdom_nrlen[ij],1.0,
922
&(Xt[0][0]),pairdom_nrlen[ij],&(evecs[0][0]),pairdom_nrlen[ij],
923
0.0,&(W[ij][0][0]),pairdom_nrlen[ij]);
926
fprintf(outfile, "\n\tW Transformation Matrix for ij = %d:\n", ij);
927
print_mat(W[ij],pairdom_len[ij],pairdom_nrlen[ij],outfile);
930
/* build the orbital energy list */
931
eps_vir[ij] = init_array(pairdom_nrlen[ij]);
932
for(i=0; i < pairdom_nrlen[ij]; i++)
933
eps_vir[ij][i] = evals[i]; /* virtual orbital energies */
936
fprintf(outfile, "\n\tVirtual orbital Energies for ij = %d:\n", ij);
937
for(i=0; i < pairdom_nrlen[ij]; i++)
938
fprintf(outfile, "%d %20.12f\n", i, eps_vir[ij][i]);
965
local.eps_occ = eps_occ;
966
local.eps_vir = eps_vir;
967
local.pairdom_nrlen = pairdom_nrlen;
969
local.weak_pair_energy = 0.0;
975
void local_done(void)
977
int i, ij, h, nocc, nvir, natom;
984
psio_write_entry(CC_INFO, "Local Cutoff", (char *) &local.cutoff,
986
psio_write_entry(CC_INFO, "Local Domain Length", (char *) local.domain_len,
988
psio_write_entry(CC_INFO, "Local Pair Domain Length", (char *) local.pairdom_len,
989
nocc*nocc*sizeof(int));
990
psio_write_entry(CC_INFO, "Local Pair Domain NR Length", (char *) local.pairdom_nrlen,
991
nocc*nocc*sizeof(int));
992
psio_write_entry(CC_INFO, "Local Weak Pairs", (char *) local.weak_pairs,
993
nocc*nocc*sizeof(int));
994
psio_write_entry(CC_INFO, "Local Occupied Orbital Energies", (char *) local.eps_occ,
995
nocc*sizeof(double));
998
for(i=0; i<nocc; i++)
999
psio_write(CC_INFO, "Local Domains", (char *) local.domain[i],
1000
natom*sizeof(int), next, &next);
1002
for(ij=0; ij<nocc*nocc; ij++)
1003
psio_write(CC_INFO, "Local Pair Domains", (char *) local.pairdomain[ij],
1004
natom*sizeof(int), next, &next);
1006
for(ij=0; ij < nocc*nocc; ij++)
1007
psio_write(CC_INFO, "Local Virtual Orbital Energies", (char *) local.eps_vir[ij],
1008
local.pairdom_nrlen[ij]*sizeof(double), next, &next);
1010
for(ij=0; ij < nocc*nocc; ij++)
1011
psio_write(CC_INFO, "Local Transformation Matrix (W)", (char *) local.W[ij][0],
1012
local.pairdom_len[ij]*local.pairdom_nrlen[ij]*sizeof(double), next, &next);
1014
for(ij=0; ij < nocc*nocc; ij++)
1015
psio_write(CC_INFO, "Local Residual Vector (V)", (char *) local.V[ij][0],
1016
nvir*local.pairdom_len[ij]*sizeof(double), next, &next);
1018
if(params.ref == 0 || params.ref == 1) {
1019
for(h=0; h < moinfo.nirreps; h++)
1020
if(moinfo.orbspi[h] && moinfo.virtpi[h]) free_block(moinfo.C[h]);
1024
free_int_matrix(local.pairdomain, nocc*nocc);
1025
free_int_matrix(local.domain, nocc);
1026
for(i=0; i < nocc*nocc; i++) {
1027
free_block(local.W[i]);
1028
free_block(local.V[i]);
1029
free(local.eps_vir[i]);
1033
free(local.eps_vir);
1035
free(local.aostart);
1038
free(local.eps_occ);
1039
free(local.domain_len);
1040
free(local.pairdom_len);
1041
free(local.pairdom_nrlen);
1042
free(local.weak_pairs);
added
removed
src/bin/ccsort/local.c
