1 /*----------------------------------------------------------------------
2 demo.c
3
4 Sample program of LIBERI
5
6 Coded by TOYODA MAsayuki, June 2009.
7 ----------------------------------------------------------------------*/
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include <time.h>
12 #include "eri.h"
13
14 #define NMESH 8192
15
16
integral(ERI_t * solver,ERI_Orbital_t p[4],int nitr,double I4[3][2],double * sec)17 static void integral(
18 ERI_t *solver,
19 ERI_Orbital_t p[4],
20 int nitr,
21 double I4[3][2], /* (OUT) integral */
22 double *sec
23 )
24 {
25 int lmax_gl, itr;
26 double *glF, *glG;
27 double cx12, cx34, R[3], a1[3], a2[3], a3[3], a4[3];
28 double c1[3], c2[3], c3[3], c4[3];
29 double Rd[3], scr;
30 clock_t t1, t2;
31
32 lmax_gl = ERI_lmax_gl(solver);
33
34 /* allocate workspace memory */
35 glF = (double*)malloc( ERI_Size_of_GLF(solver) );
36 glG = (double*)malloc( ERI_Size_of_GLF(solver) );
37
38 c1[0] = p[0].c[0]; c1[1] = p[0].c[1]; c1[2] = p[0].c[2];
39 c2[0] = p[1].c[0]; c2[1] = p[1].c[1]; c2[2] = p[1].c[2];
40 c3[0] = p[2].c[0]; c3[1] = p[2].c[1]; c3[2] = p[2].c[2];
41 c4[0] = p[3].c[0]; c4[1] = p[3].c[1]; c4[2] = p[3].c[2];
42
43 /* expansion center */
44 cx12 = ERI_Center_r2(&p[0], &p[1]);
45 cx34 = ERI_Center_r2(&p[2], &p[3]);
46
47 scr = ERI_NOSCREEN;
48
49 ERI_Overlap(solver, glF, NULL, &p[0], &p[1], cx12);
50 ERI_Overlap(solver, glG, NULL, &p[2], &p[3], cx34);
51
52 /* R */
53 ERI_Coordinate_Transform(R, a1, a2, a3, a4, c1, c2, c3, c4, cx12, cx34);
54
55 t1 = clock();
56 for (itr=0; itr<nitr; itr++) {
57 Rd[0] = R[0];
58 Rd[1] = R[1];
59 Rd[2] = R[2];
60 ERI_Integral_GL(solver, I4[0], glF, glG,
61 Rd[0], Rd[1], Rd[2], scr, lmax_gl);
62
63 Rd[0] = R[0]+1.0;
64 Rd[1] = R[1];
65 Rd[2] = R[2];
66 ERI_Integral_GL(solver, I4[1], glF, glG,
67 Rd[0], Rd[1], Rd[2], scr, lmax_gl);
68
69 Rd[0] = R[0];
70 Rd[1] = R[1]+0.4;
71 Rd[2] = R[2];
72 ERI_Integral_GL(solver, I4[2], glF, glG,
73 Rd[0], Rd[1], Rd[2], scr, lmax_gl);
74 }
75 t2 = clock();
76
77 *sec = (double)(t2-t1)/(double)CLOCKS_PER_SEC;
78
79 /* release workspace memory */
80 free(glF);
81 free(glG);
82 }
83
84
85
86
integral_pp(ERI_t * solver,ERI_Orbital_t p[4],int nitr,double I4[3][2],double * sec)87 static void integral_pp(
88 ERI_t *solver,
89 ERI_Orbital_t p[4],
90 int nitr,
91 double I4[3][2],
92 double *sec
93 )
94 {
95 int lmax_gl, ngl, itr, jmax_gl;
96 double *glF, *glG;
97 double cx12, cx34, R[3], a1[3], a2[3], a3[3], a4[3];
98 double c1[3], c2[3], c3[3], c4[3], scr, Rd[3], Rd_t[3], Rd_p[3];
99 double tmp_I4[6];
100 clock_t t1, t2;
101
102 int numR;
103 double *mul_gc, *prej, *preY;
104 int *mul_j2, mul_n;
105
106 int num_minimalR, *minimalR;
107
108 int i,j;
109
110 lmax_gl = ERI_lmax_gl(solver);
111 jmax_gl = lmax_gl*lmax_gl;
112 ngl = ERI_ngl(solver);
113 scr = ERI_NOSCREEN;
114 numR = 3;
115
116 /* allocate workspace memory */
117 glF = (double*)malloc( ERI_Size_of_GLF(solver) );
118 glG = (double*)malloc( ERI_Size_of_GLF(solver) );
119 mul_j2 = (int*)malloc(sizeof(int)*jmax_gl*jmax_gl*jmax_gl);
120 mul_gc = (double*)malloc(sizeof(double)*jmax_gl*jmax_gl*jmax_gl);
121 mul_n = (int*)malloc(sizeof(double)*jmax_gl*jmax_gl);
122 prej = (double*)malloc(sizeof(double)*ngl*lmax_gl*numR);
123 preY = (double*)malloc(sizeof(double)*numR*jmax_gl);
124 minimalR= (int*)malloc(sizeof(int)*numR);
125
126 c1[0] = p[0].c[0]; c1[1] = p[0].c[1]; c1[2] = p[0].c[2];
127 c2[0] = p[1].c[0]; c2[1] = p[1].c[1]; c2[2] = p[1].c[2];
128 c3[0] = p[2].c[0]; c3[1] = p[2].c[1]; c3[2] = p[2].c[2];
129 c4[0] = p[3].c[0]; c4[1] = p[3].c[1]; c4[2] = p[3].c[2];
130
131 /* expansion center */
132 cx12 = ERI_Center_r2(&p[0], &p[1]);
133 cx34 = ERI_Center_r2(&p[2], &p[3]);
134
135 /* R */
136 ERI_Coordinate_Transform(R, a1, a2, a3, a4, c1, c2, c3, c4, cx12, cx34);
137
138 ERI_Overlap(solver, glF, NULL, &p[0], &p[1], cx12);
139 ERI_Overlap(solver, glG, NULL, &p[2], &p[3], cx34);
140
141 Rd [0] = R[0];
142 Rd_t[0] = R[1];
143 Rd_p[0] = R[2];
144
145 Rd [1] = R[0]+1.0;
146 Rd_t[1] = R[1];
147 Rd_p[1] = R[2];
148
149 Rd [2] = R[0];
150 Rd_t[2] = R[1]+0.4;
151 Rd_p[2] = R[2];
152
153 ERI_Integral_GL_PrejY(solver, Rd, Rd_t, Rd_p, numR,
154 scr, prej, preY, mul_j2, mul_gc, mul_n, minimalR, &num_minimalR);
155
156 t1 = clock();
157 for (itr=0; itr<nitr; itr++) {
158 ERI_Integral_GL_Post(solver, tmp_I4, glF, glG, numR,
159 prej, preY, mul_j2, mul_gc, mul_n, minimalR, num_minimalR);
160 }
161 t2 = clock();
162
163 *sec = (double)(t2-t1)/(double)CLOCKS_PER_SEC;
164
165 I4[0][0] = tmp_I4[0];
166 I4[0][1] = 0.0;
167 I4[1][0] = tmp_I4[1];
168 I4[1][1] = 0.0;
169 I4[2][0] = tmp_I4[2];
170 I4[2][1] = 0.0;
171
172 /* release workspace memory */
173 free(glF);
174 free(glG);
175
176 free(mul_j2);
177 free(mul_gc);
178 free(mul_n);
179 free(prej);
180 free(preY);
181 free(minimalR);
182 }
183
184
185
186
main(void)187 int main(void)
188 {
189 int i;
190 double r, I4[3][2];
191 const double rmax = 300.0;
192 double s1, s2, s3;
193
194 ERI_t *eri;
195 ERI_Orbital_t p[4];
196
197 double xr[NMESH], gto_s[NMESH];
198
199 /* radial mesh and GTO */
200 for (i=0; i<NMESH; i++) {
201 r = rmax*(double)i/(double)NMESH;
202 xr[i] = r;
203 gto_s[i] = exp(-r*r);
204 }
205
206 /* prepare orbital informations */
207 p[0].fr = gto_s; p[0].xr = xr; p[0].ngrid = NMESH; p[0].l = 1; p[0].m = -1;
208 p[1].fr = gto_s; p[1].xr = xr; p[1].ngrid = NMESH; p[1].l = 0; p[1].m = 0;
209 p[2].fr = gto_s; p[2].xr = xr; p[2].ngrid = NMESH; p[2].l = 0; p[2].m = 0;
210 p[3].fr = gto_s; p[3].xr = xr; p[3].ngrid = NMESH; p[3].l = 0; p[3].m = 0;
211
212 /* positions */
213 p[0].c[0] = 0.2; p[0].c[1] = 0.2; p[0].c[2] = 0.0;
214 p[1].c[0] = 0.2; p[1].c[1] = -0.2; p[1].c[2] = 0.0;
215 p[2].c[0] = -0.2; p[2].c[1] = 0.3; p[2].c[2] = 0.0;
216 p[3].c[0] = -0.2; p[3].c[1] = -0.3; p[3].c[2] = 0.0;
217
218 /* initialize LIBERI */
219 //eri = ERI_Init(15, 8, 1024, 100, ERI_SH_COMPLEX, NULL);
220 eri = ERI_Init(15, 8, 1024, 100, ERI_SH_REAL, NULL);
221 if (NULL == eri) {
222 fprintf(stderr, "*** Failed to initialize LIBERI.\n");
223 return -1;
224 }
225
226 /* calculate ERI */
227 integral(eri, p, 100, I4, &s1);
228 printf(" %10.6f %10.6f\n", I4[0][0], I4[0][1]);
229 printf(" %10.6f %10.6f\n", I4[1][0], I4[1][1]);
230 printf(" %10.6f %10.6f\n", I4[2][0], I4[2][1]);
231 printf(" TIME = %10.6f\n", s1);
232 printf("\n");
233
234 integral_pp(eri, p, 100, I4, &s2);
235 printf(" %10.6f %10.6f\n", I4[0][0], I4[0][1]);
236 printf(" %10.6f %10.6f\n", I4[1][0], I4[1][1]);
237 printf(" %10.6f %10.6f\n", I4[2][0], I4[2][1]);
238 printf(" TIME = %10.6f\n", s2);
239 printf("\n");
240
241 /* release LIBERI */
242 ERI_Free(eri);
243
244 return 0;
245 }
246
247
248 /* EOF */
249