Molecular Dynamics Simulation - Science topic

I add some data about my research

Shuo-Lin Weng

Hi! No need to worry about repeating. In Amber, you can periodically save restart files by using the ntwr flag. Ever n steps, a restart file will be created which you can then use to continue your simulation.

Amber manual

Every ntwr steps during dynamics, the “restrt” file will be written, ensuring that recovery from a crash will not be so painful. No matter what the value of ntwr, a restrt file will be written at the end of the run, i.e., after nstlim steps (for dynamics) or maxcyc steps (for minimization). If ntwr < 0, a unique copy of the file, “restrt_”, is written every abs(ntwr) steps. This option is useful if for example one wants to run free energy perturbations from multiple starting points or save a series of restrt files for minimization. Default = nstlim.

To avoid saving water, you need to change the ntwprt flag.

The number of atoms to include in trajectory files (mdcrd and mdvel). This flag can be used to decrease the size of the these files, by including only the first part of the system, which is usually of greater interest (for instance, one might include only the solute and not the solvent). If ntwprt = 0, all atoms will be included. = 0 (default) Include all atoms of the system when writing trajectories. > 0 Include only atoms 1 to ntwprt when writing trajectories.

You may find this helpful:

How did you decide how many ions to pack in a simulation box ?

Yes and no. There are minor differences which can cause your results to be false. Reducing the grid box would definitely decrease the total volume so that the interactions will increase just like in a more concentrated situation. However, when your ligand and/or receptor (target molecule) goes outside the walls of grid box, it comes back from the other side (Leaving the grid from right, entering from left etc...). So in that sense, you create a new collision probability for the atoms to interact. Let's assume you have a tiny drug that you want to bind to your large protein to inhibit/suppress its activity. In terms of theoretical chemistry, it may not have any chance to attack and residue on the right side of that protein. As you limit the grid-box, some right parts of that protein will enter from the left side and your drug will be able to interact that side as well which is not supposed to happen. So if you want to increase the concentration go ahead and do it. Please do not decrease the grid-box volume to satisfy that need since it will yield wrong outcomes.

Jobin Thomas this is the pdbqt file of the docked protein-ligand complex.

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 :)

If you don’t do it already I would take care that you calculate the free energy from a bunch of frames that represent different conformations of your structure. This will give you (hopefully) a more realistic estimation. good luck

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?

It depends on the algorithm to perform the time evolution, but usually this is a sign of too large of a time step. Energy conservation is not normally enforced in the numerical evolution. Remember further that DFT is a non linear problem so instability might be more subtle. As first test I would reduce the time step or check the time evolution algorithm.

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

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

Go through following articles:

1 )

2)

3)

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

The origin corresponds to the point where the X and Y coordinate axes are zero. The average structure is just that: The structure that results from averaging the XYZ Cartesian positions of your protein. The variance that explains the PCA analysis is based on how close or far the Cartesian positions of the protein are from the average structure.

Samta,

1- If you want to generate the AM1-bcc charge for the ligand, the Gaussian optimization could be skipped; check the ligand structure after adding the H-atoms, and make sure everything is correct.

2-If you want to generate RESP charges, use this keyword "opt hf/6-31g(d) pop=mk iop(6/33=2,6/42=6) scf=tight test" and optimize your structure using Gaussian.

compute f1 poly CH2 force/tally BOR

variable fx equal c_f1[1]

variable fy equal c_f1[2]

variable fz equal c_f1[3]

fix ff poly ave/correlate 1 100 100 c_f1[1] c_f1[2] c_f1[3] type auto file force.dat ave running

I written these lines to extract force but im getting error

ERROR: Unrecognized compute style 'CH2' (src/modify.cpp:1282)

Last command: compute f1 poly CH2 force/tally BOR

where I did mistake Can you please help me out.

Read the original paper ( ) which will give you an idea regarding the coarse-graining and its basis. Four heavy atoms to one bead in going from an all-atom model to coarse grain is the convention usually in Martini.

Thank you!....rpalma@tectijuana.mx

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.

Hi Piyali Mukherjee

Depending on the polymer, you could approximate the density to 1g/ml and with the volume of your box (6^3 nm^3) calculate the mass and then the amount of water molecules.

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!

i used it its pretty good thanks Leon Sulfierry

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.

@min zang

Were you using NAMD? For your MD

I remember that I solved this by using gromacs instead of NAMD and I've never met the same problem again.

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 usually use the Yasara-Structure application to carry out molecular dynamics simulations up to 100 ns, which only takes up to 6 days. Additionally, I use a PC with 16 cores

I opened up an account using the institutional email but it still brings the error

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.

useGalaxy

webgro

Thank you Ayaz Anwar .

and I calculated PCA using GROMACS software and visualized the 2D projection of PC1 and PC2.

By conventional methods the answer is no. But if you used a coupled method of MD and Machine Learning, the answer will be yes.

Add following to config file

CUDASOAintegrate on

Command to run:

namd3 +p1 +setcpuaffinity +idlepoll +devices 0 configfile.namd

You can try K-means for your clustering analysis with different K values and use the inflection-point method to identify the most suitable K values.

Kishore Kumar Sriramoju thank you very much for your advice.

you can use NAMD energy.

"... 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

Hello,

If in each of the replicas, you have different particle numbers (# of ligands) and you plan to use that variable as an exchange variable in the replica exchange algorithm, then I think it would not work in the ensemble that usually people use in MD (NPT).

You have to simulate your system in the grand canonical ensemble (uPT) and control the exchange probability of your ligands with the ligand reservoir to keep your chemical potential fixed.

I think this is not so simple and I don't remember a package that simulates such an ensemble in a simple/direct way.

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

The steps are explained in detail in the Amber Tutorial:

Jyri Kimari

i have attached this conversation link to the NAMD Mailing List with some details to make sure whether this method is usable or not.

As previously suggested, the PDB is a textual file. You can open it and make all the modifications you need. Be careful not to misalign the columns to avoid disrupting subsequent calculations.

If you have Gromacs installed, you can renumber the PDB file using the command: gmx editconf -f file.pdb -resnr starting_number -o output.pdb

Replace file.pdb with the name of your PDB file, starting_number with the desired starting number for residue numbering, and output.pdb with the desired output file name.

The error is related to the file of the complex. Have you tried remaking the prmtop file of the complex?

Dear Dr. Punet Kumar , You can use network pharmacological analysis to get your target protein.

I tried filling the form thrice. Still couldn't get the download link. Also can you give more details on which upper part because the option for "already filled the form can't be seen".

RMSD (root mean square deviation) represents the deviation of a coordinate set from a reference set. In molecular dynamics, the reference set is represented by the starting point of the simulations and so the calculation of RMSD allows us to observe how these coordinates evolve with respect to the starting point.

On the other hand, RMSF (root mean square fluctuation) does not indicate positional differences in the time, but rather measures the individual flexibility of residues, thus how much a particular residue moves during a simulation.

So basically in MD you plot RMSD vs time and RMSF vs residue number!

Just change the seed number in each MD.

This is a cutting processing modelled by LAMMPS.

In LAMMPS, you can use MOLECULE keyword to define water molecular and then use create_atoms to create given number of water molecular in the system.

Before using MOLECULE keyword, you need define a water moelcular TEMPLATE.

For DFT calculations you could use Quantum Expresso or CASTEP software, which are free to use (CASTEP is free for academic use).

CASTEP is also capable of running AIMD, but if you want only MD simulations i would suggest CHARMM or DL_POLY (https://www.scd.stfc.ac.uk/Pages/DL_POLY.aspx), which are also free for academic use.

Hope this helps!

Based on your steps, I think there is one more step you need. The cleaved and uncleaved Na-MMT needs NVT, NPT respectively. After that, adsoprtion can be done again with NVT for cleaved, NPT for uncleaved.

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:

Here I prepare to make a gromacs simulation on gold nanoparticles (I refer the webpage, http://sobereva.com/153) and DNA base cytosine , but when grompp, it prompts the error like this ,

==============================================

ERROR 1 [file topol.top, line 14]: Atoms in the .top are not numbered consecutively from 1 (rather, atomnr = 1, while at->nr = 3)

==================================================

I don't know how to solve. I want to consult google, but I couldn't open google in C-N. here I can only open bing but bing couldn't offer more help for the error. I don't know how to solve. Thank you !

Hello ,

how can I use NAMD for the simulation of Simulation for adsorption of CO2 molecules onto the MOF for methanol production?

There is the VMD Diffusion Coefficient plugin, which computes mean square displacement and self-diffusion coefficients.

But keep in mind that dynamical properties can be affected by the electron CP dynamics (depending on your fictitious electron mass) and that long simulation lengths are necessary to converge these calculations.

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.

Where did you install M3GNet?

For which software? LAMMPS or QATK?

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.

Alisa Gorislav Thank you !

Dear Soykan,

Thank you for the answer.

Could you please share the details of the setup if that is OK with you.

Thank you

The reason for atoms moving too fast even in 2fs as time steps may be due to two issues. First, the number of atoms in the file is so large that the simulation could not be processed. Second the distance between atoms is so short that reducing the time step to even 2fs is not enough for the system to calculate the required parameters.

You can try using coarse-grained or modifying the molecules' structure to better fit the simulation.

Principal Component Analysis (PCA) is indispensable in molecular dynamics (MD) simulations for several purposes, such as simplifying, interpreting, and extracting crucial information from the complex structural datasets generated during simulation.

MD simulations produce high-dimensional data, detailing the positions and movements of numerous atoms. PCA excels at identifying the essential motions within a system.

Proteins and other macromolecules undergo critical motions for their function, such as domain rearrangements or ligand binding events. By focusing on the principal dimensions (typically PC1, PC2, and PC3), researchers can pinpoint the dominant conformational changes, thereby enhancing the understanding of the system's dynamics.

Furthermore, PCA contributes to the visualization of conformational changes. Throughout simulations, biomolecules transition between different structural states, and PCA facilitates projecting these changes onto a lower-dimensional space.

This not only simplifies the visualization of complex motions but also enables the identification and interpretation of distinct conformational states.

You can just add the following to your script and change the first time step accordingly.

set inputname output.restart [put output name corresponding to your output]

bincoordinates output.restart.coor

binvelocities output.restart.vel

extendedSystem output.restart.xsc

firsttimestep x [put in the last timestep you obtained from your logfile]

numsteps y [put the last timestep you hope to achieve]

Make sure the firsttimestep in your initial script which read 0 is hashed out so that it utilizes the new firsttimestep you added from your previous run.

Also hash out the temperature in the restart script (not the set temperature). It would look something like this depending on how you formulated your simulation script.

# -------Temperature control----

set temperature z [write temperature set by you in the previous run]

#temperature $temperature; # initial temperature

Hash out any minimization or reinitvels that you might have added in the initial script.

This is all you need to restart the simulation from the last timestep recorded. Don't forget to record the logs and the outputs for the same under different names. Otherwise, your previous results might get overwritten. Hope this helps!

thank you @tanmaykumar varma, i added '#' symbol and it worked!

Okay! I'll do this.

By free version of Autodoc we can't do MDS

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

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

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.

Heme parameters are already present in charmm36. You can prepare protein bound to heme easily in gromacs.

Dear Amol,

See the following link:

Ajay Yadav i tried to modify RNA to DNA by using discovery studio by selecting deoxyribose. However when i select Show sequence, I can still see Uracil instead of thymine. Nikhil Maroli Ajay Yadav Do you have any suggestions how to convert structure obtained by RNA composer to DNA aptamer.