- Posts: 4
Temperature of Martini Cholesterol Is Low
- da294
-
Topic Author
- Offline
- Fresh Boarder
Less
More
9 years 11 months ago #3831
by da294
Temperature of Martini Cholesterol Is Low was created by da294
Hello,
I have been working with multi-component Martini bilayers, including DPPC, DUPC and cholesterol. I couple the water and each lipid type separately to the temperature bath using the v-rescale thermostat, tau_t= 1 and ref_t = 295, as suggested elsewhere. When I use g_energy to report the temperatures of the various molecules, I get that DPPC, DUPC and W all are 295. However, the temperature of cholesterol is ~289-290. This happens in a variety of systems I have looked at. Why is that happening? Is it due to the comparatively fewer motional degrees of freedom of cholesterol compared to the lipids? My mdp file is copied below.
-David
_________________________________________________________________________________________
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 02 feb 2013 by DdJ
;
; for use with GROMACS 4.5/4.6
;
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable
; with dt=40 fs, some (especially rings and polarizable water) require 20-30 fs.
; Note that time steps of 40 fs and larger may create local heating or
; cooling in your system. Although the use of a heat bath will globally
; remove this effect, it is advised to check consistency of
; your results for somewhat smaller time steps in the range 20-30 fs.
; Time steps exceeding 40 fs should not be used; time steps smaller
; than 20 fs are also not required unless specifically stated in the itp file.
integrator = md
dt = 0.02
nsteps = 1250000000
nstcomm = 10
comm-grps = W DPPC DUPC CHOL
nstxout = 250000
nstvout = 250000
nstfout = 0
nstlog = 250000
nstenergy = 250000
nstxtcout = 0
xtc_precision = 0
xtc-grps =
energygrps = DPPC DUPC CHOL W
; NEIGHBOURLIST and MARTINI
; Due to the use of shifted potentials, the noise generated
; from particles leaving/entering the neighbour list is not so large,
; even when large time steps are being used. In practice, once every
; ten steps works fine with a neighborlist cutoff that is equal to the
; non-bonded cutoff (1.2 nm). However, to improve energy conservation
; or to avoid local heating/cooling, you may increase the update frequency
; and/or enlarge the neighbourlist cut-off (to 1.4 nm). The latter option
; is computationally less expensive and leads to improved energy conservation
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.2
; MARTINI and NONBONDED
; Standard cut-off schemes are used for the non-bonded interactions
; in the Martini model: LJ interactions are shifted to zero in the
; range 0.9-1.2 nm, and electrostatic interactions in the range 0.0-1.2 nm.
; The treatment of the non-bonded cut-offs is considered to be part of
; the force field parameterization, so we recommend not to touch these
; values as they will alter the overall balance of the force field.
; In principle you can include long range electrostatics through the use
; of PME, which could be more realistic in certain applications
; Please realize that electrostatic interactions in the Martini model are
; not considered to be very accurate to begin with, especially as the
; screening in the system is set to be uniform across the system with
; a screening constant of 15. When using PME, please make sure your
; system properties are still reasonable.
;
; With the polarizable water model, the relative electrostatic screening
; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
; apolar solvent. The polarizable water itself will perform the explicit screening
; in aqueous environment. In this case, the use of PME is more realistic.
;
; For use in combination with the Verlet-pairlist algorithm implemented
; in Gromacs 4.6 a straight cutoff in combination with the potential
; modifiers can be used. Although this will change the potential shape,
; preliminary results indicate that forcefield properties do not change a lot
; when the LJ cutoff is reduced to 1.1 nm. Be sure to test the effects for
; your particular system. The advantage is a gain of speed of 50-100%.
coulombtype = Shift ;Reaction_field (for use with Verlet-pairlist) ;PME (especially with polarizable water)
rcoulomb_switch = 0.0
rcoulomb = 1.2
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = Shift ;cutoff (for use with Verlet-pairlist)
rvdw_switch = 0.9
rvdw = 1.2 ;1.1 (for use with Verlet-pairlist)
;cutoff-scheme = verlet
;coulomb-modifier = Potential-shift
;vdw-modifier = Potential-shift
;epsilon_rf = 0 ; epsilon_rf = 0 really means epsilon_rf = infinity
;verlet-buffer-drift = 0.005
; MARTINI and TEMPERATURE/PRESSURE
; normal temperature and pressure coupling schemes can be used.
; It is recommended to couple individual groups in your system separately.
; Good temperature control can be achieved with the velocity rescale (V-rescale)
; thermostat using a coupling constant of the order of 1 ps. Even better
; temperature control can be achieved by reducing the temperature coupling
; constant to 0.1 ps, although with such tight coupling (approaching
; the time step) one can no longer speak of a weak-coupling scheme.
; We therefore recommend a coupling time constant of at least 0.5 ps.
; The Berendsen thermostat is less suited since it does not give
; a well described thermodynamic ensemble.
;
; Pressure can be controlled with the Parrinello-Rahman barostat,
; with a coupling constant in the range 4-8 ps and typical compressibility
; in the order of 10-4 - 10-5 bar-1. Note that, for equilibration purposes,
; the Berendsen thermostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the pressure
; coupling should be done semiisotropic.
tcoupl = v-rescale
tc-grps = DPPC DUPC CHOL W
tau_t = 1.0 1.0 1.0 1.0
ref_t = 295 295 295 295
Pcoupl = Berendsen
Pcoupltype = semiisotropic
tau_p = 4.0 4.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 5e-5 5e-5
ref_p = 1.0 1.0
gen_vel = yes
gen_temp = 295
gen_seed = -1
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
unconstrained_start = no
lincs_order = 4
lincs_warnangle = 30
I have been working with multi-component Martini bilayers, including DPPC, DUPC and cholesterol. I couple the water and each lipid type separately to the temperature bath using the v-rescale thermostat, tau_t= 1 and ref_t = 295, as suggested elsewhere. When I use g_energy to report the temperatures of the various molecules, I get that DPPC, DUPC and W all are 295. However, the temperature of cholesterol is ~289-290. This happens in a variety of systems I have looked at. Why is that happening? Is it due to the comparatively fewer motional degrees of freedom of cholesterol compared to the lipids? My mdp file is copied below.
-David
_________________________________________________________________________________________
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 02 feb 2013 by DdJ
;
; for use with GROMACS 4.5/4.6
;
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable
; with dt=40 fs, some (especially rings and polarizable water) require 20-30 fs.
; Note that time steps of 40 fs and larger may create local heating or
; cooling in your system. Although the use of a heat bath will globally
; remove this effect, it is advised to check consistency of
; your results for somewhat smaller time steps in the range 20-30 fs.
; Time steps exceeding 40 fs should not be used; time steps smaller
; than 20 fs are also not required unless specifically stated in the itp file.
integrator = md
dt = 0.02
nsteps = 1250000000
nstcomm = 10
comm-grps = W DPPC DUPC CHOL
nstxout = 250000
nstvout = 250000
nstfout = 0
nstlog = 250000
nstenergy = 250000
nstxtcout = 0
xtc_precision = 0
xtc-grps =
energygrps = DPPC DUPC CHOL W
; NEIGHBOURLIST and MARTINI
; Due to the use of shifted potentials, the noise generated
; from particles leaving/entering the neighbour list is not so large,
; even when large time steps are being used. In practice, once every
; ten steps works fine with a neighborlist cutoff that is equal to the
; non-bonded cutoff (1.2 nm). However, to improve energy conservation
; or to avoid local heating/cooling, you may increase the update frequency
; and/or enlarge the neighbourlist cut-off (to 1.4 nm). The latter option
; is computationally less expensive and leads to improved energy conservation
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.2
; MARTINI and NONBONDED
; Standard cut-off schemes are used for the non-bonded interactions
; in the Martini model: LJ interactions are shifted to zero in the
; range 0.9-1.2 nm, and electrostatic interactions in the range 0.0-1.2 nm.
; The treatment of the non-bonded cut-offs is considered to be part of
; the force field parameterization, so we recommend not to touch these
; values as they will alter the overall balance of the force field.
; In principle you can include long range electrostatics through the use
; of PME, which could be more realistic in certain applications
; Please realize that electrostatic interactions in the Martini model are
; not considered to be very accurate to begin with, especially as the
; screening in the system is set to be uniform across the system with
; a screening constant of 15. When using PME, please make sure your
; system properties are still reasonable.
;
; With the polarizable water model, the relative electrostatic screening
; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
; apolar solvent. The polarizable water itself will perform the explicit screening
; in aqueous environment. In this case, the use of PME is more realistic.
;
; For use in combination with the Verlet-pairlist algorithm implemented
; in Gromacs 4.6 a straight cutoff in combination with the potential
; modifiers can be used. Although this will change the potential shape,
; preliminary results indicate that forcefield properties do not change a lot
; when the LJ cutoff is reduced to 1.1 nm. Be sure to test the effects for
; your particular system. The advantage is a gain of speed of 50-100%.
coulombtype = Shift ;Reaction_field (for use with Verlet-pairlist) ;PME (especially with polarizable water)
rcoulomb_switch = 0.0
rcoulomb = 1.2
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = Shift ;cutoff (for use with Verlet-pairlist)
rvdw_switch = 0.9
rvdw = 1.2 ;1.1 (for use with Verlet-pairlist)
;cutoff-scheme = verlet
;coulomb-modifier = Potential-shift
;vdw-modifier = Potential-shift
;epsilon_rf = 0 ; epsilon_rf = 0 really means epsilon_rf = infinity
;verlet-buffer-drift = 0.005
; MARTINI and TEMPERATURE/PRESSURE
; normal temperature and pressure coupling schemes can be used.
; It is recommended to couple individual groups in your system separately.
; Good temperature control can be achieved with the velocity rescale (V-rescale)
; thermostat using a coupling constant of the order of 1 ps. Even better
; temperature control can be achieved by reducing the temperature coupling
; constant to 0.1 ps, although with such tight coupling (approaching
; the time step) one can no longer speak of a weak-coupling scheme.
; We therefore recommend a coupling time constant of at least 0.5 ps.
; The Berendsen thermostat is less suited since it does not give
; a well described thermodynamic ensemble.
;
; Pressure can be controlled with the Parrinello-Rahman barostat,
; with a coupling constant in the range 4-8 ps and typical compressibility
; in the order of 10-4 - 10-5 bar-1. Note that, for equilibration purposes,
; the Berendsen thermostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the pressure
; coupling should be done semiisotropic.
tcoupl = v-rescale
tc-grps = DPPC DUPC CHOL W
tau_t = 1.0 1.0 1.0 1.0
ref_t = 295 295 295 295
Pcoupl = Berendsen
Pcoupltype = semiisotropic
tau_p = 4.0 4.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 5e-5 5e-5
ref_p = 1.0 1.0
gen_vel = yes
gen_temp = 295
gen_seed = -1
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
unconstrained_start = no
lincs_order = 4
lincs_warnangle = 30
Please Log in or Create an account to join the conversation.
- initram
-
- Offline
- Fresh Boarder
Less
More
- Posts: 13
9 years 8 months ago #3990
by initram
Replied by initram on topic Temperature of Martini Cholesterol Is Low
Hi,
I see something similar during equilibration of cholesteryl esters.
In my case g_energy reports 300.688K although the ref_t=304K.
Would it help if I use Berendsen thermostat instead of v-rescale?
My mdp parameters are:
integrator = md
dt = 0.02
nsteps = 5000000
nstcomm = 100
comm-grps =
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 1000 ; Output frequency for energies to log file
nstenergy = 100 ; Output frequency for energies to energy file
nstxtcout = 1000 ; Output frequency for .xtc file
xtc_precision = 100
xtc-grps =
energygrps = CO
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.4
coulombtype = Shift ;Reaction_field (for use with Verlet-pairlist) ;PME (especially with polarizable water)
rcoulomb_switch = 0.0
rcoulomb = 1.2
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = Shift ;cutoff (for use with Verlet-pairlist)
rvdw_switch = 0.9
rvdw = 1.2 ;1.1 (for use with Verlet-pairlist)
tcoupl = v-rescale
tc-grps = CO
tau_t = 0.5
ref_t = 304
pcoupl = no
gen_vel = yes
gen_temp = 304
gen_seed = 173529
constraints = none
constraint_algorithm = Lincs
continuation = no
lincs_order = 4
lincs_warnangle = 30
I see something similar during equilibration of cholesteryl esters.
In my case g_energy reports 300.688K although the ref_t=304K.
Would it help if I use Berendsen thermostat instead of v-rescale?
My mdp parameters are:
integrator = md
dt = 0.02
nsteps = 5000000
nstcomm = 100
comm-grps =
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 1000 ; Output frequency for energies to log file
nstenergy = 100 ; Output frequency for energies to energy file
nstxtcout = 1000 ; Output frequency for .xtc file
xtc_precision = 100
xtc-grps =
energygrps = CO
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.4
coulombtype = Shift ;Reaction_field (for use with Verlet-pairlist) ;PME (especially with polarizable water)
rcoulomb_switch = 0.0
rcoulomb = 1.2
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = Shift ;cutoff (for use with Verlet-pairlist)
rvdw_switch = 0.9
rvdw = 1.2 ;1.1 (for use with Verlet-pairlist)
tcoupl = v-rescale
tc-grps = CO
tau_t = 0.5
ref_t = 304
pcoupl = no
gen_vel = yes
gen_temp = 304
gen_seed = 173529
constraints = none
constraint_algorithm = Lincs
continuation = no
lincs_order = 4
lincs_warnangle = 30
Please Log in or Create an account to join the conversation.
- xavier
-
- Offline
- Admin
Less
More
- Posts: 416
9 years 8 months ago #4017
by xavier
Replied by xavier on topic Temperature of Martini Cholesterol Is Low
Hi David,
We have observed cholesterol to have a slightly lower temperature than the rest of the membrane. As you may know some people in the lab have worked on an improved model of cholesterol that addresses a few issues of the old model. It is about to be finalised so hopefully more people can use it.
Also It is not clear to me is the temperature given in the output reflects entirely the use of constrains on this molecules. There are many and potentially redundant so the determination of the number of degrees of freedom taken aways might be tricky.
It is also not clear if one should use a separate temperature bath for cholesterol! You might artificially pump heat into your system that will be taken away somewhere else but another bath.
It is also not clear what cholesterol temperature in a lipid mixture means. In real life would the cholesterol have a lower temperature due to lower degrees of mobility? These are actually very interesting questions.
We have observed cholesterol to have a slightly lower temperature than the rest of the membrane. As you may know some people in the lab have worked on an improved model of cholesterol that addresses a few issues of the old model. It is about to be finalised so hopefully more people can use it.
Also It is not clear to me is the temperature given in the output reflects entirely the use of constrains on this molecules. There are many and potentially redundant so the determination of the number of degrees of freedom taken aways might be tricky.
It is also not clear if one should use a separate temperature bath for cholesterol! You might artificially pump heat into your system that will be taken away somewhere else but another bath.
It is also not clear what cholesterol temperature in a lipid mixture means. In real life would the cholesterol have a lower temperature due to lower degrees of mobility? These are actually very interesting questions.
da294 wrote: Hello,
I have been working with multi-component Martini bilayers, including DPPC, DUPC and cholesterol. I couple the water and each lipid type separately to the temperature bath using the v-rescale thermostat, tau_t= 1 and ref_t = 295, as suggested elsewhere. When I use g_energy to report the temperatures of the various molecules, I get that DPPC, DUPC and W all are 295. However, the temperature of cholesterol is ~289-290. This happens in a variety of systems I have looked at. Why is that happening? Is it due to the comparatively fewer motional degrees of freedom of cholesterol compared to the lipids? My mdp file is copied below.
-David
_________________________________________________________________________________________
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 02 feb 2013 by DdJ
;
; for use with GROMACS 4.5/4.6
;
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable
; with dt=40 fs, some (especially rings and polarizable water) require 20-30 fs.
; Note that time steps of 40 fs and larger may create local heating or
; cooling in your system. Although the use of a heat bath will globally
; remove this effect, it is advised to check consistency of
; your results for somewhat smaller time steps in the range 20-30 fs.
; Time steps exceeding 40 fs should not be used; time steps smaller
; than 20 fs are also not required unless specifically stated in the itp file.
integrator = md
dt = 0.02
nsteps = 1250000000
nstcomm = 10
comm-grps = W DPPC DUPC CHOL
nstxout = 250000
nstvout = 250000
nstfout = 0
nstlog = 250000
nstenergy = 250000
nstxtcout = 0
xtc_precision = 0
xtc-grps =
energygrps = DPPC DUPC CHOL W
; NEIGHBOURLIST and MARTINI
; Due to the use of shifted potentials, the noise generated
; from particles leaving/entering the neighbour list is not so large,
; even when large time steps are being used. In practice, once every
; ten steps works fine with a neighborlist cutoff that is equal to the
; non-bonded cutoff (1.2 nm). However, to improve energy conservation
; or to avoid local heating/cooling, you may increase the update frequency
; and/or enlarge the neighbourlist cut-off (to 1.4 nm). The latter option
; is computationally less expensive and leads to improved energy conservation
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.2
; MARTINI and NONBONDED
; Standard cut-off schemes are used for the non-bonded interactions
; in the Martini model: LJ interactions are shifted to zero in the
; range 0.9-1.2 nm, and electrostatic interactions in the range 0.0-1.2 nm.
; The treatment of the non-bonded cut-offs is considered to be part of
; the force field parameterization, so we recommend not to touch these
; values as they will alter the overall balance of the force field.
; In principle you can include long range electrostatics through the use
; of PME, which could be more realistic in certain applications
; Please realize that electrostatic interactions in the Martini model are
; not considered to be very accurate to begin with, especially as the
; screening in the system is set to be uniform across the system with
; a screening constant of 15. When using PME, please make sure your
; system properties are still reasonable.
;
; With the polarizable water model, the relative electrostatic screening
; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
; apolar solvent. The polarizable water itself will perform the explicit screening
; in aqueous environment. In this case, the use of PME is more realistic.
;
; For use in combination with the Verlet-pairlist algorithm implemented
; in Gromacs 4.6 a straight cutoff in combination with the potential
; modifiers can be used. Although this will change the potential shape,
; preliminary results indicate that forcefield properties do not change a lot
; when the LJ cutoff is reduced to 1.1 nm. Be sure to test the effects for
; your particular system. The advantage is a gain of speed of 50-100%.
coulombtype = Shift ;Reaction_field (for use with Verlet-pairlist) ;PME (especially with polarizable water)
rcoulomb_switch = 0.0
rcoulomb = 1.2
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = Shift ;cutoff (for use with Verlet-pairlist)
rvdw_switch = 0.9
rvdw = 1.2 ;1.1 (for use with Verlet-pairlist)
;cutoff-scheme = verlet
;coulomb-modifier = Potential-shift
;vdw-modifier = Potential-shift
;epsilon_rf = 0 ; epsilon_rf = 0 really means epsilon_rf = infinity
;verlet-buffer-drift = 0.005
; MARTINI and TEMPERATURE/PRESSURE
; normal temperature and pressure coupling schemes can be used.
; It is recommended to couple individual groups in your system separately.
; Good temperature control can be achieved with the velocity rescale (V-rescale)
; thermostat using a coupling constant of the order of 1 ps. Even better
; temperature control can be achieved by reducing the temperature coupling
; constant to 0.1 ps, although with such tight coupling (approaching
; the time step) one can no longer speak of a weak-coupling scheme.
; We therefore recommend a coupling time constant of at least 0.5 ps.
; The Berendsen thermostat is less suited since it does not give
; a well described thermodynamic ensemble.
;
; Pressure can be controlled with the Parrinello-Rahman barostat,
; with a coupling constant in the range 4-8 ps and typical compressibility
; in the order of 10-4 - 10-5 bar-1. Note that, for equilibration purposes,
; the Berendsen thermostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the pressure
; coupling should be done semiisotropic.
tcoupl = v-rescale
tc-grps = DPPC DUPC CHOL W
tau_t = 1.0 1.0 1.0 1.0
ref_t = 295 295 295 295
Pcoupl = Berendsen
Pcoupltype = semiisotropic
tau_p = 4.0 4.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 5e-5 5e-5
ref_p = 1.0 1.0
gen_vel = yes
gen_temp = 295
gen_seed = -1
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
unconstrained_start = no
lincs_order = 4
lincs_warnangle = 30
Please Log in or Create an account to join the conversation.
Time to create page: 0.103 seconds