Gromacs - Science topic

if you are talking about individual 3 sets of the same system. The energy should be the same. But in the case of REMD, it will not be.

Besides standard analyses like RMSD, RMSF, and PCA, other useful analyses for molecular docking and dynamics include decomposing binding free energy, analyzing hydrogen bonds, studying protein motion via PCA, investigating hydrophobic interactions, examining water dynamics, constructing residue interaction networks, identifying allosteric sites, analyzing binding pathways, assessing solvent accessible surface area, and mapping electrostatic potentials. These analyses deepen understanding of molecular interactions and enhance interpretation of simulation results.

You may find this helpful:

Thanks for your response Hamed Zahraee , Is the on the fly approach you mentioned can be used to constrain the polymer surface area at different intervals

I would recommend trying out different analysis techniques to figure out which one gives you the appropriate time and spatial resolution. For my part I often feel like raw RMS values at too coarse grained and should be accompanied by structural views of representative conformations. Really depends on the thing you like to get across though.

One thing to keep in mind is that a „calculated average“ structure is typically not something that you would encounter in nature. it might not even be a conformation that’s working with your force field. If you want to do docking with a calculated average structure I would find the frame that holds the structure that is closest to this calculated average. good luck :)

Dear Sree Haryini ,

Here is suggestion from me:

1. Reduce the number of frames: As you suggested, using a smaller portion of the trajectory can significantly reduce computation time. However, make sure that the selected portion is representative of the system's behavior.

2. Increase the interval: Instead of analyzing every frame, you can increase the interval between frames. For example, you can analyze every 10th or 20th frame instead of every 5th frame. This will reduce the number of calculations needed.

3. Parallelize the computation: Calculate by using multiple core, you can parallelize the computation to speed up the analysis. Many MMPBSA tools support parallel computing. If you use the tool by Valdes you can add mpirun -np X and X can be replace with number of processors. Example mpirun -np 16. But you also need to consider the RAM available.

4. Calculate MMGBSA - As suggested by Martin Rosellen calculate GBSA is indeed will reduce the workload for calculation. I have try before that calculate the GBSA take 2 hours while PBSA take almost 1 days in the same trajectory.

Which FF you are using?

Nowadays, Forcefields inherently supports the phosphorylation of various reisudes, you donot need to add them.

You can read further on this one here:

For adding your desired residue, you need to follow steps outlined here:

Just a quick fix, instead of using pTHR, just name it to THX or some 3 digit letter, alongwith follow the steps of how to add residue in a forece field.

Histidin can be either delta or epsilon protonated, and it is user's responsibility to check the protonation state. Before runing the gmx pdb2gmx, you should have set this protonation state and change the residue name to HIE (epsilon protonated) or HID (delta protonated). Although, pdb2gmx can approximatly guess the protonation state which are fairley good.

Since, you ignore the Hydrogen atoms, pdb2gmx will rebuilt all those h-atoms based on the residue entery of the forcefield.

OXT seems to be the terminal atom. Is Val961 is the terminal residue? Do you have proper TER record in your pdb?

You can visualize final .gro file that you generated as output in VMD and then by clicking load data into molecule you have to select .xtc file and thus you can visualize trajectory for simulation run time.

It seems that the gmx_mpi -merge command does not exist in gromacs. You must have missed something before "-merge"

Did you find the solution Divya Jhinjharia ?

And we also require the force field for metalloproteins that is compatible with small molecules.

thank you so much

Hi Priyamvada P.,

To calculate the percentage of a specific secondary structure, such as an α-helix, through counting the number of occurrences of that structure and dividing it by the total number of secondary structures created, you can follow these steps:

1. Open the file "scount.xvg" using a text editor such as Notepad++.

2. Look for the section that displays the secondary structures.

3. Count the number of occurrences of the specific structure you are interested in, such as α-helix.

4. Count the total number of secondary structures.

5. Divide the number of occurrences of the specific structure by the total number of secondary structures.

6. Multiply the result by 100 to obtain the percentage.

Organic inorganic hybrid zeolite structures with methylene groups can be simulated using GROMACS and appropriate force field parameters for bonded and non-bonded interactions. MFI-type zeolite-specific force fields can be used to simulate adsorption in other zeolite types, including FER-, TON-, MTW-, and DON-type zeolites. Morphology-dependent MFI-type intergrowth structures in zeolite crystals lead to distinct internal and external molecular diffusion barriers, with barriers originating from mismatches in structure and pore alignment.

The task for computer modeling of zeolite structures with embedded metal ions raises the problem of the choice of the calculation method (force-field approach, semi-empirical or non-empirical quantum chemical methods, or density functional theory), as well as the problem of model choice (cluster models, embedding cluster models or periodical calculations). The formation of cationic positions of monovalent metal ions in zeolites are determined by distribution of aluminum through zeolite lattice. It is generally considered that there are preferential lattice sites for Al localization exist in a high-silica MFI zeolite structure, but the list of favored sites and the energies of Al stabilization depend strongly on the computational approach.

Generally, The force field approach accurately describes interatomic interactions within silicas, aluminophosphates, and zeolites, and is a valuable tool for modeling chemically complex zeolites, there are two main factors in the destabilization of cation position and increase of Lewis acidity " sterically hindrance of the cation adsorption due to structural peculiarities of the site and the decrease of adsorption interaction which could be associated with distant mutual siting of two aluminum charge-compensating ions in the zeolite lattice”.

I am struggling with the same problem. Agnivo Gosai can you elaborate?

Many thanks Jd Tamucci for your response, very helpful indeed.

Hi, it is basically a hit and trial thing, where you test different lag times and see at what lag-time, the convergence is achieved and explains the kinetic of the system.

There is a discussion at this page:

and here is a tutorial, where they discuss how to build, have a look and hopefully you will get an idea of how to select lag time:

Hope it will be helpful.

You can compile as-

Hi,

You can download DMF coordinates from PubChem or Zinc database. PubChem provides SDF format files. After downloading, you can convert the SDF format to PDB format using OpenBabel or other source code. Then, you can generate an RTP file for DMF.

Check the GROMACS tutorial at http://www.mdtutorials.com/gmx/complex/02_topology.html.

Here they provide the procedure to generate an rtp parameter for a new molecule based on the force field.

Hope this will help you!

Ayaz Anwar make a single pdb file for every combination and try using protein-only option.

Jonathan Purnomo,

For gromacs users open-sourced tool easy to use specifically binding free energy calculations in a single run for each frame of the Simulation.

I got it, thanks

I appreciate your curiosity regarding the advantages of running Gromacs on Windows. While both Windows and Ubuntu are viable platforms, the choice often depends on user preference, familiarity, and specific requirements. Windows, being widely used, provides accessibility for a broader user base, especially those accustomed to its interface. Additionally, some researchers may find it convenient to integrate Gromacs seamlessly into their existing Windows workflows. It's essential to recognize that the suitability of an operating system for simulations can vary based on individual needs and project specifications.

I never tried it on GROMACS, but the protonation state can also be changed using H++ server, which is a step in Ambertools MD simulation.

Please follow this link: PCA and FEL through gromacs (dlcompbiobiotech.com)

Do you select the protein as the reference group for the calculation when you use that command? What do you mean with "not satisfactory"? With respect to what?

1. Create an index file with the specific segment as a separate group. You can use the make_ndx command for this. If you want to create a group for residues 1 to 20:

$echo -e "r 1-20\nq" | gmx make_ndx -f md.gro -o index.ndx

This command will create a new group with residues 1 to 10 and save the index file as index.ndx.

2. Calculate the RDF using the gmx rdf command and specify the new group as one of the groups for which to calculate the RDF:

$echo -e "Group1\nGroup2" | gmx rdf -s md.tpr -f md.xtc -n index.ndx -o rdf.xvg

To calculate the function fs(k,t) in GROMACS, you can follow these steps based on the information from the search results:

1. Velocity Distribution : The result is multiplied by the standard deviation of the velocity distribution "square root of (kT/mi)".

2. Total Energy : The resulting total energy is not an exact correspondence due to various factors.

3. Degrees of Freedom : The total number of degrees of freedom (Nf) and Boltzmann's constant (k) play a role in determining the coupling strength.

4. Mass Parameter : The coupling strength is influenced by the constant Q (mass parameter of the reservoir) in the Nosé-Hoover equations of motion.

5. Temperature Coupling : GROMACS primarily operates in a constant temperature (NVT) ensemble for most calculations.

6. Integrators : GROMACS uses leap-frog and velocity Verlet integrators for accurate integration with temperature and pressure coupling.

For a detailed calculation of fs(k,t) in GROMACS, it is recommended to refer to the official GROMACS documentation and possibly seek further guidance from forums or experts in the field for specific implementation details [1][2][3][4][5].

[2] https://manual.gromacs.org/current/reference-manual/algorithms/molecular-dynamics.html

[5] https://www.researchgate.net/post/How-to-calculate-Incoherent-Intermediate-Scattering-Function-with-dl-poly

I hope someone finds this useful, I have myself been struggling with the same problem. My solution to this is looking at the chemical time and ns/day data, which is constantly updated.

To check this go to MONITOR -> double click running MD job -> check the last entry (status -running).

Double click the last entry for a 100ns default Simulation 100000 chemical time means 100ns this might give a rough idea of percentage job progress.

Going by the ns/day data can give you a rough idea of simulation speed like 11ns/day for 100ns simulation will finish on 10th day of the run applied.

The output you get is mostly defined in the mdp file you use to compile the binary tpr file. You probably have set to zero these quantities: nstxout, nstvout, and nstfout.

Change those quantities in the mdp file (i.e. set them > 0) depending on if you want positions (nstxout), velocities (nstvout), and/or forces (nstfout). Be aware that the output is going to be very bulky, as this is a lossless format. Take also a look at this link under the output definitions here: https://manual.gromacs.org/documentation/current/user-guide/mdp-options.html#output-control

Take a look here: https://manual.gromacs.org/current/onlinehelp/gmx-msd.html

gmx make_ndx -f original.tpr -o index_res10-50.ndx

> r 10-50

> q

gmx grompp -f minimization.mdp -c equilibrated_system.gro -p topol.top -o res10-50.tpr

gmx trjconv -f production.xtc -s original.tpr -o production_res10-50.xtc -pbc mol -n index_res10-50.ndx

> r_10-50

gmx hbond -f production_res10-50.xtc -s res10-50.tpr ...

Hope this helps.

Hello RG Community in the above instance subsequent to the ffTK first tab blanks will be drawn for psf, par etc., accordingly working through all tabs depending will function best:).

I am using Amber in WSL with RTX 3050 mobile. I am getting above 50 ns/day speed for MD calculations. Applying a small tweak to the clock speed of the GPU via overclocking has helped me increase the calculation speed up to 56 ns/day.

"... Is it primarily electrostatic , or does it calculate the number of H-bonds based on the angle and distance criteria?..."

There is no contradiction between these two. Indeed, the nature of the H-Bonds in GROMACS is mainly electrostatic, the H-Bond is so to say “caused” by an interplay of electrostatic interactions between donor, acceptor and H atoms but also donor-acceptor LJ interactions. On the other hand, the algorithm “recognizes” the existence of a H-Bond applying geometric criteria (angles and distances).

I hope this helps

Paulo Netz

You need to specify -nojump as you convert the trajectory file

e.g.

"Option -pbc sets the type of periodic boundary condition treatment:

...

nojump checks if atoms jump across the box and then puts them back. This has the effect that all molecules will remain whole (provided they were whole in the initial conformation). Note that this ensures a continuous trajectory but molecules may diffuse out of the box. The starting configuration for this procedure is taken from the structure file, if one is supplied, otherwise it is the first frame.

..."

also

Looks like your protein drifted into the periodic boundary box, making part of it to jump to the other side of the box.

You need to specify NoJump before exporting the coordinates into the viewer program. Consult the manual of the specific MD analysis program you use for the exact synthax to do this. e.g. https://docs.mdanalysis.org/stable/documentation_pages/transformations/nojump.html

thank you Emmanuel Onah for your detailed answer. i used the following script for energy minimization. I will add the nstxout-compressed to see whether it produces the xtc file or not. see confly I have tried using .trr file but it contains only one frame. although my log file contains information for about 473 steps.

; minim.mdp - used as input into grompp to generate em.tpr ; Parameters describing what to do, when to stop and what to save integrator = steep ; Algorithm (steep = steepest descent minimization) emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm emstep = 0.01 ; Minimization step size nsteps = 50000 ; Maximum number of (minimization) steps to perform ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions nstlist = 1 ; Frequency to update the neighbor list and long range forces cutoff-scheme = Verlet ; Buffered neighbor searching ns_type = grid ; Method to determine neighbor list (simple, grid) coulombtype = PME ; Treatment of long range electrostatic interactions rcoulomb = 1.0 ; Short-range electrostatic cut-off rvdw = 1.0 ; Short-range Van der Waals cut-off pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions

I am facing the same issue as well, .gro files and .cpt files are not getting generated

Can you please help me with the structure factor calculation of ionic liquids using gromacs?

Thank you

Hello! From my experience, it could be a bug of Gromacs, which should have been resolved with version 2023.

Perhaps the error you encountered is due to some issue with the xtc or tpr file. Please check that both are correct; you might want to try fitting the trajectory and then rerun the command.

Does gmx covar work correctly with the same parameters?

Dear S

MD simulations are a complex and often computationally expensive technique. It is important to have a good understanding of the underlying principles and be familiar with the chosen software and force field before starting a simulation. Don't hesitate to consult with experts and use available resources to learn more about MD simulations and their applications in nanoscale systems.

Here are some additional resources that you may find helpful:

Nor Farrah Wahidah Ridzwan I am so very thankful for your time

The error message "Fatal error: Topology include file 'LIG.prm' not found" in GROMACS typically indicates that GROMACS is unable to locate the parameter (.prm) file for a ligand (LIG) during the addition of ions step. This error usually occurs when the topology file for the ligand is not specified or cannot be found in the working directory.

Here are some steps you can take to resolve this issue:

After checking and addressing these points, try running the GROMACS simulation again. If the issue persists, you may need to provide more specific information about your system setup, force field, and the contents of the topology and parameter files for further assistance.

Thank you all for your time and suggestions. Unfortunately non of the above mentioned recommendations resolved the issue. I have done some other modifications, but they also did not resolved the outcome.

Hi Mohammad,

I don't think it's possible to determine the bead types easily.

First look at the martini DB (https://mad.ibcp.fr/force_fields).

If you cannot find it there, you need to know if there is any experimental or All-atom simulation data regarding this molecule. There are some approaches, but at the end, you have to show your model works (through a series of tests including free energy validation)

1- I suggest you go through this paper first:

"Martini 3 Coarse-Grained Force Field: Small Molecules"

2- If you can break it apart and find Martini itp files for some building blocks of it, I suggest you start there. (consulting this link: https://github.com/ricalessandri/Martini3-small-molecules/blob/main/tutorials/building_block_table.pdf)

This is called "Parametrizing a new molecule based on known fragments" You can find more info here, in the first section:

3- If you have the all-atom simulation data, you can check the second section of the link above.

4- There are some other tools such as VOTCA (https://www.votca.org/csg/input_files.html) But I have not used them.

Please tell us how it went.

Hi Muhammad,

I think Nanomodeler 2.0 (https://www.nanomodeler.it/) in this link provides both top and gro files. According to their documentation, these are based on MARTINI FF.

Please let us know how it went

You should choose the pose with the highest score (-10.1) since it has a significant score difference with the second pose.

Docking algorithms like Autodock Vina evaluate multiple factors when calculating the docking score, including van der Waals interactions, electrostatic interactions, and hydrogen bonding. While the number of hydrogen bonds can be an important indicator of ligand binding, it's not the sole determinant of the docking score. Other factors such as the strength and geometry of the hydrogen bonds, as well as non-polar interactions, can also contribute to the overall score.

Sathyanarayan Balaji

1. Clean your pdb using charmm-gui (remove extra chains and water molecules)

2. Then try to run apo (protein only) simulation using gromacs

(Choose charmm27 ff)

if it goes well then try with ligand (generate topology for ligands from ligpargen or swiss param).

You can try

for easy execution

Look here:

Okay! I'll do this.

In the terminal, enter the production outputs directory and type:

vmd file.gro file.xtc

$ cd ~

sudo apt-get install libnotify4 xdg-utils libxcb-dri3-0 libgbm1 libnss3 libgtk-3-0 libatspi2.0-0 g++

sudo apt -y install ncurses-term

$ cd ~/Downloads

wget https://registrationcenter-download.intel.com/akdlm/IRC_NAS/163da6e4-56eb-4948-aba3-debcec61c064/l_BaseKit_p_2024.0.1.46_offline.sh

sudo sh ./l_BaseKit_p_2024.0.1.46_offline.sh

$ apt install ocl-icd-libopencl1

$ cd ~/Downloads

mkdir neo

cd neo

wget https://github.com/intel/intel-graphics-compiler/releases/download/igc-1.0.15136.4/intel-igc-core_1.0.15136.4_amd64.deb

wget https://github.com/intel/intel-graphics-compiler/releases/download/igc-1.0.15136.4/intel-igc-opencl_1.0.15136.4_amd64.deb

wget https://github.com/intel/compute-runtime/releases/download/23.35.27191.9/intel-level-zero-gpu-dbgsym_1.3.27191.9_amd64.ddeb

wget https://github.com/intel/compute-runtime/releases/download/23.35.27191.9/intel-level-zero-gpu_1.3.27191.9_amd64.deb

wget https://github.com/intel/compute-runtime/releases/download/23.35.27191.9/intel-opencl-icd-dbgsym_23.35.27191.9_amd64.ddeb

wget https://github.com/intel/compute-runtime/releases/download/23.35.27191.9/intel-opencl-icd_23.35.27191.9_amd64.deb

wget https://github.com/intel/compute-runtime/releases/download/23.35.27191.9/libigdgmm12_22.3.11.ci17747749_amd64.deb

$ wget https://github.com/intel/compute-runtime/releases/download/23.35.27191.9/ww35.sum

sha256sum -c ww35.sum

$ sudo dpkg -i *.deb

# cd ~/

sudo git clone --recurse-submodules https://github.com/codeplaysoftware/syclacademy.git

$ cd ~/Downloads

wget ftp://ftp.gromacs.org/gromacs/gromacs-2023.3.tar.gz

tar xfz gromacs-2023.3.tar.gz

cd gromacs-2023.3

mkdir build

cd build

$ cd ~/syclacademy/External/Catch2/single_include/

cp -r catch2 ~/Downloads/gromacs-2023.3/build

cd ~/Downloads/gromacs-2023.3/build

ls

cd ~/syclacademy/Code_Exercises/Exercise_01_Compiling_with_SYCL

cp -r source.cpp ~/Downloads/gromacs-2023.3/build

cd ~/Downloads/gromacs-2023.3/build

ls

# withou exit the build directory activate the environment of one api using the following comand:

$ source /opt/intel/oneapi/setvars.sh

icpx -fsycl -I ~/syclacademy/External/Catch2/single_include -o a.out source.cpp

ls

# this should create a file called a.out in the build directory

cd ~/

sudo apt install gcc cmake libpomp-dev hwloc libhwloc-dev

cd ~/Downloads/gromacs-2023.3/build

$ sudo su

$ source /opt/intel/oneapi/setvars.sh

$ cmake .. -DGMX_BUILD_OWN_FFTW=ON -DREGRESSIONTEST_DOWNLOAD=ON -DCMAKE_C_COMPILER=icx -DCMAKE_CXX_COMPILER=icpx -DGMX_GPU=SYCL -DGMX_HWLOC=ON -DCMAKE_INSTALL_PREFIX=/opt/gromacs-23

$ make -j 8

$ make check

$ make install

1. open the work directory

2. run as root

3. open the intel environment with: $ source /opt/intel/oneapi/setvars.sh

4. activate gromacs with: $ source /opt/gromacs-23/bin/GMXRC

Can you provide the code here so that we can tell you what parameter to change for push. In NAMD, changing the sign of velocity you can do the push or pull.

This means that during the simulation, your molecule drifted into the periodic boundary box. See here https://manual.gromacs.org/current/user-guide/terminology.html#suggested-workflow to learn how to fix this in gromacs

It is not clear what you want to do. But gmx genrestr can generate position restraints (the default value is 1000 kJ/mol nm^2, you can change this though). Moreover, pdb2gmx also generates restraint file by default. You can use it by #ifdef statement.

The solution is also provided in the error:

"Either there is a mistake in your chain, or it includes nonstandard residue

names that have not yet been added to the residuetypes.dat file in the GROMACS library directory."

Most of the forcefields have parameters for regular/standard (e.g. STY) modified amino acid residues, and they can recognize the modification. If your amino acid is not among them, you need to add the modified residue in the file "residuetypes.dat". You can open this file and see all the standard residues are written there. Create a backup copy, and then edit your desired residue.

your initial coordinates are creating this issue, you need to check for bad contacts/clashes and minimize the structure before doing NVT or NPT.

Hi Larwubah Kollie,

You can use the following link to download amber99sb-ildn force field parameters.

Citra Hasanah Please, send us a message here (https://groups.google.com/g/gmx_mmpbsa) with the error information and we'll help you through the installation process

Ankur Chaudhuri How do you solve the problem in VDM?

* the 2 would be the periodicity the 180 the Phase Shift

Sorry, I do not have any experience with STaGE or even Si simulation. I already wrote that it is a lot easier to work with already created/tested topology rather than creating your own. If you want to learn then you need to follow the manual of this software, that will be the place to start.

Plus, again, you can create your own topology and match it against the topology files from that paper. It will give you confidence that the procedure you followed is good.

AAC is actually the ligand that I'm studying

ACN stands for acetonitrile

I must mention that the .itp file was included in virtualchemistry.org

This works for me;

gmx mdrun -s md.tpr -cpi md.cpt -deffnm md

md is the file you were runnig. E.g if your tpr is AB.tpr it mean ypur cpt is AB.cpt; then type gmx mdrun -s AB.tpr -cpi AB.cpt -deffnm AB

Bibo Feng

Can you please provide the tutorial pdf as I am not able to access via the link you mentioned in the question? I am starting with REMD simulations and struggling to find some good documentation.

Thanks in advance.

thanks. but since the linked files on the tutorial are updated to v 5 and of v. 4 is no longer available I wonder if i can continue v 4 with v.5. any idea?

To convert x-axis from ps to ns

write following formula (IN BOLD) at

PLOT---> AXIS PROPERTIES---> (submenu) TICK LEBELS---> AXIS TRANFORM $t/1000

Yes, I got it

Thank you Noor Alomari ,

I got the dihedral angles using the VMD module after pressing 4. But i want to analyze the dihedral angles (PHI and PSI), throughout the simulation.

Please guide something.

Aashish Bhatt Thank you for your response, let me give a try

Hello Shradheya Raj Rajeshwar Gupta and Vivi Richard , I am also facing the same problem. The only solution I found so far is to use a job scheduler like SLURM, But, currently I am facing problem in correctly installing it. Have you find any solution overcome this without use of a Job scheduler?

Please suggest.

Thanks