← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cchbar/Wabei_RHF_FT2_a.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/cchbar/Wabei_RHF_FT2_a.c
1
 
#include <stdio.h>
2
 
#include <stdlib.h>
3
 
#include <libdpd/dpd.h>
4
 
#include <libqt/qt.h>
5
 
#define EXTERN
6
 
#include "globals.h"
7
 
 
8
 
/* Wabei_RHF_FT2_a(): Computes the following contribution to the Wabei HBAR
9
 
** matrix elements for RHF references:
10
 
**
11
 
** Z_aeib = <am|ef> [ 2 t_mi^fb - t_mi^bf ] - <mf|ae> t_mi^fb
12
 
**
13
 
** This is the final term of the Wabef HBAR element's contribution to Wabei.
14
 
** [cf. Gauss and Stanton, JCP 103, 3561-3577 (1995).]
15
 
** 
16
 
** TDC, March 2004
17
 
*/
18
 
 
19
 
void Wabei_RHF_FT2_a(void)
20
 
{
21
 
  int h, nirreps; 
22
 
  int Gam, Ga, Gm;
23
 
  int Gef, Gmf, Ge, Gf, Gae;
24
 
  int a, A, e, E, m, M, f, FF;
25
 
  int mf, am , ef;
26
 
  int *virtpi, *vir_off, *occpi, *occ_off;
27
 
  double ***W;
28
 
  double value;
29
 
  dpdbuf4 F, T2, Z;
30
 
 
31
 
  nirreps = moinfo.nirreps;
32
 
  virtpi = moinfo.virtpi;
33
 
  vir_off = moinfo.vir_off;
34
 
  occpi = moinfo.occpi;
35
 
  occ_off = moinfo.occ_off;
36
 
 
37
 
  /* Zaeib <-- - <mf|ae> t_mi^fb */
38
 
  /**** this contraction still requires storage of a full <ia|bc> set *****/
39
 
  dpd_buf4_init(&F, CC_FINTS, 0, 10, 5, 10, 5, 0, "F <ia|bc>");
40
 
  dpd_buf4_init(&T2, CC_TAMPS, 0, 10, 10, 10, 10, 0, "tIAjb");
41
 
  dpd_buf4_init(&Z, CC_TMP2, 0, 5, 10, 5, 10, 0, "Z(AE,ib)");
42
 
  dpd_contract444(&F, &T2, &Z, 1, 1, -1, 0);
43
 
  dpd_buf4_close(&T2);
44
 
  dpd_buf4_close(&F);
45
 
  dpd_buf4_close(&Z);
46
 
 
47
 
  /* Zaeib <-- <am|ef> [ 2 t_mi^fb - t_mi^bf ] */
48
 
 
49
 
   dpd_buf4_init(&Z, CC_TMP2, 0, 5, 10, 5, 10, 0, "Z(AE,ib)");
50
 
  dpd_buf4_init(&F, CC_FINTS, 0, 11, 5, 11, 5, 0, "F <ai|bc>");
51
 
  dpd_buf4_init(&T2, CC_TAMPS, 0, 10, 10, 10, 10, 0, "2 tIAjb - tIBja");
52
 
  for(h=0; h < nirreps; h++) {
53
 
    dpd_buf4_mat_irrep_init(&T2, h);
54
 
    dpd_buf4_mat_irrep_rd(&T2, h);
55
 
  }
56
 
 
57
 
  W = (double ***) malloc(nirreps * sizeof(double **));
58
 
 
59
 
  for(Ga=0; Ga < nirreps; Ga++) {
60
 
 
61
 
    /* allocate space for the input integrals */
62
 
    for(Gm=0; Gm < nirreps; Gm++) {
63
 
      Gam = Ga^Gm;
64
 
      F.matrix[Gam] = dpd_block_matrix(occpi[Gm], F.params->coltot[Gam]);
65
 
    }
66
 
 
67
 
    /* allocate space for the reordered integrals, W[e][mf] and target, Z[e][ib] */
68
 
    for(Ge=0; Ge < nirreps; Ge++) {
69
 
      Gae = Ga^Ge; Gmf = Gae;
70
 
      W[Ge] = dpd_block_matrix(virtpi[Ge], Z.params->coltot[Gmf]);
71
 
      Z.matrix[Gae] = dpd_block_matrix(virtpi[Ge], Z.params->coltot[Gmf]);
72
 
    }
73
 
 
74
 
    for(a=0; a < virtpi[Ga]; a++) {
75
 
      A = vir_off[Ga] + a;
76
 
 
77
 
      for(Gm=0; Gm < nirreps; Gm++) {
78
 
        Gam = Ga^Gm;
79
 
        Gef = Gam;
80
 
 
81
 
        /* read all <am|ef> integrals for given orbital index a --> F[m][ef] */
82
 
        dpd_buf4_mat_irrep_rd_block(&F, Gam, F.params->start13[Gam][A], occpi[Gm]);
83
 
      }
84
 
 
85
 
      for(Ge=0; Ge < nirreps; Ge++) {
86
 
        Gae = Ga^Ge; Gmf = Gae;
87
 
 
88
 
        /* sort F[m][ef] --> W[e][mf] */
89
 
        for(e=0; e < virtpi[Ge]; e++) {
90
 
          E = vir_off[Ge] + e;
91
 
 
92
 
          for(Gm=0; Gm < nirreps; Gm++) {
93
 
            Gf = Gmf^Gm;
94
 
            Gam = Ga^Gm;
95
 
 
96
 
            for(m=0; m < occpi[Gm]; m++) {
97
 
              M = occ_off[Gm] + m;
98
 
 
99
 
              for(f=0; f < virtpi[Gf]; f++) {
100
 
                FF = vir_off[Gf] + f;
101
 
 
102
 
                mf = Z.params->colidx[M][FF];
103
 
                am = F.params->rowidx[A][M];
104
 
                ef = F.params->colidx[E][FF];
105
 
 
106
 
                W[Ge][e][mf] = F.matrix[Gam][am-F.params->start13[Gam][A]][ef];
107
 
              } /* f */
108
 
            } /* m */
109
 
          } /* Gm */
110
 
        } /* e */
111
 
 
112
 
        /* read all existing Z(ae,ib) for a given orbital index a -->  Z[e][ib] */
113
 
        dpd_buf4_mat_irrep_rd_block(&Z, Gae, Z.params->start13[Gae][A], virtpi[Ge]);
114
 
 
115
 
        /* contract W[e][mf] * T[mf][ib] -> Z[e][ib] */
116
 
        if(virtpi[Ge] && Z.params->coltot[Gmf])
117
 
          C_DGEMM('n','n', virtpi[Ge], Z.params->coltot[Gmf], Z.params->coltot[Gmf], 1.0,
118
 
                  W[Ge][0], Z.params->coltot[Gmf], T2.matrix[Gmf][0], Z.params->coltot[Gmf],
119
 
                  1.0, Z.matrix[Gae][0], Z.params->coltot[Gmf]);
120
 
 
121
 
        /* write out the new Z(ae,ib) */
122
 
        dpd_buf4_mat_irrep_wrt_block(&Z, Gae, Z.params->start13[Gae][A], virtpi[Ge]);
123
 
 
124
 
      } /* Ge */
125
 
 
126
 
    } /* a */
127
 
 
128
 
    for(Gm=0; Gm < nirreps; Gm++) {
129
 
      Gam = Ga^Gm;
130
 
      dpd_free_block(F.matrix[Gam], occpi[Gm], F.params->coltot[Gam]);
131
 
    }
132
 
 
133
 
    for(Ge=0; Ge < nirreps; Ge++) {
134
 
      Gae = Ga^Ge; Gmf = Gae;
135
 
      dpd_free_block(W[Ge], virtpi[Ge], Z.params->coltot[Gmf]);
136
 
      dpd_free_block(Z.matrix[Gae], virtpi[Ge], Z.params->coltot[Gmf]);
137
 
    }
138
 
 
139
 
  } /* Ga */
140
 
 
141
 
  for(h=0; h < nirreps; h++) dpd_buf4_mat_irrep_close(&T2, h);
142
 
  dpd_buf4_close(&T2);
143
 
  dpd_buf4_sort_axpy(&Z, CC_HBAR, prqs, 11, 5, "WAbEi (Ei,Ab)", 1);
144
 
  dpd_buf4_close(&Z);
145
 
  dpd_buf4_close(&F);
146
 
 
147
 
 
148
 
}