Changeset f72dc1 for pcp/src

Timestamp:

Apr 22, 2008, 9:02:23 AM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

0f2c43

Parents:

ce81a3a

Message:

CalculateIonLocalForce(): added XC force

Force from XC potential was missing so far, is now included if ion's are CoreCorrected. In order to do that CalculateXCPotentialNoRT() has to be called beforehand.

File:

• pcp/src/pseudo.c (modified) (4 diffs)

pcp/src/pseudo.c

@@ -1068 +1068 @@
 }
 
-/** Calculates the local force acting on ion.
+/** Calculates the local force acting on each ion.
  * \param *P Problem at hand
  */
@@ -1080 +1080 @@
   struct RunStruct *R = &P->R;
   struct LatticeLevel *Lev0 = R->Lev0;
+  struct LatticeLevel *LevS = R->LevS;
   struct Density *Dens = Lev0->Dens;
+  struct  fft_plan_3d *plan = L->plan; 
+  fftw_complex *work = Dens->DensityCArray[TempDensity];
+  fftw_complex *HGC = Dens->DensityCArray[HGDensity];
+  fftw_real *HGCR = (fftw_real *)HGC;
   double ForceFac = L->Volume;
   double force, *dsum;
   fftw_complex cv;
+  
+  // precalc V_XC(r) and fft
+  //debug(P,"precalc V_XC\n");
+  SetArrayToDouble0((double *)HGCR, Dens->TotalSize*2);
+  CalculateXCPotentialNoRT(P, HGCR);
+  fft_3d_real_to_complex(plan, LevS->LevelNo, FFTNF1, HGC, work);
+  //debug(P,"precalc V_XC finished\n");
   for (is=0; is < I->Max_Types; is++) {
     SetArrayToDouble0(I->FTemp, NDIM*I->I[is].Max_IonsOfType);
@@ -1096 +1108 @@
                 c_re(cv) = c_re(PP->VCoulombc[g2])*PP->FacGauss[is][g2]+PP->phi_ps_loc[is][g2]*c_re(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
                 c_im(cv) = c_im(PP->VCoulombc[g2])*PP->FacGauss[is][g2]-PP->phi_ps_loc[is][g2]*c_im(Dens->DensityCArray[TotalDensity][Index])*R->HGcFactor;
+        if (I->I[is].corecorr == CoreCorrected) {
+          // this summand was missing before, is not thoroughly tested
+          //debug(P,"V_XC summand in CalculateIonLocalForce()\n");
+          c_re(cv) +=  (PP->formfCore[is][g2]*c_re(HGC[Index]))*R->HGcFactor;
+          c_im(cv) += -(PP->formfCore[is][g2]*c_im(HGC[Index]))*R->HGcFactor;
+          //fprintf(stderr, "(%i) FormfCore %lg\tHGC %lg+i%lg\n",P->Par.me, PP->formfCore[is][g2],c_re(HGC[Index]),c_im(HGC[Index]));
+        }
                 force = c_re(PP->ei[is][ia][g2])*c_im(cv)+c_im(PP->ei[is][ia][g2])*c_re(cv);
-                dsum[0] += 2.*(-G[0]*force);
-                dsum[1] += 2.*(-G[1]*force);
-                dsum[2] += 2.*(-G[2]*force);
+                dsum[0] += 2.*(-G[0]*force); //2.* ... why? Nowhere from derivative!
+                dsum[1] += 2.*(-G[1]*force); //2.*
+                dsum[2] += 2.*(-G[2]*force); //2.*
       }
     }
@@ -1105 +1124 @@
     for (ia=0;ia< I->I[is].Max_IonsOfType; ia++) 
       for (i=0; i<NDIM; i++)
-                I->I[is].FIonL[i+NDIM*ia] *= ForceFac;
+                                I->I[is].FIonL[i+NDIM*ia] *= ForceFac;
   }
 }