ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
5a1e020bf01c8a3454dd9159fe41de4ec4dbffa2
lammps/src/USER-CGDNA/pair_oxdna2_dh.cpp
Oliver Henrich 8afed61db1 Upgrade to oxDNA2
2017-03-28 18:16:36 +01:00

554 lines
16 KiB
C++
Raw Blame History

/* ----------------------------------------------------------------------
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 (University of Strathclyde, Glasgow)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "pair_oxdna2_dh.h"
#include "mf_oxdna.h"
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "integrate.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "atom_vec_ellipsoid.h"
#include "math_extra.h"
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace MFOxdna;
/* ---------------------------------------------------------------------- */
PairOxdna2Dh::PairOxdna2Dh(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
writedata = 1;
}
/* ---------------------------------------------------------------------- */
PairOxdna2Dh::~PairOxdna2Dh()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(qeff_dh_pf);
memory->destroy(kappa_dh);
memory->destroy(b_dh);
memory->destroy(cut_dh_ast);
memory->destroy(cutsq_dh_ast);
memory->destroy(cut_dh_c);
memory->destroy(cutsq_dh_c);
}
}
/* ----------------------------------------------------------------------
compute vector COM-sugar-phosphate backbone interaction site in oxDNA2
------------------------------------------------------------------------- */
void PairOxdna2Dh::compute_interaction_sites(double e1[3],
double e2[3], double r[3])
{
double d_cs_x=-0.34, d_cs_y=+0.3408;
r[0] = d_cs_x*e1[0] + d_cs_y*e2[0];
r[1] = d_cs_x*e1[1] + d_cs_y*e2[1];
r[2] = d_cs_x*e1[2] + d_cs_y*e2[2];
}
/* ----------------------------------------------------------------------
compute function for oxDNA pair interactions
s=sugar-phosphate backbone site, b=base site, st=stacking site
------------------------------------------------------------------------- */
void PairOxdna2Dh::compute(int eflag, int vflag)
{
double delf[3],delta[3],deltb[3]; // force, torque increment;;
double rtmp_s[3],delr[3];
double evdwl,fpair,factor_lj;
double r,rsq,rinv;
// vectors COM-backbone sites in lab frame
double ra_cs[3],rb_cs[3];
// quaternions and Cartesian unit vectors in lab frame
double *qa,ax[3],ay[3],az[3];
double *qb,bx[3],by[3],bz[3];
double *special_lj = force->special_lj;
double **x = atom->x;
double **f = atom->f;
double **torque = atom->torque;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
int *alist,*blist,*numneigh,**firstneigh;
AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
int a,b,ia,ib,anum,bnum,atype,btype;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
anum = list->inum;
alist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over pair interaction neighbours of my atoms
for (ia = 0; ia < anum; ia++) {
a = alist[ia];
atype = type[a];
qa=bonus[a].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
// vector COM-backbone site a
compute_interaction_sites(ax,ay,ra_cs);
rtmp_s[0] = x[a][0] + ra_cs[0];
rtmp_s[1] = x[a][1] + ra_cs[1];
rtmp_s[2] = x[a][2] + ra_cs[2];
blist = firstneigh[a];
bnum = numneigh[a];
for (ib = 0; ib < bnum; ib++) {
b = blist[ib];
factor_lj = special_lj[sbmask(b)]; // = 0 for nearest neighbours
b &= NEIGHMASK;
btype = type[b];
qb=bonus[b].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
// vector COM-backbone site b
compute_interaction_sites(bx,by,rb_cs);
// vector backbone site b to a
delr[0] = rtmp_s[0] - x[b][0] - rb_cs[0];
delr[1] = rtmp_s[1] - x[b][1] - rb_cs[1];
delr[2] = rtmp_s[2] - x[b][2] - rb_cs[2];
rsq = delr[0]*delr[0] + delr[1]*delr[1] + delr[2]*delr[2];
if (rsq <= cutsq_dh_c[atype][btype]) {
r = sqrt(rsq);
rinv = 1.0/r;
if (r <= cut_dh_ast[atype][btype]) {
fpair = qeff_dh_pf[atype][btype] * exp(-kappa_dh[atype][btype] * r) *
(kappa_dh[atype][btype] + rinv) * rinv * rinv;
if (eflag) {
evdwl = qeff_dh_pf[atype][btype] * exp(-kappa_dh[atype][btype]*r) * rinv;
}
}
else {
fpair = 2.0 * b_dh[atype][btype] * (cut_dh_c[atype][btype] - r) * rinv;
if (eflag) {
evdwl = b_dh[atype][btype] * (r - cut_dh_c[atype][btype]) * (r - cut_dh_c[atype][btype]);
}
}
// knock out nearest-neighbour interaction between adjacent backbone sites
fpair *= factor_lj;
evdwl *= factor_lj;
delf[0] = delr[0] * fpair;
delf[1] = delr[1] * fpair;
delf[2] = delr[2] * fpair;
// apply force and torque to each of 2 atoms
if (newton_pair || 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_pair || 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_pair,
evdwl,0.0,fpair,delr[0],delr[1],delr[2]);
}
}
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairOxdna2Dh::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(kappa_dh,n+1,n+1,"pair:kappa_dh");
memory->create(qeff_dh_pf,n+1,n+1,"pair:qeff_dh_pf");
memory->create(b_dh,n+1,n+1,"pair:b_dh");
memory->create(cut_dh_ast,n+1,n+1,"pair:cut_dh_ast");
memory->create(cutsq_dh_ast,n+1,n+1,"pair:cutsq_dh_ast");
memory->create(cut_dh_c,n+1,n+1,"pair:cut_dh_c");
memory->create(cutsq_dh_c,n+1,n+1,"pair:cutsq_dh_c");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairOxdna2Dh::settings(int narg, char **arg)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairOxdna2Dh::coeff(int narg, char **arg)
{
int count;
if (narg != 5) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/dh");
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi);
force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi);
count = 0;
double T, rhos_dh_one, qeff_dh_one;
T = force->numeric(FLERR,arg[2]);
rhos_dh_one = force->numeric(FLERR,arg[3]);
qeff_dh_one = force->numeric(FLERR,arg[4]);
double lambda_dh_one, kappa_dh_one, qeff_dh_pf_one;
double b_dh_one, cut_dh_ast_one, cut_dh_c_one;
// Debye length and inverse Debye length
/*
NOTE:
The numerical factor is the Debye length in s.u.
lambda(T = 300 K = 0.1) =
sqrt(eps_0 * eps_r * k_B * T/(2 * N_A * e^2 * 1000 mol/m^3))
* 1/oxDNA_energy_unit
(see B. Snodin et al., J. Chem. Phys. 142, 234901 (2015).)
We use
eps_0 = vacuum permittivity = 8.854187817e-12 F/m
eps_r = relative permittivity of water = 80
k_B = Boltzmann constant = 1.3806485279e-23 J/K
T = absolute temperature = 300 K
N_A = Avogadro constant = 6.02214085774e23 / mol
e = elementary charge = 1.6021766208e-19 C
oxDNA_length_unit = 8.518e-10 m
*/
lambda_dh_one = 0.3616455075438555*sqrt(T/0.1/rhos_dh_one);
kappa_dh_one = 1.0/lambda_dh_one;
// prefactor in DH interaction containing qeff^2
/*
NOTE:
The numerical factor is
qeff_dh_pf = e^2/(4 * pi * eps_0 * eps_r)
* 1/(oxDNA_energy_unit * oxDNA_length_unit)
(see B. Snodin et al., J. Chem. Phys. 142, 234901 (2015).)
In addition to the above units we use
oxDNA_energy_unit = 4.142e-20 J
*/
qeff_dh_pf_one = 0.08173808693529228*qeff_dh_one*qeff_dh_one;
// smoothing parameters - determined through continuity and differentiability
cut_dh_ast_one = 3.0*lambda_dh_one;
b_dh_one = -(exp(-cut_dh_ast_one/lambda_dh_one) * qeff_dh_pf_one * qeff_dh_pf_one *
(cut_dh_ast_one + lambda_dh_one) * (cut_dh_ast_one + lambda_dh_one))/
(-4.0 * cut_dh_ast_one * cut_dh_ast_one * cut_dh_ast_one *
lambda_dh_one * lambda_dh_one * qeff_dh_pf_one);
cut_dh_c_one = cut_dh_ast_one * (qeff_dh_pf_one*cut_dh_ast_one +
3.0*qeff_dh_pf_one * lambda_dh_one)/
(qeff_dh_pf_one * (cut_dh_ast_one+lambda_dh_one));
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
kappa_dh[i][j] = kappa_dh_one;
qeff_dh_pf[i][j] = qeff_dh_pf_one;
b_dh[i][j] = b_dh_one;
cut_dh_ast[i][j] = cut_dh_ast_one;
cut_dh_c[i][j] = cut_dh_c_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients in oxdna/dh");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairOxdna2Dh::init_style()
{
int irequest;
// request regular neighbor lists
irequest = neighbor->request(this,instance_me);
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use regular
------------------------------------------------------------------------- */
void PairOxdna2Dh::init_list(int id, NeighList *ptr)
{
if (id == 0) list = ptr;
if (id > 0) error->all(FLERR,"Respa not supported");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairOxdna2Dh::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
error->all(FLERR,"Coefficient mixing not defined in oxDNA");
}
if (offset_flag) {
error->all(FLERR,"Offset not supported in oxDNA");
}
kappa_dh[j][i] = kappa_dh[i][j];
qeff_dh_pf[j][i] = qeff_dh_pf[i][j];
b_dh[j][i] = b_dh[i][j];
cut_dh_ast[j][i] = cut_dh_ast[i][j];
cut_dh_c[j][i] = cut_dh_c[i][j];
cutsq_dh_ast[i][j] = cut_dh_ast[i][j]*cut_dh_ast[i][j];
cutsq_dh_ast[j][i] = cutsq_dh_ast[i][j];
cutsq_dh_c[i][j] = cut_dh_c[i][j]*cut_dh_c[i][j];
cutsq_dh_c[j][i] = cutsq_dh_c[i][j];
// set the master list distance cutoff
return cut_dh_c[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairOxdna2Dh::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i,j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp);
if (setflag[i][j]) {
fwrite(&kappa_dh[i][j],sizeof(double),1,fp);
fwrite(&qeff_dh_pf[i][j],sizeof(double),1,fp);
fwrite(&b_dh[i][j],sizeof(double),1,fp);
fwrite(&cut_dh_ast[i][j],sizeof(double),1,fp);
fwrite(&cut_dh_c[i][j],sizeof(double),1,fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairOxdna2Dh::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i,j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
if (setflag[i][j]) {
if (me == 0) {
fread(&kappa_dh[i][j],sizeof(double),1,fp);
fread(&qeff_dh_pf[i][j],sizeof(double),1,fp);
fread(&b_dh[i][j],sizeof(double),1,fp);
fread(&cut_dh_ast[i][j],sizeof(double),1,fp);
fread(&cut_dh_c[i][j],sizeof(double),1,fp);
}
MPI_Bcast(&kappa_dh[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&qeff_dh_pf[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&b_dh[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_dh_ast[i][j],1,MPI_DOUBLE,0,world);
MPI_Bcast(&cut_dh_c[i][j],1,MPI_DOUBLE,0,world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairOxdna2Dh::write_restart_settings(FILE *fp)
{
fwrite(&offset_flag,sizeof(int),1,fp);
fwrite(&mix_flag,sizeof(int),1,fp);
fwrite(&tail_flag,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairOxdna2Dh::read_restart_settings(FILE *fp)
{
int me = comm->me;
if (me == 0) {
fread(&offset_flag,sizeof(int),1,fp);
fread(&mix_flag,sizeof(int),1,fp);
fread(&tail_flag,sizeof(int),1,fp);
}
MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void PairOxdna2Dh::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp,"%d\
%g %g\
%g %g %g\
\n",i,
kappa_dh[i][i],qeff_dh_pf[i][i],
b_dh[i][i],cut_dh_ast[i][i],cut_dh_c[i][i]);
}
/* ----------------------------------------------------------------------
proc 0 writes all pairs to data file
------------------------------------------------------------------------- */
void PairOxdna2Dh::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d\
%g %g\
%g %g %g\
\n",i,j,
kappa_dh[i][j],qeff_dh_pf[i][j],
b_dh[i][j],cut_dh_ast[i][j],cut_dh_c[i][j]);
}
/* ---------------------------------------------------------------------- */
void *PairOxdna2Dh::extract(const char *str, int &dim)
{
dim = 2;
if (strcmp(str,"kappa_dh") == 0) return (void *) kappa_dh;
if (strcmp(str,"qeff_dh_pf") == 0) return (void *) qeff_dh_pf;
if (strcmp(str,"b_dh") == 0) return (void *) b_dh;
if (strcmp(str,"cut_dh_ast") == 0) return (void *) cut_dh_ast;
if (strcmp(str,"cut_dh_c") == 0) return (void *) cut_dh_c;
return NULL;
}
Reference in New Issue View Git Blame Copy Permalink