Skip to content

Commit

Permalink
Merge pull request #339 from msoroush/issue335
Browse files Browse the repository at this point in the history 
Fix to the issue #338 and issue #335
  • Loading branch information
@GregorySchwing
GregorySchwing authored Apr 12, 2021
2 parents 385ac79 + a90923d commit f5bdcb3
Show file tree
Hide file tree
Showing 10 changed files with 49 additions and 40 deletions.
26 changes: 7 additions & 19 deletions src/CalculateEnergy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -710,9 +710,7 @@ 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)) {
if (currentAxes.InRcut(distSq, trialPos, t, trialMol.GetCoords(), *partner, box)) {
inter[t] += forcefield.particles->CalcEn(distSq,
kind.AtomKind(partIndex),
kind.AtomKind(*partner), 1.0);
Expand Down Expand Up @@ -1066,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
Expand Down Expand Up @@ -1102,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) {
Expand Down Expand Up @@ -1136,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]),
Expand Down Expand Up @@ -1173,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));
Expand Down Expand Up @@ -1207,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]),
Expand Down Expand Up @@ -1244,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));
Expand All @@ -1269,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;

Expand Down Expand Up @@ -1301,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;

Expand Down
6 changes: 6 additions & 0 deletions src/FFExp6.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down
6 changes: 6 additions & 0 deletions src/FFParticle.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down
6 changes: 6 additions & 0 deletions src/FFShift.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down
6 changes: 6 additions & 0 deletions src/FFSwitch.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down
6 changes: 6 additions & 0 deletions src/FFSwitchMartini.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
9 changes: 3 additions & 6 deletions src/cbmc/DCCrankShaftAng.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(),
*partner, box)) {
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);
Expand Down Expand Up @@ -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;
if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(),
*partner, box)) {
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));
Expand Down Expand Up @@ -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;
if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(),
*partner, box)) {
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));
Expand Down
9 changes: 3 additions & 6 deletions src/cbmc/DCCrankShaftDih.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(),
*partner, box)) {
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);
Expand Down Expand Up @@ -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;
if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(),
*partner, box)) {
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));
Expand Down Expand Up @@ -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;
if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(),
*partner, box)) {
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));
Expand Down
6 changes: 3 additions & 3 deletions src/cbmc/DCHedronCycle.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down
9 changes: 3 additions & 6 deletions src/cbmc/DCRotateOnAtom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(),
*partner, box)) {
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);
Expand Down Expand Up @@ -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;
if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(),
*partner, box)) {
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));
Expand Down Expand Up @@ -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;
if(data->axes.InRcut(distSq, trialPos, t, mol.GetCoords(),
*partner, box)) {
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));
Expand Down

0 comments on commit f5bdcb3

Please sign in to comment.