Changeset 87b8ed for pcp

Timestamp:

Apr 22, 2008, 12:19:38 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

cc9c36

Parents:

34b70c

Message:

CalculateChemicalShieldingByReciprocalCurrentDensity(): "Some" untested magnetic forces are calculated

File:

• pcp/src/perturbed.c (modified) (1 diff)

pcp/src/perturbed.c

@@ -3248 +3248 @@
   }
   
-  // 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]
-  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));
+  if (R->MaxOuterStep > 0) { // if we do MD, calculate magnetic force with undiagonalised B fields  
+    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
+        // Finally use the magnetic moment in order to calculate the magnetic force
+        RMat33Vec3(x, Lat->ReciBasis, &(I->I[it].R[NDIM*ion]));
+        for (d=0;d<NDIM;d++)
+          n[d] = (int)(x[d]/(2.*PI)*(double)N[d]);  // round to next nearest mesh point
+//          n[d] = (int)(I->I[it].R[NDIM*ion+d]/Lat->RealBasisQ[d]*(double)N[d]);
+        for (d=0;d<NDIM;d++) { // index of induced magnetic field
+          I->I[it].FMagnetic[d+ion*NDIM] = 0.;
+          for (j=0;j<NDIM;j++) {// we to sum over all external field components
+            //fprintf(stderr,"(%i) Calculating magnetic force component %i over field component %i of ion (type %i, nr %i)\n", P->Par.me, d, j, it, ion);
+            I->I[it].FMagnetic[d+ion*NDIM] += - I->I[it].moment[ion][d] * FirstDiscreteDerivative(P, Lev0, sigma_real[d+NDIM*j], n, d, P->Par.me_comm_ST_Psi)*P->R.BField[j];
+          }
+        }
       }
-    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);
+    }
   }