diff --git a/abfe/abfe_alchemy_md.py b/abfe/abfe_alchemy_md.py index c577397..77890d5 100644 --- a/abfe/abfe_alchemy_md.py +++ b/abfe/abfe_alchemy_md.py @@ -134,18 +134,18 @@ def _parse_args(): help='Treatment of the non-bonded interactions of SC internal and SC-CC cross regions' ) # Boresch-style restraint options - boresch_grp = parser.add_argument_group('Boresch restraints') + boresch_grp = parser.add_argument_group('Boresch-style restraints') boresch_grp.add_argument( '-bc', '--bond_const', default=10.0, type=float, - help='Weight of bond restraints during ABFE' + help='Force constant for bond restraint, in kcal/mol/A^2' ) boresch_grp.add_argument( '-ac', '--angle_const', default=100.0, type=float, - help='Weight of angle restraints during ABFE' + help='Force constant for angle restraint, in kcal/mol/radian^2' ) boresch_grp.add_argument( '-dc', '--dihedral_const', default=100.0, type=float, - help='Weight of dihedral restraints during ABFE' + help='Force constant for dihedral restraint, in kcal/mol' ) boresch_grp.add_argument( '-rs', '--restraint_seed', type=int, default=None, @@ -311,7 +311,9 @@ def _add_boresch_restraint(args, prmtop, lig_masks, rec_masks, lig_idx, rst_refs .replace('TIMASK1', str(lig.timask1))\ .replace('TIMASK2', str(lig.timask2))\ .replace('SCMASK1', str(lig.scmask1))\ - .replace('SCMASK2', str(lig.scmask2)) + .replace('SCMASK2', str(lig.scmask2))\ + .replace('RESTRAINT_WT', str(args.restraint_wt))\ + .replace('RESTRAINTMASK', str(args.restraintmask)) for n, (tempi, temp0) in enumerate(zip(args.heat_temps[:-1], args.heat_temps[1:]), 1): lig_md_steps['heat-%s'%str(n)] = SOLVATED_HEAT\ .replace('CLAMBDA', str(_lambda))\ @@ -499,10 +501,10 @@ def _add_boresch_restraint(args, prmtop, lig_masks, rec_masks, lig_idx, rst_refs ) temperature = args.temperature kr = args.bond_const - kphi = args.dihedral_const ktheta = args.angle_const + kphi = args.dihedral_const rst_correction = analytic.restraint_correction( - T=temperature, ktheta=ktheta, theta0=theta0, kphi=kphi, kr=kr, r0=r0 + T=temperature, r0=r0, theta0=theta0, kr=kr, ktheta=ktheta, kphi=kphi ) results = np.array([[rst_correction]]).T df = pd.DataFrame(results, columns=['Total (kcal/mol)']) @@ -533,7 +535,9 @@ def _add_boresch_restraint(args, prmtop, lig_masks, rec_masks, lig_idx, rst_refs .replace('TIMASK1', str(comp.timask1))\ .replace('TIMASK2', str(comp.timask2))\ .replace('SCMASK1', str(comp.scmask1))\ - .replace('SCMASK2', str(comp.scmask2)) + .replace('SCMASK2', str(comp.scmask2))\ + .replace('RESTRAINT_WT', str(args.restraint_wt))\ + .replace('RESTRAINTMASK', str(args.restraintmask)) for n, (tempi, temp0) in enumerate(zip(args.heat_temps[:-1], args.heat_temps[1:]), 1): comp_md_steps['heat-%s'%str(n)] = COMPLEX_HEAT\ .replace('CLAMBDA', str(_lambda))\ @@ -700,7 +704,9 @@ def _add_boresch_restraint(args, prmtop, lig_masks, rec_masks, lig_idx, rst_refs comp_md_steps['min'] = COMPLEX_RESTRAINT_MIN\ .replace('CLAMBDA', _lambda)\ .replace('CUT', str(args.cutoff))\ - .replace('GTI_ADD_SC', str(args.gti_add_sc)) + .replace('GTI_ADD_SC', str(args.gti_add_sc))\ + .replace('RESTRAINT_WT', str(args.restraint_wt))\ + .replace('RESTRAINTMASK', str(args.restraintmask)) for n, (tempi, temp0) in enumerate(zip(args.heat_temps[:-1], args.heat_temps[1:]), 1): comp_md_steps['heat-%s'%str(n)] = COMPLEX_RESTRAINT_HEAT\ .replace('CLAMBDA', str(_lambda))\ diff --git a/abfe/md/amber_alchemy_mdin.py b/abfe/md/amber_alchemy_mdin.py index 9a1a8e4..e36ec83 100644 --- a/abfe/md/amber_alchemy_mdin.py +++ b/abfe/md/amber_alchemy_mdin.py @@ -16,15 +16,16 @@ &cntrl imin = 1, ntmin = 2, - maxcyc = 10000, - ncyc = 10000, - drms = 0.01, - drms = 0.0001, + maxcyc = 200, + ncyc = 200, + drms = 0.1, + dx0 = 0.005, ntb = 1, cut = CUT, - ntr = 0, - restraint_wt = 0.0, + ntr = 1, + restraint_wt = RESTRAINT_WT, + restraintmask = "RESTRAINTMASK", ntp = 0, pres0 = 1.0, @@ -48,7 +49,7 @@ scmask2 = SCMASK2, scalpha = 0.2, scbeta = 50.0, - + gti_lam_sch = 1, gti_ele_sc = 1, gti_vdw_sc = 1, @@ -105,7 +106,7 @@ scmask2 = SCMASK2, scalpha = 0.2, scbeta = 50.0, - + gti_lam_sch = 1, gti_ele_sc = 1, gti_vdw_sc = 1, @@ -172,7 +173,7 @@ scmask2 = SCMASK2, scalpha = 0.2, scbeta = 50.0, - + gti_lam_sch = 1, gti_ele_sc = 1, gti_vdw_sc = 1, @@ -230,7 +231,7 @@ scmask2 = SCMASK2, scalpha = 0.2, scbeta = 50.0, - + gti_lam_sch = 1, gti_ele_sc = 1, gti_vdw_sc = 1, @@ -253,15 +254,16 @@ &cntrl imin = 1, ntmin = 2, - maxcyc = 10000, - ncyc = 10000, - drms = 0.01, - drms = 0.0001, + maxcyc = 200, + ncyc = 200, + drms = 0.1, + dx0 = 0.005, ntb = 1, cut = CUT, - ntr = 0, - restraint_wt = 0.0, + ntr = 1, + restraint_wt = RESTRAINT_WT, + restraintmask = "RESTRAINTMASK", ntp = 0, pres0 = 1.0, @@ -285,7 +287,7 @@ scmask2 = SCMASK2, scalpha = 0.2, scbeta = 50.0, - + gti_lam_sch = 1, gti_ele_sc = 1, gti_vdw_sc = 1, @@ -313,7 +315,7 @@ ntx = 1, nstlim = 10000, dt = DT, - vlimit = 10, + !vlimit = 10, ntb = 1, cut = CUT, @@ -349,7 +351,7 @@ scmask2 = SCMASK2, scalpha = 0.2, scbeta = 50.0, - + gti_lam_sch = 1, gti_ele_sc = 1, gti_vdw_sc = 1, @@ -382,7 +384,7 @@ ntx = 5, nstlim = 10000, dt = DT, - vlimit = 10, + !vlimit = 10, ntb = 2, cut = CUT, @@ -481,7 +483,7 @@ scmask2 = SCMASK2, scalpha = 0.2, scbeta = 50.0, - + gti_lam_sch = 1, gti_ele_sc = 1, gti_vdw_sc = 1, @@ -514,15 +516,16 @@ &cntrl imin = 1, ntmin = 2, - maxcyc = 10000, - ncyc = 10000, - drms = 0.01, - drms = 0.0001, + maxcyc = 200, + ncyc = 200, + drms = 0.1, + dx0 = 0.005, ntb = 1, cut = CUT, - ntr = 0, - restraint_wt = 0.0, + ntr = 1, + restraint_wt = RESTRAINT_WT, + restraintmask = "RESTRAINTMASK", ntp = 0, pres0 = 1.0, @@ -558,7 +561,7 @@ ntx = 1, nstlim = 10000, dt = DT, - vlimit = 20, + !vlimit = 20, ntb = 1, cut = CUT, @@ -611,7 +614,7 @@ ntx = 5, nstlim = 10000, dt = DT, - vlimit = 20, + !vlimit = 20, ntb = 2, cut = CUT, diff --git a/abfe/md/amber_mdin.py b/abfe/md/amber_mdin.py index f8ee60b..e60dced 100644 --- a/abfe/md/amber_mdin.py +++ b/abfe/md/amber_mdin.py @@ -17,8 +17,6 @@ ntmin = 1, maxcyc = 10000, ncyc = 5000, - dx0 = 0.01, - drms = 0.0001, cut = CUT, ntr = 1, @@ -33,8 +31,6 @@ ntmin = 1, maxcyc = 10000, ncyc = 5000, - dx0 = 0.01, - drms = 0.0001, cut = CUT, ntr = 0, @@ -102,7 +98,7 @@ restraintmask = "RESTRAINTMASK", ntp = 1, - pres0 = 1.0, + pres0 = 1.0, taup = 2.0, barostat = 2, @@ -159,8 +155,6 @@ ntmin = 1, maxcyc = 10000, ncyc = 5000, - dx0 = 0.01, - drms = 0.0001, cut = CUT, ntr = 1, @@ -175,8 +169,6 @@ ntmin = 1, maxcyc = 10000, ncyc = 5000, - dx0 = 0.01, - drms = 0.0001, cut = CUT, ntr = 0, diff --git a/abfe/md/analytic.py b/abfe/md/analytic.py index 20bc0d9..b4a6676 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, ktheta, theta0, kphi, kr, r0): +def restraint_correction(T, r0, theta0, kr, ktheta, kphi): """Do Boresch-style restraint analytic correction, refer to Calculation of Standard Binding Free Energies: Aromatic Molecules @@ -25,20 +25,21 @@ def restraint_correction(T, ktheta, theta0, kphi, kr, r0): Args: T (float): system simulation temperature - ktheta (float): force constant for theta - theta0 (float): equilibrium value for theta - kphi (float): force constant for phi - kr (float): force constant for r - r0 (float): equilibrium value for r + r0 (float): equilibrium value for distance r + theta0 (float): equilibrium value for angle theta + kr (float): force constant for distance r + ktheta (float): force constant for angle theta + kphi (float): force constant for dihedral phi Returns: FtC0 + Fr (float): Boresch analytic restraint correction value """ - ktheta = float(ktheta)*2.0 + T = float(T) + r0 = float(r0) theta0 = float(theta0) - kphi = float(kphi)*2.0 kr = float(kr)*2.0 - r0 = float(r0) + ktheta = float(ktheta)*2.0 + kphi = float(kphi)*2.0 V = 1660. # in A^3 kB = 0.0019872041 # Boltzmann's constant (kcal/mol/K) @@ -61,6 +62,7 @@ def restraint_correction_schrodinger(T, r0, alpha0, theta0, gamma0, beta0, phi0, J Chem. Inf. Model. 2023, 63, 3171-3185 https://pubs.acs.org/doi/10.1021/acs.jcim.3c00013 + https://chemrxiv.org/engage/chemrxiv/article-details/63a23f6116e9a872d32f81ef Args: T (float): system simulation temperature @@ -97,22 +99,15 @@ def restraint_correction_schrodinger(T, r0, alpha0, theta0, gamma0, beta0, phi0, Z_ang_theta_r = sqrt(pi/(beta*Kang))*exp(-1./(4.*beta*Kang))*sin(theta0/180.*pi) # Eq. 5 Z_dihed_r = sqrt(pi/(beta*Kdihed))*erf(pi*sqrt(beta*Kdihed)) # Eq. 6 - # result = -Z_ang_theta_r*Z_ang_alpha_r*kB*T*log(Z_dist_r*(Z_dihed_r**3)/(8.*(pi**2)*V)) # Eq. 3 - result = -kB*T*log(Z_dist_r*Z_ang_theta_r*Z_ang_alpha_r*(Z_dihed_r**3)/(8.*(pi**2)*V)) # Eq. 3 + # result = -Z_ang_theta_r*Z_ang_alpha_r*kB*T*log(Z_dist_r*(Z_dihed_r**3)/(8.*(pi**2)*V)) # JCIM Eq. 3 + result = -kB*T*log(Z_dist_r*Z_ang_theta_r*Z_ang_alpha_r*(Z_dihed_r**3)/(8.*(pi**2)*V)) # ChemRxiv Eq. 3 return result if __name__ == '__main__': - """Just one test for restraint_correction function. + """Just do test for restraint_correction function. """ correction = restraint_correction( - T=300., ktheta=100., theta0=100.470, kphi=100., kr=10., r0=4.350 - ) - print('The reference value is 9.65 kcal/mol, while this calculated value is %.2f kcal/mol.' % correction) - - correction = restraint_correction_schrodinger( - T=300, r0=4.350, alpha0=75.020, theta0=100.470, - gamma0=79.359, beta0=172.460, phi0=5.991, - Kr=10., Kang=100., Kdihed=100. + T=298.15, ktheta=100., theta0=85.370, kphi=100., kr=10., r0=7.230, ) - print('The reference value is 9.65 kcal/mol, while this calculated value is %.2f kcal/mol.' % correction) \ No newline at end of file + print('The reference value is 8.99 kcal/mol, while this calculated value is %.2f kcal/mol.' % correction) \ No newline at end of file