Merge pull request #3461 from julient31/llg-correction

Correction of some LL features in the SPIN package

doc/src/Howto_spins.rst

@@ -30,9 +30,11 @@ can be coupled to another Langevin thermostat applied to the atoms
 using :doc:`fix langevin <fix_langevin>` in order to simulate
 thermostatted spin-lattice systems.
 
-The magnetic Gilbert damping can also be applied using :doc:`fix langevin/spin <fix_langevin_spin>`. It allows to either dissipate
-the thermal energy of the Langevin thermostat, or to perform a
-relaxation of the magnetic configuration toward an equilibrium state.
+The magnetic damping can also be applied
+using :doc:`fix langevin/spin <fix_langevin_spin>`.
+It allows to either dissipate the thermal energy of the Langevin
+thermostat, or to perform a relaxation of the magnetic configuration
+toward an equilibrium state.
 
 The command :doc:`fix setforce/spin <fix_setforce>` allows to set the
 components of the magnetic precession vectors (while erasing and
@@ -52,9 +54,11 @@ All the computed magnetic properties can be output by two main
 commands. The first one is :doc:`compute spin <compute_spin>`, that
 enables to evaluate magnetic averaged quantities, such as the total
 magnetization of the system along x, y, or z, the spin temperature, or
-the magnetic energy. The second command is :doc:`compute property/atom <compute_property_atom>`. It enables to output all the
-per atom magnetic quantities. Typically, the orientation of a given
-magnetic spin, or the magnetic force acting on this spin.
+the magnetic energy. The second command
+is :doc:`compute property/atom <compute_property_atom>`.
+It enables to output all the per atom magnetic quantities. Typically,
+the orientation of a given magnetic spin, or the magnetic force
+acting on this spin.
 
 ----------
 

doc/src/fix_langevin_spin.rst

@@ -40,7 +40,7 @@ the following stochastic differential equation:
   \times\left( \vec{\omega}_{i} \times\vec{s}_{i} \right) \right)
 
 with :math:`\lambda` the transverse damping, and :math:`\eta` a random vector.
-This equation is referred to as the stochastic Landau-Lifshitz-Gilbert (sLLG)
+This equation is referred to as the stochastic Landau-Lifshitz (sLL)
 equation.
 
 The components of :math:`\eta` are drawn from a Gaussian probability
@@ -49,7 +49,7 @@ the external thermostat T (in K in metal units).
 
 More details about this implementation are reported in :ref:`(Tranchida) <Tranchida2>`.
 
-Note: due to the form of the sLLG equation, this fix has to be defined just
+Note: due to the form of the sLL equation, this fix has to be defined just
 before the nve/spin fix (and after all other magnetic fixes).
 As an example:
 

doc/src/min_modify.rst

@@ -24,7 +24,7 @@ Syntax
          inf = max force component across all 3-vectors
          max = max force norm across all 3-vectors
        *alpha_damp* value = damping
-         damping = fictitious Gilbert damping for spin minimization (adim)
+         damping = fictitious magnetic damping for spin minimization (adim)
        *discrete_factor* value = factor
          factor = discretization factor for adaptive spin timestep (adim)
        *integrator* value = *eulerimplicit* or *verlet*
@@ -109,9 +109,9 @@ norm is replaced by the spin-torque norm.
 
 Keywords *alpha_damp* and *discrete_factor* only make sense when
 a :doc:`min_spin <min_spin>` command is declared.
-Keyword *alpha_damp* defines an analog of a magnetic Gilbert
-damping. It defines a relaxation rate toward an equilibrium for
-a given magnetic system.
+Keyword *alpha_damp* defines an analog of a magnetic damping.
+It defines a relaxation rate toward an equilibrium for a given
+magnetic system.
 Keyword *discrete_factor* defines a discretization factor for the
 adaptive timestep used in the *spin* minimization.
 See :doc:`min_spin <min_spin>` for more information about those

doc/src/min_spin.rst

@@ -39,7 +39,7 @@ timestep, according to:
 
    \frac{d \vec{s}_{i}}{dt} = \lambda\, \vec{s}_{i} \times\left( \vec{\omega}_{i} \times\vec{s}_{i} \right)
 
-with :math:`\lambda` a damping coefficient (similar to a Gilbert
+with :math:`\lambda` a damping coefficient (similar to a magnetic
 damping). :math:`\lambda` can be defined by setting the
 *alpha_damp* keyword with the :doc:`min_modify <min_modify>` command.
 

examples/SPIN/test_problems/README

@@ -44,7 +44,9 @@ directory.
   energy versus temperature. Comparison to the Langevin function
   results (computed by the python script).
   Note: This example is a reworked version of a test problem 
-  provided by Martin Kroger (ETHZ).  
+  provided by Martin Kroger (ETHZ). 
+  Note 2: Two versions of this test are deployed, one at low 
+  damping (0.01) and one at large damping (1.0).
 
 - validation_nve:
   simulates a small assembly of magnetic atoms (54). The atoms are 

examples/SPIN/test_problems/run_all.sh

@@ -12,13 +12,19 @@ cd validation_damped_precession/
 rm res_lammps.dat res_llg.dat
 cd ..
 
-# test 3: langevin, damping and Zeeman
-cd validation_langevin_precession/
+# test 3: langevin, damping and Zeeman, low damping (0.01)
+cd validation_langevin_precession_d0.01/
 ./run-test-prec.sh
 rm average_spin test-prec-spin.in res_lammps.dat res_langevin.dat
 cd ..
 
-# test 4: NVE run, test Etot preservation
+# test 4: langevin, damping and Zeeman, large damping (1.0)
+cd validation_langevin_precession_d1.0/
+./run-test-prec.sh
+rm average_spin test-prec-spin.in res_lammps.dat res_langevin.dat
+cd ..
+
+# test 5: NVE run, test Etot preservation
 cd validation_nve/
 ./run-test-nve.sh
 rm nve_spin_lattice.pdf res_lammps.dat

examples/SPIN/test_problems/validation_langevin_precession_d1.0/langevin.py

@@ -0,0 +1,22 @@
+#!/usr/bin/env python3
+
+import numpy as np, pylab, tkinter
+import matplotlib.pyplot as plt
+import mpmath as mp
+
+mub=5.78901e-5          # Bohr magneton (eV/T)
+kb=8.617333262145e-5    # Boltzman constant (eV/K)
+g=2.0                   # Lande factor (adim)
+
+Hz=10.0                # mag. field (T)
+
+#Definition of the Langevin function
+def func(t):
+    return mp.coth(g*mub*Hz/(kb*t))-1.0/(g*mub*Hz/(kb*t))
+
+npoints=200
+ti=0.01
+tf=20.0
+for i in range (0,npoints): 
+    temp=ti+i*(tf-ti)/npoints
+    print('%lf %lf %lf' % (temp,func(temp),-g*mub*Hz*func(temp)))

examples/SPIN/test_problems/validation_langevin_precession_d1.0/plot_precession.py

@@ -0,0 +1,37 @@
+#!/usr/bin/env python3
+
+#Program fitting the exchange interaction
+#Model curve: Bethe-Slater function
+import numpy as np, pylab, tkinter
+import matplotlib.pyplot as plt
+from scipy.optimize import curve_fit
+from decimal import *
+import sys, string, os
+
+
+argv = sys.argv
+if len(argv) != 3:
+  print("Syntax: ./plot_precession.py res_lammps.dat res_langevin.dat")
+  sys.exit()
+
+lammps_file = sys.argv[1]
+langevin_file = sys.argv[2]
+
+T_lmp,S_lmp,E_lmp = np.loadtxt(lammps_file, skiprows=0, usecols=(0,2,3),unpack=True)
+T_lan,S_lan,E_lan = np.loadtxt(langevin_file, skiprows=0, usecols=(0,1,2),unpack=True)
+
+plt.figure()
+plt.subplot(211)
+plt.ylabel('<Sz>')
+plt.plot(T_lmp, S_lmp, 'b-', label='LAMMPS')
+plt.plot(T_lan, S_lan, 'r--', label='Langevin')
+
+plt.subplot(212)
+plt.ylabel('E (in eV)')
+plt.plot(T_lmp, E_lmp, 'b-', label='LAMMPS')
+plt.plot(T_lan, E_lan, 'r--', label='Langevin')
+
+plt.xlabel('T (in K)')
+pylab.xlim([0,20.0])
+plt.legend()
+plt.show()

examples/SPIN/test_problems/validation_langevin_precession_d1.0/run-test-prec.sh

@@ -0,0 +1,33 @@
+#!/bin/bash
+
+tempi=0.0
+tempf=20.0
+
+rm res_*.dat
+
+# compute Lammps 
+N=20
+for (( i=0; i<$N; i++ ))
+do
+  temp="$(echo "$tempi+$i*($tempf-$tempi)/$N" | bc -l)"
+  sed s/temperature/${temp}/g test-prec-spin.template > \
+    test-prec-spin.in
+  
+  # test standard Lammps 
+  ./../../../../src/lmp_serial -in test-prec-spin.in 
+
+  # test spin/kk with Kokkos Lammps
+  # mpirun -np 1 ../../../../src/lmp_kokkos_mpi_only \
+  #   -k on -sf kk -in test-prec-spin.in
+  
+  Hz="$(tail -n 1 average_spin | awk -F " " '{print $3}')"
+  sz="$(tail -n 1 average_spin | awk -F " " '{print $5}')"
+  en="$(tail -n 1 average_spin | awk -F " " '{print $6}')"
+  echo $temp $Hz $sz $en >> res_lammps.dat
+done
+
+# compute Langevin
+python3 langevin.py > res_langevin.dat
+
+# plot results
+python3 plot_precession.py res_lammps.dat res_langevin.dat

examples/SPIN/test_problems/validation_langevin_precession_d1.0/test-prec-spin.template

@@ -0,0 +1,48 @@
+#LAMMPS in.run 
+ 
+units           metal
+atom_style      spin
+# atom_style      spin/kk
+atom_modify     map array
+boundary        p p p
+
+# read_data     singlespin.data
+
+lattice         sc 3.0
+region          box block 0.0 1.0 0.0 1.0 0.0 1.0
+create_box      1 box
+create_atoms    1 box
+
+mass            1 1.0
+set             type 1 spin 1.0 0.0 0.0 1.0
+
+# defines a pair/style for neighbor list, but do not use it
+pair_style      spin/exchange 4.0
+pair_coeff      * * exchange 1.0 0.0 0.0 1.0
+
+group bead      type 1  
+ 
+variable        H equal 10.0
+variable        Kan equal 0.0
+variable        Temperature equal temperature
+variable        RUN equal 1000000
+
+fix             1 all nve/spin lattice no
+fix             2 all precession/spin zeeman ${H} 0.0 0.0 1.0 anisotropy ${Kan} 0.0 0.0 1.0
+fix_modify      2 energy yes
+# fix             3 all langevin/spin ${Temperature} 0.01 12345
+fix             3 all langevin/spin ${Temperature} 1.0 12345
+
+compute         compute_spin all spin
+compute         outsp all property/atom spx spy spz sp
+compute         magsz all reduce ave c_outsp[3]
+
+thermo          50000
+thermo_style    custom step time temp vol pe c_compute_spin[5] etotal  
+
+variable        magnetic_energy equal c_compute_spin[5]
+
+fix             avespin all ave/time 1 ${RUN} ${RUN} v_Temperature v_H v_Kan c_magsz v_magnetic_energy file average_spin
+
+timestep        0.1
+run             ${RUN}

src/SPIN/compute_spin.cpp

@@ -43,7 +43,7 @@ using namespace MathConst;
 /* ---------------------------------------------------------------------- */
 
 ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) :
-  Compute(lmp, narg, arg), pair(nullptr), spin_pairs(nullptr)
+  Compute(lmp, narg, arg), lockprecessionspin(nullptr), pair(nullptr), spin_pairs(nullptr)
 {
   if ((narg != 3) && (narg != 4)) error->all(FLERR,"Illegal compute compute/spin command");
 
@@ -68,7 +68,8 @@ ComputeSpin::ComputeSpin(LAMMPS *lmp, int narg, char **arg) :
 ComputeSpin::~ComputeSpin()
 {
   memory->destroy(vector);
-  delete [] spin_pairs;
+  delete[] spin_pairs;
+  delete[] lockprecessionspin;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -96,7 +97,7 @@ void ComputeSpin::init()
     else npairs = hybrid->nstyles;
     for (int i = 0; i<npairs; i++) {
       if (force->pair_match("^spin",0,i)) {
-        npairspin ++;
+        npairspin++;
       }
     }
   }
@@ -135,14 +136,18 @@ void ComputeSpin::init()
     }
   }
 
-  // ptrs FixPrecessionSpin classes
+  // set ptrs for fix precession/spin styles
 
-  int iforce;
-  for (iforce = 0; iforce < modify->nfix; iforce++) {
-    if (utils::strmatch(modify->fix[iforce]->style,"^precession/spin")) {
-      precession_spin_flag = 1;
-      lockprecessionspin = dynamic_cast<FixPrecessionSpin *>(modify->fix[iforce]);
-    }
+  auto precfixes = modify->get_fix_by_style("^precession/spin");
+  nprecspin = precfixes.size();
+
+  if (nprecspin > 0) {
+    lockprecessionspin = new FixPrecessionSpin *[nprecspin];
+    precession_spin_flag = 1;
+
+    int i = 0;
+    for (auto &ifix : precfixes)
+      lockprecessionspin[i++] = dynamic_cast<FixPrecessionSpin *>(ifix);
   }
 }
 
@@ -191,7 +196,9 @@ void ComputeSpin::compute_vector()
         // update magnetic precession energies
 
         if (precession_spin_flag) {
-          magenergy += lockprecessionspin->emag[i];
+          for (int k = 0; k < nprecspin; k++) {
+            magenergy += lockprecessionspin[k]->emag[i];
+          }
         }
 
         // update magnetic pair interactions

src/SPIN/compute_spin.h

@@ -40,7 +40,8 @@ class ComputeSpin : public Compute {
 
   // pointers to magnetic fixes
 
-  class FixPrecessionSpin *lockprecessionspin;
+  int nprecspin;
+  class FixPrecessionSpin **lockprecessionspin;
 
   // pointers to magnetic pair styles
 

src/SPIN/fix_langevin_spin.cpp

@@ -109,7 +109,7 @@ void FixLangevinSpin::init()
   double hbar = force->hplanck/MY_2PI;  // eV/(rad.THz)
   double kb = force->boltz;             // eV/K
 
-  D = (alpha_t*gil_factor*kb*temp);
+  D = (alpha_t*(1.0+(alpha_t)*(alpha_t))*kb*temp);
   D /= (hbar*dts);
   sigma = sqrt(2.0*D);
 }

src/compute_property_atom.cpp

@@ -152,7 +152,6 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) :
         error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
       pack_choice[i] = &ComputePropertyAtom::pack_spx;
     } else if (strcmp(arg[iarg],"spy") == 0) {
-        error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
       if (!atom->sp_flag)
         error->all(FLERR,"Compute property/atom {} is not available", arg[iarg]);
       pack_choice[i] = &ComputePropertyAtom::pack_spy;