[gmx-users] Virtual sites parameters
Dear GROMACS user, I have been trying to define parameters for virtual sites in CHARMM lipid, and I would like to know if there is something simillar to the [ bondtypes ] for virtual sites. Specifically the way I specify the parameters now is by having lines in the topology such as (for a 3fout virtual sites): [ virtual_sites3 ] ; aiajakal functc0c1c2 40393751 4 -0.26909 -0.26909 4.5519 But then I have to do that for every virtual sites in the molecule, even if the parameters are the same. What I would like to be able to do is: [ ] ; aiajakal functc0c1c2 H C CC 4 -0.26909 -0.26909 4.5519 and then [ virtual_sites3 ] ; aiajakal functc0c1c2 40393751 4 Where atom 40 is of type H and atom 39, 37 and 51 of type C and should be are directive similar to [ bondtypes ] but for virtual site. Does this directive exist ? I had no luck checking the manual. Best regards, Bastien Loubet -- View this message in context: http://gromacs.5086.x6.nabble.com/Virtual-sites-parameters-tp5009258.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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: pullx file content with pull_geometry = position
I update the post to say that we tried to make a rerun of the trajectory on our local machine and the output of the pullx.xvg file (created with the option -px of mdrun) is: # This file was created Wed Jun 5 12:15:49 2013 # by the following command: # mdrun -s ../conf10.tpr -rerun ../conf10.xtc -px pullx_conf10.xvg -pf pullf_conf10.xvg -nt 6 -v # # mdrun is part of G R O M A C S: # # Gravel Rubs Often Many Awfully Cauterized Sores # @title Pull COM @xaxis label Time (ps) @yaxis label Position (nm) @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 X @ s1 legend 0 Y @ s2 legend 0 Z @ s3 legend 1 dX @ s4 legend 1 dY @ s5 legend 1 dZ 0. 3.47835 5.10501 11.644 -0.378345 1.13499 1.73604 20. 3.47852 5.12016 11.6554 -0.446518 1.07184 1.69255 40. 3.49803 5.12013 11.6612 -0.350032 1.12387 1.75475 (...) ** The output of the pullx.xvg file is indeed the distance between the two pull groups, but this is very different numbers than the one obtained during the simulation on a cluster (see above post). So the question stand what are the numbers given in the pullx.xvg file created during the original simulation ? -- View this message in context: http://gromacs.5086.x6.nabble.com/pullx-file-content-with-pull-geometry-position-tp5008797p5008850.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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] pullx file content with pull_geometry = position
Dear gromacs user, I run a simulation where I restrain two groups in a specific position with respect to one another: one is part of a protein (index group name: pull_group) the other is just one ion (index pull_group name r_60022). I have the following parameters in my mdp file: * ; Pull code pull= umbrella pull_geometry = position ; pull_dim= Y Y Y pull_start = no pull_ngroups= 1 pull_group0 = pull_group pull_group1 = r_60022 pull_init1 = -.4110 1.0700 1.7760 pull_rate1 = 0.0 pull_vec1 = -0.264 0.315 0.912 pull_k1 = 1000 ; kJ mol^-1 nm^-2 pull_nstxout= 100 ; every 0.2 ps pull_nstfout= 100 ; every 0.2 ps ** So the ion should be restrained around the position pull_init1 with respect to the center of mass of the protein. Here is a sample of the obtained pullx.xvg file: ** @title Pull COM @xaxis label Time (ps) @yaxis label Position (nm) @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 X @ s1 legend 0 Y @ s2 legend 0 Z @ s3 legend 1 dX @ s4 legend 1 dY @ s5 legend 1 dZ 0. 3.43102 5.2156 11.5301 -0.531018 1.0244 2.13988 0.2000 3.43096 5.21543 11.53 -0.544007 1.02157 2.14952 0.4000 3.43094 5.2146 11.5304 -0.551138 1.00177 2.12238 0.6000 3.43074 5.21354 11.5312 -0.559285 1.01417 2.13649 (...) 20. 3.45331 5.22215 11.51 -0.646908 0.9390662.11768 (...) 40. 3.45141 5.21115 11.5685 -0.689742 0.9789082.14589 (...) * As far as I understood the manual and after searching the forum I though that the first column is the time, the next three columns the position of the pulled group (r_60022) at the time and the last three columns the distance in the x y z direction between the reference group (pull_group, bad name choice I know) and the pulled group (r_60022). If that was true it would suggest that the ion is pulled more than three A away from the center of the umbrella potential. After checking the trajectory visually and by calculating the distance between pull_group and r_60022 using g_dist on the trajectory we can see it is not the case. Here is a sample of the g_dist output: @title Distance @xaxis label Time (ps) @yaxis label Distance (nm) @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 |d| @ s1 legend d\sx\N @ s2 legend d\sy\N @ s3 legend d\sz\N 0.0002.1085827 -0.37829451.13508751.7362576 20.0002.0527425 -0.44645331.07192131.6927538 40.0002.1132109 -0.34997841.12395621.7549639 (...) *** Which oscillate correctly around the pull_init1 vector as defined in the mdp file. More importantly these results are different from the results in the pullx.xvg file (at least the results I expected to have). Can somebody clarify what is actually in the pullx.xvg file ? Best, Bastien Loubet -- View this message in context: http://gromacs.5086.x6.nabble.com/pullx-file-content-with-pull-geometry-position-tp5008797.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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] Lipids and virtual site in 4.6.1
Dear GROMACS user, Last week I was watching the GTC talk from Eric Lindahl, and I noticed his insistence on virtual site to accelerate simulation (up to 5ps time step). As a matter of fact our group recently tried to use virtual site on CHARMM36 POPC bilayer and the result where less than stellar. We first did a simulation of a POPC bilayer comparing area per lipid with and without virtual site using GROMACS 4.5.5. We obtained an area per lipid around 0.62 nm^2 per lipid with no virtual site but with virtual site we obtained an area per lipid of more than 0.7 nm^2. After that we found that there was a bug report for the virial calculation with virtual site in version 4.5: http://redmine.gromacs.org/issues/908 The issue being marked as resolved for 4.6 and we tried our hand at both the beta3 and 4.6.1 version of it. The result where better but not perfect: we get 0.65~0.66 nm^2 per lipid. The problem here is that in order to increase the time step we need all hydrogen to be virtual sites, so this include the lipids hydrogen and not just the protein. However if the area per lipid is wrong the stress profile is likely to be changed, which will affect transmembrane proteins. So we do not feel confident about using virtual sites in our lipid simulation. I signal it here because there is a possibility that this is still a bug, If somebody could confirm/infirm that it would be nice. Best, Bastien -- View this message in context: http://gromacs.5086.n6.nabble.com/Lipids-and-virtual-site-in-4-6-1-tp5007063.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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: Lipids and virtual site in 4.6.1
Actually we did not increase the time step between the simulations. The plan was to increase it after we checked the basic properties of the membrane, but we never got to that point. I will try to put something together for redmine soon. Thanks for the answer, Bastien -- View this message in context: http://gromacs.5086.n6.nabble.com/Lipids-and-virtual-site-in-4-6-1-tp5007063p5007077.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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] virtual sites with CHARMM lipids
Dear GROMACS users, I have been working on a CHARMM36 POPC bilayer and in order to accelerate the simulation we want to use virtual sites for the hydrogen atoms in the lipid. We mainly want to do it for the carbon tails of the lipids. I have checked the different type of virtual sites in the GROMACS manual (3fd, 3fda, 3out, N) but none of them seems to correspond to the CH2 groups in the carbon chains. It would require an out of plane virtual site with fixed length (so using normalized vector). Can you help me with that ? Thanks, Bastien Loubet -- View this message in context: http://gromacs.5086.n6.nabble.com/virtual-sites-with-CHARMM-lipids-tp5003962.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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: Possible bug in the temperature calculation from rerun
Dear gmxusers,
Thank you for your answers so far. I have run a few more test rerun and here
is what I got:
-The number of degree of freedom are the same whatever the topology I run
with, both in the log file and from the grompp output.
-Comparing the trajectories created from rerun with gmxdump and diff show
only minor differences, such as atoms positioned on one side or the other of
the simulation box. The velocities are strictly the same however. Here is an
example:
8836288,8836290c8836288,8836290
x[394728]={ 1.32183e+02, 3.15341e+00, 5.17869e+00}
x[394729]={ 1.32054e+02, 3.10839e+00, 5.14599e+00}
x[394730]={ 1.32314e+02, 3.17612e+00, 5.13366e+00}
---
x[394728]={-7.17163e-04, 3.15341e+00, 5.17869e+00}
x[394729]={-1.29196e-01, 3.10839e+00, 5.14599e+00}
x[394730]={ 1.29883e-01, 3.17612e+00, 5.13366e+00}
-If I use another topology where I only have the interaction with water
not set to zero (still with electrostatic), I get a temperature of ~326K
intermediate between the desired value of 325K and the value I got with all
the interaction shutdown of ~327K.
-Turning the temperature/pressure coupling on or off do not change the
results.
I think there is a problem in the calculation of the kinetic energy that
translate to the temperature. However setting all the interaction but the
electrostatic to zero should not modify the kinetic energy calculated from
'mdrun -rerun'.
I do not know how to proceed from here, maybe I should post my mdp top and
log files here or on the developers forum ?
Have a good day,
Bastien Loubet
--
View this message in context:
http://gromacs.5086.n6.nabble.com/Possible-bug-in-the-temperature-calculation-from-rerun-tp5001194p5001508.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
--
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] Possible bug in the temperature calculation from rerun
Dear gmx users, We recently got a problem with the rerun feature of mdrun, and we request your help in order to help to solve it. We have run a simulation of a large POPC membrane using the coarse grained Martini force fields. From these simulations we obtained a trr trajectory file which contains both the position and the velocity of each martini particle every 20 ps. From this trajectory we rerun the simulation using the -rerun option of the mdrun program but with a different topology for which we have set all the bonded and non bonded interaction to zero. Specifically we have set all the force constant to zero in the martini_v2.P.itp and the martini_v2.0_lipids.itp file, the aim is to calculate the virial due to the electrostatic forces alone. The problem is the following: from the edr file we get from the rerun we extract the mean value of the energies using g_energy. For the edr file obtained during the simulation (and not the rerun) we obtained for the temperature and the kinetic energy, using also a 20ps time interval: Energy Average Err.Est. RMSD Tot-Drift --- Temperature 324.998 0.0049 0.413669 -0.00743801 (K) Kinetic En. 1.9525e+06 292485.21 -44.6815 (kJ/mol) This is to be expected because we have a Nose-Hover temperature coupling of 325 K. The kinetic energy is also equals to the number of degree of freedom times half the Boltzmann constant times the temperature. However for the rerun, with no bonded and non-bonded interactions, we get: Energy Average Err.Est. RMSD Tot-Drift --- Temperature 327.35 0.0053 0.415978 -0.00514402 (K) Kinetic En. 1.96663e+06 322499.08 -30.8973 (kJ/mol) The kinetic energy is again equals to the number of degree of freedom times half the Boltzmann constant times the temperature, however the temperature is 2 K off ! This is surprising as only the velocities and the masses of the particle should enter the kinetic energy. The velocities are taken from the trr trajectory file and we have checked using the Linux diff utility that our topology have the same masses than the original one. Another cause can be that we are doing the rerun on a local machine while the original simulation was run on a cluster, however we believe that this kind of numerical error cannot be responsible for a 2 degree Kelvin mistake. Furthermore we have rerun the simulation using the original topology and we obtained: Energy Average Err.Est. RMSD Tot-Drift --- Temperature 324.998 0.005 0.411827 -0.00761929 (K) Kinetic En. 1.9525e+06 302474.15 -45.7788 (kJ/mol) The slightly different value are certainly due to the previously mentioned fact that the reruns were not run on the same machine as the simulation and numerical rounding errors. This point out that it is really the topology that changes the value of the temperature. However we really have no masses differences between the two topologies and the velocities are draw from the same trr trajectory file. Any thought on this problem would be welcome. Cheers, Bastien Loubet -- View this message in context: http://gromacs.5086.n6.nabble.com/Possible-bug-in-the-temperature-calculation-from-rerun-tp5001194.html Sent from the GROMACS Users Forum mailing list archive at Nabble.com. -- 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: Possible bug in the temperature calculation from rerun
Justin Lemkul wrote On 9/21/12 8:29 AM, Bastien Loubet wrote: Dear gmx users, We recently got a problem with the rerun feature of mdrun, and we request your help in order to help to solve it. We have run a simulation of a large POPC membrane using the coarse grained Martini force fields. From these simulations we obtained a trr trajectory file which contains both the position and the velocity of each martini particle every 20 ps. From this trajectory we rerun the simulation using the -rerun option of the mdrun program but with a different topology for which we have set all the bonded and non bonded interaction to zero. Specifically we have set all the force constant to zero in the martini_v2.P.itp and the martini_v2.0_lipids.itp file, the aim is to calculate the virial due to the electrostatic forces alone. The problem is the following: from the edr file we get from the rerun we extract the mean value of the energies using g_energy. For the edr file obtained during the simulation (and not the rerun) we obtained for the temperature and the kinetic energy, using also a 20ps time interval: Energy Average Err.Est. RMSD Tot-Drift --- Temperature 324.998 0.0049 0.413669 -0.00743801 (K) Kinetic En. 1.9525e+06 292485.21 -44.6815 (kJ/mol) This is to be expected because we have a Nose-Hover temperature coupling of 325 K. The kinetic energy is also equals to the number of degree of freedom times half the Boltzmann constant times the temperature. However for the rerun, with no bonded and non-bonded interactions, we get: Energy Average Err.Est. RMSD Tot-Drift --- Temperature 327.35 0.0053 0.415978 -0.00514402 (K) Kinetic En. 1.96663e+06 322499.08 -30.8973 (kJ/mol) The kinetic energy is again equals to the number of degree of freedom times half the Boltzmann constant times the temperature, however the temperature is 2 K off ! This is surprising as only the velocities and the masses of the particle should enter the kinetic energy. The velocities are taken from the trr trajectory file and we have checked using the Linux diff utility that our topology have the same masses than the original one. Another cause can be that we are doing the rerun on a local machine while the original simulation was run on a cluster, however we believe that this kind of numerical error cannot be responsible for a 2 degree Kelvin mistake. Furthermore we have rerun the simulation using the original topology and we obtained: Energy Average Err.Est. RMSD Tot-Drift --- Temperature 324.998 0.005 0.411827 -0.00761929 (K) Kinetic En. 1.9525e+06 302474.15 -45.7788 (kJ/mol) The slightly different value are certainly due to the previously mentioned fact that the reruns were not run on the same machine as the simulation and numerical rounding errors. This point out that it is really the topology that changes the value of the temperature. However we really have no masses differences between the two topologies and the velocities are draw from the same trr trajectory file. Any thought on this problem would be welcome. Another thing to consider is the fact that an .edr file is accurate over all simulation steps, while reruns only do calculations present in the frames supplied. If you are using frames every 20 ps, that may only represent less than 1% of the total data collected in the simulation, depending on the value of dt. -Justin -- Justin A. Lemkul, Ph.D. Research Scientist 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@ 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-request@ . * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists Alright, but we did make a rerun with the same trajectory and the original topology and it gives the third results in my original post: the 2 Kelvin difference is not there. So it is really the different topology file that give the temperature difference I believe. -- View this message in context: http://gromacs.5086.n6.nabble.com/Possible-bug-in-the-temperature-calculation-from-rerun
