4
** A program to transform one- and two-electron integrals from the
5
** symmetry-orbital basis to the molecular-orbital basis.
7
** This code replaces the original transqt code developed initially
8
** in 1995 by TDC, CDS, and JTF. This version is designed to take
9
** advantage of libdpd's ability to handle easily four-index
10
** quantities, including symmetry. This version requires
11
** significantly less disk space (ca. 1/2) than the original code,
12
** and is often much faster because of its reduced I/O requirements.
14
** This version of the code can do RHF, ROHF, and UHF transformations
15
** that are compatible with all the coupled cluster codes, including
18
** Remaining tasks to achieve full replacement of transqt v.1:
19
** (1) Add reordering arrays needed for DETCI.
20
** (2) Add partial transforms for MP2 and MP2-R12.
21
** (3) Replace the backtransformation. (I want to do this with
22
** symmetry, though, so there's no hurry here.)
31
#include <libipv1/ip_lib.h>
32
#include <libciomr/libciomr.h>
33
#include <libpsio/psio.h>
34
#include <libchkpt/chkpt.h>
35
#include <libiwl/iwl.h>
37
#include <libdpd/dpd.h>
41
/* Function prototypes */
42
void init_io(int argc, char *argv[]);
45
void get_params(void);
46
void get_moinfo(void);
50
int **cacheprep_rhf(int level, int *cachefiles);
51
void cachedone_rhf(int **);
53
int file_build_presort(dpdfile4 *, int, double, long int, int,
54
int, double *, double *, double *, double *, int);
55
void transtwo_rhf(void);
56
void transone(int m, int n, double *input, double *output, double **C, int nc, int *order, int *ioff);
57
void semicanonical_fock(void);
59
main(int argc, char *argv[])
61
int nso, nmo, ntri_so, ntri_mo, nirreps;
62
int **cachelist, *cachefiles;
65
double *H, *D, *F, *oei;
66
double *H_a, *H_b, *D_a, *D_b, *F_a, *F_b;
67
double ***C, ***C_a, ***C_b;
77
if(params.semicanonical) semicanonical_fock();
83
ntri_so = nso*(nso+1)/2;
84
ntri_mo = nmo*(nmo+1)/2;
85
nirreps = moinfo.nirreps;
87
cachefiles = init_int_array(PSIO_MAXUNIT);
88
cachelist = cacheprep_rhf(params.cachelev, cachefiles); /* really just a placeholder */
90
dpd_init(0, nirreps, params.memory, 0, cachefiles, cachelist,
91
NULL, 2, moinfo.sopi, moinfo.sosym, moinfo.mopi, moinfo.mosym);
93
/*** Starting one-electron transforms and presort ***/
95
/* For the one-electron integral transform, the full MO list is best */
96
if(params.ref == 0 || params.ref == 1) {
97
C = (double ***) malloc(1 * sizeof(double **));
98
chkpt_init(PSIO_OPEN_OLD);
99
C[0] = chkpt_rd_scf();
103
C_a = (double ***) malloc(1 * sizeof(double **));
104
C_b = (double ***) malloc(1 * sizeof(double **));
105
chkpt_init(PSIO_OPEN_OLD);
106
C_a[0] = chkpt_rd_alpha_scf();
107
C_b[0] = chkpt_rd_beta_scf();
111
/* build the frozen-core density (RHF) */
112
offset = init_int_array(nirreps);
113
for(h=1; h < nirreps; h++)
114
offset[h] = offset[h-1] + moinfo.sopi[h-1];
116
if(params.ref == 0 || params.ref == 1) { /* RHF/ROHF */
117
D = init_array(ntri_so);
118
for(h=0; h < nirreps; h++)
119
for(p=offset[h]; p < offset[h]+moinfo.sopi[h]; p++)
120
for(q=offset[h]; q <=p; q++) {
122
for(i=offset[h]; i < offset[h]+moinfo.frdocc[h]; i++)
123
D[pq] += C[0][p][i] * C[0][q][i];
125
if(params.print_lvl > 2) {
126
fprintf(outfile, "\n\tFrozen-core density (SO):\n");
127
print_array(D, nso, outfile);
131
D_a = init_array(ntri_so);
132
D_b = init_array(ntri_so);
133
for(h=0; h < nirreps; h++)
134
for(p=offset[h]; p < offset[h]+moinfo.sopi[h]; p++)
135
for(q=offset[h]; q <=p; q++) {
137
for(i=offset[h]; i < offset[h]+moinfo.frdocc[h]; i++) {
138
D_a[pq] += C_a[0][p][i] * C_a[0][q][i];
139
D_b[pq] += C_b[0][p][i] * C_b[0][q][i];
142
if(params.print_lvl > 2) {
143
fprintf(outfile, "\n\tAlpha Frozen-core density (SO):\n");
144
print_array(D_a, nso, outfile);
145
fprintf(outfile, "\n\tBeta Frozen-core density (SO):\n");
146
print_array(D_b, nso, outfile);
152
/* pre-sort the SO-basis two-electron integrals and generate the fzc operator(s) */
153
if(params.ref == 0 || params.ref == 1)
154
F = init_array(ntri_so);
156
F_a = init_array(ntri_so);
157
F_b = init_array(ntri_so);
161
if(params.print_lvl) {
162
fprintf(outfile, "\n\tPresorting SO-basis two-electron integrals.\n");
165
psio_open(PSIF_SO_PRESORT, 0);
166
dpd_file4_init(&I, PSIF_SO_PRESORT, 0, 3, 3, "SO Ints (pq,rs)");
167
if(params.ref == 0 || params.ref == 1)
168
file_build_presort(&I, PSIF_SO_TEI, params.tolerance, params.memory,
169
!params.delete_tei, moinfo.nfzc, D, NULL, F, NULL, params.ref);
171
file_build_presort(&I, PSIF_SO_TEI, params.tolerance, params.memory,
172
!params.delete_tei, moinfo.nfzc, D_a, D_b, F_a, F_b, params.ref);
174
psio_close(PSIF_SO_PRESORT, 1);
175
timer_off("presort");
177
/* read the bare one-electron integrals */
178
oei = init_array(ntri_so);
179
H = init_array(ntri_so);
180
stat = iwl_rdone(PSIF_OEI, PSIF_SO_T, H, ntri_so, 0, 0, outfile);
181
stat = iwl_rdone(PSIF_OEI, PSIF_SO_V, oei, ntri_so, 0, 0, outfile);
182
for(pq=0; pq < ntri_so; pq++)
185
/* add the remaining one-electron terms to the fzc operator(s) */
186
if(params.ref == 0 || params.ref == 1) {
187
for(pq=0; pq < ntri_so; pq++)
191
for(pq=0; pq < ntri_so; pq++) {
197
/* compute the frozen-core energy and write it to the chkpt file*/
199
if(params.ref == 0 || params.ref == 1) { /* RHF/ROHF */
200
for(p=0; p < nso; p++) {
202
efzc += D[pq] * (H[pq] + F[pq]);
203
for(q=0; q < p; q++) {
205
efzc += 2.0 * D[pq] * (H[pq] + F[pq]);
210
for(p=0; p < nso; p++) {
212
efzc += 0.5 * D_a[pq] * (H[pq] + F_a[pq]);
213
efzc += 0.5 * D_b[pq] * (H[pq] + F_b[pq]);
214
for(q=0; q < p; q++) {
216
efzc += D_a[pq] * (H[pq] + F_a[pq]);
217
efzc += D_b[pq] * (H[pq] + F_b[pq]);
221
if(params.print_lvl) {
222
fprintf(outfile, "\tFrozen-core energy = %20.15f\n", efzc);
225
chkpt_init(PSIO_OPEN_OLD);
229
/* transform the bare one-electron integrals */
230
if(params.ref == 0 || params.ref == 1) {
231
transone(nso, nmo, H, oei, C[0], nmo, moinfo.pitzer2qt, ioff);
232
if(params.print_lvl > 2) {
233
fprintf(outfile, "\n\tOne-electron integrals (MO basis):\n");
234
print_array(oei, nmo, outfile);
236
iwl_wrtone(PSIF_OEI, PSIF_MO_OEI, ntri_mo, oei);
240
transone(nso, nmo, H, oei, C_a[0], nmo, moinfo.pitzer2qt_A, ioff);
241
if(params.print_lvl > 2) {
242
fprintf(outfile, "\n\tAlpha one-electron integrals (MO basis):\n");
243
print_array(oei, nmo, outfile);
245
iwl_wrtone(PSIF_OEI, PSIF_MO_A_OEI, ntri_mo, oei);
248
transone(nso, nmo, H, oei, C_b[0], nmo, moinfo.pitzer2qt_B, ioff);
249
if(params.print_lvl > 2) {
250
fprintf(outfile, "\n\tBeta one-electron integrals (MO basis):\n");
251
print_array(oei, nmo, outfile);
253
iwl_wrtone(PSIF_OEI, PSIF_MO_B_OEI, ntri_mo, oei);
256
/* transform the frozen-core operator */
257
if(params.ref == 0 || params.ref == 1) { /* RHF/ROHF */
258
transone(nso, nmo, F, oei, C[0], nmo, moinfo.pitzer2qt, ioff);
259
if(params.print_lvl > 2) {
260
fprintf(outfile, "\n\tFrozen-core operator (MO basis):\n");
261
print_array(oei, nmo, outfile);
263
iwl_wrtone(PSIF_OEI, PSIF_MO_FZC, ntri_mo, oei);
268
transone(nso, nmo, F_a, oei, C_a[0], nmo, moinfo.pitzer2qt_A, ioff);
269
if(params.print_lvl > 2) {
270
fprintf(outfile, "\n\tAlpha frozen-core operator (MO basis):\n");
271
print_array(oei, nmo, outfile);
273
iwl_wrtone(PSIF_OEI, PSIF_MO_A_FZC, ntri_mo, oei);
276
transone(nso, nmo, F_b, oei, C_b[0], nmo, moinfo.pitzer2qt_B, ioff);
277
if(params.print_lvl > 2) {
278
fprintf(outfile, "\n\tBeta frozen-core operator (MO basis):\n");
279
print_array(oei, nmo, outfile);
281
iwl_wrtone(PSIF_OEI, PSIF_MO_B_FZC, ntri_mo, oei);
286
if(params.ref == 0 || params.ref == 1) {
303
/*** One-electron transforms complete ***/
305
/*** Starting two-electron transforms ***/
307
if(params.ref == 0 || params.ref == 1) transtwo_rhf();
310
/*** Two-electron transforms complete ***/
314
cachedone_rhf(cachelist);
322
exit(PSI_RETURN_SUCCESS);
325
char *gprgid(void) { char *prgid = "TRANSQT"; return (prgid); }
327
void init_io(int argc, char *argv[])
330
extern char *gprgid(void);
332
int num_extra_args = 0;
334
extra_args = (char **) malloc(argc*sizeof(char *));
336
params.print_lvl = 1;
337
for (i=1; i<argc; i++) {
338
if (!strcmp(argv[i], "--quiet"))
339
params.print_lvl = 0;
341
extra_args[num_extra_args++] = argv[i];
344
psi_start(num_extra_args, extra_args, 0);
346
progid = (char *) malloc(strlen(gprgid())+2);
347
sprintf(progid, ":%s", gprgid());
351
if(params.print_lvl) tstart(outfile);
354
psio_open(CC_INFO, PSIO_OPEN_NEW);
359
if(params.print_lvl) {
360
fprintf(outfile, "\n");
361
fprintf(outfile,"\t**************************************************\n");
362
fprintf(outfile,"\t* TRANSQT: Program to transform integrals from *\n");
363
fprintf(outfile,"\t* the SO basis to the MO basis. *\n");
364
fprintf(outfile,"\t* *\n");
365
fprintf(outfile,"\t* Daniel, David, & Justin *\n");
366
fprintf(outfile,"\t**************************************************\n");
367
fprintf(outfile, "\n");
373
psio_close(CC_INFO,1);
375
if(params.print_lvl) tstop(outfile);
382
ioff = init_int_array(IOFF_MAX);
383
for(i=1; i < IOFF_MAX; i++)
384
ioff[i] = ioff[i-1] + i;
added
removed
src/bin/transqt2/transqt.c
