[gmx-users] Re: GPU-based workstation
I have nothing to contribute regarding the CPUs, but note that last week Nvidia introduced the GTX 780 GPU (a Titan lite), and this week it is expected to introduce the GTX 770 one (probably an overclocked 680). Therefore, disregarding possible price issues, it seems that the GTX 680 may not be the best option. Date: Mon, 27 May 2013 14:14:51 +0400 From: James Starlight jmsstarli...@gmail.commailto:jmsstarli...@gmail.com Subject: Re: Re: Re: [gmx-users] GPU-based workstation To: Discussion list for GROMACS users gmx-users@gromacs.orgmailto:gmx-users@gromacs.org Message-ID: CAALQopymnZm5r+3q4LK_ncM=vots3fvukbeky4k3pk27n2x...@mail.gmail.commailto:CAALQopymnZm5r+3q4LK_ncM=vots3fvukbeky4k3pk27n2x...@mail.gmail.com Content-Type: text/plain; charset=ISO-8859-1 On Nvidia benchmarks I've found suggestions of using of the two 6 cores CPU for systems with the 2 GPU. Assuming that I'll be using two 680 GTX cards with 256 bits and 4gb ram (not a profesional nvidia cards like TESLA) what CPU's could give me the best performance- 1 i7 of 8 cores or 2 Xeons e5 with 6 cores ? Does it meaningful to use 2 separate CPU's with several nodes each for the 2 GPU's ? James -- 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] GPU
Message: 4 Date: Mon, 11 Jun 2012 15:54:39 +1000 From: Mark Abraham mark.abra...@anu.edu.au Subject: Re: [gmx-users] GPU To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: 4fd5881f.3040...@anu.edu.au Content-Type: text/plain; charset=ISO-8859-1; format=flowed On 11/06/2012 2:32 AM, ifat shub wrote: Hi, If I understand correctly, currently the Gromacs GPU acceleration does not support energy minimization. Is this so? Are there any plans to include it in the 4.6 version or in a later one (i.e. to allow, say, integrator = steep or cg in mdrun-gpu)? I would find such options extremely useful. EM is normally so quick that it's not worth putting much effort into accelerating it, compared to the CPU-months that are spent doing subsequent MD. Mark Currently, my main use of Gromacs entails running multiple minimizations on an ensemble of states. Moreover, these states are not obtained using molecular dynamics but rather using the Concoord algorithm. Therefore, for me the bottleneck is not md but rather minimizations (specifically, cg) and so their acceleration on GPUs would be very advantageous. If such usage is not totally idiosyncratic, I hope the development team would reconsider GPU accelerating also minimizations. I suspect this would not be technically too complex given the work already done on dynamics. Thanks, Ehud Schreiber. -- 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] a request regarding pdb2gmx and specbond.dat
Hi, The mechanism for specifying special bonds for the pdb2gmx program, using the specbond.dat file (manual section 5.6.7), seems to me not general enough. The characteristic length of the bond is specifiable, but the range of lengths for which a bond is assumed is hardcoded as +-10% of the characteristic length. I, for example, am perturbing structures using the Concoord method, and lengths of disulfide bonds sometimes vary by more than 10%. I therefore suggest that two values would be added at the end of the line declaring a special bond, specifying the minimum and maximum lengths. These two values would be optional, so as not to break the existing arrangement, which would remain as the default. This should not require more than a few extra lines of code and would be, in my opinion, very useful. What do you say? Thanks, Ehud Schreiber. -- 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] trjconv -dump problem
Dear Gromacs users, I minimized a protein structure 1IARcompleted_WT.pdb, getting, among others, the files 1IARcompleted_WT_minimized.trr and 1IARcompleted_WT_minimized_potential_energy.xvg . Looking at the latter file showed the last frame to be at 217 ps: . . @ s0 legend Potential 0.00 -31997.519531 0.00 -33810.406250 200.00 -69850.609375 217.00 -69898.031250 I wanted to extract only this last frame from the .trr file, so used trjconv -f 1IARcompleted_WT_minimized.trr -o 1IARcompleted_WT_minimized_217.trr -dump 217 However, this seems to have produced a file with the t = 200 ps conformation, though the dump parameter was recorded, as trjconv output was: . . Option Type Value Description -- . . -dumptime 217 Dump frame nearest specified time (ps) . . Will write trr: Trajectory in portable xdr format trn version: GMX_trn_file (single precision) Reading frame 2 time 200.000 Dumping frame at t= 200 ps Reading frame 3 time 217.000 . . Also, using -dump 200 gave an identical file to the above. Any idea why the expected timeframe isn't reproduced? I'm using gromacs 4.5.3. Thanks, Ehud Schreiber. -- 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] FW: trjconv -dump problem
Dear Gromacs users, I minimized a protein structure 1IARcompleted_WT.pdb, getting, among others, the files 1IARcompleted_WT_minimized.trr and 1IARcompleted_WT_minimized_potential_energy.xvg . Looking at the latter file showed the last frame to be at 217 ps: . . @title Gromacs Energies @xaxis label Time (ps) @yaxis label (kJ/mol) . . @ s0 legend Potential 0.00 -31997.519531 0.00 -33810.406250 200.00 -69850.609375 217.00 -69898.031250 I wanted to extract only this last frame from the .trr file, so used trjconv -f 1IARcompleted_WT_minimized.trr -o 1IARcompleted_WT_minimized_217.trr -dump 217 However, this seems to have produced a file with the t = 200 ps conformation, though the dump parameter was recorded, as trjconv output was: . . Option Type Value Description -- . . -dumptime 217 Dump frame nearest specified time (ps) . . Will write trr: Trajectory in portable xdr format trn version: GMX_trn_file (single precision) Reading frame 2 time 200.000 Dumping frame at t= 200 ps Reading frame 3 time 217.000 . . Also, using -dump 200 gave an identical file to the above. Any idea why the expected timeframe isn't reproduced? I'm using gromacs 4.5.3. Thanks, Ehud Schreiber. -- 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: trjconv -dump problem
Dear Mark (or anybody else interested), The .trr file does include the final (t = 217) frame - first, trjconv said: Reading frame 3 time 217.000 and second this is verified by converting the whole trajectory to .gro: trjconv -f 1IARcompleted_WT_minimized.trr -s 1IARcompleted_WT_minimized.tpr -o 1IARcompleted_WT_minimized_path.gro It therefore seems the behavior is a bug, as the last frame is there, and is needed, especially in minimizations. An output in a .gro format is not sufficient for a further minimization because of its limited accuracy format of coordinates. Also, my mdrun did not produce a checkpoint file (I'm not sure whether because I didn't ask to or because the run was shorter than 15 minutes). Thanks, Ehud Schreiber. -- Message: 2 Date: Tue, 14 Feb 2012 21:58:10 +1100 From: Mark Abraham mark.abra...@anu.edu.au Subject: Re: [gmx-users] FW: trjconv -dump problem To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: 4f3a3e42.6000...@anu.edu.au Content-Type: text/plain; charset=iso-8859-1 On 14/02/2012 7:59 PM, Ehud Schreiber wrote: Dear Gromacs users, I minimized a protein structure 1IARcompleted_WT.pdb, getting, among others, the files 1IARcompleted_WT_minimized.trr and 1IARcompleted_WT_minimized_potential_energy.xvg . Looking at the latter file showed the last frame to be at 217 ps: . . @title Gromacs Energies @xaxis label Time (ps) @yaxis label (kJ/mol) . . @ s0 legend Potential 0.00 -31997.519531 0.00 -33810.406250 200.00 -69850.609375 217.00 -69898.031250 I wanted to extract only this last frame from the .trr file, so used trjconv -f 1IARcompleted_WT_minimized.trr -o 1IARcompleted_WT_minimized_217.trr -dump 217 However, this seems to have produced a file with the t = 200 ps conformation, though the dump parameter was recorded, as trjconv output was: . . Option Type Value Description -- . . -dumptime 217 Dump frame nearest specified time (ps) . . Will write trr: Trajectory in portable xdr format trn version: GMX_trn_file (single precision) Reading frame 2 time 200.000 Dumping frame at t= 200 ps Reading frame 3 time 217.000 . . Also, using -dump 200 gave an identical file to the above. Any idea why the expected timeframe isn't reproduced? I'm using gromacs 4.5.3. Your nstxout parameter means not every frame is written. Prior to the implementation of checkpointing, the final frame was written to the .trr regardless of nstxout, but that no longer occurs. The final frame is in your checkpoint file, and you can use that anywhere you might use a coordinate file - including trjconv to get a simple coordinate file from it. 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] N-Acetylglucosamine (NAG) in OPLS-AA
Hi, I am interested in the structure of a protein where some asparagine residues are N-glycosylated by N-Acetylglucosamine (NAG). As a first step I tried to model NAG itself (see e.g. http://xray.bmc.uu.se/hicup/NAG/index.html) in the OPLS-AA force field using implicit solvation; this has involved a surprisingly (for me) large amount of tinkering. Here I'd like to present the changes I needed to make to the oplsaa.ff/* force field files in order to accomplish this; I would like to elicit comments, and hopefully this will also be useful to others. 1) Added to aminoacids.rtp : [ NAG ] ; N-Acetylglucosamine [ atoms ] C1opls_195 +0.3651 HC1opls_196 +0.1001 O1opls_187 -0.7001 HO1opls_188 +0.4351 C5opls_183 +0.1701 HC5opls_185 +0.0301 O5opls_180 -0.4001 C2opls_229 +0.1402 HC2opls_140 +0.0602 N2opls_238 -0.5002 HN2opls_241 +0.3002 C3opls_158 +0.2053 HC3opls_140 +0.0603 O3opls_154 -0.6833 HO3opls_155 +0.4183 C4opls_158 +0.2054 HC4opls_140 +0.0604 O4opls_154 -0.6834 HO4opls_155 +0.4184 C6opls_157 +0.1455 HC61opls_140 +0.0605 HC62opls_140 +0.0605 O6opls_154 -0.6835 HO6opls_155 +0.4185 C7opls_235 +0.5006 O7opls_236 -0.5006 C8opls_135 -0.1807 HC81opls_140 +0.0607 HC82opls_140 +0.0607 HC83opls_140 +0.0607 [ bonds ] C1 HC1 C1 O1 C1 C2 C1 O5 O1 HO1 C2 HC2 C2 N2 C2 C3 N2 HN2 N2 C7 C3 HC3 C3 O3 C3 C4 O3 HO3 C4 HC4 C4 O4 C4 C5 O4 HO4 C5 HC5 C5 O5 C5 C6 C6 HC61 C6 HC62 C6 O6 O6 HO6 C7 O7 C7 C8 C8 HC81 C8 HC82 C8 HC83 2) Added to aminoacids.hdb : NAG12 1 5 HC1 C1 O5 O1 C2 1 5 HC2 C2 C1 N2 C3 1 5 HC3 C3 C2 O3 C4 1 5 HC4 C4 C3 O4 C5 1 5 HC5 C5 C4 C6 O5 1 2 HO1 O1 C1 C2 1 2 HO3 O3 C3 C4 1 2 HO4 O4 C4 C5 1 2 HO6 O6 C6 C5 1 1 HN2 N2 C2 C7 2 6 HC6 C6 C5 O6 3 4 HC8 C8 C7 N2 3) Added to ffbonded.itp : [ angletypes ] ; i jk func th0 cth CO CT N 1 109.700669.440 ; copied from: CT CT N ; ALAJACS 94, 2657 [ dihedraltypes ] ; i j k l func coefficients HC CO OH HO 3 0.94140 2.82420 0.0 -3.76560 0.0 0.0 ; copied from: HC CT OH HO ; alcohols AA CT CO OH HO 3 -0.44350 3.83255 0.72801 -4.11705 0.0 0.0 ; copied from: CT CT OH HO ; alcohols AA HC CO CT HC 3 0.62760 1.88280 0.0 -2.51040 0.0 0.0 ; copied from: HC CT CT HC ; hydrocarbon *new* 11/99 HC CO CT N 3 0.97069 2.91206 0.0 -3.88275 0.0 0.0 ; copied from: HC CT CT N ; N-ethylformamide HC CT CO OH 3 0.97905 2.93716 0.0 -3.91622 0.0 0.0 ; copied from: HC CT CT OH ; alcohols, ethers AA C N CT CO 3 -4.70700 2.92044 1.78656 0.0 0.0 0.0 ; copied from: C N CT CT ; N-ethylformamide N CT CT OH 3 9.89307 -4.71746 3.67774 -8.85335 0.0 0.0 ; copied from : N CT_2 CT OH ; Chi for Ser Thr N CT CO OH 3 9.89307 -4.71746 3.67774 -8.85335 0.0 0.0 ; copied from : N CT_2 CT OH ; Chi for Ser Thr N CT CO OS 3 2.92880 -1.46440 0.20920 -1.67360 0.0 0.0 ; copied from: N3 CT CT OS ; Guess for lipids 4) Added to gbsa.itp : opls_180 0.152 1 1.0800.1535 0.85 ; copied from opls_154 ; opls_183 0.18 1 1.2760.1900.72 ; copied from opls_158 ; opls_185 0.11 10.1250.85 ; copied from opls_140 ; opls_187 0.152 1 1.0800.1535 0.85 ; copied from opls_154 ; opls_188 0.11 10.1050.85 ; copied from opls_155 ; opls_195 0.18 1 1.2760.1900.72 ; copied from opls_158 ; opls_196 0.11 10.1250.85 ; copied from opls_140 ; opls_229 0.18 1 1.2760.1900.72 ; copied from opls_158 ; So, are my choices reasonable? Could this have been done in some easier fashion? Thanks, Ehud
[gmx-users] RE: How to obtain structures with large RMSD?
Hi, You can try the non-gromacs tool Concoord: http://www.mpibpc.mpg.de/groups/de_groot/concoord/ http://151.100.52.62/pdf/046_amadei_proteins_1997.pdf It was used e.g. in the CC/PBSA algorithm http://ccpbsa.biologie.uni-erlangen.de/ccpbsa/ http://www.nature.com/nmeth/journal/v6/n1/full/nmeth0109-3.html -- Message: 1 Date: Mon, 10 Oct 2011 10:04:30 -0500 From: Liu, Liang liu4...@gmail.com Subject: [gmx-users] How to obtain structures with large RMSD? To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: CABxJRoJ5MgcYx02c5CqjtbU8HhEdh1=cyu1eggzzqeot0_+...@mail.gmail.com Content-Type: text/plain; charset=iso-8859-1 Hi, all, I am trying to use Gromacs to obtain structural ensembles around native structures (PDB structures). However the simulated structures are always very close to the initial one, with RMSD 0.2. I am wondering how to obtain large-RMSD structures? Thanks. -- Best, Liang Liu -- 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] opls-aa, implicit solvent and GPUs
Hi again, Quoting from the current Gromacs GPUs webpage: Force Fields: Supported FF are Amber, CHARMM. GROMOS and OPLS-AA are not supported. Also, it states Forcefields that do not use combination rules for Lennard-Jones interactions are not supported yet. And I don't know whether opls-aa is such a forcefield. So is or isn't it currently possible to accelerate opls-aa on GPUs? Looking at the Roadmap/GROMACS 4.6 webpage, on the up side: Major features (almost) ready to be merged into the 4.6 branch: New GPU non-bonded kernels which work with almost all features of GROMACS. While on the downside: Features currently not supported by the new GPU and SSE kernels: Implicit solvent (but this will still by supported on the GPU through OpenMM) All this led me to my original question. Thanks for your time (and for updating the Roadmap webpage...), Ehud. -- Message: 2 Date: Thu, 22 Sep 2011 10:55:56 -0400 From: Justin A. Lemkul jalem...@vt.edu Subject: Re: [gmx-users] opls-aa, implicit solvent and GPUs To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: 4e7b4c7c.1010...@vt.edu Content-Type: text/plain; charset=windows-1252; format=flowed Ehud Schreiber wrote: Hi, I'm currently utilizing gromacs using the opls-aa forcefield in an implicit-solvent setting. Can any of the developers roughly predict in what version, and when, will GPU acceleration be applicable for such work? I am afraid the gromacs 4.6 roadmap webpage left me confused. You can do these simulations on GPU with the current version, or several of the previous releases. -Justin By the way, the main roadmap webpage is totally outdated and needs updating and cleaning... Thanks, Ehud. -- 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] opls-aa, implicit solvent and GPUs
Hi, I'm currently utilizing gromacs using the opls-aa forcefield in an implicit-solvent setting. Can any of the developers roughly predict in what version, and when, will GPU acceleration be applicable for such work? I am afraid the gromacs 4.6 roadmap webpage left me confused. By the way, the main roadmap webpage is totally outdated and needs updating and cleaning... Thanks, Ehud. -- 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] pdb2gmx, disulfide bonds and PDB CONECT lines
Dear Gromacs users, I am perturbing a PDB structure in a certain way, and then wish to obtain a topology for Gromacs minimization, simulation and analysis. The structure contains some disulfide bonds; however, pdb2gmx guesses them incorrectly as the distances between the Sulfur atoms change. Even using the -ss option is not enough, as some bonds are not even prompted for interactive selection. The PDB structure contains some CONECT lines specifying the disulfide bonds, but those are just ignored by pdb2gmx. Can I convince pdb2gmx to use the CONECT lines or in any other way to create a topology with the desired disulfide bonds? Preferably the solution should be non-interactive as I handle multiple perturbed structures. The perturbations use an existing utility so must be performed on the PDB level. Thanks, Ehud Schreiber. -- 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] pdb2gmx with input topology
Dear Gromacs users, I wish to convert PDB files to gromacs format while specifying the topology; however, in the pdb2gmx program, a .top file can only be an output, not an input. The need arises as I have an NMR PDB structure (including hydrogens) which I rather violently perturb into multiple PDB structures using some non-gromacs means. When I convert those structures to .gro using pdb2gmx, sometimes the protonation state of Histidine is guessed incorrectly by the program (the hydrogen moves between HD1 and HE2). Does anybody have a suggestion what can I do? I need the process to be automatic and so can't use the interactive -his option. Thanks, Ehud. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] RE: g_energy inconsistent results
Dear Mark Abraham (and anyone else interested), I have implemented your suggestion, changing in the status.mdp file to integrator = md and adding continuation = yes (this is the new name of the unconstrained_start parameter); however, this did not have any effect - the energy remained as before, different from the one obtained during the minimization. I have then encountered another (perhaps related) issue. I thought that the problem may arise from the combination of the generalized Born and all-vs-all settings. I have therefore changed in the minimization to rgbradii = rlist = rcoulomb = rvdw = 100 (from = 0). As my system is far smaller than 100 nm, I expected these parameters to provide also an all-vs-all setting (even if non-optimized). Nevertheless, the resulting structure differed from the previous minimized one (rmsd ~ 0.005 nm, delta energy ~ a few kJ/mol). Can this arise from rounding effects only or does this signify a problem? I'm not sure what rgbradii does, but the manual states that it must equal rlist. In addition, changing also status.mdp to have all radii = 100 and running on the new optimized output again does not yield the same energy as during the optimization. Does all of this make any sense to you? Thank, Ehud. Message: 6 Date: Tue, 08 Mar 2011 23:37:12 +1100 From: Mark Abraham mark.abra...@anu.edu.au Subject: Re: [gmx-users] g_energy inconsistent results To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: 4d7622f8.7050...@anu.edu.au Content-Type: text/plain; charset=iso-8859-1 On 8/03/2011 9:44 PM, Ehud Schreiber wrote: Dear Gromacs users, I am working with version 4.5.3, using the opls-aa forcefield in an implicit solvent, all-vs-all setting: pdb2gmx -ter -ff oplsaa -water none -f file.pdb I am energy-minimizing structures in 3 stages (steep, cg and l-bfgs). The last stage is the following: grompp -f em3.mdp -p topol.top -c em2.gro -t em2.trr -o em3.tpr -po em3.mdout.mdp mdrun -nice 0 -v -pd -deffnm em3 g_energy -s em3.tpr -f em3.edr -o em3.potential_energy.xvg where the mdp file is: ;;; em3.mdp ;;; integrator = l-bfgs nsteps = 5 implicit_solvent = GBSA gb_algorithm = Still sa_algorithm = Ace-approximation pbc = no rgbradii = 0 ns_type = simple nstlist = 0 rlist= 0 coulombtype = cut-off rcoulomb = 0 vdwtype = cut-off rvdw = 0 nstcalcenergy= 1 nstenergy= 1000 emtol= 0 ;;; The last line in the em3.potential_energy.xvg file should give the (potential) energy of the minimized structure em3.gro . I wish also to compute the potential energy of .gro files in general, not necessarily obtained from a simulation. For that, I prepared a .mdp file for a degenerate energy minimization, having 0 steps, designed just to give the status of the file: Zero-step EM does not calculate the initial energy because it is not useful for gradient-based energy minimization. I don't recall the details, but perhaps the first EM step is reported as step zero. Instead, you should use zero-step MD (with unconstrained_start = yes), or (for multiple single points) mdrun -rerun. You will not necessarily reproduce the g_energy energies with anything, because the energy is dependent on the state of the neighbour lists. If nstenergy is a multiple of nstlist, then those energies should be fairly reproducible. I have updated the grompp source to provide a note to the user to warn against zero-step EM. Mark ;;; status.mdp ;;; integrator = l-bfgs nsteps = 0 implicit_solvent = GBSA gb_algorithm = Still sa_algorithm = Ace-approximation pbc = no rgbradii = 0 ns_type = simple nstlist = 0 rlist= 0 coulombtype = cut-off rcoulomb = 0 vdwtype = cut-off rvdw = 0 nstcalcenergy= 1 nstenergy= 1 emtol= 0 ;;; The only changes from the former .mdp file are in nsteps and nstenergy. However, when I run this potential energy status run on em3.gro itself, grompp -f status.mdp -p topol.top -c em3.gro -o status.tpr -po status.mdout.mdp mdrun -nice 0 -v -pd -deffnm status g_energy -s status.tpr -f status.edr -o status.potential_energy.xvg and look at the (single) energy line in status.potential_energy.xvg I find that the energy does not agree with the one obtained during minimization (it's higher by some tens of kJ/mol). What am I doing wrong? How should one reliably find the energy of a given .gro file? Moreover, when changing in status.mdp to integrator = steep, the results also change dramatically -- why should the algorithm matter if no steps
[gmx-users] g_energy inconsistent results
Dear Gromacs users, I am working with version 4.5.3, using the opls-aa forcefield in an implicit solvent, all-vs-all setting: pdb2gmx -ter -ff oplsaa -water none -f file.pdb I am energy-minimizing structures in 3 stages (steep, cg and l-bfgs). The last stage is the following: grompp -f em3.mdp -p topol.top -c em2.gro -t em2.trr -o em3.tpr -po em3.mdout.mdp mdrun -nice 0 -v -pd -deffnm em3 g_energy -s em3.tpr -f em3.edr -o em3.potential_energy.xvg where the mdp file is: ;;; em3.mdp ;;; integrator = l-bfgs nsteps = 5 implicit_solvent = GBSA gb_algorithm = Still sa_algorithm = Ace-approximation pbc = no rgbradii = 0 ns_type = simple nstlist = 0 rlist= 0 coulombtype = cut-off rcoulomb = 0 vdwtype = cut-off rvdw = 0 nstcalcenergy= 1 nstenergy= 1000 emtol= 0 ;;; The last line in the em3.potential_energy.xvg file should give the (potential) energy of the minimized structure em3.gro . I wish also to compute the potential energy of .gro files in general, not necessarily obtained from a simulation. For that, I prepared a .mdp file for a degenerate energy minimization, having 0 steps, designed just to give the status of the file: ;;; status.mdp ;;; integrator = l-bfgs nsteps = 0 implicit_solvent = GBSA gb_algorithm = Still sa_algorithm = Ace-approximation pbc = no rgbradii = 0 ns_type = simple nstlist = 0 rlist= 0 coulombtype = cut-off rcoulomb = 0 vdwtype = cut-off rvdw = 0 nstcalcenergy= 1 nstenergy= 1 emtol= 0 ;;; The only changes from the former .mdp file are in nsteps and nstenergy. However, when I run this potential energy status run on em3.gro itself, grompp -f status.mdp -p topol.top -c em3.gro -o status.tpr -po status.mdout.mdp mdrun -nice 0 -v -pd -deffnm status g_energy -s status.tpr -f status.edr -o status.potential_energy.xvg and look at the (single) energy line in status.potential_energy.xvg I find that the energy does not agree with the one obtained during minimization (it's higher by some tens of kJ/mol). What am I doing wrong? How should one reliably find the energy of a given .gro file? Moreover, when changing in status.mdp to integrator = steep, the results also change dramatically - why should the algorithm matter if no steps are performed and the initial structure is explored? Thanks, Ehud. -- 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] Simulated Annealing parameters
Dear fellow GROMACS users, I am using the simulated annealing option in gromacs 4.5.2 (in an implicit solvent all-vs-all setting, by the way), i.e. using the annealing, annealing_npoints, annealing_time, annealing_temp parameters in the mdp file. I also specify the tcoupl, tc_grps and tau_t parameters for the temperature coupling, as well as the gen_vel, gen_temp, gen_seed ones for initial velocity generation. The perplexing behavior I'm encountering is that grompp requires me to specify also the ref_t parameter although it would seem to be overridden by the annealing parameters mentioned above. Am I doing anything wrong? If not, the correct behavior would be to warn when this parameter is specified rather than to demand it (and then presumably ignore it). Thanks in advance for your thoughts on the matter, Ehud Schreiber. -- 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] Generalized Born in GROMACS manual
Dear GROMACS users,
I have looked into some of the literature in order to understand better
implicit solvent models, in particular the generalized Born formalism. I
believe that section 3.18 of the v4.5 GROMACS manual, dealing with the
above, is in need of some corrections:
1) In eq. 3.135, the sum should run over all pairs i,j, in
particular allowing i=j; there should be a corresponding factor of 1/2,
and, depending on the conventions used, perhaps also a minus sign.
See e.g. eq. 5 in Still et al., J. Am. Chem. Soc. 112 (1990),
http://pubs.acs.org/doi/pdf/10.1021/ja00172a038.
2) In eqs. 3.137 and 3.138, x should be replaced by x_{ij}.
3) It might be better to add explicitly that \xi(x) = 1/\sqrt{x^2 +
exp(-x^2/4)}, and that \xi(x) ~ 1/x for x 1 and \xi(x) ~ 1 for 0 = x
1.
What do you think?
--
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: Gibbs free energy of binding
Hi Moshen, I think everybody agrees that a full calculation such as Free Energy Perturbation is the accurate, if difficult and lengthy, approach. The entropic effects usually cannot simply be ignored. All I tried to say was that there are approximation schemes for these (see the reference below). Still, I would trust such approximations only when computing binding Delta Delta G between two close variants (e.g. a wild type protein and a one residue mutation) such that most entropic contributions would tend to cancel. Ehud. -- Date: Thu, 21 Oct 2010 22:45:23 +0330 From: mohsen ramezanpour ramezanpour.moh...@gmail.com Subject: Re: [gmx-users] RE: Gibbs free energy of binding reading your idea: it seems to me I can't ignore entropy contribution because my simulation is at room tempreture. Really I couldn't understand what can I do! I am working at room tempreture and I want to estimate binding free energy(delta G),can I ignore entropy in this simulation and calculate binding free energy by the method that I said in my last email? what do you think? thank in advance for your guid On Thu, Oct 21, 2010 at 12:15 PM, David van der Spoel sp...@xray.bmc.uu.sewrote: On 2010-10-21 10.39, Ehud Schreiber wrote: Actually, I believe that using the energy difference, Delta E, as an approximation to the free energy difference, Delta G, is a valid approach (which I'm considering myself). The entropic contribution to Delta G, namely -T Delta S, may be less prominent than Delta E. In addition, Delta S can be approximated by various means - see e.g. Doig Sternberg 1995. I understand that such an approach is utilized in the Accelrys Discovery Studio. Obviously, this is an approximation that might be too crude for some applications. -- 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: Gibbs free energy of binding
Actually, I believe that using the energy difference, Delta E, as an approximation to the free energy difference, Delta G, is a valid approach (which I'm considering myself). The entropic contribution to Delta G, namely -T Delta S, may be less prominent than Delta E. In addition, Delta S can be approximated by various means - see e.g. Doig Sternberg 1995. I understand that such an approach is utilized in the Accelrys Discovery Studio. Obviously, this is an approximation that might be too crude for some applications. What do you think? -- On Oct 21, 2010, at 09:25 , Sander Pronk wrote: Hi Mohsen, The mean energy difference is only one component of the free energy difference. Before you go any further I'd suggest reading a good book on molecular simulations, like 'Understanding Molecular Simulations' by Frenkel and Smit. There's a good reason free energy calculations cover over half of that book. Sander On Oct 21, 2010, at 09:18 , mohsen ramezanpour wrote: Dear Justin If I do two MD simulations for a short time in the same conditions(of course separately for protein and drug) and calculate total energy of each one and sum them with each other as E1 as nonbonding free energy of system. then a MD simulation for Protein-drug system in the same condition and calculate it's total energy too as E2 as bound system . what does (E1-E2)mean? I think it is binding free energy,Is not it? in the other hand when we are working on NPT ensamble it means Gibbs free energy is the main energy and our total energy is equal to Gibbs free energy. Then,what is the problem? -- 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: Energy Minimization
Hi, I would allow flexible water by uncommenting the line define = -DFLEXIBLE and trying again; reducing the time step might also be required. Message: 4 Date: Mon, 18 Oct 2010 16:26:16 -0400 From: Sai Pooja saipo...@gmail.com Subject: [gmx-users] Energy Minimization To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: aanlktinomj-hzw2u+g6sa=mxqlc7q6ibwwiqj+mx6...@mail.gmail.com Content-Type: text/plain; charset=iso-8859-1 Hi all, I am using CHARMM forcefield with tables for Protein SOL interactions for alanine dipeptide. Tables supplied: table.xvg, table_Protein_SOL.xvg, tablep.xvg (Standard tables for 6-12 interactions used) Combination rule changed from '2' to '1' in forcefield.itp file. In the energy minimization step, using mdrun the following problem is encountered: Polak-Ribiere Conjugate Gradients: Tolerance (Fmax) = 1.0e+00 Number of steps=1 F-max = 2.98523e+10 on atom 4 F-Norm= 1.32072e+09 step -1: Water molecule starting at atom 1149 can not be settled. Check for bad contacts and/or reduce the timestep if appropriate. Wrote pdb files with previous and current coordinates title = Energy Minimization ; Title of run ; The following line tell the program the standard locations where to find certain files cpp = /lib/cpp ; Preprocessor ; Define can be used to control processes ;define = -DFLEXIBLE define = -DPOSRES ; Parameters describing what to do, when to stop and what to save integrator = cg; Algorithm (steep = steepest descent minimization) emtol = 1.0 ; Stop minimization when the maximum force 1.0 kJ/mol nsteps = 1 ; Maximum number of (minimization) steps to perform nstenergy = 1 ; Write energies to disk every nstenergy steps energygrps = Protein SOL energygrp_table = Protein SOL ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions ns_type = grid ; Method to determine neighbor list (simple, grid) coulombtype = User ; Treatment of long range electrostatic interactions rcoulomb= 1.0 ; long range electrostatic cut-off rvdw= 1.0 ; long range Van der Waals cut-off constraints = none ; Bond types to replace by constraints pbc = xyz ; Periodic Boundary Conditions (yes/no) Any suggestions? Pooja -- Quaerendo Invenietis-Seek and you shall discover. -- next part -- An HTML attachment was scrubbed... URL: http://lists.gromacs.org/pipermail/gmx-users/attachments/20101018/355f3a 45/attachment.html -- -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! End of gmx-users Digest, Vol 78, Issue 132 ** -- 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] best parameters for energy determination with GBSA
Dear GROMACS users, I try to compute energies of a small peptide-domain complex and of its constituents using implicit solvent, so that I can get the binding Delta E. I need the accuracy to be high enough so that I can reliably measure the Delta Delta E between the wildtype peptide and a mutation of it. I would like to get advice as to the best setting and parameters to use; are the ones given below the correct and best ones (most in an mdp file, some run control parameters of GROMACS commands)? Are other defaults in need to be changed? As the GBSA method necessitates using the inexact cut-off interactions, and my system is rather small, I consider using no cut-off at all, i.e. an all-vs.-all approach. For that I understand the parameters should be nstlist=0 ns_type=simple -pd ; in mdrun coulombtype=cut-off rcoulomb=0 vdwtype=cut-off rvdw=0 rlist=0 Then, there are the implicit solvent parameters themselves: implicit_solvent=GBSA gb_algorithm= Still sa_algorithm=Ace-approximation rgbradii=0 Finally, regarding the simulation box, I see no need for periodic boundary conditions, and I pick some big margin around the system (as there are no water molecules, no need to spare): pbc=no -d 5.0 ; in editconf What do you say? Thanks, Ehud. -- 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] potential energy in implicit solvent simulations
Dear GROMACS users, I am conducting energy minimization with an implicit solvent, using GROMACS version 4.5.1 at double precision. The minimization provides the minimal potential energy I'm interested in, but I also wanted to see the contributions to it. I therefore used g_energy_d and picked the potential energy itself (no. 10) together with what I thought are its constituents (no. 1-9). However, those contributions do not sum to the potential energy. There is also a term #Surf*SurfTen (no. 30) but it is printed as zero. Can you explain this situation? The relevant part of g_energy_d's output for a certain example is the following: # This file was created Thu Sep 16 16:51:26 2010 # by the following command: # g_energy_d -s em2.tpr -f em2.edr -o em2.potential_energy.xvg # # g_energy_d is part of G R O M A C S: # # GROningen Mixture of Alchemy and Childrens' Stories # @title Gromacs Energies @xaxis label Time (ps) @yaxis label (kJ/mol), (bar nm) @TYPE xy @ view 0.15, 0.15, 0.75, 0.85 @ legend on @ legend box on @ legend loctype view @ legend 0.78, 0.8 @ legend length 2 @ s0 legend Bond @ s1 legend Angle @ s2 legend Proper Dih. @ s3 legend Ryckaert-Bell. @ s4 legend GB Polarization @ s5 legend LJ-14 @ s6 legend Coulomb-14 @ s7 legend LJ (SR) @ s8 legend Coulomb (SR) @ s9 legend Potential @ s10 legend #Surf*SurfTen 0.00 179.387112 1063.328043 64.818404 2220.374294 440.027560 3580.261354 12299.937527 -3164.489764 -44972.543752 -27990.3807600.00 -- 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] grompp_d and mdrun_d issues with GBSA and normal mode analysis
Dear GROMACS users, I am encountering a couple of issues when trying to perform normal mode analysis in an implicit solvent (GBSA) setting. I am using version 4.5.1 with double precision; unfortunately I do not have the single precision version installed for comparison. The starting point is a small protein which was energy minimized in two stages, also in double precision and with GBSA, producing em2 files. The mdp file used for the normal mode analysis is the following: - nm.mdp - integrator = nm nsteps = 1 implicit_solvent = GBSA gb_algorithm = Still ; the default rgbradii = 1.0 ; must be equal to rlist rlist= 1.0 coulombtype = cut-off rcoulomb = 1.0 vdwtype = cut-off rvdw = 1.0 --- The command used is grompp_d -f nm.mdp -p topol.top -c em2.gro -t em2.trr -o nm.tpr which ran fine except for the following note: NOTE 1 [file nm.mdp]: You are using a plain Coulomb cut-off, which might produce artifacts. You might want to consider using PME electrostatics. Then, I tried to run mdrun_d -v -nice 0 -deffnm nm However, the command seems to be stuck, and the following serious warning is produced: Warning: 1-4 interaction between 64 and 71 at distance 3.114 which is larger than the 1-4 table size 2.000 nm These are ignored for the rest of the simulation This usually means your system is exploding, if not, you should increase table-extension in your mdp file or with user tables increase the table size Maximum force: 3.48709e+06 Maximum force probably not small enough to ensure that you are in an energy well. Be aware that negative eigenvalues may occur when the resulting matrix is diagonalized. The em2.gro is definitely not having such a large distance between atoms 64 and 71: 6GLN CB 64 1.879 1.895 2.423 . . 6GLNOE1 71 2.047 1.715 2.642 which are only about 0.346 nm apart. Also, the energy minimization, which used the same options, mutatis mutandis, -- em2.mdp -- integrator = cg nsteps = 1000 nstcgsteep = 40 implicit_solvent = GBSA gb_algorithm = Still ; the default rgbradii = 1.0 ; must be equal to rlist nstlist = 10 rlist= 1.0 coulombtype = cut-off rcoulomb = 1.0 vdwtype = cut-off rvdw = 1.0 nstenergy= 10 - converged to machine precision having a much smaller maximum force: Stepsize too small, or no change in energy. Converged to machine precision, but not to the requested precision Fmax 10 Polak-Ribiere Conjugate Gradients converged to machine precision in 0 steps, but did not reach the requested Fmax 10. Potential Energy = -2.71868395521758e+04 Maximum force = 4.63436548427712e+02 on atom 584 Norm of force = 1.25296847735530e+02 The other problem arises when I try to follow the grompp_d note above and change in nm.mdp to coulombtype = pme I then have a fatal error: ERROR 1 [file nm.mdp]: With GBSA, coulombtype must be equal to Cut-off Am I doing something wrong or are there still some problems in the new GBSA option of version 4.5.1? Thanks, Ehud Schreiber. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] RE: [nonbond_params] section , choice of parameters
Hi Amir,
Are you sure you're giving the parameters indeed as C_6 and C_{12}?
According to the combination rule in the [ defaults ] section, you may
be giving the two parameters directly as sigma and epsilon; sigma=0 (as
well as epsilon=0) would give an identically zero potential . See
sections 5.3.3 and 5.7.1 in the manual (version 4.0).
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] RE: RE: [nonbond_params] section , choice of parameters
Hi again,
Yes, you can simply specify sigma and epsilon, but indeed in your case
you do not want to do so but rather to specify C_6 and C_{12} with the
former being zero to avoid attraction.
Your .top file should include a [ defaults ] section, perhaps through an
#include of a .itp file. Perhaps yours is in ffamber99.itp (I don't have
it installed so I can't tell). In this section there are 5 parameters;
the second is a combination rule which must equal 1 for the two
numbers to represent indeed C_6 and C_{12} (see table 5.3, eq. 5.1-5.4,
and in page 112).
In short, check whether the second parameter in the [ defaults ] section
of ffamber99.itp is 1, if not that's probably the problem. Of course, if
you change this you'll screw up all the other parameters, so it seems
you could either change the force field or create your copy of the files
with all the relevant LJ parameters changed to the combination rule=1
convention.
From: Ehud Schreiber
Sent: Wednesday, July 28, 2010 3:01 PM
To: 'gmx-users@gromacs.org'
Subject: RE: [nonbond_params] section , choice of parameters
Hi Amir,
Are you sure you're giving the parameters indeed as C_6 and C_{12}?
According to the combination rule in the [ defaults ] section, you may
be giving the two parameters directly as sigma and epsilon; sigma=0 (as
well as epsilon=0) would give an identically zero potential . See
sections 5.3.3 and 5.7.1 in the manual (version 4.0).
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at http://www.gromacs.org/search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] HA1/HA2 chirality of GLY
Dear GROMACS users, Is there a convention about the chirality of the HA1 and HA2 hydrogen atoms in Glycine (relative to the C-alpha)? Glycine is non-chiral, but if one tries to mutate it to another L amino acid, one of these hydrogens should be picked as the analog of a sidechain. Which one is it? As most PDB files do not contain hydrogens, this question seems to be about the behavior of pdb2gmx when creating a .gro file with an all-atom force field (such as opls-aa). Thanks, Ehud Schreiber. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] large forces and monstrous water molecules in energy minimization step
Dear GROMACS users, When trying to simulate a pair of interacting proteins in water, I have encountered problems that ultimately resulted in the simulation crashing. I then tried to simplify the system as far as possible while retaining the problem; I now believe the problem (or at least a part of it) lies in the energy minimization step (the first molecular dynamics one). Specifically, the forces encountered during this step are very large, and some water molecules (which are supposed to be rigid) become giant and misshapen. In more details: 1) I use GROMACS version 4.0.7, single precision, on a server with two Intel x86-64 processors and the redhat 5.4 linux OS. 2) I created a PDB file, called NaCl.pdb, with only two atoms, actually Na+ and Cl- ions separated by their distance in the lattice of salt: HET NA A 1 1 HET CL A 1 1 HETNAM NA SODIUM ION HETNAM CL CHLORIDE ION FORMUL 1 NANA 1+ FORMUL 2 CLCL 1- HETATM1 NANA A 1 -1.410.0 0.01.00 0.0 NA HETATM2 CLCL A 1 +1.410.0 0.01.00 0.0 CL 3) I use the tip3p water model: pdb2gmx -f NaCl.pdb -water tip3p 4) I create the box: editconf -f conf.gro -bt dodecahedron -d 3.3 -o box.gro 5) I add water using spc216, creating the saltwater.gro file (which seems O.K. by inspection): genbox -cp box.gro -cs spc216.gro -p topol.top -o saltwater.gro 6) I create the energy minimization parameter file em.mdp: --em.mdp-- integrator = steep nsteps = 200 nstlist = 10 rlist = 1.0 coulombtype = pme rcoulomb= 1.0 vdwtype = Cut-off rvdw= 1.0 nstenergy = 10 -- 7) I prepare the em.tpr file for the energy minimization run: grompp -f em.mdp -p topol.top -c saltwater.gro -o em.tpr 8) I run the energy minimization step: mdrun -v -deffnm em 9) Looking at the em.log file I see that this step converged to machine precision but did not have maximal force 10: ... Enabling SPC water optimization for 7561 molecules. ... Max number of connections per atom is 2 Total number of connections is 30244 Max number of graph edges per atom is 2 Total number of graph edges is 30244 Going to use C-settle (7561 waters) wo = 0.33, wh =0.33, wohh = 3, rc = 0.075695, ra = 0.0390588 rb = 0.0195294, rc2 = 0.15139, rone = 1, dHH = 0.15139, dOH = 0.09572 ... Stepsize too small, or no change in energy. Converged to machine precision, but not to the requested precision Fmax 10 ... Steepest Descents converged to machine precision in 36 steps, but did not reach the requested Fmax 10. Potential Energy = -3.4678925e+05 Maximum force = 6.4623531e+05 on atom 11052 Norm of force = 5.4643726e+03 10) Looking at the em.gro file I see one monstrous water molecule (no. 3686); e.g., it has |HW2-OW| = 3.384876 nm, while the normal distance is about 0.1 nm. Its HW2 atom (no. 11054) is close to another water molecule (no. 5849), e.g., 0.047 nm from the latter's HW2 atom (no. 17543): ... 3686SOL OW11052 4.348 3.778 -0.629 3686SOLHW111053 5.360 2.601 0.505 3686SOLHW211054 6.518 1.650 0.861 ... 5849SOL OW17541 6.525 1.698 0.900 5849SOLHW117542 6.606 1.649 0.918 5849SOLHW217543 6.481 1.648 0.832 ... 11) During the simulation, several files called stepnnl.pdb were produced for problematic steps, where nn=11,15,19 and l=b,c. For example, the file step19c.pdb indeed shows a problematic water molecule no. 3686, while step19b.pdb does not. Likewise, the earlier step11c.pdb shows a problematic water molecule no. 3266 while step11b.pdb seems proper. The stdout/stderr of mdrun contains warnings like the following: ... t = 0.019 ps: Water molecule starting at atom 11052 can not be settled. Check for bad contacts and/or reduce the timestep. ... 12) Reducing the neighbor list update time, i.e., setting nstlist = 1, does not produce any change. 13) Trying to use the conjugate gradient integrator instead of steepest descent, i.e., setting integrator = cg, is even worse - the running crashes: ... Step 16, Epot=-1.258771e+35, Fnorm= nan, Fmax= inf (atom 14493) Segmentation fault Exit 139 So, am I doing something wrong? How can I avoid these problems? Thanks, Ehud Schreiber. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
