Changeset 8786c3 for pcp/src

Timestamp:

Apr 22, 2008, 11:18:34 AM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

1d77026

Parents:

d2f1b1

Message:

CalculateChemicalShieldingByReciprocalCurrentDensity(): TestReciprocalCurrentDensity() outsourced

Part of CalculateChemicalShieldingByReciprocalCurrentDensity() was shifted to new function TestReciprocalCurrentDensity() that tested the reciprocal current density (J(G)=0)

Location:

pcp/src

Files:

• perturbed.c (modified) (2 diffs)
• perturbed.h (modified) (1 diff)

pcp/src/perturbed.c ¶

@@ -3062 +3062 @@
   CurrentDensity[8] = (fftw_real *) Dens0->DensityArray[CurrentDensity8];
 
-  if (P->Call.out[StepLeaderOut]) fprintf(stderr,"(%i) Checking J_{ij} (G=0) = 0 for each i,j ... \n",P->Par.me);
+  // inverse Fourier transform current densities
+  if (P->Call.out[LeaderOut]) fprintf(stderr,"(%i) Transforming and checking J_{ij} (G=0) = 0 for each i,j ... \n",P->Par.me);
   for (in=0;in<NDIM*NDIM;in++) {
     CalculateOneDensityC(Lat, R->LevS, Dens0, CurrentDensity[in], CurrentDensityC[in], factorDC);
-    tmp = sqrt(CurrentDensityC[in][0].re*CurrentDensityC[in][0].re+CurrentDensityC[in][0].im*CurrentDensityC[in][0].im);
-    if (GArray[0].GSq < MYEPSILON) {
-      if (in % NDIM == 0) fprintf(stderr,"(%i) ",P->Par.me); 
-      if (tmp > MYEPSILON) {
-        fprintf(stderr,"J_{%i,%i} = |%e + i%e| < %e ? (%e)\t", in / NDIM, in%NDIM, CurrentDensityC[in][0].re, CurrentDensityC[in][0].im, MYEPSILON, tmp - MYEPSILON);
-      } else { 
-        fprintf(stderr,"J_{%i,%i} ok\t", in / NDIM, in%NDIM);
-      }
-      if (in % NDIM == (NDIM-1)) fprintf(stderr,"\n");
-    }   
-  }
-
+    TestReciprocalCurrent(P, CurrentDensityC[in], GArray, in);
+  }
+
+  // linking pointers to the arrays
   for (in=0;in<NDIM*NDIM;in++) {
     LockDensityArray(Dens0,in,real); // Psi1R
@@ -3397 +3390 @@
 
   gsl_matrix_complex_free(H);
+}
+
+/** Test if G=0-component of reciprocal current is 0.
+ * In most cases we do not reach a numerical sensible zero as in MYEPSILON and remain satisfied as long
+ * as the integrated current density is very small (e.g. compared to single entries in the current density array)
+ * \param *P Problem at hand
+ * \param *CurrentC pointer to reciprocal current density
+ * \param *GArray pointer to array with G vectors
+ * \param in index of current component
+ * \sa TestCurrent() these two tests are equivalent and follow by fourier transformation
+ */
+void TestReciprocalCurrent(struct Problem *P, const fftw_complex *CurrentC, struct OneGData *GArray, int in)
+{
+  double tmp;
+  tmp = sqrt(CurrentC[0].re*CurrentC[0].re+CurrentC[0].im*CurrentC[0].im);
+  if ((P->Call.out[LeaderOut]) && (GArray[0].GSq < MYEPSILON)) {
+    if (in % NDIM == 0) fprintf(stderr,"(%i) ",P->Par.me); 
+    if (tmp > MYEPSILON) {
+      fprintf(stderr,"J_{%i,%i} = |%e + i%e| < %e ? (%e)\t", in / NDIM, in%NDIM, CurrentC[0].re, CurrentC[0].im, MYEPSILON, tmp - MYEPSILON);
+    } else { 
+      fprintf(stderr,"J_{%i,%i} ok\t", in / NDIM, in%NDIM);
+    }
+    if (in % NDIM == (NDIM-1)) fprintf(stderr,"\n");
+  }   
 }
 

pcp/src/perturbed.h ¶

@@ -49 +49 @@
 void TestSawtooth(struct Problem *P, const int index);
 void TestCurrent(struct Problem *P, const int index);
+void TestReciprocalCurrent(struct Problem *P, const fftw_complex *currentC, struct OneGData *GArray, int in);
 
 #endif /*PERTURBED_H_*/