6
#include<libipv1/ip_lib.h>
8
#include<libciomr/libciomr.h>
9
#include<libint/libint.h>
15
#include"taylor_fm_eval.h"
20
/*-----------------------------------------------------------------------------------------
21
This function computes some derivatives of integrals over GIAO Gaussians with respect to
22
E and B fields (at E=0, B=0, i.e. in absence of external fields) and writes them out.
25
d2h/dBdE = -0.5 i <mu| (AB x r) r |nu>
27
d h/dB = 0.5 <mu| L + i (AB x r) h |nu>
30
d h/dE = - <mu| r |nu> -- simply electric dipole integral, already computed, hence
32
-----------------------------------------------------------------------------------------*/
33
extern "C" void giao_oe_deriv()
37
init_Taylor_Fm_Eval(BasisSet.max_am*4+1,UserOptions.cutoff);
40
/*--- allocate room for the one-e matrices ---*/
41
double** D2HDBXDEX = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
42
double** D2HDBXDEY = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
43
double** D2HDBXDEZ = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
44
double** D2HDBYDEX = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
45
double** D2HDBYDEY = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
46
double** D2HDBYDEZ = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
47
double** D2HDBZDEX = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
48
double** D2HDBZDEY = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
49
double** D2HDBZDEZ = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
50
double** DHDBX = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
51
double** DHDBY = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
52
double** DHDBZ = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
53
double** DSDBX = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
54
double** DSDBY = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
55
double** DSDBZ = block_matrix(BasisSet.num_ao,BasisSet.num_ao);
57
/*--- allocate storage for shell blocks of one electron integrals ---*/
58
int dimension = ioff[BasisSet.max_am];
59
double** d2HdBxdExtemp = block_matrix(dimension,dimension);
60
double** d2HdBxdEytemp = block_matrix(dimension,dimension);
61
double** d2HdBxdEztemp = block_matrix(dimension,dimension);
62
double** d2HdBydExtemp = block_matrix(dimension,dimension);
63
double** d2HdBydEytemp = block_matrix(dimension,dimension);
64
double** d2HdBydEztemp = block_matrix(dimension,dimension);
65
double** d2HdBzdExtemp = block_matrix(dimension,dimension);
66
double** d2HdBzdEytemp = block_matrix(dimension,dimension);
67
double** d2HdBzdEztemp = block_matrix(dimension,dimension);
68
double** dHdBxtemp = block_matrix(dimension,dimension);
69
double** dHdBytemp = block_matrix(dimension,dimension);
70
double** dHdBztemp = block_matrix(dimension,dimension);
71
double** dSdBxtemp = block_matrix(dimension,dimension);
72
double** dSdBytemp = block_matrix(dimension,dimension);
73
double** dSdBztemp = block_matrix(dimension,dimension);
75
/* Note the "+1" on the row dimension -- we are computing kinetic energy
76
integrals with one extra quantum in the bra function here */
77
double** OIX = block_matrix(BasisSet.max_am+1,BasisSet.max_am+2);
78
double** OIY = block_matrix(BasisSet.max_am+1,BasisSet.max_am+2);
79
double** OIZ = block_matrix(BasisSet.max_am+1,BasisSet.max_am+2);
80
/* same here -- need to compute integrals with extra quantum in the bra */
81
int indmax_0 = (BasisSet.max_am-1)*BasisSet.max_am*BasisSet.max_am+1;
82
int indmax_p1 = (BasisSet.max_am)*(BasisSet.max_am+1)*(BasisSet.max_am+1)+1;
83
double*** AI0 = init_box(indmax_p1,indmax_0,2*BasisSet.max_am+3);
85
for (int si=0; si<BasisSet.num_shells; si++){
86
int am_i = BasisSet.shells[si].am-1;
87
int ni = ioff[BasisSet.shells[si].am];
88
struct coordinates A = Molecule.centers[BasisSet.shells[si].center-1];
89
int ioffset = BasisSet.shells[si].fao - 1;
91
int iym = am_i+2; // Note +2 instead of +1 here -- we'll be computing nuclear attraction
92
// integrals with extra quantum in bra!
95
for (int sj=0; sj<BasisSet.num_shells; sj++){
96
int nj = ioff[BasisSet.shells[sj].am];
97
int am_j = BasisSet.shells[sj].am-1;
98
struct coordinates B = Molecule.centers[BasisSet.shells[sj].center-1];
99
int joffset = BasisSet.shells[sj].fao - 1;
104
struct shell_pair* sp = &(BasisSet.shell_pairs[si][sj]);
105
struct coordinates AB;
109
double ab2 = AB.x * AB.x;
113
/*--- zero the temporary storage for accumulating contractions ---*/
114
memset(d2HdBxdExtemp[0],0,sizeof(double)*dimension*dimension);
115
memset(d2HdBxdEytemp[0],0,sizeof(double)*dimension*dimension);
116
memset(d2HdBxdEztemp[0],0,sizeof(double)*dimension*dimension);
117
memset(d2HdBydExtemp[0],0,sizeof(double)*dimension*dimension);
118
memset(d2HdBydEytemp[0],0,sizeof(double)*dimension*dimension);
119
memset(d2HdBydEztemp[0],0,sizeof(double)*dimension*dimension);
120
memset(d2HdBzdExtemp[0],0,sizeof(double)*dimension*dimension);
121
memset(d2HdBzdEytemp[0],0,sizeof(double)*dimension*dimension);
122
memset(d2HdBzdEztemp[0],0,sizeof(double)*dimension*dimension);
123
memset(dHdBxtemp[0],0,sizeof(double)*dimension*dimension);
124
memset(dHdBytemp[0],0,sizeof(double)*dimension*dimension);
125
memset(dHdBztemp[0],0,sizeof(double)*dimension*dimension);
126
memset(dSdBxtemp[0],0,sizeof(double)*dimension*dimension);
127
memset(dSdBytemp[0],0,sizeof(double)*dimension*dimension);
128
memset(dSdBztemp[0],0,sizeof(double)*dimension*dimension);
130
/*--- contract by primitives here ---*/
131
for (int i = 0; i < BasisSet.shells[si].n_prims; i++) {
132
double a1 = sp->a1[i];
133
double inorm = sp->inorm[i];
134
for (int j = 0; j < BasisSet.shells[sj].n_prims; j++) {
135
double a2 = sp->a2[j];
136
double gam = sp->gamma[i][j];
137
double jnorm = sp->jnorm[j];
138
struct coordinates PA, PB;
139
PA.x = sp->PA[i][j][0];
140
PA.y = sp->PA[i][j][1];
141
PA.z = sp->PA[i][j][2];
142
PB.x = sp->PB[i][j][0];
143
PB.y = sp->PB[i][j][1];
144
PB.z = sp->PB[i][j][2];
145
double oog = 1.0/gam;
146
double over_pf = exp(-a1*a2*ab2*oog)*sqrt(M_PI*oog)*M_PI*oog*inorm*jnorm;
148
OI_OSrecurs(OIX,OIY,OIZ,PA,PB,gam,am_i+1,am_j+2);
150
/*--- create all am components of si ---*/
152
for(int ii = 0; ii <= am_i; ii++){
154
for(int jj = 0; jj <= ii; jj++){
157
/*--- create all am components of sj ---*/
159
for(int kk = 0; kk <= am_j; kk++){
161
for(int ll = 0; ll <= kk; ll++){
165
double x00 = OIX[l1][l2];
166
double y00 = OIY[m1][m2];
167
double z00 = OIZ[n1][n2];
169
double x01 = OIX[l1][l2+1];
170
double y01 = OIY[m1][m2+1];
171
double z01 = OIZ[n1][n2+1];
173
double x10 = OIX[l1+1][l2];
174
double y10 = OIY[m1+1][m2];
175
double z10 = OIZ[n1+1][n2];
177
double x11 = OIX[l1+1][l2+1];
178
double y11 = OIY[m1+1][m2+1];
179
double z11 = OIZ[n1+1][n2+1];
184
double x1 = x01 + x00*B.x;
185
double y1 = y01 + y00*B.y;
186
double z1 = z01 + z00*B.z;
187
double x2 = x11 + x10*B.x + x01*A.x + x00*A.x*B.x;
188
double y2 = y11 + y10*B.y + y01*A.y + y00*A.y*B.y;
189
double z2 = z11 + z10*B.z + z01*A.z + z00*A.z*B.z;
191
double dSdB_pfac = 0.5*over_pf;
195
dSdBxtemp[ai][aj] += dSdB_pfac * (AB.y * z - AB.z * y);
196
dSdBytemp[ai][aj] += dSdB_pfac * (AB.z * x - AB.x * z);
197
dSdBztemp[ai][aj] += dSdB_pfac * (AB.x * y - AB.y * x);
199
/* dSdBxtemp[ai][aj] += 2.0*dSdB_pfac * x;
200
dSdBytemp[ai][aj] += 2.0*dSdB_pfac * y;
201
dSdBztemp[ai][aj] += 2.0*dSdB_pfac * z;*/
203
double d2BE_pfac = -0.5*over_pf;
204
double qxx = x2*y0*z0;
205
double qxy = x1*y1*z0;
206
double qxz = x1*y0*z1;
208
double qyy = x0*y2*z0;
209
double qyz = x0*y1*z1;
212
double qzz = x0*y0*z2;
213
d2HdBxdExtemp[ai][aj] += d2BE_pfac * ( AB.y * qzx - AB.z * qyx);
214
d2HdBxdEytemp[ai][aj] += d2BE_pfac * ( AB.y * qzy - AB.z * qyy);
215
d2HdBxdEztemp[ai][aj] += d2BE_pfac * ( AB.y * qzz - AB.z * qyz);
216
d2HdBydExtemp[ai][aj] += d2BE_pfac * ( AB.z * qxx - AB.x * qzx);
217
d2HdBydEytemp[ai][aj] += d2BE_pfac * ( AB.z * qxy - AB.x * qzy);
218
d2HdBydEztemp[ai][aj] += d2BE_pfac * ( AB.z * qxz - AB.x * qzz);
219
d2HdBzdExtemp[ai][aj] += d2BE_pfac * ( AB.x * qyx - AB.y * qxx);
220
d2HdBzdEytemp[ai][aj] += d2BE_pfac * ( AB.x * qyy - AB.y * qxy);
221
d2HdBzdEztemp[ai][aj] += d2BE_pfac * ( AB.x * qyz - AB.y * qxz);
223
double nx = -2.0*a2*x01;
225
nx += l2*OIX[l1][l2-1];
226
double ny = -2.0*a2*y01;
228
ny += m2*OIY[m1][m2-1];
229
double nz = -2.0*a2*z01;
231
nz += n2*OIZ[n1][n2-1];
233
double dB_pfac = -0.5*over_pf;
234
dHdBxtemp[ai][aj] += dB_pfac * x00 * ( y01 * nz - z01 * ny);
235
dHdBytemp[ai][aj] += dB_pfac * y00 * ( z01 * nx - x01 * nz);
236
dHdBztemp[ai][aj] += dB_pfac * z00 * ( x01 * ny - y01 * nx);
238
double tx00 = a2*(2*l2+1)*OIX[l1][l2] - 2*a2*a2*OIX[l1][l2+2];
240
tx00 -= 0.5*l2*(l2-1)*OIX[l1][l2-2];
241
double ty00 = a2*(2*m2+1)*OIY[m1][m2] - 2*a2*a2*OIY[m1][m2+2];
243
ty00 -= 0.5*m2*(m2-1)*OIY[m1][m2-2];
244
double tz00 = a2*(2*n2+1)*OIZ[n1][n2] - 2*a2*a2*OIZ[n1][n2+2];
246
tz00 -= 0.5*n2*(n2-1)*OIZ[n1][n2-2];
248
double tx10 = a2*(2*l2+1)*(OIX[l1+1][l2] + A.x*OIX[l1][l2]) - 2*a2*a2*(OIX[l1+1][l2+2] + A.x*OIX[l1][l2+2]);
250
tx10 -= 0.5*l2*(l2-1)*(OIX[l1+1][l2-2] + A.x*OIX[l1][l2-2]);
251
double ty10 = a2*(2*m2+1)*(OIY[m1+1][m2] + A.y*OIY[m1][m2]) - 2*a2*a2*(OIY[m1+1][m2+2] + A.y*OIY[m1][m2+2]);
253
ty10 -= 0.5*m2*(m2-1)*(OIY[m1+1][m2-2] + A.y*OIY[m1][m2-2]);
254
double tz10 = a2*(2*n2+1)*(OIZ[n1+1][n2] + A.z*OIZ[n1][n2]) - 2*a2*a2*(OIZ[n1+1][n2+2] + A.z*OIZ[n1][n2+2]);
256
tz10 -= 0.5*n2*(n2-1)*(OIZ[n1+1][n2-2] + A.z*OIZ[n1][n2-2]);
258
double xT = tx10*y0*z0 + x1*ty00*z0 + x1*y0*tz00;
259
double yT = tx00*y1*z0 + x0*ty10*z0 + x0*y1*tz00;
260
double zT = tx00*y0*z1 + x0*ty00*z1 + x0*y0*tz10;
262
dB_pfac = 0.5*over_pf;
263
dHdBxtemp[ai][aj] += dB_pfac * ( AB.y * zT - AB.z * yT);
264
dHdBytemp[ai][aj] += dB_pfac * ( AB.z * xT - AB.x * zT);
265
dHdBztemp[ai][aj] += dB_pfac * ( AB.x * yT - AB.y * xT);
267
/* dHdBxtemp[ai][aj] += 2.0*dB_pfac * (tx00*y0*z0 + x0*ty00*z0+x0*y0*tz00);
268
dHdBztemp[ai][aj] += 2.0*dB_pfac * (tx00*y0*z0 + x0*ty00*z0+x0*y0*tz00);*/
270
/*double T = tx00*y0*z0 + x0*ty00*z0 + x0*y0*tz00;
271
dHdBytemp[ai][aj] += 2.0*dB_pfac*T;*/
278
} /*--- end cartesian components for (si,sj) with primitives (i,j) ---*/
280
for(int atom=0;atom<Molecule.num_atoms;atom++) {
281
struct coordinates PC;
282
PC.x = sp->P[i][j][0] - Molecule.centers[atom].x;
283
PC.y = sp->P[i][j][1] - Molecule.centers[atom].y;
284
PC.z = sp->P[i][j][2] - Molecule.centers[atom].z;
285
AI_OSrecurs(AI0,PA,PB,PC,gam,am_i+1,am_j);
287
double dB_pfac = 0.5 * over_pf * (-Molecule.centers[atom].Z_nuc);
290
for(int ii = 0; ii <= am_i; ii++){
292
for(int jj = 0; jj <= ii; jj++){
295
int iind_000 = n1*izm + m1*iym + l1*ixm;
296
int iind_100 = n1*izm + m1*iym + (l1+1)*ixm;
297
int iind_010 = n1*izm + (m1+1)*iym + l1*ixm;
298
int iind_001 = (n1+1)*izm + m1*iym + l1*ixm;
299
/*--- create all am components of sj ---*/
301
for(int kk = 0; kk <= am_j; kk++){
303
for(int ll = 0; ll <= kk; ll++){
307
int jind = n2*jzm + m2*jym + l2*jxm;
309
double V = AI0[iind_000][jind][0];
310
double xV = AI0[iind_100][jind][0] + A.x * V;
311
double yV = AI0[iind_010][jind][0] + A.y * V;
312
double zV = AI0[iind_001][jind][0] + A.z * V;
314
dHdBxtemp[ai][aj] += dB_pfac * ( AB.y * zV - AB.z * yV);
315
dHdBytemp[ai][aj] += dB_pfac * ( AB.z * xV - AB.x * zV);
316
dHdBztemp[ai][aj] += dB_pfac * ( AB.x * yV - AB.y * xV);
318
/* dHdBxtemp[ai][aj] += dB_pfac * xV;
319
dHdBytemp[ai][aj] += dB_pfac * yV;
320
dHdBztemp[ai][aj] += dB_pfac * zV;*/
327
} /*--- end cartesian components for (si,sj) with primitives (i,j) ---*/
331
} /*--- end primitive contraction ---*/
333
/*--- Normalize the contracted integrals ---*/
334
double* ptr1 = GTOs.bf_norm[am_i];
335
double* ptr2 = GTOs.bf_norm[am_j];
336
for(int i=0; i<ni; i++) {
337
double norm1 = ptr1[i];
338
for(int j=0; j<nj; j++) {
339
double norm12 = norm1*ptr2[j];
340
d2HdBxdExtemp[i][j] *= norm12;
341
d2HdBxdEytemp[i][j] *= norm12;
342
d2HdBxdEztemp[i][j] *= norm12;
343
d2HdBydExtemp[i][j] *= norm12;
344
d2HdBydEytemp[i][j] *= norm12;
345
d2HdBydEztemp[i][j] *= norm12;
346
d2HdBzdExtemp[i][j] *= norm12;
347
d2HdBzdEytemp[i][j] *= norm12;
348
d2HdBzdEztemp[i][j] *= norm12;
349
dHdBxtemp[i][j] *= norm12;
350
dHdBytemp[i][j] *= norm12;
351
dHdBztemp[i][j] *= norm12;
352
dSdBxtemp[i][j] *= norm12;
353
dSdBytemp[i][j] *= norm12;
354
dSdBztemp[i][j] *= norm12;
358
for(int i=0;i<ni;i++)
359
for(int j=0;j<nj;j++) {
360
int ii = ioffset + i;
361
int jj = joffset + j;
362
D2HDBXDEX[ii][jj] = d2HdBxdExtemp[i][j];
363
D2HDBXDEY[ii][jj] = d2HdBxdEytemp[i][j];
364
D2HDBXDEZ[ii][jj] = d2HdBxdEztemp[i][j];
365
D2HDBYDEX[ii][jj] = d2HdBydExtemp[i][j];
366
D2HDBYDEY[ii][jj] = d2HdBydEytemp[i][j];
367
D2HDBYDEZ[ii][jj] = d2HdBydEztemp[i][j];
368
D2HDBZDEX[ii][jj] = d2HdBzdExtemp[i][j];
369
D2HDBZDEY[ii][jj] = d2HdBzdEytemp[i][j];
370
D2HDBZDEZ[ii][jj] = d2HdBzdEztemp[i][j];
371
DHDBX[ii][jj] = dHdBxtemp[i][j];
372
DHDBY[ii][jj] = dHdBytemp[i][j];
373
DHDBZ[ii][jj] = dHdBztemp[i][j];
374
DSDBX[ii][jj] = dSdBxtemp[i][j];
375
DSDBY[ii][jj] = dSdBytemp[i][j];
376
DSDBZ[ii][jj] = dSdBztemp[i][j];
379
} /*--- This shell pair is done ---*/
381
/*--- flush it all away ---*/
386
free_box(AI0,indmax_p1,indmax_0);
387
dimension=BasisSet.num_ao*BasisSet.num_ao;
388
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_XX, dimension,D2HDBXDEX[0]);
389
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_XY, dimension,D2HDBXDEY[0]);
390
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_XZ, dimension,D2HDBXDEZ[0]);
391
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_YX, dimension,D2HDBYDEX[0]);
392
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_YY, dimension,D2HDBYDEY[0]);
393
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_YZ, dimension,D2HDBYDEZ[0]);
394
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_ZX, dimension,D2HDBZDEX[0]);
395
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_ZY, dimension,D2HDBZDEY[0]);
396
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_D2HDBDE_ZZ, dimension,D2HDBZDEZ[0]);
397
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DHDB_X, dimension,DHDBX[0]);
398
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DHDB_Y, dimension,DHDBY[0]);
399
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DHDB_Z, dimension,DHDBZ[0]);
400
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DSDB_X, dimension,DSDBX[0]);
401
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DSDB_Y, dimension,DSDBY[0]);
402
iwl_wrtone(IOUnits.itapOEInt_Misc, PSIF_AO_DSDB_Z, dimension,DSDBZ[0]);
403
if (UserOptions.print_lvl >= PRINT_OEI) {
404
fprintf(outfile," -dS/dB_x AO integrals:\n\n");
405
print_mat(DSDBX,BasisSet.num_ao,BasisSet.num_ao,outfile);
406
fprintf(outfile," -dS/dB_y AO integrals:\n\n");
407
print_mat(DSDBY,BasisSet.num_ao,BasisSet.num_ao,outfile);
408
fprintf(outfile," -dS/dB_z AO integrals:\n\n");
409
print_mat(DSDBZ,BasisSet.num_ao,BasisSet.num_ao,outfile);
410
fprintf(outfile," -dh/dB_x AO integrals:\n\n");
411
print_mat(DHDBX,BasisSet.num_ao,BasisSet.num_ao,outfile);
412
fprintf(outfile," -dh/dB_y AO integrals:\n\n");
413
print_mat(DHDBY,BasisSet.num_ao,BasisSet.num_ao,outfile);
414
fprintf(outfile," -dh/dB_z AO integrals:\n\n");
415
print_mat(DHDBZ,BasisSet.num_ao,BasisSet.num_ao,outfile);
416
/* fprintf(outfile," -d2h/dB_x dE_x AO integrals:\n\n");
417
print_mat(D2HDBXDEX,BasisSet.num_ao,BasisSet.num_ao,outfile);
418
fprintf(outfile," -d2h/dB_x dE_y AO integrals:\n\n");
419
print_mat(D2HDBXDEY,BasisSet.num_ao,BasisSet.num_ao,outfile);
420
fprintf(outfile," -d2h/dB_x dE_z AO integrals:\n\n");
421
print_mat(D2HDBXDEZ,BasisSet.num_ao,BasisSet.num_ao,outfile);
422
fprintf(outfile," -d2h/dB_y dE_x AO integrals:\n\n");
423
print_mat(D2HDBYDEX,BasisSet.num_ao,BasisSet.num_ao,outfile);
424
fprintf(outfile," -d2h/dB_y dE_y AO integrals:\n\n");
425
print_mat(D2HDBYDEY,BasisSet.num_ao,BasisSet.num_ao,outfile);
426
fprintf(outfile," -d2h/dB_y dE_z AO integrals:\n\n");
427
print_mat(D2HDBYDEZ,BasisSet.num_ao,BasisSet.num_ao,outfile);
428
fprintf(outfile," -d2h/dB_z dE_x AO integrals:\n\n");
429
print_mat(D2HDBZDEX,BasisSet.num_ao,BasisSet.num_ao,outfile);
430
fprintf(outfile," -d2h/dB_z dE_y AO integrals:\n\n");
431
print_mat(D2HDBZDEY,BasisSet.num_ao,BasisSet.num_ao,outfile);
432
fprintf(outfile," -d2h/dB_z dE_z AO integrals:\n\n");
433
print_mat(D2HDBZDEZ,BasisSet.num_ao,BasisSet.num_ao,outfile);*/
434
fprintf(outfile,"\n");
437
free_block(d2HdBxdExtemp);
438
free_block(d2HdBxdEytemp);
439
free_block(d2HdBxdEztemp);
440
free_block(d2HdBydExtemp);
441
free_block(d2HdBydEytemp);
442
free_block(d2HdBydEztemp);
443
free_block(d2HdBzdExtemp);
444
free_block(d2HdBzdEytemp);
445
free_block(d2HdBzdEztemp);
446
free_block(dHdBxtemp);
447
free_block(dHdBytemp);
448
free_block(dHdBztemp);
449
free_block(dSdBxtemp);
450
free_block(dSdBytemp);
451
free_block(dSdBztemp);
453
free_block(D2HDBXDEX);
454
free_block(D2HDBXDEY);
455
free_block(D2HDBXDEZ);
456
free_block(D2HDBYDEX);
457
free_block(D2HDBYDEY);
458
free_block(D2HDBYDEZ);
459
free_block(D2HDBZDEX);
460
free_block(D2HDBZDEY);
461
free_block(D2HDBZDEZ);
470
free_Taylor_Fm_Eval();