Commit JT 010318 (before APS)

This commit is contained in:
julient31 2018-03-01 09:03:20 -07:00
parent 7990826ca2
commit e4c6c46b4c
10 changed files with 6143 additions and 24 deletions

View File

@ -23,8 +23,7 @@ set group all spin/random 11 2.50
pair_style hybrid/overlay pair/spin/exchange 6.0 pair/spin/me 4.5
pair_coeff * * pair/spin/exchange exchange 6.0 -0.01575 0.0 1.965
#pair_coeff * * pair/spin/me me 4.5 0.000109 1.0 1.0 1.0
pair_coeff * * pair/spin/me me 4.5 0.00109 1.0 1.0 1.0
pair_coeff * * pair/spin/me me 4.5 0.000109 1.0 1.0 1.0
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20

View File

@ -22,11 +22,9 @@ mass 1 58.93
set group all spin/random 31 1.72
velocity all create 100 4928459 rot yes dist gaussian
#pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0 pair/spin/soc/neel 4.0
pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0
pair_coeff * * eam/alloy ../examples/SPIN/cobalt/Co_PurjaPun_2012.eam.alloy Co
pair_coeff * * eam/alloy ../examples/SPIN/cobalt_fcc/Co_PurjaPun_2012.eam.alloy Co
pair_coeff * * pair/spin/exchange exchange 4.0 0.0446928 0.003496 1.4885
#pair_coeff * * pair/spin/soc/neel neel 4.0 0.003330282 0.864159 2.13731
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
@ -55,5 +53,5 @@ thermo 10
#dump 1 all custom 50 dump_cobalt.lammpstrj type x y z spx spy spz
run 10
run 500

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,32 @@
#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 *
print("Loop begin")
#Definition of the Bethe-Slater function
def func(x,a,b,c):
return 4*a*((x/c)**2)*(1-b*(x/c)**2)*np.exp(-(x/c)**2)
#Exchange coeff table (data to fit)
rdata, Jdata = np.loadtxt('exchange_hcp_co.dat', usecols=(0,1), unpack=True)
plt.plot(rdata, Jdata, 'b-', label='data')
#Perform the fit
popt, pcov = curve_fit(func, rdata, Jdata, bounds=(0, [500.,5.,5.]))
plt.plot(rdata, func(rdata, *popt), 'r--', label='fit')
#Print the fitted params
print("Parameters: a={:.10} (in meV), b={:.10} (adim), c={:.10} (in Ang)".format(*popt))
#Ploting the result
plt.xlabel('r_ij')
pylab.xlim([0,6.5])
plt.ylabel('J_ij')
plt.legend()
plt.show()
print("Loop end")

View File

@ -0,0 +1,9 @@
2.25569176882662 73.37931034482759
2.3817863397548162 47.99999999999999
2.4518388791593697 34.39080459770115
2.507880910683012 31.816091954022987
2.5359019264448337 28.137931034482747
2.5779334500875657 25.011494252873554
2.6339754816112086 19.126436781609186
2.760070052539404 13.241379310344826
3.5446584938704033 6.068965517241367

View File

@ -0,0 +1,60 @@
# fcc cobalt in a 3d periodic box
clear
units metal
atom_style spin
dimension 3
boundary p p p
# check why?
atom_modify map array
#lattice hcp 2.5071 2.5071 4.0695
lattice hcp 2.5071
region box block 0.0 5.0 0.0 5.0 0.0 5.0
create_box 1 box
create_atoms 1 box
# setting mass, mag. moments, and interactions for cobalt
mass 1 58.93
set group all spin/random 31 1.72
velocity all create 100 4928459 rot yes dist gaussian
#pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0 pair/spin/soc/neel 4.0
pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0
pair_coeff * * eam/alloy ../examples/SPIN/cobalt_hcp/Co_PurjaPun_2012.eam.alloy Co
pair_coeff * * pair/spin/exchange exchange 4.0 0.3593 1.135028015e-05 1.064568567
#pair_coeff * * pair/spin/soc/neel neel 4.0 0.0048 0.234 1.168 2.6905 0.705 0.652
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
fix 1 all force/spin zeeman 0.01 0.0 0.0 1.0
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all integration/spin lattice yes
timestep 0.0001
compute out_mag all compute/spin
compute out_pe all pe
compute out_ke all ke
compute out_temp all temp
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]
variable mag_force equal f_1
thermo_style custom step time v_magnorm v_emag temp etotal
thermo 10
dump 1 all custom 50 dump_cobalt.lammpstrj type x y z spx spy spz
run 1000

View File

@ -48,8 +48,8 @@ velocity all create 200 4928459 rot yes dist gaussian
#Magneto-mechanic interactions for bulk fcc Cobalt
#pair_style pair/spin/exchange 4.0
#pair_style eam/alloy
#pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0 pair/spin/soc/neel 4.0
pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0
pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0 pair/spin/soc/neel 4.0
#pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0
pair_coeff * * eam/alloy ../examples/SPIN/dev/Co_PurjaPun_2012.eam.alloy Co
#pair_coeff * * ../Co_PurjaPun_2012.eam.alloy Co
@ -69,7 +69,7 @@ pair_coeff * * pair/spin/exchange exchange 4.0 0.0446928 0.003496 1.4885
#type i and j | interaction type | cutoff | K1 (eV) | K2 (adim) | K3 (Ang) (for SOC)
#pair_coeff * * pair/spin/soc/neel neel 4.0 0.003330282 0.864159 2.13731
#pair_coeff * * pair/spin/soc/neel neel 4.0 0.0048 0.234 1.168 2.6905 0.705 0.652
pair_coeff * * pair/spin/soc/neel neel 4.0 0.0048 0.234 1.168 2.6905 0.705 0.652
#pair_coeff * * pair/spin/soc/neel neel 4.0 0.0 0.864159 2.13731
@ -81,13 +81,13 @@ neigh_modify every 10 check yes delay 20
#Magnetic field fix
#Type | Intensity (T or eV) | Direction
fix 1 all force/spin zeeman 10.0 0.0 0.0 1.0
fix 1 all force/spin zeeman 0.0 0.0 0.0 1.0
#fix 1 all force/spin anisotropy 0.01 0.0 0.0 1.0
#Fix Langevin spins (merging damping and temperature)
#Temp | Alpha_trans | Alpha_long | Seed
#fix 2 all langevin/spin 0.0 0.1 0.0 21
fix 2 all langevin/spin 0.0 0.1 21
#fix 2 all langevin/spin 0.0 0.1 21
fix 2 all langevin/spin 0.0 0.0 21
#Magnetic integration fix
fix 3 all integration/spin lattice yes
@ -116,8 +116,8 @@ variable emag equal c_out_mag[5]
variable tmag equal c_out_mag[6]
thermo 50
thermo_style custom step time v_magx v_magy v_magy v_magnorm v_emag v_tmag
#thermo_style custom step time pe ke v_emag
#thermo_style custom step time v_magx v_magy v_magy v_magnorm v_emag v_tmag
thermo_style custom step time pe ke v_emag etotal
thermo_modify format float %20.15g
#Dump the positions and spin directions of magnetic particles (vmd format)
@ -125,5 +125,5 @@ thermo_modify format float %20.15g
#Running the simulations for N timesteps
#run 10
run 100
run 20000

View File

@ -40,7 +40,6 @@
#include "pair_spin_me.h"
#include "pair_spin_soc_dmi.h"
#include "pair_spin_soc_neel.h"
#include "pair_spin_soc_landau.h"
#include "respa.h"
#include "update.h"

View File

@ -220,11 +220,12 @@ void PairSpinSocNeel::compute_soc_neel(int i, int j, double rsq, double eij[3],
q1ij *= exp(-ra);
q2ij = (-2.0*q1ij/9.0);
pq1x = -(scalar_eij_si*scalar_eij_si - scalar_si_sj/3.0)*spj[0]/3.0;
pq1y = -(scalar_eij_si*scalar_eij_si - scalar_si_sj/3.0)*spj[1]/3.0;
pq1z = -(scalar_eij_si*scalar_eij_si - scalar_si_sj/3.0)*spj[2]/3.0;
double scalar_eij_si_2 = scalar_eij_si*scalar_eij_si;
pq1x = -(scalar_eij_si_2*scalar_eij_si_2 - scalar_si_sj/3.0)*spj[0]/3.0;
pq1y = -(scalar_eij_si_2*scalar_eij_si_2 - scalar_si_sj/3.0)*spj[1]/3.0;
pq1z = -(scalar_eij_si_2*scalar_eij_si_2 - scalar_si_sj/3.0)*spj[2]/3.0;
double pqt1 = scalar_eij_sj*scalar_eij_sj-scalar_si_sj/3.0;
double pqt1 = (scalar_eij_sj*scalar_eij_sj-scalar_si_sj/3.0);
pq1x += pqt1*(2.0*scalar_eij_si*eij[0] - spj[0]/3.0);
pq1y += pqt1*(2.0*scalar_eij_si*eij[1] - spj[1]/3.0);
pq1z += pqt1*(2.0*scalar_eij_si*eij[2] - spj[2]/3.0);
@ -233,7 +234,6 @@ void PairSpinSocNeel::compute_soc_neel(int i, int j, double rsq, double eij[3],
pq1y *= q1ij;
pq1z *= q1ij;
double scalar_eij_si_2 = scalar_eij_si*scalar_eij_si;
double scalar_eij_sj_3 = scalar_eij_sj*scalar_eij_sj*scalar_eij_sj;
pq2x = 3.0*scalar_eij_si_2*scalar_eij_sj*eij[0] + scalar_eij_sj_3*eij[0];
pq2y = 3.0*scalar_eij_si_2*scalar_eij_sj*eij[1] + scalar_eij_sj_3*eij[1];
@ -286,8 +286,6 @@ void PairSpinSocNeel::compute_soc_mech_neel(int i, int j, double rsq, double eij
dgij = 1.0-ra-g2[itype][jtype]*ra*(2.0-ra);
dgij *= 8.0*g_mech*rr*exp(-ra);
// this 2.0*gij/drij is in Beaujouan calc. but not recovered yet
//double pdt1 = dgij*scalar_eij_si*scalar_eij_sj;
double pdt1 = (dgij-2.0*gij/drij)*scalar_eij_si*scalar_eij_sj;
pdt1 -= scalar_si_sj*dgij/3.0;
double pdt2 = scalar_eij_sj*gij/drij;
@ -323,13 +321,20 @@ void PairSpinSocNeel::compute_soc_mech_neel(int i, int j, double rsq, double eij
double scalar_eij_si_3 = scalar_eij_si*scalar_eij_si*scalar_eij_si;
double scalar_eij_sj_3 = scalar_eij_sj*scalar_eij_sj*scalar_eij_sj;
double scalar_si_sj_2 = scalar_si_sj*scalar_si_sj;
double pqt3 = 2.0*scalar_eij_si*scalar_eij_sj_3/drij;
/* double pqt3 = 2.0*scalar_eij_si*scalar_eij_sj_3/drij;
double pqt4 = 2.0*scalar_eij_sj*scalar_eij_si_3/drij;
double pqt5 = -2.0*scalar_si_sj*scalar_eij_si/(3.0*drij);
double pqt6 = -2.0*scalar_si_sj*scalar_eij_sj/(3.0*drij);
// pq1x += q1ij*(spi[0]*(pqt3+pqt6) + spj[0]*(pqt4+));
pq1x += q1ij*(pqt3*spi[0]+pqt4*spj[0]+pqt5*spi[0]+pqt6*spi[0]);
pq1y += q1ij*(pqt3*spi[1]+pqt4*spj[1]+pqt5*spi[1]+pqt6*spj[1]);
pq1z += q1ij*(pqt3*spi[2]+pqt4*spj[2]+pqt5*spi[2]+pqt6*spj[2]);
*/
double pqt3 = 2.0*scalar_eij_si*(scalar_eij_sj_2-scalar_si_sj/3.0)/drij;
double pqt4 = 2.0*scalar_eij_sj*(scalar_eij_si_2-scalar_si_sj/3.0)/drij;
pq1x += q1ij*(pqt3*spi[0] + pqt4*spj[0]);
pq1y += q1ij*(pqt3*spi[1] + pqt4*spj[1]);
pq1z += q1ij*(pqt3*spi[2] + pqt4*spj[2]);
double pqt7 = 4.0*scalar_eij_si_2*scalar_eij_sj_2/drij;
double pqt8 = 2.0*scalar_si_sj_2*scalar_eij_sj/(3.0*drij);
double pqt9 = 2.0*scalar_si_sj_2*scalar_eij_si/(3.0*drij);
@ -337,11 +342,22 @@ void PairSpinSocNeel::compute_soc_mech_neel(int i, int j, double rsq, double eij
pq1y -= q1ij*(pqt7 + pqt8 + pqt9)*eij[1];
pq1z -= q1ij*(pqt7 + pqt8 + pqt9)*eij[2];
/*
double pqt3 = 2.0*scalar_eij_si*(scalar_eij_sj_2-scalar_si_sj/3.0)/drij;
double pqt4 = 2.0*scalar_eij_sj*(scalar_eij_si_2-scalar_si_sj/3.0)/drij;
pq1x += q1ij*(pqt3*spi[0] + pqt4*spj[0]);
pq1y += q1ij*(pqt3*spi[1] + pqt4*spj[1]);
pq1z += q1ij*(pqt3*spi[2] + pqt4*spj[2]);
*/
//double scalar_eij_si_3 = scalar_eij_si*scalar_eij_si*scalar_eij_si;
//double scalar_eij_sj_3 = scalar_eij_sj*scalar_eij_sj*scalar_eij_sj;
double pqt10 = scalar_eij_sj*scalar_eij_si_3;
double pqt11 = scalar_eij_si*scalar_eij_sj_3;
pq2x = dq2ij*(pqt10 + pqt11)*eij[0];
pq2y = dq2ij*(pqt10 + pqt11)*eij[1];
pq2z = dq2ij*(pqt10 + pqt11)*eij[2];
double pqt12 = scalar_eij_si_3/drij;
double pqt13 = scalar_eij_sj_3/drij;
double pqt14 = 3.0*scalar_eij_sj*scalar_eij_si_2/drij;