← Back to branch summary

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

« back to all changes in this revision

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