diff --git a/doc/src/Commands_bond.rst b/doc/src/Commands_bond.rst index 0cb1b723b8..dd56938773 100644 --- a/doc/src/Commands_bond.rst +++ b/doc/src/Commands_bond.rst @@ -46,6 +46,7 @@ OPT. * :doc:`oxdna2/fene ` * :doc:`oxrna2/fene ` * :doc:`quartic (o) ` + * :doc:`special ` * :doc:`table (o) ` .. _angle: diff --git a/doc/src/bond_special.rst b/doc/src/bond_special.rst new file mode 100644 index 0000000000..f7dc43a1b2 --- /dev/null +++ b/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 ` 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 ` 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 ` 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 ` command as in the example above, or in +the data file or restart files read by the :doc:`read_data ` +or :doc:`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 ` 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 +` 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 ` 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 ` must be used, and 1-4 interactions must +have weights of unity or :doc:`special_bonds dihedral yes ` +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 ` 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 ` doc +page for more info. + +This bond style requires the use of a :doc:`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 ` is declared, an error will be issued. + +Related commands +"""""""""""""""" + +:doc:`bond_coeff `, :doc:`special_bonds ` + +**Default:** none diff --git a/src/USER-MISC/README b/src/USER-MISC/README index e10d4c7904..e085676711 100644 --- a/src/USER-MISC/README +++ b/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 diff --git a/src/USER-MISC/bond_special.cpp b/src/USER-MISC/bond_special.cpp new file mode 100644 index 0000000000..6c734d482e --- /dev/null +++ b/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 +#include +#include +#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; +} + + diff --git a/src/USER-MISC/bond_special.h b/src/USER-MISC/bond_special.h new file mode 100644 index 0000000000..389990b1d7 --- /dev/null +++ b/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 +#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. + +*/