One thing I noticed is that you have 600 benzene molecules in the topology file, but the .gro file only has one benzene molecule defined. Did you create a new .gro file with more than one benzene?

I can't edit my old comment now, but after digging up the paper, I can see that there are some differences in the cutoffs and other options in your md.mdp file. The paper used a cutoff of 1.5, you used 1.3 nm, the paper used Nose-Hoover thermostat, you used modified Berendsen, and there are some other minute differences. Would it be possible for you to post a link to the actual .gro file that you used for the NVT or NPT steps, so that everything is clear?

right, yes, I can do that. The thing I posted was my initial gro file, which I editconf-ed, then I applied insert-molecule to get the 600 molecules.

initial 600 molecule gro file: ctxt.io/2/AACgtax-FA, 600 molecule gro file after NVT: ctxt.io/2/AACgtax-Eg, 600 molecule gro file after NPT: ctxt.io/2/AACgpV-XEg

I have also added my npt and nvt equilibration mdp files, if they can be of any help

@ShoubhikRMaiti, would it make a difference if I used a different thermostat? Also, in GROMACS 2020.5, we cant use rvdw != rcoulomb anymore...

From the NVT .gro file that you have supplied, the starting cubic box has 5.0 nm length. However, in the paper, they started with a 89.0511 nm^3 box, which corresponds to a length of 4.465 nm. Aside from that, they constrained all bonds with LINCS, not just the h-bonds as you have done. You could try changing these to see what you get. Even if after all this, the result does not match, then it could well be due to rvdw != rcoulomb, because the simulation is often quite sensitive to those parameters. Then you may have to use an older version of GROMACS to exactly reproduce the data.

About the thermostat, I am not sure it would cause that much of a difference, but I would use what they used in the paper, just to be sure. Because your goal here is to reproduce the data.

@ShoubhikRMaiti, I agree about the 5nm. I simply did that to ensure that all the molecules fit in the box, since i was going to equilibrate anyway with NPT later. But I will change the bond-constraints

@megamence Yes, the barostat should ideally bring the unit cell dimensions to the correct value. I am not sure why it didn't. A couple of things that you could try—(1) set both rvdw and rcoulomb to 1.5 instead of 1.3, (2) Do the NPT equilibiration step with Berendsen barostat and then use Parinello-Rahman in the MD step.

@megamence Also, what seems odd to me is that GROMACS 2020 does not allow rvdw != rcoulomb. I don't know everything about gromacs, but one thing I have seen is that some option combinations don't work. For example, if you use Nose-Hoover with Parinello-Rahman, you have to use md-vv integrator. Maybe changing the integrator (or something else) will allow you to set two different values to rvdw and rcoulomb. You will have to search the manual for this I guess.

Another thing, there is a bug in GROMACS when Nose-Hoover and Parinello-Rahman are used together which was fixed only in the 2021 version. So you could try the newer version.

@megamence I know that I am saying a lot of things at the same time here. But to be honest, I have no idea why you are getting the error, so I am throwing out my best guesses in the hope that one of them works.

Thanks for your suggestions @ShoubhikRMaiti! I really appreciate your time. So I am doing Nose-Hoover+Parrinello-Rahman, but for some reason, I am getting a ton of LINCS errors...

However, after doing the all-bonds constraint on my molecules, my density has risen from 823 to 862! thanks for catching that!

Hi @Xivi76

@megamence What are the LINCS errors that you are getting?