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

Merge pull request #1957 from david-castillo/restrain_lbound 

add moving bond restraints and a lower bound only bond restraint
This commit is contained in:
Axel Kohlmeyer 2020-04-15 14:45:07 -04:00 committed by GitHub
parent 6bac08322c c5f6b5c687
commit 4c870be433
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 164 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
doc/src/fix_restrain.rst
View File

@ -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

2
src/KOKKOS/min_kokkos.cpp
View File

@ -636,7 +636,7 @@ double MinKokkos::fnorm_inf()
double MinKokkos::fnorm_max()
{
double local_norm_max = 0.0;
{
// local variables for lambda capture

133
src/fix_restrain.cpp
View File

@ -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 {

4
src/fix_restrain.h
View File

@ -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);
};