updated protocols

freeenergylab · Aug 12, 2024 · 6f5fdc6 · 6f5fdc6
1 parent dfa85b7
commit 6f5fdc6
Show file tree

Hide file tree

Showing 4 changed files with 64 additions and 68 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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)
+    print('The reference value is 8.99 kcal/mol, while this calculated value is %.2f kcal/mol.' % correction)