- Posts: 3
Lipid bilayer + Peptide simulation
- prithvi
-
Topic Author
- Offline
- Fresh Boarder
Less
More
7 years 7 months ago #5866
by prithvi
Lipid bilayer + Peptide simulation was created by prithvi
Dear MartiniFF experts,
I am trying to perform Martini CG simulations containing peptide helices embedded into POPC bilayer with Gromacs 5.1.2. I have prepared the system using the insane.py script provided on the Martini website. I am running the system using the mdp file (martini_v2.x_new.mdp) obtained from Martini website. The entries of the mdp file are as follows -
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 15 Jul 2015 by DdJ
;
; for use with GROMACS 5
; For a thorough comparison of different mdp options in combination with the Martini force field, see:
; D.H. de Jong et al., Martini straight: boosting performance using a shorter cutoff and GPUs, submitted.
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable with dt=40 fs,
; however better energy conservation is achieved using a
; 20-30 fs timestep.
; Time steps smaller than 20 fs are not required unless specifically stated in the itp file.
integrator = md
dt = 0.02
nsteps = 200000000
nstcomm = 100
comm-grps = System
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 10000
nstenergy = 1000
nstxout-compressed = 10000
compressed-x-precision = 100
compressed-x-grps =
energygrps = Protein POPC W
; NEIGHBOURLIST and MARTINI
; To achieve faster simulations in combination with the Verlet-neighborlist
; scheme, Martini can be simulated with a straight cutoff. In order to
; do so, the cutoff distance is reduced 1.1 nm.
; Neighborlist length should be optimized depending on your hardware setup:
; updating ever 20 steps should be fine for classic systems, while updating
; every 30-40 steps might be better for GPU based systems.
; The Verlet neighborlist scheme will automatically choose a proper neighborlist
; length, based on a energy drift tolerance.
;
; Coulomb interactions can alternatively be treated using a reaction-field,
; giving slightly better properties.
; Please realize that electrostVatic 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.
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
coulombtype = cutoff
coulomb-modifier = Potential-shift-verlet
rcoulomb = 1.1
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
; 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 10e-4 - 10e-5 bar-1. Note that, for equilibration purposes,
; the Berendsen barostat 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 = Protein POPC W
tau_t = 1.0 1.0 1.0
ref_t = 310 310 310
Pcoupl = parrinello-rahman
Pcoupltype = semiisotropic
tau_p = 12.0 ;12.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 3e-4 3e-4
ref_p = 1.0 1.0
gen_vel = no
gen_temp = 320
gen_seed = 473529
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
As you may observe that I am running the simulation for 4 microseconds. The system runs perfectly fine till 4 microsecond. But at the end of 4 microseconds I see that the POPC bilayer containing the peptide has moved till to the end (along Z-axis) of the simulation box. I also tried changing comm-grps = System to comm-grps = Protein POPC W. The bilayer still moved to the end of the box. Is it normal or is something not right?
Thanking you,
Prithvi Raj Pandey
I am trying to perform Martini CG simulations containing peptide helices embedded into POPC bilayer with Gromacs 5.1.2. I have prepared the system using the insane.py script provided on the Martini website. I am running the system using the mdp file (martini_v2.x_new.mdp) obtained from Martini website. The entries of the mdp file are as follows -
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 15 Jul 2015 by DdJ
;
; for use with GROMACS 5
; For a thorough comparison of different mdp options in combination with the Martini force field, see:
; D.H. de Jong et al., Martini straight: boosting performance using a shorter cutoff and GPUs, submitted.
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable with dt=40 fs,
; however better energy conservation is achieved using a
; 20-30 fs timestep.
; Time steps smaller than 20 fs are not required unless specifically stated in the itp file.
integrator = md
dt = 0.02
nsteps = 200000000
nstcomm = 100
comm-grps = System
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 10000
nstenergy = 1000
nstxout-compressed = 10000
compressed-x-precision = 100
compressed-x-grps =
energygrps = Protein POPC W
; NEIGHBOURLIST and MARTINI
; To achieve faster simulations in combination with the Verlet-neighborlist
; scheme, Martini can be simulated with a straight cutoff. In order to
; do so, the cutoff distance is reduced 1.1 nm.
; Neighborlist length should be optimized depending on your hardware setup:
; updating ever 20 steps should be fine for classic systems, while updating
; every 30-40 steps might be better for GPU based systems.
; The Verlet neighborlist scheme will automatically choose a proper neighborlist
; length, based on a energy drift tolerance.
;
; Coulomb interactions can alternatively be treated using a reaction-field,
; giving slightly better properties.
; Please realize that electrostVatic 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.
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
coulombtype = cutoff
coulomb-modifier = Potential-shift-verlet
rcoulomb = 1.1
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
; 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 10e-4 - 10e-5 bar-1. Note that, for equilibration purposes,
; the Berendsen barostat 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 = Protein POPC W
tau_t = 1.0 1.0 1.0
ref_t = 310 310 310
Pcoupl = parrinello-rahman
Pcoupltype = semiisotropic
tau_p = 12.0 ;12.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 3e-4 3e-4
ref_p = 1.0 1.0
gen_vel = no
gen_temp = 320
gen_seed = 473529
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
As you may observe that I am running the simulation for 4 microseconds. The system runs perfectly fine till 4 microsecond. But at the end of 4 microseconds I see that the POPC bilayer containing the peptide has moved till to the end (along Z-axis) of the simulation box. I also tried changing comm-grps = System to comm-grps = Protein POPC W. The bilayer still moved to the end of the box. Is it normal or is something not right?
Thanking you,
Prithvi Raj Pandey
Please Log in or Create an account to join the conversation.
- peterkroon
-
- Offline
- Gold Boarder
Less
More
- Posts: 210
7 years 7 months ago #5870
by peterkroon
Replied by peterkroon on topic Lipid bilayer + Peptide simulation
It sounds to me it's simple diffusion across the periodic boundary conditions, and as such completely normal and expected.
Please Log in or Create an account to join the conversation.
- mnmelo
-
- Offline
- Admin
Less
More
- Posts: 89
7 years 7 months ago #5871
by mnmelo
Replied by mnmelo on topic Lipid bilayer + Peptide simulation
Hi Prithvi,
I'm not sure I agree with Peter. If your membrane starts at the center of the box in z, 4 µs does not seem enough to get so much drift (assuming your box has a sensible height, of at least 8 nm, and also that system construction with insane.py yielded something reasonable and roughly homogeneously solvated).
In any case, COM motion removal should definitely dampen such drift and it is intriguing that it didn't.
I can't see what might cause this without having a look at the actual data. Perhaps you can post the starting structure or a plot of the bilayer COM z as a function of time for both systems?
A couple of notes on COM motion removal:
Let us know,
Manel
I'm not sure I agree with Peter. If your membrane starts at the center of the box in z, 4 µs does not seem enough to get so much drift (assuming your box has a sensible height, of at least 8 nm, and also that system construction with insane.py yielded something reasonable and roughly homogeneously solvated).
In any case, COM motion removal should definitely dampen such drift and it is intriguing that it didn't.
I can't see what might cause this without having a look at the actual data. Perhaps you can post the starting structure or a plot of the bilayer COM z as a function of time for both systems?
A couple of notes on COM motion removal:
- COM motion removal only periodically resets the COM velocities. It doesn't put your structure back, which means small drift can still accumulate;
- Groups that are momentum-coupled should have their COM motion removed together. That's the case for your bilayer and peptides. This means that your comm-grps should be something like W POPC_or_Protein (you'll have to construct this last group with gmx make_ndx). The same goes for the thermostatting.
Let us know,
Manel
Please Log in or Create an account to join the conversation.
- prithvi
-
Topic Author
- Offline
- Fresh Boarder
Less
More
- Posts: 3
7 years 7 months ago - 7 years 7 months ago #5875
by prithvi
Replied by prithvi on topic Lipid bilayer + Peptide simulation
Dear Manel,
Thank you for your reply. The box dimension along Z is 10 nm. I performed
another simulation with the peptide and lipid in one index group and
coupling the COM and thermostat for this group. The modified mdp file is as
follows (the non-wat group refers to the combined index)
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 15 Jul 2015 by DdJ
;
; for use with GROMACS 5
; For a thorough comparison of different mdp options in combination with
the Martini force field, see:
; D.H. de Jong et al., Martini straight: boosting performance using a
shorter cutoff and GPUs, submitted.
title =3D Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable with dt=3D40 fs,
; however better energy conservation is achieved using a
; 20-30 fs timestep.
; Time steps smaller than 20 fs are not required unless specifically stated
in the itp file.
integrator =3D md
dt =3D 0.02
nsteps =3D 200000000 ;500000
nstcomm =3D 100
comm-grps =3D non-wat W
nstxout =3D 0
nstvout =3D 0
nstfout =3D 0
nstlog =3D 10000
nstenergy =3D 1000
nstxout-compressed =3D 10000
compressed-x-precision =3D 100
compressed-x-grps =3D
energygrps =3D Protein POPC W
; NEIGHBOURLIST and MARTINI
; To achieve faster simulations in combination with the Verlet-neighborlist
; scheme, Martini can be simulated with a straight cutoff. In order to
; do so, the cutoff distance is reduced 1.1 nm.
; Neighborlist length should be optimized depending on your hardware setup:
; updating ever 20 steps should be fine for classic systems, while updating
; every 30-40 steps might be better for GPU based systems.
; The Verlet neighborlist scheme will automatically choose a proper
neighborlist
; length, based on a energy drift tolerance.
;
; Coulomb interactions can alternatively be treated using a reaction-field,
; giving slightly better properties.
; Please realize that electrostVatic 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-dielectri=
c
; apolar solvent. The polarizable water itself will perform the explicit
screening
; in aqueous environment. In this case, the use of PME is more realistic.
cutoff-scheme =3D Verlet
nstlist =3D 20
ns_type =3D grid
pbc =3D xyz
verlet-buffer-tolerance =3D 0.005
coulombtype =3D cutoff
coulomb-modifier =3D Potential-shift-verlet
rcoulomb =3D 1.1
epsilon_r =3D 15 ; 2.5 (with polarizable water)
vdw_type =3D cutoff
vdw-modifier =3D Potential-shift-verlet
rvdw =3D 1.1
; 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 10e-4 - 10e-5 bar-1. Note that, for equilibration
purposes,
; the Berendsen barostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the
pressure
; coupling should be done semiisotropic.
tcoupl =3D v-rescale
tc-grps =3D non-wat W
tau_t =3D 1.0 1.0 ;1.0
ref_t =3D 310 310 ;310
Pcoupl =3D parrinello-rahman
Pcoupltype =3D semiisotropic
tau_p =3D 12.0 ;12.0 ;parrinello-rahman is more stable
with larger tau-p, DdJ, 20130422
compressibility =3D 3e-4 3e-4
ref_p =3D 1.0 1.0
gen_vel =3D no
gen_temp =3D 320
gen_seed =3D 473529
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
In this case too, the bilayer drifted along Z-axis. The plot for COM for
the bilayer as a function of time is in the link below for two simulations -
comm-grps = non-wat W and comm-grps = System.
drive.google.com/file/d/0B7O8OCB9pyxpMlJ...eFk/view?usp=sharing
Thank you for your reply. The box dimension along Z is 10 nm. I performed
another simulation with the peptide and lipid in one index group and
coupling the COM and thermostat for this group. The modified mdp file is as
follows (the non-wat group refers to the combined index)
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 15 Jul 2015 by DdJ
;
; for use with GROMACS 5
; For a thorough comparison of different mdp options in combination with
the Martini force field, see:
; D.H. de Jong et al., Martini straight: boosting performance using a
shorter cutoff and GPUs, submitted.
title =3D Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable with dt=3D40 fs,
; however better energy conservation is achieved using a
; 20-30 fs timestep.
; Time steps smaller than 20 fs are not required unless specifically stated
in the itp file.
integrator =3D md
dt =3D 0.02
nsteps =3D 200000000 ;500000
nstcomm =3D 100
comm-grps =3D non-wat W
nstxout =3D 0
nstvout =3D 0
nstfout =3D 0
nstlog =3D 10000
nstenergy =3D 1000
nstxout-compressed =3D 10000
compressed-x-precision =3D 100
compressed-x-grps =3D
energygrps =3D Protein POPC W
; NEIGHBOURLIST and MARTINI
; To achieve faster simulations in combination with the Verlet-neighborlist
; scheme, Martini can be simulated with a straight cutoff. In order to
; do so, the cutoff distance is reduced 1.1 nm.
; Neighborlist length should be optimized depending on your hardware setup:
; updating ever 20 steps should be fine for classic systems, while updating
; every 30-40 steps might be better for GPU based systems.
; The Verlet neighborlist scheme will automatically choose a proper
neighborlist
; length, based on a energy drift tolerance.
;
; Coulomb interactions can alternatively be treated using a reaction-field,
; giving slightly better properties.
; Please realize that electrostVatic 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-dielectri=
c
; apolar solvent. The polarizable water itself will perform the explicit
screening
; in aqueous environment. In this case, the use of PME is more realistic.
cutoff-scheme =3D Verlet
nstlist =3D 20
ns_type =3D grid
pbc =3D xyz
verlet-buffer-tolerance =3D 0.005
coulombtype =3D cutoff
coulomb-modifier =3D Potential-shift-verlet
rcoulomb =3D 1.1
epsilon_r =3D 15 ; 2.5 (with polarizable water)
vdw_type =3D cutoff
vdw-modifier =3D Potential-shift-verlet
rvdw =3D 1.1
; 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 10e-4 - 10e-5 bar-1. Note that, for equilibration
purposes,
; the Berendsen barostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the
pressure
; coupling should be done semiisotropic.
tcoupl =3D v-rescale
tc-grps =3D non-wat W
tau_t =3D 1.0 1.0 ;1.0
ref_t =3D 310 310 ;310
Pcoupl =3D parrinello-rahman
Pcoupltype =3D semiisotropic
tau_p =3D 12.0 ;12.0 ;parrinello-rahman is more stable
with larger tau-p, DdJ, 20130422
compressibility =3D 3e-4 3e-4
ref_p =3D 1.0 1.0
gen_vel =3D no
gen_temp =3D 320
gen_seed =3D 473529
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
In this case too, the bilayer drifted along Z-axis. The plot for COM for
the bilayer as a function of time is in the link below for two simulations -
comm-grps = non-wat W and comm-grps = System.
drive.google.com/file/d/0B7O8OCB9pyxpMlJ...eFk/view?usp=sharing
Last edit: 7 years 7 months ago by prithvi.
Please Log in or Create an account to join the conversation.
- prithvi
-
Topic Author
- Offline
- Fresh Boarder
Less
More
- Posts: 3
7 years 7 months ago #5876
by prithvi
Replied by prithvi on topic Lipid bilayer + Peptide simulation
In the meantime I tried running a pure POPC bilayer water system. I prepared the system using CHARMM-GUI. The system contains 128 lipid molecules in each leaflet (total 256). And while generating the system with CHARMM-GUI 32 Angstrom layer of water was chosen. Then I performed all the equilibration and production run steps suggested for systems generated from CHARMM-GUI. Following that I performed 4 microsecond simulation for this POPC bilayer system using the martini_v2.x_new.mdp file provided on Martini website. Entries for which are as follows -
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 15 Jul 2015 by DdJ
;
; for use with GROMACS 5
; For a thorough comparison of different mdp options in combination with the Martini force field, see:
; D.H. de Jong et al., Martini straight: boosting performance using a shorter cutoff and GPUs, submitted.
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable with dt=40 fs,
; however better energy conservation is achieved using a
; 20-30 fs timestep.
; Time steps smaller than 20 fs are not required unless specifically stated in the itp file.
integrator = md
dt = 0.02
nsteps = 400000000
nstcomm = 100
comm-grps = POPC W
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 20000
nstenergy = 20000
nstxout-compressed = 20000
compressed-x-precision = 100
compressed-x-grps =
energygrps = POPC W
; NEIGHBOURLIST and MARTINI
; To achieve faster simulations in combination with the Verlet-neighborlist
; scheme, Martini can be simulated with a straight cutoff. In order to
; do so, the cutoff distance is reduced 1.1 nm.
; Neighborlist length should be optimized depending on your hardware setup:
; updating ever 20 steps should be fine for classic systems, while updating
; every 30-40 steps might be better for GPU based systems.
; The Verlet neighborlist scheme will automatically choose a proper neighborlist
; length, based on a energy drift tolerance.
;
; Coulomb interactions can alternatively be treated using a reaction-field,
; giving slightly better properties.
; Please realize that electrostVatic 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.
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
coulombtype = cutoff
coulomb-modifier = Potential-shift-verlet
rcoulomb = 1.1
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
; 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 10e-4 - 10e-5 bar-1. Note that, for equilibration purposes,
; the Berendsen barostat 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 = POPC W
tau_t = 1.0 1.0
ref_t = 310 310
Pcoupl = parrinello-rahman
Pcoupltype = semiisotropic
tau_p = 12.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 3e-4 3e-4
ref_p = 1.0 1.0
gen_vel = no
gen_temp = 320
gen_seed = 473529
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
I performed two simulations, one with dt = 0.02, and the other with dt = 0.01. For the simulation with dt = 0.02, the bilayer again moved to the end of the box. But for dt = 0.01 the bilayer did not drift and stayed near the box centre. Is it necessary to use dt=0.01 for such system?
Thanking you,
Prithvi
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 15 Jul 2015 by DdJ
;
; for use with GROMACS 5
; For a thorough comparison of different mdp options in combination with the Martini force field, see:
; D.H. de Jong et al., Martini straight: boosting performance using a shorter cutoff and GPUs, submitted.
title = Martini
; TIMESTEP IN MARTINI
; Most simulations are numerically stable with dt=40 fs,
; however better energy conservation is achieved using a
; 20-30 fs timestep.
; Time steps smaller than 20 fs are not required unless specifically stated in the itp file.
integrator = md
dt = 0.02
nsteps = 400000000
nstcomm = 100
comm-grps = POPC W
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 20000
nstenergy = 20000
nstxout-compressed = 20000
compressed-x-precision = 100
compressed-x-grps =
energygrps = POPC W
; NEIGHBOURLIST and MARTINI
; To achieve faster simulations in combination with the Verlet-neighborlist
; scheme, Martini can be simulated with a straight cutoff. In order to
; do so, the cutoff distance is reduced 1.1 nm.
; Neighborlist length should be optimized depending on your hardware setup:
; updating ever 20 steps should be fine for classic systems, while updating
; every 30-40 steps might be better for GPU based systems.
; The Verlet neighborlist scheme will automatically choose a proper neighborlist
; length, based on a energy drift tolerance.
;
; Coulomb interactions can alternatively be treated using a reaction-field,
; giving slightly better properties.
; Please realize that electrostVatic 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.
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
coulombtype = cutoff
coulomb-modifier = Potential-shift-verlet
rcoulomb = 1.1
epsilon_r = 15 ; 2.5 (with polarizable water)
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
; 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 10e-4 - 10e-5 bar-1. Note that, for equilibration purposes,
; the Berendsen barostat 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 = POPC W
tau_t = 1.0 1.0
ref_t = 310 310
Pcoupl = parrinello-rahman
Pcoupltype = semiisotropic
tau_p = 12.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 3e-4 3e-4
ref_p = 1.0 1.0
gen_vel = no
gen_temp = 320
gen_seed = 473529
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs
I performed two simulations, one with dt = 0.02, and the other with dt = 0.01. For the simulation with dt = 0.02, the bilayer again moved to the end of the box. But for dt = 0.01 the bilayer did not drift and stayed near the box centre. Is it necessary to use dt=0.01 for such system?
Thanking you,
Prithvi
Please Log in or Create an account to join the conversation.
- mnmelo
-
- Offline
- Admin
Less
More
- Posts: 89
7 years 7 months ago #5904
by mnmelo
Replied by mnmelo on topic Lipid bilayer + Peptide simulation
Sorry for the delay replying.
You shouldn't need a 0.01 ps timestep.
My feeling is that there might be some GROMACS inconsistency with the new Verlet/GPU algorithm. What happens if you use the old-style mdp parameters? (with cutoff-scheme = group, rlist = 1.4 and nstlist = 10)
You shouldn't need a 0.01 ps timestep.
My feeling is that there might be some GROMACS inconsistency with the new Verlet/GPU algorithm. What happens if you use the old-style mdp parameters? (with cutoff-scheme = group, rlist = 1.4 and nstlist = 10)
Please Log in or Create an account to join the conversation.
Time to create page: 0.110 seconds