GROMACS version: gromacs/2024.3
GROMACS modification: No
Hello everyone,
I’m currently trying to build a monolayer/protein system in CHARMM-GUI using the CHARMM36 force field and simulate it with the OPC water model in GROMACS. However, CHARMM-GUI does not allow me to select the water model after choosing CHARMM36, and it defaults to TIP3P. I’ve tried to manually rename the water molecules in the output files to match the OPC water model, but I keep encountering errors when trying to prepare the simulation in GROMACS.
Here’s what I’ve done:
Renamed the TIP3P water molecules to OPC (e.g., TIP3 to OPC, OH2 to OW, H1 to HW1, H2 to HW2). I also tried renaming TIP3 to SOL instead of OPC and it didn’t work.
Used the formula a=b=0.5×const=0.147722363a = b = 0.5 \times \text{const} = 0.147722363a=b=0.5×const=0.147722363 for calculating the dummy site (which I named MW or I also tried naming it EPW) as I’ve seen both used in different OPC.itps available online.
Attempted to update the input files accordingly but received the following error:
Fatal error:
Something is wrong in the coordinate formatting of file modified_step5_input.gro.
Note that gro is fixed format (see the manual).
I’ve also tried modifying the PDB file to match the OPC water model, but I get similar errors when I try to process it in GROMACS.
My question is:
Is there a standard workflow for converting CHARMM-GUI outputs (with TIP3P) to the OPC water model for GROMACS simulations using the CHARMM36 force field? If manually modifying the files created by CHARMM-GUI is a viable solution, are there any potential issues I should be aware of that might arise later in the simulation process?
Any advice or suggestions would be greatly appreciated!
Why do you want to simulate this water model? There is a reason why CHARMM-GUI uses TIP3 water for CHARMM36m.
The easiest solution is to strip the system of all the water molecules and resolvate (gmx solvate) with the correct 4-sites model, and then point it to the OPC topology. Renaming in the .gro file is okay as long as you respect the spacing and you point to the correct topology, but adding by hand the fourth site is way more complicated than simply re-solvating.
Hi, I want to use OPC as the water model with the CHARMM36 as the literature suggests this is the optimal combination accurate description of lipid behavior at the water–air interface. Would removing the water molecules and resolvating cause an issue with the Monolayer Builder system set up from CHARMM GUI as this creates 2 monolayers separated by a water slab? I included a screenshot of how CHARMM GUI builds the system. If I resolvated, would there be a way to ensure the waters are only added between the 2 monolayers like how it is built in CHARMM GUI?
Yes, there are several ways to achieve that. The easiest is probably to just strip the system from the TIP3 water, solvate the whole system again (so also the ‘air’ interface on the top and bottom of the system), and then remove the water molecules inside the lipids and in the air volumes by creating a group with gmx select containing all the water molecules that have a z position > than the bottom leaflet lipid heads positions and < than the top leaflet lipid heads positions.
Hi, thank you so much. I’ve tried this, but I first needed to create a .gro with water in OPC format (which I haven’t been able to find available online) so I made a .gro of a single OPC water molecule by using tleap in amber to get a pdb of a water box in OPC then converting it to gro format. Is there a better way to do this?
With my OPC.gro containing 1 OPC water molecule, I tried to run: gmx_mpi solvate
This correctly added waters to the entire system (including the ‘air’ interface on the top and bottom) but it is not nearly enough waters. I’ve tried specifying the box dimensions to be larger to get more waters to be included but this hasn’t worked.
If I know I need ~13755 waters to end up in between the monolayers, is there a way to specify this when trying to solvate my system?
First, the specific water model is defined by the topology you are using rather than the .gro file. This means that the most important part is using a four-points water model to solvate the system (as OPC is four-points). You can do this by specifying the flag -cs tip4p.gro in the solvate command. You then have to check that the atom names and numbering reported in your OPC topology (the .itp file) is consistent with what GROMACS added.
I would retry with this new flag and see if you get enough water everywhere. Then you will need to remove the water from where you do not need it (which will be a combination of gmx select and gmx trjconv most likely).
I’ve been using -cs tip4p.gro with the OPC topology and confirmed that the atom names and ordering in my OPC.itp file match the water added by GROMACS. However, even when I do this, the added water is not aligned well with my system — parts of the lipid monolayers (especially at corners or edges) remain unsolvated.
Here’s what I’ve tried to resolve this:
Centered both the no-water system and the water box
Specify the dimensions when implementing tip4p.gro file with matching dimensions to the no water system or larger box dimensions
Tried omitting and specifying the -box option during solvation
Rotated both the system and water box
Despite this, there are regions where the water does not fully cover the lipids. Once I get the system solvated correctly I plan to proceed with the approach mentioned to remove the excess water that is not between the monolayers. Are there any recommendations for how to avoid this issue?
If you are adding the water molecules to the box with the lipids, then GROMACS should respect the filling of that box, i.e., you shouldn’t get any molecule that doesn’t fit in the same box as the lipids. This means that you did something like gmx solvate -cp lipid_bilayer.gro -cs tip4p.gro -o [...].
I suspect what you are having is or i) a problem with the bilayer or ii) a problem with visualization. The first case might be that the bilayer is rotated with respect to the box, but I can’t see this happen easily and would kill instantly your simulation if the lipids are rewrapped over themselves across PBCs. I don’t think this is the case.
Case ii) is simply that the representation of the lipids is not the same as that of the water, for example you have the lipids rewrapped as in an hexagonal box box the water is being added in another way. How did you generate the lipids? CHARMM-GUI? what type of box are you using (symmetry-wise)?
Thanks for the detailed response — really helpful.
I did use gmx solvate -cp lipid_bilayer.gro -cs tip4p.gro -o [...] to solvate the system. I generated the monolayer system using CHARMM-GUI’s Monolayer Builder (then removed the default TIP3P waters so I could solvate with OPC). CHARMM-GUI built the system with a rectangular box.
If this is a visualization issue, is there a fix I should apply, or will it affect the simulation? When I tried minimizing, I got this error:
Energy minimization has stopped because the force on at least one atom is not finite. This usually means atoms are overlapping. Modify the input coordinates to remove atom overlap or use soft-core potentials with the free energy code to avoid infinite forces. You could also be lucky that switching to double precision is sufficient to obtain finite forces.
Could this be related to an inconsistent setup between the lipids and water, maybe due to wrapping or misalignment during solvation?
I do not think this is the water but rather the lipids. I think in some way you rotated the box and now when the PBCs rewrap the lipids around the box you end up with lipids overlapping with other lipids and this gives rise to infinite forces.
My suggestion would be:
i) Visualize the box with VMD and draw the box vectors (pbc box on VMD terminal); how do they look? If they fit the water but not the lipids, then the lipid configurations is wrong
ii) Take again the lipids from the original CHARMM-GUI output and check that the box is alright (again the PBC box visualized in VMD is a good way to check it). If yes, then strip the water molecules again and resolvate; if not, regenerate the box as something went very very wrong with the box generation.
