source: pcp/src/pseudo.c@ 255d13

Last change on this file since 255d13 was c510a7, checked in by Frederik Heber <heber@…>, 18 years ago

char lengths of 255 and MAXDUMMYSTRING replaced with define MAXSTRINGSIZE in molecuilder and pcp
  • Property mode set to 100644
File size: 52.6 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{
[d3482a]447 char *cpiInputFileName;
[a0bcf1]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;
[c510a7]485 cpiInputFileName = (char *) Malloc(sizeof(char)*MAXSTRINGSIZE, "InitPseudoRead: *cpiInputFileName");
[a0bcf1]486 for (it = 0; it < I->Max_Types; it++) {
487 PP->lm_end[it] = I->I[it].l_max*I->I[it].l_max+1-2*I->I[it].l_loc;
488 if (PP->lm_endmax < PP->lm_end[it]) PP->lm_endmax = PP->lm_end[it];
489 //PP->phi_ps_nl[it] = (double **) Malloc(sizeof(double*)*PP->lm_end[it], "InitPseudoRead: ");
490 PP->phi_ps_nl[it] = (double **) Malloc(sizeof(double*)*PP->lm_end[it], "InitPseudoRead: ");
491 PP->ei[it] = (fftw_complex **) Malloc(sizeof(fftw_complex*)*I->I[it].Max_IonsOfType, "InitPseudoRead: ");
492 PP->expiGR[it] = (fftw_complex **) Malloc(sizeof(fftw_complex*)*I->I[it].Max_IonsOfType, "InitPseudoRead: ");
493 for (ia=0;ia<I->I[it].Max_IonsOfType;ia++) {
494 PP->ei[it][ia] = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILev0->MaxG, "InitPseudoRead: ");
495 PP->expiGR[it][ia] = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILevS->MaxG, "InitPseudoRead: ");
496 }
497 for (il=0;il<PP->lm_end[it];il++) {
498 //PP->phi_ps_nl[it][il] = (double *) Malloc(sizeof(double)*ILevS->MaxG, "InitPseudoRead: ");
499 PP->phi_ps_nl[it][il] = (double *) Malloc(sizeof(double)*ILevS->MaxG, "InitPseudoRead: ");
500 }
501 PP->wNonLoc[it] = (double *) Malloc(sizeof(double)*PP->lm_end[it], "InitPseudoRead: ");
502 sprintf(cpiInputFileName,"%s/pseudo.%02i",pseudopot_path,I->I[it].Z);
503 cpiInputFile = fopen(cpiInputFileName,"r");
504 if (cpiInputFile == NULL)
505 Error(SomeError, "Cant't find pseudopot path or file");
506 else
507 if (P->Call.out[ReadOut]) fprintf(stderr,"%s was found\n",cpiInputFileName);
508 int zeile = 1;
[961b75]509 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->zval[it], 1, critical);
510 ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, int_type, &PP->nang[it], 1, critical);
[a0bcf1]511 I->TotalZval += PP->zval[it]*(I->I[it].Max_IonsOfType);
512 PP->core[it] = (double *) Malloc(sizeof(double)*2, "InitPseudoRead: ");
513 PP->rc[it] = (double *) Malloc(sizeof(double)*2, "InitPseudoRead: ");
[961b75]514 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->core[it][0], 1, critical);
515 ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->rc[it][0], 1, critical);
516 ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &PP->core[it][1], 1, critical);
517 ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &PP->rc[it][1], 1, critical);
[a0bcf1]518 PP->rcl[it] = (double **) Malloc(sizeof(double*)*3, "InitPseudoRead: ");
519 PP->al[it] = (double **) Malloc(sizeof(double*)*3, "InitPseudoRead: ");
520 PP->bl[it] = (double **) Malloc(sizeof(double*)*3, "InitPseudoRead: ");
521 PP->FacGauss[it] = (double *) Malloc(sizeof(double)*ILev0->MaxG, "InitPseudoRead: ");
522 PP->phi_ps_loc[it] = (double *) Malloc(sizeof(double)*ILev0->MaxG, "InitPseudoRead: ");
523 I->I[it].SFactor = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILev0->MaxG, "InitPseudoRead: ");
524 for (il=0; il < 3; il++) {
525 PP->rcl[it][il] = (double *) Malloc(sizeof(double)*3, "InitPseudoRead: ");
526 PP->al[it][il] = (double *) Malloc(sizeof(double)*3, "InitPseudoRead: ");
527 PP->bl[it][il] = (double *) Malloc(sizeof(double)*3, "InitPseudoRead: ");
528 for (ib = 0; ib < 3; ib++) {
[961b75]529 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->rcl[it][il][ib], 1, critical);
530 ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->al[it][il][ib], 1, critical);
531 ParseForParameter(0,cpiInputFile, NULL, 1, 3, zeile, double_type, &PP->bl[it][il][ib], 1, critical);
[a0bcf1]532 if (PP->rcl[it][il][ib] < MYEPSILON)
533 PP->rcl[it][il][ib] = 0.0;
534 else
535 PP->rcl[it][il][ib] = 1./sqrt(PP->rcl[it][il][ib]);
536 }
537 }
[961b75]538 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &PP->mmax[it], 1, critical);
539 ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, double_type, &PP->clog[it], 1, critical);
[a0bcf1]540 if (PP->mmax[it] > PP->Mmax) PP->Mmax = PP->mmax[it];
541 PP->clog[it] = log(PP->clog[it]);
542 PP->r[it] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
543 PP->integrand2[it] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
544 PP->R[it] = (double **) Malloc(sizeof(double*)*PP->nang[it], "InitPseudoRead: ");
545 PP->v_loc[it] = (double **) Malloc(sizeof(double*)*PP->nang[it], "InitPseudoRead: ");
546 for (il=0; il < PP->nang[it]; il++) {
547 PP->R[it][il] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
548 PP->v_loc[it][il] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
549 for (j=0;j< PP->mmax[it];j++) {
[961b75]550 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &m, 1, critical);
551 ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->r[it][j], 1, critical);
552 ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &PP->R[it][il][j], 1, critical);
553 ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &PP->v_loc[it][il][j], 1, critical);
[a0bcf1]554 }
555 if (il < PP->nang[it]-1) {
[961b75]556 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &PP->mmax[it], 1, critical);
557 ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, double_type, &PP->clog[it], 1, critical);
[a0bcf1]558 PP->clog[it] = log(PP->clog[it]);
559 }
560 }
561 I->I[it].corecorr = NotCoreCorrected;
562 count = 0;
[961b75]563 count += ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &dummyr, 1, optional);
564 count += ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &dummycorewave, 1, optional);
[a0bcf1]565 //fprintf(stderr,"(%i) %lg %lg\n",P->Par.me,dummyr,dummycorewave);
566 if (count == 2) {
567 PP->r[it][0] = dummyr;
568 PP->corewave[it] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
569 PP->formfCore[it] = (double *) Malloc(sizeof(double)*ILev0->MaxG, "InitPseudoRead: ");
570 I->I[it].corecorr = CoreCorrected;
571 PP->corecorr = CoreCorrected;
572 PP->corewave[it][0] = dummycorewave/(4.*PI);
[961b75]573 ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &dummyr, 1, critical);
574 ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &dummycorewave, 1, critical);
[a0bcf1]575 for (j=1;j < PP->mmax[it]; j++) {
[961b75]576 ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->r[it][j], 1, critical);
577 ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->corewave[it][j], 1, critical);
578 ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &dummyr, 1, critical);
579 ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &dummycorewave, 1, critical);
[a0bcf1]580 PP->corewave[it][j] /= (4.*PI);
581 }
582 } else {
583 PP->corewave[it] = NULL;
584 PP->formfCore[it] = NULL;
585 }
586 fclose(cpiInputFile);
587 }
[d3482a]588 Free(cpiInputFileName, "InitPseudoRead: *cpiInputFileName");
[a0bcf1]589 PP->fnl = (fftw_complex ****) Malloc(sizeof(fftw_complex***)*(Lat->Psi.LocalNo+1), "InitPseudoRead: ");
590 for (i=0; i< Lat->Psi.LocalNo+1; i++) {
591 PP->fnl[i] = (fftw_complex ***) Malloc(sizeof(fftw_complex**)*I->Max_Types, "InitPseudoRead: ");
592 for (it=0; it < I->Max_Types; it++) {
593 PP->fnl[i][it] = (fftw_complex **) Malloc(sizeof(fftw_complex*)*PP->lm_end[it], "InitPseudoRead: ");
594 for (il =0; il < PP->lm_end[it]; il++)
595 PP->fnl[i][it][il] = (fftw_complex *) Malloc(sizeof(fftw_complex)*I->I[it].Max_IonsOfType, "InitPseudoRead: ");
596 }
597 }
598 if (fabs(I->TotalZval - P->Lat.Psi.NIDensity[Occupied]) >= MYEPSILON) {
599 if (P->Par.me_comm_ST == 0) // instead of P->Par.me, thus differing charge density output also for SpinUp/Down from both processes!
600 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);
601 Error(SomeError,"Readparameters: charged system is not implemented yet");
602 }
603 PP->CDfnl = PP->fnl[Lat->Psi.LocalNo];
604 PP->dfnl = (fftw_complex *) Malloc(sizeof(fftw_complex)*I->Max_Max_IonsOfType, "InitPseudoRead: ");
605 PP->rr = (double *) Malloc(sizeof(double)*PP->lm_endmax, "InitPseudoRead: ");
606 PP->t = (fftw_complex *) Malloc(sizeof(fftw_complex)*PP->lm_endmax, "InitPseudoRead: ");
607 PP->integrand = (double *) Malloc(sizeof(double)*PP->Mmax, "InitPseudoRead: ");
608 PP->integrand1 = (double *) Malloc(sizeof(double)*PP->Mmax, "InitPseudoRead: ");
609 for (it=0;it < I->Max_Types; it++) {
610 if (I->I[it].corecorr == CoreCorrected) {
611 for (j=0; j < PP->mmax[it]; j++) {
612 PP->integrand[j] = 4.*PI*PP->corewave[it][j]*PP->r[it][j]*PP->r[it][j]*PP->r[it][j];
613 }
614 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]));
615 }
616 }
617 PP->fac1sqrt2PI = 1./sqrt(2.*PI);
618 PP->fac1sqrtPI = 1./sqrt(PI);
619 SpeedMeasure(P, InitSimTime, StartTimeDo);
620 SpeedMeasure(P, InitLocTime, StartTimeDo);
621 CalcSG(P);
622 CalcExpiGR(P);
623 FormFacGauss(P);
624 FormFacLocPot(P);
625 SpeedMeasure(P, InitLocTime, StopTimeDo);
626 SpeedMeasure(P, InitNonLocTime, StartTimeDo);
627 FormFacNonLocPot(P);
628 SpeedMeasure(P, InitNonLocTime, StopTimeDo);
629 SpeedMeasure(P, InitLocTime, StartTimeDo);
630 CalcCoreCorrection(P, 1);
631 SpeedMeasure(P, InitLocTime, StopTimeDo);
632 SpeedMeasure(P, InitSimTime, StopTimeDo);
633}
634
635/** Updates Pseudopotentials due to change of mesh width.
636 * Calls CalcSG(), CalcExpiGR(), recalculates the form factors by calling
637 * FormFacGauss(), FormFacLocPot(), FormFacLonLocPot() and finally
638 * CalcCoreCorrection(). All speed-measured in InitLocTime respectively
639 * InitNonLocTime.
640 * \param *P Problem at hand
641 */
642inline void ChangePseudoToLevUp(struct Problem *P)
643{
644 SpeedMeasure(P, InitLocTime, StartTimeDo);
645 CalcSG(P);
646 CalcExpiGR(P);
647 FormFacGauss(P);
648 FormFacLocPot(P);
649 SpeedMeasure(P, InitLocTime, StopTimeDo);
650 SpeedMeasure(P, InitNonLocTime, StartTimeDo);
651 FormFacNonLocPot(P);
652 SpeedMeasure(P, InitNonLocTime, StopTimeDo);
653 SpeedMeasure(P, InitLocTime, StartTimeDo);
654 CalcCoreCorrection(P, 1);
655 SpeedMeasure(P, InitLocTime, StopTimeDo);
656}
657
658/** Updates Pseudopotentials due to the newly minimized wave functions.
659 * Calls CalcSG(), CalcExpiGR() and CalcCoreCorrection().
660 * \param *P Problem at hand
661 */
662inline void UpdatePseudoToNewWaves(struct Problem *P)
663{
664 SpeedMeasure(P, InitLocTime, StartTimeDo);
665 CalcSG(P);
666 CalcExpiGR(P);
667 CalcCoreCorrection(P, 0);
668 SpeedMeasure(P, InitLocTime, StopTimeDo);
669}
670
671/** Frees memory allocated from Readin.
672 * All memory is free'd that was allocated in InitPseudoRead()
673 * \param *P Problem at hand
674 * \sa InitPseudoRead()
675 */
676void RemovePseudoRead(struct Problem *P)
677{
678 int i,it,il,ia;
679 struct PseudoPot *PP = &P->PP;
680 struct Ions *I = &P->Ion;
681 struct Lattice *Lat = &P->Lat;
[64fa9e]682 Free(PP->integrand1, "RemovePseudoRead: PP->integrand1");
683 Free(PP->integrand, "RemovePseudoRead: PP->integrand");
684 Free(PP->dfnl, "RemovePseudoRead: PP->dfnl");
685 Free(PP->rr, "RemovePseudoRead: PP->rr");
686 Free(PP->t, "RemovePseudoRead: PP->t");
[a0bcf1]687 for (i=0; i< Lat->Psi.LocalNo+1; i++) {
688 for (it=0; it < I->Max_Types; it++) {
689 for (il =0; il < PP->lm_end[it]; il++)
[64fa9e]690 Free(PP->fnl[i][it][il], "RemovePseudoRead: PP->fnl[i][it][il]");
691 Free(PP->fnl[i][it], "RemovePseudoRead: PP->fnl[i][it]");
[a0bcf1]692 }
[64fa9e]693 Free(PP->fnl[i], "RemovePseudoRead: PP->fnl[i]");
[a0bcf1]694 }
695 for (it=0; it < I->Max_Types; it++) {
696 if (I->I[it].corecorr == CoreCorrected) {
[64fa9e]697 Free(PP->corewave[it], "RemovePseudoRead: PP->corewave[it]");
698 Free(PP->formfCore[it], "RemovePseudoRead: PP->formfCore[it]");
[a0bcf1]699 }
700 }
701 for (it=0; it < I->Max_Types; it++) {
702 for (il=0; il < PP->nang[it]; il++) {
[64fa9e]703 Free(PP->R[it][il], "RemovePseudoRead: PP->R[it][il]");
704 Free(PP->v_loc[it][il], "RemovePseudoRead: PP->v_loc[it][il]");
[a0bcf1]705 }
706 for (il=0; il < 3; il++) {
[64fa9e]707 Free(PP->rcl[it][il], "RemovePseudoRead: PP->rcl[it][il]");
708 Free(PP->al[it][il], "RemovePseudoRead: PP->al[it][il]");
709 Free(PP->bl[it][il], "RemovePseudoRead: PP->bl[it][il]");
[a0bcf1]710 }
711 for (il=0;il<PP->lm_end[it];il++) {
[64fa9e]712 //Free(PP->phi_ps_nl[it][il], "RemovePseudoRead: PP->phi_ps_nl[it][il]");
713 Free(PP->phi_ps_nl[it][il], "RemovePseudoRead: PP->phi_ps_nl[it][il]");
[a0bcf1]714 }
715 for (ia=0;ia<I->I[it].Max_IonsOfType;ia++) {
[64fa9e]716 Free(PP->ei[it][ia], "RemovePseudoRead: PP->ei[it][ia]");
717 Free(PP->expiGR[it][ia], "RemovePseudoRead: PP->expiGR[it][ia]");
[a0bcf1]718 }
[64fa9e]719 Free(PP->ei[it], "RemovePseudoRead: PP->ei[it]");
720 Free(PP->expiGR[it], "RemovePseudoRead: PP->expiGR[it]");
721 //Free(PP->phi_ps_nl[it], "RemovePseudoRead: PP->phi_ps_nl[it]");
722 Free(PP->phi_ps_nl[it], "RemovePseudoRead: PP->phi_ps_nl[it]");
723 Free(PP->R[it], "RemovePseudoRead: PP->R[it]");
724 Free(PP->v_loc[it], "RemovePseudoRead: PP->v_loc[it]");
725 Free(PP->r[it], "RemovePseudoRead: PP->r[it]");
726 Free(PP->integrand2[it], "RemovePseudoRead: PP->integrand2[it]");
727 Free(PP->rcl[it], "RemovePseudoRead: PP->rcl[it]");
728 Free(PP->al[it], "RemovePseudoRead: PP->al[it]");
729 Free(PP->bl[it], "RemovePseudoRead: PP->bl[it]");
730 Free(PP->FacGauss[it], "RemovePseudoRead: PP->FacGauss[it]");
731 Free(PP->phi_ps_loc[it], "RemovePseudoRead: PP->phi_ps_loc[it]");
732 Free(PP->core[it], "RemovePseudoRead: PP->core[it]");
733 Free(PP->rc[it], "RemovePseudoRead: PP->rc[it]");
734 Free(PP->wNonLoc[it], "RemovePseudoRead: PP->wNonLoc[it]");
[a0bcf1]735 }
[64fa9e]736 //Free(PP->phi_ps_nl, "RemovePseudoRead: PP->phi_ps_nl");
737 Free(PP->phi_ps_nl, "RemovePseudoRead: PP->phi_ps_nl");
738 Free(PP->R, "RemovePseudoRead: PP->R");
739 Free(PP->v_loc, "RemovePseudoRead: PP->v_loc");
740 Free(PP->r, "RemovePseudoRead: PP->r");
741 Free(PP->integrand2, "RemovePseudoRead: PP->integrand2");
742 Free(PP->rcl, "RemovePseudoRead: PP->rcl");
743 Free(PP->al, "RemovePseudoRead: PP->al");
744 Free(PP->bl, "RemovePseudoRead: PP->bl");
745 Free(PP->FacGauss, "RemovePseudoRead: PP->FacGauss");
746 Free(PP->phi_ps_loc, "RemovePseudoRead: PP->phi_ps_loc");
747 Free(PP->core, "RemovePseudoRead: PP->core");
748 Free(PP->rc, "RemovePseudoRead: PP->rc");
749 Free(PP->wNonLoc, "RemovePseudoRead: PP->wNonLoc");
750 Free(PP->ei, "RemovePseudoRead: PP->ei");
751 Free(PP->expiGR, "RemovePseudoRead: PP->expiGR");
752 Free(PP->corewave, "RemovePseudoRead: PP->corewave");
753 Free(PP->formfCore, "RemovePseudoRead: PP->formfCore");
754 Free(PP->fnl, "RemovePseudoRead: PP->fnl");
755 Free(PP->VCoulombc, "RemovePseudoRead: PP->VCoulombc");
756 Free(PP->lm_end, "RemovePseudoRead: PP->lm_end");
757 Free(PP->zval, "RemovePseudoRead: PP->zval");
758 Free(PP->nang, "RemovePseudoRead: PP->nang");
759 Free(PP->mmax, "RemovePseudoRead: PP->mmax");
760 Free(PP->clog, "RemovePseudoRead: PP->clog");
[a0bcf1]761}
762
763/* Calculate Energy and Forces */
764
765/** Calculates Gauss (self) energy term of ions \f$E_{self}\f$.
766 * \f[
767 * 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)
768 * \f]
769 * \param *P Problem at hand
770 * \param *PP PseudoPot ential structure
771 * \param *I Ions structure
772 */
773void CalculateGaussSelfEnergyNoRT(struct Problem *P, struct PseudoPot *PP, struct Ions *I)
774{
775 double SumE = 0.0;
776 int it;
777 for (it=0; it < I->Max_Types; it++) {
778 if (I->I[it].rgauss < MYEPSILON) fprintf(stderr,"CalculateGaussSelfEnergyNoRT: I->I[it].rgauss = %lg\n",I->I[it].rgauss);
779 SumE += PP->zval[it]*PP->zval[it]/I->I[it].rgauss*I->I[it].Max_IonsOfType;
780 }
781 P->Lat.E->AllTotalIonsEnergy[GaussSelfEnergy] = SumE*PP->fac1sqrt2PI;
782}
783
784/** Calculates non local pseudopotential energy \f$E_{ps,nl,i}\f$.
785 * First, calculates non-local form factors
786 * \f[
787 * f^{nl}_{i,I_s,I_a,l,m} = \sum_G \exp(-i(G R_{I_s,I_a})
788 * \phi_{I_s,l,m}^{ps,nl} (G) c_{i,G}
789 * \f]
790 * and afterwards evaluates with these
791 * \f[
792 * 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)
793 * \f]
794 * \param *P Problem at hand
795 * \param i Energy of i-th Psi
796 */
797void CalculateNonLocalEnergyNoRT(struct Problem *P, const int i)
798{
799 struct PseudoPot *PP = &P->PP;
800 struct Ions *I = &P->Ion;
801 struct Lattice *Lat = &P->Lat;
802 struct RunStruct *R = &P->R;
803 struct LatticeLevel *LevS = R->LevS;
804 int it,iot,ilm,ig,s=0;
805 double sf,enl;//,enlm;
806 fftw_complex PPdexpiGRpkg;
807 fftw_complex *dfnl = PP->dfnl;
808 fftw_complex *LocalPsi = LevS->LPsi->LocalPsi[i]; // i-th wave function coefficients
809 const double PsiFactor = Lat->Psi.LocalPsiStatus[i].PsiFactor;
810 int MinusGFactor;
811 enl=0.0;
812 for (it=0; it < I->Max_Types; it++) { // through all ion types
813 for (ilm=1; ilm <=PP->lm_end[it]; ilm++) { // over all possible m values
814 MinusGFactor = Get_pkga_MinusG(ilm-1);
815 SetArrayToDouble0((double *)dfnl, 2*I->I[it].Max_IonsOfType);
816 for (iot=0; iot < I->I[it].Max_IonsOfType; iot++) { // for each ion per type
817 s=0;
818 if (LevS->GArray[0].GSq == 0.0) {
819 dexpiGRpkg(PP, it, iot, ilm, 0, &PPdexpiGRpkg);
820 c_re(dfnl[iot]) += (c_re(LocalPsi[0])*c_re(PPdexpiGRpkg)-c_im(LocalPsi[0])*c_im(PPdexpiGRpkg));
821 c_im(dfnl[iot]) += (c_re(LocalPsi[0])*c_im(PPdexpiGRpkg)+c_im(LocalPsi[0])*c_re(PPdexpiGRpkg));
822 s++;
823 }
824 for (ig=s; ig < LevS->MaxG; ig++) {
825 dexpiGRpkg(PP, it, iot, ilm, ig, &PPdexpiGRpkg);
826 c_re(dfnl[iot]) += (c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg));
827 c_im(dfnl[iot]) += (c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg));
828 /* Minus G - due to inherent symmetry at gamma point! */
829 c_re(dfnl[iot]) += MinusGFactor*(c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg));
830 c_im(dfnl[iot]) -= MinusGFactor*(c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg));
831 }
832 }
833 // Allreduce as the coefficients are spread over many processes
834 MPI_Allreduce( dfnl, PP->fnl[i][it][ilm-1], 2*I->I[it].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
835 }
836 }
837 for (it=0; it < I->Max_Types; it++) {
838 for (ilm=1; ilm <=PP->lm_end[it]; ilm++) {
839 for (iot=0; iot < I->I[it].Max_IonsOfType; iot++) {
840 //enlm = 0.0;
841 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]);
842 //enlm = (PsiFactor*sf);
843 enl += PsiFactor*sf*PP->wNonLoc[it][ilm-1];
844 }
845 }
846 }
847 Lat->Energy[R->CurrentMin].PsiEnergy[NonLocalEnergy][i] = enl;
848}
849
850
851/** Calculates non local pseudopotential energy \f$E_{ps,nl,i}\f$ using Riemann tensor.
852 * \param *P Problem at hand
853 * \param i value
854 * \note not implemented
855 */
856void CalculateNonLocalEnergyUseRT(struct Problem *P, const int i)
857{
858 Error(SomeError, "CalculateNonLocalEnergyUseRT: Not implemented");
859}
860
861/* AddVICGP aus scp */ /* HG wird geloescht und neu gesetzt */
862
863/** Calculates Gauss, Pseudopotential, Hartreepotential and Hartree energies, also
864 * PseudoPot::VCoulombc \f$V^H (G)\f$ and Density::DensityCArray[DoubleDensityTypes#HGDensity].
865 * \f[
866 * 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})
867 * \f]
868 * \f[
869 * \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)
870 * \f]
871 * \f[
872 * 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})
873 * \f]
874 * \f[
875 * E_{Hartree} = V \sum_{G\neq 0} \frac{4\pi}{|G|^2} |\hat{n}(G)|^2
876 * \f]
877 * \f[
878 * V^H (G) = \frac{4\pi}{|G|^2} (\hat{n}(G)+\hat{n}^{Gauss}(G)) \qquad (\textnormal{section 4.3.2})
879 * \f]
880 * \param *P Problem at hand
881 * \note There are some factor 2 discrepancies to formulas due to gamma point symmetry!
882 */
883void CalculateSomeEnergyAndHGNoRT(struct Problem *P)
884{
885 struct PseudoPot *PP = &P->PP;
886 struct Ions *I = &P->Ion;
887 struct Lattice *Lat = &P->Lat;
888 struct Energy *E = Lat->E;
889 fftw_complex vp,rp,rhog,rhoe;
890 double SumEHP,Fac,SumEH,SumEG,SumEPS;
891 struct RunStruct *R = &P->R;
892 struct LatticeLevel *Lev0 = R->Lev0;
893 struct Density *Dens = Lev0->Dens;
894 int g,s=0,it,Index,i;
895 fftw_complex *HG = Dens->DensityCArray[HGDensity];
896 SetArrayToDouble0((double *)HG,2*Dens->TotalSize);
897 SumEHP = 0.0;
898 SumEH = 0.0;
899 SumEG = 0.0;
900 SumEPS = 0.0;
901 // calculates energy of local pseudopotential
902 if (Lev0->GArray[0].GSq == 0.0) {
903 Index = Lev0->GArray[0].Index;
904 c_re(vp) = 0.0;
905 c_im(vp) = 0.0;
906 for (it = 0; it < I->Max_Types; it++) {
907 c_re(vp) += (c_re(I->I[it].SFactor[0])*PP->phi_ps_loc[it][0]);
908 c_im(vp) += (c_im(I->I[it].SFactor[0])*PP->phi_ps_loc[it][0]);
909 }
910 c_re(HG[Index]) = c_re(vp);
911 SumEPS += (c_re(Dens->DensityCArray[TotalDensity][Index])*c_re(vp) +
912 c_im(Dens->DensityCArray[TotalDensity][Index])*c_im(vp))*R->HGcFactor;
913 s++;
914 }
915 for (g=s; g < Lev0->MaxG; g++) {
916 Index = Lev0->GArray[g].Index;
917 c_re(vp) = 0.0;
918 c_im(vp) = 0.0;
919 c_re(rp) = 0.0;
920 c_im(rp) = 0.0;
921 Fac = 4.*PI/(Lev0->GArray[g].GSq);
922 for (it = 0; it < I->Max_Types; it++) {
[ce81a3a]923 c_re(vp) += (c_re(I->I[it].SFactor[g])*PP->phi_ps_loc[it][g]); // V^ps,loc
[a0bcf1]924 c_im(vp) += (c_im(I->I[it].SFactor[g])*PP->phi_ps_loc[it][g]);
[ce81a3a]925 c_re(rp) += (c_re(I->I[it].SFactor[g])*PP->FacGauss[it][g]); // n^gauss
[a0bcf1]926 c_im(rp) += (c_im(I->I[it].SFactor[g])*PP->FacGauss[it][g]);
927 } // rp = n^{Gauss)(G)
[ce81a3a]928 c_re(rhog) = c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_re(rp); // n + n^gauss = n^~
[a0bcf1]929 c_im(rhog) = c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_im(rp);
930 c_re(rhoe) = c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
931 c_im(rhoe) = c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
932 // rhog = n(G) + n^{Gauss}(G), rhoe = n(G)
933 c_re(PP->VCoulombc[g]) = Fac*c_re(rhog);
934 c_im(PP->VCoulombc[g]) = -Fac*c_im(rhog);
935 c_re(HG[Index]) = c_re(vp)+Fac*c_re(rhog);
936 c_im(HG[Index]) = c_im(vp)+Fac*c_im(rhog);
937 /*if (P->first) */
938 SumEG += Fac*(c_re(rp)*c_re(rp)+c_im(rp)*c_im(rp)); // Gauss energy
939 SumEHP += Fac*(c_re(rhog)*c_re(rhog)+c_im(rhog)*c_im(rhog)); // E_ES first part
[ce81a3a]940 SumEH += Fac*(c_re(rhoe)*c_re(rhoe)+c_im(rhoe)*c_im(rhoe)); // Hartree energy
[a0bcf1]941 SumEPS += 2.*(c_re(Dens->DensityCArray[TotalDensity][Index])*c_re(vp)+
942 c_im(Dens->DensityCArray[TotalDensity][Index])*c_im(vp))*R->HGcFactor;
943 }
944 //
945 for (i=0; i<Lev0->MaxDoubleG; i++) {
946 HG[Lev0->DoubleG[2*i+1]].re = HG[Lev0->DoubleG[2*i]].re;
947 HG[Lev0->DoubleG[2*i+1]].im = -HG[Lev0->DoubleG[2*i]].im;
948 }
949 /*if (P->first) */
950 E->AllLocalDensityEnergy[GaussEnergy] = SumEG*Lat->Volume;
951 E->AllLocalDensityEnergy[PseudoEnergy] = SumEPS*Lat->Volume;
952 E->AllLocalDensityEnergy[HartreePotentialEnergy] = SumEHP*Lat->Volume;
953 E->AllLocalDensityEnergy[HartreeEnergy] = SumEH*Lat->Volume;
954 /* ReduceAllEnergy danach noch aufrufen */
955}
956
957/** Calculates \f$f^{nl}_{i,I_s,I_a,l,m}\f$ for conjugate direction Psis.
958 * \param *P Problem at hand
959 * \param *ConDir array of complex coefficients of the conjugate direction
960 * \param ***CDfnl return array [ion type, lm value, ion of type]
961 * \sa CalculateConDirHConDir() - used there
962 */
963void CalculateCDfnl(struct Problem *P, fftw_complex *ConDir, fftw_complex ***CDfnl)
964{
965 struct PseudoPot *PP = &P->PP;
966 struct Ions *I = &P->Ion;
967 struct RunStruct *R = &P->R;
968 const struct LatticeLevel *LevS = R->LevS;
969 int it,iot,ilm,ig,MinusGFactor,s;
970 fftw_complex PPdexpiGRpkg;
971 fftw_complex *dfnl = PP->dfnl;
972 for (it=0; it < I->Max_Types; it++) {
973 for (ilm=1; ilm <=PP->lm_end[it]; ilm++) {
974 MinusGFactor = Get_pkga_MinusG(ilm-1);
975 SetArrayToDouble0((double *)dfnl, 2*I->I[it].Max_IonsOfType);
976 for (iot=0; iot < I->I[it].Max_IonsOfType; iot++) {
977 s=0;
978 if (LevS->GArray[0].GSq == 0.0) {
979 dexpiGRpkg(PP, it, iot, ilm, 0, &PPdexpiGRpkg);
980 c_re(dfnl[iot]) += (c_re(ConDir[0])*c_re(PPdexpiGRpkg)-c_im(ConDir[0])*c_im(PPdexpiGRpkg));
981 c_im(dfnl[iot]) += (c_re(ConDir[0])*c_im(PPdexpiGRpkg)+c_im(ConDir[0])*c_re(PPdexpiGRpkg));
982 s++;
983 }
984 for (ig=s; ig < LevS->MaxG; ig++) {
985 dexpiGRpkg(PP, it, iot, ilm, ig, &PPdexpiGRpkg);
986 c_re(dfnl[iot]) += (c_re(ConDir[ig])*c_re(PPdexpiGRpkg)-c_im(ConDir[ig])*c_im(PPdexpiGRpkg));
987 c_im(dfnl[iot]) += (c_re(ConDir[ig])*c_im(PPdexpiGRpkg)+c_im(ConDir[ig])*c_re(PPdexpiGRpkg));
988 /* Minus G */
989 c_re(dfnl[iot]) += MinusGFactor*(c_re(ConDir[ig])*c_re(PPdexpiGRpkg)-c_im(ConDir[ig])*c_im(PPdexpiGRpkg));
990 c_im(dfnl[iot]) -= MinusGFactor*(c_re(ConDir[ig])*c_im(PPdexpiGRpkg)+c_im(ConDir[ig])*c_re(PPdexpiGRpkg));
991 }
992 }
993 MPI_Allreduce( dfnl, CDfnl[it][ilm-1], 2*I->I[it].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
994 }
995 }
996}
997
998/** Return non-local potential \f$V^{ps,nl}(G)|\Psi_{i} \rangle\f$.
999 * \a *HNL is set to zero and the non-local potential added.
1000 * \param *P Problem at hand
1001 * \param *HNL return array of real coefficients
1002 * \param ***fnl non-local form factors for specific wave function, see PseudoPot::fnl
1003 * \param PsiFactor occupation number of respective wave function
1004 * \sa CalculcateGradientNoRT() - used there in gradient calculation
1005 * CalculateConDirHConDir() - used there with conjugate directions (different non-local form factors!)
1006 */
1007void CalculateAddNLPot(struct Problem *P, fftw_complex *HNL, fftw_complex ***fnl, double PsiFactor)
1008{
1009 struct PseudoPot *PP = &P->PP;
1010 struct Ions *I = &P->Ion;
1011 struct RunStruct *R = &P->R;
1012 const struct LatticeLevel *LevS = R->LevS;
1013 int it,iot,ig,ilm;
1014 fftw_complex dexpiGR, dpkg, tt1, tt2, tt3, tt4, tt;
1015 double rr1=0, rr2=0, rr3=0, rr4=0;
1016 double *rr = PP->rr;
1017 fftw_complex *t = PP->t;
1018 SetArrayToDouble0((double *)HNL,2*LevS->MaxG);
1019 for (it=0; it < I->Max_Types; it++) { //over all types
1020 if (PP->lm_end[it] == 4) {
1021 rr1 = PP->wNonLoc[it][0];
1022 rr2 = PP->wNonLoc[it][1];
1023 rr3 = PP->wNonLoc[it][2];
1024 rr4 = PP->wNonLoc[it][3];
1025 } else {
1026 for (ilm=0; ilm < PP->lm_end[it]; ilm++) {
1027 rr[ilm] = PP->wNonLoc[it][ilm];
1028 }
1029 }
1030 for (iot =0; iot < I->I[it].Max_IonsOfType; iot++) { // over all ions of this type
1031 if (PP->lm_end[it] == 4) { // over all lm-values
1032 c_re(tt1) = -c_re(fnl[it][0][iot])*rr1;
1033 c_im(tt1) = -c_im(fnl[it][0][iot])*rr1;
1034 c_re(tt2) = -c_re(fnl[it][1][iot])*rr2;
1035 c_im(tt2) = -c_im(fnl[it][1][iot])*rr2;
1036 c_re(tt3) = -c_re(fnl[it][2][iot])*rr3;
1037 c_im(tt3) = -c_im(fnl[it][2][iot])*rr3;
1038 c_re(tt4) = -c_re(fnl[it][3][iot])*rr4;
1039 c_im(tt4) = -c_im(fnl[it][3][iot])*rr4;
1040 for (ig=0; ig < LevS->MaxG; ig++) {
1041 expiGR(PP, it, iot, ig, &dexpiGR);
1042 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);
1043 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);
1044 c_re(HNL[ig]) -= (c_re(dexpiGR)*c_re(dpkg)-c_im(dexpiGR)*c_im(dpkg))*PsiFactor;
1045 c_im(HNL[ig]) -= (c_re(dexpiGR)*c_im(dpkg)+c_im(dexpiGR)*c_re(dpkg))*PsiFactor;
1046 }
1047 } else {
1048 for (ilm=0; ilm < PP->lm_end[it]; ilm++) { // over all lm-values
1049 c_re(t[ilm]) = -c_re(fnl[it][ilm][iot])*rr[ilm];
1050 c_im(t[ilm]) = -c_im(fnl[it][ilm][iot])*rr[ilm];
1051 }
1052 if (PP->lm_end[it] > 0) {
1053 for (ig=0; ig < LevS->MaxG; ig++) {
1054 c_re(tt) = c_re(t[0])*PP->phi_ps_nl[it][0][ig];
1055 c_im(tt) = c_im(t[0])*PP->phi_ps_nl[it][0][ig];
1056 for (ilm=1; ilm < PP->lm_end[it]; ilm++) {
1057 c_re(tt) += c_re(t[ilm])*PP->phi_ps_nl[it][ilm][ig];
1058 c_im(tt) += c_im(t[ilm])*PP->phi_ps_nl[it][ilm][ig];
1059 }
1060 expiGR(PP,it,iot,ig,&dexpiGR);
1061 c_re(HNL[ig]) -= (c_re(dexpiGR)*c_re(tt)-c_im(dexpiGR)*c_im(tt))*PsiFactor;
1062 c_im(HNL[ig]) -= (c_re(dexpiGR)*c_im(tt)+c_im(dexpiGR)*c_re(tt))*PsiFactor;
1063 }
1064 }
1065 }
1066 }
1067 }
1068 if (LevS->GArray[0].GSq == 0.0)
1069 HNL[0].im = 0.0;
1070}
1071
[f72dc1]1072/** Calculates the local force acting on each ion.
[a0bcf1]1073 * \param *P Problem at hand
1074 */
1075void CalculateIonLocalForce(struct Problem *P)
1076{
[0f2c43]1077 int is,ia,g2,Index,i;
[a0bcf1]1078 double *G;
1079 struct Lattice *L = &P->Lat;
1080 struct Ions *I = &P->Ion;
1081 struct PseudoPot *PP = &P->PP;
1082 struct RunStruct *R = &P->R;
1083 struct LatticeLevel *Lev0 = R->Lev0;
[f72dc1]1084 struct LatticeLevel *LevS = R->LevS;
[a0bcf1]1085 struct Density *Dens = Lev0->Dens;
[f72dc1]1086 struct fft_plan_3d *plan = L->plan;
1087 fftw_complex *work = Dens->DensityCArray[TempDensity];
1088 fftw_complex *HGC = Dens->DensityCArray[HGDensity];
1089 fftw_real *HGCR = (fftw_real *)HGC;
[a0bcf1]1090 double ForceFac = L->Volume;
1091 double force, *dsum;
1092 fftw_complex cv;
[f72dc1]1093
1094 // precalc V_XC(r) and fft
1095 //debug(P,"precalc V_XC\n");
1096 SetArrayToDouble0((double *)HGCR, Dens->TotalSize*2);
1097 CalculateXCPotentialNoRT(P, HGCR);
1098 fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, HGC, work);
1099 //debug(P,"precalc V_XC finished\n");
[a0bcf1]1100 for (is=0; is < I->Max_Types; is++) {
1101 SetArrayToDouble0(I->FTemp, NDIM*I->I[is].Max_IonsOfType);
1102 for (ia=0;ia< I->I[is].Max_IonsOfType; ia++) {
1103 dsum = &I->FTemp[ia*NDIM];
[0f2c43]1104 g2 = 0;
1105 if (fabs(Lev0->GArray[g2].GSq) < MYEPSILON)
1106 g2++;
1107 for (; g2 < Lev0->MaxG; g2++) {
[a0bcf1]1108 Index = Lev0->GArray[g2].Index;
1109 G = Lev0->GArray[g2].G;
[0f2c43]1110 //debug(P,"Lokal Pseudopotential");
1111 c_re(cv) = PP->phi_ps_loc[is][g2]*c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
1112 c_im(cv) = -PP->phi_ps_loc[is][g2]*c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
1113 //debug(P,"Coulombpotential");
1114 c_re(cv) += c_re(PP->VCoulombc[g2])*PP->FacGauss[is][g2];
1115 c_im(cv) += c_im(PP->VCoulombc[g2])*PP->FacGauss[is][g2];
[f72dc1]1116 if (I->I[is].corecorr == CoreCorrected) {
1117 // this summand was missing before, is not thoroughly tested
1118 //debug(P,"V_XC summand in CalculateIonLocalForce()\n");
1119 c_re(cv) += (PP->formfCore[is][g2]*c_re(HGC[Index]))*R->HGcFactor;
1120 c_im(cv) += -(PP->formfCore[is][g2]*c_im(HGC[Index]))*R->HGcFactor;
1121 //fprintf(stderr, "(%i) FormfCore %lg\tHGC %lg+i%lg\n",P->Par.me, PP->formfCore[is][g2],c_re(HGC[Index]),c_im(HGC[Index]));
1122 }
[a0bcf1]1123 force = c_re(PP->ei[is][ia][g2])*c_im(cv)+c_im(PP->ei[is][ia][g2])*c_re(cv);
[f72dc1]1124 dsum[0] += 2.*(-G[0]*force); //2.* ... why? Nowhere from derivative!
1125 dsum[1] += 2.*(-G[1]*force); //2.*
1126 dsum[2] += 2.*(-G[2]*force); //2.*
[a0bcf1]1127 }
1128 }
1129 MPI_Allreduce( I->FTemp, I->I[is].FIonL, NDIM*I->I[is].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
1130 for (ia=0;ia< I->I[is].Max_IonsOfType; ia++)
1131 for (i=0; i<NDIM; i++)
[f72dc1]1132 I->I[is].FIonL[i+NDIM*ia] *= ForceFac;
[a0bcf1]1133 }
1134}
1135
1136/** Calculates the non-local force acting on ion.
1137 * \param *P Problem at hand
1138 */
1139void CalculateIonNonLocalForce(struct Problem *P)
1140{
1141 double *G;
1142 double fnl[NDIM], AllLocalfnl[NDIM], AllUpfnl[NDIM], AllDownfnl[NDIM], AllTotalfnl[NDIM];
1143 double Fac;
1144 struct Lattice *L = &P->Lat;
1145 struct Psis *Psi = &L->Psi;
1146 struct Ions *I = &P->Ion;
1147 struct PseudoPot *PP = &P->PP;
1148 struct RunStruct *R = &P->R;
1149 struct LatticeLevel *LevS = R->LevS;
1150 int MinusGFactor;
1151 fftw_complex PPdexpiGRpkg, sf_a, sf_s, *LocalPsi;
1152 int is,ia,ilm,ig,i,s,d;
1153 MPI_Status status;
1154 for (is=0; is < I->Max_Types; is++) {
1155 SetArrayToDouble0(I->I[is].FIonNL, NDIM*I->I[is].Max_IonsOfType);
1156 for (ia=0;ia< I->I[is].Max_IonsOfType; ia++) {
1157 for (d=0; d<NDIM; d++)
1158 fnl[d] = 0.0;
1159 for (i=0; i < L->Psi.LocalNo; i++) if (Psi->LocalPsiStatus[i].PsiType == P->R.CurrentMin) {
1160 LocalPsi = LevS->LPsi->LocalPsi[i];
1161 for (ilm=1; ilm <=PP->lm_end[is]; ilm++) {
1162 MinusGFactor = Get_pkga_MinusG(ilm-1);
1163 Fac = Psi->LocalPsiStatus[i].PsiFactor*PP->wNonLoc[is][ilm-1];
1164 c_re(sf_s) = c_re(PP->fnl[i][is][ilm-1][ia]);
1165 c_im(sf_s) = c_im(PP->fnl[i][is][ilm-1][ia]);
1166 s=0;
1167 if (LevS->GArray[0].GSq == 0.0) {
1168 s++;
1169 }
1170 for (ig=s; ig < LevS->MaxG; ig++) {
1171 dexpiGRpkg(PP, is, ia, ilm, ig, &PPdexpiGRpkg);
1172 G = LevS->GArray[ig].G;
1173 for (d=0; d <NDIM; d++) {
1174 c_re(sf_a) = (c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]);
1175 c_im(sf_a) = (c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]);
1176 /* Minus G */
1177 c_re(sf_a) -= MinusGFactor*(c_re(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]-c_im(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]);
1178 c_im(sf_a) += MinusGFactor*(c_re(LocalPsi[ig])*c_im(PPdexpiGRpkg)*G[d]+c_im(LocalPsi[ig])*c_re(PPdexpiGRpkg)*G[d]);
1179 fnl[d] += Fac*(c_re(sf_a)*c_im(sf_s)-c_im(sf_a)*c_re(sf_s));
1180 }
1181 }
1182 }
1183 }
1184 MPI_Allreduce( fnl, AllLocalfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
1185 switch (Psi->PsiST) {
1186 case SpinDouble:
1187 MPI_Allreduce(AllLocalfnl, AllTotalfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
1188 break;
1189 case SpinUp:
1190 MPI_Allreduce(AllLocalfnl, AllUpfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
1191 MPI_Sendrecv(AllUpfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag3,
1192 AllDownfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag4,
1193 P->Par.comm_STInter, &status );
1194 for (d=0; d< NDIM; d++)
1195 AllTotalfnl[d] = AllUpfnl[d] + AllDownfnl[d];
1196 break;
1197 case SpinDown:
1198 MPI_Allreduce(AllLocalfnl, AllDownfnl, NDIM, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
1199 MPI_Sendrecv(AllDownfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag4,
1200 AllUpfnl, NDIM, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag3,
1201 P->Par.comm_STInter, &status );
1202 for (d=0; d < NDIM; d++)
1203 AllTotalfnl[d] = AllUpfnl[d] + AllDownfnl[d];
1204 break;
1205 }
1206 for (d=0; d<NDIM;d++)
1207 I->I[is].FIonNL[d+NDIM*ia] -= AllTotalfnl[d]*2.0;
1208 }
1209 }
1210}
Note: See TracBrowser for help on using the repository browser.