Dear Mark and Richard, Thank you so much for your help.
I used shake tolerance of 1e-10 and it solved the issue. I am getting very good energy conservation. I would have never guessed that shake tolerance has such a huge impact! I also experimented LINC with the following options: constraint-algorithm = lincs continuation = no ; DO apply constraints to the start configuration lincs-order = 6 lincs-iter = 2 lincs-warnangle = 30 I was able to stop the temperature from decreasing however the there was again substantial total energy drift. Many thanks, Alireza > Message: 5 > Date: Tue, 25 Jun 2013 13:55:08 +0200 > From: Mark Abraham <mark.j.abra...@gmail.com> > Subject: Re: [gmx-users] Settles vs. 3 normal constraints (energy > conservation problem) > To: Discussion list for GROMACS users <gmx-users@gromacs.org> > Message-ID: > < > camnumaro6hv6+u+p6pvmckbfqpihb3iwtt_1hdajerm_hye...@mail.gmail.com> > Content-Type: text/plain; charset=ISO-8859-1 > > Sure, worth trying. > > Mark > > On Tue, Jun 25, 2013 at 10:07 AM, Broadbent, Richard > <richard.broadben...@imperial.ac.uk> wrote: > > For my system reducing shake-tol greatly improves the energy conservation > > generally 1.0e-7 is the largest I would use. However if you want very > good > > energy conservation 1.0e-9 or lower might be needed. > > > > This effect might only be for my system but I think it might help here > too > > > > Richard > > > > On 25/06/2013 08:45, "Mark Abraham" <mark.j.abra...@gmail.com> wrote: > > > >>On Tue, Jun 25, 2013 at 1:34 AM, S. Alireza Bagherzadeh > >><s.a.bagherzade...@gmail.com> wrote: > >>> Dear All, > >>> > >>> I have a box of 3073 tip4p water molecules. I do a 250ps nvt, then 250 > >>>ps > >>> npt and finally a 1 ns nve (production run). > >>> > >>> I used the opls forcefield and I copied the tip4p.itp to my working > >>> directory (in order to be able to make changes). > >>> > >>> In one case I used the [ settles ] directive to constraint water > >>>molecules. > >>> .top file: > >>> > >>> ; Topology for 3087 TIPT4P waters > >>> > >>> #include > >>> > > >>>"/global/software/gromacs/gromacs-4.5.5/share/gromacs/top/oplsaa.ff/force > >>>field.itp" > >>> ;SOL > >>> ;----------------------------------------------- > >>> ; water topology - liquid phase > >>> #include "./tip4p.itp" > >>> ;------------------------------------------------ > >>> [ system ] > >>> ; Name > >>> A box of 216 tip4p for protocol testing > >>> > >>> [ molecules ] > >>> ; Compound #mols > >>> SOL 3073 > >>> > >>> > >>> .itp file: > >>> ; Note the strange order of atoms to make it faster in gromacs. > >>> ; > >>> [ moleculetype ] > >>> ; molname nrexcl > >>> SOL 2 > >>> > >>> [ atoms ] > >>> ; id at type res nr residu name at name cg nr charge > >>> 1 opls_113 1 SOL OW 1 0.0 > >>> 2 opls_114 1 SOL HW1 1 0.52 > >>> 3 opls_114 1 SOL HW2 1 0.52 > >>> 4 opls_115 1 SOL MW 1 -1.04 > >>> > >>> #ifndef FLEXIBLE > >>> [ settles ] > >>> ; OW funct doh dhh > >>> 1 1 0.09572 0.15139 > >>> #else > >>> [ bonds ] > >>> ; i j funct length force.c. > >>> 1 2 1 0.09572 502416.0 0.09572 502416.0 > >>> 1 3 1 0.09572 502416.0 0.09572 502416.0 > >>> > >>> [ angles ] > >>> ; i j k funct angle force.c. > >>> 2 1 3 1 104.52 628.02 104.52 628.02 > >>> #endif > >>> > >>> [ exclusions ] > >>> 1 2 3 4 > >>> 2 1 3 4 > >>> 3 1 2 4 > >>> 4 1 2 3 > >>> > >>> ; The position of the virtual site is computed as follows: > >>> ; > >>> ; O > >>> ; > >>> ; D > >>> ; > >>> ; H H > >>> ; > >>> ; const = distance (OD) / [ cos (angle(DOH)) * distance (OH) ] > >>> ; 0.015 nm / [ cos (52.26 deg) * 0.09572 nm ] > >>> > >>> ; Vsite pos x4 = x1 + a*(x2-x1) + b*(x3-x1) > >>> > >>> [ virtual_sites3 ] > >>> ; Vsite from funct a b > >>> 4 1 2 3 1 0.128012065 0.128012065 > >>> > >>> .mdp file: > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> ; Run control > >>> integrator = md ; Leap-frog algorithm > >>> tinit = 0 ; starting time [ps] > >>> dt = 0.001 ; time step [ps] > >>> nsteps = 250000 ; number of steps > >>> nstcomm = 100 ; frequency for center of mass > >>>motion > >>> removal [steps] > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> ; Output control > >>> nstxout = 0 ; frequency to write coordinates to > >>> output trajectory file [steps] > >>> nstvout = 0 ; frequency to write velocities to > >>> output trajectory file [steps] > >>> nstfout = 0 ; frequency to write forces to > >>>output > >>> trajectory file [steps] > >>> nstlog = 500 ; frequency to write energies to > log > >>> file [steps] > >>> nstenergy = 500 ; frequency to write energies to > >>>energy > >>> file [steps] > >>> nstxtcout = 000 ; frequency to write coordinates to > >>>xtc > >>> trajectory [steps] > >>> xtc-precision = 1000 ; precision to write to xtc > >>>trajectory > >>> [real] > >>> xtc_grps = SOL > >>> energygrps = SOL > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> ; Neighborsearching and short-range nonbonded interactions > >>> nstlist = 10 ; frequency to update the neighbor > >>>list > >>> [steps] > >>> ns_type = grid ; (grid / simple) search for > >>> neighboring list > >>> pbc = xyz ; priodic boundary conditions (xyz > >>>/ no > >>> / xy) > >>> rlist = 1.7 ; cut-off distance for the > >>>short-range > >>> neighbor list [nm] > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> ; Electrostatics > >>> coulombtype = PME-Switch > >>> rcoulomb_switch = 1.3 ; where to switch the Coulomb > >>>potential > >>> [nm] > >>> rcoulomb = 1.5 ; distance for the Coulomb cut-off > >>>[nm] > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> ; van der Waals > >>> vdw-type = shift ; (the LJ normal out at rvdw_switch > >>>to > >>> reach zero at rvdw) > >>> rvdw-switch = 1.3 ; where to strat switching the LJ > >>> potential [nm] > >>> rvdw = 1.5 ; cut-off distance for vdw > potenrial > >>> [nm] > >>> DispCorr = EnerPres ; (Apply long range dispersion > >>> corrections for Energy and Pressure) > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> ; EWALD/PME/PPPM parameters > >>> fourierspacing = 0.12 ; maximum grid spacing for the FFT > >>>when > >>> using PPPM or PME [nm] > >>> pme_order = 6 ; interpolation order for PME > >>> ewald_rtol = 1e-06 ; relative strength of the > >>> Ewald-shifted direct potential at rcoulomb > >>> epsilon_surface = 0 ; dipole correction to the Ewald > >>> summation in 3d > >>> optimize_fft = no ; optimal FFT plan for the grid at > >>> startup (yes / no) > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> ; Temperature coupling > >>> tcoupl = no > >>> ;tcoupl = nose-hoover > >>> nh-chain-length = 1 > >>> tc_grps = system > >>> tau_t = 0.2 ; time constatn for coupling [ps], > 1 > >>> for each group > >>> ref_t = 300 ; refernce temperature [K], 1 for > >>>each > >>> group > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> ; Pressure coupling > >>> Pcoupl = no > >>> ;Pcoupl = Parrinello-Rahman > >>> tau_p = 0.5 ; time constant for coupling > >>> compressibility = 4.5e-05 > >>> ref_p = 1.0 ; reference pressure for coupling > >>>[bar] > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> ; Velocity generation > >>> gen_vel = no ; generating velocities according to > >>> maxwell distribution? > >>> gen_temp = 300 ; [K] > >>> gen_seed = -1 ; initialize random generator based > >>>on > >>> the process ID number [integer] > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> ; Bonds > >>> constraints = all-angles > >>> constraint-algorithm = shake > >>> continuation = no ; DO apply constraints to the start > >>> configuration > >>> shake-tol = 1e-5 > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > > >>>;------------------------------------------------------------------------ > >>>-----------------------------; > >>> > >>> This case I obtain a very good energy conservation. > >>> > >>> However, in another case where I change settles to 3 normal > constraints, > >>> there is a rather substantial energy drift (~2% per nano second) and > >>>also > >>> the T does not remain constant and it decreases (~8 K per ns). > >> > >>Hmm, that doesn't sound good. Guessing wildly, SHAKE and the v-site > >>are not properly implemented? > >> > >>> .itp file: > >>> ; Note the strange order of atoms to make it faster in gromacs. > >>> ; > >>> [ moleculetype ] > >>> ; molname nrexcl > >>> SOL 2 > >>> > >>> [ atoms ] > >>> ; id at type res nr residu name at name cg nr charge > >>> 1 opls_113 1 SOL OW 1 0.0 > >>> 2 opls_114 1 SOL HW1 1 0.52 > >>> 3 opls_114 1 SOL HW2 1 0.52 > >>> 4 opls_115 1 SOL MW 1 -1.04 > >>> > >>> #ifndef FLEXIBLE > >>> [ constraints ] > >>> ; ai aj funct b0 > >>> 1 2 1 0.09572 > >>> 1 3 1 0.09572 > >>> 2 3 1 0.15139 > >>> > >>> #else > >>> [ bonds ] > >>> ; i j funct length force.c. > >>> 1 2 1 0.09572 502416.0 0.09572 502416.0 > >>> 1 3 1 0.09572 502416.0 0.09572 502416.0 > >>> > >>> [ angles ] > >>> ; i j k funct angle force.c. > >>> 2 1 3 1 104.52 628.02 104.52 628.02 > >>> #endif > >>> > >>> [ exclusions ] > >>> 1 2 3 4 > >>> 2 1 3 4 > >>> 3 1 2 4 > >>> 4 1 2 3 > >>> > >>> ; The position of the virtual site is computed as follows: > >>> ; > >>> ; O > >>> ; > >>> ; D > >>> ; > >>> ; H H > >>> ; > >>> ; const = distance (OD) / [ cos (angle(DOH)) * distance (OH) ] > >>> ; 0.015 nm / [ cos (52.26 deg) * 0.09572 nm ] > >>> > >>> ; Vsite pos x4 = x1 + a*(x2-x1) + b*(x3-x1) > >>> > >>> [ virtual_sites3 ] > >>> ; Vsite from funct a b > >>> 4 1 2 3 1 0.128012065 0.128012065 > >>> > >>> > >>> The .mdp and .top files are the same. > >>> > >>> > >>> I was wondering if there is any difference between [settles] and 3 > >>>normal [ > >>> constraints]? > >> > >>I believe not, but Berk's probably the only expert on the > >>implementation, here... > >> > >>> * I am asking this question because in my real simulation, the above is > >>>a > >>> simplified version where I can check for energy conservation, I need to > >>>use > >>> water in two different molecule types and if I use [ settles ] for both > >>>of > >>> them I get this error: > >>> > >>> Fatal error: > >>> The [molecules] section of your topology specifies more than one block > >>>of > >>> a [moleculetype] with a [settles] block. Only one such is allowed. If > >>>you > >>> are trying to partition your solvent into different *groups* (e.g. for > >>> freezing, T-coupling, etc.) then you are using the wrong approach. > Index > >>> files specify groups. Otherwise, you may wish to change the least-used > >>> block of molecules with SETTLE constraints into 3 normal constraints. > >>> > >>> > >>> I appreciate being advised on what I might be doing wrong. > >> > >>I can't see anything wrong. I would suggest trying LINCS, which is > >>more likely to be better tested with v-sites. Please report back, > >>either way! :-) > >> > >>Mark > >>-- > >>gmx-users mailing list gmx-users@gromacs.org > >>http://lists.gromacs.org/mailman/listinfo/gmx-users > >>* Please search the archive at > >>http://www.gromacs.org/Support/Mailing_Lists/Search before posting! > >>* Please don't post (un)subscribe requests to the list. Use the > >>www interface or send it to gmx-users-requ...@gromacs.org. > >>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > > > -- > > gmx-users mailing list gmx-users@gromacs.org > > http://lists.gromacs.org/mailman/listinfo/gmx-users > > * Please search the archive at > http://www.gromacs.org/Support/Mailing_Lists/Search before posting! > > * Please don't post (un)subscribe requests to the list. Use the > > www interface or send it to gmx-users-requ...@gromacs.org. > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > > ------------------------------ > > -- > gmx-users mailing list > gmx-users@gromacs.org > http://lists.gromacs.org/mailman/listinfo/gmx-users > Please search the archive at > http://www.gromacs.org/Support/Mailing_Lists/Search before posting! > > End of gmx-users Digest, Vol 110, Issue 128 > ******************************************* > -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
