jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Commit JT 042018, new spin/precession and pair/spin (peru virtual)

This commit is contained in:
julient31 2018-04-20 11:35:32 -06:00
parent ce80d1a3ea
commit 1b8669c620
27 changed files with 579 additions and 810 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

0
doc/src/Eqs/pair_spin_soc_dmi_interaction.jpg → doc/src/Eqs/pair_spin_dmi_interaction.jpg
View File

Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 7.0 KiB

After

Width:  |  Height:  |  Size: 7.0 KiB

0
doc/src/Eqs/pair_spin_soc_dmi_interaction.tex → doc/src/Eqs/pair_spin_dmi_interaction.tex
View File

22
doc/src/fix_precession_spin.txt
View File

@ -14,11 +14,11 @@ fix ID group precession/spin style args :pre
ID, group are documented in "fix"_fix.html command :ulb,l
precession/spin = style name of this fix command :l
style = {zeeman} or {aniso} :l
style = {zeeman} or {anisotropy} :l
{zeeman} args = H x y z
H = intensity of the magnetic field (in Tesla)
x y z = vector direction of the field
{aniso} args = K x y z
{anisotropy} args = K x y z
K = intensity of the magnetic anisotropy (in eV)
x y z = vector direction of the anisotropy :pre
:ule
@ -26,7 +26,8 @@ style = {zeeman} or {aniso} :l
[Examples:]
fix 1 all precession/spin zeeman 0.1 0.0 0.0 1.0
fix 1 all precession/spin aniso 0.001 0.0 0.0 1.0 :pre
fix 1 all precession/spin anisotropy 0.001 0.0 0.0 1.0
fix 1 all precession/spin zeeman 0.1 0.0 0.0 1.0 anisotropy 0.001 0.0 0.0 1.0 :pre
[Description:]
@ -41,7 +42,7 @@ magnetic field:
with mu0 the vacuum permeability, muB the Bohr magneton (muB = 5.788 eV/T
in metal units).
Style {aniso} is used to simulate an easy axis or an easy plane
Style {anisotropy} is used to simulate an easy axis or an easy plane
for the magnetic spins in the defined group:
:c,image(Eqs/force_spin_aniso.jpg)
@ -52,10 +53,23 @@ If K>0, an easy axis is defined, and if K<0, an easy plane is defined.
In both cases, the choice of (x y z) imposes the vector direction for the force.
Only the direction of the vector is important; it's length is ignored.
Both styles can be combined within one single command line.
:line
[Restart, fix_modify, output, run start/stop, minimize info:]
By default, the energy associated to this fix is not added to the potential
energy of the system.
The "fix_modify"_fix_modify.html {energy} option is supported by this fix
to add this magnetic potential energy to the potential energy of the system,
fix 1 all precession/spin zeeman 1.0 0.0 0.0 1.0
fix_modify 1 energy yes :pre
This fix computes a global scalar which can be accessed by various
"output commands"_Section_howto.html#howto_15.
No information about this fix is written to "binary restart
files"_restart.html.

8
doc/src/pair_spin_dmi.txt
View File

@ -6,11 +6,11 @@
:line
pair_style pair/spin/dmi command :h3
pair_style spin/dmi command :h3
[Syntax:]
pair_style pair/spin/dmi cutoff :pre
pair_style spin/dmi cutoff :pre
cutoff = global cutoff pair (distance in metal units) :ulb,l
@ -18,13 +18,13 @@ cutoff = global cutoff pair (distance in metal units) :ulb,l
[Examples:]
pair_style pair/spin/dmi 4.0
pair_style spin/dmi 4.0
pair_coeff * * dmi 2.6 0.001 1.0 0.0 0.0
pair_coeff 1 2 dmi 4.0 0.00109 0.0 0.0 1.0 :pre
[Description:]
Style {pair/spin/dmi} computes the Dzyaloshinskii-Moriya (DM) interaction
Style {spin/dmi} computes the Dzyaloshinskii-Moriya (DM) interaction
between pairs of magnetic spins:
:c,image(Eqs/pair_spin_dmi_interaction.jpg)

8
doc/src/pair_spin_exchange.txt
View File

@ -6,11 +6,11 @@
:line
pair_style pair/spin/exchange command :h3
pair_style spin/exchange command :h3
[Syntax:]
pair_style pair/spin/exchange cutoff :pre
pair_style spin/exchange cutoff :pre
cutoff = global cutoff pair (distance in metal units) :ulb,l
@ -18,13 +18,13 @@ cutoff = global cutoff pair (distance in metal units) :ulb,l
[Examples:]
pair_style pair/spin/exchange 4.0
pair_style spin/exchange 4.0
pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885
pair_coeff 1 2 exchange 6.0 -0.01575 0.0 1.965 :pre
[Description:]
Style {pair/spin/exchange} computes the exchange interaction between
Style {spin/exchange} computes the exchange interaction between
pairs of magnetic spins:
:c,image(Eqs/pair_spin_exchange_interaction.jpg)

10
doc/src/pair_spin_me.txt
View File

@ -6,11 +6,11 @@
:line
pair_style pair/spin/me command :h3
pair_style spin/me command :h3
[Syntax:]
pair_style pair/spin/me cutoff :pre
pair_style spin/me cutoff :pre
cutoff = global cutoff pair (distance in metal units) :ulb,l
@ -18,12 +18,12 @@ cutoff = global cutoff pair (distance in metal units) :ulb,l
[Examples:]
pair_style pair/spin/me 4.5
pair_coeff * * pair/spin/me me 4.5 0.00109 1.0 1.0 1.0 :pre
pair_style spin/me 4.5
pair_coeff * * me 4.5 0.00109 1.0 1.0 1.0 :pre
[Description:]
Style {pair/spin/me} computes a magneto-electric interaction between
Style {spin/me} computes a magneto-electric interaction between
pairs of magnetic spins. According to the derivation reported in
"(Katsura)"_#Katsura1, this interaction is defined as:

8
doc/src/pair_spin_neel.txt
View File

@ -6,11 +6,11 @@
:line
pair_style pair/spin/neel command :h3
pair_style spin/neel command :h3
[Syntax:]
pair_style pair/spin/neel cutoff :pre
pair_style spin/neel cutoff :pre
cutoff = global cutoff pair (distance in metal units) :ulb,l
@ -18,13 +18,13 @@ cutoff = global cutoff pair (distance in metal units) :ulb,l
[Examples:]
pair_style pair/spin/neel 4.0
pair_style spin/neel 4.0
pair_coeff * * neel 4.0 0.0048 0.234 1.168 2.6905 0.705 0.652
pair_coeff 1 2 neel 4.0 0.0048 0.234 1.168 0.0 0.0 1.0 :pre
[Description:]
Style {pair/spin/neel} computes the Neel pair anisotropy model
Style {spin/neel} computes the Neel pair anisotropy model
between pairs of magnetic spins:
:c,image(Eqs/pair_spin_dmi_interaction.jpg)

8
examples/SPIN/bfo/in.spin.bfo
View File

@ -21,16 +21,16 @@ mass 1 1.0
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_style hybrid/overlay spin/exchange 6.0 spin/me 4.5
pair_coeff * * spin/exchange exchange 6.0 -0.01575 0.0 1.965
pair_coeff * * 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
fix 1 all precession/spin anisotropy 0.0000033 0.0 0.0 1.0
fix 2 all langevin/spin 0.0 0.1 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice no
timestep 0.0002

21
examples/SPIN/cobalt_fcc/in.spin.cobalt_fcc
View File

@ -19,17 +19,20 @@ create_atoms 1 box
mass 1 58.93
set group all spin/random 31 1.72
#set group all spin/random 31 1.72
set group all spin 1.72 0.0 0.0 1.0
velocity all create 100 4928459 rot yes dist gaussian
pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0
pair_style hybrid/overlay eam/alloy spin/exchange 4.0
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 * * spin/exchange exchange 4.0 0.0446928 0.003496 1.4885
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
fix 1 all precession/spin zeeman 1.0 0.0 0.0 1.0
fix_modify 1 energy yes
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
@ -42,16 +45,20 @@ compute out_pe all pe
compute out_ke all ke
compute out_temp all temp
thermo_style custom f_1
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_style custom step time v_magnorm v_emag temp etotal
thermo 10
thermo_style custom step time f_1 v_magx v_magy v_magnorm v_emag temp etotal
thermo 50
compute outsp all property/atom spx spy spz sp fmx fmy fmz
dump 100 all custom 1 dump_cobalt_fcc.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3]
run 100
run 1000

8
examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp
View File

@ -22,11 +22,11 @@ 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/neel 4.0
pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0
#pair_style hybrid/overlay eam/alloy spin/exchange 4.0 spin/neel 4.0
pair_style hybrid/overlay eam/alloy 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/neel neel 4.0 0.0048 0.234 1.168 2.6905 0.705 0.652
pair_coeff * * spin/exchange exchange 4.0 0.3593 1.135028015e-05 1.064568567
#pair_coeff * * spin/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

4
examples/SPIN/iron/in.spin.iron
View File

@ -22,9 +22,9 @@ mass 1 55.845
set group all spin/random 31 2.2
velocity all create 100 4928459 rot yes dist gaussian
pair_style hybrid/overlay eam/alloy pair/spin/exchange 3.5
pair_style hybrid/overlay eam/alloy spin/exchange 3.5
pair_coeff * * eam/alloy ../examples/SPIN/iron/Fe_Mishin2006.eam.alloy Fe
pair_coeff * * pair/spin/exchange exchange 3.4 0.02726 0.2171 1.841
pair_coeff * * spin/exchange exchange 3.4 0.02726 0.2171 1.841
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20

4
examples/SPIN/nickel/in.spin.nickel
View File

@ -23,9 +23,9 @@ set group all spin/random 31 0.63
#set group all spin 0.63 0.0 0.0 1.0
velocity all create 100 4928459 rot yes dist gaussian
pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0
pair_style hybrid/overlay eam/alloy spin/exchange 4.0
pair_coeff * * eam/alloy ../examples/SPIN/nickel/Ni99.eam.alloy Ni
pair_coeff * * pair/spin/exchange exchange 4.0 0.50 0.2280246862 1.229983475
pair_coeff * * spin/exchange exchange 4.0 0.50 0.2280246862 1.229983475
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20

120
src/SPIN/fix_nve_spin.cpp
View File

@ -58,17 +58,18 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
rsec(NULL), stack_head(NULL), stack_foot(NULL),
backward_stacks(NULL), forward_stacks(NULL),
lockpairspin(NULL)
pair(NULL), spin_pairs(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal fix/NVE/spin command");
time_integrate = 1;
sector_flag = NONE;
lattice_flag = 1;
nlocal_max = 0;
npairs = 0;
npairspin = 0;
// checking if map array or hash is defined
@ -109,9 +110,8 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) :
// initialize the magnetic interaction flags
magpair_flag = 0;
magprecession_flag = 0;
zeeman_flag = aniso_flag = 0;
pair_spin_flag = 0;
precession_spin_flag = 0;
maglangevin_flag = 0;
tdamp_flag = temp_flag = 0;
@ -126,7 +126,7 @@ FixNVESpin::~FixNVESpin()
memory->destroy(stack_foot);
memory->destroy(forward_stacks);
memory->destroy(backward_stacks);
delete lockpairspin;
delete [] spin_pairs;
}
/* ---------------------------------------------------------------------- */
@ -151,28 +151,60 @@ void FixNVESpin::init()
dtf = 0.5 * update->dt * force->ftm2v;
dts = 0.25 * update->dt;
// ptrs on PairSpin classes
// set ptrs on Pair/Spin styles
lockpairspin = new PairSpin(lmp);
magpair_flag = lockpairspin->init_pair();
if (magpair_flag != 0 && magpair_flag != 1)
error->all(FLERR,"Incorrect value of magpair_flag");
// loop 1: obtain # of Pairs, and # of Pair/Spin styles
if (force->pair_match("spin",0,0)) { // only one Pair/Spin style
pair = force->pair_match("spin",0,0);
npairs = pair->instance_total;
npairspin = 1;
} else if (force->pair_match("spin",0,1)) { // more than one Pair/Spin style
pair = force->pair_match("spin",0,1);
npairs = pair->instance_total;
for (int i = 0; i<npairs; i++) {
if (force->pair_match("spin",0,i)) {
npairspin ++;
}
}
}
// init length of vector of ptrs to Pair/Spin styles
if (npairspin > 0) {
spin_pairs = new PairSpin*[npairspin];
}
// loop 2: fill vector with ptrs to Pair/Spin styles
int count = 0;
if (npairspin == 1) {
count = 1;
spin_pairs[0] = (PairSpin *) force->pair_match("spin",0,0);
} else if (npairspin > 1) {
for (int i = 0; i<npairs; i++) {
if (force->pair_match("spin",0,i)) {
spin_pairs[count] = (PairSpin *) force->pair_match("spin",0,i);
count++;
}
}
}
if (count != npairspin)
error->all(FLERR,"Incorrect number of spin pairs");
if (npairspin >= 1) pair_spin_flag = 1;
// ptrs FixPrecessionSpin classes
int iforce;
for (iforce = 0; iforce < modify->nfix; iforce++) {
if (strstr(modify->fix[iforce]->style,"precession/spin")) {
magprecession_flag = 1;
precession_spin_flag = 1;
lockprecessionspin = (FixPrecessionSpin *) modify->fix[iforce];
}
}
if (magprecession_flag) {
if (lockprecessionspin->zeeman_flag == 1) zeeman_flag = 1;
if (lockprecessionspin->aniso_flag == 1) aniso_flag = 1;
}
// ptrs on the FixLangevinSpin class
for (iforce = 0; iforce < modify->nfix; iforce++) {
@ -307,13 +339,13 @@ void FixNVESpin::initial_integrate(int vflag)
i = backward_stacks[i];
}
}
} else if (sector_flag == 0) { // serial seq. update
} else if (sector_flag == 0) { // serial seq. update
comm->forward_comm(); // comm. positions of ghost atoms
for (int i = 0; i < nlocal; i++){ // advance quarter s for nlocal-1
for (int i = 0; i < nlocal; i++){ // advance quarter s for nlocal-1
ComputeInteractionsSpin(i);
AdvanceSingleSpin(i);
}
for (int i = nlocal-2; i >= 0; i--){ // advance quarter s for nlocal-1
for (int i = nlocal-1; i >= 0; i--){ // advance quarter s for nlocal-1
ComputeInteractionsSpin(i);
AdvanceSingleSpin(i);
}
@ -381,9 +413,10 @@ void FixNVESpin::ComputeInteractionsSpin(int i)
double **sp = atom->sp;
double **fm = atom->fm;
int eflag = 1;
int vflag = 0;
int pair_compute_flag = 1;
// force computation for spin i
@ -394,20 +427,21 @@ void FixNVESpin::ComputeInteractionsSpin(int i)
fmi[0] = fmi[1] = fmi[2] = 0.0;
// evaluate magnetic pair interactions
// update magnetic pair interactions
if (magpair_flag) lockpairspin->compute_pair_single_spin(i,fmi);
// evaluate magnetic precession interactions
if (magprecession_flag) { // magnetic precession
if (zeeman_flag) { // zeeman
lockprecessionspin->compute_zeeman(i,fmi);
}
if (aniso_flag) { // aniso
lockprecessionspin->compute_anisotropy(i,spi,fmi);
if (pair_spin_flag) {
for (int k = 0; k < npairspin; k++) {
spin_pairs[k]->compute_single_pair(i,fmi);
}
}
// update magnetic precession interactions
if (precession_spin_flag) {
lockprecessionspin->compute_single_precession(i,spi,fmi);
}
// update langevin damping and random force
if (maglangevin_flag) { // mag. langevin
if (tdamp_flag) { // transverse damping
@ -418,7 +452,7 @@ void FixNVESpin::ComputeInteractionsSpin(int i)
}
}
// replace the magnetic force fm[i] by its new value
// replace the magnetic force fm[i] by its new value fmi
fm[i][0] = fmi[0];
fm[i][1] = fmi[1];
@ -445,10 +479,18 @@ void FixNVESpin::sectoring()
const double rsy = subhi[1] - sublo[1];
const double rsz = subhi[2] - sublo[2];
// temp
//const double rv = 2.0;
//const double rv = lockpairspinexchange->cut_spin_exchange_global;
const double rv = lockpairspin->larger_cutoff;
// extract larger cutoff from PairSpin styles
double rv, cutoff;
rv = cutoff = 0.0;
int dim = 0;
for (int i = 0; i < npairspin ; i++) {
cutoff = *((double *) spin_pairs[i]->extract("cut",dim));
rv = MAX(rv,cutoff);
}
if (rv == 0.0)
error->all(FLERR,"Illegal sectoring operation");
double rax = rsx/rv;
double ray = rsy/rv;
@ -509,7 +551,7 @@ void FixNVESpin::AdvanceSingleSpin(int i)
double **sp = atom->sp;
double **fm = atom->fm;
double dtfm,msq,scale,fm2,fmsq,sp2,spsq,energy,dts2;
double cp[3],g[3];
double cp[3],g[3];
cp[0] = cp[1] = cp[2] = 0.0;
g[0] = g[1] = g[2] = 0.0;

14
src/SPIN/fix_nve_spin.h
View File

@ -67,18 +67,22 @@ friend class PairSpin;
int nlocal_max; // max value of nlocal (for lists size)
int magpair_flag; // magnetic pair flags
int magprecession_flag; // magnetic precession flags
int zeeman_flag, aniso_flag;
int pair_spin_flag; // magnetic pair flags
int precession_spin_flag; // magnetic precession flags
int maglangevin_flag; // magnetic langevin flags
int tdamp_flag, temp_flag;
// pointers to magnetic interaction classes
// pointers to magnetic fixes
class PairSpin *lockpairspin;
class FixPrecessionSpin *lockprecessionspin;
class FixLangevinSpin *locklangevinspin;
// pointers to magnetic pair styles
int npairs, npairspin; // # of pairs, and # of spin pairs
class Pair *pair;
class PairSpin **spin_pairs; // vector of spin pairs
// sectoring variables
int nsectors;

64
src/SPIN/fix_precession_spin.cpp
View File

@ -43,14 +43,12 @@ using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
enum{ZEEMAN,ANISOTROPY};
enum{CONSTANT,EQUAL};
/* ---------------------------------------------------------------------- */
FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
if (narg < 7) error->all(FLERR,"Illegal precession/spin command");
// magnetic interactions coded for cartesian coordinates
@ -75,26 +73,28 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lm
Kax = Kay = Kaz = 0.0;
zeeman_flag = aniso_flag = 0;
if (strcmp(arg[3],"zeeman") == 0) {
if (narg != 8) error->all(FLERR,"Illegal precession/spin command");
style = ZEEMAN;
zeeman_flag = 1;
H_field = force->numeric(FLERR,arg[4]);
nhx = force->numeric(FLERR,arg[5]);
nhy = force->numeric(FLERR,arg[6]);
nhz = force->numeric(FLERR,arg[7]);
magfieldstyle = CONSTANT;
} else if (strcmp(arg[3],"anisotropy") == 0) {
if (narg != 8) error->all(FLERR,"Illegal precession/spin command");
style = ANISOTROPY;
aniso_flag = 1;
Ka = force->numeric(FLERR,arg[4]);
nax = force->numeric(FLERR,arg[5]);
nay = force->numeric(FLERR,arg[6]);
naz = force->numeric(FLERR,arg[7]);
} else error->all(FLERR,"Illegal precession/spin command");
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg],"zeeman") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix precession/spin command");
zeeman_flag = 1;
H_field = force->numeric(FLERR,arg[iarg+1]);
nhx = force->numeric(FLERR,arg[iarg+2]);
nhy = force->numeric(FLERR,arg[iarg+3]);
nhz = force->numeric(FLERR,arg[iarg+4]);
iarg += 5;
} else if (strcmp(arg[iarg],"anisotropy") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix precession/spin command");
aniso_flag = 1;
Ka = force->numeric(FLERR,arg[iarg+1]);
nax = force->numeric(FLERR,arg[iarg+2]);
nay = force->numeric(FLERR,arg[iarg+3]);
naz = force->numeric(FLERR,arg[iarg+4]);
iarg += 5;
} else error->all(FLERR,"Illegal precession/spin command");
}
degree2rad = MY_PI/180.0;
time_origin = update->ntimestep;
@ -171,7 +171,7 @@ void FixPrecessionSpin::setup(int vflag)
void FixPrecessionSpin::post_force(int vflag)
{
// update gravity due to variables
// update mag field with time (potential improvement)
if (varflag != CONSTANT) {
modify->clearstep_compute();
@ -209,7 +209,19 @@ void FixPrecessionSpin::post_force(int vflag)
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixPrecessionSpin::compute_single_precession(int i, double spi[3], double fmi[3])
{
if (zeeman_flag) {
compute_zeeman(i,fmi);
}
if (aniso_flag) {
compute_anisotropy(i,spi,fmi);
}
}
@ -244,12 +256,12 @@ void FixPrecessionSpin::post_force_respa(int vflag, int ilevel, int iloop)
void FixPrecessionSpin::set_magneticprecession()
{
if (style == ZEEMAN) {
if (zeeman_flag) {
hx = H_field*nhx;
hy = H_field*nhy;
hz = H_field*nhz;
}
if (style == ANISOTROPY) {
if (aniso_flag) {
Kax = 2.0*Ka*nax;
Kay = 2.0*Ka*nay;
Kaz = 2.0*Ka*naz;

7
src/SPIN/fix_precession_spin.h
View File

@ -47,9 +47,10 @@ class FixPrecessionSpin : public Fix {
virtual void post_force_respa(int, int, int);
double compute_scalar();
int zeeman_flag, aniso_flag;
void compute_zeeman(int, double [3]);
void compute_anisotropy(int, double [3], double [3]);
int zeeman_flag, aniso_flag;
void compute_single_precession(int, double *, double *);
void compute_zeeman(int, double *);
void compute_anisotropy(int, double *, double *);
protected:
int style; // style of the magnetic precession

526
src/SPIN/pair_spin.cpp
View File

@ -52,47 +52,16 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairSpin::PairSpin(LAMMPS *lmp) : Pair(lmp),
lockfixnvespin(NULL), pair_keyword(NULL), pair_spin_keywords(NULL),
exchange_spin_styles(NULL), dmi_spin_styles(NULL),
neel_spin_styles(NULL), me_spin_styles(NULL)
PairSpin::PairSpin(LAMMPS *lmp) : Pair(lmp)
{
hbar = force->hplanck/MY_2PI;
single_enable = 0;
no_virial_fdotr_compute = 1;
lattice_flag = 0;
// init # of Pair/Spin styles
nspinstyle = 0; // # of PairSpin styles
nexchangespinstyle = 0; // # of PairSpinExchange styles
ndmispinstyle = 0; // # of PairSpinDmi styles
nneelspinstyle = 0; // # of PairSpinNeel styles
nmespinstyle = 0; // # of PairSpinMe styles
// init larger Pair/Spin style cutoff
larger_cutoff = 0.0;
}
/* ---------------------------------------------------------------------- */
PairSpin::~PairSpin()
{
if (nspinstyle) {
for (int m = 0; m < nspinstyle; m++) {
delete [] pair_spin_keywords[m];
}
}
delete [] pair_keyword;
delete [] exchange_spin_styles;
delete [] dmi_spin_styles;
delete [] neel_spin_styles;
delete [] me_spin_styles;
delete [] pair_spin_keywords;
}
PairSpin::~PairSpin() {}
/* ----------------------------------------------------------------------
global settings
@ -110,488 +79,6 @@ void PairSpin::settings(int narg, char **arg)
}
/* ----------------------------------------------------------------------
global compute, defined in Pair/Spin subclasses
------------------------------------------------------------------------- */
void PairSpin::compute(int eflag, int vflag) {}
/* ----------------------------------------------------------------------
compute all Pair/Spin interactions for atom ii
------------------------------------------------------------------------- */
void PairSpin::compute_pair_single_spin(int ii, double fmi[3])
{
const int nlocal = atom->nlocal;
int *type = atom->type;
double **x = atom->x;
double **sp = atom->sp;
double xi[3], rij[3], eij[3];
double spi[3], spj[3];
int iexchange, idmi, ineel, ime;
int i,j,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq, rd, inorm;
iexchange = idmi = ineel = ime = 0;
for (int ipair=0; ipair < nspinstyle; ipair++) {
if (strstr(pair_spin_keywords[ipair],"pair/spin/exchange")) {
inum = exchange_spin_styles[iexchange]->list->inum;
ilist = exchange_spin_styles[iexchange]->list->ilist;
numneigh = exchange_spin_styles[iexchange]->list->numneigh;
firstneigh = exchange_spin_styles[iexchange]->list->firstneigh;
i = ilist[ii];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
itype = type[ii];
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
inorm = 1.0/sqrt(rsq);
exchange_spin_styles[iexchange]->compute_exchange(i,j,rsq,fmi,spi,spj);
}
iexchange++;
}
if (strstr(pair_spin_keywords[ipair],"pair/spin/dmi")) {
inum = dmi_spin_styles[idmi]->list->inum;
ilist = dmi_spin_styles[idmi]->list->ilist;
numneigh = dmi_spin_styles[idmi]->list->numneigh;
firstneigh = dmi_spin_styles[idmi]->list->firstneigh;
i = ilist[ii];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
itype = type[ii];
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
inorm = 1.0/sqrt(rsq);
dmi_spin_styles[idmi]->compute_dmi(i,j,rsq,fmi,spi,spj);
}
idmi++;
}
if (strstr(pair_spin_keywords[ipair],"pair/spin/neel")) {
inum = neel_spin_styles[ineel]->list->inum;
ilist = neel_spin_styles[ineel]->list->ilist;
numneigh = neel_spin_styles[ineel]->list->numneigh;
firstneigh = neel_spin_styles[ineel]->list->firstneigh;
i = ilist[ii];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
itype = type[ii];
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
inorm = 1.0/sqrt(rsq);
eij[0] = rij[0]*inorm;
eij[1] = rij[1]*inorm;
eij[2] = rij[2]*inorm;
neel_spin_styles[ineel]->compute_neel(i,j,rsq,eij,fmi,spi,spj);
}
ineel++;
}
if (strstr(pair_spin_keywords[ipair],"pair/spin/me")) {
inum = me_spin_styles[ime]->list->inum;
ilist = me_spin_styles[ime]->list->ilist;
numneigh = me_spin_styles[ime]->list->numneigh;
firstneigh = me_spin_styles[ime]->list->firstneigh;
i = ilist[ii];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
itype = type[ii];
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
inorm = 1.0/sqrt(rsq);
eij[0] = rij[0]*inorm;
eij[1] = rij[1]*inorm;
eij[2] = rij[2]*inorm;
me_spin_styles[ime]->compute_me(i,j,rsq,eij,fmi,spi,spj);
}
ime++;
}
}
}
/* ----------------------------------------------------------------------
called from Fix/NVE/Spin
initialize the # and the lists of Pair/Spin styles
------------------------------------------------------------------------- */
int PairSpin::init_pair()
{
int nsub = 0;
// getting the pair style
pair_keyword = new char[strlen(force->pair_style) + 1];
strcpy(pair_keyword,force->pair_style);
// searching for number of PairSpin styles
int temp_npair = 0;
nspinstyle = 0;
pair = pair_spin_match("spin",0,nsub);
// init lists of PairSpin styles
exchange_spin_styles = new PairSpinExchange*[nspinstyle];
dmi_spin_styles = new PairSpinDmi*[nspinstyle];
neel_spin_styles = new PairSpinNeel*[nspinstyle];
me_spin_styles = new PairSpinMe*[nspinstyle];
// init lists of PairSpin names
pair_spin_keywords = new char*[nspinstyle];
nexchangespinstyle = 0;
ndmispinstyle = 0;
nneelspinstyle = 0;
nmespinstyle = 0;
// loop to define lists of Pair/Spin styles
int ispin = 0;
if (strstr(pair_keyword,"spin")) {
// Pair/Spin/Exchange style
if (strstr(pair_keyword,"pair/spin/exchange")) {
int n = strlen(pair_keyword) + 1;
pair_spin_keywords[ispin] = new char[n];
strcpy(pair_spin_keywords[ispin],pair_keyword);
exchange_spin_styles[nexchangespinstyle] = (PairSpinExchange *) force->pair;
nexchangespinstyle++;
ispin++;
}
// Pair/Spin/Dmi style
if (strstr(pair_keyword,"pair/spin/dmi")) {
int n = strlen(pair_keyword) + 1;
pair_spin_keywords[ispin] = new char[n];
strcpy(pair_spin_keywords[ispin],pair_keyword);
dmi_spin_styles[ndmispinstyle] = (PairSpinDmi *) force->pair;
ndmispinstyle++;
ispin++;
}
// Pair/Spin/Neel style
if (strstr(pair_keyword,"pair/spin/neel")) {
int n = strlen(pair_keyword) + 1;
pair_spin_keywords[ispin] = new char[n];
strcpy(pair_spin_keywords[ispin],pair_keyword);
neel_spin_styles[nneelspinstyle] = (PairSpinNeel *) force->pair;
nneelspinstyle++;
ispin++;
}
// list Pair/Spin/Me styles
if (strstr(pair_keyword,"pair/spin/me")) {
int n = strlen(pair_keyword) + 1;
pair_spin_keywords[ispin] = new char[n];
strcpy(pair_spin_keywords[ispin],pair_keyword);
me_spin_styles[nmespinstyle] = (PairSpinMe *) force->pair;
nmespinstyle++;
ispin++;
}
} else if (strstr(pair_keyword,"hybrid/overlay")) { // if hybrid/overlay
PairHybrid *lockhybrid = (PairHybrid *) force->pair;
for (int i =0; i < lockhybrid->nstyles; i++) {
// error checks
if (strcmp(lockhybrid->keywords[i],"hybrid") == 0)
error->all(FLERR,"Pair style hybrid cannot have hybrid as an argument");
if (strcmp(lockhybrid->keywords[i],"none") == 0)
error->all(FLERR,"Pair style hybrid cannot have none as an argument");
// list Pair/Spin/Exchange styles
if (strstr(lockhybrid->keywords[i],"pair/spin/exchange")) {
int n = strlen(lockhybrid->keywords[i]) + 1;
pair_spin_keywords[ispin] = new char[n];
strcpy(pair_spin_keywords[ispin],lockhybrid->keywords[i]);
exchange_spin_styles[nexchangespinstyle] = (PairSpinExchange *) lockhybrid->styles[i];
nexchangespinstyle++;
ispin++;
}
// list Pair/Spin/Dmi styles
if (strstr(lockhybrid->keywords[i],"pair/spin/dmi")) {
int n = strlen(lockhybrid->keywords[i]) + 1;
pair_spin_keywords[ispin] = new char[n];
strcpy(pair_spin_keywords[ispin],lockhybrid->keywords[i]);
dmi_spin_styles[ndmispinstyle] = (PairSpinDmi *) lockhybrid->styles[i];
ndmispinstyle++;
ispin++;
}
// list Pair/Spin/Neel styles
if (strstr(lockhybrid->keywords[i],"pair/spin/neel")) {
int n = strlen(lockhybrid->keywords[i]) + 1;
pair_spin_keywords[ispin] = new char[n];
strcpy(pair_spin_keywords[ispin],lockhybrid->keywords[i]);
neel_spin_styles[nneelspinstyle] = (PairSpinNeel *) lockhybrid->styles[i];
nneelspinstyle++;
ispin++;
}
// list Pair/Spin/Me styles
if (strstr(lockhybrid->keywords[i],"pair/spin/me")) {
int n = strlen(lockhybrid->keywords[i]) + 1;
pair_spin_keywords[ispin] = new char[n];
strcpy(pair_spin_keywords[ispin],lockhybrid->keywords[i]);
me_spin_styles[nmespinstyle] = (PairSpinMe *) lockhybrid->styles[i];
nmespinstyle++;
ispin++;
}
}
} else if (strstr(pair_keyword,"hybrid")) { // no hybrid style with PairSpin
error->all(FLERR,"Pair/Spin styles need hybrid/overlay Pair style");
} else error->all(FLERR,"Wrong arguments in PairSpin style");
if (ispin != nspinstyle)
error->all(FLERR,"Wrong number of PairSpin styles");
if ((nexchangespinstyle + ndmispinstyle + nneelspinstyle + nmespinstyle) != nspinstyle)
error->all(FLERR,"Wrong number of PairSpin styles");
if (strstr(pair_keyword,"spin") && nspinstyle != 1)
error->all(FLERR,"Wrong number of PairSpin styles");
int mag_pair_flag = 0;
if (nspinstyle >= 1) mag_pair_flag = 1;
// get larger_cutoff
larger_cutoff = larger_spin_cutoff();
if (mag_pair_flag == 1 && larger_cutoff == 0.0)
error->all(FLERR,"Wrong arguments for PairSpin styles");
return mag_pair_flag;
}
/* ----------------------------------------------------------------------
get the larger Pair/Spin style cutoff for the sectoring operation
------------------------------------------------------------------------- */
double PairSpin::larger_spin_cutoff()
{
int iexchange, idmi, ineel, ime;
double local_cutoff = 0.0;
iexchange = idmi = ineel = ime = 0;
for (int ipair=0; ipair < nspinstyle; ipair++) {
if (strstr(pair_spin_keywords[ipair],"pair/spin/exchange")) {
if (local_cutoff < exchange_spin_styles[iexchange]->cut_spin_exchange_global)
local_cutoff = exchange_spin_styles[iexchange]->cut_spin_exchange_global;
iexchange++;
}
if (strstr(pair_spin_keywords[ipair],"pair/spin/dmi")) {
if (local_cutoff < dmi_spin_styles[idmi]->cut_spin_dmi_global)
local_cutoff = dmi_spin_styles[idmi]->cut_spin_dmi_global;
idmi++;
}
if (strstr(pair_spin_keywords[ipair],"pair/spin/neel")) {
if (local_cutoff < neel_spin_styles[ineel]->cut_spin_neel_global)
local_cutoff = neel_spin_styles[ineel]->cut_spin_neel_global;
ineel++;
}
if (strstr(pair_spin_keywords[ipair],"pair/spin/me")) {
if (local_cutoff < me_spin_styles[ime]->cut_spin_me_global)
local_cutoff = me_spin_styles[ime]->cut_spin_me_global;
ime++;
}
}
if ((iexchange + idmi + ineel + ime) != nspinstyle)
error->all(FLERR,"Wrong number of PairSpin styles");
return local_cutoff;
}
/* ---------------------------------------------------------------------- */
Pair *PairSpin::pair_spin_match(const char *word, int exact, int nsub)
{
int iwhich,count;
if (exact && strcmp(pair_keyword,word) == 0) return pair;
else if (!exact && strstr(pair_keyword,word)) return pair;
else if (strstr(pair_keyword,"hybrid/overlay")) {
PairHybrid *lockhybrid = (PairHybrid *) force->pair;
count = 0;
for (int i = 0; i < lockhybrid->nstyles; i++)
if ((exact && strcmp(lockhybrid->keywords[i],word) == 0) ||
(!exact && strstr(lockhybrid->keywords[i],word))) {
iwhich = i;
count++;
nspinstyle = count;
if (nsub == count) return lockhybrid->styles[iwhich];
}
if (count == 1) return lockhybrid->styles[iwhich];
} else if (strstr(pair_keyword,"hybrid")) {
PairHybrid *lockhybrid = (PairHybrid *) force->pair;
count = 0;
for (int i = 0; i < lockhybrid->nstyles; i++)
if ((exact && strcmp(lockhybrid->keywords[i],word) == 0) ||
(!exact && strstr(lockhybrid->keywords[i],word))) {
iwhich = i;
count++;
nspinstyle = count;
if (nsub == count) return lockhybrid->styles[iwhich];
}
if (count == 1) return lockhybrid->styles[iwhich];
}
return NULL;
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSpin::allocate()
{
allocated = 1;
char *pair_keyword = new char[strlen(force->pair_style) + 1];
}
/* ----------------------------------------------------------------------
set coeffs for one or more type spin pairs (only one for now)
------------------------------------------------------------------------- */
void PairSpin::coeff(int narg, char **arg)
{
const double hbar = force->hplanck/MY_2PI;
if (!allocated) allocate();
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
@ -610,13 +97,4 @@ void PairSpin::init_style()
}
if (ifix == modify->nfix)
error->all(FLERR,"pair/spin style requires nve/spin");
// get the lattice_flag from nve/spin
for (int i = 0; i < modify->nfix; i++) {
if (strcmp(modify->fix[i]->style,"nve/spin") == 0) {
lockfixnvespin = (FixNVESpin *) modify->fix[i];
lattice_flag = lockfixnvespin->lattice_flag;
}
}
}

54
src/SPIN/pair_spin.h
View File

@ -11,22 +11,8 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------
Contributing authors: Julien Tranchida (SNL)
Aidan Thompson (SNL)
Please cite the related publication:
Tranchida, J., Plimpton, S. J., Thibaudeau, P., & Thompson, A. P. (2018).
Massively parallel symplectic algorithm for coupled magnetic spin dynamics
and molecular dynamics. arXiv preprint arXiv:1801.10233.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(pair/spin,PairSpin)
#else
#ifndef LMP_PAIR_SPIN_H
#define LMP_PAIR_SPIN_H
@ -40,46 +26,18 @@ friend class FixNVESpin;
PairSpin(class LAMMPS *);
virtual ~PairSpin();
virtual void settings(int, char **);
virtual void coeff(int, char **);
virtual void coeff(int, char **) {}
virtual void init_style();
virtual double init_one(int, int) {return 0.0;}
virtual void compute(int, int);
virtual void *extract(const char *, int &) {return NULL;}
// functions called from fix/nve/spin
// used to evaluate force on spin i inside integration loops
int init_pair(); // init call of PairSpin styles
double larger_cutoff; // larger Pair/Spin style cutoff
double larger_spin_cutoff(); // compute larger_cutoff
void compute_pair_single_spin(int, double *); // compute pairs for one atom
class Pair *pair; // unused try class
int nspinstyle; // # of magnetic pairs
class Pair *pair_spin_match(const char *, int, int); // initializing nspinstyle
char *pair_keyword; // main pair style
char **pair_spin_keywords; // list of Pair/Spin style names
// # and lists of Pair/Spin style classes
int nexchangespinstyle; // # of exchange pairs
class PairSpinExchange **exchange_spin_styles; // list of Pair/Spin/Exchange2 classes
int ndmispinstyle; // # of dmi pairs
class PairSpinDmi **dmi_spin_styles; // list of Pair/Spin/Dmi2 classes
int nneelspinstyle; // # of dmi pairs
class PairSpinNeel **neel_spin_styles; // list of Pair/Spin/Dmi2 classes
int nmespinstyle; // # of me pairs
class PairSpinMe **me_spin_styles; // list of Pair/Spin/Me2 classes
virtual void compute(int, int) {}
virtual void compute_single_pair(int, double *) {}
protected:
int lattice_flag; // if 0 spins only, 1 if spin-lattice
double hbar; // Planck constant (eV.ps.rad-1)
class FixNVESpin *lockfixnvespin; // ptr to FixNVESpin2 for setups
double hbar; // Planck constant (eV.ps.rad-1)
virtual void allocate();
virtual void allocate() {}
};
}

148
src/SPIN/pair_spin_dmi.cpp
View File

@ -101,7 +101,7 @@ void PairSpinDmi::settings(int narg, char **arg)
void PairSpinDmi::coeff(int narg, char **arg)
{
const double hbar = force->hplanck/MY_2PI;
// const double hbar = force->hplanck/MY_2PI;
if (!allocated) allocate();
@ -165,16 +165,6 @@ void PairSpinDmi::init_style()
}
if (ifix == modify->nfix)
error->all(FLERR,"pair/spin style requires nve/spin");
// get the lattice_flag from nve/spin
for (int i = 0; i < modify->nfix; i++) {
if (strcmp(modify->fix[i]->style,"nve/spin") == 0) {
lockfixnvespin = (FixNVESpin *) modify->fix[i];
lattice_flag = lockfixnvespin->lattice_flag;
}
}
}
/* ----------------------------------------------------------------------
@ -189,23 +179,32 @@ double PairSpinDmi::init_one(int i, int j)
return cut_spin_dmi_global;
}
/* ----------------------------------------------------------------------
extract the larger cutoff
------------------------------------------------------------------------- */
void *PairSpinDmi::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"cut") == 0) return (void *) &cut_spin_dmi_global;
return NULL;
}
/* ---------------------------------------------------------------------- */
void PairSpinDmi::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
int i,j,ii,jj,inum,jnum;
double evdwl, ecoul;
double xi[3], rij[3];
double spi[3], spj[3];
double fi[3], fmi[3];
double cut_spin_dmi_2;
double rsq, rd, inorm;
double rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
cut_spin_dmi_2 = cut_spin_dmi_global*cut_spin_dmi_global;
double **x = atom->x;
double **f = atom->f;
@ -221,21 +220,26 @@ void PairSpinDmi::compute(int eflag, int vflag)
firstneigh = list->firstneigh;
// dmi computation
// loop over neighbors of my atoms
// loop over all atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
// loop on neighbors
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
@ -253,23 +257,11 @@ void PairSpinDmi::compute(int eflag, int vflag)
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
inorm = 1.0/sqrt(rsq);
rij[0] *= inorm;
rij[1] *= inorm;
rij[2] *= inorm;
itype = type[i];
jtype = type[j];
// compute magnetic and mechanical components of soc_dmi
cut_spin_dmi_2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];
if (rsq <= cut_spin_dmi_2){
compute_dmi(i,j,rsq,fmi,spi,spj);
if (lattice_flag) {
compute_dmi_mech(i,j,fi,spi,spj);
}
}
compute_dmi(i,j,rsq,fmi,spi,spj);
compute_dmi_mech(i,j,fi,spi,spj);
f[i][0] += fi[0];
f[i][1] += fi[1];
@ -284,15 +276,11 @@ void PairSpinDmi::compute(int eflag, int vflag)
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
if (eflag) {
if (rsq <= cut_spin_dmi_2) {
evdwl -= spi[0]*fmi[0];
evdwl -= spi[1]*fmi[1];
evdwl -= spi[2]*fmi[2];
evdwl *= hbar;
} else evdwl = 0.0;
}
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]);
@ -300,14 +288,75 @@ void PairSpinDmi::compute(int eflag, int vflag)
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairSpinDmi::compute_single_pair(int ii, double fmi[3])
{
const int nlocal = atom->nlocal;
int *type = atom->type;
double **x = atom->x;
double **sp = atom->sp;
double xi[3], rij[3];
double spi[3], spj[3];
int iexchange, idmi, ineel, ime;
int i,j,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq, rd, inorm;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
i = ilist[ii];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
itype = type[ii];
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
inorm = 1.0/sqrt(rsq);
compute_dmi(i,j,rsq,fmi,spi,spj);
}
}
/* ----------------------------------------------------------------------
compute the dmi interaction between spin i and spin j
------------------------------------------------------------------------- */
void PairSpinDmi::compute_dmi(int i, int j, double rsq, double fmi[3], double spi[3], double spj[3])
void PairSpinDmi::compute_dmi(int i, int j, double rsq, double fmi[3], double spi[3], double spj[3])
{
int *type = atom->type;
int itype, jtype;
@ -317,15 +366,16 @@ void PairSpinDmi::compute_dmi(int i, int j, double rsq, double fmi[3], double s
double local_cut2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];
if (rsq <= local_cut2) {
double dmix,dmiy,dmiz;
double dmix, dmiy, dmiz;
dmix = DM[itype][jtype]*v_dmx[itype][jtype];
dmiy = DM[itype][jtype]*v_dmy[itype][jtype];
dmiz = DM[itype][jtype]*v_dmz[itype][jtype];
dmix = DM[itype][jtype] * v_dmx[itype][jtype];
dmiy = DM[itype][jtype] * v_dmy[itype][jtype];
dmiz = DM[itype][jtype] * v_dmz[itype][jtype];
fmi[0] -= (spj[1]*dmiz - spj[2]*dmiy);
fmi[1] -= (spj[2]*dmix - spj[0]*dmiz);
fmi[2] -= (spj[0]*dmiy - spj[1]*dmix);
fmi[0] += spj[1]*dmiz-spj[2]*dmiy;
fmi[1] += spj[2]*dmix-spj[0]*dmiz;
fmi[2] += spj[0]*dmiy-spj[1]*dmix;
}
}

5
src/SPIN/pair_spin_dmi.h
View File

@ -23,7 +23,7 @@
#ifdef PAIR_CLASS
PairStyle(pair/spin/dmi,PairSpinDmi)
PairStyle(spin/dmi,PairSpinDmi)
#else
@ -42,8 +42,11 @@ class PairSpinDmi : public PairSpin {
void coeff(int, char **);
void init_style();
double init_one(int, int);
void *extract(const char *, int &);
void compute(int, int);
void compute_single_pair(int, double *);
void compute_dmi(int, int, double, double *, double *, double *);
void compute_dmi_mech(int, int, double *, double *, double *);

118
src/SPIN/pair_spin_exchange.cpp
View File

@ -45,10 +45,13 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairSpinExchange::PairSpinExchange(LAMMPS *lmp) : PairSpin(lmp)
PairSpinExchange::PairSpinExchange(LAMMPS *lmp) : PairSpin(lmp),
lockfixnvespin(NULL)
{
hbar = force->hplanck/MY_2PI;
single_enable = 0;
no_virial_fdotr_compute = 1;
lattice_flag = 0;
}
/* ---------------------------------------------------------------------- */
@ -99,8 +102,6 @@ void PairSpinExchange::settings(int narg, char **arg)
void PairSpinExchange::coeff(int narg, char **arg)
{
const double hbar = force->hplanck/MY_2PI;
if (!allocated) allocate();
// check if args correct
@ -185,16 +186,27 @@ double PairSpinExchange::init_one(int i, int j)
return cut_spin_exchange_global;
}
/* ----------------------------------------------------------------------
extract the larger cutoff
------------------------------------------------------------------------- */
void *PairSpinExchange::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"cut") == 0) return (void *) &cut_spin_exchange_global;
return NULL;
}
/* ---------------------------------------------------------------------- */
void PairSpinExchange::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double evdwl,ecoul;
double xi[3], rij[3];
double xi[3], rij[3], eij[3];
double fi[3], fmi[3];
double cut_ex_2,cut_spin_exchange_global2;
double rsq,rd,inorm;
double rsq, inorm;
int *ilist,*jlist,*numneigh,**firstneigh;
double spi[3],spj[3];
@ -202,7 +214,7 @@ void PairSpinExchange::compute(int eflag, int vflag)
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
cut_spin_exchange_global2 = cut_spin_exchange_global*cut_spin_exchange_global;
double **x = atom->x;
double **f = atom->f;
double **fm = atom->fm;
@ -229,12 +241,14 @@ void PairSpinExchange::compute(int eflag, int vflag)
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
itype = type[i];
// loop on neighbors
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
@ -248,22 +262,19 @@ void PairSpinExchange::compute(int eflag, int vflag)
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
inorm = 1.0/sqrt(rsq);
rij[0] *= inorm;
rij[1] *= inorm;
rij[2] *= inorm;
itype = type[i];
jtype = type[j];
// compute exchange interaction
eij[0] = inorm*rij[0];
eij[1] = inorm*rij[1];
eij[2] = inorm*rij[2];
cut_ex_2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
// compute exchange interaction
if (rsq <= cut_ex_2) {
compute_exchange(i,j,rsq,fmi,spi,spj);
if (lattice_flag) {
compute_exchange_mech(i,j,rsq,rij,fi,spi,spj);
compute_exchange(i,j,rsq,fmi,spi,spj);
if (lattice_flag) {
compute_exchange_mech(i,j,rsq,eij,fi,spi,spj);
}
}
@ -285,8 +296,9 @@ void PairSpinExchange::compute(int eflag, int vflag)
if (rsq <= cut_ex_2) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
} else evdwl = 0.0;
}
}
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]);
@ -297,6 +309,71 @@ void PairSpinExchange::compute(int eflag, int vflag)
}
/* ----------------------------------------------------------------------
update the pair interactions fmi acting on the spin ii
------------------------------------------------------------------------- */
void PairSpinExchange::compute_single_pair(int ii, double fmi[3])
{
const int nlocal = atom->nlocal;
int *type = atom->type;
double **x = atom->x;
double **sp = atom->sp;
double local_cut2;
double xi[3], rij[3];
double spi[3], spj[3];
int iexchange, idmi, ineel, ime;
int i,j,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
i = ilist[ii];
itype = type[i];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
local_cut2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
if (rsq <= local_cut2) {
compute_exchange(i,j,rsq,fmi,spi,spj);
}
}
}
/* ----------------------------------------------------------------------
compute the magnetic exchange interaction between spin i and spin j
------------------------------------------------------------------------- */
@ -321,6 +398,7 @@ void PairSpinExchange::compute_exchange(int i, int j, double rsq, double fmi[3],
fmi[0] += Jex*spj[0];
fmi[1] += Jex*spj[1];
fmi[2] += Jex*spj[2];
}
}

13
src/SPIN/pair_spin_exchange.h
View File

@ -23,7 +23,7 @@
#ifdef PAIR_CLASS
PairStyle(pair/spin/exchange,PairSpinExchange)
PairStyle(spin/exchange,PairSpinExchange)
#else
@ -42,8 +42,11 @@ class PairSpinExchange : public PairSpin {
void coeff(int, char **);
void init_style();
double init_one(int, int);
void *extract(const char *, int &);
void compute(int, int);
void compute_single_pair(int, double *);
void compute_exchange(int, int, double, double *, double *, double *);
void compute_exchange_mech(int, int, double, double *, double *, double *, double *);
@ -52,13 +55,17 @@ class PairSpinExchange : public PairSpin {
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
double cut_spin_exchange_global; // global neel cutoff distance
double cut_spin_exchange_global; // global exchange cutoff distance
protected:
double **J1_mag, **J1_mech; // exchange coeffs Jij
double **J1_mag; // exchange coeffs in eV
double **J1_mech; // mech exchange coeffs in
double **J2, **J3; // J1 in eV, J2 adim, J3 in Ang
double **cut_spin_exchange; // cutoff distance exchange
int lattice_flag; // flag for mech force computation
class FixNVESpin *lockfixnvespin; // ptr to FixNVESpin for setups
void allocate();
};

105
src/SPIN/pair_spin_me.cpp
View File

@ -163,16 +163,6 @@ void PairSpinMe::init_style()
}
if (ifix == modify->nfix)
error->all(FLERR,"pair/spin style requires nve/spin");
// get the lattice_flag from nve/spin
for (int i = 0; i < modify->nfix; i++) {
if (strcmp(modify->fix[i]->style,"nve/spin") == 0) {
lockfixnvespin = (FixNVESpin *) modify->fix[i];
lattice_flag = lockfixnvespin->lattice_flag;
}
}
}
/* ----------------------------------------------------------------------
@ -187,6 +177,18 @@ double PairSpinMe::init_one(int i, int j)
return cut_spin_me_global;
}
/* ----------------------------------------------------------------------
extract the larger cutoff
------------------------------------------------------------------------- */
void *PairSpinMe::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"cut") == 0) return (void *) &cut_spin_me_global;
return NULL;
}
/* ---------------------------------------------------------------------- */
void PairSpinMe::compute(int eflag, int vflag)
@ -197,14 +199,12 @@ void PairSpinMe::compute(int eflag, int vflag)
double spi[3], spj[3];
double fi[3], fj[3];
double fmi[3], fmj[3];
double cut_me_2, cut_spin_me_global2;
double rsq, rd, inorm;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
cut_spin_me_global2 = cut_spin_me_global*cut_spin_me_global;
double **x = atom->x;
double **f = atom->f;
@ -265,13 +265,8 @@ void PairSpinMe::compute(int eflag, int vflag)
// compute me interaction
cut_me_2 = cut_spin_me[itype][jtype]*cut_spin_me[itype][jtype];
if (rsq <= cut_me_2){
compute_me(i,j,rsq,rij,fmi,spi,spj);
if (lattice_flag) {
compute_me_mech(i,j,fi,spi,spj);
}
}
compute_me(i,j,rsq,rij,fmi,spi,spj);
compute_me_mech(i,j,fi,spi,spj);
f[i][0] += fi[0];
f[i][1] += fi[1];
@ -288,11 +283,9 @@ void PairSpinMe::compute(int eflag, int vflag)
}
if (eflag) {
if (rsq <= cut_me_2) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
} else evdwl = 0.0;
}
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]);
@ -305,6 +298,70 @@ void PairSpinMe::compute(int eflag, int vflag)
/* ---------------------------------------------------------------------- */
void PairSpinMe::compute_single_pair(int ii, double fmi[3])
{
const int nlocal = atom->nlocal;
int *type = atom->type;
double **x = atom->x;
double **sp = atom->sp;
double xi[3], rij[3], eij[3];
double spi[3], spj[3];
int iexchange, idmi, ineel, ime;
int i,j,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq, rd, inorm;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
i = ilist[ii];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
eij[0] = eij[1] = eij[2] = 0.0;
jlist = firstneigh[i];
jnum = numneigh[i];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
itype = type[ii];
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
inorm = 1.0/sqrt(rsq);
eij[0] = inorm*rij[0];
eij[1] = inorm*rij[1];
eij[2] = inorm*rij[2];
compute_me(i,j,rsq,eij,fmi,spi,spj);
}
}
/* ---------------------------------------------------------------------- */
void PairSpinMe::compute_me(int i, int j, double rsq, double eij[3], double fmi[3], double spi[3], double spj[3])
{
int *type = atom->type;

5
src/SPIN/pair_spin_me.h
View File

@ -23,7 +23,7 @@
#ifdef PAIR_CLASS
PairStyle(pair/spin/me,PairSpinMe)
PairStyle(spin/me,PairSpinMe)
#else
@ -42,8 +42,11 @@ class PairSpinMe : public PairSpin {
void coeff(int, char **);
void init_style();
double init_one(int, int);
void *extract(const char *, int &);
void compute(int, int);
void compute_single_pair(int, double *);
void compute_me(int, int, double, double *, double *, double *, double *);
void compute_me_mech(int, int, double *, double *, double *);

103
src/SPIN/pair_spin_neel.cpp
View File

@ -174,16 +174,6 @@ void PairSpinNeel::init_style()
}
if (ifix == modify->nfix)
error->all(FLERR,"pair/spin style requires nve/spin");
// get the lattice_flag from nve/spin
for (int i = 0; i < modify->nfix; i++) {
if (strcmp(modify->fix[i]->style,"nve/spin") == 0) {
lockfixnvespin = (FixNVESpin *) modify->fix[i];
lattice_flag = lockfixnvespin->lattice_flag;
}
}
}
/* ----------------------------------------------------------------------
@ -198,6 +188,17 @@ double PairSpinNeel::init_one(int i, int j)
return cut_spin_neel_global;
}
/* ----------------------------------------------------------------------
extract the larger cutoff
------------------------------------------------------------------------- */
void *PairSpinNeel::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"cut") == 0) return (void *) &cut_spin_neel_global;
return NULL;
}
/* ---------------------------------------------------------------------- */
void PairSpinNeel::compute(int eflag, int vflag)
@ -207,14 +208,12 @@ void PairSpinNeel::compute(int eflag, int vflag)
double xi[3], rij[3], eij[3];
double spi[3], spj[3];
double fi[3], fmi[3];
double cut_spin_neel_2, cut_spin_neel_global2;
double rsq, rd, inorm;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
cut_spin_neel_global2 = cut_spin_neel_global*cut_spin_neel_global;
double **x = atom->x;
double **f = atom->f;
@ -273,14 +272,9 @@ void PairSpinNeel::compute(int eflag, int vflag)
// compute magnetic and mechanical components of soc_neel
cut_spin_neel_2 = cut_spin_neel[itype][jtype]*cut_spin_neel[itype][jtype];
if (rsq <= cut_spin_neel_2) {
compute_neel(i,j,rsq,eij,fmi,spi,spj);
if (lattice_flag) {
compute_neel_mech(i,j,rsq,eij,fi,spi,spj);
}
}
compute_neel(i,j,rsq,eij,fmi,spi,spj);
compute_neel_mech(i,j,rsq,eij,fi,spi,spj);
f[i][0] += fi[0];
f[i][1] += fi[1];
f[i][2] += fi[2];
@ -296,11 +290,9 @@ void PairSpinNeel::compute(int eflag, int vflag)
}
if (eflag) {
if (rsq <= cut_spin_neel_2) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
} else evdwl = 0.0;
}
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fi[0],fi[1],fi[2],rij[0],rij[1],rij[2]);
@ -313,6 +305,67 @@ void PairSpinNeel::compute(int eflag, int vflag)
/* ---------------------------------------------------------------------- */
void PairSpinNeel::compute_single_pair(int ii, double fmi[3])
{
const int nlocal = atom->nlocal;
int *type = atom->type;
double **x = atom->x;
double **sp = atom->sp;
double xi[3], rij[3], eij[3];
double spi[3], spj[3];
int iexchange, idmi, ineel, ime;
int i,j,jj,inum,jnum,itype,jtype;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq, rd, inorm;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
i = ilist[ii];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
eij[0] = eij[1] = eij[2] = 0.0;
jlist = firstneigh[i];
jnum = numneigh[i];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
itype = type[ii];
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
rij[0] = x[j][0] - xi[0];
rij[1] = x[j][1] - xi[1];
rij[2] = x[j][2] - xi[2];
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
inorm = 1.0/sqrt(rsq);
compute_neel(i,j,rsq,eij,fmi,spi,spj);
}
}
/* ---------------------------------------------------------------------- */
void PairSpinNeel::compute_neel(int i, int j, double rsq, double eij[3], double fmi[3], double spi[3], double spj[3])
{
int *type = atom->type;

5
src/SPIN/pair_spin_neel.h
View File

@ -23,7 +23,7 @@
#ifdef PAIR_CLASS
PairStyle(pair/spin/neel,PairSpinNeel)
PairStyle(spin/neel,PairSpinNeel)
#else
@ -42,8 +42,11 @@ class PairSpinNeel : public PairSpin {
void coeff(int, char **);
void init_style();
double init_one(int, int);
void *extract(const char *, int &);
void compute(int, int);
void compute_single_pair(int, double *);
void compute_neel(int, int, double, double *, double *, double *, double *);
void compute_neel_mech(int, int, double, double *, double *, double *, double *);

1
src/pair_hybrid.h
View File

@ -31,7 +31,6 @@ class PairHybrid : public Pair {
friend class FixOMP;
friend class Force;
friend class Respa;
friend class PairSpin;
friend class Info;
public:
PairHybrid(class LAMMPS *);