← Back to branch summary

~ubuntu-branches/ubuntu/precise/psicode/precise

« back to all changes in this revision

Viewing changes to src/bin/dboc/dboc.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-06-07 16:49:57 UTC
  • mfrom: (2.1.2 hardy)
  • Revision ID: james.westby@ubuntu.com-20080607164957-8pifvb133yjlkagn
Tags: 3.3.0-3
http://bugs.debian.org/464867
* debian/rules (DEB_MAKE_CHECK_TARGET): Do not abort test suite on
  failures.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Set ${bindir} to /usr/lib/psi.
* debian/rules (install/psi3): Move psi3 file to /usr/bin.
* debian/patches/07_464867_move_executables.dpatch: New patch, add
  /usr/lib/psi to the $PATH, so that the moved executables are found.
  (closes: #464867)
* debian/patches/00list: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/dboc/dboc.cc
20
20
#include <libbasis/overlap.h>
21
21
#include <masses.h>
22
22
#include <physconst.h>
 
23
#include "defines.h"
23
24
#include "molecule.h"
24
25
#include "moinfo.h"
25
26
#include "params.h"
 
27
#include "hfwfn.h"
26
28
 
27
29
#if defined HAVE_DECL_SETENV && !HAVE_DECL_SETENV
28
30
  extern int setenv(const char *, const char *, int);
37
39
static double eval_dboc();
38
40
extern void print_intro();
39
41
extern void print_params();
 
42
extern void print_geom();
40
43
extern "C" char *gprgid();
41
44
extern void setup_geoms();
42
 
extern double eval_derwfn_overlap();
 
45
extern double eval_derwfn_overlap(bool symm_coord);
43
46
extern void read_moinfo();
44
47
 
45
48
/*--- Global structures ---*/
48
51
Molecule_t Molecule;
49
52
MOInfo_t MOInfo;
50
53
Params_t Params;
51
 
BasisSet* BasisDispP;
52
 
BasisSet* BasisDispM;
53
 
 
54
 
#define MAX_GEOM_STRING 20
 
54
 
 
55
BasisSet* BasisSets[MAX_NUM_DISP];  // Array of pointers to basis set objects with current coordinate displaced by +delta, -delta, +2delta, and -2delta, respectively
 
56
HFWavefunction* HFVectors[MAX_NUM_DISP];  // Array of pointers to HF wavefunctions for displacements by +delta, -delta, +2delta, and -2delta, respectively
 
57
char* CI_Vector_Labels[MAX_NUM_DISP] = {
 
58
  "R0-delta CI Vector",
 
59
  "R0+delta CI Vector",
 
60
  "R0-2delta CI Vector",
 
61
  "R0+2delta CI Vector"
 
62
  };
 
63
 
 
64
const int MAX_GEOM_STRING=20;
55
65
 
56
66
int main(int argc, char *argv[])
57
67
{
62
72
  FILE *geometry;
63
73
 
64
74
  init_io(argc,argv);
65
 
  print_intro();
66
75
  read_molecule();
67
76
  parsing();
 
77
  print_intro();
 
78
  print_geom();
68
79
  print_params();
 
80
#if USE_MOINFO
69
81
  read_moinfo();
 
82
#endif
70
83
  setup_geoms();
71
84
  double E_dboc = eval_dboc();
72
85
  fprintf(outfile,"  E(DBOC) = %25.15lf a.u.\n",E_dboc);
73
 
  fprintf(outfile,"  E(DBOC) = %25.10lf cm^{-1}\n\n",E_dboc*_hartree2wavenumbers);
 
86
  fprintf(outfile,"  E(DBOC) = %25.5lf cm^{-1}\n\n",E_dboc*_hartree2wavenumbers);
74
87
  exit_io();
75
88
  exit(0);
76
89
}
108
121
  Params.delta = 0.0005;
109
122
  errcod = ip_data(":DBOC:DISPLACEMENT","%lf",&Params.delta,0);
110
123
 
 
124
  // how much memory to use?
 
125
  long max_memory;
 
126
  fndcor(&max_memory,infile,outfile);
 
127
  Params.max_memory = static_cast<size_t>(max_memory);
 
128
  Params.memory = Params.max_memory;
 
129
 
 
130
  // number of compute threads
 
131
  int num_threads = 1;
 
132
  errcod = ip_data("NUM_THREADS","%d",&num_threads,0);
 
133
  Params.num_threads = num_threads;
 
134
 
111
135
  Params.print_lvl = 1;
112
136
  errcod = ip_data("PRINT","%d",&Params.print_lvl,0);
113
137
 
116
140
    Params.ncoord = 0;
117
141
    ip_count(":DBOC:COORDS",&Params.ncoord,0);
118
142
    if (Params.ncoord == 0)
119
 
      done("Keyword COORDS should be a vector of 2-element vectors");
 
143
      done("Keyword COORDS should be a vector of 3-element vectors");
120
144
    Params.coords = new Params_t::Coord_t[Params.ncoord];
121
145
    for(int coord=0; coord<Params.ncoord; coord++) {
122
146
      int nelem = 0;
123
147
      ip_count(":DBOC:COORDS",&nelem,1,coord);
124
 
      if (nelem != 2)
125
 
        done("Keyword COORDS should be a vector of 2-element vectors");
 
148
      if (nelem != 3)
 
149
        done("Keyword COORDS should be a vector of 3-element vectors");
126
150
      errcod = ip_data(":DBOC:COORDS","%d",&Params.coords[coord].index,2,coord,0);
127
151
      if (Params.coords[coord].index < 0 || Params.coords[coord].index >= Molecule.natom*3)
128
152
        done("Keyword COORDS contains an out-of-bounds index");
129
153
      errcod = ip_data(":DBOC:COORDS","%lf",&Params.coords[coord].coeff,2,coord,1);
 
154
      char *tmpstr;
 
155
      errcod = ip_string(":DBOC:COORDS",&tmpstr,2,coord,2);
 
156
      if (!strcmp(tmpstr,"SYMM"))
 
157
        Params.coords[coord].symm = true;
 
158
      else if (!strcmp(tmpstr,"NONSYMM"))
 
159
        Params.coords[coord].symm = false;
 
160
      else
 
161
        done("Keyword COORDS contains an illegal symmetry specifier");
 
162
      free(tmpstr);
 
163
      Params.coords[coord].atom = Params.coords[coord].index/3;
 
164
      Params.coords[coord].xyz = Params.coords[coord].index%3;
130
165
    }
131
166
  }
132
167
  else {
133
 
    Params.ncoord = 0;
134
 
    Params.coords = NULL;
 
168
    Params.ncoord = 3*Molecule.nuniques;
 
169
    Params.coords = new Params_t::Coord_t[Params.ncoord];
 
170
    int coord = 0;
 
171
    for(int ua=0; ua<Molecule.nuniques; ua++) {
 
172
      int atom = Molecule.ua2a[ua];
 
173
      for(int xyz=0; xyz<3; xyz++) {
 
174
        Params.coords[coord].index = 3*atom + xyz;
 
175
        Params.coords[coord].atom = atom;
 
176
        Params.coords[coord].xyz = xyz;
 
177
        Params.coords[coord].coeff = (double)Molecule.ua_degen[ua];
 
178
        //
 
179
        // Figure out whether the cartesian displacement is symmetric under the symmetry operations of the stabilizer of the molecule
 
180
        // (i.e. intersection of the stabilizers of the nuclei, i.e. the subgroup which leaves the molecule invariant)
 
181
        // If it isn't -- then plus and minus displacements along this coordinate are equivalent
 
182
        //
 
183
        Params.coords[coord].symm = false;
 
184
        for(int G=0; G<Molecule.nirreps; G++) {
 
185
          int G_in_stab = true;
 
186
          for(int a=0; a<Molecule.natom; a++)
 
187
            if ((Molecule.ict[G][a]-1) != a)  // this operation isn't in the stabilizer
 
188
              G_in_stab = false;
 
189
          if (G_in_stab) {
 
190
            int Gxyz = (int) Molecule.cartrep[G][3*xyz+xyz];   // character of G in the basis of the displacement
 
191
            if (Gxyz == -1)
 
192
              Params.coords[coord].symm = true;
 
193
          }
 
194
        }
 
195
        coord++;
 
196
      }
 
197
    }
135
198
  }
136
199
 
137
200
  int isotopes_given = ip_exist(":DBOC:ISOTOPES",0);
151
214
  }
152
215
 
153
216
  Params.disp_per_coord = 2;
 
217
  errcod = ip_data(":DBOC:DISP_PER_COORD","%d",&Params.disp_per_coord,0);
 
218
  if (Params.disp_per_coord != 2 && Params.disp_per_coord != 4)
 
219
    throw std::runtime_error("dboc.cc:parsing() -- disp_per_coord must either be 2 or 4");
154
220
}
155
221
 
156
222
/*--- Open chkpt file and grab molecule info ---*/
158
224
{
159
225
  chkpt_init(PSIO_OPEN_OLD);
160
226
  Molecule.natom = chkpt_rd_natom();
 
227
  Molecule.nuniques = chkpt_rd_num_unique_atom();
161
228
  Molecule.geom = chkpt_rd_geom();
162
229
  Molecule.zvals = chkpt_rd_zvals();
 
230
  Molecule.nirreps = chkpt_rd_nirreps();
 
231
  Molecule.ua2a = chkpt_rd_ua2a();
 
232
  Molecule.ict = chkpt_rd_ict();
 
233
  Molecule.cartrep = chkpt_rd_cartrep();
163
234
  chkpt_close();
164
235
 
165
 
  fprintf(outfile, "  -Reference Geometry:\n");
166
 
  for(int i=0; i < Molecule.natom; i++) {
167
 
    fprintf(outfile, "\n   %1.0f ", Molecule.zvals[i]);
168
 
    for(int j=0; j < 3; j++)
169
 
      fprintf(outfile, "%20.10f  ", Molecule.geom[i][j]);
 
236
  // Compute degeneracy of each unique atom
 
237
  Molecule.ua_degen = init_int_array(Molecule.nuniques);
 
238
  for(int ua=0; ua<Molecule.nuniques; ua++) {
 
239
    int a = Molecule.ua2a[ua];
 
240
    int num_identity_maps = 0;
 
241
    for(int g=0; g<Molecule.nirreps; g++)
 
242
      if (a == (Molecule.ict[g][a]-1))
 
243
        num_identity_maps++;
 
244
    Molecule.ua_degen[ua] = Molecule.nirreps/num_identity_maps;
170
245
  }
171
 
  fprintf(outfile, "\n\n");
172
246
}
173
247
 
174
 
double eval_dboc()
 
248
 
 
249
// Returns an array of atomic masses, in a.u. (not in a.m.u.)
 
250
double* get_atomic_masses()
175
251
{
176
252
  double* atomic_mass;
177
253
 
178
254
  //
179
255
  // Convert isotope labels into atomic masses (in a.m.u.)
180
256
  //
 
257
 
181
258
  // Check number of atoms vs. number of isotope labels given in input
182
259
  if (Params.isotopes) {
183
260
    if (Molecule.natom != Params.nisotope)
184
 
      done("Number of atoms in molecule does not match the number of isotope labels in input.dat");
 
261
      done("Number of atoms in molecule does not match the number of isotope labels in input file");
185
262
    else {
186
263
      atomic_mass = new double[Molecule.natom];
187
264
      for(int atom=0; atom<Molecule.natom; atom++) {
207
284
  for(int atom=0; atom<Molecule.natom; atom++)
208
285
    atomic_mass[atom] /= _au2amu;
209
286
 
 
287
  return atomic_mass;
 
288
}
 
289
 
 
290
// Remove PSI3 prefixed file
 
291
void
 
292
remove_psi3_file(char* name)
 
293
{
 
294
  char* execstr = (char*) malloc( sizeof(char) * (strlen("/bin/rm -f ") +
 
295
                                                  strlen(psi_file_prefix) +
 
296
                                                  strlen(name) + 2)
 
297
                                );
 
298
  sprintf(execstr,"/bin/rm -f %s.%s",psi_file_prefix,name);
 
299
  system(execstr);
 
300
  free(execstr);
 
301
}
 
302
 
 
303
// Removes files left by DETCI/DETCAS
 
304
void
 
305
clean_detci_mess()
 
306
{
 
307
  remove_psi3_file("file14.dat");
 
308
  remove_psi3_file("thetas.dat");
 
309
  remove_psi3_file("diis.dat");
 
310
  remove_psi3_file("orbs.dat");
 
311
}
 
312
 
 
313
//
 
314
// This function will run psi for the first displaced point along the coordinate,
 
315
// save the basis set information, and save CI vector, if necessary
 
316
//
 
317
// NOTE: It's assumed that the first displacement is always by -delta
 
318
//
 
319
void run_psi_firstdisp(int disp)
 
320
{
 
321
  system("psiclean");
 
322
  if (!strcmp(Params.wfn,"DETCI") || !strcmp(Params.wfn,"DETCAS")) {
 
323
    clean_detci_mess();
 
324
  }
 
325
  char *inputcmd = new char[80];
 
326
  sprintf(inputcmd,"input --geomdat %d",disp);
 
327
  int errcod = system(inputcmd);
 
328
  if (errcod) {
 
329
    done("input failed");
 
330
  }
 
331
  errcod = system("psi3 --dboc --noinput --messy");
 
332
  if (errcod) {
 
333
    done("psi3 failed");
 
334
  }
 
335
  delete[] inputcmd;
 
336
 
 
337
  //
 
338
  // Read in the "-delta" displaced basis
 
339
  //
 
340
  chkpt_init(PSIO_OPEN_OLD);
 
341
  BasisSets[MinusDelta] = new BasisSet(PSIF_CHKPT);
 
342
  // rotate the geometry back to the original frame and change shell centers
 
343
  double **rref_m = chkpt_rd_rref();
 
344
  double **geom_m = chkpt_rd_geom();
 
345
  chkpt_close();
 
346
  double **geom_m_ref = block_matrix(Molecule.natom,3);
 
347
  mmult(geom_m,0,rref_m,0,geom_m_ref,0,Molecule.natom,3,3,0);
 
348
  free_block(geom_m);
 
349
  free_block(rref_m);
 
350
  for(int a=0; a<Molecule.natom; a++)
 
351
    BasisSets[MinusDelta]->set_center(a,geom_m_ref[a]);
 
352
    
 
353
  // Read in the "-delta" displaced HF wavefunction
 
354
  HFVectors[MinusDelta] = new HFWavefunction();
 
355
 
 
356
  // For CI method rename the saved wave function
 
357
  if (!strcmp(Params.wfn,"DETCI") || !strcmp(Params.wfn,"DETCAS")) {
 
358
    SlaterDetVector *vec;
 
359
    slaterdetvector_read(PSIF_CIVECT,"CI vector",&vec);
 
360
    slaterdetvector_write(PSIF_CIVECT,CI_Vector_Labels[0],vec);
 
361
    slaterdetvector_delete_full(vec);
 
362
 
 
363
    clean_detci_mess();
 
364
  }
 
365
}
 
366
 
 
367
void run_psi_otherdisp(int disp)
 
368
{
 
369
  int disp_coord = (disp-1)%Params.disp_per_coord;
 
370
  int coord = (disp-1)/Params.disp_per_coord;
 
371
  int symm = Params.coords[coord].symm;
 
372
 
 
373
  // If this is an odd displacement means it's a plus displacement. If plus and minus
 
374
  // displacements are equivalent for this coord (symm == true) then an equivalent
 
375
  // displacement has JUST been computed. No need to compute the wave function
 
376
  // again -- just run input, get the rref, and save the wave function.
 
377
  if (!symm || (symm && disp_coord%2 == 0)) {
 
378
    char *inputcmd = new char[80];
 
379
    sprintf(inputcmd,"input --geomdat %d",disp);
 
380
    int errcod = system(inputcmd);
 
381
    if (errcod) {
 
382
      done("input failed");
 
383
    }
 
384
    errcod = system("psi3 --dboc --noinput --messy");
 
385
    if (errcod) {
 
386
      done("psi3 failed");
 
387
    }
 
388
    delete[] inputcmd;
 
389
 
 
390
    // Read in the new displaced HF wavefunction
 
391
    HFVectors[disp_coord] = new HFWavefunction();
 
392
  }
 
393
  else {
 
394
 
 
395
    // Read in the previous displaced HF wavefunction
 
396
    HFVectors[disp_coord] = new HFWavefunction();
 
397
 
 
398
    // only need to update rref
 
399
    char *inputcmd = new char[80];
 
400
    sprintf(inputcmd,"input --savemos --geomdat %d",disp);
 
401
    int errcod = system(inputcmd);
 
402
    if (errcod) {
 
403
      done("input failed");
 
404
    }
 
405
    delete[] inputcmd;
 
406
 
 
407
    chkpt_init(PSIO_OPEN_OLD);
 
408
    double** rref_p = chkpt_rd_rref();
 
409
    chkpt_close();
 
410
    HFVectors[disp_coord]->set_rref(rref_p);
 
411
    free_block(rref_p);
 
412
  }
 
413
 
 
414
  // For CI method rename the saved wave function
 
415
  if (!strcmp(Params.wfn,"DETCI") || !strcmp(Params.wfn,"DETCAS")) {
 
416
    SlaterDetVector *vec;
 
417
    slaterdetvector_read(PSIF_CIVECT,"CI vector",&vec);
 
418
    slaterdetvector_write(PSIF_CIVECT,CI_Vector_Labels[disp_coord],vec);
 
419
    slaterdetvector_delete_full(vec);
 
420
 
 
421
    clean_detci_mess();
 
422
  }
 
423
}
 
424
 
 
425
//
 
426
// This function creates basis set objects for each displacement
 
427
//
 
428
void init_basissets(Params_t::Coord_t* coord)
 
429
{
 
430
  // BasisSetP1
 
431
  int atom = coord->index/3;
 
432
  int xyz = coord->index%3;
 
433
  double AplusD[3];
 
434
  for(int i=0; i<3; i++)
 
435
    AplusD[i] = BasisSets[MinusDelta]->get_center(atom, i);
 
436
  AplusD[xyz] += 2.0*Params.delta;
 
437
  BasisSets[PlusDelta] = new BasisSet(*BasisSets[MinusDelta]);
 
438
  BasisSets[PlusDelta]->set_center(atom,AplusD);
 
439
  
 
440
  if (Params.disp_per_coord == 4) {
 
441
    // BasisSetP2
 
442
    double Aplus2D[3];
 
443
    for(int i=0; i<3; i++)
 
444
      Aplus2D[i] = AplusD[i];
 
445
    Aplus2D[xyz] += Params.delta;
 
446
    BasisSets[Plus2Delta] = new BasisSet(*BasisSets[MinusDelta]);
 
447
    BasisSets[Plus2Delta]->set_center(atom,Aplus2D);
 
448
 
 
449
    // BasisSetM2
 
450
    double Aminus2D[3];
 
451
    for(int i=0; i<3; i++)
 
452
      Aminus2D[i] = Aplus2D[i];
 
453
    Aminus2D[xyz] -= 4.0*Params.delta;
 
454
    BasisSets[Minus2Delta] = new BasisSet(*BasisSets[MinusDelta]);
 
455
    BasisSets[Minus2Delta]->set_center(atom,Aminus2D);
 
456
  }
 
457
}
 
458
 
 
459
//
 
460
// This function destroys basis set objects for each displacement
 
461
//
 
462
void delete_basissets(Params_t::Coord_t* coord)
 
463
{
 
464
  BasisSets[MinusDelta]->~BasisSet();
 
465
  BasisSets[PlusDelta]->~BasisSet();
 
466
  if (Params.disp_per_coord == 4) {
 
467
    BasisSets[Minus2Delta]->~BasisSet();
 
468
    BasisSets[Plus2Delta]->~BasisSet();
 
469
  }
 
470
}
 
471
 
 
472
//
 
473
// This function destroys HFWavefunction objects for each displacement
 
474
//
 
475
void delete_hfwfns(Params_t::Coord_t* coord)
 
476
{
 
477
  HFVectors[MinusDelta]->~HFWavefunction();
 
478
  HFVectors[PlusDelta]->~HFWavefunction();
 
479
  if (Params.disp_per_coord == 4) {
 
480
    HFVectors[Minus2Delta]->~HFWavefunction();
 
481
    HFVectors[Plus2Delta]->~HFWavefunction();
 
482
  }
 
483
}
 
484
 
 
485
double eval_dboc()
 
486
{
 
487
 
 
488
  double* atomic_mass = get_atomic_masses();  
210
489
  const int ndisp = Params.disp_per_coord * Params.ncoord;
211
490
  double E_dboc = 0.0;
212
491
 
213
 
  for(int disp=1; disp<=ndisp; disp++) {
214
 
    Params_t::Coord_t* coord = &(Params.coords[(disp-1)/2]);
 
492
  for(int disp=1; disp<=ndisp;) {
 
493
    int current_coord = (disp-1)/Params.disp_per_coord;
 
494
    Params_t::Coord_t* coord = &(Params.coords[current_coord]);
215
495
    int atom = coord->index/3;
216
496
    int xyz = coord->index%3;
 
497
    const bool symm = coord->symm;
217
498
 
218
 
    char *inputcmd = new char[80];
219
 
#if USE_INPUT_S
220
 
    sprintf(inputcmd,"input --geomdat %d --noreorient --nocomshift",disp);
221
 
#else
222
 
    sprintf(inputcmd,"input --geomdat %d",disp);
223
 
#endif
224
 
    int errcod = system(inputcmd);
225
 
    if (errcod) {
226
 
      done("input failed");
227
 
    }
228
 
    errcod = system("psi3 --dboc --noinput --messy");
229
 
    if (errcod) {
230
 
      done("psi3 failed");
231
 
    }
 
499
    run_psi_firstdisp(disp);
232
500
    disp++;
233
 
 
234
 
    //
235
 
    // Read in the "-" displaced basis
236
 
    //
237
 
    chkpt_init(PSIO_OPEN_OLD);
238
 
    BasisDispM = new BasisSet(PSIF_CHKPT);
239
 
    // rotate the geometry back to the original frame and change shell centers
240
 
    double **rref_m = chkpt_rd_rref();
241
 
    double **geom_m = chkpt_rd_geom();
242
 
    chkpt_close();
243
 
    double **geom_m_ref = block_matrix(Molecule.natom,3);
244
 
    mmult(geom_m,0,rref_m,0,geom_m_ref,0,Molecule.natom,3,3,0);
245
 
    free_block(geom_m);
246
 
    free_block(rref_m);
247
 
    for(int a=0; a<Molecule.natom; a++)
248
 
      BasisDispM->set_center(a,geom_m_ref[a]);
249
 
 
250
 
    // For CI method rename the saved wave function
251
 
    if (!strcmp(Params.wfn,"DETCI") || !strcmp(Params.wfn,"DETCAS")) {
252
 
      SlaterDetVector *vec;
253
 
      slaterdetvector_read(PSIF_CIVECT,"CI vector",&vec);
254
 
      slaterdetvector_write(PSIF_CIVECT,"Old CI vector",vec);
255
 
      slaterdetvector_delete_full(vec);
256
 
    }
257
 
 
258
 
    // obtain "+"-displaced basis by shiting the "active" center in the
259
 
    // "-"-displaced center by twice the displacement size
260
 
    double AplusD[3];
261
 
    AplusD[0] = geom_m_ref[atom][0];
262
 
    AplusD[1] = geom_m_ref[atom][1];
263
 
    AplusD[2] = geom_m_ref[atom][2];
264
 
    AplusD[xyz] += 2*Params.delta;
265
 
    BasisDispP = new BasisSet(*BasisDispM);
266
 
    BasisDispP->set_center(atom,AplusD);
267
 
 
268
 
#if USE_INPUT_S
269
 
    sprintf(inputcmd,"input --savemos --geomdat %d --noreorient --nocomshift",disp);
270
 
#else
271
 
    sprintf(inputcmd,"input --savemos --geomdat %d",disp);
272
 
#endif
273
 
    errcod = system(inputcmd);
274
 
    if (errcod) {
275
 
      done("input failed");
276
 
    }
277
 
    errcod = system("psi3 --dboc --noinput --messy");
278
 
    if (errcod) {
279
 
      done("psi3 failed");
280
 
    }
281
 
    delete[] inputcmd;
282
 
 
283
 
    double S = eval_derwfn_overlap();
284
 
    double del2 = (S-1.0)/(2.0*Params.delta*Params.delta);
 
501
    for(int i=1;i<Params.disp_per_coord; i++) {
 
502
      run_psi_otherdisp(disp);
 
503
      disp++;
 
504
    }
 
505
    
 
506
    init_basissets(coord);
 
507
 
 
508
    // The - sign comes from the integration by parts
 
509
    double del2 = (-1.0)*eval_derwfn_overlap(symm);
285
510
    int Z = (int)Molecule.zvals[atom];
286
511
    // mass of nucleus = atomic mass - mass of electrons
287
512
    double nuclear_mass = atomic_mass[atom]  - Z;
288
513
    double E_i = -del2/(2.0*nuclear_mass);
289
514
    // multiply by the degeneracy factor
290
515
    E_i *= coord->coeff;
291
 
    if (Params.print_lvl > PrintLevels::print_intro) {
292
 
      fprintf(outfile,"  +- wave function overlap = %25.15lf\n",S);
293
 
      fprintf(outfile,"  <d2/dx2>                 = %25.15lf\n",del2);
294
 
      fprintf(outfile,"  DBOC contribution from cartesian degree of freedom %d = %20.10lf a.u.\n\n",
295
 
              coord->index+1,E_i);
 
516
    if (Params.print_lvl >= PrintLevels::print_params) {
 
517
      char xyz = 'x';  xyz = (coord->xyz == 1) ? 'y' : xyz; xyz = (coord->xyz == 2) ? 'z' : xyz;
 
518
      fprintf(outfile,"  <d2/dx2>                    = %25.15lf\n", del2);
 
519
      fprintf(outfile,"  DBOC contribution           = %20.10lf a.u.\n\n", E_i);
296
520
      fflush(outfile);
297
521
    }
298
522
    E_dboc += E_i;
303
527
      psio_close(PSIF_CIVECT,0);
304
528
    }
305
529
    // it's safe to do psiclean now
306
 
    errcod = system("psiclean");
 
530
    int errcod = system("psiclean");
307
531
    if (errcod) {
308
532
      done("psiclean");
309
533
    }
310
534
 
 
535
    delete_basissets(coord);
 
536
    delete_hfwfns(coord);
 
537
 
311
538
  }
312
539
 
313
540
  return E_dboc;
346
573
  // Psi modules called by dboc should write to a different output file
347
574
  // reset the value of PSI_OUTPUT for the duration of this run
348
575
  orig_psi_output_env = getenv("PSI_OUTPUT");
 
576
  char* ofname = (char*) malloc(strlen(psi_ofname())+1+strlen("dboc.findif.out"));
 
577
  sprintf(ofname, "%s.dboc.findif.out", psi_ofname());
349
578
#if HAVE_PUTENV
350
 
  char* tmpstr2 = strdup("PSI_OUTPUT=dboc.findif.out");
 
579
  char* tmpstr2 = (char *) malloc(12+strlen(ofname));
 
580
  sprintf(tmpstr2, "PSI_OUTPUT=%s", ofname);
351
581
  putenv(tmpstr2);
352
582
#elif HAVE_SETENV
353
 
  setenv("PSI_OUTPUT","dboc.findif.out",1);
 
583
  setenv("PSI_OUTPUT",ofname,1);
354
584
#else
355
585
#error "Have neither putenv nor setenv. Something must be very broken on this system."
356
586
#endif