← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/detci/check_energy.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/detci/check_energy.cc
 
1
/*! 
 
2
  \file
 
3
  \ingroup DETCI
 
4
  \brief Check the SCF energy
 
5
*/
 
6
#include <cstdio>
 
7
#include <cmath>
 
8
 
 
9
namespace psi { namespace detci {
 
10
 
 
11
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
 
12
 
 
13
extern int *ioff ;
 
14
 
 
15
void scf_energy(double *H, double *TE, double *energy_1, double *energy_2, 
 
16
      double *energy_e, int *docc, int *frozen_docc, int fzc_flag, 
 
17
      int nirreps, int *reorder, int *opi);
 
18
 
 
19
/*!
 
20
** check_energy(): check the SCF energy by calculating it from the two-electr.
 
21
**    integrals in the MO basis
 
22
** 
 
23
** \param  H         =  lwr tri of one-electron integrals matrix (MO basis)
 
24
** \param twoel_ints =  two electron integrals (lexically indexed, MO basis)
 
25
** \param   nocc     =  num occupied orbitals (assume closed shell case) and
 
26
**                      exclude frozen core
 
27
** \param   escf     =  scf energy to compare to
 
28
** \param   enuc     =  nuclear repulsion energy
 
29
** \param   efzc     =  frozen core energy
 
30
** \param   nirreps  =  number of irreps 
 
31
** \param   reorder  =  reordering array for Pitzer->CI ordering
 
32
** \param   opi      =  orbs per irrep in Pitzer ordering
 
33
** \param   outfile  =  file to write output to
 
34
**
 
35
** Returns: the computed SCF energy
 
36
** \ingroup DETCI
 
37
*/
 
38
double check_energy(double *H, double *twoel_ints, int *docc, int *frozen_docc,
 
39
      int fzc_flag, double escf, double enuc, double efzc, 
 
40
      int nirreps, int *reorder, int *opi, int print_lvl, FILE *outfile)
 
41
{
 
42
   double energy_1 ;     /* one-electron energy */
 
43
   double energy_2 ;     /* two-electron energy */
 
44
   double energy_e ;     /* total electronic energy */
 
45
 
 
46
   scf_energy(H, twoel_ints, &energy_1, &energy_2, &energy_e, docc,
 
47
      frozen_docc, fzc_flag, nirreps, reorder, opi);
 
48
 
 
49
   if (print_lvl) {
 
50
     fprintf(outfile,"\nCheck SCF Energy from 1- and 2-electron integrals\n\n");
 
51
     fprintf(outfile,"SCF Energy (ref):          %16.10lf\n", escf) ;
 
52
     fprintf(outfile,"Nuclear repulsion energy:  %16.10lf\n", enuc) ;
 
53
     fprintf(outfile,"One-electron energy:       %16.10lf\n", energy_1) ;
 
54
     fprintf(outfile,"Two-electron energy:       %16.10lf\n", energy_2) ;
 
55
     fprintf(outfile,"Frozen core energy:        %16.10lf\n", efzc) ;
 
56
     fprintf(outfile,"Total electronic energy:   %16.10lf\n", energy_e+efzc) ;
 
57
     fprintf(outfile,"Total SCF energy:          %16.10lf\n", enuc + 
 
58
        energy_e + efzc) ;
 
59
    
 
60
     if (fabs(enuc + efzc + energy_e - escf) > 0.00000001) {
 
61
        fprintf(outfile, 
 
62
           "\n*** Calculated Energy Differs from SCF Energy in CHKPT ! ***\n") ;
 
63
        }
 
64
   }
 
65
 
 
66
   return(enuc+efzc+energy_e); 
 
67
}   
 
68
 
 
69
 
 
70
/*!
 
71
** scf_energy(): Function calculates the SCF energy from the one- and
 
72
**   two-electron integrals in MO form (closed-shell case).
 
73
**
 
74
** David Sherrill, Sept 1993
 
75
**
 
76
** \param H        = Matrix of one-electron integrals in MO basis (lwr triangle)
 
77
** \param TE       = Two-electron integrals in MO basis, stored in 
 
78
**                   ijkl-indexed array
 
79
** \param energy_1 = pointer to hold one-electron energy
 
80
** \param energy_2 = pointer to hold two-electron energy
 
81
** \param energy_e = pointer to hold total electronic energy (sum of two 
 
82
**                   terms above)
 
83
** \param docc     = array of doubly-occupied orbitals per irrep
 
84
** \param frozen_docc = array of frozen doubly-occupied orbitals per irrep
 
85
** \param fzc_flag = remove explicit consideration of frozen core orbitals ?
 
86
** \param nirreps  = number of irreps 
 
87
** \param reorder  = reordering array Pitzer->CI order
 
88
** \param opi      = orbitals per irrep
 
89
**
 
90
** Returns: none
 
91
**
 
92
** \ingroup DETCI
 
93
*/
 
94
 
 
95
void scf_energy(double *H, double *TE, double *energy_1, double *energy_2, 
 
96
      double *energy_e, int *docc, int *frozen_docc, int fzc_flag, 
 
97
      int nirreps, int *reorder, int *opi)
 
98
{
 
99
   int irrep, irrep2, d, d2, offset, offset2, ndoc, ndoc2, nfzc, nfzc2, totfzc;
 
100
   int i, j;
 
101
   int ii, jj, iijj, ij, ijij, iiii;
 
102
 
 
103
   *energy_1 = *energy_2 = *energy_e = 0.0;
 
104
 
 
105
   totfzc=0;
 
106
   if (fzc_flag) {
 
107
     for (irrep=0; irrep<nirreps; irrep++) {
 
108
       totfzc += frozen_docc[irrep];
 
109
     }
 
110
   }
 
111
 
 
112
   for (irrep=0,offset=0; irrep<nirreps; irrep++) {
 
113
     if (irrep>0) offset += opi[irrep-1];
 
114
     ndoc = docc[irrep];
 
115
     if (fzc_flag) {
 
116
       nfzc = frozen_docc[irrep];
 
117
       ndoc -= nfzc;
 
118
     }
 
119
     else nfzc=0;
 
120
     for (d=offset+nfzc; d<ndoc+offset+nfzc; d++) {
 
121
       i = reorder[d]-totfzc;
 
122
       ii = ioff[i] + i;
 
123
       iiii = ioff[ii] + ii;
 
124
       *energy_1 += 2.0 * H[ii];
 
125
       *energy_2 += TE[iiii];
 
126
       
 
127
       for (irrep2=0,offset2=0; irrep2<=irrep; irrep2++) {
 
128
         if (irrep2>0) offset2 += opi[irrep2-1];
 
129
         ndoc2 = docc[irrep2];
 
130
         if (fzc_flag) {
 
131
           nfzc2 = frozen_docc[irrep2];
 
132
           ndoc2 -= nfzc2;
 
133
         }
 
134
         else nfzc2=0;
 
135
         
 
136
         for (d2=offset2+nfzc2; d2<ndoc2+offset2+nfzc2 && d2<d; d2++) {
 
137
           j = reorder[d2]-totfzc;
 
138
           jj = ioff[j] + j;
 
139
           iijj = INDEX(ii,jj);
 
140
           ij = INDEX(i,j);
 
141
           ijij = ioff[ij] + ij;
 
142
           *energy_2 += 4.0 * TE[iijj] - 2.0 * TE[ijij];
 
143
         }
 
144
       }
 
145
     }
 
146
   }
 
147
   
 
148
   *energy_e = *energy_1 + *energy_2;
 
149
}
 
150
 
 
151
}} // namespace psi::detci
 
152