[gmx-users] Osmotic pressure

[gozde ergin]

Hi Justin,

I corrected the nonbonded settings as your suggestion ; 

; Bond parameters
continuation         = no   
constraint_algorithm = shake     
constraints          = h-bonds 
shake_tol            = 0.0001
cutoff-scheme = Verlet
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
; Neighborsearching
ns_type     = grid            
nstlist     = 5                 
rlist       = 1.2               
rcoulomb    = 1.2             
rvdw        = 1.2              

However my osmotic pressure result is still far from the Luo&Roux 2010 study.
I have the same number of Na Cl ions and water molecules. I run the simulation for 15ns.
System size is the same, water model and force field are the same as Roux study. 
Basically I should get the same result.
Here is the result,

Osmotic pressure for 5M : My result 217 bar
                                          Roux study ~300 bar

Any suggestions would be appreciated.

> On 15 Apr 2016, at 11:50, Justin Lemkul <jalemkul at vt.edu> wrote:
> 
> 
> 
> On 4/15/16 5:20 AM, gozde ergin wrote:
>> Dear all,
>> 
>> I simulate the NaCl solution to estimate the osmotic pressure. My salt concentrations are 3,4 and 5M. I apply flat-bottom restraint to the molecules.
>> I use CHARMM36 ff with NBFIX correction.
>> 
>> After the simulation I extract the z-coordinates of restraint ions and use the equation P = F/A , F= k(zi-zwall), A=area to estimate the osmotic pressure.
>> Actually I try to get the similar results as Lou&Roux 2010 study.
>> 
>> However my osmotic pressure results are on the line of ideal solution osmotic pressure as shown equation in below;
>> 
>> P = cRT (Van’t Hoff equation)
>> 
>> but not the similar result with experiments.
>> 
> 
> What values do you actually get?  How do they compare with the values from the Roux paper?  Are you remembering to divide by 2*area in your calculation (since you have two walls)?
> 
>> Is there anyone here that get the same trend as me for osmotic pressure calculation? Or is there something that I miss?
>> 
>> 
>> Here is my .mdp file;
>> 
>> define                   = -DPOSRES
>> integrator               = md
>> dt                       = 0.002
>> nsteps                   = 2500000   ;
>> ; Output control
>> nstxout                  = 2000
>> nstvout                  = 2000
>> nstlog                   = 2000
>> nstenergy                = 2000
>> ; Bond parameters
>> continuation         = no        ;
>> constraint_algorithm = shake     ; h
>> constraints          = all-bonds ; a
>> shake_tol            = 0.0001
>> ; Neighborsearching
>> ns_type     = grid              ;
>> nstlist     = 5                 ;
>> rlist       = 1.1               ;
>> rcoulomb    = 1.1               ;
>> rvdw        = 1.1               ;
> 
> These nonbonded settings are wrong.  The values for CHARMM36 are well established and you should not deviate from them.
> 
> http://www.gromacs.org/Documentation/Terminology/Force_Fields/CHARMM <http://www.gromacs.org/Documentation/Terminology/Force_Fields/CHARMM>
> 
> -Justin
> 
>> ; Electrostatics
>> coulombtype     = PME           ;
>> pme_order       = 4             ;
>> fourierspacing  = 0.16          ;
>> 
>> tcoupl                   = berendsen
>> tc-grps                  = System
>> tau_t                    = 1.0
>> ref_t                    = 300
>> ; Pressure coupling is on
>> pcoupl              = Berendsen     ; Pressure coupling on in NPT, also weak coupling
>> pcoupltype          = semiisotropic     ; uniform scaling of x-y-z box vectors
>> tau_p               = 2.0 2.0          ; time constant, in ps
>> ref_p               = 1.0 1.0          ; reference pressure (in bar)
>> compressibility     = 0 4.5e-5        ; isothermal compressibility, bar^-1
>> refcoord_scaling    = com
>> ; Periodic boundary conditions
>> pbc     = xyz                   ; 3-D PBC
>> ; Dispersion correction
>> DispCorr    = EnerPres          ; account for cut-off vdW scheme
>> ; Velocity generation
>> gen_vel     = yes               ; Velocity generation is on
>> gen_temp    = 300               ; temperature for velocity generation
>> gen_seed    = -1                ; random seed
>> ; COM motion removal
>> ; These options remove COM motion of the system
>> nstcomm         = 10
>> comm-mode       = Linear
>> comm-grps       = System
>> 
>> 
> 
> -- 
> ==================================================
> 
> Justin A. Lemkul, Ph.D.
> Ruth L. Kirschstein NRSA Postdoctoral Fellow
> 
> Department of Pharmaceutical Sciences
> School of Pharmacy
> Health Sciences Facility II, Room 629
> University of Maryland, Baltimore
> 20 Penn St.
> Baltimore, MD 21201
> 
> jalemkul at outerbanks.umaryland.edu <mailto:jalemkul at outerbanks.umaryland.edu> | (410) 706-7441
> http://mackerell.umaryland.edu/~jalemkul <http://mackerell.umaryland.edu/~jalemkul>
> 
> ==================================================
> -- 
> Gromacs Users mailing list
> 
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List <http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List> before posting!
> 
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists <http://www.gromacs.org/Support/Mailing_Lists>
> 
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users <https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users> or send a mail to gmx-users-request at gromacs.org <mailto:gmx-users-request at gromacs.org>.