source: pcp/src/energy.c@ c510a7

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

added some comment-out fprintf lines to CalculateGapEnergy()
  • Property mode set to 100644
File size: 36.0 KB
Line 
1/** \file energy.c
2 * Energy calculation, especially kinetic and electronic.
3 *
4 * Contains routines for output of calculated values EnergyOutput(),
5 * initialization of the Energy structure InitEnergyArray() or inital calculations
6 * InitPsiEnergyCalculation(),
7 * updating UpdateEnergyArray() and UpdatePsiEnergyCalculation() and of course
8 * calculation of the various parts: Density CalculateDensityEnergy(), ionic
9 * CalculateIonsEnergy(), electronic CalculatePsiEnergy(), kinetic with RiemannTensor
10 * CalculateKineticEnergyUseRT() (calling EnergyGNSP()) and without
11 * CalculateKineticEnergyNoRT() (using EnergyLapNSP()). CalculateGapEnergy() evaluates
12 * would-be addition from unoccupied states to total energy.
13 * Smaller functions such as EnergyAllReduce() gather intermediate results from all
14 * processes and sum up to total value.
15 * \note Some others are stored in excor.c (exchange energies) and pseudo.c (density
16 * energies)
17 *
18 *
19 Project: ParallelCarParrinello
20 \author Jan Hamaekers
21 \date 2000
22
23 File: energy.c
24 $Id: energy.c,v 1.68 2007/02/09 09:13:48 foo Exp $
25*/
26
27#include <stdlib.h>
28#include <stdio.h>
29#include <stddef.h>
30#include <math.h>
31#include <string.h>
32
33#include "data.h"
34#include "errors.h"
35#include "helpers.h"
36#include "energy.h"
37#include "riemann.h"
38#include "excor.h"
39#include "myfft.h"
40#include "mymath.h"
41//#include "output.h"
42#include "run.h"
43
44/** Calculates kinetic energy for a given wave function in reciprocal base.
45 * \param *P Problem at hand
46 * \param *Lev LatticeLevel structure
47 * \param *LPsiDatA complex wave function coefficients \f$c_{i,G}\f$ (reciprocal base)
48 * \param *T kinetic eigenvalue, herein result is additionally stored
49 * \param Factor Additional factor for return value (such as occupation number \f$f_i\f$ )
50 * \return \f$E_{kin} = \frac{1}{2} f_i \sum_G |G|^2 |c_{i,G}|^2\f$
51 * \sa used in CalculateKineticEnergyNoRT()
52 */
53static double EnergyGNSP(const struct Problem *P, const struct LatticeLevel *Lev, fftw_complex *LPsiDatA, double *T, const double Factor) {
54 double LocalSP=0.0,PsiSP;
55 int i,s = 0;
56 struct OneGData *G = Lev->GArray;
57 if (Lev->GArray[0].GSq == 0.0) { // only for continuity (see other functions)
58 LocalSP += G[0].GSq*LPsiDatA[0].re*LPsiDatA[0].re;
59 s++;
60 }
61 /// Performs complex product for each reciprocal grid vector times twice this vector's norm.
62 for (i=s; i < Lev->MaxG; i++) {
63 LocalSP += G[i].GSq*2.*(LPsiDatA[i].re*LPsiDatA[i].re+LPsiDatA[i].im*LPsiDatA[i].im); // factor 2 due to gamma point symmetry
64 }
65 /// Gathers partial results by summation from all other coefficients sharing processes.
66 MPI_Allreduce ( &LocalSP, &PsiSP, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
67 /// Multiplies by \a Factor and returns result.
68 *T = Factor*PsiSP;
69 return(Factor*PsiSP);
70}
71
72/** Calculates kinetic energy for a given wave function in reciprocal base in case of Riemann tensor use.
73 * \param *P Problem at hand
74 * \param *RT RiemannTensor structure
75 * \param *src replaces reciprocal grid vectors
76 * \param *Lap replaces complex Psi wave function coefficients
77 * \param *T kinetic eigenvalue, herein result is additionally stored
78 * \param Factor Additional factor for return value
79 * \sa used in CalculateKineticEnergyUseRT(), compare with EnergyGNSP()
80 */
81static double EnergyLapNSP(const struct Problem *P, const struct RiemannTensor *RT, fftw_complex *src, fftw_complex *Lap, double *T, const double Factor) {
82 double LocalSP=0.0,PsiSP;
83 int i;
84 int s = 0;
85 struct LatticeLevel *Lev = RT->LevS;
86 if (Lev->GArray[0].GSq == 0.0) {
87 LocalSP += src[0].re*Lap[0].re+src[0].im*Lap[0].im;
88 s++;
89 }
90 for (i=s; i < Lev->MaxG; i++) {
91 LocalSP += 2.*(src[i].re*Lap[i].re+src[i].im*Lap[i].im);
92 }
93 MPI_Allreduce ( &LocalSP, &PsiSP, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
94 *T = Factor*PsiSP;
95 return(Factor*PsiSP);
96}
97
98/** Calculcates Kinetic energy without Riemann tensor.
99 * \param *P Problem at hand
100 * \param p local psi wave function number
101 * \sa EnergyGNSP()
102 */
103static void CalculateKineticEnergyNoRT(struct Problem *P, const int p)
104{
105 struct Lattice *Lat = &P->Lat;
106 struct Energy *E = Lat->E;
107 struct RunStruct *R = &P->R;
108 struct LatticeLevel *LevS = R->LevS;
109 struct LevelPsi *LPsi = LevS->LPsi;
110 struct Psis *Psi = &Lat->Psi;
111 const double PsiFactor = Psi->LocalPsiStatus[p].PsiFactor;
112 fftw_complex *psi;
113 psi = LPsi->LocalPsi[p];
114 E->PsiEnergy[KineticEnergy][p] = EnergyGNSP(P, LevS, psi, &Psi->AddData[p].T, 0.5*PsiFactor);
115}
116
117/** Calculcates Kinetic energy with Riemann tensor.
118 * CalculateRiemannLaplaceS() is called and Lattice->RT.RTLaplaceS used instead
119 * of Lattice::Psi.
120 * \param *P Problem at hand
121 * \param p local psi wave function number
122 * \sa EnergyLapNSP()
123 */
124static void CalculateKineticEnergyUseRT(struct Problem *P, const int p)
125{
126 struct Lattice *Lat = &P->Lat;
127 struct Energy *E = Lat->E;
128 struct RunStruct *R = &P->R;
129 struct LatticeLevel *LevS = R->LevS;
130 struct LevelPsi *LPsi = LevS->LPsi;
131 struct Psis *Psi = &Lat->Psi;
132 const double PsiFactor = Psi->LocalPsiStatus[p].PsiFactor;
133 fftw_complex *psi;
134 psi = LPsi->LocalPsi[p];
135 CalculateRiemannLaplaceS(Lat, &Lat->RT, psi, Lat->RT.RTLaplaceS);
136 E->PsiEnergy[KineticEnergy][p] = EnergyLapNSP(P, &Lat->RT, psi, Lat->RT.RTLaplaceS, &Psi->AddData[p].T, 0.5*PsiFactor);
137 //fprintf(stderr,"CalculateKineticEnergyUseRT: PsiEnergy[KineticEnergy] = %lg\n",E->PsiEnergy[KineticEnergy][p]);
138}
139
140/** Sums up calculated energies from all processes, calculating Energy::TotalEnergy.
141 * Via MPI_Allreduce is Energy::AllLocalPsiEnergy summed up and the sum
142 * redistributed back to all processes for the various spin cases. In
143 * case of SpinType#SpinUp or SpinType#SpinDown the "other" energy (e.g. down
144 * for an up process) is via MPI_SendRecv communicated, and lateron
145 * Energy#AllTotalPsiEnergy as the sum of two states generated in each process.<br>
146 * Finally, the Energy::TotalEnergy[0] is the sum of Energy::AllTotalPsiEnergy[],
147 * Energy::AllTotalDensityEnergy[], - Energy::AllTotalIonsEnergy[GaussSelfEnergy]
148 * and Energy::AllTotalIonsEnergy[EwaldEnergy].
149 * Again, as in UpdateEnergyArray(), this is split up completely for PsiTypeTag#UnOccupied
150 * and PsiTypeTag#Occupied.
151 * \param *P Problem at hand
152 */
153void EnergyAllReduce(struct Problem *P)
154{
155 struct Psis *Psi = &P->Lat.Psi;
156 struct RunStruct *R = &P->R;
157 struct Lattice *Lat = &P->Lat;
158 struct Energy *E = P->Lat.E;
159 struct OnePsiElement *OnePsiB;
160 MPI_Status status;
161 int i,j,m;
162 double lambda, overlap, lambda2, overlap2; //, part1, part2, tmp;
163 MPI_Allreduce(E->AllLocalDensityEnergy, E->AllTotalDensityEnergy, MAXALLDENSITYENERGY, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi );
164 switch (Psi->PsiST) {
165 case SpinDouble:
166 MPI_Allreduce(E->AllLocalPsiEnergy, E->AllTotalPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
167 break;
168 case SpinUp:
169 MPI_Allreduce(E->AllLocalPsiEnergy, E->AllUpPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
170 MPI_Sendrecv(E->AllUpPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag1,
171 E->AllDownPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag2,
172 P->Par.comm_STInter, &status );
173 for (i=0; i< MAXALLPSIENERGY; i++)
174 E->AllTotalPsiEnergy[i] = E->AllUpPsiEnergy[i] + E->AllDownPsiEnergy[i];
175 break;
176 case SpinDown:
177 MPI_Allreduce(E->AllLocalPsiEnergy, E->AllDownPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT);
178 MPI_Sendrecv(E->AllDownPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag2,
179 E->AllUpPsiEnergy, MAXALLPSIENERGY, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag1,
180 P->Par.comm_STInter, &status );
181 for (i=0; i< MAXALLPSIENERGY; i++)
182 E->AllTotalPsiEnergy[i] = E->AllUpPsiEnergy[i] + E->AllDownPsiEnergy[i];
183 break;
184 }
185 switch (R->CurrentMin) {
186 case Occupied:
187 E->TotalEnergy[0] = E->AllTotalPsiEnergy[KineticEnergy];
188 E->TotalEnergy[0] += E->AllTotalPsiEnergy[NonLocalEnergy]
189 ;
190 E->TotalEnergy[0] += E->AllTotalDensityEnergy[CorrelationEnergy]+E->AllTotalDensityEnergy[ExchangeEnergy];
191 E->TotalEnergy[0] += E->AllTotalDensityEnergy[PseudoEnergy];
192 E->TotalEnergy[0] += E->AllTotalDensityEnergy[HartreePotentialEnergy];
193 E->TotalEnergy[0] -= E->AllTotalIonsEnergy[GaussSelfEnergy];
194 E->TotalEnergy[0] += E->AllTotalIonsEnergy[EwaldEnergy];
195// fprintf(stderr,"(%i) Tot %1.5f ---\t Diff %1.5f\t ---\t Kin %1.5f\tNLErg %1.5f ---\t Corr %1.5f\tXC %1.5f\t Pseudo %1.5f\tHP %1.5f\n", P->Par.me,
196// E->TotalEnergy[0],(E->TotalEnergy[1]-E->TotalEnergy[0])/fabs(E->TotalEnergy[0]),
197// E->AllTotalPsiEnergy[KineticEnergy],E->AllTotalPsiEnergy[NonLocalEnergy],
198// E->AllTotalDensityEnergy[CorrelationEnergy],E->AllTotalDensityEnergy[ExchangeEnergy],
199// E->AllTotalDensityEnergy[PseudoEnergy],E->AllTotalDensityEnergy[HartreePotentialEnergy]);
200 break;
201 case UnOccupied:
202 E->TotalEnergy[0] = E->AllTotalPsiEnergy[KineticEnergy]+E->AllTotalPsiEnergy[NonLocalEnergy];
203 E->TotalEnergy[0] += E->AllTotalDensityEnergy[CorrelationEnergy]+E->AllTotalDensityEnergy[ExchangeEnergy];
204 E->TotalEnergy[0] += E->AllTotalDensityEnergy[PseudoEnergy]+E->AllTotalDensityEnergy[HartreePotentialEnergy];
205 break;
206 case Perturbed_P0:
207 case Perturbed_P1:
208 case Perturbed_P2:
209 case Perturbed_RxP0:
210 case Perturbed_RxP1:
211 case Perturbed_RxP2:
212 E->TotalEnergy[0] = E->AllTotalPsiEnergy[Perturbed1_0Energy]+E->AllTotalPsiEnergy[Perturbed0_1Energy];
213 E->parts[0] = E->AllTotalPsiEnergy[Perturbed1_0Energy]+E->AllTotalPsiEnergy[Perturbed0_1Energy];
214 //fprintf(stderr,"(%i) summed single: AllTotalPsiEnergy[total] = %lg + %lg = %lg\n",P->Par.me,E->AllTotalPsiEnergy[Perturbed1_0Energy],E->AllTotalPsiEnergy[Perturbed0_1Energy], E->TotalEnergy[0]);
215 //E->TotalEnergy[0] = 0.;
216 m = -1;
217 E->parts[1] = E->parts[2] = 0.;
218 for (i=0; i < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; i++) { // go through all wave functions: \sum_k
219 OnePsiB = &Psi->AllPsiStatus[i]; // grab OnePsiB
220 if (OnePsiB->PsiType == R->CurrentMin) { // drop all but occupied types
221 m++; // increase m if it is occupied wave function
222 lambda = overlap = 0.;
223 if (OnePsiB->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) { // local?
224 //fprintf(stderr,"(%i) Lambda[%i] = %lg\n",P->Par.me, m,Lat->Psi.AddData[i].Lambda);
225 lambda = Lat->Psi.AddData[OnePsiB->MyLocalNo].Lambda * Lat->Psi.LocalPsiStatus[OnePsiB->MyLocalNo].PsiFactor;
226 for (j=0;j<Psi->NoOfPsis;j++) { // over all Psis: \sum_l
227 //lambda -= Lat->Psi.lambda[j][i - Psi->TypeStartIndex[R->CurrentMin]]*GradSP(P,R->LevS,R->LevS->LPsi->LocalPsi[j + Psi->TypeStartIndex[R->CurrentMin]],R->LevS->LPsi->LocalPsi[i]);
228 overlap -= Lat->Psi.lambda[j][m]*Lat->Psi.Overlap[j][m];
229 }
230 // send results to other processes
231 }
232 //fprintf(stderr,"(%i) before MPI: \tlambda[%i] = %lg\t overlap[%i] = %lg\n",P->Par.me,m,lambda,m,overlap);
233 MPI_Allreduce( &lambda, &lambda2, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT); // lambda is just local
234 MPI_Allreduce( &overlap, &overlap2, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_PsiT); // we just summed up over local m's!
235 //fprintf(stderr,"(%i) after MPI: \t lambda[%i] = %lg\t overlap[%i] = %lg\n",P->Par.me, OnePsiB->MyGlobalNo,lambda2, m, overlap2);
236 E->TotalEnergy[0] += lambda2;
237 E->TotalEnergy[0] += overlap2;
238 E->parts[1] += lambda2;
239 E->parts[2] += overlap2;
240 }
241 }
242
243 switch(Psi->PsiST) {
244 default:
245 lambda = overlap = 0.;
246 break;
247 case SpinUp:
248 MPI_Sendrecv(&E->parts[1], 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag1,
249 &lambda, 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag2, P->Par.comm_STInter, &status );
250 E->TotalEnergy[0] += lambda;
251 MPI_Sendrecv(&E->parts[2],1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag3,
252 &overlap, 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag4, P->Par.comm_STInter, &status );
253 E->TotalEnergy[0] += overlap;
254 break;
255 case SpinDown:
256 MPI_Sendrecv(&E->parts[1], 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag2,
257 &lambda, 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag1, P->Par.comm_STInter, &status );
258 E->TotalEnergy[0] += lambda;
259 MPI_Sendrecv(&E->parts[2],1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag4,
260 &overlap, 1, MPI_DOUBLE, P->Par.me_comm_ST, EnergyTag3, P->Par.comm_STInter, &status );
261 E->TotalEnergy[0] += overlap;
262 break;
263 }
264 //fprintf(stderr,"(%i) summed single: lambda[total] = %lg\t overlap[total] = %lg\n",P->Par.me,E->parts[1]+lambda, E->parts[2]+overlap);
265 break;
266 default:
267 break;
268 }
269}
270
271/** Initializes Energy Array.
272 * Arrays are malloc'd and set to zero, then InitExchangeCorrelationEnergy()
273 * called.
274 * \param *P Problem at hand
275 */
276void InitEnergyArray(struct Problem *P)
277{
278 struct Lattice *Lat = &P->Lat;
279 struct Energy *E = Lat->E;
280 struct Psis *Psi = &Lat->Psi;
281 int i, type, oldtype = P->R.CurrentMin;
282 for (type=Occupied;type<Extra;type++) {
283 SetCurrentMinState(P, type);
284 for (i=0; i<MAXALLPSIENERGY; i++) {
285 Lat->Energy[type].PsiEnergy[i] = (double *)
286 Malloc(sizeof(double)*Psi->LocalNo,"InitEnergyCalculations: PsiKineticEnergy");
287 SetArrayToDouble0(Lat->Energy[type].PsiEnergy[i], Psi->LocalNo);
288 }
289 SetArrayToDouble0(E->AllLocalDensityEnergy, MAXALLDENSITYENERGY);
290 SetArrayToDouble0(E->AllTotalDensityEnergy, MAXALLDENSITYENERGY);
291 SetArrayToDouble0(E->AllTotalIonsEnergy, MAXALLIONSENERGY);
292 SetArrayToDouble0(E->TotalEnergy, MAXOLD);
293 }
294 InitExchangeCorrelationEnergy(P, &P->ExCo);
295 SetCurrentMinState(P, oldtype);
296}
297
298/** Shifts stored total energy values up to MAXOLD.
299 * Energy::PsiEnergy, Energy::AllLocalDensityEnergy, Energy::AllTotalDensityEnergy
300 * are all set to zero. Energy::TotalEnergy is shifted by one to make
301 * place for next value. This is done either for occupied or for the gap energies,
302 * stated by \a IncType.
303 * \param *P Problem at hand
304 */
305void UpdateEnergyArray(struct Problem *P)
306{
307 struct Lattice *Lat = &P->Lat;
308 struct RunStruct *R = &P->R;
309 struct Energy *E = Lat->E;
310 int i;
311 for (i=0; i<MAXALLPSIENERGY; i++)
312 E->PsiEnergy[i][R->OldActualLocalPsiNo] = 0.0;
313 SetArrayToDouble0(E->AllLocalDensityEnergy, MAXALLDENSITYENERGY);
314 SetArrayToDouble0(E->AllTotalDensityEnergy, MAXALLDENSITYENERGY);
315 for (i=MAXOLD-1; i > 0; i--)
316 Lat->Energy[R->CurrentMin].TotalEnergy[i] = Lat->Energy[R->CurrentMin].TotalEnergy[i-1];
317}
318
319/** Calculates ion energy.
320 * CalculateGaussSelfEnergyNoRT() is called.
321 * \param *P Problem at hand
322 * \warning CalculateIonsUseRT not implemented
323 */
324void CalculateIonsEnergy(struct Problem *P)
325{
326 struct Lattice *Lat = &P->Lat;
327 switch (Lat->RT.ActualUse) {
328 case inactive:
329 case standby:
330 CalculateGaussSelfEnergyNoRT(P, &P->PP, &P->Ion);
331 break;
332 case active:
333 Error(SomeError, "CalculateIonsUseRT: Not implemented");
334 break;
335 }
336}
337
338/** Calculates density energy.
339 * Depending on RT::ActualUse CalculateXCEnergyNoRT(), CalculateSomeEnergyAndHGNoRT()
340 * and CalculateGapEnergy() are called within SpeedMeasure() brackets.
341 * Afterwards, CalculateGradientNoRT() for RunStruct::ActualLocalPsiNo.
342 * \param *P Problem at hand
343 * \param UseInitTime whether SpeedMeasure uses InitLocTime (Init) or LocTime (Update)
344 * \warning CalculateSomeEnergyAndHGUseRT not implemented
345 * \warning call EnergyAllReduce() afterwards!
346 */
347void CalculateDensityEnergy(struct Problem *P, const int UseInitTime)
348{
349 struct Lattice *Lat = &P->Lat;
350 struct RunStruct *R = &P->R;
351 struct LatticeLevel *LevS = R->LevS;
352 switch (Lat->RT.ActualUse) {
353 case inactive:
354 case standby:
355 SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StartTimeDo);
356 CalculateXCEnergyNoRT(P);
357 CalculateSomeEnergyAndHGNoRT(P);
358 SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StopTimeDo);
359 CalculateGradientNoRT(P,
360 LevS->LPsi->LocalPsi[R->ActualLocalPsiNo],
361 Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].PsiFactor,
362 &Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda,
363 P->PP.fnl[R->ActualLocalPsiNo],
364 P->Grad.GradientArray[ActualGradient],
365 P->Grad.GradientArray[OldActualGradient],
366 P->Grad.GradientArray[HcGradient]);
367 break;
368 case active:
369 CalculateXCEnergyUseRT(P);
370 /* FIXME - Hier fehlt noch was */
371 Error(SomeError, "CalculateSomeEnergyAndHGUseRT: Not implemented");
372 break;
373 }
374 /* Achtung danach noch EnergyAllReduce aufrufen */
375}
376
377/** Calculation of an additional total energy from unoccupied states.
378 * By including additional (unoccupied) orbitals into the minimisation process,
379 * which also stand orthogonal on all others, however are minimised in the field
380 * of occupied ones (without altering it in a self-consistent manner), a first
381 * approximation to the band gap energy can be obtained. The density part of this
382 * virtual energy is stored in the PsiTypeTag#UnOccupied part of Lattice#Energy.
383 * Note that the band gap is approximated via the eigenvalues of the additional
384 * orbitals. This energy here is purely fictious.
385 * \param *P Problem at hand
386 * \note No RiemannTensor use so far implemented. And mostly copy&paste from
387 * CalculateXCEnergyNoRT() and CalculateSomeEnergyAndHGNoRT().
388 * \bug Understand CoreCorrection part in \f$E_{XC}\f$
389 */
390void CalculateGapEnergy(struct Problem *P) {
391 struct Lattice *Lat = &P->Lat;
392 struct Psis *Psi = &Lat->Psi;
393 struct Energy *E = Lat->E;
394 struct PseudoPot *PP = &P->PP;
395 struct RunStruct *R = &P->R;
396 struct LatticeLevel *Lev0 = R->Lev0;
397 struct LatticeLevel *LevS = R->LevS;
398 struct Density *Dens = Lev0->Dens;
399 struct ExCor *EC = &P->ExCo;
400 struct Ions *I = &P->Ion;
401 double SumEc = 0;
402 double rs, p = 0.0, pUp = 0.0, pDown = 0.0, zeta, rsUp, rsDown, SumEx=0.0;
403 double Factor = R->XCEnergyFactor/Lev0->MaxN;
404 fftw_complex vp,rhoe,rhoe_gap;
405 fftw_complex rp,rhog;
406 double SumEHP,Fac,SumEH,SumEG,SumEPS;
407 int g,s=0,it,Index;
408 int i;
409 fftw_complex *HG = Dens->DensityCArray[HGDensity];
410 SetArrayToDouble0((double *)HG,2*Dens->TotalSize);
411
412 //if (isnan(HG[0].re)) { fprintf(stderr,"WARNGING: CalculateGapEnergy(): HG[%i]_%i = NaN!\n", 0, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
413 SpeedMeasure(P, GapTime, StartTimeDo);
414 // === Calculate exchange and correlation energy ===
415 for (i = 0; i < Dens->LocalSizeR; i++) { // for each node in radial mesh
416 // put (corecorrected) densities in p, pUp, pDown
417 p = Dens->DensityArray[TotalDensity][i];
418 //if (isnan(p)) { fprintf(stderr,"WARNGING: CalculateGapEnergy(): Dens->DensityArray[TotalUpDensity][%i]_%i = NaN!\n", i, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
419 switch (Psi->PsiST) {
420 case SpinDouble:
421 pUp = 0.5*p;
422 pDown = 0.5*p;
423 break;
424 case SpinUp:
425 case SpinDown:
426 pUp = Dens->DensityArray[TotalUpDensity][i];
427 pDown = Dens->DensityArray[TotalDownDensity][i];
428 break;
429 default:
430 ;
431 }
432 if ((p < 0) || (pUp < 0) || (pDown < 0)) {
433 /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
434 p = 0.0;
435 pUp = 0.0;
436 pDown = 0.0;
437 }
438 rs = Calcrs(EC,p);
439 rsUp = Calcrs(EC,pUp);
440 rsDown = Calcrs(EC,pDown);
441 zeta = CalcZeta(EC,pUp,pDown);
442 // Calculation with the densities and summation
443 SumEx += CalcSEXr(EC, rsUp, rsDown, pUp, pDown);
444 SumEc += CalcSECr(EC, rs, zeta, p);
445 }
446 E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc;
447 E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*SumEx;
448
449 // === Calculate electrostatic and local pseudopotential energy ===
450 SumEHP = 0.0;
451 SumEH = 0.0;
452 SumEG = 0.0;
453 SumEPS = 0.0;
454 // calculates energy of local pseudopotential
455 if (Lev0->GArray[0].GSq == 0.0) {
456 Index = Lev0->GArray[0].Index;
457 c_re(vp) = 0.0;
458 c_im(vp) = 0.0;
459 for (it = 0; it < I->Max_Types; it++) {
460 c_re(vp) += (c_re(I->I[it].SFactor[0])*PP->phi_ps_loc[it][0]);
461 c_im(vp) += (c_im(I->I[it].SFactor[0])*PP->phi_ps_loc[it][0]);
462 }
463 //if (isnan(c_re(Dens->DensityCArray[GapDensity][Index]))) { fprintf(stderr,"WARNGING: CalculateGapEnergy(): Dens->DensityArray[GapDensity][%i]_%i.re = NaN!\n", Index, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
464 SumEPS += (c_re(Dens->DensityCArray[GapDensity][Index])*c_re(vp) +
465 c_im(Dens->DensityCArray[GapDensity][Index])*c_im(vp))*R->HGcFactor;
466 s++;
467 }
468 for (g=s; g < Lev0->MaxG; g++) {
469 Index = Lev0->GArray[g].Index;
470 c_re(vp) = 0.0;
471 c_im(vp) = 0.0;
472 c_re(rp) = 0.0;
473 c_im(rp) = 0.0;
474 Fac = 2.*PI/(Lev0->GArray[g].GSq); // 4*.. -> 2*...: 1/2 here ! (due to lacking dependence of first n^{tot} (r) on gap psis)
475 for (it = 0; it < I->Max_Types; it++) {
476 c_re(vp) += (c_re(I->I[it].SFactor[g])*PP->phi_ps_loc[it][g]);
477 c_im(vp) += (c_im(I->I[it].SFactor[g])*PP->phi_ps_loc[it][g]);
478 c_re(rp) += (c_re(I->I[it].SFactor[g])*PP->FacGauss[it][g]);
479 c_im(rp) += (c_im(I->I[it].SFactor[g])*PP->FacGauss[it][g]);
480 }
481 //if (isnan(c_re(Dens->DensityCArray[TotalDensity][Index]))) { fprintf(stderr,"WARNGING: CalculateGapEnergy(): Dens->DensityArray[TotalDensity][%i]_%i.re = NaN!\n", Index, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
482 c_re(rhog) = c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_re(rp);
483 c_im(rhog) = c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_im(rp);
484 c_re(rhoe) = c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
485 c_im(rhoe) = c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
486 //if (isnan(c_re(Dens->DensityCArray[GapDensity][Index]))) { fprintf(stderr,"WARNGING: CalculateGapEnergy(): Dens->DensityArray[GapDensity][%i]_%i.re = NaN!\n", Index, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
487 c_re(rhoe_gap) = c_re(Dens->DensityCArray[GapDensity][Index])*R->HGcFactor;
488 c_im(rhoe_gap) = c_im(Dens->DensityCArray[GapDensity][Index])*R->HGcFactor;
489 //c_re(PP->VCoulombc[g]) = Fac*c_re(rhog);
490 //c_im(PP->VCoulombc[g]) = -Fac*c_im(rhog);
491 c_re(HG[Index]) = c_re(vp) + Fac * (c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_re(rp));
492 c_im(HG[Index]) = c_im(vp) + Fac * (c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_im(rp));
493 //if (isnan(c_re(HG[Index]))) { fprintf(stderr,"WARNGING: CalculateGapEnergy(): HG[%i]_%i.re = NaN!\n", Index, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
494 //if (isnan(c_im(HG[Index]))) { fprintf(stderr,"WARNGING: CalculateGapEnergy(): HG[%i]_%i.im = NaN!\n", Index, R->ActualLocalPsiNo); Error(SomeError, "NaN-Fehler!"); }
495 /*if (P->first) */
496 //SumEG += Fac*(c_re(rp)*c_re(rp)+c_im(rp)*c_im(rp)); // Coulomb energy
497 // changed: took out gaussian part! rhog -> rhoe
498 SumEHP += 2*Fac*(c_re(rhoe)*c_re(rhoe_gap)+c_im(rhoe)*c_im(rhoe_gap)); // E_ES first part
499 //SumEH += Fac*(c_re(rhoe)*c_re(rhoe)+c_im(rhoe)*c_im(rhoe));
500 SumEPS += 2.*(c_re(Dens->DensityCArray[GapDensity][Index])*c_re(vp)+
501 c_im(Dens->DensityCArray[GapDensity][Index])*c_im(vp))*R->HGcFactor;
502 }
503 //
504 for (i=0; i<Lev0->MaxDoubleG; i++) {
505 HG[Lev0->DoubleG[2*i+1]].re = HG[Lev0->DoubleG[2*i]].re;
506 HG[Lev0->DoubleG[2*i+1]].im = -HG[Lev0->DoubleG[2*i]].im;
507 }
508 /*if (P->first) */
509 //E->AllLocalDensityEnergy[GaussEnergy] = SumEG*Lat->Volume;
510 E->AllLocalDensityEnergy[PseudoEnergy] = SumEPS*Lat->Volume; // local pseudopotential
511 E->AllLocalDensityEnergy[HartreePotentialEnergy] = SumEHP*Lat->Volume; // Gauss n^Gauss and Hartree n potential
512 //E->AllLocalDensityEnergy[HartreeEnergy] = SumEH*Lat->Volume;
513 SpeedMeasure(P, GapTime, StopTimeDo);
514
515 // === Calculate Gradient ===
516 //if (isnan(P->R.LevS->LPsi->LocalPsi[P->R.ActualLocalPsiNo][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateGapEnergy(): ActualLocalPsi_%i[%i] = NaN!\n", P->Par.me, P->R.ActualLocalPsiNo, 0); Error(SomeError, "NaN-Fehler!"); }
517 CalculateGradientNoRT(P,
518 LevS->LPsi->LocalPsi[R->ActualLocalPsiNo],
519 Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].PsiFactor,
520 &Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda,
521 P->PP.fnl[R->ActualLocalPsiNo],
522 P->Grad.GradientArray[ActualGradient],
523 P->Grad.GradientArray[OldActualGradient],
524 P->Grad.GradientArray[HcGradient]);
525 //if (isnan(P->Grad.GradientArray[ActualGradient][0].re)) { fprintf(stderr,"(%i) WARNING in CalculateGapEnergy(): ActualGradient_%i[%i] = NaN!\n", P->Par.me, P->R.ActualLocalPsiNo, 0); Error(SomeError, "NaN-Fehler!"); }
526}
527
528/** Calculates the energy of a wave function psi.
529 * Calculates kinetc and non local energies, taking care of possible
530 * Riemann tensor usage.
531 * \param *P Problem at hand
532 * \param p number of wave function Psi
533 * \param UseInitTime Whether speed measuring uses InitLocTime (Init) or LocTime (Update)
534 * \sa EnergyGNSP(), EnergyLapNSP()
535 * \warning call EnergyAllReduce() afterwards!
536 */
537void CalculatePsiEnergy(struct Problem *P, const int p, const int UseInitTime)
538{
539 struct Lattice *Lat = &(P->Lat);
540 switch (Lat->RT.ActualUse) {
541 case inactive:
542 case standby:
543 SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StartTimeDo);
544 CalculateKineticEnergyNoRT(P, p);
545 SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StopTimeDo);
546 SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StartTimeDo);
547 CalculateNonLocalEnergyNoRT(P, p);
548 SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StopTimeDo);
549 break;
550 case active:
551 SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StartTimeDo);
552 CalculateKineticEnergyUseRT(P, p);
553 SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StopTimeDo);
554 SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StartTimeDo);
555 CalculateNonLocalEnergyUseRT(P, p);
556 SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StopTimeDo);
557 break;
558 }
559 /* Achtung danach noch EnergyAllReduce aufrufen */
560}
561
562/** Initialization of wave function energy calculation.
563 * For every wave function CalculatePsiEnergy() is called,
564 * Energy::AllLocalPsiEnergy is set to sum of these.
565 * \param *P Problem at hand
566 * \param type minimisation type PsiTypeTag to calculate psi energy for
567 * \warning call EnergyAllReduce() afterwards!
568 */
569void InitPsiEnergyCalculation(struct Problem *P, enum PsiTypeTag type)
570{
571 struct Lattice *Lat = &(P->Lat);
572 int p,i, oldtype = P->R.CurrentMin;
573 SetCurrentMinState(P,type);
574 for (p=0; p < Lat->Psi.LocalNo; p++) {
575 if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
576 CalculatePsiEnergy(P, p, 1);
577 }
578 for (i=0; i< MAXALLPSIENERGY; i++) {
579 Lat->E->AllLocalPsiEnergy[i] = 0.0;
580 for (p=0; p < Lat->Psi.LocalNo; p++)
581 if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
582 Lat->E->AllLocalPsiEnergy[i] += Lat->E->PsiEnergy[i][p];
583 }
584 SetCurrentMinState(P,oldtype);
585 /* Achtung danach noch EnergyAllReduce aufrufen */
586}
587
588/** Updating wave function energy.
589 * CalculatePsiEnergy() is called for RunStruct::OldActualLocalPsiNo,
590 * Energy::AllLocalPsiEnergy is set to sum of all Psi energies.
591 * \param *P Problem at hand
592 * \warning call EnergyAllReduce() afterwards!
593 */
594void UpdatePsiEnergyCalculation(struct Problem *P)
595{
596 struct Lattice *Lat = &(P->Lat);
597 struct RunStruct *R = &P->R;
598 int p = R->OldActualLocalPsiNo, i;
599 if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
600 CalculatePsiEnergy(P, p, 0);
601 for (i=0; i< MAXALLPSIENERGY; i++) {
602 Lat->E->AllLocalPsiEnergy[i] = 0.0;
603 for (p=0; p < P->Lat.Psi.LocalNo; p++)
604 if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
605 Lat->E->AllLocalPsiEnergy[i] += Lat->E->PsiEnergy[i][p];
606 }
607 /* Achtung danach noch EnergyAllReduce aufrufen */
608}
609
610/** Output of energies.
611 * Prints either Time, PsiEnergy(Up/Down), DensityEnergy, IonsEnergy, KineticEnergy,
612 * TotalEnergy and Temperature to screen or<br>
613 * Time, TotalEnergy, NonLocalEnergy, CorrelationEnergy, ExchangeEnergy,
614 * PseudoEnergy, HartreePotentialEnergy, GaussSelfEnergy, EwaldEnergy,
615 * KineticEnergy and TotalEnergy+KineticEnergy to the energy file
616 * \param *P Problem at hand
617 * \param doout 1 - print stuff to stderr, 0 - put compact form into energy file
618 */
619void EnergyOutput(struct Problem *P, int doout)
620{
621 struct Lattice *Lat = &P->Lat;
622 struct RunStruct *R = &P->R;
623 struct LatticeLevel *Lev0 = R->Lev0;
624 struct LatticeLevel *LevS = R->LevS;
625 struct Ions *I = &P->Ion;
626 struct Psis *Psi = &Lat->Psi;
627 struct FileData *F = &P->Files;
628 FILE *eout = NULL;
629 FILE *convergence = NULL;
630 int i, it;
631 if (P->Call.out[LeaderOut] && (P->Par.me == 0) && doout) {
632 fprintf(stderr, "Time T: :\t%e\n",P->R.t);
633 fprintf(stderr, "PsiEnergy :");
634 for (i=0; i<MAXALLPSIENERGY; i++) {
635 fprintf(stderr," %e",Lat->Energy[Occupied].AllTotalPsiEnergy[i]);
636 }
637 fprintf(stderr,"\n");
638 if (Psi->PsiST != SpinDouble) {
639 fprintf(stderr, "PsiEnergyUp :");
640 for (i=0; i<MAXALLPSIENERGY; i++) {
641 fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllUpPsiEnergy[i]);
642 }
643 fprintf(stderr,"\n");
644 fprintf(stderr, "PsiEnergyDown:");
645 for (i=0; i<MAXALLPSIENERGY; i++) {
646 fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllDownPsiEnergy[i]);
647 }
648 fprintf(stderr,"\n");
649 }
650 fprintf(stderr, "DensityEnergy:");
651 for (i=0; i<MAXALLDENSITYENERGY; i++) {
652 fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllTotalDensityEnergy[i]);
653 }
654 fprintf(stderr,"\n");
655 fprintf(stderr, "IonsEnergy :");
656 for (i=0; i<MAXALLIONSENERGY; i++) {
657 fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllTotalIonsEnergy[i]);
658 }
659 fprintf(stderr,"\n");
660 fprintf(stderr, "TotalEnergy :\t%e\n",Lat->Energy[Occupied].TotalEnergy[0]);
661 if (R->DoPerturbation)
662 fprintf(stderr, "PerturbedEnergy :\t%e\t%e\t%e\t%e\t%e\t%e\n",Lat->Energy[Perturbed_P0].TotalEnergy[0], Lat->Energy[Perturbed_P1].TotalEnergy[0], Lat->Energy[Perturbed_P2].TotalEnergy[0], Lat->Energy[Perturbed_RxP0].TotalEnergy[0], Lat->Energy[Perturbed_RxP1].TotalEnergy[0], Lat->Energy[Perturbed_RxP2].TotalEnergy[0]);
663 if (R->DoUnOccupied) {
664 fprintf(stderr, "GapEnergy :\t%e\t%e\t%e\t%e\n",Lat->Energy[UnOccupied].AllTotalPsiEnergy[KineticEnergy],Lat->Energy[UnOccupied].AllTotalPsiEnergy[NonLocalEnergy], Lat->Energy[UnOccupied].AllTotalDensityEnergy[PseudoEnergy],Lat->Energy[UnOccupied].AllTotalDensityEnergy[HartreePotentialEnergy]);
665 fprintf(stderr, "TotalGapEnergy :\t%e\n",Lat->Energy[UnOccupied].TotalEnergy[0]);
666 fprintf(stderr, "BandGap :\t%e\n", Lat->Energy[UnOccupied].bandgap);
667 }
668 fprintf(stderr, "IonKinEnergy :\t%e\n",P->Ion.EKin);
669 fprintf(stderr, "Temperature :\t%e\n",P->Ion.ActualTemp);
670 }
671 if (P->Files.MeOutMes == 0 || doout == 0) return;
672 eout = F->EnergyFile;
673 if (eout != NULL) {
674 if (R->DoUnOccupied) {
675 fprintf(eout, "%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n",
676 P->R.t,
677 Lat->Energy[Occupied].TotalEnergy[0], Lat->Energy[UnOccupied].TotalEnergy[0],
678 Lat->Energy[Occupied].AllTotalPsiEnergy[KineticEnergy], Lat->Energy[Occupied].AllTotalPsiEnergy[NonLocalEnergy],
679 Lat->Energy[UnOccupied].AllTotalPsiEnergy[KineticEnergy]+Lat->Energy[UnOccupied].AllTotalPsiEnergy[NonLocalEnergy],
680 Lat->Energy[Occupied].AllTotalDensityEnergy[CorrelationEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[ExchangeEnergy],
681 Lat->Energy[Occupied].AllTotalDensityEnergy[PseudoEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[HartreePotentialEnergy],
682 Lat->Energy[UnOccupied].AllTotalDensityEnergy[PseudoEnergy]+Lat->Energy[UnOccupied].AllTotalDensityEnergy[HartreePotentialEnergy],
683 -Lat->Energy[Occupied].AllTotalIonsEnergy[GaussSelfEnergy],
684 Lat->Energy[Occupied].AllTotalIonsEnergy[EwaldEnergy],
685 P->Ion.EKin,
686 Lat->Energy[Occupied].TotalEnergy[0]+P->Ion.EKin);
687 } else {
688 fprintf(eout, "%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n",
689 P->R.t,
690 Lat->Energy[Occupied].TotalEnergy[0],
691 Lat->Energy[Occupied].AllTotalPsiEnergy[KineticEnergy], Lat->Energy[Occupied].AllTotalPsiEnergy[NonLocalEnergy],
692 Lat->Energy[Occupied].AllTotalDensityEnergy[CorrelationEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[ExchangeEnergy],
693 Lat->Energy[Occupied].AllTotalDensityEnergy[PseudoEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[HartreePotentialEnergy],
694 -Lat->Energy[Occupied].AllTotalIonsEnergy[GaussSelfEnergy],
695 Lat->Energy[Occupied].AllTotalIonsEnergy[EwaldEnergy],
696 P->Ion.EKin,
697 Lat->Energy[Occupied].TotalEnergy[0]+P->Ion.EKin);
698 }
699 fflush(eout);
700 }
701
702 // Output for testing convergence
703 if (!P->Par.me) {
704 if (R->DoPerturbation) {
705 if ((convergence = fopen("pcp.convergence.csv","r")) != NULL) { // only prepend with header if file doesn't exist yet
706 fclose(convergence);
707 convergence = fopen("pcp.convergence.csv","a");
708 } else {
709 convergence = fopen("pcp.convergence.csv","w");
710 if (convergence != NULL) {
711 fprintf(convergence,"Ecut\ta.u.\tMaxG\tN0\tN1\tN2\tChi00\tChi11\tChi22\tP0\tP1\tP2\tRxP0\tRxP1\tRxP2\tTE\t");
712 for (it=0; it < I->Max_Types; it++) { // over all types
713 fprintf(convergence,"Sigma00_%s\tSigma11_%s\tSigma22_%s\t",I->I[it].Symbol,I->I[it].Symbol,I->I[it].Symbol);
714 fprintf(convergence,"Sigma00_%s_rezi\tSigma11_%s_rezi\tSigma22_%s_rezi\t",I->I[it].Symbol,I->I[it].Symbol,I->I[it].Symbol);
715 }
716 fprintf(convergence,"Volume\tN^3\tSawStart\t");
717 fprintf(convergence,"Perturbedrxp-0\tPerturbedrxp-1\tPerturbedrxp-2\t");
718 fprintf(convergence,"Perturbedp-0\tPerturbedp-1\tPerturbedp-2\n");
719 }
720 }
721 if (convergence != NULL) {
722 fprintf(convergence,"%e\t%e\t%i\t%i\t%i\t%i\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t",
723 Lat->ECut,
724 sqrt(Lat->RealBasisSQ[0]), LevS->TotalAllMaxG, Lev0->Plan0.plan->N[0], Lev0->Plan0.plan->N[1], Lev0->Plan0.plan->N[2],
725 I->I[0].chi[0+0*NDIM], I->I[0].chi[1+1*NDIM], I->I[0].chi[2+2*NDIM],
726 Lat->Energy[Perturbed_P0].TotalEnergy[0], Lat->Energy[Perturbed_P1].TotalEnergy[0], Lat->Energy[Perturbed_P2].TotalEnergy[0],
727 Lat->Energy[Perturbed_RxP0].TotalEnergy[0], Lat->Energy[Perturbed_RxP1].TotalEnergy[0],Lat->Energy[Perturbed_RxP2].TotalEnergy[0],
728 Lat->Energy[Occupied].TotalEnergy[0]);
729 for (it=0; it < I->Max_Types; it++) { // over all types
730 //fprintf(convergence,"%e\t%e\t%e\t%e\t%e\t%e\t",
731 fprintf(convergence,"%e\t%e\t%e\t",
732 //I->I[it].sigma[0][0+0*NDIM], I->I[it].sigma[0][1+1*NDIM], I->I[it].sigma[0][2+2*NDIM],
733 I->I[it].sigma_rezi[0][0+0*NDIM], I->I[it].sigma_rezi[0][1+1*NDIM], I->I[it].sigma_rezi[0][2+2*NDIM]);
734 }
735 fprintf(convergence,"%e\t%i\t%e\t",Lat->Volume, Lev0->Plan0.plan->N[0]*Lev0->Plan0.plan->N[1]*Lev0->Plan0.plan->N[2],
736 P->Lat.SawtoothStart);
737 fprintf(convergence,"%e\t%e\t%e\t",Lat->E[Perturbed_RxP0].AllTotalPsiEnergy[Perturbed1_0Energy],Lat->E[Perturbed_RxP1].AllTotalPsiEnergy[Perturbed1_0Energy],Lat->E[Perturbed_RxP2].AllTotalPsiEnergy[Perturbed1_0Energy]);
738 fprintf(convergence,"%e\t%e\t%e\n",Lat->E[Perturbed_P0].AllTotalPsiEnergy[Perturbed1_0Energy],Lat->E[Perturbed_P1].AllTotalPsiEnergy[Perturbed1_0Energy],Lat->E[Perturbed_P2].AllTotalPsiEnergy[Perturbed1_0Energy]);
739 fclose(convergence);
740 }
741 } else {
742 if ((convergence = fopen("pcp.convergence.csv","r")) != NULL) { // only prepend with header if file doesn't exist yet
743 fclose(convergence);
744 convergence = fopen("pcp.convergence.csv","a");
745 } else {
746 convergence = fopen("pcp.convergence.csv","w");
747 if (convergence != NULL)
748 fprintf(convergence,"Ecut\ta.u.\tMaxG\tN0\tN1\tN2\tTE\tVolume\tN^3\tSawStart\n");
749 }
750 if (convergence != NULL) {
751 fprintf(convergence,"%e\t%e\t%i\t%i\t%i\t%i\t%e\t%e\t%i\t%e\n",
752 Lat->ECut,
753 sqrt(Lat->RealBasisSQ[0]), LevS->TotalAllMaxG, Lev0->Plan0.plan->N[0], Lev0->Plan0.plan->N[1], Lev0->Plan0.plan->N[2],
754 Lat->Energy[Occupied].TotalEnergy[0],
755 Lat->Volume, Lev0->Plan0.plan->N[0]*Lev0->Plan0.plan->N[1]*Lev0->Plan0.plan->N[2],
756 P->Lat.SawtoothStart);
757 fclose(convergence);
758 }
759 }
760 }
761}
Note: See TracBrowser for help on using the repository browser.