# Install/unInstall package files in LAMMPS

# mode = 0/1/2 for uninstall/install/update

mode=$1

# arg1 = file, arg2 = file it depends on

# enforce using portable C locale

LC_ALL=C

export LC_ALL

action () {

  if (test $mode = 0) then

    rm -f ../$1

  elif (! cmp -s $1 ../$1) then

    if (test -z "$2" || test -e ../$2) then

      cp $1 ..

      if (test $mode = 2) then

        echo "  updating src/$1"

      fi

    fi

  elif (test -n "$2") then

    if (test ! -e ../$2) then

      rm -f ../$1

    fi

  fi

}

# list of files with dependcies

action  bond_oxdna_fene.cpp bond_fene.h

action  bond_oxdna_fene.h bond_fene.h

action  fix_nve_dotc_langevin.cpp atom_vec_ellipsoid.h

action  fix_nve_dotc_langevin.h atom_vec_ellipsoid.h

action  fix_nve_dot.cpp atom_vec_ellipsoid.h

action  fix_nve_dot.h atom_vec_ellipsoid.h

action  mf_oxdna.h atom_vec_ellipsoid.h

action  pair_oxdna_coaxstk.cpp atom_vec_ellipsoid.h

action  pair_oxdna_coaxstk.h atom_vec_ellipsoid.h

action  pair_oxdna_excv.cpp atom_vec_ellipsoid.h

action  pair_oxdna_excv.h atom_vec_ellipsoid.h

action  pair_oxdna_hbond.cpp atom_vec_ellipsoid.h

action  pair_oxdna_hbond.h atom_vec_ellipsoid.h

action  pair_oxdna_stk.cpp atom_vec_ellipsoid.h

action  pair_oxdna_stk.h atom_vec_ellipsoid.h

action  pair_oxdna_xstk.cpp atom_vec_ellipsoid.h

action  pair_oxdna_xstk.h atom_vec_ellipsoid.h

This package contains a LAMMPS implementation of coarse-grained

models of DNA, which can be used to model sequence-specific

DNA strands.

See the doc pages and [1,2] for the individual bond and pair styles. 

The packages contains also a new Langevin-type rigid-body integrator,

which has also its own doc page and is explained in [3].

[1] T. Ouldridge, A. Louis, J. Doye, "Structural, mechanical, 

and thermodynamic properties of a coarse-grained DNA model",

J. Chem. Phys. 134, 085101 (2011).

[2] T.E. Ouldridge, Coarse-grained modelling of DNA and DNA 

self-assembly, DPhil. University of Oxford (2011).

[3] R. Davidchack, T. Ouldridge, M. Tretyakov, "New Langevin and 

gradient thermostats for rigid body dynamics", J. Chem. Phys. 142, 

144114 (2015).

Example input and data files can be found in 

/examples/USER/cgdna/examples/duplex1/ and /duplex2/.

A simple python setup tool which creates single straight or helical DNA 

strands as well as DNA duplexes and arrays of duplexes can be found in 

/examples/USER/cgdna/util/.

A technical report with more information on the model, the structure 

of the input and data file, the setup tool and the performance of the 

LAMMPS-implementation of oxDNA can be found in

/doc/src/PDF/USER-CGDNA-overview.pdf.

IMPORTANT NOTE: This package can only be used if LAMMPS is compiled

with the MOLECULE and ASPHERE packages.  These should be included 

in the LAMMPS build by typing "make yes-asphere yes-molecule" prior 

to the usual compilation (see the "Including/excluding packages" 

section of the LAMMPS manual).

The creator of this package is:

Dr Oliver Henrich

University of Edinburgh, UK

ohenrich@ph.ed.ac.uk

o.henrich@epcc.ed.ac.uk

--------------------------------------------------------------------------

Bond styles provided by this package:

bond_oxdna_fene.cpp:  backbone connectivity, a modified FENE potential

Pair styles provided by this package:

pair_oxdna_excv.cpp:  excluded volume interaction between the nucleotides

pair_oxdna_stk.cpp:  stacking interaction between consecutive nucleotides

                     on the same strand

pair_oxdna_hbond.cpp:  hydrogen-bonding interaction between complementary

                       nucleotides on different strands, e.g. A-T and C-G

pair_oxdna_xstk.cpp:  cross-stacking interaction between nucleotides

pair_oxdna_coaxstk.cpp:  coaxial stacking interaction between nucleotides

Fixes provided by this package:

fix_nve_dotc_langevin.cpp:  fix for Langevin-type rigid body integrator "C"

                            in above Ref. [3] 

fix_nve_dot.cpp:  NVE-type rigid body integrator without noise

/* ----------------------------------------------------------------------

   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: Oliver Henrich (EPCC, University of Edinburgh)

------------------------------------------------------------------------- */

#include <math.h>

#include <stdlib.h>

#include "bond_oxdna_fene.h"

#include "atom.h"

#include "neighbor.h"

#include "domain.h"

#include "comm.h"

#include "update.h"

#include "force.h"

#include "memory.h"

#include "error.h"

#include "atom_vec_ellipsoid.h"

#include "math_extra.h"

using namespace LAMMPS_NS;

/* ---------------------------------------------------------------------- */

BondOxdnaFene::BondOxdnaFene(LAMMPS *lmp) : Bond(lmp)

{

}

/* ---------------------------------------------------------------------- */

BondOxdnaFene::~BondOxdnaFene()

{

  if (allocated) {

    memory->destroy(setflag);

    memory->destroy(k);

    memory->destroy(Delta);

    memory->destroy(r0);

  }

}

/* ----------------------------------------------------------------------

   compute function for oxDNA FENE-bond interaction

   s=sugar-phosphate backbone site, b=base site, st=stacking site

------------------------------------------------------------------------- */

void BondOxdnaFene::compute(int eflag, int vflag)

{

  int a,b,in,type;

  double delf[3],delta[3],deltb[3]; // force, torque increment;;

  double delr[3],ebond,fbond;

  double rsq,Deltasq,rlogarg;

  double r,rr0,rr0sq;

  // distances COM-backbone site

  double d_cs=-0.24;

  // vectors COM-backbone site in lab frame

  double ra_cs[3],rb_cs[3];

  double *qa,ax[3],ay[3],az[3];

  double *qb,bx[3],by[3],bz[3];

  double **x = atom->x;

  double **f = atom->f;

  double **torque = atom->torque;

  AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");

  AtomVecEllipsoid::Bonus *bonus = avec->bonus;

  int **bondlist = neighbor->bondlist;

  int nbondlist = neighbor->nbondlist;

  int nlocal = atom->nlocal;

  int newton_bond = force->newton_bond;

  ebond = 0.0;

  if (eflag || vflag) ev_setup(eflag,vflag);

  else evflag = 0;

  // loop over FENE bonds

  for (in = 0; in < nbondlist; in++) {

    a = bondlist[in][1];

    b = bondlist[in][0];

    type = bondlist[in][2];

    qa=bonus[a].quat;

    MathExtra::q_to_exyz(qa,ax,ay,az);

    qb=bonus[b].quat;

    MathExtra::q_to_exyz(qb,bx,by,bz);

    // vector COM-backbone site a and b

    ra_cs[0] = d_cs*ax[0];

    ra_cs[1] = d_cs*ax[1];

    ra_cs[2] = d_cs*ax[2];

    rb_cs[0] = d_cs*bx[0];

    rb_cs[1] = d_cs*bx[1];

    rb_cs[2] = d_cs*bx[2];

    // vector backbone site b to a

    delr[0] = x[a][0] + ra_cs[0] - x[b][0] - rb_cs[0];

    delr[1] = x[a][1] + ra_cs[1] - x[b][1] - rb_cs[1];

    delr[2] = x[a][2] + ra_cs[2] - x[b][2] - rb_cs[2];

    rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];

    r = sqrt(rsq);

    rr0 = r - r0[type];

    rr0sq = rr0*rr0;

    Deltasq = Delta[type] * Delta[type];

    rlogarg = 1.0 - rr0sq/Deltasq;

    // if r -> Delta, then rlogarg < 0.0 which is an error

    // issue a warning and reset rlogarg = epsilon

    // if r > 2*Delta something serious is wrong, abort

    if (rlogarg < 0.1) {

      char str[128];

      sprintf(str,"FENE bond too long: " BIGINT_FORMAT " "

              TAGINT_FORMAT " " TAGINT_FORMAT " %g",

              update->ntimestep,atom->tag[a],atom->tag[b],r);

      error->warning(FLERR,str,0);

      if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond");

    }

    fbond = -k[type]*rr0/rlogarg/Deltasq/r;

    delf[0] = delr[0]*fbond;

    delf[1] = delr[1]*fbond;

    delf[2] = delr[2]*fbond;

    // energy

    if (eflag) {

      ebond = -0.5 * k[type]*log(rlogarg);

    }

    // apply force and torque to each of 2 atoms

    if (newton_bond || a < nlocal) {

      f[a][0] += delf[0];

      f[a][1] += delf[1];

      f[a][2] += delf[2];

      MathExtra::cross3(ra_cs,delf,delta);

      torque[a][0] += delta[0];

      torque[a][1] += delta[1];

      torque[a][2] += delta[2];

    }

    if (newton_bond || b < nlocal) {

      f[b][0] -= delf[0];

      f[b][1] -= delf[1];

      f[b][2] -= delf[2];

      MathExtra::cross3(rb_cs,delf,deltb);

      torque[b][0] -= deltb[0];

      torque[b][1] -= deltb[1];

      torque[b][2] -= deltb[2];

    }

    // increment energy and virial

    if (evflag) ev_tally(a,b,nlocal,newton_bond,ebond,fbond,delr[0],delr[1],delr[2]);

  }

}

/* ---------------------------------------------------------------------- */

void BondOxdnaFene::allocate()

{

  allocated = 1;

  int n = atom->nbondtypes;

  memory->create(k,n+1,"bond:k");

  memory->create(Delta,n+1,"bond:Delta");

  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 type

------------------------------------------------------------------------- */

void BondOxdnaFene::coeff(int narg, char **arg)

{

  if (narg != 4) error->all(FLERR,"Incorrect args for bond coefficients in oxdna_fene");

  if (!allocated) allocate();

  int ilo,ihi;

  force->bounds(FLERR,arg[0],atom->nbondtypes,ilo,ihi);

  double k_one = force->numeric(FLERR,arg[1]);

  double Delta_one = force->numeric(FLERR,arg[2]);

  double r0_one = force->numeric(FLERR,arg[3]);

  int count = 0;

  for (int i = ilo; i <= ihi; i++) {

    k[i] = k_one;

    Delta[i] = Delta_one;

    r0[i] = r0_one;

    setflag[i] = 1;

    count++;

  }

  if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients in oxdna_fene");

}

/* ----------------------------------------------------------------------

   set special_bond settings and check if valid

------------------------------------------------------------------------- */

void BondOxdnaFene::init_style()

{

  /* special bonds have to be lj = 0 1 1 and coul = 1 1 1 to exclude

     the ss excluded volume interaction between nearest neighbours   */

  force->special_lj[1] = 0.0;

  force->special_lj[2] = 1.0;

  force->special_lj[3] = 1.0;

  force->special_coul[1] = 1.0;

  force->special_coul[2] = 1.0;

  force->special_coul[3] = 1.0;

  fprintf(screen,"Finding 1-2 1-3 1-4 neighbors ...\n"

                 " Special bond factors lj:   %-10g %-10g %-10g\n"

                 " Special bond factors coul: %-10g %-10g %-10g\n",

                 force->special_lj[1],force->special_lj[2],force->special_lj[3],

                 force->special_coul[1],force->special_coul[2],force->special_coul[3]);

  if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 ||

      force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0)

  {

    if (comm->me == 0)

      error->warning(FLERR,"Use special bonds lj = 0,1,1 and coul = 1,1,1 with bond style oxdna_fene");

  }

}

/* ---------------------------------------------------------------------- */

double BondOxdnaFene::equilibrium_distance(int i)

{

  return r0[i];

}

/* ----------------------------------------------------------------------

   proc 0 writes to restart file

------------------------------------------------------------------------- */

void BondOxdnaFene::write_restart(FILE *fp)

{

  fwrite(&k[1],sizeof(double),atom->nbondtypes,fp);

  fwrite(&Delta[1],sizeof(double),atom->nbondtypes,fp);

  fwrite(&r0[1],sizeof(double),atom->nbondtypes,fp);

}

/* ----------------------------------------------------------------------

   proc 0 reads from restart file, bcasts

------------------------------------------------------------------------- */

void BondOxdnaFene::read_restart(FILE *fp)

{

  allocate();

  if (comm->me == 0) {

    fread(&k[1],sizeof(double),atom->nbondtypes,fp);

    fread(&Delta[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(&Delta[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 BondOxdnaFene::write_data(FILE *fp)

{

  for (int i = 1; i <= atom->nbondtypes; i++)

    fprintf(fp,"%d %g %g %g\n",i,k[i],r0[i],Delta[i]);

}

/* ---------------------------------------------------------------------- */

double BondOxdnaFene::single(int type, double rsq, int i, int j,

                        double &fforce)

{

  double r = sqrt(rsq);

  double rr0 = r - r0[type];

  double rr0sq = rr0*rr0;

  double Deltasq = Delta[type] * Delta[type];

  double rlogarg = 1.0 - rr0sq/Deltasq;

  // if r -> Delta, then rlogarg < 0.0 which is an error

  // issue a warning and reset rlogarg = epsilon

  // if r > 2*Delta something serious is wrong, abort

  if (rlogarg < 0.1) {

    char str[128];

    sprintf(str,"FENE bond too long: " BIGINT_FORMAT " %g",

            update->ntimestep,sqrt(rsq));

    error->warning(FLERR,str,0);

    if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond");

  }

  double eng = -0.5 * k[type]*log(rlogarg);

  fforce = -k[type]*rr0/rlogarg/Deltasq/r;

  return eng;

}

/* -*- 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: Oliver Henrich (EPCC, University of Edinburgh)

------------------------------------------------------------------------- */

#ifdef BOND_CLASS

BondStyle(oxdna_fene,BondOxdnaFene)

#else

#ifndef LMP_BOND_OXDNA_FENE_H

#define LMP_BOND_OXDNA_FENE_H

#include "bond.h"

namespace LAMMPS_NS {

class BondOxdnaFene : public Bond {

 public:

  BondOxdnaFene(class LAMMPS *);

  virtual ~BondOxdnaFene();

  virtual void compute(int, int);

  void coeff(int, char **);

  void init_style();

  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 *k,*Delta,*r0; // FENE

  void allocate();

};

}

#endif

#endif

/* ERROR/WARNING messages:

W: FENE bond too long: %ld %d %d %g

A FENE bond has stretched dangerously far.  It's interaction strength

will be truncated to attempt to prevent the bond from blowing up.

E: Bad FENE bond

Two atoms in a FENE bond have become so far apart that the bond cannot

be computed.

E: Incorrect args for bond coefficients

Self-explanatory.  Check the input script or data file.

W: Use special bonds = 0,1,1 with bond style oxdna

Most FENE models need this setting for the special_bonds command.

W: FENE bond too long: %ld %g

A FENE bond has stretched dangerously far.  It's interaction strength

will be truncated to attempt to prevent the bond from blowing up.

*/

/* ----------------------------------------------------------------------

   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: Oliver Henrich (EPCC, University of Edinburgh)

------------------------------------------------------------------------- */

#include <math.h>

#include <stdio.h>

#include <string.h>

#include "fix_nve_dot.h"

#include "math_extra.h"

#include "atom.h"

#include "atom_vec_ellipsoid.h"

#include "force.h"

#include "update.h"

#include "memory.h"

#include "error.h"

using namespace LAMMPS_NS;

using namespace FixConst;

using namespace MathExtra;

#define INERTIA 0.2          // moment of inertia prefactor for ellipsoid

/* ---------------------------------------------------------------------- */

FixNVEDot::FixNVEDot(LAMMPS *lmp, int narg, char **arg) :

  FixNVE(lmp, narg, arg) {}

/* ---------------------------------------------------------------------- */

void FixNVEDot::init()

{

  avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");

  if (!avec)

    error->all(FLERR,"Compute nve/dot requires atom style ellipsoid");

  // check that all particles are finite-size ellipsoids

  // no point particles allowed, spherical is OK

  int *ellipsoid = atom->ellipsoid;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  for (int i = 0; i < nlocal; i++)

    if (mask[i] & groupbit)

      if (ellipsoid[i] < 0)

        error->one(FLERR,"Fix nve/dot requires extended particles");

  FixNVE::init();

}

/* ---------------------------------------------------------------------- */

void FixNVEDot::initial_integrate(int vflag)

{

  double *shape,*quat;

  double fquat[4],conjqm[4],inertia[3];

  AtomVecEllipsoid::Bonus *bonus = avec->bonus;

  int *ellipsoid = atom->ellipsoid;

  double **x = atom->x;

  double **v = atom->v;

  double **f = atom->f;

  double **angmom = atom->angmom;

  double **torque = atom->torque;

  double *rmass = atom->rmass;

  int *mask = atom->mask;

  int nlocal = atom->nlocal;

  if (igroup == atom->firstgroup) nlocal = atom->nfirst;

  // set timestep here since dt may have changed or come via rRESPA

  dt = update->dt;

  dthlf = 0.5 * dt;

  for (int i = 0; i < nlocal; i++)

    if (mask[i] & groupbit) {

      dthlfm = dthlf / rmass[i];

      quat = bonus[ellipsoid[i]].quat;

      shape = bonus[ellipsoid[i]].shape;

      // update momentum by 1/2 step

      v[i][0] += dthlfm * f[i][0];

      v[i][1] += dthlfm * f[i][1];

      v[i][2] += dthlfm * f[i][2];

      // update position by full step

      x[i][0] += dt * v[i][0];

      x[i][1] += dt * v[i][1];

      x[i][2] += dt * v[i][2];

      // convert angular momentum and torque in space frame into

      // quaternion 4-momentum and 1/2 of 4-torque in body frame

      vec3_to_vec4(quat,angmom[i],conjqm);

      conjqm[0] *= 2.0;

      conjqm[1] *= 2.0;

      conjqm[2] *= 2.0;

      conjqm[3] *= 2.0;

      vec3_to_vec4(quat,torque[i],fquat);

      // update quaternion 4-momentum by 1/2 step

      conjqm[0] += dt * fquat[0];

      conjqm[1] += dt * fquat[1];

      conjqm[2] += dt * fquat[2];

      conjqm[3] += dt * fquat[3];

      // principal moments of inertia

      inertia[0] = INERTIA*rmass[i] * (shape[1]*shape[1]+shape[2]*shape[2]);

      inertia[1] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[2]*shape[2]);

      inertia[2] = INERTIA*rmass[i] * (shape[0]*shape[0]+shape[1]*shape[1]);

      // rotate quaternion and quaternion 4-momentum by full step

      no_squish_rotate(3,conjqm,quat,inertia,dthlf);

      no_squish_rotate(2,conjqm,quat,inertia,dthlf);

      no_squish_rotate(1,conjqm,quat,inertia,dt);

      no_squish_rotate(2,conjqm,quat,inertia,dthlf);

      no_squish_rotate(3,conjqm,quat,inertia,dthlf);

      qnormalize(quat);

      // convert quaternion 4-momentum in body frame back to angular momentum in space frame

      vec4_to_vec3(quat,conjqm,angmom[i]);

      angmom[i][0] *= 0.5;

      angmom[i][1] *= 0.5;

      angmom[i][2] *= 0.5;

    }

}

/* ---------------------------------------------------------------------- */

void FixNVEDot::final_integrate()

{

  double *quat;

  double fquat[4],conjqm[4];

  double conjqm_dot_quat;

  AtomVecEllipsoid::Bonus *bonus = avec->bonus;

  int *ellipsoid = atom->ellipsoid;

  double **v = atom->v;

  double **f = atom->f;

  double **angmom = atom->angmom;

  double **torque = atom->torque;

  double *rmass = atom->rmass;