2
* Revision 1.2 2004/05/03 04:32:40 crawdad
3
* Major mods based on merge with stable psi-3-2-1 release. Note that this
4
* version has not been fully tested and some scf-optn test cases do not run
5
* correctly beccause of changes in mid-March 2004 to optking.
8
/* Revision 1.1.1.1.12.2 2004/04/21 15:45:07 evaleev
9
/* Modified DIIS algorithm for RHF and ROHF to work in OSO basis rather than in
10
/* AO basis, to avoid difficulties of transforming between MO and AO bases
11
/* when linear dependencies are present.
13
/* Revision 1.1.1.1.12.1 2004/04/06 21:29:05 crawdad
14
/* Corrections to the RHF/ROHF DIIS algorithm, which was simply incorrect.
15
/* The backtransformation of the DIIS error vectors to the AO basis was not
16
/* mathematically right.
19
/* Revision 1.1.1.1 2000/02/04 22:52:29 evaleev
20
/* Started PSI 3 repository
22
/* Revision 1.3 1999/11/02 23:55:55 localpsi
23
/* Shawn Brown - (11/2/99) Modified to the code in a few major ways.
25
/* 1. Added the capability to do UHF. All of the features available with the
26
/* other refrences have been added for UHF.
28
/* 2. For UHF, I had to alter the structure of file30. (See cleanup.c for a
29
/* map) This entailed adding a pointer array right after the header in the SCF
30
/* section of file30 that pointed to all of the data for the SCF caclulation.
31
/* Functions were added to libfile30 to account for this and they are
32
/* incorporated in this code.
34
/* 3. Updated and fixed all of the problems associated with my previous
35
/* guessing code. The code no longer uses OPENTYPE to specify the type of
36
/* occupation. The keword REFERENCE and MULTP can now be used to indicate any
37
/* type of calculation. (e.g. ROHF with MULTP of 1 is an open shell singlet
38
/* ROHF calculation) This code was moved to occ_fun.c. The code can also
39
/* guess at any multplicity in a highspin case, provided enough electrons.
41
/* Revision 1.2 1999/08/17 19:04:13 evaleev
42
/* Changed the default symmetric orthogonalization to the canonical
43
/* orthogonalization. Now, if near-linear dependencies in the basis are found,
44
/* eigenvectors of the overlap matrix with eigenvalues less than 1E-6 will be
45
/* left out. This will lead to num_mo != num_so, i.e. SCF eigenvector is no
46
/* longer a square matrix. Had to rework some routines in libfile30, and add some.
47
/* The progrem prints out a warning if near-linear dependencies are found. TRANSQT
48
/* and a whole bunch of other codes has to be fixed to work with such basis sets.
50
/* Revision 1.1.1.1 1999/04/12 16:59:25 evaleev
51
/* Added a version of CSCF that can work with CINTS.
55
static char *rcsid = "$Id: diis.c 2455 2004-05-03 04:32:41Z crawdad $";
63
static double *btemp, **bold, **bmat;
64
static struct diis_mats {
65
double ***fock_c; /* Closed-shell Fock matrix in AO basis */
66
double ***fock_o; /* Open-shell Fock matrix in AO basis */
67
double ***error; /* Error matrix in AO basis */
70
void diis(scr1,scr2,scr3,c1,c2,cim,newci)
72
double **scr1, **scr2, **scr3,*c1,*c2;
82
double etemp, dotp, norm, determ, etempo;
87
double **C; /* AO->MO transform = P^-1 * C */
88
double **X, **S; /* scratch matrix */
92
bmat = (double **) block_matrix(ndiis+1,ndiis+1);
93
bold = (double **) block_matrix(ndiis,ndiis);
94
btemp = (double *) init_array(ndiis+1);
96
diism = (struct diis_mats *) malloc(sizeof(struct diis_mats)*ndiis);
98
for(i=0; i < ndiis ; i++) {
99
diism[i].fock_c = (double ***) malloc(sizeof(double **)*num_ir);
101
diism[i].fock_o = (double ***) malloc(sizeof(double **)*num_ir);
102
diism[i].error = (double ***) malloc(sizeof(double **)*num_ir);
103
for(j=0; j < num_ir ; j++) {
104
if(nn=scf_info[j].num_so) {
105
num_mo = scf_info[j].num_mo;
106
diism[i].fock_c[j] = (double **) block_matrix(nn,nn);
107
if(iopen) diism[i].fock_o[j] = (double **) block_matrix(nn,nn);
108
diism[i].error[j] = (double **) block_matrix(num_mo,num_mo);
114
scale = 1.0 + dampsv;
120
for (i=0; i < last ; i++) {
121
diism[i] = diism[i+1];
126
/* Debugging stuff */
128
errcod = ip_boolean("DIIS_PRINT",&diis_print,0);
129
if(iter == 1 && diis_print)
130
ffile(&diis_out,"diis_out.dat",0);
132
/* save ao fock matrices in fock_save */
135
for (m=0; m < num_ir ; m++) {
139
tri_to_sq(s->fock_pac,d->fock_c[m],nn);
140
if(iopen) tri_to_sq(s->fock_open,d->fock_o[m],nn);
142
/* form error matrix in mo basis */
143
mmult(s->cmat,1,d->fock_c[m],0,scr1,0,num_mo,nn,nn,0);
144
mmult(scr1,0,s->cmat,0,scr2,0,num_mo,nn,num_mo,0);
146
mmult(s->cmat,1,d->fock_o[m],0,scr1,0,num_mo,nn,nn,0);
147
mmult(scr1,0,s->cmat,0,scr3,0,num_mo,nn,num_mo,0);
150
for (i=0; i < num_mo; i++) {
151
occi = s->occ_num[i];
152
for (j=0; j <= i ; j++ ) {
153
occj = s->occ_num[j];
155
if (occi == 0.0 && occj != 0.0 ) {
156
scr1[i][j]= scr2[i][j];
157
scr1[j][i]= scr2[i][j];
160
scr1[i][j]=scr1[j][i]=0.0;
164
if ((occi == 1.0 || occi == 0.5) && occj == 2.0 )
165
etemp = scr2[i][j]-0.5*scr3[i][j];
166
else if(occi == 0.0) {
167
if (occj == 2.0) etemp = scr2[i][j];
168
else if (occj != 0.0) etemp = 0.5*scr3[i][j];
172
scr1[i][j] = scr1[j][i] = etemp;
175
if (occi != 2.0 && occi != 0.0 && occj == 2.0 )
176
etemp = cim*(c1[m]*scr2[i][j]-c2[m]*scr3[i][j]);
177
else if(occi == 0.0) {
178
if (occj == 2.0) etemp = scr2[i][j];
179
else if (occj != 0.0) etemp = cim*scr3[i][j];
183
scr1[i][j] = scr1[j][i] = etemp;
188
/* transform the error matrix into the OSO basis */
189
X = block_matrix(num_mo, num_mo);
190
C_DGEMM('n', 'n', num_mo, num_mo, num_mo, 1.0, s->ucmat[0], num_mo, scr1[0], nsfmax, 0.0, X[0], num_mo);
191
C_DGEMM('n', 't', num_mo, num_mo, num_mo, 1.0, X[0], num_mo, s->ucmat[0], num_mo, 0.0, d->error[m][0], num_mo);
194
for(i=0; i < num_mo ; i++) {
195
for(j=0; j <= i ; j++) {
196
etemp=fabs(scr1[i][j]);
197
diiser = MAX0(diiser,etemp);
203
/* then set up b matrix */
206
for (i=0; i < last ; i++) {
207
for (j=0; j <= i ; j++) {
208
bold[i][j]=bold[j][i]=bold[i+1][j+1];
212
for (i=0; i <= last ; i++) {
214
for (m=0; m < num_ir ; m++) {
218
sdot(diism[i].error[m],diism[last].error[m],num_mo,&dotp);
222
bold[i][last]=bold[last][i] = etemp;
227
norm = 1.0/bold[0][0];
228
for (i=1; i <= last+1 ; i++) {
229
bmat[i][0]=bmat[0][i] = -1.0;
231
for (j=1; j <= i ; j++) {
232
bmat[i][j]=bmat[j][i] = bold[i-1][j-1]*norm;
233
if(i==j) bmat[i][j] *= scale;
238
fprintf(diis_out,"\nBMAT for iter %d",iter);
239
print_mat(bmat,col,col,diis_out);
243
flin(bmat,btemp,col,1,&determ);
245
/* test for poorly conditioned equations */
246
while (fabs(determ) < 1.0e-19 && try < last) {
253
norm=1.0/bold[try][try];
254
for (i=1; i <= ndiis-try ; i++) {
255
bmat[i][0]=bmat[0][i] = -1.0;
256
for (j=1; j <= i ; j++) {
257
bmat[i][j]=bmat[j][i]=bold[i+try-1][j+try-1]*norm;
258
if(i==j) bmat[i][j] *= scale;
264
fprintf(diis_out,"\nCorrected BMAT for iter %d",iter);
265
print_mat(bmat,col,col,diis_out);
268
flin(bmat,btemp,col,1,&determ);
271
if(fabs(determ) < 10.0e-20) {
272
printf(" try %d no good\n",try);
276
if(iter >= it_diis) {
277
for (m=0; m < num_ir ; m++) {
280
for (i=ij=0; i < nn ; i++) {
281
for (j=0; j <= i ; j++,ij++) {
285
for (k=try; k < last+1 ; k++) {
286
if(iopen) etempo += btemp[kk]*diism[k].fock_o[m][i][j];
287
etemp += btemp[kk]*diism[k].fock_c[m][i][j];
290
if(iopen) s->fock_open[ij] = etempo;
291
s->fock_pac[ij] = etemp;
added
removed
src/bin/cscf/diis.c
