← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/response/build_A_ROHF.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/response/build_A_ROHF.cc
 
1
/*! \file
 
2
    \ingroup RESPONSE
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <libdpd/dpd.h>
 
6
#include <psifiles.h>
 
7
#include "MOInfo.h"
 
8
#include "Params.h"
 
9
#define EXTERN
 
10
#include "globals.h"
 
11
 
 
12
namespace psi { namespace response {
 
13
 
 
14
/* BUILD_A_ROHF(): Construct the molecular orbital Hessian, A, for
 
15
** ROHF orbitals. At the moment we're actually building all symmetry
 
16
** blocks of A, though for the orbital Z-vector equations we really
 
17
** only need the totally symmetric components.
 
18
**
 
19
** A(EM,AI) = 2<MI|EA> - <IM|EA> - <ME|IA> + del_MI fEA - del_EA fMI
 
20
**
 
21
** A(em,ai) = 2<mi|ea> - <im|ea> - <me|ia> + del_mi fea - del_ea fmi
 
22
**
 
23
** A(EM,ai) = 2<Mi|Ea> + del_Ma fei
 
24
**
 
25
** */
 
26
 
 
27
void build_A_ROHF(void)
 
28
{
 
29
  int h, nirreps, e, m, a, i, em, ai, E, M, A, I;
 
30
  int Esym, Msym, Asym, Isym;
 
31
  int *virtpi, *openpi, *occpi, *occ_off, *vir_off;
 
32
  int *qt_occ, *qt_vir; /* Spatial orbital translators */
 
33
  dpdfile2 fIJ, fij, fAB, fab, fIA, fia;
 
34
  dpdbuf4 Amat, Amat2, D, C;
 
35
 
 
36
  nirreps = moinfo.nirreps;
 
37
  occpi = moinfo.occpi; openpi = moinfo.openpi; virtpi = moinfo.virtpi;
 
38
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
 
39
  qt_occ = moinfo.qt_occ; qt_vir = moinfo.qt_vir;
 
40
 
 
41
  /* Two-electron integral contributions */
 
42
  dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 0, "D <ij|ab>");
 
43
  dpd_buf4_sort(&D, PSIF_MO_HESS, rpsq, 11, 11, "A(EM,AI)");
 
44
  dpd_buf4_close(&D);
 
45
  dpd_buf4_init(&Amat, PSIF_MO_HESS, 0, 11, 11, 11, 11, 0, "A(EM,AI)");
 
46
  dpd_buf4_sort(&Amat, CC_TMP0, psrq, 11, 11, "D <im|ea> (ei,am)");
 
47
  dpd_buf4_scm(&Amat, 2.0);
 
48
  dpd_buf4_copy(&Amat, CC_TMP0, "A(EM,ai)");
 
49
  dpd_buf4_init(&D, CC_TMP0, 0, 11, 11, 11, 11, 0, "D <im|ea> (ei,am)");
 
50
  dpd_buf4_axpy(&D, &Amat, -1.0);
 
51
  dpd_buf4_close(&D);
 
52
  dpd_buf4_init(&C, CC_CINTS, 0, 10, 10, 10, 10, 0, "C <ia|jb>");
 
53
  dpd_buf4_sort(&C, CC_TMP0, qpsr, 11, 11, "C <ai|bj>");
 
54
  dpd_buf4_close(&C);
 
55
  dpd_buf4_init(&C, CC_TMP0, 0, 11, 11, 11, 11, 0, "C <ai|bj>");
 
56
  dpd_buf4_axpy(&C, &Amat, -1.0);
 
57
  dpd_buf4_close(&C);
 
58
  dpd_buf4_copy(&Amat, CC_TMP0, "A(em,ai)");
 
59
  dpd_buf4_close(&Amat);
 
60
 
 
61
  /* Fock matrix contributions */
 
62
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
 
63
  dpd_file2_mat_init(&fIJ);
 
64
  dpd_file2_mat_rd(&fIJ);
 
65
  dpd_file2_init(&fij, CC_OEI, 0, 0, 0, "fij");
 
66
  dpd_file2_mat_init(&fij);
 
67
  dpd_file2_mat_rd(&fij);
 
68
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
 
69
  dpd_file2_mat_init(&fAB);
 
70
  dpd_file2_mat_rd(&fAB);
 
71
  dpd_file2_init(&fab, CC_OEI, 0, 1, 1, "fab");
 
72
  dpd_file2_mat_init(&fab);
 
73
  dpd_file2_mat_rd(&fab);
 
74
  dpd_file2_init(&fIA, CC_OEI, 0, 0, 1, "fIA");
 
75
  dpd_file2_mat_init(&fIA);
 
76
  dpd_file2_mat_rd(&fIA);
 
77
  dpd_file2_init(&fia, CC_OEI, 0, 0, 1, "fia");
 
78
  dpd_file2_mat_init(&fia);
 
79
  dpd_file2_mat_rd(&fia);
 
80
 
 
81
  dpd_buf4_init(&Amat, PSIF_MO_HESS, 0, 11, 11, 11, 11, 0, "A(EM,AI)");
 
82
  
 
83
  for(h=0; h < nirreps; h++) {
 
84
    dpd_buf4_mat_irrep_init(&Amat, h);
 
85
    dpd_buf4_mat_irrep_rd(&Amat, h);
 
86
 
 
87
    for(em=0; em < Amat.params->rowtot[h]; em++) {
 
88
      e = Amat.params->roworb[h][em][0];
 
89
      m = Amat.params->roworb[h][em][1];
 
90
      E = fAB.params->rowidx[e];
 
91
      M = fIJ.params->rowidx[m];
 
92
      Esym = fAB.params->psym[e];
 
93
      Msym = fIJ.params->psym[m];
 
94
      for(ai=0; ai < Amat.params->coltot[h]; ai++) {
 
95
        a = Amat.params->colorb[h][ai][0];
 
96
        i = Amat.params->colorb[h][ai][1];
 
97
        A = fAB.params->colidx[a];
 
98
        I = fIJ.params->colidx[i];
 
99
        Asym = fAB.params->qsym[a];
 
100
        Isym = fIJ.params->qsym[i];
 
101
 
 
102
        if((M==I) && (Esym==Asym))
 
103
          Amat.matrix[h][em][ai] += fAB.matrix[Esym][E][A];
 
104
        if((E==A) && (Msym==Isym))
 
105
          Amat.matrix[h][em][ai] -= fIJ.matrix[Msym][M][I];
 
106
 
 
107
        /* Check to see if these virtual indices actually
 
108
           correspond to open-shell orbitals --- if so, set this
 
109
           element to zero */
 
110
        if((E >= (virtpi[Esym] - openpi[Esym])) ||
 
111
           (A >= (virtpi[Asym] - openpi[Asym])) )
 
112
          Amat.matrix[h][em][ai] = 0.0;
 
113
      }
 
114
    }
 
115
 
 
116
    dpd_buf4_mat_irrep_wrt(&Amat, h);
 
117
    dpd_buf4_mat_irrep_close(&Amat, h);
 
118
  }
 
119
 
 
120
  dpd_buf4_close(&Amat);
 
121
 
 
122
  dpd_buf4_init(&Amat, CC_TMP0, 0, 11, 11, 11, 11, 0, "A(em,ai)");
 
123
  
 
124
  for(h=0; h < nirreps; h++) {
 
125
    dpd_buf4_mat_irrep_init(&Amat, h);
 
126
    dpd_buf4_mat_irrep_rd(&Amat, h);
 
127
 
 
128
    for(em=0; em < Amat.params->rowtot[h]; em++) {
 
129
      e = Amat.params->roworb[h][em][0];
 
130
      m = Amat.params->roworb[h][em][1];
 
131
      E = fab.params->rowidx[e];
 
132
      M = fij.params->rowidx[m];
 
133
      Esym = fab.params->psym[e];
 
134
      Msym = fij.params->psym[m];
 
135
      for(ai=0; ai < Amat.params->coltot[h]; ai++) {
 
136
        a = Amat.params->colorb[h][ai][0];
 
137
        i = Amat.params->colorb[h][ai][1];
 
138
        A = fab.params->colidx[a];
 
139
        I = fij.params->colidx[i];
 
140
        Asym = fab.params->qsym[a];
 
141
        Isym = fij.params->qsym[i];
 
142
 
 
143
        if((M==I) && (Esym==Asym))
 
144
          Amat.matrix[h][em][ai] += fab.matrix[Esym][E][A];
 
145
        if((E==A) && (Msym==Isym))
 
146
          Amat.matrix[h][em][ai] -= fij.matrix[Msym][M][I];
 
147
 
 
148
        /* Check to see if these occupied indices actually
 
149
           correspond to open-shell orbitals --- if so, set this
 
150
           element to zero */
 
151
        if((M >= (occpi[Msym] - openpi[Msym])) ||
 
152
           (I >= (occpi[Isym] - openpi[Isym])) )
 
153
          Amat.matrix[h][em][ai] = 0.0;
 
154
      }
 
155
    }
 
156
 
 
157
    dpd_buf4_mat_irrep_wrt(&Amat, h);
 
158
    dpd_buf4_mat_irrep_close(&Amat, h);
 
159
  }
 
160
 
 
161
  dpd_buf4_close(&Amat);
 
162
 
 
163
  dpd_buf4_init(&Amat, CC_TMP0, 0, 11, 11, 11, 11, 0, "A(EM,ai)");
 
164
 
 
165
  for(h=0; h < nirreps; h++) {
 
166
    dpd_buf4_mat_irrep_init(&Amat, h);
 
167
    dpd_buf4_mat_irrep_rd(&Amat, h);
 
168
 
 
169
    for(em=0; em < Amat.params->rowtot[h]; em++) {
 
170
      e = Amat.params->roworb[h][em][0];
 
171
      m = Amat.params->roworb[h][em][1];
 
172
      Esym = Amat.params->psym[e];
 
173
      Msym = Amat.params->qsym[m];
 
174
      E = e - vir_off[Esym];
 
175
      M = m - occ_off[Msym];
 
176
      for(ai=0; ai < Amat.params->coltot[h]; ai++) {
 
177
        a = Amat.params->colorb[h][ai][0];
 
178
        i = Amat.params->colorb[h][ai][1];
 
179
        Asym = Amat.params->rsym[a];
 
180
        Isym = Amat.params->ssym[i];
 
181
        A = a - vir_off[Asym];
 
182
        I = i - occ_off[Isym];
 
183
 
 
184
        /* This comparison is somewhat tricky.  The algebraic
 
185
           expression for the Fock matrix shift here is:
 
186
 
 
187
           A(EM,ai) += delta(M,a) f(E,i)(beta)
 
188
 
 
189
           The Kronecker Delta is actually a comparison between
 
190
           the *spatial* orbitals associated with M, and A.
 
191
           Hence we have to compare the spatial orbital
 
192
           translation of the the two absolute orbital indices. */
 
193
        if((qt_occ[m] == qt_vir[a]) && (Esym==Isym))
 
194
          Amat.matrix[h][em][ai] += fia.matrix[Isym][I][E];
 
195
 
 
196
        /* Check to see if these occupied and virtual indices
 
197
           actually correspond to open-shell orbitals --- if so,
 
198
           set this element to zero */
 
199
        if((E >= (virtpi[Esym] - openpi[Esym])) ||
 
200
           (I >= (occpi[Isym] - openpi[Isym])) )
 
201
          Amat.matrix[h][em][ai] = 0.0;
 
202
      }
 
203
    }
 
204
 
 
205
    dpd_buf4_mat_irrep_wrt(&Amat, h);
 
206
    dpd_buf4_mat_irrep_close(&Amat, h);
 
207
  }
 
208
  dpd_buf4_sort(&Amat, CC_TMP0, rspq, 11, 11, "A(em,AI)");
 
209
  dpd_buf4_close(&Amat);
 
210
 
 
211
  dpd_file2_mat_close(&fIJ);
 
212
  dpd_file2_close(&fIJ);
 
213
  dpd_file2_mat_close(&fij);
 
214
  dpd_file2_close(&fij);
 
215
  dpd_file2_mat_close(&fAB);
 
216
  dpd_file2_close(&fAB);
 
217
  dpd_file2_mat_close(&fab);
 
218
  dpd_file2_close(&fab);
 
219
  dpd_file2_mat_close(&fIA);
 
220
  dpd_file2_close(&fIA);
 
221
  dpd_file2_mat_close(&fia);
 
222
  dpd_file2_close(&fia);
 
223
 
 
224
  /* Now sum all three A-matrix components and divide by 2 */
 
225
  dpd_buf4_init(&Amat, PSIF_MO_HESS, 0, 11, 11, 11, 11, 0, "A(EM,AI)");
 
226
  dpd_buf4_copy(&Amat, CC_TMP0, "A(EM,AI)");
 
227
  dpd_buf4_init(&Amat2, CC_TMP0, 0, 11, 11, 11, 11, 0, "A(em,ai)");
 
228
  dpd_buf4_axpy(&Amat2, &Amat, 1.0);
 
229
  dpd_buf4_close(&Amat2);
 
230
  dpd_buf4_init(&Amat2, CC_TMP0, 0, 11, 11, 11, 11, 0, "A(EM,ai)");
 
231
  dpd_buf4_axpy(&Amat2, &Amat, 1.0);
 
232
  dpd_buf4_close(&Amat2);
 
233
  dpd_buf4_init(&Amat2, CC_TMP0, 0, 11, 11, 11, 11, 0, "A(em,AI)");
 
234
  dpd_buf4_axpy(&Amat2, &Amat, 1.0);
 
235
  dpd_buf4_close(&Amat2);
 
236
  dpd_buf4_scm(&Amat, 0.5);
 
237
  dpd_buf4_close(&Amat);
 
238
}
 
239
 
 
240
 
 
241
}} // namespace psi::response