We value your privacy
How can I do MMPBSA analysis on Gromacs trajectory using the MMPBSA tool of amber software packge?
Nov 28, 2012
Not sure it is possible, since these two codes use different output formats for their trajectories. An alternative, I think that you can convert the gromacs trajectory files (xtc or trr) into AMBER trajectory files with VMD. You will also need to convert your GROMACS topology file into the AMBER file format.
I would like to calculate auto-correlation function for tilt angle for the lipids in a bilayer from molecular dynamic simulation. I have the angles for all lipids from each frames for all subsequent frames but I need a code to calculate the auto-correlation function. Where can I find or get the code apart from writing my own?
How and with what program can I calculate the MD simulation of the proteins modified by polyGlu-tails of different lengths. Will suitable the Gromaks to solve such problem? Thanks for the help!
Nov 1, 2013
In addition to Hiqmet's response. It is also depend (a little) of the force field you want to use to model your poly Glu tails. Did you read the literature to see if somebody has do similar simulations?
I am doing MD simulation of ag-ab complexes using GROMACS and want to know the eletrostatic potential of a surface atom in a coordinate. Please tell me if there is any method to calculate the same that would be helpful to me.
Your question is not clear, but I think you mean how to calculate the electrostatic potential at a certain coordinate R originating from certain atom r_i. If this is the case, you must take the vector difference (R - r_i), then calculate the length of that vector, call that d_i. The electrostatic potential is (1/(4 Pi e0))*(q_i/d_i). Where e0 is the dielectric constant, and q_i is the point charge assigned to atom i by GROMACS. Then you must add over all atoms you want (sum over index i).
I tried to do MD refinement for a protein-membrane Bilayer complex by using Gromacs. In manual they mentioned only distance restraints between two atoms with in the molecule but here it is a complex with two molecules. Anyone could provide references to this topic is really appreciation.
One method is to rotate the bilayer so it is in the x-y plane . A harmonic restraint can then be applied in the z direction. From the Gromacs 4.6.1 manual, section 4.3.1
"Using three different force constants the position restraints can be turned on or off in each spatial dimension; this means that atoms can be harmonically restrained to a plane or a line. Position restraints are applied to a special ﬁxed list of atoms. Such a list is usually generated by the pdb2gmx program"
APBS calculate electrostatic energy by using the implicit solvent model. So, you do not need to include water molecules explicitly. APBS uses dielectric constant of water to calculate the electrostatic potential. Even in absence of any water molecules, obtained electrostatic potentials of solute is assumed to be similar to that of the solute present in the bulk water. If water molecule will be present explicitly, certainly it will affect the electrostatic potential, because APBS will consider water molecules as a solute similar to the protein/ligand.
Does the distance between two frozen atoms remain constant throughout the trajectory? If so, the atoms are actually frozen and just translated during the trajectory.
If this is not the case, there are a few ways to freeze atoms in g09 using cartesian coordinates. Although I have no experience in ADMP calculations in g09, you may try to see if any of these options work for your purpose:
"AddRedundant is synonymous with ModRedundant. This option may be used for job types other than optimizations."
2) try changing Geom=ModRedundant with Geom=ReadOptimize or Geom=ReadFreeze: you can find more on this keyword in the Manual (http://www.gaussian.com/g_tech/g_ur/k_geom.htm). In my experience, Geom=ReadFreeze works with ONIOM calculations, even though it's said to be deprecated in the Manual, while Geom=ReadOptimize doesn't work (at least in G09 rev A.02).
The AMBER force field can be used with GROMACS (and it is included in the standard distribution), so you do not need to use the AMBER package if all you need is the same force field. It is usually considered among the best force fields for protein molecular dynamics.
The GROMOS force fields are generally considered (from what I've heard from many, many users: I didn't compare them myself) to be not as good as AMBER, but they can be useful nonetheless.
What is critical however is the *parametrization* of your ligand. It has to be consistent with the force field you use and it has to be of good quality. Therefore your choice of force field will depend on how you can parametrize your ligand. Do you already have a known, good topology for the ligand? If yes, use the force field compatible with that topology. If not, look around to create an AMBER or GROMOS compatible ligand topology (I am unaware of tools able to create OPLS-friendly topologies but they may well exist).
Can anyone suggest a web based server or free software which is comparable to GastroPlus (GastroPlusTM, a physiologically-based simulation tool that predicts the absorption, pharmacokinetics, pharmacodynamics, and drug-drug interactions for drugs administered through intravenous, oral, ocular, or pulmonary routes)?
I have four proteins A, B, C and D. A is a transmembrane receptor and its the substrate that gets phosphorylated in its cytosolic domain. B is the kinase, therefore, Mg++ and ATP dependent. B also depends on C and D for its activity. There might as well be a domain swapping between B and C. Is it at all possible to build a reliable model of the ternary complex of BCD? I only have the sequence data. Direct structural information of these proteins is not available.
This would be a tall order. Like Felix said, you're best hope would be to find a model of a homologous multi-protein complex that you could use as a template for building a predicted model. If such a model does not exist, you could NOT perform de novo predictions of how a multi-protein complex assembles, even in the case where a reasonable model for the structure of each subunit exists.
A fruitful approach for modeling multi-protein complexes that cannot be studied directly by X-ray crystallography has been to perform docking studies in cryo-EM images of multiprotein complexes. The minimum criteria you need are high-resolution cyro-EM, and structural models for each subunit.
I've been told to use Xeno View to build my polymer and assign the Force Field parameters. The problem is that I want to use my own set of parameters and it seems that the program can only read parameters from a default list of frc file. Can you suggest a software that lets you build a polymer from a monomer, assign the atom types, add water and save the result in a topology file containing the parameters I need?
Commercial -- These tools allow you to either build or import individual molecules, assign force-fields, build complex molecular systems (initial condition) from the individuals, and then automate the running of simulations.
Of these, I have experience with the Scienomics software. It is very user friendly, has excellent molecular construction tools including a polymer builder and a crystal builder, great force-field support, and interfaces with a lot of different open source MD and Monte-Carlo programs, LAMMPS, GROMACS, etc.
I am just new to gromacs. Can anyone please help me in correcting these charges for the itp file, i.e. how to adjust them? These have been generated by prodrg, I have read the article for correcting it but couldn't make it. At further steps the error is generated while running the energy pressure and temperature.
May 2, 2013
You can use SWISSPARAM for generating topology for your ligands. But it uses charmm 27 forcefield. Its mentioned on the site how to use this topology with gromacs. Hope this helps you.
Energy minimztion probelm in Gromacs.
Apr 28, 2013
When I run the energy minimization step for Groamcs it gives the following error:
Double precision normally gives you higher accuracy.
You might need to increase your constraint accuracy, or turn
off constraints alltogether (set constraints = none in mdp file)
Steepest Descents converged to machine precision in 43 steps,
but did not reach the requested Fmax < 1000.
Potential Energy = -1.4074464e+06
Maximum force = 8.7960381e+05 on atom 2201
Norm of force = 3.9092959e+03
writing lowest energy coordinates.
How can I correct this?
First, you have to try to understand why that is happening.
1) Are you sure that the protein is well built (i.e. no bonded atoms inside aromatic rings, clashes, etc?
2) What are you minimizing? Protein, DNA? Small molecules?
Please carefully check the structure. Open the failed minimized GRO file and check for distortions or any information that can help you to evaluate the initial structure.
Besides steepest descent method, you can try minimizing the structure by conjugated gradient method (CG), a more thorough method. Simply change integrator to cg (inside minim.mdp):
integrator = cg
I am planning to make a molecular dynamics simulation by embedding my protein into nuclear membrane, to study interactions in active site. But I couldn't find any solution to put lipids different for nuclear membrane rather than cellular membrane. If possible I am interested in building system with Maestro to use Desmond.
Apr 24, 2013
I think the differences between cellular and nuclear membrane is only for their lipid compositions. So in the simulations point of view, you can use the same approach to embedded the MP in the membranes.
To construct your models of membrane and insert your membrane protein in it, you could use CHARMM-GUI (http://www.charmm-gui.org/)
Hi look at these workstation from Dell....
They may be costly but service is very good.......
May be you can also try IBM workstations
Gromacs run problem
Apr 17, 2013
I want to do 3ns simulations using gromacs on an i3 laptop. My problem is that I have to face the load shedding problem. After every 2 hours, the light goes for 1 hour, sometimes the light goes for more than an hour and a half and the laptop battery diminishes. I have heard that the graph will be effected if I shut the lid when only 10% battery is left.
Is there any solution for this problem? Can I split my simulations into different parts so that they are not effected by the light problem? If so, how can I do this? Or is there any other way to solve my problem?
yes, in Gromacs mdrun will save your run every 15 mins. please look into documentation for mdrun "-cpt" property of mdrun will help you to restart your simulation, when power goes off or PC crashes.
just append -cpt to your mdrun command
all the best
If you want to use GPU for simulation, then buy Nvidia Tesla. Cards like GTX or Quadro is nice for testing and debugging purpose. Simulations may run continuously for many days and weeks, which may affect the hardware. GTX or Quadro are not made for heavy continuous use, while simulation is long and computationally expensive process.
I would like to build a polymer membrane and do MD studies. Is there any free software available to build a polymer membrane and do MD simulation? In the future I want to study the interaction of polymer membranes with biomolecules. What is the procedure to build a model in gromacs?
You can calculate inner product of covariance matrix by g_anaeig tool using -inpr option. You can also calculate overlap in covariance matrix using same tool.
This article may be helpful.
I think you are a little bit confused. All MD production runs are unrestrained, i.e. the atoms spatial position is allowed to move.
However, in all MD simulations, the first steps (1-10 ns) must be restrained in order to allow a correct temperature (NVT) and pressure (NPT) equilibration. Only then, the restrains can be removed and the system is allowed to equilibrate on a MD production run.
Mount your GROMACS system, do an energy minimization EM) followed by a short NVT to equilibrate pressure (10-100 ps) and a restrained NPT (2-10 ns) to allow pressure equilibration. Only after this you can remove restrains and perform an unrestricted MD run.
Is it possible to perform a molecular dynamics simulation (GROMACS, AMBER ...) at a pH other than 7?
May 23, 2013
I want to do a molecular dynamics simulation with my protein. In vitro experiments were performed at pH 8. Is it possible to alter the parameters of the molecular dynamics properties with explicit solvent at pH=8? Or do I have to alter the protonation of my protein in a way that it mimics a pH of 8?
My experimental colleagues have no (little) difficulty in adjusting the pH of their solutions, as pH is a macroscopic parameter, like temperature and pressure. At the microscopic level of MD simulations, the macroscopic pH will be reflected in your choice of protonation states for titratable residues. It is generally not possible to provide much in the way of computational insights into pH-dependence of reactions. The pH is an ensemble average value and you are typically simulating systems with a copy number of one.
If your concern is about a particular residue in the active site, say a histidine, which is much less likely to be protonated at pH 8 than at pH 6.5, then you will need to run simulations with it singly protonated and others with it doubly protonated. Even at pH 8, where the probability of protonation is reduced by an order of magnitude from pH 7, remember that there are probably 10^15 copies of your protein in a small vial and millions of those histidines will remain protonated even at high pH.
I'm studying stability of certain enzymes with their mesophilic and thermophilic homologues using MD Simulation methods. Using the information obtained from the study I'm trying to design point mutations which would help in increasing stability. Is it possible to suggest some appropriate methods of studying such mutants quantitatively. e.g. in terms of melting temperatures or something like that?
I want to simulate DNA in water. All gromacs force fields do not support DNA. So I chose the amber force fields from Gromacs package.
I want to ask, are all the remaining steps the same for DNA simulations like the lyzome example?
I have done up to this step
pdb2gmx -ff amber99sb -f DNAAmber.pdb -o DNA2.pdb -p DNA.top -water spce -ignh
editconf -bt octahedron -f DNA2.pdb -o Complex.pdb -d 1.0
genbox -cp Complex.pdb -cs spc216.gro -o Complex_b4ion.pdb -p DNA.top
grompp -f em.mdp -c Complex_b4ion.pdb -p DNA.top -o Complex_b4ion.tpr
genion -s Complex_b4ion.tpr -o Complex_b4em.pdb –nname Na –nn 8 –g trp_ion.log
For the rest of the steps, can I continue from: http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/05_EM.html
I am stuck at this step
Yes, parmbsc0 is still not included in new version. So, you have to add it yourself.
You can use following command to add 0.150 M NaCl :
genion -s Complex_b4ion.tpr -o Complex_b4em.pdb -conc 0.150 -neutral
It will add sufficient number of Na+ and Cl- ions such that ion concentration will be 0.150 M and whole molecular system will be neutral.
I am working with an enzyme which has Zinc in its catalytic site.
I want to further continue my study using molecular dynamics simulation, as the zinc is involved in the catalysis process. I tried the simulation with dummy atoms approach but the mass gets distributed between dummy atom (DZ) and zinc, and has an effect on catalysis.
So if anyone has any idea about zinc parameter please let me know.
There are parameters for Zn in the CHARMM, AMBER, and GROMOS force fields. There is no zinc in the OPLS force field files that come with GROMACS, so I don't know about it. If you are really interested in the reaction, you will have to do QM/MM and will then not need the parameters since the zinc will be in the QM region.
g_rms indicates: Each structure from a trajectory (-f) is compared to a reference structure. The reference structure is taken from the structure file (-s).
So you need to put the crystal structure under the flag -s and the structure that you want to compare with the crystal structure under the flag -f.
Then, when you select the groups, the first that the tool asks for is the group for least square fit (so, the structures under the -f flag) and the second is the group for RMSD calculation (structure under the -s flag).
For better information you can see section 8.9 of Gromacs Manual.
Using NAMD, I want to simulate aquaporin channel in POPC membrane. To induce the water flow, pressure difference across the membrane is needed. I am stuck with this point. Could you please suggest how to do it?
Actually using PBC You cant just use higher pressure on one side of membrane. Recently I heared a conference talk, where guy did it by putting two membranes in system - obtained a kind of pseudo-liposome where water pressure between two membranes was higher than outside.
I want to simulate a protein using PDB file extracted from Protein Data Bank. After inspection, I noted that the protein atoms in the PDB has ANISOU atoms. In my understanding, this type of atoms has two coordinates as they are flexible (correct me if I am wrong).
My question is, can I run a molecular dynamics on this PDB file? Do I need to do any pretreatment?
Jul 9, 2013
you have to actually remove the ANISOU atoms data from the PDB file before going for simulation.
You can use the following python script to remove them
I want to simulate Piroxicam (NSAID drug) molecules with GROMACS. I have the GROMACS topology file generated by PRODRG server. As mentioned here http://pubs.acs.org/doi/abs/10.1021/ci100335w the topology generated by PRODRG server contains numerous inaccuracies. I want to validate this topology. How can I do this?
I think the easiest way is try to construct a gromacs topology file using pdb2gmx , using your PRODRG generated file as input and selecting one of the gromos forcefield options ( as the paper that you cite is based on this ff). If works, will mean that the topology generated with PRODRG is valid for run a simulation with gromos.
But without doing this, and only reading the abstract of this paper, I think that you are not going to be able, at elast not in a good way, but in the paper are proposing several points as “best practices” for parametrization of small molecules under the GROMOS force fields, try to take advantage of them.
The higher, in absolute value, of delta G, the better the energy of the protein-ligand complex. However, it doesn't correlate with the activity presenting different ligands. A priori, what should correlate, is if you have a ligand A, and the ligand B, the difference in whose delta G (AAG a-b = | dG(a) - dG(b)| ) should correlate better with the activity. No values a priori can be told, as depending of the size of the binding site and the size of your ligand, more contacts, hence, better or worse dG
If the complex protein-ligand is stable, then, molecular dynamics is required. What you need to do is monitorize the RMSD of your ligand, but it won't give so much information.
Essentially, you need to combine dG and RMSD information both together. You could have very good binders that present a good energetic profile (dG) but are less stable (RMSD high) and viceversa. At the end, both descriptors could give you a hint about if the ligand is stable and energetically well stabilized by the protein
I am a new user of GOMACS. My system contains Polymer and Water. I used gromos53a6 force field and spc water model and it ran nicely. But now I want to use TIP3P water model with the same force field (gromos53a6), in that case, I am getting a error message - "Fatal error: Atomtype OWT3 not found".
In addition, I added OWT3 and HW at the bottom of the .atp file and also corrected nonbonded.itp file of the corresponding force field. But I couldn't get rid of such kind of error.
I need the complete time series of a protein dynamics (hopefully also provided by someone wanting to collaborate). The main aim of our group is to extract topological descriptors of the corresponding protein contact network to be compared with energetic variables of the MD.
The adaptive intermolecular reactive bond order (AIREBO) potential is able to model covalent C-C bond breaking with associated changes in hybridization of atomic orbitals. In addition, the AIREBO potential reproduces the bonding characteristics (such as binding energy and bond length) and forces associated with the rotation about dihedral angles for carbon double bonds. However, If you are trying mechanical study of carbon family, the potential file should be slightly modified or it predict large strength before fracture.
The average structure calculated during analysis would be the average of all the atomic positions/co-ordinates during the length of simulation. But, I wanted to know how are the bond angles treated for computation of the average structure? Since, we use the structure for further analysis too. Are the bond lengths distorted or are they retained while the average is being computed? Also, what are the co-ordinates used for this computation?
@Tarak: Thank you for the explanation. So, if my force field parameters are correct, then the average structure computed via g_rmsf fucntion of gromacs should give me a structure which obeys the standard chemical conformations, i.e. the bond angles and bond lengths? Also, could you please guide so as to how to check if the average structure is correct as obtained?
The MD simulations using gromacs 96 force field was done on a compound obtained after semi-flexible docking and flexible docking. The protein was also simulated in water. Comparing the docking results with the original protein structure the total energy was high, so should I conclude that after docking the structure was unstable?
It's rather difficult to say. It could be yes, but it could also mean that the binding pose and the sidechain orientation, although good for the ligand, it's not ideal for the system. It will depend on the system, because it is to be expected an energetic barrier, but it cannot be calculated like that. Only with metadynamics or free energy profiles.
Like all theoretical models, they are good enough to start with, but if you need accuracy then you need to measure it yourself. If you are unsure about relative reliability of various computational methods, it's common to start with the most popular (so even when it's wrong it can be compared to everyone else who used the same software), and then get some experimental data to refine your estimate. Just remember, all models are wrong, but some are useful.
I am using Gaussian for molecular optimization and energy calculation but now I want to learn how to do the same for dimers. In the case of dimers, orientation and different types of interaction play a major role in optimization of system, so how to deal with all these parameters?
Jun 26, 2013
don't be to affraid of studying these systems. Thanks to DFT-D3 it is much easier nowadays then some years ago. In the beginnig I would suggest you prepare some input structures with different intermolecular orientations (angle, slide, distance), this will help you to see if you end up in the same minimum and to compare the energies. Regarding solvation, why not starting in vacuo to get a feeling for these calculations and later include (maybe different) solvent models? HF will not do a good job for vdw / dispersive interactions, here you need at least MP2 in a reasonable basis set including CP, but forget about that with G..
With DFT you can use BLYP with Grimmes Dispersion correction and Density Fitting to make it faster and stepwise increase the basis set size. Counterpoise correction is not recommended for DFT-D, use the QZVP basis set for the final single points, this will give you essentially no BSSE. If you want to give a try to a different software you can try Turbomole, where each system will take some days on a single core...
I don't know what software you are using but it sounds like there is an atom with an unfilled valence (i.e. not enough bonds attached to it, but also hasn't been assigned a negative charge). Or it might be that it is expecting a specific functinality on the atom but sees a hydrogen instead (i.e. it expects methyl group but there is a hydrogen instead.
It all depends on the specific software you are using and the test it is performing that your molecule is failing, but certainly at first glance it appears as though part of you molecule is incorrect in some way.
I have seen in a lot of papers the probability distribution of a particular eigenvector and comparison with Gaussian PC for the same set of analyses. They have projected the histogram of the fluctuations. I was trying to get the same for my set, I do have the rmsf projection along individual eigenvectors, say 1-8 but I was little confused with the histogram of the fluctuations as to how to obtain those.
Also, for obtaining the Gaussian PC from a list of eigenvectors, is there a straightforward measure (which I have obviously missed reading) to find which one is Gaussian in my set, or does one have to go through hit and trial along with the eigenvalue spectrum? Say, choose PC 100 and plot the projection and obtain the histogram and then see.
sorry for the delay in answering your question. I had a very busy schedule in the past weeks. I took a look to the paper you mentioned and now I have a better understanding of what you need.
Firstly, about producing histograms from the fluctuations around some projection over principal coordinates you have, of course, a number of possible solutions. Usually I produce histograms with Octave (a free version of Matlab). Say that you have an ASCII text file ("data.txt") containing 2 rows of numbers. The first being time, the second being the rmsf of the projection along some eigenvalue. These are the basic commands to plot and save an histogram (in what follows ">>>" means the Octave prompt). I will assume that you launch Octave from the same directory where you have the "data.txt" file.
>>> a = load("./data.txt");
>>> [h, x] = hist(a(:,2),50,1.0);
>>> M= [x' h'];
>>> save hist.dat M;
The command "hist" creates an histogram using data contained in the second row of the matrix a [a(:,2)], with 50 bins and normalizing to 1. I suggest you to refer to Octave user guide about the hist command.
The plot command just plots for you the histogram.
The following commands are used to save your histogram into a file named "hist.dat", which is a 2-column ASCII file with the first row being the value of the rmsf associated to the n-th bin (n = 50 in the example), while the second row contains the value of the histogram in that bin.
About the "Gaussianity" of the distributions, I think that there is no other way than trial-and-error. There are different ways to establish if the distribution is Guassian or not. One possibility is to fit the histogram with a Gaussian as both Haiguang and I suggested before. Here I suggest you a solution using Octave (not gnuplot as I wrote in the previous comment).
This is the idea: if your variable has a Gaussian distribution
p(x) = N exp[-(x-x0)^2/2s^2]
where N is the normalization constant, x0 the average value of x and s^2 the standard deviation, then the natural logarithm of p(x) is a second order polynomial
ln(p(x)) = A + Bx + Cx^2
where A, B and C trivially depend upon N, x0 and s^2. What you can do, after having produced the histogram (and so the x and h arrays) with the previous Octave commands is to do a polynomial fit or order 2 of ln(p), recalculate the value of the distribution with the fitting parameters and obtain a chi-square value. The closer the chi-square is to 0, the closer your distribution is to a Gaussian one. Here are the Octave commands:
>>> z = polyfit(x, log(h + 1.0e-15), 2);
>>> y = poly ], which allows one to treat different portions of a molecular system with different ab initio many-body methods. With CIM you choose the portion that is chemically active to be treated with the highest level ab initio many-body method, while the rest that is not chemically active, and will not play a significant role in determining the correlation energy, is treated with a lower level method.
A good example of how this can be used for real systems is P.M. Kozlowski, M. Kumar, P. Piecuch, W. Li, N.P. Bauman, J.A. Hansen, P. Lodowski, and M. Jaworska, \The Cobalt-Methyl Bond Dissociation in Methyl-cobalamin: New Benchmark Analysis Based on Density Functional Theory and Completely Renormalized Coupled-Cluster Calculations," J. Chem. Theory Comput. 8, 1870-1894 (2012). There are some nice illustrations that demonstrate how the chemically active portion is chosen, and most importantly, it is chosen not by the user, thus making it easier to use.
I know this is different from ONIOM, but the thinking for choosing the higher-level and lower-level portions are the same. Which portions are chemically important and which are mainly spectators, if you will. Also, if you're interested, CIM will be freely available in GAMESS most likely in the next version that is released.
There is a very simple way of articulating this convention. Objects roll down-hill, the negative sign just says that: Force is in the direction of DECREASING potential energy; this convention makes potential energy diagrams intuitive!
The protein is of 543 amino acids, its MD simulation is done using Gromacs, the problem is occurring when the duration is set for 1000 PS and contemporary steps. The problem occurs in 457 step showing "molecules too much far in the trajectory". Could any one suggest how to overcome this issue?
Actually I think you can do it, why not. The only problem is to find a suitable force field to reproduce hematite accurately. It is a pretty common iron oxide, so you will probably find something in the literature.
Regarding the surfaces, yes, you should create them with a visualization tool or with a drawing tool. The second one is easier of course. I recommend you the free software GDIS, or the commercial code Material Studio if you have access to it.
In terms of how the force fields compare w.r.t structure and dynamics, rather than speed factors, there have been a few recent articles that compared a number of the more popular force fields and compared them against experiment.
How can I measure the formation of water mediated hydrogen bonds formed between a ligand and a protein complex during simulation. I'm using gromacs 4.5v. Is there any specific way by which we can calculate the number of hydrogen bonds?
Water mediated hydrogen bonds are very difficult to monitor in a molecular dynamics simulation. 1) They are much less stable (around 4 kcal.mol⁻¹) than protein-protein hydrogen-bonds (~8 kcal.mol⁻¹), and 2) as they are more labile, they break more easily and are rapidly replaced by another water molecule. You have to select a particular residue and all water molecules to find out how many HB are established, but you will not find which specific water molecules (and how many) are involved). If your water molecules are relatively stable and 'trapped' between ligand and protein, maybe you can check them. If not, is more accurate to chech ligand-protein interactions (with LigPlot software) that focus only in water molecules...
For further informations on g_hbond, check: "Van der Spoel, D., Van Maaren, P. J., Larsson, P., & Tîmneanu, N. (2006). Thermodynamics of hydrogen bonding in hydrophilic and hydrophobic media. J Phys Chem B, 110(9), 4393–8. doi:10.1021/jp0572535"
I began to run a MD simulation using the GROMACS program. In the first step, pdb2gmx, it was asked to select the force filed. How do I select the right force field?
Dec 19, 2012
Many force fields can be used to study protein-ligand interactions, for example CHARMM with CGenff, AMBER with GAFF, OPLS, etc.). The choice of the force field for your simulations *should not* depend on the MD code. It is particularly true if the program you want to use accepts different potentials (such as GROMACS, NAMD). To choose the better force field for your simulations, you will need to answer some questions and also read the literature !!! For example:
- Does my ligand or my protein have been already studied with MD simulations. if yes, with which force field ? If not, can I used (say GAFF and AMBER parameters to model accurately my ligand and my protein, respectively?) If yes, good !! , Try these force fields in your simulations and see if you obtain accurate results with them. If not, try different combination of force fields . If all the parameters available in literature are not able to model accurately your protein-ligand interactions, you will need to parameterize your own force field. It is a long and complicated task !!!
Check GROMACS (www.gromacs.org). In my profile, you will find an article about membrane protein simulation (DOI: 10.1021/ct300083m), more specifically ABC transporters (P-glycoprotein) on a POPC membrane. This paper have all the steps/methods you will need.
Any help on membrane building or any other issue, please ask.
I want to simulate dobule stranded DNA using Gromacs. I read that I should use an amber forcefield for DNA simulation and I have crossed a website explaining some parameters that need to be changed in the Gromacs package to become compatible with this amber force filed. Is it possible to use the pdb2gmx command on PDF file for DNA (taken from pdb data bank) and choose amber forcefield for that? Are those changes necessary? And if pdb2gmx does not give error, how can I verify the correctness of the procedure?
Make sure that the double strands do not come apart. Some DNA simulations will require constraints to keep the two strands together.
DSSP instalation error.
Dec 7, 2012
I want to do secondary structure analysis using DSSP. I am getting an error during the execution.
error command: Failed to execute command: /usr/local/bin/dssp -na dd02brnL dd1YGM3g > /dev/null 2> /dev/null. . Please help me out.
I will be glad to learn from anyone who has had experience (or knowledge) of how to estimate configurational entropy from MD trajectories. In particular, I will like to obtain an estimate from cluster analyses obtained from trajectories from different force fields.
Dec 3, 2012
I have recently written a tool to calculate angular spread of single residues from MD trajectories In fact the output is a "flexibility spectrum", where the flexibility of each residue is used to characterize the given trajectory. The tool can also distinguish local fluctuations from conformational changes. Let me know if you want to try it.
I have a question about Material Studio. Actually I'm a beginner to deal with this software. I would like to use forcite calculation to do my simulation. I have a problem with this calculation. I would like to use COMPASS as a force field but one of my compounds (Boron) can't be found in this force field. Then, I tried to change the force field using Universal Force Field and it works. My questions are:
1. What is the different between using COMPASS and Universal Force Field for MD simulation?
2. If I want to use COMPASS Force Field, it means I have to add the compound (Boron) to the database. How to add the database of this compound? Is there any equation that we have to add?
3. Do you have any suggestions on how to do this simulation? Should I change to another method (GULP or the others)?
U must first install ambertools (its fortunately free!). Then make your .prmtop and .crd files using tleap or xleap (included in Ambertools). then use this file to convert .prmtop and .crd files to .top and .gro files
How to retrieve the structures for the energy minima from the generated Gibbs free energy landscapes in Gromacs?
Nov 28, 2012
I have retrieved the Gibbs free energy landscapes for my principal components and located the energy minima. Can anyone suggest how to retrieve the structures from the specific energy clusters on the free energy landscape? I have attached the paper which I am using as reference.
I want to simulate a modeled protein which is around 7 angstrom from native and I want to see whether it goes towards the native with time. What could be the parameters so that the protein does not change much?
Nov 28, 2012
As you probably known folding of small protein with classical explicit MD simulations (without use enhanced sampling approach, such as T-REMD, metadynmics, etc..) take a (very) long time and need a lot computational resources.
For me, It normal that RMSF of C and N terminal are highers, since generally these parts are more flexible than other part of the protein. To conclude I will suggest you to do more refinements of your modeled protein and carry out other MD.
Based on what criteria should we select different temperature levels while simulating a point mutation to study the folding pattern?
Nov 23, 2012
I want to perform MD simulations for two point mutations of my protein at two or three different temperature levels and compare the folding at different temperatures. I want to know whether there is any criteria for selecting different temperatures. Any suggestions?
You could predict the temperature stability change of your mutated protein and choose temperatures for simulation below and above the predicted melting temperature.
Tools that predict the melting temperature change according to mutations are:
Between the ligand and protein, unfortunately gromacs don't have any tool. But, you can download LigPLOT (or LigPLOT+ with GUI), a simple program that quantifies the number of non-bonded interactions and hydrogen-bonds between any ligand and the protein.
I don't have expertise on SYBYL, but in GROMACS you cannot break chemical bonds (as seen in a chemical reaction). The only way you can simulate a similar effect is a multi-step quantum calculation with GAUSSIAN or employing a QM/MM simulation, with explicit solvent that you can change. But, in QM, solvent effects can be difficult to simulate.
Dear Mr. Penumutchu, hear are my answer for your questions. Hope these answer satisfies you.. :)
1. If you conduct a MD simulation for protein complex after docking, you may predict its stability based on the curve of simulation time versus potential energy of the complex. This complex stability calculation thus very important to obtain the best ligands or lead compounds.
2. To answer this question maybe you have to see the software manual. Usually, to analyze the MD simulation results you have to draw curve between potential energy of the complex versus simulation time or RMSD/RMSF (root mean square deviation/fluctuation) versus simulation time.
3. You could compare the RMSD/RMSF curve between simulation result and experimental result from protein data bank (PDB) and then you could conclude is your simulation result valid or not.
I am doing molecular dynamic simulation of one compound to see the binding mode in its protein target, but the absolute configuration of one chiral center of this compound changed during the dynamic simulation. How can I fix its absolute configuration to be unchanged?
HumanNet scores were designed to indicate pairs of genes that are more or less likely to share the same GO functional annotation (GO "biological process"). For an example of construction of one of the individual features: for CX, or co-expression, lots of microarray experiments in different conditions, calculating the pearson correlation between all pairs of genes for their expression over different conditions, and using regression to get a score from those correlation values related to GO functional annotation. The individual features were integrated based on their power in predicting GO functional annotation, resulting in an "integrated" score, the last column in the dataset.
I am trying a MD simulation on fructose (as a practice) using VMD and NAMD. I am having difficulties in generating the *.PSF file. I am always getting an error saying some atoms can’t be parameterized. Can you please help me to fix this issue?
I am using toppar CHARMM36 all-atom carbohydrate force field files.
Hi, I personally think glycam force field is better for carbohydrates. and u can get some help in Amber mailing list
Can anyone suggest any tools or softwares which would aid in visualizing the eigenvectors obtained from PCA of gromacs trajectories.
Mar 19, 2013
I want to visualize the eigen vectors individually and compare the relative motion. Initially i tried the dynamite server , but it is down from past few months. so, if anyone is using a similar or better tool please let me know
You can get trajectory for the given eigenvector through g_analyze using "-extr " option and frame number in this trajectory is control by "-nframes" options. Then, you can use VMD to visualize the trajectory, which will give you the motion along the given eigenvector.
Could you please help me out how to do MD simulation using Gromacs for a protein which has Threonine residue being phosphorylated (TPO). Ending up with Error "Residue TPO not found in Residue topology File."
Mar 19, 2013
For proteins containing the phosphorylated residues, you need to edit the forcefield so that the new residues can be recognized by the forcefield or u can download either one of this forcefield (ffG43a1p.tar.gz/gromos43a1p-4.5.1.tgz)
If the metal has any kind of activity in the active site, as far as I know, any docking program can calculate that. In contrast, in the case of the metal has only "structural effects" (such as sustain the ligand which participates in the reaction), you can do it doing some tricks as change the name to the metal from X to Zn (in the case of GOLD, that works) or create a new atom type. In one of the works where I collaborate, my colleges use a "huge proton" with 6 valences... it's chemically meaningful but can model the structural effect of the metal.
Getting error in autodock execution for autogrid running.
Mar 17, 2013
When I am running autogrid run, I am getting error that I cannot open or find Grid parameter file, I am using Windows 7.
Thanks @Rajnikant for ur kind response, I have some more doubts,
I am facing problem when performing docking of metal with protein, can any body suggest me how to make metal pdb file, when i am running autogrid of complex, i am getting error as unknown ligand atom type x, add parameters for it to the parameter library first. what should i do now,
I am using Win 7,I have one more doubt like, do we need cygwin for autodock 4.2, or not, as I have cygwin and while execution i am following the path from directly the scrippt research institute instead of from cygwin folder.
I have saved my input files in user in one folder, do I need to keep it in C drive itself the one which I am giving to program pathname for autogrid. exe selection at time of execution.
I have found that it is the spatial fourier transform of gs(r,t). But gs(r,t) is calculated at a particular time. So what time should I select? And then for the fourier transform which k value should I use? Since I am not familiar with fourier transformation can you suggest me which algorithm should I use?
Either you can do the Fourier transform of this quantity or you can directly calculate intermediate scattering function by the formula attached here. The <> averaging is the ensemble average which is generally being done from the time trajectory given that your system is ergodic. Basically, one can theoretically derive those expression from the spatial fourier transform of pair-pair correlation function, gs(r,t).
If you want to follow the direct path (calculating from gs(r,t)), then you can choose the time interval 0,5, 10, 15 or any regular timestep interval upto a large value - lets say 100000 timesteps. Please remember you have to have the trajectory snapshots which are consistent with these interval.
I want to have comprehensive information on the force fields present in gromacs, in the PDB2GMX command.
Any article on comparative analysis of ff will be welcomed.
Also, info. on which ff to use on a particular system and its explanation would be useful.
Mar 6, 2013
GROMACS contains a lot of force fields for many types of biomolecular simulations (i.e. protein, detergents, lipids, ionic solutions, etc.) so it is difficult to give you a definitive response. Be more precise in your question. For protein force fields you can see the reference attached:
In case you second question, there is no clear/definitive response since it depends a lot of what do you want to do. A good start is to (always) read the literature about previous MD works of similar system.
I am trying to run a quench molecular dynamics on a molecule. This molecule consists of Cyclopentadienyl rings with Fe together with a (CH)2 chain. I am using Materials Studio 6.1 (Accelrys) version program package.
The problem is, when I click on the run bottom then it is showing "the dynamics can't integrate the equation of motion for a particle of zero (or small) mass". But if I cut the Cyclopentadienyl rings together with Fe then there is no problem. That means there is no problem with the rest of the (CH)2 chain.
The problem that you are facing is that you have a cyclopentadienyl ring which has a dummy atom representing the Pi system. This is not recognised by Forcite or Discover as a valid particle for molecular dynamics because it has no mass associated with it and hence the program gives an error.
I think that it is pretty hard to study MD of organometallics using classical simulations techniques. I did a quick search for "organometallics" on
This turned up many hits but publications were nearly all using DMol3. I guess you could potentially delete the Pi system and create bonds between the cyclopentadienyl fragment and the metal atom but I am not convinced this would give good results (picking a forcefield might be challenging!)
I have read a couple of tutorials and am trying to incorporate the basic analysis tools available in gromacs like g_rmsd, g_rmsf, g_rms, g_covar, etc. But as I have recently started, I am confused about how to proceed in an appropriate direction. I am interested in studying the allosteric movements in my protein. Any help for me to channelize properly would be really kind.
Feb 27, 2013
Sorry but your question is not very precise. A good start is to read the literature about similar system as yours Probably you have ligand(s) in your system (correct ?), so you could examine the influence of the ligand on the protein dynamic against the apo protein. For instance, by examining
- the stability of the protein vs. time : g_rmsd
- the fluctuation of the protein domains w/o ligand : g_rmsf
- the protein domain motions w/o ligand : g_covar
- Interaction between protein residues (in the pocket ) - ligand : g_hbonds and g_saltbr
Hi, a simulation lenght depends from what information you want to get from your MD.
Generally a MD simulation should be carried out untill the property you are looking for reaches convergence. For example, if you are interested in generating a set of snapshots for MMPBSA analysis, the property you might look at can be the protein RMSD, so you can run your simulation untill the RMSD vs time plot reaches a plateau. If you start from a well defined crystal structure, generally you'll reach convergence in a relatively short time. Otherwise, if you have to test different ligands docked in the same protein structure, simulation time will depend on how the protein need to "adapt" to the different ligands.
Differently, if you need or want to observe large protein motion following a particular "perturbation" (i.e. activation following the binding of an agonists to an inactive state of the protein), simulation time might be decidedly longer.
By different time scale such as 5000ps, 10000ps etc. for structural analysis by different server.
Feb 23, 2013
You don't need a webserver to do that, you can use the g_angle, g_hbonds tools of GROMACS with the appropriate snapshots and index file of your MD trajectory, since you did your MD with GROMACS. Go to the gromacs website and start the helps of these tools.
Additionnally, i don't think you did 5000 ns or 10000 ns of MD ;)
I want to get the psi and phi angles on a particular residue and plot them on ramachandran plot. I went through some tutorials and was informed to use ndx file for that. But can someone give me an actual command to do so.
Mar 8, 2013
I think this is what u are looking for
To extract the angles for a single residue, e.g. ILE-32, type:
How to retrieve a phosphorylated form of a protein insilico to perform molecular dynamics studies.
Feb 23, 2013
In my protein, tyrosine kinase (TYR) is phosphorylated. I have the structure of the normal form. Is there a server, tool, or any other method to replace the TYR with phospho-TYR? I want to use this phosphorylated form for my MD studies.
I am guessing GROMACS does not recognize your phosphotyrosine modification because its forcefield parameters are not defined. You can download Phosphotyrosine forcefield parameters for GROMACS here:
Its a modified GROMOS96 43a1 forcefield that includes phosphoserine, phosphothreonine, and phosphotyrosine parameters
Briefly below are the steps you need to take.
1. untar the downloaded file in your current directory where you setup your simulations.
2. You can use a variety of tools like the ones suggested by others to add the phosphotyrosine residue into your PDB file
3. To let GROMACS recognize the residue correctly, the important thing is to make sure the phosphotyrosine residue in your input PDB file has its residue name as PTR and the atoms are named as listed below and strictly in that order.
NOTE: If the phosphotyrosine residue you inserted using VMD or any other tool has a different atom naming convention, you can manually edit the PDB file and change them to the ones listed above.
NOTE: Explore the downloaded forcefield directory to find PTR residue listing in the ffG43a1p.rtp file.
4. using this modified PDB file as your input run pdb2gmx -ff ffG43a1p <including the other usual command line options)
You should now be able to interactively choose the forcefield ffG43a1p to include phosphotyrosine in your simulation box.
Hope this helps
Begin with a fairly sparse k-point grid, or even Gamma, and increase the energy cut-off in steps of 20% from a low value until the total energy of the model system stops decreasing on each step by, say 0.1 meV. Then, increase the grid density monotonically until the total energy stops changing by the same amount as chosen previously for the convergence criterion, and that's it.
I am working on RNA simulations and having a problem. When I start simulation from full stretched sequence, it never goes to native structure (a T-loop). It is a very small RNA (20 nt), and it is never showing any base-pairing during simulation. I am using amber99sb.
Can anybody suggest to me something which could be helpful for my work?
Feb 14, 2013
I think the problem is your force field. The AMBER99SB force is not optimized for DNA and RNA simulations . It is mainly (but not only) used for proteins. For simulations of RNA/DNA, I would suggest you to use the recent bsc0χOL3 force field (J Phys Chem B 116(33):9899-916 (2012), PMID 22809319)
For more information, you can also read for instance the following papers
Effect of Guanine to Inosine Substitution on Stability of Canonical DNA and RNA Duplexes: Molecular Dynamics Thermodynamics Integration Study
Miroslav Krepl, Michal Otyepka, Pavel Banáš, and Jiří Šponer
The Journal of Physical Chemistry B
2013 117 (6), 1872-1879 DOI: 10.1021/jp311180u
The Nuclear Magnetic Resonance of CCCC RNA Reveals a Right-Handed Helix, and Revised Parameters for AMBER Force Field Torsions Improve Structural Predictions from Molecular Dynamics
Jason D. Tubbs, David E. Condon, Scott D. Kennedy, Melanie Hauser, Philip C. Bevilacqua, and Douglas H. Turner
2013 52 (6), 996-1010
Performance of Molecular Mechanics Force Fields for RNA Simulations: Stability of UUCG and GNRA Hairpins J. Chem. Theory Comput., 2010, 6 (12), pp 3836–3849
There is no easy answer for this question. The amount of time required for a simulation depends a lot upon the properties that you want to study, and how long the system will need to provide you with sufficiently high-quality data that your answers will be reliable. A time appropriate for one system may be completely insufficient or overkill for another simulation of a different material or under different state conditions.
There are a few basic rules of thumb that you can follow:
(1) You need to an equilibrate a system until there is no discernible trace of how you "built" the system in the data that you're collecting.
(2) Your "production runs" should be long enough to reduce the uncertainty in your your statistical averaging to an "acceptable" level. (You can specify what that level of uncertainty is. The uncertainty in the averaging can be predicted using classical tools, or can be estimated with methods such as those developed by Flyvbjerg and Petersen.)
In my opinion Schrodinger is the best option if you can afford it. However, if your budget is limited than ICM (Molsoft) is one of the best options. It costs only around 3000$ and can achieve better docking results than Glide. Other good options are Molegro's MVD for docking and MDM for QSAR. MVD is also CUDA ported. I personally tested it and the CUDA version achieve really good results and 17x speed up!
I am new to simulations and gromacs. I have docked the protein-ligand complex and now want to simulate it using Gromacs. I am using the best docked conformation of protein-ligand structure. Should I remove hydrogens from it, as hydrogens are causing problems, and ignore the missing residues and generate the topology of protein ( and generate the ligand topology by PRODRUG server)? Is this the correct method?
That is no missing residue. Simply, in GROMACS, atom names must match. For instance, if the acidic proton in serine is named H, but in the database is named HA, they don't match and give that specific error.
Open the aminoacid database (aminoacids.rtp) of the correspondent FF and check the atom (probably HA). Then, go to your protein and see if the correspondent proton is also HA:
- If not, change it to HA;
- If missing, you have a protein protonation problem. Re-protonate the protein.
But I faced a problem in the shrinking phase. I have done up to 14 cycles of the shrinking phase using inflate gro. I got only 2.14513126335007 nm^2 area per lipid. The protein was in the membrane up to 13 cycles but after the 14th cycle it came away from the lipid bilayer. I was unable to find the error. Please help me in this regard.
It is possible, but will depend on the flexibility of your receptor. With MD, the whole structure will adapt towards the perturbation made by each specific ligand. However, it can be difficult to calculate binding energies. One way is through Gibbs Free Energy: if you simulate ligand and protein alone, and then the complex ligand-protein, you can calculate the difference between energies and extrapolate the binding energy. Another thing that you can search is for evidences of 'induced-fit' mechanism between ligand and receptor.
I have a protein (crystal structure) complexed with two molecules of the same ligand (inhibitor). I am trying to build a pharmacophore model for searching potential hits. But the problem is that two molecules are bound in same binding pocket but in different conformation (image attached). Both the molecules are interacting with hydrophobic interaction. The binding pocket consists of exclusively hydrophobic residues. Could some one please suggest me how can I develop a pharmacophore model by using these? Shall I straightforwardly consider any one of the ligand molecules? I was thinking of merging two similar pharmacophore features to make one point. But the problem is that the binding conformation of both the inhibitor molecules is different.
Also, if I want to do docking of several inhibitors in the same site before going for pharmacophore modelling (for building common pharmacophore model), how should I do this?
Please help me in solving this. Thanking you in anticipation.
Given your ligand and the pocket size, there may not be a simple answer or a straightforward way use existing pharmacophore generating software. You have to decide if you want to find a ligand that is of similar size, where two of them will bind next to each other. Or do you want to find a single ligand that spans the whole binding site region.
An initial place to start is to try to find similar ligand based on simple molecular descriptors (e.g. via ShaEP, Feature trees, MOE’s MACCS). Carefully selected descriptors might allow you to find a few similar compounds that would bind and stack the way the native ligand does.
However, I would recommend starting even more basic, and form some hypotheses that you can computationally test. A better understanding of the system may help you develop a suitable approach for searching for new ligands. By investigating the system, you will develop some chemical intuition concerning it. Some basic questions I would have for you system are
1) Could the ligands adopt different binding poses under physiological conditions?
2) What are the ligand-protein binding strengths of each ligand alone and when together?
3) What is the ligand-ligand interaction energy?
4) How dynamic is the binding pocket?
5) What are the more dynamically stable ligand-protein interactions?
6) How important are all of those sulfur atoms (if that is what the yellow atoms are), and what changes occur if they are replaced with oxygen atoms?
7) What forces are important for binding? How large is the molecule’s dipole moment?
8) How conformationally labile is the molecule?
Insight into many of these questions can be gain through good MD simulations. Such an investigation would be time consuming, but potentially rewarding.