Sure, worth trying. Mark
On Tue, Jun 25, 2013 at 10:07 AM, Broadbent, Richard <[email protected]> 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" <[email protected]> wrote: > >>On Tue, Jun 25, 2013 at 1:34 AM, S. Alireza Bagherzadeh >><[email protected]> 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 [email protected] >>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 [email protected]. >>* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > -- > gmx-users mailing list [email protected] > 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 [email protected]. > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing list [email protected] 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 [email protected]. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
