Dear Paul 

You can control the temperature in gromacs. Just use the following lines in .mdp file. 

You can adjust parameters in your case.

; SIMULATED ANNEALING

; Type of annealing for each temperature group (no/single/periodic)

annealing = single

; Number of time points to use for specifying annealing in each group

annealing-npoints = 5

; List of times at the annealing points for each group

annealing-time =0 5000 10000 15000 20000

; Temp. at each annealing point, for each group.

annealing-temp = 200 250 300 350 400

An immediate increase from 0 -> 300 K suggests you did not appropriately generate velocities at the outset of the run and should be expected with ref-t = 300. The thermostat will try to maintain the temperature you specify, which is precisely what it is doing. To equilibrate at a given temperature, you should generate velocities appropriately. But if your goal is truly is heating, you must use simulated annealing, as suggested above. There is no other way to gradually heat the system.

Hi,

I want only heating of the protein. Actually I tried simulated annealing  with

annealing = single

annealing_npoints = 2

annealing_time = 0 10

annealing_temp = 0 300

but  temperature  never reaches to 300 K . It  increases upto 100 K  only . Should I use higher temperature upper limit ? Why does my system never reach temperature value as mentioned in the .mdp file ?

Provide a full .mdp file. Snippets are not enough to tell what is going on. Heating only the protein is extremely artificial; what are you trying to accomplish by doing so? If you really want to do this, you have to use multiple tc-grps and multiple heating series in the simulated annealing setup (you can keep groups at constant temperature this way while others change). But this sort of temperature gradient is completely non-physical.

Hi,

  My main goal is to equilibriate the protein structure in step by step manner and finally to run for a long time in NVT condition. I took the crystal structure , added solvent and ions to neutralize the system and then optimized the structure keeping all the protein heavy atom fixed without any temperature and pressure coupling (we call this step 'min-prot-fixed'). Then I heated the system from 0 K to 300 K keeping protein heavy atoms fixed ('heat-prot-fixed').After this we relaxed the constrained by constraining all protein heavy atoms by 25 kcal/mol/angstrom^2('min-prot-25') and again minimized the system('min-prot-25') followed by heating from 0 K to 300 K without T and P coupling ('heat-prot-25').Then equilibriation was done in NPT condition with temperature and pressure coupling('NPT-25') followed by gradually relaxing the constraints to 10('NPT-10') then 5('NPT-5') and completely free('NPT-free').Then I ran it for a longer time like 100 ns in NVT condition.

In the step heat-prot-fixed and heat-prot-25 I need to heat the system from 0 K to 300 K .I tried with two methods as I mentioned before producing a problem in each of the cases.In the case where I used simulated annealing the problem is that temperature never reaches to 300 K rather it stops its' increament at 100 K .Using v-rescale coupling temperature starts at 0 K and reaches to 300 K and stabilizes but the increament happen very rapidly.I want a slow increament of temperature and  extent of this increament should be modulated by user.I used gen-vel = 0 to start temperature at 0 K.If I use gen-vel = 300 K which is suggested by many of the tutorials starting temperature become 300 K which I don't want.Here is the full .mdp file for the step heat-prot-fixed step with simulated annealing

title = Ubiquitin system heat protein fixed equilibration

; Run parameters

integrator = md ; leap-frog integrator

nsteps = 10000 ; 2 * 50000 = 100 ps

dt = 0.002 ; 2 fs

freezegrps = protein-noh

freezedim = Y Y Y

; Output control

nstxout = 500 ; save coordinates every 1.0 ps

nstvout = 500 ; save velocities every 1.0 ps

nstenergy = 1 ; save energies every 1/500 ps

nstlog = 500 ; update log file every 1.0 ps

; 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

cutoff-scheme = Verlet

ns_type = grid ; search neighboring grid cells

nstlist = 10 ; 20 fs, largely irrelevant with Verlet

rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)

rvdw = 1.0 ; 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.16 ; grid spacing for FFT

; Temperature coupling is on

annealing = single

annealing_npoints = 2

annealing_time = 0 10

annealing_temp = 0 300

; 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 = 0 ; temperature for Maxwell distribution

gen_seed = -1 ; generate a random seed

Here is the full .mdp file for the step heat-prot-fixed step with v-rescale tcouple

title = Ubiquitin system heat protein fixed equilibration

; Run parameters

integrator = md ; leap-frog integrator

nsteps = 10000 ; 2 * 50000 = 100 ps

dt = 0.002 ; 2 fs

freezegrps = protein-noh

freezedim = Y Y Y

; Output control

nstxout = 500 ; save coordinates every 1.0 ps

nstvout = 500 ; save velocities every 1.0 ps

nstenergy = 1 ; save energies every 1/500 ps

nstlog = 500 ; update log file every 1.0 ps

; 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

cutoff-scheme = Verlet

ns_type = grid ; search neighboring grid cells

nstlist = 10 ; 20 fs, largely irrelevant with Verlet

rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)

rvdw = 1.0 ; 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.16 ; grid spacing for FFT

; Temperature coupling is on

tcoupl = v-rescale ; modified Berendsen thermostat

tc-grps = all ; two coupling groups - more accurate

tau_t = 0.4 ; 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

; Dispersion correction

DispCorr = EnerPres ; account for cut-off vdW scheme

; Velocity generationgen_vel = yes ; assign velocities from Maxwell distribution

gen_temp = 0 ; temperature for Maxwell distribution

gen_seed = -1 ; generate a random seed

Your preparation protocol is far more complex than is necessary to reasonably prepare the system. But do as you see fit. You're freezing the protein, which means it cannot have velocities, so the output temperature is some artificial representation of a system that is being annealed and otherwise perturbed. You're doing more harm than good here. Freezing is a severe perturbation, and it actually inhibits equilibration rather than helping it. You'll start in some quasi-equilibrated, non-equilibrium state in which the protein is frozen (literally, since its velocities are zero) in a comparatively warm bath of solvent. This is completely unphysical.

A normal protocol might be:

1. Energy minimization (restraints here again inhibit your progress)

2. Heating from 0-300 K under NVT conditions, with restraints on the protein (not freezegrps)

3. Restrained equilibration at 300 K under an NPT ensemble, conducted in several phases if you wish to step down the restraints

4. Production

Dear Justin,

I follwed the protocol you suggested above but the same problem is generating i.e. temperature is not increasing to 300 K at heating step.I performed energy minimization and then started heating the system from 0-300 K under NVT conditions with 25 kcal/mol/A^2 restraints on all protein heavy atoms.After performing this step I found that temperature never reach 300 K , it increases upto 100 K only.Here is the full .mdp file

title = ubiquitin heating under NVT equilibration

define = -DPOSRES ; position restrain the protein

; Run parameters

integrator = md ; leap-frog integrator

nsteps = 10000 ; 2 * 10000 = 20000 fs = 20ps

dt = 0.002 ; 2 fs

; Output control

nstxout = 500 ; save coordinates every 1.0 ps

nstvout = 500 ; save velocities every 1.0 ps

nstenergy = 1 ; save energies every 1.0 ps

nstlog = 1 ; update log file every 1.0 ps

; 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

cutoff-scheme = Verlet

ns_type = grid ; search neighboring grid cells

nstlist = 10 ; 20 fs, largely irrelevant with Verlet

rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)

rvdw = 1.0 ; 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.16 ; grid spacing for FFT

; Temperature coupling is on

annealing = single

annealing_npoints = 5

annealing_time = 0 4 8 12 20

annealing_temp = 0 75 150 300 300

; 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 = 0 ; temperature for Maxwell distribution

gen_seed = -1 ; generate a random seed

Your .mdp file does not set tcoupl, so by default it is off. You need to use some thermostat in order to regulate temperature. I'm surprised the run even worked without it.

Dear Justin ,

I modified the .mdp file and here is the full modified .mdp file

title = ubiquitin heating under NVT equilibration

define = -DPOSRES ; position restrain the protein

; Run parameters

integrator = md ; leap-frog integrator

nsteps = 10000 ; 2 * 50000000 = 100 ns

dt = 0.002 ; 2 fs

; Output control

nstxout = 500 ; save coordinates every 1.0 ps

nstvout = 500 ; save velocities every 1.0 ps

nstenergy = 1 ; save energies every 1.0 ps

nstlog = 1 ; update log file every 1.0 ps

; 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

cutoff-scheme = Verlet

ns_type = grid ; search neighboring grid cells

nstlist = 10 ; 20 fs, largely irrelevant with Verlet

rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)

rvdw = 1.0 ; 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.16 ; grid spacing for FFT

; Temperature coupling is on

tcoupl = v-rescale ; modified Berendsen thermostat

tc-grps = system

tau_t = 0.4 ; time constant, in ps

ref_t = 300 ; reference temperature, one for each group, in K

annealing = single

annealing_npoints = 5

annealing_time = 0 4 8 12 20

annealing_temp = 0 75 150 300 300

; 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 = 0 ; temperature for Maxwell distribution

gen_seed = -1 ; generate a random seed

Now temperature increases from 0 K to 300 K but I found some abnormal bumps and discontinuity  in the temperature vs. time plot. I have attached one plot for the above mentioned .mdp file. I changed tau_t from 0.4 to 0.1 but still the same problem was coming.

No idea why that would happen. You're not stopping and re-starting the job in any way, are you?  The sudden spikes in temperature suggest discontinuities.

It is a continuous 20 ps simulation and I am not interrupting in between.I am saving energy,temperature etc. at every timestep i.e. at the interval of 2 fs (in normal tutorial it was suggested that save energy and all the output after a picosecond interval and if I do that my picture will be clear removing all the bumps and discontinuities).I have attached temperature plot for another protein system (apo-azurin , pdb id = 1E65) and I found same abnormalities.I have attached one zoomed in figure of the discontinuous region and you will find that for a single time I got multiple temperature value at that region.Did you find bumps and discontinuity ever in your heating step with your protein if you save temperature at each of the timestep ?

title = Simulated Annealing

define = -DPOSRES -DPOSRES_LIPID ; restrain protein and lipid P8

; Run parameters

integrator = md ; leap-frog integrator

dt = 0.002 ; 2 fs

nsteps = 250000 ; 500 ps

; Bond parameters

continuation = no ; starting up

constraints = all-bonds ; constrain all bond lengths

constraint-algorithm = lincs ; holonomic constraints

lincs-iter = 1 ; accuracy of LINCS

lincs-order = 4 ; also related to accuracy

; Output control

nstxout = 1000 ; save coordinates every 2 ps

nstvout = 1000 ; save velocities every 2 ps

nstfout = 1000 ; save forces every 2 ps

nstenergy = 1000 ; save energies every 2 ps

; Neighborsearching

nstlist = 5 ; 10 fs

ns_type = grid ; search neighboring grid cells

rlist = 1.2 ; short-range neighborlist cutoff (nm)

rcoulomb = 1.2 ; short-range electrostatic cutoff (nm)

rvdw = 1.2 ; short-range van der Waals cutoff (nm)

; 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 in three groups

Tcoupl = Berendsen ; Weak coupling

tc_grps = Protein DPPC SOL_CL ; three coupling groups - more accurate

tau_t = 0.1 0.1 0.1 ; time constant, in ps

ref_t = 323 323 323 ; reference temperature, one for each group, in K

; Pressure coupling

Pcoupl = Berendsen ; Weak coupling

Pcoupltype = semiisotropic ; uniform scaling of x-y vectors, independent z

ref_p = 1.0 1.0 ; reference pressure, x-y, z (in bar)

compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility, bar^-1

; Generate velocities is off

gen_vel = no ; Velocity generation is off

; Periodic boundary conditions are on in all directions

pbc = xyz ; 3-D PBC

; Dispersion correction

DispCorr = EnerPres ; account for cut-off vdW scheme

; COM motion removal

; These options remove motion of the protein/bilayer relative to the solvent/ions

nstcomm = 1

comm-mode = Linear

comm-grps = Protein_DPPC SOL_CL

; Simulated annealing

annealing = single single single ; single sequence of points for each T-coupling group

annealing_npoints = 2 2 2 ; two points - start and end temperatures

annealing_time = 0 500 0 500 0 500 ; time frame of heating - heat over period of 500 ps

annealing_temp = 0 323 0 323 0 323 ; start and end temperatures

Hi Justin Lemkul , I have a system in which I have a fixed silicon substrate on which I am placing a water cube. I have been able to equilibrate the water to give contact angles and make it roll on the surface using a gradient of chemical potential on the surface by adjusting lennard-jones potentials values.

I am currently trying to evaporate a droplet while equilibration is occurring, but I realised through gmx energy commands that the substrate temperature, though coupled with v-rescale and coupling constant as 0 (very strong), the temperature of my substrate has always remained 0K in all of my previous simulations. I think it is because of freezegroups on the substrate that I have imposed in all 3 dimensions. How do I facilitate the transfer of heat from the fixed solid substrate to the droplet? How do I actually have heating of the substrate? And again, how should I successfully execute a system in which I have a droplet which is un-coupled but placed on a substrate which is getting heated from temperature T1 to T2 to have evaporation?

Thanks,

Raj

Raj Narayanan I'm not sure you can even do that. The use of freezegrps means every collision is inelastic and the substrate never has any velocities, thus no kinetic energy or temperature. You certainly can't freeze this entity. It sounds to me like what you want is to turn temperature coupling off for this group (tau_t = -1) but I've never attempted anything like what you're trying. You may find better features in software like LAMMPS which is more geared towards the materials community.

Thanks a lot @Justin Lemkul.

Hi Justin Lemkul ,

Could I use position restraints instead of a freezegrps? Would it allow me to give a temperature to a substrate?