4
** Return values of one and two-electron integrals
7
** Center for Computational Quantum Chemistry
8
** University of Georgia
10
** Updated 3/18/95 to exclude frozen virtual orbitals.
11
** Updated 3/28/95 to exclude frozen core orbitals.
16
#include <libiwl/iwl.h>
17
#include <libciomr/libciomr.h>
24
#define MIN0(a,b) (((a)<(b)) ? (a) : (b))
25
#define MAX0(a,b) (((a)>(b)) ? (a) : (b))
30
int i, j, ij, k, l, kl, ijkl, ijij;
31
int nmotri, nmotri_full;
33
extern void check_energy(double *H, double *twoel_ints, int *docc,
34
int *frozen_docc, int fzc_flag, double escf, double enuc, double efzc,
35
int nirreps, int *reorder, int *opi, FILE *outfile);
37
double *tmp_onel_ints;
39
/* allocate memory for one and two electron integrals */
40
nmotri_full = (CalcInfo.nmo * (CalcInfo.nmo + 1)) / 2;
41
nmotri = (CalcInfo.num_ci_orbs * (CalcInfo.num_ci_orbs + 1)) / 2 ;
42
CalcInfo.onel_ints = (double *) init_array(nmotri) ;
43
CalcInfo.twoel_ints = (double *) init_array(nmotri * (nmotri + 1) / 2);
44
CalcInfo.maxK = (double *) init_array(CalcInfo.num_ci_orbs);
46
/* replace old iwl_rdone call with the new one which filters
47
* note we must always filter the one-elec ints in the new scheme,
48
* although the two-elec ints only need filtering for first derivs
52
/* iwl_rdone(Parameters.oei_file, nmotri, CalcInfo.onel_ints,
53
&(CalcInfo.efzc), Parameters.oei_erase); */
55
tmp_onel_ints = init_array(nmotri_full);
56
iwl_rdone(Parameters.oei_file, PSIF_MO_FZC, tmp_onel_ints, nmotri_full,
57
Parameters.oei_erase, (Parameters.print_lvl>4), outfile);
58
filter(tmp_onel_ints, CalcInfo.onel_ints, ioff, CalcInfo.nmo,
59
CalcInfo.num_fzc_orbs, CalcInfo.num_fzv_orbs);
62
iwl_rdtwo(Parameters.tei_file, CalcInfo.twoel_ints, ioff, CalcInfo.nmo,
63
Parameters.filter_ints ? CalcInfo.num_fzc_orbs : 0,
64
Parameters.filter_ints ? CalcInfo.num_fzv_orbs : 0,
65
(Parameters.print_lvl>4), outfile);
67
/* Determine maximum K integral for use in averaging the diagonal */
68
/* Hamiltonian matrix elements over spin-coupling set */
69
if (Parameters.hd_ave) {
70
for(i=0; i<CalcInfo.num_ci_orbs; i++)
71
for(j=0; j<CalcInfo.num_ci_orbs; j++) {
72
/* if (i==j) continue; */
73
ij = ioff[MAX0(i,j)] + MIN0(i,j);
75
value = CalcInfo.twoel_ints[ijij];
76
if(value > CalcInfo.maxK[i]) CalcInfo.maxK[i] = value;
78
for(i=0; i<CalcInfo.num_ci_orbs; i++) {
79
if(CalcInfo.maxK[i] > CalcInfo.maxKlist)
80
CalcInfo.maxKlist = CalcInfo.maxK[i];
81
if (Parameters.print_lvl > 4)
82
fprintf(outfile,"maxK[%d] = %lf\n",i, CalcInfo.maxK[i]);
86
if (Parameters.print_lvl > 4) {
87
fprintf(outfile, "\nOne-electron integrals\n") ;
88
for (i=0, ij=0; i<CalcInfo.num_ci_orbs; i++) {
89
for (j=0; j<=i; j++, ij++) {
90
fprintf(outfile, "h(%d)(%d) = %11.7lf\n", i, j,
91
CalcInfo.onel_ints[ij]) ;
94
fprintf(outfile, "\n") ;
97
if (Parameters.print_lvl > 4) {
98
fprintf(outfile, "\nmaxKlist = %lf\n",CalcInfo.maxKlist);
99
fprintf(outfile, "\nTwo-electron integrals\n");
100
for (i=0; i<CalcInfo.num_ci_orbs; i++) {
101
for (j=0; j<=i; j++) {
102
ij = ioff[MAX0(i,j)] + MIN0(i,j) ;
103
for (k=0; k<=i; k++) {
104
for (l=0; l<=k; l++) {
105
kl = ioff[MAX0(k,l)] + MIN0(k,l) ;
106
ijkl = ioff[MAX0(ij,kl)] + MIN0(ij,kl) ;
107
fprintf(outfile, "%2d %2d %2d %2d (%4d) = %10.6lf\n",
108
i, j, k, l, ijkl, CalcInfo.twoel_ints[ijkl]);
115
if (Parameters.print_lvl) {
116
check_energy(CalcInfo.onel_ints, CalcInfo.twoel_ints, CalcInfo.docc,
117
CalcInfo.frozen_docc, Parameters.fzc, CalcInfo.escf, CalcInfo.enuc,
118
CalcInfo.efzc, CalcInfo.nirreps, CalcInfo.reorder,
119
CalcInfo.orbs_per_irr, outfile);
126
double get_onel(int i, int j)
133
value = CalcInfo.onel_ints[ij] ;
138
value = CalcInfo.onel_ints[ij] ;
141
return(CalcInfo.onel_ints[ij]) ;
145
double get_twoel(int i, int j, int k, int l)
149
ij = ioff[MAX0(i,j)] ;
151
kl = ioff[MAX0(k,l)] ;
153
ijkl = ioff[MAX0(ij,kl)] ;
154
ijkl += MIN0(ij,kl) ;
157
return(CalcInfo.twoel_ints[ijkl]) ;
163
** tf_onel_ints(): Function lumps together one-electron contributions
164
** so that h'_{ij} = h_{ij} - 1/2 SUM_k (ik|kj)
165
** The term h' arises in the calculation of sigma1 and sigma2 via
166
** equation (20) of Olsen, Roos, et. al. JCP 1988
169
void tf_onel_ints(int printflag, FILE *outfile)
171
int i, j, k, ij, ik, kj, ikkj ;
173
double *tei, *teptr ;
177
/* set up some shorthand notation (speed up access) */
178
nbf = CalcInfo.num_ci_orbs ;
179
tei = CalcInfo.twoel_ints ;
180
ntri = (nbf * (nbf + 1)) / 2;
182
/* ok, new special thing for CASSCF...if there are *no* excitations
183
into restricted orbitals, and if Parameters.fci=TRUE, then we
184
do *not* want to sum over the restricted virts in h' or else
185
we would need to account for RAS-out-of-space contributions
186
(requiring fci=false).
188
if (Parameters.fci && (nbf > Parameters.ras3_lvl) &&
189
Parameters.ras34_max == 0)
190
nbf = Parameters.ras3_lvl;
192
/* allocate space for the new array */
193
CalcInfo.tf_onel_ints = init_array(ntri) ;
195
/* fill up the new array */
196
for (i=0,ij=0; i<nbf; i++)
197
for (j=0; j<=i; j++) {
198
tval = CalcInfo.onel_ints[ij] ;
200
for (k=0; k<nbf; k++) {
201
ik = ioff[MAX0(i,k)] + MIN0(i,k) ;
202
kj = ioff[MAX0(k,j)] + MIN0(k,j) ;
203
ikkj = ioff[ik] + kj ;
205
tval -= 0.5 * (*teptr) ;
208
CalcInfo.tf_onel_ints[ij++] = tval ;
210
/* print if necessary */
212
fprintf(outfile, "\nh' matrix\n") ;
213
print_array(CalcInfo.tf_onel_ints, nbf, outfile) ;
214
fprintf(outfile, "\n") ;
221
** form_gmat(): Form the g matrix necessary for restriction to the RAS
222
** subspaces (i.e. to eliminate contributions of out-of-space terms).
223
** See equations (28-29) in Olsen, Roos, et. al. JCP 1988
226
void form_gmat(int printflag, FILE *outfile)
231
int i, j, k, ij, ii, ik, kj, ikkj, iiij ;
234
/* set up some shorthand notation (speed up access) */
235
nbf = CalcInfo.num_ci_orbs ;
236
oei = CalcInfo.onel_ints ;
237
tei = CalcInfo.twoel_ints ;
239
/* allocate space for the new array */
240
/* CalcInfo.gmat = init_matrix(nbf, nbf) ; */
241
/* why not use init_blockmatix here? */
242
CalcInfo.gmat = (double **) malloc (nbf * sizeof(double *));
243
CalcInfo.gmat[0] = (double *) malloc (nbf * nbf * sizeof(double));
244
for (i=1; i<nbf; i++) {
245
CalcInfo.gmat[i] = CalcInfo.gmat[i-1] + nbf;
248
/* fill up the new array */
249
for (i=0; i<nbf; i++) {
250
for (j=i+1; j<nbf; j++) {
253
for (k=0; k<i; k++) {
256
ikkj = ioff[kj] + ik ;
259
CalcInfo.gmat[i][j] = tval ;
263
for (i=0, ij=0; i<nbf; i++) {
264
for (j=0; j<=i; j++,ij++) {
266
for (k=0; k<i; k++) {
268
kj = ioff[MAX0(k,j)] + MIN0(k,j) ;
269
ikkj = ioff[ik] + kj ;
273
iiij = ioff[ii] + ij ;
274
if (i==j) tval -= 0.5 * tei[iiij] ;
275
else tval -= tei[iiij] ;
276
CalcInfo.gmat[i][j] = tval ;
281
fprintf(outfile, "\ng matrix\n") ;
282
print_mat(CalcInfo.gmat, nbf, nbf, outfile) ;
283
fprintf(outfile, "\n") ;
added
removed
src/bin/detci/ints.c
