← Back to branch summary

~ubuntu-branches/ubuntu/vivid/psicode/vivid

« back to all changes in this revision

Viewing changes to src/bin/ccsort/cphf_B.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:
src/bin/ccsort/cphf_B.cc
 
1
/*! \file
 
2
    \ingroup CCSORT
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#include <cstdio>
 
6
#include <cstdlib>
 
7
#include <cstring>
 
8
#include <libqt/qt.h>
 
9
#include <libdpd/dpd.h>
 
10
#include <libciomr/libciomr.h>
 
11
#include <psifiles.h>
 
12
#include "MOInfo.h"
 
13
#define EXTERN
 
14
#include "globals.h"
 
15
 
 
16
namespace psi { namespace ccsort {
 
17
 
 
18
void cphf_B(const char *cart)
 
19
{
 
20
  int irrep, row, a, i, asym, isym, num_ai, info, *ipiv;
 
21
  double *vector, *vector2, polar;
 
22
  dpdbuf4 A;
 
23
  dpdfile2 L;
 
24
 
 
25
  psio_open(PSIF_MO_HESS, 1);
 
26
  dpd_buf4_init(&A, PSIF_MO_HESS, 0, 11, 11, 11, 11, 0, "B(AI,BJ)");
 
27
 
 
28
  if (!strcmp(cart, "X")) {
 
29
    irrep = moinfo.irrep_x;
 
30
 
 
31
    /* sort L elements into a single vector for lineq solver */
 
32
    dpd_file2_init(&L, CC_OEI, irrep, 0, 1, "L_X_IA");
 
33
    dpd_file2_mat_init(&L);
 
34
    dpd_file2_mat_rd(&L);
 
35
    num_ai = A.params->rowtot[irrep];
 
36
    vector = init_array(num_ai);
 
37
    for(row=0; row < num_ai; row++) {
 
38
      a = A.params->roworb[irrep][row][0];
 
39
      i = A.params->roworb[irrep][row][1];
 
40
      asym = A.params->psym[a];
 
41
      isym = A.params->qsym[i];
 
42
      vector[row] = L.matrix[isym][i-A.params->qoff[isym]][a-A.params->poff[asym]];
 
43
    }
 
44
    dpd_file2_mat_close(&L);
 
45
    dpd_file2_close(&L);
 
46
 
 
47
    /* grab current irrep of MO Hessian */
 
48
    dpd_buf4_mat_irrep_init(&A, irrep);
 
49
    dpd_buf4_mat_irrep_rd(&A, irrep);
 
50
 
 
51
 
 
52
    /* solve CPHF equations */
 
53
    ipiv = init_int_array(num_ai);
 
54
    info = C_DGESV(num_ai, 1, &(A.matrix[irrep][0][0]), num_ai, ipiv, vector, num_ai);
 
55
    if(info) {
 
56
      fprintf(outfile, "CCSORT: cphf_B: Error in C_DGESV.  Info = %d.  Exiting.\n", info);
 
57
      exit(PSI_RETURN_FAILURE);
 
58
    }
 
59
    free(ipiv);
 
60
 
 
61
    dpd_buf4_mat_irrep_close(&A, irrep);
 
62
 
 
63
    /* sort CPHF solution to DPD format */
 
64
    dpd_file2_init(&L, CC_OEI, irrep, 1, 0, "CPHF Ub_X_AI");
 
65
    dpd_file2_mat_init(&L);
 
66
    for(row=0; row < num_ai; row++) {
 
67
      a = A.params->roworb[irrep][row][0];
 
68
      i = A.params->roworb[irrep][row][1];
 
69
      asym = A.params->psym[a];
 
70
      isym = A.params->qsym[i];
 
71
      L.matrix[asym][a-A.params->poff[asym]][i-A.params->qoff[isym]] = vector[row];
 
72
    }
 
73
    dpd_file2_mat_wrt(&L);
 
74
    dpd_file2_close(&L);
 
75
  }
 
76
 
 
77
  if (!strcmp(cart, "Y")) {
 
78
    irrep = moinfo.irrep_y;
 
79
 
 
80
    /* sort L elements into a single vector for lineq solver */
 
81
    dpd_file2_init(&L, CC_OEI, irrep, 0, 1, "L_Y_IA");
 
82
    dpd_file2_mat_init(&L);
 
83
    dpd_file2_mat_rd(&L);
 
84
    num_ai = A.params->rowtot[irrep];
 
85
    vector = init_array(num_ai);
 
86
    for(row=0; row < num_ai; row++) {
 
87
      a = A.params->roworb[irrep][row][0];
 
88
      i = A.params->roworb[irrep][row][1];
 
89
      asym = A.params->psym[a];
 
90
      isym = A.params->qsym[i];
 
91
      vector[row] = L.matrix[isym][i-A.params->qoff[isym]][a-A.params->poff[asym]];
 
92
    }
 
93
    dpd_file2_mat_close(&L);
 
94
    dpd_file2_close(&L);
 
95
 
 
96
    /* grab current irrep of MO Hessian */
 
97
    dpd_buf4_mat_irrep_init(&A, irrep);
 
98
    dpd_buf4_mat_irrep_rd(&A, irrep);
 
99
 
 
100
 
 
101
    /* solve CPHF equations */
 
102
    ipiv = init_int_array(num_ai);
 
103
    info = C_DGESV(num_ai, 1, &(A.matrix[irrep][0][0]), num_ai, ipiv, vector, num_ai);
 
104
    if(info) {
 
105
      fprintf(outfile, "CCSORT: cphf_B: Error in C_DGESV.  Info = %d.  Exiting.\n", info);
 
106
      exit(PSI_RETURN_FAILURE);
 
107
    }
 
108
    free(ipiv);
 
109
 
 
110
    dpd_buf4_mat_irrep_close(&A, irrep);
 
111
 
 
112
    /* sort CPHF solution to DPD format */
 
113
    dpd_file2_init(&L, CC_OEI, irrep, 1, 0, "CPHF Ub_Y_AI");
 
114
    dpd_file2_mat_init(&L);
 
115
    for(row=0; row < num_ai; row++) {
 
116
      a = A.params->roworb[irrep][row][0];
 
117
      i = A.params->roworb[irrep][row][1];
 
118
      asym = A.params->psym[a];
 
119
      isym = A.params->qsym[i];
 
120
      L.matrix[asym][a-A.params->poff[asym]][i-A.params->qoff[isym]] = vector[row];
 
121
    }
 
122
    dpd_file2_mat_wrt(&L);
 
123
    dpd_file2_close(&L);
 
124
  }
 
125
 
 
126
  if (!strcmp(cart, "Z")) {
 
127
    irrep = moinfo.irrep_z;
 
128
 
 
129
    /* sort L elements into a single vector for lineq solver */
 
130
    dpd_file2_init(&L, CC_OEI, irrep, 0, 1, "L_Z_IA");
 
131
    dpd_file2_mat_init(&L);
 
132
    dpd_file2_mat_rd(&L);
 
133
    num_ai = A.params->rowtot[irrep];
 
134
    vector = init_array(num_ai);
 
135
    for(row=0; row < num_ai; row++) {
 
136
      a = A.params->roworb[irrep][row][0];
 
137
      i = A.params->roworb[irrep][row][1];
 
138
      asym = A.params->psym[a];
 
139
      isym = A.params->qsym[i];
 
140
      vector[row] = L.matrix[isym][i-A.params->qoff[isym]][a-A.params->poff[asym]];
 
141
    }
 
142
    dpd_file2_mat_close(&L);
 
143
    dpd_file2_close(&L);
 
144
 
 
145
    /*
 
146
      vector2 = init_array(num_ai);
 
147
      for(row=0; row < num_ai; row++) vector2[row] = vector[row];
 
148
    */
 
149
 
 
150
    /* grab current irrep of MO Hessian */
 
151
    dpd_buf4_mat_irrep_init(&A, irrep);
 
152
    dpd_buf4_mat_irrep_rd(&A, irrep);
 
153
 
 
154
    /* solve CPHF equations */
 
155
    ipiv = init_int_array(num_ai);
 
156
    info = C_DGESV(num_ai, 1, &(A.matrix[irrep][0][0]), num_ai, ipiv, vector, num_ai);
 
157
    if(info) {
 
158
      fprintf(outfile, "CCSORT: cphf_B: Error in C_DGESV.  Info = %d.  Exiting.\n", info);
 
159
      exit(PSI_RETURN_FAILURE);
 
160
    }
 
161
    free(ipiv);
 
162
 
 
163
    dpd_buf4_mat_irrep_close(&A, irrep);
 
164
 
 
165
    /*
 
166
      polar = 0.0;
 
167
      for(row=0; row < num_ai; row++) polar += vector2[row] * vector[row];
 
168
      polar *= 4.0;
 
169
      fprintf(outfile, "\talpha_zz = %20.12f\n", polar);
 
170
      free(vector2);
 
171
    */
 
172
 
 
173
    /* sort CPHF solution to DPD format */
 
174
    dpd_file2_init(&L, CC_OEI, irrep, 1, 0, "CPHF Ub_Z_AI");
 
175
    dpd_file2_mat_init(&L);
 
176
    for(row=0; row < num_ai; row++) {
 
177
      a = A.params->roworb[irrep][row][0];
 
178
      i = A.params->roworb[irrep][row][1];
 
179
      asym = A.params->psym[a];
 
180
      isym = A.params->qsym[i];
 
181
      L.matrix[asym][a-A.params->poff[asym]][i-A.params->qoff[isym]] = vector[row];
 
182
    }
 
183
    dpd_file2_mat_wrt(&L);
 
184
    dpd_file2_close(&L);
 
185
  }
 
186
 
 
187
  dpd_buf4_close(&A);
 
188
  if (!strcmp(cart,"Z"))
 
189
    psio_close(PSIF_MO_HESS, 0);
 
190
  else
 
191
    psio_close(PSIF_MO_HESS, 1);
 
192
}
 
193
 
 
194
}} // namespace psi::ccsort