add gromos bond style

doc/src/Eqs/bond_gromos.tex

@@ -0,0 +1,10 @@
+\documentclass[12pt]{article}
+\pagestyle{empty}
+
+\begin{document}
+
+$$
+   E = K (r^2 - r_0^2)^2
+$$
+
+\end{document}

doc/src/Section_commands.txt

@@ -1111,6 +1111,7 @@ KOKKOS, o = USER-OMP, t = OPT.
 "class2 (ko)"_bond_class2.html,
 "fene (iko)"_bond_fene.html,
 "fene/expand (o)"_bond_fene_expand.html,
+"gromos (o)"_bond_gromos.html,
 "harmonic (ko)"_bond_harmonic.html,
 "morse (o)"_bond_morse.html,
 "nonlinear (o)"_bond_nonlinear.html,

doc/src/bond_gromos.txt

@@ -0,0 +1,73 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+bond_style gromos command :h3
+bond_style gromos/omp command :h3
+
+[Syntax:]
+
+bond_style gromos :pre
+
+[Examples:]
+
+bond_style gromos
+bond_coeff 5 80.0 1.2 :pre
+
+[Description:]
+
+The {gromos} bond style uses the potential
+
+:c,image(Eqs/bond_gromos.jpg)
+
+where r0 is the equilibrium bond distance.  Note that the usual 1/4
+factor is included in K.
+
+The following coefficients must be defined for each bond type via the
+"bond_coeff"_bond_coeff.html command as in the example above, or in
+the data file or restart files read by the "read_data"_read_data.html
+or "read_restart"_read_restart.html commands:
+
+K (energy/distance^4)
+r0 (distance) :ul
+
+:line
+
+Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
+functionally the same as the corresponding style without the suffix.
+They have been optimized to run faster, depending on your available
+hardware, as discussed in "Section 5"_Section_accelerate.html
+of the manual.  The accelerated styles take the same arguments and
+should produce the same results, except for round-off and precision
+issues.
+
+These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
+USER-OMP and OPT packages, respectively.  They are only enabled if
+LAMMPS was built with those packages.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info.
+
+You can specify the accelerated styles explicitly in your input script
+by including their suffix, or you can use the "-suffix command-line
+switch"_Section_start.html#start_6 when you invoke LAMMPS, or you can
+use the "suffix"_suffix.html command in your input script.
+
+See "Section 5"_Section_accelerate.html of the manual for
+more instructions on how to use the accelerated styles effectively.
+
+:line
+
+[Restrictions:]
+
+This bond style can only be used if LAMMPS was built with the
+MOLECULE package.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info on packages.
+
+[Related commands:]
+
+"bond_coeff"_bond_coeff.html, "delete_bonds"_delete_bonds.html
+
+[Default:] none

doc/src/bonds.txt

@@ -8,6 +8,7 @@ Bond Styles :h1
    bond_class2
    bond_fene
    bond_fene_expand
+   bond_gromos
    bond_harmonic
    bond_harmonic_shift
    bond_harmonic_shift_cut

doc/src/lammps.book

@@ -515,7 +515,7 @@ pair_zero.html
 bond_class2.html
 bond_fene.html
 bond_fene_expand.html
-bond_oxdna.html
+bond_gromos.html
 bond_harmonic.html
 bond_harmonic_shift.html
 bond_harmonic_shift_cut.html
@@ -523,6 +523,7 @@ bond_hybrid.html
 bond_morse.html
 bond_none.html
 bond_nonlinear.html
+bond_oxdna.html
 bond_quartic.html
 bond_table.html
 bond_zero.html

src/.gitignore

@@ -185,6 +185,8 @@
 /bond_fene.h
 /bond_fene_expand.cpp
 /bond_fene_expand.h
+/bond_gromos.cpp
+/bond_gromos.h
 /bond_harmonic.cpp
 /bond_harmonic.h
 /bond_harmonic_shift.cpp

src/MOLECULE/bond_gromos.cpp

@@ -0,0 +1,210 @@
+/* ----------------------------------------------------------------------
+   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: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include "bond_gromos.h"
+#include "atom.h"
+#include "neighbor.h"
+#include "domain.h"
+#include "comm.h"
+#include "force.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+BondGromos::BondGromos(LAMMPS *lmp) : Bond(lmp)
+{
+  reinitflag = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+BondGromos::~BondGromos()
+{
+  if (allocated && !copymode) {
+    memory->destroy(setflag);
+    memory->destroy(k);
+    memory->destroy(r0);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondGromos::compute(int eflag, int vflag)
+{
+  int i1,i2,n,type;
+  double delx,dely,delz,ebond,fbond;
+
+  ebond = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = 0;
+
+  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];
+
+    delx = x[i1][0] - x[i2][0];
+    dely = x[i1][1] - x[i2][1];
+    delz = x[i1][2] - x[i2][2];
+
+    const double rsq = delx*delx + dely*dely + delz*delz;
+    const double dr = rsq - r0[type]*r0[type];
+    const double kdr = k[type]*dr;
+
+    // force & energy
+
+    fbond = -4.0 * kdr;
+    if (eflag) ebond = kdr;
+
+    // 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 BondGromos::allocate()
+{
+  allocated = 1;
+  int n = atom->nbondtypes;
+
+  memory->create(k,n+1,"bond:k");
+  memory->create(r0,n+1,"bond:r0");
+
+  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 BondGromos::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 k_one = force->numeric(FLERR,arg[1]);
+  double r0_one = force->numeric(FLERR,arg[2]);
+
+  int count = 0;
+  for (int i = ilo; i <= ihi; i++) {
+    k[i] = k_one;
+    r0[i] = r0_one;
+    setflag[i] = 1;
+    count++;
+  }
+
+  if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   return an equilbrium bond length
+------------------------------------------------------------------------- */
+
+double BondGromos::equilibrium_distance(int i)
+{
+  return r0[i];
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 writes out coeffs to restart file
+------------------------------------------------------------------------- */
+
+void BondGromos::write_restart(FILE *fp)
+{
+  fwrite(&k[1],sizeof(double),atom->nbondtypes,fp);
+  fwrite(&r0[1],sizeof(double),atom->nbondtypes,fp);
+}
+
+/* ----------------------------------------------------------------------
+   proc 0 reads coeffs from restart file, bcasts them
+------------------------------------------------------------------------- */
+
+void BondGromos::read_restart(FILE *fp)
+{
+  allocate();
+
+  if (comm->me == 0) {
+    fread(&k[1],sizeof(double),atom->nbondtypes,fp);
+    fread(&r0[1],sizeof(double),atom->nbondtypes,fp);
+  }
+  MPI_Bcast(&k[1],atom->nbondtypes,MPI_DOUBLE,0,world);
+  MPI_Bcast(&r0[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 BondGromos::write_data(FILE *fp)
+{
+  for (int i = 1; i <= atom->nbondtypes; i++)
+    fprintf(fp,"%d %g %g\n",i,k[i],r0[i]);
+}
+
+/* ---------------------------------------------------------------------- */
+
+double BondGromos::single(int type, double rsq, int i, int j,
+                        double &fforce)
+{
+  double dr = rsq - r0[type]*r0[type];
+  fforce = -4.0*k[type] * dr;
+  return k[type]*dr;
+}
+
+/* ----------------------------------------------------------------------
+    Return ptr to internal members upon request.
+------------------------------------------------------------------------ */
+void *BondGromos::extract( char *str, int &dim )
+{
+  dim = 1;
+  if( strcmp(str,"kappa")==0) return (void*) k;
+  if( strcmp(str,"r0")==0) return (void*) r0;
+  return NULL;
+}

src/MOLECULE/bond_gromos.h

@@ -0,0 +1,58 @@
+/* -*- 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 BOND_CLASS
+
+BondStyle(gromos,BondGromos)
+
+#else
+
+#ifndef LMP_BOND_GROMOS_H
+#define LMP_BOND_GROMOS_H
+
+#include <stdio.h>
+#include "bond.h"
+
+namespace LAMMPS_NS {
+
+class BondGromos : public Bond {
+ public:
+  BondGromos(class LAMMPS *);
+  virtual ~BondGromos();
+  virtual 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 &);
+  virtual void *extract(char *, int &);
+
+ protected:
+  double *k,*r0;
+
+  virtual void allocate();
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Incorrect args for bond coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+*/

src/USER-OMP/bond_gromos_omp.cpp

@@ -0,0 +1,129 @@
+/* ----------------------------------------------------------------------
+   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: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#include "bond_gromos_omp.h"
+#include "atom.h"
+#include "comm.h"
+#include "force.h"
+#include "neighbor.h"
+#include "domain.h"
+
+#include <math.h>
+
+#include "suffix.h"
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+BondGromosOMP::BondGromosOMP(class LAMMPS *lmp)
+  : BondGromos(lmp), ThrOMP(lmp,THR_BOND)
+{
+  suffix_flag |= Suffix::OMP;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void BondGromosOMP::compute(int eflag, int vflag)
+{
+
+  if (eflag || vflag) {
+    ev_setup(eflag,vflag);
+  } else evflag = 0;
+
+  const int nall = atom->nlocal + atom->nghost;
+  const int nthreads = comm->nthreads;
+  const int inum = neighbor->nbondlist;
+
+#if defined(_OPENMP)
+#pragma omp parallel default(none) shared(eflag,vflag)
+#endif
+  {
+    int ifrom, ito, tid;
+
+    loop_setup_thr(ifrom, ito, tid, inum, nthreads);
+    ThrData *thr = fix->get_thr(tid);
+    thr->timer(Timer::START);
+    ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);
+
+    if (inum > 0) {
+      if (evflag) {
+        if (eflag) {
+          if (force->newton_bond) eval<1,1,1>(ifrom, ito, thr);
+          else eval<1,1,0>(ifrom, ito, thr);
+        } else {
+          if (force->newton_bond) eval<1,0,1>(ifrom, ito, thr);
+          else eval<1,0,0>(ifrom, ito, thr);
+        }
+      } else {
+        if (force->newton_bond) eval<0,0,1>(ifrom, ito, thr);
+        else eval<0,0,0>(ifrom, ito, thr);
+      }
+    }
+    thr->timer(Timer::BOND);
+    reduce_thr(this, eflag, vflag, thr);
+  } // end of omp parallel region
+}
+
+template <int EVFLAG, int EFLAG, int NEWTON_BOND>
+void BondGromosOMP::eval(int nfrom, int nto, ThrData * const thr)
+{
+  int i1,i2,n,type;
+  double delx,dely,delz,ebond,fbond;
+
+  const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
+  dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0];
+  const int3_t * _noalias const bondlist = (int3_t *) neighbor->bondlist[0];
+  const int nlocal = atom->nlocal;
+  ebond = 0.0;
+
+  for (n = nfrom; n < nto; n++) {
+    i1 = bondlist[n].a;
+    i2 = bondlist[n].b;
+    type = bondlist[n].t;
+
+    delx = x[i1].x - x[i2].x;
+    dely = x[i1].y - x[i2].y;
+    delz = x[i1].z - x[i2].z;
+
+    const double rsq = delx*delx + dely*dely + delz*delz;
+    const double dr = rsq - r0[type]*r0[type];
+    const double kdr = k[type]*dr;
+
+    // force & energy
+
+    fbond = -4.0 * kdr;
+
+    if (EFLAG) ebond = kdr;
+
+    // apply force to each of 2 atoms
+
+    if (NEWTON_BOND || i1 < nlocal) {
+      f[i1].x += delx*fbond;
+      f[i1].y += dely*fbond;
+      f[i1].z += delz*fbond;
+    }
+
+    if (NEWTON_BOND || i2 < nlocal) {
+      f[i2].x -= delx*fbond;
+      f[i2].y -= dely*fbond;
+      f[i2].z -= delz*fbond;
+    }
+
+    if (EVFLAG) ev_tally_thr(this,i1,i2,nlocal,NEWTON_BOND,
+                             ebond,fbond,delx,dely,delz,thr);
+  }
+}

src/USER-OMP/bond_gromos_omp.h

@@ -0,0 +1,46 @@
+/* -*- 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: Axel Kohlmeyer (Temple U)
+------------------------------------------------------------------------- */
+
+#ifdef BOND_CLASS
+
+BondStyle(gromos/omp,BondGromosOMP)
+
+#else
+
+#ifndef LMP_BOND_GROMOS_OMP_H
+#define LMP_BOND_GROMOS_OMP_H
+
+#include "bond_gromos.h"
+#include "thr_omp.h"
+
+namespace LAMMPS_NS {
+
+class BondGromosOMP : public BondGromos, public ThrOMP {
+
+ public:
+  BondGromosOMP(class LAMMPS *lmp);
+  virtual void compute(int, int);
+
+ private:
+  template <int EVFLAG, int EFLAG, int NEWTON_BOND>
+  void eval(int ifrom, int ito, ThrData * const thr);
+};
+
+}
+
+#endif
+#endif