diff --git a/doc/src/bond_fene.rst b/doc/src/bond_fene.rst index be7775489a..24b3a4e1e2 100644 --- a/doc/src/bond_fene.rst +++ b/doc/src/bond_fene.rst @@ -51,7 +51,7 @@ in the same form as in pair style :doc:`nm/cut `. The bond energy is th .. math:: - E = -0.5 K r_0^2 \ln \left[ 1 - \left(\frac{r}{R_0}\right)^2\right] + \frac{E_0}{(n-m)} \left[ m \left(\frac{r_0}{r}\right)^n - n \left(\frac{r_0}{r}\right)^m \right] + E = -0.5 K R_0^2 \ln \left[ 1 - \left(\frac{r}{R_0}\right)^2\right] + \frac{E_0}{(n-m)} \left[ m \left(\frac{r_0}{r}\right)^n - n \left(\frac{r_0}{r}\right)^m \right] Similar to the *fene* style, the generalized Lennard-Jones is cut off at the potential minimum, :math:`r_0`, to be repulsive only. The following diff --git a/src/CLASS2/pair_lj_class2_coul_cut.cpp b/src/CLASS2/pair_lj_class2_coul_cut.cpp index d05eb9acae..3f66f1550f 100644 --- a/src/CLASS2/pair_lj_class2_coul_cut.cpp +++ b/src/CLASS2/pair_lj_class2_coul_cut.cpp @@ -32,6 +32,7 @@ using namespace MathConst; PairLJClass2CoulCut::PairLJClass2CoulCut(LAMMPS *lmp) : Pair(lmp) { + born_matrix_enable = 1; writedata = 1; centroidstressflag = CENTROID_SAME; } @@ -470,6 +471,37 @@ double PairLJClass2CoulCut::single(int i, int j, int itype, int jtype, double rs /* ---------------------------------------------------------------------- */ +void PairLJClass2CoulCut::born_matrix(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &dupair, + double &du2pair) +{ + double rinv, r2inv, r3inv, r7inv, r8inv; + double du_lj, du2_lj, du_coul, du2_coul; + + double *q = atom->q; + double qqrd2e = force->qqrd2e; + + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r7inv = r3inv * r3inv * rinv; + r8inv = r7inv * rinv; + + // Reminder: lj1[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 9.0); + // Reminder: lj2[i][j] = 18.0 * epsilon[i][j] * pow(sigma[i][j], 6.0); + du_lj = r7inv * (lj2[itype][jtype] - lj1[itype][jtype] * r3inv); + du2_lj = r8inv * (10 * lj1[itype][jtype] * r3inv - 7 * lj2[itype][jtype]); + + // Reminder: qqrd2e converts q^2/r to energy w/ dielectric constant + du_coul = -qqrd2e * q[i] * q[j] * r2inv; + du2_coul = 2.0 * qqrd2e * q[i] * q[j] * r3inv; + + dupair = factor_lj * du_lj + factor_coul * du_coul; + du2pair = factor_lj * du2_lj + factor_coul * du2_coul; +} + +/* ---------------------------------------------------------------------- */ + void *PairLJClass2CoulCut::extract(const char *str, int &dim) { dim = 2; diff --git a/src/CLASS2/pair_lj_class2_coul_cut.h b/src/CLASS2/pair_lj_class2_coul_cut.h index ca5d6830f5..c4e9ba2ec9 100644 --- a/src/CLASS2/pair_lj_class2_coul_cut.h +++ b/src/CLASS2/pair_lj_class2_coul_cut.h @@ -40,6 +40,7 @@ class PairLJClass2CoulCut : public Pair { void write_data(FILE *) override; void write_data_all(FILE *) override; double single(int, int, int, int, double, double, double, double &) override; + void born_matrix(int, int, int, int, double, double, double, double &, double &) override; void *extract(const char *, int &) override; protected: diff --git a/src/EXTRA-MOLECULE/bond_fene_nm.cpp b/src/EXTRA-MOLECULE/bond_fene_nm.cpp index 5451479881..59f60379bd 100644 --- a/src/EXTRA-MOLECULE/bond_fene_nm.cpp +++ b/src/EXTRA-MOLECULE/bond_fene_nm.cpp @@ -28,7 +28,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondFENENM::BondFENENM(LAMMPS *lmp) : BondFENE(lmp), nn(nullptr), mm(nullptr) {} +BondFENENM::BondFENENM(LAMMPS *lmp) : BondFENE(lmp), nn(nullptr), mm(nullptr) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -270,6 +273,27 @@ double BondFENENM::single(int type, double rsq, int /*i*/, int /*j*/, double &ff /* ---------------------------------------------------------------------- */ +void BondFENENM::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du, double &du2) +{ + double r = sqrt(rsq); + double r0sq = r0[type] * r0[type]; + double rlogarg = 1.0 - rsq / r0sq; + + // Contribution from the attractive term + du = k[type] * r / rlogarg; + du2 = k[type] * (1.0 + rsq / r0sq) / (rlogarg * rlogarg); + + // Contribution from the repulsive Lennard-Jones term + if (rsq < sigma[type] * sigma[type]) { + double prefactor = epsilon[type] * nn[type] * mm[type] / (nn[type] - mm[type]); + du += prefactor * (pow(sigma[type] / r, mm[type]) - pow(sigma[type] / r, nn[type])) / r; + du2 += prefactor * ((nn[type] + 1.0) * pow(sigma[type] / r, nn[type]) - + (mm[type] + 1.0) * pow(sigma[type] / r, mm[type])) / rsq; + } +} + +/* ---------------------------------------------------------------------- */ + void *BondFENENM::extract(const char *str, int &dim) { dim = 1; diff --git a/src/EXTRA-MOLECULE/bond_fene_nm.h b/src/EXTRA-MOLECULE/bond_fene_nm.h index 07da20e425..4e9899648b 100644 --- a/src/EXTRA-MOLECULE/bond_fene_nm.h +++ b/src/EXTRA-MOLECULE/bond_fene_nm.h @@ -35,6 +35,7 @@ class BondFENENM : public BondFENE { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, double, int, int, double &) override; + void born_matrix(int, double, int, int, double &, double &) override; void *extract(const char *, int &) override; protected: diff --git a/src/EXTRA-MOLECULE/bond_nonlinear.cpp b/src/EXTRA-MOLECULE/bond_nonlinear.cpp index 873c54e917..a2955b7d2e 100644 --- a/src/EXTRA-MOLECULE/bond_nonlinear.cpp +++ b/src/EXTRA-MOLECULE/bond_nonlinear.cpp @@ -28,7 +28,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondNonlinear::BondNonlinear(LAMMPS *lmp) : Bond(lmp) {} +BondNonlinear::BondNonlinear(LAMMPS *lmp) : Bond(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -207,6 +210,21 @@ double BondNonlinear::single(int type, double rsq, int /*i*/, int /*j*/, /* ---------------------------------------------------------------------- */ +void BondNonlinear::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du, double &du2) +{ + double r = sqrt(rsq); + double dr = r - r0[type]; + double drsq = dr * dr; + double lamdasq = lamda[type] * lamda[type]; + double denom = lamdasq - drsq; + double denomsq = denom * denom; + + du = 2.0 * epsilon[type] * lamdasq * dr / denomsq; + du2 = 2.0 * epsilon[type] * lamdasq * (lamdasq + 3.0 * drsq)/ (denomsq * denom); +} + +/* ---------------------------------------------------------------------- */ + void *BondNonlinear::extract(const char *str, int &dim) { dim = 1; diff --git a/src/EXTRA-MOLECULE/bond_nonlinear.h b/src/EXTRA-MOLECULE/bond_nonlinear.h index 3023add89d..f092643662 100644 --- a/src/EXTRA-MOLECULE/bond_nonlinear.h +++ b/src/EXTRA-MOLECULE/bond_nonlinear.h @@ -35,6 +35,7 @@ class BondNonlinear : public Bond { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, double, int, int, double &) override; + void born_matrix(int, double, int, int, double &, double &) override; void *extract(const char *, int &) override; protected: diff --git a/src/EXTRA-PAIR/pair_nm_cut.cpp b/src/EXTRA-PAIR/pair_nm_cut.cpp index 18b5810abc..3156ec224d 100644 --- a/src/EXTRA-PAIR/pair_nm_cut.cpp +++ b/src/EXTRA-PAIR/pair_nm_cut.cpp @@ -36,6 +36,7 @@ using namespace MathConst; PairNMCut::PairNMCut(LAMMPS *lmp) : Pair(lmp) { + born_matrix_enable = 1; writedata = 1; } @@ -416,6 +417,25 @@ double PairNMCut::single(int /*i*/, int /*j*/, int itype, int jtype, /* ---------------------------------------------------------------------- */ +void PairNMCut::born_matrix(int /*i*/, int /*j*/, int itype, int jtype, double rsq, + double /*factor_coul*/, double factor_lj, double &dupair, + double &du2pair) +{ + double r = sqrt(rsq); + double prefactor = e0nm[itype][jtype]*nm[itype][jtype]; + + double du = prefactor * + (r0m[itype][jtype]/pow(r,mm[itype][jtype]) - r0n[itype][jtype]/pow(r,nn[itype][jtype])) / r; + double du2 = prefactor * + (r0n[itype][jtype]*(nn[itype][jtype] + 1.0) / pow(r,nn[itype][jtype]) - + r0m[itype][jtype]*(mm[itype][jtype] + 1.0) / pow(r,mm[itype][jtype])) / rsq; + + dupair = factor_lj * du; + du2pair = factor_lj * du2; +} + +/* ---------------------------------------------------------------------- */ + void *PairNMCut::extract(const char *str, int &dim) { dim = 2; diff --git a/src/EXTRA-PAIR/pair_nm_cut.h b/src/EXTRA-PAIR/pair_nm_cut.h index ab91747252..a9a800af23 100644 --- a/src/EXTRA-PAIR/pair_nm_cut.h +++ b/src/EXTRA-PAIR/pair_nm_cut.h @@ -40,6 +40,7 @@ class PairNMCut : public Pair { void write_data(FILE *) override; void write_data_all(FILE *) override; double single(int, int, int, int, double, double, double, double &) override; + void born_matrix(int, int, int, int, double, double, double, double &, double &) override; void *extract(const char *, int &) override; protected: diff --git a/src/EXTRA-PAIR/pair_nm_cut_coul_cut.cpp b/src/EXTRA-PAIR/pair_nm_cut_coul_cut.cpp index adc6d5a058..531f0e615f 100644 --- a/src/EXTRA-PAIR/pair_nm_cut_coul_cut.cpp +++ b/src/EXTRA-PAIR/pair_nm_cut_coul_cut.cpp @@ -37,6 +37,7 @@ using namespace MathConst; PairNMCutCoulCut::PairNMCutCoulCut(LAMMPS *lmp) : Pair(lmp) { + born_matrix_enable = 1; writedata = 1; } @@ -481,6 +482,38 @@ double PairNMCutCoulCut::single(int i, int j, int itype, int jtype, /* ---------------------------------------------------------------------- */ +void PairNMCutCoulCut::born_matrix(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &dupair, + double &du2pair) +{ + double r, rinv, r2inv, r3inv; + double du_lj, du2_lj, du_coul, du2_coul; + + double *q = atom->q; + double qqrd2e = force->qqrd2e; + + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r = sqrt(rsq); + + double prefactor = e0nm[itype][jtype]*nm[itype][jtype]; + + du_lj = prefactor * + (r0m[itype][jtype]/pow(r,mm[itype][jtype]) - r0n[itype][jtype]/pow(r,nn[itype][jtype])) / r; + du2_lj = prefactor * + (r0n[itype][jtype]*(nn[itype][jtype] + 1.0) / pow(r,nn[itype][jtype]) - + r0m[itype][jtype]*(mm[itype][jtype] + 1.0) / pow(r,mm[itype][jtype])) / rsq; + + du_coul = -qqrd2e * q[i] * q[j] * r2inv; + du2_coul = 2.0 * qqrd2e * q[i] * q[j] * r3inv; + + dupair = factor_lj * du_lj + factor_coul * du_coul; + du2pair = factor_lj * du2_lj + factor_coul * du2_coul; +} + +/* ---------------------------------------------------------------------- */ + void *PairNMCutCoulCut::extract(const char *str, int &dim) { dim = 2; diff --git a/src/EXTRA-PAIR/pair_nm_cut_coul_cut.h b/src/EXTRA-PAIR/pair_nm_cut_coul_cut.h index ff51920c63..f2254f7fd3 100644 --- a/src/EXTRA-PAIR/pair_nm_cut_coul_cut.h +++ b/src/EXTRA-PAIR/pair_nm_cut_coul_cut.h @@ -41,6 +41,7 @@ class PairNMCutCoulCut : public Pair { void write_data(FILE *) override; void write_data_all(FILE *) override; double single(int, int, int, int, double, double, double, double &) override; + void born_matrix(int, int, int, int, double, double, double, double &, double &) override; void *extract(const char *, int &) override; protected: diff --git a/src/MOLECULE/angle_harmonic.cpp b/src/MOLECULE/angle_harmonic.cpp index 4d0683b9bc..e9f1c528ef 100644 --- a/src/MOLECULE/angle_harmonic.cpp +++ b/src/MOLECULE/angle_harmonic.cpp @@ -35,6 +35,7 @@ static constexpr double SMALL = 0.001; AngleHarmonic::AngleHarmonic(LAMMPS *_lmp) : Angle(_lmp) { + born_matrix_enable = 1; k = nullptr; theta0 = nullptr; } @@ -266,6 +267,35 @@ double AngleHarmonic::single(int type, int i1, int i2, int i3) return tk * dtheta; } +/* ---------------------------------------------------------------------- */ + +void AngleHarmonic::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 theta = acos(c); + + double dtheta = theta - theta0[type]; + du = -2 * k[type] * dtheta / sin(theta); + du2 = 2 * k[type] * (sin(theta) - dtheta * cos(theta)) / pow(sin(theta), 3); +} + /* ---------------------------------------------------------------------- return ptr to internal members upon request ------------------------------------------------------------------------ */ diff --git a/src/MOLECULE/angle_harmonic.h b/src/MOLECULE/angle_harmonic.h index c50df2ab4e..284a5316ef 100644 --- a/src/MOLECULE/angle_harmonic.h +++ b/src/MOLECULE/angle_harmonic.h @@ -35,6 +35,7 @@ class AngleHarmonic : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override; void *extract(const char *, int &) override; protected: diff --git a/src/MOLECULE/bond_fene.cpp b/src/MOLECULE/bond_fene.cpp index 0f279fccb6..7f55c89296 100644 --- a/src/MOLECULE/bond_fene.cpp +++ b/src/MOLECULE/bond_fene.cpp @@ -30,6 +30,13 @@ using MathConst::MY_CUBEROOT2; /* ---------------------------------------------------------------------- */ +BondFENE::BondFENE(LAMMPS *_lmp) : Bond(_lmp) +{ + born_matrix_enable = 1; +} + +/* ---------------------------------------------------------------------- */ + BondFENE::~BondFENE() { if (allocated && !copymode) { @@ -262,6 +269,28 @@ double BondFENE::single(int type, double rsq, int /*i*/, int /*j*/, double &ffor /* ---------------------------------------------------------------------- */ +void BondFENE::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du, double &du2) +{ + double r = sqrt(rsq); + double r0sq = r0[type] * r0[type]; + double rlogarg = 1.0 - rsq / r0sq; + + // Contribution from the attractive term + du = k[type] * r / rlogarg; + du2 = k[type] * (1.0 + rsq / r0sq) / (rlogarg * rlogarg); + + // Contribution from the repulsive Lennard-Jones term + if (rsq < MY_CUBEROOT2 * sigma[type] * sigma[type]) { + double sr2 = sigma[type] * sigma[type] / rsq; + double sr6 = sr2 * sr2 * sr2; + + du += 48.0 * epsilon[type] * sr6 * (0.5 - sr6) / r; + du2 += 48.0 * epsilon[type] * sr6 * (13.0 * sr6 - 3.5) / rsq; + } +} + +/* ---------------------------------------------------------------------- */ + void *BondFENE::extract(const char *str, int &dim) { dim = 1; diff --git a/src/MOLECULE/bond_fene.h b/src/MOLECULE/bond_fene.h index eb9e981a32..1628e7a2a2 100644 --- a/src/MOLECULE/bond_fene.h +++ b/src/MOLECULE/bond_fene.h @@ -26,7 +26,7 @@ namespace LAMMPS_NS { class BondFENE : public Bond { public: - BondFENE(class LAMMPS *_lmp) : Bond(_lmp) {} + BondFENE(class LAMMPS *); ~BondFENE() override; void compute(int, int) override; void coeff(int, char **) override; @@ -36,6 +36,7 @@ class BondFENE : public Bond { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, double, int, int, double &) override; + void born_matrix(int, double, int, int, double &, double &) override; void *extract(const char *, int &) override; protected: diff --git a/src/MOLECULE/bond_morse.cpp b/src/MOLECULE/bond_morse.cpp index e9555698da..b7020b0441 100644 --- a/src/MOLECULE/bond_morse.cpp +++ b/src/MOLECULE/bond_morse.cpp @@ -30,7 +30,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondMorse::BondMorse(LAMMPS *_lmp) : Bond(_lmp) {} +BondMorse::BondMorse(LAMMPS *_lmp) : Bond(_lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -209,6 +212,18 @@ double BondMorse::single(int type, double rsq, int /*i*/, int /*j*/, double &ffo /* ---------------------------------------------------------------------- */ +void BondMorse::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du, double &du2) +{ + double r = sqrt(rsq); + double dr = r - r0[type]; + double ralpha = exp(-alpha[type] * dr); + + du = 2.0 * d0[type] * alpha[type] * (1.0 - ralpha) * ralpha; + du2 = -2.0 * d0[type] * alpha[type] * alpha[type] * (1.0 - 2.0 * ralpha) * ralpha; +} + +/* ---------------------------------------------------------------------- */ + void *BondMorse::extract(const char *str, int &dim) { dim = 1; diff --git a/src/MOLECULE/bond_morse.h b/src/MOLECULE/bond_morse.h index 5515ea023e..4d0e8e359c 100644 --- a/src/MOLECULE/bond_morse.h +++ b/src/MOLECULE/bond_morse.h @@ -35,6 +35,7 @@ class BondMorse : public Bond { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, double, int, int, double &) override; + void born_matrix(int, double, int, int, double &, double &) override; void *extract(const char *, int &) override; protected: diff --git a/src/MOLECULE/dihedral_multi_harmonic.cpp b/src/MOLECULE/dihedral_multi_harmonic.cpp index 11213b59cd..8e6685cac9 100644 --- a/src/MOLECULE/dihedral_multi_harmonic.cpp +++ b/src/MOLECULE/dihedral_multi_harmonic.cpp @@ -36,6 +36,7 @@ static constexpr double SMALL = 0.001; DihedralMultiHarmonic::DihedralMultiHarmonic(LAMMPS *_lmp) : Dihedral(_lmp) { writedata = 1; + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -322,3 +323,87 @@ void DihedralMultiHarmonic::write_data(FILE *fp) for (int i = 1; i <= atom->ndihedraltypes; i++) fprintf(fp, "%d %g %g %g %g %g\n", i, a1[i], a2[i], a3[i], a4[i], a5[i]); } + +/* ---------------------------------------------------------------------- */ + +void DihedralMultiHarmonic::born_matrix(int nd, int i1, int i2, int i3, int i4, + double &du, double &du2) +{ + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + double sb1, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2; + double b2mag, b3mag2, b3mag, ctmp, r12c1, c1mag, r12c2; + double c2mag, sc1, sc2, s12, c; + double sin2; + + double **x = atom->x; + int **dihedrallist = neighbor->dihedrallist; + + int type = dihedrallist[nd][4]; + + // 1st bond + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + // 2nd bond + 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; + + // 3rd bond + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + + // c0 calculation + sb1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z); + sb3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3; + + // 1st and 2nd angle + b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x * vb2x + vb1y * vb2y + vb1z * vb2z; + r12c1 = 1.0 / (b1mag * b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm * vb3x + vb2ym * vb3y + vb2zm * vb3z; + r12c2 = 1.0 / (b2mag * b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + sin2 = MAX(1.0 - c1mag * c1mag, 0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0 / sc1; + + sin2 = MAX(1.0 - c2mag * c2mag, 0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0 / sc2; + + s12 = sc1 * sc2; + c = (c0 + c1mag * c2mag) * s12; + + // error check + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4); + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + du = a2[type] + c * (2.0 * a3[type] + c * (3.0 * a4[type] + c * 4.0 * a5[type])); + du2 = 2.0 * a3[type] + 6.0 * c * (a4[type] + 2.0 * a5[type] * c); +} diff --git a/src/MOLECULE/dihedral_multi_harmonic.h b/src/MOLECULE/dihedral_multi_harmonic.h index 4780028181..45f1112b4f 100644 --- a/src/MOLECULE/dihedral_multi_harmonic.h +++ b/src/MOLECULE/dihedral_multi_harmonic.h @@ -33,6 +33,7 @@ class DihedralMultiHarmonic : public Dihedral { void write_restart(FILE *) override; void read_restart(FILE *) override; void write_data(FILE *) override; + void born_matrix(int, int, int, int, int, double &, double &) override; protected: double *a1, *a2, *a3, *a4, *a5; diff --git a/src/MOLECULE/dihedral_opls.cpp b/src/MOLECULE/dihedral_opls.cpp index af34cd24db..eced454d68 100644 --- a/src/MOLECULE/dihedral_opls.cpp +++ b/src/MOLECULE/dihedral_opls.cpp @@ -37,6 +37,7 @@ static constexpr double SMALLER = 0.00001; DihedralOPLS::DihedralOPLS(LAMMPS *_lmp) : Dihedral(_lmp) { writedata = 1; + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -332,3 +333,101 @@ void DihedralOPLS::write_data(FILE *fp) for (int i = 1; i <= atom->ndihedraltypes; i++) fprintf(fp, "%d %g %g %g %g\n", i, 2.0 * k1[i], 2.0 * k2[i], 2.0 * k3[i], 2.0 * k4[i]); } + +/* ----------------------------------------------------------------------*/ + +void DihedralOPLS::born_matrix(int nd, int i1, int i2, int i3, int i4, + double &du, double &du2) +{ + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + double sb1, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2; + double b2mag, b3mag2, b3mag, ctmp, r12c1, c1mag, r12c2; + double c2mag, sc1, sc2, s12, c; + double cx, cy, cz, cmag, dx, phi, si, sin2; + + double **x = atom->x; + int **dihedrallist = neighbor->dihedrallist; + + int type = dihedrallist[nd][4]; + + // 1st bond + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + // 2nd bond + 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; + + // 3rd bond + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + + // c0 calculation + sb1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z); + sb3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3; + + // 1st and 2nd angle + b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x * vb2x + vb1y * vb2y + vb1z * vb2z; + r12c1 = 1.0 / (b1mag * b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm * vb3x + vb2ym * vb3y + vb2zm * vb3z; + r12c2 = 1.0 / (b2mag * b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + sin2 = MAX(1.0 - c1mag * c1mag, 0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0 / sc1; + + sin2 = MAX(1.0 - c2mag * c2mag, 0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0 / sc2; + + s12 = sc1 * sc2; + c = (c0 + c1mag * c2mag) * s12; + + cx = vb1y * vb2z - vb1z * vb2y; + cy = vb1z * vb2x - vb1x * vb2z; + cz = vb1x * vb2y - vb1y * vb2x; + cmag = sqrt(cx * cx + cy * cy + cz * cz); + dx = (cx * vb3x + cy * vb3y + cz * vb3z) / cmag / b3mag; + + // error check + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4); + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + phi = acos(c); + if (dx < 0.0) phi *= -1.0; + si = sin(phi); + if (fabs(si) < SMALLER) si = SMALLER; + + du = k1[type] - 2.0 * k2[type] * sin(2.0 * phi) / si + 3.0 * k3[type] * sin(3.0 * phi) / si + - 4.0 * k4[type] * sin(4.0 * phi) / si; + du2 = (4.0 * k2[type] * si * cos(2.0 * phi) - 2.0 * k2[type] * sin(2.0 * phi) + - 9.0 * k3[type] * si * cos(3.0 * phi) + 3.0 * k3[type] * sin(3.0 * phi) + + 16.0 * k4[type] * si * cos(4.0 * phi) - 4.0 * k4[type] * sin(4.0 * phi)) / (si * si * si); +} diff --git a/src/MOLECULE/dihedral_opls.h b/src/MOLECULE/dihedral_opls.h index c6286b3b6f..82180cf323 100644 --- a/src/MOLECULE/dihedral_opls.h +++ b/src/MOLECULE/dihedral_opls.h @@ -33,6 +33,7 @@ class DihedralOPLS : public Dihedral { void write_restart(FILE *) override; void read_restart(FILE *) override; void write_data(FILE *) override; + void born_matrix(int, int, int, int, int, double &, double &) override; protected: double *k1, *k2, *k3, *k4; diff --git a/src/pair_buck_coul_cut.cpp b/src/pair_buck_coul_cut.cpp index 4f85ccd430..369b77474f 100644 --- a/src/pair_buck_coul_cut.cpp +++ b/src/pair_buck_coul_cut.cpp @@ -36,6 +36,7 @@ using namespace MathConst; PairBuckCoulCut::PairBuckCoulCut(LAMMPS *lmp) : Pair(lmp) { + born_matrix_enable = 1; writedata = 1; } @@ -475,6 +476,43 @@ double PairBuckCoulCut::single(int i, int j, int itype, int jtype, double rsq, d /* ---------------------------------------------------------------------- */ +void PairBuckCoulCut::born_matrix(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &dupair, + double &du2pair) +{ + double rinv, r2inv, r3inv, r6inv, r7inv, r8inv, r, rexp; + double du_lj, du2_lj, du_coul, du2_coul; + + double *q = atom->q; + double qqrd2e = force->qqrd2e; + + r = sqrt(rsq); + rexp = exp(-r*rhoinv[itype][jtype]); + + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r2inv * r2inv * r2inv; + r7inv = r6inv * rinv; + r8inv = r6inv * r2inv; + + // Reminder: buck1[itype][jtype] = a[itype][jtype]/rho[itype][jtype]; + // Reminder: buck2[itype][jtype] = 6.0*c[itype][jtype]; + + du_lj = buck2[itype][jtype] * r7inv - buck1[itype][jtype] * rexp; + du2_lj = (buck1[itype][jtype] / rho[itype][jtype]) * rexp - 7 * buck2[itype][jtype] * r8inv; + + // Reminder: qqrd2e converts q^2/r to energy w/ dielectric constant + + du_coul = -qqrd2e * q[i] * q[j] * r2inv; + du2_coul = 2.0 * qqrd2e * q[i] * q[j] * r3inv; + + dupair = factor_lj * du_lj + factor_coul * du_coul; + du2pair = factor_lj * du2_lj + factor_coul * du2_coul; +} + +/* ---------------------------------------------------------------------- */ + void *PairBuckCoulCut::extract(const char *str, int &dim) { dim = 2; diff --git a/src/pair_buck_coul_cut.h b/src/pair_buck_coul_cut.h index b5f19def38..3b13cb46b3 100644 --- a/src/pair_buck_coul_cut.h +++ b/src/pair_buck_coul_cut.h @@ -40,6 +40,7 @@ class PairBuckCoulCut : public Pair { void write_data(FILE *) override; void write_data_all(FILE *) override; double single(int, int, int, int, double, double, double, double &) override; + void born_matrix(int, int, int, int, double, double, double, double &, double &) override; void *extract(const char *, int &) override; protected: diff --git a/src/pair_coul_cut.cpp b/src/pair_coul_cut.cpp index 9246e535fe..f0a749baa8 100644 --- a/src/pair_coul_cut.cpp +++ b/src/pair_coul_cut.cpp @@ -30,6 +30,7 @@ using namespace LAMMPS_NS; PairCoulCut::PairCoulCut(LAMMPS *lmp) : Pair(lmp) { + born_matrix_enable = 1; writedata = 1; } @@ -329,6 +330,31 @@ double PairCoulCut::single(int i, int j, int /*itype*/, int /*jtype*/, double rs /* ---------------------------------------------------------------------- */ +void PairCoulCut::born_matrix(int i, int j, int /*itype*/, int /*jtype*/, double rsq, + double factor_coul, double /*factor_lj*/, double &dupair, + double &du2pair) +{ + double rinv, r2inv, r3inv; + double du_coul, du2_coul; + + double *q = atom->q; + double qqrd2e = force->qqrd2e; + + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + + // Reminder: qqrd2e converts q^2/r to energy w/ dielectric constant + + du_coul = -qqrd2e * q[i] * q[j] * r2inv; + du2_coul = 2.0 * qqrd2e * q[i] * q[j] * r3inv; + + dupair = factor_coul * du_coul; + du2pair = factor_coul * du2_coul; +} + +/* ---------------------------------------------------------------------- */ + void *PairCoulCut::extract(const char *str, int &dim) { dim = 2; diff --git a/src/pair_coul_cut.h b/src/pair_coul_cut.h index dc3a3df64f..0890a7da33 100644 --- a/src/pair_coul_cut.h +++ b/src/pair_coul_cut.h @@ -40,6 +40,7 @@ class PairCoulCut : public Pair { void write_data(FILE *) override; void write_data_all(FILE *) override; double single(int, int, int, int, double, double, double, double &) override; + void born_matrix(int, int, int, int, double, double, double, double &, double &) override; void *extract(const char *, int &) override; protected: diff --git a/src/pair_coul_debye.cpp b/src/pair_coul_debye.cpp index 0f375f6e00..285b742074 100644 --- a/src/pair_coul_debye.cpp +++ b/src/pair_coul_debye.cpp @@ -26,7 +26,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairCoulDebye::PairCoulDebye(LAMMPS *lmp) : PairCoulCut(lmp) {} +PairCoulDebye::PairCoulDebye(LAMMPS *lmp) : PairCoulCut(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -177,3 +180,29 @@ double PairCoulDebye::single(int i, int j, int /*itype*/, int /*jtype*/, phicoul = force->qqrd2e * atom->q[i]*atom->q[j] * rinv * screening; return factor_coul*phicoul; } + +/* ---------------------------------------------------------------------- */ + +void PairCoulDebye::born_matrix(int i, int j, int /*itype*/, int /*jtype*/, double rsq, + double factor_coul, double /*factor_lj*/, double &dupair, + double &du2pair) +{ + double r, rinv, r2inv, r3inv, screening; + double du_coul, du2_coul; + + double *q = atom->q; + double qqrd2e = force->qqrd2e; + + r = sqrt(rsq); + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + screening = exp(-kappa*r); + + // Reminder: qqrd2e converts q^2/r to energy w/ dielectric constant + du_coul = -qqrd2e * q[i] * q[j] * r2inv * (1 + kappa * r) * screening; + du2_coul = qqrd2e * q[i] * q[j] * r3inv * (2 + 2 * kappa * r + kappa * kappa * rsq) * screening; + + dupair = factor_coul * du_coul; + du2pair = factor_coul * du2_coul; +} diff --git a/src/pair_coul_debye.h b/src/pair_coul_debye.h index 0f28efb46a..1d8b5760fa 100644 --- a/src/pair_coul_debye.h +++ b/src/pair_coul_debye.h @@ -32,6 +32,7 @@ class PairCoulDebye : public PairCoulCut { void write_restart_settings(FILE *) override; void read_restart_settings(FILE *) override; double single(int, int, int, int, double, double, double, double &) override; + void born_matrix(int, int, int, int, double, double, double, double &, double &) override; protected: double kappa; diff --git a/src/pair_lj_cut_coul_cut.cpp b/src/pair_lj_cut_coul_cut.cpp index bc5f5ed93e..5438e82377 100644 --- a/src/pair_lj_cut_coul_cut.cpp +++ b/src/pair_lj_cut_coul_cut.cpp @@ -32,6 +32,7 @@ using namespace MathConst; PairLJCutCoulCut::PairLJCutCoulCut(LAMMPS *lmp) : Pair(lmp) { + born_matrix_enable = 1; writedata = 1; } @@ -460,6 +461,35 @@ double PairLJCutCoulCut::single(int i, int j, int itype, int jtype, double rsq, /* ---------------------------------------------------------------------- */ +void PairLJCutCoulCut::born_matrix(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &dupair, + double &du2pair) +{ + double rinv, r2inv, r3inv, r6inv, du, du2; + double du_lj, du2_lj, du_coul, du2_coul; + + double *q = atom->q; + double qqrd2e = force->qqrd2e; + + r2inv = 1.0 / rsq; + rinv = sqrt(r2inv); + r3inv = r2inv * rinv; + r6inv = r2inv * r2inv * r2inv; + + // Reminder: lj1 = 48*e*s^12, lj2 = 24*e*s^6 + + du_lj = r6inv * rinv * (lj2[itype][jtype] - lj1[itype][jtype] * r6inv); + du2_lj = r6inv * r2inv * (13 * lj1[itype][jtype] * r6inv - 7 * lj2[itype][jtype]); + + du_coul = -qqrd2e * q[i] * q[j] * r2inv; + du2_coul = 2.0 * qqrd2e * q[i] * q[j] * r3inv; + + dupair = factor_lj * du_lj + factor_coul * du_coul; + du2pair = factor_lj * du2_lj + factor_coul * du2_coul; +} + +/* ---------------------------------------------------------------------- */ + void *PairLJCutCoulCut::extract(const char *str, int &dim) { dim = 2; diff --git a/src/pair_lj_cut_coul_cut.h b/src/pair_lj_cut_coul_cut.h index b4e1bfb5c9..fdf54d9340 100644 --- a/src/pair_lj_cut_coul_cut.h +++ b/src/pair_lj_cut_coul_cut.h @@ -40,6 +40,7 @@ class PairLJCutCoulCut : public Pair { void write_data(FILE *) override; void write_data_all(FILE *) override; double single(int, int, int, int, double, double, double, double &) override; + void born_matrix(int, int, int, int, double, double, double, double &, double &) override; void *extract(const char *, int &) override; protected: