From a110b1d47564963ffd68249817cec722142cae7c Mon Sep 17 00:00:00 2001 From: Germain Clavier Date: Sun, 6 Feb 2022 20:00:55 +0100 Subject: [PATCH] 1) Added born_matrix method to bond_harmonic, angle_cosine_squared and dihedral_nharmonic 2) Changed compute_born_matrix so it can takes arguments pair bond angle dihedral to output only selected elements. Also moved flags to creator method and added error/warning messages --- src/EXTRA-COMPUTE/compute_born_matrix.cpp | 164 ++++++++++++++-------- src/EXTRA-MOLECULE/dihedral_nharmonic.cpp | 61 ++++++++ src/EXTRA-MOLECULE/dihedral_nharmonic.h | 1 + src/MOLECULE/angle_cosine_squared.cpp | 30 ++++ src/MOLECULE/angle_cosine_squared.h | 1 + src/MOLECULE/bond_harmonic.cpp | 15 ++ src/MOLECULE/bond_harmonic.h | 1 + 7 files changed, 215 insertions(+), 58 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.cpp b/src/EXTRA-COMPUTE/compute_born_matrix.cpp index 682120626f..70b6535c7a 100644 --- a/src/EXTRA-COMPUTE/compute_born_matrix.cpp +++ b/src/EXTRA-COMPUTE/compute_born_matrix.cpp @@ -43,13 +43,26 @@ using namespace LAMMPS_NS; #define BIG 1000000000 -static int const albe[21][2] = { +// This table is used to pick the 3d rij vector indices used to +// compute the 6 indices long Voigt stress vector +static int const sigma_albe[6][2] = { + {0, 0}, // s11 + {1, 1}, // s22 + {2, 2}, // s33 + {1, 2}, // s44 + {0, 2}, // s55 + {0, 1}, // s66 +}; + +// This table is used to pick the correct indices from the Voigt +// stress vector to compute the Cij matrix (21 terms, see doc) contribution +static int const C_albe[21][2] = { {0, 0}, // C11 {1, 1}, // C22 {2, 2}, // C33 - {1, 2}, // C44 - {0, 2}, // C55 - {0, 1}, // C66 + {3, 3}, // C44 + {4, 4}, // C55 + {5, 5}, // C66 {0, 1}, // C12 {0, 2}, // C13 {0, 3}, // C14 @@ -67,6 +80,8 @@ static int const albe[21][2] = { {4, 5} // C56 }; +// This table is used to pick the 3d rij vector indices used to +// compute the 21 indices long Cij matrix static int const albemunu[21][4] = { {0, 0, 0, 0}, // C11 {1, 1, 1, 1}, // C22 @@ -103,6 +118,75 @@ ComputeBornMatrix::ComputeBornMatrix(LAMMPS *lmp, int narg, char **arg) : Comput // Error check + if (narg < 3) error->all(FLERR,"Illegal compute pressure command"); + + pairflag = 0; + bondflag = 0; + angleflag = 0; + dihedflag = 0; + impflag = 0; + kspaceflag = 0; + if (narg == 3) { + pairflag = 1; + bondflag = 1; + angleflag = 1; + dihedflag = 1; + impflag = 1; + kspaceflag = 1; + } else { + int iarg = 3; + while (iarg < narg) { + if (strcmp(arg[iarg],"pair") == 0) pairflag = 1; + else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1; + else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1; + else if (strcmp(arg[iarg],"dihedral") == 0) dihedflag = 1; + else if (strcmp(arg[iarg],"improper") == 0) impflag = 1; + else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1; + else error->all(FLERR,"Illegal compute pressure command"); + ++iarg; + } + } + + // Error check + + if (pairflag) { + if (force->pair) { + if (force->pair->born_matrix_enable == 0) { + if (comm->me == 0) error->warning(FLERR, "Pair style does not support compute born/matrix"); + } + } + } + if (bondflag) { + if (force->bond) { + if (force->bond->born_matrix_enable == 0) { + if (comm->me == 0) error->warning(FLERR, "Bond style does not support compute born/matrix"); + } + } + } + if (angleflag) { + if (force->angle) { + if (force->angle->born_matrix_enable == 0) { + if (comm->me == 0) error->warning(FLERR, "Angle style does not support compute born/matrix"); + } + } + } + if (dihedflag) { + if (force->dihedral) { + if (force->dihedral->born_matrix_enable == 0) { + if (comm->me == 0) error->warning(FLERR, "Dihedral style does not support compute born/matrix"); + } + } + } + if (impflag) { + if (force->improper) { + if (force->improper->born_matrix_enable == 0) { + if (comm->me == 0) error->warning(FLERR, "Improper style does not support compute born/matrix"); + } + } + } + if (kspaceflag) { + error->warning(FLERR, "KSPACE contribution not supported by compute born/matrix"); + } // Initialize some variables values_local = values_global = vector = nullptr; @@ -138,42 +222,16 @@ void ComputeBornMatrix::init() // Timestep Value dt = update->dt; - pairflag = 0; - bondflag = 0; - angleflag = 0; - dihedflag = 0; - impflag = 0; - kspaceflag = 0; + // IF everything works fine, + // this is to be removed + // + // if (force->pair) pairflag = 1; + // if (force->bond) bondflag = 1; + // if (force->angle) angleflag = 1; + // if (force->dihedral) dihedflag = 1; + // if (force->improper) impflag = 1; + // if (force->kspace) kspaceflag = 1; - if (force->pair) pairflag = 1; - if (force->bond) bondflag = 1; - if (force->angle) angleflag = 1; - if (force->dihedral) dihedflag = 1; - if (force->improper) impflag = 1; - if (force->kspace) kspaceflag = 1; - - // Error check - - if (comm->me == 0) { - if (pairflag && (force->pair->born_matrix_enable == 0)) { - error->all(FLERR, "Pair style does not support compute born/matrix"); - } - if (bondflag && (force->bond->born_matrix_enable == 0)) { - error->warning(FLERR, "Bond style does not support compute born/matrix"); - } - if (angleflag && (force->angle->born_matrix_enable == 0)) { - error->warning(FLERR, "Angle style does not support compute born/matrix"); - } - if (dihedflag && (force->dihedral->born_matrix_enable == 0)) { - error->warning(FLERR, "Dihedral style does not support compute born/matrix"); - } - if (impflag && (force->improper->born_matrix_enable == 0)) { - error->warning(FLERR, "Improper style does not support compute born/matrix"); - } - if (kspaceflag) { - error->warning(FLERR, "KSPACE contribution not supported by compute born/matrix"); - } - } // need an occasional half neighbor list int irequest = neighbor->request((void *) this); @@ -204,7 +262,6 @@ void ComputeBornMatrix::compute_vector() // Compute Born contribution if (pairflag) compute_pairs(); - // For now these functions are commented if (bondflag) compute_bonds(); if (angleflag) compute_angles(); if (dihedflag) compute_dihedrals(); @@ -331,11 +388,9 @@ void ComputeBornMatrix::compute_pairs() all atoms in interaction must be known to proc if bond is deleted or turned off (type <= 0) do not count or count contribution - COMMENTED FOR NOW ---------------------------------------------------------------------- */ void ComputeBornMatrix::compute_bonds() { - /* ---------------------------------------------------------------------- int i,m,n,nb,atom1,atom2,imol,iatom,btype,ivar; tagint tagprev; double dx,dy,dz,rsq; @@ -426,7 +481,6 @@ void ComputeBornMatrix::compute_bonds() } m += 21; } -------------------------------------------------------------------------- */ } /* ---------------------------------------------------------------------- @@ -436,11 +490,9 @@ void ComputeBornMatrix::compute_bonds() all atoms in interaction must be known to proc if bond is deleted or turned off (type <= 0) do not count or count contribution - COMMENTED FOR NOW ---------------------------------------------------------------------- */ void ComputeBornMatrix::compute_angles() { - /* ---------------------------------------------------------------------- int i,m,n,na,atom1,atom2,atom3,imol,iatom,atype,ivar; tagint tagprev; double delx1,dely1,delz1,delx2,dely2,delz2; @@ -568,8 +620,8 @@ void ComputeBornMatrix::compute_angles() c = 0; d = 0; for (i = 0; i<6; i++) { - a = albe[i][0]; - b = albe[i][1]; + a = sigma_albe[i][0]; + b = sigma_albe[i][1]; dcos[i] = cost*(del1[a]*del2[b]+del1[b]*del2[a]*r1r2inv - del1[a]*del1[b]*rsq1inv - del2[a]*del2[b]*rsq2inv); } @@ -578,8 +630,8 @@ void ComputeBornMatrix::compute_angles() b = albemunu[i][1]; c = albemunu[i][2]; d = albemunu[i][3]; - e = albe[i][0]; - f = albe[i][1]; + e = C_albe[i][0]; + f = C_albe[i][1]; d2lncos[i] = 2*(del1[a]*del1[b]*del1[c]*del1[d]*rsq1inv*rsq1inv + del2[a]*del2[b]*del2[c]*del2[d]*rsq2inv*rsq2inv) - (del1[a]*del2[b]+del1[b]*del2[a]) * @@ -592,7 +644,6 @@ void ComputeBornMatrix::compute_angles() } m+=21; } -------------------------------------------------------------------------- */ } /* ---------------------------------------------------------------------- @@ -601,12 +652,10 @@ void ComputeBornMatrix::compute_angles() all atoms in interaction must be in group all atoms in interaction must be known to proc if flag is set, compute requested info about dihedral - COMMENTED FOR NOW ------------------------------------------------------------------------- */ void ComputeBornMatrix::compute_dihedrals() { - /* ---------------------------------------------------------------------- int i,m,n,nd,atom1,atom2,atom3,atom4,imol,iatom,dtype,ivar; tagint tagprev; double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; @@ -781,8 +830,8 @@ void ComputeBornMatrix::compute_dihedrals() e = 0; f = 0; for (i = 0; i<6; i++) { - a = albe[i][0]; - b = albe[i][1]; + a = sigma_albe[i][0]; + b = sigma_albe[i][1]; dmm[i] = 2*(b2sq*b1[a]*b1[b]+b1sq*b2[a]*b2[b] - b1b2*(b1[a]*b2[b]+b1[b]*b2[a])); dnn[i] = 2*(b3sq*b2[a]*b2[b]+b2sq*b3[a]*b3[b] - b2b3*(b2[a]*b3[b]+b2[b]*b3[a])); dmn[i] = b1b2*(b2[a]*b3[b]+b2[b]*b3[a]) + b2b3*(b1[a]*b2[b]+b1[b]*b2[a]) @@ -794,8 +843,8 @@ void ComputeBornMatrix::compute_dihedrals() b = albemunu[i][1]; c = albemunu[i][2]; d = albemunu[i][3]; - e = albe[i][0]; - f = albe[i][1]; + e = C_albe[i][0]; + f = C_albe[i][1]; d2mm[i] = 4*(b1[a]*b1[b]*b2[c]*b2[d] + b1[c]*b1[d]*b2[a]*b2[b]) - 8*(b1[a]*b2[b]+b1[b]*b2[a])*(b1[c]*b2[d]+b1[d]*b2[c]); d2nn[i] = 4*(b2[a]*b2[b]*b3[c]*b3[d] + b2[c]*b2[d]*b3[a]*b3[b]) @@ -810,12 +859,11 @@ void ComputeBornMatrix::compute_dihedrals() + ra2inv*ra2inv*dmm[e]*dmm[f] - ra2inv*d2mm[i] + rb2inv*rb2inv*dnn[e]*dnn[f] - - rb2inv*d2nn[i] ); + - rb2inv*d2nn[i]); values_local[m+i] += dudih*d2cos[i] + du2dih*dcos[e]*dcos[f]; } } } m+=21; } -------------------------------------------------------------------------- */ } diff --git a/src/EXTRA-MOLECULE/dihedral_nharmonic.cpp b/src/EXTRA-MOLECULE/dihedral_nharmonic.cpp index 41a3d66ba6..be2cb2fcc5 100644 --- a/src/EXTRA-MOLECULE/dihedral_nharmonic.cpp +++ b/src/EXTRA-MOLECULE/dihedral_nharmonic.cpp @@ -39,6 +39,7 @@ DihedralNHarmonic::DihedralNHarmonic(LAMMPS *lmp) : Dihedral(lmp) { writedata = 1; a = nullptr; + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -336,3 +337,63 @@ void DihedralNHarmonic::write_data(FILE *fp) } } + +/* ----------------------------------------------------------------------*/ + +void DihedralNHarmonic::born_matrix(int nd, int i1, int i2, int i3, int i4, + double &dudih, double &du2dih) { + int i,type; + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; + double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,rginv,ra2inv,rb2inv,rabinv; + double c,s,kf; + + int **dihedrallist = neighbor->dihedrallist; + double **x = atom->x; + + type = dihedrallist[nd][4]; + + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + + // c,s calculation + + ax = vb1y*vb2zm - vb1z*vb2ym; + ay = vb1z*vb2xm - vb1x*vb2zm; + az = vb1x*vb2ym - vb1y*vb2xm; + bx = vb3y*vb2zm - vb3z*vb2ym; + by = vb3z*vb2xm - vb3x*vb2zm; + bz = vb3x*vb2ym - vb3y*vb2xm; + + rasq = ax*ax + ay*ay + az*az; + rbsq = bx*bx + by*by + bz*bz; + + ra2inv = rb2inv = 0.0; + if (rasq > 0) ra2inv = 1.0/rasq; + if (rbsq > 0) rb2inv = 1.0/rbsq; + rabinv = sqrt(ra2inv*rb2inv); + + c = (ax*bx + ay*by + az*bz)*rabinv; + + dudih = 0.0; + du2dih = 0.0; + for (i = 1; i < nterms[type]; i++) { + dudih += this->a[type][i]*i*pow(c,i-1); + } + + for (i = 2; i < nterms[type]; i++) { + du2dih += this->a[type][i]*i*(i-1)*pow(c, i-2); + } +} diff --git a/src/EXTRA-MOLECULE/dihedral_nharmonic.h b/src/EXTRA-MOLECULE/dihedral_nharmonic.h index ff41bd7ef5..c8a6245857 100644 --- a/src/EXTRA-MOLECULE/dihedral_nharmonic.h +++ b/src/EXTRA-MOLECULE/dihedral_nharmonic.h @@ -33,6 +33,7 @@ class DihedralNHarmonic : public Dihedral { void write_restart(FILE *); void read_restart(FILE *); void write_data(FILE *); + void born_matrix(int /*dtype*/, int, int, int, int, double&, double&); protected: int *nterms; diff --git a/src/MOLECULE/angle_cosine_squared.cpp b/src/MOLECULE/angle_cosine_squared.cpp index 4351095e56..5fe4768f6e 100644 --- a/src/MOLECULE/angle_cosine_squared.cpp +++ b/src/MOLECULE/angle_cosine_squared.cpp @@ -40,6 +40,7 @@ AngleCosineSquared::AngleCosineSquared(LAMMPS *lmp) : Angle(lmp) { k = nullptr; theta0 = nullptr; + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -263,3 +264,32 @@ double AngleCosineSquared::single(int type, int i1, int i2, int i3) double tk = k[type] * dcostheta; return tk*dcostheta; } + +/* ---------------------------------------------------------------------- */ + +void AngleCosineSquared::born_matrix(int type, int i1, int i2, int i3, double& du, double& du2) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + double dcostheta = c - cos(theta0[type]); + double tk = k[type] * dcostheta; + du2 = 2*k[type]; + du = du2*dcostheta; +} diff --git a/src/MOLECULE/angle_cosine_squared.h b/src/MOLECULE/angle_cosine_squared.h index d36aac64dd..58c13d519b 100644 --- a/src/MOLECULE/angle_cosine_squared.h +++ b/src/MOLECULE/angle_cosine_squared.h @@ -35,6 +35,7 @@ class AngleCosineSquared : public Angle { void read_restart(FILE *); void write_data(FILE *); virtual double single(int, int, int, int); + void born_matrix(int type, int i1, int i2, int i3, double& du, double& du2); protected: double *k, *theta0; diff --git a/src/MOLECULE/bond_harmonic.cpp b/src/MOLECULE/bond_harmonic.cpp index 32988d8d0a..39a9a41334 100644 --- a/src/MOLECULE/bond_harmonic.cpp +++ b/src/MOLECULE/bond_harmonic.cpp @@ -31,6 +31,7 @@ using namespace LAMMPS_NS; BondHarmonic::BondHarmonic(LAMMPS *lmp) : Bond(lmp) { reinitflag = 1; + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -201,6 +202,20 @@ double BondHarmonic::single(int type, double rsq, int /*i*/, int /*j*/, return rk*dr; } + +/* ---------------------------------------------------------------------- */ + +void BondHarmonic::born_matrix(int type, double rsq, int /*i*/, int /*j*/, + double &du, double& du2) +{ + double r = sqrt(rsq); + double dr = r - r0[type]; + du2 = 0.0; + du = 0.0; + du2 = 2*k[type]; + if (r > 0.0) du = du2*dr; +} + /* ---------------------------------------------------------------------- Return ptr to internal members upon request. ------------------------------------------------------------------------ */ diff --git a/src/MOLECULE/bond_harmonic.h b/src/MOLECULE/bond_harmonic.h index eb6c5e062e..7e9f6807ec 100644 --- a/src/MOLECULE/bond_harmonic.h +++ b/src/MOLECULE/bond_harmonic.h @@ -35,6 +35,7 @@ class BondHarmonic : public Bond { virtual void read_restart(FILE *); void write_data(FILE *); double single(int, double, int, int, double &); + void born_matrix(int, double, int, int, double &, double &); virtual void *extract(const char *, int &); protected: