Questions about softcore potential for vdW and electrostatic interactions

[freeenergylab]

trafficstars

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

[IAlibay]

Hi @freeenergylab - thanks for opening this. Apologies for the slow reply, we have mostly been out of office, we'll get back to you as soon as we can.

[freeenergylab]

Dose anyone help give me some comments about this question? Thanks!

[freeenergylab]

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.

It should be one typo error in the SI of this JCTC paper. The correct version which your code was based on should be like: . Is it right?

[freeenergylab]

@IAlibay @hannahbaumann @dotsdl @jameseastwood Hi guys,

For RBFE hybrid protocol, why not implement softcore potential for electrostatic interactions as defined in the above mentioned JCTC paper? If so, are there some obstacles on implementing it?

Looking forward to your answers.

Thanks.

Pengfei

[IAlibay]

Hi @freeenergylab - just letting you know what we're not ignoring you. I'm OOO today and will get back to you on softcore electrostatics early next week. The other question you have, might take a bit longer since I think there's some historical context from the Chodera lab.

[freeenergylab]

@IAlibay Thanks for your reply. Looking forward to your next comments.

[IAlibay]

Hi @freeenergylab

So for now I'll stick to answering your second question, i.e. "why have we not implemented electrostatic sofcore scaling". There's a lot of detail here so I'll offer a brief (and somewhat flawed) answer and we can possibly discuss this more (see my other message re: invite to the OpenFE slack).

NonbondedForce (which handles PME in OpenMM) allows for lambda scaling, but not altering the functional form. Thus to implement a softcored version one would either need to: 1) implement a new force as an OpenMM plugin at the C++ layer, or 2) create a combination of custom nonbonded forces which corrects the NonbondedForce scaling. The second option is somewhat easier, but has a performance impact and can be messy when dealing with reciprocal space contributions. The former is a rather large undertaking.

For now, the solution of not softcoring partial charges and just having them being annihilated through the NonbondedForce machinery is something we've considered as an imperfect but "suitable" temporary approach.

Long term, handling softcore PME is something we would love to do, particularly so that we can do LJ-PME for non-homogenous systems. However, it's not something we've been able to invest limited resources towards. We definitely would love it if folks in the community had spare capacity to deal with this!

[IAlibay]

For your first question about the softcore formalism, we need a bit of time to chat with some folks from the OpenFE TAC.

[freeenergylab]

Hi @freeenergylab

So for now I'll stick to answering your second question, i.e. "why have we not implemented electrostatic sofcore scaling". There's a lot of detail here so I'll offer a brief (and somewhat flawed) answer and we can possibly discuss this more (see my other message re: invite to the OpenFE slack).

NonbondedForce (which handles PME in OpenMM) allows for lambda scaling, but not altering the functional form. Thus to implement a softcored version one would either need to: 1) implement a new force as an OpenMM plugin at the C++ layer, or 2) create a combination of custom nonbonded forces which corrects the NonbondedForce scaling. The second option is somewhat easier, but has a performance impact and can be messy when dealing with reciprocal space contributions. The former is a rather large undertaking.

For now, the solution of not softcoring partial charges and just having them being annihilated through the NonbondedForce machinery is something we've considered as an imperfect but "suitable" temporary approach.

Long term, handling softcore PME is something we would love to do, particularly so that we can do LJ-PME for non-homogenous systems. However, it's not something we've been able to invest limited resources towards. We definitely would love it if folks in the community had spare capacity to deal with this!

I got your points. Thanks. Let me think a little bit.