I want to run simulations with Amber force fields that have inconsistent treatment of 1-4 interactions. The protein FF and GAFF have the usual fudgeLJ=0.5 and fudgeQQ=0.8333. However, GLYCAM has no scaling of 1-4 interactions. I think I have found a way to work around this, but would like to double-check that I understand things correctly.
To avoid the fudgeLJ and fudgeQQ applying to the sugar molecules, can I simply set nrexcl=2 and leave the [ pairs ] section empty? That way the 1-4 pairs are never excluded in the first place, so don’t have to be added back in (with scaling applied). Or is there something I’m missing?
Obviously this wouldn’t work if scaling was still required, but by a different amount (such as with Lipid21). Would the developers ever be willing to allow fudgeLJ and fudgeQQ to be set at the molecule level instead of the parameter level, similar to nrexcl? That way scaling for such molecules could be set to 0 with pair interactions explicitly defined instead. One of the great things about GROMACS is its flexibility. Such a change would make it much easier to combine (supposedly compatible) FFs like those from Amber. Multiple scaling definitions could trigger warnings rather than errors to still discourage it being used inappropriately.
In fact it was quite trivial, you can try out set_pairs_fudgeas described briefly in the documentation, and see if that solves the problem with GLYCAM - I should probably add a few utility functions to make working with glycans easier.
However, I’ve seen posts in this forum from GROMACS developers saying that [ pairs_nb ] shouldn’t be used in this way, so I was looking for an alternative method.
Gromologist looks very interesting. Thanks for adding the new function. I’ll give it a try and do a few tests to check that the [ pairs ] section is behaving as you are suggesting under the hood.
I used Gromologist’s compare_topologies_by_energy function and checked that the implementation gives essentially the same LJ and Coulomb energies with default fudge factors, you can check whether that’s the case for your system:
I checked the energy of a single glucose molecule at step 0 with scaling factors of 1 vs the default Amber scaling factors with the explicit pair definitions using a) [ pairs_nb ] and b) [ pairs ] function type 2. Both gave me the exact same energies, so both methods appear to work well despite not being intended for this purpose (but for free energy sims).
Strangely enough, using nrexcl=2 and removing everything from [ pairs ] gave me slightly different energies. The LJ energy varied at the 3rd sig fig, whereas the Coulomb energy agreed until the 6th sig fig. So, I’ll be falling back on the free energy methods of modifying the pairs.
Thanks for your help milosz.wieczor. I will definitely check out Gromologist and figure out how to implement it into my work flow.
By the way, as I said I haven’t really worked much with GLYCAM, so if other elements of the workflow require custom scripts or weird workarounds, I’m always happy to implement them as extra functions, and perhaps make that into a tutorial.