ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Commit JT 051719

Browse Source 
- removed qsymp pair style
- cleaned doc (pair/spin/diole and kspace_style)
- cleaned kspace .cpp/h files
This commit is contained in:
julient31
2019-05-17 15:04:14 -06:00
parent 41872e37e6
commit fbb78e7b78
16 changed files with 294 additions and 1073 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 38 KiB

42
doc/src/Eqs/pair_spin_dipole.tex Normal file
View File

@ -0,0 +1,42 @@
\documentclass[preview]{standalone}
\usepackage{varwidth}
\usepackage[utf8x]{inputenc}
\usepackage{amsmath,amssymb,graphics,bm,setspace}
\begin{document}
\begin{varwidth}{50in}
\begin{equation}
\mathcal{H}_{\rm long}=
-\frac{\mu_{0} \left( \mu_B\right)^2}{4\pi}
\sum_{i,j,i\neq j}^{N}
\frac{g_i g_j}{r_{ij}^3}
\Big(3
\left(\bm{e}_{ij}\cdot \bm{s}_{i}\right)
\left(\bm{e}_{ij}\cdot \bm{s}_{j}\right)
-\bm{s}_i\cdot\bm{s}_j \Big)
\nonumber
\end{equation}
\begin{equation}
\bm{\omega}_i =
\frac{\mu_0 (\mu_B)^2}{4\pi\hbar}\sum_{j}
\frac{g_i g_j}{r_{ij}^3}
\, \Big(
3\,(\bm{e}_{ij}\cdot\bm{s}_{j})\bm{e}_{ij}
-\bm{s}_{j} \Big) \nonumber
\end{equation}
\begin{equation}
\bm{F}_i =
\frac{3\, \mu_0 (\mu_B)^2}{4\pi} \sum_j
\frac{g_i g_j}{r_{ij}^4}
\Big[\big( (\bm{s}_i\cdot\bm{s}_j)
-5(\bm{e}_{ij}\cdot\bm{s}_i)
(\bm{e}_{ij}\cdot\bm{s}_j)\big) \bm{e}_{ij}+
\big(
(\bm{e}_{ij}\cdot\bm{s}_i)\bm{s}_j+
(\bm{e}_{ij}\cdot\bm{s}_j)\bm{s}_i
\big)
\Big]
\nonumber
\end{equation}
\end{varwidth}
\end{document}

24
doc/src/kspace_style.txt
View File

@ -20,6 +20,10 @@ style = {none} or {ewald} or {ewald/disp} or {ewald/omp} or {pppm} or {pppm/cg}
accuracy = desired relative error in forces
{ewald/omp} value = accuracy
accuracy = desired relative error in forces
{ewald/dipole} value = accuracy
accuracy = desired relative error in forces
{ewald/dipole/spin} value = accuracy
accuracy = desired relative error in forces
{pppm} value = accuracy
accuracy = desired relative error in forces
{pppm/cg} values = accuracy (smallq)
@ -47,6 +51,10 @@ style = {none} or {ewald} or {ewald/disp} or {ewald/omp} or {pppm} or {pppm/cg}
accuracy = desired relative error in forces
{pppm/stagger} value = accuracy
accuracy = desired relative error in forces
{pppm/dipole} value = accuracy
accuracy = desired relative error in forces
{pppm/dipole/spin} value = accuracy
accuracy = desired relative error in forces
{msm} value = accuracy
accuracy = desired relative error in forces
{msm/cg} value = accuracy (smallq)
@ -105,9 +113,13 @@ The {ewald/disp} style adds a long-range dispersion sum option for
but in a more efficient manner than the {ewald} style. The 1/r^6
capability means that Lennard-Jones or Buckingham potentials can be
used without a cutoff, i.e. they become full long-range potentials.
The {ewald/disp} style can also be used with point-dipoles
"(Toukmaji)"_#Toukmaji and is currently the only kspace solver in
LAMMPS with this capability.
The {ewald/dipole} style adds long-range standard Ewald summations
for dipole-dipole interactions.
The {ewald/dipole/spin} style adds long-range standard Ewald
summations for magnetic dipole-dipole interactions between
magnetic spins.
:line
@ -128,6 +140,12 @@ The optional {smallq} argument defines the cutoff for the absolute
charge value which determines whether a particle is considered charged
or not. Its default value is 1.0e-5.
The {pppm/dipole} style invokes a particle-particle particle-mesh solver
for dipole-dipole interactions.
The {pppm/dipole/spin} style invokes a particle-particle particle-mesh solver
for magnetic dipole-dipole interactions between magnetic spins.
The {pppm/tip4p} style is identical to the {pppm} style except that it
adds a charge at the massless 4th site in each TIP4P water molecule.
It should be used with "pair styles"_pair_style.html with a

214
doc/src/pair_spin_dipole.txt
View File

@ -8,15 +8,13 @@
pair_style spin/dipole/cut command :h3
pair_style spin/dipole/long command :h3
pair_style spin/dipole/long/qsymp command :h3
[Syntax:]
pair_style spin/dipole/cut cutoff
pair_style spin/dipole/long cutoff
pair_style spin/dipole/long/qsymp cutoff :pre
pair_style spin/dipole/long cutoff :pre
cutoff = global cutoff for Magnetic dipole energy and forces
cutoff = global cutoff for magnetic dipole energy and forces
(optional) (distance units) :ulb,l
:ule
@ -31,179 +29,31 @@ pair_coeff * * 1.0 1.0
pair_coeff 2 3 1.0 1.0 2.5 4.0 scale 0.5
pair_coeff 2 3 1.0 1.0 2.5 4.0 :pre
pair_style spin/dipole/long/qsymp 10.0
pair_coeff * * 1.0 1.0
pair_coeff 2 3 1.0 1.0 2.5 4.0 :pre
[Description:]
Style {spin/dipole/cut} computes a short-range dipole-dipole
interactions between pairs of magnetic particles that each
interaction between pairs of magnetic particles that each
have a magnetic spin.
The magnetic dipole-dipole interactions are computed by the
following formulas for the energy (E), force
(F), and torque (T) between particles I and J.
following formulas for the magnetic energy, magnetic precession
vector omega and mechanical force between particles I and J.
:c,image(Eqs/pair_dipole.jpg)
:c,image(Eqs/pair_spin_dipole.jpg)
where qi and qj are the charges on the two particles, pi and pj are
the dipole moment vectors of the two particles, r is their separation
distance, and the vector r = Ri - Rj is the separation vector between
the two particles. Note that Eqq and Fqq are simply Coulombic energy
and force, Fij = -Fji as symmetric forces, and Tij != -Tji since the
torques do not act symmetrically. These formulas are discussed in
"(Allen)"_#Allen2 and in "(Toukmaji)"_#Toukmaji2.
where si and sj are the spin on two magnetic particles,
r is their separation distance, and the vector e = (Ri - Rj)/|Ri - Rj|
is the direction vector between the two particles.
Also note, that in the code, all of these terms (except Elj) have a
C/epsilon prefactor, the same as the Coulombic term in the LJ +
Coulombic pair styles discussed "here"_pair_lj.html. C is an
energy-conversion constant and epsilon is the dielectric constant
which can be set by the "dielectric"_dielectric.html command. The
same is true of the equations that follow for other dipole pair
styles.
Style {lj/sf/dipole/sf} computes "shifted-force" interactions between
pairs of particles that each have a charge and/or a point dipole
moment. In general, a shifted-force potential is a (slightly) modified
potential containing extra terms that make both the energy and its
derivative go to zero at the cutoff distance; this removes
(cutoff-related) problems in energy conservation and any numerical
instability in the equations of motion "(Allen)"_#Allen2. Shifted-force
interactions for the Lennard-Jones (E_LJ), charge-charge (Eqq),
charge-dipole (Eqp), dipole-charge (Epq) and dipole-dipole (Epp)
potentials are computed by these formulas for the energy (E), force
(F), and torque (T) between particles I and J:
:c,image(Eqs/pair_dipole_sf.jpg)
:c,image(Eqs/pair_dipole_sf2.jpg)
where epsilon and sigma are the standard LJ parameters, r_c is the
cutoff, qi and qj are the charges on the two particles, pi and pj are
the dipole moment vectors of the two particles, r is their separation
distance, and the vector r = Ri - Rj is the separation vector between
the two particles. Note that Eqq and Fqq are simply Coulombic energy
and force, Fij = -Fji as symmetric forces, and Tij != -Tji since the
torques do not act symmetrically. The shifted-force formula for the
Lennard-Jones potential is reported in "(Stoddard)"_#Stoddard. The
original (non-shifted) formulas for the electrostatic potentials,
forces and torques can be found in "(Price)"_#Price2. The shifted-force
electrostatic potentials have been obtained by applying equation 5.13
of "(Allen)"_#Allen2. The formulas for the corresponding forces and
torques have been obtained by applying the 'chain rule' as in appendix
C.3 of "(Allen)"_#Allen2.
If one cutoff is specified in the pair_style command, it is used for
both the LJ and Coulombic (q,p) terms. If two cutoffs are specified,
they are used as cutoffs for the LJ and Coulombic (q,p) terms
respectively. This pair style also supports an optional {scale} keyword
as part of a pair_coeff statement, where the interactions can be
scaled according to this factor. This scale factor is also made available
for use with fix adapt.
Style {lj/cut/dipole/long} computes long-range point-dipole
interactions as discussed in "(Toukmaji)"_#Toukmaji2. Dipole-dipole,
dipole-charge, and charge-charge interactions are all supported, along
with the standard 12/6 Lennard-Jones interactions, which are computed
with a cutoff. A "kspace_style"_kspace_style.html must be defined to
use this pair style. Currently, only "kspace_style
ewald/disp"_kspace_style.html support long-range point-dipole
interactions.
Style {lj/long/dipole/long} also computes point-dipole interactions as
discussed in "(Toukmaji)"_#Toukmaji2. Long-range dipole-dipole,
dipole-charge, and charge-charge interactions are all supported, along
with the standard 12/6 Lennard-Jones interactions. LJ interactions
can be cutoff or long-ranged.
For style {lj/long/dipole/long}, if {flag_lj} is set to {long}, no
cutoff is used on the LJ 1/r^6 dispersion term. The long-range
portion is calculated by using the "kspace_style
ewald_disp"_kspace_style.html command. The specified LJ cutoff then
determines which portion of the LJ interactions are computed directly
by the pair potential versus which part is computed in reciprocal
space via the Kspace style. If {flag_lj} is set to {cut}, the LJ
interactions are simply cutoff, as with "pair_style
lj/cut"_pair_lj.html. If {flag_lj} is set to {off}, LJ interactions
are not computed at all.
If {flag_coul} is set to {long}, no cutoff is used on the Coulombic or
dipole interactions. The long-range portion is calculated by using
{ewald_disp} of the "kspace_style"_kspace_style.html command. If
{flag_coul} is set to {off}, Coulombic and dipole interactions are not
computed at all.
Atoms with dipole moments should be integrated using the "fix
nve/sphere update dipole"_fix_nve_sphere.html or the "fix
nvt/sphere update dipole"_fix_nvt_sphere.html command to rotate the
dipole moments. The {omega} option on the "fix
langevin"_fix_langevin.html command can be used to thermostat the
rotational motion. The "compute temp/sphere"_compute_temp_sphere.html
command can be used to monitor the temperature, since it includes
rotational degrees of freedom. The "atom_style
hybrid dipole sphere"_atom_style.html command should be used since
it defines the point dipoles and their rotational state.
The magnitude and orientation of the dipole moment for each particle
can be defined by the "set"_set.html command or in the "Atoms" section
of the data file read in by the "read_data"_read_data.html command.
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, or by mixing as described below:
epsilon (energy units)
sigma (distance units)
cutoff1 (distance units)
cutoff2 (distance units) :ul
The latter 2 coefficients are optional. If not specified, the global
LJ and Coulombic cutoffs specified in the pair_style command are used.
If only one cutoff is specified, it is used as the cutoff for both LJ
and Coulombic interactions for this type pair. If both coefficients
are specified, they are used as the LJ and Coulombic cutoffs for this
type pair.
Style {spin/dipole/long} computes long-range magnetic dipole-dipole
interaction.
A "kspace_style"_kspace_style.html must be defined to
use this pair style. Currently, "kspace_style
ewald/dipole/spin"_kspace_style.html and "kspace_style
pppm/dipole/spin"_kspace_style.html support long-range magnetic
dipole-dipole interactions.
:line
Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
functionally the same as the corresponding style without the suffix.
They have been optimized to run faster, depending on your available
hardware, as discussed on the "Speed packages"_Speed_packages.html doc
page. The accelerated styles take the same arguments and should
produce the same results, except for round-off and precision issues.
These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
USER-OMP and OPT packages, respectively. They are only enabled if
LAMMPS was built with those packages. See the "Build
package"_Build_package.html doc page for more info.
You can specify the accelerated styles explicitly in your input script
by including their suffix, or you can use the "-suffix command-line
switch"_Run_options.html when you invoke LAMMPS, or you can use the
"suffix"_suffix.html command in your input script.
See the "Speed packages"_Speed_packages.html doc page for more
instructions on how to use the accelerated styles effectively.
:line
[Mixing, shift, table, tail correction, restart, rRESPA info]:
For atom type pairs I,J and I != J, the epsilon and sigma coefficients
and cutoff distances for this pair style can be mixed. The default
mix value is {geometric}. See the "pair_modify" command for details.
For atom type pairs I,J and I != J, the A, sigma, d1, and d2
coefficients and cutoff distance for this pair style can be mixed. A
is an energy value mixed like a LJ epsilon. D1 and d2 are distance
values and are mixed like sigma. The default mix value is
{geometric}. See the "pair_modify" command for details.
This pair style does not support the "pair_modify"_pair_modify.html
shift option for the energy of the Lennard-Jones portion of the pair
interaction; such energy goes to zero at the cutoff by construction.
The "pair_modify"_pair_modify.html table option is not relevant
for this pair style.
@ -215,28 +65,20 @@ This pair style writes its information to "binary restart
files"_restart.html, so pair_style and pair_coeff commands do not need
to be specified in an input script that reads a restart file.
This pair style can only be used via the {pair} keyword of the
"run_style respa"_run_style.html command. It does not support the
{inner}, {middle}, {outer} keywords.
[Restrictions:]
The {lj/cut/dipole/cut}, {lj/cut/dipole/long}, and
{lj/long/dipole/long} styles are part of the DIPOLE package. They are
only enabled if LAMMPS was built with that package. See the "Build
package"_Build_package.html doc page for more info.
The {spin/dipole/cut} and {spin/dipole/long} styles are part of
the SPIN package. They are only enabled if LAMMPS was built with that
package. See the "Build package"_Build_package.html doc page for more
info.
The {lj/sf/dipole/sf} style is part of the USER-MISC package. It is
only enabled if LAMMPS was built with that package. See the "Build
package"_Build_package.html doc page for more info.
Using dipole pair styles with {electron} "units"_units.html is not
Using dipole/spin pair styles with {electron} "units"_units.html is not
currently supported.
[Related commands:]
"pair_coeff"_pair_coeff.html, "set"_set.html, "read_data"_read_data.html,
"fix nve/sphere"_fix_nve_sphere.html, "fix nvt/sphere"_fix_nvt_sphere.html
"pair_coeff"_pair_coeff.html, "kspace_style"_kspace_style.html
"fix nve/spin"_fix_nve_spin.html
[Default:] none
@ -245,13 +87,3 @@ currently supported.
:link(Allen2)
[(Allen)] Allen and Tildesley, Computer Simulation of Liquids,
Clarendon Press, Oxford, 1987.
:link(Toukmaji2)
[(Toukmaji)] Toukmaji, Sagui, Board, and Darden, J Chem Phys, 113,
10913 (2000).
:link(Stoddard)
[(Stoddard)] Stoddard and Ford, Phys Rev A, 8, 1504 (1973).
:link(Price2)
[(Price)] Price, Stone and Alderton, Mol Phys, 52, 987 (1984).

9
examples/SPIN/dipole_spin/in.spin.iron_ewald
View File

@ -21,12 +21,16 @@ mass 1 55.845
set group all spin 2.2 -1.0 0.0 0.0
velocity all create 100 4928459 rot yes dist gaussian
pair_style hybrid/overlay eam/alloy spin/exchange 3.5
pair_style hybrid/overlay eam/alloy spin/exchange 3.5 spin/dipole/long 8.0
pair_coeff * * eam/alloy Fe_Mishin2006.eam.alloy Fe
pair_coeff * * spin/exchange exchange 3.4 0.02726 0.2171 1.841
pair_coeff * * spin/dipole/long 8.0
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
kspace_style ewald/dipole/spin 1.0e-4
fix 1 all precession/spin cubic 0.001 0.0005 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0
fix_modify 1 energy yes
fix 2 all langevin/spin 0.0 0.0 21
@ -55,6 +59,3 @@ compute outsp all property/atom spx spy spz sp fmx fmy fmz
dump 100 all custom 1 dump_iron.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3]
run 2000
# min_style spin
# min_modify alpha_damp 1.0 discrete_factor 10
# minimize 1.0e-16 1.0e-16 10000 10000

7
examples/SPIN/dipole_spin/in.spin.iron_pppm
View File

@ -21,12 +21,17 @@ mass 1 55.845
set group all spin 2.2 -1.0 0.0 0.0
velocity all create 100 4928459 rot yes dist gaussian
pair_style hybrid/overlay eam/alloy spin/exchange 3.5
pair_style hybrid/overlay eam/alloy spin/exchange 3.5 spin/dipole/long 8.0
pair_coeff * * eam/alloy Fe_Mishin2006.eam.alloy Fe
pair_coeff * * spin/exchange exchange 3.4 0.02726 0.2171 1.841
pair_coeff * * spin/dipole/long 8.0
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
kspace_style pppm/dipole/spin 1.0e-4
kspace_modify compute yes
fix 1 all precession/spin cubic 0.001 0.0005 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0
fix_modify 1 energy yes
fix 2 all langevin/spin 0.0 0.0 21

12
src/.gitignore vendored
View File

@ -166,6 +166,10 @@
/pair_spin.h
/pair_spin_dmi.cpp
/pair_spin_dmi.h
/pair_spin_dipole_cut.cpp
/pair_spin_dipole_cut.h
/pair_spin_dipole_long.cpp
/pair_spin_dipole_long.h
/pair_spin_exchange.cpp
/pair_spin_exchange.h
/pair_spin_magelec.cpp
@ -427,6 +431,10 @@
/ewald.h
/ewald_cg.cpp
/ewald_cg.h
/ewald_dipole.cpp
/ewald_dipole.h
/ewald_dipole_spin.cpp
/ewald_dipole_spin.h
/ewald_disp.cpp
/ewald_disp.h
/ewald_n.cpp
@ -1026,6 +1034,10 @@
/pppm.h
/pppm_cg.cpp
/pppm_cg.h
/pppm_dipole.cpp
/pppm_dipole.h
/pppm_dipole_spin.cpp
/pppm_dipole_spin.h
/pppm_disp.cpp
/pppm_disp.h
/pppm_disp_tip4p.cpp

1
src/KSPACE/ewald_dipole.cpp
View File

@ -94,7 +94,6 @@ void EwaldDipole::init()
error->all(FLERR,"Cannot use EwaldDipole with 2d simulation");
if (!atom->mu) error->all(FLERR,"Kspace style requires atom attribute mu");
//if (!atom->q_flag) error->all(FLERR,"Kspace style requires atom attribute q");
if (dipoleflag && strcmp(update->unit_style,"electron") == 0)
error->all(FLERR,"Cannot (yet) use 'electron' units with dipoles");

5
src/KSPACE/ewald_dipole_spin.cpp
View File

@ -51,8 +51,6 @@ EwaldDipoleSpin::EwaldDipoleSpin(LAMMPS *lmp) :
spinflag = 1;
hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
//mub = 5.78901e-5; // in eV/T
//mu_0 = 1.2566370614e-6; // in T.m/A
mub = 9.274e-4; // in A.Ang^2
mu_0 = 785.15; // in eV/Ang/A^2
mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV.Ang^3
@ -437,7 +435,6 @@ void EwaldDipoleSpin::compute(int eflag, int vflag)
// compute field for torque calculation
//partial2 = exprl*sfacrl_all[k] + expim*sfacim_all[k];
partial_peratom = exprl*sfacrl_all[k] + expim*sfacim_all[k];
tk[i][0] += partial2*eg[k][0];
tk[i][1] += partial2*eg[k][1];
@ -456,12 +453,10 @@ void EwaldDipoleSpin::compute(int eflag, int vflag)
// (for per-atom energy and virial calc.)
if (evflag_atom) {
//partial_peratom = exprl*sfacrl_all[k] + expim*sfacim_all[k];
if (eflag_atom) eatom[i] += mudotk*ug[k]*partial_peratom;
if (vflag_atom)
for (j = 0; j < 6; j++)
vatom[i][j] += (ug[k]*mudotk*vg[k][j]*partial_peratom - vcik[j]);
//vatom[i][j] += ug[k] * (vg[k][j]*partial_peratom - vcik[j]);
}
}
}

4
src/KSPACE/pppm_dipole_spin.cpp
View File

@ -70,8 +70,6 @@ PPPMDipoleSpin::PPPMDipoleSpin(LAMMPS *lmp) :
spinflag = 1;
hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
//mub = 5.78901e-5; // in eV/T
//mu_0 = 1.2566370614e-6; // in T.m/A
mub = 9.274e-4; // in A.Ang^2
mu_0 = 785.15; // in eV/Ang/A^2
mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV.Ang^3
@ -157,7 +155,7 @@ void PPPMDipoleSpin::init()
int itmp = 0;
double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
// probably not the correct extract here
// check the correct extract here
if (p_cutoff == NULL)
error->all(FLERR,"KSpace style is incompatible with Pair style");
cutoff = *p_cutoff;

141
src/SPIN/pair_spin_dipole_cut.cpp
View File

@ -58,7 +58,7 @@ lockfixnvespin(NULL)
mub = 9.274e-4; // in A.Ang^2
mu_0 = 785.15; // in eV/Ang/A^2
mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV.Ang^3
mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV
//mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV
mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz
}
@ -115,8 +115,8 @@ void PairSpinDipoleCut::settings(int narg, char **arg)
void PairSpinDipoleCut::coeff(int narg, char **arg)
{
if (!allocated) allocate();
if (narg < 1 || narg > 3)
if (narg != 3)
error->all(FLERR,"Incorrect args in pair_style command");
int ilo,ihi,jlo,jhi;
@ -216,24 +216,21 @@ void *PairSpinDipoleCut::extract(const char *str, int &dim)
void PairSpinDipoleCut::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double rinv,r2inv,r3inv,rsq;
double evdwl,ecoul;
double xi[3],rij[3],eij[3];
double spi[4],spj[4],fi[3],fmi[3];
double local_cut2;
int *ilist,*jlist,*numneigh,**firstneigh;
double rinv,r2inv,r3inv,rsq,local_cut2,evdwl,ecoul;
double xi[3],rij[3],eij[3],spi[4],spj[4],fi[3],fmi[3];
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double **x = atom->x;
double **f = atom->f;
double **fm = atom->fm;
double **sp = atom->sp;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
@ -320,70 +317,90 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
/* ----------------------------------------------------------------------
update the pair interaction fmi acting on the spin ii
adding 1/r (for r in [0,rc]) contribution to the pair
removing erf(r)/r (for r in [0,rc]) from the kspace force
------------------------------------------------------------------------- */
void PairSpinDipoleCut::compute_single_pair(int ii, double fmi[3])
{
int i,j,jj,jnum,itype,jtype;
double rsq,rinv,r2inv,r3inv;
double xi[3],rij[3],eij[3];
double spi[4],spj[4];
double local_cut2;
int j,jnum,itype,jtype,ntypes;
int *ilist,*jlist,*numneigh,**firstneigh;
double rsq,rinv,r2inv,r3inv,local_cut2;
double xi[3],rij[3],eij[3],spi[4],spj[4];
int k,locflag;
int *type = atom->type;
double **x = atom->x;
double **sp = atom->sp;
int *type = atom->type;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// computation of the exchange interaction
// loop over neighbors of atom i
i = ilist[ii];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
spi[3] = sp[i][3];
jlist = firstneigh[i];
jnum = numneigh[i];
itype = type[i];
// check if interaction applies to type of ii
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
itype = type[ii];
ntypes = atom->ntypes;
locflag = 0;
k = 1;
while (k <= ntypes) {
if (k <= itype) {
if (setflag[k][itype] == 1) {
locflag =1;
break;
}
k++;
} else if (k > itype) {
if (setflag[itype][k] == 1) {
locflag =1;
break;
}
k++;
} else error->all(FLERR,"Wrong type number");
}
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
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];
rinv = 1.0/sqrt(rsq);
eij[0] = rij[0]*rinv;
eij[1] = rij[1]*rinv;
eij[2] = rij[2]*rinv;
// if interaction applies to type ii,
// locflag = 1 and compute pair interaction
if (locflag == 1) {
local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype];
xi[0] = x[ii][0];
xi[1] = x[ii][1];
xi[2] = x[ii][2];
spi[0] = sp[ii][0];
spi[1] = sp[ii][1];
spi[2] = sp[ii][2];
spi[3] = sp[ii][3];
jlist = firstneigh[ii];
jnum = numneigh[ii];
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
if (rsq < local_cut2) {
r2inv = 1.0/rsq;
r3inv = r2inv*rinv;
// compute dipolar interaction
compute_dipolar(i,j,eij,fmi,spi,spj,r3inv);
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
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];
rinv = 1.0/sqrt(rsq);
eij[0] = rij[0]*rinv;
eij[1] = rij[1]*rinv;
eij[2] = rij[2]*rinv;
local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype];
if (rsq < local_cut2) {
r2inv = 1.0/rsq;
r3inv = r2inv*rinv;
// compute dipolar interaction
compute_dipolar(ii,j,eij,fmi,spi,spj,r3inv);
}
}
}
}
@ -424,11 +441,11 @@ void PairSpinDipoleCut::compute_dipolar_mech(int i, int j, double eij[3],
sjeij = spj[0]*eij[0] + spj[1]*eij[1] + spj[2]*eij[2];
bij = sisj - 5.0*sieij*sjeij;
pre = mub2mu0*gigjri4;
pre = 3.0*mub2mu0*gigjri4;
fi[0] += pre * (eij[0] * bij + (sjeij*spi[0] + sieij*spj[0]));
fi[1] += pre * (eij[1] * bij + (sjeij*spi[1] + sieij*spj[1]));
fi[2] += pre * (eij[2] * bij + (sjeij*spi[2] + sieij*spj[2]));
fi[0] -= pre * (eij[0] * bij + (sjeij*spi[0] + sieij*spj[0]));
fi[1] -= pre * (eij[1] * bij + (sjeij*spi[1] + sieij*spj[1]));
fi[2] -= pre * (eij[2] * bij + (sjeij*spi[2] + sieij*spj[2]));
}
/* ----------------------------------------------------------------------

198
src/SPIN/pair_spin_dipole_long.cpp
View File

@ -60,8 +60,6 @@ lockfixnvespin(NULL)
lattice_flag = 0;
hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
//mub = 5.78901e-5; // in eV/T
//mu_0 = 1.2566370614e-6; // in T.m/A
mub = 9.274e-4; // in A.Ang^2
mu_0 = 785.15; // in eV/Ang/A^2
mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV.Ang^3
@ -120,18 +118,14 @@ void PairSpinDipoleLong::coeff(int narg, char **arg)
{
if (!allocated) allocate();
// check if args correct
if (strcmp(arg[2],"long") != 0)
error->all(FLERR,"Incorrect args in pair_style command");
if (narg < 1 || narg > 4)
if (narg != 3)
error->all(FLERR,"Incorrect args in pair_style command");
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
double spin_long_cut_one = force->numeric(FLERR,arg[3]);
double spin_long_cut_one = force->numeric(FLERR,arg[2]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
@ -181,15 +175,10 @@ void PairSpinDipoleLong::init_style()
// insure use of KSpace long-range solver, set g_ewald
//if (force->kspace == NULL)
// error->all(FLERR,"Pair style requires a KSpace style");
//g_ewald = force->kspace->g_ewald;
// test case
g_ewald = 0.1;
if (force->kspace == NULL)
error->all(FLERR,"Pair style requires a KSpace style");
g_ewald = force->kspace->g_ewald;
}
/* ----------------------------------------------------------------------
@ -312,7 +301,6 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
if (rsq < local_cut2) {
r2inv = 1.0/rsq;
//rinv = sqrt(r2inv);
r = sqrt(rsq);
grij = g_ewald * r;
@ -327,8 +315,6 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
compute_long(i,j,eij,bij,fmi,spi,spj);
compute_long_mech(i,j,eij,bij,fmi,spi,spj);
//compute_long(i,j,rij,bij,fmi,spi,spj);
//compute_long_mech(i,j,rij,bij,fmi,spi,spj);
}
// force accumulation
@ -368,95 +354,117 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
void PairSpinDipoleLong::compute_single_pair(int ii, double fmi[3])
{
int i,j,jj,jnum,itype,jtype;
double r,rinv,r2inv,rsq;
double grij,expm2,t,erfc;
double bij[4];
double xi[3],rij[3],eij[3];
double spi[4],spj[4];
double local_cut2;
double pre1,pre2,pre3;
//int i,j,jj,jnum,itype,jtype;
int j,jj,jnum,itype,jtype,ntypes;
int k,locflag;
int *ilist,*jlist,*numneigh,**firstneigh;
double r,rinv,r2inv,rsq,grij,expm2,t,erfc;
double local_cut2,pre1,pre2,pre3;
double bij[4],xi[3],rij[3],eij[3],spi[4],spj[4];
int *type = atom->type;
double **x = atom->x;
double **sp = atom->sp;
double **fm_long = atom->fm_long;
int *type = atom->type;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// check if interaction applies to type of ii
pre1 = 2.0 * g_ewald / MY_PIS;
pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS;
pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS;
// computation of the exchange interaction
// loop over neighbors of atom i
i = ilist[ii];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
spi[3] = sp[i][3];
jlist = firstneigh[i];
jnum = numneigh[i];
itype = type[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
fmi[0] = fmi[1] = fmi[2] = 0.0;
bij[0] = bij[1] = bij[2] = bij[3] = 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];
rinv = 1.0/sqrt(rsq);
eij[0] = rij[0]*rinv;
eij[1] = rij[1]*rinv;
eij[2] = rij[2]*rinv;
local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype];
if (rsq < local_cut2) {
r2inv = 1.0/rsq;
//rinv = sqrt(r2inv);
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
bij[0] = erfc * rinv;
bij[1] = (bij[0] + pre1*expm2) * r2inv;
bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv;
bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
compute_long(i,j,eij,bij,fmi,spi,spj);
}
itype = type[ii];
ntypes = atom->ntypes;
locflag = 0;
k = 1;
while (k <= ntypes) {
if (k <= itype) {
if (setflag[k][itype] == 1) {
locflag =1;
break;
}
k++;
} else if (k > itype) {
if (setflag[itype][k] == 1) {
locflag =1;
break;
}
k++;
} else error->all(FLERR,"Wrong type number");
}
// adding the kspace components to fm
// if interaction applies to type ii,
// locflag = 1 and compute pair interaction
//printf("test fm before: %g, %g, %g \n",fmi[0],fmi[1],fmi[2]);
//printf("test fm_long: %g, %g, %g \n",fm_long[i][0],fm_long[i][1],fm_long[i][2]);
fmi[0] += fm_long[i][0];
fmi[1] += fm_long[i][1];
fmi[2] += fm_long[i][2];
//printf("test fm after: %g, %g, %g \n",fmi[0],fmi[1],fmi[2]);
if (locflag == 1) {
pre1 = 2.0 * g_ewald / MY_PIS;
pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS;
pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS;
// computation of the exchange interaction
// loop over neighbors of atom i
//i = ilist[ii];
xi[0] = x[ii][0];
xi[1] = x[ii][1];
xi[2] = x[ii][2];
spi[0] = sp[ii][0];
spi[1] = sp[ii][1];
spi[2] = sp[ii][2];
spi[3] = sp[ii][3];
jlist = firstneigh[ii];
jnum = numneigh[ii];
//itype = type[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
fmi[0] = fmi[1] = fmi[2] = 0.0;
bij[0] = bij[1] = bij[2] = bij[3] = 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];
rinv = 1.0/sqrt(rsq);
eij[0] = rij[0]*rinv;
eij[1] = rij[1]*rinv;
eij[2] = rij[2]*rinv;
local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype];
if (rsq < local_cut2) {
r2inv = 1.0/rsq;
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
bij[0] = erfc * rinv;
bij[1] = (bij[0] + pre1*expm2) * r2inv;
bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv;
bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
compute_long(ii,j,eij,bij,fmi,spi,spj);
}
}
// adding the kspace components to fm
fmi[0] += fm_long[ii][0];
fmi[1] += fm_long[ii][1];
fmi[2] += fm_long[ii][2];
}
}
/* ----------------------------------------------------------------------

2
src/SPIN/pair_spin_dipole_long.h
View File

@ -89,7 +89,7 @@ E: Pair dipole/long requires atom attributes q, mu, torque
The atom style defined does not have these attributes.
E: Cannot (yet) use 'electron' units with dipoles
E: Can only use 'metal' units with spins
This feature is not yet supported.

606
src/SPIN/pair_spin_dipole_long_qsymp.cpp
View File

@ -1,606 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
www.cs.sandia.gov/~sjplimp/lammps.html
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------
Contributing authors: Julien Tranchida (SNL)
Stan Moore (SNL)
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "pair_spin_dipole_long_qsymp.h"
#include "atom.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "fix_nve_spin.h"
#include "force.h"
#include "kspace.h"
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "error.h"
#include "update.h"
using namespace LAMMPS_NS;
using namespace MathConst;
#define EWALD_F 1.12837917
#define EWALD_P 0.3275911
#define A1 0.254829592
#define A2 -0.284496736
#define A3 1.421413741
#define A4 -1.453152027
#define A5 1.061405429
/* ---------------------------------------------------------------------- */
PairSpinDipoleLongQsymp::PairSpinDipoleLongQsymp(LAMMPS *lmp) : PairSpin(lmp),
lockfixnvespin(NULL)
{
single_enable = 0;
ewaldflag = pppmflag = spinflag = 1;
respa_enable = 0;
no_virial_fdotr_compute = 1;
lattice_flag = 0;
hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
//mub = 5.78901e-5; // in eV/T
//mu_0 = 1.2566370614e-6; // in T.m/A
mub = 9.274e-4; // in A.Ang^2
mu_0 = 785.15; // in eV/Ang/A^2
mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV.Ang^3
//mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV
mub2mu0 = mub * mub * mu_0 / (4.0*MY_PI); // in eV
mub2mu0hbinv = mub2mu0 / hbar; // in rad.THz
}
/* ----------------------------------------------------------------------
free all arrays
------------------------------------------------------------------------- */
PairSpinDipoleLongQsymp::~PairSpinDipoleLongQsymp()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cut_spin_long);
memory->destroy(cutsq);
}
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2)
error->all(FLERR,"Incorrect args in pair_style command");
if (strcmp(update->unit_style,"metal") != 0)
error->all(FLERR,"Spin simulations require metal unit style");
cut_spin_long_global = force->numeric(FLERR,arg[0]);
// reset cutoffs that have been explicitly set
if (allocated) {
int i,j;
for (i = 1; i <= atom->ntypes; i++) {
for (j = i+1; j <= atom->ntypes; j++) {
if (setflag[i][j]) {
cut_spin_long[i][j] = cut_spin_long_global;
}
}
}
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::coeff(int narg, char **arg)
{
if (!allocated) allocate();
// check if args correct
if (strcmp(arg[2],"long") != 0)
error->all(FLERR,"Incorrect args in pair_style command");
if (narg < 1 || narg > 4)
error->all(FLERR,"Incorrect args in pair_style command");
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
double spin_long_cut_one = force->numeric(FLERR,arg[3]);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
setflag[i][j] = 1;
cut_spin_long[i][j] = spin_long_cut_one;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::init_style()
{
if (!atom->sp_flag)
error->all(FLERR,"Pair spin requires atom/spin style");
// need a full neighbor list
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
// checking if nve/spin is a listed fix
int ifix = 0;
while (ifix < modify->nfix) {
if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
ifix++;
}
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;
}
}
// insure use of KSpace long-range solver, set g_ewald
if (force->kspace == NULL)
error->all(FLERR,"Pair style requires a KSpace style");
g_ewald = force->kspace->g_ewald;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairSpinDipoleLongQsymp::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
cut_spin_long[j][i] = cut_spin_long[i][j];
return cut_spin_long_global;
}
/* ----------------------------------------------------------------------
extract the larger cutoff if "cut" or "cut_coul"
------------------------------------------------------------------------- */
void *PairSpinDipoleLongQsymp::extract(const char *str, int &dim)
{
if (strcmp(str,"cut") == 0) {
dim = 0;
return (void *) &cut_spin_long_global;
} else if (strcmp(str,"cut_coul") == 0) {
dim = 0;
return (void *) &cut_spin_long_global;
} else if (strcmp(str,"ewald_order") == 0) {
ewald_order = 0;
ewald_order |= 1<<1;
ewald_order |= 1<<3;
dim = 0;
return (void *) &ewald_order;
} else if (strcmp(str,"ewald_mix") == 0) {
dim = 0;
return (void *) &mix_flag;
}
return NULL;
}
/* ---------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double r,rinv,r2inv,rsq;
double grij,expm2,t,erfc,erf;
double sjdotr,gigj;
double bij[4];
double cij[4];
double evdwl,ecoul;
double xi[3],rij[3];
double spi[4],spj[4],fi[3],fmi[3];
double fmx_erf_s,fmy_erf_s,fmz_erf_s;
double local_cut2;
double pre1,pre2,pre3;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
double **fm = atom->fm;
double **fm_long = atom->fm_long;
double **sp = atom->sp;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
pre1 = 2.0 * g_ewald / MY_PIS;
pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS;
pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS;
// computation of the exchange interaction
// loop over atoms and their neighbors
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];
spi[3] = sp[i][3];
itype = type[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
evdwl = 0.0;
fi[0] = fi[1] = fi[2] = 0.0;
fmi[0] = fmi[1] = fmi[2] = 0.0;
fmx_erf_s = fmy_erf_s = fmz_erf_s = 0.0;
bij[0] = bij[1] = bij[2] = bij[3] = 0.0;
cij[0] = cij[1] = cij[2] = cij[3] = 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];
local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype];
if (rsq < local_cut2) {
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
erf = 1.0 - t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2;
// evaluating erfc for mech. force and energy calc.
bij[0] = erfc * rinv;
bij[1] = (bij[0] + pre1*expm2) * r2inv;
bij[2] = (3.0*bij[1] + pre2*expm2) * r2inv;
bij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
compute_long(i,j,rij,bij,fmi,spi,spj);
compute_long_mech(i,j,rij,bij,fmi,spi,spj);
// evaluating erf comp. for fm_kspace correction
cij[0] = erf * rinv;
cij[1] = (bij[0] + pre1*expm2) * r2inv;
cij[2] = (3.0*bij[1] + pre2*expm2) * r2inv;
//cij[3] = (5.0*bij[2] + pre3*expm2) * r2inv;
gigj = spi[3] * spj[3];
sjdotr = spj[0]*rij[0] + spj[1]*rij[1] + spj[2]*rij[2];
// evaluating short-range correction to the kspace part on [0,rc]
fmx_erf_s += gigj * (cij[2] * sjdotr * rij[0] - cij[1] * spj[0]);
fmy_erf_s += gigj * (cij[2] * sjdotr * rij[1] - cij[1] * spj[1]);
fmz_erf_s += gigj * (cij[2] * sjdotr * rij[2] - cij[1] * spj[2]);
}
// force accumulation
f[i][0] += fi[0] * mub2mu0;
f[i][1] += fi[1] * mub2mu0;
f[i][2] += fi[2] * mub2mu0;
fm[i][0] += fmi[0] * mub2mu0hbinv;
fm[i][1] += fmi[1] * mub2mu0hbinv;
fm[i][2] += fmi[2] * mub2mu0hbinv;
// correction of the fm_kspace
fm_long[i][0] -= (mub2mu0hbinv * fmx_erf_s);
fm_long[i][1] -= (mub2mu0hbinv * fmy_erf_s);
fm_long[i][2] -= (mub2mu0hbinv * fmz_erf_s);
if (newton_pair || j < nlocal) {
f[j][0] -= fi[0];
f[j][1] -= fi[1];
f[j][2] -= fi[2];
}
if (eflag) {
if (rsq <= local_cut2) {
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]);
}
}
}
/* ----------------------------------------------------------------------
update the pair interaction fmi acting on the spin ii
adding 1/r (for r in [0,rc]) contribution to the pair
removing erf(r)/r (for r in [0,rc]) from the kspace force
------------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::compute_single_pair(int ii, double fmi[3])
{
int i,j,jj,jnum,itype,jtype;
double rinv,r2inv,r3inv,rsq;
double sjdotr,sjdotrr3inv;
double b1,b2,gigj;
double bij[4],xi[3],rij[3];
double spi[4],spj[4];
double local_cut2;
int *ilist,*jlist,*numneigh,**firstneigh;
double fmx_s,fmy_s,fmz_s;
double **x = atom->x;
double **sp = atom->sp;
double **fm_long = atom->fm_long;
int *type = atom->type;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// computation of the exchange interaction
// loop over neighbors of atom i
i = ilist[ii];
xi[0] = x[i][0];
xi[1] = x[i][1];
xi[2] = x[i][2];
spi[3] = sp[i][3];
jlist = firstneigh[i];
jnum = numneigh[i];
itype = type[i];
fmx_s = fmy_s = fmz_s = 0.0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = type[j];
spj[0] = sp[j][0];
spj[1] = sp[j][1];
spj[2] = sp[j][2];
spj[3] = sp[j][3];
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];
local_cut2 = cut_spin_long[itype][jtype]*cut_spin_long[itype][jtype];
// evaluating full dipolar interaction on [0,rc]
if (rsq < local_cut2) {
r2inv = 1.0/rsq;
rinv = sqrt(r2inv);
r3inv = r2inv*rinv;
gigj = spi[3] * spj[3];
sjdotr = spj[0]*rij[0] + spj[1]*rij[1] + spj[2]*rij[2];
sjdotrr3inv = 3.0 * sjdotr * r3inv;
fmx_s += mub2mu0hbinv * gigj * (sjdotrr3inv * rij[0] - sp[i][0]/r3inv);
fmy_s += mub2mu0hbinv * gigj * (sjdotrr3inv * rij[1] - sp[i][1]/r3inv);
fmz_s += mub2mu0hbinv * gigj * (sjdotrr3inv * rij[2] - sp[i][2]/r3inv);
}
}
// adding truncated kspace force and short-range full force
fmi[0] += (fmx_s + fm_long[i][0]);
fmi[1] += (fmy_s + fm_long[i][1]);
fmi[2] += (fmz_s + fm_long[i][2]);
}
/* ----------------------------------------------------------------------
compute dipolar interaction between spins i and j
------------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::compute_long(int i, int j, double rij[3],
double bij[4], double fmi[3], double spi[4], double spj[4])
{
double sjdotr;
double b1,b2,gigj;
gigj = spi[3] * spj[3];
sjdotr = spj[0]*rij[0] + spj[1]*rij[1] + spj[2]*rij[2];
b1 = bij[1];
b2 = bij[2];
fmi[0] += mub2mu0hbinv * gigj * (b2 * sjdotr *rij[0] - b1 * spj[0]);
fmi[1] += mub2mu0hbinv * gigj * (b2 * sjdotr *rij[1] - b1 * spj[1]);
fmi[2] += mub2mu0hbinv * gigj * (b2 * sjdotr *rij[2] - b1 * spj[2]);
}
/* ----------------------------------------------------------------------
compute the mechanical force due to the dipolar interaction between
atom i and atom j
------------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::compute_long_mech(int i, int j, double rij[3],
double bij[4], double fi[3], double spi[3], double spj[3])
{
double sdots,sidotr,sjdotr,b2,b3;
double g1,g2,g1b2_g2b3,gigj;
gigj = spi[3] * spj[3];
sdots = spi[0]*spj[0] + spi[1]*spj[1] + spi[2]*spj[2];
sidotr = spi[0]*rij[0] + spi[1]*rij[1] + spi[2]*rij[2];
sjdotr = spj[0]*rij[0] + spj[1]*rij[1] + spj[2]*rij[2];
b2 = bij[2];
b3 = bij[3];
g1 = sdots;
g2 = -sidotr*sjdotr;
g1b2_g2b3 = g1*b2 + g2*b3;
fi[0] += gigj * (rij[0] * g1b2_g2b3 +
b2 * (sjdotr*spi[0] + sidotr*spj[0]));
fi[1] += gigj * (rij[1] * g1b2_g2b3 +
b2 * (sjdotr*spi[1] + sidotr*spj[1]));
fi[2] += gigj * (rij[2] * g1b2_g2b3 +
b2 * (sjdotr*spi[2] + sidotr*spj[2]));
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cut_spin_long,n+1,n+1,"pair/spin/long:cut_spin_long");
memory->create(cutsq,n+1,n+1,"pair/spin/long:cutsq");
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&cut_spin_long[i][j],sizeof(int),1,fp);
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++) {
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&cut_spin_long[i][j],sizeof(int),1,fp);
}
MPI_Bcast(&cut_spin_long[i][j],1,MPI_INT,0,world);
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::write_restart_settings(FILE *fp)
{
fwrite(&cut_spin_long_global,sizeof(double),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairSpinDipoleLongQsymp::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
fread(&cut_spin_long_global,sizeof(double),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
}
MPI_Bcast(&cut_spin_long_global,1,MPI_DOUBLE,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
}

100
src/SPIN/pair_spin_dipole_long_qsymp.h
View File

@ -1,100 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
www.cs.sandia.gov/~sjplimp/lammps.html
Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(spin/dipole/long/qsymp,PairSpinDipoleLongQsymp)
#else
#ifndef LMP_PAIR_SPIN_DIPOLE_LONG_QSYMP_H
#define LMP_PAIR_SPIN_DIPOLE_LONG_QSYMP_H
#include "pair_spin.h"
namespace LAMMPS_NS {
class PairSpinDipoleLongQsymp : public PairSpin {
public:
double cut_coul;
double **sigma;
PairSpinDipoleLongQsymp(class LAMMPS *);
virtual ~PairSpinDipoleLongQsymp();
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
void init_style();
void *extract(const char *, int &);
void compute(int, int);
void compute_single_pair(int, double *);
void compute_long(int, int, double *, double *, double *,
double *, double *);
void compute_long_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_long_global; // global long cutoff distance
protected:
double hbar; // reduced Planck's constant
double mub; // Bohr's magneton
double mu_0; // vacuum permeability
double mub2mu0; // prefactor for mech force
double mub2mu0hbinv; // prefactor for mag force
double **cut_spin_long; // cutoff distance long
double g_ewald;
int ewald_order;
int lattice_flag; // flag for mech force computation
class FixNVESpin *lockfixnvespin; // ptr for setups
void allocate();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Incorrect args in pair_style command
Self-explanatory.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair dipole/long requires atom attributes q, mu, torque
The atom style defined does not have these attributes.
E: Cannot (yet) use 'electron' units with dipoles
This feature is not yet supported.
E: Pair style requires a KSpace style
No kspace style is defined.
*/

2
src/SPIN/pair_spin_exchange.cpp
View File

@ -75,7 +75,7 @@ PairSpinExchange::~PairSpinExchange()
void PairSpinExchange::settings(int narg, char **arg)
{
if (narg < 1 || narg > 2)
if (narg < 1 || narg > 7)
error->all(FLERR,"Incorrect number of args in pair_style pair/spin command");
if (strcmp(update->unit_style,"metal") != 0)