2
/*************************************************************************/
5
/* Written by Edward Seidl (a NON-hog) */
7
/* Parts liberally ripped off from HFS group SCF */
8
/* code written in FORTRAN (Boo!!!) */
10
/* modified April 17, 1991 to use new input format developed by */
12
/**************************************************************************/
14
/* Description of input */
17
/* This is a character string to be included in the output. */
18
/* There is no default. */
21
/* This is the type of wavefunction which is ultimately desired. */
22
/* The default is SCF. */
24
/* DIRECT_SCF = boolean */
25
/* Flag to request the direct formation of the Fock matrix */
27
/* OPENTYPE = string */
28
/* This specifies the state desired. It can be one of NONE */
29
/* (for a closed shell singlet), SINGLET (for an open shell */
30
/* singlet), HIGHSPIN (for any high spin open shell system), */
31
/* TWOCON (for a two configuration singlet), or SPECIAL. */
32
/* If SPECIAL is given, then alpha and beta coupling */
33
/* coefficients must be given with the ALPHA and BETA keywords. */
34
/* The default is NONE. */
36
/* DOCC = integer_vector */
37
/* This gives the number of doubly occupied orbitals in each */
38
/* irreducible representation. There is no default. */
40
/* SOCC = integer_vector */
41
/* This gives the number of singly occupied orbitals in each */
42
/* irreducible representation. If OPENTYPE = NONE this defaults */
43
/* to the zero vector. Otherwise, there is no default. */
45
/* DERTYPE = string */
46
/* This specifies the order of derivative that is to even- */
47
/* tually be done. It is used by the scf program to */
48
/* determine if certain files are to be written and it is */
49
/* also used to determine the default convergence of the */
50
/* wavefunction. The default is FIRST. */
52
/* MAXITER = integer */
53
/* This gives the maximum number of iterations. The default */
56
/* CONVERGENCE = integer */
57
/* The convergence criterion is 10**(-integer). The default is */
58
/* 7 if both DERTYPE = NONE and WFN = SCF are given and 10 */
61
/* LEVELSHIFT = real */
62
/* This specifies the level shift. The default is 1.0. */
65
/* There are also a large number of less commonly used input */
66
/* parameters. If you do not understand what the following */
67
/* options mean, then make sure that they do not appear in your */
68
/* input. The defaults will work in the overwhelming majority */
69
/* of cases. These are specified with the following keywords: */
72
/* REORDER = boolean */
73
/* The molecular orbitals will be reordered if this is */
74
/* true, in which case, the MOORDER parameter must be */
75
/* present. The default is false. */
77
/* MOORDER = integer_vector */
78
/* This specifies a molecular orbital reordering vector. */
79
/* It will only be used if REORDER = YES. This vector */
80
/* contains first the ordering for the orbitals in the */
81
/* first irreducible representation and then the second */
82
/* and so on. The first orbital of each irreducible */
83
/* representation is numbered 1. There is no default. */
85
/* ALPHA = real_vector */
86
/* If OPENTYPE = SPECIAL, then this parameter gives the */
87
/* alpha coupling coefficients. The number of elements in */
88
/* this vector is MM(MM+1)/2, where MM is the number of */
89
/* irreducible representations containing singly occupied */
90
/* molecular orbitals. There is no default. */
92
/* BETA = real_vector */
93
/* If OPENTYPE = SPECIAL, then this parameter gives the */
94
/* beta coupling coefficients. The number of elements in */
95
/* this vector is MM(MM+1)/2, where MM is the number of */
96
/* irreducible representations containing singly occupied */
97
/* molecular orbitals. There is no default. */
99
/* RESTART = boolean */
100
/* The calculation will restart from the old wavefunction */
101
/* if RESTART is true. If the old wavefunction does not */
102
/* exist, then the cscf program will generate its own ini- */
103
/* tial guess automatically. Possible values for this */
104
/* parameter are TRUE, YES, 1, FALSE, NO, and 0. The */
105
/* default is true. */
107
/* IPRINT = integer */
108
/* This is a print option. The default is 0. */
110
/* ROTATE = boolean */
111
/* The molecular orbitals will not be rotated if this is */
112
/* false. The rotation only affects the virtual orbitals */
113
/* for open shell systems. This parameter must be true */
114
/* for correlated gradients and it must be false for */
115
/* second and higher derivatives. The default is false if */
116
/* WFN = SCF and true otherwise. */
119
/* This determines whether diis will be used. The default is */
120
/* false for OPENTYPE = TWOCON and true otherwise. */
122
/* NDIIS = integer */
123
/* This gives the number of error matrices to use in the diis */
124
/* procedure. The default is 6 for closed shell, 4 for open */
125
/* shell, and 3 for tcscf. */
127
/* DIISSTART = integer */
128
/* This gives the first iteration for which DIIS will be */
129
/* used. The default is 0. */
131
/* DIISDAMP = real */
132
/* This gives the damping factor for the diis procedure. The */
133
/* default is 0.0 for closed shell, 0.02 for open shell, and */
134
/* 0.01 for tcscf. */
137
/* This is used in tcscf to determine how often the ci */
138
/* coefficients are recalculated. A small number (~0.25) */
139
/* will cause them to be recalculated nearly every scf */
140
/* iteration. The default is 0.5. */
143
/* FOCK_TYPE = integer */
144
/* Only used for tcscf and open shell calculations. */
145
/* If FOCK_TYPE = 0 use a simple form for fock_eff */
146
/* " " = 1 use form of fock_eff suitable for */
147
/* high-spin cases */
148
/* " " > 1 experimental fock matrices */
150
/**************************************************************************/
154
/*-------------------------------------------------------------------------
155
READ THIS FIRST: This code is a hack of the original CSCF which employed
156
symmetric orthogonalization procedure. This code uses the canonical
157
orthogonalization procedure. To decrease the amount of
158
rewriting I had to do I left the initialization part (init_scf and
159
init_scf2) intact. Thus, although I have to allocate more space than
160
I need, I should be fine otherwise.
162
- Edward Valeev, August'99
163
-------------------------------------------------------------------------*/
166
static char *rcsid = "$Id: cscf.c,v 1.21.4.1 2003/12/31 01:59:54 crawdad Exp $";
168
#include "includes.h"
170
#include <libipv1/ip_lib.h>
171
#include <libpsio/psio.h>
172
#include <libchkpt/chkpt.h>
173
#include <libqt/qt.h>
175
void print_initial_vec();
176
extern void write_scf_matrices(void);
183
char *prog_name="CSCF3.0: An SCF program written in C";
184
char *output="APPEND ";
186
ip_value_t *ipvalue=NULL;
187
int errcod, orthog_only, mo_print;
190
errcod = psi_start(argc-1, argv+1, 0);
191
if (errcod != PSI_RETURN_SUCCESS)
192
exit(PSI_RETURN_FAILURE);
196
fprintf(outfile,"\n%13c------------------------------------------\n",' ');
197
fprintf(outfile,"\n%16c%s\n",' ',prog_name);
198
fprintf(outfile,"\n%14cWritten by too many people to mention here\n",' ');
199
fprintf(outfile,"\n%13c------------------------------------------\n",' ');
203
itap33 = PSIF_SO_TEI;
213
itapDSCF = PSIF_DSCF;
214
itap92 = PSIF_SO_PKSUPER1;
215
itap93 = PSIF_SO_PKSUPER2;
217
/* JPK 6/1/00 integral accuracy: dynamic(default)=1, static=0 */
219
eri_cutoff = 1.0E-14;
220
ip_boolean("DYN_ACC",&dyn_acc,0);
224
/* CDS 3/6/02 add flag to do only orthogonalization */
225
errcod = ip_string("WFN",&wfn,0);
226
if (strcmp(wfn,"DETCAS")==0) orthog_only = 1;
227
else orthog_only = 0;
228
ip_boolean("ORTHOG_ONLY",&orthog_only,0);
231
/* open integrals file(s) */
235
/* STB (6/30/99) - Function added because in order to initialize things
236
one must know whether you are doing UHF or restricted */
238
chkpt_init(PSIO_OPEN_OLD);
242
/* initialize some constants and arrays */
248
/* read input.dat, get occupations, flags, etc. */
252
/* we can't just orthogonalize the orbitals if there aren't any */
253
if (inflg != 1) orthog_only = 0;
255
/* set up other useful arrays */
259
/* get one electron integrals */
264
for (i=0; i < num_ir; i++) {
267
fprintf(outfile,"\nsmat for irrep %s\n",s->irrep_label);
268
print_array(s->smat,nn,outfile);
269
fprintf(outfile,"\ntmat for irrep %s\n",s->irrep_label);
270
print_array(s->tmat,nn,outfile);
271
fprintf(outfile,"\nhmat for irrep %s\n",s->irrep_label);
272
print_array(s->hmat,nn,outfile);
277
/* form S-1/2 matrix sahalf */
282
for (i=0; i < num_ir ; i++) {
285
fprintf(outfile,"\nsahalf for irrep %s\n",s->irrep_label);
286
print_mat(s->sahalf,nn,s->num_mo,outfile);
291
/* if no initial vector, form one from core hamiltonian */
293
if (inflg == 2) form_vec();
296
/* guess or designate orbital occupations*/
299
/* orthogonalize old vector and form first density matrix */
307
/* Print out the first vector */
311
/* if we are only orthogonalizing, then quit here */
313
fprintf(outfile, "Only orbital orthogonalization has been performed\n");
315
errcod = ip_boolean("PRINT_MOS",&mo_print,0);
316
if (mo_print) print_mos();
317
write_scf_matrices();
321
exit(PSI_RETURN_SUCCESS);
328
/*--- Prepare alpha and beta eigenvectors from
329
the core Hamiltonian guess ---*/
337
/* Decide how to form the Fock matrix */
341
/* read in tei's and form supermatrix */
346
Pmat.key = strdup("P-supermatrix");
347
Pmat.bufpos = PSIO_ZERO;
349
PKmat.key = strdup("PK-supermatrix");
350
PKmat.bufpos = PSIO_ZERO;
351
psio_open(Pmat.unit,PSIO_OPEN_NEW);
352
psio_open(PKmat.unit,PSIO_OPEN_NEW);
357
/* form the Fock matrix directly */
358
/*check for rohf singlet...doesn't work direct*/
359
if(!uhf && singlet) {
360
fprintf(outfile,"\n rohf open shell singlet doesn't work direct\n");
361
fprintf(outfile," remove 'direct_scf = true' from input\n");
362
fprintf(stderr,"rohf open shell singlet doesn't work direct\n");
363
fprintf(stderr,"remove 'direct_scf = true' from input\n");
366
exit(PSI_RETURN_FAILURE);
370
if(dyn_acc) fprintf(outfile,"\n Using inexpensive integrals");
378
if(twocon) scf_iter_2();
379
else if(uhf) uhf_iter();
388
void print_initial_vec()
393
for(irrep=0;irrep<num_ir;irrep++) {
394
s = &scf_info[irrep];
398
fprintf(outfile,"\n Initial alpha vector for irrep %s\n",s->irrep_label);
399
print_mat(spin_info[0].scf_spin[irrep].cmat,nso,nmo,outfile);
400
fprintf(outfile,"\n Initial beta vector for irrep %s\n",s->irrep_label);
401
print_mat(spin_info[1].scf_spin[irrep].cmat,nso,nmo,outfile);
404
fprintf(outfile,"\nInitial vector for irrep %s\n",s->irrep_label);
405
print_mat(s->cmat,nso,nmo,outfile);
added
removed
src/bin/cscf/cscf.c
