← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/bin/cints/Tools/ccinfo.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/cints/Tools/ccinfo.c
 
1
#include<stdio.h>
 
2
#include<stdlib.h>
 
3
#include<math.h>
 
4
#include<libciomr/libciomr.h>
 
5
#include<libpsio/psio.h>
 
6
#include<libdpd/dpd.h>
 
7
#include<libchkpt/chkpt.h>
 
8
#include<libint/libint.h>
 
9
#include<libqt/qt.h>
 
10
#include <psifiles.h>
 
11
#include <ccfiles.h>
 
12
 
 
13
#include"moinfo.h"
 
14
#include"moinfo_corr.h"
 
15
 
 
16
#include"defines.h"
 
17
#define EXTERN
 
18
#include"global.h"
 
19
#include"small_fns.h"
 
20
 
 
21
static int *cachefiles, **cachelist;
 
22
 
 
23
void init_ccinfo()
 
24
{
 
25
  int h, ij, ab, row, col, row_offset, col_offset;
 
26
  int nactive, nocc, nvirt, nao;
 
27
  int *occpi, *occ_sym;
 
28
  int *virtpi, *vir_sym;
 
29
  double **A;
 
30
  double **T2_MO, **T2_AO, **T2;
 
31
  dpdbuf4 tau;
 
32
 
 
33
  nao = BasisSet.num_ao;
 
34
 
 
35
  /* grab some basic data from CC_INFO */
 
36
  psio_open(CC_INFO, PSIO_OPEN_OLD);
 
37
  psio_read_entry(CC_INFO, "No. of Active Orbitals", (char *) &(nactive),
 
38
                  sizeof(int));
 
39
  occpi = init_int_array(Symmetry.nirreps);
 
40
  virtpi = init_int_array(Symmetry.nirreps);
 
41
  occ_sym = init_int_array(nactive);
 
42
  vir_sym = init_int_array(nactive);
 
43
  psio_read_entry(CC_INFO, "Active Occ Orbs Per Irrep",
 
44
                  (char *) occpi, sizeof(int)*Symmetry.nirreps);
 
45
  psio_read_entry(CC_INFO, "Active Virt Orbs Per Irrep",
 
46
                  (char *) virtpi, sizeof(int)*Symmetry.nirreps);
 
47
  psio_read_entry(CC_INFO, "Active Occ Orb Symmetry",
 
48
                  (char *) occ_sym, sizeof(int)*Symmetry.nirreps);
 
49
  psio_read_entry(CC_INFO, "Active Virt Orb Symmetry",
 
50
                  (char *) vir_sym, sizeof(int)*Symmetry.nirreps);
 
51
  psio_close(CC_INFO, 1);
 
52
 
 
53
  nocc = 0; nvirt = 0;
 
54
  for(h=0; h < Symmetry.nirreps; h++) {
 
55
    nocc += occpi[h];  nvirt += virtpi[h];
 
56
  }
 
57
 
 
58
  cachefiles = init_int_array(PSIO_MAXUNIT);
 
59
  /* assuming no caching for DPD here */
 
60
  cachelist = init_int_matrix(12,12);
 
61
 
 
62
  init_moinfo();
 
63
  init_moinfo_corr();
 
64
 
 
65
  /* open the T-amplitude file for r/w */
 
66
  psio_open(CC_TAMPS, PSIO_OPEN_OLD);
 
67
 
 
68
  /*--- Initialize DPD library ---*/
 
69
  dpd_init(1, Symmetry.nirreps, UserOptions.memory, 0, cachefiles, cachelist, NULL, 
 
70
           2, occpi, occ_sym, virtpi, vir_sym);
 
71
 
 
72
  /* Grab the MO-basis T2's provided by ccenergy */
 
73
  T2_MO = block_matrix(nocc*nocc,nvirt*nvirt);
 
74
  dpd_buf4_init(&tau, CC_TAMPS, 0, 0, 5, 0, 5, 0, "tauIjAb");
 
75
  row_offset = col_offset = 0;
 
76
  for(h=0,row=0,col=0; h < Symmetry.nirreps; h++) {
 
77
    dpd_buf4_mat_irrep_init(&tau, h);
 
78
    dpd_buf4_mat_irrep_rd(&tau, h);
 
79
 
 
80
    if(h) { 
 
81
      row_offset += tau.params->rowtot[h-1];  
 
82
      col_offset += tau.params->coltot[h-1];
 
83
    }
 
84
 
 
85
    for(row=0; row < tau.params->rowtot[h]; row++) {
 
86
      for(col=0; col < tau.params->coltot[h]; col++) {
 
87
        T2_MO[row_offset+row][col_offset+col] = tau.matrix[h][row][col];
 
88
      }
 
89
    }
 
90
 
 
91
  dpd_buf4_mat_irrep_close(&tau, h);
 
92
  }
 
93
  dpd_buf4_close(&tau);
 
94
 
 
95
  /* half-transform the T2's */
 
96
  A = block_matrix(nvirt, nao);
 
97
  T2_AO = block_matrix(nocc*nocc, nao*nao);
 
98
  for(ij=0; ij < nocc*nocc; ij++) {
 
99
 
 
100
 
 
101
    C_DGEMM('n','n',nvirt,nao,nvirt,1.0,&(T2_MO[ij][0]),nvirt,
 
102
            &(MOInfo.scf_evec_uocc[0][0][0]),nao,0.0,A[0],nao);
 
103
 
 
104
    C_DGEMM('t','n',nao,nao,nvirt,1.0,&(MOInfo.scf_evec_uocc[0][0][0]),nao,
 
105
            A[0],nao,0.0,&(T2_AO[ij][0]),nao);
 
106
  }
 
107
  free_block(A);
 
108
  free_block(T2_MO);
 
109
 
 
110
  /* transpose T2 to AO*AO x MO*MO */
 
111
  T2 = block_matrix(nao*nao,nocc*nocc);
 
112
  for(ij=0; ij < nocc*nocc; ij++)
 
113
    for(ab=0; ab < nao*nao; ab++) 
 
114
      T2[ab][ij] = T2_AO[ij][ab];
 
115
  free_block(T2_AO);
 
116
 
 
117
  CCInfo.T2_s = T2;
 
118
  CCInfo.T2_t = block_matrix(nao*nao,nocc*nocc);
 
119
  CCInfo.nocc = nocc;
 
120
  CCInfo.nvirt = nvirt;
 
121
 
 
122
  free(occpi);  free(occ_sym);
 
123
  free(virtpi); free(vir_sym);
 
124
      
 
125
  return;
 
126
}
 
127
 
 
128
 
 
129
void cleanup_ccinfo()
 
130
{
 
131
  int h, ij, ab, row, col, row_offset, col_offset;
 
132
  int nocc, nvirt, nao;
 
133
  double **T2, **A, **T2_MO;
 
134
  dpdbuf4 tau;
 
135
 
 
136
  nocc = CCInfo.nocc;
 
137
  nvirt = CCInfo.nvirt;
 
138
  nao = BasisSet.num_ao;
 
139
 
 
140
  free_block(CCInfo.T2_s);
 
141
 
 
142
  /* transpose target T2 back to MO*MO x AO*AO */
 
143
  T2 = block_matrix(nocc*nocc, nao*nao); 
 
144
  for(ij=0; ij < nocc*nocc; ij++) {
 
145
    for(ab=0; ab < nao*nao; ab++) {
 
146
       T2[ij][ab] = CCInfo.T2_t[ab][ij];
 
147
    }
 
148
  }
 
149
  free_block(CCInfo.T2_t);
 
150
 
 
151
  /* transform T2 back to MO basis */
 
152
  A = block_matrix(nao,nvirt);
 
153
  T2_MO = block_matrix(nocc*nocc,nvirt*nvirt);
 
154
  for(ij=0; ij < nocc*nocc; ij++) {
 
155
    C_DGEMM('n','t',nao,nvirt,nao,1.0,T2[ij],nao,
 
156
            MOInfo.scf_evec_uocc[0][0],nao,0.0,A[0],nvirt);
 
157
    C_DGEMM('n','n',nvirt,nvirt,nao,1.0,MOInfo.scf_evec_uocc[0][0],nao,
 
158
            A[0],nvirt,0.0,T2_MO[ij],nvirt);
 
159
  }
 
160
  free_block(A);
 
161
  free_block(T2);
 
162
 
 
163
  /* Put T2 back into the DPD buffer */
 
164
  dpd_buf4_init(&tau, CC_TAMPS, 0, 0, 5, 0, 5, 0, "New tIjAb");
 
165
  row_offset = col_offset = 0;
 
166
  for(h=0; h < Symmetry.nirreps; h++) {
 
167
    dpd_buf4_mat_irrep_init(&tau, h);
 
168
    dpd_buf4_mat_irrep_rd(&tau, h);
 
169
 
 
170
    if(h) {
 
171
      row_offset += tau.params->rowtot[h-1];
 
172
      col_offset += tau.params->coltot[h-1];
 
173
    }
 
174
 
 
175
    for(row=0; row < tau.params->rowtot[h]; row++) {
 
176
      for(col=0; col < tau.params->coltot[h]; col++) {
 
177
         tau.matrix[h][row][col] += T2_MO[row+row_offset][col+col_offset];
 
178
      }
 
179
    }
 
180
    dpd_buf4_mat_irrep_wrt(&tau, h);
 
181
    dpd_buf4_mat_irrep_close(&tau, h);
 
182
  }
 
183
  dpd_buf4_close(&tau);
 
184
 
 
185
  free_block(T2_MO);
 
186
 
 
187
  free_int_matrix(cachelist,12);
 
188
  free(cachefiles);
 
189
  dpd_close(1);
 
190
 
 
191
  /* close the T-amplitude file */
 
192
  psio_close(CC_TAMPS, 1);
 
193
 
 
194
  return;
 
195
}
 
196
 
 
197