Adding a fix setforce/spin to the SPIN package

doc/src/Howto_spins.txt

@@ -10,7 +10,7 @@ Documentation"_ld - "LAMMPS Commands"_lc :c
 Magnetic spins :h3
 
 The magnetic spin simulations are enabled by the SPIN package, whose
-implementation is detailed in "Tranchida"_#Tranchida7.
+implementation is detailed in "Tranchida"_#Tranchida.
 
 The model represents the simulation of atomic magnetic spins coupled
 to lattice vibrations. The dynamics of those magnetic spins can be used
@@ -36,13 +36,28 @@ A Langevin thermostat can be applied to those magnetic spins using
 "fix langevin/spin"_fix_langevin_spin.html. Typically, this thermostat
 can be coupled to another Langevin thermostat applied to the atoms
 using "fix langevin"_fix_langevin.html in order to simulate
-thermostatted spin-lattice system.
+thermostatted spin-lattice systems.
 
 The magnetic Gilbert damping can also be applied using "fix
 langevin/spin"_fix_langevin_spin.html. 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 "fix setforce/spin"_fix_setforce.html allows to set the 
+components of the magnetic precession vectors (while erasing and 
+replacing the previously computed magnetic precession vectors on 
+the atom). 
+This command can be used to freeze the magnetic moment of certain 
+atoms in the simulation by zeroing their precession vector. 
+
+The command "fix nve/spin"_fix_nve_spin.html can be used to
+perform a symplectic integration of the combined dynamics of spins 
+and atomic motions.
+
+The minimization style "min/spin"_min_spin.html can be applied
+to the spins to perform a minimization of the spin configuration. 
+
+
 All the computed magnetic properties can be output by two main
 commands. The first one is "compute spin"_compute_spin.html, that
 enables to evaluate magnetic averaged quantities, such as the total
@@ -54,6 +69,6 @@ magnetic spin, or the magnetic force acting on this spin.
 
 :line
 
-:link(Tranchida7)
+:link(Tranchida)
 [(Tranchida)] Tranchida, Plimpton, Thibaudeau and Thompson,
-arXiv preprint arXiv:1801.10233, (2018).
+Journal of Computational Physics, 372, 406-425, (2018).

doc/src/fix_langevin_spin.txt

@@ -99,4 +99,4 @@ integration fix (e.g. {fix nve/spin}).
 
 :link(Tranchida2)
 [(Tranchida)] Tranchida, Plimpton, Thibaudeau and Thompson,
-Journal of Computational Physics, (2018).
+Journal of Computational Physics, 372, 406-425, (2018).

doc/src/fix_nve_spin.txt

@@ -73,4 +73,4 @@ instead of "array" is also valid.
 
 :link(Tranchida1)
 [(Tranchida)] Tranchida, Plimpton, Thibaudeau and Thompson,
-Journal of Computational Physics, (2018).
+Journal of Computational Physics, 372, 406-425, (2018).

doc/src/fix_setforce.txt

@@ -8,6 +8,7 @@
 
 fix setforce command :h3
 fix setforce/kk command :h3
+fix setforce/spin command :h3
 
 [Syntax:]
 
@@ -27,6 +28,7 @@ keyword = {region}  :l
 
 fix freeze indenter setforce 0.0 0.0 0.0
 fix 2 edge setforce NULL 0.0 0.0
+fix 1 edge setforce/spin 0.0 0.0 0.0
 fix 2 edge setforce NULL 0.0 v_oscillate :pre
 
 [Description:]
@@ -65,6 +67,19 @@ to it.
 
 :line
 
+Style {spin} suffix sets the components of the magnetic precession 
+vectors instead of the mechanical forces. This also erases all 
+previously computed magnetic precession vectors on the atom, though 
+additional magnetic fixes could add new forces.
+
+This command can be used to freeze the magnetic moment of certain 
+atoms in the simulation by zeroing their precession vector. 
+
+All options defined above remain valid, they just apply to the magnetic 
+precession vectors instead of the forces.
+
+:line
+
 Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
 functionally the same as the corresponding style without the suffix.
 They have been optimized to run faster, depending on your available
@@ -117,7 +132,10 @@ forces to any value besides zero when performing a minimization.  Use
 the "fix addforce"_fix_addforce.html command if you want to apply a
 non-zero force to atoms during a minimization.
 
-[Restrictions:] none
+[Restrictions:] 
+
+The fix {setforce/spin} only makes sense when LAMMPS was built with the
+SPIN package.
 
 [Related commands:]
 

doc/src/pair_spin_dmi.txt

@@ -88,4 +88,4 @@ package"_Build_package.html doc page for more info.
 Physical Review B, 88(18), 184422. (2013).
 :link(Tranchida5)
 [(Tranchida)] Tranchida, Plimpton, Thibaudeau and Thompson,
-Journal of Computational Physics, (2018).
+Journal of Computational Physics, 372, 406-425, (2018).

doc/src/pair_spin_exchange.txt

@@ -95,4 +95,4 @@ package"_Build_package.html doc page for more info.
 
 :link(Tranchida3)
 [(Tranchida)] Tranchida, Plimpton, Thibaudeau and Thompson,
-Journal of Computational Physics, (2018).
+Journal of Computational Physics, 372, 406-425, (2018).

doc/src/pair_spin_magelec.txt

@@ -70,4 +70,4 @@ package"_Build_package.html doc page for more info.
 
 :link(Tranchida4)
 [(Tranchida)] Tranchida, Plimpton, Thibaudeau, and Thompson,
-Journal of Computational Physics, (2018).
+Journal of Computational Physics, 372, 406-425, (2018).

doc/src/pair_spin_neel.txt

@@ -80,4 +80,4 @@ package"_Build_package.html doc page for more info.
 
 :link(Tranchida6)
 [(Tranchida)] Tranchida, Plimpton, Thibaudeau and Thompson,
-Journal of Computational Physics, (2018).
+Journal of Computational Physics, 372, 406-425, (2018).

examples/SPIN/setforce_spin/in.spinmin.setforce

@@ -0,0 +1,59 @@
+
+units 		metal
+dimension 	3
+boundary 	f f f
+atom_style 	spin
+
+# necessary for the serial algorithm (sametag)
+atom_modify 	map array 
+
+lattice 	sc 3.0
+region 		box block 0.0 10.0 0.0 10.0 0.0 1.0
+create_box 	2 box
+region		reg1 block 0.0 10.0 	0.0 5.0 	0.0 1.0
+region		reg2 block 0.0 10.0 	6.0 10.0 	0.0 1.0
+create_atoms 	1 region reg1 
+create_atoms 	2 region reg2
+
+# setting mass, mag. moments, and interactions for bcc iron
+
+mass		1 55.845
+mass		2 55.845
+set 		region reg1 spin 2.2 0.0 0.0 1.0
+set 		region reg2 spin/random 31 2.2
+
+group 		fixed_spin region reg1
+
+pair_style 	hybrid/overlay spin/exchange 3.1 spin/dmi 3.1
+pair_coeff 	* * spin/exchange exchange 3.1 -0.01593 0.06626915552 1.211
+pair_coeff	* * spin/dmi dmi 3.1 0.12e-03 0.0 0.0 1.0
+
+neighbor 	0.1 bin
+neigh_modify 	every 10 check yes delay 20
+
+fix 		1 all precession/spin zeeman 0.0 0.0 0.0 1.0 anisotropy 5e-05 0.0 0.0 1.0 
+fix_modify	1 energy yes
+fix 		2 fixed_spin setforce/spin 0.0 0.0 0.0
+fix 		3 all langevin/spin 0.0 0.1 21
+fix		4 all nve/spin lattice no 
+
+timestep	0.0001
+
+compute 	out_mag    all spin
+variable 	magx      equal c_out_mag[1]
+variable 	magy      equal c_out_mag[2]
+variable 	magz      equal c_out_mag[3]
+variable 	magnorm   equal c_out_mag[4]
+variable 	emag      equal c_out_mag[5]
+variable 	tmag      equal c_out_mag[6]
+
+thermo          1000
+thermo_style    custom step time v_magx v_magz v_magnorm v_tmag etotal
+thermo_modify   format float %20.15g
+
+compute 	outsp all property/atom spx spy spz sp fmx fmy fmz
+dump 		1 all custom 1000 dump.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3] c_outsp[5] c_outsp[6] c_outsp[7]
+
+min_style	spin
+min_modify 	alpha_damp 1.0 discrete_factor 20.0
+minimize        1.0e-16 1.0e-16 50000 1000

src/SPIN/fix_nve_spin.cpp

@@ -31,9 +31,10 @@
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
-#include "fix_precession_spin.h"
 #include "fix_nve_spin.h"
+#include "fix_precession_spin.h"
 #include "fix_langevin_spin.h"
+#include "fix_setforce_spin.h"
 #include "force.h"
 #include "math_vector.h"
 #include "math_extra.h"
@@ -88,7 +89,6 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) :
   npairs = 0;
   npairspin = 0;
 
-
   // checking if map array or hash is defined
 
   if (atom->map_style == 0)
@@ -132,6 +132,7 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) :
   precession_spin_flag = 0;
   maglangevin_flag = 0;
   tdamp_flag = temp_flag = 0;
+  setforce_spin_flag = 0;
 
 }
 
@@ -238,6 +239,15 @@ void FixNVESpin::init()
    if (locklangevinspin->temp_flag == 1) temp_flag = 1;
   }
 
+  // ptrs FixSetForceSpin classes
+
+  for (iforce = 0; iforce < modify->nfix; iforce++) {
+    if (strstr(modify->fix[iforce]->style,"setforce/spin")) {
+      setforce_spin_flag = 1;
+      locksetforcespin = (FixSetForceSpin *) modify->fix[iforce];
+    }
+  }
+
   // setting the sector variables/lists
 
   nsectors = 0;
@@ -250,13 +260,10 @@ void FixNVESpin::init()
   // init. size of stacking lists (sectoring)
 
   nlocal_max = atom->nlocal;
-  stack_head = memory->grow(stack_head,nsectors,"NVE/spin:stack_head");
-  stack_foot = memory->grow(stack_foot,nsectors,"NVE/spin:stack_foot");
-  backward_stacks = memory->grow(backward_stacks,nlocal_max,"NVE/spin:backward_stacks");
-  forward_stacks = memory->grow(forward_stacks,nlocal_max,"NVE/spin:forward_stacks");
-  if (nlocal_max == 0)
-    error->all(FLERR,"Incorrect value of nlocal_max");
-
+  memory->grow(stack_head,nsectors,"NVE/spin:stack_head");
+  memory->grow(stack_foot,nsectors,"NVE/spin:stack_foot");
+  memory->grow(backward_stacks,nlocal_max,"NVE/spin:backward_stacks");
+  memory->grow(forward_stacks,nlocal_max,"NVE/spin:forward_stacks");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -389,8 +396,8 @@ void FixNVESpin::pre_neighbor()
 
   if (nlocal_max < nlocal) {                    // grow linked lists if necessary
     nlocal_max = nlocal;
-    backward_stacks = memory->grow(backward_stacks,nlocal_max,"NVE/spin:backward_stacks");
-    forward_stacks = memory->grow(forward_stacks,nlocal_max,"NVE/spin:forward_stacks");
+    memory->grow(backward_stacks,nlocal_max,"NVE/spin:backward_stacks");
+    memory->grow(forward_stacks,nlocal_max,"NVE/spin:forward_stacks");
   }
 
   for (int j = 0; j < nsectors; j++) {
@@ -462,6 +469,12 @@ void FixNVESpin::ComputeInteractionsSpin(int i)
     }
   }
 
+  // update setforce of magnetic interactions
+
+  if (setforce_spin_flag) {
+    locksetforcespin->single_setforce_spin(i,fmi);
+  }
+
   // replace the magnetic force fm[i] by its new value fmi
 
   fm[i][0] = fmi[0];
@@ -569,33 +582,33 @@ void FixNVESpin::AdvanceSingleSpin(int i)
   energy = (sp[i][0]*fm[i][0])+(sp[i][1]*fm[i][1])+(sp[i][2]*fm[i][2]);
   dts2 = dts*dts;
 
-  cp[0] = fm[i][1]*sp[i][2]-fm[i][2]*sp[i][1];
-  cp[1] = fm[i][2]*sp[i][0]-fm[i][0]*sp[i][2];
-  cp[2] = fm[i][0]*sp[i][1]-fm[i][1]*sp[i][0];
+  cp[0] = fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1];
+  cp[1] = fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2];
+  cp[2] = fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0];
 
-  g[0] = sp[i][0]+cp[0]*dts;
-  g[1] = sp[i][1]+cp[1]*dts;
-  g[2] = sp[i][2]+cp[2]*dts;
+  g[0] = sp[i][0] + cp[0]*dts;
+  g[1] = sp[i][1] + cp[1]*dts;
+  g[2] = sp[i][2] + cp[2]*dts;
 
-  g[0] += (fm[i][0]*energy-0.5*sp[i][0]*fm2)*0.5*dts2;
-  g[1] += (fm[i][1]*energy-0.5*sp[i][1]*fm2)*0.5*dts2;
-  g[2] += (fm[i][2]*energy-0.5*sp[i][2]*fm2)*0.5*dts2;
+  g[0] += (fm[i][0]*energy - 0.5*sp[i][0]*fm2)*0.5*dts2;
+  g[1] += (fm[i][1]*energy - 0.5*sp[i][1]*fm2)*0.5*dts2;
+  g[2] += (fm[i][2]*energy - 0.5*sp[i][2]*fm2)*0.5*dts2;
 
-  g[0] /= (1+0.25*fm2*dts2);
-  g[1] /= (1+0.25*fm2*dts2);
-  g[2] /= (1+0.25*fm2*dts2);
+  g[0] /= (1.0 + 0.25*fm2*dts2);
+  g[1] /= (1.0 + 0.25*fm2*dts2);
+  g[2] /= (1.0 + 0.25*fm2*dts2);
 
   sp[i][0] = g[0];
   sp[i][1] = g[1];
   sp[i][2] = g[2];
 
   // renormalization (check if necessary)
 
-  msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
-  scale = 1.0/sqrt(msq);
-  sp[i][0] *= scale;
-  sp[i][1] *= scale;
-  sp[i][2] *= scale;
+  //msq = g[0]*g[0] + g[1]*g[1] + g[2]*g[2];
+  //scale = 1.0/sqrt(msq);
+  //sp[i][0] *= scale;
+  //sp[i][1] *= scale;
+  //sp[i][2] *= scale;
 
   // comm. sp[i] to atoms with same tag (for serial algo)
 

src/SPIN/fix_nve_spin.h

@@ -58,14 +58,16 @@ friend class PairSpin;
   int nlocal_max;			// max value of nlocal (for lists size)
 
   int pair_spin_flag;			// magnetic pair flags
-  int precession_spin_flag;			// magnetic precession flags
+  int precession_spin_flag;		// magnetic precession flags
   int maglangevin_flag;			// magnetic langevin flags
   int tdamp_flag, temp_flag;
+  int setforce_spin_flag;
 
   // pointers to magnetic fixes
 
   class FixPrecessionSpin *lockprecessionspin;
   class FixLangevinSpin *locklangevinspin;
+  class FixSetForceSpin *locksetforcespin;
 
   // pointers to magnetic pair styles