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

Merge pull request #3668 from evoyiatzis/master

Browse Source 
Implementation of analytical expressions for Born matrix
This commit is contained in:
Axel Kohlmeyer
2023-03-13 17:19:05 -04:00
committed by GitHub
parent 7b07fe04c8 d5f6c7e0be
commit 6bec23e7db
29 changed files with 528 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/bond_fene.rst
View File

@ -51,7 +51,7 @@ in the same form as in pair style :doc:`nm/cut <pair_nm>`. 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

32
src/CLASS2/pair_lj_class2_coul_cut.cpp
View File

@ -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;

1
src/CLASS2/pair_lj_class2_coul_cut.h
View File

@ -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:

26
src/EXTRA-MOLECULE/bond_fene_nm.cpp
View File

@ -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;

1
src/EXTRA-MOLECULE/bond_fene_nm.h
View File

@ -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:

20
src/EXTRA-MOLECULE/bond_nonlinear.cpp
View File

@ -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;

1
src/EXTRA-MOLECULE/bond_nonlinear.h
View File

@ -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:

20
src/EXTRA-PAIR/pair_nm_cut.cpp
View File

@ -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;

1
src/EXTRA-PAIR/pair_nm_cut.h
View File

@ -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:

33
src/EXTRA-PAIR/pair_nm_cut_coul_cut.cpp
View File

@ -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;

1
src/EXTRA-PAIR/pair_nm_cut_coul_cut.h
View File

@ -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:

30
src/MOLECULE/angle_harmonic.cpp
View File

@ -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
------------------------------------------------------------------------ */

1
src/MOLECULE/angle_harmonic.h
View File

@ -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:

29
src/MOLECULE/bond_fene.cpp
View File

@ -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;

3
src/MOLECULE/bond_fene.h
View File

@ -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:

17
src/MOLECULE/bond_morse.cpp
View File

@ -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;

1
src/MOLECULE/bond_morse.h
View File

@ -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:

85
src/MOLECULE/dihedral_multi_harmonic.cpp
View File

@ -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);
}

1
src/MOLECULE/dihedral_multi_harmonic.h
View File

@ -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;

99
src/MOLECULE/dihedral_opls.cpp
View File

@ -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);
}

1
src/MOLECULE/dihedral_opls.h
View File

@ -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;

38
src/pair_buck_coul_cut.cpp
View File

@ -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;

1
src/pair_buck_coul_cut.h
View File

@ -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:

26
src/pair_coul_cut.cpp
View File

@ -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;

1
src/pair_coul_cut.h
View File

@ -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:

31
src/pair_coul_debye.cpp
View File

@ -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;
}

1
src/pair_coul_debye.h
View File

@ -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;

30
src/pair_lj_cut_coul_cut.cpp
View File

@ -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;

1
src/pair_lj_cut_coul_cut.h
View File

@ -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: