← Back to branch summary

~ubuntu-branches/ubuntu/precise/psicode/precise

« back to all changes in this revision

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

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