<p style="border:0; padding:0; margin:0; font-family:'Gulim'; font-size:10pt; cursor: text;"><p style="font-family: Gulim; font-size: 10pt; border: 0px; padding: 0px; margin: 0px; cursor: text;">Hi Justin and all,&nbsp;</p><blockquote style="font-family: Gulim; font-size: 12px; border-left-style: solid; border-left-width: 2px; margin: 0pt 0pt 0pt 0.8em; padding-left: 1em;"><div style="font-family:arial,돋움;line-height:1.5"><div style="margin-top:5px;"><pre>&gt;     &gt;
&gt;     &gt; Thanks for the input, I thought a time step of 2 fs is small enough?
&gt;     &gt;
&gt;
&gt;     Not always for unconstrained bonds. The largest stable time step is
&gt;     determined by the size of the period of the fastest oscillation, which is
&gt;     vibrations of bonds to hydrogen. 2fs stretches the friendship there, and
&gt;     you would certainly want to equilibrate thoroughly at a smaller time step
&gt;     if you are doing messy things like freezing nanotubes and using pressure
&gt;     coupling (which would not be my go-to method without a clear demonstration
&gt;     that it produces valid results). More generally, equilibration of a
&gt;     not-quite-right structure works better with a smaller timestep, because
&gt;     that copes better with the artificially large forces you get...
&gt;
&gt; Thanks for the very informative input, this should serve me well as a rule of
&gt; thumb for other systems as well.
&gt;
&gt;
&gt;     &gt; Should I change my constraint to constraints = h-bonds instead?
&gt;     &gt;
&gt;
&gt;     That would be normal.
&gt;
&gt;
&gt;     &gt; Or,
&gt;     &gt;
&gt;     &gt;
&gt;     &gt; Should I use:
&gt;     &gt;
&gt;     &gt; dt = 0.001        ; 1 fs
&gt;     &gt;
&gt;     &gt; constraints = none
&gt;     &gt;
&gt;     &gt; constraint-algorithm = shake
&gt;     &gt;
&gt;
&gt;     Plausible, but you want to run the final simulation with constraints for
&gt;     the higher ns/day you obtain...
&gt;
&gt; I think I got a little confused here.
&gt; For energy minimisation, I should impose no constraints to obtain the most
&gt; realistic starting structure for NVT.
&gt; As for NVT, one should use constraints to achieve a much quicker equilibration,
&gt; is that correct?
&gt;

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.</pre><pre><br></pre></div></div></blockquote><font face="monospace">Yes, I do want to use constraints in later calculations, so I will do the two-step minimisation as mentioned.<br></font><font face="monospace"><br></font><span style="font-family: monospace; font-size: 12px; line-height: 18px; white-space: pre;">Isn't setting water to flexible effectively calls for [ settles ], which means the calculation will be run with constraints, <br></span><font face="monospace">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.<br></font><font face="monospace"><br></font><span style="font-family: monospace;">Always a good thing to be cautious about the distorted geo. that cannot be constrained during MD can be the root of all evil.</span><font face="monospace"><span style="font-size: 12px; line-height: 18px; white-space: pre;"><br></span></font><blockquote style="font-family: Gulim; font-size: 12px; border-left-style: solid; border-left-width: 2px; margin: 0pt 0pt 0pt 0.8em; padding-left: 1em;"><div style="font-family:arial,돋움;line-height:1.5"><div style="margin-top:5px;"><pre>
&gt; For instance (TIPS3P.itp as shown below), to achieve a higher ns/day in my NVT
&gt; run, I should have constraints = none in my mdp setting (using [bonds] and
&gt; [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).
</pre><pre><br></pre></div></div></blockquote><font face="monospace">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.<br><br></font><font face="monospace"><span style="font-size: 12px; line-height: 18px; white-space: pre;">According to the manual, SHAKE cannot be used with energy minimisation. </span></font><br><font face="monospace"><span style="font-size: 12px; line-height: 18px; white-space: pre;">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?<br></span></font><blockquote style="font-size: 12px; border-left-style: solid; border-left-width: 2px; margin: 0pt 0pt 0pt 0.8em; padding-left: 1em;"><div style="font-family: arial, 돋움; line-height: 1.5;"><div style="margin-top: 5px;"></div></div></blockquote><font face="monospace"><span style="font-size: 12px; line-height: 18px; white-space: pre;">Can both LINCS and SHAKE be turned off during minimisation without constraints (constraints = none, i.e. using [ bonds ] and [ angles ] for water and hydronium), </span></font><br><font face="monospace"><span style="font-size: 12px; line-height: 18px; white-space: pre;">and then constraint-algorithm = SHAKE for the second minimisation with constraints (i.e. using [ settles ] and [ constraints ] for water and hydronium, respectively)?<br></span></font><blockquote style="font-size: 12px; border-left-style: solid; border-left-width: 2px; margin: 0pt 0pt 0pt 0.8em; padding-left: 1em;"><div style="font-family: arial, 돋움; line-height: 1.5;"><div style="margin-top: 5px;"></div></div></blockquote><blockquote style="font-family: Gulim; font-size: 12px; border-left-style: solid; border-left-width: 2px; margin: 0pt 0pt 0pt 0.8em; padding-left: 1em;"><div style="font-family:arial,돋움;line-height:1.5"><div style="margin-top:5px;"><pre><br></pre><pre>&gt; or have define = -DFLEXIBLE in my .mdp setting to achieve a better energy
&gt; minimisation (using SETTLES)?
&gt;

See above.

&gt; [ moleculetype ]
&gt; ; molname      nrexcl
&gt; SOL             2
&gt;
&gt; [ atoms ]
&gt; ; id   at type  res nr  residu name     at name         cg nr   charge
&gt; 1       OT      1       SOL              OW             1       -0.834
&gt; 2       HT      1       SOL             HW1             1        0.417
&gt; 3       HT      1       SOL             HW2             1        0.417
&gt;
&gt; #ifndef FLEXIBLE
&gt; [ settles ]
&gt; ; i       j     funct   length
&gt; 1         1     0.09572 0.15139
&gt;
&gt; [ exclusions ]
&gt; 1 2 3
&gt; 2 1 3
&gt; 3 1 2
&gt; #else
&gt; [ bonds ]
&gt; ; i     j       funct   length  force.c.
&gt; 1       2       1       0.09572 376560.0 0.09572        376560.0
&gt; 1       3       1       0.09572 376560.0 0.09572        376560.0
&gt;
&gt; [ angles ]
&gt; ; i      j      k       funct   angle   force.c.
&gt; 2        1      3       1       104.52  460.24  104.52  460.24
&gt; #endif
&gt;
&gt;
&gt; I would appreciate it if anyone could point me to the right direction, once and for all.
&gt;
&gt; For constraints in a system with tips3p.itp, and hydronium.itp, how does one define "FLEXIBLE" or "CONSTRAINTS"?
&gt;

define = -DFLEXIBLE in the .mdp file.  See the manual.
</pre></div></div></blockquote><font face="monospace">I ought to write the question differently, I was meant to ask, how does one tell GROMACS not to impose constraints on both water and hydronium.&nbsp;<br><span style="font-size: 12px; line-height: 18px; white-space: pre;"><br></span></font><blockquote style="font-family: Gulim; font-size: 12px; border-left-style: solid; border-left-width: 2px; margin: 0pt 0pt 0pt 0.8em; padding-left: 1em;"><div style="font-family:arial,돋움;line-height:1.5"><div style="margin-top:5px;"><pre><span style="font-family: arial, 돋움; line-height: 1.5;">&gt; Below is the .itp for the hydronium model:</span></pre><pre>&gt;
&gt; [ moleculetype ]
&gt; ; molname   nrexcl
&gt; H3O          2
&gt;
&gt; [ atoms ]
&gt; ; id    at type     res nr  residu name at name  cg nr  charge   mass
&gt; 1       H3O-O        1      H3O          OW       1      -0.59   15.9994
&gt; 2       H3O-H        1      H3O         H31       1       0.53   1.008
&gt; 3       H3O-H        1      H3O         H32       1       0.53   1.008
&gt; 4       H3O-H        1      H3O         H33       1       0.53   1.008
&gt;
&gt; ; use #ifdef FLEXIBLE or #ifdef CONSTRAINTS
&gt; #ifdef FLEXIBLE
&gt; [ bonds ]
&gt; ;i  j   funct   length (nm)  force.c. (kJ/mol)
&gt; 1   2   1       0.102      ; Hydronium OH bond
&gt; 1   3   1       0.102      ; Wolf used 1000? check this
&gt; 1   4   1       0.102      ; Hub used 400,000
&gt;
&gt; [ angles ]
&gt; ; i      j      k       funct   angle   force.c.
&gt; 2        1      3       1       112     ; Hub used 400 force
&gt; 4        1      3       1       112
&gt; 2        1      4       1       112
&gt;
&gt; #else
&gt; [ constraints ]
&gt; ; i    funct   doh    dhh
&gt; 1 2      2       0.102
&gt; 1 3      2       0.102
&gt; 1 4      2       0.102
&gt; 2 3      2       0.169124
&gt; 3 4      2       0.169124
&gt; 2 4      2       0.169124
&gt; #endif
&gt;
&gt; [ exclusions ]
&gt; 1  2  3  4
&gt; 2  1  3  4
&gt; 3  1  2  4
&gt; 4  1  2  3
&gt;
&gt;
&gt; I tried the two approaches, by either having define = -DFLEXIBLE in my .mdp file, or #define FLEXIBLE in my topol.top file,
&gt;

Either using the "define" keyword or explicitly using #define in the topology 
will trigger the same behavior.  The more flexible approach is to use .mdp 
settings, lest you perhaps forget one day and end up doing the run incorrectly 
because you forgot you hacked the topology!
</pre></div></div></blockquote><font face="monospace">Thanks for the tip! Yes, I do make changes more regularly on .mdp, less so on the topology; one can forget to amend the topology for that matter.&nbsp;<br><br>Thank you for your input Justin!<br><br><span style="font-size: 12px; line-height: 18px; white-space: pre;">Regards,</span><br><span style="font-size: 12px; line-height: 18px; white-space: pre;">Kester<br></span></font><blockquote style="font-family: Gulim; font-size: 12px; border-left-style: solid; border-left-width: 2px; margin: 0pt 0pt 0pt 0.8em; padding-left: 1em;"><div style="font-family:arial,돋움;line-height:1.5"><div style="margin-top:5px;"><pre><br></pre><pre>-Justin

&gt; and received this error: Incorrect number of parameters - found 1, expected 2 or 4 for Bond.
&gt;
&gt; Setting the LD random seed to 3234307754
&gt; Generated 207 of the 210 non-bonded parameter combinations
&gt; Generating 1-4 interactions: fudge = 0.5
&gt; Generated 210 of the 210 1-4 parameter combinations
&gt;
&gt;
&gt; Thank you for your time and patience.
&gt;
&gt; Regards,
&gt;
&gt; Kester
&gt;
&gt;
&gt;     Mark
&gt;
&gt;     Cheers,
&gt;     &gt;
&gt;     &gt; Kester
&gt;     &gt;
&gt;     &gt;
&gt;     &gt; You should start by using a time step more appropriate to your use of
&gt;     &gt; constraints = none, or vice versa.
&gt;     &gt;
&gt;     &gt; Mark
&gt;     &gt;
&gt;     &gt; On Wed, Oct 22, 2014 at 4:26 AM, Kester Wong  wrote:
&gt;     &gt;
&gt;     &gt; &gt; Dear all,
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; I am intrigued by the error that I have received this morning, as follows:
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; The charge group starting at atom 37920 moved more than the distance
&gt;     &gt; &gt; allowed by the domain decomposition (1.400000) in direction Y
&gt;     &gt; &gt;
&gt;     &gt; &gt; distance out of cell 7.144197
&gt;     &gt; &gt;
&gt;     &gt; &gt; Old coordinates:   54.499   39.697   14.058
&gt;     &gt; &gt;
&gt;     &gt; &gt; New coordinates:   54.499   39.697   14.058
&gt;     &gt; &gt;
&gt;     &gt; &gt; Old cell boundaries in direction Y:   19.468   32.846
&gt;     &gt; &gt;
&gt;     &gt; &gt; New cell boundaries in direction Y:   18.667   32.553
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; If the charge group moved, if any, then shouldn't the new coordinates be
&gt;     &gt; &gt; showing different values?
&gt;     &gt; &gt;
&gt;     &gt; &gt; It was a new NPT calculation that I continued from a 20ns NVT run, so it
&gt;     &gt; &gt; should be equilibrated well enough, for a simple water on graphene system.
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; Using GROMACS 5.0, I was unable to start the NPT with Berendsen
&gt;     &gt; &gt; thermostat, so I had to go with Parrinello-Rahman with a tau-P at 5.0 (my
&gt;     &gt; &gt; tau-T was 0.5, following the general rule of thumb, that tau-P should be
&gt;     &gt; &gt; larger than tau-T to allow thermal degrees of freedom to equilibrate faster
&gt;     &gt; &gt; than the box size).
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; Should I reduce tau-P?
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; Below is my NPT.mdp parameters:
&gt;     &gt; &gt;
&gt;     &gt; &gt; ; RUN CONTROL PARAMETERS
&gt;     &gt; &gt;
&gt;     &gt; &gt; integrator               = md
&gt;     &gt; &gt;
&gt;     &gt; &gt; tinit                    = 20000.000
&gt;     &gt; &gt;
&gt;     &gt; &gt; dt                       = 0.002    ;  2 fs
&gt;     &gt; &gt;
&gt;     &gt; &gt; nsteps                   = 5000000  ; 10,000 ps
&gt;     &gt; &gt;
&gt;     &gt; &gt; comm-mode                = Linear
&gt;     &gt; &gt;
&gt;     &gt; &gt; nstcomm                  = 10
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; emtol                    = 0.001
&gt;     &gt; &gt;
&gt;     &gt; &gt; niter                    = 30
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; ; OUTPUT CONTROL OPTIONS
&gt;     &gt; &gt;
&gt;     &gt; &gt; nstxout                  = 5000    ; save coordinates every 8 ps
&gt;     &gt; &gt;
&gt;     &gt; &gt; nstvout                  = 5000    ; save velocities every 8 ps
&gt;     &gt; &gt;
&gt;     &gt; &gt; nstfout                  = 0
&gt;     &gt; &gt;
&gt;     &gt; &gt; nstlog                   = 5000    ; update log file every 8 ps
&gt;     &gt; &gt;
&gt;     &gt; &gt; nstenergy                = 5000    ; save energies every 8 ps
&gt;     &gt; &gt;
&gt;     &gt; &gt; ;nstxtcout is replaced with nstxout-compressed
&gt;     &gt; &gt;
&gt;     &gt; &gt; nstxout-compressed       = 5000    ; no. of steps between writing
&gt;     &gt; &gt; coordinates using lossy compression
&gt;     &gt; &gt;
&gt;     &gt; &gt; ;xtc-precision is replaced with compressed-x-precision
&gt;     &gt; &gt;
&gt;     &gt; &gt; compressed-x-precision   = 3000    ; precision with which to write to the
&gt;     &gt; &gt; compressed trajectory file
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; ; Periodic boundary conditions: xyz (default), no (vacuum)
&gt;     &gt; &gt;
&gt;     &gt; &gt; cutoff-scheme            = Verlet ; default for GROMACS-5.0
&gt;     &gt; &gt;
&gt;     &gt; &gt; nstlist                  = 1      ; A smaller nstlist may reduce LINCS
&gt;     &gt; &gt; warnings
&gt;     &gt; &gt;
&gt;     &gt; &gt; ns-type                  = grid   ;
&gt;     &gt; &gt;
&gt;     &gt; &gt; pbc                      = xyz
&gt;     &gt; &gt;
&gt;     &gt; &gt; ;periodic_molecules       = yes  ; cannot be used with SHAKE
&gt;     &gt; &gt;
&gt;     &gt; &gt; ; nblist cut-off
&gt;     &gt; &gt;
&gt;     &gt; &gt; rlist                    = 1.40   ; short-range neighbourlist cutoff
&gt;     &gt; &gt;
&gt;     &gt; &gt; nstcalclr                = 1
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; coulombtype              = PME
&gt;     &gt; &gt;
&gt;     &gt; &gt; coulomb-modifier         = Potential-shift-Verlet
&gt;     &gt; &gt;
&gt;     &gt; &gt; r_coulomb                = 1.40  ; short-range electrostatic cutoff
&gt;     &gt; &gt;
&gt;     &gt; &gt; rcoulomb-switch          = 0     ; used when potential-switch
&gt;     &gt; &gt; (coulomb-modifier) is specified
&gt;     &gt; &gt;
&gt;     &gt; &gt; epsilon-r                = 1     ; default value=1
&gt;     &gt; &gt;
&gt;     &gt; &gt; pme_order                = 4
&gt;     &gt; &gt;
&gt;     &gt; &gt; fourierspacing           = 0.12
&gt;     &gt; &gt;
&gt;     &gt; &gt; fourier_nx               = 0
&gt;     &gt; &gt;
&gt;     &gt; &gt; fourier_ny               = 0
&gt;     &gt; &gt;
&gt;     &gt; &gt; fourier_nz               = 0
&gt;     &gt; &gt;
&gt;     &gt; &gt; ewald_rtol               = 1e-05
&gt;     &gt; &gt;
&gt;     &gt; &gt; ewald_geometry           = 3d
&gt;     &gt; &gt;
&gt;     &gt; &gt; epsilon_surface          = 0      ; Use 0 when simulation is done with PME
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; ;vdw-type = switch is replaced with the following two flags in GROMACS-5.0
&gt;     &gt; &gt;
&gt;     &gt; &gt; vdwtype                  = Cut-off
&gt;     &gt; &gt;
&gt;     &gt; &gt; vdw_modifier             = Potential-switch
&gt;     &gt; &gt;
&gt;     &gt; &gt; ; cut-off lengths
&gt;     &gt; &gt;
&gt;     &gt; &gt; rvdw-switch              = 1.00
&gt;     &gt; &gt; rvdw                     = 1.40   ; short-range van der Waals cutoff
&gt;     &gt; &gt;
&gt;     &gt; &gt; constraints              = none
&gt;     &gt; &gt; constraint-algorithm     = SHAKE
&gt;     &gt; &gt; continuation             = yes    ; starting from NVT run
&gt;     &gt; &gt; Shake-SOR                = no
&gt;     &gt; &gt; shake-tol                = 0.0001
&gt;     &gt; &gt; ;lincs-order              = 4     ; accuracy of LINCS
&gt;     &gt; &gt; ;lincs-iter               = 24    ; also related to accuracy
&gt;     &gt; &gt; ; Lincs will write a warning to the stderr if in one step a bond
&gt;     &gt; &gt; ; rotates over more degrees than
&gt;     &gt; &gt; ;lincs-warnangle          = 30
&gt;     &gt; &gt;
&gt;     &gt; &gt; Tcoupl                   = v-rescale
&gt;     &gt; &gt; nh-chain-length          = 1    ; the leapfrog md integrator only supports
&gt;     &gt; &gt; 1
&gt;     &gt; &gt; ; Groups to couple separately
&gt;     &gt; &gt; tc-grps                  = GRA SOL
&gt;     &gt; &gt; ; Time constant (ps) and reference temperature (K)
&gt;     &gt; &gt; tau-t                    = 0.5 0.5
&gt;     &gt; &gt; ref-t                    = 300 300
&gt;     &gt; &gt; nsttcouple               = 20  ; added by k-wong ; not specified in
&gt;     &gt; &gt; Yanbin's MDP
&gt;     &gt; &gt;
&gt;     &gt; &gt; Pcoupl                   = Parrinello-Rahman
&gt;     &gt; &gt; Pcoupltype               = isotropic
&gt;     &gt; &gt; ; Time constant (ps), compressibility (1/bar) and reference P (bar)
&gt;     &gt; &gt; tau-p                    = 5.0  ; time constant in ps
&gt;     &gt; &gt; compressibility          = 4.5e-5  ; isothermal compressibility of water,
&gt;     &gt; &gt; bar^-1
&gt;     &gt; &gt; ref-p                    = 1.0  ; reference pressure, in bar
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; ; Non-equilibrium MD stuff
&gt;     &gt; &gt; freezegrps               = GRA
&gt;     &gt; &gt; freezedim                = Y Y Y
&gt;     &gt; &gt;
&gt;     &gt; &gt; gen_vel                  = no
&gt;     &gt; &gt; ;gen_temp                 = 300
&gt;     &gt; &gt; ;gen_seed                 = 173529
&gt;     &gt; &gt;
&gt;     &gt; &gt;
&gt;     &gt; &gt; Thanks in advance!

&gt;     &gt; &gt; Regards,
&gt;     &gt; &gt;
&gt;     &gt; &gt; Kester

</pre></div></div></blockquote></p>
<img src='http://mail.ibs.re.kr/checkread/ODg2NjMy/Z214LXVzZXJzQGdyb21hY3Mub3Jn/' width='1px' height='1px' />