Commit JT 111319

#LAMMPS in.run 

units           metal

atom_style      spin

atom_modify     map array

boundary        f f f 

read_data	two_spins.data

pair_style      spin/exchange 3.1

pair_coeff	* * exchange 3.1 11.254 0.0 1.0

group bead      type 1  

variable        H equal 0.0

variable        Kan equal 0.0

variable        Temperature equal 0.0 

variable        RUN equal 30000

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

compute		out_mag    all spin

compute		out_pe     all pe

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]

thermo_style    custom step time v_magx v_magy v_magz v_emag pe etotal

thermo          10

timestep	0.0001

run             ${RUN}

#!/usr/bin/env python3

import numpy as np , pylab, tkinter

import math

import matplotlib.pyplot as plt

import mpmath as mp

hbar=0.658212           # Planck's constant (eV.fs/rad)

J0=0.05                 # per-neighbor exchange interaction (eV)

S1 = np.array([1.0, 0.0, 0.0])

S2 = np.array([0.0, 1.0, 0.0])

alpha=0.01              # damping coefficient

pi=math.pi

N=30000                 # number of timesteps

dt=0.1                  # timestep (fs)

# Rodrigues rotation formula

def rotation_matrix(axis, theta):

  """

  Return the rotation matrix associated with counterclockwise

  rotation about the given axis by theta radians

  """

  axis = np.asarray(axis)

  a = math.cos(theta / 2.0)

  b, c, d = -axis * math.sin(theta / 2.0)

  aa, bb, cc, dd = a * a, b * b, c * c, d * d

  bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d

  return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],

      [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],

      [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])

# calculating precession field of spin Sr

def calc_rot_vector(Sr,Sf):

  rot = (J0/hbar)*(Sf-alpha*np.cross(Sf,Sr))/(1.0+alpha**2)

  return rot

# second-order ST decomposition as implemented in LAMMPS 

for t in range (0,N):

  # advance s1 by dt/4

  wf1 = calc_rot_vector(S1,S2)

  theta=dt*np.linalg.norm(wf1)*0.25

  axis=wf1/np.linalg.norm(wf1)

  S1 = np.dot(rotation_matrix(axis, theta), S1)

  # advance s2 by dt/2

  wf2 = calc_rot_vector(S2,S1)

  theta=dt*np.linalg.norm(wf2)*0.5

  axis=wf2/np.linalg.norm(wf2)

  S2 = np.dot(rotation_matrix(axis, theta), S2)

  # advance s1 by dt/2

  wf1 = calc_rot_vector(S1,S2)

  theta=dt*np.linalg.norm(wf1)*0.5

  axis=wf1/np.linalg.norm(wf1)

  S1 = np.dot(rotation_matrix(axis, theta), S1)

  # advance s2 by dt/2

  wf2 = calc_rot_vector(S2,S1)

  theta=dt*np.linalg.norm(wf2)*0.5

  axis=wf2/np.linalg.norm(wf2)

  S2 = np.dot(rotation_matrix(axis, theta), S2)

  # advance s1 by dt/4

  wf1 = calc_rot_vector(S1,S2)

  theta=dt*np.linalg.norm(wf1)*0.25

  axis=wf1/np.linalg.norm(wf1)

  S1 = np.dot(rotation_matrix(axis, theta), S1)

  # calc. average magnetization

  Sm = (S1+S2)*0.5

  # calc. energy

  # en = -hbar*(np.dot(S1,wf1)+np.dot(S2,wf2))

  en = -2.0*J0*(np.dot(S1,S2))

  # print res. in ps for comparison with LAMMPS

  print(t*dt/1000.0,Sm[0],Sm[1],Sm[2],en)

  # print(t*dt/1000.0,S1[0],S2[0],S1[1],S2[1],S1[2],S2[2],en)

#!/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_llg.dat")

  sys.exit()

lammps_file = sys.argv[1]

llg_file = sys.argv[2]

t_lmp,Sx_lmp,Sy_lmp,Sz_lmp,e_lmp = np.loadtxt(lammps_file,skiprows=0, usecols=(1,2,3,4,5),unpack=True)

t_llg,Sx_llg,Sy_llg,Sz_llg,e_llg = np.loadtxt(llg_file,skiprows=0, usecols=(0,1,2,3,4),unpack=True)

plt.figure()

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(412)

plt.ylabel('Sy')

plt.plot(t_lmp, Sy_lmp, 'b-', label='LAMMPS')

plt.plot(t_llg, Sy_llg, 'r--', label='LLG')

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('E (eV)')

plt.plot(t_lmp, e_lmp, 'b-', label='LAMMPS')

plt.plot(t_llg, e_llg, 'r--', label='LLG')

plt.xlabel('time (in ps)')

plt.legend()

plt.show()