How do i decouple multiple species lets say positive ions and negative ions during free energy calculations ?
For decoupling single species A, you can specify in : “couple-moltype = A”
But what should the parameters be if i want to decouple species A and B ?
I tried “couple-moltype = A B” , but that gave me an error.
I tried “couple-moltype = A " and “couple-moltype = B”, that also gave me an error.
The user manual instruction is not that clear on decoupling multiple species.
" If you want to decouple one of several copies of a molecule, you need to copy and rename the molecule definition in the topology” - Gromacs Manual
The code is not generally designed to handle multiple species being transformed. The only valid input into couple-moltype is a single [moleculetype] from the topology. If you want to transform multiple entities, they must be unified into one [moleculetype] and that name specified for couple-moltype. This may involve a lot of topology hacking, depending on what you’re doing.
This refers to the case in which there are many of one species and you only want to transform one molecule, e.g. the hydration free energy of a water molecule in a box of pure water. You can’t specify SOL as couple-moltype because that applies to everything; you would have to create a duplicate [moleculetype] with identical topology definition but a different name.
Ah I see, Currently I am trying to calculate free energy calcuations of ionic liquid(solute) in ionic liquids(solvent). Now, the solute ionic liquid has cation and anion species which has separate topologies. So the problem now is how do I treat these as a single species ? Can you refer me to any help page or documentation regarding that “topology hacking” you mentioned above.
Yes thats correct, thats how I do it, if it were a single species.
There’s no documentation because this isn’t something that’s normally done. But if you have A and B in separate topologies, you have [moleculetype] definitions for each. You need to merge them, which then requires renumbering all of the interactions involving the second species. So rather than having A with 1…NA and B numbered from 1…NB, you will have A+B where B has been added to A and renumbered from (1+NA)…(NA+NB). You must be very careful to renumber appropriately, otherwise the topology will make no sense. There are no GROMACS utilities to do this unless you create .rtp entries and generate a topology with pdb2gmx -merge.
I have a similar question: I want to decouple a trimeric coiled coil peptide, the topology for which was prepared with CHARMM-GUI. It is represented as 3 different instances in the topology and manually rewriting the numbers would take a couple of days or even more. I am trying to understand whether a subsequent decoupling of each of the three helices could be a correct approach? The system is a protein with restrained z-axis at the interface of a membrane and water.
Decoupling biomolecules is unlikely to ever converge or be numerically stable. I just posted this to another thread, but if you want a free energy of solvation of a protein (or complex), the typical approach is MM/PBSA, not a λ-dependent transformation.