source: pcp/src/test.c@ 460e95

Last change on this file since 460e95 was a0bcf1, checked in by Frederik Heber <heber@…>, 18 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: 13.4 KB
RevLine 
[a0bcf1]1void CalculatePerturbedEnergy(struct Problem *P, int l, int DoGradient)
2{
3 struct Lattice *Lat = &P->Lat;
4 struct Psis *Psi = &Lat->Psi;
5 struct Energy *E = Lat->E;
6 struct PseudoPot *PP = &P->PP;
7 struct RunStruct *R = &P->R;
8 struct LatticeLevel *LevS = R->LevS;
9 int state = R->CurrentMin;
10 int l_normal = Psi->TypeStartIndex[Occupied] + (l - Psi->TypeStartIndex[state]); // offset l to \varphi_l^{(0)}
11 int ActNum = l - Psi->TypeStartIndex[state] + Psi->TypeStartIndex[1] * Psi->LocalPsiStatus[l].my_color_comm_ST_Psi;
12 int g, i, m, j;
13 double lambda, Lambda;
14 double RElambda10, RELambda10;
15 double RElambda01, RELambda01;
16 double IMlambda10, IMLambda10;
17 double IMlambda01, IMLambda01;
18 fftw_complex *source = LevS->LPsi->LocalPsi[l];
19 fftw_complex *grad = P->Grad.GradientArray[ActualGradient];
20 fftw_complex *Hc_grad = P->Grad.GradientArray[HcGradient];
21 fftw_complex *H1c_grad = P->Grad.GradientArray[H1cGradient];
22 fftw_complex *TempPsi_0 = H1c_grad;
23 fftw_complex *TempPsi_1 = P->Grad.GradientArray[TempGradient];
24 fftw_complex *varphi_1, *varphi_0;
25 struct OnePsiElement *OnePsiB, *LOnePsiB;
26 fftw_complex *LPsiDatB=NULL;
27 int ElementSize = (sizeof(fftw_complex) / sizeof(double)), RecvSource;
28 MPI_Status status;
29
30 // ============ Calculate H^(0) psi^(1) =============================
31 SetArrayToDouble0((double *)Hc_grad,2*R->InitLevS->MaxG);
32 ApplyTotalHamiltonian(P,source,Hc_grad, PP->fnl[l], 1, 1);
33
34 // ============ ENERGY FUNCTIONAL Evaluation PART 1 ================
35 varphi_0 = LevS->LPsi->LocalPsi[l_normal];
36 varphi_1 = LevS->LPsi->LocalPsi[l];
37 switch (state) {
38 case Perturbed_P0:
39 CalculatePerturbationOperator_P(P,varphi_0,TempPsi_0,0,0); // \nabla_0 | \varphi_l^{(0)} \rangle
40 CalculatePerturbationOperator_P(P,varphi_1,TempPsi_1,0,0); // \nabla_0 | \varphi_l^{(1)} \rangle
41 break;
42 case Perturbed_P1:
43 CalculatePerturbationOperator_P(P,varphi_0,TempPsi_0,1,0); // \nabla_1 | \varphi_l^{(0)} \rangle
44 CalculatePerturbationOperator_P(P,varphi_1,TempPsi_1,1,0); // \nabla_1 | \varphi_l^{(1)} \rangle
45 break;
46 case Perturbed_P2:
47 CalculatePerturbationOperator_P(P,varphi_0,TempPsi_0,2,0); // \nabla_1 | \varphi_l^{(0)} \rangle
48 CalculatePerturbationOperator_P(P,varphi_1,TempPsi_1,2,0); // \nabla_1 | \varphi_l^{(1)} \rangle
49 break;
50 case Perturbed_RxP0:
51 CalculatePerturbationOperator_RxP(P,varphi_0,TempPsi_0,l,0,0); // r \times \nabla | \varphi_l^{(0)} \rangle
52 CalculatePerturbationOperator_RxP(P,varphi_1,TempPsi_1,l,0,0); // r \times \nabla | \varphi_l^{(1)} \rangle
53 //CalculatePerturbationOperator_R(P,varphi_0,TempPsi_0,0,l); // r \times \nabla | \varphi_l^{(0)} \rangle
54 //CalculatePerturbationOperator_R(P,varphi_1,TempPsi_1,0,l); // r \times \nabla | \varphi_l^{(1)} \rangle
55 break;
56 case Perturbed_RxP1:
57 CalculatePerturbationOperator_RxP(P,varphi_0,TempPsi_0,l,1,0); // r \times \nabla | \varphi_l^{(0)} \rangle
58 CalculatePerturbationOperator_RxP(P,varphi_1,TempPsi_1,l,1,0); // r \times \nabla | \varphi_l^{(1)} \rangle
59 //CalculatePerturbationOperator_R(P,varphi_0,TempPsi_0,1,l); // r \times \nabla | \varphi_l^{(0)} \rangle
60 //CalculatePerturbationOperator_R(P,varphi_1,TempPsi_1,1,l); // r \times \nabla | \varphi_l^{(1)} \rangle
61 break;
62 case Perturbed_RxP2:
63 CalculatePerturbationOperator_RxP(P,varphi_0,TempPsi_0,l,2,0); // r \times \nabla | \varphi_l^{(0)} \rangle
64 CalculatePerturbationOperator_RxP(P,varphi_1,TempPsi_1,l,2,0); // r \times \nabla | \varphi_l^{(1)} \rangle
65 //CalculatePerturbationOperator_R(P,varphi_0,TempPsi_0,2,l); // r \times \nabla | \varphi_l^{(0)} \rangle
66 //CalculatePerturbationOperator_R(P,varphi_1,TempPsi_1,2,l); // r \times \nabla | \varphi_l^{(1)} \rangle
67 break;
68 default:
69 fprintf(stderr,"(%i) CalculatePerturbedEnergy called whilst not within perturbation run: CurrentMin = %i !\n",P->Par.me, R->CurrentMin);
70 break;
71 }
72
73 // ============ GRADIENT and EIGENVALUE Evaluation Part 1==============
74 lambda = 0.0;
75 if ((DoGradient) && (grad != NULL)) {
76 g = 0;
77 if (LevS->GArray[0].GSq == 0.0) {
78 lambda += Hc_grad[0].re*source[0].re;
79 grad[0].re = -(Hc_grad[0].re + TempPsi_0[0].re);
80 grad[0].im = -(Hc_grad[0].im + TempPsi_0[0].im);
81 g++;
82 }
83 for (;g<LevS->MaxG;g++) {
84 lambda += 2.*(Hc_grad[g].re*source[g].re + Hc_grad[g].im*source[g].im);
85 grad[g].re = -(Hc_grad[g].re + TempPsi_0[g].re);
86 grad[g].im = -(Hc_grad[g].im + TempPsi_0[g].im);
87 }
88
89 m = -1;
90 for (j=0; j < Psi->MaxPsiOfType+P->Par.Max_me_comm_ST_PsiT; j++) { // go through all wave functions
91 OnePsiB = &Psi->AllPsiStatus[j]; // grab OnePsiB
92 if (OnePsiB->PsiType == state) { // drop all but the ones of current min state
93 m++; // increase m if it is type-specific wave function
94 if (OnePsiB->my_color_comm_ST_Psi == P->Par.my_color_comm_ST_Psi) // local?
95 LOnePsiB = &Psi->LocalPsiStatus[OnePsiB->MyLocalNo];
96 else
97 LOnePsiB = NULL;
98 if (LOnePsiB == NULL) { // if it's not local ... receive it from respective process into TempPsi
99 RecvSource = OnePsiB->my_color_comm_ST_Psi;
100 MPI_Recv( LevS->LPsi->TempPsi, LevS->MaxG*ElementSize, MPI_DOUBLE, RecvSource, PerturbedTag, P->Par.comm_ST_PsiT, &status );
101 LPsiDatB=LevS->LPsi->TempPsi;
102 } else { // .. otherwise send it to all other processes (Max_me... - 1)
103 for (i=0;i<P->Par.Max_me_comm_ST_PsiT;i++)
104 if (i != OnePsiB->my_color_comm_ST_Psi)
105 MPI_Send( LevS->LPsi->LocalPsi[OnePsiB->MyLocalNo], LevS->MaxG*ElementSize, MPI_DOUBLE, i, PerturbedTag, P->Par.comm_ST_PsiT);
106 LPsiDatB=LevS->LPsi->LocalPsi[OnePsiB->MyLocalNo];
107 } // LPsiDatB is now set to the coefficients of OnePsi either stored or MPI_Received
108
109 g = 0;
110 if (LevS->GArray[0].GSq == 0.0) { // perform the summation
111 grad[0].re += Lat->Psi.lambda[ActNum][m]*LPsiDatB[0].re;
112 grad[0].im += Lat->Psi.lambda[ActNum][m]*LPsiDatB[0].im;
113 g++;
114 }
115 for (;g<LevS->MaxG;g++) {
116 grad[g].re += Lat->Psi.lambda[ActNum][m]*LPsiDatB[g].re;
117 grad[g].im += Lat->Psi.lambda[ActNum][m]*LPsiDatB[g].im;
118 }
119 }
120 }
121 } else {
122 g = 0;
123 if (LevS->GArray[0].GSq == 0.0) {
124 lambda += Hc_grad[0].re*source[0].re;
125 g++;
126 }
127 for (;g<LevS->MaxG;g++)
128 lambda += 2.*(Hc_grad[g].re*source[g].re + Hc_grad[g].im*source[g].im);
129 }
130 MPI_Allreduce ( &lambda, &Lambda, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
131 //fprintf(stderr,"(%i) Lambda[%i] = %lg\n",P->Par.me, l, Lambda);
132 Lat->Psi.AddData[l].Lambda = Lambda;
133
134 // ============ ENERGY FUNCTIONAL Evaluation PART 2 ================
135 // varphi_1 jas negative symmetry, returning TemPsi_0 from CalculatePerturbedOperator also, thus real part of scalar product
136 // "-" due to purely imaginary wave function is on left hand side, thus becomes complex conjugated: i -> -i
137 // (-i goes into pert. op., "-" remains when on right hand side)
138 IMlambda10 = GradImSP(P,LevS,varphi_1,TempPsi_0) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
139 IMlambda01 = -GradImSP(P,LevS,varphi_0,TempPsi_1) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
140 RElambda10 = GradSP(P,LevS,varphi_1,TempPsi_0) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
141 RElambda01 = -GradSP(P,LevS,varphi_0,TempPsi_1) * sqrt(Psi->LocalPsiStatus[l].PsiFactor * Psi->LocalPsiStatus[l_normal].PsiFactor);
142
143 MPI_Allreduce ( &RElambda10, &RELambda10, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
144 MPI_Allreduce ( &RElambda01, &RELambda01, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
145 MPI_Allreduce ( &IMlambda10, &IMLambda10, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
146 MPI_Allreduce ( &IMlambda01, &IMLambda01, 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);
147
148 E->PsiEnergy[Perturbed1_0Energy][l] = RELambda10;
149 E->PsiEnergy[Perturbed0_1Energy][l] = RELambda01;
150 if (P->Par.me == 0) {
151 fprintf(stderr,"RE.Lambda10[%i] = %lg\t RE.Lambda01[%i] = %lg\n", l, RELambda10, l, RELambda01);
152 fprintf(stderr,"IM.Lambda10[%i] = %lg\t IM.Lambda01[%i] = %lg\n", l, IMLambda10, l, IMLambda01);
153 }
154}
155
156void CalculatePerturbationOperator_RxP(struct Problem *P, fftw_complex *source, fftw_complex *dest, int l, int index, int symmetry)
157{
158 struct RunStruct *R = &P->R;
159 struct LatticeLevel *Lev0 = R->Lev0;
160 struct LatticeLevel *LevS = R->LevS;
161 struct Lattice *Lat = &P->Lat;
162 struct fft_plan_3d *plan = Lat->plan;
163 struct Density *Dens0 = Lev0->Dens;
164 fftw_complex *tempdestRC = Dens0->DensityCArray[CurrentDensity0];
165 fftw_real *tempdestR = (fftw_real *) tempdestRC;
166 fftw_complex *tempdestRC2 = Dens0->DensityCArray[CurrentDensity1];
167 fftw_real *tempdestR2 = (fftw_real *) tempdestRC2;
168 fftw_complex *work = Dens0->DensityCArray[TempDensity];
169 fftw_complex *PsiC = (fftw_complex *) Dens0->DensityCArray[ActualPsiDensity];
170 fftw_real *PsiCR = (fftw_real *) PsiC;
171 fftw_real *RealPhiR = (fftw_real *) Dens0->DensityArray[TempDensity];
172 fftw_real *RealPhiR2 = (fftw_real *) Dens0->DensityArray[Temp2Density];
173 fftw_complex *posfac, *destsnd, *destrcv, *destsnd2, *destrcv2;
174 double x[NDIM], fac[NDIM], *WCentre;
175 int n[NDIM], N0, n0, g, Index, pos, iS, i0;
176
177 // init pointers and values
178 int myPE = P->Par.me_comm_ST_Psi;
179 int N[NDIM], NUp[NDIM];
180 N[0] = LevS->Plan0.plan->N[0];
181 N[1] = LevS->Plan0.plan->N[1];
182 N[2] = LevS->Plan0.plan->N[2];
183 NUp[0] = LevS->NUp[0];
184 NUp[1] = LevS->NUp[1];
185 NUp[2] = LevS->NUp[2];
186 double FFTFactor = 1./LevS->MaxN;
187 N0 = LevS->Plan0.plan->local_nx;
188 WCentre = Lat->Psi.AddData[l].WannierCentre;
189
190 // blow up source coefficients
191 SetArrayToDouble0((double *)tempdestRC ,Dens0->TotalSize*2);
192 SetArrayToDouble0((double *)tempdestRC2,Dens0->TotalSize*2);
193 SetArrayToDouble0((double *)RealPhiR ,Dens0->TotalSize*2);
194 SetArrayToDouble0((double *)RealPhiR2,Dens0->TotalSize*2);
195 SetArrayToDouble0((double *)PsiC,Dens0->TotalSize*2);
196 for (g=0; g<LevS->MaxG; g++) {
197 Index = LevS->GArray[g].Index;
198 posfac = &LevS->PosFactorUp[LevS->MaxNUp*g];
199 destrcv = &tempdestRC[LevS->MaxNUp*Index];
200 destrcv2 = &tempdestRC2[LevS->MaxNUp*Index];
201 for (pos=0; pos<LevS->MaxNUp; pos++) {
202 destrcv [pos].re = LevS->GArray[g].G[cross(index,1)]*(( source[g].im)*posfac[pos].re-(-source[g].re)*posfac[pos].im);
203 destrcv [pos].im = LevS->GArray[g].G[cross(index,1)]*(( source[g].im)*posfac[pos].im+(-source[g].re)*posfac[pos].re);
204 destrcv2[pos].re = LevS->GArray[g].G[cross(index,3)]*(( source[g].im)*posfac[pos].re-(-source[g].re)*posfac[pos].im);
205 destrcv2[pos].im = LevS->GArray[g].G[cross(index,3)]*(( source[g].im)*posfac[pos].im+(-source[g].re)*posfac[pos].re);
206 }
207 }
208 for (g=0; g<LevS->MaxDoubleG; g++) {
209 destsnd = &tempdestRC [LevS->DoubleG[2*g]*LevS->MaxNUp];
210 destrcv = &tempdestRC [LevS->DoubleG[2*g+1]*LevS->MaxNUp];
211 destsnd2 = &tempdestRC2[LevS->DoubleG[2*g]*LevS->MaxNUp];
212 destrcv2 = &tempdestRC2[LevS->DoubleG[2*g+1]*LevS->MaxNUp];
213 for (pos=0; pos<LevS->MaxNUp; pos++) {
214 destrcv [pos].re = destsnd [pos].re;
215 destrcv [pos].im = -destsnd [pos].im;
216 destrcv2[pos].re = destsnd2[pos].re;
217 destrcv2[pos].im = -destsnd2[pos].im;
218 }
219 }
220
221 // fourier transform blown up wave function
222 fft_3d_complex_to_real(plan,LevS->LevelNo, FFTNFUp, tempdestRC , work);
223 DensityRTransformPos(LevS,tempdestR ,RealPhiR );
224 fft_3d_complex_to_real(plan,LevS->LevelNo, FFTNFUp, tempdestRC2, work);
225 DensityRTransformPos(LevS,tempdestR2,RealPhiR2);
226
227 //fft_Psi(P,source,RealPhiR ,cross(index,1),6);
228 //fft_Psi(P,source,RealPhiR2,cross(index,3),6);
229
230 // for every point on the real grid multiply with component of position vector
231 for (n0=0; n0<N0; n0++)
232 for (n[1]=0; n[1]<N[1]; n[1]++)
233 for (n[2]=0; n[2]<N[2]; n[2]++) {
234 n[0] = n0 + N0 * myPE;
235 fac[0] = (double)(n[0])/(double)((N[0]));
236 fac[1] = (double)(n[1])/(double)((N[1]));
237 fac[2] = (double)(n[2])/(double)((N[2]));
238 RMat33Vec3(x,Lat->RealBasis,fac);
239 iS = n[2] + N[2]*(n[1] + N[1]*n0); // mind splitting of x axis due to multiple processes
240 i0 = n[2]*NUp[2]+N[2]*NUp[2]*(n[1]*NUp[1]+N[1]*NUp[1]*n0*NUp[0]);
241 //PsiCR[iS] = ((double)n[cross(index,0)]/(double)N[cross(index,0)]*Lat->RealBasisQ[cross(index,0)] - WCentre[cross(index,0)])*RealPhiR[i0] - ((double)n[cross(index,2)]/(double)N[cross(index,2)]*Lat->RealBasisQ[cross(index,2)] - WCentre[cross(index,2)])*RealPhiR2[i0];
242 PsiCR[iS] =
243 MinImageConv(Lat,x[cross(index,0)],WCentre[cross(index,0)],cross(index,0)) * RealPhiR [i0]
244 - MinImageConv(Lat,x[cross(index,2)],WCentre[cross(index,2)],cross(index,2)) * RealPhiR2[i0];
245 //tmp += truedist(Lat,x[index_r],WCentre[index_r],index_r) * RealPhiR[i0];
246 //tmp += sawtooth(P,truedist(Lat,x[index_r],WCentre[index_r],index_r), index_r)*RealPhiR[i0];
247 //(Fehler mit falschem Ort ist vor dieser Stelle!): ueber result = RealPhiR[i0] * (x[index_r]) * RealPhiR[i0]; gecheckt
248 }
249
250 // inverse fourier transform
251 fft_3d_real_to_complex(plan,LevS->LevelNo, FFTNF1, PsiC, work);
252
253 // copy to destination array
254 for (g=0; g<LevS->MaxG; g++) {
255 Index = LevS->GArray[g].Index;
256 dest[g].re = ( PsiC[Index].re)*FFTFactor;
257 dest[g].im = ( PsiC[Index].im)*FFTFactor;
258 if (LevS->GArray[g].GSq == 0)
259 dest[g].im = 0; // imaginary of G=0 is zero
260 }
261}
Note: See TracBrowser for help on using the repository browser.