3
\brief Enter brief description of file here
9
#include<libipv1/ip_lib.h>
10
#include<libiwl/iwl.h>
11
#include<libciomr/libciomr.h>
12
#include<libint/libint.h>
24
/*!---------------------------------------------------------------------------
25
This function computes AO overlap, dipole and quadrupole moment integrals,
26
and nabla integrals and writes them out to disk
27
---------------------------------------------------------------------------*/
30
struct coordinates PA, PB, AB, A, B;
31
struct shell_pair *sp;
32
struct unique_shell_pair *usp;
33
register int i, j, k, l, ii, jj, kk, ll;
38
int l1, l2, m1, m2, n1, n2;
39
int ioffset, joffset ;
53
double **stemp, **mxtemp, **mytemp, **mztemp, **qxxtemp, **qxytemp, **qxztemp,
54
**qyytemp, **qyztemp, **qzztemp, **nxtemp, **nytemp, **nztemp,
56
double *S, *MX, *MY, *MZ, *QXX, *QXY, *QXZ, *QYY, *QYZ, *QZZ, *NX, *NY, *NZ;
57
double inorm, jnorm, over_pf;
58
double *ptr1, *ptr2, norm1, norm12;
59
double **OIX, **OIY, **OIZ;
63
/*--- allocate room for the one-e matrices ---*/
64
dimension = ioff[BasisSet.num_ao];
65
S = init_array(dimension);
66
MX = init_array(dimension);
67
MY = init_array(dimension);
68
MZ = init_array(dimension);
69
QXX = init_array(dimension);
70
QXY = init_array(dimension);
71
QXZ = init_array(dimension);
72
QYY = init_array(dimension);
73
QYZ = init_array(dimension);
74
QZZ = init_array(dimension);
75
NX = init_array(dimension);
76
NY = init_array(dimension);
77
NZ = init_array(dimension);
79
/*--- allocate storage for shell blocks of one electron integrals ---*/
80
dimension = ioff[BasisSet.max_am];
81
stemp = block_matrix(dimension,dimension);
82
mxtemp = block_matrix(dimension,dimension);
83
mytemp = block_matrix(dimension,dimension);
84
mztemp = block_matrix(dimension,dimension);
85
qxxtemp = block_matrix(dimension,dimension);
86
qxytemp = block_matrix(dimension,dimension);
87
qxztemp = block_matrix(dimension,dimension);
88
qyytemp = block_matrix(dimension,dimension);
89
qyztemp = block_matrix(dimension,dimension);
90
qzztemp = block_matrix(dimension,dimension);
91
nxtemp = block_matrix(dimension,dimension);
92
nytemp = block_matrix(dimension,dimension);
93
nztemp = block_matrix(dimension,dimension);
95
/* Note the "+1" -- we are computing dipole and quadrupole moment integrals here */
96
OIX = block_matrix(BasisSet.max_am+1,BasisSet.max_am+1);
97
OIY = block_matrix(BasisSet.max_am+1,BasisSet.max_am+1);
98
OIZ = block_matrix(BasisSet.max_am+1,BasisSet.max_am+1);
100
C = UserOptions.origin;
101
if(UserOptions.print_lvl >= PRINT_OEI)
102
fprintf(outfile, " Reference point for elec. mom. ints. = (%5.3f, %5.3f, %5.3f)\n", C.x, C.y, C.z);
104
for (si=0; si<BasisSet.num_shells; si++){
105
am_i = BasisSet.shells[si].am-1;
106
ni = ioff[BasisSet.shells[si].am];
107
A = Molecule.centers[BasisSet.shells[si].center-1];
108
ioffset = BasisSet.shells[si].fao - 1;
109
for (sj=0; sj<=si; sj++){
110
nj = ioff[BasisSet.shells[sj].am];
111
am_j = BasisSet.shells[sj].am-1;
112
B = Molecule.centers[BasisSet.shells[sj].center-1];
113
joffset = BasisSet.shells[sj].fao - 1;
115
sp = &(BasisSet.shell_pairs[si][sj]);
123
/*--- zero the temporary storage for accumulating contractions ---*/
124
memset(stemp[0],0,sizeof(double)*dimension*dimension);
125
memset(mxtemp[0],0,sizeof(double)*dimension*dimension);
126
memset(mytemp[0],0,sizeof(double)*dimension*dimension);
127
memset(mztemp[0],0,sizeof(double)*dimension*dimension);
128
memset(qxxtemp[0],0,sizeof(double)*dimension*dimension);
129
memset(qxytemp[0],0,sizeof(double)*dimension*dimension);
130
memset(qxztemp[0],0,sizeof(double)*dimension*dimension);
131
memset(qyytemp[0],0,sizeof(double)*dimension*dimension);
132
memset(qyztemp[0],0,sizeof(double)*dimension*dimension);
133
memset(qzztemp[0],0,sizeof(double)*dimension*dimension);
134
memset(nxtemp[0],0,sizeof(double)*dimension*dimension);
135
memset(nytemp[0],0,sizeof(double)*dimension*dimension);
136
memset(nztemp[0],0,sizeof(double)*dimension*dimension);
138
/*--- contract by primitives here ---*/
139
for (i = 0; i < BasisSet.shells[si].n_prims; i++) {
141
inorm = sp->inorm[i];
142
for (j = 0; j < BasisSet.shells[sj].n_prims; j++) {
144
gam = sp->gamma[i][j];
145
jnorm = sp->jnorm[j];
146
PA.x = sp->PA[i][j][0];
147
PA.y = sp->PA[i][j][1];
148
PA.z = sp->PA[i][j][2];
149
PB.x = sp->PB[i][j][0];
150
PB.y = sp->PB[i][j][1];
151
PB.z = sp->PB[i][j][2];
153
over_pf = exp(-a1*a2*ab2*oog)*sqrt(M_PI*oog)*M_PI*oog*inorm*jnorm;
155
OI_OSrecurs(OIX,OIY,OIZ,PA,PB,gam,am_i+1,am_j+1);
157
/*--- create all am components of si ---*/
159
for(ii = 0; ii <= am_i; ii++){
161
for(jj = 0; jj <= ii; jj++){
164
/*--- create all am components of sj ---*/
166
for(kk = 0; kk <= am_j; kk++){
168
for(ll = 0; ll <= kk; ll++){
172
x00 = OIX[l1][l2]; y00 = OIY[m1][m2]; z00 = OIZ[n1][n2];
173
x01 = OIX[l1][l2+1]; y01 = OIY[m1][m2+1]; z01 = OIZ[n1][n2+1];
174
x10 = OIX[l1+1][l2]; y10 = OIY[m1+1][m2]; z10 = OIZ[n1+1][n2];
175
x11 = OIX[l1+1][l2+1]; y11 = OIY[m1+1][m2+1]; z11 = OIZ[n1+1][n2+1];
176
stemp[ai][aj] += over_pf*x00*y00*z00;
177
/*--- electrons have negative charge ---*/
178
mxtemp[ai][aj] -= over_pf*(x01+x00*(B.x-C.x))*y00*z00;
179
mytemp[ai][aj] -= over_pf*x00*(y01+y00*(B.y-C.y))*z00;
180
mztemp[ai][aj] -= over_pf*x00*y00*(z01+z00*(B.z-C.z));
181
qxxtemp[ai][aj] -= over_pf*(x11 + x10*(B.x-C.x) + x01*(A.x-C.x)
182
+ x00*(A.x-C.x)*(B.x-C.x))*y00*z00;
183
qyytemp[ai][aj] -= over_pf*(y11 + y10*(B.y-C.y) + y01*(A.y-C.y)
184
+ y00*(A.y-C.y)*(B.y-C.y))*x00*z00;
185
qzztemp[ai][aj] -= over_pf*(z11 + z10*(B.z-C.z) + z01*(A.z-C.z)
186
+ z00*(A.z-C.z)*(B.z-C.z))*x00*y00;
187
qxytemp[ai][aj] -= over_pf*(x01+x00*(B.x-C.x))*(y01+y00*(B.y-C.y))*z00;
188
qxztemp[ai][aj] -= over_pf*(x01+x00*(B.x-C.x))*y00*(z01+z00*(B.z-C.z));
189
qyztemp[ai][aj] -= over_pf*x00*(y01+y00*(B.y-C.y))*(z01+z00*(B.z-C.z));
193
nx += l2*OIX[l1][l2-1];
194
nxtemp[ai][aj] += nx*y00*z00*over_pf;
197
ny += m2*OIY[m1][m2-1];
198
nytemp[ai][aj] += x00*ny*z00*over_pf;
201
nz += n2*OIZ[n1][n2-1];
202
nztemp[ai][aj] += x00*y00*nz*over_pf;
209
} /*--- end cartesian components for (si,sj) with primitives (i,j) ---*/
211
} /*--- end primitive contraction ---*/
213
/*--- Normalize the contracted integrals ---*/
214
ptr1 = GTOs.bf_norm[am_i];
215
ptr2 = GTOs.bf_norm[am_j];
216
for(i=0; i<ni; i++) {
218
for(j=0; j<nj; j++) {
219
norm12 = norm1*ptr2[j];
220
stemp[i][j] *= norm12;
221
mxtemp[i][j] *= norm12;
222
mytemp[i][j] *= norm12;
223
mztemp[i][j] *= norm12;
224
qxxtemp[i][j] *= norm12;
225
qxytemp[i][j] *= norm12;
226
qxztemp[i][j] *= norm12;
227
qyytemp[i][j] *= norm12;
228
qyztemp[i][j] *= norm12;
229
qzztemp[i][j] *= norm12;
230
nxtemp[i][j] *= norm12;
231
nytemp[i][j] *= norm12;
232
nztemp[i][j] *= norm12;
238
ij = INDEX(ioffset + i, joffset + j);
240
MX[ij] = mxtemp[i][j];
241
MY[ij] = mytemp[i][j];
242
MZ[ij] = mztemp[i][j];
243
QXX[ij] = qxxtemp[i][j];
244
QXY[ij] = qxytemp[i][j];
245
QXZ[ij] = qxztemp[i][j];
246
QYY[ij] = qyytemp[i][j];
247
QYZ[ij] = qyztemp[i][j];
248
QZZ[ij] = qzztemp[i][j];
249
NX[ij] = nxtemp[i][j];
250
NY[ij] = nytemp[i][j];
251
NZ[ij] = nztemp[i][j];
254
} /*--- This shell pair is done ---*/
256
/*--- flush it all away ---*/
261
dimension=ioff[BasisSet.num_ao];
262
iwl_wrtone(IOUnits.itapS_AO, PSIF_AO_S, dimension,S);
263
iwl_wrtone(IOUnits.itapMX_AO,PSIF_AO_MX,dimension,MX);
264
iwl_wrtone(IOUnits.itapMY_AO,PSIF_AO_MY,dimension,MY);
265
iwl_wrtone(IOUnits.itapMZ_AO,PSIF_AO_MZ,dimension,MZ);
266
iwl_wrtone(IOUnits.itapQXX_AO,PSIF_AO_QXX,dimension,QXX);
267
iwl_wrtone(IOUnits.itapQXY_AO,PSIF_AO_QXY,dimension,QXY);
268
iwl_wrtone(IOUnits.itapQXZ_AO,PSIF_AO_QXZ,dimension,QXZ);
269
iwl_wrtone(IOUnits.itapQYY_AO,PSIF_AO_QYY,dimension,QYY);
270
iwl_wrtone(IOUnits.itapQYZ_AO,PSIF_AO_QYZ,dimension,QYZ);
271
iwl_wrtone(IOUnits.itapQZZ_AO,PSIF_AO_QZZ,dimension,QZZ);
272
iwl_wrtone(IOUnits.itapNablaX_AO,PSIF_AO_NablaX,dimension,NX);
273
iwl_wrtone(IOUnits.itapNablaY_AO,PSIF_AO_NablaY,dimension,NY);
274
iwl_wrtone(IOUnits.itapNablaZ_AO,PSIF_AO_NablaZ,dimension,NZ);
275
if (UserOptions.print_lvl >= PRINT_OEI) {
276
fprintf(outfile," -Overlap AO integrals:\n\n");
277
print_array(S,BasisSet.num_ao,outfile);
278
fprintf(outfile," -mu(x) AO integrals:\n\n");
279
print_array(MX,BasisSet.num_ao,outfile);
280
fprintf(outfile," -mu(y) AO integrals:\n\n");
281
print_array(MY,BasisSet.num_ao,outfile);
282
fprintf(outfile," -mu(z) AO integrals:\n\n");
283
print_array(MZ,BasisSet.num_ao,outfile);
284
fprintf(outfile," -q(xx) AO integrals:\n\n");
285
print_array(QXX,BasisSet.num_ao,outfile);
286
fprintf(outfile," -q(xy) AO integrals:\n\n");
287
print_array(QXY,BasisSet.num_ao,outfile);
288
fprintf(outfile," -q(xz) AO integrals:\n\n");
289
print_array(QXZ,BasisSet.num_ao,outfile);
290
fprintf(outfile," -q(yy) AO integrals:\n\n");
291
print_array(QYY,BasisSet.num_ao,outfile);
292
fprintf(outfile," -q(yz) AO integrals:\n\n");
293
print_array(QYZ,BasisSet.num_ao,outfile);
294
fprintf(outfile," -q(zz) AO integrals:\n\n");
295
print_array(QZZ,BasisSet.num_ao,outfile);
296
fprintf(outfile," -Nabla_x AO integrals:\n\n");
297
print_array(NX,BasisSet.num_ao,outfile);
298
fprintf(outfile," -Nabla_y AO integrals:\n\n");
299
print_array(NY,BasisSet.num_ao,outfile);
300
fprintf(outfile," -Nabla_z AO integrals:\n\n");
301
print_array(NZ,BasisSet.num_ao,outfile);
302
fprintf(outfile,"\n");
added
removed
src/bin/cints/OEProp_Ints/moment_ints.cc
