Changeset cc9c36

Timestamp:

Apr 22, 2008, 12:29:35 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

9bdd86

Parents:

87b8ed

Message:

CalculateChemicalShielding() renamed to CalculateMagneticMoment() and (small) changes to code.

The routine does already almost all that's wanted for moments, only prefactors needed change. The rest was just renaming [sS]igma to [mM]oment.

Location:

pcp/src

Files:

• perturbed.c (modified) (10 diffs)
• perturbed.h (modified) (1 diff)
• run.c (modified) (1 diff)

pcp/src/perturbed.c

@@ -3059 +3059 @@
   fftw_real *sigma_real[NDIM_NDIM];
   double sigma,Sigma;
-  int it, ion, in, dex, g, Index, i;
+  double x[NDIM];
+  int it, g, ion, in, dex, Index, i, j, d;
+  int n[NDIM];
   int *N = Lev0->Plan0.plan->N;
   //const double FFTfactor = 1.;///Lev0->MaxN;
@@ -3305 +3307 @@
 
 
-/** Calculates the chemical shielding tensor at the positions of the nuclei.
- * The chemical shielding tensor at position R is defined as the proportionality factor between the induced and
- * the externally applied field.
+/** Calculates the magnetic moment at the positions of the nuclei.
+ * The magnetic moment at position R is defined as
  * \f[
- *  \sigma_{ij} (R) = \frac{\delta B_j^{ind} (R)}{\delta B_i^{ext}} 
- *  = \frac{\mu_0}{4 \pi} \int d^3 r' \left ( \frac{r'-r}{| r' - r |^3} \times J_i (r') \right )_j 
+ *  m_{ij} (R) = \frac{1}{2} \int d^3 r' \left ( (r'-R) \times J_i (r') \right )_j 
  * \f]
  * One after another for each nuclear position is the tensor evaluated and the result printed
@@ -3316 +3316 @@
  * \param *P Problem at hand
  * \sa CalculateMagneticSusceptibility() - similar calculation, yet without translation to ion centers.
- * \warning This routine is out-dated due to being numerically unstable because of the singularity which is not
- * considered carefully, recommendend replacement is CalculateChemicalShieldingByReciprocalCurrentDensity().
  */
-void CalculateChemicalShielding(struct Problem *P)
+void CalculateMagneticMoment(struct Problem *P)
 {
   struct RunStruct *R = &P->R;
@@ -3326 +3324 @@
   struct Density *Dens0 = R->Lev0->Dens;
   struct Ions *I = &P->Ion;
-  double sigma[NDIM*NDIM],Sigma[NDIM*NDIM];
+  double moment[NDIM*NDIM],Moment[NDIM*NDIM];
   fftw_real *CurrentDensity[NDIM*NDIM];
   int it, ion, in, dex, i0, n[NDIM], n0, i;//, *NUp;
   double r[NDIM], fac[NDIM], X[NDIM];
   const double discrete_factor = Lat->Volume/Lev0->MaxN;
-  double eta, delta_sigma, S, A, iso;
+  double eta, delta_moment, S, A, iso;
   const int myPE =  P->Par.me_comm_ST_Psi;
   int N[NDIM];
@@ -3338 +3336 @@
   N[2] = Lev0->Plan0.plan->N[2];
   const int N0 = Lev0->Plan0.plan->local_nx;
-  FILE *SigmaFile;
-  char suffixsigma[255];
+  FILE *MomentFile;
+  char suffixmoment[255];
   time_t seconds;  
   time(&seconds); // get current time
   
+   if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) Integrating current density to evaluate magnetic moment\n", P->Par.me);
+ 
   // set pointers onto current density
   CurrentDensity[0] = (fftw_real *) Dens0->DensityArray[CurrentDensity0];
@@ -3357 +3357 @@
   for (it=0; it < I->Max_Types; it++) {  // integration over all types
     for (ion=0; ion < I->I[it].Max_IonsOfType; ion++) { // and each ion of type
-      if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) Shielding Tensor for Ion %i of element %s \\sigma_ij = \n",P->Par.me, ion, I->I[it].Name);
+      if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) Magnetic dipole moment Tensor for Ion %i of element %s \\moment_ij = \n",P->Par.me, ion, I->I[it].Name);
       for (in=0; in<NDIM; in++) {// index i of vector component in integrand
         for (dex=0; dex<NDIM; dex++) { // index j of B component derivation in current density tensor
-          sigma[in+dex*NDIM] = 0.;
+          moment[in+dex*NDIM] = 0.;
     
           for(n0=0;n0<N0;n0++) // do the integration over real space
@@ -3373 +3373 @@
                 i0 = n[2]+N[2]*(n[1]+N[1]*(n0));  // the index of current density must match LocalSizeR!
                 //z = MinImageConv(Lat,r, I->I[it].R[NDIM*ion],in);  // "in" always is missing third component in cross product
-                sigma[in+dex*NDIM] += (X[cross(in,0)] * CurrentDensity[dex*NDIM+cross(in,1)][i0] - X[cross(in,2)] * CurrentDensity[dex*NDIM+cross(in,3)][i0]);
+                moment[in+dex*NDIM] += (X[cross(in,0)] * CurrentDensity[dex*NDIM+cross(in,1)][i0] - X[cross(in,2)] * CurrentDensity[dex*NDIM+cross(in,3)][i0]);
                 //if (it == 0 && ion == 0) fprintf(stderr,"(%i) moment[%i][%i] += (%e * %e - %e * %e) = %e\n", P->Par.me, in, dex, x,CurrentDensity[dex*NDIM+cross(in,1)][i0],y,CurrentDensity[dex*NDIM+cross(in,3)][i0],moment[in+dex*NDIM]); 
               }
-          sigma[in+dex*NDIM] *= -mu0*discrete_factor/(4.*PI); // due to summation instead of integration
-          MPI_Allreduce ( &sigma[in+dex*NDIM], &Sigma[in+dex*NDIM], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);  // sum "LocalSize to TotalSize"
-          I->I[it].sigma[ion][in+dex*NDIM] = Sigma[in+dex*NDIM];          
-          if (P->Call.out[ValueOut]) fprintf(stderr," %e", Sigma[in+dex*NDIM]);
+          //moment[in+dex*NDIM] *= -mu0*discrete_factor/(4.*PI); // due to summation instead of integration
+          moment[in+dex*NDIM] *= 1./2.*discrete_factor; // due to summation instead of integration
+          MPI_Allreduce ( &moment[in+dex*NDIM], &Moment[in+dex*NDIM], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);  // sum "LocalSize to TotalSize"
+          I->I[it].moment[ion][in+dex*NDIM] = Moment[in+dex*NDIM];          
+          if (P->Call.out[ValueOut]) fprintf(stderr," %e", Moment[in+dex*NDIM]);
         }
         if (P->Call.out[ValueOut]) fprintf(stderr,"\n");
@@ -3386 +3387 @@
       for (in=0;in<NDIM;in++)
         for (dex=0;dex<NDIM;dex++)
-          gsl_matrix_complex_set(H,in,dex,gsl_complex_rect((Sigma[in+dex*NDIM]+Sigma[dex+in*NDIM])/2.,0));
+          gsl_matrix_complex_set(H,in,dex,gsl_complex_rect((Moment[in+dex*NDIM]+Moment[dex+in*NDIM])/2.,0));
       // output tensor to file
       if (P->Par.me == 0) {
-        sprintf(&suffixsigma[0], ".sigma_i%i_%s.L%i.csv", ion, I->I[it].Symbol, Lev0->LevelNo);
-        OpenFile(P, &SigmaFile, suffixsigma, "a", P->Call.out[ReadOut]);
-        fprintf(SigmaFile,"# chemical shielding tensor sigma[01,02,03,10,11,12,20,21,22], seed %i, config %s, run on %s", R->Seed, P->Files.default_path, ctime(&seconds));
-        fprintf(SigmaFile,"%lg\t", P->Lat.ECut/(Lat->LevelSizes[0]*Lat->LevelSizes[0]));
+        sprintf(&suffixmoment[0], ".moment_i%i_%s.L%i.csv", ion, I->I[it].Symbol, Lev0->LevelNo);
+        OpenFile(P, &MomentFile, suffixmoment, "a", P->Call.out[ReadOut]);
+        fprintf(MomentFile,"# magnetic tensor moment[01,02,03,10,11,12,20,21,22], seed %i, config %s, run on %s", R->Seed, P->Files.default_path, ctime(&seconds));
+        fprintf(MomentFile,"%lg\t", Lev0->ECut/4.);  // ECut is in Rydberg
         for (in=0;in<NDIM*NDIM;in++)
-          fprintf(SigmaFile,"%e\t", Sigma[in]);
-        fprintf(SigmaFile,"\n");
-        fclose(SigmaFile);
+          fprintf(MomentFile,"%e\t", Moment[in]);
+        fprintf(MomentFile,"\n");
+        fclose(MomentFile);
       }
-      // diagonalize sigma
+      // diagonalize moment
       gsl_vector *eval = gsl_vector_alloc(NDIM);
       gsl_eigen_herm_workspace *w = gsl_eigen_herm_alloc(NDIM);
@@ -3414 +3415 @@
       iso = 0;
       for (i=0;i<NDIM;i++) {
-        I->I[it].sigma[ion][i] = gsl_vector_get(eval,i);
-        iso += Sigma[i+i*NDIM]/3.;
+        I->I[it].moment[ion][i] = gsl_vector_get(eval,i);
+        iso += Moment[i+i*NDIM]/3.;
       } 
       eta = (gsl_vector_get(eval,1)-gsl_vector_get(eval,0))/(gsl_vector_get(eval,2)-iso);
-      delta_sigma = gsl_vector_get(eval,2) - 0.5*(gsl_vector_get(eval,0)+gsl_vector_get(eval,1));
-      S = (delta_sigma*delta_sigma)*(1+1./3.*eta*eta);
+      delta_moment = gsl_vector_get(eval,2) - 0.5*(gsl_vector_get(eval,0)+gsl_vector_get(eval,1));
+      S = (delta_moment*delta_moment)*(1+1./3.*eta*eta);
       A = 0.;
       for (i=0;i<NDIM;i++) {
         in = cross(i,0);
         dex = cross(i,1);
-        A += pow(-1,i)*pow(0.5*(Sigma[in+dex*NDIM]-Sigma[dex+in*NDIM]),2);
+        A += pow(-1,i)*pow(0.5*(Moment[in+dex*NDIM]-Moment[dex+in*NDIM]),2);
       }
       if (P->Call.out[ValueOut]) {
@@ -3431 +3432 @@
           fprintf(stderr,"\t%lg",gsl_vector_get(eval,i));
         fprintf(stderr,"\nshielding  : %e\n", iso);
-        fprintf(stderr,"anisotropy : %e\n", delta_sigma);
+        fprintf(stderr,"anisotropy : %e\n", delta_moment);
         fprintf(stderr,"asymmetry  : %e\n", eta);
         fprintf(stderr,"S          : %e\n", S);

pcp/src/perturbed.h

@@ -31 +31 @@
 void ApplyTotalHamiltonian(struct Problem *P, const fftw_complex *source, fftw_complex *dest, fftw_complex ***fnl, const double PsiFactor, const int first);
 void CalculateMagneticSusceptibility (struct Problem *P);
-void CalculateChemicalShielding(struct Problem *P);
+void CalculateMagneticMoment(struct Problem *P);
 void CalculateChemicalShieldingByReciprocalCurrentDensity(struct Problem *P);
 void CalculateChemicalShieldingbyDESolver(struct Problem *P);

pcp/src/run.c

@@ -1138 +1138 @@
           CalculateMagneticSusceptibility(P);
           debug(P,"Normal calculation of shielding over R-space");
-          CalculateChemicalShielding(P);
+          CalculateMagneticMoment(P);
           CalculateChemicalShieldingByReciprocalCurrentDensity(P);
           SpeedMeasure(P, CurrDensTime, StopTimeDo);