Commit JT 022120

- added message for only one precession/spin (+doc)
- added a per pair/spin class emag table

doc/src/fix_precession_spin.rst

@@ -82,6 +82,7 @@ function for the same parameters.
 
 .. image:: JPG/zeeman_langevin.jpg
    :align: center
+   :width: 600
 
 The temperature effects are accounted for by connecting the spin
 :math:`i` to a thermal bath using a Langevin thermostat (see
@@ -159,11 +160,15 @@ No information about this fix is written to :doc:`binary restart files <restart>
 Restrictions
 """"""""""""
 
-
 The *precession/spin* style is part of the SPIN package.  This style
 is only enabled if LAMMPS was built with this package, and if the
 atom\_style "spin" was declared.  See the :doc:`Build package <Build_package>` doc page for more info.
 
+The *precession/spin* style can only be declared once. If more
+than one precession type (for example combining an anisotropy and a Zeeman interactions)
+has to be declared, they have to be chained in the same command
+line (as shown in the examples above). 
+
 Related commands
 """"""""""""""""
 

src/SPIN/compute_spin.cpp

@@ -69,6 +69,22 @@ void ComputeSpin::init()
   hbar = force->hplanck/MY_2PI;
   kb = force->boltz;
 
+  // loop 1: obtain # of Pairs, and # of Pair/Spin styles
+
+  if (force->pair_match("spin",0,0)) {        // only one Pair/Spin style
+    pair = force->pair_match("spin",0,0);
+    npairs = pair->instance_total;
+    npairspin = 1;
+  } else if (force->pair_match("spin",0,1)) { // more than one Pair/Spin style
+    pair = force->pair_match("spin",0,1);
+    npairs = pair->instance_total;
+    for (int i = 0; i<npairs; i++) {
+      if (force->pair_match("spin",0,i)) {
+        npairspin ++;
+      }
+    }
+  }
+
   // init length of vector of ptrs to Pair/Spin styles
 
   if (npairspin > 0) {
@@ -166,7 +182,7 @@ void ComputeSpin::compute_vector()
 
         if (pair_spin_flag) {
           for (int k = 0; k < npairspin; k++) {
-            // spin_pairs[k]->compute_single_pair(i,fmi);
+            magenergy += spin_pairs[k]->emag[i];
           }
         }
 

src/SPIN/fix_precession_spin.cpp

@@ -197,6 +197,15 @@ void FixPrecessionSpin::init()
         error->all(FLERR,"Illegal precession/spin command");
   }
 
+  // check that fix precession/spin is only declared once
+
+  int iprec = 0;
+  for (int iforce = 0; iforce < modify->nfix; iforce++)
+    if (strstr(modify->fix[iforce]->style,"precession/spin")) iprec++;
+  if (iprec > 1)
+    error->all(FLERR,"precession/spin command can only be declared once");
+
+
   varflag = CONSTANT;
   if (magfieldstyle != CONSTANT) varflag = EQUAL;
 

src/SPIN/pair_spin.cpp

@@ -29,6 +29,7 @@
 #include "fix.h"
 #include "force.h"
 #include "math_const.h"
+#include "memory.h"
 #include "modify.h"
 #include "neighbor.h"
 #include "neigh_request.h"
@@ -98,4 +99,9 @@ void PairSpin::init_style()
   if (ifix >=0)
     lattice_flag = ((FixNVESpin *) modify->fix[ifix])->lattice_flag;
 
+  // test emag list storing mag energies
+  // init. size of energy stacking lists
+
+  nlocal_max = atom->nlocal;
+  memory->grow(emag,nlocal_max,"pair/spin:emag");
 }

src/SPIN/pair_spin.h

@@ -31,6 +31,10 @@ friend class FixNVESpin;
 
   virtual void compute(int, int) {}
   virtual void compute_single_pair(int, double *) {}
+  
+  // test emag list storing mag energies
+  int nlocal_max;                       // max value of nlocal (for size of lists)
+  double *emag;                         // energy list
 
  protected:
   double hbar;                          // Planck constant (eV.ps.rad-1)

src/SPIN/pair_spin_dipole_cut.cpp

@@ -64,6 +64,9 @@ PairSpinDipoleCut::~PairSpinDipoleCut()
     memory->destroy(setflag);
     memory->destroy(cut_spin_long);
     memory->destroy(cutsq);
+  
+    // test emag list storing mag energies
+    memory->destroy(emag);
   }
 }
 
@@ -185,6 +188,13 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // test emag list storing mag energies
+  // checking size of emag
+  if (nlocal_max < nlocal) {                    // grow emag lists if necessary
+    nlocal_max = nlocal;
+    memory->grow(emag,nlocal_max,"pair/spin:emag");
+  }
+
   // computation of the exchange interaction
   // loop over atoms and their neighbors
 

src/SPIN/pair_spin_dipole_long.cpp

@@ -69,6 +69,9 @@ PairSpinDipoleLong::~PairSpinDipoleLong()
     memory->destroy(setflag);
     memory->destroy(cut_spin_long);
     memory->destroy(cutsq);
+  
+    // test emag list storing mag energies
+    memory->destroy(emag);
   }
 }
 
@@ -212,6 +215,13 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // test emag list storing mag energies
+  // checking size of emag
+  if (nlocal_max < nlocal) {                    // grow emag lists if necessary
+    nlocal_max = nlocal;
+    memory->grow(emag,nlocal_max,"pair/spin:emag");
+  }
+
   pre1 = 2.0 * g_ewald / MY_PIS;
   pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS;
   pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS;
@@ -221,16 +231,20 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
 
   for (ii = 0; ii < inum; ii++) {
     i = ilist[ii];
+    itype = type[i];
+    
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
     xi[0] = x[i][0];
     xi[1] = x[i][1];
     xi[2] = x[i][2];
-    jlist = firstneigh[i];
-    jnum = numneigh[i];
     spi[0] = sp[i][0];
     spi[1] = sp[i][1];
     spi[2] = sp[i][2];
     spi[3] = sp[i][3];
-    itype = type[i];
+  
+    // test emag list storing mag energies
+    emag[i] = 0.0;
 
     for (jj = 0; jj < jnum; jj++) {
       j = jlist[jj];
@@ -294,9 +308,9 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
 
       if (eflag) {
         if (rsq <= local_cut2) {
-          evdwl -= spi[0]*fmi[0] + spi[1]*fmi[1] +
-            spi[2]*fmi[2];
+          evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
           evdwl *= 0.5*hbar;
+          emag[i] += evdwl;
         }
       } else evdwl = 0.0;
 

src/SPIN/pair_spin_dmi.cpp

@@ -53,6 +53,9 @@ PairSpinDmi::~PairSpinDmi()
     memory->destroy(vmech_dmy);
     memory->destroy(vmech_dmz);
     memory->destroy(cutsq);
+  
+    // test emag list storing mag energies
+    memory->destroy(emag);
   }
 }
 
@@ -191,6 +194,13 @@ void PairSpinDmi::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // test emag list storing mag energies
+  // checking size of emag
+  if (nlocal_max < nlocal) {                    // grow emag lists if necessary
+    nlocal_max = nlocal;
+    memory->grow(emag,nlocal_max,"pair/spin:emag");
+  }
+
   // dmi computation
   // loop over all atoms
 
@@ -206,7 +216,9 @@ void PairSpinDmi::compute(int eflag, int vflag)
     spi[0] = sp[i][0];
     spi[1] = sp[i][1];
     spi[2] = sp[i][2];
-
+  
+    // test emag list storing mag energies
+    emag[i] = 0.0;
 
     // loop on neighbors
 
@@ -260,6 +272,7 @@ void PairSpinDmi::compute(int eflag, int vflag)
       if (eflag) {
         evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
         evdwl *= 0.5*hbar;
+        emag[i] += evdwl;
       } else evdwl = 0.0;
 
       if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,

src/SPIN/pair_spin_exchange.cpp

@@ -50,6 +50,9 @@ PairSpinExchange::~PairSpinExchange()
     memory->destroy(J2);
     memory->destroy(J3);
     memory->destroy(cutsq); // to be implemented
+  
+    // test emag list storing mag energies
+    memory->destroy(emag);
   }
 }
 
@@ -176,6 +179,13 @@ void PairSpinExchange::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // test emag list storing mag energies
+  // checking size of emag
+  if (nlocal_max < nlocal) {                    // grow emag lists if necessary
+    nlocal_max = nlocal;
+    memory->grow(emag,nlocal_max,"pair/spin:emag");
+  }
+
   // computation of the exchange interaction
   // loop over atoms and their neighbors
 
@@ -191,6 +201,9 @@ void PairSpinExchange::compute(int eflag, int vflag)
     spi[0] = sp[i][0];
     spi[1] = sp[i][1];
     spi[2] = sp[i][2];
+  
+    // test emag list storing mag energies
+    emag[i] = 0.0;
 
     // loop on neighbors
 
@@ -243,6 +256,10 @@ void PairSpinExchange::compute(int eflag, int vflag)
       if (eflag) {
         evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
         evdwl *= 0.5*hbar;
+        // printf("test ex energy: %g \n",evdwl);
+        // evdwl = -0.5*compute_energy(i,j,rsq,spi,spj);
+        // printf("test ex energy: %g \n",evdwl);
+        emag[i] += evdwl;
         // evdwl *= hbar;
       } else evdwl = 0.0;
 
@@ -385,6 +402,29 @@ void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double ei
   fi[2] -= Jex_mech*eij[2];
 }
 
+/* ----------------------------------------------------------------------
+   compute energy of spin pair i and j
+------------------------------------------------------------------------- */
+
+// double PairSpinExchange::compute_energy(int i, int j, double rsq, double spi[3], double spj[3])
+// {
+//   int *type = atom->type;
+//   int itype, jtype;
+//   double Jex, ra;
+//   double energy = 0.0;
+//   itype = type[i];
+//   jtype = type[j];
+// 
+//   Jex = J1_mech[itype][jtype];
+//   ra = rsq/J3[itype][jtype]/J3[itype][jtype];
+//   Jex = 4.0*Jex*ra;
+//   Jex *= (1.0-J2[itype][jtype]*ra);
+//   Jex *= exp(-ra);
+// 
+//   energy = Jex*(spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]);
+//   return energy;
+// }
+
 /* ----------------------------------------------------------------------
    allocate all arrays
 ------------------------------------------------------------------------- */

src/SPIN/pair_spin_exchange.h

@@ -39,6 +39,8 @@ class PairSpinExchange : public PairSpin {
   void compute_exchange(int, int, double, double *, double *);
   void compute_exchange_mech(int, int, double, double *, double *, double *, double *);
 
+  // double compute_energy(int , int , double , double *, double *);
+
   void write_restart(FILE *);
   void read_restart(FILE *);
   void write_restart_settings(FILE *);

src/SPIN/pair_spin_magelec.cpp

@@ -51,6 +51,9 @@ PairSpinMagelec::~PairSpinMagelec()
     memory->destroy(v_mey);
     memory->destroy(v_mez);
     memory->destroy(cutsq); // to be deteled
+  
+    // test emag list storing mag energies
+    memory->destroy(emag);
   }
 }
 
@@ -185,6 +188,13 @@ void PairSpinMagelec::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // test emag list storing mag energies
+  // checking size of emag
+  if (nlocal_max < nlocal) {                    // grow emag lists if necessary
+    nlocal_max = nlocal;
+    memory->grow(emag,nlocal_max,"pair/spin:emag");
+  }
+
   // magneto-electric computation
   // loop over atoms and their neighbors
 
@@ -200,6 +210,9 @@ void PairSpinMagelec::compute(int eflag, int vflag)
     spi[0] = sp[i][0];
     spi[1] = sp[i][1];
     spi[2] = sp[i][2];
+  
+    // test emag list storing mag energies
+    emag[i] = 0.0;
 
     // loop on neighbors
 
@@ -252,6 +265,7 @@ void PairSpinMagelec::compute(int eflag, int vflag)
       if (eflag) {
         evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
         evdwl *= 0.5*hbar;
+        emag[i] += evdwl;
       } else evdwl = 0.0;
 
       if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,

src/SPIN/pair_spin_neel.cpp

@@ -54,6 +54,9 @@ PairSpinNeel::~PairSpinNeel()
     memory->destroy(q2);
     memory->destroy(q3);
     memory->destroy(cutsq); // to be deleted
+    
+    // test emag list storing mag energies
+    memory->destroy(emag);
   }
 }
 
@@ -190,6 +193,13 @@ void PairSpinNeel::compute(int eflag, int vflag)
   numneigh = list->numneigh;
   firstneigh = list->firstneigh;
 
+  // test emag list storing mag energies
+  // checking size of emag
+  if (nlocal_max < nlocal) {                    // grow emag lists if necessary
+    nlocal_max = nlocal;
+    memory->grow(emag,nlocal_max,"pair/spin:emag");
+  }
+
   // computation of the neel interaction
   // loop over atoms and their neighbors
 
@@ -206,6 +216,9 @@ void PairSpinNeel::compute(int eflag, int vflag)
     spi[1] = sp[i][1];
     spi[2] = sp[i][2];
 
+    // test emag list storing mag energies
+    emag[i] = 0.0;
+    
     // loop on neighbors
 
     for (jj = 0; jj < jnum; jj++) {
@@ -262,6 +275,7 @@ void PairSpinNeel::compute(int eflag, int vflag)
         // evdwl = (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
         evdwl = compute_neel_energy(i,j,rsq,eij,spi,spj);
         evdwl *= 0.5*hbar;
+        emag[i] += evdwl;
       } else evdwl = 0.0;
 
       if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,