source: pcp/src/grad.c@ fdbf0c

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

config.h is included in each and every file. After trying to compile on JUMP (with xlc).
  • Property mode set to 100644
File size: 58.9 KB
Line 
1/** \file grad.c
2 * Gradient Calculation.
3 *
4 * The most important structure is Gradient with its Gradient::GradientArray.
5 * It contains the various stages of the calculated conjugate direction, beginning
6 * from the Psis::LocalPsiStatus. Subsequently, CalculateCGGradient(), CalculatePreConGrad(),
7 * CalculateConDir() and CalculateLineSearch() are called which use earlier stages to
8 * calculate the latter ones: ActualGradient, PreConDirGradient, ConDirGradient and OldConDirGradient,
9 * do the one-dimensional minimisation and finally update the coefficients of the wave function,
10 * whilst these ought to remain orthogonal.
11 *
12 * Initialization of the Gradient structure and arrays is done in InitGradients(), deallocation
13 * in RemoveGradients().
14 *
15 * Small functions calculate the scalar product between two gradient directions GradSP(), find
16 * the minimal angle in the line search CalculateDeltaI() or evaluate complex expressions
17 * CalculateConDirHConDir(), down to basic calculation of the gradient of the wave function
18 * CalculateGradientNoRT(). CalculateHamiltonian() sets up the Hamiltonian for the Kohn-Sham
19 * orbitals and calculates the eigenvalues.
20 *
21 * Other functions update the coefficients after the position of the ion has changed -
22 * CalcAlphaForUpdateWaveAfterIonMove() and UpdateWaveAfterIonMove().
23 *
24 *
25 Project: ParallelCarParrinello
26 \author Jan Hamaekers
27 \date 2000
28
29 File: grad.c
30 $Id: grad.c,v 1.90 2007-03-29 13:38:04 foo Exp $
31*/
32
33#ifdef HAVE_CONFIG_H
34#include <config.h>
35#endif
36
37#include <stdlib.h>
38#include <stdio.h>
39#include <stddef.h>
40#include <math.h>
41#include <string.h>
42#include <gsl/gsl_matrix.h>
43#include <gsl/gsl_eigen.h>
44#include <gsl/gsl_heapsort.h>
45#include <gsl/gsl_complex.h>
46#include <gsl/gsl_complex_math.h>
47#include <gsl/gsl_sort_vector.h>
48#include <gsl/gsl_blas.h>
49
50#include "mymath.h"
51#include "data.h"
52#include "errors.h"
53#include "energy.h"
54#include "helpers.h"
55#include "myfft.h"
56#include "density.h"
57#include "gramsch.h"
58#include "perturbed.h"
59#include "pseudo.h"
60#include "grad.h"
61#include "run.h"
62#include "wannier.h"
63
64/** Initialization of Gradient Calculation.
65 * Gradient::GradientArray[GraSchGradient] is set onto "extra" Psis::LocalPsi (needed
66 * for orthogonal conjugate gradient directions!),
67 * for the remaining GradientTypes memory is allocated and zerod.
68 * \param *P Problem at hand
69 * \param *Grad Gradient structure
70 */
71void InitGradients(struct Problem *P, struct Gradient *Grad)
72{
73 struct Lattice *Lat = &P->Lat;
74 struct RunStruct *R = &P->R;
75 struct LatticeLevel *LevS = R->LevS;
76 struct LatticeLevel *ILevS = R->InitLevS;
77 struct Psis *Psi = &Lat->Psi;
78 int i;
79 Grad->GradientArray[GraSchGradient] = LevS->LPsi->LocalPsi[Psi->LocalNo];
80 for (i=1; i<MaxGradientArrayTypes; i++) {
81 Grad->GradientArray[i] = (fftw_complex *) Malloc(sizeof(fftw_complex)*ILevS->MaxG, "InitGradients");
82 SetArrayToDouble0((double*)Grad->GradientArray[i], 2*ILevS->MaxG);
83 }
84}
85
86/** Deallocates memory used in Gradient::GradientArray.
87 * Gradient::GradientArray[GraSchGradient] is set to null.
88 * \param *P Problem at hand
89 * \param *Grad Gradient structure
90 */
91void RemoveGradients(struct Problem *P, struct Gradient *Grad)
92{
93 int i;
94 P=P;
95 Grad->GradientArray[GraSchGradient] = NULL;
96 for (i=1; i<MaxGradientArrayTypes; i++)
97 Free(Grad->GradientArray[i], "RemoveGradients: Grad->GradientArray[i]");
98}
99
100/** Calculates the scalar product of two wave functions with equal symmetry.
101 * \param *P Problem at hand
102 * \param *Lev LatticeLevel structure
103 * \param *LPsiDatA first array of complex wave functions coefficients
104 * \param *LPsiDatB second array of complex wave functions coefficients
105 * \return complex scalar product of the two wave functions.
106 * \warning MPI_Allreduce has to be done still, remember coefficients are shared among processes!
107 */
108double GradSP(struct Problem *P, struct LatticeLevel *Lev, const fftw_complex *LPsiDatA, const fftw_complex *LPsiDatB) {
109 double LocalSP=0.0;
110 int i = 0;
111 if (Lev->GArray[0].GSq == 0.0) {
112 LocalSP += LPsiDatA[0].re * LPsiDatB[0].re;
113 i++;
114 }
115 for (; i < Lev->MaxG; i++) {
116 LocalSP += 2.*(LPsiDatA[i].re * LPsiDatB[i].re + LPsiDatA[i].im * LPsiDatB[i].im); // factor 2 because of the symmetry arising from gamma point!
117 }
118 return(LocalSP);
119}
120
121/** Calculates the scalar product of two wave functions with unequal symmetry.
122 * \param *P Problem at hand
123 * \param *Lev LatticeLevel structure
124 * \param *LPsiDatA first array of complex wave functions coefficients
125 * \param *LPsiDatB second array of complex wave functions coefficients
126 * \return complex scalar product of the two wave functions.
127 * \warning MPI_Allreduce has to be done still, remember coefficients are shared among processes!
128 */
129double GradImSP(struct Problem *P, struct LatticeLevel *Lev, fftw_complex *LPsiDatA, fftw_complex *LPsiDatB) {
130 double LocalSP=0.0;
131 int i = 0;
132 if (Lev->GArray[0].GSq == 0.0) {
133 i++;
134 }
135 for (; i < Lev->MaxG; i++) {
136 LocalSP += 2.*(LPsiDatA[i].re * LPsiDatB[i].im - LPsiDatA[i].im * LPsiDatB[i].re); // factor 2 because of the symmetry arising from gamma point!
137 }
138 return(LocalSP);
139}
140
141/** Calculation of gradient of wave function(without RiemannTensor).
142 * We want to calculate the variational derivative \f$\frac{\delta}{\delta \psi^*_j} E[{\psi^*_j}]\f$,
143 * which we need for the minimisation of the energy functional. In the plane wave basis this evaluates to
144 * \f[
145 * \langle \chi_G| -\frac{1}{2}\nabla^2 + \underbrace{V^H + V^{ps,loc} + V^{XC}}_{V^{loc}} + V^{ps,nl} | \psi_i \rangle
146 * \qquad (section 4.3)
147 * \f]
148 * Therefore, DensityType#HGDensity is back-transformed, exchange correlation potential added up (CalculateXCPotentialNoRT()),
149 * and vector multiplied with the wave function Psi (per node R), transformation into reciprocal grid space and we are back
150 * in the plane wave basis set as desired. CalculateAddNLPot() adds the non-local potential. Finally, the kinetic term
151 * is added up (simple in reciprocal space). The energy eigenvalue \a *Lambda is computed by summation over all reciprocal grid
152 * vectors and in the meantime also the gradient direction (negative derivative) stored in \a *grad while the old values
153 * are moved to \a *oldgrad. At last, \a *Lambda is summed up among all processes.
154 * \f[
155 * \varphi_i^{(m)} (G) = \frac{1}{2} |G|^2 c_{i,G} + \underbrace{V^H (G) + V^{ps,loc} (G)}_{+V^{HG} (G)} + V^{ps,nl} (G) | \psi_i (G) \rangle
156 * \qquad (\textnormal{end of section 4.5.2})
157 * \f]
158 * \param *P Problem at hand, contains Lattice and RunStruct
159 * \param *source source wave function coefficients
160 * \param PsiFactor occupation number, see Psis::PsiFactor
161 * \param *Lambda eigenvalue \f$\lambda_i = \langle \psi_i |h|\psi_i \rangle\f$, see OnePsiElementAddData::Lambda
162 * \param ***fnl non-local form factors over ion type, ion of type and reciprocal grid vector, see PseudoPot::fnl
163 * \param *grad return array for (negative) gradient coefficients \f$\varphi_i^{(m)} (G)\f$
164 * \param *oldgrad return array for (negative) gradient coefficients \f$\varphi_i^{(m-1)} (G)\f$
165 * \param *Hc_grad return array for complex coefficients of \f$H | \psi_i (G) \rangle\f$, see GradientArray#HcGradient
166 * \sa CalculateConDirHConDir() - similar calculation for another explanation.
167 */
168void CalculateGradientNoRT(struct Problem *P, fftw_complex *source, const double PsiFactor, double *Lambda, fftw_complex ***fnl, fftw_complex *grad, fftw_complex *oldgrad, fftw_complex *Hc_grad)
169{
170 struct Lattice *Lat = &P->Lat;
171 struct RunStruct *R = &P->R;
172 struct LatticeLevel *Lev0 = R->Lev0;
173 struct LatticeLevel *LevS = R->LevS;
174 struct Density *Dens0 = Lev0->Dens;
175 struct fft_plan_3d *plan = Lat->plan;
176 int Index, g;
177 fftw_complex *work = Dens0->DensityCArray[TempDensity];
178 fftw_real *HGcR = Dens0->DensityArray[HGcDensity];
179 fftw_complex *HGcRC = (fftw_complex*)HGcR;
180 fftw_complex *HGC = Dens0->DensityCArray[HGDensity];
181 fftw_real *HGCR = (fftw_real *)HGC;
182 fftw_complex *PsiC = Dens0->DensityCArray[ActualPsiDensity];
183 fftw_real *PsiCR = (fftw_real *)PsiC;
184 int nx,ny,nz,iS,i0,s;
185 const int Nx = LevS->Plan0.plan->local_nx;
186 const int Ny = LevS->Plan0.plan->N[1];
187 const int Nz = LevS->Plan0.plan->N[2];
188 const int NUpx = LevS->NUp[0];
189 const int NUpy = LevS->NUp[1];
190 const int NUpz = LevS->NUp[2];
191 double lambda;
192 const double HGcRCFactor = 1./LevS->MaxN;
193 SpeedMeasure(P, LocTime, StartTimeDo);
194 // back-transform HGDensity
195 //if (isnan(HGC[0].re)) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): HGC[%i]_%i.re = NaN!\n", 0, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
196 fft_3d_complex_to_real(plan, Lev0->LevelNo, FFTNF1, HGC, work);
197 // evaluate exchange potential with this density, add up onto HGCR
198 //if (isnan(HGCR[0])) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): HGcR[:%i:]_%i = NaN!\n", 0, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
199 CalculateXCPotentialNoRT(P, HGCR); // add V^{xc} or \epsilon^{xc} upon V^H + V^{ps}
200 for (nx=0;nx<Nx;nx++)
201 for (ny=0;ny<Ny;ny++)
202 for (nz=0;nz<Nz;nz++) {
203 i0 = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
204 iS = nz+Nz*(ny+Ny*nx);
205 HGcR[iS] = HGCR[i0]*PsiCR[i0]; /* Matrix Vector Mult */
206 //if (isnan(PsiCR[i0])) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): PsiCR[%i]_%i = NaN!\n", i0, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
207 //if (isnan(HGCR[i0])) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): HGCR[%i]_%i = NaN!\n", i0, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
208 }
209 fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, HGcRC, work);
210 SpeedMeasure(P, LocTime, StopTimeDo);
211 /* NonLocalPP */
212 //SpeedMeasure(P, NonLocTime, StartTimeDo);
213 CalculateAddNLPot(P, Hc_grad, fnl, PsiFactor); // also resets Hc_grad beforehand
214 //SetArrayToDouble0((double *)Hc_grad,2*LevS->MaxG);
215 SpeedMeasure(P, NonLocTime, StopTimeDo);
216 /* Lambda und Grad */
217 for (g=0;g<LevS->MaxG;g++) {
218 Index = LevS->GArray[g].Index; /* FIXME - factoren */
219 Hc_grad[g].re += PsiFactor*(HGcRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[g].GSq*source[g].re);
220 Hc_grad[g].im += PsiFactor*(HGcRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[g].GSq*source[g].im);
221 //if (isnan(source[g].re)) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): source_%i(%i) = NaN!\n", R->ActualLocalPsiNo, g); Error(SomeError, "NaN-Fehler!"); }
222 }
223 if (grad != NULL) {
224 lambda = 0.0;
225 s = 0;
226 if (LevS->GArray[0].GSq == 0.0) {
227 lambda += Hc_grad[0].re*source[0].re + Hc_grad[0].im*source[0].im;
228 oldgrad[0].re = grad[0].re; // store away old preconditioned steepest decent gradient
229 oldgrad[0].im = grad[0].im;
230 grad[0].re = -Hc_grad[0].re;
231 grad[0].im = -Hc_grad[0].im;
232 s++;
233 }
234 for (g=s;g<LevS->MaxG;g++) {
235 lambda += 2.*(Hc_grad[g].re*source[g].re + Hc_grad[g].im*source[g].im);
236 oldgrad[g].re = grad[g].re; // store away old preconditioned steepest decent gradient
237 oldgrad[g].im = grad[g].im;
238 grad[g].re = -Hc_grad[g].re;
239 grad[g].im = -Hc_grad[g].im;
240 //if (isnan(Hc_grad[g].re)) { fprintf(stderr,"WARNGING: CalculateGradientNoRT(): Hc_grad_%i(%i) = NaN!\n", R->ActualLocalPsiNo, g); Error(SomeError, "NaN-Fehler!"); }
241 }
242 MPI_Allreduce ( &lambda, Lambda, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
243 //fprintf(stderr,"(%i) Lambda = %e\n",P->Par.me, *Lambda);
244 }
245}
246
247
248/** Calculates the Hamiltonian matrix in Kohn-Sham basis set and outputs to file.
249 * In a similiar manner to TestGramSch() all necessary coefficients are retrieved from the
250 * other processes and \f$\langle Psi_i | H | Psi_j \rangle\f$ = \f$\lambda_{i,j}\f$ calculated.
251 * \param *P Problem at hand
252 * \sa CalculateGradientNoRT() - same calculation (of \f$\lambda_i\f$), yet only for diagonal elements
253 */
254void CalculateHamiltonian(struct Problem *P) {
255 struct Lattice *Lat = &P->Lat;
256 struct FileData *F = &P->Files;
257 struct RunStruct *R = &P->R;
258 struct Psis *Psi = &Lat->Psi;
259 struct LatticeLevel *Lev0 = R->Lev0;
260 struct LatticeLevel *LevS = R->LevS;
261 //struct LevelPsi *LPsi = LevS->LPsi;
262 struct Density *Dens0 = Lev0->Dens;
263 int Index;
264 struct fft_plan_3d *plan = Lat->plan;
265 fftw_complex *work = Dens0->DensityCArray[TempDensity];
266 fftw_real *HGcR = Dens0->DensityArray[HGcDensity];
267 fftw_complex *HGcRC = (fftw_complex*)HGcR;
268 fftw_complex *HGC = Dens0->DensityCArray[HGDensity];
269 fftw_real *HGCR = (fftw_real *)HGC;
270 fftw_complex *PsiC = Dens0->DensityCArray[ActualPsiDensity];
271 fftw_real *PsiCR = (fftw_real *)PsiC;
272 fftw_complex ***fnl;
273 fftw_complex *Hc_grad = P->Grad.GradientArray[HcGradient];
274 int nx,ny,nz,iS,i0,s;
275 const int Nx = LevS->Plan0.plan->local_nx;
276 const int Ny = LevS->Plan0.plan->N[1];
277 const int Nz = LevS->Plan0.plan->N[2];
278 const int NUpx = LevS->NUp[0];
279 const int NUpy = LevS->NUp[1];
280 const int NUpz = LevS->NUp[2];
281 double lambda, Lambda;
282 const double HGcRCFactor = 1./LevS->MaxN;
283 double pot_xc, Pot_xc;
284 int PsiFactorA = 1;
285 int PsiFactorB = 1;
286 int NoPsis = 0;
287 int NoOfPsis = 0;
288
289 if (P->Call.out[LeaderOut])
290 fprintf(stderr, "(%i)Setting up Hamiltonian for Level %d...",P->Par.me, LevS->LevelNo);
291 if (P->Call.out[StepLeaderOut])
292 fprintf(stderr, "\n(%i) H = \n",P->Par.me);
293
294 ApplyTotalHamiltonian(P,LevS->LPsi->LocalPsi[0],Hc_grad, P->PP.fnl[0], 1., 1); // update HGDensity->
295
296 // go through all pairs of orbitals
297 int i,k,l,m,j=0,RecvSource;
298 MPI_Status status;
299 struct OnePsiElement *OnePsiB, *OnePsiA, *LOnePsiA, *LOnePsiB;
300 int ElementSize = (sizeof(fftw_complex) / sizeof(double));
301 fftw_complex *LPsiDatA=NULL, *LPsiDatB;
302
303 l = -1;
304 switch (Psi->PsiST) {
305 case SpinDouble:
306 NoOfPsis = NoPsis = Psi->GlobalNo[PsiMaxNoDouble];
307 NoOfPsis += Psi->GlobalNo[PsiMaxAdd];
308 break;
309 case SpinUp:
310 NoOfPsis = NoPsis = Psi->GlobalNo[PsiMaxNoUp];
311 NoOfPsis += Psi->GlobalNo[PsiMaxAdd];
312 break;
313 case SpinDown:
314 NoOfPsis = NoPsis = Psi->GlobalNo[PsiMaxNoDown];
315 NoOfPsis += Psi->GlobalNo[PsiMaxAdd];
316 break;
317 };
318 gsl_matrix *H = gsl_matrix_calloc(NoOfPsis,NoOfPsis);
319 for (i=0; i < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; i++) { // go through all wave functions
320 //fprintf(stderr,"(%i) GlobalNo: %d\tLocalNo: %d\n", P->Par.me,Psi->AllPsiStatus[i].MyGlobalNo,Psi->AllPsiStatus[i].MyLocalNo);
321 OnePsiA = &Psi->AllPsiStatus[i]; // grab OnePsiA
322 if (OnePsiA->PsiType <= UnOccupied) { // drop the extra and perturbed ones
323 l++;
324 if (OnePsiA->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // local?
325 LOnePsiA = &Psi->LocalPsiStatus[OnePsiA->MyLocalNo];
326 else
327 LOnePsiA = NULL;
328 if (LOnePsiA == NULL) { // if it's not local ... receive it from respective process into TempPsi
329 RecvSource = OnePsiA->my_color_comm_ST_Psi;
330 MPI_Recv( LevS->LPsi->TempPsi, LevS->MaxG*ElementSize, MPI_DOUBLE, RecvSource, HamiltonianTag, P->Par.comm_ST_PsiT, &status );
331 LPsiDatA=LevS->LPsi->TempPsi;
332 } else { // .. otherwise send it to all other processes (Max_me... - 1)
333 for (k=0;k<P->Par.Max_me_comm_ST_PsiT;k++)
334 if (k != OnePsiA->my_color_comm_ST_Psi)
335 MPI_Send( LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo], LevS->MaxG*ElementSize, MPI_DOUBLE, k, HamiltonianTag, P->Par.comm_ST_PsiT);
336 LPsiDatA=LevS->LPsi->LocalPsi[OnePsiA->MyLocalNo];
337 } // LPsiDatA is now set to the coefficients of OnePsi either stored or MPI_Received
338
339 if (LOnePsiA != NULL) { // if fnl (!) are locally available!
340 // Trick is: fft wave function by moving it to higher level and there take only each ...th element into account
341 PsiFactorA = 1; //Psi->LocalPsiStatus[OnePsiA->MyLocalNo].PsiFactor;
342 CalculateOneDensityR(Lat, LevS, Dens0, LPsiDatA, Dens0->DensityArray[ActualDensity], R->FactorDensityR*PsiFactorA, 1);
343 // note: factor is not used when storing result in DensityCArray[ActualPsiDensity] in CalculateOneDensityR()!
344 for (nx=0;nx<Nx;nx++)
345 for (ny=0;ny<Ny;ny++)
346 for (nz=0;nz<Nz;nz++) {
347 i0 = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
348 iS = nz+Nz*(ny+Ny*nx);
349 HGcR[iS] = HGCR[i0]*PsiCR[i0]; /* Matrix Vector Mult */
350 }
351 // V^{HG} (R) + V^{XC} (R) | \Psi_b (R) \rangle
352 fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, HGcRC, work);
353 // V^{HG} (G) + V^{XC} (G) | \Psi_b (G) \rangle
354 /* NonLocalPP */
355 fnl = P->PP.fnl[OnePsiA->MyLocalNo]; // fnl only locally available!
356 CalculateAddNLPot(P, Hc_grad, fnl, PsiFactorA); // sets Hc_grad to zero beforehand
357 // V^{ps,nl} (G)| \Psi_b (G) + [V^{HG} (G) + V^{XC} (G)] | \Psi_b (G) \rangle
358 /* Lambda und Grad */
359 s = 0;
360 if (LevS->GArray[0].GSq == 0.0) {
361 Index = LevS->GArray[0].Index; /* FIXME - factoren */
362 Hc_grad[0].re += PsiFactorA*(HGcRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[0].GSq*LPsiDatA[0].re);
363 Hc_grad[0].im += PsiFactorA*(HGcRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[0].GSq*LPsiDatA[0].im);
364 s++;
365 }
366 for (;s<LevS->MaxG;s++) {
367 Index = LevS->GArray[s].Index; /* FIXME - factoren */
368 Hc_grad[s].re += PsiFactorA*(HGcRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[s].GSq*LPsiDatA[s].re);
369 Hc_grad[s].im += PsiFactorA*(HGcRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[s].GSq*LPsiDatA[s].im);
370 }
371 // 1/2 |G|^2 + V^{ps,nl} (G) + V^{HG} (G) + V^{XC} (G) | \Psi_b (G) \rangle
372 } else {
373 SetArrayToDouble0((double *)Hc_grad,2*LevS->MaxG);
374 }
375
376 m = -1;
377 // former border of following loop: Psi->GlobalNo[Psi->PsiST+1]+Psi->GlobalNo[Psi->PsiST+4]+P->Par.Max_me_comm_ST_PsiT
378 for (j=0; j < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; j++) { // go through all wave functions
379 OnePsiB = &Psi->AllPsiStatus[j]; // grab OnePsiA
380 if (OnePsiB->PsiType <= UnOccupied) { // drop the extra ones
381 m++;
382 if (OnePsiB->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // local?
383 LOnePsiB = &Psi->LocalPsiStatus[OnePsiB->MyLocalNo];
384 else
385 LOnePsiB = NULL;
386 if (LOnePsiB == NULL) { // if it's not local ... receive it from respective process into TempPsi
387 RecvSource = OnePsiB->my_color_comm_ST_Psi;
388 MPI_Recv( LevS->LPsi->TempPsi, LevS->MaxG*ElementSize, MPI_DOUBLE, RecvSource, HamiltonianTag, P->Par.comm_ST_PsiT, &status );
389 LPsiDatB=LevS->LPsi->TempPsi;
390 } else { // .. otherwise send it to all other processes (Max_me... - 1)
391 for (k=0;k<P->Par.Max_me_comm_ST_PsiT;k++)
392 if (k != OnePsiB->my_color_comm_ST_Psi)
393 MPI_Send( LevS->LPsi->LocalPsi[OnePsiB->MyLocalNo], LevS->MaxG*ElementSize, MPI_DOUBLE, k, HamiltonianTag, P->Par.comm_ST_PsiT);
394 LPsiDatB=LevS->LPsi->LocalPsi[OnePsiB->MyLocalNo];
395 } // LPsiDatB is now set to the coefficients of OnePsi either stored or MPI_Received
396
397 PsiFactorB = 1; //Psi->LocalPsiStatus[OnePsiB->MyLocalNo].PsiFactor;
398 lambda = 0; //PsiFactorB * GradSP(P, LevS, LPsiDatB, Hc_grad);
399 Lambda = 0;
400 s=0;
401 if (LevS->GArray[0].GSq == 0.0) {
402 lambda += GSL_REAL(gsl_complex_mul(gsl_complex_conjugate(convertComplex(LPsiDatB[s])), convertComplex(Hc_grad[s])));
403 s++;
404 }
405 for (; s < LevS->MaxG; s++) {
406 lambda += GSL_REAL(gsl_complex_mul(gsl_complex_conjugate(convertComplex(LPsiDatB[s])), convertComplex(Hc_grad[s])));
407 lambda += GSL_REAL(gsl_complex_mul(convertComplex(LPsiDatB[s]), gsl_complex_conjugate(convertComplex(Hc_grad[s])))); // c(-G) = c^\ast (G) due to gamma point symmetry !
408 } // due to the symmetry the resulting lambda is both real and symmetric in l,m (i,j) !
409 lambda *= PsiFactorB;
410 // \langle \Psi_a | 1/2 |G|^2 + V^{ps,nl} (G) + V^{HG} (G) + V^{XC} (G) | \Psi_b (G) \rangle
411 MPI_Allreduce ( &lambda, &Lambda, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST); // gather from all(!) (coeffs- and Psis-sharing) groups
412
413 //if (P->Par.me == 0) {
414 // and Output
415 if (P->Call.out[StepLeaderOut])
416 fprintf(stderr,"%2.5lg\t",Lambda);
417 if (P->Par.me == 0) fprintf(F->HamiltonianFile, "%lg\t", Lambda);
418 //}
419 if (m < NoOfPsis && l < NoOfPsis) {
420 //fprintf(stderr,"Index (%i/%i,%i/%i): %e \n",l,NoOfPsis,m,NoOfPsis, Lambda);
421 gsl_matrix_set(H,l,m,Lambda);
422 } else
423 fprintf(stderr,"Index (%i,%i) out of range %i\n",l,m,NoOfPsis);
424 }
425 if ((P->Par.me == 0) && (j == NoOfPsis-1))
426 fprintf(F->HamiltonianFile, "\n");
427 if ((P->Call.out[StepLeaderOut]) && (j == NoOfPsis-1))
428 fprintf(stderr,"\n");
429 }
430 }
431 }
432 if (P->Par.me == 0) {
433 fprintf(F->HamiltonianFile, "\n");
434 fflush(F->HamiltonianFile);
435 }
436 if ((P->Call.out[LeaderOut]) && (!P->Call.out[StepLeaderOut])) {
437 fprintf(stderr,"finished.\n");
438 }
439
440 // storing H matrix for latter use in perturbed calculation
441 for (i=0;i<Psi->NoOfPsis;i++)
442 for (j=0;j<Psi->NoOfPsis;j++)
443 Psi->lambda[i][j] = gsl_matrix_get(H,i,j);
444
445 // retrieve EW from H
446 gsl_vector *eval = gsl_vector_alloc(NoOfPsis);
447 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc(NoOfPsis);
448 gsl_eigen_symm(H, eval, w);
449 gsl_eigen_symm_free(w);
450 gsl_matrix_free(H);
451
452 // print eigenvalues
453 if(P->Call.out[LeaderOut]) {
454 fprintf(stderr,"(%i) Unsorted Eigenvalues:", P->Par.me);
455 for (i=0;i<NoOfPsis;i++)
456 fprintf(stderr,"\t%lg",gsl_vector_get(eval,i));
457 fprintf(stderr,"\n");
458 }
459 gsl_sort_vector(eval); // sort eigenvalues
460 // print eigenvalues
461 if(P->Call.out[LeaderOut]) {
462 fprintf(stderr,"(%i) All Eigenvalues:", P->Par.me);
463 for (i=0;i<NoOfPsis;i++)
464 fprintf(stderr,"\t%lg",gsl_vector_get(eval,i));
465 fprintf(stderr,"\n");
466 fprintf(stderr,"(%i) KS Eigenvalues:", P->Par.me);
467 for (i=0;i<NoOfPsis;i++)
468 fprintf(stderr,"\t%lg",Psi->AddData[i].Lambda);
469 fprintf(stderr,"\n");
470 }
471
472 // print difference between highest occupied and lowest unoccupied
473 if (P->Lat.Psi.GlobalNo[PsiMaxAdd] != 0) {
474 Lat->Energy[UnOccupied].bandgap = gsl_vector_get(eval, Psi->NoOfPsis) - gsl_vector_get(eval, Psi->NoOfPsis-1);
475 //Lat->Energy[UnOccupied].homolumo = gsl_vector_get(eval, Psi->NoOfPsis/2 + 1) - gsl_vector_get(eval, (Psi->NoOfPsis-1)/2);
476 if (P->Call.out[NormalOut])
477 fprintf(stderr,"(%i) Band gap: \\epsilon_{%d+%d}(%d+1) - \\epsilon_{%d+%d}(%d) = %lg Ht\n", P->Par.me, NoPsis,P->Lat.Psi.GlobalNo[PsiMaxAdd],NoPsis, NoPsis, P->Lat.Psi.GlobalNo[PsiMaxAdd], NoPsis, Lat->Energy[UnOccupied].bandgap);
478 //fprintf(stderr,"(%i) HOMO-LUMO gap (HOMO is: %i, LUMO is %i): %lg Ht\n", P->Par.me, Psi->NoOfPsis/2 + 1, (Psi->NoOfPsis-1)/2, Lat->Energy[UnOccupied].homolumo);
479 }
480 // now, calculate \int v_xc (r) n(r) d^3r
481 SetArrayToDouble0((double *)HGCR,2*Dens0->TotalSize);
482 CalculateXCPotentialNoRT(P, HGCR); /* FIXME - kann man halbieren */
483 pot_xc = 0;
484 Pot_xc = 0;
485 for (i=0;i<Dens0->LocalSizeR;i++)
486 pot_xc += HGCR[i]*Dens0->DensityArray[TotalDensity][i];
487 MPI_Allreduce ( &pot_xc, &Pot_xc, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
488 if (P->Call.out[StepLeaderOut]) {
489 fprintf(stderr,"(%i) Exchange-Correlation potential energy: \\int v_xc[n] (r) n(r) d^3r = %lg\n", P->Par.me, HGcRCFactor*Pot_xc*Lat->Volume);
490 }
491 gsl_vector_free(eval);
492}
493
494
495/** Calculation of direction of steepest descent.
496 * The Gradient has been evaluated in CalculateGradient(). Now the SD-direction
497 * \f$\zeta_i = - (H-\lambda_i)\psi_i\f$ is stored in
498 * Gradient::GradientArray[GradientTypes#GraSchGradient].
499 * Afterwards, the "extra" local Psi (which is the above gradient) is orthogonalized
500 * under exclusion of the to be minimialized Psi to ascertain Orthonormality of this
501 * gradient with respect to all other states/orbits.
502 * \param *P Problem at hand, contains Lattice and RunStruct
503 * \param *source source wave function coefficients, \f$\psi_i\f$
504 * \param *grad gradient coefficients, see CalculateGradient(), on return contains conjugate gradient coefficients, \f$\varphi_i^{(m)}\f$
505 * \param *Lambda eigenvalue \f$\lambda_i\f$ = \f$\langle \psi_i^{(m)}|H|\psi_i^{(m)} \rangle\f$
506 * \warning coefficients in \a *grad are overwritten
507 */
508static void CalculateCGGradient(struct Problem *P, fftw_complex *source, fftw_complex *grad, double Lambda)
509{
510 struct Lattice *Lat = &P->Lat;
511 struct RunStruct *R = &P->R;
512 struct Psis *Psi = &Lat->Psi;
513 struct LatticeLevel *LevS = R->LevS;
514 fftw_complex *GradOrtho = P->Grad.GradientArray[GraSchGradient];
515 int g;
516 int ElementSize = (sizeof(fftw_complex) / sizeof(double));
517 for (g=0;g<LevS->MaxG;g++) {
518 GradOrtho[g].re = grad[g].re+Lambda*source[g].re;
519 GradOrtho[g].im = grad[g].im+Lambda*source[g].im;
520 }
521 // include the ExtraPsi (which is the GraSchGradient!)
522 SetGramSchExtraPsi(P, Psi, NotOrthogonal);
523 // exclude the minimised Psi
524 SetGramSchActualPsi(P, Psi, NotUsedToOrtho);
525 SpeedMeasure(P, GramSchTime, StartTimeDo);
526 // makes conjugate gradient orthogonal to all other orbits
527 //fprintf(stderr,"CalculateCGGradient: GramSch() for extra orbital\n");
528 GramSch(P, LevS, Psi, Orthogonalize);
529 SpeedMeasure(P, GramSchTime, StopTimeDo);
530 SetGramSchActualPsi(P, Psi, IsOrthonormal);
531 memcpy(grad, GradOrtho, ElementSize*LevS->MaxG*sizeof(double));
532}
533
534/** Preconditioning for the gradient calculation.
535 * A suitable matrix for preconditioning is [TPA89]
536 * \f[
537 * M_{G,G'} = \delta_{G,G'} \frac{27+18x+12x^2+8x^3}{27+18x+12x^2+8x^3+16x^4} \qquad (5.3)
538 * \f]
539 * where \f$x= \frac{\langle \xi_G|-\frac{1}{2}\nabla^2|\xi_G \rangle}
540 * {\langle \psi_i^{(m)}|-\frac{1}{2}\nabla^2|\psi_i^{(m)} \rangle}\f$.
541 * (This simple form is due to the diagonal dominance of the kinetic term in Hamiltonian, and
542 * x is much simpler to evaluate in reciprocal space!)
543 * With this one the higher plane wave states are being transformed such that they
544 * nearly equal the error vector and converge at an equal rate.
545 * So that the new steepest descent direction is \f$M \tilde{\zeta}_i\f$.
546 * Afterwards, the new direction is orthonormalized via GramSch() and the "extra"
547 * local Psi in a likewise manner to CalculateCGGradient().
548 *
549 * \param *P Problem at hand, contains LatticeLevel and RunStruct
550 * \param *grad array of steepest descent coefficients, \f$\varphi_i^{(m)}\f$
551 * \param *PCgrad return array for preconditioned steepest descent coefficients, \f$\tilde{\varphi}_i^{(m)}\f$
552 * \param T kinetic eigenvalue \f$\langle \psi_i | - \frac{1}{2} \nabla^2 | \psi_i \rangle\f$
553 */
554static void CalculatePreConGrad(struct Problem *P, fftw_complex *grad, fftw_complex *PCgrad, double T)
555{
556 struct RunStruct *R = &P->R;
557 struct LatticeLevel *LevS = R->LevS;
558 int g;
559 double x, K, dK;
560 struct Psis *Psi = &P->Lat.Psi;
561 fftw_complex *PCOrtho = P->Grad.GradientArray[GraSchGradient];
562 int ElementSize = (sizeof(fftw_complex) / sizeof(double));
563 if (fabs(T) < MYEPSILON) T = 1;
564 for (g=0;g<LevS->MaxG;g++) {
565 x = LevS->GArray[g].GSq;
566 x /= T;
567 K = (((8*x+12)*x +18)*x+27);
568 dK = ((((16*x+8)*x +12)*x+18)*x+27);
569 K /= dK;
570 c_re(PCOrtho[g]) = K*c_re(grad[g]);
571 c_im(PCOrtho[g]) = K*c_im(grad[g]);
572 }
573 SetGramSchExtraPsi(P, Psi, NotOrthogonal);
574 SpeedMeasure(P, GramSchTime, StartTimeDo);
575 // preconditioned direction is orthogonalized
576 //fprintf(stderr,"CalculatePreConGrad: GramSch() for extra orbital\n");
577 GramSch(P, LevS, Psi, Orthogonalize);
578 SpeedMeasure(P, GramSchTime, StopTimeDo);
579 memcpy(PCgrad, PCOrtho, ElementSize*LevS->MaxG*sizeof(double));
580}
581
582/** Calculation of conjugate direction.
583 * The new conjugate direction is according to the algorithm of Fletcher-Reeves:
584 * \f[
585 * \varphi^{(m)}_i =
586 * \tilde{\eta}^{(m)}_i + \frac{\langle \tilde{\eta}^{(m)}_i | \tilde{\varphi}_i^{(m)} \rangle}{\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle} \varphi_i^{(m-1)}
587 * \f]
588 * Instead we use the algorithm of Polak-Ribiere:
589 * \f[
590 * \varphi^{(m)}_i =
591 * \tilde{\eta}^{(m)}_i + \frac{\langle \tilde{\eta}^{(m)}_i | \bigl (\tilde{\varphi}_i^{(m)} \rangle - \tilde{\varphi}_i^{(m-1)} \rangle \bigr)}{\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle} \varphi_i^{(m-1)}
592 * \f]
593 * with an additional automatic restart:
594 * \f[
595 * \varphi^{(m)}_i = \tilde{\eta}^{(m)}_i \quad \textnormal{if}\ \frac{\langle \tilde{\eta}^{(m)}_i | \bigl (\tilde{\varphi}_i^{(m)} \rangle - \tilde{\varphi}_i^{(m-1)} \rangle \bigr)}{\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle} <= 0
596 * \f]
597 * That's why ConDirGradient and OldConDirGradient are stored, the scalar product - summed up among all processes -
598 * is saved in OnePsiAddData::Gamma. The conjugate direction is again orthonormalized by GramSch() and stored
599 * in \a *ConDir.
600 * \param *P Problem at hand, containing LatticeLevel and RunStruct
601 * \param step minimisation step m
602 * \param *GammaDiv scalar product \f$\langle \tilde{\eta}^{(m-1)}_i | \tilde{\varphi}_i^{(m-1)} \rangle\f$, on calling expects (m-1), on return contains (m)
603 * \param *grad array for coefficients of steepest descent direction \f$\tilde{\zeta}_i^{(m)}\f$
604 * \param *oldgrad array for coefficients of steepest descent direction of \a step-1 \f$\tilde{\zeta}_i^{(m-1)}\f$
605 * \param *PCgrad array for coefficients of preconditioned steepest descent direction \f$\tilde{\eta}_i^{(m)}\f$
606 * \param *ConDir return array for coefficients of conjugate array of \a step, \f$\tilde{\varphi}_i^{(m)}\f$ = \f$\frac{\tilde{\varphi}_i^{(m)}} {||\tilde{\varphi}_i^{(m)}||}\f$
607 * \param *ConDir_old array for coefficients of conjugate array of \a step-1, \f$\tilde{\varphi}_i^{(m-1)}\f$
608 */
609static void CalculateConDir(struct Problem *P, int step, double *GammaDiv,
610 fftw_complex *grad,
611 fftw_complex *oldgrad,
612 fftw_complex *PCgrad,
613 fftw_complex *ConDir, fftw_complex *ConDir_old) {
614 struct RunStruct *R = &P->R;
615 struct LatticeLevel *LevS = R->LevS;
616 struct Psis *Psi = &P->Lat.Psi;
617 fftw_complex *Ortho = P->Grad.GradientArray[GraSchGradient];
618 int ElementSize = (sizeof(fftw_complex) / sizeof(double));
619 int g;
620 double S[3], Gamma, GammaDivOld = *GammaDiv;
621 double LocalSP[2], PsiSP[2];
622 LocalSP[0] = GradSP(P, LevS, PCgrad, grad);
623 LocalSP[1] = GradSP(P, LevS, PCgrad, oldgrad);
624
625 MPI_Allreduce ( LocalSP, PsiSP, 2, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
626 *GammaDiv = S[0] = PsiSP[0];
627 S[1] = GammaDivOld;
628 S[2] = PsiSP[1];
629 if (step) { // only in later steps is the scalar product used, but always condir stored in oldcondir and Ortho (working gradient)
630 if (fabs(S[1]) < MYEPSILON) {
631 if (fabs((S[0]-S[2])) < MYEPSILON)
632 Gamma = 1.0;
633 else
634 Gamma = 0.0;
635 } else
636 Gamma = (S[0]-S[2])/S[1]; // the final CG beta coefficient (Polak-Ribiere)
637 if (Gamma < 0) { // Polak-Ribiere with automatic restart: gamma = max {0, gamma}
638 Gamma = 0.;
639 if (P->Call.out[StepLeaderOut]) fprintf(stderr,"(%i) Polak-Ribiere: Restarting iteration by setting gamma = 0.\n", P->Par.me);
640 }
641 for (g=0; g < LevS->MaxG; g++) {
642 c_re(ConDir[g]) = c_re(PCgrad[g]) + Gamma*c_re(ConDir_old[g]);
643 c_im(ConDir[g]) = c_im(PCgrad[g]) + Gamma*c_im(ConDir_old[g]);
644 c_re(ConDir_old[g]) = c_re(ConDir[g]);
645 c_im(ConDir_old[g]) = c_im(ConDir[g]);
646 c_re(Ortho[g]) = c_re(ConDir[g]);
647 c_im(Ortho[g]) = c_im(ConDir[g]);
648 }
649 } else {
650 Gamma = 0.0;
651 for (g=0; g < LevS->MaxG; g++) {
652 c_re(ConDir[g]) = c_re(PCgrad[g]);
653 c_im(ConDir[g]) = c_im(PCgrad[g]);
654 c_re(ConDir_old[g]) = c_re(ConDir[g]);
655 c_im(ConDir_old[g]) = c_im(ConDir[g]);
656 c_re(Ortho[g]) = c_re(ConDir[g]);
657 c_im(Ortho[g]) = c_im(ConDir[g]);
658 }
659 }
660 // orthonormalize
661 SetGramSchExtraPsi(P, Psi, NotOrthogonal);
662 SpeedMeasure(P, GramSchTime, StartTimeDo);
663 //fprintf(stderr,"CalculateConDir: GramSch() for extra orbital\n");
664 GramSch(P, LevS, Psi, Orthonormalize);
665 SpeedMeasure(P, GramSchTime, StopTimeDo);
666 memcpy(ConDir, Ortho, ElementSize*LevS->MaxG*sizeof(double));
667}
668
669/** Calculation of \f$\langle \widetilde{\widetilde{\varphi}}_i^{(m)}|H|\widetilde{\widetilde{\varphi}}_i^{(m)} \rangle\f$.
670 * Evaluates most of the parts of the second derivative of the energy functional \f$\frac{\delta^2 E}{\delta \Theta^2}|_{\Theta=0}\f$,
671 * most importantly the energy expection value of the conjugate direction.
672 * \f[
673 * \langle \widetilde{\widetilde{\varphi}}_i^{(m)} (G) | -\frac{1}{2} \nabla^2 + \underbrace{V^H (G) + V^{ps,loc} (G)}_{=V^{HG}(G)} + V^{ps,nl} | \widetilde{\widetilde{\varphi}} (G) \rangle
674 * \f]
675 * The conjugate direction coefficients \f$\varphi_i^{(m)} (G)\f$ are fftransformed to \f$\varphi_i^{(m)} (R)\f$,
676 * times the DoubleDensityTypes#HGDensity \f$V^{HG} (R) \f$we have the first part,
677 * then \f$\rho (R) = 2\frac{f_i}{V} \Re (\varphi_i^{(m)} (R) \cdot \psi_i (R))\f$is evaluated,
678 * HGDensity is fftransformed to \f$V^{HG} | \varphi_i^{(m)} (G) \rangle\f$,
679 * CalculateAddNLPot() adds non-local potential \f$V^{ps,nl} (G)\f$ to it,
680 * kinetic term is added in the reciprocal base as well \f$-\frac{1}{2} |G|^2\f$,
681 * and at last scalar product between it and the conjugate direction, and done.
682 *
683 * Meanwhile the fftransformed ConDir density \f$\rho (R)\f$ is used for the ddEddt0s \f$A_H\f$, \f$A_{XC}\f$.
684 * In case of RunStruct::DoCalcCGGauss in the following expression additionally \f$\hat{n}^{Gauss}\f$
685 * is evaluted, otherwise left out:
686 * \f[
687 * A^H = \frac{4\pi}{V} \sum_G \frac{|\frac{\rho(G)}{V^2 N} + \hat{n}^{Gauss}|}{|G|^2}
688 * \qquad (\textnormal{section 5.1, line search})
689 * \f]
690 * \param *P Problem at hand
691 * \param *ConDir conjugate direction \f$\tilde{\tilde{\varphi}}_i^{(m)}\f$
692 * \param PsiFactor occupation number, see Psis::PsiFactor
693 * \param *ConDirHConDir first return value: energy expectation value for conjugate direction
694 * \param *HartreeddEddt0 second return value: second derivative of Hartree energy \f$A_H\f$
695 * \param *XCddEddt0 third return value: second derivative of exchange correlation energy \f$A_{XC}\f$
696 * \warning The MPI_Allreduce for the scalar product in the end has not been done!
697 * \sa CalculateLineSearch() - used there
698 */
699void CalculateConDirHConDir(struct Problem *P, fftw_complex *ConDir, const double PsiFactor,
700 double *ConDirHConDir, double *HartreeddEddt0, double *XCddEddt0)
701{
702 struct Lattice *Lat = &P->Lat;
703 struct PseudoPot *PP = &P->PP;
704 struct RunStruct *R = &P->R;
705 struct LatticeLevel *Lev0 = R->Lev0;
706 struct LatticeLevel *LevS = R->LevS;
707 struct Density *Dens0 = Lev0->Dens;
708 int Index, g, i, pos, s=0;
709 struct fft_plan_3d *plan = Lat->plan;
710 fftw_complex *work = Dens0->DensityCArray[TempDensity];
711 fftw_real *PsiCD = (fftw_real *)work;
712 fftw_complex *PsiCDC = work;
713 fftw_real *PsiR = (fftw_real *)Dens0->DensityCArray[ActualPsiDensity];
714 fftw_complex *tempdestRC = (fftw_complex *)Dens0->DensityArray[TempDensity];
715 fftw_real *HGCDR = Dens0->DensityArray[HGcDensity];
716 fftw_complex *HGCDRC = (fftw_complex*)HGCDR;
717 fftw_complex *HGC = Dens0->DensityCArray[HGDensity];
718 fftw_real *HGCR = (fftw_real *)HGC;
719 fftw_complex *posfac, *destpos, *destRCS, *destRCD;
720 const double HartreeFactorR = 2.*R->FactorDensityR*PsiFactor;
721 double HartreeFactorC = 4.*PI*Lat->Volume;
722 double HartreeFactorCGauss = 1./Lev0->MaxN; //R->FactorDensityC*R->HGcFactor/Lev0->MaxN; // same expression as first and second factor cancel each other
723 int nx,ny,nz,iS,i0;
724 const int Nx = LevS->Plan0.plan->local_nx;
725 const int Ny = LevS->Plan0.plan->N[1];
726 const int Nz = LevS->Plan0.plan->N[2];
727 const int NUpx = LevS->NUp[0];
728 const int NUpy = LevS->NUp[1];
729 const int NUpz = LevS->NUp[2];
730 const double HGcRCFactor = 1./LevS->MaxN;
731 fftw_complex *HConDir = P->Grad.GradientArray[TempGradient];
732 int DoCalcCGGauss = R->DoCalcCGGauss;
733 fftw_complex rp, rhog;
734 int it;
735 struct Ions *I = &P->Ion;
736
737 /* ConDir H ConDir */
738 // (G)->(R): retrieve density from wave functions into PsiCD in the usual manner (yet within the same level!) (see CalculateOneDensityR())
739 SetArrayToDouble0((double *)tempdestRC, Dens0->TotalSize*2);
740 for (i=0;i<LevS->MaxG;i++) {
741 Index = LevS->GArray[i].Index;
742 posfac = &LevS->PosFactorUp[LevS->MaxNUp*i];
743 destpos = &tempdestRC[LevS->MaxNUp*Index];
744 for (pos=0; pos < LevS->MaxNUp; pos++) {
745 destpos[pos].re = ConDir[i].re*posfac[pos].re-ConDir[i].im*posfac[pos].im;
746 destpos[pos].im = ConDir[i].re*posfac[pos].im+ConDir[i].im*posfac[pos].re;
747 }
748 }
749 for (i=0; i<LevS->MaxDoubleG; i++) {
750 destRCS = &tempdestRC[LevS->DoubleG[2*i]*LevS->MaxNUp];
751 destRCD = &tempdestRC[LevS->DoubleG[2*i+1]*LevS->MaxNUp];
752 for (pos=0; pos < LevS->MaxNUp; pos++) {
753 destRCD[pos].re = destRCS[pos].re;
754 destRCD[pos].im = -destRCS[pos].im;
755 }
756 }
757 fft_3d_complex_to_real(plan, LevS->LevelNo, FFTNFUp, tempdestRC, work);
758 DensityRTransformPos(LevS,(fftw_real*)tempdestRC, PsiCD);
759 // HGCR contains real(!) HGDensity from CalculateGradientNoRT() !
760 for (nx=0;nx<Nx;nx++)
761 for (ny=0;ny<Ny;ny++)
762 for (nz=0;nz<Nz;nz++) {
763 i0 = nz*NUpz+Nz*NUpz*(ny*NUpy+Ny*NUpy*nx*NUpx);
764 iS = nz+Nz*(ny+Ny*nx);
765 HGCDR[iS] = HGCR[i0]*PsiCD[i0]; /* Matrix Vector Mult */
766 }
767 // now we have V^HG + V^XC | \varphi \rangle
768 /*
769 Hartree & Exchange Correlaton ddt0
770 */
771 for (i=0; i<Dens0->LocalSizeR; i++)
772 PsiCD[i] *= HartreeFactorR*PsiR[i]; // is rho(r) = PsiFactor times \psi_i is in PsiR times 1/V times \varphi_i^{(m)}
773 *XCddEddt0 = CalculateXCddEddt0NoRT(P, PsiCD); // calcs A^{XC}
774 fft_3d_real_to_complex(plan, Lev0->LevelNo, FFTNF1, PsiCDC, tempdestRC);
775 s = 0;
776 *HartreeddEddt0 = 0.0;
777 if (Lev0->GArray[0].GSq == 0.0)
778 s++;
779 if (DoCalcCGGauss) { // this part is wrong, n^Gauss is not part of HartreeddEddt0 (only appears in 2*(<phi_i| H |phi_i> - ...))
780 for (g=s; g < Lev0->MaxG; g++) {
781 Index = Lev0->GArray[g].Index;
782 rp.re = 0.0; rp.im = 0.0;
783 for (it = 0; it < I->Max_Types; it++) {
784 c_re(rp) += (c_re(I->I[it].SFactor[g])*PP->FacGauss[it][g]);
785 c_im(rp) += (c_im(I->I[it].SFactor[g])*PP->FacGauss[it][g]);
786 }
787 c_re(rhog) = c_re(PsiCDC[Index])*HartreeFactorCGauss+c_re(rp);
788 c_im(rhog) = c_im(PsiCDC[Index])*HartreeFactorCGauss+c_im(rp);
789 if (Lev0->GArray[g].GSq < MYEPSILON) fprintf(stderr,"CalculateConDirHConDir: Lev0->GArray[g].GSq = %lg",Lev0->GArray[g].GSq);
790 *HartreeddEddt0 += (c_re(rhog)*c_re(rhog)+c_im(rhog)*c_im(rhog))/(Lev0->GArray[g].GSq);
791 }
792 } else {
793 HartreeFactorC *= HartreeFactorCGauss*HartreeFactorCGauss;
794 for (g=s; g < Lev0->MaxG; g++) {
795 Index = Lev0->GArray[g].Index;
796 if (Lev0->GArray[g].GSq < MYEPSILON) fprintf(stderr,"CalculateConDirHConDir: Lev0->GArray[g].GSq = %lg",Lev0->GArray[g].GSq);
797 *HartreeddEddt0 += (c_re(PsiCDC[Index])*c_re(PsiCDC[Index])+c_im(PsiCDC[Index])*c_im(PsiCDC[Index]))/(Lev0->GArray[g].GSq);
798 }
799 }
800 *HartreeddEddt0 *= HartreeFactorC;
801 /* fertig mit ddEddt0 */
802
803 fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, HGCDRC, work);
804 CalculateCDfnl(P, ConDir, PP->CDfnl); // calculate needed non-local form factors
805 CalculateAddNLPot(P, HConDir, PP->CDfnl, PsiFactor); // add non-local potential, sets to zero beforehand
806 //SetArrayToDouble0((double *)HConDir,2*LevS->MaxG);
807 // now we have V^{HG} + V^{XC} + V^{ps,nl} | \varphi \rangle
808 for (g=0;g<LevS->MaxG;g++) { // factor by one over volume and add kinetic part
809 Index = LevS->GArray[g].Index; /* FIXME - factoren */
810 HConDir[g].re += PsiFactor*(HGCDRC[Index].re*HGcRCFactor + 0.5*LevS->GArray[g].GSq*ConDir[g].re);
811 HConDir[g].im += PsiFactor*(HGCDRC[Index].im*HGcRCFactor + 0.5*LevS->GArray[g].GSq*ConDir[g].im);
812 }
813 // now, T^{kin} + V^{HG} + V^{XC} + V^{ps,nl} | \varphi \rangle
814 *ConDirHConDir = GradSP(P, LevS, ConDir, HConDir);
815 // finally, \langle \varphi | V^{HG} + V^{XC} + V^{ps,nl} | \varphi \rangle which is returned
816}
817
818/** Find the minimal \f$\Theta_{min}\f$ for the line search.
819 * We find it as
820 * \f[
821 * \widetilde{\Theta}_{min} = \frac{1}{2} \tan^{-1} \Bigl ( -\frac{ \frac{\delta E}{\delta \Theta}_{|\Theta =0} }
822 * { \frac{1}{2} \frac{\delta^2 E}{\delta^2 \Theta}_{|\Theta =0} } \Bigr )
823 * \f]
824 * Thus, with given \f$E(0)\f$, \f$\frac{\delta E}{\delta \Theta}_{|\Theta =0}\f$ and \f$\frac{\delta^2 E}{\delta^2 \Theta}_{|\Theta =0}\f$
825 * \f$\Theta_{min}\f$is calculated and the unknowns A, B and C (see CalculateLineSearch()) obtained. So that our
826 * approximating formula \f$\widetilde{E}(\Theta) = A + B cos(2\Theta) + CSin(2\Theta)\f$ is fully known and we
827 * can evaluate the energy functional and its derivative at \f$\Theta_{min} \f$.
828 * \param E0 energy eigenvalue at \f$\Theta = 0\f$: \f$A= E(0) - B\f$
829 * \param dEdt0 first derivative of energy eigenvalue at \f$\Theta = 0\f$: \f$C = \frac{1}{2} \frac{\delta E}{\delta \Theta}|_{\Theta=0}\f$
830 * \param ddEddt0 first derivative of energy eigenvalue at \f$\Theta = 0\f$: \f$B = -\frac{1}{4} \frac{\delta^2 E}{\delta \Theta^2}|_{\Theta=0}\f$
831 * \param *E returned energy functional at \f$\Theta_{min}\f$
832 * \param *dE returned first derivative energy functional at \f$\Theta_{min}\f$
833 * \param *ddE returned second derivative energy functional at \f$\Theta_{min}\f$
834 * \param *dcos cosine of the found minimal angle
835 * \param *dsin sine of the found minimal angle
836 * \return \f$\Theta_{min}\f$
837 * \sa CalculateLineSearch() - used there
838 */
839double CalculateDeltaI(double E0, double dEdt0, double ddEddt0, double *E, double *dE, double *ddE, double *dcos, double *dsin)
840{
841 double delta, dcos2, dsin2, A, B, Eavg;
842 if (fabs(ddEddt0) < MYEPSILON) fprintf(stderr,"CalculateDeltaI: ddEddt0 = %lg",ddEddt0);
843 delta = atan(-2.*dEdt0/ddEddt0)/2.0; //.5*
844 *dcos = cos(delta);
845 *dsin = sin(delta);
846 dcos2 = cos(delta*2.0); //2.*
847 dsin2 = sin(delta*2.0); //2.*
848 A = -ddEddt0/4.;
849 B = dEdt0/2.;
850 //fprintf(stderr,"CalculateDeltaI: A = %lg, B = %lg\n", A, B);
851 Eavg = E0-A;
852 *E = Eavg + A*dcos2 + B*dsin2;
853 *dE = -2.*A*dsin2+2.*B*dcos2;
854 *ddE = -4.*A*dcos2-4.*B*dsin2;
855 if (*ddE < 0) { // if it's a maximum (indicated by positive second derivative at found theta)
856 if (delta < 0) { // go the minimum - the two are always separated by PI/2 by the ansatz
857 delta += PI/2.0; ///2.
858 } else {
859 delta -= PI/2.0; ///2.
860 }
861 *dcos = cos(delta);
862 *dsin = sin(delta);
863 dcos2 = cos(delta*2.0); //2.*
864 dsin2 = sin(delta*2.0); //2.*
865 *E = Eavg + A*dcos2 + B*dsin2;
866 *dE = -2.*A*dsin2+2.*B*dcos2;
867 *ddE = -4.*A*dcos2-4.*B*dsin2;
868 }
869 return(delta);
870}
871
872
873/** Line search: One-dimensional minimisation along conjugate direction.
874 * We want to minimize the energy functional along the calculated conjugate direction.
875 * States \f$\psi_i^{(m)}(\Theta) = \widetilde{\psi}_i^{(m)} \cos(\Theta) + \widetilde{\widetilde{\varphi}}_i^{(m)}
876 * \sin(\Theta)\f$ are normalized and orthogonal for every \f$\Theta \in {\cal R}\f$. Thus we
877 * want to find the minimal \f$\Theta_{min}\f$. Intenting to evaluate the energy functional as
878 * seldom as possible, we approximate it by:
879 * \f$\widetilde{E}(\Theta) = A + B cos(2\Theta) + CSin(2\Theta)\f$
880 * where the functions shall match the energy functional and its first and second derivative at
881 * \f$\Theta = 0\f$, fixing the unknowns A, B and C.\n
882 * \f$\frac{\delta E}{\delta \Theta}_{|\Theta =0} =
883 * 2 \Re \bigl ( \langle \widetilde{\widetilde{\varphi}}_i^{(m)}|H|\psi_i^{(m)} \rangle \bigr )\f$
884 * and
885 * \f$\frac{\delta^2 E}{\delta^2 \Theta}_{|\Theta =0} =
886 * 2\bigl( \langle \widetilde{\widetilde{\varphi}}_i^{(m)}|H|\widetilde{\widetilde{\varphi}}_i^{(m)} \rangle
887 * - \langle \psi_i^{(m)}|H|\psi_i^{(m)} \rangle
888 * + \underbrace{\int_\Omega \nu_H[\rho](r)\rho(r)dr}_{A_H}
889 * + \underbrace{\int_\Omega \rho(r)^2 \frac{\delta V^{xc}}{\delta n}(r) dr}_{A_{xc}}\f$. \n
890 * \f$A_H\f$ and \f$A_{xc}\f$ can be calculated from the Hartree respectively the exchange
891 * correlation energy.\n
892 * Basically, this boils down to using CalculateConDirHConDir() (we already have OnePsiAddData::Lambda),
893 * in evaluation of the energy functional and its derivatives as shown in the above formulas
894 * and calling with these CalculateDeltaI() which finds the above \f$\Theta_{min} \f$. \n
895 * The coefficients of the current wave function are recalculated, thus we have the new wave function,
896 * the total energy and its derivatives printed to stderr.\n
897 *
898 * \param *P Problem at hand
899 * \param *source current minimised wave function
900 * \param *ConDir conjugate direction, \f$\tilde{\tilde{\varphi}}_i^{(m)}\f$
901 * \param *Hc_grad complex coefficients of \f$H | \psi_i (G) \rangle\f$, see GradientArray#HcGradient
902 * \param x if NULL do approximative line search by second order taylor expansion, otherwise use given
903 * parameter between -M_PI...M_PI. Is simply put through to CalculateLineSearch()
904 * \param PsiFactor occupation number
905 */
906static void CalculateLineSearch(struct Problem *P,
907 fftw_complex *source,
908 fftw_complex *ConDir,
909 fftw_complex *Hc_grad, const double PsiFactor, const double *x)
910{
911 struct Lattice *Lat = &P->Lat;
912 struct RunStruct *R = &P->R;
913 struct LatticeLevel *LevS = R->LevS;
914 struct FileData *F = &P->Files;
915 //struct LatticeLevel *Lev0 = R->Lev0;
916 //struct Density *Dens = Lev0->Dens;
917 //struct Psis *Psi = &P->Lat.Psi;
918 struct Energy *En = Lat->E;
919 //int ElementSize = (sizeof(fftw_complex) / sizeof(double));
920 double dEdt0, ddEddt0, ConDirHConDir, Eavg, A, B;//, sourceHsource;
921 double E0, E, dE, ddE, delta, dcos, dsin, dcos2, dsin2;
922 double EI, dEI, ddEI, deltaI, dcosI, dsinI;
923 double HartreeddEddt0, XCddEddt0;
924 double d[4],D[4], Diff;
925 //double factorDR = R->FactorDensityR;
926 int g, i;
927 //int ActNum;
928
929 // ========= dE / dt | 0 ============
930 //fprintf(stderr,"(%i) |ConDir| = %lg, |Hc_grad| = %lg\n", P->Par.me, sqrt(GradSP(P, LevS, ConDir, ConDir)), sqrt(GradSP(P, LevS, Hc_grad, Hc_grad)));
931 d[0] = 2.*GradSP(P, LevS, ConDir, Hc_grad);
932
933 // ========== ddE / ddt | 0 =========
934 CalculateConDirHConDir(P, ConDir, PsiFactor, &d[1], &d[2], &d[3]);
935
936 MPI_Allreduce ( &d, &D, 4, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
937 dEdt0 = D[0];
938 for (i=MAXOLD-1; i > 0; i--)
939 En->dEdt0[i] = En->dEdt0[i-1];
940 En->dEdt0[0] = dEdt0;
941 ConDirHConDir = D[1];
942 HartreeddEddt0 = D[2];
943 XCddEddt0 = D[3];
944 ddEddt0 = 0.0;
945 switch (R->CurrentMin) {
946 case Occupied:
947 ddEddt0 = HartreeddEddt0+XCddEddt0; // these terms drop out in unoccupied variation due to lacking dependence of n^{tot}
948 case UnOccupied:
949 ddEddt0 += 2.*(ConDirHConDir - Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda);
950 break;
951 default:
952 ddEddt0 = 0.0;
953 break;
954 }
955 //if (R->CurrentMin <= UnOccupied)
956 //fprintf(stderr,"ConDirHConDir = %lg\tLambda = %lg\t Hartree = %lg\t XC = %lg\n",ConDirHConDir,Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda,HartreeddEddt0,XCddEddt0);
957
958 for (i=MAXOLD-1; i > 0; i--)
959 En->ddEddt0[i] = En->ddEddt0[i-1];
960 En->ddEddt0[0] = ddEddt0;
961 E0 = En->TotalEnergy[0];
962
963 if (x == NULL) {
964 // delta
965 //if (isnan(E0)) { fprintf(stderr,"(%i) WARNING in CalculateLineSearch(): E0_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
966 //if (isnan(dEdt0)) { fprintf(stderr,"(%i) WARNING in CalculateLineSearch(): dEdt0_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
967 //if (isnan(ddEddt0)) { fprintf(stderr,"(%i) WARNING in CalculateLineSearch(): ddEddt0_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
968 deltaI = CalculateDeltaI(E0, dEdt0, ddEddt0,
969 &EI, &dEI, &ddEI, &dcosI, &dsinI);
970 } else {
971 deltaI = *x;
972 dcosI = cos(deltaI);
973 dsinI = sin(deltaI);
974 dcos2 = cos(deltaI*2.0); //2.*
975 dsin2 = sin(deltaI*2.0); //2.*
976 A = -ddEddt0/4.;
977 B = dEdt0/2.;
978 //fprintf(stderr,"CalculateDeltaI: A = %lg, B = %lg\n", A, B);
979 Eavg = E0-A;
980 EI = Eavg + A*dcos2 + B*dsin2;
981 dEI = -2.*A*dsin2+2.*B*dcos2;
982 ddEI = -4.*A*dcos2-4.*B*dsin2;
983 }
984 delta = deltaI; E = EI; dE = dEI; ddE = ddEI; dcos = dcosI; dsin = dsinI;
985 if (R->ScanPotential >= R->ScanAtStep) {
986 delta = (R->ScanPotential-(R->ScanAtStep-1))*R->ScanDelta;
987 dcos = cos(delta);
988 dsin = sin(delta);
989 //fprintf(stderr,"(%i) Scaning Potential form with delta %lg\n", P->Par.me, delta);
990 }
991 // shift energy delta values
992 for (i=MAXOLD-1; i > 0; i--) {
993 En->delta[i] = En->delta[i-1];
994 En->ATE[i] = En->ATE[i-1];
995 }
996 // store new one
997 En->delta[0] = delta;
998 En->ATE[0] = E;
999
1000 if (R->MinStep != 0)
1001 Diff = fabs(En->TotalEnergy[1] - E0)/(En->TotalEnergy[1] - E0)*fabs((E0 - En->ATE[1])/E0);
1002 else
1003 Diff = 0.;
1004
1005 // rotate local Psi and ConDir according to angle Theta
1006 for (g = 0; g < LevS->MaxG; g++) { // Here all coefficients are updated for the new found wave function
1007 //if (isnan(ConDir[g].re)) { fprintf(stderr,"WARNGING: CalculateLineSearch(): ConDir_%i(%i) = NaN!\n", R->ActualLocalPsiNo, g); Error(SomeError, "NaN-Fehler!"); }
1008 c_re(source[g]) = c_re(source[g])*dcos + c_re(ConDir[g])*dsin;
1009 c_im(source[g]) = c_im(source[g])*dcos + c_im(ConDir[g])*dsin;
1010 }
1011 // print to stderr
1012 if (P->Call.out[StepLeaderOut] && x == NULL) {
1013 if (1)
1014 switch (R->CurrentMin) {
1015 case UnOccupied:
1016 fprintf(stderr, "(%i,%i,%i)S(%i,%i,%i):\tTGE: %e\tATGE: %e\t Diff: %e\t --- d: %e\tdEdt0: %e\tddEddt0: %e\n",P->Par.my_color_comm_ST,P->Par.me_comm_ST, P->Par.me_comm_ST_PsiT, R->MinStep, R->ActualLocalPsiNo, R->PsiStep, Lat->Energy[UnOccupied].TotalEnergy[0], E, Diff, delta, dEdt0, ddEddt0);
1017 break;
1018 default:
1019 fprintf(stderr, "(%i,%i,%i)S(%i,%i,%i):\tTE: %e\tATE: %e\t Diff: %e\t --- d: %e\tdEdt0: %e\tddEddt0: %e\n",P->Par.my_color_comm_ST,P->Par.me_comm_ST, P->Par.me_comm_ST_PsiT, R->MinStep, R->ActualLocalPsiNo, R->PsiStep, Lat->Energy[Occupied].TotalEnergy[0], E, Diff, delta, dEdt0, ddEddt0);
1020 break;
1021 }
1022 }
1023 if ((P->Par.me == 0) && (R->ScanFlag || (R->ScanPotential < R->ScanAtStep))) {
1024 fprintf(F->MinimisationFile, "%i\t%i\t%i\t%e\t%e\t%e\t%e\t%e\n",R->MinStep, R->ActualLocalPsiNo, R->PsiStep, E0, E, delta, dEdt0, ddEddt0);
1025 fflush(F->MinimisationFile);
1026 }
1027}
1028
1029/** Find the new minimised wave function.
1030 * Note, the gradient wave the function is expected in GradientTypes#ActualGradient. Then,
1031 * subsequently calls (performing thereby a conjugate gradient minimisation algorithm
1032 * of Fletchers&Reeves):
1033 * - CalculateCGGradient()\n
1034 * find direction of steepest descent: GradientTypes#ActualGradient
1035 * - CalculatePreConGrad()\n
1036 * precondition steepest descent direction: GradientTypes#PreConGradient
1037 * - CalculateConDir()\n
1038 * calculate conjugate direction: GradientTypes#ConDirGradient and GradientTypes#OldConDirGradient
1039 * - CalculateLineSearch()\n
1040 * perform the minimisation in this one-dimensional subspace: LevelPsi::LocalPsi[RunStruct::ActualLocalPsiNo]
1041 * \note CalculateGradient() is not re-done.
1042 * \param *P Problem at hand
1043 * \param x if NULL do approximative line search by second order taylor expansion, otherwise use given paramter between -M_PI...M_PI
1044 */
1045void CalculateNewWave(struct Problem *P, double *x)
1046{
1047 struct Lattice *Lat = &P->Lat;
1048 struct RunStruct *R = &P->R;
1049 struct LatticeLevel *LevS = R->LevS;
1050 const int i = R->ActualLocalPsiNo;
1051 const int ConDirStep = R->PsiStep;
1052 if (((!R->ScanPotential) || (R->ScanPotential < R->ScanAtStep)) && (Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent != 1)) { // only evaluate CG once
1053 //if (isnan(P->Grad.GradientArray[ActualGradient][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): ActualGradient_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
1054 //fprintf(stderr,"(%i) Evaluating ConDirGradient from SD direction...\n", P->Par.me);
1055 CalculateCGGradient(P, LevS->LPsi->LocalPsi[i],
1056 P->Grad.GradientArray[ActualGradient],
1057 Lat->Psi.AddData[i].Lambda);
1058 //fprintf(stderr,"CalculateCGGradient: %lg\n", GradSP(P, LevS, P->Grad.GradientArray[ActualGradient], P->Grad.GradientArray[ActualGradient]));
1059 CalculatePreConGrad(P, P->Grad.GradientArray[ActualGradient],
1060 P->Grad.GradientArray[PreConGradient],
1061 Lat->Psi.AddData[i].T);
1062 //if (isnan(P->Grad.GradientArray[PreConGradient][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): PreConGradient_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
1063 //fprintf(stderr,"CalculatePreConGrad: %lg\n", GradSP(P, LevS, P->Grad.GradientArray[PreConGradient], P->Grad.GradientArray[PreConGradient]));
1064 CalculateConDir(P, ConDirStep, &Lat->Psi.AddData[i].Gamma,
1065 P->Grad.GradientArray[ActualGradient],
1066 P->Grad.GradientArray[OldActualGradient],
1067 P->Grad.GradientArray[PreConGradient],
1068 P->Grad.GradientArray[ConDirGradient],
1069 P->Grad.GradientArray[OldConDirGradient]);
1070 //if (isnan(P->Grad.GradientArray[ConDirGradient][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): ConDirGradient_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
1071 //fprintf(stderr,"CalculateConDir: %lg\n", GradSP(P, LevS, P->Grad.GradientArray[ConDirGradient], P->Grad.GradientArray[ConDirGradient]));
1072 } //else
1073 //TestForOrth(P,R->LevS,P->Grad.GradientArray[GraSchGradient]);
1074 //fprintf(stderr,"(%i) DoBrent[%i] = %i\n", P->Par.me, R->ActualLocalPsiNo, Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].DoBrent);
1075 if (R->ScanPotential) R->ScanPotential++; // increment by one
1076 CalculateLineSearch(P,
1077 LevS->LPsi->LocalPsi[i],
1078 P->Grad.GradientArray[ConDirGradient],
1079 P->Grad.GradientArray[HcGradient],
1080 Lat->Psi.LocalPsiStatus[i].PsiFactor, x);
1081 //if (isnan(LevS->LPsi->LocalPsi[i][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateNewWave(): {\\wildetilde \\varphi}_%i[%i] = NaN!\n", P->Par.me, i, 0); Error(SomeError, "NaN-Fehler!"); }
1082 SetGramSchExtraPsi(P, &Lat->Psi, NotUsedToOrtho);
1083 //TestForOrth(P,R->LevS,P->Grad.GradientArray[GraSchGradient]);
1084}
1085
1086
1087/** Calculates the parameter alpha for UpdateWaveAfterIonMove()
1088 * Takes ratios of differences of old and older positions with old and current,
1089 * respectively with old and older. Alpha is then the ratio between the two:
1090 * \f[
1091 * \alpha = \frac{v_x[old] \cdot v_x[now] + v_y[old] \cdot v_y[now] + v_z[old] \cdot v_z[now]}{v_x[old]^2 + v_y[old]^2 + v_z[old]^2}
1092 * \f]
1093 * Something like a rate of change in the ion velocity.
1094 * \param *P Problem at hand
1095 * \return alpha
1096 */
1097static double CalcAlphaForUpdateWaveAfterIonMove(struct Problem *P)
1098{
1099 struct Ions *I = &P->Ion;
1100 struct RunStruct *Run = &P->R;
1101 double alpha=0;
1102 double Sum1=0;
1103 double Sum2=0;
1104 int is,ia;
1105 double *R,*R_old,*R_old_old;
1106 Run->AlphaStep++;
1107 if (Run->AlphaStep <= 1) return(0.0);
1108 for (is=0; is < I->Max_Types; is++) {
1109 for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
1110 //if (I->I[is].IMT[ia] == MoveIon) { // Even FixedIon moves, only not by other's forces
1111 R = &I->I[is].R[NDIM*ia];
1112 R_old = &I->I[is].R_old[NDIM*ia];
1113 R_old_old = &I->I[is].R_old_old[NDIM*ia];
1114 Sum1 += (R_old[0]-R_old_old[0])*(R[0]-R_old[0]);
1115 Sum1 += (R_old[1]-R_old_old[1])*(R[1]-R_old[1]);
1116 Sum1 += (R_old[2]-R_old_old[2])*(R[2]-R_old[2]);
1117 Sum2 += (R_old[0]-R_old_old[0])*(R_old[0]-R_old_old[0]);
1118 Sum2 += (R_old[1]-R_old_old[1])*(R_old[1]-R_old_old[1]);
1119 Sum2 += (R_old[2]-R_old_old[2])*(R_old[2]-R_old_old[2]);
1120 //}
1121 }
1122 }
1123 if (fabs(Sum2) > MYEPSILON) {
1124 alpha = Sum1/Sum2;
1125 } else {
1126 fprintf(stderr,"(%i) Sum2 <= MYEPSILON\n", P->Par.me);
1127 alpha = 0.0;
1128 }
1129 //if (isnan(alpha)) { fprintf(stderr,"(%i) CalcAlphaForUpdateWaveAfterIonMove: alpha is NaN!\n", P->Par.me); Error(SomeError, "NaN-Fehler!");}
1130 return (alpha);
1131}
1132
1133/** Updates Wave functions after UpdateIonsR.
1134 * LocalPsi coefficients are shifted by the difference between LevelPsi::LocalPsi and LevelPsi::OldLocalPsi
1135 * multiplied by alpha, new value stored in old one.
1136 * Afterwards all local wave functions are set to NotOrthogonal and then orthonormalized via
1137 * GramSch().
1138 *
1139 * \param *P Problem at hand
1140 * \sa CalcAlphaForUpdateWaveAfterIonMove() - alpha is calculated here
1141 */
1142double UpdateWaveAfterIonMove(struct Problem *P)
1143#if 0
1144static double CalculateDeltaII(double E0, double dEdt0, double E1, double deltastep, double *E, double *dE, double *ddE, double *dcos, double *dsin)
1145{
1146 double delta, dcos2, dsin2, A, B, Eavg;
1147 if (fabs(1-cos(2*deltastep) < MYEPSILON) fprintf(stderr,"CalculateDeltaII: (1-cos(2*deltastep) = %lg",(1-cos(2*deltastep));
1148 A = (E0-E1+0.5*dEdt0*sin(2*deltastep))/(1-cos(2*deltastep));
1149 B = 0.5*dEdt0;
1150 Eavg = E0-A;
1151 if (fabs(A) < MYEPSILON) fprintf(stderr,"CalculateDeltaII: A = %lg",A);
1152 delta = 0.5*atan(B/A);
1153 *dcos = cos(delta);
1154 *dsin = sin(delta);
1155 dcos2 = cos(2.*delta);
1156 dsin2 = sin(2.*delta);
1157 *E = Eavg + A*dcos2 + B*dsin2;
1158 *dE = -2.*A*dsin2+2.*B*dcos2;
1159 *ddE = -4.*A*dcos2-4.*B*dsin2;
1160 if (*ddE < 0) {
1161 if (delta < 0) {
1162 delta += PI/2.;
1163 } else {
1164 delta -= PI/2.;
1165 }
1166 *dcos = cos(delta);
1167 *dsin = sin(delta);
1168 dcos2 = cos(2.*delta);
1169 dsin2 = sin(2.*delta);
1170 *E = Eavg + A*dcos2 + B*dsin2;
1171 *dE = -2.*A*dsin2+2.*B*dcos2;
1172 *ddE = -4.*A*dcos2-4.*B*dsin2;
1173 }
1174 return(delta);
1175}
1176#endif
1177{
1178 struct RunStruct *R = &P->R;
1179 struct LatticeLevel *LevS = R->LevS;
1180 struct Lattice *Lat = &P->Lat;
1181 struct Psis *Psi = &Lat->Psi;
1182 double alpha=CalcAlphaForUpdateWaveAfterIonMove(P);
1183 int i,g,p,ResetNo=0;
1184 fftw_complex *source, *source_old;
1185 struct OnePsiElement *OnePsi = NULL;
1186 //fprintf(stderr, "%lg\n", alpha);
1187 if (P->R.AlphaStep > 1 && alpha != 0.0) {
1188 for (i=Psi->TypeStartIndex[Occupied]; i < Psi->TypeStartIndex[Perturbed_P0]; i++) {
1189 source = LevS->LPsi->LocalPsi[i];
1190 source_old = LevS->LPsi->OldLocalPsi[i];
1191 for (g=0; g < LevS->MaxG; g++) {
1192 c_re(source[g]) = c_re(source[g])+alpha*(c_re(source[g])-c_re(source_old[g]));
1193 c_im(source[g]) = c_im(source[g])+alpha*(c_im(source[g])-c_im(source_old[g]));
1194 c_re(source_old[g])=c_re(source[g]);
1195 c_im(source_old[g])=c_im(source[g]);
1196 }
1197 }
1198 } else {
1199 for (i=Psi->TypeStartIndex[Occupied]; i < Psi->TypeStartIndex[Perturbed_P0]; i++) {
1200 source = LevS->LPsi->LocalPsi[i];
1201 source_old = LevS->LPsi->OldLocalPsi[i];
1202 for (g=0; g < LevS->MaxG; g++) {
1203 c_re(source_old[g])=c_re(source[g]);
1204 c_im(source_old[g])=c_im(source[g]);
1205 }
1206 }
1207 }
1208 //fprintf(stderr, "(%i) alpha = %lg\n", P->Par.me, alpha);
1209 if (alpha != 0.0) {
1210 for (p=0; p< Psi->LocalNo; p++) {
1211 Psi->LocalPsiStatus[p].PsiGramSchStatus = (Psi->LocalPsiStatus[p].PsiType == Occupied) ? (int)NotOrthogonal : (int)NotUsedToOrtho;
1212 //if (R->CurrentMin > UnOccupied)
1213 //fprintf(stderr,"(%i) Setting L-Status of %i to %i\n",P->Par.me,p, Psi->LocalPsiStatus[p].PsiGramSchStatus);
1214 }
1215 ResetNo=0;
1216 for (i=0; i < P->Par.Max_me_comm_ST_PsiT; i++) {
1217 for (p=ResetNo; p < ResetNo+Psi->AllLocalNo[i]-1; p++) {
1218 Psi->AllPsiStatus[p].PsiGramSchStatus = (Psi->AllPsiStatus[p].PsiType == Occupied) ? (int)NotOrthogonal : (int)NotUsedToOrtho;
1219 //if (R->CurrentMin > UnOccupied)
1220 //fprintf(stderr,"(%i) Setting A-Status of %i to %i\n",P->Par.me,p, Psi->AllPsiStatus[p].PsiGramSchStatus);
1221 }
1222 ResetNo += Psi->AllLocalNo[i];
1223 OnePsi = &Psi->AllPsiStatus[ResetNo-1]; // extra wave function is set to NotUsed
1224 //if (R->CurrentMin > UnOccupied) fprintf(stderr,"Setting A-Status of %i to %i\n",ResetNo-1, NotUsedToOrtho);
1225 OnePsi->PsiGramSchStatus = (int)NotUsedToOrtho;
1226 }
1227 SpeedMeasure(P, InitGramSchTime, StartTimeDo);
1228 //fprintf(stderr,"(%i) UpdateWaveAfterIonMove: ResetGramSch() for %i orbitals\n",P->Par.me, p);
1229 GramSch(P, LevS, Psi, Orthonormalize);
1230 SpeedMeasure(P, InitGramSchTime, StartTimeDo);
1231 }
1232 return(alpha);
1233}
1234
Note: See TracBrowser for help on using the repository browser.