Questions about softcore potential for vdW and electrostatic interactions

Dear OpenFE team developers,

` def _nonbonded_custom(self, v2): """ Get a part of the nonbonded energy expression when there is no cutoff.

    Returns
    -------
    sterics_energy_expression : str
        The energy expression for U_sterics
    electrostatics_energy_expression : str
        The energy expression for electrostatics
    """
    # Soft-core Lennard-Jones
    if v2:
        sterics_energy_expression = "U_sterics = select(step(r - r_LJ), 4*epsilon*x*(x-1.0), U_sterics_quad);"
        sterics_energy_expression += f"U_sterics_quad = Force*(((r - r_LJ)^2)/2 - (r - r_LJ)) + U_sterics_cut;"
        sterics_energy_expression += f"U_sterics_cut = 4*epsilon*((sigma/r_LJ)^6)*(((sigma/r_LJ)^6) - 1.0);"
        sterics_energy_expression += f"Force = -4*epsilon*((-12*sigma^12)/(r_LJ^13) + (6*sigma^6)/(r_LJ^7));"
        sterics_energy_expression += f"x = (sigma/r)^6;"
        sterics_energy_expression += f"r_LJ = softcore_alpha*((26/7)*(sigma^6)*lambda_sterics_deprecated)^(1/6);"
        sterics_energy_expression += f"lambda_sterics_deprecated = new_interaction*(1.0 - lambda_sterics_insert) + old_interaction*lambda_sterics_delete;"
    else:
        sterics_energy_expression = "U_sterics = 4*epsilon*x*(x-1.0); x = (sigma/reff_sterics)^6;"

    return sterics_energy_expression

It seems that the defined U_sterics_quad here was not consistent with the Eq. 1.3 of S1 in Supporting Information of JCTC 8, (2012): 2373-2382.

By the way, actually the alchemical transformation designed here is two-step protocol (usually named stepwise protocol in RBFE calculations). If using softcore potential for electrostatic interactions as defined in the above mentioned JCTC paper, the one-step protocol (also named concerted protocol) can be implemented? If not, are there some obstacles on implementing it?

Thanks for your reply in advance.

Pengfei

OpenFreeEnergy / openfe

Questions about softcore potential for vdW and electrostatic interactions #886