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

Corrected energy for capped force potential

Browse Source
This commit is contained in:
Oliver Henrich
2024-05-15 15:24:36 +01:00
parent 99b3233564
commit 59ce8c966c
1 changed files with 26 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -167,6 +167,7 @@ void BondOxdnaFene::compute(int eflag, int vflag)
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
const double rlogarg_min = 0.1;
ebond = 0.0;
ev_init(eflag, vflag);
@ -219,14 +220,34 @@ void BondOxdnaFene::compute(int eflag, int vflag)
Deltasq = Delta[type] * Delta[type];
rlogarg = 1.0 - rr0sq / Deltasq;
// if r -> Delta, then rlogarg < 0.0 which is an error
// issue a warning and reset rlogarg = epsilon
// if r > 2*Delta something serious is wrong, abort
// energy
if (eflag) {
ebond = -0.5 * k[type] * log(rlogarg);
}
if (rlogarg < 0.1) {
if (rlogarg < rlogarg_min) {
// if r-r0 -> Delta, then rlogarg < 0.0 which is an error
// issue warning and reset rlogarg = rlogarg_min to cap force to
// F_max = F(r_max) = F(r_min) = F(r_0 +/- Delta*sqrt(1-rlogarg_min))
error->warning(FLERR, "FENE bond too long: {} {} {} {}", update->ntimestep, atom->tag[a],
atom->tag[b], r);
rlogarg = 0.1;
rlogarg = rlogarg_min;
// energy of capped force potential if r > r_max and r < r_min
if (eflag) {
// if overstretched E(r) = E(r_max) + F_max * (r-r_max)
if (r > r0[type]) {
ebond = -0.5 * k[type] * log(rlogarg) +
k[type] * sqrt(1.0-rlogarg) / rlogarg / Delta[type] *
(r - r0[type] - Delta[type] * sqrt(1.0-rlogarg));
}
// if overcompressed E(r) = E(r_min) + F_max * (r_min - r)
if (r < r0[type]) {
ebond = -0.5 * k[type] * log(rlogarg) +
k[type] * sqrt(1.0-rlogarg) / rlogarg / Delta[type] *
(r0[type] - Delta[type] * sqrt(1.0-rlogarg) - r);
}
}
}
fbond = -k[type] * rr0 / rlogarg / Deltasq / r;
@ -234,10 +255,6 @@ void BondOxdnaFene::compute(int eflag, int vflag)
delf[1] = delr[1] * fbond;
delf[2] = delr[2] * fbond;
// energy
if (eflag) { ebond = -0.5 * k[type] * log(rlogarg); }
// apply force and torque to each of 2 atoms
if (newton_bond || a < nlocal) {