GROMACS version: gromacs-2022.2
GROMACS modification: No
Hi All,
(maybe more relevant to Dr. Christian Blau)
It has been long I did not do MD simulations. I came across an old Bioexcel webinar about density mediated simulation by Christian Blau. It seems to be a very helpful program. I gave a try to fit alphafold protein model in cryo map density mrc format (only single protein map density extracted from a full map which has more than one protein). It does not seem to fit very fine in the map density as some helices are not within the map density. In brief I am doing these steps:
gmx pdb2gmx -f af.pdb -o af_processed.gro # select charmm27
gmx editconf -f af_processed.gro -o af_center.gro -c
gmx grompp -c af_center.gro -f densguidsim.mdp -p topol.top -o densguidsim.tpr -maxwarn 2
gmx mdrun -v -s densguidsim.tpr -o out.trr
for the grommp step I used maxwarn because I have net charge in protein and am generating velocities. My mdp settings are below, which probably is not great for density mediated simulations; it is mainly a NPT simulation mdp set up but at the end I simply inserted density guided simulation mdp settings. The reason I included NPT settings because if I don’t use the settings below the bond parameters my simulation crashes within 300 ps (although I haven’t tried to figure out exactly what other mdp settings are required apart from the density mediated simulation settings so that it does not crash)
I have following questions:
a) for how long one should usually run? I got rough fitting in less than 1 ns, shall I keep it running more for finer fitting?
b) are there any good mdp settings to start with (apart from the section specifically related to density mediated simulation)
c) in the webinar it is suggested that “Pick periodic image closest to density center!”. I am not really sure what does it mean?
d) in the webinar there are nice movies how the model fits into the map. When I open my pbc whole corrected trr trajectory along with mrc file in chimerax I see my protein and map are very far apart. I can bring them close manually using mouse but I am not sure how Christian’s protein model and map appear to be starting close to each other :) please let me know.
This is mdp:
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000000 ; 2 * 50000 = 100 ps
dt = 0.002 ; in ps and = 2 fs
; Output control
nstxout = 5000 ; save coordinates every 10.0 ps
nstvout = 5000 ; save velocities every 10.0 ps
nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps
; Bond parameters
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings hello
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein ; two coupling groups - more accurate
tau_t = 0.1
ref_t = 300
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = yes ; Velocity generation is off
;;;;;;;;;;;;;;;;; density guided simulation settings
density-guided-simulation-active = true
density-guided-simulation-group = Protein
density-guided-simulation-similarity-measure = relative-entropy;inner-product
density-guided-simulation-atom-spreading-weight = mass ; unity
density-guided-simulation-force-constant = 100
density-guided-simulation-gaussian-transform-spreading-width = 0.2
density-guided-simulation-gaussian-transform-spreading-range-in-multiples-of-width = 4
density-guided-simulation-reference-density-filename = my.mrc
density-guided-simulation-nst = 20
density-guided-simulation-normalize-densities = true
density-guided-simulation-adaptive-force-scaling = false
density-guided-simulation-adaptive-force-scaling-time-constant = 100 ; in ps,
thank you
best regards,
Dave