← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cceom/precondition.c

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck
  • Date: 2008-06-07 16:49:57 UTC
  • mfrom: (2.1.2 hardy)
  • Revision ID: james.westby@ubuntu.com-20080607164957-8pifvb133yjlkagn
Tags: 3.3.0-3
http://bugs.debian.org/464867
* debian/rules (DEB_MAKE_CHECK_TARGET): Do not abort test suite on
  failures.
* debian/rules (DEB_CONFIGURE_EXTRA_FLAGS): Set ${bindir} to /usr/lib/psi.
* debian/rules (install/psi3): Move psi3 file to /usr/bin.
* debian/patches/07_464867_move_executables.dpatch: New patch, add
  /usr/lib/psi to the $PATH, so that the moved executables are found.
  (closes: #464867)
* debian/patches/00list: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/cceom/precondition.c
 
1
#include <stdlib.h>
1
2
#include <stdio.h>
2
3
#include <math.h>
 
4
#include <string.h>
 
5
#include <libqt/qt.h>
3
6
#define EXTERN
4
7
#include "globals.h"
5
8
 
123
126
 
124
127
  /* Local correlation decs */
125
128
  int nso, nocc, nvir;
126
 
  int *pairdom_len, *pairdom_nrlen, *weak_pairs;
127
 
  double ***V, ***W, *eps_occ, **eps_vir;
128
129
  double *T1tilde, *T1bar, **T2tilde, **T2bar, **X1, **X2;
 
130
  psio_address next;
129
131
 
130
132
  C_irr = RIA->my_irrep;
131
133
  nirreps = RIA->params->nirreps;
132
134
 
133
 
  if(params.local && local.filter_singles) {
 
135
  if(params.local) {
134
136
 
135
137
    nso = local.nso;
136
138
    nocc = local.nocc;
137
139
    nvir = local.nvir;
138
 
    V = local.V;
139
 
    W = local.W;
140
 
    eps_occ = local.eps_occ;
141
 
    eps_vir = local.eps_vir;
142
 
    pairdom_len = local.pairdom_len;
143
 
    pairdom_nrlen = local.pairdom_nrlen;
144
 
    weak_pairs = local.weak_pairs;
145
 
 
146
 
    dpd_file2_mat_init(RIA);
147
 
    dpd_file2_mat_rd(RIA);
148
 
 
149
 
    for(i=0; i < nocc; i++) {
150
 
      ii = i * nocc +i;
151
 
 
152
 
      if(!pairdom_len[ii]) {
153
 
        fprintf(outfile, "\n\tlocal_filter_T1: Pair ii = [%d] is zero-length, which makes no sense.\n",ii);
154
 
        exit(2);
155
 
      }
156
 
 
157
 
      T1tilde = init_array(pairdom_len[ii]);
158
 
      T1bar = init_array(pairdom_nrlen[ii]);
159
 
 
160
 
      /* Transform the virtuals to the redundant projected virtual basis */
161
 
      C_DGEMV('t', nvir, pairdom_len[ii], 1.0, &(V[ii][0][0]), pairdom_len[ii], 
162
 
              &(RIA->matrix[0][i][0]), 1, 0.0, &(T1tilde[0]), 1);
163
 
 
164
 
      /* Transform the virtuals to the non-redundant virtual basis */
165
 
      C_DGEMV('t', pairdom_len[ii], pairdom_nrlen[ii], 1.0, &(W[ii][0][0]), pairdom_nrlen[ii], 
166
 
              &(T1tilde[0]), 1, 0.0, &(T1bar[0]), 1);
167
 
 
168
 
      for(a=0; a < pairdom_nrlen[ii]; a++) {
169
 
        tval = eval + eps_occ[i] - eps_vir[ii][a];
170
 
        if(fabs(tval) > 0.0001) T1bar[a] /= tval;
171
 
      }
172
 
 
173
 
      /* Transform the new T1's to the redundant projected virtual basis */
174
 
      C_DGEMV('n', pairdom_len[ii], pairdom_nrlen[ii], 1.0, &(W[ii][0][0]), pairdom_nrlen[ii],
175
 
              &(T1bar[0]), 1, 0.0, &(T1tilde[0]), 1);
176
 
 
177
 
      /* Transform the new T1's to the MO basis */
178
 
      C_DGEMV('n', nvir, pairdom_len[ii], 1.0, &(V[ii][0][0]), pairdom_len[ii], 
179
 
              &(T1tilde[0]), 1, 0.0, &(RIA->matrix[0][i][0]), 1);
180
 
 
181
 
      free(T1bar);
182
 
      free(T1tilde);
183
 
    }
184
 
 
185
 
    dpd_file2_mat_wrt(RIA);
186
 
    dpd_file2_mat_close(RIA);
 
140
 
 
141
    local.pairdom_len = init_int_array(nocc*nocc);
 
142
    local.pairdom_nrlen = init_int_array(nocc*nocc);
 
143
    local.eps_occ = init_array(nocc);
 
144
    local.weak_pairs = init_int_array(nocc*nocc);
 
145
    psio_read_entry(CC_INFO, "Local Pair Domain Length", (char *) local.pairdom_len,
 
146
                    nocc*nocc*sizeof(int));
 
147
    psio_read_entry(CC_INFO, "Local Pair Domain Length (Non-redundant basis)", (char *) local.pairdom_nrlen,
 
148
                    nocc*nocc*sizeof(int));
 
149
    psio_read_entry(CC_INFO, "Local Occupied Orbital Energies", (char *) local.eps_occ,
 
150
                    nocc*sizeof(double));
 
151
    psio_read_entry(CC_INFO, "Local Weak Pairs", (char *) local.weak_pairs,
 
152
                    nocc*nocc*sizeof(int));
 
153
 
 
154
    local.W = (double ***) malloc(nocc * nocc * sizeof(double **));
 
155
    local.V = (double ***) malloc(nocc * nocc * sizeof(double **));
 
156
    local.eps_vir = (double **) malloc(nocc * nocc * sizeof(double *));
 
157
    next = PSIO_ZERO;
 
158
    for(ij=0; ij < nocc*nocc; ij++) {
 
159
      local.eps_vir[ij] = init_array(local.pairdom_nrlen[ij]);
 
160
      psio_read(CC_INFO, "Local Virtual Orbital Energies", (char *) local.eps_vir[ij],
 
161
                local.pairdom_nrlen[ij]*sizeof(double), next, &next);
 
162
    }
 
163
    next = PSIO_ZERO;
 
164
    for(ij=0; ij < nocc*nocc; ij++) {
 
165
      local.V[ij] = block_matrix(nvir,local.pairdom_len[ij]);
 
166
      psio_read(CC_INFO, "Local Residual Vector (V)", (char *) local.V[ij][0],
 
167
                nvir*local.pairdom_len[ij]*sizeof(double), next, &next);
 
168
    }
 
169
    next = PSIO_ZERO;
 
170
    for(ij=0; ij < nocc*nocc; ij++) {
 
171
      local.W[ij] = block_matrix(local.pairdom_len[ij],local.pairdom_nrlen[ij]);
 
172
      psio_read(CC_INFO, "Local Transformation Matrix (W)", (char *) local.W[ij][0],
 
173
                local.pairdom_len[ij]*local.pairdom_nrlen[ij]*sizeof(double), next, &next);
 
174
    }
 
175
 
 
176
    if(local.filter_singles) {
 
177
 
 
178
      dpd_file2_mat_init(RIA);
 
179
      dpd_file2_mat_rd(RIA);
 
180
 
 
181
      for(i=0; i < nocc; i++) {
 
182
        ii = i * nocc +i;
 
183
 
 
184
        if(!local.pairdom_len[ii]) {
 
185
          fprintf(outfile, "\n\tlocal_filter_T1: Pair ii = [%d] is zero-length, which makes no sense.\n",ii);
 
186
          exit(2);
 
187
        }
 
188
 
 
189
        T1tilde = init_array(local.pairdom_len[ii]);
 
190
        T1bar = init_array(local.pairdom_nrlen[ii]);
 
191
 
 
192
        /* Transform the virtuals to the redundant projected virtual basis */
 
193
        C_DGEMV('t', nvir, local.pairdom_len[ii], 1.0, &(local.V[ii][0][0]), local.pairdom_len[ii], 
 
194
                &(RIA->matrix[0][i][0]), 1, 0.0, &(T1tilde[0]), 1);
 
195
 
 
196
        /* Transform the virtuals to the non-redundant virtual basis */
 
197
        C_DGEMV('t', local.pairdom_len[ii], local.pairdom_nrlen[ii], 1.0, &(local.W[ii][0][0]), local.pairdom_nrlen[ii], 
 
198
                &(T1tilde[0]), 1, 0.0, &(T1bar[0]), 1);
 
199
 
 
200
        for(a=0; a < local.pairdom_nrlen[ii]; a++) {
 
201
          tval = eval + local.eps_occ[i] - local.eps_vir[ii][a];
 
202
          if(fabs(tval) > 0.0001) T1bar[a] /= tval;
 
203
        }
 
204
 
 
205
        /* Transform the new T1's to the redundant projected virtual basis */
 
206
        C_DGEMV('n', local.pairdom_len[ii], local.pairdom_nrlen[ii], 1.0, &(local.W[ii][0][0]), local.pairdom_nrlen[ii],
 
207
                &(T1bar[0]), 1, 0.0, &(T1tilde[0]), 1);
 
208
 
 
209
        /* Transform the new T1's to the MO basis */
 
210
        C_DGEMV('n', nvir, local.pairdom_len[ii], 1.0, &(local.V[ii][0][0]), local.pairdom_len[ii], 
 
211
                &(T1tilde[0]), 1, 0.0, &(RIA->matrix[0][i][0]), 1);
 
212
 
 
213
        free(T1bar);
 
214
        free(T1tilde);
 
215
      }
 
216
 
 
217
      dpd_file2_mat_wrt(RIA);
 
218
      dpd_file2_mat_close(RIA);
 
219
    }
187
220
  }
188
221
  else {
189
222
    dpd_file2_mat_init(RIA);
206
239
 
207
240
  if(params.local) {
208
241
 
209
 
    nso = local.nso;
210
 
    nocc = local.nocc;
211
 
    nvir = local.nvir;
212
 
    V = local.V;
213
 
    W = local.W;
214
 
    eps_occ = local.eps_occ;
215
 
    eps_vir = local.eps_vir;
216
 
    pairdom_len = local.pairdom_len;
217
 
    pairdom_nrlen = local.pairdom_nrlen;
218
 
    weak_pairs = local.weak_pairs;
219
 
 
220
242
    dpd_buf4_mat_irrep_init(RIjAb, 0);
221
243
    dpd_buf4_mat_irrep_rd(RIjAb, 0);
222
244
 
228
250
    for(i=0,ij=0; i < nocc; i++) {
229
251
      for(j=0; j < nocc; j++,ij++) {
230
252
 
231
 
        if(!weak_pairs[ij]) {
 
253
        if(!local.weak_pairs[ij]) {
232
254
 
233
255
          /* Transform the virtuals to the redundant projected virtual basis */
234
 
          C_DGEMM('t', 'n', pairdom_len[ij], nvir, nvir, 1.0, &(V[ij][0][0]), pairdom_len[ij],
 
256
          C_DGEMM('t', 'n', local.pairdom_len[ij], nvir, nvir, 1.0, &(local.V[ij][0][0]), local.pairdom_len[ij],
235
257
                  &(RIjAb->matrix[0][ij][0]), nvir, 0.0, &(X1[0][0]), nvir);
236
 
          C_DGEMM('n', 'n', pairdom_len[ij], pairdom_len[ij], nvir, 1.0, &(X1[0][0]), nvir,
237
 
                  &(V[ij][0][0]), pairdom_len[ij], 0.0, &(T2tilde[0][0]), nso);
 
258
          C_DGEMM('n', 'n', local.pairdom_len[ij], local.pairdom_len[ij], nvir, 1.0, &(X1[0][0]), nvir,
 
259
                  &(local.V[ij][0][0]), local.pairdom_len[ij], 0.0, &(T2tilde[0][0]), nso);
238
260
 
239
261
          /* Transform the virtuals to the non-redundant virtual basis */
240
 
          C_DGEMM('t', 'n', pairdom_nrlen[ij], pairdom_len[ij], pairdom_len[ij], 1.0, 
241
 
                  &(W[ij][0][0]), pairdom_nrlen[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
242
 
          C_DGEMM('n', 'n', pairdom_nrlen[ij], pairdom_nrlen[ij], pairdom_len[ij], 1.0, 
243
 
                  &(X2[0][0]), nso, &(W[ij][0][0]), pairdom_nrlen[ij], 0.0, &(T2bar[0][0]), nvir);
 
262
          C_DGEMM('t', 'n', local.pairdom_nrlen[ij], local.pairdom_len[ij], local.pairdom_len[ij], 1.0, 
 
263
                  &(local.W[ij][0][0]), local.pairdom_nrlen[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
 
264
          C_DGEMM('n', 'n', local.pairdom_nrlen[ij], local.pairdom_nrlen[ij], local.pairdom_len[ij], 1.0, 
 
265
                  &(X2[0][0]), nso, &(local.W[ij][0][0]), local.pairdom_nrlen[ij], 0.0, &(T2bar[0][0]), nvir);
244
266
 
245
267
          /* Divide the new amplitudes by the denominators */
246
 
          for(a=0; a < pairdom_nrlen[ij]; a++) {
247
 
            for(b=0; b < pairdom_nrlen[ij]; b++) {
248
 
              tval = eval + eps_occ[i]+eps_occ[j]-eps_vir[ij][a]-eps_vir[ij][b];
 
268
          for(a=0; a < local.pairdom_nrlen[ij]; a++) {
 
269
            for(b=0; b < local.pairdom_nrlen[ij]; b++) {
 
270
              tval = eval + local.eps_occ[i]+local.eps_occ[j]-local.eps_vir[ij][a]-local.eps_vir[ij][b];
249
271
              if(fabs(tval) > 0.0001) T2bar[a][b] /= tval;
250
272
            }
251
273
          }
252
274
 
253
275
          /* Transform the new T2's to the redundant virtual basis */
254
 
          C_DGEMM('n', 'n', pairdom_len[ij], pairdom_nrlen[ij], pairdom_nrlen[ij], 1.0, 
255
 
                  &(W[ij][0][0]), pairdom_nrlen[ij], &(T2bar[0][0]), nvir, 0.0, &(X1[0][0]), nvir);
256
 
          C_DGEMM('n','t', pairdom_len[ij], pairdom_len[ij], pairdom_nrlen[ij], 1.0, 
257
 
                  &(X1[0][0]), nvir, &(W[ij][0][0]), pairdom_nrlen[ij], 0.0, &(T2tilde[0][0]), nso);
 
276
          C_DGEMM('n', 'n', local.pairdom_len[ij], local.pairdom_nrlen[ij], local.pairdom_nrlen[ij], 1.0, 
 
277
                  &(local.W[ij][0][0]), local.pairdom_nrlen[ij], &(T2bar[0][0]), nvir, 0.0, &(X1[0][0]), nvir);
 
278
          C_DGEMM('n','t', local.pairdom_len[ij], local.pairdom_len[ij], local.pairdom_nrlen[ij], 1.0, 
 
279
                  &(X1[0][0]), nvir, &(local.W[ij][0][0]), local.pairdom_nrlen[ij], 0.0, &(T2tilde[0][0]), nso);
258
280
 
259
281
          /* Transform the new T2's to the MO basis */
260
 
          C_DGEMM('n', 'n', nvir, pairdom_len[ij], pairdom_len[ij], 1.0, 
261
 
                  &(V[ij][0][0]), pairdom_len[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
262
 
          C_DGEMM('n', 't', nvir, nvir, pairdom_len[ij], 1.0, &(X2[0][0]), nso,
263
 
                  &(V[ij][0][0]), pairdom_len[ij], 0.0, &(RIjAb->matrix[0][ij][0]), nvir);
 
282
          C_DGEMM('n', 'n', nvir, local.pairdom_len[ij], local.pairdom_len[ij], 1.0, 
 
283
                  &(local.V[ij][0][0]), local.pairdom_len[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
 
284
          C_DGEMM('n', 't', nvir, nvir, local.pairdom_len[ij], 1.0, &(X2[0][0]), nso,
 
285
                  &(local.V[ij][0][0]), local.pairdom_len[ij], 0.0, &(RIjAb->matrix[0][ij][0]), nvir);
264
286
        }
265
287
        else  /* This must be a neglected weak pair; force it to zero */
266
288
          memset((void *) RIjAb->matrix[0][ij], 0, nvir*nvir*sizeof(double));
274
296
 
275
297
    dpd_buf4_mat_irrep_wrt(RIjAb, 0);
276
298
    dpd_buf4_mat_irrep_close(RIjAb, 0);
 
299
 
 
300
    /* Free Local Memory */
 
301
    for(i=0; i < nocc*nocc; i++) {
 
302
      free_block(local.W[i]);
 
303
      free_block(local.V[i]);
 
304
      free(local.eps_vir[i]);
 
305
    }
 
306
    free(local.W);
 
307
    free(local.V);
 
308
    free(local.eps_vir);
 
309
 
 
310
    free(local.eps_occ);
 
311
    free(local.pairdom_len);
 
312
    free(local.pairdom_nrlen);
 
313
    free(local.weak_pairs);
 
314
 
277
315
  }
278
316
  else {
279
317