[gmx-users] system cooling down when runing NVE

[Thanh Binh NGUYEN]

   Dear Gromacs experts,
> I'm just a newbaby in Gromacs, and hence I have a lot of problem  
> when  running this program. I try to run NVE simulation of a  
> protein. First,  I run an NTP ensemble, follow by NVT and finally,  
> NVE. In NPT and NVT, T remains as constant, however, in NVE   
> simulation, the temperature drops gradually instead of maintain   
> approximately constant. Could you give me any advice to solve this   
> problem?
> Thank you very much.
> Regards,
> Nguyen
>
> P.S: below is my mdp file for three simulations:
> NPT ensemble:
> title = DEN2pro_MD
> define = -DFLEXIBLE
> constraints = all-bonds
> integrator = md
> dt = 0.002 ; ps !
> nsteps = 1000000 ; total 2000 ps.
> nstcomm = 1
> nstxout = 5000 ; output coordinates every 10 ps
> nstvout = 5000 ; output velocities every 10 ps
> nstfout = 0
> nstlog = 10
> nstenergy = 10
> nstlist = 10
> ns_type = grid
> rlist = 0.9
> coulombtype = PME
> rcoulomb = 0.9
> rvdw = 1.0
> fourierspacing = 0.12
> fourier_nx = 0
> fourier_ny = 0
> fourier_nz = 0
> pme_order = 6
> ewald_rtol = 1e-5
> optimize_fft = yes
> ; Berendsen temperature coupling is on in four groups
> Tcoupl = berendsen
> tau_t = 0.1 0.1
> tc_grps = protein non-protein
> ref_t = 300 300
> ; Pressure coupling is on
> Pcoupl = berendsen
> pcoupltype = isotropic
> tau_p = 1.0
> compressibility = 4.5e-5
> ref_p = 1.0
> ; Generate velocites is on at 300 K.
> gen_vel = yes
> gen_temp = 300.0
> gen_seed = 173529
>
>
>
> NVT ensemble
> title       = DEN2pro_MD NVT equilibration
> define      = -DPOSRES  ; position restrain the protein and ligand
> ; Run parameters
> integrator  = md        ; leap-frog integrator
> nsteps      = 100000     ; 2 * 100000 = 200 ps
> dt          = 0.002     ; 2 fs
> ; Output control
> nstxout     = 5000       ; save coordinates every 10 ps
> nstvout     = 5000       ; save velocities every 10 ps
> nstenergy   = 100       ; save energies every 0.2 ps
> nstlog      = 100       ; update log file every 0.2 ps
> energygrps  = Protein
> ; Bond parameters
> continuation    = no            ; first dynamics run
> constraint_algorithm = lincs    ; holonomic constraints
> constraints     = all-bonds     ; all bonds (even heavy atom-H  
> bonds)  constrained
> lincs_iter      = 1             ; accuracy of LINCS
> lincs_order     = 4             ; also related to accuracy
> ; Neighborsearching
> ns_type     = grid      ; search neighboring grid cells
> nstlist     = 5         ; 10 fs
> rlist       = 0.9       ; short-range neighborlist cutoff (in nm)
> rcoulomb    = 0.9       ; short-range electrostatic cutoff (in nm)
> rvdw        = 1.4       ; short-range van der Waals cutoff (in nm)
> ; Electrostatics
> coulombtype     = PME       ; Particle Mesh Ewald for long-range   
> electrostatics
> pme_order       = 4         ; cubic interpolation
> fourierspacing  = 0.12      ; grid spacing for FFT
> ; Temperature coupling is on
> tcoupl      = V-rescale                     ; modified Berendsen thermostat
> tc-grps     = Protein Non-protein	    ; two coupling groups - more accurate
> tau_t       = 0.1   0.1                     ; time constant, in ps
> ref_t       = 300   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
> ; Dispersion correction
> DispCorr    = EnerPres  ; account for cut-off vdW scheme
> ; Velocity generation
> gen_vel     = yes       ; assign velocities from Maxwell distribution
> gen_temp    = 300       ; temperature for Maxwell distribution
> gen_seed    = -1        ; generate a random seed
>
>
> NVE ensemble:
> title       = DEN2pro_MD NVE equilibration
> define      = -DPOSRES  ; position restrain the protein and ligand
> ; Run parameters
> integrator  = md-vv     ; leap-frog integrator
> nsteps      = 1000000   ; 2 * 1000000 = 2000 ps
> dt          = 0.002     ; 2 fs
> ; Output control
> nstxout     = 5000      ; save coordinates every 10 ps
> nstvout     = 5000      ; save velocities every 10 ps
> nstenergy   = 100       ; save energies every 0.2 ps
> nstlog      = 100       ; update log file every 0.2 ps
> energygrps  = Protein
> ; Bond parameters
> continuation    = no            ; first dynamics run
> constraint_algorithm = lincs    ; holonomic constraints
> constraints     = all-bonds     ; all bonds (even heavy atom-H  
> bonds)  constrained
> lincs_iter      = 1             ; accuracy of LINCS
> lincs_order     = 4             ; also related to accuracy
> ; Neighborsearching
> ns_type     = grid      ; search neighboring grid cells
> nstlist     = 5         ; 10 fs
> rlist       = 1.6       ; short-range neighborlist cutoff (in nm)
> rcoulomb    = 0.9       ; short-range electrostatic cutoff (in nm)
> rlistlong   = 1.6
> ;treatment of van der waals interactions
> vdwtype = Shift
> rvdw 	= 0.95 		; short-range van der Waals cutoff (in nm)
> rvdw-switch = 0.9
> ; Electrostatics
> coulombtype     = PME-Switch; Particle Mesh Ewald for long-range   
> electrostatics
> pme_order       = 4         ; cubic interpolation
> fourierspacing  = 0.12      ; grid spacing for FFT
> ; Temperature coupling is on
> tcoupl      = no                     	    ; no temperature coupling in NVE
> ; Pressure coupling is off
> pcoupl      = no        ; no pressure coupling in NVE
> ; Periodic boundary conditions
> pbc         = xyz       ; 3-D PBC
> ; Dispersion correction
> DispCorr    = EnerPres  ; account for cut-off vdW scheme
> ; Velocity generation
> gen_vel     = yes       ; assign velocities from Maxwell distribution
> gen_temp    = 300       ; temperature for Maxwell distribution
> gen_seed    = -1        ; generate a random seed
>
>