From 8735077012556128a1b52affe3847011148743ff Mon Sep 17 00:00:00 2001 From: lipengfei Date: Mon, 14 Oct 2024 22:11:19 +0800 Subject: [PATCH 1/2] fix analytic bug --- abfe/abfe_alchemy_md.py | 15 ++++++++++++++- abfe/md/analytic.py | 21 +++++++++++++++++---- 2 files changed, 31 insertions(+), 5 deletions(-) diff --git a/abfe/abfe_alchemy_md.py b/abfe/abfe_alchemy_md.py index 5678e48..2397960 100644 --- a/abfe/abfe_alchemy_md.py +++ b/abfe/abfe_alchemy_md.py @@ -502,13 +502,26 @@ def _add_boresch_restraint(args, prmtop, lig_masks, rec_masks, lig_idx, rst_refs temperature = args.temperature kr = args.bond_const ktheta = args.angle_const + kalpha = args.angle_const kphi = args.dihedral_const + kbeta = args.dihedral_const + kgamma = args.dihedral_const rst_correction = analytic.restraint_correction( - T=temperature, r0=r0, theta0=theta0, kr=kr, ktheta=ktheta, kphi=kphi + T=temperature, r0=r0, theta0=theta0, + kr=kr, ktheta=ktheta, kalpha=kalpha, kphi=kphi, kbeta=kbeta, kgamma=kgamma ) results = np.array([[rst_correction]]).T df = pd.DataFrame(results, columns=['Total (kcal/mol)']) df.to_csv('results.csv', float_format='%.2f') + # Just dump out the correction results from schrodinger for comparison. + _rst_correction = analytic.restraint_correction_schrodinger( + T=temperature, r0=r0, alpha0=alpha0, theta0=theta0, + gamma0=gamma0, beta0=beta0, phi0=phi0, + Kr=args.bond_const, Kang=args.angle_const, Kdihed=args.dihedral_const + ) + _results = np.array([[_rst_correction]]).T + _df = pd.DataFrame(_results, columns=['Total (kcal/mol)']) + _df.to_csv('results_schrodinger.csv', float_format='%.2f') #========================================================================== # 2.2 VDW #========================================================================== diff --git a/abfe/md/analytic.py b/abfe/md/analytic.py index b4a6676..4756789 100644 --- a/abfe/md/analytic.py +++ b/abfe/md/analytic.py @@ -12,7 +12,7 @@ """ from math import pi, log, sin, sqrt, exp, erf -def restraint_correction(T, r0, theta0, kr, ktheta, kphi): +def restraint_correction(T, r0, theta0, kr, ktheta, kalpha, kphi, kbeta, kgamma): """Do Boresch-style restraint analytic correction, refer to Calculation of Standard Binding Free Energies: Aromatic Molecules @@ -29,7 +29,10 @@ def restraint_correction(T, r0, theta0, kr, ktheta, kphi): theta0 (float): equilibrium value for angle theta kr (float): force constant for distance r ktheta (float): force constant for angle theta + kalpha (float): force constant for angle alpha kphi (float): force constant for dihedral phi + kbeta (float): force constant for dihedral beta + kgamma (float): force constant for dihedral gamma Returns: FtC0 + Fr (float): Boresch analytic restraint correction value @@ -39,13 +42,16 @@ def restraint_correction(T, r0, theta0, kr, ktheta, kphi): theta0 = float(theta0) kr = float(kr)*2.0 ktheta = float(ktheta)*2.0 + kalpha = float(kalpha)*2.0 kphi = float(kphi)*2.0 + kbeta = float(kbeta)*2.0 + kgamma = float(kgamma)*2.0 V = 1660. # in A^3 kB = 0.0019872041 # Boltzmann's constant (kcal/mol/K) FtC0 = -(kB*T)*log(r0**2*sin(theta0/180.0*pi)\ *sqrt((2*pi*kB*T)**3/(kr*ktheta*kphi))/V) # Eq. 38 - Fr = -(kB*T)*log(1.0/(8*pi**2)*(2*pi*kB*T/kr)**1.5) # Eq. 40 + Fr = -(kB*T)*log(1.0/(8*pi**2)*sqrt((2*pi*kB*T)**3/(kalpha*kbeta*kgamma))) # Eq. 40 print('Translational: %.2f kcal/mol' % FtC0) print('Rotational: %.2f kcal/mol' % Fr) @@ -108,6 +114,13 @@ def restraint_correction_schrodinger(T, r0, alpha0, theta0, gamma0, beta0, phi0, """Just do test for restraint_correction function. """ correction = restraint_correction( - T=298.15, ktheta=100., theta0=85.370, kphi=100., kr=10., r0=7.230, + T=298.15, r0=4.470, theta0=101.230, + kr=10, ktheta=100, kalpha=100, kphi=100, kbeta=100, kgamma=100 ) - print('The reference value is 8.99 kcal/mol, while this calculated value is %.2f kcal/mol.' % correction) \ No newline at end of file + print('The reference value is 11.62 kcal/mol, while this calculated restraint_correction is %.2f kcal/mol.' % correction) + correction = restraint_correction_schrodinger( + T=298.15, r0=4.470, alpha0=72.140, theta0=101.230, + gamma0=77.887, beta0=-137.094, phi0=9.048, + Kr=10, Kang=100, Kdihed=100 + ) + print('The reference value is 11.65 kcal/mol, while this calculated restraint_correction_schrodinger is %.2f kcal/mol.' % correction) \ No newline at end of file From 8d8eddc5c2cb76f349b448d3af07edd4556a2698 Mon Sep 17 00:00:00 2001 From: lipengfei Date: Tue, 15 Oct 2024 10:11:48 +0800 Subject: [PATCH 2/2] updates on dG values --- testing/abfe_testing/dG_components_results.csv | 5 ++--- testing/abfe_testing/dG_results.csv | 5 ++--- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/testing/abfe_testing/dG_components_results.csv b/testing/abfe_testing/dG_components_results.csv index 44eae69..b7d6546 100644 --- a/testing/abfe_testing/dG_components_results.csv +++ b/testing/abfe_testing/dG_components_results.csv @@ -1,4 +1,3 @@ ,MBAR_dG_solvated_vdw,MBAR_dG_solvated_cc,TI_dG_solvated_vdw,TI_dG_solvated_cc,MBAR_dG_complex_vdw,MBAR_dG_complex_cc,MBAR_dG_complex_restraint,MBAR_dG_analytic,TI_dG_complex_vdw,TI_dG_complex_cc,TI_dG_complex_restraint,TI_dG_analytic -bnz,1.06,0.00,1.04,0.00,7.92,0.00,-6.11,9.57,7.75,0.00,-6.71,9.57 -,MBAR_dG_solvated_vdw,MBAR_dG_solvated_cc,TI_dG_solvated_vdw,TI_dG_solvated_cc,MBAR_dG_complex_vdw,MBAR_dG_complex_cc,MBAR_dG_complex_restraint,MBAR_dG_analytic,TI_dG_complex_vdw,TI_dG_complex_cc,TI_dG_complex_restraint,TI_dG_analytic -phn,5.18,0.00,5.12,0.00,9.66,0.00,-6.06,9.61,9.74,0.00,-6.27,9.61 +phn,5.18,0.00,5.12,0.00,9.66,0.00,-6.06,11.63,9.74,0.00,-6.27,11.63 +bnz,1.06,0.00,1.04,0.00,7.92,0.00,-6.11,11.64,7.75,0.00,-6.71,11.64 diff --git a/testing/abfe_testing/dG_results.csv b/testing/abfe_testing/dG_results.csv index 6a36f1e..96218be 100644 --- a/testing/abfe_testing/dG_results.csv +++ b/testing/abfe_testing/dG_results.csv @@ -1,4 +1,3 @@ ,MBAR_dG,MBAR_dG_std,TI_dG,TI_dG_std -bnz,-3.40,0.08,-3.85,0.10 -,MBAR_dG,MBAR_dG_std,TI_dG,TI_dG_std -phn,-0.94,0.10,-1.28,0.10 +phn,1.08,0.10,0.74,0.10 +bnz,-1.33,0.08,-1.78,0.10