GROMACS version: 2020.4
GROMACS modification: No
Hello everyone,
is there an easy way to duplicate the entries in a topology file when adding a second identical molecule?
My system contains two identical polymer chains A and B with 337 atoms each. In my topology file, there’s entries for chain A and in the [ molecules ] section I’ve specified nmols = 2 to account for chain B. This works fine for simulations.
The issue arises when I run gmx energy on the system, I don’t get any interaction energies for chain B. I’m assuming this is because there are no entries in the topology file for its atoms (indices 338-674).
Any tips would be appreciated.
Cheers,
Ben
Hi,
According to what you wrote the topology file looks correctly defined.
Now the two molecules have the same name, the one you used in the directive [ molecules ], even if you had 2 different name in the coordinate file. If you differ mol name in the coordinate file, grompp will give a warning telling that the name in topology dominates the name in the coordinate files.
If you want to specify energy group in the mdp file, you have to generate an index file where you defined chain A and chain B and give such an index file as input to grompp. If this is the problem, you can rerun to extract the interaction energy.
Alessandra
Hi Alessandra,
thanks for responding.
Yes, the two molecules have the same name, so I’m not getting any warning from grompp. In fact, I haven’t had any issues running simulations or analysing them with Gromac’s inbuilt tools. I’ve just named them ‘chainA’ and ‘chainB’ here to explain the issue.
I’ve defined chainA and chainB (i.e. atoms 1-337 and 338-674 respectively) in an index file, specified the two groups in the .mdp run input file, and given these to grompp in preparation for gmx mdrun -rerun. The issue is that I get interaction energies for chainA (i.e. between chainA-solvent) but interaction energy = 0 for anything involving chainB. What I’ve assumed is that this is because atoms 338-674 aren’t specified explicitly in the topology file, but perhaps the issue is something else.
Ben
Hi,
that sounds strange. The topology description is used to describe both molecule. Can you upload the topol.top?
Maybe you define twice chainA? or most probably
what I have suggested (with rerun) is not working.
If you visualize the trajectory and if you do not see water molecules or chainA penetrating/overlapping into chainB, then the parameters are specify also for chainB.
Alessandra
Hi Alessandra,
unfortunately as a new user, I can’t upload files. I’ll include relevant snippets below:
From the index file (I’ve omitted all the atoms, as well as where the solvent etc has been defined)…
[ chain_A ]
1 2 ... 337 338
[ chain_B ]
339 340 ... 675 676
From the topology file, chainA begins with…
[ moleculetype ]
;name nrexcl
car 3
[ atoms ]
; nr type resi res atom cgnr charge mass ; qtot bond_type
1 Cg 1 G4X C4 1 0.207370 12.01000 ; qtot 0.207
and ends with…
338 Ho 16 DGX HO1 338 0.413269 1.00800 ; qtot -8.000
… and chainB is accounted for at the bottom of the topology file with…
[ molecules ]
; Compound nmols
car 2
From the run input .mdp file…
;Energy groups
energygrps = chain_A chain_B
And at the end of the log file from mdrun:
Epot (kJ/mol) Coul-SR LJ-SR
chain_A-chain_A -6.97257e+03 -5.38362e+01
chain_A-chain_B 0.00000e+00 0.00000e+00
chain_A-rest -7.57088e+03 -9.24845e+02
chain_B-chain_B 7.67618e+03 5.40154e+02
chain_B-rest 0.00000e+00 0.00000e+00
rest-rest -3.65684e+06 6.06919e+05
Visualising the trajectory, it’s clear chainB’s parameters have been specified as chainB behaves exactly like chainA (as expected). So the issue of no interaction energy being returned must either be an error in how I’ve set up and included the energygrps in the run input files, or something going wrong when the calculation of interaction energy is getting the parameters from the topology file.
Thanks for your help with this,
Ben
Hi,
your topology looks fine and from the information you provide I do not see problem in the setting.
It looks like a bug in printing out the energy. Did you try to run gmx energy on the re-run ener.edr? Did you get the same energies?
\Alessandra
Hi Alessandra,
OK, good to know all the files look fine.
I tried running gmx energy on both the re-run .edr file, and well as re-simulating the system with the energygroups specified and running gmx energy on the .edr output of that, and in both cases the interaction energy involving chainB was zero.
Thanks for your help with this.
Cheers,
Ben