1 /*----------------------------------------------------------------------
2 exx_index.c
3
4 Indexing
5
6 Coded by M. Toyoda, 10/NOV/2009
7 ----------------------------------------------------------------------*/
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include "exx.h"
12 #include "exx_vector.h"
13 #include "exx_index.h"
14
15
16 /*----------------------------------------------------------------------
17 EXX_Index_Cmp_Sites
18
19 This defines relation between two off-site sites.
20 SITE1 is given by atom-site {iatom1} at the central cell, and SITE2
21 by atom-site {iatom2} at a cell {icell}.
22
23 SITE1 is SMALLER than SITE2 if
24 ic0<icell OR (ic0==icell AND iatom1<iatom2)),
25 SITE1 is LARGER than SITE2 if
26 ic0>icell OR (ic0==icell AND iatom1>iatom2)),
27 and SITE1 is EQUAL to SITE2 if
28 ic0==icell AND iatom1==iatom2.
29
30 This returns -1, +1, or 0 if SITE1 is smaller than, larger than, or
31 equal to SITE2, respectively.
32 ----------------------------------------------------------------------*/
EXX_Index_Cmp_Sites(int iatom1,int iatom2,int icell,int nshell)33 int EXX_Index_Cmp_Sites(
34 int iatom1,
35 int iatom2,
36 int icell,
37 int nshell
38 )
39 {
40 int ncd, ncell, ic0;
41
42 ncd = 2*nshell+1;
43 ncell = ncd*ncd*ncd;
44 ic0 = (ncell-1)/2;
45
46 if (icell>ic0) {
47 return -1;
48 } else if (icell==ic0) {
49 if (iatom1<iatom2) {
50 return -1;
51 } else if (iatom1==iatom2) {
52 return 0;
53 }
54 }
55
56 return 1;
57 }
58
59
60
61 /*----------------------------------------------------------------------
62 EXX_Index_Cmp_OP
63
64 This defines relation between two OPs.
65
66 SITE1 is given by atom-site {iatom1} at the central cell, and SITE2
67 by atom-site {iatom2} at a cell {icell}.
68
69 SITE1 is SMALLER than SITE2 if
70 ic0<icell OR (ic0==icell AND iatom1<iatom2)),
71 SITE1 is LARGER than SITE2 if
72 ic0>icell OR (ic0==icell AND iatom1>iatom2)),
73 and SITE1 is EQUAL to SITE2 if
74 ic0==icell AND iatom1==iatom2.
75
76 This returns -1, +1, or 0 if SITE1 is smaller than, larger than, or
77 equal to SITE2, respectively.
78 ----------------------------------------------------------------------*/
EXX_Index_Cmp_OP(int iatom1,int iatom2,int iatom3,int iatom4,int icell1,int icell2,int icell3,int nshell_op,int nshell_ep)79 int EXX_Index_Cmp_OP(
80 int iatom1,
81 int iatom2,
82 int iatom3,
83 int iatom4,
84 int icell1,
85 int icell2,
86 int icell3,
87 int nshell_op,
88 int nshell_ep
89 )
90 {
91 int ncd_op, ncell_op, ic0_op;
92 int ncd_ep, ncell_ep, ic0_ep;
93
94 ncd_op = 2*nshell_op+1;
95 ncell_op = ncd_op*ncd_op*ncd_op;
96 ic0_op = (ncell_op-1)/2;
97
98 ncd_ep = 2*nshell_ep+1;
99 ncell_ep = ncd_ep*ncd_ep*ncd_ep;
100 ic0_ep = (ncell_ep-1)/2;
101
102 if (ic0_ep>icell3) {
103 return +1;
104 } else if (ic0_ep<icell3) {
105 return -1;
106 } else {
107 /* icell3==ic0_ep */
108 if (iatom1>iatom3) {
109 return +1;
110 } else if (iatom1<iatom3) {
111 return -1;
112 } else {
113 /* iatom1==iatom3 */
114 if (icell1>icell2) {
115 return +1;
116 } else if (icell1<icell2) {
117 return -1;
118 } else {
119 /* icell1==icell2 */
120 if (iatom2>iatom4) {
121 return +1;
122 } else if (iatom2<iatom4) {
123 return -1;
124 } else {
125 /* iatom2==iatom4 */
126 return 0;
127 }
128 }
129 }
130 }
131 fprintf(stderr, "***ERROR: %s (%d)\n", __FILE__, __LINE__);
132 abort();
133 return 0;
134 }
135
136
137
138
139 /*----------------------------------------------------------------------
140 EXX_Index_Cell
141
142
143 ----------------------------------------------------------------------*/
EXX_Index_Cell2XYZ(int icell,int nshell,int c[3])144 void EXX_Index_Cell2XYZ(
145 int icell,
146 int nshell,
147 int c[3]
148 )
149 {
150 int ncd;
151 ncd = 2*nshell+1;
152 c[0] = icell%ncd - nshell;
153 c[1] = (icell/ncd)%ncd - nshell;
154 c[2] = (icell/ncd/ncd)%ncd - nshell;
155 }
156
157
EXX_Index_XYZ2Cell(int c[3],int nshell)158 int EXX_Index_XYZ2Cell(
159 int c[3],
160 int nshell
161 )
162 {
163 int ncd;
164 ncd = 2*nshell+1;
165
166 if (abs(c[0])>nshell || abs(c[1])>nshell || abs(c[2])>nshell) {
167 return -1;
168 }
169
170 return (c[0]+nshell) + (c[1]+nshell)*ncd + (c[2]+nshell)*ncd*ncd;
171 }
172
173
EXX_Index_Cell2Cartesian(int icell,int nshell,const double pvec[9],double c[3])174 void EXX_Index_Cell2Cartesian(
175 int icell,
176 int nshell,
177 const double pvec[9],
178 double c[3]
179 )
180 {
181 int ncd;
182 double cf[3];
183 ncd = 2*nshell+1;
184 cf[0] = (double)(icell%ncd - nshell);
185 cf[1] = (double)((icell/ncd)%ncd - nshell);
186 cf[2] = (double)((icell/ncd/ncd)%ncd - nshell);
187 EXX_Vector_F2C(c, cf, pvec);
188 }
189
190
191 /*----------------------------------------------------------------------
192 EXX_Index_OP_NShell
193
194 Minimun number of shells for overlaping pairs (OPs).
195
196 IN:
197 natom : number of atoms in the unit cell
198 atom_rc : confinement length of atoms
199 atom_v : position of atoms in fractional cooridnate
200 pvec : primitive translational vectors
201
202 OUT:
203 return : minimum number of shells for OPs.
204 If something wrong happens, return value is -1.
205
206 Note:
207 {atom_v} must be given in fractional coordinate. For example, if
208 the atom is at the body-center of the unit cell, then the {atom_v}
209 for the atom should be {0.5, 0.5, 0.5}.
210 ----------------------------------------------------------------------*/
EXX_Index_OP_NShell(int natom,const double * atom_rc,const double * atom_v,const double * pvec)211 int EXX_Index_OP_NShell(
212 int natom,
213 const double *atom_rc,
214 const double *atom_v,
215 const double *pvec
216 )
217 {
218 int i, nop, nop0;
219
220 /* number of OPs wthin the central cell */
221 nop0 = EXX_Index_OP(natom, 0, atom_rc, atom_v, pvec, NULL, NULL, NULL);
222
223 /* increase nshell untill the number of OPs converges */
224 for (i=1; i<100; i++) {
225 nop = EXX_Index_OP(natom, i, atom_rc, atom_v, pvec, NULL, NULL, NULL);
226 if (nop0 == nop) { return i-1; } /* converged */
227 nop0 = nop;
228 }
229
230 return -1; /* no convergence */
231 }
232
233
234
235
236 /*----------------------------------------------------------------------
237 EXX_Index_EP_NShell
238
239 Minimun number of shells for exchanging pairs (EPs).
240
241 IN:
242 natom : number of atoms in the unit cell
243 atom_rc : confinement length of atoms
244 atom_v : position of atoms in fractional cooridnate
245 pvec : primitive translational vectors
246 rc_cut : truncation length
247
248 OUT:
249 return : minimum number of shells for OPs.
250 If something wrong happens, return value is -1.
251
252 Note:
253 {atom_v} must be given in fractional coordinate. For example, if
254 the atom is at the body-center of the unit cell, then the {atom_v}
255 for the atom should be {0.5, 0.5, 0.5}.
256 ----------------------------------------------------------------------*/
EXX_Index_EP_NShell(int natom,const double * atom_rc,const double * atom_v,const double * pvec,double rc_cut)257 int EXX_Index_EP_NShell(
258 int natom,
259 const double *atom_rc,
260 const double *atom_v,
261 const double *pvec,
262 double rc_cut
263 )
264 {
265 int i, nep, nep0;
266
267 /* number of OPs wthin the central cell */
268 nep0 = EXX_Index_EP(natom, 0, atom_rc, atom_v, pvec, rc_cut,
269 NULL, NULL, NULL);
270
271 /* increase nshell untill the number of OPs converges */
272 for (i=1; i<100; i++) {
273 nep = EXX_Index_EP(natom, i, atom_rc, atom_v, pvec, rc_cut,
274 NULL, NULL, NULL);
275 if (nep0 == nep) { return i-1; } /* converged */
276 nep0 = nep;
277 }
278
279 return -1; /* no convergence */
280 }
281
282
283
284
285 /*----------------------------------------------------------------------
286 EXX_Index_OP
287
288 List up OPs.
289 Array for {op_atom1}, {op_atom2}, and {op_cell} must be of enough
290 length. One can know the required minimum length for the arrays by
291 calling this function with specifing NULL for {op_atom1}, {op_atom2},
292 and {op_cell}.
293
294 A typical usage of this function is as follows:
295 ---
296 nshell = EXX_Index_OP_NShell(natom, atom_rc, atom_v, pvec);
297 nop = CCM_OP_Enum(natom, nshell, atom_rc, atom_v, pvec,
298 NULL, NULL, NULL);
299 op_atom1 = (int*)malloc(sizeof(int)*nop);
300 op_atom2 = (int*)malloc(sizeof(int)*nop);
301 op_cell = (int*)malloc(sizeof(int)*nop);
302 CCM_OP_Enum(natom, nshell, atom_rc, atom_v, pvec,
303 op_atom1, op_atom2, op_cell);
304 ---
305
306 IN:
307 natom : number of atoms in the unit cell.
308 nshell : number of shells to be considered.
309 atom_rc : confinement length of atoms.
310 atom_v : atomic positions in fractional coordinates.
311 pvec : primitive translational vectors.
312
313 OUT:
314 op_atom1 : array for atom1 of OPs (optional)
315 op_atom2 : array for atom2 of OPs (optional)
316 op_cell : array for cells of OPs (optional)
317 return : number of OPs
318 ----------------------------------------------------------------------*/
EXX_Index_OP(int natom,int nshell,const double * atom_rc,const double * atom_v,const double * pvec,int * op_atm1,int * op_atm2,int * op_cell)319 int EXX_Index_OP(
320 int natom,
321 int nshell,
322 const double *atom_rc,
323 const double *atom_v,
324 const double *pvec,
325 int *op_atm1, /* OUT, optional */
326 int *op_atm2, /* OUT, optional */
327 int *op_cell /* OUT, optional */
328 )
329 {
330 int iatm1, iatm2, icell, ncell, ncd, nop, ic0;
331 double c1[3], c2[3], rc1, rc2, cc[3], x, y, z, d;
332
333 /* number of cells to consider */
334 ncd = 2*nshell+1;
335 ncell = ncd*ncd*ncd;
336 ic0 = (ncell-1)/2;
337 nop = 0;
338
339 for (iatm1=0; iatm1<natom; iatm1++) {
340 /* atom 1 in the central cell */
341 EXX_Vector_F2C(c1, &atom_v[3*iatm1], pvec);
342 /* frac2cart(c1, &atom_v[3*iatm1], pvec, ic0, nshell); */
343 rc1 = atom_rc[iatm1];
344
345 for (icell=ic0; icell<ncell; icell++) {
346 /* lattice cell */
347
348 for (iatm2=0; iatm2<natom; iatm2++) {
349 if (icell==ic0 && iatm2 < iatm1 ) { continue; }
350 /* atom 2 in the lattice cell */
351 EXX_Vector_F2C_Offsite(c2, &atom_v[3*iatm2], pvec, icell, nshell);
352 /* frac2cart(c2, &atom_v[3*iatm2], pvec, icell, nshell); */
353 rc2 = atom_rc[iatm2];
354
355 d = EXX_Vector_Distance(c1, c2);
356 if (d > rc1+rc2) { continue; }
357
358 if (op_atm1) { op_atm1[nop] = iatm1; }
359 if (op_atm2) { op_atm2[nop] = iatm2; }
360 if (op_cell) { op_cell[nop] = icell; }
361 nop++;
362
363 } /* loop of iatm2 */
364 } /* loop of icell */
365 } /* loop of iatm1 */
366
367 return nop;
368 }
369
370
371
372
373
374 /*----------------------------------------------------------------------
375 EXX_Index_EQ
376
377 List up EQs.
378
379 IN:
380 natom : number of atoms in the unit cell.
381 nshell : number of shells to be considered.
382 atom_rc : confinement length of atoms.
383 atom_v : atomic positions in fractional coordinates.
384 pvec : primitive translational vectors.
385 rc_cut : truncation length
386
387 OUT:
388 eq_atom1 : array for atom1 of EQs (optional)
389 eq_atom2 : array for atom2 of EQs (optional)
390 eq_atom3 : array for atom3 of EQs (optional)
391 eq_atom4 : array for atom4 of EQs (optional)
392 eq_cell1 : array for cell1 of EQs (optional)
393 eq_cell2 : array for cell2 of EQs (optional)
394 eq_cell3 : array for cell3 of EQs (optional)
395 return : number of EQs
396 ----------------------------------------------------------------------*/
EXX_Index_EQ(int natom,int nshell,const double * atom_rc,const double * atom_v,const double * pvec,double rc_cut,int * eq_atom1,int * eq_atom2,int * eq_atom3,int * eq_atom4,int * eq_cell1,int * eq_cell2,int * eq_cell3)397 int EXX_Index_EQ(
398 int natom,
399 int nshell,
400 const double *atom_rc,
401 const double *atom_v,
402 const double *pvec,
403 double rc_cut,
404 int *eq_atom1, /* OUT, optional */
405 int *eq_atom2, /* OUT, optional */
406 int *eq_atom3, /* OUT, optional */
407 int *eq_atom4, /* OUT, optional */
408 int *eq_cell1, /* OUT, optional */
409 int *eq_cell2, /* OUT, optional */
410 int *eq_cell3 /* OUT, optional */
411 )
412 {
413 int iatm1, iatm2, iatm3, iatm4, icell1, icell2, icell3;
414 int neq, ncell, ncd, nshell_op, ncell_op, ncd_op, ic0_op;
415 int iscr;
416 double x, y, z, rc1, rc2, rc3, rc4, rc12, rc34, cx12, cx34;
417 double c1[3], c2[3], c3[3], c4[3], c12[3], c34[3], c34_off[3];
418 double cc1[3], cc2[3], cc3[3];
419 double d12, d34, d;
420
421 /* cells for OP */
422 nshell_op = EXX_Index_OP_NShell(natom, atom_rc, atom_v, pvec);
423 ncd_op = 2*nshell_op+1;
424 ncell_op = ncd_op*ncd_op*ncd_op;
425 ic0_op = (ncell_op-1)/2;
426
427 /* cells for EQ */
428 ncd = 2*nshell+1;
429 ncell = ncd*ncd*ncd;
430
431 neq = 0;
432
433 /* bare/screneed */
434 iscr = (rc_cut>1e-10);
435
436 for (iatm1=0; iatm1<natom; iatm1++) {
437 /* atom 1 in the central cell */
438 EXX_Vector_F2C(c1, &atom_v[3*iatm1], pvec);
439 /* frac2cart(c1, &atom_v[3*iatm1], pvec, ic0_op, nshell_op); */
440 rc1 = atom_rc[iatm1];
441
442 for (icell1=0; icell1<ncell_op; icell1++) {
443 /* lattice cell */
444
445 for (iatm2=0; iatm2<natom; iatm2++) {
446 /* atom 2 in the lattice cell */
447 EXX_Vector_F2C_Offsite(c2, &atom_v[3*iatm2], pvec, icell1, nshell_op);
448 /* frac2cart(c2, &atom_v[3*iatm2], pvec, icell1, nshell_op); */
449 rc2 = atom_rc[iatm2];
450
451 /* d12 = distance(c1, c2); */
452 d12 = EXX_Vector_Distance(c1, c2);
453 if (d12 > rc1+rc2) { continue; }
454
455 EXX_Vector_PAO_Overlap(rc1, rc2, d12, &rc12, &cx12);
456 /* slo_pair(rc1, rc2, d12, &rc12, &cx12); */
457 c12[0] = cx12*c1[0] + (1.0-cx12)*c2[0];
458 c12[1] = cx12*c1[1] + (1.0-cx12)*c2[1];
459 c12[2] = cx12*c1[2] + (1.0-cx12)*c2[2];
460
461 for (iatm3=0; iatm3<natom; iatm3++) {
462 /* atom 3 in the central cell */
463 EXX_Vector_F2C(c3, &atom_v[3*iatm3], pvec);
464 /* frac2cart(c3, &atom_v[3*iatm3], pvec, ic0_op, nshell_op); */
465 rc3 = atom_rc[iatm3];
466
467 for (icell2=0; icell2<ncell_op; icell2++) {
468 /* lattice cell */
469
470 for (iatm4=0; iatm4<natom; iatm4++) {
471 /* atom 4 in the lattice cell */
472 EXX_Vector_F2C_Offsite(c4, &atom_v[3*iatm4], pvec,
473 icell2, nshell_op);
474 /* frac2cart(c4, &atom_v[3*iatm4], pvec, icell2, nshell_op); */
475 rc4 = atom_rc[iatm4];
476
477 /* d34 = distance(c3, c4); */
478 d34 = EXX_Vector_Distance(c3, c4);
479 if (d34 > rc3+rc4) { continue; }
480
481 EXX_Vector_PAO_Overlap(rc3, rc4, d34, &rc34, &cx34);
482 /* slo_pair(rc3, rc4, d34, &rc34, &cx34); */
483 c34[0] = cx34*c3[0] + (1.0-cx34)*c4[0];
484 c34[1] = cx34*c3[1] + (1.0-cx34)*c4[1];
485 c34[2] = cx34*c3[2] + (1.0-cx34)*c4[2];
486
487 for (icell3=0; icell3<ncell; icell3++) {
488 /* lattice cell */
489 x = (double)( icell3%ncd - nshell );
490 y = (double)( (icell3/ncd)%ncd - nshell );
491 z = (double)( (icell3/ncd/ncd)%ncd - nshell );
492 cc3[0] = pvec[0]*x + pvec[3]*y + pvec[6]*z;
493 cc3[1] = pvec[1]*x + pvec[4]*y + pvec[7]*z;
494 cc3[2] = pvec[2]*x + pvec[5]*y + pvec[8]*z;
495
496 c34_off[0] = c34[0] + cc3[0];
497 c34_off[1] = c34[1] + cc3[1];
498 c34_off[2] = c34[2] + cc3[2];
499
500 if (iscr) {
501 /* d = distance(c12, c34_off); */
502 d = EXX_Vector_Distance(c12, c34_off);
503 if (d > rc_cut + rc12 + rc34) { continue; }
504 }
505
506 if (eq_atom1) { eq_atom1[neq] = iatm1; }
507 if (eq_atom2) { eq_atom2[neq] = iatm2; }
508 if (eq_atom3) { eq_atom3[neq] = iatm3; }
509 if (eq_atom4) { eq_atom4[neq] = iatm4; }
510 if (eq_cell1) { eq_cell1[neq] = icell1; }
511 if (eq_cell2) { eq_cell2[neq] = icell2; }
512 if (eq_cell3) { eq_cell3[neq] = icell3; }
513
514 neq++;
515 }
516 } /* iatm4 */
517 } /* icell2 */
518 } /* iatm3 */
519 } /* iatm2 */
520 } /* icell1 */
521 } /* iatm1 */
522
523 return neq;
524 }
525
526
527
528
529 /*----------------------------------------------------------------------
530 EXX_Index_EP
531
532 List up EQs.
533
534 IN:
535 natom : number of atoms in the unit cell.
536 nshell : number of shells to be considered.
537 atom_rc : confinement length of atoms.
538 atom_v : atomic positions in fractional coordinates.
539 pvec : primitive translational vectors.
540 rc_cut : truncation length
541
542 OUT:
543 ep_atom1 : array for atom1 of EPs (optional)
544 ep_atom2 : array for atom2 of EPs (optional)
545 ep_cell : array for cells of EPs (optional)
546 return : number of EPs
547 ----------------------------------------------------------------------*/
EXX_Index_EP(int natom,int nshell,const double * atom_rc,const double * atom_v,const double * pvec,double rc_cut,int * ep_atom1,int * ep_atom2,int * ep_cell)548 int EXX_Index_EP(
549 int natom,
550 int nshell,
551 const double *atom_rc,
552 const double *atom_v,
553 const double *pvec,
554 double rc_cut,
555 int *ep_atom1, /* OUT, optional */
556 int *ep_atom2, /* OUT, optional */
557 int *ep_cell /* OUT, optional */
558 )
559 {
560 int iatm1, iatm2, iatm3, iatm4, icell1, icell2, icell3;
561 int nshell_op, ncd_op, ncell_op, ic0_op, ncd, ncell, iflag, nep;
562 int iscr;
563 double x, y, z, c1[3], c2[3], c3[3], c4[3], cc3[3], c12[3], c34[3];
564 double rc1, rc2, rc3, rc4, d12, d34, rc12, rc34, cx12, cx34;
565 double c34_off[3], d;
566
567 /* cells for OP */
568 nshell_op = EXX_Index_OP_NShell(natom, atom_rc, atom_v, pvec);
569 ncd_op = 2*nshell_op+1;
570 ncell_op = ncd_op*ncd_op*ncd_op;
571 ic0_op = (ncell_op-1)/2;
572
573 /* cells for EQ */
574 ncd = 2*nshell+1;
575 ncell = ncd*ncd*ncd;
576
577 nep = 0;
578
579 /* bare/screen */
580 iscr = (rc_cut>1e-10);
581
582 for (iatm1=0; iatm1<natom; iatm1++) {
583 /* atom 1 in the central cell */
584 EXX_Vector_F2C(c1, &atom_v[3*iatm1], pvec);
585 /* frac2cart(c1, &atom_v[3*iatm1], pvec, ic0_op, nshell_op); */
586 rc1 = atom_rc[iatm1];
587
588 for (iatm3=0; iatm3<natom; iatm3++) {
589 /* atom 3 in the central cell */
590 EXX_Vector_F2C(c3, &atom_v[3*iatm3], pvec);
591 /* frac2cart(c3, &atom_v[3*iatm3], pvec, ic0_op, nshell_op); */
592 rc3 = atom_rc[iatm3];
593
594 for (icell3=0; icell3<ncell; icell3++) {
595 /* lattice cell */
596 x = (double)( icell3%ncd - nshell );
597 y = (double)( (icell3/ncd)%ncd - nshell );
598 z = (double)( (icell3/ncd/ncd)%ncd - nshell );
599 cc3[0] = pvec[0]*x + pvec[3]*y + pvec[6]*z;
600 cc3[1] = pvec[1]*x + pvec[4]*y + pvec[7]*z;
601 cc3[2] = pvec[2]*x + pvec[5]*y + pvec[8]*z;
602
603 iflag = 0;
604
605 for (icell1=0; icell1<ncell_op && 0==iflag; icell1++) {
606 /* lattice cell */
607
608 for (iatm2=0; iatm2<natom && 0==iflag; iatm2++) {
609 /* atom 2 in the lattice cell */
610 EXX_Vector_F2C_Offsite(c2, &atom_v[3*iatm2], pvec,
611 icell1, nshell_op);
612 /* frac2cart(c2, &atom_v[3*iatm2], pvec, icell1, nshell_op); */
613 rc2 = atom_rc[iatm2];
614
615 /* d12 = distance(c1, c2); */
616 d12 = EXX_Vector_Distance(c1, c2);
617 if (d12 > rc1+rc2) { continue; }
618
619 EXX_Vector_PAO_Overlap(rc1, rc2, d12, &rc12, &cx12);
620 /* slo_pair(rc1, rc2, d12, &rc12, &cx12); */
621 c12[0] = cx12*c1[0] + (1.0-cx12)*c2[0];
622 c12[1] = cx12*c1[1] + (1.0-cx12)*c2[1];
623 c12[2] = cx12*c1[2] + (1.0-cx12)*c2[2];
624
625 for (icell2=0; icell2<ncell_op && 0==iflag; icell2++) {
626 /* lattice cell */
627
628 for (iatm4=0; iatm4<natom && 0==iflag; iatm4++) {
629 /* atom 4 in the lattice cell */
630 EXX_Vector_F2C_Offsite(c4, &atom_v[3*iatm4], pvec,
631 icell2, nshell_op);
632 /* frac2cart(c4, &atom_v[3*iatm4], pvec, icell2, nshell_op); */
633 rc4 = atom_rc[iatm4];
634
635 /* d34 = distance(c3, c4); */
636 d34 = EXX_Vector_Distance(c3, c4);
637 if (d34 > rc3+rc4) { continue; }
638
639 if (iscr) {
640 EXX_Vector_PAO_Overlap(rc3, rc4, d34, &rc34, &cx34);
641 /* slo_pair(rc3, rc4, d34, &rc34, &cx34); */
642 c34[0] = cx34*c3[0] + (1.0-cx34)*c4[0];
643 c34[1] = cx34*c3[1] + (1.0-cx34)*c4[1];
644 c34[2] = cx34*c3[2] + (1.0-cx34)*c4[2];
645
646 c34_off[0] = c34[0] + cc3[0];
647 c34_off[1] = c34[1] + cc3[1];
648 c34_off[2] = c34[2] + cc3[2];
649
650 /* d = distance(c12, c34_off); */
651 d = EXX_Vector_Distance(c12, c34_off);
652 if (d > rc_cut + rc12 + rc34) { continue; }
653 }
654 iflag = 1;
655 } /* iatm4 */
656 } /* icell2 */
657 } /* iatm2 */
658 } /* icell1 */
659
660 if (iflag) {
661 if (ep_atom1) { ep_atom1[nep] = iatm1; }
662 if (ep_atom2) { ep_atom2[nep] = iatm3; }
663 if (ep_cell) { ep_cell[nep] = icell3; }
664 nep++;
665 }
666
667 } /* icell3*/
668 } /* iatm2 */
669 } /* iatm1 */
670
671 return nep;
672 }
673
674
675
EXX_Index_NQ_Full(int natom,int nshell,const double * atom_rc,const double * atom_v,const double * pvec,double rc_cut)676 int EXX_Index_NQ_Full(
677 int natom,
678 int nshell,
679 const double *atom_rc,
680 const double *atom_v,
681 const double *pvec,
682 double rc_cut
683 )
684 {
685 int iatm1, iatm2, iatm3, iatm4, icell1, icell2, icell3, ic0;
686 int nq, ncell, ncd, nshell_op, ncell_op, ncd_op, ic0_op;
687 int iscr;
688 double x, y, z, rc1, rc2, rc3, rc4, rc12, rc34, cx12, cx34;
689 double c1[3], c2[3], c3[3], c4[3], c12[3], c34[3], c34_off[3];
690 double cc1[3], cc2[3], cc3[3];
691 double d12, d34, d;
692
693 /* cells for OP */
694 nshell_op = EXX_Index_OP_NShell(natom, atom_rc, atom_v, pvec);
695 ncd_op = 2*nshell_op+1;
696 ncell_op = ncd_op*ncd_op*ncd_op;
697 ic0_op = (ncell_op-1)/2;
698
699 /* cells for EQ */
700 ncd = 2*nshell+1;
701 ncell = ncd*ncd*ncd;
702 ic0 = (ncell-1)/2;
703
704 nq = 0;
705
706 /* bear/screened */
707 iscr = (rc_cut>1e-10);
708
709 for (iatm1=0; iatm1<natom; iatm1++) {
710 /* atom 1 in the central cell */
711 EXX_Vector_F2C(c1, &atom_v[3*iatm1], pvec);
712 rc1 = atom_rc[iatm1];
713
714 for (icell1=0; icell1<ncell_op; icell1++) {
715 /* lattice cell */
716
717 for (iatm2=0; iatm2<natom; iatm2++) {
718 /* atom 2 in the lattice cell */
719
720 EXX_Vector_F2C_Offsite(c2, &atom_v[3*iatm2], pvec, icell1, nshell_op);
721 rc2 = atom_rc[iatm2];
722
723 /* d12 = distance(c1, c2); */
724 d12 = EXX_Vector_Distance(c1, c2);
725 if (d12 > rc1+rc2) { continue; }
726
727 EXX_Vector_PAO_Overlap(rc1, rc2, d12, &rc12, &cx12);
728 c12[0] = cx12*c1[0] + (1.0-cx12)*c2[0];
729 c12[1] = cx12*c1[1] + (1.0-cx12)*c2[1];
730 c12[2] = cx12*c1[2] + (1.0-cx12)*c2[2];
731
732 for (iatm3=0; iatm3<natom; iatm3++) {
733 /* atom 3 in the central cell */
734 EXX_Vector_F2C(c3, &atom_v[3*iatm3], pvec);
735 rc3 = atom_rc[iatm3];
736
737 for (icell2=0; icell2<ncell_op; icell2++) {
738 /* lattice cell */
739
740 for (iatm4=0; iatm4<natom; iatm4++) {
741
742 /* atom 4 in the lattice cell */
743 EXX_Vector_F2C_Offsite(c4, &atom_v[3*iatm4], pvec,
744 icell2, nshell_op);
745 rc4 = atom_rc[iatm4];
746
747 /* d34 = distance(c3, c4); */
748 d34 = EXX_Vector_Distance(c3, c4);
749 if (d34 > rc3+rc4) { continue; }
750
751 EXX_Vector_PAO_Overlap(rc3, rc4, d34, &rc34, &cx34);
752 c34[0] = cx34*c3[0] + (1.0-cx34)*c4[0];
753 c34[1] = cx34*c3[1] + (1.0-cx34)*c4[1];
754 c34[2] = cx34*c3[2] + (1.0-cx34)*c4[2];
755
756 for (icell3=0; icell3<ncell; icell3++) {
757 /* lattice cell */
758 if (iscr) {
759 x = (double)( icell3%ncd - nshell );
760 y = (double)( (icell3/ncd)%ncd - nshell );
761 z = (double)( (icell3/ncd/ncd)%ncd - nshell );
762 cc3[0] = pvec[0]*x + pvec[3]*y + pvec[6]*z;
763 cc3[1] = pvec[1]*x + pvec[4]*y + pvec[7]*z;
764 cc3[2] = pvec[2]*x + pvec[5]*y + pvec[8]*z;
765
766 c34_off[0] = c34[0] + cc3[0];
767 c34_off[1] = c34[1] + cc3[1];
768 c34_off[2] = c34[2] + cc3[2];
769
770 d = EXX_Vector_Distance(c12, c34_off);
771 if (d > rc_cut + rc12 + rc34) { continue; }
772 }
773
774 nq++;
775 }
776 } /* iatm4 */
777 } /* icell2 */
778 } /* iatm3 */
779 } /* iatm2 */
780 } /* icell1 */
781 } /* iatm1 */
782
783 return nq;
784 }
785
786
787
EXX_Index_NQ_Reduced(int natom,int nshell,const double * atom_rc,const double * atom_v,const double * pvec,double rc_cut)788 int EXX_Index_NQ_Reduced(
789 int natom,
790 int nshell,
791 const double *atom_rc,
792 const double *atom_v,
793 const double *pvec,
794 double rc_cut
795 )
796 {
797 int iatm1, iatm2, iatm3, iatm4, icell1, icell2, icell3, ic0;
798 int nq, ncell, ncd, nshell_op, ncell_op, ncd_op, ic0_op;
799 int iscr;
800 double x, y, z, rc1, rc2, rc3, rc4, rc12, rc34, cx12, cx34;
801 double c1[3], c2[3], c3[3], c4[3], c12[3], c34[3], c34_off[3];
802 double cc1[3], cc2[3], cc3[3];
803 double d12, d34, d;
804 int f;
805
806 /* cells for OP */
807 nshell_op = EXX_Index_OP_NShell(natom, atom_rc, atom_v, pvec);
808 ncd_op = 2*nshell_op+1;
809 ncell_op = ncd_op*ncd_op*ncd_op;
810 ic0_op = (ncell_op-1)/2;
811
812 /* cells for EQ */
813 ncd = 2*nshell+1;
814 ncell = ncd*ncd*ncd;
815 ic0 = (ncell-1)/2;
816
817 nq = 0;
818
819 /* bare/screened */
820 iscr = (rc_cut>1e-10);
821
822 for (iatm1=0; iatm1<natom; iatm1++) {
823 /* atom 1 in the central cell */
824 EXX_Vector_F2C(c1, &atom_v[3*iatm1], pvec);
825 rc1 = atom_rc[iatm1];
826
827 for (icell1=0; icell1<ncell_op; icell1++) {
828 /* lattice cell */
829
830 for (iatm2=0; iatm2<natom; iatm2++) {
831 /* atom 2 in the lattice cell */
832 f = EXX_Index_Cmp_Sites(iatm1, iatm2, icell1, nshell_op);
833 if (f>0) { continue; }
834
835 EXX_Vector_F2C_Offsite(c2, &atom_v[3*iatm2], pvec, icell1, nshell_op);
836 rc2 = atom_rc[iatm2];
837
838 /* d12 = distance(c1, c2); */
839 d12 = EXX_Vector_Distance(c1, c2);
840 if (d12 > rc1+rc2) { continue; }
841
842 EXX_Vector_PAO_Overlap(rc1, rc2, d12, &rc12, &cx12);
843 c12[0] = cx12*c1[0] + (1.0-cx12)*c2[0];
844 c12[1] = cx12*c1[1] + (1.0-cx12)*c2[1];
845 c12[2] = cx12*c1[2] + (1.0-cx12)*c2[2];
846
847 for (iatm3=0; iatm3<natom; iatm3++) {
848 /* atom 3 in the central cell */
849 EXX_Vector_F2C(c3, &atom_v[3*iatm3], pvec);
850 rc3 = atom_rc[iatm3];
851
852 for (icell2=0; icell2<ncell_op; icell2++) {
853 /* lattice cell */
854
855 for (iatm4=0; iatm4<natom; iatm4++) {
856 f = EXX_Index_Cmp_Sites(iatm3, iatm4, icell2, nshell_op);
857 if (f>0) { continue; }
858
859 /* atom 4 in the lattice cell */
860 EXX_Vector_F2C_Offsite(c4, &atom_v[3*iatm4], pvec,
861 icell2, nshell_op);
862 rc4 = atom_rc[iatm4];
863
864 /* d34 = distance(c3, c4); */
865 d34 = EXX_Vector_Distance(c3, c4);
866 if (d34 > rc3+rc4) { continue; }
867
868 EXX_Vector_PAO_Overlap(rc3, rc4, d34, &rc34, &cx34);
869 c34[0] = cx34*c3[0] + (1.0-cx34)*c4[0];
870 c34[1] = cx34*c3[1] + (1.0-cx34)*c4[1];
871 c34[2] = cx34*c3[2] + (1.0-cx34)*c4[2];
872
873 for (icell3=0; icell3<ncell; icell3++) {
874 f = EXX_Index_Cmp_OP(iatm1, iatm2, iatm3, iatm4,
875 icell1, icell2, icell3, nshell_op, nshell);
876 if (f>0) { continue; }
877
878 if (iscr) {
879 /* lattice cell */
880 x = (double)( icell3%ncd - nshell );
881 y = (double)( (icell3/ncd)%ncd - nshell );
882 z = (double)( (icell3/ncd/ncd)%ncd - nshell );
883 cc3[0] = pvec[0]*x + pvec[3]*y + pvec[6]*z;
884 cc3[1] = pvec[1]*x + pvec[4]*y + pvec[7]*z;
885 cc3[2] = pvec[2]*x + pvec[5]*y + pvec[8]*z;
886
887 c34_off[0] = c34[0] + cc3[0];
888 c34_off[1] = c34[1] + cc3[1];
889 c34_off[2] = c34[2] + cc3[2];
890
891 d = EXX_Vector_Distance(c12, c34_off);
892 if (d > rc_cut + rc12 + rc34) { continue; }
893 }
894
895 nq++;
896 }
897 } /* iatm4 */
898 } /* icell2 */
899 } /* iatm3 */
900 } /* iatm2 */
901 } /* icell1 */
902 } /* iatm1 */
903
904 return nq;
905 }
906
907
908
909