Re: [gmx-users] charge group moved too far between two domain decomposition step

2014-10-23 Thread Justin Lemkul



On 10/23/14 1:47 AM, Kester Wong wrote:

Thanks for the input Mark,



 Hi Mark,


 Thanks for the input, I thought a time step of 2 fs is small enough?


Not always for unconstrained bonds. The largest stable time step is
determined by the size of the period of the fastest oscillation, which is
vibrations of bonds to hydrogen. 2fs stretches the friendship there, and
you would certainly want to equilibrate thoroughly at a smaller time step
if you are doing messy things like freezing nanotubes and using pressure
coupling (which would not be my go-to method without a clear demonstration
that it produces valid results). More generally, equilibration of a
not-quite-right structure works better with a smaller timestep, because
that copes better with the artificially large forces you get...

Thanks for the very informative input, this should serve me well as a rule of
thumb for other systems as well.



 Should I change my constraint to constraints = h-bonds instead?


That would be normal.


 Or,


 Should I use:

 dt = 0.001; 1 fs

 constraints = none

 constraint-algorithm = shake


Plausible, but you want to run the final simulation with constraints for
the higher ns/day you obtain...

I think I got a little confused here.
For energy minimisation, I should impose no constraints to obtain the most
realistic starting structure for NVT.
As for NVT, one should use constraints to achieve a much quicker equilibration,
is that correct?



Constraints generally do not have any real effect on how fast the target 
properties reach an equilibrated state.  One thing to note - if you're going to 
run with constraints, you should minimize with constraints.  If you have 
problems during EM, setting the water to flexible can help, but you should 
minimize again with constraints, otherwise you often end up with distorted 
geometries that cannot be constrained during dynamics.



For instance (TIPS3P.itp as shown below), to achieve a higher ns/day in my NVT
run, I should have constraints = none in my mdp setting (using [bonds] and
[angles]);


If you want greater performance in terms of ns/day, you need a larger time step, 
in which case you do want constraints.  Now, the thing you have to consider is 
that if there is a [settles] block, constraints = none has no real effect 
since the SETTLE algorithm will be used (constraints present in the topology 
always override constraints = none in the .mdp file).



or have define = -DFLEXIBLE in my .mdp setting to achieve a better energy
minimisation (using SETTLES)?



See above.


[ moleculetype ]
; molname  nrexcl
SOL 2

[ atoms ]
; id   at type  res nr  residu name at name cg nr   charge
1   OT  1   SOL  OW 1   -0.834
2   HT  1   SOL HW1 10.417
3   HT  1   SOL HW2 10.417

#ifndef FLEXIBLE
[ settles ]
; i   j funct   length
1 1 0.09572 0.15139

[ exclusions ]
1 2 3
2 1 3
3 1 2
#else
[ bonds ]
; i j   funct   length  force.c.
1   2   1   0.09572 376560.0 0.09572376560.0
1   3   1   0.09572 376560.0 0.09572376560.0

[ angles ]
; i  j  k   funct   angle   force.c.
21  3   1   104.52  460.24  104.52  460.24
#endif


I would appreciate it if anyone could point me to the right direction, once and 
for all.

For constraints in a system with tips3p.itp, and hydronium.itp, how does one define 
FLEXIBLE or CONSTRAINTS?



define = -DFLEXIBLE in the .mdp file.  See the manual.


Below is the .itp for the hydronium model:

[ moleculetype ]
; molname   nrexcl
H3O  2

[ atoms ]
; idat type res nr  residu name at name  cg nr  charge   mass
1   H3O-O1  H3O  OW   1  -0.59   15.9994
2   H3O-H1  H3O H31   1   0.53   1.008
3   H3O-H1  H3O H32   1   0.53   1.008
4   H3O-H1  H3O H33   1   0.53   1.008

; use #ifdef FLEXIBLE or #ifdef CONSTRAINTS
#ifdef FLEXIBLE
[ bonds ]
;i  j   funct   length (nm)  force.c. (kJ/mol)
1   2   1   0.102  ; Hydronium OH bond
1   3   1   0.102  ; Wolf used 1000? check this
1   4   1   0.102  ; Hub used 400,000

[ angles ]
; i  j  k   funct   angle   force.c.
21  3   1   112 ; Hub used 400 force
41  3   1   112
21  4   1   112

#else
[ constraints ]
; ifunct   dohdhh
1 2  2   0.102
1 3  2   0.102
1 4  2   0.102
2 3  2   0.169124
3 4  2   0.169124
2 4  2   0.169124
#endif

[ exclusions ]
1  2  3  4
2  1  3  4
3  1  2  4
4  1  2  3


I tried the two approaches, by either having define = -DFLEXIBLE in 

Re: [gmx-users] charge group moved too far between two domain decomposition step

2014-10-23 Thread Kester Wong
Hi Justin and all, 
  Thanks for the input, I thought a time step of 2 fs is small enough?
 

 Not always for unconstrained bonds. The largest stable time step is
 determined by the size of the period of the fastest oscillation, which is
 vibrations of bonds to hydrogen. 2fs stretches the friendship there, and
 you would certainly want to equilibrate thoroughly at a smaller time step
 if you are doing messy things like freezing nanotubes and using pressure
 coupling (which would not be my go-to method without a clear demonstration
 that it produces valid results). More generally, equilibration of a
 not-quite-right structure works better with a smaller timestep, because
 that copes better with the artificially large forces you get...

 Thanks for the very informative input, this should serve me well as a rule of
 thumb for other systems as well.


  Should I change my constraint to constraints = h-bonds instead?
 

 That would be normal.


  Or,
 
 
  Should I use:
 
  dt = 0.001; 1 fs
 
  constraints = none
 
  constraint-algorithm = shake
 

 Plausible, but you want to run the final simulation with constraints for
 the higher ns/day you obtain...

 I think I got a little confused here.
 For energy minimisation, I should impose no constraints to obtain the most
 realistic starting structure for NVT.
 As for NVT, one should use constraints to achieve a much quicker equilibration,
 is that correct?


Constraints generally do not have any real effect on how fast the target 
properties reach an equilibrated state.  One thing to note - if you're going to 
run with constraints, you should minimize with constraints.  If you have 
problems during EM, setting the water to flexible can help, but you should 
minimize again with constraints, otherwise you often end up with distorted 
geometries that cannot be constrained during dynamics.Yes, I do want to use constraints in later calculations, so I will do the two-step minimisation as mentioned.Isn't setting water to flexible effectively calls for [ settles ], which means the calculation will be run with constraints, shouldn't the [ bonds ] and [ angles ] restraints be used instead? For this reason, I asked "For constraints in a system with tips3p.ipt, and hydronium.itp, how does one define "FLEXIBLE" or "CONSTRAINTS"" in the previous email.Always a good thing to be cautious about the distorted geo. that cannot be constrained during MD can be the root of all evil.
 For instance (TIPS3P.itp as shown below), to achieve a higher ns/day in my NVT
 run, I should have constraints = none in my mdp setting (using [bonds] and
 [angles]);

If you want greater performance in terms of ns/day, you need a larger time step, 
in which case you do want constraints.  Now, the thing you have to consider is 
that if there is a [settles] block, "constraints = none" has no real effect 
since the SETTLE algorithm will be used (constraints present in the topology 
always override "constraints = none" in the .mdp file).
This is exactly what I wanted to find out, it is clear to me now that constraints = none will have no effect since GROMACS will use what ever constraints defined in the topology (if any), as stated in the manual. Thank you for that.According to the manual, SHAKE cannot be used with energy minimisation. How does one minimise the water (settles) and ions system without using SHAKE, given that the maxwarn-angle for LINCS is 90 degree and the H-O-H angle of hydronium (H3O) is 112 degree?Can both LINCS and SHAKE be turned off during minimisation without constraints (constraints = none, i.e. using [ bonds ] and [ angles ] for water and hydronium), and then constraint-algorithm = SHAKE for the second minimisation with constraints (i.e. using [ settles ] and [ constraints ] for water and hydronium, respectively)? or have define = -DFLEXIBLE in my .mdp setting to achieve a better energy
 minimisation (using SETTLES)?


See above.

 [ moleculetype ]
 ; molname  nrexcl
 SOL 2

 [ atoms ]
 ; id   at type  res nr  residu name at name cg nr   charge
 1   OT  1   SOL  OW 1   -0.834
 2   HT  1   SOL HW1 10.417
 3   HT  1   SOL HW2 10.417

 #ifndef FLEXIBLE
 [ settles ]
 ; i   j funct   length
 1 1 0.09572 0.15139

 [ exclusions ]
 1 2 3
 2 1 3
 3 1 2
 #else
 [ bonds ]
 ; i j   funct   length  force.c.
 1   2   1   0.09572 376560.0 0.09572376560.0
 1   3   1   0.09572 376560.0 0.09572376560.0

 [ angles ]
 ; i  j  k   funct   angle   force.c.
 21  3   1   104.52  460.24  104.52  460.24
 #endif


 I would appreciate it if anyone 

Re: [gmx-users] charge group moved too far between two domain decomposition step

2014-10-23 Thread Justin Lemkul



On 10/23/14 10:19 PM, Kester Wong wrote:

Hi Justin and all,

 
  Thanks for the input, I thought a time step of 2 fs is small enough?
 

 Not always for unconstrained bonds. The largest stable time step is
 determined by the size of the period of the fastest oscillation, 
which is
 vibrations of bonds to hydrogen. 2fs stretches the friendship there, 
and
 you would certainly want to equilibrate thoroughly at a smaller time 
step
 if you are doing messy things like freezing nanotubes and using 
pressure
 coupling (which would not be my go-to method without a clear 
demonstration
 that it produces valid results). More generally, equilibration of a
 not-quite-right structure works better with a smaller timestep, 
because
 that copes better with the artificially large forces you get...

 Thanks for the very informative input, this should serve me well as a 
rule of
 thumb for other systems as well.


  Should I change my constraint to constraints = h-bonds instead?
 

 That would be normal.


  Or,
 
 
  Should I use:
 
  dt = 0.001; 1 fs
 
  constraints = none
 
  constraint-algorithm = shake
 

 Plausible, but you want to run the final simulation with constraints 
for
 the higher ns/day you obtain...

 I think I got a little confused here.
 For energy minimisation, I should impose no constraints to obtain the most
 realistic starting structure for NVT.
 As for NVT, one should use constraints to achieve a much quicker 
equilibration,
 is that correct?


Constraints generally do not have any real effect on how fast the target
properties reach an equilibrated state.  One thing to note - if you're 
going to
run with constraints, you should minimize with constraints.  If you have
problems during EM, setting the water to flexible can help, but you should
minimize again with constraints, otherwise you often end up with distorted
geometries that cannot be constrained during dynamics.


Yes, I do want to use constraints in later calculations, so I will do the
two-step minimisation as mentioned.

Isn't setting water to flexible effectively calls for [ settles ], which means
the calculation will be run with constraints,
shouldn't the [ bonds ] and [ angles ] restraints be used instead? For this
reason, I asked For constraints in a system with tips3p.ipt, and hydronium.itp,
how does one define FLEXIBLE or CONSTRAINTS in the previous email.



Note that the topology uses different constructs.  One is #ifndef, the other is 
#ifdef.  With #ifndef, it says if FLEXIBLE is not defined, use SETTLE, in 
contrast to the more direct use of #ifdef.



Always a good thing to be cautious about the distorted geo. that cannot be
constrained during MD can be the root of all evil.

 For instance (TIPS3P.itp as shown below), to achieve a higher ns/day in 
my NVT
 run, I should have constraints = none in my mdp setting (using [bonds] and
 [angles]);

If you want greater performance in terms of ns/day, you need a larger time 
step,
in which case you do want constraints.  Now, the thing you have to consider 
is
that if there is a [settles] block, constraints = none has no real effect
since the SETTLE algorithm will be used (constraints present in the topology
always override constraints = none in the .mdp file).


This is exactly what I wanted to find out, it is clear to me now that
constraints = none will have no effect since GROMACS will use what ever
constraints defined in the topology (if any), as stated in the manual. Thank you
for that.



LINCS is more stable than SHAKE, anyway, so it is preferable.


According to the manual, SHAKE cannot be used with energy minimisation.
How does one minimise the water (settles) and ions system without using SHAKE,
given that the maxwarn-angle for LINCS is 90 degree and the H-O-H angle of
hydronium (H3O) is 112 degree?



Bond rotations and topology angles are unrelated.


Can both LINCS and SHAKE be turned off during minimisation without constraints
(constraints = none, i.e. using [ bonds ] and [ angles ] for water and 
hydronium),
and then constraint-algorithm = SHAKE for the second minimisation with
constraints (i.e. using [ settles ] and [ constraints ] for water and hydronium,
respectively)?



Not if the statement in the manual is accurate, but you can constrain with 
LINCS.



 or have define = -DFLEXIBLE in my .mdp setting to achieve a better energy
 minimisation (using SETTLES)?


See above.

 [ moleculetype ]
 ; molname  nrexcl
 SOL 2

 [ atoms ]
 ; id   at type  res nr  residu name at name cg nr   charge
 1   OT  1   SOL  OW 1   

Re: [gmx-users] charge group moved too far between two domain decomposition step

2014-10-22 Thread Mark Abraham
On Wed, Oct 22, 2014 at 6:37 AM, Kester Wong kester2...@ibs.re.kr wrote:

 Hi Mark,


 Thanks for the input, I thought a time step of 2 fs is small enough?


Not always for unconstrained bonds. The largest stable time step is
determined by the size of the period of the fastest oscillation, which is
vibrations of bonds to hydrogen. 2fs stretches the friendship there, and
you would certainly want to equilibrate thoroughly at a smaller time step
if you are doing messy things like freezing nanotubes and using pressure
coupling (which would not be my go-to method without a clear demonstration
that it produces valid results). More generally, equilibration of a
not-quite-right structure works better with a smaller timestep, because
that copes better with the artificially large forces you get...


 Should I change my constraint to constraints = h-bonds instead?


That would be normal.


 Or,


 Should I use:

 dt = 0.001; 1 fs

 constraints = none

 constraint-algorithm = shake


Plausible, but you want to run the final simulation with constraints for
the higher ns/day you obtain...

Mark

Cheers,

 Kester


 You should start by using a time step more appropriate to your use of
 constraints = none, or vice versa.

 Mark

 On Wed, Oct 22, 2014 at 4:26 AM, Kester Wong  wrote:

  Dear all,
 
 
  I am intrigued by the error that I have received this morning, as follows:
 
 
 
  The charge group starting at atom 37920 moved more than the distance
  allowed by the domain decomposition (1.40) in direction Y
 
  distance out of cell 7.144197
 
  Old coordinates:   54.499   39.697   14.058
 
  New coordinates:   54.499   39.697   14.058
 
  Old cell boundaries in direction Y:   19.468   32.846
 
  New cell boundaries in direction Y:   18.667   32.553
 
 
  If the charge group moved, if any, then shouldn't the new coordinates be
  showing different values?
 
  It was a new NPT calculation that I continued from a 20ns NVT run, so it
  should be equilibrated well enough, for a simple water on graphene system.
 
 
  Using GROMACS 5.0, I was unable to start the NPT with Berendsen
  thermostat, so I had to go with Parrinello-Rahman with a tau-P at 5.0 (my
  tau-T was 0.5, following the general rule of thumb, that tau-P should be
  larger than tau-T to allow thermal degrees of freedom to equilibrate faster
  than the box size).
 
 
  Should I reduce tau-P?
 
 
 
 
  Below is my NPT.mdp parameters:
 
  ; RUN CONTROL PARAMETERS
 
  integrator   = md
 
  tinit= 2.000
 
  dt   = 0.002;  2 fs
 
  nsteps   = 500  ; 10,000 ps
 
  comm-mode= Linear
 
  nstcomm  = 10
 
 
  emtol= 0.001
 
  niter= 30
 
 
  ; OUTPUT CONTROL OPTIONS
 
  nstxout  = 5000; save coordinates every 8 ps
 
  nstvout  = 5000; save velocities every 8 ps
 
  nstfout  = 0
 
  nstlog   = 5000; update log file every 8 ps
 
  nstenergy= 5000; save energies every 8 ps
 
  ;nstxtcout is replaced with nstxout-compressed
 
  nstxout-compressed   = 5000; no. of steps between writing
  coordinates using lossy compression
 
  ;xtc-precision is replaced with compressed-x-precision
 
  compressed-x-precision   = 3000; precision with which to write to the
  compressed trajectory file
 
 
  ; Periodic boundary conditions: xyz (default), no (vacuum)
 
  cutoff-scheme= Verlet ; default for GROMACS-5.0
 
  nstlist  = 1  ; A smaller nstlist may reduce LINCS
  warnings
 
  ns-type  = grid   ;
 
  pbc  = xyz
 
  ;periodic_molecules   = yes  ; cannot be used with SHAKE
 
  ; nblist cut-off
 
  rlist= 1.40   ; short-range neighbourlist cutoff
 
  nstcalclr= 1
 
 
  coulombtype  = PME
 
  coulomb-modifier = Potential-shift-Verlet
 
  r_coulomb= 1.40  ; short-range electrostatic cutoff
 
  rcoulomb-switch  = 0 ; used when potential-switch
  (coulomb-modifier) is specified
 
  epsilon-r= 1 ; default value=1
 
  pme_order= 4
 
  fourierspacing   = 0.12
 
  fourier_nx   = 0
 
  fourier_ny   = 0
 
  fourier_nz   = 0
 
  ewald_rtol   = 1e-05
 
  ewald_geometry   = 3d
 
  epsilon_surface  = 0  ; Use 0 when simulation is done with PME
 
 
 
  ;vdw-type = switch is replaced with the following two flags in GROMACS-5.0
 
  vdwtype  = Cut-off
 
  vdw_modifier = Potential-switch
 
  ; cut-off lengths
 
  rvdw-switch  = 1.00
  rvdw = 1.40   ; short-range van der Waals cutoff
 
  constraints  = none
  constraint-algorithm = SHAKE
  continuation = yes; starting from NVT run
  Shake-SOR  

Re: [gmx-users] charge group moved too far between two domain decomposition step

2014-10-22 Thread Kester Wong
Thanks for the input Mark,
 Hi Mark,


 Thanks for the input, I thought a time step of 2 fs is small enough?


Not always for unconstrained bonds. The largest stable time step is
determined by the size of the period of the fastest oscillation, which is
vibrations of bonds to hydrogen. 2fs stretches the friendship there, and
you would certainly want to equilibrate thoroughly at a smaller time step
if you are doing messy things like freezing nanotubes and using pressure
coupling (which would not be my go-to method without a clear demonstration
that it produces valid results). More generally, equilibration of a
not-quite-right structure works better with a smaller timestep, because
that copes better with the artificially large forces you get...
Thanks for the very informative input, this should serve me well as a rule of thumb for other systems as well.
 Should I change my constraint to constraints = h-bonds instead?


That would be normal.


 Or,


 Should I use:

 dt = 0.001; 1 fs

 constraints = none

 constraint-algorithm = shake


Plausible, but you want to run the final simulation with constraints for
the higher ns/day you obtain...
I think I got a little confused here.For energy minimisation, I should impose no constraints to obtain the most realistic starting structure for NVT.As for NVT, one should use constraints to achieve a much quicker equilibration, is that correct?For instance (TIPS3P.itp as shown below), to achieve a higher ns/day in my NVT run, I should have constraints = none in my mdp setting (using [bonds] and [angles]);or have define = -DFLEXIBLE in my .mdp setting to achieve a better energy minimisation (using SETTLES)?[ moleculetype ]; molname   nrexclSOL   2[ atoms ]; id  at type res nr residu name   at name cg nr  charge1OT   1SOL   OW   1-0.8342HT   1SOL   HW1   10.4173HT   1SOL   HW2   10.417#ifndef FLEXIBLE[ settles ]; ij   funct  length1 1   0.09572 0.15139[ exclusions ]1 2 32 1 33 1 2#else[ bonds ]; i   jfunct  length force.c.1210.09572 376560.0 0.09572376560.01310.09572 376560.0 0.09572376560.0[ angles ]; i   j   kfunct  angle  force.c.21   31104.52 460.24 104.52 460.24#endifI would appreciate it if anyone could point me to the right direction, once and for all. For constraints in a system with tips3p.itp, and hydronium.itp, how does one define "FLEXIBLE" or "CONSTRAINTS"?Below is the .itp for the hydronium model:[ moleculetype ]
; molname   nrexcl
H3O  2

[ atoms ]
; idat type res nr  residu name at name  cg nr  charge   mass
1   H3O-O1  H3O  OW   1  -0.59   15.9994
2   H3O-H1  H3O H31   1   0.53   1.008
3   H3O-H1  H3O H32   1   0.53   1.008
4   H3O-H1  H3O H33   1   0.53   1.008

; use #ifdef FLEXIBLE or #ifdef CONSTRAINTS
#ifdef FLEXIBLE
[ bonds ]
;i  j   funct   length (nm)  force.c. (kJ/mol)
1   2   1   0.102  ; Hydronium OH bond
1   3   1   0.102  ; Wolf used 1000? check this
1   4   1   0.102  ; Hub used 400,000

[ angles ]
; i  j  k   funct   angle   force.c.
21  3   1   112 ; Hub used 400 force
41  3   1   112
21  4   1   112

#else
[ constraints ]
; ifunct   dohdhh
1 2  2   0.102
1 3  2   0.102
1 4  2   0.102
2 3  2   0.169124
3 4  2   0.169124
2 4  2   0.169124
#endif

[ exclusions ]
1  2  3  4
2  1  3  4
3  1  2  4
4  1  2  3
I tried the two approaches, by either having define = -DFLEXIBLE in my .mdp file, or #define FLEXIBLE in my topol.top file,and received this error: Incorrect number of parameters - found 1, expected 2 or 4 for Bond.Setting the LD random seed to 3234307754
Generated 207 of the 210 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 0.5
Generated 210 of the 210 1-4 parameter combinations
Thank you for your time and patience.Regards,Kester
Mark

Cheers,

 Kester


 You should start by using a time step more appropriate to your use of
 constraints = none, or vice versa.

 Mark

 On Wed, Oct 22, 2014 at 4:26 AM, Kester Wong  wrote:

  Dear all,
 
 
  I am intrigued by the error that I have received this morning, as follows:
 
 
 
  The charge group starting at atom 37920 moved more than the distance
  allowed by the domain decomposition (1.40) in direction Y
 
  distance out of cell 7.144197
 
  Old coordinates:   54.499   39.697   14.058
 
  New coordinates:   54.499   39.697   14.058
 
  Old cell boundaries in direction Y:   19.468   32.846
 
  New cell boundaries in 

Re: [gmx-users] charge group moved too far between two domain decomposition step

2014-10-21 Thread Kester Wong
Hi Mark,Thanks for the input, I thought a time step of 2 fs is small enough?Should I change my constraint to constraints = h-bonds instead?Or,Should I use:dt = 0.001; 1 fsconstraints = noneconstraint-algorithm = shakeCheers,Kester
You should start by using a time step more appropriate to your use of
constraints = none, or vice versa.

Mark

On Wed, Oct 22, 2014 at 4:26 AM, Kester Wong  wrote:

 Dear all,


 I am intrigued by the error that I have received this morning, as follows:



 The charge group starting at atom 37920 moved more than the distance
 allowed by the domain decomposition (1.40) in direction Y

 distance out of cell 7.144197

 Old coordinates:   54.499   39.697   14.058

 New coordinates:   54.499   39.697   14.058

 Old cell boundaries in direction Y:   19.468   32.846

 New cell boundaries in direction Y:   18.667   32.553


 If the charge group moved, if any, then shouldn't the new coordinates be
 showing different values?

 It was a new NPT calculation that I continued from a 20ns NVT run, so it
 should be equilibrated well enough, for a simple water on graphene system.


 Using GROMACS 5.0, I was unable to start the NPT with Berendsen
 thermostat, so I had to go with Parrinello-Rahman with a tau-P at 5.0 (my
 tau-T was 0.5, following the general rule of thumb, that tau-P should be
 larger than tau-T to allow thermal degrees of freedom to equilibrate faster
 than the box size).


 Should I reduce tau-P?




 Below is my NPT.mdp parameters:

 ; RUN CONTROL PARAMETERS

 integrator   = md

 tinit= 2.000

 dt   = 0.002;  2 fs

 nsteps   = 500  ; 10,000 ps

 comm-mode= Linear

 nstcomm  = 10


 emtol= 0.001

 niter= 30


 ; OUTPUT CONTROL OPTIONS

 nstxout  = 5000; save coordinates every 8 ps

 nstvout  = 5000; save velocities every 8 ps

 nstfout  = 0

 nstlog   = 5000; update log file every 8 ps

 nstenergy= 5000; save energies every 8 ps

 ;nstxtcout is replaced with nstxout-compressed

 nstxout-compressed   = 5000; no. of steps between writing
 coordinates using lossy compression

 ;xtc-precision is replaced with compressed-x-precision

 compressed-x-precision   = 3000; precision with which to write to the
 compressed trajectory file


 ; Periodic boundary conditions: xyz (default), no (vacuum)

 cutoff-scheme= Verlet ; default for GROMACS-5.0

 nstlist  = 1  ; A smaller nstlist may reduce LINCS
 warnings

 ns-type  = grid   ;

 pbc  = xyz

 ;periodic_molecules   = yes  ; cannot be used with SHAKE

 ; nblist cut-off

 rlist= 1.40   ; short-range neighbourlist cutoff

 nstcalclr= 1


 coulombtype  = PME

 coulomb-modifier = Potential-shift-Verlet

 r_coulomb= 1.40  ; short-range electrostatic cutoff

 rcoulomb-switch  = 0 ; used when potential-switch
 (coulomb-modifier) is specified

 epsilon-r= 1 ; default value=1

 pme_order= 4

 fourierspacing   = 0.12

 fourier_nx   = 0

 fourier_ny   = 0

 fourier_nz   = 0

 ewald_rtol   = 1e-05

 ewald_geometry   = 3d

 epsilon_surface  = 0  ; Use 0 when simulation is done with PME



 ;vdw-type = switch is replaced with the following two flags in GROMACS-5.0

 vdwtype  = Cut-off

 vdw_modifier = Potential-switch

 ; cut-off lengths

 rvdw-switch  = 1.00
 rvdw = 1.40   ; short-range van der Waals cutoff

 constraints  = none
 constraint-algorithm = SHAKE
 continuation = yes; starting from NVT run
 Shake-SOR= no
 shake-tol= 0.0001
 ;lincs-order  = 4 ; accuracy of LINCS
 ;lincs-iter   = 24; also related to accuracy
 ; Lincs will write a warning to the stderr if in one step a bond
 ; rotates over more degrees than
 ;lincs-warnangle  = 30

 Tcoupl   = v-rescale
 nh-chain-length  = 1; the leapfrog md integrator only supports
 1
 ; Groups to couple separately
 tc-grps  = GRA SOL
 ; Time constant (ps) and reference temperature (K)
 tau-t= 0.5 0.5
 ref-t= 300 300
 nsttcouple   = 20  ; added by k-wong ; not specified in
 Yanbin's MDP

 Pcoupl   = Parrinello-Rahman
 Pcoupltype   = isotropic
 ; Time constant (ps), compressibility (1/bar) and reference P