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

Commit JT 072618 

- improvements documentation (dmi and exchange)
- correction error cross product in pair_spin_dmi.cpp
- implementation mech. part in pair_spin_dmi.cpp
- correction in all pairs: init_one for [j][i] couples
- correction in atom_vec_spin.cpp: index error in read_data
- some improvements in pair_spin_dmi.cpp and pair_spin_magelec.cpp
This commit is contained in:
julient31 2018-07-26 08:36:30 -06:00 committed by Axel Kohlmeyer
parent ae0979e1ad
commit 28993d9823
20 changed files with 272 additions and 157 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/src/Eqs/pair_spin_dmi_forces.jpg Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 13 KiB

14
doc/src/Eqs/pair_spin_dmi_forces.tex Normal file
View File

@ -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}

BIN
doc/src/Eqs/pair_spin_dmi_interaction.jpg
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 6.2 KiB

After

Width:  |  Height:  |  Size: 7.7 KiB

2
doc/src/Eqs/pair_spin_dmi_interaction.tex
View File

@ -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

BIN
doc/src/Eqs/pair_spin_exchange_forces.jpg
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 13 KiB

After

Width:  |  Height:  |  Size: 13 KiB

10
doc/src/Eqs/pair_spin_exchange_forces.tex
View File

@ -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}

BIN
doc/src/Eqs/pair_spin_exchange_interaction.jpg
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 6.5 KiB

After

Width:  |  Height:  |  Size: 5.8 KiB

2
doc/src/Eqs/pair_spin_exchange_interaction.tex
View File

@ -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}

41
doc/src/pair_spin_dmi.txt
View File

@ -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).

28
doc/src/pair_spin_exchange.txt
View File

@ -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:]

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

@ -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

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

@ -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

6
src/SPIN/atom_vec_spin.cpp
View File

@ -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]);

20
src/SPIN/compute_spin.cpp
View File

@ -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");

116
src/SPIN/pair_spin_dmi.cpp
View File

@ -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);
}
}

15
src/SPIN/pair_spin_dmi.h
View File

@ -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();
};

52
src/SPIN/pair_spin_exchange.cpp
View File

@ -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];
}
/* ----------------------------------------------------------------------

88
src/SPIN/pair_spin_magelec.cpp
View File

@ -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];
vx = v_mex[itype][jtype];
vy = v_mey[itype][jtype];
vz = v_mez[itype][jtype];
if (rsq <= local_cut2) {
double meix,meiy,meiz;
double rx, ry, rz;
double vx, vy, vz;
meix = vy*eij[2] - vz*eij[1];
meiy = vz*eij[0] - vx*eij[2];
meiz = vx*eij[1] - vy*eij[0];
rx = eij[0];
ry = eij[1];
rz = eij[2];
meix *= ME[itype][jtype];
meiy *= ME[itype][jtype];
meiz *= ME[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 *= 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;
}
/* ---------------------------------------------------------------------- */

10
src/SPIN/pair_spin_neel.cpp
View File

@ -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;
}

4
src/SPIN/pair_spin_neel.h
View File

@ -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