source: pcp/src/energy.c@ 774ae8

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

-initial commit
-Minimum set of files needed from ESPACK SVN repository
-Switch to three tantamount package parts instead of all relating to pcp (as at some time Ralf's might find inclusion as well)
  • Property mode set to 100644
File size: 34.2 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 SpeedMeasure(P, GapTime, StartTimeDo);
413 // === Calculate exchange and correlation energy ===
414 for (i = 0; i < Dens->LocalSizeR; i++) { // for each node in radial mesh
415 // put (corecorrected) densities in p, pUp, pDown
416 p = Dens->DensityArray[TotalDensity][i];
417 switch (Psi->PsiST) {
418 case SpinDouble:
419 pUp = 0.5*p;
420 pDown = 0.5*p;
421 break;
422 case SpinUp:
423 case SpinDown:
424 pUp = Dens->DensityArray[TotalUpDensity][i];
425 pDown = Dens->DensityArray[TotalDownDensity][i];
426 break;
427 default:
428 ;
429 }
430 if ((p < 0) || (pUp < 0) || (pDown < 0)) {
431 /*fprintf(stderr,"index %i pc %g p %g\n",i,Dens->DensityArray[CoreWaveDensity][i],p);*/
432 p = 0.0;
433 pUp = 0.0;
434 pDown = 0.0;
435 }
436 rs = Calcrs(EC,p);
437 rsUp = Calcrs(EC,pUp);
438 rsDown = Calcrs(EC,pDown);
439 zeta = CalcZeta(EC,pUp,pDown);
440 // Calculation with the densities and summation
441 SumEx += CalcSEXr(EC, rsUp, rsDown, pUp, pDown);
442 SumEc += CalcSECr(EC, rs, zeta, p);
443 }
444 E->AllLocalDensityEnergy[CorrelationEnergy] = Factor*SumEc;
445 E->AllLocalDensityEnergy[ExchangeEnergy] = Factor*SumEx;
446
447 // === Calculate electrostatic and local pseudopotential energy ===
448 SumEHP = 0.0;
449 SumEH = 0.0;
450 SumEG = 0.0;
451 SumEPS = 0.0;
452 // calculates energy of local pseudopotential
453 if (Lev0->GArray[0].GSq == 0.0) {
454 Index = Lev0->GArray[0].Index;
455 c_re(vp) = 0.0;
456 c_im(vp) = 0.0;
457 for (it = 0; it < I->Max_Types; it++) {
458 c_re(vp) += (c_re(I->I[it].SFactor[0])*PP->phi_ps_loc[it][0]);
459 c_im(vp) += (c_im(I->I[it].SFactor[0])*PP->phi_ps_loc[it][0]);
460 }
461 SumEPS += (c_re(Dens->DensityCArray[GapDensity][Index])*c_re(vp) +
462 c_im(Dens->DensityCArray[GapDensity][Index])*c_im(vp))*R->HGcFactor;
463 s++;
464 }
465 for (g=s; g < Lev0->MaxG; g++) {
466 Index = Lev0->GArray[g].Index;
467 c_re(vp) = 0.0;
468 c_im(vp) = 0.0;
469 c_re(rp) = 0.0;
470 c_im(rp) = 0.0;
471 Fac = 2.*PI/(Lev0->GArray[g].GSq); // 4*.. -> 2*...: 1/2 here ! (due to lacking dependence of first n^{tot} (r) on gap psis)
472 for (it = 0; it < I->Max_Types; it++) {
473 c_re(vp) += (c_re(I->I[it].SFactor[g])*PP->phi_ps_loc[it][g]);
474 c_im(vp) += (c_im(I->I[it].SFactor[g])*PP->phi_ps_loc[it][g]);
475 c_re(rp) += (c_re(I->I[it].SFactor[g])*PP->FacGauss[it][g]);
476 c_im(rp) += (c_im(I->I[it].SFactor[g])*PP->FacGauss[it][g]);
477 }
478 c_re(rhog) = c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_re(rp);
479 c_im(rhog) = c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_im(rp);
480 c_re(rhoe) = c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
481 c_im(rhoe) = c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
482 c_re(rhoe_gap) = c_re(Dens->DensityCArray[GapDensity][Index])*R->HGcFactor;
483 c_im(rhoe_gap) = c_im(Dens->DensityCArray[GapDensity][Index])*R->HGcFactor;
484 //c_re(PP->VCoulombc[g]) = Fac*c_re(rhog);
485 //c_im(PP->VCoulombc[g]) = -Fac*c_im(rhog);
486 c_re(HG[Index]) = c_re(vp) + Fac * (c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_re(rp));
487 c_im(HG[Index]) = c_im(vp) + Fac * (c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor+c_im(rp));
488 /*if (P->first) */
489 //SumEG += Fac*(c_re(rp)*c_re(rp)+c_im(rp)*c_im(rp)); // Coulomb energy
490 // changed: took out gaussian part! rhog -> rhoe
491 SumEHP += 2*Fac*(c_re(rhoe)*c_re(rhoe_gap)+c_im(rhoe)*c_im(rhoe_gap)); // E_ES first part
492 //SumEH += Fac*(c_re(rhoe)*c_re(rhoe)+c_im(rhoe)*c_im(rhoe));
493 SumEPS += 2.*(c_re(Dens->DensityCArray[GapDensity][Index])*c_re(vp)+
494 c_im(Dens->DensityCArray[GapDensity][Index])*c_im(vp))*R->HGcFactor;
495 }
496 //
497 for (i=0; i<Lev0->MaxDoubleG; i++) {
498 HG[Lev0->DoubleG[2*i+1]].re = HG[Lev0->DoubleG[2*i]].re;
499 HG[Lev0->DoubleG[2*i+1]].im = -HG[Lev0->DoubleG[2*i]].im;
500 }
501 /*if (P->first) */
502 //E->AllLocalDensityEnergy[GaussEnergy] = SumEG*Lat->Volume;
503 E->AllLocalDensityEnergy[PseudoEnergy] = SumEPS*Lat->Volume; // local pseudopotential
504 E->AllLocalDensityEnergy[HartreePotentialEnergy] = SumEHP*Lat->Volume; // Gauss n^Gauss and Hartree n potential
505 //E->AllLocalDensityEnergy[HartreeEnergy] = SumEH*Lat->Volume;
506 SpeedMeasure(P, GapTime, StopTimeDo);
507
508 // === Calculate Gradient ===
509 CalculateGradientNoRT(P,
510 LevS->LPsi->LocalPsi[R->ActualLocalPsiNo],
511 Lat->Psi.LocalPsiStatus[R->ActualLocalPsiNo].PsiFactor,
512 &Lat->Psi.AddData[R->ActualLocalPsiNo].Lambda,
513 P->PP.fnl[R->ActualLocalPsiNo],
514 P->Grad.GradientArray[ActualGradient],
515 P->Grad.GradientArray[OldActualGradient],
516 P->Grad.GradientArray[HcGradient]);
517}
518
519/** Calculates the energy of a wave function psi.
520 * Calculates kinetc and non local energies, taking care of possible
521 * Riemann tensor usage.
522 * \param *P Problem at hand
523 * \param p number of wave function Psi
524 * \param UseInitTime Whether speed measuring uses InitLocTime (Init) or LocTime (Update)
525 * \sa EnergyGNSP(), EnergyLapNSP()
526 * \warning call EnergyAllReduce() afterwards!
527 */
528void CalculatePsiEnergy(struct Problem *P, const int p, const int UseInitTime)
529{
530 struct Lattice *Lat = &(P->Lat);
531 switch (Lat->RT.ActualUse) {
532 case inactive:
533 case standby:
534 SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StartTimeDo);
535 CalculateKineticEnergyNoRT(P, p);
536 SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StopTimeDo);
537 SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StartTimeDo);
538 CalculateNonLocalEnergyNoRT(P, p);
539 SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StopTimeDo);
540 break;
541 case active:
542 SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StartTimeDo);
543 CalculateKineticEnergyUseRT(P, p);
544 SpeedMeasure(P, (UseInitTime == 1 ? InitLocTime : LocTime), StopTimeDo);
545 SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StartTimeDo);
546 CalculateNonLocalEnergyUseRT(P, p);
547 SpeedMeasure(P, (UseInitTime == 1 ? InitNonLocTime : NonLocTime), StopTimeDo);
548 break;
549 }
550 /* Achtung danach noch EnergyAllReduce aufrufen */
551}
552
553/** Initialization of wave function energy calculation.
554 * For every wave function CalculatePsiEnergy() is called,
555 * Energy::AllLocalPsiEnergy is set to sum of these.
556 * \param *P Problem at hand
557 * \param type minimisation type PsiTypeTag to calculate psi energy for
558 * \warning call EnergyAllReduce() afterwards!
559 */
560void InitPsiEnergyCalculation(struct Problem *P, enum PsiTypeTag type)
561{
562 struct Lattice *Lat = &(P->Lat);
563 int p,i, oldtype = P->R.CurrentMin;
564 SetCurrentMinState(P,type);
565 for (p=0; p < Lat->Psi.LocalNo; p++) {
566 if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
567 CalculatePsiEnergy(P, p, 1);
568 }
569 for (i=0; i< MAXALLPSIENERGY; i++) {
570 Lat->E->AllLocalPsiEnergy[i] = 0.0;
571 for (p=0; p < Lat->Psi.LocalNo; p++)
572 if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
573 Lat->E->AllLocalPsiEnergy[i] += Lat->E->PsiEnergy[i][p];
574 }
575 SetCurrentMinState(P,oldtype);
576 /* Achtung danach noch EnergyAllReduce aufrufen */
577}
578
579/** Updating wave function energy.
580 * CalculatePsiEnergy() is called for RunStruct::OldActualLocalPsiNo,
581 * Energy::AllLocalPsiEnergy is set to sum of all Psi energies.
582 * \param *P Problem at hand
583 * \warning call EnergyAllReduce() afterwards!
584 */
585void UpdatePsiEnergyCalculation(struct Problem *P)
586{
587 struct Lattice *Lat = &(P->Lat);
588 struct RunStruct *R = &P->R;
589 int p = R->OldActualLocalPsiNo, i;
590 if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
591 CalculatePsiEnergy(P, p, 0);
592 for (i=0; i< MAXALLPSIENERGY; i++) {
593 Lat->E->AllLocalPsiEnergy[i] = 0.0;
594 for (p=0; p < P->Lat.Psi.LocalNo; p++)
595 if (P->Lat.Psi.LocalPsiStatus[p].PsiType == P->R.CurrentMin)
596 Lat->E->AllLocalPsiEnergy[i] += Lat->E->PsiEnergy[i][p];
597 }
598 /* Achtung danach noch EnergyAllReduce aufrufen */
599}
600
601/** Output of energies.
602 * Prints either Time, PsiEnergy(Up/Down), DensityEnergy, IonsEnergy, KineticEnergy,
603 * TotalEnergy and Temperature to screen or<br>
604 * Time, TotalEnergy, NonLocalEnergy, CorrelationEnergy, ExchangeEnergy,
605 * PseudoEnergy, HartreePotentialEnergy, GaussSelfEnergy, EwaldEnergy,
606 * KineticEnergy and TotalEnergy+KineticEnergy to the energy file
607 * \param *P Problem at hand
608 * \param doout 1 - print stuff to stderr, 0 - put compact form into energy file
609 */
610void EnergyOutput(struct Problem *P, int doout)
611{
612 struct Lattice *Lat = &P->Lat;
613 struct RunStruct *R = &P->R;
614 struct LatticeLevel *Lev0 = R->Lev0;
615 struct LatticeLevel *LevS = R->LevS;
616 struct Ions *I = &P->Ion;
617 struct Psis *Psi = &Lat->Psi;
618 struct FileData *F = &P->Files;
619 FILE *eout = NULL;
620 FILE *convergence = NULL;
621 int i, it;
622 if (P->Call.out[LeaderOut] && (P->Par.me == 0) && doout) {
623 fprintf(stderr, "Time T: :\t%e\n",P->R.t);
624 fprintf(stderr, "PsiEnergy :");
625 for (i=0; i<MAXALLPSIENERGY; i++) {
626 fprintf(stderr," %e",Lat->Energy[Occupied].AllTotalPsiEnergy[i]);
627 }
628 fprintf(stderr,"\n");
629 if (Psi->PsiST != SpinDouble) {
630 fprintf(stderr, "PsiEnergyUp :");
631 for (i=0; i<MAXALLPSIENERGY; i++) {
632 fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllUpPsiEnergy[i]);
633 }
634 fprintf(stderr,"\n");
635 fprintf(stderr, "PsiEnergyDown:");
636 for (i=0; i<MAXALLPSIENERGY; i++) {
637 fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllDownPsiEnergy[i]);
638 }
639 fprintf(stderr,"\n");
640 }
641 fprintf(stderr, "DensityEnergy:");
642 for (i=0; i<MAXALLDENSITYENERGY; i++) {
643 fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllTotalDensityEnergy[i]);
644 }
645 fprintf(stderr,"\n");
646 fprintf(stderr, "IonsEnergy :");
647 for (i=0; i<MAXALLIONSENERGY; i++) {
648 fprintf(stderr,"\t%e",Lat->Energy[Occupied].AllTotalIonsEnergy[i]);
649 }
650 fprintf(stderr,"\n");
651 fprintf(stderr, "TotalEnergy :\t%e\n",Lat->Energy[Occupied].TotalEnergy[0]);
652 if (R->DoPerturbation)
653 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]);
654 if (R->DoUnOccupied) {
655 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]);
656 fprintf(stderr, "TotalGapEnergy :\t%e\n",Lat->Energy[UnOccupied].TotalEnergy[0]);
657 fprintf(stderr, "BandGap :\t%e\n", Lat->Energy[UnOccupied].bandgap);
658 }
659 fprintf(stderr, "IonKinEnergy :\t%e\n",P->Ion.EKin);
660 fprintf(stderr, "Temperature :\t%e\n",P->Ion.ActualTemp);
661 }
662 if (P->Files.MeOutMes == 0 || doout == 0) return;
663 eout = F->EnergyFile;
664 if (eout != NULL) {
665 if (R->DoUnOccupied) {
666 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",
667 P->R.t,
668 Lat->Energy[Occupied].TotalEnergy[0], Lat->Energy[UnOccupied].TotalEnergy[0],
669 Lat->Energy[Occupied].AllTotalPsiEnergy[KineticEnergy], Lat->Energy[Occupied].AllTotalPsiEnergy[NonLocalEnergy],
670 Lat->Energy[UnOccupied].AllTotalPsiEnergy[KineticEnergy]+Lat->Energy[UnOccupied].AllTotalPsiEnergy[NonLocalEnergy],
671 Lat->Energy[Occupied].AllTotalDensityEnergy[CorrelationEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[ExchangeEnergy],
672 Lat->Energy[Occupied].AllTotalDensityEnergy[PseudoEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[HartreePotentialEnergy],
673 Lat->Energy[UnOccupied].AllTotalDensityEnergy[PseudoEnergy]+Lat->Energy[UnOccupied].AllTotalDensityEnergy[HartreePotentialEnergy],
674 -Lat->Energy[Occupied].AllTotalIonsEnergy[GaussSelfEnergy],
675 Lat->Energy[Occupied].AllTotalIonsEnergy[EwaldEnergy],
676 P->Ion.EKin,
677 Lat->Energy[Occupied].TotalEnergy[0]+P->Ion.EKin);
678 } else {
679 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",
680 P->R.t,
681 Lat->Energy[Occupied].TotalEnergy[0],
682 Lat->Energy[Occupied].AllTotalPsiEnergy[KineticEnergy], Lat->Energy[Occupied].AllTotalPsiEnergy[NonLocalEnergy],
683 Lat->Energy[Occupied].AllTotalDensityEnergy[CorrelationEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[ExchangeEnergy],
684 Lat->Energy[Occupied].AllTotalDensityEnergy[PseudoEnergy], Lat->Energy[Occupied].AllTotalDensityEnergy[HartreePotentialEnergy],
685 -Lat->Energy[Occupied].AllTotalIonsEnergy[GaussSelfEnergy],
686 Lat->Energy[Occupied].AllTotalIonsEnergy[EwaldEnergy],
687 P->Ion.EKin,
688 Lat->Energy[Occupied].TotalEnergy[0]+P->Ion.EKin);
689 }
690 fflush(eout);
691 }
692
693 // Output for testing convergence
694 if (!P->Par.me) {
695 if (R->DoPerturbation) {
696 if ((convergence = fopen("pcp.convergence.csv","r")) != NULL) { // only prepend with header if file doesn't exist yet
697 fclose(convergence);
698 convergence = fopen("pcp.convergence.csv","a");
699 } else {
700 convergence = fopen("pcp.convergence.csv","w");
701 if (convergence != NULL) {
702 fprintf(convergence,"Ecut\ta.u.\tMaxG\tN0\tN1\tN2\tChi00\tChi11\tChi22\tP0\tP1\tP2\tRxP0\tRxP1\tRxP2\tTE\t");
703 for (it=0; it < I->Max_Types; it++) { // over all types
704 fprintf(convergence,"Sigma00_%s\tSigma11_%s\tSigma22_%s\t",I->I[it].Symbol,I->I[it].Symbol,I->I[it].Symbol);
705 fprintf(convergence,"Sigma00_%s_rezi\tSigma11_%s_rezi\tSigma22_%s_rezi\t",I->I[it].Symbol,I->I[it].Symbol,I->I[it].Symbol);
706 }
707 fprintf(convergence,"Volume\tN^3\tSawStart\t");
708 fprintf(convergence,"Perturbedrxp-0\tPerturbedrxp-1\tPerturbedrxp-2\t");
709 fprintf(convergence,"Perturbedp-0\tPerturbedp-1\tPerturbedp-2\n");
710 }
711 }
712 if (convergence != NULL) {
713 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",
714 Lat->ECut,
715 Lat->RealBasisQ[0], LevS->TotalAllMaxG, Lev0->Plan0.plan->N[0], Lev0->Plan0.plan->N[1], Lev0->Plan0.plan->N[2],
716 I->I[0].chi[0+0*NDIM], I->I[0].chi[1+1*NDIM], I->I[0].chi[2+2*NDIM],
717 Lat->Energy[Perturbed_P0].TotalEnergy[0], Lat->Energy[Perturbed_P1].TotalEnergy[0], Lat->Energy[Perturbed_P2].TotalEnergy[0],
718 Lat->Energy[Perturbed_RxP0].TotalEnergy[0], Lat->Energy[Perturbed_RxP1].TotalEnergy[0],Lat->Energy[Perturbed_RxP2].TotalEnergy[0],
719 Lat->Energy[Occupied].TotalEnergy[0]);
720 for (it=0; it < I->Max_Types; it++) { // over all types
721 //fprintf(convergence,"%e\t%e\t%e\t%e\t%e\t%e\t",
722 fprintf(convergence,"%e\t%e\t%e\t",
723 //I->I[it].sigma[0][0+0*NDIM], I->I[it].sigma[0][1+1*NDIM], I->I[it].sigma[0][2+2*NDIM],
724 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]);
725 }
726 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],
727 P->Lat.SawtoothStart);
728 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]);
729 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]);
730 fclose(convergence);
731 }
732 } else {
733 if ((convergence = fopen("pcp.convergence.csv","r")) != NULL) { // only prepend with header if file doesn't exist yet
734 fclose(convergence);
735 convergence = fopen("pcp.convergence.csv","a");
736 } else {
737 convergence = fopen("pcp.convergence.csv","w");
738 if (convergence != NULL)
739 fprintf(convergence,"Ecut\ta.u.\tMaxG\tN0\tN1\tN2\tTE\tVolume\tN^3\tSawStart\n");
740 }
741 if (convergence != NULL) {
742 fprintf(convergence,"%e\t%e\t%i\t%i\t%i\t%i\t%e\t%e\t%i\t%e\n",
743 Lat->ECut,
744 Lat->RealBasisQ[0], LevS->TotalAllMaxG, Lev0->Plan0.plan->N[0], Lev0->Plan0.plan->N[1], Lev0->Plan0.plan->N[2],
745 Lat->Energy[Occupied].TotalEnergy[0],
746 Lat->Volume, Lev0->Plan0.plan->N[0]*Lev0->Plan0.plan->N[1]*Lev0->Plan0.plan->N[2],
747 P->Lat.SawtoothStart);
748 fclose(convergence);
749 }
750 }
751 }
752}
Note: See TracBrowser for help on using the repository browser.