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

Merge branch 'master' into collected-small-changes

This commit is contained in:
Axel Kohlmeyer 2019-10-07 18:54:33 +02:00
parent 6580700943 edda8a4d07
commit 8c063343a1
70 changed files with 3774 additions and 262 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
cmake/Modules/Packages/KOKKOS.cmake
View File

@ -1,15 +1,21 @@
if(PKG_KOKKOS)
set(LAMMPS_LIB_KOKKOS_SRC_DIR ${LAMMPS_LIB_SOURCE_DIR}/kokkos)
set(LAMMPS_LIB_KOKKOS_BIN_DIR ${LAMMPS_LIB_BINARY_DIR}/kokkos)
option(EXTERNAL_KOKKOS "Build against external kokkos library")
if(EXTERNAL_KOKKOS)
find_package(Kokkos REQUIRED)
list(APPEND LAMMPS_LINK_LIBS Kokkos::kokkos)
else()
set(LAMMPS_LIB_KOKKOS_SRC_DIR ${LAMMPS_LIB_SOURCE_DIR}/kokkos)
set(LAMMPS_LIB_KOKKOS_BIN_DIR ${LAMMPS_LIB_BINARY_DIR}/kokkos)
add_subdirectory(${LAMMPS_LIB_KOKKOS_SRC_DIR} ${LAMMPS_LIB_KOKKOS_BIN_DIR})
set(Kokkos_INCLUDE_DIRS ${LAMMPS_LIB_KOKKOS_SRC_DIR}/core/src
${LAMMPS_LIB_KOKKOS_SRC_DIR}/containers/src
${LAMMPS_LIB_KOKKOS_SRC_DIR}/algorithms/src
${LAMMPS_LIB_KOKKOS_BIN_DIR})
include_directories(${Kokkos_INCLUDE_DIRS})
list(APPEND LAMMPS_LINK_LIBS kokkos)
endif()
add_definitions(-DLMP_KOKKOS)
add_subdirectory(${LAMMPS_LIB_KOKKOS_SRC_DIR} ${LAMMPS_LIB_KOKKOS_BIN_DIR})
set(Kokkos_INCLUDE_DIRS ${LAMMPS_LIB_KOKKOS_SRC_DIR}/core/src
${LAMMPS_LIB_KOKKOS_SRC_DIR}/containers/src
${LAMMPS_LIB_KOKKOS_SRC_DIR}/algorithms/src
${LAMMPS_LIB_KOKKOS_BIN_DIR})
include_directories(${Kokkos_INCLUDE_DIRS})
list(APPEND LAMMPS_LINK_LIBS kokkos)
set(KOKKOS_PKG_SOURCES_DIR ${LAMMPS_SOURCE_DIR}/KOKKOS)
set(KOKKOS_PKG_SOURCES ${KOKKOS_PKG_SOURCES_DIR}/kokkos.cpp

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

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 14 KiB

15
doc/src/Eqs/norm_inf.tex Normal file
View File

@ -0,0 +1,15 @@
\documentclass[preview]{standalone}
\usepackage{varwidth}
\usepackage[utf8x]{inputenc}
\usepackage{amsmath, amssymb, graphics, setspace}
\begin{document}
\begin{varwidth}{50in}
\begin{equation}
|| \vec{F} ||_{inf}
= {\rm max}\left(|F_1^1|, |F_1^2|, |F_1^3| \cdots,
|F_N^1|, |F_N^2|, |F_N^3|\right)
\nonumber
\end{equation}
\end{varwidth}
\end{document}

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

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 9.2 KiB

15
doc/src/Eqs/norm_max.tex Normal file
View File

@ -0,0 +1,15 @@
\documentclass[preview]{standalone}
\usepackage{varwidth}
\usepackage[utf8x]{inputenc}
\usepackage{amsmath, amssymb, graphics, setspace}
\begin{document}
\begin{varwidth}{50in}
\begin{equation}
% \left| \left| \vec{F} \right| \right|_2
|| \vec{F} ||_{max}
= {\rm max}\left(||\vec{F}_1||, \cdots, ||\vec{F}_N||\right)
\nonumber
\end{equation}
\end{varwidth}
\end{document}

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

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 5.9 KiB

15
doc/src/Eqs/norm_two.tex Normal file
View File

@ -0,0 +1,15 @@
\documentclass[preview]{standalone}
\usepackage{varwidth}
\usepackage[utf8x]{inputenc}
\usepackage{amsmath, amssymb, graphics, setspace}
\begin{document}
\begin{varwidth}{50in}
\begin{equation}
% \left| \left| \vec{F} \right| \right|_2
|| \vec{F} ||_{2}
= \sqrt{\vec{F}_1+ \cdots + \vec{F}_N}
\nonumber
\end{equation}
\end{varwidth}
\end{document}

16
doc/src/Errors_messages.txt
View File

@ -4791,6 +4791,22 @@ Self-explanatory. :dd
This fix option cannot be used with point particles. :dd
{Fix langevin gjf and respa are not compatible} :dt
Self-explanatory. :dd
{Fix langevin gjf cannot have period equal to dt/2} :dt
If the period is equal to dt/2 then division by zero will happen. :dd
{Fix langevin gjf should come before fix nve} :dt
Self-explanatory. :dd
{Fix langevin gjf with tbias is not yet implemented with kokkos} :dt
This option is not yet available. :dd
{Fix langevin omega is not yet implemented with kokkos} :dt
This option is not yet available. :dd

4
doc/src/Errors_warnings.txt
View File

@ -248,6 +248,10 @@ included one or more of the following: kspace, triclinic, a hybrid
pair style, an eam pair style, or no "single" function for the pair
style. :dd
{Fix langevin gjf using random gaussians is not implemented with kokkos} :dt
This will most likely cause errors in kinetic fluctuations.
{Fix property/atom mol or charge w/out ghost communication} :dt
A model typically needs these properties defined for ghost atoms. :dd

17
doc/src/fix_bond_react.txt
View File

@ -266,7 +266,7 @@ either 'none' or 'charges.' Further details are provided in the
discussion of the 'update_edges' keyword. The fourth optional section
begins with the keyword 'Constraints' and lists additional criteria
that must be satisfied in order for the reaction to occur. Currently,
there is one type of constraint available, as discussed below.
there are two types of constraints available, as discussed below.
A sample map file is given below:
@ -300,14 +300,23 @@ Equivalences :pre
:line
Any number of additional constraints may be specified in the
Constraints section of the map file. Currently there is one type of
additional constraint, of type 'distance', whose syntax is as follows:
Constraints section of the map file. The constraint of type 'distance'
has syntax as follows:
distance {ID1} {ID2} {rmin} {rmax} :pre
where 'distance' is the required keyword, {ID1} and {ID2} are
pre-reaction atom IDs, and these two atoms must be separated by a
distance between {rmin} and {rmax} for the reaction to occur. This
distance between {rmin} and {rmax} for the reaction to occur.
The constraint of type 'angle' has the following syntax:
angle {ID1} {ID2} {ID3} {amin} {amax} :pre
where 'angle' is the required keyword, {ID1}, {ID2} and {ID3} are
pre-reaction atom IDs, and these three atoms must form an angle
between {amin} and {amax} for the reaction to occur (where {ID2} is
the central atom). Angles must be specified in degrees. This
constraint can be used to enforce a certain orientation between
reacting molecules.

58
doc/src/fix_langevin.txt
View File

@ -24,9 +24,10 @@ keyword = {angmom} or {omega} or {scale} or {tally} or {zero} :l
{angmom} value = {no} or factor
{no} = do not thermostat rotational degrees of freedom via the angular momentum
factor = do thermostat rotational degrees of freedom via the angular momentum and apply numeric scale factor as discussed below
{gjf} value = {no} or {yes}
{gjf} value = {no} or {vfull} or {vhalf}
{no} = use standard formulation
{yes} = use Gronbech-Jensen/Farago formulation
{vfull} = use Gronbech-Jensen/Farago formulation
{vhalf} = use 2GJ formulation
{omega} value = {no} or {yes}
{no} = do not thermostat rotational degrees of freedom via the angular velocity
{yes} = do thermostat rotational degrees of freedom via the angular velocity
@ -217,6 +218,10 @@ the particles. As described below, this energy can then be printed
out or added to the potential energy of the system to monitor energy
conservation.
NOTE: this accumulated energy does NOT include kinetic energy removed
by the {zero} flag. LAMMPS will print a warning when both options are
active.
The keyword {zero} can be used to eliminate drift due to the
thermostat. Because the random forces on different atoms are
independent, they do not sum exactly to zero. As a result, this fix
@ -232,29 +237,24 @@ The keyword {gjf} can be used to run the "Gronbech-Jensen/Farago
described in the papers cited below, the purpose of this method is to
enable longer timesteps to be used (up to the numerical stability
limit of the integrator), while still producing the correct Boltzmann
distribution of atom positions. It is implemented within LAMMPS, by
changing how the random force is applied so that it is composed of
the average of two random forces representing half-contributions from
the previous and current time intervals.
distribution of atom positions.
In common with all methods based on Verlet integration, the
discretized velocities generated by this method in conjunction with
velocity-Verlet time integration are not exactly conjugate to the
positions. As a result the temperature (computed from the discretized
velocities) will be systematically lower than the target temperature,
by a small amount which grows with the timestep. Nonetheless, the
distribution of atom positions will still be consistent with the
The current implementation provides the user with the option to output
the velocity in one of two forms: {vfull} or {vhalf}, which replaces
the outdated option {yes}. The {gjf} option {vfull} outputs the on-site
velocity given in "Gronbech-Jensen/Farago"_#Gronbech-Jensen; this velocity
is shown to be systematically lower than the target temperature by a small
amount, which grows quadratically with the timestep.
The {gjf} option {vhalf} outputs the 2GJ half-step velocity given in
"Gronbech Jensen/Gronbech-Jensen"_#2Gronbech-Jensen; this velocity is shown
to not have any linear statistical errors for any stable time step.
An overview of statistically correct Boltzmann and Maxwell-Boltzmann
sampling of true on-site and true half-step velocities is given in
"Gronbech-Jensen_#1Gronbech-Jensen.
Regardless of the choice of output velocity, the sampling of the configurational
distribution of atom positions is the same, and linearly consistent with the
target temperature.
As an example of using the {gjf} keyword, for molecules containing C-H
bonds, configurational properties generated with dt = 2.5 fs and tdamp
= 100 fs are indistinguishable from dt = 0.5 fs. Because the velocity
distribution systematically decreases with increasing timestep, the
method should not be used to generate properties that depend on the
velocity distribution, such as the velocity auto-correlation function
(VACF). In this example, the velocity distribution at dt = 2.5fs
generates an average temperature of 220 K, instead of 300 K.
:line
Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
@ -312,7 +312,10 @@ This fix can ramp its target temperature over multiple runs, using the
This fix is not invoked during "energy minimization"_minimize.html.
[Restrictions:] none
[Restrictions:]
For {gjf} do not choose damp=dt/2. {gjf} is not compatible
with run_style respa.
[Related commands:]
@ -335,5 +338,10 @@ types, tally = no, zero = no, gjf = no.
: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)
(2013); Gronbech-Jensen, Hayre, and Farago, Comp Phys Comm, 185, 524 (2014)
:link(2Gronbech-Jensen)
[(Gronbech-Jensen)] Gronbech Jensen and Gronbech-Jensen, Mol Phys, 117, 2511 (2019)
:link(1Gronbech-Jensen)
[(Gronbech-Jensen)] Gronbech-Jensen, Mol Phys (2019); https://doi.org/10.1080/00268976.2019.1662506

2
doc/src/fix_langevin_spin.txt
View File

@ -50,7 +50,7 @@ As an example:
fix 1 all precession/spin zeeman 0.01 0.0 0.0 1.0
fix 2 all langevin/spin 300.0 0.01 21
fix 3 all nve/spin lattice yes :pre
fix 3 all nve/spin lattice moving :pre
is correct, but defining a force/spin command after the langevin/spin command
would give an error message.

20
doc/src/fix_nve_spin.txt
View File

@ -15,22 +15,26 @@ fix ID group-ID nve/spin keyword values :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
nve/spin = style name of this fix command :l
keyword = {lattice} :l
{lattice} value = {no} or {yes} :pre
{lattice} value = {moving} or {frozen}
moving = integrate both spin and atomic degress of freedom
frozen = integrate spins on a fixed lattice :pre
:ule
[Examples:]
fix 3 all nve/spin lattice yes
fix 1 all nve/spin lattice no :pre
fix 3 all nve/spin lattice moving
fix 1 all nve/spin lattice frozen :pre
[Description:]
Perform a symplectic integration for the spin or spin-lattice system.
The {lattice} keyword defines if the spins are integrated on a lattice
of fixed atoms (lattice = no), or if atoms are moving (lattice = yes).
By default (lattice = yes), a spin-lattice integration is performed.
of fixed atoms (lattice = frozen), or if atoms are moving
(lattice = moving).
The first case corresponds to a spin dynamics calculation, and
the second to a spin-lattice calculation.
By default a spin-lattice integration is performed (lattice = moving).
The {nve/spin} fix applies a Suzuki-Trotter decomposition to
the equations of motion of the spin lattice system, following the scheme:
@ -63,7 +67,9 @@ instead of "array" is also valid.
"atom_style spin"_atom_style.html, "fix nve"_fix_nve.html
[Default:] none
[Default:]
The option default is lattice = moving.
:line

58
doc/src/min_modify.txt
View File

@ -13,11 +13,15 @@ min_modify command :h3
min_modify keyword values ... :pre
one or more keyword/value pairs may be listed :ulb,l
keyword = {dmax} or {line} or {alpha_damp} or {discrete_factor}
keyword = {dmax} or {line} or {norm} or {alpha_damp} or {discrete_factor}
{dmax} value = max
max = maximum distance for line search to move (distance units)
{line} value = {backtrack} or {quadratic} or {forcezero}
backtrack,quadratic,forcezero = style of linesearch to use
{line} value = {backtrack} or {quadratic} or {forcezero} or {spin_cubic} or {spin_none}
backtrack,quadratic,forcezero,spin_cubic,spin_none = style of linesearch to use
{norm} value = {two} or {max}
two = Euclidean two-norm (length of 3N vector)
inf = max force component across all 3-vectors
max = max force norm across all 3-vectors
{alpha_damp} value = damping
damping = fictitious Gilbert damping for spin minimization (adim)
{discrete_factor} value = factor
@ -69,6 +73,28 @@ difference of two large values (energy before and energy after) and
that difference may be smaller than machine epsilon even if atoms
could move in the gradient direction to reduce forces further.
The choice of a norm can be modified for the min styles {cg}, {sd},
{quickmin}, {fire}, {spin}, {spin/cg} and {spin/lbfgs} using
the {norm} keyword.
The default {two} norm computes the 2-norm (Euclidean length) of the
global force vector:
:c,image(Eqs/norm_two.jpg)
The {max} norm computes the length of the 3-vector force
for each atom (2-norm), and takes the maximum value of those across
all atoms
:c,image(Eqs/norm_max.jpg)
The {inf} norm takes the maximum component across the forces of
all atoms in the system:
:c,image(Eqs/norm_inf.jpg)
For the min styles {spin}, {spin/cg} and {spin/lbfgs}, the force
norm is replaced by the spin-torque norm.
Keywords {alpha_damp} and {discrete_factor} only make sense when
a "min_spin"_min_spin.html command is declared.
Keyword {alpha_damp} defines an analog of a magnetic Gilbert
@ -77,10 +103,24 @@ a given magnetic system.
Keyword {discrete_factor} defines a discretization factor for the
adaptive timestep used in the {spin} minimization.
See "min_spin"_min_spin.html for more information about those
quantities.
Default values are {alpha_damp} = 1.0 and {discrete_factor} = 10.0.
quantities.
[Restrictions:] none
The choice of a line search algorithm for the {spin/cg} and
{spin/lbfgs} styles can be specified via the {line} keyword.
The {spin_cubic} and {spin_none} only make sense when one of those
two minimization styles is declared.
The {spin_cubic} performs the line search based on a cubic interpolation
of the energy along the search direction. The {spin_none} keyword
deactivates the line search procedure.
The {spin_none} is a default value for {line} keyword for both {spin/lbfgs}
and {spin/cg}. Convergence of {spin/lbfgs} can be more robust if
{spin_cubic} line search is used.
[Restrictions:]
For magnetic GNEB calculations, only {spin_none} value for {line} keyword can be used
when styles {spin/cg} and {spin/lbfgs} are employed.
See "neb/spin"_neb_spin.html for more explanation.
[Related commands:]
@ -88,4 +128,8 @@ Default values are {alpha_damp} = 1.0 and {discrete_factor} = 10.0.
[Default:]
The option defaults are dmax = 0.1 and line = quadratic.
The option defaults are dmax = 0.1, line = quadratic and norm = two.
For the {spin}, {spin/cg} and {spin/lbfgs} styles, the
option defaults are alpha_damp = 1.0, discrete_factor = 10.0,
line = spin_none, and norm = euclidean.

53
doc/src/min_spin.txt
View File

@ -6,14 +6,19 @@
:line
min_style spin command :h3
min_style spin/cg command :h3
min_style spin/lbfgs command :h3
[Syntax:]
min_style spin :pre
min_style spin
min_style spin/cg
min_style spin/lbfgs :pre
[Examples:]
min_style spin :pre
min_style spin/lbfgs
min_modify line spin_cubic discrete_factor 10.0 :pre
[Description:]
@ -46,10 +51,39 @@ definition of this timestep.
{discrete_factor} can be defined with the "min_modify"_min_modify.html
command.
NOTE: The {spin} style replaces the force tolerance by a torque
tolerance. See "minimize"_minimize.html for more explanation.
Style {spin/cg} defines an orthogonal spin optimization
(OSO) combined to a conjugate gradient (CG) algorithm.
The "min_modify"_min_modify.html command can be used to
couple the {spin/cg} to a line search procedure, and to modify the
discretization factor {discrete_factor}.
By default, style {spin/cg} does not employ the line search procedure
and uses the adaptive time-step technique in the same way as style {spin}.
[Restrictions:]
Style {spin/lbfgs} defines an orthogonal spin optimization
(OSO) combined to a limited-memory Broyden-Fletcher-Goldfarb-Shanno
(L-BFGS) algorithm.
By default, style {spin/lbfgs} does not employ line search procedure.
If the line search procedure is not used then the discrete factor defines
the maximum root mean squared rotation angle of spins by equation {pi/(5*Kappa)}.
The default value for Kappa is 10.
The {spin_cubic} line search can improve the convergence of the
{spin/lbfgs} algorithm.
The "min_modify"_min_modify.html command can be used to
activate the line search procedure, and to modify the
discretization factor {discrete_factor}.
For more information about styles {spin/cg} and {spin/lbfgs},
see their implementation reported in "(Ivanov)"_#Ivanov1.
NOTE: All the {spin} styles replace the force tolerance by a torque
tolerance. See "minimize"_minimize.html for more explanation.
NOTE: The {spin/cg} and {spin/lbfgs} styles can be used
for magnetic NEB calculations only if the line search procedure
is deactivated. See "neb/spin"_neb_spin.html for more explanation.
[Restrictions:]
This minimization procedure is only applied to spin degrees of
freedom for a frozen lattice configuration.
@ -61,5 +95,10 @@ freedom for a frozen lattice configuration.
[Default:]
The option defaults are {alpha_damp} = 1.0 and {discrete_factor} =
10.0.
The option defaults are {alpha_damp} = 1.0, {discrete_factor} =
10.0, {line} = spin_none and {norm} = euclidean.
:line
:link(Ivanov1)
[(Ivanov)] Ivanov, Uzdin, Jonsson. arXiv preprint arXiv:1904.02669, (2019).

18
doc/src/min_style.txt
View File

@ -11,7 +11,7 @@ min_style command :h3
min_style style :pre
style = {cg} or {cg/kk} or {hftn} or {sd} or {quickmin} or {fire} or {spin} :ul
style = {cg} or {hftn} or {sd} or {quickmin} or {fire} or {spin} or {spin/cg} or {spin/lbfgs} :ul
[Examples:]
@ -64,11 +64,25 @@ a minimization.
Style {spin} is a damped spin dynamics with an adaptive
timestep.
See the "min/spin"_min_spin.html doc page for more information.
Style {spin/cg} uses an orthogonal spin optimization (OSO)
combined to a conjugate gradient (CG) approach to minimize spin
configurations.
Style {spin/lbfgs} uses an orthogonal spin optimization (OSO)
combined to a limited-memory Broyden-Fletcher-Goldfarb-Shanno
(LBFGS) approach to minimize spin configurations.
See the "min/spin"_min_spin.html doc page for more information
about the {spin}, {spin/cg} and {spin/lbfgs} styles.
Either the {quickmin} and {fire} styles are useful in the context of
nudged elastic band (NEB) calculations via the "neb"_neb.html command.
Either the {spin}, {spin/cg} and {spin/lbfgs} styles are useful
in the context of magnetic geodesic nudged elastic band (GNEB) calculations
via the "neb/spin"_neb_spin.html command.
NOTE: The damped dynamic minimizers use whatever timestep you have
defined via the "timestep"_timestep.html command. Often they will
converge more quickly if you use a timestep about 10x larger than you

11
doc/src/minimize.txt
View File

@ -104,12 +104,13 @@ the line search fails because the step distance backtracks to 0.0
the number of outer iterations or timesteps exceeds {maxiter}
the number of total force evaluations exceeds {maxeval} :ul
NOTE: the "minimization style"_min_style.html {spin} replaces
NOTE: the "minimization style"_min_style.html {spin},
{spin/cg}, and {spin/lbfgs} replace
the force tolerance {ftol} by a torque tolerance.
The minimization procedure stops if the 2-norm (length) of the
global torque vector (defined as the cross product between the
spins and their precession vectors omega) is less than {ftol},
or if any of the other criteria are met.
The minimization procedure stops if the 2-norm (length) of the torque vector on atom
(defined as the cross product between the
atomic spin and its precession vectors omega) is less than {ftol},
or if any of the other criteria are met. Torque have the same units as the energy.
NOTE: You can also use the "fix halt"_fix_halt.html command to specify
a general criterion for exiting a minimization, that is a calculation

18
doc/src/neb_spin.txt
View File

@ -59,9 +59,9 @@ 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.
NOTE: As explained below, a GNEB calculation performs a
minimization across all the replicas. One of the "spin"_min_spin.html
style minimizers 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.
@ -170,10 +170,11 @@ 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 spin dynamics
"min_style"_min_style.html, using the {spin} style (see
"min_spin"_min_spin.html for more information).
minimization procedure. To enable
this, you must first define a
"min_style"_min_style.html, using either the {spin},
{spin/cg}, or {spin/lbfgs} style (see
"min_spin"_min_spin.html for more information).
The other styles cannot be used, since they relax the lattice
degrees of freedom instead of the spins.
@ -358,6 +359,9 @@ 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.
For magnetic GNEB calculations, only {spin_none} value for {line} keyword can be used
when styles {spin/cg} and {spin/lbfgs} are employed.
:line
[Related commands:]

5
doc/src/pair_spin_dipole.txt
View File

@ -25,9 +25,8 @@ pair_coeff * * 10.0
pair_coeff 2 3 8.0 :pre
pair_style spin/dipole/long 9.0
pair_coeff * * 1.0 1.0
pair_coeff 2 3 1.0 1.0 2.5 4.0 scale 0.5
pair_coeff 2 3 1.0 1.0 2.5 4.0 :pre
pair_coeff * * 10.0
pair_coeff 2 3 6.0 :pre
[Description:]

6
doc/utils/sphinx-config/false_positives.txt
View File

@ -275,6 +275,7 @@ Broadwell
Broglie
brownian
brownw
Broyden
Bryantsev
Btarget
btype
@ -988,6 +989,7 @@ gmask
Gmask
gneb
GNEB
Goldfarb
googlemail
Gordan
GPa
@ -1405,6 +1407,7 @@ Laupretre
lavenderblush
lawngreen
lB
lbfgs
lbl
LBtype
lcbop
@ -2041,6 +2044,7 @@ Orsi
ortho
orthonormal
orthorhombic
oso
ot
Otype
Ouldridge
@ -2272,6 +2276,7 @@ qoffload
qopenmp
qoverride
qtb
quadratically
quadrupolar
Quant
quartic
@ -2513,6 +2518,7 @@ setvel
sfftw
Sg
Shan
Shanno
shapex
shapey
shapez

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

@ -32,7 +32,7 @@ neigh_modify every 10 check yes delay 20
fix 1 all precession/spin anisotropy 0.0000033 0.0 0.0 1.0
fix 2 all langevin/spin 0.0 0.1 21
fix 3 all nve/spin lattice no
fix 3 all nve/spin lattice frozen
timestep 0.0002

2
examples/SPIN/cobalt_fcc/in.spin.cobalt_fcc
View File

@ -35,7 +35,7 @@ fix_modify 1 energy yes
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice moving
timestep 0.0001
# compute and output options

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

@ -37,7 +37,7 @@ neigh_modify every 10 check yes delay 20
fix 1 all precession/spin anisotropy 0.01 0.0 0.0 1.0
#fix 2 all langevin/spin 0.0 0.0 21
fix 2 all langevin/spin 0.0 0.1 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice moving
timestep 0.0001

2
examples/SPIN/dipole_spin/in.spin.iron_dipole_cut
View File

@ -33,7 +33,7 @@ fix 1 all precession/spin cubic 0.001 0.0005 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1
fix_modify 1 energy yes
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice moving
timestep 0.0001
# compute and output options

2
examples/SPIN/dipole_spin/in.spin.iron_dipole_ewald
View File

@ -35,7 +35,7 @@ fix 1 all precession/spin cubic 0.001 0.0005 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1
fix_modify 1 energy yes
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice moving
timestep 0.0001
# compute and output options

2
examples/SPIN/dipole_spin/in.spin.iron_dipole_pppm
View File

@ -36,7 +36,7 @@ fix 1 all precession/spin cubic 0.001 0.0005 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1
fix_modify 1 energy yes
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice moving
timestep 0.0001
# compute and output options

2
examples/SPIN/iron/in.spin.iron
View File

@ -33,7 +33,7 @@ neigh_modify every 10 check yes delay 20
fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice moving
timestep 0.0001
# compute and output options

2
examples/SPIN/iron/in.spin.iron_cubic
View File

@ -31,7 +31,7 @@ fix 1 all precession/spin cubic 0.001 0.0005 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1
fix_modify 1 energy yes
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice moving
timestep 0.0001
# compute and output options

2
examples/SPIN/nickel/in.spin.nickel
View File

@ -33,7 +33,7 @@ neigh_modify every 10 check yes delay 20
fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice moving
timestep 0.0001
# compute and output options

2
examples/SPIN/nickel/in.spin.nickel_cubic
View File

@ -35,7 +35,7 @@ fix 1 all precession/spin cubic -0.0001 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.
fix_modify 1 energy yes
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice moving
timestep 0.0001
# compute and output options

2
examples/SPIN/read_restart/in.spin.read_data
View File

@ -20,7 +20,7 @@ neigh_modify every 1 check no delay 0
fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice moving
timestep 0.0001
# define outputs and computes

2
examples/SPIN/read_restart/in.spin.restart
View File

@ -24,7 +24,7 @@ neigh_modify every 1 check no delay 0
fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
fix 3 all nve/spin lattice moving
timestep 0.0001
# define outputs

2
examples/SPIN/read_restart/in.spin.write_restart
View File

@ -29,7 +29,7 @@ neigh_modify every 10 check yes delay 20
fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0
fix 2 all langevin/spin 100.0 0.01 21
fix 3 all nve/spin lattice no
fix 3 all nve/spin lattice frozen
timestep 0.0001
# compute and output options

2
examples/SPIN/setforce_spin/in.spinmin.setforce
View File

@ -35,7 +35,7 @@ fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0 anisotropy 5e-05 0.0 0.0 1.0
fix_modify 1 energy yes
fix 2 fixed_spin setforce/spin 0.0 0.0 0.0
fix 3 all langevin/spin 0.0 0.1 21
fix 4 all nve/spin lattice no
fix 4 all nve/spin lattice frozen
timestep 0.0001

54
examples/SPIN/spinmin/in.spinmin_cg.bfo Normal file
View File

@ -0,0 +1,54 @@
# bfo in a 3d periodic box
units metal
dimension 3
boundary p p f
atom_style spin
# necessary for the serial algorithm (sametag)
atom_modify map array
lattice sc 3.96
region box block 0.0 34.0 0.0 34.0 0.0 1.0
create_box 1 box
create_atoms 1 box
# setting mass, mag. moments, and interactions for bcc iron
mass 1 1.0
set group all spin/random 11 2.50
pair_style hybrid/overlay spin/exchange 6.0 spin/magelec 4.5 spin/dmi 4.5
pair_coeff * * spin/exchange exchange 6.0 -0.01575 0.0 1.965
# pair_coeff * * spin/magelec magelec 4.5 0.000109 1.0 1.0 1.0
pair_coeff * * spin/magelec magelec 4.5 0.00109 1.0 1.0 1.0
pair_coeff * * spin/dmi dmi 4.5 0.00005 1.0 1.0 1.0
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
fix 1 all precession/spin anisotropy 0.0000033 0.0 0.0 1.0
fix_modify 1 energy yes
timestep 0.0001
compute out_mag all spin
compute out_pe all pe
compute out_ke all ke
compute out_temp all temp
variable magz equal c_out_mag[3]
variable magnorm equal c_out_mag[4]
variable emag equal c_out_mag[5]
variable tmag equal c_out_mag[6]
thermo 100
thermo_style custom step time v_magnorm v_emag v_tmag temp etotal
thermo_modify format float %20.15g
compute outsp all property/atom spx spy spz sp fmx fmy fmz
dump 1 all custom 50 dump.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3] c_outsp[4] c_outsp[5] c_outsp[6] c_outsp[7]
min_style spin/cg
# min_modify line spin_none discrete_factor 10.0
minimize 1.0e-10 1.0e-10 10000 10000

55
examples/SPIN/spinmin/in.spinmin_lbfgs.bfo Normal file
View File

@ -0,0 +1,55 @@
# bfo in a 3d periodic box
units metal
dimension 3
boundary p p f
atom_style spin
# necessary for the serial algorithm (sametag)
atom_modify map array
lattice sc 3.96
region box block 0.0 34.0 0.0 34.0 0.0 1.0
create_box 1 box
create_atoms 1 box
# setting mass, mag. moments, and interactions for bcc iron
mass 1 1.0
set group all spin/random 11 2.50
pair_style hybrid/overlay spin/exchange 6.0 spin/magelec 4.5 spin/dmi 4.5
pair_coeff * * spin/exchange exchange 6.0 -0.01575 0.0 1.965
#pair_coeff * * spin/magelec magelec 4.5 0.000109 1.0 1.0 1.0
pair_coeff * * spin/magelec magelec 4.5 0.00109 1.0 1.0 1.0
pair_coeff * * spin/dmi dmi 4.5 0.00005 1.0 1.0 1.0
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
fix 1 all precession/spin anisotropy 0.0000033 0.0 0.0 1.0
fix_modify 1 energy yes
timestep 0.0001
compute out_mag all spin
compute out_pe all pe
compute out_ke all ke
compute out_temp all temp
variable magz equal c_out_mag[3]
variable magnorm equal c_out_mag[4]
variable emag equal c_out_mag[5]
variable tmag equal c_out_mag[6]
thermo 50
thermo_style custom step time v_magnorm v_emag v_tmag temp etotal
thermo_modify format float %20.15g
compute outsp all property/atom spx spy spz sp fmx fmy fmz
dump 1 all custom 50 dump.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3] c_outsp[4] c_outsp[5] c_outsp[6] c_outsp[7]
min_style spin/lbfgs
# min_modify line spin_cubic discrete_factor 10.0
min_modify norm max
minimize 1.0e-15 1.0e-10 10000 1000

13
examples/gjf/README.md Normal file
View File

@ -0,0 +1,13 @@
# LAMMPS GJF-2GJ THERMOSTAT EXAMPLE
## GJF-2GJ THERMOSTAT
This directory contains the ingredients to run an NVT simulation using the GJF-2GJ thermostat.
Example:
```
NP=4 #number of processors
mpirun -np $NP lmp_mpi -in.gjf.vhalf
```
## Required LAMMPS packages: MOLECULE package

886
examples/gjf/argon.lmp Normal file
View File

@ -0,0 +1,886 @@
LAMMPS description
864 atoms
0 bonds
0 angles
0 dihedrals
0 impropers
1 atom types
0 bond types
0 angle types
0 dihedral types
0 improper types
0.0000000 32.146000 xlo xhi
0.0000000 32.146000 ylo yhi
0.0000000 32.146000 zlo zhi
Atoms
1 1 1 0.0000000 0.0000000 2.6790000 2.6790000
2 2 1 0.0000000 0.0000000 2.6790000 8.0360000
3 3 1 0.0000000 0.0000000 2.6790000 13.3940000
4 4 1 0.0000000 0.0000000 2.6790000 18.7520000
5 5 1 0.0000000 0.0000000 2.6790000 24.1090000
6 6 1 0.0000000 0.0000000 2.6790000 29.4670000
7 7 1 0.0000000 0.0000000 8.0360000 2.6790000
8 8 1 0.0000000 0.0000000 8.0360000 8.0360000
9 9 1 0.0000000 0.0000000 8.0360000 13.3940000
10 10 1 0.0000000 0.0000000 8.0360000 18.7520000
11 11 1 0.0000000 0.0000000 8.0360000 24.1090000
12 12 1 0.0000000 0.0000000 8.0360000 29.4670000
13 13 1 0.0000000 0.0000000 13.3940000 2.6790000
14 14 1 0.0000000 0.0000000 13.3940000 8.0360000
15 15 1 0.0000000 0.0000000 13.3940000 13.3940000
16 16 1 0.0000000 0.0000000 13.3940000 18.7520000
17 17 1 0.0000000 0.0000000 13.3940000 24.1090000
18 18 1 0.0000000 0.0000000 13.3940000 29.4670000
19 19 1 0.0000000 0.0000000 18.7520000 2.6790000
20 20 1 0.0000000 0.0000000 18.7520000 8.0360000
21 21 1 0.0000000 0.0000000 18.7520000 13.3940000
22 22 1 0.0000000 0.0000000 18.7520000 18.7520000
23 23 1 0.0000000 0.0000000 18.7520000 24.1090000
24 24 1 0.0000000 0.0000000 18.7520000 29.4670000
25 25 1 0.0000000 0.0000000 24.1090000 2.6790000
26 26 1 0.0000000 0.0000000 24.1090000 8.0360000
27 27 1 0.0000000 0.0000000 24.1090000 13.3940000
28 28 1 0.0000000 0.0000000 24.1090000 18.7520000
29 29 1 0.0000000 0.0000000 24.1090000 24.1090000
30 30 1 0.0000000 0.0000000 24.1090000 29.4670000
31 31 1 0.0000000 0.0000000 29.4670000 2.6790000
32 32 1 0.0000000 0.0000000 29.4670000 8.0360000
33 33 1 0.0000000 0.0000000 29.4670000 13.3940000
34 34 1 0.0000000 0.0000000 29.4670000 18.7520000
35 35 1 0.0000000 0.0000000 29.4670000 24.1090000
36 36 1 0.0000000 0.0000000 29.4670000 29.4670000
37 37 1 0.0000000 5.3580000 2.6790000 2.6790000
38 38 1 0.0000000 5.3580000 2.6790000 8.0360000
39 39 1 0.0000000 5.3580000 2.6790000 13.3940000
40 40 1 0.0000000 5.3580000 2.6790000 18.7520000
41 41 1 0.0000000 5.3580000 2.6790000 24.1090000
42 42 1 0.0000000 5.3580000 2.6790000 29.4670000
43 43 1 0.0000000 5.3580000 8.0360000 2.6790000
44 44 1 0.0000000 5.3580000 8.0360000 8.0360000
45 45 1 0.0000000 5.3580000 8.0360000 13.3940000
46 46 1 0.0000000 5.3580000 8.0360000 18.7520000
47 47 1 0.0000000 5.3580000 8.0360000 24.1090000
48 48 1 0.0000000 5.3580000 8.0360000 29.4670000
49 49 1 0.0000000 5.3580000 13.3940000 2.6790000
50 50 1 0.0000000 5.3580000 13.3940000 8.0360000
51 51 1 0.0000000 5.3580000 13.3940000 13.3940000
52 52 1 0.0000000 5.3580000 13.3940000 18.7520000
53 53 1 0.0000000 5.3580000 13.3940000 24.1090000
54 54 1 0.0000000 5.3580000 13.3940000 29.4670000
55 55 1 0.0000000 5.3580000 18.7520000 2.6790000
56 56 1 0.0000000 5.3580000 18.7520000 8.0360000
57 57 1 0.0000000 5.3580000 18.7520000 13.3940000
58 58 1 0.0000000 5.3580000 18.7520000 18.7520000
59 59 1 0.0000000 5.3580000 18.7520000 24.1090000
60 60 1 0.0000000 5.3580000 18.7520000 29.4670000
61 61 1 0.0000000 5.3580000 24.1090000 2.6790000
62 62 1 0.0000000 5.3580000 24.1090000 8.0360000
63 63 1 0.0000000 5.3580000 24.1090000 13.3940000
64 64 1 0.0000000 5.3580000 24.1090000 18.7520000
65 65 1 0.0000000 5.3580000 24.1090000 24.1090000
66 66 1 0.0000000 5.3580000 24.1090000 29.4670000
67 67 1 0.0000000 5.3580000 29.4670000 2.6790000
68 68 1 0.0000000 5.3580000 29.4670000 8.0360000
69 69 1 0.0000000 5.3580000 29.4670000 13.3940000
70 70 1 0.0000000 5.3580000 29.4670000 18.7520000
71 71 1 0.0000000 5.3580000 29.4670000 24.1090000
72 72 1 0.0000000 5.3580000 29.4670000 29.4670000
73 73 1 0.0000000 10.7150000 2.6790000 2.6790000
74 74 1 0.0000000 10.7150000 2.6790000 8.0360000
75 75 1 0.0000000 10.7150000 2.6790000 13.3940000
76 76 1 0.0000000 10.7150000 2.6790000 18.7520000
77 77 1 0.0000000 10.7150000 2.6790000 24.1090000
78 78 1 0.0000000 10.7150000 2.6790000 29.4670000
79 79 1 0.0000000 10.7150000 8.0360000 2.6790000
80 80 1 0.0000000 10.7150000 8.0360000 8.0360000
81 81 1 0.0000000 10.7150000 8.0360000 13.3940000
82 82 1 0.0000000 10.7150000 8.0360000 18.7520000
83 83 1 0.0000000 10.7150000 8.0360000 24.1090000
84 84 1 0.0000000 10.7150000 8.0360000 29.4670000
85 85 1 0.0000000 10.7150000 13.3940000 2.6790000
86 86 1 0.0000000 10.7150000 13.3940000 8.0360000
87 87 1 0.0000000 10.7150000 13.3940000 13.3940000
88 88 1 0.0000000 10.7150000 13.3940000 18.7520000
89 89 1 0.0000000 10.7150000 13.3940000 24.1090000
90 90 1 0.0000000 10.7150000 13.3940000 29.4670000
91 91 1 0.0000000 10.7150000 18.7520000 2.6790000
92 92 1 0.0000000 10.7150000 18.7520000 8.0360000
93 93 1 0.0000000 10.7150000 18.7520000 13.3940000
94 94 1 0.0000000 10.7150000 18.7520000 18.7520000
95 95 1 0.0000000 10.7150000 18.7520000 24.1090000
96 96 1 0.0000000 10.7150000 18.7520000 29.4670000
97 97 1 0.0000000 10.7150000 24.1090000 2.6790000
98 98 1 0.0000000 10.7150000 24.1090000 8.0360000
99 99 1 0.0000000 10.7150000 24.1090000 13.3940000
100 100 1 0.0000000 10.7150000 24.1090000 18.7520000
101 101 1 0.0000000 10.7150000 24.1090000 24.1090000
102 102 1 0.0000000 10.7150000 24.1090000 29.4670000
103 103 1 0.0000000 10.7150000 29.4670000 2.6790000
104 104 1 0.0000000 10.7150000 29.4670000 8.0360000
105 105 1 0.0000000 10.7150000 29.4670000 13.3940000
106 106 1 0.0000000 10.7150000 29.4670000 18.7520000
107 107 1 0.0000000 10.7150000 29.4670000 24.1090000
108 108 1 0.0000000 10.7150000 29.4670000 29.4670000
109 109 1 0.0000000 16.0730000 2.6790000 2.6790000
110 110 1 0.0000000 16.0730000 2.6790000 8.0360000
111 111 1 0.0000000 16.0730000 2.6790000 13.3940000
112 112 1 0.0000000 16.0730000 2.6790000 18.7520000
113 113 1 0.0000000 16.0730000 2.6790000 24.1090000
114 114 1 0.0000000 16.0730000 2.6790000 29.4670000
115 115 1 0.0000000 16.0730000 8.0360000 2.6790000
116 116 1 0.0000000 16.0730000 8.0360000 8.0360000
117 117 1 0.0000000 16.0730000 8.0360000 13.3940000
118 118 1 0.0000000 16.0730000 8.0360000 18.7520000
119 119 1 0.0000000 16.0730000 8.0360000 24.1090000
120 120 1 0.0000000 16.0730000 8.0360000 29.4670000
121 121 1 0.0000000 16.0730000 13.3940000 2.6790000
122 122 1 0.0000000 16.0730000 13.3940000 8.0360000
123 123 1 0.0000000 16.0730000 13.3940000 13.3940000
124 124 1 0.0000000 16.0730000 13.3940000 18.7520000
125 125 1 0.0000000 16.0730000 13.3940000 24.1090000
126 126 1 0.0000000 16.0730000 13.3940000 29.4670000
127 127 1 0.0000000 16.0730000 18.7520000 2.6790000
128 128 1 0.0000000 16.0730000 18.7520000 8.0360000
129 129 1 0.0000000 16.0730000 18.7520000 13.3940000
130 130 1 0.0000000 16.0730000 18.7520000 18.7520000
131 131 1 0.0000000 16.0730000 18.7520000 24.1090000
132 132 1 0.0000000 16.0730000 18.7520000 29.4670000
133 133 1 0.0000000 16.0730000 24.1090000 2.6790000
134 134 1 0.0000000 16.0730000 24.1090000 8.0360000
135 135 1 0.0000000 16.0730000 24.1090000 13.3940000
136 136 1 0.0000000 16.0730000 24.1090000 18.7520000
137 137 1 0.0000000 16.0730000 24.1090000 24.1090000
138 138 1 0.0000000 16.0730000 24.1090000 29.4670000
139 139 1 0.0000000 16.0730000 29.4670000 2.6790000
140 140 1 0.0000000 16.0730000 29.4670000 8.0360000
141 141 1 0.0000000 16.0730000 29.4670000 13.3940000
142 142 1 0.0000000 16.0730000 29.4670000 18.7520000
143 143 1 0.0000000 16.0730000 29.4670000 24.1090000
144 144 1 0.0000000 16.0730000 29.4670000 29.4670000
145 145 1 0.0000000 21.4310000 2.6790000 2.6790000
146 146 1 0.0000000 21.4310000 2.6790000 8.0360000
147 147 1 0.0000000 21.4310000 2.6790000 13.3940000
148 148 1 0.0000000 21.4310000 2.6790000 18.7520000
149 149 1 0.0000000 21.4310000 2.6790000 24.1090000
150 150 1 0.0000000 21.4310000 2.6790000 29.4670000
151 151 1 0.0000000 21.4310000 8.0360000 2.6790000
152 152 1 0.0000000 21.4310000 8.0360000 8.0360000
153 153 1 0.0000000 21.4310000 8.0360000 13.3940000
154 154 1 0.0000000 21.4310000 8.0360000 18.7520000
155 155 1 0.0000000 21.4310000 8.0360000 24.1090000
156 156 1 0.0000000 21.4310000 8.0360000 29.4670000
157 157 1 0.0000000 21.4310000 13.3940000 2.6790000
158 158 1 0.0000000 21.4310000 13.3940000 8.0360000
159 159 1 0.0000000 21.4310000 13.3940000 13.3940000
160 160 1 0.0000000 21.4310000 13.3940000 18.7520000
161 161 1 0.0000000 21.4310000 13.3940000 24.1090000
162 162 1 0.0000000 21.4310000 13.3940000 29.4670000
163 163 1 0.0000000 21.4310000 18.7520000 2.6790000
164 164 1 0.0000000 21.4310000 18.7520000 8.0360000
165 165 1 0.0000000 21.4310000 18.7520000 13.3940000
166 166 1 0.0000000 21.4310000 18.7520000 18.7520000
167 167 1 0.0000000 21.4310000 18.7520000 24.1090000
168 168 1 0.0000000 21.4310000 18.7520000 29.4670000
169 169 1 0.0000000 21.4310000 24.1090000 2.6790000
170 170 1 0.0000000 21.4310000 24.1090000 8.0360000
171 171 1 0.0000000 21.4310000 24.1090000 13.3940000
172 172 1 0.0000000 21.4310000 24.1090000 18.7520000
173 173 1 0.0000000 21.4310000 24.1090000 24.1090000
174 174 1 0.0000000 21.4310000 24.1090000 29.4670000
175 175 1 0.0000000 21.4310000 29.4670000 2.6790000
176 176 1 0.0000000 21.4310000 29.4670000 8.0360000
177 177 1 0.0000000 21.4310000 29.4670000 13.3940000
178 178 1 0.0000000 21.4310000 29.4670000 18.7520000
179 179 1 0.0000000 21.4310000 29.4670000 24.1090000
180 180 1 0.0000000 21.4310000 29.4670000 29.4670000
181 181 1 0.0000000 26.7880000 2.6790000 2.6790000
182 182 1 0.0000000 26.7880000 2.6790000 8.0360000
183 183 1 0.0000000 26.7880000 2.6790000 13.3940000
184 184 1 0.0000000 26.7880000 2.6790000 18.7520000
185 185 1 0.0000000 26.7880000 2.6790000 24.1090000
186 186 1 0.0000000 26.7880000 2.6790000 29.4670000
187 187 1 0.0000000 26.7880000 8.0360000 2.6790000
188 188 1 0.0000000 26.7880000 8.0360000 8.0360000
189 189 1 0.0000000 26.7880000 8.0360000 13.3940000
190 190 1 0.0000000 26.7880000 8.0360000 18.7520000
191 191 1 0.0000000 26.7880000 8.0360000 24.1090000
192 192 1 0.0000000 26.7880000 8.0360000 29.4670000
193 193 1 0.0000000 26.7880000 13.3940000 2.6790000
194 194 1 0.0000000 26.7880000 13.3940000 8.0360000
195 195 1 0.0000000 26.7880000 13.3940000 13.3940000
196 196 1 0.0000000 26.7880000 13.3940000 18.7520000
197 197 1 0.0000000 26.7880000 13.3940000 24.1090000
198 198 1 0.0000000 26.7880000 13.3940000 29.4670000
199 199 1 0.0000000 26.7880000 18.7520000 2.6790000
200 200 1 0.0000000 26.7880000 18.7520000 8.0360000
201 201 1 0.0000000 26.7880000 18.7520000 13.3940000
202 202 1 0.0000000 26.7880000 18.7520000 18.7520000
203 203 1 0.0000000 26.7880000 18.7520000 24.1090000
204 204 1 0.0000000 26.7880000 18.7520000 29.4670000
205 205 1 0.0000000 26.7880000 24.1090000 2.6790000
206 206 1 0.0000000 26.7880000 24.1090000 8.0360000
207 207 1 0.0000000 26.7880000 24.1090000 13.3940000
208 208 1 0.0000000 26.7880000 24.1090000 18.7520000
209 209 1 0.0000000 26.7880000 24.1090000 24.1090000
210 210 1 0.0000000 26.7880000 24.1090000 29.4670000
211 211 1 0.0000000 26.7880000 29.4670000 2.6790000
212 212 1 0.0000000 26.7880000 29.4670000 8.0360000
213 213 1 0.0000000 26.7880000 29.4670000 13.3940000
214 214 1 0.0000000 26.7880000 29.4670000 18.7520000
215 215 1 0.0000000 26.7880000 29.4670000 24.1090000
216 216 1 0.0000000 26.7880000 29.4670000 29.4670000
217 217 1 0.0000000 2.6790000 5.3580000 2.6790000
218 218 1 0.0000000 2.6790000 5.3580000 8.0360000
219 219 1 0.0000000 2.6790000 5.3580000 13.3940000
220 220 1 0.0000000 2.6790000 5.3580000 18.7520000
221 221 1 0.0000000 2.6790000 5.3580000 24.1090000
222 222 1 0.0000000 2.6790000 5.3580000 29.4670000
223 223 1 0.0000000 2.6790000 10.7150000 2.6790000
224 224 1 0.0000000 2.6790000 10.7150000 8.0360000
225 225 1 0.0000000 2.6790000 10.7150000 13.3940000
226 226 1 0.0000000 2.6790000 10.7150000 18.7520000
227 227 1 0.0000000 2.6790000 10.7150000 24.1090000
228 228 1 0.0000000 2.6790000 10.7150000 29.4670000
229 229 1 0.0000000 2.6790000 16.0730000 2.6790000
230 230 1 0.0000000 2.6790000 16.0730000 8.0360000
231 231 1 0.0000000 2.6790000 16.0730000 13.3940000
232 232 1 0.0000000 2.6790000 16.0730000 18.7520000
233 233 1 0.0000000 2.6790000 16.0730000 24.1090000
234 234 1 0.0000000 2.6790000 16.0730000 29.4670000
235 235 1 0.0000000 2.6790000 21.4310000 2.6790000
236 236 1 0.0000000 2.6790000 21.4310000 8.0360000
237 237 1 0.0000000 2.6790000 21.4310000 13.3940000
238 238 1 0.0000000 2.6790000 21.4310000 18.7520000
239 239 1 0.0000000 2.6790000 21.4310000 24.1090000
240 240 1 0.0000000 2.6790000 21.4310000 29.4670000
241 241 1 0.0000000 2.6790000 26.7880000 2.6790000
242 242 1 0.0000000 2.6790000 26.7880000 8.0360000
243 243 1 0.0000000 2.6790000 26.7880000 13.3940000
244 244 1 0.0000000 2.6790000 26.7880000 18.7520000
245 245 1 0.0000000 2.6790000 26.7880000 24.1090000
246 246 1 0.0000000 2.6790000 26.7880000 29.4670000
247 247 1 0.0000000 2.6790000 32.1460000 2.6790000
248 248 1 0.0000000 2.6790000 32.1460000 8.0360000
249 249 1 0.0000000 2.6790000 32.1460000 13.3940000
250 250 1 0.0000000 2.6790000 32.1460000 18.7520000
251 251 1 0.0000000 2.6790000 32.1460000 24.1090000
252 252 1 0.0000000 2.6790000 32.1460000 29.4670000
253 253 1 0.0000000 8.0360000 5.3580000 2.6790000
254 254 1 0.0000000 8.0360000 5.3580000 8.0360000
255 255 1 0.0000000 8.0360000 5.3580000 13.3940000
256 256 1 0.0000000 8.0360000 5.3580000 18.7520000
257 257 1 0.0000000 8.0360000 5.3580000 24.1090000
258 258 1 0.0000000 8.0360000 5.3580000 29.4670000
259 259 1 0.0000000 8.0360000 10.7150000 2.6790000
260 260 1 0.0000000 8.0360000 10.7150000 8.0360000
261 261 1 0.0000000 8.0360000 10.7150000 13.3940000
262 262 1 0.0000000 8.0360000 10.7150000 18.7520000
263 263 1 0.0000000 8.0360000 10.7150000 24.1090000
264 264 1 0.0000000 8.0360000 10.7150000 29.4670000
265 265 1 0.0000000 8.0360000 16.0730000 2.6790000
266 266 1 0.0000000 8.0360000 16.0730000 8.0360000
267 267 1 0.0000000 8.0360000 16.0730000 13.3940000
268 268 1 0.0000000 8.0360000 16.0730000 18.7520000
269 269 1 0.0000000 8.0360000 16.0730000 24.1090000
270 270 1 0.0000000 8.0360000 16.0730000 29.4670000
271 271 1 0.0000000 8.0360000 21.4310000 2.6790000
272 272 1 0.0000000 8.0360000 21.4310000 8.0360000
273 273 1 0.0000000 8.0360000 21.4310000 13.3940000
274 274 1 0.0000000 8.0360000 21.4310000 18.7520000
275 275 1 0.0000000 8.0360000 21.4310000 24.1090000
276 276 1 0.0000000 8.0360000 21.4310000 29.4670000
277 277 1 0.0000000 8.0360000 26.7880000 2.6790000
278 278 1 0.0000000 8.0360000 26.7880000 8.0360000
279 279 1 0.0000000 8.0360000 26.7880000 13.3940000
280 280 1 0.0000000 8.0360000 26.7880000 18.7520000
281 281 1 0.0000000 8.0360000 26.7880000 24.1090000
282 282 1 0.0000000 8.0360000 26.7880000 29.4670000
283 283 1 0.0000000 8.0360000 32.1460000 2.6790000
284 284 1 0.0000000 8.0360000 32.1460000 8.0360000
285 285 1 0.0000000 8.0360000 32.1460000 13.3940000
286 286 1 0.0000000 8.0360000 32.1460000 18.7520000
287 287 1 0.0000000 8.0360000 32.1460000 24.1090000
288 288 1 0.0000000 8.0360000 32.1460000 29.4670000
289 289 1 0.0000000 13.3940000 5.3580000 2.6790000
290 290 1 0.0000000 13.3940000 5.3580000 8.0360000
291 291 1 0.0000000 13.3940000 5.3580000 13.3940000
292 292 1 0.0000000 13.3940000 5.3580000 18.7520000
293 293 1 0.0000000 13.3940000 5.3580000 24.1090000
294 294 1 0.0000000 13.3940000 5.3580000 29.4670000
295 295 1 0.0000000 13.3940000 10.7150000 2.6790000
296 296 1 0.0000000 13.3940000 10.7150000 8.0360000
297 297 1 0.0000000 13.3940000 10.7150000 13.3940000
298 298 1 0.0000000 13.3940000 10.7150000 18.7520000
299 299 1 0.0000000 13.3940000 10.7150000 24.1090000
300 300 1 0.0000000 13.3940000 10.7150000 29.4670000
301 301 1 0.0000000 13.3940000 16.0730000 2.6790000
302 302 1 0.0000000 13.3940000 16.0730000 8.0360000
303 303 1 0.0000000 13.3940000 16.0730000 13.3940000
304 304 1 0.0000000 13.3940000 16.0730000 18.7520000
305 305 1 0.0000000 13.3940000 16.0730000 24.1090000
306 306 1 0.0000000 13.3940000 16.0730000 29.4670000
307 307 1 0.0000000 13.3940000 21.4310000 2.6790000
308 308 1 0.0000000 13.3940000 21.4310000 8.0360000
309 309 1 0.0000000 13.3940000 21.4310000 13.3940000
310 310 1 0.0000000 13.3940000 21.4310000 18.7520000
311 311 1 0.0000000 13.3940000 21.4310000 24.1090000
312 312 1 0.0000000 13.3940000 21.4310000 29.4670000
313 313 1 0.0000000 13.3940000 26.7880000 2.6790000
314 314 1 0.0000000 13.3940000 26.7880000 8.0360000
315 315 1 0.0000000 13.3940000 26.7880000 13.3940000
316 316 1 0.0000000 13.3940000 26.7880000 18.7520000
317 317 1 0.0000000 13.3940000 26.7880000 24.1090000
318 318 1 0.0000000 13.3940000 26.7880000 29.4670000
319 319 1 0.0000000 13.3940000 32.1460000 2.6790000
320 320 1 0.0000000 13.3940000 32.1460000 8.0360000
321 321 1 0.0000000 13.3940000 32.1460000 13.3940000
322 322 1 0.0000000 13.3940000 32.1460000 18.7520000
323 323 1 0.0000000 13.3940000 32.1460000 24.1090000
324 324 1 0.0000000 13.3940000 32.1460000 29.4670000
325 325 1 0.0000000 18.7520000 5.3580000 2.6790000
326 326 1 0.0000000 18.7520000 5.3580000 8.0360000
327 327 1 0.0000000 18.7520000 5.3580000 13.3940000
328 328 1 0.0000000 18.7520000 5.3580000 18.7520000
329 329 1 0.0000000 18.7520000 5.3580000 24.1090000
330 330 1 0.0000000 18.7520000 5.3580000 29.4670000
331 331 1 0.0000000 18.7520000 10.7150000 2.6790000
332 332 1 0.0000000 18.7520000 10.7150000 8.0360000
333 333 1 0.0000000 18.7520000 10.7150000 13.3940000
334 334 1 0.0000000 18.7520000 10.7150000 18.7520000
335 335 1 0.0000000 18.7520000 10.7150000 24.1090000
336 336 1 0.0000000 18.7520000 10.7150000 29.4670000
337 337 1 0.0000000 18.7520000 16.0730000 2.6790000
338 338 1 0.0000000 18.7520000 16.0730000 8.0360000
339 339 1 0.0000000 18.7520000 16.0730000 13.3940000
340 340 1 0.0000000 18.7520000 16.0730000 18.7520000
341 341 1 0.0000000 18.7520000 16.0730000 24.1090000
342 342 1 0.0000000 18.7520000 16.0730000 29.4670000
343 343 1 0.0000000 18.7520000 21.4310000 2.6790000
344 344 1 0.0000000 18.7520000 21.4310000 8.0360000
345 345 1 0.0000000 18.7520000 21.4310000 13.3940000
346 346 1 0.0000000 18.7520000 21.4310000 18.7520000
347 347 1 0.0000000 18.7520000 21.4310000 24.1090000
348 348 1 0.0000000 18.7520000 21.4310000 29.4670000
349 349 1 0.0000000 18.7520000 26.7880000 2.6790000
350 350 1 0.0000000 18.7520000 26.7880000 8.0360000
351 351 1 0.0000000 18.7520000 26.7880000 13.3940000
352 352 1 0.0000000 18.7520000 26.7880000 18.7520000
353 353 1 0.0000000 18.7520000 26.7880000 24.1090000
354 354 1 0.0000000 18.7520000 26.7880000 29.4670000
355 355 1 0.0000000 18.7520000 32.1460000 2.6790000
356 356 1 0.0000000 18.7520000 32.1460000 8.0360000
357 357 1 0.0000000 18.7520000 32.1460000 13.3940000
358 358 1 0.0000000 18.7520000 32.1460000 18.7520000
359 359 1 0.0000000 18.7520000 32.1460000 24.1090000
360 360 1 0.0000000 18.7520000 32.1460000 29.4670000
361 361 1 0.0000000 24.1090000 5.3580000 2.6790000
362 362 1 0.0000000 24.1090000 5.3580000 8.0360000
363 363 1 0.0000000 24.1090000 5.3580000 13.3940000
364 364 1 0.0000000 24.1090000 5.3580000 18.7520000
365 365 1 0.0000000 24.1090000 5.3580000 24.1090000
366 366 1 0.0000000 24.1090000 5.3580000 29.4670000
367 367 1 0.0000000 24.1090000 10.7150000 2.6790000
368 368 1 0.0000000 24.1090000 10.7150000 8.0360000
369 369 1 0.0000000 24.1090000 10.7150000 13.3940000
370 370 1 0.0000000 24.1090000 10.7150000 18.7520000
371 371 1 0.0000000 24.1090000 10.7150000 24.1090000
372 372 1 0.0000000 24.1090000 10.7150000 29.4670000
373 373 1 0.0000000 24.1090000 16.0730000 2.6790000
374 374 1 0.0000000 24.1090000 16.0730000 8.0360000
375 375 1 0.0000000 24.1090000 16.0730000 13.3940000
376 376 1 0.0000000 24.1090000 16.0730000 18.7520000
377 377 1 0.0000000 24.1090000 16.0730000 24.1090000
378 378 1 0.0000000 24.1090000 16.0730000 29.4670000
379 379 1 0.0000000 24.1090000 21.4310000 2.6790000
380 380 1 0.0000000 24.1090000 21.4310000 8.0360000
381 381 1 0.0000000 24.1090000 21.4310000 13.3940000
382 382 1 0.0000000 24.1090000 21.4310000 18.7520000
383 383 1 0.0000000 24.1090000 21.4310000 24.1090000
384 384 1 0.0000000 24.1090000 21.4310000 29.4670000
385 385 1 0.0000000 24.1090000 26.7880000 2.6790000
386 386 1 0.0000000 24.1090000 26.7880000 8.0360000
387 387 1 0.0000000 24.1090000 26.7880000 13.3940000
388 388 1 0.0000000 24.1090000 26.7880000 18.7520000
389 389 1 0.0000000 24.1090000 26.7880000 24.1090000
390 390 1 0.0000000 24.1090000 26.7880000 29.4670000
391 391 1 0.0000000 24.1090000 32.1460000 2.6790000
392 392 1 0.0000000 24.1090000 32.1460000 8.0360000
393 393 1 0.0000000 24.1090000 32.1460000 13.3940000
394 394 1 0.0000000 24.1090000 32.1460000 18.7520000
395 395 1 0.0000000 24.1090000 32.1460000 24.1090000
396 396 1 0.0000000 24.1090000 32.1460000 29.4670000
397 397 1 0.0000000 29.4670000 5.3580000 2.6790000
398 398 1 0.0000000 29.4670000 5.3580000 8.0360000
399 399 1 0.0000000 29.4670000 5.3580000 13.3940000
400 400 1 0.0000000 29.4670000 5.3580000 18.7520000
401 401 1 0.0000000 29.4670000 5.3580000 24.1090000
402 402 1 0.0000000 29.4670000 5.3580000 29.4670000
403 403 1 0.0000000 29.4670000 10.7150000 2.6790000
404 404 1 0.0000000 29.4670000 10.7150000 8.0360000
405 405 1 0.0000000 29.4670000 10.7150000 13.3940000
406 406 1 0.0000000 29.4670000 10.7150000 18.7520000
407 407 1 0.0000000 29.4670000 10.7150000 24.1090000
408 408 1 0.0000000 29.4670000 10.7150000 29.4670000
409 409 1 0.0000000 29.4670000 16.0730000 2.6790000
410 410 1 0.0000000 29.4670000 16.0730000 8.0360000
411 411 1 0.0000000 29.4670000 16.0730000 13.3940000
412 412 1 0.0000000 29.4670000 16.0730000 18.7520000
413 413 1 0.0000000 29.4670000 16.0730000 24.1090000
414 414 1 0.0000000 29.4670000 16.0730000 29.4670000
415 415 1 0.0000000 29.4670000 21.4310000 2.6790000
416 416 1 0.0000000 29.4670000 21.4310000 8.0360000
417 417 1 0.0000000 29.4670000 21.4310000 13.3940000
418 418 1 0.0000000 29.4670000 21.4310000 18.7520000
419 419 1 0.0000000 29.4670000 21.4310000 24.1090000
420 420 1 0.0000000 29.4670000 21.4310000 29.4670000
421 421 1 0.0000000 29.4670000 26.7880000 2.6790000
422 422 1 0.0000000 29.4670000 26.7880000 8.0360000
423 423 1 0.0000000 29.4670000 26.7880000 13.3940000
424 424 1 0.0000000 29.4670000 26.7880000 18.7520000
425 425 1 0.0000000 29.4670000 26.7880000 24.1090000
426 426 1 0.0000000 29.4670000 26.7880000 29.4670000
427 427 1 0.0000000 29.4670000 32.1460000 2.6790000
428 428 1 0.0000000 29.4670000 32.1460000 8.0360000
429 429 1 0.0000000 29.4670000 32.1460000 13.3940000
430 430 1 0.0000000 29.4670000 32.1460000 18.7520000
431 431 1 0.0000000 29.4670000 32.1460000 24.1090000
432 432 1 0.0000000 29.4670000 32.1460000 29.4670000
433 433 1 0.0000000 2.6790000 2.6790000 5.3580000
434 434 1 0.0000000 2.6790000 2.6790000 10.7150000
435 435 1 0.0000000 2.6790000 2.6790000 16.0730000
436 436 1 0.0000000 2.6790000 2.6790000 21.4310000
437 437 1 0.0000000 2.6790000 2.6790000 26.7880000
438 438 1 0.0000000 2.6790000 2.6790000 32.1460000
439 439 1 0.0000000 2.6790000 8.0360000 5.3580000
440 440 1 0.0000000 2.6790000 8.0360000 10.7150000
441 441 1 0.0000000 2.6790000 8.0360000 16.0730000
442 442 1 0.0000000 2.6790000 8.0360000 21.4310000
443 443 1 0.0000000 2.6790000 8.0360000 26.7880000
444 444 1 0.0000000 2.6790000 8.0360000 32.1460000
445 445 1 0.0000000 2.6790000 13.3940000 5.3580000
446 446 1 0.0000000 2.6790000 13.3940000 10.7150000
447 447 1 0.0000000 2.6790000 13.3940000 16.0730000
448 448 1 0.0000000 2.6790000 13.3940000 21.4310000
449 449 1 0.0000000 2.6790000 13.3940000 26.7880000
450 450 1 0.0000000 2.6790000 13.3940000 32.1460000
451 451 1 0.0000000 2.6790000 18.7520000 5.3580000
452 452 1 0.0000000 2.6790000 18.7520000 10.7150000
453 453 1 0.0000000 2.6790000 18.7520000 16.0730000
454 454 1 0.0000000 2.6790000 18.7520000 21.4310000
455 455 1 0.0000000 2.6790000 18.7520000 26.7880000
456 456 1 0.0000000 2.6790000 18.7520000 32.1460000
457 457 1 0.0000000 2.6790000 24.1090000 5.3580000
458 458 1 0.0000000 2.6790000 24.1090000 10.7150000
459 459 1 0.0000000 2.6790000 24.1090000 16.0730000
460 460 1 0.0000000 2.6790000 24.1090000 21.4310000
461 461 1 0.0000000 2.6790000 24.1090000 26.7880000
462 462 1 0.0000000 2.6790000 24.1090000 32.1460000
463 463 1 0.0000000 2.6790000 29.4670000 5.3580000
464 464 1 0.0000000 2.6790000 29.4670000 10.7150000
465 465 1 0.0000000 2.6790000 29.4670000 16.0730000
466 466 1 0.0000000 2.6790000 29.4670000 21.4310000
467 467 1 0.0000000 2.6790000 29.4670000 26.7880000
468 468 1 0.0000000 2.6790000 29.4670000 32.1460000
469 469 1 0.0000000 8.0360000 2.6790000 5.3580000
470 470 1 0.0000000 8.0360000 2.6790000 10.7150000
471 471 1 0.0000000 8.0360000 2.6790000 16.0730000
472 472 1 0.0000000 8.0360000 2.6790000 21.4310000
473 473 1 0.0000000 8.0360000 2.6790000 26.7880000
474 474 1 0.0000000 8.0360000 2.6790000 32.1460000
475 475 1 0.0000000 8.0360000 8.0360000 5.3580000
476 476 1 0.0000000 8.0360000 8.0360000 10.7150000
477 477 1 0.0000000 8.0360000 8.0360000 16.0730000
478 478 1 0.0000000 8.0360000 8.0360000 21.4310000
479 479 1 0.0000000 8.0360000 8.0360000 26.7880000
480 480 1 0.0000000 8.0360000 8.0360000 32.1460000
481 481 1 0.0000000 8.0360000 13.3940000 5.3580000
482 482 1 0.0000000 8.0360000 13.3940000 10.7150000
483 483 1 0.0000000 8.0360000 13.3940000 16.0730000
484 484 1 0.0000000 8.0360000 13.3940000 21.4310000
485 485 1 0.0000000 8.0360000 13.3940000 26.7880000
486 486 1 0.0000000 8.0360000 13.3940000 32.1460000
487 487 1 0.0000000 8.0360000 18.7520000 5.3580000
488 488 1 0.0000000 8.0360000 18.7520000 10.7150000
489 489 1 0.0000000 8.0360000 18.7520000 16.0730000
490 490 1 0.0000000 8.0360000 18.7520000 21.4310000
491 491 1 0.0000000 8.0360000 18.7520000 26.7880000
492 492 1 0.0000000 8.0360000 18.7520000 32.1460000
493 493 1 0.0000000 8.0360000 24.1090000 5.3580000
494 494 1 0.0000000 8.0360000 24.1090000 10.7150000
495 495 1 0.0000000 8.0360000 24.1090000 16.0730000
496 496 1 0.0000000 8.0360000 24.1090000 21.4310000
497 497 1 0.0000000 8.0360000 24.1090000 26.7880000
498 498 1 0.0000000 8.0360000 24.1090000 32.1460000
499 499 1 0.0000000 8.0360000 29.4670000 5.3580000
500 500 1 0.0000000 8.0360000 29.4670000 10.7150000
501 501 1 0.0000000 8.0360000 29.4670000 16.0730000
502 502 1 0.0000000 8.0360000 29.4670000 21.4310000
503 503 1 0.0000000 8.0360000 29.4670000 26.7880000
504 504 1 0.0000000 8.0360000 29.4670000 32.1460000
505 505 1 0.0000000 13.3940000 2.6790000 5.3580000
506 506 1 0.0000000 13.3940000 2.6790000 10.7150000
507 507 1 0.0000000 13.3940000 2.6790000 16.0730000
508 508 1 0.0000000 13.3940000 2.6790000 21.4310000
509 509 1 0.0000000 13.3940000 2.6790000 26.7880000
510 510 1 0.0000000 13.3940000 2.6790000 32.1460000
511 511 1 0.0000000 13.3940000 8.0360000 5.3580000
512 512 1 0.0000000 13.3940000 8.0360000 10.7150000
513 513 1 0.0000000 13.3940000 8.0360000 16.0730000
514 514 1 0.0000000 13.3940000 8.0360000 21.4310000
515 515 1 0.0000000 13.3940000 8.0360000 26.7880000
516 516 1 0.0000000 13.3940000 8.0360000 32.1460000
517 517 1 0.0000000 13.3940000 13.3940000 5.3580000
518 518 1 0.0000000 13.3940000 13.3940000 10.7150000
519 519 1 0.0000000 13.3940000 13.3940000 16.0730000
520 520 1 0.0000000 13.3940000 13.3940000 21.4310000
521 521 1 0.0000000 13.3940000 13.3940000 26.7880000
522 522 1 0.0000000 13.3940000 13.3940000 32.1460000
523 523 1 0.0000000 13.3940000 18.7520000 5.3580000
524 524 1 0.0000000 13.3940000 18.7520000 10.7150000
525 525 1 0.0000000 13.3940000 18.7520000 16.0730000
526 526 1 0.0000000 13.3940000 18.7520000 21.4310000
527 527 1 0.0000000 13.3940000 18.7520000 26.7880000
528 528 1 0.0000000 13.3940000 18.7520000 32.1460000
529 529 1 0.0000000 13.3940000 24.1090000 5.3580000
530 530 1 0.0000000 13.3940000 24.1090000 10.7150000
531 531 1 0.0000000 13.3940000 24.1090000 16.0730000
532 532 1 0.0000000 13.3940000 24.1090000 21.4310000
533 533 1 0.0000000 13.3940000 24.1090000 26.7880000
534 534 1 0.0000000 13.3940000 24.1090000 32.1460000
535 535 1 0.0000000 13.3940000 29.4670000 5.3580000
536 536 1 0.0000000 13.3940000 29.4670000 10.7150000
537 537 1 0.0000000 13.3940000 29.4670000 16.0730000
538 538 1 0.0000000 13.3940000 29.4670000 21.4310000
539 539 1 0.0000000 13.3940000 29.4670000 26.7880000
540 540 1 0.0000000 13.3940000 29.4670000 32.1460000
541 541 1 0.0000000 18.7520000 2.6790000 5.3580000
542 542 1 0.0000000 18.7520000 2.6790000 10.7150000
543 543 1 0.0000000 18.7520000 2.6790000 16.0730000
544 544 1 0.0000000 18.7520000 2.6790000 21.4310000
545 545 1 0.0000000 18.7520000 2.6790000 26.7880000
546 546 1 0.0000000 18.7520000 2.6790000 32.1460000
547 547 1 0.0000000 18.7520000 8.0360000 5.3580000
548 548 1 0.0000000 18.7520000 8.0360000 10.7150000
549 549 1 0.0000000 18.7520000 8.0360000 16.0730000
550 550 1 0.0000000 18.7520000 8.0360000 21.4310000
551 551 1 0.0000000 18.7520000 8.0360000 26.7880000
552 552 1 0.0000000 18.7520000 8.0360000 32.1460000
553 553 1 0.0000000 18.7520000 13.3940000 5.3580000
554 554 1 0.0000000 18.7520000 13.3940000 10.7150000
555 555 1 0.0000000 18.7520000 13.3940000 16.0730000
556 556 1 0.0000000 18.7520000 13.3940000 21.4310000
557 557 1 0.0000000 18.7520000 13.3940000 26.7880000
558 558 1 0.0000000 18.7520000 13.3940000 32.1460000
559 559 1 0.0000000 18.7520000 18.7520000 5.3580000
560 560 1 0.0000000 18.7520000 18.7520000 10.7150000
561 561 1 0.0000000 18.7520000 18.7520000 16.0730000
562 562 1 0.0000000 18.7520000 18.7520000 21.4310000
563 563 1 0.0000000 18.7520000 18.7520000 26.7880000
564 564 1 0.0000000 18.7520000 18.7520000 32.1460000
565 565 1 0.0000000 18.7520000 24.1090000 5.3580000
566 566 1 0.0000000 18.7520000 24.1090000 10.7150000
567 567 1 0.0000000 18.7520000 24.1090000 16.0730000
568 568 1 0.0000000 18.7520000 24.1090000 21.4310000
569 569 1 0.0000000 18.7520000 24.1090000 26.7880000
570 570 1 0.0000000 18.7520000 24.1090000 32.1460000
571 571 1 0.0000000 18.7520000 29.4670000 5.3580000
572 572 1 0.0000000 18.7520000 29.4670000 10.7150000
573 573 1 0.0000000 18.7520000 29.4670000 16.0730000
574 574 1 0.0000000 18.7520000 29.4670000 21.4310000
575 575 1 0.0000000 18.7520000 29.4670000 26.7880000
576 576 1 0.0000000 18.7520000 29.4670000 32.1460000
577 577 1 0.0000000 24.1090000 2.6790000 5.3580000
578 578 1 0.0000000 24.1090000 2.6790000 10.7150000
579 579 1 0.0000000 24.1090000 2.6790000 16.0730000
580 580 1 0.0000000 24.1090000 2.6790000 21.4310000
581 581 1 0.0000000 24.1090000 2.6790000 26.7880000
582 582 1 0.0000000 24.1090000 2.6790000 32.1460000
583 583 1 0.0000000 24.1090000 8.0360000 5.3580000
584 584 1 0.0000000 24.1090000 8.0360000 10.7150000
585 585 1 0.0000000 24.1090000 8.0360000 16.0730000
586 586 1 0.0000000 24.1090000 8.0360000 21.4310000
587 587 1 0.0000000 24.1090000 8.0360000 26.7880000
588 588 1 0.0000000 24.1090000 8.0360000 32.1460000
589 589 1 0.0000000 24.1090000 13.3940000 5.3580000
590 590 1 0.0000000 24.1090000 13.3940000 10.7150000
591 591 1 0.0000000 24.1090000 13.3940000 16.0730000
592 592 1 0.0000000 24.1090000 13.3940000 21.4310000
593 593 1 0.0000000 24.1090000 13.3940000 26.7880000
594 594 1 0.0000000 24.1090000 13.3940000 32.1460000
595 595 1 0.0000000 24.1090000 18.7520000 5.3580000
596 596 1 0.0000000 24.1090000 18.7520000 10.7150000
597 597 1 0.0000000 24.1090000 18.7520000 16.0730000
598 598 1 0.0000000 24.1090000 18.7520000 21.4310000
599 599 1 0.0000000 24.1090000 18.7520000 26.7880000
600 600 1 0.0000000 24.1090000 18.7520000 32.1460000
601 601 1 0.0000000 24.1090000 24.1090000 5.3580000
602 602 1 0.0000000 24.1090000 24.1090000 10.7150000
603 603 1 0.0000000 24.1090000 24.1090000 16.0730000
604 604 1 0.0000000 24.1090000 24.1090000 21.4310000
605 605 1 0.0000000 24.1090000 24.1090000 26.7880000
606 606 1 0.0000000 24.1090000 24.1090000 32.1460000
607 607 1 0.0000000 24.1090000 29.4670000 5.3580000
608 608 1 0.0000000 24.1090000 29.4670000 10.7150000
609 609 1 0.0000000 24.1090000 29.4670000 16.0730000
610 610 1 0.0000000 24.1090000 29.4670000 21.4310000
611 611 1 0.0000000 24.1090000 29.4670000 26.7880000
612 612 1 0.0000000 24.1090000 29.4670000 32.1460000
613 613 1 0.0000000 29.4670000 2.6790000 5.3580000
614 614 1 0.0000000 29.4670000 2.6790000 10.7150000
615 615 1 0.0000000 29.4670000 2.6790000 16.0730000
616 616 1 0.0000000 29.4670000 2.6790000 21.4310000
617 617 1 0.0000000 29.4670000 2.6790000 26.7880000
618 618 1 0.0000000 29.4670000 2.6790000 32.1460000
619 619 1 0.0000000 29.4670000 8.0360000 5.3580000
620 620 1 0.0000000 29.4670000 8.0360000 10.7150000
621 621 1 0.0000000 29.4670000 8.0360000 16.0730000
622 622 1 0.0000000 29.4670000 8.0360000 21.4310000
623 623 1 0.0000000 29.4670000 8.0360000 26.7880000
624 624 1 0.0000000 29.4670000 8.0360000 32.1460000
625 625 1 0.0000000 29.4670000 13.3940000 5.3580000
626 626 1 0.0000000 29.4670000 13.3940000 10.7150000
627 627 1 0.0000000 29.4670000 13.3940000 16.0730000
628 628 1 0.0000000 29.4670000 13.3940000 21.4310000
629 629 1 0.0000000 29.4670000 13.3940000 26.7880000
630 630 1 0.0000000 29.4670000 13.3940000 32.1460000
631 631 1 0.0000000 29.4670000 18.7520000 5.3580000
632 632 1 0.0000000 29.4670000 18.7520000 10.7150000
633 633 1 0.0000000 29.4670000 18.7520000 16.0730000
634 634 1 0.0000000 29.4670000 18.7520000 21.4310000
635 635 1 0.0000000 29.4670000 18.7520000 26.7880000
636 636 1 0.0000000 29.4670000 18.7520000 32.1460000
637 637 1 0.0000000 29.4670000 24.1090000 5.3580000
638 638 1 0.0000000 29.4670000 24.1090000 10.7150000
639 639 1 0.0000000 29.4670000 24.1090000 16.0730000
640 640 1 0.0000000 29.4670000 24.1090000 21.4310000
641 641 1 0.0000000 29.4670000 24.1090000 26.7880000
642 642 1 0.0000000 29.4670000 24.1090000 32.1460000
643 643 1 0.0000000 29.4670000 29.4670000 5.3580000
644 644 1 0.0000000 29.4670000 29.4670000 10.7150000
645 645 1 0.0000000 29.4670000 29.4670000 16.0730000
646 646 1 0.0000000 29.4670000 29.4670000 21.4310000
647 647 1 0.0000000 29.4670000 29.4670000 26.7880000
648 648 1 0.0000000 29.4670000 29.4670000 32.1460000
649 649 1 0.0000000 0.0000000 5.3580000 5.3580000
650 650 1 0.0000000 0.0000000 5.3580000 10.7150000
651 651 1 0.0000000 0.0000000 5.3580000 16.0730000
652 652 1 0.0000000 0.0000000 5.3580000 21.4310000
653 653 1 0.0000000 0.0000000 5.3580000 26.7880000
654 654 1 0.0000000 0.0000000 5.3580000 32.1460000
655 655 1 0.0000000 0.0000000 10.7150000 5.3580000
656 656 1 0.0000000 0.0000000 10.7150000 10.7150000
657 657 1 0.0000000 0.0000000 10.7150000 16.0730000
658 658 1 0.0000000 0.0000000 10.7150000 21.4310000
659 659 1 0.0000000 0.0000000 10.7150000 26.7880000
660 660 1 0.0000000 0.0000000 10.7150000 32.1460000
661 661 1 0.0000000 0.0000000 16.0730000 5.3580000
662 662 1 0.0000000 0.0000000 16.0730000 10.7150000
663 663 1 0.0000000 0.0000000 16.0730000 16.0730000
664 664 1 0.0000000 0.0000000 16.0730000 21.4310000
665 665 1 0.0000000 0.0000000 16.0730000 26.7880000
666 666 1 0.0000000 0.0000000 16.0730000 32.1460000
667 667 1 0.0000000 0.0000000 21.4310000 5.3580000
668 668 1 0.0000000 0.0000000 21.4310000 10.7150000
669 669 1 0.0000000 0.0000000 21.4310000 16.0730000
670 670 1 0.0000000 0.0000000 21.4310000 21.4310000
671 671 1 0.0000000 0.0000000 21.4310000 26.7880000
672 672 1 0.0000000 0.0000000 21.4310000 32.1460000
673 673 1 0.0000000 0.0000000 26.7880000 5.3580000
674 674 1 0.0000000 0.0000000 26.7880000 10.7150000
675 675 1 0.0000000 0.0000000 26.7880000 16.0730000
676 676 1 0.0000000 0.0000000 26.7880000 21.4310000
677 677 1 0.0000000 0.0000000 26.7880000 26.7880000
678 678 1 0.0000000 0.0000000 26.7880000 32.1460000
679 679 1 0.0000000 0.0000000 32.1460000 5.3580000
680 680 1 0.0000000 0.0000000 32.1460000 10.7150000
681 681 1 0.0000000 0.0000000 32.1460000 16.0730000
682 682 1 0.0000000 0.0000000 32.1460000 21.4310000
683 683 1 0.0000000 0.0000000 32.1460000 26.7880000
684 684 1 0.0000000 0.0000000 32.1460000 32.1460000
685 685 1 0.0000000 5.3580000 5.3580000 5.3580000
686 686 1 0.0000000 5.3580000 5.3580000 10.7150000
687 687 1 0.0000000 5.3580000 5.3580000 16.0730000
688 688 1 0.0000000 5.3580000 5.3580000 21.4310000
689 689 1 0.0000000 5.3580000 5.3580000 26.7880000
690 690 1 0.0000000 5.3580000 5.3580000 32.1460000
691 691 1 0.0000000 5.3580000 10.7150000 5.3580000
692 692 1 0.0000000 5.3580000 10.7150000 10.7150000
693 693 1 0.0000000 5.3580000 10.7150000 16.0730000
694 694 1 0.0000000 5.3580000 10.7150000 21.4310000
695 695 1 0.0000000 5.3580000 10.7150000 26.7880000
696 696 1 0.0000000 5.3580000 10.7150000 32.1460000
697 697 1 0.0000000 5.3580000 16.0730000 5.3580000
698 698 1 0.0000000 5.3580000 16.0730000 10.7150000
699 699 1 0.0000000 5.3580000 16.0730000 16.0730000
700 700 1 0.0000000 5.3580000 16.0730000 21.4310000
701 701 1 0.0000000 5.3580000 16.0730000 26.7880000
702 702 1 0.0000000 5.3580000 16.0730000 32.1460000
703 703 1 0.0000000 5.3580000 21.4310000 5.3580000
704 704 1 0.0000000 5.3580000 21.4310000 10.7150000
705 705 1 0.0000000 5.3580000 21.4310000 16.0730000
706 706 1 0.0000000 5.3580000 21.4310000 21.4310000
707 707 1 0.0000000 5.3580000 21.4310000 26.7880000
708 708 1 0.0000000 5.3580000 21.4310000 32.1460000
709 709 1 0.0000000 5.3580000 26.7880000 5.3580000
710 710 1 0.0000000 5.3580000 26.7880000 10.7150000
711 711 1 0.0000000 5.3580000 26.7880000 16.0730000
712 712 1 0.0000000 5.3580000 26.7880000 21.4310000
713 713 1 0.0000000 5.3580000 26.7880000 26.7880000
714 714 1 0.0000000 5.3580000 26.7880000 32.1460000
715 715 1 0.0000000 5.3580000 32.1460000 5.3580000
716 716 1 0.0000000 5.3580000 32.1460000 10.7150000
717 717 1 0.0000000 5.3580000 32.1460000 16.0730000
718 718 1 0.0000000 5.3580000 32.1460000 21.4310000
719 719 1 0.0000000 5.3580000 32.1460000 26.7880000
720 720 1 0.0000000 5.3580000 32.1460000 32.1460000
721 721 1 0.0000000 10.7150000 5.3580000 5.3580000
722 722 1 0.0000000 10.7150000 5.3580000 10.7150000
723 723 1 0.0000000 10.7150000 5.3580000 16.0730000
724 724 1 0.0000000 10.7150000 5.3580000 21.4310000
725 725 1 0.0000000 10.7150000 5.3580000 26.7880000
726 726 1 0.0000000 10.7150000 5.3580000 32.1460000
727 727 1 0.0000000 10.7150000 10.7150000 5.3580000
728 728 1 0.0000000 10.7150000 10.7150000 10.7150000
729 729 1 0.0000000 10.7150000 10.7150000 16.0730000
730 730 1 0.0000000 10.7150000 10.7150000 21.4310000
731 731 1 0.0000000 10.7150000 10.7150000 26.7880000
732 732 1 0.0000000 10.7150000 10.7150000 32.1460000
733 733 1 0.0000000 10.7150000 16.0730000 5.3580000
734 734 1 0.0000000 10.7150000 16.0730000 10.7150000
735 735 1 0.0000000 10.7150000 16.0730000 16.0730000
736 736 1 0.0000000 10.7150000 16.0730000 21.4310000
737 737 1 0.0000000 10.7150000 16.0730000 26.7880000
738 738 1 0.0000000 10.7150000 16.0730000 32.1460000
739 739 1 0.0000000 10.7150000 21.4310000 5.3580000
740 740 1 0.0000000 10.7150000 21.4310000 10.7150000
741 741 1 0.0000000 10.7150000 21.4310000 16.0730000
742 742 1 0.0000000 10.7150000 21.4310000 21.4310000
743 743 1 0.0000000 10.7150000 21.4310000 26.7880000
744 744 1 0.0000000 10.7150000 21.4310000 32.1460000
745 745 1 0.0000000 10.7150000 26.7880000 5.3580000
746 746 1 0.0000000 10.7150000 26.7880000 10.7150000
747 747 1 0.0000000 10.7150000 26.7880000 16.0730000
748 748 1 0.0000000 10.7150000 26.7880000 21.4310000
749 749 1 0.0000000 10.7150000 26.7880000 26.7880000
750 750 1 0.0000000 10.7150000 26.7880000 32.1460000
751 751 1 0.0000000 10.7150000 32.1460000 5.3580000
752 752 1 0.0000000 10.7150000 32.1460000 10.7150000
753 753 1 0.0000000 10.7150000 32.1460000 16.0730000
754 754 1 0.0000000 10.7150000 32.1460000 21.4310000
755 755 1 0.0000000 10.7150000 32.1460000 26.7880000
756 756 1 0.0000000 10.7150000 32.1460000 32.1460000
757 757 1 0.0000000 16.0730000 5.3580000 5.3580000
758 758 1 0.0000000 16.0730000 5.3580000 10.7150000
759 759 1 0.0000000 16.0730000 5.3580000 16.0730000
760 760 1 0.0000000 16.0730000 5.3580000 21.4310000
761 761 1 0.0000000 16.0730000 5.3580000 26.7880000
762 762 1 0.0000000 16.0730000 5.3580000 32.1460000
763 763 1 0.0000000 16.0730000 10.7150000 5.3580000
764 764 1 0.0000000 16.0730000 10.7150000 10.7150000
765 765 1 0.0000000 16.0730000 10.7150000 16.0730000
766 766 1 0.0000000 16.0730000 10.7150000 21.4310000
767 767 1 0.0000000 16.0730000 10.7150000 26.7880000
768 768 1 0.0000000 16.0730000 10.7150000 32.1460000
769 769 1 0.0000000 16.0730000 16.0730000 5.3580000
770 770 1 0.0000000 16.0730000 16.0730000 10.7150000
771 771 1 0.0000000 16.0730000 16.0730000 16.0730000
772 772 1 0.0000000 16.0730000 16.0730000 21.4310000
773 773 1 0.0000000 16.0730000 16.0730000 26.7880000
774 774 1 0.0000000 16.0730000 16.0730000 32.1460000
775 775 1 0.0000000 16.0730000 21.4310000 5.3580000
776 776 1 0.0000000 16.0730000 21.4310000 10.7150000
777 777 1 0.0000000 16.0730000 21.4310000 16.0730000
778 778 1 0.0000000 16.0730000 21.4310000 21.4310000
779 779 1 0.0000000 16.0730000 21.4310000 26.7880000
780 780 1 0.0000000 16.0730000 21.4310000 32.1460000
781 781 1 0.0000000 16.0730000 26.7880000 5.3580000
782 782 1 0.0000000 16.0730000 26.7880000 10.7150000
783 783 1 0.0000000 16.0730000 26.7880000 16.0730000
784 784 1 0.0000000 16.0730000 26.7880000 21.4310000
785 785 1 0.0000000 16.0730000 26.7880000 26.7880000
786 786 1 0.0000000 16.0730000 26.7880000 32.1460000
787 787 1 0.0000000 16.0730000 32.1460000 5.3580000
788 788 1 0.0000000 16.0730000 32.1460000 10.7150000
789 789 1 0.0000000 16.0730000 32.1460000 16.0730000
790 790 1 0.0000000 16.0730000 32.1460000 21.4310000
791 791 1 0.0000000 16.0730000 32.1460000 26.7880000
792 792 1 0.0000000 16.0730000 32.1460000 32.1460000
793 793 1 0.0000000 21.4310000 5.3580000 5.3580000
794 794 1 0.0000000 21.4310000 5.3580000 10.7150000
795 795 1 0.0000000 21.4310000 5.3580000 16.0730000
796 796 1 0.0000000 21.4310000 5.3580000 21.4310000
797 797 1 0.0000000 21.4310000 5.3580000 26.7880000
798 798 1 0.0000000 21.4310000 5.3580000 32.1460000
799 799 1 0.0000000 21.4310000 10.7150000 5.3580000
800 800 1 0.0000000 21.4310000 10.7150000 10.7150000
801 801 1 0.0000000 21.4310000 10.7150000 16.0730000
802 802 1 0.0000000 21.4310000 10.7150000 21.4310000
803 803 1 0.0000000 21.4310000 10.7150000 26.7880000
804 804 1 0.0000000 21.4310000 10.7150000 32.1460000
805 805 1 0.0000000 21.4310000 16.0730000 5.3580000
806 806 1 0.0000000 21.4310000 16.0730000 10.7150000
807 807 1 0.0000000 21.4310000 16.0730000 16.0730000
808 808 1 0.0000000 21.4310000 16.0730000 21.4310000
809 809 1 0.0000000 21.4310000 16.0730000 26.7880000
810 810 1 0.0000000 21.4310000 16.0730000 32.1460000
811 811 1 0.0000000 21.4310000 21.4310000 5.3580000
812 812 1 0.0000000 21.4310000 21.4310000 10.7150000
813 813 1 0.0000000 21.4310000 21.4310000 16.0730000
814 814 1 0.0000000 21.4310000 21.4310000 21.4310000
815 815 1 0.0000000 21.4310000 21.4310000 26.7880000
816 816 1 0.0000000 21.4310000 21.4310000 32.1460000
817 817 1 0.0000000 21.4310000 26.7880000 5.3580000
818 818 1 0.0000000 21.4310000 26.7880000 10.7150000
819 819 1 0.0000000 21.4310000 26.7880000 16.0730000
820 820 1 0.0000000 21.4310000 26.7880000 21.4310000
821 821 1 0.0000000 21.4310000 26.7880000 26.7880000
822 822 1 0.0000000 21.4310000 26.7880000 32.1460000
823 823 1 0.0000000 21.4310000 32.1460000 5.3580000
824 824 1 0.0000000 21.4310000 32.1460000 10.7150000
825 825 1 0.0000000 21.4310000 32.1460000 16.0730000
826 826 1 0.0000000 21.4310000 32.1460000 21.4310000
827 827 1 0.0000000 21.4310000 32.1460000 26.7880000
828 828 1 0.0000000 21.4310000 32.1460000 32.1460000
829 829 1 0.0000000 26.7880000 5.3580000 5.3580000
830 830 1 0.0000000 26.7880000 5.3580000 10.7150000
831 831 1 0.0000000 26.7880000 5.3580000 16.0730000
832 832 1 0.0000000 26.7880000 5.3580000 21.4310000
833 833 1 0.0000000 26.7880000 5.3580000 26.7880000
834 834 1 0.0000000 26.7880000 5.3580000 32.1460000
835 835 1 0.0000000 26.7880000 10.7150000 5.3580000
836 836 1 0.0000000 26.7880000 10.7150000 10.7150000
837 837 1 0.0000000 26.7880000 10.7150000 16.0730000
838 838 1 0.0000000 26.7880000 10.7150000 21.4310000
839 839 1 0.0000000 26.7880000 10.7150000 26.7880000
840 840 1 0.0000000 26.7880000 10.7150000 32.1460000
841 841 1 0.0000000 26.7880000 16.0730000 5.3580000
842 842 1 0.0000000 26.7880000 16.0730000 10.7150000
843 843 1 0.0000000 26.7880000 16.0730000 16.0730000
844 844 1 0.0000000 26.7880000 16.0730000 21.4310000
845 845 1 0.0000000 26.7880000 16.0730000 26.7880000
846 846 1 0.0000000 26.7880000 16.0730000 32.1460000
847 847 1 0.0000000 26.7880000 21.4310000 5.3580000
848 848 1 0.0000000 26.7880000 21.4310000 10.7150000
849 849 1 0.0000000 26.7880000 21.4310000 16.0730000
850 850 1 0.0000000 26.7880000 21.4310000 21.4310000
851 851 1 0.0000000 26.7880000 21.4310000 26.7880000
852 852 1 0.0000000 26.7880000 21.4310000 32.1460000
853 853 1 0.0000000 26.7880000 26.7880000 5.3580000
854 854 1 0.0000000 26.7880000 26.7880000 10.7150000
855 855 1 0.0000000 26.7880000 26.7880000 16.0730000
856 856 1 0.0000000 26.7880000 26.7880000 21.4310000
857 857 1 0.0000000 26.7880000 26.7880000 26.7880000
858 858 1 0.0000000 26.7880000 26.7880000 32.1460000
859 859 1 0.0000000 26.7880000 32.1460000 5.3580000
860 860 1 0.0000000 26.7880000 32.1460000 10.7150000
861 861 1 0.0000000 26.7880000 32.1460000 16.0730000
862 862 1 0.0000000 26.7880000 32.1460000 21.4310000
863 863 1 0.0000000 26.7880000 32.1460000 26.7880000
864 864 1 0.0000000 26.7880000 32.1460000 32.1460000

20
examples/gjf/ff-argon.lmp Normal file
View File

@ -0,0 +1,20 @@
#############################
#Atoms types - mass - charge#
#############################
#@ 1 atom types #!THIS LINE IS NECESSARY DON'T SPEND HOURS FINDING THAT OUT!#
variable Ar equal 1
#############
#Atom Masses#
#############
mass ${Ar} 39.903
###########################
#Pair Potentials - Tersoff#
###########################
pair_style lj/cubic
pair_coeff * * 0.0102701 3.42

23
examples/gjf/in.gjf.vfull Normal file
View File

@ -0,0 +1,23 @@
# GJF-2GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
include ff-argon.lmp
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
fix nve all nve
fix lang all langevin 10 10 1 26488 gjf vfull
thermo 200
run 50000

23
examples/gjf/in.gjf.vhalf Normal file
View File

@ -0,0 +1,23 @@
# GJF-2GJ thermostat
units metal
atom_style full
boundary p p p
read_data argon.lmp
include ff-argon.lmp
velocity all create 10 2357 mom yes dist gaussian
neighbor 1 bin
timestep 0.1
fix nve all nve
fix lang all langevin 10 10 1 26488 gjf vhalf
thermo 200
run 50000

4
src/.gitignore vendored
View File

@ -161,6 +161,10 @@
/fix_setforce_spin.h
/min_spin.cpp
/min_spin.h
/min_spin_cg.cpp
/min_spin_cg.h
/min_spin_lbfgs.cpp
/min_spin_lbfgs.h
/neb_spin.cpp
/neb_spin.h
/pair_spin.cpp

166
src/KOKKOS/fix_langevin_kokkos.cpp
View File

@ -61,7 +61,6 @@ FixLangevinKokkos<DeviceType>::FixLangevinKokkos(LAMMPS *lmp, int narg, char **a
k_ratio.template modify<LMPHostType>();
if(gjfflag){
nvalues = 3;
grow_arrays(atomKK->nmax);
atom->add_callback(0);
// initialize franprev to zero
@ -69,8 +68,12 @@ FixLangevinKokkos<DeviceType>::FixLangevinKokkos(LAMMPS *lmp, int narg, char **a
franprev[i][0] = 0.0;
franprev[i][1] = 0.0;
franprev[i][2] = 0.0;
lv[i][0] = 0.0;
lv[i][1] = 0.0;
lv[i][2] = 0.0;
}
k_franprev.template modify<LMPHostType>();
k_lv.template modify<LMPHostType>();
}
if(zeroflag){
k_fsumall = tdual_double_1d_3n("langevin:fsumall");
@ -94,6 +97,7 @@ FixLangevinKokkos<DeviceType>::~FixLangevinKokkos()
memoryKK->destroy_kokkos(k_ratio,ratio);
memoryKK->destroy_kokkos(k_flangevin,flangevin);
if(gjfflag) memoryKK->destroy_kokkos(k_franprev,franprev);
if(gjfflag) memoryKK->destroy_kokkos(k_lv,lv);
memoryKK->destroy_kokkos(k_tforce,tforce);
}
@ -107,6 +111,10 @@ void FixLangevinKokkos<DeviceType>::init()
error->all(FLERR,"Fix langevin omega is not yet implemented with kokkos");
if(ascale)
error->all(FLERR,"Fix langevin angmom is not yet implemented with kokkos");
if(gjfflag && tbiasflag)
error->all(FLERR,"Fix langevin gjf + tbias is not yet implemented with kokkos");
if(gjfflag && tbiasflag)
error->warning(FLERR,"Fix langevin gjf + kokkos is not implemented with random gaussians");
// prefactors are modified in the init
k_gfactor1.template modify<LMPHostType>();
@ -121,6 +129,40 @@ void FixLangevinKokkos<DeviceType>::grow_arrays(int nmax)
memoryKK->grow_kokkos(k_franprev,franprev,nmax,3,"langevin:franprev");
d_franprev = k_franprev.template view<DeviceType>();
h_franprev = k_franprev.template view<LMPHostType>();
memoryKK->grow_kokkos(k_lv,lv,nmax,3,"langevin:lv");
d_lv = k_lv.template view<DeviceType>();
h_lv = k_lv.template view<LMPHostType>();
}
/* ---------------------------------------------------------------------- */
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::initial_integrate(int vflag)
{
atomKK->sync(execution_space,datamask_read);
atomKK->modified(execution_space,datamask_modify);
v = atomKK->k_v.view<DeviceType>();
f = atomKK->k_f.view<DeviceType>();
int nlocal = atomKK->nlocal;
if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst;
FixLangevinKokkosInitialIntegrateFunctor<DeviceType> functor(this);
Kokkos::parallel_for(nlocal,functor);
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixLangevinKokkos<DeviceType>::initial_integrate_item(int i) const
{
if (mask[i] & groupbit) {
f(i,0) /= gjfa;
f(i,1) /= gjfa;
f(i,2) /= gjfa;
v(i,0) = d_lv(i,0);
v(i,1) = d_lv(i,1);
v(i,2) = d_lv(i,2);
}
}
/* ---------------------------------------------------------------------- */
@ -140,6 +182,7 @@ void FixLangevinKokkos<DeviceType>::post_force(int vflag)
k_gfactor2.template sync<DeviceType>();
k_ratio.template sync<DeviceType>();
if(gjfflag) k_franprev.template sync<DeviceType>();
if(gjfflag) k_lv.template sync<DeviceType>();
boltz = force->boltz;
dt = update->dt;
@ -177,7 +220,7 @@ void FixLangevinKokkos<DeviceType>::post_force(int vflag)
atomKK->sync(temperature->execution_space,temperature->datamask_read);
temperature->compute_scalar();
temperature->remove_bias_all(); // modifies velocities
// if temeprature compute is kokkosized host-devcie comm won't be needed
// if temeprature compute is kokkosized host-device comm won't be needed
atomKK->modified(temperature->execution_space,temperature->datamask_modify);
atomKK->sync(execution_space,temperature->datamask_modify);
}
@ -481,6 +524,7 @@ void FixLangevinKokkos<DeviceType>::post_force(int vflag)
// set modify flags for the views modified in post_force functor
if (gjfflag) k_franprev.template modify<DeviceType>();
if (gjfflag) k_lv.template modify<DeviceType>();
if (tallyflag) k_flangevin.template modify<DeviceType>();
// set total force to zero
@ -550,6 +594,10 @@ FSUM FixLangevinKokkos<DeviceType>::post_force_item(int i) const
}
if (Tp_GJF) {
d_lv(i,0) = gjfsib*v(i,0);
d_lv(i,1) = gjfsib*v(i,1);
d_lv(i,2) = gjfsib*v(i,2);
fswap = 0.5*(fran[0]+d_franprev(i,0));
d_franprev(i,0) = fran[0];
fran[0] = fswap;
@ -560,15 +608,15 @@ FSUM FixLangevinKokkos<DeviceType>::post_force_item(int i) const
d_franprev(i,2) = fran[2];
fran[2] = fswap;
fdrag[0] *= gjffac;
fdrag[1] *= gjffac;
fdrag[2] *= gjffac;
fran[0] *= gjffac;
fran[1] *= gjffac;
fran[2] *= gjffac;
f(i,0) *= gjffac;
f(i,1) *= gjffac;
f(i,2) *= gjffac;
fdrag[0] *= gjfa;
fdrag[1] *= gjfa;
fdrag[2] *= gjfa;
fran[0] *= gjfa;
fran[1] *= gjfa;
fran[2] *= gjfa;
f(i,0) *= gjfa;
f(i,1) *= gjfa;
f(i,2) *= gjfa;
}
f(i,0) += fdrag[0] + fran[0];
@ -576,6 +624,17 @@ FSUM FixLangevinKokkos<DeviceType>::post_force_item(int i) const
f(i,2) += fdrag[2] + fran[2];
if (Tp_TALLY) {
if (Tp_GJF){
fdrag[0] = gamma1*d_lv(i,0)/gjfsib/gjfsib;
fdrag[1] = gamma1*d_lv(i,1)/gjfsib/gjfsib;
fdrag[2] = gamma1*d_lv(i,2)/gjfsib/gjfsib;
fswap = (2*fran[0]/gjfa - d_franprev(i,0))/gjfsib;
fran[0] = fswap;
fswap = (2*fran[1]/gjfa - d_franprev(i,1))/gjfsib;
fran[1] = fswap;
fswap = (2*fran[2]/gjfa - d_franprev(i,2))/gjfsib;
fran[2] = fswap;
}
d_flangevin(i,0) = fdrag[0] + fran[0];
d_flangevin(i,1) = fdrag[1] + fran[1];
d_flangevin(i,2) = fdrag[2] + fran[2];
@ -719,9 +778,10 @@ double FixLangevinKokkos<DeviceType>::compute_energy_item(int i) const
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::end_of_step()
{
if (!tallyflag) return;
if (!tallyflag && !gjfflag) return;
v = atomKK->k_v.template view<DeviceType>();
f = atomKK->k_f.template view<DeviceType>();
mask = atomKK->k_mask.template view<DeviceType>();
atomKK->sync(execution_space,V_MASK | MASK_MASK);
@ -733,9 +793,81 @@ void FixLangevinKokkos<DeviceType>::end_of_step()
FixLangevinKokkosTallyEnergyFunctor<DeviceType> tally_functor(this);
Kokkos::parallel_reduce(nlocal,tally_functor,energy_onestep);
if (gjfflag){
if (rmass.data()) {
FixLangevinKokkosEndOfStepFunctor<DeviceType,1> functor(this);
Kokkos::parallel_for(nlocal,functor);
} else {
mass = atomKK->k_mass.view<DeviceType>();
FixLangevinKokkosEndOfStepFunctor<DeviceType,0> functor(this);
Kokkos::parallel_for(nlocal,functor);
}
}
energy += energy_onestep*update->dt;
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixLangevinKokkos<DeviceType>::end_of_step_item(int i) const {
double tmp[3];
if (mask[i] & groupbit) {
const double dtfm = force->ftm2v * 0.5 * dt / mass[type[i]];
tmp[0] = v(i,0);
tmp[1] = v(i,1);
tmp[2] = v(i,2);
if (!osflag){
v(i,0) = d_lv(i,0);
v(i,1) = d_lv(i,1);
v(i,2) = d_lv(i,2);
} else {
v(i,0) = 0.5 * gjfsib * gjfsib * (v(i,0) + dtfm * f(i,0) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,0) - d_franprev(i,0)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,0);
v(i,1) = 0.5 * gjfsib * gjfsib * (v(i,1) + dtfm * f(i,1) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,0) - d_franprev(i,1)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,1);
v(i,2) = 0.5 * gjfsib * gjfsib * (v(i,2) + dtfm * f(i,2) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,0) - d_franprev(i,2)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,2);
}
d_lv(i,0) = tmp[0];
d_lv(i,1) = tmp[1];
d_lv(i,2) = tmp[2];
}
}
template<class DeviceType>
KOKKOS_INLINE_FUNCTION
void FixLangevinKokkos<DeviceType>::end_of_step_rmass_item(int i) const
{
double tmp[3];
if (mask[i] & groupbit) {
const double dtfm = force->ftm2v * 0.5 * dt / rmass[i];
tmp[0] = v(i,0);
tmp[1] = v(i,1);
tmp[2] = v(i,2);
if (!osflag){
v(i,0) = d_lv(i,0);
v(i,1) = d_lv(i,1);
v(i,2) = d_lv(i,2);
} else {
v(i,0) = 0.5 * gjfsib * gjfsib * (v(i,0) + dtfm * f(i,0) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,0) - d_franprev(i,0)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,0);
v(i,1) = 0.5 * gjfsib * gjfsib * (v(i,1) + dtfm * f(i,1) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,1) - d_franprev(i,1)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,1);
v(i,2) = 0.5 * gjfsib * gjfsib * (v(i,2) + dtfm * f(i,2) / gjfa) +
dtfm * 0.5 * (gjfsib * d_flangevin(i,2) - d_franprev(i,2)) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * d_lv(i,2);
}
d_lv(i,0) = tmp[0];
d_lv(i,1) = tmp[1];
d_lv(i,2) = tmp[2];
}
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
------------------------------------------------------------------------- */
@ -743,10 +875,15 @@ void FixLangevinKokkos<DeviceType>::end_of_step()
template<class DeviceType>
void FixLangevinKokkos<DeviceType>::copy_arrays(int i, int j, int delflag)
{
for (int m = 0; m < nvalues; m++)
h_franprev(j,m) = h_franprev(i,m);
h_franprev(j,0) = h_franprev(i,0);
h_franprev(j,1) = h_franprev(i,1);
h_franprev(j,2) = h_franprev(i,2);
h_lv(j,0) = h_lv(i,0);
h_lv(j,1) = h_lv(i,1);
h_lv(j,2) = h_lv(i,2);
k_franprev.template modify<LMPHostType>();
k_lv.template modify<LMPHostType>();
}
@ -765,6 +902,7 @@ void FixLangevinKokkos<DeviceType>::cleanup_copy()
tforce = NULL;
gjfflag = 0;
franprev = NULL;
lv = NULL;
id = style = NULL;
vatom = NULL;
}

59
src/KOKKOS/fix_langevin_kokkos.h
View File

@ -56,6 +56,9 @@ namespace LAMMPS_NS {
template<class DeviceType>
class FixLangevinKokkos;
template <class DeviceType>
class FixLangevinKokkosInitialIntegrateFunctor;
template<class DeviceType,int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO>
class FixLangevinKokkosPostForceFunctor;
@ -72,6 +75,7 @@ namespace LAMMPS_NS {
void cleanup_copy();
void init();
void initial_integrate(int);
void post_force(int);
void reset_dt();
void grow_arrays(int);
@ -79,6 +83,12 @@ namespace LAMMPS_NS {
double compute_scalar();
void end_of_step();
KOKKOS_INLINE_FUNCTION
void initial_integrate_item(int) const;
KOKKOS_INLINE_FUNCTION
void initial_integrate_rmass_item(int) const;
template<int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO>
KOKKOS_INLINE_FUNCTION
@ -90,14 +100,25 @@ namespace LAMMPS_NS {
KOKKOS_INLINE_FUNCTION
double compute_energy_item(int) const;
KOKKOS_INLINE_FUNCTION
void end_of_step_item(int) const;
KOKKOS_INLINE_FUNCTION
void end_of_step_rmass_item(int) const;
private:
class CommKokkos *commKK;
typename ArrayTypes<DeviceType>::t_float_1d rmass;
typename ArrayTypes<DeviceType>::t_float_1d mass;
typename ArrayTypes<DeviceType>::tdual_double_2d k_franprev;
typename ArrayTypes<DeviceType>::t_double_2d d_franprev;
HAT::t_double_2d h_franprev;
typename ArrayTypes<DeviceType>::tdual_double_2d k_lv;
typename ArrayTypes<DeviceType>::t_double_2d d_lv;
HAT::t_double_2d h_lv;
typename ArrayTypes<DeviceType>::tdual_double_2d k_flangevin;
typename ArrayTypes<DeviceType>::t_double_2d d_flangevin;
HAT::t_double_2d h_flangevin;
@ -130,6 +151,21 @@ namespace LAMMPS_NS {
};
template <class DeviceType>
struct FixLangevinKokkosInitialIntegrateFunctor {
typedef DeviceType device_type ;
FixLangevinKokkos<DeviceType> c;
FixLangevinKokkosInitialIntegrateFunctor(FixLangevinKokkos<DeviceType>* c_ptr):
c(*c_ptr) {c.cleanup_copy();};
KOKKOS_INLINE_FUNCTION
void operator()(const int i) const {
c.initial_integrate_item(i);
}
};
template <class DeviceType,int Tp_TSTYLEATOM, int Tp_GJF, int Tp_TALLY,
int Tp_BIAS, int Tp_RMASS, int Tp_ZERO>
struct FixLangevinKokkosPostForceFunctor {
@ -207,6 +243,21 @@ namespace LAMMPS_NS {
update += source;
}
};
template <class DeviceType, int RMass>
struct FixLangevinKokkosEndOfStepFunctor {
typedef DeviceType device_type ;
FixLangevinKokkos<DeviceType> c;
FixLangevinKokkosEndOfStepFunctor(FixLangevinKokkos<DeviceType>* c_ptr):
c(*c_ptr) {c.cleanup_copy();}
KOKKOS_INLINE_FUNCTION
void operator()(const int i) const {
if (RMass) c.end_of_step_rmass_item(i);
else c.end_of_step_item(i);
}
};
}
#endif
@ -231,4 +282,12 @@ E: Fix langevin variable returned negative temperature
Self-explanatory.
E: Fix langevin gjf with tbias is not yet implemented with kokkos
This option is not yet available.
W: Fix langevin gjf using random gaussians is not implemented with kokkos
This will most likely cause errors in kinetic fluctuations.
*/

5
src/SPIN/fix_nve_spin.cpp
View File

@ -91,12 +91,17 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) :
// defining lattice_flag
// changing the lattice option, from (yes,no) -> (moving,frozen)
// for now, (yes,no) still works (to avoid user's confusions).
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg],"lattice") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix/NVE/spin command");
if (strcmp(arg[iarg+1],"no") == 0) lattice_flag = 0;
else if (strcmp(arg[iarg+1],"frozen") == 0) lattice_flag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) lattice_flag = 1;
else if (strcmp(arg[iarg+1],"moving") == 0) lattice_flag = 1;
else error->all(FLERR,"Illegal fix/NVE/spin command");
iarg += 2;
} else error->all(FLERR,"Illegal fix/NVE/spin command");

2
src/SPIN/fix_nve_spin.h
View File

@ -52,7 +52,7 @@ friend class PairSpin;
double dtv, dtf, dts; // velocity, force, and spin timesteps
int nlocal_max; // max value of nlocal (for lists size)
int nlocal_max; // max value of nlocal (for size of lists)
int pair_spin_flag; // magnetic pair flags
int long_spin_flag; // magnetic long-range flag

59
src/SPIN/min_spin.cpp
View File

@ -41,15 +41,15 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
MinSpin::MinSpin(LAMMPS *lmp) : Min(lmp) {}
MinSpin::MinSpin(LAMMPS *lmp) : Min(lmp) {
alpha_damp = 1.0;
discrete_factor = 10.0;
}
/* ---------------------------------------------------------------------- */
void MinSpin::init()
{
alpha_damp = 1.0;
discrete_factor = 10.0;
Min::init();
dts = dt = update->dt;
@ -77,12 +77,12 @@ void MinSpin::setup_style()
int MinSpin::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"alpha_damp") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
if (narg < 2) error->all(FLERR,"Illegal min_modify command");
alpha_damp = force->numeric(FLERR,arg[1]);
return 2;
}
if (strcmp(arg[0],"discrete_factor") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
if (narg < 2) error->all(FLERR,"Illegal min_modify command");
discrete_factor = force->numeric(FLERR,arg[1]);
return 2;
}
@ -116,7 +116,7 @@ void MinSpin::reset_vectors()
int MinSpin::iterate(int maxiter)
{
bigint ntimestep;
double fmdotfm;
double fmdotfm,fmsq;
int flag,flagall;
for (int iter = 0; iter < maxiter; iter++) {
@ -130,7 +130,7 @@ int MinSpin::iterate(int maxiter)
// optimize timestep accross processes / replicas
// need a force calculation for timestep optimization
energy_force(0);
if (iter == 0) energy_force(0);
dts = evaluate_dt();
// apply damped precessional dynamics to the spins
@ -163,8 +163,13 @@ int MinSpin::iterate(int maxiter)
// magnetic torque tolerance criterion
// sync across replicas if running multi-replica minimization
fmdotfm = fmsq = 0.0;
if (update->ftol > 0.0) {
fmdotfm = fmnorm_sqr();
if (normstyle == MAX) fmsq = max_torque(); // max torque norm
else if (normstyle == INF) fmsq = inf_torque(); // inf torque norm
else if (normstyle == TWO) fmsq = total_torque(); // Euclidean torque 2-norm
else error->all(FLERR,"Illegal min_modify command");
fmdotfm = fmsq*fmsq;
if (update->multireplica == 0) {
if (fmdotfm < update->ftol*update->ftol) return FTOL;
} else {
@ -242,7 +247,7 @@ void MinSpin::advance_spins(double dts)
double **sp = atom->sp;
double **fm = atom->fm;
double tdampx,tdampy,tdampz;
double msq,scale,fm2,energy,dts2;
double fm2,energy,dts2;
double cp[3],g[3];
dts2 = dts*dts;
@ -288,37 +293,3 @@ void MinSpin::advance_spins(double dts)
// because no need for simplecticity
}
}
/* ----------------------------------------------------------------------
compute and return ||mag. torque||_2^2
------------------------------------------------------------------------- */
double MinSpin::fmnorm_sqr()
{
int nlocal = atom->nlocal;
double tx,ty,tz;
double **sp = atom->sp;
double **fm = atom->fm;
// calc. magnetic torques
double local_norm2_sqr = 0.0;
for (int i = 0; i < nlocal; i++) {
tx = (fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]);
ty = (fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]);
tz = (fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]);
local_norm2_sqr += tx*tx + ty*ty + tz*tz;
}
// no extra atom calc. for spins
if (nextra_atom)
error->all(FLERR,"extra atom option not available yet");
double norm2_sqr = 0.0;
MPI_Allreduce(&local_norm2_sqr,&norm2_sqr,1,MPI_DOUBLE,MPI_SUM,world);
return norm2_sqr;
}

1
src/SPIN/min_spin.h
View File

@ -35,7 +35,6 @@ class MinSpin : public Min {
int iterate(int);
double evaluate_dt();
void advance_spins(double);
double fmnorm_sqr();
private:

649
src/SPIN/min_spin_cg.cpp Normal file
View File

@ -0,0 +1,649 @@
/* ----------------------------------------------------------------------
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: Aleksei Ivanov (University of Iceland)
Julien Tranchida (SNL)
Please cite the related publication:
Ivanov, A. V., Uzdin, V. M., & Jónsson, H. (2019). Fast and Robust
Algorithm for the Minimisation of the Energy of Spin Systems. arXiv
preprint arXiv:1904.02669.
------------------------------------------------------------------------- */
#include <mpi.h>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "min_spin_cg.h"
#include "universe.h"
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "force.h"
#include "update.h"
#include "output.h"
#include "timer.h"
#include "error.h"
#include "memory.h"
#include "modify.h"
#include "math_special.h"
#include "math_const.h"
#include "universe.h"
using namespace LAMMPS_NS;
using namespace MathConst;
static const char cite_minstyle_spin_cg[] =
"min_style spin/cg command:\n\n"
"@article{ivanov2019fast,\n"
"title={Fast and Robust Algorithm for the Minimisation of the Energy of "
"Spin Systems},\n"
"author={Ivanov, A. V and Uzdin, V. M. and J{\'o}nsson, H.},\n"
"journal={arXiv preprint arXiv:1904.02669},\n"
"year={2019}\n"
"}\n\n";
// EPS_ENERGY = minimum normalization for energy tolerance
#define EPS_ENERGY 1.0e-8
#define DELAYSTEP 5
/* ---------------------------------------------------------------------- */
MinSpinCG::MinSpinCG(LAMMPS *lmp) :
Min(lmp), g_old(NULL), g_cur(NULL), p_s(NULL), sp_copy(NULL)
{
if (lmp->citeme) lmp->citeme->add(cite_minstyle_spin_cg);
nlocal_max = 0;
// nreplica = number of partitions
// ireplica = which world I am in universe
nreplica = universe->nworlds;
ireplica = universe->iworld;
use_line_search = 0; // no line search as default option for CG
discrete_factor = 10.0;
}
/* ---------------------------------------------------------------------- */
MinSpinCG::~MinSpinCG()
{
memory->destroy(g_old);
memory->destroy(g_cur);
memory->destroy(p_s);
if (use_line_search)
memory->destroy(sp_copy);
}
/* ---------------------------------------------------------------------- */
void MinSpinCG::init()
{
local_iter = 0;
der_e_cur = 0.0;
der_e_pr = 0.0;
Min::init();
// warning if line_search combined to gneb
if ((nreplica >= 1) && (linestyle != 4) && (comm->me == 0))
error->warning(FLERR,"Line search incompatible gneb");
// set back use_line_search to 0 if more than one replica
if (linestyle == 3 && nreplica == 1){
use_line_search = 1;
}
else{
use_line_search = 0;
}
dts = dt = update->dt;
last_negative = update->ntimestep;
// allocate tables
nlocal_max = atom->nlocal;
memory->grow(g_old,3*nlocal_max,"min/spin/cg:g_old");
memory->grow(g_cur,3*nlocal_max,"min/spin/cg:g_cur");
memory->grow(p_s,3*nlocal_max,"min/spin/cg:p_s");
if (use_line_search)
memory->grow(sp_copy,nlocal_max,3,"min/spin/cg:sp_copy");
}
/* ---------------------------------------------------------------------- */
void MinSpinCG::setup_style()
{
double **v = atom->v;
int nlocal = atom->nlocal;
// check if the atom/spin style is defined
if (!atom->sp_flag)
error->all(FLERR,"min spin/cg requires atom/spin style");
for (int i = 0; i < nlocal; i++)
v[i][0] = v[i][1] = v[i][2] = 0.0;
}
/* ---------------------------------------------------------------------- */
int MinSpinCG::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"discrete_factor") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
discrete_factor = force->numeric(FLERR,arg[1]);
return 2;
}
return 0;
}
/* ----------------------------------------------------------------------
set current vector lengths and pointers
called after atoms have migrated
------------------------------------------------------------------------- */
void MinSpinCG::reset_vectors()
{
// atomic dof
// size sp is 4N vector
nvec = 4 * atom->nlocal;
if (nvec) spvec = atom->sp[0];
nvec = 3 * atom->nlocal;
if (nvec) fmvec = atom->fm[0];
if (nvec) xvec = atom->x[0];
if (nvec) fvec = atom->f[0];
}
/* ----------------------------------------------------------------------
minimization via orthogonal spin optimisation
------------------------------------------------------------------------- */
int MinSpinCG::iterate(int maxiter)
{
int nlocal = atom->nlocal;
bigint ntimestep;
double fmdotfm,fmsq;
int flag, flagall;
double **sp = atom->sp;
double der_e_cur_tmp = 0.0;
if (nlocal_max < nlocal) {
local_iter = 0;
nlocal_max = nlocal;
memory->grow(g_old,3*nlocal_max,"min/spin/cg:g_old");
memory->grow(g_cur,3*nlocal_max,"min/spin/cg:g_cur");
memory->grow(p_s,3*nlocal_max,"min/spin/cg:p_s");
if (use_line_search)
memory->grow(sp_copy,nlocal_max,3,"min/spin/cg:sp_copy");
}
for (int iter = 0; iter < maxiter; iter++) {
if (timer->check_timeout(niter))
return TIMEOUT;
ntimestep = ++update->ntimestep;
niter++;
// optimize timestep accross processes / replicas
// need a force calculation for timestep optimization
if (use_line_search) {
// here we need to do line search
if (local_iter == 0){
calc_gradient();
}
calc_search_direction();
der_e_cur = 0.0;
for (int i = 0; i < 3 * nlocal; i++)
der_e_cur += g_cur[i] * p_s[i];
MPI_Allreduce(&der_e_cur,&der_e_cur_tmp,1,MPI_DOUBLE,MPI_SUM,world);
der_e_cur = der_e_cur_tmp;
if (update->multireplica == 1) {
MPI_Allreduce(&der_e_cur_tmp,&der_e_cur,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
}
for (int i = 0; i < nlocal; i++)
for (int j = 0; j < 3; j++)
sp_copy[i][j] = sp[i][j];
eprevious = ecurrent;
der_e_pr = der_e_cur;
calc_and_make_step(0.0, 1.0, 0);
}
else{
// here we don't do line search
// if gneb calc., nreplica > 1
// then calculate gradients and advance spins
// of intermediate replicas only
calc_gradient();
calc_search_direction();
advance_spins();
neval++;
eprevious = ecurrent;
ecurrent = energy_force(0);
}
// energy tolerance criterion
// only check after DELAYSTEP elapsed since velocties reset to 0
// sync across replicas if running multi-replica minimization
if (update->etol > 0.0 && ntimestep-last_negative > DELAYSTEP) {
if (update->multireplica == 0) {
if (fabs(ecurrent-eprevious) <
update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY))
return ETOL;
} else {
if (fabs(ecurrent-eprevious) <
update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY))
flag = 0;
else flag = 1;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
if (flagall == 0) return ETOL;
}
}
// magnetic torque tolerance criterion
// sync across replicas if running multi-replica minimization
fmdotfm = fmsq = 0.0;
if (update->ftol > 0.0) {
if (normstyle == MAX) fmsq = max_torque(); // max torque norm
else if (normstyle == INF) fmsq = inf_torque(); // inf torque norm
else if (normstyle == TWO) fmsq = total_torque(); // Euclidean torque 2-norm
else error->all(FLERR,"Illegal min_modify command");
fmdotfm = fmsq*fmsq;
if (update->multireplica == 0) {
if (fmdotfm < update->ftol*update->ftol) return FTOL;
} else {
if (fmdotfm < update->ftol*update->ftol) flag = 0;
else flag = 1;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
if (flagall == 0) return FTOL;
}
}
// output for thermo, dump, restart files
if (output->next == ntimestep) {
timer->stamp();
output->write(ntimestep);
timer->stamp(Timer::OUTPUT);
}
}
return MAXITER;
}
/* ----------------------------------------------------------------------
calculate gradients
---------------------------------------------------------------------- */
void MinSpinCG::calc_gradient()
{
int nlocal = atom->nlocal;
double **sp = atom->sp;
double **fm = atom->fm;
double hbar = force->hplanck/MY_2PI;
double factor;
if (use_line_search)
factor = hbar;
else factor = evaluate_dt();
// loop on all spins on proc.
for (int i = 0; i < nlocal; i++) {
g_cur[3 * i + 0] = (fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]) * factor;
g_cur[3 * i + 1] = -(fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]) * factor;
g_cur[3 * i + 2] = (fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]) * factor;
}
}
/* ----------------------------------------------------------------------
search direction:
The Fletcher-Reeves conj. grad. method
See Jorge Nocedal and Stephen J. Wright 'Numerical
Optimization' Second Edition, 2006 (p. 121)
---------------------------------------------------------------------- */
void MinSpinCG::calc_search_direction()
{
int nlocal = atom->nlocal;
double g2old = 0.0;
double g2 = 0.0;
double beta = 0.0;
double g2_global = 0.0;
double g2old_global = 0.0;
double factor = 1.0;
// for multiple replica do not move end points
if (nreplica > 1)
if (ireplica == 0 || ireplica == nreplica - 1)
factor = 0.0;
if (local_iter == 0 || local_iter % 5 == 0){ // steepest descent direction
for (int i = 0; i < 3 * nlocal; i++) {
p_s[i] = -g_cur[i] * factor;
g_old[i] = g_cur[i] * factor;
}
} else { // conjugate direction
for (int i = 0; i < 3 * nlocal; i++) {
g2old += g_old[i] * g_old[i];
g2 += g_cur[i] * g_cur[i];
}
// now we need to collect/broadcast beta on this replica
// need to check what is beta for GNEB
MPI_Allreduce(&g2,&g2_global,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&g2old,&g2old_global,1,MPI_DOUBLE,MPI_SUM,world);
// Sum over all replicas. Good for GNEB.
if (nreplica > 1) {
g2 = g2_global * factor;
g2old = g2old_global * factor;
MPI_Allreduce(&g2,&g2_global,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
MPI_Allreduce(&g2old,&g2old_global,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
}
if (fabs(g2_global) < 1.0e-60) beta = 0.0;
else beta = g2_global / g2old_global;
// calculate conjugate direction
for (int i = 0; i < 3 * nlocal; i++) {
p_s[i] = (beta * p_s[i] - g_cur[i]) * factor;
g_old[i] = g_cur[i] * factor;
}
}
local_iter++;
}
/* ----------------------------------------------------------------------
rotation of spins along the search direction
---------------------------------------------------------------------- */
void MinSpinCG::advance_spins()
{
int nlocal = atom->nlocal;
double **sp = atom->sp;
double rot_mat[9]; // exponential of matrix made of search direction
double s_new[3];
// loop on all spins on proc.
for (int i = 0; i < nlocal; i++) {
rodrigues_rotation(p_s + 3 * i, rot_mat);
// rotate spins
vm3(rot_mat, sp[i], s_new);
for (int j = 0; j < 3; j++) sp[i][j] = s_new[j];
}
}
/* ----------------------------------------------------------------------
calculate 3x3 matrix exponential using Rodrigues' formula
(R. Murray, Z. Li, and S. Shankar Sastry,
A Mathematical Introduction to
Robotic Manipulation (1994), p. 28 and 30).
upp_tr - vector x, y, z so that one calculate
U = exp(A) with A= [[0, x, y],
[-x, 0, z],
[-y, -z, 0]]
------------------------------------------------------------------------- */
void MinSpinCG::rodrigues_rotation(const double *upp_tr, double *out)
{
double theta,A,B,D,x,y,z;
double s1,s2,s3,a1,a2,a3;
if (fabs(upp_tr[0]) < 1.0e-40 &&
fabs(upp_tr[1]) < 1.0e-40 &&
fabs(upp_tr[2]) < 1.0e-40){
// if upp_tr is zero, return unity matrix
for(int k = 0; k < 3; k++){
for(int m = 0; m < 3; m++){
if (m == k) out[3 * k + m] = 1.0;
else out[3 * k + m] = 0.0;
}
}
return;
}
theta = sqrt(upp_tr[0] * upp_tr[0] +
upp_tr[1] * upp_tr[1] +
upp_tr[2] * upp_tr[2]);
A = cos(theta);
B = sin(theta);
D = 1.0 - A;
x = upp_tr[0]/theta;
y = upp_tr[1]/theta;
z = upp_tr[2]/theta;
// diagonal elements of U
out[0] = A + z * z * D;
out[4] = A + y * y * D;
out[8] = A + x * x * D;
// off diagonal of U
s1 = -y * z *D;
s2 = x * z * D;
s3 = -x * y * D;
a1 = x * B;
a2 = y * B;
a3 = z * B;
out[1] = s1 + a1;
out[3] = s1 - a1;
out[2] = s2 + a2;
out[6] = s2 - a2;
out[5] = s3 + a3;
out[7] = s3 - a3;
}
/* ----------------------------------------------------------------------
out = vector^T x m,
m -- 3x3 matrix , v -- 3-d vector
------------------------------------------------------------------------- */
void MinSpinCG::vm3(const double *m, const double *v, double *out)
{
for(int i = 0; i < 3; i++){
out[i] = 0.0;
for(int j = 0; j < 3; j++) out[i] += *(m + 3 * j + i) * v[j];
}
}
/* ----------------------------------------------------------------------
advance spins
------------------------------------------------------------------------- */
void MinSpinCG::make_step(double c, double *energy_and_der)
{
double p_scaled[3];
int nlocal = atom->nlocal;
double rot_mat[9]; // exponential of matrix made of search direction
double s_new[3];
double **sp = atom->sp;
double der_e_cur_tmp = 0.0;
for (int i = 0; i < nlocal; i++) {
// scale the search direction
for (int j = 0; j < 3; j++) p_scaled[j] = c * p_s[3 * i + j];
// calculate rotation matrix
rodrigues_rotation(p_scaled, rot_mat);
// rotate spins
vm3(rot_mat, sp[i], s_new);
for (int j = 0; j < 3; j++) sp[i][j] = s_new[j];
}
ecurrent = energy_force(0);
calc_gradient();
neval++;
der_e_cur = 0.0;
for (int i = 0; i < 3 * nlocal; i++) {
der_e_cur += g_cur[i] * p_s[i];
}
MPI_Allreduce(&der_e_cur,&der_e_cur_tmp,1,MPI_DOUBLE,MPI_SUM,world);
der_e_cur = der_e_cur_tmp;
if (update->multireplica == 1) {
MPI_Allreduce(&der_e_cur_tmp,&der_e_cur,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
}
energy_and_der[0] = ecurrent;
energy_and_der[1] = der_e_cur;
}
/* ----------------------------------------------------------------------
Calculate step length which satisfies approximate Wolfe conditions
using the cubic interpolation
------------------------------------------------------------------------- */
int MinSpinCG::calc_and_make_step(double a, double b, int index)
{
double e_and_d[2] = {0.0,0.0};
double alpha,c1,c2,c3;
double **sp = atom->sp;
int nlocal = atom->nlocal;
make_step(b,e_and_d);
ecurrent = e_and_d[0];
der_e_cur = e_and_d[1];
index++;
if (adescent(eprevious,e_and_d[0]) || index == 5){
MPI_Bcast(&b,1,MPI_DOUBLE,0,world);
for (int i = 0; i < 3 * nlocal; i++) {
p_s[i] = b * p_s[i];
}
return 1;
}
else {
double r,f0,f1,df0,df1;
r = b - a;
f0 = eprevious;
f1 = ecurrent;
df0 = der_e_pr;
df1 = der_e_cur;
c1 = -2.0*(f1-f0)/(r*r*r)+(df1+df0)/(r*r);
c2 = 3.0*(f1-f0)/(r*r)-(df1+2.0*df0)/(r);
c3 = df0;
// f(x) = c1 x^3 + c2 x^2 + c3 x^1 + c4
// has minimum at alpha below. We do not check boundaries.
alpha = (-c2 + sqrt(c2*c2 - 3.0*c1*c3))/(3.0*c1);
MPI_Bcast(&alpha,1,MPI_DOUBLE,0,world);
if (alpha < 0.0) alpha = r/2.0;
for (int i = 0; i < nlocal; i++) {
for (int j = 0; j < 3; j++) sp[i][j] = sp_copy[i][j];
}
calc_and_make_step(0.0, alpha, index);
}
return 0;
}
/* ----------------------------------------------------------------------
Approximate descent
------------------------------------------------------------------------- */
int MinSpinCG::adescent(double phi_0, double phi_j){
double eps = 1.0e-6;
if (phi_j<=phi_0+eps*fabs(phi_0))
return 1;
else
return 0;
}
/* ----------------------------------------------------------------------
evaluate max timestep
---------------------------------------------------------------------- */
double MinSpinCG::evaluate_dt()
{
double dtmax;
double fmsq;
double fmaxsqone,fmaxsqloc,fmaxsqall;
int nlocal = atom->nlocal;
double **fm = atom->fm;
// finding max fm on this proc.
fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0;
for (int i = 0; i < nlocal; i++) {
fmsq = fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2];
fmaxsqone = MAX(fmaxsqone,fmsq);
}
// finding max fm on this replica
fmaxsqloc = fmaxsqone;
MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world);
// finding max fm over all replicas, if necessary
// this communicator would be invalid for multiprocess replicas
fmaxsqall = fmaxsqloc;
if (update->multireplica == 1) {
fmaxsqall = fmaxsqloc;
MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld);
}
if (fmaxsqall == 0.0)
error->all(FLERR,"Incorrect fmaxsqall calculation");
// define max timestep by dividing by the
// inverse of max frequency by discrete_factor
dtmax = MY_2PI/(discrete_factor*sqrt(fmaxsqall));
return dtmax;
}

71
src/SPIN/min_spin_cg.h Normal file
View File

@ -0,0 +1,71 @@
/* -*- 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 MINIMIZE_CLASS
MinimizeStyle(spin/cg, MinSpinCG)
#else
#ifndef LMP_MIN_SPIN_CG_H
#define LMP_MIN_SPIN_CG_H
#include "min.h"
namespace LAMMPS_NS {
class MinSpinCG: public Min {
public:
MinSpinCG(class LAMMPS *);
virtual ~MinSpinCG();
void init();
void setup_style();
void reset_vectors();
int modify_param(int, char **);
int iterate(int);
private:
int local_iter; // for neb
int nlocal_max; // max value of nlocal (for size of lists)
int use_line_search; // use line search or not.
int ireplica,nreplica; // for neb
double dt; // global timestep
double dts; // spin timestep
double discrete_factor; // factor for spin timestep evaluation
double der_e_cur; // current derivative along search dir.
double der_e_pr; // previous derivative along search dir.
double *spvec; // variables for atomic dof, as 1d vector
double *fmvec; // variables for atomic dof, as 1d vector
double *g_old; // gradient vector at previous step
double *g_cur; // current gradient vector
double *p_s; // search direction vector
double **sp_copy; // copy of the spins
void advance_spins();
void calc_gradient();
void calc_search_direction();
void vm3(const double *, const double *, double *);
void rodrigues_rotation(const double *, double *);
void make_step(double, double *);
int calc_and_make_step(double, double, int);
int adescent(double, double);
double evaluate_dt();
double maximum_rotation(double *);
bigint last_negative;
};
}
#endif
#endif

754
src/SPIN/min_spin_lbfgs.cpp Normal file
View File

@ -0,0 +1,754 @@
/* ----------------------------------------------------------------------
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: Aleksei Ivanov (University of Iceland)
Julien Tranchida (SNL)
Please cite the related publication:
Ivanov, A. V., Uzdin, V. M., & Jónsson, H. (2019). Fast and Robust
Algorithm for the Minimisation of the Energy of Spin Systems. arXiv
preprint arXiv:1904.02669.
------------------------------------------------------------------------- */
#include <mpi.h>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "min_spin_lbfgs.h"
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "force.h"
#include "update.h"
#include "output.h"
#include "timer.h"
#include "error.h"
#include "memory.h"
#include "modify.h"
#include "math_special.h"
#include "math_const.h"
#include "universe.h"
using namespace LAMMPS_NS;
using namespace MathConst;
static const char cite_minstyle_spin_lbfgs[] =
"min_style spin/lbfgs command:\n\n"
"@article{ivanov2019fast,\n"
"title={Fast and Robust Algorithm for the Minimisation of the Energy of "
"Spin Systems},\n"
"author={Ivanov, A. V and Uzdin, V. M. and J{\'o}nsson, H.},\n"
"journal={arXiv preprint arXiv:1904.02669},\n"
"year={2019}\n"
"}\n\n";
// EPS_ENERGY = minimum normalization for energy tolerance
#define EPS_ENERGY 1.0e-8
#define DELAYSTEP 5
/* ---------------------------------------------------------------------- */
MinSpinLBFGS::MinSpinLBFGS(LAMMPS *lmp) :
Min(lmp), g_old(NULL), g_cur(NULL), p_s(NULL), rho(NULL), ds(NULL), dy(NULL), sp_copy(NULL)
{
if (lmp->citeme) lmp->citeme->add(cite_minstyle_spin_lbfgs);
nlocal_max = 0;
// nreplica = number of partitions
// ireplica = which world I am in universe
nreplica = universe->nworlds;
ireplica = universe->iworld;
use_line_search = 0; // no line search as default option for LBFGS
maxepsrot = MY_2PI / (100.0);
}
/* ---------------------------------------------------------------------- */
MinSpinLBFGS::~MinSpinLBFGS()
{
memory->destroy(g_old);
memory->destroy(g_cur);
memory->destroy(p_s);
memory->destroy(ds);
memory->destroy(dy);
memory->destroy(rho);
if (use_line_search)
memory->destroy(sp_copy);
}
/* ---------------------------------------------------------------------- */
void MinSpinLBFGS::init()
{
num_mem = 3;
local_iter = 0;
der_e_cur = 0.0;
der_e_pr = 0.0;
Min::init();
// warning if line_search combined to gneb
if ((nreplica >= 1) && (linestyle != 4) && (comm->me == 0))
error->warning(FLERR,"Line search incompatible gneb");
// set back use_line_search to 0 if more than one replica
if (linestyle == 3 && nreplica == 1){
use_line_search = 1;
}
else{
use_line_search = 0;
}
last_negative = update->ntimestep;
// allocate tables
nlocal_max = atom->nlocal;
memory->grow(g_old,3*nlocal_max,"min/spin/lbfgs:g_old");
memory->grow(g_cur,3*nlocal_max,"min/spin/lbfgs:g_cur");
memory->grow(p_s,3*nlocal_max,"min/spin/lbfgs:p_s");
memory->grow(rho,num_mem,"min/spin/lbfgs:rho");
memory->grow(ds,num_mem,3*nlocal_max,"min/spin/lbfgs:ds");
memory->grow(dy,num_mem,3*nlocal_max,"min/spin/lbfgs:dy");
if (use_line_search)
memory->grow(sp_copy,nlocal_max,3,"min/spin/lbfgs:sp_copy");
}
/* ---------------------------------------------------------------------- */
void MinSpinLBFGS::setup_style()
{
double **v = atom->v;
int nlocal = atom->nlocal;
// check if the atom/spin style is defined
if (!atom->sp_flag)
error->all(FLERR,"min spin/lbfgs requires atom/spin style");
for (int i = 0; i < nlocal; i++)
v[i][0] = v[i][1] = v[i][2] = 0.0;
}
/* ---------------------------------------------------------------------- */
int MinSpinLBFGS::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"discrete_factor") == 0) {
if (narg < 2) error->all(FLERR,"Illegal min_modify command");
double discrete_factor;
discrete_factor = force->numeric(FLERR,arg[1]);
maxepsrot = MY_2PI / (10 * discrete_factor);
return 2;
}
return 0;
}
/* ----------------------------------------------------------------------
set current vector lengths and pointers
called after atoms have migrated
------------------------------------------------------------------------- */
void MinSpinLBFGS::reset_vectors()
{
// atomic dof
// size sp is 4N vector
nvec = 4 * atom->nlocal;
if (nvec) spvec = atom->sp[0];
nvec = 3 * atom->nlocal;
if (nvec) fmvec = atom->fm[0];
if (nvec) xvec = atom->x[0];
if (nvec) fvec = atom->f[0];
}
/* ----------------------------------------------------------------------
minimization via damped spin dynamics
------------------------------------------------------------------------- */
int MinSpinLBFGS::iterate(int maxiter)
{
int nlocal = atom->nlocal;
bigint ntimestep;
double fmdotfm,fmsq;
int flag, flagall;
double **sp = atom->sp;
double der_e_cur_tmp = 0.0;
if (nlocal_max < nlocal) {
nlocal_max = nlocal;
local_iter = 0;
memory->grow(g_old,3*nlocal_max,"min/spin/lbfgs:g_old");
memory->grow(g_cur,3*nlocal_max,"min/spin/lbfgs:g_cur");
memory->grow(p_s,3*nlocal_max,"min/spin/lbfgs:p_s");
memory->grow(rho,num_mem,"min/spin/lbfgs:rho");
memory->grow(ds,num_mem,3*nlocal_max,"min/spin/lbfgs:ds");
memory->grow(dy,num_mem,3*nlocal_max,"min/spin/lbfgs:dy");
if (use_line_search)
memory->grow(sp_copy,nlocal_max,3,"min/spin/lbfgs:sp_copy");
}
for (int iter = 0; iter < maxiter; iter++) {
if (timer->check_timeout(niter))
return TIMEOUT;
ntimestep = ++update->ntimestep;
niter++;
// optimize timestep accross processes / replicas
// need a force calculation for timestep optimization
if (use_line_search) {
// here we need to do line search
if (local_iter == 0){
eprevious = ecurrent;
ecurrent = energy_force(0);
calc_gradient();
}
calc_search_direction();
der_e_cur = 0.0;
for (int i = 0; i < 3 * nlocal; i++)
der_e_cur += g_cur[i] * p_s[i];
MPI_Allreduce(&der_e_cur,&der_e_cur_tmp,1,MPI_DOUBLE,MPI_SUM,world);
der_e_cur = der_e_cur_tmp;
if (update->multireplica == 1) {
MPI_Allreduce(&der_e_cur_tmp,&der_e_cur,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
}
for (int i = 0; i < nlocal; i++)
for (int j = 0; j < 3; j++)
sp_copy[i][j] = sp[i][j];
eprevious = ecurrent;
der_e_pr = der_e_cur;
calc_and_make_step(0.0, 1.0, 0);
}
else{
// here we don't do line search
// but use cutoff rotation angle
// if gneb calc., nreplica > 1
// then calculate gradients and advance spins
// of intermediate replicas only
eprevious = ecurrent;
ecurrent = energy_force(0);
calc_gradient();
calc_search_direction();
advance_spins();
neval++;
}
// energy tolerance criterion
// only check after DELAYSTEP elapsed since velocties reset to 0
// sync across replicas if running multi-replica minimization
if (update->etol > 0.0 && ntimestep-last_negative > DELAYSTEP) {
if (update->multireplica == 0) {
if (fabs(ecurrent-eprevious) <
update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY))
return ETOL;
} else {
if (fabs(ecurrent-eprevious) <
update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY))
flag = 0;
else flag = 1;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
if (flagall == 0) return ETOL;
}
}
// magnetic torque tolerance criterion
// sync across replicas if running multi-replica minimization
fmdotfm = fmsq = 0.0;
if (update->ftol > 0.0) {
if (normstyle == MAX) fmsq = max_torque(); // max torque norm
else if (normstyle == INF) fmsq = inf_torque(); // inf torque norm
else if (normstyle == TWO) fmsq = total_torque(); // Euclidean torque 2-norm
else error->all(FLERR,"Illegal min_modify command");
fmdotfm = fmsq*fmsq;
if (update->multireplica == 0) {
if (fmdotfm < update->ftol*update->ftol) return FTOL;
} else {
if (fmdotfm < update->ftol*update->ftol) flag = 0;
else flag = 1;
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld);
if (flagall == 0) return FTOL;
}
}
// output for thermo, dump, restart files
if (output->next == ntimestep) {
timer->stamp();
output->write(ntimestep);
timer->stamp(Timer::OUTPUT);
}
}
return MAXITER;
}
/* ----------------------------------------------------------------------
calculate gradients
---------------------------------------------------------------------- */
void MinSpinLBFGS::calc_gradient()
{
int nlocal = atom->nlocal;
double **sp = atom->sp;
double **fm = atom->fm;
double hbar = force->hplanck/MY_2PI;
// loop on all spins on proc.
for (int i = 0; i < nlocal; i++) {
g_cur[3 * i + 0] = (fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]) * hbar;
g_cur[3 * i + 1] = -(fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]) * hbar;
g_cur[3 * i + 2] = (fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]) * hbar;
}
}
/* ----------------------------------------------------------------------
search direction:
Limited-memory BFGS.
See Jorge Nocedal and Stephen J. Wright 'Numerical
Optimization' Second Edition, 2006 (p. 177)
---------------------------------------------------------------------- */
void MinSpinLBFGS::calc_search_direction()
{
int nlocal = atom->nlocal;
double dyds = 0.0;
double sq = 0.0;
double yy = 0.0;
double yr = 0.0;
double beta = 0.0;
double dyds_global = 0.0;
double sq_global = 0.0;
double yy_global = 0.0;
double yr_global = 0.0;
int m_index = local_iter % num_mem; // memory index
int c_ind = 0;
double *q;
double *alpha;
double factor;
double scaling = 1.0;
// for multiple replica do not move end points
if (nreplica > 1) {
if (ireplica == 0 || ireplica == nreplica - 1) {
factor = 0.0;
}
else factor = 1.0;
}else{
factor = 1.0;
}
if (local_iter == 0){ // steepest descent direction
//if no line search then calculate maximum rotation
if (use_line_search == 0)
scaling = maximum_rotation(g_cur);
for (int i = 0; i < 3 * nlocal; i++) {
p_s[i] = -g_cur[i] * factor * scaling;;
g_old[i] = g_cur[i] * factor;
for (int k = 0; k < num_mem; k++){
ds[k][i] = 0.0;
dy[k][i] = 0.0;
}
}
for (int k = 0; k < num_mem; k++)
rho[k] = 0.0;
} else {
dyds = 0.0;
for (int i = 0; i < 3 * nlocal; i++) {
ds[m_index][i] = p_s[i];
dy[m_index][i] = g_cur[i] - g_old[i];
dyds += ds[m_index][i] * dy[m_index][i];
}
MPI_Allreduce(&dyds, &dyds_global, 1, MPI_DOUBLE, MPI_SUM, world);
if (nreplica > 1) {
dyds_global *= factor;
dyds = dyds_global;
MPI_Allreduce(&dyds, &dyds_global, 1,MPI_DOUBLE,MPI_SUM,universe->uworld);
}
if (fabs(dyds_global) > 1.0e-60) rho[m_index] = 1.0 / dyds_global;
else rho[m_index] = 1.0e60;
if (rho[m_index] < 0.0){
local_iter = 0;
return calc_search_direction();
}
q = (double *) calloc(3*nlocal, sizeof(double));
alpha = (double *) calloc(num_mem, sizeof(double));
// set the q vector
for (int i = 0; i < 3 * nlocal; i++) {
q[i] = g_cur[i];
}
// loop over last m indecies
for(int k = num_mem - 1; k > -1; k--) {
// this loop should run from the newest memory to the oldest one.
c_ind = (k + m_index + 1) % num_mem;
// dot product between dg and q
sq = 0.0;
for (int i = 0; i < 3 * nlocal; i++) {
sq += ds[c_ind][i] * q[i];
}
MPI_Allreduce(&sq,&sq_global,1,MPI_DOUBLE,MPI_SUM,world);
if (nreplica > 1) {
sq_global *= factor;
sq = sq_global;
MPI_Allreduce(&sq,&sq_global,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
}
// update alpha
alpha[c_ind] = rho[c_ind] * sq_global;
// update q
for (int i = 0; i < 3 * nlocal; i++) {
q[i] -= alpha[c_ind] * dy[c_ind][i];
}
}
// dot product between dg with itself
yy = 0.0;
for (int i = 0; i < 3 * nlocal; i++) {
yy += dy[m_index][i] * dy[m_index][i];
}
MPI_Allreduce(&yy,&yy_global,1,MPI_DOUBLE,MPI_SUM,world);
if (nreplica > 1) {
yy_global *= factor;
yy = yy_global;
MPI_Allreduce(&yy,&yy_global,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
}
// calculate now search direction
double devis = rho[m_index] * yy_global;
if (fabs(devis) > 1.0e-60) {
for (int i = 0; i < 3 * nlocal; i++) {
p_s[i] = factor * q[i] / devis;
}
}else{
for (int i = 0; i < 3 * nlocal; i++) {
p_s[i] = factor * q[i] * 1.0e60;
}
}
for (int k = 0; k < num_mem; k++){
// this loop should run from the oldest memory to the newest one.
if (local_iter < num_mem) c_ind = k;
else c_ind = (k + m_index + 1) % num_mem;
// dot product between p and da
yr = 0.0;
for (int i = 0; i < 3 * nlocal; i++) {
yr += dy[c_ind][i] * p_s[i];
}
MPI_Allreduce(&yr,&yr_global,1,MPI_DOUBLE,MPI_SUM,world);
if (nreplica > 1) {
yr_global *= factor;
yr = yr_global;
MPI_Allreduce(&yr,&yr_global,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
}
beta = rho[c_ind] * yr_global;
for (int i = 0; i < 3 * nlocal; i++) {
p_s[i] += ds[c_ind][i] * (alpha[c_ind] - beta);
}
}
if (use_line_search == 0)
scaling = maximum_rotation(p_s);
for (int i = 0; i < 3 * nlocal; i++) {
p_s[i] = - factor * p_s[i] * scaling;
g_old[i] = g_cur[i] * factor;
}
free(q);
free(alpha);
}
local_iter++;
}
/* ----------------------------------------------------------------------
rotation of spins along the search direction
---------------------------------------------------------------------- */
void MinSpinLBFGS::advance_spins()
{
int nlocal = atom->nlocal;
double **sp = atom->sp;
double rot_mat[9]; // exponential of matrix made of search direction
double s_new[3];
// loop on all spins on proc.
for (int i = 0; i < nlocal; i++) {
rodrigues_rotation(p_s + 3 * i, rot_mat);
// rotate spins
vm3(rot_mat, sp[i], s_new);
for (int j = 0; j < 3; j++) sp[i][j] = s_new[j];
}
}
/* ----------------------------------------------------------------------
calculate 3x3 matrix exponential using Rodrigues' formula
(R. Murray, Z. Li, and S. Shankar Sastry,
A Mathematical Introduction to
Robotic Manipulation (1994), p. 28 and 30).
upp_tr - vector x, y, z so that one calculate
U = exp(A) with A= [[0, x, y],
[-x, 0, z],
[-y, -z, 0]]
------------------------------------------------------------------------- */
void MinSpinLBFGS::rodrigues_rotation(const double *upp_tr, double *out)
{
double theta,A,B,D,x,y,z;
double s1,s2,s3,a1,a2,a3;
if (fabs(upp_tr[0]) < 1.0e-40 &&
fabs(upp_tr[1]) < 1.0e-40 &&
fabs(upp_tr[2]) < 1.0e-40){
// if upp_tr is zero, return unity matrix
for(int k = 0; k < 3; k++){
for(int m = 0; m < 3; m++){
if (m == k) out[3 * k + m] = 1.0;
else out[3 * k + m] = 0.0;
}
}
return;
}
theta = sqrt(upp_tr[0] * upp_tr[0] +
upp_tr[1] * upp_tr[1] +
upp_tr[2] * upp_tr[2]);
A = cos(theta);
B = sin(theta);
D = 1 - A;
x = upp_tr[0]/theta;
y = upp_tr[1]/theta;
z = upp_tr[2]/theta;
// diagonal elements of U
out[0] = A + z * z * D;
out[4] = A + y * y * D;
out[8] = A + x * x * D;
// off diagonal of U
s1 = -y * z *D;
s2 = x * z * D;
s3 = -x * y * D;
a1 = x * B;
a2 = y * B;
a3 = z * B;
out[1] = s1 + a1;
out[3] = s1 - a1;
out[2] = s2 + a2;
out[6] = s2 - a2;
out[5] = s3 + a3;
out[7] = s3 - a3;
}
/* ----------------------------------------------------------------------
out = vector^T x m,
m -- 3x3 matrix , v -- 3-d vector
------------------------------------------------------------------------- */
void MinSpinLBFGS::vm3(const double *m, const double *v, double *out)
{
for(int i = 0; i < 3; i++){
out[i] = 0.0;
for(int j = 0; j < 3; j++)
out[i] += *(m + 3 * j + i) * v[j];
}
}
void MinSpinLBFGS::make_step(double c, double *energy_and_der)
{
double p_scaled[3];
int nlocal = atom->nlocal;
double rot_mat[9]; // exponential of matrix made of search direction
double s_new[3];
double **sp = atom->sp;
double der_e_cur_tmp = 0.0;
for (int i = 0; i < nlocal; i++) {
// scale the search direction
for (int j = 0; j < 3; j++) p_scaled[j] = c * p_s[3 * i + j];
// calculate rotation matrix
rodrigues_rotation(p_scaled, rot_mat);
// rotate spins
vm3(rot_mat, sp[i], s_new);
for (int j = 0; j < 3; j++) sp[i][j] = s_new[j];
}
ecurrent = energy_force(0);
calc_gradient();
neval++;
der_e_cur = 0.0;
for (int i = 0; i < 3 * nlocal; i++) {
der_e_cur += g_cur[i] * p_s[i];
}
MPI_Allreduce(&der_e_cur,&der_e_cur_tmp, 1, MPI_DOUBLE, MPI_SUM, world);
der_e_cur = der_e_cur_tmp;
if (update->multireplica == 1) {
MPI_Allreduce(&der_e_cur_tmp,&der_e_cur,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
}
energy_and_der[0] = ecurrent;
energy_and_der[1] = der_e_cur;
}
/* ----------------------------------------------------------------------
Calculate step length which satisfies approximate Wolfe conditions
using the cubic interpolation
------------------------------------------------------------------------- */
int MinSpinLBFGS::calc_and_make_step(double a, double b, int index)
{
double e_and_d[2] = {0.0,0.0};
double alpha,c1,c2,c3;
double **sp = atom->sp;
int nlocal = atom->nlocal;
make_step(b,e_and_d);
ecurrent = e_and_d[0];
der_e_cur = e_and_d[1];
index++;
if (adescent(eprevious,e_and_d[0]) || index == 5){
MPI_Bcast(&b,1,MPI_DOUBLE,0,world);
for (int i = 0; i < 3 * nlocal; i++) {
p_s[i] = b * p_s[i];
}
return 1;
}
else{
double r,f0,f1,df0,df1;
r = b - a;
f0 = eprevious;
f1 = ecurrent;
df0 = der_e_pr;
df1 = der_e_cur;
c1 = -2.0*(f1-f0)/(r*r*r)+(df1+df0)/(r*r);
c2 = 3.0*(f1-f0)/(r*r)-(df1+2.0*df0)/(r);
c3 = df0;
// f(x) = c1 x^3 + c2 x^2 + c3 x^1 + c4
// has minimum at alpha below. We do not check boundaries.
alpha = (-c2 + sqrt(c2*c2 - 3.0*c1*c3))/(3.0*c1);
MPI_Bcast(&alpha,1,MPI_DOUBLE,0,world);
if (alpha < 0.0) alpha = r/2.0;
for (int i = 0; i < nlocal; i++) {
for (int j = 0; j < 3; j++) sp[i][j] = sp_copy[i][j];
}
calc_and_make_step(0.0, alpha, index);
}
return 0;
}
/* ----------------------------------------------------------------------
Approximate descent
------------------------------------------------------------------------- */
int MinSpinLBFGS::adescent(double phi_0, double phi_j){
double eps = 1.0e-6;
if (phi_j<=phi_0+eps*fabs(phi_0))
return 1;
else
return 0;
}
double MinSpinLBFGS::maximum_rotation(double *p)
{
double norm2,norm2_global,scaling,alpha;
int nlocal = atom->nlocal;
int ntotal = 0;
norm2 = 0.0;
for (int i = 0; i < 3 * nlocal; i++) norm2 += p[i] * p[i];
MPI_Allreduce(&norm2,&norm2_global,1,MPI_DOUBLE,MPI_SUM,world);
if (nreplica > 1) {
norm2 = norm2_global;
MPI_Allreduce(&norm2,&norm2_global,1,MPI_DOUBLE,MPI_SUM,universe->uworld);
}
MPI_Allreduce(&nlocal,&ntotal,1,MPI_INT,MPI_SUM,world);
if (nreplica > 1) {
nlocal = ntotal;
MPI_Allreduce(&nlocal,&ntotal,1,MPI_INT,MPI_SUM,universe->uworld);
}
scaling = (maxepsrot * sqrt((double) ntotal / norm2_global));
if (scaling < 1.0) alpha = scaling;
else alpha = 1.0;
return alpha;
}

72
src/SPIN/min_spin_lbfgs.h Normal file
View File

@ -0,0 +1,72 @@
/* -*- 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 MINIMIZE_CLASS
MinimizeStyle(spin/lbfgs, MinSpinLBFGS)
#else
#ifndef LMP_MIN_SPIN_LBFGS_H
#define LMP_MIN_SPIN_LBFGS_H
#include "min.h"
namespace LAMMPS_NS {
class MinSpinLBFGS: public Min {
public:
MinSpinLBFGS(class LAMMPS *);
virtual ~MinSpinLBFGS();
void init();
void setup_style();
int modify_param(int, char **);
void reset_vectors();
int iterate(int);
private:
int local_iter; // for neb
int use_line_search; // use line search or not.
int nlocal_max; // max value of nlocal (for size of lists)
int ireplica,nreplica; // for neb
double der_e_cur; // current derivative along search dir.
double der_e_pr; // previous derivative along search dir.
double maxepsrot;
double *spvec; // variables for atomic dof, as 1d vector
double *fmvec; // variables for atomic dof, as 1d vector
double *g_old; // gradient vector at previous step
double *g_cur; // current gradient vector
double *p_s; // search direction vector
void advance_spins();
void calc_gradient();
void calc_search_direction();
void vm3(const double *, const double *, double *);
void rodrigues_rotation(const double *, double *);
void make_step(double, double *);
int calc_and_make_step(double, double, int);
int adescent(double, double);
double maximum_rotation(double *);
double *rho; // estimation of curvature
double **ds; // change in rotation matrix between two iterations, da
double **dy; // change in gradients between two iterations, dg
double **sp_copy; // copy of the spins
int num_mem; // number of stored steps
bigint last_negative;
};
}
#endif
#endif

10
src/SPIN/neb_spin.cpp
View File

@ -43,6 +43,8 @@
#include "timer.h"
#include "memory.h"
#include "error.h"
#include "math_const.h"
#include "utils.h"
using namespace LAMMPS_NS;
@ -186,8 +188,8 @@ void NEBSpin::run()
if (update->minimize->searchflag)
error->all(FLERR,"NEBSpin requires damped dynamics minimizer");
if (strcmp(update->minimize_style,"spin") != 0)
error->all(FLERR,"NEBSpin requires spin minimizer");
if (!utils::strmatch(update->minimize_style,"^spin"))
error->all(FLERR,"NEBSpin requires a spin minimizer");
// setup regular NEBSpin minimization
@ -243,6 +245,8 @@ void NEBSpin::run()
timer->init();
timer->barrier_start();
// if(ireplica != 0 && ireplica != nreplica -1)
while (update->minimize->niter < n1steps) {
update->minimize->run(nevery);
print_status();
@ -639,7 +643,7 @@ int NEBSpin::initial_rotation(double *spi, double *sploc, double fraction)
kcrossy = kz*spix - kx*spiz;
kcrossz = kx*spiy - ky*spix;
kdots = kx*spix + ky*spiz + kz*spiz;
kdots = kx*spix + ky*spiy + kz*spiz;
omega = acos(sidotsf);
omega *= fraction;

2
src/SPIN/pair_spin_dipole_cut.cpp
View File

@ -323,7 +323,7 @@ void PairSpinDipoleCut::compute(int eflag, int vflag)
void PairSpinDipoleCut::compute_single_pair(int ii, double fmi[3])
{
int j,jnum,itype,jtype,ntypes;
int *ilist,*jlist,*numneigh,**firstneigh;
int *jlist,*numneigh,**firstneigh;
double rsq,rinv,r2inv,r3inv,local_cut2;
double xi[3],rij[3],eij[3],spi[4],spj[4];

5
src/SPIN/pair_spin_dipole_long.cpp
View File

@ -354,10 +354,9 @@ void PairSpinDipoleLong::compute(int eflag, int vflag)
void PairSpinDipoleLong::compute_single_pair(int ii, double fmi[3])
{
//int i,j,jj,jnum,itype,jtype;
int j,jj,jnum,itype,jtype,ntypes;
int k,locflag;
int *ilist,*jlist,*numneigh,**firstneigh;
int *jlist,*numneigh,**firstneigh;
double r,rinv,r2inv,rsq,grij,expm2,t,erfc;
double local_cut2,pre1,pre2,pre3;
double bij[4],xi[3],rij[3],eij[3],spi[4],spj[4];
@ -367,7 +366,6 @@ void PairSpinDipoleLong::compute_single_pair(int ii, double fmi[3])
double **sp = atom->sp;
double **fm_long = atom->fm_long;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
@ -405,7 +403,6 @@ void PairSpinDipoleLong::compute_single_pair(int ii, double fmi[3])
// computation of the exchange interaction
// loop over neighbors of atom i
//i = ilist[ii];
xi[0] = x[ii][0];
xi[1] = x[ii][1];
xi[2] = x[ii][2];

8
src/SPIN/pair_spin_dmi.cpp
View File

@ -313,7 +313,7 @@ void PairSpinDmi::compute(int eflag, int vflag)
if (eflag) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
evdwl *= 0.5*hbar;
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
@ -427,9 +427,9 @@ void PairSpinDmi::compute_dmi(int i, int j, double eij[3], double fmi[3], double
dmiy = eij[2]*v_dmx[itype][jtype] - eij[0]*v_dmz[itype][jtype];
dmiz = eij[0]*v_dmy[itype][jtype] - eij[1]*v_dmx[itype][jtype];
fmi[0] -= (dmiy*spj[2] - dmiz*spj[1]);
fmi[1] -= (dmiz*spj[0] - dmix*spj[2]);
fmi[2] -= (dmix*spj[1] - dmiy*spj[0]);
fmi[0] -= 2.0*(dmiy*spj[2] - dmiz*spj[1]);
fmi[1] -= 2.0*(dmiz*spj[0] - dmix*spj[2]);
fmi[2] -= 2.0*(dmix*spj[1] - dmiy*spj[0]);
}
/* ----------------------------------------------------------------------

10
src/SPIN/pair_spin_exchange.cpp
View File

@ -296,7 +296,7 @@ void PairSpinExchange::compute(int eflag, int vflag)
if (eflag) {
evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
evdwl *= hbar;
evdwl *= 0.5*hbar;
} else evdwl = 0.0;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
@ -405,9 +405,9 @@ void PairSpinExchange::compute_exchange(int i, int j, double rsq, double fmi[3],
Jex *= (1.0-J2[itype][jtype]*ra);
Jex *= exp(-ra);
fmi[0] += Jex*spj[0];
fmi[1] += Jex*spj[1];
fmi[2] += Jex*spj[2];
fmi[0] += 2.0*Jex*spj[0];
fmi[1] += 2.0*Jex*spj[1];
fmi[2] += 2.0*Jex*spj[2];
}
/* ----------------------------------------------------------------------
@ -504,7 +504,7 @@ void PairSpinExchange::read_restart(FILE *fp)
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(&J3[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);

6
src/SPIN/pair_spin_neel.cpp
View File

@ -643,10 +643,8 @@ void PairSpinNeel::allocate()
memory->create(q3,n+1,n+1,"pair/spin/soc/neel:q3");
memory->create(cutsq,n+1,n+1,"pair/spin/soc/neel:cutsq");
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
@ -694,11 +692,11 @@ void PairSpinNeel::read_restart(FILE *fp)
fread(&g1[i][j],sizeof(double),1,fp);
fread(&g1_mech[i][j],sizeof(double),1,fp);
fread(&g2[i][j],sizeof(double),1,fp);
fread(&g2[i][j],sizeof(double),1,fp);
fread(&g3[i][j],sizeof(double),1,fp);
fread(&q1[i][j],sizeof(double),1,fp);
fread(&q1_mech[i][j],sizeof(double),1,fp);
fread(&q2[i][j],sizeof(double),1,fp);
fread(&q2[i][j],sizeof(double),1,fp);
fread(&q3[i][j],sizeof(double),1,fp);
fread(&cut_spin_neel[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&g1[i][j],1,MPI_DOUBLE,0,world);

112
src/USER-MISC/fix_bond_react.cpp
View File

@ -35,6 +35,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
#include "molecule.h"
#include "group.h"
#include "citeme.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
@ -42,6 +43,7 @@ Contributing Author: Jacob Gissinger (jacob.gissinger@colorado.edu)
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
static const char cite_fix_bond_react[] =
"fix bond/react:\n\n"
@ -57,7 +59,7 @@ static const char cite_fix_bond_react[] =
#define BIG 1.0e20
#define DELTA 16
#define MAXGUESS 20 // max # of guesses allowed by superimpose algorithm
#define MAXCONARGS 5 // max # of arguments for any type of constraint
#define MAXCONARGS 7 // max # of arguments for any type of constraint + rxnID
// various statuses of superimpose algorithm:
// ACCEPT: site successfully matched to pre-reacted template
@ -68,6 +70,9 @@ static const char cite_fix_bond_react[] =
// RESTORE: restore mode, load most recent restore point
enum{ACCEPT,REJECT,PROCEED,CONTINUE,GUESSFAIL,RESTORE};
// types of available reaction constraints
enum{DISTANCE,ANGLE};
/* ---------------------------------------------------------------------- */
FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
@ -94,6 +99,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
global_freq = 1;
extvector = 0;
rxnID = 0;
nconstraints = 0;
status = PROCEED;
nxspecial = NULL;
@ -169,8 +175,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(limit_duration,nreacts,"bond/react:limit_duration");
memory->create(stabilize_steps_flag,nreacts,"bond/react:stabilize_steps_flag");
memory->create(update_edges_flag,nreacts,"bond/react:update_edges_flag");
memory->create(nconstraints,nreacts,"bond/react:nconstraints");
memory->create(constraints,nreacts,MAXCONARGS,"bond/react:constraints");
memory->create(constraints,1,MAXCONARGS,"bond/react:constraints");
memory->create(iatomtype,nreacts,"bond/react:iatomtype");
memory->create(jatomtype,nreacts,"bond/react:jatomtype");
memory->create(ibonding,nreacts,"bond/react:ibonding");
@ -188,7 +193,6 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
max_rxn[i] = INT_MAX;
stabilize_steps_flag[i] = 0;
update_edges_flag[i] = 0;
nconstraints[i] = 0;
// set default limit duration to 60 timesteps
limit_duration[i] = 60;
reaction_count[i] = 0;
@ -1138,6 +1142,22 @@ void FixBondReact::superimpose_algorithm()
glove[myjbonding-1][1] = created[lcl_inst][1][rxnID];
glove_counter++;
// special case, only two atoms in reaction templates
// then: bonding onemol_nxspecials guaranteed to be equal, and either 0 or 1
if (glove_counter == onemol->natoms) {
tagint local_atom1 = atom->map(glove[myibonding-1][1]);
tagint local_atom2 = atom->map(glove[myjbonding-1][1]);
if ( (nxspecial[local_atom1][0] == onemol_nxspecial[myibonding-1][0] &&
nxspecial[local_atom2][0] == nxspecial[local_atom1][0]) &&
(nxspecial[local_atom1][0] == 0 ||
xspecial[local_atom1][0] == atom->tag[local_atom2]) &&
check_constraints() ) {
status = ACCEPT;
glove_ghostcheck();
} else
status = REJECT;
}
avail_guesses = 0;
for (int i = 0; i < max_natoms; i++)
@ -1617,21 +1637,50 @@ evaluate constraints: return 0 if any aren't satisfied
int FixBondReact::check_constraints()
{
tagint atom1,atom2;
tagint atom1,atom2,atom3;
double delx,dely,delz,rsq;
double delx1,dely1,delz1,delx2,dely2,delz2;
double rsq1,rsq2,r1,r2,c;
double **x = atom->x;
for (int i = 0; i < nconstraints[rxnID]; i++) {
if (constraints[rxnID][0] == 0) { // 'distance' type
atom1 = atom->map(glove[(int) constraints[rxnID][1]-1][1]);
atom2 = atom->map(glove[(int) constraints[rxnID][2]-1][1]);
delx = x[atom1][0] - x[atom2][0];
dely = x[atom1][1] - x[atom2][1];
delz = x[atom1][2] - x[atom2][2];
domain->minimum_image(delx,dely,delz); // ghost location fix
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < constraints[rxnID][3] || rsq > constraints[rxnID][4]) return 0;
for (int i = 0; i < nconstraints; i++) {
if (constraints[i][0] == rxnID) {
if (constraints[i][1] == DISTANCE) {
atom1 = atom->map(glove[(int) constraints[i][2]-1][1]);
atom2 = atom->map(glove[(int) constraints[i][3]-1][1]);
delx = x[atom1][0] - x[atom2][0];
dely = x[atom1][1] - x[atom2][1];
delz = x[atom1][2] - x[atom2][2];
domain->minimum_image(delx,dely,delz); // ghost location fix
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < constraints[i][4] || rsq > constraints[i][5]) return 0;
} else if (constraints[i][1] == ANGLE) {
atom1 = atom->map(glove[(int) constraints[i][2]-1][1]);
atom2 = atom->map(glove[(int) constraints[i][3]-1][1]);
atom3 = atom->map(glove[(int) constraints[i][4]-1][1]);
// 1st bond
delx1 = x[atom1][0] - x[atom2][0];
dely1 = x[atom1][1] - x[atom2][1];
delz1 = x[atom1][2] - x[atom2][2];
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
r1 = sqrt(rsq1);
// 2nd bond
delx2 = x[atom3][0] - x[atom2][0];
dely2 = x[atom3][1] - x[atom2][1];
delz2 = x[atom3][2] - x[atom2][2];
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
r2 = sqrt(rsq2);
// angle (cos and sin)
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
c /= r1*r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
if (acos(c) < constraints[i][5] || acos(c) > constraints[i][6]) return 0;
}
}
}
return 1;
@ -2757,13 +2806,20 @@ void FixBondReact::read(int myrxn)
if (strspn(line," \t\n\r") == strlen(line)) continue;
if (strstr(line,"edgeIDs")) sscanf(line,"%d",&nedge);
else if (strstr(line,"equivalences")) sscanf(line,"%d",&nequivalent);
else if (strstr(line,"equivalences")) {
sscanf(line,"%d",&nequivalent);
if (nequivalent != onemol->natoms)
error->one(FLERR,"Bond/react: Number of equivalences in map file must "
"equal number of atoms in reaction templates");
}
else if (strstr(line,"customIDs")) sscanf(line,"%d",&ncustom);
else if (strstr(line,"deleteIDs")) sscanf(line,"%d",&ndelete);
else if (strstr(line,"constraints")) sscanf(line,"%d",&nconstraints[myrxn]);
else if (strstr(line,"constraints")) sscanf(line,"%d",&nconstr);
else break;
}
memory->grow(constraints,nconstraints+nconstr,MAXCONARGS,"bond/react:constraints");
//count = NULL;
// grab keyword and skip next line
@ -2874,18 +2930,28 @@ void FixBondReact::Constraints(char *line, int myrxn)
double tmp[MAXCONARGS];
int n = strlen("distance") + 1;
char *constraint_type = new char[n];
for (int i = 0; i < nconstraints[myrxn]; i++) {
for (int i = 0; i < nconstr; i++) {
readline(line);
sscanf(line,"%s",constraint_type);
constraints[nconstraints][0] = myrxn;
if (strcmp(constraint_type,"distance") == 0) {
constraints[myrxn][0] = 0; // 0 = 'distance' ...maybe use another enum eventually
constraints[nconstraints][1] = DISTANCE;
sscanf(line,"%*s %lg %lg %lg %lg",&tmp[0],&tmp[1],&tmp[2],&tmp[3]);
constraints[myrxn][1] = tmp[0];
constraints[myrxn][2] = tmp[1];
constraints[myrxn][3] = tmp[2]*tmp[2]; // using square of distance
constraints[myrxn][4] = tmp[3]*tmp[3];
constraints[nconstraints][2] = tmp[0];
constraints[nconstraints][3] = tmp[1];
constraints[nconstraints][4] = tmp[2]*tmp[2]; // using square of distance
constraints[nconstraints][5] = tmp[3]*tmp[3];
} else if (strcmp(constraint_type,"angle") == 0) {
constraints[nconstraints][1] = ANGLE;
sscanf(line,"%*s %lg %lg %lg %lg %lg",&tmp[0],&tmp[1],&tmp[2],&tmp[3],&tmp[4]);
constraints[nconstraints][2] = tmp[0];
constraints[nconstraints][3] = tmp[1];
constraints[nconstraints][4] = tmp[2];
constraints[nconstraints][5] = tmp[3]/180.0 * MY_PI;
constraints[nconstraints][6] = tmp[4]/180.0 * MY_PI;
} else
error->one(FLERR,"Bond/react: Illegal constraint type in 'Constraints' section of map file");
nconstraints++;
}
delete [] constraint_type;
}

4
src/USER-MISC/fix_bond_react.h
View File

@ -64,7 +64,7 @@ class FixBondReact : public Fix {
int custom_exclude_flag;
int *stabilize_steps_flag;
int *update_edges_flag;
int *nconstraints;
int nconstraints;
double **constraints;
int status;
int *groupbits;
@ -108,7 +108,7 @@ class FixBondReact : public Fix {
int *ibonding,*jbonding;
int *closeneigh; // indicates if bonding atoms of a rxn are 1-2, 1-3, or 1-4 neighbors
int nedge,nequivalent,ncustom,ndelete; // number of edge, equivalent, custom atoms in mapping file
int nedge,nequivalent,ncustom,ndelete,nconstr; // # edge, equivalent, custom atoms in mapping file
int attempted_rxn; // there was an attempt!
int *local_rxn_count;
int *ghostly_rxn_count;

347
src/fix_langevin.cpp
View File

@ -14,6 +14,9 @@
/* ----------------------------------------------------------------------
Contributing authors: Carolyn Phillips (U Mich), reservoir energy tally
Aidan Thompson (SNL) GJF formulation
Charles Sievers & Niels Gronbech-Jensen (UC Davis)
updated GJF formulation and included
statistically correct 2GJ velocity
------------------------------------------------------------------------- */
#include "fix_langevin.h"
@ -35,6 +38,7 @@
#include "memory.h"
#include "error.h"
#include "group.h"
#include "utils.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -50,7 +54,7 @@ enum{CONSTANT,EQUAL,ATOM};
FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
gjfflag(0), gfactor1(NULL), gfactor2(NULL), ratio(NULL), tstr(NULL),
flangevin(NULL), tforce(NULL), franprev(NULL), id_temp(NULL), random(NULL)
flangevin(NULL), tforce(NULL), franprev(NULL), lv(NULL), id_temp(NULL), random(NULL)
{
if (narg < 7) error->all(FLERR,"Illegal fix langevin command");
@ -96,6 +100,7 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
oflag = 0;
tallyflag = 0;
zeroflag = 0;
osflag = 0;
int iarg = 7;
while (iarg < narg) {
@ -106,8 +111,15 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
iarg += 2;
} else if (strcmp(arg[iarg],"gjf") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix langevin command");
if (strcmp(arg[iarg+1],"no") == 0) gjfflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) gjfflag = 1;
if (strcmp(arg[iarg+1],"no") == 0) {gjfflag = 0; osflag = 0;}
else if (strcmp(arg[iarg+1],"vfull") == 0) {
gjfflag = 1;
osflag = 1;
}
else if (strcmp(arg[iarg+1],"vhalf") == 0) {
gjfflag = 1;
osflag = 0;
}
else error->all(FLERR,"Illegal fix langevin command");
iarg += 2;
} else if (strcmp(arg[iarg],"omega") == 0) {
@ -153,6 +165,7 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
flangevin = NULL;
flangevin_allocated = 0;
franprev = NULL;
lv = NULL;
tforce = NULL;
maxatom1 = maxatom2 = 0;
@ -161,7 +174,6 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
// no need to set peratom_flag, b/c data is for internal use only
if (gjfflag) {
nvalues = 3;
grow_arrays(atom->nmax);
atom->add_callback(0);
@ -172,6 +184,9 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) :
franprev[i][0] = 0.0;
franprev[i][1] = 0.0;
franprev[i][2] = 0.0;
lv[i][0] = 0.0;
lv[i][1] = 0.0;
lv[i][2] = 0.0;
}
}
@ -192,6 +207,7 @@ FixLangevin::~FixLangevin()
if (gjfflag) {
memory->destroy(franprev);
memory->destroy(lv);
atom->delete_callback(id,0);
}
}
@ -201,6 +217,7 @@ FixLangevin::~FixLangevin()
int FixLangevin::setmask()
{
int mask = 0;
if (gjfflag) mask |= INITIAL_INTEGRATE;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= END_OF_STEP;
@ -212,6 +229,21 @@ int FixLangevin::setmask()
void FixLangevin::init()
{
if (gjfflag){
if (t_period*2 == update->dt)
error->all(FLERR,"Fix langevin gjf cannot have t_period equal to dt/2");
// warn if any integrate fix comes after this one
int before = 1;
int flag = 0;
for (int i = 0; i < modify->nfix; i++) {
if (strcmp(id,modify->fix[i]->id) == 0) before = 0;
else if ((modify->fmask[i] && utils::strmatch(modify->fix[i]->style,"^nve")) && before) flag = 1;
}
if (flag && comm->me == 0)
error->all(FLERR,"Fix langevin gjf should come before fix nve");
}
if (oflag && !atom->sphere_flag)
error->all(FLERR,"Fix langevin omega requires atom style sphere");
if (ascale && !atom->ellipsoid_flag)
@ -261,9 +293,14 @@ void FixLangevin::init()
if (!atom->rmass) {
for (int i = 1; i <= atom->ntypes; i++) {
gfactor1[i] = -atom->mass[i] / t_period / force->ftm2v;
gfactor2[i] = sqrt(atom->mass[i]) *
sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v;
if (gjfflag)
gfactor2[i] = sqrt(atom->mass[i]) *
sqrt(2.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v;
else
gfactor2[i] = sqrt(atom->mass[i]) *
sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
force->ftm2v;
gfactor1[i] *= 1.0/ratio[i];
gfactor2[i] *= 1.0/sqrt(ratio[i]);
}
@ -275,14 +312,60 @@ void FixLangevin::init()
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
if (gjfflag) gjffac = 1.0/(1.0+update->dt/2.0/t_period);
if (utils::strmatch(update->integrate_style,"^respa") && gjfflag)
error->all(FLERR,"Fix langevin gjf and respa are not compatible");
if (gjfflag) gjfa = (1.0-update->dt/2.0/t_period)/(1.0+update->dt/2.0/t_period);
if (gjfflag) gjfsib = sqrt(1.0+update->dt/2.0/t_period);
}
/* ---------------------------------------------------------------------- */
void FixLangevin::setup(int vflag)
{
if (gjfflag){
double dtfm;
double dt = update->dt;
double **v = atom->v;
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = 0.5 * dt / rmass[i];
v[i][0] -= dtfm * f[i][0];
v[i][1] -= dtfm * f[i][1];
v[i][2] -= dtfm * f[i][2];
if (tbiasflag)
temperature->remove_bias(i,v[i]);
v[i][0] /= gjfa*gjfsib*gjfsib;
v[i][1] /= gjfa*gjfsib*gjfsib;
v[i][2] /= gjfa*gjfsib*gjfsib;
if (tbiasflag)
temperature->restore_bias(i,v[i]);
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = 0.5 * dt / mass[type[i]];
v[i][0] -= dtfm * f[i][0];
v[i][1] -= dtfm * f[i][1];
v[i][2] -= dtfm * f[i][2];
if (tbiasflag)
temperature->remove_bias(i,v[i]);
v[i][0] /= gjfa*gjfsib*gjfsib;
v[i][1] /= gjfa*gjfsib*gjfsib;
v[i][2] /= gjfa*gjfsib*gjfsib;
if (tbiasflag)
temperature->restore_bias(i,v[i]);
}
}
}
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else {
@ -290,6 +373,61 @@ void FixLangevin::setup(int vflag)
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
if (gjfflag){
double dtfm;
double dt = update->dt;
double **f = atom->f;
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = 0.5 * dt / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
lv[i][0] = v[i][0];
lv[i][1] = v[i][1];
lv[i][2] = v[i][2];
}
//
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = 0.5 * dt / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
lv[i][0] = v[i][0];
lv[i][1] = v[i][1];
lv[i][2] = v[i][2];
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixLangevin::initial_integrate(int /* vflag */)
{
double **v = atom->v;
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
f[i][0] /= gjfa;
f[i][1] /= gjfa;
f[i][2] /= gjfa;
v[i][0] = lv[i][0];
v[i][1] = lv[i][1];
v[i][2] = lv[i][2];
}
}
/* ---------------------------------------------------------------------- */
@ -304,7 +442,7 @@ void FixLangevin::post_force(int /*vflag*/)
if (tstyle == ATOM)
if (gjfflag)
if (tallyflag)
if (tallyflag || osflag)
if (tbiasflag == BIAS)
if (rmass)
if (zeroflag) post_force_templated<1,1,1,1,1,1>();
@ -511,7 +649,10 @@ void FixLangevin::post_force_templated()
if (Tp_TSTYLEATOM) tsqrt = sqrt(tforce[i]);
if (Tp_RMASS) {
gamma1 = -rmass[i] / t_period / ftm2v;
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
if (Tp_GJF)
gamma2 = sqrt(rmass[i]) * sqrt(2.0*boltz/t_period/dt/mvv2e) / ftm2v;
else
gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
gamma1 *= 1.0/ratio[type[i]];
gamma2 *= 1.0/sqrt(ratio[type[i]]) * tsqrt;
} else {
@ -519,9 +660,16 @@ void FixLangevin::post_force_templated()
gamma2 = gfactor2[type[i]] * tsqrt;
}
fran[0] = gamma2*(random->uniform()-0.5);
fran[1] = gamma2*(random->uniform()-0.5);
fran[2] = gamma2*(random->uniform()-0.5);
if (Tp_GJF){
fran[0] = gamma2*random->gaussian();
fran[1] = gamma2*random->gaussian();
fran[2] = gamma2*random->gaussian();
}
else{
fran[0] = gamma2*(random->uniform()-0.5);
fran[1] = gamma2*(random->uniform()-0.5);
fran[2] = gamma2*(random->uniform()-0.5);
}
if (Tp_BIAS) {
temperature->remove_bias(i,v[i]);
@ -539,6 +687,16 @@ void FixLangevin::post_force_templated()
}
if (Tp_GJF) {
if (Tp_BIAS)
temperature->remove_bias(i,v[i]);
lv[i][0] = gjfsib*v[i][0];
lv[i][1] = gjfsib*v[i][1];
lv[i][2] = gjfsib*v[i][2];
if (Tp_BIAS)
temperature->restore_bias(i,v[i]);
if (Tp_BIAS)
temperature->restore_bias(i,lv[i]);
fswap = 0.5*(fran[0]+franprev[i][0]);
franprev[i][0] = fran[0];
fran[0] = fswap;
@ -549,32 +707,44 @@ void FixLangevin::post_force_templated()
franprev[i][2] = fran[2];
fran[2] = fswap;
fdrag[0] *= gjffac;
fdrag[1] *= gjffac;
fdrag[2] *= gjffac;
fran[0] *= gjffac;
fran[1] *= gjffac;
fran[2] *= gjffac;
f[i][0] *= gjffac;
f[i][1] *= gjffac;
f[i][2] *= gjffac;
fdrag[0] *= gjfa;
fdrag[1] *= gjfa;
fdrag[2] *= gjfa;
fran[0] *= gjfa;
fran[1] *= gjfa;
fran[2] *= gjfa;
f[i][0] *= gjfa;
f[i][1] *= gjfa;
f[i][2] *= gjfa;
}
f[i][0] += fdrag[0] + fran[0];
f[i][1] += fdrag[1] + fran[1];
f[i][2] += fdrag[2] + fran[2];
if (Tp_TALLY) {
flangevin[i][0] = fdrag[0] + fran[0];
flangevin[i][1] = fdrag[1] + fran[1];
flangevin[i][2] = fdrag[2] + fran[2];
}
if (Tp_ZERO) {
fsum[0] += fran[0];
fsum[1] += fran[1];
fsum[2] += fran[2];
}
if (Tp_TALLY) {
if (Tp_GJF){
fdrag[0] = gamma1*lv[i][0]/gjfsib/gjfsib;
fdrag[1] = gamma1*lv[i][1]/gjfsib/gjfsib;
fdrag[2] = gamma1*lv[i][2]/gjfsib/gjfsib;
fswap = (2*fran[0]/gjfa - franprev[i][0])/gjfsib;
fran[0] = fswap;
fswap = (2*fran[1]/gjfa - franprev[i][1])/gjfsib;
fran[1] = fswap;
fswap = (2*fran[2]/gjfa - franprev[i][2])/gjfsib;
fran[2] = fswap;
}
flangevin[i][0] = fdrag[0] + fran[0];
flangevin[i][1] = fdrag[1] + fran[1];
flangevin[i][2] = fdrag[2] + fran[2];
}
}
}
@ -754,18 +924,72 @@ void FixLangevin::angmom_thermostat()
void FixLangevin::end_of_step()
{
if (!tallyflag) return;
if (!tallyflag && !gjfflag) return;
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double dtfm;
double dt = update->dt;
double *mass = atom->mass;
double *rmass = atom->rmass;
double **f = atom->f;
int *type = atom->type;
energy_onestep = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2];
if (tallyflag){
if (gjfflag){
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (tbiasflag)
temperature->remove_bias(i, lv[i]);
energy_onestep += flangevin[i][0]*lv[i][0] + flangevin[i][1]*lv[i][1] +
flangevin[i][2]*lv[i][2];
if (tbiasflag)
temperature->restore_bias(i, lv[i]);
}
}
else
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2];
}
if (gjfflag){
double tmp[3];
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
tmp[0] = v[i][0];
tmp[1] = v[i][1];
tmp[2] = v[i][2];
if (!osflag){
v[i][0] = lv[i][0];
v[i][1] = lv[i][1];
v[i][2] = lv[i][2];
}
else{
if (atom->rmass) {
dtfm = force->ftm2v * 0.5 * dt / rmass[i];
} else {
dtfm = force->ftm2v * 0.5 * dt / mass[type[i]];
}
v[i][0] = 0.5 * gjfsib*gjfsib*(v[i][0] + dtfm * f[i][0] / gjfa) +
dtfm * 0.5 * (gjfsib * flangevin[i][0] - franprev[i][0]) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * lv[i][0];
v[i][1] = 0.5 * gjfsib*gjfsib*(v[i][1] + dtfm * f[i][1] / gjfa) +
dtfm * 0.5 * (gjfsib * flangevin[i][1] - franprev[i][1]) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * lv[i][1];
v[i][2] = 0.5 * gjfsib*gjfsib*(v[i][2] + dtfm * f[i][2] / gjfa) +
dtfm * 0.5 * (gjfsib * flangevin[i][2] - franprev[i][2]) +
(gjfsib * gjfa * 0.5 + dt * 0.25 / t_period / gjfsib) * lv[i][2];
}
lv[i][0] = tmp[0];
lv[i][1] = tmp[1];
lv[i][2] = tmp[2];
}
}
energy += energy_onestep*update->dt;
}
@ -831,11 +1055,25 @@ double FixLangevin::compute_scalar()
if (update->ntimestep == update->beginstep) {
energy_onestep = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2];
energy = 0.5*energy_onestep*update->dt;
if (!gjfflag){
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
energy_onestep += flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2];
energy = 0.5*energy_onestep*update->dt;
}
else{
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit){
if (tbiasflag)
temperature->remove_bias(i, lv[i]);
energy_onestep += flangevin[i][0]*lv[i][0] + flangevin[i][1]*lv[i][1] +
flangevin[i][2]*lv[i][2];
if (tbiasflag)
temperature->restore_bias(i, lv[i]);
}
energy = -0.5*energy_onestep*update->dt;
}
}
// convert midstep energy back to previous fullstep energy
@ -867,8 +1105,8 @@ void *FixLangevin::extract(const char *str, int &dim)
double FixLangevin::memory_usage()
{
double bytes = 0.0;
if (gjfflag) bytes += atom->nmax*3 * sizeof(double);
if (tallyflag) bytes += atom->nmax*3 * sizeof(double);
if (gjfflag) bytes += atom->nmax*6 * sizeof(double);
if (tallyflag || osflag) bytes += atom->nmax*3 * sizeof(double);
if (tforce) bytes += atom->nmax * sizeof(double);
return bytes;
}
@ -880,6 +1118,7 @@ double FixLangevin::memory_usage()
void FixLangevin::grow_arrays(int nmax)
{
memory->grow(franprev,nmax,3,"fix_langevin:franprev");
memory->grow(lv,nmax,3,"fix_langevin:lv");
}
/* ----------------------------------------------------------------------
@ -888,8 +1127,12 @@ void FixLangevin::grow_arrays(int nmax)
void FixLangevin::copy_arrays(int i, int j, int /*delflag*/)
{
for (int m = 0; m < nvalues; m++)
franprev[j][m] = franprev[i][m];
franprev[j][0] = franprev[i][0];
franprev[j][1] = franprev[i][1];
franprev[j][2] = franprev[i][2];
lv[j][0] = lv[i][0];
lv[j][1] = lv[i][1];
lv[j][2] = lv[i][2];
}
/* ----------------------------------------------------------------------
@ -898,8 +1141,14 @@ void FixLangevin::copy_arrays(int i, int j, int /*delflag*/)
int FixLangevin::pack_exchange(int i, double *buf)
{
for (int m = 0; m < nvalues; m++) buf[m] = franprev[i][m];
return nvalues;
int n = 0;
buf[n++] = franprev[i][0];
buf[n++] = franprev[i][1];
buf[n++] = franprev[i][2];
buf[n++] = lv[i][0];
buf[n++] = lv[i][1];
buf[n++] = lv[i][2];
return n;
}
/* ----------------------------------------------------------------------
@ -908,6 +1157,12 @@ int FixLangevin::pack_exchange(int i, double *buf)
int FixLangevin::unpack_exchange(int nlocal, double *buf)
{
for (int m = 0; m < nvalues; m++) franprev[nlocal][m] = buf[m];
return nvalues;
int n = 0;
franprev[nlocal][0] = buf[n++];
franprev[nlocal][1] = buf[n++];
franprev[nlocal][2] = buf[n++];
lv[nlocal][0] = buf[n++];
lv[nlocal][1] = buf[n++];
lv[nlocal][2] = buf[n++];
return n;
}

19
src/fix_langevin.h
View File

@ -31,6 +31,7 @@ class FixLangevin : public Fix {
int setmask();
void init();
void setup(int);
virtual void initial_integrate(int);
virtual void post_force(int);
void post_force_respa(int, int, int);
virtual void end_of_step();
@ -46,7 +47,7 @@ class FixLangevin : public Fix {
int unpack_exchange(int, double *);
protected:
int gjfflag,oflag,tallyflag,zeroflag,tbiasflag;
int gjfflag,nvalues,osflag,oflag,tallyflag,zeroflag,tbiasflag;
int flangevin_allocated;
double ascale;
double t_start,t_stop,t_period,t_target;
@ -54,7 +55,7 @@ class FixLangevin : public Fix {
double energy,energy_onestep;
double tsqrt;
int tstyle,tvar;
double gjffac;
double gjfa, gjfsib; //gjf a and gjf sqrt inverse b
char *tstr;
class AtomVecEllipsoid *avec;
@ -63,7 +64,7 @@ class FixLangevin : public Fix {
double **flangevin;
double *tforce;
double **franprev;
int nvalues;
double **lv; //half step velocity
char *id_temp;
class Compute *temperature;
@ -139,6 +140,18 @@ E: Fix_modify temperature ID does not compute temperature
The compute ID assigned to the fix must compute temperature.
E: Fix langevin gjf cannot have period equal to dt/2
If the period is equal to dt/2 then division by zero will happen.
E: Fix langevin gjf should come before fix nve
Self-explanatory
E: Fix langevin gjf and respa are not compatible
Self-explanatory
W: Group for fix_modify temp != fix group
The fix_modify command is specifying a temperature computation that

143
src/min.cpp
View File

@ -42,10 +42,12 @@
#include "output.h"
#include "thermo.h"
#include "timer.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
@ -54,6 +56,7 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
dmax = 0.1;
searchflag = 0;
linestyle = 1;
normstyle = TWO;
elist_global = elist_atom = NULL;
vlist_global = vlist_atom = NULL;
@ -659,6 +662,15 @@ void Min::modify_params(int narg, char **arg)
if (strcmp(arg[iarg+1],"backtrack") == 0) linestyle = 0;
else if (strcmp(arg[iarg+1],"quadratic") == 0) linestyle = 1;
else if (strcmp(arg[iarg+1],"forcezero") == 0) linestyle = 2;
else if (strcmp(arg[iarg+1],"spin_cubic") == 0) linestyle = 3;
else if (strcmp(arg[iarg+1],"spin_none") == 0) linestyle = 4;
else error->all(FLERR,"Illegal min_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"norm") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command");
if (strcmp(arg[iarg+1],"two") == 0) normstyle = TWO;
else if (strcmp(arg[iarg+1],"max") == 0) normstyle = MAX;
else if (strcmp(arg[iarg+1],"inf") == 0) normstyle = INF;
else error->all(FLERR,"Illegal min_modify command");
iarg += 2;
} else {
@ -822,6 +834,137 @@ double Min::fnorm_inf()
return norm_inf;
}
/* ----------------------------------------------------------------------
compute and return ||force||_max (inf norm per-vector)
------------------------------------------------------------------------- */
double Min::fnorm_max()
{
int i,n;
double fdotf,*fatom;
double local_norm_max = 0.0;
for (i = 0; i < nvec; i+=3) {
fdotf = fvec[i]*fvec[i]+fvec[i+1]*fvec[i+1]+fvec[i+2]*fvec[i+2];
local_norm_max = MAX(fdotf,local_norm_max);
}
if (nextra_atom) {
for (int m = 0; m < nextra_atom; m++) {
fatom = fextra_atom[m];
n = extra_nlen[m];
for (i = 0; i < n; i+=3)
fdotf = fvec[i]*fvec[i]+fvec[i+1]*fvec[i+1]+fvec[i+2]*fvec[i+2];
local_norm_max = MAX(fdotf,local_norm_max);
}
}
double norm_max = 0.0;
MPI_Allreduce(&local_norm_max,&norm_max,1,MPI_DOUBLE,MPI_MAX,world);
if (nextra_global)
for (i = 0; i < n; i+=3)
fdotf = fvec[i]*fvec[i]+fvec[i+1]*fvec[i+1]+fvec[i+2]*fvec[i+2];
norm_max = MAX(fdotf,norm_max);
return norm_max;
}
/* ----------------------------------------------------------------------
compute and return sum_i||mag. torque_i||_2 (in eV)
------------------------------------------------------------------------- */
double Min::total_torque()
{
double fmsq,ftotsqone,ftotsqall;
int nlocal = atom->nlocal;
double hbar = force->hplanck/MY_2PI;
double tx,ty,tz;
double **sp = atom->sp;
double **fm = atom->fm;
fmsq = ftotsqone = ftotsqall = 0.0;
for (int i = 0; i < nlocal; i++) {
tx = fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1];
ty = fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2];
tz = fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0];
fmsq = tx*tx + ty*ty + tz*tz;
ftotsqone += fmsq;
}
// summing all fmsqtot on this replica
MPI_Allreduce(&ftotsqone,&ftotsqall,1,MPI_DOUBLE,MPI_SUM,world);
// multiply it by hbar so that units are in eV
return sqrt(ftotsqall) * hbar;
}
/* ----------------------------------------------------------------------
compute and return max_i ||mag. torque components|| (in eV)
------------------------------------------------------------------------- */
double Min::inf_torque()
{
double fmsq,fmaxsqone,fmaxsqall;
int nlocal = atom->nlocal;
double hbar = force->hplanck/MY_2PI;
double tx,ty,tz;
double **sp = atom->sp;
double **fm = atom->fm;
fmsq = fmaxsqone = fmaxsqall = 0.0;
for (int i = 0; i < nlocal; i++) {
tx = fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1];
ty = fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2];
tz = fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0];
fmaxsqone = MAX(fmaxsqone,tx*tx);
fmaxsqone = MAX(fmaxsqone,ty*ty);
fmaxsqone = MAX(fmaxsqone,tz*tz);
}
// finding max fm on this replica
fmaxsqall = fmaxsqone;
MPI_Allreduce(&fmaxsqone,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,world);
// multiply it by hbar so that units are in eV
return sqrt(fmaxsqall) * hbar;
}
/* ----------------------------------------------------------------------
compute and return max_i ||mag. torque_i|| (in eV)
------------------------------------------------------------------------- */
double Min::max_torque()
{
double fmsq,fmaxsqone,fmaxsqall;
int nlocal = atom->nlocal;
double hbar = force->hplanck/MY_2PI;
double tx,ty,tz;
double **sp = atom->sp;
double **fm = atom->fm;
fmsq = fmaxsqone = fmaxsqall = 0.0;
for (int i = 0; i < nlocal; i++) {
tx = fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1];
ty = fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2];
tz = fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0];
fmsq = tx*tx + ty*ty + tz*tz;
fmaxsqone = MAX(fmaxsqone,fmsq);
}
// finding max fm on this replica
fmaxsqall = fmaxsqone;
MPI_Allreduce(&fmaxsqone,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,world);
// multiply it by hbar so that units are in eV
return sqrt(fmaxsqall) * hbar;
}
/* ----------------------------------------------------------------------
possible stop conditions
------------------------------------------------------------------------- */

22
src/min.h
View File

@ -39,8 +39,16 @@ class Min : protected Pointers {
virtual bigint memory_usage() {return 0;}
void modify_params(int, char **);
virtual int modify_param(int, char **) {return 0;}
virtual double fnorm_sqr();
virtual double fnorm_inf();
double fnorm_sqr();
double fnorm_inf();
double fnorm_max();
enum{TWO,MAX,INF};
// methods for spin minimizers
double total_torque();
double inf_torque();
double max_torque();
virtual void init_style() {}
virtual void setup_style() = 0;
@ -56,8 +64,11 @@ class Min : protected Pointers {
int virial_style; // compute virial explicitly or implicitly
int external_force_clear; // clear forces locally or externally
double dmax; // max dist to move any atom in one step
int linestyle; // 0 = backtrack, 1 = quadratic, 2 = forcezero
double dmax; // max dist to move any atom in one step
int linestyle; // 0 = backtrack, 1 = quadratic, 2 = forcezero
// 3 = spin_cubic, 4 = spin_none
int normstyle; // TWO, MAX or INF flag for force norm evaluation
int nelist_global,nelist_atom; // # of PE,virial computes to check
int nvlist_global,nvlist_atom;
@ -104,9 +115,6 @@ class Min : protected Pointers {
virtual double energy_force(int);
virtual void force_clear();
double compute_force_norm_sqr();
double compute_force_norm_inf();
void ev_setup();
void ev_set(bigint);

16
src/min_cg.cpp
View File

@ -14,6 +14,7 @@
#include "min_cg.h"
#include <mpi.h>
#include <cmath>
#include "error.h"
#include "update.h"
#include "output.h"
#include "timer.h"
@ -35,7 +36,7 @@ MinCG::MinCG(LAMMPS *lmp) : MinLineSearch(lmp) {}
int MinCG::iterate(int maxiter)
{
int i,m,n,fail,ntimestep;
double beta,gg,dot[2],dotall[2];
double beta,gg,dot[2],dotall[2],fmax;
double *fatom,*gatom,*hatom;
// nlimit = max # of CG iterations before restarting
@ -90,6 +91,7 @@ int MinCG::iterate(int maxiter)
dot[0] += fvec[i]*fvec[i];
dot[1] += fvec[i]*g[i];
}
if (nextra_atom)
for (m = 0; m < nextra_atom; m++) {
fatom = fextra_atom[m];
@ -98,6 +100,7 @@ int MinCG::iterate(int maxiter)
for (i = 0; i < n; i++) {
dot[0] += fatom[i]*fatom[i];
dot[1] += fatom[i]*gatom[i];
fmax = MAX(fmax,fatom[i]*fatom[i]);
}
}
MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world);
@ -107,7 +110,16 @@ int MinCG::iterate(int maxiter)
dotall[1] += fextra[i]*gextra[i];
}
if (dotall[0] < update->ftol*update->ftol) return FTOL;
fmax = 0.0;
if (normstyle == MAX) { // max force norm
fmax = fnorm_max();
if (fmax < update->ftol*update->ftol) return FTOL;
} else if (normstyle == INF) { // infinite force norm
fmax = fnorm_inf();
if (fmax < update->ftol*update->ftol) return FTOL;
} else if (normstyle == TWO) { // Euclidean force 2-norm
if (dotall[0] < update->ftol*update->ftol) return FTOL;
} else error->all(FLERR,"Illegal min_modify command");
// update new search direction h from new f = -Grad(x) and old g
// this is Polak-Ribieri formulation

8
src/min_fire.cpp
View File

@ -16,6 +16,7 @@
#include <cmath>
#include "universe.h"
#include "atom.h"
#include "error.h"
#include "force.h"
#include "update.h"
#include "output.h"
@ -80,7 +81,7 @@ void MinFire::reset_vectors()
int MinFire::iterate(int maxiter)
{
bigint ntimestep;
double vmax,vdotf,vdotfall,vdotv,vdotvall,fdotf,fdotfall;
double vmax,vdotf,vdotfall,vdotv,vdotvall,fdotf,fdotfloc,fdotfall;
double scale1,scale2;
double dtvone,dtv,dtf,dtfm;
int flag,flagall;
@ -250,7 +251,10 @@ int MinFire::iterate(int maxiter)
// sync across replicas if running multi-replica minimization
if (update->ftol > 0.0) {
fdotf = fnorm_sqr();
if (normstyle == MAX) fdotf = fnorm_max(); // max force norm
else if (normstyle == INF) fdotf = fnorm_inf(); // inf force norm
else if (normstyle == TWO) fdotf = fnorm_sqr(); // Euclidean force 2-norm
else error->all(FLERR,"Illegal min_modify command");
if (update->multireplica == 0) {
if (fdotf < update->ftol*update->ftol) return FTOL;
} else {

4
src/min_hftn.cpp
View File

@ -22,6 +22,7 @@
#include <cmath>
#include <cstring>
#include "atom.h"
#include "error.h"
#include "fix_minimize.h"
#include "modify.h"
#include "output.h"
@ -112,6 +113,9 @@ void MinHFTN::init()
{
Min::init();
if (normstyle == MAX)
error->all(FLERR,"Incorrect min_modify option");
for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) {
if (_daExtraGlobal[i] != NULL)
delete [] _daExtraGlobal[i];

8
src/min_quickmin.cpp
View File

@ -16,6 +16,7 @@
#include <cmath>
#include "universe.h"
#include "atom.h"
#include "error.h"
#include "force.h"
#include "update.h"
#include "output.h"
@ -75,7 +76,7 @@ void MinQuickMin::reset_vectors()
int MinQuickMin::iterate(int maxiter)
{
bigint ntimestep;
double vmax,vdotf,vdotfall,fdotf,fdotfall,scale;
double vmax,vdotf,vdotfall,fdotf,fdotfloc,fdotfall,scale;
double dtvone,dtv,dtf,dtfm;
int flag,flagall;
@ -215,7 +216,10 @@ int MinQuickMin::iterate(int maxiter)
// sync across replicas if running multi-replica minimization
if (update->ftol > 0.0) {
fdotf = fnorm_sqr();
if (normstyle == MAX) fdotfloc = fnorm_max(); // max force norm
else if (normstyle == INF) fdotfloc = fnorm_inf(); // inf force norm
else if (normstyle == TWO) fdotfloc = fnorm_sqr(); // Euclidean force 2-norm
else error->all(FLERR,"Illegal min_modify command");
if (update->multireplica == 0) {
if (fdotf < update->ftol*update->ftol) return FTOL;
} else {

8
src/min_sd.cpp
View File

@ -13,6 +13,7 @@
#include "min_sd.h"
#include <cmath>
#include "error.h"
#include "update.h"
#include "output.h"
#include "timer.h"
@ -34,7 +35,7 @@ MinSD::MinSD(LAMMPS *lmp) : MinLineSearch(lmp) {}
int MinSD::iterate(int maxiter)
{
int i,m,n,fail,ntimestep;
double fdotf;
double fdotf,fdotfloc;
double *fatom,*hatom;
// initialize working vectors
@ -77,7 +78,10 @@ int MinSD::iterate(int maxiter)
// force tolerance criterion
fdotf = fnorm_sqr();
if (normstyle == MAX) fdotf = fnorm_max(); // max force norm
else if (normstyle == INF) fdotf = fnorm_inf(); // infinite force norm
else if (normstyle == TWO) fdotf = fnorm_sqr(); // Euclidean force 2-norm
else error->all(FLERR,"Illegal min_modify command");
if (fdotf < update->ftol*update->ftol) return FTOL;
// set new search direction h to f = -Grad(x)

3
src/neighbor.cpp
View File

@ -1941,6 +1941,7 @@ int Neighbor::decide()
conservative shrink procedure:
compute distance each of 8 corners of box has moved since last reneighbor
reduce skin distance by sum of 2 largest of the 8 values
if reduced skin distance is negative, set to zero
new trigger = 1/2 of reduced skin distance
for orthogonal box, only need 2 lo/hi corners
for triclinic, need all 8 corners since deformations can displace all 8
@ -1962,6 +1963,7 @@ int Neighbor::check_distance()
delz = bboxhi[2] - boxhi_hold[2];
delta2 = sqrt(delx*delx + dely*dely + delz*delz);
delta = 0.5 * (skin - (delta1+delta2));
if (delta < 0.0) delta = 0.0;
deltasq = delta*delta;
} else {
domain->box_corners();
@ -1975,6 +1977,7 @@ int Neighbor::check_distance()
else if (delta > delta2) delta2 = delta;
}
delta = 0.5 * (skin - (delta1+delta2));
if (delta < 0.0) delta = 0.0;
deltasq = delta*delta;
}
} else deltasq = triggersq;