forked from lijiext/lammps
Commit Julien 10/19/2017
- New files for the pair interactions - New files for the documentation - Spin orbit coupling via Neel approach
This commit is contained in:
julient31
parent
54832a8fe4
commit
72d9795d7f
|
|
@ -0,0 +1,63 @@
|
|||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
compute spin command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
compute ID group-ID spin :pre
|
||||
|
||||
ID, group-ID are documented in "compute"_compute.html command
|
||||
spin = style name of this compute command
|
||||
|
||||
[Examples:]
|
||||
|
||||
compute out_mag all spin
|
||||
|
||||
[Description:]
|
||||
|
||||
Define a computation that calculates the magnetic quantities for a system
|
||||
of atoms having spins.
|
||||
|
||||
This compute calculates 6 magnetic quantities.
|
||||
The three first ones are the x,y and z cooredinates of the total magnetization.
|
||||
The fourth quantity is the norm of the total magnetization.
|
||||
The fifth one is referred to as the spin temperature, according
|
||||
to the work of "(Nurdin)"_#Nurdin1.
|
||||
|
||||
The sixth quantity is the magnetic energy.
|
||||
|
||||
|
||||
A simple way to output the results of the compute spin
|
||||
calculation to a file is to use the "fix ave/time"_fix_ave_time.html
|
||||
command, for example:
|
||||
|
||||
compute mag all compute/spin
|
||||
fix outmag all ave/time 1 1 50 c_mag[2] c_mag[3] c_mag[4] c_mag[7] file mag.dat
|
||||
|
||||
[Output info:]
|
||||
|
||||
The array values are "intensive". The array values will be in
|
||||
metal units ("units"_units.html).
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
The {spin} compute is part of the SPIN package.
|
||||
This compute is only enabled if LAMMPS was built with this package.
|
||||
See the "Making LAMMPS"_Section_start.html#start_3 section for more info.
|
||||
|
||||
[Related commands:] none
|
||||
|
||||
|
||||
[Default:] none
|
||||
|
||||
:line
|
||||
|
||||
:link(Nurdin1)
|
||||
[(Nurdin)] Nurdin and Schotte Phys Rev E, 61(4), 3579 (2000)
|
||||
|
||||
|
|
@ -0,0 +1,82 @@
|
|||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
fix force/spin command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
fix ID group force/spin style args :pre
|
||||
|
||||
ID, group are documented in "fix"_fix.html command :ulb,l
|
||||
force/spin = style name of this fix command :l
|
||||
style = {zeeman} or {aniso} :l
|
||||
{zeeman} args = H x y z
|
||||
H = intensity of the magnetic field (in Tesla)
|
||||
x y z = vector direction of the field
|
||||
{aniso} args = K x y z
|
||||
K = intensity of the magnetic anisotropy (in eV)
|
||||
x y z = vector direction of the anisotropy
|
||||
:ule
|
||||
|
||||
[Examples:]
|
||||
|
||||
fix 1 all force/spin zeeman 0.1 0.0 0.0 1.0
|
||||
fix 1 all force/spin aniso 0.001 0.0 0.0 1.0
|
||||
|
||||
|
||||
[Description:]
|
||||
|
||||
Impose a force torque to each magnetic spin in the group.
|
||||
|
||||
Style {zeeman} is used for the simulation of the interaction
|
||||
between the magnetic spins in the defined group and an external
|
||||
magnetic field.
|
||||
|
||||
Style {aniso} is used to simulate an easy axis or an easy plane
|
||||
for the magnetic spins in the defined group.
|
||||
If K>0, an easy axis is defined, and if K<0, an easy plane is
|
||||
defined.
|
||||
|
||||
In both cases, x y z impose is the vector direction for the force.
|
||||
Only the direction of the vector is important; it's length is ignored.
|
||||
|
||||
:line
|
||||
|
||||
|
||||
[Restart, fix_modify, output, run start/stop, minimize info:]
|
||||
|
||||
No information about this fix is written to "binary restart
|
||||
files"_restart.html.
|
||||
|
||||
The "fix_modify"_fix_modify.html {energy} option is supported by this
|
||||
fix to add the gravitational potential energy of the system to the
|
||||
system's potential energy as part of "thermodynamic
|
||||
output"_thermo_style.html.
|
||||
|
||||
The "fix_modify"_fix_modify.html {respa} option is supported by this
|
||||
fix. This allows to set at which level of the "r-RESPA"_run_style.html
|
||||
integrator the fix is adding its forces. Default is the outermost level.
|
||||
|
||||
This fix computes a global scalar which can be accessed by various
|
||||
"output commands"_Section_howto.html#howto_15. This scalar is the
|
||||
gravitational potential energy of the particles in the defined field,
|
||||
namely mass * (g dot x) for each particles, where x and mass are the
|
||||
particles position and mass, and g is the gravitational field. The
|
||||
scalar value calculated by this fix is "extensive".
|
||||
|
||||
No parameter of this fix can be used with the {start/stop} keywords of
|
||||
the "run"_run.html command. This fix is not invoked during "energy
|
||||
minimization"_minimize.html.
|
||||
|
||||
[Restrictions:] none
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"atom_style sphere"_atom_style.html, "fix addforce"_fix_addforce.html
|
||||
|
||||
[Default:] none
|
||||
|
|
@ -0,0 +1,65 @@
|
|||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
fix integration/spin command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
fix ID group-ID integration/spin style :pre
|
||||
|
||||
ID, group-ID are documented in "fix"_fix.html command :ulb,l
|
||||
integration/spin = style name of this fix command :l
|
||||
style = {serial} or {mpi} :l
|
||||
{serial} value = factor
|
||||
factor = do thermostat rotational degrees of freedom via the angular momentum and apply numeric scale factor as discussed below :pre
|
||||
{mpi}
|
||||
:ule
|
||||
|
||||
[Examples:]
|
||||
|
||||
fix 3 all integration/spin serial
|
||||
fix 1 all integration/spin mpi
|
||||
|
||||
[Description:]
|
||||
|
||||
A Suzuki-Trotter decomposition is applied to the integration of the
|
||||
spin-lattice system.
|
||||
|
||||
:line
|
||||
|
||||
The {style} value defines if a serial of a parallel
|
||||
algorithm has to be used for the integration.
|
||||
|
||||
The parallel algorithm uses a sectoring method as
|
||||
described in .
|
||||
|
||||
:line
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
This fix style can only be used if LAMMPS was built with the
|
||||
SPIN package. See the "Making LAMMPS"_Section_start.html#start_3
|
||||
section for more info on packages.
|
||||
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"fix nve"_fix_nve.html, "pair spin"_pair_spin.html,
|
||||
"compute spin"_compute_spin.html, "fix langevin spin"_fix_langevin_spin.html,
|
||||
|
||||
|
||||
[Default:] none
|
||||
|
||||
:line
|
||||
|
||||
:link(Davidchack2)
|
||||
[(Davidchack)] R.L Davidchack, T.E. Ouldridge, M.V. Tretyakov. J. Chem. Phys. 142, 144114 (2015).
|
||||
:link(Miller2)
|
||||
[(Miller)] T. F. Miller III, M. Eleftheriou, P. Pattnaik, A. Ndirango, G. J. Martyna, J. Chem. Phys., 116, 8649-8659 (2002).
|
||||
:link(Dunweg3)
|
||||
[(Dunweg)] B. Dunweg, W. Paul, Int. J. Mod. Phys. C, 2, 817-27 (1991).
|
||||
|
|
@ -0,0 +1,91 @@
|
|||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
fix langevin/spin command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
fix ID group-ID langevin/spin T Tdamp Ldamp seed :pre
|
||||
|
||||
ID, group-ID are documented in "fix"_fix.html command :ulb,l
|
||||
langevin/spin = style name of this fix command :l
|
||||
T = desired temperature of the bath (temperature units, K in metal units) :l
|
||||
Tdamp = transverse magnetic damping parameter (adim) :l
|
||||
Ldamp = longitudinal magnetic damping parameter (adim) :l
|
||||
seed = random number seed to use for white noise (positive integer) :l
|
||||
:ule
|
||||
|
||||
[Examples:]
|
||||
|
||||
fix 2 all langevin/spin 300.0 0.01 0.0 21 :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
Apply a Langevin thermostat as described in "(Mayergoyz)"_#Mayergoyz1
|
||||
to the magnetic spins associated to the atoms.
|
||||
Used with "fix integration spin"_fix_integration_spin.html, this
|
||||
command performs Brownian dynamics (BD), apply a random torque
|
||||
and a transverse dissipation to each spin:
|
||||
|
||||
Rand torque =
|
||||
Transverse dmping =
|
||||
D is proportional to sqrt(Kb T m / (hbar dt damp)) :pre
|
||||
|
||||
|
||||
Note: The random # {seed} must be a positive integer. A Marsaglia random
|
||||
number generator is used. Each processor uses the input seed to
|
||||
generate its own unique seed and its own stream of random numbers.
|
||||
Thus the dynamics of the system will not be identical on two runs on
|
||||
different numbers of processors.
|
||||
|
||||
:line
|
||||
|
||||
[Restart, fix_modify, output, run start/stop, minimize info:]
|
||||
|
||||
No information about this fix is written to "binary restart
|
||||
files"_restart.html. Because the state of the random number generator
|
||||
is not saved in restart files, this means you cannot do "exact"
|
||||
restarts with this fix, where the simulation continues on the same as
|
||||
if no restart had taken place. However, in a statistical sense, a
|
||||
restarted simulation should produce the same behavior.
|
||||
|
||||
This fix is not invoked during "energy minimization"_minimize.html.
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
The {langevin/spin} fix is part of the SPIN package.
|
||||
These styles are only enabled if LAMMPS was built with this package.
|
||||
See the "Making LAMMPS"_Section_start.html#start_3 section for more info.
|
||||
|
||||
The numerical integration has to be performed with {fix/integration/spin}
|
||||
when {langevin/spin} is enabled.
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"fix integration spin"_fix_integration_spin.html,
|
||||
"pair spin"_pair_spin.html
|
||||
|
||||
[Default:]
|
||||
|
||||
The option defaults are angmom = no, omega = no, scale = 1.0 for all
|
||||
types, tally = no, zero = no, gjf = no.
|
||||
|
||||
:line
|
||||
|
||||
:link(Mayergoyz1)
|
||||
Mayergoyz, Bertotti, and Serpico (2009). Elsevier (2009)
|
||||
|
||||
|
||||
|
||||
:link(Schneider1)
|
||||
[(Schneider)] Schneider and Stoll, Phys Rev B, 17, 1302 (1978).
|
||||
|
||||
:link(Gronbech-Jensen)
|
||||
[(Gronbech-Jensen)] Gronbech-Jensen and Farago, Mol Phys, 111, 983
|
||||
(2013); Gronbech-Jensen, Hayre, and Farago, Comp Phys Comm,
|
||||
185, 524 (2014)
|
||||
|
|
@ -0,0 +1,90 @@
|
|||
<script type="text/javascript"
|
||||
src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
|
||||
</script>
|
||||
<script type="text/x-mathjax-config">
|
||||
MathJax.Hub.Config({ TeX: { equationNumbers: {autoNumber: "AMS"} } });
|
||||
</script>
|
||||
|
||||
|
||||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Section_commands.html#comm)
|
||||
|
||||
:line
|
||||
|
||||
pair_style pair/spin/exchange command :h3
|
||||
pair_style pair/spin/dmi command :h3
|
||||
pair_style pair/spin/me command :h3
|
||||
pair_style pair/spin/soc command :h3
|
||||
|
||||
|
||||
[Syntax:]
|
||||
|
||||
pair_style pair/spin/exchange cutoff
|
||||
pair_style pair/spin/dmi cutoff
|
||||
pair_style pair/spin/me cutoff
|
||||
pair_style pair/spin/soc/neel cutoff :pre
|
||||
|
||||
cutoff = global cutoff pair (distance units) :ulb,l
|
||||
|
||||
:ule
|
||||
|
||||
[Examples:]
|
||||
|
||||
pair_style pair/spin/exchange 4.0
|
||||
pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885
|
||||
pair_coeff 1 2 exchange 4.0 0.0446928 0.003496 1.4885 :pre
|
||||
|
||||
pair_style pair/spin/dmi 4.0
|
||||
pair_coeff * * dmi 2.6 0.001 0.0 0.0 1.0 :pre
|
||||
|
||||
pair_style pair/spin/soc/neel 4.0
|
||||
pair_coeff * * neel 4.0 -0.003330282 0.864159 2.13731 :pre
|
||||
|
||||
|
||||
[Description:]
|
||||
|
||||
Style {pair/spin} computes interactions between pairs of particles
|
||||
that each have a spin.
|
||||
|
||||
|
||||
The {pair/spin} style compute the short-range spin-spin interactions.
|
||||
|
||||
|
||||
\begin\{equation\}
|
||||
J_{i=0} = \sum_i \alpha
|
||||
\end\{equation\}
|
||||
$$\sum_{i=0}^n i^2 = \frac{(n^2+n)(2n+1)}{6}$$
|
||||
|
||||
|
||||
:c,image(Eqs/pair_lj.jpg)
|
||||
|
||||
The {lj/cut} and {lj/cut/coul/long} pair styles support the use of the
|
||||
{inner}, {middle}, and {outer} keywords of the "run_style
|
||||
respa"_run_style.html command, meaning the pairwise forces can be
|
||||
|
||||
:line
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
All the {pair/spin} styles are part of the SPIN package.
|
||||
These styles are only enabled if LAMMPS was built with this package.
|
||||
See the "Making LAMMPS"_Section_start.html#start_3 section for more info.
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"eam"
|
||||
|
||||
"pair_coeff"_pair_coeff.html
|
||||
|
||||
[Default:] none
|
||||
|
||||
:line
|
||||
|
||||
:link(Beaujouan1)
|
||||
[(Beaujouan)] Beaujouan, Thibaudeau, and Barreteau, Phys Rev B, 86(17), 174409 (2012)
|
||||
|
||||
:link(Perera2)
|
||||
[(Perera)] Perera, Eisenbach, Nicholson, Stocks, and Landau, Phys Rev B, 93(6), 060402 (2016)
|
||||
Binary file not shown.
|
|
@ -48,75 +48,95 @@ create_atoms 1 box
|
|||
#######Settings########
|
||||
#######################
|
||||
|
||||
#isolating 1 atom into a group
|
||||
group single_spin id 10
|
||||
|
||||
#Setting one or more properties of one or more atoms.
|
||||
#Setting mass
|
||||
mass 1 1.0
|
||||
#set group all mass 1.0
|
||||
mass 1 58.93
|
||||
#set group all mass 58.93
|
||||
|
||||
#Setting spins orientation and moment
|
||||
set group all spin/random 11 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
|
||||
#set group single_spin spin/random 11 1.72
|
||||
|
||||
#velocity all create 600 4928459 rot yes dist gaussian
|
||||
|
||||
#Magneto-mechanic interactions for bulk fcc Cobalt
|
||||
pair_style hybrid/overlay eam/alloy pair/spin 4.0
|
||||
#pair_style pair/spin 4.0
|
||||
#pair_style eam/alloy
|
||||
#pair_style hybrid/overlay eam/alloy pair/spin/soc 4.0
|
||||
pair_style hybrid/overlay eam/alloy pair/spin/exchange 4.0 pair/spin/soc/neel 4.0
|
||||
#pair_style hybrid/overlay eam/alloy pair/spin 4.0 pair/spin/soc 4.0
|
||||
|
||||
pair_coeff * * eam/alloy ../Co_PurjaPun_2012.eam.alloy Co
|
||||
#pair_coeff * * ../Co_PurjaPun_2012.eam.alloy Co
|
||||
|
||||
#pair_style pair/spin 4.0
|
||||
#type i and j | interaction type | cutoff | J1 (eV) | J2 (adim) | J3 (Ang) (for Exchange)
|
||||
pair_coeff * * pair/spin exchange 4.0 0.0446928 0.003496 1.4885
|
||||
#pair_coeff * * exchange 4.0 0.0446928 0.003496 1.4885
|
||||
#pair_coeff * * exchange 4.0 0.0 0.003496 1.4885
|
||||
pair_coeff * * pair/spin/exchange exchange 4.0 0.0446928 0.003496 1.4885
|
||||
#pair_coeff * * pair/spin exchange 4.0 0.0 0.003496 1.4885
|
||||
|
||||
#type i and j | interaction type | cutoff | Int (eV) | [dx,dy,dz] (for DMI and ME)
|
||||
#pair_coeff * * dmi 2.6 0.001 0.0 0.0 1.0
|
||||
#pair_coeff * * me 2.6 0.01 1.0 1.0 1.0
|
||||
|
||||
#type i and j | interaction type | cutoff | K1 (eV) | K2 (adim) | K3 (Ang) (for SOC)
|
||||
pair_coeff * * pair/spin/soc/neel neel 4.0 -0.003330282 0.864159 2.13731
|
||||
|
||||
#Define a skin distance, update neigh list every
|
||||
#neighbor 1.0 bin
|
||||
#neigh_modify every 10 check yes delay 20
|
||||
neighbor 0.0 bin
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 check no delay 0
|
||||
|
||||
#Magnetic field fix
|
||||
#Type | Intensity (T or eV) | Direction
|
||||
fix 1 all force/spin zeeman 0.0 0.0 0.0 1.0
|
||||
#fix 1 all force/spin anisotropy 0.01 0.0 0.0 1.0
|
||||
#fix 1 all force/spin anisotropy -0.1 0.0 0.0 1.0
|
||||
#fix 1 all force/spin anisotropy 0.1 0.0 1.0 0.0
|
||||
|
||||
#Fix Langevin spins (merging damping and temperature)
|
||||
#Temp | Alpha_trans | Alpha_long | Seed
|
||||
#fix 2 all langevin/spin 0.0 0.1 0.0 21
|
||||
#fix 2 all langevin/spin 1.0 0.1 0.0 21
|
||||
fix 2 all langevin/spin 0.0 0.0 0.0 21
|
||||
fix 2 all langevin/spin 300.0 0.1 0.0 21
|
||||
|
||||
#Magnetic integration fix
|
||||
fix 3 all integration/spin serial
|
||||
|
||||
#compute real time, total magnetization, magnetic energy, and spin temperature
|
||||
#Iteration | Time | Mx | My | Mz | |M| | Em | Tm
|
||||
compute mag all compute/spin
|
||||
fix outmag all ave/time 1 1 10 c_mag[1] c_mag[2] c_mag[3] c_mag[4] c_mag[5] c_mag[6] c_mag[7] file mag_Co_nodamp.dat format %20.16g
|
||||
|
||||
#compute mechanical energy
|
||||
compute mech all pe
|
||||
fix outmech all ave/time 1 1 10 c_mech file mech_Co_nodamp.dat format %20.16g
|
||||
|
||||
#compute kinetic energy
|
||||
compute kinetic all ke
|
||||
fix outke all ave/time 1 1 10 c_kinetic file kinetic_Co_nodamp.dat format %20.16g
|
||||
#compute mag all compute/spin
|
||||
#fix outmag all ave/time 1 1 50 c_mag[1] c_mag[2] c_mag[3] c_mag[4] c_mag[5] c_mag[6] c_mag[7] file mag_VSRSV.dat format %20.16g
|
||||
|
||||
#Setting the timestep for the simulation (in ps)
|
||||
timestep 0.0001
|
||||
timestep 0.0002
|
||||
|
||||
##################
|
||||
#######run########
|
||||
##################
|
||||
|
||||
thermo 500
|
||||
#fix_modify 1 energy yes
|
||||
|
||||
compute out_mag all compute/spin
|
||||
compute out_pe all pe
|
||||
compute out_ke all ke
|
||||
|
||||
variable magnorm equal c_out_mag[5]
|
||||
variable emag equal c_out_mag[6]
|
||||
variable tmag equal c_out_mag[7]
|
||||
variable mag_force equal f_1
|
||||
|
||||
thermo 500
|
||||
thermo_style custom step time v_emag c_out_pe c_out_ke v_mag_force v_magnorm v_tmag etotal
|
||||
#fix out_vals all ave/time 1 1 50 step v_emag file temp_lattice_VSRSV.dat format %20.16g
|
||||
|
||||
#Dump the positions and spin directions of magnetic particles (vmd format)
|
||||
dump 1 all custom 10 dump_spin_MM.lammpstrj type x y z spx spy spz
|
||||
dump 1 all custom 100 dump_VSRSV.lammpstrj type x y z spx spy spz
|
||||
|
||||
#Running the simulations for N timesteps
|
||||
run 100000
|
||||
run 1000000
|
||||
#run 1
|
||||
|
||||
|
|
|
|||
|
|
@ -100,14 +100,14 @@ void ComputeSpin::compute_vector()
|
|||
mag[0] += sp[i][0];
|
||||
mag[1] += sp[i][1];
|
||||
mag[2] += sp[i][2];
|
||||
magenergy += sp[i][0]*fm[i][0];
|
||||
magenergy += sp[i][1]*fm[i][1];
|
||||
magenergy += sp[i][2]*fm[i][2];
|
||||
magenergy -= sp[i][0]*fm[i][0];
|
||||
magenergy -= sp[i][1]*fm[i][1];
|
||||
magenergy -= 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]*sp[i][0]+sp[i][1]*sp[i][1]+sp[i][2]*sp[i][2];
|
||||
tempdenom += sp[i][0]*fm[i][0]+fm[i][1]*sp[i][1]+sp[i][2]*fm[i][2];
|
||||
countsp++;
|
||||
}
|
||||
}
|
||||
|
|
@ -124,15 +124,16 @@ void ComputeSpin::compute_vector()
|
|||
magtot[0] *= scale;
|
||||
magtot[1] *= scale;
|
||||
magtot[2] *= scale;
|
||||
magtot[3] = sqrt(square(magtot[0])+square(magtot[1])+square(magtot[2]));
|
||||
spintemperature = hbar*tempnumtot/2.0/kb/tempdenomtot;
|
||||
magtot[3] = sqrt((magtot[0]*magtot[0])+(magtot[1]*magtot[1])+(magtot[2]*magtot[2]));
|
||||
spintemperature = hbar*tempnumtot;
|
||||
spintemperature /= (2.0*kb*tempdenomtot);
|
||||
|
||||
vector[0] = invoked_vector*update->dt;
|
||||
vector[1] = magtot[0];
|
||||
vector[2] = magtot[1];
|
||||
vector[3] = magtot[2];
|
||||
vector[4] = magtot[3];
|
||||
vector[5] = -0.5*magenergytot*hbar;
|
||||
vector[5] = magenergytot*hbar;
|
||||
vector[6] = spintemperature;
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -65,10 +65,12 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a
|
|||
magfieldstyle = CONSTANT;
|
||||
|
||||
H_field = 0.0;
|
||||
Hx = Hy = Hz = 0.0;
|
||||
nhx = nhy = nhz = 0.0;
|
||||
hx = hy = hz = 0.0;
|
||||
Ka = 0.0;
|
||||
nax = nay = naz = 0.0;
|
||||
Kax = Kay = Kaz = 0.0;
|
||||
|
||||
|
||||
zeeman_flag = aniso_flag = 0;
|
||||
|
||||
if (strcmp(arg[3],"zeeman") == 0) {
|
||||
|
|
@ -76,25 +78,25 @@ FixForceSpin::FixForceSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, a
|
|||
style = ZEEMAN;
|
||||
zeeman_flag = 1;
|
||||
H_field = force->numeric(FLERR,arg[4]);
|
||||
Hx = force->numeric(FLERR,arg[5]);
|
||||
Hy = force->numeric(FLERR,arg[6]);
|
||||
Hz = force->numeric(FLERR,arg[7]);
|
||||
nhx = force->numeric(FLERR,arg[5]);
|
||||
nhy = force->numeric(FLERR,arg[6]);
|
||||
nhz = force->numeric(FLERR,arg[7]);
|
||||
magfieldstyle = CONSTANT;
|
||||
} else if (strcmp(arg[3],"anisotropy") == 0) {
|
||||
if (narg != 8) error->all(FLERR,"Illegal force/spin command");
|
||||
style = ANISOTROPY;
|
||||
aniso_flag = 1;
|
||||
Ka = force->numeric(FLERR,arg[4]);
|
||||
Kax = force->numeric(FLERR,arg[5]);
|
||||
Kay = force->numeric(FLERR,arg[6]);
|
||||
Kaz = force->numeric(FLERR,arg[7]);
|
||||
nax = force->numeric(FLERR,arg[5]);
|
||||
nay = force->numeric(FLERR,arg[6]);
|
||||
naz = force->numeric(FLERR,arg[7]);
|
||||
} else error->all(FLERR,"Illegal force/spin command");
|
||||
|
||||
degree2rad = MY_PI/180.0;
|
||||
time_origin = update->ntimestep;
|
||||
|
||||
eflag = 0;
|
||||
emag = 0.0;
|
||||
emag = 0.0;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -191,22 +193,23 @@ void FixForceSpin::post_force(int vflag)
|
|||
fmi[0] = fmi[1] = fmi[2] = 0.0;
|
||||
if (zeeman_flag) {
|
||||
compute_zeeman(i,fmi);
|
||||
emag -= sp[i][0] * fmi[0];
|
||||
emag -= sp[i][1] * fmi[1];
|
||||
emag -= sp[i][2] * fmi[2];
|
||||
}
|
||||
if (aniso_flag) {
|
||||
spi[0] = sp[i][0];
|
||||
spi[1] = sp[i][1];
|
||||
spi[2] = sp[i][2];
|
||||
compute_anisotropy(i,spi,fmi);
|
||||
emag -= 0.5 * sp[i][0] * fmi[0];
|
||||
emag -= 0.5 * sp[i][1] * fmi[1];
|
||||
emag -= 0.5 * sp[i][2] * fmi[2];
|
||||
}
|
||||
fm[i][0] += fmi[0];
|
||||
fm[i][1] += fmi[1];
|
||||
fm[i][2] += fmi[2];
|
||||
|
||||
emag -= mumag[i]*spi[0]*fmi[0];
|
||||
emag -= mumag[i]*spi[1]*fmi[1];
|
||||
emag -= mumag[i]*spi[2]*fmi[2];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -214,19 +217,19 @@ void FixForceSpin::post_force(int vflag)
|
|||
void FixForceSpin::compute_zeeman(int i, double *fmi)
|
||||
{
|
||||
double *mumag = atom->mumag;
|
||||
fmi[0] += mumag[i]*xmag;
|
||||
fmi[1] += mumag[i]*ymag;
|
||||
fmi[2] += mumag[i]*zmag;
|
||||
fmi[0] += mumag[i]*hx;
|
||||
fmi[1] += mumag[i]*hy;
|
||||
fmi[2] += mumag[i]*hz;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixForceSpin::compute_anisotropy(int i, double * spi, double *fmi)
|
||||
{
|
||||
double scalar = Kax*spi[0] + Kay*spi[1] + Kaz*spi[2];
|
||||
fmi[0] += scalar*xmag;
|
||||
fmi[1] += scalar*ymag;
|
||||
fmi[2] += scalar*zmag;
|
||||
double scalar = nax*spi[0] + nay*spi[1] + naz*spi[2];
|
||||
fmi[0] += scalar*Kax;
|
||||
fmi[1] += scalar*Kay;
|
||||
fmi[2] += scalar*Kaz;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -240,14 +243,14 @@ void FixForceSpin::post_force_respa(int vflag, int ilevel, int iloop)
|
|||
void FixForceSpin::set_magneticforce()
|
||||
{
|
||||
if (style == ZEEMAN) {
|
||||
xmag = H_field*Hx;
|
||||
ymag = H_field*Hy;
|
||||
zmag = H_field*Hz;
|
||||
hx = H_field*nhx;
|
||||
hy = H_field*nhy;
|
||||
hz = H_field*nhz;
|
||||
}
|
||||
if (style == ANISOTROPY) {
|
||||
xmag = 2.0*Ka*Kax;
|
||||
ymag = 2.0*Ka*Kay;
|
||||
zmag = 2.0*Ka*Kaz;
|
||||
Kax = 2.0*Ka*nax;
|
||||
Kay = 2.0*Ka*nay;
|
||||
Kaz = 2.0*Ka*naz;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -257,13 +260,14 @@ void FixForceSpin::set_magneticforce()
|
|||
|
||||
double FixForceSpin::compute_scalar()
|
||||
{
|
||||
double emag_all=0.0;
|
||||
// only sum across procs one time
|
||||
if (eflag == 0) {
|
||||
MPI_Allreduce(&emag,&emag_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
eflag = 1;
|
||||
}
|
||||
|
||||
emag *= hbar;
|
||||
|
||||
emag_all *= hbar;
|
||||
|
||||
return emag_all;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -44,13 +44,12 @@ class FixForceSpin : public Fix {
|
|||
protected:
|
||||
int style; // style of the magnetic force
|
||||
|
||||
double xmag, ymag, zmag; // temp. force variables
|
||||
double degree2rad;
|
||||
double hbar;
|
||||
int ilevel_respa;
|
||||
int time_origin;
|
||||
int eflag;
|
||||
double emag, emag_all;
|
||||
double emag;
|
||||
|
||||
int varflag;
|
||||
int magfieldstyle;
|
||||
|
|
@ -59,11 +58,13 @@ class FixForceSpin : public Fix {
|
|||
|
||||
// zeeman field intensity and direction
|
||||
double H_field;
|
||||
double Hx, Hy, Hz;
|
||||
double nhx, nhy, nhz;
|
||||
double hx, hy, hz; // temp. force variables
|
||||
|
||||
// magnetic anisotropy intensity and direction
|
||||
double Ka;
|
||||
double Kax, Kay, Kaz;
|
||||
double nax, nay, naz;
|
||||
double Kax, Kay, Kaz; // temp. force variables
|
||||
|
||||
// temp. spin variables
|
||||
double *spi, *fmi;
|
||||
|
|
|
|||
|
|
@ -37,6 +37,11 @@
|
|||
#include "pair.h"
|
||||
#include "pair_hybrid.h"
|
||||
#include "pair_spin.h"
|
||||
#include "pair_spin_exchange.h"
|
||||
#include "pair_spin_me.h"
|
||||
#include "pair_spin_soc_dmi.h"
|
||||
#include "pair_spin_soc_neel.h"
|
||||
#include "pair_spin_soc_landau.h"
|
||||
#include "respa.h"
|
||||
#include "update.h"
|
||||
|
||||
|
|
@ -74,13 +79,16 @@ FixIntegrationSpin::FixIntegrationSpin(LAMMPS *lmp, int narg, char **arg) :
|
|||
if (extra == SPIN && !atom->mumag_flag)
|
||||
error->all(FLERR,"Fix integration/spin requires spin attribute mumag");
|
||||
magpair_flag = 0;
|
||||
exch_flag = dmi_flag = me_flag = 0;
|
||||
soc_flag = 0;
|
||||
exch_flag = 0;
|
||||
magforce_flag = 0;
|
||||
zeeman_flag = aniso_flag = 0;
|
||||
maglangevin_flag = 0;
|
||||
tdamp_flag = temp_flag = 0;
|
||||
|
||||
lockpairspin = NULL;
|
||||
lockpairspinexchange = NULL;
|
||||
lockpairspinsocneel = NULL;
|
||||
lockforcespin = NULL;
|
||||
locklangevinspin = NULL;
|
||||
}
|
||||
|
|
@ -95,6 +103,7 @@ FixIntegrationSpin::~FixIntegrationSpin(){
|
|||
memory->destroy(rsec);
|
||||
memory->destroy(seci);
|
||||
|
||||
memory->destroy(rij);
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fmi);
|
||||
|
|
@ -118,32 +127,39 @@ void FixIntegrationSpin::init()
|
|||
//FixNVE::init();
|
||||
|
||||
// set timesteps
|
||||
dtv = update->dt;
|
||||
dtv = 0.5 * update->dt;
|
||||
dtf = 0.5 * update->dt * force->ftm2v;
|
||||
dts = update->dt;
|
||||
dts = 0.5 * update->dt;
|
||||
|
||||
memory->create(xi,3,"integrations:xi");
|
||||
memory->create(sec,3,"integrations:sec");
|
||||
memory->create(rsec,3,"integrations:rsec");
|
||||
memory->create(seci,3,"integrations:seci");
|
||||
|
||||
memory->create(rij,3,"integrations:xi");
|
||||
memory->create(spi,3,"integrations:spi");
|
||||
memory->create(spj,3,"integrations:spj");
|
||||
memory->create(fmi,3,"integrations:fmi");
|
||||
memory->create(fmj,3,"integrations:fmj");
|
||||
|
||||
|
||||
if (strstr(force->pair_style,"pair/spin")) {
|
||||
magpair_flag = 1;
|
||||
lockpairspin = (PairSpin *) force->pair;
|
||||
if (strstr(force->pair_style,"pair/spin/exchange")) {
|
||||
exch_flag = 1;
|
||||
lockpairspinexchange = (PairSpinExchange *) force->pair;
|
||||
} else if (strstr(force->pair_style,"pair/spin/soc")) {
|
||||
soc_flag = 1;
|
||||
lockpairspinsocneel = (PairSpinSocNeel *) force->pair;
|
||||
} else if (strstr(force->pair_style,"hybrid/overlay")) {
|
||||
PairHybrid *lockhybrid = (PairHybrid *) force->pair;
|
||||
int nhybrid_styles = lockhybrid->nstyles;
|
||||
int ipair;
|
||||
for (ipair = 0; ipair < nhybrid_styles; ipair++) {
|
||||
if (strstr(lockhybrid->keywords[ipair],"pair/spin")) {
|
||||
magpair_flag = 1;
|
||||
lockpairspin = (PairSpin *) lockhybrid->styles[ipair];
|
||||
if (strstr(lockhybrid->keywords[ipair],"pair/spin/exchange")) {
|
||||
exch_flag = 1;
|
||||
lockpairspinexchange = (PairSpinExchange *) lockhybrid->styles[ipair];
|
||||
} else if (strstr(lockhybrid->keywords[ipair],"pair/spin/soc")) {
|
||||
soc_flag = 1;
|
||||
lockpairspinsocneel = (PairSpinSocNeel *) lockhybrid->styles[ipair];
|
||||
}
|
||||
}
|
||||
} else error->all(FLERR,"Illegal fix integration/spin command");
|
||||
|
|
@ -165,12 +181,10 @@ void FixIntegrationSpin::init()
|
|||
}
|
||||
}
|
||||
|
||||
// set flags for the different magnetic interactionsi
|
||||
if (magpair_flag) {
|
||||
if (lockpairspin->exch_flag == 1) exch_flag = 1;
|
||||
if (lockpairspin->dmi_flag == 1) dmi_flag = 1;
|
||||
if (lockpairspin->me_flag == 1) me_flag = 1;
|
||||
}
|
||||
// set flags for the different magnetic interactions
|
||||
// if (magpair_flag) {
|
||||
// if (lockpairspinexchange->exch_flag == 1) exch_flag = 1;
|
||||
// }
|
||||
|
||||
if (magforce_flag) {
|
||||
if (lockforcespin->zeeman_flag == 1) zeeman_flag = 1;
|
||||
|
|
@ -206,16 +220,130 @@ void FixIntegrationSpin::initial_integrate(int vflag)
|
|||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
|
||||
// printf("before mechmag, spin 1: [%g,%g,%g] \n",f[1][0],f[1][1],f[1][2]);
|
||||
|
||||
// compute and add magneto-mech. force
|
||||
// if (magpair_flag == 1) {
|
||||
// if (exch_flag) lockpairspin->compute_magnetomech(0,vflag);
|
||||
// }
|
||||
#define VSRSV_TEST
|
||||
#if defined VSRSV_TEST
|
||||
|
||||
// printf("after mechmag, spin 1: [%g,%g,%g] \n",f[1][0],f[1][1],f[1][2]);
|
||||
// update half v for all particles
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) dtfm = dtf / rmass[i];
|
||||
else dtfm = dtf / mass[type[i]];
|
||||
v[i][0] += dtfm * f[i][0];
|
||||
v[i][1] += dtfm * f[i][1];
|
||||
v[i][2] += dtfm * f[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
// update half v all particles
|
||||
// update s for all particles
|
||||
if (extra == SPIN) {
|
||||
if (mpi_flag == 1) {
|
||||
int nseci;
|
||||
// mpi seq. update
|
||||
for (int j = 0; j < nsectors; j++) {
|
||||
comm->forward_comm();
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
for (int j = nsectors-1; j >= 0; j--) {
|
||||
comm->forward_comm();
|
||||
for (int i = nlocal-1; i >= 0; i--) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
} else if (mpi_flag == 0) {
|
||||
// serial seq. update
|
||||
// advance quarter s for nlocal-1
|
||||
for (int i = 0; i < nlocal-1; i++){
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
}
|
||||
// advance half s for 1
|
||||
ComputeInteractionsSpin(nlocal-1);
|
||||
AdvanceSingleSpin(nlocal-1,dts,sp,fm);
|
||||
// advance quarter s for nlocal-1
|
||||
for (int i = nlocal-2; i >= 0; i--){
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
}
|
||||
} else error->all(FLERR,"Illegal fix integration/spin command");
|
||||
}
|
||||
|
||||
// update x for all particles
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
x[i][0] += 2.0 * dtv * v[i][0];
|
||||
x[i][1] += 2.0 * dtv * v[i][1];
|
||||
x[i][2] += 2.0 * dtv * v[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
// update half s for all particles
|
||||
if (extra == SPIN) {
|
||||
if (mpi_flag == 1) {
|
||||
int nseci;
|
||||
// mpi seq. update
|
||||
for (int j = 0; j < nsectors; j++) {
|
||||
comm->forward_comm();
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
for (int j = nsectors-1; j >= 0; j--) {
|
||||
comm->forward_comm();
|
||||
for (int i = nlocal-1; i >= 0; i--) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
} else if (mpi_flag == 0) {
|
||||
// serial seq. update
|
||||
// advance quarter s for nlocal-2 particles
|
||||
for (int i = nlocal-1; i >= 1; i--){
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
}
|
||||
// advance half s for 1
|
||||
ComputeInteractionsSpin(0);
|
||||
AdvanceSingleSpin(0,dts,sp,fm);
|
||||
// advance quarter s for nlocal-2 particles
|
||||
for (int i = 1; i < nlocal; i++){
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
}
|
||||
} else error->all(FLERR,"Illegal fix integration/spin command");
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
//#define VRSRV
|
||||
#if defined VRSRV
|
||||
// update half v for all particles
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) dtfm = dtf / rmass[i];
|
||||
|
|
@ -229,16 +357,17 @@ void FixIntegrationSpin::initial_integrate(int vflag)
|
|||
// update half x for all particles
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
x[i][0] += 0.5 * dtv * v[i][0];
|
||||
x[i][1] += 0.5 * dtv * v[i][1];
|
||||
x[i][2] += 0.5 * dtv * v[i][2];
|
||||
x[i][0] += dtv * v[i][0];
|
||||
x[i][1] += dtv * v[i][1];
|
||||
x[i][2] += dtv * v[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// update s for all particles
|
||||
if (extra == SPIN) {
|
||||
if (mpi_flag == 1) {
|
||||
int nseci;
|
||||
// mpi seq. update spins for all particles
|
||||
// mpi seq. update
|
||||
for (int j = 0; j < nsectors; j++) {
|
||||
comm->forward_comm();
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
|
|
@ -248,7 +377,7 @@ void FixIntegrationSpin::initial_integrate(int vflag)
|
|||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
for (int j = nsectors-1; j >= 0; j--) {
|
||||
|
|
@ -260,18 +389,18 @@ void FixIntegrationSpin::initial_integrate(int vflag)
|
|||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
} else if (mpi_flag == 0) {
|
||||
// serial seq. update spins for all particles
|
||||
// serial seq. update
|
||||
for (int i = 0; i < nlocal; i++){
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
for (int i = nlocal-1; i >= 0; i--){
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
} else error->all(FLERR,"Illegal fix integration/spin command");
|
||||
}
|
||||
|
|
@ -279,11 +408,75 @@ void FixIntegrationSpin::initial_integrate(int vflag)
|
|||
// update half x for all particles
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
x[i][0] += 0.5 * dtv * v[i][0];
|
||||
x[i][1] += 0.5 * dtv * v[i][1];
|
||||
x[i][2] += 0.5 * dtv * v[i][2];
|
||||
x[i][0] += dtv * v[i][0];
|
||||
x[i][1] += dtv * v[i][1];
|
||||
x[i][2] += dtv * v[i][2];
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
//#define RSVSR_TEST
|
||||
#if defined RSVSR_TEST
|
||||
|
||||
// update half x for all particles
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
x[i][0] += dtv * v[i][0];
|
||||
x[i][1] += dtv * v[i][1];
|
||||
x[i][2] += dtv * v[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
// update half s for all particles
|
||||
if (extra == SPIN) {
|
||||
if (mpi_flag == 1) {
|
||||
int nseci;
|
||||
// mpi seq. update
|
||||
for (int j = 0; j < nsectors; j++) {
|
||||
comm->forward_comm();
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
for (int j = nsectors-1; j >= 0; j--) {
|
||||
comm->forward_comm();
|
||||
for (int i = nlocal-1; i >= 0; i--) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
} else if (mpi_flag == 0) {
|
||||
// serial seq. update
|
||||
// advance quarter s for nlocal-2 particles
|
||||
for (int i = 0; i < nlocal-2; i++){
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
}
|
||||
// advance half s for nlocal-1
|
||||
ComputeInteractionsSpin(nlocal-1);
|
||||
AdvanceSingleSpin(nlocal-1,dts,sp,fm);
|
||||
// advance quarter s for nlocal-2 particles
|
||||
for (int i = nlocal-2; i >= 0; i--){
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
}
|
||||
} else error->all(FLERR,"Illegal fix integration/spin command");
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
|
@ -302,17 +495,17 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
|
|||
const int newton_pair = force->newton_pair;
|
||||
|
||||
// add test here
|
||||
if (magpair_flag) {
|
||||
inum = lockpairspin->list->inum;
|
||||
ilist = lockpairspin->list->ilist;
|
||||
numneigh = lockpairspin->list->numneigh;
|
||||
firstneigh = lockpairspin->list->firstneigh;
|
||||
if (exch_flag) {
|
||||
inum = lockpairspinexchange->list->inum;
|
||||
ilist = lockpairspinexchange->list->ilist;
|
||||
numneigh = lockpairspinexchange->list->numneigh;
|
||||
firstneigh = lockpairspinexchange->list->firstneigh;
|
||||
}
|
||||
|
||||
double xtmp,ytmp,ztmp;
|
||||
double rsq,rd,delx,dely,delz;
|
||||
double cut_ex_2, cut_dmi_2, cut_me_2;
|
||||
cut_ex_2 = cut_dmi_2 = cut_me_2 = 0.0;
|
||||
double rsq, rd;
|
||||
double delx, dely, delz;
|
||||
double temp_cut, cut_2, inorm;
|
||||
temp_cut = cut_2 = inorm = 0.0;
|
||||
|
||||
int eflag = 1;
|
||||
int vflag = 0;
|
||||
|
|
@ -341,43 +534,43 @@ void FixIntegrationSpin::ComputeInteractionsSpin(int ii)
|
|||
for (int jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
itype = type[ii];
|
||||
jtype = type[j];
|
||||
|
||||
spj[0] = sp[j][0];
|
||||
spj[1] = sp[j][1];
|
||||
spj[2] = sp[j][2];
|
||||
|
||||
delx = xi[0] - x[j][0];
|
||||
dely = xi[1] - x[j][1];
|
||||
delz = xi[2] - x[j][2];
|
||||
delx = x[j][0] - xi[0];
|
||||
dely = x[j][1] - xi[1];
|
||||
delz = x[j][2] - xi[2];
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
itype = type[ii];
|
||||
jtype = type[j];
|
||||
inorm = 1.0/sqrt(rsq);
|
||||
|
||||
if (magpair_flag) { // mag. pair inter.
|
||||
double temp_cut;
|
||||
if (exch_flag) { // exchange
|
||||
temp_cut = lockpairspin->cut_spin_exchange[itype][jtype];
|
||||
cut_ex_2 = temp_cut*temp_cut;
|
||||
if (rsq <= cut_ex_2) {
|
||||
lockpairspin->compute_exchange(i,j,rsq,fmi,fmj,spi,spj);
|
||||
}
|
||||
}
|
||||
if (dmi_flag) { // dmi
|
||||
temp_cut = lockpairspin->cut_spin_dmi[itype][jtype];
|
||||
cut_dmi_2 = temp_cut*temp_cut;
|
||||
if (rsq <= cut_dmi_2) {
|
||||
lockpairspin->compute_dmi(i,j,fmi,fmj,spi,spj);
|
||||
}
|
||||
}
|
||||
if (me_flag) { // me
|
||||
temp_cut = lockpairspin->cut_spin_me[itype][jtype];
|
||||
cut_me_2 = temp_cut*temp_cut;
|
||||
if (rsq <= cut_me_2) {
|
||||
lockpairspin->compute_me(i,j,fmi,fmj,spi,spj);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
rij[0] = delx*inorm;
|
||||
rij[1] = dely*inorm;
|
||||
rij[2] = delz*inorm;
|
||||
|
||||
temp_cut = 0.0;
|
||||
|
||||
if (exch_flag) { // exchange
|
||||
temp_cut = lockpairspinexchange->cut_spin_exchange[itype][jtype];
|
||||
cut_2 = temp_cut*temp_cut;
|
||||
if (rsq <= cut_2) {
|
||||
lockpairspinexchange->compute_exchange(i,j,rsq,fmi,fmj,spi,spj);
|
||||
}
|
||||
}
|
||||
|
||||
if (soc_flag) { // soc
|
||||
temp_cut = lockpairspinsocneel->cut_soc_neel[itype][jtype];
|
||||
cut_2 = temp_cut*temp_cut;
|
||||
if (rsq <= cut_2) {
|
||||
lockpairspinsocneel->compute_soc_neel(i,j,rsq,rij,fmi,fmj,spi,spj);
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
if (magforce_flag) { // mag. forces
|
||||
if (zeeman_flag) { // zeeman
|
||||
lockforcespin->compute_zeeman(i,fmi);
|
||||
|
|
@ -421,7 +614,7 @@ void FixIntegrationSpin::sectoring()
|
|||
const double rsy = subhi[1] - sublo[1];
|
||||
const double rsz = subhi[2] - sublo[2];
|
||||
|
||||
const double rv = lockpairspin->cut_spin_pair_global;
|
||||
const double rv = lockpairspinexchange->cut_spin_exchange_global;
|
||||
|
||||
double rax = rsx/rv;
|
||||
double ray = rsy/rv;
|
||||
|
|
@ -474,9 +667,9 @@ int FixIntegrationSpin::coords2sector(double *xi)
|
|||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixIntegrationSpin::AdvanceSingleSpin(int i, double dts, double **sp, double **fm)
|
||||
void FixIntegrationSpin::AdvanceSingleSpin(int i, double dtl, double **sp, double **fm)
|
||||
{
|
||||
double dtfm,msq,scale,fm2,fmsq,sp2,spsq,energy;
|
||||
double dtfm,msq,scale,fm2,fmsq,sp2,spsq,energy,dts2;
|
||||
double cp[3],g[3];
|
||||
|
||||
cp[0] = cp[1] = cp[2] = 0.0;
|
||||
|
|
@ -484,22 +677,23 @@ void FixIntegrationSpin::AdvanceSingleSpin(int i, double dts, double **sp, doubl
|
|||
fm2 = (fm[i][0]*fm[i][0])+(fm[i][1]*fm[i][1])+(fm[i][2]*fm[i][2]);
|
||||
fmsq = sqrt(fm2);
|
||||
energy = (sp[i][0]*fm[i][0])+(sp[i][1]*fm[i][1])+(sp[i][2]*fm[i][2]);
|
||||
|
||||
dts2 = dtl*dtl;
|
||||
|
||||
cp[0] = fm[i][1]*sp[i][2]-fm[i][2]*sp[i][1];
|
||||
cp[1] = fm[i][2]*sp[i][0]-fm[i][0]*sp[i][2];
|
||||
cp[2] = fm[i][0]*sp[i][1]-fm[i][1]*sp[i][0];
|
||||
|
||||
g[0] = sp[i][0]+cp[0]*dts;
|
||||
g[1] = sp[i][1]+cp[1]*dts;
|
||||
g[2] = sp[i][2]+cp[2]*dts;
|
||||
g[0] = sp[i][0]+cp[0]*dtl;
|
||||
g[1] = sp[i][1]+cp[1]*dtl;
|
||||
g[2] = sp[i][2]+cp[2]*dtl;
|
||||
|
||||
g[0] += (fm[i][0]*energy-0.5*sp[i][0]*fm2)*0.5*dts*dts;
|
||||
g[1] += (fm[i][1]*energy-0.5*sp[i][1]*fm2)*0.5*dts*dts;
|
||||
g[2] += (fm[i][2]*energy-0.5*sp[i][2]*fm2)*0.5*dts*dts;
|
||||
g[0] += (fm[i][0]*energy-0.5*sp[i][0]*fm2)*0.5*dts2;
|
||||
g[1] += (fm[i][1]*energy-0.5*sp[i][1]*fm2)*0.5*dts2;
|
||||
g[2] += (fm[i][2]*energy-0.5*sp[i][2]*fm2)*0.5*dts2;
|
||||
|
||||
g[0] /= (1+0.25*fm2*dts*dts);
|
||||
g[1] /= (1+0.25*fm2*dts*dts);
|
||||
g[2] /= (1+0.25*fm2*dts*dts);
|
||||
g[0] /= (1+0.25*fm2*dts2);
|
||||
g[1] /= (1+0.25*fm2*dts2);
|
||||
g[2] /= (1+0.25*fm2*dts2);
|
||||
|
||||
sp[i][0] = g[0];
|
||||
sp[i][1] = g[1];
|
||||
|
|
@ -523,6 +717,8 @@ void FixIntegrationSpin::final_integrate()
|
|||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double **sp = atom->sp;
|
||||
double **fm = atom->fm;
|
||||
double *rmass = atom->rmass;
|
||||
double *mass = atom->mass;
|
||||
int nlocal = atom->nlocal;
|
||||
|
|
@ -530,13 +726,11 @@ void FixIntegrationSpin::final_integrate()
|
|||
int *type = atom->type;
|
||||
int *mask = atom->mask;
|
||||
|
||||
// compute and add magneto-mech. force
|
||||
// if (magpair_flag == 1) {
|
||||
// if (exch_flag) lockpairspin->compute_magnetomech(0,0);
|
||||
// }
|
||||
#define VSRSV_TEST
|
||||
#if defined VSRSV_TEST
|
||||
|
||||
// update half v for all particles
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
for (int i = nlocal-1; i >= 0; i--) {
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) dtfm = dtf / rmass[i];
|
||||
else dtfm = dtf / mass[type[i]];
|
||||
|
|
@ -545,4 +739,93 @@ void FixIntegrationSpin::final_integrate()
|
|||
v[i][2] += dtfm * f[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
//#define VRSRV
|
||||
#if defined VRSRV
|
||||
// update half v for all particles
|
||||
for (int i = nlocal-1; i >= 0; i--) {
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) dtfm = dtf / rmass[i];
|
||||
else dtfm = dtf / mass[type[i]];
|
||||
v[i][0] += dtfm * f[i][0];
|
||||
v[i][1] += dtfm * f[i][1];
|
||||
v[i][2] += dtfm * f[i][2];
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
//#define RSVSR_TEST
|
||||
#if defined RSVSR_TEST
|
||||
|
||||
// update v for all particles
|
||||
for (int i = nlocal-1; i >= 0; i--) {
|
||||
if (mask[i] & groupbit) {
|
||||
if (rmass) dtfm = dtf / rmass[i];
|
||||
else dtfm = dtf / mass[type[i]];
|
||||
v[i][0] += 2.0 * dtfm * f[i][0];
|
||||
v[i][1] += 2.0 * dtfm * f[i][1];
|
||||
v[i][2] += 2.0 * dtfm * f[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
// update half s for all particles
|
||||
if (extra == SPIN) {
|
||||
if (mpi_flag == 1) {
|
||||
int nseci;
|
||||
// mpi seq. update
|
||||
for (int j = 0; j < nsectors; j++) {
|
||||
comm->forward_comm();
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
for (int j = nsectors-1; j >= 0; j--) {
|
||||
comm->forward_comm();
|
||||
for (int i = nlocal-1; i >= 0; i--) {
|
||||
xi[0] = x[i][0];
|
||||
xi[1] = x[i][1];
|
||||
xi[2] = x[i][2];
|
||||
nseci = coords2sector(xi);
|
||||
if (j != nseci) continue;
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,dts,sp,fm);
|
||||
}
|
||||
}
|
||||
} else if (mpi_flag == 0) {
|
||||
// serial seq. update
|
||||
// advance quarter s for nlocal-2 particles
|
||||
for (int i = nlocal-1; i >= 1; i--){
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
}
|
||||
// advance half s for nlocal-1
|
||||
ComputeInteractionsSpin(0);
|
||||
AdvanceSingleSpin(0,dts,sp,fm);
|
||||
// advance quarter s for nlocal-2 particles
|
||||
for (int i = 1; i < nlocal-1; i++){
|
||||
ComputeInteractionsSpin(i);
|
||||
AdvanceSingleSpin(i,0.5*dts,sp,fm);
|
||||
}
|
||||
} else error->all(FLERR,"Illegal fix integration/spin command");
|
||||
}
|
||||
|
||||
// update half x for all particles
|
||||
for (int i = nlocal-1; i >= 0; i--) {
|
||||
if (mask[i] & groupbit) {
|
||||
x[i][0] += dtv * v[i][0];
|
||||
x[i][1] += dtv * v[i][1];
|
||||
x[i][2] += dtv * v[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -50,20 +50,23 @@ class FixIntegrationSpin : public Fix {
|
|||
|
||||
// mag. interaction flags
|
||||
int magpair_flag;
|
||||
int exch_flag, dmi_flag, me_flag;
|
||||
int soc_flag;
|
||||
int exch_flag;
|
||||
int magforce_flag;
|
||||
int zeeman_flag, aniso_flag;
|
||||
int maglangevin_flag;
|
||||
int tdamp_flag, temp_flag;
|
||||
|
||||
// pointers to mag. interaction classes
|
||||
// pointers to interaction classes
|
||||
class PairHybrid *lockhybrid;
|
||||
class PairSpin *lockpairspin;
|
||||
class PairSpinExchange *lockpairspinexchange;
|
||||
class PairSpinSocNeel *lockpairspinsocneel;
|
||||
class FixForceSpin *lockforcespin;
|
||||
class FixLangevinSpin *locklangevinspin;
|
||||
|
||||
// temporary variables
|
||||
double *xi;
|
||||
double *xi, *rij;
|
||||
double *spi, *spj;
|
||||
double *fmi, *fmj;
|
||||
|
||||
|
|
|
|||
|
|
@ -129,16 +129,17 @@ void FixLangevinSpin::init()
|
|||
}
|
||||
if (flag_force >= flag_lang) error->all(FLERR,"Fix langevin/spin should come after all other spin fixes");
|
||||
|
||||
memory->create(spi,3,"pair:spi");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(spi,3,"langevin:spi");
|
||||
memory->create(fmi,3,"langevin:fmi");
|
||||
|
||||
gil_factor = 1.0/(1.0+(alpha_t)*(alpha_t));
|
||||
dts = update->dt;
|
||||
Gil_factor = 1.0/(1.0+(alpha_t)*(alpha_t));
|
||||
|
||||
double hbar = force->hplanck/MY_2PI; //eV/(rad.THz)
|
||||
double kb = force->boltz;
|
||||
D = (MY_2PI*Gil_factor*kb*temp)/hbar/dts;
|
||||
sigma = sqrt(D);
|
||||
double hbar = force->hplanck/MY_2PI; // eV/(rad.THz)
|
||||
double kb = force->boltz; // eV/K
|
||||
D = (MY_2PI*alpha_t*gil_factor*kb*temp);
|
||||
D /= (hbar*dts);
|
||||
sigma = sqrt(2.0*D);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -155,7 +156,7 @@ void FixLangevinSpin::setup(int vflag)
|
|||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
/*
|
||||
void FixLangevinSpin::post_force(int vflag)
|
||||
{
|
||||
double **sp = atom->sp;
|
||||
|
|
@ -163,11 +164,6 @@ void FixLangevinSpin::post_force(int vflag)
|
|||
int *mask = atom->mask;
|
||||
const int nlocal = atom->nlocal;
|
||||
|
||||
double sx, sy, sz;
|
||||
double fmx, fmy, fmz;
|
||||
double cpx, cpy, cpz;
|
||||
double rx, ry, rz;
|
||||
|
||||
// add the damping to the effective field of each spin
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
|
|
@ -194,6 +190,7 @@ void FixLangevinSpin::post_force(int vflag)
|
|||
}
|
||||
|
||||
}
|
||||
*/
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void FixLangevinSpin::add_tdamping(double *spi, double *fmi)
|
||||
|
|
@ -213,25 +210,27 @@ void FixLangevinSpin::add_temperature(double *fmi)
|
|||
{
|
||||
//#define GAUSSIAN_R
|
||||
#if defined GAUSSIAN_R
|
||||
// drawing gausian random dist
|
||||
// drawing gaussian random dist
|
||||
double rx = sigma*random->gaussian();
|
||||
double ry = sigma*random->gaussian();
|
||||
double rz = sigma*random->gaussian();
|
||||
#else
|
||||
double rx = sigma*(random->uniform() - 0.5);
|
||||
double ry = sigma*(random->uniform() - 0.5);
|
||||
double rz = sigma*(random->uniform() - 0.5);
|
||||
double rx = sigma*(-1.0+2.0*random->uniform());
|
||||
double ry = sigma*(-1.0+2.0*random->uniform());
|
||||
double rz = sigma*(-1.0+2.0*random->uniform());
|
||||
#endif
|
||||
|
||||
// printf("test Gaussian vals: %g \n",rx);
|
||||
|
||||
// adding the random field
|
||||
fmi[0] += rx;
|
||||
fmi[1] += ry;
|
||||
fmi[2] += rz;
|
||||
|
||||
// adding Gilbert's prefactor
|
||||
fmi[0] *= Gil_factor;
|
||||
fmi[1] *= Gil_factor;
|
||||
fmi[2] *= Gil_factor;
|
||||
fmi[0] *= gil_factor;
|
||||
fmi[1] *= gil_factor;
|
||||
fmi[2] *= gil_factor;
|
||||
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -31,7 +31,7 @@ class FixLangevinSpin : public Fix {
|
|||
int setmask();
|
||||
void init();
|
||||
void setup(int);
|
||||
virtual void post_force(int);
|
||||
// virtual void post_force(int);
|
||||
void post_force_respa(int, int, int);
|
||||
// add transverse damping and temperature
|
||||
void add_tdamping(double *, double *);
|
||||
|
|
@ -41,11 +41,13 @@ class FixLangevinSpin : public Fix {
|
|||
|
||||
protected:
|
||||
double *spi, *fmi;
|
||||
// transverse and longitudinal damping coeff.
|
||||
double alpha_t, alpha_l;
|
||||
// timestep, temperature, noise intensity
|
||||
double dts,temp,D,sigma;
|
||||
double Gil_factor;
|
||||
double alpha_t; // mag. transverse damping coeff.
|
||||
double alpha_l; // mag. longitudinal damping coeff.
|
||||
double dts; // timestep
|
||||
double temp; // spin bath temperature
|
||||
double D,sigma; // bath intensity var.
|
||||
|
||||
double gil_factor; // Gilbert's prefactor
|
||||
|
||||
char *id_temp;
|
||||
class Compute *temperature;
|
||||
|
|
|
|||
|
|
@ -0,0 +1,484 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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)
|
||||
Aidan Thompson (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "pair_hybrid.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "pair_spin_exchange.h"
|
||||
#include "update.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpinExchange::PairSpinExchange(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
hbar = force->hplanck/MY_2PI;
|
||||
|
||||
newton_pair_spin = 0; // no newton pair for now
|
||||
// newton_pair = 0;
|
||||
|
||||
single_enable = 0;
|
||||
exch_flag = 0;
|
||||
exch_mech_flag = 0;
|
||||
|
||||
no_virial_fdotr_compute = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpinExchange::~PairSpinExchange()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
|
||||
memory->destroy(cut_spin_exchange);
|
||||
memory->destroy(J1_mag);
|
||||
memory->destroy(J1_mech);
|
||||
memory->destroy(J2);
|
||||
memory->destroy(J3);
|
||||
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fi);
|
||||
memory->destroy(fj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
memory->destroy(rij);
|
||||
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinExchange::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double evdwl,ecoul;
|
||||
double xtmp,ytmp,ztmp;
|
||||
double fix,fiy,fiz,fjx,fjy,fjz;
|
||||
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
|
||||
double cut_ex_2,cut_spin_exchange_global2;
|
||||
double rsq,rd;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
cut_spin_exchange_global2 = cut_spin_exchange_global*cut_spin_exchange_global;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double **fm = atom->fm;
|
||||
double *mumag = atom->mumag;
|
||||
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;
|
||||
|
||||
// pair spin computations
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
spi[0] = sp[i][0];
|
||||
spi[1] = sp[i][1];
|
||||
spi[2] = sp[i][2];
|
||||
|
||||
// loop on neighbors
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
spj[0] = sp[j][0];
|
||||
spj[1] = sp[j][1];
|
||||
spj[2] = sp[j][2];
|
||||
|
||||
evdwl = 0.0;
|
||||
|
||||
fi[0] = fi[1] = fi[2] = 0.0;
|
||||
fj[0] = fj[1] = fj[2] = 0.0;
|
||||
fmi[0] = fmi[1] = fmi[2] = 0.0;
|
||||
fmj[0] = fmj[1] = fmj[2] = 0.0;
|
||||
rij[0] = rij[1] = rij[2] = 0.0;
|
||||
|
||||
rij[0] = x[j][0] - xtmp;
|
||||
rij[1] = x[j][1] - ytmp;
|
||||
rij[2] = x[j][2] - ztmp;
|
||||
|
||||
// square of inter-atomic distance
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
double inorm = 1.0/sqrt(rsq);
|
||||
rij[0] *= inorm;
|
||||
rij[1] *= inorm;
|
||||
rij[2] *= inorm;
|
||||
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
// exchange interaction
|
||||
if (exch_flag) {
|
||||
cut_ex_2 = cut_spin_exchange[itype][jtype]*cut_spin_exchange[itype][jtype];
|
||||
if (rsq <= cut_ex_2) {
|
||||
compute_exchange(i,j,rsq,fmi,fmj,spi,spj);
|
||||
compute_exchange_mech(i,j,rsq,rij,fi,fj,spi,spj);
|
||||
}
|
||||
}
|
||||
|
||||
f[i][0] += fi[0];
|
||||
f[i][1] += fi[1];
|
||||
f[i][2] += fi[2];
|
||||
fm[i][0] += fmi[0];
|
||||
fm[i][1] += fmi[1];
|
||||
fm[i][2] += fmi[2];
|
||||
|
||||
// if (newton_pair || j < nlocal) {
|
||||
if (newton_pair_spin) {
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
f[j][2] += fj[2];
|
||||
fm[j][0] += fmj[0];
|
||||
fm[j][1] += fmj[1];
|
||||
fm[j][2] += fmj[2];
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq <= cut_ex_2) {
|
||||
evdwl -= spi[0]*fmi[0];
|
||||
evdwl -= spi[1]*fmi[1];
|
||||
evdwl -= 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]);
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinExchange::compute_exchange(int i, int j, double rsq, double *fmi, double *fmj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
double Jex, ra;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
ra = rsq/J3[itype][jtype]/J3[itype][jtype];
|
||||
Jex = 4.0*J1_mag[itype][jtype]*ra;
|
||||
Jex *= (1.0-J2[itype][jtype]*ra);
|
||||
Jex *= exp(-ra);
|
||||
|
||||
fmi[0] += 0.5*Jex*spj[0];
|
||||
fmi[1] += 0.5*Jex*spj[1];
|
||||
fmi[2] += 0.5*Jex*spj[2];
|
||||
|
||||
fmj[0] += 0.5*Jex*spi[0];
|
||||
fmj[1] += 0.5*Jex*spi[1];
|
||||
fmj[2] += 0.5*Jex*spi[2];
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double *rij, double *fi, double *fj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
double Jex, Jex_mech, ra, rr, iJ3;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
Jex = J1_mech[itype][jtype];
|
||||
iJ3 = 1.0/(J3[itype][jtype]*J3[itype][jtype]);
|
||||
|
||||
ra = rsq*iJ3;
|
||||
rr = sqrt(rsq)*iJ3;
|
||||
|
||||
Jex_mech = 1.0-ra-J2[itype][jtype]*ra*(2.0-ra);
|
||||
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];
|
||||
|
||||
fj[0] -= Jex_mech*rij[0];
|
||||
fj[1] -= Jex_mech*rij[1];
|
||||
fj[2] -= Jex_mech*rij[2];
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinExchange::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_exchange,n+1,n+1,"pair:cut_spin_exchange");
|
||||
memory->create(J1_mag,n+1,n+1,"pair:J1_mag");
|
||||
memory->create(J1_mech,n+1,n+1,"pair:J1_mech");
|
||||
memory->create(J2,n+1,n+1,"pair:J2");
|
||||
memory->create(J3,n+1,n+1,"pair:J3");
|
||||
|
||||
memory->create(spi,3,"pair:spi");
|
||||
memory->create(spj,3,"pair:spj");
|
||||
memory->create(fi,3,"pair:fi");
|
||||
memory->create(fj,3,"pair:fj");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(fmj,3,"pair:fmj");
|
||||
memory->create(rij,3,"pair:rij");
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinExchange::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 1 || narg > 2)
|
||||
error->all(FLERR,"Incorrect number of args in pair_style pair/spin command");
|
||||
|
||||
if (strcmp(update->unit_style,"metal") != 0)
|
||||
error->all(FLERR,"Spin simulations require metal unit style");
|
||||
|
||||
cut_spin_exchange_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_exchange[i][j] = cut_spin_exchange_global;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type spin pairs (only one for now)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinExchange::coeff(int narg, char **arg)
|
||||
{
|
||||
const double hbar = force->hplanck/MY_2PI;
|
||||
|
||||
if (!allocated) allocate();
|
||||
|
||||
// set exch_mech_flag to 1 if magneto-mech simulation
|
||||
if (strstr(force->pair_style,"pair/spin")) {
|
||||
exch_mech_flag = 0;
|
||||
} else if (strstr(force->pair_style,"hybrid/overlay")) {
|
||||
exch_mech_flag = 1;
|
||||
} else error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
|
||||
if (strcmp(arg[2],"exchange")==0){
|
||||
if (narg != 7) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
exch_flag = 1;
|
||||
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
|
||||
const double rij = force->numeric(FLERR,arg[3]);
|
||||
const double j1 = (force->numeric(FLERR,arg[4]));
|
||||
const double j2 = force->numeric(FLERR,arg[5]);
|
||||
const double j3 = force->numeric(FLERR,arg[6]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
cut_spin_exchange[i][j] = rij;
|
||||
J1_mag[i][j] = j1/hbar;
|
||||
if (exch_mech_flag) {
|
||||
J1_mech[i][j] = j1;
|
||||
} else {
|
||||
J1_mech[i][j] = 0.0;
|
||||
}
|
||||
J2[i][j] = j2;
|
||||
J3[i][j] = j3;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
} else error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinExchange::init_style()
|
||||
{
|
||||
if (!atom->sp_flag || !atom->mumag_flag)
|
||||
error->all(FLERR,"Pair spin requires atom attributes sp, mumag");
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
// check this half/full request
|
||||
#define FULLNEI
|
||||
#if defined FULLNEI
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairSpinExchange::init_one(int i, int j)
|
||||
{
|
||||
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
|
||||
return cut_spin_exchange_global;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinExchange::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]) {
|
||||
if (exch_flag){
|
||||
fwrite(&J1_mag[i][j],sizeof(double),1,fp);
|
||||
fwrite(&J1_mech[i][j],sizeof(double),1,fp);
|
||||
fwrite(&J2[i][j],sizeof(double),1,fp);
|
||||
fwrite(&J3[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_spin_exchange[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinExchange::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(&J1_mag[i][j],sizeof(double),1,fp);
|
||||
fread(&J1_mech[i][j],sizeof(double),1,fp);
|
||||
fread(&J2[i][j],sizeof(double),1,fp);
|
||||
fread(&J2[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_spin_exchange[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&J1_mag[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&J1_mech[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&J2[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&J3[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_spin_exchange[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinExchange::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_spin_exchange_global,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinExchange::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&cut_spin_exchange_global,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&cut_spin_exchange_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
|
@ -0,0 +1,89 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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(pair/spin/exchange,PairSpinExchange)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_SPIN_EXCHANGE_H
|
||||
#define LMP_PAIR_SPIN_EXCHANGE_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairSpinExchange : public Pair {
|
||||
public:
|
||||
PairSpinExchange(class LAMMPS *);
|
||||
virtual ~PairSpinExchange();
|
||||
virtual void compute(int, int);
|
||||
void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
void compute_exchange(int, int, double, double *, double *,double *, double *);
|
||||
void compute_exchange_mech(int, int, double, double *, double *, double *,double *, double *);
|
||||
|
||||
int exch_flag; // mag. exchange flag
|
||||
int exch_mech_flag; // mech. exchange flags
|
||||
|
||||
double cut_spin_exchange_global; // global exchange cutoff
|
||||
double **cut_spin_exchange; // cutoff distance exchange
|
||||
|
||||
protected:
|
||||
int newton_pair_spin;
|
||||
double hbar;
|
||||
|
||||
double **J1_mag, **J1_mech; // exchange coeffs Jij
|
||||
double **J2, **J3; // J1 in eV, J2 adim, J3 in Ang
|
||||
|
||||
double *spi, *spj; // temp. spin vals. in compute
|
||||
double *fi, *fj; // temp. mech. forces in compute
|
||||
double *fmi, *fmj; // temp. mag. forces in compute
|
||||
double *rij; // norm. inter atomic vectors
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Incorrect args in pair_spin command
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Spin simulations require metal unit style
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair spin requires atom attributes sp, mumag
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,460 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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)
|
||||
Aidan Thompson (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "pair_hybrid.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "pair_spin_me.h"
|
||||
#include "update.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpinMe::PairSpinMe(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
hbar = force->hplanck/MY_2PI;
|
||||
|
||||
newton_pair_spin = 0; // no newton pair for now
|
||||
// newton_pair = 0;
|
||||
|
||||
single_enable = 0;
|
||||
me_flag = 0;
|
||||
|
||||
no_virial_fdotr_compute = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpinMe::~PairSpinMe()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
|
||||
memory->destroy(cut_spin_me);
|
||||
memory->destroy(ME);
|
||||
memory->destroy(v_mex);
|
||||
memory->destroy(v_mey);
|
||||
memory->destroy(v_mez);
|
||||
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fi);
|
||||
memory->destroy(fj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
memory->destroy(rij);
|
||||
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinMe::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double evdwl,ecoul;
|
||||
double xi,yi,zi;
|
||||
double fix,fiy,fiz,fjx,fjy,fjz;
|
||||
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
|
||||
double cut_me_2,cut_spin_me_global2;
|
||||
double rsq,rd;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
cut_spin_me_global2 = cut_spin_me_global*cut_spin_me_global;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double **fm = atom->fm;
|
||||
double *mumag = atom->mumag;
|
||||
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;
|
||||
|
||||
// magneto-electric computation
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xi = x[i][0];
|
||||
yi = x[i][1];
|
||||
zi = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
spi[0] = sp[i][0];
|
||||
spi[1] = sp[i][1];
|
||||
spi[2] = sp[i][2];
|
||||
|
||||
// loop on neighbors
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
spj[0] = sp[j][0];
|
||||
spj[1] = sp[j][1];
|
||||
spj[2] = sp[j][2];
|
||||
|
||||
evdwl = 0.0;
|
||||
|
||||
fi[0] = fi[1] = fi[2] = 0.0;
|
||||
fj[0] = fj[1] = fj[2] = 0.0;
|
||||
fmi[0] = fmi[1] = fmi[2] = 0.0;
|
||||
fmj[0] = fmj[1] = fmj[2] = 0.0;
|
||||
rij[0] = rij[1] = rij[2] = 0.0;
|
||||
|
||||
rij[0] = x[j][0] - xi;
|
||||
rij[1] = x[j][1] - yi;
|
||||
rij[2] = x[j][2] - zi;
|
||||
|
||||
// square of inter-atomic distance
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
double inorm = 1.0/sqrt(rsq);
|
||||
rij[0] *= inorm;
|
||||
rij[1] *= inorm;
|
||||
rij[2] *= inorm;
|
||||
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
// me interaction
|
||||
if (me_flag){
|
||||
cut_me_2 = cut_spin_me[itype][jtype]*cut_spin_me[itype][jtype];
|
||||
if (rsq <= cut_me_2){
|
||||
compute_me(i,j,fmi,fmj,spi,spj);
|
||||
}
|
||||
}
|
||||
|
||||
f[i][0] += fi[0];
|
||||
f[i][1] += fi[1];
|
||||
f[i][2] += fi[2];
|
||||
fm[i][0] += fmi[0];
|
||||
fm[i][1] += fmi[1];
|
||||
fm[i][2] += fmi[2];
|
||||
|
||||
// if (newton_pair || j < nlocal) {
|
||||
if (newton_pair_spin) {
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
f[j][2] += fj[2];
|
||||
fm[j][0] += fmj[0];
|
||||
fm[j][1] += fmj[1];
|
||||
fm[j][2] += fmj[2];
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq <= cut_me_2) {
|
||||
evdwl -= spi[0]*fmi[0];
|
||||
evdwl -= spi[1]*fmi[1];
|
||||
evdwl -= 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]);
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinMe::compute_me(int i, int j, double *fmi, double *fmj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
double **sp = atom->sp;
|
||||
double **x = atom->x;
|
||||
double meix,meiy,meiz;
|
||||
double rx, ry, rz, inorm;
|
||||
|
||||
rx = x[j][0] - x[i][0];
|
||||
ry = x[j][1] - x[i][1];
|
||||
rz = x[j][2] - x[i][2];
|
||||
inorm = 1.0/sqrt(rx*rx+ry*ry+rz*rz);
|
||||
rx *= inorm;
|
||||
ry *= inorm;
|
||||
rz *= inorm;
|
||||
|
||||
meix = v_mey[itype][jtype]*rz - v_mez[itype][jtype]*ry;
|
||||
meiy = v_mez[itype][jtype]*rx - v_mex[itype][jtype]*rz;
|
||||
meiz = v_mex[itype][jtype]*ry - v_mey[itype][jtype]*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;
|
||||
|
||||
fmj[0] -= spi[1]*meiz - spi[2]*meiy;
|
||||
fmj[1] -= spi[2]*meix - spi[0]*meiz;
|
||||
fmj[2] -= spi[0]*meiy - spi[1]*meix;
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinMe::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_me,n+1,n+1,"pair:cut_spin_me");
|
||||
memory->create(ME,n+1,n+1,"pair:ME");
|
||||
memory->create(v_mex,n+1,n+1,"pair:ME_vector_x");
|
||||
memory->create(v_mey,n+1,n+1,"pair:ME_vector_y");
|
||||
memory->create(v_mez,n+1,n+1,"pair:ME_vector_z");
|
||||
|
||||
memory->create(spi,3,"pair:spi");
|
||||
memory->create(spj,3,"pair:spj");
|
||||
memory->create(fi,3,"pair:fi");
|
||||
memory->create(fj,3,"pair:fj");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(fmj,3,"pair:fmj");
|
||||
memory->create(rij,3,"pair:rij");
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinMe::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 1 || narg > 2)
|
||||
error->all(FLERR,"Incorrect number of args in pair_style pair/spin command");
|
||||
|
||||
if (strcmp(update->unit_style,"metal") != 0)
|
||||
error->all(FLERR,"Spin simulations require metal unit style");
|
||||
|
||||
cut_spin_me_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_me[i][j] = cut_spin_me_global;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type spin pairs (only one for now)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinMe::coeff(int narg, char **arg)
|
||||
{
|
||||
const double hbar = force->hplanck/MY_2PI;
|
||||
|
||||
if (!allocated) allocate();
|
||||
|
||||
if (strcmp(arg[2],"me")==0) {
|
||||
if (narg != 8) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
me_flag = 1;
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
|
||||
const double rij = force->numeric(FLERR,arg[3]);
|
||||
const double me = (force->numeric(FLERR,arg[4]))/hbar;
|
||||
double mex = force->numeric(FLERR,arg[5]);
|
||||
double mey = force->numeric(FLERR,arg[6]);
|
||||
double mez = force->numeric(FLERR,arg[7]);
|
||||
|
||||
double inorm = 1.0/(mex*mex+mey*mey+mez*mez);
|
||||
mex *= inorm;
|
||||
mey *= inorm;
|
||||
mez *= inorm;
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
cut_spin_me[i][j] = rij;
|
||||
ME[i][j] = me;
|
||||
v_mex[i][j] = mex;
|
||||
v_mey[i][j] = mey;
|
||||
v_mez[i][j] = mez;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
} else error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinMe::init_style()
|
||||
{
|
||||
if (!atom->sp_flag || !atom->mumag_flag)
|
||||
error->all(FLERR,"Pair spin requires atom attributes sp, mumag");
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
// check this half/full request
|
||||
#define FULLNEI
|
||||
#if defined FULLNEI
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairSpinMe::init_one(int i, int j)
|
||||
{
|
||||
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
|
||||
return cut_spin_me_global;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinMe::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]) {
|
||||
if (me_flag) {
|
||||
fwrite(&ME[i][j],sizeof(double),1,fp);
|
||||
fwrite(&v_mex[i][j],sizeof(double),1,fp);
|
||||
fwrite(&v_mey[i][j],sizeof(double),1,fp);
|
||||
fwrite(&v_mez[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_spin_me[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinMe::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(&ME[i][j],sizeof(double),1,fp);
|
||||
fread(&v_mex[i][j],sizeof(double),1,fp);
|
||||
fread(&v_mey[i][j],sizeof(double),1,fp);
|
||||
fread(&v_mez[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_spin_me[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&ME[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&v_mex[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&v_mey[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&v_mez[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_spin_me[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinMe::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_spin_me_global,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinMe::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&cut_spin_me_global,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&cut_spin_me_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
|
@ -0,0 +1,87 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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(pair/spin/me,PairSpinMe)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_SPIN_ME_H
|
||||
#define LMP_PAIR_SPIN_ME_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairSpinMe : public Pair {
|
||||
public:
|
||||
PairSpinMe(class LAMMPS *);
|
||||
virtual ~PairSpinMe();
|
||||
virtual void compute(int, int);
|
||||
void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
void compute_me(int, int, double *, double *, double *, double *);
|
||||
|
||||
int me_flag; // me flag
|
||||
|
||||
double cut_spin_me_global; // global me cutoff
|
||||
double **cut_spin_me; // me cutoff distance
|
||||
|
||||
protected:
|
||||
int newton_pair_spin;
|
||||
double hbar;
|
||||
|
||||
double **ME; // me coeff in eV
|
||||
double **v_mex, **v_mey, **v_mez;// me direction
|
||||
|
||||
double *spi, *spj; // temp. spin vals. in compute
|
||||
double *fi, *fj; // temp. mech. forces in compute
|
||||
double *fmi, *fmj; // temp. mag. forces in compute
|
||||
double *rij; // norm. inter atomic vectors
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Incorrect args in pair_spin command
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Spin simulations require metal unit style
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair spin requires atom attributes sp, mumag
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,444 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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)
|
||||
Aidan Thompson (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "pair_hybrid.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "pair_spin_soc_dmi.h"
|
||||
#include "update.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpinSocDmi::PairSpinSocDmi(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
hbar = force->hplanck/MY_2PI;
|
||||
|
||||
newton_pair_spin = 0; // no newton pair for now
|
||||
// newton_pair = 0;
|
||||
|
||||
single_enable = 0;
|
||||
dmi_flag = 0;
|
||||
|
||||
no_virial_fdotr_compute = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpinSocDmi::~PairSpinSocDmi()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
|
||||
memory->destroy(cut_spin_dmi);
|
||||
memory->destroy(DM);
|
||||
memory->destroy(v_dmx);
|
||||
memory->destroy(v_dmy);
|
||||
memory->destroy(v_dmz);
|
||||
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fi);
|
||||
memory->destroy(fj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
memory->destroy(rij);
|
||||
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocDmi::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double evdwl,ecoul;
|
||||
double xi,yi,zi;
|
||||
double fix,fiy,fiz,fjx,fjy,fjz;
|
||||
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
|
||||
double cut_dmi_2;
|
||||
double rsq,rd;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
cut_dmi_2 = cut_spin_dmi_global*cut_spin_dmi_global;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double **fm = atom->fm;
|
||||
double *mumag = atom->mumag;
|
||||
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;
|
||||
|
||||
// dmi computation
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xi = x[i][0];
|
||||
yi = x[i][1];
|
||||
zi = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
spi[0] = sp[i][0];
|
||||
spi[1] = sp[i][1];
|
||||
spi[2] = sp[i][2];
|
||||
|
||||
// loop on neighbors
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
spj[0] = sp[j][0];
|
||||
spj[1] = sp[j][1];
|
||||
spj[2] = sp[j][2];
|
||||
|
||||
evdwl = 0.0;
|
||||
|
||||
fi[0] = fi[1] = fi[2] = 0.0;
|
||||
fj[0] = fj[1] = fj[2] = 0.0;
|
||||
fmi[0] = fmi[1] = fmi[2] = 0.0;
|
||||
fmj[0] = fmj[1] = fmj[2] = 0.0;
|
||||
rij[0] = rij[1] = rij[2] = 0.0;
|
||||
|
||||
rij[0] = x[j][0] - xi;
|
||||
rij[1] = x[j][1] - yi;
|
||||
rij[2] = x[j][2] - zi;
|
||||
|
||||
// square of inter-atomic distance
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
double inorm = 1.0/sqrt(rsq);
|
||||
rij[0] *= inorm;
|
||||
rij[1] *= inorm;
|
||||
rij[2] *= inorm;
|
||||
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
// dm interaction
|
||||
if (dmi_flag){
|
||||
cut_dmi_2 = cut_spin_dmi[itype][jtype]*cut_spin_dmi[itype][jtype];
|
||||
if (rsq <= cut_dmi_2){
|
||||
compute_dmi(i,j,fmi,fmj,spi,spj);
|
||||
}
|
||||
}
|
||||
|
||||
f[i][0] += fi[0];
|
||||
f[i][1] += fi[1];
|
||||
f[i][2] += fi[2];
|
||||
fm[i][0] += fmi[0];
|
||||
fm[i][1] += fmi[1];
|
||||
fm[i][2] += fmi[2];
|
||||
|
||||
// if (newton_pair || j < nlocal) {
|
||||
if (newton_pair_spin) {
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
f[j][2] += fj[2];
|
||||
fm[j][0] += fmj[0];
|
||||
fm[j][1] += fmj[1];
|
||||
fm[j][2] += fmj[2];
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq <= cut_dmi_2) {
|
||||
evdwl -= spi[0]*fmi[0];
|
||||
evdwl -= spi[1]*fmi[1];
|
||||
evdwl -= 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]);
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinSocDmi::compute_dmi(int i, int j, double *fmi, double *fmj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
double dmix,dmiy,dmiz;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
dmix = DM[itype][jtype]*v_dmx[itype][jtype];
|
||||
dmiy = DM[itype][jtype]*v_dmy[itype][jtype];
|
||||
dmiz = DM[itype][jtype]*v_dmz[itype][jtype];
|
||||
|
||||
fmi[0] += spj[1]*dmiz-spj[2]*dmiy;
|
||||
fmi[1] += spj[2]*dmix-spj[0]*dmiz;
|
||||
fmi[2] += spj[0]*dmiy-spj[1]*dmix;
|
||||
|
||||
fmj[0] -= spi[1]*dmiz-spi[2]*dmiy;
|
||||
fmj[1] -= spi[2]*dmix-spi[0]*dmiz;
|
||||
fmj[2] -= spi[0]*dmiy-spi[1]*dmix;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocDmi::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_dmi,n+1,n+1,"pair:cut_spin_dmi");
|
||||
memory->create(DM,n+1,n+1,"pair:DM");
|
||||
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(spi,3,"pair:spi");
|
||||
memory->create(spj,3,"pair:spj");
|
||||
memory->create(fi,3,"pair:fi");
|
||||
memory->create(fj,3,"pair:fj");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(fmj,3,"pair:fmj");
|
||||
memory->create(rij,3,"pair:rij");
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocDmi::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 1 || narg > 2)
|
||||
error->all(FLERR,"Incorrect number of args in pair/spin/dmi command");
|
||||
|
||||
if (strcmp(update->unit_style,"metal") != 0)
|
||||
error->all(FLERR,"Spin simulations require metal unit style");
|
||||
|
||||
cut_spin_dmi_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_dmi[i][j] = cut_spin_dmi_global;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type spin pairs (only one for now)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocDmi::coeff(int narg, char **arg)
|
||||
{
|
||||
const double hbar = force->hplanck/MY_2PI;
|
||||
|
||||
if (!allocated) allocate();
|
||||
|
||||
if (strcmp(arg[2],"dmi")==0) {
|
||||
if (narg != 8) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
dmi_flag = 1;
|
||||
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
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;
|
||||
double dmx = force->numeric(FLERR,arg[5]);
|
||||
double dmy = force->numeric(FLERR,arg[6]);
|
||||
double dmz = force->numeric(FLERR,arg[7]);
|
||||
|
||||
double inorm = 1.0/(dmx*dmx+dmy*dmy+dmz*dmz);
|
||||
dmx *= inorm;
|
||||
dmy *= inorm;
|
||||
dmz *= inorm;
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
cut_spin_dmi[i][j] = rij;
|
||||
DM[i][j] = dm;
|
||||
v_dmx[i][j] = dmx;
|
||||
v_dmy[i][j] = dmy;
|
||||
v_dmz[i][j] = dmz;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
} else error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocDmi::init_style()
|
||||
{
|
||||
if (!atom->sp_flag || !atom->mumag_flag)
|
||||
error->all(FLERR,"Pair spin requires atom attributes sp, mumag");
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
// check this half/full request
|
||||
#define FULLNEI
|
||||
#if defined FULLNEI
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairSpinSocDmi::init_one(int i, int j)
|
||||
{
|
||||
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
|
||||
return cut_spin_dmi_global;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocDmi::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]) {
|
||||
if (dmi_flag) {
|
||||
fwrite(&DM[i][j],sizeof(double),1,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(&cut_spin_dmi[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocDmi::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(&DM[i][j],sizeof(double),1,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(&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(&cut_spin_dmi[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocDmi::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_spin_dmi_global,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocDmi::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&cut_spin_dmi_global,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&cut_spin_dmi_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
|
@ -0,0 +1,87 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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(pair/spin/soc/dmi,PairSpinSocDmi)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_SPIN_SOC_DMI_H
|
||||
#define LMP_PAIR_SPIN_SOC_DMI_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairSpinSocDmi : public Pair {
|
||||
public:
|
||||
PairSpinSocDmi(class LAMMPS *);
|
||||
virtual ~PairSpinSocDmi();
|
||||
virtual void compute(int, int);
|
||||
void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
void compute_dmi(int, int, double *, double *, double *, double *);
|
||||
|
||||
int dmi_flag; // dmi flag
|
||||
|
||||
double cut_spin_dmi_global; // short range pair cutoff
|
||||
double **cut_spin_dmi; // cutoff distance dmi
|
||||
|
||||
protected:
|
||||
int newton_pair_spin;
|
||||
double hbar;
|
||||
|
||||
double **DM; // dmi coeff in eV
|
||||
double **v_dmx, **v_dmy, **v_dmz;// dmi direction
|
||||
|
||||
double *spi, *spj; // temp. spin vals. in compute
|
||||
double *fi, *fj; // temp. mech. forces in compute
|
||||
double *fmi, *fmj; // temp. mag. forces in compute
|
||||
double *rij; // norm. inter atomic vectors
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Incorrect args in pair_spin command
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Spin simulations require metal unit style
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair spin requires atom attributes sp, mumag
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,498 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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)
|
||||
Aidan Thompson (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "pair_hybrid.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "pair_spin_soc_landau.h"
|
||||
#include "update.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpinSocLandau::PairSpinSocLandau(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
hbar = force->hplanck/MY_2PI;
|
||||
|
||||
newton_pair_spin = 0; // no newton pair for now
|
||||
// newton_pair = 0;
|
||||
|
||||
single_enable = 0;
|
||||
soc_neel_flag = 0;
|
||||
|
||||
no_virial_fdotr_compute = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpinSocLandau::~PairSpinSocLandau()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
|
||||
memory->destroy(cut_soc_neel);
|
||||
memory->destroy(K1);
|
||||
memory->destroy(K1_mech);
|
||||
memory->destroy(K2);
|
||||
memory->destroy(K3);
|
||||
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fi);
|
||||
memory->destroy(fj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
memory->destroy(rij);
|
||||
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocLandau::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double evdwl,ecoul;
|
||||
double xi,yi,zi;
|
||||
double fix,fiy,fiz,fjx,fjy,fjz;
|
||||
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
|
||||
double cut_soc_neel_2,cut_soc_global2;
|
||||
double rsq,rd;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
cut_soc_global2 = cut_soc_global*cut_soc_global;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double **fm = atom->fm;
|
||||
double *mumag = atom->mumag;
|
||||
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;
|
||||
|
||||
// pair spin computations
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xi = x[i][0];
|
||||
yi = x[i][1];
|
||||
zi = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
spi[0] = sp[i][0];
|
||||
spi[1] = sp[i][1];
|
||||
spi[2] = sp[i][2];
|
||||
|
||||
// loop on neighbors
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
spj[0] = sp[j][0];
|
||||
spj[1] = sp[j][1];
|
||||
spj[2] = sp[j][2];
|
||||
|
||||
evdwl = 0.0;
|
||||
|
||||
fi[0] = fi[1] = fi[2] = 0.0;
|
||||
fj[0] = fj[1] = fj[2] = 0.0;
|
||||
fmi[0] = fmi[1] = fmi[2] = 0.0;
|
||||
fmj[0] = fmj[1] = fmj[2] = 0.0;
|
||||
rij[0] = rij[1] = rij[2] = 0.0;
|
||||
|
||||
rij[0] = x[j][0] - xi;
|
||||
rij[1] = x[j][1] - yi;
|
||||
rij[2] = x[j][2] - zi;
|
||||
|
||||
// square of inter-atomic distance
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
double inorm = 1.0/sqrt(rsq);
|
||||
rij[0] *= inorm;
|
||||
rij[1] *= inorm;
|
||||
rij[2] *= inorm;
|
||||
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
// compute mag. and mech. components of soc
|
||||
cut_soc_neel_2 = cut_soc_neel[itype][jtype]*cut_soc_neel[itype][jtype];
|
||||
if (rsq <= cut_soc_neel_2) {
|
||||
compute_soc_neel(i,j,rsq,rij,fmi,fmj,spi,spj);
|
||||
compute_soc_mech_neel(i,j,rsq,rij,fi,fj,spi,spj);
|
||||
}
|
||||
|
||||
f[i][0] += fi[0];
|
||||
f[i][1] += fi[1];
|
||||
f[i][2] += fi[2];
|
||||
fm[i][0] += fmi[0];
|
||||
fm[i][1] += fmi[1];
|
||||
fm[i][2] += fmi[2];
|
||||
|
||||
// if (newton_pair || j < nlocal) {
|
||||
if (newton_pair_spin) {
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
f[j][2] += fj[2];
|
||||
fm[j][0] += fmj[0];
|
||||
fm[j][1] += fmj[1];
|
||||
fm[j][2] += fmj[2];
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq <= cut_soc_neel_2) {
|
||||
evdwl -= spi[0]*fmi[0];
|
||||
evdwl -= spi[1]*fmi[1];
|
||||
evdwl -= 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]);
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinSocLandau::compute_soc_neel(int i, int j, double rsq, double *rij, double *fmi, double *fmj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
double Kij, Kij_3, ra, scalar;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
ra = rsq/K3[itype][jtype]/K3[itype][jtype];
|
||||
Kij = 4.0*K1[itype][jtype]*ra;
|
||||
Kij *= (1.0-K2[itype][jtype]*ra);
|
||||
Kij *= exp(-ra);
|
||||
|
||||
scalar = rij[0]*spj[0]+rij[1]*spj[1]+rij[2]*spj[2];
|
||||
Kij_3 = Kij/3.0;
|
||||
|
||||
fmi[0] += Kij*scalar*rij[0]-Kij_3*spj[0];
|
||||
fmi[1] += Kij*scalar*rij[1]-Kij_3*spj[1];
|
||||
fmi[2] += Kij*scalar*rij[2]-Kij_3*spj[2];
|
||||
|
||||
fmj[0] -= Kij*scalar*rij[0]+Kij_3*spi[0];
|
||||
fmj[1] -= Kij*scalar*rij[1]+Kij_3*spi[1];
|
||||
fmj[2] -= Kij*scalar*rij[2]+Kij_3*spi[2];
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinSocLandau::compute_soc_mech_neel(int i, int j, double rsq, double *rij, double *fi, double *fj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
double scalar_si_sj, scalar_rij_si, scalar_rij_sj;
|
||||
double K_mech, Kij, dKij, ra, rr, drij, iK3;
|
||||
double t1, t2, t3;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
scalar_si_sj = spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2];
|
||||
scalar_rij_si = rij[0]*spi[0]+rij[1]*spi[1]+rij[2]*spi[2];
|
||||
scalar_rij_sj = rij[0]*spj[0]+rij[1]*spj[1]+rij[2]*spj[2];
|
||||
|
||||
K_mech = K1_mech[itype][jtype];
|
||||
iK3 = 1.0/(K3[itype][jtype]*K3[itype][jtype]);
|
||||
|
||||
drij = sqrt(rsq);
|
||||
ra = rsq*iK3;
|
||||
rr = drij*iK3;
|
||||
|
||||
Kij *= (1.0-K2[itype][jtype]*ra);
|
||||
Kij *= 4.0*K_mech*ra*exp(-ra);
|
||||
|
||||
dKij = 1.0-ra-K2[itype][jtype]*ra*(2.0-ra);
|
||||
dKij *= 8.0*K_mech*rr*exp(-ra);
|
||||
|
||||
t1 = (dKij-2.0*Kij/drij)*scalar_rij_si*scalar_rij_sj;
|
||||
t1 -= scalar_si_sj*dKij/3.0;
|
||||
t2 = scalar_rij_sj*Kij/drij;
|
||||
t3 = scalar_rij_si*Kij/drij;
|
||||
|
||||
fi[0] += t1*rij[0]+t2*spi[0]+t3*spj[0];
|
||||
fi[1] += t1*rij[1]+t2*spi[1]+t3*spj[1];
|
||||
fi[2] += t1*rij[2]+t2*spi[2]+t3*spj[2];
|
||||
|
||||
fj[0] -= t1*rij[0]-t2*spi[0]-t3*spj[0];
|
||||
fj[1] -= t1*rij[1]-t2*spi[1]-t3*spj[1];
|
||||
fj[2] -= t1*rij[2]-t2*spi[2]-t3*spj[2];
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocLandau::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_soc_neel,n+1,n+1,"pair:cut_soc_neel");
|
||||
memory->create(K1,n+1,n+1,"pair:K1");
|
||||
memory->create(K1_mech,n+1,n+1,"pair:K1_mech");
|
||||
memory->create(K2,n+1,n+1,"pair:K2");
|
||||
memory->create(K3,n+1,n+1,"pair:K3");
|
||||
|
||||
memory->create(spi,3,"pair:spi");
|
||||
memory->create(spj,3,"pair:spj");
|
||||
memory->create(fi,3,"pair:fi");
|
||||
memory->create(fj,3,"pair:fj");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(fmj,3,"pair:fmj");
|
||||
memory->create(rij,3,"pair:rij");
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocLandau::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 1 || narg > 2)
|
||||
error->all(FLERR,"Incorrect number of args in pair_style pair/spin command");
|
||||
|
||||
if (strcmp(update->unit_style,"metal") != 0)
|
||||
error->all(FLERR,"Spin simulations require metal unit style");
|
||||
|
||||
cut_soc_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_soc_neel[i][j] = cut_soc_global;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type spin pairs (only one for now)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocLandau::coeff(int narg, char **arg)
|
||||
{
|
||||
const double hbar = force->hplanck/MY_2PI;
|
||||
|
||||
if (!allocated) allocate();
|
||||
|
||||
// set mech_flag to 1 if magneto-mech simulation
|
||||
//no longer correct: can be hybrid without magneto-mech
|
||||
if (strstr(force->pair_style,"pair/spin")) {
|
||||
mech_flag = 0;
|
||||
} else if (strstr(force->pair_style,"hybrid/overlay")) {
|
||||
mech_flag = 1;
|
||||
} else error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
|
||||
if (strcmp(arg[2],"neel")==0){
|
||||
if (narg != 7) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
soc_neel_flag = 1;
|
||||
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
|
||||
const double rij = force->numeric(FLERR,arg[3]);
|
||||
const double k1 = (force->numeric(FLERR,arg[4]));
|
||||
const double k2 = force->numeric(FLERR,arg[5]);
|
||||
const double k3 = force->numeric(FLERR,arg[6]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
cut_soc_neel[i][j] = rij;
|
||||
K1[i][j] = k1/hbar;
|
||||
if (mech_flag) {
|
||||
K1_mech[i][j] = k1;
|
||||
} else {
|
||||
K1_mech[i][j] = 0.0;
|
||||
}
|
||||
K2[i][j] = k2;
|
||||
K3[i][j] = k3;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
} else error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocLandau::init_style()
|
||||
{
|
||||
if (!atom->sp_flag || !atom->mumag_flag)
|
||||
error->all(FLERR,"Pair spin requires atom attributes sp, mumag");
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
// check this half/full request
|
||||
#define FULLNEI
|
||||
#if defined FULLNEI
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairSpinSocLandau::init_one(int i, int j)
|
||||
{
|
||||
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
|
||||
return cut_soc_global;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocLandau::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]) {
|
||||
if (soc_neel_flag){
|
||||
fwrite(&K1[i][j],sizeof(double),1,fp);
|
||||
fwrite(&K1_mech[i][j],sizeof(double),1,fp);
|
||||
fwrite(&K2[i][j],sizeof(double),1,fp);
|
||||
fwrite(&K3[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_soc_neel[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocLandau::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(&K1[i][j],sizeof(double),1,fp);
|
||||
fread(&K1_mech[i][j],sizeof(double),1,fp);
|
||||
fread(&K2[i][j],sizeof(double),1,fp);
|
||||
fread(&K2[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_soc_neel[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&K1[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&K1_mech[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&K2[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&K3[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_soc_neel[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocLandau::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_soc_global,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocLandau::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&cut_soc_global,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&cut_soc_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
|
@ -0,0 +1,89 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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(pair/spin/soc/landau,PairSpinSocLandau)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_SPIN_SOC_LANDAU_H
|
||||
#define LMP_PAIR_SPIN_SOC_LANDAU_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairSpinSocLandau : public Pair {
|
||||
public:
|
||||
PairSpinSocLandau(class LAMMPS *);
|
||||
virtual ~PairSpinSocLandau();
|
||||
virtual void compute(int, int);
|
||||
void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
void compute_soc_neel(int, int, double, double *, double *, double *,double *, double *);
|
||||
void compute_soc_mech_neel(int, int, double, double *, double *, double *,double *, double *);
|
||||
|
||||
int soc_neel_flag; // soc neel flag
|
||||
int mech_flag; // mech calc. flag
|
||||
|
||||
double cut_soc_global;
|
||||
double **cut_soc_neel; // cutoff distance exchange
|
||||
|
||||
protected:
|
||||
int newton_pair_spin;
|
||||
double hbar;
|
||||
|
||||
double **K1, **K1_mech; // exchange coeffs Kij
|
||||
double **K2, **K3; // K1 in eV, K2 adim, K3 in Ang
|
||||
|
||||
double *spi, *spj; // temp. spin vals. in compute
|
||||
double *fi, *fj; // temp. mech. forces in compute
|
||||
double *fmi, *fmj; // temp. mag. forces in compute
|
||||
double *rij; // norm. inter atomic vectors
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Incorrect args in pair_spin command
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Spin simulations require metal unit style
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair spin requires atom attributes sp, mumag
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
*/
|
||||
|
|
@ -0,0 +1,498 @@
|
|||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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)
|
||||
Aidan Thompson (SNL)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "pair_hybrid.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neigh_request.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "pair_spin_soc_neel.h"
|
||||
#include "update.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpinSocNeel::PairSpinSocNeel(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
hbar = force->hplanck/MY_2PI;
|
||||
|
||||
newton_pair_spin = 0; // no newton pair for now
|
||||
// newton_pair = 0;
|
||||
|
||||
single_enable = 0;
|
||||
soc_neel_flag = 0;
|
||||
|
||||
no_virial_fdotr_compute = 1;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairSpinSocNeel::~PairSpinSocNeel()
|
||||
{
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
|
||||
memory->destroy(cut_soc_neel);
|
||||
memory->destroy(K1);
|
||||
memory->destroy(K1_mech);
|
||||
memory->destroy(K2);
|
||||
memory->destroy(K3);
|
||||
|
||||
memory->destroy(spi);
|
||||
memory->destroy(spj);
|
||||
memory->destroy(fi);
|
||||
memory->destroy(fj);
|
||||
memory->destroy(fmi);
|
||||
memory->destroy(fmj);
|
||||
memory->destroy(rij);
|
||||
|
||||
memory->destroy(cutsq);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocNeel::compute(int eflag, int vflag)
|
||||
{
|
||||
int i,j,ii,jj,inum,jnum,itype,jtype;
|
||||
double evdwl,ecoul;
|
||||
double xi,yi,zi;
|
||||
double fix,fiy,fiz,fjx,fjy,fjz;
|
||||
double fmix,fmiy,fmiz,fmjx,fmjy,fmjz;
|
||||
double cut_soc_neel_2,cut_soc_global2;
|
||||
double rsq,rd;
|
||||
int *ilist,*jlist,*numneigh,**firstneigh;
|
||||
|
||||
evdwl = ecoul = 0.0;
|
||||
if (eflag || vflag) ev_setup(eflag,vflag);
|
||||
else evflag = vflag_fdotr = 0;
|
||||
cut_soc_global2 = cut_soc_global*cut_soc_global;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
double **fm = atom->fm;
|
||||
double *mumag = atom->mumag;
|
||||
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;
|
||||
|
||||
// pair spin computations
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xi = x[i][0];
|
||||
yi = x[i][1];
|
||||
zi = x[i][2];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
spi[0] = sp[i][0];
|
||||
spi[1] = sp[i][1];
|
||||
spi[2] = sp[i][2];
|
||||
|
||||
// loop on neighbors
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
spj[0] = sp[j][0];
|
||||
spj[1] = sp[j][1];
|
||||
spj[2] = sp[j][2];
|
||||
|
||||
evdwl = 0.0;
|
||||
|
||||
fi[0] = fi[1] = fi[2] = 0.0;
|
||||
fj[0] = fj[1] = fj[2] = 0.0;
|
||||
fmi[0] = fmi[1] = fmi[2] = 0.0;
|
||||
fmj[0] = fmj[1] = fmj[2] = 0.0;
|
||||
rij[0] = rij[1] = rij[2] = 0.0;
|
||||
|
||||
rij[0] = x[j][0] - xi;
|
||||
rij[1] = x[j][1] - yi;
|
||||
rij[2] = x[j][2] - zi;
|
||||
|
||||
// square of inter-atomic distance
|
||||
rsq = rij[0]*rij[0] + rij[1]*rij[1] + rij[2]*rij[2];
|
||||
double inorm = 1.0/sqrt(rsq);
|
||||
rij[0] *= inorm;
|
||||
rij[1] *= inorm;
|
||||
rij[2] *= inorm;
|
||||
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
// compute mag. and mech. components of soc
|
||||
cut_soc_neel_2 = cut_soc_neel[itype][jtype]*cut_soc_neel[itype][jtype];
|
||||
if (rsq <= cut_soc_neel_2) {
|
||||
compute_soc_neel(i,j,rsq,rij,fmi,fmj,spi,spj);
|
||||
compute_soc_mech_neel(i,j,rsq,rij,fi,fj,spi,spj);
|
||||
}
|
||||
|
||||
f[i][0] += fi[0];
|
||||
f[i][1] += fi[1];
|
||||
f[i][2] += fi[2];
|
||||
fm[i][0] += fmi[0];
|
||||
fm[i][1] += fmi[1];
|
||||
fm[i][2] += fmi[2];
|
||||
|
||||
// if (newton_pair || j < nlocal) {
|
||||
if (newton_pair_spin) {
|
||||
f[j][0] += fj[0];
|
||||
f[j][1] += fj[1];
|
||||
f[j][2] += fj[2];
|
||||
fm[j][0] += fmj[0];
|
||||
fm[j][1] += fmj[1];
|
||||
fm[j][2] += fmj[2];
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
if (rsq <= cut_soc_neel_2) {
|
||||
evdwl -= spi[0]*fmi[0];
|
||||
evdwl -= spi[1]*fmi[1];
|
||||
evdwl -= 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]);
|
||||
}
|
||||
}
|
||||
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinSocNeel::compute_soc_neel(int i, int j, double rsq, double *rij, double *fmi, double *fmj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
double Kij, Kij_3, ra, scalar;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
ra = rsq/K3[itype][jtype]/K3[itype][jtype];
|
||||
Kij = 4.0*K1[itype][jtype]*ra;
|
||||
Kij *= (1.0-K2[itype][jtype]*ra);
|
||||
Kij *= exp(-ra);
|
||||
|
||||
scalar = rij[0]*spj[0]+rij[1]*spj[1]+rij[2]*spj[2];
|
||||
Kij_3 = Kij/3.0;
|
||||
|
||||
fmi[0] += Kij*scalar*rij[0]-Kij_3*spj[0];
|
||||
fmi[1] += Kij*scalar*rij[1]-Kij_3*spj[1];
|
||||
fmi[2] += Kij*scalar*rij[2]-Kij_3*spj[2];
|
||||
|
||||
fmj[0] -= Kij*scalar*rij[0]+Kij_3*spi[0];
|
||||
fmj[1] -= Kij*scalar*rij[1]+Kij_3*spi[1];
|
||||
fmj[2] -= Kij*scalar*rij[2]+Kij_3*spi[2];
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
void PairSpinSocNeel::compute_soc_mech_neel(int i, int j, double rsq, double *rij, double *fi, double *fj, double *spi, double *spj)
|
||||
{
|
||||
int *type = atom->type;
|
||||
int itype, jtype;
|
||||
double scalar_si_sj, scalar_rij_si, scalar_rij_sj;
|
||||
double K_mech, Kij, dKij, ra, rr, drij, iK3;
|
||||
double t1, t2, t3;
|
||||
itype = type[i];
|
||||
jtype = type[j];
|
||||
|
||||
scalar_si_sj = spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2];
|
||||
scalar_rij_si = rij[0]*spi[0]+rij[1]*spi[1]+rij[2]*spi[2];
|
||||
scalar_rij_sj = rij[0]*spj[0]+rij[1]*spj[1]+rij[2]*spj[2];
|
||||
|
||||
K_mech = K1_mech[itype][jtype];
|
||||
iK3 = 1.0/(K3[itype][jtype]*K3[itype][jtype]);
|
||||
|
||||
drij = sqrt(rsq);
|
||||
ra = rsq*iK3;
|
||||
rr = drij*iK3;
|
||||
|
||||
Kij *= (1.0-K2[itype][jtype]*ra);
|
||||
Kij *= 4.0*K_mech*ra*exp(-ra);
|
||||
|
||||
dKij = 1.0-ra-K2[itype][jtype]*ra*(2.0-ra);
|
||||
dKij *= 8.0*K_mech*rr*exp(-ra);
|
||||
|
||||
t1 = (dKij-2.0*Kij/drij)*scalar_rij_si*scalar_rij_sj;
|
||||
t1 -= scalar_si_sj*dKij/3.0;
|
||||
t2 = scalar_rij_sj*Kij/drij;
|
||||
t3 = scalar_rij_si*Kij/drij;
|
||||
|
||||
fi[0] += t1*rij[0]+t2*spi[0]+t3*spj[0];
|
||||
fi[1] += t1*rij[1]+t2*spi[1]+t3*spj[1];
|
||||
fi[2] += t1*rij[2]+t2*spi[2]+t3*spj[2];
|
||||
|
||||
fj[0] -= t1*rij[0]-t2*spi[0]-t3*spj[0];
|
||||
fj[1] -= t1*rij[1]-t2*spi[1]-t3*spj[1];
|
||||
fj[2] -= t1*rij[2]-t2*spi[2]-t3*spj[2];
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocNeel::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_soc_neel,n+1,n+1,"pair:cut_soc_neel");
|
||||
memory->create(K1,n+1,n+1,"pair:K1");
|
||||
memory->create(K1_mech,n+1,n+1,"pair:K1_mech");
|
||||
memory->create(K2,n+1,n+1,"pair:K2");
|
||||
memory->create(K3,n+1,n+1,"pair:K3");
|
||||
|
||||
memory->create(spi,3,"pair:spi");
|
||||
memory->create(spj,3,"pair:spj");
|
||||
memory->create(fi,3,"pair:fi");
|
||||
memory->create(fj,3,"pair:fj");
|
||||
memory->create(fmi,3,"pair:fmi");
|
||||
memory->create(fmj,3,"pair:fmj");
|
||||
memory->create(rij,3,"pair:rij");
|
||||
|
||||
memory->create(cutsq,n+1,n+1,"pair:cutsq");
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocNeel::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg < 1 || narg > 2)
|
||||
error->all(FLERR,"Incorrect number of args in pair_style pair/spin command");
|
||||
|
||||
if (strcmp(update->unit_style,"metal") != 0)
|
||||
error->all(FLERR,"Spin simulations require metal unit style");
|
||||
|
||||
cut_soc_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_soc_neel[i][j] = cut_soc_global;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type spin pairs (only one for now)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocNeel::coeff(int narg, char **arg)
|
||||
{
|
||||
const double hbar = force->hplanck/MY_2PI;
|
||||
|
||||
if (!allocated) allocate();
|
||||
|
||||
// set mech_flag to 1 if magneto-mech simulation
|
||||
//no longer correct: can be hybrid without magneto-mech
|
||||
if (strstr(force->pair_style,"pair/spin")) {
|
||||
mech_flag = 0;
|
||||
} else if (strstr(force->pair_style,"hybrid/overlay")) {
|
||||
mech_flag = 1;
|
||||
} else error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
|
||||
if (strcmp(arg[2],"neel")==0){
|
||||
if (narg != 7) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
soc_neel_flag = 1;
|
||||
|
||||
int ilo,ihi,jlo,jhi;
|
||||
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
|
||||
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
|
||||
|
||||
const double rij = force->numeric(FLERR,arg[3]);
|
||||
const double k1 = (force->numeric(FLERR,arg[4]));
|
||||
const double k2 = force->numeric(FLERR,arg[5]);
|
||||
const double k3 = force->numeric(FLERR,arg[6]);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo,i); j <= jhi; j++) {
|
||||
cut_soc_neel[i][j] = rij;
|
||||
K1[i][j] = k1/hbar;
|
||||
if (mech_flag) {
|
||||
K1_mech[i][j] = k1;
|
||||
} else {
|
||||
K1_mech[i][j] = 0.0;
|
||||
}
|
||||
K2[i][j] = k2;
|
||||
K3[i][j] = k3;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
if (count == 0) error->all(FLERR,"Incorrect args in pair_style command");
|
||||
} else error->all(FLERR,"Incorrect args in pair_style command");
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocNeel::init_style()
|
||||
{
|
||||
if (!atom->sp_flag || !atom->mumag_flag)
|
||||
error->all(FLERR,"Pair spin requires atom attributes sp, mumag");
|
||||
|
||||
neighbor->request(this,instance_me);
|
||||
|
||||
// check this half/full request
|
||||
#define FULLNEI
|
||||
#if defined FULLNEI
|
||||
int irequest = neighbor->request(this,instance_me);
|
||||
neighbor->requests[irequest]->half = 0;
|
||||
neighbor->requests[irequest]->full = 1;
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairSpinSocNeel::init_one(int i, int j)
|
||||
{
|
||||
|
||||
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
|
||||
|
||||
return cut_soc_global;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocNeel::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]) {
|
||||
if (soc_neel_flag){
|
||||
fwrite(&K1[i][j],sizeof(double),1,fp);
|
||||
fwrite(&K1_mech[i][j],sizeof(double),1,fp);
|
||||
fwrite(&K2[i][j],sizeof(double),1,fp);
|
||||
fwrite(&K3[i][j],sizeof(double),1,fp);
|
||||
fwrite(&cut_soc_neel[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocNeel::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(&K1[i][j],sizeof(double),1,fp);
|
||||
fread(&K1_mech[i][j],sizeof(double),1,fp);
|
||||
fread(&K2[i][j],sizeof(double),1,fp);
|
||||
fread(&K2[i][j],sizeof(double),1,fp);
|
||||
fread(&cut_soc_neel[i][j],sizeof(double),1,fp);
|
||||
}
|
||||
MPI_Bcast(&K1[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&K1_mech[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&K2[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&K3[i][j],1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&cut_soc_neel[i][j],1,MPI_DOUBLE,0,world);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocNeel::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_soc_global,sizeof(double),1,fp);
|
||||
fwrite(&offset_flag,sizeof(int),1,fp);
|
||||
fwrite(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairSpinSocNeel::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
fread(&cut_soc_global,sizeof(double),1,fp);
|
||||
fread(&offset_flag,sizeof(int),1,fp);
|
||||
fread(&mix_flag,sizeof(int),1,fp);
|
||||
}
|
||||
MPI_Bcast(&cut_soc_global,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
|
||||
}
|
||||
|
|
@ -0,0 +1,89 @@
|
|||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
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(pair/spin/soc/neel,PairSpinSocNeel)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_SPIN_SOC_NEEL_H
|
||||
#define LMP_PAIR_SPIN_SOC_NEEL_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairSpinSocNeel : public Pair {
|
||||
public:
|
||||
PairSpinSocNeel(class LAMMPS *);
|
||||
virtual ~PairSpinSocNeel();
|
||||
virtual void compute(int, int);
|
||||
void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
|
||||
void compute_soc_neel(int, int, double, double *, double *, double *,double *, double *);
|
||||
void compute_soc_mech_neel(int, int, double, double *, double *, double *,double *, double *);
|
||||
|
||||
int soc_neel_flag; // soc neel flag
|
||||
int mech_flag; // mech calc. flag
|
||||
|
||||
double cut_soc_global;
|
||||
double **cut_soc_neel; // cutoff distance exchange
|
||||
|
||||
protected:
|
||||
int newton_pair_spin;
|
||||
double hbar;
|
||||
|
||||
double **K1, **K1_mech; // exchange coeffs Kij
|
||||
double **K2, **K3; // K1 in eV, K2 adim, K3 in Ang
|
||||
|
||||
double *spi, *spj; // temp. spin vals. in compute
|
||||
double *fi, *fj; // temp. mech. forces in compute
|
||||
double *fmi, *fmj; // temp. mag. forces in compute
|
||||
double *rij; // norm. inter atomic vectors
|
||||
|
||||
void allocate();
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Incorrect args in pair_spin command
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Spin simulations require metal unit style
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Incorrect args for pair coefficients
|
||||
|
||||
Self-explanatory. Check the input script or data file.
|
||||
|
||||
E: Pair spin requires atom attributes sp, mumag
|
||||
|
||||
The atom style defined does not have these attributes.
|
||||
|
||||
*/
|
||||
|
|
@ -312,6 +312,7 @@ void Verlet::run(int n)
|
|||
timer->stamp(Timer::PAIR);
|
||||
}
|
||||
|
||||
|
||||
if (atom->molecular) {
|
||||
if (force->bond) force->bond->compute(eflag,vflag);
|
||||
if (force->angle) force->angle->compute(eflag,vflag);
|
||||
|
|
|
|||
|
|
@ -70,7 +70,7 @@ proc enable_trace {} {
|
|||
trace variable vmd_frame([molinfo top]) w vmd_draw_spin
|
||||
}
|
||||
|
||||
set molid [mol addfile {/home/jtranch/Documents/lammps/src/dump_spin_MM.lammpstrj} type {lammpstrj} autobonds off first 0 last -1 step 1 waitfor all]
|
||||
set molid [mol addfile {/home/jtranch/Documents/lammps/src/dump_VSRSV.lammpstrj} type {lammpstrj} autobonds off first 0 last -1 step 1 waitfor all]
|
||||
scale by 0.5
|
||||
animate style Loop
|
||||
enable_trace
|
||||
|
|
|
|||
Loading…
Reference in New Issue