← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/cints/Default_Deriv2/deriv2.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/cints/Default_Deriv2/deriv2.cc
 
1
/*! \file deriv2.cc
 
2
    \ingroup CINTS
 
3
    \brief Driver for second derivative integrals.
 
4
*/
 
5
#include<cstdio>
 
6
#include<cstdlib>
 
7
 
 
8
#include<libipv1/ip_lib.h>
 
9
#include<cmath>
 
10
#include<libciomr/libciomr.h>
 
11
#include<libpsio/psio.h>
 
12
#include<libint/libint.h>
 
13
#include<libiwl/iwl.h>
 
14
#include<libqt/qt.h>
 
15
#include<psifiles.h>
 
16
#include <Tools/prints.h>
 
17
 
 
18
#include"defines.h"
 
19
#define EXTERN
 
20
#include"global.h"
 
21
#include <stdexcept>
 
22
#include"moinfo.h"
 
23
#include"compute_scf_opdm.h"
 
24
#include"read_gen_opdm.h"
 
25
#include"enuc_deriv2.h"
 
26
#include"oe_deriv2.h"
 
27
#include"te_deriv2_scf.h"
 
28
#include"file11.h"
 
29
 
 
30
#define DO_NUC 1
 
31
#define DO_OE 1
 
32
#define DO_TE 1
 
33
namespace psi {
 
34
  namespace CINTS {
 
35
  //! Driver for integral second derivatives.
 
36
  void deriv2(void)
 
37
  {
 
38
  int i, j, ij, coord, ntri, natom, size;
 
39
  char *label;
 
40
  double **tmpmat, **my_grad, *outbuf;
 
41
  double **C;
 
42
 
 
43
  /*--- Hessian in the canonical frame ---*/
 
44
  Hess = block_matrix(Molecule.num_atoms*3,Molecule.num_atoms*3);
 
45
 
 
46
  /* Derivative Fock and overlap matrices */
 
47
  F = (double ***) malloc(Molecule.num_atoms * 3 * sizeof(double **));
 
48
  S = (double ***) malloc(Molecule.num_atoms * 3 * sizeof(double **));
 
49
  HDS = (double ***) malloc(Molecule.num_atoms * 3 * sizeof(double **));
 
50
  for(i=0; i < Molecule.num_atoms*3; i++) {
 
51
    F[i] = block_matrix(BasisSet.num_ao, BasisSet.num_ao);
 
52
    S[i] = block_matrix(BasisSet.num_ao, BasisSet.num_ao);
 
53
    HDS[i] = block_matrix(BasisSet.num_ao, BasisSet.num_ao);
 
54
  }
 
55
 
 
56
  init_moinfo();
 
57
  compute_scf_opdm();
 
58
 
 
59
#if DO_NUC
 
60
  enuc_deriv2();
 
61
#endif
 
62
#if DO_OE
 
63
  oe_deriv2();
 
64
#endif
 
65
 
 
66
 
 
67
 
 
68
  /* Multiply the F's by 2 -- accounts for orbital population factor */
 
69
  for(coord=0; coord < Molecule.num_atoms*3; coord++) {
 
70
    for(i=0; i < BasisSet.num_ao; i++)
 
71
      for(j=0; j < BasisSet.num_ao; j++)
 
72
        F[coord][i][j] *= 2.0;
 
73
  }
 
74
 
 
75
#if DO_TE
 
76
  te_deriv2_scf();
 
77
  /* te_deriv2_scf_symm(); */
 
78
#endif
 
79
 
 
80
  /* divide the whole thing by 2 */
 
81
  for(coord=0; coord < Molecule.num_atoms*3; coord++) {
 
82
    for(i=0; i < BasisSet.num_ao; i++)
 
83
      for(j=0; j < BasisSet.num_ao; j++)
 
84
        F[coord][i][j] *= 0.5;
 
85
  }
 
86
 
 
87
  /* symmetrize the F's and S's */
 
88
  for(coord=0; coord < Molecule.num_atoms*3; coord++) {
 
89
    if (UserOptions.print_lvl >= PRINT_OEDERIV) {
 
90
      fprintf(outfile, "AO-basis Overlap Derivs (Pre-Symm) (%d)", coord);
 
91
      print_mat(S[coord],BasisSet.num_ao,BasisSet.num_ao,outfile);
 
92
    }
 
93
 
 
94
    for(i=0; i < BasisSet.num_ao; i++) {
 
95
      for(j=0; j <= i; j++) {
 
96
        if(i!=j) {
 
97
          F[coord][i][j] = F[coord][j][i] = 0.5 * (F[coord][i][j] + F[coord][j][i]);
 
98
          S[coord][i][j] = S[coord][j][i] = 0.5 * (S[coord][i][j] + S[coord][j][i]);
 
99
        }
 
100
      }
 
101
    }
 
102
 
 
103
    if (UserOptions.print_lvl >= PRINT_OEDERIV) {
 
104
      fprintf(outfile, "AO-basis Overlap Derivs (%d)", coord);
 
105
      print_mat(S[coord],BasisSet.num_ao,BasisSet.num_ao,outfile);
 
106
    }
 
107
  }
 
108
 
 
109
  /* Transform Fock and Overlap derivatives to the MO basis */
 
110
  tmpmat = block_matrix(MOInfo.num_mo, BasisSet.num_ao);
 
111
  for(i=0; i < Molecule.num_atoms*3; i++) {
 
112
    C_DGEMM('n','n',MOInfo.num_mo,BasisSet.num_ao,BasisSet.num_ao,1.0,
 
113
            &(MOInfo.scf_evec[0][0][0]),BasisSet.num_ao,&(F[i][0][0]),BasisSet.num_ao,
 
114
            0.0,&(tmpmat[0][0]),BasisSet.num_ao);
 
115
    C_DGEMM('n','t',MOInfo.num_mo,MOInfo.num_mo,BasisSet.num_ao,1.0,
 
116
            &(tmpmat[0][0]),BasisSet.num_ao,&(MOInfo.scf_evec[0][0][0]),BasisSet.num_ao,
 
117
            0.0,&(F[i][0][0]),BasisSet.num_ao);
 
118
 
 
119
    C_DGEMM('n','n',MOInfo.num_mo,BasisSet.num_ao,BasisSet.num_ao,1.0,
 
120
            &(MOInfo.scf_evec[0][0][0]),BasisSet.num_ao,&(S[i][0][0]),BasisSet.num_ao,
 
121
            0.0,&(tmpmat[0][0]),BasisSet.num_ao);
 
122
    C_DGEMM('n','t',MOInfo.num_mo,MOInfo.num_mo,BasisSet.num_ao,1.0,
 
123
            &(tmpmat[0][0]),BasisSet.num_ao,&(MOInfo.scf_evec[0][0][0]),BasisSet.num_ao,
 
124
            0.0,&(S[i][0][0]),BasisSet.num_ao);
 
125
  }
 
126
  free_block(tmpmat);
 
127
 
 
128
  /* write derivative Fock and overlap derivatives to disk */
 
129
  ntri = MOInfo.num_mo * (MOInfo.num_mo + 1)/2;
 
130
  outbuf = init_array(ntri);
 
131
  label = (char *) malloc(PSIO_KEYLEN * sizeof(char));
 
132
  for(i=0; i < PSIO_KEYLEN; i++) label[i] = '\0';
 
133
  for(coord=0; coord < Molecule.num_atoms*3; coord++) {
 
134
 
 
135
    for(i=0, ij=0; i < MOInfo.num_mo; i++) {
 
136
      for(j=0; j <= i; j++, ij++) { /* Lower triangles only */
 
137
        outbuf[ij] = F[coord][i][j];
 
138
      }
 
139
    }
 
140
 
 
141
    sprintf(label, "MO-basis Fock Derivs (%d)", coord);
 
142
    iwl_wrtone(PSIF_OEI, label, ntri, outbuf);
 
143
    if (UserOptions.print_lvl >= PRINT_OEDERIV) {
 
144
      fprintf(outfile,"  -%s\n",label);
 
145
      print_mat(F[coord],MOInfo.num_mo,MOInfo.num_mo,outfile);
 
146
    }
 
147
    for(i=0; i < PSIO_KEYLEN; i++) label[i] = '\0';
 
148
  }
 
149
  for(coord=0; coord < Molecule.num_atoms*3; coord++) {
 
150
 
 
151
    for(i=0, ij=0; i < MOInfo.num_mo; i++) {
 
152
      for(j=0; j <= i; j++, ij++) { /* Lower triangles only */
 
153
        outbuf[ij] = S[coord][i][j];
 
154
      }
 
155
    }
 
156
 
 
157
    sprintf(label, "MO-basis Overlap Derivs (%d)", coord);
 
158
    iwl_wrtone(PSIF_OEI, label, ntri, outbuf);
 
159
    if (UserOptions.print_lvl >= PRINT_OEDERIV) {
 
160
      fprintf(outfile,"  -%s\n",label);
 
161
      print_mat(S[coord],MOInfo.num_mo,MOInfo.num_mo,outfile);
 
162
    }
 
163
    for(i=0; i < PSIO_KEYLEN; i++) label[i] = '\0';
 
164
  }
 
165
  /* Write half-differentiated overlap integrals to disk */
 
166
  size = BasisSet.num_ao * BasisSet.num_ao;
 
167
  for(coord=0; coord < Molecule.num_atoms*3; coord++) {
 
168
 
 
169
    sprintf(label, "AO-basis Half-Diff Overlap (%d)", coord);
 
170
    psio_open(PSIF_OEI, PSIO_OPEN_OLD);
 
171
    psio_write_entry(PSIF_OEI, label, (char *) HDS[coord][0], size*sizeof(double));
 
172
    psio_close(PSIF_OEI, 1);
 
173
    if (UserOptions.print_lvl >= PRINT_OEDERIV) {
 
174
      fprintf(outfile,"  -%s\n",label);
 
175
      print_mat(HDS[coord],BasisSet.num_ao,BasisSet.num_ao,outfile);
 
176
    }
 
177
    for(i=0; i < PSIO_KEYLEN; i++) label[i] = '\0';
 
178
  }
 
179
  free(outbuf);
 
180
  free(label);
 
181
 
 
182
  /* print_atommat("Skeleton contribution to the molecular Hessian (a.u.)",Hess); */
 
183
  natom = Molecule.num_atoms;
 
184
  psio_open(PSIF_DERINFO, PSIO_OPEN_NEW);
 
185
  psio_write_entry(PSIF_DERINFO, "Skeleton Hessian", (char *) Hess[0], natom*3*natom*3*sizeof(double));
 
186
  psio_close(PSIF_DERINFO, 1);
 
187
 
 
188
  for(i=0; i < Molecule.num_atoms*3; i++) {
 
189
    free_block(F[i]);
 
190
    free_block(S[i]);
 
191
    free_block(HDS[i]);
 
192
  }
 
193
  free(F);  free(S);  free(HDS);
 
194
  free_block(Hess);
 
195
  cleanup_moinfo();
 
196
 
 
197
  return;
 
198
}
 
199
 
 
200
};};