← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccdensity/build_A_UHF.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2006-09-10 14:01:33 UTC
  • Revision ID: james.westby@ubuntu.com-20060910140133-ib2j86trekykfsfv
Tags: upstream-3.2.3
Import upstream version 3.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/ccdensity/build_A_UHF.c
 
1
#include <libdpd/dpd.h>
 
2
#define EXTERN
 
3
#include "globals.h"
 
4
 
 
5
/* BUILD_A_UHF(): Construct the components of the molecular orbital
 
6
** Hessian, A, for UHF orbitals. At the moment we're actually building
 
7
** all symmetry blocks of A, though for the orbital Z-vector equations
 
8
** we really only need the totally symmetric components.
 
9
**
 
10
** A(AI,JB) = delta_IJ f_AB - delta_AB f_IJ + <IJ||AB> - <JA||IB> 
 
11
**
 
12
** A(ai,jb) = delta_ij f_ab - delta_ab f_ij + <ij||ab> - <ja||ib> 
 
13
**
 
14
** A(AI,jb) = 2<Ij|Ab>
 
15
**
 
16
** TDC, January 2003
 
17
*/
 
18
 
 
19
void build_A_UHF(void)
 
20
{
 
21
  int h, nirreps, row, col;
 
22
  int a, i, b, j;
 
23
  int A, I, B, J;
 
24
  int Asym, Isym, Bsym, Jsym;
 
25
  dpdfile2 fIJ, fij, fAB, fab, fIA, fia;
 
26
  dpdbuf4 Amat, D, C;
 
27
 
 
28
  nirreps = moinfo.nirreps;
 
29
 
 
30
  dpd_file2_init(&fIJ, CC_OEI, 0, 0, 0, "fIJ");
 
31
  dpd_file2_mat_init(&fIJ);
 
32
  dpd_file2_mat_rd(&fIJ);
 
33
  dpd_file2_init(&fij, CC_OEI, 0, 2, 2, "fij");
 
34
  dpd_file2_mat_init(&fij);
 
35
  dpd_file2_mat_rd(&fij);
 
36
  dpd_file2_init(&fAB, CC_OEI, 0, 1, 1, "fAB");
 
37
  dpd_file2_mat_init(&fAB);
 
38
  dpd_file2_mat_rd(&fAB);
 
39
  dpd_file2_init(&fab, CC_OEI, 0, 3, 3, "fab");
 
40
  dpd_file2_mat_init(&fab);
 
41
  dpd_file2_mat_rd(&fab);
 
42
  dpd_file2_init(&fIA, CC_OEI, 0, 0, 1, "fIA");
 
43
  dpd_file2_mat_init(&fIA);
 
44
  dpd_file2_mat_rd(&fIA);
 
45
  dpd_file2_init(&fia, CC_OEI, 0, 2, 3, "fia");
 
46
  dpd_file2_mat_init(&fia);
 
47
  dpd_file2_mat_rd(&fia);
 
48
 
 
49
  dpd_buf4_init(&D, CC_DINTS, 0, 0, 5, 0, 5, 1, "D <IJ|AB>");
 
50
  dpd_buf4_sort(&D, CC_MISC, rpsq, 21, 21, "A(AI,BJ)");
 
51
  dpd_buf4_close(&D);
 
52
  dpd_buf4_init(&D, CC_CINTS, 0, 20, 20, 20, 20, 0, "C <IA||JB>");
 
53
  dpd_buf4_sort_axpy(&D, CC_MISC, qrsp, 21, 21, "A(AI,BJ)", -1.0);
 
54
  dpd_buf4_close(&D);
 
55
 
 
56
  dpd_buf4_init(&Amat, CC_MISC, 0, 21, 21, 21, 21, 0, "A(AI,BJ)");
 
57
  for(h=0; h < nirreps; h++) {
 
58
    dpd_buf4_mat_irrep_init(&Amat, h);
 
59
    dpd_buf4_mat_irrep_rd(&Amat, h);
 
60
 
 
61
    for(row=0; row < Amat.params->rowtot[h]; row++) {
 
62
      a = Amat.params->roworb[h][row][0];
 
63
      i = Amat.params->roworb[h][row][1];
 
64
      A = fAB.params->rowidx[a];
 
65
      I = fIJ.params->rowidx[i];
 
66
      Asym = fAB.params->psym[a];
 
67
      Isym = fIJ.params->psym[i];
 
68
      for(col=0; col < Amat.params->coltot[h]; col++) {
 
69
        b = Amat.params->colorb[h][col][0];
 
70
        j = Amat.params->colorb[h][col][1];
 
71
        B = fAB.params->colidx[b];
 
72
        J = fIJ.params->colidx[j];
 
73
        Bsym = fAB.params->qsym[b];
 
74
        Jsym = fIJ.params->qsym[j];
 
75
 
 
76
        if((I==J) && (Asym==Bsym))
 
77
          Amat.matrix[h][row][col] += fAB.matrix[Asym][A][B];
 
78
        if((A==B) && (Isym==Jsym))
 
79
          Amat.matrix[h][row][col] -= fIJ.matrix[Isym][I][J];
 
80
 
 
81
      }
 
82
    }
 
83
 
 
84
    dpd_buf4_mat_irrep_wrt(&Amat, h);
 
85
    dpd_buf4_mat_irrep_close(&Amat, h);
 
86
  }
 
87
  dpd_buf4_close(&Amat);
 
88
 
 
89
  dpd_buf4_init(&D, CC_DINTS, 0, 10, 15, 10, 15, 1, "D <ij|ab>");
 
90
  dpd_buf4_sort(&D, CC_MISC, rpsq, 31, 31, "A(ai,bj)");
 
91
  dpd_buf4_close(&D);
 
92
  dpd_buf4_init(&D, CC_CINTS, 0, 30, 30, 30, 30, 0, "C <ia||jb>");
 
93
  dpd_buf4_sort_axpy(&D, CC_MISC, qrsp, 31, 31, "A(ai,bj)", -1.0);
 
94
  dpd_buf4_close(&D);
 
95
 
 
96
  dpd_buf4_init(&Amat, CC_MISC, 0, 31, 31, 31, 31, 0, "A(ai,bj)");
 
97
  for(h=0; h < nirreps; h++) {
 
98
    dpd_buf4_mat_irrep_init(&Amat, h);
 
99
    dpd_buf4_mat_irrep_rd(&Amat, h);
 
100
 
 
101
    for(row=0; row < Amat.params->rowtot[h]; row++) {
 
102
      a = Amat.params->roworb[h][row][0];
 
103
      i = Amat.params->roworb[h][row][1];
 
104
      A = fab.params->rowidx[a];
 
105
      I = fij.params->rowidx[i];
 
106
      Asym = fab.params->psym[a];
 
107
      Isym = fij.params->psym[i];
 
108
      for(col=0; col < Amat.params->coltot[h]; col++) {
 
109
        b = Amat.params->colorb[h][col][0];
 
110
        j = Amat.params->colorb[h][col][1];
 
111
        B = fab.params->colidx[b];
 
112
        J = fij.params->colidx[j];
 
113
        Bsym = fab.params->qsym[b];
 
114
        Jsym = fij.params->qsym[j];
 
115
 
 
116
        if((I==J) && (Asym==Bsym))
 
117
          Amat.matrix[h][row][col] += fab.matrix[Asym][A][B];
 
118
        if((A==B) && (Isym==Jsym))
 
119
          Amat.matrix[h][row][col] -= fij.matrix[Isym][I][J];
 
120
 
 
121
      }
 
122
    }
 
123
 
 
124
    dpd_buf4_mat_irrep_wrt(&Amat, h);
 
125
    dpd_buf4_mat_irrep_close(&Amat, h);
 
126
  }
 
127
  dpd_buf4_close(&Amat);
 
128
 
 
129
  dpd_buf4_init(&D, CC_DINTS, 0, 22, 28, 22, 28, 0, "D <Ij|Ab>");
 
130
  dpd_buf4_sort(&D, CC_MISC, rpsq, 21, 31, "A(AI,bj)");
 
131
  dpd_buf4_close(&D);
 
132
  dpd_buf4_init(&Amat, CC_MISC, 0, 21, 31, 21, 31, 0, "A(AI,bj)");
 
133
  dpd_buf4_scm(&Amat, 2);
 
134
  dpd_buf4_close(&Amat);
 
135
 
 
136
  dpd_file2_mat_close(&fIJ);
 
137
  dpd_file2_close(&fIJ);
 
138
  dpd_file2_mat_close(&fij);
 
139
  dpd_file2_close(&fij);
 
140
  dpd_file2_mat_close(&fAB);
 
141
  dpd_file2_close(&fAB);
 
142
  dpd_file2_mat_close(&fab);
 
143
  dpd_file2_close(&fab);
 
144
  dpd_file2_mat_close(&fIA);
 
145
  dpd_file2_close(&fIA);
 
146
  dpd_file2_mat_close(&fia);
 
147
  dpd_file2_close(&fia);
 
148
}
 
149