source: pcp/src/pseudo.c@ 88e890

Last change on this file since 88e890 was 961b75, checked in by Frederik Heber <heber@…>, 18 years ago

ParseForParameter() has new parameter repetition, which reads the repetition-th occurence of the given keyword
  • Property mode set to 100644
File size: 51.2 KB
RevLine 
[a0bcf1]1/** \file pseudo.c
2 * PseudoPotentials to approximate core densities for plane waves.
3 * Note: Pseudopotentials replace the electron-core-interaction potential \f$V^{El-K}\f$.
4 *
5 * Herein are various functions dealing with the Pseudopotential approximation, which
6 * calculate local and non-local form factors FormFacGauss(), FormFacLocPot(), FormFacNonLocPot(),
7 * energies CalculateNonLocalEnergyNoRT(), CalculateIonNonLocalForce(), CalculateSomeEnergyAndHGNoRT()
8 * and forces CalculateIonLocalForce(), CalculateIonNonLocalForce(),
9 * also small functions for evaluating often used expressions: CalcExpiGR(), dexpiGRpkg(), Get_pkga_MinusG(),
10 * CalculateCDfnl(), CalcSG(), CalcCoreCorrection() - some of these are updated in case of new waves or
11 * finer mesh by calling UpdatePseudoToNewWaves(), ChangePseudoToLevUp(),
12 * for allocation and parameter readin InitPseudoRead() and memory deallocation
13 * RemovePseudoRead().
14 *
15 * Missing so far are: CalculateAddNLPot()
16 *
17 Project: ParallelCarParrinello
18 \author Jan Hamaekers
19 \date 2000
20
21 File: pseudo.c
22 $Id: pseudo.c,v 1.51 2007-03-29 13:38:47 foo Exp $
23*/
24
25#include <stdlib.h>
26#include <stdio.h>
27#include <stddef.h>
28#include <math.h>
29#include <string.h>
30
31
32// use double precision fft when we have it
33#ifdef HAVE_CONFIG_H
34#include <config.h>
35#endif
36
37#ifdef HAVE_DFFTW_H
38#include "dfftw.h"
39#else
40#include "fftw.h"
41#endif
42
43#include "data.h"
44#include "errors.h"
45#include "helpers.h"
46#include "ions.h"
47#include "init.h"
48#include "mymath.h"
49#include "myfft.h"
50#include "pseudo.h"
51#include "run.h"
52
53/** Calculates structure factor \f$S_{I_s} (G)\f$.
54 * \f[
55 * S_{I_s} (G) = \sum_{I_a} \exp(i G\cdot R_{I_s,I_a}) \qquad (4.14)
56 * \f]
57 * \param *P Problem at hand
58 */
59static void CalcSG(struct Problem *P)
60{
61 struct PseudoPot *PP = &P->PP;
62 struct Ions *I = &P->Ion;
63 struct RunStruct *R = &P->R;
64 struct LatticeLevel *Lev0 = R->Lev0;
65 fftw_complex dS;
66 double dSPGRi;
67 int it,iot,g;
68 struct OneGData *G = R->Lev0->GArray;
69 for (it=0; it < I->Max_Types; it++)
70 for (g=0; g < Lev0->MaxG; g++) {
71 c_re(dS) = 0;
72 c_im(dS) = 0;
73 for (iot=0; iot < I->I[it].Max_IonsOfType; iot++) {
74 dSPGRi = RSP3(G[g].G,&I->I[it].R[NDIM*iot]);
75 c_re(PP->ei[it][iot][g]) = cos(dSPGRi);
76 c_im(PP->ei[it][iot][g]) = -sin(dSPGRi);
77 c_re(dS) += c_re(PP->ei[it][iot][g]);
78 c_im(dS) += c_im(PP->ei[it][iot][g]);
79 }
80 c_re(I->I[it].SFactor[g]) = c_re(dS);
81 c_im(I->I[it].SFactor[g]) = c_im(dS);
82 }
83}
84
85
86/** Calculates \f$\exp(-i(G\cdot R)) = \cos(GR) - i\cdot\sin(GR)\f$.
87 * Fills PseudoPot::expiGR array.
88 * \param *P Problem at hand
89 * \sa expiGR() for element retrieval
90 */
91static void CalcExpiGR(struct Problem *P)
92{
93 struct PseudoPot *PP = &P->PP;
94 struct Ions *I = &P->Ion;
95 struct RunStruct *R = &P->R;
96 struct LatticeLevel *LevS = R->LevS;
97 double dSPGRi;
98 int it,iot,g;
99 struct OneGData *G = LevS->GArray;
100 for (it=0; it < I->Max_Types; it++)
101 for (iot=0; iot < I->I[it].Max_IonsOfType; iot++)
102 for (g=0; g < LevS->MaxG; g++) {
103 dSPGRi = RSP3(G[g].G,&I->I[it].R[NDIM*iot]);
104 c_re(PP->expiGR[it][iot][g]) = cos(dSPGRi);
105 c_im(PP->expiGR[it][iot][g]) = -sin(dSPGRi);
106 }
107}
108
109/** Calculates \f$\exp(-i(G)R_{I_s,I_a})\f$.
110 * \param *PP PseudoPot ential structure
111 * \param it ion type \f$I_s\f$
112 * \param iot ion number within certain type (ion of type) \f$I_a\f$
113 * \param g reciprocal grid vector from GArray
114 * \param *res complex return value
115 * \sa CalcExpiGR() fills the array
116 */
117#ifdef HAVE_INLINE
118inline static void expiGR(struct PseudoPot *PP, int it, int iot, int g, fftw_complex *res) {
119#else
120static void expiGR(struct PseudoPot *PP, int it, int iot, int g, fftw_complex *res) {
121#endif
122 c_re(res[0]) = c_re(PP->expiGR[it][iot][g]);
123 c_im(res[0]) = c_im(PP->expiGR[it][iot][g]);
124}
125
126/** Calculates \f$\exp(-i(G)R_{I_s,I_a})\phi_{I_s,l,m}^{ps,nl} (G)\f$.
127 * \param *PP PseudoPot ential structure
128 * \param it ion type \f$I_s\f$
129 * \param iot ion number within certain type (ion of type), \f$I_a\f$
130 * \param ilm m value for ion type
131 * \param g reciprocal grid vector from GArray
132 * \param *res complex return value
133 */
134#ifdef HAVE_INLINE
135inline static void dexpiGRpkg(struct PseudoPot *PP, int it, int iot, int ilm, int g, fftw_complex *res) {
136#else
137static void dexpiGRpkg(struct PseudoPot *PP, int it, int iot, int ilm, int g, fftw_complex *res) {
138#endif
139 expiGR(PP, it, iot, g, res);
140 c_re(res[0]) = c_re(res[0])*PP->phi_ps_nl[it][ilm-1][g];
141 c_im(res[0]) = -c_im(res[0])*PP->phi_ps_nl[it][ilm-1][g];
142}
143
144/** Calculates Gaussian form factor \f$\phi^{gauss}_{I_s}\f$.
145 * \f$I_s\f$ ion types, V volume of cell, \f$r_I^{Gauss}\f$ is defined by the
146 * core charge within via \f$n^{Gauss}\f$, G reciprocal grid vector
147 * \f[
148 * \phi^{gauss}_{I_s} = - \frac{Z_{I_s}}{V} \exp(-\frac{1}{4} (r^{Gauss}_{I_s})^2 |G|^2)
149 * \qquad (4.23)
150 * \f]
151 * \param *P Problem at hand
152 */
153static void FormFacGauss(struct Problem *P)
154{
155 struct PseudoPot *PP = &P->PP;
156 struct Ions *I = &P->Ion;
157 struct Lattice *Lat = &P->Lat;
158 struct RunStruct *R = &P->R;
159 struct LatticeLevel *Lev0 = R->Lev0;
160 double dFac;
161 int it, g;
162 struct OneGData *G = Lev0->GArray;
163 for (it=0; it < I->Max_Types; it++) { // for each all types ...
164 dFac = 0.25*I->I[it].rgauss*I->I[it].rgauss;
165 for (g=0; g < Lev0->MaxG; g++) { // for each reciprocal grid vector ...
166 PP->FacGauss[it][g] = -PP->zval[it]/Lat->Volume*exp(-dFac*G[g].GSq);
167 }
168 }
169}
170
171/** Calculates local PseudoPot form factor \f$\phi_{I_s}^{ps,loc} (G)\f$.
172 * \f$I_s\f$ ion types, \f$I_a\f$ ion number of type, V volume of cell,
173 * \f$r_I^{Gauss}\f$ is defined by the core charge within via \f$n^{Gauss}\f$,
174 * G reciprocal grid vector
175 * \f[
176 * \phi_{I_s}^{ps,loc} (G) = \frac{4}{\pi} \int_0^\infty r^2 j_0(r|G|)
177 * \Bigr ( \nu_{I_s}^{loc} (r) + \frac{Z_{I_s}}{r}erf\bigr ( \frac{r}
178 * {r_{I_s}^{Gauss}} \bigl )\Bigl )
179 * \qquad (4.15)
180 * \f]
181 * \param *P Problem at hand
182 * \note profiling says, sin()/cos() is biggest cpu hog
183 */
184static void FormFacLocPot(struct Problem *P)
185{
186 struct PseudoPot *PP = &P->PP;
187 struct Ions *I = &P->Ion;
188 struct Lattice *Lat = &P->Lat;
189 struct RunStruct *R = &P->R;
190 struct LatticeLevel *Lev0 = R->Lev0;
191 double vv,dint,darg;
192 int I_s, ll, ir, g; // I_s = iontype, ir = ion radial pos, s = counter for g (0 treated specially), g index for Grid vectors GArray
193 struct OneGData *G = Lev0->GArray;
194 for (I_s=0; I_s < I->Max_Types; I_s++) {
195 ll = I->I[I_s].l_loc-1;
196 for (ir = 0; ir < PP->mmax[I_s]; ir++) { //fill integrand1 function array with its values: Z/r * erf(r/R_g)
197 if (fabs(PP->r[I_s][ir]) < MYEPSILON) fprintf(stderr,"FormFacLocPot: PP->r[it][ir] = %lg\n",PP->r[I_s][ir]);
198 if (fabs(I->I[I_s].rgauss) < MYEPSILON) fprintf(stderr,"FormFacLocPot: I->I[it].rgauss = %lg\n",I->I[I_s].rgauss);
199 vv = -PP->zval[I_s]/PP->r[I_s][ir] * derf(PP->r[I_s][ir]/I->I[I_s].rgauss);
200 PP->integrand1[ir] = PP->v_loc[I_s][ll][ir]-vv;
201 }
202 g=0;
203 if (Lev0->GArray[0].GSq == 0.0) { // if grid vectors contain zero vector, treat special
204 for (ir = 0; ir < PP->mmax[I_s]; ir++) { // fill the to be integrated function array with its values
205 PP->integrand[ir] = PP->integrand1[ir]*PP->r[I_s][ir]*PP->r[I_s][ir]*PP->r[I_s][ir];
206 }
207 dint = Simps(PP->mmax[I_s],PP->integrand,PP->clog[I_s]); // do the integration
208 PP->phi_ps_loc[I_s][0] = dint*4.*PI/Lat->Volume;
209 g++;
210 }
211 for (; g < Lev0->MaxG; g++) { // pre-calculate for the rest of all the grid vectors
212 for (ir = 0; ir < PP->mmax[I_s]; ir++) { // r^3 * j_0 * dummy1
213 darg = PP->r[I_s][ir]*sqrt(G[g].GSq);
214 if (fabs(darg) < MYEPSILON) fprintf(stderr,"FormFacLocPot: darg = %lg\n",darg);
215 PP->integrand[ir] = PP->integrand1[ir]*sin(darg)/darg*PP->r[I_s][ir]*PP->r[I_s][ir]*PP->r[I_s][ir];
216 }
217 dint = Simps(PP->mmax[I_s],PP->integrand,PP->clog[I_s]);
218 PP->phi_ps_loc[I_s][g] = dint*4.*PI/Lat->Volume;
219 }
220 }
221}
222
223/** Determines whether its a symmetric or antisymmetric non-local pseudopotential form factor pkga[][value][].
224 * \param value pgga l value
225 * \return 1 - symmetric (0,4,5,6,7,8) , -1 - antisymmetric (1,2,3)
226 */
227inline static int Get_pkga_MinusG(const int value)
228{
229 switch (value) {
230 case 0: return(1);
231 case 1:
232 case 2:
233 case 3: return(-1);
234 case 4:
235 case 5:
236 case 6:
237 case 7:
238 case 8: return(1);
239 default:
240 Error(SomeError,"Get_pkga_MinusG: Unknown Value");
241 }
242 return 0;
243}
244
245/** Calculates non-local PseudoPot form factor \f$\phi_{I_s,l,m}^{ps,nl} (G)\f$ and non-local weight \f$w^{nl}_{I_s,l}\f$.
246 * \f$I_s\f$ ion types, \f$I_a\f$ ion number of type, V volume of cell,
247 * \f$j_l\f$ Bessel function of l-th order,
248 * G reciprocal grid vector
249 * \f[
250 * \phi_{I_s,l,m}^{ps,nl} (G) = \sqrt{\frac{4\pi}{2l+1}} \int_0^\infty r^2 j_l(|G|r)
251 * \Delta \nu^{nl}_{I_s} (r) R_{I_s,l}(r) y_{lm} (\vartheta_G, \varphi_G) dr
252 * \qquad (4.17)
253 * \f]
254 * \f[
255 * w^{nl}_{I_s,l} = \frac{4\pi}{V} (2l+1) \Bigr ( \int^\infty_0 r^2 R_{I_s,l} (r)
256 * \Delta \nu^{nl}_{I_s,l} (r) R_{I_s,l} (r) dr \Bigl )^{-1}
257 * \qquad (4.18)
258 * \f]
259 * \param *P Problem at hand
260 * \note see e.g. MathWorld for the spherical harmonics in cartesian coordinates
261 */
262static void FormFacNonLocPot(struct Problem *P)
263{
264 int I_s,ir,j,ll,il,ml,g,i;
265 double delta_v,arg,ag,fac,cosx,cosy,cosz,dint,sqrt3=sqrt(3.0);
266 struct PseudoPot *PP = &P->PP;
267 struct Ions *I = &P->Ion;
268 struct Lattice *Lat = &P->Lat;
269 struct RunStruct *R = &P->R;
270 struct LatticeLevel *LevS = R->LevS;
271 struct OneGData *G = LevS->GArray;
272 for (I_s=0; I_s < I->Max_Types; I_s++) { // through all ion types
273 ml = -1; // m value
274 for (il=0; il < I->I[I_s].l_max; il++) { // through all possible l values
275 ll = I->I[I_s].l_loc-1; // index for the local potential
276 if (il != ll) {
277 // calculate non-local weights
278 for (ir = 0; ir < PP->mmax[I_s]; ir++) {
279 delta_v = PP->v_loc[I_s][il][ir]-PP->v_loc[I_s][ll][ir]; // v_l - v_l_loc
280 PP->integrand2[I_s][ir] = delta_v*PP->R[I_s][il][ir]*PP->r[I_s][ir]*PP->r[I_s][ir]; // for form factor
281 PP->integrand1[ir] = delta_v*PP->R[I_s][il][ir]*PP->R[I_s][il][ir]*PP->r[I_s][ir]; // for weight
282 }
283 dint = Simps(PP->mmax[I_s],PP->integrand1,PP->clog[I_s]);
284 if (fabs(dint) < MYEPSILON) fprintf(stderr,"FormFacNonLocPot: dint[%d] = %lg\n",I_s,dint);
285 if (il+1 < 4) {
286 for(j =1; j <= 2*(il+1)-1; j++) {
287 PP->wNonLoc[I_s][ml+j]=(2.0*(il+1)-1.0)*4.*PI/Lat->Volume/dint; // store weight
288 }
289 }
290 // calculate non-local form factors
291 for (g=0; g < LevS->MaxG; g++) {
292 arg = sqrt(G[g].GSq);
293 switch (il) { // switch on the order of the needed bessel function, note: only implemented till second order
294 case 0:
295 if (arg == 0.0) { // arg = 0 ...
296 for (ir=0; ir < PP->mmax[I_s]; ir++) {
297 PP->integrand1[ir] = PP->integrand2[I_s][ir];
298 }
299 } else { // ... or evaluate J_0
300 for (ir=0; ir < PP->mmax[I_s]; ir++) {
301 fac=arg*PP->r[I_s][ir];
302 if (fabs(fac) < MYEPSILON) fprintf(stderr,"FormFacNonLocPot: fac[%d][%d] = %lg\n",I_s,ir,fac);
303 fac=sin(fac)/fac; // CHECKED: j_0 (x)
304 PP->integrand1[ir] = fac*PP->integrand2[I_s][ir];
305 }
306 }
307 dint = Simps(PP->mmax[I_s],PP->integrand1,PP->clog[I_s]);
308 PP->phi_ps_nl[I_s][ml+1][g]=dint; // m = 0
309 break;
310 case 1:
311 if (arg == 0.0) { // arg = 0 ...
312 for (j=1; j <= 3; j++) {
313 PP->phi_ps_nl[I_s][ml+j][g]=0.0;
314 }
315 } else { // ... or evaluate J_1
316 for (ir=0; ir < PP->mmax[I_s]; ir++) {
317 fac=arg*PP->r[I_s][ir];
318 if (fabs(fac) < MYEPSILON) fprintf(stderr,"FormFacNonLocPot: fac[%d][%d] = %lg\n",I_s,ir,fac);
319 fac=(sin(fac)/fac-cos(fac))/fac; // CHECKED: j_1(x)
320 PP->integrand1[ir] = fac*PP->integrand2[I_s][ir];
321 }
322 dint = Simps(PP->mmax[I_s],PP->integrand1,PP->clog[I_s]);
323 PP->phi_ps_nl[I_s][ml+1][g]=dint*G[g].G[0]/arg; // m = -1
324 PP->phi_ps_nl[I_s][ml+2][g]=dint*G[g].G[1]/arg; // m = +1
325 PP->phi_ps_nl[I_s][ml+3][g]=dint*G[g].G[2]/arg; // m = 0
326 }
327 break;
328 case 2:
329 if (arg == 0.0) { // arg = 0 ...
330 for (j=1; j <= 5; j++) {
331 PP->phi_ps_nl[I_s][ml+j][g]=0.0;
332 }
333 } else { // ... or evaluate J_2
334 for (ir=0; ir < PP->mmax[I_s]; ir++) {
335 ag=arg*PP->r[I_s][ir];
336 if (fabs(ag) < MYEPSILON) fprintf(stderr,"FormFacNonLocPot: ag[%d][%d] = %lg\n",I_s,ir,ag);
337 fac=(3.0/(ag*ag)-1.0)*sin(ag)-3.0/ag*cos(ag);
338 fac=fac/ag; // CHECKED: j_2(x)
339 PP->integrand1[ir] = fac*PP->integrand2[I_s][ir];
340 }
341 dint = Simps(PP->mmax[I_s],PP->integrand1,PP->clog[I_s]);
342 if (fabs(arg) < MYEPSILON) fprintf(stderr,"FormFacNonLocPot: arg = %lg\n",arg);
343 cosx = G[g].G[0]/arg;
344 cosy = G[g].G[1]/arg;
345 cosz = G[g].G[2]/arg;
346 PP->phi_ps_nl[I_s][ml+1][g]=dint*0.5*(3.0*cosz*cosz-1.0); // m = 0
347 PP->phi_ps_nl[I_s][ml+2][g]=dint*sqrt3*cosz*cosx; // m = +1
348 PP->phi_ps_nl[I_s][ml+3][g]=dint*sqrt3*cosz*cosy; // m = -1
349 PP->phi_ps_nl[I_s][ml+4][g]=dint*sqrt3*cosx*cosy; // m = -2
350 PP->phi_ps_nl[I_s][ml+5][g]=dint*sqrt3/2.0*(cosx*cosx-cosy*cosy); // m = +2
351 }
352 break;
353 default:
354 Error(SomeError, "FormFacNonLocPot: Order of sperhical Bessel function not implemented!");
355 break;
356 }
357 // phi_ps_nl_a removed, not needed, just a mem hog
358// for (j=1; j<=2*(il+1)-1;j++) {
359// PP->phi_ps_nl_a[I_s][ml+j][g] = PP->phi_ps_nl[I_s][ml+j][g];
360// }
361 }
362 ml += (2*(il+1)-1);
363 }
364 }
365 for (i=0; i < I->I[I_s].l_max*I->I[I_s].l_max-2*I->I[I_s].l_loc+1; i++)
366 if(P->Call.out[ReadOut])
367 fprintf(stderr,"wnl:(%i,%i) = %g\n",I_s,i,PP->wNonLoc[I_s][i]);
368 }
369}
370
371/** Calculates partial core form factors \f$\phi^{pc}_{I_s}(G)\f$ and \f$\hat{n}^{pc} (G)\f$.
372 * \f[
373 * \phi^{pc}_{I_s}(G) = \frac{4\pi}{V} \int_0^\infty r^2 j_0(|G|r) n^{pc}_{I_s} (r) dr \qquad (4.21)
374 * \f]
375 * \f[
376 * \hat{n}^{pc} = \sum_{I_s} S(G) \phi^{pc}_{I_s}(G)
377 * \f]
378 * \param *P Problem at hand
379 * \param Create 1 - calculate form factors, 0 - don't
380 */
381static void CalcCoreCorrection(struct Problem *P, const int Create)
382{
383 int i, s = 0;
384 int it,g,Index,ir;
385 double arg;
386 struct PseudoPot *PP = &P->PP;
387 struct Ions *I = &P->Ion;
388 struct Lattice *Lat = &P->Lat;
389 struct RunStruct *R = &P->R;
390 struct LatticeLevel *Lev0 = R->Lev0;
391 struct OneGData *G = Lev0->GArray;
392 struct Density *Dens = Lev0->Dens;
393 struct fft_plan_3d *plan = Lat->plan;
394 fftw_complex *work = (fftw_complex *)Dens->DensityCArray[TempDensity];
395 fftw_complex *CoreWave = (fftw_complex *)Dens->DensityArray[CoreWaveDensity];
396 // if desired the form factors are calculated
397 if (Create) {
398 for (it=0; it < I->Max_Types; it++) {
399 if (I->I[it].corecorr == CoreCorrected) {
400 s=0;
401 if (Lev0->GArray[0].GSq == 0.0) {
402 for (ir=0; ir < PP->mmax[it]; ir++) {
403 PP->integrand[ir] = 4.*PI*PP->corewave[it][ir]*PP->r[it][ir]*PP->r[it][ir]*PP->r[it][ir]/Lat->Volume;
404 }
405 PP->formfCore[it][0] = Simps(PP->mmax[it], PP->integrand, PP->clog[it]);
406 s++;
407 }
408 for (g=s; g < Lev0->MaxG; g++) {
409 for (ir=0; ir < PP->mmax[it]; ir++) {
410 arg=PP->r[it][ir]*sqrt(G[g].GSq);
411 if (fabs(arg) < MYEPSILON) fprintf(stderr,"CalcCoreCorrection: arg = %lg\n",arg);
412 PP->integrand1[ir] = PP->integrand[ir]*sin(arg)/arg;
413 }
414 PP->formfCore[it][g] = Simps(PP->mmax[it], PP->integrand1, PP->clog[it]);
415 }
416 }
417 }
418 }
419 // here the density is evaluated, note again that it has to be enfolded due to the finer resolution
420 SetArrayToDouble0((double *)CoreWave,2*Dens->TotalSize);
421 for (g=0; g < Lev0->MaxG; g++) {
422 Index = Lev0->GArray[g].Index;
423 for (it = 0; it < I->Max_Types; it++) {
424 if (I->I[it].corecorr == CoreCorrected) {
425 c_re(CoreWave[Index]) += PP->formfCore[it][g]*c_re(I->I[it].SFactor[g]);
426 c_im(CoreWave[Index]) += PP->formfCore[it][g]*c_im(I->I[it].SFactor[g]);
427 }
428 }
429 }
430 for (i=0; i<Lev0->MaxDoubleG; i++) {
431 CoreWave[Lev0->DoubleG[2*i+1]].re = CoreWave[Lev0->DoubleG[2*i]].re;
432 CoreWave[Lev0->DoubleG[2*i+1]].im = -CoreWave[Lev0->DoubleG[2*i]].im;
433 }
434 fft_3d_complex_to_real(plan, Lev0->LevelNo, FFTNF1, CoreWave, work);
435}
436
437/** Reads data for Pseudopotentials.
438 * PseudoPot structure is allocated and intialized, depending on the given Z values in the
439 * config file, corresponding pseudopotential files are read (generated from FHMD98 package),
440 * basically calls ChangePseudoToLevUp() (only it's implemented by hand).
441 * \param *P Problem at hand
442 * \param *pseudopot_path Path to pseudopot files
443 * \sa ParseForParameter() for parsing, RemovePseudoRead() for deallocation
444 */
445void InitPseudoRead(struct Problem *P, char *pseudopot_path)
446{
447 char cpiInputFileName[255];
448 FILE *cpiInputFile;
449 int i,it,ib,il,m,j,ia;
450 int count = 0;
451 double dummycorewave = 0., dummyr = 0.;
452 struct PseudoPot *PP = &P->PP;
453 struct Ions *I = &P->Ion;
454 struct Lattice *Lat = &P->Lat;
455 struct RunStruct *R = &P->R;
456 struct LatticeLevel *ILev0 = R->InitLev0;
457 struct LatticeLevel *ILevS = R->InitLevS;
458 PP->Mmax = 0;
459 PP->core = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
460 PP->rc = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
461 PP->rcl = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
462 PP->al = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
463 PP->bl = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
464 PP->mmax = (int *) Malloc(sizeof(int)*I->Max_Types, "InitPseudoRead: ");
465 PP->clog = (double *) Malloc(sizeof(double)*I->Max_Types, "InitPseudoRead: ");
466 PP->R = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
467 PP->r = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
468 PP->integrand2 = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
469 PP->v_loc = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
470 PP->zval = (double *)Malloc(sizeof(double)*I->Max_Types, "InitPseudoRead: ");
471 PP->nang = (int *) Malloc(sizeof(int)*I->Max_Types, "InitPseudoRead: ");
472 PP->FacGauss = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
473 PP->phi_ps_loc = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
474 PP->wNonLoc = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
475 //PP->phi_ps_nl = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
476 PP->phi_ps_nl = (double ***) Malloc(sizeof(double**)*I->Max_Types, "InitPseudoRead: ");
477 PP->lm_end = (int *) Malloc(sizeof(int)*I->Max_Types, "InitPseudoRead: ");
478 PP->ei = (fftw_complex ***) Malloc(sizeof(fftw_complex**)*I->Max_Types, "InitPseudoRead: ");
479 PP->expiGR = (fftw_complex ***) Malloc(sizeof(fftw_complex**)*I->Max_Types, "InitPseudoRead: ");
480 PP->VCoulombc = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILev0->MaxG, "InitPseudoRead: ");
481 PP->corewave = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
482 PP->formfCore = (double **) Malloc(sizeof(double*)*I->Max_Types, "InitPseudoRead: ");
483 PP->lm_endmax = 0;
484 PP->corecorr = NotCoreCorrected;
485 for (it = 0; it < I->Max_Types; it++) {
486 PP->lm_end[it] = I->I[it].l_max*I->I[it].l_max+1-2*I->I[it].l_loc;
487 if (PP->lm_endmax < PP->lm_end[it]) PP->lm_endmax = PP->lm_end[it];
488 //PP->phi_ps_nl[it] = (double **) Malloc(sizeof(double*)*PP->lm_end[it], "InitPseudoRead: ");
489 PP->phi_ps_nl[it] = (double **) Malloc(sizeof(double*)*PP->lm_end[it], "InitPseudoRead: ");
490 PP->ei[it] = (fftw_complex **) Malloc(sizeof(fftw_complex*)*I->I[it].Max_IonsOfType, "InitPseudoRead: ");
491 PP->expiGR[it] = (fftw_complex **) Malloc(sizeof(fftw_complex*)*I->I[it].Max_IonsOfType, "InitPseudoRead: ");
492 for (ia=0;ia<I->I[it].Max_IonsOfType;ia++) {
493 PP->ei[it][ia] = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILev0->MaxG, "InitPseudoRead: ");
494 PP->expiGR[it][ia] = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILevS->MaxG, "InitPseudoRead: ");
495 }
496 for (il=0;il<PP->lm_end[it];il++) {
497 //PP->phi_ps_nl[it][il] = (double *) Malloc(sizeof(double)*ILevS->MaxG, "InitPseudoRead: ");
498 PP->phi_ps_nl[it][il] = (double *) Malloc(sizeof(double)*ILevS->MaxG, "InitPseudoRead: ");
499 }
500 PP->wNonLoc[it] = (double *) Malloc(sizeof(double)*PP->lm_end[it], "InitPseudoRead: ");
501 sprintf(cpiInputFileName,"%s/pseudo.%02i",pseudopot_path,I->I[it].Z);
502 cpiInputFile = fopen(cpiInputFileName,"r");
503 if (cpiInputFile == NULL)
504 Error(SomeError, "Cant't find pseudopot path or file");
505 else
506 if (P->Call.out[ReadOut]) fprintf(stderr,"%s was found\n",cpiInputFileName);
507 int zeile = 1;
[961b75]508 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->zval[it], 1, critical);
509 ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, int_type, &PP->nang[it], 1, critical);
[a0bcf1]510 I->TotalZval += PP->zval[it]*(I->I[it].Max_IonsOfType);
511 PP->core[it] = (double *) Malloc(sizeof(double)*2, "InitPseudoRead: ");
512 PP->rc[it] = (double *) Malloc(sizeof(double)*2, "InitPseudoRead: ");
[961b75]513 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->core[it][0], 1, critical);
514 ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->rc[it][0], 1, critical);
515 ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &PP->core[it][1], 1, critical);
516 ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &PP->rc[it][1], 1, critical);
[a0bcf1]517 PP->rcl[it] = (double **) Malloc(sizeof(double*)*3, "InitPseudoRead: ");
518 PP->al[it] = (double **) Malloc(sizeof(double*)*3, "InitPseudoRead: ");
519 PP->bl[it] = (double **) Malloc(sizeof(double*)*3, "InitPseudoRead: ");
520 PP->FacGauss[it] = (double *) Malloc(sizeof(double)*ILev0->MaxG, "InitPseudoRead: ");
521 PP->phi_ps_loc[it] = (double *) Malloc(sizeof(double)*ILev0->MaxG, "InitPseudoRead: ");
522 I->I[it].SFactor = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILev0->MaxG, "InitPseudoRead: ");
523 for (il=0; il < 3; il++) {
524 PP->rcl[it][il] = (double *) Malloc(sizeof(double)*3, "InitPseudoRead: ");
525 PP->al[it][il] = (double *) Malloc(sizeof(double)*3, "InitPseudoRead: ");
526 PP->bl[it][il] = (double *) Malloc(sizeof(double)*3, "InitPseudoRead: ");
527 for (ib = 0; ib < 3; ib++) {
[961b75]528 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->rcl[it][il][ib], 1, critical);
529 ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->al[it][il][ib], 1, critical);
530 ParseForParameter(0,cpiInputFile, NULL, 1, 3, zeile, double_type, &PP->bl[it][il][ib], 1, critical);
[a0bcf1]531 if (PP->rcl[it][il][ib] < MYEPSILON)
532 PP->rcl[it][il][ib] = 0.0;
533 else
534 PP->rcl[it][il][ib] = 1./sqrt(PP->rcl[it][il][ib]);
535 }
536 }
[961b75]537 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &PP->mmax[it], 1, critical);
538 ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, double_type, &PP->clog[it], 1, critical);
[a0bcf1]539 if (PP->mmax[it] > PP->Mmax) PP->Mmax = PP->mmax[it];
540 PP->clog[it] = log(PP->clog[it]);
541 PP->r[it] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
542 PP->integrand2[it] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
543 PP->R[it] = (double **) Malloc(sizeof(double*)*PP->nang[it], "InitPseudoRead: ");
544 PP->v_loc[it] = (double **) Malloc(sizeof(double*)*PP->nang[it], "InitPseudoRead: ");
545 for (il=0; il < PP->nang[it]; il++) {
546 PP->R[it][il] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
547 PP->v_loc[it][il] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
548 for (j=0;j< PP->mmax[it];j++) {
[961b75]549 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &m, 1, critical);
550 ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->r[it][j], 1, critical);
551 ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &PP->R[it][il][j], 1, critical);
552 ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &PP->v_loc[it][il][j], 1, critical);
[a0bcf1]553 }
554 if (il < PP->nang[it]-1) {
[961b75]555 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &PP->mmax[it], 1, critical);
556 ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, double_type, &PP->clog[it], 1, critical);
[a0bcf1]557 PP->clog[it] = log(PP->clog[it]);
558 }
559 }
560 I->I[it].corecorr = NotCoreCorrected;
561 count = 0;
[961b75]562 count += ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &dummyr, 1, optional);
563 count += ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &dummycorewave, 1, optional);
[a0bcf1]564 //fprintf(stderr,"(%i) %lg %lg\n",P->Par.me,dummyr,dummycorewave);
565 if (count == 2) {
566 PP->r[it][0] = dummyr;
567 PP->corewave[it] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
568 PP->formfCore[it] = (double *) Malloc(sizeof(double)*ILev0->MaxG, "InitPseudoRead: ");
569 I->I[it].corecorr = CoreCorrected;
570 PP->corecorr = CoreCorrected;
571 PP->corewave[it][0] = dummycorewave/(4.*PI);
[961b75]572 ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &dummyr, 1, critical);
573 ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &dummycorewave, 1, critical);
[a0bcf1]574 for (j=1;j < PP->mmax[it]; j++) {
[961b75]575 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->r[it][j], 1, critical);
576 ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->corewave[it][j], 1, critical);
577 ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &dummyr, 1, critical);
578 ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &dummycorewave, 1, critical);
[a0bcf1]579 PP->corewave[it][j] /= (4.*PI);
580 }
581 } else {
582 PP->corewave[it] = NULL;
583 PP->formfCore[it] = NULL;
584 }
585 fclose(cpiInputFile);
586 }
587 PP->fnl = (fftw_complex ****) Malloc(sizeof(fftw_complex***)*(Lat->Psi.LocalNo+1), "InitPseudoRead: ");
588 for (i=0; i< Lat->Psi.LocalNo+1; i++) {
589 PP->fnl[i] = (fftw_complex ***) Malloc(sizeof(fftw_complex**)*I->Max_Types, "InitPseudoRead: ");
590 for (it=0; it < I->Max_Types; it++) {
591 PP->fnl[i][it] = (fftw_complex **) Malloc(sizeof(fftw_complex*)*PP->lm_end[it], "InitPseudoRead: ");
592 for (il =0; il < PP->lm_end[it]; il++)
593 PP->fnl[i][it][il] = (fftw_complex *) Malloc(sizeof(fftw_complex)*I->I[it].Max_IonsOfType, "InitPseudoRead: ");
594 }
595 }
596 if (fabs(I->TotalZval - P->Lat.Psi.NIDensity[Occupied]) >= MYEPSILON) {
597 if (P->Par.me_comm_ST == 0) // instead of P->Par.me, thus differing charge density output also for SpinUp/Down from both processes!
598 fprintf(stderr, "TotalZval %i != NIDensity[0] %g eps (%g >= %g)\n",I->TotalZval, P->Lat.Psi.NIDensity[Occupied], fabs(I->TotalZval - P->Lat.Psi.NIDensity[Occupied]),MYEPSILON);
599 Error(SomeError,"Readparameters: charged system is not implemented yet");
600 }
601 PP->CDfnl = PP->fnl[Lat->Psi.LocalNo];
602 PP->dfnl = (fftw_complex *) Malloc(sizeof(fftw_complex)*I->Max_Max_IonsOfType, "InitPseudoRead: ");
603 PP->rr = (double *) Malloc(sizeof(double)*PP->lm_endmax, "InitPseudoRead: ");
604 PP->t = (fftw_complex *) Malloc(sizeof(fftw_complex)*PP->lm_endmax, "InitPseudoRead: ");
605 PP->integrand = (double *) Malloc(sizeof(double)*PP->Mmax, "InitPseudoRead: ");
606 PP->integrand1 = (double *) Malloc(sizeof(double)*PP->Mmax, "InitPseudoRead: ");
607 for (it=0;it < I->Max_Types; it++) {
608 if (I->I[it].corecorr == CoreCorrected) {
609 for (j=0; j < PP->mmax[it]; j++) {
610 PP->integrand[j] = 4.*PI*PP->corewave[it][j]*PP->r[it][j]*PP->r[it][j]*PP->r[it][j];
611 }
612 if(P->Call.out[ReadOut]) fprintf(stderr,"Ion[%i] is core corrected: core charge = %g\n", it, Simps(PP->mmax[it], PP->integrand, PP->clog[it]));
613 }
614 }
615 PP->fac1sqrt2PI = 1./sqrt(2.*PI);
616 PP->fac1sqrtPI = 1./sqrt(PI);
617 SpeedMeasure(P, InitSimTime, StartTimeDo);
618 SpeedMeasure(P, InitLocTime, StartTimeDo);
619 CalcSG(P);
620 CalcExpiGR(P);
621 FormFacGauss(P);
622 FormFacLocPot(P);
623 SpeedMeasure(P, InitLocTime, StopTimeDo);
624 SpeedMeasure(P, InitNonLocTime, StartTimeDo);
625 FormFacNonLocPot(P);
626 SpeedMeasure(P, InitNonLocTime, StopTimeDo);
627 SpeedMeasure(P, InitLocTime, StartTimeDo);
628 CalcCoreCorrection(P, 1);
629 SpeedMeasure(P, InitLocTime, StopTimeDo);
630 SpeedMeasure(P, InitSimTime, StopTimeDo);
631}
632
633/** Updates Pseudopotentials due to change of mesh width.
634 * Calls CalcSG(), CalcExpiGR(), recalculates the form factors by calling
635 * FormFacGauss(), FormFacLocPot(), FormFacLonLocPot() and finally
636 * CalcCoreCorrection(). All speed-measured in InitLocTime respectively
637 * InitNonLocTime.
638 * \param *P Problem at hand
639 */
640inline void ChangePseudoToLevUp(struct Problem *P)
641{
642 SpeedMeasure(P, InitLocTime, StartTimeDo);
643 CalcSG(P);
644 CalcExpiGR(P);
645 FormFacGauss(P);
646 FormFacLocPot(P);
647 SpeedMeasure(P, InitLocTime, StopTimeDo);
648 SpeedMeasure(P, InitNonLocTime, StartTimeDo);
649 FormFacNonLocPot(P);
650 SpeedMeasure(P, InitNonLocTime, StopTimeDo);
651 SpeedMeasure(P, InitLocTime, StartTimeDo);
652 CalcCoreCorrection(P, 1);
653 SpeedMeasure(P, InitLocTime, StopTimeDo);
654}
655
656/** Updates Pseudopotentials due to the newly minimized wave functions.
657 * Calls CalcSG(), CalcExpiGR() and CalcCoreCorrection().
658 * \param *P Problem at hand
659 */
660inline void UpdatePseudoToNewWaves(struct Problem *P)
661{
662 SpeedMeasure(P, InitLocTime, StartTimeDo);
663 CalcSG(P);
664 CalcExpiGR(P);
665 CalcCoreCorrection(P, 0);
666 SpeedMeasure(P, InitLocTime, StopTimeDo);
667}
668
669/** Frees memory allocated from Readin.
670 * All memory is free'd that was allocated in InitPseudoRead()
671 * \param *P Problem at hand
672 * \sa InitPseudoRead()
673 */
674void RemovePseudoRead(struct Problem *P)
675{
676 int i,it,il,ia;
677 struct PseudoPot *PP = &P->PP;
678 struct Ions *I = &P->Ion;
679 struct Lattice *Lat = &P->Lat;
[64fa9e]680 Free(PP->integrand1, "RemovePseudoRead: PP->integrand1");
681 Free(PP->integrand, "RemovePseudoRead: PP->integrand");
682 Free(PP->dfnl, "RemovePseudoRead: PP->dfnl");
683 Free(PP->rr, "RemovePseudoRead: PP->rr");
684 Free(PP->t, "RemovePseudoRead: PP->t");
[a0bcf1]685 for (i=0; i< Lat->Psi.LocalNo+1; i++) {
686 for (it=0; it < I->Max_Types; it++) {
687 for (il =0; il < PP->lm_end[it]; il++)
[64fa9e]688 Free(PP->fnl[i][it][il], "RemovePseudoRead: PP->fnl[i][it][il]");
689 Free(PP->fnl[i][it], "RemovePseudoRead: PP->fnl[i][it]");
[a0bcf1]690 }
[64fa9e]691 Free(PP->fnl[i], "RemovePseudoRead: PP->fnl[i]");
[a0bcf1]692 }
693 for (it=0; it < I->Max_Types; it++) {
694 if (I->I[it].corecorr == CoreCorrected) {
[64fa9e]695 Free(PP->corewave[it], "RemovePseudoRead: PP->corewave[it]");
696 Free(PP->formfCore[it], "RemovePseudoRead: PP->formfCore[it]");
[a0bcf1]697 }
698 }
699 for (it=0; it < I->Max_Types; it++) {
700 for (il=0; il < PP->nang[it]; il++) {
[64fa9e]701 Free(PP->R[it][il], "RemovePseudoRead: PP->R[it][il]");
702 Free(PP->v_loc[it][il], "RemovePseudoRead: PP->v_loc[it][il]");
[a0bcf1]703 }
704 for (il=0; il < 3; il++) {
[64fa9e]705 Free(PP->rcl[it][il], "RemovePseudoRead: PP->rcl[it][il]");
706 Free(PP->al[it][il], "RemovePseudoRead: PP->al[it][il]");
707 Free(PP->bl[it][il], "RemovePseudoRead: PP->bl[it][il]");
[a0bcf1]708 }
709 for (il=0;il<PP->lm_end[it];il++) {
[64fa9e]710 //Free(PP->phi_ps_nl[it][il], "RemovePseudoRead: PP->phi_ps_nl[it][il]");
711 Free(PP->phi_ps_nl[it][il], "RemovePseudoRead: PP->phi_ps_nl[it][il]");
[a0bcf1]712 }
713 for (ia=0;ia<I->I[it].Max_IonsOfType;ia++) {
[64fa9e]714 Free(PP->ei[it][ia], "RemovePseudoRead: PP->ei[it][ia]");
715 Free(PP->expiGR[it][ia], "RemovePseudoRead: PP->expiGR[it][ia]");
[a0bcf1]716 }
[64fa9e]717 Free(PP->ei[it], "RemovePseudoRead: PP->ei[it]");
718 Free(PP->expiGR[it], "RemovePseudoRead: PP->expiGR[it]");
719 //Free(PP->phi_ps_nl[it], "RemovePseudoRead: PP->phi_ps_nl[it]");
720 Free(PP->phi_ps_nl[it], "RemovePseudoRead: PP->phi_ps_nl[it]");
721 Free(PP->R[it], "RemovePseudoRead: PP->R[it]");
722 Free(PP->v_loc[it], "RemovePseudoRead: PP->v_loc[it]");
723 Free(PP->r[it], "RemovePseudoRead: PP->r[it]");
724 Free(PP->integrand2[it], "RemovePseudoRead: PP->integrand2[it]");
725 Free(PP->rcl[it], "RemovePseudoRead: PP->rcl[it]");
726 Free(PP->al[it], "RemovePseudoRead: PP->al[it]");
727 Free(PP->bl[it], "RemovePseudoRead: PP->bl[it]");
728 Free(PP->FacGauss[it], "RemovePseudoRead: PP->FacGauss[it]");
729 Free(PP->phi_ps_loc[it], "RemovePseudoRead: PP->phi_ps_loc[it]");
730 Free(PP->core[it], "RemovePseudoRead: PP->core[it]");
731 Free(PP->rc[it], "RemovePseudoRead: PP->rc[it]");
732 Free(PP->wNonLoc[it], "RemovePseudoRead: PP->wNonLoc[it]");
[a0bcf1]733 }
[64fa9e]734 //Free(PP->phi_ps_nl, "RemovePseudoRead: PP->phi_ps_nl");
735 Free(PP->phi_ps_nl, "RemovePseudoRead: PP->phi_ps_nl");
736 Free(PP->R, "RemovePseudoRead: PP->R");
737 Free(PP->v_loc, "RemovePseudoRead: PP->v_loc");
738 Free(PP->r, "RemovePseudoRead: PP->r");
739 Free(PP->integrand2, "RemovePseudoRead: PP->integrand2");
740 Free(PP->rcl, "RemovePseudoRead: PP->rcl");
741 Free(PP->al, "RemovePseudoRead: PP->al");
742 Free(PP->bl, "RemovePseudoRead: PP->bl");
743 Free(PP->FacGauss, "RemovePseudoRead: PP->FacGauss");
744 Free(PP->phi_ps_loc, "RemovePseudoRead: PP->phi_ps_loc");
745 Free(PP->core, "RemovePseudoRead: PP->core");
746 Free(PP->rc, "RemovePseudoRead: PP->rc");
747 Free(PP->wNonLoc, "RemovePseudoRead: PP->wNonLoc");
748 Free(PP->ei, "RemovePseudoRead: PP->ei");
749 Free(PP->expiGR, "RemovePseudoRead: PP->expiGR");
750 Free(PP->corewave, "RemovePseudoRead: PP->corewave");
751 Free(PP->formfCore, "RemovePseudoRead: PP->formfCore");
752 Free(PP->fnl, "RemovePseudoRead: PP->fnl");
753 Free(PP->VCoulombc, "RemovePseudoRead: PP->VCoulombc");
754 Free(PP->lm_end, "RemovePseudoRead: PP->lm_end");
755 Free(PP->zval, "RemovePseudoRead: PP->zval");
756 Free(PP->nang, "RemovePseudoRead: PP->nang");
757 Free(PP->mmax, "RemovePseudoRead: PP->mmax");
758 Free(PP->clog, "RemovePseudoRead: PP->clog");
[a0bcf1]759}
760
761/* Calculate Energy and Forces */
762
763/** Calculates Gauss (self) energy term of ions \f$E_{self}\f$.
764 * \f[
765 * E_{self} = \sum_{I_s,I_a} \frac{1}{\sqrt{2\pi}} \frac{Z_{I_s}^2}{r_{I_s}^{Gauss}} \qquad (4.10 last part)
766 * \f]
767 * \param *P Problem at hand
768 * \param *PP PseudoPot ential structure
769 * \param *I Ions structure
770 */
771void CalculateGaussSelfEnergyNoRT(struct Problem *P, struct PseudoPot *PP, struct Ions *I)
772{
773 double SumE = 0.0;
774 int it;
775 for (it=0; it < I->Max_Types; it++) {
776 if (I->I[it].rgauss < MYEPSILON) fprintf(stderr,"CalculateGaussSelfEnergyNoRT: I->I[it].rgauss = %lg\n",I->I[it].rgauss);
777 SumE += PP->zval[it]*PP->zval[it]/I->I[it].rgauss*I->I[it].Max_IonsOfType;
778 }
779 P->Lat.E->AllTotalIonsEnergy[GaussSelfEnergy] = SumE*PP->fac1sqrt2PI;
780}
781
782/** Calculates non local pseudopotential energy \f$E_{ps,nl,i}\f$.
783 * First, calculates non-local form factors
784 * \f[
785 * f^{nl}_{i,I_s,I_a,l,m} = \sum_G \exp(-i(G R_{I_s,I_a})
786 * \phi_{I_s,l,m}^{ps,nl} (G) c_{i,G}
787 * \f]
788 * and afterwards evaluates with these
789 * \f[
790 * E_{ps,nl,i} = f_{i} \sum_{I_s,I_a} \sum_{l,m} w^{nl}_{I_s,l} | f^{nl}_{i,I_s,I_a,l,m} |^2 \qquad (4.16)
791 * \f]
792 * \param *P Problem at hand
793 * \param i Energy of i-th Psi
794 */
795void CalculateNonLocalEnergyNoRT(struct Problem *P, const int i)
796{
797 struct PseudoPot *PP = &P->PP;
798 struct Ions *I = &P->Ion;
799 struct Lattice *Lat = &P->Lat;
800 struct RunStruct *R = &P->R;
801 struct LatticeLevel *LevS = R->LevS;
802 int it,iot,ilm,ig,s=0;
803 double sf,enl;//,enlm;
804 fftw_complex PPdexpiGRpkg;
805 fftw_complex *dfnl = PP->dfnl;
806 fftw_complex *LocalPsi = LevS->LPsi->LocalPsi[i]; // i-th wave function coefficients
807 const double PsiFactor = Lat->Psi.LocalPsiStatus[i].PsiFactor;
808 int MinusGFactor;
809 enl=0.0;
810 for (it=0; it < I->Max_Types; it++) { // through all ion types
811 for (ilm=1; ilm <=PP->lm_end[it]; ilm++) { // over all possible m values
812 MinusGFactor = Get_pkga_MinusG(ilm-1);
813 SetArrayToDouble0((double *)dfnl, 2*I->I[it].Max_IonsOfType);
814 for (iot=0; iot < I->I[it].Max_IonsOfType; iot++) { // for each ion per type
815 s=0;
816 if (LevS->GArray[0].GSq == 0.0) {
817 dexpiGRpkg(PP, it, iot, ilm, 0, &PPdexpiGRpkg);
818 c_re(dfnl[iot]) += (c_re(LocalPsi[0])*c_re(PPdexpiGRpkg)-c_im(LocalPsi[0])*c_im(PPdexpiGRpkg));
819 c_im(dfnl[iot]) += (c_re(LocalPsi[0])*c_im(PPdexpiGRpkg)+c_im(LocalPsi[0])*c_re(PPdexpiGRpkg));
820 s++;
821 }
822 for (ig=s; ig < LevS->MaxG; ig++) {
823 dexpiGRpkg(PP, it, iot, ilm, ig, &PPdexpiGRpkg);
824 c_re(dfnl[iot]) += (c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg));
825 c_im(dfnl[iot]) += (c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg));
826 /* Minus G - due to inherent symmetry at gamma point! */
827 c_re(dfnl[iot]) += MinusGFactor*(c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg));
828 c_im(dfnl[iot]) -= MinusGFactor*(c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg));
829 }
830 }
831 // Allreduce as the coefficients are spread over many processes
832 MPI_Allreduce( dfnl, PP->fnl[i][it][ilm-1], 2*I->I[it].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
833 }
834 }
835 for (it=0; it < I->Max_Types; it++) {
836 for (ilm=1; ilm <=PP->lm_end[it]; ilm++) {
837 for (iot=0; iot < I->I[it].Max_IonsOfType; iot++) {
838 //enlm = 0.0;
839 sf = c_re(PP->fnl[i][it][ilm-1][iot])*c_re(PP->fnl[i][it][ilm-1][iot])+c_im(PP->fnl[i][it][ilm-1][iot])*c_im(PP->fnl[i][it][ilm-1][iot]);
840 //enlm = (PsiFactor*sf);
841 enl += PsiFactor*sf*PP->wNonLoc[it][ilm-1];
842 }
843 }
844 }
845 Lat->Energy[R->CurrentMin].PsiEnergy[NonLocalEnergy][i] = enl;
846}
847
848
849/** Calculates non local pseudopotential energy \f$E_{ps,nl,i}\f$ using Riemann tensor.
850 * \param *P Problem at hand
851 * \param i value
852 * \note not implemented
853 */
854void CalculateNonLocalEnergyUseRT(struct Problem *P, const int i)
855{
856 Error(SomeError, "CalculateNonLocalEnergyUseRT: Not implemented");
857}
858
859/* AddVICGP aus scp */ /* HG wird geloescht und neu gesetzt */
860
861/** Calculates Gauss, Pseudopotential, Hartreepotential and Hartree energies, also
862 * PseudoPot::VCoulombc \f$V^H (G)\f$ and Density::DensityCArray[DoubleDensityTypes#HGDensity].
863 * \f[
864 * E_{Gauss} = V \sum_{G\neq 0} \frac{4\pi}{|G|^2} |\overbrace{\sum_{I_s} S_{I_s} (G) \phi^{Gauss}_{I_s} (G)}^{\hat{n}^{Gauss}(G)}|^2 \qquad (\textnormal{just before 4.22})
865 * \f]
866 * \f[
867 * \widetilde{E}_{ps,loc} = 2V \sum_G \sum_{I_s} S_{I_s} (G) \phi^{ps,loc}_{I_s} (G) \hat{n}^\ast (G) \quad (4.13)
868 * \f]
869 * \f[
870 * E_{Hpot} = V \sum_{G\neq 0} \frac{4\pi}{|G|^2} |\hat{n}(G)+\hat{n}^{Gauss}(G)|^2 \qquad (\textnormal{first part 4.10})
871 * \f]
872 * \f[
873 * E_{Hartree} = V \sum_{G\neq 0} \frac{4\pi}{|G|^2} |\hat{n}(G)|^2
874 * \f]
875 * \f[
876 * V^H (G) = \frac{4\pi}{|G|^2} (\hat{n}(G)+\hat{n}^{Gauss}(G)) \qquad (\textnormal{section 4.3.2})
877 * \f]
878 * \param *P Problem at hand
879 * \note There are some factor 2 discrepancies to formulas due to gamma point symmetry!
880 */
881void CalculateSomeEnergyAndHGNoRT(struct Problem *P)
882{
883 struct PseudoPot *PP = &P->PP;
884 struct Ions *I = &P->Ion;
885 struct Lattice *Lat = &P->Lat;
886 struct Energy *E = Lat->E;
887 fftw_complex vp,rp,rhog,rhoe;
888 double SumEHP,Fac,SumEH,SumEG,SumEPS;
889 struct RunStruct *R = &P->R;
890 struct LatticeLevel *Lev0 = R->Lev0;
891 struct Density *Dens = Lev0->Dens;
892 int g,s=0,it,Index,i;
893 fftw_complex *HG = Dens->DensityCArray[HGDensity];
894 SetArrayToDouble0((double *)HG,2*Dens->TotalSize);
895 SumEHP = 0.0;
896 SumEH = 0.0;
897 SumEG = 0.0;
898 SumEPS = 0.0;
899 // calculates energy of local pseudopotential
900 if (Lev0->GArray[0].GSq == 0.0) {
901 Index = Lev0->GArray[0].Index;
902 c_re(vp) = 0.0;
903 c_im(vp) = 0.0;
904 for (it = 0; it < I->Max_Types; it++) {
905 c_re(vp) += (c_re(I->I[it].SFactor[0])*PP->phi_ps_loc[it][0]);
906 c_im(vp) += (c_im(I->I[it].SFactor[0])*PP->phi_ps_loc[it][0]);
907 }
908 c_re(HG[Index]) = c_re(vp);
909 SumEPS += (c_re(Dens->DensityCArray[TotalDensity][Index])*c_re(vp) +
910 c_im(Dens->DensityCArray[TotalDensity][Index])*c_im(vp))*R->HGcFactor;
911 s++;
912 }
913 for (g=s; g < Lev0->MaxG; g++) {
914 Index = Lev0->GArray[g].Index;
915 c_re(vp) = 0.0;
916 c_im(vp) = 0.0;
917 c_re(rp) = 0.0;
918 c_im(rp) = 0.0;
919 Fac = 4.*PI/(Lev0->GArray[g].GSq);
920 for (it = 0; it < I->Max_Types; it++) {
921 c_re(vp) += (c_re(I->I[it].SFactor[g])*PP->phi_ps_loc[it][g]);
922 c_im(vp) += (c_im(I->I[it].SFactor[g])*PP->phi_ps_loc[it][g]);
923 c_re(rp) += (c_re(I->I[it].SFactor[g])*PP->FacGauss[it][g]);
924 c_im(rp) += (c_im(I->I[it].SFactor[g])*PP->FacGauss[it][g]);
925 } // rp = n^{Gauss)(G)
926 c_re(rhog) = c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_re(rp);
927 c_im(rhog) = c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_im(rp);
928 c_re(rhoe) = c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
929 c_im(rhoe) = c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
930 // rhog = n(G) + n^{Gauss}(G), rhoe = n(G)
931 c_re(PP->VCoulombc[g]) = Fac*c_re(rhog);
932 c_im(PP->VCoulombc[g]) = -Fac*c_im(rhog);
933 c_re(HG[Index]) = c_re(vp)+Fac*c_re(rhog);
934 c_im(HG[Index]) = c_im(vp)+Fac*c_im(rhog);
935 /*if (P->first) */
936 SumEG += Fac*(c_re(rp)*c_re(rp)+c_im(rp)*c_im(rp)); // Gauss energy
937 SumEHP += Fac*(c_re(rhog)*c_re(rhog)+c_im(rhog)*c_im(rhog)); // E_ES first part
938 SumEH += Fac*(c_re(rhoe)*c_re(rhoe)+c_im(rhoe)*c_im(rhoe));
939 SumEPS += 2.*(c_re(Dens->DensityCArray[TotalDensity][Index])*c_re(vp)+
940 c_im(Dens->DensityCArray[TotalDensity][Index])*c_im(vp))*R->HGcFactor;
941 }
942 //
943 for (i=0; i<Lev0->MaxDoubleG; i++) {
944 HG[Lev0->DoubleG[2*i+1]].re = HG[Lev0->DoubleG[2*i]].re;
945 HG[Lev0->DoubleG[2*i+1]].im = -HG[Lev0->DoubleG[2*i]].im;
946 }
947 /*if (P->first) */
948 E->AllLocalDensityEnergy[GaussEnergy] = SumEG*Lat->Volume;
949 E->AllLocalDensityEnergy[PseudoEnergy] = SumEPS*Lat->Volume;
950 E->AllLocalDensityEnergy[HartreePotentialEnergy] = SumEHP*Lat->Volume;
951 E->AllLocalDensityEnergy[HartreeEnergy] = SumEH*Lat->Volume;
952 /* ReduceAllEnergy danach noch aufrufen */
953}
954
955/** Calculates \f$f^{nl}_{i,I_s,I_a,l,m}\f$ for conjugate direction Psis.
956 * \param *P Problem at hand
957 * \param *ConDir array of complex coefficients of the conjugate direction
958 * \param ***CDfnl return array [ion type, lm value, ion of type]
959 * \sa CalculateConDirHConDir() - used there
960 */
961void CalculateCDfnl(struct Problem *P, fftw_complex *ConDir, fftw_complex ***CDfnl)
962{
963 struct PseudoPot *PP = &P->PP;
964 struct Ions *I = &P->Ion;
965 struct RunStruct *R = &P->R;
966 const struct LatticeLevel *LevS = R->LevS;
967 int it,iot,ilm,ig,MinusGFactor,s;
968 fftw_complex PPdexpiGRpkg;
969 fftw_complex *dfnl = PP->dfnl;
970 for (it=0; it < I->Max_Types; it++) {
971 for (ilm=1; ilm <=PP->lm_end[it]; ilm++) {
972 MinusGFactor = Get_pkga_MinusG(ilm-1);
973 SetArrayToDouble0((double *)dfnl, 2*I->I[it].Max_IonsOfType);
974 for (iot=0; iot < I->I[it].Max_IonsOfType; iot++) {
975 s=0;
976 if (LevS->GArray[0].GSq == 0.0) {
977 dexpiGRpkg(PP, it, iot, ilm, 0, &PPdexpiGRpkg);
978 c_re(dfnl[iot]) += (c_re(ConDir[0])*c_re(PPdexpiGRpkg)-c_im(ConDir[0])*c_im(PPdexpiGRpkg));
979 c_im(dfnl[iot]) += (c_re(ConDir[0])*c_im(PPdexpiGRpkg)+c_im(ConDir[0])*c_re(PPdexpiGRpkg));
980 s++;
981 }
982 for (ig=s; ig < LevS->MaxG; ig++) {
983 dexpiGRpkg(PP, it, iot, ilm, ig, &PPdexpiGRpkg);
984 c_re(dfnl[iot]) += (c_re(ConDir[ig])*c_re(PPdexpiGRpkg)-c_im(ConDir[ig])*c_im(PPdexpiGRpkg));
985 c_im(dfnl[iot]) += (c_re(ConDir[ig])*c_im(PPdexpiGRpkg)+c_im(ConDir[ig])*c_re(PPdexpiGRpkg));
986 /* Minus G */
987 c_re(dfnl[iot]) += MinusGFactor*(c_re(ConDir[ig])*c_re(PPdexpiGRpkg)-c_im(ConDir[ig])*c_im(PPdexpiGRpkg));
988 c_im(dfnl[iot]) -= MinusGFactor*(c_re(ConDir[ig])*c_im(PPdexpiGRpkg)+c_im(ConDir[ig])*c_re(PPdexpiGRpkg));
989 }
990 }
991 MPI_Allreduce( dfnl, CDfnl[it][ilm-1], 2*I->I[it].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
992 }
993 }
994}
995
996/** Return non-local potential \f$V^{ps,nl}(G)|\Psi_{i} \rangle\f$.
997 * \a *HNL is set to zero and the non-local potential added.
998 * \param *P Problem at hand
999 * \param *HNL return array of real coefficients
1000 * \param ***fnl non-local form factors for specific wave function, see PseudoPot::fnl
1001 * \param PsiFactor occupation number of respective wave function
1002 * \sa CalculcateGradientNoRT() - used there in gradient calculation
1003 * CalculateConDirHConDir() - used there with conjugate directions (different non-local form factors!)
1004 */
1005void CalculateAddNLPot(struct Problem *P, fftw_complex *HNL, fftw_complex ***fnl, double PsiFactor)
1006{
1007 struct PseudoPot *PP = &P->PP;
1008 struct Ions *I = &P->Ion;
1009 struct RunStruct *R = &P->R;
1010 const struct LatticeLevel *LevS = R->LevS;
1011 int it,iot,ig,ilm;
1012 fftw_complex dexpiGR, dpkg, tt1, tt2, tt3, tt4, tt;
1013 double rr1=0, rr2=0, rr3=0, rr4=0;
1014 double *rr = PP->rr;
1015 fftw_complex *t = PP->t;
1016 SetArrayToDouble0((double *)HNL,2*LevS->MaxG);
1017 for (it=0; it < I->Max_Types; it++) { //over all types
1018 if (PP->lm_end[it] == 4) {
1019 rr1 = PP->wNonLoc[it][0];
1020 rr2 = PP->wNonLoc[it][1];
1021 rr3 = PP->wNonLoc[it][2];
1022 rr4 = PP->wNonLoc[it][3];
1023 } else {
1024 for (ilm=0; ilm < PP->lm_end[it]; ilm++) {
1025 rr[ilm] = PP->wNonLoc[it][ilm];
1026 }
1027 }
1028 for (iot =0; iot < I->I[it].Max_IonsOfType; iot++) { // over all ions of this type
1029 if (PP->lm_end[it] == 4) { // over all lm-values
1030 c_re(tt1) = -c_re(fnl[it][0][iot])*rr1;
1031 c_im(tt1) = -c_im(fnl[it][0][iot])*rr1;
1032 c_re(tt2) = -c_re(fnl[it][1][iot])*rr2;
1033 c_im(tt2) = -c_im(fnl[it][1][iot])*rr2;
1034 c_re(tt3) = -c_re(fnl[it][2][iot])*rr3;
1035 c_im(tt3) = -c_im(fnl[it][2][iot])*rr3;
1036 c_re(tt4) = -c_re(fnl[it][3][iot])*rr4;
1037 c_im(tt4) = -c_im(fnl[it][3][iot])*rr4;
1038 for (ig=0; ig < LevS->MaxG; ig++) {
1039 expiGR(PP, it, iot, ig, &dexpiGR);
1040 c_re(dpkg) = PP->phi_ps_nl[it][0][ig]*c_re(tt1)+PP->phi_ps_nl[it][1][ig]*c_re(tt2)+PP->phi_ps_nl[it][2][ig]*c_re(tt3)+PP->phi_ps_nl[it][3][ig]*c_re(tt4);
1041 c_im(dpkg) = PP->phi_ps_nl[it][0][ig]*c_im(tt1)+PP->phi_ps_nl[it][1][ig]*c_im(tt2)+PP->phi_ps_nl[it][2][ig]*c_im(tt3)+PP->phi_ps_nl[it][3][ig]*c_im(tt4);
1042 c_re(HNL[ig]) -= (c_re(dexpiGR)*c_re(dpkg)-c_im(dexpiGR)*c_im(dpkg))*PsiFactor;
1043 c_im(HNL[ig]) -= (c_re(dexpiGR)*c_im(dpkg)+c_im(dexpiGR)*c_re(dpkg))*PsiFactor;
1044 }
1045 } else {
1046 for (ilm=0; ilm < PP->lm_end[it]; ilm++) { // over all lm-values
1047 c_re(t[ilm]) = -c_re(fnl[it][ilm][iot])*rr[ilm];
1048 c_im(t[ilm]) = -c_im(fnl[it][ilm][iot])*rr[ilm];
1049 }
1050 if (PP->lm_end[it] > 0) {
1051 for (ig=0; ig < LevS->MaxG; ig++) {
1052 c_re(tt) = c_re(t[0])*PP->phi_ps_nl[it][0][ig];
1053 c_im(tt) = c_im(t[0])*PP->phi_ps_nl[it][0][ig];
1054 for (ilm=1; ilm < PP->lm_end[it]; ilm++) {
1055 c_re(tt) += c_re(t[ilm])*PP->phi_ps_nl[it][ilm][ig];
1056 c_im(tt) += c_im(t[ilm])*PP->phi_ps_nl[it][ilm][ig];
1057 }
1058 expiGR(PP,it,iot,ig,&dexpiGR);
1059 c_re(HNL[ig]) -= (c_re(dexpiGR)*c_re(tt)-c_im(dexpiGR)*c_im(tt))*PsiFactor;
1060 c_im(HNL[ig]) -= (c_re(dexpiGR)*c_im(tt)+c_im(dexpiGR)*c_re(tt))*PsiFactor;
1061 }
1062 }
1063 }
1064 }
1065 }
1066 if (LevS->GArray[0].GSq == 0.0)
1067 HNL[0].im = 0.0;
1068}
1069
1070/** Calculates the local force acting on ion.
1071 * \param *P Problem at hand
1072 */
1073void CalculateIonLocalForce(struct Problem *P)
1074{
1075 int is,ia,g2,Index,s,i;
1076 double *G;
1077 struct Lattice *L = &P->Lat;
1078 struct Ions *I = &P->Ion;
1079 struct PseudoPot *PP = &P->PP;
1080 struct RunStruct *R = &P->R;
1081 struct LatticeLevel *Lev0 = R->Lev0;
1082 struct Density *Dens = Lev0->Dens;
1083 double ForceFac = L->Volume;
1084 double force, *dsum;
1085 fftw_complex cv;
1086 for (is=0; is < I->Max_Types; is++) {
1087 SetArrayToDouble0(I->FTemp, NDIM*I->I[is].Max_IonsOfType);
1088 for (ia=0;ia< I->I[is].Max_IonsOfType; ia++) {
1089 dsum = &I->FTemp[ia*NDIM];
1090 s = 0;
1091 if (Lev0->GArray[0].GSq == 0.0)
1092 s++;
1093 for (g2=s; g2 < Lev0->MaxG; g2++) {
1094 Index = Lev0->GArray[g2].Index;
1095 G = Lev0->GArray[g2].G;
1096 c_re(cv) = c_re(PP->VCoulombc[g2])*PP->FacGauss[is][g2]+PP->phi_ps_loc[is][g2]*c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
1097 c_im(cv) = c_im(PP->VCoulombc[g2])*PP->FacGauss[is][g2]-PP->phi_ps_loc[is][g2]*c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
1098 force = c_re(PP->ei[is][ia][g2])*c_im(cv)+c_im(PP->ei[is][ia][g2])*c_re(cv);
1099 dsum[0] += 2.*(-G[0]*force);
1100 dsum[1] += 2.*(-G[1]*force);
1101 dsum[2] += 2.*(-G[2]*force);
1102 }
1103 }
1104 MPI_Allreduce( I->FTemp, I->I[is].FIonL, NDIM*I->I[is].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
1105 for (ia=0;ia< I->I[is].Max_IonsOfType; ia++)
1106 for (i=0; i<NDIM; i++)
1107 I->I[is].FIonL[i+NDIM*ia] *= ForceFac;
1108 }
1109}
1110
1111/** Calculates the non-local force acting on ion.
1112 * \param *P Problem at hand
1113 */
1114void CalculateIonNonLocalForce(struct Problem *P)
1115{
1116 double *G;
1117 double fnl[NDIM], AllLocalfnl[NDIM], AllUpfnl[NDIM], AllDownfnl[NDIM], AllTotalfnl[NDIM];
1118 double Fac;
1119 struct Lattice *L = &P->Lat;
1120 struct Psis *Psi = &L->Psi;
1121 struct Ions *I = &P->Ion;
1122 struct PseudoPot *PP = &P->PP;
1123 struct RunStruct *R = &P->R;
1124 struct LatticeLevel *LevS = R->LevS;
1125 int MinusGFactor;
1126 fftw_complex PPdexpiGRpkg, sf_a, sf_s, *LocalPsi;
1127 int is,ia,ilm,ig,i,s,d;
1128 MPI_Status status;
1129 for (is=0; is < I->Max_Types; is++) {
1130 SetArrayToDouble0(I->I[is].FIonNL, NDIM*I->I[is].Max_IonsOfType);
1131 for (ia=0;ia< I->I[is].Max_IonsOfType; ia++) {
1132 for (d=0; d<NDIM; d++)
1133 fnl[d] = 0.0;
1134 for (i=0; i < L->Psi.LocalNo; i++) if (Psi->LocalPsiStatus[i].PsiType == P->R.CurrentMin) {
1135 LocalPsi = LevS->LPsi->LocalPsi[i];
1136 for (ilm=1; ilm <=PP->lm_end[is]; ilm++) {
1137 MinusGFactor = Get_pkga_MinusG(ilm-1);
1138 Fac = Psi->LocalPsiStatus[i].PsiFactor*PP->wNonLoc[is][ilm-1];
1139 c_re(sf_s) = c_re(PP->fnl[i][is][ilm-1][ia]);
1140 c_im(sf_s) = c_im(PP->fnl[i][is][ilm-1][ia]);
1141 s=0;
1142 if (LevS->GArray[0].GSq == 0.0) {
1143 s++;
1144 }
1145 for (ig=s; ig < LevS->MaxG; ig++) {
1146 dexpiGRpkg(PP, is, ia, ilm, ig, &PPdexpiGRpkg);
1147 G = LevS->GArray[ig].G;
1148 for (d=0; d <NDIM; d++) {
1149 c_re(sf_a) = (c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]);
1150 c_im(sf_a) = (c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]);
1151 /* Minus G */
1152 c_re(sf_a) -= MinusGFactor*(c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]);
1153 c_im(sf_a) += MinusGFactor*(c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]);
1154 fnl[d] += Fac*(c_re(sf_a)*c_im(sf_s)-c_im(sf_a)*c_re(sf_s));
1155 }
1156 }
1157 }
1158 }
1159 MPI_Allreduce( fnl, AllLocalfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
1160 switch (Psi->PsiST) {
1161 case SpinDouble:
1162 MPI_Allreduce(AllLocalfnl, AllTotalfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
1163 break;
1164 case SpinUp:
1165 MPI_Allreduce(AllLocalfnl, AllUpfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
1166 MPI_Sendrecv(AllUpfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag3,
1167 AllDownfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag4,
1168 P->Par.comm_STInter, &status );
1169 for (d=0; d< NDIM; d++)
1170 AllTotalfnl[d] = AllUpfnl[d] + AllDownfnl[d];
1171 break;
1172 case SpinDown:
1173 MPI_Allreduce(AllLocalfnl, AllDownfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
1174 MPI_Sendrecv(AllDownfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag4,
1175 AllUpfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag3,
1176 P->Par.comm_STInter, &status );
1177 for (d=0; d < NDIM; d++)
1178 AllTotalfnl[d] = AllUpfnl[d] + AllDownfnl[d];
1179 break;
1180 }
1181 for (d=0; d<NDIM;d++)
1182 I->I[is].FIonNL[d+NDIM*ia] -= AllTotalfnl[d]*2.0;
1183 }
1184 }
1185}
Note: See TracBrowser for help on using the repository browser.