lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into sdpd to resolve merge conflicts 

# Conflicts:
#	cmake/CMakeLists.txt
#	src/Makefile
This commit is contained in:
Axel Kohlmeyer 2018-11-08 16:32:02 -05:00
parent bbb371134d 3681bc853d
commit e0955f6434
71 changed files with 9459 additions and 528 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
cmake/CMakeLists.txt
View File

@ -171,8 +171,9 @@ set(DEFAULT_PACKAGES ASPHERE BODY CLASS2 COLLOID COMPRESS DIPOLE GRANULAR
USER-BOCS USER-CGDNA USER-MESO USER-CGSDK USER-COLVARS USER-DIFFRACTION
USER-DPD USER-DRUDE USER-EFF USER-FEP USER-H5MD USER-LB USER-MANIFOLD
USER-MEAMC USER-MGPT USER-MISC USER-MOFFF USER-MOLFILE USER-NETCDF
USER-PHONON USER-PTM USER-QTB USER-REAXC USER-SCAFACOS USER-SDPD USER-SMD
USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK USER-QUIP USER-QMMM)
USER-PHONON USER-PLUMED USER-PTM USER-QTB USER-REAXC USER-SCAFACOS
USER-SDPD USER-SMD USER-SMTBQ USER-SPH USER-TALLY USER-UEF USER-VTK
USER-QUIP USER-QMMM)
set(ACCEL_PACKAGES USER-OMP KOKKOS OPT USER-INTEL GPU)
set(OTHER_PACKAGES CORESHELL QEQ)
foreach(PKG ${DEFAULT_PACKAGES})
@ -193,6 +194,8 @@ if(PKG_MEAM OR PKG_USER-H5MD OR PKG_USER-QMMM OR PKG_USER-SCAFACOS)
enable_language(C)
endif()
include_directories(${LAMMPS_SOURCE_DIR})
# do MPI detection after language activation, if MPI for these language is required
find_package(MPI QUIET)
option(BUILD_MPI "Build MPI version" ${MPI_FOUND})
@ -526,6 +529,32 @@ if(PKG_USER-SCAFACOS)
include_directories(${SCAFACOS_INCLUDE_DIRS})
endif()
if(PKG_USER-PLUMED)
find_package(GSL REQUIRED)
option(DOWNLOAD_PLUMED "Download Plumed (instead of using the system's one)" OFF)
if(DOWNLOAD_PLUMED)
include(ExternalProject)
ExternalProject_Add(plumed_build
URL https://github.com/plumed/plumed2/releases/download/v2.4.3/plumed-src-2.4.3.tgz
URL_MD5 b1be7c48971627febc11c61b70767fc5
BUILD_IN_SOURCE 1
CONFIGURE_COMMAND <SOURCE_DIR>/configure --prefix=<INSTALL_DIR>
$<$<BOOL:${BUILD_SHARED_LIBS}>:--with-pic> )
ExternalProject_get_property(plumed_build INSTALL_DIR)
set(PLUMED_INSTALL_DIR ${INSTALL_DIR})
list(APPEND LAMMPS_DEPS plumed_build)
list(APPEND LAMMPS_LINK_LIBS ${PLUMED_INSTALL_DIR}/lib/plumed/obj/kernel.o
${PLUMED_INSTALL_DIR}/lib/plumed/obj/PlumedStatic.o ${GSL_LIBRARIES} ${CMAKE_DL_LIBS})
set(PLUMED_INCLUDE_DIRS "${PLUMED_INSTALL_DIR}/include")
else()
find_package(PkgConfig REQUIRED)
pkg_check_modules(PLUMED plumed REQUIRED)
include(${PLUMED_LIBDIR}/plumed/src/lib/Plumed.cmake.static)
list(APPEND LAMMPS_LINK_LIBS ${PLUMED_LOAD})
endif()
include_directories(${PLUMED_INCLUDE_DIRS})
endif()
if(PKG_USER-MOLFILE)
add_library(molfile INTERFACE)
target_include_directories(molfile INTERFACE ${LAMMPS_LIB_SOURCE_DIR}/molfile)
@ -1180,7 +1209,6 @@ set(LAMMPS_STYLE_HEADERS_DIR ${CMAKE_CURRENT_BINARY_DIR}/styles)
GenerateStyleHeaders(${LAMMPS_STYLE_HEADERS_DIR})
include_directories(${LAMMPS_SOURCE_DIR})
include_directories(${LAMMPS_STYLE_HEADERS_DIR})
######################################

24
cmake/pkgconfig/liblammps.pc.in
View File

@ -1,9 +1,29 @@
# pkg-config file for lammps
# https://people.freedesktop.org/~dbn/pkg-config-guide.html
# Usage: cc `pkg-config --cflags --libs liblammps` -o myapp myapp.c
# after you added @CMAKE_INSTALL_FULL_LIBDIR@/pkg-config to PKG_CONFIG_PATH,
# Add the directory where lammps.pc got installed to your PKG_CONFIG_PATH
# e.g. export PKG_CONFIG_PATH=@CMAKE_INSTALL_FULL_LIBDIR@/pkgconfig
# Use this on commandline with:
# c++ `pkg-config --cflags --libs lammps` -o myapp myapp.cpp
# Use this in a Makefile:
# myapp: myapp.cpp
# $(CC) `pkg-config --cflags --libs lammps` -o $@ $<
# Use this in autotools:
# configure.ac:
# PKG_CHECK_MODULES([LAMMPS], [lammps])
# Makefile.am:
# myapp_CFLAGS = $(LAMMPS_CFLAGS)
# myapp_LDADD = $(LAMMPS_LIBS)
# Use this in CMake:
# CMakeLists.txt:
# find_package(PkgConfig)
# pkg_check_modules(LAMMPS IMPORTED_TARGET lammps)
# target_link_libraries(<lib> PkgConfig::LAMMPS)
prefix=@CMAKE_INSTALL_PREFIX@
libdir=@CMAKE_INSTALL_FULL_LIBDIR@
includedir=@CMAKE_INSTALL_FULL_INCLUDEDIR@

57
doc/src/Build_extras.txt
View File

@ -41,6 +41,7 @@ This is the list of packages that may require additional steps.
"USER-ATC"_#user-atc,
"USER-AWPMD"_#user-awpmd,
"USER-COLVARS"_#user-colvars,
"USER-PLUMED" _#user-plumed,
"USER-H5MD"_#user-h5md,
"USER-INTEL"_#user-intel,
"USER-MOLFILE"_#user-molfile,
@ -712,6 +713,62 @@ a corresponding Makefile.lammps.machine file.
:line
USER-PLUMED package :h4,link(user-plumed)
[CMake build]:
[Traditional make]:
Before building LAMMPS with this package, you must first build
PLUMED. We recommending building PLUMED separately to LAMMPS using
the instructions that can be found at http://plumed.github.io/doc-master/user-doc/html/_installation.html.
Before compiling LAMMPS you can then install the fix plumed command
and compile LAMMPS in the usual manner:
make yes-user-plumed
make machine :pre
Once this compilation completes you should be able to run LAMMPS in the usual
way. When running LAMMPS with an input script that contains a fix
plumed command LAMMPS will try to call the PLUMED runtime library. PLUMED
must therefore be available in your path if LAMMPS is compiled in this way.
On some machines it is not possible to call runtime libraries in the way described
above. When compiling on these machines it is thus better to statically link
PLUMED when compiling LAMMPS. To do this you must either download a PLUMED
tarball from http://www.plumed.org/get-it or clone it using
git clone https://github.com/plumed/plumed2.git. If you download the tarball
unpack it in the /lib/plumed directory. Similarly if you clone
it clone it to the /lib/plumed directory as if there is a version of PLUMED within
this directory LAMMPS will always try to statically link the version of PLUMED
that this directory contains instead of dynamically linking the library.
Once you have downloaded PLUMED into /lib/plumed you must again build the code
here by following the instructions that can be found at
http://plumed.github.io/doc-master/user-doc/html/_installation.html.
You can statically link PLUMED manually and if you want to access the full
range of PLUMED functionalities this is what you should do. If you only want the
basic range of functionalities, however, (i.e. no user contributed modules) then
you can download and compile PLUMED in one step from the lammps/src dir, using a
command like like those below:
make lib-plumed # print help message
make lib-plumed args="-b" # download and build the latest stable version of PLUMED
These commands will simply invoke the lib/plumed/Install.py script with
args specified. Furthermore, once the script has completed you should
have a compiled version of PLUMED. With this built you can install/un-install
PLUMED and build LAMMPS in the usual manner:
make yes-user-plumed
make machine :pre
make no-user-plumed
make machine :pre
:line
USER-H5MD package :h4,link(user-h5md)
To build with this package you must have the HDF5 software package

BIN
doc/src/PDF/colvars-refman-lammps.pdf
View File

Binary file not shown.

28
doc/src/Packages_details.txt
View File

@ -1201,6 +1201,34 @@ examples/USER/colvars :ul
:line
USER-PLUMED package :link(USER-PLUMED),h4
[Contents:]
The fix plumed command allows you to use the plugin for molecular
dynamics PLUMED to analyse and bias your LAMMPS trajectory on the fly.
In practise PLUMED is called from within the lammps input script by using
the "fix plumed _fix_plumed.html command.
[Authors:] The PLUMED library is written and maintained by
Massimilliano Bonomi, Giovanni Bussi, Carlo Camiloni and
Gareth Tribello.
[Install:]
This package has "specific installation
instructions"_Build_extras.html#gpu on the "Build
extras"_Build_extras.html doc page.
[Supporting info:]
src/USER-PLUMED/README
lib/plumed/README
"fix plumed "_fix_plumed.html
examples/USER/plumed :ul
:line
USER-DIFFRACTION package :link(PKG-USER-DIFFRACTION),h4
[Contents:]

4
doc/src/fix_bond_break.txt
View File

@ -137,8 +137,8 @@ doc page for more info.
[Related commands:]
"fix bond/create"_fix_bond_create.html, "fix
bond/swap"_fix_bond_swap.html, "dump local"_dump.html,
"special_bonds"_special_bonds.html
bond/react"_fix_bond_react.html, "fix bond/swap"_fix_bond_swap.html,
"dump local"_dump.html, "special_bonds"_special_bonds.html
[Default:]

4
doc/src/fix_bond_create.txt
View File

@ -232,8 +232,8 @@ doc page for more info.
[Related commands:]
"fix bond/break"_fix_bond_break.html, "fix
bond/swap"_fix_bond_swap.html, "dump local"_dump.html,
"special_bonds"_special_bonds.html
bond/react"_fix_bond_react.html, "fix bond/swap"_fix_bond_swap.html,
"dump local"_dump.html, "special_bonds"_special_bonds.html
[Default:]

95
doc/src/fix_bond_react.txt
View File

@ -24,11 +24,11 @@ common_keyword = {stabilization} :l
{stabilization} values = {no} or {yes} {group-ID} {xmax}
{no} = no reaction site stabilization
{yes} = perform reaction site stabilization
{group-ID} = user-assigned ID for all non-reacting atoms (group created internally)
{group-ID} = user-assigned prefix for the dynamic group of non-reacting atoms
{xmax} = xmax value that is used by an internally created "nve/limit"_fix_nve_limit.html integrator :pre
react = mandatory argument indicating new reaction specification :l
react-ID = user-assigned name for the reaction :l
react-group-ID = only atoms in this group are available for the reaction :l
react-group-ID = only atoms in this group are considered for the reaction :l
Nevery = attempt reaction every this many steps :l
Rmin = bonding pair atoms must be separated by more than Rmin to initiate reaction (distance units) :l
Rmax = bonding pair atoms must be separated by less than Rmax to initiate reaction (distance units) :l
@ -41,14 +41,18 @@ react = mandatory argument indicating new reaction specification :l
fraction = initiate reaction with this probability if otherwise eligible
seed = random number seed (positive integer)
{stabilize_steps} value = timesteps
timesteps = number of timesteps to apply internally created nve/limit.html :pre
timesteps = number of timesteps to apply internally created nve/limit.html
{update_edges} value = {none} or {charges} :l
none = do not update topology near the edges of reaction templates
charges = update atomic charges of all atoms in reaction templates
custom = force the update of user-specified atomic charges :pre
:ule
[Examples:]
molecule mol1 pre_reacted_topology.txt
molecule mol2 post_reacted_topology.txt
fix 5 all bond/react stabilization no react myrxn1 all 1 0 3.25 mol1 mol2 map_file.txt :pre
fix 5 all bond/react react myrxn1 all 1 0 3.25 mol1 mol2 map_file.txt :pre
molecule mol1 pre_reacted_rxn1.txt
molecule mol2 post_reacted_rxn1.txt
@ -57,7 +61,7 @@ molecule mol4 post_reacted_rxn2.txt
fix 5 all bond/react stabilization yes nvt_grp .03 &
react myrxn1 all 1 0 3.25 mol1 mol2 map_file_rxn1.txt prob 0.50 12345 &
react myrxn2 all 1 0 2.75 mol3 mol4 map_file_rxn2.txt prob 0.25 12345
fix 6 nvt_grp nvt temp 300 300 100 # set thermostat after bond/react :pre
fix 6 nvt_grp_REACT nvt temp 300 300 100 # set thermostat after bond/react :pre
[Description:]
@ -99,19 +103,29 @@ involved in any new reactions. The {xmax} value keyword should
typically be set to the maximum distance that non-reacting atoms move
during the simulation.
The group-ID set using the {stabilization} keyword should be a
previously unused group-ID. It cannot be specified as 'all'. The fix
bond/react command creates a "dynamic group"_group.html of this name
that includes all non-reacting atoms. This dynamic group-ID should
then be used by a subsequent system-wide time integrator such as nvt,
npt, or nve, as shown in the second example above. It is currently
necessary to place the time integration command after the fix
bond/react command due to the internal dynamic grouping performed by
fix bond/react.
The group-ID set using the {stabilization} keyword can be an existing
static group or a previously-unused group-ID. It cannot be specified
as 'all'. If the group-ID is previously unused, the fix bond/react
command creates a "dynamic group"_group.html that is initialized to
include all atoms. If the group-ID is that of an existing static
group, the group is used as the parent group of new,
internally-created dynamic group. In both cases, this new dynamic
group is named by appending '_REACT' to the group-ID, e.g.
nvt_grp_REACT. By specifying an existing group, you may thermostat
constant-topology parts of your system separately. The dynamic group
contains only non-reacting atoms at a given timestep, and therefore
should be used by a subsequent system-wide time integrator such as
nvt, npt, or nve, as shown in the second example above. The time
integration command should be placed after the fix bond/react command
due to the internal dynamic grouping performed by fix bond/react.
NOTE: The internally created group currently applies to all atoms in
the system, i.e. you should generally not have a separate thermostat
which acts on the 'all' group.
NOTE: If the group-ID is an existing static group, react-group-IDs
should also be specified as this static group, or a subset.
NOTE: If the group-ID is previously unused, the internally created
group applies to all atoms in the system, i.e. you should generally
not have a separate thermostat which acts on the 'all' group, or any
other group.
The following comments pertain to each {react} argument:
@ -155,7 +169,17 @@ Some atoms in the pre-reacted template that are not reacting may have
missing topology with respect to the simulation. For example, the
pre-reacted template may contain an atom that would connect to the
rest of a long polymer chain. These are referred to as edge atoms, and
are also specified in the map file.
are also specified in the map file. When the pre-reaction template
contains edge atoms, not all atoms, bonds, charges, etc. specified in
the reaction templates will be updated. Specifically, topology that
involves only atoms that are 'too near' to template edges will not be
updated. The definition of 'too near the edge' depends on which
interactions are defined in the simulation. If the simulation has
defined dihedrals, atoms within two bonds of edge atoms are considered
'too near the edge.' If the simulation defines angles, but not
dihedrals, atoms within one bond of edge atoms are considered 'too
near the edge.' If just bonds are defined, only edge atoms are
considered 'too near the edge.'
Note that some care must be taken when a building a molecule template
for a given simulation. All atom types in the pre-reacted template
@ -178,23 +202,30 @@ A discussion of correctly handling this is also provided on the
The map file is a text document with the following format:
A map file has a header and a body. The header of map file the
contains one mandatory keyword and one optional keyword. The mandatory
keyword is 'equivalences' and the optional keyword is 'edgeIDs':
contains one mandatory keyword and two optional keywords. The mandatory
keyword is 'equivalences' and the optional keywords are 'edgeIDs' and
'customIDs':
N {equivalences} = # of atoms N in the reaction molecule templates
N {edgeIDs} = # of edge atoms N in the pre-reacted molecule template :pre
N {edgeIDs} = # of edge atoms N in the pre-reacted molecule template
N {customIDs} = # of atoms N that are specified for a custom update :pre
The body of the map file contains two mandatory sections and one
optional section. The first mandatory section begins with the keyword
The body of the map file contains two mandatory sections and two
optional sections. The first mandatory section begins with the keyword
'BondingIDs' and lists the atom IDs of the bonding atom pair in the
pre-reacted molecule template. The second mandatory section begins
with the keyword 'Equivalences' and lists a one-to-one correspondence
between atom IDs of the pre- and post-reacted templates. The first
column is an atom ID of the pre-reacted molecule template, and the
second column is the corresponding atom ID of the post-reacted
molecule template. The optional section begins with the keyword
molecule template. The first optional section begins with the keyword
'EdgeIDs' and lists the atom IDs of edge atoms in the pre-reacted
molecule template.
molecule template. The second optional section begins with the keyword
'Custom Edges' and allows for forcing the update of a specific atom's
atomic charge. The first column is the ID of an atom near the edge of
the pre-reacted molecule template, and the value of the second column
is either 'none' or 'charges.' Further details are provided in the
discussion of the 'update_edges' keyword.
A sample map file is given below:
@ -255,6 +286,18 @@ The {stabilize_steps} keyword allows for the specification of how many
timesteps a reaction site is stabilized before being returned to the
overall system thermostat.
The {update_edges} keyword can increase the number of atoms whose
atomic charges are updated, when the pre-reaction template contains
edge atoms. When the value is set to 'charges,' all atoms' atomic
charges are updated to those specified by the post-reaction template,
including atoms near the edge of reaction templates. When the value is
set to 'custom,' an additional section must be included in the map
file that specifies whether to update charges, on a per-atom basis.
The format of this section is detailed above. Listing a pre-reaction
atom ID with a value of 'charges' will force the update of the atom's
charge, even if it is near a template edge. Atoms not near a template
edge are unaffected by this setting.
In order to produce the most physical behavior, this 'reaction site
equilibration time' should be tuned to be as small as possible while
retaining stability for a given system or reaction step. After a
@ -323,7 +366,7 @@ bond/break"_fix_bond_break.html, "fix bond/swap"_fix_bond_swap.html,
[Default:]
The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60
The option defaults are stabilization = no, prob = 1.0, stabilize_steps = 60, update_edges = none
:line

6
doc/src/fix_nve_dot.txt
View File

@ -36,8 +36,8 @@ The command is equivalent to the "fix nve"_fix_nve.html.
The particles are always considered to have a finite size.
An example input file can be found in /examples/USER/cgdna/examples/duplex1/.
A technical report with more information on this integrator can be found
"here"_PDF/USER-CGDNA-overview.pdf.
Further details of the implementation and stability of the integrator are contained in "(Henrich)"_#Henrich3.
The preprint version of the article can be found "here"_PDF/USER-CGDNA.pdf.
:line
@ -59,3 +59,5 @@ See the "Build package"_Build_package.html doc page for more info.
[(Davidchack)] R.L Davidchack, T.E. Ouldridge, and M.V. Tretyakov. J. Chem. Phys. 142, 144114 (2015).
:link(Miller1)
[(Miller)] T. F. Miller III, M. Eleftheriou, P. Pattnaik, A. Ndirango, G. J. Martyna, J. Chem. Phys., 116, 8649-8659 (2002).
:link(Henrich3)
[(Henrich)] O. Henrich, Y. A. Gutierrez-Fosado, T. Curk, T. E. Ouldridge, Eur. Phys. J. E 41, 57 (2018).

6
doc/src/fix_nve_dotc_langevin.txt
View File

@ -114,8 +114,8 @@ The scale factor after the {angmom} keyword gives the ratio of the rotational to
the translational friction coefficient.
An example input file can be found in /examples/USER/cgdna/examples/duplex2/.
A technical report with more information on this integrator can be found
"here"_PDF/USER-CGDNA-overview.pdf.
Further details of the implementation and stability of the integrators are contained in "(Henrich)"_#Henrich4.
The preprint version of the article can be found "here"_PDF/USER-CGDNA.pdf.
:line
@ -139,3 +139,5 @@ See the "Build package"_Build_package.html doc page for more info.
[(Miller)] T. F. Miller III, M. Eleftheriou, P. Pattnaik, A. Ndirango, G. J. Martyna, J. Chem. Phys., 116, 8649-8659 (2002).
:link(Dunweg3)
[(Dunweg)] B. Dunweg, W. Paul, Int. J. Mod. Phys. C, 2, 817-27 (1991).
:link(Henrich4)
[(Henrich)] O. Henrich, Y. A. Gutierrez-Fosado, T. Curk, T. E. Ouldridge, Eur. Phys. J. E 41, 57 (2018).

117
doc/src/fix_plumed.txt Normal file
View File

@ -0,0 +1,117 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Commands_all.html)
:line
fix plumed command :h3
[Syntax:]
fix ID group-ID plumed keyword value ... :pre
ID, group-ID are documented in "fix"_fix.html command :ulb,l
plumed = style name of this fix command :l
keyword = {plumedfile} or {outfile} :l
{plumedfile} arg = name of PLUMED input file to use (default: NULL)
{outfile} arg = name of file on which to write the PLUMED log (default: NULL) :pre
:ule
[Examples:]
fix pl all plumed all plumed plumedfile plumed.dat outfile p.log
[Description:]
This fix instructs LAMMPS to call the PLUMED library, which allows one
to perform various forms of trajectory analysis on the fly and to also
use methods such as umbrella sampling and metadynamics to enhance the
sampling of phase space.
The documentation included here only describes the fix plumed command.
This command is LAMMPS specific whereas most of the functionality
implemented in PLUMED will work with a range of MD codes and also when
PLUMED is used as a stand alone code. The full documentation for PLUMED
is available at "this website"_http://www.plumed.org/documentation
The PLUMED library is developed at
"https://github.com/plumed/plumed2"_https://github.com/plumed/plumed2 A
detailed discussion of the code can be found in "(PLUMED)"_#PLUMED.
There are some example scripts for using this package with LAMMPS in the
examples/USER/plumed directory.
:line
The command to call PLUMED above is reasonably self explanatory. Within
the input file for lammps the user is required to specify the input file
for PLUMED and a file on which to output the PLUMED log. The user must
specify both of these arguments every time PLUMED is to be used.
Furthermore, the fix plumed command should appear in the LAMMPS input
file after the relevant input paramters (e.g. the timestep) have been
set.
The {group-ID} entry is ignored. LAMMPS will always pass all the atoms
to PLUMED and there can only be one instance of the plumed fix at a
time. The plumed fix communicates the minimum amount of information
required and the PLUMED supports multiple, completely independent
collective variables, multiple independent biases and multiple
independent forms of analysis. There is thus really no restriction in
functionality by only allowing only one plumed fix in the LAMMPS input.
The {plumedfile} keyword allows the user to specify the name of the
PLUMED input file. Instructions as to what should be included in a
plumed input file can be found in the "documentation for
PLUMED"_http://www.plumed.org/documentation.
The {outfile} keyword allows the user to specify the name of a file on
which to output the PLUMED log. This log file normally just parots the
information that is contained in the input file. The names of the files
on which the results from the various analyses that have been performed
using PLUMED will be specified by the user in the PLUMED input file.
[Restart, fix_modify, output, run start/stop, minimize info:]
When performing a restart of a calculation that involves PLUMED you must
include a RESTART command in the PLUMED input file as detailed in the
"PLUMED documentation"_http://www.plumed.org/documentation. When the
restart command is found in the PLUMED input PLUMED will append to the
files that were generated in the run that was performed previously.
Furthermore, any history dependent bias potentials that were accumulated
in previous calculations will be read in when the restart command is
included in the PLUMED input.
The "fix_modify"_fix_modify.html {energy} option is not supported by
this fix.
Nothing is computed by this fix that can be accessed by any of the
"output commands"_Howto_output.html within LAMMPS. All the quantities
of interest can be output by commands that are native to PLUMED,
however.
[Restrictions:]
This fix is part of the USER-PLUMED package. It is only enabled if
LAMMPS was built with that package. See the "Build
package"_Build_package.html doc page for more info.
There can only be one plumed fix active at a time. Since the interface
communicates only the minimum amount of information and since the PLUMED
module itself can handle an arbitrary number of analysis and biasing
methods, this is not a limitation of functionality.
[Related commands:]
"fix smd"_fix_smd.html
"fix colvars"_fix_colvars.html
[Default:]
The default options are plumedfile = NULL and outfile = NULL
:line
:link(PLUMED)
[(PLUMED)] G.A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comp. Phys. Comm 185, 604 (2014)

1
doc/src/fixes.txt
View File

@ -118,6 +118,7 @@ Fixes :h1
fix_phonon
fix_pimd
fix_planeforce
fix_plumed
fix_poems
fix_pour
fix_precession_spin

3
doc/src/group.txt
View File

@ -225,8 +225,7 @@ atomfile-style variable. The variable is evaluated and atoms whose
per-atom values are 0.0, are removed from the dynamic group. If the {property}
keyword is used, the per-atom property name must be a previously defined
per-atom property. The per-atom property is evaluated and atoms whose
values are 0.0 are removed from the dynamic group, otherwise they
are added to the group.
values are 0.0 are removed from the dynamic group.
The assignment of atoms to a dynamic group is done at the beginning of
each run and on every timestep that is a multiple of {N}, which is the

1
doc/src/lammps.book
View File

@ -338,6 +338,7 @@ fix_orient.html
fix_phonon.html
fix_pimd.html
fix_planeforce.html
fix_plumed.html
fix_poems.html
fix_pour.html
fix_precession_spin.html

2
examples/USER/misc/bond_react/nylon,6-6_melt/in.large_nylon_melt
View File

@ -36,7 +36,7 @@ fix myrxns all bond/react stabilization yes statted_grp .03 &
react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map
# stable at 800K
fix 1 statted_grp nvt temp 800 800 100
fix 1 statted_grp_REACT nvt temp 800 800 100
# in order to customize behavior of reacting atoms,
# you can use the internally created 'bond_react_MASTER_group', like so:

2
examples/USER/misc/bond_react/tiny_nylon/in.tiny_nylon
View File

@ -36,7 +36,7 @@ fix myrxns all bond/react stabilization yes statted_grp .03 &
react rxn1 all 1 0.0 2.9 mol1 mol2 rxn1_stp1_map &
react rxn2 all 1 0.0 5.0 mol3 mol4 rxn1_stp2_map
fix 1 statted_grp nvt temp 300 300 100
fix 1 statted_grp_REACT nvt temp 300 300 100
fix 4 bond_react_MASTER_group temp/rescale 1 300 300 10 1

6531
examples/USER/plumed/data.peptide Normal file
View File

File diff suppressed because it is too large Load Diff

41
examples/USER/plumed/in.peptide-plumed Normal file
View File

@ -0,0 +1,41 @@
# Solvated 5-mer peptide
units real
atom_style full
pair_style lj/charmm/coul/long 8.0 10.0 10.0
bond_style harmonic
angle_style charmm
dihedral_style charmm
improper_style harmonic
kspace_style pppm 0.0001
read_data data.peptide
neighbor 2.0 bin
neigh_modify delay 5
timestep 2.0
group peptide type <= 12
group one id 2 4 5 6
group two id 80 82 83 84
group ref id 37
group colvar union one two ref
fix 1 all nvt temp 275.0 275.0 100.0 tchain 1
fix 2 all plumed plumedfile plumed.dat outfile p.log
fix 2a ref setforce 0.0 0.0 0.0
fix 4 all shake 0.0001 10 100 b 4 6 8 10 12 14 18 a 31
#dump 1 colvar custom 1 dump.colvar.lammpstrj id xu yu zu fx fy fz
#dump_modify 1 sort id
thermo_style custom step temp etotal pe ke epair ebond f_2
thermo 10
variable step equal step
variable pe equal pe
run 101

162
examples/USER/plumed/log.24Oct18.peptide-plumed.g++.1 Normal file
View File

@ -0,0 +1,162 @@
LAMMPS (24 Oct 2018)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:87)
using 1 OpenMP thread(s) per MPI task
# Solvated 5-mer peptide
units real
atom_style full
pair_style lj/charmm/coul/long 8.0 10.0 10.0
bond_style harmonic
angle_style charmm
dihedral_style charmm
improper_style harmonic
kspace_style pppm 0.0001
read_data data.peptide
orthogonal box = (36.8402 41.0137 29.7681) to (64.2116 68.3851 57.1395)
1 by 1 by 1 MPI processor grid
reading atoms ...
2004 atoms
reading velocities ...
2004 velocities
scanning bonds ...
3 = max bonds/atom
scanning angles ...
6 = max angles/atom
scanning dihedrals ...
14 = max dihedrals/atom
scanning impropers ...
1 = max impropers/atom
reading bonds ...
1365 bonds
reading angles ...
786 angles
reading dihedrals ...
207 dihedrals
reading impropers ...
12 impropers
4 = max # of 1-2 neighbors
7 = max # of 1-3 neighbors
14 = max # of 1-4 neighbors
18 = max # of special neighbors
neighbor 2.0 bin
neigh_modify delay 5
timestep 2.0
group peptide type <= 12
84 atoms in group peptide
group one id 2 4 5 6
4 atoms in group one
group two id 80 82 83 84
4 atoms in group two
group ref id 37
1 atoms in group ref
group colvar union one two ref
9 atoms in group colvar
fix 1 all nvt temp 275.0 275.0 100.0 tchain 1
fix 2 all plumed plumedfile plumed.dat outfile p.log
fix 2a ref setforce 0.0 0.0 0.0
fix 4 all shake 0.0001 10 100 b 4 6 8 10 12 14 18 a 31
19 = # of size 2 clusters
6 = # of size 3 clusters
3 = # of size 4 clusters
640 = # of frozen angles
#dump 1 colvar custom 1 dump.colvar.lammpstrj id xu yu zu fx fy fz
#dump_modify 1 sort id
thermo_style custom step temp etotal pe ke epair ebond f_2
thermo 10
variable step equal step
variable pe equal pe
run 101
PPPM initialization ...
using 12-bit tables for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.268725
grid = 15 15 15
stencil order = 5
estimated absolute RMS force accuracy = 0.0228209
estimated relative force accuracy = 6.87243e-05
using double precision FFTs
3d grid and FFT values/proc = 10648 3375
Neighbor list info ...
update every 1 steps, delay 5 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 5 5 5
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/charmm/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
SHAKE stats (type/ave/delta) on step 0
4 1.111 1.44264e-05
6 0.996998 7.26967e-06
8 1.08 1.32536e-05
10 1.111 1.22749e-05
12 1.08 1.11767e-05
14 0.96 0
18 0.957206 4.37979e-05
31 104.519 0.00396029
Per MPI rank memory allocation (min/avg/max) = 18.7 | 18.7 | 18.7 Mbytes
Step Temp TotEng PotEng KinEng E_pair E_bond f_2
0 282.10052 -5237.458 -6372.3766 1134.9186 -6442.768 16.557152 0
10 276.9783 -5234.3057 -6348.6171 1114.3114 -6421.6171 17.024361 0.47785504
20 279.08532 -5226.4036 -6349.1917 1122.7881 -6441.0169 20.764378 0.52605302
30 282.32141 -5222.3866 -6358.1939 1135.8073 -6448.9785 22.945165 0.65106011
40 276.34173 -5218.7623 -6330.5128 1111.7504 -6423.7566 15.655345 0.23795099
50 286.12741 -5215.9248 -6367.0439 1151.1192 -6449.2655 17.420975 0.42646205
60 273.01449 -5217.7381 -6316.1026 1098.3646 -6406.4709 21.800931 0.92327815
70 274.67549 -5221.0246 -6326.0716 1105.047 -6409.7721 19.41235 0.0016975896
80 273.74824 -5224.7613 -6326.0778 1101.3165 -6418.5055 19.206793 0.48550348
90 284.32594 -5229.195 -6373.0667 1143.8717 -6461.3467 21.124789 0.5468014
SHAKE stats (type/ave/delta) on step 100
4 1.111 2.06868e-06
6 0.996999 2.09521e-06
8 1.08 1.10835e-06
10 1.111 2.46599e-06
12 1.08 8.86314e-07
14 0.959999 0
18 0.9572 9.14098e-06
31 104.52 0.000760401
100 270.40648 -5234.9604 -6322.8327 1087.8723 -6417.73 19.666404 0.0094784372
101 270.99811 -5235.8295 -6326.082 1090.2525 -6418.8974 17.285816 0.086681332
Loop time of 2.12948 on 1 procs for 101 steps with 2004 atoms
Performance: 8.196 ns/day, 2.928 hours/ns, 47.429 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 1.757 | 1.757 | 1.757 | 0.0 | 82.51
Bond | 0.0052233 | 0.0052233 | 0.0052233 | 0.0 | 0.25
Kspace | 0.14772 | 0.14772 | 0.14772 | 0.0 | 6.94
Neigh | 0.16455 | 0.16455 | 0.16455 | 0.0 | 7.73
Comm | 0.0083704 | 0.0083704 | 0.0083704 | 0.0 | 0.39
Output | 0.00031424 | 0.00031424 | 0.00031424 | 0.0 | 0.01
Modify | 0.044411 | 0.044411 | 0.044411 | 0.0 | 2.09
Other | | 0.001851 | | | 0.09
Nlocal: 2004 ave 2004 max 2004 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 11134 ave 11134 max 11134 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 707961 ave 707961 max 707961 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 707961
Ave neighs/atom = 353.274
Ave special neighs/atom = 2.34032
Neighbor list builds = 8
Dangerous builds = 0
Total wall time: 0:00:02

162
examples/USER/plumed/log.24Oct18.peptide-plumed.g++.4 Normal file
View File

@ -0,0 +1,162 @@
LAMMPS (24 Oct 2018)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (../comm.cpp:87)
using 1 OpenMP thread(s) per MPI task
# Solvated 5-mer peptide
units real
atom_style full
pair_style lj/charmm/coul/long 8.0 10.0 10.0
bond_style harmonic
angle_style charmm
dihedral_style charmm
improper_style harmonic
kspace_style pppm 0.0001
read_data data.peptide
orthogonal box = (36.8402 41.0137 29.7681) to (64.2116 68.3851 57.1395)
1 by 2 by 2 MPI processor grid
reading atoms ...
2004 atoms
reading velocities ...
2004 velocities
scanning bonds ...
3 = max bonds/atom
scanning angles ...
6 = max angles/atom
scanning dihedrals ...
14 = max dihedrals/atom
scanning impropers ...
1 = max impropers/atom
reading bonds ...
1365 bonds
reading angles ...
786 angles
reading dihedrals ...
207 dihedrals
reading impropers ...
12 impropers
4 = max # of 1-2 neighbors
7 = max # of 1-3 neighbors
14 = max # of 1-4 neighbors
18 = max # of special neighbors
neighbor 2.0 bin
neigh_modify delay 5
timestep 2.0
group peptide type <= 12
84 atoms in group peptide
group one id 2 4 5 6
4 atoms in group one
group two id 80 82 83 84
4 atoms in group two
group ref id 37
1 atoms in group ref
group colvar union one two ref
9 atoms in group colvar
fix 1 all nvt temp 275.0 275.0 100.0 tchain 1
fix 2 all plumed plumedfile plumed.dat outfile p.log
fix 2a ref setforce 0.0 0.0 0.0
fix 4 all shake 0.0001 10 100 b 4 6 8 10 12 14 18 a 31
19 = # of size 2 clusters
6 = # of size 3 clusters
3 = # of size 4 clusters
640 = # of frozen angles
#dump 1 colvar custom 1 dump.colvar.lammpstrj id xu yu zu fx fy fz
#dump_modify 1 sort id
thermo_style custom step temp etotal pe ke epair ebond f_2
thermo 10
variable step equal step
variable pe equal pe
run 101
PPPM initialization ...
using 12-bit tables for long-range coulomb (../kspace.cpp:321)
G vector (1/distance) = 0.268725
grid = 15 15 15
stencil order = 5
estimated absolute RMS force accuracy = 0.0228209
estimated relative force accuracy = 6.87243e-05
using double precision FFTs
3d grid and FFT values/proc = 4312 960
Neighbor list info ...
update every 1 steps, delay 5 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 5 5 5
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/charmm/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
SHAKE stats (type/ave/delta) on step 0
4 1.111 1.44264e-05
6 0.996998 7.26967e-06
8 1.08 1.32536e-05
10 1.111 1.22749e-05
12 1.08 1.11767e-05
14 0.96 0
18 0.957206 4.37979e-05
31 104.519 0.00396029
Per MPI rank memory allocation (min/avg/max) = 15.65 | 15.86 | 16.05 Mbytes
Step Temp TotEng PotEng KinEng E_pair E_bond f_2
0 282.10052 -5237.458 -6372.3766 1134.9186 -6442.768 16.557152 0
10 276.9783 -5234.3057 -6348.6171 1114.3114 -6421.6171 17.024361 0.47785504
20 279.08532 -5226.4036 -6349.1917 1122.7881 -6441.0169 20.764378 0.52605302
30 282.32141 -5222.3866 -6358.1939 1135.8073 -6448.9785 22.945165 0.65106011
40 276.34173 -5218.7623 -6330.5128 1111.7504 -6423.7566 15.655345 0.23795099
50 286.12741 -5215.9248 -6367.0439 1151.1192 -6449.2655 17.420975 0.42646205
60 273.01449 -5217.7381 -6316.1026 1098.3646 -6406.4709 21.800931 0.92327815
70 274.67549 -5221.0246 -6326.0716 1105.047 -6409.7721 19.41235 0.0016975896
80 273.74824 -5224.7613 -6326.0778 1101.3165 -6418.5055 19.206793 0.48550348
90 284.32594 -5229.195 -6373.0667 1143.8717 -6461.3466 21.124789 0.5468014
SHAKE stats (type/ave/delta) on step 100
4 1.111 2.06868e-06
6 0.996999 2.09521e-06
8 1.08 1.10835e-06
10 1.111 2.46599e-06
12 1.08 8.86314e-07
14 0.959999 0
18 0.9572 9.14098e-06
31 104.52 0.000760401
100 270.40648 -5234.9604 -6322.8327 1087.8723 -6417.73 19.666404 0.009478437
101 270.99811 -5235.8295 -6326.082 1090.2525 -6418.8974 17.285816 0.086681332
Loop time of 1.16767 on 4 procs for 101 steps with 2004 atoms
Performance: 14.947 ns/day, 1.606 hours/ns, 86.497 timesteps/s
97.8% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.84633 | 0.86651 | 0.88617 | 1.6 | 74.21
Bond | 0.0010614 | 0.0027132 | 0.004288 | 3.0 | 0.23
Kspace | 0.095008 | 0.1162 | 0.13491 | 4.3 | 9.95
Neigh | 0.07834 | 0.078424 | 0.078516 | 0.0 | 6.72
Comm | 0.03314 | 0.033299 | 0.033426 | 0.1 | 2.85
Output | 0.00036979 | 0.00076199 | 0.0019338 | 0.0 | 0.07
Modify | 0.063471 | 0.064011 | 0.065312 | 0.3 | 5.48
Other | | 0.005751 | | | 0.49
Nlocal: 501 ave 512 max 492 min
Histogram: 1 0 0 1 0 1 0 0 0 1
Nghost: 6570.25 ave 6604 max 6529 min
Histogram: 1 0 0 1 0 0 0 0 1 1
Neighs: 176990 ave 181122 max 173551 min
Histogram: 1 1 0 0 0 0 1 0 0 1
Total # of neighbors = 707961
Ave neighs/atom = 353.274
Ave special neighs/atom = 2.34032
Neighbor list builds = 8
Dangerous builds = 0
Total wall time: 0:00:01

4
examples/USER/plumed/plumed.dat Normal file
View File

@ -0,0 +1,4 @@
UNITS LENGTH=A ENERGY=kcal/mol
dd: DISTANCE ATOMS=45,48
RESTRAINT ARG=dd KAPPA=2000 AT=6.0
PRINT ARG=dd FILE=colvar

103
examples/USER/plumed/reference/colvar Normal file
View File

@ -0,0 +1,103 @@
#! FIELDS time dd
0.000000 5.985554
0.002000 6.002880
0.004000 6.022015
0.006000 6.029922
0.008000 6.020103
0.010000 5.996906
0.012000 5.972734
0.014000 5.960079
0.016000 5.963714
0.018000 5.978140
0.020000 5.991813
0.022000 5.995155
0.024000 5.987021
0.026000 5.975340
0.028000 5.971456
0.030000 5.981945
0.032000 6.003550
0.034000 6.024743
0.036000 6.032990
0.038000 6.022936
0.040000 6.000131
0.042000 5.977800
0.044000 5.968692
0.046000 5.977224
0.048000 5.996934
0.050000 6.014800
0.052000 6.019586
0.054000 6.008803
0.056000 5.989809
0.058000 5.974484
0.060000 5.971140
0.062000 5.979074
0.064000 5.989379
0.066000 5.991356
0.068000 5.980176
0.070000 5.960625
0.072000 5.944401
0.074000 5.942614
0.076000 5.958402
0.078000 5.984574
0.080000 6.007964
0.082000 6.017667
0.084000 6.011795
0.086000 5.998304
0.088000 5.989405
0.090000 5.993275
0.092000 6.008545
0.094000 6.025183
0.096000 6.031186
0.098000 6.020651
0.100000 5.997952
0.102000 5.975230
0.104000 5.964757
0.106000 5.971150
0.108000 5.988568
0.110000 6.004676
0.112000 6.008731
0.114000 5.998481
0.116000 5.981406
0.118000 5.969615
0.120000 5.971827
0.122000 5.987658
0.124000 6.007888
0.126000 6.020477
0.128000 6.018377
0.130000 6.004046
0.132000 5.987682
0.134000 5.980338
0.136000 5.986534
0.138000 6.001303
0.140000 6.013589
0.142000 6.013717
0.144000 6.000028
0.146000 5.980283
0.148000 5.966836
0.150000 5.968670
0.152000 5.985459
0.154000 6.007612
0.156000 6.022374
0.158000 6.022034
0.160000 6.008851
0.162000 5.993355
0.164000 5.987212
0.166000 5.995452
0.168000 6.013111
0.170000 6.028386
0.172000 6.030387
0.174000 6.016468
0.176000 5.994191
0.178000 5.976616
0.180000 5.973983
0.182000 5.987185
0.184000 6.007275
0.186000 6.021338
0.188000 6.020837
0.190000 6.006955
0.192000 5.989433
0.194000 5.979796
0.196000 5.983601
0.198000 5.996921
0.200000 6.009310
0.202000 6.011114

57
examples/USER/plumed/reference/p.log Normal file
View File

@ -0,0 +1,57 @@
PLUMED: PLUMED is starting
PLUMED: Version: 2.4.2 (git: Unknown) compiled on Jul 11 2018 at 19:09:03
PLUMED: Please cite this paper when using PLUMED [1]
PLUMED: For further information see the PLUMED web page at http://www.plumed.org
PLUMED: Root: /Users/gareth/MD_code/lammps-permanent/lammps/lib/plumed/plumed2-2.4.2/
PLUMED: For installed feature, see /Users/gareth/MD_code/lammps-permanent/lammps/lib/plumed/plumed2-2.4.2//src/config/config.txt
PLUMED: Molecular dynamics engine: LAMMPS
PLUMED: Precision of reals: 8
PLUMED: Running over 1 node
PLUMED: Number of threads: 1
PLUMED: Cache line size: 512
PLUMED: Number of atoms: 2004
PLUMED: File suffix:
PLUMED: FILE: plumed.dat
PLUMED: Action UNITS
PLUMED: with label @0
PLUMED: length: A
PLUMED: energy: kcal/mol
PLUMED: time: ps
PLUMED: charge: e
PLUMED: mass: amu
PLUMED: using physical units
PLUMED: inside PLUMED, Boltzmann constant is 0.001987
PLUMED: Action DISTANCE
PLUMED: with label dd
PLUMED: between atoms 45 48
PLUMED: using periodic boundary conditions
PLUMED: Action RESTRAINT
PLUMED: with label @2
PLUMED: with arguments dd
PLUMED: added component to this action: @2.bias
PLUMED: at 6.000000
PLUMED: with harmonic force constant 2000.000000
PLUMED: and linear force constant 0.000000
PLUMED: added component to this action: @2.force2
PLUMED: Action PRINT
PLUMED: with label @3
PLUMED: with stride 1
PLUMED: with arguments dd
PLUMED: on file colvar
PLUMED: with format %f
PLUMED: END FILE: plumed.dat
PLUMED: Timestep: 0.002000
PLUMED: KbT has not been set by the MD engine
PLUMED: It should be set by hand where needed
PLUMED: Relevant bibliography:
PLUMED: [1] Tribello, Bonomi, Branduardi, Camilloni, and Bussi, Comput. Phys. Commun. 185, 604 (2014)
PLUMED: Please read and cite where appropriate!
PLUMED: Finished setup
PLUMED: Cycles Total Average Minumum Maximum
PLUMED: 1 0.020354 0.020354 0.020354 0.020354
PLUMED: 1 Prepare dependencies 102 0.000256 0.000003 0.000001 0.000006
PLUMED: 2 Sharing data 102 0.010002 0.000098 0.000078 0.000546
PLUMED: 3 Waiting for data 102 0.001398 0.000014 0.000011 0.000072
PLUMED: 4 Calculating (forward loop) 102 0.001797 0.000018 0.000013 0.000058
PLUMED: 5 Applying (backward loop) 102 0.002666 0.000026 0.000022 0.000062
PLUMED: 6 Update 102 0.001126 0.000011 0.000007 0.000055

122
lib/colvars/colvar.cpp
View File

@ -426,8 +426,8 @@ int colvar::init_custom_function(std::string const &conf)
if (x.size() != value_evaluators.size()) {
cvm::error("Error: based on custom function type, expected "
+ cvm::to_str(x.size()) + " scalar expressions, but "
+ cvm::to_str(value_evaluators.size() + " were found.\n"));
+ cvm::to_str(x.size()) + " scalar expressions, but "
+ cvm::to_str(value_evaluators.size()) + " were found.\n");
return INPUT_ERROR;
}
@ -455,36 +455,42 @@ int colvar::init_grid_parameters(std::string const &conf)
}
lower_boundary.type(value());
upper_boundary.type(value());
upper_wall.type(value());
if (is_enabled(f_cv_scalar)) {
if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {
enable(f_cv_lower_boundary);
}
std::string lw_conf, uw_conf;
if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0, parse_silent)) {
cvm::log("Warning: lowerWallConstant and lowerWall are deprecated, "
"please define a harmonicWalls bias instead.\n");
lower_wall.type(value());
get_keyval(conf, "lowerWall", lower_wall, lower_boundary);
lw_conf = std::string("\n\
lowerWallConstant "+cvm::to_str(lower_wall_k*width*width)+"\n\
lowerWalls "+cvm::to_str(lower_wall)+"\n");
}
if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) {
enable(f_cv_upper_boundary);
}
if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0, parse_silent)) {
cvm::log("Warning: upperWallConstant and upperWall are deprecated, "
"please define a harmonicWalls bias instead.\n");
std::string lw_conf, uw_conf;
if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0,
parse_silent)) {
cvm::log("Reading legacy options lowerWall and lowerWallConstant: "
"consider using a harmonicWalls restraint.\n");
lower_wall.type(value());
if (!get_keyval(conf, "lowerWall", lower_wall, lower_boundary)) {
cvm::log("Warning: lowerWall will need to be "
"defined explicitly in the next release.\n");
}
lw_conf = std::string("\n\
lowerWallConstant "+cvm::to_str(lower_wall_k*width*width)+"\n\
lowerWalls "+cvm::to_str(lower_wall)+"\n");
}
if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0,
parse_silent)) {
cvm::log("Reading legacy options upperWall and upperWallConstant: "
"consider using a harmonicWalls restraint.\n");
upper_wall.type(value());
get_keyval(conf, "upperWall", upper_wall, upper_boundary);
if (!get_keyval(conf, "upperWall", upper_wall, upper_boundary)) {
cvm::log("Warning: upperWall will need to be "
"defined explicitly in the next release.\n");
}
uw_conf = std::string("\n\
upperWallConstant "+cvm::to_str(upper_wall_k*width*width)+"\n\
upperWalls "+cvm::to_str(upper_wall)+"\n");
@ -677,6 +683,7 @@ template<typename def_class_name> int colvar::init_components_type(std::string c
if (cvcp != NULL) {
cvcs.push_back(cvcp);
cvcp->check_keywords(def_conf, def_config_key);
cvcp->config_key = def_config_key;
if (cvm::get_error()) {
cvm::error("Error: in setting up component \""+
std::string(def_config_key)+"\".\n", INPUT_ERROR);
@ -1022,13 +1029,13 @@ int colvar::calc()
int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
{
colvarproxy *proxy = cvm::main()->proxy;
int error_code = COLVARS_OK;
if (cvm::debug())
cvm::log("Calculating colvar \""+this->name+"\", components "+
cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");
colvarproxy *proxy = cvm::main()->proxy;
int error_code = COLVARS_OK;
error_code |= check_cvc_range(first_cvc, num_cvcs);
if (error_code != COLVARS_OK) {
return error_code;
@ -1059,9 +1066,10 @@ int colvar::collect_cvc_data()
if (cvm::debug())
cvm::log("Calculating colvar \""+this->name+"\"'s properties.\n");
colvarproxy *proxy = cvm::main()->proxy;
int error_code = COLVARS_OK;
if (cvm::step_relative() > 0) {
if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
// Total force depends on Jacobian derivative from previous timestep
// collect_cvc_total_forces() uses the previous value of jd
error_code |= collect_cvc_total_forces();
@ -1069,6 +1077,10 @@ int colvar::collect_cvc_data()
error_code |= collect_cvc_values();
error_code |= collect_cvc_gradients();
error_code |= collect_cvc_Jacobians();
if (proxy->total_forces_same_step()){
// Use Jacobian derivative from this timestep
error_code |= collect_cvc_total_forces();
}
error_code |= calc_colvar_properties();
if (cvm::debug())
@ -1394,6 +1406,13 @@ int colvar::calc_colvar_properties()
vr.reset(); // (already 0; added for clarity)
}
// Special case of a repeated timestep (eg. multiple NAMD "run" statements)
// revert values of the extended coordinate and velocity prior to latest integration
if (cvm::step_relative() == prev_timestep) {
xr = prev_xr;
vr = prev_vr;
}
// report the restraint center as "value"
x_reported = xr;
v_reported = vr;
@ -1450,7 +1469,6 @@ cvm::real colvar::update_forces_energy()
// extended variable if there is one
if (is_enabled(f_cv_extended_Lagrangian)) {
if (cvm::debug()) {
cvm::log("Updating extended-Lagrangian degree of freedom.\n");
}
@ -1502,6 +1520,11 @@ cvm::real colvar::update_forces_energy()
ft_reported = f_ext;
}
// backup in case we need to revert this integration timestep
// if the same MD timestep is re-run
prev_xr = xr;
prev_vr = vr;
// leapfrog: starting from x_i, f_i, v_(i-1/2)
vr += (0.5 * dt) * f_ext / ext_mass;
// Because of leapfrog, kinetic energy at time i is approximate
@ -1638,18 +1661,26 @@ int colvar::set_cvc_flags(std::vector<bool> const &flags)
}
void colvar::update_active_cvc_square_norm()
{
active_cvc_square_norm = 0.0;
for (size_t i = 0; i < cvcs.size(); i++) {
if (cvcs[i]->is_enabled()) {
active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;
}
}
}
int colvar::update_cvc_flags()
{
// Update the enabled/disabled status of cvcs if necessary
if (cvc_flags.size()) {
n_active_cvcs = 0;
active_cvc_square_norm = 0.;
for (size_t i = 0; i < cvcs.size(); i++) {
cvcs[i]->set_enabled(f_cvc_active, cvc_flags[i]);
if (cvcs[i]->is_enabled()) {
n_active_cvcs++;
active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;
}
}
if (!n_active_cvcs) {
@ -1657,12 +1688,49 @@ int colvar::update_cvc_flags()
return COLVARS_ERROR;
}
cvc_flags.clear();
update_active_cvc_square_norm();
}
return COLVARS_OK;
}
int colvar::update_cvc_config(std::vector<std::string> const &confs)
{
cvm::log("Updating configuration for colvar \""+name+"\"");
if (confs.size() != cvcs.size()) {
return cvm::error("Error: Wrong number of CVC config strings. "
"For those CVCs that are not being changed, try passing "
"an empty string.", INPUT_ERROR);
}
int error_code = COLVARS_OK;
int num_changes = 0;
for (size_t i = 0; i < cvcs.size(); i++) {
if (confs[i].size()) {
std::string conf(confs[i]);
cvm::increase_depth();
error_code |= cvcs[i]->colvar::cvc::init(conf);
error_code |= cvcs[i]->check_keywords(conf,
cvcs[i]->config_key.c_str());
cvm::decrease_depth();
num_changes++;
}
}
if (num_changes == 0) {
cvm::log("Warning: no changes were applied through modifycvcs; "
"please check that its argument is a list of strings.\n");
}
update_active_cvc_square_norm();
return error_code;
}
// ******************** METRIC FUNCTIONS ********************
// Use the metrics defined by \link cvc \endlink objects

20
lib/colvars/colvar.h
View File

@ -171,8 +171,12 @@ protected:
// Options for extended_lagrangian
/// Restraint center
colvarvalue xr;
/// Previous value of the restraint center;
colvarvalue prev_xr;
/// Velocity of the restraint center
colvarvalue vr;
/// Previous velocity of the restraint center
colvarvalue prev_vr;
/// Mass of the restraint center
cvm::real ext_mass;
/// Restraint force constant
@ -352,6 +356,9 @@ public:
/// \brief Updates the flags in the CVC objects, and their number
int update_cvc_flags();
/// \brief Modify the configuration of CVCs (currently, only base class data)
int update_cvc_config(std::vector<std::string> const &confs);
protected:
/// \brief Number of CVC objects with an active flag
size_t n_active_cvcs;
@ -359,10 +366,17 @@ protected:
/// Sum of square coefficients for active cvcs
cvm::real active_cvc_square_norm;
/// Update the sum of square coefficients for active cvcs
void update_active_cvc_square_norm();
/// \brief Absolute timestep number when this colvar was last updated
int prev_timestep;
public:
/// \brief Return the number of CVC objects defined
inline size_t num_cvcs() const { return cvcs.size(); }
/// \brief Return the number of CVC objects with an active flag (as set by update_cvc_flags)
inline size_t num_active_cvcs() const { return n_active_cvcs; }
@ -371,21 +385,21 @@ public:
///
/// Handles correctly symmetries and periodic boundary conditions
cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
///
/// Handles correctly symmetries and periodic boundary conditions
colvarvalue dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
///
/// Handles correctly symmetries and periodic boundary conditions
colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to wrap a value into a standard interval

4
lib/colvars/colvaratoms.cpp
View File

@ -21,6 +21,8 @@ cvm::atom::atom()
{
index = -1;
id = -1;
mass = 1.0;
charge = 1.0;
reset_data();
}
@ -395,7 +397,7 @@ int cvm::atom_group::parse(std::string const &group_conf)
}
// NOTE: calls to add_atom() and/or add_atom_id() are in the proxy-implemented function
cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
parse_error |= cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
}
}

2
lib/colvars/colvaratoms.h
View File

@ -90,7 +90,7 @@ public:
/// Destructor
~atom();
/// Set mutable data (everything except id and mass) to zero; update mass
/// Set mutable data (everything except id and mass) to zero
inline void reset_data()
{
pos = cvm::atom_pos(0.0);

58
lib/colvars/colvarbias_abf.cpp
View File

@ -564,6 +564,8 @@ int colvarbias_abf::replica_share() {
return COLVARS_OK;
}
void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool append)
{
std::string samples_out_name = prefix + ".count";
@ -572,10 +574,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::ostream *samples_os =
cvm::proxy->output_stream(samples_out_name, mode);
if (!samples_os) {
cvm::error("Error opening ABF samples file " + samples_out_name + " for writing");
return;
}
if (!samples_os) return;
samples->write_multicol(*samples_os);
cvm::proxy->close_output_stream(samples_out_name);
@ -583,10 +582,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
if (num_variables() > 2) {
std::string samples_dx_out_name = prefix + ".count.dx";
std::ostream *samples_dx_os = cvm::proxy->output_stream(samples_dx_out_name, mode);
if (!samples_os) {
cvm::error("Error opening samples file " + samples_dx_out_name + " for writing");
return;
}
if (!samples_os) return;
samples->write_opendx(*samples_dx_os);
*samples_dx_os << std::endl;
cvm::proxy->close_output_stream(samples_dx_out_name);
@ -594,10 +590,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::ostream *gradients_os =
cvm::proxy->output_stream(gradients_out_name, mode);
if (!gradients_os) {
cvm::error("Error opening ABF gradient file " + gradients_out_name + " for writing");
return;
}
if (!gradients_os) return;
gradients->write_multicol(*gradients_os);
cvm::proxy->close_output_stream(gradients_out_name);
@ -609,20 +602,14 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::string pmf_out_name = prefix + ".pmf";
std::ostream *pmf_os = cvm::proxy->output_stream(pmf_out_name, mode);
if (!pmf_os) {
cvm::error("Error opening pmf file " + pmf_out_name + " for writing");
return;
}
if (!pmf_os) return;
pmf->write_multicol(*pmf_os);
// In dimension higher than 2, dx is easier to handle and visualize
if (num_variables() > 2) {
std::string pmf_dx_out_name = prefix + ".pmf.dx";
std::ostream *pmf_dx_os = cvm::proxy->output_stream(pmf_dx_out_name, mode);
if (!pmf_dx_os) {
cvm::error("Error opening pmf file " + pmf_dx_out_name + " for writing");
return;
}
if (!pmf_dx_os) return;
pmf->write_opendx(*pmf_dx_os);
*pmf_dx_os << std::endl;
cvm::proxy->close_output_stream(pmf_dx_out_name);
@ -639,10 +626,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::ostream *z_samples_os =
cvm::proxy->output_stream(z_samples_out_name, mode);
if (!z_samples_os) {
cvm::error("Error opening eABF z-histogram file " + z_samples_out_name + " for writing");
return;
}
if (!z_samples_os) return;
z_samples->write_multicol(*z_samples_os);
cvm::proxy->close_output_stream(z_samples_out_name);
@ -651,10 +635,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::ostream *z_gradients_os =
cvm::proxy->output_stream(z_gradients_out_name, mode);
if (!z_gradients_os) {
cvm::error("Error opening eABF z-gradient file " + z_gradients_out_name + " for writing");
return;
}
if (!z_gradients_os) return;
z_gradients->write_multicol(*z_gradients_os);
cvm::proxy->close_output_stream(z_gradients_out_name);
}
@ -672,10 +653,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::ostream *czar_gradients_os =
cvm::proxy->output_stream(czar_gradients_out_name, mode);
if (!czar_gradients_os) {
cvm::error("Error opening CZAR gradient file " + czar_gradients_out_name + " for writing");
return;
}
if (!czar_gradients_os) return;
czar_gradients->write_multicol(*czar_gradients_os);
cvm::proxy->close_output_stream(czar_gradients_out_name);
@ -688,20 +666,14 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::string czar_pmf_out_name = prefix + ".czar.pmf";
std::ostream *czar_pmf_os = cvm::proxy->output_stream(czar_pmf_out_name, mode);
if (!czar_pmf_os) {
cvm::error("Error opening CZAR pmf file " + czar_pmf_out_name + " for writing");
return;
}
if (!czar_pmf_os) return;
czar_pmf->write_multicol(*czar_pmf_os);
// In dimension higher than 2, dx is easier to handle and visualize
if (num_variables() > 2) {
std::string czar_pmf_dx_out_name = prefix + ".czar.pmf.dx";
std::ostream *czar_pmf_dx_os = cvm::proxy->output_stream(czar_pmf_dx_out_name, mode);
if (!czar_pmf_dx_os) {
cvm::error("Error opening CZAR pmf file " + czar_pmf_dx_out_name + " for writing");
return;
}
if (!czar_pmf_dx_os) return;
czar_pmf->write_opendx(*czar_pmf_dx_os);
*czar_pmf_dx_os << std::endl;
cvm::proxy->close_output_stream(czar_pmf_dx_out_name);
@ -854,3 +826,9 @@ std::istream & colvarbias_abf::read_state_data(std::istream& is)
return is;
}
int colvarbias_abf::write_output_files()
{
write_gradients_samples(output_prefix);
return COLVARS_OK;
}

6
lib/colvars/colvarbias_abf.h
View File

@ -136,13 +136,13 @@ private:
/// Write human-readable FE gradients and sample count, and DX file in dim > 2
void write_gradients_samples(const std::string &prefix, bool append = false);
void write_last_gradients_samples(const std::string &prefix, bool append = false);
/// Read human-readable FE gradients and sample count (if not using restart)
void read_gradients_samples();
std::istream& read_state_data(std::istream&);
std::ostream& write_state_data(std::ostream&);
virtual std::istream& read_state_data(std::istream&);
virtual std::ostream& write_state_data(std::ostream&);
virtual int write_output_files();
};
#endif

22
lib/colvars/colvarbias_restraint.cpp
View File

@ -996,11 +996,16 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
if ((lower_walls.size() > 0) && (upper_walls.size() > 0)) {
for (i = 0; i < num_variables(); i++) {
if (lower_walls[i] >= upper_walls[i]) {
cvm::error("Error: one upper wall, "+
cvm::to_str(upper_walls[i])+
", is not higher than the lower wall, "+
cvm::to_str(lower_walls[i])+".\n",
INPUT_ERROR);
return cvm::error("Error: one upper wall, "+
cvm::to_str(upper_walls[i])+
", is not higher than the lower wall, "+
cvm::to_str(lower_walls[i])+".\n",
INPUT_ERROR);
}
if (variables(i)->dist2(lower_walls[i], upper_walls[i]) < 1.0e-12) {
return cvm::error("Error: lower wall and upper wall are equal "
"in the domain of the variable \""+
variables(i)->name+"\".\n", INPUT_ERROR);
}
}
if (lower_wall_k * upper_wall_k == 0.0) {
@ -1279,13 +1284,16 @@ cvm::real colvarbias_restraint_linear::energy_difference(std::string const &conf
cvm::real colvarbias_restraint_linear::restraint_potential(size_t i) const
{
return force_k / variables(i)->width * (variables(i)->value() - colvar_centers[i]);
return force_k / variables(i)->width * (variables(i)->value() -
colvar_centers[i]).sum();
}
colvarvalue const colvarbias_restraint_linear::restraint_force(size_t i) const
{
return -1.0 * force_k / variables(i)->width;
colvarvalue dummy(variables(i)->value());
dummy.set_ones();
return -1.0 * force_k / variables(i)->width * dummy;
}

25
lib/colvars/colvarcomp.cpp
View File

@ -43,16 +43,27 @@ colvar::cvc::cvc(std::string const &conf)
int colvar::cvc::init(std::string const &conf)
{
int error_code = COLVARS_OK;
if (cvm::debug())
cvm::log("Initializing cvc base object.\n");
get_keyval(conf, "name", this->name, this->name);
std::string const old_name(name);
if (name.size() > 0) {
// Temporary description until child object is initialized
description = "cvc " + name;
} else {
description = "uninitialized cvc";
cvm::log("Updating configuration for component \""+name+"\"");
}
if (get_keyval(conf, "name", name, name)) {
if (name.size() > 0) {
description = "cvc \"" + name + "\" of type " + function_type;
} else {
description = "unnamed cvc";
}
if ((name != old_name) && (old_name.size() > 0)) {
cvm::error("Error: cannot rename component \""+old_name+
"\" after initialization (new name = \""+name+"\")",
INPUT_ERROR);
name = old_name;
}
}
get_keyval(conf, "componentCoeff", sup_coeff, sup_coeff);
@ -78,7 +89,7 @@ int colvar::cvc::init(std::string const &conf)
if (cvm::debug())
cvm::log("Done initializing cvc base object.\n");
return error_code;
return cvm::get_error();
}

118
lib/colvars/colvarcomp.h
View File

@ -82,6 +82,9 @@ public:
/// this variable definition should be set within the constructor.
std::string function_type;
/// Keyword used in the input to denote this CVC
std::string config_key;
/// \brief Coefficient in the polynomial combination (default: 1.0)
cvm::real sup_coeff;
/// \brief Exponent in the polynomial combination (default: 1)
@ -834,49 +837,70 @@ protected:
cvm::real r0;
/// \brief "Cutoff vector" for anisotropic calculation
cvm::rvector r0_vec;
/// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be
/// used
/// \brief Whether r/r0 or \vec{r}*\vec{1/r0_vec} should be used
bool b_anisotropic;
/// Integer exponent of the function numerator
int en;
/// Integer exponent of the function denominator
int ed;
/// \brief If true, group2 will be treated as a single atom
/// (default: loop over all pairs of atoms in group1 and group2)
bool b_group2_center_only;
/// \brief If true, group2 will be treated as a single atom, stored in this
/// accessory group
cvm::atom_group *group2_center;
/// Tolerance for the pair list
cvm::real tolerance;
/// Frequency of update of the pair list
int pairlist_freq;
/// Pair list
bool *pairlist;
public:
/// Constructor
coordnum(std::string const &conf);
coordnum();
virtual ~coordnum() {}
~coordnum();
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
template<bool b_gradients>
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
/// coordination number \param exp_num \i n exponent \param exp_den
/// \i m exponent \param A1 atom \param A2 atom
static cvm::real switching_function(cvm::real const &r0,
int const &exp_num, int const &exp_den,
cvm::atom &A1, cvm::atom &A2);
template<bool b_gradients>
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
/// vector of different cutoffs in the three directions \param
/// exp_num \i n exponent \param exp_den \i m exponent \param A1
/// atom \param A2 atom
static cvm::real switching_function(cvm::rvector const &r0_vec,
int const &exp_num, int const &exp_den,
cvm::atom &A1, cvm::atom &A2);
virtual cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;
virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
enum {
ef_null = 0,
ef_gradients = 1,
ef_anisotropic = (1<<8),
ef_use_pairlist = (1<<9),
ef_rebuild_pairlist = (1<<10)
};
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), where x = |A1-A2|/r0 \param r0, r0_vec "cutoff" for
/// the coordination number (scalar or vector depending on user choice)
/// \param en Numerator exponent \param ed Denominator exponent \param First
/// atom \param Second atom \param pairlist_elem pointer to pair flag for
/// this pair \param tolerance A pair is defined as having a larger
/// coordination than this number
template<int flags>
static cvm::real switching_function(cvm::real const &r0,
cvm::rvector const &r0_vec,
int en,
int ed,
cvm::atom &A1,
cvm::atom &A2,
bool **pairlist_elem,
cvm::real tolerance);
/// Main workhorse function
template<int flags> int compute_coordnum();
};
@ -887,7 +911,8 @@ class colvar::selfcoordnum
: public colvar::cvc
{
protected:
/// First atom group
/// Selected atoms
cvm::atom_group *group1;
/// \brief "Cutoff" for isotropic calculation (default)
cvm::real r0;
@ -895,22 +920,18 @@ protected:
int en;
/// Integer exponent of the function denominator
int ed;
cvm::real tolerance;
int pairlist_freq;
bool *pairlist;
public:
/// Constructor
selfcoordnum(std::string const &conf);
selfcoordnum();
virtual ~selfcoordnum() {}
~selfcoordnum();
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
template<bool b_gradients>
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
/// coordination number \param exp_num \i n exponent \param exp_den
/// \i m exponent \param A1 atom \param A2 atom
static cvm::real switching_function(cvm::real const &r0,
int const &exp_num, int const &exp_den,
cvm::atom &A1, cvm::atom &A2);
virtual cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;
@ -918,6 +939,9 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
/// Main workhorse function
template<int flags> int compute_selfcoordnum();
};
@ -947,26 +971,6 @@ public:
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
template<bool b_gradients>
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
/// coordination number \param exp_num \i n exponent \param exp_den
/// \i m exponent \param A1 atom \param A2 atom
static cvm::real switching_function(cvm::real const &r0,
int const &exp_num, int const &exp_den,
cvm::atom &A1, cvm::atom &A2);
/*
template<bool b_gradients>
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
/// vector of different cutoffs in the three directions \param
/// exp_num \i n exponent \param exp_den \i m exponent \param A1
/// atom \param A2 atom
static cvm::real switching_function(cvm::rvector const &r0_vec,
int const &exp_num, int const &exp_den,
cvm::atom &A1, cvm::atom &A2);
*/
virtual cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;

521
lib/colvars/colvarcomp_coordnums.cpp
View File

@ -18,45 +18,36 @@
template<bool calculate_gradients>
template<int flags>
cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
int const &en,
int const &ed,
cvm::rvector const &r0_vec,
int en,
int ed,
cvm::atom &A1,
cvm::atom &A2)
cvm::atom &A2,
bool **pairlist_elem,
cvm::real pairlist_tol)
{
cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
cvm::real const l2 = diff.norm2()/(r0*r0);
// Assume en and ed are even integers, and avoid sqrt in the following
int const en2 = en/2;
int const ed2 = ed/2;
cvm::real const xn = cvm::integer_power(l2, en2);
cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
if (calculate_gradients) {
cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
A1.grad += (-1.0)*dFdl2*dl2dx;
A2.grad += dFdl2*dl2dx;
if ((flags & ef_use_pairlist) && !(flags & ef_rebuild_pairlist)) {
bool const within = **pairlist_elem;
(*pairlist_elem)++;
if (!within) {
return 0.0;
}
}
return func;
}
cvm::rvector const r0sq_vec(r0_vec.x*r0_vec.x,
r0_vec.y*r0_vec.y,
r0_vec.z*r0_vec.z);
template<bool calculate_gradients>
cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
int const &en,
int const &ed,
cvm::atom &A1,
cvm::atom &A2)
{
cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
cvm::rvector const scal_diff(diff.x/((flags & ef_anisotropic) ?
r0_vec.x : r0),
diff.y/((flags & ef_anisotropic) ?
r0_vec.y : r0),
diff.z/((flags & ef_anisotropic) ?
r0_vec.z : r0));
cvm::real const l2 = scal_diff.norm2();
// Assume en and ed are even integers, and avoid sqrt in the following
@ -65,22 +56,45 @@ cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
cvm::real const xn = cvm::integer_power(l2, en2);
cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
//The subtraction and division stretches the function back to the range of [0,1] from [pairlist_tol,1]
cvm::real const func = (((1.0-xn)/(1.0-xd)) - pairlist_tol) / (1.0-pairlist_tol);
if (calculate_gradients) {
cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
(2.0/(r0_vec.y*r0_vec.y))*diff.y,
(2.0/(r0_vec.z*r0_vec.z))*diff.z);
if (flags & ef_rebuild_pairlist) {
//Particles just outside of the cutoff also are considered if they come near.
**pairlist_elem = (func > (-pairlist_tol * 0.5)) ? true : false;
(*pairlist_elem)++;
}
//If the value is too small, we need to exclude it, rather than let it contribute to the sum or the gradients.
if (func < 0)
return 0;
if (flags & ef_gradients) {
//This is the old, completely correct expression for dFdl2:
//cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) -
// func*ed2*(xd/l2))*(-1.0);
//This can become:
//cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2)*(1.0-xn)/(1.0-xn) -
// func*ed2*(xd/l2))*(-1.0);
//Recognizing that func = (1.0-xn)/(1.0-xd), we can group together the "func" and get a version of dFdl2 that is 0
//when func=0, which lets us skip this gradient calculation when func=0.
cvm::real const dFdl2 = func * ((ed2*xd/((1.0-xd)*l2)) - (en2*xn/((1.0-xn)*l2)));
cvm::rvector const dl2dx((2.0/((flags & ef_anisotropic) ? r0sq_vec.x :
r0*r0)) * diff.x,
(2.0/((flags & ef_anisotropic) ? r0sq_vec.y :
r0*r0)) * diff.y,
(2.0/((flags & ef_anisotropic) ? r0sq_vec.z :
r0*r0)) * diff.z);
A1.grad += (-1.0)*dFdl2*dl2dx;
A2.grad += dFdl2*dl2dx;
}
return func;
}
colvar::coordnum::coordnum(std::string const &conf)
: cvc(conf), b_anisotropic(false), b_group2_center_only(false)
: cvc(conf), b_anisotropic(false), group2_center(NULL), pairlist(NULL)
{
function_type = "coordnum";
x.type(colvarvalue::type_scalar);
@ -90,23 +104,26 @@ colvar::coordnum::coordnum(std::string const &conf)
if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) {
cvm::error("Error: group1 and group2 share a common atom (number: " +
cvm::to_str(atom_number) + ")\n");
cvm::to_str(atom_number) + ")\n", INPUT_ERROR);
return;
}
if (group1->b_dummy) {
cvm::error("Error: only group2 is allowed to be a dummy atom\n");
cvm::error("Error: only group2 is allowed to be a dummy atom\n",
INPUT_ERROR);
return;
}
bool const b_isotropic = get_keyval(conf, "cutoff", r0,
cvm::real(4.0 * cvm::unit_angstrom()));
if (get_keyval(conf, "cutoff3", r0_vec, cvm::rvector(4.0 * cvm::unit_angstrom(),
4.0 * cvm::unit_angstrom(),
4.0 * cvm::unit_angstrom()))) {
if (get_keyval(conf, "cutoff3", r0_vec,
cvm::rvector(4.0 * cvm::unit_angstrom(),
4.0 * cvm::unit_angstrom(),
4.0 * cvm::unit_angstrom()))) {
if (b_isotropic) {
cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" at the same time.\n",
cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" "
"at the same time.\n",
INPUT_ERROR);
return;
}
@ -135,86 +152,178 @@ colvar::coordnum::coordnum(std::string const &conf)
cvm::log("Warning: only minimum-image distances are used by this variable.\n");
}
bool b_group2_center_only = false;
get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy);
if (b_group2_center_only) {
if (!group2_center) {
group2_center = new cvm::atom_group();
group2_center->add_atom(cvm::atom());
}
}
get_keyval(conf, "tolerance", tolerance, 0.0);
if (tolerance > 0) {
get_keyval(conf, "pairListFrequency", pairlist_freq, 100);
if ( ! (pairlist_freq > 0) ) {
cvm::error("Error: non-positive pairlistfrequency provided.\n",
INPUT_ERROR);
return; // and do not allocate the pairlists below
}
if (b_group2_center_only) {
pairlist = new bool[group1->size()];
}
else {
pairlist = new bool[group1->size() * group2->size()];
}
}
}
colvar::coordnum::coordnum()
: b_anisotropic(false), b_group2_center_only(false)
: b_anisotropic(false), group2_center(NULL), pairlist(NULL)
{
function_type = "coordnum";
x.type(colvarvalue::type_scalar);
}
void colvar::coordnum::calc_value()
colvar::coordnum::~coordnum()
{
x.real_value = 0.0;
if (pairlist != NULL) {
delete [] pairlist;
}
if (group2_center != NULL) {
delete group2_center;
}
}
if (b_group2_center_only) {
// create a fake atom to hold the group2 com coordinates
cvm::atom group2_com_atom;
group2_com_atom.pos = group2->center_of_mass();
template<int compute_flags> int colvar::coordnum::compute_coordnum()
{
if (group2_center) {
(*group2_center)[0].pos = group2->center_of_mass();
group2_center->calc_required_properties();
}
cvm::atom_group *group2p = group2_center ? group2_center : group2;
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, group2_com_atom);
} else {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
x.real_value += switching_function<false>(r0, en, ed, *ai1, group2_com_atom);
bool const use_pairlist = (pairlist != NULL);
bool const rebuild_pairlist = (pairlist != NULL) &&
(cvm::step_relative() % pairlist_freq == 0);
bool *pairlist_elem = use_pairlist ? pairlist : NULL;
cvm::atom_iter ai1 = group1->begin(), ai2 = group2p->begin();
cvm::atom_iter const ai1_end = group1->end();
cvm::atom_iter const ai2_end = group2p->end();
if (b_anisotropic) {
if (use_pairlist) {
if (rebuild_pairlist) {
int const flags = compute_flags | ef_anisotropic | ef_use_pairlist |
ef_rebuild_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
} else {
int const flags = compute_flags | ef_anisotropic | ef_use_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
}
} else { // if (use_pairlist) {
int const flags = compute_flags | ef_anisotropic;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
NULL, 0.0);
}
}
}
} else {
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, *ai2);
if (use_pairlist) {
if (rebuild_pairlist) {
int const flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
} else {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
x.real_value += switching_function<false>(r0, en, ed, *ai1, *ai2);
} else {
int const flags = compute_flags | ef_use_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
}
} else { // if (use_pairlist) {
int const flags = compute_flags;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
NULL, 0.0);
}
}
}
}
if (compute_flags & ef_gradients) {
if (group2_center) {
group2->set_weighted_gradient((*group2_center)[0].grad);
}
}
return COLVARS_OK;
}
void colvar::coordnum::calc_value()
{
x.real_value = 0.0;
if (is_enabled(f_cvc_gradient)) {
compute_coordnum<ef_gradients>();
} else {
compute_coordnum<ef_null>();
}
}
void colvar::coordnum::calc_gradients()
{
if (b_group2_center_only) {
// create a fake atom to hold the group2 com coordinates
cvm::atom group2_com_atom;
group2_com_atom.pos = group2->center_of_mass();
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
switching_function<true>(r0_vec, en, ed, *ai1, group2_com_atom);
} else {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
switching_function<true>(r0, en, ed, *ai1, group2_com_atom);
}
group2->set_weighted_gradient(group2_com_atom.grad);
} else {
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
switching_function<true>(r0_vec, en, ed, *ai1, *ai2);
}
} else {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
switching_function<true>(r0, en, ed, *ai1, *ai2);
}
}
}
// Gradients are computed by calc_value() if f_cvc_gradients is enabled
}
@ -235,7 +344,7 @@ simple_scalar_dist_functions(coordnum)
// h_bond member functions
colvar::h_bond::h_bond(std::string const &conf)
: cvc(conf)
: cvc(conf)
{
if (cvm::debug())
cvm::log("Initializing h_bond object.\n");
@ -307,13 +416,24 @@ colvar::h_bond::~h_bond()
void colvar::h_bond::calc_value()
{
x.real_value = colvar::coordnum::switching_function<false>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
int const flags = coordnum::ef_null;
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
x.real_value =
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
(*atom_groups[0])[0],
(*atom_groups[0])[1],
NULL, 0.0);
}
void colvar::h_bond::calc_gradients()
{
colvar::coordnum::switching_function<true>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
int const flags = coordnum::ef_gradients;
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
(*atom_groups[0])[0],
(*atom_groups[0])[1],
NULL, 0.0);
}
@ -328,7 +448,7 @@ simple_scalar_dist_functions(h_bond)
colvar::selfcoordnum::selfcoordnum(std::string const &conf)
: cvc(conf)
: cvc(conf), pairlist(NULL)
{
function_type = "selfcoordnum";
x.type(colvarvalue::type_scalar);
@ -353,36 +473,115 @@ colvar::selfcoordnum::selfcoordnum(std::string const &conf)
if (!is_enabled(f_cvc_pbc_minimum_image)) {
cvm::log("Warning: only minimum-image distances are used by this variable.\n");
}
get_keyval(conf, "tolerance", tolerance, 0.0);
if (tolerance > 0) {
get_keyval(conf, "pairListFrequency", pairlist_freq, 100);
if ( ! (pairlist_freq > 0) ) {
cvm::error("Error: non-positive pairlistfrequency provided.\n",
INPUT_ERROR);
return;
}
pairlist = new bool[(group1->size()-1) * (group1->size()-1)];
}
}
colvar::selfcoordnum::selfcoordnum()
: pairlist(NULL)
{
function_type = "selfcoordnum";
x.type(colvarvalue::type_scalar);
}
colvar::selfcoordnum::~selfcoordnum()
{
if (pairlist != NULL) {
delete [] pairlist;
}
}
template<int compute_flags> int colvar::selfcoordnum::compute_selfcoordnum()
{
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
bool const use_pairlist = (pairlist != NULL);
bool const rebuild_pairlist = (pairlist != NULL) &&
(cvm::step_relative() % pairlist_freq == 0);
bool *pairlist_elem = use_pairlist ? pairlist : NULL;
size_t i = 0, j = 0;
size_t const n = group1->size();
// Always isotropic (TODO: enable the ellipsoid?)
if (use_pairlist) {
if (rebuild_pairlist) {
int const flags = compute_flags | coordnum::ef_use_pairlist |
coordnum::ef_rebuild_pairlist;
for (i = 0; i < n - 1; i++) {
for (j = i + 1; j < n; j++) {
x.real_value +=
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
(*group1)[i],
(*group1)[j],
&pairlist_elem,
tolerance);
}
}
} else {
int const flags = compute_flags | coordnum::ef_use_pairlist;
for (i = 0; i < n - 1; i++) {
for (j = i + 1; j < n; j++) {
x.real_value +=
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
(*group1)[i],
(*group1)[j],
&pairlist_elem,
tolerance);
}
}
}
} else { // if (use_pairlist) {
int const flags = compute_flags | coordnum::ef_null;
for (i = 0; i < n - 1; i++) {
for (j = i + 1; j < n; j++) {
x.real_value +=
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
(*group1)[i],
(*group1)[j],
&pairlist_elem,
tolerance);
}
}
}
return COLVARS_OK;
}
void colvar::selfcoordnum::calc_value()
{
x.real_value = 0.0;
for (size_t i = 0; i < group1->size() - 1; i++) {
for (size_t j = i + 1; j < group1->size(); j++) {
x.real_value += colvar::coordnum::switching_function<false>(r0, en, ed, (*group1)[i], (*group1)[j]);
}
if (is_enabled(f_cvc_gradient)) {
compute_selfcoordnum<coordnum::ef_gradients>();
} else {
compute_selfcoordnum<coordnum::ef_null>();
}
}
void colvar::selfcoordnum::calc_gradients()
{
for (size_t i = 0; i < group1->size() - 1; i++) {
for (size_t j = i + 1; j < group1->size(); j++) {
colvar::coordnum::switching_function<true>(r0, en, ed, (*group1)[i], (*group1)[j]);
}
}
// Gradients are computed by calc_value() if f_cvc_gradients is enabled
}
void colvar::selfcoordnum::apply_force(colvarvalue const &force)
{
if (!group1->noforce) {
@ -394,6 +593,7 @@ void colvar::selfcoordnum::apply_force(colvarvalue const &force)
simple_scalar_dist_functions(selfcoordnum)
// groupcoordnum member functions
colvar::groupcoordnum::groupcoordnum(std::string const &conf)
: distance(conf), b_anisotropic(false)
@ -415,7 +615,7 @@ colvar::groupcoordnum::groupcoordnum(std::string const &conf)
if (b_scale) {
cvm::error("Error: cannot specify \"scale\" and "
"\"scale3\" at the same time.\n");
"\"scale3\" at the same time.\n");
return;
}
b_anisotropic = true;
@ -453,95 +653,56 @@ colvar::groupcoordnum::groupcoordnum()
}
template<bool calculate_gradients>
cvm::real colvar::groupcoordnum::switching_function(cvm::real const &r0,
int const &en,
int const &ed,
cvm::atom &A1,
cvm::atom &A2)
{
cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
cvm::real const l2 = diff.norm2()/(r0*r0);
// Assume en and ed are even integers, and avoid sqrt in the following
int const en2 = en/2;
int const ed2 = ed/2;
cvm::real const xn = cvm::integer_power(l2, en2);
cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
if (calculate_gradients) {
cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
A1.grad += (-1.0)*dFdl2*dl2dx;
A2.grad += dFdl2*dl2dx;
}
return func;
}
#if 0 // AMG: I don't think there's any reason to support anisotropic,
// and I don't have those flags below in calc_value, but
// if I need them, I'll also need to uncomment this method
template<bool calculate_gradients>
cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
int const &en,
int const &ed,
cvm::atom &A1,
cvm::atom &A2)
{
cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
cvm::real const l2 = scal_diff.norm2();
// Assume en and ed are even integers, and avoid sqrt in the following
int const en2 = en/2;
int const ed2 = ed/2;
cvm::real const xn = cvm::integer_power(l2, en2);
cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
if (calculate_gradients) {
cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
(2.0/(r0_vec.y*r0_vec.y))*diff.y,
(2.0/(r0_vec.z*r0_vec.z))*diff.z);
A1.grad += (-1.0)*dFdl2*dl2dx;
A2.grad += dFdl2*dl2dx;
}
return func;
}
#endif
void colvar::groupcoordnum::calc_value()
{
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
// create fake atoms to hold the com coordinates
cvm::atom group1_com_atom;
cvm::atom group2_com_atom;
group1_com_atom.pos = group1->center_of_mass();
group2_com_atom.pos = group2->center_of_mass();
x.real_value = coordnum::switching_function<false>(r0, en, ed,
group1_com_atom, group2_com_atom);
if (b_anisotropic) {
int const flags = coordnum::ef_anisotropic;
x.real_value = coordnum::switching_function<flags>(r0, r0_vec, en, ed,
group1_com_atom,
group2_com_atom,
NULL, 0.0);
} else {
int const flags = coordnum::ef_null;
x.real_value = coordnum::switching_function<flags>(r0, r0_vec, en, ed,
group1_com_atom,
group2_com_atom,
NULL, 0.0);
}
}
void colvar::groupcoordnum::calc_gradients()
{
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
cvm::atom group1_com_atom;
cvm::atom group2_com_atom;
group1_com_atom.pos = group1->center_of_mass();
group2_com_atom.pos = group2->center_of_mass();
coordnum::switching_function<true>(r0, en, ed, group1_com_atom, group2_com_atom);
if (b_anisotropic) {
int const flags = coordnum::ef_gradients | coordnum::ef_anisotropic;
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
group1_com_atom,
group2_com_atom,
NULL, 0.0);
} else {
int const flags = coordnum::ef_gradients;
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
group1_com_atom,
group2_com_atom,
NULL, 0.0);
}
group1->set_weighted_gradient(group1_com_atom.grad);
group2->set_weighted_gradient(group2_com_atom.grad);
}

2
lib/colvars/colvarcomp_distances.cpp
View File

@ -939,7 +939,7 @@ colvar::rmsd::rmsd(std::string const &conf)
bool b_Jacobian_derivative = true;
if (atoms->fitting_group != NULL && b_Jacobian_derivative) {
cvm::log("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: "
cvm::log("The option \"fittingGroup\" (alternative group for fitting) was enabled: "
"Jacobian derivatives of the RMSD will not be calculated.\n");
b_Jacobian_derivative = false;
}

31
lib/colvars/colvarmodule.cpp
View File

@ -139,7 +139,7 @@ int colvarmodule::read_config_file(char const *config_filename)
// read the config file into a string
std::string conf = "";
std::string line;
while (colvarparse::getline_nocomments(config_s, line)) {
while (parse->read_config_line(config_s, line)) {
// Delete lines that contain only white space after removing comments
if (line.find_first_not_of(colvarparse::white_space) != std::string::npos)
conf.append(line+"\n");
@ -159,11 +159,12 @@ int colvarmodule::read_config_string(std::string const &config_str)
// strip the comments away
std::string conf = "";
std::string line;
while (colvarparse::getline_nocomments(config_s, line)) {
while (parse->read_config_line(config_s, line)) {
// Delete lines that contain only white space after removing comments
if (line.find_first_not_of(colvarparse::white_space) != std::string::npos)
conf.append(line+"\n");
}
return parse_config(conf);
}
@ -191,6 +192,12 @@ int colvarmodule::parse_config(std::string &conf)
{
extra_conf.clear();
// Check that the input has matching braces
if (colvarparse::check_braces(conf, 0) != COLVARS_OK) {
return cvm::error("Error: unmatched curly braces in configuration.\n",
INPUT_ERROR);
}
// Parse global options
if (catch_input_errors(parse_global_params(conf))) {
return get_error();
@ -235,6 +242,12 @@ int colvarmodule::parse_config(std::string &conf)
}
std::string const & colvarmodule::get_config() const
{
return parse->get_config();
}
int colvarmodule::append_new_config(std::string const &new_conf)
{
extra_conf += new_conf;
@ -246,9 +259,13 @@ int colvarmodule::parse_global_params(std::string const &conf)
{
colvarmodule *cvm = cvm::main();
std::string index_file_name;
if (parse->get_keyval(conf, "indexFile", index_file_name)) {
cvm->read_index_file(index_file_name.c_str());
{
std::string index_file_name;
size_t pos = 0;
while (parse->key_lookup(conf, "indexFile", &index_file_name, &pos)) {
cvm->read_index_file(index_file_name.c_str());
index_file_name.clear();
}
}
if (parse->get_keyval(conf, "smp", proxy->b_smp_active, proxy->b_smp_active)) {
@ -1073,10 +1090,10 @@ colvarmodule::~colvarmodule()
int colvarmodule::reset()
{
parse->init();
cvm::log("Resetting the Collective Variables module.\n");
parse->init();
// Iterate backwards because we are deleting the elements as we go
for (std::vector<colvarbias *>::reverse_iterator bi = biases.rbegin();
bi != biases.rend();

8
lib/colvars/colvarmodule.h
View File

@ -131,7 +131,7 @@ public:
/// Module-wide error state
/// see constants at the top of this file
protected:
private:
static int errorCode;
@ -274,6 +274,9 @@ public:
/// \brief Parse a "clean" config string (no comments)
int parse_config(std::string &conf);
/// Get the configuration string read so far (includes comments)
std::string const & get_config() const;
// Parse functions (setup internal data based on a string)
/// Allow reading from Windows text files using using std::getline
@ -296,6 +299,9 @@ public:
private:
/// Configuration string read so far by the module (includes comments)
std::string config_string;
/// Auto-generated configuration during parsing (e.g. to implement
/// back-compatibility)
std::string extra_conf;

55
lib/colvars/colvarparse.cpp
View File

@ -43,9 +43,10 @@ template<typename TYPE> bool colvarparse::_get_keyval_scalar_(std::string const
}
} while (b_found);
if (found_count > 1)
cvm::log("Warning: found more than one instance of \""+
std::string(key)+"\".\n");
if (found_count > 1) {
cvm::error("Error: found more than one instance of \""+
std::string(key)+"\".\n", INPUT_ERROR);
}
if (data.size()) {
std::istringstream is(data);
@ -98,9 +99,10 @@ bool colvarparse::_get_keyval_scalar_string_(std::string const &conf,
}
} while (b_found);
if (found_count > 1)
cvm::log("Warning: found more than one instance of \""+
std::string(key)+"\".\n");
if (found_count > 1) {
cvm::error("Error: found more than one instance of \""+
std::string(key)+"\".\n", INPUT_ERROR);
}
if (data.size()) {
std::istringstream is(data);
@ -162,9 +164,10 @@ template<typename TYPE> bool colvarparse::_get_keyval_vector_(std::string const
}
} while (b_found);
if (found_count > 1)
cvm::log("Warning: found more than one instance of \""+
std::string(key)+"\".\n");
if (found_count > 1) {
cvm::error("Error: found more than one instance of \""+
std::string(key)+"\".\n", INPUT_ERROR);
}
if (data.size()) {
std::istringstream is(data);
@ -319,9 +322,10 @@ bool colvarparse::get_keyval(std::string const &conf,
}
} while (b_found);
if (found_count > 1)
cvm::log("Warning: found more than one instance of \""+
std::string(key)+"\".\n");
if (found_count > 1) {
cvm::error("Error: found more than one instance of \""+
std::string(key)+"\".\n", INPUT_ERROR);
}
if (data.size()) {
if ( (data == std::string("on")) ||
@ -535,6 +539,19 @@ int colvarparse::check_keywords(std::string &conf, char const *key)
}
std::istream & colvarparse::read_config_line(std::istream &is,
std::string &line)
{
cvm::getline(is, line);
config_string += line+'\n';
size_t const comment = line.find('#');
if (comment != std::string::npos) {
line.erase(comment);
}
return is;
}
std::istream & colvarparse::getline_nocomments(std::istream &is,
std::string &line)
{
@ -607,7 +624,7 @@ bool colvarparse::key_lookup(std::string const &conf,
}
// check that there are matching braces between here and the end of conf
bool const b_not_within_block = brace_check(conf, pos);
bool const b_not_within_block = (check_braces(conf, pos) == COLVARS_OK);
bool const b_isolated = (b_isolated_left && b_isolated_right &&
b_not_within_block);
@ -781,19 +798,15 @@ std::istream & operator>> (std::istream &is, colvarparse::read_block const &rb)
}
bool colvarparse::brace_check(std::string const &conf,
int colvarparse::check_braces(std::string const &conf,
size_t const start_pos)
{
size_t brace_count = 0;
int brace_count = 0;
size_t brace = start_pos;
while ( (brace = conf.find_first_of("{}", brace)) != std::string::npos) {
while ((brace = conf.find_first_of("{}", brace)) != std::string::npos) {
if (conf[brace] == '{') brace_count++;
if (conf[brace] == '}') brace_count--;
brace++;
}
if (brace_count != 0)
return false;
else
return true;
return (brace_count != 0) ? INPUT_ERROR : COLVARS_OK;
}

24
lib/colvars/colvarparse.h
View File

@ -24,7 +24,7 @@
/// need to parse input inherit from this
class colvarparse {
private:
protected:
/// \brief List of legal keywords for this object: this is updated
/// by each call to colvarparse::get_keyval() or
@ -47,7 +47,7 @@ private:
/// \brief Remove all the values from the config string
void strip_values(std::string &conf);
/// \brief Configuration string of the object
/// \brief Configuration string of the object (includes comments)
std::string config_string;
public:
@ -72,7 +72,7 @@ public:
}
/// Set a new config string for this object
inline void init(const std::string& conf)
inline void init(std::string const &conf)
{
if (! config_string.size()) {
init();
@ -80,7 +80,8 @@ public:
}
}
inline const std::string& get_config()
/// Get the configuration string (includes comments)
inline std::string const & get_config() const
{
return config_string;
}
@ -284,14 +285,19 @@ public:
std::string *data = NULL,
size_t *save_pos = NULL);
/// \brief Reads a configuration line, adds it to config_string, and returns
/// the stream \param is Input stream \param s String that will hold the
/// configuration line, with comments stripped
std::istream & read_config_line(std::istream &is, std::string &line);
/// \brief Works as std::getline() but also removes everything
/// between a comment character and the following newline
static std::istream & getline_nocomments(std::istream &is,
std::string &s);
static std::istream & getline_nocomments(std::istream &is, std::string &s);
/// Check if the content of the file has matching braces
bool brace_check(std::string const &conf,
size_t const start_pos = 0);
/// \brief Check if the content of a config string has matching braces
/// \param conf The configuration string \param start_pos Start the count
/// from this position
static int check_braces(std::string const &conf, size_t const start_pos);
};

20
lib/colvars/colvarproxy.cpp
View File

@ -289,13 +289,23 @@ colvarproxy_smp::colvarproxy_smp()
omp_lock_state = NULL;
#if defined(_OPENMP)
if (smp_thread_id() == 0) {
omp_lock_state = reinterpret_cast<void *>(new omp_lock_t);
omp_init_lock(reinterpret_cast<omp_lock_t *>(omp_lock_state));
}
#endif
}
colvarproxy_smp::~colvarproxy_smp() {}
colvarproxy_smp::~colvarproxy_smp()
{
#if defined(_OPENMP)
if (smp_thread_id() == 0) {
if (omp_lock_state) {
delete reinterpret_cast<omp_lock_t *>(omp_lock_state);
}
}
#endif
}
int colvarproxy_smp::smp_enabled()
@ -499,6 +509,14 @@ char const *colvarproxy_script::script_obj_to_str(unsigned char *obj)
}
std::vector<std::string> colvarproxy_script::script_obj_to_str_vector(unsigned char *obj)
{
cvm::error("Error: trying to print a script object without a scripting "
"language interface.\n", BUG_ERROR);
return std::vector<std::string>();
}
int colvarproxy_script::run_force_callback()
{
return COLVARS_NOT_IMPLEMENTED;

3
lib/colvars/colvarproxy.h
View File

@ -461,6 +461,9 @@ public:
/// Convert a script object (Tcl or Python call argument) to a C string
virtual char const *script_obj_to_str(unsigned char *obj);
/// Convert a script object (Tcl or Python call argument) to a vector of strings
virtual std::vector<std::string> script_obj_to_str_vector(unsigned char *obj);
/// Pointer to the scripting interface object
/// (does not need to be allocated in a new interface)
colvarscript *script;

2
lib/colvars/colvars_version.h
View File

@ -1,5 +1,5 @@
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2018-04-29"
#define COLVARS_VERSION "2018-10-16"
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars

35
lib/colvars/colvarscript.cpp
View File

@ -82,6 +82,14 @@ int colvarscript::run(int objc, unsigned char *const objv[])
int error_code = COLVARS_OK;
// If command is found in map, execute it
std::string const cmd_key("cv_"+cmd);
if (comm_str_map.count(cmd_key) > 0) {
error_code |= (*(comm_fns[comm_str_map[cmd_key]]))(
reinterpret_cast<void *>(this), objc, objv);
return error_code;
}
if (cmd == "colvar") {
if (objc < 3) {
result = "Missing parameters\n" + help_string();
@ -295,11 +303,12 @@ int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[])
}
if (subcmd == "delete") {
size_t i;
for (i = 0; i < cv->biases.size(); i++) {
while (cv->biases.size() > 0) {
size_t i = cv->biases.size()-1;
cvm::log("Warning: before deleting colvar " + cv->name
+ ", deleting parent bias " + cv->biases[i]->name);
delete cv->biases[i];
}
cv->biases.clear();
// colvar destructor is tasked with the cleanup
delete cv;
// TODO this could be done by the destructors
@ -373,6 +382,23 @@ int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[])
return COLVARS_OK;
}
if (subcmd == "modifycvcs") {
if (objc < 4) {
result = "cvcflags: missing parameter: vector of strings";
return COLVARSCRIPT_ERROR;
}
std::vector<std::string> const confs(proxy->script_obj_to_str_vector(objv[3]));
cvm::increase_depth();
int res = cv->update_cvc_config(confs);
cvm::decrease_depth();
if (res != COLVARS_OK) {
result = "Error setting CVC flags";
return COLVARSCRIPT_ERROR;
}
result = "0";
return COLVARS_OK;
}
if ((subcmd == "get") || (subcmd == "set") || (subcmd == "state")) {
return proc_features(cv, objc, objv);
}
@ -547,6 +573,8 @@ std::string colvarscript::help_string() const
Managing the Colvars module:\n\
configfile <file name> -- read configuration from a file\n\
config <string> -- read configuration from the given string\n\
getconfig -- get the module's configuration string\n\
resetindexgroups -- clear the index groups loaded so far\n\
reset -- delete all internal configuration\n\
delete -- delete this Colvars module instance\n\
version -- return version of Colvars code\n\
@ -579,6 +607,7 @@ Accessing collective variables:\n\
colvar <name> gettotalforce -- return total force of colvar <name>\n\
colvar <name> getconfig -- return config string of colvar <name>\n\
colvar <name> cvcflags <fl> -- enable or disable cvcs according to 0/1 flags\n\
colvar <name> modifycvcs <str> -- pass new config strings to each CVC\n\
colvar <name> get <f> -- get the value of the colvar feature <f>\n\
colvar <name> set <f> <val> -- set the value of the colvar feature <f>\n\
\n\

19
lib/colvars/colvarscript.h
View File

@ -71,8 +71,10 @@ public:
cv_help,
cv_version,
cv_config,
cv_getconfig,
cv_configfile,
cv_reset,
cv_resetindexgroups,
cv_delete,
cv_list,
cv_list_biases,
@ -254,6 +256,23 @@ extern "C" {
return COLVARSCRIPT_ERROR;
)
CVSCRIPT(cv_getconfig,
"Get the module's configuration string read so far",
0, 0,
{ },
script->set_str_result(cvm::main()->get_config());
return COLVARS_OK;
)
CVSCRIPT(cv_resetindexgroups,
"Clear the index groups loaded so far, allowing to replace them",
0, 0,
{ },
cvm::main()->index_group_names.clear();
cvm::main()->index_groups.clear();
return COLVARS_OK;
)
CVSCRIPT(cv_addenergy,
"Add an energy to the MD engine",
1, 1,

34
lib/colvars/colvarvalue.cpp
View File

@ -379,6 +379,40 @@ void colvarvalue::set_random()
}
void colvarvalue::set_ones(cvm::real assigned_value)
{
size_t ic;
switch (this->type()) {
case colvarvalue::type_scalar:
this->real_value = assigned_value;
break;
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
this->rvector_value.x = assigned_value;
this->rvector_value.y = assigned_value;
this->rvector_value.z = assigned_value;
break;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
this->quaternion_value.q0 = assigned_value;
this->quaternion_value.q1 = assigned_value;
this->quaternion_value.q2 = assigned_value;
this->quaternion_value.q3 = assigned_value;
break;
case colvarvalue::type_vector:
for (ic = 0; ic < this->vector1d_value.size(); ic++) {
this->vector1d_value[ic] = assigned_value;
}
break;
case colvarvalue::type_notset:
default:
undef_op();
break;
}
}
void colvarvalue::undef_op() const
{
cvm::error("Error: Undefined operation on a colvar of type \""+

43
lib/colvars/colvarvalue.h
View File

@ -187,6 +187,12 @@ public:
return std::sqrt(this->norm2());
}
/// Sum of the components of this colvarvalue (if more than one dimension)
cvm::real sum() const;
/// Return a colvarvalue object of the same type and all components set to 1
colvarvalue ones() const;
/// Square distance between this \link colvarvalue \endlink and another
cvm::real dist2(colvarvalue const &x2) const;
@ -272,17 +278,21 @@ public:
/// Get a single colvarvalue out of elements of the vector
colvarvalue const get_elem(int const i_begin, int const i_end, Type const vt) const;
/// Get a single colvarvalue out of elements of the vector
colvarvalue const get_elem(int const icv) const;
/// Set elements of the vector from a single colvarvalue (uses the rank of x
/// to compute the length)
void set_elem(int const icv, colvarvalue const &x);
/// Set elements of the vector from a single colvarvalue
void set_elem(int const i_begin, int const i_end, colvarvalue const &x);
/// Make each element a random number in N(0,1)
void set_random();
/// Get a single colvarvalue out of elements of the vector
colvarvalue const get_elem(int const icv) const;
/// Set elements of the vector from a single colvarvalue
void set_elem(int const icv, colvarvalue const &x);
/// Make each element equal to the given argument
void set_ones(cvm::real assigned_value = 1.0);
/// Get a scalar number out of an element of the vector
cvm::real operator [] (int const i) const;
@ -683,6 +693,29 @@ inline cvm::real colvarvalue::norm2() const
}
inline cvm::real colvarvalue::sum() const
{
switch (value_type) {
case colvarvalue::type_scalar:
return (this->real_value);
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
return (this->rvector_value).x + (this->rvector_value).y +
(this->rvector_value).z;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
return (this->quaternion_value).q0 + (this->quaternion_value).q1 +
(this->quaternion_value).q2 + (this->quaternion_value).q3;
case colvarvalue::type_vector:
return (this->vector1d_value).sum();
case colvarvalue::type_notset:
default:
return 0.0;
}
}
inline cvm::real colvarvalue::dist2(colvarvalue const &x2) const
{
colvarvalue::check_types(*this, x2);

2
lib/latte/Install.py
View File

@ -78,7 +78,6 @@ def geturl(url,fname):
if which('curl') != None:
cmd = 'curl -L -o "%s" %s' % (fname,url)
print(cmd)
try:
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
success = True
@ -87,7 +86,6 @@ def geturl(url,fname):
if not success and which('wget') != None:
cmd = 'wget -O "%s" %s' % (fname,url)
print(cmd)
try:
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
success = True

4
lib/plumed/.gitignore vendored Normal file
View File

@ -0,0 +1,4 @@
/plumed2*
/includelink
/liblink
/plumed-*

206
lib/plumed/Install.py Normal file
View File

@ -0,0 +1,206 @@
#!/usr/bin/env python
# Install.py tool to download, unpack, build, and link to the plumed2 library
# used to automate the steps described in the README file in this dir
from __future__ import print_function
import sys,os,re,subprocess,hashlib
# help message
help = """
Syntax from src dir: make lib-plumed args="-b"
or: make lib-plumed args="-b -v 2.4.3"
or: make lib-plumed args="-p /usr/local/plumed2-2.4.3"
Syntax from lib dir: python Install.py -b -v 2.4.3
or: python Install.py -b
or: python Install.py -p /usr/local/plumed2-2.4.3
specify one or more options, order does not matter
-b = download and build the plumed2 library
-p = specify folder of existing plumed2 installation
-v = set version of plumed2 to download and build (default: 2.4.3)
Example:
make lib-plumed args="-b" # download/build in lib/plumed/plumed2
make lib-plumed args="-p $HOME/plumed-2.4.3" # use existing Plumed2 installation in $HOME/plumed-2.4.3
"""
# settings
version = "2.4.3"
# known checksums for different PLUMED versions. used to validate the download.
checksums = { \
'2.4.2' : '88188743a6e03ef076e5377d03ebb0e7', \
'2.4.3' : 'b1be7c48971627febc11c61b70767fc5', \
'2.5b' : 'e341bdef469be1da058b8a0b97a3db22', \
}
#checksums = { \
# '2.4.2' : '0f66f24b4c763ae8b2f39574113e9935', \
# '2.4.3' : 'dc38de0ffd59d13950d8f1ef1ce05574', \
# }
# print error message or help
def error(str=None):
if not str: print(help)
else: print("ERROR",str)
sys.exit()
# expand to full path name
# process leading '~' or relative path
def fullpath(path):
return os.path.abspath(os.path.expanduser(path))
def which(program):
def is_exe(fpath):
return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
fpath, fname = os.path.split(program)
if fpath:
if is_exe(program):
return program
else:
for path in os.environ["PATH"].split(os.pathsep):
path = path.strip('"')
exe_file = os.path.join(path, program)
if is_exe(exe_file):
return exe_file
return None
def geturl(url,fname):
success = False
if which('curl') != None:
cmd = 'curl -L -o "%s" %s' % (fname,url)
try:
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
success = True
except subprocess.CalledProcessError as e:
print("Calling curl failed with: %s" % e.output.decode('UTF-8'))
if not success and which('wget') != None:
cmd = 'wget -O "%s" %s' % (fname,url)
try:
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
success = True
except subprocess.CalledProcessError as e:
print("Calling wget failed with: %s" % e.output.decode('UTF-8'))
if not success:
error("Failed to download source code with 'curl' or 'wget'")
return
def checkmd5sum(md5sum,fname):
with open(fname,'rb') as fh:
m = hashlib.md5()
while True:
data = fh.read(81920)
if not data:
break
m.update(data)
fh.close()
return m.hexdigest() == md5sum
# parse args
args = sys.argv[1:]
nargs = len(args)
if nargs == 0: error()
homepath = "."
buildflag = False
pathflag = False
suffixflag = False
linkflag = True
iarg = 0
while iarg < nargs:
if args[iarg] == "-v":
if iarg+2 > nargs: error()
version = args[iarg+1]
iarg += 2
elif args[iarg] == "-p":
if iarg+2 > nargs: error()
plumedpath = fullpath(args[iarg+1])
pathflag = True
iarg += 2
elif args[iarg] == "-b":
buildflag = True
iarg += 1
else: error()
homepath = fullpath(homepath)
if (pathflag):
if not os.path.isdir(plumedpath): error("Plumed2 path does not exist")
homedir = plumedpath
if (buildflag and pathflag):
error("Cannot use -b and -p flag at the same time")
if (not buildflag and not pathflag):
error("Have to use either -b or -p flag")
# download and unpack plumed2 tarball
if buildflag:
url = "https://github.com/plumed/plumed2/releases/download/v%s/plumed-src-%s.tgz" % (version,version)
filename = "plumed-src-%s.tar.gz" %version
#url = "https://github.com/plumed/plumed2/archive/v%s.tar.gz" % version
#filename = "v%s.tar.gz" %version
print("Downloading plumed ...")
geturl(url,filename)
# verify downloaded archive integrity via md5 checksum, if known.
if version in checksums:
if not checkmd5sum(checksums[version],filename):
error("Checksum for plumed2 library does not match")
print("Unpacking plumed2 source tarball ...")
if os.path.exists("%s/plumed-%s" % (homepath,version)):
cmd = 'rm -rf "%s/plumed-%s"' % (homepath,version)
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
#if os.path.exists("%s/plumed2-%s" % (homepath,version)):
# cmd = 'rm -rf "%s/plumed2-%s"' % (homepath,version)
# subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
if os.path.exists("%s/plumed2" % (homepath)):
cmd = 'rm -rf "%s/plumed2"' % (homepath)
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
cmd = 'cd "%s"; tar -xzvf %s' % (homepath,filename)
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
os.remove("%s/%s" % (homepath,filename))
# build plumed
if buildflag:
print("Building plumed ...")
cmd = 'cd %s/plumed-%s; ./configure --prefix=%s/plumed2 --enable-static-patch ; make ; make install' % (homepath,version,homepath)
#cmd = 'cd %s/plumed2-%s; ./configure --prefix=%s/plumed2 --enable-static-patch ; make ; make install' % (homepath,version,homepath)
txt = subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
print(txt.decode('UTF-8'))
#
# create 2 links in lib/plumed to plumed2 installation dir
if linkflag:
print("Creating links to plumed2 include and lib files")
if os.path.isfile("includelink") or os.path.islink("includelink"):
os.remove("includelink")
if os.path.isfile("liblink") or os.path.islink("liblink"):
os.remove("liblink")
cmd = 'ln -s "%s/plumed2/include" includelink' % homepath
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
cmd = 'ln -s "%s/plumed2/lib" liblink' % homepath
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
if os.path.isfile("Makefile.lammps.static"):
print("Creating Makefile.lammps")
cmd = 'cat liblink/plumed/src/lib/Plumed.inc.static Makefile.lammps.static > Makefile.lammps'
subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)

5
lib/plumed/Makefile.lammps.static Normal file
View File

@ -0,0 +1,5 @@
# Settings that the LAMMPS build will import when this package library is used
plumed_SYSINC =
plumed_SYSLIB = $(PLUMED_LOAD)
plumed_SYSPATH =

56
lib/plumed/README Normal file
View File

@ -0,0 +1,56 @@
This directory contains links to the PLUMED library which is required
to use the PLUMED package and its fix plumed command in a
LAMMPS input script. PLUMED should only be downloaded into this directory if
you wish to statically link the library. If you wish to link PLUMED as
a dynamic library (as we recommend) then you can compile and build PLUMED
separately to LAMMPS. To use PLUMED in conjuction with LAMMPS you then simply
need to ensure that the PLUMED library is in your path at runtime.
More info about the PLUMED library can be found at http://www.plumed.org.
You can type "make lib-plumed" from the src directory to see help on
how to download, build and statically link PLUMED via make commands, or you can
do the same thing by typing "python Install.py" from within this
directory. Alternatively you can download and build PLUMED manually by following the instructions
below.
-----------------
Instructions:
1. Download PLUMED either as a tarball from
http://www.plumed.org/get-it
or clone it using git clone https://github.com/plumed/plumed2.git.
If you download the tarball
unpack it in unpack it in this /lib/plumed directory.
Similarly if you clone it clone it to the /lib/plumed
directory.
2. Compile PLUMED from within its home directory. In the
simplest cases this be done by issuing the commands
% ./configure
% make
More detailed instructions can be found at
http://plumed.github.io/doc-master/user-doc/html/_installation.html
3. There is no need to install PLUMED if you only wish
to use it from LAMMPS. You should thus only run
make install if you want to use PLUMED as a stand-alone
code or from some other code. To install it you can
run the following commands:
a) install under the default /usr/local
% sudo make install
b) install under a user-writeable location by first
changing the PREFIX variable when running the
configure command file, then
% make install
-----------------
When these steps are complete you can build LAMMPS
with the PLUMED package installed:
% cd lammps/src
% make yes-user-plumed
% make mpi (or whatever target you wish)

2
src/.gitignore vendored
View File

@ -59,6 +59,8 @@
/colvarproxy_lammps_version.h
/fix_colvars.cpp
/fix_colvars.h
/fix_plumed.cpp
/fix_plumed.h
/dump_molfile.cpp
/dump_molfile.h
/molfile_interface.cpp

9
src/Makefile
View File

@ -61,21 +61,22 @@ PACKUSER = user-atc user-awpmd user-bocs user-cgdna user-cgsdk user-colvars \
user-diffraction user-dpd user-drude user-eff user-fep user-h5md \
user-intel user-lb user-manifold user-meamc user-meso \
user-mgpt user-misc user-mofff user-molfile \
user-netcdf user-omp user-phonon user-ptm user-qmmm user-qtb \
user-quip user-reaxc user-scafacos user-smd user-smtbq \
user-netcdf user-omp user-phonon user-plumed user-ptm user-qmmm \
user-qtb user-quip user-reaxc user-scafacos user-smd user-smtbq \
user-sdpd user-sph user-tally user-uef user-vtk
PACKLIB = compress gpu kim kokkos latte meam message mpiio mscg poems \
python reax voronoi \
user-atc user-awpmd user-colvars user-h5md user-lb user-molfile \
user-netcdf user-qmmm user-quip user-scafacos user-smd user-vtk
user-netcdf user-plumed user-qmmm user-quip user-scafacos \
user-smd user-vtk
PACKSYS = compress mpiio python user-lb
PACKINT = gpu kokkos meam message poems reax user-atc user-awpmd user-colvars
PACKEXT = kim latte mscg voronoi \
user-h5md user-molfile user-netcdf user-qmmm user-quip \
user-h5md user-molfile user-netcdf user-plumed user-qmmm user-quip \
user-smd user-vtk
PACKALL = $(PACKAGE) $(PACKUSER)

3
src/SPIN/fix_nve_spin.cpp
View File

@ -60,8 +60,11 @@ static const char cite_fix_nve_spin[] =
"dynamics and molecular dynamics},\n"
"author={Tranchida, J and Plimpton, SJ and Thibaudeau, P and Thompson, AP},\n"
"journal={Journal of Computational Physics},\n"
"volume={372},\n"
"pages={406-425},\n"
"year={2018},\n"
"publisher={Elsevier}\n"
"doi={10.1016/j.jcp.2018.06.042}\n"
"}\n\n";
enum{NONE};

2
src/USER-COLVARS/colvarproxy_lammps_version.h
View File

@ -1,5 +1,5 @@
#ifndef COLVARPROXY_VERSION
#define COLVARPROXY_VERSION "2018-04-29"
#define COLVARPROXY_VERSION "2018-08-29"
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars

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

@ -78,6 +78,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
fix1 = NULL;
fix2 = NULL;
fix3 = NULL;
if (narg < 8) error->all(FLERR,"Illegal fix bond/react command: "
"too few arguments");
@ -163,6 +164,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(seed,nreacts,"bond/react:seed");
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(iatomtype,nreacts,"bond/react:iatomtype");
memory->create(jatomtype,nreacts,"bond/react:jatomtype");
memory->create(ibonding,nreacts,"bond/react:ibonding");
@ -178,6 +180,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
fraction[i] = 1;
seed[i] = 12345;
stabilize_steps_flag[i] = 0;
update_edges_flag[i] = 0;
// set default limit duration to 60 timesteps
limit_duration[i] = 60;
reaction_count[i] = 0;
@ -249,6 +252,14 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
limit_duration[rxn] = force->numeric(FLERR,arg[iarg+1]);
stabilize_steps_flag[rxn] = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"update_edges") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix bond/react command: "
"'update_edges' has too few arguments");
if (strcmp(arg[iarg+1],"none") == 0) update_edges_flag[rxn] = 0;
else if (strcmp(arg[iarg+1],"charges") == 0) update_edges_flag[rxn] = 1;
else if (strcmp(arg[iarg+1],"custom") == 0) update_edges_flag[rxn] = 2;
else error->all(FLERR,"Illegal value for 'update_edges' keyword'");
iarg += 2;
} else error->all(FLERR,"Illegal fix bond/react command: unknown keyword");
}
}
@ -263,15 +274,23 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(reverse_equiv,max_natoms,2,nreacts,"bond/react:reverse_equiv");
memory->create(edge,max_natoms,nreacts,"bond/react:edge");
memory->create(landlocked_atoms,max_natoms,nreacts,"bond/react:landlocked_atoms");
memory->create(custom_edges,max_natoms,nreacts,"bond/react:custom_edges");
for (int j = 0; j < nreacts; j++)
for (int i = 0; i < max_natoms; i++) edge[i][j] = 0;
for (int i = 0; i < max_natoms; i++) {
edge[i][j] = 0;
if (update_edges_flag[j] == 1) custom_edges[i][j] = 1;
else custom_edges[i][j] = 0;
}
// read all superimpose files afterward
// read all map files afterward
for (int i = 0; i < nreacts; i++) {
open(files[i]);
onemol = atom->molecules[unreacted_mol[i]];
twomol = atom->molecules[reacted_mol[i]];
if (onemol->natoms != twomol->natoms)
error->all(FLERR,"Post-reacted template must contain the same "
"number of atoms as the pre-reacted template");
get_molxspecials();
read(i);
fclose(fp);
@ -347,6 +366,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
id_fix1 = NULL;
id_fix2 = NULL;
id_fix3 = NULL;
statted_id = NULL;
custom_exclude_flag = 0;
}
/* ---------------------------------------------------------------------- */
@ -370,6 +392,7 @@ FixBondReact::~FixBondReact()
memory->destroy(edge);
memory->destroy(equivalences);
memory->destroy(reverse_equiv);
memory->destroy(custom_edges);
memory->destroy(nevery);
memory->destroy(cutsq);
@ -379,6 +402,7 @@ FixBondReact::~FixBondReact()
memory->destroy(seed);
memory->destroy(limit_duration);
memory->destroy(stabilize_steps_flag);
memory->destroy(update_edges_flag);
memory->destroy(iatomtype);
memory->destroy(jatomtype);
@ -418,11 +442,15 @@ FixBondReact::~FixBondReact()
// check nfix in case all fixes have already been deleted
if (id_fix1 == NULL && modify->nfix) modify->delete_fix(id_fix1);
delete [] id_fix1;
if (id_fix3 == NULL && modify->nfix) modify->delete_fix(id_fix3);
delete [] id_fix3;
}
if (id_fix2 == NULL && modify->nfix) modify->delete_fix(id_fix2);
delete [] id_fix2;
delete [] statted_id;
delete [] guess_branch;
delete [] pioneer_count;
}
@ -453,61 +481,126 @@ void FixBondReact::post_constructor()
int ifix = modify->find_fix(id_fix2);
if (ifix == -1) {
char **newarg = new char*[8];
char **newarg = new char*[7];
newarg[0] = (char *) "bond_react_props_internal";
newarg[1] = (char *) "all"; // group ID is ignored
newarg[2] = (char *) "property/atom";
newarg[3] = (char *) "i_limit_tags";
newarg[4] = (char *) "i_statted_tags";
newarg[5] = (char *) "i_react_tags";
newarg[6] = (char *) "ghost";
newarg[7] = (char *) "yes";
modify->add_fix(8,newarg);
fix2 = modify->fix[modify->nfix-1];
newarg[4] = (char *) "i_react_tags";
newarg[5] = (char *) "ghost";
newarg[6] = (char *) "yes";
modify->add_fix(7,newarg);
delete [] newarg;
}
// create master_group if not already existing
if (group->find(master_group) == -1) {
group->find_or_create(master_group);
char **newarg;
newarg = new char*[5];
newarg[0] = master_group;
newarg[1] = (char *) "dynamic";
newarg[2] = (char *) "all";
newarg[3] = (char *) "property";
newarg[4] = (char *) "limit_tags";
group->assign(5,newarg);
delete [] newarg;
}
// on to statted_tags (system-wide thermostat)
// intialize per-atom statted_flags to 1
// (only if not already initialized by restart)
// NOTE: limit_tags and react_tags automaticaly intitialized to zero (unless read from restart)
if (fix2->restart_reset != 1) {
int flag;
int index = atom->find_custom("statted_tags",flag);
int *i_statted_tags = atom->ivector[index];
for (int i = 0; i < atom->nlocal; i++)
i_statted_tags[i] = 1;
}
group->find_or_create(master_group);
char **newarg;
newarg = new char*[5];
newarg[0] = master_group;
newarg[1] = (char *) "dynamic";
newarg[2] = (char *) "all";
newarg[3] = (char *) "property";
newarg[4] = (char *) "limit_tags";
group->assign(5,newarg);
delete [] newarg;
if (stabilization_flag == 1) {
// create exclude_group if not already existing
if (group->find(exclude_group) == -1) {
int igroup = group->find(exclude_group);
// create exclude_group if not already existing, or use as parent group if static
if (igroup == -1 || group->dynamic[igroup] == 0) {
// create stabilization per-atom property
len = strlen("bond_react_stabilization_internal") + 1;
id_fix3 = new char[len];
strcpy(id_fix3,"bond_react_stabilization_internal");
ifix = modify->find_fix(id_fix3);
if (ifix == -1) {
char **newarg = new char*[6];
newarg[0] = (char *) id_fix3;
newarg[1] = (char *) "all"; // group ID is ignored
newarg[2] = (char *) "property/atom";
newarg[3] = (char *) "i_statted_tags";
newarg[4] = (char *) "ghost";
newarg[5] = (char *) "yes";
modify->add_fix(6,newarg);
fix3 = modify->fix[modify->nfix-1];
delete [] newarg;
}
len = strlen("statted_tags") + 1;
statted_id = new char[len];
strcpy(statted_id,"statted_tags");
// if static group exists, use as parent group
// also, rename dynamic exclude_group by appending '_REACT'
char *exclude_PARENT_group;
int n = strlen(exclude_group) + 1;
exclude_PARENT_group = new char[n];
strcpy(exclude_PARENT_group,exclude_group);
n += strlen("_REACT");
delete [] exclude_group;
exclude_group = new char[n];
strcpy(exclude_group,exclude_PARENT_group);
strcat(exclude_group,"_REACT");
group->find_or_create(exclude_group);
char **newarg;
newarg = new char*[5];
newarg[0] = exclude_group;
newarg[1] = (char *) "dynamic";
newarg[2] = (char *) "all";
if (igroup == -1) newarg[2] = (char *) "all";
else newarg[2] = (char *) exclude_PARENT_group;
newarg[3] = (char *) "property";
newarg[4] = (char *) "statted_tags";
group->assign(5,newarg);
delete [] newarg;
}
delete [] exclude_PARENT_group;
// on to statted_tags (system-wide thermostat)
// intialize per-atom statted_flags to 1
// (only if not already initialized by restart)
if (fix3->restart_reset != 1) {
int flag;
int index = atom->find_custom("statted_tags",flag);
int *i_statted_tags = atom->ivector[index];
for (int i = 0; i < atom->nlocal; i++)
i_statted_tags[i] = 1;
}
} else {
// sleeping code, for future capabilities
custom_exclude_flag = 1;
// first we have to find correct fix group reference
int n = strlen("GROUP_") + strlen(exclude_group) + 1;
char *fix_group = new char[n];
strcpy(fix_group,"GROUP_");
strcat(fix_group,exclude_group);
int ifix = modify->find_fix(fix_group);
Fix *fix = modify->fix[ifix];
delete [] fix_group;
// this returns names of corresponding property
int unused;
char * idprop;
idprop = (char *) fix->extract("property",unused);
if (idprop == NULL)
error->all(FLERR,"Exclude group must be a per-atom property group");
len = strlen(idprop) + 1;
statted_id = new char[len];
strcpy(statted_id,idprop);
// intialize per-atom statted_tags to 1
// need to correct for smooth restarts
//int flag;
//int index = atom->find_custom(statted_id,flag);
//int *i_statted_tags = atom->ivector[index];
//for (int i = 0; i < atom->nlocal; i++)
// i_statted_tags[i] = 1;
}
// let's create a new nve/limit fix to limit newly reacted atoms
len = strlen("bond_react_MASTER_nve_limit") + 1;
@ -526,40 +619,6 @@ void FixBondReact::post_constructor()
fix1 = modify->fix[modify->nfix-1];
delete [] newarg;
}
}
// currently must redefine dynamic groups so they are updated at proper time
// -> should double check as to why
int must_redefine_groups = 1;
if (must_redefine_groups) {
group->find_or_create(master_group);
char **newarg;
newarg = new char*[5];
newarg[0] = master_group;
newarg[1] = (char *) "dynamic";
newarg[2] = (char *) "all";
newarg[3] = (char *) "property";
newarg[4] = (char *) "limit_tags";
group->assign(5,newarg);
delete [] newarg;
}
if (stabilization_flag == 1) {
if (must_redefine_groups) {
group->find_or_create(exclude_group);
char **newarg;
newarg = new char*[5];
newarg[0] = exclude_group;
newarg[1] = (char *) "dynamic";
newarg[2] = (char *) "all";
newarg[3] = (char *) "property";
newarg[4] = (char *) "statted_tags";
group->assign(5,newarg);
delete [] newarg;
}
}
}
@ -1804,8 +1863,11 @@ void FixBondReact::limit_bond(int limit_bond_mode)
int index1 = atom->find_custom("limit_tags",flag);
int *i_limit_tags = atom->ivector[index1];
int index2 = atom->find_custom("statted_tags",flag);
int *i_statted_tags = atom->ivector[index2];
int *i_statted_tags;
if (stabilization_flag == 1) {
int index2 = atom->find_custom(statted_id,flag);
i_statted_tags = atom->ivector[index2];
}
int index3 = atom->find_custom("react_tags",flag);
int *i_react_tags = atom->ivector[index3];
@ -1813,7 +1875,7 @@ void FixBondReact::limit_bond(int limit_bond_mode)
for (int i = 0; i < temp_limit_num; i++) {
// update->ntimestep could be 0. so add 1 throughout
i_limit_tags[atom->map(temp_limit_glove[i])] = update->ntimestep + 1;
i_statted_tags[atom->map(temp_limit_glove[i])] = 0;
if (stabilization_flag == 1) i_statted_tags[atom->map(temp_limit_glove[i])] = 0;
i_react_tags[atom->map(temp_limit_glove[i])] = rxnID;
}
@ -1834,8 +1896,11 @@ void FixBondReact::unlimit_bond()
int index1 = atom->find_custom("limit_tags",flag);
int *i_limit_tags = atom->ivector[index1];
int index2 = atom->find_custom("statted_tags",flag);
int *i_statted_tags = atom->ivector[index2];
int *i_statted_tags;
if (stabilization_flag == 1) {
int index2 = atom->find_custom(statted_id,flag);
i_statted_tags = atom->ivector[index2];
}
int index3 = atom->find_custom("react_tags",flag);
int *i_react_tags = atom->ivector[index3];
@ -1845,7 +1910,7 @@ void FixBondReact::unlimit_bond()
// first '1': indexing offset, second '1': for next step
if (i_limit_tags[i] != 0 && (update->ntimestep + 1 - i_limit_tags[i]) > limit_duration[i_react_tags[i]]) { // + 1
i_limit_tags[i] = 0;
i_statted_tags[i] = 1;
if (stabilization_flag == 1) i_statted_tags[i] = 1;
i_react_tags[i] = 0;
}
}
@ -2022,13 +2087,13 @@ void FixBondReact::update_everything()
}
// update charges and types of landlocked atoms
// here, add check for charge instead of requiring it
for (int i = 0; i < update_num_mega; i++) {
rxnID = update_mega_glove[0][i];
twomol = atom->molecules[reacted_mol[rxnID]];
for (int j = 0; j < twomol->natoms; j++) {
int jj = equivalences[j][1][rxnID]-1;
if (landlocked_atoms[j][rxnID] == 1 && atom->map(update_mega_glove[jj+1][i]) >= 0 &&
if ((landlocked_atoms[j][rxnID] == 1 || custom_edges[jj][rxnID] == 1) &&
atom->map(update_mega_glove[jj+1][i]) >= 0 &&
atom->map(update_mega_glove[jj+1][i]) < nlocal) {
type[atom->map(update_mega_glove[jj+1][i])] = twomol->type[j];
if (twomol->qflag && atom->q_flag) {
@ -2470,6 +2535,7 @@ void FixBondReact::read(int myrxn)
if (strstr(line,"edgeIDs")) sscanf(line,"%d",&nedge);
else if (strstr(line,"equivalences")) sscanf(line,"%d",&nequivalent);
else if (strstr(line,"customIDs")) sscanf(line,"%d",&ncustom);
else break;
}
@ -2482,7 +2548,7 @@ void FixBondReact::read(int myrxn)
// loop over sections of superimpose file
int equivflag = 0, edgeflag = 0, bondflag = 0;
int equivflag = 0, edgeflag = 0, bondflag = 0, customedgesflag = 0;
while (strlen(keyword)) {
if (strcmp(keyword,"BondingIDs") == 0) {
bondflag = 1;
@ -2496,6 +2562,9 @@ void FixBondReact::read(int myrxn)
} else if (strcmp(keyword,"Equivalences") == 0) {
equivflag = 1;
Equivalences(line, myrxn);
} else if (strcmp(keyword,"Custom Edges") == 0) {
customedgesflag = 1;
CustomEdges(line, myrxn);
} else error->one(FLERR,"Unknown section in superimpose file");
parse_keyword(1,line,keyword);
@ -2505,6 +2574,12 @@ void FixBondReact::read(int myrxn)
// error check
if (bondflag == 0 || equivflag == 0)
error->all(FLERR,"Superimpose file missing BondingIDs or Equivalences section\n");
if (update_edges_flag[myrxn] == 2 && customedgesflag == 0)
error->all(FLERR,"Map file must have a Custom Edges section when using 'update_edges custom'\n");
if (update_edges_flag[myrxn] != 2 && customedgesflag == 1)
error->all(FLERR,"Specify 'update_edges custom' to include Custom Edges section in map file\n");
}
void FixBondReact::EdgeIDs(char *line, int myrxn)
@ -2535,6 +2610,26 @@ void FixBondReact::Equivalences(char *line, int myrxn)
}
}
void FixBondReact::CustomEdges(char *line, int myrxn)
{
// 0 for 'none', 1 for 'charges'
int tmp;
int n = MAX(strlen("none"),strlen("charges")) + 1;
char *edgemode = new char[n];
for (int i = 0; i < ncustom; i++) {
readline(line);
sscanf(line,"%d %s",&tmp,edgemode);
if (strcmp(edgemode,"none") == 0)
custom_edges[tmp-1][myrxn] = 0;
else if (strcmp(edgemode,"charges") == 0)
custom_edges[tmp-1][myrxn] = 1;
else
error->one(FLERR,"Illegal value in 'Custom Edges' section of map file");
}
delete [] edgemode;
}
void FixBondReact::open(char *file)
{
fp = fopen(file,"r");

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

@ -57,7 +57,9 @@ class FixBondReact : public Fix {
double **cutsq,*fraction;
tagint lastcheck;
int stabilization_flag;
int custom_exclude_flag;
int *stabilize_steps_flag;
int *update_edges_flag;
int status;
int *groupbits;
@ -76,9 +78,9 @@ class FixBondReact : public Fix {
class Molecule *onemol; // pre-reacted molecule template
class Molecule *twomol; // post-reacted molecule template
Fix *fix1; // nve/limit used to relax reaction sites
Fix *fix2; // properties/atom used to indicate 1) indicate relaxing atoms
// 2) system-wide thermostat
// 3) to which 'react' atom belongs
Fix *fix2; // properties/atom used to indicate 1) relaxing atoms
// 2) to which 'react' atom belongs
Fix *fix3; // property/atom used for system-wide thermostat
class RanMars **random;
class NeighList *list;
@ -87,6 +89,8 @@ class FixBondReact : public Fix {
char *nve_limit_xmax; // indicates max distance allowed to move when relaxing
char *id_fix1; // id of internally created fix nve/limit
char *id_fix2; // id of internally created fix per-atom properties
char *id_fix3; // id of internally created 'stabilization group' per-atom property fix
char *statted_id; // name of 'stabilization group' per-atom property
char *master_group; // group containing relaxing atoms from all fix rxns
char *exclude_group; // group for system-wide thermostat
@ -97,7 +101,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; // number of edge, equivalent atoms in mapping file
int nedge,nequivalent,ncustom; // number of edge, equivalent, custom atoms in mapping file
int attempted_rxn; // there was an attempt!
int *local_rxn_count;
int *ghostly_rxn_count;
@ -111,6 +115,7 @@ class FixBondReact : public Fix {
int ***equivalences; // relation between pre- and post-reacted templates
int ***reverse_equiv; // re-ordered equivalences
int **landlocked_atoms; // all atoms at least three bonds away from edge atoms
int **custom_edges; // atoms in molecule templates with incorrect valences
int **nxspecial,**onemol_nxspecial,**twomol_nxspecial; // full number of 1-4 neighbors
tagint **xspecial,**onemol_xspecial,**twomol_xspecial; // full 1-4 neighbor list
@ -132,6 +137,7 @@ class FixBondReact : public Fix {
void read(int);
void EdgeIDs(char *,int);
void Equivalences(char *,int);
void CustomEdges(char *,int);
void make_a_guess ();
void neighbor_loop();

65
src/USER-PLUMED/Install.sh Executable file
View File

@ -0,0 +1,65 @@
# Install/unInstall package files in LAMMPS
# mode = 0/1/2 for uninstall/install/update
mode=$1
# enforce using portable C locale
LC_ALL=C
export LC_ALL
# arg1 = file, arg2 = file it depends on
action () {
if (test $mode = 0) then
rm -f ../$1
elif (! cmp -s $1 ../$1) then
if (test -z "$2" || test -e ../$2) then
cp $1 ..
if (test $mode = 2) then
echo " updating src/$1"
fi
fi
elif (test -n "$2") then
if (test ! -e ../$2) then
rm -f ../$1
fi
fi
}
# all package files with no dependencies
for file in *.cpp *.h; do
test -f ${file} && action $file
done
# edit 2 Makefile.package files to include/exclude package info
if (test $1 = 1) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*plumed[^ \t]* //' ../Makefile.package
sed -i -e 's|^PKG_INC =[ \t]*|&-I../../lib/plumed/includelink |' ../Makefile.package
sed -i -e 's|^PKG_SYSINC =[ \t]*|&$(plumed_SYSINC) |' ../Makefile.package
sed -i -e 's|^PKG_SYSLIB =[ \t]*|&$(plumed_SYSLIB) |' ../Makefile.package
sed -i -e 's|^PKG_SYSPATH =[ \t]*|&$(plumed_SYSPATH) |' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*plumed.*$/d' ../Makefile.package.settings
# multiline form needed for BSD sed on Macs
sed -i -e '4 i \
include ..\/..\/lib\/plumed\/Makefile.lammps
' ../Makefile.package.settings
fi
elif (test $1 = 0) then
if (test -e ../Makefile.package) then
sed -i -e 's/[^ \t]*plumed[^ \t]* //' ../Makefile.package
fi
if (test -e ../Makefile.package.settings) then
sed -i -e '/^include.*plumed.*$/d' ../Makefile.package.settings
fi
fi

70
src/USER-PLUMED/README Normal file
View File

@ -0,0 +1,70 @@
This package implements the "fix plumed" command, which can be used
in a LAMMPS input script.
The fix allows enhanced sampling methods such as umbrella sampling and
metadynamics to be used. Furthermore, PLUMED can be used to perform a
wide range of analyses on trajectories on the fly as they are generated.
The package uses the "PLUMED" library, whose source code is not included
in the LAMMPS source code distribution. The files in the USER-PLUMED package
folder implement an interface between LAMMPS and PLUMED, that are written
and maintained by Gareth Tribello (gareth.tribello@gmail.com).
PLUMED must instead be downloaded and compiled separately to LAMMPS. This building
and compiling of PLUMED can be done before or after the building of LAMMPS as LAMMPS
can call PLUMED as a dynamic library. There is also the possibility to link PLUEMD
statically. In this case a copy of PLUMED must be downloaded into the lib/plumed
directory. This copy of PLUMED will then always be linked into the code at compile
time.
However you decide to link PLUMED (statically or dynamically) you must run the command:
make yes-user-plumed
before compiling LAMMPS in order to enable the module. In addition, if you have chosen to
link PLUMED dynamically you must ensure that PLUMED is in your
PATH when running a LAMMPS calculation that takes advantage of PLUMED. If
PLUMED is linked as a runtime library and if PLUMED is not in the PATH an error will be returned whenever LAMMPS encounters
the fix plumed command in its input. To be clear, however, a LAMMPS executable that was dynamically linked with PLUMED will run
even if PLUMED is not in the path if as long as the input does not contain a fix
plumed command.
If you wish to statically link PLUMED you must download PLUMED to the /lib/plumed directory before compiling LAMMPS. You can
download a tar ball into that directory or you can clone the plumed2 repository from github there. Once you have created a
directory containing a distribution of PLUMED within /lib/plumed you then must build PLUMED within that directory by issuing
the usual commands. It is worth noting that we have provided a script that will download and build PLUMED for you with
a minimal set of options. To run this script you need to issue the following command:
make lib-plumed args="-b"
in the src directory.
More info about the PLUMED library can be found at:
www.plumed.org
and in the reference articles:
PLUMED2: New feathers for an old bird
G.A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi,
Comp. Phys. Comm 185, 604 (2014)
https://doi.org/10.1016/j.cpc.2013.09.018
PLUMED: a portable plugin for free energy calculations with molecular dynamics
M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi, P. Raiteri, D. Donadio, F. Marinelli, F. Pietrucci, R.A. Broglia and M. Parrinello
Comp. Phys. Comm. 180, 1961 (2009)
https://doi.org/10.1016/j.cpc.2009.05.011
Instructions explaining how to use PLUMED and LAMMPS in tandem can be found on the PLUMED website, which also gives
numerous example scripts for PLUMED as well as citations to articles that dcoment the various methods that are
implemented within PLUMED.
There are also example scripts for using this package in the folder
examples/USER/plumed, as well as on the GitHub page for PLUMED.
Please contact Gareth Tribello (gareth.tribello@gmail.com) for questions
regarding this package.
---------------------------------
Version: 2016-12-22

559
src/USER-PLUMED/fix_plumed.cpp Normal file
View File

@ -0,0 +1,559 @@
/* ----------------------------------------------------------------------
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: Gareth Tribello (Queens U, Belfast)
Pablo Piaggi (EPFL)
------------------------------------------------------------------------- */
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "atom.h"
#include "comm.h"
#include "update.h"
#include "force.h"
#include "respa.h"
#include "domain.h"
#include "error.h"
#include "group.h"
#include "fix_plumed.h"
#include "universe.h"
#include "compute.h"
#include "modify.h"
#include "pair.h"
/*
Use statically linked C++ interface to plumed
*/
#define __PLUMED_WRAPPER_CXX 1
#include "plumed/wrapper/Plumed.h"
/* -------------------------------------------------------------------- */
using namespace LAMMPS_NS;
using namespace FixConst;
#define INVOKED_SCALAR 1
FixPlumed::FixPlumed(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
p(NULL), nlocal(0), gatindex(NULL), masses(NULL), charges(NULL),
id_pe(NULL), id_press(NULL)
{
if (!atom->tag_enable)
error->all(FLERR,"Fix plumed requires atom tags");
if (atom->tag_consecutive() == 0)
error->all(FLERR,"Fix plumed requires consecutive atom IDs");
if (igroup != 0 && comm->me == 0)
error->warning(FLERR,"Fix group for fix plumed is not 'all'. "
"Group will be ignored.");
p=new PLMD::Plumed;
// Check API version
int api_version;
p->cmd("getApiVersion",&api_version);
if (api_version > 6)
error->all(FLERR,"Incompatible API version for PLUMED in fix plumed");
// If the -partition option is activated then enable
// inter-partition communication
if (universe->existflag == 1) {
int me;
MPI_Comm inter_comm;
MPI_Comm_rank(world,&me);
// Change MPI_COMM_WORLD to universe->uworld which seems more appropriate
MPI_Comm_split(universe->uworld,me,0,&inter_comm);
p->cmd("GREX setMPIIntracomm",&world);
if (me == 0) {
// The inter-partition communicator is only defined for the root in
// each partition (a.k.a. world). This is due to the way in which
// it is defined inside plumed.
p->cmd("GREX setMPIIntercomm",&inter_comm);
}
p->cmd("GREX init",NULL);
}
// The general communicator is independent of the existence of partitions,
// if there are partitions, world is defined within each partition,
// whereas if partitions are not defined then world is equal to
// MPI_COMM_WORLD.
p->cmd("setMPIComm",&world);
// Set up units
// LAMMPS units wrt kj/mol - nm - ps
// Set up units
if(strcmp(update->unit_style,"lj") == 0) {
// LAMMPS units lj
p->cmd("setNaturalUnits");
} else {
// Conversion factor from LAMMPS energy units to kJ/mol (units of PLUMED)
double energyUnits=1.0;
// LAMMPS units real :: kcal/mol;
if (strcmp(update->unit_style,"real") == 0) {
energyUnits=4.184;
// LAMMPS units metal :: eV;
} else if (strcmp(update->unit_style,"metal") == 0) {
energyUnits=96.48530749925792;
// LAMMPS units si :: Joule;
} else if (strcmp(update->unit_style,"si") == 0) {
energyUnits=0.001;
// LAMMPS units cgs :: erg;
} else if (strcmp(update->unit_style,"cgs") == 0) {
energyUnits=6.0221418e13;
// LAMMPS units electron :: Hartree;
} else if (strcmp(update->unit_style,"electron") == 0) {
energyUnits=2625.5257;
} else error->all(FLERR,"Fix plumed cannot handle your choice of units");
// Conversion factor from LAMMPS length units to nm (units of PLUMED)
double lengthUnits=0.1/force->angstrom;
// Conversion factor from LAMMPS time unit to ps (units of PLUMED)
double timeUnits=0.001/force->femtosecond;
p->cmd("setMDEnergyUnits",&energyUnits);
p->cmd("setMDLengthUnits",&lengthUnits);
p->cmd("setMDTimeUnits",&timeUnits);
}
// Read fix parameters:
int next=0;
for (int i=3;i<narg;++i) {
if (!strcmp(arg[i],"outfile")) {
next=1;
} else if (next==1) {
if (universe->existflag == 1) {
// Each replica writes an independent log file
// with suffix equal to the replica id
char str_num[32], logFile[1024];
sprintf(str_num,".%d",universe->iworld);
strncpy(logFile,arg[i],1024-32);
strcat(logFile,str_num);
p->cmd("setLogFile",logFile);
next=0;
} else {
// partition option not used
p->cmd("setLogFile",arg[i]);
next=0;
}
} else if (!strcmp(arg[i],"plumedfile")) {
next=2;
} else if (next==2) {
p->cmd("setPlumedDat",arg[i]);
next=0;
} else error->all(FLERR,"Syntax error - use 'fix <fix-ID> plumed "
"plumedfile plumed.dat outfile plumed.out' ");
}
if (next==1) error->all(FLERR,"missing argument for outfile option");
if (next==2) error->all(FLERR,"missing argument for plumedfile option");
p->cmd("setMDEngine","LAMMPS");
if (atom->natoms > MAXSMALLINT)
error->all(FLERR,"Fix plumed can only handle up to 2.1 billion atoms");
natoms=int(atom->natoms);
p->cmd("setNatoms",&natoms);
double dt=update->dt;
p->cmd("setTimestep",&dt);
virial_flag=1;
thermo_virial=1;
scalar_flag = 1;
// This is the real initialization:
p->cmd("init");
// Define compute to calculate potential energy
id_pe = new char[7];
id_pe = (char *) "plmd_pe";
char **newarg = new char*[3];
newarg[0] = id_pe;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pe";
modify->add_compute(3,newarg);
delete [] newarg;
int ipe = modify->find_compute(id_pe);
c_pe = modify->compute[ipe];
// Define compute to calculate pressure tensor
id_press = new char[9];
id_press = (char *) "plmd_press";
newarg = new char*[5];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = (char *) "NULL";
newarg[4] = (char *) "virial";
modify->add_compute(5,newarg);
delete [] newarg;
int ipress = modify->find_compute(id_press);
c_press = modify->compute[ipress];
for (int i = 0; i < modify->nfix; i++) {
const char * const check_style = modify->fix[i]->style;
// There must be only one
if (strcmp(check_style,"plumed") == 0)
error->all(FLERR,"There must be only one instance of fix plumed");
// Avoid conflict with fixes that define internal pressure computes.
// See comment in the setup method
if ((strncmp(check_style,"nph",3) == 0) ||
(strncmp(check_style,"npt",3) == 0) ||
(strncmp(check_style,"rigid/nph",9) == 0) ||
(strncmp(check_style,"rigid/npt",9) == 0) ||
(strncmp(check_style,"msst",4) == 0) ||
(strncmp(check_style,"nphug",5) == 0) ||
(strncmp(check_style,"ipi",3) == 0) ||
(strncmp(check_style,"press/berendsen",15) == 0) ||
(strncmp(check_style,"qbmsst",6) == 0))
error->all(FLERR,"Fix plumed must be defined before any other fixes, "
"that compute pressure internally");
}
}
FixPlumed::~FixPlumed()
{
delete p;
modify->delete_compute(id_pe);
modify->delete_compute(id_press);
delete[] masses;
delete[] charges;
delete[] gatindex;
}
int FixPlumed::setmask()
{
// set with a bitmask how and when apply the force from plumed
int mask = 0;
mask |= POST_FORCE;
mask |= THERMO_ENERGY;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
void FixPlumed::init()
{
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
// This avoids nan pressure if compute_pressure is called
// in a setup method
for (int i=0;i<6;i++) virial[i] = 0.;
}
void FixPlumed::setup(int vflag)
{
// Here there is a crucial issue connected to constant pressure
// simulations. The fix_nh will call the compute_pressure inside
// the setup method, that is executed once and for all at the
// beginning of the simulation. Since our fix has a contribution
// to the virial, when this happens the variable virial must have
// been calculated. In other words, the setup method of fix_plumed
// has to be executed first. This creates a race condition with the
// setup method of fix_nh. This is why in the constructor I check if
// nh fixes have already been called.
if (strcmp(update->integrate_style,"verlet") == 0)
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}
void FixPlumed::min_setup(int vflag)
{
// This has to be checked.
// For instance it might have problems with fix_box_relax
post_force(vflag);
}
void FixPlumed::post_force(int /* vflag */)
{
int update_gatindex=0;
if (natoms != int(atom->natoms))
error->all(FLERR,"Fix plumed does not support simulations with varying "
"numbers of atoms");
// Try to find out if the domain decomposition has been updated:
if (nlocal != atom->nlocal) {
if (charges) delete [] charges;
if (masses) delete [] masses;
if (gatindex) delete [] gatindex;
nlocal=atom->nlocal;
gatindex=new int [nlocal];
masses=new double [nlocal];
charges=new double [nlocal];
update_gatindex=1;
} else {
for (int i=0;i<nlocal;i++) {
if (gatindex[i]!=atom->tag[i]-1) {
update_gatindex=1;
break;
}
}
}
MPI_Allreduce(MPI_IN_PLACE,&update_gatindex,1,MPI_INT,MPI_SUM,world);
// In case it has been updated, rebuild the local mass/charges array
// and tell plumed about the change:
if (update_gatindex) {
for (int i=0;i<nlocal;i++) gatindex[i]=atom->tag[i]-1;
// Get masses
if (atom->rmass_flag) {
for (int i=0;i<nlocal;i++) masses[i]=atom->rmass[i];
} else {
for (int i=0;i<nlocal;i++) masses[i]=atom->mass[atom->type[i]];
}
// Get charges
if (atom->q_flag) {
for (int i=0;i<nlocal;i++) charges[i]=atom->q[i];
} else {
for (int i=0;i<nlocal;i++) charges[i]=0.0;
}
p->cmd("setAtomsNlocal",&nlocal);
p->cmd("setAtomsGatindex",gatindex);
}
// set up local virial/box. plumed uses full 3x3 matrices
double plmd_virial[3][3];
for (int i=0;i<3;i++) for (int j=0;j<3;j++) plmd_virial[i][j]=0.0;
double box[3][3];
for (int i=0;i<3;i++) for (int j=0;j<3;j++) box[i][j]=0.0;
box[0][0]=domain->h[0];
box[1][1]=domain->h[1];
box[2][2]=domain->h[2];
box[2][1]=domain->h[3];
box[2][0]=domain->h[4];
box[1][0]=domain->h[5];
// Make initial of virial of this fix zero
// The following line is very important, otherwise
// the compute pressure will include
for (int i=0;i<6;++i) virial[i] = 0.;
// local variable with timestep:
if (update->ntimestep > MAXSMALLINT)
error->all(FLERR,"Fix plumed can only handle up to 2.1 billion timesteps");
int step=int(update->ntimestep);
// pass all pointers to plumed:
p->cmd("setStep",&step);
p->cmd("setPositions",&atom->x[0][0]);
p->cmd("setBox",&box[0][0]);
p->cmd("setForces",&atom->f[0][0]);
p->cmd("setMasses",&masses[0]);
p->cmd("setCharges",&charges[0]);
p->cmd("getBias",&bias);
// Pass virial to plumed
// If energy is needed virial_plmd is equal to Lammps' virial
// If energy is not needed virial_plmd is initialized to zero
// In the first case the virial will be rescaled and an extra term will be added
// In the latter case only an extra term will be added
p->cmd("setVirial",&plmd_virial[0][0]);
p->cmd("prepareCalc");
plumedNeedsEnergy=0;
p->cmd("isEnergyNeeded",&plumedNeedsEnergy);
// Pass potential energy and virial if needed
double *virial_lmp;
if (plumedNeedsEnergy) {
// Error if tail corrections are included
if (force->pair && force->pair->tail_flag && comm->me == 0)
error->warning(FLERR,"Tail corrections to the pair potential included."
" The energy cannot be biased correctly in this case."
" Remove the tail corrections by removing the"
" command: pair_modify tail yes");
// compute the potential energy
double pot_energy = 0.;
c_pe->compute_scalar();
pot_energy = c_pe->scalar;
// Divide energy by number of processes
// Plumed wants it this way
int nprocs;
MPI_Comm_size(world,&nprocs);
pot_energy /= nprocs;
p->cmd("setEnergy",&pot_energy);
// Compute pressure due to the virial (no kinetic energy term!)
c_press->compute_vector();
virial_lmp = c_press->vector;
// Check if pressure is finite
if (!std::isfinite(virial_lmp[0]) || !std::isfinite(virial_lmp[1])
|| !std::isfinite(virial_lmp[2]) || !std::isfinite(virial_lmp[3])
|| !std::isfinite(virial_lmp[4]) || !std::isfinite(virial_lmp[5]))
error->all(FLERR,"Non-numeric virial - Plumed cannot work with that");
// Convert pressure to virial per number of MPI processes
// From now on all virials are divided by the number of MPI processes
double nktv2p = force->nktv2p;
double inv_volume;
if (domain->dimension == 3) {
inv_volume = 1.0 / (domain->xprd * domain->yprd * domain->zprd);
} else {
inv_volume = 1.0 / (domain->xprd * domain->yprd);
}
for (int i=0;i<6;i++) virial_lmp[i] /= (inv_volume * nktv2p * nprocs);
// Convert virial from lammps to plumed representation
plmd_virial[0][0]=-virial_lmp[0];
plmd_virial[1][1]=-virial_lmp[1];
plmd_virial[2][2]=-virial_lmp[2];
plmd_virial[0][1]=-virial_lmp[3];
plmd_virial[0][2]=-virial_lmp[4];
plmd_virial[1][2]=-virial_lmp[5];
} else {
virial_lmp = new double[6];
for (int i=0;i<6;i++) virial_lmp[i] = 0.;
}
// do the real calculation:
p->cmd("performCalc");
// retransform virial to lammps representation and assign it to this
// fix's virial. Plumed is giving back the full virial and therefore
// we have to subtract the initial virial i.e. virial_lmp.
// The vector virial contains only the contribution added by plumed.
// The calculation of the pressure will be done by a compute pressure
// and will include this contribution.
virial[0] = -plmd_virial[0][0]-virial_lmp[0];
virial[1] = -plmd_virial[1][1]-virial_lmp[1];
virial[2] = -plmd_virial[2][2]-virial_lmp[2];
virial[3] = -plmd_virial[0][1]-virial_lmp[3];
virial[4] = -plmd_virial[0][2]-virial_lmp[4];
virial[5] = -plmd_virial[1][2]-virial_lmp[5];
// Ask for the computes in the next time step
// such that the virial and energy are tallied.
// This should be changed to something that triggers the
// calculation only if plumed needs it.
c_pe->addstep(update->ntimestep+1);
c_press->addstep(update->ntimestep+1);
}
void FixPlumed::post_force_respa(int vflag, int ilevel, int /* iloop */)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
void FixPlumed::min_post_force(int vflag)
{
post_force(vflag);
}
void FixPlumed::reset_dt()
{
error->all(FLERR,"Cannot change the time step when fix plumed is active");
}
double FixPlumed::compute_scalar()
{
return bias;
}
int FixPlumed::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"pe") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
modify->delete_compute(id_pe); delete [] id_pe;
int n = strlen(arg[1]) + 1;
id_pe = new char[n];
strcpy(id_pe,arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify potential energy ID");
c_pe = modify->compute[icompute];
if (c_pe->peflag == 0)
error->all(FLERR,"Fix_modify plmd_pe ID does not compute potential energy");
if (c_pe->igroup != 0 && comm->me == 0)
error->warning(FLERR,"Potential for fix PLUMED is not for group all");
return 2;
} else if (strcmp(arg[0],"press") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
modify->delete_compute(id_press); delete [] id_press;
int n = strlen(arg[1]) + 1;
id_press = new char[n];
strcpy(id_press,arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");
c_press = modify->compute[icompute];
if (c_press->pressflag == 0)
error->all(FLERR,"Fix_modify pressure ID does not compute pressure");
if (c_press->igroup != 0 && comm->me == 0)
error->warning(FLERR,"Virial for fix PLUMED is not for group all");
return 2;
}
return 0;
}
double FixPlumed::memory_usage()
{
return double((8+8+4)*atom->nlocal);
}

68
src/USER-PLUMED/fix_plumed.h Normal file
View File

@ -0,0 +1,68 @@
/* -*- 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 FIX_CLASS
FixStyle(plumed,FixPlumed)
#else
#ifndef LMP_FIX_PLUMED_H
#define LMP_FIX_PLUMED_H
#include "fix.h"
// forward declaration
namespace PLMD {
class Plumed;
}
namespace LAMMPS_NS {
class FixPlumed : public Fix {
public:
FixPlumed(class LAMMPS *, int, char **);
~FixPlumed();
int setmask();
void init();
void setup(int);
void min_setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
double compute_scalar();
void reset_dt();
int modify_param(int narg, char **arg);
double memory_usage();
private:
PLMD::Plumed *p; // pointer to plumed object
int nlocal; // number of atoms local to this process
int natoms; // total number of atoms
int *gatindex; // array of atom indexes local to this process
double *masses; // array of masses for local atoms
double *charges; // array of charges for local atoms
int nlevels_respa; // this is something to enable respa
double bias; // output bias potential
class Compute *c_pe; // Compute for the energy
class Compute *c_press; // Compute for the pressure
int plumedNeedsEnergy; // Flag to trigger calculation of the
// energy and virial
char *id_pe, *id_press; // ID for potential energy and pressure compute
};
};
#endif
#endif

1
src/VORONOI/compute_voronoi_atom.cpp
View File

@ -19,6 +19,7 @@
#include <cmath>
#include <cstring>
#include <cstdlib>
#include "voro++.hh"
#include "compute_voronoi_atom.h"
#include "atom.h"
#include "group.h"

7
src/VORONOI/compute_voronoi_atom.h
View File

@ -21,7 +21,12 @@ ComputeStyle(voronoi/atom,ComputeVoronoi)
#define LMP_COMPUTE_VORONOI_H
#include "compute.h"
#include "voro++.hh"
namespace voro {
class container;
class container_poly;
class voronoicell_neighbor;
}
namespace LAMMPS_NS {

2
src/comm_brick.cpp
View File

@ -1477,8 +1477,10 @@ void CommBrick::free_multi()
void *CommBrick::extract(const char *str, int &dim)
{
dim = 0;
if (strcmp(str,"localsendlist") == 0) {
int i, iswap, isend;
dim = 1;
if (!localsendlist)
memory->create(localsendlist,atom->nlocal,"comm:localsendlist");
else

23
src/compute_adf.cpp
View File

@ -46,14 +46,17 @@ enum{DEGREE, RADIAN, COSINE};
ComputeADF::ComputeADF(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
ilo(NULL), ihi(NULL), jlo(NULL), jhi(NULL), klo(NULL), khi(NULL),
iatomcount(NULL), iatomcountall(NULL),
hist(NULL), histall(NULL),
iatomflag(NULL),
jatomflag(NULL), rcutinnerj(NULL), rcutouterj(NULL),
katomflag(NULL), rcutinnerk(NULL), rcutouterk(NULL),
maxjatom(NULL), numjatom(NULL), neighjatom(NULL),
maxkatom(NULL), numkatom(NULL), neighkatom(NULL),
maxjkatom(NULL), numjkatom(NULL), neighjkatom(NULL), bothjkatom(NULL)
rcutinnerj(NULL), rcutinnerk(NULL),
rcutouterj(NULL), rcutouterk(NULL),
list(NULL),
iatomcount(NULL), iatomcountall(NULL), iatomflag(NULL),
maxjatom(NULL), maxkatom(NULL),
numjatom(NULL), numkatom(NULL),
neighjatom(NULL),neighkatom(NULL),
jatomflag(NULL), katomflag(NULL),
maxjkatom(NULL), numjkatom(NULL),
neighjkatom(NULL), bothjkatom(NULL), delrjkatom(NULL)
{
int nargsperadf = 7;
@ -358,9 +361,9 @@ void ComputeADF::init_list(int /*id*/, NeighList *ptr)
void ComputeADF::compute_array()
{
int i,j,k,m,ii,jj,jatom,katom,jk,jjj,kkk;
int inum,jnum,itype,jtype,ibin,ihisto;
int inum,jnum,itype,jtype,ibin;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *ilist,*jlist,*klist,*numneigh,**firstneigh;
int *ilist,*jlist,*numneigh,**firstneigh;
double factor_lj,factor_coul;
double delr1[3],delr2[3],rinv1,rinv2,rinv12,cs,theta;
@ -394,11 +397,9 @@ void ComputeADF::compute_array()
double **x = atom->x;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double *special_coul = force->special_coul;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];

2
src/compute_angle_local.cpp
View File

@ -38,7 +38,7 @@ enum{THETA,ENG,VARIABLE};
ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
bstyle(NULL), vstr(NULL), vvar(NULL), tstr(NULL), vlocal(NULL), alocal(NULL)
bstyle(NULL), vvar(NULL), tstr(NULL), vstr(NULL), vlocal(NULL), alocal(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal compute angle/local command");

2
src/compute_bond_local.cpp
View File

@ -39,7 +39,7 @@ enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,VARIABLE};
ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
bstyle(NULL), vstr(NULL), vvar(NULL), dstr(NULL), vlocal(NULL), alocal(NULL)
bstyle(NULL), vvar(NULL), dstr(NULL), vstr(NULL), vlocal(NULL), alocal(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal compute bond/local command");

2
src/compute_dihedral_local.cpp
View File

@ -40,7 +40,7 @@ enum{PHI,VARIABLE};
ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
bstyle(NULL), vstr(NULL), vvar(NULL), pstr(NULL), vlocal(NULL), alocal(NULL)
bstyle(NULL), vvar(NULL), pstr(NULL), vstr(NULL), vlocal(NULL), alocal(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal compute dihedral/local command");

2
src/deprecated.cpp
View File

@ -36,7 +36,7 @@ static void writemsg(LAMMPS *lmp, const char *msg, int abend=1)
/* ---------------------------------------------------------------------- */
void Deprecated::command(int narg, char **arg)
void Deprecated::command(int /* narg */, char ** /* arg */)
{
if (strcmp(input->command,"DEPRECATED") == 0) {
writemsg(lmp,"\nCommand 'DEPRECATED' is a dummy command\n\n",0);

4
tools/msi2lmp/src/MakeLists.c
View File

@ -363,7 +363,7 @@ int count_dihedrals()
if (i != k) {
for (ll=0; ll < atoms[k].no_connect; ll++) {
l = atoms[k].conn_no[ll];
if (l != j) n++;
if ((l != j) && (i != l)) n++
}
}
}
@ -391,7 +391,7 @@ void build_dihedrals_list()
if (i != k) {
for (ll=0; ll < atoms[k].no_connect; ll++) {
l = atoms[k].conn_no[ll];
if (l != j) {
if ((l != j) && (i != l)) {
dihedrals[n ].type = 0;
dihedrals[n ].members[0] = i;
dihedrals[n ].members[1] = j;

2
tools/msi2lmp/src/msi2lmp.c
View File

@ -2,6 +2,8 @@
*
* msi2lmp.exe
*
* v3.9.9 AK- Teach msi2lmp to not generate dihedrals with identical 1-4 atoms
*
* v3.9.8 AK- Improved whitespace handling in parsing topology and force
* field files to avoid bogus warnings about type name truncation
*

2
tools/msi2lmp/src/msi2lmp.h
View File

@ -36,7 +36,7 @@
# include <stdio.h>
#define MSI2LMP_VERSION "v3.9.8 / 06 Oct 2016"
#define MSI2LMP_VERSION "v3.9.9 / 05 Nov 2018"
#define PI_180 0.01745329251994329576