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

Implement born_matrix in bond_fene.cpp

Browse Source
This commit is contained in:
Evangelos Voyiatzis
2023-02-23 18:15:43 +02:00
committed by GitHub
parent 5fb11e3f06
commit bb75ed5071
1 changed files with 29 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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;