mirror of https://github.com/lammps/lammps.git
Commit JT 040519
- initial rotation with Rodrigues' formula + exception - worked on neb_spin documentation - removed pair/spin warning for neb/spin
This commit is contained in:
julient31
parent
e45e92b1cb
commit
da16a7e50b
Binary file not shown.
|
After Width: | Height: | Size: 8.2 KiB |
|
|
@ -0,0 +1,16 @@
|
|||
\documentclass[preview]{standalone}
|
||||
\usepackage{varwidth}
|
||||
\usepackage[utf8x]{inputenc}
|
||||
\usepackage{amsmath, amssymb, graphics, setspace}
|
||||
|
||||
\begin{document}
|
||||
\begin{varwidth}{50in}
|
||||
\begin{equation}
|
||||
\vec{k}_i =
|
||||
\frac{\vec{m}_i^I \times \vec{m}_i^F}{\left|\vec{m}_i^I
|
||||
\times \vec{m}_i^F\right|}
|
||||
%&{\rm ~if~}& \vec{m}_i^I \times \vec{m}_i^F
|
||||
, \nonumber
|
||||
\end{equation}
|
||||
\end{varwidth}
|
||||
\end{document}
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 20 KiB |
|
|
@ -0,0 +1,16 @@
|
|||
\documentclass[preview]{standalone}
|
||||
\usepackage{varwidth}
|
||||
\usepackage[utf8x]{inputenc}
|
||||
\usepackage{amsmath, amssymb, graphics, setspace}
|
||||
|
||||
\begin{document}
|
||||
\begin{varwidth}{50in}
|
||||
\begin{equation}
|
||||
\vec{m}_i^{\nu} = \vec{m}_i^{I} \cos(\omega_i^{\nu})
|
||||
+ (\vec{k}_i \times \vec{m}_i^{I}) \sin(\omega_i^{\nu})
|
||||
+ (1.0-\cos(\omega_i^{\nu})) \vec{k}_i (\vec{k}_i\cdot
|
||||
\vec{m}_i^{I})
|
||||
, \nonumber
|
||||
\end{equation}
|
||||
\end{varwidth}
|
||||
\end{document}
|
||||
|
|
@ -0,0 +1,430 @@
|
|||
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
|
||||
|
||||
:link(lws,http://lammps.sandia.gov)
|
||||
:link(ld,Manual.html)
|
||||
:link(lc,Commands_all.html)
|
||||
|
||||
:line
|
||||
|
||||
neb command :h3
|
||||
|
||||
[Syntax:]
|
||||
|
||||
neb/spin etol ttol N1 N2 Nevery file-style arg keyword :pre
|
||||
|
||||
etol = stopping tolerance for energy (energy units) :ulb,l
|
||||
ttol = stopping tolerance for torque ( units) :l
|
||||
N1 = max # of iterations (timesteps) to run initial NEB :l
|
||||
N2 = max # of iterations (timesteps) to run barrier-climbing NEB :l
|
||||
Nevery = print replica energies and reaction coordinates every this many timesteps :l
|
||||
file-style = {final} or {each} or {none} :l
|
||||
{final} arg = filename
|
||||
filename = file with initial coords for final replica
|
||||
coords for intermediate replicas are linearly interpolated
|
||||
between first and last replica
|
||||
{each} arg = filename
|
||||
filename = unique filename for each replica (except first)
|
||||
with its initial coords
|
||||
{none} arg = no argument all replicas assumed to already have
|
||||
their initial coords :pre
|
||||
keyword = {verbose}
|
||||
:ule
|
||||
|
||||
[Examples:]
|
||||
|
||||
neb/spin 0.1 0.0 1000 500 50 final coords.final
|
||||
neb/spin 0.0 0.001 1000 500 50 each coords.initial.$i
|
||||
neb/spin 0.0 0.001 1000 500 50 none verbose :pre
|
||||
|
||||
[Description:]
|
||||
|
||||
Perform a geodesic nudged elastic band (GNEB) calculation using multiple
|
||||
replicas of a system. Two or more replicas must be used; the first
|
||||
and last are the end points of the transition path.
|
||||
|
||||
GNEB is a method for finding both the spin configurations and height
|
||||
of the energy barrier associated with a transition state, e.g.
|
||||
spins to perform a collective rotation from one energy basin to
|
||||
another.
|
||||
The implementation in LAMMPS follows the discussion in the
|
||||
following paper: "(Bessarab)"_#Bessarab.
|
||||
|
||||
Each replica runs on a partition of one or more processors. Processor
|
||||
partitions are defined at run-time using the "-partition command-line
|
||||
switch"_Run_options.html. Note that if you have MPI installed, you
|
||||
can run a multi-replica simulation with more replicas (partitions)
|
||||
than you have physical processors, e.g you can run a 10-replica
|
||||
simulation on just one or two processors. You will simply not get the
|
||||
performance speed-up you would see with one or more physical
|
||||
processors per replica. See the "Howto replica"_Howto_replica.html
|
||||
doc page for further discussion.
|
||||
|
||||
NOTE: As explained below, a GNEB calculation performs a damped dynamics
|
||||
minimization across all the replicas. The "spin"_min_spin.html
|
||||
style minimizer has to be defined in your input script.
|
||||
|
||||
When a GNEB calculation is performed, it is assumed that each replica
|
||||
is running the same system, though LAMMPS does not check for this.
|
||||
I.e. the simulation domain, the number of magnetic atoms, the
|
||||
interaction potentials, and the starting configuration when the neb
|
||||
command is issued should be the same for every replica.
|
||||
|
||||
In a GNEB calculation each replica is connected to other replicas by
|
||||
inter-replica nudging forces. These forces are imposed by the "fix
|
||||
neb/spin"_fix_neb_spin.html command, which must be used in conjunction
|
||||
with the neb command.
|
||||
The group used to define the fix neb/spin command defines the
|
||||
GNEB magnetic atoms which are the only ones that inter-replica springs
|
||||
are applied to.
|
||||
If the group does not include all magnetic atoms, then non-GNEB
|
||||
magnetic atoms have no inter-replica springs and the torques they feel
|
||||
and their precessional motion is computed in the usual way due only
|
||||
to other magnetic atoms within their replica.
|
||||
Conceptually, the non-GNEB atoms provide a background force field for
|
||||
the GNEB atoms.
|
||||
Their magnetic spins can be allowed to precess during the GNEB
|
||||
minimization procedure.
|
||||
|
||||
The initial spin configuration for each of the replicas can be
|
||||
specified in different manners via the {file-style} setting, as
|
||||
discussed below. Only atomic spins whose initial coordinates should
|
||||
differ from the current configuration need to be specified.
|
||||
|
||||
Conceptually, the initial and final configurations for the first
|
||||
replica should be states on either side of an energy barrier.
|
||||
|
||||
As explained below, the initial configurations of intermediate
|
||||
replicas can be spin coordinates interpolated in a linear fashion
|
||||
between the first and last replicas. This is often adequate for
|
||||
simple transitions. For more complex transitions, it may lead to slow
|
||||
convergence or even bad results if the minimum energy path (MEP, see
|
||||
below) of states over the barrier cannot be correctly converged to
|
||||
from such an initial path. In this case, you will want to generate
|
||||
initial states for the intermediate replicas that are geometrically
|
||||
closer to the MEP and read them in.
|
||||
|
||||
###################################################################
|
||||
|
||||
:line
|
||||
|
||||
For a {file-style} setting of {final}, a filename is specified which
|
||||
contains atomic and spin coordinates for zero or more atoms, in the
|
||||
format described below.
|
||||
For each atom that appears in the file, the new coordinates are
|
||||
assigned to that atom in the final replica. Each intermediate replica
|
||||
also assigns a new spin to that atom in an interpolated manner.
|
||||
This is done by using the current direction of the spin at the starting
|
||||
point and the read-in direction as the final point.
|
||||
The angular distance between them is calculated, and the new direction
|
||||
is assigned to be a fraction of the angular distance.
|
||||
|
||||
NOTE: The "angular distance" between the starting and final point is
|
||||
evaluated in the geodesic sense, as described in "(Bessarab)"_#Bessarab.
|
||||
|
||||
NOTE: The angular interpolation between the starting and final point
|
||||
is achieved using Rodrigues formula:
|
||||
|
||||
:c,image(Eqs/neb_spin_rodrigues_formula.jpg)
|
||||
|
||||
with m_i^I is the initial spin configuration for the spin i,
|
||||
where the rotation and k_i is defined as:
|
||||
|
||||
:c,image(Eqs/neb_spin_k.jpg)
|
||||
|
||||
The distance between them is calculated, and the new position
|
||||
is assigned to be a fraction of the distance. E.g. if there are 10
|
||||
replicas, the 2nd replica will assign a position that is 10% of the
|
||||
distance along a line between the starting and final point, and the
|
||||
9th replica will assign a position that is 90% of the distance along
|
||||
the line. Note that for this procedure to produce consistent
|
||||
coordinates across all the replicas, the current coordinates need to
|
||||
be the same in all replicas. LAMMPS does not check for this, but
|
||||
invalid initial configurations will likely result if it is not the
|
||||
case.
|
||||
|
||||
NOTE: The "distance" between the starting and final point is
|
||||
calculated in a minimum-image sense for a periodic simulation box.
|
||||
This means that if the two positions are on opposite sides of a box
|
||||
(periodic in that dimension), the distance between them will be small,
|
||||
because the periodic image of one of the atoms is close to the other.
|
||||
Similarly, even if the assigned position resulting from the
|
||||
interpolation is outside the periodic box, the atom will be wrapped
|
||||
back into the box when the NEB calculation begins.
|
||||
|
||||
For a {file-style} setting of {each}, a filename is specified which is
|
||||
assumed to be unique to each replica. This can be done by using a
|
||||
variable in the filename, e.g.
|
||||
|
||||
variable i equal part
|
||||
neb 0.0 0.001 1000 500 50 each coords.initial.$i :pre
|
||||
|
||||
which in this case will substitute the partition ID (0 to N-1) for the
|
||||
variable I, which is also effectively the replica ID. See the
|
||||
"variable"_variable.html command for other options, such as using
|
||||
world-, universe-, or uloop-style variables.
|
||||
|
||||
Each replica (except the first replica) will read its file, formatted
|
||||
as described below, and for any atom that appears in the file, assign
|
||||
the specified coordinates to its atom. The various files do not need
|
||||
to contain the same set of atoms.
|
||||
|
||||
For a {file-style} setting of {none}, no filename is specified. Each
|
||||
replica is assumed to already be in its initial configuration at the
|
||||
time the neb command is issued. This allows each replica to define
|
||||
its own configuration by reading a replica-specific data or restart or
|
||||
dump file, via the "read_data"_read_data.html,
|
||||
"read_restart"_read_restart.html, or "read_dump"_read_dump.html
|
||||
commands. The replica-specific names of these files can be specified
|
||||
as in the discussion above for the {each} file-style. Also see the
|
||||
section below for how a NEB calculation can produce restart files, so
|
||||
that a long calculation can be restarted if needed.
|
||||
|
||||
NOTE: None of the {file-style} settings change the initial
|
||||
configuration of any atom in the first replica. The first replica
|
||||
must thus be in the correct initial configuration at the time the neb
|
||||
command is issued.
|
||||
|
||||
:line
|
||||
|
||||
A NEB calculation proceeds in two stages, each of which is a
|
||||
minimization procedure, performed via damped dynamics. To enable
|
||||
this, you must first define a damped dynamics
|
||||
"min_style"_min_style.html, such as {quickmin} or {fire}. The {cg},
|
||||
{sd}, and {hftn} styles cannot be used, since they perform iterative
|
||||
line searches in their inner loop, which cannot be easily synchronized
|
||||
across multiple replicas.
|
||||
|
||||
The minimizer tolerances for energy and force are set by {etol} and
|
||||
{ftol}, the same as for the "minimize"_minimize.html command.
|
||||
|
||||
A non-zero {etol} means that the NEB calculation will terminate if the
|
||||
energy criterion is met by every replica. The energies being compared
|
||||
to {etol} do not include any contribution from the inter-replica
|
||||
nudging forces, since these are non-conservative. A non-zero {ftol}
|
||||
means that the NEB calculation will terminate if the force criterion
|
||||
is met by every replica. The forces being compared to {ftol} include
|
||||
the inter-replica nudging forces.
|
||||
|
||||
The maximum number of iterations in each stage is set by {N1} and
|
||||
{N2}. These are effectively timestep counts since each iteration of
|
||||
damped dynamics is like a single timestep in a dynamics
|
||||
"run"_run.html. During both stages, the potential energy of each
|
||||
replica and its normalized distance along the reaction path (reaction
|
||||
coordinate RD) will be printed to the screen and log file every
|
||||
{Nevery} timesteps. The RD is 0 and 1 for the first and last replica.
|
||||
For intermediate replicas, it is the cumulative distance (normalized
|
||||
by the total cumulative distance) between adjacent replicas, where
|
||||
"distance" is defined as the length of the 3N-vector of differences in
|
||||
atomic coordinates, where N is the number of NEB atoms involved in the
|
||||
transition. These outputs allow you to monitor NEB's progress in
|
||||
finding a good energy barrier. {N1} and {N2} must both be multiples
|
||||
of {Nevery}.
|
||||
|
||||
In the first stage of NEB, the set of replicas should converge toward
|
||||
a minimum energy path (MEP) of conformational states that transition
|
||||
over a barrier. The MEP for a transition is defined as a sequence of
|
||||
3N-dimensional states, each of which has a potential energy gradient
|
||||
parallel to the MEP itself. The configuration of highest energy along
|
||||
a MEP corresponds to a saddle point. The replica states will also be
|
||||
roughly equally spaced along the MEP due to the inter-replica nudging
|
||||
force added by the "fix neb"_fix_neb.html command.
|
||||
|
||||
In the second stage of NEB, the replica with the highest energy is
|
||||
selected and the inter-replica forces on it are converted to a force
|
||||
that drives its atom coordinates to the top or saddle point of the
|
||||
barrier, via the barrier-climbing calculation described in
|
||||
"(HenkelmanB)"_#HenkelmanB. As before, the other replicas rearrange
|
||||
themselves along the MEP so as to be roughly equally spaced.
|
||||
|
||||
When both stages are complete, if the NEB calculation was successful,
|
||||
the configurations of the replicas should be along (close to) the MEP
|
||||
and the replica with the highest energy should be an atomic
|
||||
configuration at (close to) the saddle point of the transition. The
|
||||
potential energies for the set of replicas represents the energy
|
||||
profile of the transition along the MEP.
|
||||
|
||||
:line
|
||||
|
||||
A few other settings in your input script are required or advised to
|
||||
perform a NEB calculation. See the NOTE about the choice of timestep
|
||||
at the beginning of this doc page.
|
||||
|
||||
An atom map must be defined which it is not by default for "atom_style
|
||||
atomic"_atom_style.html problems. The "atom_modify
|
||||
map"_atom_modify.html command can be used to do this.
|
||||
|
||||
The minimizers in LAMMPS operate on all atoms in your system, even
|
||||
non-NEB atoms, as defined above. To prevent non-NEB atoms from moving
|
||||
during the minimization, you should use the "fix
|
||||
setforce"_fix_setforce.html command to set the force on each of those
|
||||
atoms to 0.0. This is not required, and may not even be desired in
|
||||
some cases, but if those atoms move too far (e.g. because the initial
|
||||
state of your system was not well-minimized), it can cause problems
|
||||
for the NEB procedure.
|
||||
|
||||
The damped dynamics "minimizers"_min_style.html, such as {quickmin}
|
||||
and {fire}), adjust the position and velocity of the atoms via an
|
||||
Euler integration step. Thus you must define an appropriate
|
||||
"timestep"_timestep.html to use with NEB. As mentioned above, NEB
|
||||
will often converge more quickly if you use a timestep about 10x
|
||||
larger than you would normally use for dynamics simulations.
|
||||
|
||||
:line
|
||||
|
||||
Each file read by the neb command containing atomic coordinates used
|
||||
to initialize one or more replicas must be formatted as follows.
|
||||
|
||||
The file can be ASCII text or a gzipped text file (detected by a .gz
|
||||
suffix). The file can contain initial blank lines or comment lines
|
||||
starting with "#" which are ignored. The first non-blank, non-comment
|
||||
line should list N = the number of lines to follow. The N successive
|
||||
lines contain the following information:
|
||||
|
||||
ID1 x1 y1 z1
|
||||
ID2 x2 y2 z2
|
||||
...
|
||||
IDN xN yN zN :pre
|
||||
|
||||
The fields are the atom ID, followed by the x,y,z coordinates. The
|
||||
lines can be listed in any order. Additional trailing information on
|
||||
the line is OK, such as a comment.
|
||||
|
||||
Note that for a typical NEB calculation you do not need to specify
|
||||
initial coordinates for very many atoms to produce differing starting
|
||||
and final replicas whose intermediate replicas will converge to the
|
||||
energy barrier. Typically only new coordinates for atoms
|
||||
geometrically near the barrier need be specified.
|
||||
|
||||
Also note there is no requirement that the atoms in the file
|
||||
correspond to the NEB atoms in the group defined by the "fix
|
||||
neb"_fix_neb.html command. Not every NEB atom need be in the file,
|
||||
and non-NEB atoms can be listed in the file.
|
||||
|
||||
:line
|
||||
|
||||
Four kinds of output can be generated during a NEB calculation: energy
|
||||
barrier statistics, thermodynamic output by each replica, dump files,
|
||||
and restart files.
|
||||
|
||||
When running with multiple partitions (each of which is a replica in
|
||||
this case), the print-out to the screen and master log.lammps file
|
||||
contains a line of output, printed once every {Nevery} timesteps. It
|
||||
contains the timestep, the maximum force per replica, the maximum
|
||||
force per atom (in any replica), potential gradients in the initial,
|
||||
final, and climbing replicas, the forward and backward energy
|
||||
barriers, the total reaction coordinate (RDT), and the normalized
|
||||
reaction coordinate and potential energy of each replica.
|
||||
|
||||
The "maximum force per replica" is the two-norm of the 3N-length force
|
||||
vector for the atoms in each replica, maximized across replicas, which
|
||||
is what the {ftol} setting is checking against. In this case, N is
|
||||
all the atoms in each replica. The "maximum force per atom" is the
|
||||
maximum force component of any atom in any replica. The potential
|
||||
gradients are the two-norm of the 3N-length force vector solely due to
|
||||
the interaction potential i.e. without adding in inter-replica
|
||||
forces.
|
||||
|
||||
The "reaction coordinate" (RD) for each replica is the two-norm of the
|
||||
3N-length vector of distances between its atoms and the preceding
|
||||
replica's atoms, added to the RD of the preceding replica. The RD of
|
||||
the first replica RD1 = 0.0; the RD of the final replica RDN = RDT,
|
||||
the total reaction coordinate. The normalized RDs are divided by RDT,
|
||||
so that they form a monotonically increasing sequence from zero to
|
||||
one. When computing RD, N only includes the atoms being operated on by
|
||||
the fix neb command.
|
||||
|
||||
The forward (reverse) energy barrier is the potential energy of the
|
||||
highest replica minus the energy of the first (last) replica.
|
||||
|
||||
Supplementary information for all replicas can be printed out to the
|
||||
screen and master log.lammps file by adding the verbose keyword. This
|
||||
information include the following. The "path angle" (pathangle) for
|
||||
the replica i which is the angle between the 3N-length vectors (Ri-1 -
|
||||
Ri) and (Ri+1 - Ri) (where Ri is the atomic coordinates of replica
|
||||
i). A "path angle" of 180 indicates that replicas i-1, i and i+1 are
|
||||
aligned. "angletangrad" is the angle between the 3N-length tangent
|
||||
vector and the 3N-length force vector at image i. The tangent vector
|
||||
is calculated as in "(HenkelmanA)"_#HenkelmanA for all intermediate
|
||||
replicas and at R2 - R1 and RM - RM-1 for the first and last replica,
|
||||
respectively. "anglegrad" is the angle between the 3N-length energy
|
||||
gradient vector of replica i and that of replica i+1. It is not
|
||||
defined for the final replica and reads nan. gradV is the norm of the
|
||||
energy gradient of image i. ReplicaForce is the two-norm of the
|
||||
3N-length force vector (including nudging forces) for replica i.
|
||||
MaxAtomForce is the maximum force component of any atom in replica i.
|
||||
|
||||
When a NEB calculation does not converge properly, the supplementary
|
||||
information can help understanding what is going wrong. For instance
|
||||
when the path angle becomes acute, the definition of tangent used in
|
||||
the NEB calculation is questionable and the NEB cannot may diverge
|
||||
"(Maras)"_#Maras2.
|
||||
|
||||
|
||||
When running on multiple partitions, LAMMPS produces additional log
|
||||
files for each partition, e.g. log.lammps.0, log.lammps.1, etc. For a
|
||||
NEB calculation, these contain the thermodynamic output for each
|
||||
replica.
|
||||
|
||||
If "dump"_dump.html commands in the input script define a filename
|
||||
that includes a {universe} or {uloop} style "variable"_variable.html,
|
||||
then one dump file (per dump command) will be created for each
|
||||
replica. At the end of the NEB calculation, the final snapshot in
|
||||
each file will contain the sequence of snapshots that transition the
|
||||
system over the energy barrier. Earlier snapshots will show the
|
||||
convergence of the replicas to the MEP.
|
||||
|
||||
Likewise, "restart"_restart.html filenames can be specified with a
|
||||
{universe} or {uloop} style "variable"_variable.html, to generate
|
||||
restart files for each replica. These may be useful if the NEB
|
||||
calculation fails to converge properly to the MEP, and you wish to
|
||||
restart the calculation from an intermediate point with altered
|
||||
parameters.
|
||||
|
||||
There are 2 Python scripts provided in the tools/python directory,
|
||||
neb_combine.py and neb_final.py, which are useful in analyzing output
|
||||
from a NEB calculation. Assume a NEB simulation with M replicas, and
|
||||
the NEB atoms labeled with a specific atom type.
|
||||
|
||||
The neb_combine.py script extracts atom coords for the NEB atoms from
|
||||
all M dump files and creates a single dump file where each snapshot
|
||||
contains the NEB atoms from all the replicas and one copy of non-NEB
|
||||
atoms from the first replica (presumed to be identical in other
|
||||
replicas). This can be visualized/animated to see how the NEB atoms
|
||||
relax as the NEB calculation proceeds.
|
||||
|
||||
The neb_final.py script extracts the final snapshot from each of the M
|
||||
dump files to create a single dump file with M snapshots. This can be
|
||||
visualized to watch the system make its transition over the energy
|
||||
barrier.
|
||||
|
||||
To illustrate, here are images from the final snapshot produced by the
|
||||
neb_combine.py script run on the dump files produced by the two
|
||||
example input scripts in examples/neb. Click on them to see a larger
|
||||
image.
|
||||
|
||||
:image(JPG/hop1_small.jpg,JPG/hop1.jpg)
|
||||
:image(JPG/hop2_small.jpg,JPG/hop2.jpg)
|
||||
|
||||
:line
|
||||
|
||||
[Restrictions:]
|
||||
|
||||
This command can only be used if LAMMPS was built with the SPIN
|
||||
package. See the "Build package"_Build_package.html doc
|
||||
page for more info.
|
||||
|
||||
:line
|
||||
|
||||
[Related commands:]
|
||||
|
||||
"min/spin"_min_spin.html, "fix neb/spin"_fix_neb_spin.html
|
||||
|
||||
[Default:]
|
||||
|
||||
none
|
||||
|
||||
:line
|
||||
|
||||
:link(Bessarab)
|
||||
[(Bessarab)] Bessarab, Uzdin, Jonsson, Comp Phys Comm, 196,
|
||||
335-347 (2015).
|
||||
|
|
@ -7,11 +7,10 @@ atom_style spin
|
|||
# necessary for the serial algorithm (sametag)
|
||||
atom_modify map array
|
||||
|
||||
read_data initial.iron_spin
|
||||
|
||||
# setting mass, mag. moments, and interactions for bcc iron
|
||||
# (mass not necessary for fixed lattice calculation)
|
||||
|
||||
read_data initial.iron_spin
|
||||
mass 1 55.845
|
||||
|
||||
pair_style spin/exchange 3.5
|
||||
|
|
|
|||
|
|
@ -7,11 +7,10 @@ atom_style spin
|
|||
# necessary for the serial algorithm (sametag)
|
||||
atom_modify map array
|
||||
|
||||
read_data initial.skyrmion
|
||||
|
||||
# setting mass, mag. moments, and interactions for bcc iron
|
||||
# (mass not necessary for fixed lattice calculation)
|
||||
|
||||
read_data initial.skyrmion
|
||||
mass 1 55.845
|
||||
|
||||
pair_style hybrid/overlay spin/exchange 3.1 spin/dmi 3.1
|
||||
|
|
|
|||
|
|
@ -60,20 +60,20 @@ class FixNEB_spin : public Fix {
|
|||
double **spprev,**spnext,**fmnext;
|
||||
double **springF;
|
||||
double **tangent;
|
||||
double **xsend,**xrecv; // coords to send/recv to/from other replica
|
||||
double **fsend,**frecv; // coords to send/recv to/from other replica
|
||||
double **spsend,**sprecv; // sp to send/recv to/from other replica
|
||||
double **fmsend,**fmrecv; // fm to send/recv to/from other replica
|
||||
tagint *tagsend,*tagrecv; // ditto for atom IDs
|
||||
double **xsend,**xrecv; // coords to send/recv to/from other replica
|
||||
double **fsend,**frecv; // coords to send/recv to/from other replica
|
||||
double **spsend,**sprecv; // sp to send/recv to/from other replica
|
||||
double **fmsend,**fmrecv; // fm to send/recv to/from other replica
|
||||
tagint *tagsend,*tagrecv; // ditto for atom IDs
|
||||
|
||||
// info gathered from all procs in my replica
|
||||
double **xsendall,**xrecvall; // coords to send/recv to/from other replica
|
||||
double **fsendall,**frecvall; // force to send/recv to/from other replica
|
||||
double **spsendall,**sprecvall; // sp to send/recv to/from other replica
|
||||
double **fmsendall,**fmrecvall; // fm to send/recv to/from other replica
|
||||
tagint *tagsendall,*tagrecvall; // ditto for atom IDs
|
||||
// info gathered from all procs in my replica
|
||||
double **xsendall,**xrecvall; // coords to send/recv to/from other replica
|
||||
double **fsendall,**frecvall; // force to send/recv to/from other replica
|
||||
double **spsendall,**sprecvall; // sp to send/recv to/from other replica
|
||||
double **fmsendall,**fmrecvall; // fm to send/recv to/from other replica
|
||||
tagint *tagsendall,*tagrecvall; // ditto for atom IDs
|
||||
|
||||
int *counts,*displacements; // used for MPI_Gather
|
||||
int *counts,*displacements; // used for MPI_Gather
|
||||
|
||||
double geodesic_distance(double *, double *);
|
||||
void inter_replica_comm();
|
||||
|
|
@ -97,16 +97,16 @@ E: Potential energy ID for fix neb does not exist
|
|||
|
||||
Self-explanatory.
|
||||
|
||||
E: Too many active NEB atoms
|
||||
E: Too many active GNEB atoms
|
||||
|
||||
UNDOCUMENTED
|
||||
|
||||
E: Too many atoms for NEB
|
||||
E: Too many atoms for GNEB
|
||||
|
||||
UNDOCUMENTED
|
||||
|
||||
U: Atom count changed in fix neb
|
||||
|
||||
This is not allowed in a NEB calculation.
|
||||
This is not allowed in a GNEB calculation.
|
||||
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -111,6 +111,8 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
|
|||
double **sp = atom->sp;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
int temp_flag,rot_flag;
|
||||
temp_flag = rot_flag = 0;
|
||||
int ii = 0;
|
||||
double spinit[3],spfinal[3];
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
|
|
@ -123,7 +125,7 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
|
|||
spfinal[2] = buf_final[ii+2];
|
||||
|
||||
// interpolate intermediate spin states
|
||||
|
||||
|
||||
if (fraction == 0.0) {
|
||||
sp[i][0] = spinit[0];
|
||||
sp[i][1] = spinit[1];
|
||||
|
|
@ -133,7 +135,8 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
|
|||
sp[i][1] = spfinal[1];
|
||||
sp[i][2] = spfinal[2];
|
||||
} else {
|
||||
initial_rotation(spinit,spfinal,fraction);
|
||||
temp_flag = initial_rotation(spinit,spfinal,fraction);
|
||||
rot_flag = MAX(temp_flag,rot_flag);
|
||||
sp[i][0] = spfinal[0];
|
||||
sp[i][1] = spfinal[1];
|
||||
sp[i][2] = spfinal[2];
|
||||
|
|
@ -141,6 +144,14 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
|
|||
|
||||
ii += 3;
|
||||
}
|
||||
|
||||
// warning message if one or more couples (spi,spf) were aligned
|
||||
// this breaks Rodrigues' formula, and an arbitrary rotation
|
||||
// vector has to be chosen
|
||||
|
||||
if ((rot_flag > 0) && (comm->me == 0))
|
||||
error->warning(FLERR,"arbitrary initial rotation of one or more spin(s)");
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
|
@ -494,6 +505,8 @@ void NEB_spin::readfile(char *file, int flag)
|
|||
|
||||
int ncount = 0;
|
||||
|
||||
int temp_flag,rot_flag;
|
||||
temp_flag = rot_flag = 0;
|
||||
int nread = 0;
|
||||
while (nread < nlines) {
|
||||
nchunk = MIN(nlines-nread,CHUNK);
|
||||
|
|
@ -566,7 +579,8 @@ void NEB_spin::readfile(char *file, int flag)
|
|||
sp[m][1] = spfinal[1];
|
||||
sp[m][2] = spfinal[2];
|
||||
} else {
|
||||
initial_rotation(spinit,spfinal,fraction);
|
||||
temp_flag = initial_rotation(spinit,spfinal,fraction);
|
||||
rot_flag = MAX(temp_flag,rot_flag);
|
||||
sp[m][0] = spfinal[0];
|
||||
sp[m][1] = spfinal[1];
|
||||
sp[m][2] = spfinal[2];
|
||||
|
|
@ -588,6 +602,13 @@ void NEB_spin::readfile(char *file, int flag)
|
|||
nread += nchunk;
|
||||
}
|
||||
|
||||
// warning message if one or more couples (spi,spf) were aligned
|
||||
// this breaks Rodrigues' formula, and an arbitrary rotation
|
||||
// vector has to be chosen
|
||||
|
||||
if ((rot_flag > 0) && (comm->me == 0))
|
||||
error->warning(FLERR,"arbitrary initial rotation of one or more spin(s)");
|
||||
|
||||
// check that all atom IDs in file were found by a proc
|
||||
|
||||
if (flag == 0) {
|
||||
|
|
@ -621,69 +642,24 @@ void NEB_spin::readfile(char *file, int flag)
|
|||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
initial configuration of spin sploc using Rodrigues' formula
|
||||
initial configuration of intermediate spins using Rodrigues' formula
|
||||
interpolates between initial (spi) and final (stored in sploc)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction)
|
||||
int NEB_spin::initial_rotation(double *spi, double *sploc, double fraction)
|
||||
{
|
||||
|
||||
// implementing initial rotation using atan2
|
||||
// this may not be a sufficient routine, need more accurate verifications
|
||||
// no interpolation for initial and final replica
|
||||
|
||||
// interpolation only for intermediate replica
|
||||
if (fraction == 0.0 || fraction == 1.0) return 0;
|
||||
|
||||
if (fraction == 0.0 || fraction == 1.0) return;
|
||||
|
||||
// initial, final and inter ang. values
|
||||
|
||||
//double itheta,iphi,ftheta,fphi,ktheta,kphi;
|
||||
//double spix,spiy,spiz,spfx,spfy,spfz;
|
||||
//double spkx,spky,spkz,iknorm;
|
||||
|
||||
//spix = spi[0];
|
||||
//spiy = spi[1];
|
||||
//spiz = spi[2];
|
||||
|
||||
//spfx = sploc[0];
|
||||
//spfy = sploc[1];
|
||||
//spfz = sploc[2];
|
||||
|
||||
//iphi = itheta = fphi = ftheta = 0.0;
|
||||
|
||||
//iphi = acos(spiz);
|
||||
//if (sin(iphi) != 0.0)
|
||||
// itheta = acos(spix/sin(iphi));
|
||||
|
||||
//fphi = acos(spfz);
|
||||
//if (sin(fphi) != 0.0)
|
||||
// ftheta = acos(spfx/sin(fphi));
|
||||
|
||||
//kphi = iphi + fraction*(fphi-iphi);
|
||||
//ktheta = itheta + fraction*(ftheta-itheta);
|
||||
|
||||
//spkx = cos(ktheta)*sin(kphi);
|
||||
//spky = sin(ktheta)*sin(kphi);
|
||||
//spkz = cos(kphi);
|
||||
|
||||
//double knormsq = spkx*spkx + spky*spky + spkz*spkz;
|
||||
//if (knormsq != 0.0)
|
||||
// iknorm = 1.0/sqrt(knormsq);
|
||||
|
||||
//spkx *= iknorm;
|
||||
//spky *= iknorm;
|
||||
//spkz *= iknorm;
|
||||
|
||||
//sploc[0] = spkx;
|
||||
//sploc[1] = spky;
|
||||
//sploc[2] = spkz;
|
||||
|
||||
int rot_flag = 0;
|
||||
double kx,ky,kz;
|
||||
double spix,spiy,spiz,spfx,spfy,spfz;
|
||||
double kcrossx,kcrossy,kcrossz,knormsq;
|
||||
double kdots;
|
||||
double spkx,spky,spkz;
|
||||
double sdot,omega,iknorm,isnorm;
|
||||
double sidotsf,omega,iknorm,isnorm;
|
||||
|
||||
spix = spi[0];
|
||||
spiy = spi[1];
|
||||
|
|
@ -698,43 +674,73 @@ void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction)
|
|||
kz = spix*spfy - spiy*spfx;
|
||||
|
||||
knormsq = kx*kx+ky*ky+kz*kz;
|
||||
|
||||
if (knormsq != 0.0) {
|
||||
iknorm = 1.0/sqrt(knormsq);
|
||||
kx *= iknorm;
|
||||
ky *= iknorm;
|
||||
kz *= iknorm;
|
||||
sidotsf = spix*spfx + spiy*spfy + spiz*spfz;
|
||||
|
||||
// if knormsq == 0.0, init and final spins are aligned
|
||||
// Rodrigues' formula breaks, needs to define another axis k
|
||||
|
||||
if (knormsq == 0.0) {
|
||||
if (sidotsf > 0.0) { // spins aligned and in same direction
|
||||
return 0;
|
||||
} else if (sidotsf < 0.0) { // spins aligned and in opposite directions
|
||||
|
||||
// defining a rotation axis
|
||||
// first guess, k = spi x [100]
|
||||
// second guess, k = spi x [010]
|
||||
|
||||
if (spiy*spiy + spiz*spiz != 0.0) { // spin not along [100]
|
||||
kx = 0.0;
|
||||
ky = spiz;
|
||||
kz = -spiy;
|
||||
} else if (spix*spix + spiz*spiz != 0.0) { // spin not along [010]
|
||||
kx = -spiz;
|
||||
ky = 0.0;
|
||||
kz = spix;
|
||||
} else error->all(FLERR,"Incorrect initial rotation operation");
|
||||
rot_flag = 1;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
kcrossx = ky*spiz - kz*spiy;
|
||||
kcrossy = kz*spix - kx*spiz;
|
||||
kcrossz = kx*spiy - ky*spix;
|
||||
|
||||
kdots = kx*spix + ky*spiz + kz*spiz;
|
||||
sdot = spix*spfx + spiy*spfy + spiz*spfz;
|
||||
|
||||
omega = acos(sdot);
|
||||
omega = acos(sidotsf);
|
||||
omega *= fraction;
|
||||
|
||||
spkx = spix*cos(omega) + kcrossx*sin(omega);
|
||||
spky = spiy*cos(omega) + kcrossy*sin(omega);
|
||||
spkz = spiz*cos(omega) + kcrossz*sin(omega);
|
||||
// applying Rodrigues' formula
|
||||
|
||||
spkx = spix*cos(omega);
|
||||
spky = spiy*cos(omega);
|
||||
spkz = spiz*cos(omega);
|
||||
|
||||
spkx += kcrossx*sin(omega);
|
||||
spky += kcrossy*sin(omega);
|
||||
spkz += kcrossz*sin(omega);
|
||||
|
||||
spkx += kx*kdots*(1.0-cos(omega));
|
||||
spky += ky*kdots*(1.0-cos(omega));
|
||||
spkz += kz*kdots*(1.0-cos(omega));
|
||||
|
||||
// normalizing resulting spin vector
|
||||
|
||||
isnorm = 1.0/sqrt(spkx*spkx+spky*spky+spkz*spkz);
|
||||
if (isnorm == 0.0)
|
||||
error->all(FLERR,"Incorrect rotation operation");
|
||||
error->all(FLERR,"Incorrect initial rotation operation");
|
||||
|
||||
spkx *= isnorm;
|
||||
spky *= isnorm;
|
||||
spkz *= isnorm;
|
||||
|
||||
// returns rotated spin
|
||||
|
||||
sploc[0] = spkx;
|
||||
sploc[1] = spky;
|
||||
sploc[2] = spkz;
|
||||
|
||||
return rot_flag;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -57,7 +57,7 @@ class NEB_spin : protected Pointers {
|
|||
double *fmaxatomInRepl; // force on an image
|
||||
|
||||
void readfile(char *, int);
|
||||
void initial_rotation(double *, double *, double);
|
||||
int initial_rotation(double *, double *, double);
|
||||
void open(char *);
|
||||
void print_status();
|
||||
};
|
||||
|
|
|
|||
|
|
@ -171,10 +171,11 @@ void PairSpinDmi::init_style()
|
|||
int ifix = 0;
|
||||
while (ifix < modify->nfix) {
|
||||
if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
|
||||
if (strcmp(modify->fix[ifix]->style,"neb/spin") == 0) break;
|
||||
ifix++;
|
||||
}
|
||||
if ((ifix == modify->nfix) && (comm->me == 0))
|
||||
error->warning(FLERR,"Using pair/spin style without nve/spin");
|
||||
error->warning(FLERR,"Using pair/spin style without nve/spin or neb/spin");
|
||||
|
||||
// get the lattice_flag from nve/spin
|
||||
|
||||
|
|
|
|||
|
|
@ -162,7 +162,7 @@ void PairSpinExchange::init_style()
|
|||
ifix++;
|
||||
}
|
||||
if ((ifix == modify->nfix) && (comm->me == 0))
|
||||
error->warning(FLERR,"Using pair/spin style without nve/spin");
|
||||
error->warning(FLERR,"Using pair/spin style without nve/spin or neb/spin");
|
||||
|
||||
// get the lattice_flag from nve/spin
|
||||
|
||||
|
|
|
|||
|
|
@ -164,10 +164,11 @@ void PairSpinMagelec::init_style()
|
|||
int ifix = 0;
|
||||
while (ifix < modify->nfix) {
|
||||
if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
|
||||
if (strcmp(modify->fix[ifix]->style,"neb/spin") == 0) break;
|
||||
ifix++;
|
||||
}
|
||||
if ((ifix == modify->nfix) && (comm->me == 0))
|
||||
error->warning(FLERR,"Using pair/spin style without nve/spin");
|
||||
error->warning(FLERR,"Using pair/spin style without nve/spin or neb/spin");
|
||||
|
||||
// get the lattice_flag from nve/spin
|
||||
|
||||
|
|
|
|||
|
|
@ -171,10 +171,11 @@ void PairSpinNeel::init_style()
|
|||
int ifix = 0;
|
||||
while (ifix < modify->nfix) {
|
||||
if (strcmp(modify->fix[ifix]->style,"nve/spin") == 0) break;
|
||||
if (strcmp(modify->fix[ifix]->style,"neb/spin") == 0) break;
|
||||
ifix++;
|
||||
}
|
||||
if ((ifix == modify->nfix) && (comm->me == 0))
|
||||
error->warning(FLERR,"Using pair/spin style without nve/spin");
|
||||
error->warning(FLERR,"Using pair/spin style without nve/spin or neb/spin");
|
||||
|
||||
// get the lattice_flag from nve/spin
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue