← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/dboc/roci.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/dboc/roci.cc
 
1
/*! \file
 
2
 \ingroup DBOC
 
3
 \brief Enter brief description of file here 
 
4
 */
1
5
#include <iostream>
2
 
#include <stdio.h>
3
 
#include <stdlib.h>
4
 
#include <math.h>
5
 
extern "C" {
 
6
#include <cstdio>
 
7
#include <cstdlib>
 
8
#include <cmath>
6
9
#include <libciomr/libciomr.h>
7
10
#include <libqt/qt.h>
8
11
#include <libqt/slaterdset.h>
9
12
#include <psifiles.h>
10
 
}
11
13
#include "moinfo.h"
12
14
#include "float.h"
13
15
#include "linalg.h"
14
16
#include "mo_overlap.h"
15
17
#include "hfwfn.h"
16
18
 
 
19
// Set to 1 to try using unit MO overlap between + and - displacements
 
20
#define TRY_UNIT_OVERLAP 0
 
21
 
17
22
using namespace std;
18
 
 
19
 
// Wrap a,b indices into one composite index assuming S2 symmetry
 
23
namespace psi {
 
24
  namespace dboc {
 
25
    
 
26
    // Wrap a,b indices into one composite index assuming S2 symmetry
20
27
#define INDEX2(a,b) ((a) > (b)) ? ( (((a)*(a+1)) >> 1) + (b) ) : ( (((b)*(b+1)) >> 1) + (a) )
21
 
// Wrap a>=b indices into one composite index assuming S2 symmetry
 
28
    // Wrap a>=b indices into one composite index assuming S2 symmetry
22
29
#define INDEX2_ORD(a,b) ( (((a)*(a+1))/2) + (b) )
23
 
// Wrap a>=b>=c indices into one composite index assuming S3 symmetry
 
30
    // Wrap a>=b>=c indices into one composite index assuming S3 symmetry
24
31
#define INDEX3_ORD(a,b,c) ( ((a)*(((a)+4)*((a)-1)+6)/6) + (((b)*(b+1))/2) + (c) )
25
 
 
26
 
extern MOInfo_t MOInfo;
27
 
extern FILE *outfile;
28
 
extern char *CI_Vector_Labels[MAX_NUM_DISP];
29
 
extern HFWavefunction* HFVectors[MAX_NUM_DISP];
30
 
extern void done(const char *);
31
 
extern void mo_maps(short int**, short int**);
32
 
 
33
 
double eval_roci_derwfn_overlap(DisplacementIndex LDisp, DisplacementIndex RDisp)
34
 
{
35
 
  // Read in CI vectors
36
 
  SlaterDetVector *vecm, *vecp;
37
 
  slaterdetvector_read(PSIF_CIVECT,CI_Vector_Labels[RDisp],&vecm);
38
 
  slaterdetvector_read(PSIF_CIVECT,CI_Vector_Labels[LDisp],&vecp);
39
 
 
40
 
  int nfzc = vecm->sdset->alphastrings->nfzc;
 
32
    
 
33
    extern "C" FILE *outfile;
 
34
 
 
35
    extern char *CI_Vector_Labels[MAX_NUM_DISP];
 
36
    extern HFWavefunction* HFVectors[MAX_NUM_DISP];
 
37
    extern void done(const char *);
 
38
    extern void mo_maps(short int**, short int**);
 
39
 
 
40
    double eval_roci_derwfn_overlap(DisplacementIndex LDisp, DisplacementIndex RDisp)
 
41
    {
 
42
      // Read in CI vectors
 
43
      SlaterDetVector *vecm, *vecp;
 
44
      slaterdetvector_read(PSIF_CIVECT,CI_Vector_Labels[RDisp],&vecm);
 
45
      slaterdetvector_read(PSIF_CIVECT,CI_Vector_Labels[LDisp],&vecp);
 
46
 
 
47
      int nfzc = vecm->sdset->alphastrings->nfzc;
41
48
#if USE_MOINFO
42
 
  int nalpha = MOInfo.nalpha;
43
 
  int nbeta = MOInfo.nbeta;
44
 
#else
45
 
  int nalpha = HFVectors[LDisp]->nalpha();
46
 
  int nbeta = HFVectors[LDisp]->nbeta();
47
 
#endif
48
 
  int nact_a = nalpha - nfzc;
49
 
  int nact_b = nbeta - nfzc;
50
 
 
51
 
  FLOAT **CSC_full = eval_S_alpha(LDisp,RDisp);
52
 
  FLOAT **CSC_a = create_matrix(nalpha,nalpha);
53
 
  FLOAT **CSC_b = create_matrix(nbeta,nbeta);
54
 
  int *tmpintvec = new int[nalpha];
55
 
 
56
 
  // Compute overlap between strings for alpha spin case
57
 
  StringSet *ssetm;
58
 
  ssetm = vecm->sdset->alphastrings;
59
 
  short int* fzc_occ = ssetm->fzc_occ;
60
 
  int nstr_a = ssetm->size;
61
 
  FLOAT **S_a = create_matrix(nstr_a,nstr_a);
62
 
  // Assume the order of strings is the same for - and + displacements
63
 
  for(int jp=0; jp<nstr_a; jp++) {
64
 
    String *str_j = &ssetm->strings[jp];
65
 
    for(int im=0; im<nstr_a; im++) {
66
 
      String *str_i = &ssetm->strings[im];
67
 
 
68
 
      for(int j=0;j<nact_a;j++)
69
 
        for(int i=0;i<nact_a;i++)
70
 
          CSC_a[j+nfzc][i+nfzc] = CSC_full[str_j->occ[j]][str_i->occ[i]];
71
 
 
72
 
      // frozen orbitals need to be mapped to pitzer order manually
73
 
      for(int i=0; i<nfzc; i++) {
74
 
        int ii = fzc_occ[i];
75
 
        for(int j=0; j<nact_a; j++)
76
 
          CSC_a[j+nfzc][i] = CSC_full[str_j->occ[j]][i];
77
 
      }
78
 
 
79
 
      for(int j=0; j<nfzc; j++) {
80
 
        int jj = fzc_occ[j];
81
 
        for(int i=0; i<nact_a; i++)
82
 
          CSC_a[j][i+nfzc] = CSC_full[jj][str_i->occ[i]];
83
 
      }
84
 
 
85
 
      for(int i=0;i<nfzc;i++) {
86
 
        int ii = fzc_occ[i];
87
 
        for(int j=0;j<nfzc;j++)
88
 
          CSC_a[i][j] = CSC_full[ii][fzc_occ[j]];
89
 
      }
90
 
 
91
 
      FLOAT sign;
92
 
      lu_decom(CSC_a, nalpha, tmpintvec, &sign);
93
 
      FLOAT deter1 = 1.0;
94
 
      for(int i=0;i<nalpha;i++)
95
 
        deter1 *= CSC_a[i][i];
96
 
 
97
 
      S_a[jp][im] = sign*deter1;
98
 
    }
99
 
  }
100
 
 
101
 
  // Compute overlap between strings for beta spin case
102
 
  ssetm = vecm->sdset->betastrings;
103
 
  int nstr_b = ssetm->size;
104
 
  FLOAT **S_b = create_matrix(nstr_b,nstr_b);
105
 
  // Assume the order of strings is the same for - and + displacements
106
 
  for(int jp=0; jp<nstr_b; jp++) {
107
 
    String *str_j = &ssetm->strings[jp];
108
 
    for(int im=0; im<nstr_b; im++) {
109
 
      String *str_i = &ssetm->strings[im];
110
 
 
111
 
      for(int j=0;j<nact_b;j++)
112
 
        for(int i=0;i<nact_b;i++)
113
 
          CSC_b[j+nfzc][i+nfzc] = CSC_full[str_j->occ[j]][str_i->occ[i]];
114
 
 
115
 
      // frozen orbitals need to be mapped to pitzer order manually
116
 
      for(int i=0; i<nfzc; i++) {
117
 
        int ii = fzc_occ[i];
118
 
        for(int j=0; j<nact_b; j++)
119
 
          CSC_b[j+nfzc][i] = CSC_full[str_j->occ[j]][ii];
120
 
      }
121
 
 
122
 
      for(int j=0; j<nfzc; j++) {
123
 
        int jj = fzc_occ[j];
124
 
        for(int i=0; i<nact_b; i++)
125
 
          CSC_b[j][i+nfzc] = CSC_full[jj][str_i->occ[i]];
126
 
      }
127
 
 
128
 
      for(int i=0;i<nfzc;i++) {
129
 
        int ii = fzc_occ[i];
130
 
        for(int j=0;j<nfzc;j++)
131
 
          CSC_b[i][j] = CSC_full[ii][fzc_occ[j]];
132
 
      }
133
 
 
134
 
      FLOAT sign;
135
 
      lu_decom(CSC_b, nbeta, tmpintvec, &sign);
136
 
      FLOAT deter1 = 1.0;
137
 
      for(int i=0;i<nbeta;i++)
138
 
        deter1 *= CSC_b[i][i];
139
 
 
140
 
      S_b[jp][im] = sign*deter1;
141
 
    }
142
 
  }
143
 
 
144
 
  // Evaluate total overlap in the highest available precision
145
 
  int ndets = vecm->size;
146
 
  FLOAT S_tot = 0.0;
147
 
  for(int I=0; I<ndets; I++) {
148
 
    SlaterDet *detI = vecm->sdset->dets + I;
149
 
    int Istra = detI->alphastring;
150
 
    int Istrb = detI->betastring;
151
 
    FLOAT cI = vecm->coeffs[I];
152
 
 
153
 
    for(int J=0; J<ndets; J++) {
154
 
      SlaterDet *detJ = vecp->sdset->dets + J;
155
 
      int Jstra = detJ->alphastring;
156
 
      int Jstrb = detJ->betastring;
157
 
      FLOAT cJ = vecp->coeffs[J];
158
 
      
159
 
      FLOAT S = S_a[Jstra][Istra] * S_b[Jstrb][Istrb];
160
 
 
161
 
      FLOAT contrib = cI * S * cJ;
162
 
      S_tot += cI * S * cJ;
163
 
      /* fprintf(outfile,"  %3d %3d %+15.10Le", I, J, cI);
164
 
      fprintf(outfile," %+15.10Le", cJ);
165
 
      fprintf(outfile," %+25.15Le", S);
166
 
      fprintf(outfile," %+25.15Le", contrib);
167
 
      fprintf(outfile," %+25.15Le\n", S_tot); */
168
 
 
169
 
    }
170
 
  }
171
 
  
172
 
 
173
 
  // Cleanup
174
 
  slaterdetvector_delete_full(vecm);
175
 
  slaterdetvector_delete_full(vecp);
176
 
  delete[] tmpintvec;
177
 
  delete_matrix(CSC_a);
178
 
  delete_matrix(CSC_full);
179
 
  delete_matrix(CSC_b);
180
 
  delete_matrix(S_a);
181
 
  delete_matrix(S_b);
182
 
  double S_tot_double = (double) S_tot;
183
 
  return fabs(S_tot_double);
184
 
}
185
 
 
 
49
      extern MOInfo_t MOInfo;
 
50
      int nalpha = MOInfo.nalpha;
 
51
      int nbeta = MOInfo.nbeta;
 
52
#else
 
53
      int nalpha = HFVectors[LDisp]->nalpha();
 
54
      int nbeta = HFVectors[LDisp]->nbeta();
 
55
#endif
 
56
      int nact_a = nalpha - nfzc;
 
57
      int nact_b = nbeta - nfzc;
 
58
 
 
59
#if !TRY_UNIT_OVERLAP
 
60
      FLOAT **CSC_full = eval_S_alpha(LDisp,RDisp);
 
61
#else
 
62
      const int num_mo = HFVectors[LDisp]->num_mo();
 
63
      FLOAT **CSC_full = create_matrix(num_mo, num_mo);
 
64
      for(int i=0; i<num_mo; ++i)
 
65
        CSC_full[i][i] = 1.0;
 
66
#endif
 
67
      FLOAT **CSC_a = create_matrix(nalpha,nalpha);
 
68
      FLOAT **CSC_b = create_matrix(nbeta,nbeta);
 
69
      int *tmpintvec = new int[nalpha];
 
70
 
 
71
      // Compute overlap between strings for alpha spin case
 
72
      StringSet *ssetm;
 
73
      ssetm = vecm->sdset->alphastrings;
 
74
      short int* fzc_occ = ssetm->fzc_occ;
 
75
      int nstr_a = ssetm->size;
 
76
      FLOAT **S_a = create_matrix(nstr_a,nstr_a);
 
77
      // Assume the order of strings is the same for - and + displacements
 
78
      for(int jp=0; jp<nstr_a; jp++) {
 
79
        String *str_j = &ssetm->strings[jp];
 
80
        for(int im=0; im<nstr_a; im++) {
 
81
          String *str_i = &ssetm->strings[im];
 
82
 
 
83
          for(int j=0;j<nact_a;j++)
 
84
            for(int i=0;i<nact_a;i++)
 
85
              CSC_a[j+nfzc][i+nfzc] = CSC_full[str_j->occ[j]][str_i->occ[i]];
 
86
 
 
87
          // frozen orbitals need to be mapped to pitzer order manually
 
88
          for(int i=0; i<nfzc; i++) {
 
89
            int ii = fzc_occ[i];
 
90
            for(int j=0; j<nact_a; j++)
 
91
              CSC_a[j+nfzc][i] = CSC_full[str_j->occ[j]][ii];
 
92
          }
 
93
 
 
94
          for(int j=0; j<nfzc; j++) {
 
95
            int jj = fzc_occ[j];
 
96
            for(int i=0; i<nact_a; i++)
 
97
              CSC_a[j][i+nfzc] = CSC_full[jj][str_i->occ[i]];
 
98
          }
 
99
 
 
100
          for(int i=0;i<nfzc;i++) {
 
101
            int ii = fzc_occ[i];
 
102
            for(int j=0;j<nfzc;j++)
 
103
              CSC_a[i][j] = CSC_full[ii][fzc_occ[j]];
 
104
          }
 
105
 
 
106
          FLOAT sign;
 
107
          lu_decom(CSC_a, nalpha, tmpintvec, &sign);
 
108
          FLOAT deter1 = 1.0;
 
109
          for(int i=0;i<nalpha;i++)
 
110
          deter1 *= CSC_a[i][i];
 
111
 
 
112
          S_a[jp][im] = sign*deter1;
 
113
        }
 
114
      }
 
115
 
 
116
      // Compute overlap between strings for beta spin case
 
117
      ssetm = vecm->sdset->betastrings;
 
118
      int nstr_b = ssetm->size;
 
119
      FLOAT **S_b = create_matrix(nstr_b,nstr_b);
 
120
      // Assume the order of strings is the same for - and + displacements
 
121
      for(int jp=0; jp<nstr_b; jp++) {
 
122
        String *str_j = &ssetm->strings[jp];
 
123
        for(int im=0; im<nstr_b; im++) {
 
124
          String *str_i = &ssetm->strings[im];
 
125
 
 
126
          for(int j=0;j<nact_b;j++)
 
127
            for(int i=0;i<nact_b;i++)
 
128
              CSC_b[j+nfzc][i+nfzc] = CSC_full[str_j->occ[j]][str_i->occ[i]];
 
129
 
 
130
          // frozen orbitals need to be mapped to pitzer order manually
 
131
          for(int i=0; i<nfzc; i++) {
 
132
            int ii = fzc_occ[i];
 
133
            for(int j=0; j<nact_b; j++)
 
134
              CSC_b[j+nfzc][i] = CSC_full[str_j->occ[j]][ii];
 
135
          }
 
136
 
 
137
          for(int j=0; j<nfzc; j++) {
 
138
            int jj = fzc_occ[j];
 
139
            for(int i=0; i<nact_b; i++)
 
140
              CSC_b[j][i+nfzc] = CSC_full[jj][str_i->occ[i]];
 
141
          }
 
142
 
 
143
          for(int i=0;i<nfzc;i++) {
 
144
            int ii = fzc_occ[i];
 
145
            for(int j=0;j<nfzc;j++)
 
146
              CSC_b[i][j] = CSC_full[ii][fzc_occ[j]];
 
147
          }
 
148
 
 
149
          FLOAT sign;
 
150
          lu_decom(CSC_b, nbeta, tmpintvec, &sign);
 
151
          FLOAT deter1 = 1.0;
 
152
          for(int i=0;i<nbeta;i++)
 
153
          deter1 *= CSC_b[i][i];
 
154
 
 
155
          S_b[jp][im] = sign*deter1;
 
156
        }
 
157
      }
 
158
 
 
159
      // Evaluate total overlap in the highest available precision
 
160
      int ndets = vecm->size;
 
161
      FLOAT S_tot = 0.0;
 
162
      for(int I=0; I<ndets; I++) {
 
163
        SlaterDet *detI = vecm->sdset->dets + I;
 
164
        int Istra = detI->alphastring;
 
165
        int Istrb = detI->betastring;
 
166
        FLOAT cI = vecm->coeffs[I];
 
167
 
 
168
        for(int J=0; J<ndets; J++) {
 
169
          SlaterDet *detJ = vecp->sdset->dets + J;
 
170
          int Jstra = detJ->alphastring;
 
171
          int Jstrb = detJ->betastring;
 
172
          FLOAT cJ = vecp->coeffs[J];
 
173
 
 
174
          FLOAT S = S_a[Jstra][Istra] * S_b[Jstrb][Istrb];
 
175
 
 
176
          FLOAT contrib = cI * S * cJ;
 
177
          S_tot += cI * S * cJ;
 
178
          /* fprintf(outfile,"  %3d %3d %+15.10Le", I, J, cI);
 
179
           fprintf(outfile," %+15.10Le", cJ);
 
180
           fprintf(outfile," %+25.15Le", S);
 
181
           fprintf(outfile," %+25.15Le", contrib);
 
182
           fprintf(outfile," %+25.15Le\n", S_tot); */
 
183
 
 
184
        }
 
185
      }
 
186
 
 
187
      // Cleanup
 
188
      slaterdetvector_delete_full(vecm);
 
189
      slaterdetvector_delete_full(vecp);
 
190
      delete[] tmpintvec;
 
191
      delete_matrix(CSC_a);
 
192
      delete_matrix(CSC_full);
 
193
      delete_matrix(CSC_b);
 
194
      delete_matrix(S_a,nstr_a,nstr_a);
 
195
      delete_matrix(S_b,nstr_b,nstr_b);
 
196
      double S_tot_double = (double) S_tot;
 
197
      return fabs(S_tot_double);
 
198
    }
 
199
 
 
200
  }} // namespace psi::dboc