Langevin dynamics of protein in vacuum

User discussions
1

GROMACS version: 2023-2
GROMACS modification: No
Dear All,
I am trying to implement Langevin dynamics for an in-vacuum protein system with PBC and large cutoffs for Coulomb and VdW interactions (no PME), and OPLS-AA as ff. I have two questions:

  1. When using the recommended value tau_t = 2 ps, which here defines the inverse friction constant, the recorded temperature (using gmx energy) stabilizes at ~265 K, while my set temperature (in mdp file) is 300K. I only made the recorded and set temperature match when lowering tau_t to 0.2. Although I am still seeing the predicted behavior for my protein (compaction, about in a slower fashion), I am not sure about how to interpret the temperature change. In addition, I would expect tau_t for an in-vacuum or gas phase system to be higher, as the inverse friction constant of water should be lower (higher friction constant). How could this be interpreted?
  2. Using Langevin as both integrator and thermostat has been the only way I have got my vacuum and PBC protein system to work with a commonly used force field like OPLS. Using other integrators and thermostats have resulted in the system being very unstable, especially at high temperatures. Has Langevin dynamics been used for in-vacuum protein systems? I have not found published articles about this.
    Below is my nvt.mdp file with tau_t = 0.2 ps:
    title = OPLSAA NVT equilibration
    ; Run parameters
    integrator = sd ;
    nsteps = 1000000 ;
    dt = 0.001 ;
    ; Mode for center of mass motion removal
    comm-mode = Linear ; remove center of mass
    ; Output control
    nstxout = 0 ; save coordinates every 1.0 ps
    nstvout = 0 ; save velocities every 1.0 ps
    nstfout = 0 ; nstvout, and nstfout
    nstenergy = 1000 ; save energies every 1.0 ps
    nstlog = 1000 ; update log file every 1.0 ps
    nstxout-compressed = 1000 ; save compressed coordinates every 1 ps (default 5000)
    compressed-x-grps = System ; save the whole system
    ; Bond parameters
    continuation = no ; first dynamics run
    constraint_algorithm = lincs ; holonomic constraints
    constraints = h-bonds ; bonds involving H are constrained
    lincs_iter = 2 ; accuracy of LINCS
    lincs_order = 4 ; also related to accuracy
    ; Nonbonded settings
    cutoff-scheme = Verlet ; Buffered neighbor searching
    ns_type = grid ; search neighboring grid cells
    nstlist = 10 ; 20 fs, largely irrelevant with Verlet
    rcoulomb = 333.3 ; short-range electrostatic cutoff (in nm)
    rvdw = 333.3 ; short-range van der Waals cutoff (in nm)
    DispCorr = no ; account for cut-off vdW scheme
    rlist = 333.3 ; short-range neighbour list cut-off
    ; Electrostatics
    coulombtype = cutoff ; Particle Mesh Ewald for long-range electrostatics
    pme_order = 4 ; cubic interpolation
    fourierspacing = 0.16 ; grid spacing for FFT
    ; Temperature coupling is on
    tc-grps = System ; two coupling groups - more accurate
    tau_t = 0.2 ; time constant, in ps
    ref_t = 300 ; reference temperature, one for each group, in K
    ; Pressure coupling is off
    pcoupl = no ; no pressure coupling in NVT
    ; Periodic boundary conditions
    pbc = xyz ; 3-D PBC
    ; Velocity generation
    gen_vel = yes ; assign velocities from Maxwell distribution
    gen_temp = 300 ; temperature for Maxwell distribution
    gen_seed = -1 ; generate a random seed
2

I don’t see a reason for using Langevin Dynamics. Having said that, both LD and normal MD should not results in such large temperature deviations. I don’t see anything strange in your MD settings, so I have no clue what could be the cause of this.

3

Thank you for your response. Could you expand on why there is no reason for using LD here?