4
#include <libciomr/libciomr.h>
5
#include <libint/libint.h>
6
#include <libchkpt/chkpt.h>
7
#include <libiwl/iwl.h>
11
FILE *infile, *outfile;
12
char *psi_file_prefix;
14
void init_io(int argc, char *argv[]);
18
int main(int argc, char *argv[])
20
int iter, s, t, A, k, l, m, p, q, inew, iold, max, max_col, phase_ok, phase_chk;
21
int i, j, ij, am, atom, shell_length, offset, stat;
22
int nirreps, nao, nmo, nso, natom, nshell, noei, nocc, errcod, nfzc;
23
int *stype, *snuc, *aostart, *aostop, *ao2atom, *l_length;
24
int *clsdpi, *openpi, *orbspi, *dummy, *order, *frdocc;
25
double *ss, **S, **scf, **u, **Ctmp, **C, *evals, **MO_S, **scf_old, **X;
27
int *orb_order, *orb_boolean, puream;
28
double P, PiiA, Pst, Pss, Ptt, Ast, Bst, AB;
29
double Uss, Utt, Ust, Uts, Cks, Ckt, **U, **V, **VV, **F;
30
double cos4a, alpha, alphamax, alphalast, conv, norm;
38
chkpt_init(PSIO_OPEN_OLD);
42
natom = chkpt_rd_natom();
43
nshell = chkpt_rd_nshell();
44
stype = chkpt_rd_stype();
45
snuc = chkpt_rd_snuc();
46
u = chkpt_rd_usotao();
47
nirreps = chkpt_rd_nirreps();
48
clsdpi = chkpt_rd_clsdpi();
49
openpi = chkpt_rd_openpi();
50
orbspi = chkpt_rd_orbspi();
52
evals = chkpt_rd_evals();
55
/* A couple of error traps */
57
fprintf(outfile, "\n\tError: localization is only valid in C1 symmetry!\n");
58
exit(PSI_RETURN_FAILURE);
61
fprintf(outfile, "\n\tError: localization available for closed-shells only!\n");
62
exit(PSI_RETURN_FAILURE);
65
/* Frozen orbital info */
66
frdocc = get_frzcpi();
70
/* Compute the length of each AM block */
71
ip_boolean("PUREAM", &(puream), 0);
72
l_length = init_int_array(LIBINT_MAX_AM);
74
for(l=1; l < (LIBINT_MAX_AM); l++) {
75
if(puream) l_length[l] = 2 * l + 1;
76
else l_length[l] = l_length[l-1] + l + 1;
79
/* Set up the atom->AO and AO->atom lookup arrays */
80
aostart = init_int_array(natom);
81
aostop = init_int_array(natom);
82
for(i=0,atom=-1,offset=0; i < nshell; i++) {
84
shell_length = l_length[am];
86
if(atom != snuc[i]-1) {
87
if(atom != -1) aostop[atom] = offset-1;
89
aostart[atom] = offset;
92
offset += shell_length;
94
aostop[atom] = offset-1;
96
ao2atom = init_int_array(nso);
97
for(i=0; i < natom; i++)
98
for(j=aostart[i]; j <= aostop[i]; j++) ao2atom[j] = i;
100
/* Get the overlap integrals -- these should be identical to AO S */
101
noei = nso*(nso+1)/2;
102
ss = init_array(noei);
103
stat = iwl_rdone(PSIF_OEI,PSIF_SO_S,ss,noei,0,0,outfile);
104
S = block_matrix(nso,nso);
105
for(i=0,ij=0; i < nso; i++)
106
for(j=0; j <= i; j++,ij++) {
107
S[i][j] = S[j][i] = ss[ij];
111
/* Compute nocc --- closed-shells only */
112
for(i=0,nocc=0; i < nirreps; i++) nocc += clsdpi[i];
114
fprintf(outfile, "\tNumber of doubly occupied orbitals: %d\n\n", nocc);
116
fprintf(outfile, "\tIter Pop. Localization Max. Rotation Angle Conv\n");
117
fprintf(outfile, "\t------------------------------------------------------------\n");
119
U = block_matrix(nocc, nocc);
120
V = block_matrix(nocc, nocc);
121
for(i=0; i < nocc; i++) V[i][i] = 1.0;
122
VV = block_matrix(nocc, nocc);
124
for(iter=0; iter < 100; iter++) {
127
for(i=nfzc; i < nocc; i++) {
128
for(A=0; A < natom; A++) {
131
for(l=aostart[A]; l <= aostop[A]; l++)
132
for(k=0; k < nso; k++)
133
PiiA += C[k][i] * C[l][i] * S[k][l];
139
/* Compute 2x2 rotations for Pipek-Mezey localization */
141
for(s=nfzc; s < nocc; s++) {
142
for(t=nfzc; t < s; t++) {
146
for(A=0; A < natom; A++) {
148
Pst = Pss = Ptt = 0.0;
150
for(l=aostart[A]; l <= aostop[A]; l++) {
151
for(k=0; k < nso; k++) {
152
Pst += 0.5 * (C[k][s] * C[l][t] +
153
C[l][s] * C[k][t]) * S[k][l];
155
Pss += C[k][s] * C[l][s] * S[k][l];
157
Ptt += C[k][t] * C[l][t] * S[k][l];
161
Ast += Pst * Pst - 0.25 * (Pss - Ptt) * (Pss - Ptt);
162
Bst += Pst * (Pss - Ptt);
166
/* Compute the rotation angle */
167
AB = Ast * Ast + Bst * Bst;
169
cos4a = -Ast/sqrt(AB);
170
alpha = 0.25 * acos(cos4a) * (Bst > 0 ? 1 : -1);
174
/* Keep up with the maximum 2x2 rotation angle */
175
alphamax = (fabs(alpha) > alphamax ? alpha : alphamax);
182
/* Now do the rotation */
183
for(k=0; k < nso; k++) {
186
C[k][s] = Uss * Cks + Ust * Ckt;
187
C[k][t] = Uts * Cks + Utt * Ckt;
190
zero_mat(U, nocc, nocc);
191
for(i=0; i < nocc; i++) U[i][i] = 1.0;
198
zero_mat(VV, nocc, nocc);
199
for(i=0; i < nocc; i++) {
200
for(j=0; j < nocc; j++) {
201
for(k=0; k < nocc; k++) {
202
VV[i][j] += V[i][k] * U[j][k];
207
for(i=0; i < nocc; i++)
208
for(j=0; j < nocc; j++)
214
conv = fabs(alphamax) - fabs(alphalast);
215
fprintf(outfile, "\t%4d %20.10f %20.10f %4.3e\n", iter, P, alphamax, conv);
216
if((iter > 2) && ((fabs(conv) < 1e-12) || alphamax == 0.0)) break;
217
alphalast = alphamax;
223
/* print_mat(V, nocc, nocc, outfile); */
225
/* Transform occupied orbital eigenvalues */
226
F = block_matrix(nocc, nocc);
227
for(i=0; i < nocc; i++)
228
for(j=0; j < nocc; j++)
229
for(k=0; k < nocc; k++)
230
F[i][j] += V[k][i] * evals[k] * V[k][j];
233
fprintf(outfile, "\nTransformed Orbital Energies:\n");
234
print_mat(F, nocc, nocc, outfile);
237
/* Compute a reordering array based on the diagonal elements of F */
238
orb_order = init_int_array(nocc);
239
orb_boolean = init_int_array(nocc);
240
for(i=0; i < nocc; i++) { orb_order[i] = 0; orb_boolean[i] = 0; }
242
for(i=0,max=0; i < nocc; i++) /* First, find the overall maximum */
243
if(fabs(F[i][i]) > fabs(F[max][max])) max = i;
245
orb_order[0] = max; orb_boolean[max] = 1;
247
for(i=1; i < nocc; i++) {
249
while(orb_boolean[max]) max++; /* Find an unused max */
250
for(j=0; j < nocc; j++)
251
if((fabs(F[j][j]) >= fabs(F[max][max])) && !orb_boolean[j]) max = j;
252
orb_order[i] = max; orb_boolean[max] = 1;
256
for(i=0; i < nocc; i++) fprintf(outfile, "%d %d\n", i, orb_order[i]);
260
fprintf(outfile, "\n\tPipek-Mezey Localized MO's (before sort):\n");
261
print_mat(C, nso, nmo, outfile);
264
/* Now reorder the localized MO's according to F */
265
Ctmp = block_matrix(nso,nocc);
266
for(i=0; i < nocc; i++)
267
for(j=0; j < nso; j++) Ctmp[j][i] = C[j][i];
269
for(i=0; i < nocc; i++) {
271
for(j=0; j < nso; j++) C[j][i] = Ctmp[j][iold];
272
evals[i] = F[iold][iold];
277
errcod = ip_boolean("PRINT_MOS", &(print), 0);
279
fprintf(outfile, "\n\tPipek-Mezey Localized MO's (after sort):\n");
280
print_mat(C, nso, nmo, outfile);
283
/* Check MO normalization */
285
for(i=0; i < nmo; i++) {
287
for(j=0; j < nso; j++)
288
for(k=0; k < nso; k++) {
289
norm += C[j][i] * C[k][i] * S[j][k];
292
fprintf(outfile, "norm[%d] = %20.10f\n", i, norm);
296
/* correct orbital phases for amplitude restarts */
297
chkpt_init(PSIO_OPEN_OLD);
298
scf_old = chkpt_rd_local_scf();
300
if (scf_old != NULL) {
301
MO_S = block_matrix(nmo, nmo);
302
X = block_matrix(nso, nmo);
303
C_DGEMM('n','n',nso, nmo, nso, 1, &(S[0][0]), nso, &(C[0][0]), nmo,
305
C_DGEMM('t','n',nmo, nmo, nso, 1, &(scf_old[0][0]), nmo, &(X[0][0]), nmo,
306
0, &(MO_S[0][0]), nmo);
310
fprintf(outfile, "Approximate Overlap Matrix\n");
311
print_mat(MO_S, nmo, nmo, outfile);
314
for(p=0; p < nmo; p++) {
316
for(q=0; q < nmo; q++) {
317
if(fabs(MO_S[p][q]) > max) {
318
max = fabs(MO_S[p][q]); max_col = q;
321
if(max_col != p) phase_ok = 0;
324
chkpt_init(PSIO_OPEN_OLD);
325
chkpt_wt_phase_check(phase_ok);
328
for(p=0; p < nmo; p++) {
329
if(MO_S[p][p] < 0.0) {
330
for(q=0; q < nso; q++)
341
/* Write the new MO's to chkpt */
342
chkpt_init(PSIO_OPEN_OLD);
344
chkpt_wt_local_scf(C);
350
fprintf(outfile, "\n\tLocalization of occupied orbitals complete.\n");
353
exit(PSI_RETURN_SUCCESS);
356
void init_io(int argc, char * argv[])
358
extern char *gprgid();
361
progid = (char *) malloc(strlen(gprgid())+2);
362
sprintf(progid, ":%s",gprgid());
364
psi_start(argc-1,argv+1,0);
365
ip_cwk_add(":INPUT"); /* for checking puream keyword */
375
fprintf(outfile, "\n");
376
fprintf(outfile, "\t\t\t**************************\n");
377
fprintf(outfile, "\t\t\t* *\n");
378
fprintf(outfile, "\t\t\t* LOCAL *\n");
379
fprintf(outfile, "\t\t\t* *\n");
380
fprintf(outfile, "\t\t\t**************************\n");
381
fprintf(outfile, "\n");
393
char *prgid = "LOCAL";
added
removed
src/bin/localize/local.c
