I am trying to follow this tutorial.
https://degrootlab.github.io/pmx/tutorials/protein_mut.html
After all computations up to Non-equilibrium MD are done, I got this plot. I repeated the same calculations and got very similar result which is different from the one on the website. The only change I made was decrease the time step to 1fs. I have adjusted all other parameters so that the data collection was recorded every 20ps. What could be the potential problem I made during this tutorial? Also, I noticed the force constant on one of the carbon atom (CG 100) was very high during the minimization. However, I do not think this is the cause of the issue here. My DG is negative -13kJ/mol in both try.
Thank you very much for your help in advance.
if you changed the time step, did you also change nsteps and delta-lambda?
Thank you very much for the response. Yes, I did change that to 1.25e-5 for 1 ps time step with 80ps runs.
Here is the part in r_nonequil. mdp
; RUN CONTROL
;----------------------------------------------------
integrator = sd ; stochastic leap-frog integrator
nsteps = 80000 ; 1 * 40000 fs = 80 ps
dt = 0.001 ; 1 fs
comm-mode = Linear ; remove center of mass translation
nstcomm = 1 ; frequency for center of mass motion removal
; OUTPUT CONTROL
;----------------------------------------------------
nstxout = 20000 ; save coordinates to .trr every 20 ps
nstvout = 20000 ; save velocities to .trr every 20 ps
nstfout = 20000 ; save forces to .trr every 20 ps
nstxout-compressed = 20000 ; xtc compressed trajectory output every 20 ps
compressed-x-precision = 1000 ; precision with which to write to the compressed trajectory file
nstlog = 10000 ; update log file every 10 ps
nstenergy = 2000 ; save energies every 2 ps
nstcalcenergy = 1 ; calculate energies at every step
; FREE ENERGY
;----------------------------------------------------
free-energy = yes
init-lambda = 1 ; start from state B
delta-lambda = -1.25e-05 ; complete transition in this number of steps
sc-coul = yes ; use soft-core also for coulombic interactions
sc-alpha = 0.3 ; soft-core
sc-sigma = 0.25 ;
sc-power = 1 ;
nstdhdl = 1 ; write to dhdl at each step
yes, that looks OK. Did you also try with the unmidified mdp file? What results do you then get?
It could be that the tutorial used an older/different pmx version. It is best to compare ddG, because dG depends on the topology construction. Could you complete the thermodynamic cycle with the tripeptide and check if ddG is the same as in the tutorial? Also, did you use the same ff as in the tutorial?
I am sorry for the slow response. I just could not run the unmodefied mdp because it crashed several time during the equilibrium run. Also, I just finished run for the unfolded DG part. I will post both result (They are off by approximately -23 or -24kJ/mol, respectively). However, the DDG is approximately +16kJ/mol which is within the error suggested by the tutorial. So, the issue I had maybe relate to update(?) in the GROMACS FF(amber99sb-star-ildn-mut.ff/forcefield.itp). I will try to repeat this unfolded part again to make sure I get the similar number. Thank you very much for your suggestions.
Ok, it seems that all is as expected: ddG matches the tutorial result. The difference in dG is likely due to the changes in the topology construction for the mutations between pmx versions