For my 10 ns MD simulations for my protein-ligand complexes, I first created .prmtop topology and .inpcrd coordinate files in tleap with a 12 A water box (using the solvateoct command) and added ions to neutralize the complex. Then using pmemd, I ran minimization, heating, and equilibration steps similar to some of the online Amber tutorials I found before I did 10 ns of production and created .rst files and .mdcrd files. But I noticed that for some of my complexes, the RMSD didn't converge completely after 10 ns because there were very sharp sudden peaks in my RMSD graphs that I created with cpptraj. I've attached a few images of my RMSD graphs here to show the difference. I've been looking through Amber tutorials and the manual for months, but am still having issues with knowing the right settings for my MD simulation, so any advice for a beginner would be much appreciated!
Here are the settings that I used for my input files:
min.in:
&cntrl imin=1, ntx=1, ntmin=1, maxcyc=50000, ncyc=10000, ntpr=100, ntf=1, ntc=1, ntb=1, cut=10.0 &end
heat.in
&cntrl imin=0, irest=0, ntx=1, ntb=1, cut=10, ntr=1, ntc=2, ntf=2, tempi=0.0, temp0=300.0, ioutfm=1, iwrap=1, ntt=3, gamma_ln=2.0, nstlim=50000, dt=0.002, ntpr=1250, ntwx=1250, ntwr=2500 / Keep protein fixed with weak restraints 100.0 RES 1 208 / END END
equil.in
&cntrl imin=0, irest=1, ntx=5, ntb=2, ntp=1, pres0=1.0, taup=2.0, cut=10, ntr=0, ntc=2, ntf=2, tempi=300.0, temp0=300.0, ntt=3, ioutfm=1, iwrap=1, gamma_ln=2.0, nstlim=50000, dt=0.002, ntpr=1250, ntwx=1250, ntwr=2500, / Keep protein fixed with weak restraints 100.0 RES 1 208 / &end
prod.in:
NPT production &cntrl imin=0, ntx=5, irest=1, ntpr=1250, ntwr=2500, ntwx=1250, ntf=2, ntc=2, ntb=2, ioutfm=1, iwrap=1, cut=10.0, nstlim=500000, dt=0.002, temp0=300, ntt=3, gamma_ln=2.0, ntp=1, pres0=1.0, taup=2.0, / Keep protein fixed with weak restraints 5.0 RES 1 208 / &end
Here are the commands that I used to run the simulation:
pmemd.MPI -O -i min.in -p complex12solvated.prmtop -c complex12solvated.inpcrd -o min12.out -r min12.rst -x min12.mdcrd -inf min12.mdinfo -l min12.logfile
pmemd.MPI -O -i heat.in -p complex12solvated.prmtop -c min12.rst -o heat12.out -r heat12.rst -ref min12.rst -x heat12.mdcrd -inf heat12.mdinfo -l heat12.logfile
pmemd.MPI -O -i equil.in -p complex12solvated.prmtop -c heat12.rst -o equil12.out -r equil12.rst -ref heat12.rst -x equil12.mdcrd -inf equil12.mdinfo -l equil12.logfile
pmemd.MPI -O -i prod.in -p complex12solvated.prmtop -c equil12.rst -o prod12.out -r prod12.rst -ref equil12.rst -x prod12.mdcrd -inf prod12.mdinfo -l prod12.logfile
pmemd.MPI -O -i prod.in -p complex12solvated.prmtop -c prod12.rst -o prod12-extend2ns.out -r prod12-extend2ns.rst -ref prod12.rst -x prod12-extend2ns.mdcrd -inf prod12-extend2ns.mdinfo -l prod12-extend2ns.logfile
pmemd.MPI -O -i prod.in -p complex12solvated.prmtop -c prod12-extend2ns.rst -o prod12-extend3ns.out -r prod12-extend3ns.rst -ref prod12-extend2ns.rst -x prod12-extend3ns.mdcrd -inf prod12-extend3ns.mdinfo -l prod12-extend3ns.logfile
pmemd.MPI -O -i prod.in -p complex12solvated.prmtop -c prod12-extend3ns.rst -o prod12-extend4ns.out -r prod12-extend4ns.rst -ref prod12-extend3ns.rst -x prod12-extend4ns.mdcrd -inf prod12-extend4ns.mdinfo -l prod12-extend4ns.logfile
pmemd.MPI -O -i prod.in -p complex12solvated.prmtop -c prod12-extend4ns.rst -o prod12-extend5ns.out -r prod12-extend5ns.rst -ref prod12-extend4ns.rst -x prod12-extend5ns.mdcrd -inf prod12-extend5ns.mdinfo -l prod12-extend5ns.logfile
pmemd.MPI -O -i prod.in -p complex12solvated.prmtop -c prod12-extend5ns.rst -o prod12-extend6ns.out -r prod12-extend6ns.rst -ref prod12-extend5ns.rst -x prod12-extend6ns.mdcrd -inf prod12-extend6ns.mdinfo -l prod12-extend6ns.logfile
pmemd.MPI -O -i prod.in -p complex12solvated.prmtop -c prod12-extend6ns.rst -o prod12-extend7ns.out -r prod12-extend7ns.rst -ref prod12-extend6ns.rst -x prod12-extend7ns.mdcrd -inf prod12-extend7ns.mdinfo -l prod12-extend7ns.logfile
pmemd.MPI -O -i prod.in -p complex12solvated.prmtop -c prod12-extend7ns.rst -o prod12-extend8ns.out -r prod12-extend8ns.rst -ref prod12-extend7ns.rst -x prod12-extend8ns.mdcrd -inf prod12-extend8ns.mdinfo -l prod12-extend8ns.logfile
pmemd.MPI -O -i prod.in -p complex12solvated.prmtop -c prod12-extend8ns.rst -o prod12-extend9ns.out -r prod12-extend9ns.rst -ref prod12-extend8ns.rst -x prod12-extend9ns.mdcrd -inf prod12-extend9ns.mdinfo -l prod12-extend9ns.logfile
pmemd.MPI -O -i prod.in -p complex12solvated.prmtop -c prod12-extend9ns.rst -o prod12-extend10ns.out -r prod12-extend10ns.rst -ref prod12-extend9ns.rst -x prod12-extend10ns.mdcrd -inf prod12-extend10ns.mdinfo -l prod12-extend10ns.logfile
To answer your question in short - you cannot.
Convergence of a molecular dynamics simulation has been discussed nicely in the following thread and I advice you to go over it.
https://www.researchgate.net/post/What_do_you_understand_by_convergence_in_molecular_dynamics
Since you used the term "some of your ligands", I assume you have multiple ligands bound to a protein and "some of the ligands" have high RMSD variations ( I cannot see your images). So, those ligands may be not in the "best" starting conformation. Did you use a docking tool to get these initial conformations or a PDB structure?.
Aravinda Munasinghe Thank you for your response and for the link! I've fixed the pictures in my original question. My pictures look similar to the graph in this question actually (https://www.researchgate.net/post/Does_anyone_have_experience_with_unconventional_RMSD) , except I'm using Amber instead of Gromacs so I don't know how to fix this issue.
I'm a little confused from what I've read in the link you posted. So is looking at RMSD convergence not meaningful at all then? My advisor said that this was important for me to make sure that the RMSD of the carbon atoms for my simualtions converged before I could run MMPBSA calculations, and some of the Amber tutorials ( http://ambermd.org/tutorials/advanced/tutorial3/section1.htm ) also looked at the rmsd of the backbone atoms, so I thought this was important?
For my simulations, I used one receptor with 25 different ligands. So I have 25 separate protein-ligand complexes that I am individually running 35 ns MD simulations with using pmemd. Each protein only has one ligand. To get the RMSD of the carbon atoms for each complex, I used cpptraj to analyze the 35 one nanosecond mdcrd files that were created for that complex.
Here are all the commands that I entered in cpptraj for complex 3, for example:
parm complex3solvated.prmtop
list parm
parminfo 0
trajin prod3.mdcrd
trajin prod3-extend2ns.mdcrd
trajin prod3-extend3ns.mdcrd
trajin prod3-extend4ns.mdcrd
trajin prod3-extend5ns.mdcrd
trajin prod3-extend6ns.mdcrd
trajin prod3-extend7ns.mdcrd
trajin prod3-extend8ns.mdcrd
trajin prod3-extend9ns.mdcrd
trajin prod3-extend10ns.mdcrd
trajin prod3-extend11ns.mdcrd
trajin prod3-extend12ns.mdcrd
trajin prod3-extend13ns.mdcrd
trajin prod3-extend14ns.mdcrd
trajin prod3-extend15ns.mdcrd
trajin prod3-extend16ns.mdcrd
trajin prod3-extend17ns.mdcrd
trajin prod3-extend18ns.mdcrd
trajin prod3-extend19ns.mdcrd
trajin prod3-extend20ns.mdcrd
trajin prod3-extend21ns.mdcrd
trajin prod3-extend22ns.mdcrd
trajin prod3-extend23ns.mdcrd
trajin prod3-extend24ns.mdcrd
trajin prod3-extend25ns.mdcrd
trajin prod3-extend26ns.mdcrd
trajin prod3-extend27ns.mdcrd
trajin prod3-extend28ns.mdcrd
trajin prod3-extend29ns.mdcrd
trajin prod3-extend30ns.mdcrd
trajin prod3-extend31ns.mdcrd
trajin prod3-extend32ns.mdcrd
trajin prod3-extend33ns.mdcrd
trajin prod3-extend34ns.mdcrd
trajin prod3-extend35ns.mdcrd
list trajin
trajout 35ns_com3_CACN.trj
strip :WAT
reference complex3solvated.inpcrd
rms reference out 35nscomplex3.rmsd :1-207@CA
rms reference out 35nsligand3.rmsd :208
rms reference out 35nscomplex3.gnu :1-207@CA
rms reference out 35nsligand3.gnu :208
list actions
run
quit
I'm a little confused by how to know if my ligands are not in a good starting conformation; could you explain this a little more please? I first downloaded my receptor protein from the PDB website and prepared it using the Schrodinger Maestro software and then did virtual screening of a chemical library of 10,000+ compounds using the GOLD docking software and picked ~100 ligands with top-scoring positions, and then I used Schrodinger Prime MM-GBSA calculations and picked the 25 protein-ligand complexes with the most negative binding energy for the MD simulations.
Yes, sudden jumping you see in RMSD is more likely due to periodic boundary conditions. In CPPTRAJ you can align your trajectory wrt the initial frame. That should technically help you with this rapid jumping in RMSD.
If this is not the case, you should load this trajectory in VMD and see whats going on in that time period.
And about auto image in cpptraj,
autoimage < protein_selection >
https://amber-md.github.io/cpptraj/CPPTRAJ.xhtml#magicparlabel-4231
about ligands initial conformation - Docking is a good starting point. But remember, computational models just doing what they are designed to do. So, just cz it throwing out some structures doesn't mean they are the best structures. That's why most ppl run MD on docking outputs. Different docking tools can give you different structures (based on parameters they use) for the same system. Its up to you to decide which structures make chemically more sense or employ additional constrains to screen them properly (as you do with MMGBSA).
- Badges
- Science topic