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

Revert the sloppy intraBonded PR #474

Merged
merged 2 commits into from
Nov 29, 2022
Merged

Revert the sloppy intraBonded PR #474

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7a371ce
resolve conflicts
GregorySchwing Nov 7, 2022
2864c5b
ran clang formatting
GregorySchwing Nov 29, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion lib/Endian.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,6 @@
*/

#include <stdlib.h>
#include <stdint.h>

#define bswap_64(x) \
((((x)&0xff00000000000000ull) >> 56) | (((x)&0x00ff000000000000ull) >> 40) | \
Expand Down
58 changes: 15 additions & 43 deletions src/CalculateEnergy.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,9 +90,8 @@ SystemPotential CalculateEnergy::SystemTotal() {
// system intra
for (uint b = 0; b < BOX_TOTAL; ++b) {
GOMC_EVENT_START(1, GomcProfileEvent::EN_BOX_INTRA);
double bondEnergy[4] = {0};
double bondEn = 0.0, anglEn = 0.0, diheEn = 0.0, nonbondEn = 0.0,
correction = 0.0;
double bondEnergy[2] = {0};
double bondEn = 0.0, nonbondEn = 0.0, correction = 0.0;
MoleculeLookup::box_iterator thisMol = molLookup.BoxBegin(b);
MoleculeLookup::box_iterator end = molLookup.BoxEnd(b);
std::vector<uint> molID;
Expand All @@ -104,23 +103,18 @@ SystemPotential CalculateEnergy::SystemTotal() {

#ifdef _OPENMP
#pragma omp parallel for default(none) private(bondEnergy) shared(b, molID) \
reduction(+:bondEn, anglEn, diheEn, nonbondEn, correction)
reduction(+:bondEn, nonbondEn, correction)
#endif
for (int i = 0; i < (int)molID.size(); i++) {
// calculate nonbonded energy
MoleculeIntra(molID[i], b, bondEnergy);
bondEn += bondEnergy[0];
anglEn += bondEnergy[1];
diheEn += bondEnergy[2];
nonbondEn += bondEnergy[3];
nonbondEn += bondEnergy[1];
// calculate correction term of electrostatic interaction
correction += calcEwald->MolCorrection(molID[i], b);
}

pot.boxEnergy[b].intraBond = bondEn + anglEn + diheEn;
pot.boxEnergy[b].bond = bondEn;
pot.boxEnergy[b].angle = anglEn;
pot.boxEnergy[b].dihedral = diheEn;
pot.boxEnergy[b].intraBond = bondEn;
pot.boxEnergy[b].intraNonbond = nonbondEn;
// calculate self term of electrostatic interaction
pot.boxEnergy[b].self = calcEwald->BoxSelf(b);
Expand Down Expand Up @@ -869,26 +863,26 @@ Intermolecular CalculateEnergy::MoleculeTailVirChange(const uint box,
void CalculateEnergy::MoleculeIntra(const uint molIndex, const uint box,
double *bondEn) const {
GOMC_EVENT_START(1, GomcProfileEvent::EN_MOL_INTRA);
bondEn[0] = 0.0, bondEn[1] = 0.0, bondEn[2] = 0.0, bondEn[3] = 0.0;
bondEn[0] = 0.0, bondEn[1] = 0.0;

MoleculeKind &molKind = mols.kinds[mols.kIndex[molIndex]];
// *2 because we'll be storing inverse bond vectors
XYZArray bondVec(molKind.bondList.count * 2);

BondVectors(bondVec, molKind, molIndex, box);
MolBond(bondEn[0], molKind, bondVec, molIndex, box);
MolAngle(bondEn[1], molKind, bondVec, box);
MolDihedral(bondEn[2], molKind, bondVec, box);
MolNonbond(bondEn[3], molKind, molIndex, box);
MolNonbond_1_4(bondEn[3], molKind, molIndex, box);
MolNonbond_1_3(bondEn[3], molKind, molIndex, box);
MolAngle(bondEn[0], molKind, bondVec, box);
MolDihedral(bondEn[0], molKind, bondVec, box);
MolNonbond(bondEn[1], molKind, molIndex, box);
MolNonbond_1_4(bondEn[1], molKind, molIndex, box);
MolNonbond_1_3(bondEn[1], molKind, molIndex, box);
GOMC_EVENT_STOP(1, GomcProfileEvent::EN_MOL_INTRA);
}

// used in molecule exchange for calculating bonded and intraNonbonded energy
Energy CalculateEnergy::MoleculeIntra(cbmc::TrialMol const &mol) const {
GOMC_EVENT_START(1, GomcProfileEvent::EN_MOL_INTRA);
double bondEn = 0.0, anglEn = 0.0, diheEn = 0.0, intraNonbondEn = 0.0;
double bondEn = 0.0, intraNonbondEn = 0.0;
// *2 because we'll be storing inverse bond vectors
const MoleculeKind &molKind = mol.GetKind();
uint count = molKind.bondList.count;
Expand All @@ -897,35 +891,13 @@ Energy CalculateEnergy::MoleculeIntra(cbmc::TrialMol const &mol) const {

BondVectors(bondVec, mol, bondExist, molKind);
MolBond(bondEn, mol, bondVec, bondExist, molKind);
MolAngle(anglEn, mol, bondVec, bondExist, molKind);
MolDihedral(diheEn, mol, bondVec, bondExist, molKind);
MolAngle(bondEn, mol, bondVec, bondExist, molKind);
MolDihedral(bondEn, mol, bondVec, bondExist, molKind);
MolNonbond(intraNonbondEn, mol, molKind);
MolNonbond_1_4(intraNonbondEn, mol, molKind);
MolNonbond_1_3(intraNonbondEn, mol, molKind);
GOMC_EVENT_STOP(1, GomcProfileEvent::EN_MOL_INTRA);
return Energy(bondEn, anglEn, diheEn, bondEn + anglEn + diheEn,
intraNonbondEn, 0.0, 0.0, 0.0, 0.0, 0.0);
}

// used in molecule exchange for calculating bonded and intraNonbonded energy
Energy CalculateEnergy::UpdateBondAngleDihe(cbmc::TrialMol const &mol) const {
GOMC_EVENT_START(1, GomcProfileEvent::EN_MOL_INTRA);
double bondEn = 0.0, anglEn = 0.0, diheEn = 0.0, intraNonbondEn = 0.0;
// *2 because we'll be storing inverse bond vectors
const MoleculeKind &molKind = mol.GetKind();
uint count = molKind.bondList.count;
XYZArray bondVec(count * 2);
std::vector<bool> bondExist(count * 2, false);

BondVectors(bondVec, mol, bondExist, molKind);
MolBond(bondEn, mol, bondVec, bondExist, molKind);
MolAngle(anglEn, mol, bondVec, bondExist, molKind);
MolDihedral(diheEn, mol, bondVec, bondExist, molKind);
MolNonbond(intraNonbondEn, mol, molKind);
MolNonbond_1_4(intraNonbondEn, mol, molKind);
MolNonbond_1_3(intraNonbondEn, mol, molKind);
GOMC_EVENT_STOP(1, GomcProfileEvent::EN_MOL_INTRA);
return Energy(bondEn, anglEn, diheEn, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
return Energy(bondEn, intraNonbondEn, 0.0, 0.0, 0.0, 0.0, 0.0);
}

void CalculateEnergy::BondVectors(XYZArray &vecs, MoleculeKind const &molKind,
Expand Down
3 changes: 0 additions & 3 deletions src/CalculateEnergy.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -143,9 +143,6 @@ class CalculateEnergy {
// used in molecule exchange for calculating bonded and intraNonbonded energy
Energy MoleculeIntra(cbmc::TrialMol const &mol) const;

// used in molecule exchange for calculating bonded and intraNonbonded energy
Energy UpdateBondAngleDihe(cbmc::TrialMol const &mol) const;

//! Calculates Nonbonded 1_3 intramolecule energy of a full molecule
// for Martini forcefield
double IntraEnergy_1_3(const double distSq, const uint atom1,
Expand Down
7 changes: 4 additions & 3 deletions src/Clock.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,8 @@ along with this program, also can be found at
#endif

struct Clock {
Clock() : lastTime(0.0), stepsPerOut(0), prevStep(0), firstStep(0),
lastStep(0) {}
Clock()
: lastTime(0.0), stepsPerOut(0), prevStep(0), firstStep(0), lastStep(0) {}
void Init(const ulong steps, const ulong totSt, const ulong startStep) {
stepsPerOut = steps;
prevStep = startStep;
Expand Down Expand Up @@ -119,7 +119,8 @@ inline void Clock::CompletionTime(uint &day, uint &hr, uint &min) {
#if defined(__linux__) || defined(__APPLE__)
speed = (double)(prevStep - firstStep) / (lastTime - strt);
#elif (_WIN32) || (__CYGWIN__)
speed = (double)(prevStep - firstStep) / ((double)(lastTime - strt) / CLOCKS_PER_SEC);
speed = (double)(prevStep - firstStep) /
((double)(lastTime - strt) / CLOCKS_PER_SEC);
#endif
ulong rem = lastStep - prevStep;
ulong sec = rem / speed;
Expand Down
6 changes: 0 additions & 6 deletions src/ConsoleOutput.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -304,9 +304,6 @@ void ConsoleOutput::PrintEnergy(const uint box, Energy const &en,

printElement(en.total, elementWidth);
printElement(en.intraBond, elementWidth);
printElement(en.bond, elementWidth);
printElement(en.angle, elementWidth);
printElement(en.dihedral, elementWidth);
printElement(en.intraNonbond, elementWidth);
printElement(en.inter, elementWidth);
printElement(en.tailCorrection, elementWidth);
Expand All @@ -325,9 +322,6 @@ void ConsoleOutput::PrintEnergyTitle() {

printElement("TOTAL", elementWidth);
printElement("INTRA(B)", elementWidth);
printElement("BOND(B)", elementWidth);
printElement("ANGLE(B)", elementWidth);
printElement("DIHEDRAL(B)", elementWidth);
printElement("INTRA(NB)", elementWidth);
printElement("INTER(LJ)", elementWidth);
printElement("LRC", elementWidth);
Expand Down
142 changes: 65 additions & 77 deletions src/EnergyTypes.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,23 +77,14 @@ struct Intermolecular {
class Energy {
public:
Energy()
: bond(0.0), angle(0.0), dihedral(0.0), intraBond(0.0), intraNonbond(0.0),
inter(0.0), tailCorrection(0.0), total(0.0), real(0.0), recip(0.0),
self(0.0), correction(0.0), totalElect(0.0) {}
// For CBMC
Energy(double intraBond, double nonbond, double inter, double real,
double recip, double self, double correc)
: bond(0.0), angle(0.0), dihedral(0.0), intraBond(intraBond),
intraNonbond(nonbond), inter(inter), tailCorrection(0.0), total(0.0),
real(real), recip(recip), self(self), correction(correc),
totalElect(0.0) {}
Energy(double bond, double angle, double dihedral, double intraBond,
double nonbond, double inter, double real, double recip, double self,
double correc)
: bond(bond), angle(angle), dihedral(dihedral), intraBond(bond),
intraNonbond(nonbond), inter(inter), tailCorrection(0.0), total(0.0),
real(real), recip(recip), self(self), correction(correc),
: intraBond(0.0), intraNonbond(0.0), inter(0.0), tailCorrection(0.0),
total(0.0), real(0.0), recip(0.0), self(0.0), correction(0.0),
totalElect(0.0) {}
Energy(double bond, double nonbond, double inter, double real, double recip,
double self, double correc)
: intraBond(bond), intraNonbond(nonbond), inter(inter),
tailCorrection(0.0), total(0.0), real(real), recip(recip), self(self),
correction(correc), totalElect(0.0) {}

// VALUE SETTERS
double Total() {
Expand All @@ -108,9 +99,6 @@ class Energy {
}

void Zero() {
bond = 0.0;
angle = 0.0;
dihedral = 0.0;
intraBond = 0.0;
intraNonbond = 0.0;
inter = 0.0;
Expand Down Expand Up @@ -138,15 +126,12 @@ class Energy {

// private:
// MEMBERS
double bond, angle, dihedral, intraBond, intraNonbond, inter, tailCorrection,
total, real, recip, self, correction, totalElect;
double intraBond, intraNonbond, inter, tailCorrection, total, real, recip,
self, correction, totalElect;
};

inline Energy &Energy::operator-=(Energy const &rhs) {
inter -= rhs.inter;
bond -= rhs.bond;
angle -= rhs.angle;
dihedral -= rhs.dihedral;
intraBond -= rhs.intraBond;
intraNonbond -= rhs.intraNonbond;
tailCorrection -= rhs.tailCorrection;
Expand All @@ -162,9 +147,6 @@ inline Energy &Energy::operator-=(Energy const &rhs) {

inline Energy &Energy::operator+=(Energy const &rhs) {
inter += rhs.inter;
bond += rhs.bond;
angle += rhs.angle;
dihedral += rhs.dihedral;
intraBond += rhs.intraBond;
intraNonbond += rhs.intraNonbond;
tailCorrection += rhs.tailCorrection;
Expand All @@ -180,9 +162,6 @@ inline Energy &Energy::operator+=(Energy const &rhs) {

inline Energy &Energy::operator*=(double const &rhs) {
inter *= rhs;
bond *= rhs;
angle *= rhs;
dihedral *= rhs;
intraBond *= rhs;
intraNonbond *= rhs;
tailCorrection *= rhs;
Expand Down Expand Up @@ -436,26 +415,38 @@ inline bool SystemPotential::ComparePotentials(SystemPotential &other) {
std::cout << "difference: "
<< totalEnergy.intraBond - other.totalEnergy.intraBond
<< std::endl;
if (totalEnergy.bond != other.totalEnergy.bond) {
std::cout << "my bond: " << totalEnergy.bond
<< " other bond: " << other.totalEnergy.bond << std::endl;
std::cout << "difference: " << totalEnergy.bond - other.totalEnergy.bond
<< std::endl;
}
if (totalEnergy.angle != other.totalEnergy.angle) {
std::cout << "my angle: " << totalEnergy.angle
<< " other angle: " << other.totalEnergy.angle << std::endl;
std::cout << "difference: " << totalEnergy.angle - other.totalEnergy.angle
<< std::endl;
}
if (totalEnergy.dihedral != other.totalEnergy.dihedral) {
std::cout << "my dihedral: " << totalEnergy.dihedral
<< " other dihedral: " << other.totalEnergy.dihedral
<< std::endl;
std::cout << "difference: "
<< totalEnergy.dihedral - other.totalEnergy.dihedral
<< std::endl;
}
returnVal = false;
}
if (totalEnergy.intraNonbond != other.totalEnergy.intraNonbond) {
std::cout << "my intraNonbond: " << totalEnergy.intraNonbond
<< " other intraNonbond: " << other.totalEnergy.intraNonbond
<< std::endl;
std::cout << "difference: "
<< totalEnergy.intraNonbond - other.totalEnergy.intraNonbond
<< std::endl;
returnVal = false;
}
if (totalEnergy.tailCorrection != other.totalEnergy.tailCorrection) {
std::cout << "my LRC: " << totalEnergy.tailCorrection
<< " other LRC: " << other.totalEnergy.tailCorrection
<< std::endl;
std::cout << "difference: "
<< totalEnergy.tailCorrection - other.totalEnergy.tailCorrection
<< std::endl;
returnVal = false;
}
if (totalEnergy.real != other.totalEnergy.real) {
std::cout << "my real: " << totalEnergy.real
<< " other real: " << other.totalEnergy.real << std::endl;
std::cout << "difference: " << totalEnergy.real - other.totalEnergy.real
<< std::endl;
returnVal = false;
}
if (totalEnergy.recip != other.totalEnergy.recip) {
std::cout << "my recip: " << totalEnergy.recip
<< " other recip: " << other.totalEnergy.recip << std::endl;
std::cout << "difference: " << totalEnergy.recip - other.totalEnergy.recip
<< std::endl;
returnVal = false;
}
if (totalEnergy.intraNonbond != other.totalEnergy.intraNonbond) {
Expand Down Expand Up @@ -522,33 +513,30 @@ inline bool SystemPotential::ComparePotentials(SystemPotential &other) {
}

#ifndef NDEBUG
inline std::ostream &operator<<(std::ostream &out, Energy &en) {
// Save existing settings for the ostream
std::streamsize ss = out.precision();
std::ios_base::fmtflags ff = out.flags();

en.Total();
en.TotalElect();

out << std::setprecision(6) << std::fixed;
out << "\tTotal: " << en.total << " IntraB: " << en.intraBond
<< " Bond: " << en.bond << " Angle: " << en.angle
<< " Dihedral: " << en.dihedral << " IntraNB: " << en.intraNonbond
<< " Inter: " << en.inter
<< " Tail Correction: " << en.tailCorrection;

if (en.totalElect != 0.0) {
out << std::endl
<< "\tTotal Electric: " << en.totalElect << " Real: " << en.real
<< " Recip: " << en.recip << " Self: " << en.self
<< " Correction: " << en.correction;

// Restore ostream settings to prior value
out << std::setprecision(ss);
out.flags(ff);
}
return out;
}
inline std::ostream &operator<<(std::ostream &out, Energy &en) {
// Save existing settings for the ostream
std::streamsize ss = out.precision();
std::ios_base::fmtflags ff = out.flags();

en.Total();
en.TotalElect();

out << std::setprecision(6) << std::fixed;
out << "\tTotal: " << en.total << " IntraB: " << en.intraBond
<< " IntraNB: " << en.intraNonbond << " Inter: " << en.inter
<< " LRC: " << en.tailCorrection;
if (en.totalElect != 0.0) {
out << std::endl
<< "\tTotal Electric: " << en.totalElect << " Real: " << en.real
<< " Recip: " << en.recip << " Self: " << en.self
<< " Correction: " << en.correction;

// Restore ostream settings to prior value
out << std::setprecision(ss);
out.flags(ff);
}
return out;
}
#endif

#endif
Loading