[gmx-users] cut-off g_hbond
Hi all I am using g_hbond to calculate the number of hydrogen bonds in an alcohol system. I am using the following command: g_hbond -f traj.trr -nonitacc I am slightly confused about the cut off r. Using the above command, is the default cut-off r (0.35nm) for the donor-acceptor distance or the hydrogen-acceptor ? Cheers Gavin -- 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
Re: [gmx-users] cut-off g_hbond
Hi Justin Thanks, that's what I thought. Is there a need to define a hydrogen-acceptor distance as well. I read in a few articles that this was the case, usually 0.25nm ? Cheers Gavin Justin Lemkul wrote: On 10/5/12 7:20 AM, Gavin Melaugh wrote: Hi all I am using g_hbond to calculate the number of hydrogen bonds in an alcohol system. I am using the following command: g_hbond -f traj.trr -nonitacc I am slightly confused about the cut off r. Using the above command, is the default cut-off r (0.35nm) for the donor-acceptor distance or the hydrogen-acceptor ? This is controlled by the -da switch. The default is donor-acceptor; if one uses g_hbond -noda, then the cutoff is applied to the hydrogen-acceptor distance. -Justin -- 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
Re: [gmx-users] cut-off g_hbond
Is there any need to use the r2 option the ? Justin Lemkul wrote: On 10/5/12 7:39 AM, Gavin Melaugh wrote: Hi Justin Thanks, that's what I thought. Is there a need to define a hydrogen-acceptor distance as well. I read in a few articles that this was the case, usually 0.25nm ? One can define hydrogen bonds in a number of ways, but a D-A distance of 0.35 nm and an H-A distance of 0.25 nm are basically equivalent, given that O-H and N-H bonds are about 0.1 nm. -Justin -- 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
Re: [gmx-users] cut-off g_hbond
O.K thanks Justin Lemkul wrote: On 10/5/12 8:15 AM, Gavin Melaugh wrote: Is there any need to use the r2 option the ? That option is not relevant here. It is undocumented in the current release, but the help description has been updated for the next release: http://redmine.gromacs.org/projects/gromacs/repository/revisions/1fd5b54ee7e1c3889a9b110690b36374136a2cbf/diff/src/tools/gmx_hbond.c -Justin -- 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] g_hbond
Hi I have an alcohol system and I want to calculate the number of H bonds during the trajectory. My atom type labels are C2-C2-C2-CO-OH-HO. CO denotes carbon bonded to oxygen, OH denotes alcohol oxygen, and HO denotes alcohol hydrogen. In g_hbond how do I specify my groups to consider for H bonding. If I uses group 0 which is the whole system it says that I have twice as many acceptors as donors, which shouldn't be the case. Cheers Gavin -- 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] wham histograms
Hi all Are the histograms from histo.xvg (output of g_wham) the biased or unbiased distributions? Cheers Gavin -- 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] extra factor in PMF
Hi all According to some references 1) J. Chem. Phys, 128, 044106, (2008) 2) Comparison of Methods to Compute the Potential of Mean Force Daniel Trzesniak, Anna-Pitschna E. Kunz, Wilfred F. van Gunsteren Prof. Article first published online: 28 NOV 2006 DOI: 10.1002/cphc.200600527 There is an extra factor that has to be taken into account when calculating the PMF. This extra factor (2kTln(r)) results from the transformation from the 3N Cartesian coordinates to 3N internal coordinates. In the case of a reaction coordinate defined by a distance r between two species the Jacobian (r*2sin(theta)) of the transformation for 3D to 1D leads to an extra term in the PMF. Therefore the true PMF should be PMF(true) = PMF(g_wham) + 2kTln(r) This extra factor is only significant at lower values of the reaction coordinate. My question is; Does g_wham take this into account ? Cheers Gavin -- 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] DD load balancing is limited by minimum cell size in dimension Z
Hi all What does the following note mean in the log file DD load balancing is limited by minimum cell size in dimension Z Is is purely a performance related issue ? Cheers Gavin -- 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] genbox
Hi all I have a system of 40 solute molecules in a solvent of 480 crown ether molecules. I am not trying to insert 100 methane molecules into this relaxed and well equilibrated structure using genbox. There are clearly visible cavities in the fluid but genbox only alllows the insertion of 8 methane molecules. How can I circumvent this problem ? The command I use is genbox -cp test.gro -ci methane.gro -nmol 100 -try 5 -p combined.top where combined.top includes all three itp files, and test.gro is the initial solute and solvent configuration. Cheers Gavin -- 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] fatal error
Hi all I am trying to relax the density of an alcohol system that I have just set up, but after a certain period of time I get the following error Fatal error: The X-size of the box (6.800588) times the triclinic skew factor (1.00) is smaller than the number of DD cells (4) times the smallest allowed cell size (1.70) I have never had this error before when I ran simulations on much more processors for the alkane analogues of the alcohol. What is the source of this error ? Cheers Gavin -- 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] bug in g_msd
Hi all Is there a bug with g_msd when using the -mol flag. I have my own program that calculates the MSD. If I compare it with the gromacs utility for one system the curves are the exact same, however when I compare with another system th curves are very different. Someone else mentioned something similar a few days ago, so I was just wondering if this was the case. Cheers Gavin -- 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] [Fwd: bug in g_msd]
---BeginMessage--- Hi all Is there a bug with g_msd when using the -mol flag. I have my own program that calculates the MSD. If I compare it with the gromacs utility for one system the curves are the exact same, however when I compare with another system th curves are very different. Someone else mentioned something similar a few days ago, so I was just wondering if this was the case. Cheers Gavin ---End Message--- -- 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] g_msd segmentation fault
Dear All I have a system of 40 solute molecules in 480 crown ether solvent molecules. When I ran the msd analysis on the solvent molecules using the following comand. g_msd_login_d -f traj.trr -s topol.tpr -mol -o 1_12_400K_sol_msd.xvg I get a segmentation fault as follows Select a group: 4 Selected 4: 'SOL' Split group of 16800 atoms into 480 molecules trn version: GMX_trn_file (double precision) Reading frame 800 time 4.000 Segmentation fault However when I perform the same analysis on the solute molecules it runs to completion. I have checked the configuration of the system at frame 800 and everything seems to be fine. I have also analysed the energy and there seems to be no problem. Has anybody any idea of what might be happening? Cheers Gavin -- 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] folding of coordinates in trajectory files
Hi All I was wondering. I assume gromacs writes folded coordinates in the trajectory files. If so does it use rx(i) = rx(i) -boxl *nint(rx(i) / boxl) to fold the coordinates at each step? and where does it define the origin (0,0,0) Cheers Gavin -- 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
Re: [gmx-users] Error note
Hi Justin I saw this and had a question. How important is it for a charge group to be neutral? In one of my systems I model a solvent of crown ethers using an all atom model. The smallest neutral unit comprises a CH2-O-CH2 (7 atoms). I have used this as a charge group to and get no warnings. Cheers Gavin Justin A. Lemkul wrote: RAMYA NAGA wrote: Dear friends, iam getting this note while doing pressure coupling of protein-ligand complex. can anybody help me how to handle this?? The largest charge group contains 11 atoms. Since atoms only see each other when the centers of geometry of the charge groups they belong to are within the cut-off distance, too large charge groups can lead to serious cut-off artifacts. For efficiency and accuracy, charge group should consist of a few atoms. For all-atom force fields use: CH3, CH2, CH, NH2, NH, OH, CO2, CO, etc. The error message explains what is wrong and how it should be fixed. Charge groups should be small. You have one that is very large. -Justin -- 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
Re: [gmx-users] Error note
I am using 4.5.5 Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin I saw this and had a question. How important is it for a charge group to be neutral? It isn't. Conventional use dictates that the charge group bear an integral charge. With PME, this is not necessary. Many force fields do not use charge groups at all (i.e., single-atom charge groups). In one of my systems I model a solvent of crown ethers using an all atom model. The smallest neutral unit comprises a CH2-O-CH2 (7 atoms). I have used this as a charge group to and get no warnings. The only implication I can think of would be in neighbor searching. If the group is large, then short-range forces may not be calculated accurately. The warning from grompp was a recent addition to the code; if you're using an older version you may not have triggered it. -Justin -- 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] pull-code
Hi all I am returning to a query I had a few weeks ago regarding a discrepancy between two free energy curves. One calculated using umbrella sampling, the other calculated via the reversible work theorem from the RDF. There is sufficient sampling of the dynamics in the RDF so this method is viable. Anyway in the pull-code I use pull_geometry = dist and pull_dim=Y Y Y. The free energy curve from the pull-code method does not give me a minimum at the zero value of the order parameter whereas the RDF method does. Someone said before about double counting of positive distances at small values of the order parameter and therefore information is lost at very small distances. Is this correct? I am slightly concerned that my curves are not giving me the correct information involving a very important state in my reaction coordinate. Also when this dist restraint (which cannot be negative) is implemented are there issues with the normalisation of the histograms from g_wham? Cheers Gavin -- 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
Re: [gmx-users] pull-code
Hi Thomas Many thanks for the reply again. At larger distances the two curves match up quite well. The curve from the reversible work theorem is better behaved and smoother but this could be solely due to statistics. I am slightly confused about your statement If the small circle moves between 0 and any value 0 everything should be fine. Do you not mean 0 and any value 0 ? Cheers Gavin Thomas Schlesier wrote: Hi Gavin, if i remember correctly it was a system about pulling a ligand from a binding pocket? To make the system simpler we have a big circle and in the middle a small circle. And we assume that the potential minimum for the interaction between both circles is when the small cirlce is in the middle of the large circle. Now we do the Umbrella sampling. For a window which is centered at a distance which is sligthly greater then 0, we will get problems. Assume small circle is sligthly shifted to the right. And the other windows are also in this dircetion. (- reaction coordinate goes from zero to the right dircetion) If the small circle moves between 0 and any value 0 everythig should be fine. But if the small circle moves to the left, we will also get a positive distance. Problem is from the above defined reaction coordinate it should be a negative distance. So we are counting the positive distances too much. To check this, you could use *g_dist* to calculate the distance for both molecules for the problematic windows. Then project the resulting vector onto your reaction coordinate. Then you should see the crossings between the right and left side. How do the two free energy curves compare for larger distances, where you can be sure, that you do not have this 'crossing problem'? Greetings Thomas - Hi all I am returning to a query I had a few weeks ago regarding a discrepancy between two free energy curves. One calculated using umbrella sampling, the other calculated via the reversible work theorem from the RDF. There is sufficient sampling of the dynamics in the RDF so this method is viable. Anyway in the pull-code I use pull_geometry = dist and pull_dim=Y Y Y. The free energy curve from the pull-code method does not give me a minimum at the zero value of the order parameter whereas the RDF method does. Someone said before about double counting of positive distances at small values of the order parameter and therefore information is lost at very small distances. Is this correct? I am slightly concerned that my curves are not giving me the correct information involving a very important state in my reaction coordinate. Also when this dist restraint (which cannot be negative) is implemented are there issues with the normalisation of the histograms from g_wham? Cheers Gavin -- 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
Re: [gmx-users] pull-code
Hi Thomas I am sorry to bother you but if you could answer a few questions I have about how the pull-code works with respect to my system I would really appreciate it. My system is a liquid and I am trying to pull one substituent of one liquid molecule in a certain region of another liquid molecule using as I said earlier pull_geometry =dist and pull_dim =YYY. Also pull_start = no Say I have a window with pull_init=0 1) At the very start of the simulation the pullcode calculates the vector between the two groups? Is there anything particularly significant about this initial vector? Is this distance vector recalculated at every step? 2) Does it compare the distance vector an time t with that at time 0.? 3) Given the initial vector can this vector change (i.e. dierection) or does the distance between the two groups vary only along this vector (i.e. in a line) ? Cheers Gavin Thomas Schlesier wrote: Hi Gavin, if i remember correctly it was a system about pulling a ligand from a binding pocket? To make the system simpler we have a big circle and in the middle a small circle. And we assume that the potential minimum for the interaction between both circles is when the small cirlce is in the middle of the large circle. Now we do the Umbrella sampling. For a window which is centered at a distance which is sligthly greater then 0, we will get problems. Assume small circle is sligthly shifted to the right. And the other windows are also in this dircetion. (- reaction coordinate goes from zero to the right dircetion) If the small circle moves between 0 and any value 0 everythig should be fine. But if the small circle moves to the left, we will also get a positive distance. Problem is from the above defined reaction coordinate it should be a negative distance. So we are counting the positive distances too much. To check this, you could use *g_dist* to calculate the distance for both molecules for the problematic windows. Then project the resulting vector onto your reaction coordinate. Then you should see the crossings between the right and left side. How do the two free energy curves compare for larger distances, where you can be sure, that you do not have this 'crossing problem'? Greetings Thomas - Hi all I am returning to a query I had a few weeks ago regarding a discrepancy between two free energy curves. One calculated using umbrella sampling, the other calculated via the reversible work theorem from the RDF. There is sufficient sampling of the dynamics in the RDF so this method is viable. Anyway in the pull-code I use pull_geometry = dist and pull_dim=Y Y Y. The free energy curve from the pull-code method does not give me a minimum at the zero value of the order parameter whereas the RDF method does. Someone said before about double counting of positive distances at small values of the order parameter and therefore information is lost at very small distances. Is this correct? I am slightly concerned that my curves are not giving me the correct information involving a very important state in my reaction coordinate. Also when this dist restraint (which cannot be negative) is implemented are there issues with the normalisation of the histograms from g_wham? Cheers Gavin -- 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
Re: [gmx-users] RDF(PMF) and Umbrella sampling
Hi Justin Again, many thanks for the reply. So when the COM distance changes sign, what effect does that have on the distribution of the COM distance about the mean value for that window i.e. If say my ref dist in 0 nm and the umbrella sampling allows the distance to sample distances say at 0.02 nm to -0.02nm. What happens to negative values? Obviously they are not counted as negative in the distribution or else it would be centred at zero/ Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Thanks very much. One last question. What do you mean when you say COM reference distance is changing signs? I thought the COM distance was the absolute distance between the two groups and therefore cannot be negative? The pull code deals in vectors. Signs can change. The use of distance as a geometry is perhaps somewhat misleading. -Justin Cheers Gavin Dariush Mohammadyani wrote: Hi Gavin, A question arose for me: why did you consider the (rate = 0)? Dariush On Fri, Jan 6, 2012 at 11:47 AM, Gavin Melaugh gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Hi Justin Just a quick clarification regarding my previous point. With geometry = distance, and pull_dim =Y Y Y . Is the pull_group sampling all dimensions equally (or without prejudice) about pull_init ? And iN your first reply what did you mean about by straight pull ? Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Thanks for the reply. I wanted my pulling to be free in all directions, that is in the liquid state with no defined reaction coordinate i.e not along a specific axis. This is why I used geometry = distance. Would you agree with this approach? I suppose there is an argument that can be made for a more free approach such as this one, but you're going to get the artifact you observed the instant your pull group moves past a zero COM distance. Whether or not this is a significant problem is something you'll have to determine. -Justin By free I mean. The absolute distance between the COG of the ref group and that of the pull group. Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Dear all I have a query regarding umbrella sampling simulations that I have carried out to study a dynamical process of a guest inserting into a host. I always get get a wall tending off to infinity at or just before the zero distance between the two species. The process I describe, for one system in particular, happens readily and I have compared the PMF from a non constrained simulation (via the RDF and reversible work theorem) and the same PMF from a set of umbrella sampling simulations. They agree quite well but in the non constrained simulation I get a minimum practically at zero whereas for the umbrella sampling the minimum is shifted and there is an infinite wall close to zero. This wall is not present from the reversible work theorem. Why the infinite wall? Why does the black histogram not centre around zero. Is this an artefact of the umbrella technique? Please see attached the profile from the umbrella sampling technique, and the corresponding histograms. What's happening is the COM reference distance is changing signs, so you get an artifact. The distance geometry is relatively inflexible and is only suitable for straight pulls of continuously increasing or continuously decreasing COM distance. You should try using the position geometry instead. There are some notes that you may find useful in my tutorial: http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/05a_pull_tips.html -Justin Here is an excerpt from one of the umbrella mdp files. pull= umbrella pull_geometry = distance pull_dim = Y Y Y pull_start = no pull_ngroups = 1 pull_group0 = cage_1 pull_group1 = tail pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1 pull_nstxout = 150 pull_nstfout = 150 Cheers Gavin -- gmx-users mailing listgmx-users@gromacs.org mailto:gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support
Re: [gmx-users] RDF(PMF) and Umbrella sampling
Hi Thomas Many Thanks for the reply. I have attached the histograms I get for the umbrella sampling. My concern is the black histogram in which I have set the ref dist to be zero. The distances in the first few windows are very small as I really wanted this difficult region to be sampled well and therefore there might be some sampling in the negative distances. Does the picture here fit with your explanation below ? Cheers Gavin Thomas Schlesier wrote: I think that for the histogram all contribution with a negative sign would be add to the contributions from positive distances. If your distribution is be a perfect gaussian with zero mean, you would end up with half a gaussian with double high (for positive distances) and zero for negative distances. For this case it would be easy to correct the histogram. But if the histogram isn't a perfect gaussian, you couldn't say anything. Date: Tue, 10 Jan 2012 10:47:43 + From: Gavin Melaughgmelaug...@qub.ac.uk Subject: Re: [gmx-users] RDF(PMF) and Umbrella sampling To: Discussion list for GROMACS usersgmx-users@gromacs.org Message-ID:4f0c174f.2050...@qub.ac.uk Content-Type: text/plain; charset=ISO-8859-1 Hi Justin Again, many thanks for the reply. So when the COM distance changes sign, what effect does that have on the distribution of the COM distance about the mean value for that window i.e. If say my ref dist in 0 nm and the umbrella sampling allows the distance to sample distances say at 0.02 nm to -0.02nm. What happens to negative values? Obviously they are not counted as negative in the distribution or else it would be centred at zero/ Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Thanks very much. One last question. What do you mean when you say COM reference distance is changing signs? I thought the COM distance was the absolute distance between the two groups and therefore cannot be negative? The pull code deals in vectors. Signs can change. The use of distance as a geometry is perhaps somewhat misleading. -Justin Cheers Gavin Dariush Mohammadyani wrote: Hi Gavin, A question arose for me: why did you consider the (rate = 0)? Dariush On Fri, Jan 6, 2012 at 11:47 AM, Gavin Melaughgmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Hi Justin Just a quick clarification regarding my previous point. With geometry = distance, and pull_dim =Y Y Y . Is the pull_group sampling all dimensions equally (or without prejudice) about pull_init ? And iN your first reply what did you mean about by straight pull ? Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Thanks for the reply. I wanted my pulling to be free in all directions, that is in the liquid state with no defined reaction coordinate i.e not along a specific axis. This is why I used geometry = distance. Would you agree with this approach? I suppose there is an argument that can be made for a more free approach such as this one, but you're going to get the artifact you observed the instant your pull group moves past a zero COM distance. Whether or not this is a significant problem is something you'll have to determine. -Justin By free I mean. The absolute distance between the COG of the ref group and that of the pull group. Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Dear all I have a query regarding umbrella sampling simulations that I have carried out to study a dynamical process of a guest inserting into a host. I always get get a wall tending off to infinity at or just before the zero distance between the two species. The process I describe, for one system in particular, happens readily and I have compared the PMF from a non constrained simulation (via the RDF and reversible work theorem) and the same PMF from a set of umbrella sampling simulations. They agree quite well but in the non constrained simulation I get a minimum practically at zero whereas for the umbrella sampling the minimum is shifted and there is an infinite wall close to zero. This wall is not present from the reversible work theorem. Why the infinite wall? Why does the black histogram not centre around zero. Is this an artefact of the umbrella technique? Please see attached the profile from the umbrella sampling technique, and the corresponding histograms. What's happening is the COM reference distance is changing signs, so you get
[gmx-users] RDF(PMF) and Umbrella sampling
Dear all I have a query regarding umbrella sampling simulations that I have carried out to study a dynamical process of a guest inserting into a host. I always get get a wall tending off to infinity at or just before the zero distance between the two species. The process I describe, for one system in particular, happens readily and I have compared the PMF from a non constrained simulation (via the RDF and reversible work theorem) and the same PMF from a set of umbrella sampling simulations. They agree quite well but in the non constrained simulation I get a minimum practically at zero whereas for the umbrella sampling the minimum is shifted and there is an infinite wall close to zero. This wall is not present from the reversible work theorem. Why the infinite wall? Why does the black histogram not centre around zero. Is this an artefact of the umbrella technique? Please see attached the profile from the umbrella sampling technique, and the corresponding histograms. Here is an excerpt from one of the umbrella mdp files. pull= umbrella pull_geometry = distance pull_dim = Y Y Y pull_start = no pull_ngroups = 1 pull_group0 = cage_1 pull_group1 = tail pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1 pull_nstxout = 150 pull_nstfout = 150 Cheers Gavin inline: histo.pnginline: profile.png-- 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
Re: [gmx-users] RDF(PMF) and Umbrella sampling
Hi Justin Thanks for the reply. I wanted my pulling to be free in all directions, that is in the liquid state with no defined reaction coordinate i.e not along a specific axis. This is why I used geometry = distance. Would you agree with this approach? By free I mean. The absolute distance between the COG of the ref group and that of the pull group. Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Dear all I have a query regarding umbrella sampling simulations that I have carried out to study a dynamical process of a guest inserting into a host. I always get get a wall tending off to infinity at or just before the zero distance between the two species. The process I describe, for one system in particular, happens readily and I have compared the PMF from a non constrained simulation (via the RDF and reversible work theorem) and the same PMF from a set of umbrella sampling simulations. They agree quite well but in the non constrained simulation I get a minimum practically at zero whereas for the umbrella sampling the minimum is shifted and there is an infinite wall close to zero. This wall is not present from the reversible work theorem. Why the infinite wall? Why does the black histogram not centre around zero. Is this an artefact of the umbrella technique? Please see attached the profile from the umbrella sampling technique, and the corresponding histograms. What's happening is the COM reference distance is changing signs, so you get an artifact. The distance geometry is relatively inflexible and is only suitable for straight pulls of continuously increasing or continuously decreasing COM distance. You should try using the position geometry instead. There are some notes that you may find useful in my tutorial: http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/05a_pull_tips.html -Justin Here is an excerpt from one of the umbrella mdp files. pull= umbrella pull_geometry = distance pull_dim = Y Y Y pull_start = no pull_ngroups = 1 pull_group0 = cage_1 pull_group1 = tail pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1 pull_nstxout = 150 pull_nstfout = 150 Cheers Gavin -- 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
Re: [gmx-users] RDF(PMF) and Umbrella sampling
Hi Justin Just a quick clarification regarding my previous point. With geometry = distance, and pull_dim =Y Y Y . Is the pull_group sampling all dimensions equally (or without prejudice) about pull_init ? And iN your first reply what did you mean about by straight pull ? Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Thanks for the reply. I wanted my pulling to be free in all directions, that is in the liquid state with no defined reaction coordinate i.e not along a specific axis. This is why I used geometry = distance. Would you agree with this approach? I suppose there is an argument that can be made for a more free approach such as this one, but you're going to get the artifact you observed the instant your pull group moves past a zero COM distance. Whether or not this is a significant problem is something you'll have to determine. -Justin By free I mean. The absolute distance between the COG of the ref group and that of the pull group. Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Dear all I have a query regarding umbrella sampling simulations that I have carried out to study a dynamical process of a guest inserting into a host. I always get get a wall tending off to infinity at or just before the zero distance between the two species. The process I describe, for one system in particular, happens readily and I have compared the PMF from a non constrained simulation (via the RDF and reversible work theorem) and the same PMF from a set of umbrella sampling simulations. They agree quite well but in the non constrained simulation I get a minimum practically at zero whereas for the umbrella sampling the minimum is shifted and there is an infinite wall close to zero. This wall is not present from the reversible work theorem. Why the infinite wall? Why does the black histogram not centre around zero. Is this an artefact of the umbrella technique? Please see attached the profile from the umbrella sampling technique, and the corresponding histograms. What's happening is the COM reference distance is changing signs, so you get an artifact. The distance geometry is relatively inflexible and is only suitable for straight pulls of continuously increasing or continuously decreasing COM distance. You should try using the position geometry instead. There are some notes that you may find useful in my tutorial: http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/05a_pull_tips.html -Justin Here is an excerpt from one of the umbrella mdp files. pull= umbrella pull_geometry = distance pull_dim = Y Y Y pull_start = no pull_ngroups = 1 pull_group0 = cage_1 pull_group1 = tail pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1 pull_nstxout = 150 pull_nstfout = 150 Cheers Gavin -- 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
Re: [gmx-users] RDF(PMF) and Umbrella sampling
Hi Justin Thanks very much. One last question. What do you mean when you say COM reference distance is changing signs? I thought the COM distance was the absolute distance between the two groups and therefore cannot be negative? Cheers Gavin Dariush Mohammadyani wrote: Hi Gavin, A question arose for me: why did you consider the (rate = 0)? Dariush On Fri, Jan 6, 2012 at 11:47 AM, Gavin Melaugh gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Hi Justin Just a quick clarification regarding my previous point. With geometry = distance, and pull_dim =Y Y Y . Is the pull_group sampling all dimensions equally (or without prejudice) about pull_init ? And iN your first reply what did you mean about by straight pull ? Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Thanks for the reply. I wanted my pulling to be free in all directions, that is in the liquid state with no defined reaction coordinate i.e not along a specific axis. This is why I used geometry = distance. Would you agree with this approach? I suppose there is an argument that can be made for a more free approach such as this one, but you're going to get the artifact you observed the instant your pull group moves past a zero COM distance. Whether or not this is a significant problem is something you'll have to determine. -Justin By free I mean. The absolute distance between the COG of the ref group and that of the pull group. Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Dear all I have a query regarding umbrella sampling simulations that I have carried out to study a dynamical process of a guest inserting into a host. I always get get a wall tending off to infinity at or just before the zero distance between the two species. The process I describe, for one system in particular, happens readily and I have compared the PMF from a non constrained simulation (via the RDF and reversible work theorem) and the same PMF from a set of umbrella sampling simulations. They agree quite well but in the non constrained simulation I get a minimum practically at zero whereas for the umbrella sampling the minimum is shifted and there is an infinite wall close to zero. This wall is not present from the reversible work theorem. Why the infinite wall? Why does the black histogram not centre around zero. Is this an artefact of the umbrella technique? Please see attached the profile from the umbrella sampling technique, and the corresponding histograms. What's happening is the COM reference distance is changing signs, so you get an artifact. The distance geometry is relatively inflexible and is only suitable for straight pulls of continuously increasing or continuously decreasing COM distance. You should try using the position geometry instead. There are some notes that you may find useful in my tutorial: http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/05a_pull_tips.html -Justin Here is an excerpt from one of the umbrella mdp files. pull= umbrella pull_geometry = distance pull_dim = Y Y Y pull_start = no pull_ngroups = 1 pull_group0 = cage_1 pull_group1 = tail pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1 pull_nstxout = 150 pull_nstfout = 150 Cheers Gavin -- gmx-users mailing listgmx-users@gromacs.org mailto: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 mailto:gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- Kind Regards, Dariush Mohammadyani Department of Structural Biology University of Pittsburgh School of Medicine Biomedical Science Tower 3 3501 Fifth Avenue Pittsburgh, PA 15261 USA -- 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
[gmx-users] [Fwd: h-bonds constraints]
---BeginMessage--- Hi I want to run an NPT simulation with all h-bonds constrained. How does grompp identify the Hydrogen atoms given that forcefield labels like HA, HC, HE are used. Is it the mass? Many Thanks Gavin ---End Message--- -- 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] h-bonds constraints
Hi I want to run an NPT simulation with all h-bonds constrained. How does grompp identify the Hydrogen atoms given that forcefield labels like HA, HC, HE are used. Is it the mass? Many Thanks Gavin -- 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] discrepancy trjconv and gmxdump
Hi all I have generated a gro file from a traj.trr file using trjconv. When I use gmxdump on the same traj.trr file to output a generic format history file it seems that there is a discrepancy in the coordinates of some atoms in a particular frame. Essentially I output frame 5 using trjconv, then I take the corresponding frame from the history file and covert it to gro format. I then run vimdiff on the two files. There are sign differences with the coordiante components. The atoms in question are at the edge of the periodic box. I was just wondering why there is a discrepancy in the two utilities? Many thanks Gavin -- 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
Re: [gmx-users] discrepancy trjconv and gmxdump
Hi Mark My Apologies. Here is more info. gmxdump -f traj.trr history From history I take the coordinates at 250 ps and convert to gro file say test.gro trjconv -f traj.trr -dump 250 -o frame5.gro I then compare frame5.gro with test.gro Please find attached an excerpts of both files (in one file) and pay attention to the sign difference of the x coordinate at atom s 38 50 51 59 Cheers Gavin Mark Abraham wrote: On 2/11/2011 10:23 PM, Gavin Melaugh wrote: Hi all I have generated a gro file from a traj.trr file using trjconv. When I use gmxdump on the same traj.trr file to output a generic format history file it seems that there is a discrepancy in the coordinates of some atoms in a particular frame. Essentially I output frame 5 using trjconv, then I take the corresponding frame from the history file and covert it to gro format. I then run vimdiff on the two files. There are sign differences with the coordiante components. The atoms in question are at the edge of the periodic box. I was just wondering why there is a discrepancy in the two utilities? Without actual command lines that produced output and actual examples of what you think is anomalous we can't say. Mark History trajdump 240 240 1CGE CA1 0.287 4.485 6.3451CGE CA1 0.287 4.485 6.345 1CGE CB2 0.381 4.539 6.4421CGE CB2 0.381 4.539 6.442 1CGE HC3 0.182 4.477 6.3731CGE HC3 0.182 4.477 6.373 1CGE CA4 0.504 4.564 6.3911CGE CA4 0.504 4.564 6.391 1CGE CB5 0.551 4.515 6.2641CGE CB5 0.551 4.515 6.264 1CGE HC6 0.570 4.621 6.4461CGE HC6 0.570 4.621 6.446 1CGE CA7 0.467 4.440 6.1811CGE CA7 0.467 4.440 6.181 1CGE CB8 0.328 4.437 6.2181CGE CB8 0.328 4.437 6.218 1CGE HC9 0.506 4.388 6.0961CGE HC9 0.506 4.388 6.096 1CGE CU 10 0.336 4.595 6.5791CGE CU 10 0.336 4.595 6.579 1CGE NU 11 0.232 4.544 6.6311CGE NU 11 0.232 4.544 6.631 1CGE CH 12 0.161 4.574 6.7571CGE CH 12 0.161 4.574 6.757 1CGE CU 13 0.692 4.527 6.2171CGE CU 13 0.692 4.527 6.217 1CGE NU 14 0.754 4.485 6.1191CGE NU 14 0.754 4.485 6.119 1CGE CH 15 0.897 4.530 6.1271CGE CH 15 0.897 4.530 6.127 1CGE CU 16 0.231 4.395 6.1221CGE CU 16 0.231 4.395 6.122 1CGE NU 17 0.107 4.398 6.1581CGE NU 17 0.107 4.398 6.158 1CGE CH 18 0.002 4.359 6.0701CGE CH 18 0.002 4.359 6.070 1CGE CA 19 0.633 5.058 5.9601CGE CA 19 0.633 5.058 5.960 1CGE CB 20 0.690 4.968 6.0571CGE CB 20 0.690 4.968 6.057 1CGE HC 21 0.662 5.025 5.8541CGE HC 21 0.662 5.025 5.854 1CGE CA 22 0.677 4.992 6.1931CGE CA 22 0.677 4.992 6.193 1CGE CB 23 0.577 5.090 6.2311CGE CB 23 0.577 5.090 6.231 1CGE HC 24 0.714 4.916 6.2601CGE HC 24 0.714 4.916 6.260 1CGE CA 25 0.509 5.175 6.1401CGE CA 25 0.509 5.175 6.140 1CGE CB 26 0.536 5.158 5.9961CGE CB 26 0.536 5.158 5.996 1CGE HC 27 0.434 5.238 6.1781CGE HC 27 0.434 5.238 6.178 1CGE CU 28 0.773 4.848 6.0111CGE CU 28 0.773 4.848 6.011 1CGE NU 29 0.847 4.767 6.0831CGE NU 29 0.847 4.767 6.083 1CGE CH 30 0.922 4.651 6.0421CGE CH 30 0.922 4.651 6.041 1CGE CU 31 0.555 5.129 6.3741CGE CU 31 0.555 5.129 6.374 1CGE NU 32 0.626 5.070 6.4671CGE NU 32 0.626 5.070 6.467 1CGE CH 33 0.609 5.112 6.6071CGE CH 33 0.609 5.112 6.607 1CGE CU 34 0.466 5.249 5.9001CGE CU 34 0.466 5.249 5.900 1CGE NU 35 0.501 5.268 5.7861CGE NU 35 0.501 5.268 5.786 1CGE CH 36 0.430 5.358 5.6981CGE CH 36 0.430 5.358 5.698 1CGE CA 37 0.030 5.088 5.9021CGE CA
Re: [gmx-users] discrepancy trjconv and gmxdump
Actually Mark, I may have made a very trivial error. Forget about it for now. Cheers Gavin Gavin Melaugh wrote: Hi Mark My Apologies. Here is more info. gmxdump -f traj.trr history From history I take the coordinates at 250 ps and convert to gro file say test.gro trjconv -f traj.trr -dump 250 -o frame5.gro I then compare frame5.gro with test.gro Please find attached an excerpts of both files (in one file) and pay attention to the sign difference of the x coordinate at atom s 38 50 51 59 Cheers Gavin Mark Abraham wrote: On 2/11/2011 10:23 PM, Gavin Melaugh wrote: Hi all I have generated a gro file from a traj.trr file using trjconv. When I use gmxdump on the same traj.trr file to output a generic format history file it seems that there is a discrepancy in the coordinates of some atoms in a particular frame. Essentially I output frame 5 using trjconv, then I take the corresponding frame from the history file and covert it to gro format. I then run vimdiff on the two files. There are sign differences with the coordiante components. The atoms in question are at the edge of the periodic box. I was just wondering why there is a discrepancy in the two utilities? Without actual command lines that produced output and actual examples of what you think is anomalous we can't say. Mark -- 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] (no subject)
Hi all A very quick question. If I want to constrain all bonds in my system, is it just a case of setting constraints = all-bonds. Do I have to select the constraint algorithm. Here is my mpd file, I got no errors but its seems that this might be too simple. Note there is no relevance in the name umbrella here. Cheers Gavin umbrella1.mdp Description: application/mdp -- 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] bond constraints
Hi all Sorry forgot to give a subject in the first e-mail. A very quick question. If I want to constrain all bonds in my system, is it just a case of setting constraints = all-bonds. Do I have to select the constraint algorithm. Here is my mdp file, I got no errors but its seems that this might be too simple. Note there is no relevance in the name umbrella here. Cheers Gavin -- 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
Re: [gmx-users] atom types
Hi Justin Again thanks for the reply. I am not disagreeing with you but If I don't include a [pairs] directive in the topology file (with gen_pairs =yes), then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log file. When I include the [pair s] directive then both types of interaction are written to the log file. Therefore does gen_pairs= yes + [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to fudge LJ and QQ? Thanks Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Sorry to labour on this but: I don't quite understand what you mean when you say that nonbonded pair interactions are not Coulombic. Surely nonbonded charged atoms interact with each other, when close enough? (or by Nonbonded pair interactions do you explcitly mean 1-4,1-5, etc.) What are the 1-4 Coulombic interactions generated by, if not by gen_pairs =yes in my case? I have read this section of the manual loads and though I had a comprehensive understanding of it, but now I am confused again. Sure, there are nonbonded interactions for 1-4, 1-5, etc. But the purpose of [pairtype] generation is for LJ terms only. They are special 1-4 interactions between different atomtypes. Look at any force field for which [pairtypes] are listed - they have only C6 and C12 terms. Charges are not used for these calculations, but they are applied later during the MD using normal Coulombic equations and FudgeQQ. -Justin Thanks Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin I have checked the tpr file. Now it seems to assign the the two type of CHs as the same atom type, but at the same time with the specified charge from the [atoms] directive, as I expected. Concerning 1-4 interactions and gen_pairs =yes, my concern is this; from the pair list and using gen_pairs = yes, does grompp then take the 1-4 coulombic interaction for CH from the [atomtypes] directive (as the meaning of gen_pairs =yes)? Or does it assign the charge based on the atom index in the pair list? Charges are irrelevant for generation of pair interactions. Nonbonded pair interactions are LJ, not Coulombic. You will certainly have 1-4 Coulombic interactions, but they are not generated by gen_pairs. See manual section 5.3.4. -Justin Many Thanks Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi all A very quick question. I have an atom-type labelled CH in the atom-types with a particular charge, and in the atom list I assign some of these specific atoms with zero charge as below. When I generate 1,4 interactions using gen_pairs =yes, what charge for the CH type does it use? Does gromacs assign the CH with the different charge as a new atom type. Charges set in [atomtypes] are not used. The zero charge is assigned. Verify this by using gmxdump on your .tpr file. -Justin ;Parameter level [defaults] ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ 1 3 yes0.5 0.5 [atomtypes] ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 HC 1.008000 0.115000 A 0.242000 0.125520 CU 13.019000 0.265000 A 0.35 0.334720 NU 14.007000 -0.597000 A 0.325000 0.711280 CH 13.019000 0.332000 A 0.385000 0.334720 C3 15.035000 0.00 A 0.391000 0.669440 C2 14.027000 0.00 A 0.390500 0.493712 ;Molecular level [moleculetype] ; name nrexcl isotridecylcage 3 [atoms] . 72 CH 1 CGECH 24 0.3320 13.0190 73 C2 1 CGEC2 25 0. 14.0270 74 C2 1 CGEC2 25 0. 14.0270 75 C2 1 CGEC2 25 0. 14.0270 76 C2 1 CGEC2 26 0. 14.0270 77 C2 1 CGEC2 26 0. 14.0270 78 C2 1 CGEC2 26 0. 14.0270 79 C2 1 CGEC2 27 0. 14.0270 80 C2 1 CGEC2 27 0. 14.0270 81 C2 1 CGEC2 27 0. 14.0270 82 C2 1 CGEC2 28 0. 14.0270 83 CH 1 CGECH 28 0. 13.0190 84 C3 1 CGEC3 29 0. 15.0350 Many Thanks Gavin -- 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
Re: [gmx-users] atom types
Hi Mark Thanks for the reply. I am currently reading that section of the manual and, unless I am completely mistaken, it seems to vindicate what I am saying. Extra Lennard-Jones and electrostatic interactions between pairs of atoms in a molecule can be added in the [pairs] section of a molecule definition. In my [atom types] directive I have atomtype, charge mass, sigma and epsilon etc. All nonbonding parameters are then calculated according to the combination rule (in my case 3). 1-4 interactions are then calculated based on the information in [pairs] directive (all atoms are three bond away). I just have the atom indices of each pair in this directive therefore with gen_pairs = yes, the interaction parameters between each pair (which are 1-4) are calculated based on Fudge LJ and Fudge QQ (which are both 0.5 in my case). All of this in conjunction with nrexcl =3. Or am I completely wrong? In my set up then, are 1-4 Coulombic interactions determined by the pair list and fudge QQ? Many Thanks Gavin Mark Abraham wrote: On 02/08/11, *Gavin Melaugh * gmelaug...@qub.ac.uk wrote: Hi Justin Again thanks for the reply. I am not disagreeing with you but If I don't include a [pairs] directive in the topology file (with gen_pairs =yes), then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log file. When I include the [pair s] directive then both types of interaction are written to the log file. Therefore does gen_pairs= yes + [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to fudge LJ and QQ? Does manual section 5.3.4 answer your question? Mark Thanks Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Sorry to labour on this but: I don't quite understand what you mean when you say that nonbonded pair interactions are not Coulombic. Surely nonbonded charged atoms interact with each other, when close enough? (or by Nonbonded pair interactions do you explcitly mean 1-4,1-5, etc.) What are the 1-4 Coulombic interactions generated by, if not by gen_pairs =yes in my case? I have read this section of the manual loads and though I had a comprehensive understanding of it, but now I am confused again. Sure, there are nonbonded interactions for 1-4, 1-5, etc. But the purpose of [pairtype] generation is for LJ terms only. They are special 1-4 interactions between different atomtypes. Look at any force field for which [pairtypes] are listed - they have only C6 and C12 terms. Charges are not used for these calculations, but they are applied later during the MD using normal Coulombic equations and FudgeQQ. -Justin Thanks Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin I have checked the tpr file. Now it seems to assign the the two type of CHs as the same atom type, but at the same time with the specified charge from the [atoms] directive, as I expected. Concerning 1-4 interactions and gen_pairs =yes, my concern is this; from the pair list and using gen_pairs = yes, does grompp then take the 1-4 coulombic interaction for CH from the [atomtypes] directive (as the meaning of gen_pairs =yes)? Or does it assign the charge based on the atom index in the pair list? Charges are irrelevant for generation of pair interactions. Nonbonded pair interactions are LJ, not Coulombic. You will certainly have 1-4 Coulombic interactions, but they are not generated by gen_pairs. See manual section 5.3.4. -Justin Many Thanks Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi all A very quick question. I have an atom-type labelled CH in the atom-types with a particular charge, and in the atom list I assign some of these specific atoms with zero charge as below. When I generate 1,4 interactions using gen_pairs =yes, what charge for the CH type does it use? Does gromacs assign the CH with the different charge as a new atom type. Charges set in [atomtypes] are not used. The zero charge is assigned. Verify this by using gmxdump on your .tpr file. -Justin ;Parameter level [defaults] ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ 1 3 yes0.5 0.5 [atomtypes] ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 HC 1.008000 0.115000 A 0.242000 0.125520 CU 13.019000 0.265000 A 0.35 0.334720 NU 14.007000 -0.597000 A 0.325000 0.711280 CH 13.019000 0.332000 A 0.385000 0.334720 C3 15.035000 0.00 A 0.391000 0.669440 C2 14.027000 0.00 A 0.390500 0.493712 ;Molecular level [moleculetype
Re: [gmx-users] atom types
Yes I think the example vindicates what I am saying as well. I suppose I the contradiction ( I'll call it the point of confusion) you refer to is perhaps when Justin (who is always more than helpful) said that Charges are irrelevant for generation of pair interactions. Nonbonded pair interactions are LJ, not Coulombic. You will certainly have 1-4 Coulombic interactions, but they are not generated by gen_pairs. See manual section 5.3.4. My sequence of 1-4 interaction generation should go like this I suppose: e.g. [pairs] 3 6 no parameters present therefore get from [pairtypes] directive. no [pairtypes] directive therefore get from [non_bonded parameters] directive as gen pairs = yes again no [non_bonded parameters] directive. Therefore generate 1,4 interaction parameters based on the normal sigma and epsilon values (comb rule 3) present in [atomtypes] directive, in accordance with fudge LJ and QQ. My point is, then in conclusion, that in this way surely the 1,4 electrostatic interactions are determined by the pair list and in my case gen_pairs = yes no? Many Thanks Gavin Mark Abraham wrote: On 02/08/11, *Gavin Melaugh * gmelaug...@qub.ac.uk wrote: Hi Mark Thanks for the reply. I am currently reading that section of the manual and, unless I am completely mistaken, it seems to vindicate what I am saying. Extra Lennard-Jones and electrostatic interactions between pairs of atoms in a molecule can be added in the [pairs] section of a molecule definition. In my [atom types] directive I have atomtype, charge mass, sigma and epsilon etc. All nonbonding parameters are then calculated according to the combination rule (in my case 3). 1-4 interactions are then calculated based on the information in [pairs] directive (all atoms are three bond away). I just have the atom indices of each pair in this directive therefore with gen_pairs = yes, the interaction parameters between each pair (which are 1-4) are calculated based on Fudge LJ and Fudge QQ (which are both 0.5 in my case). All of this in conjunction with nrexcl =3. That will generate parameters for the interactions listed in [pairs] that do not have corresponding [pairtypes]. FudgeLJ and [nonbond_params] are used in such generation, per other parts of 5.7. Or am I completely wrong? In my set up then, are 1-4 Coulombic interactions determined by the pair list and fudge QQ? If the contradiction you think exists is this one... On 02/08/11, *Gavin Melaugh * gmelaug...@qub.ac.uk wrote: Hi Justin Again thanks for the reply. I am not disagreeing with you but If I don't include a [pairs] directive in the topology file (with gen_pairs =yes), then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log file. When I include the [pair s] directive then both types of interaction are written to the log file. Therefore does gen_pairs= yes + [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to fudge LJ and QQ? ... then 5.3.4 indicates that the presence of a [pairs] directive will generate the 1,4 output fields. The parameters for that output are taken from [pairtypes]. If gen-pairs=yes then the parameters are generated, else some warning/error occurs. The example in 5.7.1 has some more explanation about the use of the fudge parameters. Mark -- 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] to clear up the confusion regarding fudge LJ and QQ
Hi Justin and Mark Fudge LJ is only used when gen_pairs =yes. Fudge QQ is used regardless. I suppose a more concise and sensible question is the following: Are the atoms involved in the 1,4 Coulombic interactions specified in the [pairs] directive? If not how else are they generated. Cheers Gavin Gavin Melaugh wrote: Yes I think the example vindicates what I am saying as well. I suppose I the contradiction ( I'll call it the point of confusion) you refer to is perhaps when Justin (who is always more than helpful) said that Charges are irrelevant for generation of pair interactions. Nonbonded pair interactions are LJ, not Coulombic. You will certainly have 1-4 Coulombic interactions, but they are not generated by gen_pairs. See manual section 5.3.4. My sequence of 1-4 interaction generation should go like this I suppose: e.g. [pairs] 3 6 no parameters present therefore get from [pairtypes] directive. no [pairtypes] directive therefore get from [non_bonded parameters] directive as gen pairs = yes again no [non_bonded parameters] directive. Therefore generate 1,4 interaction parameters based on the normal sigma and epsilon values (comb rule 3) present in [atomtypes] directive, in accordance with fudge LJ and QQ. My point is, then in conclusion, that in this way surely the 1,4 electrostatic interactions are determined by the pair list and in my case gen_pairs = yes no? Many Thanks Gavin Mark Abraham wrote: On 02/08/11, *Gavin Melaugh * gmelaug...@qub.ac.uk wrote: Hi Mark Thanks for the reply. I am currently reading that section of the manual and, unless I am completely mistaken, it seems to vindicate what I am saying. Extra Lennard-Jones and electrostatic interactions between pairs of atoms in a molecule can be added in the [pairs] section of a molecule definition. In my [atom types] directive I have atomtype, charge mass, sigma and epsilon etc. All nonbonding parameters are then calculated according to the combination rule (in my case 3). 1-4 interactions are then calculated based on the information in [pairs] directive (all atoms are three bond away). I just have the atom indices of each pair in this directive therefore with gen_pairs = yes, the interaction parameters between each pair (which are 1-4) are calculated based on Fudge LJ and Fudge QQ (which are both 0.5 in my case). All of this in conjunction with nrexcl =3. That will generate parameters for the interactions listed in [pairs] that do not have corresponding [pairtypes]. FudgeLJ and [nonbond_params] are used in such generation, per other parts of 5.7. Or am I completely wrong? In my set up then, are 1-4 Coulombic interactions determined by the pair list and fudge QQ? If the contradiction you think exists is this one... On 02/08/11, *Gavin Melaugh * gmelaug...@qub.ac.uk wrote: Hi Justin Again thanks for the reply. I am not disagreeing with you but If I don't include a [pairs] directive in the topology file (with gen_pairs =yes), then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log file. When I include the [pair s] directive then both types of interaction are written to the log file. Therefore does gen_pairs= yes + [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to fudge LJ and QQ? ... then 5.3.4 indicates that the presence of a [pairs] directive will generate the 1,4 output fields. The parameters for that output are taken from [pairtypes]. If gen-pairs=yes then the parameters are generated, else some warning/error occurs. The example in 5.7.1 has some more explanation about the use of the fudge parameters. Mark -- 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
Re: [gmx-users] atom types
Hi Justin I see that we may have got our wires crossed from the off. Consider the [pairs] directive, which determines which atoms interact in a 1,4 manner. Consider two atoms listed the [pairs] directive. From the point of the Coulombic interaction between these two atoms I suppose my original question should have been: Does mdrun, when calculating the Coulombic potential between these two atoms, use the charges assigned to the atoms in the [atomtypes] directive or [atoms] directive ? Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Yes I think the example vindicates what I am saying as well. I suppose I the contradiction ( I'll call it the point of confusion) you refer to is perhaps when Justin (who is always more than helpful) said that Charges are irrelevant for generation of pair interactions. Nonbonded pair interactions are LJ, not Coulombic. You will certainly have 1-4 Coulombic interactions, but they are not generated by gen_pairs. See manual section 5.3.4. Charges *are* irrelevant - the information is not used when generating pairs, which I thought was the original question. The charge information is used during MD, when the pair list tells mdrun which atoms interact in what way. Then you get Coul. 1-4 terms. Perhaps I missed your point, but this whole thread started as which charge is used to generate pairs? The answer is still none. My sequence of 1-4 interaction generation should go like this I suppose: e.g. [pairs] 3 6 no parameters present therefore get from [pairtypes] directive. no [pairtypes] directive therefore get from [non_bonded parameters] directive as gen pairs = yes again no [non_bonded parameters] directive. Therefore generate 1,4 interaction parameters based on the normal sigma and epsilon values (comb rule 3) present in [atomtypes] directive, in accordance with fudge LJ and QQ. My point is, then in conclusion, that in this way surely the 1,4 electrostatic interactions are determined by the pair list and in my case gen_pairs = yes no? Yes. -Justin Many Thanks Gavin Mark Abraham wrote: On 02/08/11, *Gavin Melaugh * gmelaug...@qub.ac.uk wrote: Hi Mark Thanks for the reply. I am currently reading that section of the manual and, unless I am completely mistaken, it seems to vindicate what I am saying. Extra Lennard-Jones and electrostatic interactions between pairs of atoms in a molecule can be added in the [pairs] section of a molecule definition. In my [atom types] directive I have atomtype, charge mass, sigma and epsilon etc. All nonbonding parameters are then calculated according to the combination rule (in my case 3). 1-4 interactions are then calculated based on the information in [pairs] directive (all atoms are three bond away). I just have the atom indices of each pair in this directive therefore with gen_pairs = yes, the interaction parameters between each pair (which are 1-4) are calculated based on Fudge LJ and Fudge QQ (which are both 0.5 in my case). All of this in conjunction with nrexcl =3. That will generate parameters for the interactions listed in [pairs] that do not have corresponding [pairtypes]. FudgeLJ and [nonbond_params] are used in such generation, per other parts of 5.7. Or am I completely wrong? In my set up then, are 1-4 Coulombic interactions determined by the pair list and fudge QQ? If the contradiction you think exists is this one... On 02/08/11, *Gavin Melaugh * gmelaug...@qub.ac.uk wrote: Hi Justin Again thanks for the reply. I am not disagreeing with you but If I don't include a [pairs] directive in the topology file (with gen_pairs =yes), then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log file. When I include the [pair s] directive then both types of interaction are written to the log file. Therefore does gen_pairs= yes + [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to fudge LJ and QQ? ... then 5.3.4 indicates that the presence of a [pairs] directive will generate the 1,4 output fields. The parameters for that output are taken from [pairtypes]. If gen-pairs=yes then the parameters are generated, else some warning/error occurs. The example in 5.7.1 has some more explanation about the use of the fudge parameters. Mark -- 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
Re: [gmx-users] atom types
I hope this doesn't come across as stupid, or worse insolent. But what is the point in stating the charges in the atom type section then ? Gavin Mark Abraham wrote: On 02/08/11, *Gavin Melaugh * gmelaug...@qub.ac.uk wrote: Hi Justin I see that we may have got our wires crossed from the off. Consider the [pairs] directive, which determines which atoms interact in a 1,4 manner. Consider two atoms listed the [pairs] directive. From the point of the Coulombic interaction between these two atoms I suppose my original question should have been: Does mdrun, when calculating the Coulombic potential between these two atoms, use the charges assigned to the atoms in the [atomtypes] directive or [atoms] directive ? The charges from [atoms] are used. The charges in [atomtypes] are never used in any forcefield currently used with GROMACS. It might not be possible for grompp to ever use them, but I'd have to check the code for that. Mark Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Yes I think the example vindicates what I am saying as well. I suppose I the contradiction ( I'll call it the point of confusion) you refer to is perhaps when Justin (who is always more than helpful) said that Charges are irrelevant for generation of pair interactions. Nonbonded pair interactions are LJ, not Coulombic. You will certainly have 1-4 Coulombic interactions, but they are not generated by gen_pairs. See manual section 5.3.4. Charges *are* irrelevant - the information is not used when generating pairs, which I thought was the original question. The charge information is used during MD, when the pair list tells mdrun which atoms interact in what way. Then you get Coul. 1-4 terms. Perhaps I missed your point, but this whole thread started as which charge is used to generate pairs? The answer is still none. My sequence of 1-4 interaction generation should go like this I suppose: e.g. [pairs] 3 6 no parameters present therefore get from [pairtypes] directive. no [pairtypes] directive therefore get from [non_bonded parameters] directive as gen pairs = yes again no [non_bonded parameters] directive. Therefore generate 1,4 interaction parameters based on the normal sigma and epsilon values (comb rule 3) present in [atomtypes] directive, in accordance with fudge LJ and QQ. My point is, then in conclusion, that in this way surely the 1,4 electrostatic interactions are determined by the pair list and in my case gen_pairs = yes no? Yes. -Justin Many Thanks Gavin Mark Abraham wrote: On 02/08/11, *Gavin Melaugh * gmelaug...@qub.ac.uk wrote: Hi Mark Thanks for the reply. I am currently reading that section of the manual and, unless I am completely mistaken, it seems to vindicate what I am saying. Extra Lennard-Jones and electrostatic interactions between pairs of atoms in a molecule can be added in the [pairs] section of a molecule definition. In my [atom types] directive I have atomtype, charge mass, sigma and epsilon etc. All nonbonding parameters are then calculated according to the combination rule (in my case 3). 1-4 interactions are then calculated based on the information in [pairs] directive (all atoms are three bond away). I just have the atom indices of each pair in this directive therefore with gen_pairs = yes, the interaction parameters between each pair (which are 1-4) are calculated based on Fudge LJ and Fudge QQ (which are both 0.5 in my case). All of this in conjunction with nrexcl =3. That will generate parameters for the interactions listed in [pairs] that do not have corresponding [pairtypes]. FudgeLJ and [nonbond_params] are used in such generation, per other parts of 5.7. Or am I completely wrong? In my set up then, are 1-4 Coulombic interactions determined by the pair list and fudge QQ? If the contradiction you think exists is this one... On 02/08/11, *Gavin Melaugh * gmelaug...@qub.ac.uk wrote: Hi Justin Again thanks for the reply. I am not disagreeing with you but If I don't include a [pairs] directive in the topology file (with gen_pairs =yes), then there are no 1-4 LJ nor 1-4 Coulombic energies written in the log file. When I include the [pair s] directive then both types of interaction are written to the log file. Therefore does gen_pairs= yes + [pairs] directive generate 1,4 LJ and 1,4 Coulomb according to fudge LJ and QQ? ... then 5.3.4 indicates that the presence of a [pairs] directive will generate the 1,4 output fields. The parameters for that output are taken from [pairtypes]. If gen-pairs=yes then the parameters are generated, else some warning/error occurs. The example in 5.7.1 has some more explanation about the use of the fudge parameters. Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo
[gmx-users] atom types
Hi all A very quick question. I have an atom-type labelled CH in the atom-types with a particular charge, and in the atom list I assign some of these specific atoms with zero charge as below. When I generate 1,4 interactions using gen_pairs =yes, what charge for the CH type does it use? Does gromacs assign the CH with the different charge as a new atom type. ;Parameter level [defaults] ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ 1 3 yes0.5 0.5 [atomtypes] ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 HC 1.008000 0.115000 A 0.242000 0.125520 CU 13.019000 0.265000 A 0.35 0.334720 NU 14.007000 -0.597000 A 0.325000 0.711280 CH 13.019000 0.332000 A 0.385000 0.334720 C3 15.035000 0.00 A 0.391000 0.669440 C2 14.027000 0.00 A 0.390500 0.493712 ;Molecular level [moleculetype] ; name nrexcl isotridecylcage 3 [atoms] . 72 CH 1 CGECH 24 0.3320 13.0190 73 C2 1 CGEC2 25 0. 14.0270 74 C2 1 CGEC2 25 0. 14.0270 75 C2 1 CGEC2 25 0. 14.0270 76 C2 1 CGEC2 26 0. 14.0270 77 C2 1 CGEC2 26 0. 14.0270 78 C2 1 CGEC2 26 0. 14.0270 79 C2 1 CGEC2 27 0. 14.0270 80 C2 1 CGEC2 27 0. 14.0270 81 C2 1 CGEC2 27 0. 14.0270 82 C2 1 CGEC2 28 0. 14.0270 83 CH 1 CGECH 28 0. 13.0190 84 C3 1 CGEC3 29 0. 15.0350 Many Thanks Gavin -- 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
Re: [gmx-users] atom types
Hi Justin I have checked the tpr file. Now it seems to assign the the two type of CHs as the same atom type, but at the same time with the specified charge from the [atoms] directive, as I expected. Concerning 1-4 interactions and gen_pairs =yes, my concern is this; from the pair list and using gen_pairs = yes, does grompp then take the 1-4 coulombic interaction for CH from the [atomtypes] directive (as the meaning of gen_pairs =yes)? Or does it assign the charge based on the atom index in the pair list? Many Thanks Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi all A very quick question. I have an atom-type labelled CH in the atom-types with a particular charge, and in the atom list I assign some of these specific atoms with zero charge as below. When I generate 1,4 interactions using gen_pairs =yes, what charge for the CH type does it use? Does gromacs assign the CH with the different charge as a new atom type. Charges set in [atomtypes] are not used. The zero charge is assigned. Verify this by using gmxdump on your .tpr file. -Justin ;Parameter level [defaults] ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ 1 3 yes0.5 0.5 [atomtypes] ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 HC 1.008000 0.115000 A 0.242000 0.125520 CU 13.019000 0.265000 A 0.35 0.334720 NU 14.007000 -0.597000 A 0.325000 0.711280 CH 13.019000 0.332000 A 0.385000 0.334720 C3 15.035000 0.00 A 0.391000 0.669440 C2 14.027000 0.00 A 0.390500 0.493712 ;Molecular level [moleculetype] ; name nrexcl isotridecylcage 3 [atoms] . 72 CH 1 CGECH 24 0.3320 13.0190 73 C2 1 CGEC2 25 0. 14.0270 74 C2 1 CGEC2 25 0. 14.0270 75 C2 1 CGEC2 25 0. 14.0270 76 C2 1 CGEC2 26 0. 14.0270 77 C2 1 CGEC2 26 0. 14.0270 78 C2 1 CGEC2 26 0. 14.0270 79 C2 1 CGEC2 27 0. 14.0270 80 C2 1 CGEC2 27 0. 14.0270 81 C2 1 CGEC2 27 0. 14.0270 82 C2 1 CGEC2 28 0. 14.0270 83 CH 1 CGECH 28 0. 13.0190 84 C3 1 CGEC3 29 0. 15.0350 Many Thanks Gavin -- 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
Re: [gmx-users] atom types
Hi Justin Sorry to labour on this but: I don't quite understand what you mean when you say that nonbonded pair interactions are not Coulombic. Surely nonbonded charged atoms interact with each other, when close enough? (or by Nonbonded pair interactions do you explcitly mean 1-4,1-5, etc.) What are the 1-4 Coulombic interactions generated by, if not by gen_pairs =yes in my case? I have read this section of the manual loads and though I had a comprehensive understanding of it, but now I am confused again. Thanks Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin I have checked the tpr file. Now it seems to assign the the two type of CHs as the same atom type, but at the same time with the specified charge from the [atoms] directive, as I expected. Concerning 1-4 interactions and gen_pairs =yes, my concern is this; from the pair list and using gen_pairs = yes, does grompp then take the 1-4 coulombic interaction for CH from the [atomtypes] directive (as the meaning of gen_pairs =yes)? Or does it assign the charge based on the atom index in the pair list? Charges are irrelevant for generation of pair interactions. Nonbonded pair interactions are LJ, not Coulombic. You will certainly have 1-4 Coulombic interactions, but they are not generated by gen_pairs. See manual section 5.3.4. -Justin Many Thanks Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi all A very quick question. I have an atom-type labelled CH in the atom-types with a particular charge, and in the atom list I assign some of these specific atoms with zero charge as below. When I generate 1,4 interactions using gen_pairs =yes, what charge for the CH type does it use? Does gromacs assign the CH with the different charge as a new atom type. Charges set in [atomtypes] are not used. The zero charge is assigned. Verify this by using gmxdump on your .tpr file. -Justin ;Parameter level [defaults] ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ 1 3 yes0.5 0.5 [atomtypes] ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 HC 1.008000 0.115000 A 0.242000 0.125520 CU 13.019000 0.265000 A 0.35 0.334720 NU 14.007000 -0.597000 A 0.325000 0.711280 CH 13.019000 0.332000 A 0.385000 0.334720 C3 15.035000 0.00 A 0.391000 0.669440 C2 14.027000 0.00 A 0.390500 0.493712 ;Molecular level [moleculetype] ; name nrexcl isotridecylcage 3 [atoms] . 72 CH 1 CGECH 24 0.3320 13.0190 73 C2 1 CGEC2 25 0. 14.0270 74 C2 1 CGEC2 25 0. 14.0270 75 C2 1 CGEC2 25 0. 14.0270 76 C2 1 CGEC2 26 0. 14.0270 77 C2 1 CGEC2 26 0. 14.0270 78 C2 1 CGEC2 26 0. 14.0270 79 C2 1 CGEC2 27 0. 14.0270 80 C2 1 CGEC2 27 0. 14.0270 81 C2 1 CGEC2 27 0. 14.0270 82 C2 1 CGEC2 28 0. 14.0270 83 CH 1 CGECH 28 0. 13.0190 84 C3 1 CGEC3 29 0. 15.0350 Many Thanks Gavin -- 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] lattice Coulomb energy between energy groups
Hi all I have set up a simulation with PBC using PME and specify an energy group enclusion involving some virtual sites and a hydrocarbon chain. i.e. hydrocarbon chain is one energy group and the virtual sites constitute the other group. Grompp generates the following warning Can not exclude the lattice Coulomb energy between energy groups I have checked the mailing list and the manual but there doesn't seem to be any concise answers. Given the fact that both groups are not charged is it O.K. to ignore this warning. Many Thanks Gavin -- 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
Re: [gmx-users] lattice Coulomb energy between energy groups
Hi Justin Thanks for the reply. I take it then, if I use the maxwarn option, it will still apply the exclusion for the short-range non-bonding interactions. Cheers Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi all I have set up a simulation with PBC using PME and specify an energy group enclusion involving some virtual sites and a hydrocarbon chain. i.e. hydrocarbon chain is one energy group and the virtual sites constitute the other group. Grompp generates the following warning Can not exclude the lattice Coulomb energy between energy groups I have checked the mailing list and the manual but there doesn't seem to be any concise answers. Given the fact that both groups are not charged is it O.K. to ignore this warning. Probably. The warning indicates that one cannot use energygrp_excl to exclude components of the PME mesh; such exclusions only apply to short-range interactions that can be decomposed pairwise. Since the PME mesh is dependent on all the entities in the system, there is not way to simply calculate the mesh, except for some interactions between groups A and B. If there are no charges on groups A and B, you should be fine. -Justin -- 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] convergence
Hi all I have generated a PMF curve over 15 ns. Does g_wham have a facility whereby I can calculate the PMF over say 7 ns, to check for convergence. There doesn't seem to be anything in the manual. Many thanks Gavin -- 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: convergence
Please ignore my last question, I have found the answer using the g_wham -h option Cheers Gavin Gavin Melaugh wrote: Hi all I have generated a PMF curve over 15 ns. Does g_wham have a facility whereby I can calculate the PMF over say 7 ns, to check for convergence. There doesn't seem to be anything in the manual. Many thanks Gavin -- 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] error bars g_wham
Hi all I have read the manual and the recent JCTC paper on g_wham, and I was wondering how to actually get the error bars on the profile.xvg file outputted from g_wham. Many Thanks Gavin -- 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
Re: [gmx-users] Re: distance/direction PMF
Hi Justin I think what is happening is that grompp is returning the projection of the current vector between the two reference groups onto the stated vector in the mdp file. i.e the dot product of the two vectors. The resulting PMF plot must therefore be a function of this varying cross product. I say varying as I don't think the dynamics constrain the pull group alone that vector alone. When you change the vector in the mdp file to be that of the new vector connecting the two groups grompp gives the actual distance between the two groups. Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Yeah I set pull_init =0, and pull_start = yes, and it still gave a distance at start as 0.78nm. It must be doing something funny because with pull_vec, because when I use pull_geometry = distance and pull_dim = Y Y Y, the distance 0.815nm is returned as the distance at start, which is the actual distance between the two groups. I have no idea why this is happening. If you'd like me to troubleshoot further and see if I can get to the bottom of this, please send me (off-list is OK): 1. Coordinate file 2. Topology file(s) 3. The .mdp file that is giving the weird result -Justin Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Sorry maybe I was unclear before with my first point; the specified distance in the mdp file is 0.78nm, which grompp returns back as the ref distance. The distance at the start that grompp returns is also 0.78nm, even though the actual distance is 0.815nm. That's why I was asking what this distance actually means, as it is not the absolute distance between the reference groups. As I said before, you're telling grompp to make that the reference distance. You say the specified distance in the mdp file is 0.78 nm, therefore grompp is doing what it's told. If that's not the desired distance, then don't set it as such. Have you tried the combination I suggested in the last message? Does it give an initial distance of 0.815 nm? Also as regards to your second point; what if in the initial configurations the two reference groups do not lie along the stated vector? Hard to say, but probably you'll get some very high forces at the outset of the simulation as the simulation as the umbrella potential tries to force the pulled group to conform to the restraint specifications. Whether or not that negatively impacts the ultimate result of the simulation is also an important consideration. If you're pulling along a given vector, you should start with your reference coordinates as close to the desired location as possible. -Justin Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Justin I did a short test on a particular window with the specified vector with pull_init =0.78nm. I then grompp(ed) the final configuration and it gave me a distance of 0.78 as the initial ditsance, fair enough. However when I viewed the final line from g_dist the absolute distance or modulus was 0.815 nm, which agreed with the distance I calculated using a molecular configuration editor. My question is therefore this, when using pull_geometry =direction with a specified vector, How should one interpret the initial distance provided by grompp? and how does Gromacs deal with the fact that the vector between the two point changes during the run ? You're telling grompp that the initial distance is 0.78 nm, so it's spitting that back out. If that's not the true distance, then specify the correct quantity :) The distance should be correctly detected with: pull_start = yes pull_init = 0 As for the second question, the vector doesn't change over the simulation, but the position of pull_group1 along that vector will. That's what the umbrella potential is doing - it allows for harmonic oscillation along any of the dimensions specified by pull_vec/pull_dim (depending on the settings). -Justin -- 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] fatal error
Hi all Why is pull geometry direction not supported in g_wham ? I got the following error. I did this to compare with the same simulations but with pull geometry = distance. Fatal error: Pull geometry direction not supported -- 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] distance/direction PMF
Hi all Could someone please get back to me on this. I have ran two sets of umbrella sampling simulations 1- Using pull_geometry = distance, and pull_dim = Y Y Y 2- Using pull_geometry = direction, with pull_vec = 0.462808 0.494125 0.735968 In both cases I wish to calculate the PMF for taking a host out of a guest. I specify the distance between the COM of the two groups in the respective mdp files. The two curves I get are completely different which has leads me to ask a couple of questions: 1) Does the vector change as the dynamics evolves (i.e is it just relative between the two points) or does it remain fixed in space? 2) Using umbrella sampling in the second case does the relative position of the pull group (guest) fluctuate about the distance along the specified vector? Note- I couldn't run g_wham with pull_geometry = direction in my current version of gromacs, so I ran it in 4.5.3 but I had to use the -if pullf-files.xvg option. Some help would be appreciated as I have to be sure of what the program is actually doing. Many thanks Gavin -- 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
Re: [gmx-users] Re: distance/direction PMF
Justin I did a short test on a particular window with the specified vector with pull_init =0.78nm. I then grompp(ed) the final configuration and it gave me a distance of 0.78 as the initial ditsance, fair enough. However when I viewed the final line from g_dist the absolute distance or modulus was 0.815 nm, which agreed with the distance I calculated using a molecular configuration editor. My question is therefore this, when using pull_geometry =direction with a specified vector, How should one interpret the initial distance provided by grompp? and how does Gromacs deal with the fact that the vector between the two point changes during the run ? Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Many thanks for the reply I know technically I am doing two different things. However due to my starting conifgurations in both cases I would expect similar results. Could you confirm the following points 1) In pull_geometry=distance, and pull_dim =YYY; is the absolute distance between the COM of both groups in all directions considered ? Yes. 2) In pull_geometry = direction, with a specified vector; I assume that it is the distance between the two groups along this vector that is considered? Should be. Compare the output of g_dist with what grompp reports as the initial distance. 3) Is the specified vector taken to be the vector connecting the two points? Therefore are the desired distances in the mdp files essentially the magnitudes of this vector at each window. Per the documentation, grompp normalizes the vector, so no, not directly. But if you have proper settings for pull_init/pull_start, then the reference distance should be correctly calculated along the specified vector for each window. -Justin -- 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
Re: [gmx-users] Re: distance/direction PMF
Hi Justin Sorry maybe I was unclear before with my first point; the specified distance in the mdp file is 0.78nm, which grompp returns back as the ref distance. The distance at the start that grompp returns is also 0.78nm, even though the actual distance is 0.815nm. That's why I was asking what this distance actually means, as it is not the absolute distance between the reference groups. Also as regards to your second point; what if in the initial configurations the two reference groups do not lie along the stated vector? Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Justin I did a short test on a particular window with the specified vector with pull_init =0.78nm. I then grompp(ed) the final configuration and it gave me a distance of 0.78 as the initial ditsance, fair enough. However when I viewed the final line from g_dist the absolute distance or modulus was 0.815 nm, which agreed with the distance I calculated using a molecular configuration editor. My question is therefore this, when using pull_geometry =direction with a specified vector, How should one interpret the initial distance provided by grompp? and how does Gromacs deal with the fact that the vector between the two point changes during the run ? You're telling grompp that the initial distance is 0.78 nm, so it's spitting that back out. If that's not the true distance, then specify the correct quantity :) The distance should be correctly detected with: pull_start = yes pull_init = 0 As for the second question, the vector doesn't change over the simulation, but the position of pull_group1 along that vector will. That's what the umbrella potential is doing - it allows for harmonic oscillation along any of the dimensions specified by pull_vec/pull_dim (depending on the settings). -Justin -- 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
Re: [gmx-users] Re: distance/direction PMF
Yeah I set pull_init =0, and pull_start = yes, and it still gave a distance at start as 0.78nm. It must be doing something funny because with pull_vec, because when I use pull_geometry = distance and pull_dim = Y Y Y, the distance 0.815nm is returned as the distance at start, which is the actual distance between the two groups. Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Sorry maybe I was unclear before with my first point; the specified distance in the mdp file is 0.78nm, which grompp returns back as the ref distance. The distance at the start that grompp returns is also 0.78nm, even though the actual distance is 0.815nm. That's why I was asking what this distance actually means, as it is not the absolute distance between the reference groups. As I said before, you're telling grompp to make that the reference distance. You say the specified distance in the mdp file is 0.78 nm, therefore grompp is doing what it's told. If that's not the desired distance, then don't set it as such. Have you tried the combination I suggested in the last message? Does it give an initial distance of 0.815 nm? Also as regards to your second point; what if in the initial configurations the two reference groups do not lie along the stated vector? Hard to say, but probably you'll get some very high forces at the outset of the simulation as the simulation as the umbrella potential tries to force the pulled group to conform to the restraint specifications. Whether or not that negatively impacts the ultimate result of the simulation is also an important consideration. If you're pulling along a given vector, you should start with your reference coordinates as close to the desired location as possible. -Justin Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Justin I did a short test on a particular window with the specified vector with pull_init =0.78nm. I then grompp(ed) the final configuration and it gave me a distance of 0.78 as the initial ditsance, fair enough. However when I viewed the final line from g_dist the absolute distance or modulus was 0.815 nm, which agreed with the distance I calculated using a molecular configuration editor. My question is therefore this, when using pull_geometry =direction with a specified vector, How should one interpret the initial distance provided by grompp? and how does Gromacs deal with the fact that the vector between the two point changes during the run ? You're telling grompp that the initial distance is 0.78 nm, so it's spitting that back out. If that's not the true distance, then specify the correct quantity :) The distance should be correctly detected with: pull_start = yes pull_init = 0 As for the second question, the vector doesn't change over the simulation, but the position of pull_group1 along that vector will. That's what the umbrella potential is doing - it allows for harmonic oscillation along any of the dimensions specified by pull_vec/pull_dim (depending on the settings). -Justin -- 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] energy conservation
Hi all I have run a very long NVT simulation 500ns of one molecule (220 atoms) to try to assess the fluctuations in a particular event. When I analyse the energies it seems that the conserved quantity drifts significantly. Why would this be? Has it got to do with the thermostat I am using? Please find below the mdp file. title = Pull test cpp = include = define = integrator = md nsteps = 3 dt = 0.001 nstxout = 15 nstvout = 15 nstlog = 15 nstenergy = 15000 nstfout = 15 pbc = no nstlist = 0 ns_type = simple vdwtype = cut-off rlist = 0 rvdw_switch = 0 rvdw= 0 coulombtype = cut-off rcoulomb= 0 tcoupl = nose-hoover tc_grps = system tau_t = 0.1 ref_t = 400 gen_vel = no gen_temp= constraints = none comm_mode = angular Cheers Gavin -- 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] pull_vec
Hi all In umbrella sampling simulations: If pull_geometry = direction, and pull_vec is obviously a specified vector, does the distance between the two groups only vary along that vector? i.e Is the umbrella potential only defined along that vector? Many Thanks Gavin -- 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] Energy groups
Hi all If I have a topology with set up below (excerpts). And define energy groups in the .mdp file say 'virtsites' and 'endgroup' in correspondence with the index file. Will defining energygrp_exl 'virtsites' and 'endgroup' prevent nonbonding interactions between any atom in either of these groups. Also will all other nonbonding interactions stated in the topology remain the same? Many Thanks Gavin ;Parameter level [defaults] ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ 1 3 yes0.5 0.5 [atomtypes] ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 HC 1.008000 0.115000 A 0.242000 0.125520 CU 13.019000 0.265000 A 0.35 0.334720 NU 14.007000 -0.597000 A 0.325000 0.711280 CH 13.019000 0.332000 A 0.385000 0.334720 C3 15.035000 0.00 A 0.391000 0.669440 C2 14.027000 0.00 A 0.390500 0.493712 VS 0.0 0.0V 0.0 0.0 [nonbond_params] ;i j funct sigma epsilon VS C3 1 0.10.03153 -- 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] exclusions
Hi all
I am cheking the non binary version topol.tpr form grompp and I was
wondering what the following format for the exclusion section signifies?
excls:
nr=228
nra=1920
excls[0][0..12]={0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 16}
excls[1][13..25]={0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 15}
Also I take it, that this is a result of nrexcl =3
Cheers
Gavin
--
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
Re: [gmx-users] exclusions
Regarding my earlier post- If someone could confirm.
Taking into account that atom indices start from zero.
Then excls[0] = exclusions centred around atom 1.
[0..12] = there are 13 exclusions each one labelled form 0-12.
I then assume that you read the exclusions with respect to atom 0
(actual index =1)
In the next set do you read the exclusion with respect to atom 1 (actual
index =2).
and so on
Many thanks in advance
Gavin
Gavin Melaugh wrote:
Hi all
I am cheking the non binary version topol.tpr form grompp and I was
wondering what the following format for the exclusion section signifies?
excls:
nr=228
nra=1920
excls[0][0..12]={0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 16}
excls[1][13..25]={0, 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 15}
Also I take it, that this is a result of nrexcl =3
Cheers
Gavin
--
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] energy group error
Hi All I have set up a topology file and with virtual sites (exerpts below), whereby my virtual sites only interact with the C3 atoms. I am trying to construct a PMF wrt one particular C3. This atom I do not want to interact with the virtual sites at all. I thought the best way to do this was via energy exclusion groups. So I set up some groups in the index file (see below), and defined the energy groups like so energygrps = vsites vsitex energygrp_excl = vsites vsitex but I get the following error Fatal error: atoms 415 and 416 in charge group 89 of molecule type 'name' are in different energy groups I have no charge group 89 in my molecule. Can someone please suggest where I might be going wrong. ;Parameter level [defaults] ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ 1 3 yes0.5 0.5 [atomtypes] ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 HC 1.008000 0.115000 A 0.242000 0.125520 CU 13.019000 0.265000 A 0.35 0.334720 NU 14.007000 -0.597000 A 0.325000 0.711280 CH 13.019000 0.332000 A 0.385000 0.334720 C3 15.035000 0.00 A 0.391000 0.669440 C2 14.027000 0.00 A 0.390500 0.493712 VS 0.0 0.0V 0.0 0.0 [nonbond_params] ;i j funct sigma epsilon VS C3 1 0.10.03153 [ cage_1 ] 123456789 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 75 89 103 117 131 145 159 173 187 201 215 [ tail ] 416 [ vsites ] 229 230 231 232 461 462 463 464 [ vsitex ] 416 -- 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
Re: [gmx-users] energy group error
Hi again To rephrase. If I have two atoms constructing a charge group, is it not possible to have one of these atoms constructing an energy group on its own? Cheers Gavin Gavin Melaugh wrote: Hi All I have set up a topology file and with virtual sites (exerpts below), whereby my virtual sites only interact with the C3 atoms. I am trying to construct a PMF wrt one particular C3. This atom I do not want to interact with the virtual sites at all. I thought the best way to do this was via energy exclusion groups. So I set up some groups in the index file (see below), and defined the energy groups like so energygrps = vsites vsitex energygrp_excl = vsites vsitex but I get the following error Fatal error: atoms 415 and 416 in charge group 89 of molecule type 'name' are in different energy groups I have no charge group 89 in my molecule. Can someone please suggest where I might be going wrong. ;Parameter level [defaults] ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ 1 3 yes0.5 0.5 [atomtypes] ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 HC 1.008000 0.115000 A 0.242000 0.125520 CU 13.019000 0.265000 A 0.35 0.334720 NU 14.007000 -0.597000 A 0.325000 0.711280 CH 13.019000 0.332000 A 0.385000 0.334720 C3 15.035000 0.00 A 0.391000 0.669440 C2 14.027000 0.00 A 0.390500 0.493712 VS 0.0 0.0V 0.0 0.0 [nonbond_params] ;i j funct sigma epsilon VS C3 1 0.10.03153 [ cage_1 ] 123456789 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 75 89 103 117 131 145 159 173 187 201 215 [ tail ] 416 [ vsites ] 229 230 231 232 461 462 463 464 [ vsitex ] 416 -- 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] charge groups
Hi all Do atoms belonging to the same charge group have to be indexed consecutively in the topology file or can you have the following. [atoms] 1 CA 1 CGECA 1 -0.1150 12.0110 2 CB 1 CGECB 1 0. 12.0110 3 CA 1 CGECA 2 -0.1150 12.0110 4 CB 1 CGE :wq CB 2 0. 12.0110 5 CA 1 CGECA 3 -0.1150 12.0110 6 CB 1 CGECB 3 0. 12.0110 7 HC 1 CGEHC 1 0.1150 1.0080 8 CU 1 CGECU 13 0.2650 13.0190 9 NU 1 CGENU 13 -0.5970 14.0070 10 HC 1 CGEHC 2 0.1150 1.0080 11 CU 1 CGECU 14 0.2650 13.0190 12 NU 1 CGENU 14 -0.5970 14.0070 13 HC 1 CGEHC 3 0.1150 1.0080 Cheers Gavin -- 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
Re: [gmx-users] charge groups
As you can see atomnr 7 HC should belong to the same charge group (1) as
atomnr 12. From the tpr file atomnr 7 (in this file 6) seems to have
its own charge group (3).
[atoms]
; atomnr type resnr residue namecgnr charge mass
1 CA 1 CGECA 1 -0.1150 12.0110
2 CB 1 CGECB 1 0. 12.0110
3 CA 1 CGECA 2 -0.1150 12.0110
4 CB 1 CGECB 2 0. 12.0110
5 CA 1 CGECA 3 -0.1150 12.0110
6 CB 1 CGECB 3 0. 12.0110
7 HC 1 CGEHC 1 0.1150 1.0080
8 CU 1 CGECU 13 0.2650 13.0190
9 NU 1 CGENU 13 -0.5970 14.0070
10 HC 1 CGEHC 2 0.1150 1.0080
11 CU 1 CGECU 14 0.2650 13.0190
12 NU 1 CGENU 14 -0.5970 14.0070
13 HC 1 CGEHC 3 0.1150 1.0080
14 CU 1 CGECU 15 0.2650 13.0190
15 NU 1 CGENU 15 -0.5970 14.0070
cgs:
nr=112
cgs[0]={0..1}
cgs[1]={2..3}
cgs[2]={4..5}
cgs[3]={6..6}
cgs[4]={7..8}
cgs[5]={9..9}
cgs[6]={10..11}
cgs[7]={12..12}
cgs[8]={13..14}
cgs[9]={15..16}
cgs[10]={17..18}
cgs[11]={19..20}
cgs[12]={21..21}
cgs[13]={22..23}
cgs[14]={24..24}
cgs[15]={25..26}
cgs[16]={27..27}
cgs[17]={28..29}
cgs[18]={30..31}
cgs[19]={32..33}
cgs[20]={34..35}
cgs[21]={36..36}
cgs[22]={37..38}
cgs[23]={39..39}
cgs[24]={40..41}
cgs[25]={42..42}
cgs[26]={43..44}
cgs[27]={45..46}
Justin A. Lemkul wrote:
Gavin Melaugh wrote:
Hi Justin
The two files processed.top and topol.top are the exact same. What is
worrying me is that when I look at the topol.tpr from gmxdump, my charge
groups do not seem to be specified the way I stated in the topology
file.
Can you provide a relevant snippet of the gmxdump output to
demonstrate this?
Consecutive numbering is not stated as a requirement in the manual,
and if it is in the code, then it needs to be either documented as
being necessary, and/or an enhancement request needs to be filed on
redmine for non-consecutive charge group numbering.
-Justin
Gavin
Justin A. Lemkul wrote:
Gavin Melaugh wrote:
Hi all
Do atoms belonging to the same charge group have to be indexed
consecutively in the topology file or can you have the following.
[atoms]
1 CA 1 CGECA 1 -0.1150 12.0110
2 CB 1 CGECB 1 0. 12.0110
3 CA 1 CGECA 2 -0.1150 12.0110
4 CB 1 CGE :wq
CB 2 0. 12.0110
Assuming you fix this, then yes, splitting charge groups can work.
Use grompp -pp to see the post-processed topology to make sure
everything came out as expected.
-Justin
5 CA 1 CGECA 3 -0.1150 12.0110
6 CB 1 CGECB 3 0. 12.0110
7 HC 1 CGEHC 1 0.1150 1.0080
8 CU 1 CGECU 13 0.2650 13.0190
9 NU 1 CGENU 13 -0.5970 14.0070
10 HC 1 CGEHC 2 0.1150 1.0080
11 CU 1 CGECU 14 0.2650 13.0190
12 NU 1 CGENU 14 -0.5970 14.0070
13 HC 1 CGEHC 3 0.1150 1.0080
Cheers
Gavin
--
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
Re: [gmx-users] charge groups
O.K cheers
I assume that this is irrelevant in vacuum with no cut-offs
Gavin
Justin A. Lemkul wrote:
Gavin Melaugh wrote:
As you can see atomnr 7 HC should belong to the same charge group (1) as
atomnr 12. From the tpr file atomnr 7 (in this file 6) seems to have
its own charge group (3).
[atoms]
; atomnr type resnr residue namecgnr charge mass
1 CA 1 CGECA 1 -0.1150 12.0110
2 CB 1 CGECB 1 0. 12.0110
3 CA 1 CGECA 2 -0.1150 12.0110
4 CB 1 CGECB 2 0. 12.0110
5 CA 1 CGECA 3 -0.1150 12.0110
6 CB 1 CGECB 3 0. 12.0110
7 HC 1 CGEHC 1 0.1150 1.0080
8 CU 1 CGECU 13 0.2650 13.0190
9 NU 1 CGENU 13 -0.5970 14.0070
10 HC 1 CGEHC 2 0.1150 1.0080
11 CU 1 CGECU 14 0.2650 13.0190
12 NU 1 CGENU 14 -0.5970 14.0070
13 HC 1 CGEHC 3 0.1150 1.0080
14 CU 1 CGECU 15 0.2650 13.0190
15 NU 1 CGENU 15 -0.5970 14.0070
cgs:
nr=112
cgs[0]={0..1}
cgs[1]={2..3}
cgs[2]={4..5}
cgs[3]={6..6}
cgs[4]={7..8}
cgs[5]={9..9}
cgs[6]={10..11}
cgs[7]={12..12}
cgs[8]={13..14}
cgs[9]={15..16}
cgs[10]={17..18}
cgs[11]={19..20}
cgs[12]={21..21}
cgs[13]={22..23}
cgs[14]={24..24}
cgs[15]={25..26}
cgs[16]={27..27}
cgs[17]={28..29}
cgs[18]={30..31}
cgs[19]={32..33}
cgs[20]={34..35}
cgs[21]={36..36}
cgs[22]={37..38}
cgs[23]={39..39}
cgs[24]={40..41}
cgs[25]={42..42}
cgs[26]={43..44}
cgs[27]={45..46}
Then it seems pretty clear that one cannot have discontinuous charge
groups. I will find somewhere to put this fact in the manual. It may
be a useful enhancement request. If you need these discontinuous
charge groups in your topology, you'll have to reorganize it such that
all grouped atoms are consecutive.
-Justin
--
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
Re: [gmx-users] charge groups
Also. Obviously I arranged the charge groups like this so that they
would be neutral, I just didn't realise that the atoms comprising the
charge groups had to be in sequence.
Justin A. Lemkul wrote:
Gavin Melaugh wrote:
As you can see atomnr 7 HC should belong to the same charge group (1) as
atomnr 12. From the tpr file atomnr 7 (in this file 6) seems to have
its own charge group (3).
[atoms]
; atomnr type resnr residue namecgnr charge mass
1 CA 1 CGECA 1 -0.1150 12.0110
2 CB 1 CGECB 1 0. 12.0110
3 CA 1 CGECA 2 -0.1150 12.0110
4 CB 1 CGECB 2 0. 12.0110
5 CA 1 CGECA 3 -0.1150 12.0110
6 CB 1 CGECB 3 0. 12.0110
7 HC 1 CGEHC 1 0.1150 1.0080
8 CU 1 CGECU 13 0.2650 13.0190
9 NU 1 CGENU 13 -0.5970 14.0070
10 HC 1 CGEHC 2 0.1150 1.0080
11 CU 1 CGECU 14 0.2650 13.0190
12 NU 1 CGENU 14 -0.5970 14.0070
13 HC 1 CGEHC 3 0.1150 1.0080
14 CU 1 CGECU 15 0.2650 13.0190
15 NU 1 CGENU 15 -0.5970 14.0070
cgs:
nr=112
cgs[0]={0..1}
cgs[1]={2..3}
cgs[2]={4..5}
cgs[3]={6..6}
cgs[4]={7..8}
cgs[5]={9..9}
cgs[6]={10..11}
cgs[7]={12..12}
cgs[8]={13..14}
cgs[9]={15..16}
cgs[10]={17..18}
cgs[11]={19..20}
cgs[12]={21..21}
cgs[13]={22..23}
cgs[14]={24..24}
cgs[15]={25..26}
cgs[16]={27..27}
cgs[17]={28..29}
cgs[18]={30..31}
cgs[19]={32..33}
cgs[20]={34..35}
cgs[21]={36..36}
cgs[22]={37..38}
cgs[23]={39..39}
cgs[24]={40..41}
cgs[25]={42..42}
cgs[26]={43..44}
cgs[27]={45..46}
Then it seems pretty clear that one cannot have discontinuous charge
groups. I will find somewhere to put this fact in the manual. It may
be a useful enhancement request. If you need these discontinuous
charge groups in your topology, you'll have to reorganize it such that
all grouped atoms are consecutive.
-Justin
--
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
Re: [gmx-users] virtual sites
Hi Sikandar A couple of questions regarding the virtual sites. 1) Do I have to number the virtual site in accordance with the atom indices of the rest of the molecule? 2) Is the parameters for the virtual site declared in the atomtypes directive? Cheers Gavin Sikandar Mashayak wrote: in doing so .. by default all pair interactions with virtual sites would result in zero forces except those between atoms defined in [nonbond_params] On Thu, May 5, 2011 at 11:57 AM, Sikandar Mashayak symasha...@gmail.com mailto:symasha...@gmail.com wrote: hey another approach to do this without using energy group exclusion is to define non-bonded interactions parameter explicitly between atoms in ffnonbonded.itp file. You can specify sigma and epsilon to be zero in virtual sites atoms definition and specify individual pair interactions parameters using non-bonded interactions like following ... ; virtual site VS1 00 0 D 0 0 VS2 00 0 D 0 0 VS3 00 0 D 0 0 [ nonbond_params ] VS1 C 1 1.0 0.25 VS2 C 1 1.0 0.25 VS2 C 1 1.0 0.25 -- sikandar On Thu, May 5, 2011 at 11:47 AM, Gavin Melaugh gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Hi Justin I do not intend to have charges on the sites. All I want is; when a CH3 group gets close to the site it feels a repulsive force. I have calculated a sigma and epsilon value for this interaction. Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin I am reading the manual at the moment. I want to include some virtual sites in my molecule so that only surrounding CH3s atom type C3 interact with then. All other atoms I don't want to interact with them. Do I create energy groups in the index file called say virtual sites and exclusions, and list all the indices of the atom types that I don't want to interact with the virtual site in one group and all the virtual sites in another. e.g [virtual sites] 17 18 19 20 [virtsite_exclus] 1 2 3 4 5 6 7 8 9 . In a general sense, yes, that's the right approach. Note that if any of these sites is charged and/or you're using PME, then this whole exclusion thing goes out the window, as has been discussed several times in recent days. Using energygrp_excl applies only to short-range nonbonded interactions. If you need complete exclusion, you may have to look into tabulated potentials if this is the case. -Justin -- gmx-users mailing listgmx-users@gromacs.org mailto: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 mailto:gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- 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
Re: [gmx-users] virtual sites set up for topology file
Hi I am trying to alter a topology to include 3 virtual sites and I have a few queries, the answers to which are not obvious form the manual. Do I declare the virtual sites in the atomtypes directive like so ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 VS1 0.00.0V 0.0 0.0 VS2 0.00.0V 0.0 0.0 VS3 0.00.0V 0.0 0.0 Do I give them an index number in the atoms directive e.g. [atoms] ; atomnr type resnr residue namecgnr charge mass 1 CA 1 CGECA 1 -0.1150 12.0110 2 CB 1 CGECB 1 0. 12.0110 3 CA 1 CGECA 2 -0.1150 12.0110 4 CB 1 CGECB 2 0. 12.0110 5 CA 1 CGECA 3 -0.1150 12.0110 .. .. 229 VS1 2 VIRVS1 85 0. 0. 230 VS2 3 VIRVS2 86 0. 0. 231 VS3 4 VIRVS3 87 0. 0. Then if I want to set up a virtual site between the COG of 3 atoms do I do it in the following way [virtual_sitesn] ;site COG of three hydrogens at window ;site i jkfunc 229740 581 230 1025 551 231 2837 521 Cheers Gavin Gavin Melaugh wrote: Hi Sikandar A couple of questions regarding the virtual sites. 1) Do I have to number the virtual site in accordance with the atom indices of the rest of the molecule? 2) Is the parameters for the virtual site declared in the atomtypes directive? Cheers Gavin Sikandar Mashayak wrote: in doing so .. by default all pair interactions with virtual sites would result in zero forces except those between atoms defined in [nonbond_params] On Thu, May 5, 2011 at 11:57 AM, Sikandar Mashayak symasha...@gmail.com mailto:symasha...@gmail.com wrote: hey another approach to do this without using energy group exclusion is to define non-bonded interactions parameter explicitly between atoms in ffnonbonded.itp file. You can specify sigma and epsilon to be zero in virtual sites atoms definition and specify individual pair interactions parameters using non-bonded interactions like following ... ; virtual site VS1 00 0 D 0 0 VS2 00 0 D 0 0 VS3 00 0 D 0 0 [ nonbond_params ] VS1 C 1 1.0 0.25 VS2 C 1 1.0 0.25 VS2 C 1 1.0 0.25 -- sikandar On Thu, May 5, 2011 at 11:47 AM, Gavin Melaugh gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Hi Justin I do not intend to have charges on the sites. All I want is; when a CH3 group gets close to the site it feels a repulsive force. I have calculated a sigma and epsilon value for this interaction. Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin I am reading the manual at the moment. I want to include some virtual sites in my molecule so that only surrounding CH3s atom type C3 interact with then. All other atoms I don't want to interact with them. Do I create energy groups in the index file called say virtual sites and exclusions, and list all the indices of the atom types that I don't want to interact with the virtual site in one group and all the virtual sites in another. e.g [virtual sites] 17 18 19 20 [virtsite_exclus] 1 2 3 4 5 6 7 8 9 . In a general sense, yes, that's the right approach. Note that if any of these sites is charged and/or you're using PME, then this whole exclusion thing goes out the window, as has been discussed several times in recent days. Using energygrp_excl applies only to short-range nonbonded interactions. If you need complete exclusion, you may have to look into tabulated potentials if this is the case. -Justin -- gmx-users mailing listgmx-users@gromacs.org mailto:gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx
Re: [gmx-users] virtual sites set up for topology file
Hi Having tried to run grompp using the data below I keep getting the following error Fatal error: Unknown vsiten type 7 Does anyone know why this might be? Gavin Gavin Melaugh wrote: Hi I am trying to alter a topology to include 3 virtual sites and I have a few queries, the answers to which are not obvious form the manual. Do I declare the virtual sites in the atomtypes directive like so ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 VS1 0.00.0V 0.0 0.0 VS2 0.00.0V 0.0 0.0 VS3 0.00.0V 0.0 0.0 Do I give them an index number in the atoms directive e.g. [atoms] ; atomnr type resnr residue namecgnr charge mass 1 CA 1 CGECA 1 -0.1150 12.0110 2 CB 1 CGECB 1 0. 12.0110 3 CA 1 CGECA 2 -0.1150 12.0110 4 CB 1 CGECB 2 0. 12.0110 5 CA 1 CGECA 3 -0.1150 12.0110 .. .. 229 VS1 2 VIRVS1 85 0. 0. 230 VS2 3 VIRVS2 86 0. 0. 231 VS3 4 VIRVS3 87 0. 0. Then if I want to set up a virtual site between the COG of 3 atoms do I do it in the following way [virtual_sitesn] ;site COG of three hydrogens at window ;site i jkfunc 229740 581 230 1025 551 231 2837 521 Cheers Gavin Gavin Melaugh wrote: Hi Sikandar A couple of questions regarding the virtual sites. 1) Do I have to number the virtual site in accordance with the atom indices of the rest of the molecule? 2) Is the parameters for the virtual site declared in the atomtypes directive? Cheers Gavin Sikandar Mashayak wrote: in doing so .. by default all pair interactions with virtual sites would result in zero forces except those between atoms defined in [nonbond_params] On Thu, May 5, 2011 at 11:57 AM, Sikandar Mashayak symasha...@gmail.com mailto:symasha...@gmail.com wrote: hey another approach to do this without using energy group exclusion is to define non-bonded interactions parameter explicitly between atoms in ffnonbonded.itp file. You can specify sigma and epsilon to be zero in virtual sites atoms definition and specify individual pair interactions parameters using non-bonded interactions like following ... ; virtual site VS1 00 0 D 0 0 VS2 00 0 D 0 0 VS3 00 0 D 0 0 [ nonbond_params ] VS1 C 1 1.0 0.25 VS2 C 1 1.0 0.25 VS2 C 1 1.0 0.25 -- sikandar On Thu, May 5, 2011 at 11:47 AM, Gavin Melaugh gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Hi Justin I do not intend to have charges on the sites. All I want is; when a CH3 group gets close to the site it feels a repulsive force. I have calculated a sigma and epsilon value for this interaction. Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin I am reading the manual at the moment. I want to include some virtual sites in my molecule so that only surrounding CH3s atom type C3 interact with then. All other atoms I don't want to interact with them. Do I create energy groups in the index file called say virtual sites and exclusions, and list all the indices of the atom types that I don't want to interact with the virtual site in one group and all the virtual sites in another. e.g [virtual sites] 17 18 19 20 [virtsite_exclus] 1 2 3 4 5 6 7 8 9 . In a general sense, yes, that's the right approach. Note that if any of these sites is charged and/or you're using PME, then this whole exclusion thing goes out the window, as has been discussed several times in recent days. Using energygrp_excl applies only to short-range nonbonded interactions. If you need complete exclusion, you may have to look into tabulated potentials
Re: [gmx-users] virtual sites set up for topology file
Hi Justin Thanks for the reply. To create a virtual site at the centre of geometry of 3 atoms, according to the manual, do you not say: [virtual_sitesn], the index of the site, the index of the three atoms and then the function type 1 which determines that it is COG. Or as I now realise. State [virtual_site3], site index atom indices, and function 1. The fact that there are no parameters then by default must mean it is COG. Is this correct? Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Having tried to run grompp using the data below I keep getting the following error Fatal error: Unknown vsiten type 7 Does anyone know why this might be? Your [virtual_sites] directive is not correct, either in its name (no such thing as virtual_sitesn - the n should be replaced by a digit indicating the type) or in the contents. See manual section 5.2.2. -Justin Gavin Gavin Melaugh wrote: Hi I am trying to alter a topology to include 3 virtual sites and I have a few queries, the answers to which are not obvious form the manual. Do I declare the virtual sites in the atomtypes directive like so ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 VS1 0.00.0V 0.0 0.0 VS2 0.00.0V 0.0 0.0 VS3 0.00.0V 0.0 0.0 Do I give them an index number in the atoms directive e.g. [atoms] ; atomnr type resnr residue namecgnr charge mass 1 CA 1 CGECA 1 -0.1150 12.0110 2 CB 1 CGECB 1 0. 12.0110 3 CA 1 CGECA 2 -0.1150 12.0110 4 CB 1 CGECB 2 0. 12.0110 5 CA 1 CGECA 3 -0.1150 12.0110 .. .. 229 VS1 2 VIRVS1 85 0. 0. 230 VS2 3 VIRVS2 86 0. 0. 231 VS3 4 VIRVS3 87 0. 0. Then if I want to set up a virtual site between the COG of 3 atoms do I do it in the following way [virtual_sitesn] ;site COG of three hydrogens at window ;site i jkfunc 229740 581 230 1025 551 231 2837 521 Cheers Gavin Gavin Melaugh wrote: Hi Sikandar A couple of questions regarding the virtual sites. 1) Do I have to number the virtual site in accordance with the atom indices of the rest of the molecule? 2) Is the parameters for the virtual site declared in the atomtypes directive? Cheers Gavin Sikandar Mashayak wrote: in doing so .. by default all pair interactions with virtual sites would result in zero forces except those between atoms defined in [nonbond_params] On Thu, May 5, 2011 at 11:57 AM, Sikandar Mashayak symasha...@gmail.com mailto:symasha...@gmail.com wrote: hey another approach to do this without using energy group exclusion is to define non-bonded interactions parameter explicitly between atoms in ffnonbonded.itp file. You can specify sigma and epsilon to be zero in virtual sites atoms definition and specify individual pair interactions parameters using non-bonded interactions like following ... ; virtual site VS1 00 0 D 0 0 VS2 00 0 D 0 0 VS3 00 0 D 0 0 [ nonbond_params ] VS1 C 1 1.0 0.25 VS2 C 1 1.0 0.25 VS2 C 1 1.0 0.25 -- sikandar On Thu, May 5, 2011 at 11:47 AM, Gavin Melaugh gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Hi Justin I do not intend to have charges on the sites. All I want is; when a CH3 group gets close to the site it feels a repulsive force. I have calculated a sigma and epsilon value for this interaction. Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin I am reading the manual at the moment. I want to include some virtual sites in my molecule so that only surrounding CH3s atom type C3 interact with then. All other atoms I don't want to interact with them. Do I create energy groups in the index file called say virtual sites and exclusions, and list all the indices of the atom types that I don't want to interact
Re: [gmx-users] virtual sites set up for topology file
Justin I have tried this but I am now getting different errors. I take it that: I specify the virtual sites in the atomtypes directive as I have seen from examples? I index the virtual sites in the atoms directive in accordance with the rest of the molecule. atom numbers go from 1-228, therefore I label the 3 virtual sites 229 to231. The error I get now is Atom index (229) in virtual_sites3 out of bounds (1-228). This probably means that you have inserted topology section virtual_sites3 in a part belonging to a different molecule than you intended to. In that case move the virtual_sites3 section to the right molecule. Do I have to have the virtual sites in the gro file also? This doesn't make sense Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin Thanks for the reply. To create a virtual site at the centre of geometry of 3 atoms, according to the manual, do you not say: [virtual_sitesn], the index of the site, the index of the three atoms and then the function type 1 which determines that it is COG. OK, so that would seem to be right from Table 5.6, but it's not discussed anywhere else, so I suspect that it may be a feature that either got broken along the way, or for some other reason doesn't work, since it's not :) Or as I now realise. State [virtual_site3], site index atom indices, and function 1. The fact that there are no parameters then by default must mean it is COG. Is this correct? I don't know if [virtual_sites3] can be specified without any parameters. It seems to me that this shouldn't be correct. It never hurts to try. -Justin -- 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
Re: [gmx-users] virtual sites set up for topology file
Yeah I see your point about the types. With regard to the initial configuration state I would have assumed that gromacs knew the initial position of the virtual site when I stated that it was to be at the COG of the 3 atoms in the [virtual_sitesn] directive. Justin A. Lemkul wrote: Gavin Melaugh wrote: Justin I have tried this but I am now getting different errors. I take it that: I specify the virtual sites in the atomtypes directive as I have seen from examples? Virtual sites are included in all the force fields already, but if you want some custom name, then yes, include them in a new [atomtypes] directive. I see no reason to create three distinct, but identical, types as you have. I index the virtual sites in the atoms directive in accordance with the rest of the molecule. atom numbers go from 1-228, therefore I label the 3 virtual sites 229 to231. The error I get now is Atom index (229) in virtual_sites3 out of bounds (1-228). This probably means that you have inserted topology section virtual_sites3 in a part belonging to a different molecule than you intended to. In that case move the virtual_sites3 section to the right molecule. Do I have to have the virtual sites in the gro file also? This doesn't make sense Yes, you need their coordinates as part of the initial state. -Justin -- 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
Re: [gmx-users] virtual sites set up for topology file
O.k Cheers Justin A. Lemkul wrote: Gavin Melaugh wrote: Yeah I see your point about the types. With regard to the initial configuration state I would have assumed that gromacs knew the initial position of the virtual site when I stated that it was to be at the COG of the 3 atoms in the [virtual_sitesn] directive. I believe this information is only used during mdrun to construct the position after forces have been applied and coordinates updated. You can run a simple protein through pdb2gmx with -vsite hydrogens to see how all of this should be put together. -Justin -- 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
Re: [gmx-users] virtual sites set up for topology file
Hi Mark Cheers.Could you perhaps shed some insight into format of [virtual_siten] when trying to define the VS at the COG of three atoms. It is not obvious from the manual Cheers Gavin Mark Abraham wrote: On 6/05/2011 10:17 PM, Gavin Melaugh wrote: Yeah I see your point about the types. With regard to the initial configuration state I would have assumed that gromacs knew the initial position of the virtual site when I stated that it was to be at the COG of the 3 atoms in the [virtual_sitesn] directive. Yes it does, but (IIRC) the implementation requires that each of the constructing atoms and the virtual atom(s) are distinct atoms. mdrun does the book-keeping to know which atoms are sensible at different stages of the integration. So the input coordinate file needs to have such atoms there, but I rather suspect their coordinates will get ignored - that's easily tested, of course. Mark Justin A. Lemkul wrote: Gavin Melaugh wrote: Justin I have tried this but I am now getting different errors. I take it that: I specify the virtual sites in the atomtypes directive as I have seen from examples? Virtual sites are included in all the force fields already, but if you want some custom name, then yes, include them in a new [atomtypes] directive. I see no reason to create three distinct, but identical, types as you have. I index the virtual sites in the atoms directive in accordance with the rest of the molecule. atom numbers go from 1-228, therefore I label the 3 virtual sites 229 to231. The error I get now is Atom index (229) in virtual_sites3 out of bounds (1-228). This probably means that you have inserted topology section virtual_sites3 in a part belonging to a different molecule than you intended to. In that case move the virtual_sites3 section to the right molecule. Do I have to have the virtual sites in the gro file also? This doesn't make sense Yes, you need their coordinates as part of the initial state. -Justin -- 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
Re: [gmx-users] virtual sites set up for topology file
Hi Mark Many thanks to you and all your colleagues for replying. This has worked, at least there are now no errors. Where the manual is incorrect is that it leads you to believe that you state; [virtual _siten] ;index of VS index of atoms for COG func 229 8 11 151 However the function type should be stated before the atom indices comprising the COG ;index of VS func indices 229 1 8 11 15 Cheers Gavin Mark Abraham wrote: On 6/05/2011 10:39 PM, Gavin Melaugh wrote: Hi Mark Cheers.Could you perhaps shed some insight into format of [virtual_siten] when trying to define the VS at the COG of three atoms. It is not obvious from the manual Indeed, it's undocumented - but I think Sikander's experience from my discussion in this thread earlier this month http://lists.gromacs.org/pipermail/gmx-users/2011-May/060994.html should point the way for you. Mark Cheers Gavin Mark Abraham wrote: On 6/05/2011 10:17 PM, Gavin Melaugh wrote: Yeah I see your point about the types. With regard to the initial configuration state I would have assumed that gromacs knew the initial position of the virtual site when I stated that it was to be at the COG of the 3 atoms in the [virtual_sitesn] directive. Yes it does, but (IIRC) the implementation requires that each of the constructing atoms and the virtual atom(s) are distinct atoms. mdrun does the book-keeping to know which atoms are sensible at different stages of the integration. So the input coordinate file needs to have such atoms there, but I rather suspect their coordinates will get ignored - that's easily tested, of course. Mark Justin A. Lemkul wrote: Gavin Melaugh wrote: Justin I have tried this but I am now getting different errors. I take it that: I specify the virtual sites in the atomtypes directive as I have seen from examples? Virtual sites are included in all the force fields already, but if you want some custom name, then yes, include them in a new [atomtypes] directive. I see no reason to create three distinct, but identical, types as you have. I index the virtual sites in the atoms directive in accordance with the rest of the molecule. atom numbers go from 1-228, therefore I label the 3 virtual sites 229 to231. The error I get now is Atom index (229) in virtual_sites3 out of bounds (1-228). This probably means that you have inserted topology section virtual_sites3 in a part belonging to a different molecule than you intended to. In that case move the virtual_sites3 section to the right molecule. Do I have to have the virtual sites in the gro file also? This doesn't make sense Yes, you need their coordinates as part of the initial state. -Justin -- 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
Re: [gmx-users] virtual sites set up for topology file
Am I correct in saying now, from the following topology (exerpts), that the virtual site VS will only interact with C3? I guess I don't have to give the atom indices of this interaction in the pair list which I use only for 1_4 interactions? Can I use sigma and epsilon in the nonbond_params directive like in atomtypes. ;Parameter level [defaults] ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ 1 3 yes0.5 0.5 [atomtypes] ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 HC 1.008000 0.115000 A 0.242000 0.125520 CU 13.019000 0.265000 A 0.35 0.334720 NU 14.007000 -0.597000 A 0.325000 0.711280 CH 13.019000 0.332000 A 0.385000 0.334720 C3 15.035000 0.00 A 0.390500 0.732200 C2 14.027000 0.00 A 0.390500 0.493712 VS 0.0 0.0V 0.0 0.0 [nonbond_params] ;ij func sigmaepsilon VS C31 0.1 0.03153 -- 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
Re: [gmx-users] virtual sites set up for topology file
Cheers Sikandar I take it that because combination rule 3 (provide sigma and epsilon) is stated gromacs assumes that all values in nonbonding parameters are sigma and epsilon. I know this tp be tru for the atomtypes but does it filter down to all intermolecular interactions. Cheers Gavin Sikandar Mashayak wrote: Yes, since sigma and epsilon are zero for VS then interactions between and VS-VS and VS-(any other atom) would result in zero force. Since you explicitly define the interaction between VS-C3, combination rule won't be used and C6(VS_C3) and C12(VS_C3) will be computed as per sigma(VS_C3) epsilon(VS_C3) you defined under [nonbond_params] cheers sikandar On Fri, May 6, 2011 at 9:03 AM, Gavin Melaugh gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Am I correct in saying now, from the following topology (exerpts), that the virtual site VS will only interact with C3? I guess I don't have to give the atom indices of this interaction in the pair list which I use only for 1_4 interactions? Can I use sigma and epsilon in the nonbond_params directive like in atomtypes. ;Parameter level [defaults] ; nbfunccomb-rule gen-pairsfudgeLJ fudgeQQ 1 3 yes0.5 0.5 [atomtypes] ;type mass charge ptype sigma(nm) epsilon(kjmol-1) CB 12.011000 0.00 A 0.355000 0.292880 CA 12.011000 -0.115000 A 0.355000 0.292880 HC 1.008000 0.115000 A 0.242000 0.125520 CU 13.019000 0.265000 A 0.35 0.334720 NU 14.007000 -0.597000 A 0.325000 0.711280 CH 13.019000 0.332000 A 0.385000 0.334720 C3 15.035000 0.00 A 0.390500 0.732200 C2 14.027000 0.00 A 0.390500 0.493712 VS 0.0 0.0V 0.0 0.0 [nonbond_params] ;ij func sigmaepsilon VS C31 0.1 0.03153 -- gmx-users mailing listgmx-users@gromacs.org mailto: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 mailto:gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- 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] nrexcl=3
Hi all A very quick question. From the manual the definition of nrexcl =3 is excluding nonbonded interactions between atoms that are no further than 3 bonds away. In say butane does this means that atoms 1 and 4 do not interact, or I take it they interact according to the pair list. Gavin -- 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
Hi All Is it possible to have a virtual site interact with only specific atoms and not interact at all with everything else? Cheers Gavin -- 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
Re: [gmx-users] virtual sites
Hi Sikandar Cheers. How do you actually define the energy groups? Gavin Sikandar Mashayak wrote: yes, you can make virtual sites interact with only specific sites by using Energy Exclusion between energy groups. This can be done by defining energy groups for virtual sites and other atoms, then exclude or include the non-bonded interactions between them accordingly... On Thu, May 5, 2011 at 7:30 AM, Gavin Melaugh gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Hi All Is it possible to have a virtual site interact with only specific atoms and not interact at all with everything else? Cheers Gavin -- gmx-users mailing listgmx-users@gromacs.org mailto: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 mailto:gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- 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
Re: [gmx-users] virtual sites
Hi Justin I am reading the manual at the moment. I want to include some virtual sites in my molecule so that only surrounding CH3s atom type C3 interact with then. All other atoms I don't want to interact with them. Do I create energy groups in the index file called say virtual sites and exclusions, and list all the indices of the atom types that I don't want to interact with the virtual site in one group and all the virtual sites in another. e.g [virtual sites] 17 18 19 20 [virtsite_exclus] 1 2 3 4 5 6 7 8 9 . Cheers Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Sikandar Cheers. How do you actually define the energy groups? Read in the manual about the .mdp option energygrps and apply custom index groups as necessary. -Justin Gavin Sikandar Mashayak wrote: yes, you can make virtual sites interact with only specific sites by using Energy Exclusion between energy groups. This can be done by defining energy groups for virtual sites and other atoms, then exclude or include the non-bonded interactions between them accordingly... On Thu, May 5, 2011 at 7:30 AM, Gavin Melaugh gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Hi All Is it possible to have a virtual site interact with only specific atoms and not interact at all with everything else? Cheers Gavin -- gmx-users mailing listgmx-users@gromacs.org mailto: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 mailto:gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- 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
Re: [gmx-users] virtual sites
Hi Sikandar Many Thanks for the reply. When you define energygrps in the mdp file, does that mean that only these groups are written to the energy file. I take it that all iother information regarding the energy is written as well i.e if you left this section blank Cheers Gavin Sikandar Mashayak wrote: you can do it as following 1. create index file make_indx -f conf.gro -n index.ndx ( select VS and any other atoms you want lets say OTHERS) 2. in .mdp define energygrps VS OTHERS 3. exclude interactions by specifying energygrp_excl VS OTHERS in .mdp file for more you can refer to manual.. I hope it helps -- sikandar On Thu, May 5, 2011 at 11:00 AM, Gavin Melaugh gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Hi Sikandar Cheers. How do you actually define the energy groups? Gavin Sikandar Mashayak wrote: yes, you can make virtual sites interact with only specific sites by using Energy Exclusion between energy groups. This can be done by defining energy groups for virtual sites and other atoms, then exclude or include the non-bonded interactions between them accordingly... On Thu, May 5, 2011 at 7:30 AM, Gavin Melaugh gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk mailto:gmelaug...@qub.ac.uk wrote: Hi All Is it possible to have a virtual site interact with only specific atoms and not interact at all with everything else? Cheers Gavin -- gmx-users mailing listgmx-users@gromacs.org mailto:gmx-users@gromacs.org mailto:gmx-users@gromacs.org mailto: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 mailto:gmx-users-requ...@gromacs.org mailto:gmx-users-requ...@gromacs.org mailto:gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org mailto: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 mailto:gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- 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
Re: [gmx-users] virtual sites
Hi Justin I do not intend to have charges on the sites. All I want is; when a CH3 group gets close to the site it feels a repulsive force. I have calculated a sigma and epsilon value for this interaction. Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin I am reading the manual at the moment. I want to include some virtual sites in my molecule so that only surrounding CH3s atom type C3 interact with then. All other atoms I don't want to interact with them. Do I create energy groups in the index file called say virtual sites and exclusions, and list all the indices of the atom types that I don't want to interact with the virtual site in one group and all the virtual sites in another. e.g [virtual sites] 17 18 19 20 [virtsite_exclus] 1 2 3 4 5 6 7 8 9 . In a general sense, yes, that's the right approach. Note that if any of these sites is charged and/or you're using PME, then this whole exclusion thing goes out the window, as has been discussed several times in recent days. Using energygrp_excl applies only to short-range nonbonded interactions. If you need complete exclusion, you may have to look into tabulated potentials if this is the case. -Justin -- 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] free energy
Hi all I have performed a PMF simulation of taking part of a molecule out of the cavity of a host using umbrella sampling. The free energy curve suggests that the guest prefers to be outside the host as the bound state is higher in energy, or the free energy difference to go in is positive. However when I view the trajectories for each window it appears that there is always more bound states than unbound. This leads me to doubt my free energy profile? Cheers Gavin -- 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
Re: [gmx-users] free energy
Hi Chris My windows are restrained obviously using the force constant in the mdp file. The trajectories that I have viewed are those of the individual biased sampling windows. I have not put on the unbiased simulations yet. There is also the following issue: The simulations involve two identical molecules containing hydrocarbon chains. I calculate the PMF to take a specific hydrocarbon chain of one molecule out of a specific site on the neighbouring molecule. When this guest chain goes out, other chains can go in (intramolecular or intermolecular). Will this affect the profile? Gavin chris.ne...@utoronto.ca wrote: Your windows are restrained. The PMF that you get out of WHAM is a representation of the relative sampling after removing the umbrella biases. Sounds like you are saying that you look at the still-biased trajectories and you see different a different distribution of states than you do in the de-biased PMF. not sure what the problem is here. Perhaps go back to some review literature on US. Now, if you saw more bound than unbound in unrestrained simulations, then that's a different story, but that doesn't appear to be what you are talking about. Chris. -- original message -- Hi all I have performed a PMF simulation of taking part of a molecule out of the cavity of a host using umbrella sampling. The free energy curve suggests that the guest prefers to be outside the host as the bound state is higher in energy, or the free energy difference to go in is positive. However when I view the trajectories for each window it appears that there is always more bound states than unbound. This leads me to doubt my free energy profile? Cheers Gavin -- 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
Re: [gmx-users] free energy
The current simulations are currently in vacuum. Does the following mdp file seem ok. Note that this is a production run using the final configuration after a lengthy equilibration. Also I was thinking about trying to prevent the other chains from entering the specific site; is there a a way to implement this in gromacs. Gavin chris.ne...@utoronto.ca wrote: yes it will, and so will it affect the profile if water molecules or ions go in when both chains are absent. You'll need to determine what question you are trying to answer and also think pretty hard about what your PMF really means in the context of this system. Chris. -- original message -- Hi Chris My windows are restrained obviously using the force constant in the mdp file. The trajectories that I have viewed are those of the individual biased sampling windows. I have not put on the unbiased simulations yet. There is also the following issue: The simulations involve two identical molecules containing hydrocarbon chains. I calculate the PMF to take a specific hydrocarbon chain of one molecule out of a specific site on the neighbouring molecule. When this guest chain goes out, other chains can go in (intramolecular or intermolecular). Will this affect the profile? Gavin umbrella81.mdp Description: application/mdp -- 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
Re: [gmx-users] Umbrella Sampling
Hi I assume that the order of the file names in the tpr-files.dat and pullx-files.dat is irrelevant for g_wham . Cheers Gavin -- 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
Re: [gmx-users] Umbrella Sampling
Hi Justin Yeah I know the tpr files must be in corresponding order with pullx.xvg files. What I meant was should they be in order of distance. i.e say If my windows go from 0 to 1.0 nm with windows every 0.1nm, could I list the files in any order or does it have to like 0.1 0.2 0.3 0.4 Gavin Justin A. Lemkul wrote: From g_wham -h: The tpr and pullx files must be in corresponding order, i.e. the first tpr created the first pullx, etc. -Justin Gavin Melaugh wrote: Hi I assume that the order of the file names in the tpr-files.dat and pullx-files.dat is irrelevant for g_wham . Cheers Gavin -- 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
Re: [gmx-users] Umbrella Sampling
Cheers Chris Is the best way to check for convergence; to keep adding in more histograms until the curves converge. Also your comment 'don't remove any data', do you mean to keep histograms that are not so good. Gavin chris.ne...@utoronto.ca wrote: your comment: which should be centred around 0.80nm is flawed, as i mentioned earlier. also, it is not g_wham that is sensitive but the convergence and sampling of phase space that is sensitive. don`t remove any data. do evaluate your convergence. without convergence measures, a pmf is worse than useless. chris. -- original message -- Cheers Chris If I remove the red histogram (the first of the wider distributions), which should be centred around 0.80nm but is actually centred around 0.78 nm; and add in some more histograms with higher force constants the profile changes slightly. It seems that g_wham is very sensitive to these subtleties. How do I know which curve is correct? I have about six such curves that differ slightly in this manner. Gavin chris.ne...@utoronto.ca wrote: [Hide Quoted Text] looks fine to me, no need to do that extra sampling that I suggested since it appears that you already did this -- benefits of seeing real data ;). If you want to understand why your histograms are not always centered at r0 (note that this is just fine) then you should read more about US, WHAM, and how to bias/debias the data for US (I am sure that there are textbooks around that explain this). The only case in which all of your histograms will be centered at their respective r0 is when the underlying PMF is exactly flat. Chris. -- original message -- Hi Chris many thanks again for the advise. I have, or at least I thought have sampled my barrier region to death, but as I say some histograms may not be centred around r0. I will proceed with what you suggest. Please find attached a picture of the histograms, the corresponding profile, and a sample mdp file that I use. -- 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] Umbrella sampling
Hi All I have generated several PMF curves for the one system using umbrella sampling. In the first part of the curve (barrier region) I use a high force constant with small intervals between the windows. The latter part of the curve I use a lower force constant with larger window spacing. Anyway I have a few issues that I need clarifying: 1 - Can you have too much overlap between windows? 2 - Does the distribution at each window have to centered around the desired r0? (If not does this affect the free energy?) 3- If you over sample one particular window, will it affect the curve? Many thanks GAvin -- 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
Re: [gmx-users] Umbrella sampling
Hi Xavier Thanks for the reply. With respect to your answer of my first query. What if you had two windows practically on top of each other, but one was not supposed to be there. e.g A window with r0 of 0.80 nm and centred at 0.78 nm and a window with r0 of 0.78 nm centred at 0.78nm. Gavin XAvier Periole wrote: On Mar 31, 2011, at 11:53 AM, Gavin Melaugh wrote: Hi All I have generated several PMF curves for the one system using umbrella sampling. In the first part of the curve (barrier region) I use a high force constant with small intervals between the windows. The latter part of the curve I use a lower force constant with larger window spacing. Anyway I have a few issues that I need clarifying: 1 - Can you have too much overlap between windows? no, there no such a thing of too much overlap :)) You could even put two identical windows with same 100% overlap ... no problem. 2 - Does the distribution at each window have to centered around the desired r0? (If not does this affect the free energy?) The deviation of the distribution from the r0 is what dictates the profile. The more away from the disired r0 the higher the free energy of the system. 3- If you over sample one particular window, will it affect the curve? There is no such a thing of over sampling ... the only thing you can have is not enough sampling. Many thanks GAvin -- 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 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
Re: [gmx-users] Umbrella sampling
Sorry I am not sure that I follow. Will the window with r0 =0.80 giving the distribution centred around 0.78nm not drive my free energy profile up. If I remove this window prior to running g_wham the free energy goes down. Should I increase the force constant so that the mean of the window is 0.80nm (bearing in mind that this is near the barrier region). Gavin XAvier Periole wrote: You can present the data differently: you have two windows at 0.78 nm giving different distribution. That indicates these windows are not converged. Does not mean that the others (0.80 nm) are converged :)) On Mar 31, 2011, at 12:20 PM, Gavin Melaugh wrote: Hi Xavier Thanks for the reply. With respect to your answer of my first query. What if you had two windows practically on top of each other, but one was not supposed to be there. e.g A window with r0 of 0.80 nm and centred at 0.78 nm and a window with r0 of 0.78 nm centred at 0.78nm. Gavin XAvier Periole wrote: On Mar 31, 2011, at 11:53 AM, Gavin Melaugh wrote: Hi All I have generated several PMF curves for the one system using umbrella sampling. In the first part of the curve (barrier region) I use a high force constant with small intervals between the windows. The latter part of the curve I use a lower force constant with larger window spacing. Anyway I have a few issues that I need clarifying: 1 - Can you have too much overlap between windows? no, there no such a thing of too much overlap :)) You could even put two identical windows with same 100% overlap ... no problem. 2 - Does the distribution at each window have to centered around the desired r0? (If not does this affect the free energy?) The deviation of the distribution from the r0 is what dictates the profile. The more away from the disired r0 the higher the free energy of the system. 3- If you over sample one particular window, will it affect the curve? There is no such a thing of over sampling ... the only thing you can have is not enough sampling. Many thanks GAvin -- 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 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 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
Re: [gmx-users] Umbrella sampling
Thanks Justin for the reply Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Sorry I am not sure that I follow. Will the window with r0 =0.80 giving the distribution centred around 0.78nm not drive my free energy profile up. If I remove this window prior to running g_wham the free energy goes down. Should I increase the force constant so that the mean of the window is 0.80nm (bearing in mind that this is near the barrier region). If you have an incomplete or otherwise discontinuous free energy profile, then you won't get a correct result, but it's not simply due to oversampling one region. It's that the oversampling results in undersampling another region. Increasing the force constant for the window centered around 0.80 nm should work. -Justin Gavin XAvier Periole wrote: You can present the data differently: you have two windows at 0.78 nm giving different distribution. That indicates these windows are not converged. Does not mean that the others (0.80 nm) are converged :)) On Mar 31, 2011, at 12:20 PM, Gavin Melaugh wrote: Hi Xavier Thanks for the reply. With respect to your answer of my first query. What if you had two windows practically on top of each other, but one was not supposed to be there. e.g A window with r0 of 0.80 nm and centred at 0.78 nm and a window with r0 of 0.78 nm centred at 0.78nm. Gavin XAvier Periole wrote: On Mar 31, 2011, at 11:53 AM, Gavin Melaugh wrote: Hi All I have generated several PMF curves for the one system using umbrella sampling. In the first part of the curve (barrier region) I use a high force constant with small intervals between the windows. The latter part of the curve I use a lower force constant with larger window spacing. Anyway I have a few issues that I need clarifying: 1 - Can you have too much overlap between windows? no, there no such a thing of too much overlap :)) You could even put two identical windows with same 100% overlap ... no problem. 2 - Does the distribution at each window have to centered around the desired r0? (If not does this affect the free energy?) The deviation of the distribution from the r0 is what dictates the profile. The more away from the disired r0 the higher the free energy of the system. 3- If you over sample one particular window, will it affect the curve? There is no such a thing of over sampling ... the only thing you can have is not enough sampling. Many thanks GAvin -- 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 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 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] Umbrella Sampling
Hi Chris many thanks again for the advise. I have, or at least I thought have sampled my barrier region to death, but as I say some histograms may not be centred around r0. I will proceed with what you suggest. Please find attached a picture of the histograms, the corresponding profile, and a sample mdp file that I use. groprofile.agr.gz Description: GNU Zip compressed data hist.agr.gz Description: GNU Zip compressed data -- 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
Re: [gmx-users] Umbrella sampling
Here is also the sample .mdp file chris.ne...@utoronto.ca wrote: can you show us your mdp and the pmf and the histograms for the data that you put into wham? It's a lot easier to diagnose with the real data. In the general case where umbrellas are spaced equally along your reaction coordinate, sampling overlap between umbrellas will always decrease anywhere the PMF is concave down, such as your barrier region. I would suggest adding windows at every 0.005 nm spanning the region that you consider to have a sampling problem (e.g. 0.76 to 0.83 ?) and use a much stronger force constant here (e.g. 2-3 times as strong as you used for umbrellas with 0.02 nm spacing). You have identified a problem, so my suggestion is to not bother fiddling with adding one extra replica but sample that region to death with strong force constants. I presume that this will not impact your overall efficiency too much if the reaction coordinate is long overall. Finally, you will need to evaluate the convergence of your PMF overall and perhaps of this region in particular, especially if you want to know the dG to cross it or between two low energy states on either side of it. Chris. -- original message -- Sorry I am not sure that I follow. Will the window with r0 =0.80 giving the distribution centred around 0.78nm not drive my free energy profile up. If I remove this window prior to running g_wham the free energy goes down. Should I increase the force constant so that the mean of the window is 0.80nm (bearing in mind that this is near the barrier region). Gavin XAvier Periole wrote: You can present the data differently: you have two windows at 0.78 nm giving different distribution. That indicates these windows are not converged. Does not mean that the others (0.80 nm) are converged :)) On Mar 31, 2011, at 12:20 PM, Gavin Melaugh wrote: Hi Xavier Thanks for the reply. With respect to your answer of my first query. What if you had two windows practically on top of each other, but one was not supposed to be there. e.g A window with r0 of 0.80 nm and centred at 0.78 nm and a window with r0 of 0.78 nm centred at 0.78nm. Gavin XAvier Periole wrote: On Mar 31, 2011, at 11:53 AM, Gavin Melaugh wrote: Hi All I have generated several PMF curves for the one system using umbrella sampling. In the first part of the curve (barrier region) I use a high force constant with small intervals between the windows. The latter part of the curve I use a lower force constant with larger window spacing. Anyway I have a few issues that I need clarifying: 1 - Can you have too much overlap between windows? no, there no such a thing of too much overlap :)) You could even put two identical windows with same 100% overlap ... no problem. 2 - Does the distribution at each window have to centered around the desired r0? (If not does this affect the free energy?) The deviation of the distribution from the r0 is what dictates the profile. The more away from the disired r0 the higher the free energy of the system. 3- If you over sample one particular window, will it affect the curve? There is no such a thing of over sampling ... the only thing you can have is not enough sampling. Many thanks GAvin -- gmx-users mailing listgmx-users at 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-request at gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users at 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-request at gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists umbrella114.mdp Description: application/mdp -- 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
Re: [gmx-users] Umbrella Sampling
Cheers Chris If I remove the red histogram (the first of the wider distributions), which should be centred around 0.80nm but is actually centred around 0.78 nm; and add in some more histograms with higher force constants the profile changes slightly. It seems that g_wham is very sensitive to these subtleties. How do I know which curve is correct? I have about six such curves that differ slightly in this manner. Gavin chris.ne...@utoronto.ca wrote: looks fine to me, no need to do that extra sampling that I suggested since it appears that you already did this -- benefits of seeing real data ;). If you want to understand why your histograms are not always centered at r0 (note that this is just fine) then you should read more about US, WHAM, and how to bias/debias the data for US (I am sure that there are textbooks around that explain this). The only case in which all of your histograms will be centered at their respective r0 is when the underlying PMF is exactly flat. Chris. -- original message -- Hi Chris many thanks again for the advise. I have, or at least I thought have sampled my barrier region to death, but as I say some histograms may not be centred around r0. I will proceed with what you suggest. Please find attached a picture of the histograms, the corresponding profile, and a sample mdp file that I use. -- next part -- A non-text attachment was scrubbed... Name: groprofile.agr.gz Type: application/x-gzip Size: 3805 bytes Desc: not available Url : http://lists.gromacs.org/pipermail/gmx-users/attachments/20110331/2a5b0bd6/groprofile.agr.gz -- next part -- A non-text attachment was scrubbed... Name: hist.agr.gz Type: application/x-gzip Size: 25250 bytes Desc: not available Url : http://lists.gromacs.org/pipermail/gmx-users/attachments/20110331/2a5b0bd6/hist.agr.gz -- modifiedhist.xvg.gz Description: GNU Zip compressed data modified.agr.gz Description: GNU Zip compressed data -- 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.xvg format
Hi all I have a very quick question. The output for umbrella pulling simulations (pullx.xvg) often appears to have an irregular format, in that the columns don't seem to be perfectly aligned (see below). I have never had a problem with this as g_wham never reported any errors. I was just wondering is this an issue? Cheers Gavin 8631.6006 1.05024 3.941 3.80205 0.0919597 0.147718 -0.0215892 8631.9004 1.03157 3.94253 3.80054 0.0290956 0.175955 -0.0447413 8632.2002 1.01474 3.96139 3.81348 0.0594728 0.160099 -0.0482623 8632.5000 1.0197 3.95987 3.82872 0.0515857 0.140793 -0.0359526 8632.8008 1.03466 3.94717 3.83522 0.0775257 0.156349 -0.0268523 8633.1006 1.03959 3.94529 3.83996 0.0630389 0.140153 -0.0315862 8633.4004 1.0372 3.95998 3.82624 0.0489344 0.1468 -0.0239522 8633.7002 1.02653 3.95994 3.82712 0.0632881 0.142731 -0.0574301 8634. 1.01196 3.95079 3.81532 0.0594033 0.128616 -0.057459 8634.3008 1.00748 3.93597 3.82667 0.1461230.114964 -0.0386781 8634.6006 1.00209 3.91869 3.83553 0.0928916 0.119439 -0.0338013 8634.9004 0.9937043.92392 3.8084 0.0669909 0.152878-0.0255926 8635.2002 1.00164 3.94898 3.79449 0.0824905 0.170253 -0.0137497 8635.5000 1.0168 3.95829 3.78047 0.0581190.158457 -0.0219142 8635.8008 1.03037 3.95717 3.7914 0.0950175 0.145733 -0.0685743 8636.1006 1.0152 3.94208 3.78807 0.0705052 0.147887 -0.0117258 8636.4004 1.00647 3.93244 3.79199 0.0921640.130635 -0.00548112 8636.7002 1.02261 3.91966 3.76668 0.0743383 0.160227 -0.0491515 8637. 1.0466 3.91742 3.73989 0.0766793 0.134678 -0.09385 8637.3008 1.04931 3.93179 3.75026 0.1267890.115951 -0.0487055 8637.6006 1.04238 3.94086 3.73288 0.1274040.150352 0.00932272 8637.9004 1.02266 3.92861 3.74975 0.0541594 0.164885 -0.0608739 8638.2002 1.01286 3.93057 3.77248 0.14425 0.10085 -0.00612713 8638.5000 1.01337 3.94557 3.79916 0.11815 0.0991607 0.00891958 8638.8008 0.9870623.93091 3.81182 0.0363212 0.168897-0.0292177 8639.1006 0.9953363.93407 3.82284 0.122306 0.0993002 -0.00813846 8639.4004 1.00739 3.95131 3.82143 0.0971151 0.124606 -0.0229913 8639.7002 1.01473 3.98067 3.82504 0.1302830.117168 -0.0369952 8640. 1.02179 3.97666 3.80322 0.1076740.13277 -0.0346299 8640.3008 1.02445 3.96855 3.79499 0.0598816 0.185437 -0.0311454 8640.6006 1.03611 3.94857 3.79149 0.1051130.122864 0.00836338 8640.9004 1.02776 3.94737 3.79431 0.09896 0.093913-0.0594025 8641.2002 0.9931783.93066 3.80311 0.0987073 0.146735-0.0481467 8641.5000 0.9816533.91943 3.7912 0.0527064 0.155903-0.0340922 8641.8008 0.99508 3.92219 3.80137 0.1122350.08475 -0.0709295 8642.1006 1.0068 3.91384 3.78044 0.0627984 0.156858 -0.0764401 8642.4004 1.00386 3.93476 3.80116 0.0806420.0993683 -0.050131 8642.7002 0.9932483.90789 3.80662 0.0890923 0.148829-0.0430867 8643. 0.9925173.87576 3.79073 0.091716 0.134674-0.0131645 8643.3008 0.9952343.89502 3.7805 0.0704844 0.156837-0.0469381 8643.6006 1.01357 3.91115 3.75715 0.0796926 0.128354 -0.0687171 8643.9004 1.01555 3.90522 3.75396 0.0917188 0.140592 -0.0575233 -- 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
Re: [gmx-users] pullx.xvg format
Hi Justin Thanks for the reply. I expanded the terminal and it didn't make a difference. Do you still think it's O.K? Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi all I have a very quick question. The output for umbrella pulling simulations (pullx.xvg) often appears to have an irregular format, in that the columns don't seem to be perfectly aligned (see below). I have never had a problem with this as g_wham never reported any errors. I was just wondering is this an issue? It shouldn't be. Data are parsed based on line breaks. I imagine this output might look normal if you expand the size of your terminal to get everything on one line. -Justin Cheers Gavin 8631.6006 1.05024 3.941 3.80205 0.0919597 0.147718 -0.0215892 8631.9004 1.03157 3.94253 3.80054 0.0290956 0.175955 -0.0447413 8632.2002 1.01474 3.96139 3.81348 0.0594728 0.160099 -0.0482623 8632.5000 1.0197 3.95987 3.82872 0.0515857 0.140793 -0.0359526 8632.8008 1.03466 3.94717 3.83522 0.0775257 0.156349 -0.0268523 8633.1006 1.03959 3.94529 3.83996 0.0630389 0.140153 -0.0315862 8633.4004 1.0372 3.95998 3.82624 0.0489344 0.1468 -0.0239522 8633.7002 1.02653 3.95994 3.82712 0.0632881 0.142731 -0.0574301 8634. 1.01196 3.95079 3.81532 0.0594033 0.128616 -0.057459 8634.3008 1.00748 3.93597 3.82667 0.146123 0.114964 -0.0386781 8634.6006 1.00209 3.91869 3.83553 0.0928916 0.119439 -0.0338013 8634.9004 0.9937043.92392 3.8084 0.0669909 0.152878-0.0255926 8635.2002 1.00164 3.94898 3.79449 0.0824905 0.170253 -0.0137497 8635.5000 1.0168 3.95829 3.78047 0.058119 0.158457 -0.0219142 8635.8008 1.03037 3.95717 3.7914 0.0950175 0.145733 -0.0685743 8636.1006 1.0152 3.94208 3.78807 0.0705052 0.147887 -0.0117258 8636.4004 1.00647 3.93244 3.79199 0.092164 0.130635 -0.00548112 8636.7002 1.02261 3.91966 3.76668 0.0743383 0.160227 -0.0491515 8637. 1.0466 3.91742 3.73989 0.0766793 0.134678 -0.09385 8637.3008 1.04931 3.93179 3.75026 0.126789 0.115951 -0.0487055 8637.6006 1.04238 3.94086 3.73288 0.127404 0.150352 0.00932272 8637.9004 1.02266 3.92861 3.74975 0.0541594 0.164885 -0.0608739 8638.2002 1.01286 3.93057 3.77248 0.14425 0.10085 -0.00612713 8638.5000 1.01337 3.94557 3.79916 0.11815 0.0991607 0.00891958 8638.8008 0.9870623.93091 3.81182 0.0363212 0.168897-0.0292177 8639.1006 0.9953363.93407 3.82284 0.122306 0.0993002 -0.00813846 8639.4004 1.00739 3.95131 3.82143 0.0971151 0.124606 -0.0229913 8639.7002 1.01473 3.98067 3.82504 0.130283 0.117168 -0.0369952 8640. 1.02179 3.97666 3.80322 0.1076740.13277 -0.0346299 8640.3008 1.02445 3.96855 3.79499 0.0598816 0.185437 -0.0311454 8640.6006 1.03611 3.94857 3.79149 0.105113 0.122864 0.00836338 8640.9004 1.02776 3.94737 3.79431 0.09896 0.093913 -0.0594025 8641.2002 0.9931783.93066 3.80311 0.0987073 0.146735-0.0481467 8641.5000 0.9816533.91943 3.7912 0.0527064 0.155903-0.0340922 8641.8008 0.99508 3.92219 3.80137 0.1122350.08475 -0.0709295 8642.1006 1.0068 3.91384 3.78044 0.0627984 0.156858 -0.0764401 8642.4004 1.00386 3.93476 3.80116 0.080642 0.0993683 -0.050131 8642.7002 0.9932483.90789 3.80662 0.0890923 0.148829-0.0430867 8643. 0.9925173.87576 3.79073 0.091716 0.134674-0.0131645 8643.3008 0.9952343.89502 3.7805 0.0704844 0.156837-0.0469381 8643.6006 1.01357 3.91115 3.75715 0.0796926 0.128354 -0.0687171 8643.9004 1.01555 3.90522 3.75396 0.0917188 0.140592 -0.0575233 -- 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] sum of largest charge group radii
Hi I have installed the new version of gromacs and so far I only use it as a check when running grompp as it gnerates more information about cut off charge group radii etc. For a system of 256 hexane molecules in a cubic box I get the following error. WARNING 1 [file pbc.mdp]: The sum of the two largest charge group radii (5.231524) is larger than rlist (1.70) This doesn't even make sense since my as system itself fits in a 3.2 nm cubic box. I am pretty sure my topology is correct. I have seen on the mailing list that there have been problems with this but they have all been for bigger molecules. Gavin -- 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
Re: [gmx-users] sum of largest charge group radii
Hi Justin By broken do you mean split at the edge of the periodic box? If so then yes. Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi I have installed the new version of gromacs and so far I only use it as a check when running grompp as it gnerates more information about cut off charge group radii etc. For a system of 256 hexane molecules in a cubic box I get the following error. WARNING 1 [file pbc.mdp]: The sum of the two largest charge group radii (5.231524) is larger than rlist (1.70) This doesn't even make sense since my as system itself fits in a 3.2 nm cubic box. I am pretty sure my topology is correct. I have seen on the mailing list that there have been problems with this but they have all been for bigger molecules. Are any of the molecules broken? That was the exact case just a couple of days ago: http://lists.gromacs.org/pipermail/gmx-users/2011-January/057980.html -Justin Gavin -- 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
Re: [gmx-users] sum of largest charge group radii
Then is this a bug?, or do I have to have a box in which all the molecules are whole. Does that not defeat the purpose of having a periodic box? Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin By broken do you mean split at the edge of the periodic box? If so then yes. Then that's your problem. -Justin Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi I have installed the new version of gromacs and so far I only use it as a check when running grompp as it gnerates more information about cut off charge group radii etc. For a system of 256 hexane molecules in a cubic box I get the following error. WARNING 1 [file pbc.mdp]: The sum of the two largest charge group radii (5.231524) is larger than rlist (1.70) This doesn't even make sense since my as system itself fits in a 3.2 nm cubic box. I am pretty sure my topology is correct. I have seen on the mailing list that there have been problems with this but they have all been for bigger molecules. Are any of the molecules broken? That was the exact case just a couple of days ago: http://lists.gromacs.org/pipermail/gmx-users/2011-January/057980.html -Justin Gavin -- 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
Re: [gmx-users] sum of largest charge group radii
Yeah I did that and there was no error, cheers. I did the same thing a week a go for a united atom model of hexane, I did not get the warning then even though the molecules were still broken. Peculiar? Justin A. Lemkul wrote: Gavin Melaugh wrote: Then is this a bug?, or do I have to have a box in which all the molecules are whole. Does that not defeat the purpose of having a periodic box? In a sense, maybe this is a grompp bug, in that it does not account for periodicity. There is no harm in using this run input file, though, which you can verify by making the molecules whole and then re-running grompp. You shouldn't see the same error. -Justin Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi Justin By broken do you mean split at the edge of the periodic box? If so then yes. Then that's your problem. -Justin Gavin Justin A. Lemkul wrote: Gavin Melaugh wrote: Hi I have installed the new version of gromacs and so far I only use it as a check when running grompp as it gnerates more information about cut off charge group radii etc. For a system of 256 hexane molecules in a cubic box I get the following error. WARNING 1 [file pbc.mdp]: The sum of the two largest charge group radii (5.231524) is larger than rlist (1.70) This doesn't even make sense since my as system itself fits in a 3.2 nm cubic box. I am pretty sure my topology is correct. I have seen on the mailing list that there have been problems with this but they have all been for bigger molecules. Are any of the molecules broken? That was the exact case just a couple of days ago: http://lists.gromacs.org/pipermail/gmx-users/2011-January/057980.html -Justin Gavin -- 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] Pull code
Dear All Sorry wrong subject title in previous post. Can someone please tell me how to generated the plot of mean force having ran the pull code at several distances using constraints? Gavin -- 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
