Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Cleanup: Force field #7406

Merged
merged 21 commits into from
May 10, 2024
Merged

Cleanup: Force field #7406

Show file tree
Hide file tree
Changes from 14 commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
2ed6b55
fix unsigned int to int comparison
AnnaBruenisholz Apr 17, 2024
fd8a0f3
revert previous mistake
AnnaBruenisholz Apr 17, 2024
2f6834c
declaration & init together, sinthetaSq in [0, 1]
AnnaBruenisholz Apr 22, 2024
6e56c25
using std::swap
AnnaBruenisholz Apr 22, 2024
4f3cd8f
use that sinThetaSq in [0,1]
AnnaBruenisholz Apr 22, 2024
446845b
declare & init at same time
AnnaBruenisholz Apr 22, 2024
df28453
use knowledge about target range
AnnaBruenisholz Apr 22, 2024
20d0bf9
use std::clamp
AnnaBruenisholz Apr 22, 2024
be8c773
use std::max
AnnaBruenisholz Apr 22, 2024
2d21acc
numerically more stable trigonometrics
AnnaBruenisholz Apr 30, 2024
786a0c5
numerically more stable trigonometrics
AnnaBruenisholz Apr 30, 2024
3b9c2f7
numerically more stable trigonometrics
AnnaBruenisholz Apr 30, 2024
0864889
range based for-loop
AnnaBruenisholz Apr 30, 2024
8432eb9
actually do the assignement
AnnaBruenisholz Apr 30, 2024
a486b39
Update Code/ForceField/MMFF/Params.h
AnnaBruenisholz May 8, 2024
f56e8f0
implemented suggested changes
AnnaBruenisholz May 8, 2024
e9b778f
Revert "implemented suggested changes"
AnnaBruenisholz May 10, 2024
09a73d2
auto typing
AnnaBruenisholz May 10, 2024
fbf1176
remove old comment
AnnaBruenisholz May 10, 2024
a627d7a
revert to numerically more stable expression
AnnaBruenisholz May 10, 2024
bd70aae
now correctly formatted
AnnaBruenisholz May 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
3 changes: 1 addition & 2 deletions Code/ForceField/MMFF/DistanceConstraint.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,8 +92,7 @@ void DistanceConstraintContrib::getGrad(double *pos, double *grad) const {
double *end1Coords = &(pos[3 * d_end1Idx]);
double *end2Coords = &(pos[3 * d_end2Idx]);
for (unsigned int i = 0; i < 3; ++i) {
double dGrad;
dGrad =
double dGrad =
preFactor * (end1Coords[i] - end2Coords[i]) / std::max(dist, 1.0e-8);
grad[3 * d_end1Idx + i] += dGrad;
grad[3 * d_end2Idx + i] -= dGrad;
Expand Down
82 changes: 33 additions & 49 deletions Code/ForceField/MMFF/Params.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@
// written the code for std::map, so the two methods
// can be toggled defining RDKIT_MMFF_PARAMS_USE_STD_MAP

//#define RDKIT_MMFF_PARAMS_USE_STD_MAP 1
// #define RDKIT_MMFF_PARAMS_USE_STD_MAP 1

namespace ForceFields {
namespace MMFF {
Expand All @@ -43,13 +43,7 @@ const double MDYNE_A_TO_KCAL_MOL = 143.9325;
inline bool isDoubleZero(const double x) {
return ((x < 1.0e-10) && (x > -1.0e-10));
}
inline void clipToOne(double &x) {
if (x > 1.0) {
x = 1.0;
} else if (x < -1.0) {
x = -1.0;
}
}
inline void clipToOne(double &x) { x = std::clamp(x, -1.0, 1.0); }

//! class to store MMFF atom type equivalence levels
class RDKIT_FORCEFIELD_EXPORT MMFFDef {
Expand Down Expand Up @@ -178,9 +172,8 @@ class RDKIT_FORCEFIELD_EXPORT MMFFDefCollection {
*/
const MMFFDef *operator()(const unsigned int atomType) const {
#ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
std::map<const unsigned int, MMFFDef>::const_iterator res;
res = d_params.find(atomType);

std::map<const unsigned int, MMFFDef>::const_iterator res =
AnnaBruenisholz marked this conversation as resolved.
Show resolved Hide resolved
d_params.find(atomType);
return ((res != d_params.end()) ? &((*res).second) : NULL);
#else
return ((atomType && (atomType <= d_params.size()))
Expand All @@ -206,8 +199,8 @@ class RDKIT_FORCEFIELD_EXPORT MMFFPropCollection {
*/
const MMFFProp *operator()(const unsigned int atomType) const {
#ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
std::map<const unsigned int, MMFFProp>::const_iterator res;
res = d_params.find(atomType);
std::map<const unsigned int, MMFFProp>::const_iterator res =
AnnaBruenisholz marked this conversation as resolved.
Show resolved Hide resolved
d_params.find(atomType);

return ((res != d_params.end()) ? &((*res).second) : NULL);
#else
Expand Down Expand Up @@ -239,9 +232,8 @@ class RDKIT_FORCEFIELD_EXPORT MMFFPBCICollection {
*/
const MMFFPBCI *operator()(const unsigned int atomType) const {
#ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
std::map<const unsigned int, MMFFPBCI>::const_iterator res;
res = d_params.find(atomType);

std::map<const unsigned int, MMFFPBCI>::const_iterator res =
d_params.find(atomType);
return ((res != d_params.end()) ? &((*res).second) : NULL);
#else
return ((atomType && (atomType <= d_params.size()))
Expand Down Expand Up @@ -279,9 +271,9 @@ class RDKIT_FORCEFIELD_EXPORT MMFFChgCollection {
}
#ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
std::map<const unsigned int,
std::map<const unsigned int, MMFFChg>>::const_iterator res1;
std::map<const unsigned int, MMFFChg>>::const_iterator res1 =
d_params[bondType].find(canIAtomType);
std::map<const unsigned int, MMFFChg>::const_iterator res2;
res1 = d_params[bondType].find(canIAtomType);
if (res1 != d_params[bondType].end()) {
res2 = ((*res1).second).find(canJAtomType);
if (res2 != ((*res1).second).end()) {
Expand All @@ -291,10 +283,8 @@ class RDKIT_FORCEFIELD_EXPORT MMFFChgCollection {
#else
std::pair<std::vector<std::uint8_t>::const_iterator,
std::vector<std::uint8_t>::const_iterator>
bounds;

bounds =
std::equal_range(d_iAtomType.begin(), d_iAtomType.end(), canIAtomType);
bounds = std::equal_range(d_iAtomType.begin(), d_iAtomType.end(),
AnnaBruenisholz marked this conversation as resolved.
Show resolved Hide resolved
canIAtomType);
if (bounds.first != bounds.second) {
bounds = std::equal_range(
d_jAtomType.begin() + (bounds.first - d_iAtomType.begin()),
Expand Down Expand Up @@ -350,11 +340,10 @@ class RDKIT_FORCEFIELD_EXPORT MMFFBondCollection {
std::map<const unsigned int,
std::map<const unsigned int,
std::map<const unsigned int, MMFFBond>>>::const_iterator
res1;
res1 = d_params.find(bondType);
std::map<const unsigned int,
std::map<const unsigned int, MMFFBond>>::const_iterator res2;
std::map<const unsigned int, MMFFBond>::const_iterator res3;
res1 = d_params.find(bondType);
if (res1 != d_params.end()) {
res2 = ((*res1).second).find(canAtomType);
if (res2 != ((*res1).second).end()) {
Expand All @@ -367,9 +356,8 @@ class RDKIT_FORCEFIELD_EXPORT MMFFBondCollection {
#else
std::pair<std::vector<std::uint8_t>::const_iterator,
std::vector<std::uint8_t>::const_iterator>
bounds;
bounds =
std::equal_range(d_iAtomType.begin(), d_iAtomType.end(), canAtomType);
bounds = std::equal_range(d_iAtomType.begin(), d_iAtomType.end(),
canAtomType);
if (bounds.first != bounds.second) {
bounds = std::equal_range(
d_jAtomType.begin() + (bounds.first - d_iAtomType.begin()),
Expand Down Expand Up @@ -420,9 +408,9 @@ class RDKIT_FORCEFIELD_EXPORT MMFFBndkCollection {
}
#ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
std::map<const unsigned int,
std::map<const unsigned int, MMFFBond>>::const_iterator res1;
std::map<const unsigned int, MMFFBond>>::const_iterator res1 =
d_params.find(canAtomicNum);
std::map<const unsigned int, MMFFBond>::const_iterator res2;
res1 = d_params.find(canAtomicNum);
if (res1 != d_params.end()) {
res2 = ((*res1).second).find(canNbrAtomicNum);
if (res2 != ((*res1).second).end()) {
Expand All @@ -432,9 +420,8 @@ class RDKIT_FORCEFIELD_EXPORT MMFFBndkCollection {
#else
std::pair<std::vector<std::uint8_t>::const_iterator,
std::vector<std::uint8_t>::const_iterator>
bounds;
bounds = std::equal_range(d_iAtomicNum.begin(), d_iAtomicNum.end(),
canAtomicNum);
bounds = std::equal_range(d_iAtomicNum.begin(), d_iAtomicNum.end(),
canAtomicNum);
if (bounds.first != bounds.second) {
bounds = std::equal_range(
d_jAtomicNum.begin() + (bounds.first - d_iAtomicNum.begin()),
Expand Down Expand Up @@ -477,9 +464,8 @@ class RDKIT_FORCEFIELD_EXPORT MMFFHerschbachLaurieCollection {
#ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
std::map<const unsigned int,
std::map<const unsigned int, MMFFHerschbachLaurie>>::const_iterator
res1;
res1 = d_params.find(canIRow);
std::map<const unsigned int, MMFFHerschbachLaurie>::const_iterator res2;
res1 = d_params.find(canIRow);
if (res1 != d_params.end()) {
res2 = ((*res1).second).find(canJRow);
if (res2 != ((*res1).second).end()) {
Expand All @@ -489,8 +475,7 @@ class RDKIT_FORCEFIELD_EXPORT MMFFHerschbachLaurieCollection {
#else
std::pair<std::vector<std::uint8_t>::const_iterator,
std::vector<std::uint8_t>::const_iterator>
bounds;
bounds = std::equal_range(d_iRow.begin(), d_iRow.end(), canIRow);
bounds = std::equal_range(d_iRow.begin(), d_iRow.end(), canIRow);
if (bounds.first != bounds.second) {
bounds = std::equal_range(
d_jRow.begin() + (bounds.first - d_iRow.begin()),
Expand Down Expand Up @@ -525,8 +510,8 @@ class RDKIT_FORCEFIELD_EXPORT MMFFCovRadPauEleCollection {
*/
const MMFFCovRadPauEle *operator()(const unsigned int atomicNum) const {
#ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
std::map<const unsigned int, MMFFCovRadPauEle>::const_iterator res;
res = d_params.find(atomicNum);
std::map<const unsigned int, MMFFCovRadPauEle>::const_iterator res =
d_params.find(atomicNum);

return ((res != d_params.end()) ? &((*res).second) : NULL);
#else
Expand Down Expand Up @@ -586,9 +571,9 @@ class RDKIT_FORCEFIELD_EXPORT MMFFAngleCollection {
unsigned int canIAtomType = (*mmffDef)(iAtomType)->eqLevel[iter];
unsigned int canKAtomType = (*mmffDef)(kAtomType)->eqLevel[iter];
if (canIAtomType > canKAtomType) {
unsigned int temp = canKAtomType;
canKAtomType = canIAtomType;
canIAtomType = temp;
unsigned int temp = canIAtomType;
canIAtomType = canKAtomType;
canKAtomType = temp;
Comment on lines +535 to +537
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This could have been a std::swap(canIAtomType, canKAtomType);

}
res1 = d_params.find(angleType);
if (res1 != d_params.end()) {
Expand Down Expand Up @@ -618,10 +603,11 @@ class RDKIT_FORCEFIELD_EXPORT MMFFAngleCollection {
unsigned int canIAtomType = (*mmffDef)(iAtomType)->eqLevel[iter];
unsigned int canKAtomType = (*mmffDef)(kAtomType)->eqLevel[iter];
if (canIAtomType > canKAtomType) {
unsigned int temp = canKAtomType;
canKAtomType = canIAtomType;
canIAtomType = temp;
unsigned int temp = canIAtomType;
canIAtomType = canKAtomType;
canKAtomType = temp;
Comment on lines +562 to +564
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

std::swap here too

}

bounds = std::equal_range(
d_iAtomType.begin() + (jBounds.first - d_jAtomType.begin()),
d_iAtomType.begin() + (jBounds.second - d_jAtomType.begin()),
Expand Down Expand Up @@ -692,15 +678,14 @@ class RDKIT_FORCEFIELD_EXPORT MMFFStbnCollection {
std::map<const unsigned int,
std::map<const unsigned int,
std::map<const unsigned int, MMFFStbn>>>>::
const_iterator res1;
const_iterator res1 = d_params.find(canStretchBendType);
std::map<const unsigned int,
std::map<const unsigned int,
std::map<const unsigned int, MMFFStbn>>>::const_iterator
res2;
std::map<const unsigned int,
std::map<const unsigned int, MMFFStbn>>::const_iterator res3;
std::map<const unsigned int, MMFFStbn>::const_iterator res4;
res1 = d_params.find(canStretchBendType);
if (res1 != d_params.end()) {
res2 = ((*res1).second).find(canIAtomType);
if (res2 != ((*res1).second).end()) {
Expand Down Expand Up @@ -1104,9 +1089,8 @@ class RDKIT_FORCEFIELD_EXPORT MMFFVdWCollection {
*/
const MMFFVdW *operator()(const unsigned int atomType) const {
#ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
std::map<const unsigned int, MMFFVdW>::const_iterator res;
res = d_params.find(atomType);

std::map<const unsigned int, MMFFVdW>::const_iterator res =
d_params.find(atomType);
return (res != d_params.end() ? &((*res).second) : NULL);
#else
std::pair<std::vector<std::uint8_t>::const_iterator,
Expand Down
3 changes: 2 additions & 1 deletion Code/ForceField/MMFF/PositionConstraint.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include <cmath>
#include <ForceField/ForceField.h>
#include <RDGeneral/Invariant.h>
#include <algorithm>

namespace ForceFields {
namespace MMFF {
Expand All @@ -39,7 +40,7 @@ double PositionConstraintContrib::getEnergy(double *pos) const {
RDGeom::Point3D p(pos[3 * d_atIdx], pos[3 * d_atIdx + 1],
pos[3 * d_atIdx + 2]);
double dist = (p - d_pos0).length();
double distTerm = (dist > d_maxDispl) ? dist - d_maxDispl : 0.0;
double distTerm = std::max(dist - d_maxDispl, 0.0);
double res = 0.5 * d_forceConstant * distTerm * distTerm;

return res;
Expand Down
3 changes: 1 addition & 2 deletions Code/ForceField/MMFF/StretchBend.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,8 +106,7 @@ void StretchBendContrib::getGrad(double *pos, double *grad) const {
double cosTheta = p12.dotProduct(p32);
clipToOne(cosTheta);
double sinThetaSq = 1.0 - cosTheta * cosTheta;
double sinTheta =
std::max(((sinThetaSq > 0.0) ? sqrt(sinThetaSq) : 0.0), 1.0e-8);
double sinTheta = std::max(sqrt(sinThetaSq), 1.0e-8);
AnnaBruenisholz marked this conversation as resolved.
Show resolved Hide resolved
double angleTerm = RAD2DEG * acos(cosTheta) - d_theta0;
double distTerm = RAD2DEG * (d_forceConstants.first * (dist1 - d_restLen1) +
d_forceConstants.second * (dist2 - d_restLen2));
Expand Down
30 changes: 7 additions & 23 deletions Code/ForceField/MMFF/testMMFFForceField.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@
#include <sstream>
#include <string>
#include "testMMFFForceField.h"
#include <utility>

#define FCON_TOLERANCE 0.05
#define ENERGY_TOLERANCE 0.025
Expand Down Expand Up @@ -181,46 +182,29 @@ bool sortTorsionInstanceVec(TorsionInstance *a, TorsionInstance *b) {

void fixBondStretchInstance(BondStretchInstance *bondStretchInstance) {
if (bondStretchInstance->iAtomType > bondStretchInstance->jAtomType) {
unsigned int temp;
temp = bondStretchInstance->iAtomType;
bondStretchInstance->iAtomType = bondStretchInstance->jAtomType;
bondStretchInstance->jAtomType = temp;
std::swap(bondStretchInstance->iAtomType, bondStretchInstance->jAtomType);
}
}

void fixAngleBendInstance(AngleBendInstance *angleBendInstance) {
if (angleBendInstance->iAtomType > angleBendInstance->kAtomType) {
unsigned int temp;
temp = angleBendInstance->iAtomType;
angleBendInstance->iAtomType = angleBendInstance->kAtomType;
angleBendInstance->kAtomType = temp;
std::swap(angleBendInstance->iAtomType, angleBendInstance->kAtomType);
}
}

void fixOopBendInstance(OopBendInstance *oopBendInstance) {
if (oopBendInstance->iAtomType > oopBendInstance->kAtomType) {
unsigned int temp;
temp = oopBendInstance->iAtomType;
oopBendInstance->iAtomType = oopBendInstance->kAtomType;
oopBendInstance->kAtomType = temp;
std::swap(oopBendInstance->iAtomType, oopBendInstance->kAtomType);
}
}

void fixTorsionInstance(TorsionInstance *torsionInstance) {
if (torsionInstance->jAtomType > torsionInstance->kAtomType) {
unsigned int temp;
temp = torsionInstance->jAtomType;
torsionInstance->jAtomType = torsionInstance->kAtomType;
torsionInstance->kAtomType = temp;
temp = torsionInstance->iAtomType;
torsionInstance->iAtomType = torsionInstance->lAtomType;
torsionInstance->lAtomType = temp;
std::swap(torsionInstance->jAtomType, torsionInstance->kAtomType);
std::swap(torsionInstance->iAtomType, torsionInstance->lAtomType);
} else if (torsionInstance->jAtomType == torsionInstance->kAtomType) {
unsigned int temp;
if (torsionInstance->iAtomType > torsionInstance->lAtomType) {
temp = torsionInstance->iAtomType;
torsionInstance->iAtomType = torsionInstance->lAtomType;
torsionInstance->lAtomType = temp;
std::swap(torsionInstance->iAtomType, torsionInstance->lAtomType);
}
}
}
Expand Down
3 changes: 1 addition & 2 deletions Code/ForceField/UFF/AngleBend.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,8 +162,7 @@ void AngleBendContrib::getGrad(double *pos, double *grad) const {
double cosTheta = r[0].dotProduct(r[1]);
clipToOne(cosTheta);
double sinThetaSq = 1.0 - cosTheta * cosTheta;
double sinTheta =
std::max(((sinThetaSq > 0.0) ? sqrt(sinThetaSq) : 0.0), 1.0e-8);
double sinTheta = std::max(sqrt(sinThetaSq), 1.0e-8);
AnnaBruenisholz marked this conversation as resolved.
Show resolved Hide resolved

// std::cerr << "GRAD: " << cosTheta << " (" << acos(cosTheta)<< "), ";
// std::cerr << sinTheta << " (" << asin(sinTheta)<< ")" << std::endl;
Expand Down
3 changes: 1 addition & 2 deletions Code/ForceField/UFF/DistanceConstraint.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,8 +92,7 @@ void DistanceConstraintContrib::getGrad(double *pos, double *grad) const {
double *end1Coords = &(pos[3 * d_end1Idx]);
double *end2Coords = &(pos[3 * d_end2Idx]);
for (unsigned int i = 0; i < 3; ++i) {
double dGrad;
dGrad =
double dGrad =
preFactor * (end1Coords[i] - end2Coords[i]) / std::max(dist, 1.0e-8);
grad[3 * d_end1Idx + i] += dGrad;
grad[3 * d_end2Idx + i] -= dGrad;
Expand Down
7 changes: 3 additions & 4 deletions Code/ForceField/UFF/Inversion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,12 +168,11 @@ void InversionContrib::getGrad(double *pos, double *grad) const {
double cosY = n.dotProduct(rJL);
clipToOne(cosY);
double sinYSq = 1.0 - cosY * cosY;
double sinY = std::max(((sinYSq > 0.0) ? sqrt(sinYSq) : 0.0), 1.0e-8);
double sinY = std::max(sqrt(sinYSq), 1.0e-8);
AnnaBruenisholz marked this conversation as resolved.
Show resolved Hide resolved
double cosTheta = rJI.dotProduct(rJK);
clipToOne(cosTheta);
double sinThetaSq = std::max(1.0 - cosTheta * cosTheta, 1.0e-8);
double sinTheta =
std::max(((sinThetaSq > 0.0) ? sqrt(sinThetaSq) : 0.0), 1.0e-8);
double sinThetaSq = 1.0 - cosTheta * cosTheta;
double sinTheta = std::max(sqrt(sinThetaSq), 1.0e-8);
// sin(2 * W) = 2 * sin(W) * cos(W) = 2 * cos(Y) * sin(Y)
double dE_dW = -d_forceConstant * (d_C1 * cosY - 4.0 * d_C2 * cosY * sinY);
RDGeom::Point3D t1 = rJL.crossProduct(rJK);
Expand Down
9 changes: 2 additions & 7 deletions Code/ForceField/UFF/Params.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
#include <string>
#include <cmath>
#include <map>
#include <algorithm>

#ifndef M_PI
#define M_PI 3.14159265358979323846
Expand All @@ -28,13 +29,7 @@ const double RAD2DEG = 180.0 / M_PI;
inline bool isDoubleZero(const double x) {
return ((x < 1.0e-10) && (x > -1.0e-10));
}
inline void clipToOne(double &x) {
if (x > 1.0) {
x = 1.0;
} else if (x < -1.0) {
x = -1.0;
}
}
inline void clipToOne(double &x) { x = std::clamp(x, -1.0, 1.0); }

//! class to store UFF parameters for bond stretching
class RDKIT_FORCEFIELD_EXPORT UFFBond {
Expand Down
3 changes: 2 additions & 1 deletion Code/ForceField/UFF/PositionConstraint.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
// of the RDKit source tree.
//
#include "PositionConstraint.h"
#include <algorithm>
#include <cmath>
#include <ForceField/ForceField.h>
#include <RDGeneral/Invariant.h>
Expand Down Expand Up @@ -39,7 +40,7 @@ double PositionConstraintContrib::getEnergy(double *pos) const {
RDGeom::Point3D p(pos[3 * d_atIdx], pos[3 * d_atIdx + 1],
pos[3 * d_atIdx + 2]);
double dist = (p - d_pos0).length();
double distTerm = (dist > d_maxDispl) ? dist - d_maxDispl : 0.0;
double distTerm = std::max(dist - d_maxDispl, 0.0);
double res = 0.5 * d_forceConstant * distTerm * distTerm;

return res;
Expand Down