~ubuntu-branches/ubuntu/karmic/psicode/karmic

« back to all changes in this revision

Viewing changes to src/bin/cscf/guess.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
/*--------------------------------------------------------
2
 
 
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.
9
 
 
10
 
---------------------------------------------------------*/
11
 
#define EXTERN
12
 
#include "includes.h"
13
 
#include "common.h"
14
 
#include <libipv1/ip_lib.h>
15
 
 
16
 
void guess()
17
 
{
18
 
  struct symm *s;
19
 
  int i,j,bool,k,m,l,o;
20
 
  int errcod;
21
 
  int nc,no,nh,nn;
22
 
  int netmp=0;
23
 
  int size;
24
 
  int optri;
25
 
  char *guess_opt,*occ_fix_str;
26
 
  reftype reftmp;
27
 
  int mguess;
28
 
 
29
 
  /*
30
 
   ip_cwk_clear();
31
 
   ip_cwk_add(":DEFAULT");
32
 
   ip_cwk_add(":SCF");
33
 
  */
34
 
 
35
 
/*--- Set occupations to zero----*/
36
 
  for(i=0; i < num_ir; i++){
37
 
      scf_info[i].nclosed=0;
38
 
      scf_info[i].nopen=0;
39
 
    }
40
 
 
41
 
/* ---- Get reference type for comparison for restarting ----*/
42
 
  if(inflg == 1) reftmp = chkpt_rd_ref();
43
 
 
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() ) {
48
 
        occ_read();
49
 
  }
50
 
  else {
51
 
      
52
 
          errcod = ip_exist("DOCC",0);
53
 
          if(errcod) {
54
 
              fprintf(outfile,"\n  Using DOCC and SOCC to \n");
55
 
              fprintf(outfile,"  determine occupations\n\n");
56
 
              
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);
62
 
                }
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);
67
 
              }
68
 
              
69
 
              if(iopen || uhf) {
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);
75
 
                  }
76
 
                  
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);
82
 
                  }
83
 
              }
84
 
              
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 = 
89
 
                                ip_data("SOCC","%d"
90
 
                                        ,&scf_info[i].nopen,1,i);
91
 
                  if(iopen || uhf) errcod = 
92
 
                                ip_data("HOCC","%d"
93
 
                                        ,&scf_info[i].nhalf,1,i);
94
 
              }
95
 
             
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;
101
 
              if(netmp != nelec){
102
 
                 fprintf(outfile,"\n  DOCC and SOCC have the wrong number of");
103
 
                 fprintf(outfile,"\n  in them.  Should have %d total electrons"
104
 
                         ,nelec);
105
 
                 fprintf(outfile,"\n  and there are %d present\n\n\n",netmp);
106
 
                 exit(PSI_RETURN_FAILURE);
107
 
              }
108
 
          }
109
 
          else if(inflg == 1 && (reftmp == refnum) )
110
 
              occ_read();
111
 
              
112
 
          else{
113
 
              fprintf(outfile, "  Using core guess to determine occupations\n\n"); 
114
 
              
115
 
/* sort orbitals into two arrays, one with the symmetry label and 
116
 
   one with the eigenvalues */
117
 
              
118
 
              sortev();
119
 
              
120
 
/* calculate the occupations */
121
 
              
122
 
              occ_calc();
123
 
          }
124
 
  }
125
 
/* STB - 7/10/00 for DFT to ship the eigenvector */
126
 
/* calculate the number of total closed and open shells */
127
 
          if(uhf){
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
133
 
                      + scf_info[i].nopen;
134
 
              }
135
 
          }
136
 
          else{
137
 
              for(i=0; i < num_ir; i++){
138
 
                  n_closed += scf_info[i].nclosed;
139
 
              }
140
 
          }
141
 
/* output occupations to outfile */
142
 
          occ_out();
143
 
          
144
 
/* Setup occupation data for the calculation */
145
 
 
146
 
      n_open = 0;
147
 
 
148
 
/* Convert occupation data to symmetry info */
149
 
      
150
 
      for(i=0; i < num_ir; i++) {
151
 
          s=&scf_info[i];
152
 
          
153
 
          if (s->nopen|| s->nhalf) n_open++;
154
 
          if (nn = s->num_so) {
155
 
              nc = s->nclosed;
156
 
              no = s->nopen;
157
 
              nh = s->nhalf;
158
 
              
159
 
              if (  (nn < nc + no + nh)
160
 
                    ||(nc < 0)
161
 
                    ||(no < 0)
162
 
                    ||(nh < 0)) {
163
 
                  const char* fmt
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);
168
 
                } 
169
 
              
170
 
              if(uhf){
171
 
                  spin_info[0].scf_spin[i].noccup = nc+no;
172
 
                  spin_info[1].scf_spin[i].noccup = nc;
173
 
                  if(nh != 0){
174
 
                      fprintf(outfile,"\nCannot use HOCC with UHF\n");
175
 
                      exit(PSI_RETURN_FAILURE);
176
 
                  }
177
 
              }
178
 
              
179
 
              if(iopen) {
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]);
185
 
              }
186
 
              if(twocon) {
187
 
                  s->pmat2 = (double *) init_array(ioff[nn]);
188
 
                  s->pmato2 = (double *) init_array(ioff[nn]);
189
 
              }
190
 
              if(uhf){
191
 
                  for(j=0; j < nc+no; j++)
192
 
                      spin_info[0].scf_spin[i].occ_num[j] = 1.0;
193
 
                  
194
 
                  for(j=0;j < nc; j++)
195
 
                      spin_info[1].scf_spin[i].occ_num[j] = 1.0;
196
 
 
197
 
              }
198
 
              else{
199
 
                  for (j=0; j < nc ; j++) {
200
 
                      s->occ_num[j] = 2.0;
201
 
                  }
202
 
                  for (j=nc; j < nc+no ; j++) {
203
 
                      s->occ_num[j] = 1.0;
204
 
                  }
205
 
                  if(nh) s->occ_num[nc+no] = nh*0.5;
206
 
              }
207
 
          }
208
 
      }
209
 
      
210
 
      optri = n_open*(n_open+1)/2;
211
 
      
212
 
      if (iopen) {
213
 
          int mm1=1,mm2=1;
214
 
          if(optri == 0){
215
 
              fprintf(outfile,"\nNot an open shell molecule.\n");
216
 
              fprintf(outfile,"Re-check opentype!!!!\n\n");
217
 
              exit(PSI_RETURN_FAILURE);
218
 
          }
219
 
          
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;
224
 
          
225
 
          alpha = (double *) init_array(ioff[optri]);
226
 
          beta = (double *) init_array(ioff[optri]);
227
 
          
228
 
          if (twocon) {
229
 
              if(n_open == 2) {
230
 
                  alpha[0] = 0.0;
231
 
                  alpha[1] = 0.0;
232
 
                  alpha[2] = 0.0;
233
 
                  beta[0] = 0.0;
234
 
                  beta[1] = 1.0;
235
 
                  beta[2] = 0.0;
236
 
              }
237
 
              else {
238
 
                  fprintf(outfile,
239
 
                          "this program cannot handle same symmetry\n");
240
 
                  fprintf(outfile," tcscf. try SCFX\n");
241
 
                  exit(PSI_RETURN_FAILURE);
242
 
              }
243
 
          }
244
 
          else if(singlet) {
245
 
              if(n_open == 2) {
246
 
                  alpha[0] = 0.0;
247
 
                  alpha[0] = 0.0;
248
 
                  alpha[1] = 0.0;
249
 
                  alpha[2] = 0.0;
250
 
                  beta[0] = 1.0;
251
 
                  beta[1] = -3.0;
252
 
                  beta[2] = 1.0;
253
 
              }
254
 
              else {
255
 
                  fprintf(outfile,
256
 
                          "this program cannot handle same symmetry\n");
257
 
                  fprintf(outfile," singlets. try SCFX\n");
258
 
                  exit(PSI_RETURN_FAILURE);
259
 
              }
260
 
          }
261
 
          else if(hsos) {
262
 
              for(i=0; i < optri ; i++) {
263
 
                  alpha[i]=0.0;
264
 
                  beta[i]=1.0;
265
 
              }
266
 
          }
267
 
          else {
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);
271
 
                  beta[i] = -beta[i];
272
 
              }
273
 
          }
274
 
          for (i=0; i < optri; i++) {
275
 
              if(iter == 0 || print & 1)
276
 
                  fprintf(outfile,"        %d  %d       %f     %f\n",mm1,mm2,
277
 
                          alpha[i],-beta[i]);
278
 
              mm2++;
279
 
              if (mm2 > mm1) {
280
 
                  mm1++;
281
 
                  mm2 = 1;
282
 
              }
283
 
          }
284
 
      }
285
 
}       
286
 
 
287