Re: [gmx-users] Distance for window spacing in US?
So in umbrella sampling, does it really matter if the vector connecting the COMs of protein and ligand is NOT parallel with the vector of the pulling force (although the pull rate is 0)? I guess as long as the force in each sampling window is along the same direction, we can in theory get potential of mean force along any direction. Thanks, Yun On Wed, Feb 20, 2013 at 1:46 PM, Justin Lemkul jalem...@vt.edu wrote: On 2/20/13 11:01 AM, Yun Shi wrote: Hi Justin, I was not able to find relevant posts in the archive. Any special keyword for searching? I'm afraid I can't think of anything magical off the top of my head, no. But I wonder if it would introduce any additional error to deltaG calculations if I restrain the dissociation pathway to x axis, using pull_vecl = 1 0 0. Would the program then forcing the velocities along y and z axes to be zero? No. If you choose the restraint direction to only coincide with x, then restraint forces are only applied along x. That does not mean there can be no motion along y and z. You have to decide on a suitable dimension or set of dimensions based on the geometry of your system. -Justin -- Justin A. Lemkul, Ph.D. Research Scientist Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Protein unfolded after COM pulling
A general question: In COM pulling (pull two molecules away), does the weakest inter-molecular interaction always break first? Or is it the interaction aligned on the pulling direction break first? In other words, does every part of the molecule feel the pulling force? Or the part close to the COM feel stronger force that the part of molecule that is far away from its COM? Thanks, Yun On Mon, Feb 4, 2013 at 6:35 PM, Justin Lemkul jalem...@vt.edu wrote: On 2/4/13 9:32 PM, Yun Shi wrote: Hi all, I am pulling one monomer of a tetrameric protein away from the other three monomers with pull_k1 = 1 and pull_rate1 = 0.0005. As shown in the attached picture, the green monomer (being pulled) has unfolded while two regions on it are still interacting with the other three monomers at the end of my pulling simulations (13 ns, so COM distance has been pulled 6.5 nm). Attachments to the list are prohibited; if you have files or images to share you must provide a URL. If I want to calculate deltaG of this process, can I use distance restraints on the green monomer being pulled to maintain its folded structure? Otherwise, I will have to construct a larger box or pull even slower than 0.5 nm/ns. Yes, you probably need distance restraints. Complex restraints often fail with domain decomposition, so you may need to use particle decomposition instead. -Justin -- Justin A. Lemkul, Ph.D. Research Scientist Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Many energygrps to output
I guess I will do mdrun -rerun 400 times then. Thanks, Yun On Thu, Feb 7, 2013 at 9:06 PM, Bogdan Costescu bcoste...@gmail.com wrote: On Thu, Feb 7, 2013 at 8:13 PM, Yun Shi yunsh...@gmail.com wrote: So instead of making an index file with 399 groups of each residue in A and typing in rerun.mdp file 400 group names as energygrps = A1 A2 A3 A4 ... A399 B, and issuing g_energy command 399 times, There can be maximum 256 groups defined at once (Chapter 3.3 in manual), so you can't do this anyway... Cheers, Bogdan -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] amber to gromacs error
On Sun, Dec 9, 2012 at 9:46 AM, Albert mailmd2...@gmail.com wrote: hello: I am using the command: acpype.py -p prmtop -x S13.rst to convert Amber system into Gromacs system, but it failed when I try to generate .tpr file: WARNING 1 [file prmtop_GMX.top, line 19]: Too few parameters on line (source file toppush.c, line 300) WARNING 2 [file prmtop_GMX.top, line 22]: Too few parameters on line (source file toppush.c, line 300) Generated 820 of the 820 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 820 of the 820 1-4 parameter combinations --- Program grompp, VERSION 4.5.5 Source code file: toppush.c, line: 1166 Fatal error: Atomtype 3C not found For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors --- This Doesn't Suck, It's a Black Hole ! (K.A. Feenstra) Not sure about too few But it seems you have some atoms with type 3C that is not included (defined) in top files. My ligand contains both glycam FF and GAFF parameters from Amber. It is quite strange when I open the generated .top file because there is even no include amber99SB.ff/forcefiled.itp line Here is the .top file I got: My resulting .top and .gro files worked well without these include ... stuff, and I believe all the atom types and vdw parameters, bond, angle, and dihedral parameters, as well as atomic charges have already been included in the resulting .top file. ; prmtop_GMX.top created by acpype (Rev: 7234) on Sun Dec 9 18:35:28 2012 [ defaults ] ; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ 1 2 yes 0.5 0.8333 [ atomtypes ] ;name bond_type mass charge ptype sigma epsilon Amb N3 N3 0.0 0.0 A 3.25000e-01 7.11280e-01 ; 1.82 0.1700 HH 0.0 0.0 A 1.06908e-01 6.56888e-02 ; 0.60 0.0157 CX CX 0.0 0.0 A 3.39967e-01 4.57730e-01 ; 1.91 0.1094 HP HP 0.0 0.0 A 1.95998e-01 6.56888e-02 ; 1.10 0.0157 CT CT 0.0 0.0 A 3.39967e-01 4.57730e-01 ; 1.91 0.1094 HC HC 0.0 0.0 A 2.64953e-01 6.56888e-02 ; 1.49 0.0157 CC 0.0 0.0 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860 OO 0.0 0.0 A 2.95992e-01 8.78640e-01 ; 1.66 0.2100 NN 0.0 0.0 A 3.25000e-01 7.11280e-01 ; 1.82 0.1700 H1 H1 0.0 0.0 A 2.47135e-01 6.56888e-02 ; 1.39 0.0157 3C 3C 0.0 0.0 A 3.39967e-01 4.57730e-01 ; 1.91 0.1094 OH OH 0.0 0.0 A 3.06647e-01 8.80314e-01 ; 1.72 0.2104 HO HO 0.0 0.0 A 0.0e+00 0.0e+00 ; 0.00 0. 2C 2C 0.0 0.0 A 3.39967e-01 4.57730e-01 ; 1.91 0.1094 SS 0.0 0.0 A 3.56359e-01 1.04600e+00 ; 2.00 0.2500 CA CA 0.0 0.0 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860 HA HA 0.0 0.0 A 2.59964e-01 6.27600e-02 ; 1.46 0.0150 CO CO 0.0 0.0 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860 O2 O2 0.0 0.0 A 2.95992e-01 8.78640e-01 ; 1.66 0.2100 C* C* 0.0 0.0 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860 CW CW 0.0 0.0 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860 H4 H4 0.0 0.0 A 2.51055e-01 6.27600e-02 ; 1.41 0.0150 NA NA 0.0 0.0 A 3.25000e-01 7.11280e-01 ; 1.82 0.1700 CN CN 0.0 0.0 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860 CB CB 0.0 0.0 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860 C8 C8 0.0 0.0 A 3.39967e-01 4.57730e-01 ; 1.91 0.1094 CC CC 0.0 0.0 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860 NB NB 0.0 0.0 A 3.25000e-01 7.11280e-01 ; 1.82 0.1700 CR CR 0.0 0.0 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860 H5 H5 0.0 0.0 A 2.42146e-01 6.27600e-02 ; 1.36 0.0150 N2 N2 0.0 0.0 A 3.25000e-01 7.11280e-01 ; 1.82 0.1700 CV CV 0.0 0.0 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860 SH SH 0.0 0.0 A 3.56359e-01 1.04600e+00 ; 2.00 0.2500 HS HS 0.0 0.0 A 1.06908e-01 6.56888e-02 ; 0.60 0.0157 Ho Ho 0.0 0.0 A 3.56359e-02 1.25520e-01 ; 0.20
Re: [gmx-users] How to avoid adding ions close to ligand
I did hope the ions will move out eventually. But after my ~70ns of conventional MD (with duplicate MD runs and the protein as a dimer with identical sequence), they were still there in the binding site. So I assume it would be much better to start without any salt ions beside my ligand. Could anyone suggest a way around this? Thanks, Yun On Mon, Nov 26, 2012 at 12:39 PM, David van der Spoel sp...@xray.bmc.uu.se wrote: On 2012-11-26 21:28, Yun Shi wrote: Hi everyone, I am doing conventional MD of a protein-ligand system with a mobile loop as part of the binding site. Presumably, the positive Arg side chain on the mobile loop will eventually move towards the negative carboxylic group on my ligand. But I found the addition of NaCl (0.15 M conc.) had some effect on this movement, since the random addition could put Na+ or Cl- ions between the mobile loop and my ligand. I tried generating a index containing only SOL far not close to my ligand, but apparently genion requires a continuous solvent group. So is there any other way to achieve this? Trying different numbers for -seed option seems inefficient and is dependent on luck. Thanks, Yun Maybe your assumption is wrong? Run a long MD simulation and you will find out. -- David van der Spoel, Ph.D., Professor of Biology Dept. of Cell Molec. Biol., Uppsala University. Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. sp...@xray.bmc.uu.sehttp://folding.bmc.uu.se -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] How to delete atoms during MD runs?
Hi everyone, I am running MD of an enzyme containing co-factor FADH2. After the enzyme stabilize, I want to reduce the FADH2 to FAD, which for me is to delete two H atoms and change the topology file and .gro file correspondingly. However, I want to continue the MD run with the state (velocities, temperature, pressure...) right before the reduction (deletion of 2 H atoms), which requires .cpt file. How could I match the .cpt file with the new structure that has 2 H atoms deleted? Thanks for any suggestion or way-around, Yun -- 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] How to delete atoms during MD runs?
Hi Mark, I just want to look into the conformational change of the enzyme after reduction of the co-factor FAD. So instead of feed into a reducing molecule and do some complicated QM simulation, I wonder if this way would be better. I understand that this is not physical. But what could be an alternative way to look at the system evolution upon reduction of the co-factor? Thanks, Yun On Sat, Sep 8, 2012 at 3:17 PM, Mark Abraham mark.abra...@anu.edu.auwrote: On 9/09/2012 8:09 AM, Yun Shi wrote: Hi everyone, I am running MD of an enzyme containing co-factor FADH2. After the enzyme stabilize, I want to reduce the FADH2 to FAD, which for me is to delete two H atoms and change the topology file and .gro file correspondingly. However, I want to continue the MD run with the state (velocities, temperature, pressure...) right before the reduction (deletion of 2 H atoms), which requires .cpt file. How could I match the .cpt file with the new structure that has 2 H atoms deleted? You can't. There's no sense in which the two would be continuous, and as far as I can see little virtue in seeking it. It might be possible to convert the .cpt to a .trr and form the correct subset, projecting the missing velocities properly, etc. which could be given to grompp with a new .top... but what's the point? 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
Re: [gmx-users] rvdw and dispersion correction for amber99sb
In the original paper [Comparison of multiple Amber force fields and development of improved protein backbone parameters], 0.8 nm was used; but I did find recent papers using 1.0 nm for rvdw. And I myself used 1.4 nm, all without obvious abnormality. On Mon, Aug 20, 2012 at 10:13 PM, Mark Abraham mark.abra...@anu.edu.au wrote: On 21/08/2012 2:30 PM, Yun Shi wrote: Hello all, I am simulating protein-ligand complex with amber99sb force field in TIP3P water. What would be a reasonable value rvdw? I saw someone uses 1.0, and I did not find any abnormality when using 1.4 nm for rvdw, but where can I find the right reference? It's cited in the GROMACS manual. Of course, you read it, and reports of similar studies others did with the force field before choosing it,right? :-) And is this force field parameterized with or without dispersion correction? (Better to turn it on or off?) Generally, one should follow the parametrization strategy, unless/until someone shows that observables relevant to yours are unaffected or improved by a change. Mark -- 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 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
Re: [gmx-users] how to use the triplets index file generated by g_hbond
Hi Justin, Thank you for the perl script. But I wonder now why there is a difference between what shows in the hbmap.xpm file and the summary_HBmap.dat generated. For example, I found many occurrences of the second to last H-bond in the 20hbmap.xpm file attached (visualized by GIMP), but in the summary_HBmap.dat file there showed only 0.1 % Exist. of this specific H-bond. Anything wrong with my interpretation? Please see attached the files I used/ generated. Thanks, Yun On Sat, May 5, 2012 at 10:43 AM, Justin A. Lemkul jalem...@vt.edu wrote: On 5/5/12 1:22 PM, Yun Shi wrote: Hello all, I have used g_hbond with -hbn option to generate a .ndx file that has Acceptor - Donor - Hydrogen in each line of the last index group. But I wonder I could I use this index file to monitor each hydrogen bond specified by these triplets along the trajectory? Depending on what you're trying to do, you can map the triplets in hbond.ndx with the hbmap.xpm file to obtain, for instance, the % of total frames that a particular hydrogen bond exists. I've done something along those lines with a script I call plot_hbmap.pl, available from: http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justin/**scripts.htmlhttp://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/scripts.html Please note the usage information in the header of the file. Also, I understand that the -ins in no longer available. But how can we analyze those water-mediated hydrogen bonds? Doing something like this will require careful use of index groups and external scripting. Presumably, if you have two groups (let's call them A and B), if group A forms hydrogen bonds to a certain group of SOL residues, and if group B does the same, with some of those SOL molecules overlapping in some given frames, you can make an argument for a water-mediated hydrogen bond. Determining the existence of the hydrogen bonds is easy with g_hbond and index groups. Determining simultaneous occurrence requires an external script, and can make use of hbmap.xpm (to determine which frames contain the hydrogen bonds) and hbond.ndx (to determine which groups are participating in those hydrogen bonds). -Justin -- ==**== Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.**vt.edu/Pages/Personal/justinhttp://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ==**== -- 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 attachment: 20hbmap.xpm summary_HBmap.dat Description: Binary data mod_hb.ndx Description: Binary data -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] how to use the triplets index file generated by g_hbond
Hello all, I have used g_hbond with -hbn option to generate a .ndx file that has Acceptor - Donor - Hydrogen in each line of the last index group. But I wonder I could I use this index file to monitor each hydrogen bond specified by these triplets along the trajectory? Also, I understand that the -ins in no longer available. But how can we analyze those water-mediated hydrogen bonds? Thanks for any suggestion, Yun -- 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 g_select bound water molecules
Hello all, Can anyone tell me how to get in help examples in g_select by typing some commands? Anyway, I want to select bound waters between my ligand and protein using -select 'resname SOL within 0.5 ...'. Any idea? Thanks for any suggestion. Yun -- 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 to g_select bound water molecules
Never mind. I got it. On Mon, Apr 2, 2012 at 4:20 PM, Yun Shi yunsh...@gmail.com wrote: Hello all, Can anyone tell me how to get in help examples in g_select by typing some commands? Anyway, I want to select bound waters between my ligand and protein using -select 'resname SOL within 0.5 ...'. Any idea? Thanks for any suggestion. Yun -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Re: how to g_select bound water molecules
OK. I have a protein ab and ligand lig, and I want to select the bound water between them. So I used g_select -f trajout.xtc -s mdrun.tpr -n index.ndx -select 'boundwaters resname SOL and within 0.4 of group ab and within 0.4 of group lig' -on boundwaters.ndx Then I can use the ndx file to do other things. The hard part is what after -select On Mon, Apr 2, 2012 at 5:11 PM, Dallas Warren dallas.war...@monash.eduwrote: For future reference of others that may struggle with g_select, it would be a very good idea for you to put on the emailing list how you achieved it. ** ** Catch ya, Dr. Dallas Warren Medicinal Chemistry and Drug Action Monash Institute of Pharmaceutical Sciences, Monash University 381 Royal Parade, Parkville VIC 3010 dallas.war...@monash.edu +61 3 9903 9304 - When the only tool you own is a hammer, every problem begins to resemble a nail. ** ** *From:* gmx-users-boun...@gromacs.org [mailto: gmx-users-boun...@gromacs.org] *On Behalf Of *Yun Shi *Sent:* Tuesday, 3 April 2012 9:58 AM *To:* Discussion list for GROMACS users *Subject:* [gmx-users] Re: how to g_select bound water molecules ** ** Never mind. I got it. On Mon, Apr 2, 2012 at 4:20 PM, Yun Shi yunsh...@gmail.com wrote: Hello all, Can anyone tell me how to get in help examples in g_select by typing some commands? Anyway, I want to select bound waters between my ligand and protein using -select 'resname SOL within 0.5 ...'. Any idea? Thanks for any suggestion. Yun ** ** -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Scale down/ up velocities when coupling temperature and pressure?
Hi all, I am just wondering about how GROMACS works when coupling temperature and pressure. Assuming the simplest coupling method, like Berendsen, does GROMACS just scale down/ up kinetic energy (KE) for each particle in the system, such as KE1/KE2 = T1/T2 ? Similarly for pressure, KE1/KE2 = P1/P2? So GROMACS scales the square root of velocities for every particle after each time step? - BTW, the number of degrees of freedom = 3N - Nc - Ncom according to the manual. For conventional MD, is Nc equal to the total number of constrained bonds (like all bonds in gromos ff, and all bonds involving hydrogens in amber99sb ff)? Thanks for any explanation. Yun -- 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] sudden jumps in RMSD etc.
Hi all, I am doing duplicate MD simulations with a protein-ligand system. After processing one trajectory by trjconv with the optioin -pbc nojump, I still find abrupt jumps (on the scale of nm) in RMSDs and COM distances. Then I tried -pbc mol -ur compact, which did not work. And then -fit progressive on protein atoms, which again did not completely eliminate those jumps in protein atom RMSDs. I wonder if I have not used the right option? Thanks for any advice, Yun -- 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 is a central structure
Hello everyone, I am using g_cluster with gromos method to do some clustering, and by default, -cl writes the central structure of each cluster obtained. So I wonder what 'central structure' mean? Assuming that I cluster based on RMSD values relative to the starting conformation, and that I then get a cluster with 3 structures (RMSD = 0.1, 0.2, 0.3 nm, respectively), should the structure with RMSD = 0.2 nm be the central structure? What if there are 4 structures in this cluster, with RMSD = 0.1, 0.2, 0.3, 0.4 nm, respectively? Thanks for any explanation, Yun -- 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] nrexcl = 2 or 3?
Hi Mark, I do not quite understand. For example, in amber ff, 1-4 interactions (except the part from dihedral interactions) are calculated according to non-bonded parameters and then scaled by 1/2 or 5/6. When setting nrexcl = 3, which is the default, aren't 1 - 4 interactions excluded from using non-bonded parameters? Thanks, Yun On Wed, Dec 14, 2011 at 4:33 PM, Mark Abraham mark.abra...@anu.edu.auwrote: On 15/12/11, *lq z *lqz@gmail.com wrote: Dear GMXers, I need to use charmm ff (same question for amber ffs), and am confused with what value I should give to nrexcl in [ moleculetype ], 2 or 3? Whatever value pdb2gmx generates for CHARMM, probably 3. According to the manual (v 4.5), nrexcl = 3 stands for excluding non-bonded interactions between atoms that are no further than 3 bonds away. In charmm ff, 1-4 interactions are needed and it is no further than 3 bonds away, ... but 1-4 interactions are regarded as bonded interactions. so I should use nrexcl = 2. However, it is 3 in charmm27.ff/ folder. Yes, CHARMM has special parameters for some 1-4 interactions, which are dealt with seperately. By the way, I did a simple test with a dihedral a-b-c-d. Only a and d have charges. If I use nrexcl =3, no electrostatics energy is reported while the value matches with the one calculated by hand when I set nrexcl to be 2. As expected. Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users 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] MD simulation of Glycoproteion
2011/12/14 陈应广 525342...@qq.com ** Dear gromacs users I used Gromacs in order to get a MD simulation of Glycoproteion.now I have got the Glycoproteion's PDB file,When I want to MD by GMX,it gave a Warnning:Fatal error: Residue not found in residue topology database. And I know it was because of the force field,I want to know which force field should be choose in GMX,And where I can get the force field?? *Any suggestions? * I guess you have to check the residue names, especially those glyco parts, which are probably not in the force field you chose. ** That's all. 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 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 edit a .xvg file?
Hi everyone, I wonder if there is any tool similar to trjconv that can be used to edit a .xvg file? I just want to extract selected data points (like taking data every 10 ps when the original .xvg contains data every 1 ps) and to concatenate .xvg files according to the time of data. Thanks, Yun -- 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] how to edit a .xvg file?
Thanks for the reply. In my case, they are the pullf and pullx files from a COM pulling simulation. So does this mean writing a script would be the best way? Thanks, Yun == == -- Message: 5 Date: Wed, 23 Nov 2011 19:04:33 -0500 From: Justin A. Lemkul jalem...@vt.edu Subject: Re: [gmx-users] how to edit a .xvg file? To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: 4ecd8a11.9030...@vt.edu Content-Type: text/plain; charset=ISO-8859-1; format=flowed Yun Shi wrote: Hi everyone, I wonder if there is any tool similar to trjconv that can be used to edit a .xvg file? I just want to extract selected data points (like taking data every 10 ps when the original .xvg contains data every 1 ps) and to concatenate .xvg files according to the time of data. Several options: 1. Use trjconv to reduce output frequency. 2. Use simple scripting or *nix operations to manipulate the resulting .xvg files. Skipping lines in output is trivial for a language like Perl, and concatenation can be done with simple *nix utilities like 'cat' and others. -Justin -- -- Message: 6 Date: Thu, 24 Nov 2011 13:23:02 +1100 From: Mark Abraham mark.abra...@anu.edu.au Subject: Re: [gmx-users] how to edit a .xvg file? To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: 4ecdaa86.7050...@anu.edu.au Content-Type: text/plain; charset=ISO-8859-1; format=flowed On 24/11/2011 10:31 AM, Yun Shi wrote: Hi everyone, I wonder if there is any tool similar to trjconv that can be used to edit a .xvg file? I just want to extract selected data points (like taking data every 10 ps when the original .xvg contains data every 1 ps) and to concatenate .xvg files according to the time of data. It's usually easiest to re-create the file using the -dt option of the tool. Note that the command that created the .xvg file is recorded in a comment in the file. Mark -- -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Re: how to edit a .xvg file?
Sorry for this question. The bash script turned out to be a one-liner. Yun -- 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] the force constant in constant speed umbrella pulling
Hi everyone, I am just wondering about the mathematical and physical meaning of this force constant when pulling a ligand from its receptor. So I can imagine a dummy atom is linked to the ligand via a spring, and it is moving away from the receptor at 1 nm/ns with the spring force constant 1000 kJ/mol/nm^2. After it moved a very small distance, for example 0.0001 nm, which would be the change in the length of the spring, there will be a corresponding force 0.1 kJ/mol/nm, which is 0.166 pN. So is this force only directly exerted on the COM of the ligand? Or it's been distributed across every atom within the ligand? From this imagined model, is seems the force constant does not have as much effect on the quality (i.e. better quality means closer to a reversible thermodynamic process) of the pulling process as the pulling rate. But I wonder, would reducing force constant by 50% have the same effect on decreasing computational load as reducing pulling rate by 50%? Supposing a set pulled distance 5 nm is required, which one should save more computational time? Thanks for any suggestion. Yun -- 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] picking out resident water molecule?
Hi everyone, I am doing MD simulation with a protein-ligand system, and I want to pick out the water molecules (their residue numbers or coordinates in any frame) that simultaneously contact (within 0.4 nm range for heavy atoms) the ligand and the protein, so that I could plot the lifetime of these resident water molecules within a trajectory. Does anyone know how to achieve this within GROMACS? Thanks for any advice. Yun -- 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: sudden drop of minimal periodic distance
Hi Tsjerk, Thanks for the advice. After I did the nojump conversion of trajectories, there were still quite a few frames with min. dist. under 2.0 nm, but all of them were above 1.4 nm, which I set for rvdw. This is not ideal, but I wonder if it's still OK to use this trajectory? Thanks, Yun Message-ID: cabze1sg24mu7zsh7asvmfainrjsmvgkomg3uy-ixftokwvr...@mail.gmail.com Content-Type: text/plain; charset=iso-8859-1 Hi Yun, Make sure to remove jumps from the trajectory (trjconv -pbc nojump) before using g_mindist. Also visually check a frame that is reported to have closed contacts. Hope it helps, Tsjerk On Nov 10, 2011 1:45 AM, Yun Shi yunsh...@gmail.com wrote: Sorry, I just found that even if I use a dodecahedron box with -d 1.2 nm, the min periodic image dist still dropped abruptly to 0.172 or something like this after around 35 ns or 30 ns (different trajectory with same topology). So I wonder if this is just inevitable and we should live with it? Thanks, Yun On Wed, Nov 9, 2011 at 4:18 PM, Yun Shi yunsh...@gmail.com wrote: Hello everyone, I am ... -- 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 -- next part -- An HTML attachment was scrubbed... URL: http://lists.gromacs.org/pipermail/gmx-users/attachments/2010/7b9ed4a5/attachment-0001.html -- -- 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] where is Coul-LR?
So would it be reasonable to set rcoulomb = 2 or even 3 nm when rerunning a trajectory? I am looking at a ligand-antibody system, and I guess the long-range electrostatic interactions will not be small. A trick proposed by Nicolas in the mailing list during 2007 is to set charges to 0.00 for everything except the atoms we are interested in. This would make Coul-LR: = Coul.-recip., which is calculated by PME (more accurate than just setting rcoulomb to 3nm?) and we can read from .edr file. So I wonder if this method is theoretically sound? Thanks, Yun -- 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] sudden drop of minimal periodic distance
Hello everyone, I am using g_mindist with -pi option to look at the minimal distance between periodic images of my protein-ligand system. It appears that after a certain amount of time (12 ns or 30 ns or ...), there would be a sudden drop of min distance from well above 2 nm to around 0.15 nm. I wonder if this is caused by the octahedral box I used in the MD simulation. But should I still be able to keep the part of trajectory before the sudden drop as valid so that I could do some analysis with my system? Thanks, Yun -- 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: sudden drop of minimal periodic distance
Sorry, I just found that even if I use a dodecahedron box with -d 1.2 nm, the min periodic image dist still dropped abruptly to 0.172 or something like this after around 35 ns or 30 ns (different trajectory with same topology). So I wonder if this is just inevitable and we should live with it? Thanks, Yun On Wed, Nov 9, 2011 at 4:18 PM, Yun Shi yunsh...@gmail.com wrote: Hello everyone, I am using g_mindist with -pi option to look at the minimal distance between periodic images of my protein-ligand system. It appears that after a certain amount of time (12 ns or 30 ns or ...), there would be a sudden drop of min distance from well above 2 nm to around 0.15 nm. I wonder if this is caused by the octahedral box I used in the MD simulation. But should I still be able to keep the part of trajectory before the sudden drop as valid so that I could do some analysis with my system? Thanks, Yun -- 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] where is Coul-LR?
Hello all, I understand that setting rcoulomb rlist should give me Coul-LR from the .edr file. But I set rcoulomb = rlist since PME was used to calculate long range electrostatic interactions, and when I tried g_energy, I only have: 58 Coul-SR:Protein-LIG 59 LJ-SR:Protein-LIG 60 LJ-LR:Protein-LIG 61 Coul-14:Protein-LIG 62 LJ-14:Protein-LIG So where is Coul-LR? How should I get it? Thanks, Yun -- 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] fudge QQ values
Hello everyone, I just have a quick question. So the fudgeQQ value for amber99sb force field in gromacs is 0.8333, but I wonder if we should use 0.83 or 0.8333? I quite different simulation results by manipulating this value. In my antibody - ligand system, 0.8333 showed ligand not stable, dropping out of the binding site, while 0.8333 showed the same ligand fairly stable. The trajectories are about 50 ns long. Thanks for any suggestion. Yun -- 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] Thioester bond problem
Hi Alberto, I used a stupid method to deal with this kind of non-standard moiety. You can use AmberTools to parameterize your thioester, together with your proteins, and then use acpype to convert they topology and coordinate files to gromacs format. But I am not sure if you want to use amber force field. Yun -- 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] position restraints on heavy atoms or all?
Hello all, I am using amber99SB to model an antibody with organic ligands. I know we could choose to restrain all atoms or only heavy atoms during the equilibration. But I wonder if this really matters for my system. As far as I know, equilibration is aimed at getting the temperature and pressure right, mainly for the solvent. And as long as we are going to do a 'real' MD production run without restraints, it should not matter too much how we restrained the antibody and ligands? I am not if I am right about this. Any suggestions? Thanks, Yun -- 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] would trjcat correctly delete extra frames?
Hi all, I am doing simulations on cluster piece by piece with -maxh and -noappend options of mdrun. However, one piece crushed way before approaching the max hours for unknown reasons. As a result, the part0004.trr file contains a couple of frames ahead of part0004.cpt file, since the most .cpt was only written about 5 or 10 minutes before the crush. So I plan to continue with the part0004.cpt file, and I expect the part0005.trr file I will get would have its initial 3 or 4 frames duplicated (only # steps and time) with the last 3 or 4 frames in part0004.trr file. And when I concatenate them together with trjcat, would the last 3 or 4 frames in part0004.trr file be automatically deleted? and only the initial 3 or 4 frames in part0005.trr be included in the resulting .trr file? Thanks for any suggestion! Yun -- 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: reasons for slow computation?
Hi Mark, I am not quite sure. While the details are different, I guess it is still valid to compare some statistical properties, right? For example, I should still be able to compare the average COM distance of these two ligands over long trajectories, and make the conclusion that one is better than the other since the distance is smaller. So should I be able to compare any averaged, or say, clustered properties? Thanks, Yun On 11/10/2011 1:40 PM, Yun Shi wrote: Hi Justin, I guess you are right, that some processors on that cluster appear to be much slower than others. More likely is that the mapping of processes to processors is faulty, as Justin said. But I am still wondering that, would the difference in initial maximum inter charge-group distances (0.451 nm vs 0.450 nm) and minimum initial size of DD gird (0.620nm vs 0.618nm) make the two simulations NOT comparable? Because I want to compare the interactions of two different ligands with the same protein. That's just a detail of how the simulation is being implemented in parallel. There are a vast number of algorithmically equivalent ways to do this, and the initial conditions determine which is chosen. Since your initial conditions differ, so do the details. Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] reasons for slow computation?
Hi all, I am doing MD simulation on two almost identical protein-ligand systems with GROMACS4.5.4 and amber99SB force fields. Almost every single parameter I used for this two systems are the same (I literally copied the mdp files for both), except that one has 65235 atoms while the other has 65205 atoms, which is caused by the difference of 10 water molecules. However, the 65235-atom one is significantly slower than the 65205-atom one, the computational speed (performance) ratio is roughly 1:3. I am really confused, since the only differences between them that I can think of are: 1. The configuration of ligands, that is, alpha sugar vs beta sugar 2. number of water molecules, 19477 vs 19467 3. octahedral box size, ~726.119 nm^3 vs ~726.322 nm^3 4. probably the initial velocities, since I am used gen_seed= -1 5. the cpu specs, as I am using a cluster (I specified procs=72 for both system and I assume they were assigned to different cpus randomly), which consists of Intel Xeon E5450 quad-core processor, running at 3.0 GHz, and Intel Xeon X5650 6-core processors running at 2.67 GHz How come the performance could have a 1:3 difference? Thanks for any suggestion! Yun -- 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: reasons for slow computation?
And another difference I noticed from .log files are: .. Initial maximum inter charge-group distances: two-body bonded interactions: 0.451 nm, LJ-14, atoms 3586 4808 multi-body bonded interactions: 0.451 nm, Proper Dih., atoms 3586 4808 Minimum cell size due to bonded interactions: 0.496 nm Optimizing the DD grid for 48 cells with a minimum initial size of 0.620 nm The maximum allowed number of cells is: X 12 Y 12 Z 12 Domain decomposition grid 8 x 3 x 2, separate PME nodes 24 PME domain decomposition: 8 x 3 x 1 Interleaving PP and PME nodes This is a particle-particle only node Domain decomposition nodeid 0, coordinates 0 0 0 Using two step summing over 21 groups of on average 2.3 processes ... for the slow one, while for the fast one: Initial maximum inter charge-group distances: two-body bonded interactions: 0.450 nm, LJ-14, atoms 2064 3014 multi-body bonded interactions: 0.450 nm, Proper Dih., atoms 2064 3014 Minimum cell size due to bonded interactions: 0.494 nm Optimizing the DD grid for 48 cells with a minimum initial size of 0.618 nm The maximum allowed number of cells is: X 12 Y 12 Z 12 Domain decomposition grid 8 x 3 x 2, separate PME nodes 24 PME domain decomposition: 8 x 3 x 1 Interleaving PP and PME nodes This is a particle-particle only node Domain decomposition nodeid 0, coordinates 0 0 0 Using two step summing over 18 groups of on average 2.7 processes Are these differences likely to generate huge performance difference? Thanks, Yun -- 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] not guaranteed to be binary identical ?
Hi all, I am doing my MD simulations piece by piece, with -maxh and -noappend options, so that I can link pieces of trajectories together afterward. But as I did part0001 with 48 cores, and part0002 with 72 cores, the log file told me that: #nodes mismatch, current program: 72 checkpoint file: 48 #PME-nodes mismatch, current program: -1 checkpoint file: 16 Gromacs binary or parallel settings not identical to previous run. Continuation is exact, but is not guaranteed to be binary identical. I do not quite understand what it really means by binary identical? Would this affect my results? In other words, would using different # cores for each part give me different conformational assembles from using the same # cores each part? Thanks for any help! Yun -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Re: reasons for slow computation?
Hi Justin, I guess you are right, that some processors on that cluster appear to be much slower than others. But I am still wondering that, would the difference in initial maximum inter charge-group distances (0.451 nm vs 0.450 nm) and minimum initial size of DD gird (0.620nm vs 0.618nm) make the two simulations NOT comparable? Because I want to compare the interactions of two different ligands with the same protein. Thanks, Yun -- 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] acpype generated different tip3p water paramters
Hi Alan, So is acpype using a conversion factor of 4.184 for dihedral force constant? I found some dihedral constants as 0.156 in the amber format, which should be 0.156*4.184=0.652704 in gromacs unit. However, acpype gave a force constant of 0.65084 after conversion, which is slightly off. I wonder if this is OK, and I suspect it might be that I used rdparm to check the amber format value, which only gives 3 decimals for force constants. Thanks for the reply, Yun Hi Yun, ACPYPE is working fine. What happens here is I choose the reproduce the exact values one sees in AMBER. Now why GMX tip3p file choose a different value, I don't know. Nevertheless, it's pretty simple to put whatever value you want there if you think you need. Alan -- 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] different tip3p water parameters between AMBER and GROMACS
Hi all, I found that the 99SB force field in AMBERTOOLS1.5 and GROMACS4.5.4 have different force constants for bond and angle parameters, after converting them to the same unit. So for those in AMBER [ moleculetype ] ; molname nrexcl ; TIP3P model WAT 2 [ atoms ] ; nr type resnr residue atom cgnr charge mass 1 OW 1 WAT O 1 -0.834 16.0 2 HW 1 WATH1 1 0.4171.00800 3 HW 1 WATH2 1 0.4171.00800 #ifdef FLEXIBLE [ bonds ] ; i j funct length force.c. 1 2 1 0.09572 462750.4 0.09572 462750.4 1 3 1 0.09572 462750.4 0.09572 462750.4 [ angles ] ; i j k funct angle force.c. 2 1 3 1 104.520836.800 104.520836.800 #else [ settles ] ; i j funct length 1 1 0.09572 0.15139 [ exclusions ] 1 2 3 2 1 3 3 1 2 #endif while for GROMACS, that is, in the amber99sb.ff/tip3p.itp file, [ moleculetype ] ; molname nrexcl SOL 2 [ atoms ] ; id at type res nr res name at name cg nr chargemass 1 OW 1 SOL OW 1 -0.83416.0 2 HW 1 SOL HW1 1 0.417 1.00800 3 HW 1 SOL HW2 1 0.417 1.00800 #ifndef FLEXIBLE [ settles ] ; OWfunct doh dhh 1 1 0.09572 0.15139 [ exclusions ] 1 2 3 2 1 3 3 1 2 #else [ bonds ] ; i j funct length force_constant 1 2 1 0.09572 502416.0 0.09572502416.0 1 3 1 0.09572 502416.0 0.09572502416.0 [ angles ] ; i j k funct angle force_constant 2 1 3 1 104.52 628.02 104.52 628.02 #endif So I wonder why different parameters for tip3p water are used in GMX, even for the save force field ? Thanks for any suggestion, Yun -- 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] acpype generated different tip3p water paramters
Hi all, I just noted that the tip3p water converted from amber format to gromacs format is [ moleculetype ] ; molname nrexcl ; TIP3P model WAT 2 [ atoms ] ; nr type resnr residue atom cgnr charge mass 1 OW 1 WAT O 1 -0.834 16.0 2 HW 1 WATH1 1 0.4171.00800 3 HW 1 WATH2 1 0.4171.00800 #ifdef FLEXIBLE [ bonds ] ; i j funct length force.c. 1 2 1 0.09572 462750.4 0.09572 462750.4 1 3 1 0.09572 462750.4 0.09572 462750.4 [ angles ] ; i j k funct angle force.c. 2 1 3 1 104.520836.800 104.520836.800 #else [ settles ] ; i j funct length 1 1 0.09572 0.15139 [ exclusions ] 1 2 3 2 1 3 3 1 2 #endif while in amber99sb.ff/tip3p.itp. it is moleculetype ] ; molname nrexcl SOL 2 [ atoms ] ; id at type res nr res name at name cg nr chargemass 1 OW 1 SOL OW 1 -0.83416.0 2 HW 1 SOL HW1 1 0.417 1.00800 3 HW 1 SOL HW2 1 0.417 1.00800 #ifndef FLEXIBLE [ settles ] ; OWfunct doh dhh 1 1 0.09572 0.15139 [ exclusions ] 1 2 3 2 1 3 3 1 2 #else [ bonds ] ; i j funct length force_constant 1 2 1 0.09572 502416.0 0.09572502416.0 1 3 1 0.09572 502416.0 0.09572502416.0 [ angles ] ; i j k funct angle force_constant 2 1 3 1 104.52 628.02 104.52 628.02 #endif So it seems that the force_constants for O-H bond and H-O-H angle are different? Does this mean amber and gromacs use different parameters for tip3p water? or it's just acpype is not working right? Thanks, Yun -- 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] different sets of fudgeQQ and fudgeLJ
Hi Mark, I don't quite understand what it follows only one [defaults] section can exist in an entire topology. Then how to specify different fudge values for different subsets of non-bonded interactions? When defining new atom types, should I always use 'new' atomtypes names? For example, if the atom type H2 already exists for part A, then I should use something different like H2B to define similar atomtypes in part B? But if I can use only one [ defaults ] section, even within different .itp files under the same .top file, how could I tell the program to apply different fudge values to H2 as defined in A.itp and H2B as defined in B.itp? Thanks, Yun Date: Sun, 18 Sep 2011 09:25:38 +1000 From: Mark Abraham mark.abra...@anu.edu.au Subject: Re: [gmx-users] different sets of fudgeQQ and fudgeLJ To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: 4e752c72.4040...@anu.edu.au Content-Type: text/plain; charset=ISO-8859-1; format=flowed On 18/09/2011 8:58 AM, Yun Shi wrote: Hi all, I want to apply different values of LJ and QQ scaling factors for two interacting molecules A and B. Since I already have the .itp files for each molecule, should I just add something like: [ defaults ] ; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ 12 yes 0.5 0.8333 at the very beginning of A.itp and B.itp respectively (the actual values are different from this eg.)? I wonder if .itp file format was designed for this kind of purpose? :) I guess we could even define different sets of atom types for A and B, right? As you will see in the examples in chapter 5, only one [defaults] section can exist in an entire topology. Further, the fudge values only apply to a subset of non-bonded interactions. For the kind of calculation you have in mind, you will need to define new atom types for the modified interactions. Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Re: different sets of fudgeQQ and fudgeLJ
Hi Mark and Justin, I think I should be more specific here. So i.e., I want to study the interaction between a protein receptor and a carbohydrate ligand with MD simulation, and I plan to use ff99sb for protein while glycam06 for carbohydrate. Since the two force fields are parameterized using different scaling factors, I should use different fudge values, although these values will not DIRECTLY affect inter-molecular interactions. So how to set this up, if it is possible, within .top file and .itp files? Thanks, Yun -- 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] different sets of fudgeQQ and fudgeLJ
Hi all, I want to apply different values of LJ and QQ scaling factors for two interacting molecules A and B. Since I already have the .itp files for each molecule, should I just add something like: [ defaults ] ; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ 12 yes 0.5 0.8333 at the very beginning of A.itp and B.itp respectively (the actual values are different from this eg.)? I wonder if .itp file format was designed for this kind of purpose? :) I guess we could even define different sets of atom types for A and B, right? Thanks, Yun -- 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: amb2gmx.pl to convert GLYCAM topology
Hi Alan, For example, in the Glycam_06g.dat file, you can find: OH-CG-CG-OS 1 -1.10 0.0-1 So this dihedral parameter has a force constant of -1.10, and this is what I mean by GLYCAM force field assigns negative force constants to some dihedrals. I did try the GMX45 approach, and using the conversion factor 4.184, I got a difference of about 0.05%. I am not sure if this is caused not setting step size in the sander minimization. Regards, Yun Date: Tue, 13 Sep 2011 12:03:10 +0100 From: Alan alanwil...@gmail.com Subject: Re: [gmx-users] Re: amb2gmx.pl to convert GLYCAM topology To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: caeznbznq98pbhjdr1smpogdtz_pgfqsk40izlfpq7-nizyw...@mail.gmail.com Content-Type: text/plain; charset=utf-8 Hi Yun, Have you read http://ambermd.org/formats.html? In particular, this note: NOTE: *the atom numbers in the following arrays that describe bonds, angles, and dihedrals are coordinate array indexes for runtime speed. The true atom number equals the absolute value of the number divided by three, plus one. In the case of the dihedrals, if the fourth atom is negative, this implies that the dihedral is an improper. If the third atom is negative, this implies that the end group interations are to be ignored. End group interactions are ignored, for example, in dihedrals of various ring systems (to prevent double counting of 1-4 interactions) and in multiterm dihedrals. * I may be failing to understand what you mean by GLYCAM force field assigns negative force constants to some dihedrals. Anyway, since GMX 4.5 can go without RB convertions, you can do this: acpype -x disac.inpcrd -p disac.prmtop --gmx45 If you have sander, you can do just one step of EM and compare against one step EM with GMX. Do the proper conversions and Energies diff should be 0.001%. Cheers, Alan On 12 September 2011 21:21, Yun Shi yunsh...@gmail.com wrote: Hi all, I am not a CS person, but I did find something in acpype.py as . if phase in [0, 180]: properDihedralsGmx45.append([ item[0].atoms, phaseRaw, kPhi, period]) if not self.gmx45: if kPhi 0: V[period] = 2 * kPhi * cal if period == 1: C[0] += 0.5 * V[period] if phase == 0: C[1] -= 0.5 * V[period] else: C[1] += 0.5 * V[period] elif period == 2: .. kPhi here seems to be the dihedral force constant, and it seems if kPhi 0, no value will be assigned to C[0], C[1], C[2] ... I wonder if the negative dihedral force constants problem could be solved by changing 'kPhi 0' to 'kPhi != 0' for acpype? Thanks, Yun -- 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 -- Alan Wilter SOUSA da SILVA, DSc Bioinformatician, UniProt - PANDA, EMBL-EBI CB10 1SD, Hinxton, Cambridge, UK +44 1223 49 4588 -- next part -- An HTML attachment was scrubbed... URL: http://lists.gromacs.org/pipermail/gmx-users/attachments/20110913/f9bc31d1/attachment-0001.html -- Message: 5 Date: Tue, 13 Sep 2011 16:46:44 +0530 From: om prakash ombioi...@gmail.com Subject: [gmx-users] Unsubscribe me Please To: gmx-users@gromacs.org Message-ID: CAM5rxXNDPYi_aEEd+7efAmZoPex8B4Um5FwgJKg6hfm=t9y...@mail.gmail.com Content-Type: text/plain; charset=iso-8859-1 Unsubscribe me Please -- Om Prakash Sharma Ph.D Scholar DIT JRF Centre for Bioinformatics Pondicherry University Pondicherry-605014 -- next part -- -- 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: amb2gmx.pl to convert GLYCAM topology
Hi all, I am not a CS person, but I did find something in acpype.py as . if phase in [0, 180]: properDihedralsGmx45.append([item[0].atoms, phaseRaw, kPhi, period]) if not self.gmx45: if kPhi 0: V[period] = 2 * kPhi * cal if period == 1: C[0] += 0.5 * V[period] if phase == 0: C[1] -= 0.5 * V[period] else: C[1] += 0.5 * V[period] elif period == 2: .. kPhi here seems to be the dihedral force constant, and it seems if kPhi 0, no value will be assigned to C[0], C[1], C[2] ... I wonder if the negative dihedral force constants problem could be solved by changing 'kPhi 0' to 'kPhi != 0' for acpype? Thanks, Yun -- 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] amb2gmx.pl to convert GLYCAM topology
Hi all,
I understand this problem has been discussed before, but it seems no
conclusion has been drawn.
GLYCAM force field assigns negative force constants to some dihedrals, and
when amb2gmx.pl was used to convert prmtop file to gromacs top file, these
negative values seem to be ignored. Some people proposed that we change the
code in amb2gmx.pl, that is:
...
# get all force constants for each line of a dihedral #
my $lines = $i -1 +$numijkl;
for(my $j=$i;$j=$lines;$j++){
my $period = abs($pn{$j});
if($pk{$j}0) {
$V[$period] = 2*$pk{$j}*$cal/$idivf{$j};
}
...
the $pk{$j}0 is modified to $pk{$j}!=0.
Others suggest to modify the original prmtop file, that is, to remove the
negative signs, and correspondingly, change the phase shift from 0 to 180.
Then amb2gmx.pl could be used to correctly convert the topology.
I am wondering if the first approach has been validated, since the second
one seems complicated and laborious to carry out.
Thanks for any advice,
Yun
--
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: amb2gmx.pl to convert GLYCAM topology
Hi Alan,
I am not sure if my acpype version is not updated. But I did try, and it
behaved the same as amb2gmx.pl for dihedrals.
Yun
Or why not trying acpype?
Cheers,
Alan
On 9 September 2011 07:37, Mark Abraham mark.abra...@anu.edu.au wrote:
On 9/09/2011 4:21 PM, Yun Shi wrote:
Hi all,
I understand this problem has been discussed before, but it seems no
conclusion has been drawn.
Someone needs to do some work and report back :-)
GLYCAM force field assigns negative force constants to some dihedrals, and
when amb2gmx.pl was used to convert prmtop file to gromacs top file, these
negative values seem to be ignored. Some people proposed that we change
the
code in amb2gmx.pl, that is:
...
# get all force constants for each line of a dihedral #
my $lines = $i -1 +$numijkl;
for(my $j=$i;$j=$lines;$j++){
my $period = abs($pn{$j});
if($pk{$j}0) {
$V[$period] = 2*$pk{$j}*$cal/$idivf{$j};
}
...
the $pk{$j}0 is modified to $pk{$j}!=0.
Others suggest to modify the original prmtop file, that is, to remove the
negative signs, and correspondingly, change the phase shift from 0 to 180.
Then amb2gmx.pl could be used to correctly convert the topology.
I am wondering if the first approach has been validated, since the second
one seems complicated and laborious to carry out.
Seems like a straightforward job for regular expression replacement using
sed/perl/python/whatever. It might even be a one-liner.
Mark
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
Alan Wilter SOUSA da SILVA, DSc
Bioinformatician, UniProt - PANDA, EMBL-EBI
CB10 1SD, Hinxton, Cambridge, UK
+44 1223 49 4588
-- next part --
An HTML attachment was scrubbed...
URL:
http://lists.gromacs.org/pipermail/gmx-users/attachments/20110909/cce09225/attachment.html
--
--
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] missing parameters in auto-generated topology file?
Hi all, I used pdb2gmx and selected amber99sb for generation of itp files of a normal peptide within GROMACS 4.5.4. But I saw that all the bonds, angles, and dihedral parameters (c0, c1, c2 ...) were not present in the itp file, while only funct is defined. It seems the same thing happens with opls ff. So I wonder if those parameters are automatically generated according to corresponding atom types when implementing grompp? What would happen if I have exotic species and some atomic sequences are not defined in the [ angletypes ] or [ dihedraltypes ] section of amber99sb.ff/ffboned.itp ? In the other case, when I used amb2gmx.pl to convert glycam force field-parametized prmtop file to gromacs top file, I saw something like: [ dihedrals ] ;i j k l func C0 ... C5 19 18 20 21 3 0.75312 2.25936 0.0 -3.01248 0.0 0.0 ; 41 118 19 3 0.20920 0.62760 0.0 -0.83680 0.0 0.0 ; Although I understand '3' in the fifth column indicates a Ryckaert-Bellemans-potential function, but what parameters are those numbers followed (0.75312 2.25936 0.0-3.01248 0.0 0.0) correspond to? Thanks for any help! Yun -- 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] gd_29 or gd_41 ?
Hi, For a H-NL-CH1-CH2 (H1-N-CA-C) dihedral angle in a N-terminal MET, why would pdb2gmx automatically assign gd_29 ? In /gromos53a6.ff/ffbonded.itp, it appears: #define gd_29 0.000 3.77 3 ; -C,CHn,SI- 0.9 ; ... #define gd_41 0.000 3.77 6 ; -CHn-NT- 0.9 ; ... Thanks, Yun Shi -- 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: gd_29 or gd_41 ?
similarly, why would CA-CB-CG-OD2/1 in ASP be gd_40 instead of gd_34 ? #define gd_34 0.000 5.92 3 ; -CHn,SI-CHn-1.4 ; ... #define gd_40 0.0001.0 6 ; -CHn-C,NR(ring), CR1- 0.24 ; Something wrong with my interpretation of the comments after each definition? Thanks, Yun Shi On Wed, Aug 24, 2011 at 3:00 PM, Yun Shi yunsh...@gmail.com wrote: Hi, For a H-NL-CH1-CH2 (H1-N-CA-C) dihedral angle in a N-terminal MET, why would pdb2gmx automatically assign gd_29 ? In /gromos53a6.ff/ffbonded.itp, it appears: #define gd_29 0.000 3.77 3 ; -C,CHn,SI- 0.9 ; ... #define gd_41 0.000 3.77 6 ; -CHn-NT- 0.9 ; ... Thanks, Yun Shi -- 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] sugar dihedral parameters
Hi, In the ATB-distributed 53a6_carbo_new.dat, there are more than one set of parameters specified for some dihedral angles. For example, O5-C5-C6-O6 has a 1-fold term as gd_5, and a 3-fold term gd_37; C3-C2-C1-O5 has a 2-fold term as gd_17, and a 3-fold term as gd_34. So which one should I use? Thanks, Yun Shi -- 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] parameters for CH1 -- S -- CH2 ?
Hi all, I have a small molecule with a part like: CH1 --- CH1 --- S --- CH2 --- CH1 | CH1 I first tried PRODRG, and it turned out PRODRG assigned bond, angle, and dihedral parameters according what are already present in the gromos53a6.ff/ffbonded.itp file. Since this file does not contain information regarding this CH1 --- S --- CH2 pattern, PRODRG just assigned bond (gb_31),angle (ga_4) as if my molecule has CH3 --- S --- CH2 pattern, and dihedral parameters for CH1 -- CH1 -- S -- CH2 as if it was CH1 -- CH2 -- S -- CH2. So I guess I have to do this manually. But as no corresponding parameters defined in ffbonded.itp file, should I derive these parameters via chemical intuition and try validate them in some way? Or is there an easier alternative way to do this? Thanks for any suggestion! Yun -- 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] validation of ligand parameters
Hi all, I have thought of three experimental values to validate the parameters I assigned for a oligosaccharide, namely, free enthalpy of solvation, nOe effect, and some J3-couplings. For free enthalpy of solvation, I wondered why delta G was used to indicate enthalpy? Would what be measured is just the heat, i.e., delta H ? For nOe effect, should I use g_rmsdist to calculate 1/r3 and 1/r6 averaged distances? Should I correlate these values to corresponding nOe initial build-up rate or just a normal steady-state NOE? For J-couplings, I cannot find a tool within GROMACS to do this. The command g_chi only computes NMR 3J coupling constants for amino acid backbond and sidechain atoms? Thanks for any suggestion! Yun Shi -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] To modify PRODRG-generated .itp file?
Hi all, In addition to modifying atomic charges to the .itp file generated by PRODRG, I also want to change the vdw C12 and C6 parameters of selected 1-4 LJ interactions and modify some 1-5 LJ interactions (or any non-bonded LJ interactions) with the vdw C12 and C6 parameters I prefer. So should I just make changes to the [ pairs ] section of the .itp file? For example, change ; ai aj fuc0, c1, ... 1 4 1 to ; ai aj fuc0, c1, ... 1 4 1 0.004724 1.84e-6 How can I make sure it will override the parameters in the [ pairtypes ] section of the ffnonbonded.itp file? And I should make these values (c0 for C6, c1 for C12) in accordance with the units in [ pairtypes ] section, right? Moreover, I want to change some dihedral parameters. But I saw the .itp file generated by PRODRG as something like [ dihedrals ] ; ai aj ak al fuc0, c1, m, ... 10 5 4 3 1 0.05.9 3 0.05.9 3 so 0.0 is the phase shift, 5.9 is the force constant, 3 is the multiplicity, but why these values appeared twice in a row? When I modify them, should I also write them twice as this format? Will it override dihedral parameters that are already in the ffbonded.itp file? Thanks very much, Yun -- 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] LINCS with amber99SB ?
Hi all, I want to do MD simulation with amber99SB force field, which I found originally developed using SHAKE algorithm to constrain all bonds involving hydrogen atoms. But it seems to me that SHAKE is still not supported with domain decomposition in GROMACS4.5.4, and we can only use LINCS for bond constraints if we want to run MD simulations in parallel. So what should I use in simulating with amber99SB in parallel? LINCS constraining allbonds? LINCS constraining hbonds only? Or just run with SHAKE in one node? I would be very grateful if someone could tell me the relative accuracy (or if they are all good enough) and computational costs of these options above. And maybe I should ask this question in the developers mailing list, but would GROMACS support SHAKE in parallel in the near future, like in a 4.5.5 version? Thanks a lot, Yun Shi -- 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: LINCS with amber99SB ?
Thank you for the reply. Just to clarify. In the original gromos 53a6 paper, all bonds were constrained, while in amber99sb paper only bonds involving H atoms were constrained. So as long as I stick to the same groups of bonds constrained as in the development of respective forcefields, it is always OK to use LINCS rather than SHAKE, right? And for those bonds not constrained when using amber99sb force field, would GROMACS automatically apply the harmonic bond stretching functional form as specified in the [ bonds ] section of corresponding .itp files? Regards, Yun Yun Shi wrote: Hi all, I want to do MD simulation with amber99SB force field, which I found originally developed using SHAKE algorithm to constrain all bonds involving hydrogen atoms. But it seems to me that SHAKE is still not supported with domain decomposition in GROMACS4.5.4, and we can only use LINCS for bond constraints if we want to run MD simulations in parallel. So what should I use in simulating with amber99SB in parallel? LINCS constraining allbonds? LINCS constraining hbonds only? Or just run with SHAKE in one node? You should be able to use LINCS with either all-bonds or hbonds. Running your simulation in serial, even for a relatively small system, will likely be painfully slow. I would be very grateful if someone could tell me the relative accuracy (or if they are all good enough) and computational costs of these options above. Either should be fine. LINCS is actually somewhat more stable than SHAKE. And maybe I should ask this question in the developers mailing list, but would GROMACS support SHAKE in parallel in the near future, like in a 4.5.5 version? There will not be any new features in 4.5.5, and the 4.6 list of new features is already basically agreed upon. Gromacs is undergoing major code changes that limit the amount of new features that will come out soon. You can file a feature request on redmine.gromacs.org, but do not expect to see any activity on the issue until probably version 5.0. -Justin -- == == Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing 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: Re: [gmx-users] charge assignment for the anomeric carbon in 53a6 ff
Hi Justin, Thanks a lot for the replies. I wonder what are the newer versions you indicated, but I find one as: *A reoptimized GROMOS force field for hexopyranose-based carbohydrates accounting for the relative free energies of ring conformers, anomers, epimers, hydroxymethyl rotamers, and glycosidic linkage conformers*. This new version 56ACARBO is claimed to be nearly equivalent to 53A6 for non-carbohydrate systems, and I find, at least, atomic charges for C1 and C5 look more chemically reasonable. But changes for a couple of 1-4 van der Waals parameter C12 were made to 53a6 force field, and additional atom types were added. So could you tell me which file contain those LJ parameters with respect to different atom types? Alternatively, would it be reasonably good enough to use your Lysozyme-JZ4 strategy, i.e., to build initial topology of carbohydrate molecules by PRODRG, and then modify the atomic charges according to the 56ACARBO paper? Thanks, Yun Shi Yun Shi wrote: Hi all, I am doing MD simulation of some carbohydrate-protein complex with this 53a6 force-field. I noted that in any oligosaccharide, the charge assigned for anomeric carbon is 0.232 while C5 (i.e. in a non-terminal mannopyranose) is 0.376. This is kind of counter-intuitive for me since the anomeric carbon C1 is linked to two oxygen atoms O5 and O1 while C5 is linked to only one oxygen atom O5, and both C1 and C5 are CH type atom. And in other force filed like GLYCAM06, C1 always has a larger positive partial charge than C5. I understand that GROMOS force fields are derived in quite different ways, but I still wonder if it makes sense to switch the charges on C1 and C5? No, but the charges used may not be the best ones anyway. There are newer revisions of 53A6 that are much better. It seems these carbohydrate parameters did not change in the new 54a7 ff, which was tested against nOe and J-coupling data rather than free enthalpy of hydration and solvation, which really puzzles me. 54A7 was principally a revision of protein atom types and one lipid type; nothing in the sugar parameters was modified. -Justin -- == == Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] charge assignment for the anomeric carbon in 53a6 ff
Thanks again. Since I am running MD simulation on a public cluster, I cannot modify atomtypes.atp and ffnonbonded.itp on that cluster. But I could get every parameter ready and prepare the xxx.tpr files locally, and pass them to the cluster. So the practical question is, does the tpr file contain all those C12 and C6 LJ parameters already? Or mdrun need to retrieve these parameters according to the atom types from corresponding files within the gromos53a6.ff folder? Regards, Yun Shi Yun Shi wrote: Hi Justin, Thanks a lot for the replies. I wonder what are the newer versions you indicated, but I find one as: /A reoptimized GROMOS force field for hexopyranose-based carbohydrates accounting for the relative free energies of ring conformers, anomers, epimers, hydroxymethyl rotamers, and glycosidic linkage conformers/ . This new version 56A_CARBO is claimed to be nearly equivalent to 53A6 for non-carbohydrate systems, and I find, at least, atomic charges for C1 and C5 look more chemically reasonable. But changes for a couple of 1-4 van der Waals parameter C12 were made to 53a6 force field, and additional atom types were added. So could you tell me which file contain those LJ parameters with respect to different atom types? New atom types need to be added to atomtypes.atp and their parameters defined in ffnonbonded.itp. Alternatively, would it be reasonably good enough to use your Lysozyme-JZ4 strategy, i.e., to build initial topology of carbohydrate molecules by PRODRG, and then modify the atomic charges according to the 56A_CARBO paper? Maybe, but check the bonded parameters carefully. Sugars are particularly sensitive to incorrect angle and dihedral parameters, especially. -Justin -- == == Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] charge assignment for the anomeric carbon in 53a6 ff
Hi all, I am doing MD simulation of some carbohydrate-protein complex with this 53a6 force-field. I noted that in any oligosaccharide, the charge assigned for anomeric carbon is 0.232 while C5 (i.e. in a non-terminal mannopyranose) is 0.376. This is kind of counter-intuitive for me since the anomeric carbon C1 is linked to two oxygen atoms O5 and O1 while C5 is linked to only one oxygen atom O5, and both C1 and C5 are CH type atom. And in other force filed like GLYCAM06, C1 always has a larger positive partial charge than C5. I understand that GROMOS force fields are derived in quite different ways, but I still wonder if it makes sense to switch the charges on C1 and C5? It seems these carbohydrate parameters did not change in the new 54a7 ff, which was tested against nOe and J-coupling data rather than free enthalpy of hydration and solvation, which really puzzles me. Thank you for any explanation and suggestion. Yun Shi -- 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] setting vdw cutoff with specific force-filed?
Hi Justin and Mark, Thank you very much for the reply. I was using table 7 (Normal van der Waals Parameters) to calculate non-bonded vdw interactions that are not between third neighbors, such as CH1 carbons between different chains in a biomolecular system. Anything wrong here? I understand that it is the force that dictates the MD evolution, and I calculated in this case as F = 12 * 9.85^2 / 1.5^13 = 5.98 kJ/mol/nm for the repulsion term. The force from different directions on a atom in a homogeneous system would cancel each other to some extend, but what about the energy arises from this interaction? Would this considerably affect the calculation of, say, binding energy of a ligand to a receptor from thermodynamic integration or pulling simulation? Besides, the GROMOS 53a6 paper used triple range scheme for calculations of nonbonded interactions, and I guess it was rlist = 0.8 nm while rvdw = rcoulomb = 1.4 nm. So is this considered to be accurate enough in calculating free enthalpies of solvation since we know the interactions between 0.8 and 1.4 nm were calculated every 5 steps? The paper also used reaction-field instead PME to account for long-range electrostatic interactions. I heard some people argue that PME would be more accurate and it seemed to be utilized more often even in gromacs tutorials. So does this mean certain accuracy could be achieved by using triple range scheme and reaction-field together because the errors they incur respectively somehow cancel out each other? Thanks a lot, Yun Shi On 04/08/11, Justin A. Lemkul jalem...@vt.edu wrote: Yun Shi wrote: Hi all, I am working with GROMOS 53a6 ff in GROMACS 4.5, and I assume a Lennard-Jones interaction function was used for short-range vdw interactions. From the reference paper /A Biomolecular Force Field Based on the Free Enthalpy of Hydration and Solvation: The GROMOS Force-Field Parameter Sets 53A5 and 53A6/, I found that for example, when rvdw = 1.5nm, the repulsion term of the interaction between two CH1 type atoms (C12ij = 9.85^2) can be calculated as 9.850*9.850 / (1.5^12) = 0.747786 kJ/mol. So I wonder if this value is considered to be small enough to be ignored. You should pay attention to the column headings in table 7 so that you can compute the contribution correctly. However, the magnitude of the energy of any particular interaction is not really of any concern. The evolution of the system depends on the *forces*, and it is likely that the sum of the forces on any atom from all its repulsion interactions from atoms that are (say) 1.4nm to 1.5nm away is very close to zero, except in highly non-homogeneous spatial distributions of particles. In any case, the sum of that contribution will be much smaller than the other contributions. Mark In addition, it seems not until 5 nm does the dispersion term become larger than the repulsion term in this case, so would turning on Dispersion Correction between, say 1.5 to 5 nm introduce more errors than turning it off? You should use the cutoff described the authors of the force field, in this case rvdw=1.4. Unless you can demonstrate that by using a different value you can achieve superior results, stick with the specifics of parameterization. I have never seen ill effects of setting rvdw=1.4 and using dispersion correction with this force field. -Justin -- == == Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Re: setting vdw cutoff with specific force-filed?
Hi Justin, I got it now. During the 10fs, even water molecules with a speed of 500 m/s only have a 0.005 nm displacement, which is far less than than 0.9 nm or 1.4 nm. Thanks again! Yun Yun Shi wrote: Hi Justin and Mark, Thank you very much for the reply. I was using table 7 (Normal van der Waals Parameters) to calculate non-bonded vdw interactions that are not between third neighbors, such as CH1 carbons between different chains in a biomolecular system. Anything wrong here? I understand that it is the force that dictates the MD evolution, and I calculated in this case as F = 12 * 9.85^2 / 1.5^13 = 5.98 kJ/mol/nm for the repulsion term. The force from different directions on a atom in a homogeneous system would cancel each other to some extend, but what about the energy arises from this interaction? Would this considerably affect the calculation of, say, binding energy of a ligand to a receptor from thermodynamic integration or pulling simulation? Check your units and the column headings of Table 7. Plugging in 9.85 as the energy will give you a wildly inflated result. The C12 parameters listed are actually square roots and listed as 10^-3. I think you will find the resulting energies and forces are vastly smaller for a simple interaction between two atoms. Besides, the GROMOS 53a6 paper used triple range scheme for calculations of nonbonded interactions, and I guess it was rlist = 0.8 nm while rvdw = rcoulomb = 1.4 nm. So is this considered to be accurate enough in calculating free enthalpies of solvation since we know the interactions between 0.8 and 1.4 nm were calculated every 5 steps? There is no need to update the neighbor list every single step. Typically, water is the fastest-diffusing molecule in the system, but it will generally not have a dramatic displacement on the scale of 10 fs or so. The paper also used reaction-field instead PME to account for long-range electrostatic interactions. I heard some people argue that PME would be more accurate and it seemed to be utilized more often even in gromacs tutorials. So does this mean certain accuracy could be achieved by using triple range scheme and reaction-field together because the errors they incur respectively somehow cancel out each other? PME is substantially more accurate. Using it also requires rlist=rcoulomb, so the exact details of the Gromos96 derivation may be somewhat outdated. Typical settings for Gromos96 would be something like: rlist = 0.9 rcoulomb = 0.9 rvdw = 1.4 nstlist = 5 coulombtype = PME Note that the value of rcoulomb and rlist can vary a bit as a consequence of PME. -Justin Thanks a lot, Yun Shi On 04/08/11, Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu wrote: Yun Shi wrote: Hi all, I am working with GROMOS 53a6 ff in GROMACS 4.5, and I assume a Lennard-Jones interaction function was used for short-range vdw interactions. From the reference paper /A Biomolecular Force Field Based on the Free Enthalpy of Hydration and Solvation: The GROMOS Force-Field Parameter Sets 53A5 and 53A6/, I found that for example, when rvdw = 1.5nm, the repulsion term of the interaction between two CH1 type atoms (C12ij = 9.85^2) can be calculated as 9.850*9.850 / (1.5^12) = 0.747786 kJ/mol. So I wonder if this value is considered to be small enough to be ignored. You should pay attention to the column headings in table 7 so that you can compute the contribution correctly. However, the magnitude of the energy of any particular interaction is not really of any concern. The evolution of the system depends on the *forces*, and it is likely that the sum of the forces on any atom from all its repulsion interactions from atoms that are (say) 1.4nm to 1.5nm away is very close to zero, except in highly non-homogeneous spatial distributions of particles. In any case, the sum of that contribution will be much smaller than the other contributions. Mark In addition, it seems not until 5 nm does the dispersion term become larger than the repulsion term in this case, so would turning on Dispersion Correction between, say 1.5 to 5 nm introduce more errors than turning it off? You should use the cutoff described the authors of the force field, in this case rvdw=1.4. Unless you can demonstrate that by using a different value you can achieve superior results, stick with the specifics of parameterization. I have never seen ill effects of setting rvdw=1.4 and using dispersion correction with this force field. -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
[gmx-users] setting vdw cutoff with specific force-filed?
Hi all, I am working with GROMOS 53a6 ff in GROMACS 4.5, and I assume a Lennard-Jones interaction function was used for short-range vdw interactions. From the reference paper *A Biomolecular Force Field Based on the Free Enthalpy of Hydration and Solvation: The GROMOS Force-Field Parameter Sets 53A5 and 53A6*, I found that for example, when rvdw = 1.5nm, the repulsion term of the interaction between two CH1 type atoms (C12ij = 9.85^2) can be calculated as 9.850*9.850 / (1.5^12) = 0.747786 kJ/mol. So I wonder if this value is considered to be small enough to be ignored. In addition, it seems not until 5 nm does the dispersion term become larger than the repulsion term in this case, so would turning on Dispersion Correction between, say 1.5 to 5 nm introduce more errors than turning it off? Any suggestion would be appreciated! Thanks, Yun Shi -- 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 protonate a gmx trajectory?
Hi there, I am using GROMACS4.5.4 and the GROMOS 53A6 united-atom force field to do MD with my protein-ligand system. I was trying to add aliphatic hydrogens to the .gro and .trr or .xtc trajectory file with g_protonate, but the Fatal error was: Library file in current dir nor not found ffgmx2in defalt directories. Then re-name the directory in ../share/top/gmx2.ff to ffgmx2, and this time g_protonate seems to be able to find those .hdb files. But another error follows as: Reading frames from gro file 'Protein in water', 66743 atoms. Reading frame 0 time0.000 Opening force field file ffgmx2/aminoacids.hdb Opening force field file ffgmx2/atomtypes.atp Atomtype 1 Opening force field file ffgmx2/aminoacids.n.tdb Opening force field file ffgmx2/aminoacids.c.tdb Segmentation fault I searched for 'Segmentation fault', nothing relevant to g_protonate came up. So I wonder if anyone could kindly help me out? Or is there an easy alternative way to protonate the trajectory file? Thanks a lot for your attention. Yun -- 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
