diff --git a/doc/src/Eqs/pair_spin_dmi_forces.jpg b/doc/src/Eqs/pair_spin_dmi_forces.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..fba6a91cbf6fdb6ee6736073a114d08c907d6d3c
Binary files /dev/null and b/doc/src/Eqs/pair_spin_dmi_forces.jpg differ
diff --git a/doc/src/Eqs/pair_spin_dmi_forces.tex b/doc/src/Eqs/pair_spin_dmi_forces.tex
new file mode 100644
index 0000000000000000000000000000000000000000..1c0c246db4c70852212f86c9a7ab7e6b518c323c
--- /dev/null
+++ b/doc/src/Eqs/pair_spin_dmi_forces.tex
@@ -0,0 +1,14 @@
+\documentclass[preview]{standalone}
+\usepackage{varwidth}
+\usepackage[utf8x]{inputenc}
+\usepackage{amsmath,amssymb,amsthm,bm}
+\begin{document}
+\begin{varwidth}{50in}
+ \begin{equation}
+ \vec{\omega}_i = -\frac{1}{\hbar} \sum_{j}^{Neighb} \vec{s}_{j}\times \left(\vec{e}_{ij}\times \vec{D} \right)
+ ~~{\rm and}~~
+ \vec{F}_i = -\sum_{j}^{Neighb} \frac{1}{r_{ij}} \vec{D} \times \left( \vec{s}_{i}\times \vec{s}_{j} \right)
+ , \nonumber
+ \end{equation}
+\end{varwidth}
+\end{document}
diff --git a/doc/src/Eqs/pair_spin_dmi_interaction.jpg b/doc/src/Eqs/pair_spin_dmi_interaction.jpg
index 1d15b2199add75eefecb424f75ea2dbd363a1afe..3eb24c67e3b26f97c362fd05982ead5fb0f5c150 100644
Binary files a/doc/src/Eqs/pair_spin_dmi_interaction.jpg and b/doc/src/Eqs/pair_spin_dmi_interaction.jpg differ
diff --git a/doc/src/Eqs/pair_spin_dmi_interaction.tex b/doc/src/Eqs/pair_spin_dmi_interaction.tex
index 5a5a776f030088efe17fcbe5e131e4686e8de2a3..79f63a333a6c823258f5d9012fa40561e00fb9cb 100644
--- a/doc/src/Eqs/pair_spin_dmi_interaction.tex
+++ b/doc/src/Eqs/pair_spin_dmi_interaction.tex
@@ -5,7 +5,7 @@
\begin{document}
\begin{varwidth}{50in}
\begin{equation}
- \bm{H}_{dm} = -\sum_{{ i,j}=1,i\neq j}^{N}
+ \bm{H}_{dm} = \sum_{{ i,j}=1,i\neq j}^{N}
\left( \vec{e}_{ij} \times \vec{D} \right)
\cdot\left(\vec{s}_{i}\times \vec{s}_{j}\right),
\nonumber
diff --git a/doc/src/Eqs/pair_spin_exchange_forces.jpg b/doc/src/Eqs/pair_spin_exchange_forces.jpg
index 2b0469bf4db84885fd8b81679cdd9571ba9e3574..b312a9ccdaaf44b7647f2c444c7d79a740cea88c 100644
Binary files a/doc/src/Eqs/pair_spin_exchange_forces.jpg and b/doc/src/Eqs/pair_spin_exchange_forces.jpg differ
diff --git a/doc/src/Eqs/pair_spin_exchange_forces.tex b/doc/src/Eqs/pair_spin_exchange_forces.tex
index 088696d5efbc8f44b7b9c5f87ca8984984927f2f..ac5ef682f3ca9d9532ff6d258cb1012529058636 100644
--- a/doc/src/Eqs/pair_spin_exchange_forces.tex
+++ b/doc/src/Eqs/pair_spin_exchange_forces.tex
@@ -5,10 +5,12 @@
\begin{document}
\begin{varwidth}{50in}
\begin{equation}
- \vec{F}^{i} = \sum_{j}^{Neighbor} \frac{\partial {J} \left(r_{ij} \right)}{
- \partial r_{ij}} \left( \vec{s}_{i}\cdot \vec{s}_{j} \right) \vec{r}_{ij}
- ~~{\rm and}~~ \vec{\omega}^{i} = \frac{1}{\hbar} \sum_{j}^{Neighbor} {J}
- \left(r_{ij} \right)\,\vec{s}_{j} \nonumber
+ \vec{\omega}_{i} = \frac{1}{\hbar} \sum_{j}^{Neighb} {J}
+ \left(r_{ij} \right)\,\vec{s}_{j}
+ ~~{\rm and}~~
+ \vec{F}_{i} = \sum_{j}^{Neighb} \frac{\partial {J} \left(r_{ij} \right)}{
+ \partial r_{ij}} \left( \vec{s}_{i}\cdot \vec{s}_{j} \right) \vec{e}_{ij}
+ \nonumber
\end{equation}
\end{varwidth}
\end{document}
diff --git a/doc/src/Eqs/pair_spin_exchange_interaction.jpg b/doc/src/Eqs/pair_spin_exchange_interaction.jpg
index d51524d27ca413ec4ded189107c2219b061f7d7e..c70d8a6554003e37472d648c1944a13ee9d3f859 100644
Binary files a/doc/src/Eqs/pair_spin_exchange_interaction.jpg and b/doc/src/Eqs/pair_spin_exchange_interaction.jpg differ
diff --git a/doc/src/Eqs/pair_spin_exchange_interaction.tex b/doc/src/Eqs/pair_spin_exchange_interaction.tex
index 6e598f75ac04b4318c8d5679c032a75c5e82aeec..f20b3e5740ce2db9cb86a70910ca2b80f263fe4f 100644
--- a/doc/src/Eqs/pair_spin_exchange_interaction.tex
+++ b/doc/src/Eqs/pair_spin_exchange_interaction.tex
@@ -5,7 +5,7 @@
\begin{document}
\begin{varwidth}{50in}
\begin{equation}
- \bm{H}_{exchange} ~=~ -\sum_{i,j,i\neq j}^{N} {J} \left(r_{ij} \right)\, \vec{s}_{i}\cdot \vec{s}_{j} \nonumber
+ \bm{H}_{ex} ~=~ -\sum_{i,j,i\neq j}^{N} {J} \left(r_{ij} \right)\, \vec{s}_{i}\cdot \vec{s}_{j} \nonumber
\end{equation}
\end{varwidth}
\end{document}
diff --git a/doc/src/pair_spin_dmi.txt b/doc/src/pair_spin_dmi.txt
index 24877906f3de032b60b88db2cdc11ca3d7140f49..a1d0cb2265b9f58a46d33072c9e819acb364c00f 100644
--- a/doc/src/pair_spin_dmi.txt
+++ b/doc/src/pair_spin_dmi.txt
@@ -25,24 +25,46 @@ pair_coeff 1 2 dmi 4.0 0.00109 0.0 0.0 1.0 :pre
[Description:]
Style {spin/dmi} computes the Dzyaloshinskii-Moriya (DM) interaction
-between pairs of magnetic spins:
+between pairs of magnetic spins.
+According to the expression reported in "(Rohart)"_#Rohart, one has
+the following DM energy:
:c,image(Eqs/pair_spin_dmi_interaction.jpg)
where si and sj are two neighboring magnetic spins of two particles,
-eij = (ri - rj)/|ri-rj| is the normalized separation vector between the
-two particles, and D is the DM vector defining the intensity and the
-sign of the interaction.
+eij = (ri - rj)/|ri-rj| is the unit vector between sites i and j,
+and D is the DM vector defining the intensity (in eV) and the direction
+of the interaction.
-Examples and more explanations about this interaction and its parametrization are
-reported in "(Tranchida)"_#Tranchida5.
+In "(Rohart)"_#Rohart, D is defined as the direction normal to the film oriented
+from the high spin-orbit layer to the magnetic ultrathin film.
-From this DM interaction, each spin i will be submitted to a magnetic torque
-omega and its associated atom to a force F (for spin-lattice calculations only).
+The application of a spin-lattice Poisson bracket to this energy (as described
+in "(Tranchida)"_#Tranchida5) allows to derive a magnetic torque omega, and a
+mechanical force F (for spin-lattice calculations only) for each magnetic
+particle i:
+
+:c,image(Eqs/pair_spin_dmi_forces.jpg)
More details about the derivation of these torques/forces are reported in
"(Tranchida)"_#Tranchida5.
+For the {spin/dmi} pair style, the following coefficients must be defined for
+each pair of atoms types via the "pair_coeff"_pair_coeff.html command as in
+the examples above, or in the data file or restart files read by the
+"read_data"_read_data.html or "read_restart"_read_restart.html commands, and
+set in the following order:
+
+rc (distance units)
+|D| (energy units)
+Dx, Dy, Dz (direction of D) :ul
+
+Note that rc is the radius cutoff of the considered DM interaction, |D| is
+the norm of the DM vector (in eV), and Dx, Dy and Dz define its direction.
+
+None of those coefficients is optional. If not specified, the {spin/dmi}
+pair style cannot be used.
+
:line
[Restrictions:]
@@ -61,6 +83,9 @@ See the "Making LAMMPS"_Section_start.html#start_3 section for more info.
:line
+:link(Rohart)
+[(Rohart)] Rohart and Thiaville,
+Physical Review B, 88(18), 184422. (2013).
:link(Tranchida5)
[(Tranchida)] Tranchida, Plimpton, Thibaudeau and Thompson,
Journal of Computational Physics, (2018).
diff --git a/doc/src/pair_spin_exchange.txt b/doc/src/pair_spin_exchange.txt
index ad3357cb5eebcfab4c12567752fbf95711ecb6a1..a4f37a98eb4e9a205f8c08842d70ef3b6035276d 100644
--- a/doc/src/pair_spin_exchange.txt
+++ b/doc/src/pair_spin_exchange.txt
@@ -32,18 +32,15 @@ pairs of magnetic spins:
where si and sj are two neighboring magnetic spins of two particles,
rij = ri - rj is the inter-atomic distance between the two particles,
and J(rij) is a function defining the intensity and the sign of the exchange
-interaction.
-
-This function is defined as:
+interaction for different neighboring shells. This function is defined as:
:c,image(Eqs/pair_spin_exchange_function.jpg)
where a, b and d are the three constant coefficients defined in the associated
-"pair_coeff" command.
+"pair_coeff" command (see below for more explanations).
The coefficients a, b, and d need to be fitted so that the function above matches with
the value of the exchange interaction for the N neighbor shells taken into account.
-
Examples and more explanations about this function and its parametrization are reported
in "(Tranchida)"_#Tranchida3.
@@ -54,11 +51,30 @@ such as:
:c,image(Eqs/pair_spin_exchange_forces.jpg)
-with h the Planck constant (in metal units).
+with h the Planck constant (in metal units), and eij = (ri - rj)/|ri-rj| the unit
+vector between sites i and j.
More details about the derivation of these torques/forces are reported in
"(Tranchida)"_#Tranchida3.
+For the {spin/exchange} pair style, the following coefficients must be defined
+for each pair of atoms types via the "pair_coeff"_pair_coeff.html command as in
+the examples above, or in the data file or restart files read by the
+"read_data"_read_data.html or "read_restart"_read_restart.html commands, and
+set in the following order:
+
+rc (distance units)
+a (energy units)
+b (adim parameter)
+d (distance units) :ul
+
+Note that rc is the radius cutoff of the considered exchange interaction,
+and a, b and d are the three coefficients performing the parametrization
+of the function J(rij) defined above.
+
+None of those coefficients is optional. If not specified, the
+{spin/exchange} pair style cannot be used.
+
:line
[Restrictions:]
diff --git a/examples/SPIN/bfo/in.spin.bfo b/examples/SPIN/bfo/in.spin.bfo
index 2442b12b72294d2b138831565a868a65c5395af5..de23ba87ba2b825c129b88fd55fcfcb33703eb01 100644
--- a/examples/SPIN/bfo/in.spin.bfo
+++ b/examples/SPIN/bfo/in.spin.bfo
@@ -21,9 +21,11 @@ mass 1 1.0
set group all spin/random 11 2.50
-pair_style hybrid/overlay spin/exchange 6.0 spin/magelec 4.5
+#pair_style hybrid/overlay spin/exchange 6.0 spin/magelec 4.5
+pair_style hybrid/overlay spin/exchange 6.0 spin/magelec 4.5 spin/dmi 4.5
pair_coeff * * spin/exchange exchange 6.0 -0.01575 0.0 1.965
pair_coeff * * spin/magelec magelec 4.5 0.000109 1.0 1.0 1.0
+pair_coeff * * spin/dmi dmi 4.5 0.00005 1.0 1.0 1.0
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
@@ -44,10 +46,11 @@ 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 50
+#thermo_style custom step time v_magnorm v_emag temp etotal
+thermo_style custom step time v_magnorm pe ke v_emag temp etotal
+thermo 10
compute outsp all property/atom spx spy spz sp fmx fmy fmz
dump 100 all custom 1 dump_bfo.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3]
-run 5000
+run 2000
diff --git a/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp b/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp
index 4a42ec419aeb1a02603d86cb748d8964b439ea06..35aa1df86c4a2c2bdfff40f3e35f032939311b90 100644
--- a/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp
+++ b/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp
@@ -19,8 +19,8 @@ create_atoms 1 box
mass 1 58.93
-#set group all spin/random 31 1.72
-set group all spin 1.72 0.0 0.0 1.0
+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 spin/exchange 4.0 spin/neel 4.0
@@ -29,11 +29,11 @@ pair_coeff * * eam/alloy Co_PurjaPun_2012.eam.alloy Co
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
-fix 1 all precession/spin zeeman 1.0 0.0 0.0 1.0
+#fix 1 all precession/spin zeeman 1.0 0.0 0.0 1.0
+fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
@@ -56,4 +56,4 @@ thermo 10
compute outsp all property/atom spx spy spz sp fmx fmy fmz
dump 100 all custom 1 dump_cobalt_hcp.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3]
-run 2000
+run 20000
diff --git a/src/SPIN/atom_vec_spin.cpp b/src/SPIN/atom_vec_spin.cpp
index 8b47eff624c27535f6d1cb297a792de3aadf99ea..6460a6185fb72d0635b623dbd7e8a2e7a2d27ba6 100644
--- a/src/SPIN/atom_vec_spin.cpp
+++ b/src/SPIN/atom_vec_spin.cpp
@@ -816,9 +816,9 @@ void AtomVecSpin::data_atom(double *coord, imageint imagetmp, char **values)
x[nlocal][2] = coord[2];
sp[nlocal][3] = atof(values[2]);
- sp[nlocal][0] = atof(values[5]);
- sp[nlocal][1] = atof(values[6]);
- sp[nlocal][2] = atof(values[7]);
+ sp[nlocal][0] = atof(values[6]);
+ sp[nlocal][1] = atof(values[7]);
+ sp[nlocal][2] = atof(values[8]);
double inorm = 1.0/sqrt(sp[nlocal][0]*sp[nlocal][0] +
sp[nlocal][1]*sp[nlocal][1] +
sp[nlocal][2]*sp[nlocal][2]);
diff --git a/src/SPIN/compute_spin.cpp b/src/SPIN/compute_spin.cpp
index b67f62d53d61cdf8fce10d1a09a07bb59939aa92..dc16190c98763053d76fb84e5a423eca570f860d 100644
--- a/src/SPIN/compute_spin.cpp
+++ b/src/SPIN/compute_spin.cpp
@@ -105,16 +105,16 @@ void ComputeSpin::compute_vector()
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (atom->sp_flag) {
- mag[0] += sp[i][0];
- mag[1] += sp[i][1];
- mag[2] += sp[i][2];
- magenergy -= (sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]);
- tx = sp[i][1]*fm[i][2]-sp[i][2]*fm[i][1];
- ty = sp[i][2]*fm[i][0]-sp[i][0]*fm[i][2];
- tz = sp[i][0]*fm[i][1]-sp[i][1]*fm[i][0];
- tempnum += tx*tx+ty*ty+tz*tz;
- tempdenom += sp[i][0]*fm[i][0]+fm[i][1]*sp[i][1]+sp[i][2]*fm[i][2];
- countsp++;
+ mag[0] += sp[i][0];
+ mag[1] += sp[i][1];
+ mag[2] += sp[i][2];
+ magenergy -= (sp[i][0]*fm[i][0] + sp[i][1]*fm[i][1] + sp[i][2]*fm[i][2]);
+ tx = sp[i][1]*fm[i][2]-sp[i][2]*fm[i][1];
+ ty = sp[i][2]*fm[i][0]-sp[i][0]*fm[i][2];
+ tz = sp[i][0]*fm[i][1]-sp[i][1]*fm[i][0];
+ tempnum += tx*tx+ty*ty+tz*tz;
+ tempdenom += sp[i][0]*fm[i][0]+fm[i][1]*sp[i][1]+sp[i][2]*fm[i][2];
+ countsp++;
}
}
else error->all(FLERR,"Compute compute/spin requires atom/spin style");
diff --git a/src/SPIN/pair_spin_dmi.cpp b/src/SPIN/pair_spin_dmi.cpp
index b792969c5dcd5bc4c5d238143de90455c8cf885d..08e2c63e7f6ba147c6062da17f92a68b2a5d090a 100644
--- a/src/SPIN/pair_spin_dmi.cpp
+++ b/src/SPIN/pair_spin_dmi.cpp
@@ -65,6 +65,9 @@ PairSpinDmi::~PairSpinDmi()
memory->destroy(v_dmx);
memory->destroy(v_dmy);
memory->destroy(v_dmz);
+ memory->destroy(vmech_dmx);
+ memory->destroy(vmech_dmy);
+ memory->destroy(vmech_dmz);
memory->destroy(cutsq);
}
}
@@ -118,7 +121,7 @@ void PairSpinDmi::coeff(int narg, char **arg)
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
const double rij = force->numeric(FLERR,arg[3]);
- const double dm = (force->numeric(FLERR,arg[4]))/hbar;
+ const double dm = (force->numeric(FLERR,arg[4]));
double dmx = force->numeric(FLERR,arg[5]);
double dmy = force->numeric(FLERR,arg[6]);
double dmz = force->numeric(FLERR,arg[7]);
@@ -133,9 +136,12 @@ void PairSpinDmi::coeff(int narg, char **arg)
for (int j = MAX(jlo,i); j <= jhi; j++) {
cut_spin_dmi[i][j] = rij;
DM[i][j] = dm;
- v_dmx[i][j] = dmx * dm;
- v_dmy[i][j] = dmy * dm;
- v_dmz[i][j] = dmz * dm;
+ v_dmx[i][j] = dmx * dm / hbar;
+ v_dmy[i][j] = dmy * dm / hbar;
+ v_dmz[i][j] = dmz * dm / hbar;
+ vmech_dmx[i][j] = dmx * dm;
+ vmech_dmy[i][j] = dmy * dm;
+ vmech_dmz[i][j] = dmz * dm;
setflag[i][j] = 1;
count++;
}
@@ -187,9 +193,17 @@ void PairSpinDmi::init_style()
double PairSpinDmi::init_one(int i, int j)
{
-
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+ DM[j][i] = DM[i][j];
+ v_dmx[j][i] = v_dmx[i][j];
+ v_dmy[j][i] = v_dmy[i][j];
+ v_dmz[j][i] = v_dmz[i][j];
+ vmech_dmx[j][i] = vmech_dmx[i][j];
+ vmech_dmy[j][i] = vmech_dmy[i][j];
+ vmech_dmz[j][i] = vmech_dmz[i][j];
+ cut_spin_dmi[j][i] = cut_spin_dmi[i][j];
+
return cut_spin_dmi_global;
}
@@ -210,7 +224,8 @@ void PairSpinDmi::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double evdwl, ecoul;
- double xi[3], rij[3], eij[3];
+ double xi[3], eij[3];
+ double delx,dely,delz;
double spi[3], spj[3];
double fi[3], fmi[3];
double local_cut2;
@@ -264,20 +279,17 @@ void PairSpinDmi::compute(int eflag, int vflag)
spj[2] = sp[j][2];
evdwl = 0.0;
-
fi[0] = fi[1] = fi[2] = 0.0;
fmi[0] = fmi[1] = fmi[2] = 0.0;
- rij[0] = rij[1] = rij[2] = 0.0;
- eij[0] = eij[1] = eij[2] = 0.0;
- 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];
+ delx = xi[0] - x[j][0];
+ dely = xi[1] - x[j][1];
+ delz = xi[2] - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
inorm = 1.0/sqrt(rsq);
- eij[0] = rij[0]*inorm;
- eij[1] = rij[1]*inorm;
- eij[2] = rij[2]*inorm;
+ eij[0] = -inorm*delx;
+ eij[1] = -inorm*dely;
+ eij[2] = -inorm*delz;
local_cut2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];
@@ -286,7 +298,7 @@ void PairSpinDmi::compute(int eflag, int vflag)
if (rsq <= local_cut2) {
compute_dmi(i,j,eij,fmi,spj);
if (lattice_flag) {
- compute_dmi_mech(fi);
+ compute_dmi_mech(i,j,rsq,eij,fi,spi,spj);
}
}
@@ -309,7 +321,7 @@ void PairSpinDmi::compute(int eflag, int vflag)
} 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]);
+ evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
}
}
@@ -325,8 +337,8 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3])
double **x = atom->x;
double **sp = atom->sp;
double local_cut2;
-
- double xi[3], rij[3], eij[3];
+ double xi[3], eij[3];
+ double delx,dely,delz;
double spj[3];
int i,j,jnum,itype,jtype;
@@ -358,14 +370,14 @@ void PairSpinDmi::compute_single_pair(int ii, double fmi[3])
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];
+ delx = xi[0] - x[j][0];
+ dely = xi[1] - x[j][1];
+ delz = xi[2] - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
inorm = 1.0/sqrt(rsq);
- eij[0] = rij[0]*inorm;
- eij[1] = rij[1]*inorm;
- eij[2] = rij[2]*inorm;
+ eij[0] = -inorm*delx;
+ eij[1] = -inorm*dely;
+ eij[2] = -inorm*delz;
local_cut2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];
@@ -390,23 +402,45 @@ void PairSpinDmi::compute_dmi(int i, int j, double eij[3], double fmi[3], double
jtype = type[j];
dmix = eij[1]*v_dmz[itype][jtype] - eij[2]*v_dmy[itype][jtype];
- dmiy = eij[2]*v_dmx[itype][jtype] - eij[2]*v_dmz[itype][jtype];
+ dmiy = eij[2]*v_dmx[itype][jtype] - eij[0]*v_dmz[itype][jtype];
dmiz = eij[0]*v_dmy[itype][jtype] - eij[1]*v_dmx[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);
}
/* ----------------------------------------------------------------------
compute the mechanical force due to the dmi interaction between atom i and atom j
------------------------------------------------------------------------- */
-void PairSpinDmi::compute_dmi_mech(double fi[3])
+void PairSpinDmi::compute_dmi_mech(int i, int j, double rsq, double eij[3],
+ double fi[3], double spi[3], double spj[3])
{
- fi[0] += 0.0;
- fi[1] += 0.0;
- fi[2] += 0.0;
+ int *type = atom->type;
+ int itype, jtype;
+ double dmix,dmiy,dmiz;
+ itype = type[i];
+ jtype = type[j];
+ double csx,csy,csz,cdmx,cdmy,cdmz,irij;
+
+ irij = 1.0/sqrt(rsq);
+
+ dmix = vmech_dmx[itype][jtype];
+ dmiy = vmech_dmy[itype][jtype];
+ dmiz = vmech_dmz[itype][jtype];
+
+ csx = (spi[1]*spj[2] - spi[2]*spj[1]);
+ csy = (spi[2]*spj[0] - spi[0]*spj[2]);
+ csz = (spi[0]*spj[1] - spi[1]*spj[0]);
+
+ cdmx = (dmiy*csz - dmiz*csy);
+ cdmy = (dmiz*csx - dmix*csz);
+ cdmz = (dmix*csy - dmiy*csz);
+
+ fi[0] += irij*cdmx;
+ fi[1] += irij*cdmy;
+ fi[2] += irij*cdmz;
}
/* ----------------------------------------------------------------------
@@ -428,6 +462,9 @@ void PairSpinDmi::allocate()
memory->create(v_dmx,n+1,n+1,"pair:DM_vector_x");
memory->create(v_dmy,n+1,n+1,"pair:DM_vector_y");
memory->create(v_dmz,n+1,n+1,"pair:DM_vector_z");
+ memory->create(vmech_dmx,n+1,n+1,"pair:DMmech_vector_x");
+ memory->create(vmech_dmy,n+1,n+1,"pair:DMmech_vector_y");
+ memory->create(vmech_dmz,n+1,n+1,"pair:DMmech_vector_z");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
@@ -451,6 +488,9 @@ void PairSpinDmi::write_restart(FILE *fp)
fwrite(&v_dmx[i][j],sizeof(double),1,fp);
fwrite(&v_dmy[i][j],sizeof(double),1,fp);
fwrite(&v_dmz[i][j],sizeof(double),1,fp);
+ fwrite(&vmech_dmx[i][j],sizeof(double),1,fp);
+ fwrite(&vmech_dmy[i][j],sizeof(double),1,fp);
+ fwrite(&vmech_dmz[i][j],sizeof(double),1,fp);
fwrite(&cut_spin_dmi[i][j],sizeof(double),1,fp);
}
}
@@ -478,12 +518,18 @@ void PairSpinDmi::read_restart(FILE *fp)
fread(&v_dmx[i][j],sizeof(double),1,fp);
fread(&v_dmy[i][j],sizeof(double),1,fp);
fread(&v_dmz[i][j],sizeof(double),1,fp);
+ fread(&vmech_dmx[i][j],sizeof(double),1,fp);
+ fread(&vmech_dmy[i][j],sizeof(double),1,fp);
+ fread(&vmech_dmz[i][j],sizeof(double),1,fp);
fread(&cut_spin_dmi[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&DM[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&v_dmx[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&v_dmy[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&v_dmz[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&vmech_dmx[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&vmech_dmy[i][j],1,MPI_DOUBLE,0,world);
+ MPI_Bcast(&vmech_dmz[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_spin_dmi[i][j],1,MPI_DOUBLE,0,world);
}
}
diff --git a/src/SPIN/pair_spin_dmi.h b/src/SPIN/pair_spin_dmi.h
index a309f0c8d52ade48d42cffc26254e5a5c6839988..68e42e879dfe15577d6bf73fe8f80a5593530dfd 100644
--- a/src/SPIN/pair_spin_dmi.h
+++ b/src/SPIN/pair_spin_dmi.h
@@ -38,22 +38,23 @@ class PairSpinDmi : public PairSpin {
void compute_single_pair(int, double *);
void compute_dmi(int, int, double *, double *, double *);
- void compute_dmi_mech(double *);
+ void compute_dmi_mech(int, int, double, double *, double *, double *, double *);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
- double cut_spin_dmi_global; // short range pair cutoff
+ double cut_spin_dmi_global; // short range pair cutoff
protected:
- double **DM; // dmi coeff in eV
- double **v_dmx, **v_dmy, **v_dmz; // dmi direction
- double **cut_spin_dmi; // cutoff distance dmi
+ double **DM; // dmi coeff in eV
+ double **v_dmx, **v_dmy, **v_dmz; // dmi direction
+ double **vmech_dmx, **vmech_dmy, **vmech_dmz; // dmi mech direction
+ double **cut_spin_dmi; // cutoff distance dmi
- int lattice_flag; // flag for mech force computation
- class FixNVESpin *lockfixnvespin; // ptr to FixNVESpin for setups
+ 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_exchange.cpp b/src/SPIN/pair_spin_exchange.cpp
index 1b7b36b6dbac6d6cb69ef566e1bd71132c3b4d3b..cc074bb97d9b9dd91d2bfcc8e378d7cd56a7ba5e 100644
--- a/src/SPIN/pair_spin_exchange.cpp
+++ b/src/SPIN/pair_spin_exchange.cpp
@@ -65,7 +65,7 @@ PairSpinExchange::~PairSpinExchange()
memory->destroy(J1_mech);
memory->destroy(J2);
memory->destroy(J3);
- memory->destroy(cutsq); // to be deleted
+ memory->destroy(cutsq); // to be implemented
}
}
@@ -134,8 +134,8 @@ void PairSpinExchange::coeff(int narg, char **arg)
count++;
}
}
- if (count == 0)
- error->all(FLERR,"Incorrect args in pair_style command");
+
+ if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
}
/* ----------------------------------------------------------------------
@@ -183,6 +183,12 @@ double PairSpinExchange::init_one(int i, int j)
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+ J1_mag[j][i] = J1_mag[i][j];
+ J1_mech[j][i] = J1_mech[i][j];
+ J2[j][i] = J2[i][j];
+ J3[j][i] = J3[i][j];
+ cut_spin_exchange[j][i] = cut_spin_exchange[i][j];
+
return cut_spin_exchange_global;
}
@@ -203,7 +209,8 @@ void PairSpinExchange::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double evdwl, ecoul;
- double xi[3], rij[3], eij[3];
+ double xi[3], eij[3];
+ double delx,dely,delz;
double spi[3], spj[3];
double fi[3], fmi[3];
double local_cut2;
@@ -255,18 +262,17 @@ void PairSpinExchange::compute(int eflag, int vflag)
spj[2] = sp[j][2];
evdwl = 0.0;
-
fi[0] = fi[1] = fi[2] = 0.0;
fmi[0] = fmi[1] = fmi[2] = 0.0;
- 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];
+ delx = xi[0] - x[j][0];
+ dely = xi[1] - x[j][1];
+ delz = xi[2] - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
inorm = 1.0/sqrt(rsq);
- eij[0] = inorm*rij[0];
- eij[1] = inorm*rij[1];
- eij[2] = inorm*rij[2];
+ eij[0] = -inorm*delx;
+ eij[1] = -inorm*dely;
+ eij[2] = -inorm*delz;
local_cut2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
@@ -298,7 +304,7 @@ void PairSpinExchange::compute(int eflag, int vflag)
} 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]);
+ evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
}
}
@@ -317,8 +323,8 @@ void PairSpinExchange::compute_single_pair(int ii, double fmi[3])
double **x = atom->x;
double **sp = atom->sp;
double local_cut2;
-
double xi[3], rij[3];
+ double delx,dely,delz;
double spj[3];
int i,j,jnum,itype,jtype;
@@ -351,15 +357,14 @@ void PairSpinExchange::compute_single_pair(int ii, double fmi[3])
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];
+ delx = xi[0] - x[j][0];
+ dely = xi[1] - x[j][1];
+ delz = xi[2] - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
if (rsq <= local_cut2) {
compute_exchange(i,j,rsq,fmi,spj);
}
-
}
}
@@ -390,7 +395,8 @@ void PairSpinExchange::compute_exchange(int i, int j, double rsq, double fmi[3],
compute the mechanical force due to the exchange interaction between atom i and atom j
------------------------------------------------------------------------- */
-void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double rij[3], double fi[3], double spi[3], double spj[3])
+void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double eij[3],
+ double fi[3], double spi[3], double spj[3])
{
int *type = atom->type;
int itype, jtype;
@@ -408,9 +414,9 @@ void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double ri
Jex_mech *= 8.0*Jex*rr*exp(-ra);
Jex_mech *= (spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]);
- fi[0] -= Jex_mech*rij[0];
- fi[1] -= Jex_mech*rij[1];
- fi[2] -= Jex_mech*rij[2];
+ fi[0] -= Jex_mech*eij[0];
+ fi[1] -= Jex_mech*eij[1];
+ fi[2] -= Jex_mech*eij[2];
}
/* ----------------------------------------------------------------------
diff --git a/src/SPIN/pair_spin_magelec.cpp b/src/SPIN/pair_spin_magelec.cpp
index 315b691d17ea38a265a16587ee1191ef7ca0f7cd..6bc1f71947333ac405e29d6342696165b47379ab 100644
--- a/src/SPIN/pair_spin_magelec.cpp
+++ b/src/SPIN/pair_spin_magelec.cpp
@@ -187,8 +187,14 @@ void PairSpinMagelec::init_style()
double PairSpinMagelec::init_one(int i, int j)
{
+ if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
- if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+ ME[j][i] = ME[i][j];
+ ME_mech[j][i] = ME_mech[i][j];
+ v_mex[j][i] = v_mex[i][j];
+ v_mey[j][i] = v_mey[i][j];
+ v_mez[j][i] = v_mez[i][j];
+ cut_spin_magelec[j][i] = cut_spin_magelec[i][j];
return cut_spin_magelec_global;
}
@@ -211,7 +217,8 @@ void PairSpinMagelec::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double evdwl, ecoul;
- double xi[3], rij[3], eij[3];
+ double xi[3], eij[3];
+ double delx,dely,delz;
double spi[3], spj[3];
double fi[3], fmi[3];
double local_cut2;
@@ -263,18 +270,17 @@ void PairSpinMagelec::compute(int eflag, int vflag)
spj[2] = sp[j][2];
evdwl = 0.0;
-
fi[0] = fi[1] = fi[2] = 0.0;
fmi[0] = fmi[1] = fmi[2] = 0.0;
- 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];
+ delx = xi[0] - x[j][0];
+ dely = xi[1] - x[j][1];
+ delz = xi[2] - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
inorm = 1.0/sqrt(rsq);
- eij[0] = inorm*rij[0];
- eij[1] = inorm*rij[1];
- eij[2] = inorm*rij[2];
+ eij[0] = -inorm*delx;
+ eij[1] = -inorm*dely;
+ eij[2] = -inorm*delz;
local_cut2 = cut_spin_magelec[itype][jtype]*cut_spin_magelec[itype][jtype];
@@ -301,12 +307,12 @@ void PairSpinMagelec::compute(int eflag, int vflag)
}
if (eflag) {
- evdwl = (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
+ 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]);
+ evdwl,ecoul,fi[0],fi[1],fi[2],delx,dely,delz);
}
}
@@ -322,8 +328,8 @@ void PairSpinMagelec::compute_single_pair(int ii, double fmi[3])
double **x = atom->x;
double **sp = atom->sp;
double local_cut2;
-
- double xi[3], rij[3], eij[3];
+ double xi[3], eij[3];
+ double delx,dely,delz;
double spj[3];
int i,j,jnum,itype,jtype;
@@ -342,8 +348,6 @@ void PairSpinMagelec::compute_single_pair(int ii, double fmi[3])
xi[1] = x[i][1];
xi[2] = x[i][2];
- eij[0] = eij[1] = eij[2] = 0.0;
-
jlist = firstneigh[i];
jnum = numneigh[i];
@@ -358,14 +362,14 @@ void PairSpinMagelec::compute_single_pair(int ii, double fmi[3])
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];
+ delx = xi[0] - x[j][0];
+ dely = xi[1] - x[j][1];
+ delz = xi[2] - x[j][2];
+ rsq = delx*delx + dely*dely + delz*delz;
inorm = 1.0/sqrt(rsq);
- eij[0] = inorm*rij[0];
- eij[1] = inorm*rij[1];
- eij[2] = inorm*rij[2];
+ eij[0] = -inorm*delx;
+ eij[1] = -inorm*dely;
+ eij[2] = -inorm*delz;
if (rsq <= local_cut2) {
compute_magelec(i,j,rsq,eij,fmi,spj);
@@ -380,36 +384,26 @@ void PairSpinMagelec::compute_magelec(int i, int j, double rsq, double eij[3], d
{
int *type = atom->type;
int itype, jtype;
+ double meix,meiy,meiz;
+ double vx,vy,vz;
itype = type[i];
jtype = type[j];
- double local_cut2 = cut_spin_magelec[itype][jtype]*cut_spin_magelec[itype][jtype];
-
- if (rsq <= local_cut2) {
- double meix,meiy,meiz;
- double rx, ry, rz;
- double vx, vy, vz;
-
- rx = eij[0];
- ry = eij[1];
- rz = eij[2];
-
- vx = v_mex[itype][jtype];
- vy = v_mey[itype][jtype];
- vz = v_mez[itype][jtype];
+ vx = v_mex[itype][jtype];
+ vy = v_mey[itype][jtype];
+ vz = v_mez[itype][jtype];
- meix = vy*rz - vz*ry;
- meiy = vz*rx - vx*rz;
- meiz = vx*ry - vy*rx;
+ meix = vy*eij[2] - vz*eij[1];
+ meiy = vz*eij[0] - vx*eij[2];
+ meiz = vx*eij[1] - vy*eij[0];
- meix *= ME[itype][jtype];
- meiy *= ME[itype][jtype];
- meiz *= ME[itype][jtype];
+ meix *= ME[itype][jtype];
+ meiy *= ME[itype][jtype];
+ meiz *= ME[itype][jtype];
- fmi[0] += spj[1]*meiz - spj[2]*meiy;
- fmi[1] += spj[2]*meix - spj[0]*meiz;
- fmi[2] += spj[0]*meiy - spj[1]*meix;
- }
+ fmi[0] += spj[1]*meiz - spj[2]*meiy;
+ fmi[1] += spj[2]*meix - spj[0]*meiz;
+ fmi[2] += spj[0]*meiy - spj[1]*meix;
}
/* ---------------------------------------------------------------------- */
diff --git a/src/SPIN/pair_spin_neel.cpp b/src/SPIN/pair_spin_neel.cpp
index 0daafad7567d32ff841a2bfd5bf0a5e400d1eeeb..55f537cf4fe1c4b2af7de3dc3a1945cca14b4f0e 100644
--- a/src/SPIN/pair_spin_neel.cpp
+++ b/src/SPIN/pair_spin_neel.cpp
@@ -193,8 +193,16 @@ void PairSpinNeel::init_style()
double PairSpinNeel::init_one(int i, int j)
{
-
- if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+ if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
+
+ g1[j][i] = g1[i][j];
+ g1_mech[j][i] = g1_mech[i][j];
+ g2[j][i] = g2[i][j];
+ g3[j][i] = g3[i][j];
+ q1[j][i] = q1[i][j];
+ q1_mech[j][i] = q1_mech[i][j];
+ q2[j][i] = q2[i][j];
+ q3[j][i] = q3[i][j];
return cut_spin_neel_global;
}
diff --git a/src/SPIN/pair_spin_neel.h b/src/SPIN/pair_spin_neel.h
index 934d4a93ad5f2bdd37a53da2403d9c7d7c049c76..f60d7d2dca0e6ac72ea3c0111698f79d60ea2f37 100644
--- a/src/SPIN/pair_spin_neel.h
+++ b/src/SPIN/pair_spin_neel.h
@@ -51,9 +51,9 @@ class PairSpinNeel : public PairSpin {
// pseudo-dipolar and pseudo-quadrupolar coeff.
- double **g1, **g1_mech; // exchange coeffs gij
+ double **g1, **g1_mech; // neel coeffs gij
double **g2, **g3; // g1 in eV, g2 adim, g3 in Ang
- double **q1, **q1_mech; // exchange coeffs qij
+ double **q1, **q1_mech; // neel coeffs qij
double **q2, **q3; // q1 in eV, q2 adim, q3 in Ang
double **cut_spin_neel; // cutoff distance exchange