[gmx-users] Pressure coupling problem

[Fabian Casteblanco]

Thank you Justin and Peter for your responses.  I tried extending the time
on the npt equilibration.  It helped but not much.  My final pressure after
the MD run was about 1.15 bar compared to ref_p which was set to 1 bar.
Peter, I will try to analyze the potential error using RMSD and Drift.

Thanks.

On Mon, Apr 11, 2011 at 4:55 PM, Fabian Casteblanco <
fabian.castebla...@gmail.com> 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.
>
> 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  ;als

Re: [gmx-users] Pressure coupling problem

[Peter C. Lai]

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
> > 

[gmx-users] Pressure coupling problem

[Fabian Casteblanco]

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.

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
;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
;Temperatur

Re: [gmx-users] Pressure coupling problem

[Justin A. Lemkul]

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
;Neighborh

[gmx-users] Pressure coupling problem

[Fabian Casteblanco]

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.

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
;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
;Temp

Re: [gmx-users] Pressure Coupling Problem

[Lucio Montero]

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,  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 wrote:

  [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

[Joe Joe]

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,  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  >wrote:
>
> [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] Pressure Coupling Problem

[chris . neale]

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 wrote:

[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 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

[Joe Joe]

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 wrote:

> 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
> the system (425K atoms) is quite stable. The  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.
>
> ___
> 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

[gmx-users] Pressure Coupling Problem

[Chris Neale]

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
the system (425K atoms) is quite stable. The  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.

___
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

[Joe Joe]

On Thu, Apr 9, 2009 at 6:36 AM, Justin A. Lemkul  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  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

[Justin A. Lemkul]

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 
 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)?  What are your box dimensions doing?  To get that kind of sudden 
change in density, your box dimensions would have to expand astronomically?


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?  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?


-Justin


;
;   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 t

Re: [gmx-users] Pressure Coupling Problem

[Joe Joe]

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  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 nstl

Re: [gmx-users] Pressure Coupling Problem

[Joe Joe]

On Wed, Apr 8, 2009 at 11:31 AM, Roland Schulz  wrote:

>
>
> On Wed, Apr 8, 2009 at 7:53 AM, Joe Joe  wrote:
>
>> HI Chris,
>>
>> On Tue, Apr 7, 2009 at 9:31 PM,  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
>

Yes, my point was that when electrostatics = PME then Gromacs ignores  the
rcoulomb-switch parameter.

>
>
>
>>
>>>
>>> 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  

Re: [gmx-users] Pressure Coupling Problem

[Joe Joe]

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  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>> 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, >> 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)
>>

Re: [gmx-users] Pressure Coupling Problem

[Joe Joe]

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,  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,  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,  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   

Re: [gmx-users] Pressure Coupling Problem

[Roland Schulz]

On Wed, Apr 8, 2009 at 7:53 AM, Joe Joe  wrote:

> HI Chris,
>
> On Tue, Apr 7, 2009 at 9:31 PM,  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 

Re: [gmx-users] Pressure Coupling Problem

[Justin A. Lemkul]

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, > 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, 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 
;   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 co

[gmx-users] Pressure Coupling Problem

[chris . neale]

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,  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,  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 algori

Re: [gmx-users] Pressure Coupling Problem

[Joe Joe]

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,  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,  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 

[gmx-users] Pressure Coupling Problem

[chris . neale]

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,  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_

Re: [gmx-users] Pressure Coupling Problem

[Joe Joe]

HI Chris,

On Tue, Apr 7, 2009 at 9:31 PM,  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 

[gmx-users] Pressure Coupling Problem

[Joe Joe]

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
; Pressure coupling
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
; Scaling of reference coordinates, No, All or COM
refcoord_scaling = No
; Random seed for Andersen thermostat
andersen_seed= 815131
___
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] Pressure Coupling Problem

[chris . neale]

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

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.


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
; Pressure coupling
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
; Scaling of reference coordinates, No, All or COM
refcoord_scaling = No
; Random seed for Andersen thermostat
andersen_seed= 815131
---