diff --git a/cmake/Modules/Packages/KOKKOS.cmake b/cmake/Modules/Packages/KOKKOS.cmake index cc1e051629..4134fed597 100644 --- a/cmake/Modules/Packages/KOKKOS.cmake +++ b/cmake/Modules/Packages/KOKKOS.cmake @@ -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 diff --git a/doc/src/Eqs/norm_inf.jpg b/doc/src/Eqs/norm_inf.jpg new file mode 100644 index 0000000000..42a2afb3d2 Binary files /dev/null and b/doc/src/Eqs/norm_inf.jpg differ diff --git a/doc/src/Eqs/norm_inf.tex b/doc/src/Eqs/norm_inf.tex new file mode 100644 index 0000000000..b97a9f88f7 --- /dev/null +++ b/doc/src/Eqs/norm_inf.tex @@ -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} diff --git a/doc/src/Eqs/norm_max.jpg b/doc/src/Eqs/norm_max.jpg new file mode 100644 index 0000000000..c10db9a531 Binary files /dev/null and b/doc/src/Eqs/norm_max.jpg differ diff --git a/doc/src/Eqs/norm_max.tex b/doc/src/Eqs/norm_max.tex new file mode 100644 index 0000000000..3b2198bdf0 --- /dev/null +++ b/doc/src/Eqs/norm_max.tex @@ -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} diff --git a/doc/src/Eqs/norm_two.jpg b/doc/src/Eqs/norm_two.jpg new file mode 100644 index 0000000000..5554de32f9 Binary files /dev/null and b/doc/src/Eqs/norm_two.jpg differ diff --git a/doc/src/Eqs/norm_two.tex b/doc/src/Eqs/norm_two.tex new file mode 100644 index 0000000000..d428081a49 --- /dev/null +++ b/doc/src/Eqs/norm_two.tex @@ -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} diff --git a/doc/src/Errors_messages.txt b/doc/src/Errors_messages.txt index c131b10ec6..33593d4d53 100644 --- a/doc/src/Errors_messages.txt +++ b/doc/src/Errors_messages.txt @@ -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 diff --git a/doc/src/Errors_warnings.txt b/doc/src/Errors_warnings.txt index 9f346ba8c1..fbd857f162 100644 --- a/doc/src/Errors_warnings.txt +++ b/doc/src/Errors_warnings.txt @@ -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 diff --git a/doc/src/fix_bond_react.txt b/doc/src/fix_bond_react.txt index 4f4506944f..fb8ed95afb 100644 --- a/doc/src/fix_bond_react.txt +++ b/doc/src/fix_bond_react.txt @@ -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. diff --git a/doc/src/fix_langevin.txt b/doc/src/fix_langevin.txt index 861eed4a6f..07d8b274aa 100644 --- a/doc/src/fix_langevin.txt +++ b/doc/src/fix_langevin.txt @@ -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 diff --git a/doc/src/fix_langevin_spin.txt b/doc/src/fix_langevin_spin.txt index e4065adad5..e5dccc5e57 100644 --- a/doc/src/fix_langevin_spin.txt +++ b/doc/src/fix_langevin_spin.txt @@ -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. diff --git a/doc/src/fix_nve_spin.txt b/doc/src/fix_nve_spin.txt index 7b382bb6ad..8e6284639b 100644 --- a/doc/src/fix_nve_spin.txt +++ b/doc/src/fix_nve_spin.txt @@ -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 diff --git a/doc/src/min_modify.txt b/doc/src/min_modify.txt index 434f7d05e6..c23ca67852 100644 --- a/doc/src/min_modify.txt +++ b/doc/src/min_modify.txt @@ -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. diff --git a/doc/src/min_spin.txt b/doc/src/min_spin.txt index 2a0f0e5397..f5bdf428bb 100644 --- a/doc/src/min_spin.txt +++ b/doc/src/min_spin.txt @@ -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). diff --git a/doc/src/min_style.txt b/doc/src/min_style.txt index 7c8aa0ae29..8816d62ffa 100644 --- a/doc/src/min_style.txt +++ b/doc/src/min_style.txt @@ -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 diff --git a/doc/src/minimize.txt b/doc/src/minimize.txt index b109235ecc..da203c7d2d 100644 --- a/doc/src/minimize.txt +++ b/doc/src/minimize.txt @@ -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 diff --git a/doc/src/neb_spin.txt b/doc/src/neb_spin.txt index e72ec63b06..217606289f 100644 --- a/doc/src/neb_spin.txt +++ b/doc/src/neb_spin.txt @@ -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:] diff --git a/doc/src/pair_spin_dipole.txt b/doc/src/pair_spin_dipole.txt index d903a3f07d..36bcd25d5e 100644 --- a/doc/src/pair_spin_dipole.txt +++ b/doc/src/pair_spin_dipole.txt @@ -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:] diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index fd4612bf78..13088139f9 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -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 diff --git a/examples/SPIN/bfo/in.spin.bfo b/examples/SPIN/bfo/in.spin.bfo index e3c88b0f06..2cd9200121 100644 --- a/examples/SPIN/bfo/in.spin.bfo +++ b/examples/SPIN/bfo/in.spin.bfo @@ -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 diff --git a/examples/SPIN/cobalt_fcc/in.spin.cobalt_fcc b/examples/SPIN/cobalt_fcc/in.spin.cobalt_fcc index ea98eeba94..9193faa798 100644 --- a/examples/SPIN/cobalt_fcc/in.spin.cobalt_fcc +++ b/examples/SPIN/cobalt_fcc/in.spin.cobalt_fcc @@ -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 diff --git a/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp b/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp index 3f34838553..b9ede5f09c 100644 --- a/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp +++ b/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp @@ -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 diff --git a/examples/SPIN/dipole_spin/in.spin.iron_dipole_cut b/examples/SPIN/dipole_spin/in.spin.iron_dipole_cut index a409fe0563..34f7fea0d3 100644 --- a/examples/SPIN/dipole_spin/in.spin.iron_dipole_cut +++ b/examples/SPIN/dipole_spin/in.spin.iron_dipole_cut @@ -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 diff --git a/examples/SPIN/dipole_spin/in.spin.iron_dipole_ewald b/examples/SPIN/dipole_spin/in.spin.iron_dipole_ewald index 58b44b55fe..f694bc5ddc 100644 --- a/examples/SPIN/dipole_spin/in.spin.iron_dipole_ewald +++ b/examples/SPIN/dipole_spin/in.spin.iron_dipole_ewald @@ -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 diff --git a/examples/SPIN/dipole_spin/in.spin.iron_dipole_pppm b/examples/SPIN/dipole_spin/in.spin.iron_dipole_pppm index 28d7e4a4bc..4175038ade 100644 --- a/examples/SPIN/dipole_spin/in.spin.iron_dipole_pppm +++ b/examples/SPIN/dipole_spin/in.spin.iron_dipole_pppm @@ -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 diff --git a/examples/SPIN/iron/in.spin.iron b/examples/SPIN/iron/in.spin.iron index bb1b0e1b4d..3468575493 100644 --- a/examples/SPIN/iron/in.spin.iron +++ b/examples/SPIN/iron/in.spin.iron @@ -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 diff --git a/examples/SPIN/iron/in.spin.iron_cubic b/examples/SPIN/iron/in.spin.iron_cubic index d4703a2959..859d9df0fa 100644 --- a/examples/SPIN/iron/in.spin.iron_cubic +++ b/examples/SPIN/iron/in.spin.iron_cubic @@ -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 diff --git a/examples/SPIN/nickel/in.spin.nickel b/examples/SPIN/nickel/in.spin.nickel index 0ed2fac410..caa1c940ae 100644 --- a/examples/SPIN/nickel/in.spin.nickel +++ b/examples/SPIN/nickel/in.spin.nickel @@ -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 diff --git a/examples/SPIN/nickel/in.spin.nickel_cubic b/examples/SPIN/nickel/in.spin.nickel_cubic index 3c97b284ae..76ea23689a 100644 --- a/examples/SPIN/nickel/in.spin.nickel_cubic +++ b/examples/SPIN/nickel/in.spin.nickel_cubic @@ -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 diff --git a/examples/SPIN/read_restart/in.spin.read_data b/examples/SPIN/read_restart/in.spin.read_data index a450421699..e788ecf67e 100644 --- a/examples/SPIN/read_restart/in.spin.read_data +++ b/examples/SPIN/read_restart/in.spin.read_data @@ -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 diff --git a/examples/SPIN/read_restart/in.spin.restart b/examples/SPIN/read_restart/in.spin.restart index 39157fdac4..ccce25b254 100644 --- a/examples/SPIN/read_restart/in.spin.restart +++ b/examples/SPIN/read_restart/in.spin.restart @@ -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 diff --git a/examples/SPIN/read_restart/in.spin.write_restart b/examples/SPIN/read_restart/in.spin.write_restart index 42f07fd316..c127101093 100644 --- a/examples/SPIN/read_restart/in.spin.write_restart +++ b/examples/SPIN/read_restart/in.spin.write_restart @@ -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 diff --git a/examples/SPIN/setforce_spin/in.spinmin.setforce b/examples/SPIN/setforce_spin/in.spinmin.setforce index 10d4df66ed..822768e0ef 100644 --- a/examples/SPIN/setforce_spin/in.spinmin.setforce +++ b/examples/SPIN/setforce_spin/in.spinmin.setforce @@ -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 diff --git a/examples/SPIN/spinmin/in.spinmin_cg.bfo b/examples/SPIN/spinmin/in.spinmin_cg.bfo new file mode 100644 index 0000000000..9d57399a56 --- /dev/null +++ b/examples/SPIN/spinmin/in.spinmin_cg.bfo @@ -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 diff --git a/examples/SPIN/spinmin/in.spinmin_lbfgs.bfo b/examples/SPIN/spinmin/in.spinmin_lbfgs.bfo new file mode 100644 index 0000000000..56cd6b8fae --- /dev/null +++ b/examples/SPIN/spinmin/in.spinmin_lbfgs.bfo @@ -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 diff --git a/examples/gjf/README.md b/examples/gjf/README.md new file mode 100644 index 0000000000..e6886cb2dd --- /dev/null +++ b/examples/gjf/README.md @@ -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 diff --git a/examples/gjf/argon.lmp b/examples/gjf/argon.lmp new file mode 100644 index 0000000000..00214b4c54 --- /dev/null +++ b/examples/gjf/argon.lmp @@ -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 + diff --git a/examples/gjf/ff-argon.lmp b/examples/gjf/ff-argon.lmp new file mode 100644 index 0000000000..b6f7bc931a --- /dev/null +++ b/examples/gjf/ff-argon.lmp @@ -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 + diff --git a/examples/gjf/in.gjf.vfull b/examples/gjf/in.gjf.vfull new file mode 100644 index 0000000000..19420e22ca --- /dev/null +++ b/examples/gjf/in.gjf.vfull @@ -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 + + diff --git a/examples/gjf/in.gjf.vhalf b/examples/gjf/in.gjf.vhalf new file mode 100644 index 0000000000..74e2089595 --- /dev/null +++ b/examples/gjf/in.gjf.vhalf @@ -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 + + diff --git a/src/.gitignore b/src/.gitignore index 2d79018b28..5848874d94 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -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 diff --git a/src/KOKKOS/fix_langevin_kokkos.cpp b/src/KOKKOS/fix_langevin_kokkos.cpp index 651f790a25..426bcb49e3 100644 --- a/src/KOKKOS/fix_langevin_kokkos.cpp +++ b/src/KOKKOS/fix_langevin_kokkos.cpp @@ -61,7 +61,6 @@ FixLangevinKokkos::FixLangevinKokkos(LAMMPS *lmp, int narg, char **a k_ratio.template modify(); if(gjfflag){ - nvalues = 3; grow_arrays(atomKK->nmax); atom->add_callback(0); // initialize franprev to zero @@ -69,8 +68,12 @@ FixLangevinKokkos::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(); + k_lv.template modify(); } if(zeroflag){ k_fsumall = tdual_double_1d_3n("langevin:fsumall"); @@ -94,6 +97,7 @@ FixLangevinKokkos::~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::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(); @@ -121,6 +129,40 @@ void FixLangevinKokkos::grow_arrays(int nmax) memoryKK->grow_kokkos(k_franprev,franprev,nmax,3,"langevin:franprev"); d_franprev = k_franprev.template view(); h_franprev = k_franprev.template view(); + memoryKK->grow_kokkos(k_lv,lv,nmax,3,"langevin:lv"); + d_lv = k_lv.template view(); + h_lv = k_lv.template view(); +} + +/* ---------------------------------------------------------------------- */ + +template +void FixLangevinKokkos::initial_integrate(int vflag) +{ + atomKK->sync(execution_space,datamask_read); + atomKK->modified(execution_space,datamask_modify); + + v = atomKK->k_v.view(); + f = atomKK->k_f.view(); + int nlocal = atomKK->nlocal; + if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst; + + FixLangevinKokkosInitialIntegrateFunctor functor(this); + Kokkos::parallel_for(nlocal,functor); +} + +template +KOKKOS_INLINE_FUNCTION +void FixLangevinKokkos::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::post_force(int vflag) k_gfactor2.template sync(); k_ratio.template sync(); if(gjfflag) k_franprev.template sync(); + if(gjfflag) k_lv.template sync(); boltz = force->boltz; dt = update->dt; @@ -177,7 +220,7 @@ void FixLangevinKokkos::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::post_force(int vflag) // set modify flags for the views modified in post_force functor if (gjfflag) k_franprev.template modify(); + if (gjfflag) k_lv.template modify(); if (tallyflag) k_flangevin.template modify(); // set total force to zero @@ -550,6 +594,10 @@ FSUM FixLangevinKokkos::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::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::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::compute_energy_item(int i) const template void FixLangevinKokkos::end_of_step() { - if (!tallyflag) return; + if (!tallyflag && !gjfflag) return; v = atomKK->k_v.template view(); + f = atomKK->k_f.template view(); mask = atomKK->k_mask.template view(); atomKK->sync(execution_space,V_MASK | MASK_MASK); @@ -733,9 +793,81 @@ void FixLangevinKokkos::end_of_step() FixLangevinKokkosTallyEnergyFunctor tally_functor(this); Kokkos::parallel_reduce(nlocal,tally_functor,energy_onestep); + if (gjfflag){ + if (rmass.data()) { + FixLangevinKokkosEndOfStepFunctor functor(this); + Kokkos::parallel_for(nlocal,functor); + } else { + mass = atomKK->k_mass.view(); + FixLangevinKokkosEndOfStepFunctor functor(this); + Kokkos::parallel_for(nlocal,functor); + } + } + energy += energy_onestep*update->dt; } +template +KOKKOS_INLINE_FUNCTION +void FixLangevinKokkos::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 +KOKKOS_INLINE_FUNCTION +void FixLangevinKokkos::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::end_of_step() template void FixLangevinKokkos::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(); + k_lv.template modify(); } @@ -765,6 +902,7 @@ void FixLangevinKokkos::cleanup_copy() tforce = NULL; gjfflag = 0; franprev = NULL; + lv = NULL; id = style = NULL; vatom = NULL; } diff --git a/src/KOKKOS/fix_langevin_kokkos.h b/src/KOKKOS/fix_langevin_kokkos.h index 140fea81d6..a6d467dfd7 100644 --- a/src/KOKKOS/fix_langevin_kokkos.h +++ b/src/KOKKOS/fix_langevin_kokkos.h @@ -56,6 +56,9 @@ namespace LAMMPS_NS { template class FixLangevinKokkos; + template + class FixLangevinKokkosInitialIntegrateFunctor; + template 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 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::t_float_1d rmass; + typename ArrayTypes::t_float_1d mass; typename ArrayTypes::tdual_double_2d k_franprev; typename ArrayTypes::t_double_2d d_franprev; HAT::t_double_2d h_franprev; + typename ArrayTypes::tdual_double_2d k_lv; + typename ArrayTypes::t_double_2d d_lv; + HAT::t_double_2d h_lv; + typename ArrayTypes::tdual_double_2d k_flangevin; typename ArrayTypes::t_double_2d d_flangevin; HAT::t_double_2d h_flangevin; @@ -130,6 +151,21 @@ namespace LAMMPS_NS { }; + template + struct FixLangevinKokkosInitialIntegrateFunctor { + typedef DeviceType device_type ; + FixLangevinKokkos c; + + FixLangevinKokkosInitialIntegrateFunctor(FixLangevinKokkos* c_ptr): + c(*c_ptr) {c.cleanup_copy();}; + + KOKKOS_INLINE_FUNCTION + void operator()(const int i) const { + c.initial_integrate_item(i); + } + }; + + template struct FixLangevinKokkosPostForceFunctor { @@ -207,6 +243,21 @@ namespace LAMMPS_NS { update += source; } }; + + template + struct FixLangevinKokkosEndOfStepFunctor { + typedef DeviceType device_type ; + FixLangevinKokkos c; + + FixLangevinKokkosEndOfStepFunctor(FixLangevinKokkos* 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. + */ diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp index b1b466b5a4..9b4f1916ae 100644 --- a/src/SPIN/fix_nve_spin.cpp +++ b/src/SPIN/fix_nve_spin.cpp @@ -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"); diff --git a/src/SPIN/fix_nve_spin.h b/src/SPIN/fix_nve_spin.h index 1e3b87c1fb..5871f721be 100644 --- a/src/SPIN/fix_nve_spin.h +++ b/src/SPIN/fix_nve_spin.h @@ -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 diff --git a/src/SPIN/min_spin.cpp b/src/SPIN/min_spin.cpp index 553526ea8e..7315aca056 100644 --- a/src/SPIN/min_spin.cpp +++ b/src/SPIN/min_spin.cpp @@ -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; -} - diff --git a/src/SPIN/min_spin.h b/src/SPIN/min_spin.h index fbc624a9cc..f2df81e58c 100644 --- a/src/SPIN/min_spin.h +++ b/src/SPIN/min_spin.h @@ -35,7 +35,6 @@ class MinSpin : public Min { int iterate(int); double evaluate_dt(); void advance_spins(double); - double fmnorm_sqr(); private: diff --git a/src/SPIN/min_spin_cg.cpp b/src/SPIN/min_spin_cg.cpp new file mode 100644 index 0000000000..95bbcf437b --- /dev/null +++ b/src/SPIN/min_spin_cg.cpp @@ -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 +#include +#include +#include +#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; +} diff --git a/src/SPIN/min_spin_cg.h b/src/SPIN/min_spin_cg.h new file mode 100644 index 0000000000..0eed7a61e6 --- /dev/null +++ b/src/SPIN/min_spin_cg.h @@ -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 diff --git a/src/SPIN/min_spin_lbfgs.cpp b/src/SPIN/min_spin_lbfgs.cpp new file mode 100644 index 0000000000..f86bdd5d48 --- /dev/null +++ b/src/SPIN/min_spin_lbfgs.cpp @@ -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 +#include +#include +#include +#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; +} diff --git a/src/SPIN/min_spin_lbfgs.h b/src/SPIN/min_spin_lbfgs.h new file mode 100644 index 0000000000..cead605b32 --- /dev/null +++ b/src/SPIN/min_spin_lbfgs.h @@ -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 diff --git a/src/SPIN/neb_spin.cpp b/src/SPIN/neb_spin.cpp index 2bcfe6573a..559fd1cb49 100644 --- a/src/SPIN/neb_spin.cpp +++ b/src/SPIN/neb_spin.cpp @@ -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; diff --git a/src/SPIN/pair_spin_dipole_cut.cpp b/src/SPIN/pair_spin_dipole_cut.cpp index fb0ec8bb20..ebcc59d499 100644 --- a/src/SPIN/pair_spin_dipole_cut.cpp +++ b/src/SPIN/pair_spin_dipole_cut.cpp @@ -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]; diff --git a/src/SPIN/pair_spin_dipole_long.cpp b/src/SPIN/pair_spin_dipole_long.cpp index 670fccfce2..bb42fb70b4 100644 --- a/src/SPIN/pair_spin_dipole_long.cpp +++ b/src/SPIN/pair_spin_dipole_long.cpp @@ -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]; diff --git a/src/SPIN/pair_spin_dmi.cpp b/src/SPIN/pair_spin_dmi.cpp index afd2deaa6a..3fee84d5fc 100644 --- a/src/SPIN/pair_spin_dmi.cpp +++ b/src/SPIN/pair_spin_dmi.cpp @@ -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]); } /* ---------------------------------------------------------------------- diff --git a/src/SPIN/pair_spin_exchange.cpp b/src/SPIN/pair_spin_exchange.cpp index cc28018ad0..7a54eba9d7 100644 --- a/src/SPIN/pair_spin_exchange.cpp +++ b/src/SPIN/pair_spin_exchange.cpp @@ -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); diff --git a/src/SPIN/pair_spin_neel.cpp b/src/SPIN/pair_spin_neel.cpp index 355ba20f39..017682593a 100644 --- a/src/SPIN/pair_spin_neel.cpp +++ b/src/SPIN/pair_spin_neel.cpp @@ -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); diff --git a/src/USER-MISC/fix_bond_react.cpp b/src/USER-MISC/fix_bond_react.cpp index 6b61a7b84d..5bf6ed68b9 100644 --- a/src/USER-MISC/fix_bond_react.cpp +++ b/src/USER-MISC/fix_bond_react.cpp @@ -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; } diff --git a/src/USER-MISC/fix_bond_react.h b/src/USER-MISC/fix_bond_react.h index e5452cb226..eda26f129d 100644 --- a/src/USER-MISC/fix_bond_react.h +++ b/src/USER-MISC/fix_bond_react.h @@ -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; diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp index 1e86a90218..f805fddb4b 100644 --- a/src/fix_langevin.cpp +++ b/src/fix_langevin.cpp @@ -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; } diff --git a/src/fix_langevin.h b/src/fix_langevin.h index 4b5570ac2e..868b71a44d 100644 --- a/src/fix_langevin.h +++ b/src/fix_langevin.h @@ -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 diff --git a/src/min.cpp b/src/min.cpp index b6bed01c75..b57dd44d4f 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -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 ------------------------------------------------------------------------- */ diff --git a/src/min.h b/src/min.h index 835e6b3ed5..61f9ce0bda 100644 --- a/src/min.h +++ b/src/min.h @@ -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); diff --git a/src/min_cg.cpp b/src/min_cg.cpp index ff318d23da..80dde25f51 100644 --- a/src/min_cg.cpp +++ b/src/min_cg.cpp @@ -14,6 +14,7 @@ #include "min_cg.h" #include #include +#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 diff --git a/src/min_fire.cpp b/src/min_fire.cpp index b7f853afd2..e0cc2437d4 100644 --- a/src/min_fire.cpp +++ b/src/min_fire.cpp @@ -16,6 +16,7 @@ #include #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 { diff --git a/src/min_hftn.cpp b/src/min_hftn.cpp index dcfed10fef..63432aab63 100644 --- a/src/min_hftn.cpp +++ b/src/min_hftn.cpp @@ -22,6 +22,7 @@ #include #include #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]; diff --git a/src/min_quickmin.cpp b/src/min_quickmin.cpp index 04e72b046d..ef649b4dac 100644 --- a/src/min_quickmin.cpp +++ b/src/min_quickmin.cpp @@ -16,6 +16,7 @@ #include #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 { diff --git a/src/min_sd.cpp b/src/min_sd.cpp index e4e159003f..8541b0ccdf 100644 --- a/src/min_sd.cpp +++ b/src/min_sd.cpp @@ -13,6 +13,7 @@ #include "min_sd.h" #include +#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) diff --git a/src/neighbor.cpp b/src/neighbor.cpp index d38aed08c0..2c77a13258 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -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;