1
/*--------------------------------------------------------
3
guess.c: Function that reads the guessing
4
parameters from input and then either
5
uses them to form an initial guess or
6
calculates an initial guess from
7
multiplicity and the charge using a
8
diagonalization of the core guess.
10
---------------------------------------------------------*/
14
#include <libipv1/ip_lib.h>
25
char *guess_opt,*occ_fix_str;
31
ip_cwk_add(":DEFAULT");
35
/*--- Set occupations to zero----*/
36
for(i=0; i < num_ir; i++){
37
scf_info[i].nclosed=0;
41
/* ---- Get reference type for comparison for restarting ----*/
42
if(inflg == 1) reftmp = chkpt_rd_ref();
44
/*---If first calculation either read occupations from DOCC or form an initial guess---*/
45
/* check checkpoint flag to see if occupations in ckhpt file should be forced;
46
this is needed to override DOCC vectors in finite difference calculations */
47
if ( chkpt_rd_override_occ() ) {
52
errcod = ip_exist("DOCC",0);
54
fprintf(outfile,"\n Using DOCC and SOCC to \n");
55
fprintf(outfile," determine occupations\n\n");
57
errcod = ip_count("DOCC",&size,0);
58
if(errcod == IPE_OK && size != num_ir) {
59
fprintf(outfile,"\n DOCC array is the wrong size\n");
60
fprintf(outfile," is %d, should be %d\n",size,num_ir);
61
exit(PSI_RETURN_FAILURE);
63
if(errcod != IPE_OK && !iopen) {
64
fprintf(outfile,"\n try adding some electrons buddy!\n");
65
fprintf(outfile," need DOCC\n");
66
exit(PSI_RETURN_FAILURE);
70
errcod = ip_count("SOCC",&size,0);
71
if(errcod == IPE_OK && size != num_ir) {
72
fprintf(outfile,"\n SOCC array is the wrong size\n");
73
fprintf(outfile," is %d, should be %d\n",size,num_ir);
74
exit(PSI_RETURN_FAILURE);
77
errcod = ip_count("HOCC",&size,0);
78
if(errcod == IPE_OK && size != num_ir) {
79
fprintf(outfile,"\n HOCC array is the wrong size\n");
80
fprintf(outfile," is %d, should be %d\n",size,num_ir);
81
exit(PSI_RETURN_FAILURE);
85
for (i=0; i < num_ir; i++) {
86
errcod = ip_data("DOCC","%d"
87
,&scf_info[i].nclosed,1,i);
88
if(iopen || uhf) errcod =
90
,&scf_info[i].nopen,1,i);
91
if(iopen || uhf) errcod =
93
,&scf_info[i].nhalf,1,i);
96
/* STB (10/29/99) Make sure that DOCC and SOCC have the right
97
number of electrons in them */
98
for(i = 0;i<num_ir;i++)
99
netmp += (2*scf_info[i].nclosed)
100
+ scf_info[i].nopen + scf_info[i].nhalf;
102
fprintf(outfile,"\n DOCC and SOCC have the wrong number of");
103
fprintf(outfile,"\n in them. Should have %d total electrons"
105
fprintf(outfile,"\n and there are %d present\n\n\n",netmp);
106
exit(PSI_RETURN_FAILURE);
109
else if(inflg == 1 && (reftmp == refnum) )
113
fprintf(outfile, " Using core guess to determine occupations\n\n");
115
/* sort orbitals into two arrays, one with the symmetry label and
116
one with the eigenvalues */
120
/* calculate the occupations */
125
/* STB - 7/10/00 for DFT to ship the eigenvector */
126
/* calculate the number of total closed and open shells */
128
for(i=0;i<num_ir;i++){
129
b_elec += scf_info[i].nclosed;
130
spin_info[1].scf_spin[i].nclosed = scf_info[i].nclosed;
131
a_elec += scf_info[i].nclosed+scf_info[i].nopen;
132
spin_info[0].scf_spin[i].nclosed = scf_info[i].nclosed
137
for(i=0; i < num_ir; i++){
138
n_closed += scf_info[i].nclosed;
141
/* output occupations to outfile */
144
/* Setup occupation data for the calculation */
148
/* Convert occupation data to symmetry info */
150
for(i=0; i < num_ir; i++) {
153
if (s->nopen|| s->nhalf) n_open++;
154
if (nn = s->num_so) {
159
if ( (nn < nc + no + nh)
164
= "cscf: invalid number of electrons in irrep %d\n";
165
fprintf(stderr,fmt,i);
166
fprintf(outfile,fmt,i);
167
exit(PSI_RETURN_FAILURE);
171
spin_info[0].scf_spin[i].noccup = nc+no;
172
spin_info[1].scf_spin[i].noccup = nc;
174
fprintf(outfile,"\nCannot use HOCC with UHF\n");
175
exit(PSI_RETURN_FAILURE);
180
s->fock_eff = (double *) init_array(ioff[nn]);
181
s->fock_open = (double *) init_array(ioff[nn]);
182
s->pmato = (double *) init_array(ioff[nn]);
183
s->dpmato = (double *) init_array(ioff[nn]);
184
s->gmato = (double *) init_array(ioff[nn]);
187
s->pmat2 = (double *) init_array(ioff[nn]);
188
s->pmato2 = (double *) init_array(ioff[nn]);
191
for(j=0; j < nc+no; j++)
192
spin_info[0].scf_spin[i].occ_num[j] = 1.0;
195
spin_info[1].scf_spin[i].occ_num[j] = 1.0;
199
for (j=0; j < nc ; j++) {
202
for (j=nc; j < nc+no ; j++) {
205
if(nh) s->occ_num[nc+no] = nh*0.5;
210
optri = n_open*(n_open+1)/2;
215
fprintf(outfile,"\nNot an open shell molecule.\n");
216
fprintf(outfile,"Re-check opentype!!!!\n\n");
217
exit(PSI_RETURN_FAILURE);
220
if(iter == 0 || print & 1){
221
fprintf(outfile,"\n open-shell energy coeffs\n");
222
fprintf(outfile," open shell pair alpha beta\n");}
223
optri = (n_open*(n_open+1))/2;
225
alpha = (double *) init_array(ioff[optri]);
226
beta = (double *) init_array(ioff[optri]);
239
"this program cannot handle same symmetry\n");
240
fprintf(outfile," tcscf. try SCFX\n");
241
exit(PSI_RETURN_FAILURE);
256
"this program cannot handle same symmetry\n");
257
fprintf(outfile," singlets. try SCFX\n");
258
exit(PSI_RETURN_FAILURE);
262
for(i=0; i < optri ; i++) {
268
for (i=0; i < optri ; i++) {
269
errcod = ip_data("ALPHA","%lf",&alpha[i],1,i);
270
errcod = ip_data("BETA","%lf",&beta[i],1,i);
274
for (i=0; i < optri; i++) {
275
if(iter == 0 || print & 1)
276
fprintf(outfile," %d %d %f %f\n",mm1,mm2,