forked from lijiext/lammps
Merge branch 'master' into polymorphic-update
This commit is contained in:
Axel Kohlmeyer
commit
21e9db6bdf
.github
cmake/Modules/Packages
doc
lib/colvars
src
|
|
@ -22,7 +22,6 @@ src/SPIN/* @julient31
|
|||
src/USER-CGDNA/* @ohenrich
|
||||
src/USER-CGSDK/* @akohlmey
|
||||
src/USER-COLVARS/* @giacomofiorin
|
||||
src/USER-DPD/* @timattox
|
||||
src/USER-INTEL/* @wmbrownintel
|
||||
src/USER-MANIFOLD/* @Pakketeretet2
|
||||
src/USER-MEAMC/* @martok
|
||||
|
|
@ -119,6 +118,7 @@ tools/emacs/* @HaoZeke
|
|||
|
||||
# cmake
|
||||
cmake/* @junghans @rbberger
|
||||
cmake/Modules/Packages/USER-COLVARS.cmake @junghans @rbberger @giacomofiorin
|
||||
|
||||
# python
|
||||
python/* @rbberger
|
||||
|
|
|
|||
|
|
@ -12,6 +12,8 @@ if(COLVARS_LEPTON)
|
|||
if(NOT BUILD_SHARED_LIBS)
|
||||
install(TARGETS lepton EXPORT LAMMPS_Targets LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR})
|
||||
endif()
|
||||
# Change the define below to LEPTON_BUILDING_SHARED_LIBRARY when linking Lepton as a DLL with MSVC
|
||||
target_compile_definitions(lepton PRIVATE -DLEPTON_BUILDING_STATIC_LIBRARY)
|
||||
set_target_properties(lepton PROPERTIES OUTPUT_NAME lammps_lepton${LAMMPS_MACHINE})
|
||||
target_include_directories(lepton PRIVATE ${LEPTON_DIR}/include)
|
||||
endif()
|
||||
|
|
@ -27,6 +29,8 @@ target_link_libraries(lammps PRIVATE colvars)
|
|||
|
||||
if(COLVARS_LEPTON)
|
||||
target_link_libraries(lammps PRIVATE lepton)
|
||||
target_compile_options(colvars PRIVATE -DLEPTON)
|
||||
target_compile_definitions(colvars PRIVATE -DLEPTON)
|
||||
# Disable the line below when linking Lepton as a DLL with MSVC
|
||||
target_compile_definitions(colvars PRIVATE -DLEPTON_USE_STATIC_LIBRARIES)
|
||||
target_include_directories(colvars PUBLIC ${LEPTON_DIR}/include)
|
||||
endif()
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
.TH LAMMPS "19 March 2020" "2020-03-19"
|
||||
.TH LAMMPS "15 April 2020" "2020-04-15"
|
||||
.SH NAME
|
||||
.B LAMMPS
|
||||
\- Molecular Dynamics Simulator.
|
||||
|
|
|
|||
|
|
@ -142,24 +142,40 @@ new owning processors.
|
|||
|
||||
.. note::
|
||||
|
||||
The simulation box size/shape can be changed by arbitrarily
|
||||
large amounts by this command. This is not a problem, except that the
|
||||
The simulation box size/shape can be changed by arbitrarily large
|
||||
amounts by this command. This is not a problem, except that the
|
||||
mapping of processors to the simulation box is not changed from its
|
||||
initial 3d configuration; see the :doc:`processors <processors>`
|
||||
command. Thus, if the box size/shape changes dramatically, the
|
||||
mapping of processors to the simulation box may not end up as optimal
|
||||
as the initial mapping attempted to be.
|
||||
mapping of processors to the simulation box may not end up as
|
||||
optimal as the initial mapping attempted to be. You may wish to
|
||||
re-balance the atoms by using the :doc:`balance <balance>` command
|
||||
if that is the case.
|
||||
|
||||
.. note::
|
||||
|
||||
Because the keywords used in this command are applied one at a
|
||||
time to the simulation box and the atoms in it, care must be taken
|
||||
with triclinic cells to avoid exceeding the limits on skew after each
|
||||
transformation in the sequence. If skew is exceeded before the final
|
||||
transformation this can be avoided by changing the order of the
|
||||
sequence, or breaking the transformation into two or more smaller
|
||||
transformations. For more information on the allowed limits for box
|
||||
skew see the discussion on triclinic boxes on :doc:`Howto triclinic <Howto_triclinic>` doc page.
|
||||
You cannot use this command after reading a restart file (and
|
||||
before a run is performed) if the restart file stored per-atom
|
||||
information from a fix and any of the specified keywords change the
|
||||
box size or shape or boundary conditions. This is because atoms
|
||||
may be moved to new processors and the restart info will not
|
||||
migrate with them. LAMMPS will generate an error if this could
|
||||
happen. Only the *ortho* and *triclinic* keywords do not trigger
|
||||
this error. One solution is to perform a "run 0" command before
|
||||
using the change_box command. This clears the per-atom restart
|
||||
data, whether it has been re-assigned to a new fix or not.
|
||||
|
||||
.. note::
|
||||
|
||||
Because the keywords used in this command are applied one at a time
|
||||
to the simulation box and the atoms in it, care must be taken with
|
||||
triclinic cells to avoid exceeding the limits on skew after each
|
||||
transformation in the sequence. If skew is exceeded before the
|
||||
final transformation this can be avoided by changing the order of
|
||||
the sequence, or breaking the transformation into two or more
|
||||
smaller transformations. For more information on the allowed
|
||||
limits for box skew see the discussion on triclinic boxes on
|
||||
:doc:`Howto triclinic <Howto_triclinic>` doc page.
|
||||
|
||||
----------
|
||||
|
||||
|
|
|
|||
|
|
@ -17,10 +17,18 @@ Syntax
|
|||
|
||||
.. parsed-literal::
|
||||
|
||||
*bond* args = atom1 atom2 Kstart Kstop r0
|
||||
*bond* args = atom1 atom2 Kstart Kstop r0start (r0stop)
|
||||
atom1,atom2 = IDs of 2 atoms in bond
|
||||
Kstart,Kstop = restraint coefficients at start/end of run (energy units)
|
||||
r0 = equilibrium bond distance (distance units)
|
||||
r0start = equilibrium bond distance at start of run (distance units)
|
||||
r0stop = equilibrium bond distance at end of run (optional) (distance units). If not
|
||||
specified it is assumed to be equal to r0start
|
||||
*lbond* args = atom1 atom2 Kstart Kstop r0start (r0stop)
|
||||
atom1,atom2 = IDs of 2 atoms in bond
|
||||
Kstart,Kstop = restraint coefficients at start/end of run (energy units)
|
||||
r0start = equilibrium bond distance at start of run (distance units)
|
||||
r0stop = equilibrium bond distance at end of run (optional) (distance units). If not
|
||||
specified it is assumed to be equal to r0start
|
||||
*angle* args = atom1 atom2 atom3 Kstart Kstop theta0
|
||||
atom1,atom2,atom3 = IDs of 3 atoms in angle, atom2 = middle atom
|
||||
Kstart,Kstop = restraint coefficients at start/end of run (energy units)
|
||||
|
|
@ -38,6 +46,7 @@ Examples
|
|||
.. code-block:: LAMMPS
|
||||
|
||||
fix holdem all restrain bond 45 48 2000.0 2000.0 2.75
|
||||
fix holdem all restrain lbond 45 48 2000.0 2000.0 2.75
|
||||
fix holdem all restrain dihedral 1 2 3 4 2000.0 2000.0 120.0
|
||||
fix holdem all restrain bond 45 48 2000.0 2000.0 2.75 dihedral 1 2 3 4 2000.0 2000.0 120.0
|
||||
fix texas_holdem all restrain dihedral 1 2 3 4 0.0 2000.0 120.0 dihedral 1 2 3 5 0.0 2000.0 -120.0 dihedral 1 2 3 6 0.0 2000.0 0.0
|
||||
|
|
@ -141,6 +150,29 @@ is included in :math:`K`.
|
|||
|
||||
----------
|
||||
|
||||
The *lbond* keyword applies a lower bound bond restraint to the specified atoms
|
||||
using the same functional form used by the :doc:`bond_style harmonic <bond_harmonic>` command if the distance between
|
||||
the atoms is smaller than the equilibrium bond distance and 0 otherwise. The potential associated with
|
||||
the restraint is
|
||||
|
||||
.. math::
|
||||
|
||||
E = K (r - r_0)^2 ,if\;r < r_0
|
||||
|
||||
.. math::
|
||||
|
||||
E = 0 \qquad\quad\quad ,if\;r \ge r_0
|
||||
|
||||
with the following coefficients:
|
||||
|
||||
* :math:`K` (energy/distance\^2)
|
||||
* :math:`r_0` (distance)
|
||||
|
||||
:math:`K` and :math:`r_0` are specified with the fix. Note that the usual 1/2 factor
|
||||
is included in :math:`K`.
|
||||
|
||||
----------
|
||||
|
||||
The *angle* keyword applies an angle restraint to the specified atoms
|
||||
using the same functional form used by the :doc:`angle_style harmonic <angle_harmonic>` command. The potential associated with
|
||||
the restraint is
|
||||
|
|
|
|||
|
|
@ -61,7 +61,7 @@ ifeq ($(COLVARS_LEPTON),no)
|
|||
LEPTON_INCFLAGS =
|
||||
COLVARS_OBJS = $(COLVARS_SRCS:.cpp=.o)
|
||||
else
|
||||
LEPTON_INCFLAGS = -Ilepton/include -DLEPTON
|
||||
LEPTON_INCFLAGS = -Ilepton/include -DLEPTON -DLEPTON_USE_STATIC_LIBRARIES
|
||||
COLVARS_OBJS = $(COLVARS_SRCS:.cpp=.o) $(LEPTON_SRCS:.cpp=.o)
|
||||
endif
|
||||
|
||||
|
|
@ -82,4 +82,20 @@ Makefile.deps: $(COLVARS_SRCS)
|
|||
done
|
||||
|
||||
include Makefile.deps
|
||||
|
||||
# Exceptions to pattern rule above for Lepton objects
|
||||
|
||||
lepton/src/CompiledExpression.o: lepton/src/CompiledExpression.cpp
|
||||
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
|
||||
lepton/src/ExpressionProgram.o: lepton/src/ExpressionProgram.cpp
|
||||
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
|
||||
lepton/src/ExpressionTreeNode.o: lepton/src/ExpressionTreeNode.cpp
|
||||
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
|
||||
lepton/src/Operation.o: lepton/src/Operation.cpp
|
||||
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
|
||||
lepton/src/ParsedExpression.o: lepton/src/ParsedExpression.cpp
|
||||
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
|
||||
lepton/src/Parser.o: lepton/src/Parser.cpp
|
||||
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
|
||||
|
||||
include Makefile.lepton.deps # Hand-generated
|
||||
|
|
|
|||
|
|
@ -15,6 +15,7 @@
|
|||
Contributing authors: Axel Kohlmeyer (Temple U),
|
||||
Ryan S. Elliott (UMN)
|
||||
Ellad B. Tadmor (UMN)
|
||||
Ronald Miller (Carleton U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -57,8 +58,10 @@
|
|||
|
||||
#include "kim_interactions.h"
|
||||
#include <cstring>
|
||||
#include <cstdio>
|
||||
#include <string>
|
||||
#include <sstream>
|
||||
#include <vector>
|
||||
#include "error.h"
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
|
|
@ -80,6 +83,8 @@ extern "C" {
|
|||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
#define MAXLINE 1024
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void KimInteractions::command(int narg, char **arg)
|
||||
|
|
@ -223,8 +228,27 @@ void KimInteractions::do_setup(int narg, char **arg)
|
|||
for (int j=0; j < sim_lines; ++j) {
|
||||
KIM_SimulatorModel_GetSimulatorFieldLine(
|
||||
simulatorModel,sim_model_idx,j,&sim_value);
|
||||
input->one(sim_value);
|
||||
}
|
||||
char strbuf[MAXLINE];
|
||||
char * strword;
|
||||
strcpy(strbuf,sim_value);
|
||||
strword = strtok(strbuf," \t");
|
||||
if (0==strcmp(strword,"KIM_MATCH_PAIRS")) {
|
||||
// Notes regarding the KIM_MATCH_PAIRS command
|
||||
// * This is an INTERNAL command.
|
||||
// * It is intended for use only by KIM Simulator Models.
|
||||
// * It is not possible to use this command outside of the context
|
||||
// of the kim_interactions command and KIM Simulator Models.
|
||||
// * The command performs a transformation from symbolic
|
||||
// string-based atom types to lammps numeric atom types for
|
||||
// the pair_coeff settings.
|
||||
// * The command is not documented fully as it is expected to be
|
||||
// temporary. Eventually it should be replaced by a more
|
||||
// comprehensive symbolic types support in lammps.
|
||||
KIM_MATCH_PAIRS(sim_value);
|
||||
} else {
|
||||
input->one(sim_value);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -263,6 +287,73 @@ void KimInteractions::do_setup(int narg, char **arg)
|
|||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void KimInteractions::KIM_MATCH_PAIRS(char const *const input_line) const
|
||||
{
|
||||
char strbuf[MAXLINE];
|
||||
strcpy(strbuf,input_line);
|
||||
char *cmd, *filename;
|
||||
cmd = strtok(strbuf," \t");
|
||||
filename = strtok(NULL," \t");
|
||||
|
||||
FILE *fp;
|
||||
fp = fopen(filename,"r");
|
||||
if (fp == NULL) {
|
||||
error->one(FLERR,"Parameter file not found");
|
||||
}
|
||||
|
||||
std::vector<char *> species;
|
||||
for (int i = 0; i < atom->ntypes; ++i)
|
||||
{
|
||||
char *str;
|
||||
str = strtok(NULL," \t");
|
||||
if (str == NULL)
|
||||
error->one(FLERR,"Incorrect args for KIM_MATCH_PAIRS command");
|
||||
species.push_back(str);
|
||||
}
|
||||
|
||||
char line[MAXLINE],*ptr;
|
||||
int n, eof = 0;
|
||||
|
||||
while (1) {
|
||||
if (comm->me == 0) {
|
||||
ptr = fgets(line,MAXLINE,fp);
|
||||
if (ptr == NULL) {
|
||||
eof = 1;
|
||||
fclose(fp);
|
||||
} else n = strlen(line) + 1;
|
||||
}
|
||||
MPI_Bcast(&eof,1,MPI_INT,0,world);
|
||||
if (eof) break;
|
||||
MPI_Bcast(&n,1,MPI_INT,0,world);
|
||||
MPI_Bcast(line,n,MPI_CHAR,0,world);
|
||||
|
||||
char *species1, *species2, *the_rest;
|
||||
ptr = line;
|
||||
species1 = strtok(ptr," \t");
|
||||
species2 = strtok(NULL," \t");
|
||||
the_rest = strtok(NULL,"\n");
|
||||
|
||||
for (int type_a = 0; type_a < atom->ntypes; ++type_a) {
|
||||
for (int type_b = type_a; type_b < atom->ntypes; ++type_b) {
|
||||
if(((strcmp(species[type_a],species1) == 0) &&
|
||||
(strcmp(species[type_b],species2) == 0))
|
||||
||
|
||||
((strcmp(species[type_b],species1) == 0) &&
|
||||
(strcmp(species[type_a],species2) == 0))
|
||||
) {
|
||||
char pair_command[MAXLINE];
|
||||
sprintf(pair_command,"pair_coeff %i %i %s",type_a+1,type_b+1,
|
||||
the_rest);
|
||||
input->one(pair_command);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
fclose(fp);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int KimInteractions::species_to_atomic_no(std::string const species) const
|
||||
{
|
||||
if (species == "H") return 1;
|
||||
|
|
|
|||
|
|
@ -15,6 +15,7 @@
|
|||
Contributing authors: Axel Kohlmeyer (Temple U),
|
||||
Ryan S. Elliott (UMN)
|
||||
Ellad B. Tadmor (UMN)
|
||||
Ronald Miller (Carleton U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
|
@ -76,6 +77,7 @@ class KimInteractions : protected Pointers {
|
|||
private:
|
||||
void do_setup(int, char **);
|
||||
int species_to_atomic_no(std::string const species) const;
|
||||
void KIM_MATCH_PAIRS(char const *const input_line) const;
|
||||
void kim_interactions_log_delimiter(std::string const begin_end) const;
|
||||
};
|
||||
|
||||
|
|
|
|||
|
|
@ -636,7 +636,7 @@ double MinKokkos::fnorm_inf()
|
|||
|
||||
double MinKokkos::fnorm_max()
|
||||
{
|
||||
|
||||
|
||||
double local_norm_max = 0.0;
|
||||
{
|
||||
// local variables for lambda capture
|
||||
|
|
|
|||
|
|
@ -31,7 +31,6 @@
|
|||
#include "neigh_request.h"
|
||||
#include "math_special.h"
|
||||
#include "pair_dpd_fdt_energy.h"
|
||||
#include "utils.h"
|
||||
|
||||
#include <vector> // std::vector<>
|
||||
#include <algorithm> // std::max
|
||||
|
|
@ -256,10 +255,6 @@ void FixRX::post_constructor()
|
|||
int nUniqueSpecies = 0;
|
||||
bool match;
|
||||
|
||||
for (int i = 0; i < modify->nfix; i++)
|
||||
if (utils::strmatch(modify->fix[i]->style,"^property/atom") == 0)
|
||||
error->all(FLERR,"fix rx cannot be combined with fix property/atom");
|
||||
|
||||
char **tmpspecies = new char*[maxspecies];
|
||||
int tmpmaxstrlen = 0;
|
||||
for(int jj=0; jj < maxspecies; jj++)
|
||||
|
|
|
|||
|
|
@ -46,9 +46,6 @@ void ChangeBox::command(int narg, char **arg)
|
|||
if (domain->box_exist == 0)
|
||||
error->all(FLERR,"Change_box command before simulation box is defined");
|
||||
if (narg < 2) error->all(FLERR,"Illegal change_box command");
|
||||
if (modify->nfix_restart_peratom)
|
||||
error->all(FLERR,"Cannot change_box after "
|
||||
"reading restart file with per-atom info");
|
||||
|
||||
if (comm->me == 0 && screen) fprintf(screen,"Changing box ...\n");
|
||||
|
||||
|
|
@ -174,6 +171,21 @@ void ChangeBox::command(int narg, char **arg)
|
|||
|
||||
if (nops == 0) error->all(FLERR,"Illegal change_box command");
|
||||
|
||||
// move_atoms = 1 if need to move atoms to new procs after box changes
|
||||
// anything other than ORTHO or TRICLINIC may cause atom movement
|
||||
|
||||
int move_atoms = 0;
|
||||
for (int m = 0; m < nops; m++) {
|
||||
if (ops[m].style != ORTHO || ops[m].style != TRICLINIC) move_atoms = 1;
|
||||
}
|
||||
|
||||
// error if moving atoms and there is stored per-atom restart state
|
||||
// disallowed b/c restart per-atom fix info will not move with atoms
|
||||
|
||||
if (move_atoms && modify->nfix_restart_peratom)
|
||||
error->all(FLERR,"Change_box parameter not allowed after "
|
||||
"reading restart file with per-atom info");
|
||||
|
||||
// read options from end of input line
|
||||
|
||||
options(narg-iarg,&arg[iarg]);
|
||||
|
|
@ -350,6 +362,10 @@ void ChangeBox::command(int narg, char **arg)
|
|||
if (domain->triclinic) domain->lamda2x(atom->nlocal);
|
||||
}
|
||||
|
||||
// done if don't need to move atoms
|
||||
|
||||
if (!move_atoms) return;
|
||||
|
||||
// move atoms back inside simulation box and to new processors
|
||||
// use remap() instead of pbc()
|
||||
// in case box moved a long distance relative to atoms
|
||||
|
|
|
|||
|
|
@ -34,7 +34,7 @@ using namespace LAMMPS_NS;
|
|||
using namespace FixConst;
|
||||
using namespace MathConst;
|
||||
|
||||
enum{BOND,ANGLE,DIHEDRAL};
|
||||
enum{BOND,LBOUND,ANGLE,DIHEDRAL};
|
||||
|
||||
#define TOLERANCE 0.05
|
||||
#define SMALL 0.001
|
||||
|
|
@ -45,7 +45,7 @@ enum{BOND,ANGLE,DIHEDRAL};
|
|||
FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg),
|
||||
rstyle(NULL), mult(NULL), ids(NULL), kstart(NULL), kstop(NULL), target(NULL),
|
||||
cos_target(NULL), sin_target(NULL)
|
||||
deqstart(NULL), deqstop(NULL), cos_target(NULL), sin_target(NULL)
|
||||
{
|
||||
if (narg < 4) error->all(FLERR,"Illegal fix restrain command");
|
||||
|
||||
|
|
@ -72,6 +72,8 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
|
|||
memory->grow(kstart,maxrestrain,"restrain:kstart");
|
||||
memory->grow(kstop,maxrestrain,"restrain:kstop");
|
||||
memory->grow(target,maxrestrain,"restrain:target");
|
||||
memory->grow(deqstart,maxrestrain,"restrain:deqstart");
|
||||
memory->grow(deqstop,maxrestrain,"restrain:deqstop");
|
||||
memory->grow(cos_target,maxrestrain,"restrain:cos_target");
|
||||
memory->grow(sin_target,maxrestrain,"restrain:sin_target");
|
||||
}
|
||||
|
|
@ -83,8 +85,29 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
|
|||
ids[nrestrain][1] = force->inumeric(FLERR,arg[iarg+2]);
|
||||
kstart[nrestrain] = force->numeric(FLERR,arg[iarg+3]);
|
||||
kstop[nrestrain] = force->numeric(FLERR,arg[iarg+4]);
|
||||
target[nrestrain] = force->numeric(FLERR,arg[iarg+5]);
|
||||
iarg += 6;
|
||||
deqstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]);
|
||||
if (iarg+6 == narg) {
|
||||
deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+5]);
|
||||
iarg += 6;
|
||||
} else {
|
||||
deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+6]);
|
||||
iarg += 7;
|
||||
}
|
||||
} else if (strcmp(arg[iarg],"lbound") == 0) {
|
||||
if (iarg+6 > narg) error->all(FLERR,"Illegal fix restrain command");
|
||||
rstyle[nrestrain] = LBOUND;
|
||||
ids[nrestrain][0] = force->inumeric(FLERR,arg[iarg+1]);
|
||||
ids[nrestrain][1] = force->inumeric(FLERR,arg[iarg+2]);
|
||||
kstart[nrestrain] = force->numeric(FLERR,arg[iarg+3]);
|
||||
kstop[nrestrain] = force->numeric(FLERR,arg[iarg+4]);
|
||||
deqstart[nrestrain] = force->numeric(FLERR,arg[iarg+5]);
|
||||
if (iarg+6 == narg) {
|
||||
deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+5]);
|
||||
iarg += 6;
|
||||
} else {
|
||||
deqstop[nrestrain] = force->numeric(FLERR,arg[iarg+6]);
|
||||
iarg += 7;
|
||||
}
|
||||
} else if (strcmp(arg[iarg],"angle") == 0) {
|
||||
if (iarg+7 > narg) error->all(FLERR,"Illegal fix restrain command");
|
||||
rstyle[nrestrain] = ANGLE;
|
||||
|
|
@ -139,6 +162,8 @@ FixRestrain::~FixRestrain()
|
|||
memory->destroy(kstart);
|
||||
memory->destroy(kstop);
|
||||
memory->destroy(target);
|
||||
memory->destroy(deqstart);
|
||||
memory->destroy(deqstop);
|
||||
memory->destroy(cos_target);
|
||||
memory->destroy(sin_target);
|
||||
}
|
||||
|
|
@ -192,11 +217,13 @@ void FixRestrain::post_force(int /*vflag*/)
|
|||
energy = 0.0;
|
||||
|
||||
ebond = 0.0;
|
||||
elbound = 0.0;
|
||||
eangle = 0.0;
|
||||
edihed = 0.0;
|
||||
|
||||
for (int m = 0; m < nrestrain; m++)
|
||||
if (rstyle[m] == BOND) restrain_bond(m);
|
||||
else if (rstyle[m] == LBOUND) restrain_lbound(m);
|
||||
else if (rstyle[m] == ANGLE) restrain_angle(m);
|
||||
else if (rstyle[m] == DIHEDRAL) restrain_dihedral(m);
|
||||
}
|
||||
|
|
@ -233,6 +260,7 @@ void FixRestrain::restrain_bond(int m)
|
|||
double delta = update->ntimestep - update->beginstep;
|
||||
if (delta != 0.0) delta /= update->endstep - update->beginstep;
|
||||
double k = kstart[m] + delta * (kstop[m] - kstart[m]);
|
||||
double deq = deqstart[m] + delta * (deqstop[m] - deqstart[m]);
|
||||
|
||||
i1 = atom->map(ids[m][0]);
|
||||
i2 = atom->map(ids[m][1]);
|
||||
|
|
@ -269,7 +297,7 @@ void FixRestrain::restrain_bond(int m)
|
|||
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
r = sqrt(rsq);
|
||||
dr = r - target[m];
|
||||
dr = r - deq;
|
||||
rk = k * dr;
|
||||
|
||||
// force & energy
|
||||
|
|
@ -277,7 +305,7 @@ void FixRestrain::restrain_bond(int m)
|
|||
if (r > 0.0) fbond = -2.0*rk/r;
|
||||
else fbond = 0.0;
|
||||
|
||||
ebond += rk*dr;
|
||||
ebond += rk*dr;
|
||||
energy += rk*dr;
|
||||
|
||||
// apply force to each of 2 atoms
|
||||
|
|
@ -295,6 +323,94 @@ void FixRestrain::restrain_bond(int m)
|
|||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
apply harmonic lower-bound bond restraints
|
||||
---------------------------------------------------------------------- */
|
||||
|
||||
void FixRestrain::restrain_lbound(int m)
|
||||
{
|
||||
int i1,i2;
|
||||
double delx,dely,delz,fbond;
|
||||
double rsq,r,dr,rk;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int nlocal = atom->nlocal;
|
||||
int newton_bond = force->newton_bond;
|
||||
|
||||
double delta = update->ntimestep - update->beginstep;
|
||||
if (delta != 0.0) delta /= update->endstep - update->beginstep;
|
||||
double k = kstart[m] + delta * (kstop[m] - kstart[m]);
|
||||
double deq = deqstart[m] + delta * (deqstop[m] - deqstart[m]);
|
||||
|
||||
i1 = atom->map(ids[m][0]);
|
||||
i2 = atom->map(ids[m][1]);
|
||||
|
||||
// newton_bond on: only processor owning i2 computes restraint
|
||||
// newton_bond off: only processors owning either of i1,i2 computes restraint
|
||||
|
||||
if (newton_bond) {
|
||||
if (i2 == -1 || i2 >= nlocal) return;
|
||||
if (i1 == -1) {
|
||||
char str[128];
|
||||
sprintf(str,
|
||||
"Restrain atoms %d %d missing on proc %d at step " BIGINT_FORMAT,
|
||||
ids[m][0],ids[m][1],
|
||||
comm->me,update->ntimestep);
|
||||
error->one(FLERR,str);
|
||||
}
|
||||
} else {
|
||||
if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal)) return;
|
||||
if (i1 == -1 || i2 == -1) {
|
||||
char str[128];
|
||||
sprintf(str,
|
||||
"Restrain atoms %d %d missing on proc %d at step " BIGINT_FORMAT,
|
||||
ids[m][0],ids[m][1],
|
||||
comm->me,update->ntimestep);
|
||||
error->one(FLERR,str);
|
||||
}
|
||||
}
|
||||
|
||||
delx = x[i1][0] - x[i2][0];
|
||||
dely = x[i1][1] - x[i2][1];
|
||||
delz = x[i1][2] - x[i2][2];
|
||||
domain->minimum_image(delx,dely,delz);
|
||||
|
||||
rsq = delx*delx + dely*dely + delz*delz;
|
||||
r = sqrt(rsq);
|
||||
dr = r - deq;
|
||||
rk = k * dr;
|
||||
|
||||
// force & energy
|
||||
|
||||
if (dr < 0) {
|
||||
if (r > 0.0) fbond = -2.0*rk/r;
|
||||
else fbond = 0.0;
|
||||
|
||||
elbound += rk*dr;
|
||||
energy += rk*dr;
|
||||
} else {
|
||||
fbond = 0.0;
|
||||
|
||||
elbound += 0.0;
|
||||
energy += 0.0;
|
||||
}
|
||||
|
||||
// apply force to each of 2 atoms
|
||||
|
||||
if (newton_bond || i1 < nlocal) {
|
||||
f[i1][0] += delx*fbond;
|
||||
f[i1][1] += dely*fbond;
|
||||
f[i1][2] += delz*fbond;
|
||||
}
|
||||
|
||||
if (newton_bond || i2 < nlocal) {
|
||||
f[i2][0] -= delx*fbond;
|
||||
f[i2][1] -= dely*fbond;
|
||||
f[i2][2] -= delz*fbond;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
apply harmonic angle restraints
|
||||
---------------------------------------------------------------------- */
|
||||
|
|
@ -655,9 +771,12 @@ double FixRestrain::compute_vector(int n)
|
|||
MPI_Allreduce(&ebond,&ebond_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
return ebond_all;
|
||||
} else if (n == 1) {
|
||||
MPI_Allreduce(&elbound,&elbound_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
return elbound_all;
|
||||
} else if (n == 3) {
|
||||
MPI_Allreduce(&eangle,&eangle_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
return eangle_all;
|
||||
} else if (n == 2) {
|
||||
} else if (n == 4) {
|
||||
MPI_Allreduce(&edihed,&edihed_all,1,MPI_DOUBLE,MPI_SUM,world);
|
||||
return edihed_all;
|
||||
} else {
|
||||
|
|
|
|||
|
|
@ -44,14 +44,16 @@ class FixRestrain : public Fix {
|
|||
int *rstyle;
|
||||
int *mult;
|
||||
int **ids;
|
||||
double *kstart,*kstop,*target;
|
||||
double *kstart,*kstop,*deqstart,*deqstop,*target;
|
||||
double *cos_target,*sin_target;
|
||||
double energy,energy_all;
|
||||
double ebond,ebond_all;
|
||||
double elbound,elbound_all;
|
||||
double eangle,eangle_all;
|
||||
double edihed,edihed_all;
|
||||
|
||||
void restrain_bond(int);
|
||||
void restrain_lbound(int);
|
||||
void restrain_angle(int);
|
||||
void restrain_dihedral(int);
|
||||
};
|
||||
|
|
|
|||
|
|
@ -1 +1 @@
|
|||
#define LAMMPS_VERSION "19 Mar 2020"
|
||||
#define LAMMPS_VERSION "15 Apr 2020"
|
||||
|
|
|
|||
Loading…
Reference in New Issue