3
#include <libciomr/libciomr.h>
5
#include <libiwl/iwl.h>
10
#include "yoshimine.h"
13
#define INDEX(i,j) ((i>j) ? (ioff[(i)]+(j)) : (ioff[(j)]+(i)))
15
/*--------------------------------------------------------
16
This routine transforms (|r12,T1|)-type integrals
17
from AO basis into (ip|[r12,T2]|jq)-type MO integrals
18
(p and q are general MO indices).
20
*******************************************************
21
NOTE: Full transfromation is done here! No frozen core
22
Nevermind the variables names
23
*******************************************************
24
--------------------------------------------------------*/
27
static void make_arrays(double ****Cdocc,
28
int **active_docc, int **active_virt,
29
int **focact, int **locact, int **fvract, int **lvract,
30
int **occ, int **vir, int **ioff3, int nmo);
32
void transform_two_mp2r12a_t(void)
35
int psym, qsym, rsym, ssym,pqsym;
36
int pfirst, qfirst, rfirst, sfirst;
37
int plast, qlast, rlast, slast;
39
int ksym, lsym, klsym;
40
int kfirst, lfirst, ifirst, jfirst;
41
int klast, llast, ilast, jlast;
47
int *sopi, *clsdpi,*orbspi;
48
int *first_so, *last_so, *first, *last;
49
int *active_virt, *active_docc;
50
int *focact, *locact, *fvract, *lvract;
54
long int maxcor, maxcord;
56
double ***Cdocc, ***C;
62
struct yoshimine YBuffP, YBuffJ;
63
int max_buckets, print_lvl, first_tmp_file;
64
int presort_file, keep_presort, jfile, keep_half_tf, mfile;
69
fprintf(outfile,"\n\tThis system seems to be open-shell.\n");
70
fprintf(outfile,"MP2R12/A restricted transformations are available only ");
71
fprintf(outfile,"for closed-shell molecules.\n");
73
exit(PSI_RETURN_FAILURE);
76
ntri = moinfo.noeints;
79
nirreps = moinfo.nirreps;
81
orbspi = moinfo.orbspi;
82
clsdpi = moinfo.clsdpi;
83
first_so = moinfo.first_so;
84
last_so = moinfo.last_so;
87
max_buckets = params.max_buckets;
88
print_lvl = params.print_lvl;
89
first_tmp_file = params.first_tmp_file;
90
presort_file = params.presort_file;
91
keep_presort = params.keep_presort;
93
keep_half_tf = params.keep_half_tf;
95
tolerance = params.tolerance;
96
print_integrals = params.print_te_ints;
97
maxcor = params.maxcor;
98
maxcord = params.maxcord;
102
/** Pre-sort the two-electron Integrals **/
103
if (params.print_lvl) {
104
fprintf(outfile, "\n\tPre-sorting two-electron ints...\n\n");
108
yosh_init(&YBuffP, ntri, ntri, maxcor, maxcord, max_buckets,
109
first_tmp_file, tolerance, outfile);
112
fprintf(outfile, "\tPresort");
113
yosh_print(&YBuffP, outfile);
114
fprintf(outfile, "\n");
118
yosh_init_buckets(&YBuffP);
120
/*------------------------------------------------------------------------
121
Notice that the third to last argument of this call to yosh_rdtwo is 0.
122
It is because [r12,T1] integrals are not symmetric with respect to
123
a permutation of bra and ket.
124
------------------------------------------------------------------------*/
125
yosh_rdtwo(&YBuffP, params.src_tei_file, params.delete_src_tei, orbspi, nirreps, ioff, 0,
126
0, moinfo.fzc_density,
127
moinfo.fzc_operator, 0, (print_lvl > 5), outfile);
129
yosh_close_buckets(&YBuffP, 0);
131
yosh_sort(&YBuffP, presort_file, 0, ioff, NULL, nso, ntri,
132
0, 1, 0, 0, 1, (print_lvl > 5), outfile);
134
yosh_done(&YBuffP); /* Pre-transform complete */
136
/* no frozen core here */
140
for(h=0; h < nirreps; h++) {
142
nvirt += orbspi[h] - clsdpi[h];
146
"\tEntering MP2R12/A two-electron integral transformation...\n\n");
148
iwl_buf_init(&PBuff, presort_file, tolerance, 1, 1);
149
yosh_init(&YBuffJ, ndocc*nmo, ntri, maxcor, maxcord, max_buckets,
150
first_tmp_file, tolerance, outfile);
152
yosh_init_buckets(&YBuffJ);
154
fprintf(outfile, "\tHalf-transform");
155
yosh_print(&YBuffJ, outfile);
156
fprintf(outfile, "\n");
159
A = block_matrix(nso,nso);
160
B = block_matrix(nso,nso);
161
P_block = init_array(ntri);
163
make_arrays(&Cdocc, &active_docc, &active_virt,
164
&focact, &locact, &fvract, &lvract, &occ, &vir, &ioff3, nmo);
166
for(psym=0; psym < nirreps; psym++) {
167
pfirst = first_so[psym];
168
plast = last_so[psym];
169
for(p=pfirst; p <= plast; p++) {
170
for(qsym=0; qsym <= psym; qsym++) {
171
qfirst = first_so[qsym];
172
qlast = last_so[qsym];
174
for(q=qfirst; (q<=qlast) && (q <= p); q++) {
177
zero_arr(P_block,ntri);
178
iwl_buf_rd(&PBuff, pq, P_block, ioff, ioff, 0, 0, outfile);
180
for(rsym=0; rsym < nirreps; rsym++) {
181
rfirst = first_so[rsym];
182
rlast = last_so[rsym];
183
kfirst = focact[rsym];
184
klast = locact[rsym];
186
sfirst = first_so[ssym];
187
slast = last_so[ssym];
188
lfirst = first[ssym];
190
for(r=rfirst,R=0; r <= rlast; r++,R++) {
191
for(s=sfirst,S=0; s <= slast; s++,S++) {
193
/*------------------------------------
194
(pq|[r12,T1]|rs) = (pq|[r12,T1]|sr)
195
------------------------------------*/
197
A[R][S] = P_block[rs];
200
mmult(A,0,C[ssym],0,B,0,sopi[rsym],sopi[ssym],
202
mmult(Cdocc[rsym],1,B,0,A,0,active_docc[rsym],
203
sopi[rsym],orbspi[ssym],0);
205
yosh_wrt_arr_mp2r12a(&YBuffJ, p, q, pq, pqsym, A,
206
rsym, focact, locact, first, last,
208
(print_lvl >4), outfile);
215
/*---------------------------------------
216
Now the integrals are (kl|[r12,T2]|pq)
217
---------------------------------------*/
219
fprintf(outfile, "\tSorting half-transformed integrals...\n");
221
iwl_buf_close(&PBuff, keep_presort);
223
yosh_close_buckets(&YBuffJ,0);
224
yosh_sort(&YBuffJ, jfile, 0, ioff3, ioff, nso, ntri, 0, 1, 1, nmo,
225
0, (print_lvl > 5), outfile);
228
fprintf(outfile, "\tFinished half-transform...\n");
229
fprintf(outfile, "\tWorking on second half...\n");
232
iwl_buf_init(&JBuff, jfile, tolerance, 1, 1);
233
iwl_buf_init(&MBuff, mfile, tolerance, 0, 0);
237
for(ksym=0; ksym < nirreps; ksym++) {
238
kfirst = (focact[ksym] < 0) ? -1 : occ[focact[ksym]];
239
klast = (locact[ksym] < 0) ? -2 : occ[locact[ksym]];
240
for(k=kfirst; k <= klast; k++) {
241
for(lsym=0; lsym < nirreps; lsym++) {
242
lfirst = first[lsym];
245
for(l=lfirst; l <= llast; l++) {
248
zero_arr(J_block, ntri);
249
iwl_buf_rd(&JBuff, kl, J_block, ioff3, ioff, 1, 0, outfile);
251
for(psym=0; psym < nirreps; psym++) {
252
pfirst = first_so[psym];
253
plast = last_so[psym];
254
ifirst = focact[psym];
255
ilast = locact[psym];
257
qfirst = first_so[qsym];
258
qlast = last_so[qsym];
259
jfirst = first[qsym];
261
for(p=pfirst,P=0; p <= plast; p++,P++) {
262
for(q=qfirst,Q=0; q <= qlast; q++,Q++) {
263
/*-------------------------------------
264
(kl|[r12,T2]|pq) = -(kl|[r12,T2]|qp)
265
-------------------------------------*/
267
A[P][Q] = J_block[ioff[p]+q];
269
A[P][Q] = (-1.0)*J_block[ioff[q]+p];
272
mmult(A,0,C[qsym],0,B,0,sopi[psym],sopi[qsym],
274
mmult(Cdocc[psym],1,B,0,A,0,active_docc[psym],
275
sopi[psym],orbspi[qsym],0);
276
iwl_buf_wrt_mp2r12a(&MBuff, k, l, kl, klsym, A, psym,
277
focact, locact, first, last,
278
occ, 0, ioff3,print_integrals,outfile);
285
/*---------------------------------------------------------------
286
NOTE: Written integrals are of [r12,T2] operator, not [r12,T1]
287
It doesn't really matter but is good to remember!
288
---------------------------------------------------------------*/
293
iwl_buf_close(&JBuff, keep_half_tf);
294
/* Need to flush last buffer, else it's not written to disk */
295
iwl_buf_flush(&MBuff, 1);
296
iwl_buf_close(&MBuff, 1);
298
fprintf(outfile, "\n\tTransformation finished.\n");
300
fprintf(outfile, "\tTwo-electron integrals written to file%d.\n",mfile);
306
void make_arrays(double ****Cdocc,
307
int **active_docc, int **active_virt,
308
int **focact, int **locact, int **fvract, int **lvract,
309
int **occ, int **vir, int **ioff3, int nmo)
311
int h, i, offset, count;
312
int first_offset, last_offset;
315
/* active_docc[] and active_uocc[] are the number of active doubly-occupied
316
and unoccupied orbitals, respectively, in each irrep. */
317
*active_docc = init_int_array(moinfo.nirreps);
318
*active_virt = init_int_array(moinfo.nirreps);
319
for(h=0; h < moinfo.nirreps; h++) {
320
(*active_docc)[h] = moinfo.clsdpi[h];
321
(*active_virt)[h] = moinfo.virtpi[h];
324
/* focact[] and locact[] supply the first and last orbital indices (Pitzer
325
ordering) for the occupied orbitals in each irrep. */
326
*focact = init_int_array(moinfo.nirreps);
327
*locact = init_int_array(moinfo.nirreps);
328
for(h=0; h < moinfo.nirreps; h++) {
333
last_offset = moinfo.clsdpi[0] - 1;
334
(*focact)[0] = first_offset;
335
(*locact)[0] = last_offset;
336
for(h=1; h < moinfo.nirreps; h++) {
337
first_offset += moinfo.orbspi[h-1];
338
last_offset += moinfo.virtpi[h-1]+moinfo.orbspi[h]-moinfo.virtpi[h];
339
if((*active_docc)[h]) {
340
(*focact)[h] = first_offset;
341
(*locact)[h] = last_offset;
345
/* frvact[] and lrvact[] supply the first and last orbital indices (Pitzer
346
ordering) for the virtual orbitals in each irrep. */
347
*fvract = init_int_array(moinfo.nirreps);
348
*lvract = init_int_array(moinfo.nirreps);
349
for(h=0; h < moinfo.nirreps; h++) {
353
first_offset = moinfo.clsdpi[0];
354
last_offset = moinfo.orbspi[0] - 1;
355
(*fvract)[0] = first_offset;
356
(*lvract)[0] = last_offset;
357
for(h=1; h < moinfo.nirreps; h++) {
358
first_offset += moinfo.virtpi[h-1] + moinfo.clsdpi[h];
359
last_offset += moinfo.orbspi[h];
360
if((*active_virt)[h]) {
361
(*fvract)[h] = first_offset;
362
(*lvract)[h] = last_offset;
366
/* Construct occupied and virtual Pitzer -> QTS ordering arrays for
367
occupied (occ[]) and virtual (vir[]) orbitals */
368
*occ = init_int_array(moinfo.nmo);
369
*vir = init_int_array(moinfo.nmo);
370
for(i=0; i< moinfo.nmo; i++) {
377
for(h=0; h < moinfo.nirreps; h++) {
379
offset += moinfo.orbspi[h-1];
380
for(i=offset; i < (offset+moinfo.clsdpi[h]); i++) {
387
for(h=0; h < moinfo.nirreps; h++) {
389
offset += moinfo.orbspi[h-1];
390
for(i=offset+moinfo.clsdpi[h]; i < (offset+moinfo.orbspi[h]); i++) {
395
/* Generate ioff3 array. This array gives the row offset for an
396
ndocc x nmo matrix */
397
*ioff3 = init_int_array(MAXIOFF3);
398
for(i=0; i < MAXIOFF3; i++) {
402
/* Organize MOs for docc set only */
403
if(params.print_mos) {
404
fprintf(outfile,"\n\tSCF Eigenvectors (Occupied Set):\n");
405
fprintf(outfile, "\t--------------------------------\n");
407
*Cdocc = (double ***) malloc(moinfo.nirreps * sizeof(double **));
408
for(h=0; h < moinfo.nirreps; h++) {
409
if((*active_docc)[h]) {
410
(*Cdocc)[h] = block_matrix(moinfo.sopi[h],(*active_docc)[h]);
412
for(p=moinfo.first_so[h]; p <= moinfo.last_so[h]; p++) {
414
for(q=(*focact)[h]; q <= (*locact)[h]; q++) {
416
(*Cdocc)[h][row][col] = moinfo.scf_vector[p][q];
419
if(params.print_mos) {
420
fprintf(outfile,"\n\tDoubly Occupied Orbitals for Irrep %s\n",
422
print_mat((*Cdocc)[h],moinfo.sopi[h],(*active_docc)[h],
added
removed
src/bin/transqt/transform_two_mp2r12a_t.c
