Changeset 34b70c for pcp/src/perturbed.c

Timestamp:

Apr 22, 2008, 12:18:49 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

87b8ed

Parents:

acd467

Message:

CalculateChemicalShieldingByReciprocalCurrentDensity(): NICS map written to file

fabs of sigma_real is taken (only positive values can be displayed), then B (sigma is basically B_Induced) is diagonalised and absolute values are delivered to OutputVis().

File:

• pcp/src/perturbed.c (modified) (2 diffs)

pcp/src/perturbed.c

@@ -3046 +3046 @@
   struct Lattice *Lat = &P->Lat;
   struct LatticeLevel *Lev0 = R->Lev0;
+  struct FileData *F = &P->Files;
   struct Ions *I = &P->Ion;
   struct Density *Dens0 = Lev0->Dens;
@@ -3263 +3264 @@
       sigma_real[in][i] = gsl_vector_get(eval,in);
   }
-  // Output of diagonalized magnetic field densities for each direction
-  for (i=0;i<NDIM;i++)
-    OutputVis(P, sigma_real[i]);    
+   
+  if (Lev0->LevelNo == 0) {
+    if (!P->Par.me && P->Call.out[NormalOut]) fprintf(stderr,"(%i)Output of NICS map ...\n", P->Par.me);
+    // Output of magnetic field densities for each direction
+    //for (i=0;i<NDIM*NDIM;i++)
+    //  OutputVis(P, sigma_real[i]);    
+    // Diagonalizing the tensor "field" B_ij [r]
+   if (P->Call.out[ValueOut])  fprintf(stderr,"(%i) Diagonalizing B_ij [r] ... \n", P->Par.me);
+    for (i=0; i< Dens0->LocalSizeR; i++) {
+      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
+          //fprintf(stderr,"(%i) Setting B_(%i,%i)[%i] ... \n", P->Par.me, in,dex,i);
+          gsl_matrix_complex_set(H,in,dex,gsl_complex_rect((sigma_real[in+dex*NDIM][i]+sigma_real[dex+in*NDIM][i])/2.,0.));
+        }
+      gsl_eigen_herm(H, eval, w);
+      gsl_sort_vector(eval);  // sort eigenvalues
+      for (in=0;in<NDIM;in++)
+        sigma_real[in][i] = gsl_vector_get(eval,in);
+    }
+  }
+
+  // now absolute the B values (as density scales them by log) and output
+  if (F->DoOutNICS) {
+    for (i=0; i< Dens0->LocalSizeR; i++)
+      for (in=0;in<NDIM;in++)
+        sigma_real[in][i] = fabs(sigma_real[in][i]);
+    // Output of diagonalized magnetic field densities for each direction
+    for (i=0;i<NDIM;i++)
+      OutputVis(P, sigma_real[i]);    
+  }
   for (i=0;i<NDIM*NDIM;i++)
     UnLockDensityArray(Dens0,i,real);  // sigma_imag/real free