Can I use a customnonbondedforce to manipulate the interactions in an alchemical system?

[Nithishwer]

Hi,

I am trying to do GCMC water sampling on my alchemical system, where a ligand is partially turned off. I was wondering if I can use openmmtools's alchemy module along with grand for this. This is how grand manipulates water interactions when doing the GCMC/NCMC process:

    def customiseForces(self):
        """
        Create a CustomNonbondedForce to handle water-water interactions and modify the original NonbondedForce
        to ignore water interactions
        """
        #  Need to make sure that the electrostatics are handled using PME (for now)
        if self.nonbonded_force.getNonbondedMethod() != openmm.NonbondedForce.PME:
            self.raiseError("Currently only supporting PME for long range electrostatics")

        # Define the energy expression for the softcore sterics
        energy_expression = ("U;"
                             "U = (lambda^soft_a) * 4 * epsilon * x * (x-1.0);"  # Softcore energy
                             "x = (sigma/reff)^6;"  # Define x as sigma/r(effective)
                             # Calculate effective distance
                             "reff = sigma*((soft_alpha*(1.0-lambda)^soft_b + (r/sigma)^soft_c))^(1/soft_c);"
                             # Define combining rules
                             "sigma = 0.5*(sigma1+sigma2); epsilon = sqrt(epsilon1*epsilon2); lambda = lambda1*lambda2")

        # Create a customised sterics force
        custom_sterics = openmm.CustomNonbondedForce(energy_expression)
        # Add necessary particle parameters
        custom_sterics.addPerParticleParameter("sigma")
        custom_sterics.addPerParticleParameter("epsilon")
        custom_sterics.addPerParticleParameter("lambda")
        # Assume that the system is periodic (for now)
        custom_sterics.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffPeriodic)
        # Transfer properties from the original force
        custom_sterics.setUseSwitchingFunction(self.nonbonded_force.getUseSwitchingFunction())
        custom_sterics.setCutoffDistance(self.nonbonded_force.getCutoffDistance())
        custom_sterics.setSwitchingDistance(self.nonbonded_force.getSwitchingDistance())
        self.nonbonded_force.setUseDispersionCorrection(False)
        custom_sterics.setUseLongRangeCorrection(self.nonbonded_force.getUseDispersionCorrection())
        # Set softcore parameters
        custom_sterics.addGlobalParameter('soft_alpha', 0.5)
        custom_sterics.addGlobalParameter('soft_a', 1)
        custom_sterics.addGlobalParameter('soft_b', 1)
        custom_sterics.addGlobalParameter('soft_c', 6)

        # Get a list of all water and non-water atom IDs
        water_atom_ids = []
        for resid, residue in enumerate(self.topology.residues()):
            if resid in self.water_resids:
                for atom in residue.atoms():
                    water_atom_ids.append(atom.index)

        # Copy all steric interactions into the custom force, and remove them from the original force
        for atom_idx in range(self.nonbonded_force.getNumParticles()):
            # Get atom parameters
            [charge, sigma, epsilon] = self.nonbonded_force.getParticleParameters(atom_idx)

            # Make sure that sigma is not equal to zero
            if np.isclose(sigma._value, 0.0):
                sigma = 1.0 * unit.angstrom

            # Add particle to the custom force (with lambda=1 for now)
            custom_sterics.addParticle([sigma, epsilon, 1.0])

            # Disable steric interactions in the original force by setting epsilon=0 (keep the charges for PME purposes)
            self.nonbonded_force.setParticleParameters(atom_idx, charge, sigma, abs(0))

        # Copy over all exceptions into the new force as exclusions
        # Exceptions between non-water atoms will be excluded here, and handled by the NonbondedForce
        # If exceptions (other than ignored interactions) are found involving water atoms, we have a problem
        for exception_idx in range(self.nonbonded_force.getNumExceptions()):
            [i, j, chargeprod, sigma, epsilon] = self.nonbonded_force.getExceptionParameters(exception_idx)

            # If epsilon is greater than zero, this is a non-zero exception, which must be checked
            if epsilon > 0.0 * unit.kilojoule_per_mole:
                if i in water_atom_ids or j in water_atom_ids:
                    self.raiseError("Non-zero exception interaction found involving water atoms ({} & {}). grand is"
                                    " not currently able to support this".format(i, j))

            # Add this to the list of exclusions
            custom_sterics.addExclusion(i, j)

        # Add the custom force to the system
        self.system.addForce(custom_sterics)
        self.custom_nb_force = custom_sterics

        return None

Is there any obvious reasons why using grand with openmmtools alchemy is not a good idea? I have asked a question over at grand too. Any thoughts/suggestions are appreciated!