[gmx-users] Errors from alchemical_analysis.py from analyzing .xvg files
Hello, This request is related to python programs used for analyzing gromacs outputs, not gromacs itself. If there is a more appropriate forum for this question, let me know and I apologize for posting here. I have performed a bunch of thermodynamic-integration simulations for liquid water in Gromacs, and am trying to analyze the .xvg files with the python script mentioned in the ethanol solvation tutorial: http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy (I think the tutorial is out of date - the python script is now called alchemical_analysis.py and I think it is called with slightly different flags). It may be relevant that I'm using Windows 10 and installed Anaconda2 (64 bit) with the gui installer and then pymbar through Anaconda2 before downloading alchemical-analysis-master. In a Windows command prompt, I entered: >python alchemical_analysis.py -d C:\\ -t 200 -p 1 -v > results This is the result: Traceback (most recent call last): File "alchemical_analysis.py", line 1224, in main() File "alchemical_analysis.py", line 1166, in main nsnapshots, lv, dhdlt, u_klt = parser_gromacs.readDataGromacs(P) File "C:\<...>\alchemical-analysis-master\alchemical_analysis\parser_gromacs.py", line 187, in readDataGromacs fs = sorted(fs, key=F.sortedHelper) File "C:\<...>\alchemical-analysis-master\alchemical_analysis\parser_gromacs.py", line 48, in sortedHelper self.state = l[0] = int(l[0]) # Will be of use for selective MBAR analysis. IndexError: list index out of range where <...> is just the directory to which I downloaded alchemical-analysis-master. I don't understand the error and cannot find help for it online. In the directory where I have the .xvg files, all the files are simply named integers: 0.xvg, 1.xvg, 2.xvg, 20.xvg, and there are no other files in that folder. I don't want to tinker with the python code and accidentally mess something up. Am I calling the script correctly? Could it be the operating system (Windows vs. Linux)? I might be able to get on a Linux machine if it would help. Thanks very much, N. H. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Free energy of ice crystals (Thermodynamic Integration)
Hello, I would like to perform some free energy calculations of ice Ic (cubic ice) at various temperatures. However, I'm having difficulty finding information on how Gromacs can be used to do this. The goal is to gradually turn off VDW and Coulomb interactions, just as one would do in a liquid, but also turn on virtual spring-like potentials to create an Einstein crystal, which is where I'm having trouble. I know you can fix atom locations, and considering bond potentials are spring-like, my best idea is to create a virtual copy of each molecule (not virtual sites, but just atoms that have no charge and don't interact with other atoms/molecules), fix it in the crystal lattice, and create bonds between corresponding atoms with an equilibrium distance of 0. But this means I have to simulate twice as many atoms, and more importantly I don't know how to couple only some of a molecule's bonds with lambda. Am I even on the right track, or is there a more straightforward way to do this? Ar e there any tutorials for free energies of solids that I could take a look at? (I couldn't find any) Thanks for your help! N. H. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Free Energy of Liquid Water
You were right, it's much faster when I only calculate free energy for one molecule. My simulation has 1000 water molecules in it, but I created a new topology file and gave one molecule a different name (but identical properties) and now it runs the same simulation in about 36% more time, which is much more manageable. Thank you very much for the help! I don't think there is a problem with my topology or mdp files now. Regards, Nathan H. - Original Message - From: "Michael Shirts" To: "Discussion list for GROMACS users" Sent: Wednesday, October 7, 2015 11:14:29 PM Subject: Re: [gmx-users] Free Energy of Liquid Water If ALL the particles are changing with the free energy coupling parameter, then GROMACS will slow down quite a bit. If only one molecule is changing, then it shouldn't be that much slower (20-30% slower?) But without .mdp and.top files, it's hard to say exactly what is happening. On Wed, Oct 7, 2015 at 11:08 PM, Nathan K Houtz wrote: > Thanks for explaining, (and Professor Farais de Moura as well) that makes > sense. > > On a different note, I've noticed a huge difference in performance between my > equillibration run and my actual simulation, and I'm wondering if this is > normal. For comparison, I ran a couple small (meaningless) simulations where > the only difference was the collection of free energy information and Gromacs > slows down by a factor of 39. I had been under the impression that I would be > able to do each simulation in just under an hour but now it looks like it > will take the better part of two days, which is too much. I'll just have to > run shorter or fewer simulations if I can't improve that, but I'm hoping I'm > just doing something wrong. I'm essentially using modified versions of the > .mdp file I found for this tutorial: > http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy. > I changed things like turning off the pressure coupling, changing the > temperature, etc. The settings for free energy calculations are unchanged. Is > the calculation of dh/dl re > ally that expensive, or is this abnormal? > > Thank you for the help as always. Regards, > Nathan H. > > - Original Message - > From: "Michael Shirts" > To: "Discussion list for GROMACS users" > Cc: "gromacs org gmx-users" > Sent: Tuesday, October 6, 2015 3:05:38 PM > Subject: Re: [gmx-users] Free Energy of Liquid Water > > For a pure fluid, G = N \mu. And \mu = (dG/dN)_(T,P). So you only > need to change one molecule to ideal gas to get the change in free > energy. The free energy of transfer of water from liquid to gas is > indeed the free energy of solvation of one water molecule in bath of > water. So there's a reason why you're just finding tutorials of > solvation free energies! > > On Thu, Oct 1, 2015 at 10:44 PM, Nathan K Houtz wrote: >> Hi everyone, >> >> I would like to use Gromacs to do Thermodynamic Integration (TI) from liquid >> water (TIP4P model) to an ideal gas, to find the relative free energy. To do >> this, I believe one generally integrates above the critical point by >> increasing the temperature above the critical temperature and then relaxing >> the pressure until the system is a diffuse gas. The mdp options >> documentation is helpful, and I went through an ethanol solvation tutorial, >> but there doesn't appear to be a "pressure-lambda" or a "volume-lambda" >> option that I could use to do the second part. How can I get Gromacs to >> calculate the dh/dl derivative while relaxing the pressure? >> >> In addition, all of the tutorials I found for thermodynamic integration were >> for finding solvation free energies. The coulomb and VDW forces are >> essentially changed from "completely on" to "completely off". But in my >> case, I'd like to change the temperature and pressure between two nonzero >> values. I don't want to begin my simulation at 0K and 0atm, but lambda >> *must* go from 0 to 1. How can I define both starting and ending points for >> the temperature and pressure (or volume, or density)? >> >> Thanks for your help! >> Nathan H. >> -- >> Gromacs Users mailing list >> >> * Please search the archive at >> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! >> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >> >> * For (un)subscribe requests visit >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a >> mail to gmx-users-requ...
Re: [gmx-users] Free Energy of Liquid Water
Thanks for explaining, (and Professor Farais de Moura as well) that makes sense. On a different note, I've noticed a huge difference in performance between my equillibration run and my actual simulation, and I'm wondering if this is normal. For comparison, I ran a couple small (meaningless) simulations where the only difference was the collection of free energy information and Gromacs slows down by a factor of 39. I had been under the impression that I would be able to do each simulation in just under an hour but now it looks like it will take the better part of two days, which is too much. I'll just have to run shorter or fewer simulations if I can't improve that, but I'm hoping I'm just doing something wrong. I'm essentially using modified versions of the .mdp file I found for this tutorial: http://www.alchemistry.org/wiki/GROMACS_4.6_example:_Direct_ethanol_solvation_free_energy. I changed things like turning off the pressure coupling, changing the temperature, etc. The settings for free energy calculations are unchanged. Is the calculation of dh/dl re ally that expensive, or is this abnormal? Thank you for the help as always. Regards, Nathan H. - Original Message - From: "Michael Shirts" To: "Discussion list for GROMACS users" Cc: "gromacs org gmx-users" Sent: Tuesday, October 6, 2015 3:05:38 PM Subject: Re: [gmx-users] Free Energy of Liquid Water For a pure fluid, G = N \mu. And \mu = (dG/dN)_(T,P). So you only need to change one molecule to ideal gas to get the change in free energy. The free energy of transfer of water from liquid to gas is indeed the free energy of solvation of one water molecule in bath of water. So there's a reason why you're just finding tutorials of solvation free energies! On Thu, Oct 1, 2015 at 10:44 PM, Nathan K Houtz wrote: > Hi everyone, > > I would like to use Gromacs to do Thermodynamic Integration (TI) from liquid > water (TIP4P model) to an ideal gas, to find the relative free energy. To do > this, I believe one generally integrates above the critical point by > increasing the temperature above the critical temperature and then relaxing > the pressure until the system is a diffuse gas. The mdp options documentation > is helpful, and I went through an ethanol solvation tutorial, but there > doesn't appear to be a "pressure-lambda" or a "volume-lambda" option that I > could use to do the second part. How can I get Gromacs to calculate the dh/dl > derivative while relaxing the pressure? > > In addition, all of the tutorials I found for thermodynamic integration were > for finding solvation free energies. The coulomb and VDW forces are > essentially changed from "completely on" to "completely off". But in my case, > I'd like to change the temperature and pressure between two nonzero values. I > don't want to begin my simulation at 0K and 0atm, but lambda *must* go from 0 > to 1. How can I define both starting and ending points for the temperature > and pressure (or volume, or density)? > > Thanks for your help! > Nathan H. > -- > Gromacs Users mailing list > > * Please search the archive at > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a > mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Free Energy of Liquid Water
Thank you, Professor Farais de Moura. That method of integrating to an ideal gas makes sense. However, I'm having trouble figuring out how to deal with the temperature. I hope my question isn't too basic, but I can't find any examples online and I know I am misunderstanding how it works. I want to integrate from about 200K to 1000K or so (above the critical temperature, which I believe is somewhere in the 700's of Kelvin for TIP4P water). To do this, I thought I would set the reference temperature to 1000, and increment temperature-lambda from 0.2 to 1, thinking that temperature-lambda would scale the absolute temperature. But after running a the simulations, it's evident that temperature-lambda does not affect the thermostat. Should I set the reference temperature to different temperatures for each run? What does the temperature-lambda affect in that case? Thanks for your help! Regards, Nathan H. - Original Message - From: "André Farias de Moura" To: "Discussion list for GROMACS users" Sent: Friday, October 2, 2015 9:26:28 AM Subject: Re: [gmx-users] Free Energy of Liquid Water Apart from stability/convergence issues, I guess that turning off all intermolecular interactions should take you to the ideal gas straightforwardly, but in a different (P,T) point as compared to your target. But if you managed to alchemically turn water into an ideal gas, then you just need to apply standard free change for an ideal gas along a (P,T) process to achieve your target state. (I have not found the reference, but I read a paper doing just that with Monte Carlo simulations a few years ago) you should be able to track the conversion of water into an ideal gas by means of the g(r) profiles, which should change from the typical TIP4P profiles to g(r)=1 for all distances ranging from zero to half of the smallest box length. On Thu, Oct 1, 2015 at 11:44 PM, Nathan K Houtz wrote: > Hi everyone, > > I would like to use Gromacs to do Thermodynamic Integration (TI) from > liquid water (TIP4P model) to an ideal gas, to find the relative free > energy. To do this, I believe one generally integrates above the critical > point by increasing the temperature above the critical temperature and then > relaxing the pressure until the system is a diffuse gas. The mdp options > documentation is helpful, and I went through an ethanol solvation tutorial, > but there doesn't appear to be a "pressure-lambda" or a "volume-lambda" > option that I could use to do the second part. How can I get Gromacs to > calculate the dh/dl derivative while relaxing the pressure? > > In addition, all of the tutorials I found for thermodynamic integration > were for finding solvation free energies. The coulomb and VDW forces are > essentially changed from "completely on" to "completely off". But in my > case, I'd like to change the temperature and pressure between two nonzero > values. I don't want to begin my simulation at 0K and 0atm, but lambda > *must* go from 0 to 1. How can I define both starting and ending points for > the temperature and pressure (or volume, or density)? > > Thanks for your help! > Nathan H. > -- > Gromacs Users mailing list > > * Please search the archive at > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before > posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > send a mail to gmx-users-requ...@gromacs.org. > -- _ Prof. Dr. André Farias de Moura Department of Chemistry Federal University of São Carlos São Carlos - Brazil phone: +55-16-3351-8090 -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Free Energy of Liquid Water
Hi everyone, I would like to use Gromacs to do Thermodynamic Integration (TI) from liquid water (TIP4P model) to an ideal gas, to find the relative free energy. To do this, I believe one generally integrates above the critical point by increasing the temperature above the critical temperature and then relaxing the pressure until the system is a diffuse gas. The mdp options documentation is helpful, and I went through an ethanol solvation tutorial, but there doesn't appear to be a "pressure-lambda" or a "volume-lambda" option that I could use to do the second part. How can I get Gromacs to calculate the dh/dl derivative while relaxing the pressure? In addition, all of the tutorials I found for thermodynamic integration were for finding solvation free energies. The coulomb and VDW forces are essentially changed from "completely on" to "completely off". But in my case, I'd like to change the temperature and pressure between two nonzero values. I don't want to begin my simulation at 0K and 0atm, but lambda *must* go from 0 to 1. How can I define both starting and ending points for the temperature and pressure (or volume, or density)? Thanks for your help! Nathan H. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] chemical potential differences between two component liquids and crystal structures in gromacs
Hello everyone, I am interested in calculating the change in chemical potential between a metastable liquid, consisting or solvent and solute, and the crystal structure of the pure solute. For instance, one system I am looking at is tetrolic acid in CCL4 as my metastable liquid, and one of the polymorphs of tetrolic acid as the final crystal structure. I have looked through a couple tutorials on free energy calculations using gromacs, but I get the impression that this type of calculation may be less routine than say a solvation free energy which I came across. Can anyone tell me if this type of calculation is possible in gromacs, and where else I could be looking for directions as to how to proceed? Thank you, Nathan -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] simulations of ice beginning to spin
Thanks for your reply. How do I get gromacs to recognize that I want a triclinic box? I thought it was simply defined in the .gro file (as VMD recognizes). I know I can use gmx editconf to modify the gro file to define a triclinic box, but do I need anything in my input file? Thanks, Nathan - Original Message - From: "Vitaly V. Chaban" To: "gmx-users" Sent: Sunday, August 23, 2015 10:36:12 AM Subject: Re: [gmx-users] simulations of ice beginning to spin I suggest you to make a movie... Maybe, indeed, something is wrong with triclinic box definition (?) This not a rotation problem, since your crystal is outside the box at the end. Thus, the box which VMD drew is not the box which gromacs used during dynamics. On Fri, Aug 21, 2015 at 9:13 PM, Nathan K Houtz wrote: > Hello Eric, > > Thanks for your reply. I actually had made a mistake there, so that was a > good catch. I used "Linear" instead of "Angular" to constrain movement, so > gromacs actually wasn't fixing the rotation. However, I fixed it and redid > the simulations and I'm still getting very nearly the same results. Does > anything else come to mind? > > Thanks for your help, > Nathan Houtz > > - Original Message - > From: "Eric Smoll" > To: gmx-us...@gromacs.org > Sent: Thursday, August 20, 2015 5:53:18 PM > Subject: Re: [gmx-users] simulations of ice beginning to spin > > Hello Nathan, > > Did you set comm-mode and nstcomm > <http://manual.gromacs.org/online/mdp_opt.html#run> in your mdp file? > > Best, > Eric > > On Thu, Aug 20, 2015 at 3:37 PM, Nathan K Houtz wrote: > >> Hello, >> >> I'm sorry to submit a query twice, but I have not received a reply since >> last week. The original message is below. As the subject line states, my >> simulations are developing a rotation that I cannot explain. I was >> wondering if the box size could be an issue, so I redid both simulations in >> NPT with the Parrinello-Rahman Barostat. It did not solve the rotation >> problem. Both outputs looked very similar to the results from NVT. In >> addition, the corners of the cubic simulation still became disordered. >> That's the only update I have. Please see below for further details. Thanks! >> >> Regards, >> Nathan >> >> - Original Message - >> From: "Nathan K Houtz" >> To: "gromacs org gmx-users" >> Sent: Thursday, August 13, 2015 10:57:43 PM >> Subject: simulations of ice beginning to spin >> >> Hello, >> >> I am simulating two kinds of ice: ice Ih and ice Ic (cubic ice). Both >> simulations seem to have developed some rotation spontaneously and I'd like >> to know if I can control that somehow. I'm also simulating one in a >> triclinic box but the output gro file appears cubic. I'm not sure if I >> should be concerned about that or not. Lastly, the Ice Ic structure begins >> to fall apart during my simulation and I don't think it should. I'd >> appreciate help on any of these matters. Here are some more details: >> >> Ice Ih is done in a triclinic box, and originally looks like this: >> http://imgur.com/VQg9xQz. I use a rigid TIP4P/Ice water model, >> constrained by shake, and simulate it in NVT for 1,000,000 time steps (2 >> ns) at a temperature of 217K and a density of 0.920 g/cm3. At the end, it >> looks like this: http://imgur.com/8uzddRH. First of all, I'm wondering if >> gromacs has correctly applied the triclinic box, as it appears that >> periodic boundary conditions have turned it into a cube. VMD is showing the >> box described by the .gro file. Secondly, it has been rotated, at least a >> few degrees, about the z axis. It did not noticeably rotate about the x or >> y axes. I'm not sure how to explain how the simulation would develop any >> angular momentum. In my .mdp file, I have gromacs get rid of angular >> momentum every 100 steps, but I don't think it should develop any in the >> first place. Any advice? >> >> The other simulation, ice ic, didn't go quite as well. Here is the >> original orientation, viewed from one of the corners: >> http://imgur.com/i0ncpju (I realize the orthographic projections make it >> harder to see the actual structure, but they make it easier to see unique >> patterns from various angles, which I'm trying to use to determine how it >> has rotated). This one was also simulated for 1,000,000 timesteps (2 ns) at >> 217K with TIP4P/Ice constrained by shake, but at a density of 0.931 g/cm3. >> The resulting structure is this: http://imgur.com/Lwl
Re: [gmx-users] simulations of ice beginning to spin
Hello Eric, Thanks for your reply. I actually had made a mistake there, so that was a good catch. I used "Linear" instead of "Angular" to constrain movement, so gromacs actually wasn't fixing the rotation. However, I fixed it and redid the simulations and I'm still getting very nearly the same results. Does anything else come to mind? Thanks for your help, Nathan Houtz - Original Message - From: "Eric Smoll" To: gmx-us...@gromacs.org Sent: Thursday, August 20, 2015 5:53:18 PM Subject: Re: [gmx-users] simulations of ice beginning to spin Hello Nathan, Did you set comm-mode and nstcomm <http://manual.gromacs.org/online/mdp_opt.html#run> in your mdp file? Best, Eric On Thu, Aug 20, 2015 at 3:37 PM, Nathan K Houtz wrote: > Hello, > > I'm sorry to submit a query twice, but I have not received a reply since > last week. The original message is below. As the subject line states, my > simulations are developing a rotation that I cannot explain. I was > wondering if the box size could be an issue, so I redid both simulations in > NPT with the Parrinello-Rahman Barostat. It did not solve the rotation > problem. Both outputs looked very similar to the results from NVT. In > addition, the corners of the cubic simulation still became disordered. > That's the only update I have. Please see below for further details. Thanks! > > Regards, > Nathan > > - Original Message - > From: "Nathan K Houtz" > To: "gromacs org gmx-users" > Sent: Thursday, August 13, 2015 10:57:43 PM > Subject: simulations of ice beginning to spin > > Hello, > > I am simulating two kinds of ice: ice Ih and ice Ic (cubic ice). Both > simulations seem to have developed some rotation spontaneously and I'd like > to know if I can control that somehow. I'm also simulating one in a > triclinic box but the output gro file appears cubic. I'm not sure if I > should be concerned about that or not. Lastly, the Ice Ic structure begins > to fall apart during my simulation and I don't think it should. I'd > appreciate help on any of these matters. Here are some more details: > > Ice Ih is done in a triclinic box, and originally looks like this: > http://imgur.com/VQg9xQz. I use a rigid TIP4P/Ice water model, > constrained by shake, and simulate it in NVT for 1,000,000 time steps (2 > ns) at a temperature of 217K and a density of 0.920 g/cm3. At the end, it > looks like this: http://imgur.com/8uzddRH. First of all, I'm wondering if > gromacs has correctly applied the triclinic box, as it appears that > periodic boundary conditions have turned it into a cube. VMD is showing the > box described by the .gro file. Secondly, it has been rotated, at least a > few degrees, about the z axis. It did not noticeably rotate about the x or > y axes. I'm not sure how to explain how the simulation would develop any > angular momentum. In my .mdp file, I have gromacs get rid of angular > momentum every 100 steps, but I don't think it should develop any in the > first place. Any advice? > > The other simulation, ice ic, didn't go quite as well. Here is the > original orientation, viewed from one of the corners: > http://imgur.com/i0ncpju (I realize the orthographic projections make it > harder to see the actual structure, but they make it easier to see unique > patterns from various angles, which I'm trying to use to determine how it > has rotated). This one was also simulated for 1,000,000 timesteps (2 ns) at > 217K with TIP4P/Ice constrained by shake, but at a density of 0.931 g/cm3. > The resulting structure is this: http://imgur.com/Lwl8UUw. I'm not sure > I've got the exact corresponding orientation in this view here, but I > believe it is. It's looking down the x-axis this time. This simulation got > rotated much worse. Here are the before and after views looking down the > z-axis, just like I showed for the ice Ih above: http://imgur.com/a/LRQR3. > With such a large number of timesteps, I cannot output enough frames in the > trajectory file to view this as a "smooth" movie and ide > ntify exactly when and why these rotations happen. It's more like a very > long and rapid slide show. For the ice Ic, at least, it is possible that > the whole box did not rotate, but all of the molecules simply reoriented > themselves. But since at least part of the box appears to still be in the > ice Ic structure, I'm not sure why they would do that. Ice Ic is only > metastable at any temperature, so I would assume that if it falls apart, > the structure would not come back. This is another issue. I would like to > still have cubic ice structure at the end of this simulation, but clearly > two of the corn
Re: [gmx-users] simulations of ice beginning to spin
Hello, I'm sorry to submit a query twice, but I have not received a reply since last week. The original message is below. As the subject line states, my simulations are developing a rotation that I cannot explain. I was wondering if the box size could be an issue, so I redid both simulations in NPT with the Parrinello-Rahman Barostat. It did not solve the rotation problem. Both outputs looked very similar to the results from NVT. In addition, the corners of the cubic simulation still became disordered. That's the only update I have. Please see below for further details. Thanks! Regards, Nathan - Original Message ----- From: "Nathan K Houtz" To: "gromacs org gmx-users" Sent: Thursday, August 13, 2015 10:57:43 PM Subject: simulations of ice beginning to spin Hello, I am simulating two kinds of ice: ice Ih and ice Ic (cubic ice). Both simulations seem to have developed some rotation spontaneously and I'd like to know if I can control that somehow. I'm also simulating one in a triclinic box but the output gro file appears cubic. I'm not sure if I should be concerned about that or not. Lastly, the Ice Ic structure begins to fall apart during my simulation and I don't think it should. I'd appreciate help on any of these matters. Here are some more details: Ice Ih is done in a triclinic box, and originally looks like this: http://imgur.com/VQg9xQz. I use a rigid TIP4P/Ice water model, constrained by shake, and simulate it in NVT for 1,000,000 time steps (2 ns) at a temperature of 217K and a density of 0.920 g/cm3. At the end, it looks like this: http://imgur.com/8uzddRH. First of all, I'm wondering if gromacs has correctly applied the triclinic box, as it appears that periodic boundary conditions have turned it into a cube. VMD is showing the box described by the .gro file. Secondly, it has been rotated, at least a few degrees, about the z axis. It did not noticeably rotate about the x or y axes. I'm not sure how to explain how the simulation would develop any angular momentum. In my .mdp file, I have gromacs get rid of angular momentum every 100 steps, but I don't think it should develop any in the first place. Any advice? The other simulation, ice ic, didn't go quite as well. Here is the original orientation, viewed from one of the corners: http://imgur.com/i0ncpju (I realize the orthographic projections make it harder to see the actual structure, but they make it easier to see unique patterns from various angles, which I'm trying to use to determine how it has rotated). This one was also simulated for 1,000,000 timesteps (2 ns) at 217K with TIP4P/Ice constrained by shake, but at a density of 0.931 g/cm3. The resulting structure is this: http://imgur.com/Lwl8UUw. I'm not sure I've got the exact corresponding orientation in this view here, but I believe it is. It's looking down the x-axis this time. This simulation got rotated much worse. Here are the before and after views looking down the z-axis, just like I showed for the ice Ih above: http://imgur.com/a/LRQR3. With such a large number of timesteps, I cannot output enough frames in the trajectory file to view this as a "smooth" movie and ide ntify exactly when and why these rotations happen. It's more like a very long and rapid slide show. For the ice Ic, at least, it is possible that the whole box did not rotate, but all of the molecules simply reoriented themselves. But since at least part of the box appears to still be in the ice Ic structure, I'm not sure why they would do that. Ice Ic is only metastable at any temperature, so I would assume that if it falls apart, the structure would not come back. This is another issue. I would like to still have cubic ice structure at the end of this simulation, but clearly two of the corners have disastrously broken up. However, what looks like most of the molecules appear to still retain the structure. Do you think that it could be a problem with my box size, or the orientation? I'd like to know what things I might try to prevent the disintegration of the cubic structure. Thanks for any help. Regards, Nathan -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] simulations of ice beginning to spin
Hello, I am simulating two kinds of ice: ice Ih and ice Ic (cubic ice). Both simulations seem to have developed some rotation spontaneously and I'd like to know if I can control that somehow. I'm also simulating one in a triclinic box but the output gro file appears cubic. I'm not sure if I should be concerned about that or not. Lastly, the Ice Ic structure begins to fall apart during my simulation and I don't think it should. I'd appreciate help on any of these matters. Here are some more details: Ice Ih is done in a triclinic box, and originally looks like this: http://imgur.com/VQg9xQz. I use a rigid TIP4P/Ice water model, constrained by shake, and simulate it in NVT for 1,000,000 time steps (2 ns) at a temperature of 217K and a density of 0.920 g/cm3. At the end, it looks like this: http://imgur.com/8uzddRH. First of all, I'm wondering if gromacs has correctly applied the triclinic box, as it appears that periodic boundary conditions have turned it into a cube. VMD is showing the box described by the .gro file. Secondly, it has been rotated, at least a few degrees, about the z axis. It did not noticeably rotate about the x or y axes. I'm not sure how to explain how the simulation would develop any angular momentum. In my .mdp file, I have gromacs get rid of angular momentum every 100 steps, but I don't think it should develop any in the first place. Any advice? The other simulation, ice ic, didn't go quite as well. Here is the original orientation, viewed from one of the corners: http://imgur.com/i0ncpju (I realize the orthographic projections make it harder to see the actual structure, but they make it easier to see unique patterns from various angles, which I'm trying to use to determine how it has rotated). This one was also simulated for 1,000,000 timesteps (2 ns) at 217K with TIP4P/Ice constrained by shake, but at a density of 0.931 g/cm3. The resulting structure is this: http://imgur.com/Lwl8UUw. I'm not sure I've got the exact corresponding orientation in this view here, but I believe it is. It's looking down the x-axis this time. This simulation got rotated much worse. Here are the before and after views looking down the z-axis, just like I showed for the ice Ih above: http://imgur.com/a/LRQR3. With such a large number of timesteps, I cannot output enough frames in the trajectory file to view this as a "smooth" movie and ide ntify exactly when and why these rotations happen. It's more like a very long and rapid slide show. For the ice Ic, at least, it is possible that the whole box did not rotate, but all of the molecules simply reoriented themselves. But since at least part of the box appears to still be in the ice Ic structure, I'm not sure why they would do that. Ice Ic is only metastable at any temperature, so I would assume that if it falls apart, the structure would not come back. This is another issue. I would like to still have cubic ice structure at the end of this simulation, but clearly two of the corners have disastrously broken up. However, what looks like most of the molecules appear to still retain the structure. Do you think that it could be a problem with my box size, or the orientation? I'd like to know what things I might try to prevent the disintegration of the cubic structure. Thanks for any help. Regards, Nathan -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] trouble diagnosing a simulation that's blowing up (shake not converging --> segmentation fault)
Hi, Thanks Mark. I'm sorry though, I don't think I understand what you mean. I thought the [system] section of the topology file was just a name for the system. How should I imply an order for the atoms? I did double check that the ordering in atoms is the same as it is in my .gro file (grommp -c input). 'OW', 'HW', 'HW', then 'DW' - Original Message - From: "Mark Abraham" To: gmx-us...@gromacs.org Sent: Friday, July 10, 2015 7:24:49 PM Subject: Re: [gmx-users] trouble diagnosing a simulation that's blowing up (shake not converging --> segmentation fault) Hi, That could only happen if your grompp -c input used a different atom ordering from that implied by your [system] section of your .top. grompp warns about this, but maybe you didn't notice... (Or the structure you loaded into VMD is a mismatch to the trajectory, so its heuristics for guessing where bonds actually are get double-crossed.) Mark On Sat, Jul 11, 2015 at 1:15 AM Nathan K Houtz wrote: > Actually, I think I found the problem. When i looked at the pdb files in > vmd the first time I missed it, but a colleague had a hunch and I found > that he was right! I think that shake is constraining atoms to the wrong > molecules. Here's a screenshot I took where you can see what I mean: > http://imgur.com/1bgThDv > > I'm still not sure how to fix it though. The molecule numbering in my .gro > file is correct, and there are only 4 atoms per molecule (or 3 plus the > dummy). How can I tell shake to look at each molecule individually? Or have > I made a mistake somewhere? Here are my most recent files: > > .gro file: (only part of it -- there is a character limit. this is the > first and last ten molecules with the box dimensions at the bottom) > http://textuploader.com/e898 > > topology file: > http://textuploader.com/e89g > > minimization: > http://textuploader.com/e8rc > > .mdp file: > http://textuploader.com/e896 > > log file from a failed simulation: > http://textuploader.com/e80y > > Thanks for your help! If you do still want the full .gro file to run the > simulation, I can look for a different file host. > Nathan > > > - Original Message - > From: "Justin Lemkul" > To: gmx-us...@gromacs.org > Sent: Friday, July 10, 2015 5:06:34 PM > Subject: Re: [gmx-users] trouble diagnosing a simulation that's blowing up > (shake not converging --> segmentation fault) > > > > On 7/9/15 9:23 PM, Nathan K Houtz wrote: > > Thanks for your explanations, Dr. Lemkul. > > > > I had already corrected a couple of the things you suggested. Gromacs > won't actually let me run with Nose-Hoover and Parrinello-Rahman together > (or at least, it gives a warning not to do that and stops). I'd like to run > in NVT anyway, so I had set Pcoupl to 'no'. I had also increased the > cutoffs and my nstlist is 20. > > > > I have now also changed the tcoupl to v-rescale, but unfortunately that > alone didn't help. So I'm not sure exactly what I can inspect to find the > source of my error. I know that temperature, pressure, and total energy are > all warning signs of bad things if they misbehave, but they were all fairly > constant throughout the simulation. (This is for the flexible-model > simulation:) The starting energy for 1700 water molecules at 120K was > -1.08744e+05, and it finished at -1.06047e+05 with no significant > variations on any step. Pressure and temperature were also fine. I > attempted to look at the results visually in vmd but I might have done > something wrong because the .trr file has some error in it and vmd crashes. > But the starting geometry after minimization looks fine: none of the > molecules were moved out of their position in the crystal. > > > > For the constrained molecules, the energy stays about the same (about > -1e+5) for the first 8 timesteps and then quickly blows up to +2e15 at step > 12 before it obviously fails. I get shake warnings from step 0 though. In > vmd, the first 8 steps look reasonable (checking the .pdb files that > gromacs outputs when there are warnings) and all the molecules seem to hold > their places in the crystal, but then at step 9 suddenly some of the > molecules become misshapen with very long or very short bonds. Of course it > goes downhill from there, but I can't figure out what's causing the > problems since it's clearly happening from the very beginning (according to > the shake warnings). > > > > What other things could I look at to troubleshoot my problem? > > > > Not sure, but at least the failure happens fast. Upload all your input > files > somewhere and I'll
Re: [gmx-users] trouble diagnosing a simulation that's blowing up (shake not converging --> segmentation fault)
Actually, I think I found the problem. When i looked at the pdb files in vmd the first time I missed it, but a colleague had a hunch and I found that he was right! I think that shake is constraining atoms to the wrong molecules. Here's a screenshot I took where you can see what I mean:http://imgur.com/1bgThDv I'm still not sure how to fix it though. The molecule numbering in my .gro file is correct, and there are only 4 atoms per molecule (or 3 plus the dummy). How can I tell shake to look at each molecule individually? Or have I made a mistake somewhere? Here are my most recent files: .gro file: (only part of it -- there is a character limit. this is the first and last ten molecules with the box dimensions at the bottom) http://textuploader.com/e898 topology file: http://textuploader.com/e89g minimization: http://textuploader.com/e8rc .mdp file: http://textuploader.com/e896 log file from a failed simulation: http://textuploader.com/e80y Thanks for your help! If you do still want the full .gro file to run the simulation, I can look for a different file host. Nathan - Original Message - From: "Justin Lemkul" To: gmx-us...@gromacs.org Sent: Friday, July 10, 2015 5:06:34 PM Subject: Re: [gmx-users] trouble diagnosing a simulation that's blowing up (shake not converging --> segmentation fault) On 7/9/15 9:23 PM, Nathan K Houtz wrote: > Thanks for your explanations, Dr. Lemkul. > > I had already corrected a couple of the things you suggested. Gromacs won't > actually let me run with Nose-Hoover and Parrinello-Rahman together (or at > least, it gives a warning not to do that and stops). I'd like to run in NVT > anyway, so I had set Pcoupl to 'no'. I had also increased the cutoffs and my > nstlist is 20. > > I have now also changed the tcoupl to v-rescale, but unfortunately that alone > didn't help. So I'm not sure exactly what I can inspect to find the source of > my error. I know that temperature, pressure, and total energy are all warning > signs of bad things if they misbehave, but they were all fairly constant > throughout the simulation. (This is for the flexible-model simulation:) The > starting energy for 1700 water molecules at 120K was -1.08744e+05, and it > finished at -1.06047e+05 with no significant variations on any step. Pressure > and temperature were also fine. I attempted to look at the results visually > in vmd but I might have done something wrong because the .trr file has some > error in it and vmd crashes. But the starting geometry after minimization > looks fine: none of the molecules were moved out of their position in the > crystal. > > For the constrained molecules, the energy stays about the same (about -1e+5) > for the first 8 timesteps and then quickly blows up to +2e15 at step 12 > before it obviously fails. I get shake warnings from step 0 though. In vmd, > the first 8 steps look reasonable (checking the .pdb files that gromacs > outputs when there are warnings) and all the molecules seem to hold their > places in the crystal, but then at step 9 suddenly some of the molecules > become misshapen with very long or very short bonds. Of course it goes > downhill from there, but I can't figure out what's causing the problems since > it's clearly happening from the very beginning (according to the shake > warnings). > > What other things could I look at to troubleshoot my problem? > Not sure, but at least the failure happens fast. Upload all your input files somewhere and I'll take a few minutes to see if I can spot anything. -Justin > Thanks for your help, > Nathan > > - Original Message - > From: "Justin Lemkul" > To: gmx-us...@gromacs.org > Sent: Thursday, July 9, 2015 7:39:47 AM > Subject: Re: [gmx-users] trouble diagnosing a simulation that's blowing up > (shake not converging --> segmentation fault) > > > > On 7/8/15 8:06 PM, Nathan K Houtz wrote: >> Hello, >> >> I deleted the email and can't respond to my last reply directly - sorry! I >> got this response from Mark Abraham: >> >> Hi, Try doing some EM and initial equilibration with no constraints at all, >> perhaps? Mark >> >> I tried commenting out the shake commands, and got a short (5000 step) >> simulation to run just fine without blowing up. Before, I would get shake >> warnings from the first few steps and a segmentation fault around step 13 or >> 14. I would like to be able to simulate with rigid molecules, though. Why >> would the simulation work with flexible molecules but not rigid ones? >> > > Flexible water allows weird geometry, which is probably coming up and causing > your constraint algorithm to fail. I'd insp
Re: [gmx-users] trouble diagnosing a simulation that's blowing up (shake not converging --> segmentation fault)
Thanks for your explanations, Dr. Lemkul. I had already corrected a couple of the things you suggested. Gromacs won't actually let me run with Nose-Hoover and Parrinello-Rahman together (or at least, it gives a warning not to do that and stops). I'd like to run in NVT anyway, so I had set Pcoupl to 'no'. I had also increased the cutoffs and my nstlist is 20. I have now also changed the tcoupl to v-rescale, but unfortunately that alone didn't help. So I'm not sure exactly what I can inspect to find the source of my error. I know that temperature, pressure, and total energy are all warning signs of bad things if they misbehave, but they were all fairly constant throughout the simulation. (This is for the flexible-model simulation:) The starting energy for 1700 water molecules at 120K was -1.08744e+05, and it finished at -1.06047e+05 with no significant variations on any step. Pressure and temperature were also fine. I attempted to look at the results visually in vmd but I might have done something wrong because the .trr file has some error in it and vmd crashes. But the starting geometry after minimization looks fine: none of the molecules were moved out of their position in the crystal. For the constrained molecules, the energy stays about the same (about -1e+5) for the first 8 timesteps and then quickly blows up to +2e15 at step 12 before it obviously fails. I get shake warnings from step 0 though. In vmd, the first 8 steps look reasonable (checking the .pdb files that gromacs outputs when there are warnings) and all the molecules seem to hold their places in the crystal, but then at step 9 suddenly some of the molecules become misshapen with very long or very short bonds. Of course it goes downhill from there, but I can't figure out what's causing the problems since it's clearly happening from the very beginning (according to the shake warnings). What other things could I look at to troubleshoot my problem? Thanks for your help, Nathan - Original Message - From: "Justin Lemkul" To: gmx-us...@gromacs.org Sent: Thursday, July 9, 2015 7:39:47 AM Subject: Re: [gmx-users] trouble diagnosing a simulation that's blowing up (shake not converging --> segmentation fault) On 7/8/15 8:06 PM, Nathan K Houtz wrote: > Hello, > > I deleted the email and can't respond to my last reply directly - sorry! I > got this response from Mark Abraham: > > Hi, Try doing some EM and initial equilibration with no constraints at all, > perhaps? Mark > > I tried commenting out the shake commands, and got a short (5000 step) > simulation to run just fine without blowing up. Before, I would get shake > warnings from the first few steps and a segmentation fault around step 13 or > 14. I would like to be able to simulate with rigid molecules, though. Why > would the simulation work with flexible molecules but not rigid ones? > Flexible water allows weird geometry, which is probably coming up and causing your constraint algorithm to fail. I'd inspect the outcome carefully. Just because it runs doesn't mean it's right. > Also, in the example .mdp file for tip4p water, there is the (outdated) > option, 'unconstrained-start', which is now 'continuation'. I got errors when > trying to make the input .tpr file when I attempted to set that option to > 'yes'. The warning said it was because I want Gromacs to generate velocities > to start the simulation, which is incompatible with that command. Is there > another way I can try to start the simulation unconstrained? Or would you > suggest another idea to fix my shake warnings? > What this setting says is "have the constraints already been solved (continuation = yes) or should mdrun constrain the starting configuration itself (continuation = no)." The .mdp file has a number of weird settings. I would never use Nose-Hoover and Parrinello-Rahman when generating velocities; that's likely to be very unstable. See if a more forgiving thermostat and/or barostat resolves the issue, e.g. tcoupl = v-rescale and pcoupl = Berendsen. Also note that the cutoffs are very short (probably because they want to use a small box, but note that it *does* affect the physics) and nstlist = 1 is totally unnecessary. It doesn't hurt your physics, but it's a major waste of performance. -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Pl
Re: [gmx-users] trouble diagnosing a simulation that's blowing up (shake not converging --> segmentation fault)
Hello, I deleted the email and can't respond to my last reply directly - sorry! I got this response from Mark Abraham: Hi, Try doing some EM and initial equilibration with no constraints at all, perhaps? Mark I tried commenting out the shake commands, and got a short (5000 step) simulation to run just fine without blowing up. Before, I would get shake warnings from the first few steps and a segmentation fault around step 13 or 14. I would like to be able to simulate with rigid molecules, though. Why would the simulation work with flexible molecules but not rigid ones? Also, in the example .mdp file for tip4p water, there is the (outdated) option, 'unconstrained-start', which is now 'continuation'. I got errors when trying to make the input .tpr file when I attempted to set that option to 'yes'. The warning said it was because I want Gromacs to generate velocities to start the simulation, which is incompatible with that command. Is there another way I can try to start the simulation unconstrained? Or would you suggest another idea to fix my shake warnings? Thanks very much for your help, Nathan - Original Message - From: "Nathan K Houtz" To: gmx-us...@gromacs.org Sent: Monday, July 6, 2015 9:08:33 PM Subject: trouble diagnosing a simulation that's blowing up (shake not converging --> segmentation fault) Hello, I'm attempting to simulate ice ih in a triclinic box with a rigid tip4p water model. As the subject says, I get warnings (almost immediately) that "shake failed to converge in 1000 steps" and then eventually a segmentation fault. Gromacs documentation suggests that this is a result of my simulation "blowing up." I tried decreasing the timestep down to 0.0001 ps, lowering the temperature, and minimizing further (down to 100 kJ/mol/nm) but nothing worked. I mostly used modified versions of the files available on this page: http://www.sklogwiki.org/SklogWiki/index.php/GROMACS_files_for_the_TIP4P/2005_model except the .gro file, which is my own. I did comment out the pressure coupling algorithms because I got a note that the barostat was unsuitable for equillibration, but I want NVT anyway. The temperature is 120.0 K and the coulomb type is PME. I've double checked and am satisfied with the interaction parameters in the topology file, so I'm out of ideas on where my mistake migh! t be. Here are the commands and files I used, in case it helps: > gmx grompp -f minim.mdp -c water_5100.gro -p water_topol.top -o em.tpr > gmx mdrun -v -deffnm em > gmx grompp -f water_5100.mdp -c em.gro -p water_topol.top -o water.tpr > gmx mdrun -deffnm water Here is the minim.mdp file: http://textuploader.com/iggv The water_topol.top file: http://textuploader.com/igpd And the water_5100.mdp file: http://textuploader.com/igzj You'll notice a slight mistake in the naming convention. There are 1700 molecules in the simulation, but everything is named something like '..._5100', even though tip4p water has a virtual site (making 4 particles). My .gro file does in fact have 6800 atoms, not 5100, I just neglected to rename the files because it's easier. Unfortunately I can't upload the .gro file because it's too big. I realize the problem could be there but I won't attach a large file unless it's specifically requested. Here's a sample: Triclinic Frozen Water Simulation 6800 1WATER OW1 0.261 0.000 0.000 0. 0. 0. 1WATER HW2 0.261 0.000 0.082 0. 0. 0. 1WATER HW3 0.343 -0.000 -0.027 0. 0. 0. 1WATER DW4 0.272 0.000 0.007 0. 0. 0. 2WATER OW5 0.261 0.000 0.736 0. 0. 0. 2WATER HW6 0.261 0.000 0.818 0. 0. 0. 2WATER HW7 0.343 0.000 0.709 0. 0. 0. 2WATER DW8 0.272 0.000 0.743 0. 0. 0. 3WATER OW9 0.261 0.000 1.472 0. 0. 0. 3WATER HW 10 0.261 0.000 1.554 0. 0. 0. 3WATER HW 11 0.343 0.000 1.445 0. 0. 0. 3WATER DW 12 0.272 0.000 1.479 0. 0. 0. 4WATER OW 13 0.261 0.000 2.208 0. 0. 0. 4WATER HW 14 0.261 0.000 2.290 0. 0. 0. 4WATER HW 15 0.343 0.000 2.181 0. 0. 0. 4WATER DW 16 0.272 0.000 2.215 0. 0. 0. 5WATER OW 17 0.261 0.000 2.944 0. 0. 0. 5WATER HW 18 0.261 0.000 3.026 0. 0. 0. 5WATER HW 19 0.343 0.000 2.917 0. 0. 0. 5WATER DW 20 0.272 0.000 2.951 0. 0. 0. 6WATER OW 21 -0.130 0.677 0.000 0. 0. 0. 6WATER HW 22 -0.130 0.677
[gmx-users] trouble diagnosing a simulation that's blowing up (shake not converging --> segmentation fault)
Hello, I'm attempting to simulate ice ih in a triclinic box with a rigid tip4p water model. As the subject says, I get warnings (almost immediately) that "shake failed to converge in 1000 steps" and then eventually a segmentation fault. Gromacs documentation suggests that this is a result of my simulation "blowing up." I tried decreasing the timestep down to 0.0001 ps, lowering the temperature, and minimizing further (down to 100 kJ/mol/nm) but nothing worked. I mostly used modified versions of the files available on this page: http://www.sklogwiki.org/SklogWiki/index.php/GROMACS_files_for_the_TIP4P/2005_model except the .gro file, which is my own. I did comment out the pressure coupling algorithms because I got a note that the barostat was unsuitable for equillibration, but I want NVT anyway. The temperature is 120.0 K and the coulomb type is PME. I've double checked and am satisfied with the interaction parameters in the topology file, so I'm out of ideas on where my mistake migh! t be. Here are the commands and files I used, in case it helps: > gmx grompp -f minim.mdp -c water_5100.gro -p water_topol.top -o em.tpr > gmx mdrun -v -deffnm em > gmx grompp -f water_5100.mdp -c em.gro -p water_topol.top -o water.tpr > gmx mdrun -deffnm water Here is the minim.mdp file: http://textuploader.com/iggv The water_topol.top file: http://textuploader.com/igpd And the water_5100.mdp file: http://textuploader.com/igzj You'll notice a slight mistake in the naming convention. There are 1700 molecules in the simulation, but everything is named something like '..._5100', even though tip4p water has a virtual site (making 4 particles). My .gro file does in fact have 6800 atoms, not 5100, I just neglected to rename the files because it's easier. Unfortunately I can't upload the .gro file because it's too big. I realize the problem could be there but I won't attach a large file unless it's specifically requested. Here's a sample: Triclinic Frozen Water Simulation 6800 1WATER OW1 0.261 0.000 0.000 0. 0. 0. 1WATER HW2 0.261 0.000 0.082 0. 0. 0. 1WATER HW3 0.343 -0.000 -0.027 0. 0. 0. 1WATER DW4 0.272 0.000 0.007 0. 0. 0. 2WATER OW5 0.261 0.000 0.736 0. 0. 0. 2WATER HW6 0.261 0.000 0.818 0. 0. 0. 2WATER HW7 0.343 0.000 0.709 0. 0. 0. 2WATER DW8 0.272 0.000 0.743 0. 0. 0. 3WATER OW9 0.261 0.000 1.472 0. 0. 0. 3WATER HW 10 0.261 0.000 1.554 0. 0. 0. 3WATER HW 11 0.343 0.000 1.445 0. 0. 0. 3WATER DW 12 0.272 0.000 1.479 0. 0. 0. 4WATER OW 13 0.261 0.000 2.208 0. 0. 0. 4WATER HW 14 0.261 0.000 2.290 0. 0. 0. 4WATER HW 15 0.343 0.000 2.181 0. 0. 0. 4WATER DW 16 0.272 0.000 2.215 0. 0. 0. 5WATER OW 17 0.261 0.000 2.944 0. 0. 0. 5WATER HW 18 0.261 0.000 3.026 0. 0. 0. 5WATER HW 19 0.343 0.000 2.917 0. 0. 0. 5WATER DW 20 0.272 0.000 2.951 0. 0. 0. 6WATER OW 21 -0.130 0.677 0.000 0. 0. 0. 6WATER HW 22 -0.130 0.677 0.082 0. 0. 0. 6WATER HW 23 -0.048 0.677 -0.027 0. 0. 0. 6WATER DW 24 -0.119 0.677 0.007 0. 0. 0. . . . . . . 1696WATER OW 6781 1.825 3.160 0.276 0. 0. 0. 1696WATER HW 6782 1.784 3.232 0.304 0. 0. 0. 1696WATER HW 6783 1.907 3.161 0.304 0. 0. 0. 1696WATER DW 6784 1.830 3.170 0.283 0. 0. 0. 1697WATER OW 6785 1.825 3.160 1.012 0. 0. 0. 1697WATER HW 6786 1.784 3.232 1.040 0. 0. 0. 1697WATER HW 6787 1.907 3.161 1.040 0. 0. 0. 1697WATER DW 6788 1.830 3.170 1.019 0. 0. 0. 1698WATER OW 6789 1.825 3.160 1.748 0. 0. 0. 1698WATER HW 6790 1.784 3.232 1.776 0. 0. 0. 1698WATER HW 6791 1.907 3.161 1.776 0. 0. 0. 1698WATER DW 6792 1.830 3.170 1.755 0. 0. 0. 1699WATER OW 6793 1.825 3.160 2.484 0. 0. 0. 1699WATER HW 6794 1.784 3.232 2.512 0. 0. 0. 1699WATER HW 6795 1.907 3.161 2.512 0. 0. 0. 1699WATER DW 6796 1.830 3.170 2.491 0. 0. 0. 1700WATER OW 6797 1.825 3.160 3.220 0. 0. 0. 1700WATER HW 6798 1.784 3.232 3.248 0. 0. 0. 1700WATER HW 6799 1.907 3.161 3.248 0. 0. 0. 1700W
[gmx-users] Density of pure tetrolic acid does not agree with experimental data
Hello Gromacs users, I tried to create a 6.5 nm cubic box of tetrolic acid (otherwise known as 2-butynoic acid), but the density is too low. According to this, http://www.chemspider.com/Chemical-Structure.61810.html, the density should be about 0.964 g/mL. Given the molecular weight of tetrolic acid, this corresponds to right around 1900 molecules in the 6.5 nm box I prepared. Unfortunately, when I attempted to fill the box, Gromacs only found room for 1676 molecules, resulting in a density of 0.852 g/mL. I know that models are not 100% accurate, but I'm worried about a difference that big. Should I expect the system to condense when I minimize, and just shrink the box afterwards? Or do you think there is another problem? Thanks for your help, Nathan -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Electrostatic force cutoffs
Hello all, I'm wondering if gromacs will allow me to use either a wolf electrostatic cutoff method or a damp shifted force cutoff method, similar to what is described by C. J. Fennell (2006). (I think the paper is: http://scitation.aip.org/content/aip/journal/jcp/124/23/10.1063/1.2206581) If so, are there any tutorials or how-to's on how to implement them? If it matters, I am attempting to simulate tetrolic acid (also called 2-butynoic acid) in a solution of water (Tip3p). Thanks very much for your help! N.H. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Electrostatic force cutoffs
Hello all, I'm wondering if gromacs will allow me to use either a wolf electrostatic cutoff method or a damp shifted force cutoff method, similar to what is described by C. J. Fennell (2006). (I think the paper is: http://scitation.aip.org/content/aip/journal/jcp/124/23/10.1063/1.2206581) If so, are there any tutorials or how-to's on how to implement them? If it matters, I am attempting to simulate tetrolic acid (also called 2-butynoic acid) in a solution of water (Tip3p). Thanks very much for your help! N.H. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Simulating Multiple Solute Particles
Dr. Lemkul, Thank you again - your tutorials are very helpful for me. I want to get through your Virtual Sites tutorial but I got stuck near the beginning. Gromacs thinks the minim.mdp file is for an earlier version but the introduction page said it's for 5.0 and later. I'm using gromacs 5.0.2, so I thought I should be fine. I copied this command: > gmx grompp -f minim.mdp -c box.pdb -p topol.top -o min.tpr after having made my box, of course, and got the following error: >Command line: > gmx grompp -f minim.mdp -c box.pdb -p topol.top -o min.tpr > > >NOTE 1 [file minim.mdp, line 19]: > minim.mdp did not specify a value for the .mdp option "cutoff-scheme". > Probably it was first intended for use with GROMACS before 4.6. In 4.6, > the Verlet scheme was introduced, but the group scheme was still the > default. The default is now the Verlet scheme, so you will observe > different behaviour. > >Ignoring obsolete mdp entry 'title' > >ERROR 1 [file minim.mdp]: > With Verlet lists only full pbc or pbc=xy with walls is supported > > >ERROR 2 [file minim.mdp]: > With Verlet lists nstlist should be larger than 0 > > >NOTE 2 [file minim.mdp]: > With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note > that with the Verlet scheme, nstlist has no effect on the accuracy of > your simulation. I tried to figure out how I could modify it myself but I don't really know what I'm doing. The sample mdp file (http://manual.gromacs.org/online/mdp.html) is much longer than the one you linked (http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/vsites/Files/minim.mdp) and I'm not sure which fields are necessary. Is the file outdated, or have I made a mistake? Regards, N.H. - Original Message - From: "Justin Lemkul" To: gmx-us...@gromacs.org Sent: Monday, November 10, 2014 7:40:58 AM Subject: Re: [gmx-users] Simulating Multiple Solute Particles On 11/10/14 3:07 AM, Nathan K Houtz wrote: > Thank you again, Dr. Lemkul! I went with the first option because my advisor > wants me to use tip3p water. Unfortunately, I immediately got myself stuck on > the next step, but this time I think the solution may not be as trivial. The > command to output em.tpr seems to work without any errors. The command, > Well, strictly speaking, either option works because the parameters of TIP3P are not force field-dependent, just translated between different conventions for nonbonded parameters, but that's just a pedantic issue. >> gmx mdrun -deffnm em > > gives me a lot of errors. (It stops after 1000). Here's a small sample: > >> Steepest Descents: >>Tolerance (Fmax) = 1.0e+03 >>Number of steps=5 >> >> Step 1, time 0.001 (ps) LINCS WARNING >> relative constraint deviation after LINCS: >> rms 0.373893, max 0.619462 (between atoms 1747 and 1748) >> bonds that rotated more than 30 degrees: >> atom 1 atom 2 angle previous, current, constraint length >> 4 5 74.80.1130 0.0501 0.1204 >> 11 12 75.00.1132 0.0510 0.1204 >> 18 19 74.90.1129 0.0497 0.1204 >> 25 26 75.60.1126 0.0488 0.1204 >> 32 33 74.90.1132 0.0513 0.1204 >> 39 40 75.70.1129 0.0503 0.1204 > > etc. There are 7 atoms (the methyl group is a united atom, otherwise there > would be 10) in each TTA molecule, so you'll notice that it's complaining > about the bond between the 4th and 5th atoms on every one of them. My best > guess is that my topology file has insufficient or incorrect constraints, > allowing this particular bond to rotate. The entire molecule is supposed to > be rigid. I did constrain all the distances, but not the angles or dihedrals. > One of the example topology files I used did the same: constrained distances > but not angles or dihedrals. I was skeptical because it seems to allow some > degrees of freedom but I assumed gromacs somehow took care of it. Even in > this instance, it doesn't make sense to me that only one of the bonds is > problematic if none of them are constrained. So, I'm unsure of what exactly > is being left out. The bond in question is between the double bonded C=O. My > molecule looks like this: > > H O >\ // > H - C ≡ C - C >/ \ > H O-H > > My topology file is available from a previous message in this thread (minus > the last force field correction), linked here from the archive: > https://www.mail-archive.com/gromacs.org_gmx-users@maillist.sys.kth.se/msg07776.html > > If anyone knows directl
Re: [gmx-users] Simulating Multiple Solute Particles
Thank you again, Dr. Lemkul! I went with the first option because my advisor wants me to use tip3p water. Unfortunately, I immediately got myself stuck on the next step, but this time I think the solution may not be as trivial. The command to output em.tpr seems to work without any errors. The command, > gmx mdrun -deffnm em gives me a lot of errors. (It stops after 1000). Here's a small sample: >Steepest Descents: > Tolerance (Fmax) = 1.0e+03 > Number of steps=5 > >Step 1, time 0.001 (ps) LINCS WARNING >relative constraint deviation after LINCS: >rms 0.373893, max 0.619462 (between atoms 1747 and 1748) >bonds that rotated more than 30 degrees: > atom 1 atom 2 angle previous, current, constraint length > 4 5 74.80.1130 0.0501 0.1204 > 11 12 75.00.1132 0.0510 0.1204 > 18 19 74.90.1129 0.0497 0.1204 > 25 26 75.60.1126 0.0488 0.1204 > 32 33 74.90.1132 0.0513 0.1204 > 39 40 75.70.1129 0.0503 0.1204 etc. There are 7 atoms (the methyl group is a united atom, otherwise there would be 10) in each TTA molecule, so you'll notice that it's complaining about the bond between the 4th and 5th atoms on every one of them. My best guess is that my topology file has insufficient or incorrect constraints, allowing this particular bond to rotate. The entire molecule is supposed to be rigid. I did constrain all the distances, but not the angles or dihedrals. One of the example topology files I used did the same: constrained distances but not angles or dihedrals. I was skeptical because it seems to allow some degrees of freedom but I assumed gromacs somehow took care of it. Even in this instance, it doesn't make sense to me that only one of the bonds is problematic if none of them are constrained. So, I'm unsure of what exactly is being left out. The bond in question is between the double bonded C=O. My molecule looks like this: H O \ // H - C ≡ C - C / \ H O-H My topology file is available from a previous message in this thread (minus the last force field correction), linked here from the archive: https://www.mail-archive.com/gromacs.org_gmx-users@maillist.sys.kth.se/msg07776.html If anyone knows directly how to fix the error, I would greatly appreciate it, but helping me understand the error would also be very useful: Does "relative constraint deviation" mean that a constraint has been broken, or that there is no constraint where there should be one? Also, when it's saying the bonds are rotating more than 30 degrees, is that the bond angle changing or does it mean that a dihedral is changing by that much? Thanks again in advance - I really appreciate the help! N.H. - Original Message - From: "Justin Lemkul" To: gmx-us...@gromacs.org Sent: Tuesday, November 4, 2014 7:45:39 AM Subject: Re: [gmx-users] Simulating Multiple Solute Particles On 11/4/14 12:45 AM, Nathan K Houtz wrote: > Thanks for your reply. However, I'm still confused. I thought that the > command: > > #include "oplsaa.ff/tip3p.itp" > > is a call to the opls-aa force field. If this is not the correct way to > include the force field parameters, how should I do that? > That line does not call a force field, it calls a topology for TIP3P, which makes use of OPLS-AA parameters (well, more correctly TIP3P parameters as translated into OPLS-AA atom types and functional form). Your options are either (1) #include "oplsaa.ff/forcefield.itp" at the start of the topology and remove the [defaults] directive in your TETR topology or (2) translate TIP3P into atom types that are self-contained in your topology. -Justin -- == Justin A. Lemkul, Ph.D. Ruth L. Kirschstein NRSA Postdoctoral Fellow Department of Pharmaceutical Sciences School of Pharmacy Health Sciences Facility II, Room 629 University of Maryland, Baltimore 20 Penn St. Baltimore, MD 21201 jalem...@outerbanks.umaryland.edu | (410) 706-7441 http://mackerell.umaryland.edu/~jalemkul == -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Simulating Multiple Solute Particles
Thanks for your reply. However, I'm still confused. I thought that the command: #include "oplsaa.ff/tip3p.itp" is a call to the opls-aa force field. If this is not the correct way to include the force field parameters, how should I do that? Regards, Nathan Houtz - Original Message - From: "Justin Lemkul" To: gmx-us...@gromacs.org Sent: Monday, November 3, 2014 12:44:52 PM Subject: Re: [gmx-users] Simulating Multiple Solute Particles On 11/3/14 12:50 AM, Nathan K Houtz wrote: > Sorry everybody, I accidentally unsubscribed from the mailing list and missed > my last response. Dr. Lemkul posted a message but I can't reply directly to > it, so here's what he said: > >> Your topology is constructed incorrectly, but without seeing it in its >> entirety, it is impossible to say. Please post the file for download >> somewhere and provide a link, otherwise copy and paste its entire contents >> into a reply if the message will be small enough. >> >> >> -Justin >> >> -- >> == >> >> Justin A. Lemkul, Ph.D. >> Ruth L. Kirschstein NRSA Postdoctoral Fellow >> >> Department of Pharmaceutical Sciences >> School of Pharmacy >> Health Sciences Facility II, Room 629 >> University of Maryland, Baltimore >> 20 Penn St. >> Baltimore, MD 21201 >> >> jalem...@outerbanks.umaryland.edu | (410) 706-7441 >> http://mackerell.umaryland.edu/~jalemkul > > Thank you Justin! I think it's easier to just paste it here. This is my > topol.top file: > >> [ defaults ] >> ; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ >> 1 2 yes 0.5 0.5 >> >> [ atomtypes ] >> ;nameat.num mass charge ptype sigma epsilon >> CH3 15 15.035000.000 A 0.39100 0.66944 >> ;united atom Methyl [C-CH3] >> C 12 12.011000.000 A 0.35000 0.27614 ;C >> [R3CCOO-] (first carbon) >>CM 12 12.011000.000 A 0.37500 0.43932 >> ;carboxyl C [R-(C=O)-OH] >> O8 15.999400.000 A 0.30800 0.71128 >> ;carboxyl O [C-O-H] >>OM8 15.999400.000 A 0.29600 0.87864 >> ;carbonyl O [C=O] >> H11.007940.000 A 0.0 0.0 >> ;carboxyl H [C-O-H] >>OW8 15.999400.000 A 0.31660 0.65000 >>HW11.007940.000 A 0.0 0.0 >> >> [ moleculetype ] >> ; Namenrexcl >> TETR3 >> >> [ atoms ] >> ; nrtype resnr residuatomcgnr charge mass >> 1 CH3 1TETR C1 1 0 15.03500 >> 2 C 1TETR C2 1 0 12.01100 >> 3 C 1TETR C3 1 0 12.01100 >> 4 CM 1TETR C4 1 0.38 12.01100 >> 5 OM 1TETR O1 1 -0.38 15.99940 >> 6 O 1TETR O2 1 -0.40 15.99940 >> 7 H 1TETR H4 1 0.401.00794 >> >> [constraints] >> ; i j funct distance >> 1 2 10.14550 >> 2 3 10.11780 >> 3 4 10.14420 >> 4 5 10.12040 >> 4 6 10.13100 >> 6 7 10.89000 >> >> ; Include SPC water topology >> #include "oplsaa.ff/tip3p.itp" >> Here's your problem. Look at the contents of tip3p.itp: [ atoms ] ; idat type res nr residu name at name cg nr charge 1 opls_111 1 SOL OW 1 -0.834 2 opls_112 1 SOL HW1 10.417 3 opls_112 1 SOL HW2 10.417 It uses OPLS-AA atom types, but you haven't defined what those are because nothing above calls the OPLS-AA force field. Hence, grompp dies because it does not know how to assign nonbonded parameters to opls_111 or opls_112. >> [ system ] >> Tetrolic Acid in water >> >> [ molecules ] >> Tetrolic Acid 500 >> SOL 5481 > > I made the file myself, so I hope I didn't make any syntax errors. I > discovered with the .gro file that the spacing is very important since > gromacs doesn't
Re: [gmx-users] Simulating Multiple Solute Particles
Sorry everybody, I accidentally unsubscribed from the mailing list and missed my last response. Dr. Lemkul posted a message but I can't reply directly to it, so here's what he said: >Your topology is constructed incorrectly, but without seeing it in its >entirety, it is impossible to say. Please post the file for download somewhere >and provide a link, otherwise copy and paste its entire contents into a reply >if the message will be small enough. > > >-Justin > >-- >== > >Justin A. Lemkul, Ph.D. >Ruth L. Kirschstein NRSA Postdoctoral Fellow > >Department of Pharmaceutical Sciences >School of Pharmacy >Health Sciences Facility II, Room 629 >University of Maryland, Baltimore >20 Penn St. >Baltimore, MD 21201 > >jalem...@outerbanks.umaryland.edu | (410) 706-7441 >http://mackerell.umaryland.edu/~jalemkul Thank you Justin! I think it's easier to just paste it here. This is my topol.top file: >[ defaults ] >; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ > 1 2 yes 0.5 0.5 > >[ atomtypes ] >;name at.num mass charge ptype sigma epsilon > CH3 15 15.035000.000 A 0.39100 0.66944 > ;united atom Methyl [C-CH3] >C 12 12.011000.000 A 0.35000 0.27614 ;C > [R3CCOO-] (first carbon) > CM 12 12.011000.000 A 0.37500 0.43932 > ;carboxyl C [R-(C=O)-OH] >O8 15.999400.000 A 0.30800 0.71128 > ;carboxyl O [C-O-H] > OM8 15.999400.000 A 0.29600 0.87864 > ;carbonyl O [C=O] >H11.007940.000 A 0.0 0.0 > ;carboxyl H [C-O-H] > OW8 15.999400.000 A 0.31660 0.65000 > HW11.007940.000 A 0.0 0.0 > >[ moleculetype ] >; Namenrexcl > TETR3 > >[ atoms ] >; nrtype resnr residuatomcgnr charge mass > 1 CH3 1TETR C1 1 0 15.03500 > 2 C 1TETR C2 1 0 12.01100 > 3 C 1TETR C3 1 0 12.01100 > 4 CM 1TETR C4 1 0.38 12.01100 > 5 OM 1TETR O1 1 -0.38 15.99940 > 6 O 1TETR O2 1 -0.40 15.99940 > 7 H 1TETR H4 1 0.401.00794 > >[constraints] >; i j funct distance >1 2 10.14550 >2 3 10.11780 >3 4 10.14420 >4 5 10.12040 >4 6 10.13100 >6 7 10.89000 > >; Include SPC water topology >#include "oplsaa.ff/tip3p.itp" > >[ system ] >Tetrolic Acid in water > >[ molecules ] >Tetrolic Acid500 >SOL 5481 I made the file myself, so I hope I didn't make any syntax errors. I discovered with the .gro file that the spacing is very important since gromacs doesn't use delimiters but allows only a fixed number of characters for each entry. Is the .top file the same way? I just used a template from some example file and tried to make the numbers fit appropriately under the headings. If it's not correct, why didn't gromacs complain when I solvated it? Doesn't it read in the topology file then as well? I'm subscribed again so I won't miss a reply. Thanks again for the guidance! N. H. - Original Message - From: "Nathan K Houtz" To: gmx-us...@gromacs.org Sent: Sunday, October 26, 2014 10:49:35 PM Subject: Re: [gmx-users] Simulating Multiple Solute Particles Thanks for your help. And actually my version of gromacs does not have genbox anymore. (http://www.gromacs.org/Documentation/How-tos/Tool_Changes_for_5.0 at the bottom: GENBOX - This tool has been split to gmx solvate and gmx insert-molecules.) But it's no big deal, insert-molecules seemed to work just fine! Here are the commands I used: >>gmx insert-molecules -f tetrolic_acid.gro -ci tetrolic_acid.gro -nmol 499 -o >>tetrolic500.gro >>gmx solvate -cp tetrolic500.gro -cs spc216.gro -o tetrolic_solv.gro -p >>topol.top where tetrolic_acid.gro and topol.top are files I created myself to define the molecule. I checked the output in vmd and it looks like exactly what I want. But now I'm stuck on another step and I get a confusing error. I tried to minimize the energy via: >>gmx grompp -f minim.mdp -c tetrolic_solv.gro -p top
Re: [gmx-users] Simulating Multiple Solute Particles
Thanks for your help. And actually my version of gromacs does not have genbox anymore. (http://www.gromacs.org/Documentation/How-tos/Tool_Changes_for_5.0 at the bottom: GENBOX - This tool has been split to gmx solvate and gmx insert-molecules.) But it's no big deal, insert-molecules seemed to work just fine! Here are the commands I used: >>gmx insert-molecules -f tetrolic_acid.gro -ci tetrolic_acid.gro -nmol 499 -o >>tetrolic500.gro >>gmx solvate -cp tetrolic500.gro -cs spc216.gro -o tetrolic_solv.gro -p >>topol.top where tetrolic_acid.gro and topol.top are files I created myself to define the molecule. I checked the output in vmd and it looks like exactly what I want. But now I'm stuck on another step and I get a confusing error. I tried to minimize the energy via: >>gmx grompp -f minim.mdp -c tetrolic_solv.gro -p topol.top -o em.tpr I stole the minim.mdp file from a tutorial (http://bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme/Files/minim.mdp). I don't know how to properly create an mdp file from scratch, but judging from the comments, it seems reasonable and I don't see anything that must be specific to the case in the tutorial. Anyway, I don't think it caused my error. The output for that command is: Fatal error: Atomtype opls_111 not found For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors The link is not helpful, however. Googling exactly the error turns up very few results. The problem seems to be either with my topology file or with the itp file for water, which is referenced in my topology file like this: #include "oplsaa.ff/tip3p.itp" right after the constraints section. I also tried <#include "oplsaa.ff/spc.itp"> but it came up with a nearly identical error: Fatal error: Atomtype opls_116 not found For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors The closest thing to help I found by searching was this blog post: http://www.somewhereville.com/?p=114 However, his initial spc.itp file looks much different than mine. I think he must just have an older version of gromacs or something. Does anyone else know how I could do the energy minimization? Thanks again for all your help! N.H. - Original Message - From: "Mark Abraham" To: "Discussion list for GROMACS users" Sent: Saturday, October 25, 2014 6:10:21 AM Subject: Re: [gmx-users] Simulating Multiple Solute Particles On Sat, Oct 25, 2014 at 7:23 AM, Nathan K Houtz wrote: > Hello, I apologize for any ignorance but I'm quite new to gromacs and am > confused about a few things. > > I want to run some simulations of a small molecule, tetrolic acid > (CH3-C≡C-COOH) dissolved in a box of water. I'm doing this to see if a > united atom approximation on the methyl group is sufficiently good for > saving a bit of computation time. I have created a .gro file and a .top > file and attempted to put it into gromacs but I don't know how to model > more than one molecule of TTA. All of the tutorials I've looked up involve > solvating enormous proteins and they only do one at a time. I'd like to > simulate hundreds of TTA molecules with an appropriate amount of water > molecules. I attempted to do this by specifying the numbers of TTA and H2O > in my topology file, but when I tried to solvate it, gromacs removed the > line specifying how many waters and replaced it with its own number (216) > and made the box way too small to fit any more TTA. I feel like this may be > trivial, but I'm stuck. I really appreciate anyone's help. > You can adapt http://www.gromacs.org/Documentation/How-tos/Mixed_Solvents to do this. > After I figure this out, I'm going to try do dissolve tetrolic acid in > chloroform, then ethanol as well. I suspect gromacs doesn't have a built in > .gro or .top files for either of those (certainly not chloroform). Would > the solvate command work just as well if I created my own? > Yes, but perhaps gmx solvate -p won't quite do what you want. Do equilbrate that box decently; garbage in -> garbage out. Mark > Thanks in advance! > N.H. > -- > Gromacs Users mailing list > > * Please search the archive at > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before > posting! > > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists > > * For (un)subscribe requests visit > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or > send a mail to gmx-users-requ...@gromacs.org. > -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can
[gmx-users] Simulating Multiple Solute Particles
Hello, I apologize for any ignorance but I'm quite new to gromacs and am confused about a few things. I want to run some simulations of a small molecule, tetrolic acid (CH3-C≡C-COOH) dissolved in a box of water. I'm doing this to see if a united atom approximation on the methyl group is sufficiently good for saving a bit of computation time. I have created a .gro file and a .top file and attempted to put it into gromacs but I don't know how to model more than one molecule of TTA. All of the tutorials I've looked up involve solvating enormous proteins and they only do one at a time. I'd like to simulate hundreds of TTA molecules with an appropriate amount of water molecules. I attempted to do this by specifying the numbers of TTA and H2O in my topology file, but when I tried to solvate it, gromacs removed the line specifying how many waters and replaced it with its own number (216) and made the box way too small to fit any more TTA. I feel like this may be trivial, but I'm stuck. I really appreciate anyone's help. After I figure this out, I'm going to try do dissolve tetrolic acid in chloroform, then ethanol as well. I suspect gromacs doesn't have a built in .gro or .top files for either of those (certainly not chloroform). Would the solvate command work just as well if I created my own? Thanks in advance! N.H. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
