diff --git a/cpp/common/molecular/MolWorld_sp3.h b/cpp/common/molecular/MolWorld_sp3.h index 6759db6a..3e34fcba 100644 --- a/cpp/common/molecular/MolWorld_sp3.h +++ b/cpp/common/molecular/MolWorld_sp3.h @@ -3435,7 +3435,21 @@ void runGlobalOptimization(int Findex, std::vector* a, std::vector* dE_dLamda){ + + + + + + + + + + + + +// Free Energy + +void MD_with_mexican_hat(int nbStep, int nbMDSteps, int nbMDStepsThermalize, double* lamda, double d, double T, double gamma, double dt, std::vector* dE_dLamda){ int atom = 0; double f, r, inv_c; @@ -3445,10 +3459,9 @@ void MD_with_mexican_hat(int nbStep, int nbMDSteps, int nbMDStepsThermalize, dou for(int i = 0; i < nbMDStepsThermalize; i++){ r = ffl.apos[atom].x; inv_c = 1/lamda[L]; - dE_dLamda[L][i] = - 4.0 * d * inv_c * r * r * inv_c * inv_c * ( r * r * inv_c * inv_c - 1); f = - 4.0 * d * inv_c * r * inv_c * ( r * inv_c * r * inv_c - 1); ffl.fapos[atom].set(f, 0, 0); - ffl.move_atom_Langevin( atom, dt_default, 10000.0, go.gamma_damp, go.T_target ); + ffl.move_atom_Langevin( atom, dt, 10000.0, gamma, T ); } for(int i = 0; i < nbMDSteps; i++){ @@ -3457,12 +3470,12 @@ void MD_with_mexican_hat(int nbStep, int nbMDSteps, int nbMDStepsThermalize, dou dE_dLamda[L][i] = - 4.0 * d * inv_c * r * r * inv_c * inv_c * ( r * r * inv_c * inv_c - 1); f = - 4.0 * d * inv_c * r * inv_c * ( r * inv_c * r * inv_c - 1); ffl.fapos[atom].set(f, 0, 0); - ffl.move_atom_Langevin( atom, dt_default, 10000.0, go.gamma_damp, go.T_target ); + ffl.move_atom_Langevin( atom, dt, 10000.0, gamma, T ); } } } -void JE_init(int nbInits, int nbMDSteps, int nbMDStepsThermalize, double lamda_init, double d, double* position){ +void JE_init(int nbInits, int nbMDSteps, int nbMDStepsThermalize, double lamda_init, double d, double T, double gamma, double dt, double* position){ int atom = 0; double f, r, inv_c; int n = nbMDSteps/nbInits; @@ -3475,7 +3488,7 @@ void JE_init(int nbInits, int nbMDSteps, int nbMDStepsThermalize, double lamda_i inv_c = 1/lamda_init; f = - 4.0 * d * inv_c * r * inv_c * ( r * inv_c * r * inv_c - 1); ffl.fapos[atom].set(f, 0, 0); - ffl.move_atom_Langevin( atom, dt_default, 10000.0, go.gamma_damp, go.T_target ); + ffl.move_atom_Langevin( atom, dt, 10000.0, gamma, T ); } for(int i = 0; i < nbMDSteps; i++){ @@ -3483,7 +3496,7 @@ void JE_init(int nbInits, int nbMDSteps, int nbMDStepsThermalize, double lamda_i inv_c = 1/lamda_init; f = - 4.0 * d * inv_c * r * inv_c * ( r * inv_c * r * inv_c - 1); ffl.fapos[atom].set(f, 0, 0); - ffl.move_atom_Langevin( atom, dt_default, 10000.0, go.gamma_damp, go.T_target ); + ffl.move_atom_Langevin( atom, dt, 10000.0, gamma, T ); n--; if(n == 0){ position[L] = r; @@ -3495,7 +3508,7 @@ void JE_init(int nbInits, int nbMDSteps, int nbMDStepsThermalize, double lamda_i } -void JE_propagation(double initial_position, double* lamda, int nbSteps, double d, double* W, int nbMD = 1000){ +void JE_propagation(double initial_position, double* lamda, int nbSteps, double d, double T, double gamma, double dt, double* W, int nbMD = 1000){ int atom = 0; double r, inv_c; int n = 1000; @@ -3507,15 +3520,13 @@ void JE_propagation(double initial_position, double* lamda, int nbSteps, double for(int L = 1; L < nbSteps; L++){ inv_c = 1/lamda[L]; - printf("r = %f\n", r); for (int t = 0; t < nbMD; t++) { f = -4.0 * d * inv_c * r * inv_c * (r * inv_c * r * inv_c - 1); ffl.fapos[atom].set(f, 0, 0); - ffl.move_atom_Langevin(atom, dt_default, 10000.0, go.gamma_damp, go.T_target); + ffl.move_atom_Langevin(atom, dt, 10000.0, gamma, T); r = ffl.apos[atom].x; } - if(L>1)E_prev = E; E = d * ( r * r * inv_c * inv_c - 1) * ( r * r * inv_c * inv_c - 1); @@ -3524,8 +3535,63 @@ void JE_propagation(double initial_position, double* lamda, int nbSteps, double } } +void thermodynamic_integration(int nbStep, int nMDSteps, double d_lamda, std::vector* dE_dLamda, double* TI, double* sigmaTI){ + std::vector Fs(nbStep, 0.0); + std::vector sigmaF(nbStep, 0.0); + for (int L = 0; L < nbStep; L++) + { + for (int i = 0; i < nMDSteps; i++) + { + Fs[L] += dE_dLamda[L][i]; + } + Fs[L] = Fs[L] / nMDSteps; + + double sigma = 0; + for(int i = 0; i < nMDSteps; i++){ + sigma += (dE_dLamda[L][i] - Fs[L]) * (dE_dLamda[L][i] - Fs[L]); + } + sigmaF[L] = std::sqrt(sigma / (nMDSteps - 1)); + + if (L > 0) + { + TI[L] = TI[L - 1] + 0.5 * (Fs[L] + Fs[L - 1]) * d_lamda; + + double variance_contribution = 0.25 * d_lamda * d_lamda * (sigmaF[L] * sigmaF[L] + sigmaF[L - 1] * sigmaF[L - 1]); + sigmaTI[L] = std::sqrt(sigmaTI[L - 1] * sigmaTI[L - 1] + variance_contribution); + } + else + { + TI[L] = 0.0; + sigmaTI[L] = 0.0; + } + printf("TI=%f, sigmaTI=%f\n", TI[L], sigmaTI[L]); + } + +} -double compute_Free_energy(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.05){ +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 gamma = 1/(dt*tdamp); + double d = 0.1; + + std::vector> dE_dLamda(nbStep, std::vector(nMDSteps)); + + std::vector lamda(nbStep); + double d_lamda = (lamda2 - lamda1) / (nbStep - 1); + for (int i = 0; i < nbStep; ++i) { + lamda[i] = lamda1 + i * d_lamda; + } + + MD_with_mexican_hat(nbStep, nMDSteps, nEQsteps, lamda.data(), d, T, gamma, dt, dE_dLamda.data()); + + std::vector TI(nbStep, 0.0); + std::vector sigmaTI(nbStep, 0.0); + thermodynamic_integration(nbStep, nMDSteps, d_lamda, dE_dLamda.data(), TI.data(), sigmaTI.data()); + + double deltaF = TI[nbStep-1]; + return deltaF; +} + +double compute_Free_energy(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){ // set global parameters double gamma = 1/(dt*tdamp); // = go.gamma_damp go.gamma_damp = gamma; @@ -3534,7 +3600,10 @@ double compute_Free_energy(double lamda1, double lamda2, int n, int *dc, int nbS go.bExploring = true; bConstrains = true; bFreeEnergyCalc = true; -/* + //return mexican_hat_TI(lamda1, lamda2, nbStep, nMDSteps, nEQsteps, tdamp, T, dt); + + +#ifdef mexican_hat_JE int nbInits = 1000; std::vector> W(nbInits, std::vector(nMDSteps)); std::vector initial_positions(nbInits, 0.0); @@ -3549,12 +3618,12 @@ int nbInits = 1000; } - JE_init(nbInits, nMDSteps, nbMDStepsThermalize, lamda1, d, initial_positions.data()); + JE_init(nbInits, nMDSteps, nbMDStepsThermalize, lamda1, d, T, gamma, dt, initial_positions.data()); std::vector boltzman_factor(nbStep, 0.0); std::vector sigma_boltzman_factor(nbStep, 0.0); for (int i=0; i F_JE(nbStep, 0.0); std::vector sigma_JE(nbStep, 0.0); @@ -3593,37 +3662,10 @@ std::vector Ref(nbStep, 0.0); file << lamda[i] << " " << (F_JE[i]) << " " << sigma_JE[i] << " " << Ref[i] << std::endl; } file.close(); -*/ - // Mexican hat potential (work) -/* double d = 0.1; - - std::vector> dE_dLamda(nbStep, std::vector(nMDSteps))+; - int nbMDStepsThermalize = 1000; - std::vector lamda(nbStep); - double d_lamda = (lamda2 - lamda1) / (nbStep - 1); - for (int i = 0; i < nbStep; ++i) { - lamda[i] = lamda1 + i * d_lamda; - } +#endif - MD_with_mexican_hat(nbStep, nMDSteps, nbMDStepsThermalize, lamda.data(), d, dE_dLamda.data()); - - std::vector Fs(nbStep, 0.0); - std::vector TI(nbStep, 0.0); - for(int L = 0; L < nbStep; L++){ - for(int i = 0; i < nMDSteps; i++){ - Fs[L] += dE_dLamda[L][i]; - } - Fs[L] = Fs[L]/nMDSteps; - if(L > 0) - TI[L] = TI[L-1] + Fs[L] * d_lamda; - else - TI[L] = Fs[L] * d_lamda; - printf("%f\n", TI[L]); - } - double deltaF = TI[nbStep-1]; - return deltaF; - */ + /* @@ -3769,8 +3811,8 @@ double constant = 3 * const_kB * T / (1.2 * 1.2 * ((double)ffl.natoms-1)); printf("lamda=%f, lamda_bar=%f, TI=%f, Ref=%f, T_measured=%f\n", lamda[L], lamda_bar, TI[L], Ref[L], T_measured); } Ref[0] = 0.0; - -/* +*/ +#ifdef three_atoms_problem_TI //////////////////////////// Three atom problem //////////////////////////////////// //std::vector> Energy(nbStep, std::vector(nMDSteps)); std::vector> dE_dLamda(nbStep, std::vector(nMDSteps)); @@ -3838,9 +3880,9 @@ ffl.apos[1].x += d_lamda*0.5; \ // printf("ffl.apos[0]=(%f, %f, %f), ffl.apos[1]=(%f, %f, %f), lamda=%f, lamda_bar=%f, TI=%f, T_measured=%f\n", ffl.apos[0].x, ffl.apos[0].y, ffl.apos[0].z, ffl.apos[1].x, ffl.apos[1].y, ffl.apos[1].z, lamda[L], lamda_bar, TI[L], T_measured); } -*/ -/* + + double deltaF = TI[nbStep - 1]; // save lamda and TI into dat file with headers "lamda" and "FreeEnergy" @@ -3852,7 +3894,8 @@ ffl.apos[1].x += d_lamda*0.5; \ } file.close(); - printf("deltaF=%f\n", deltaF);*/ + printf("deltaF=%f\n", deltaF); +#endif return 0; } }; diff --git a/cpp/common/molecular/constrains.h b/cpp/common/molecular/constrains.h index 46f55a9a..80898dee 100644 --- a/cpp/common/molecular/constrains.h +++ b/cpp/common/molecular/constrains.h @@ -99,7 +99,7 @@ struct DistConstr{ double l = d.norm(); double f,E; E = spring( l, ls, ks, flim, f ); - printf( "ks=(%g,%g) ls(%g,%g) l %g f %g E %g \n", ks.x,ks.y, ls.x,ls.y, l,f,E ); + //f = ks.x*l; E= 0.5*ks.x*l*l; d.mul(f/l); // fs[ias.b].sub(d); diff --git a/cpp/tests/test_free_energy.cpp b/cpp/tests/test_free_energy.cpp index 76edc830..7fe59913 100644 --- a/cpp/tests/test_free_energy.cpp +++ b/cpp/tests/test_free_energy.cpp @@ -1,41 +1,91 @@ #include +#include "globals.h" #include "testUtils.h" #include "MolWorld_sp3.h" -// def init( -// xyz_name =None, -// surf_name =None, -// smile_name=None, -// constr_name=None, -// sElementTypes = "data/ElementTypes.dat", -// sAtomTypes = "data/AtomTypes.dat", -// sBondTypes = "data/BondTypes.dat", -// sAngleTypes = "data/AngleTypes.dat", -// sDihedralTypes = "data/DihedralTypes.dat", -// bMMFF=True, -// bEpairs=False, -// nPBC=(1,1,0), -// gridStep=0.1, -// bUFF=False, -// b141=True, -// bSimple=False, -// bConj=True, -// bCumulene=True -// ): -// global glob_bMMFF -// glob_bMMFF = bMMFF -// nPBC=np.array(nPBC,dtype=np.int32) -// return lib.init( cstr(xyz_name), cstr(surf_name), cstr(smile_name), cstr(constr_name), bMMFF, bEpairs, bUFF, b141, bSimple, bConj, bCumulene, nPBC, gridStep, cstr(sElementTypes), cstr(sAtomTypes), cstr(sBondTypes), cstr(sAngleTypes), cstr(sDihedralTypes) ) -int main(int argc, char **argv) { - // ::testing::InitGoogleTest(&argc, argv); - MolWorld_sp3 W; + + +char dir_cpp[1024]; +void init(MolWorld_sp3& W){ + W.smile_name = nullptr; + W.xyz_name = "../common_resources/xyz/pyridine"; + W.surf_name = nullptr; + W.constr_name= nullptr; + 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); + + W.initParams( path_ElementTypes, path_AtomTypes, path_BondTypes, path_AngleTypes, path_DihedralTypes); + + + printf("INITTIALIZING\n"); W.init(); +} + +TEST(MolWorld_sp3, initParams) { + MolWorld_sp3 W; + 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); + + W.initParams( path_ElementTypes, path_AtomTypes, path_BondTypes, path_AngleTypes, path_DihedralTypes); + ASSERT_GT(W.params.atypes.size(), 0); +} + + + +TEST(MolWorld_sp3, mexican_hat_TI ){ + MolWorld_sp3 W; + init(W); + double initial = 0.5; + double final = 4.0; + int N = 1; + int colectiveVariable[N] = {0}; + double E = W.compute_Free_energy(initial, final, N, colectiveVariable, 100, 100000, 10000, 100, 100, 0.5); + printf("E = %f\n", E); + EXPECT_DOUBLE_EQ(E, 0.08241649205729196); +} + +int main(int argc, char **argv) { + ::testing::InitGoogleTest(&argc, argv); char* current_dir = getcwd(NULL, 0); - printf("%s", current_dir); - return 0;//RUN_ALL_TESTS(); + char* cpp_pos = strstr(current_dir, "/cpp"); + if(cpp_pos) {*(cpp_pos+4) = '\0';} + else { "Error: cannot find /cpp in current directory -> Cannot run tests"; exit(1); } + snprintf(dir_cpp, sizeof(dir_cpp), "%s", current_dir); + free(current_dir); + + return RUN_ALL_TESTS(); } diff --git a/examples/tMMFF/run.py b/examples/tMMFF/run.py index c0265a69..891e3d8e 100755 --- a/examples/tMMFF/run.py +++ b/examples/tMMFF/run.py @@ -29,12 +29,12 @@ def scanPlot( nscan = 1000, span=(0.0,8.0), dir=(0.0,0.0,1.0), p0=(0.0,0.0,0.0), nbStep = 100 nMDsteps = 100000 nEQsteps = 10000 -t_damp = 20 -T = 10 -dt = 0.05 +t_damp = 100 +T = 100 +dt = 0.5 mmff.setVerbosity(verbosity=3) -mmff.init( xyz_name="data/enthropic_spring_10", bMMFF=True ) +mmff.init( xyz_name="data/H2", bMMFF=True ) collectiveVariable = np.array([0], dtype=np.int32) E = mmff.compute_Free_energy(0.5, 4.0, collectiveVariable, nbStep=nbStep, nMDsteps=nMDsteps, nEQsteps=nEQsteps, t_damp=t_damp, T=T, dt=dt) print("E=", E) diff --git a/examples/tMMFF/run_sample.py b/examples/tMMFF/run_sample.py deleted file mode 100755 index 569a5f1b..00000000 --- a/examples/tMMFF/run_sample.py +++ /dev/null @@ -1,171 +0,0 @@ -import sys -import numpy as np -import os -import matplotlib.pyplot as plt - -sys.path.append("../../") -from pyBall import atomicUtils as au -from pyBall import MMFF as mmff - -def numDeriv( xs, Es): - dx = xs[1]-xs[0] - Fs = (Es[2:]-Es[:-2])/(2*dx) - return -Fs - -def cos_half( a ): - E = (1-np.cos(a*0.5))*0.5 - F = np.sin(a*0.5)*0.25 - return E,F - -#======== Body - -''' -# ----- Angle -#mmff.sample_evalAngleCos( xs, lmin=1, lmax=1, kmin=1, kmax=1, flim=1e+300, Es=None, Fs=None) -xs = np.linspace(-np.pi,np.pi,100) -#Es,Fs = mmff.sample_evalAngleCos( xs, ang0=np.pi*0.00 ) -#print("Fs",Fs) -#Es,Fs = mmff.sample_evalAngleCos( xs, ang0=np.pi*0.25 ) -#Es,Fs = mmff.sample_evalAngleCos( xs, ang0=np.pi*0.50 ) -#Es,Fs = mmff.sample_evalAngleCos( xs, ang0=np.pi*0.75 ) -#Es,Fs = mmff.sample_evalAngleCos( xs, ang0=np.pi*1.00 ) - -Es,Fs = mmff.sample_evalAngleCosHalf( xs, ang0=np.pi*0.00 ) -#Es,Fs = mmff.sample_evalAngleCosHalf( xs, ang0=np.pi*0.1 ) -#Es,Fs = mmff.sample_evalAngleCosHalf( xs, ang0=np.pi*0.25 ) -#Es,Fs = mmff.sample_evalAngleCosHalf( xs, ang0=np.pi*0.50 ) -#Es,Fs = mmff.sample_evalAngleCosHalf ( xs, ang0=np.pi*0.75 ) -#Es,Fs = mmff.sample_evalAngleCosHalf( xs, ang0=np.pi*0.9 ) -#Es,Fs = mmff.sample_evalAngleCosHalf( xs, ang0=np.pi*1.00 ) -Fnum = numDeriv(xs,Es) #ratios = Fs[1:-1]/Fnum ;print("ratios", ratios) -#Es,Fs = cos_half( xs + np.pi*0.5 ) -xs/=np.pi -plt.figure(); plt.plot(xs, Es, label="E"); plt.plot(xs, Fs, label="F_ana"); plt.plot(xs[1:-1], Fnum, label="F_num"); plt.grid(); plt.legend() -plt.show() -''' - - -''' -# ----- PiPiAlignment -xs = np.linspace(-np.pi,np.pi,1000) -Es,Fs = mmff.sample_evalPiAling( xs, ang0=np.pi*0.00 ) -Fnum = numDeriv(xs,Es) #ratios = Fs[1:-1]/Fnum ;print("ratios", ratios) -xs/=np.pi -plt.figure(); plt.plot(xs, Es, label="E"); plt.plot(xs, Fs, label="F_ana"); plt.plot(xs[1:-1], Fnum, label="F_num"); plt.grid(); plt.legend() -plt.show() -''' - -''' -# ----- Bond -xs = np.linspace(0.1,3.0,100) -Es,Fs = mmff.sample_evalBond( xs) -Fnum = numDeriv(xs,Es) #ratios = Fs[1:-1]/Fnum ;print("ratios", ratios) -plt.figure(); plt.plot(xs, Es, label="E"); plt.plot(xs, Fs, label="F_ana"); plt.plot(xs[1:-1], Fnum, label="F_num"); plt.grid(); plt.legend() -plt.show() -''' - -''' -# ----- evalAtom -xs = np.linspace(-3.0,3.0,1000) -mmff.init( xyz_name="data/polymer-2_new" ) -Es,Fs = mmff.sample_evalAtom( xs, ia=1) -Fnum = numDeriv(xs,Es) #ratios = Fs[1:-1]/Fnum ;print("ratios", ratios) -plt.figure(); plt.plot(xs, Es, label="E"); plt.plot(xs, Fs, label="F_ana"); plt.plot(xs[1:-1], Fnum, label="F_num"); plt.grid(); plt.legend() -plt.show() -''' - -''' -# ----- getLJQH -xs = np.linspace(0.1,3.0,1000) -Es,Fs = mmff.sample_getLJQH( xs) -Fnum = numDeriv(xs,Es) #ratios = Fs[1:-1]/Fnum ;print("ratios", ratios) -plt.figure(); plt.plot(xs, Es, label="E"); plt.plot(xs, Fs, label="F_ana"); plt.plot(xs[1:-1], Fnum, label="F_num"); plt.grid(); plt.legend() -plt.show() -''' - -''' -# ----- evalLJQs_ng4_PBC_atom_omp, evalLJQs_ng4_atom_omp,... -xs = np.linspace(-3.0,3.0,100) -mmff.init( xyz_name="data/pyridine" ) -#Es,Fs = mmff.sample_evalLJQs_ng4_PBC_atom_omp ( xs) -#Es,Fs = mmff.sample_evalLJQs_ng4_atom_omp ( xs) -Es,Fs = mmff.sample_evalLJQs_PBC_atom_omp ( xs) -#Es,Fs = mmff.sample_evalLJQs_atom_omp ( xs) -#Es,Fs = mmff.sample_addForce_Tricubic ( xs) -Fnum = numDeriv(xs,Es) #ratios = Fs[1:-1]/Fnum ;print("ratios", ratios) -plt.figure(); plt.plot(xs, Es, label="E"); plt.plot(xs, Fs, label="F_ana"); plt.plot(xs[1:-1], Fnum, label="F_num"); plt.grid(); plt.legend() -plt.show() -''' - -''' -# ----- Surf -xs = np.linspace(0.0,3.0,100) -mmff.init( xyz_name="data/pyridine", surf_name="data/NaCl_1x1_L2" ) -#Es,Fs = mmff.sample_addForce_Tricubic( xs) -#Es,Fs = mmff.sample_addForce(xs) -#Es,Fs = mmff.sample_evalMorsePBC_sym( xs) -#Es,Fs = mmff.sample_springbound( xs) -Fnum = numDeriv(xs,Es) #ratios = Fs[1:-1]/Fnum ;print("ratios", ratios) -plt.figure(); plt.plot(xs, Es, label="E"); plt.plot(xs, Fs, label="F_ana"); plt.plot(xs[1:-1], Fnum, label="F_num"); plt.grid(); plt.legend() -plt.show() -''' - -''' -# ----- Constrains -xs = np.linspace(0.0,3.0,100) -mmff.init( xyz_name="data/nHexadecan_dicarboxylic", constr_name="hexan-dicarboxylic.cons" ) -Es,Fs = mmff.sample_applyConstr( xs ) -Fnum = numDeriv(xs,Es) #ratios = Fs[1:-1]/Fnum ;print("ratios", ratios) -plt.figure(); plt.plot(xs, Es, label="E"); plt.plot(xs, Fs, label="F_ana"); plt.plot(xs[1:-1], Fnum, label="F_num"); plt.grid(); plt.legend() -plt.show() -''' - -''' -# ----- SetSwitches+run_omp -N = 1000 -xs = np.linspace(-3.0,3.0,N) -mmff.init( xyz_name="data/polymer-2_new", surf_name="data/NaCl_1x1_L2" ) -#mmff.init( xyz_name="data/O2")# ) -mmff.setSwitches(CheckInvariants=-1, PBC_nonBond=1, PBC_evalAtom=1, NonBonded=1, MMFF=1, doBonds=1, Angles=1, PiSigma=1, PiPiI=1, bNonBondNeighs=1, bSurfAtoms=1, bGridFF=-1, bTricubic=1, bConstrZ=-1, bConstrains=-1) -Es,Fs = mmff.sample_movementOfAtom(xs, ia=1) -Fnum = numDeriv(xs,Es) #ratios = Fs[1:-1]/Fnum ;print("ratios", ratios) -plt.figure(); plt.plot(xs, Es, label="E"); plt.plot(xs, Fs, label="F_ana"); plt.plot(xs[1:-1], Fnum, label="F_num"); plt.grid(); plt.legend() -plt.show() -''' -'''''' -# ----- lvecs -N = 10 -xs = np.linspace(1.0,3.0,N) -mmff.init( xyz_name="data/polymer-2_new", surf_name="data/NaCl_1x1_L2" ) -mmff.setSwitches(CheckInvariants=-1, PBC_nonBond=1, PBC_evalAtom=1, NonBonded=1, MMFF=1, doBonds=1, Angles=1, PiSigma=1, PiPiI=1, bNonBondNeighs=1, bSurfAtoms=-1, bGridFF=-1, bTricubic=-1, bConstrZ=-1, bConstrains=-1) -Es, Fs = mmff.sample_lvecs( xs ) -Fnum = numDeriv(xs,Es) #ratios = Fs[1:-1]/Fnum ;print("ratios", ratios) -plt.figure(); plt.plot(xs, Es, label="E"); plt.plot(xs, Fs, label="F_ana"); plt.plot(xs[1:-1], Fnum, label="F_num"); plt.grid(); plt.legend() -#print(Fs) -plt.show() - - -'''constr_name="hexan-dicarboxylic.cons", -# ----- DistConstr -#mmff.sample_DistConstr( xs, lmin=1, lmax=1, kmin=1, kmax=1, flim=1e+300, Es=None, Fs=None) -xs = np.linspace(0.0,3.0,100) -Es,Fs = mmff.sample_DistConstr( xs, lmin=0.9, lmax=1.2, flim=0.5 ) # ;print("Fs",Fs) -plt.figure(); plt.plot(xs, Es,'.-', label="E"); plt.plot(xs, Fs, label="F_ana"); plt.plot(xs[1:-1], numDeriv(xs,Es), label="F_num"); plt.grid(); plt.legend() -''' - -''' -# ----- SplineConstr -dx = 1.5 -x0 = 0.5 -Eps = np.array( [1.0, 0.0,-1.0,-0.5,-0.2,-0.1] ) -xp = (np.array(range(len(Eps)))-1)*dx + x0 -xs = np.linspace(0.0,6.0,100) -Es,Fs = mmff.sample_SplineConstr( xs, Eps, x0=x0, dx=dx ) # ;print("Fs",Fs) -plt.figure(); -plt.plot(xp, Eps, 'o-k', lw=0.2, label="Eps"); -plt.plot(xs, Es,'.-', label="E"); -plt.plot(xs, Fs,'-', label="F_ana"); -plt.plot(xs[1:-1],numDeriv(xs,Es),':', label="F_num"); -plt.grid(); plt.legend() - -plt.show()''' \ No newline at end of file