GROMACS version: 2020
I am trying to run the Energy minmization of the Ribosome (PDB ID 6i0y) with a Protein in the exit channel. I have freezed the Ribosome atoms using freezegrps in mdp file and left the protein free.
But once I run the energy minimization I result in
very high potential energy in +7.65 *e+07 range which is extermely high and
did not reach the condition Fmax<1000.
Using this state I cannot run the Production MD simulation
I tried to removing peptide and running it still only with ribosome the potential energy is very high. May I know how to resolve this issue. Any input would be really helpful.
Some additional questions:
- Does freezing of atoms still influence in the potential Energy?
- Can I eliminate the forces of particular atoms which are in boundary, will it help in solving this issue.
Yes, frozen atoms still contribute to the energy - they are simply not moved during geometry optimization. If you’re freezing atoms that clash with other ones, your potential energy will stay high because the procedure won’t be able to resolve the clashes.
But usually Gromacs gives you extra information on where the highest forces occur, so the reasonable thing is to visualize the system and see what’s happening around there. With many automatic processing tools (and even deep learning models), you can get side or main chains going through aromatic rings or other crazy arrangements, so be sure to go back to the structure and resolve any major issues it might have.
But my focus is to freeze the Ribosome group atoms. Since just freezing still calculates the forces but to avoid the force calculation I was thinking to include energygrp-excl parameter this way it would not calculate the forces and result in high potential energy. Rather only the forces will be calculated from the atoms which are not freezed and not included in energygrp-excl But I am not able to set it in Gromacs 2023.
I tried checking with GROMACS 2019 as few posts suggested but results in issue of installation 403 Forbidden Regression_Test * http://gerrit.gromacs.org/download/regressiontests-2019.tar.gz
May I know is this right way and how to fix this issue of fixing specific atoms, not resulting high potential energy.
Don’t freeze them. Apply a very strong position restraint. I’m not sure that energygrp-excl works with the Verlet scheme. It needed the group scheme, which was removed (unless I am mistaken).
Hi @jalemkul I understand that Position restraints are used to limit the movement of atoms while still allowing them some freedom to vibrate. This is typically implemented by adding a harmonic potential to restrain the atoms around their initial positions.
But my goal is to fix the high potential energy and reduce time both, If I freeze the atoms fully and also eliminate the forces so that they dont result in high potential energy and also help me in not wasting my time calculating those forces.
Freezing does not skip the force calculation and therefore does not save time. It just skips the update of positions. That’s why you have unsolvable high forces. Again I emphasize the only legitimate way to do this in GROMACS is with position restraints.
Okay I get it. I will restrain it but along with that what if I add the if condition in the source code to not calculate the forces of certain atoms based on their atom id. I think this way I will be able to not only fix atoms but also reduce time.
In most systems, the majority of time is spent on calculating the interactions involving solvent anyway. If you want to focus on a small subsystem of the ribosome and save computational time, it might be more reasonable to extract a subset of the full complex and (optionally) restrain the interface with the rest of the ribosome, given it’s sufficiently far away. This has its own challenges and approximations, but for certain purposes the increase in sampling time might outweigh the error resulting from reduced mobility.
Thank you, Yeah I will try cropping the ribosome further to small system which i want to focus.
But, I see that restraining the subset of the atoms does not fully fix the positions there is still movement in the atoms I see during the simulation. Am I doing something wrong let me know:
I need to fix the ribosome atoms fully except for protein in exit channel, I tried using freezegrps but it was leading to very high potential energy, so my goal was to remove the forces on those ribosome atoms for which I am looking into source code still, which is not resolved.
As alternative I found position restraints will help me in fixing the position by applying the counter force so tried including in the topology.top file but did not work as mentioned below
#include "topol_Protein_chain_2.itp"
#ifdef POSRES
#include "posre_Protein_chain_2.itp"
#endif
#include "topol_Protein_chain_7.itp"
#include "topol_RNA_chain_A.itp"
#ifdef POSRES
#include "posre_RNA_chain_A.itp"
#endif
- I run the Energy Minimization on freezegrps led to high potential energy in e+14 and Fmax did not converge.
- so I tried with restraints on ribosome atoms but no restraints on protein chain 7 mentioned above. Potential energy resulted to e+07 which is still high, Fmax did converge.
But the Ribosome atoms displaced from the original position which i wanted to be fixed.
Continuing with the output of Energy Minimization then running nvt.
Lastly, I ran the Production MD. This displaced the ribosome atoms further with increased RMSD in ribosome atoms alone.
Am i doing something wrong?
- So position restraints does not fully fix the atoms position, they still jiggle and also displace from original position?
@jalemkul @milosz.wieczor
Restraints and constraints work very differently on the theoretical level. The former, as you say, just allow for some wiggle a tiny bit around the reference position (with the amount of wiggling controlled by the force constant(s). The latter require special treatment in the algorithms (integration, thermo- and barostats etc), as they remove some degrees of freedom.
Again, whenever you refer to high energy in minimization, I suggest you look at the highest-force atom indicated by mdrun, and check it visually. This will help you understand where the clashes are located, and whether they are acceptable given your strategy to put restraints.
If parts of your system drift away even when they were supposed to be held by position restraints, check the .log to make sure there’s an energy contribution coming from position restraints. If not, that might indicate you failed to turn them on correctly (e.g. because of a misspelled or commented-out directive in .mdp).
Dear Milosz,
I faced a similar error but for zeolite. As I run minimization the mdrun specifies the number of atoms with the highest potential energy. I checked them on the structure but I saw no specific disorder. What should I do about that? I attached the file showing atoms with high energy in Blue color. All of them are in the boundry.
It should be mentioned that all the structure consists of just Si and O in SIO4 form.
Thanks
The fact that they are on the interface suggests that maybe your topology does not feature bonds across the periodic boundary, or perhaps your box size was changed and does not exactly fit the size of the system. Hard to tell without further info but I’d first investigate these 2 options (in this order).
Thanks for the reply.
What do you mean by " your topology does not feature bonds across the periodic boundary"? What should I do to fix that? I also checked with larger box sizes but it still shows the errors. Also, there are lots of problems with NPT and NVT levels (segmentation fault, core dumped)
Thanks again
Well, first of all what is your idea about the periodic boundary treatment? Should it be a 3D/2D/1D crystal, or a single slab with extra vacuum around it?
Dear Milosz
The system is in 3D. I would like to study the diffusion of water in zeolite. in the first step, I must equilibrate the zeolite framework in a pure state (without water). The system is really simple. Also, I got the forcefield parameters from an article published in Nature. However, I face lots of strange errors. The issue that the instabilities are in the interface of the box is more strange to me.
Thanks again for your attention
The parameters are likely fine, what I’m asking is how you generate the topology. In the framework, neighboring atoms are connected by chemical bonds, and these are listed in the topology. I’m simply asking if the atoms with, say, the highest values of the Y-coordinate are covalently linked to atoms with the lowest values of the Y-coordinate through the periodic boundary (and similarly for X and Z).
I assume you understand how PBC is applied in molecular systems, but you can go to the VMD > Graphical Representations > Periodic tab and turn on neighboring periodic images to understand my point.
If the atoms that are neighbors through the periodic boundary don’t have covalent bonds, they will create clashes as they’d only interact through nonbonded interactions. Then, during minimization box size remains fixed, so there’s little the algorithm can do to relieve the strain.
Dear Milosz,
Thank you for explaining the concept of Periodic Boundary Conditions (PBC) to me. In the picture I uploaded, you can see that if the oxygen atom (red) at the interface connects to the silicon atom (yellow), a network is formed (SiO4 network). Because the structure of SiO3 in the interface is incomplete. However, I’m unsure how this connection is defined under PBC. I also suspect that the “inconsistent shift” errors might be related to this disorder in the structure. I believe I need to include a specific setting in the .mdp file that is tailored for zeolite structures.
What do you suggest?
The thing is, VMD will not visualize these connections anyway, but if you want to model a continuous solid state without any breaking points or irregularities, these bonds must be included in your topology (.itp/.top) files. You can check the atom numbers in VMD and verify if they are present in the topology file under the [ bonds ] section.
I have never really worked on material science so I don’t know which tools are available to handle it properly, but if I remember correctly, Gromacs’ gmx pdb2gmx cannot currently add bonds through periodic boundaries.
I could add a function in my Gromologist tool to automatically detect atoms that are close in the periodic space but not in “real” space, and add a bond in the topology between them - that would actually be useful for e.g. handling periodic DNA systems; let me know if that’s of interest, and perhaps send a sample topology file to test/validate.
You’re also right that there’s an .mdp setting to enable the handling of periodic molecules in the simulation, you should set periodic-molecules = yes.
Thanks again for your help,
You are right. I think pdb2gmx is not a good tool for such complex structures. I also checked periodic-molecules = yes before but it didn’t work.
I have worked on polymer systems for years and this is the first time I’m working on Zeolites.
I think I have to use other tools to make topology. That’s very kind of you if you suggest any other tools to check the peridic structure.
Thanks again
As I said, if you want me to take a look and try to implement it, feel free to send me the (even most basic) structure + topology to milosz.wieczor@irbbarcelona.org.
Gromologist has subroutines for adding bonds in the topology, the only question is how to easily identify pairs of atoms that should be bonded across the PBC boundary. I could add a criterion saying “pairs between type A and type B whose PBC-distance is less than cutoff d and regular distance is more than half of the box size”.
Still, I don’t know very well what your topology looks like, so maybe there’s a more general way to get there with the basic steps.
