diff --git a/doc/src/Eqs/pair_spin_soc_dmi_interaction.jpg b/doc/src/Eqs/pair_spin_dmi_interaction.jpg similarity index 100% rename from doc/src/Eqs/pair_spin_soc_dmi_interaction.jpg rename to doc/src/Eqs/pair_spin_dmi_interaction.jpg diff --git a/doc/src/Eqs/pair_spin_soc_dmi_interaction.tex b/doc/src/Eqs/pair_spin_dmi_interaction.tex similarity index 100% rename from doc/src/Eqs/pair_spin_soc_dmi_interaction.tex rename to doc/src/Eqs/pair_spin_dmi_interaction.tex diff --git a/doc/src/fix_precession_spin.txt b/doc/src/fix_precession_spin.txt index 9ef7a7230c..4133d7dd57 100644 --- a/doc/src/fix_precession_spin.txt +++ b/doc/src/fix_precession_spin.txt @@ -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. diff --git a/doc/src/pair_spin_dmi.txt b/doc/src/pair_spin_dmi.txt index e2987ce1df..9ad74d567e 100644 --- a/doc/src/pair_spin_dmi.txt +++ b/doc/src/pair_spin_dmi.txt @@ -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) diff --git a/doc/src/pair_spin_exchange.txt b/doc/src/pair_spin_exchange.txt index e6938afd40..90a75705c7 100644 --- a/doc/src/pair_spin_exchange.txt +++ b/doc/src/pair_spin_exchange.txt @@ -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) diff --git a/doc/src/pair_spin_me.txt b/doc/src/pair_spin_me.txt index bc4b6b840b..54c6d43300 100644 --- a/doc/src/pair_spin_me.txt +++ b/doc/src/pair_spin_me.txt @@ -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: diff --git a/doc/src/pair_spin_neel.txt b/doc/src/pair_spin_neel.txt index 84893fecb1..5c82b6ef6f 100644 --- a/doc/src/pair_spin_neel.txt +++ b/doc/src/pair_spin_neel.txt @@ -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) diff --git a/examples/SPIN/bfo/in.spin.bfo b/examples/SPIN/bfo/in.spin.bfo index dfeb48eb08..55cc53446d 100644 --- a/examples/SPIN/bfo/in.spin.bfo +++ b/examples/SPIN/bfo/in.spin.bfo @@ -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 diff --git a/examples/SPIN/cobalt_fcc/in.spin.cobalt_fcc b/examples/SPIN/cobalt_fcc/in.spin.cobalt_fcc index 45d9d983bd..c1d7e4f903 100644 --- a/examples/SPIN/cobalt_fcc/in.spin.cobalt_fcc +++ b/examples/SPIN/cobalt_fcc/in.spin.cobalt_fcc @@ -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 diff --git a/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp b/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp index 24d19a2252..127e24bf2b 100644 --- a/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp +++ b/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp @@ -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 diff --git a/examples/SPIN/iron/in.spin.iron b/examples/SPIN/iron/in.spin.iron index 83301309e0..0f40a5af21 100644 --- a/examples/SPIN/iron/in.spin.iron +++ b/examples/SPIN/iron/in.spin.iron @@ -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 diff --git a/examples/SPIN/nickel/in.spin.nickel b/examples/SPIN/nickel/in.spin.nickel index 0b1faa5319..2fa16b8a1b 100644 --- a/examples/SPIN/nickel/in.spin.nickel +++ b/examples/SPIN/nickel/in.spin.nickel @@ -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 diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp index c4dca269ad..2d3f73cd77 100644 --- a/src/SPIN/fix_nve_spin.cpp +++ b/src/SPIN/fix_nve_spin.cpp @@ -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; ipair_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; ipair_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; diff --git a/src/SPIN/fix_nve_spin.h b/src/SPIN/fix_nve_spin.h index 89824a9bd1..686c391299 100644 --- a/src/SPIN/fix_nve_spin.h +++ b/src/SPIN/fix_nve_spin.h @@ -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; diff --git a/src/SPIN/fix_precession_spin.cpp b/src/SPIN/fix_precession_spin.cpp index c5fde83fe6..67984a6009 100644 --- a/src/SPIN/fix_precession_spin.cpp +++ b/src/SPIN/fix_precession_spin.cpp @@ -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; diff --git a/src/SPIN/fix_precession_spin.h b/src/SPIN/fix_precession_spin.h index 71155ddd2b..5e1047ff67 100644 --- a/src/SPIN/fix_precession_spin.h +++ b/src/SPIN/fix_precession_spin.h @@ -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 diff --git a/src/SPIN/pair_spin.cpp b/src/SPIN/pair_spin.cpp index a40f2a86b7..d0df3868df 100755 --- a/src/SPIN/pair_spin.cpp +++ b/src/SPIN/pair_spin.cpp @@ -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; - } - } } diff --git a/src/SPIN/pair_spin.h b/src/SPIN/pair_spin.h index 8fe189801d..e71f2eb117 100755 --- a/src/SPIN/pair_spin.h +++ b/src/SPIN/pair_spin.h @@ -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() {} }; } diff --git a/src/SPIN/pair_spin_dmi.cpp b/src/SPIN/pair_spin_dmi.cpp index 5c1ab13065..99b41f77c8 100755 --- a/src/SPIN/pair_spin_dmi.cpp +++ b/src/SPIN/pair_spin_dmi.cpp @@ -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; } } diff --git a/src/SPIN/pair_spin_dmi.h b/src/SPIN/pair_spin_dmi.h index ff6faa2a46..da9f8d1494 100755 --- a/src/SPIN/pair_spin_dmi.h +++ b/src/SPIN/pair_spin_dmi.h @@ -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 *); diff --git a/src/SPIN/pair_spin_exchange.cpp b/src/SPIN/pair_spin_exchange.cpp index c899fba40a..145d49fe85 100755 --- a/src/SPIN/pair_spin_exchange.cpp +++ b/src/SPIN/pair_spin_exchange.cpp @@ -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]; + } } diff --git a/src/SPIN/pair_spin_exchange.h b/src/SPIN/pair_spin_exchange.h index ab22ef9255..66750743bb 100755 --- a/src/SPIN/pair_spin_exchange.h +++ b/src/SPIN/pair_spin_exchange.h @@ -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(); }; diff --git a/src/SPIN/pair_spin_me.cpp b/src/SPIN/pair_spin_me.cpp index e617470c8a..ee32127c23 100755 --- a/src/SPIN/pair_spin_me.cpp +++ b/src/SPIN/pair_spin_me.cpp @@ -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; diff --git a/src/SPIN/pair_spin_me.h b/src/SPIN/pair_spin_me.h index 05e6719710..58ea6b3eda 100755 --- a/src/SPIN/pair_spin_me.h +++ b/src/SPIN/pair_spin_me.h @@ -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 *); diff --git a/src/SPIN/pair_spin_neel.cpp b/src/SPIN/pair_spin_neel.cpp index 205ee736ed..a34931af85 100755 --- a/src/SPIN/pair_spin_neel.cpp +++ b/src/SPIN/pair_spin_neel.cpp @@ -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; diff --git a/src/SPIN/pair_spin_neel.h b/src/SPIN/pair_spin_neel.h index 970844a2a7..114e1843de 100755 --- a/src/SPIN/pair_spin_neel.h +++ b/src/SPIN/pair_spin_neel.h @@ -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 *); diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h index c01f57026b..202847b028 100644 --- a/src/pair_hybrid.h +++ b/src/pair_hybrid.h @@ -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 *);