Re: [gmx-users] g_hbond results inconsistent between v4.0.7 and v4.5.1
Hi All: Has there been further resolution on the difference between g_hbond in v 4.0 and 4.5 (I'm running 4.5.3)? For counting water-water (tip4p) h-bonds in a box of 506 waters I get the following nhbdist output, which appears to neglect pbc since there are such a large number of unbound donor-H (and the other columns are consequently too low): @title Number of donor-H with N HBs @xaxis label Time (ps) @yaxis label N @TYPE xy @ view 0.15, 0.15, 0.75, 0.85 @ legend on @ legend box on @ legend loctype view @ legend 0.78, 0.8 @ legend length 2 @ s0 legend 0 HBs @ s1 legend 1 HB @ s2 legend 2 HBs @ s3 legend 3 HBs @ s4 legend Total 0.0e+00 706 174 132 0 438 1.0e-01 708 163 140 1 446 2.0e-01 707 168 137 0 442 3.0e-01 712 161 136 3 442 4.0e-01 713 165 134 0 433 5.0e-01 709 169 133 1 438 6.0e-01 705 169 137 1 446 7.0e-01 699 169 142 2 459 8.0e-01 710 167 133 2 439 9.0e-01 700 189 122 1 436 1.0e+00 707 169 136 0 441 1.1e+00 705 163 144 0 451 1.2e+00 707 171 133 1 440 1.3e+00 709 158 145 0 448 1.4e+00 718 155 138 1 434 1.5e+00 711 164 137 0 438 1.6e+00 713 160 139 0 438 1.7e+00 702 181 129 0 439 1.8e+00 707 162 142 1 449 1.9e+00 706 167 139 0 445 2.0e+00 709 168 135 0 438 2.1e+00 708 167 136 1 442 2.2e+00 708 171 133 0 437 2.3e+00 709 162 140 1 445 2.4e+00 698 171 142 1 458 2.5e+00 701 177 133 1 446 2.6e+00 706 177 129 0 435 2.7e+00 708 178 126 0 430 2.8e+00 698 187 127 0 441 2.9e+00 700 190 120 2 436 3.0e+00 703 176 131 2 444 3.1e+00 702 181 127 2 441 3.2e+00 694 184 134 0 452 I've tried playing with different index assignments, and I get what appears to be a reasonable output for 0 HBs by making the two groups appear different, excluding the dummy atom from one. Although, because the program then double counts, the other columns don't convey the intended information. Robin -- Robin C. Underwood Chemistry Department 560 Oval Drive West Lafayette, IN 47907 -- gmx-users mailing listgmx-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] Note on oscillation period / time step
GMX Users: When I run grompp on methanol in water in v. 4.5.3 I get the following Note: NOTE 1 [file methanol.top, line 73]: The bond in molecule-type MET between atoms 5 OA and 6 HO has an estimated oscillational period of 9.0e-03 ps, which is less than 10 times the time step of 1.0e-03 ps. Maybe you forgot to change the constraints mdp option. Here is my mdp file: ; User Robin ; Wed Dec.1 2010 ; Input file ; title = Yo cpp = /lib/cpp integrator = md dt = 0.001; ps ! nsteps = 110 ; total 1000 ps=1 ns nstcomm = 10 ; (x) coordinates, (v) is velocities, (f) is forces ; This affects traj.trr ; these numbers shouldn't be too small (100,000), takes up disk space nstxout = 1000 nstvout = 1000 nstfout = 1000 ; Output frequency for energies nstlog = 100 nstenergy = 100 ; This affects .xtc file for movies this should be a fairly small number nstxtcout = 100 xtc-precision = 1000 ; parameters for neighbors nstlist = 10 ns_type = grid rlist = 0.9 ; electrostatics coulombtype = pme pme_order = 4 ewald_rtol = 1e-5 ewald_geometry = 3d epsilon_surface = 0 fourierspacing = 0.1 fourier_nx = 0.0 fourier_ny = 0.0 fourier_nz = 0.0 rcoulomb= 0.9 optimize_fft= yes Dispcorr= EnerPres vdwtype = cutoff rvdw_switch = 0.85 rvdw= 0.9 ; Berendsen temperature coupling is on in two groups Tcoupl = v_rescale tc-grps = system tau_t = 0.1 ref_t = 300 ; Energy monitoring energygrps = MET SOL ; Isotropic pressure coupling is now on Pcoupl = berendsen Pcoupltype = isotropic tau_p = 0.5 compressibility = 4.5e-5 ref_p = 1.0 ; Generate velocites is off at 300 K. gen_vel = no gen_temp= 300.0 gen_seed= 173529 Here is my top file: ; ; ; Include forcefield parameters #include ffoplsaa.itp ; ; Methanol, Jorgensen et al. JACS 118 pp. 11225 (1996) ; [ moleculetype ] ; name nrexcl MET 3 [ atoms ] ; nrtype resnr residuatomcgnrcharge mass 1 opls_157 1 METC 1 0.145 2 opls_156 1 METH 1 0.04 3 opls_156 1 METH 1 0.04 4 opls_156 1 METH 1 0.04 5 opls_154 1 MET OA 1 -0.683 6 opls_155 1 MET HO 1 0.418 [ bonds ] ; aiaj funct c0 c1 1 2 1 1 3 1 1 4 1 1 5 1 5 6 1 [ pairs ] ; i j funct 2 6 3 6 4 6 [ angles ] ; aiajak funct c0c1 ; H3 2 1 5 1 3 1 5 1 4 1 5 1 4 1 3 1 4 1 2 1 [ system ] ; Name MET in water [ molecules ] ; Compound#mols MET 1 SOL 503 Do I need to add some sort of constraint on my methanol molecule? Thanks, Robin -- Robin C. Underwood Chemistry Department 560 Oval Drive West Lafayette, IN 47907 -- gmx-users mailing listgmx-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] Re: Note on oscillation period / time step
Thanks for all the responses. After rechecking my top and manually inputting the [ ...types ] fields, the note remains. It is not clear to me whether I need to modify my top file or add some additional constraint (either in the top or mdp). Is it OK to have a bond oscillation period of this magnitude relative to my timestep? I've not seen this note in previous versions. Robin Javier Cerezo wrote: Hello. It looks like you're missing the bonded-interactions parameters. [ bonds ] ; aiaj funct *R0 Kb *12 1 * XX* 1 3 1 * XX* 1 4 1 * XX* 1 5 1 * XX* 5 6 1 * XX* (also angles...) Is it implicit somewhere? I don't know if it is automatic and if so, from where the parameters are taken. If it is by using atoms names and the ffbonded.itp from oplsaa, then atom name OA is missing in the itp file. I think that, according to the code (grompp.c), this note may arise if the force constants are zero. Atom types are interpolated from ffbonded.itp and assigned. If the problem were in the bonded parameters, a very different fatal error arises. There is no problem here. -Justin Javier El 02/12/10 17:11, Robin C. Underwood escribió: GMX Users: When I run grompp on methanol in water in v. 4.5.3 I get the following Note: NOTE 1 [file methanol.top, line 73]: The bond in molecule-type MET between atoms 5 OA and 6 HO has an estimated oscillational period of 9.0e-03 ps, which is less than 10 times the time step of 1.0e-03 ps. Maybe you forgot to change the constraints mdp option. Here is my mdp file: ;User Robin ;Wed Dec.1 2010 ;Input file ; title = Yo cpp = /lib/cpp integrator = md dt = 0.001 ; ps ! nsteps = 110 ; total 1000 ps=1 ns nstcomm = 10 ; (x) coordinates, (v) is velocities, (f) is forces ; This affects traj.trr ; these numbers shouldn't be too small (100,000), takes up disk space nstxout = 1000 nstvout = 1000 nstfout = 1000 ; Output frequency for energies nstlog = 100 nstenergy = 100 ; This affects .xtc file for movies this should be a fairly small number nstxtcout = 100 xtc-precision = 1000 ; parameters for neighbors nstlist = 10 ns_type = grid rlist = 0.9 ; electrostatics coulombtype = pme pme_order= 4 ewald_rtol = 1e-5 ewald_geometry = 3d epsilon_surface = 0 fourierspacing = 0.1 fourier_nx = 0.0 fourier_ny = 0.0 fourier_nz = 0.0 rcoulomb= 0.9 optimize_fft= yes Dispcorr = EnerPres vdwtype = cutoff rvdw_switch = 0.85 rvdw= 0.9 ; Berendsen temperature coupling is on in two groups Tcoupl = v_rescale tc-grps = system tau_t = 0.1 ref_t = 300 ; Energy monitoring energygrps = MET SOL ; Isotropic pressure coupling is now on Pcoupl = berendsen Pcoupltype = isotropic tau_p = 0.5 compressibility = 4.5e-5 ref_p = 1.0 ; Generate velocites is off at 300 K. gen_vel = no gen_temp= 300.0 gen_seed= 173529 Here is my top file: ; ; ; Include forcefield parameters #include ffoplsaa.itp ; ; Methanol, Jorgensen et al. JACS 118 pp. 11225 (1996) ; [ moleculetype ] ; name nrexcl MET 3 [ atoms ] ; nrtype resnr residuatomcgnrcharge mass 1 opls_157 1 METC 1 0.145 2 opls_156 1 METH 1 0.04 3 opls_156 1 METH 1 0.04 4 opls_156 1 METH 1 0.04 5 opls_154 1 MET OA 1 -0.683 6 opls_155 1 MET HO 1 0.418 [ bonds ] ; aiaj funct c0 c1 12 1 13 1 14 1 15 1 56 1 [ pairs ] ; i j funct 26 36 46 [ angles ] ; aiajak funct c0c1 ; H3 2 1 5 1 3 1 5 1 4 1 5 1 4 1 31 4 1 2 1 [ system ] ; Name MET in water [ molecules ] ; Compound#mols MET 1 SOL 503 Do I need to add some sort of constraint on my methanol molecule? Thanks, Robin -- Javier CEREZO BASTIDA Estudiante de Doctorado - Dpto. Química-Física Universidad de Murcia 30100 MURCIA (España) Tlf.(+34
[gmx-users] g_hbond modification
I would like to modify the g_hbond code (or preferably, know how to implement g_hbond if it is already capable) to count the number of non-hydrogen bound water OH donors in a simulation. I define a non-hydrogen bound donor as one that does not meet the distance requirement, or one that may meet the distance requirement for a hydrogen bond, but does not meet the angle requirement. An inferred value of non-hydrogen bound OH donors is not exact because there is not a rigorously defined maximum value for the number of hydrogen bonds per water, (e.g. a single OH donor may qualify as part of a hydrogen bond at two adjacent water's O-acceptor sites, making the maximum possible number of hydrogen bonds for the water of this particular OH donor greater than 4). Is there a way to implement g_hbond to do this? If not, any specific information on what to consider, and how to modify the g_hbond source code to count non-hydrogen bound water OH donors is greatly appreciated. Robin -- Robin C. Underwood Chemistry Department 560 Oval Drive West Lafayette, IN 47907 -- gmx-users mailing listgmx-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] Re: g_hbond modification
Thanks for the prompt reply Justin, but I'm still not convinced that I can extract the information I seek from the g_hbond program without modification. I realize this wasn't clear in my initial message: I would like to know how many OH donors are not participating in ANY hydrogen bonds. If I calculate the number of non-hydrogen bound OH donors as the difference between the number of pairs that meet the distance requirement, and the actual number of hydrogen bonds (column 3 minus column 2 in hbnum), then I would be over-counting, as this does not exclude pairs with OH donors that are also included in a hydrogen bond with another O acceptor. Robin Quoting Robin C. Underwood rcund...@purdue.edu: I would like to modify the g_hbond code (or preferably, know how to implement g_hbond if it is already capable) to count the number of non-hydrogen bound water OH donors in a simulation. I define a non-hydrogen bound donor as one that does not meet the distance requirement, or one that may meet the distance requirement for a hydrogen bond, but does not meet the angle requirement. An inferred value of non-hydrogen bound OH donors is not exact because there is not a rigorously defined maximum value for the number of hydrogen bonds per water, (e.g. a single OH donor may qualify as part of a hydrogen bond at two adjacent water's O-acceptor sites, making the maximum possible number of hydrogen bonds for the water of this particular OH donor greater than 4). Is there a way to implement g_hbond to do this? If not, any specific information on what to consider, and how to modify the g_hbond source code to count non-hydrogen bound water OH donors is greatly appreciated. Robin -- Robin C. Underwood Chemistry Department 560 Oval Drive West Lafayette, IN 47907 -- Robin C. Underwood Chemistry Department 560 Oval Drive West Lafayette, IN 47907 -- gmx-users mailing listgmx-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] how to choose rlist and cutoffs
Greetings, My system consists of large non-polar spheres in tip4p water, and the user potential option is implemented to define vdw interactions. Because the radius of the sphere is on the order of 1 nm, I am using a cutoff of 1.9 nm for the vdw portion. I would like to save on computational expense by using a cutoff for electrostatics (non-user defined) and choosing a much shorter electrostatics cutoff, around 0.9 nm. Because I must also decrease my rlist value to 0.9, I am wondering how this may effect the interaction energy between the sphere and water. In particular, if my rlist is not much larger than the radius of the sphere, or less than that radius, will I be neglecting important energetic information with rlist=0.9 (even though the electrostatic interaction between the sphere and water is practically 0)? I did see section 4.6.3 in the Gromacs manual, which leads me to think that it is reasonable to decrease rlist because my vdw interactions should still be computed, but I want to be sure this is correct. Robin -- Robin C. Underwood Chemistry Department 560 Oval Drive West Lafayette, IN 47907 -- gmx-users mailing listgmx-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] vapor phase pmf correction
Hello Gromacs users, I have a simple question about whether it is possible to implement a kinetic energy correction to vapor phase potential of mean force calculations using g_wham. I have used the umbrella pulling method (outlined in the umbrella sampling tutorial by Justin on the Gromacs tutorial web page) but for a very simple system- pulling apart two methanes. When I run g_wham, my resulting pmf graph is shifted uniformly by about -13 kJ/mol from achieving a baseline of zero, and has a deeper minimum than I would expect based on the literature. Thus, I would like to use g_wham on only the potential energy portion, excluding the kinetic portion from my calculation. Is it possible to do this? I have read the entropic effects section 6 in the manual, but do not understand the statistical mechanical origin of this equation 6.1. I assume r represents the pull parameter, in my case, distance between the COM of the two groups. If so, such a correction does not shift my pmf to a baseline of 0. Robin -- Robin C. Underwood Chemistry Department 560 Oval Drive West Lafayette, IN 47907 -- gmx-users mailing listgmx-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
