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
