Well, when I set the box, I used 'editconf' command to rescale the box to have the right density which was ~997. After simulation, I got the following:
Energy Average RMSD Fluct. Drift Tot-Drift ---------------------------------------------------------------------------- --- Density (SI) 956.765 150.266 130.596 0.257477 257.477 Payman -----Original Message----- From: [email protected] [mailto:[email protected]] On Behalf Of David van der Spoel Sent: June 8, 2009 11:29 AM To: Discussion list for GROMACS users Subject: Re: [gmx-users] Energies in simulation Payman Pirzadeh wrote: > Here is the content of .itp file which I developed: > > ; This is an itp file to describe water's six-site model by H. Nada and J.P. > J. M. Van der Eerden, J. Chem. Phys. Vol.118, no.16, pp7401-7413 (2003) > ; This model is a combination of TIP4P and TIP5P. It has three LJ sites and > 3 Coulomb sites > ; O-H bond length is 0.980A, HOH angle is 108.00degrees, LOL angle is 111.00 > degrees, O-M and O-L are about 0.230A and 0.8892A respectively > > [ defaults ] > ; non-bondedtype combrule genpairs FudgeLJ > FudgeQQ N > 1 2 NO > > [ atomtypes ] > ;name mass charge ptype c6 c12 > OW 15.9994 0.0 A 0.3115 0.714845562 > HW 1.00800 0.477 A 0.0673 0.11541 > MW 0.000 -0.866 D 0.00 0.00 > LW 0.00 -0.044 D 0.00 0.00 > > [ moleculetype ] > ;molname nrexcl > SOL 1 > > [ atoms ] > ; nr atomtype resnr residuename atom cgnr charge > 1 OW 1 SOL OW 1 0.0 > 2 HW 1 SOL HW1 1 0.477 > 3 HW 1 SOL HW2 1 0.477 > 4 MW 1 SOL MW 1 -0.866 > 5 LW 1 SOL LP1 1 -0.044 > 6 LW 1 SOL LP2 1 -0.044 > > [ settles ] > ; OW function doh dhh > 1 1 0.0980 0.15856 > > [ dummies3 ] > ; These set of parameters are for M site which can be easily calculated > using TIP4P calculations from tip4p.itp > ; So, it will be described as dummy site 3: r(v)= r(i) + a*(r(i)-r(j)) + > b*(r(i)-r(k)) > ; const = |OH|/{|OH|*cos(HOH/2)} => Due to vector algebra a=b=const/2. > Remember that OM is in the same direction of OH bonds. > ; Remember this site is in the same plane of OH bonds; so, its function 1 > ; > ; site from function a b > 4 1 2 3 1 0.199642536 0.199642536 > > ; Now we define the position of L sites which can be obtained from tip5p.itp > ; So, it will be described as dummy site 3out: r(v) = r(i) + a*(r(i)-r(j)) + > b*(r(i)-r(k)) + c*(r(ij)Xr(ik)) > ; const1 = {|OL|*cos(LOL/2)}/{|OH|*cos(HOH/2)} => Due to vector algebra > |a|=|b|=const/2. since the lone pairs are in opposite direction of OH bonds, > a minus sign is added. This part is similar to M site. > ; const2 = {|OL|*sin(LOL/2)}/{|OH|*|OH|*sin(HOH)} => The denominator is the > magnitude of vector product of OH bonds. > ; This sites are tetrahedral sites; so, its function 4 > ; > ; site from function a b c > 5 1 2 3 4 -0.437172388 -0.437172388 > 8.022961206 > 6 1 2 3 4 -0.437172388 -0.437172388 > -8.022961206 > > [ exclusions ] > 1 2 3 4 5 6 > 2 1 3 4 5 6 > 3 1 2 4 5 6 > 4 1 2 3 5 6 > 5 1 2 3 4 6 > 6 1 2 3 4 5 > > > And here is the mdp file which I used for the simulation run: > cpp = cpp > include = -I../top > define = > > ; Run control > > integrator = md > dt = 0.001 ;1 fs > nsteps = 1000000 ;10 ns > comm_mode = linear > nstcomm = 1 > > ;Output control > > nstxout = 5000 > nstlog = 5000 > nstenergy = 5000 > nstxtcout = 1000 > xtc_grps = > energygrps = > > ; Neighbour Searching > > nstlist = 10 > ns_type = grid > rlist = 0.9 > > ; Electrostatistics > > coulombtype = PME > rcoulomb = 0.9 > epsilon_r = 1 > > ; Vdw > > vdwtype = cut-off > rvdw = 1.2 > DispCorr = EnerPres > > ;Ewald > > fourierspacing = 0.12 > pme_order = 4 > ewald_rtol = 1e-6 > optimize_fft = yes > > ; Temperature coupling > > tcoupl = Berendsen > tc-grps = System > tau_t = 0.1 > ref_t = 300 > > ; Pressure Coupling > > Pcoupl = Berendsen > Pcoupltype = isotropic > tau_p = 1.0 > compressibility = 5.5e-5 > ref_p = 1.0 > gen_vel = yes > > > The expected Potential energy is supposed to be around -41.5kJ/mol while my > potential is around -22.2kJ/mol. I calculated the energies by g_energy > command. > And do yo have the right density? > Payman > > > -----Original Message----- > From: [email protected] [mailto:[email protected]] > On Behalf Of David van der Spoel > Sent: June 8, 2009 11:13 AM > To: Discussion list for GROMACS users > Subject: Re: [gmx-users] Energies in simulation > > Payman Pirzadeh wrote: >> To the best of my knowledge no, but how can I check that? >> > A. read the original paper: is your topology correct? Are the simulation > parameters the same? > > B. post the itp file here and mdp file and specify energy and expected > energy. How about energy units? > >> -----Original Message----- >> From: [email protected] [mailto:[email protected]] >> On Behalf Of David van der Spoel >> Sent: June 8, 2009 11:06 AM >> To: Discussion list for GROMACS users >> Subject: Re: [gmx-users] Energies in simulation >> >> Payman Pirzadeh wrote: >>> Hi, >>> >>> I am using my own water model which I developed its .itp file. When >>> simulation is done after 1ns and energy is kinetic and potential >>> energies are analyzed, the kinetic value is almost OK, but the potential >>> energy is almost half of the value reported in literature and another MD >>> code that I am currently using. I double-checked the parameters I gave >>> in the .itp with TIP4P and TIP5P to make sure everything is correct in >>> format and unit. But I can not figure out the problem. Any ideas? >>> >> Is there any self-energy involved (i.e. a monomer energy that yo have to >> subtract)? >>> >>> >>> Payman >>> >>> >>> >>> >>> >>> >>> ------------------------------------------------------------------------ >>> >>> _______________________________________________ >>> gmx-users mailing list [email protected] >>> http://lists.gromacs.org/mailman/listinfo/gmx-users >>> Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php >> > > -- David. ________________________________________________________________________ David van der Spoel, PhD, Professor of Biology Dept. of Cell and Molecular Biology, Uppsala University. Husargatan 3, Box 596, 75124 Uppsala, Sweden phone: 46 18 471 4205 fax: 46 18 511 755 [email protected] [email protected] http://folding.bmc.uu.se ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ _______________________________________________ gmx-users mailing list [email protected] http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php _______________________________________________ gmx-users mailing list [email protected] http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/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/mailing_lists/users.php
