1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include "openmx_common.h"
5 #include "exx.h"
6 #include "exx_debug.h"
7 #include "exx_log.h"
8
9
EXX_Debug_Copy_DM(int MaxN,double ***** CDM,EXX_t * exx,dcomplex **** exx_CDM,int symbrk)10 void EXX_Debug_Copy_DM(
11 int MaxN,
12 double *****CDM,
13 EXX_t *exx,
14 dcomplex ****exx_CDM,
15 int symbrk
16 )
17 {
18 int spin, i, j, myrank, nproc, iproc, ib;
19 int GA_AN, GB_AN, MA_AN, LB_AN;
20 int iep, nep, iatm, nb, nb1, nb2, nbmax;
21 const int *ep_atom1, *ep_atom2;
22 double dm1, dm2, diff;
23 double w, du, dd;
24
25 int nbuf;
26 double *buffer;
27 MPI_Comm comm;
28
29 nep = EXX_Number_of_EP(exx);
30 ep_atom1 = EXX_Array_EP_Atom1(exx);
31 ep_atom2 = EXX_Array_EP_Atom2(exx);
32
33 /* MPI information */
34 comm = g_exx_mpicomm;
35 MPI_Comm_rank(comm, &myrank);
36 MPI_Comm_size(comm, &nproc);
37
38 /* max number of basis */
39 nbmax = 0;
40 for (i=1; i<=atomnum; i++) {
41 nb = Spe_Total_CNO[WhatSpecies[i]];
42 if (nb>nbmax) { nbmax=nb; }
43 }
44
45 /* allocate buffer */
46 nbuf = nep*nbmax*nbmax;
47 buffer = (double*)malloc(sizeof(double)*nbuf);
48
49 for (spin=0; spin<=SpinP_switch; spin++){
50
51 /* clear buffer */
52 for (i=0; i<nbuf; i++) { buffer[i] = 0.0; }
53
54 /* copy DM to buffer */
55 for (iep=0; iep<nep; iep++) {
56 GA_AN = ep_atom1[iep]+1;
57 GB_AN = ep_atom2[iep]+1;
58
59 nb1 = Spe_Total_CNO[WhatSpecies[GA_AN]];
60 nb2 = Spe_Total_CNO[WhatSpecies[GB_AN]];
61
62 /* find MA_AN */
63 for (MA_AN=1; MA_AN<=Matomnum; MA_AN++) {
64 if (GA_AN==M2G[MA_AN]) { break; }
65 }
66 if (GA_AN != M2G[MA_AN]) { continue; }
67
68 /* find LB_AN */
69 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
70 if (GB_AN==natn[GA_AN][LB_AN]) { break; }
71 }
72 if (GB_AN!=natn[GA_AN][LB_AN]) { continue; }
73
74 for (i=0; i<nb1; i++){
75 for (j=0; j<nb2; j++){
76 ib = iep*nbmax*nbmax + i*nbmax + j;
77 buffer[ib] = CDM[spin][MA_AN][LB_AN][i][j];
78 }
79 }
80 } /* loop of iep */
81
82 /* all-to-all communication */
83 MPI_Allreduce(buffer, buffer, nbuf, MPI_DOUBLE, MPI_SUM, comm);
84
85 /* copy buffer to EXX-DM */
86 for (iep=0; iep<nep; iep++) {
87 GA_AN = ep_atom1[iep]+1;
88 GB_AN = ep_atom2[iep]+1;
89 nb1 = Spe_Total_CNO[WhatSpecies[GA_AN]];
90 nb2 = Spe_Total_CNO[WhatSpecies[GB_AN]];
91
92 /* symmetry breaking factor */
93 if (SpinP_switch==1 && symbrk) {
94 du = InitN_USpin[GA_AN]+ InitN_USpin[GB_AN];
95 dd = InitN_DSpin[GA_AN]+ InitN_DSpin[GB_AN];
96 w = (spin==0) ? (2.0*du/(dd + du)) : (2.0*dd/(dd+du));
97 } else {
98 w = 1.0;
99 }
100
101 for (i=0; i<nb1; i++){
102 for (j=0; j<nb2; j++){
103 ib = iep*nbmax*nbmax + i*nbmax + j;
104 exx_CDM[spin][iep][i][j].r = w*buffer[ib];
105 exx_CDM[spin][iep][i][j].i = 0.0;
106 }
107 }
108 } /* loop of iep */
109 } /* loop of spin */
110
111 free(buffer);
112 }
113
114
115
EXX_Initial_DM(EXX_t * exx,dcomplex **** exx_CDM)116 void EXX_Initial_DM(
117 EXX_t *exx,
118 dcomplex ****exx_CDM
119 )
120 {
121 int i, iep, nep, GA_AN, GB_AN, nb, spin;
122 double den;
123 const int *ep_atom1, *ep_atom2;
124
125 nep = EXX_Number_of_EP(exx);
126 ep_atom1 = EXX_Array_EP_Atom1(exx);
127 ep_atom2 = EXX_Array_EP_Atom2(exx);
128
129 for (iep=0; iep<nep; iep++) {
130 GA_AN = ep_atom1[iep]+1;
131 GB_AN = ep_atom2[iep]+1;
132 if (GA_AN != GB_AN) { continue; }
133 nb = Spe_Total_CNO[WhatSpecies[GA_AN]];
134
135 den = InitN_USpin[GA_AN]/nb;
136 for (i=0; i<nb; i++){
137 exx_CDM[spin][iep][i][i].r = den;
138 exx_CDM[spin][iep][i][i].i = 0.0;
139 }
140
141 if (0<SpinP_switch) {
142 den = InitN_DSpin[GA_AN]/nb;
143 for (i=0; i<nb; i++){
144 exx_CDM[spin][iep][i][i].r = den;
145 exx_CDM[spin][iep][i][i].i = 0.0;
146 }
147 }
148 } /* loop of iep */
149 }
150
151
152
153 #if 0
154 static find_ep(EXX_t *exx, int iatm1, int iatm2, int Rx, int Ry, int Rz)
155 {
156 int nep, iep, nshell, ncd, ia1, ia2, icell, ix, iy, iz;
157 const int *ep_atom1, *ep_atom2, *ep_cell;
158
159 nep = EXX_Number_of_EP(exx);
160 nshell = EXX_Number_of_EP_Shells(exx);
161 ncd = 2*nshell+1;
162
163 ep_atom1 = EXX_Array_EP_Atom1(exx);
164 ep_atom2 = EXX_Array_EP_Atom2(exx);
165 ep_cell = EXX_Array_EP_Cell(exx);
166
167 for (iep=0; iep<nep; iep++) {
168 ia1 = ep_atom1[iep];
169 ia1 = ep_atom2[iep];
170 icell = ep_cell[iep];
171 ix = icell%ncd - nshell;
172 iy = (icell/ncd)%ncd - nshell;
173 iz = (icell/ncd/ncd)%ncd - nshell;
174 if (ia1==iatm1 && ia2== iatm2
175 && ix==Rx && iy==Ry && iz==Rz) { return iep; }
176 }
177
178 return -1;
179 }
180 #endif
181
182
EXX_Debug_Check_DM(EXX_t * exx,dcomplex **** exx_DM,double ***** DM)183 void EXX_Debug_Check_DM(
184 EXX_t *exx,
185 dcomplex ****exx_DM,
186 double *****DM
187 )
188 {
189 int i, j, spin, MA_AN, GA_AN, Anum, LB_AN, GB_AN;
190 int nb1, nb2;
191 int ia1, ia2, l1, l2, l3, iep, Rn;
192 double err, maxerr, x, y;
193
194 maxerr = 0.0;
195
196 for (spin=0; spin<=SpinP_switch; spin++){
197 for (MA_AN=1; MA_AN<=Matomnum; MA_AN++){
198 GA_AN = M2G[MA_AN];
199 nb1 = Spe_Total_CNO[WhatSpecies[GA_AN]];
200 for (LB_AN=0; LB_AN<=FNAN[GA_AN]; LB_AN++){
201 GB_AN = natn[GA_AN][LB_AN];
202 nb2 = Spe_Total_CNO[WhatSpecies[GB_AN]];
203 Rn = ncn[GA_AN][LB_AN];
204
205 l1 = atv_ijk[Rn][1];
206 l2 = atv_ijk[Rn][2];
207 l3 = atv_ijk[Rn][3];
208
209 /*iep = find_ep(exx, GA_AN-1, GB_AN-1, l1, l2, l3);*/
210 iep = EXX_Find_EP(exx, GA_AN-1, GB_AN-1, l1, l2, l3);
211
212 if (-1==iep) { EXX_ERROR("iep not found"); }
213
214 for (i=0; i<nb1; i++){
215 for (j=0; j<nb2; j++){
216 x = DM[spin][MA_AN][LB_AN][i][j];
217 y = exx_DM[spin][iep][i][j].r;
218 err = fabs(x-y);
219 if (err>1e-8) {
220 EXX_WARN("err too large");
221 fprintf(stderr, " DM= %20.12f\n",
222 DM[spin][MA_AN][LB_AN][i][j]);
223 fprintf(stderr, " XDM= %20.12f %20.12f\n",
224 exx_DM[spin][iep][i][j].r, exx_DM[spin][iep][i][j].i);
225 }
226 if (err>maxerr) { maxerr = err; }
227 }
228 }
229 } /* loop of LB_AN */
230 } /* loop of MA_AN */
231 } /* loop of spin */
232
233 if (maxerr < 1e-10) {
234 fprintf(stderr, "DM-Check passed (maxerr=%20.12f)\n", maxerr);
235 } else {
236 fprintf(stderr, "DM-Check failed (maxerr=%20.12f)\n", maxerr);
237 }
238 }
239