← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cis/mp2.c

  • 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/cis/mp2.c
 
1
#include <math.h>
 
2
#include <libdpd/dpd.h>
 
3
#include <psifiles.h>
 
4
#define EXTERN
 
5
#include "globals.h"
 
6
 
 
7
void mp2(void)
 
8
{
 
9
  int iter, h, nirreps, row, col;
 
10
  double energy, conv, rms, value;
 
11
  dpdfile2 F;
 
12
  dpdbuf4 D, T2, newT2, Z;
 
13
 
 
14
  nirreps = moinfo.nirreps;
 
15
 
 
16
  if(params.ref == 0) { /** RHF **/
 
17
 
 
18
    /* build initial guess amplitudes */
 
19
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
20
    dpd_buf4_copy(&D, CC_MISC, "MP2 tIjAb");
 
21
    dpd_buf4_close(&D);
 
22
 
 
23
    dpd_buf4_init(&T2, CC_MISC, 0, 0, 5, 0, 5, 0, "MP2 tIjAb");
 
24
    if(params.local) local_filter_T2(&T2);
 
25
    else {
 
26
      dpd_buf4_init(&D, CC_DENOM, 0, 0, 5, 0, 5, 0, "dIjAb");
 
27
      dpd_buf4_dirprd(&D, &T2);
 
28
      dpd_buf4_close(&D);
 
29
    }
 
30
 
 
31
    dpd_buf4_copy(&T2, CC_MISC, "New MP2 tIjAb");
 
32
    dpd_buf4_close(&T2);
 
33
 
 
34
    dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D 2<ij|ab> - <ij|ba>");
 
35
    dpd_buf4_init(&T2, CC_MISC, 0, 0, 5, 0, 5, 0, "MP2 tIjAb");
 
36
    energy = dpd_buf4_dot(&D, &T2);
 
37
    dpd_buf4_close(&T2);
 
38
    dpd_buf4_close(&D);
 
39
 
 
40
    if(params.local) {
 
41
      fprintf(outfile, "\n\tSolving for LMP2 wave function:\n");
 
42
      fprintf(outfile,   "\t-------------------------------\n");
 
43
      fprintf(outfile, "\titer = %d  LMP2 Energy = %20.14f\n", 0, energy);
 
44
    }
 
45
    else {
 
46
      fprintf(outfile, "\n\tSolving for MP2 wave function:\n");
 
47
      fprintf(outfile,   "\t-------------------------------\n");
 
48
      fprintf(outfile, "\titer = %d  MP2 Energy = %20.14f\n", 0, energy);
 
49
    }
 
50
 
 
51
    conv = 0;
 
52
    for(iter=1; iter < params.maxiter; iter++) {
 
53
 
 
54
      dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
55
      dpd_buf4_copy(&D, CC_MISC, "New MP2 tIjAb Increment");
 
56
      dpd_buf4_close(&D);
 
57
 
 
58
      dpd_buf4_init(&newT2, CC_MISC, 0, 0, 5, 0, 5, 0, "New MP2 tIjAb Increment");
 
59
      dpd_buf4_init(&T2, CC_MISC, 0, 0, 5, 0, 5, 0, "MP2 tIjAb");
 
60
 
 
61
      dpd_file2_init(&F, CC_OEI, 0, 0, 0, "fIJ");
 
62
      dpd_contract424(&T2, &F, &newT2, 1, 0, 1, -1, 1);
 
63
      dpd_contract244(&F, &T2, &newT2, 0, 0, 0, -1, 1);
 
64
      dpd_file2_close(&F);
 
65
 
 
66
      dpd_file2_init(&F, CC_OEI, 0, 1, 1, "fAB");
 
67
      dpd_contract244(&F, &T2, &newT2, 1, 2, 1, 1, 1);
 
68
      dpd_contract424(&T2, &F, &newT2, 3, 1, 0, 1, 1);
 
69
      dpd_file2_close(&F);
 
70
 
 
71
      dpd_buf4_close(&T2);
 
72
 
 
73
      if(params.local) {
 
74
        local_filter_T2(&newT2);
 
75
      }
 
76
      else {
 
77
        dpd_buf4_init(&D, CC_DENOM, 0, 0, 5, 0, 5, 0, "dIjAb");
 
78
        dpd_buf4_dirprd(&D, &newT2);
 
79
        dpd_buf4_close(&D);
 
80
      }
 
81
 
 
82
      dpd_buf4_close(&newT2);
 
83
 
 
84
      dpd_buf4_init(&newT2, CC_MISC, 0, 0, 5, 0, 5, 0, "New MP2 tIjAb");
 
85
      dpd_buf4_init(&T2, CC_MISC, 0, 0, 5, 0, 5, 0, "New MP2 tIjAb Increment");
 
86
      dpd_buf4_axpy(&T2, &newT2, 1);
 
87
      dpd_buf4_close(&T2);
 
88
 
 
89
      dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D 2<ij|ab> - <ij|ba>");
 
90
      energy = dpd_buf4_dot(&D, &newT2);
 
91
      dpd_buf4_close(&D);
 
92
      dpd_buf4_close(&newT2);
 
93
 
 
94
      dpd_buf4_init(&newT2, CC_MISC, 0, 0, 5, 0, 5, 0, "New MP2 tIjAb");
 
95
      dpd_buf4_init(&T2, CC_MISC, 0, 0, 5, 0, 5, 0, "MP2 tIjAb");
 
96
      rms = 0.0;
 
97
      for(h=0; h < nirreps; h++) {
 
98
        dpd_buf4_mat_irrep_init(&newT2, h);
 
99
        dpd_buf4_mat_irrep_rd(&newT2, h);
 
100
        dpd_buf4_mat_irrep_init(&T2, h);
 
101
        dpd_buf4_mat_irrep_rd(&T2, h);
 
102
 
 
103
        for(row=0; row < T2.params->rowtot[h]; row++)
 
104
          for(col=0; col < T2.params->coltot[h]; col++) {
 
105
            value = newT2.matrix[h][row][col] - T2.matrix[h][row][col];
 
106
            rms += value * value;
 
107
          }
 
108
 
 
109
        dpd_buf4_mat_irrep_close(&T2, h);
 
110
        dpd_buf4_mat_irrep_close(&newT2, h);
 
111
      }
 
112
      dpd_buf4_close(&T2);
 
113
      dpd_buf4_close(&newT2);
 
114
      rms = sqrt(rms);
 
115
 
 
116
      if(params.local) {
 
117
        fprintf(outfile, "\titer = %d   LMP2 Energy = %20.14f   RMS = %4.3e\n", iter, energy, rms);
 
118
      }
 
119
      else {
 
120
        fprintf(outfile, "\titer = %d   MP2 Energy = %20.14f   RMS = %4.3e\n", iter, energy, rms);
 
121
      }
 
122
 
 
123
      if(rms < params.convergence) {
 
124
        conv = 1;
 
125
        fprintf(outfile, "\n\tMP2 iterations converged.\n\n");
 
126
        break;
 
127
      }
 
128
      else {
 
129
        dpd_buf4_init(&T2, CC_MISC, 0, 0, 5, 0, 5, 0, "New MP2 tIjAb");
 
130
        dpd_buf4_copy(&T2, CC_MISC, "MP2 tIjAb");
 
131
        dpd_buf4_close(&T2);
 
132
      }
 
133
    }
 
134
 
 
135
    if(!conv) {
 
136
      fprintf(outfile, "\n\tMP2 iterative procedure failed.\n");
 
137
      exit(PSI_RETURN_FAILURE);
 
138
    }
 
139
 
 
140
    /* spin adapt the final amplitudes */
 
141
    dpd_buf4_init(&T2, CC_MISC, 0, 0, 5, 0, 5, 0, "MP2 tIjAb");
 
142
    dpd_buf4_sort(&T2, CC_TMP0, pqsr, 0, 5, "MP2 tIjbA");
 
143
    dpd_buf4_copy(&T2, CC_MISC, "MP2 2 tIjAb - tIjbA");
 
144
    dpd_buf4_close(&T2);
 
145
    dpd_buf4_init(&T2, CC_MISC, 0, 0, 5, 0, 5, 0, "MP2 2 tIjAb - tIjbA");
 
146
    dpd_buf4_scm(&T2, 2);
 
147
    dpd_buf4_init(&Z, CC_TMP0, 0, 0, 5, 0, 5, 0, "MP2 tIjbA");
 
148
    dpd_buf4_axpy(&Z, &T2, -1);
 
149
    dpd_buf4_close(&Z);
 
150
    dpd_buf4_close(&T2);
 
151
  }
 
152
  else if(params.ref == 2) { /** UHF **/
 
153
    dpd_buf4_init(&D, CC_DINTS, 0, 2, 7, 2, 7, 0, "D <IJ||AB> (I>J,A>B)");
 
154
    dpd_buf4_copy(&D, CC_MISC, "MP2 tIJAB");
 
155
    dpd_buf4_close(&D);
 
156
    dpd_buf4_init(&T2, CC_MISC, 0, 2, 7, 2, 7, 0, "MP2 tIJAB");
 
157
    dpd_buf4_init(&D, CC_DENOM, 0, 1, 6, 1, 6, 0, "dIJAB");
 
158
    dpd_buf4_dirprd(&D, &T2);
 
159
    dpd_buf4_close(&D);
 
160
    dpd_buf4_close(&T2);
 
161
 
 
162
    dpd_buf4_init(&D, CC_DINTS, 0, 12, 17, 12, 17, 0, "D <ij||ab> (i>j,a>b)");
 
163
    dpd_buf4_copy(&D, CC_MISC, "MP2 tijab");
 
164
    dpd_buf4_close(&D);
 
165
    dpd_buf4_init(&T2, CC_MISC, 0, 12, 17, 12, 17, 0, "MP2 tijab");
 
166
    dpd_buf4_init(&D, CC_DENOM, 0, 11, 16, 11, 16, 0, "dijab");
 
167
    dpd_buf4_dirprd(&D, &T2);
 
168
    dpd_buf4_close(&D);
 
169
    dpd_buf4_close(&T2);
 
170
 
 
171
    dpd_buf4_init(&D, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
 
172
    dpd_buf4_copy(&D, CC_MISC, "MP2 tIjAb");
 
173
    dpd_buf4_close(&D);
 
174
    dpd_buf4_init(&T2, CC_MISC, 0, 22, 28, 22, 28, 0, "MP2 tIjAb");
 
175
    dpd_buf4_init(&D, CC_DENOM, 0, 22, 28, 22, 28, 0, "dIjAb");
 
176
    dpd_buf4_dirprd(&D, &T2);
 
177
    dpd_buf4_close(&D);
 
178
    dpd_buf4_close(&T2);
 
179
  }
 
180
}
 
181