← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/optking/disp_freq_energy_cart.cc

  • 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:
src/bin/optking/disp_freq_energy_cart.cc
1
 
/* DISP_FREQ_ENERGY_CART makes symmetry-adapted cartesian coordinates and displaces
2
 
   along them to setup frequencies by energy points */
3
 
 
4
 
#if HAVE_CMATH
5
 
# include <cmath>
6
 
#else
7
 
# include <math.h>
8
 
#endif
9
 
 
10
 
extern "C" {
11
 
#include <stdio.h>
 
1
/*! \file
 
2
    \ingroup OPTKING
 
3
    \brief 
 
4
 DISP_FREQ_ENERGY_CART makes symmetry-adapted cartesian coordinates and displaces
 
5
   along them to setup frequencies by energy points
 
6
Formulas for finite differences
 
7
3-point - diagonal
 
8
  O(1/h^2): [f(1,0) + f(-1,0) - 2f(0,0)]/(h^2)
 
9
   [which is same as : [f(2,0) + f(-2,0) - 2f(0,0)]/(4h^2)
 
10
3-point - off-diagonal
 
11
  O(1/h^2): old-way: [f(1,1) - f(1,-1) - f(-1,1) + f(-1,-1)]/(4h^2)
 
12
  O(1/h^2): new-way: [f(1,1)+f(-1,-1)+2f(0,0) -f(1,0) -f(-1,0) -f(0,1) -f(0,-1)]/(2h^2)
 
13
 
 
14
5-point formula - diagonal
 
15
  O(1/h^4): [-f(-2,0) + 16f(-1,0) + 16f(1,0) - f(2,0) - 30f(0,0)] / (12h^2)
 
16
5-point formula - off-diagonal
 
17
  O(1/h^4): old-way [ 1f(-2,-2) - 8f(-1,-2) + 8f(+1,-2) - 1f(+2,-2) - 8f(-2,-1)
 
18
  + 64f(-1,-1) - 64f(+1,-1) + 8f(+2,-1) + 8f(-2,+1) - 64f(-1,+1) + 64f(+1,+1)
 
19
 - 8f(+2,+1) - 1f(-2,+2) + 8f(-1,+2) - 8f(+1,+2) + 1f(+2,+2) / (144h^2)
 
20
  O(1/h^4): new-way [ -1f(-1,-2) - 1f(-2,-1) + 9f(-1,-1) - 1f(+1,-1)
 
21
    - 1f(-1,1) + 9f(+1,+1) - 1f(+2,+1) - 1f(1,2)
 
22
    + 1f(-2,0) - 7f(-1,0)  - 7f(+1,0) + 1f(+2,0)
 
23
    + 1f(0,-2) - 7f(0,-1)  - 7f(0,+1) + 1f(0,+2) + 12f(0,0)]/(12h^2)
 
24
*/
 
25
 
 
26
#include <cmath>
 
27
#include <cstdio>
12
28
#include <libchkpt/chkpt.h>
13
 
#include <stdlib.h>
14
 
#include <string.h>
15
 
#include <ctype.h>
 
29
#include <cstdlib>
 
30
#include <cstring>
 
31
#include <cctype>
16
32
#include <libciomr/libciomr.h>
17
33
#include <libqt/qt.h>
18
34
#include <libipv1/ip_lib.h>
19
35
#include <physconst.h>
20
36
#include <libpsio/psio.h>
21
37
#include <psifiles.h>
22
 
}
23
38
 
24
39
#define EXTERN
25
40
#include "opt.h"
29
44
#include "salc.h"
30
45
#include "bond_lengths.h"
31
46
 
 
47
namespace psi { namespace optking {
 
48
 
32
49
extern int get_irrep_xyz( double **cartrep, int xyz);
33
50
extern int check_coordinates(int natom, double *coord, double *masses, double *Zvals,
34
51
    int *ndisp_small, double ***disp_small);
171
188
    }
172
189
  }
173
190
  nconstraints = cnt+1;
174
 
  /*
175
 
  if (print) {
176
 
    fprintf(outfile,"Normalized constraints\n");
177
 
    print_mat(constraints,nconstraints,3*natom,outfile);
178
 
  }
179
 
  */
180
191
 
181
192
  /* Orthogonalize rotations and translations exactly-is this ever necessary?*/
182
193
  constraints_ortho = block_matrix(nconstraints,3*natom);
187
198
    for (j=0; j<3*natom; ++j)
188
199
      constraints[i][j] = constraints_ortho[i][j];
189
200
  free_block(constraints_ortho);
190
 
  /*
191
 
  if (print) {
192
 
    fprintf(outfile,"Orthogonal constraints\n");
193
 
    print_mat(constraints,nconstraints,3*natom,outfile);
194
 
  }
195
 
  */
196
201
 
197
202
  /**** Form symmetry-adapted cartesian vectors ****/
198
203
  salc_orig = (double ***) malloc(nirreps*sizeof(double **));
290
295
  if (print) {
291
296
    for (irrep=0; irrep < nirreps; ++irrep)
292
297
      if (nsalc[irrep]) {
293
 
        fprintf(outfile,"Projected Cartesian Displacements, irrep %d\n", irrep);
 
298
        fprintf(outfile,"Symmetry Adapted Cartesian Coordinates, irrep %d\n", irrep);
294
299
        print_mat(salc[irrep], nsalc[irrep], 3*natom, outfile);
295
300
      }
296
301
  }
298
303
  check_coordinates(natom,coord,masses,Zvals,nsalc,salc);
299
304
 
300
305
  /* determine number of displacements */
301
 
  /* 1/h^2 * (f(1,0) + f(-1,0) + 2f(0,0) ) */
302
306
  ndisp_all = 0;
303
307
  ndisp = init_int_array(nirreps);
304
308
  for (irrep=0; irrep<nirreps; ++irrep) {
319
323
 
320
324
    /* off-diagonal displacements */
321
325
    if (optinfo.points == 3)
322
 
      ndisp[irrep] += 4 * nsalc[irrep] * (nsalc[irrep] - 1) / 2;
 
326
      ndisp[irrep] += 2 * nsalc[irrep] * (nsalc[irrep] - 1) / 2;
323
327
    else if (optinfo.points == 5)
324
 
      ndisp[irrep] += 16 * nsalc[irrep] * (nsalc[irrep] - 1) / 2;
 
328
      ndisp[irrep] += 8 * nsalc[irrep] * (nsalc[irrep] - 1) / 2;
325
329
 
326
330
    ndisp_all += ndisp[irrep];
327
331
  }
381
385
    /* add off-diagonal displacements */
382
386
 
383
387
    for (i=0; i<nsalc[irrep]; ++i) {
384
 
      for (j=i+1; j<nsalc[irrep]; ++j) {
 
388
      for (j=0; j<i; ++j) {
385
389
        if (optinfo.points == 3) {
386
390
          for (k=0; k < 3*natom; ++k) {
387
391
            geoms[irrep][cnt+0][k] += ( + salc[irrep][i][k] + salc[irrep][j][k] )
388
392
              * optinfo.disp_size / sqrt(masses[k]);
389
 
            geoms[irrep][cnt+1][k] += ( + salc[irrep][i][k] - salc[irrep][j][k] )
390
 
              * optinfo.disp_size / sqrt(masses[k]);
391
 
            geoms[irrep][cnt+2][k] += ( - salc[irrep][i][k] + salc[irrep][j][k] )
392
 
              * optinfo.disp_size / sqrt(masses[k]);
393
 
            geoms[irrep][cnt+3][k] += ( - salc[irrep][i][k] - salc[irrep][j][k] )
 
393
            geoms[irrep][cnt+1][k] += ( - salc[irrep][i][k] - salc[irrep][j][k] )
394
394
              * optinfo.disp_size / sqrt(masses[k]);
395
395
          }
396
 
          cnt += 4;
 
396
          cnt += 2;
397
397
        }
398
398
        else if (optinfo.points == 5) {
399
399
          for (k=0; k < 3*natom; ++k) {
400
 
            geoms[irrep][cnt+0][k] += ( - 2.0 * salc[irrep][i][k] - 2.0 * salc[irrep][j][k] )
401
 
              * optinfo.disp_size / sqrt(masses[k]);
402
 
            geoms[irrep][cnt+1][k] += ( - 1.0 * salc[irrep][i][k] - 2.0 * salc[irrep][j][k] )
403
 
              * optinfo.disp_size / sqrt(masses[k]);
404
 
            geoms[irrep][cnt+2][k] += ( + 1.0 * salc[irrep][i][k] - 2.0 * salc[irrep][j][k] )
405
 
              * optinfo.disp_size / sqrt(masses[k]);
406
 
            geoms[irrep][cnt+3][k] += ( + 2.0 * salc[irrep][i][k] - 2.0 * salc[irrep][j][k] )
407
 
              * optinfo.disp_size / sqrt(masses[k]);
408
 
            geoms[irrep][cnt+4][k] += ( - 2.0 * salc[irrep][i][k] - 1.0 * salc[irrep][j][k] )
409
 
              * optinfo.disp_size / sqrt(masses[k]);
410
 
            geoms[irrep][cnt+5][k] += ( - 1.0 * salc[irrep][i][k] - 1.0 * salc[irrep][j][k] )
411
 
              * optinfo.disp_size / sqrt(masses[k]);
412
 
            geoms[irrep][cnt+6][k] += ( + 1.0 * salc[irrep][i][k] - 1.0 * salc[irrep][j][k] )
413
 
              * optinfo.disp_size / sqrt(masses[k]);
414
 
            geoms[irrep][cnt+7][k] += ( + 2.0 * salc[irrep][i][k] - 1.0 * salc[irrep][j][k] )
415
 
              * optinfo.disp_size / sqrt(masses[k]);
416
 
            geoms[irrep][cnt+8][k] += ( - 2.0 * salc[irrep][i][k] + 1.0 * salc[irrep][j][k] )
417
 
              * optinfo.disp_size / sqrt(masses[k]);
418
 
            geoms[irrep][cnt+9][k] += ( - 1.0 * salc[irrep][i][k] + 1.0 * salc[irrep][j][k] )
419
 
              * optinfo.disp_size / sqrt(masses[k]);
420
 
            geoms[irrep][cnt+10][k] += ( + 1.0 * salc[irrep][i][k] + 1.0 * salc[irrep][j][k] )
421
 
              * optinfo.disp_size / sqrt(masses[k]);
422
 
            geoms[irrep][cnt+11][k] += ( + 2.0 * salc[irrep][i][k] + 1.0 * salc[irrep][j][k] )
423
 
              * optinfo.disp_size / sqrt(masses[k]);
424
 
            geoms[irrep][cnt+12][k] += ( - 2.0 * salc[irrep][i][k] + 2.0 * salc[irrep][j][k] )
425
 
              * optinfo.disp_size / sqrt(masses[k]);
426
 
            geoms[irrep][cnt+13][k] += ( - 1.0 * salc[irrep][i][k] + 2.0 * salc[irrep][j][k] )
427
 
              * optinfo.disp_size / sqrt(masses[k]);
428
 
            geoms[irrep][cnt+14][k] += ( + 1.0 * salc[irrep][i][k] + 2.0 * salc[irrep][j][k] )
429
 
              * optinfo.disp_size / sqrt(masses[k]);
430
 
            geoms[irrep][cnt+15][k] += ( + 2.0 * salc[irrep][i][k] + 2.0 * salc[irrep][j][k] )
 
400
            geoms[irrep][cnt+0][k] += ( - 1.0 * salc[irrep][i][k] - 2.0 * salc[irrep][j][k] )
 
401
              * optinfo.disp_size / sqrt(masses[k]);
 
402
            geoms[irrep][cnt+1][k] += ( - 2.0 * salc[irrep][i][k] - 1.0 * salc[irrep][j][k] )
 
403
              * optinfo.disp_size / sqrt(masses[k]);
 
404
            geoms[irrep][cnt+2][k] += ( - 1.0 * salc[irrep][i][k] - 1.0 * salc[irrep][j][k] )
 
405
              * optinfo.disp_size / sqrt(masses[k]);
 
406
            geoms[irrep][cnt+3][k] += ( + 1.0 * salc[irrep][i][k] - 1.0 * salc[irrep][j][k] )
 
407
              * optinfo.disp_size / sqrt(masses[k]);
 
408
            geoms[irrep][cnt+4][k] += ( - 1.0 * salc[irrep][i][k] + 1.0 * salc[irrep][j][k] )
 
409
              * optinfo.disp_size / sqrt(masses[k]);
 
410
            geoms[irrep][cnt+5][k] += ( + 1.0 * salc[irrep][i][k] + 1.0 * salc[irrep][j][k] )
 
411
              * optinfo.disp_size / sqrt(masses[k]);
 
412
            geoms[irrep][cnt+6][k] += ( + 2.0 * salc[irrep][i][k] + 1.0 * salc[irrep][j][k] )
 
413
              * optinfo.disp_size / sqrt(masses[k]);
 
414
            geoms[irrep][cnt+7][k] += ( + 1.0 * salc[irrep][i][k] + 2.0 * salc[irrep][j][k] )
431
415
              * optinfo.disp_size / sqrt(masses[k]);
432
416
          }
433
 
          cnt += 16;
 
417
          cnt += 8;
434
418
        }
435
419
      }
436
420
    }
437
421
  }
438
422
  free(masses);
439
423
 
 
424
  if(optinfo.external_energies){
 
425
    int counter = 1;
 
426
    FILE *fp_dispcart = fopen("dispcart", "w");
 
427
    // ACS (11/06) Print the displaced and reference geometry for external programs
 
428
    // Use 1 based counting, and remember that the first "displacement" is the reference itself
 
429
    fprintf(fp_dispcart,"%d\n",counter++);
 
430
    for(int n = 0; n < carts.get_natom(); n++){
 
431
      fprintf(fp_dispcart,"%16.10lf %16.10lf %16.10lf\n",coord[3*n],coord[3*n+1],coord[3*n+2]);
 
432
    }
 
433
 
 
434
    for (irrep=0; irrep<nirreps; ++irrep) {
 
435
      for(int disp = 0; disp< ndisp[irrep]; disp++){
 
436
        fprintf(fp_dispcart,"%d\n",counter++);
 
437
        for(int n = 0; n < carts.get_natom(); n++){
 
438
          fprintf(fp_dispcart,"%16.10lf %16.10lf %16.10lf\n",
 
439
            geoms[irrep][disp][3*n],geoms[irrep][disp][3*n+1],geoms[irrep][disp][3*n+2]);
 
440
        }
 
441
      }
 
442
    }
 
443
    fclose(fp_dispcart);
 
444
  }
 
445
 
440
446
  /**** write out info to PSIF_OPTKING, geoms and coordinates in irrep order ****/
441
447
  open_PSIF();
442
448
 
454
460
    if (ndisp[irrep]) {
455
461
      psio_write(PSIF_OPTKING, "OPT: Displaced geometries",
456
462
          (char *) &(geoms[irrep][0][0]), ndisp[irrep]*3*natom*sizeof(double), next, &next);
457
 
fprintf(outfile,"Displaced Geometries of irrep %d.\n",irrep);
458
 
print_mat(geoms[irrep],ndisp[irrep],3*natom,outfile);
 
463
//fprintf(outfile,"Displaced Geometries of irrep %d.\n",irrep);
 
464
//print_mat(geoms[irrep],ndisp[irrep],3*natom,outfile);
459
465
    }
460
466
    free_block(geoms[irrep]);
461
467
  }
488
494
  disp_e = init_array(ndisp_all);
489
495
  optinfo.micro_iteration = 0;
490
496
 
491
 
  /*
492
 
  ndisp_all = 8;
493
 
  ndisp[1] = 0;
494
 
  ndisp[2] = 0;
495
 
  ndisp[3] = 0;
496
 
  nsalc_all = 4;
497
 
  nsalc[1] = 0;
498
 
  nsalc[2] = 0;
499
 
  nsalc[3] = 0;
500
 
  */
501
 
 
502
497
  psio_write_entry(PSIF_OPTKING, "OPT: Displaced energies",
503
498
      (char *) &(disp_e[0]), ndisp_all*sizeof(double));
504
499
  psio_write_entry(PSIF_OPTKING, "OPT: Reference geometry",
533
528
  return(ndisp_all);
534
529
}
535
530
 
 
531
}} /* namespace psi::optking */
 
532