[gmx-users] Re: Steered MD
First of all, sorry for the late answer. More generally it is said the the rupture force depends logarithmically on the loading rate (velocity times spring constant). - So the rup.force should also depend logarithmically on the spring constant (in the simple bell modell). Think the reason, why most people vary the velocity is, it is easier. If you change the spring constant you also change the fluctuations in the force (higher constant - higher fluctuations) and it becomes way more difficult to determine a rupture force. For the dependence of the spring constant comes one paper to my mind: S. Izrailev, S. Stepaniants, M. Balsera, Y. Oono, and K. Schulten Biophysical Journal Volume 72 April 1997 1568-1581 Molecular Dynamics Study of Unbinding of the Avidin-Biotin Complex This might be a good starting point to find more recent publications on this topic. Hope this helps you. Greetings Thomas Am 09.10.2013 15:59, schrieb gmx-users-requ...@gromacs.org: Dear Gromacs Users, I am trying to look for some references regarding the SMD. I found one which tells about logarithmically dependency between the Rupture force (maximum pulling force) obtained from SMD and the pulling rate. Just wondering whether you are aware (or tested) the dependency between rupture force and harmonic spring constant? Steven -- 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: Umbrella Sampling _ pulled ion
Think your .mdp-file looks reasonable. If you are totaly unsure you could determine the PMF of two water molecules. As a reference one can use the radial distribution function g(r) and calculate the PMF as V_pmf(r) = -kT LN(g(r)) One side note: Since the force constant you are using is relative large you will need a relative small distance increment between the individual umbrella windows. Greatings Thomas Am 25.07.2013 13:56, schrieb gmx-users-requ...@gromacs.org: After running for 100 ps, I visualized the pullf_umbrella.xvg, in this plot I found the F value around 100 kJ/mol/ns. But force constant which I had set in md_US.mdp file was 4000 KJ/mol/ns. Does this difference show me that the US has not done properly? ;define??? ??? = -DPOSRES ??? integrator? = md??? ; leap-frog integrator nsteps? =10 ; 1 * 10 = 100 ps dt? = 0.001 ; 1 fs nstcomm = 10 tinit?? = 0 ; Output control nstxout = 5000?? ; save coordinates every 100 ps nstvout = 5000?? ; save velocities every 100 ps nstenergy?? = 1000?? ; save energies every 2 ps nstfout = 5000 nstxtcout?? = 5000??? ??? ??? ??? ; every 10 ps continuation??? = yes??? ; first dynamics run constraint_algorithm = lincs??? ; holonomic constraints constraints = h-bonds ; all bonds (even heavy atom-H bonds) constrained lincs_iter? = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid? ; search neighboring grid cells nstlist = 5 ; 10 fs rlist?? = 1.2?? ; short-range neighborlist cutoff (in nm) rlistlong?? = 1.4 rcoulomb??? = 1.2?? ; short-range electrostatic cutoff (in nm) rvdw??? = 1.2?? ; short-range van der Waals cutoff (in nm) vdwtype = switch rvdw_switch = 1.0 ; Electrostatics coulombtype = PME?? ; Particle Mesh Ewald for long-range electrostatics pme_order?? = 4 ; cubic interpolation fourierspacing? = 0.12? ; grid spacing for FFT fourier_nx? = 0 fourier_ny? = 0 fourier_nz? = 0 ewald_rtol? = 1e-5 optimize_fft??? = yes ; Temperature coupling is on tcoupl? = Nose-Hoover ; modified Berendsen thermostat tc-grps = Protein_POPC Water_Ces_CL??? ??? ; two coupling groups - more accurate tau_t?? = 0.5?? 0.5?? ??? ; time constant, in ps ref_t?? = 310?? 310??? ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl? = Parrinello-Rahman??? ; no pressure coupling in NVT pcoupltype? = semiisotropic tau_p?? = 4 ref_p?? = 1.01325 1.01325 compressibility = 4.5e-5 4.5e-5 ; Periodic boundary conditions pbc = xyz?? ; 3-D PBC ; Long-range dispersion correction DispCorr??? = EnerPres ; Pull code pull??? = umbrella pull_geometry?? = position pull_dim??? = N N Y pull_start? = yes pull_ngroups??? = 1 pull_group0 = POPC pull_group1 = Ces_ion pull_init1? = 0.0 0.0 0.0 pull_rate1? = 0.0 0.0 0.00 pull_vec1??? = 0 0 1 pull_k1 = 4000? ; kJ mol^-1 nm^-2 pull_nstxout??? = 1000? ; every 2 ps pull_nstfout??? = 1000? ; every 2 ps ; Velocity generation gen_vel = no?? ; assign velocities from Maxwell distribution refcoord_scaling??? = com Would you please let me know your suggestions ? -- 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: snapshot
Hi, the following: trjconv -s md.tpr -f md.trr -o beta.gro -dt 250 -sep should work. -dt 250 ; write output every 250 ps -sep ; to write each frame (but the output frequency gets overwritten by '-dt 250') the additional use of '-dump x' might be the problem, since with '-dump x' GMX writes only the frame at time x Greetings Thomas Am 15.07.2013 15:03, schrieb gmx-users-requ...@gromacs.org: Hi, The following command I used for the snapshots but I'm getting one frame time 0 time 500. what are the changes I need to do in the command to get snapshots in equal intervals of time for ex: 250ps. My production MD trajectory was 15ns. g_trjconv -s md.tpr -f md.trr -o beta.gro -skip 1 -dt 0 -dump 0 -sep Thanks in Advance. Rama On Mon, Jul 15, 2013 at 2:02 AM, Dr. Vitaly Chaban vvcha...@gmail.comwrote: look for either -dt or -skip. Dr. Vitaly V. Chaban On Mon, Jul 15, 2013 at 2:03 AM, Ramaramkishn...@gmail.com wrote: Hi, How to get a snapshots in equal intervals of time (250ps) from production MD trajectory. I'm using -sep , -t0, -timestep but output came only one .gro file. -- 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: Umbrella Sampling settings
But you need each of these lines for both cases (SMD and US). Probably one could skip two lines and use the default values, but it's better to set them manually. See below for comments (comments are under the related entry): Thanks for your reply. But when I don't understand why these extra lines are needed to set when are not advantageous practically! :-( Sincerely, Shima - Original Message - From: Justin Lemkul jalem...@vt.edu To: Shima Arasteh shima_arasteh2...@yahoo.com; Discussion list for GROMACS users gmx-users@gromacs.org Cc: Sent: Friday, July 12, 2013 1:37 AM Subject: Re: [gmx-users] Umbrella Sampling settings On 7/11/13 4:21 PM, Shima Arasteh wrote: Hi, I want to run Umbrella Sampling on my system. In initial configurations, an ion is located in center of the window. Some mdp file settings for running US, as I found in US tutorial are : ; Pull code pull? ? ? ? ? ? = umbrella how do you want to pull (umbrella / constant force / ...) pull_geometry? = distance second part for 'how to pull' ; see the manual pull_dim? ? ? ? = N N Y in which direction do you want to pull pull_start? ? ? = yes should the initial distance of 'pull_group0' and 'pull_group1' be added to 'pull_init1' pull_ngroups? ? = 1 number of pulled groups pull_group0? ? = Chain_B referene group pull_group1? ? = Chain_A pulled group pull_init1? ? ? = 0 place where origin of the potential is place pull_rate1? ? ? = 0.0 pulling veloity (default value could be 0, but if you want to pull with a zero-velocity it's safer to give this value manually pull_k1? ? ? ? = 4000? ? ? ; kJ mol^-1 nm^-2 force constant for umbrella / force for constant force pull_nstxout? ? = 1000? ? ? ; every 2 ps output coordinates pull_nstfout? ? = 1000? ? ? ; every 2 ps output forces So, you see each line is sensible, there are no 'extra' lines. Greetings Thomas But I'd like to know which lines are specifically for US? Because in this step, no group is supposed to be pulled but there are some lines written here related to pulling! All of them are related to umbrella sampling.? Pulling (steered MD) and umbrella sampling simply use common parts of the pull code in Gromacs because US requires a restraint potential.? Whether or not that restraint potential induces net displacement (steering, i.e. non-zero pull_rate) or not (zero pull rate, restrain to a given set of conditions) is the only difference.? Both processes require reference and pull groups, geometry information, etc. -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] Re: Umbrella Sampling settings
In GROMACS groups are called via the *.ndx file (default: index.ndx) Be aware that 'pull_dim' determines in which diretions (x,y,z) the umbrella potential acts. So use N N Y , if you want that the ion can move freely (onsidering the pull) in the xy-plane and Y Y Y if you want to also restrit the movement in the xy-plane. Am 12.07.2013 17:32, schrieb gmx-users-requ...@gromacs.org: Allright. As I said earlier, my system is a lipid bilayer. A channel is inserted in it and I want to run US on this system. An ion is considered in center of the each window, the reaction coordinate is set to z,? so the group which is pulled is an ion, and my ref group would be COM of the protein. But I don't know what statement is supposed to write in mdp settings exactly: ; Pull code pull??? = umbrella pull_geometry?? = position pull_dim??? = N N Y pull_start? = yes pull_ngroups??? = 1 pull_group0 = COM of protein pull_group1 = ion pull_init1? = 0 pull_rate1? = 0.0 pull_k1 = 4000? ; kJ mol^-1 nm^-2 pull_nstxout??? = 1000? ; every 2 ps pull_nstfout??? = 1000? ; every 2 ps IN fact, to implement such settings, how I make the US understand to get the COM of protein as the ref group and the proposed ion as the pulled group? Would you please give me any suggestions? Thanks for all your time and consideration. Sincerely, Shima - Original Message - From: Justin Lemkuljalem...@vt.edu To: Discussion list for GROMACS usersgmx-users@gromacs.org Cc: Sent: Friday, July 12, 2013 1:41 AM Subject: Re: [gmx-users] Umbrella Sampling settings On 7/11/13 5:10 PM, Shima Arasteh wrote: Thanks for your reply. But when I don't understand why these extra lines are needed to set when are not advantageous practically!:-( There's nothing extra.? Everything here has a functional purpose. -Justin Sincerely, Shima - Original Message - From: Justin Lemkuljalem...@vt.edu To: Shima Arastehshima_arasteh2...@yahoo.com; Discussion list for GROMACS usersgmx-users@gromacs.org Cc: Sent: Friday, July 12, 2013 1:37 AM Subject: Re: [gmx-users] Umbrella Sampling settings On 7/11/13 4:21 PM, Shima Arasteh wrote: Hi, I want to run Umbrella Sampling on my system. In initial configurations, an ion is located in center of the window. Some mdp file settings for running US, as I found in US tutorial are : ; Pull code pull? ? ? ? ? ? = umbrella pull_geometry?? = distance pull_dim? ? ? ? = N N Y pull_start? ? ? = yes pull_ngroups? ? = 1 pull_group0? ?? = Chain_B pull_group1? ?? = Chain_A pull_init1? ? ? = 0 pull_rate1? ? ? = 0.0 pull_k1? ? ? ?? = 4000? ? ? ; kJ mol^-1 nm^-2 pull_nstxout? ? = 1000? ? ? ; every 2 ps pull_nstfout? ? = 1000? ? ? ; every 2 ps But I'd like to know which lines are specifically for US? Because in this step, no group is supposed to be pulled but there are some lines written here related to pulling! All of them are related to umbrella sampling.? Pulling (steered MD) and umbrella sampling simply use common parts of the pull code in Gromacs because US requires a restraint potential.? Whether or not that restraint potential induces net displacement (steering, i.e. non-zero pull_rate) or not (zero pull rate, restrain to a given set of conditions) is the only difference.? Both processes require reference and pull groups, geometry information, etc. -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] Re: Umbrella sampling- force vs time plots
Generally: Using a higher force constant and / or pulling velocity drives the system faster out of equilibrium, which results in higher rupture forces. Varying the force constant has two effects. The softer the potential is, the larger are the fluctuations in the coordinates but the lower are the fluctuations in the force (see Paper: Biophysical Journal Volume 72 April 1997 1568-1581). So a lower force constant gives 'nicer looking' force-extension curves. Problem with a low force constant is that it's harder to detect intermediates in the system (especially if they have small energy barriers). I think the was a paper from Grubmüller which showed this in a picture. But i don't remember the title and only believe it was from Grubmüller. But one this effect also in the PMF obtained from Jarzynski's equality - here a too low force constant (lower than the curvature of the PMF) gives deviations from the exact result (G. Hummer, A. Szabo; PNAS 107 (50): 21441-21446, 2010). Hello all, I am trying to understand the force vs time plots using Gromacs' umbrella sampling method. I am trying to pull a short polymer chain from the interior of a micelle and see what the PMF looks like. I use the following parameters to run the pulling simulation for 500ps to pull the polymer over a distance of 5nm: pull=umbrella pull_geometry=direction pull_vec1=1 0 0 pull_start=yes pull_ngroups=1 pull_group0=surf pull_group1=poly pull_rate1=0.01 pull_k1=1000 After the simulation, pullf.xvg plot I obtained is a linearly increasing plot with time and similar result when pull_rate1=0.001 nm per ps. I am not The slope and the shape of the force-extension curve (extension = v*t) depends on the underlying potential surface. A linear increasing force corresponds (in first approximation) to an harmonic potential. In experiments where they use sometimes long linker molecules, one observes a force which increases steeper with increasing extension - this corresponds to a Worm-Like-Chain. sure if this is right. My question is, on what basis do we select the optimum pull_rate1 and pull_k1 for a particular system? Or is it just a choice of parameters as long as the system does not deform? How does an ideal force-time plot look like and does the choice of pull_k1 affect the histogram? It appears, the entire procedure depends on the choice of input of these two variables. I would greatly appreciate if someone can explain this concept. Thanks a lot. Andy Greetings Thomas -- 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: Is non-linear data output/storage possible?
I see the same problem as Justin. The real problem lies here: mdrun -nt 4 -v -s 10step.tpr -c nonlinear.gro -e nonlinear.edr -x nonlinear.xtc -o nonlinear.trr -g nonlinear.log -cpi nonlinear.cpt -cpo nonlinear.cpt -rdd 2.5 Due to the '-cpi nonlinear.cpt' mdrun thinks the the time now the end-time of the longer simulation. But in the grompp step, you tell GMX thaat the time starts at 0 and should go till the shorter end-time. You could use Justin's approach. Or skip the '-cpi nonlinear.cpt' in the mdrun step. Problem with the later is, that you'll end with many trajectories which all start at t=0. One thing which is strange: Why do you use in the second grompp/mdrun step, the *.tpr file of the previous run as input for the coordinates? Would think it will be better to use the final coordinates of the previous run and not the initial coordinates. (Being confused) Greetings Thomas Am 02.07.2013 20:51, schrieb gmx-users-requ...@gromacs.org: Hi, So the commands are: grompp -f 10step.mdp -p dppc_bilayer.top -o 10step.tpr -c dppc_bilayer.gro mdrun -nt 4 -s 10step.tpr -c nonlinear.gro -e nonlinear.edr -x nonlinear.xtc -o nonlinear.trr -g nonlinear.log -cpo nonlinear.cpt -rdd 2.5 grompp -f 100step.mdp -c 10step.tpr -p dppc_bilayer.top -o 100step.tpr -t nonlinear.cpt mdrun -nt 4 -s 100step.tpr -c nonlinear.gro -e nonlinear.edr -x nonlinear.xtc -o nonlinear.trr -g nonlinear.log -cpi nonlinear.cpt -cpo nonlinear.cpt -rdd 2.5 Up till this point the output is what I would expect. However now I try to use the 10step mdp again in order to save it every ten time steps - grompp -f 10step.mdp -c 100step.tpr -p dppc_bilayer.top -o 10step.tpr -t nonlinear.cpt mdrun -nt 4 -v -s 10step.tpr -c nonlinear.gro -e nonlinear.edr -x nonlinear.xtc -o nonlinear.trr -g nonlinear.log -cpi nonlinear.cpt -cpo nonlinear.cpt -rdd 2.5 And at this mdrun I get WARNING: This run will generate roughly 186617909559164928 Mb of data starting mdrun 'DPPC BILAYER' 100 steps, infinite ps (continuing from step 1000, 20.0 ps). This is the saving every 10 steps mdp file, the hundred step is identical except steps is a 1000 and nstxout is 100. ; VARIOUS PREPROCESSING OPTIONS = title= Martini cpp = /usr/bin/cpp ; RUN CONTROL PARAMETERS = integrator = md ; start time and timestep in ps = tinit= 0.0 This is your problem. You're no longer at t = 0. The checkpoint file specifies some time in the future, so grompp thinks you're going backwards in time and thinks you're going to generate infinite data by doing so. Then mdrun gets confused, because it thinks you're trying to go backwards as well, from step 1000 to step 100. Proper use of tinit and init_step in the .mdp file should alleviate this problem. -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] Re: Pulling
If you use periodic boundary conditions, there is no need that the protein stays at one side of the box. For the pulling simulation: Read the chapters 6.4 (explains the pull-code) and 7.3.21 (explains the mdp-paramters). Additional information you an also get from Justin Tutorial for Umbrella Sampling (http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/). Greetings Thomas Am 30.06.2013 01:07, schrieb gmx-users-requ...@gromacs.org: Dear users, I have a sytem including protein, lipids, and water. My protein is in center of the box. Now i want it stays at one side of the box. Which tool or command should i use to pull the protein to a any mong mu???n location ? Thankful for any help ! Thu -- 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: umbrella sampling for two polymer interaction
Sorry, my writing was not really excat. If you use a reference group, the force / potential acting on the pulled group is always relative to the reference group. If you use 'pull_geometry = distance' the origin potential is always in the distance 'pull_init1' along the vector from the reference group to the pulled group. (assuming the pulling velocity is zero, else there would be an additional contribution). Imagine an MD step. Both molecules have moved. 1) Calculate new origin of the umbrella potential 2) Calculate the forces 3) Let only these forces act which are allowed from 'pull_dim = X Y Z' Some cases: Y Y Y - Move the pulled group to the origin of umbrella potential Y N N - let the force act only along the X-axis N N N - let no force act 4) do next MD step and start again at (1) So if one component is N the pulled group can move freely along this axis (relative to the reference group). Just imagine the step of the pull-code as a step which happens after a MD step and remember that both step are independed from each other (ok, both use the coordinates from the previous step, but they are n that sense independend from each other as two following MD steps in a normal simulation). This means, that in the MD step, both groups can move freely, as there is now umbrella potential. But in the pull-code step, we switch the potential on. Sorry confused. In AMBER I remember I did once methane-methane interaction,just distance based umbrella sampling. But there I did not provide any direction. So should it not be N N N in Gromacs if I want to allow them freely move in xyz. With 'pull_gemoetry' you tell GROMACS how it should setup the pulling potential. With 'distance' it don't needs an explizit pull_vector, since GROMACS uses the vector connecting both groups. Am 31.05.2013 21:04, schrieb gmx-users-requ...@gromacs.org: Dear Thomas, Thanks a lot again for great reply. Please clarify this also, If one uses only 'Y N N' B would only move along the x-axis due to the pull, but could move freely in the yz-plane You never want to use the pull-code and 'pull_dim = N N N' This would mean that there is no force acting between your two groups. Then one could have skipped using the pull-code... So keeping Y N N will allow free movement in yz plane. So if I want A B to move freely in xyz but just keep them separated by some distance with spring const (just like two balloons tied to each other flying in air). Sorry confused. In AMBER I remember I did once methane-methane interaction, just distance based umbrella sampling. But there I did not provide any direction. So should it not be N N N in Gromacs if I want to allow them freely move in xyz. thanks, On Fri, May 31, 2013 at 8:53 PM, Thomas Schlesierschl...@uni-mainz.dewrote: For comments to your questions see below. More general: (somewhat longer than i wanted. Hope you find some answers here) Imagine two interacting particles A and B which are alinged to the x-axis. We take A as the reference group, B as pulled group and put the origin of the umbrella potential on top of B (pull_start=yes). Simulation starts - A and B moves. Pull-code step: From A we calculate the new position of the umbrella potential, this is unequal to B, since B moved and our reference point move. Now we have a force acting on B, 'pull_dim' controls in which directions the force acts. With 'Y Y Y' B is pulled towards the origin of the umbrella potential (and with this to the position it should have relative to A). If one uses only 'Y N N' B would only move along the x-axis due to the pull, but could move freely in the yz-plane. In the end one would get a structure where A and B have the right distance in the x-axis but are miles away from each other in the yz-plane. Now imagine we pull B away from A. Since MD simulations normally separete the movement of the center of mass of the system, it would look like A and B would move away from a middle point. (Interchanging A and B should give the same results). Think in your case (doing umbrella sampling) the mdp-file you suggested would be most appropiate (with 'pull_start=yes' and 'pull_ngroups=1'). This gives you a potential which fixes the distance between the two proteins. One thing you should be aware is that if you restrain the distance in 3d, you have to account for entropic effects (see also the GROMACS manual). If you restrain the system only in one direction, these don't arise. Think this is the reason why one sees some work with umbrella sampling were the restraint works only in one direction. Am 31.05.2013 17:20, schriebgmx-users-requ...@gromacs.org: Dear Thomas, Thanks a lot for your time and nice explanation. I was not able to get specially the pull_start flag but now its quite clear. I feel sorry, that should be pull_dim = N N N in my case. Also I will be much thankful if you please can help me to make understand following: STOP!!! You never want to use the pull-code and 'pull_dim = N
[gmx-users] Re: umbrella sampling for two polymer interaction
Look also into the manual. But the tutorial is a nice place to start. For further comments see below: Dear Lloyd, I have read that but my system is different regards, On Thu, May 30, 2013 at 8:28 PM, lloyd riggslloyd.ri...@gmx.ch wrote: Dear Jiom, Look at justines tutorial, there's example pull .mdp. Stephan Watkins I want to do Umbrella sampling between two different polymers (A and B) interacting with each other with starting configuration separated by some distance and I am trying to bring them closer. I have some queries regarding pull inputs: (this is for to run a umbrella sampling at some distance) pull = umbrella pull_geometry = distance pull_dim = Y Y Y pull_start = ??? pull_ngroups = 2? pull_group0 = polymer_B pull_group1 = polymer_A pull_init1 = 0 pull_rate1 = 0.0 please suggest for following: 1) pull_dim I have set to Y Y Y: Is this correct I do not want to make it interact with some directional vector I don't really understand the question. If you use 'pull_dim = Y Y Y' the pulling potential acts in all 3 dimensions, if you use 'pull_dim = Y N N' the pulling potential acts only in the X direction and your two groups an move freely in the YZ-plane. 2) Which should be group0 or group1, in other words should I pull both together or how I should decide which one should be reference and which to be pulled as both are different polymers? Depends on what you want to do. Easiest way would be define one polymer a group0/reference group and the other as group1/pulled group. System shouldn't care about which polymer is which group. If you do a pulling simulation, there can be reason for chosing the groups (protein = reference , ligand = pulled group, since we want to pull it away) 3) And also what should be pull_ngroups because if there is no reference group then it should be 2 Better use a reference group - pull_ngroups = 1 You don't want to pull in absolute coordinates, when your system can rotate... 4) I am not able to understand pull_start option with pull_init1. In this case if it is set to yes and 0.0 respectively then does that mean this combination is equivalent to pull_start = No if I just assume pull_init1 does not have any default value (which is 0.0); not existing From the setup which you have written above: polymer_B is the reference group. the origin of the pulling potential is at 'pull_init1' (a length) along the vector which goes from polymer_B to polymer_A (sine you use pull_geometry = distance). If you set pull_init1=0 and pull_start=no, polymer_A will crash into polymer_B (since the origin of the umbrella potential is directly at the center of mass of polymer_B). If you set pull_init1=0 and pull_start=yes, GROMACS adds the distance between polymer_B and A to pull_init1 (- so pull_init1 is now greater than 0.0). Now the origin of the umbrella potential is at the center of mass of polymer_A. - A is restrainted to a certian distance of B. 5) Also finally where are upper and lower bounds defined. pull_k1 = 1000 is harmonic applied to some equilibrium distance value. How this distance is taken by the programme (or it is just the starting distance taken between two groups) and what are the ± values defined. (say in AMBER I define r1,r2,r3,r4; where r2=r3 which is assumed equilibrium value and r1 is lower and r4 is upper value which defines shape of potential) The umbrella potential is a simple harmonic potential (so no fancy stuff as in AMBER) with V = 1/2 k x^2 where x is the violation of the equilibrium distance. For your setup V = 1/2 (pull_init1 - distance(B-A))^2 where distance(B-A) means the distance between both polymers. Greetings Thomas -- 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: umbrella sampling for two polymer interaction
For comments to your questions see below. More general: (somewhat longer than i wanted. Hope you find some answers here) Imagine two interacting particles A and B which are alinged to the x-axis. We take A as the reference group, B as pulled group and put the origin of the umbrella potential on top of B (pull_start=yes). Simulation starts - A and B moves. Pull-code step: From A we calculate the new position of the umbrella potential, this is unequal to B, since B moved and our reference point move. Now we have a force acting on B, 'pull_dim' controls in which directions the force acts. With 'Y Y Y' B is pulled towards the origin of the umbrella potential (and with this to the position it should have relative to A). If one uses only 'Y N N' B would only move along the x-axis due to the pull, but could move freely in the yz-plane. In the end one would get a structure where A and B have the right distance in the x-axis but are miles away from each other in the yz-plane. Now imagine we pull B away from A. Since MD simulations normally separete the movement of the center of mass of the system, it would look like A and B would move away from a middle point. (Interchanging A and B should give the same results). Think in your case (doing umbrella sampling) the mdp-file you suggested would be most appropiate (with 'pull_start=yes' and 'pull_ngroups=1'). This gives you a potential which fixes the distance between the two proteins. One thing you should be aware is that if you restrain the distance in 3d, you have to account for entropic effects (see also the GROMACS manual). If you restrain the system only in one direction, these don't arise. Think this is the reason why one sees some work with umbrella sampling were the restraint works only in one direction. Am 31.05.2013 17:20, schrieb gmx-users-requ...@gromacs.org: Dear Thomas, Thanks a lot for your time and nice explanation. I was not able to get specially the pull_start flag but now its quite clear. I feel sorry, that should be pull_dim = N N N in my case. Also I will be much thankful if you please can help me to make understand following: STOP!!! You never want to use the pull-code and 'pull_dim = N N N' This would mean that there is no force acting between your two groups. Then one could have skipped using the pull-code... 1) If you do a pulling simulation, there can be reason for chosing the groups (protein = reference , ligand = pulled group, since we want to pull it away) This indeed is correct but I am not able to get depth of this. I mean to say lets keep ligand as a reference and protein as pulled group then yes it sounds stupid but I am not able to provide a reason myself why we can not keep ligand as reference and pull protein rather !! Think this setup should also work. For some simple systems i imagine it should give identical results. For complex system i would also think so. But i can't comment on this with actual expirience. The dimer systems which i investigated were symmetric... 2) 3) And also what should be pull_ngroups because if there is no reference group then it should be 2 Better use a reference group - pull_ngroups = 1 You don't want to pull in absolute coordinates, when your system can rotate.. I am not able to understand this part. Can you please provide some example so that it makes easier to understand this Imagine only a single protein which you want to unfold. In an equilibrium simulation the protein can freely rotate in the box. If we use the N-terminus as the reference group and the C-terminus a the pulled group, the origin of the umbrella potential will always be updated and will account for movement of the N-terminus (reference group). If one would pull in absolute coordinates, one would need to give the position of the umbrella potential in absolute space. The molecule can move, but the origin of the potential will always stay fixed at one place. Think in the end this would be equal to an position restraint of said group. If one would want to restrain the distance of two groups in such a way, one would need two umbrella potentials. But since these would be equal to two position restraints, there would be no coupling between the two. I mean, if both groups move around but would have the same distance it should be ok since the distance is fine. But both umbrella potential would pull each group back to the initial position. Much thanks again, regards, Jiom On Fri, May 31, 2013 at 1:21 PM, Thomas Schlesierschl...@uni-mainz.dewrote: Look also into the manual. But the tutorial is a nice place to start. For further comments see below: Dear Lloyd, I have read that but my system is different regards, On Thu, May 30, 2013 at 8:28 PM, lloyd riggslloyd.ri...@gmx.ch wrote: Dear Jiom, Look at justines tutorial, there's example pull .mdp. Stephan Watkins I want to do Umbrella sampling between two different
[gmx-users] Re: Running Pull Code
The three steps (EM, NVT and NPT) are to equilibrate the system. How much time these steps need depends on the system. But i would assume a ouple of nanosecounds are sufficient for most systems. You could look into the literature, how long other people equilibrate systems which are similar to ours. If the system is equilibrated, you an start to perform the pulling simulation to obtain the individual structure for the later umbrella sampling. Greetings Thomas Am 17.05.2013 07:46, schrieb gmx-users-requ...@gromacs.org: Hi, I have a system composed of POPC/peptide/water/ions. I aim to study ion conduction through the peptide using umbrella sampling. I built the system and ran EM, NVT, NPT successfully, but have not run md yet. I' d like to know if the system is required of passing a few nanoseconds md? Or I might be able to go to Umbrella Sampling straight after NPT? As I studied in Justin's tutorial, running pull code is done after some typical steps of every simulation ( EM, NVT, NPT). But I dont know if is correct generally for other systems as well? Would you please give me any suggestions? Thanks in advance. Sincerely, Shima? -- 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] Performance (GMX4.6.1): MPI vs Threads
Dear all, if one performs a parallel calculation on a single node / computer with more than 1 core, is there a speed difference between MPI and Threads? Problem is for gromacs 4.6.1 i'm facing problems to link it to OpenMpi. A serial version (without threads) worked, so i think i should be able to compile a version which supports threads. Since we only run calculations on a single node on our cluster (no infini-band), i'm wondering if the programm with threads would be sufficient in our case. Greetings Thomas -- 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: umbrella sampling
Hi, it seems that you've only coupled your protein to the thermostat, but not the solvent, hence the error message. Generally one would couple both domains of the protein to one thermostat and the solvent (inluding ions) to another thermostat. Side note: If you want to use the WHAM method for constructing the PMF from the simulations, you need to use an harmonic potential instead of constraints. This method is named 'umbrella sampling'. Using constraints and integrating the constraint force to obtain the PMF is called 'thermodynamic integration'. Greetings Thomas Am 13.05.2013 10:33, schrieb gmx-users-requ...@gromacs.org: HI, I would like to compute an umbrella sampling simulation for o protein divided in two domain, with Center of mass pulling using as constraint between the two domains. And the constraint is applied instead of a harmonic potential I create a pull.mdp file with this option for the pull: title = Umbrella pulling simulation ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 25; 500 ps nstcomm = 10 ; Output parameters nstxout = 5000 ; every 10 ps nstvout = 5000 nstfout = 500 nstxtcout = 500 ; every 1 ps nstenergy = 500 ; Bond parameters constraint_algorithm= lincs constraints = all-bonds continuation= yes ; continuing from NPT ; Single-range cutoff scheme nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb= 1.4 rvdw= 1.4 ; PME electrostatics parameters coulombtype = PME fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 optimize_fft= yes *; Berendsen temperature coupling is on in two groups Tcoupl = Nose-Hoover tc_grps = domain_1 domain_2 tau_t = 0.50.5 ref_t = 310310* ; Pressure coupling is on Pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 1.0 compressibility = 4.5e-5 ref_p = 1.0 refcoord_scaling = com ; Generate velocities is off gen_vel = no ; Periodic boundary conditions are on in all directions pbc = xyz ; Long-range dispersion correction DispCorr= EnerPres *; Pull code pull= constraint pull_geometry = distance ; simple distance increase pull_dim= Y Y Y pull_start = yes ; define initial COM distance 0 pull_ngroups= 2 pull_group0 = domain_1 pull_group1 = domain_2 pull_rate1 = 0.01 ; 0.01 nm per ps = 10 nm per ns pull_k1 = 1000 ; kJ mol^-1 nm^-2* I don't know if I use the good option for compute this king of umbrella sampling. When I run the grompp command: *grompp -f md_pull.mdp -c min1.gro -p topol.top -n index.ndx -o pull.tpr* I optain this error: * Fatal error: 108147 atoms are not part of any of the T-Coupling groups* I am little confuse with the mdp file option. I don't know which name I should clarify for T-Coupling groups ( in green) ?? If some one can help me . THanks a lot. Nawel -- 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: RDF - water and protein
Hi, the manual mentions that with the option '-surf' no normalization is done. So it's normal the RDF will inrease with larger distances, since the further you'll go away from the protein, the bigger the spherical shell is (from which the RDF for distance r is calculated) and the more water molecules will be in this shell. Because of this using the '-nopbc' option shouldn't change anything. One thing which i find rather strange, is that for really high distances the RDF goes down. Would have thought that with pbc the RDF should increase all the time. But i have no idea to perform a sensible normalization afterwards. Greetings Thomas Am 13.05.2013 17:24, schrieb gmx-users-requ...@gromacs.org: Dear Gmx Users, I run long simulation of my protein with 50 small molecules in water. I calculated the RDF (Protein - Water) using -surf mol and -rdf mol_com. Please, take a look at my plot: http://speedy.sh/tmJbD/rdf-P-W.png Could you please, explain me why the second peak is so high? Shall I use option -nopbc ? How can I obtain the density [kg/m3] versus the distance? I wish to compare it with RDF of Protein-Ligands so wish to use some normalization for them to be comparable - peaks will be visible. Steven -- 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: PMF calculation between protein and ligand
If you refer to Justin's tutorial, you can use it for every system, where the reaction coordinate is the distane between two parts of the system. You probably need to make some little adjustments, but in the end it doesn't matter if you pull two proteins apart, or protein+ligand, 2 water molecules, or even streth a protein. Greetings Thomas Am 02.05.2013 08:59, schrieb gmx-users-requ...@gromacs.org: Sir I have query as to how to go ahead for potential mean force (PMF) calculation between protein and ligand. As per the umbrella sampling protocol you have provided us in gromacs tutotrial it says it is for protein molecules... -- 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: how is the pulling force measured
I never looked into the code, but i understanded it the following way (for pull_geometry = distane, direction, position). For the pulling the reference group is held fixed and the force is applied to the pulled group. Since the reference group an move during the rest of the simulation (all steps not involving pulling), the pull-stuff is updated in every step (meaning if you use 'pull_geometry = position' you're pulling to a point relative to the reference group. If the referene group moves, the point to where you pull also moves - so pulling is always relative to the reference group, which leads to the fact the reference group is fixed for the pulling). If you simulate with 'pbc' you should separate the movement of the center of mass (com) of the system. If you look into a video of the trajectory (vmd or similar program), it will look like that both groups will be pulled. But this is from the separation of the com-movement. If you want a 'nice' movie you can use 'g_trjconv' (i think) to define a group which is always in the enter of the box. One important side note: If you pull only in one dimension (x,y or z), there will be no pulling-force acting on the other two directions. In your case your pulled-group could move freely in the xz-plane. In my opinion it's safer to pull in all three dimensions. Greetings Thomas Am 30.04.2013 09:37, schrieb gmx-users-requ...@gromacs.org: Hi Alex, I read the manual but I got confused with the details. My pull_geometry = position with pulling direction along y. I am trying to dissociate two interacting proteins. So I set the last amino acid of one protein as my reference group and the last amino acid of the other protein as the pulled group. Now, I don't know exactly how to interpret what reference group and pulled group mean, and what is their difference; since in the manual it says 'there is no difference in treatment of the reference and pulled group (except with the cylinder geometry).' Also, it says that measured force output is the force of the pulled group. So basically, I am trying to visualize how the system works. Using the pull_geometry stated above, are the two groups (reference and pulled) attached to two virtual springs? Will this statement be correct? Force was applied by moving the clamped ends of the two springs (i = 1, 2) in opposite directions with constant velocity v to positions Zi(t) = zi(0) + vt, where zi(0) is the initial position of the end of the molecule. The forces fi(t) at the two ends were measured at every time step using the relationship fi(t) = k (zi(t) ??? Zi(t)). I hope my questions make sense. Thank you. -- 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: gmx-users Digest, Vol 108, Issue 165
TL;DR version (longer version below): Due to the stochastic nature of SMD (and pulling experiments in general) it is quite natural that the results for different simulations will not be excatly the same. I would say you are fine. Thing is, if you do ligand unbinding SMD simulations the unbinding process is a stochastic process and the results differ somewhat, most importantly the rupture force (highest force before the rupture event) will differ between the simulations and you'll get a distributions of rupture forces. In a simple two-state model, the slope of force vs time, should be the same for each simulation before and after the rupture event (but probably different slopes before and after the rupture), only the time when the rupture event happens will change. If your system is more complex (say a couple of bonds must be broken before the ligand unbinds), the slope before the rupture may vary between some simulations (in some simulations some of the bonds break earlier than in other simulatons. In http://pubs.acs.org/doi/full/10.1021/jp3115644 i looked into a system which can be described by a two-state model. The results from one typical trajectory are shown, but mostly i discuss the averaged results. If one compares the single trajectory and the averaged, one clearly sees that the trajectories varies around the rupture event, but at the start of the simulations the system behaves almost equal for each simulation. Hope this helps. Greetings Thomas Am 26.04.2013 12:00, schrieb gmx-users-requ...@gromacs.org: Dear Users, I am running my puling simulations of ligand with constant velocity. First I minimize and equilibrate my system: grompp -f minim.mdp -c Solvions.gro -p topol.top -o em.tpr mdrun -s em.tpr -deffnm em grompp -f nvt298US.mdp -n index.ndx -c em.gro -p topol.top -o nvt298.tpr mdrun -s nvt298.tpr -deffnm nvt298 grompp -f npt298US.mdp -n index.ndx -c nvt298.gro -t nvt298.cpt -p topol.top -o npt298.tpr mdrun -s npt298.tpr -deffnm npt298 Then I run 10 pulling simulations with the same mdp file: grompp -f pull.mdp -c npt298.gro -p topol.top -n index.ndx -t npt298.cpt -o pull_1.tpr mpiexec mdrun -s pull_1.tpr -deffnm pull_1 ... grompp -f pull.mdp -c npt298.gro -p topol.top -n index.ndx -t npt298.cpt -o pull_10.tpr mpiexec mdrun -s pull_10.tpr -deffnm pull_10 I get 3 different (but similar) profiles (Force vs time) with 10 simulations as some of them produce exactly the same results... In another system with the same methodology I get 10 similar but different profiles. I am wondering why in this case only 3 types are possible... Shall I try grompp without -t npt.cpt ? Steven -- 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: SMD - reproducibility
Think i now understand your question. Forget what i wrote before. I could imagine the the 'grompp -t npt.cpt' part is a problem. If the simulations would be numerical reproducible, one should get the same results. As they are not, the results will differ somewhat (would think the more, the longer you simulate). But two different simulations would be more equal to each other, than two simulations which start with different velocity distributions for the particles. If you're interested in an stochastic analysis of your system (meaning simulations which are not equal - performing many pulling experiments in reality, one would also have many different starting points) you could do two things: 1) Run a look npt simulation, and use different frames to start the SMD simulations. From each frame the particles should have a different velocity distribution and the results of the SMD simulations should also differ. (Depending on how many SMD simulations you want to perform, this might get expensive, since the starting frame for SMD should be separated by more then a few ps.) 2) Dump the 'npt.cpt' file and randomly determine new velocities for each particle at the start of each SMD simulation. Since each simulation has a different velocity distribution, the SMD simulation wont be the same. This approach has only one weak point. Due to assigning new random velocities you destroy the thermal equilibrium of the system. But if the system was well equilibrated before, this distrubance should only be small and after the first 100-200 ps of the SMD simulaton the system is in thermal equilibrium. If the complete SMD simulation is much longer (couple of ns), the interesting stuff would happen longer after the inital simulation time with the destroyed equilibrium. Hope this helps Thomas Am 26.04.2013 12:00, schrieb gmx-users-requ...@gromacs.org: Dear Users, I am running my puling simulations of ligand with constant velocity. First I minimize and equilibrate my system: grompp -f minim.mdp -c Solvions.gro -p topol.top -o em.tpr mdrun -s em.tpr -deffnm em grompp -f nvt298US.mdp -n index.ndx -c em.gro -p topol.top -o nvt298.tpr mdrun -s nvt298.tpr -deffnm nvt298 grompp -f npt298US.mdp -n index.ndx -c nvt298.gro -t nvt298.cpt -p topol.top -o npt298.tpr mdrun -s npt298.tpr -deffnm npt298 Then I run 10 pulling simulations with the same mdp file: grompp -f pull.mdp -c npt298.gro -p topol.top -n index.ndx -t npt298.cpt -o pull_1.tpr mpiexec mdrun -s pull_1.tpr -deffnm pull_1 ... grompp -f pull.mdp -c npt298.gro -p topol.top -n index.ndx -t npt298.cpt -o pull_10.tpr mpiexec mdrun -s pull_10.tpr -deffnm pull_10 I get 3 different (but similar) profiles (Force vs time) with 10 simulations as some of them produce exactly the same results... In another system with the same methodology I get 10 similar but different profiles. I am wondering why in this case only 3 types are possible... Shall I try grompp without -t npt.cpt ? Steven -- 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: SMD - reproducibility
Don't know. One idea i have: Take a flexible and a relative rigid system and perform simulations with the same starting conditions (- using -t *.cpt). I would imagine that for the flexible system the trajectories start earlier to deviate, since more stuff could happen (system is more flexible - greater configurational space). For the rigid system the configurational space is smaller, so the probability is higher to always follow the same trajectory if one starts with a predefinded velocity and direction. But don't know if this is true, but it's the first thing which comes to my mind. Am 26.04.2013 14:01, schrieb gmx-users-requ...@gromacs.org: Thanks for this. I think option 2 is more reasonable. However, still do not know why I get sometimes 3 types of profiels and sometimes 10 for 10 SMD simulations... On Fri, Apr 26, 2013 at 12:46 PM, Thomas Schlesierschl...@uni-mainz.dewrote: Think i now understand your question. Forget what i wrote before. I could imagine the the 'grompp -t npt.cpt' part is a problem. If the simulations would be numerical reproducible, one should get the same results. As they are not, the results will differ somewhat (would think the more, the longer you simulate). But two different simulations would be more equal to each other, than two simulations which start with different velocity distributions for the particles. If you're interested in an stochastic analysis of your system (meaning simulations which are not equal - performing many pulling experiments in reality, one would also have many different starting points) you could do two things: 1) Run a look npt simulation, and use different frames to start the SMD simulations. From each frame the particles should have a different velocity distribution and the results of the SMD simulations should also differ. (Depending on how many SMD simulations you want to perform, this might get expensive, since the starting frame for SMD should be separated by more then a few ps.) 2) Dump the 'npt.cpt' file and randomly determine new velocities for each particle at the start of each SMD simulation. Since each simulation has a different velocity distribution, the SMD simulation wont be the same. This approach has only one weak point. Due to assigning new random velocities you destroy the thermal equilibrium of the system. But if the system was well equilibrated before, this distrubance should only be small and after the first 100-200 ps of the SMD simulaton the system is in thermal equilibrium. If the complete SMD simulation is much longer (couple of ns), the interesting stuff would happen longer after the inital simulation time with the destroyed equilibrium. Hope this helps Thomas Am 26.04.2013 12:00, schriebgmx-users-requ...@gromacs.org: Dear Users, I am running my puling simulations of ligand with constant velocity. First I minimize and equilibrate my system: grompp -f minim.mdp -c Solvions.gro -p topol.top -o em.tpr mdrun -s em.tpr -deffnm em grompp -f nvt298US.mdp -n index.ndx -c em.gro -p topol.top -o nvt298.tpr mdrun -s nvt298.tpr -deffnm nvt298 grompp -f npt298US.mdp -n index.ndx -c nvt298.gro -t nvt298.cpt -p topol.top -o npt298.tpr mdrun -s npt298.tpr -deffnm npt298 Then I run 10 pulling simulations with the same mdp file: grompp -f pull.mdp -c npt298.gro -p topol.top -n index.ndx -t npt298.cpt -o pull_1.tpr mpiexec mdrun -s pull_1.tpr -deffnm pull_1 ... grompp -f pull.mdp -c npt298.gro -p topol.top -n index.ndx -t npt298.cpt -o pull_10.tpr mpiexec mdrun -s pull_10.tpr -deffnm pull_10 I get 3 different (but similar) profiles (Force vs time) with 10 simulations as some of them produce exactly the same results... In another system with the same methodology I get 10 similar but different profiles. I am wondering why in this case only 3 types are possible... Shall I try grompp without -t npt.cpt ? Steven -- 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: Is it possible to pull 2 different groups with respect to a third group as reference?
Yes, from 4.5.x manual: pull_ngroups: (1) The number of pull groups, not including the reference group. [...] Just set 'pull_ngroups = 2' and then make entries for pull_group1 - pull_vec1 ... pull_group2 - pull_vec2 ... and so on... greetings thomas Am 12.04.2013 12:00, schrieb gmx-users-requ...@gromacs.org: I have 3 groups on the system, and I want to pull 2 groups with respect to the 3rd one. Is it possible in GROMACS? Thanking you, Saikat -- 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: fail to pull
Must get the bus, so only short answer. You could try to use constraint pulling instead of an umbrella potential. Then the ligand should move 1nm in 1ns. And you could see is the setup is ok, or if it would be better to pull into another direction... greetings thomas Am 07.04.2013 18:05, schrieb gmx-users-requ...@gromacs.org: Hard to tell. Does your ligand have a suitable exit pathway exactly aligned along the x-axis? Have you tried increasing the pull rate? How long is the simulation? I don't even see nsteps in the above .mdp file. How about increasing the force constant? Is the vector connecting the COM of the entire protein and the COM of the ligand suitable for describing the exit pathway? -Justin Hello Justin: thanks a lot for kind rely. Yes, I adjust the conformation of whole protein/ligand so that it can exist from X-axies. I only show part of the .mdp file so some of then are not shown. ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 50; 500 ps nstcomm = 10 ; Output parameters nstxout = 5000 ; every 10 ps nstvout = 5000 nstfout = 1000 nstxtcout = 1000 ; every 1 ps nstenergy = 1000 probably I should consider use part of the protein such as residues around binding pocket as COM for protein instead of whole protein? I applied for 1ns with rate pull_rate1= 0.001, so at then end of pulling, the distance for COMprotein and COM ligand should be 10A. Probably this is too short for whole protein as COM? Let me clear one thing up first. 1 ns of pulling with a 0.001 nm/ps pull rate will not necessarily cause the ligand to be displaced by 1 nm. The particle pulling the virtual spring will be displaced by 1 nm, but the ligand will only move as a function of this applied force and the restoring forces (i.e. interactions between the ligand and protein). Choosing a more suitable reference group and running the simulation for longer will produce the desired result. -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] SMD : pulling both primes of DNA
The easiest solution would be using 'pull_geometry = distance' in all three dimensions. Then you can be sure that both groups are pulled together. Small remark: One group would be fixed for the pulling, and the second group gets pulled towards the first group. So if you want to have both groups really close, this setup would be the easiest solution. If you really to pull both groups to a point in the middle between them, i think you would need 'pull_geometry = position': Use some part of the DNA as the reference group and put the origin of the pulling (via 'pull_init') to the middle point. Then make two pull-groups + 'pull_vec1' or 'pull_vec2', from each DNA part, which you want to pull to the middle point. In principle this setup should work, but you can run into problems if the structure of the DNA changes and the middle point, to where you pull both groups, moves somewhat around. hard to explain... I would stick to the first approach... What pull-mode did you use? You could use 'pull = constraint', measure the distance between both groups and then pull with a given pulling velocity so long till both groups should meet. Since you use a constraint instead of a restraining potential, the groups should be definitly pulled together (or some other parameters in the .mdp file are wrong). Greetings Thomas Am 27.03.2013 20:16, schrieb gmx-users-requ...@gromacs.org: Hello All, My simulation system is composed of DNA and want to pull both of the primes at the same time towards each other. I have tried every possible set of parameters in pull code...and i am not able to pull them together. it's like ... the DNA molecule is aligned along Y-axis and if want to pull both ends at the same time that means pulling in -Y and +Y direction at the same time. if anyone successfully implemented such ...simulations... Please provide me with some idea about this problem. Thank you Raghav -- 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 with large pulling distance (larger than half of the box size)
Look for pull_geometry = direction_periodic This should solve the problem. Greetings Thomas Am 20.03.2013 12:00, schrieb gmx-users-requ...@gromacs.org: Dear all, I want to use Umbrella sampling method to calculate the potential of mean force. Unfortunately, the distance of my two groups is very large (about 10nm, larger than the half of simulation box size 16nm). Without increasing the box size (too large!!), is there any ideas to solve this problem? In another way, is it possible to modify the pulling code (pull.c) to force the umbrella sampling method only calculate the distance of the groups in the box but not the periodic image? -- 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] Treating electrostatics and van der Waals interactions differently
Looking into the manual, i find under the 'energy monitor group': 'Mutual interactions between all energy monitor groups are compiled during the simulation. This is done separately for Lennard-Jones and Coulomb terms. In principle up to 256 groups could be defined, but that would lead to 256x256 items! Better use this concept sparingly. All non-bonded interactions between pairs of energy monitor groups can be excluded (see sec. 7.3). Pairs of particles from excluded pairs of energy monitor groups are not put into the pair list. This can result in a significant speedup for simulations where interactions within or between parts of the system are not required.' But i found nowhere an option to tell GROMACS to handle the exclusions for LJ and Coulomb interactions differently. Don't know if it's possible. For standard exclusions it seems that LJ and Coulomb interaction get treated equally - if one excludes one, one also excludes the other... One idea: But only possible if one can define two bonds, between the same atoms. I don't know if this is possible. One idea for a work-around would be the following (but it will slow the simulation). Exclude all 1-2 and 1-3 interactions. To get the Coulomb interacions back, construct table interactions for bonds. GMX manual 4.5.3: Look for '4.2.13 Tabulated interaction functions' and '6.7 Tabulated interaction functions' It seems that one can construct tabulated interactions also for bonds. Potential looks like: V = k * f(r) where f(r) is a cubic spline, think here one could use the values for the coulomb-interaction for the non-bonded table. For k use the product of the two charges. In the topology, one would need to write both atoms and this force-constant (our product of charges), similar to normal bonds. If two bonds between the same atoms are not possible, one could define tabulated interactions for the bonds + Coulomb interactions ... but this would be a lot of work. Probably it would be easier, to hack GMX, that it ignores Exclusions for Coulomb interactions. Or another persons has hopefully a better idea. Greetings Thomas Hi all, I am attempting to simulate a system that has strong ionic character so I would like to treat the electrostatics and van der Waals interactions separately. For example I would like to include all pairs of atoms in the electrostatics calculation but I would like to exclude 1-2 and 1-3 neighbors in the van der Waals calculation. Is this possible to do in GROMACS? If so, how might this be accomplished? Thanks in advance for your help. Jeff Woodford Assistant Professor of Chemistry Missouri Western State University -- 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: Temperature coupling problem: reference T is, different to given T in mdp file
Hi, you use a relative high tau_t: tcoupl = v-rescale tc-grps = CO2 PARA PAR2 tau_t = 1.0 1.0 1.0 ref_t = 410 410 410 In my simulations i use a value of 0.1 don't know if this would help, but it was the only thing which catched my eye, as i read you .mdp parameters. Greetings Thomas Am 22.01.2013 15:36, schrieb gmx-users-requ...@gromacs.org: Dear gmx-users, 1) I am running simulation of two paracetamol molecules dissolved in supercritical CO2 in NVT ensemble. The problem is that the resulting temperature of the system is larger than the temperature specified in mdp file by ~ 20 K (i set ref_t = 410 K in mdp file). Berendsen and v-rescale thermostats behave similarly. The output from g_energy: Statistics over 250001 steps [ 0. through 500. ps ], 1 data sets All statistics are over 2501 points Energy Average Err.Est. RMSD Tot-Drift --- T-System 430.3421.1 14.1796 6.30923 (K) 2) When i run the same simulation but set 3 different subsystems (all CO2 molecules, the first paracetamol molecule, the second paracetamol molecule) coupled to different bathes with the same temperature, the resulting temperatures do not match the specified temperature = 410 K. The output from g_energy: This approach is unsound. There are insufficient degrees of freedom in single molecules to justify distinct thermostats. Statistics over 250001 steps [ 0. through 500. ps ], 4 data sets All statistics are over 2501 points Energy Average Err.Est. RMSD Tot-Drift --- Temperature428.821.2 14.356 -3.39055 (K) T-CO2 427.432 1.2 14.6848 -5.24 (K) T-PARA452.434 3.7 91.2873 20.656 (K) T-PAR2454.211 10 93.7436 37.8375 (K) I would appreciate very much if someone gives me a hint what could be wrong. Considering all data points will include early frames that are not necessarily anywhere close to being equilibrated. You could try block averaging to see if the results converge to the correct result over the latter part of the simulation, or at least see if the system is approaching the target temperature. Some information about system: There are 200 CO2 molecules, cubic box: 4 x 4 x 4 nm3. I set coulombtype = Cut-off just for testing purposes, PME gives similar results. It would be better to diagnose the system using PME. Even if the results are similar, rounding errors are compounded when using plain cutoffs, which render the results inherently unreliable. -Justin Some mdp control options (for the 2) case): ; RUN CONTROL PARAMETERS = integrator = md dt = 0.002 nsteps = 25 comm-mode = Linear nstcomm = 1 ; NEIGHBORSEARCHING PARAMETERS = nstlist = 10 ns_type = grid pbc = xyz rlist = 1.1 ; OPTIONS FOR ELECTROSTATICS AND VDW = coulombtype = Cut-off rcoulomb = 1.1 vdw_type = Shift rvdw = 1.0 ; OPTIONS FOR WEAK COUPLING ALGORITHMS = tcoupl = v-rescale tc-grps = CO2 PARA PAR2 tau_t = 1.0 1.0 1.0 ref_t = 410 410 410 ; GENERATE VELOCITIES FOR STARTUP RUN = gen_vel = yes gen_temp = 410 gen_seed = 473529 ; OPTIONS FOR BONDS = constraints = hbonds constraint_algorithm = lincs unconstrained_start = no shake_tol = 0.1 lincs_order = 4 lincs_warnangle = 30 morse = no lincs_iter = 2 -- 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: MD Simulation of Calcium with protein
Am 15.01.2013 12:52, schrieb gmx-users-requ...@gromacs.org: On 1/15/13 5:49 AM, Devika N T wrote: HI I would like to know the protocol to be followed for performing MD simulation for calcium with protein (Calmodulin) Can I follow the same protocol which is followed for a protein? Probably, as long as the force field you choose can deal with calcium (most can). The only special thing to do would be to couple the Ca2+ ions to the protein for the purpose of temperature coupling, which requires a custom index group. -Justin @Justin: Why would you couple the calcium ion with the protein and not with the water / solvent? Just being interested. Thought commonly one couples the ions with the solvent. Thomas -- 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: MD Simulation of Calcium with protein
Am 15.01.2013 14:17, schrieb gmx-users-requ...@gromacs.org: On 1/15/13 8:12 AM, Thomas Schlesier wrote: Am 15.01.2013 12:52, schriebgmx-users-requ...@gromacs.org: On 1/15/13 5:49 AM, Devika N T wrote: HI I would like to know the protocol to be followed for performing MD simulation for calcium with protein (Calmodulin) Can I follow the same protocol which is followed for a protein? Probably, as long as the force field you choose can deal with calcium (most can). The only special thing to do would be to couple the Ca2+ ions to the protein for the purpose of temperature coupling, which requires a custom index group. -Justin @Justin: Why would you couple the calcium ion with the protein and not with the water / solvent? Just being interested. Thought commonly one couples the ions with the solvent. If the ion is structurally bound to the protein, I think it makes more sense to couple it to the protein rather than the solvent since the dynamics of the two species are more tightly linked. It's not a free-floating ion like others that might be in solution. -Justin This makes sense. Thanks for the answer. Thomas -- 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] Question about 'mdrun -nosum'
Dear all, i have a small question regarding the '-nosum' option of 'mdrun'. The manual states: For a global thermostat and/or barostat the temperature and/or pressure will also only be updated every nstlist steps. With this option the energy file will not contain averages and fluctuations over all integration steps. Second sentence is clear to me. But the first sentence give me some thoughts. I think this would introduce some errors, but what is about the magnitude of these? I would expected that the errors are in the same order of the errors in the forces due to the neighborlist search (error due to the fact that one uses 'nstlist=5' instead of 'nstlist=1') and would be more or less negilegible, if one doesn't use a very larger 'nstlist'-value. But to be on the save side i wanted to ask. Greetings Thomas -- 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: Pulling ion - US
I would also use the same residue from the pulling for the US. One thing you should be aware of is the pulling dimension: Now you have the pull-code only ativated for the z-direction. If you use this still in the US the ion can move freely in the xy-plane (freely in the sense of what is possible from the surrounding). One extreme example: The ion bounds to the protein and somehow (don't ask, think this wont happen in reality) and diffses away (in the xy-plane) after a long time it's so far away from the protein that the are no interaction with the protein and the ion interacts only with the surrounding water. Now you don't measure with the US potential the interaction of the ion with the protein, but the free diffusion of the ion. This case wont happen, since the probability that the ion unbounds itself from the protein goes down to the cellar. But i hope you get the idea what the gerneral problem is. If you the pull-dim in each direction this problem wouldn't occur, since the movement of the ion is also restrained in the xy-plane. Am 10.12.2012 21:40, schrieb gmx-users-requ...@gromacs.org: Would you also specify in each US window specific residue instead of the whole protein? Sreven On Mon, Dec 10, 2012 at 2:47 PM, Steven Neumanns.neuman...@gmail.com wrote: On Mon, Dec 10, 2012 at 2:11 PM, Justin Lemkuljalem...@vt.edu wrote: On 12/10/12 9:01 AM, Steven Neumann wrote: Dear Gmx Users, I am pulling away cation from the protein glutamic acid residue with: pull= umbrella pull_geometry = distance ; simple distance increase pull_dim= N N Y pull_start = yes ; define initial COM distance 0 pull_ngroups= 1 pull_group0 = Protein pull_group1 = NA pull_rate1 = 0.01 pull_k1 = 500 ; kJ mol^-1 nm^-2 I tried different pulling rates and simulation time to pull it 3 nm away. I tried pull rate of 0.1; 0.01 and 0.001. The interaction is so strong that the force reaches 600 kJ/mol/nm2 and they do not become separated - with position restraints protein looses its secondary structure and is draged by the ion - they do not become separated. Would you suggest constant force pulling in this case? Then I will extract initial coordinates for US windows. Can I use then US with harmonic potential in windows then and WHAM? You can generate coordinates in any way you wish. I would think that, regardless of the pull method, setting pull_group0 to the actual residue to which the ion is coordinated would be significantly more effective than pulling with respect to the entire protein, though it seems rather strange that the dissociation of an ion would cause a protein to unfold. A stronger force constant in pull_k1 may also help. -Justin Thank you Justin. That indeed helped. Steven -- 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] help about opls-aa for thiophene
Have a look there: http://virtualchemistry.org/molecules/110-02-1/index.php virtualchemistry.org is a really nice site (from David van der Spoel, and others i think), which has many paramters for solvents for the GAFF and OPLS force field. And also Physical properties for these. Greetings Thomas C4H4S. Thiophene is common compound . it seems oplss-aa does not have the parameters for it (such as dihedral angle). Any expert of opls-aa forcefield can help with ff parameters for thiophene? Thanks very much! Steven -- 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: help about opls-aa for thiophene
They have also the complete force field parameters under: OPLS Auxiliary topologyoplsaaff.itp if there are more paramters then in the oplsaa.ff directiory, then they have probably developed these parameters. They give this paper as a reference for the calculations. Probably something is mentioned there? Carl Caleman, Paul J. van Maaren, Minyan Hong, Jochen S. Hub, Luciano T. Costa and David van der Spoel Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant, J. Chem. Theor. Comput., 8, 61-74 (2012).DOI Greetings Thomas PS: If you reply with the whole digest, please delete all the stuff from the other questions. This makes it far easier to follow one question. Am 06.12.2012 18:21, schrieb gmx-users-requ...@gromacs.org: Dear Tomas, Thanks for your information! I look at the website you mentioned: http://virtualchemistry.org/molecules/110-02-1/index.php The *top file is available and the atom opls-aa types are assigned on *top file. But I can not find these parameters of dihedral angles from the default gmx (VERSION 4.5.5) directory of oplsaa.ff or the opls-aa (2005). For example: the parameters for dihedral of s-cw-cs-cs is missing Do you know where i can find these parameters? what versioni of opls-aa they are using? Thanks, Steven -- 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: force vs time plot
This should be in the pullf.xvg (time and then the forces). Am 24.11.2012 19:55, schrieb gmx-users-requ...@gromacs.org: Hi to all gromacs users, I am trying to run an umbrella sampling and i am getting the initial conformations by pulling simulations but i want to check the simulation through the force vs time plot to see if my complex did or did not separate, so how can i get this plot? Thank you for your time Paula -- 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: Bonds - force constant for given Beads
A good start might be: Phys. Chem. Chem. Phys., 2011, 13, 10437–10448 This paper is about hybrid-models (mixing CG and AA). But they discuss 'boltzmann inversion' and 'force matching', which are both methods to obtain CG-potentials. Since they use small molecules it focusses on nonbonded interactions, but one can probably transfer the methods to bonded interactions. Greetings Thomas Am 26.11.2012 18:45, schrieb gmx-users-requ...@gromacs.org: the distance distribution should be given by the Boltzmann factor of the potential energy function between the beads, assigning V(r)=0 for the most probable distance in your histogram. that's what you get when you take a molecule in vacuum and for instance you compare the dihedral distribution with the dihedral potential, the distribution is just exp[-U(theta)/RT], except maybe for an additive constant. you should be aware that the distribution may change appreciably depending on the environment, so this approach may be tricky: you may include implicit solvent effects on your bonded parameters and then you would end up with a forcefield that counts twice the solvent effects on the internal structure if you add explicit solvent in your model system. I hope it helps. Andre On Mon, Nov 26, 2012 at 2:06 PM, Steven Neumanns.neuman...@gmail.comwrote: Dear Gmx Users, I am planning to build coarse grained model based on the all atom simulation. I created (using VMD) beads representing 2-4 atoms of my protein chain. I want to extract bonded parameters. The equilibrium lenght for bonds (between specified beads) would be the average over the equilibrium from all atom simulation using g_dist between Centre of Mass of group of atoms belonging to given bead. My question: How can I extract the force constant for the bonds from all atom simulation between those beads? Thank you, Steven -- 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: Strange form of RDF curve
The rdf should not depend on the choice 'which is group A and which is group B'! Not if the system is well equilibrated and not if you consider only a single snapshot (in that case the rdf looks like garbage if the system is not really huge, but the RDF(A-B) and RDF(B-A) must be the same). I mean the the RDF is the radial distribution function. A distribution about distances, and in distances one has no information about direction, or from which point one goes to another. In an A and B system one has N_A * N_B distances, no matter if A or B is the reference! Going to the manual the given equation for the RDF is g(r) = 1/rho_B 1/N_A sum_A sum_B delta(r_AB - r) / (4pi r^2) Since after the sums it is clear that there is no difference between A and B one could argue that 1/rho_B 1/N_A is not equal to 1/rho_A 1/N_B looking into the picture which defines rho one recognises that rho=N/V - g_AB(r) = g_BA(r) Am 20.11.2012 12:00, schrieb gmx-users-requ...@gromacs.org: the opening statement from the statistical mechanical standpoint means a thorough sampling has been attained, so I do agree that the reference group choice would matter for an every day situation. Thoroughness is not attainable for most complex systems we might be interested in. On Tue, Nov 20, 2012 at 1:21 AM, Justin Lemkuljalem...@vt.edu wrote: On 11/19/12 10:04 PM, Andr? Farias de Moura wrote: from the statistical thermodynamics standpoint, rdf must be the same for both choices of reference group, i.e.,solute-solvent and solvent-solute must yield exactly the same rdf, the only difference being expected for the cumulative numbers, which depend on the particles number density. The reason why rdf's must be the same is the fact that rdf are connected to the pmf between the particles, and the pmf do not depend upon the choice of reference group, just on the reaction coordinate connecting the groups. Doesn't this explanation assume that the system is homogeneous, or at least, that the solute (reference group) of interest exhaustively samples configurations within the solvent? It's not intuitive to me why one would expect an Arg-water RDF to be the same as a Water-Arg RDF when the Arg residue is only part of a 582-residue protein, which means a considerable portion of the simulation unit cell contains neither the Arg of interest nor water, a fact that significantly impacts the binning as shown in the manual. The actual g_rdf command is essential in this regard, as well, depending on how the RDF is being calculated. One can produce wildly different results by changing only a a single element of the command. -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] Force vs distance plot in pulling simulation?
But be aware that the force depends on the pulling velocity. If you perform the simulation with two different pulling velocities you'll get two different forces for each distance. The easiest way to get the force as a function of the distance (without the bias of the pulling velocity) would be thermodynamik intergration. This is similar to the umbrella sampling. You set up some windows with different distances. But instead of restraining the distance with an umbrella potential (as in umbrella sampling) you use a constraint to fix the distance and than can determine the constraint force. Greetings Thomas Am 17.11.2012 16:44, schrieb gmx-users-requ...@gromacs.org: Hi Erik and Justin Thanks, writing such a script is easy. The point of it all would be to be able to map the magnitude of the pulling force to what I see happen in the pulling simulation. How else would you get an understanding of what the pulling force means? Thanks /PK The only solution is to write a simple script that parses out the columns you want. -Justin I don't see the point though. Except for checking implementation of the pull code. -- 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: Question about scaling
Am 13.11.2012 06:16, schrieb gmx-users-requ...@gromacs.org: Dear all, i did some scaling tests for a cluster and i'm a little bit clueless about the results. So first the setup: Cluster: Saxonid 6100, Opteron 6272 16C 2.100GHz, Infiniband QDR GROMACS version: 4.0.7 and 4.5.5 Compiler: GCC 4.7.0 MPI: Intel MPI 4.0.3.008 FFT-library: ACML 5.1.0 fma4 System: 895 spce water molecules this is a somewhat small system I would say. Simulation time: 750 ps (0.002 fs timestep) Cut-off: 1.0 nm but with long-range correction ( DispCorr = EnerPres ; PME (standard settings) - but in each case no extra CPU solely for PME) V-rescale thermostat and Parrinello-Rahman barostat I get the following timings (seconds), whereas is calculated as the time which would be needed for 1 CPU (so if a job on 2 CPUs took X s the time would be 2 * X s). These timings were taken from the *.log file, at the end of the 'real cycle and time accounting' - section. Timings: gmx-version 1cpu2cpu4cpu 4.0.7 422333843540 4.5.5 378032552878 Do you mean CPUs or CPU cores? Are you using the IB network or are you running single-node? Meant number of cores and all cores are on the same node. I'm a little bit clueless about the results. I always thought, that if i have a non-interacting system and double the amount of CPUs, i You do use PME, which means a global interaction of all charges. would get a simulation which takes only half the time (so the times as defined above would be equal). If the system does have interactions, i would lose some performance due to communication. Due to node imbalance there could be a further loss of performance. Keeping this in mind, i can only explain the timings for version 4.0.7 2cpu - 4cpu (2cpu a little bit faster, since going to 4cpu leads to more communication - loss of performance). All the other timings, especially that 1cpu takes in each case longer than the other cases, i do not understand. Probalby the system is too small and / or the simulation time is too short for a scaling test. But i would assume that the amount of time to setup the simulation would be equal for all three cases of one GROMACS-version. Only other explaination, which comes to my mind, would be that something went wrong during the installation of the programs? You might want to take a closer look at the timings in the md.log output files, this will give you a clue where the bottleneck is, and also tell you about the communication-computation ratio. Best, Carsten Please, can somebody enlighten me? Here are the timings from the log-file: #cores: 1 2 4 (all cores are on the same node) Computing: Domain decomp. 41.747.8up DD comm. load 0.0 0.0 - Comm. coord. 17.830.5up Neighbor search614.1 355.4 323.7 down Force 2401.6 1968.7 1676.0 down Wait + Comm. F 15.131.4up PME mesh 596.3 710.4 639.1 - Write traj.1.2 0.8 0.6 down Update 49.744.037.6down Constraints79.370.460.0down Comm. energies 3.2 5.3 up Rest 38.327.125.4down Total 3780.5 3254.6 2877.5 down PME redist. X/F133.0 120.5 down PME spread/gather 511.3 465.7 396.8 down PME 3D-FFT 59.488.9102.2 up PME solve 25.222.218.9down The two calculations-parts for which the most time is saved for going parallel are: 1) Forces 2) Neighbor search (ok, going from 2cores to 4cores does not make a big differences, but from 1core to 2 or 4 saves much time) Is there any good explains for this time saving? I would have thought that the system has a set number of interaction and one has to calculate all these interactions. If i divide the set in 2 or 4 smaller sets, the number of interactions shouldn't change and so the calculation time shouldn't change? Or is something fancy in the algorithm, which reducces the time spent for calling up the arrays if the calculation is for a smaller set of interactions? -- 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: Question about scaling
Sorry for reposting, but forgot one comment and added it now below: Am 13.11.2012 06:16, schrieb gmx-users-request at gromacs.org: Dear all, i did some scaling tests for a cluster and i'm a little bit clueless about the results. So first the setup: Cluster: Saxonid 6100, Opteron 6272 16C 2.100GHz, Infiniband QDR GROMACS version: 4.0.7 and 4.5.5 Compiler: GCC 4.7.0 MPI: Intel MPI 4.0.3.008 FFT-library: ACML 5.1.0 fma4 System: 895 spce water molecules this is a somewhat small system I would say. Simulation time: 750 ps (0.002 fs timestep) Cut-off: 1.0 nm but with long-range correction ( DispCorr = EnerPres ; PME (standard settings) - but in each case no extra CPU solely for PME) V-rescale thermostat and Parrinello-Rahman barostat I get the following timings (seconds), whereas is calculated as the time which would be needed for 1 CPU (so if a job on 2 CPUs took X s the time would be 2 * X s). These timings were taken from the *.log file, at the end of the 'real cycle and time accounting' - section. Timings: gmx-version1cpu2cpu4cpu 4.0.7 422333843540 4.5.5 378032552878 Do you mean CPUs or CPU cores? Are you using the IB network or are you running single-node? Meant number of cores and all cores are on the same node. I'm a little bit clueless about the results. I always thought, that if i have a non-interacting system and double the amount of CPUs, i You do use PME, which means a global interaction of all charges. would get a simulation which takes only half the time (so the times as defined above would be equal). If the system does have interactions, i would lose some performance due to communication. Due to node imbalance there could be a further loss of performance. Keeping this in mind, i can only explain the timings for version 4.0.7 2cpu - 4cpu (2cpu a little bit faster, since going to 4cpu leads to more communication - loss of performance). All the other timings, especially that 1cpu takes in each case longer than the other cases, i do not understand. Probalby the system is too small and / or the simulation time is too short for a scaling test. But i would assume that the amount of time to setup the simulation would be equal for all three cases of one GROMACS-version. Only other explaination, which comes to my mind, would be that something went wrong during the installation of the programs? You might want to take a closer look at the timings in the md.log output files, this will give you a clue where the bottleneck is, and also tell you about the communication-computation ratio. Best, Carsten Please, can somebody enlighten me? Here are the timings from the log-file (for GMX 4.5.5): #cores: 1 2 4 (all cores are on the same node) Computing: Domain decomp. 41.747.8up DD comm. load 0.0 0.0 - Comm. coord. 17.830.5up Neighbor search614.1 355.4 323.7 down Force 2401.6 1968.7 1676.0 down Wait + Comm. F 15.131.4up PME mesh 596.3 710.4 639.1 - Write traj.1.2 0.8 0.6 down Update 49.744.037.6down Constraints79.370.460.0down Comm. energies 3.2 5.3 up Rest 38.327.125.4down Total 3780.5 3254.6 2877.5 down PME redist. X/F133.0 120.5 down PME spread/gather 511.3 465.7 396.8 down PME 3D-FFT 59.488.9102.2 up PME solve 25.222.218.9down The two calculations-parts for which the most time is saved for going parallel are: 1) Forces 2) Neighbor search (ok, going from 2cores to 4cores does not make a big differences, but from 1core to 2 or 4 saves much time) For GMX 4.0.7 ist looks similar, whereas the difference between 2 and 4 cores is not so high as for GMX 4.5.5 Is there any good explains for this time saving? I would have thought that the system has a set number of interaction and one has to calculate all these interactions. If i divide the set in 2 or 4 smaller sets, the number of interactions shouldn't change and so the calculation time shouldn't change? Or is something fancy in the algorithm, which reducces the time spent for calling up the arrays if the calculation is for a smaller set of interactions? -- 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] Question about scaling
Dear all, i did some scaling tests for a cluster and i'm a little bit clueless about the results. So first the setup: Cluster: Saxonid 6100, Opteron 6272 16C 2.100GHz, Infiniband QDR GROMACS version: 4.0.7 and 4.5.5 Compiler: GCC 4.7.0 MPI: Intel MPI 4.0.3.008 FFT-library: ACML 5.1.0 fma4 System: 895 spce water molecules Simulation time: 750 ps (0.002 fs timestep) Cut-off: 1.0 nm but with long-range correction ( DispCorr = EnerPres ; PME (standard settings) - but in each case no extra CPU solely for PME) V-rescale thermostat and Parrinello-Rahman barostat I get the following timings (seconds), whereas is calculated as the time which would be needed for 1 CPU (so if a job on 2 CPUs took X s the time would be 2 * X s). These timings were taken from the *.log file, at the end of the 'real cycle and time accounting' - section. Timings: gmx-version 1cpu2cpu4cpu 4.0.7 422333843540 4.5.5 378032552878 I'm a little bit clueless about the results. I always thought, that if i have a non-interacting system and double the amount of CPUs, i would get a simulation which takes only half the time (so the times as defined above would be equal). If the system does have interactions, i would lose some performance due to communication. Due to node imbalance there could be a further loss of performance. Keeping this in mind, i can only explain the timings for version 4.0.7 2cpu - 4cpu (2cpu a little bit faster, since going to 4cpu leads to more communication - loss of performance). All the other timings, especially that 1cpu takes in each case longer than the other cases, i do not understand. Probalby the system is too small and / or the simulation time is too short for a scaling test. But i would assume that the amount of time to setup the simulation would be equal for all three cases of one GROMACS-version. Only other explaination, which comes to my mind, would be that something went wrong during the installation of the programs... Please, can somebody enlighten me? Greetings Thomas -- 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] GROMACS with different gcc and FFT versions but one unique *tpr file
Dear all, i have access to a cluster on which GROMACS is compiled with a different version of GCC and a different FFT libary (compared to the local machine). Will this affect simulationns if i prepare the *.tpr on the local machine and run the simulation on the cluster and the local machine? Sorry if this is a dumb question. I could imagine that the two simulations will be not numerical identical due to the different FFT libaries, but how strong this effect is and what else could happen i have no idea... Greetings Thomas -- 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=constraint gives zero forces
See the SHAKE algorithm in the manual. Especially equation (3.98) - G_i = \sum_k \lambda_k * (\partial \sigma_k)/(\partial r_i) where G_i is the constraint force \sigma_k are the equations for the constraints and \lambda_k is the lagrange multiplier 'Understanding molecular simulation' (D. Frenkel, B. Smit) gives these equations for bond-constraints as \sigma = r_ij^2 - d_ij^2 where r_ij is the real distance, and d_ij to one to which one wants to constrain the bond length from this \lambda_k would correspond to a force constant. Greetings Thomas Am 10.10.2012 18:41, schrieb gmx-users-requ...@gromacs.org: How can there be forces for holonomic constraints? Is this described by an equation in the manual? Just because there are values in the pullf.xvg file does not mean that these values are forces. If they are forces, what is the force constant and what is the equation that defines this force? Thank you, Chris. -- original message -- But for GMX 4.0.7 there are forces in the pullf.xvg. The forces which arise rom the contraint the hold the two groups fixed. I use them for thermodynamic integration... I use the following mdp-parameters, probably this gives you an idea what you might make different: ; AFM OPTIONS pull= constraint pull_geometry = distance pull_dim= Y Y Y pull_start = yes pull_nstfout= 10 pull_nstxout= 50 pull_ngroups= 1 pull_group0 = REF pull_group1 = ZUG pull_rate1 = 0.000 pull_init1 = 0.000 pull_constr_tol = 1e-06 But as Chris said, better post the complete mdp-file (probably for both GMX 4.0.7 and 4.5.x) so we can comment on them. Grettings Thomas -- 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=constraint gives zero forces
But for GMX 4.0.7 there are forces in the pullf.xvg. The forces which arise rom the contraint the hold the two groups fixed. I use them for thermodynamic integration... I use the following mdp-parameters, probably this gives you an idea what you might make different: ; AFM OPTIONS pull= constraint pull_geometry = distance pull_dim= Y Y Y pull_start = yes pull_nstfout= 10 pull_nstxout= 50 pull_ngroups= 1 pull_group0 = REF pull_group1 = ZUG pull_rate1 = 0.000 pull_init1 = 0.000 pull_constr_tol = 1e-06 But as Chris said, better post the complete mdp-file (probably for both GMX 4.0.7 and 4.5.x) so we can comment on them. Grettings Thomas Am 09.10.2012 04:00, schrieb gmx-users-requ...@gromacs.org: Please post your entire .mdp file and a snip of the output in your pullf and pullc files. (Your initial post on this topic was also missing these, although the text reads as if you intended to include them). I'll note that there are no forces when using constraints, so the fact that you get zero forces for a constrained run is not really surprising. I guess the thing is that it doesn't work to keep atoms in place for you, which we can help you with if you post more details. Chris. -- original message -- Following up on this post. I've tried the same runs using version 4.0.7, which gave immediate segmentation faults. Not sure if this is a clue or a trivial consequence of switching versions, but there it is. Any other ideas why the pullf output just contains zeros? Cheers, Alex -- 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=constraint gives zero forces
For me it seems that the problem with the constraints iteration is only relevant if both groups are connected with contraints over the rest of the system. It seems that your pull-code parameters are essentially the same like mine (apart from the fact that you have more groups). In my case the system consisted of a dimer, and the pull-constraint acted between both monomers (which had internal constraints) and i had never any problems when the seperation between the two molecules was reasonable. But nevertheless if the system breaks if you exchange the constraints to normal bonds, the is somewhere other a (probably additional) problem. Is the system stable, when you do not use the pull-code? Am 09.10.2012 14:26, schrieb gmx-users-requ...@gromacs.org: Hi Erik, Yes - I had a feeling I'd read that in the manual somewhere, which was why I tried replacing such constraints with bonds. As I wrote above, the LINCS warnings disappeared but the segfault remained. Alex 2012/10/9 Erik Marklund [via GROMACS]ml-node+s5086n5001821...@n6.nabble.com: Hi, Do you know there are issues with using pull=constraint on molecules that have constrained bonds? It's mentioned in the manual somewhere. Erik 9 okt 2012 kl. 11.39 skrev alex.bjorling: Sorry - forgot to mention that before crashing, the run with all other constraints removed produces a single line of pullf output: 0. -812.401-4002.84482.04 1951.47 138.953 -1806.55 -601.0072644.79 447.018 1768.6 -214.64 -199.829-2746.97 1177.7 476.39 288.535 -559.274123.08 114.493 851.86 550.558 As Thomas Schlesier mentions here, http://gromacs.5086.n6.nabble.com/pull-constraint-gives-zero-forces-tp5001817.html, the pullf output apparently contains the forces necessary to enforce the constraints. / Alex alex.bjorling wrote Thanks guys, Fixing the bug, recompiling and trying again results in a segfault like with version 4.0.7. I interpret this as GROMACS working fine, and suppose that there's something wrong with my input. Will continue this thread here and would really appreciate any ideas on how to proceed. Before the segfault, I get a bunch of LINCS warnings, all concerning atoms that have other constraints in the topology. Manually replacing these by stiff bonds in the itp gets rid of the LINCS warnings, but still produces an immediate segfault. The complete mdp follows. (Chris: previously posted via the web forum - it seems then you have to click the nabble link to see it). Cheers, Alex 50_constr3.mdp: ** pull= constraint pull_geometry = position pull_dim= Y Y Y pull_start = yes pull_group0 = BB pull_nstxout= 1000 pull_nstfout= 1000 pull_ngroups= 7 pull_constr_tol = 1e-6 pull_group1 = group1 pull_init1 = 0.0 0.0 0.0 pull_rate1 = 0.0 0.0 0.0 pull_group2 = group2 pull_init2 = 0.0 0.0 0.0 pull_rate2 = 0.0 0.0 0.0 pull_group3 = group3 pull_init3 = 0.0 0.0 0.0 pull_rate3 = 0.0 0.0 0.0 pull_group4 = group4 pull_init4 = 0.0 0.0 0.0 pull_rate4 = 0.0 0.0 0.0 pull_group5 = group5 pull_init5 = 0.0 0.0 0.0 pull_rate5 = 0.0 0.0 0.0 pull_group6 = group6 pull_init6 = 0.0 0.0 0.0 pull_rate6 = 0.0 0.0 0.0 pull_group7 = group7 pull_init7 = 0.0 0.0 0.0 pull_rate7 = 0.0 0.0 0.0 ; ; STANDARD MD INPUT OPTIONS FOR MARTINI 2.0 or 2.1 ; ; TIMESTEP IN MARTINI ; Most simulations are numerically stable ; with dt=40 fs, some (especially rings) require 20-30 fs. ; Note that time steps of 40 fs and larger may create local heating or ; cooling in your system. Although the use of a heat bath will globally ; remove this effect, it is advised to check consistency of ; your results for somewhat smaller time steps in the range 20-30 fs. ; Time steps exceeding 40 fs should not be used; time steps smaller ; than 20 fs are also not required. ;define = -DPOSRES integrator = md tinit= 0.0 dt = 0.02 nsteps = 250 ; 50 ns nstcomm = 1 comm-grps = nstxout = 0 nstvout = 0 nstfout = 0 nstlog = 1000 nstenergy= 100 nstxtcout= 1000 xtc_precision= 100 xtc-grps = energygrps = Protein ; NEIGHBOURLIST and MARTINI ; Due to the use of shifted potentials, the noise generated ; from particles leaving/entering the neighbour list is not so large, ; even when large time steps are being used. In practice, once every ; ten steps works fine with a neighborlist cutoff that is equal to the ; non-bonded cutoff (1.2 nm). However, to improve energy conservation ; or to avoid local heating/cooling, you may increase the update frequency ; and/or enlarge the neighbourlist cut-off (to 1.4 nm). The latter option ; is computationally
[gmx-users] RE: Re: Binding Energy to Binding affinity (Kd)
You could use 'constraint' pull-mode instead of the 'umbrella' mode. Than the distance would change gradually and you won't observe the fluctuations in the distance. greetings thomas Am 04.10.2012 16:58, schrieb gmx-users-requ...@gromacs.org: On 10/4/12 10:52 AM, jiang wrote: Justin Lemkul wrote On 10/2/12 4:39 AM, Du Jiangfeng (BIOCH) wrote: Hi Justin, I used ~20 windows to sample ~2 nm pulling. I notice that the distance between the complex being increased during the pulling but not gradually. At the distance of 0-1nm, there are 70 snapshots (the distance sometime increased sometimes decreased). At the distance of 1-2nm, there are only 30 snapshots (the distance kept increasing always). At the distance more than 3nm, the distance increased as 0.3nm of each snapshot, is it normal and reliable? I will assume you are using a harmonic potential (umbrella) to do the pulling. In this case, your observations are totally normal. When two species interact strongly, it is harder to pull them apart, thus the spring extends further to induce a larger force before displacement occurs. As the restoring forces are overcome, it is easier to move the pulled group through solution, so it makes more steady progress as the molecules are separated. Hi Justin, it is right I am using umbrella pulling. Now here is another hurdle in front of me: How to select the snapshots for umbrella samples? Since the distance between two groups went higher or lower at the beginning of the pulling. For example, during the pulling simulation, the distance changes like: 0.46 0.42 0.46 0.43 0.44 0.42 0.45 0.44 0.43 0.45 0.44 0.45 0.43 0.44 0.44 0.54 0.52 0.63 0.65 0.72 0.8 0.92 1.2 1.5 (nm) . I suppose it doesn't matter which snapshots to be chosen, as long as the snapshots can indicate a good spacing, the PMF result always should be same, right? You need reasonable spacing and sufficient sampling in each window to allow for proper overlap of the umbrella potentials. -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] what's the difference between gen_seed and ld_seed?
Am 23.08.2012 23:11, schrieb gmx-users-requ...@gromacs.org: hello : I am a little bit confused about the difference between gen_seed and ld_seed. I checked the manual, it is said: gen_seed used to initialize random generator for random velocities, when gen_seed is set to -1, the seed is calculated from the process ID number. This is often used coupled with gen_vel which is Generate velocities in grompp according to a Maxwell distribution at temperaturegen_temp [K], with random seed gen_seed. This is only meaningful with integratormd. As indicated in Jonhn E.Kerrigan's tutorial, gen_seed=-1 is always turned on in NVT, NPT and MD production step. However, in Justin's Lysozyme in Water tutorial, this option is only turned on in NVT and the following NPT and MD production were turned off. Yes because you undo your equilibration when you reassign initial velocities using gen_seed in NPT and MD production (NPT) (md/md-vv integrator). I am just wondering which option would be better for our simulations? How about ld_seed? here is the statement from manual: ld_seed: (1993) [integer] used to initialize random generator for thermal noise for stochastic and Brownian dynamics. When ld_seed is set to -1, the seed is calculated from the process ID. When running BD or SD on multiple processors, each processor uses a seed equal to ld_seed plus the processor number. when should we turn this option on? As stated, you use this option for use with the bd or sd integrators. Also the 'v-rescale' thermostat uses the ld_seed. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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] deltaG from PMF
As a first step, i would shift all curves so, that the energy of the minium is for all plots at the same (aribarity) value. The minimum should be the point which has sampled the best. If you shift then all values, it should be easier to spot differences between the plots. And probably make the analysis for every 10ns slices 10-20, 20-30, ... then it's easier to see if you have a drift or fluctuations. If you identify problematic regions you could do there additional umbrella simulations, probably with higher force constants, in order to sample that region better. In theory the PMF should converge if you wait long enough, till the system equilibrates under the external restraints (how long this takes? nobody knows). On could probably wait a long time, big question here is what is the biggest error you would tolerate. On the other hand one can't simulate forever... Look for what other people used for simulation-times for umbrella sampling (i have the impression that 50ns is rather long, but i could fool me here) and what they did to estimate if the calculations are converged. Then decide, what to do. Think nobody here would / will tell you this or this error is ok ... but if you do something reasonable, which is in accord what others do / did it should be ok. Am 24.08.2012 14:59, schrieb gmx-users-requ...@gromacs.org: I did g_wham calculation as you suggested for: 10-30 ns 30-50 ns 10-50 ns In some cases of PMF curves the drift is 0.3 kcal/mol but in some cases it is over 1.5 kcal/mol so huge drifts. In terms of huge drifts what would you suggest ? Steven On Tue, Aug 21, 2012 at 5:48 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 12:46 PM, Steven Neumann wrote: Thanks Thomas. Justin, could you please comment on this? I agree with everything Thomas has said. There's not really anything to say. -Justin Steven On Tue, Aug 21, 2012 at 5:37 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Am 21.08.2012 18:22, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or similar. But 'g_wham -h' prints some info about how to use the error analysis Greetings Thomas Steven Am 21.08.2012 17:25,schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ;
[gmx-users] Re: deltaG from PMF
Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Am 21.08.2012 17:25, schrieb gmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkul jalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 318 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed ; These options remove COM motion of the system ; Pull code pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 500 ; kJ mol^-1 nm^-2 pull_nstxout= 1000; every 2 ps pull_nstfout= 1000 ; every 2 ps Based on these settings you're allowing grompp to set the reference distance to whatever it finds in the coordinate file. It seems clear to me that the sampling indicates what I said before - you have an energy minimum somewhere other than where you started with. What that state corresponds to relative to what you think is going on is for you to decide based on the nature of your system. Whatever is occurring at 0.6 nm of COM separation is of particular interest, since the energy minimum is so distinct. So based on this the deltaG will correspond to -5.22 as the initial state was taken at 0.4 nm corresponding to 0 kcal/mol as the moment corresponding to the
[gmx-users] deltaG from PMF
Am 21.08.2012 18:22, schrieb gmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:49 PM, Thomas Schlesierschl...@uni-mainz.de wrote: Since your simulations of the individual windows are about 50 ns, you could first dismiss the first 10 ns for equilibration, and then perform the WHAM analysis for 10-30 ns and 30-50 ns. If everything is fine, you should see no drift. If you want to have more data for the analysis you could also use 5ns ; 5-27.5ns and 27.5-50ns. From the PMF it seems that the equilibrium state should be around 0.6 nm. To be sure, you can perform a normal simulation (without any restraints) from you initial starting window (~0.4nm) and a window near the minima (~0.6nm). Then after the equilibration phase, look at the distribution of the distance along the reaction coordinate. If in both cases the maximum is at ~0.6nm, this should be the 'true' equilibrium state of the system (instead of the first window of the PMF calculation) and i would measure \Delta G from this point. Greetings Thomas Thanks Thomas for this but finally I realised that my first configuration corresponds to 0.6 nm which is the minima so I take the free energy difference based on this value and plateau. I want also to calculate error bars. Would you do this: Final PMF curve for 10-50 ns Error bars from: g_wham -b 3 -e 4 g_wham -b 5 -e 6 Think this approach would be good to see if you have any drifts. But for error bars there is something implemented in 'g_wham'. But i never used it, since for my system umbrella sampling is not really applicable, only TI. So i can't comment on it, if there is anything one should be aware of, or similar. But 'g_wham -h' prints some info about how to use the error analysis Greetings Thomas Steven Am 21.08.2012 17:25, schriebgmx-users-requ...@gromacs.org: On Tue, Aug 21, 2012 at 4:21 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:18 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 4:13 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 11:09 AM, Steven Neumann wrote: On Tue, Aug 21, 2012 at 3:48 PM, Justin Lemkuljalem...@vt.edu wrote: On 8/21/12 10:42 AM, Steven Neumann wrote: Please see the example of the plot from US simulations and WHAM: http://speedy.sh/Ecr3A/PMF.JPG First grompp of frame 0 corresponds to 0.8 nm - this is what was shown by grompp at the end. The mdp file: ; Run parameters define = -DPOSRES_T integrator = md; leap-frog integrator nsteps = 2500 ; 100ns dt = 0.002 ; 2 fs tinit = 0 nstcomm = 10 ; Output control nstxout = 0 ; save coordinates every 100 ps nstvout = 0 ; save velocities every nstxtcout = 5; every 10 ps nstenergy = 1000 ; Bond parameters continuation= no ; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs vdwtype = Switch rvdw-switch = 1.0 rlist = 1.4 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.4 ; short-range electrostatic cutoff (in nm) rvdw= 1.2 ; short-range van der Waals cutoff (in nm) ewald_rtol = 1e-5 ; relative strenght of the Ewald-shifted potential rcoulomb ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.12 ; grid spacing for FFT fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 318 318 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 318 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed ; These options remove COM motion of the system ; Pull code pull
[gmx-users] Re: angle constraints
Thanks for the clarification. Forgot that virtual sites have no mass. With this info it is clear why a setup with 2 'normal' and one dummy atom would not work. Am 26.07.2012 01:42, schrieb gmx-users-requ...@gromacs.org: On 26/07/2012 4:12 AM, Broadbent, Richard wrote: Virtual sites are by definition have no mass. If you simply ignore the mass of the carbon the molecule will be too light and its translational momentum will therefore be too small meaning it will move too quickly. If you place half the mass of the carbon on each oxygen the moment of inertia will be wrong and the molecule will spin too slowly. All correct so far. In practice you have to decide what you want to loose or if a balance between the two is better. Not true, as illustrated by the link I gave earlier in the thread, which nobody seems to have read/understood. One needs at least two massive particles to describe the available degrees of freedom of a linear molecule, and using exactly two side-steps the angle constraint issue. Each must have half the total mass of CO2 and the distance between them is chosen to reproduce the moment of inertia. These will not be in suitable positions to have non-bonded interactions, of course. Then three (massless) virtual sites are constructed from those two, and these are the only ones that have the non-bonded interactions. Richard On 25/07/2012 14:44, Thomas Schlesierschl...@uni-mainz.de wrote: Ok, read the topic about the acetonitril. But i'm somewhat clueless: Why is the following setup wrong: Use 2 particles as normal atoms. Put the third as a dummy in between. Give each particle its 'normal' mass? I would assume that this system should have the right mass and moment of inertia, due to the fact the all individual masses and the positions one the particles would be correct. The virtual site so constructed cannot have mass, so this cannot be an accurate model. Only idea i have, why this setup could be flawed, would be that the third particle does only interact indirectly through the other two particles (i mean, virtual site interacts normally with all othe particles, but the force which would act on the dummy get redistributed to the other particles)... and then it's mass does not come into play, since it new position is determined only by the other two particles. so the complete molecule would move with a reduced mass?!? Still not an accurate model - you'd have a CO2 with three sites and mass only at two of the sites, so either the mass or moment of intertia must be wrong. Mark Can anyone comment on this? greetings thomas On 25/07/2012 10:08 PM, Thomas Schlesier wrote: What you have done there looks very strange... easiest wy would be: define the two oxygens as normal atoms (1,2), give them a bondlength twotimes the C-O bond length define the carbon as a dummy (3), while you construct it's position from both oxygens with a=0.5 one thing i don't know is how to handle the mass: 1) give both oxygen half of the system mass 2) give all atoms their normal mass would tend to (2) One should want to get both the total mass and the moment of inertia correct... http://lists.gromacs.org/pipermail/gmx-users/2003-September/007095.html. Mark greetings thomas Am 25.07.2012 13:15, schrieb gmx-users-request at gromacs.org: How to choose the positions of the dummy atoms while constraining the angle for a linear triatomic molecule? The topology for a such molecule ( af for example CO2 ) is as follows [ moleculetype ] ; Namenrexcl CO2 2 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeBchargeB massB ; residue 503 CO rtp CO q 0.0 1 D1503 CO D1 1 0 21.90158 ; qtot 0 2 D2503 CO D2 2 0 21.90158 ; qtot 0 3 CE503 CO CO 30.7 0.0 ; qtot 0.7 4 OE503 COOC1 4 -0.35 0.0 ; qtot 0.35 5 OE503 COOC2 5 -0.35 0.0 ; qtot 0.35 [ constraints ] ; ai aj funct b0 1 2 1 0.2000 [ dummies2 ] ; aiajak funct a 3 1 2 1 0.0170 4 1 2 1 0.1000 5 1 2 1 0.2170 [ exclusions ] 3 4 5 4 5 3 5 4 3 The .rtp file for CO2 [ CO ] [ atoms ] D1 D1 0. 1 D2 D2 0. 2 COCE 0.70003 OC1 OE -0.3500 4 OC2 OE -0.35005 [ bonds ] CO OC1 CO OC2 Can anyone please check above file parts whether I'm doing correct or not ? -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before
[gmx-users] Re: angle constraints
What you have done there looks very strange... easiest wy would be: define the two oxygens as normal atoms (1,2), give them a bondlength twotimes the C-O bond length define the carbon as a dummy (3), while you construct it's position from both oxygens with a=0.5 one thing i don't know is how to handle the mass: 1) give both oxygen half of the system mass 2) give all atoms their normal mass would tend to (2) greetings thomas Am 25.07.2012 13:15, schrieb gmx-users-requ...@gromacs.org: How to choose the positions of the dummy atoms while constraining the angle for a linear triatomic molecule? The topology for a such molecule ( af for example CO2 ) is as follows [ moleculetype ] ; Namenrexcl CO2 2 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeBchargeB massB ; residue 503 CO rtp CO q 0.0 1 D1503 CO D1 1 0 21.90158 ; qtot 0 2 D2503 CO D2 2 0 21.90158 ; qtot 0 3 CE503 CO CO 30.7 0.0 ; qtot 0.7 4 OE503 COOC1 4 -0.35 0.0 ; qtot 0.35 5 OE503 COOC2 5 -0.35 0.0 ; qtot 0.35 [ constraints ] ; ai aj funct b0 1 2 1 0.2000 [ dummies2 ] ; aiajak funct a 3 1 2 1 0.0170 4 1 2 1 0.1000 5 1 2 1 0.2170 [ exclusions ] 3 4 5 4 5 3 5 4 3 The .rtp file for CO2 [ CO ] [ atoms ] D1 D1 0. 1 D2 D2 0. 2 COCE 0.7000 3 OC1 OE -0.35004 OC2 OE -0.35005 [ bonds ] CO OC1 CO OC2 Can anyone please check above file parts whether I'm doing correct or not ? -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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] angle constraints
Ok, read the topic about the acetonitril. But i'm somewhat clueless: Why is the following setup wrong: Use 2 particles as normal atoms. Put the third as a dummy in between. Give each particle its 'normal' mass? I would assume that this system should have the right mass and moment of inertia, due to the fact the all individual masses and the positions one the particles would be correct. Only idea i have, why this setup could be flawed, would be that the third particle does only interact indirectly through the other two particles (i mean, virtual site interacts normally with all othe particles, but the force which would act on the dummy get redistributed to the other particles)... and then it's mass does not come into play, since it new position is determined only by the other two particles. so the complete molecule would move with a reduced mass?!? Can anyone comment on this? greetings thomas On 25/07/2012 10:08 PM, Thomas Schlesier wrote: What you have done there looks very strange... easiest wy would be: define the two oxygens as normal atoms (1,2), give them a bondlength twotimes the C-O bond length define the carbon as a dummy (3), while you construct it's position from both oxygens with a=0.5 one thing i don't know is how to handle the mass: 1) give both oxygen half of the system mass 2) give all atoms their normal mass would tend to (2) One should want to get both the total mass and the moment of inertia correct... http://lists.gromacs.org/pipermail/gmx-users/2003-September/007095.html. Mark greetings thomas Am 25.07.2012 13:15, schrieb gmx-users-request at gromacs.org: How to choose the positions of the dummy atoms while constraining the angle for a linear triatomic molecule? The topology for a such molecule ( af for example CO2 ) is as follows [ moleculetype ] ; Namenrexcl CO2 2 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeBchargeB massB ; residue 503 CO rtp CO q 0.0 1 D1503 CO D1 1 0 21.90158 ; qtot 0 2 D2503 CO D2 2 0 21.90158 ; qtot 0 3 CE503 CO CO 30.7 0.0 ; qtot 0.7 4 OE503 COOC1 4 -0.35 0.0 ; qtot 0.35 5 OE503 COOC2 5 -0.35 0.0 ; qtot 0.35 [ constraints ] ; ai aj funct b0 1 2 1 0.2000 [ dummies2 ] ; aiajak funct a 3 1 2 1 0.0170 4 1 2 1 0.1000 5 1 2 1 0.2170 [ exclusions ] 3 4 5 4 5 3 5 4 3 The .rtp file for CO2 [ CO ] [ atoms ] D1 D1 0. 1 D2 D2 0. 2 COCE 0.70003 OC1 OE -0.3500 4 OC2 OE -0.35005 [ bonds ] CO OC1 CO OC2 Can anyone please check above file parts whether I'm doing correct or not ? -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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: angle constraints
As others said: type 2 virtual site check chapter 4.7 and 5.2.2 in the manual (version 4.5.x). Greetings Thomas Am 24.07.2012 12:00, schrieb gmx-users-requ...@gromacs.org: Sorry I should mention it at the very beginning that I have a linear molecule and the angle is to be constrained at 180 degree. So what is the best way of constraining the angle for the linear molecule ? On Tue, Jul 24, 2012 at 2:06 PM, Mark Abrahammark.abra...@anu.edu.au wrote: On 24/07/2012 6:07 PM, Broadbent, Richard wrote: An Angle constraint amounts to fixing a triangle. To do this you need to constrain the distances between all the atoms as you know the of the bonds 6064, 6063 and 6063, 6065 and the angle between the two bonds it is a trivial geometry problem calculate the length of the third side of the triangle (6064,6065). However, as you are attempting to hold these atoms in a straight line I would suggest that a type 2 virtual site might (depending on your system) be a better idea. Indeed, a much better idea. Mark Richard On 24/07/2012 07:21, tarak karmakartarak20...@gmail.com wrote: Oh ! Thnaks I saw that table, the angle_restrain option is there but not constraints . Anyway if suppose, I fix the distance between the two terminal atoms of the molecule, the angle will eventually be fixed at a particular given value. Is that the logic ? Actually I searched for this problem so many times but didn't get proper clue; in one of those mails I saw someone has dealt with some dummy atoms. I could not able to digest that logic. On Tue, Jul 24, 2012 at 11:01 AM, Mark Abrahammark.abra...@anu.edu.au wrote: On 24/07/2012 3:21 PM, tarak karmakar wrote: Dear All, I am constraining one angle in my protein sample by incorporating [ constraints ] block in topology file as [ constraints ] ; index1 index2 index3 funct angle 6064 6063 6065 1 180.0 while doing that its showing the following error Program grompp, VERSION 4.5.5 Source code file: topdirs.c, line: 174 Fatal error: Invalid constraints type 6065 For more information and tips for troubleshooting, please check the GROMACS website athttp://www.gromacs.org/Documentation/Errors As you will see in table 5.6, this is not a valid option for [constraints] - you can only specify bonds. You will need to construct a triangle of bond constraints. Mark Then I rechecked the angle block, that specific angle is there in that angle section, part of it as follows 6039 6057 6058 1 6039 6057 6059 1 6058 6057 6059 1 6064 6063 6065 1 6067 6066 6068 1 6067 6066 6069 1 6068 6066 6069 1 6071 6070 6072 1 6071 6070 6073 1 6072 6070 6073 1 [ constraints ] ; index1 index2 index3 funct angle 6064 6063 6065 1 180.0 [ dihedrals ] ; aiajakal functc0c1 c2c3c4 -- gmx-users mailing listgmx-us...@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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 togmx-users-requ...@gromacs.org. * Can't post? Readhttp://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-us...@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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 togmx-users-requ...@gromacs.org. * Can't post? Readhttp://www.gromacs.org/Support/Mailing_Lists -- Tarak Karmakar Molecular Simulation Lab. Chemistry and Physics of Materials Unit Jawaharlal Nehru Centre for Advanced Scientific Research Jakkur P. O. Bangalore - 560 064 Karnataka, INDIA Ph. (lab) : +91-80-22082809 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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] Final state not reached in pulling simulation
It could be possible tht you do not pull into the 'right' direction. if there is another group between 'GTP' and 'Residue' you will get clashes and 'Residue' won't move further (could be a water molecule, or some other part of 'GTP'). If this happens you should observe an increase in the force due to the umbrella potential. If the problems are due to waters molecules which block the 'pathway' you could just delete them. If another group is in the way, you might want to change the pull-vector (and if lucky find the right one). But don't know what would be the best strategy in this case. Maybe you can look into what docking-people do, seems to me that your simulation is related to what they do (but myself has absolute no knowledge about docking simulations). greetings thomas Am 12.07.2012 17:26, schrieb gmx-users-requ...@gromacs.org: Hello all, I m performing a pulling simulation on my Protein-Mg-GTP complex. I have considered pulling between the GTP and a residue of protein. The pull code in the .mdp file im using is as follows: ; Pull code pull= umbrella pull_geometry = distance ; simple distance increase pull_dim= N N Y pull_start = yes ; define initial COM distance 0 pull_ngroups= 1 pull_group0 = GTP pull_group1 = Residue pull_rate1 = -0.5 ; 0.5 nm per ps = .05 nm per ns pull_k1 = 1 ; kJ mol^-1 nm^-2 The initial distance between GTP and the residue was 7 A and the desired one was 3A. After the completion of run (10ns), I could get a trajectory where the final distance was still 4.25 A. I tried to continue the simulation for another 10ns with the same value for pull_k1 parameter and one by increasing the value to 100,000 also. In both of the case, the trajectories showed the distance stabilized near _4.25 A only. Can anyone please tell me the reason behind it? What should I do, so that I could get the desired distance ? Any suggestion and help is welcome !!! Thanks, Neeru Sharma -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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] specifying the direction of Pull in US
If you want to pull along a vector which connects to groups, the easiest way is to run 'g_dist' over your starting *.gro file. this measures the distance and the vector connecting both groups. From GROMACS-4.0.x you don't need to normalise the vector. So you can directly use this vector. greetings thomas Thanks for ur suggestion Justin, I'm facing trouble in setting that vector, actually I cant figure out how can i set up a vector. Is there any easier way with which i can set up a vector. Thanks -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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] Nucleic acid simulation
Heard that RNA/DNA system could be a little trickier than proteins due to the many negative charges. I found somewhen a nice article about RNA simulations in general. Probably some questions you have / will have are answered there: A short guide for molecular dynamics simulations of RNA systems Yaser Hashem, Pascal Auffinger Methods 47 (2009) 187–197 greetings thomas Am 04.07.2012 16:52, schrieb gmx-users-requ...@gromacs.org: On 7/4/12 6:45 AM, Ravi Raja Merugu wrote: Hello every one, Im interesting in performing a MD for Protein - RNA complex , Can any one suggest a good tutorial. Such systems do not differ significantly from simulations of simple proteins in water, since pdb2gmx can produce topologies for both proteins and nucleic acids. The remaining workflow is basically the same. -Justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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] Pulling ligand - Different Profiles (Force vs time)
As a side note: The rupture process is a stochastic process, so a single rupture force is meaningless, since it is a distributed property. So you need to do many simulations to get the distribution / average rupture force. It that same like equilibrium properties, one doesn't determine them from only a single frame from a trajectory. Am 27.06.2012 15:52, schrieb gmx-users-requ...@gromacs.org: On 6/27/12 9:36 AM, Steven Neumann wrote: On Wed, Jun 27, 2012 at 1:51 PM, Justin A. Lemkuljalem...@vt.edu wrote: On 6/27/12 7:48 AM, Steven Neumann wrote: Dear Gmx Users, I obtained a protein-ligand complex from 100ns simulation. Now I am pulling my ligand away from the protein after the energy minimzation in water and equilibration of 100ps (two coupling baths: Protein, LIG_Water_and_ions). Then I proceed my pulling : grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr mdrun -s pull.tpr -deffnm pull title = Umbrella pulling simulation define = -DPOSRES ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 50; 1 ns nstcomm = 10 ; Output parameters nstxout = 0 nstvout = 0 nstfout = 500 nstxtcout = 1000 ; every 1 ps nstenergy = 500 ; Bond parameters constraint_algorithm= lincs constraints = all-bonds continuation= yes ; continuing from NPT ; Single-range cutoff scheme nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb= 1.4 rvdw= 1.2 vdwtype = Switch rvdw-switch = 1.0 ; PME electrostatics parameters coulombtype = PME fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein LIG_Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 298 298 ; reference temperature, one for each group, in K ; Pressure coupling is on Pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 2.0 compressibility = 4.5e-5 ref_p = 1.0 ; Generate velocities is off gen_vel = no ; Periodic boundary conditions are on in all directions pbc = xyz ; Long-range dispersion correction DispCorr= EnerPres ; Pull code pull= umbrella pull_geometry = distance ; simple distance increase pull_dim= N N Y pull_start = yes ; define initial COM distance 0 pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG pull_rate1 = 0.004 ; 0.004 nm per ps = 4 nm per ns pull_k1 = 500 ; kJ mol^-1 nm^-2 I run 3 pulling simulations with the same mdp and I obtain 3 different profiles (Force vs time). Then I used 2xlonger pulling with the same pulling distance and I run 3 simulations again. Each time I obtain different profile. Can anyone explain me this? I am using velocities from npt simulation as above (gen_vel = no and continuation = yes) so I presume the output should be similar. Please, advice. I assume you're passing a checkpoint file to grompp? If you're relying on velocities from the .gro file, they are of insufficient precision to guarantee proper continuation. Thank you Justin. I am using according to your tutorial: grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr mdrun -s pull.tpr -deffnm pull Would you suggest: grompp -f pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr mdrun -s pull.tpr -cpi npt.cpt -deffnm pull ?? No, I would not, especially if the NPT run uses position restraints, in which case the two phases are different. I missed the command line in the earlier message. What you are doing makes sense. Profiles do not vary slightly - the maximum pulling force (breaking point) varies from 290 to 500 kJ/mol nm2 which is really a lot. Consult the points below and watch your trajectories. If you're getting different forces, your ligands are experiencing different interactions. SMD is a path-dependent, non-equilibrium process. Good sampling and a justifiable path are key. -Justin Small variations are inherent in any simulation set, and in the case of pulling, small changes (though intentional) are the basis for Jarzynski's method. In any case, all MD simulations are chaotic and so it depends on what your definition of different is in the context of whether or not there are meaningful changes imparted through the course of each simulation. Also note that in the absence of the -reprod flag, the same .tpr file may result in a slightly
[gmx-users] Force constant - units
You have 1mol of your system. conversition factor for kJ/(mol*nm) - pN is approx 1.661 Am 29.06.2012 10:33, schrieb gmx-users-requ...@gromacs.org: Dear Gmx Users, How to recalculate the force constant from the harmonic potential: 1 [pN/A] into [kJ/mol nm2] ? Where is the [mol] here? Thanks, Steven -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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] PMF trails off to infinity.
think you encounter the problem, that you construct your pmf from a 3d simulation and project it onto 1d, but do no correction. For TI (if you constrain the distance in all three directions) the pmf is given by V_pmf(r) = - \int [ F_c + 2/(beta*r) ] dr with F_c the constraint force and \beta = 1/(k_B*T). for umbrella sampling one should need the same factor -2/beta \int 1/r dr if one restraints the system in 3 directions. Since 'g_wham' uses no *.tpr it can not know the value of 'pull_dim' one should introduce the factor afterwards. There could also the problem that the dummy-atom interacts with the ligand, or other clashes, but i don't think so, because the force looks too fine. If the ligand would crash into something one should see a greatly increasing force. like around 1000-2500 in your force-profile, when the protein still helds the ligand in the binding pocket. Additional comment: 'pull_geometry = distance' is a really bad idea for pmf generation if one has an actual distance of 0 and the system is not isotropic, since distance only knows distances and has no ideas about directrions: place a reference particle at the origin of a box put a second particle on top of it pull along x pull along -x in both cases you get the distance r=|x|=|-x| if the system is isotropic, everything is fine if system is anisotropic you get problems on the nice side, this problem should affect simulations where one pulls a ligand away from a binding pocket, only for the windows which is centered the nearest to 0. for mapping an ion-channel with a reference group in the middle of the channel it's (/ it should be) far more worse, (if the system is not isotropic). Somewhere in the archive is a longer explainition, i discussed the topic with 2-3 different people... Am 03.07.2012 11:34, schrieb gmx-users-requ...@gromacs.org: Basically what I'm doing is pulling a ligand out of a protein towards a dummy atom, which has a mass of 1 and no charge. I've attached the a portion .mdp files for both the smd portion and the umbrella sampling. I know that the ligand gets very close possibly crashing into the dummy atom. So from what you're saying, I'm thinking this might be the source of the problem. - Laura ## smd.mdp## ; Generate velocites is on at 300 K. gen_vel = yes gen_temp = 300.0 gen_seed = 123456 ;COM PULLING ; Pull type: no, umbrella, constraint or constant_force pull = umbrella ; Pull geometry: distance, direction, cylinder or position pull_geometry = distance ; Select components for the pull vector. default: Y Y Y pull_dim = Y Y Y ; Cylinder radius for dynamic reaction force groups (nm) pull_r1 = 1 ; Switch from r1 to r0 in case of dynamic reaction force pull_r0 = 1.5 pull_constr_tol = 1e-06 pull_start = yes pull_nstxout = 10 pull_nstfout = 1 ; Number of pull groups pull_ngroups = 1 ; Group name, weight (default all 1), vector, init, rate (nm/ps), kJ/(mol*nm^2) pull_group0 =JJJ pull_weights0 = pull_pbcatom0 = 0 pull_group1 = CPZ pull_weights1 = pull_pbcatom1 = 0 pull_vec1 = pull_init1 = pull_rate1 = -0.0006 pull_k1 = 800 pull_kB1 = ### um.mdp ## ;COM PULLING ; Pull type: no, umbrella, constraint or constant_force pull = umbrella ; Pull geometry: distance, direction, cylinder or position pull_geometry = distance ; Select components for the pull vector. default: Y Y Y pull_dim = Y Y Y ; Cylinder radius for dynamic reaction force groups (nm) pull_r1 = 1 ; Switch from r1 to r0 in case of dynamic reaction force pull_r0 = 1.5 pull_constr_tol = 1e-06 pull_start = yes pull_nstxout = 100 pull_nstfout = 100 ; Number of pull groups pull_ngroups = 1 ; Group name, weight (default all 1), vector, init, rate (nm/ps), kJ/(mol*nm^2) pull_group0 =JJJ pull_weights0 = pull_pbcatom0 = 0 pull_group1 = CPZ pull_weights1 = pull_pbcatom1 = 0 pull_vec1 = pull_init1 = 0 pull_rate1 = 0 pull_k1 = 1000 pull_kB1 = On 07/02/2012 04:38 PM, Justin A. Lemkul wrote: On 7/2/12 4:30 PM, Laura Kingsley wrote: Just for clarification, the PMF is read from right to left and the force profile is read from left to right. The dramatic change in the magnitude and sign of the force, coupled with the steady increase in PMF, indicates to me that some elements of your system are crashing into one another. In the absence of an accompanying explanation of what you're doing (description of system, procedure with .mdp parameters, etc) that's the best I can offer. -Justin On 07/02/2012 04:27 PM, Kingsley, Laura J wrote: Here is a link to both the PMF profile and the force profile: http://s1064.photobucket.com/albums/u370/laurakingsley/?action=viewcurrent=pull_fig.jpg -Original Message- From:gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] On Behalf Of Justin A. Lemkul Sent: Monday, July 02, 2012 3:56 PM To: Discussion list for GROMACS users Subject: Re: [gmx-users] PMF trails off to infinity. On 7/2/12 3:53 PM, Laura Kingsley
[gmx-users] Umbrella - Force constant
It also depends in some cases strongly on the system. I have a two-state system in which both states are rather narrow (doing a normal pulling simulation, the end-to-end-distance seems nearly constant). In these two regions one could use small force constants. but both state are seperated by a distance of more then 0.5 nm. So if one would sample the transition state (between both states) accurately one would need really high force constants... Was too lazy to determine the right force constant for every region, so used in the end TI (thermodynamic integration), which one could describe as umbrella sampling with constraints. Am 03.07.2012 16:49, schrieb gmx-users-requ...@gromacs.org: On Tue, Jul 3, 2012 at 2:04 PM, Justin A. Lemkuljalem...@vt.edu wrote: On 7/3/12 8:41 AM, Steven Neumann wrote: Dear Gmx Users, Do you know or can you suggest some results based on the comparison of the force constant in Umbrell Sampling? Any literature? That would be lovely, but I've never seen such a thing. One could probably write a book with all the test cases that would be required. My gut tells me that you can't generalize too much in terms of pulling simulations - the approach depends on what is being pulled (small molecule, peptide, large protein), what the medium is (water, membrane, etc), and what the interacting partner is (protein surface, ion channel, binding pocket). As far as I understand when you use the same staring coordinates (from the same pulling simulation) for windows but you just change the force constant (e.g. from 500 to 2000 kJ/mol nm2) you should increase number of windows (for f=2000) as smaller force constant will cover wider neigboruring distances - that makes sense. I am curious whether the final result will be the same? I guess with stronger force it will converge faster but more windows are required. is it the only one difference? Without a systematic comparison, it's hard to say, but in theory if one samples sufficiently and has good overlap between neighboring windows, the results should converge to the same answer. If someone knows of some applicable literature that has done such comparisons, please post a reference. I'd love to see it. Most SMD and US methodology is written with hand-waving explanations as to what the authors did and why it worked, and I have a suspicion that most reviewers don't have a better idea so they can't refute such claims. -Justin Thanks Justin. The only thing I found: Beno?t Roux, The calculation of the potential of mean force using computer simulations, Computer Physics Communications, Volume 91, Issues 1?3, 2 September 1995, Pages 275-282, ISSN 0010-4655, 10.1016/0010-4655(95)00053-I. (http://www.sciencedirect.com/science/article/pii/001046559500053I): Furthermore, the convergence properties of umbrella sampling calculations may be exploited more effectively using WHAM. Generating short umbrella sampling simulations for a large number of narrow windows is computationally more advantageous than generating longer simulations with a smaller number of wider windows (...) If Nw simulations are used to cover the whole range L, the force constant K of the umbrella sampling potential must be chosen to insure a proper overlap between the adjacent windows, i.e., each window should cover a range of dL = L/Nw and the value of K should be on the order of kBT/dL^2 based on the magnitude of the rms fluctuations. It follows that the total simulation time Trot needed to generate the Nw windows varies as ,~ L2/Nw. Thus, it is more advantageous to run short umbrella sampling simulations for a large number of narrow windows (the simulation time required to prepare and equilibrate the windows, which is also important, is ignored in this simple analysis) If anyone would find something, please post. that is an interesting gap in free energy calculations. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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: gmx-users Digest, Vol 99, Issue 12
ok, had still the 4.0.x version in mind, there the .tpr-files were not used. could be that the factor was introduced in 4.5.x, but i don't think so. If there would be steric crashes, which would account for the strong increase in the PMF, one should see them in the force profile. for more information on the factor in regard to TI: E. Paci, G. Ciccotti, M. Ferrario; Chem. Phys. Letters 176, 6: 581-587, 1991 On 7/3/12 2:50 PM, Thomas Schlesier wrote: think you encounter the problem, that you construct your pmf from a 3d simulation and project it onto 1d, but do no correction. For TI (if you constrain the distance in all three directions) the pmf is given by V_pmf(r) = - \int [ F_c + 2/(beta*r) ] dr with F_c the constraint force and \beta = 1/(k_B*T). for umbrella sampling one should need the same factor -2/beta \int 1/r dr if one restraints the system in 3 directions. Since 'g_wham' uses no *.tpr it can not know the value of 'pull_dim' one should introduce the factor afterwards. The input to g_wham (at least in modern Gromacs versions) is a list of .tpr files (passed to the -it flag) and a list of pullx or pullf files (with -ix or -if). How is g_wham ignorant of these pull settings? -Justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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] Residue 'DUM' not found in residue topology database
'pdb2gmx' is mainly a tool to get the topology for (bigger) molecules from a *.pdb file. 'pdbgmx' is expecting an entry in the *.rtp file for the dummy atoms. You could add one in that file... (see further below) But if the dummys only interact via nonbonded interactions it would be more convient to write the topology by hand. In the *.top file, you would need after the protein a new entry with [ moleculetype ] and [ atoms ] But you should also have this nonbonded interaction in the *nb.itp and file of your force field and the atom definition in the *.atp file. If there are also bonded parameters the new *.rtp entry would be better. For this look at the definition of an aminoacid, or a single molecule and write something similar for your dummy atoms. Essential you would need the atoms definitions and charges, and then an entry for which atoms is bonded to which. But you should also have these interactions (nonbonded and bonded) in the *nb.itp and *bon.itp files of your force field and the atom definition in the *.atp file. The basic informations for both ways are in chapter 5 of the manual. Hope this gives some hints. Greetings Thomas Schlesier Am 20.06.2012 16:18, schrieb gmx-users-requ...@gromacs.org: Hi everybody, I try to use GROMACS for my protein where I added a layer of dummy atoms simulating the membrane around it. But now when I want to call pdb2gmx I always get the error: Fatal error: Residue 'DUM' not found in residue topology database I understand the error, that there is no entry for the dummy residue in the .tpr file. But is there a possibility to add it there. Because the file don't look like that I can just write DUM in it. my pdb file looks like this: ATOM 2597 CB LEU A 313 12.816 -29.877 -23.547 1.00 55.31 C ATOM 2598 CG LEU A 313 14.261 -29.695 -23.059 1.00 54.06 C ATOM 2599 CD1 LEU A 313 15.169 -30.640 -23.866 1.00 53.26 C ATOM 2600 CD2 LEU A 313 14.757 -28.257 -23.179 1.00 54.46 C TER HETATM 3076 DUM DUM A 314 -4.000 -49.000 6.000 1.00 1.70 HETATM 3077 DUM DUM A 314 -4.000 -49.000 7.000 1.00 1.70 HETATM 3078 DUM DUM A 314 -4.000 -48.000 6.000 1.00 1.70 HETATM 3079 DUM DUM A 314 -4.000 -48.000 7.000 1.00 1.70 . . . . HETATM 36213 O HOH A 527 25.281 -35.299 13.147 1.00 37.03 O HETATM 36214 O HOH A 528 46.452 -27.955 -20.733 1.00 19.63 O HETATM 36215 O HOH A 530 32.439 -23.614 -26.720 1.00 27.12 O TER Best regards, Eva -- 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] Residue 'DUM' not found in residue topology database
Yes, meant that file (didn't know the name for GMX 4.5.x, i use mainly 4.0.7). You can use: editconf -f *.pdb -o *.gro to convert the *.pdb file into a *.gro file. Since the *.gro file has only informations about: atom names residue name coordinates while the *.pdb has also additional informations. Only problem could be the atom and residue names. But if i'm not mistaken, the *.pdb file needs the 'right' atom-names and residue names (the same atom names, which are also used in the *.rtp file) in order that 'pdb2gmx' can do its work (- producing the *.top). If you need the topology for the protein, you could delete all non-protein atoms and use 'pdb2gmx' to get it. Then later write the data for the dummy atoms into the *.top. But Justin hit an important point: You need parameters for the dummy - dummy dummy - protein / water interactions, because GROMACS has none, and would give you the next error. You can look into the literature if someone made parameters for this, probably you are lucky. Else you need to determine the parameters yourself, which is not so easy. Again chapter 5 would be your friend for how force fields work. And looking into coarse-graining papers could be a good idea, since they do something similar. Hope this helps and wish you luck for finding / getting parameters. Thomas Am 20.06.2012 18:13, schrieb gmx-users-requ...@gromacs.org: With the *nb.itp file you ment the ffnonbonded.itp file, right? But there is one thing: Of course you are right, I can modify the .top file and add my Dummy atom there. But because I got this error I mentioned, the pdb2gmx didn't produce a .gro file. Bests, Eva 'pdb2gmx' is mainly a tool to get the topology for (bigger) molecules from a *.pdb file. 'pdbgmx' is expecting an entry in the *.rtp file for the dummy atoms. You could add one in that file... (see further below) But if the dummys only interact via nonbonded interactions it would be more convient to write the topology by hand. In the *.top file, you would need after the protein a new entry with [ moleculetype ] and [ atoms ] But you should also have this nonbonded interaction in the *nb.itp and file of your force field and the atom definition in the *.atp file. If there are also bonded parameters the new *.rtp entry would be better. For this look at the definition of an aminoacid, or a single molecule and write something similar for your dummy atoms. Essential you would need the atoms definitions and charges, and then an entry for which atoms is bonded to which. But you should also have these interactions (nonbonded and bonded) in the *nb.itp and *bon.itp files of your force field and the atom definition in the *.atp file. The basic informations for both ways are in chapter 5 of the manual. Hope this gives some hints. Greetings Thomas Schlesier Am 20.06.2012 16:18, schriebgmx-users-requ...@gromacs.org: Hi everybody, I try to use GROMACS for my protein where I added a layer of dummy atoms simulating the membrane around it. But now when I want to call pdb2gmx I always get the error: Fatal error: Residue 'DUM' not found in residue topology database I understand the error, that there is no entry for the dummy residue in the .tpr file. But is there a possibility to add it there. Because the file don't look like that I can just write DUM in it. my pdb file looks like this: ATOM 2597 CB LEU A 313 12.816 -29.877 -23.547 1.00 55.31 C ATOM 2598 CG LEU A 313 14.261 -29.695 -23.059 1.00 54.06 C ATOM 2599 CD1 LEU A 313 15.169 -30.640 -23.866 1.00 53.26 C ATOM 2600 CD2 LEU A 313 14.757 -28.257 -23.179 1.00 54.46 C TER HETATM 3076 DUM DUM A 314 -4.000 -49.000 6.000 1.00 1.70 HETATM 3077 DUM DUM A 314 -4.000 -49.000 7.000 1.00 1.70 HETATM 3078 DUM DUM A 314 -4.000 -48.000 6.000 1.00 1.70 HETATM 3079 DUM DUM A 314 -4.000 -48.000 7.000 1.00 1.70 . . . . HETATM 36213 O HOH A 527 25.281 -35.299 13.147 1.00 37.03 O HETATM 36214 O HOH A 528 46.452 -27.955 -20.733 1.00 19.63 O HETATM 36215 O HOH A 530 32.439 -23.614 -26.720 1.00 27.12 O TER Best regards, Eva -- gmx-users mailing listgmx-us...@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 togmx-users-requ...@gromacs.org. Can't post? Readhttp://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
[gmx-users] vsite as interaction site for COM of benzene
Hi all, i have a more conceptional question, for using vsites as interaction-centers for coarse-grained particles: First the simple case: I want to simulate one benzene molecule (atomistic - aa) in coarse-grained (cg-) benzene (each benzene molecule as a single particle). For the cg-cg interaction of the benzene i have calculated a table potential. But for the aa-cg interaction I'm a little bit clueless: I want to use the COM of the aa-benzene as an interaction-site. But there is no vsite, which i could construct from 6 atoms. Side-note: The vsite would not interact with the with the atoms, only with cg-benzenes So there are two ideas: 1) Using one '3out'-vsite, which is construted from 3 non-neighbour atoms (if we had mesitylene instead of benzene, it would be the three C-atoms to which only hydrogens are bound). The the force from the aa-cg interaction would be distributed to these three atoms and i would hope that the bond-constraints would do 'the rest' (i.e. they would mimick that the whole aa-benzene molecule interacts with the cg-particles) 2) Using two '3out'-vsites. First viste the same as in (1), the second comes from the set of the other three C-atoms. These two vsites, should be at the same position (more or less, for a perfect energy-minimized benzene molecule it would be so; but i would assume the error is only minimal). Since the potential is pair-additive, i would use for this table potential only half of the potential, so that the potential of both particles would add up to the 'true' potential. The forces would be calculated from the 'half-potential'. Big questions is now: Are (1) and (2) equivalent? (2) Seems like the 'right way' to do. But (1) would be easier to implement :) But the big question is, if the bond-constraints would correct the approximation, that only 3 out of 6 atoms have a direct interaction with the cg-particles. Side note: The real system is a little more complex. There are fracments of molecules where i could not really construct two vsites which would have more or less the same position. In this case, method (1) would really help. Any ideas? Greetings Thomas -- 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] vsite as interaction site for COM of benzene
Dear Richard, thanks for the idea, i will try this. Relating to your side note: I perform pulling simulations of a calixarene-catenane dimer (both together 600 atoms) and i mesitylene as a solvent. To access longer time-scales I wanted to coarse-grain the solvent only. My supervisor said that the mesitylene rotates on such a short time-scale, that it would be ok, so coarse-grain it as one particle. But i also thought 3 particles would be an option. But since it is easier to coarse-grain the mesitylene as one particle, i started with that, at this stage, my first interest to get a somewhat reasonable simulation, which runs, since it is my first project involving coarse-graining and vsites. If the results are not good enough, i still can make the system more complex :) but it is easier to start with the simple case. I'm now at the to construct the AA-CG interaction. For this is started with an atomistic mesitylene in coarse-grained mesitylene (strictly speaking it's still an CG-CG interaction, but the single mesitylene has an atomistic representation, as later the calixarene). For the question here i used benzene instead of mesitylene, since it is smaller but the concept is the same. Dear Thomas, Whilst in principle a constraint should do what you are asking, you might be better off using multiple vsites to solve the problem so that there is an algebraic as apposed to numerical mapping of the forces. My immediate thought is to use two type 3 (or if necessary type 3fd) virtual sites to define two points near the centre of the benzene. You could then place your tabulated potential on a type 2 virtual site between these two. I believe that should algebraically link all the atomic sites together making the process more stable and it should be quite easy to transfer it to other systems. On a separate note why if you are coarse graining the benzenes to a single point particle are you then maintaining all the atoms? Could you go a stage further and map the atoms onto 3 non interacting point masses? That would maintain angular momentum etc. whilst reducing the computational overhead. Although it becomes very complex to implement if your benzene's are in the backbone of a polymer. Hope that helps, Richard On 08/06/2012 18:18, Thomas Schlesierschlesi at uni-mainz.de http://lists.gromacs.org/mailman/listinfo/gmx-users wrote: / Hi all, // // i have a more conceptional question, for using vsites as // interaction-centers for coarse-grained particles: // // First the simple case: // I want to simulate one benzene molecule (atomistic - aa) in // coarse-grained (cg-) benzene (each benzene molecule as a single // particle). For the cg-cg interaction of the benzene i have calculated a // table potential. // // But for the aa-cg interaction I'm a little bit clueless: // I want to use the COM of the aa-benzene as an interaction-site. But // there is no vsite, which i could construct from 6 atoms. // Side-note: The vsite would not interact with the with the atoms, only // with cg-benzenes // // So there are two ideas: // // 1) Using one '3out'-vsite, which is construted from 3 non-neighbour // atoms (if we had mesitylene instead of benzene, it would be the three // C-atoms to which only hydrogens are bound). The the force from the aa-cg // interaction would be distributed to these three atoms and i would hope // that the bond-constraints would do 'the rest' (i.e. they would mimick // that the whole aa-benzene molecule interacts with the cg-particles) // // 2) Using two '3out'-vsites. First viste the same as in (1), the second // comes from the set of the other three C-atoms. These two vsites, should // be at the same position (more or less, for a perfect energy-minimized // benzene molecule it would be so; but i would assume the error is only // minimal). Since the potential is pair-additive, i would use for this // table potential only half of the potential, so that the potential of // both particles would add up to the 'true' potential. The forces would be // calculated from the 'half-potential'. // // Big questions is now: // Are (1) and (2) equivalent? // (2) Seems like the 'right way' to do. // But (1) would be easier to implement :) But the big question is, if the // bond-constraints would correct the approximation, that only 3 out of 6 // atoms have a direct interaction with the cg-particles. // // Side note: // The real system is a little more complex. There are fracments of // molecules where i could not really construct two vsites which would have // more or less the same position. In this case, method (1) would really help. // // Any ideas? // Greetings // Thomas / -- 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
[gmx-users] Re: boundaries in PMF
One comment: If the channel is horizontal orientated, and the COM of MOL lies in the middle, you have two directions to go out of the channel: two the left side (quasi negative distance) and to the right side (quasi positive distance). What happens with 'pull_geometry=distance' is, that only the distance between MOL and Na matters. So you get the same (positive) distance if you go to the left or go to the right. The resulting PMF will be for only half of the channel (since the are no negative distances), but it will be evaluated from going into both directions (since the distances for left and right are the same). To distingiush if you go left or right, you would need to use 'pull_geometry=direction/position'. If you are 100% sure, that the system is symmetric and that the COM of MOL is in the middle of the channel, then your results are right. But if the COM of MOL is a little more to the right or left (relative to the middle of the channel) you would get problems. Since then the left and right direction are not equal. One thing you could check is. Make a analysis with g_wham on the windows which are on the left side of the center and one for the right side of the center. If everything is converged and the system is symmetric + COM of MOL is in the middle of the channel, both PMF should look the same. If they differ one or more of the above things are not correct. Am 31.05.2012 22:03, schrieb gmx-users-requ...@gromacs.org: OK, however, just one point about your last comment: I suspect this is why g_wham is finding a range of values only equal to half of your expected reaction coordinate; it is considering only the positive displacement along the reaction coordinate. It seems like all the channel is explored, not only one half. If g_wham was only considering the positive displacement I should see a profile for just one half of the channel, shouldn?t I? and I can see a profile typical for the entire channel. Rebeca. Date: Thu, 31 May 2012 15:42:06 -0400 From:jalem...@vt.edu To:gmx-users@gromacs.org Subject: Re: [gmx-users] Re: boundaries in PMF On 5/31/12 3:37 PM, Rebeca Garc?a Fandi?o wrote: Hi, the center of mass of my channel is at the middle of the ion channel, and it is a symmetric system, so I suppose these results would be OK. Anyway, I will check the options you propose. If you are sampling regions above and below the channel/membrane, then the distance geometry is not appropriate; you need either direction or (perhaps the most flexible option) position. There are a number of useful discussions on such topics in the list archive and in my umbrella sampling tutorial for you to consider. You can likely create a series of .tpr files from new .mdp files with correct options to run the analysis. I suspect this is why g_wham is finding a range of values only equal to half of your expected reaction coordinate; it is considering only the positive displacement along the reaction coordinate. -Justin Thanks a lot for all your comments!! Best wishes, Rebeca. Date: Thu, 31 May 2012 20:08:26 +0200 From:schl...@uni-mainz.de To:gmx-users@gromacs.org Subject: [gmx-users] Re: boundaries in PMF Where is the center of mass of reference group (MOL) located? It seems that the COM is near the middle of the ion channel. Since you use 'pull_geometry=distance', g_wham will look only for the distance between 'MOL' and 'Na' and that leads to problem. If the com of 'MOL' sits in the center of the channel (which is around 5nm long), one has 2.5nm in both directions. Since g_wham uses only the distance you get the PMF for half of the channel, but with the data of both parts. If the channel would be symmetric and the com of 'MOL' would lie exactly in the middle of the channel, this could be ok. But if even one of both assumptions is wrong, the results would be errorous. A better approach would be to use 'pull_geometry=direction', since the you define a vector along which the windows lie and do not have the problem that the distance is in both directions (along positive and negative vector) the same. Only problem could be that g_wham supports 'pull_geometry=direction' only from gromacs 4.5.x (don't know this, since instead of umbrella smapling i use thermodynamik integration, where one uses constraints (instead of restraints) and integrates the constraint-force; but the conceptual problem with 'distance/direction' is the same). Another approach (with 'pull_geometry=distance') would be to use a reference group which is just outside of the channel (like going some steps away from the channel, along the vector which goes through the channel). If there is only
[gmx-users] Re: boundaries in PMF
Where is the center of mass of reference group (MOL) located? It seems that the COM is near the middle of the ion channel. Since you use 'pull_geometry=distance', g_wham will look only for the distance between 'MOL' and 'Na' and that leads to problem. If the com of 'MOL' sits in the center of the channel (which is around 5nm long), one has 2.5nm in both directions. Since g_wham uses only the distance you get the PMF for half of the channel, but with the data of both parts. If the channel would be symmetric and the com of 'MOL' would lie exactly in the middle of the channel, this could be ok. But if even one of both assumptions is wrong, the results would be errorous. A better approach would be to use 'pull_geometry=direction', since the you define a vector along which the windows lie and do not have the problem that the distance is in both directions (along positive and negative vector) the same. Only problem could be that g_wham supports 'pull_geometry=direction' only from gromacs 4.5.x (don't know this, since instead of umbrella smapling i use thermodynamik integration, where one uses constraints (instead of restraints) and integrates the constraint-force; but the conceptual problem with 'distance/direction' is the same). Another approach (with 'pull_geometry=distance') would be to use a reference group which is just outside of the channel (like going some steps away from the channel, along the vector which goes through the channel). If there is only water, it would be bad, because then the reference group would be move away. But then on could use the entry and exit of the channel as a reference point for two simulations. In the case that the entry is the reference group, the PMF would be ill defined near the entry, but from the simulation with the exit as reference you would get the right PMF for the entry region and vice versa. Greetings Thomas Am 31.05.2012 19:39, schrieb gmx-users-requ...@gromacs.org: Thanks a lot for your quick answer. The mdp file I have used is copied below. What is strange is that when I look at the *gro files for the different windows (100 windows in total), i. e: window 1: 8704Na Na56458 5.134 5.085 5.824 window 50: 8704Na Na56458 5.053 5.081 3.459 window 100: 8704Na Na56458 4.990 5.042 0.951 you can see that the z-coordinate goes from 0.951 to 5.824 nm I should have a total distance in the x-axis of about 5 nm, and however, the boundaries calculated by g_wham are Determined boundaries to 0.35 to 2.603290 Can you see anything in the mdp file which is causing this behaviour...? Thanks again for your help! MDP FILE USED: title = Umbrella pulling simulation define = -DPOSRES_B define = ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 50 ; 1 ns nstcomm = 10 ; Output parameters nstxout = 5000 ; every 10 ps nstvout = 5000 nstfout = 5000 nstxtcout = 5000 ; every 10 ps nstenergy = 5000 ; Bond parameters constraint_algorithm = lincs constraints = all-bonds continuation = yes ; Single-range cutoff scheme nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb = 1.4 rvdw = 1.4 ; PME electrostatics parameters coulombtype = PME fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 optimize_fft = yes ;Berendsen temperature coupling is on Tcoupl = V-rescale tau_t = 0.1 0.1 0.1 tc-grps = MOL dop WAT_Cl_Na ref_t = 300 300 300 ; Pressure coupling Pcoupl = Parrinello-Rahman Pcoupltype = Semiisotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) tau-p = 1.0 1.0 compressibility = 4.6E-5 4.6E-5 ref-p = 1.0 1.0 ; Generate velocities is off gen_vel = no ; Periodic boundary conditions are on in all directions pbc = xyz ; Long-range dispersion correction DispCorr = EnerPres ; Pull code pull = umbrella pull_geometry = distance pull_dim = N N Y pull_start = yes pull_ngroups = 1 pull_group0 = MOL pull_group1 = Na pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1000 ; kJ mol^-1 nm^-2 pull_nstxout = 1000 ; every 2 ps pull_nstfout = 1000 ; every 2 ps Date: Thu, 31 May 2012 13:25:26 -0400 From:jalem...@vt.edu To:gmx-users@gromacs.org Subject: Re: [gmx-users] boundaries in PMF On 5/31/12 1:20 PM, Rebeca Garc?a Fandi?o wrote: Hi, I am trying to calculate a PMF for an ion along a channel. Everything went OK, but when I used g_wham I got a profile with strange dimensions for the x-axis. What are the boundaries g_wham is using for calculating the units of x-axis? The values are based on the restraint distances along the reaction coordinate. I have used: g_wham -it tpr-files.dat -if pullf-files.dat -o file_output.xvg -hist file_histo_output.xvg -unit kCal -cycl yes (version 4.0.7) and the units I got in the x-axis are determined by the boundaries: Determined boundaries to 0.35 to 2.603290 Which are these units? nm? Yes. The z coordinate for my ion explores at least 5 nm!! I am a bit confused about
[gmx-users] Re: Justin umbrella sampling tutorial......
I think all the answers to your question are in the tutorial. Probably read first the lysozyme tutorial and then the umbrella tutorial again. But here is a more general answer: Normally you have two types of simulations: preperation (which is also equilibration) production and the need of restraints depends on what you want to achive. So if you want to calculate some equilibrium property of a protein in water you do first a preperation simulation to equilibrate the system (NVE, NVT or NPT - depending for what ensemble you want the property). Normally during this you put position restraints to the protein backbone, so that the protein structure does not gets disturbed during the part where you equilibrate the water. but if your protein is fairly stable / rigid, you don't need these restraints. Then you do the production simulation. Normally without restraints, since you want to know the property for a protein in equilibrium. But you could also do this with the same position restraints, in the extrem case, when you freeze the protein, you would get the property for a special protein geometry (the one to which ou restrain the coordinates). - So it really depends on what you want to do!!! For the umbrella tutorial: As far as i remember, Justin used the restrains for one of the chains (B), so to conserve the structure of said chain, because when one pull chain A away this would have an influence to the structure of chain B (they interact with each other). BUT: one could also dismiss these restrains during the umbrella simulations, since it is expected the structure of chain B depends on the distance to chain A (since they interact). You could also restrain only that part of chain B, which does not interact directly with chain A, to conserve a part of the structure of B (so you would be in the middle between the two extreme cases above). In all three cases you get a different answer (like above for the equilibrium property). One could argue which is the better option, but i think all three have their right to exist, since the answer to this also depends on the initial question. Best would be to do all three cases, but mostly one does not have so much time. So one must stick to one. Which leads to: First think about your question. Then think about how to answer it (think often there will be different possibilities / methods). Then decide what to do... Since the answer you get (may) depend on the method you use, you should justify your decission. Hope this helps greetings thomas Thank you Justin for these correct explanation Its really clear my lot of queries.. For the tutorial, NPT is conducted with restraints on all protein heavy atoms. The production runs are conducted by restraining only one chain for practical reasons. These is my question ; If we are doing NPT with restraints on all protein atom and production run by conducted by restraining only one chain for protein... means NPT and productin run mdp files should be different , Where these information in mdp files??? It my request to you please tell me why these mdp files are almost same in parameter .. (Where you mentioned in mdp for npt to do restraints on all proteins heavy atoms and for production md restraining only one chain ..) I don't understand what more I can say beyond what I already have. -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] Re: Wierd results from Umbrella sampling
The picture of pullx.xvg is not really helpful. The reason is that in pullx.xvg the second column (after the time) is the position of the reference group (which naturally moves; but it's not the thing in which wie are interested) and the third column is the distance from the reference group to the pulled group (this is the one which interestes you). Another thing you could try is to pull in all 3 dimensions, if you pull only in the z-direction your pulled group can move freely in the x-y plane and wouldn't be affected by the umbrella potential. But i think for the histogram only the direction in which one pulls are accounted for, so the wide histograms are somewhat strange. For the problems with the constraints i have absolute no idea, never had this error message. Even don't know if the problems you get from there follow you into the pull-mode simulation. Does the constraint simulation run without the pull-code, or do you get there the same problem? greetings thomas I increased the potential to 2350 instead of 1000 in my md_umbrella.mdp (the other parameters kept unchanged), then the movement of my protein was less than the previous one, which we can see from pullx.xvg (attached). But it still not reach to the designed region. So I am going to increase the potential higher again. It still does not work if I switched on the LINCS. The error is that the number of constraint is more than the number of freedom. I am running another trail with high potential (let's say 3000, ok?) and dt=2 ft, which value is definitely not recommended by Martini developers. With Regards, Jiangfeng. Message: 3 Date: Tue, 22 May 2012 18:53:54 +0200 From: Thomas Schlesierschl...@uni-mainz.de Subject: [gmx-users] Re: Wierd results from Umbrella sampling, (Justin A. Lemkul) To:gmx-users@gromacs.org Message-ID:4fbbc4a2.9020...@uni-mainz.de Content-Type: text/plain; charset=ISO-8859-1; format=flowed I never worked with the MARTINI (or other coarse-grained) force field, but this in the umbrella.mdp title = Umbrella pulling simulation integrator = md dt = 0.019 looks suspicious. The dt is about an order of magnitude greater than one uses in normal (bond-)constrainted md-simulations. I know that in coarse-grained simulations the potentials are smoother than in atomistic force fields and one therefore could use a higher timestep, but i don't know if you can go so high. Anyway you should look for the umbrella potential. If the force constant and the timestep are both too high you could get problems: Assume a particle in a harmonic well. If the force constant is high and the timestep too high, you wont sample the harmonic well, but 'jump' each time from one side to the other (high force + great timestep - long movement). So i would test it first with a timestep of 0.002 (same as you used for pulling sim). -- Message: 4 Date: Tue, 22 May 2012 19:15:41 +0200 From: Justin A. Lemkuljalem...@vt.edu Subject: Re: [gmx-users] Re: Wierd results from Umbrella sampling, (Justin A. Lemkul) To: Discussion list for GROMACS usersgmx-users@gromacs.org Message-ID:4fbbc9bd.8070...@vt.edu Content-Type: text/plain; charset=ISO-8859-1; format=flowed On 5/22/12 6:53 PM, Thomas Schlesier wrote: I never worked with the MARTINI (or other coarse-grained) force field, but this in the umbrella.mdp title = Umbrella pulling simulation integrator = md dt = 0.019 looks suspicious. The dt is about an order of magnitude greater than one uses in normal (bond-)constrainted md-simulations. I know that in coarse-grained simulations the potentials are smoother than in atomistic force fields and one therefore could use a higher timestep, but i don't know if you can go so high. Anyway you should look for the umbrella potential. If the force constant and the timestep are both too high you could get problems: Assume a particle in a harmonic well. If the force constant is high and the timestep too high, you wont sample the harmonic well, but 'jump' each time from one side to the other (high force + great timestep - long movement). So i would test it first with a timestep of 0.002 (same as you used for pulling sim). Testing this would be good, though the MARTINI developers routinely report timesteps of 20-40 fs as normal use. I have never obtained stable trajectories above 20 fs, but I also do not do much coarse graining. It could very well be that the pull code is not stable with such an integration step. -Justin Am 22.05.2012 18:37, schrieb gmx-users-requ...@gromacs.org: Thank you for your reply, Thomas. Under your explanation, I am well understood of the terms: pull_k and pull_rate. Here I upload both md_umbrella.mdp and md_pulling.mdp. I have to mention that it is coarse gained system with Martini force field. At the same time, i am going to run a simulation without pull code but with LINCS constraint, and also run another system with a huge pull_k but without
[gmx-users] Re: Wierd results from Umbrella, sampling (Justin A. Lemkul)
Think it would be best to show the .mdp file, else we can only guess what the parameters are. From the histogram it looks like that the force constant of the restraining potential is too low, since the histograms are really wide, but pull_k1=1000 is a 'normal' value. On thing which confueses me is you said that fluctuations from g_dist are about 0.4nm, but in the histogram i looks like the distance fluctuates about 1nm or even more. for this be sure, that you compare the right distances - if you pull only in z-direction, the only look into the z-direction from the output from g_dist. Since the x, and y-directions would be unaffected by the umbrella potential. for the error message with the constraints: what happens if you run the system with constraints but without the pull-code? for pull_k10 and pull_rate1=0: if you're pulling with an umbrella potential pull_k1 defines the force constant of the potential (hook'sches law). Let's assume you put a spring to your molecule with pull_rate1=0, it means that the other end of the spring doesn't move, and the spring restrains the position of the molecule (molecule cn't diffuse away). in umbrella sampling you want to restrian of molecule (or part of it) relative to another one / part - so pull_rate1=0 with pull_rate10, it's like you pull the spring away for the molecule, spring gets strechted - wants to relax - molecule follows spring (like that what happens in afm pulling) Am 22.05.2012 15:40, schrieb gmx-users-requ...@gromacs.org: Hi Justin, As for the maximum of 0.4 nm fluctuation of my pulled protein, I used g_dist to calculate the distance between my pulled protein and stable part in a window, where the distance is treated as the fluctuation of the protein along z-axis. Well, from pullx.xvg, the position moved a lot (3nm for instance.) As for the windows simulation, I didn't apply constraint but only the internal constraint in the itp file. I still don't understand why it have to do constraint? why not give a fully freedom to run the simulation? In the md_umbrella.mdp, I set pull_k1=1000; pull_rate1= 0.0, but apparently, I am confused with pull_k10 combined with pull_rate1=0. In my mdp, i set none to LINCS, because if I use all-bonds, an error of 1099 constraints but degrees of freedom is only1074 occurs. Actually, there is no any window with a designed distance. Here I attach the histo.xvg, profile.xvg. Thank you with regards, Jiangfeng. On 5/22/12 9:36 AM, Du Jiangfeng (BIOCH) wrote: Dear Justin, Based on your questions to my simulation, I posted here yesterday hopefully it was the correct way to reply in this forum. You've still replied to the entire digest message (which I've cut out); please make sure to keep replies free of superfluous posts in the future. The archive is already pretty hopeless, but let's not make it worse:) In this morning I got a list of new windows of umbrella sampling, the overlap is sufficient enough, but I saw another problem: In the histogram figure, the base of peak covers the distance of 2 nm instead of 0.2 nm., that's horrible! However, when I checked back to the simulation results of each window, the fluctuation of my pulled protein is only 0.4nm in maximum. So the base of peak shouldn't cover such long distance, right? If the peaks aren't corresponding to the desired restraint distances, then there are several potential problems: 1. Your restraints aren't set up the way you think they are (check grompp output and .tpr file contents to be sure) 2. Your restraints are ineffectual (in which case you may need to revisit the force constant) I can't determine from your description what's going on. What do you mean by a maximum of 0.4 nm fluctuation? In what quantity? What do the contents of pullx.xvg show you for the problematic window, and for that matter, the others? Are any of them producing the desired restraint distances? Again I will ask you to share an image of the histogram and PMF profile; these would be very helpful to see. -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] Re: Wierd results from Umbrella sampling, (Justin A. Lemkul)
I never worked with the MARTINI (or other coarse-grained) force field, but this in the umbrella.mdp title = Umbrella pulling simulation integrator = md dt = 0.019 looks suspicious. The dt is about an order of magnitude greater than one uses in normal (bond-)constrainted md-simulations. I know that in coarse-grained simulations the potentials are smoother than in atomistic force fields and one therefore could use a higher timestep, but i don't know if you can go so high. Anyway you should look for the umbrella potential. If the force constant and the timestep are both too high you could get problems: Assume a particle in a harmonic well. If the force constant is high and the timestep too high, you wont sample the harmonic well, but 'jump' each time from one side to the other (high force + great timestep - long movement). So i would test it first with a timestep of 0.002 (same as you used for pulling sim). Am 22.05.2012 18:37, schrieb gmx-users-requ...@gromacs.org: Thank you for your reply, Thomas. Under your explanation, I am well understood of the terms: pull_k and pull_rate. Here I upload both md_umbrella.mdp and md_pulling.mdp. I have to mention that it is coarse gained system with Martini force field. At the same time, i am going to run a simulation without pull code but with LINCS constraint, and also run another system with a huge pull_k but without LINCS. Hope I could get some interesting information. With Regards, Jiangfeng. Message: 3 Date: Tue, 22 May 2012 16:36:47 +0200 From: Thomas Schlesierschl...@uni-mainz.de Subject: [gmx-users] Re: Wierd results from Umbrella, sampling (Justin A. Lemkul) To:gmx-users@gromacs.org Message-ID:4fbba47f.5010...@uni-mainz.de Content-Type: text/plain; charset=ISO-8859-1; format=flowed Think it would be best to show the .mdp file, else we can only guess what the parameters are. From the histogram it looks like that the force constant of the restraining potential is too low, since the histograms are really wide, but pull_k1=1000 is a 'normal' value. On thing which confueses me is you said that fluctuations from g_dist are about 0.4nm, but in the histogram i looks like the distance fluctuates about 1nm or even more. for this be sure, that you compare the right distances - if you pull only in z-direction, the only look into the z-direction from the output from g_dist. Since the x, and y-directions would be unaffected by the umbrella potential. for the error message with the constraints: what happens if you run the system with constraints but without the pull-code? for pull_k10 and pull_rate1=0: if you're pulling with an umbrella potential pull_k1 defines the force constant of the potential (hook'sches law). Let's assume you put a spring to your molecule with pull_rate1=0, it means that the other end of the spring doesn't move, and the spring restrains the position of the molecule (molecule cn't diffuse away). in umbrella sampling you want to restrian of molecule (or part of it) relative to another one / part - so pull_rate1=0 with pull_rate10, it's like you pull the spring away for the molecule, spring gets strechted - wants to relax - molecule follows spring (like that what happens in afm pulling) Am 22.05.2012 15:40, schriebgmx-users-requ...@gromacs.org: Hi Justin, As for the maximum of 0.4 nm fluctuation of my pulled protein, I used g_dist to calculate the distance between my pulled protein and stable part in a window, where the distance is treated as the fluctuation of the protein along z-axis. Well, from pullx.xvg, the position moved a lot (3nm for instance.) As for the windows simulation, I didn't apply constraint but only the internal constraint in the itp file. I still don't understand why it have to do constraint? why not give a fully freedom to run the simulation? In the md_umbrella.mdp, I set pull_k1=1000; pull_rate1= 0.0, but apparently, I am confused with pull_k10 combined with pull_rate1=0. In my mdp, i set none to LINCS, because if I use all-bonds, an error of 1099 constraints but degrees of freedom is only1074 occurs. Actually, there is no any window with a designed distance. Here I attach the histo.xvg, profile.xvg. Thank you with regards, Jiangfeng. On 5/22/12 9:36 AM, Du Jiangfeng (BIOCH) wrote: Dear Justin, Based on your questions to my simulation, I posted here yesterday hopefully it was the correct way to reply in this forum. You've still replied to the entire digest message (which I've cut out); please make sure to keep replies free of superfluous posts in the future. The archive is already pretty hopeless, but let's not make it worse:) In this morning I got a list of new windows of umbrella sampling, the overlap is sufficient enough, but I saw another problem: In the histogram figure, the base of peak covers the distance of 2 nm instead of 0.2 nm., that's horrible! However, when I checked
[gmx-users] Re: Proper 1-octanol box preparation
Hi, relating to the picture see: http://www.gromacs.org/Documentation/Terminology/Periodic_Boundary_Conditions It's just a matter of periodic boundary conditions, you don't have two, but just one cluster of octanol. What i don't understand is, that your box doesn't collapse (becoming smaller) if you perfom the NPT simulation. Only thing is that i would search the value of compressibility of octanol, but this couldn't be the origin of the problem. Greetings Thomas Am 03.05.2012 12:23, schrieb gmx-users-requ...@gromacs.org: Hi GMX users! I prepared a box of 125 1-octanol molecules using genconf -f octanol_single_molecule.gro -nbox 5 5 5 and I tried to equilibrate this system in NPT ensemble to get proper density and nice cube box, similar to octanol configuration that one can download from http://www.gromacs.org/Downloads/User_contributions/Molecule_topologies. But I always got something completely different. Some cluster of the molecules at the bottom of the simulation box, far away from dense cube configuration. You can see picture here: https://picasaweb.google.com/okroutil/Science#5738216101659210130. I tried different thermo- and barostats, different coupling times, compressibility, first NVT then NPT equil., but with no success. I also tried more chaotic initial configuration generated with genbox -ci ...gro - nmol 125 -box 5 5 5, no success. So my question is: is there some simulation protocol or some setting (higher pressure at the beginning of equilibration?), any trick suitable for this type of solvent? Till now I worked only with water and surfaces, so this is new area for me... Thank you very much for any answer and have a nice day! Ondrej Kroutil (Faculty of Health and Social Studies, South Bohemian University, Czech Republic) P.S.: I used gaff ff and topologies downloaded from http://virtualchemistry.org/molecules/111-87-5/index.php and this NPT eq. input: integrator = md dt = 0.002 nsteps = 50 comm_mode= linear nstcomm = 1000 nstxout = 0 nstxtcout= 100 nstvout = 0 nstfout = 0 nstlog = 500 nstlist = 10 ns_type = grid rlist= 1.4 coulombtype = PME rcoulomb = 1.4 rvdw = 1.4 constraints = none constraint_algorithm = lincs ;shake_tol = 0.1 lincs_iter = 1 fourierspacing = 0.1 pme_order = 4 ewald_rtol = 1e-5 ewald_geometry = 3d optimize_fft= yes ; Nose-Hoover temperature coupling Tcoupl = berendsen tau_t = 1 tc_grps= system ref_t = 298. ; No Pressure Pcoupl = berendsen pcoupltype = isotropic tau_p = 0.5 compressibility = 4.6e-5 ref_p = 1.0 ; OTHER periodic_molecules = no pbc = xyz gen_vel = yes gen_temp= 298.15 gen_seed= -1-- Ond??ej Kroutil -- 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] Extract pullF from trajectory
Am 27.04.2012 15:31, schrieb gmx-users-requ...@gromacs.org: Hi Gmx Users, I run umbrella sampling simulation and for one window I lost my files: pullf.xvg and pullx.xvg. Is there any way to extract it based on the trajectory from this window? Steven Should be possible, but you will be limited to the output-rate of your trajectory: You should be able to extract the position of the pulled (PULL) and reference (REF) group from the trajectory. From this you can recalculate the pullx.xvg (think it's position of REF, then vector connecting both groups). From the position of REF and the input for the pulling, you can calculate the position of the origin of the umbrella potential. Then just calculate the distance of the PULL to this origin to get the 'extension' - multiply by the force constant and you get the force. If you had 'pull_geometry=distance' it's trivial: (1) 'pull_init' is distance REF to origin (2) 'g_dist' to get distance from REF to PULL (1)-(2) gives distance of PULL to origin * force constant - force For the other options it's also straight forward: in the end you want to know the distance of the umbrella origin - PULL * from input you know the distance REF - origin * 'g_dist' gives you the distance REF - PULL Thing to note: * Depending on the using of 'pull_start' you need to modify 'pull_init' If you have 'pull_start=yes', just use 'g_dist' to get the initial distance of REF to PULL, add this to 'pull_init' * If you pull only in one dimension, just only use this dimension and not just the whole vector Greetings Thomas -- 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] -com switch in g_rdf
Using the option '-com' results in a RDF were the reference group is the center of mass (COM) of the first group. - In a system with pure benzene and using the whole system as a group, this reference point would be somewhere in the middle of the box. And choosing this group twice, would result in a RDF of the distances of every benzene atom to the reference point (somewhere in the middle of the box). What you are looking is the option: -rdf: -rdf enum atomRDF type: atom, mol_com, mol_cog, res_com or res_cog whereas 'atom' is the standard option (calculating RDF between every atom of both groups) For calculating the benzene COM-COM RDF you would need the option '-rdf mol_com' and choosing the index-group which consists of all benzene atoms twice. Hope this helps Greetings Thomas --- Hi, If you have time, I have a somewhat embarrassing question to ask. It is embarrassing because I feel like I should be able to understand or search for the answer, but I have not had success. So I apologize if this is a silly question (and it probably is). On the mailing list, I have read this post from several years ago: http://lists.gromacs.org/pipermail/gmx-users/2006-May/021743.html The poster asks if one can use the -com option in g_rdf to calculate the (benzene center-of-mass)-(benzene center-of-mass) radial distribution function. Dr. van der Spoel responds and says that the -com option only computes COM of the central group. What you want is not implemented. My question is, what is the central group? g_rdf asks the user to select two entries for computing the RDF. Typically I have been doing atom-atom RDFs. This is straightforward; if I have a system of water with atoms OW and HW, and I want to compute the OW-HW RDF, I just need to make two index file entries: one with all of the OW atoms, and one with all of the HW atoms. But what happens if I have a system of benzene and I want to compute the (benzene center-of-mass)-(benzene center-of-mass). This is not possible in the current implementation, I guess, because the index file does not know about molecules (only about atoms). My question is, what if I selected all of the atoms in my pure benzene system, and used this selection for both selections as prompted by g_rdf. What would this compute if I used the -com switch? Since the index file does not know about molecules, it wouldn't compute the center mass of each benzene molecule. But my question is, of what would it compute the center of mass? Thank you! Andrew DeYoung Carnegie Mellon University -- 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
AFAIK this factor isn't included in g_wham. Think i tested this some years ago with an older GROMACS version. If one looks into the 'gmx_wham.c' and searches for 'ln' one finds it only in the comments for the gaussian random numbers. Greetings Thomas 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] Re: g_wham problem with negative COM differences
Am 18.04.2012 12:00, schrieb gmx-users-requ...@gromacs.org: Send gmx-users mailing list submissions to gmx-users@gromacs.org To subscribe or unsubscribe via the World Wide Web, visit http://lists.gromacs.org/mailman/listinfo/gmx-users or, via email, send a message with subject or body 'help' to gmx-users-requ...@gromacs.org You can reach the person managing the list at gmx-users-ow...@gromacs.org When replying, please edit your Subject line so it is more specific than Re: Contents of gmx-users digest... Today's Topics: 1. Re: g_wham problem with negative COM differences (Anni Kauko) 2. Error- Atom not found in residue seq nr while adding atom (aiswarya pawar) -- Message: 1 Date: Wed, 18 Apr 2012 12:19:08 +0300 From: Anni Kaukoaka...@sbc.su.se Subject: Re: [gmx-users] g_wham problem with negative COM differences To: gmx-users@gromacs.org Message-ID: calxandpfxzgyhxdvbsqmsrfgy7-man2f_-vdkbpkrgtzq-t...@mail.gmail.com Content-Type: text/plain; charset=iso-8859-1 Anni Kauko wrote: Date: Wed, 11 Apr 2012 08:38:05 -0400 From: Justin A. Lemkuljalem...@vt.edumailto:jalem...@vt.edu Subject: Re: [gmx-users] g_wham problem with negative COM differences To: Discussion list for GROMACS usersgmx-users@gromacs.org mailto:gmx-users@gromacs.org Message-ID:4f857b2d.3050...@vt.edu mailto:4f857b2d.3050...@vt.edu Content-Type: text/plain; charset=ISO-8859-1; format=flowed Anni Kauko wrote: Hi! I try to perform pmf calculations for case where a peptide shifts through the membrane. My COM differences should vary from 2.3 to -2.5. My problem is that g_wham plots negative COM difference as they would be positive. In pullx-files the COM differences are treated correctly (look below). My peptide is not symmetric, so profile curves are not symmetric, so loosing the sign for COM difference screws my profile curve completely. I did not manage to find any pre-existing answers to this problem from internet. First datalines from pullx files: (sorry for strange file names...) pull_umbr_0.xvg: 0. 6.26031 2.27369 pullz_umbr_23.xvg: 0. 6.09702 0.0233141 pullz_umbr_50.xvg: 0. 6.02097 -2.50088 g_wham command: g_wham -b 5000 -it tpr_files.dat -ix pullz_files.dat -o profile_test.xvg -hist histo_test.xvg -unit kCal My pull code: pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = POPC_POPS ; reference group is bilayer pull_group1 = C-alpha__r_92-94 ; group that is actually pulled pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1000 ; kJ mol-1 nm-2 Your problem stems from the use of distance geometry. This method assumes the sign along the reaction coordinate does not change, i.e. always positive or always negative. If the sign changes, this simple method fails. You should be using something like position to allow for a vector to be specified. Perhaps you can reconstruct the PMF by separately analyzing the positive restraint distances and negative restraint distances (note here that distance really refers to a vector quantity, and thus it can have a sign), or otherwise create new .tpr files using position geometry, though I don't know if g_wham will accept them or not. -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.eduhttp://vt.edu | (540) 231-9080 tel:%28540%29%20231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin Thank's! I managed to solve my g_wham problem by doing two things: 1. New tpr-files with proper pull code for g_wham. 2. I also needed to modify signs of pullf values: If value for pullx distance was negative, I reversed the sign of corresponding pullf value. I did that by my own script. The new pull code: ; Pull code pull= umbrella pull_geometry = direction pull_vec1 = 0 0 1 pull_start = yes pull_ngroups= 1 pull_group0 = POPC_POPS ; reference group is bilayer pull_group1 = C-alpha__r_92-94 ; group that is actually pulled pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1000 ; kJ mol-1 nm-2 I am little bit confuced, why I needed to
[gmx-users] g_wham problem with negative COM differences
New try, think last time message didn't reached the list :( Original-Nachricht Betreff: Re: g_wham problem with negative COM differences Datum: Fri, 13 Apr 2012 14:55:15 +0200 Von: Thomas Schlesier schl...@uni-mainz.de An: gmx-users@gromacs.org Anni Kauko wrote: Date: Wed, 11 Apr 2012 08:38:05 -0400 From: Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu Subject: Re: [gmx-users] g_wham problem with negative COM differences To: Discussion list for GROMACS users gmx-users@gromacs.org mailto:gmx-users@gromacs.org Message-ID: 4f857b2d.3050...@vt.edu mailto:4f857b2d.3050...@vt.edu Content-Type: text/plain; charset=ISO-8859-1; format=flowed Anni Kauko wrote: Hi! I try to perform pmf calculations for case where a peptide shifts through the membrane. My COM differences should vary from 2.3 to -2.5. My problem is that g_wham plots negative COM difference as they would be positive. In pullx-files the COM differences are treated correctly (look below). My peptide is not symmetric, so profile curves are not symmetric, so loosing the sign for COM difference screws my profile curve completely. I did not manage to find any pre-existing answers to this problem from internet. First datalines from pullx files: (sorry for strange file names...) pull_umbr_0.xvg: 0. 6.26031 2.27369 pullz_umbr_23.xvg: 0. 6.09702 0.0233141 pullz_umbr_50.xvg: 0. 6.02097 -2.50088 g_wham command: g_wham -b 5000 -it tpr_files.dat -ix pullz_files.dat -o profile_test.xvg -hist histo_test.xvg -unit kCal My pull code: pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = POPC_POPS ; reference group is bilayer pull_group1 = C-alpha__r_92-94 ; group that is actually pulled pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1000 ; kJ mol-1 nm-2 Your problem stems from the use of distance geometry. This method assumes the sign along the reaction coordinate does not change, i.e. always positive or always negative. If the sign changes, this simple method fails. You should be using something like position to allow for a vector to be specified. Perhaps you can reconstruct the PMF by separately analyzing the positive restraint distances and negative restraint distances (note here that distance really refers to a vector quantity, and thus it can have a sign), or otherwise create new .tpr files using position geometry, though I don't know if g_wham will accept them or not. -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu http://vt.edu | (540) 231-9080 tel:%28540%29%20231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin Thank's! I managed to solve my g_wham problem by doing two things: 1. New tpr-files with proper pull code for g_wham. 2. I also needed to modify signs of pullf values: If value for pullx distance was negative, I reversed the sign of corresponding pullf value. I did that by my own script. The new pull code: ; Pull code pull= umbrella pull_geometry = direction pull_vec1 = 0 0 1 pull_start = yes pull_ngroups= 1 pull_group0 = POPC_POPS ; reference group is bilayer pull_group1 = C-alpha__r_92-94 ; group that is actually pulled pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1000 ; kJ mol-1 nm-2 I am little bit confuced, why I needed to tweak signes of pullf values. But like that I got the curve that resembles two half curve made for positive and negative pullx distances separately. That curve also makes sense from biochemical point of view. Such changes do not seem appropriate to me. If you change the sign of the pulling force, you change the implication of what that value means. What happens if you run your simulations with the new (more appropriate) .mdp file? Do the forces have the same magnitude, but opposite sign? I don't think that the problem is so easy fixed. I had with another person about a month ago a lengthy discussion on the list. The problem is the following: If you use 'distance' as pull_geometry the position of the minimum of the umbrella potential is determined by the vector connecting the ref- and pulled group, but there is no information about
[gmx-users] Re: g_wham problem with negative COM differences
Anni Kauko wrote: Date: Wed, 11 Apr 2012 08:38:05 -0400 From: Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu Subject: Re: [gmx-users] g_wham problem with negative COM differences To: Discussion list for GROMACS users gmx-users@gromacs.org mailto:gmx-users@gromacs.org Message-ID: 4f857b2d.3050...@vt.edu mailto:4f857b2d.3050...@vt.edu Content-Type: text/plain; charset=ISO-8859-1; format=flowed Anni Kauko wrote: Hi! I try to perform pmf calculations for case where a peptide shifts through the membrane. My COM differences should vary from 2.3 to -2.5. My problem is that g_wham plots negative COM difference as they would be positive. In pullx-files the COM differences are treated correctly (look below). My peptide is not symmetric, so profile curves are not symmetric, so loosing the sign for COM difference screws my profile curve completely. I did not manage to find any pre-existing answers to this problem from internet. First datalines from pullx files: (sorry for strange file names...) pull_umbr_0.xvg: 0. 6.26031 2.27369 pullz_umbr_23.xvg: 0. 6.09702 0.0233141 pullz_umbr_50.xvg: 0. 6.02097 -2.50088 g_wham command: g_wham -b 5000 -it tpr_files.dat -ix pullz_files.dat -o profile_test.xvg -hist histo_test.xvg -unit kCal My pull code: pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = POPC_POPS ; reference group is bilayer pull_group1 = C-alpha__r_92-94 ; group that is actually pulled pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1000 ; kJ mol-1 nm-2 Your problem stems from the use of distance geometry. This method assumes the sign along the reaction coordinate does not change, i.e. always positive or always negative. If the sign changes, this simple method fails. You should be using something like position to allow for a vector to be specified. Perhaps you can reconstruct the PMF by separately analyzing the positive restraint distances and negative restraint distances (note here that distance really refers to a vector quantity, and thus it can have a sign), or otherwise create new .tpr files using position geometry, though I don't know if g_wham will accept them or not. -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu http://vt.edu | (540) 231-9080 tel:%28540%29%20231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin Thank's! I managed to solve my g_wham problem by doing two things: 1. New tpr-files with proper pull code for g_wham. 2. I also needed to modify signs of pullf values: If value for pullx distance was negative, I reversed the sign of corresponding pullf value. I did that by my own script. The new pull code: ; Pull code pull= umbrella pull_geometry = direction pull_vec1 = 0 0 1 pull_start = yes pull_ngroups= 1 pull_group0 = POPC_POPS ; reference group is bilayer pull_group1 = C-alpha__r_92-94 ; group that is actually pulled pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1000 ; kJ mol-1 nm-2 I am little bit confuced, why I needed to tweak signes of pullf values. But like that I got the curve that resembles two half curve made for positive and negative pullx distances separately. That curve also makes sense from biochemical point of view. Such changes do not seem appropriate to me. If you change the sign of the pulling force, you change the implication of what that value means. What happens if you run your simulations with the new (more appropriate) .mdp file? Do the forces have the same magnitude, but opposite sign? I don't think that the problem is so easy fixed. I had with another person about a month ago a lengthy discussion on the list. The problem is the following: If you use 'distance' as pull_geometry the position of the minimum of the umbrella potential is determined by the vector connecting the ref- and pulled group, but there is no information about the direction. Assume this (1D): reference particle at origin (0) minimum of umbrella potential (1) - via pull_init1 = 1 pulled particle
[gmx-users] pull-code
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
[gmx-users] Umbrella Pulling
As far as i remember for PBC the distance between the reference and the pulled group are relevant, and not the distance between the reference group and the virtual spring (place where the pulling potential is zero). Looking into the code of GROMACS-4.0.7 backs this. If the distance between the reference group and spring would be relevant, one should see first an increasing force, then when pbc matters, the force should decrease, because we pull now in the other direction. Looking into the code of GROMACS-4.5.4, it seems there is an error message for pbc violations: gmx_fatal(FARGS,Distance of pull group %d (%f nm) is larger than 0.49 times the box size (%f),g,sqrt(dr2),max_dist2); pull.c; line 329 One thing is that the force constant for the pulling is quite small, so probably you need more force the move the ligand? (But i can't really comment on this, because i can't find your force vs time plot, so i don't know how large the force actually is). Greetings Thomas Steven Neumann wrote: Hello Justin, As you recommended I run longer pulling of my ligand. I pull the ligand which is on the top of my protein so the top of the protein is pulled so that the whole protein rotated app. 30 degrees - the low part of the protein came out of the box due to the rotation. After app 600 ps the ligand does not move any more and the force applied increase linearly. As I saw the trajectory, ligand does not collide with a periodic image but it does not move. Please, see attacheg plot force vs time. Can you explain why the force increase linearly and ligand is not pulled any more? This is my mdp file for pulling: Is the pulled distance greater than half of the box vector in z? That's the only reason I can see things going haywire. With simple distance geometry, there are some limitations like that. Rotation itself is not a problem. You're separating two objects; there's no clear orientation dependence in that case. -Justin As you can see the pulling distance is 5 nm. The lenght of my box in Z direction is 12 nm. Thus, it is not. Any other ideas? No, it's not. You're doing 1 ns of simulation (contrary to the comment, 50 * 0.002 = 1000 ps) and thus pulling at 10 nm per ns you are trying to pull 10 nm in this simulation. So once your molecule hits 6 nm (right at about 600 ps from your plot and the expected pull rate), you're violating the criterion that I suspected. -Justin title = Umbrella pulling simulation ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 50; 500 ps nstcomm = 10 ; Output parameters nstxout = 5000 ; every 10 ps nstvout = 5000 nstfout = 500 nstxtcout = 500 ; every 1 ps nstenergy = 500 ; Bond parameters constraint_algorithm= lincs constraints = all-bonds continuation= yes ; continuing from NPT ; Single-range cutoff scheme nstlist = 5 ns_type = grid rlist = 0.9 rcoulomb= 0.9 rvdw= 0.9 ; PME electrostatics parameters coulombtype = PME fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 optimize_fft= yes ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc_grps = Protein_LIG Water_and_ions ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 298 298 ; reference temperature, one for each group, in K ; Pressure coupling is on Pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 1.0 compressibility = 4.5e-5 ref_p = 1.0 ; Generate velocities is off gen_vel = no ; Periodic boundary conditions are on in all directions pbc = xyz ; Long-range dispersion correction DispCorr= EnerPres ; Pull code pull= umbrella pull_geometry = distance ; simple distance increase pull_dim= N N Y pull_start = yes ; define initial COM distance 0 pull_ngroups= 1 pull_group0 = Protein pull_group1 = LIG182 pull_rate1 = 0.01 ; 0.01 nm per ps = 10 nm per ns pull_k1 = 200 ; kJ mol^-1 nm^-2 Would you recommend position restrained of the
[gmx-users] pull-code
Yes you are right, should be between 0 and 0. Do you have a window for a distance equal 0? This window should behave similar to the RDF-analysis. Because there are no directions. Or to reformulate the problem. We make an umbrella window for a distance of 1. If particle stays there everything is fine. If particle moves to 0, it should be also fine (particle sees a force of k*1). If paticle moves to -1, it should see a force of k*2, but instead, the distance is 0 - no force. If you have the umbrella window centered at 0 this problem vansihes - if particle move it sees always a force. But one thing gives me headaches. I don't have this problem in my pulling simulation, because the distance between my reference and pulled group can not become zero: But concerning the reaction coordinate it will have a similar flaw like the RDF i think: It doesn't matter in wish direction the particle moves (left or right) due to the distance we would always say it moves along the reaction coordinate. In reality it moves sometimes in the negative direction of the reaction coordinate, but we always say it's a positve distance - so positive value on the reaction coordinate. For an isotropic system this would not matter, but for system which we have a anisotrop reaction coordinate it should matter. Greetings Thomas Am 17.02.2012 17:07, schrieb gmx-users-requ...@gromacs.org: 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 value0 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
[gmx-users] pull-code
Yes, the pulling is in all direction. But i tried to explain it in one dimension. Lol, send last mail only to Gavin. The real problem (i think) is the following: We have our pulled group 1nm right from the reference group, and the umbrella potential is also centered the (we used pull_start=yes, or pull_init=1). The umbrella force is zero. Now we move the particle 2nm to the left, to we have the same distance, but into the left direction. Since the distance is the same, the umbrella force is zero. On its way the umbrella force increases, maximum when the pulled group is on top of the reference group, and then decreases back to zero. But we wanted a potential, where it doesn't matter if the particle moves to the left or right, force should always increase. But with pull_geometry=distance we can't get it. If the distances between reference and pulled group are big enough we don't run into this problem. -- Hello, I was following your argument but then doesnt pull_dim=Y Y Y mean it pulls in all directions at once? I use pull_dim=N N Y , or just 1 pull direction and it works. I might be wrong, but would be interesting if you got it to work like that for a small molecule. Stephan Watkins Original-Nachricht Datum: Fri, 17 Feb 2012 16:34:22 +0100 Von: Thomas Schlesier schlesi at uni-mainz.de An: gmx-users at gromacs.org Betreff: [gmx-users] pull-code 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 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@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: Npt equilibration of the membrane-mimicking CCl4 layer
AFAIK vdwradii.dat is only used for *genbox*. For the actual simulation the force field paramters of CCl4 will be used. One thing which you could check is what is the compressibility of CCl4 (the value you use reminds me as that of water) and try with this. I do not know if the is a protocol for what to do if you have an biphasic system with different compressibilities. But this is only an idea, do not if it would solve the problem. greetings thomas Mark, I'm using exact all parameters wich I found in different experimental work. By the way reducing of integration step to 1fs provide me with better equilibration of the Ccl4 system ( I've being obtained stabile system during 3 ns) but I had a problems during ntp equilibration when I inserted test peptide into that pre-equilibrated Ccl4 system and made two surrounded layer of water. My system was quickly eqxanded on Z-dimension and slightly shrinked on X. I think that such problem could be due to some problems with the vdw radius value for CCl4. E.g I didnt find this value in the vdwradii.dat file. James 2012/2/16 Mark Abraham mark.abra...@anu.edu.au On 16/02/2012 1:45 AM, James Starlight wrote: Mark, I've used that dimensions in accordance to some literature where the same membrane-mimicking simulation were performed. I've tried to rise cutoffs Don't, that breaks your model physics and makes it even more likely you will encounter problems with the system dimensions becoming too small for the cut-off! and dicrease integration step but my system have been stil crashed during npt. I'm using pcoupl= Parrinello-Rahman wich I've found in the KALP tutorial because I have not found the same npt example file in the Biphastic tutorial :) So you're following some other work and not copying their equilibration protocol and/or model physics? Could you advise me another p_coup algorithm for my Ccl4 system? There's only two choices available. Manual 3.4.9 specifically warns against one of them for equilibration. What is there to say? You should be sure to construct a simple case and get the model physics validated. For the moment, forget about all the stuff where you were struggling to insert more CCl4 into a box with CCl4 (probably creating a far-from-equilibrium starting configuration). Don't try to learn to run on stilts while shaving. Learn to shave, then to walk on stilts, then to run, then start combining them. gmx-users Digest, Vol 94, Issue 95 Mark James -- Forwarded message -- From: Mark Abraham mark.abra...@anu.edu.au Date: 2012/2/15 Subject: Re: [gmx-users] Npt equilibration of the membrane-mimicking CCl4 layer To: Discussion list for GROMACS users gmx-users@gromacs.org On 15/02/2012 4:45 PM, James Starlight wrote: Mark, due to hight density the volume of my system have been slightly increased and during NPT phase I've obtained error Fatal error: One of the box vectors has become shorter than twice the cut-off length or box_yy-|box_zy| or box_zz has become smaller than the cut-off. I'm using 0.9 for electrostatic and 1.4 for vdw cutofs and the dimensions of my box was 6.5 3 3 on the initial step and 6.6 3 3.3 before crush :) I want prevent such expansion of my system by increasing of pressure and/ or compressibility but I have not found exact sollution yet. Your system is dangerously small for those cut-offs if your initial density is not correct for your model physics. Your y and z dimensions only just contain a full cut-off sphere. You should also make sure you are following the advice about choice of P-coupling algorithm in manual 3.4.9, and consider using a very small integration time step. I remain unconvinced by this thread that you have generated a starting configuration that does not have atomic clashes. Mark James 2012/2/14 Mark Abraham mark.abra...@anu.edu.au On 14/02/2012 11:01 PM, James Starlight wrote: This also was solved by the some extra minimisation steps. I've forced with another problem :D During npt equilibration my system have slightly expanded so my desired volume and density were perturbed. I've noticed the below options in npt wich could help me ref_p= 1 1 compressibility = 4.5e-5 i'm using this compressibility value because I'm modelling the lipid-like environment so I think that I must increase pressure. Could you remind me the dependence of pressure from density and volume for liquids ? :) Your forcefield, simulation cell contents and .mdp settings will determine the equilibrium density. Whether you need to do anything depends on whether you've made a statistically significant post-equilibration measurement of your average density. Haphazardly increasing the reference pressure for the coupling will
[gmx-users] Npt equilibration of the membrane-mimicking CCl4 layer
I assume that you energy minimisd the system, but still have atomic clashes? One thing which helped me in a similar case, was a short simulation at low temperature with a really small timestep (about 3-5 magnitudes smaller than the normal timestep). With this the atoms which clashes move away from each other, but due to the very small timestep, they are not fast as rockets and shouldn't lead to an exploding system. Then when most clashes are resolved you can continue with a normal equilibration. Greetings Thomas It seems that I've fixed that problem by reduce vdv radii for Cl during defining of my box Eventually I've obtained box with the desired density than I've delete vdvradii.dat for my wor dir by when I've launched equilibration I've oibtained Fatal error: Too many LINCS warnings (1598) If you know what you are doing you can adjust the lincs warning threshold in your mdp file I've never seen this before I'm using 1.o cutoff for pme and 1.4 for vdv my LINKS parameters are ; Bond parameters continuation= no; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints= all-bonds; all bonds (even heavy atom-H bonds) constrained lincs_iter= 1; accuracy of LINCS lincs_order= 4; also related to accuracy How I could solve it? James 2012/2/14 James Starlight jmsstarli...@gmail.com Mark, I've checked only density value with 500 molecules Ccl4 I have density that is twisely less that I need ( in accordance to the literature ). Also I've checked my box visually and found that the box is not properly tightly packed so I dont know why genbox didnt add some extra mollecules :( In other words I wounder to know if there is any way to add some extra molecules to the pre defined box to make my system more tighly packed ( to short distance between existing molecules and place new ones in the new space ) ? James 2012/2/14 Mark Abraham mark.abra...@anu.edu.au On 14/02/2012 4:57 PM, James Starlight wrote: Justin, Firstly I've created the box of desired size with only 500 molecules ( I need 1000) Than I've tried to add extra 200 molecules by means of Genbox genbox -cp super_box.gro -ci Ccl4.gro -nmol 200 -o new_solv.gro but no molecules have been added Added 0 molecules (out of 200 requested) of Cl4 ... then there are no gaps large enough to insert your molecules. Either make gaps, or check out genbox -h for advice on defining the radii. also I've tried genbox -cp super_box.gro -cp Ccl4.gro -nmol 200 -o new_solv.gro Two -cp options is not what you want, and -nmol probably only works with -ci. but system were crashed with message Reading solute configuration God Rules Over Mankind, Animals, Cosmos and Such Containing 2500 atoms in 500 residues Initialising van der waals distances... WARNING: masses and atomic (Van der Waals) radii will be determined based on residue and atom names. These numbers can deviate from the correct mass and radius of the atom type. Reading solvent configuration God Rules Over Mankind, Animals, Cosmos and Such solvent configuration contains 5 atoms in 1 residues Is there any ways to add extra mollecules to the pre defined box ? Yes - but there has to be room for them. Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/**mailman/listinfo/gmx-usershttp://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/** Support/Mailing_Lists/Searchhttp://www.gromacs.org/Support/Mailing_Lists/Searchbefore 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_Listshttp://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] Npt equilibration of the membrane-mimicking, CCl4 layer
Hi, i second Justins seond idea (creating a small box of equilibrated CCl4 and then fill the simulation box via the -cs option). Depending if you have other molecules in your system, make the simulation box a little bit bigger, because you will get some holes. In the subsequent NPT simulation these holes will vanish and your box will shrink. The first way has the problem, that if you have holes between some CCl4 molecules, no new molecule would fit there. Then you could make a short equilibration, so that the holes gather to bigger holes and fill them again with the -ci option. If you're lucky, you will reach after a very long time the required density. (In my case i used mesitylene as a solvent and after day of the automated iteration (try to put molecules into the box, equilibrate, try to put more molecules into the box, ...) i dismissed this approach.) greetings thomas James Starlight wrote: Justin, 2012/2/6 Justin A. Lemkuljalem...@vt.edumailto:jalem...@vt.edu Some simple calculations using the desired density and the box dimensions (to get the volume) will tell you exactly how many molecules you need. If you only suppose you've got a reasonable number, there are better ways to be sure ;) In accordance to my calculations I need in ~ 1000 CCl4 to obtain 1.600 g/cm^-3 density with my box dimensins By the way way I try to define such system by genbox -ci ccl4.gro -nmol 1000 -box 8.6 6.5 3 -o ccl4_box.gro I've obtain memory error :( Is there any other way to define such system with the desired size and n_mollecules with the less computation demands? Two options: 1. Add them incrementally to the box (100 at a time or so, until the desired value is reached) 2. Build a small box of solvent (a few hundred molecules), then use genbox -cs instead of -ci to use the small solvent box to create the one of desired size. -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] Re: Is there a way to omit particles with q=0, from Coulomb-/PME-calculations?
On 17/01/2012 4:55 AM, Thomas Schlesier wrote: Dear all, Is there a way to omit particles with zero charge from calculations for Coulomb-interactions or PME? In my calculations i want to coarse-grain my solvent, but the solute should be still represented by atoms. In doing so the solvent-molecules have a zero charge. I noticed that for a simulation with only the CG-solvent significant time was spent for the PME-part of the simulation. If i would simulate the complete system (atomic solute + coarse-grained solvent), i would save only time for the reduced number of particles (compared to atomistic solvent). But if i could omit the zero-charge solvent from the Coulomb-/PME-part, it would save much additional time. Is there an easy way for the omission, or would one have to hack the code? If the latter is true, how hard would it be and where do i have to look? (First idea would be to create an index-file group with all non-zero-charged particles and then run in the loops needed for Coulomb/PME only over this subset of particles.) I have only experience with Fortran and not with C++. Only other solution which comes to my mind would be to use plain cut-offs for the Coulomb-part. This would save time required for doing PME but will in turn cost time for the calculations of zeros (Coulomb-interaction for the CG-solvent). But more importantly would introduce artifacts from the plain cut-off :( Particles with zero charge are not included in neighbour lists used for calculating Coulomb interactions. The statistics in the M E G A - F L O P S A C C O U N T I N G section of the .log file will show that there is significant use of loops that do not have Coul component. So already these have no effect on half of the PME calculation. I don't know whether the grid part is similarly optimized, but you can test this yourself by comparing timing of runs with and without charged solvent. Mark Ok, i will test this. But here is the data i obtained for two simulations, one with plain cut-off and the other with PME. As one sees the simulation with plain cut-offs is much faster (by a factor of 6). --- With PME: M E G A - F L O P S A C C O U N T I N G RF=Reaction-Field FE=Free Energy SCFE=Soft-Core/Free Energy T=TabulatedW3=SPC/TIP3pW4=TIP4p (single or pairs) NF=No Forces Computing: M-Number M-Flops % Flops --- VdW(T) 1132.029152 61129.574 0.1 Outer nonbonded loop1020.997718 10209.977 0.0 Calc Weights 16725.001338 602100.048 0.6 Spread Q Bspline 356800.028544 713600.057 0.7 Gather F Bspline 356800.028544 4281600.343 4.4 3D-FFT 9936400.79491279491206.35981.6 Solve PME 18.0144001152.92211.8 NS-Pairs2210.718786 46425.095 0.0 Reset In Box1115.003345.000 0.0 CG-CoM 1115.0004463345.001 0.0 Virial 7825.000626 140850.011 0.1 Ext.ens. Update 5575.000446 301050.024 0.3 Stop-CM 5575.000446 55750.004 0.1 Calc-Ekin 5575.000892 150525.024 0.2 --- Total 97381137.440 100.0 --- D O M A I N D E C O M P O S I T I O N S T A T I S T I C S av. #atoms communicated per step for force: 2 x 94.1 Average load imbalance: 10.7 % Part of the total run time spent waiting due to load imbalance: 0.1 % R E A L C Y C L E A N D T I M E A C C O U N T I N G Computing: Nodes Number G-CyclesSeconds % --- Domain decomp. 4250 903.835 308.1 1.8 Comm. coord. 4 1251 321.930 109.7 0.6 Neighbor search4251 1955.330 666.5 3.8 Force 4 1251 696.668 237.5 1.4 Wait + Comm. F 4 1251 384.107 130.9 0.7 PME mesh 4 125143854.81814948.285.3 Write traj.4 50011.4890.5 0.0 Update 4 1251 1137.630 387.8 2.2 Comm. energies 4 1251 1074.541 366.3 2.1 Rest 41093.194 372.6 2.1 --- Total
[gmx-users] Is there a way to omit particles with, q=0, from Coulomb-/PME-calculations?
But would there be a way to optimize it further? In my real simulation i would have a charged solute and the uncharged solvent (both have nearly the same number of particles). If i could omit the uncharged solvent from the long-ranged coulomb-calculation (PME) it would save much time. Or is there a reason that some of the PME stuff is also calculated for uncharged particles? (Ok, i know that this is a rather specical system, in so far that in most md-simulations the number of uncharged particles is negligible.) Would it be probably better to move the question to the developer-list? Greetings Thomas On 17/01/2012 7:32 PM, Thomas Schlesier wrote: On 17/01/2012 4:55 AM, Thomas Schlesier wrote: Dear all, Is there a way to omit particles with zero charge from calculations for Coulomb-interactions or PME? In my calculations i want to coarse-grain my solvent, but the solute should be still represented by atoms. In doing so the solvent-molecules have a zero charge. I noticed that for a simulation with only the CG-solvent significant time was spent for the PME-part of the simulation. If i would simulate the complete system (atomic solute + coarse-grained solvent), i would save only time for the reduced number of particles (compared to atomistic solvent). But if i could omit the zero-charge solvent from the Coulomb-/PME-part, it would save much additional time. Is there an easy way for the omission, or would one have to hack the code? If the latter is true, how hard would it be and where do i have to look? (First idea would be to create an index-file group with all non-zero-charged particles and then run in the loops needed for Coulomb/PME only over this subset of particles.) I have only experience with Fortran and not with C++. Only other solution which comes to my mind would be to use plain cut-offs for the Coulomb-part. This would save time required for doing PME but will in turn cost time for the calculations of zeros (Coulomb-interaction for the CG-solvent). But more importantly would introduce artifacts from the plain cut-off :( Particles with zero charge are not included in neighbour lists used for calculating Coulomb interactions. The statistics in the M E G A -F L O P S A C C O U N T I N G section of the .log file will show that there is significant use of loops that do not have Coul component. So already these have no effect on half of the PME calculation. I don't know whether the grid part is similarly optimized, but you can test this yourself by comparing timing of runs with and without charged solvent. Mark Ok, i will test this. But here is the data i obtained for two simulations, one with plain cut-off and the other with PME. As one sees the simulation with plain cut-offs is much faster (by a factor of 6). Yes. I think I have seen this before for PME when (some grid cells) are lacking (many) charged particles. You will see that the nonbonded loops are always VdW(T) for tabulated VdW - you have no charges at all in this system and GROMACS has already optimized its choice of nonbonded loops accordingly. You would see Coul(T) + VdW(T) if your solvent had charge. It's not a meaningful test of the performance of PME vs cut-off, either, because there are no charges. Mark --- With PME: M E G A - F L O P S A C C O U N T I N G RF=Reaction-Field FE=Free Energy SCFE=Soft-Core/Free Energy T=TabulatedW3=SPC/TIP3pW4=TIP4p (single or pairs) NF=No Forces Computing: M-Number M-Flops % Flops --- VdW(T) 1132.029152 61129.574 0.1 Outer nonbonded loop1020.997718 10209.977 0.0 Calc Weights 16725.001338 602100.048 0.6 Spread Q Bspline 356800.028544 713600.057 0.7 Gather F Bspline 356800.028544 4281600.343 4.4 3D-FFT 9936400.79491279491206.35981.6 Solve PME 18.0144001152.92211.8 NS-Pairs2210.718786 46425.095 0.0 Reset In Box1115.003345.000 0.0 CG-CoM 1115.0004463345.001 0.0 Virial 7825.000626 140850.011 0.1 Ext.ens. Update 5575.000446 301050.024 0.3 Stop-CM 5575.000446 55750.004 0.1 Calc-Ekin 5575.000892 150525.024 0.2 --- Total 97381137.440 100.0 --- D O M A I N D E C O M P O S I T I O N S T A T I S T I C S av. #atoms communicated per step
[gmx-users] Re: Is there a way to omit particles with, q=0, from, Coulomb-/PME-calculations?
Thanks Carsten. Now i see the problem. Hi Thomas, Am Jan 17, 2012 um 10:29 AM schrieb Thomas Schlesier: But would there be a way to optimize it further? In my real simulation i would have a charged solute and the uncharged solvent (both have nearly the same number of particles). If i could omit the uncharged solvent from the long-ranged coulomb-calculation (PME) it would save much time. Or is there a reason that some of the PME stuff is also calculated for uncharged particles? For PME you need the Fourier-transformed charge grid and you get back the potential grid from which you interpolate the forces on the charged atoms. The charges are spread each on typically 4x4x4 (=PME order) grid points, and in this spreading only charged atoms will take part. So the spreading part (and also the force interpolation part) will become faster with less charges. However, the rest of PME (the Fourier transforms and calculations in reciprocal space) are unaffected by the number of charges. For this only the size of the whole PME grid matters. You could try to lower the number of PME grid points (enlarge fourierspacing) and at the same time enhance the PME order (to 6 for example) to keep a comparable force accuracy. You could also try to shift more load to real space, which will also lower the number of PME grid points (g_tune_pme can do that for you). But I am not shure that you can get large performance benefits from that. Best, Carsten (Ok, i know that this is a rather specical system, in so far that in most md-simulations the number of uncharged particles is negligible.) Would it be probably better to move the question to the developer-list? Greetings Thomas On 17/01/2012 7:32 PM, Thomas Schlesier wrote: On 17/01/2012 4:55 AM, Thomas Schlesier wrote: Dear all, Is there a way to omit particles with zero charge from calculations for Coulomb-interactions or PME? In my calculations i want to coarse-grain my solvent, but the solute should be still represented by atoms. In doing so the solvent-molecules have a zero charge. I noticed that for a simulation with only the CG-solvent significant time was spent for the PME-part of the simulation. If i would simulate the complete system (atomic solute + coarse-grained solvent), i would save only time for the reduced number of particles (compared to atomistic solvent). But if i could omit the zero-charge solvent from the Coulomb-/PME-part, it would save much additional time. Is there an easy way for the omission, or would one have to hack the code? If the latter is true, how hard would it be and where do i have to look? (First idea would be to create an index-file group with all non-zero-charged particles and then run in the loops needed for Coulomb/PME only over this subset of particles.) I have only experience with Fortran and not with C++. Only other solution which comes to my mind would be to use plain cut-offs for the Coulomb-part. This would save time required for doing PME but will in turn cost time for the calculations of zeros (Coulomb-interaction for the CG-solvent). But more importantly would introduce artifacts from the plain cut-off :( Particles with zero charge are not included in neighbour lists used for calculating Coulomb interactions. The statistics in the M E G A -F L O P S A C C O U N T I N G section of the .log file will show that there is significant use of loops that do not have Coul component. So already these have no effect on half of the PME calculation. I don't know whether the grid part is similarly optimized, but you can test this yourself by comparing timing of runs with and without charged solvent. Mark Ok, i will test this. But here is the data i obtained for two simulations, one with plain cut-off and the other with PME. As one sees the simulation with plain cut-offs is much faster (by a factor of 6). Yes. I think I have seen this before for PME when (some grid cells) are lacking (many) charged particles. You will see that the nonbonded loops are always VdW(T) for tabulated VdW - you have no charges at all in this system and GROMACS has already optimized its choice of nonbonded loops accordingly. You would see Coul(T) + VdW(T) if your solvent had charge. It's not a meaningful test of the performance of PME vs cut-off, either, because there are no charges. Mark --- With PME: M E G A - F L O P S A C C O U N T I N G RF=Reaction-Field FE=Free Energy SCFE=Soft-Core/Free Energy T=TabulatedW3=SPC/TIP3pW4=TIP4p (single or pairs) NF=No Forces Computing: M-Number M-Flops % Flops --- VdW(T) 1132.029152 61129.574 0.1 Outer nonbonded loop1020.997718 10209.977 0.0 Calc Weights 16725.001338 602100.048
[gmx-users] table-potential, why table from r=0-r_c+1?
Dear all, what is the reason, that the tabulated potential must go till r_c+1 (r_c = cut-off radius) and not only up to r_c? I think we only calculate the interactions till r_c and truncate the rest. So everything behind r_c would be redundant information (due to the truncation it would be set to 0). If this would be right, i could fill everything from r_c to r_c+1 with 0's, but this sounds to easy. Greetings Thomas -- 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] Is there a way to omit particles with q=0 from Coulomb-/PME-calculations?
Dear all, Is there a way to omit particles with zero charge from calculations for Coulomb-interactions or PME? In my calculations i want to coarse-grain my solvent, but the solute should be still represented by atoms. In doing so the solvent-molecules have a zero charge. I noticed that for a simulation with only the CG-solvent significant time was spent for the PME-part of the simulation. If i would simulate the complete system (atomic solute + coarse-grained solvent), i would save only time for the reduced number of particles (compared to atomistic solvent). But if i could omit the zero-charge solvent from the Coulomb-/PME-part, it would save much additional time. Is there an easy way for the omission, or would one have to hack the code? If the latter is true, how hard would it be and where do i have to look? (First idea would be to create an index-file group with all non-zero-charged particles and then run in the loops needed for Coulomb/PME only over this subset of particles.) I have only experience with Fortran and not with C++. Only other solution which comes to my mind would be to use plain cut-offs for the Coulomb-part. This would save time required for doing PME but will in turn cost time for the calculations of zeros (Coulomb-interaction for the CG-solvent). But more importantly would introduce artifacts from the plain cut-off :( Hope anyone could help. Greetings Thomas -- 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] Coarse-graining and cut-offs
Dear all, first of all, sorry to this rather conceptional question, which is not totally to GROMOACS related. But probably anyone of you can help. In my simulations I use mesitylene as a solvent. In future i want to coarse-grain the full atomic mesitylene to an effective one-particle. For this i did a NPT-simulation to obtain the RDF (with regard to the COM of the molecules) from which i want to calculate an effective interaction potential by iterative boltzmann. The full-atomistic simulations (slovate + solvent) employed a cut-off of 1.4 nm. At this distance (1.4 nm) the RDF shows it second minimum (RDF=0.97) and after 1.8 nm the RDF is around 1.0 +/- 0.15 (I got the RDF only up to a distance of 2.3 nm). Now my question is how long should be the table for the effective potential (i.e. maximum distance for an interaction)? If i would use 1.4 nm, i would have a jump in the interaction potential (which is normal due to the truncation of the cut-off). But for coulomb-interaction one would have PME or other stuff which would correct the artifects to some degree. Or would it be better to use a long interaction table (interactions at longer distances) and truncate the table at a distance after 1.4 nm where the RDF is equal to 1? Hope anybody can give some insight. Greetings Thomas -- 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
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 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
Re: [gmx-users] RDF(PMF) and Umbrella sampling
The histograms are really crowded, it would be better to plot only the black one and probably the red one, to see the it better. Two ideas which can probably solve the problem (ok, both assume, that the host has a certain shape). You said you investigate guest-insertion to a host. I would think that the guest molecule can't bind from all directions to the host. If the host is 'U' shaped then, the guest can come only from above. Idea 1: Take some atoms at the lower part from 'U' as the reference. So you do not run into the problem with with a zero ref-distance. Bad thing is, that you would start from all over again. Idea 2: we make some tricks to get negative distances. I think the big problem is that from pull_geometry=distance, you get always positive distances (really, only a idea, WHAM could make problems, i don't really know how WHAM and the GROMACS implementation of it work): Instead of using pull_geometry=distance, use pull_geometry=position. Define a plane at zero ref-distance, which is roughly orthogonal to the path which the guest can enter. In my example 'U-' - would be the plane. From the pullx.xvg you can calculate the pull-ref distance (which is always positive). Then you can add a sign to the distances, depending if the guest is above or below the plane. (Probably it is better to take only the part of the pull-ref distance which is orthogonal to the plane. Then you should also modify the forces in the same way.) From this distances you can also calculate pullf.xvg. Now you must create a *.tpr for this windows with pull_geometry=distance, instead of use pull_geometry=position (WHAM doesn't like position, as far as i know). Now start WHAM and pray that it will work. All we have done is, that we have now negative distances in pullx.xvg. I don't know if the fact that we now have negative numbers makes any problem for WHAM (from the *.tpr we would expect that we have only positive numbers). Wish you luck. greetings thomas 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.
[gmx-users] RDF: mol_com and atom at the same time
Dear all, i would like to know if it is possible to get the rdf between the center of mass of a molecule and individual atoms of said molecule? In my case i have mesitylene and i would like to calculate the RDF between the COM and the Methyl-carbon-atoms. My problem is that i would need the options '-rdf mol_com' and '-rdf atom' at the same time. The only idea i have is to use '-yescom' and make the analysis for each molecule individually and then average overall the resulting RDFs. But probably there is an easier solution? greetings thomas -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] How to list OPLS parameters
Hi all, for a publication i want to list the used OPLS parameters for the investigated molecules. In GROMACS all atomtypes are uniquely defined by the atomtype `opls_xyz`. And from the atomtypes one can deduct the bonded parameters. So it is sufficient to list only the atomtypes and probably charges. My question is now, is this definition common knowledge, or only GROMACS-intern? From the header of `ffoplsaa.atp` one sees that the first 65 atomtypes are for the OPLS-UA force field, and the names of the atomtypes are the same as in the paper which is mentioned. Concerning the OPLS-AA force field, is there the GROMACS numbering scheme also common knowledge? Greetings Thomas -- 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:umbrella sampling with pull=constraint
If you use constraints it would not be umbrella sampling, where you need to sample around a restraint structure to get the histograms for WHAM or another analysis-technic. So if you want to do umbrella sampling either make to box bigger and/or make the restraints harder. But you can also use constraints instead of restraints. But it would require another analysis-technic. Think it's called thermodynamic integration. In the end you also sample many different windows with varying distances, but instead of using WHAM you integrate the constraint force to obtain the PMF. http://www.sciencedirect.com/science/article/pii/000926149190259C One important thing to note is, that you will need more windows compared to umbrella-sampling, since in each window you sample only on distance and not around one distance (restrains allow for a change of the distance). Another thing to note is the pull_dim. For determining the PMF of an end-toend distance of similar i would turn all three dimensions to 'yes'. Imagine to parallel sheets of paper, with z perpendicular to the papers. With your setup you would only fix the z-distance between both sheets, but they can move freely in the xy-plane. Hope that helps. greetings thomas Hello, I have done umbrella sampling with pull=umbrella and I found that the pulling group has high fluctuations and sometimes moving out of the periodic box. I think that the harmonic potential is not properly applied and thus the pulling group is not retained with the specified COM distance between the reference and pulling group. Can we use pull=constraint option to retain the pulling group within the COM distance? my pulling code is as bellow with restrain at 0.5nm distance from ref group. I just want to get some idea about this pull code modification. pull= constraint pull_geometry= distance pull_dim= N N Y pull_start = no pull_ngroups= 1 pull_group0 = Chain_A pull_group1 = Chain_B pull_init1 = 0.50 pull_rate1 = 0.0 pull_k1 = 1000 pull_nstxout= 500 pull_nstfout= 500 Regards, Vijay. -- 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] g_dist error
never used the -dist option, but i think you have here a missunderstanding: t: 275 20230 SOL 62618 OW 0.341434 (nm) i think this means: at time=275 atom 62618 (which is a OW) from residue 20230 (which is a SOL) is 0.341434nm away from your protein atom. the x SOL y OW means not the distance between some SOl-atom and OW, but that the OW atom is from the residue x from the SOL group!!! And when am already specifying only one atom from protein ie say 1322. why do i get this kind of output- t: 275 20230 SOL 62618 OW 0.341434 (nm) t: 275 22019 SOL 67985 OW 0.171584 (nm) t: 276 10768 SOL 34232 OW 0.303328 (nm) t: 276 20230 SOL 62618 OW 0.325176 (nm) ... you get this output, because the programm assumes that you remember which was group 1 and so it does not need to write it in every line. since i never used the -dist option, the stuff above could be false, so short thing you can do is and look if atom 62618 is in residue 20230, and then the above should be right. greetings thomas Date: Tue, 13 Sep 2011 13:32:26 +0530 From: aiswarya pawaraiswarya.pa...@gmail.com Subject: Re: [gmx-users] g_dist error To: Discussion list for GROMACS usersgmx-users@gromacs.org Message-ID: CAEa6cRA+7W5Hed_KHtPN3cqmPD8xyE5ZJrS-HSap6-2C6h2=9...@mail.gmail.com Content-Type: text/plain; charset=iso-8859-1 g_dist -f md.xtc -s md.tpr -dist 0.5 -e 500 if the index name not given it takes the default index file, so there isnt any wrong in selecting the atoms. On Tue, Sep 13, 2011 at 1:24 PM, Mark Abrahammark.abra...@anu.edu.auwrote: On 13/09/2011 4:14 PM, aiswarya pawar wrote: Mark, the command line goes like this- g_dist -f md.xtc -s md.tpr -dist 0.5 -e 500 This command line is not using an index file. The index groups defined for the grompp that produced the .tpr are being used (which may be the default ones, depending what you did). Please copy and show the interactive input you made to g_dist after it showed the groups it knew about. the index file has group1- a_1322 ( this is just a single atom from a protein. its in sidechain) group2- a_OW ( this is water atoms) Your output is listing the time of the frame, and the residue number, residue name, atom number, and atom name of the matching atom. Apparently a water molecule can sometimes be closer than 0.2nm, and sometimes not. Mark now as per the utility it should give me and output as- t:1 1322 a 54119 OW 0.389 but am getting something different On Tue, Sep 13, 2011 at 11:24 AM, Mark Abrahammark.abra...@anu.edu.auwrote: On 13/09/2011 3:40 PM, aiswarya pawar wrote: Hi Mark, The -dist option says- print all the atoms in group 2 that are closer than a certain distance to the center of mass of group 1. That means it should give me the distance from OW to protein atom. If you choose correct groups that are correctly defined with respect to your trajectory and use a large enough distance cutoff. And when am already specifying only one atom from protein ie say 1322. why do i get this kind of output- We can't tell. We don't know what your command line is, what's in your simulation system, the contents of your index groups, or which groups you've selected for which role. Clearly something is not working properly, and our time constraints mean that we're going to assume you've done something wrong, in the absence of evidence to the contrary. Mark t: 275 20230 SOL 62618 OW 0.341434 (nm) t: 275 22019 SOL 67985 OW 0.171584 (nm) t: 276 10768 SOL 34232 OW 0.303328 (nm) t: 276 20230 SOL 62618 OW 0.325176 (nm) t: 276 22019 SOL 67985 OW 0.187259 (nm) t: 277 10768 SOL 34232 OW 0.306008 (nm) t: 277 20230 SOL 62618 OW 0.326195 (nm) t: 277 22019 SOL 67985 OW 0.181089 (nm) t: 278 10768 SOL 34232 OW 0.274507 (nm) t: 278 22019 SOL 67985 OW 0.292652 (nm) t: 279 10618 SOL 33782 OW 0.319922 (nm) t: 280 10618 SOL 33782 OW 0.330082 (nm) t: 280 22019 SOL 67985 OW 0.330203 (nm) t: 281 8273 SOL 26747 OW 0.278731 (nm) t: 281 10618 SOL 33782 OW 0.325434 (nm) t: 281 11535 SOL 36533 OW 0.200327 (nm) t: 281 17036 SOL 53036 OW 0.343946 (nm) t: 282 8273 SOL 26747 OW 0.256558 (nm) t: 282 11535 SOL 36533 OW 0.327147 (nm) t: 283 8273 SOL 26747 OW 0.165061 (nm) t: 283 10618 SOL 33782 OW 0.306578 (nm) t: 283 17191 SOL 53501 OW 0.333075 (nm) t: 284 8273 SOL 26747 OW 0.19427 (nm) t: 284 10618 SOL 33782 OW 0.321927 (nm) t: 284 17191 SOL 53501 OW 0.30832 (nm) On Tue, Sep 13, 2011 at 10:40 AM, Mark Abrahammark.abra...@anu.edu.auwrote: On 13/09/2011 2:27 PM, aiswarya.pa...@gmail.com wrote: Iam -dist option because I need the distance between two groups That is not what g_dist -dist does. Please read g_dist -h. excluding -dist gives me X,Y,Z output which I don't want. And other output which you do, but you have to use -o to get it. Read g_dist -h. And am not specifying an -o. You need to specify -o to achieve your purpose. However, as I
[gmx-users] mirroring a trajectorie
Hi all, is it possible to mirror a trajectorie? I have done pulling simulations, where i first pulled two molecules apart and later used the pulled them together (starting form the last frame of the pulling-(apart)-simulation). Now want to calculate the rmsd of structures for the path. So i would like to compared first frame of pull-apart with last frame of pull-together ... Only way i could think of, would be to write the .xtc into individual .gro-files, then order them and finally put them together to the mirrored .xtc This should work, but probably there is an easier way of doing this? Greetings Thomas -- 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: pull forces
Message: 5 Date: Tue, 15 Mar 2011 07:17:28 -0700 (PDT) From: Michael Brunsteinermbx0...@yahoo.com Subject: [gmx-users] pull forces To: gmx usersgmx-users@gromacs.org Message-ID:613152.30411...@web120517.mail.ne1.yahoo.com Content-Type: text/plain; charset=us-ascii Dear all, does anybody know what the forces that are saved in the f*xvg files from the pull-code actually are? are these: 1) the (sum of the) forces due to the normal non-bonded interactions in the system acting on the ref and the the full groups. 2) only the forces from the harmonic restraint. 3) the sum of both? if you mean the pullf.xvg it's 2. only the force from the harmonic restraint (if you're using 'pull=umbrella') i just did a test writing out the forces from the trajectory file, and looking at the relevant forces (the ones in the z-dimension in my case, as i use pulldim N N Y and constrain the groups to stay put in the x and y dimensions) it seems as if the forces on the pull group are oscillating as expected, but the forces on the reference group are close to, but not exactly(!) zero. i am not sure how to interpret this ... cheers Michael -- 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] problems with randoms seeds
hi all, i'm trying out GROMACS 4.5.3 i'm simulating (sd-integrator) a small rna hairpin in vacuum, without pbc for 100ps. System has about 380 atoms. If i do the simulation a second time (grompp + mdrun) i get identical results. first i had only *gen_seed = -1* but latter i set also *ld_seed = -1*, but nothing changes. one thing which really confuses me is the ld_seed: output from grompp: Setting the LD random seed to 1870 md.log: ld_seed = 1993 position restraints are only applied to 3atoms. has anyone an idea what happens? greetings thomas here is the mdp-file: OPLS Lysozyme NVT equilibration define = -DPOSRES ; position restrain the protein ; Run parameters integrator = sd; leap-frog integrator nsteps = 5 ; 2 * 5 = 100 ps dt = 0.002 ; 2 fs comm_mode = angular ; Output control nstxout = 100 ; save coordinates every 0.2 ps nstvout = 100 ; save velocities every 0.2 ps nstenergy = 100 ; save energies every 0.2 ps nstlog = 100 ; update log file every 0.2 ps ; Bond parameters continuation= no; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = simple; search neighboring grid cells nstlist = 0 ; 10 fs rlist = 0.0 ; short-range neighborlist cutoff (in nm) rcoulomb= 0.0 ; short-range electrostatic cutoff (in nm) rvdw= 0.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = cut-off ; Particle Mesh Ewald for long-range electrostatics ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = System; two coupling groups - more accurate tau_t = 2.0 ; time constant, in ps ref_t = 300 ; reference temperature, one for each group, in K ; Pressure coupling is off pcoupl = no; no pressure coupling in NVT ; Periodic boundary conditions pbc = no; 3-D PBC ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 300 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed ld_seed = -1 -- 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: problems with randoms seeds
ok, problem solved... used same .tpr file all the time should stop working for today :) greetings thomas On 02/22/2011 08:49 PM, Thomas Schlesier wrote: hi all, i'm trying out GROMACS 4.5.3 i'm simulating (sd-integrator) a small rna hairpin in vacuum, without pbc for 100ps. System has about 380 atoms. If i do the simulation a second time (grompp + mdrun) i get identical results. first i had only *gen_seed = -1* but latter i set also *ld_seed = -1*, but nothing changes. one thing which really confuses me is the ld_seed: output from grompp: Setting the LD random seed to 1870 md.log: ld_seed = 1993 position restraints are only applied to 3atoms. has anyone an idea what happens? greetings thomas here is the mdp-file: OPLS Lysozyme NVT equilibration define = -DPOSRES ; position restrain the protein ; Run parameters integrator = sd ; leap-frog integrator nsteps = 5 ; 2 * 5 = 100 ps dt = 0.002 ; 2 fs comm_mode = angular ; Output control nstxout = 100 ; save coordinates every 0.2 ps nstvout = 100 ; save velocities every 0.2 ps nstenergy = 100 ; save energies every 0.2 ps nstlog = 100 ; update log file every 0.2 ps ; Bond parameters continuation = no ; first dynamics run constraint_algorithm = lincs ; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = simple ; search neighboring grid cells nstlist = 0 ; 10 fs rlist = 0.0 ; short-range neighborlist cutoff (in nm) rcoulomb = 0.0 ; short-range electrostatic cutoff (in nm) rvdw = 0.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = cut-off ; Particle Mesh Ewald for long-range electrostatics ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = System ; two coupling groups - more accurate tau_t = 2.0 ; time constant, in ps ref_t = 300 ; reference temperature, one for each group, in K ; Pressure coupling is off pcoupl = no ; no pressure coupling in NVT ; Periodic boundary conditions pbc = no ; 3-D PBC ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp = 300 ; temperature for Maxwell distribution gen_seed = -1 ; generate a random seed ld_seed = -1 -- 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
