source: pcp/src/pseudo.c@ f70c2a

Last change on this file since f70c2a was e08f45, checked in by Frederik Heber <heber@…>, 17 years ago

Merge branch 'ConcaveHull' of ../espack2 into ConcaveHull

Conflicts:

molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/linkedcell.cpp
molecuilder/src/linkedcell.hpp
molecuilder/src/vector.cpp
molecuilder/src/vector.hpp
util/src/NanoCreator.c

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