Simulating a simple liquid with GROMACS

User discussions
1

GROMACS version: 2020.5
GROMACS modification: No

Dear fellow GROMACS users,

I’m new to this. I have very low experience with molecular dynamics simulations. My original expertise is COSMO-RS, a semi-explicit molecular thermodynamic theory dealing with contact statistics between molecules.

However, I am interested in running MD simulations for complex liquids so that I can track aggregation phenomena, through radial distributions functions, Kirkwood-Buff integrals and so on. Then perhaps later I can use these MD simulations as inspiration and comparison to extend COSMO-RS theory.

Since Gromacs is installed on my university server, I decided that I would try making a simple liquid cyclohexane phase equilibration (NVT, and then NPT) at ambient temperature and pression conditions. Then I could calculate some RDF, and somme self-solvation energy, using the equilibrated phase, as an exercise. I tried to get help from the various tutorials and manuals that are available online, all the “get started” or “quick and dirty” pages I could find.

But the problem is: creating an input file for such a simple system seems to be far from trivial and even with all of that help I couldn’t get what I want.

So here is how I tried to start, from scratch.

I opened my dear Avogadro software and drew/equilibrated a nice cyclohexane molecule, which I saved as a .pdb file, which seems to be a compatible format for GROMACS. Now, I learned that I need a topology file which contains the molecular topology, and also the force field parameters for my molecule.

So, I tried to input this cyclohexane.pdb file into pdb2gmx, but it didn’t work. It seems strange to me, but the simple cyclohexane molecule seems to not be found in the Gromacs database.

So after tweaking around, I eventually got to know that I can find this desired topology file on PRODRG. After asking for a token and inputting my pdb there, I got a cyclohexane.top file, apparently using some OPLS-AA force field which seems to be a relatively general, unbiased one. Nice. I am not interested in tweaking with force fields for the moment, my priority is to actually get operational, and able to run a calculation.

Now I wanted to create a starting box filled with cyclohexane molecules. I found some tutorial online and the following line got me what I want:

gmx_mpi insert-molecules -ci cyclohexane.pdb -nmol 1200 -box 5 5 5 -o chx_box.gro

Now I had my box of cyclohexane nicely stored in the chx_box.gro file. Very nice. I felt that I had some momentum there!

Then, I got to know that to equilibrate my box, I need to use the “mdrun” program. Upon reading the mdrun program, I learnt that my “chx_box.gro” is actually NOT the required input file to equilibrate. No! I need to specify the parameters of the simulation in a separate .mdp file, and I need to use the “grompp” program to create an input file that bundles the gro file, the top file, and the mdp file. Fine. I found a default one which seemed to be all-around good on the internet and saved it as example_mdp.mdp.

And then, I tried to actually bundle all of these together, with the following command line:
gmx_mpi grompp -f example_mdp.mdp -c chx_box.gro -p cyclohexane.top

But things went wrong here. I really did not expect that they would go wrong at that point. I thought that this was some trivial formatting and that I would get into trouble while running an actual simulation on the server.

The error message is:
Fatal error:
Syntax error - File cyclohexane.top, line 17
Last line read:
‘[ moleculetype ]’
Invalid order for directive moleculetype

I am out of options. Looking at the online error documentation of Gromacs didn’t help. Since I used the .top file for a single molecule, while the .gro file is for a system composed of that molecule, I guess that the problem comes from there. But here’s the thing: I have absolutely no idea about how to get a .top file for the entire initial box. How should I do this?

Best regards and thanks in advance for your help!

2

Hi,

the order of directives (e.g. #includes) is crucial in .top files; could you please show the topology file, so that I can help pin-pointing where the problem is?
Also, bear in mind that the topology file have to include the information for all the molecules in the simulation box, not just a single one.

The grompp call seems correct; I suggest you take a look at genconf: one easy way to obtain a large box of fluid is to stack smaller ones, which may be easier to produce via insert-molecules.

Cheers,
Michele

3

Dear MichelePellegrino,

Thanks for your response, and your help.

Enclosed you will find the topology file (on a weTransfer link as, because my account is new, this website’s interface does not allow me to upload files):

Do you spot anything wrong with it?

Best regards.

4

Alright, at a first glance it seems two things are missing, one at the beginning and one at the end:

  1. before [ moleculetype ] you need the [ atomtypes ] directive where you list all the atom types in the system;

  2. at the end of the file you need to include [ system ] (gives a name to the configuration you are simulating) and [ molecules ] which lists the number of molecules for each residure, e.g:
    [ molecules ]
    RE1 1234
    RE2 765

5

Hi Michele,

I don’t know how to fill the [ atomtypes ] directive. What should I put in here? Is there any example file? I didn’t find any on the Internet.

6

Hi,

there is one example in Gromacs 2021 reference manual under the section ‘5.6 Topologies’ (Reference Manual — GROMACS 2021 documentation).

Basically, you ned to provide name, atomic number, mass, charge, particle type (e.g. V: virtual or A:atom) as well as the sigma and epsilon paremeters of LJ potential. You can specify the type of combination law for vdW in the [ default ] directive before [ atomtypes ].

7

How do you calculate the Kirkwood-Buff integrals?