ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3104 from ohenrich/cg-dna

Browse Source 
CG-DNA unit test and performance enhancement
This commit is contained in:
Axel Kohlmeyer
2022-02-01 18:30:13 -05:00
committed by GitHub
parent 76c57d54c1 e8ce01079d
commit 172f88a0d5
22 changed files with 975 additions and 186 deletions
Download Patch File Download Diff File Expand all files Collapse all files

47
src/CG-DNA/bond_oxdna_fene.cpp
View File

@ -26,6 +26,7 @@
#include "atom_vec_ellipsoid.h"
#include "math_extra.h"
#include "pair.h"
#include <cmath>
@ -145,16 +146,16 @@ void BondOxdnaFene::ev_tally_xyz(int i, int j, int nlocal, int newton_bond, doub
------------------------------------------------------------------------- */
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;
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;
// 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 ra_cs[3],rb_cs[3];
// Cartesian unit vectors in lab frame
double ax[3],ay[3],az[3];
double bx[3],by[3],bz[3];
double **x = atom->x;
double **f = atom->f;
@ -172,6 +173,12 @@ void BondOxdnaFene::compute(int eflag, int vflag)
ebond = 0.0;
ev_init(eflag, vflag);
// n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv
int dim;
nx_xtrct = (double **) force->pair->extract("nx",dim);
ny_xtrct = (double **) force->pair->extract("ny",dim);
nz_xtrct = (double **) force->pair->extract("nz",dim);
// loop over FENE bonds
for (in = 0; in < nbondlist; in++) {
@ -180,10 +187,24 @@ void BondOxdnaFene::compute(int eflag, int vflag)
b = bondlist[in][0];
type = bondlist[in][2];
qa = bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa, ax, ay, az);
qb = bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb, bx, by, bz);
ax[0] = nx_xtrct[a][0];
ax[1] = nx_xtrct[a][1];
ax[2] = nx_xtrct[a][2];
ay[0] = ny_xtrct[a][0];
ay[1] = ny_xtrct[a][1];
ay[2] = ny_xtrct[a][2];
az[0] = nz_xtrct[a][0];
az[1] = nz_xtrct[a][1];
az[2] = nz_xtrct[a][2];
bx[0] = nx_xtrct[b][0];
bx[1] = nx_xtrct[b][1];
bx[2] = nx_xtrct[b][2];
by[0] = ny_xtrct[b][0];
by[1] = ny_xtrct[b][1];
by[2] = ny_xtrct[b][2];
bz[0] = nz_xtrct[b][0];
bz[1] = nz_xtrct[b][1];
bz[2] = nz_xtrct[b][2];
// vector COM-backbone site a and b
compute_interaction_sites(ax, ay, az, ra_cs);

1
src/CG-DNA/bond_oxdna_fene.h
View File

@ -40,6 +40,7 @@ class BondOxdnaFene : public Bond {
protected:
double *k, *Delta, *r0; // FENE
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
void allocate();
void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double);

30
src/CG-DNA/pair_oxdna2_coaxstk.cpp
View File

@ -118,9 +118,9 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
double ra_cs[3],ra_cst[3];
double rb_cs[3],rb_cst[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];
// Cartesian unit vectors in lab frame
double ax[3],ay[3],az[3];
double bx[3],by[3],bz[3];
double **x = atom->x;
double **f = atom->f;
@ -149,6 +149,11 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// n(x/z)_xtrct = extracted local unit vectors from oxdna_excv
int dim;
nx_xtrct = (double **) force->pair->extract("nx",dim);
nz_xtrct = (double **) force->pair->extract("nz",dim);
// loop over pair interaction neighbors of my atoms
for (ia = 0; ia < anum; ia++) {
@ -156,8 +161,9 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
a = alist[ia];
atype = type[a];
qa=bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
ax[0] = nx_xtrct[a][0];
ax[1] = nx_xtrct[a][1];
ax[2] = nx_xtrct[a][2];
// vector COM a - stacking site a
ra_cst[0] = d_cst*ax[0];
@ -180,8 +186,9 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
btype = type[b];
qb=bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
bx[0] = nx_xtrct[b][0];
bx[1] = nx_xtrct[b][1];
bx[2] = nx_xtrct[b][2];
// vector COM b - stacking site b
rb_cst[0] = d_cst*bx[0];
@ -232,6 +239,13 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
// early rejection criterium
if (f4f6t1) {
az[0] = nz_xtrct[a][0];
az[1] = nz_xtrct[a][1];
az[2] = nz_xtrct[a][2];
bz[0] = nz_xtrct[b][0];
bz[1] = nz_xtrct[b][1];
bz[2] = nz_xtrct[b][2];
cost4 = MathExtra::dot3(az,bz);
if (cost4 > 1.0) cost4 = 1.0;
if (cost4 < -1.0) cost4 = -1.0;
@ -306,7 +320,7 @@ void PairOxdna2Coaxstk::compute(int eflag, int vflag)
DF4(theta6p, a_cxst6[atype][btype], theta_cxst6_0[atype][btype], dtheta_cxst6_ast[atype][btype],
b_cxst6[atype][btype], dtheta_cxst6_c[atype][btype])*rsint;
// force, torque and virial contribution for forces between stacking sites
// force, torque and virial contribution for forces between stacking sites
fpair = 0.0;

6
src/CG-DNA/pair_oxdna2_coaxstk.h
View File

@ -46,20 +46,16 @@ class PairOxdna2Coaxstk : public Pair {
double **k_cxst, **cut_cxst_0, **cut_cxst_c, **cut_cxst_lo, **cut_cxst_hi;
double **cut_cxst_lc, **cut_cxst_hc, **b_cxst_lo, **b_cxst_hi;
double **cutsq_cxst_hc;
double **a_cxst1, **theta_cxst1_0, **dtheta_cxst1_ast;
double **b_cxst1, **dtheta_cxst1_c;
double **a_cxst4, **theta_cxst4_0, **dtheta_cxst4_ast;
double **b_cxst4, **dtheta_cxst4_c;
double **a_cxst5, **theta_cxst5_0, **dtheta_cxst5_ast;
double **b_cxst5, **dtheta_cxst5_c;
double **a_cxst6, **theta_cxst6_0, **dtheta_cxst6_ast;
double **b_cxst6, **dtheta_cxst6_c;
double **AA_cxst1, **BB_cxst1;
double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
virtual void allocate();
};

36
src/CG-DNA/pair_oxdna2_dh.cpp
View File

@ -87,10 +87,9 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
// 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;
// Cartesian unit vectors in lab frame
double ax[3],ay[3],az[3];
double bx[3],by[3],bz[3];
double **x = atom->x;
double **f = atom->f;
@ -100,6 +99,7 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
int *alist,*blist,*numneigh,**firstneigh;
double *special_lj = force->special_lj;
AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
@ -115,6 +115,12 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv
int dim;
nx_xtrct = (double **) force->pair->extract("nx",dim);
ny_xtrct = (double **) force->pair->extract("ny",dim);
nz_xtrct = (double **) force->pair->extract("nz",dim);
// loop over pair interaction neighbors of my atoms
for (ia = 0; ia < anum; ia++) {
@ -122,8 +128,15 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
a = alist[ia];
atype = type[a];
qa=bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
ax[0] = nx_xtrct[a][0];
ax[1] = nx_xtrct[a][1];
ax[2] = nx_xtrct[a][2];
ay[0] = ny_xtrct[a][0];
ay[1] = ny_xtrct[a][1];
ay[2] = ny_xtrct[a][2];
az[0] = nz_xtrct[a][0];
az[1] = nz_xtrct[a][1];
az[2] = nz_xtrct[a][2];
// vector COM-backbone site a
compute_interaction_sites(ax,ay,az,ra_cs);
@ -142,8 +155,15 @@ void PairOxdna2Dh::compute(int eflag, int vflag)
b &= NEIGHMASK;
btype = type[b];
qb=bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
bx[0] = nx_xtrct[b][0];
bx[1] = nx_xtrct[b][1];
bx[2] = nx_xtrct[b][2];
by[0] = ny_xtrct[b][0];
by[1] = ny_xtrct[b][1];
by[2] = ny_xtrct[b][2];
bz[0] = nz_xtrct[b][0];
bz[1] = nz_xtrct[b][1];
bz[2] = nz_xtrct[b][2];
// vector COM-backbone site b
compute_interaction_sites(bx,by,bz,rb_cs);

1
src/CG-DNA/pair_oxdna2_dh.h
View File

@ -45,6 +45,7 @@ class PairOxdna2Dh : public Pair {
protected:
double **qeff_dh_pf, **kappa_dh;
double **b_dh, **cut_dh_ast, **cutsq_dh_ast, **cut_dh_c, **cutsq_dh_c;
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
virtual void allocate();
};

37
src/CG-DNA/pair_oxdna_coaxstk.cpp
View File

@ -106,7 +106,6 @@ PairOxdnaCoaxstk::~PairOxdnaCoaxstk()
void PairOxdnaCoaxstk::compute(int eflag, int vflag)
{
double delf[3],delt[3],delta[3],deltb[3]; // force, torque increment;
double evdwl,fpair,finc,tpair,factor_lj;
double v1tmp[3],v2tmp[3],v3tmp[3];
@ -129,10 +128,9 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
// vectors COM-backbone site, COM-stacking site in lab frame
double ra_cs[3],ra_cst[3];
double rb_cs[3],rb_cst[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];
// Cartesian unit vectors in lab frame
double ax[3],ay[3],az[3];
double bx[3],by[3],bz[3];
double **x = atom->x;
double **f = atom->f;
@ -161,6 +159,12 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv
int dim;
nx_xtrct = (double **) force->pair->extract("nx",dim);
ny_xtrct = (double **) force->pair->extract("ny",dim);
nz_xtrct = (double **) force->pair->extract("nz",dim);
// loop over pair interaction neighbors of my atoms
for (ia = 0; ia < anum; ia++) {
@ -168,8 +172,10 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
a = alist[ia];
atype = type[a];
qa=bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
ax[0] = nx_xtrct[a][0];
ax[1] = nx_xtrct[a][1];
ax[2] = nx_xtrct[a][2];
// a(y/z) not needed here as oxDNA(1) co-linear
// vector COM a - stacking site a
ra_cst[0] = d_cst*ax[0];
@ -192,8 +198,10 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
btype = type[b];
qb=bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
bx[0] = nx_xtrct[b][0];
bx[1] = nx_xtrct[b][1];
bx[2] = nx_xtrct[b][2];
// b(y/z) not needed here as oxDNA(1) co-linear
// vector COM b - stacking site b
rb_cst[0] = d_cst*bx[0];
@ -245,6 +253,13 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
// early rejection criterium
if (f4t1) {
az[0] = nz_xtrct[a][0];
az[1] = nz_xtrct[a][1];
az[2] = nz_xtrct[a][2];
bz[0] = nz_xtrct[b][0];
bz[1] = nz_xtrct[b][1];
bz[2] = nz_xtrct[b][2];
cost4 = MathExtra::dot3(az,bz);
if (cost4 > 1.0) cost4 = 1.0;
if (cost4 < -1.0) cost4 = -1.0;
@ -380,6 +395,10 @@ void PairOxdnaCoaxstk::compute(int eflag, int vflag)
// cosphi3 and cosphi4 (=cosphi3) force and virial
if (cosphi3) {
ay[0] = ny_xtrct[a][0];
ay[1] = ny_xtrct[a][1];
ay[2] = ny_xtrct[a][2];
finc = -f2 * f4t1* f4t4 * f4t5 * f4t6 * 2.0 * f5c3 * df5c3 * factor_lj;
fpair += finc;

6
src/CG-DNA/pair_oxdna_coaxstk.h
View File

@ -47,21 +47,17 @@ class PairOxdnaCoaxstk : public Pair {
double **k_cxst, **cut_cxst_0, **cut_cxst_c, **cut_cxst_lo, **cut_cxst_hi;
double **cut_cxst_lc, **cut_cxst_hc, **b_cxst_lo, **b_cxst_hi;
double **cutsq_cxst_hc;
double **a_cxst1, **theta_cxst1_0, **dtheta_cxst1_ast;
double **b_cxst1, **dtheta_cxst1_c;
double **a_cxst4, **theta_cxst4_0, **dtheta_cxst4_ast;
double **b_cxst4, **dtheta_cxst4_c;
double **a_cxst5, **theta_cxst5_0, **dtheta_cxst5_ast;
double **b_cxst5, **dtheta_cxst5_c;
double **a_cxst6, **theta_cxst6_0, **dtheta_cxst6_ast;
double **b_cxst6, **dtheta_cxst6_c;
double **a_cxst3p, **cosphi_cxst3p_ast, **b_cxst3p, **cosphi_cxst3p_c;
double **a_cxst4p, **cosphi_cxst4p_ast, **b_cxst4p, **cosphi_cxst4p_c;
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
virtual void allocate();
};

113
src/CG-DNA/pair_oxdna_excv.cpp
View File

@ -39,6 +39,9 @@ PairOxdnaExcv::PairOxdnaExcv(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
writedata = 1;
// set comm size needed by this Pair
comm_forward = 9;
}
/* ---------------------------------------------------------------------- */
@ -47,6 +50,10 @@ PairOxdnaExcv::~PairOxdnaExcv()
{
if (allocated) {
memory->destroy(nx);
memory->destroy(ny);
memory->destroy(nz);
memory->destroy(setflag);
memory->destroy(cutsq);
@ -108,7 +115,6 @@ void PairOxdnaExcv::compute_interaction_sites(double e1[3], double /*e2*/[3],
void PairOxdnaExcv::compute(int eflag, int vflag)
{
double delf[3],delta[3],deltb[3]; // force, torque increment;
double evdwl,fpair,factor_lj;
double rtmp_s[3],rtmp_b[3];
@ -118,10 +124,10 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
// vectors COM-backbone site, COM-base site in lab frame
double ra_cs[3],ra_cb[3];
double rb_cs[3],rb_cb[3];
// Cartesian unit vectors in lab frame
double ax[3],ay[3],az[3];
double bx[3],by[3],bz[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;
@ -137,7 +143,7 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
int *ellipsoid = atom->ellipsoid;
int a,b,ia,ib,anum,bnum,atype,btype;
int n,a,b,in,ia,ib,anum,bnum,atype,btype;
evdwl = 0.0;
ev_init(eflag,vflag);
@ -147,6 +153,29 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over all local atoms, calculation of local reference frame
for (in = 0; in < atom->nlocal; in++) {
int n = alist[in];
double *qn,nx_temp[3],ny_temp[3],nz_temp[3]; // quaternion and Cartesian unit vectors in lab frame
qn=bonus[ellipsoid[n]].quat;
MathExtra::q_to_exyz(qn,nx_temp,ny_temp,nz_temp);
nx[n][0] = nx_temp[0];
nx[n][1] = nx_temp[1];
nx[n][2] = nx_temp[2];
ny[n][0] = ny_temp[0];
ny[n][1] = ny_temp[1];
ny[n][2] = ny_temp[2];
nz[n][0] = nz_temp[0];
nz[n][1] = nz_temp[1];
nz[n][2] = nz_temp[2];
}
comm->forward_comm_pair(this);
// loop over pair interaction neighbors of my atoms
for (ia = 0; ia < anum; ia++) {
@ -154,8 +183,15 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
a = alist[ia];
atype = type[a];
qa=bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
ax[0] = nx[a][0];
ax[1] = nx[a][1];
ax[2] = nx[a][2];
ay[0] = ny[a][0];
ay[1] = ny[a][1];
ay[2] = ny[a][2];
az[0] = nz[a][0];
az[1] = nz[a][1];
az[2] = nz[a][2];
// vector COM - backbone and base site a
compute_interaction_sites(ax,ay,az,ra_cs,ra_cb);
@ -179,8 +215,15 @@ void PairOxdnaExcv::compute(int eflag, int vflag)
btype = type[b];
qb=bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
bx[0] = nx[b][0];
bx[1] = nx[b][1];
bx[2] = nx[b][2];
by[0] = ny[b][0];
by[1] = ny[b][1];
by[2] = ny[b][2];
bz[0] = nz[b][0];
bz[1] = nz[b][1];
bz[2] = nz[b][2];
// vector COM - backbone and base site b
compute_interaction_sites(bx,by,bz,rb_cs,rb_cb);
@ -400,6 +443,10 @@ void PairOxdnaExcv::allocate()
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(nx,atom->nmax,3,"pair:nx");
memory->create(ny,atom->nmax,3,"pair:ny");
memory->create(nz,atom->nmax,3,"pair:nz");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(epsilon_ss,n+1,n+1,"pair:epsilon_ss");
@ -806,10 +853,58 @@ void PairOxdnaExcv::write_data_all(FILE *fp)
/* ---------------------------------------------------------------------- */
int PairOxdnaExcv::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
{
int i,j,m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = nx[j][0];
buf[m++] = nx[j][1];
buf[m++] = nx[j][2];
buf[m++] = ny[j][0];
buf[m++] = ny[j][1];
buf[m++] = ny[j][2];
buf[m++] = nz[j][0];
buf[m++] = nz[j][1];
buf[m++] = nz[j][2];
}
return m;
}
/* ---------------------------------------------------------------------- */
void PairOxdnaExcv::unpack_forward_comm(int n, int first, double *buf)
{
int i,m,last;
m = 0;
last = first + n;
for (i = first; i < last; i++) {
nx[i][0] = buf[m++];
nx[i][1] = buf[m++];
nx[i][2] = buf[m++];
ny[i][0] = buf[m++];
ny[i][1] = buf[m++];
ny[i][2] = buf[m++];
nz[i][0] = buf[m++];
nz[i][1] = buf[m++];
nz[i][2] = buf[m++];
}
}
/* ---------------------------------------------------------------------- */
void *PairOxdnaExcv::extract(const char *str, int &dim)
{
dim = 2;
if (strcmp(str,"nx") == 0) return (void *) nx;
if (strcmp(str,"ny") == 0) return (void *) ny;
if (strcmp(str,"nz") == 0) return (void *) nz;
if (strcmp(str,"epsilon_ss") == 0) return (void *) epsilon_ss;
if (strcmp(str,"sigma_ss") == 0) return (void *) sigma_ss;
if (strcmp(str,"cut_ss_ast") == 0) return (void *) cut_ss_ast;

3
src/CG-DNA/pair_oxdna_excv.h
View File

@ -41,6 +41,8 @@ class PairOxdnaExcv : public Pair {
void write_data(FILE *) override;
void write_data_all(FILE *) override;
void *extract(const char *, int &) override;
virtual int pack_forward_comm(int, int *, double *, int, int *);
virtual void unpack_forward_comm(int, int, double *);
protected:
// s=sugar-phosphate backbone site, b=base site, st=stacking site
@ -52,6 +54,7 @@ class PairOxdnaExcv : public Pair {
double **lj1_sb, **lj2_sb, **b_sb, **cut_sb_c, **cutsq_sb_c;
double **epsilon_bb, **sigma_bb, **cut_bb_ast, **cutsq_bb_ast;
double **lj1_bb, **lj2_bb, **b_bb, **cut_bb_c, **cutsq_bb_c;
double **nx, **ny, **nz; // per-atom arrays for local unit vectors
virtual void allocate();
};

30
src/CG-DNA/pair_oxdna_hbond.cpp
View File

@ -148,10 +148,9 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
double d_chb=+0.4;
// vectors COM-h-bonding site in lab frame
double ra_chb[3],rb_chb[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];
// Cartesian unit vectors in lab frame
double ax[3],ay[3],az[3];
double bx[3],by[3],bz[3];
double **x = atom->x;
double **f = atom->f;
@ -180,6 +179,12 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv
int dim;
nx_xtrct = (double **) force->pair->extract("nx",dim);
ny_xtrct = (double **) force->pair->extract("ny",dim);
nz_xtrct = (double **) force->pair->extract("nz",dim);
// loop over pair interaction neighbors of my atoms
for (ia = 0; ia < anum; ia++) {
@ -187,8 +192,9 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
a = alist[ia];
atype = type[a];
qa=bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
ax[0] = nx_xtrct[a][0];
ax[1] = nx_xtrct[a][1];
ax[2] = nx_xtrct[a][2];
ra_chb[0] = d_chb*ax[0];
ra_chb[1] = d_chb*ax[1];
@ -205,8 +211,9 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
btype = type[b];
qb=bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
bx[0] = nx_xtrct[b][0];
bx[1] = nx_xtrct[b][1];
bx[2] = nx_xtrct[b][2];
rb_chb[0] = d_chb*bx[0];
rb_chb[1] = d_chb*bx[1];
@ -265,6 +272,13 @@ void PairOxdnaHbond::compute(int eflag, int vflag)
// early rejection criterium
if (f4t3) {
az[0] = nz_xtrct[a][0];
az[1] = nz_xtrct[a][1];
az[2] = nz_xtrct[a][2];
bz[0] = nz_xtrct[b][0];
bz[1] = nz_xtrct[b][1];
bz[2] = nz_xtrct[b][2];
cost4 = MathExtra::dot3(az,bz);
if (cost4 > 1.0) cost4 = 1.0;
if (cost4 < -1.0) cost4 = -1.0;

8
src/CG-DNA/pair_oxdna_hbond.h
View File

@ -48,25 +48,19 @@ class PairOxdnaHbond : public Pair {
double **epsilon_hb, **a_hb, **cut_hb_0, **cut_hb_c, **cut_hb_lo, **cut_hb_hi;
double **cut_hb_lc, **cut_hb_hc, **b_hb_lo, **b_hb_hi, **shift_hb;
double **cutsq_hb_hc;
double **a_hb1, **theta_hb1_0, **dtheta_hb1_ast;
double **b_hb1, **dtheta_hb1_c;
double **a_hb2, **theta_hb2_0, **dtheta_hb2_ast;
double **b_hb2, **dtheta_hb2_c;
double **a_hb3, **theta_hb3_0, **dtheta_hb3_ast;
double **b_hb3, **dtheta_hb3_c;
double **a_hb4, **theta_hb4_0, **dtheta_hb4_ast;
double **b_hb4, **dtheta_hb4_c;
double **a_hb7, **theta_hb7_0, **dtheta_hb7_ast;
double **b_hb7, **dtheta_hb7_c;
double **a_hb8, **theta_hb8_0, **dtheta_hb8_ast;
double **b_hb8, **dtheta_hb8_c;
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
int seqdepflag;
virtual void allocate();

39
src/CG-DNA/pair_oxdna_stk.cpp
View File

@ -212,7 +212,6 @@ void PairOxdnaStk::ev_tally_xyz(int i, int j, int nlocal, int newton_bond,
void PairOxdnaStk::compute(int eflag, int vflag)
{
double delf[3],delta[3],deltb[3]; // force, torque increment;
double evdwl,fpair,finc,tpair;
double delr_ss[3],delr_ss_norm[3],rsq_ss,r_ss,rinv_ss;
@ -227,10 +226,9 @@ void PairOxdnaStk::compute(int eflag, int vflag)
// vectors COM-backbone site, COM-stacking site in lab frame
double ra_cs[3],ra_cst[3];
double rb_cs[3],rb_cst[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];
// Cartesian unit vectors in lab frame
double ax[3],ay[3],az[3];
double bx[3],by[3],bz[3];
double **x = atom->x;
double **f = atom->f;
@ -257,6 +255,12 @@ void PairOxdnaStk::compute(int eflag, int vflag)
evdwl = 0.0;
ev_init(eflag,vflag);
// n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv
int dim;
nx_xtrct = (double **) force->pair->extract("nx",dim);
ny_xtrct = (double **) force->pair->extract("ny",dim);
nz_xtrct = (double **) force->pair->extract("nz",dim);
// loop over stacking interaction neighbors using bond topology
for (in = 0; in < nbondlist; in++) {
@ -275,10 +279,13 @@ void PairOxdnaStk::compute(int eflag, int vflag)
// a now in 3' direction, b in 5' direction
qa=bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
qb=bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
ax[0] = nx_xtrct[a][0];
ax[1] = nx_xtrct[a][1];
ax[2] = nx_xtrct[a][2];
bx[0] = nx_xtrct[b][0];
bx[1] = nx_xtrct[b][1];
bx[2] = nx_xtrct[b][2];
// (a/b)y/z not needed here as oxDNA(1) co-linear
// vector COM a - stacking site a
ra_cst[0] = d_cst*ax[0];
@ -336,6 +343,13 @@ void PairOxdnaStk::compute(int eflag, int vflag)
// early rejection criterium
if (f1) {
az[0] = nz_xtrct[a][0];
az[1] = nz_xtrct[a][1];
az[2] = nz_xtrct[a][2];
bz[0] = nz_xtrct[b][0];
bz[1] = nz_xtrct[b][1];
bz[2] = nz_xtrct[b][2];
// theta4 angle and correction
cost4 = MathExtra::dot3(bz,az);
if (cost4 > 1.0) cost4 = 1.0;
@ -360,6 +374,13 @@ void PairOxdnaStk::compute(int eflag, int vflag)
// early rejection criterium
if (f4t5) {
ay[0] = ny_xtrct[a][0];
ay[1] = ny_xtrct[a][1];
ay[2] = ny_xtrct[a][2];
by[0] = ny_xtrct[b][0];
by[1] = ny_xtrct[b][1];
by[2] = ny_xtrct[b][2];
cost6p = MathExtra::dot3(delr_st_norm,az);
if (cost6p > 1.0) cost6p = 1.0;
if (cost6p < -1.0) cost6p = -1.0;

2
src/CG-DNA/pair_oxdna_stk.h
View File

@ -59,7 +59,7 @@ class PairOxdnaStk : public Pair {
double **b_st6, **dtheta_st6_c;
double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c;
double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c;
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
int seqdepflag;
virtual void allocate();

33
src/CG-DNA/pair_oxdna_xstk.cpp
View File

@ -111,7 +111,6 @@ PairOxdnaXstk::~PairOxdnaXstk()
void PairOxdnaXstk::compute(int eflag, int vflag)
{
double delf[3],delta[3],deltb[3]; // force, torque increment;
double evdwl,fpair,finc,tpair,factor_lj;
double delr_hb[3],delr_hb_norm[3],rsq_hb,r_hb,rinv_hb;
@ -126,10 +125,9 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
double d_chb=+0.4;
// vectors COM-h-bonding site in lab frame
double ra_chb[3],rb_chb[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];
// Cartesian unit vectors in lab frame
double ax[3],ay[3],az[3];
double bx[3],by[3],bz[3];
double **x = atom->x;
double **f = atom->f;
@ -158,6 +156,12 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv
int dim;
nx_xtrct = (double **) force->pair->extract("nx",dim);
ny_xtrct = (double **) force->pair->extract("ny",dim);
nz_xtrct = (double **) force->pair->extract("nz",dim);
// loop over pair interaction neighbors of my atoms
for (ia = 0; ia < anum; ia++) {
@ -165,8 +169,10 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
a = alist[ia];
atype = type[a];
qa=bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
ax[0] = nx_xtrct[a][0];
ax[1] = nx_xtrct[a][1];
ax[2] = nx_xtrct[a][2];
// a(y/z) not needed here as oxDNA(1) co-linear
ra_chb[0] = d_chb*ax[0];
ra_chb[1] = d_chb*ax[1];
@ -183,8 +189,10 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
btype = type[b];
qb=bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
bx[0] = nx_xtrct[b][0];
bx[1] = nx_xtrct[b][1];
bx[2] = nx_xtrct[b][2];
// b(y/z) not needed here as oxDNA(1) co-linear
rb_chb[0] = d_chb*bx[0];
rb_chb[1] = d_chb*bx[1];
@ -243,6 +251,13 @@ void PairOxdnaXstk::compute(int eflag, int vflag)
// early rejection criterium
if (f4t3) {
az[0] = nz_xtrct[a][0];
az[1] = nz_xtrct[a][1];
az[2] = nz_xtrct[a][2];
bz[0] = nz_xtrct[b][0];
bz[1] = nz_xtrct[b][1];
bz[2] = nz_xtrct[b][2];
cost4 = MathExtra::dot3(az,bz);
if (cost4 > 1.0) cost4 = 1.0;
if (cost4 < -1.0) cost4 = -1.0;

7
src/CG-DNA/pair_oxdna_xstk.h
View File

@ -47,24 +47,19 @@ class PairOxdnaXstk : public Pair {
double **k_xst, **cut_xst_0, **cut_xst_c, **cut_xst_lo, **cut_xst_hi;
double **cut_xst_lc, **cut_xst_hc, **b_xst_lo, **b_xst_hi;
double **cutsq_xst_hc;
double **a_xst1, **theta_xst1_0, **dtheta_xst1_ast;
double **b_xst1, **dtheta_xst1_c;
double **a_xst2, **theta_xst2_0, **dtheta_xst2_ast;
double **b_xst2, **dtheta_xst2_c;
double **a_xst3, **theta_xst3_0, **dtheta_xst3_ast;
double **b_xst3, **dtheta_xst3_c;
double **a_xst4, **theta_xst4_0, **dtheta_xst4_ast;
double **b_xst4, **dtheta_xst4_c;
double **a_xst7, **theta_xst7_0, **dtheta_xst7_ast;
double **b_xst7, **dtheta_xst7_c;
double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast;
double **b_xst8, **dtheta_xst8_c;
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
virtual void allocate();
};

36
src/CG-DNA/pair_oxrna2_stk.cpp
View File

@ -245,9 +245,9 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
double ra_cs[3],ra_cst[3];
double rb_cs[3],rb_cst[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];
// Cartesian unit vectors in lab frame
double ax[3],ay[3],az[3];
double bx[3],by[3],bz[3];
double **x = atom->x;
double **f = atom->f;
@ -274,6 +274,12 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
evdwl = 0.0;
ev_init(eflag,vflag);
// n(x/y/z)_xtrct = extracted local unit vectors from oxdna_excv
int dim;
nx_xtrct = (double **) force->pair->extract("nx",dim);
ny_xtrct = (double **) force->pair->extract("ny",dim);
nz_xtrct = (double **) force->pair->extract("nz",dim);
// loop over stacking interaction neighbors using bond topology
for (in = 0; in < nbondlist; in++) {
@ -290,12 +296,26 @@ void PairOxrna2Stk::compute(int eflag, int vflag)
}
// a now in 3' direction, b in 5' direction
// a now in 3' direction, b in 5' direction
qa=bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
qb=bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
ax[0] = nx_xtrct[a][0];
ax[1] = nx_xtrct[a][1];
ax[2] = nx_xtrct[a][2];
ay[0] = ny_xtrct[a][0];
ay[1] = ny_xtrct[a][1];
ay[2] = ny_xtrct[a][2];
az[0] = nz_xtrct[a][0];
az[1] = nz_xtrct[a][1];
az[2] = nz_xtrct[a][2];
bx[0] = nx_xtrct[b][0];
bx[1] = nx_xtrct[b][1];
bx[2] = nx_xtrct[b][2];
by[0] = ny_xtrct[b][0];
by[1] = ny_xtrct[b][1];
by[2] = ny_xtrct[b][2];
bz[0] = nz_xtrct[b][0];
bz[1] = nz_xtrct[b][1];
bz[2] = nz_xtrct[b][2];
// vector COM a - 5'-stacking site a
ra_cst[0] = d_cst_x_5p*ax[0] + d_cst_y_5p*ay[0];

2
src/CG-DNA/pair_oxrna2_stk.h
View File

@ -60,7 +60,7 @@ class PairOxrna2Stk : public Pair {
double **b_st10, **dtheta_st10_c;
double **a_st1, **cosphi_st1_ast, **b_st1, **cosphi_st1_c;
double **a_st2, **cosphi_st2_ast, **b_st2, **cosphi_st2_c;
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
int seqdepflag;
virtual void allocate();

29
src/CG-DNA/pair_oxrna2_xstk.cpp
View File

@ -120,9 +120,9 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
// vectors COM-h-bonding site in lab frame
double ra_chb[3],rb_chb[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];
// Cartesian unit vectors in lab frame
double ax[3],ay[3],az[3];
double bx[3],by[3],bz[3];
double **x = atom->x;
double **f = atom->f;
@ -151,6 +151,11 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// n(x/z)_xtrct = extracted local unit vectors from oxdna_excv
int dim;
nx_xtrct = (double **) force->pair->extract("nx",dim);
nz_xtrct = (double **) force->pair->extract("nz",dim);
// loop over pair interaction neighbors of my atoms
for (ia = 0; ia < anum; ia++) {
@ -158,8 +163,9 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
a = alist[ia];
atype = type[a];
qa=bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa,ax,ay,az);
ax[0] = nx_xtrct[a][0];
ax[1] = nx_xtrct[a][1];
ax[2] = nx_xtrct[a][2];
ra_chb[0] = d_chb*ax[0];
ra_chb[1] = d_chb*ax[1];
@ -176,8 +182,9 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
btype = type[b];
qb=bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb,bx,by,bz);
bx[0] = nx_xtrct[b][0];
bx[1] = nx_xtrct[b][1];
bx[2] = nx_xtrct[b][2];
rb_chb[0] = d_chb*bx[0];
rb_chb[1] = d_chb*bx[1];
@ -236,6 +243,10 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
// early rejection criterium
if (f4t3) {
az[0] = nz_xtrct[a][0];
az[1] = nz_xtrct[a][1];
az[2] = nz_xtrct[a][2];
cost7 = -1.0*MathExtra::dot3(az,delr_hb_norm);
if (cost7 > 1.0) cost7 = 1.0;
if (cost7 < -1.0) cost7 = -1.0;
@ -250,6 +261,10 @@ void PairOxrna2Xstk::compute(int eflag, int vflag)
// early rejection criterium
if (f4t7) {
bz[0] = nz_xtrct[b][0];
bz[1] = nz_xtrct[b][1];
bz[2] = nz_xtrct[b][2];
cost8 = MathExtra::dot3(bz,delr_hb_norm);
if (cost8 > 1.0) cost8 = 1.0;
if (cost8 < -1.0) cost8 = -1.0;

6
src/CG-DNA/pair_oxrna2_xstk.h
View File

@ -46,21 +46,17 @@ class PairOxrna2Xstk : public Pair {
double **k_xst, **cut_xst_0, **cut_xst_c, **cut_xst_lo, **cut_xst_hi;
double **cut_xst_lc, **cut_xst_hc, **b_xst_lo, **b_xst_hi;
double **cutsq_xst_hc;
double **a_xst1, **theta_xst1_0, **dtheta_xst1_ast;
double **b_xst1, **dtheta_xst1_c;
double **a_xst2, **theta_xst2_0, **dtheta_xst2_ast;
double **b_xst2, **dtheta_xst2_c;
double **a_xst3, **theta_xst3_0, **dtheta_xst3_ast;
double **b_xst3, **dtheta_xst3_c;
double **a_xst7, **theta_xst7_0, **dtheta_xst7_ast;
double **b_xst7, **dtheta_xst7_c;
double **a_xst8, **theta_xst8_0, **dtheta_xst8_ast;
double **b_xst8, **dtheta_xst8_c;
double **nx_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
virtual void allocate();
};