[gmx-users] Order parameter
Dear all, I am using GROMACS-4.6.5 for studying gas hydrates and related systems. A useful order parameter which indicates water to hydrate transition is the F4 order parameter where F4 is related to the H--O--O--H dihedral angle for a pair of water molecules. It is found that F4 parameter is very frequently mentioned in several works on gas hydrates done using GROMACS. I want to know whether GROMACS can calculate the value of the F4 parameter. Or is it possible for the user to define an order parameter? The g_hydorder tool in GROMACS refers to an order parameter different in definition from the F4 parameter. (P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.) Please comment. Regards, Sujith. -- 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] F4-Order-parameter
Dear all, I am using GROMACS-4.6.5 for studying gas hydrates and related systems. A useful order parameter which indicates water to hydrate transition is the F4 order parameter where F4 is related to the H--O--O--H dihedral angle for a pair of water molecules Exact expression for F4 can be found in Phys. Chem. Chem. Phys., 2015, 17, 9509-9518. The g_hydorder tool in GROMACS refers to an order parameter different in definition from the F4 parameter. (P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.) It is found that F4 parameter is very frequently used in several works on gas hydrates done using GROMACS (including the PCCP article mentioned above). I was wondering whether it is possible for any version of GROMACS to compute the value of F4 for my system (water + gas molecules). If not I guess I would have to write a program to do the same. Please comment. Regards, Sujith. -- 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] Pdb file
Hello, I guess it is not unusual for gaseous systems to have positive energy. I hope that the Fmax criteria is satisfied. Sujith. On Wed, Aug 12, 2015 at 10:30 AM, Sunil Ghimire ghimiresuni...@gmail.com wrote: I mean potential energy of system (1600 argon atoms) is positive after energy mininization . On 12 Aug 2015 16:58, Vitaly V. Chaban vvcha...@gmail.com wrote: L_BFGS On Tue, Aug 11, 2015 at 11:37 PM, Sunil Ghimire ghimiresuni...@gmail.com wrote: I created the pdb file by genconf but the system was not minimized, i got positive potential energy. What may be the problem? On 11 Aug 2015 19:38, Mark Abraham mark.j.abra...@gmail.com wrote: Hi, gmx genconf on a small box with a single argon is a useful starting point. Mark On Tue, Aug 11, 2015 at 3:16 PM Sunil Ghimire ghimiresuni...@gmail.com wrote: How can i create the .pdb file for 1600 atoms of argon ? -- 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. -- 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.
[gmx-users] Unexpected free energy values.
Dear all, I am studying the CH4-CO2-H2O mixture at low T(270K) and high P (20 bar). Forcefields used: CH4 : Single point LJ site; CO2 : EPM2 (Rigid with two virtual sites); H2O: TIP4P. I computed the free energy of solvation for a CH4 molecule in H2O in the presence of CO2 at different concentrations. The BAR method was used for the calculation. I decoupled the CH4 molecule in 20 steps with lambda varying by 0.05 each step. Experiments have revealed that CO2 and CH4 are cosolvents for each other in water and CH4 solubility will increase in presence of CO2. However, I am getting more positive value for the free energy of solvation of CH4 when the CO2 concentration is increased. Based on the experimental data one would expect that the free energy becomes less positive (more negative) since CO2 was observed to enhance CH4 solubility. The forcefield combination used here was used earlier for studying CH4-CO2-H2O system, though in a different context where the authors simulated hydrate formation. So I guess the morels used are fine. I wonder whether the problems lies in my MDP file . The free energy calculation parameters used are given below. ; Free energy control stuff free_energy = yes init_lambda = 0.0 delta_lambda = 0 foreign_lambda = 0.05 sc-alpha = 0.5 sc-power = 1.0 sc-sigma = 0.373 ; equal to sigma of the single point CH4 model. couple-moltype = CH4 couple-lambda0 = vdw couple-lambda1 = none couple-intramol = no nstdhdl = 10 Any comment will be of great help. Regards, Sujith. -- 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] Regarding-virtual-sites-tutorial
Hello Justin, Thanks for the mail. I am interested in knowing whether the particle type assigned has any role in the way it is treated by the program. For both virtual sites and the regular atoms, interaction parameters and mass are assigned. So in Gromacs, does the particle type (A or V) has any role beyond these interaction parameters. Regards, Sujith. -- 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] Regarding-virtual-sites-tutorial
Dear Justin, I am following your tutorial about CO2 topology using virtual sites. I find that in the topology provided, you have not reassigned the particle type as virtual (V) for the C and O atoms (which are treated as the virtual sites in the model). Instead, they are directly selected from the OPLSAA forcefield files, in which the particle type is always 'A'. I was wondering if the choice of particle type has any effect on the results. Please comment. Regards Sujith. -- 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] About airtual site tutorial
Dear Justin, I am following your tutorial about CO2 topology using virtual sites. I find that in the topology provided, you have not reassigned the particle type as virtual (V) for the C and O atoms (which are treated as the virtual sites in the model). Instead, they are directly selected from the OPLSAA forcefield files, in which the particle type is always 'A'. I was wondering if the choice of particle type has any effect on the results. Please comment. Regards Sujith. -- 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] Location of Bubbles.
Hello, I am observing nucleation and growth of bubbles in a solution of methane in water with methane present at various levels of supersaturation. The system consists of a total 0f 2000 molecules in a cubic box with the number of CH4 varying from 80 to 200. I am using a single point Lennard Jones potential for methane along with the TIP4P water model with short range interactions cutoff at 1.2nm. Five different simulations were performed and in each case a methane bubble was formed . I find that in all cases the bubble is located near the edge or a face of the box rather than being within the box. I wonder whether this is a coincidence or if it is an artefact. I am interested to know if someone had similar observations. Please comment. Regards, Sujith. -- 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] Location of Bubbles.
Thank you Mark and Jan. On Sat, Jul 18, 2015 at 3:59 AM, Mark Abraham mark.j.abra...@gmail.com wrote: Hi, If your box is periodic then there are no edges or faces... Also note that if you mentally divide your box in e.g. three parts in each dimension, the 27 resulting cubes have 26 that are near an original edge or face, so of course things will look like they happen near edges and faces... Mark On Sat, 18 Jul 2015 09:49 Jan Jirsák janjir...@gmail.com wrote: Hi Sujith, you can try to translate the initial configuration by, say, L/2 in all directions using trjconv. If the bubble appears again at the edge, I would suspect artifacts. If you are using Verlet cutoff scheme, you can also try to change it to group - it helped in my simulation when a droplet was pulled to the origin by a spurious force ( http://gromacs.org_gmx-users.maillist.sys.kth.narkive.com/6qnoy5o4/strange-behavior-of-water-droplet-on-a-solid-surface-with-verlet-cutoff-scheme-in-gmx-5-0-4 ). Regards, Jan 2015-07-18 8:53 GMT+02:00 sujithkakkat . sujithk...@gmail.com: Hello, I am observing nucleation and growth of bubbles in a solution of methane in water with methane present at various levels of supersaturation. The system consists of a total 0f 2000 molecules in a cubic box with the number of CH4 varying from 80 to 200. I am using a single point Lennard Jones potential for methane along with the TIP4P water model with short range interactions cutoff at 1.2nm. Five different simulations were performed and in each case a methane bubble was formed . I find that in all cases the bubble is located near the edge or a face of the box rather than being within the box. I wonder whether this is a coincidence or if it is an artefact. I am interested to know if someone had similar observations. Please comment. Regards, Sujith. -- 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. -- 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] Thermostats and MSD
Dear all, I am simulating a system of small hydrophobic solutes ( CH4) dissolved in water. Simulations are performed in the NPT simulation, and I would like to determine the diffusion coefficients from MSD plots. I have two related questions; (1) Is the thermostat going to seriously affect the diffusivity of the molecules. I am using Nose-Hoover thermostat with a time constant of 0.2 ps. Is the choice of time constant very crucial while studying diffusivity. Should I avoid thermostat and go for NVE simulations for studying diffusivity. (2) From radial distribution functions, it is found that there is association between the hydrophobic solvents. I would like to know, if in case a small cluster of the hydrophobic solutes are formed, then the whether the rotational motion of these clusters contribute to dispalcement in the MSD calculation. Please comment. Regards, Sujith. -- 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] Question on forcefields?
Dear all, Is it possible that a model for a simple molecule (say CO2) which can accurately predict macroscopic properties (like density), but still can have error in the microscopic behavior, like diffusivity or the pair interaction functions ? Please share your thoughts. Regards, Sujith. -- 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] Virtual sites and diffusion coefficients.
Hello Mark, I have tested both models in simulating a variety of properties. The original parametrization results involved predicting the critical point of CO2 and also density at various P and T. I did several NVT simulations to predict density-vs-pressure data. In this case the model with virtual interaction sites gave much better results. In addition to this, simulations performed were performed on pure CO2 at supercritical conditions and mean square displacement was calculated for the CO2 molecule. When compared to the reported ab initio dynamics results the one obtained using virtual sites gave much better agreement. The density of CO2 dissolved in H2O at very dilute concentrations was also studied using both models of CO2 and the results were in good agreement with experimental data for the model with virtual sites. In this case the one without virtual sites also did reasonably well. However, the disparity appears when considering the Measn Square displacement of dilute CO2 solution in water. Here the model without virtual sites is closer to the experimental data. This has become very confusing now. The original parametrization of he model did not involve predicting the transport properties. Also in the published results using the EPM2 model I don't find any mention of the virtual interaction sites. Regards, Sujith. On Tue, Feb 3, 2015 at 1:40 PM, Mark Abraham mark.j.abra...@gmail.com wrote: Hi, With which set of parameters are you able to reproduce the original parametrization results? Mark On Tue, Feb 3, 2015 at 8:04 AM, sujithkakkat . sujithk...@gmail.com wrote: Dear all, I tried simulating the CO2-water mixture with CO2 at a low concentration (0.003 mole fraction) at various temperatures at 1 atmos. I used the EPM2 parameters for the CO2 molecules which includes a flexible OCO angle. I did a separate simulation on the similar system with CO2 model constructed using virtual sites (following Justin Lemkul's tutorial) , with the interaction parameters same as the EPM2 model. In both cases the water model TIP4P was used. I ran 5ns simulation after equilibrium was attained for the T and P values. When plotting the mean square displacement and fitting the graph to get the diffusion coefficients, what I find is that the diffusion coefficient obtained in the case of the CO2 model using virtual sites is much higher (approximately two times) of that obtained for the EPM2 model without the virtual sites. Let me know if any one came across this before. I would like to hear from the more experienced users, what they think about this. I will provide further details on the topology and MD parameters if necessary. Regards, Sujith. -- 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.
[gmx-users] Questions about CO2.
Dear all, I went through the mails about the use of virtual sites for representing linear molecules like CO2. I also read that GROMACS cannot handle angle bending potentials at 180 degrees in an appropriate manner. My system consists of EPM2-CO2 dissollved in TIP4- H2O. I made topology files for CO2 with virtual sites and also made another topology without the use of virtual sites providing the equilibrium angle (180) and the bending force constant. Both these are given at the end. I am trying to use the EPM2 model for CO2 reported in* J. Phys. Chem. 1995,99, 12021-12024* by *Jonathan G. Harris and Kwong H. Yung* under the title Carbon Dioxide's Liquid-Vapor Coexistence Curve and Critical Properties As Predicted by a Simple Molecular Model. I have the following questions; (i) I am using gromacs 4.6.5. Is the problem related to treating angle bending potentials for linear molecules solved in later versions of gromacs? (ii) The authors Harris and Yung stress that their model is flexible in terms of bending in the abstract of their article. So if I build this model using virtual sites then I believe I am neglecting the flexibility. Also while using virtual sites the model looks more like a dimer with two point masses, instead of a a three point EPM2 model. So can I call this an EPM2 model any more.? Do you think introducing virtual sites affects the results which are expected from a simple looking EPM2 model, with three charged-Lennard-Jones point masses with a flexible angle. I appreciate (and often surprised by) the fact that a simple model like EPM2 which can be specified by a total of just nine parameters can accurately predict the properties of CO2 at supercritical conditions which otherwise may a be a very intricate condensed system to study. This prompts me to think that any small change in the model can impart an error that can get magnified for a system of thousands of molecules. (iii) I found that I am getting highly positive energy values when using the topology built using virtual sites. To add to this I am forced to use very short time steps (0.2 fs) for simulations, since with 2fs, system was unstable which I believe has to do with the positive energy and high force values. Whereas, using the topology built in the usual way without virtual sites, I get negative energy values and simulations could easily be performed using a 2fs time step (system may be more stable) . I checked the angles in this case at the final state of the system, and found them lying between 175 to 180 degrees, which I am not sure if it is all right. Is this slight deviation from 180 degree caused by the failure of GROMACS to maintain linearity? Or, is it perfectly normal since the model itself is flexible.? Please share your views. Regards, Sujith. *TOPOLOGY USING VIRTUAL SITES:* [atomtypes] ; name mass charge ptypesigma epsilon D 22.0049 0. A 0. 0. CA 0.0.6512 A 0.2757 0.2339 CO 0. -0.3256 A 0.3033 0.6695 [moleculetype] ; name nrexcl CO22 [atoms] ; nr type resnr residue atom cgnr charge mass 1 D 1 CO2 D1 10. 22.0049 2 D 1 CO2 D2 10. 22.0049 3 CA1 CO2CA 10.6512 0. 4 CO1 CO2OC11 -0.3256 0. 5 CO1 CO2OC21 -0.3256 0. [constraints] ; i j funct doc 1 2 1 0.195948 [virtual_sites2] ; i j k funct a 3 1 2 1 0.5 4 1 2 1 1.08638006 5 2 1 1 1.08638006 *TOPOLOGY WITHOUT VIRTUAL SITES:* [ defaults ] ; nbfunc comb-rule gen-pairsfudgeLJ fudgeQQ 1 2 no 1.0 1.0 [atomtypes] ;name mass charge ptype sigma epsilon CA 12.011000.6512 A 0.27570.2339 CO 15.99940 -0.3256 A 0.30330.6695 [ nonbond_params ] ; i j funct sigmaepsilon CA CO10.2892 0.3955 [moleculetype] ; name nrexcl CO22 [atoms] ; nr type resnr residue atom cgnr charge mass 1 CA1 CO2 CA 1 0.6512 12.0110 2 CO1 CO2 OC1 1 -0.3256 15.9994 3 CO1 CO2 OC2 1 -0.3256 15.9994 [angles] ; i j k funct ao ak 3 1 2 1 180 1236 [constraints] ; i j funct doc 1 2 1 0.1149 1 3 1 0.1149 [exclusions] 1 2 3 2 1 3 3 1 2 -- 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] PMF and RDF
Hi Erik and Justin, Thanks for the response. Erik, I was thinking that in the case of hydrophobic solutes, there is a higher chance of proper sampling of all points along the inter solute distance. I believe the water shell around the hydrophibic solvent can break easier than the that in the case of ions. In that case I hope standard simulations can give better solute-solute RDFs and PMF may be avoided. The article which I refered to studies hydrophobic solutes and they have done PMF calculations. Regards, Sujith. On Tue, Dec 9, 2014 at 7:41 PM, Erik Marklund erik.markl...@chem.ox.ac.uk wrote: On 9 Dec 2014, at 11:00, Erik Marklund erik.markl...@chem.ox.ac.uk wrote: Dear Sujith, Umbrella sampling does exactly that, but adds a biasing potential to sample high-energy regions of the reaction coordinate in separate simulations. The -px output when you do umbrella sampling with mdrun is indeed a sampling of distances along the reaction coordinate, which if you histogram it is the rdf. g_wham uses that data (or alternatively the force along the reaction coordinate instead of the distance) to build a PMF. My question to you is why do you think rdfs are significantly easier than pmts? To clarify what I meant, in light of what Justin wrote earlier, is that the pmf and ref can *in principle* be calculated from exactly the same data. The *conventional* way to get the pmf involves biasing the simulations to sample poorly populated regions, whereas the rdf is commonly inferred from simulations without such bias. In principle the rdf can also be obtained from biased simulations if that bias is corrected for during analysis, and similarly a pmf can be obtained without bias. As Justin said, there are many situations where the sampling is insufficient for some distances however. Erik Kind regards, Erik On 9 Dec 2014, at 06:39, sujithkakkat . sujithk...@gmail.com wrote: Dear all, I read in *Phys. Chem. Chem. Phys., 2009, 11, 10427-10437*, that the radial distribution function is directly related to Potential of mean force through RDF=exp(-PMF/kT). My question is why would someone worry about computing PMF in a simple case like interaction between two small solute molecules in water , along the intermolecular distance, when one can get the RDF between the solutes , which I believe is easier than PMF calculation. Another article *Biophysical Chemistry 101-102 (2002), 295-307 *reports PMF between solute molecules from Monte Carlo simulations. Why not just find RDF. Regards, Sujith. -- 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. -- 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] PMF and RDF
Dear all, I read in *Phys. Chem. Chem. Phys., 2009, 11, 10427-10437*, that the radial distribution function is directly related to Potential of mean force through RDF=exp(-PMF/kT). My question is why would someone worry about computing PMF in a simple case like interaction between two small solute molecules in water , along the intermolecular distance, when one can get the RDF between the solutes , which I believe is easier than PMF calculation. Another article *Biophysical Chemistry 101-102 (2002), 295-307 *reports PMF between solute molecules from Monte Carlo simulations. Why not just find RDF. Regards, Sujith. -- 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] Problems with CO2-water simulation
Dear all, I am trying to simulate a water-carbondioxide system at 250K and 30bar. The system consists of 85 CO2 and 1204 water molecules in a 3x3x4 nm box. I used EPM2 and TIP4P potentials for CO2 and water following *Energy Environ. Sci., 2012, 5, 7033-7041* where the authors used the same.\ I am facing problems with the energy potential energy values which are very high and positive ( of the order 10^7). I tried the same simulation with the OPLSAA topology of CO2 which I obtained from Dr. Justin Lemkul's tutorial on virtual site with TIP4P water. The problem still remains, potential energy values being highly positive. I am totally stuck at this point and is unable to find how to tackle the problem. Firstly, I am wondering whether the positive potential energy is something to be worried about. If so, it would be great if you could provide me some direction in which I should think, or suggest me something to read about on this issue. Given below is the parameter file I used (I have redueced the step size due to a avoid persistent blowing up problem); title = co2--in-water ; Run parameters integrator= md nsteps= 250 dt = 0.0002 ; Output control nstxout = 5000 nstvout = 5000 nstenergy= 5000 nstlog = 5000 nstxtcout = 5000 ; Bond parameters continuation = no constraint_algorithm= lincs constraints = all-bonds lincs_iter = 2 lincs_order = 4 ; Neighborsearching cutoff-scheme = Verlet ns_type= grid nstlist = 10 rcoulomb = 1.2 rvdw = 1.2 ; Electrostatics coulombtype = PME pme_order = 4 fourierspacing = 0.16 ; Temperature coupling is on tcoupl = Nose-Hoover tc-grps= CO2 SOL tau_t = 0.4 0.4 ref_t= 250 250 ; Pressure coupling is on pcoupl = Parrinello-Rahman pcoupltype= isotropic tau_p= 2.0 ref_p = 30.0 compressibility = 4.5e-5 refcoord_scaling = com ; Periodic boundary conditions pbc= xyz ; Dispersion correction DispCorr = EnerPres ; Velocity generation gen_vel = no gen_temp= 250 gen_seed= -1 The EPM2 topology file which I made is given below; [atomtypes] ; namemasscharge ptypesigmaepsilon D 22.0049 0. A 0. 0. CA 0. 0.6512 A 0.2757 0.2339 CO 0. -0.3256 A 0.3033 0.6695 [moleculetype] ; name nrexcl CO22 [atoms] ; nrtyperesnr residueatomcgnr chargemass 1 D 1 CO2D1 1 0.22.0049 2 D 1 CO2D2 1 0.22.0049 3 CA 1 CO2CA 1 0.6512 0. 4 CO 1 CO2OC1 1-0.3256 0. 5 CO 1 CO2OC2 1-0.3256 0. [constraints] ; i j functdoc 1 2 1 0.195948 [virtual_sites2] ; i j k funct a 3 1 2 1 0.5 4 1 2 1 1.08638006 5 2 1 1 1.08638006 Regards, Sujith. -- 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] Problems with CO2-water simulation
Hello Tsjerk, Thanks for the response. I find your argument intuitive. But my simulations are at 30bars and 250K. Isn't the conditions supposed to play a role in making the dissolved state metastable. My system at 250K and 30 bar is like a closed coke can in the refrigerator. I believe there are no bubbles in the coke can until we open it to release the pressure. Also, I have performed simulations with methane-water systems at similar conditions. There I found the potential energy is negative. Since , I guess CO2 to be more soluble in water, doesn't that make the positive potential energy in CO2-water case unusual. Regards, Sujith. On Mon, Dec 8, 2014 at 12:35 PM, Tsjerk Wassenaar tsje...@gmail.com wrote: Hi Sujith, There is nothing special about positive potential energy. In this case it explains why soft drinks gives bubbles: the CO2 wants to get out. In addition, if it could, it would drive the system to HCO3-/H3O+, which is an energetically more favorable way to store CO2 in water. Cheers, Tsjerk On Dec 8, 2014 6:50 AM, sujithkakkat . sujithk...@gmail.com wrote: Dear all, I am trying to simulate a water-carbondioxide system at 250K and 30bar. The system consists of 85 CO2 and 1204 water molecules in a 3x3x4 nm box. I used EPM2 and TIP4P potentials for CO2 and water following *Energy Environ. Sci., 2012, 5, 7033-7041* where the authors used the same.\ I am facing problems with the energy potential energy values which are very high and positive ( of the order 10^7). I tried the same simulation with the OPLSAA topology of CO2 which I obtained from Dr. Justin Lemkul's tutorial on virtual site with TIP4P water. The problem still remains, potential energy values being highly positive. I am totally stuck at this point and is unable to find how to tackle the problem. Firstly, I am wondering whether the positive potential energy is something to be worried about. If so, it would be great if you could provide me some direction in which I should think, or suggest me something to read about on this issue. Given below is the parameter file I used (I have redueced the step size due to a avoid persistent blowing up problem); title = co2--in-water ; Run parameters integrator= md nsteps= 250 dt = 0.0002 ; Output control nstxout = 5000 nstvout = 5000 nstenergy= 5000 nstlog = 5000 nstxtcout = 5000 ; Bond parameters continuation = no constraint_algorithm= lincs constraints = all-bonds lincs_iter = 2 lincs_order = 4 ; Neighborsearching cutoff-scheme = Verlet ns_type= grid nstlist = 10 rcoulomb = 1.2 rvdw = 1.2 ; Electrostatics coulombtype = PME pme_order = 4 fourierspacing = 0.16 ; Temperature coupling is on tcoupl = Nose-Hoover tc-grps= CO2 SOL tau_t = 0.4 0.4 ref_t= 250 250 ; Pressure coupling is on pcoupl = Parrinello-Rahman pcoupltype= isotropic tau_p= 2.0 ref_p = 30.0 compressibility = 4.5e-5 refcoord_scaling = com ; Periodic boundary conditions pbc= xyz ; Dispersion correction DispCorr = EnerPres ; Velocity generation gen_vel = no gen_temp= 250 gen_seed= -1 The EPM2 topology file which I made is given below; [atomtypes] ; namemasscharge ptypesigmaepsilon D 22.0049 0. A 0. 0. CA 0. 0.6512 A 0.2757 0.2339 CO 0. -0.3256 A 0.3033 0.6695 [moleculetype] ; name nrexcl CO22 [atoms] ; nrtyperesnr residueatomcgnr chargemass 1 D 1 CO2D1 1 0.22.0049 2 D 1 CO2D2 1 0.22.0049 3 CA 1 CO2CA 1 0.6512 0. 4 CO 1 CO2OC1 1-0.3256 0. 5 CO 1 CO2OC2 1-0.3256 0. [constraints] ; i j functdoc 1 2 1 0.195948 [virtual_sites2] ; i j k funct
[gmx-users] question about modified LJ combination rule
Dear all, I am studying a water-methane system with tip4p-Ice for water and oplsaa for methane. What I want to try is modify the intermolecular non-bonded interaction between methane and water. I just want to multiply the combination rule (geometric mean) used to obtain the mean well depth (epsilon) component by a scaling factor, without changing any of the non bonding function types. I am aware that user defined non-bonding functions are possible. However, here I don't want to change the function type, but only slightly modify the combination rule 3 , by multiplying the epsilon part with a scaling factor. Also, I find that there is no nonbond_params section in OPLSAA where I think it would have been possible to edit cross LJ terms . So, is there an easy way I can do it rather than going for user defined non-bonded functions . Is it ok to define nonbond_params in topology file with oplsaa forcefield? Please comment. Thanks, Sujith. -- 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] RDF of solvation shell molecules
Hello, I would like to know if it is possible to plot the RDF of water molecules belonging to the solvation shell of a solute. I mean, I want to know the effect of the solute on the water arrangement in the solvation shell. Is there a way to do this? Regards, Sujith. -- 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] RDF of solvation shell molecules
Hello Chandan, Thank you for the response. I have used g_rdf, but I guess it gives the RDF for the whole system. Say, if I choose water oxygen, it will plot RDF from water molecules in the whole system. But, I want the RDF for the solvation shell water molecules ,leaving the bulk water . I came across such a plot in an article where the author studied the structure of solvation shell . So do you think it is possible to do that GROMACS? Regards, Sujith. On Thu, May 22, 2014 at 9:31 PM, Chandan Choudhury iitd...@gmail.comwrote: Hi Sujith, g_rdf is the tool you need. Chandan On Thu, May 22, 2014 at 8:18 PM, sujithkakkat . sujithk...@gmail.com wrote: Hello, I would like to know if it is possible to plot the RDF of water molecules belonging to the solvation shell of a solute. I mean, I want to know the effect of the solute on the water arrangement in the solvation shell. Is there a way to do this? Regards, Sujith. -- 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. -- -- Chandan Kumar Choudhury National Chemical Laboratory, Pune India -- 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] High density after genbox
Thank you Mark, My Bad. Now fixed. Hope the simulations I did before were not affected, I had defined atomtypes with the correct mass in the .itp file used. Sujith On Wed, May 21, 2014 at 12:59 AM, Mark Abraham mark.j.abra...@gmail.comwrote: Probably garbage in, garbage out. However you're measuring the density probably depends on share/top/atommass.dat, which relies on matching atom names to infer atom types and thus masses. If your atom names don't follow its assumptions... at least some tools warn about this in the output. Did the tool you were using not do so? Mark On Tue, May 20, 2014 at 12:50 PM, sujithkakkat . sujithk...@gmail.com wrote: Hello, I am working on a water-methane system with OPLSAA for methane and TIP4P/Ice for water. What I find strange is that after solvating the methane molecule with water (tip4p/ice) using the genbox routine, the density appears to be very high at 1661.04 (g/l) . The number of water molecules is normal for the box size, and therefore should give lower density values. After equilibration steps the density gets back to normal. Earlier , I ignored this problem and went ahead with the simulations, and the end results were found to be good. For eg, the solvation free energy of methane I calculated , was very close to reported values and also the value obtained in Justin's tutorial. Another test was simulating the tip4p/ice water , which gave results, for example the RDFs very close to the reported results. Could some one comment what is going on that gives the high density? Regards, Sujith. -- 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.
[gmx-users] High density after genbox
Hello, I am working on a water-methane system with OPLSAA for methane and TIP4P/Ice for water. What I find strange is that after solvating the methane molecule with water (tip4p/ice) using the genbox routine, the density appears to be very high at 1661.04 (g/l) . The number of water molecules is normal for the box size, and therefore should give lower density values. After equilibration steps the density gets back to normal. Earlier , I ignored this problem and went ahead with the simulations, and the end results were found to be good. For eg, the solvation free energy of methane I calculated , was very close to reported values and also the value obtained in Justin's tutorial. Another test was simulating the tip4p/ice water , which gave results, for example the RDFs very close to the reported results. Could some one comment what is going on that gives the high density? Regards, Sujith. -- 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] Positive energy .
Hello Justin, I looked at the energy components. What I find is that Short Range LJ is highly positive (208522 kJ/mol) . Also the short range coulomb interaction (-2028.46 kJ/mol) , although negative looks weak (I think so , because an earlier minimization of a system with methane( oplsaa) in water (tip4p_Ice) gave a more negative value for short range coulomb. I expect CO2 to have stronger coulombic interactions with water ) . Do you think any of the cutoff values in .mdp are really bad? I tried changing them but still kept getting poor results. Regards, Sujith. On Wed, May 14, 2014 at 8:58 PM, Justin Lemkul jalem...@vt.edu wrote: On 5/14/14, 12:39 AM, sujithkakkat . wrote: Hello Justin and others, Thank you for the response . Here are mode details of the system; ( in the order ; MDP settings, itp and top files) . The system consists of 460 water molecules and 1 CO2 molecule . I am trying to find the free energy of solvation for CO2 in water at various temperatures and pressure. *MDP settings:* ; Run control integrator = steep nsteps = 2 ; EM criteria and other stuff emtol = 10 emstep = 0.01 niter= 20 nbfgscorr = 10 ; Output control nstlog = 1 nstenergy = 1 ; Neighborsearching and short-range nonbonded interactions nstlist= 1 ns_type = grid pbc = xyz rlist = 1.0 ; Electrostatics coulombtype = PME rcoulomb = 1.0 ; van der Waals vdw-type = switch rvdw-switch =0.8 rvdw = 0.9 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; Spacing for the PME/PPPM FFT grid fourierspacing = 0.12 ; EWALD/PME/PPPM parameters pme_order = 6 ewald_rtol = 1e-06 epsilon_surface = 0 optimize_fft= no ; Temperature and pressure coupling are off during EM tcoupl = no pcoupl = no ; Generate velocities to start gen_vel = no ; options for bonds constraints = all-bonds ; Type of constraint algorithm constraint-algorithm = lincs ; Do not constrain the starting configuration continuation = no ; Highest order in the expansion of the constraint coupling matrix lincs-order = 12 *--- -* *.itp for tip4p/Ice* [atomtypes] ;name mass charge ptypesigmaepsilon IW0 0.000 D 0.0 0.0 OW4 15.99940 0.000 A 0.31668 0.88211 HW1.00800 0.000 A 0.0E+00 0.0E+00 [moleculetype] ; name nrexcl water 1 [atoms] ; nr type resnr residu atom cgnr charge mass 1 OW4 1 H2OOW1 1 0 15.994 2 HW 1 H2OHW2 1 0.58971.008 3 HW 1 H2OHW3 1 0.58971.008 4 IW 1 H2OMW4 1-1.17940.0 [constraints] ;i j funct doh dhh 1 2 1 0.09572 1 3 1 0.09572 2 3 1 0.15139 [exclusions] 1 2 3 4 2 1 3 4 3 1 2 4 4 1 2 3 ; The position of the dummy is computed as follows: ; ; O ; ; D ; ; H H ; ; const = distance (OD) / [ cos (angle(DOH))* distance (OH) ] ; 0.015 nm / [ cos (52.26 deg) * 0.09572 nm] ; Dummy pos x4 = x1 + a*(x2-x1) + b*(x3-x1) [virtual_sites3] ; Dummy fromfunct a b 4 1 2 3 1 0.13458 0.13458 *--- .itp for CO2* [atomtypes] ; namemasscharge ptypesigmaepsilon D 22.0049 0. A 0. 0. CA 0. 0.6512 A 0.2757 0.2339 CO 0. -0.3256 A 0.3033 0.6695 [moleculetype] ; name nrexcl CO22 [atoms] ; nrtyperesnr residueatomcgnr chargemass 1 D 1 CO2D1 1 0.22.0049 2 D 1 CO2D2 1 0.22.0049 3 CA 1 CO2CA 1 0.6512 0. 4 CO 1 CO2OC1 1-0.3256 0. 5 CO 1 CO2OC2 1-0.3256 0. [constraints] ; i j functdoc 1 2 1 0.195948 [virtual_sites2] ; i j k funct a 3 1 2 1 0.5 4 1 2 1 1.08638006 5 2 1 1 1.08638006
[gmx-users] Positive energy .
Dear all, I am running simulations for a water-CO2 system using the OPLSAA forcefield, where I use TIP4P/Ice water model and the EPM2 model for CO2. After energy minimization I found the value of potential energy to be positive, and also the maximum force value was high (I requested Fmax 10 ). The values are as follows; Potential Energy= 1.7600962e+05 Maximum force = 6.5140656e+01 on atom 1176 Norm of force= 8.7009411e+00 Earlier when I did a simulation for gaseous CO2 using the same model, the potential energies were always positive. But I doubt whether the energy can be positive here. I want to make sure that I am not using any wrong parameters, before I proceed. Please comment. My MDP file is as given below; ; Run control integrator = steep nsteps = 2 ; EM criteria and other stuff emtol = 10 emstep = 0.01 niter= 20 nbfgscorr = 10 ; Output control nstlog = 1 nstenergy = 1 ; Neighborsearching and short-range nonbonded interactions nstlist= 1 ns_type = grid pbc = xyz rlist = 1.0 ; Electrostatics coulombtype = PME rcoulomb = 1.0 ; van der Waals vdw-type = switch rvdw-switch =0.8 rvdw = 0.9 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; Spacing for the PME/PPPM FFT grid fourierspacing = 0.12 ; EWALD/PME/PPPM parameters pme_order = 6 ewald_rtol = 1e-06 epsilon_surface = 0 optimize_fft= no ; Temperature and pressure coupling are off during EM tcoupl = no pcoupl = no ; Generate velocities to start gen_vel = no ; options for bonds constraints = all-bonds ; Type of constraint algorithm constraint-algorithm = lincs ; Do not constrain the starting configuration continuation = no ; Highest order in the expansion of the constraint coupling matrix lincs-order = 12 -- 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] Positive energy .
Hello Justin and others, Thank you for the response . Here are mode details of the system; ( in the order ; MDP settings, itp and top files) . The system consists of 460 water molecules and 1 CO2 molecule . I am trying to find the free energy of solvation for CO2 in water at various temperatures and pressure. *MDP settings:* ; Run control integrator = steep nsteps = 2 ; EM criteria and other stuff emtol = 10 emstep = 0.01 niter= 20 nbfgscorr = 10 ; Output control nstlog = 1 nstenergy = 1 ; Neighborsearching and short-range nonbonded interactions nstlist= 1 ns_type = grid pbc = xyz rlist = 1.0 ; Electrostatics coulombtype = PME rcoulomb = 1.0 ; van der Waals vdw-type = switch rvdw-switch =0.8 rvdw = 0.9 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; Spacing for the PME/PPPM FFT grid fourierspacing = 0.12 ; EWALD/PME/PPPM parameters pme_order = 6 ewald_rtol = 1e-06 epsilon_surface = 0 optimize_fft= no ; Temperature and pressure coupling are off during EM tcoupl = no pcoupl = no ; Generate velocities to start gen_vel = no ; options for bonds constraints = all-bonds ; Type of constraint algorithm constraint-algorithm = lincs ; Do not constrain the starting configuration continuation = no ; Highest order in the expansion of the constraint coupling matrix lincs-order = 12 ** *.itp for tip4p/Ice* [atomtypes] ;name mass charge ptypesigmaepsilon IW0 0.000 D 0.0 0.0 OW4 15.99940 0.000 A 0.31668 0.88211 HW1.00800 0.000 A 0.0E+00 0.0E+00 [moleculetype] ; name nrexcl water 1 [atoms] ; nr type resnr residu atom cgnr charge mass 1 OW4 1 H2OOW1 1 0 15.994 2 HW 1 H2OHW2 1 0.58971.008 3 HW 1 H2OHW3 1 0.58971.008 4 IW 1 H2OMW4 1-1.17940.0 [constraints] ;i j funct doh dhh 1 2 1 0.09572 1 3 1 0.09572 2 3 1 0.15139 [exclusions] 1 2 3 4 2 1 3 4 3 1 2 4 4 1 2 3 ; The position of the dummy is computed as follows: ; ; O ; ; D ; ; H H ; ; const = distance (OD) / [ cos (angle(DOH))* distance (OH) ] ; 0.015 nm / [ cos (52.26 deg) * 0.09572 nm] ; Dummy pos x4 = x1 + a*(x2-x1) + b*(x3-x1) [virtual_sites3] ; Dummy fromfunct a b 4 1 2 3 1 0.13458 0.13458 *---.itp for CO2* [atomtypes] ; namemasscharge ptypesigmaepsilon D 22.0049 0. A 0. 0. CA 0. 0.6512 A 0.2757 0.2339 CO 0. -0.3256 A 0.3033 0.6695 [moleculetype] ; name nrexcl CO22 [atoms] ; nrtyperesnr residueatomcgnr chargemass 1 D 1 CO2D1 1 0.22.0049 2 D 1 CO2D2 1 0.22.0049 3 CA 1 CO2CA 1 0.6512 0. 4 CO 1 CO2OC1 1-0.3256 0. 5 CO 1 CO2OC2 1-0.3256 0. [constraints] ; i j functdoc 1 2 1 0.195948 [virtual_sites2] ; i j k funct a 3 1 2 1 0.5 4 1 2 1 1.08638006 5 2 1 1 1.08638006 *-* *.top file*#include Oplsaa.ff/forcefield.itp #include Oplsaa.ff/atomtype.itp #include Oplsaa.ff/tip4p_Ice.itp #include Oplsaa.ff/co2_epm2_dummy.itp [ system ] ; Name CO2 in water [ molecules ] ; Compound#mols CO2 1 SOL 460 *-* Thank you for your time. Regards, Sujith On Tue, May 13, 2014 at 9:57 PM, Justin Lemkul jalem...@vt.edu wrote: On 5/13/14, 11:21 AM, Rasoul Nasiri wrote: Why do you think positive values are doubtful? What matters are energy differences
[gmx-users] About cutoffs for OPLSAA
Dear all, I want to study a system consisting of water, methane and CO2. I will use TIP4P/Ice , OPLSAA and EPM2 models respectively for the components . My problem is, how do I make a choice of the cutoff lengths, for a system which uses different molecular models which were parametrized individually at different cutoff values. I am not sure whether I should use the LJ and Coulomb cutoff values to the OPLSAA forcefield which I would be using. The paper about EPM2 (J.Phys.Chem. Vol. 99, N0.31, 1995) mentions cutoff values of 1nm, whereas the article (J. Chem. Phys. 122, 234511,2005) about TIP4P/Ice model mentions an LJ cutoff value of 0.85nm. To add to my confusion, in Justin's manual on free energy calculation which used OPLSAA forcefield with TIP3P water model, cutoff lengths 1.0 and 0.9nm are used for coulombs and vdW interactions, whereas, article (J. Phys. Chem. B 2012, 116, 14115-14125) which uses OPLSAA again but with a TIP4P-Ew model for water uses slghtly different values (rvdw=0.95 and coulomb = 0.85nm) . Please let me know what you think. Regards, Sujith. -- 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] about tutorial on virtual sites
Hello Justin, I found that in the virtual sites section of your tutorial, the energy values of the system large positive . You had warned in the topology file that the choice of the partial charges used are not guaranteed to work always with OPLSAA . However, the momentum of inertial values points indicate proper dynamics of the system. Is the positive potential energy an out come of using the OPLSAA forcefield ? Is it normal to have positive potential energy? Regards, Sujith. -- 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] Reproducibility in free energy calculations.
Hello Justin, I tried the free energy calculations in your tutorial at different sets of conditions. Also I repeated a simulation at the same conditions with same parameters, thrice on the same computer on same number of processors. The results in the three cases where different (4.16 +/- 0.19 , 4.36 +/-0.14 , 4.54 +/-0.15 kJ/mol ) . I am not surprised, since I thought that this might happen to a small system (241 water + 1 methane). However, I want to be sure that this is OK. I guess in cases like this I should go for an average value of many runs. Regards, Sujith. -- 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] Dubious results with NPT
Hello, This time I tried sesiisotropic scaling with Berendsen Barostat. The compressibility was set to zero in Z-direction allowing scaling only in X and Y. The NPT simulation was run for 10ns. The result did not look good and the average pressure was ~36 bar (reference pressure is 1 bar). Even more confusing was the graphs of various properties against time. Pressure vs time shows a sudden decrease in fluctuations at about 5ns and thereafter . The pressure fluctuations lowered to around ~170bar from ~650 bar. A dramatic change was observed in many properties at around this time. For eg, the short range LJ interaction was found to decrease sharply at ~5ns. Also, total energy, Kinetic, Potential energies and the temperature showed a rapid drop in the magnitude of their fluctuations from the same time. Nothing notable was observed in the system trajectory at around this time. I am totally confused, and I feel that I made some very bad choice again the MDP file given below. Please let me know what you think. Suggest something on this issue which I should read about. ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.001 nsteps = 1000 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 25000 nstvout = 25000 nstfout = 25000 nstlog = 100 nstenergy = 100 nstxtcout = 100 xtc_precision = 1000 xtc_grps = System energygrps = System ; 7.3.9 Neighbor Searching nstlist= 10 ns_type = grid pbc = xyz rlist = 0.9 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.9 ; 7.3.11 VdW vdwtype = cut-off rvdw = 1.4 DispCorr = EnerPres ; 7.3.13 Ewald fourierspacing = 0.12 pme_order = 4 ewald_rtol= 1e-5 ; 7.3.14 Temperature Coupling tcoupl = nose-hoover tc_grps = CHXSOL tau_t = 0.50.5 ref_t = 298298 ; 7.3.15 Pressure Coupling pcoupl = berendsen pcoupltype = semiisotropic tau_p = 2.0 compressibility= 4.5e-5 0.0 ref_p = 1.01.0 ; 7.3.17 Velocity Generation gen_vel = no gen_temp= 298 gen_seed = -1 ; 7.3.18 Bonds constraints = none constraint_algorithm= LINCS continuation= no lincs_order = 4 lincs_iter = 1 lincs_warnangle = 30 Regards, Sujith. On Wed, Mar 12, 2014 at 5:28 PM, Justin Lemkul jalem...@vt.edu wrote: On 3/12/14, 7:51 AM, sujithkakkat . wrote: Hello, After a while I got back to the problem posted here. This issue was the large value for average pressure(~25 bar against the reference pressure of 1 bar) in NPT simulations with parrinello rahman barostat. The system studied is cyclohexane-water system with an interface. The forcefield is Gromos96. Gromacs 4.6.5 version is used. I recollect a post from Michael Shirts which says that with systems that are heterogeneous in direction (not uniform in x y and z), any errors in the pressure may be magnified. Is this what is happening in my case? Based on your suggestions , I tried the following to analyze/solve the problem, and till now the results have not improved a bit. (i) The system was simplified. Each phases was simulated separately. The results were good, with proper values of average pressure after NPT simulation for both independent water and cyclohexane systems. However, when doing the same with similar parameters on the water-cyclohexane combined system with parrinello-rahman barostat , still gives very high average pressure (even after 10ns ) (ii) Changing the tau_p values did not help either. (I tried 1ps, 2ps, and 4ps) My question is whether, I can try the semiisotropic scaling? The manual says that it is useful in case of systems with interface. I am not sure, since none of the replies I got suggested that option, which gives the impression that it is not the best solution. But , to my limited knowledge that is the only option left to be tried to solve the issue. Worth a shot. Also try changing the barostat to Berendsen to see if that changes anything. There are a lot of interrelated issues here, possibly a problem with the barostat, the tau_p value, the type of coupling, the combination of the integrator and barostat...you see how complicated it gets, and that's only some of the possible issues. -Justin The mdp file used is, ; 7.3.3 Run Control integrator = md
Re: [gmx-users] Dubious results with NPT
Hello Mark, I guess you were asking whether I ran the simulation as a continuation to a previous 5ns run? But the results I got are from a continuous 10ns simulation. Sujith. On Fri, Mar 14, 2014 at 4:28 PM, Mark Abraham mark.j.abra...@gmail.comwrote: The simplest explanation would be that you've appended to a previous 5ns trajectory, not run a new trajectory. Check the .log file and the length of time you expected this job to run (wall and simulation). Mark On Fri, Mar 14, 2014 at 9:02 AM, sujithkakkat . sujithk...@gmail.com wrote: Hello, This time I tried sesiisotropic scaling with Berendsen Barostat. The compressibility was set to zero in Z-direction allowing scaling only in X and Y. The NPT simulation was run for 10ns. The result did not look good and the average pressure was ~36 bar (reference pressure is 1 bar). Even more confusing was the graphs of various properties against time. Pressure vs time shows a sudden decrease in fluctuations at about 5ns and thereafter . The pressure fluctuations lowered to around ~170bar from ~650 bar. A dramatic change was observed in many properties at around this time. For eg, the short range LJ interaction was found to decrease sharply at ~5ns. Also, total energy, Kinetic, Potential energies and the temperature showed a rapid drop in the magnitude of their fluctuations from the same time. Nothing notable was observed in the system trajectory at around this time. I am totally confused, and I feel that I made some very bad choice again the MDP file given below. Please let me know what you think. Suggest something on this issue which I should read about. ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.001 nsteps = 1000 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 25000 nstvout = 25000 nstfout = 25000 nstlog = 100 nstenergy = 100 nstxtcout = 100 xtc_precision = 1000 xtc_grps = System energygrps = System ; 7.3.9 Neighbor Searching nstlist= 10 ns_type = grid pbc = xyz rlist = 0.9 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.9 ; 7.3.11 VdW vdwtype = cut-off rvdw = 1.4 DispCorr = EnerPres ; 7.3.13 Ewald fourierspacing = 0.12 pme_order = 4 ewald_rtol= 1e-5 ; 7.3.14 Temperature Coupling tcoupl = nose-hoover tc_grps = CHXSOL tau_t = 0.50.5 ref_t = 298298 ; 7.3.15 Pressure Coupling pcoupl = berendsen pcoupltype = semiisotropic tau_p = 2.0 compressibility= 4.5e-5 0.0 ref_p = 1.01.0 ; 7.3.17 Velocity Generation gen_vel = no gen_temp= 298 gen_seed = -1 ; 7.3.18 Bonds constraints = none constraint_algorithm= LINCS continuation= no lincs_order = 4 lincs_iter = 1 lincs_warnangle = 30 Regards, Sujith. On Wed, Mar 12, 2014 at 5:28 PM, Justin Lemkul jalem...@vt.edu wrote: On 3/12/14, 7:51 AM, sujithkakkat . wrote: Hello, After a while I got back to the problem posted here. This issue was the large value for average pressure(~25 bar against the reference pressure of 1 bar) in NPT simulations with parrinello rahman barostat. The system studied is cyclohexane-water system with an interface. The forcefield is Gromos96. Gromacs 4.6.5 version is used. I recollect a post from Michael Shirts which says that with systems that are heterogeneous in direction (not uniform in x y and z), any errors in the pressure may be magnified. Is this what is happening in my case? Based on your suggestions , I tried the following to analyze/solve the problem, and till now the results have not improved a bit. (i) The system was simplified. Each phases was simulated separately. The results were good, with proper values of average pressure after NPT simulation for both independent water and cyclohexane systems. However, when doing the same with similar parameters on the water-cyclohexane combined system with parrinello-rahman barostat , still gives very high average pressure (even after 10ns ) (ii) Changing the tau_p values did not help either. (I tried 1ps, 2ps, and 4ps
Re: [gmx-users] Dubious results with NPT
Hello, After a while I got back to the problem posted here. This issue was the large value for average pressure(~25 bar against the reference pressure of 1 bar) in NPT simulations with parrinello rahman barostat. The system studied is cyclohexane-water system with an interface. The forcefield is Gromos96. Gromacs 4.6.5 version is used. I recollect a post from Michael Shirts which says that with systems that are heterogeneous in direction (not uniform in x y and z), any errors in the pressure may be magnified. Is this what is happening in my case? Based on your suggestions , I tried the following to analyze/solve the problem, and till now the results have not improved a bit. (i) The system was simplified. Each phases was simulated separately. The results were good, with proper values of average pressure after NPT simulation for both independent water and cyclohexane systems. However, when doing the same with similar parameters on the water-cyclohexane combined system with parrinello-rahman barostat , still gives very high average pressure (even after 10ns ) (ii) Changing the tau_p values did not help either. (I tried 1ps, 2ps, and 4ps) My question is whether, I can try the semiisotropic scaling? The manual says that it is useful in case of systems with interface. I am not sure, since none of the replies I got suggested that option, which gives the impression that it is not the best solution. But , to my limited knowledge that is the only option left to be tried to solve the issue. The mdp file used is, ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.001 nsteps = 1000 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 25000 nstvout = 25000 nstfout = 25000 nstlog = 100 nstenergy = 100 nstxtcout = 100 xtc_precision = 1000 xtc_grps = System energygrps = System ; 7.3.9 Neighbor Searching nstlist= 10 ns_type = grid pbc = xyz rlist = 0.9 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.9 ; 7.3.11 VdW vdwtype = cut-off rvdw = 1.4 DispCorr = EnerPres ; 7.3.13 Ewald fourierspacing = 0.12 pme_order = 4 ewald_rtol= 1e-5 ; 7.3.14 Temperature Coupling tcoupl = nose-hoover tc_grps = CHXSOL tau_t = 0.50.5 ref_t = 298298 ; 7.3.15 Pressure Coupling pcoupl = parrinello-rahman pcoupltype = isotropic tau_p = 2.0 compressibility= 4.5e-5 ref_p = 1.0 ; 7.3.17 Velocity Generation gen_vel = no gen_temp= 298 gen_seed = -1 ; 7.3.18 Bonds constraints = none constraint_algorithm= LINCS continuation= no lincs_order = 4 lincs_iter = 1 lincs_warnangle = 30 Please comment. Sujith. -- 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] A question about blowing up
Hello, I am following justin lemkul's tutorials for the cyclohexane-water biphasic system using Gromos96 forcefield. Things went well till NVT equilibration where I had this Blowing up issue. Shortening the time step alone did not help, and it was difficult for me to figure out what was happening from watching the trajectory. g_energy data indicated a spike on the, short range LJ contribution at around 10ps. Then I switched off the constraints (but not for water, SETTLE still used) and ran the simulation again with a small time step. Things worked well this time , the system reaching the target temperature. I want to make sure that switching off the constraints is a reasonable thing to solve Blowing up and that it would not affect in any way the accuracy of further simulations. http://www.gromacs.org/Documentation/Terminology/Blowing_Up says that constraints are exactly not the source of the Blowing up problem. However, http://www.gromacs.org/Documentation/How-tos/Steps_to_Perform_a_Simulation suggests removing constraints as a solution. So, in general, is the problem solved by removing constrains, or am I cheating? Please comment. Here is the .mdp file which solved the Blowing up issue. ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.0002 nsteps = 250 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 25000 nstvout = 25000 nstfout = 25000 nstlog = 100 nstenergy = 100 nstxtcout = 100 xtc_precision = 1000 xtc_grps = System energygrps = System ; 7.3.9 Neighbor Searching nstlist= 1 ns_type = grid pbc = xyz rlist = 0.9 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.9 ; 7.3.11 VdW vdwtype = cut-off rvdw = 1.4 DispCorr = EnerPres ; 7.3.13 Ewald fourierspacing = 0.12 pme_order = 4 ewald_rtol= 1e-5 ; 7.3.14 Temperature Coupling tcoupl = berendsen tc_grps = CHXSOL tau_t = 0.10.1 ref_t = 298298 ; 7.3.17 Velocity Generation gen_vel = yes gen_temp= 298 gen_seed = -1 ; 7.3.18 Bonds constraints = none constraint_algorithm= LINCS continuation= no lincs_order = 4 lincs_iter = 1 lincs_warnangle = 30 ~ Thanks, Sujith. -- 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] Severe error with NPT
Dear GROMACS users, I am new to GROMACS, and recently started using the version 4.6.5. I have seen a lot of NPT related issues raised earlier in this forum, but in my case the error looks much more severe. I am following Justin Lemkul's tutorial ( http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/biphasic/01_genconf.html) on Biphasic system, containing cyclohexane and water.Everything went well till the NPT simulation to bring system to reference pressure of 1 bar. Here are the details of the system , .mdp file , what I did and where I find the problem. SYSTEM : 512 cyclohexane + 2720 water molecules. CURRENT STATUS: (1)Energy minimization : energy converged. (2)NVT (tcoupl=berendsen , tau_t=0.1ps / ref_t=298K) completed and target pressure of 298 K attained. (3)NPT (tcoupl=nose-hoover, tau_t=0.5ps / pcoupl=berendsen, tau_p=1ps / ref_p = 1bar , total time=15 ns) completed with the following result: Energy Average Err.Est. RMSD Tot-Drift --- Pressure1.08924 0.67 114.66 -1.74033 (bar) Since berendsen barostat doesn't generate the true ensemble, an equilibration with pcoupl=parrinello-rahman is performed in the next step. (4) NPT (pcoupl=parrinello-rahman, tau_p=5ps, total time=5ns) PROBLEM FACED: The average pressure too high as shown below which I feel is not going to improve. Energy Average Err.Est. RMSD Tot-Drift --- Pressure23.6651 0.53144.464 -1.00634 (bar) I am aware of the fact that pressure is subject to large fluctuations in small sized systems and that this may affect the average value of the pressure. But , here the average pressure looks too large to be ignored. The pressure-vs-time graph doesn't show any upward trend, and the pressure looks like fluctuating about a central value. Here is the .mdp file. Only changes made from the previous .mdp for berendsen pressure coupling , are with pcoupl and tau_p. ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.002 nsteps = 250 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 2500 nstvout = 2500 nstfout = 2500 nstlog = 2500 nstenergy = 100 nstxtcout = 1000 xtc_precision = 1000 xtc_grps= System energygrps = System ; 7.3.9 Neighbor Searching nstlist = 1 ns_type= grid pbc = xyz rlist = 0.8 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.8 ; 7.3.11 VdW vdwtype = cut-off rvdw = 0.8 DispCorr= EnerPres ; 7.3.13 Ewald fourierspacing = 0.12 pme_order = 4 ewald_rtol = 1e-5 ; 7.3.14 Temperature Coupling tcoupl = nose-hoover tc_grps = CHXSOL tau_t= 0.50.5 ref_t= 298298 ; 7.3.15 Pressure Coupling pcoupl= parrinello-rahman pcoupltype = isotropic tau_p = 5.0 compressibility= 4.5e-5 ref_p = 1.0 ; 7.3.17 Velocity Generation gen_vel = no ; 7.3.18 Bonds constraints = all-bonds constraint_algorithm= LINCS continuation = yes lincs_order = 4 lincs_iter= 1 lincs_warnangle = 30 npt.mdp 63L, 4080C I guess there is something seriously wrong in the choice of methods/parameters in the .mdp file, which I cant figure out. Kindly go through and let me know your comments. I would be happy to give any further details. Any help would be appreciated. Regards, Sujith 107,4 Bot gromacs.org_gmx-users@maillist.sys.kth.se -- 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] Dubious results with NPT
Hello, Thank you both for the comments. I am using gromos96 forcefield . I read a little bit and as you said the nonbonded cutoff has to be higher. The tau_p=5ps was chosen , since the manual mentions that the value has to be raised by 4-5 times on going from berendsen to parrinello-rahman barostat, though I did not completely follow the reasons behind it. I will try with lower values. I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. I guess the sigma for the case of water was low and therefore the small cut-off of 0.8nm did not matter. However the case of cyclohexane alone remains to be tried. I guess Dr Vitaly was saying about using a switch/shift function. I will try the simulation with the new settings and see. Regards, Sujith. On Sun, Feb 23, 2014 at 10:37 PM, Dr. Vitaly Chaban vvcha...@gmail.comwrote: There is such thing in simulations as energy conservation... If you use vdwtype= cut-off and this cut-off happens at 0.8nm, while sigma for the largest atom is ~0.34nm, the problems are inevitable. Your cut-off should not be smaller than 0.90nm, and you need to apply a more polite method to bring pairwise energy term down to zero at this cut-off. Dr. Vitaly V. Chaban On Sun, Feb 23, 2014 at 3:07 PM, Justin Lemkul jalem...@vt.edu wrote: On 2/23/14, 8:30 AM, sujith wrote: Dear GROMACS users, I am new to GROMACS, and recently started using the version 4.6.5. I have seen a lot of NPT related issues raised earlier in this forum, but in my case the error looks much more severe. I am following Justin Lemkul's tutorial ( http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/biphasic/01_genconf.html ) on Biphasic system, containing cyclohexane and water.Everything went well till the NPT simulation to bring system to reference pressure of 1 bar. Here are the details of the system , .mdp file , what I did and where I find the problem. SYSTEM : 512 cyclohexane + 2720 water molecules. CURRENT STATUS: (1)Energy minimization : energy converged. (2)NVT (tcoupl=berendsen , tau_t=0.1ps / ref_t=298K) completed and target pressure of 298 K attained. (3)NPT (tcoupl=nose-hoover, tau_t=0.5ps / pcoupl=berendsen, tau_p=1ps / ref_p = 1bar , total time=15 ns) completed with the following result: Energy Average Err.Est. RMSD Tot-Drift --- Pressure1.08924 0.67 114.66 -1.74033 (bar) Since berendsen barostat doesn't generate the true ensemble, an equilibration with pcoupl=parrinello-rahman is performed in the next step. (4) NPT (pcoupl=parrinello-rahman, tau_p=5ps, total time=5ns) PROBLEM FACED: The average pressure too high as shown below which I feel is not going to improve. Energy Average Err.Est. RMSD Tot-Drift --- Pressure23.6651 0.53144.464 -1.00634 (bar) I am aware of the fact that pressure is subject to large fluctuations in small sized systems and that this may affect the average value of the pressure. But , here the average pressure looks too large to be ignored. The pressure-vs-time graph doesn't show any upward trend, and the pressure looks like fluctuating about a central value. Here is the .mdp file. Only changes made from the previous .mdp for berendsen pressure coupling , are with pcoupl and tau_p. ; 7.3.3 Run Control integrator = md tinit = 0 dt = 0.002 nsteps = 250 comm_mode = Linear nstcomm = 1 comm_grps = CHX SOL ; 7.3.8 Output Control nstxout = 2500 nstvout = 2500 nstfout = 2500 nstlog = 2500 nstenergy = 100 nstxtcout = 1000 xtc_precision = 1000 xtc_grps= System energygrps = System ; 7.3.9 Neighbor Searching nstlist = 1 ns_type= grid pbc = xyz rlist = 0.8 ; 7.3.10 Electrostatics coulombtype = PME rcoulomb = 0.8 ; 7.3.11 VdW vdwtype = cut-off rvdw = 0.8 What force field are you using? If it's Gromos96 like my tutorial, the value of rvdw is much too short and can lead to artifacts. Using 0.8 for
