3
\brief Enter brief description of file here
7
#include <libciomr/libciomr.h>
9
#include <libiwl/iwl.h>
14
#include "yoshimine.h"
16
namespace psi { namespace transqt {
19
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
21
/*--------------------------------------------------------
22
This routine transforms ERIs and (|r12|)-type integrals
23
from AO basis into (ip|jq)-type MO integrals (p and q
24
are general MO indices).
26
*******************************************************
27
NOTE: Full transfromation is done here! No frozen core
28
Nevermind the variables names
29
*******************************************************
30
--------------------------------------------------------*/
34
static void make_arrays(double ****Cdocc,
35
int **active_docc, int **active_virt,
36
int **focact, int **locact, int **fvract, int **lvract,
37
int **occ, int **vir, int **ioff3, int nmo);
39
void transform_two_mp2r12a_gr(void)
42
int psym, qsym, rsym, ssym,pqsym;
43
int pfirst, qfirst, rfirst, sfirst;
44
int plast, qlast, rlast, slast;
46
int ksym, lsym, klsym;
47
int kfirst, lfirst, ifirst, jfirst;
48
int klast, llast, ilast, jlast;
54
int *sopi, *clsdpi,*orbspi;
55
int *first_so, *last_so, *first, *last;
56
int *active_virt, *active_docc;
57
int *focact, *locact, *fvract, *lvract;
61
long int maxcor, maxcord;
63
double ***Cdocc, ***C;
69
struct yoshimine YBuffP, YBuffJ;
70
int max_buckets, print_lvl, first_tmp_file;
71
int presort_file, keep_presort, jfile, keep_half_tf, mfile;
76
fprintf(outfile,"\n\tThis system seems to be open-shell.\n");
77
fprintf(outfile,"MP2R12/A restricted transformations are available only ");
78
fprintf(outfile,"for closed-shell molecules.\n");
80
exit(PSI_RETURN_FAILURE);
83
ntri = moinfo.noeints;
86
nirreps = moinfo.nirreps;
88
orbspi = moinfo.orbspi;
89
clsdpi = moinfo.clsdpi;
90
first_so = moinfo.first_so;
91
last_so = moinfo.last_so;
94
max_buckets = params.max_buckets;
95
print_lvl = params.print_lvl;
96
first_tmp_file = params.first_tmp_file;
97
presort_file = params.presort_file;
98
keep_presort = params.keep_presort;
100
keep_half_tf = params.keep_half_tf;
101
mfile = params.mfile;
102
tolerance = params.tolerance;
103
print_integrals = params.print_te_ints;
104
maxcor = params.maxcor;
105
maxcord = params.maxcord;
109
/** Pre-sort the two-electron Integrals **/
110
if (params.print_lvl) {
111
fprintf(outfile, "\n\tPre-sorting two-electron ints...\n\n");
115
yosh_init(&YBuffP, ntri, ntri, maxcor, maxcord, max_buckets,
116
first_tmp_file, tolerance, outfile);
119
fprintf(outfile, "\tPresort");
120
yosh_print(&YBuffP, outfile);
121
fprintf(outfile, "\n");
125
yosh_init_buckets(&YBuffP);
127
yosh_rdtwo(&YBuffP, params.src_tei_file, params.delete_src_tei, sopi, nirreps, ioff, 0,
128
0, moinfo.fzc_density,
129
moinfo.fzc_operator, 1, (print_lvl > 5), outfile);
131
yosh_close_buckets(&YBuffP, 0);
133
yosh_sort(&YBuffP, presort_file, 0, ioff, NULL, nso, ntri,
134
0, 1, 0, 0, 1, (print_lvl > 5), outfile);
136
yosh_done(&YBuffP); /* Pre-transform complete */
138
/* no frozen core here */
142
for(h=0; h < nirreps; h++) {
144
nvirt += orbspi[h] - clsdpi[h];
148
"\tEntering MP2R12/A two-electron integral transformation...\n\n");
150
iwl_buf_init(&PBuff, presort_file, tolerance, 1, 1);
151
yosh_init(&YBuffJ, ndocc*nmo, ntri, maxcor, maxcord, max_buckets,
152
first_tmp_file, tolerance, outfile);
154
yosh_init_buckets(&YBuffJ);
156
fprintf(outfile, "\tHalf-transform");
157
yosh_print(&YBuffJ, outfile);
158
fprintf(outfile, "\n");
161
A = block_matrix(nso,nso);
162
B = block_matrix(nso,nso);
163
P_block = init_array(ntri);
165
make_arrays(&Cdocc, &active_docc, &active_virt,
166
&focact, &locact, &fvract, &lvract, &occ, &vir, &ioff3, nmo);
168
for(psym=0; psym < nirreps; psym++) {
169
pfirst = first_so[psym];
170
plast = last_so[psym];
171
for(p=pfirst; p <= plast; p++) {
172
for(qsym=0; qsym <= psym; qsym++) {
173
qfirst = first_so[qsym];
174
qlast = last_so[qsym];
176
for(q=qfirst; (q<=qlast) && (q <= p); q++) {
179
zero_arr(P_block,ntri);
180
iwl_buf_rd(&PBuff, pq, P_block, ioff, ioff, 0, 0, outfile);
182
for(rsym=0; rsym < nirreps; rsym++) {
183
rfirst = first_so[rsym];
184
rlast = last_so[rsym];
185
kfirst = focact[rsym];
186
klast = locact[rsym];
188
sfirst = first_so[ssym];
189
slast = last_so[ssym];
190
lfirst = first[ssym];
192
for(r=rfirst,R=0; r <= rlast; r++,R++) {
193
for(s=sfirst,S=0; s <= slast; s++,S++) {
195
A[R][S] = P_block[rs];
198
mmult(A,0,C[ssym],0,B,0,sopi[rsym],sopi[ssym],
200
mmult(Cdocc[rsym],1,B,0,A,0,active_docc[rsym],
201
sopi[rsym],orbspi[ssym],0);
203
yosh_wrt_arr_mp2r12a(&YBuffJ, p, q, pq, pqsym, A,
204
rsym, focact, locact, first, last,
206
(print_lvl >4), outfile);
213
fprintf(outfile, "\tSorting half-transformed integrals...\n");
215
iwl_buf_close(&PBuff, keep_presort);
217
yosh_close_buckets(&YBuffJ,0);
218
yosh_sort(&YBuffJ, jfile, 0, ioff3, ioff, nso, ntri, 0, 1, 1, nmo,
219
0, (print_lvl > 5), outfile);
222
fprintf(outfile, "\tFinished half-transform...\n");
223
fprintf(outfile, "\tWorking on second half...\n");
226
iwl_buf_init(&JBuff, jfile, tolerance, 1, 1);
227
iwl_buf_init(&MBuff, mfile, tolerance, 0, 0);
231
for(ksym=0; ksym < nirreps; ksym++) {
232
kfirst = (focact[ksym] < 0) ? -1 : occ[focact[ksym]];
233
klast = (locact[ksym] < 0) ? -2 : occ[locact[ksym]];
234
for(k=kfirst; k <= klast; k++) {
235
for(lsym=0; lsym < nirreps; lsym++) {
236
lfirst = first[lsym];
239
for(l=lfirst; l <= llast; l++) {
242
zero_arr(J_block, ntri);
243
iwl_buf_rd(&JBuff, kl, J_block, ioff3, ioff, 1, 0, outfile);
245
for(psym=0; psym < nirreps; psym++) {
246
pfirst = first_so[psym];
247
plast = last_so[psym];
248
ifirst = focact[psym];
249
ilast = locact[psym];
251
qfirst = first_so[qsym];
252
qlast = last_so[qsym];
253
jfirst = first[qsym];
255
for(p=pfirst,P=0; p <= plast; p++,P++) {
256
for(q=qfirst,Q=0; q <= qlast; q++,Q++) {
258
A[P][Q] = J_block[pq];
261
mmult(A,0,C[qsym],0,B,0,sopi[psym],sopi[qsym],
263
mmult(Cdocc[psym],1,B,0,A,0,active_docc[psym],
264
sopi[psym],orbspi[qsym],0);
265
iwl_buf_wrt_mp2r12a(&MBuff, k, l, kl, klsym, A, psym,
266
focact, locact, first, last,
267
occ, 1, ioff3,print_integrals,outfile);
277
iwl_buf_close(&JBuff, keep_half_tf);
278
/* Need to flush last buffer, else it's not written to disk */
279
iwl_buf_flush(&MBuff, 1);
280
iwl_buf_close(&MBuff, 1);
282
fprintf(outfile, "\n\tTransformation finished.\n");
284
fprintf(outfile, "\tTwo-electron integrals written to file%d.\n",mfile);
290
void make_arrays(double ****Cdocc,
291
int **active_docc, int **active_virt,
292
int **focact, int **locact, int **fvract, int **lvract,
293
int **occ, int **vir, int **ioff3, int nmo)
295
int h, i, offset, count;
296
int first_offset, last_offset;
299
/* active_docc[] and active_uocc[] are the number of active doubly-occupied
300
and unoccupied orbitals, respectively, in each irrep. */
301
*active_docc = init_int_array(moinfo.nirreps);
302
*active_virt = init_int_array(moinfo.nirreps);
303
for(h=0; h < moinfo.nirreps; h++) {
304
(*active_docc)[h] = moinfo.clsdpi[h];
305
(*active_virt)[h] = moinfo.virtpi[h];
308
/* focact[] and locact[] supply the first and last orbital indices (Pitzer
309
ordering) for the occupied orbitals in each irrep. */
310
*focact = init_int_array(moinfo.nirreps);
311
*locact = init_int_array(moinfo.nirreps);
312
for(h=0; h < moinfo.nirreps; h++) {
317
last_offset = moinfo.clsdpi[0] - 1;
318
(*focact)[0] = first_offset;
319
(*locact)[0] = last_offset;
320
for(h=1; h < moinfo.nirreps; h++) {
321
first_offset += moinfo.orbspi[h-1];
322
last_offset += moinfo.virtpi[h-1]+moinfo.orbspi[h]-moinfo.virtpi[h];
323
if((*active_docc)[h]) {
324
(*focact)[h] = first_offset;
325
(*locact)[h] = last_offset;
329
/* frvact[] and lrvact[] supply the first and last orbital indices (Pitzer
330
ordering) for the virtual orbitals in each irrep. */
331
*fvract = init_int_array(moinfo.nirreps);
332
*lvract = init_int_array(moinfo.nirreps);
333
for(h=0; h < moinfo.nirreps; h++) {
337
first_offset = moinfo.clsdpi[0];
338
last_offset = moinfo.orbspi[0] - 1;
339
(*fvract)[0] = first_offset;
340
(*lvract)[0] = last_offset;
341
for(h=1; h < moinfo.nirreps; h++) {
342
first_offset += moinfo.virtpi[h-1] + moinfo.clsdpi[h];
343
last_offset += moinfo.orbspi[h];
344
if((*active_virt)[h]) {
345
(*fvract)[h] = first_offset;
346
(*lvract)[h] = last_offset;
350
/* Construct occupied and virtual Pitzer -> QTS ordering arrays for
351
occupied (occ[]) and virtual (vir[]) orbitals */
352
*occ = init_int_array(moinfo.nmo);
353
*vir = init_int_array(moinfo.nmo);
354
for(i=0; i< moinfo.nmo; i++) {
361
for(h=0; h < moinfo.nirreps; h++) {
363
offset += moinfo.orbspi[h-1];
364
for(i=offset; i < (offset+moinfo.clsdpi[h]); i++) {
371
for(h=0; h < moinfo.nirreps; h++) {
373
offset += moinfo.orbspi[h-1];
374
for(i=offset+moinfo.clsdpi[h]; i < (offset+moinfo.orbspi[h]); i++) {
379
/* Generate ioff3 array. This array gives the row offset for an
380
ndocc x nmo matrix */
381
*ioff3 = init_int_array(MAXIOFF3);
382
for(i=0; i < MAXIOFF3; i++) {
386
/* Organize MOs for docc set only */
387
if(params.print_mos) {
388
fprintf(outfile,"\n\tSCF Eigenvectors (Occupied Set):\n");
389
fprintf(outfile, "\t--------------------------------\n");
391
*Cdocc = (double ***) malloc(moinfo.nirreps * sizeof(double **));
392
for(h=0; h < moinfo.nirreps; h++) {
393
if((*active_docc)[h]) {
394
(*Cdocc)[h] = block_matrix(moinfo.sopi[h],(*active_docc)[h]);
396
for(p=moinfo.first_so[h]; p <= moinfo.last_so[h]; p++) {
398
for(q=(*focact)[h]; q <= (*locact)[h]; q++) {
400
(*Cdocc)[h][row][col] = moinfo.scf_vector[p][q];
403
if(params.print_mos) {
404
fprintf(outfile,"\n\tDoubly Occupied Orbitals for Irrep %s\n",
406
print_mat((*Cdocc)[h],moinfo.sopi[h],(*active_docc)[h],
413
}} // end namespace psi::transqt
added
removed
src/bin/transqt/transform_two_mp2r12a_gr.cc
