1 /*----------------------------------------------------------------------
2 exx_step1.c
3 ----------------------------------------------------------------------*/
4 #include <stdio.h>
5 #include <stdlib.h>
6 #include <math.h>
7 #include <time.h>
8 #include "exx.h"
9 #include "exx_vector.h"
10 #include "exx_log.h"
11 #include "exx_step1.h"
12 #include "exx_file_overlap.h"
13 #include "eri.h"
14 #include "eri_sf.h"
15 #include "eri_def.h"
16
17
EXX_Step1(EXX_t * exx,int natom,int * atom_sp,int nspec,int * spec_nb,double * spec_rc,int ** spec_l,int ** spec_m,int ** spec_mul,double **** spec_fr,double ** spec_xr,int * spec_nmesh)18 int EXX_Step1(
19 EXX_t *exx,
20 int natom,
21 int *atom_sp, /* [natom] */
22 int nspec,
23 int *spec_nb, /* [nspec] */
24 double *spec_rc, /* [nspec] */
25 int **spec_l, /* [nspec][nbmax] */
26 int **spec_m, /* [nspec][nbmax] */
27 int **spec_mul, /* [nspec][nbmax] */
28 double ****spec_fr, /* [nspec][lmax][mulmax][nmeshmax] */
29 double **spec_xr, /* [nspec][nmeshmax] */
30 int *spec_nmesh /* [nspec] */
31 )
32 {
33 int i, j;
34 int iatom1, iatom2, icell, ispec1, ispec2, ib1, ib2, iop;
35 int nb1, nb2, nbmax, nshell;
36 int nop, nop_n, nop_m, nop_local, iop0, npp_local;
37 int lmax, lmax_gl, ngrid, ngl, ndglf, ndalp, jmax_gl;
38 int l, m, mul;
39 double x, y, z, c1[3], c2[3], rc1, rc2, d, rc12, cx12, c12[3];
40 double ar1, at1, ap1, ar2, at2, ap2;
41
42 const int *op_atom1, *op_atom2, *op_cell, *atom_nb;
43 double *fk, *gam, *alp1, *alp2, *ovlp, *ovlF;
44 double **glf; /* [nop_local][ndglf*nbmax*nbmax] */
45 const double *fr, *xr, *pvec, *atom_v;
46
47 ERI_Init_Misc info;
48 ERI_t *solver;
49
50 clock_t clk1, clk2;
51 int nstep;
52 double et;
53
54
55 int nq_full, nq_sig;
56
57 char path_cache[EXX_PATHLEN];
58
59 int myrank, nproc;
60
61 #ifdef EXX_USE_MPI
62 MPI_Comm comm;
63 MPI_Status mpistat;
64 #endif
65
66 #ifdef EXX_USE_MPI
67 comm = g_exx_mpicomm;
68 MPI_Comm_rank(comm, &myrank);
69 MPI_Comm_size(comm, &nproc);
70 #else
71 myrank = EXX_ROOT_RANK;
72 nproc = 1;
73 #endif
74
75 op_atom1 = EXX_Array_OP_Atom1(exx);
76 op_atom2 = EXX_Array_OP_Atom2(exx);
77 op_cell = EXX_Array_OP_Cell(exx);
78 atom_nb = EXX_atom_nb(exx);
79 pvec = EXX_pvec(exx);
80 atom_v = EXX_atom_v(exx);
81 nbmax = EXX_nbmax(exx);
82
83 nshell = EXX_Number_of_OP_Shells(exx);
84 nop = EXX_Number_of_OP(exx);
85 nop_n = nop / nproc;
86 nop_m = nop - nop_n*nproc;
87
88 EXX_LOG_TRACE_INTEGER("NOP", nop);
89 EXX_LOG_TRACE_INTEGER("NOP_N", nop_n);
90 EXX_LOG_TRACE_INTEGER("NOP_M", nop_m);
91
92 sprintf(path_cache, "%s/exxovlp.dat", EXX_CacheDir(exx));
93
94 #ifdef EXX_USE_MPI
95 MPI_Barrier(comm);
96 #endif
97
98 nop_local = (myrank<nop_m) ? (nop_n+1) : (nop_n);
99 iop0 = 0;
100 for (i=0; i<myrank; i++) {
101 iop0 += (i<nop_m) ? (nop_n+1) : (nop_n);
102 }
103
104 EXX_LOG_TRACE_INTEGER("IOP0" , iop0);
105 EXX_LOG_TRACE_INTEGER("NOP_N" , nop_n);
106 EXX_LOG_TRACE_INTEGER("NOP_M" , nop_m);
107 EXX_LOG_TRACE_INTEGER("NOP_LOCAL", nop_local);
108
109 /*----- initialize LIBERI -----*/
110 lmax = g_exx_liberi_lmax;
111 lmax_gl = (g_exx_liberi_lmax+1)/2;
112 ngrid = g_exx_liberi_ngrid;
113 ngl = g_exx_liberi_ngl;;
114
115 info.sbttype = ERI_SBT_LINEAR;
116 info.rmax = 80.0;
117
118 solver = ERI_Init(lmax, lmax_gl, ngrid, ngl, ERI_SH_REAL, &info);
119 if (NULL==solver) { EXX_ERROR( "failed to initialize LIBERI" ); }
120
121 #if 1
122 for (i=0; i<natom; i++) {
123 EXX_Vector_F2C(c1, &atom_v[3*i], pvec);
124 EXX_Log_Print(" atom %2d : ( %10.4f, %10.4f, %10.4f)\n",
125 i, c1[0], c1[1], c1[2]);
126 }
127 EXX_Log_Print("\n");
128 #endif
129
130 jmax_gl = lmax_gl*lmax_gl;
131
132 npp_local = 0;
133 for (i=0; i<nop_local; i++) {
134 iop = iop0 + i;
135 nb1 = atom_nb[op_atom1[iop]];
136 nb2 = atom_nb[op_atom2[iop]];
137 npp_local += nb1*nb2;
138 }
139
140 ndglf = ERI_Size_of_GLF(solver)/sizeof(double);
141 ndalp = ERI_Size_of_Alpha(solver)/sizeof(double);
142 fk = (double*)malloc( ERI_Size_of_Orbital(solver) );
143 gam = (double*)malloc( ERI_Size_of_Gamma(solver) );
144 alp1 = (double*)malloc( sizeof(double)*ndalp*nbmax );
145 alp2 = (double*)malloc( sizeof(double)*ndalp*nbmax );
146 ovlp = (double*)malloc( ERI_Size_of_Overlap(solver) );
147 ovlF = (double*)malloc( ERI_Size_of_Overlap(solver) );
148 glf = (double**)malloc( sizeof(double*)*nop_local );
149 for (i=0; i<nop_local; i++) {
150 glf[i] = (double*)malloc( sizeof(double)*ndglf*nbmax*nbmax );
151 }
152
153 clk1 = clock();
154
155 nstep = 0;
156
157 for (i=0; i<nop_local; i++) {
158 iop = iop0 + i;
159 iatom1 = op_atom1[iop];
160 iatom2 = op_atom2[iop];
161 icell = op_cell[iop];
162
163 /* atom 1 */
164 EXX_Vector_F2C(c1, &atom_v[3*iatom1], pvec);
165 ispec1 = atom_sp[iatom1];
166 rc1 = spec_rc[ispec1];
167 nb1 = spec_nb[ispec1];
168
169 /* atom 2 */
170 EXX_Vector_F2C_Offsite(c2, &atom_v[3*iatom2], pvec, icell, nshell);
171 ispec2 = atom_sp[iatom2];
172 rc2 = spec_rc[ispec2];
173 nb2 = spec_nb[ispec2];
174
175 d = EXX_Vector_Distance(c1, c2);
176 EXX_Vector_PAO_Overlap(rc1, rc2, d, &rc12, &cx12);
177 c12[0] = cx12*c1[0] + (1.0-cx12)*c2[0];
178 c12[1] = cx12*c1[1] + (1.0-cx12)*c2[1];
179 c12[2] = cx12*c1[2] + (1.0-cx12)*c2[2];
180
181 /* translation */
182 c1[0] -= c12[0]; c1[1] -= c12[1]; c1[2] -= c12[2];
183 c2[0] -= c12[0]; c2[1] -= c12[1]; c2[2] -= c12[2];
184
185 /* cartesian -> spherical */
186 EXX_Vector_C2S(c1, &ar1, &at1, &ap1);
187 EXX_Vector_C2S(c2, &ar2, &at2, &ap2);
188
189 /* pre-calculate alpha terms */
190 for (ib1=0; ib1<nb1; ib1++) {
191 l = spec_l[ispec1][ib1];
192 mul = spec_mul[ispec1][ib1];
193 m = spec_m[ispec1][ib1];
194 fr = spec_fr[ispec1][l][mul];
195 xr = spec_xr[ispec1];
196 ngrid = spec_nmesh[ispec1];
197
198 ERI_Transform_Orbital(solver, fk, fr, xr, ngrid, l);
199 ERI_LL_Gamma(solver, gam, NULL, fk, ar1);
200 ERI_LL_Alpha(solver, &alp1[ndalp*ib1], gam, at1, ap1, l, m);
201 }
202
203 for (ib2=0; ib2<nb2; ib2++) {
204 l = spec_l[ispec2][ib2];
205 mul = spec_mul[ispec2][ib2];
206 m = spec_m[ispec2][ib2];
207 fr = spec_fr[ispec2][l][mul];
208 xr = spec_xr[ispec2];
209 ngrid = spec_nmesh[ispec2];
210
211 ERI_Transform_Orbital(solver, fk, fr, xr, ngrid, l);
212 ERI_LL_Gamma(solver, gam, NULL, fk, ar2);
213 ERI_LL_Alpha(solver, &alp2[ndalp*ib2], gam, at2, ap2, l, m);
214 }
215
216 for (ib1=0; ib1<nb1; ib1++) {
217 for (ib2=0; ib2<nb2; ib2++) {
218 ERI_LL_Overlap(solver, ovlp, &alp1[ndalp*ib1], &alp2[ndalp*ib2]);
219 ERI_Transform_Overlap(solver, ovlF, ovlp);
220 ERI_GL_Interpolate(solver, &glf[i][(ib1*nbmax+ib2)*ndglf], ovlF);
221 nstep++;
222 }
223 }
224 } /* loop of iop */
225
226 clk2 = clock();
227
228 #ifdef EXX_USE_MPI
229 MPI_Barrier(comm);
230 #endif
231
232 /* save data */
233 EXX_File_Overlap_Write(ndglf, nop_local, nbmax, jmax_gl, glf, path_cache);
234
235 #ifdef EXX_USE_MPI
236 MPI_Barrier(comm);
237 #endif
238
239 free(fk);
240 free(gam);
241 free(alp1);
242 free(alp2);
243 free(ovlp);
244 free(ovlF);
245 for (i=0; i<nop_local; i++) { free( glf[i]); }
246 free(glf);
247
248 ERI_Free(solver);
249
250 et = (double)(clk2-clk1)/(double)CLOCKS_PER_SEC;
251
252 EXX_LOG_MESSAGE("REPORT OF STEP1:");
253 EXX_LOG_TRACE_INTEGER("NPROC", nproc);
254 EXX_LOG_TRACE_DOUBLE("TIME ", et);
255 EXX_LOG_TRACE_DOUBLE("AVG_TIME", et/(double)nstep);
256
257 return 0;
258 }
259
260
261