Commit JT 111519

- add README file to the benchmark examples repo
- removed comments from src/SPIN files

examples/SPIN/README

@@ -1,6 +1,10 @@
 This directory contains examples and applications of the SPIN package
 =====================================================================
 
+- the benchmark directory provides comparison between LAMMPS
+  results and a series of simple test problems (coded as python
+  scripts).
+
 - the iron, cobalt_hcp, cobalt_fcc and nickel directories provide 
   examples of spin-lattice calculations.
 
@@ -16,5 +20,5 @@ VMD. See the vmd repository for help vizualizing results with VMD.
 
 ** Note, the aim of this repository is mainly to provide users with
 examples. Better values and tuning of the magnetic and mechanical 
-interactions can be achieved for more accurate materials 
+interactions can (have to...) be achieved for more accurate materials 
 simulations. **

examples/SPIN/benchmark/README

@@ -0,0 +1,47 @@
+** The objective of the benchmark examples in this directory 
+   is the following twofold: 
+- verify the implementation of the LAMMPS' SPIN package by 
+  comparing its results to well-known analytic results, or 
+  to simple test problems,
+- provide users with simple comparisons, allowing them to 
+  better understand what is implemented in the code.
+
+The LAMMPS input file (bench-*) can be modified, as well as the
+associated python script, in order to try different comparisons.
+
+All scripts can be run by executing the shell script from its
+directory. Example:
+./run-bench-exchange.sh from the benchmarck_damped_exchange/
+directory. 
+
+** Below a brief description of the different benchmark 
+   problems:
+
+- benchmarck_damped_precession:
+  simulates the damped precession of a single spin in a magnetic
+  field. 
+  Run as: ./run-bench-prec.sh
+  Output: x, y and z components of the magnetization, and
+  magnetic energy.
+
+- benchmarck_damped_exchange:
+  simulates two spins interacting through the exchange
+  interaction. The parameters in the LAMMPS input script
+  (bench-spin-precession.in) are calibrated to match the
+  exchange definition in the python script (llg_exchange.py). 
+  Run as: ./run-bench-exchange.sh
+  Output: average magnetization resulting from the damped 
+  precession of the two interacting spins. Also plots the 
+  evolution of the magnetic energy.
+  
+- benchmarck_langevin_precession:
+  simulates a single spin in a magnetic field and in contact
+  with a thermal bath, and compares the statistical averages of 
+  the output to the analytical result of the Langevin function. 
+  Run as: ./run-bench-prec.sh
+  Output: statistical average of the z-component of the
+  magnetization (along the applied field) and of the magnetic
+  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).  

examples/SPIN/benchmark/benchmarck_damped_precession/bench-spin-precession.in

@@ -15,7 +15,6 @@ create_atoms    1 box
 mass            1 1.0
 set             type 1 spin 2.0 1.0 0.0 0.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
 
@@ -24,7 +23,7 @@ group bead      type 1
 variable        H equal 10.0
 variable        Kan equal 0.0
 variable        Temperature equal 0.0 
-variable        RUN equal 100000
+variable        Nsteps equal 500000
 
 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
@@ -45,4 +44,4 @@ thermo          100
 
 timestep	0.0001
 
-run             ${RUN}
+run             ${Nsteps}

examples/SPIN/benchmark/benchmarck_damped_precession/llg_precession.py

@@ -17,7 +17,7 @@ Bext = np.array([0.0, 0.0, 1.0])
 Sn = 2.0                # spin norm (in # of muB)
 S = np.array([1.0, 0.0, 0.0])
 
-N=100000                # number of timesteps
+N=500000                # number of timesteps
 dt=0.1                  # timestep (fs)
 
 # Rodrigues rotation formula
@@ -46,6 +46,7 @@ for t in range (0,N):
   theta=dt*np.linalg.norm(wf)
   axis=wf/np.linalg.norm(wf)
   S = np.dot(rotation_matrix(axis, theta), S)
+  en = -hbar*gyro*Sn*Bnrm*np.dot(S,Bext)
   # print res. in ps for comparison with LAMMPS
-  print(t*dt/1000.0,S[0],S[1],S[2])
+  print(t*dt/1000.0,S[0],S[1],S[2],en)
 

examples/SPIN/benchmark/benchmarck_damped_precession/plot_precession.py

@@ -15,25 +15,31 @@ if len(argv) != 3:
 lammps_file = sys.argv[1]
 llg_file = sys.argv[2]
 
-t_lmp,Sx_lmp,Sy_lmp,Sz_lmp = np.loadtxt(lammps_file, skiprows=0, usecols=(1,2,3,4),unpack=True)
-t_llg,Sx_llg,Sy_llg,Sz_llg = np.loadtxt(llg_file, skiprows=0, usecols=(0,1,2,3),unpack=True)
+t_lmp,Sx_lmp,Sy_lmp,Sz_lmp,en_lmp = np.loadtxt(lammps_file,
+        skiprows=0, usecols=(1,2,3,4,5),unpack=True)
+t_llg,Sx_llg,Sy_llg,Sz_llg,en_llg = np.loadtxt(llg_file, skiprows=0, usecols=(0,1,2,3,4),unpack=True)
 
 plt.figure()
-plt.subplot(311)
+plt.subplot(411)
 plt.ylabel('Sx')
 plt.plot(t_lmp, Sx_lmp, 'b-', label='LAMMPS')
 plt.plot(t_llg, Sx_llg, 'r--', label='LLG')
 
-plt.subplot(312)
+plt.subplot(412)
 plt.ylabel('Sy')
 plt.plot(t_lmp, Sy_lmp, 'b-', label='LAMMPS')
 plt.plot(t_llg, Sy_llg, 'r--', label='LLG')
 
-plt.subplot(313)
+plt.subplot(413)
 plt.ylabel('Sz')
 plt.plot(t_lmp, Sz_lmp, 'b-', label='LAMMPS')
 plt.plot(t_llg, Sz_llg, 'r--', label='LLG')
 
+plt.subplot(414)
+plt.ylabel('En (eV)')
+plt.plot(t_lmp, en_lmp, 'b-', label='LAMMPS')
+plt.plot(t_llg, en_llg, 'r--', label='LLG')
+
 plt.xlabel('time (in ps)')
 plt.legend()
 plt.show()

examples/SPIN/benchmark/benchmarck_langevin_exchange/bench-exchange-spin.template

@@ -1,41 +0,0 @@
-#LAMMPS in.run 
- 
-units           metal
-atom_style      spin
-atom_modify     map array
-boundary        p p p
-
-lattice         sc 3.0
-region          box block 0.0 2.0 0.0 2.0 0.0 2.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 3.1
-pair_coeff  	* * exchange 3.1 11.254 0.0 1.0
-
-variable        H equal 0.0
-variable        Kan equal 0.0
-variable        Temperature equal temperature
-variable        RUN equal 100000
-
-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             3 all langevin/spin ${Temperature} 0.01 12345
-
-thermo          500000
-thermo_style    custom step time temp vol 
-timestep        0.1
-
-compute         compute_spin all spin
-variable 	mag_energy equal c_compute_spin[5]
-variable 	AVEs equal c_compute_spin[4]
-
-fix 		avespin all ave/time 1 ${RUN} ${RUN} v_Temperature v_H v_Kan v_AVEs v_mag_energy file _av_spin
-
-run             ${RUN}
-
-shell           cat _av_spin

examples/SPIN/benchmark/benchmarck_langevin_exchange/langevin-exchange.py

@@ -1,34 +0,0 @@
-#!/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)
-J0=0.05                 # per-neighbor exchange interaction (eV)
-z=6                     # number of NN (bcc)
-g=2.0                   # Lande factor (adim)
-Hz=10.0                # mag. field (T)
-
-#Definition of the Langevin function
-def func(sm,t):
-    return mp.coth(z*J0*sm/(kb*t))-1.0/(z*J0*sm/(kb*t))
-
-npoints=200
-tolerance=1e-5
-ti=0.01
-tf=2000.0
-sz=1.0
-szg=0.5
-for i in range (0,npoints):
-    temp=ti+i*(tf-ti)/npoints
-    count=0
-    sz=1.0
-    szg=0.5
-    while (abs(sz-szg)) >= tolerance:
-        sz=szg
-        szg=func(sz,temp)
-        count+=1
-    emag=-z*J0*sz*sz
-    print('%d %lf %lf %lf %lf' % (temp,szg,sz,emag,count))

examples/SPIN/benchmark/benchmarck_langevin_exchange/plot_exchange.py

@@ -1,34 +0,0 @@
-#!/usr/bin/env python3
-
-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,2,3),unpack=True)
-
-plt.figure()
-plt.subplot(211)
-plt.ylabel('<S>')
-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)')
-plt.legend()
-plt.show()

examples/SPIN/benchmark/benchmarck_langevin_exchange/run-bench-exchange.sh

@@ -1,28 +0,0 @@
-#!/bin/bash
-
-# set initial and final temperature (K)
-tempi=0.0
-tempf=2000.0
-
-rm res_*.dat
-
-# run Lammps calculation
-N=20
-for (( i=0; i<$N; i++ ))
-do
-  temp="$(echo "$tempi+$i*($tempf-$tempi)/$N" | bc -l)"
-  sed s/temperature/${temp}/g bench-exchange-spin.template > \
-    bench-exchange-spin.in
-  ../../../../src/lmp_serial \
-    -in bench-exchange-spin.in 
-  Hz="$(tail -n 1 _av_spin | awk -F " " '{print $3}')"
-  sm="$(tail -n 1 _av_spin | awk -F " " '{print $5}')"
-  en="$(tail -n 1 _av_spin | awk -F " " '{print $6}')"
-  echo $temp $Hz $sm $en >> res_lammps.dat
-done
-
-# run Langevin function calculation
-python3 -m langevin-exchange.py > res_langevin.dat
-
-# plot comparison
-python3 -m plot_exchange.py res_lammps.dat res_langevin.dat

examples/SPIN/benchmark/benchmarck_langevin_precession/bench-prec-spin.template

@@ -33,14 +33,14 @@ fix             3 all langevin/spin ${Temperature} 0.01 12345
 
 compute         compute_spin all spin
 compute         outsp all property/atom spx spy spz sp
-compute         AVEsz all reduce ave c_outsp[3]
+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_AVEsz v_magnetic_energy file _av_spin
+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}

examples/SPIN/benchmark/benchmarck_langevin_precession/run-bench-prec.sh

@@ -14,9 +14,9 @@ do
     bench-prec-spin.in
   ./../../../../src/lmp_serial \
     -in bench-prec-spin.in 
-  Hz="$(tail -n 1 _av_spin | awk -F " " '{print $3}')"
-  sz="$(tail -n 1 _av_spin | awk -F " " '{print $5}')"
-  en="$(tail -n 1 _av_spin | awk -F " " '{print $6}')"
+  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
 

src/SPIN/compute_spin.cpp

@@ -104,7 +104,6 @@ void ComputeSpin::compute_vector()
         mag[0] += sp[i][0];
         mag[1] += sp[i][1];
         mag[2] += sp[i][2];
-        // magenergy -= 2.0*(sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]);
         magenergy -= (sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]);
         tx = sp[i][1]*fm[i][2]-sp[i][2]*fm[i][1];
         ty = sp[i][2]*fm[i][0]-sp[i][0]*fm[i][2];
@@ -129,7 +128,6 @@ 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 /= (kb*tempdenomtot);
   spintemperature /= (2.0*kb*tempdenomtot);
 
   vector[0] = magtot[0];

src/SPIN/fix_precession_spin.cpp

@@ -257,7 +257,6 @@ void FixPrecessionSpin::post_force(int /* vflag */)
 
       if (zeeman_flag) {          // compute Zeeman interaction
         compute_zeeman(i,fmi);
-        // epreci -= 2.0*hbar*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
         epreci -= hbar*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
       }
 
@@ -296,9 +295,6 @@ void FixPrecessionSpin::compute_single_precession(int i, double spi[3], double f
 void FixPrecessionSpin::compute_zeeman(int i, double fmi[3])
 {
   double **sp = atom->sp;
-  // fmi[0] += 0.5*sp[i][3]*hx;
-  // fmi[1] += 0.5*sp[i][3]*hy;
-  // fmi[2] += 0.5*sp[i][3]*hz;
   fmi[0] += sp[i][3]*hx;
   fmi[1] += sp[i][3]*hy;
   fmi[2] += sp[i][3]*hz;
@@ -309,9 +305,9 @@ void FixPrecessionSpin::compute_zeeman(int i, double fmi[3])
 void FixPrecessionSpin::compute_anisotropy(double spi[3], double fmi[3])
 {
   double scalar = nax*spi[0] + nay*spi[1] + naz*spi[2];
-  fmi[0] += 0.5*scalar*Kax;
-  fmi[1] += 0.5*scalar*Kay;
-  fmi[2] += 0.5*scalar*Kaz;
+  fmi[0] += scalar*Kax;
+  fmi[1] += scalar*Kay;
+  fmi[2] += scalar*Kaz;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -320,7 +316,7 @@ double FixPrecessionSpin::compute_anisotropy_energy(double spi[3])
 {
   double energy = 0.0;
   double scalar = nax*spi[0] + nay*spi[1] + naz*spi[2];
-  energy = 2.0*Ka*scalar*scalar;
+  energy = Ka*scalar*scalar;
   return energy; 
 }
 
@@ -358,9 +354,9 @@ void FixPrecessionSpin::compute_cubic(double spi[3], double fmi[3])
   sixy = k2ch*(nc1y*six1 + nc2y*six2 + nc3y*six3);
   sixz = k2ch*(nc1z*six1 + nc2z*six2 + nc3z*six3);
   
-  fmi[0] += 0.5*(fourx + sixx);
-  fmi[1] += 0.5*(foury + sixy);
-  fmi[2] += 0.5*(fourz + sixz);
+  fmi[0] += (fourx + sixx);
+  fmi[1] += (foury + sixy);
+  fmi[2] += (fourz + sixz);
 }
 
 /* ----------------------------------------------------------------------
@@ -379,7 +375,7 @@ double FixPrecessionSpin::compute_cubic_energy(double spi[3])
   energy = k1c*(skx*skx*sky*sky + sky*sky*skz*skz + skx*skx*skz*skz);
   energy += k2c*skx*skx*sky*sky*skz*skz;
 
-  return 2.0*energy;
+  return energy;
 }
 
 /* ---------------------------------------------------------------------- */

src/SPIN/pair_spin_exchange.cpp

@@ -242,7 +242,6 @@ 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;
         evdwl *= hbar;
       } else evdwl = 0.0;
 
@@ -352,11 +351,6 @@ void PairSpinExchange::compute_exchange(int i, int j, double rsq, double fmi[3],
   Jex *= (1.0-J2[itype][jtype]*ra);
   Jex *= exp(-ra);
 
-  // printf("Exchange : %g %g \n",Jex,Jex*hbar);
-
-  // fmi[0] += 2.0*Jex*spj[0];
-  // fmi[1] += 2.0*Jex*spj[1];
-  // fmi[2] += 2.0*Jex*spj[2];
   fmi[0] += Jex*spj[0];
   fmi[1] += Jex*spj[1];
   fmi[2] += Jex*spj[2];