Skip to content

Commit

Permalink
for some reason no difference between dsp and dsf?
Browse files Browse the repository at this point in the history
  • Loading branch information
@GregorySchwing
GregorySchwing committed Feb 22, 2024
1 parent e6ae7b2 commit 3bee7a0
Show file tree
Hide file tree
Showing 4 changed files with 62 additions and 30 deletions.
27 changes: 27 additions & 0 deletions src/Forcefield.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,4 +158,31 @@ void Forcefield::SetWolfAlphaAndWolfFactors(double rcc, double wa, uint b){
void Forcefield::SetRCutCoulomb(double rcc, uint b){
rCutCoulomb[b]=rcc;
rCutCoulombSq[b]=rcc*rcc;
}


void Forcefield::SetWolfMethod(uint wolfmethod, uint potential){

if (wolfmethod == 0){
intramoleculardsf = false;
simpleself = false;
} else if (wolfmethod == 1) {
intramoleculardsf = false;
simpleself = true;
} else if (wolfmethod == 2) {
intramoleculardsf = true;
simpleself = false;
} else {
std::cerr << "ERROR: METHOD VALID VALUES [0,1,2]" << std::endl;
exit(1);
}

if (potential == 0){
dsf = false;
} else if (potential == 1){
dsf = true;
} else {
std::cerr << "ERROR: METHOD VALID VALUES [0,1]" << std::endl;
exit(1);
}
}
1 change: 1 addition & 0 deletions src/Forcefield.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ class Forcefield {

void SetWolfAlphaAndWolfFactors(double, double, uint);
void SetRCutCoulomb(double, uint);
void SetWolfMethod(uint, uint);

FFParticle *particles; //!< For LJ/Mie energy between unbonded atoms
// for LJ, shift and switch type
Expand Down
61 changes: 32 additions & 29 deletions src/WolfCalibrationOutput.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,11 +61,11 @@ void WolfCalibrationOutput::Init(pdb_setup::Atoms const& atoms,
//WriteGraceParFile();
for (uint b = 0; b < BOXES_WITH_U_NB; ++b) {
if(wolfAlphaRangeRead[b] && wolfCutoffCoulombRangeRead[b]){
//for (uint wolfKind = 0; wolfKind < 1; ++wolfKind){
for (uint wolfKind = 0; wolfKind < WOLF_TOTAL_KINDS; ++wolfKind){
for (uint coulKind = 0; coulKind < COUL_TOTAL_KINDS; ++coulKind){
sumRelativeErrorVec[b][coulKind].resize(alphaSize[b]*cutoffCoulombSize[b]);
sumRelativeErrorVec[b][wolfKind][coulKind].resize(alphaSize[b]*cutoffCoulombSize[b]);
}
//}
}
}
}
}
Expand Down Expand Up @@ -95,7 +95,7 @@ std::string WolfCalibrationOutput::getFileName(int b, int wolfKind, int coulKind
void WolfCalibrationOutput::WriteHeader()
{
for (uint b = 0; b < BOXES_WITH_U_NB; ++b) {
for (uint wolfKind = 0; wolfKind < 1; ++wolfKind){
for (uint wolfKind = 0; wolfKind < WOLF_TOTAL_KINDS; ++wolfKind){
for (uint coulKind = 0; coulKind < COUL_TOTAL_KINDS; ++coulKind){
outF.open((getFileName(b, wolfKind, coulKind, uniqueName)+".dat").c_str(), std::ofstream::out);
if (outF.is_open()) {
Expand Down Expand Up @@ -147,7 +147,7 @@ void WolfCalibrationOutput::WriteGraceParFile()
firstRow += "TITLE SIZE 2 \n";
firstRow += "LEGEND .8,.45\n";
firstRow += "with g0\n";
for (uint wolfKind = 0; wolfKind < 1; ++wolfKind){
for (uint wolfKind = 0; wolfKind < WOLF_TOTAL_KINDS; ++wolfKind){
for (uint coulKind = 0; coulKind < COUL_TOTAL_KINDS; ++coulKind){
std::string title = WOLF_KINDS[wolfKind] + " " + COUL_KINDS[coulKind];
firstRow += "\ts";
Expand Down Expand Up @@ -180,7 +180,7 @@ void WolfCalibrationOutput::WriteGraceParFileWRcut()
firstRow += "TITLE SIZE 2 \n";
firstRow += "LEGEND .8,.45\n";
firstRow += "with g0\n";
for (uint wolfKind = 0; wolfKind < 1; ++wolfKind){
for (uint wolfKind = 0; wolfKind < WOLF_TOTAL_KINDS; ++wolfKind){
for (uint coulKind = 0; coulKind < COUL_TOTAL_KINDS; ++coulKind){
int min_rcut_index=mapWK_CK_BOX_to_bestRCutIndex[std::make_tuple(wolfKind, coulKind, b)];
double best_rcut = wolfCutoffCoulombStart[b] + min_rcut_index*wolfCutoffCoulombDelta[b];
Expand Down Expand Up @@ -210,7 +210,7 @@ void WolfCalibrationOutput::DoOutput(const ulong step) {
// return;
CalculateGrid();
for (uint b = 0; b < BOXES_WITH_U_NB; ++b) {
for (uint wolfKind = 0; wolfKind < 1; ++wolfKind){
for (uint wolfKind = 0; wolfKind < WOLF_TOTAL_KINDS; ++wolfKind){
for (uint coulKind = 0; coulKind < COUL_TOTAL_KINDS; ++coulKind){
outF.open((getFileName(b, wolfKind, coulKind, uniqueName)+".dat").c_str(), std::ofstream::out);
if (outF.is_open()) {
Expand All @@ -227,7 +227,7 @@ void WolfCalibrationOutput::DoOutput(const ulong step) {
dataRow += std::to_string(alpha);
for (uint RCutIndex = 0; RCutIndex < cutoffCoulombSize[b]; ++RCutIndex) {
dataRow += ",";
double relativeError = 100.00*((sumRelativeErrorVec[b][coulKind][GetIndex(RCutIndex, alphaIndex, b)].mean()-ewaldAvg[b].mean())/ewaldAvg[b].mean());
double relativeError = 100.00*((sumRelativeErrorVec[b][wolfKind][coulKind][GetIndex(RCutIndex, alphaIndex, b)].mean()-ewaldAvg[b].mean())/ewaldAvg[b].mean());
dataRow += std::to_string(relativeError);
}
outF << dataRow << std::endl;
Expand All @@ -242,7 +242,7 @@ void WolfCalibrationOutput::DoOutput(const ulong step) {
outF.open((uniqueName + "_WOLF_CALIBRATION_BOX_" + std::to_string(b) + "_BEST_ALPHAS.csv").c_str(), std::ofstream::out);
if (outF.is_open()) {
std::string firstRow = "WolfKind,CoulKind,Rcut,Alpha,Error\n";
for (uint wolfKind = 0; wolfKind < 1; ++wolfKind){
for (uint wolfKind = 0; wolfKind < WOLF_TOTAL_KINDS; ++wolfKind){
for (uint coulKind = 0; coulKind < COUL_TOTAL_KINDS; ++coulKind){
int min_i = 0;
int min_a_index = 0;
Expand All @@ -253,14 +253,14 @@ void WolfCalibrationOutput::DoOutput(const ulong step) {
double alpha = wolfAlphaStart[b] + alphaIndex*wolfAlphaDelta[b];
for (uint RCutIndex = 0; RCutIndex < cutoffCoulombSize[b]; ++RCutIndex) {
double rCutCoulomb = wolfCutoffCoulombStart[b] + RCutIndex*wolfCutoffCoulombDelta[b];
double err = std::abs(sumRelativeErrorVec[b][coulKind][GetIndex(RCutIndex, alphaIndex, b)].mean()-ewaldAvg[b].mean());
double err = std::abs(sumRelativeErrorVec[b][wolfKind][coulKind][GetIndex(RCutIndex, alphaIndex, b)].mean()-ewaldAvg[b].mean());
double mean1 = ewaldAvg[b].mean();
double sd1 = ewaldAvg[b].sd();
double n = ewaldAvg[b].count();

double mean2 = sumRelativeErrorVec[b][coulKind][GetIndex(RCutIndex, alphaIndex, b)].mean();
double sd2 = sumRelativeErrorVec[b][coulKind][GetIndex(RCutIndex, alphaIndex, b)].sd();
double m = sumRelativeErrorVec[b][coulKind][GetIndex(RCutIndex, alphaIndex, b)].count();
double mean2 = sumRelativeErrorVec[b][wolfKind][coulKind][GetIndex(RCutIndex, alphaIndex, b)].mean();
double sd2 = sumRelativeErrorVec[b][wolfKind][coulKind][GetIndex(RCutIndex, alphaIndex, b)].sd();
double m = sumRelativeErrorVec[b][wolfKind][coulKind][GetIndex(RCutIndex, alphaIndex, b)].count();

double t_test;
if (n == 1)
Expand Down Expand Up @@ -302,10 +302,10 @@ void WolfCalibrationOutput::DoOutput(const ulong step) {
double a = wolfAlphaStart[b] + alphaIndex*wolfAlphaDelta[b];
std::string firstRow = "";
firstRow += std::to_string(a) + "\t";
for (uint wolfKind = 0; wolfKind < 1; ++wolfKind){
for (uint wolfKind = 0; wolfKind < WOLF_TOTAL_KINDS; ++wolfKind){
for (uint coulKind = 0; coulKind < COUL_TOTAL_KINDS; ++coulKind){
int min_rcut_index = mapWK_CK_BOX_to_bestRCutIndex[std::make_tuple(wolfKind, coulKind, b)];
double min_err = 100.00*((sumRelativeErrorVec[b][coulKind][GetIndex(min_rcut_index, alphaIndex, b)].mean()-ewaldAvg[b].mean())/ewaldAvg[b].mean());
double min_err = 100.00*((sumRelativeErrorVec[b][wolfKind][coulKind][GetIndex(min_rcut_index, alphaIndex, b)].mean()-ewaldAvg[b].mean())/ewaldAvg[b].mean());
firstRow += std::to_string(min_err) + "\t";
//firstRow += std::to_string(sumRelativeError[b][i]/numSamples) + "\t";
}
Expand Down Expand Up @@ -338,20 +338,23 @@ void WolfCalibrationOutput::CalculateGrid() {
for (uint b = 0; b < BOXES_WITH_U_NB; ++b) {
//printf("EwAtStep %lu %.*e\n", step, Digs, ewaldRef.boxEnergy[b].totalElect);
ewaldAvg[b].add_value(ewaldRef.boxEnergy[b].totalElect);
for (uint coulKind = 0; coulKind < COUL_TOTAL_KINDS; ++coulKind){
statValRef.forcefield.dsf = coulKind;
for (uint alphaIndex = 0; alphaIndex < alphaSize[b]; ++alphaIndex) {
double alpha = wolfAlphaStart[b] + alphaIndex*wolfAlphaDelta[b];
for (uint RCutIndex = 0; RCutIndex < cutoffCoulombSize[b]; ++RCutIndex) {
double rCutCoulomb = wolfCutoffCoulombStart[b] + RCutIndex*wolfCutoffCoulombDelta[b];
// Wolf class has references to these forcefield values
statValRef.forcefield.SetWolfAlphaAndWolfFactors(rCutCoulomb, alpha, b);
#ifdef GOMC_CUDA
statValRef.forcefield.particles->updateWolfEwald();
#endif
SystemPotential wolfTot = calcEn.SystemTotal();
wolfTot.Total();
sumRelativeErrorVec[b][coulKind][GetIndex(RCutIndex, alphaIndex, b)].add_value(wolfTot.boxEnergy[b].totalElect);
for (uint wolfKind = 0; wolfKind < WOLF_TOTAL_KINDS; ++wolfKind){
for (uint coulKind = 0; coulKind < COUL_TOTAL_KINDS; ++coulKind){
//statValRef.forcefield.dsf = coulKind;
statValRef.forcefield.SetWolfMethod(wolfKind,coulKind);
for (uint alphaIndex = 0; alphaIndex < alphaSize[b]; ++alphaIndex) {
double alpha = wolfAlphaStart[b] + alphaIndex*wolfAlphaDelta[b];
for (uint RCutIndex = 0; RCutIndex < cutoffCoulombSize[b]; ++RCutIndex) {
double rCutCoulomb = wolfCutoffCoulombStart[b] + RCutIndex*wolfCutoffCoulombDelta[b];
// Wolf class has references to these forcefield values
statValRef.forcefield.SetWolfAlphaAndWolfFactors(rCutCoulomb, alpha, b);
#ifdef GOMC_CUDA
statValRef.forcefield.particles->updateWolfEwald();
#endif
SystemPotential wolfTot = calcEn.SystemTotal();
wolfTot.Total();
sumRelativeErrorVec[b][wolfKind][coulKind][GetIndex(RCutIndex, alphaIndex, b)].add_value(wolfTot.boxEnergy[b].totalElect);
}
}
}
}
Expand Down
3 changes: 2 additions & 1 deletion src/WolfCalibrationOutput.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,8 @@ class WolfCalibrationOutput : public OutputableBase
bool wolfAlphaRangeRead[2];
bool wolfCutoffCoulombRangeRead[2];

std::vector<Welford<double>> sumRelativeErrorVec[BOX_TOTAL][COUL_TOTAL_KINDS];
//std::vector<Welford<double>> sumRelativeErrorVec[BOX_TOTAL][COUL_TOTAL_KINDS];
std::vector<Welford<double>> sumRelativeErrorVec[BOX_TOTAL][WOLF_TOTAL_KINDS][COUL_TOTAL_KINDS];
std::map<std::tuple<int, int, int>, int> mapWK_CK_BOX_to_bestRCutIndex;

Welford<double> ewaldAvg[BOX_TOTAL];
Expand Down

0 comments on commit 3bee7a0

Please sign in to comment.