source: pcp/src/grad.c@ b76b97

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

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