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

« back to all changes in this revision

Viewing changes to src/bin/ccdensity/sortI_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:
 
1
/*! \file
 
2
    \ingroup CCDENSITY
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <libdpd/dpd.h>
 
7
#include <cmath>
 
8
#include <libciomr/libciomr.h>
 
9
#include <libiwl/iwl.h>
 
10
#include "MOInfo.h"
 
11
#include "Params.h"
 
12
#include "Frozen.h"
 
13
#define EXTERN
 
14
#include "globals.h"
 
15
 
 
16
namespace psi { namespace ccdensity {
 
17
 
 
18
/* SORTI_ROHF(): Place all the components of the ROHF Lagrangian into
 
19
** a large matrix, I (moinfo.I), which we also symmetrize by computing
 
20
** Ipq = 1/2 (Ipq + Iqp).  This matrix is later written to disk in
 
21
** dump() for subsequent backtransformation.  Note that some of the
 
22
** components of the Lagrangian computed into the IIJ, Iij, IIA, and
 
23
** Iia matrices remain non-symmetric (e.g., IIJ neq IJI).  I re-used
 
24
** my sortone.c code here, so don't let some of the variable names
 
25
** confuse you. */
 
26
 
 
27
void sortI_ROHF(void)
 
28
{
 
29
  int h, nirreps, nmo, nfzv, nfzc, nclsd, nopen;
 
30
  int row, col, i, j, I, J, a, b, A, B, p, q;
 
31
  int *occpi, *virtpi, *occ_off, *vir_off; 
 
32
  int *occ_sym, *vir_sym, *openpi;
 
33
  int *qt_occ, *qt_vir;
 
34
  double **O, chksum, value;
 
35
  dpdfile2 D;
 
36
 
 
37
  nmo = moinfo.nmo;
 
38
  nfzc = moinfo.nfzc;
 
39
  nfzv = moinfo.nfzv;
 
40
  nclsd = moinfo.nclsd;
 
41
  nopen = moinfo.nopen;
 
42
  nirreps = moinfo.nirreps;
 
43
  occpi = moinfo.occpi; virtpi = moinfo.virtpi;
 
44
  occ_off = moinfo.occ_off; vir_off = moinfo.vir_off;
 
45
  occ_sym = moinfo.occ_sym; vir_sym = moinfo.vir_sym;
 
46
  openpi = moinfo.openpi;
 
47
  qt_occ = moinfo.qt_occ; qt_vir = moinfo.qt_vir;
 
48
 
 
49
  O = block_matrix(nmo,nmo);
 
50
 
 
51
  /* Sort alpha components first */
 
52
  dpd_file2_init(&D, CC_OEI, 0, 0, 0, "I(I,J)");
 
53
  dpd_file2_mat_init(&D);
 
54
  dpd_file2_mat_rd(&D);
 
55
  for(h=0; h < nirreps; h++) {
 
56
      for(i=0; i < occpi[h]; i++) {
 
57
          I = qt_occ[occ_off[h] + i];
 
58
          for(j=0; j < occpi[h]; j++) {
 
59
              J = qt_occ[occ_off[h] + j];
 
60
              O[I][J] += D.matrix[h][i][j];
 
61
            }
 
62
        }
 
63
    }
 
64
  dpd_file2_mat_close(&D);
 
65
  dpd_file2_close(&D);
 
66
 
 
67
  dpd_file2_init(&D, CC_OEI, 0, 1, 1, "I'AB");
 
68
  dpd_file2_mat_init(&D);
 
69
  dpd_file2_mat_rd(&D);
 
70
  for(h=0; h < nirreps; h++) {
 
71
      for(a=0; a < (virtpi[h] - openpi[h]); a++) {
 
72
          A = qt_vir[vir_off[h] + a];
 
73
          for(b=0; b < (virtpi[h] - openpi[h]); b++) {
 
74
              B = qt_vir[vir_off[h] + b];
 
75
 
 
76
              O[A][B] += D.matrix[h][a][b];
 
77
            }
 
78
        }
 
79
    }
 
80
  dpd_file2_mat_close(&D);
 
81
  dpd_file2_close(&D);
 
82
 
 
83
  dpd_file2_init(&D, CC_OEI, 0, 0, 1, "I(I,A)");
 
84
  dpd_file2_mat_init(&D);
 
85
  dpd_file2_mat_rd(&D);
 
86
  for(h=0; h < nirreps; h++) {
 
87
      for(i=0; i < occpi[h]; i++) {
 
88
          I = qt_occ[occ_off[h] + i];
 
89
          for(a=0; a < (virtpi[h] - openpi[h]); a++) {
 
90
              A = qt_vir[vir_off[h] + a];
 
91
 
 
92
              O[A][I] += D.matrix[h][i][a];
 
93
              O[I][A] += D.matrix[h][i][a];
 
94
            }
 
95
        }
 
96
    }
 
97
  dpd_file2_mat_close(&D);
 
98
  dpd_file2_close(&D);
 
99
 
 
100
  /* Sort beta components */
 
101
  dpd_file2_init(&D, CC_OEI, 0, 0, 0, "I(i,j)");
 
102
  dpd_file2_mat_init(&D); 
 
103
  dpd_file2_mat_rd(&D);
 
104
  for(h=0; h < nirreps; h++) {
 
105
      for(i=0; i < (occpi[h] - openpi[h]); i++) { 
 
106
          I = qt_occ[occ_off[h] + i];
 
107
          for(j=0; j < (occpi[h] - openpi[h]); j++) {
 
108
              J = qt_occ[occ_off[h] + j];
 
109
              O[I][J] += D.matrix[h][i][j];
 
110
            }
 
111
        }
 
112
    }
 
113
  dpd_file2_mat_close(&D);
 
114
  dpd_file2_close(&D);
 
115
 
 
116
  dpd_file2_init(&D, CC_OEI, 0, 1, 1, "I'ab");
 
117
  dpd_file2_mat_init(&D);
 
118
  dpd_file2_mat_rd(&D);
 
119
  for(h=0; h < nirreps; h++) {
 
120
      for(a=0; a < virtpi[h]; a++) {
 
121
          A = qt_vir[vir_off[h] + a];
 
122
          for(b=0; b < virtpi[h]; b++) {
 
123
              B = qt_vir[vir_off[h] + b];
 
124
 
 
125
              O[A][B] += D.matrix[h][a][b];
 
126
            }
 
127
        }
 
128
    }
 
129
  dpd_file2_mat_close(&D);
 
130
  dpd_file2_close(&D);
 
131
 
 
132
  dpd_file2_init(&D, CC_OEI, 0, 0, 1, "I(i,a)");
 
133
  dpd_file2_mat_init(&D);
 
134
  dpd_file2_mat_rd(&D);
 
135
  for(h=0; h < nirreps; h++) {
 
136
      for(i=0; i < (occpi[h] - openpi[h]); i++) {
 
137
          I = qt_occ[occ_off[h] + i];
 
138
          for(a=0; a < virtpi[h]; a++) {
 
139
              A = qt_vir[vir_off[h] + a];
 
140
 
 
141
              O[A][I] += D.matrix[h][i][a];
 
142
              O[I][A] += D.matrix[h][i][a];
 
143
            }
 
144
        }
 
145
    }
 
146
  dpd_file2_mat_close(&D);
 
147
  dpd_file2_close(&D);
 
148
 
 
149
  /* Symmetrize the Lagrangian */
 
150
  for(p=0; p < (nmo-nfzv); p++) {
 
151
      for(q=0; q < p; q++) {
 
152
          value = 0.5*(O[p][q] + O[q][p]);
 
153
          O[p][q] = O[q][p] = value;
 
154
        }
 
155
    }
 
156
 
 
157
  /* Multiply the Lagrangian by -2.0 for the final energy derivative
 
158
     expression */
 
159
  for(p=0; p < (nmo-nfzv); p++) {
 
160
      for(q=0; q < (nmo-nfzv); q++) {
 
161
          O[p][q] *= -2.0;
 
162
        }
 
163
    }
 
164
 
 
165
  moinfo.I = O;
 
166
}
 
167
 
 
168
}} // namespace psi::ccdensity