1 /*----------------------------------------------------------------------
2 eri.c
3
4 High-level interface of LIBERI
5
6 Coded by TOYODA Masayuki, 19 Jun. 2009.
7 ----------------------------------------------------------------------*/
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include <assert.h>
12 #include <time.h>
13 #include "eri_def.h"
14 #include "eri.h"
15
16
17
18
coord_spherical(double x,double y,double z,double * r,double * theta,double * phi)19 static void coord_spherical(
20 double x, /* (IN) Cartesian coordinate */
21 double y,
22 double z,
23 double *r, /* (OUT) spherical coordinate */
24 double *theta,
25 double *phi
26 )
27 {
28 *r = sqrt(x*x + y*y + z*z);
29 *theta = atan2(sqrt(x*x+y*y),z);
30 *phi = atan2(y,x);
31 }
32
33
34
35
ERI_Overlap(ERI_t * solver,double * glF,double * dglF[3],const ERI_Orbital_t * orb1,const ERI_Orbital_t * orb2,double cx)36 void ERI_Overlap(
37 ERI_t *solver, /* (INOUT) ERI object */
38 double *glF, /* (OUT) F term */
39 double *dglF[3], /* (OUT) derivatives of F */
40 const ERI_Orbital_t *orb1, /* (IN) orbital 1 */
41 const ERI_Orbital_t *orb2, /* (IN) orbital 2 */
42 double cx /* (IN) expansion center */
43 )
44 {
45 int ngrid, lmax, jmax, j;
46
47 const double *fr1, *fr2, *xr1, *xr2;
48 int ngrid1, ngrid2, l1, l2, m1, m2;
49 double c[3], c1[3], c2[3];
50 double ar1, at1, ap1, ar2, at2, ap2;
51
52 double *fk1, *fk2, *gam1, *gam2, *alp1, *alp2, *p, *F;
53 double *dgam1, *dgam2, *dalp1[3], *dalp2[3], *dp[3], *dF[3];
54
55 ngrid = ERI_ngrid(solver);
56 lmax = ERI_lmax(solver);
57 jmax = lmax*lmax;
58
59 /* WF 1 */
60 fr1 = orb1->fr;
61 xr1 = orb1->xr;
62 ngrid1 = orb1->ngrid;
63 l1 = orb1->l;
64 m1 = orb1->m;
65 c1[0] = orb1->c[0];
66 c1[1] = orb1->c[1];
67 c1[2] = orb1->c[2];
68
69 /* WF 2 */
70 fr2 = orb2->fr;
71 xr2 = orb2->xr;
72 ngrid2 = orb2->ngrid;
73 l2 = orb2->l;
74 m2 = orb2->m;
75 c2[0] = orb2->c[0];
76 c2[1] = orb2->c[1];
77 c2[2] = orb2->c[2];
78
79 STEPTRACE("ERI_Overlap: in\n");
80
81 /* expansion centers */
82 c[0] = cx*c1[0] + (1.0-cx)*c2[0];
83 c[1] = cx*c1[1] + (1.0-cx)*c2[1];
84 c[2] = cx*c1[2] + (1.0-cx)*c2[2];
85
86 coord_spherical(c1[0]-c[0], c1[1]-c[1], c1[2]-c[2], &ar1, &at1, &ap1);
87 coord_spherical(c2[0]-c[0], c2[1]-c[1], c2[2]-c[2], &ar2, &at2, &ap2);
88
89 /* allocate workspace memory */
90 fk1 = (double*)malloc( ERI_Size_of_Orbital(solver) );
91 fk2 = (double*)malloc( ERI_Size_of_Orbital(solver) );
92 gam1 = (double*)malloc( ERI_Size_of_Gamma(solver) );
93 gam2 = (double*)malloc( ERI_Size_of_Gamma(solver) );
94 alp1 = (double*)malloc( ERI_Size_of_Alpha(solver) );
95 alp2 = (double*)malloc( ERI_Size_of_Alpha(solver) );
96 p = (double*)malloc( ERI_Size_of_Overlap(solver) );
97 F = (double*)malloc( ERI_Size_of_Overlap(solver) );
98
99 if (dglF != NULL) {
100 dgam1 = (double*)malloc( ERI_Size_of_Gamma(solver) );
101 dgam2 = (double*)malloc( ERI_Size_of_Gamma(solver) );
102 dalp1[0] = (double*)malloc( ERI_Size_of_Alpha(solver) );
103 dalp1[1] = (double*)malloc( ERI_Size_of_Alpha(solver) );
104 dalp1[2] = (double*)malloc( ERI_Size_of_Alpha(solver) );
105 dalp2[0] = (double*)malloc( ERI_Size_of_Alpha(solver) );
106 dalp2[1] = (double*)malloc( ERI_Size_of_Alpha(solver) );
107 dalp2[2] = (double*)malloc( ERI_Size_of_Alpha(solver) );
108 dp[0] = (double*)malloc( ERI_Size_of_Overlap(solver) );
109 dp[1] = (double*)malloc( ERI_Size_of_Overlap(solver) );
110 dp[2] = (double*)malloc( ERI_Size_of_Overlap(solver) );
111 dF[0] = (double*)malloc( ERI_Size_of_Overlap(solver) );
112 dF[1] = (double*)malloc( ERI_Size_of_Overlap(solver) );
113 dF[2] = (double*)malloc( ERI_Size_of_Overlap(solver) );
114 }
115
116 ERI_Transform_Orbital(solver, fk1, fr1, xr1, ngrid1, l1);
117 ERI_Transform_Orbital(solver, fk2, fr2, xr2, ngrid2, l2);
118
119 if (NULL == dglF) {
120 /* without derivatives */
121 ERI_LL_Gamma(solver, gam1, NULL, fk1, ar1);
122 ERI_LL_Gamma(solver, gam2, NULL, fk2, ar2);
123 ERI_LL_Alpha(solver, alp1, gam1, at1, ap1, l1, m1);
124 ERI_LL_Alpha(solver, alp2, gam2, at2, ap2, l2, m2);
125 ERI_LL_Overlap(solver, p, alp1, alp2);
126 ERI_Transform_Overlap(solver, F, p );
127 ERI_GL_Interpolate(solver, glF, F);
128 } else {
129 ERI_LL_Gamma(solver, gam1, dgam1, fk1, ar1);
130 ERI_LL_Gamma(solver, gam2, dgam2, fk2, ar2);
131 ERI_LL_Alpha_d(solver, alp1, dalp1, gam1, dgam1, ar1, at1, ap1, l1, m1);
132 ERI_LL_Alpha_d(solver, alp2, dalp2, gam2, dgam2, ar2, at2, ap2, l2, m2);
133 ERI_LL_Overlap_d(solver, p, dp, alp1, dalp1, alp2, dalp2, cx);
134 ERI_Transform_Overlap(solver, F, p );
135 ERI_Transform_Overlap(solver, dF[0], dp[0] );
136 ERI_Transform_Overlap(solver, dF[1], dp[1] );
137 ERI_Transform_Overlap(solver, dF[2], dp[2] );
138 ERI_GL_Interpolate(solver, glF, F);
139 ERI_GL_Interpolate(solver, dglF[0], dF[0]);
140 ERI_GL_Interpolate(solver, dglF[1], dF[1]);
141 ERI_GL_Interpolate(solver, dglF[2], dF[2]);
142 }
143
144 /* release workspace memory */
145 free(fk1);
146 free(fk2);
147 free(gam1);
148 free(gam2);
149 free(alp1);
150 free(alp2);
151 free(p);
152 free(F);
153
154 if (NULL != dglF) {
155 free(dgam1);
156 free(dgam2);
157 free(dalp1[0]);
158 free(dalp1[1]);
159 free(dalp1[2]);
160 free(dalp2[0]);
161 free(dalp2[1]);
162 free(dalp2[2]);
163 free(dp[0]);
164 free(dp[1]);
165 free(dp[2]);
166 free(dF[0]);
167 free(dF[1]);
168 free(dF[2]);
169 }
170
171 STEPTRACE("ERI_Overlap: out\n");
172 }
173
174
175
ERI_Integral(ERI_t * solver,double I4[2],double * dI4,const ERI_Orbital_t * orb1,const ERI_Orbital_t * orb2,const ERI_Orbital_t * orb3,const ERI_Orbital_t * orb4,double scr)176 void ERI_Integral(
177 ERI_t *solver,
178 double I4[2], /* (OUT) integral */
179 double *dI4, /* (OUT) derivatives */
180 const ERI_Orbital_t *orb1,
181 const ERI_Orbital_t *orb2,
182 const ERI_Orbital_t *orb3,
183 const ERI_Orbital_t *orb4,
184 double scr
185 )
186 {
187 int lmax_gl;
188 double *glF, *glG, *dglF[3], *dglG[3];
189 double cx12, cx34, R[3], a1[3], a2[3], a3[3], a4[3];
190 double c1[3], c2[3], c3[3], c4[3];
191 clock_t clk1, clk2, clk3;
192
193 lmax_gl = ERI_lmax_gl(solver);
194
195 STEPTRACE("ERI_Integral: in");
196
197 /* allocate workspace memory */
198 glF = (double*)malloc( ERI_Size_of_GLF(solver) );
199 glG = (double*)malloc( ERI_Size_of_GLF(solver) );
200 if (dI4) {
201 dglF[0] = (double*)malloc( ERI_Size_of_GLF(solver) );
202 dglF[1] = (double*)malloc( ERI_Size_of_GLF(solver) );
203 dglF[2] = (double*)malloc( ERI_Size_of_GLF(solver) );
204 dglG[0] = (double*)malloc( ERI_Size_of_GLF(solver) );
205 dglG[1] = (double*)malloc( ERI_Size_of_GLF(solver) );
206 dglG[2] = (double*)malloc( ERI_Size_of_GLF(solver) );
207 }
208
209 c1[0] = orb1->c[0]; c1[1] = orb1->c[1]; c1[2] = orb1->c[2];
210 c2[0] = orb2->c[0]; c2[1] = orb2->c[1]; c2[2] = orb2->c[2];
211 c3[0] = orb3->c[0]; c3[1] = orb3->c[1]; c3[2] = orb3->c[2];
212 c4[0] = orb4->c[0]; c4[1] = orb4->c[1]; c4[2] = orb4->c[2];
213
214 /* expansion center */
215 cx12 = ERI_Center_r2(orb1, orb2);
216 cx34 = ERI_Center_r2(orb3, orb4);
217
218 /* R */
219 ERI_Coordinate_Transform(R, a1, a2, a3, a4, c1, c2, c3, c4, cx12, cx34);
220
221 if (NULL==dI4) {
222 ERI_Overlap(solver, glF, NULL, orb1, orb2, cx12);
223 ERI_Overlap(solver, glG, NULL, orb3, orb4, cx34);
224 ERI_Integral_GL(solver, I4, glF, glG, R[0], R[1], R[2], scr, lmax_gl);
225 } else {
226 ERI_Overlap(solver, glF, dglF, orb1, orb2, cx12);
227 ERI_Overlap(solver, glG, dglG, orb3, orb4, cx34);
228 ERI_Integral_GL_d(solver, I4, dI4, glF, glG, dglF, dglG,
229 R[0], R[1], R[2], cx12, cx34, 1e-10, scr, lmax_gl);
230 }
231
232 /* release workspace memory */
233 free(glF);
234 free(glG);
235 if (dI4) {
236 free(dglF[0]);
237 free(dglF[1]);
238 free(dglF[2]);
239 free(dglG[0]);
240 free(dglG[1]);
241 free(dglG[2]);
242 }
243
244 STEPTRACE("ERI_Integral: out");
245 }
246
247
248
249
250 /*----------------------------------------------------------------------
251 orbital_T
252
253 kinetic energy
254 ----------------------------------------------------------------------*/
orbital_T(ERI_t * solver,const double * fr,const double * xr,int nrgrid,int l)255 static double orbital_T(
256 ERI_t *solver,
257 const double *fr,
258 const double *xr,
259 int nrgrid,
260 int l
261 )
262 {
263 int i;
264 double k, af, ke, dk;
265 double *fk;
266
267 int nkgrid;
268 const double *kmesh;
269
270 nkgrid = ERI_ngrid(solver);
271 kmesh = ERI_Mesh_Array_k(solver);
272
273 fk = (double*)malloc( ERI_Size_of_Orbital(solver) );
274
275 ERI_Transform_Orbital(solver, fk, fr, xr, nrgrid, l);
276
277 ke = 0.0;
278 for (i=0; i<nkgrid; i++) {
279 k = ERI_Mesh_k(solver, i);
280 dk = ERI_Mesh_dk(solver, i);
281 af = fk[2*i+0]*fk[2*i+0]+fk[2*i+1]*fk[2*i+1];
282 ke += 8.0*(k*k*k*k)*af*dk;
283 }
284
285 free(fk);
286
287 return ke;
288 }
289
290
291
292
293 /*----------------------------------------------------------------------
294 orbital_r2
295
296 <r^2>
297 ----------------------------------------------------------------------*/
orbital_r2(const double * fr,const double * xr,int ngrid)298 static double orbital_r2(
299 const double *fr,
300 const double *xr,
301 int ngrid
302 )
303 {
304 int ir;
305 double norm, br2k, x0, x1, r2, r, f;
306
307 norm = 0.0;
308 br2k = 0.0;
309 for (ir=0; ir<ngrid-1; ir++) {
310 x0 = xr[ir];
311 x1 = xr[ir+1];
312 r = 0.5*(x0+x1);
313 r2 = r*r;
314 f = 0.5*(fr[ir]+fr[ir+1]);
315 norm += r2*f*f*(x1-x0);
316 br2k += r2*r2*f*f*(x1-x0);
317 }
318
319 return br2k/norm;
320 }
321
322
323
324
325 /*----------------------------------------------------------------------
326 orbital_r
327
328 <r>
329 ----------------------------------------------------------------------*/
orbital_r(const double * fr,const double * xr,int ngrid)330 static double orbital_r(
331 const double *fr,
332 const double *xr,
333 int ngrid
334 )
335 {
336 int ir;
337 double norm, brk, x0, x1, r, r2, f;
338
339 norm = 0.0;
340 brk = 0.0;
341 for (ir=0; ir<ngrid-1; ir++) {
342 x0 = xr[ir];
343 x1 = xr[ir+1];
344 r = 0.5*(x0+x1);
345 r2 = r*r;
346 f = 0.5*(fr[ir]+fr[ir+1]);
347 norm += r2*f*f*(x1-x0);
348 brk += r2*r*f*f*(x1-x0);
349 }
350
351 return brk/norm;
352 }
353
354
355
356
ERI_Coordinate_Transform(double out_R[3],double out_a1[3],double out_a2[3],double out_a3[3],double out_a4[3],const double in_a1[3],const double in_a2[3],const double in_a3[3],const double in_a4[3],double x12,double x34)357 void ERI_Coordinate_Transform(
358 double out_R[3], /* (OUT) displacement in spherical coord. */
359 double out_a1[3], /* (OUT) translated center 1 in spherical coord. */
360 double out_a2[3], /* (OUT) translated center 2 in spherical coord. */
361 double out_a3[3], /* (OUT) translated center 3 in spherical coord. */
362 double out_a4[3], /* (OUT) translated center 4 in spherical coord. */
363 const double in_a1[3], /* (IN) center 1 in Cartesian coord. */
364 const double in_a2[3], /* (IN) center 2 in Cartesian coord. */
365 const double in_a3[3], /* (IN) center 3 in Cartesian coord. */
366 const double in_a4[3], /* (IN) center 4 in Cartesian coord. */
367 double x12, /* (IN) ratio */
368 double x34 /* (IN) ratio */
369 )
370 {
371 double c12[3], c34[3];
372 double a1[3], a2[3], a3[3], a4[3];
373
374 c12[0] = x12*in_a1[0] + (1.0-x12)*in_a2[0];
375 c12[1] = x12*in_a1[1] + (1.0-x12)*in_a2[1];
376 c12[2] = x12*in_a1[2] + (1.0-x12)*in_a2[2];
377 c34[0] = x34*in_a3[0] + (1.0-x34)*in_a4[0];
378 c34[1] = x34*in_a3[1] + (1.0-x34)*in_a4[1];
379 c34[2] = x34*in_a3[2] + (1.0-x34)*in_a4[2];
380
381 a1[0] = in_a1[0]-c12[0]; a1[1] = in_a1[1]-c12[1]; a1[2] = in_a1[2]-c12[2];
382 a2[0] = in_a2[0]-c12[0]; a2[1] = in_a2[1]-c12[1]; a2[2] = in_a2[2]-c12[2];
383 a3[0] = in_a3[0]-c34[0]; a3[1] = in_a3[1]-c34[1]; a3[2] = in_a3[2]-c34[2];
384 a4[0] = in_a4[0]-c34[0]; a4[1] = in_a4[1]-c34[1]; a4[2] = in_a4[2]-c34[2];
385
386 coord_spherical(a1[0], a1[1], a1[2], &out_a1[0], &out_a1[1], &out_a1[2]);
387 coord_spherical(a2[0], a2[1], a2[2], &out_a2[0], &out_a2[1], &out_a2[2]);
388 coord_spherical(a3[0], a3[1], a3[2], &out_a3[0], &out_a3[1], &out_a3[2]);
389 coord_spherical(a4[0], a4[1], a4[2], &out_a4[0], &out_a4[1], &out_a4[2]);
390
391 coord_spherical(c34[0]-c12[0], c34[1]-c12[1], c34[2]-c12[2],
392 &out_R[0], &out_R[1], &out_R[2]);
393 }
394
395
396
397
398 /*----------------------------------------------------------------------
399 ERI_Center_r2
400 ----------------------------------------------------------------------*/
ERI_Center_r2(const ERI_Orbital_t * pA,const ERI_Orbital_t * pB)401 double ERI_Center_r2(
402 const ERI_Orbital_t *pA,
403 const ERI_Orbital_t *pB
404 )
405 {
406 double d1 = orbital_r2(pA->fr, pA->xr, pA->ngrid);
407 double d2 = orbital_r2(pB->fr, pB->xr, pB->ngrid);
408
409 return d2/(d1+d2);
410 }
411
412
413
414
415 /*----------------------------------------------------------------------
416 ERI_Center_DK
417 ----------------------------------------------------------------------*/
ERI_Center_DK(ERI_t * solver,const ERI_Orbital_t * pA,const ERI_Orbital_t * pB)418 double ERI_Center_DK(
419 ERI_t *solver,
420 const ERI_Orbital_t *pA,
421 const ERI_Orbital_t *pB
422 )
423 {
424 double d1 = orbital_r2(pA->fr, pA->xr, pA->ngrid);
425 double d2 = orbital_r2(pB->fr, pB->xr, pB->ngrid);
426 double k1 = orbital_T(solver, pA->fr, pA->xr, pA->ngrid, pA->l);
427 double k2 = orbital_T(solver, pB->fr, pB->xr, pB->ngrid, pB->l);
428
429 return (d2*k1)/(d2*k1+d1*k2);
430 }
431
432
433 /* EOF */
434