Well, I tried to equilibrate the system, and I think the program has included all the distortion that the system had been through. When I checked the trajectory, the box was highly distorted as the simulation started, then returned to its cubic shape. Also, I am a bit confused by the word "strange". Well, based on manual the rvdw>rlist and it is. Also, authors in the paper have tested U with various cutoff radii from 8-11A. The results were very close! To make sure we are on the same page, the paper address is J. Chem. Phys., Vol118, No.16 pp7401 (2003).
Payman -----Original Message----- From: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] On Behalf Of David van der Spoel Sent: June 8, 2009 11:58 AM To: Discussion list for GROMACS users Subject: Re: [gmx-users] Energies in simulation Payman Pirzadeh wrote: > 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 The RMSD is very large. Are you sure this is in equilibrium? It could be that your box is exploding. Please check as well that the strange twin range cut-off that you are using is what the original authors used. > > > > Payman > > -----Original Message----- > From: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] > 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: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] >> 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: gmx-users-boun...@gromacs.org > [mailto:gmx-users-boun...@gromacs.org] >>> 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 gmx-users@gromacs.org >>>> 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 gmx-users-requ...@gromacs.org. >>>> 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 sp...@xray.bmc.uu.se sp...@gromacs.org http://folding.bmc.uu.se ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ _______________________________________________ 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/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/mailing_lists/users.php _______________________________________________ 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/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/mailing_lists/users.php