Merge pull request #1865 from danicholson/bond-style-special

Add bond style special

doc/src/Commands_bond.rst

@@ -46,6 +46,7 @@ OPT.
    * :doc:`oxdna2/fene <bond_oxdna>`
    * :doc:`oxrna2/fene <bond_oxdna>`
    * :doc:`quartic (o) <bond_quartic>`
+   * :doc:`special <bond_special>`
    * :doc:`table (o) <bond_table>`
 
 .. _angle:

doc/src/bond_special.rst

@@ -0,0 +1,106 @@
+.. index:: bond_style special
+
+bond_style special command
+=================================
+
+Syntax
+""""""
+
+
+.. code-block:: LAMMPS
+
+   bond_style special
+
+Examples
+""""""""
+
+
+.. code-block:: LAMMPS
+
+   bond_style special
+   bond_coeff 0.5 0.5
+
+Description
+"""""""""""
+
+The *special* bond style can be used to create conceptual bonds which
+effectively impose weightings on the pairwise Lennard Jones and/or
+Coulombic interactions between selected pairs of particles in the
+system.  The form of the pairwise interaction will be whatever is
+computed by the :doc:`pair_style <pair_style>` command defined for the
+system; this command defines the weightings for its two terms.
+
+This command can thus be useful to apply weightings that cannot be
+handled by the :doc:`special_bonds <special_bonds>` command, such as
+on 1-5 or 1-6 interactions.  Or it can be used to add pairwise forces
+between one or more pairs of atoms that otherwise would not be include
+in the :doc:`pair_style <pair_style>` computation.
+
+The potential for this bond style has the form
+
+.. math::
+
+   E =  w_{LJ} E_{LJ} + w_{Coul} E_{Coul}
+
+The following coefficients must be defined for each bond type via the
+:doc:`bond_coeff <bond_coeff>` command as in the example above, or in
+the data file or restart files read by the :doc:`read_data <read_data>`
+or :doc:`read_restart <read_restart>` commands:
+
+* :math:`w_{LJ}` weight (0.0 to 1.0) on pairwise Lennard-Jones interactions
+
+* :math:`w_{Coul}` weight (0.0 to 1.0) on pairwise Coulombic interactions
+
+----------
+
+Normally this bond style should be used in conjunction with one (or
+more) other bond styles which compute forces between atoms directly
+bonded to each other in a molecule.  This means the :doc:`bond_style
+hybrid <bond_hybrid>` command should be used with bond_style special
+as one of its sub-styles.
+
+Note that the same as for any other bond style, pairs of bonded atoms
+must be enumerated in the data file read by the :doc:`read_data
+<read_data>` command.  Thus if this command is used to weight all 1-5
+interactions in the system, all the 1-5 pairs of atoms must be listed
+in the "Bonds" section of the data file.
+
+This bond style imposes strict requirements on settings made with the
+:doc:`special_bonds <special_bonds>` command.  These requirements
+ensure that the new bonds created by this style do not create spurious
+1-2, 1-3, or 1-4 interactions within the molecular topology.
+
+Specifically 1-2 interactions must have weights of zero, 1-3
+interactions must either have weights of unity or :doc:`special_bonds
+angle yes <special_bonds>` must be used, and 1-4 interactions must
+have weights of unity or :doc:`special_bonds dihedral yes <special_bonds>` 
+must be used.
+
+If this command is used to create bonded interactions between
+particles that are further apart than usual (e.g. 1-5 or 1-6
+interactions), this style may require an increase in the communication
+cutoff via the :doc:`comm_modify cutoff <comm_modify>` command.  If
+LAMMPS cannot find a partner atom in a bond, an error will be issued.
+
+Restrictions
+""""""""""""
+
+This bond style can only be used if LAMMPS was built with the
+USER-MISC package.  See the :doc:`Build package <Build_package>` doc
+page for more info.
+
+This bond style requires the use of a :doc:`pair_style <pair_style>` which
+computes a pairwise additive interaction and provides the ability to
+compute interactions for individual pairs of atoms.  Manybody potentials
+are not compatible in general, but also some other pair styles are missing
+the required functionality and thus will cause an error.
+
+This command is not compatible with long-range Coulombic interactions. If a 
+`kspace_style <kspace_style>` is declared, an error will be issued.
+
+Related commands
+""""""""""""""""
+
+:doc:`bond_coeff <bond_coeff>`, :doc:`special_bonds <special_bonds>`
+
+**Default:** none

src/USER-MISC/README

@@ -26,6 +26,7 @@ angle_style dipole, Mario Orsi, orsimario at gmail.com, 10 Jan 12
 angle_style quartic, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12
 bond_style harmonic/shift, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
 bond_style harmonic/shift/cut, Carsten Svaneborg, science at zqex.dk, 8 Aug 11
+bond_style special, David Nicholson, davidanich at gmail.com, 31 Jan 2020
 compute ackland/atom, Gerolf Ziegenhain, gerolf at ziegenhain.com, 4 Oct 2007
 compute basal/atom, Christopher Barrett, cdb333 at cavs.msstate.edu, 3 Mar 2013
 compute cnp/atom, Paulo Branicio (USC), branicio at usc.edu, 31 May 2017

src/USER-MISC/bond_special.cpp

@@ -0,0 +1,222 @@
+/* ----------------------------------------------------------------------
+   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 Author: David Nicholson (MIT)
+------------------------------------------------------------------------- */
+
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+#include "bond_special.h"
+#include "atom.h"
+#include "neighbor.h"
+#include "domain.h"
+#include "comm.h"
+#include "force.h"
+#include "pair.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+BondSpecial::BondSpecial(LAMMPS *lmp) : Bond(lmp) {}
+
+/* ---------------------------------------------------------------------- */
+
+BondSpecial::~BondSpecial()
+{
+  if (allocated) {
+    memory->destroy(setflag);
+    memory->destroy(factor_lj);
+    memory->destroy(factor_coul);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondSpecial::init_style()
+{
+  if (force->pair == NULL) error->all(FLERR,"No pair style defined");
+  else if ((force->pair->single_enable == 0) || force->pair->manybody_flag)
+    error->all(FLERR,"Pair style does not support bond style special");
+
+  if (force->special_lj[1] != 0.0 || force->special_coul[1] != 0)
+    error->all(FLERR,"Invalid 1-2 setting for bond style special.");
+
+  if (force->special_angle != 1 && (force->special_lj[2] != 1.0 ||
+                                    force->special_coul[2] != 1.0))
+    error->all(FLERR,"Invalid 1-3 setting for bond style special.");
+
+  if (force->special_dihedral != 1 && (force->special_lj[3] != 1.0 ||
+                                       force->special_coul[3] != 1.0))
+    error->all(FLERR,"Invalid 1-4 setting for bond style special.");
+
+  if (force->kspace != NULL)
+    error->all(FLERR,"Bond style special is not compatible with long range "
+        "Coulombic interactions");
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondSpecial::compute(int eflag, int vflag)
+{
+  int i1,i2,i1type,i2type,n,type;
+  double delx,dely,delz,ebond,fbond;
+  double rsq;
+
+  ev_init(eflag,vflag);
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int **bondlist = neighbor->bondlist;
+  int nbondlist = neighbor->nbondlist;
+  int nlocal = atom->nlocal;
+  int newton_bond = force->newton_bond;
+
+  for (n = 0; n < nbondlist; n++) {
+    i1 = bondlist[n][0];
+    i2 = bondlist[n][1];
+    type = bondlist[n][2];
+
+    i1type = atom->type[i1];
+    i2type = atom->type[i2];
+
+    delx = x[i1][0] - x[i2][0];
+    dely = x[i1][1] - x[i2][1];
+    delz = x[i1][2] - x[i2][2];
+
+    rsq = delx*delx + dely*dely + delz*delz;
+
+    ebond = force->pair->single(i1,i2,i1type,i2type,rsq,factor_lj[type],
+        factor_coul[type],fbond);
+
+    // 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;
+    }
+
+    if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondSpecial::allocate()
+{
+  allocated = 1;
+  int n = atom->nbondtypes;
+
+  memory->create(factor_lj,n+1,"bond:factor_lj");
+  memory->create(factor_coul,n+1,"bond:factor_coul");
+
+  memory->create(setflag,n+1,"bond:setflag");
+  for (int i = 1; i <= n; i++) setflag[i] = 0;
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more types
+------------------------------------------------------------------------- */
+
+void BondSpecial::coeff(int narg, char **arg)
+{
+  if (narg != 3) error->all(FLERR,"Incorrect args for bond coefficients");
+  if (!allocated) allocate();
+
+  int ilo,ihi;
+  force->bounds(FLERR,arg[0],atom->nbondtypes,ilo,ihi);
+
+  double factor_lj_one = force->numeric(FLERR,arg[1]);
+  double factor_coul_one = force->numeric(FLERR,arg[2]);
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    factor_lj[i] = factor_lj_one;
+    factor_coul[i] = factor_coul_one;
+    setflag[i] = 1;
+    count++;
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   return an equilbrium bond length
+------------------------------------------------------------------------- */
+
+double BondSpecial::equilibrium_distance(int i)
+{
+  return 0.;
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes out coeffs to restart file
+------------------------------------------------------------------------- */
+
+void BondSpecial::write_restart(FILE *fp)
+{
+  fwrite(&factor_lj[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&factor_coul[1],sizeof(double),atom->nbondtypes,fp);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads coeffs from restart file, bcasts them
+------------------------------------------------------------------------- */
+
+void BondSpecial::read_restart(FILE *fp)
+{
+  allocate();
+
+  if (comm->me == 0) {
+    fread(&factor_lj[1],sizeof(double),atom->nbondtypes,fp);
+    fread(&factor_coul[1],sizeof(double),atom->nbondtypes,fp);
+  }
+  MPI_Bcast(&factor_lj[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&factor_coul[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+
+  for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void BondSpecial::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->nbondtypes; i++)
+    fprintf(fp,"%d %g %g\n",i,factor_lj[i],factor_coul[i]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double BondSpecial::single(int type, double rsq, int i, int j,
+                        double &fforce)
+{
+  int itype = atom->type[i];
+  int jtype = atom->type[j];
+  double ebond;
+  ebond = force->pair->single(i,j,itype,jtype,rsq,factor_lj[type],
+      factor_coul[type],fforce);
+  return ebond;
+}
+
+

src/USER-MISC/bond_special.h

@@ -0,0 +1,78 @@
+/* -*- 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.
+
+   Contributing Author: David Nicholson (MIT)
+------------------------------------------------------------------------- */
+
+#ifdef BOND_CLASS
+
+BondStyle(special,BondSpecial)
+
+#else
+
+#ifndef LMP_BOND_SPECIAL_H
+#define LMP_BOND_SPECIAL_H
+
+#include <cstdio>
+#include "bond.h"
+
+namespace LAMMPS_NS {
+
+class BondSpecial: public Bond {
+ public:
+  BondSpecial(class LAMMPS *);
+  virtual ~BondSpecial();
+  void init_style();
+  void compute(int, int);
+  void coeff(int, char **);
+  double equilibrium_distance(int);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_data(FILE *);
+  double single(int, double, int, int, double &);
+
+ protected:
+  double *factor_lj,*factor_coul;
+
+  void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Incorrect args for bond coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+E: Invalid 1-2 setting for bond style special.
+
+Bond style special must be used with zero factors for 1-2 special bonds.
+
+E: Invalid 1-3 setting for bond style special.
+
+Bond style special must be used with 1.0 factors for 1-3 special bonds or the
+angle keyword set to yes.
+
+E: Invalid 1-4 setting for bond style special.
+
+Bond style special must be used with 1.0 factors for 1-4 special bonds or the
+dihedral keyword set to yes.
+
+E: Bond style special is not compatible with long range Coulombic interactions.
+
+Self-explanatory.
+
+*/