1 /*----------------------------------------------------------------------
2 exx.c
3
4 M. Toyoda, 10 Nov. 2009.
5 ----------------------------------------------------------------------*/
6 #include <stdio.h>
7 #include <stdlib.h>
8 #include <math.h>
9 #include <string.h>
10 #include <assert.h>
11 #include "exx.h"
12 #include "exx_log.h"
13 #include "exx_vector.h"
14 #include "exx_index.h"
15
16
17 struct EXX_Struct {
18 /* unit cell */
19 int natom; /* number of atoms */
20 double *atom_rc; /* confinement length of atoms */
21 double *atom_v; /* atomic positions */
22 int *atom_nb; /* number of basis */
23 double pvec[9]; /* primitive translational vectors */
24 int nbmax; /* max number of basis */
25
26 /* parameters */
27 double w_scr; /* screening parameter */
28 double rc_cut; /* truncation length */
29
30 /* overlaping pairs */
31 int nshell_op; /* number of shells */
32 int nop; /* number of pairs */
33 int *op_atom1; /* arrays for atom1 */
34 int *op_atom2; /* arrays for atom2 */
35 int *op_cell; /* arrays for cells */
36
37 /* exchanging pairs */
38 int nshell_ep; /* number of shells */
39 int nep; /* number of pairs */
40 int *ep_atom1; /* arrays for atom1 */
41 int *ep_atom2; /* arrays for atom2 */
42 int *ep_cell; /* arrays for cells */
43
44 /* path */
45 char path_cachedir[EXX_PATHLEN];
46 };
47
48
49
50
51 int g_exx_skip1 = 0;
52 int g_exx_skip2 = 0;
53 char g_exx_cachedir[EXX_PATHLEN];
54 int g_exx_liberi_lmax = 20;
55 int g_exx_liberi_ngrid = 1024;
56 int g_exx_liberi_ngl = 100;
57 double g_exx_rc_cut = -1.0;
58 double g_exx_w_scr = -1.0;
59
60
61
62 /*----------------------------------------------------------------------
63 EXX_New
64
65 IN:
66 natom : number of atoms
67 atom_rc : confinement length of atoms
68 atom_v : atomic positions
69 pvec : primitive translational vectors
70 w_scr : screening parameter
71 rc_cut : truncation length
72
73 OUT:
74 return : pointer to EXX object
75 ----------------------------------------------------------------------*/
EXX_New(int natom,const double * atom_v,const int * atom_sp,int nspec,const double * spec_rc,const int * spec_nb,const double * pvec,double w_scr,double rc_cut,int system_type,const char * cachedir)76 EXX_t* EXX_New(
77 int natom,
78 const double *atom_v,
79 const int *atom_sp,
80 int nspec,
81 const double *spec_rc,
82 const int *spec_nb,
83 const double *pvec,
84 double w_scr,
85 double rc_cut,
86 int system_type,
87 const char *cachedir
88 )
89 {
90 int i, j, nshell_op, nshell_ep, nop, nep, nbmax;
91 EXX_t *self;
92
93 self = (EXX_t*)malloc(sizeof(EXX_t));
94
95 /* max number of basis */
96 nbmax = 0;
97 for (i=0; i<nspec; i++) {
98 if (spec_nb[i] > nbmax) { nbmax = spec_nb[i]; }
99 }
100
101 /* allocate memory */
102 self->atom_rc = (double*)malloc(sizeof(double)*natom);
103 self->atom_v = (double*)malloc(sizeof(double)*natom*3);
104 self->atom_nb = (int*)malloc(sizeof(int)*natom);
105
106 /* copy data */
107 self->natom = natom;
108 self->nbmax = nbmax;
109 for (i=0; i<natom; i++) {
110 self->atom_v[3*i+0] = atom_v[3*i+0];
111 self->atom_v[3*i+1] = atom_v[3*i+1];
112 self->atom_v[3*i+2] = atom_v[3*i+2];
113 j = atom_sp[i];
114 self->atom_rc[i] = spec_rc[j];
115 self->atom_nb[i] = spec_nb[j];
116 }
117
118 /* number of shells for OP and EP */
119 switch (system_type) {
120 case EXX_SYSTEM_PERIODIC:
121 nshell_op =
122 EXX_Index_OP_NShell(natom, self->atom_rc, self->atom_v, pvec);
123 nshell_ep =
124 EXX_Index_EP_NShell(natom, self->atom_rc, self->atom_v, pvec, rc_cut);
125 break;
126 case EXX_SYSTEM_CLUSTER:
127 nshell_op = 0;
128 nshell_ep = 0;
129 break;
130 default:
131 EXX_ERROR( "undefined system_type" );
132 break;
133 }
134
135 /* number of pais */
136 nop = EXX_Index_OP(natom, nshell_op, self->atom_rc, self->atom_v,
137 pvec, NULL, NULL, NULL);
138 nep = EXX_Index_EP(natom, nshell_ep, self->atom_rc, self->atom_v,
139 pvec, rc_cut, NULL, NULL, NULL);
140
141 /* allocate memory */
142 self->op_atom1 = (int*)malloc(sizeof(int)*nop);
143 self->op_atom2 = (int*)malloc(sizeof(int)*nop);
144 self->op_cell = (int*)malloc(sizeof(int)*nop);
145 self->ep_atom1 = (int*)malloc(sizeof(int)*nep);
146 self->ep_atom2 = (int*)malloc(sizeof(int)*nep);
147 self->ep_cell = (int*)malloc(sizeof(int)*nep);
148
149 /* enumurate pairs */
150 EXX_Index_OP(natom, nshell_op, self->atom_rc, self->atom_v, pvec,
151 self->op_atom1, self->op_atom2, self->op_cell);
152 EXX_Index_EP(natom, nshell_ep, self->atom_rc, self->atom_v, pvec, rc_cut,
153 self->ep_atom1, self->ep_atom2, self->ep_cell);
154
155 /* copy data */
156 for (i=0; i<9; i++) { self->pvec[i] = pvec[i]; }
157 self->w_scr = w_scr;
158 self->rc_cut = rc_cut;
159 self->nshell_op = nshell_op;
160 self->nop = nop;
161 self->nshell_ep = nshell_ep;
162 self->nep = nep;
163
164 strncpy(self->path_cachedir, cachedir, EXX_PATHLEN);
165
166 return self;
167 }
168
169
170
171
EXX_Free(EXX_t * self)172 void EXX_Free(EXX_t * self)
173 {
174 if (self) {
175 if (self->atom_rc) { free(self->atom_rc); }
176 if (self->atom_v) { free(self->atom_v); }
177 if (self->atom_nb) { free(self->atom_nb); }
178 if (self->op_atom1) { free(self->op_atom1); }
179 if (self->op_atom2) { free(self->op_atom2); }
180 if (self->op_cell) { free(self->op_cell); }
181 if (self->ep_atom1) { free(self->ep_atom1); }
182 if (self->ep_atom2) { free(self->ep_atom2); }
183 if (self->ep_cell) { free(self->ep_cell); }
184 free(self);
185 }
186 }
187
188
189
EXX_natom(const EXX_t * self)190 int EXX_natom(const EXX_t *self)
191 {
192 return self->natom;
193 }
194
EXX_atom_rc(const EXX_t * self)195 const double* EXX_atom_rc(const EXX_t *self)
196 {
197 return self->atom_rc;
198 }
199
EXX_atom_v(const EXX_t * self)200 const double* EXX_atom_v(const EXX_t *self)
201 {
202 return self->atom_v;
203 }
204
EXX_atom_nb(const EXX_t * self)205 const int* EXX_atom_nb(const EXX_t *self)
206 {
207 return self->atom_nb;
208 }
209
EXX_pvec(const EXX_t * self)210 const double* EXX_pvec(const EXX_t *self)
211 {
212 return self->pvec;
213 }
214
EXX_w_scr(const EXX_t * self)215 double EXX_w_scr(const EXX_t *self)
216 {
217 return self->w_scr;
218 }
219
EXX_rc_cut(const EXX_t * self)220 double EXX_rc_cut(const EXX_t *self)
221 {
222 return self->rc_cut;
223 }
224
EXX_nbmax(const EXX_t * self)225 int EXX_nbmax(const EXX_t *self)
226 {
227 return self->nbmax;
228 }
229
EXX_Number_of_OP_Shells(const EXX_t * self)230 int EXX_Number_of_OP_Shells(const EXX_t *self)
231 {
232 return self->nshell_op;
233 }
234
EXX_Number_of_OP(const EXX_t * self)235 int EXX_Number_of_OP(const EXX_t *self)
236 {
237 return self->nop;
238 }
239
EXX_Array_OP_Atom1(const EXX_t * self)240 const int* EXX_Array_OP_Atom1(const EXX_t *self)
241 {
242 return self->op_atom1;
243 }
244
EXX_Array_OP_Atom2(const EXX_t * self)245 const int* EXX_Array_OP_Atom2(const EXX_t *self)
246 {
247 return self->op_atom2;
248 }
249
EXX_Array_OP_Cell(const EXX_t * self)250 const int* EXX_Array_OP_Cell(const EXX_t *self)
251 {
252 return self->op_cell;
253 }
254
EXX_Number_of_EP_Shells(const EXX_t * self)255 int EXX_Number_of_EP_Shells(const EXX_t *self)
256 {
257 return self->nshell_ep;
258 }
259
EXX_Number_of_EP(const EXX_t * self)260 int EXX_Number_of_EP(const EXX_t *self)
261 {
262 return self->nep;
263 }
264
EXX_Array_EP_Atom1(const EXX_t * self)265 const int* EXX_Array_EP_Atom1(const EXX_t *self)
266 {
267 return self->ep_atom1;
268 }
269
EXX_Array_EP_Atom2(const EXX_t * self)270 const int* EXX_Array_EP_Atom2(const EXX_t *self)
271 {
272 return self->ep_atom2;
273 }
274
EXX_Array_EP_Cell(const EXX_t * self)275 const int* EXX_Array_EP_Cell(const EXX_t *self)
276 {
277 return self->ep_cell;
278 }
279
EXX_CacheDir(const EXX_t * self)280 const char* EXX_CacheDir(const EXX_t *self)
281 {
282 return self->path_cachedir;
283 }
284
285
286
287 /*----------------------------------------------------------------------
288 EXX_Find_OP
289
290 Find OP index whose iatom1, iatom2, and icell are equivalent with
291 those specified. If index is not found, this returns -1.
292 ----------------------------------------------------------------------*/
EXX_Find_OP(const EXX_t * self,int iatom1,int iatom2,int iR_x,int iR_y,int iR_z)293 int EXX_Find_OP(
294 const EXX_t *self,
295 int iatom1,
296 int iatom2,
297 int iR_x,
298 int iR_y,
299 int iR_z
300 )
301 {
302 int i, nop, nshell, ncd, iR_max, icell;
303 const int *op_atom1, *op_atom2, *op_cell;
304
305 nop = EXX_Number_of_OP(self);
306 nshell = EXX_Number_of_OP_Shells(self);
307 ncd = 2*nshell+1;
308
309 op_atom1 = EXX_Array_OP_Atom1(self);
310 op_atom2 = EXX_Array_OP_Atom2(self);
311 op_cell = EXX_Array_OP_Cell(self);
312
313 iR_max = abs(iR_x);
314 if (abs(iR_y)>iR_max) { iR_max = abs(iR_y); }
315 if (abs(iR_z)>iR_max) { iR_max = abs(iR_z); }
316 if (iR_max>nshell) { return -1; }
317
318 icell = (iR_x+nshell) + (iR_y+nshell)*ncd + (iR_z+nshell)*ncd*ncd;
319
320 for (i=0; i<nop; i++) {
321 if (iatom1==op_atom1[i] && iatom2==op_atom2[i] && icell==op_cell[i]) {
322 return i;
323 }
324 }
325
326 return -1;
327 }
328
329
330 /*----------------------------------------------------------------------
331 EXX_Find_EP
332
333 Find EP index whose iatom1, iatom2, and icell are equivalent with
334 those specified. If index is not found, this returns -1.
335 ----------------------------------------------------------------------*/
EXX_Find_EP(const EXX_t * self,int iatom1,int iatom2,int iR_x,int iR_y,int iR_z)336 int EXX_Find_EP(
337 const EXX_t *self,
338 int iatom1,
339 int iatom2,
340 int iR_x,
341 int iR_y,
342 int iR_z
343 )
344 {
345 int i, nep, nshell, ncd, iR_max, icell;
346 const int *ep_atom1, *ep_atom2, *ep_cell;
347
348 nep = EXX_Number_of_EP(self);
349 nshell = EXX_Number_of_EP_Shells(self);
350 ncd = 2*nshell+1;
351
352 ep_atom1 = EXX_Array_EP_Atom1(self);
353 ep_atom2 = EXX_Array_EP_Atom2(self);
354 ep_cell = EXX_Array_EP_Cell(self);
355
356 iR_max = abs(iR_x);
357 if (abs(iR_y)>iR_max) { iR_max = abs(iR_y); }
358 if (abs(iR_z)>iR_max) { iR_max = abs(iR_z); }
359 if (iR_max>nshell) { return -1; }
360
361 icell = (iR_x+nshell) + (iR_y+nshell)*ncd + (iR_z+nshell)*ncd*ncd;
362
363 for (i=0; i<nep; i++) {
364 if (iatom1==ep_atom1[i] && iatom2==ep_atom2[i] && icell==ep_cell[i]) {
365 return i;
366 }
367 }
368
369 return -1;
370 }
371
372
373
374 #if 0
375 /*---*/
376
377
378 static void quad(
379 const EXX_t *xfm,
380 int iatom1,
381 int iatom2,
382 int iatom3,
383 int iatom4,
384 int icell1,
385 int icell2,
386 int icell3,
387 int *w,
388 int q_ep1[8],
389 int q_ep2[8],
390 int q_cell[8],
391 int *mul
392 )
393 {
394 int nshell_op, ncd_op, ncell_op, ic0_op;
395 int nshell_ep, ncd_ep, ncell_ep, ic0_ep;
396 int f1, f2, f3;
397 int iR1, imR1, iR2, imR2, iRd, imRd, iRdmR1, iR1mRd;
398 int iR2pRd, imR2mRd, iR2mR1pRd, iR1mR2mRd;
399 int vR1[3], vR2[3], vRd[3], v[3];
400 int n;
401
402 nshell_op = EXX_Number_of_OP_Shells(xfm);
403 nshell_ep = EXX_Number_of_EP_Shells(xfm);
404
405 ncd_op = 2*nshell_op+1;
406 ncell_op = ncd_op*ncd_op*ncd_op;
407 ic0_op = (ncell_op-1)/2;
408
409 ncd_ep = 2*nshell_ep+1;
410 ncell_ep = ncd_ep*ncd_ep*ncd_ep;
411 ic0_ep = (ncell_ep-1)/2;
412
413 /* weight factor */
414 {
415 f1 = EXX_Index_Cmp_Sites(iatom1, iatom2, icell1, nshell_op); /* i=j ? */
416 f2 = EXX_Index_Cmp_Sites(iatom3, iatom4, icell2, nshell_op); /* k=l ? */
417
418 f3 = EXX_Index_Cmp_OP(iatom1, iatom2, iatom3, iatom4,
419 icell1, icell2, icell3, nshell_op, nshell_ep);
420
421 /* fool-proof test */
422 assert( f1 <= 0 );
423 assert( f2 <= 0 );
424 assert( f3 <= 0 );
425
426 if (0==f1 && 0==f2 && 0==f3) { *w = 8; }
427 else if (0==f1 && 0==f2 && -1==f3) { *w = 4; }
428 else if (-1==f1 && -1==f2 && -1==f3) { *w = 1; }
429 else { *w = 2; }
430 }
431
432 /* cell vectors */
433 {
434 /* IN : icell1, icell2, icell3, nshell_op, nshell_ep
435 OUT: iR1, imR1, iR2, imR2, iRd, imRd, iRdmR1, iR1mRd
436 iR2pRd, imR2mRd, iR2mR1pRd, iR1mR2mRd
437 TMP: vR1, vR2, vR3, v */
438
439 EXX_Index_Cell2XYZ(icell1, nshell_op, vR1);
440 EXX_Index_Cell2XYZ(icell2, nshell_op, vR2);
441 EXX_Index_Cell2XYZ(icell3, nshell_ep, vRd);
442
443 /* R1 */
444 v[0] = vR1[0];
445 v[1] = vR1[1];
446 v[2] = vR1[2];
447 iR1 = EXX_Index_XYZ2Cell(v, nshell_ep);
448
449 /* -R1 */
450 v[0] = -vR1[0];
451 v[1] = -vR1[1];
452 v[2] = -vR1[2];
453 imR1 = EXX_Index_XYZ2Cell(v, nshell_ep);
454
455 /* R2 */
456 v[0] = vR2[0];
457 v[1] = vR2[1];
458 v[2] = vR2[2];
459 iR2 = EXX_Index_XYZ2Cell(v, nshell_ep);
460
461 /* -R2 */
462 v[0] = -vR2[0];
463 v[1] = -vR2[1];
464 v[2] = -vR2[2];
465 imR2 = EXX_Index_XYZ2Cell(v, nshell_ep);
466
467 /* Rd */
468 v[0] = vRd[0];
469 v[1] = vRd[1];
470 v[2] = vRd[2];
471 iRd = EXX_Index_XYZ2Cell(v, nshell_ep);
472
473 /* -Rd */
474 v[0] = -vRd[0];
475 v[1] = -vRd[1];
476 v[2] = -vRd[2];
477 imRd = EXX_Index_XYZ2Cell(v, nshell_ep);
478
479 /* Rd-R1 */
480 v[0] = vRd[0]-vR1[0];
481 v[1] = vRd[1]-vR1[1];
482 v[2] = vRd[2]-vR1[2];
483 iRdmR1 = EXX_Index_XYZ2Cell(v, nshell_ep);
484
485 /* R1-Rd */
486 v[0] = vR1[0]-vRd[0];
487 v[1] = vR1[1]-vRd[1];
488 v[2] = vR1[2]-vRd[2];
489 iR1mRd = EXX_Index_XYZ2Cell(v, nshell_ep);
490
491 /* R2+Rd */
492 v[0] = vR2[0]+vRd[0];
493 v[1] = vR2[1]+vRd[1];
494 v[2] = vR2[2]+vRd[2];
495 iR2pRd = EXX_Index_XYZ2Cell(v, nshell_ep);
496
497 /* -R2-Rd */
498 v[0] = -vR2[0]-vRd[0];
499 v[1] = -vR2[1]-vRd[1];
500 v[2] = -vR2[2]-vRd[2];
501 imR2mRd = EXX_Index_XYZ2Cell(v, nshell_ep);
502
503 /* R2-R1+Rd */
504 v[0] = vR2[0]-vR1[0]+vRd[0];
505 v[1] = vR2[1]-vR1[1]+vRd[1];
506 v[2] = vR2[2]-vR1[2]+vRd[2];
507 iR2mR1pRd = EXX_Index_XYZ2Cell(v, nshell_ep);
508
509 /* R1-R2-Rd */
510 v[0] = vR1[0]-vR2[0]-vRd[0];
511 v[1] = vR1[1]-vR2[1]-vRd[1];
512 v[2] = vR1[2]-vR2[2]-vRd[2];
513 iR1mR2mRd = EXX_Index_XYZ2Cell(v, nshell_ep);
514 }
515
516 n = 0;
517
518 /* x1=(i,k,Rd), x2=(j,l,R2-R1+Rd), Rm=R1 */
519 if (iRd>=0 && iR2mR1pRd>=0) {
520 q_ep1[n] = EXX_Find_EP(xfm, iatom1, iatom3, iRd);
521 q_ep2[n] = EXX_Find_EP(xfm, iatom2, iatom4, iR2mR1pRd);
522 q_cell[n] = iR1;
523 assert(q_ep1[n] >= 0);
524 assert(q_ep2[n] >= 0);
525 assert(q_cell[n] >= 0);
526 n++;
527 }
528
529 /* x1=(j,k,Rd-R1), x2=(i,l,R2+Rd), Rm=-R1 */
530 if (iRdmR1>=0 && iR2pRd>=0) {
531 q_ep1[n] = EXX_Find_EP(xfm, iatom2, iatom3, iRdmR1);
532 q_ep2[n] = EXX_Find_EP(xfm, iatom1, iatom4, iR2pRd);
533 q_cell[n] = imR1;
534 assert(q_ep1[n] >= 0);
535 assert(q_ep2[n] >= 0);
536 assert(q_cell[n] >= 0);
537 n++;
538 }
539
540 /* x1=(i,l,R2+Rd), x2=(j,k,Rd-R1), Rm=R1 */
541 if (iR2pRd>=0 && iRdmR1>=0) {
542 q_ep1[n] = EXX_Find_EP(xfm, iatom1, iatom4, iR2pRd);
543 q_ep2[n] = EXX_Find_EP(xfm, iatom2, iatom3, iRdmR1);
544 q_cell[n] = iR1;
545 assert(q_ep1[n] >= 0);
546 assert(q_ep2[n] >= 0);
547 assert(q_cell[n] >= 0);
548 n++;
549 }
550
551 /* x1=(j,l,R2-R1+Rd), x2=(i,k,Rd), Rm=-R1 */
552 if (iR2mR1pRd>=0 && iRd>=0) {
553 q_ep1[n] = EXX_Find_EP(xfm, iatom2, iatom4, iR2mR1pRd);
554 q_ep2[n] = EXX_Find_EP(xfm, iatom1, iatom3, iRd);
555 q_cell[n] = imR1;
556 assert(q_ep1[n] >= 0);
557 assert(q_ep2[n] >= 0);
558 assert(q_cell[n] >= 0);
559 n++;
560 }
561
562 /* x1=(k,i,-Rd), x2=(l,j,R1-R2-Rd), Rm=R2 */
563 if (imR1>=0 && iR1mR2mRd>=0) {
564 q_ep1[n] = EXX_Find_EP(xfm, iatom3, iatom1, imR1);
565 q_ep2[n] = EXX_Find_EP(xfm, iatom4, iatom2, iR1mR2mRd);
566 q_cell[n] = iR2;
567 assert(q_ep1[n] >= 0);
568 assert(q_ep2[n] >= 0);
569 assert(q_cell[n] >= 0);
570 n++;
571 }
572
573 /* x1=(k,j,R1-Rd), x2=(l,i,-R2-Rd), Rm=R2 */
574 if (iR1mRd>=0 && imR2mRd>=0) {
575 q_ep1[n] = EXX_Find_EP(xfm, iatom3, iatom2, iR1mRd);
576 q_ep2[n] = EXX_Find_EP(xfm, iatom4, iatom1, imR2mRd);
577 q_cell[n] = iR2;
578 assert(q_ep1[n] >= 0);
579 assert(q_ep2[n] >= 0);
580 assert(q_cell[n] >= 0);
581 n++;
582 }
583
584 /* x1=(l,i,-R2-Rd), x2=(k,j,R1-Rd), Rm=-R2 */
585 if (imR2mRd>=0 && iR1mRd) {
586 q_ep1[n] = EXX_Find_EP(xfm, iatom4, iatom1, imR2mRd);
587 q_ep2[n] = EXX_Find_EP(xfm, iatom3, iatom2, iR1mRd);
588 q_cell[n] = imR2;
589 assert(q_ep1[n] >= 0);
590 assert(q_ep2[n] >= 0);
591 assert(q_cell[n] >= 0);
592 n++;
593 }
594
595 /* x1=(l,j,R1-R2-Rd), x2=(k,i,-Rd), Rm=-R2 */
596 if (iR1mR2mRd>=0 && imRd>=0) {
597 q_ep1[n] = EXX_Find_EP(xfm, iatom4, iatom2, iR1mR2mRd);
598 q_ep2[n] = EXX_Find_EP(xfm, iatom3, iatom1, imRd);
599 q_cell[n] = imR2;
600 assert(q_ep1[n] >= 0);
601 assert(q_ep2[n] >= 0);
602 assert(q_cell[n] >= 0);
603 n++;
604 }
605
606 *mul = n;
607 }
608
609
610
611 int EXX_Make_Quartets(
612 const EXX_t *xfm,
613 int *q_op1, /* [nq] */
614 int *q_op2, /* [nq] */
615 int *q_opd, /* [nq] */
616 int *q_wf, /* [nq] */
617 int *q_ep1, /* [nq*8] */
618 int *q_ep2, /* [nq*8] */
619 int *q_mul /* [nq] */
620 )
621 {
622 int i, f;
623 int nshell_op, nshell_ep, ncd_ep, ncell_ep, ic0_ep, nq, nop;
624 int iop1, iop2, iatom1, iatom2, iatom3, iatom4;
625 int icell1, icell2, icell3;
626 double rc1, rc2, rc3, rc4, rc12, rc34, d12, d34, d, cx12, cx34;
627 double c1[3], c2[3], c3[3], c4[3], c12[3], c34[3], c34_off[3];
628 double cc[3], cc_frac[3];
629 int cell_xyz[3];
630 int w, iep1[8], iep2[8], iRm[8], nmul;
631
632 const int *op_atom1, *op_atom2, *op_cell;
633 const double *atom_rc, *atom_v, *pvec;
634 double rc_cut, x, y, z;
635 int natom;
636 int ncd_op, ncell_op, ic0_op;
637
638 natom = EXX_natom(xfm);
639 nshell_op = EXX_Number_of_OP_Shells(xfm);
640 nshell_ep = EXX_Number_of_EP_Shells(xfm);
641 nop = EXX_Number_of_OP(xfm);
642
643 op_atom1 = EXX_Array_OP_Atom1(xfm);
644 op_atom2 = EXX_Array_OP_Atom2(xfm);
645 op_cell = EXX_Array_OP_Cell(xfm);
646
647 atom_rc = EXX_atom_rc(xfm);
648 atom_v = EXX_atom_v(xfm);
649 pvec = EXX_pvec(xfm);
650 rc_cut = EXX_rc_cut(xfm);
651
652 ncd_op = 2*nshell_op+1;
653 ncell_op = ncd_op*ncd_op*ncd_op;
654 ic0_op = (ncell_op-1)/2;
655
656 ncd_ep = 2*nshell_ep+1;
657 ncell_ep = ncd_ep*ncd_ep*ncd_ep;
658 ic0_ep = (ncell_ep-1)/2;
659
660 nq = 0;
661
662 for (iop1=0; iop1<nop; iop1++) {
663 iatom1 = op_atom1[iop1];
664 iatom2 = op_atom2[iop1];
665 icell1 = op_cell[iop1];
666
667 EXX_Vector_F2C(c1, &atom_v[3*iatom1], pvec);
668 EXX_Vector_F2C_Offsite(c2, &atom_v[3*iatom2], pvec, icell1, nshell_op);
669
670 rc1 = atom_rc[iatom1];
671 rc2 = atom_rc[iatom2];
672 d12 = EXX_Vector_Distance(c1, c2);
673 EXX_Vector_PAO_Overlap(rc1, rc2, d12, &rc12, &cx12);
674
675 c12[0] = cx12*c1[0] + (1.0-cx12)*c2[0];
676 c12[1] = cx12*c1[1] + (1.0-cx12)*c2[1];
677 c12[2] = cx12*c1[2] + (1.0-cx12)*c2[2];
678
679 for (iop2=0; iop2<nop; iop2++) {
680 iatom3 = op_atom1[iop2];
681 iatom4 = op_atom2[iop2];
682 icell2 = op_cell[iop2];
683
684 EXX_Vector_F2C(c3, &atom_v[3*iatom3], pvec);
685 EXX_Vector_F2C_Offsite(c4, &atom_v[3*iatom4], pvec, icell2, nshell_op);
686
687 rc3 = atom_rc[iatom3];
688 rc4 = atom_rc[iatom4];
689 d34 = EXX_Vector_Distance(c3, c4);
690 EXX_Vector_PAO_Overlap(rc3, rc4, d34, &rc34, &cx34);
691
692 c34[0] = cx34*c3[0] + (1.0-cx34)*c4[0];
693 c34[1] = cx34*c3[1] + (1.0-cx34)*c4[1];
694 c34[2] = cx34*c3[2] + (1.0-cx34)*c4[2];
695
696 for (icell3=ic0_ep; icell3<ncell_ep; icell3++) {
697 f = EXX_Index_Cmp_OP(iatom1, iatom2, iatom3, iatom4,
698 icell1, icell2, icell3,
699 nshell_op, nshell_ep);
700 if (f>0) { continue; }
701
702 EXX_Index_Cell2Cartesian(icell3, nshell_ep, pvec, cc);
703
704 c34_off[0] = c34[0] + cc[0];
705 c34_off[1] = c34[1] + cc[1];
706 c34_off[2] = c34[2] + cc[2];
707
708 d = EXX_Vector_Distance(c12, c34_off);
709 if (d > rc_cut + rc12 + rc34) { continue; }
710
711 quad(xfm, iatom1, iatom2, iatom3, iatom4,
712 icell1, icell2, icell3,
713 &w, iep1, iep2, iRm, &nmul);
714
715 if (q_op1) { q_op1[nq] = iop1; }
716 if (q_op2) { q_op2[nq] = iop2; }
717 if (q_opd) { q_opd[nq] = icell3; }
718 if (q_wf) { q_wf[nq] = w; }
719 for (i=0; i<nmul; i++) {
720 if (q_ep1) { q_ep1[8*nq+i] = iep1[i]; }
721 if (q_ep2) { q_ep2[8*nq+i] = iep2[i]; }
722 }
723 if (q_mul) { q_mul[nq] = nmul; }
724
725 nq++;
726 } /* icell3 */
727 } /* iop2 */
728 } /* iop1 */
729
730 return nq;
731 }
732 #endif
733
734