# Path: openmmtools/forces.py
# class MultipleForcesError(Exception):
# class NoForceFoundError(Exception):
# class RadiallySymmetricRestraintForce(utils.RestorableOpenMMObject):
# class RadiallySymmetricCentroidRestraintForce(RadiallySymmetricRestraintForce,
#                                               openmm.CustomCentroidBondForce):
# class RadiallySymmetricBondRestraintForce(RadiallySymmetricRestraintForce,
#                                           openmm.CustomBondForce):
# class HarmonicRestraintForceMixIn(object):
# class HarmonicRestraintForce(HarmonicRestraintForceMixIn,
#                              RadiallySymmetricCentroidRestraintForce):
# class HarmonicRestraintBondForce(HarmonicRestraintForceMixIn,
#                                  RadiallySymmetricBondRestraintForce):
# class FlatBottomRestraintForceMixIn(object):
# class FlatBottomRestraintForce(FlatBottomRestraintForceMixIn,
#                                RadiallySymmetricCentroidRestraintForce):
# class FlatBottomRestraintBondForce(FlatBottomRestraintForceMixIn,
#                                    RadiallySymmetricBondRestraintForce):
# class UnshiftedReactionFieldForce(openmm.CustomNonbondedForce):
# class SwitchedReactionFieldForce(openmm.CustomNonbondedForce):
# def iterate_forces(system):
# def find_forces(system, force_type, only_one=False, include_subclasses=False):
# def _compute_sphere_volume(radius):
# def _compute_harmonic_volume(radius, spring_constant, beta):
# def _compute_harmonic_radius(spring_constant, potential_energy):
#     def __init__(self, restraint_parameters, restrained_atom_indices1,
#                  restrained_atom_indices2, controlling_parameter_name,
#                  *args, **kwargs):
#     def _create_bond(self, bond_parameter_values, restrained_atom_indices1,
#                      restrained_atom_indices2):
#     def restrained_atom_indices1(self):
#     def restrained_atom_indices2(self):
#     def restraint_parameters(self):
#     def controlling_parameter_name(self):
#     def distance_at_energy(self, potential_energy):
#     def compute_standard_state_correction(self, thermodynamic_state, square_well=False,
#                                           radius_cutoff=None, energy_cutoff=None,
#                                           max_volume=None):
#     def _compute_restraint_volume(self, thermodynamic_state, square_well,
#                                   radius_cutoff, energy_cutoff):
#     def _integrate_restraint_volume(self, thermodynamic_state, square_well,
#                                     radius_cutoff, energy_cutoff):
#         def restraint_potential_func(r):
#         def integrand(r):
#     def _determine_integral_limits(self, thermodynamic_state, radius_cutoff,
#                                    energy_cutoff, potential_energy_func):
#     def __init__(self, energy_function, restraint_parameters,
#                  restrained_atom_indices1, restrained_atom_indices2,
#                  controlling_parameter_name='lambda_restraints'):
#     def restrained_atom_indices1(self):
#     def restrained_atom_indices1(self, atom_indices):
#     def restrained_atom_indices2(self):
#     def restrained_atom_indices2(self, atom_indices):
#     def _create_bond(self, bond_parameter_values, restrained_atom_indices1,
#                      restrained_atom_indices2):
#     def __init__(self, energy_function, restraint_parameters,
#                  restrained_atom_index1, restrained_atom_index2,
#                  controlling_parameter_name='lambda_restraints'):
#     def restrained_atom_indices1(self):
#     def restrained_atom_indices1(self, atom_indices):
#     def restrained_atom_indices2(self):
#     def restrained_atom_indices2(self, atom_indices):
#     def _create_bond(self, bond_parameter_values, restrained_atom_indices1, restrained_atom_indices2):
#     def __init__(self, spring_constant, *args, **kwargs):
#     def spring_constant(self):
#     def distance_at_energy(self, potential_energy):
#     def _compute_restraint_volume(self, thermodynamic_state, square_well,
#                                   radius_cutoff, energy_cutoff):
#     def __init__(self, spring_constant, well_radius, *args, **kwargs):
#     def spring_constant(self):
#     def well_radius(self):
#     def distance_at_energy(self, potential_energy):
#     def _compute_restraint_volume(self, thermodynamic_state, square_well,
#                                   radius_cutoff, energy_cutoff):
#     def _determine_integral_limits(self, thermodynamic_state, radius_cutoff,
#                                    energy_cutoff, potential_energy_func=None):
#     def __init__(self, cutoff_distance=15*unit.angstroms, switch_width=1.0*unit.angstrom,
#                  reaction_field_dielectric=78.3):
#     def from_nonbonded_force(cls, nonbonded_force, switch_width=1.0*unit.angstrom):
#     def from_system(cls, system, switch_width=1.0*unit.angstrom):
#     def __init__(self, cutoff_distance=15*unit.angstroms, switch_width=1.0*unit.angstrom,
#                  reaction_field_dielectric=78.3):
#     def from_nonbonded_force(cls, nonbonded_force, switch_width=1.0*unit.angstrom):
#     def from_system(cls, system, switch_width=1.0*unit.angstrom):


import copy
import numpy as np
import mdtraj
    import openmm
    import doctest
    from openmm import unit
    from simtk import openmm, unit
from openmmtools import forces
    Note that `AbsoluteAlchemicalFactory.create_alchemical_system()` can NOT
    handle the resulting `System` object yet since the `CustomNonbondedForce`
    are not recognized and re-coded.

    Parameters
    ----------
    reference_system : openmm.System
        The system to use as a reference. This will be modified only if
        `return_copy` is `False`.
    switch_width : openmm.unit.Quantity, default=1.0*angstrom
        Switch width for electrostatics (units of distance).
    return_copy : bool, optional, default=True
        If `True`, `reference_system` is not modified, and a copy is returned.
        Setting it to `False` speeds up the function execution but modifies
        the `reference_system` object.
    shifted : bool, optional, default=False
        If `True`, a shifted reaction-field will be used.

    Returns
    -------
    system : openmm.System
        System with reaction-field converted to c_rf = 0

    """
    if return_copy:
        system = copy.deepcopy(reference_system)
    else:
        system = reference_system

    if shifted:

# Path: openmmtools/forces.py
# class MultipleForcesError(Exception):
# class NoForceFoundError(Exception):
# class RadiallySymmetricRestraintForce(utils.RestorableOpenMMObject):
# class RadiallySymmetricCentroidRestraintForce(RadiallySymmetricRestraintForce,
#                                               openmm.CustomCentroidBondForce):
# class RadiallySymmetricBondRestraintForce(RadiallySymmetricRestraintForce,
#                                           openmm.CustomBondForce):
# class HarmonicRestraintForceMixIn(object):
# class HarmonicRestraintForce(HarmonicRestraintForceMixIn,
#                              RadiallySymmetricCentroidRestraintForce):
# class HarmonicRestraintBondForce(HarmonicRestraintForceMixIn,
#                                  RadiallySymmetricBondRestraintForce):
# class FlatBottomRestraintForceMixIn(object):
# class FlatBottomRestraintForce(FlatBottomRestraintForceMixIn,
#                                RadiallySymmetricCentroidRestraintForce):
# class FlatBottomRestraintBondForce(FlatBottomRestraintForceMixIn,
#                                    RadiallySymmetricBondRestraintForce):
# class UnshiftedReactionFieldForce(openmm.CustomNonbondedForce):
# class SwitchedReactionFieldForce(openmm.CustomNonbondedForce):
# def iterate_forces(system):
# def find_forces(system, force_type, only_one=False, include_subclasses=False):
# def _compute_sphere_volume(radius):
# def _compute_harmonic_volume(radius, spring_constant, beta):
# def _compute_harmonic_radius(spring_constant, potential_energy):
#     def __init__(self, restraint_parameters, restrained_atom_indices1,
#                  restrained_atom_indices2, controlling_parameter_name,
#                  *args, **kwargs):
#     def _create_bond(self, bond_parameter_values, restrained_atom_indices1,
#                      restrained_atom_indices2):
#     def restrained_atom_indices1(self):
#     def restrained_atom_indices2(self):
#     def restraint_parameters(self):
#     def controlling_parameter_name(self):
#     def distance_at_energy(self, potential_energy):
#     def compute_standard_state_correction(self, thermodynamic_state, square_well=False,
#                                           radius_cutoff=None, energy_cutoff=None,
#                                           max_volume=None):
#     def _compute_restraint_volume(self, thermodynamic_state, square_well,
#                                   radius_cutoff, energy_cutoff):
#     def _integrate_restraint_volume(self, thermodynamic_state, square_well,
#                                     radius_cutoff, energy_cutoff):
#         def restraint_potential_func(r):
#         def integrand(r):
#     def _determine_integral_limits(self, thermodynamic_state, radius_cutoff,
#                                    energy_cutoff, potential_energy_func):
#     def __init__(self, energy_function, restraint_parameters,
#                  restrained_atom_indices1, restrained_atom_indices2,
#                  controlling_parameter_name='lambda_restraints'):
#     def restrained_atom_indices1(self):
#     def restrained_atom_indices1(self, atom_indices):
#     def restrained_atom_indices2(self):
#     def restrained_atom_indices2(self, atom_indices):
#     def _create_bond(self, bond_parameter_values, restrained_atom_indices1,
#                      restrained_atom_indices2):
#     def __init__(self, energy_function, restraint_parameters,
#                  restrained_atom_index1, restrained_atom_index2,
#                  controlling_parameter_name='lambda_restraints'):
#     def restrained_atom_indices1(self):
#     def restrained_atom_indices1(self, atom_indices):
#     def restrained_atom_indices2(self):
#     def restrained_atom_indices2(self, atom_indices):
#     def _create_bond(self, bond_parameter_values, restrained_atom_indices1, restrained_atom_indices2):
#     def __init__(self, spring_constant, *args, **kwargs):
#     def spring_constant(self):
#     def distance_at_energy(self, potential_energy):
#     def _compute_restraint_volume(self, thermodynamic_state, square_well,
#                                   radius_cutoff, energy_cutoff):
#     def __init__(self, spring_constant, well_radius, *args, **kwargs):
#     def spring_constant(self):
#     def well_radius(self):
#     def distance_at_energy(self, potential_energy):
#     def _compute_restraint_volume(self, thermodynamic_state, square_well,
#                                   radius_cutoff, energy_cutoff):
#     def _determine_integral_limits(self, thermodynamic_state, radius_cutoff,
#                                    energy_cutoff, potential_energy_func=None):
#     def __init__(self, cutoff_distance=15*unit.angstroms, switch_width=1.0*unit.angstrom,
#                  reaction_field_dielectric=78.3):
#     def from_nonbonded_force(cls, nonbonded_force, switch_width=1.0*unit.angstrom):
#     def from_system(cls, system, switch_width=1.0*unit.angstrom):
#     def __init__(self, cutoff_distance=15*unit.angstroms, switch_width=1.0*unit.angstrom,
#                  reaction_field_dielectric=78.3):
#     def from_nonbonded_force(cls, nonbonded_force, switch_width=1.0*unit.angstrom):
#     def from_system(cls, system, switch_width=1.0*unit.angstrom):


# Path: openmmtools/forcefactories.py
import copy
import numpy as np
import mdtraj
    import openmm
    import doctest
    from openmm import unit
    from simtk import openmm, unit
from openmmtools import forces

    Note that `AbsoluteAlchemicalFactory.create_alchemical_system()` can NOT
    handle the resulting `System` object yet since the `CustomNonbondedForce`
    are not recognized and re-coded.

    Parameters
    ----------
    reference_system : openmm.System
        The system to use as a reference. This will be modified only if
        `return_copy` is `False`.
    switch_width : openmm.unit.Quantity, default=1.0*angstrom
        Switch width for electrostatics (units of distance).
    return_copy : bool, optional, default=True
        If `True`, `reference_system` is not modified, and a copy is returned.
        Setting it to `False` speeds up the function execution but modifies
        the `reference_system` object.
    shifted : bool, optional, default=False
        If `True`, a shifted reaction-field will be used.

    Returns
    -------
    system : openmm.System
        System with reaction-field converted to c_rf = 0

    """
    if return_copy:
        system = copy.deepcopy(reference_system)
    else:
        system = reference_system

    if shifted:

        force_constructor = getattr(forces, 'SwitchedReactionFieldForce')