← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/psicode/intrepid

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
Import upstream version 3.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/dboc/rccsd.cc
 
1
#include <iostream>
 
2
#include <stdio.h>
 
3
#include <stdlib.h>
 
4
#include <math.h>
 
5
extern "C" {
 
6
#include <libciomr/libciomr.h>
 
7
#include <libqt/qt.h>
 
8
#include <psifiles.h>
 
9
}
 
10
#include "moinfo.h"
 
11
#include "mo_overlap.h"
 
12
#include "float.h"
 
13
#include "linalg.h"
 
14
 
 
15
using namespace std;
 
16
 
 
17
// Wrap a,b indices into one composite index assuming S2 symmetry
 
18
#define INDEX2(a,b) ((a) > (b)) ? ( (((a)*(a+1)) >> 1) + (b) ) : ( (((b)*(b+1)) >> 1) + (a) )
 
19
// Wrap a>=b indices into one composite index assuming S2 symmetry
 
20
#define INDEX2_ORD(a,b) ( (((a)*(a+1))/2) + (b) )
 
21
// Wrap a>=b>=c indices into one composite index assuming S3 symmetry
 
22
#define INDEX3_ORD(a,b,c) ( ((a)*(((a)+4)*((a)-1)+6)/6) + (((b)*(b+1))/2) + (c) )
 
23
 
 
24
extern MOInfo_t MOInfo;
 
25
extern FILE *outfile;
 
26
extern void done(const char *);
 
27
 
 
28
double eval_rccsd_derwfn_overlap()
 
29
{
 
30
  int ndocc = MOInfo.ndocc;
 
31
  int nvirt = MOInfo.num_mo - MOInfo.ndocc;
 
32
  FLOAT **CSC_full = eval_S_alpha();
 
33
  FLOAT **CSC = create_matrix(ndocc,ndocc);
 
34
  int *tmpintvec = new int[ndocc];
 
35
 
 
36
  FLOAT **RmSp = create_matrix(ndocc,nvirt);
 
37
  FLOAT **RmDp = create_matrix(INDEX2_ORD(ndocc,0),INDEX2_ORD(nvirt,0));
 
38
  FLOAT **RmTp = create_matrix(INDEX3_ORD(ndocc,0,0),INDEX3_ORD(nvirt,0,0));
 
39
 
 
40
  //
 
41
  // Evaluate reference-reference overlap <Ref(-)|Ref(+)>
 
42
  //
 
43
  for(int i=0;i<ndocc;i++)
 
44
    for(int j=0;j<ndocc;j++)
 
45
      CSC[i][j] = CSC_full[i][j];
 
46
  //  C_DGETRF(ndocc,ndocc,&(CSC[0][0]),ndocc,tmpintvec);
 
47
  FLOAT sign;
 
48
  lu_decom(CSC, ndocc, tmpintvec, &sign);
 
49
  FLOAT deter_ref = 1.0;
 
50
  for(int i=0;i<ndocc;i++)
 
51
    deter_ref *= CSC[i][i];
 
52
 
 
53
  //
 
54
  // Evaluate all overlaps <Ref(-)|Singles(+)>
 
55
  //
 
56
  fprintf(outfile,"\n  -Overlap of Ref(-) with Singles(-):\n");
 
57
  for(int mo_i=0; mo_i<ndocc; mo_i++) {
 
58
    for(int mo_a=ndocc; mo_a<MOInfo.num_mo; mo_a++) {
 
59
 
 
60
      // Before generating the single restore the original (<Ref(-)|Ref(+)>) overlap matrix
 
61
      for(int i=0;i<ndocc;i++)
 
62
        for(int j=0;j<ndocc;j++)
 
63
          CSC[i][j] = CSC_full[i][j];
 
64
      
 
65
      // Replace ith column with ath column
 
66
      for(int p=0; p<ndocc; p++)
 
67
        CSC[p][mo_i] = CSC_full[p][mo_a];
 
68
 
 
69
      // Compute the determinant
 
70
      //      C_DGETRF(ndocc,ndocc,&(CSC[0][0]),ndocc,tmpintvec);
 
71
      lu_decom(CSC, ndocc, tmpintvec, &sign);
 
72
      FLOAT deter1 = 1.0;
 
73
      for(int i=0;i<ndocc;i++)
 
74
        deter1 *= CSC[i][i];
 
75
 
 
76
      RmSp[mo_i][mo_a-ndocc] = deter_ref * deter1;
 
77
      fprintf(outfile,"  %d %d %20.15lf\n",mo_i,mo_a,RmSp[mo_i][mo_a-ndocc]);
 
78
 
 
79
    }
 
80
  }
 
81
 
 
82
  //
 
83
  // Evaluate all overlaps <Ref(-)|IJABDoubles(+)>
 
84
  //
 
85
  fprintf(outfile,"\n  -Overlap of Ref(-) with IJABDoubles(-):\n");
 
86
  for(int mo_i=0; mo_i<ndocc; mo_i++) {
 
87
    for(int mo_j=0; mo_j<mo_i; mo_j++) {
 
88
      for(int mo_a=ndocc, a=0; mo_a<MOInfo.num_mo; mo_a++,a++) {
 
89
        for(int mo_b=ndocc, b=0; mo_b<mo_a; mo_b++,b++) {
 
90
 
 
91
          // Before generating the single restore the original (<Ref(-)|Ref(+)>) overlap matrix
 
92
          for(int i=0;i<ndocc;i++)
 
93
            for(int j=0;j<ndocc;j++)
 
94
              CSC[i][j] = CSC_full[i][j];
 
95
          
 
96
          // Replace ith column with ath column
 
97
          for(int p=0; p<ndocc; p++)
 
98
            CSC[p][mo_i] = CSC_full[p][mo_a];
 
99
          // Replace jth column with bth column
 
100
          for(int p=0; p<ndocc; p++)
 
101
            CSC[p][mo_j] = CSC_full[p][mo_b];
 
102
          
 
103
          // Compute the determinant
 
104
          //      C_DGETRF(ndocc,ndocc,&(CSC[0][0]),ndocc,tmpintvec);
 
105
          lu_decom(CSC, ndocc, tmpintvec, &sign);
 
106
          FLOAT deter1 = 1.0;
 
107
          for(int i=0;i<ndocc;i++)
 
108
            deter1 *= CSC[i][i];
 
109
          
 
110
          int ij = INDEX2_ORD(mo_i,mo_j);
 
111
          int ab = INDEX2_ORD(a,b);
 
112
          RmDp[ij][ab] = deter_ref * deter1;
 
113
          fprintf(outfile,"  %d %d %d %d %20.15lf\n",mo_i,mo_j,mo_a,mo_b,RmDp[ij][ab]);
 
114
          
 
115
        }
 
116
      }
 
117
    }
 
118
  }
 
119
 
 
120
  //
 
121
  // Evaluate all overlaps <Ref(-)|IJKABCTriples(+)>
 
122
  //
 
123
  fprintf(outfile,"\n  -Overlap of Ref(-) with IJKABCTriples(-):\n");
 
124
  for(int mo_i=0; mo_i<ndocc; mo_i++) {
 
125
    for(int mo_j=0; mo_j<mo_i; mo_j++) {
 
126
      for(int mo_k=0; mo_k<mo_j; mo_k++) {
 
127
        for(int mo_a=ndocc, a=0; mo_a<MOInfo.num_mo; mo_a++,a++) {
 
128
          for(int mo_b=ndocc, b=0; mo_b<mo_a; mo_b++,b++) {
 
129
            for(int mo_c=ndocc, c=0; mo_c<mo_b; mo_c++,c++) {
 
130
 
 
131
              // Before generating the single restore the original (<Ref(-)|Ref(+)>) overlap matrix
 
132
              for(int i=0;i<ndocc;i++)
 
133
                for(int j=0;j<ndocc;j++)
 
134
                  CSC[i][j] = CSC_full[i][j];
 
135
          
 
136
              // Replace ith column with ath column
 
137
              for(int p=0; p<ndocc; p++)
 
138
                CSC[p][mo_i] = CSC_full[p][mo_a];
 
139
              // Replace jth column with bth column
 
140
              for(int p=0; p<ndocc; p++)
 
141
                CSC[p][mo_j] = CSC_full[p][mo_b];
 
142
              // Replace jth column with bth column
 
143
              for(int p=0; p<ndocc; p++)
 
144
                CSC[p][mo_k] = CSC_full[p][mo_c];
 
145
          
 
146
              // Compute the determinant
 
147
              // C_DGETRF(ndocc,ndocc,&(CSC[0][0]),ndocc,tmpintvec);
 
148
              lu_decom(CSC, ndocc, tmpintvec, &sign);
 
149
              FLOAT deter1 = 1.0;
 
150
              for(int i=0;i<ndocc;i++)
 
151
                deter1 *= CSC[i][i];
 
152
              
 
153
              int ijk = INDEX3_ORD(mo_i,mo_j,mo_k);
 
154
              int abc = INDEX3_ORD(a,b,c);
 
155
              RmTp[ijk][abc] = deter_ref * deter1;
 
156
              fprintf(outfile,"  %d %d %d %d %d %d %20.15lf\n",mo_i,mo_j,mo_k,mo_a,mo_b,mo_c,RmTp[ijk][abc]);
 
157
              
 
158
            }
 
159
          }
 
160
        }
 
161
      }
 
162
    }
 
163
  }
 
164
 
 
165
  delete[] tmpintvec;
 
166
  delete_matrix(CSC);
 
167
  delete_matrix(CSC_full);
 
168
  delete_matrix(RmTp);
 
169
  delete_matrix(RmDp);
 
170
  delete_matrix(RmSp);
 
171
  return (double)deter_ref*deter_ref;
 
172
}
 
173