Source code for schrodinger.application.desmond.fep_schedule
"""This module handles the allocation of lambda windows for FEP simulations.
Different schemes of FEP simulations and interactions require different
schedules of lambda windows. The FepSchedule class abstracts the creation of
these schedules. This class is meant to be instantiated via its derived
classes (with mixins) or, preferably, using the get_fep_schedule function.
Copyright Schrodinger, LLC. All rights reserved.
"""
from typing import List
from typing import Optional
from typing import Tuple
from typing import Type
from schrodinger.application.desmond import constants
[docs]class FepSchedule:
"""Lambda schedule for FEP calculations.
This class should not be instantiated directly. Instead, subclass it
combining the mixins below.
Attributes
----------
n_win : int
Number of lambda windows
charge_div : float or None
Determines the fraction of lambda windows devoted to turning on
electrostatic interactions. This should be overriden by a mix-in
class or by the user.
bonded_div : float or None
Determines the fraction of lambda windows devoted to turning on
bonded interactions. This should be overriden by a mix-in class or by
the user.
"""
n_win: Optional[int] = None
charge_div: Optional[float] = None
bonded_div: Optional[float] = None
# The following three attributes are mean to hold the actual schedules.
_vdw: Optional[List[float]] = None
_coulomb: Optional[List[float]] = None
_bonded: Optional[List[float]] = None
[docs] def __init__(self, n_win: int, charge_div: Optional[float] = None) -> None:
if n_win <= 0:
raise ValueError(f'Invalid number of lambda windows: {n_win}')
assert self.charge_div is not None, \
(f'Invalid charge_div. The class f{self.__class__.__name__!r} should '
'not be instantiated directly')
assert self.bonded_div is not None, \
(f'Invalid bonded_div. The class f{self.__class__.__name__!r} should '
'not be instantiated directly')
self.n_win = n_win
if charge_div is not None:
self.charge_div = charge_div
n_coulomb = max(int(self.n_win / self.charge_div), 1) + 1
n_vdw = self.n_win - n_coulomb + 1
_coulomb = [0.0] * (self.n_win - n_coulomb)
coulomb = self._calc_charge_lambda(n_coulomb)
self._coulomb = _coulomb + coulomb
vdw_ = [1.0] * (self.n_win - n_vdw + 1)
vdw = self._calc_vdw_lambda(n_vdw)
self._vdw = vdw + vdw_
n_bonded = self.n_win - int(self.n_win / self.bonded_div)
bonded_ = [1.0] * (n_win - n_bonded - 1)
bonded = self._calc_bonded_lambda(n_bonded)
self._bonded = bonded + bonded_
@staticmethod
def _calc_vdw_lambda(n: int) -> List[float]:
"""Return lambda values for vdW potentials.
"""
def cubic_poly(x: float) -> float:
return 0.13418 * x - 0.03096 * x * x + 0.0037542 * x * x * x
return [cubic_poly(k * 8.0 / (n - 1)) for k in range(n - 1)]
@staticmethod
def _calc_charge_lambda(n: int) -> List[float]:
"""Return lambda values for Coulomb potentials.
"""
return [i / (n - 1) for i in range(n)]
@staticmethod
def _calc_bonded_lambda(n: int) -> List[float]:
"""Return lambda values for bonded potentials.
"""
return [i / (n - 1) for i in range(n)]
[docs] def get_lambda(self):
return self # For backwards-compatibility.
[docs]class FepScheduleAlchemical(FepSchedule):
@property
def vdw(self):
return [1.0] * self.n_win
@property
def coulomb(self):
return [0.0] * self.n_win
@property
def vdwA(self):
return list(reversed(self._vdw))
@property
def vdwB(self):
return self._vdw
@property
def chargeA(self):
return list(reversed(self._coulomb))
@property
def chargeB(self):
return self._coulomb
@property
def bondedA(self):
return list(reversed(self._bonded))
@property
def bondedB(self):
return self._bonded
[docs]class FepScheduleBinding(FepSchedule):
@property
def vdw(self):
return self._vdw
@property
def coulomb(self):
return self._coulomb
@property
def vdwA(self):
return [1.0] * self.n_win
@property
def vdwB(self):
return [1.0] * self.n_win
@property
def chargeA(self):
return [0.0] * self.n_win
@property
def chargeB(self):
return [0.0] * self.n_win
@property
def bondedA(self):
return [1.0] * self.n_win
@property
def bondedB(self):
return [1.0] * self.n_win
[docs]class FepSchemeMixIn:
"""Configure schedule of FEP schedule.
"""
charge_div = 5 / 2
bonded_div = float('inf')
[docs]class FepSchemeChargeMixIn(FepSchemeMixIn):
charge_div = 3 / 2
@staticmethod
def _calc_charge_lambda(n) -> List[float]:
"""Return a list of lambda values for Coulomb potentials.
Note that the computed schedule is quadratic.
"""
return [1.0 - (i / (n - 1) - 1.0)**2 for i in range(n)]
[docs]class FepSchemeQuickchargeMixIn(FepSchemeMixIn):
charge_div = 5.0
[docs]class FepSchemeSuperquickchargeMixIn(FepSchemeMixIn):
charge_div = 10.0
bonded_div = 10.0
[docs]def get_fep_schedule_class(fep_type: str = 'alchemical',
scheme: str = 'default') -> Type[FepSchedule]:
"""Return FEP schedule class.
See also
--------
get_fep_schedule
"""
base_classes = {
'alchemical': FepScheduleAlchemical,
'binding': FepScheduleBinding
}
mixin_classes = {
'default': FepSchemeMixIn,
'flexible': FepSchemeMixIn,
'charge': FepSchemeChargeMixIn,
'quickcharge': FepSchemeQuickchargeMixIn,
'superquickcharge': FepSchemeSuperquickchargeMixIn
}
class Schedule(mixin_classes[scheme], base_classes[fep_type]):
pass
return Schedule
[docs]def get_fep_schedule(n_win: int,
fep_type: str = 'alchemical',
scheme: str = 'default') -> FepSchedule:
"""Instantiate FEP schedule.
Return a class encapsulating a concrete lambda schedule for the specified
FEP simulation type and scheme.
:param n_win: Number of lambda windows.
:type n_win: int
:param fep_type: Type of FEP simulation. Can be either 'alchemical'
(default) or 'binding'.
:type fep_type: str
:param scheme: Simulation scheme. Can be one of 'default', 'flexible',
'charge', 'quickcharge', or 'superquickcharge'.
:type scheme: str
:return: Concrete schedule of the specified type, scheme, and number of
windows.
:rtype: FepSchedule
"""
FepSchedule = get_fep_schedule_class(fep_type, scheme)
return FepSchedule(n_win)
[docs]class FepScheduleAlchemicalDefault(FepSchemeMixIn, FepScheduleAlchemical):
"""Convenience class for the alchemical schedule with the default scheme.
"""
[docs]class FepScheduleAlchemicalCharge(FepSchemeChargeMixIn, FepScheduleAlchemical):
"""Convenience class for the alchemical schedule with the charge scheme.
"""
[docs]class FepScheduleBindingDefault(FepSchemeMixIn, FepScheduleBinding):
"""Convenience class for the binding schedule with the default scheme.
"""
[docs]def get_absolute_binding_num_lambda_windows(
protocol: constants.SIMULATION_PROTOCOL,
restrained: bool) -> Tuple[int, int]:
"""
Get the number of lambda windows in the absolute binding lambda
schedules. These are derived from the lambda schedules in
mmshare/data/desmond/abfep
:return: A tuple of # of complex windows, # of solvent windows
"""
if restrained:
if protocol == constants.SIMULATION_PROTOCOL.CHARGED:
return 128, 68
else:
return 80, 68
else:
if protocol == constants.SIMULATION_PROTOCOL.CHARGED:
return 108, 60
else:
return 68, 60