← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cis/diag.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/cis/diag.c
1
 
#include <stdio.h>
2
 
#include <stdlib.h>
3
 
#include <libciomr/libciomr.h>
4
 
#include <libdpd/dpd.h>
5
 
#include <libqt/qt.h>
6
 
#include <physconst.h>
7
 
#define EXTERN
8
 
#include "globals.h"
9
 
 
10
 
void diag(void)
11
 
{
12
 
  int h, i, dim, dim_A, dim_B;
13
 
  int nroot, root;
14
 
  int ai, bj, ck, c, k, C, K, Ksym;
15
 
  double **A, *eps, **v;
16
 
  char lbl[32];
17
 
  dpdbuf4 A_AA, A_BB, A_AB;
18
 
  dpdfile2 B;
19
 
  double **singlet_evals, **triplet_evals, **uhf_evals;
20
 
 
21
 
  if(params.ref == 0) { /**RHF **/
22
 
 
23
 
    singlet_evals = (double **) malloc(moinfo.nirreps * sizeof(double *));
24
 
    for(h=0; h < moinfo.nirreps; h++)
25
 
      singlet_evals[h] = init_array(params.rpi[h]);
26
 
 
27
 
    dpd_buf4_init(&A_AA, CC_MISC, 0, 11, 11, 11, 11, 0, "A(AI,BJ)");
28
 
 
29
 
    for(h=0; h < moinfo.nirreps; h++) {
30
 
      dim = A_AA.params->rowtot[h];
31
 
      eps = init_array(dim);
32
 
      v = block_matrix(dim, dim);
33
 
 
34
 
      dpd_buf4_mat_irrep_init(&A_AA, h);
35
 
      dpd_buf4_mat_irrep_rd(&A_AA, h);
36
 
 
37
 
      if(!strcmp(params.diag_method,"FULL")) 
38
 
        sq_rsp(dim, dim, A_AA.matrix[h], eps, 1, v, 1e-14);
39
 
      else if(!strcmp(params.diag_method,"DAVIDSON")) {
40
 
        nroot = david(A_AA.matrix[h], dim, params.rpi[h], eps, v, params.convergence, 0);
41
 
 
42
 
        if(nroot != params.rpi[h])
43
 
          fprintf(outfile, "\tDavidson algorithm converged only %d roots in irrep %d.\n", 
44
 
                  nroot, h);
45
 
      }
46
 
 
47
 
      dpd_buf4_mat_irrep_close(&A_AA, h);
48
 
 
49
 
      /*
50
 
        fprintf(outfile, "RHF-CIS Singlet Eigenvectors for irrep %d:\n", h);
51
 
        print_mat(v, dim, params.rpi[h], outfile);
52
 
      */
53
 
 
54
 
      /* Store the eigenvectors in DPD entries and save the eigenvalues*/
55
 
      for(root=0; root < params.rpi[h]; root++) {
56
 
 
57
 
        sprintf(lbl, "BIA(%d)[%d] singlet", root, h);
58
 
        dpd_file2_init(&B, CC_OEI, h, 0, 1, lbl);
59
 
        dpd_file2_mat_init(&B);
60
 
        for(ck=0; ck < dim; ck++) {
61
 
          c = A_AA.params->roworb[h][ck][0];
62
 
          k = A_AA.params->roworb[h][ck][1];
63
 
 
64
 
          K = B.params->rowidx[k];
65
 
          C = B.params->colidx[c];
66
 
 
67
 
          Ksym = B.params->psym[k];
68
 
 
69
 
          B.matrix[Ksym][K][C] = v[ck][root];
70
 
        }
71
 
        dpd_file2_mat_wrt(&B);
72
 
        dpd_file2_mat_close(&B);
73
 
        dpd_file2_close(&B);
74
 
 
75
 
        singlet_evals[h][root] = eps[root];
76
 
      }
77
 
 
78
 
 
79
 
      free(eps);
80
 
      free_block(v);
81
 
    }
82
 
    dpd_buf4_close(&A_AA);
83
 
 
84
 
    triplet_evals = (double **) malloc(moinfo.nirreps * sizeof(double *));
85
 
    for(h=0; h < moinfo.nirreps; h++)
86
 
      triplet_evals[h] = init_array(params.rpi[h]);
87
 
 
88
 
    dpd_buf4_init(&A_AA, CC_MISC, 0, 11, 11, 11, 11, 0, "A(AI,BJ) triplet");
89
 
 
90
 
    for(h=0; h < moinfo.nirreps; h++) {
91
 
      dim = A_AA.params->rowtot[h];
92
 
      eps = init_array(dim);
93
 
      v = block_matrix(dim, dim);
94
 
 
95
 
      dpd_buf4_mat_irrep_init(&A_AA, h);
96
 
      dpd_buf4_mat_irrep_rd(&A_AA, h);
97
 
 
98
 
      if(!strcmp(params.diag_method, "FULL"))
99
 
        sq_rsp(dim, dim, A_AA.matrix[h], eps, 1, v, 1e-14);
100
 
      else if(!strcmp(params.diag_method,"DAVIDSON")) {
101
 
        nroot = david(A_AA.matrix[h], dim, params.rpi[h], eps, v, params.convergence, 0);
102
 
 
103
 
        if(nroot != params.rpi[h])
104
 
          fprintf(outfile, "\tDavidson algorithm converged only %d roots in irrep %d.\n", 
105
 
                  nroot, h);
106
 
      }
107
 
 
108
 
      dpd_buf4_mat_irrep_close(&A_AA, h);
109
 
 
110
 
      /* Store the eigenvectors in DPD entries and save the eigenvalues*/
111
 
      for(root=0; root < params.rpi[h]; root++) {
112
 
 
113
 
        sprintf(lbl, "BIA(%d)[%d] triplet", root, h);
114
 
        dpd_file2_init(&B, CC_OEI, h, 0, 1, lbl);
115
 
        dpd_file2_mat_init(&B);
116
 
        for(ck=0; ck < dim; ck++) {
117
 
          c = A_AA.params->roworb[h][ck][0];
118
 
          k = A_AA.params->roworb[h][ck][1];
119
 
 
120
 
          K = B.params->rowidx[k];
121
 
          C = B.params->colidx[c];
122
 
 
123
 
          Ksym = B.params->psym[k];
124
 
 
125
 
          B.matrix[Ksym][K][C] = v[ck][root];
126
 
        }
127
 
        dpd_file2_mat_wrt(&B);
128
 
        dpd_file2_mat_close(&B);
129
 
        dpd_file2_close(&B);
130
 
 
131
 
        triplet_evals[h][root] = eps[root];
132
 
      }
133
 
 
134
 
      free(eps);
135
 
      free_block(v);
136
 
    }
137
 
    dpd_buf4_close(&A_AA);
138
 
 
139
 
    moinfo.singlet_evals = singlet_evals;
140
 
    moinfo.triplet_evals = triplet_evals;
141
 
 
142
 
  }
143
 
  else if(params.ref == 2) { /** UHF **/
144
 
 
145
 
    uhf_evals = (double **) malloc(moinfo.nirreps * sizeof(double *));
146
 
    for(h=0; h < moinfo.nirreps; h++)
147
 
      uhf_evals[h] = init_array(params.rpi[h]);
148
 
 
149
 
    dpd_buf4_init(&A_AA, CC_MISC, 0, 21, 21, 21, 21, 0, "A(AI,BJ)");
150
 
    dpd_buf4_init(&A_BB, CC_MISC, 0, 31, 31, 31, 31, 0, "A(ai,bj)");
151
 
    dpd_buf4_init(&A_AB, CC_MISC, 0, 21, 31, 21, 31, 0, "A(AI,bj)");
152
 
    for(h=0; h < moinfo.nirreps; h++) {
153
 
      dim_A = A_AA.params->rowtot[h];
154
 
      dim_B = A_BB.params->rowtot[h];
155
 
 
156
 
      dim = dim_A + dim_B;
157
 
 
158
 
      A = block_matrix(dim, dim);
159
 
      eps = init_array(dim);
160
 
      v = block_matrix(dim, dim);
161
 
 
162
 
      dpd_buf4_mat_irrep_init(&A_AA, h);
163
 
      dpd_buf4_mat_irrep_rd(&A_AA, h);
164
 
      for(ai=0; ai < dim_A; ai++)
165
 
        for(bj=0; bj < dim_A; bj++)
166
 
          A[ai][bj] = A_AA.matrix[h][ai][bj];
167
 
      dpd_buf4_mat_irrep_close(&A_AA, h);
168
 
 
169
 
      dpd_buf4_mat_irrep_init(&A_BB, h);
170
 
      dpd_buf4_mat_irrep_rd(&A_BB, h);
171
 
      for(ai=0; ai < dim_B; ai++)
172
 
        for(bj=0; bj < dim_B; bj++)
173
 
          A[ai+dim_A][bj+dim_A] = A_BB.matrix[h][ai][bj];
174
 
      dpd_buf4_mat_irrep_close(&A_BB, h);
175
 
 
176
 
      dpd_buf4_mat_irrep_init(&A_AB, h);
177
 
      dpd_buf4_mat_irrep_rd(&A_AB, h);
178
 
      for(ai=0; ai < dim_A; ai++)
179
 
        for(bj=0; bj < dim_B; bj++)
180
 
          A[ai][bj+dim_A] = A[bj+dim_A][ai] = A_AB.matrix[h][ai][bj];
181
 
      dpd_buf4_mat_irrep_close(&A_AB, h);
182
 
 
183
 
      if(!strcmp(params.diag_method,"FULL"))
184
 
        sq_rsp(dim, dim, A, eps, 1, v, 1e-12);
185
 
      else if(!strcmp(params.diag_method,"DAVIDSON")) {
186
 
        nroot = david(A, dim, params.rpi[h], eps, v, params.convergence, 0);
187
 
 
188
 
        if(nroot != params.rpi[h])
189
 
          fprintf(outfile, "\tDavidson algorithm converged only %d roots in irrep %d.\n", 
190
 
                  nroot, h);
191
 
      }
192
 
 
193
 
      /*
194
 
        fprintf(outfile, "UHF-CIS Eigenvectors for irrep %d:\n", h);
195
 
        print_mat(v, dim, params.rpi[h], outfile);
196
 
      */
197
 
 
198
 
      /* Store the eigenvectors in DPD entries and save the eigenvalues */
199
 
      for(root=0; root < params.rpi[h]; root++) {
200
 
        sprintf(lbl, "BIA(%d)[%d]", root, h);
201
 
        dpd_file2_init(&B, CC_OEI, h, 0, 1, lbl);
202
 
        dpd_file2_mat_init(&B);
203
 
        for(ck=0; ck < dim_A; ck++) {
204
 
          c = A_AA.params->roworb[h][ck][0];
205
 
          k = A_AA.params->roworb[h][ck][1];
206
 
 
207
 
          K = B.params->rowidx[k];
208
 
          C = B.params->colidx[c];
209
 
 
210
 
          Ksym = B.params->psym[k];
211
 
 
212
 
          B.matrix[Ksym][K][C] = v[ck][root];
213
 
        }
214
 
        dpd_file2_mat_wrt(&B);
215
 
        dpd_file2_mat_close(&B);
216
 
        dpd_file2_close(&B);
217
 
 
218
 
        sprintf(lbl, "Bia(%d)[%d]", root, h);
219
 
        dpd_file2_init(&B, CC_OEI, h, 2, 3, lbl);
220
 
        dpd_file2_mat_init(&B);
221
 
        for(ck=0; ck < dim_B; ck++) {
222
 
          c = A_BB.params->roworb[h][ck][0];
223
 
          k = A_BB.params->roworb[h][ck][1];
224
 
 
225
 
          K = B.params->rowidx[k];
226
 
          C = B.params->colidx[c];
227
 
 
228
 
          Ksym = B.params->psym[k];
229
 
 
230
 
          B.matrix[Ksym][K][C] = v[ck+dim_A][root];
231
 
        }
232
 
        dpd_file2_mat_wrt(&B);
233
 
        dpd_file2_mat_close(&B);
234
 
        dpd_file2_close(&B);
235
 
 
236
 
        uhf_evals[h][root] = eps[root];
237
 
      }
238
 
 
239
 
      free(eps);
240
 
      free_block(v);
241
 
      free_block(A);
242
 
    }
243
 
    dpd_buf4_close(&A_AA);
244
 
    dpd_buf4_close(&A_BB);
245
 
    dpd_buf4_close(&A_AB);
246
 
 
247
 
    moinfo.uhf_evals = uhf_evals;
248
 
  }
249
 
 
250
 
}