Skip to content

Commit

Permalink
Add comparison for free energy: Jarzynski eq. vs. thermodynamic integ…
Browse files Browse the repository at this point in the history
…ration.
  • Loading branch information
kocimil1 committed Feb 6, 2025
1 parent 6c18fae commit 84ce683
Show file tree
Hide file tree
Showing 14 changed files with 189 additions and 23 deletions.
3 changes: 3 additions & 0 deletions cpp/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -167,6 +167,9 @@ add_subdirectory( ${MY_SRC_DIR}/libs )
# set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
add_subdirectory( tests )

# ====== Results
add_subdirectory( results )

# ====== OpenCL

if ( WITH_OPENCL )
Expand Down
6 changes: 3 additions & 3 deletions cpp/common/math/Forces.h
Original file line number Diff line number Diff line change
Expand Up @@ -339,14 +339,14 @@ inline double getLJQH( const Vec3d& dp, Vec3d& f, const Quat4d& REQH, const doub
double E,F;
// ---- Electrostatic
const double ir2_ = 1/( r2 + R2damp );
E = 0*COULOMB_CONST*REQH.z*sqrt( ir2_ );
F = 0*E*ir2_ ;
E = COULOMB_CONST*REQH.z*sqrt( ir2_ );
F = E*ir2_ ;
// --- LJ
const double ir2 = 1/r2;
const double u2 = REQH.x*REQH.x*ir2;
const double u6 = u2*u2*u2;
const double vdW = u6*REQH.y;
const double H = 0*u6*u6* ((REQH.w<0) ? REQH.w*REQH.y : 0.0); // H-bond correction
const double H = u6*u6* ((REQH.w<0) ? REQH.w*REQH.y : 0.0); // H-bond correction
E += (u6-2.)*vdW + H ;
F += ((u6-1.)*vdW + H )*ir2*12 ;
f.set_mul( dp, -F );
Expand Down
60 changes: 40 additions & 20 deletions cpp/common/molecular/MolWorld_sp3.h
Original file line number Diff line number Diff line change
Expand Up @@ -3233,24 +3233,27 @@ void thermodynamic_integration(int nbStep, int nMDSteps, double d_lamda, std::ve

}

void store_TI(std::string filename, std::vector<double> lamda, std::vector<double> TI, std::vector<double> sigmaTI, std::vector<double> Ref = std::vector<double>(0)){
void store_TI(std::string filename, std::vector<double> lamda, std::vector<double> TI, std::vector<double> sigmaTI = std::vector<double>(0), std::vector<double> Ref = std::vector<double>(0)){
std::ofstream outfile(filename);
if (!outfile.is_open()) {
std::cerr << "Error opening file: " << filename << std::endl;
return;
}
outfile << "lambda TI sigmaTI Reference" << std::endl;
outfile << "lambda TI";
if(sigmaTI.size()>0) outfile << " sigmaTI";
if(Ref.size()>0) outfile << "Reference";
outfile << std::endl;
for (int L = 0; L < lamda.size(); L++) {
outfile << lamda[L] << " " << TI[L] << " " << sigmaTI[L];
outfile << lamda[L] << " " << TI[L];
if(sigmaTI.size()>0) outfile << " " << sigmaTI[L];
if(Ref.size() > 0) outfile << " " << Ref[L];
outfile << std::endl;
}
outfile.close();
}

double mexican_hat_TI(double lamda1, double lamda2, int nbStep = 100, int nMDSteps = 100000, int nEQsteps = 10000, double tdamp = 100.0, double T = 300, double dt = 0.5){
double mexican_hat_TI(double lamda1, double lamda2, int nbStep = 100, int nMDSteps = 100000, int nEQsteps = 10000, double tdamp = 100.0, double T = 300, double dt = 0.5, double d = 0.1){
double gamma = 1/(dt*tdamp);
double d = 0.1;

std::vector<std::vector<double>> dE_dLamda(nbStep, std::vector<double>(nMDSteps));

Expand All @@ -3265,6 +3268,8 @@ double mexican_hat_TI(double lamda1, double lamda2, int nbStep = 100, int nMDSte
std::vector<double> TI(nbStep, 0.0);
std::vector<double> sigmaTI(nbStep, 0.0);
thermodynamic_integration(nbStep, nMDSteps, d_lamda, dE_dLamda.data(), TI.data(), sigmaTI.data());

store_TI("results/TI_plot.dat", lamda, TI);

double deltaF = TI[nbStep-1];
return deltaF;
Expand Down Expand Up @@ -3319,12 +3324,13 @@ void mexican_hat_MD_JE(double *x, double *v, double lamda, int nbMDsteps, double
//printf("temperature = %f\n", v2 / const_kB);
}

double mexican_hat_JE(double lamda1, double lamda2, int nbProdSteps = 100000, int nrealization = 1000, int nEQsteps = 10000, double tdamp = 100.0, double T = 300, double dt = 0.05){
double mexican_hat_JE(double lamda1, double lamda2, int nbProdSteps = 100000, int nrealization = 1000, int nEQsteps = 10000, double tdamp = 100.0, double T = 300, double dt = 0.05, int nSampleSteps = 0, double d = 0.1){

double gamma = 1/(dt*tdamp);
double d = 0.1;

int nSampleSteps = nEQsteps;
if(!nSampleSteps){
nSampleSteps = nEQsteps;
}

// the algorithm is not stable for higher x_eq (depends on nb MD steps)
std::vector<double> expave(nbProdSteps, 0.0);
Expand Down Expand Up @@ -3445,22 +3451,27 @@ int run_omp_three_atoms_problem(int niter_max, double dt, double Fconv = 1e-6, d
return itr;
}

double three_atoms_problem_JE(double lamda1, double lamda2, int nbProdSteps = 100000, int nrealization = 1000, int nEQsteps = 10000, double tdamp = 100.0, double T = 300, double dt = 0.05){
double three_atoms_problem_JE(double lamda1, double lamda2, int nbProdSteps = 100000, int nrealization = 1000, int nEQsteps = 10000, double tdamp = 100.0, double T = 300, double dt = 0.05, int nSampleSteps = 0){
ffl.apos[0].set(-1.0,0.0,0.0);
ffl.apos[1].set( 1.0,0.0,0.0);
ffl.apos[2].set( 1.0,3.0,0.0);
for(int i=0; i<ffl.natoms; i++){
ffl.REQs[i].x = 1.4430;
ffl.REQs[i].y = 0.5; // this number corresponds to write 0.1 in AtomTypes Hydrogen EvdW
ffl.REQs[i].y = 0.316228; // this number corresponds to write 0.1 in AtomTypes Hydrogen EvdW
ffl.REQs[i].z = 0.000000;
ffl.REQs[i].w = 0.000000;
}

int nSampleSteps = nEQsteps;

if(!nSampleSteps){
nSampleSteps = nEQsteps;
}

bFreeEnergyCalc = true;
ffl.bTestThreeAtoms = true; // turn off atom 0<->1 interaction
double gamma = 1/(dt*tdamp);
go.gamma_damp = gamma;
go.T_target = T;

std::vector<double> expave(nbProdSteps, 0.0);
double E_new, E_old;
Expand Down Expand Up @@ -3528,7 +3539,7 @@ double three_atoms_problem_JE(double lamda1, double lamda2, int nbProdSteps = 10
for (int i = 0; i < nbProdSteps; ++i) {
lamda[i] = lamda1 + i * d_lamda;
}
store_JE("results/JE_plot.dat", lamda, F_JE);
store_JE("results/JE_plot_3AP.dat", lamda, F_JE);

return F_JE[nbProdSteps-1];
}
Expand All @@ -3550,6 +3561,8 @@ double three_atoms_problem_TI(double lamda1, double lamda2, int nbStep = 100, in
bool calculate_temperature = true;
bMoving = true;
double gamma = 1/(dt*tdamp);
go.gamma_damp = gamma;
go.T_target = T;

std::vector<std::vector<double>> dE_dLamda(nbStep, std::vector<double>(nMDSteps));

Expand Down Expand Up @@ -3587,7 +3600,7 @@ double three_atoms_problem_TI(double lamda1, double lamda2, int nbStep = 100, in
std::vector<double> sigmaTI(nbStep, 0.0);
thermodynamic_integration(nbStep, nMDSteps, d_lamda, dE_dLamda.data(), TI.data(), sigmaTI.data());

//store_TI("results/TI_plot.dat", lamda, TI, sigmaTI);
store_TI("results/TI_plot_3AP.dat", lamda, TI, sigmaTI);

if (verbosity>3){
if (calculate_temperature)
Expand Down Expand Up @@ -3689,22 +3702,27 @@ int run_omp_entropic_spring(int niter_max, double dt, double Fconv = 1e-6, doubl
return itr;
}

double entropic_spring_JE(double lamda1, double lamda2, int n, int *dc, int nbProdSteps = 10000, int nrealization = 100, int nEQsteps = 10000, double tdamp = 100.0, double T = 300, double dt = 0.5)
double entropic_spring_JE(double lamda1, double lamda2, int n, int *dc, int nbProdSteps = 10000, int nrealization = 100, int nEQsteps = 10000, double tdamp = 100.0, double T = 300, double dt = 0.5, int nSampleSteps = 0)
{
ffl.doAngles = false;
ffl.doPiPiI = false;
ffl.doPiPiT = false;
ffl.doPiSigma = false;
ffl.doBonds = true;
ffl.bSubtractBondNonBond = true;
bConstrains = true;

bFreeEnergyCalc = true;
double gamma = 1 / (dt * tdamp);
go.gamma_damp = gamma;
go.T_target = T;

std::vector<double> expave(nbProdSteps, 0.0);
double E_new, E_old;

int nSampleSteps = nEQsteps;

if(!nSampleSteps){
nSampleSteps = nEQsteps;
}
// define colective variable
for (int i = 0; i < n; i++){
constrs.bonds[dc[i]].ls.set(lamda1 / n);
Expand Down Expand Up @@ -3779,22 +3797,25 @@ double entropic_spring_JE(double lamda1, double lamda2, int n, int *dc, int nbPr
double constant = 3 * const_kB * T / (lamda_bar * lamda_bar * ((double)ffl.natoms - 1));
Ref[L] = Ref[L - 1] + 0.5 * constant * (lamda[L] + lamda[L - 1]) * d_lamda;
}
store_JE("results/TI_plot.dat", lamda, F_JE, sigma_JE, Ref);
store_JE("results/JE_plot_ES.dat", lamda, F_JE, sigma_JE, Ref);

return F_JE[nbProdSteps-1];
}

double entropic_spring_TI(double lamda1, double lamda2, int n, int *dc, int nbStep = 100, int nMDsteps = 100000, int nEQsteps = 10000, double tdamp = 100.0, double T = 300, double dt = 0.5)
{
if(verbosity > 2)printf("MolWorld_sp3::entropic_spring_TI(%f, %f, %d, %d, %d, %d, %d, %f, %f, %f)\n", lamda1, lamda2, n, dc[0], nbStep, nMDsteps, nEQsteps, tdamp, T, dt);
ffl.doAngles = false;
ffl.doPiPiI = false;
ffl.doPiPiT = false;
ffl.doPiSigma= false;
ffl.doBonds = true;
ffl.bSubtractBondNonBond = true;
bConstrains = true;

bFreeEnergyCalc = true;
double gamma = 1 / (dt * tdamp);
go.gamma_damp = gamma;
go.T_target = T;

// define colective variable
for (int i = 0; i < n; i++)
Expand Down Expand Up @@ -3838,7 +3859,7 @@ double entropic_spring_JE(double lamda1, double lamda2, int n, int *dc, int nbPr
Ref[L] = Ref[L - 1] + 0.5 * constant * (lamda[L] + lamda[L - 1]) * d_lamda;
}

store_TI("results/TI_plot.dat", lamda, TI, sigmaTI, Ref);
store_TI("results/TI_plot_ES.dat", lamda, TI, sigmaTI, Ref);

return TI[nbStep - 1];
}
Expand All @@ -3852,7 +3873,6 @@ double entropic_spring_JE(double lamda1, double lamda2, int n, int *dc, int nbPr
double gamma = 1 / (dt * tdamp); // = go.gamma_damp
go.gamma_damp = gamma;
go.T_target = T;
int nbInits = 1000;

go.bExploring = true;
bConstrains = true;
Expand Down
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
File renamed without changes.
6 changes: 6 additions & 0 deletions cpp/results/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
# Set the output directory explicitly
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/results)

# Add the executable for the test
add_executable(free_energy free_energy.cpp)
target_link_libraries(free_energy PRIVATE MMFF_lib)
116 changes: 116 additions & 0 deletions cpp/results/free_energy.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,116 @@
#include "globals.h"
#include "testUtils.h"

#include "MolWorld_sp3.h"



char dir_cpp[1024];
void init(MolWorld_sp3& W, const char* xyz_name, const char* constr_name=nullptr){
W.smile_name = nullptr;
W.xyz_name = xyz_name;
W.surf_name = nullptr;
W.constr_name= constr_name;
W.bMMFF = true;
W.bEpairs = false;
W.gridStep = 0.01;
W.nPBC = Vec3i({1,1,1});
W.bUFF = false;
W.b141 = true;
W.bSimple = false;
W.bConj = true;
W.bCumulene = true;
W.bGridFF = false;



char path_ElementTypes[1024];
char path_AtomTypes[1024];
char path_BondTypes[1024];
char path_AngleTypes[1024];
char path_DihedralTypes[1024];
snprintf(path_ElementTypes, sizeof(path_ElementTypes), "%s/common_resources/ElementTypes.dat", dir_cpp);
snprintf(path_AtomTypes, sizeof(path_AtomTypes), "%s/common_resources/AtomTypes.dat", dir_cpp);
snprintf(path_BondTypes, sizeof(path_BondTypes), "%s/common_resources/BondTypes.dat", dir_cpp);
snprintf(path_AngleTypes, sizeof(path_AngleTypes), "%s/common_resources/AngleTypes.dat", dir_cpp);
snprintf(path_DihedralTypes, sizeof(path_DihedralTypes), "%s/common_resources/DihedralTypes.dat", dir_cpp);

verbosity = -1;
W.initParams( path_ElementTypes, path_AtomTypes, path_BondTypes, path_AngleTypes, path_DihedralTypes);

W.init();
}

void mexican_hat(double lamda1, double lamda2, int nbStep, int nMDsteps, int nEQsteps, double tdamp, double T, double dt, int nbProdSteps, int nrealization, int nSampleSteps){
char xyz_path [1024];
snprintf(xyz_path, sizeof(xyz_path), "%s/common_resources/xyz/H2O", dir_cpp);


MolWorld_sp3 W;
init(W, xyz_path);

// run mexican hat with TI
W.mexican_hat_TI(lamda1, lamda2, nbStep, nMDsteps, nEQsteps, tdamp, T, dt);

// run mexican hat with JE
W.mexican_hat_JE(lamda1, lamda2, nbProdSteps, nrealization, nEQsteps, tdamp, T, dt, nSampleSteps);
}

void three_atoms_problem(double lamda1, double lamda2, int nbStep, int nMDsteps, int nEQsteps, double tdamp, double T, double dt, int nbProdSteps, int nrealization, int nSampleSteps){
char xyz_path [1024];
snprintf(xyz_path, sizeof(xyz_path), "%s/common_resources/xyz/H2O", dir_cpp);

MolWorld_sp3 W;
init(W, xyz_path);

// run three atom problem with TI
W.three_atoms_problem_TI(lamda1, lamda2, nbStep, nMDsteps, nEQsteps, tdamp, T, dt);
// run three atom problem with JE
W.three_atoms_problem_JE(lamda1, lamda2, nbProdSteps, nrealization, nEQsteps, tdamp, T, dt, nSampleSteps);
}

void entropic_spring(double lamda1, double lamda2, int n, int *dc, int nbStep, int nMDsteps, int nEQsteps, double tdamp, double T, double dt, int nbProdSteps, int nrealization, int nSampleSteps, int natoms){
if(!dc){
dc = new int[n];
for(int i=0; i<n; i++){
dc[i] = i;
printf("%d\n", dc[i]);
}
}


char xyz_path [1024];
snprintf(xyz_path, sizeof(xyz_path), "%s/common_resources/entropic_spring_%d", dir_cpp, natoms);

char constr_path [1024];
snprintf(constr_path, sizeof(constr_path), "%s/common_resources/entropic_spring_%d.cons", dir_cpp, natoms);

MolWorld_sp3 W;
init(W, xyz_path, constr_path);

// run entropic spring with TI
W.entropic_spring_TI(lamda1, lamda2, n, dc, nbStep, nMDsteps, nEQsteps, tdamp, T, dt);
// run entropic spring with JE
W.entropic_spring_JE(lamda1, lamda2, n, dc, nbProdSteps, nrealization, nEQsteps, tdamp, T, dt, nSampleSteps);
}



int main(int argc, char **argv)
{
printf("generate data for toy models for free energy\n");

char* current_dir = getcwd(NULL, 0);
char* cpp_pos = strstr(current_dir, "/cpp");
if(cpp_pos) {*(cpp_pos+4) = '\0';}
else { "Error: cannot find /cpp in current directory path -> Cannot run program"; exit(1); }
snprintf(dir_cpp, sizeof(dir_cpp), "%s", current_dir);
free(current_dir);



//mexican_hat(std::stod(argv[1]), std::stod(argv[2]), std::stoi(argv[3]), std::stoi(argv[4]), std::stoi(argv[5]), std::stod(argv[6]), std::stod(argv[7]), std::stod(argv[8]), std::stoi(argv[9]), std::stoi(argv[10]), std::stoi(argv[11]));
//three_atoms_problem(std::stod(argv[1]), std::stod(argv[2]), std::stoi(argv[3]), std::stoi(argv[4]), std::stoi(argv[5]), std::stod(argv[6]), std::stod(argv[7]), std::stod(argv[8]), std::stoi(argv[9]), std::stoi(argv[10]), std::stoi(argv[11]));
entropic_spring(std::stod(argv[1]), std::stod(argv[2]), 1, nullptr, std::stoi(argv[3]), std::stoi(argv[4]), std::stoi(argv[5]), std::stod(argv[6]), std::stod(argv[7]), std::stod(argv[8]), std::stoi(argv[9]), std::stoi(argv[10]), std::stoi(argv[11]), std::stoi(argv[12]));
return 0;
}
21 changes: 21 additions & 0 deletions examples/tMMFF/results/JE_plot.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
# plot file TI_plot.dat with gnuplot. In the file, there are 3 columns "lamda", "FreeEnergy" and "Reference", which are the x and 2 y values. Add legend from the header of the file.
# /bin/python3 /home/kocimil1/Documents/testing_space/Free_energy/Three_particle_problem/main.py


gnuplot -e "
set terminal pngcairo size 1920,1080 enhanced font 'Verdana,10';
set output 'results/TI_plot.png';
set title 'TI plot';
set xlabel 'lamda';
set ylabel 'Energy [eV]';
set format x '%.1f';
set mxtics 5;
set xtics auto;
set ytics auto;
plot 'results/TI_plot.dat' using 1:2:3 with yerrorbars title 'Thermodynamic integration', \
'results/TI_plot.dat' using 1:4 with lines title 'Reference';
"
# show the plot
display results/TI_plot.png

0 comments on commit 84ce683

Please sign in to comment.