Re: [gmx-users] Pressure coupling
Probably. See manual 7.3.15. Semiisotropic with some zero compressibilities is probably what you want. Mark On Sun, Jun 2, 2013 at 9:33 PM, Marcelo Vanean vanea...@gmail.com wrote: Hello everyone. I want to know if can be applied pressure coupling only in the z direction, allowing the edges x and y simulation box with fixed size. -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Pressure coupling doubt
On 29/03/2012 5:38 PM, bipin singh wrote: Hello, I have two doubts regarding pressure coupling in Gromacs: 1) When I use pcoupl=no the mdp.out shows the following ; Pressure coupling pcoupl = no Pcoupltype = Isotropic nstpcouple = -1 ; Time constant (ps), compressibility (1/bar) and reference P (bar) tau-p= 1 compressibility = ref-p= I have not used the pcoupl(=no) then why it is showing Pcoupltype=Isotropic. ... because that was either the value you input, or the default value. Neither matters, because you said not to use pressure coupling. (2) This is a silly question: What will happen if we use following option in mdp of NVT simulation: pcoupl=no compressibility= 4.5e-5 does it will affect the NVT criteria ? or It will ignore the compressibility input? It will ignore it. Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Pressure coupling and membrane-type simulations
Andrew DeYoung wrote: Hi, I am interested in doing a membrane-type simulation, in which I have all-atom membrane walls parallel to xy plane, at z = -z_0 and z = +z_0 (where z_0 is a constant). I would like to run an NPT simulation at 1 atm. What type of pressure coupling should I use? Isotropic pressure coupling requires only one input value (either compressibility or reference pressure ref_p), whereas I think that semiisotropic, anisotropic, and surface-tension pressure coupling require specification of both the compressibility and ref_p (tensors). Clearly, I should not use isotropic pressure coupling, because clearly my system is not isotropic. However, what if I do not know and cannot find in the literature the compressibility of the liquid system that I am placing between the membrane walls? If I do not know the compressibility very accurately or at all, then it seems that I cannot use semiisotropic, anisotropic, or surface-tension pressure coupling. If you have time, I would like to ask an additional question. Now suppose I know the compressibility of the liquid between the membrane. Now what pressure coupling type should I use; should I use semiisotropic, anisotropic, or surface-tension pressure coupling? Both semiisotropic and surface-tension look reasonable. In the manual (http://manual.gromacs.org/current/online/mdp_opt.html#pc), semiisotropic pressure coupling is useful for systems that are isotropic in x and y, but different in z (which is the situation I have here). Surface-tension also looks like it describes a similar situation, but it requires the specification of the surface tension of the liquid, which I do not know. I am sorry that my questions are quite vague. If you have time, do you have any general thoughts? Or can you please recommend any papers that would help me understand and choose between the pressure coupling types? I doubt I can really answer much of this, but isn't the most pressing consideration the walls themselves? If they're reasonably rigid, the compressibility of the fluid layer is largely irrelevant, isn't it? Then again, if the walls are not rigid, then they become a liability under pressure coupling as they may buckle. Seems to me this should be the principal concern. -Justin -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Pressure coupling and membrane-type simulations
On 2012-02-21 03:41:07PM -0500, Justin A. Lemkul wrote: Andrew DeYoung wrote: Hi, I am interested in doing a membrane-type simulation, in which I have all-atom membrane walls parallel to xy plane, at z = -z_0 and z = +z_0 (where z_0 is a constant). I would like to run an NPT simulation at 1 atm. What type of pressure coupling should I use? Isotropic pressure coupling requires only one input value (either compressibility or reference pressure ref_p), whereas I think that semiisotropic, anisotropic, and surface-tension pressure coupling require specification of both the compressibility and ref_p (tensors). Clearly, I should not use isotropic pressure coupling, because clearly my system is not isotropic. However, what if I do not know and cannot find in the literature the compressibility of the liquid system that I am placing between the membrane walls? If I do not know the compressibility very accurately or at all, then it seems that I cannot use semiisotropic, anisotropic, or surface-tension pressure coupling. If you have time, I would like to ask an additional question. Now suppose I know the compressibility of the liquid between the membrane. Now what pressure coupling type should I use; should I use semiisotropic, anisotropic, or surface-tension pressure coupling? Both semiisotropic and surface-tension look reasonable. In the manual (http://manual.gromacs.org/current/online/mdp_opt.html#pc), semiisotropic pressure coupling is useful for systems that are isotropic in x and y, but different in z (which is the situation I have here). Surface-tension also looks like it describes a similar situation, but it requires the specification of the surface tension of the liquid, which I do not know. I am sorry that my questions are quite vague. If you have time, do you have any general thoughts? Or can you please recommend any papers that would help me understand and choose between the pressure coupling types? I doubt I can really answer much of this, but isn't the most pressing consideration the walls themselves? If they're reasonably rigid, the compressibility of the fluid layer is largely irrelevant, isn't it? Then again, if the walls are not rigid, then they become a liability under pressure coupling as they may buckle. Seems to me this should be the principal concern. -Justin Another consideration is that the documentation states when using built-in implicit walls at z=0 and z=zbox to use anistropic pressure coupling with the x/y compressibility set to 0 otherwise the surface area will change. Maybe start with that and go from there? -- == Peter C. Lai| University of Alabama-Birmingham Programmer/Analyst | KAUL 752A Genetics, Div. of Research | 705 South 20th Street p...@uab.edu| Birmingham AL 35294-4461 (205) 690-0808 | == -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] pressure coupling
On 17/12/2011 4:58 PM, mohammad agha wrote: Dear GROMACS Users I have a warning after doing pr.mdp as followed, I read errors in GROMACS site and checked mailing list, but my warning is only step1 and after that equilibration is run normally till end. 500 steps, 15.0 ps. step 0 Step 1 Warning: pressure scaling more than 1%, mu: 1.04679 1.04679 1.04679 Step 1 Warning: pressure scaling more than 1%, mu: 1.04679 1.04679 1.04679 my pr.mdp about pressure coupling is: pcoupl= berendsen pcoupltype= isotropic; uniform scaling of box vectors tau_p= 4 ; time constant, in ps ref_p= 1.0 ; reference pressure, in bar compressibility = 4.5e-5 ;isothermal compressibility of water, bar^-1 I work in Martini Coarse-Grained and my time step is 0.03 ps. I increased tau_p but my warning increased to step 1 and step 11, I think this warning since it is only in initial steps isn't important, is my idea right? That kind of thing is normal in the very early stages of NPT equilibration when the volume is not quite right. mdrun is just warning you that there are some changes going on. Use g_energy to see those. If it bothers you, then you will want to adjust your density (i.e. number of water molecules) before starting. Mark -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] pressure coupling
Thank you very much. Yes, I checked NPT run and I saw the size and volume of box had been increased. So I increased size of box to 12. My system consists of: 151 lipids+151 ions+16200 water and antifreezawater molecules May I know your idea about this, please? Best Regards Sara From: Mark Abraham mark.abra...@anu.edu.au To: Discussion list for GROMACS users gmx-users@gromacs.org Sent: Saturday, December 17, 2011 10:22 AM Subject: Re: [gmx-users] pressure coupling On 17/12/2011 5:40 PM, mohammad agha wrote: Thank you very much. Excuse me, I did another thing, may I know is it right, Please? I did editconf before equilibration (pr.mdp) and increased my box vector from 11.2 to 12: editconf -f em3.gro -o out.gro -c -d 1.0 -bt cubic -box 12 12 12 -center 6 6 6 and I minimized my system (out.gro) and then did equilibration (pr.mdp). My system gave me warning for box sizes 11.4, 11.6 and 11.8, consequently I select 12 and pr.mdp is run without warning. Is it right? We don't know. You've been changing your box volume but haven't said anything about the contents. You need to look at how the box changes over the NPT run, like I said last time. If your generated solvent density was wrong such that you needed to add about 20% more volume, then expect to need to equilibrate for quite a while. Mark Best Regards Sara From: Mark Abraham mark.abra...@anu.edu.au To: Discussion list for GROMACS users gmx-users@gromacs.org Sent: Saturday, December 17, 2011 9:31 AM Subject: Re: [gmx-users] pressure coupling On 17/12/2011 4:58 PM, mohammad agha wrote: Dear GROMACS Users I have a warning after doing pr.mdp as followed, I read errors in GROMACS site and checked mailing list, but my warning is only step1 and after that equilibration is run normally till end. 500 steps, 15.0 ps. step 0 Step 1 Warning: pressure scaling more than 1%, mu: 1.04679 1.04679 1.04679 Step 1 Warning: pressure scaling more than 1%, mu: 1.04679 1.04679 1.04679 my pr.mdp about pressure coupling is: pcoupl = berendsen pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 4 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5 ;isothermal compressibility of water, bar^-1 I work in Martini Coarse-Grained and my time step is 0.03 ps. I increased tau_p but my warning increased to step 1 and step 11, I think this warning since it is only in initial steps isn't important, is my idea right? That kind of thing is normal in the very early stages of NPT equilibration when the volume is not quite right. mdrun is just warning you that there are some changes going on. Use g_energy to see those. If it bothers you, then you will want to adjust your density (i.e. number of water molecules) before starting. Mark -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists -- gmx-users mailing list gmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists-- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Pressure coupling problem
Fabian Casteblanco wrote: Hi, I'm still in my first few months of using Gromacs. I started by creating an *.itp and *.top file for /Ethanol/ using CHARMM force field parameters. I made the molecule and it looked fine, put 1000 molecules in a box, energy minimized it to a negative potential energy, viewed it on VMD, again looks fine. When I started running the NVT script, I set it equal to a ref_T of 298 K. It equilibrated at the temperature. Then I tried using an NPT script to equilibrate it to a ref_p of 1 bar. This is where I get the problem. The output shows the density is close to the actual experimental value of 0.789 g/cm^3. But for some reason, my pressure never gets an average of 1 bar. It keeps oscillating, which I understand is normal, but the average is always 1.3 or 1.4 bar (it seems the longer I let it run, the larger the average pressure; 1.38 for 50,000 steps,dt=0.002 and 1.45 for 75,000 steps,dt=0.002). I don't understand why the ref_p of 1 bar is not working when I run this NPT.mdp script file. My simple goal is to have 1000 molecules of ethanol using CHARMM ff parameters at 25degC and 1 bar and somewhere near the experimental density. Your equilibration period (100-150 ps) is rather short, and the systematic increase suggests that you're simply not equilibrated yet. Also bear in mind that a quantity that is prone to pressure fluctuations in the hundreds to thousands can only be so accurate. There was a very thorough discussion about the statistical significance of pressure values that are not equal to ref_p just some time ago. You may want to look through the archives to find this discussion. -Justin I would really appreciate anybody's help! I'm new to this but I'm eager to keep getting better. Thanks. _NVT SCRIPT (this works fine and takes me to 298 K)_ File Edit Options Buffers Tools Help title =CHARMM ETHANOL NVT equilibration ;define =-DPOSRES ;position restrain the protein ;Run parameters integrator =md ;leap-frog algorithm nsteps =5 ;2 * 5 = 100 ps dt =0.002 ;2fs ;Output control nstxout =100;save coordinates every 0.2 ps nstvout =100;save velocities every 0.2 ps nstenergy =100;save energies every 0.2 ps nstlog =100;update log file every 0.2 ps ;Bond parameters continuation=no ;first dynamics run constraint_algorithm=lincs ;holonomic constraints constraints =all-bonds ;all bonds (even heavy atom-H bonds)constraind lincs_iter =1 ;accuracy of LINCS lincs_order =4 ;also related to accuracy ;Neighborhood searching ns_type =grid ;search neighboring grid cells nstlist =5 ;10 fs rlist =1.0;short-range neighborlist cutoff (in nm) rcoulomb=1.0;short-range electrostatic cutoff (in nm) rvdw=1.0;short-range van der Waals cutoff (in nm) ;Electrostatics coulombtype =PME;Particle Mesh Ewald for long-range electrostat\ ;ics pme_order =4 ;cubic interpolation fourierspacing =0.16 ;grid spacing for FFT ;Temperature coupling is on tcoupl =V-rescale ;modified Berendsen thermostat tc_grps =SYSTEM ;two coupling groups - more accurate tau_t =0.1;0.1 ;time constant, in ps ref_t =298;25 ;reference temperature, one for each \ ;group, in K ;Pressure coupling is off pcoupl =no ;no pressure coupling in NVT ;Periodic boundary conditions pbc =xyz; 3-D PBC ;Dispersion correction DispCorr=EnerPres ;account for cut-off vdW scheme ;Velocity generation gen_vel =yes;assign velocities from Maxwell distribution gen_temp=25 ;temperature for Maxwell distribution gen_seed=-1 ;generate a random seed ;END _NPT SCRIPT_ File Edit Options Buffers Tools Help title =Ethanol npt equilibration ;define =-DPOSRES ;position restrain the protein ;Run parameters integrator =md ;leap-frog algorithm nsteps =5 ;2 * 5 = 100 ps dt =0.002 ;2fs ;Output control nstxout =100;save coordinates every 0.2 ps nstvout =100;save velocities every 0.2 ps nstenergy =100;save energies every 0.2 ps nstlog =100;update log file every 0.2 ps ;Bond parameters continuation=yes;Restarting after NVT constraint_algorithm=lincs ;holonomic constraints constraints =all-bonds ;all bonds (even heavy atom-H bonds)constraind lincs_iter =1 ;accuracy of LINCS lincs_order =4 ;also related to accuracy
Re: [gmx-users] Pressure coupling problem
So your density graph looks stabilized? I also tend to look for changes in box x, y, z as well since the scale of their changes is easier to track. Sometimes it helps to look at the error vs. rmsd vs total drift statistics as well for such parameters that are easier to track - again if density shows stability for 0-150ps then check your boxes, it might be shrinking or growing due to the pressure perturbation, and you can use the average rate of change of those to find your equilibration point instead of trying to do something with a -300 to 300 bar pressure smear or whatnot. On 2011-04-11 04:17:59PM -0500, Justin A. Lemkul wrote: Fabian Casteblanco wrote: Hi, I'm still in my first few months of using Gromacs. I started by creating an *.itp and *.top file for /Ethanol/ using CHARMM force field parameters. I made the molecule and it looked fine, put 1000 molecules in a box, energy minimized it to a negative potential energy, viewed it on VMD, again looks fine. When I started running the NVT script, I set it equal to a ref_T of 298 K. It equilibrated at the temperature. Then I tried using an NPT script to equilibrate it to a ref_p of 1 bar. This is where I get the problem. The output shows the density is close to the actual experimental value of 0.789 g/cm^3. But for some reason, my pressure never gets an average of 1 bar. It keeps oscillating, which I understand is normal, but the average is always 1.3 or 1.4 bar (it seems the longer I let it run, the larger the average pressure; 1.38 for 50,000 steps,dt=0.002 and 1.45 for 75,000 steps,dt=0.002). I don't understand why the ref_p of 1 bar is not working when I run this NPT.mdp script file. My simple goal is to have 1000 molecules of ethanol using CHARMM ff parameters at 25degC and 1 bar and somewhere near the experimental density. Your equilibration period (100-150 ps) is rather short, and the systematic increase suggests that you're simply not equilibrated yet. Also bear in mind that a quantity that is prone to pressure fluctuations in the hundreds to thousands can only be so accurate. There was a very thorough discussion about the statistical significance of pressure values that are not equal to ref_p just some time ago. You may want to look through the archives to find this discussion. -Justin I would really appreciate anybody's help! I'm new to this but I'm eager to keep getting better. Thanks. _NVT SCRIPT (this works fine and takes me to 298 K)_ File Edit Options Buffers Tools Help title =CHARMM ETHANOL NVT equilibration ;define =-DPOSRES ;position restrain the protein ;Run parameters integrator =md ;leap-frog algorithm nsteps =5 ;2 * 5 = 100 ps dt =0.002 ;2fs ;Output control nstxout =100;save coordinates every 0.2 ps nstvout =100;save velocities every 0.2 ps nstenergy =100;save energies every 0.2 ps nstlog =100;update log file every 0.2 ps ;Bond parameters continuation=no ;first dynamics run constraint_algorithm=lincs ;holonomic constraints constraints =all-bonds ;all bonds (even heavy atom-H bonds)constraind lincs_iter =1 ;accuracy of LINCS lincs_order =4 ;also related to accuracy ;Neighborhood searching ns_type =grid ;search neighboring grid cells nstlist =5 ;10 fs rlist =1.0;short-range neighborlist cutoff (in nm) rcoulomb=1.0;short-range electrostatic cutoff (in nm) rvdw=1.0;short-range van der Waals cutoff (in nm) ;Electrostatics coulombtype =PME;Particle Mesh Ewald for long-range electrostat\ ;ics pme_order =4 ;cubic interpolation fourierspacing =0.16 ;grid spacing for FFT ;Temperature coupling is on tcoupl =V-rescale ;modified Berendsen thermostat tc_grps =SYSTEM ;two coupling groups - more accurate tau_t =0.1;0.1 ;time constant, in ps ref_t =298;25 ;reference temperature, one for each \ ;group, in K ;Pressure coupling is off pcoupl =no ;no pressure coupling in NVT ;Periodic boundary conditions pbc =xyz; 3-D PBC ;Dispersion correction DispCorr=EnerPres ;account for cut-off vdW scheme ;Velocity generation gen_vel =yes;assign velocities from Maxwell distribution gen_temp=25 ;temperature for Maxwell distribution gen_seed=-1 ;generate a random seed ;END _NPT SCRIPT_ File Edit Options Buffers Tools Help title =Ethanol npt
Re: [gmx-users] pressure coupling not enough values (I need 2)
afsaneh maleki wrote: Hi, grompp show the error as below: ERROR: pressure coupling not enough values (I need 2) I used ref_p =1. In file.mdp How to solve this problem? Without seeing the rest of your pressure coupling settings, I can only assume that you're using semi-isotropic coupling, and thus, per the manual: Pressure coupling which is isotropic in the x and y direction, but different in the z direction. This can be useful for membrane simulations. 2 values are needed for x/y and z directions respectively. -Justin thanks in advance, Afsaneh -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
Re: [gmx-users] Pressure coupling
Sai Pooja wrote: Hi, I am aware there have been many threads on this topic and I have looked at some of them. I would still like to make this query. I am running an npt equilibration run for a solvated protein using Charmm27-Tip3p (nocmap). I ran the simulation for 1ns. The average values of pressure converge to 1.1 (ref=1.0), however, the fluctuations are not only large but range from -10^2 --- +10^2. I am pasting an excerpt from my log file here: Seems perfectly normal to me. http://www.gromacs.org/Documentation/Terminology/Pressure#Fluctuation snip Moreover, in this particular version of gromacs (git version, 20th June, 2010), RMS fluctuations are not reported in the log file. This might be nice to restore, but you can still get all of that from g_energy. -Justin Pooja -- Quaerendo Invenietis-Seek and you shall discover. -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Pressure coupling
On 2010-07-15 16.01, Justin A. Lemkul wrote: Sai Pooja wrote: Hi, I am aware there have been many threads on this topic and I have looked at some of them. I would still like to make this query. I am running an npt equilibration run for a solvated protein using Charmm27-Tip3p (nocmap). I ran the simulation for 1ns. The average values of pressure converge to 1.1 (ref=1.0), however, the fluctuations are not only large but range from -10^2 --- +10^2. I am pasting an excerpt from my log file here: Seems perfectly normal to me. http://www.gromacs.org/Documentation/Terminology/Pressure#Fluctuation snip Moreover, in this particular version of gromacs (git version, 20th June, 2010), RMS fluctuations are not reported in the log file. mdrun -nocompact This might be nice to restore, but you can still get all of that from g_energy. -Justin Pooja -- Quaerendo Invenietis-Seek and you shall discover. -- David van der Spoel, Ph.D., Professor of Biology Dept. of Cell Molec. Biol., Uppsala University. Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. sp...@xray.bmc.uu.sehttp://folding.bmc.uu.se -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] pressure coupling
Shuangxing Dai wrote: Hi, all, I am trying to do anisotropic coupling using Parrinello-Rahman. I want to reach the equilibrium state at 1000K, 0.1MPa( 1 bar). However, I find that the fluctuation is so large ( p=80.22 bar, T=1007.23 K, fluctuation is 626 bar for pressure and 11.21 K for temperature.) I run energy minimization first, then Berenderson pressure coupling for 50 ps, finally Parrinello-Rahman pressure coupling for 50 ps. Anyone has experience on anisotropic coupling and know why? My system has just 4000 atoms. This is my .mdp file: For a small system, large fluctations in the pressure are to be expected: http://www.gromacs.org/Documentation/Terminology/Pressure -Justin ;define = -DPOSRES ; RUN CONTROL PARAMETERS = integrator = sd ; start time and timestep in ps = tinit= 0 dt = 0.001 nsteps = 5 ; number of steps for center of mass motion removal = nstcomm = 100 ; OUTPUT CONTROL OPTIONS = ; Output frequency for coords (x), velocities (v) and forces (f) = nstxout = 0 nstvout = 0 nstfout = 0 ; Output frequency for energies to log file and energy file = nstlog = 10 nstenergy= 10 ; Output frequency and precision for xtc file = nstxtcout= 10 xtc-precision= 1000 ; NEIGHBORSEARCHING PARAMETERS = ; nblist update frequency = nstlist = 20 ; ns algorithm (simple or grid) = ns_type = grid ;OPTIONS FOR PRESSURE COUPLING Pcoupl = Parrinello-Rahman pcoupltype = anisotropic tau_p= 5 compressibility = 2.1645e-07 2.1645e-07 2.7322e-07 0 0 0 ref_p= 1 1 1 0 0 0 ;OPTIONS FOR TEMPERATURE COUPLING tc_grps = system tau_t= 1 ref_t= 1000 ; OPTIONS FOR BONDS = constraints = hbonds ; Type of constraint algorithm = constraint-algorithm = Lincs ; Do not constrain the start configuration = unconstrained-start = no ; Relative tolerance of shake = shake-tol= 0.0001 ; Highest order in the expansion of the constraint coupling matrix = lincs-order = 12 ; Lincs will write a warning to the stderr if in one step a bond = ; rotates over more degrees than = lincs-warnangle = 30 ; Periodic boundary conditions: xyz, no, xy pbc = xyz periodic_molecules = no ; nblist cut-off rlist= 1 ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = PME rcoulomb = 1 ; Method for doing Van der Waals vdw-type = Cut-off ; cut-off lengths rvdw = 1 ; Spacing for the PME/PPPM FFT grid fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 ; EWALD/PME/PPPM parameters pme_order= 6 ewald_rtol = 1e-4 ewald_geometry = 3d epsilon_surface = 0 optimize_fft = no I used the same procedure to reach different state ( like 0.01 K, 100K, 300K, 500K, 800K, 1000K). Only the 0.01 K got the correct average pressure and small fluctuation. T p(bar) delta p T (K) delta T 0.011.03E+002.12E+009.95E-031.21E-04 100 5.48359 71.5565 91.198 1.25218 300 2.62807 162.926 290.324 3.36594 500 16.6356 1048.76 492.616 5.45231 800 59.3245 594.773 799.314 9.18876 100080.2257 626.789 1007.23 11.2159 Thanks in advance, Shuangxing Dai -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar MILES-IGERT Trainee Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] pressure coupling
Thanks for the response. Yes, I know that large fluctuations will be expected. But why the average was not even correct? Also, how can I reach some state ( some pressre and temperature) and know/tell that the system is in equilibrium if the average is not correct with large fluctuation? I mean I do not know the criteria of Gromacs. Or maybe my procedure is wrong? Message: 1 Date: Mon, 05 Jul 2010 10:39:30 -0400 From: Justin A. Lemkul jalem...@vt.edu Subject: Re: [gmx-users] pressure coupling To: Discussion list for GROMACS users gmx-users@gromacs.org Message-ID: 4c31eea2.7070...@vt.edu Content-Type: text/plain; charset=UTF-8; format=flowed Shuangxing Dai wrote: Hi, all, I am trying to do anisotropic coupling using Parrinello-Rahman. I want to reach the equilibrium state at 1000K, 0.1MPa( 1 bar). However, I find that the fluctuation is so large ( p=80.22 bar, T=1007.23 K, fluctuation is 626 bar for pressure and 11.21 K for temperature.) I run energy minimization first, then Berenderson pressure coupling for 50 ps, finally Parrinello-Rahman pressure coupling for 50 ps. Anyone has experience on anisotropic coupling and know why? My system has just 4000 atoms. This is my .mdp file: For a small system, large fluctations in the pressure are to be expected: http://www.gromacs.org/Documentation/Terminology/Pressure -Justin ;define = -DPOSRES ; RUN CONTROL PARAMETERS = integrator = sd ; start time and timestep in ps = tinit= 0 dt = 0.001 nsteps = 5 ; number of steps for center of mass motion removal = nstcomm = 100 ; OUTPUT CONTROL OPTIONS = ; Output frequency for coords (x), velocities (v) and forces (f) = nstxout = 0 nstvout = 0 nstfout = 0 ; Output frequency for energies to log file and energy file = nstlog = 10 nstenergy= 10 ; Output frequency and precision for xtc file = nstxtcout= 10 xtc-precision= 1000 ; NEIGHBORSEARCHING PARAMETERS = ; nblist update frequency = nstlist = 20 ; ns algorithm (simple or grid) = ns_type = grid ;OPTIONS FOR PRESSURE COUPLING Pcoupl = Parrinello-Rahman pcoupltype = anisotropic tau_p= 5 compressibility = 2.1645e-07 2.1645e-07 2.7322e-07 0 0 0 ref_p= 1 1 1 0 0 0 ;OPTIONS FOR TEMPERATURE COUPLING tc_grps = system tau_t= 1 ref_t= 1000 ; OPTIONS FOR BONDS = constraints = hbonds ; Type of constraint algorithm = constraint-algorithm = Lincs ; Do not constrain the start configuration = unconstrained-start = no ; Relative tolerance of shake = shake-tol= 0.0001 ; Highest order in the expansion of the constraint coupling matrix = lincs-order = 12 ; Lincs will write a warning to the stderr if in one step a bond = ; rotates over more degrees than = lincs-warnangle = 30 ; Periodic boundary conditions: xyz, no, xy pbc = xyz periodic_molecules = no ; nblist cut-off rlist= 1 ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = PME rcoulomb = 1 ; Method for doing Van der Waals vdw-type = Cut-off ; cut-off lengths rvdw = 1 ; Spacing for the PME/PPPM FFT grid fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 ; EWALD/PME/PPPM parameters pme_order= 6 ewald_rtol = 1e-4 ewald_geometry = 3d epsilon_surface = 0 optimize_fft = no I used the same procedure to reach different state ( like 0.01 K, 100K, 300K, 500K, 800K, 1000K). Only the 0.01 K got the correct average pressure and small fluctuation. T p(bar) delta p T (K) delta T 0.01 1.03E+002.12E+009.95E-031.21E-04 100 5.48359 71.5565 91.198 1.25218 300 2.62807 162.926 290.324 3.36594 500 16.6356 1048.76 492.616 5.45231 800 59.3245 594.773 799.314 9.18876 1000 80.2257 626.789 1007.23 11.2159 Thanks in advance, Shuangxing Dai Thanks, Shuangxing Dai -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists
Re: [gmx-users] pressure coupling
Shuangxing Dai wrote: Thanks for the response. Yes, I know that large fluctuations will be expected. But why the average was not even correct? Also, how can I reach some state ( some pressre and temperature) and know/tell that the system is in equilibrium if the average is not correct with large fluctuation? I mean I do not know the criteria of Gromacs. Or maybe my procedure is wrong? You've only run 100 ps total equilibration. For an extremely high temperature, I might expect it to take quite a bit longer, especially since you're dealing with a system that will be subject to large fluctuations (due to its size). -Justin Message: 1 Date: Mon, 05 Jul 2010 10:39:30 -0400 From: Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu Subject: Re: [gmx-users] pressure coupling To: Discussion list for GROMACS users gmx-users@gromacs.org mailto:gmx-users@gromacs.org Message-ID: 4c31eea2.7070...@vt.edu mailto:4c31eea2.7070...@vt.edu Content-Type: text/plain; charset=UTF-8; format=flowed Shuangxing Dai wrote: Hi, all, I am trying to do anisotropic coupling using Parrinello-Rahman. I want to reach the equilibrium state at 1000K, 0.1MPa( 1 bar). However, I find that the fluctuation is so large ( p=80.22 bar, T=1007.23 K, fluctuation is 626 bar for pressure and 11.21 K for temperature.) I run energy minimization first, then Berenderson pressure coupling for 50 ps, finally Parrinello-Rahman pressure coupling for 50 ps. Anyone has experience on anisotropic coupling and know why? My system has just 4000 atoms. This is my .mdp file: For a small system, large fluctations in the pressure are to be expected: http://www.gromacs.org/Documentation/Terminology/Pressure -Justin ;define = -DPOSRES ; RUN CONTROL PARAMETERS = integrator = sd ; start time and timestep in ps = tinit= 0 dt = 0.001 nsteps = 5 ; number of steps for center of mass motion removal = nstcomm = 100 ; OUTPUT CONTROL OPTIONS = ; Output frequency for coords (x), velocities (v) and forces (f) = nstxout = 0 nstvout = 0 nstfout = 0 ; Output frequency for energies to log file and energy file = nstlog = 10 nstenergy= 10 ; Output frequency and precision for xtc file = nstxtcout= 10 xtc-precision= 1000 ; NEIGHBORSEARCHING PARAMETERS = ; nblist update frequency = nstlist = 20 ; ns algorithm (simple or grid) = ns_type = grid ;OPTIONS FOR PRESSURE COUPLING Pcoupl = Parrinello-Rahman pcoupltype = anisotropic tau_p= 5 compressibility = 2.1645e-07 2.1645e-07 2.7322e-07 0 0 0 ref_p= 1 1 1 0 0 0 ;OPTIONS FOR TEMPERATURE COUPLING tc_grps = system tau_t= 1 ref_t= 1000 ; OPTIONS FOR BONDS = constraints = hbonds ; Type of constraint algorithm = constraint-algorithm = Lincs ; Do not constrain the start configuration = unconstrained-start = no ; Relative tolerance of shake = shake-tol= 0.0001 ; Highest order in the expansion of the constraint coupling matrix = lincs-order = 12 ; Lincs will write a warning to the stderr if in one step a bond = ; rotates over more degrees than = lincs-warnangle = 30 ; Periodic boundary conditions: xyz, no, xy pbc = xyz periodic_molecules = no ; nblist cut-off rlist= 1 ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = PME rcoulomb = 1 ; Method for doing Van der Waals vdw-type = Cut-off ; cut-off lengths rvdw = 1 ; Spacing for the PME/PPPM FFT grid fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 ; EWALD/PME/PPPM parameters pme_order= 6 ewald_rtol = 1e-4 ewald_geometry = 3d epsilon_surface = 0 optimize_fft = no I used the same procedure to reach different state ( like 0.01 K, 100K, 300K, 500K, 800K, 1000K). Only the 0.01 K got the correct average pressure and small fluctuation. T p(bar) delta p T (K) delta T 0.01 1.03E+002.12E+009.95E-031.21E-04 100 5.48359 71.5565 91.198 1.25218 300 2.62807 162.926 290.324 3.36594 500 16.6356 1048.76 492.616 5.45231 800 59.3245 594.773 799.314 9.18876 1000 80.2257 626.789 1007.23 11.2159 Thanks in advance, Shuangxing Dai Thanks
Re: [gmx-users] Pressure coupling and cut-off
Dear Justin: Thank you for your suggestions. using shift for vdwtype helps a lot. Can I ask you another question? The density of my system is about 8% larger than the experimental value. Do you have any suggestions on how to reduce the density of the system. Or once the force field and all the parameters are fixed, the density is independent of the MD procedure? I am using NPT ensemble with the mdp I mentioned before with added vdwtype = shift. On Thu, May 14, 2009 at 4:35 PM, Justin A. Lemkul jalem...@vt.edu wrote: Yanmei Song wrote: Dear Justin: Thanks for your response. Here is the complete my .mdp file: title = pdm cpp = /lib/cpp constraints = all_bonds integrator = md dt = 0.004 ; ps ! nsteps = 250 ; total 10ns. nstcomm = 1 nstxout = 5 nstvout = 5 nstfout = 0 nstlog = 5000 nstenergy = 5000 nstxtcout = 25000 nstlist = 10 ns_type = grid pbc = xyz coulombtype = PME rlist = 1.4 rcoulomb= 1.4 rvdw= 1.4 fourierspacing = 0.20 pme_order = 4 ewald_rtol = 1e-5 ; Berendsen temperature coupling is on in one groups Tcoupl = berendsen tc_grps = PDM tau_t = 0.1 ref_t = 300 ; Energy monitoring energygrps = PDM ; Isotropic pressure coupling is now on Pcoupl = berendsen pcoupltype = isotropic ;pc-grps= PDM tau_p = 1.0 ref_p = 1.0 compressibility = 4.5e-5 ; Generate velocites is off at 300 K. gen_vel = yes gen_temp= 300.0 gen_seed= 10 The problem you're seeing could be an artifact of the shorter cutoff. Have you tried using DispCorr = EnerPres? Or what about using a Shift function for vdwtype? You might see better energy conservation in that case compared to a plain cutoff. -Justin On Thu, May 14, 2009 at 4:06 PM, Justin A. Lemkul jalem...@vt.edumailto: jalem...@vt.edu wrote: Yanmei Song wrote: Dear All: I have question about the pressure coupling. I have done a 10ns simulation with 19800 atoms for 120 large molecules using the following pressure coupling. Tcoupl = berendsen tc_grps = PDM tau_t = 0.1 ref_t = 300 Pcoupl = berendsen pcoupltype = isotropic ;pc-grps= PDM tau_p = 1.0 ref_p = 1.0 compressibility = 4.5e-5 Then I did g_energy for the last 3ns and got the results: Energy Average RMSD Fluct. Drift Tot-Drift --- Potential-98061 0 0 0.616681850.04 Temperature 303.561109.602109.602 0.000181791 0.545372 Pressure (bar) 4.48840.811109.8 -0.169835 -509.506 For such a long run the pressure drift is still too much and seem hasn't approached 1bar. Does it mean the system hasn't reach equilibrium yet. I did a similar system by using the same method. it just take 2 or 3ns to reach the equilibrium. and the pressure is around 1.01after the run. The only difference is the cutoff changing from 1.2 to 1.4. Does the cufoff of 1.4 is too large to make the system running slower. Or the pressure coupling method is not working well. Anyone can give me any suggestions? I think it will depend on the interplay of other parameters as well. Posting a complete .mdp file may be more helpful. -Justin --Yanmei Song Ph.D. Candidate Department of Chemical Engineering Arizona State University ___ gmx-users mailing listgmx-users@gromacs.org mailto:gmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org mailto:gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar Department
Re: [gmx-users] Pressure coupling and cut-off
Yanmei Song wrote: Dear Justin: Thank you for your suggestions. using shift for vdwtype helps a lot. Can I ask you another question? The density of my system is about 8% larger than the experimental value. Do you have any suggestions on how to reduce the density of the system. Or once the force field and all the parameters are fixed, the density is independent of the MD procedure? I am using NPT ensemble with the mdp I mentioned before with added vdwtype = shift. Output like density, pressure, energy conservation, etc. are a consequence of both the force field parameters and their proper use. Look into the primary literature for the force field you are using and see how well the parameters are expected to behave. A difference of 8% doesn't sound that unreasonable to me, but the primary literature will be the judge. -Justin On Thu, May 14, 2009 at 4:35 PM, Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu wrote: Yanmei Song wrote: Dear Justin: Thanks for your response. Here is the complete my .mdp file: title = pdm cpp = /lib/cpp constraints = all_bonds integrator = md dt = 0.004 ; ps ! nsteps = 250 ; total 10ns. nstcomm = 1 nstxout = 5 nstvout = 5 nstfout = 0 nstlog = 5000 nstenergy = 5000 nstxtcout = 25000 nstlist = 10 ns_type = grid pbc = xyz coulombtype = PME rlist = 1.4 rcoulomb= 1.4 rvdw= 1.4 fourierspacing = 0.20 pme_order = 4 ewald_rtol = 1e-5 ; Berendsen temperature coupling is on in one groups Tcoupl = berendsen tc_grps = PDM tau_t = 0.1 ref_t = 300 ; Energy monitoring energygrps = PDM ; Isotropic pressure coupling is now on Pcoupl = berendsen pcoupltype = isotropic ;pc-grps= PDM tau_p = 1.0 ref_p = 1.0 compressibility = 4.5e-5 ; Generate velocites is off at 300 K. gen_vel = yes gen_temp= 300.0 gen_seed= 10 The problem you're seeing could be an artifact of the shorter cutoff. Have you tried using DispCorr = EnerPres? Or what about using a Shift function for vdwtype? You might see better energy conservation in that case compared to a plain cutoff. -Justin On Thu, May 14, 2009 at 4:06 PM, Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu mailto:jalem...@vt.edu mailto:jalem...@vt.edu wrote: Yanmei Song wrote: Dear All: I have question about the pressure coupling. I have done a 10ns simulation with 19800 atoms for 120 large molecules using the following pressure coupling. Tcoupl = berendsen tc_grps = PDM tau_t = 0.1 ref_t = 300 Pcoupl = berendsen pcoupltype = isotropic ;pc-grps= PDM tau_p = 1.0 ref_p = 1.0 compressibility = 4.5e-5 Then I did g_energy for the last 3ns and got the results: Energy Average RMSD Fluct. Drift Tot-Drift --- Potential-98061 0 0 0.616681850.04 Temperature 303.561109.602109.602 0.000181791 0.545372 Pressure (bar) 4.48840.811109.8 -0.169835 -509.506 For such a long run the pressure drift is still too much and seem hasn't approached 1bar. Does it mean the system hasn't reach equilibrium yet. I did a similar system by using the same method. it just take 2 or 3ns to reach the equilibrium. and the pressure is around 1.01after the run. The only difference is the cutoff changing from 1.2 to 1.4. Does the cufoff of 1.4 is too large to make the system running slower. Or the pressure
Re: [gmx-users] Pressure coupling and cut-off
Yanmei Song wrote: Dear All: I have question about the pressure coupling. I have done a 10ns simulation with 19800 atoms for 120 large molecules using the following pressure coupling. Tcoupl = berendsen tc_grps = PDM tau_t = 0.1 ref_t = 300 Pcoupl = berendsen pcoupltype = isotropic ;pc-grps= PDM tau_p = 1.0 ref_p = 1.0 compressibility = 4.5e-5 Then I did g_energy for the last 3ns and got the results: Energy Average RMSD Fluct. Drift Tot-Drift --- Potential-98061 0 00.61668 1850.04 Temperature 303.561109.602109.602 0.000181791 0.545372 Pressure (bar) 4.48840.811109.8 -0.169835 -509.506 For such a long run the pressure drift is still too much and seem hasn't approached 1bar. Does it mean the system hasn't reach equilibrium yet. I did a similar system by using the same method. it just take 2 or 3ns to reach the equilibrium. and the pressure is around 1.01after the run. The only difference is the cutoff changing from 1.2 to 1.4. Does the cufoff of 1.4 is too large to make the system running slower. Or the pressure coupling method is not working well. Anyone can give me any suggestions? I think it will depend on the interplay of other parameters as well. Posting a complete .mdp file may be more helpful. -Justin -- Yanmei Song Ph.D. Candidate Department of Chemical Engineering Arizona State University ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Pressure coupling and cut-off
Dear Justin: Thanks for your response. Here is the complete my .mdp file: title = pdm cpp = /lib/cpp constraints = all_bonds integrator = md dt = 0.004 ; ps ! nsteps = 250 ; total 10ns. nstcomm = 1 nstxout = 5 nstvout = 5 nstfout = 0 nstlog = 5000 nstenergy = 5000 nstxtcout = 25000 nstlist = 10 ns_type = grid pbc = xyz coulombtype = PME rlist = 1.4 rcoulomb= 1.4 rvdw= 1.4 fourierspacing = 0.20 pme_order = 4 ewald_rtol = 1e-5 ; Berendsen temperature coupling is on in one groups Tcoupl = berendsen tc_grps = PDM tau_t = 0.1 ref_t = 300 ; Energy monitoring energygrps = PDM ; Isotropic pressure coupling is now on Pcoupl = berendsen pcoupltype = isotropic ;pc-grps= PDM tau_p = 1.0 ref_p = 1.0 compressibility = 4.5e-5 ; Generate velocites is off at 300 K. gen_vel = yes gen_temp= 300.0 gen_seed= 10 On Thu, May 14, 2009 at 4:06 PM, Justin A. Lemkul jalem...@vt.edu wrote: Yanmei Song wrote: Dear All: I have question about the pressure coupling. I have done a 10ns simulation with 19800 atoms for 120 large molecules using the following pressure coupling. Tcoupl = berendsen tc_grps = PDM tau_t = 0.1 ref_t = 300 Pcoupl = berendsen pcoupltype = isotropic ;pc-grps= PDM tau_p = 1.0 ref_p = 1.0 compressibility = 4.5e-5 Then I did g_energy for the last 3ns and got the results: Energy Average RMSD Fluct. Drift Tot-Drift --- Potential-98061 0 00.61668 1850.04 Temperature 303.561109.602109.602 0.000181791 0.545372 Pressure (bar) 4.48840.811109.8 -0.169835 -509.506 For such a long run the pressure drift is still too much and seem hasn't approached 1bar. Does it mean the system hasn't reach equilibrium yet. I did a similar system by using the same method. it just take 2 or 3ns to reach the equilibrium. and the pressure is around 1.01after the run. The only difference is the cutoff changing from 1.2 to 1.4. Does the cufoff of 1.4 is too large to make the system running slower. Or the pressure coupling method is not working well. Anyone can give me any suggestions? I think it will depend on the interplay of other parameters as well. Posting a complete .mdp file may be more helpful. -Justin -- Yanmei Song Ph.D. Candidate Department of Chemical Engineering Arizona State University ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php -- Yanmei Song Ph.D. Candidate Department of Chemical Engineering Arizona State University ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Pressure coupling and cut-off
Yanmei Song wrote: Dear Justin: Thanks for your response. Here is the complete my .mdp file: title = pdm cpp = /lib/cpp constraints = all_bonds integrator = md dt = 0.004 ; ps ! nsteps = 250 ; total 10ns. nstcomm = 1 nstxout = 5 nstvout = 5 nstfout = 0 nstlog = 5000 nstenergy = 5000 nstxtcout = 25000 nstlist = 10 ns_type = grid pbc = xyz coulombtype = PME rlist = 1.4 rcoulomb= 1.4 rvdw= 1.4 fourierspacing = 0.20 pme_order = 4 ewald_rtol = 1e-5 ; Berendsen temperature coupling is on in one groups Tcoupl = berendsen tc_grps = PDM tau_t = 0.1 ref_t = 300 ; Energy monitoring energygrps = PDM ; Isotropic pressure coupling is now on Pcoupl = berendsen pcoupltype = isotropic ;pc-grps= PDM tau_p = 1.0 ref_p = 1.0 compressibility = 4.5e-5 ; Generate velocites is off at 300 K. gen_vel = yes gen_temp= 300.0 gen_seed= 10 The problem you're seeing could be an artifact of the shorter cutoff. Have you tried using DispCorr = EnerPres? Or what about using a Shift function for vdwtype? You might see better energy conservation in that case compared to a plain cutoff. -Justin On Thu, May 14, 2009 at 4:06 PM, Justin A. Lemkul jalem...@vt.edu mailto:jalem...@vt.edu wrote: Yanmei Song wrote: Dear All: I have question about the pressure coupling. I have done a 10ns simulation with 19800 atoms for 120 large molecules using the following pressure coupling. Tcoupl = berendsen tc_grps = PDM tau_t = 0.1 ref_t = 300 Pcoupl = berendsen pcoupltype = isotropic ;pc-grps= PDM tau_p = 1.0 ref_p = 1.0 compressibility = 4.5e-5 Then I did g_energy for the last 3ns and got the results: Energy Average RMSD Fluct. Drift Tot-Drift --- Potential-98061 0 0 0.616681850.04 Temperature 303.561109.602109.602 0.000181791 0.545372 Pressure (bar) 4.48840.811109.8 -0.169835 -509.506 For such a long run the pressure drift is still too much and seem hasn't approached 1bar. Does it mean the system hasn't reach equilibrium yet. I did a similar system by using the same method. it just take 2 or 3ns to reach the equilibrium. and the pressure is around 1.01after the run. The only difference is the cutoff changing from 1.2 to 1.4. Does the cufoff of 1.4 is too large to make the system running slower. Or the pressure coupling method is not working well. Anyone can give me any suggestions? I think it will depend on the interplay of other parameters as well. Posting a complete .mdp file may be more helpful. -Justin -- Yanmei Song Ph.D. Candidate Department of Chemical Engineering Arizona State University ___ gmx-users mailing listgmx-users@gromacs.org mailto:gmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org mailto:gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php -- Justin A. Lemkul Ph.D. Candidate ICTAS Doctoral Scholar Department of Biochemistry Virginia Tech Blacksburg, VA jalemkul[at]vt.edu http://vt.edu | (540) 231-9080 http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin ___ gmx-users mailing listgmx-users@gromacs.org mailto:gmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests
Re: [gmx-users] Pressure Coupling Problem
I don't know the conditions you use with your system, but sometimes there are problems when you use pressure coupling in a system with positionally restrained or fixed atoms. Cheers. Lucio. From: Joe Joe Sent: Friday, April 10, 2009 12:12 PM To: Discussion list for GROMACS users Subject: Re: [gmx-users] Pressure Coupling Problem I finally figured it out. I went through every parameter step by step and it turns our I had epsilon-r set to 80. Not sure why I had that. Wish gromacs would have given me a warning (hint hint). That explains why my P.E. was 10^-5 instead of 10^-6. Thanks everyone for trying!!!. Sometime the most obvious mistake takes some time to figure out. Cheers, Ilya On Fri, Apr 10, 2009 at 8:05 AM, chris.ne...@utoronto.ca wrote: Alright, sorry that I wasn't able to help. I'm confused by some apparent contradictions in your posts and I'm not sure that I'm going to be useful to you here. Quoting http://www.gromacs.org/pipermail/gmx-users/2009-April/041173.html: No matter how much minimization I do the volume of the box expands when run it using berendsen pressure coupling with a tau_p -f .1. Quoting http://www.gromacs.org/pipermail/gmx-users/2009-April/041159.html: So I got my small water box (800 waters) to behave stably with pressure coupling after more minimization ... Good luck. Chris. -- original message -- Nope not an A/nm problem. As a simple test I take spc.gro from share/top. I reconfigure the box (i.e. editconf -f spc.gro -d 1.0 -c -bt cubic -o water_center. I then solvate with genbox, minimize and run using the mdp file I provided earlier. No matter how much minimization I do the volume of the box expands when run it using berendsen pressure coupling with a tau_p -f .1. Ilya On Thu, Apr 9, 2009 at 11:22 AM, Chris Neale chris.ne...@utoronto.cawrote: [Hide Quoted Text] So your problem with the small water box was solved simply by adding more minimization? I then suspect that all of your problems are simply related to a bad starting structure -- and by the sound of it is really is very bad. Are you sure that you don't have an angstrom / nm problem here? Chris. -- original message -- So I got my small water box (800 waters) to behave stably with pressure coupling after more minimization but I still can't get my large system to work with pressure coupling. I tried minimizing but I can never get the Fmax to be less 102, which is pretty normal for protein/water simulations of large proteins, at least from my experience. I have since run 400 ps NVT as ... ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use thewww interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Pressure Coupling Problem
So I got my small water box (800 waters) to behave stably with pressure coupling after more minimization but I still can't get my large system to work with pressure coupling. I tried minimizing but I can never get the Fmax to be less 10^2, which is pretty normal for protein/water simulations of large proteins, at least from my experience. I have since run 400 ps NVT as the system (425K atoms) is quite stable. The P.E. is 2E-05. Since I am using 4fs time steps gromacs won't let me use a tau_p less than .4. Not sure what else to do except run NVT, which is what I was going to do after I got the density equilibrated. BTW, I am using octahedral PBC, but that should not make a difference with respect to P coupling, should it? Below is my whole mdp file. As a reminder my density in the system goes from 1.0 - .1 in 10 ps with Pcoupl = Berendsen and Tau_p = .4. If I increase Tau_P then the amount of time it takes for my system to expand increases but it still expands. ; ; File 'mdout.mdp' was generated ; By user: relly (508) ; On host: master.simprota.com ; At date: Fri Mar 6 20:17:33 2009 ; ; VARIOUS PREPROCESSING OPTIONS ; Preprocessor information: use cpp syntax. ; e.g.: -I/home/joe/doe -I/home/mary/hoe include = ; e.g.: -DI_Want_Cookies -DMe_Too define = ; RUN CONTROL PARAMETERS integrator = md ; Start time and timestep in ps tinit= 0 dt = 0.004 ;nsteps = 25 nsteps = 250 ; For exact run continuation or redoing part of a run ; Part index is updated automatically on checkpointing (keeps files separate) simulation_part = 1 init_step= 0 ; mode for center of mass motion removal comm_mode= linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm_grps= system ; LANGEVIN DYNAMICS OPTIONS ; Friction coefficient (amu/ps) and random seed bd-fric = 0 ld-seed = 1993 ; ENERGY MINIMIZATION OPTIONS ; Force tolerance and initial step-size emtol= 10 emstep = 0.01 ; Max number of iterations in relax_shells niter= 20 ; Step size (ps^2) for minimization of flexible constraints fcstep = 0 ; Frequency of steepest descents steps when doing CG nstcgsteep = 1000 nbfgscorr= 10 ; TEST PARTICLE INSERTION OPTIONS rtpi = 0.05 ; OUTPUT CONTROL OPTIONS ; Output frequency for coords (x), velocities (v) and forces (f) nstxout = 12500 nstvout = 0 nstfout = 0 ; Output frequency for energies to log file and energy file nstlog = 10 nstenergy= 10 ; Output frequency and precision for xtc file nstxtcout= 250 xtc-precision= 1000 ; This selects the subset of atoms for the xtc file. You can ; select multiple groups. By default all atoms will be written. xtc-grps = protein ; Selection of energy groups energygrps = Protein SOL ; NEIGHBORSEARCHING PARAMETERS ; nblist update frequency nstlist = 5 ; ns algorithm (simple or grid) ns_type = grid ; Periodic boundary conditions: xyz, no, xy pbc = xyz periodic_molecules = no ; nblist cut-off rlist= 1.0 ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = PME rcoulomb-switch = .9 rcoulomb = 1.0 ; Relative dielectric constant for the medium and the reaction field epsilon-r= 80 epsilon_rf = 1 ; Method for doing Van der Waals vdw-type = Switch ; cut-off lengths rvdw-switch = .8 rvdw = 1.0 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; Extension of the potential lookup tables beyond the cut-off table-extension = 1 ; Seperate tables between energy group pairs energygrp_table = ; Spacing for the PME/PPPM FFT grid fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 ; EWALD/PME/PPPM parameters pme_order= 4 ewald_rtol = 1.e-05 ewald_geometry = 3d epsilon_surface = 0 optimize_fft = no ; IMPLICIT SOLVENT ALGORITHM implicit_solvent = No ; GENERALIZED BORN ELECTROSTATICS ; Algorithm for calculating Born radii gb_algorithm = Still ; Frequency of calculating the Born radii inside rlist nstgbradii = 1 ; Cutoff for Born radii calculation; the contribution from atoms ; between rlist and rgbradii is updated every
Re: [gmx-users] Pressure Coupling Problem
On Thu, Apr 9, 2009 at 6:36 AM, Justin A. Lemkul jalem...@vt.edu wrote: Joe Joe wrote: So I got my small water box (800 waters) to behave stably with pressure coupling after more minimization but I still can't get my large system to work with pressure coupling. I tried minimizing but I can never get the Fmax to be less 10^2, which is pretty normal for protein/water simulations of large proteins, at least from my experience. I have since run 400 ps NVT as the system (425K atoms) is quite stable. The P.E. is 2E-05. Since I am using 4fs time steps gromacs won't let me use a tau_p less than .4. Not sure what else to do except run NVT, which is what I was going to do after I got the density equilibrated. BTW, I am using octahedral PBC, but that should not make a difference with respect to P coupling, should it? Below is my whole mdp file. As a reminder my density in the system goes from 1.0 - .1 in 10 ps with Pcoupl = Berendsen and Tau_p = .4. If I increase Tau_P then the amount of time it takes for my system to expand increases but it still expands. This seems truly bizarre. How are you measuring the density (g_density, g_energy, etc)? Both g_energy and g_density. What are your box dimensions doing? To get that kind of sudden change in density, your box dimensions would have to expand astronomically? Yep. It's also curious that your 425K-atom system only has a PE on the order of 10^5; my systems with 100K-200K have around 10^6 - 10^7; are you sure the minimization is reasonable, and you are not simply seeing the effects of the classic blowing up problem? If that was the case would not the NVT also not behave stably? I also agree that 10^5 seems to high. Most of that should come from water though, correct? Why would the water not relax. Maybe I should just expand the box a bit and see what happens? What does your trajectory show? If you have multiple proteins or other large species present, does minimization of each component individually prior to system assembly help? I have one very large protein (antibody). I did do a gas phase minimization of the protein prior to solvation and it does not help. ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Pressure Coupling Problem
HI Chris, On Tue, Apr 7, 2009 at 9:31 PM, chris.ne...@utoronto.ca wrote: Hi Ilya, First thing that comes to mind is that it is strange to couple a coulombic switching function with PME. While this could possibly be done correctly, I doubt that it is in fact done in the way that you expect (i.e. correctly) in gromacs. In fact, I think that grompp/mdrun should probably throw an error here -- unless it is actually handled in the proper way, and a developer could help you here to figure out if you are indeed getting what you desire. coulombtype = PME rcoulomb-switch = .9 rcoulomb = 1.0 I am pretty sure gromacs ignores the rcoulomb-switch parameter in the case of PME but I will give it a try. However, it is not clear to me that this should cause a system to continuously expand. Still, you do not give very good information about what you mean by continuously expand. Can you please provide some information on that? e.g. amount of time and total volume change. My box density goes from ~1.0 to .5 in 5 ps with a compressibility of 5E-05. It goes from ~1.0 to .94 in 300 ps with a compressibility of 5E-06. In both case the slope of density(t) is negative and never levels off. Chris -- original message -- Hi I am having some pressure coupling issues. I have a fairly large protein/water system 400K+ atoms. It minimizes just fine (F 1000). If I run NVE it conserves energy with appropriate parameter settings. If I run NVT it is stable. When I turn on Pcoupl (i.e. Berendsen or Parinello Rahman), the system just continuously expands. My parameters are as follows. Any ideas? Best, Ilya ; ; File 'mdout.mdp' was generated ; By user: relly (508) ; On host: master.simprota.com ; At date: Fri Mar 6 20:17:33 2009 ; ; VARIOUS PREPROCESSING OPTIONS ; Preprocessor information: use cpp syntax. ; e.g.: -I/home/joe/doe -I/home/mary/hoe include = ; e.g.: -DI_Want_Cookies -DMe_Too define = ; RUN CONTROL PARAMETERS integrator = md ; Start time and timestep in ps tinit= 0 dt = 0.004 ;nsteps = 25 nsteps = 250 ; For exact run continuation or redoing part of a run ; Part index is updated automatically on checkpointing (keeps files separate) simulation_part = 1 init_step= 0 ; mode for center of mass motion removal comm_mode= linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm_grps= system ; OUTPUT CONTROL OPTIONS ; Output frequency for coords (x), velocities (v) and forces (f) nstxout = 0 nstvout = 0 nstfout = 0 ; Output frequency for energies to log file and energy file nstlog = 10 nstenergy= 10 ; Output frequency and precision for xtc file nstxtcout= 250 xtc-precision= 1000 ; This selects the subset of atoms for the xtc file. You can ; select multiple groups. By default all atoms will be written. xtc-grps = protein ; Selection of energy groups energygrps = ; NEIGHBORSEARCHING PARAMETERS ; nblist update frequency nstlist = 5 ; ns algorithm (simple or grid) ns_type = grid ; Periodic boundary conditions: xyz, no, xy pbc = xyz periodic_molecules = no ; nblist cut-off rlist= 1.0 ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = PME rcoulomb-switch = .9 rcoulomb = 1.0 ; Relative dielectric constant for the medium and the reaction field epsilon-r= 80 epsilon_rf = 1 ; Method for doing Van der Waals vdw-type = Switch ; cut-off lengths rvdw-switch = .9 rvdw = 1.0 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; Extension of the potential lookup tables beyond the cut-off table-extension = 1 ; Seperate tables between energy group pairs energygrp_table = ; Spacing for the PME/PPPM FFT grid fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 ; EWALD/PME/PPPM parameters pme_order= 4 ewald_rtol = 1.e-05 ewald_geometry = 3d epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS ; Temperature coupling Tcoupl = V-rescale ; Groups to couple separately tc-grps = System ; Time constant (ps) and
Re: [gmx-users] Pressure Coupling Problem
Hi Chris, When I create the topology for the 4fs timestep I use pdb2gmx -vsite h. I set up the correct constraints. I've tested it and it conserves energy in NVE. I run all he sims with constraints=all-bonds. I am now running a single water box (800 water molecules) with 1s time steps and the volume keeps blowing up. Thanks, Ilya On Wed, Apr 8, 2009 at 8:37 AM, chris.ne...@utoronto.ca wrote: Hi Ilya, If you did include the entire mdp file then you have a time step of 4 fs and no constraints (other than water). For a timestep of 2 fs, you should constrain all-bonds (or some would say at least h-bonds) and for 4 fs then you should also constrain angles involving hydrogens (need a new .itp file for this). Can you try with a 1 fs timestep and see how it goes? Still, I am surprised that everything works out at NVT, but this is certainly worth the test. Do you have other systems running fine with these mdp options in NVT? Chris. -- original message -- HI Chris, On Tue, Apr 7, 2009 at 9:31 PM, chris.neale at utoronto.ca wrote: Hi Ilya, First thing that comes to mind is that it is strange to couple a coulombic switching function with PME. While this could possibly be done correctly, I doubt that it is in fact done in the way that you expect (i.e. correctly) in gromacs. In fact, I think that grompp/mdrun should probably throw an error here -- unless it is actually handled in the proper way, and a developer could help you here to figure out if you are indeed getting what you desire. coulombtype = PME rcoulomb-switch = .9 rcoulomb = 1.0 I am pretty sure gromacs ignores the rcoulomb-switch parameter in the case of PME but I will give it a try. However, it is not clear to me that this should cause a system to continuously expand. Still, you do not give very good information about what you mean by continuously expand. Can you please provide some information on that? e.g. amount of time and total volume change. My box density goes from ~1.0 to .5 in 5 ps with a compressibility of 5E-05. It goes from ~1.0 to .94 in 300 ps with a compressibility of 5E-06. In both case the slope of density(t) is negative and never levels off. Chris -- original message -- Hi I am having some pressure coupling issues. I have a fairly large protein/water system 400K+ atoms. It minimizes just fine (F 1000). If I run NVE it conserves energy with appropriate parameter settings. If I run NVT it is stable. When I turn on Pcoupl (i.e. Berendsen or Parinello Rahman), the system just continuously expands. My parameters are as follows. Any ideas? Best, Ilya ; ; File 'mdout.mdp' was generated ; By user: relly (508) ; On host: master.simprota.com ; At date: Fri Mar 6 20:17:33 2009 ; ; VARIOUS PREPROCESSING OPTIONS ; Preprocessor information: use cpp syntax. ; e.g.: -I/home/joe/doe -I/home/mary/hoe include = ; e.g.: -DI_Want_Cookies -DMe_Too define = ; RUN CONTROL PARAMETERS integrator = md ; Start time and timestep in ps tinit= 0 dt = 0.004 ;nsteps = 25 nsteps = 250 ; For exact run continuation or redoing part of a run ; Part index is updated automatically on checkpointing (keeps files separate) simulation_part = 1 init_step= 0 ; mode for center of mass motion removal comm_mode= linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm_grps= system ; OUTPUT CONTROL OPTIONS ; Output frequency for coords (x), velocities (v) and forces (f) nstxout = 0 nstvout = 0 nstfout = 0 ; Output frequency for energies to log file and energy file nstlog = 10 nstenergy= 10 ; Output frequency and precision for xtc file nstxtcout= 250 xtc-precision= 1000 ; This selects the subset of atoms for the xtc file. You can ; select multiple groups. By default all atoms will be written. xtc-grps = protein ; Selection of energy groups energygrps = ; NEIGHBORSEARCHING PARAMETERS ; nblist update frequency nstlist = 5 ; ns algorithm (simple or grid) ns_type = grid ; Periodic boundary conditions: xyz, no, xy pbc = xyz periodic_molecules = no ; nblist cut-off rlist= 1.0 ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = PME rcoulomb-switch = .9 rcoulomb = 1.0 ; Relative dielectric constant for the medium and the reaction field epsilon-r= 80 epsilon_rf = 1 ; Method for doing Van der
Re: [gmx-users] Pressure Coupling Problem
Joe Joe wrote: Hi Chris, When I create the topology for the 4fs timestep I use pdb2gmx -vsite h. I set up the correct constraints. I've tested it and it conserves energy in NVE. I run all he sims with constraints=all-bonds. I am now running a single water box (800 water molecules) with 1s time steps and the volume keeps blowing up. In addition to what Chris has been saying about constraints, consider your pressure coupling settings themselves Pcoupl = Berendsen Pcoupltype = Isotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) tau_p= 10 compressibility = 4.5e-5 ref_p= 1.01325 A 10-ps relaxation time for a system that is not necessarily well-equilibrated is too weak, I think. Try 1.0 - 2.0 ps for tau_p. If your system is expanding rapidly in as little as 5 ps, I would think you lack appropriate pressure regulation. -Justin Thanks, Ilya On Wed, Apr 8, 2009 at 8:37 AM, chris.ne...@utoronto.ca mailto:chris.ne...@utoronto.ca wrote: Hi Ilya, If you did include the entire mdp file then you have a time step of 4 fs and no constraints (other than water). For a timestep of 2 fs, you should constrain all-bonds (or some would say at least h-bonds) and for 4 fs then you should also constrain angles involving hydrogens (need a new .itp file for this). Can you try with a 1 fs timestep and see how it goes? Still, I am surprised that everything works out at NVT, but this is certainly worth the test. Do you have other systems running fine with these mdp options in NVT? Chris. -- original message -- HI Chris, On Tue, Apr 7, 2009 at 9:31 PM, chris.neale at utoronto.ca http://utoronto.ca wrote: Hi Ilya, First thing that comes to mind is that it is strange to couple a coulombic switching function with PME. While this could possibly be done correctly, I doubt that it is in fact done in the way that you expect (i.e. correctly) in gromacs. In fact, I think that grompp/mdrun should probably throw an error here -- unless it is actually handled in the proper way, and a developer could help you here to figure out if you are indeed getting what you desire. coulombtype = PME rcoulomb-switch = .9 rcoulomb = 1.0 I am pretty sure gromacs ignores the rcoulomb-switch parameter in the case of PME but I will give it a try. However, it is not clear to me that this should cause a system to continuously expand. Still, you do not give very good information about what you mean by continuously expand. Can you please provide some information on that? e.g. amount of time and total volume change. My box density goes from ~1.0 to .5 in 5 ps with a compressibility of 5E-05. It goes from ~1.0 to .94 in 300 ps with a compressibility of 5E-06. In both case the slope of density(t) is negative and never levels off. Chris -- original message -- Hi I am having some pressure coupling issues. I have a fairly large protein/water system 400K+ atoms. It minimizes just fine (F 1000). If I run NVE it conserves energy with appropriate parameter settings. If I run NVT it is stable. When I turn on Pcoupl (i.e. Berendsen or Parinello Rahman), the system just continuously expands. My parameters are as follows. Any ideas? Best, Ilya ; ; File 'mdout.mdp' was generated ; By user: relly (508) ; On host: master.simprota.com http://master.simprota.com ; At date: Fri Mar 6 20:17:33 2009 ; ; VARIOUS PREPROCESSING OPTIONS ; Preprocessor information: use cpp syntax. ; e.g.: -I/home/joe/doe -I/home/mary/hoe include = ; e.g.: -DI_Want_Cookies -DMe_Too define = ; RUN CONTROL PARAMETERS integrator = md ; Start time and timestep in ps tinit= 0 dt = 0.004 ;nsteps = 25 nsteps = 250 ; For exact run continuation or redoing part of a run ; Part index is updated automatically on checkpointing (keeps files separate) simulation_part = 1 init_step= 0 ; mode for center of mass motion removal comm_mode= linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm_grps= system ; OUTPUT
Re: [gmx-users] Pressure Coupling Problem
On Wed, Apr 8, 2009 at 7:53 AM, Joe Joe ilcho...@gmail.com wrote: HI Chris, On Tue, Apr 7, 2009 at 9:31 PM, chris.ne...@utoronto.ca wrote: Hi Ilya, First thing that comes to mind is that it is strange to couple a coulombic switching function with PME. While this could possibly be done correctly, I doubt that it is in fact done in the way that you expect (i.e. correctly) in gromacs. In fact, I think that grompp/mdrun should probably throw an error here -- unless it is actually handled in the proper way, and a developer could help you here to figure out if you are indeed getting what you desire. coulombtype = PME rcoulomb-switch = .9 rcoulomb = 1.0 I am pretty sure gromacs ignores the rcoulomb-switch parameter in the case of PME but I will give it a try. It is indeed supported and does work correctly. But you have to set coulombtype PME-Switch. mdp options says: This is mainly useful constant energy simulations. For constant temperature simulations the advantage of improved energy conservation is usually outweighed by the small loss in accuracy of the electrostatics. Roland Chris -- original message -- Hi I am having some pressure coupling issues. I have a fairly large protein/water system 400K+ atoms. It minimizes just fine (F 1000). If I run NVE it conserves energy with appropriate parameter settings. If I run NVT it is stable. When I turn on Pcoupl (i.e. Berendsen or Parinello Rahman), the system just continuously expands. My parameters are as follows. Any ideas? Best, Ilya ; ; File 'mdout.mdp' was generated ; By user: relly (508) ; On host: master.simprota.com ; At date: Fri Mar 6 20:17:33 2009 ; ; VARIOUS PREPROCESSING OPTIONS ; Preprocessor information: use cpp syntax. ; e.g.: -I/home/joe/doe -I/home/mary/hoe include = ; e.g.: -DI_Want_Cookies -DMe_Too define = ; RUN CONTROL PARAMETERS integrator = md ; Start time and timestep in ps tinit= 0 dt = 0.004 ;nsteps = 25 nsteps = 250 ; For exact run continuation or redoing part of a run ; Part index is updated automatically on checkpointing (keeps files separate) simulation_part = 1 init_step= 0 ; mode for center of mass motion removal comm_mode= linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm_grps= system ; OUTPUT CONTROL OPTIONS ; Output frequency for coords (x), velocities (v) and forces (f) nstxout = 0 nstvout = 0 nstfout = 0 ; Output frequency for energies to log file and energy file nstlog = 10 nstenergy= 10 ; Output frequency and precision for xtc file nstxtcout= 250 xtc-precision= 1000 ; This selects the subset of atoms for the xtc file. You can ; select multiple groups. By default all atoms will be written. xtc-grps = protein ; Selection of energy groups energygrps = ; NEIGHBORSEARCHING PARAMETERS ; nblist update frequency nstlist = 5 ; ns algorithm (simple or grid) ns_type = grid ; Periodic boundary conditions: xyz, no, xy pbc = xyz periodic_molecules = no ; nblist cut-off rlist= 1.0 ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = PME rcoulomb-switch = .9 rcoulomb = 1.0 ; Relative dielectric constant for the medium and the reaction field epsilon-r= 80 epsilon_rf = 1 ; Method for doing Van der Waals vdw-type = Switch ; cut-off lengths rvdw-switch = .9 rvdw = 1.0 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; Extension of the potential lookup tables beyond the cut-off table-extension = 1 ; Seperate tables between energy group pairs energygrp_table = ; Spacing for the PME/PPPM FFT grid fourierspacing = 0.12 ; FFT grid size, when a value is 0 fourierspacing will be used fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 ; EWALD/PME/PPPM parameters pme_order= 4 ewald_rtol = 1.e-05 ewald_geometry = 3d epsilon_surface = 0 optimize_fft = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS ; Temperature coupling Tcoupl = V-rescale ; Groups to couple separately tc-grps = System ; Time constant (ps) and reference temperature (K) tau_t= 0.1 ref_t= 298.0 ;
Re: [gmx-users] Pressure Coupling Problem
Yeah I only gave a partial. Tried to remove the QM params. I do use constraints = all-bonds. On Wed, Apr 8, 2009 at 9:18 AM, chris.ne...@utoronto.ca wrote: You say I run all he sims with constraints=all-bonds, but I don't see that in the mdp options that you provided. I even put your text in a file and grepped for it just to be sure. Did you only give us a partial mdp file? Try adding this to your mdp file: constraints = all-bonds ; REMOVE_FOR_EM constraint_algorithm= lincs ; REMOVE_FOR_EM lincs-iter = 1 ; REMOVE_FOR_EM lincs-order = 6 ; REMOVE_FOR_EM see, for example, http://www.gromacs.org/pipermail/gmx-users/2008-October/037545.html http://www.gromacs.org/pipermail/gmx-users/2008-November/037673.html --- original message --- Hi Chris, When I create the topology for the 4fs timestep I use pdb2gmx -vsite h. I set up the correct constraints. I've tested it and it conserves energy in NVE. I run all he sims with constraints=all-bonds. I am now running a single water box (800 water molecules) with 1s time steps and the volume keeps blowing up. Thanks, Ilya On Wed, Apr 8, 2009 at 8:37 AM, chris.neale at utoronto.ca wrote: Hi Ilya, If you did include the entire mdp file then you have a time step of 4 fs and no constraints (other than water). For a timestep of 2 fs, you should constrain all-bonds (or some would say at least h-bonds) and for 4 fs then you should also constrain angles involving hydrogens (need a new .itp file for this). Can you try with a 1 fs timestep and see how it goes? Still, I am surprised that everything works out at NVT, but this is certainly worth the test. Do you have other systems running fine with these mdp options in NVT? Chris. -- original message -- HI Chris, On Tue, Apr 7, 2009 at 9:31 PM, chris.neale at utoronto.ca wrote: Hi Ilya, First thing that comes to mind is that it is strange to couple a coulombic switching function with PME. While this could possibly be done correctly, I doubt that it is in fact done in the way that you expect (i.e. correctly) in gromacs. In fact, I think that grompp/mdrun should probably throw an error here -- unless it is actually handled in the proper way, and a developer could help you here to figure out if you are indeed getting what you desire. coulombtype = PME rcoulomb-switch = .9 rcoulomb = 1.0 I am pretty sure gromacs ignores the rcoulomb-switch parameter in the case of PME but I will give it a try. However, it is not clear to me that this should cause a system to continuously expand. Still, you do not give very good information about what you mean by continuously expand. Can you please provide some information on that? e.g. amount of time and total volume change. My box density goes from ~1.0 to .5 in 5 ps with a compressibility of 5E-05. It goes from ~1.0 to .94 in 300 ps with a compressibility of 5E-06. In both case the slope of density(t) is negative and never levels off. Chris -- original message -- Hi I am having some pressure coupling issues. I have a fairly large protein/water system 400K+ atoms. It minimizes just fine (F 1000). If I run NVE it conserves energy with appropriate parameter settings. If I run NVT it is stable. When I turn on Pcoupl (i.e. Berendsen or Parinello Rahman), the system just continuously expands. My parameters are as follows. Any ideas? Best, Ilya ; ; File 'mdout.mdp' was generated ; By user: relly (508) ; On host: master.simprota.com ; At date: Fri Mar 6 20:17:33 2009 ; ; VARIOUS PREPROCESSING OPTIONS ; Preprocessor information: use cpp syntax. ; e.g.: -I/home/joe/doe -I/home/mary/hoe include = ; e.g.: -DI_Want_Cookies -DMe_Too define = ; RUN CONTROL PARAMETERS integrator = md ; Start time and timestep in ps tinit= 0 dt = 0.004 ;nsteps = 25 nsteps = 250 ; For exact run continuation or redoing part of a run ; Part index is updated automatically on checkpointing (keeps files separate) simulation_part = 1 init_step= 0 ; mode for center of mass motion removal comm_mode= linear ; number of steps for center of mass motion removal nstcomm = 1 ; group(s) for center of mass motion removal comm_grps= system ; OUTPUT CONTROL OPTIONS ; Output frequency for coords (x), velocities (v) and forces (f) nstxout = 0 nstvout = 0 nstfout = 0 ; Output frequency for energies to log file and energy file nstlog = 10 nstenergy= 10 ; Output frequency and precision for xtc file nstxtcout= 250 xtc-precision= 1000 ;
Re: [gmx-users] Pressure Coupling Problem
I tried .1, and 10 ps tau_p values. I guess I can try smaller values. On Wed, Apr 8, 2009 at 10:44 AM, Justin A. Lemkul jalem...@vt.edu wrote: Joe Joe wrote: Hi Chris, When I create the topology for the 4fs timestep I use pdb2gmx -vsite h. I set up the correct constraints. I've tested it and it conserves energy in NVE. I run all he sims with constraints=all-bonds. I am now running a single water box (800 water molecules) with 1s time steps and the volume keeps blowing up. In addition to what Chris has been saying about constraints, consider your pressure coupling settings themselves Pcoupl = Berendsen Pcoupltype = Isotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) tau_p= 10 compressibility = 4.5e-5 ref_p= 1.01325 A 10-ps relaxation time for a system that is not necessarily well-equilibrated is too weak, I think. Try 1.0 - 2.0 ps for tau_p. If your system is expanding rapidly in as little as 5 ps, I would think you lack appropriate pressure regulation. -Justin Thanks, Ilya On Wed, Apr 8, 2009 at 8:37 AM, chris.ne...@utoronto.ca mailto: chris.ne...@utoronto.ca wrote: Hi Ilya, If you did include the entire mdp file then you have a time step of 4 fs and no constraints (other than water). For a timestep of 2 fs, you should constrain all-bonds (or some would say at least h-bonds) and for 4 fs then you should also constrain angles involving hydrogens (need a new .itp file for this). Can you try with a 1 fs timestep and see how it goes? Still, I am surprised that everything works out at NVT, but this is certainly worth the test. Do you have other systems running fine with these mdp options in NVT? Chris. -- original message -- HI Chris, On Tue, Apr 7, 2009 at 9:31 PM, chris.neale at utoronto.ca http://utoronto.ca wrote: Hi Ilya, First thing that comes to mind is that it is strange to couple a coulombic switching function with PME. While this could possibly be done correctly, I doubt that it is in fact done in the way that you expect (i.e. correctly) in gromacs. In fact, I think that grompp/mdrun should probably throw an error here -- unless it is actually handled in the proper way, and a developer could help you here to figure out if you are indeed getting what you desire. coulombtype = PME rcoulomb-switch = .9 rcoulomb = 1.0 I am pretty sure gromacs ignores the rcoulomb-switch parameter in the case of PME but I will give it a try. However, it is not clear to me that this should cause a system to continuously expand. Still, you do not give very good information about what you mean by continuously expand. Can you please provide some information on that? e.g. amount of time and total volume change. My box density goes from ~1.0 to .5 in 5 ps with a compressibility of 5E-05. It goes from ~1.0 to .94 in 300 ps with a compressibility of 5E-06. In both case the slope of density(t) is negative and never levels off. Chris -- original message -- Hi I am having some pressure coupling issues. I have a fairly large protein/water system 400K+ atoms. It minimizes just fine (F 1000). If I run NVE it conserves energy with appropriate parameter settings. If I run NVT it is stable. When I turn on Pcoupl (i.e. Berendsen or Parinello Rahman), the system just continuously expands. My parameters are as follows. Any ideas? Best, Ilya ; ; File 'mdout.mdp' was generated ; By user: relly (508) ; On host: master.simprota.com http://master.simprota.com ; At date: Fri Mar 6 20:17:33 2009 ; ; VARIOUS PREPROCESSING OPTIONS ; Preprocessor information: use cpp syntax. ; e.g.: -I/home/joe/doe -I/home/mary/hoe include = ; e.g.: -DI_Want_Cookies -DMe_Too define = ; RUN CONTROL PARAMETERS integrator = md ; Start time and timestep in ps tinit= 0 dt = 0.004 ;nsteps = 25 nsteps = 250 ; For exact run continuation or redoing part of a run ; Part index is updated automatically on checkpointing (keeps files separate) simulation_part = 1 init_step= 0 ; mode for center of mass motion removal comm_mode= linear ; number of steps for center of mass motion
Re: [gmx-users] pressure coupling for bilayers
semiisotropic is more indicated for bilayers. pragya chohan [EMAIL PROTECTED] wrote: hi i want to rum a NPT ensemble on my bilayer+protein system. What is the best pcoupletype to use. I have seen some posts on gmx user group advising anisotriopic type with tau_p 5.0 ps cpmpressibility: 4.5e-5 4.5e-5 4.5e-5 0.0 0.0 0.0 Another suggestion is welcome. Thanks in advance Pragya _ Search from any Web page with powerful protection. Get the FREE Windows Live Toolbar Today! http://toolbar.live.com/?mkt=en-in___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php - XAvier Periole - PhD 1- Institute of Molecular Assemblies City University of New York - USA 2- Molecular Dynamics-Group University of Groningen - The Netherlands http://md.chem.rug.nl/~periole - ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Pressure Coupling
Thank you very much for the reply, and I'll try to better search the board first :) Cheers, George Dallas B. Warren wrote: Issue like this has been discussed previously, try searching the emailing list for a solution first before posting. Means you get your answer faster :-) Pressure is a macroscopic property that is being monitored and adjusted on a microscopic scale. A variation of that order is entirely normal. How much and quickly it varies depends on what type of coupling you use and the coupling constants. My own simulations oscillate between roughly 400 bar and -400 bar. Catch ya, Dr. Dallas Warren Lecturer Department of Pharmaceutical Biology and Pharmacology Victorian College of Pharmacy, Monash University 381 Royal Parade, Parkville VIC 3010 [EMAIL PROTECTED] +61 3 9903 9524 - When the only tool you own is a hammer, every problem begins to resemble a nail. ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Pressure coupling question
toma0052 wrote: Hello, I am simulating a lipid bilayer system, and am having some trouble understanding the pressure output. When I run a simulation, even when I do not add any perturbations or fix any atoms, the pressure oscillations are quite large. The temperature coupling seems fine. After about 1ps, the system is near the reference temperature, and the oscillations are only about 1 or 2 degrees. So, I expected something, roughly, similar with the pressure. For the pressure, however, my reference pressure is 1.0 bar, but the pressure in the system after a few picoseconds seems to range from -400 to 400 bar. Is this normal? Do I just need to wait longer? Does this mean that I did not run the energy minimization long enough and there are some high forces? Is there something else that I am doing wrong? In my mdp file, the pressure coupling looks like: Pcoupl = berendsen Pcoupltype = isotropic tau_p= 0.5 compressibility = 4e-5 ref_p= 1.0 I have also tried Parrinello-Rahman pressure coupling with the same result. Thanks Mike ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php this is normal. increase your system size to reduce the fluctuations. scales as sqrt(n) :(. -- David. David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group, Dept. of Cell and Molecular Biology, Uppsala University. Husargatan 3, Box 596, 75124 Uppsala, Sweden phone: 46 18 471 4205 fax: 46 18 511 755 [EMAIL PROTECTED] [EMAIL PROTECTED] http://folding.bmc.uu.se ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Re: [gmx-users] Pressure coupling appears to produce errors on the P655+ aix 5.2
Arthur Roberts wrote: Hi, All, I would like to thank Carsten Kutzner, Erik Lindahl, and David van der Spoel for their advise with my previous problem with not getting any data with mdrun. It turned out to be due to the fourier spacing and the PME order. In addition, David van der Spoel's suggestion of using a water shell was quite useful for troubleshooting. I have a new problem: When pressure coupling is used with mdrun, I get the following error: Step 209 Warning: pressure scaling more than 1%, mu: 1.02233 1.02233 1.02233 I tried lengthening the tau_p to 4, but it didn't help as recommended for a similar problem in the mail archives. is the system well minimized and equilibrated? you can try running without pressure coupling for a while. then you can try running without PME for testing only.. Details of the problem can be found here: http://cetus.mchem.washington.edu/pub/supercomputer/Pressure%20causes%20errors.html Best wishes, Art University of Washington ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php -- David. David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group, Dept. of Cell and Molecular Biology, Uppsala University. Husargatan 3, Box 596, 75124 Uppsala, Sweden phone: 46 18 471 4205 fax: 46 18 511 755 [EMAIL PROTECTED] [EMAIL PROTECTED] http://folding.bmc.uu.se ___ gmx-users mailing listgmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please don't post (un)subscribe requests to the list. Use the www interface or send it to [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
