Hello Justin, I used 0.1fs timestep and I could get dielectric constant around 56 which looks ok. To get the dielectric constant around ~75, should I reduce the timestep around 0.05fs.
I used temperature 298k for my simulation and average tmep. I got is 297.83 K. I am doing NPT and used 1.0 bar as reference pressure in my .mdp. I am getting average pressure around 2.58 bar. I run the equilibration for 1ns. Energy Average RMSD Fluct. Drift Tot-Drift ------------------------------------------------------------------------------- Temperature 297.83 10.9687 10.9643 -0.00165807 -1.07774 Pressure (bar) 2.57551 4246.51 4246.51 -0.0360116 -23.4075 Why there is change in pressure? Is there any wrong with simulation? Nilesh On Wed, May 11, 2011 4:18 pm, Justin A. Lemkul wrote: > > > Nilesh Dhumal wrote: > >> I have general question. If this parameters are not giving the proper >> results then How come I get the proper results with same parameter if I >> use 1/2 value of Kbond with the same parameters. >> > > Random chance, I imagine. Your .mdp settings are not correct and the > cutoffs are unreasonably short for any modern force field. If you're > seeking to reproduce others' work, then you have to use the same protocol. > You're not > doing that. > >> I checked the temp when I use half value of Kbond and it goes max. >> ~320. >> >> > > This would also indicate instability. You managed to get a "correct" > result for one observable from an incorrect simulation wherein the > thermodynamic observables do not correspond to what you have set. > > The way I found all of this out was first and foremost by watching the > trajectory. After about 30 ps, the water molecules freeze in place and > vibrate as if they are about to shear apart. At this exact time, the > value of epsilon (from epsilon.xvg, output by g_dipoles) plummets down, > after having steadily increased for the first 30 ps, as expected. The > temperature spikes and the pressure begins to oscillate by +/- 10^5, far > higher than what it should for a system of a few hundred water molecules. > > Watch your trajectories, make sure you're getting the right average > temperature and pressure, and that all your other energetic terms > stabilize, then analyze the result of the simulation. Otherwise, any > inferred output is meaningless. > > -Justin > > >> Nilesh >> >> >> On Wed, May 11, 2011 4:03 pm, Justin A. Lemkul wrote: >> >> >>> Nilesh Dhumal wrote: >>> >>> >>>> Hello Justin, >>>> >>>> >>>> >>>> I used 0.5 fs time step and still I got dielectric constant ~2. >>>> >>>> >>>> >>>> This is the md.mdp file I used. I checked the temperature it >>>> didn't go more than ~ 320 K. >>>> >>> Still, that's unacceptable if your ref_t is 295. The fact that the >>> temperature is still rocketing up suggests the same instability I >>> found. >>> >>>> title = cpeptide MD cpp = >>>> /usr/bin/cpp >>>> integrator = md dt = 0.001 ; ps ! >>>> nsteps = 6500000 ; total 5 ps. >>>> nstcomm = 1000 nstxout = 1000 nstvout = >>>> 1000 >>>> nstfout = 1 nstlist = 1 ns_type >>>> = >>>> grid rlist = 0.6 rcoulomb = 0.6 rvdw = >>>> 0.7 >>>> >>>> >>> These cutoffs make no sense. The paper you've cited used 0.9 nm. If >>> you're trying to reproduce previous work, use the same settings. >>> >>>> coulombtype = PME vdwtype = cut-off >>> Similarly here, you haven't applied dispersion correction, which you >>> should, per the methods used in the JCP article. >>> >>> -Justin >>> >>> >>> >>>> pbc = xyz fourierspacing = 0.12 fourier_nx = 0 >>>> fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = >>>> 1e-5 >>>> optimize_fft = yes ; Berendsen temperature coupling is on >>>> Tcoupl = v-rescale >>>> tau_t = 0.1 tc-grps =system ref_t = 295 ; Pressure coupling is >>>> on Pcoupl = parrinello-rahman >>>> pcoupltype = isotropic tau_p = 0.5 >>>> compressibility = 4.5e-5 >>>> ref_p = 1.0 ; Generate velocites is on at 300 K. >>>> gen_vel = yes gen_temp = 295.0 gen_seed = >>>> 173529 >>>> >>>> >>>> >>>> >>>> Nilesh >>>> >>>> >>>> >>>> >>>> On Tue, May 10, 2011 6:14 pm, Justin A. Lemkul wrote: >>>> >>>> >>>> >>>>> The gmx-developers list is not the forum for these types of >>>>> questions. I am replying via gmx-users, which is where the >>>>> discussion should stay. >>>>> >>>>> I took a few minutes to dig into this. My conclusion is that >>>>> your system is not stable. I would encourage you to analyze the >>>>> temperature and pressure of your systems that are giving odd >>>>> results. When I used the >>>>> topology you provided, the temperature spiked to over 600 K and >>>>> the box began to oscillate extensively as if the system were about >>>>> to explode. The resulting epsilon value was about 1.6. >>>>> >>>>> >>>>> >>>>> If I reduce the timestep to 0.5 fs, the results are much more >>>>> reasonable, and using the Ka and Kb values in the paper (not >>>>> halved!) I get a much >>>>> more sensible result, even after just 100 ps. It looks to me as >>>>> if the combination of these parameters and a 1-fs timestep is not >>>>> entirely stable. I know the original authors used dt = 1 fs, but >>>>> the point remains. >>>>> >>>>> -Justin >>>>> >>>>> >>>>> >>>>> >>>>> Nilesh Dhumal wrote: >>>>> >>>>> >>>>> >>>>>> Hello, >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> I am using flexible water model for my system. I am referring a >>>>>> paper J. Chem. Phys. 124, 024503 (2006). Author have used Amber >>>>>> type force field. >>>>>> >>>>>> >>>>>> i.e. 1/2 Kbond (r-req)**2 + 1/2 Kangle(theta-thetaeq)**2. >>>>>> >>>>>> Kbond = 443153.3808 kJ/mol nm**2 >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> Kangle= 317.5656 kJ/mol rad**2. >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> I am using olss-aa force field parameters in Gromacs VERSION >>>>>> 4.0.7. >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> I checked some papers in which author have used oplsaa force >>>>>> field in Gromacs. 1/2 factor is not in opls force field if I >>>>>> compare opls and amber. >>>>>> >>>>>> >>>>>> I didnt get the proper dielectric constant for water when I >>>>>> used the parameters reported in paper >>>>>> >>>>>> (Kbond = 443153.3808 kJ/mol nm**2,Kangle= 317.5656 kJ/mol >>>>>> rad**2) >>>>>> >>>>>> >>>>>> >>>>>> I half the value of Kb and get the proper dielectric constant >>>>>> (~80) >>>>>> for water reported in paper. If I half Kangle then I dont get >>>>>> proper value. Below are the results for the dielectric constant >>>>>> of water. Bond length is nm. >>>>>> >>>>>> >>>>>> >>>>>> Here I have done some analysis. The original value reported in >>>>>> J. >>>>>> Chem. >>>>>> Phys. 124, 024503 2006, paper are >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> Kbond = 443153.3808 kJ/mol nm**2 >>>>>> Kangle = 317.5656 kJ/mol rad**2. >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> bond length Kbond angle Kangle dielectric >>>>>> constant 0.1012 >>>>>> 443153.3808 113.24 317.5656 ~1.9 : orginal value >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> 0.1012 221576.6904 113.24 317.5656 ~80 : 1/2 >>>>>> (Kbond) >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> 0.1012 443153.3808 113.24 158.7828 ~1.58 : 1/2 >>>>>> (kangle) >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> 0.1012 221576.6904 113.24 317.5656 ~1.9 : 1/2 >>>>>> (Kbond)&(Kangle) >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> Here I pasted spc_fw.itp file used for the simulation. >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> ; Derived from parsing of runfiles/alat.top.orig >>>>>> [ defaults ] >>>>>> ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ >>>>>> ;1 3 yes 0.5 0.5 >>>>>> ; comb-rule 3 is square-root sigma, the OPLSAA version >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> [ atomtypes ] >>>>>> ; full atom descriptions are available in ffoplsaa.atp >>>>>> ; name bond_type mass charge ptype sigma >>>>>> epsilon opls_1001 OW 8 15.99940 -0.820 A >>>>>> 3.1655e-01 >>>>>> 6.503e-01 >>>>>> opls_1002 HW 1 1.00800 0.410 A 0.00e+00 >>>>>> 0.00e+00 >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> [ bondtypes ] >>>>>> ; i j func b0 kb >>>>>> OW HW 1 0.1012 443153.3808 ; J. Chem. Phys. >>>>>> (2006),124,024503 >>>>>> [ angletypes ] >>>>>> ; i j k func th0 cth >>>>>> HW OW HW 1 113.24 317.5656 ; J. Chem. Phys. >>>>>> (2006),124,024503 >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> [ moleculetype ] >>>>>> ; Name nrexcl >>>>>> WAT 3 >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> [ atoms ] >>>>>> ; nr type resnr residue atom cgnr charge >>>>>> mass typeB chargeB massB 1 opls_1001 1 WAT OW >>>>>> 1 >>>>>> -0.82 >>>>>> 15.99940 ; >>>>>> 2 opls_1002 1 WAT HW1 1 0.41 1.008 >>>>>> ; >>>>>> 3 opls_1002 1 WAT HW2 1 0.41 1.008 >>>>>> ; >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> [ bonds ] >>>>>> ; i j funct >>>>>> 1 2 1 >>>>>> 1 3 1 >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> [ angles ] >>>>>> ; i j k funct >>>>>> 2 1 3 1 >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> In water.top file, I included spc_fw.itp file. >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> ; Include water topology >>>>>> #include "spc_fw.itp" >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> I run the simulation 6.5 ns for collecting data and I have >>>>>> total 256 >>>>>> water molecules. >>>>>> >>>>>> >>>>>> Is there anything wrong with my .itp file? >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> I am using Gromacs VERSION 4.0.7. >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> Thanks >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> Nilesh >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>> -- >>>>> ======================================== >>>>> >>>>> >>>>> >>>>> >>>>> Justin A. Lemkul >>>>> Ph.D. Candidate >>>>> ICTAS Doctoral Scholar >>>>> MILES-IGERT Trainee >>>>> Department of Biochemistry >>>>> Virginia Tech >>>>> Blacksburg, VA >>>>> jalemkul[at]vt.edu | (540) 231-9080 >>>>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin >>>>> >>>>> >>>>> >>>>> >>>>> ======================================== >>>>> -- >>>>> 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/Support/Mailing_Lists/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/Support/Mailing_Lists >>>>> >>>>> >>>>> >>>>> >>>>> >>>> >>>> >>> -- >>> ======================================== >>> >>> >>> >>> Justin A. Lemkul >>> Ph.D. Candidate >>> ICTAS Doctoral Scholar >>> MILES-IGERT Trainee >>> Department of Biochemistry >>> Virginia Tech >>> Blacksburg, VA >>> jalemkul[at]vt.edu | (540) 231-9080 >>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin >>> >>> >>> >>> ======================================== >>> -- >>> 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/Support/Mailing_Lists/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/Support/Mailing_Lists >>> >>> >>> >>> >> >> >> > > -- > ======================================== > > > Justin A. Lemkul > Ph.D. Candidate > ICTAS Doctoral Scholar > MILES-IGERT Trainee > Department of Biochemistry > Virginia Tech > Blacksburg, VA > jalemkul[at]vt.edu | (540) 231-9080 > http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin > > > ======================================== > -- > 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/Support/Mailing_Lists/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/Support/Mailing_Lists > > > -- 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/Support/Mailing_Lists/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/Support/Mailing_Lists