Commit JT 052318

src/SPIN/compute_spin.cpp

@@ -132,6 +132,7 @@ void ComputeSpin::compute_vector()
   magtot[2] *= scale;
   magtot[3] = sqrt((magtot[0]*magtot[0])+(magtot[1]*magtot[1])+(magtot[2]*magtot[2]));
   spintemperature = hbar*tempnumtot;    
+  //spintemperature /= (2.0*kb*tempdenomtot);
   spintemperature /= (kb*tempdenomtot);
  
   vector[0] = magtot[0];

src/SPIN/fix_langevin_spin.cpp

@@ -134,8 +134,10 @@ void FixLangevinSpin::init()
   double hbar = force->hplanck/MY_2PI;	// eV/(rad.THz)
   double kb = force->boltz;		// eV/K
   D = (MY_2PI*alpha_t*gil_factor*kb*temp);
+  //D = (alpha_t*gil_factor*kb*temp);
   D /= (hbar*dts);
   sigma = sqrt(2.0*D);
+  //sigma = sqrt(D);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -171,9 +173,14 @@ void FixLangevinSpin::add_tdamping(double spi[3], double fmi[3])
 void FixLangevinSpin::add_temperature(double fmi[3]) 
 {
 
-  double rx = sigma*(-1.0+2.0*random->uniform());
-  double ry = sigma*(-1.0+2.0*random->uniform());
-  double rz = sigma*(-1.0+2.0*random->uniform());
+  double rx = sigma*(2.0*random->uniform() - 1.0);
+  double ry = sigma*(2.0*random->uniform() - 1.0);
+  double rz = sigma*(2.0*random->uniform() - 1.0);
+  //printf("test rd : %g \n",2.0*random->uniform() - 1.0);
+  //printf("test gaussian : %g \n", random->gaussian());
+  //double rx = sigma*(random->gaussian());
+  //double ry = sigma*(random->gaussian());
+  //double rz = sigma*(random->gaussian());
 
   // adding the random field 
 

src/SPIN/fix_nve_spin.cpp

@@ -164,6 +164,7 @@ void FixNVESpin::init()
   dtv = update->dt;
   dtf = 0.5 * update->dt * force->ftm2v;
   dts = 0.25 * update->dt;
+  npairs = npairspin = 0;
 
   // set ptrs on Pair/Spin styles
 

src/SPIN/pair_spin_exchange.cpp

@@ -173,6 +173,8 @@ void PairSpinExchange::init_style()
     }
   }
 
+  printf("test lattice flag: %d \n",lattice_flag);
+
 }
 
 /* ----------------------------------------------------------------------

src/SPIN/pair_spin_neel.cpp

@@ -285,7 +285,7 @@ void PairSpinNeel::compute(int eflag, int vflag)
 
       local_cut2 = cut_spin_neel[itype][jtype]*cut_spin_neel[itype][jtype];
 
-      // compute magnetic and mechanical components of soc_neel
+      // compute magnetic and mechanical components of neel
 
       if (rsq <= local_cut2) {
 	compute_neel(i,j,rsq,eij,fmi,spi,spj);
@@ -465,55 +465,55 @@ void PairSpinNeel::compute_neel_mech(int i, int j, double rsq, double eij[3], do
   itype = type[i];
   jtype = type[j];
   
-    double g_mech, gij, dgij;
-    double q_mech, q1ij, dq1ij;
-    double q2ij, dq2ij;
-    double pdx, pdy, pdz;
-    double pq1x, pq1y, pq1z;
-    double pq2x, pq2y, pq2z;
-    double ra, rr, drij, ig3, iq3;
-
-    drij = sqrt(rsq);
-  
-    double scalar_si_sj = spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2];
-    double scalar_eij_si = eij[0]*spi[0]+eij[1]*spi[1]+eij[2]*spi[2];
-    double scalar_eij_sj = eij[0]*spj[0]+eij[1]*spj[1]+eij[2]*spj[2];
-  
-    // pseudo-dipolar component
-    
-    g_mech = g1_mech[itype][jtype];        
-    ig3 = 1.0/(g3[itype][jtype]*g3[itype][jtype]);
-  
-    ra = rsq*ig3; 
-    rr = drij*ig3;
-  
-    gij = 4.0*g_mech*ra;
-    gij *= (1.0-g2[itype][jtype]*ra);
-    gij *= exp(-ra);
-  
-    dgij = 1.0-ra-g2[itype][jtype]*ra*(2.0-ra);
-    dgij *= 8.0*g_mech*rr*exp(-ra);
-  
-    double pdt1 = (dgij-2.0*gij/drij)*scalar_eij_si*scalar_eij_sj;
-    pdt1 -= scalar_si_sj*dgij/3.0;
-    double pdt2 = scalar_eij_sj*gij/drij;
-    double pdt3 = scalar_eij_si*gij/drij;
-    pdx = -(pdt1*eij[0] + pdt2*spi[0] + pdt3*spj[0]);
-    pdy = -(pdt1*eij[1] + pdt2*spi[1] + pdt3*spj[1]);
-    pdz = -(pdt1*eij[2] + pdt2*spi[2] + pdt3*spj[2]);
-  
-    // pseudo-quadrupolar component
-    
-    q_mech = q1_mech[itype][jtype];        
-    iq3 = 1.0/(q3[itype][jtype]*q3[itype][jtype]);
+  double g_mech, gij, dgij;
+  double q_mech, q1ij, dq1ij;
+  double q2ij, dq2ij;
+  double pdx, pdy, pdz;
+  double pq1x, pq1y, pq1z;
+  double pq2x, pq2y, pq2z;
+  double ra, rr, drij, ig3, iq3;
+
+  drij = sqrt(rsq);
+
+  double scalar_si_sj = spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2];
+  double scalar_eij_si = eij[0]*spi[0]+eij[1]*spi[1]+eij[2]*spi[2];
+  double scalar_eij_sj = eij[0]*spj[0]+eij[1]*spj[1]+eij[2]*spj[2];
+
+  // pseudo-dipolar component
   
-    ra = rsq*iq3; 
-    rr = drij*iq3;
+  g_mech = g1_mech[itype][jtype];        
+  ig3 = 1.0/(g3[itype][jtype]*g3[itype][jtype]);
+
+  ra = rsq*ig3; 
+  rr = drij*ig3;
+
+  gij = 4.0*g_mech*ra;
+  gij *= (1.0-g2[itype][jtype]*ra);
+  gij *= exp(-ra);
+
+  dgij = 1.0-ra-g2[itype][jtype]*ra*(2.0-ra);
+  dgij *= 8.0*g_mech*rr*exp(-ra);
+
+  double pdt1 = (dgij-2.0*gij/drij)*scalar_eij_si*scalar_eij_sj;
+  pdt1 -= scalar_si_sj*dgij/3.0;
+  double pdt2 = scalar_eij_sj*gij/drij;
+  double pdt3 = scalar_eij_si*gij/drij;
+  pdx = -(pdt1*eij[0] + pdt2*spi[0] + pdt3*spj[0]);
+  pdy = -(pdt1*eij[1] + pdt2*spi[1] + pdt3*spj[1]);
+  pdz = -(pdt1*eij[2] + pdt2*spi[2] + pdt3*spj[2]);
+
+  // pseudo-quadrupolar component
   
-    q1ij = 4.0*q_mech*ra;
-    q1ij *= (1.0-q2[itype][jtype]*ra);
-    q1ij *= exp(-ra);
-    q2ij = -2.0*q1ij/9.0;
+  q_mech = q1_mech[itype][jtype];        
+  iq3 = 1.0/(q3[itype][jtype]*q3[itype][jtype]);
+
+  ra = rsq*iq3; 
+  rr = drij*iq3;
+
+  q1ij = 4.0*q_mech*ra;
+  q1ij *= (1.0-q2[itype][jtype]*ra);
+  q1ij *= exp(-ra);
+  q2ij = -2.0*q1ij/9.0;
   
   dq1ij = 1.0-ra-q2[itype][jtype]*ra*(2.0-ra);
   dq1ij *= 8.0*q_mech*rr*exp(-ra);