Re: [gmx-users] Different temperatures for different groups, even with Nose-Hoover
The comment from Justin and your answer to it were pretty convincing. I think the problem is solved :)). Tom is correct to get T with g_energy on a simulation that ran with one Tcoup, you need to rerun the trajectory (trr) with a tpr in which you define two Tcoup. This should give you T for each group. XAvier. On Oct 29, 2009, at 11:11 PM, Michael Lerner wrote: On Thu, Oct 29, 2009 at 5:31 PM, XAvier Periole x.peri...@rug.nl wrote: On Oct 29, 2009, at 9:26 PM, Michael Lerner wrote: I calculated temperatures with g_traj -f md2.trr -s md2.tpr -n index.ndx -ot -ng 5 I did not know you could get the temperature throught g_traj ... the -ot option is the one giving you the temperature? Yes. Do you get the same temperature difference with g_energy? Or may be the info is not in there? May be you could rerun you trajectory using a new tpr where you define two different temp coupling groups. And see if you observe the same temperature given by g_energy. It's possible that I'm doing something wrong, but I couldn't convince g_energy to give me the temperature of different groups when I used a single thermostat for the entire system. If I use two temperature coupling groups, I get the expected results. As Justin pointed out, it looks like the error was mine for not understanding the details of g_traj. When I manually correct for the degrees of freedom, I get temperatures for Protein/Sol that are close enough to each other (i.e. within a couple of standard deviations). Thanks, -Michael Thanks, -Michael -- begin md2.mdp -- ; RUN CONTROL PARAMETERS = integrator = md ; start time and timestep in ps = tinit= 0.0 dt = 0.020 nsteps = 2550 ; mode for center of mass motion removal comm-mode= linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm-grps= System ; OUTPUT CONTROL OPTIONS = ; Output frequency for coords (x), velocities (v) and forces (f) = nstxout = 5 nstvout = 5 nstfout = 0 ; Output frequency for energies to log file and energy file = nstlog = 2500 nstenergy= 2500 ; Output frequency and precision for xtc file = nstxtcout= 2500 xtc_precision= 100 ; This selects the subset of atoms for the xtc file. You can = ; select multiple groups. By default all atoms will be written. = xtc-grps = ; Selection of energy groups = energygrps = System Protein W WF ; NEIGHBORSEARCHING PARAMETERS = ; nblist update frequency = nstlist = 5 ; ns algorithm (simple or grid) = ns_type = grid ; Periodic boundary conditions: xyz or none = pbc = xyz ; nblist cut-off = rlist= 1.2 ; OPTIONS FOR ELECTROSTATICS AND VDW = ; Method for doing electrostatics = coulombtype = Shift rcoulomb_switch = 0.0 rcoulomb = 1.2 ; Dielectric constant (DC) for cut-off or DC of reaction field = epsilon_r= 15 ; Method for doing Van der Waals = vdw_type = Shift ; cut-off lengths= rvdw_switch = 0.9 rvdw = 1.2 ; Apply long range dispersion corrections for Energy and Pressure = DispCorr = No ; Spacing for the PME/PPPM FFT grid = fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used = fourier_nx = 10 fourier_ny = 10 fourier_nz = 10 ; EWALD/PME/PPPM parameters = pme_order= 4 ewald_rtol = 1e-05 epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS = ; Temperature coupling = tcoupl = Nose-Hoover ; Groups to couple separately = tc-grps = System ; Time constant (ps) and reference temperature (K) = tau_t= 2 ref_t= 323 ; Pressure coupling = Pcoupl = Parrinello-Rahman Pcoupltype = isotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) = tau_p= 5 compressibility = 5e-5 ref_p= 1.0 ; SIMULATED ANNEALING CONTROL = annealing= no ; GENERATE VELOCITIES FOR STARTUP RUN = continuation = yes gen_vel = no ;gen_temp = 323 ;gen_seed = 473529 ; OPTIONS FOR BONDS = constraints = none ; Type of constraint algorithm = constraint_algorithm = Lincs ; Relative tolerance of shake = shake_tol= 0.0001 ; Highest order in the expansion of the constraint coupling matrix = lincs_order = 4 ; Lincs will write a warning to the stderr if in one
[gmx-users] Different temperatures for different groups, even with Nose-Hoover
Hi, I'm using GROMACS 4.0.5 to do a MARTINI simulation of a CG protein in a CG water box, looking at the diffusion constant of the protein among other things. I'm using the Nose-Hoover thermostat and Parrinello-Rahman barostat with tc-grps = System. In CHARMM, a protein-in-water simulation with an extended ensemble will show equal temperatures for the protein and the water. However, I'm finding that, with ref_t = 323, I get a SOL (both W and WF particles) temperature of 323 and a Protein temperature of 295. Is this what I should expect? The system has 62.5k particles, the box is 20 x 20 x 20, and the relevant .mdp parameters I've used are: tcoupl = Nose-Hoover tc-grps = System tau_t= 2 ref_t= 323 Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p= 5 compressibility = 5e-5 ref_p= 1.0 These parameters seem to lead to stable systems for bilayer simulations. I've included the whole .mdp file at the end in case I've forgotten to mention something relevant. I calculated temperatures with g_traj -f md2.trr -s md2.tpr -n index.ndx -ot -ng 5 Thanks, -Michael -- begin md2.mdp -- ; RUN CONTROL PARAMETERS = integrator = md ; start time and timestep in ps = tinit= 0.0 dt = 0.020 nsteps = 2550 ; mode for center of mass motion removal comm-mode= linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm-grps= System ; OUTPUT CONTROL OPTIONS = ; Output frequency for coords (x), velocities (v) and forces (f) = nstxout = 5 nstvout = 5 nstfout = 0 ; Output frequency for energies to log file and energy file = nstlog = 2500 nstenergy= 2500 ; Output frequency and precision for xtc file = nstxtcout= 2500 xtc_precision= 100 ; This selects the subset of atoms for the xtc file. You can = ; select multiple groups. By default all atoms will be written. = xtc-grps = ; Selection of energy groups = energygrps = System Protein W WF ; NEIGHBORSEARCHING PARAMETERS = ; nblist update frequency = nstlist = 5 ; ns algorithm (simple or grid) = ns_type = grid ; Periodic boundary conditions: xyz or none = pbc = xyz ; nblist cut-off = rlist= 1.2 ; OPTIONS FOR ELECTROSTATICS AND VDW = ; Method for doing electrostatics = coulombtype = Shift rcoulomb_switch = 0.0 rcoulomb = 1.2 ; Dielectric constant (DC) for cut-off or DC of reaction field = epsilon_r= 15 ; Method for doing Van der Waals = vdw_type = Shift ; cut-off lengths= rvdw_switch = 0.9 rvdw = 1.2 ; Apply long range dispersion corrections for Energy and Pressure = DispCorr = No ; Spacing for the PME/PPPM FFT grid = fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used = fourier_nx = 10 fourier_ny = 10 fourier_nz = 10 ; EWALD/PME/PPPM parameters = pme_order= 4 ewald_rtol = 1e-05 epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS = ; Temperature coupling = tcoupl = Nose-Hoover ; Groups to couple separately = tc-grps = System ; Time constant (ps) and reference temperature (K) = tau_t= 2 ref_t= 323 ; Pressure coupling = Pcoupl = Parrinello-Rahman Pcoupltype = isotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) = tau_p= 5 compressibility = 5e-5 ref_p= 1.0 ; SIMULATED ANNEALING CONTROL = annealing= no ; GENERATE VELOCITIES FOR STARTUP RUN = continuation = yes gen_vel = no ;gen_temp = 323 ;gen_seed = 473529 ; OPTIONS FOR BONDS = constraints = none ; Type of constraint algorithm = constraint_algorithm = Lincs ; Relative tolerance of shake = shake_tol= 0.0001 ; Highest order in the expansion of the constraint coupling matrix = lincs_order = 4 ; Lincs will write a warning to the stderr if in one step a bond = ; rotates over more degrees than = lincs_warnangle = 30 ; Convert harmonic bonds to morse potentials = morse= no ; NMR refinement stuff = ; Distance restraints type: No, Simple or Ensemble = disre= No ; Force weighting of pairs in one distance restraint: Equal or
Re: [gmx-users] Different temperatures for different groups, even with Nose-Hoover
On Oct 29, 2009, at 9:26 PM, Michael Lerner wrote: Hi, I'm using GROMACS 4.0.5 to do a MARTINI simulation of a CG protein in a CG water box, looking at the diffusion constant of the protein among other things. I'm using the Nose-Hoover thermostat and Parrinello-Rahman barostat with tc-grps = System. In CHARMM, a protein-in-water simulation with an extended ensemble will show equal temperatures for the protein and the water. However, I'm finding that, with ref_t = 323, I get a SOL (both W and WF particles) temperature of 323 and a Protein temperature of 295. Is this what I should expect? No. The system has 62.5k particles, the box is 20 x 20 x 20, and the relevant .mdp parameters I've used are: tcoupl = Nose-Hoover tc-grps = System tau_t= 2 ref_t= 323 Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p= 5 compressibility = 5e-5 ref_p= 1.0 These parameters seem to lead to stable systems for bilayer simulations. I've included the whole .mdp file at the end in case I've forgotten to mention something relevant. I calculated temperatures with g_traj -f md2.trr -s md2.tpr -n index.ndx -ot -ng 5 I did not know you could get the temperature throught g_traj ... the - ot option is the one giving you the temperature? Do you get the same temperature difference with g_energy? Or may be the info is not in there? May be you could rerun you trajectory using a new tpr where you define two different temp coupling groups. And see if you observe the same temperature given by g_energy. Thanks, -Michael -- begin md2.mdp -- ; RUN CONTROL PARAMETERS = integrator = md ; start time and timestep in ps = tinit= 0.0 dt = 0.020 nsteps = 2550 ; mode for center of mass motion removal comm-mode= linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm-grps= System ; OUTPUT CONTROL OPTIONS = ; Output frequency for coords (x), velocities (v) and forces (f) = nstxout = 5 nstvout = 5 nstfout = 0 ; Output frequency for energies to log file and energy file = nstlog = 2500 nstenergy= 2500 ; Output frequency and precision for xtc file = nstxtcout= 2500 xtc_precision= 100 ; This selects the subset of atoms for the xtc file. You can = ; select multiple groups. By default all atoms will be written. = xtc-grps = ; Selection of energy groups = energygrps = System Protein W WF ; NEIGHBORSEARCHING PARAMETERS = ; nblist update frequency = nstlist = 5 ; ns algorithm (simple or grid) = ns_type = grid ; Periodic boundary conditions: xyz or none = pbc = xyz ; nblist cut-off = rlist= 1.2 ; OPTIONS FOR ELECTROSTATICS AND VDW = ; Method for doing electrostatics = coulombtype = Shift rcoulomb_switch = 0.0 rcoulomb = 1.2 ; Dielectric constant (DC) for cut-off or DC of reaction field = epsilon_r= 15 ; Method for doing Van der Waals = vdw_type = Shift ; cut-off lengths= rvdw_switch = 0.9 rvdw = 1.2 ; Apply long range dispersion corrections for Energy and Pressure = DispCorr = No ; Spacing for the PME/PPPM FFT grid = fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used = fourier_nx = 10 fourier_ny = 10 fourier_nz = 10 ; EWALD/PME/PPPM parameters = pme_order= 4 ewald_rtol = 1e-05 epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS = ; Temperature coupling = tcoupl = Nose-Hoover ; Groups to couple separately = tc-grps = System ; Time constant (ps) and reference temperature (K) = tau_t= 2 ref_t= 323 ; Pressure coupling = Pcoupl = Parrinello-Rahman Pcoupltype = isotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) = tau_p= 5 compressibility = 5e-5 ref_p= 1.0 ; SIMULATED ANNEALING CONTROL = annealing= no ; GENERATE VELOCITIES FOR STARTUP RUN = continuation = yes gen_vel = no ;gen_temp = 323 ;gen_seed = 473529 ; OPTIONS FOR BONDS = constraints = none ; Type of constraint algorithm = constraint_algorithm = Lincs ; Relative tolerance of shake = shake_tol= 0.0001 ;
Re: [gmx-users] Different temperatures for different groups, even with Nose-Hoover
Michael Lerner wrote: Hi, I'm using GROMACS 4.0.5 to do a MARTINI simulation of a CG protein in a CG water box, looking at the diffusion constant of the protein among other things. I'm using the Nose-Hoover thermostat and Parrinello-Rahman barostat with tc-grps = System. In CHARMM, a protein-in-water simulation with an extended ensemble will show equal temperatures for the protein and the water. However, I'm finding that, with ref_t = 323, I get a SOL (both W and WF particles) temperature of 323 and a Protein temperature of 295. Is this what I should expect? The system has 62.5k particles, the box is 20 x 20 x 20, and the relevant .mdp parameters I've used are: tcoupl = Nose-Hoover tc-grps = System tau_t= 2 ref_t= 323 Pcoupl = Parrinello-Rahman Pcoupltype = isotropic tau_p= 5 compressibility = 5e-5 ref_p= 1.0 These parameters seem to lead to stable systems for bilayer simulations. I've included the whole .mdp file at the end in case I've forgotten to mention something relevant. I calculated temperatures with g_traj -f md2.trr -s md2.tpr -n index.ndx -ot -ng 5 I see you are using constraints = none, but are there any constraints defined in the topology? If so, as noted in g_traj -h: Option -ot plots the temperature of each group, provided velocities are present in the trajectory file. No corrections are made for constrained degrees of freedom! This implies -com. So you will get an inherently incorrect answer when using g_traj, if there are any constraints. -Justin Thanks, -Michael -- begin md2.mdp -- ; RUN CONTROL PARAMETERS = integrator = md ; start time and timestep in ps = tinit= 0.0 dt = 0.020 nsteps = 2550 ; mode for center of mass motion removal comm-mode= linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm-grps= System ; OUTPUT CONTROL OPTIONS = ; Output frequency for coords (x), velocities (v) and forces (f) = nstxout = 5 nstvout = 5 nstfout = 0 ; Output frequency for energies to log file and energy file = nstlog = 2500 nstenergy= 2500 ; Output frequency and precision for xtc file = nstxtcout= 2500 xtc_precision= 100 ; This selects the subset of atoms for the xtc file. You can = ; select multiple groups. By default all atoms will be written. = xtc-grps = ; Selection of energy groups = energygrps = System Protein W WF ; NEIGHBORSEARCHING PARAMETERS = ; nblist update frequency = nstlist = 5 ; ns algorithm (simple or grid) = ns_type = grid ; Periodic boundary conditions: xyz or none = pbc = xyz ; nblist cut-off = rlist= 1.2 ; OPTIONS FOR ELECTROSTATICS AND VDW = ; Method for doing electrostatics = coulombtype = Shift rcoulomb_switch = 0.0 rcoulomb = 1.2 ; Dielectric constant (DC) for cut-off or DC of reaction field = epsilon_r= 15 ; Method for doing Van der Waals = vdw_type = Shift ; cut-off lengths= rvdw_switch = 0.9 rvdw = 1.2 ; Apply long range dispersion corrections for Energy and Pressure = DispCorr = No ; Spacing for the PME/PPPM FFT grid = fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used = fourier_nx = 10 fourier_ny = 10 fourier_nz = 10 ; EWALD/PME/PPPM parameters = pme_order= 4 ewald_rtol = 1e-05 epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS = ; Temperature coupling = tcoupl = Nose-Hoover ; Groups to couple separately = tc-grps = System ; Time constant (ps) and reference temperature (K) = tau_t= 2 ref_t= 323 ; Pressure coupling = Pcoupl = Parrinello-Rahman Pcoupltype = isotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) = tau_p= 5 compressibility = 5e-5 ref_p= 1.0 ; SIMULATED ANNEALING CONTROL = annealing= no ; GENERATE VELOCITIES FOR STARTUP RUN = continuation = yes gen_vel = no ;gen_temp = 323 ;gen_seed = 473529 ; OPTIONS FOR BONDS = constraints = none ; Type of constraint algorithm = constraint_algorithm = Lincs ; Relative tolerance of shake = shake_tol= 0.0001 ; Highest order in
Re: [gmx-users] Different temperatures for different groups, even with Nose-Hoover
On Thu, Oct 29, 2009 at 5:38 PM, Justin A. Lemkul jalem...@vt.edu wrote:
I see you are using constraints = none, but are there any constraints
defined in the topology? If so, as noted in g_traj -h:
Option -ot plots the temperature of each group, provided velocities are
present in the trajectory file. No corrections are made for constrained
degrees of freedom! This implies -com.
So you will get an inherently incorrect answer when using g_traj, if there
are any constraints.
Oops! I actually had read some previous messages about that, but I got
confused and thought the problem was with g_energy, not g_traj. I'm not
particularly familiar with the GROMACS source code, but gmx_traj.c has this
function:
static real temp(rvec v[],real mass[],int isize,atom_id index[])
{
real ekin2=0;
int i,j;
for(i=0; iisize; i++) {
j = index[i];
ekin2 += mass[j]*norm2(v[j]);
}
return ekin2/(3*isize*BOLTZ);
}
so it looks like g_traj is just assuming that there are 3*N degrees of
freedom. My protein topology file has 304 particles and 83 lines in the
constraints section, so I think the correction should be
295.3496*((3*304)/(3*304 - 83)) = 324.9202
for my protein vs. 322.9525 for the W's and 322.9279 for the WF's.
So, they're still a little off, but doing some block averaging and looking
at standard errors, it's well within the acceptable range.
Thanks,
-Michael
-Justin
Thanks,
-Michael
-- begin md2.mdp --
; RUN CONTROL PARAMETERS =
integrator = md
; start time and timestep in ps =
tinit= 0.0
dt = 0.020
nsteps = 2550
; mode for center of mass motion removal
comm-mode= linear
; number of steps for center of mass motion removal
nstcomm = 1
; group(s) for center of mass motion removal
comm-grps= System
; OUTPUT CONTROL OPTIONS =
; Output frequency for coords (x), velocities (v) and forces (f) =
nstxout = 5
nstvout = 5
nstfout = 0
; Output frequency for energies to log file and energy file =
nstlog = 2500
nstenergy= 2500
; Output frequency and precision for xtc file =
nstxtcout= 2500
xtc_precision= 100
; This selects the subset of atoms for the xtc file. You can =
; select multiple groups. By default all atoms will be written. =
xtc-grps =
; Selection of energy groups =
energygrps = System Protein W WF
; NEIGHBORSEARCHING PARAMETERS =
; nblist update frequency =
nstlist = 5
; ns algorithm (simple or grid) =
ns_type = grid
; Periodic boundary conditions: xyz or none =
pbc = xyz
; nblist cut-off =
rlist= 1.2
; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype = Shift
rcoulomb_switch = 0.0
rcoulomb = 1.2
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon_r= 15
; Method for doing Van der Waals =
vdw_type = Shift
; cut-off lengths=
rvdw_switch = 0.9
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr = No
; Spacing for the PME/PPPM FFT grid =
fourierspacing = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used =
fourier_nx = 10
fourier_ny = 10
fourier_nz = 10
; EWALD/PME/PPPM parameters =
pme_order= 4
ewald_rtol = 1e-05
epsilon_surface = 0
optimize_fft = no
; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; Temperature coupling =
tcoupl = Nose-Hoover
; Groups to couple separately =
tc-grps = System
; Time constant (ps) and reference temperature (K) =
tau_t= 2
ref_t= 323
; Pressure coupling =
Pcoupl = Parrinello-Rahman
Pcoupltype = isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau_p= 5
compressibility = 5e-5
ref_p= 1.0
; SIMULATED ANNEALING CONTROL =
annealing= no
; GENERATE VELOCITIES FOR STARTUP RUN =
continuation = yes
gen_vel = no
;gen_temp = 323
;gen_seed = 473529
; OPTIONS FOR BONDS =
constraints = none
; Type of constraint algorithm =
constraint_algorithm = Lincs
; Relative tolerance of shake =
shake_tol= 0.0001
; Highest order in the expansion of the constraint coupling matrix =
lincs_order = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more
Re: [gmx-users] Different temperatures for different groups, even with Nose-Hoover
On Thu, Oct 29, 2009 at 5:31 PM, XAvier Periole x.peri...@rug.nl wrote: On Oct 29, 2009, at 9:26 PM, Michael Lerner wrote: I calculated temperatures with g_traj -f md2.trr -s md2.tpr -n index.ndx -ot -ng 5 I did not know you could get the temperature throught g_traj ... the -ot option is the one giving you the temperature? Yes. Do you get the same temperature difference with g_energy? Or may be the info is not in there? May be you could rerun you trajectory using a new tpr where you define two different temp coupling groups. And see if you observe the same temperature given by g_energy. It's possible that I'm doing something wrong, but I couldn't convince g_energy to give me the temperature of different groups when I used a single thermostat for the entire system. If I use two temperature coupling groups, I get the expected results. As Justin pointed out, it looks like the error was mine for not understanding the details of g_traj. When I manually correct for the degrees of freedom, I get temperatures for Protein/Sol that are close enough to each other (i.e. within a couple of standard deviations). Thanks, -Michael Thanks, -Michael -- begin md2.mdp -- ; RUN CONTROL PARAMETERS = integrator = md ; start time and timestep in ps = tinit= 0.0 dt = 0.020 nsteps = 2550 ; mode for center of mass motion removal comm-mode= linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm-grps= System ; OUTPUT CONTROL OPTIONS = ; Output frequency for coords (x), velocities (v) and forces (f) = nstxout = 5 nstvout = 5 nstfout = 0 ; Output frequency for energies to log file and energy file = nstlog = 2500 nstenergy= 2500 ; Output frequency and precision for xtc file = nstxtcout= 2500 xtc_precision= 100 ; This selects the subset of atoms for the xtc file. You can = ; select multiple groups. By default all atoms will be written. = xtc-grps = ; Selection of energy groups = energygrps = System Protein W WF ; NEIGHBORSEARCHING PARAMETERS = ; nblist update frequency = nstlist = 5 ; ns algorithm (simple or grid) = ns_type = grid ; Periodic boundary conditions: xyz or none = pbc = xyz ; nblist cut-off = rlist= 1.2 ; OPTIONS FOR ELECTROSTATICS AND VDW = ; Method for doing electrostatics = coulombtype = Shift rcoulomb_switch = 0.0 rcoulomb = 1.2 ; Dielectric constant (DC) for cut-off or DC of reaction field = epsilon_r= 15 ; Method for doing Van der Waals = vdw_type = Shift ; cut-off lengths= rvdw_switch = 0.9 rvdw = 1.2 ; Apply long range dispersion corrections for Energy and Pressure = DispCorr = No ; Spacing for the PME/PPPM FFT grid = fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used = fourier_nx = 10 fourier_ny = 10 fourier_nz = 10 ; EWALD/PME/PPPM parameters = pme_order= 4 ewald_rtol = 1e-05 epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS = ; Temperature coupling = tcoupl = Nose-Hoover ; Groups to couple separately = tc-grps = System ; Time constant (ps) and reference temperature (K) = tau_t= 2 ref_t= 323 ; Pressure coupling = Pcoupl = Parrinello-Rahman Pcoupltype = isotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) = tau_p= 5 compressibility = 5e-5 ref_p= 1.0 ; SIMULATED ANNEALING CONTROL = annealing= no ; GENERATE VELOCITIES FOR STARTUP RUN = continuation = yes gen_vel = no ;gen_temp = 323 ;gen_seed = 473529 ; OPTIONS FOR BONDS = constraints = none ; Type of constraint algorithm = constraint_algorithm = Lincs ; Relative tolerance of shake = shake_tol= 0.0001 ; Highest order in the expansion of the constraint coupling matrix = lincs_order = 4 ; Lincs will write a warning to the stderr if in one step a bond = ; rotates over more degrees than = lincs_warnangle = 30 ; Convert harmonic bonds to morse potentials = morse= no ; NMR refinement stuff = ; Distance restraints type: No, Simple or Ensemble = disre= No ; Force
Re: [gmx-users] Different temperatures for different groups, even with Nose-Hoover
To get g_energy to give you the temperature (or any other property) of a specific group you need to have defined appropriate energygrps in your mdp file. This is independent of the choice of groups you use for the temperature coupling. If you want to calculate these temperatures of the groups after you have already ran your simulations you can use the mdrun -rerun option to save you a lot of time. Tom --On Thursday, October 29, 2009 18:11:53 -0400 Michael Lerner mglerner+grom...@gmail.com wrote: On Thu, Oct 29, 2009 at 5:31 PM, XAvier Periole x.peri...@rug.nl wrote: On Oct 29, 2009, at 9:26 PM, Michael Lerner wrote: I calculated temperatures with g_traj -f md2.trr -s md2.tpr -n index.ndx -ot -ng 5 I did not know you could get the temperature throught g_traj ... the -ot option is the one giving you the temperature? Yes. Do you get the same temperature difference with g_energy? Or may be the info is not in there? May be you could rerun you trajectory using a new tpr where you define two different temp coupling groups. And see if you observe the same temperature given by g_energy. It's possible that I'm doing something wrong, but I couldn't convince g_energy to give me the temperature of different groups when I used a single thermostat for the entire system. If I use two temperature coupling groups, I get the expected results. As Justin pointed out, it looks like the error was mine for not understanding the details of g_traj. When I manually correct for the degrees of freedom, I get temperatures for Protein/Sol that are close enough to each other (i.e. within a couple of standard deviations). Thanks, -Michael Thanks, -Michael -- begin md2.mdp -- ; RUN CONTROL PARAMETERS = integrator = md ; start time and timestep in ps = tinit = 0.0 dt = 0.020 nsteps = 2550 ; mode for center of mass motion removal comm-mode = linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm-grps = System ; OUTPUT CONTROL OPTIONS = ; Output frequency for coords (x), velocities (v) and forces (f) = nstxout = 5 nstvout = 5 nstfout = 0 ; Output frequency for energies to log file and energy file = nstlog = 2500 nstenergy = 2500 ; Output frequency and precision for xtc file = nstxtcout = 2500 xtc_precision = 100 ; This selects the subset of atoms for the xtc file. You can = ; select multiple groups. By default all atoms will be written. = xtc-grps = ; Selection of energy groups = energygrps = System Protein W WF ; NEIGHBORSEARCHING PARAMETERS = ; nblist update frequency = nstlist = 5 ; ns algorithm (simple or grid) = ns_type = grid ; Periodic boundary conditions: xyz or none = pbc = xyz ; nblist cut-off = rlist = 1.2 ; OPTIONS FOR ELECTROSTATICS AND VDW = ; Method for doing electrostatics = coulombtype = Shift rcoulomb_switch = 0.0 rcoulomb = 1.2 ; Dielectric constant (DC) for cut-off or DC of reaction field = epsilon_r = 15 ; Method for doing Van der Waals = vdw_type = Shift ; cut-off lengths = rvdw_switch = 0.9 rvdw = 1.2 ; Apply long range dispersion corrections for Energy and Pressure = DispCorr = No ; Spacing for the PME/PPPM FFT grid = fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used = fourier_nx = 10 fourier_ny = 10 fourier_nz = 10 ; EWALD/PME/PPPM parameters = pme_order = 4 ewald_rtol = 1e-05 epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS = ; Temperature coupling = tcoupl = Nose-Hoover ; Groups to couple separately = tc-grps = System ; Time constant (ps) and reference temperature (K) = tau_t = 2 ref_t = 323 ; Pressure coupling = Pcoupl = Parrinello-Rahman Pcoupltype = isotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) = tau_p = 5 compressibility = 5e-5 ref_p = 1.0 ; SIMULATED ANNEALING CONTROL = annealing = no ; GENERATE VELOCITIES FOR STARTUP RUN = continuation = yes gen_vel = no ;gen_temp = 323 ;gen_seed = 473529 ; OPTIONS FOR BONDS = constraints = none ; Type of constraint algorithm = constraint_algorithm = Lincs ; Relative tolerance of shake = shake_tol = 0.0001 ; Highest order in
