From 807a31829e68e799d6682ab8b026f22a232bc07d Mon Sep 17 00:00:00 2001 From: Mohammad Soroush Date: Mon, 5 Apr 2021 17:21:51 -0400 Subject: [PATCH 1/2] Fix to the issue #338 and issue #335 Increase the tolerance to 10 degree. If sum Angles +/- 10 was close to 360, we keep the angles rigid for cyclic molecules. --- src/CalculateEnergy.cpp | 9 ++++----- src/cbmc/DCCrankShaftAng.cpp | 12 ++++++------ src/cbmc/DCCrankShaftDih.cpp | 12 ++++++------ src/cbmc/DCHedronCycle.cpp | 6 +++--- src/cbmc/DCRotateOnAtom.cpp | 12 ++++++------ 5 files changed, 25 insertions(+), 26 deletions(-) diff --git a/src/CalculateEnergy.cpp b/src/CalculateEnergy.cpp index 7f6ad5619..e6809b956 100644 --- a/src/CalculateEnergy.cpp +++ b/src/CalculateEnergy.cpp @@ -710,9 +710,8 @@ void CalculateEnergy::ParticleNonbonded(double* inter, if (trialMol.AtomExists(*partner)) { for (uint t = 0; t < trials; ++t) { double distSq; - - if (currentAxes.InRcut(distSq, trialPos, t, trialMol.GetCoords(), - *partner, box)) { + currentAxes.InRcut(distSq, trialPos, t, trialMol.GetCoords(), *partner, box); + if (forcefield.rCutSq > distSq) { inter[t] += forcefield.particles->CalcEn(distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner), 1.0); @@ -1269,7 +1268,7 @@ double CalculateEnergy::IntraEnergy_1_3(const double distSq, const uint atom1, { if(!forcefield.OneThree) return 0.0; - else if(forcefield.rCutSq < distSq) + else if(forcefield.rCutSq <= distSq) return 0.0; double eng = 0.0; @@ -1301,7 +1300,7 @@ double CalculateEnergy::IntraEnergy_1_4(const double distSq, const uint atom1, { if(!forcefield.OneFour) return 0.0; - else if(forcefield.rCutSq < distSq) + else if(forcefield.rCutSq <= distSq) return 0.0; double eng = 0.0; diff --git a/src/cbmc/DCCrankShaftAng.cpp b/src/cbmc/DCCrankShaftAng.cpp index 629c9df1b..718ffa18d 100644 --- a/src/cbmc/DCCrankShaftAng.cpp +++ b/src/cbmc/DCCrankShaftAng.cpp @@ -394,8 +394,8 @@ void DCCrankShaftAng::ParticleNonbonded1_N(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), - *partner, box)) { + data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); + if(data->ff.rCutSq > distSq) { nonbonded[t] += data->ff.particles->CalcEn(distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner), 1.0); @@ -431,8 +431,8 @@ void DCCrankShaftAng::ParticleNonbonded1_4(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), - *partner, box)) { + data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); + if(data->ff.rCutSq > distSq) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner)); @@ -468,8 +468,8 @@ void DCCrankShaftAng::ParticleNonbonded1_3(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), - *partner, box)) { + data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); + if(data->ff.rCutSq > distSq) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner)); diff --git a/src/cbmc/DCCrankShaftDih.cpp b/src/cbmc/DCCrankShaftDih.cpp index bdae2b0ca..51f9197c8 100644 --- a/src/cbmc/DCCrankShaftDih.cpp +++ b/src/cbmc/DCCrankShaftDih.cpp @@ -404,8 +404,8 @@ void DCCrankShaftDih::ParticleNonbonded1_N(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), - *partner, box)) { + data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); + if(data->ff.rCutSq > distSq) { nonbonded[t] += data->ff.particles->CalcEn(distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner), 1.0); @@ -441,8 +441,8 @@ void DCCrankShaftDih::ParticleNonbonded1_4(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), - *partner, box)) { + data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); + if(data->ff.rCutSq > distSq) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner)); @@ -478,8 +478,8 @@ void DCCrankShaftDih::ParticleNonbonded1_3(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), - *partner, box)) { + data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); + if(data->ff.rCutSq > distSq) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner)); diff --git a/src/cbmc/DCHedronCycle.cpp b/src/cbmc/DCHedronCycle.cpp index 058a7e54f..95b0f4e98 100644 --- a/src/cbmc/DCHedronCycle.cpp +++ b/src/cbmc/DCHedronCycle.cpp @@ -77,10 +77,10 @@ DCHedronCycle::DCHedronCycle(DCData* data, const mol_setup::MolKind& kind, for (uint a = 0; a < angles.size(); a++) { sumAngle += data->ff.angles->Angle(angles[a].kind); } - //If sum of angles = 2*pi = 6.283, it means they are in a plane + //If sum of angles (sumAngle +/- 10) ~ 2*pi = 6.283, it means they are in a plane //To avoid geometric conflict for flexible angle, we consider it fixed and - //let crankshaft sample it. 0.044 ~= 5 degree - bool angleInPlane = (std::abs(2.0 * M_PI - sumAngle) < 0.044); + //let crankshaft sample it. 0.1745 ~= 10 degree tolerance + bool angleInPlane = (std::abs(2.0 * M_PI - sumAngle) < 0.1745); bool constrainAngInRing = false; for (uint i = 0; i < nBonds; ++i) { diff --git a/src/cbmc/DCRotateOnAtom.cpp b/src/cbmc/DCRotateOnAtom.cpp index a48c093fc..94acdda05 100644 --- a/src/cbmc/DCRotateOnAtom.cpp +++ b/src/cbmc/DCRotateOnAtom.cpp @@ -422,8 +422,8 @@ void DCRotateOnAtom::ParticleNonbonded1_N(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), - *partner, box)) { + data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); + if(data->ff.rCutSq > distSq) { nonbonded[t] += data->ff.particles->CalcEn(distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner), 1.0); @@ -459,8 +459,8 @@ void DCRotateOnAtom::ParticleNonbonded1_4(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), - *partner, box)) { + data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); + if(data->ff.rCutSq > distSq) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner)); @@ -496,8 +496,8 @@ void DCRotateOnAtom::ParticleNonbonded1_3(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), - *partner, box)) { + data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); + if(data->ff.rCutSq > distSq) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner)); From a90923dbb92f5f6f9c6cfcccf31c2941faa62831 Mon Sep 17 00:00:00 2001 From: Mohammad Soroush Date: Mon, 12 Apr 2021 12:02:20 -0400 Subject: [PATCH 2/2] Implement general method to check distance for intra-nonbonded energies (elec and VdW) to avoid mistake in future implementations. All distance is now will be checked in force field class level. --- src/CalculateEnergy.cpp | 25 +++++++------------------ src/FFExp6.h | 6 ++++++ src/FFParticle.cpp | 6 ++++++ src/FFShift.h | 6 ++++++ src/FFSwitch.h | 6 ++++++ src/FFSwitchMartini.h | 6 ++++++ src/cbmc/DCCrankShaftAng.cpp | 9 +++------ src/cbmc/DCCrankShaftDih.cpp | 9 +++------ src/cbmc/DCRotateOnAtom.cpp | 9 +++------ 9 files changed, 46 insertions(+), 36 deletions(-) diff --git a/src/CalculateEnergy.cpp b/src/CalculateEnergy.cpp index e6809b956..4e320b649 100644 --- a/src/CalculateEnergy.cpp +++ b/src/CalculateEnergy.cpp @@ -710,8 +710,7 @@ void CalculateEnergy::ParticleNonbonded(double* inter, if (trialMol.AtomExists(*partner)) { for (uint t = 0; t < trials; ++t) { double distSq; - currentAxes.InRcut(distSq, trialPos, t, trialMol.GetCoords(), *partner, box); - if (forcefield.rCutSq > distSq) { + if (currentAxes.InRcut(distSq, trialPos, t, trialMol.GetCoords(), *partner, box)) { inter[t] += forcefield.particles->CalcEn(distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner), 1.0); @@ -1065,8 +1064,7 @@ void CalculateEnergy::MolNonbond(double & energy, for (uint i = 0; i < molKind.nonBonded.count; ++i) { uint p1 = mols.start[molIndex] + molKind.nonBonded.part1[i]; uint p2 = mols.start[molIndex] + molKind.nonBonded.part2[i]; - currentAxes.InRcut(distSq, currentCoords, p1, p2, box); - if (forcefield.rCutSq > distSq) { + if (currentAxes.InRcut(distSq, currentCoords, p1, p2, box)) { energy += forcefield.particles->CalcEn(distSq, molKind.AtomKind (molKind.nonBonded.part1[i]), molKind.AtomKind @@ -1101,8 +1099,7 @@ void CalculateEnergy::MolNonbond(double & energy, cbmc::TrialMol const &mol, uint p1 = molKind.nonBonded.part1[i]; uint p2 = molKind.nonBonded.part2[i]; if(mol.AtomExists(p1) && mol.AtomExists(p2)) { - currentAxes.InRcut(distSq, mol.GetCoords(), p1, p2, mol.GetBox()); - if (forcefield.rCutSq > distSq) { + if (currentAxes.InRcut(distSq, mol.GetCoords(), p1, p2, mol.GetBox())) { energy += forcefield.particles->CalcEn(distSq, molKind.AtomKind(p1), molKind.AtomKind(p2), 1.0); if (electrostatic) { @@ -1135,8 +1132,7 @@ void CalculateEnergy::MolNonbond_1_4(double & energy, for (uint i = 0; i < molKind.nonBonded_1_4.count; ++i) { uint p1 = mols.start[molIndex] + molKind.nonBonded_1_4.part1[i]; uint p2 = mols.start[molIndex] + molKind.nonBonded_1_4.part2[i]; - currentAxes.InRcut(distSq, currentCoords, p1, p2, box); - if (forcefield.rCutSq > distSq) { + if (currentAxes.InRcut(distSq, currentCoords, p1, p2, box)) { forcefield.particles->CalcAdd_1_4(energy, distSq, molKind.AtomKind (molKind.nonBonded_1_4.part1[i]), @@ -1172,8 +1168,7 @@ void CalculateEnergy::MolNonbond_1_4(double & energy, uint p1 = molKind.nonBonded_1_4.part1[i]; uint p2 = molKind.nonBonded_1_4.part2[i]; if(mol.AtomExists(p1) && mol.AtomExists(p2)) { - currentAxes.InRcut(distSq, mol.GetCoords(), p1, p2, mol.GetBox()); - if (forcefield.rCutSq > distSq) { + if (currentAxes.InRcut(distSq, mol.GetCoords(), p1, p2, mol.GetBox())) { forcefield.particles->CalcAdd_1_4(energy, distSq, molKind.AtomKind(p1), molKind.AtomKind(p2)); @@ -1206,8 +1201,7 @@ void CalculateEnergy::MolNonbond_1_3(double & energy, for (uint i = 0; i < molKind.nonBonded_1_3.count; ++i) { uint p1 = mols.start[molIndex] + molKind.nonBonded_1_3.part1[i]; uint p2 = mols.start[molIndex] + molKind.nonBonded_1_3.part2[i]; - currentAxes.InRcut(distSq, currentCoords, p1, p2, box); - if (forcefield.rCutSq > distSq) { + if (currentAxes.InRcut(distSq, currentCoords, p1, p2, box)) { forcefield.particles->CalcAdd_1_4(energy, distSq, molKind.AtomKind (molKind.nonBonded_1_3.part1[i]), @@ -1243,8 +1237,7 @@ void CalculateEnergy::MolNonbond_1_3(double & energy, uint p1 = molKind.nonBonded_1_3.part1[i]; uint p2 = molKind.nonBonded_1_3.part2[i]; if(mol.AtomExists(p1) && mol.AtomExists(p2)) { - currentAxes.InRcut(distSq, mol.GetCoords(), p1, p2, mol.GetBox()); - if (forcefield.rCutSq > distSq) { + if (currentAxes.InRcut(distSq, mol.GetCoords(), p1, p2, mol.GetBox())) { forcefield.particles->CalcAdd_1_4(energy, distSq, molKind.AtomKind(p1), molKind.AtomKind(p2)); @@ -1268,8 +1261,6 @@ double CalculateEnergy::IntraEnergy_1_3(const double distSq, const uint atom1, { if(!forcefield.OneThree) return 0.0; - else if(forcefield.rCutSq <= distSq) - return 0.0; double eng = 0.0; @@ -1300,8 +1291,6 @@ double CalculateEnergy::IntraEnergy_1_4(const double distSq, const uint atom1, { if(!forcefield.OneFour) return 0.0; - else if(forcefield.rCutSq <= distSq) - return 0.0; double eng = 0.0; diff --git a/src/FFExp6.h b/src/FFExp6.h index 86c20dec9..0a207068e 100644 --- a/src/FFExp6.h +++ b/src/FFExp6.h @@ -153,6 +153,9 @@ inline void FF_EXP6::Init(ff_setup::Particle const& mie, inline void FF_EXP6::CalcAdd_1_4(double& en, const double distSq, const uint kind1, const uint kind2) const { + if(forcefield.rCutSq < distSq) + return; + uint idx = FlatIndex(kind1, kind2); if(distSq < rMaxSq_1_4[idx]) { en += num::BIGNUM; @@ -174,6 +177,9 @@ inline void FF_EXP6::CalcCoulombAdd_1_4(double& en, const double distSq, const double qi_qj_Fact, const bool NB) const { + if(forcefield.rCutSq < distSq) + return; + double dist = sqrt(distSq); if(NB) en += qi_qj_Fact / dist; diff --git a/src/FFParticle.cpp b/src/FFParticle.cpp index 3d08149a0..667270b30 100644 --- a/src/FFParticle.cpp +++ b/src/FFParticle.cpp @@ -265,6 +265,9 @@ inline double FFParticle::GetRmax_1_4(const uint i, const uint j) const inline void FFParticle::CalcAdd_1_4(double& en, const double distSq, const uint kind1, const uint kind2) const { + if(forcefield.rCutSq < distSq) + return; + uint index = FlatIndex(kind1, kind2); double rRat2 = sigmaSq_1_4[index] / distSq; double rRat4 = rRat2 * rRat2; @@ -284,6 +287,9 @@ inline void FFParticle::CalcCoulombAdd_1_4(double& en, const double distSq, const double qi_qj_Fact, const bool NB) const { + if(forcefield.rCutSq < distSq) + return; + double dist = sqrt(distSq); if(NB) en += qi_qj_Fact / dist; diff --git a/src/FFShift.h b/src/FFShift.h index 8e9190c05..6fd1e9f8d 100644 --- a/src/FFShift.h +++ b/src/FFShift.h @@ -127,6 +127,9 @@ inline void FF_SHIFT::Init(ff_setup::Particle const& mie, inline void FF_SHIFT::CalcAdd_1_4(double& en, const double distSq, const uint kind1, const uint kind2) const { + if(forcefield.rCutSq < distSq) + return; + uint index = FlatIndex(kind1, kind2); double rRat2 = sigmaSq_1_4[index] / distSq; double rRat4 = rRat2 * rRat2; @@ -146,6 +149,9 @@ inline void FF_SHIFT::CalcCoulombAdd_1_4(double& en, const double distSq, const double qi_qj_Fact, const bool NB) const { + if(forcefield.rCutSq < distSq) + return; + double dist = sqrt(distSq); if(NB) en += qi_qj_Fact / dist; diff --git a/src/FFSwitch.h b/src/FFSwitch.h index ae1c7f951..c289b2e33 100644 --- a/src/FFSwitch.h +++ b/src/FFSwitch.h @@ -118,6 +118,9 @@ inline void FF_SWITCH::Init(ff_setup::Particle const& mie, inline void FF_SWITCH::CalcAdd_1_4(double& en, const double distSq, const uint kind1, const uint kind2) const { + if(forcefield.rCutSq < distSq) + return; + uint index = FlatIndex(kind1, kind2); double rCutSq_rijSq = forcefield.rCutSq - distSq; double rCutSq_rijSq_Sq = rCutSq_rijSq * rCutSq_rijSq; @@ -144,6 +147,9 @@ inline void FF_SWITCH::CalcCoulombAdd_1_4(double& en, const double distSq, const double qi_qj_Fact, const bool NB) const { + if(forcefield.rCutSq < distSq) + return; + double dist = sqrt(distSq); if(NB) en += qi_qj_Fact / dist; diff --git a/src/FFSwitchMartini.h b/src/FFSwitchMartini.h index d749e3d3e..72a483477 100644 --- a/src/FFSwitchMartini.h +++ b/src/FFSwitchMartini.h @@ -220,6 +220,9 @@ inline void FF_SWITCH_MARTINI::CalcAdd_1_4(double& en, const double distSq, const uint kind1, const uint kind2) const { + if(forcefield.rCutSq < distSq) + return; + uint index = FlatIndex(kind1, kind2); double r_2 = 1.0 / distSq; double r_4 = r_2 * r_2; @@ -253,6 +256,9 @@ inline void FF_SWITCH_MARTINI::CalcCoulombAdd_1_4(double& en, const double qi_qj_Fact, const bool NB) const { + if(forcefield.rCutSq < distSq) + return; + double dist = sqrt(distSq); if(NB) en += qi_qj_Fact / dist; diff --git a/src/cbmc/DCCrankShaftAng.cpp b/src/cbmc/DCCrankShaftAng.cpp index 718ffa18d..20fbf19c6 100644 --- a/src/cbmc/DCCrankShaftAng.cpp +++ b/src/cbmc/DCCrankShaftAng.cpp @@ -394,8 +394,7 @@ void DCCrankShaftAng::ParticleNonbonded1_N(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); - if(data->ff.rCutSq > distSq) { + if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box)) { nonbonded[t] += data->ff.particles->CalcEn(distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner), 1.0); @@ -431,8 +430,7 @@ void DCCrankShaftAng::ParticleNonbonded1_4(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); - if(data->ff.rCutSq > distSq) { + if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box)) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner)); @@ -468,8 +466,7 @@ void DCCrankShaftAng::ParticleNonbonded1_3(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); - if(data->ff.rCutSq > distSq) { + if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box)) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner)); diff --git a/src/cbmc/DCCrankShaftDih.cpp b/src/cbmc/DCCrankShaftDih.cpp index 51f9197c8..3810e3a99 100644 --- a/src/cbmc/DCCrankShaftDih.cpp +++ b/src/cbmc/DCCrankShaftDih.cpp @@ -404,8 +404,7 @@ void DCCrankShaftDih::ParticleNonbonded1_N(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); - if(data->ff.rCutSq > distSq) { + if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box)) { nonbonded[t] += data->ff.particles->CalcEn(distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner), 1.0); @@ -441,8 +440,7 @@ void DCCrankShaftDih::ParticleNonbonded1_4(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); - if(data->ff.rCutSq > distSq) { + if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box)) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner)); @@ -478,8 +476,7 @@ void DCCrankShaftDih::ParticleNonbonded1_3(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); - if(data->ff.rCutSq > distSq) { + if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box)) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner)); diff --git a/src/cbmc/DCRotateOnAtom.cpp b/src/cbmc/DCRotateOnAtom.cpp index 94acdda05..2f69e6483 100644 --- a/src/cbmc/DCRotateOnAtom.cpp +++ b/src/cbmc/DCRotateOnAtom.cpp @@ -422,8 +422,7 @@ void DCRotateOnAtom::ParticleNonbonded1_N(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); - if(data->ff.rCutSq > distSq) { + if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box)) { nonbonded[t] += data->ff.particles->CalcEn(distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner), 1.0); @@ -459,8 +458,7 @@ void DCRotateOnAtom::ParticleNonbonded1_4(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); - if(data->ff.rCutSq > distSq) { + if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box)) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner)); @@ -496,8 +494,7 @@ void DCRotateOnAtom::ParticleNonbonded1_3(cbmc::TrialMol const& mol, (std::find(atoms.begin(), atoms.end(), *partner) == atoms.end())) { for (uint t = 0; t < trials; ++t) { double distSq; - data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box); - if(data->ff.rCutSq > distSq) { + if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(), *partner, box)) { data->ff.particles->CalcAdd_1_4(nonbonded[t], distSq, kind.AtomKind(partIndex), kind.AtomKind(*partner));