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/forcefield.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
