← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/ccenergy/local.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/ccenergy/local.cc
 
1
/*! \file
 
2
    \ingroup CCENERGY
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
 
 
6
#include <cstdio>
 
7
#include <cstdlib>
 
8
#include <cstring>
 
9
#include <cmath>
 
10
#include <libipv1/ip_lib.h>
 
11
#include <libciomr/libciomr.h>
 
12
#include <libpsio/psio.h>
 
13
#include <libiwl/iwl.h>
 
14
#include <libint/libint.h>
 
15
#include <libchkpt/chkpt.h>
 
16
#include <libqt/qt.h>
 
17
#include <libdpd/dpd.h>
 
18
#include <psifiles.h>
 
19
#include "Local.h"
 
20
#include "Params.h"
 
21
#include "MOInfo.h"
 
22
#define EXTERN
 
23
#include "globals.h"
 
24
 
 
25
namespace psi { namespace ccenergy {
 
26
 
 
27
/*! 
 
28
** local_init(): Set up parameters of local excitation domains.
 
29
**
 
30
** The orbital domains constructed here are based on those described
 
31
** in Broughton and Pulay, J. Comp. Chem. 14, 736-740 (1993).  The
 
32
** localization of the occupied orbitals is done elsewhere (see the
 
33
** program "local").  Pair domains are defined as the union of pairs
 
34
** of single occupied orbital domains.  "Weak pairs", which are
 
35
** defined as pair domains whose individual occupied orbital domains
 
36
** have no atoms in common, are identified (cf. int *weak_pairs).
 
37
** 
 
38
** TDC, Jan-June 2002
 
39
*/
 
40
 
 
41
void local_init(void)
 
42
{
 
43
  int i, k, ij, nocc;
 
44
 
 
45
  chkpt_init(PSIO_OPEN_OLD);
 
46
  local.natom = chkpt_rd_natom();
 
47
  chkpt_close();
 
48
 
 
49
  local.nso = moinfo.nso;
 
50
  local.nocc = moinfo.occpi[0]; /* active doubly occupied orbitals */
 
51
  local.nvir = moinfo.virtpi[0]; /* active virtual orbitals */
 
52
  nocc = local.nocc;
 
53
 
 
54
  local.weak_pair_energy = 0.0;
 
55
 
 
56
  local.weak_pairs = init_int_array(nocc*nocc);
 
57
  psio_read_entry(CC_INFO, "Local Weak Pairs", (char *) local.weak_pairs,
 
58
                  local.nocc*local.nocc*sizeof(int));
 
59
 
 
60
  fprintf(outfile, "\tLocalization parameters ready.\n\n");
 
61
  fflush(outfile);
 
62
}
 
63
 
 
64
void local_done(void)
 
65
{
 
66
  fprintf(outfile, "\tLocal parameters free.\n");
 
67
  fflush(outfile);
 
68
}
 
69
 
 
70
void local_filter_T1(dpdfile2 *T1)
 
71
{
 
72
  int i, a, ij, ii;
 
73
  int nocc, nvir;
 
74
  double *T1tilde, *T1bar;
 
75
  psio_address next;
 
76
 
 
77
  nocc = local.nocc;
 
78
  nvir = local.nvir;
 
79
 
 
80
/*   local.weak_pairs = init_int_array(nocc*nocc); */
 
81
  local.pairdom_len = init_int_array(nocc*nocc);
 
82
  local.pairdom_nrlen = init_int_array(nocc*nocc);
 
83
  local.eps_occ = init_array(nocc);
 
84
/*   psio_read_entry(CC_INFO, "Local Weak Pairs", (char *) local.weak_pairs, */
 
85
/*                local.nocc*local.nocc*sizeof(int)); */
 
86
  psio_read_entry(CC_INFO, "Local Pair Domain Length", (char *) local.pairdom_len,
 
87
                  nocc*nocc*sizeof(int));
 
88
  psio_read_entry(CC_INFO, "Local Pair Domain NR Length", (char *) local.pairdom_nrlen,
 
89
                  nocc*nocc*sizeof(int));
 
90
  psio_read_entry(CC_INFO, "Local Occupied Orbital Energies", (char *) local.eps_occ,
 
91
                  nocc*sizeof(double));
 
92
 
 
93
  local.W = (double ***) malloc(nocc * nocc * sizeof(double **));
 
94
  local.V = (double ***) malloc(nocc * nocc * sizeof(double **));
 
95
  local.eps_vir = (double **) malloc(nocc * nocc * sizeof(double *));
 
96
  next = PSIO_ZERO;
 
97
  for(ij=0; ij < nocc*nocc; ij++) {
 
98
      local.eps_vir[ij] = init_array(local.pairdom_nrlen[ij]);
 
99
      psio_read(CC_INFO, "Local Virtual Orbital Energies", (char *) local.eps_vir[ij],
 
100
                local.pairdom_nrlen[ij]*sizeof(double), next, &next);
 
101
    }
 
102
  next = PSIO_ZERO;
 
103
  for(ij=0; ij < nocc*nocc; ij++) {
 
104
      local.V[ij] = block_matrix(nvir,local.pairdom_len[ij]);
 
105
      psio_read(CC_INFO, "Local Residual Vector (V)", (char *) local.V[ij][0],
 
106
                nvir*local.pairdom_len[ij]*sizeof(double), next, &next);
 
107
    }
 
108
  next = PSIO_ZERO;
 
109
  for(ij=0; ij < nocc*nocc; ij++) {
 
110
      local.W[ij] = block_matrix(local.pairdom_len[ij],local.pairdom_nrlen[ij]);
 
111
      psio_read(CC_INFO, "Local Transformation Matrix (W)", (char *) local.W[ij][0],
 
112
                local.pairdom_len[ij]*local.pairdom_nrlen[ij]*sizeof(double), next, &next);
 
113
    }
 
114
 
 
115
  dpd_file2_mat_init(T1);
 
116
  dpd_file2_mat_rd(T1);
 
117
 
 
118
  for(i=0; i<nocc; i++) {
 
119
    ii = i*nocc + i;  /* diagonal element of pair matrices */
 
120
 
 
121
    if(!local.pairdom_len[ii]) {
 
122
      fprintf(outfile, "\n\tlocal_filter_T1: Pair ii = [%d] is zero-length, which makes no sense.\n",ii);
 
123
      exit(PSI_RETURN_FAILURE);
 
124
    }
 
125
 
 
126
    T1tilde = init_array(local.pairdom_len[ii]);
 
127
    T1bar = init_array(local.pairdom_nrlen[ii]);
 
128
 
 
129
    /* Transform the virtuals to the redundant projected virtual basis */
 
130
    C_DGEMV('t', nvir, local.pairdom_len[ii], 1.0, &(local.V[ii][0][0]), local.pairdom_len[ii], 
 
131
            &(T1->matrix[0][i][0]), 1, 0.0, &(T1tilde[0]), 1);
 
132
 
 
133
    /* Transform the virtuals to the non-redundant virtual basis */
 
134
    C_DGEMV('t', local.pairdom_len[ii], local.pairdom_nrlen[ii], 1.0, &(local.W[ii][0][0]), local.pairdom_nrlen[ii], 
 
135
            &(T1tilde[0]), 1, 0.0, &(T1bar[0]), 1);
 
136
 
 
137
    /* Apply the denominators */
 
138
    for(a=0; a < local.pairdom_nrlen[ii]; a++)
 
139
      T1bar[a] /= (local.eps_occ[i] - local.eps_vir[ii][a]);
 
140
 
 
141
    /* Transform the new T1's to the redundant projected virtual basis */
 
142
    C_DGEMV('n', local.pairdom_len[ii], local.pairdom_nrlen[ii], 1.0, &(local.W[ii][0][0]), local.pairdom_nrlen[ii],
 
143
            &(T1bar[0]), 1, 0.0, &(T1tilde[0]), 1);
 
144
 
 
145
 
 
146
    /* Transform the new T1's to the MO basis */
 
147
    C_DGEMV('n', nvir, local.pairdom_len[ii], 1.0, &(local.V[ii][0][0]), local.pairdom_len[ii], 
 
148
            &(T1tilde[0]), 1, 0.0, &(T1->matrix[0][i][0]), 1);
 
149
 
 
150
    free(T1bar);
 
151
    free(T1tilde);
 
152
 
 
153
  }
 
154
 
 
155
  dpd_file2_mat_wrt(T1);
 
156
  dpd_file2_mat_close(T1);
 
157
 
 
158
  for(ij=0; ij<nocc*nocc; ij++) {
 
159
      free_block(local.W[ij]);
 
160
      free_block(local.V[ij]);
 
161
      free(local.eps_vir[ij]);
 
162
    }
 
163
  free(local.W);
 
164
  free(local.V);
 
165
  free(local.eps_vir);
 
166
 
 
167
  free(local.eps_occ);
 
168
  free(local.pairdom_len);
 
169
  free(local.pairdom_nrlen);
 
170
/*   free(local.weak_pairs); */
 
171
}
 
172
 
 
173
void local_filter_T2(dpdbuf4 *T2)
 
174
{
 
175
  int ij, i, j, a, b;
 
176
  int nso, nocc, nvir;
 
177
  double **X1, **X2, **T2tilde, **T2bar;
 
178
  psio_address next;
 
179
 
 
180
  nso = local.nso;
 
181
  nocc = local.nocc;
 
182
  nvir = local.nvir;
 
183
 
 
184
/*   local.weak_pairs = init_int_array(nocc*nocc); */
 
185
  local.pairdom_len = init_int_array(nocc*nocc);
 
186
  local.pairdom_nrlen = init_int_array(nocc*nocc);
 
187
  local.eps_occ = init_array(nocc);
 
188
/*   psio_read_entry(CC_INFO, "Local Weak Pairs", (char *) local.weak_pairs, */
 
189
/*                local.nocc*local.nocc*sizeof(int)); */
 
190
  psio_read_entry(CC_INFO, "Local Pair Domain Length", (char *) local.pairdom_len,
 
191
                  nocc*nocc*sizeof(int));
 
192
  psio_read_entry(CC_INFO, "Local Pair Domain NR Length", (char *) local.pairdom_nrlen,
 
193
                  nocc*nocc*sizeof(int));
 
194
  psio_read_entry(CC_INFO, "Local Occupied Orbital Energies", (char *) local.eps_occ,
 
195
                  nocc*sizeof(double));
 
196
 
 
197
  local.W = (double ***) malloc(nocc * nocc * sizeof(double **));
 
198
  local.V = (double ***) malloc(nocc * nocc * sizeof(double **));
 
199
  local.eps_vir = (double **) malloc(nocc * nocc * sizeof(double *));
 
200
  next = PSIO_ZERO;
 
201
  for(ij=0; ij < nocc*nocc; ij++) {
 
202
      local.eps_vir[ij] = init_array(local.pairdom_nrlen[ij]);
 
203
      psio_read(CC_INFO, "Local Virtual Orbital Energies", (char *) local.eps_vir[ij],
 
204
                local.pairdom_nrlen[ij]*sizeof(double), next, &next);
 
205
    }
 
206
  next = PSIO_ZERO;
 
207
  for(ij=0; ij < nocc*nocc; ij++) {
 
208
      local.V[ij] = block_matrix(nvir,local.pairdom_len[ij]);
 
209
      psio_read(CC_INFO, "Local Residual Vector (V)", (char *) local.V[ij][0],
 
210
                nvir*local.pairdom_len[ij]*sizeof(double), next, &next);
 
211
    }
 
212
  next = PSIO_ZERO;
 
213
  for(ij=0; ij < nocc*nocc; ij++) {
 
214
      local.W[ij] = block_matrix(local.pairdom_len[ij],local.pairdom_nrlen[ij]);
 
215
      psio_read(CC_INFO, "Local Transformation Matrix (W)", (char *) local.W[ij][0],
 
216
                local.pairdom_len[ij]*local.pairdom_nrlen[ij]*sizeof(double), next, &next);
 
217
    }
 
218
 
 
219
  /* Grab the MO-basis T2's */
 
220
  dpd_buf4_mat_irrep_init(T2, 0);
 
221
  dpd_buf4_mat_irrep_rd(T2, 0);
 
222
 
 
223
  X1 = block_matrix(nso,nvir);
 
224
  X2 = block_matrix(nvir,nso);
 
225
  T2tilde = block_matrix(nso,nso);
 
226
  T2bar = block_matrix(nvir, nvir);
 
227
 
 
228
  for(i=0,ij=0; i<nocc; i++) {
 
229
    for(j=0; j<nocc; j++,ij++) {
 
230
 
 
231
      if(!local.weak_pairs[ij]) {
 
232
 
 
233
        /* Transform the virtuals to the redundant projected virtual basis */
 
234
        C_DGEMM('t', 'n', local.pairdom_len[ij], nvir, nvir, 1.0, &(local.V[ij][0][0]), local.pairdom_len[ij],
 
235
                &(T2->matrix[0][ij][0]), nvir, 0.0, &(X1[0][0]), nvir);
 
236
        C_DGEMM('n', 'n', local.pairdom_len[ij], local.pairdom_len[ij], nvir, 1.0, &(X1[0][0]), nvir,
 
237
                &(local.V[ij][0][0]), local.pairdom_len[ij], 0.0, &(T2tilde[0][0]), nso);
 
238
 
 
239
        /* Transform the virtuals to the non-redundant virtual basis */
 
240
        C_DGEMM('t', 'n', local.pairdom_nrlen[ij], local.pairdom_len[ij], local.pairdom_len[ij], 1.0,
 
241
                &(local.W[ij][0][0]), local.pairdom_nrlen[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
 
242
        C_DGEMM('n', 'n', local.pairdom_nrlen[ij], local.pairdom_nrlen[ij], local.pairdom_len[ij], 1.0,
 
243
                &(X2[0][0]), nso, &(local.W[ij][0][0]), local.pairdom_nrlen[ij], 0.0, &(T2bar[0][0]), nvir);
 
244
 
 
245
        /* Divide the new amplitudes by the denominators */
 
246
        for(a=0; a < local.pairdom_nrlen[ij]; a++)
 
247
          for(b=0; b < local.pairdom_nrlen[ij]; b++)
 
248
            T2bar[a][b] /= (local.eps_occ[i] + local.eps_occ[j]
 
249
                            - local.eps_vir[ij][a] - local.eps_vir[ij][b]);
 
250
 
 
251
        /* Transform the new T2's to the redundant virtual basis */
 
252
        C_DGEMM('n', 'n', local.pairdom_len[ij], local.pairdom_nrlen[ij], local.pairdom_nrlen[ij], 1.0,
 
253
                &(local.W[ij][0][0]), local.pairdom_nrlen[ij], &(T2bar[0][0]), nvir, 0.0, &(X1[0][0]), nvir);
 
254
        C_DGEMM('n','t', local.pairdom_len[ij], local.pairdom_len[ij], local.pairdom_nrlen[ij], 1.0,
 
255
                &(X1[0][0]), nvir, &(local.W[ij][0][0]), local.pairdom_nrlen[ij], 0.0, &(T2tilde[0][0]), nso);
 
256
 
 
257
        /* Transform the new T2's to the MO basis */
 
258
        C_DGEMM('n', 'n', nvir, local.pairdom_len[ij], local.pairdom_len[ij], 1.0,
 
259
                &(local.V[ij][0][0]), local.pairdom_len[ij], &(T2tilde[0][0]), nso, 0.0, &(X2[0][0]), nso);
 
260
        C_DGEMM('n', 't', nvir, nvir, local.pairdom_len[ij], 1.0, &(X2[0][0]), nso,
 
261
                &(local.V[ij][0][0]), local.pairdom_len[ij], 0.0, &(T2->matrix[0][ij][0]), nvir);
 
262
      }
 
263
      else  /* This must be a neglected weak pair; force it to zero */
 
264
        memset((void *) T2->matrix[0][ij], 0, nvir*nvir*sizeof(double));
 
265
 
 
266
    }
 
267
  }
 
268
 
 
269
  free_block(X1);
 
270
  free_block(X2);
 
271
  free_block(T2tilde);
 
272
  free_block(T2bar);
 
273
 
 
274
  /* Write the updated MO-basis T2's to disk */
 
275
  dpd_buf4_mat_irrep_wrt(T2, 0);
 
276
  dpd_buf4_mat_irrep_close(T2, 0);
 
277
 
 
278
  for(i=0; i<nocc*nocc; i++) {
 
279
      free_block(local.W[i]);
 
280
      free_block(local.V[i]);
 
281
      free(local.eps_vir[i]);
 
282
    }
 
283
  free(local.W);
 
284
  free(local.V);
 
285
  free(local.eps_vir);
 
286
 
 
287
  free(local.eps_occ);
 
288
  free(local.pairdom_len);
 
289
  free(local.pairdom_nrlen);
 
290
/*   free(local.weak_pairs); */
 
291
}
 
292
}} // namespace psi::ccenergy