[gmx-users] Re: Unsteady Density

[Fabian Casteblanco]

Thanks Justin.

I'm currently running the MD run separately so that the average value
won't be affected by the unequilibrated values at the beginning of the
NPT run where the density started in the 400s and ran up to the ~high
700s.  I had done this before when I had nstlist at 5 but I did run
NPT for a lot shorter time last time so Im hoping this time it will be
better since last time it was still off (0.781 g/cm3 to 0.791 g.cm3).
I will also change the compressibility.  Do you think increasing the
rlist, rcoulomb, rvdw value to about 1.6, 1.8 nm might also work a bit
better?

Thanks again.

On Mon, May 2, 2011 at 2:16 PM, Fabian Casteblanco
<fabian.casteblanco at gmail.com> wrote:
> Hello,
>
> I have been trying to use CHARMM forcefield to simulate 1000 molecules
> of Methanol in a box but I seem to fall short on the density every
> time I run.  At first, I had my rlist at 1 but CHARMM recommends
> greater than 1.2-1.4 nm so I changed accordingly, then I changed it so
> that the short neighbor list would update more freqently and finally I
> even added a greater amount of steps, hoping that the pressure and
> density would settle down at some value.
>
> Below, I had taken and modified from the Lysozyme tutorial.  One thing
> I noticed is the contraint -algorithm and how the constraints are set
> to 'all-bonds'.  Is this necessary?  After running NVT (leveled off
> near 298K, seemed ok), then NPT, for 300,000 steps, the pressure is
> close to the reference 1 bar, ~1.1,1.2 bar, but I noticed the density
> is still oscillating, not by much, but it still leaves me questioning
> why its not settling around to the literature value of 0.791 g/cm3.  I
> posted the last several lines from the analysis *.xvg file.  NPT
> script is also below.
>
> Is there still more steps that I have to run for?  The density just
> seems to be oscillating too much and although a literature value of
> 0.791 g/cm3 is within the data, it doesnt seem to be settling around
> that value.  I have similar problems for 1-propanol to where my
> simulated density is under that of the literature value.  Could it be
> possible that maybe I should use different NVT, NPT algorithms
> (tcoupl, pcoupl)?
>
> I would really appreciate anyones help. I'm still in my first few
> months of learning the program.  Thanks!
>
>  591.800000  38039.125000  297.409546  -206.071747  779.428894
>  592.000000  38035.421875  294.951111  -20.557419  777.816956
>  592.200000  37982.812500  298.063324  -463.192047  776.430542
>  592.400000  37681.421875  295.430878   36.529854  776.460876
>  592.600000  37476.148438  297.193787  -28.961288  775.634644
>  592.800000  37611.574219  297.059875  -553.369263  774.096924
>  593.000000  37682.277344  288.297974  -293.567963  776.055237
>  593.200000  37147.824219  295.985901  -221.989441  781.570679
>  593.400000  37404.906250  293.938721  399.093018  788.921570
>  593.600000  37339.914062  297.464783  415.891235  793.838440
>  593.800000  37909.019531  292.707184  574.818726  796.730286
>  594.000000  37369.820312  304.095703  -119.569496  796.235962
>  594.200000  37293.445312  294.449005  -63.335049  794.578552
>  594.400000  37196.917969  295.323761  104.148239  791.275085
>  594.600000  37776.199219  298.004333   28.989655  785.244629
>  594.800000  37893.097656  302.397308  -268.065765  779.457886
>  595.000000  38123.460938  295.343109  -140.914047  775.682678
>  595.200000  38002.054688  302.233124  -91.893478  774.813660
>  595.400000  37835.652344  297.808319  441.931885  776.142090
>  595.600000  38103.964844  296.369019  490.766754  778.422302
>  595.800000  37686.902344  303.358093  -19.726471  781.057922
>  596.000000  37586.890625  293.751648  -92.387199  783.724854
>  596.200000  37408.414062  300.007385  -288.156036  785.990417
>  596.400000  37894.898438  295.720520  346.861206  788.753174
>  596.600000  37673.378906  288.311432   89.059952  789.179993
>  596.800000  36881.062500  295.265564   70.789131  791.552673
>  597.000000  37077.789062  296.261444   72.207787  793.970398
>  597.200000  37123.230469  296.511475  -232.434616  793.653564
>  597.400000  37426.152344  304.281891  -133.300797  791.791077
>  597.600000  37514.332031  305.183197  -259.364136  787.592712
>  597.800000  37393.667969  301.326111  -150.162338  785.053711
>  598.000000  37263.054688  299.811554  326.851837  785.553894
>  598.200000  37203.261719  300.688904  225.526001  788.008667
>  598.400000  37511.636719  302.506714   66.299194  790.725281
>  598.600000  37460.843750  299.512909  198.943665  793.353333
>  598.800000  37248.921875  291.575134  -93.283127  793.582703
>  599.000000  36870.480469  293.515381  211.069107  792.929443
>  599.200000  37398.097656  293.184448  -163.371140  789.150513
>  599.400000  37402.238281  297.402008   45.650658  786.303711
>  599.600000  37419.816406  298.520508  -147.688004  783.492371
>  599.800000  37497.332031  296.529877  -15.991615  782.294250
>  600.000000  37691.468750  298.839844  -68.580627  781.974182
>
>
>
> title           =Methanol npt equilibration
>
> ;Run parameters
> integrator      =md                     ;leap-frog algorithm
> nsteps  =300000         ;2 * 50000 = 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         =1              ;2 fs
> rlist                   =1.4            ;short-range neighborlist cutoff (in nm)
> rcoulomb                =1.4            ;short-range electrostatic cutoff (in nm)
> rvdw                    =1.4            ;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                 ;one coupling group
> tau_t                   =0.1                          ;time constant, in ps
> ref_t                   =298                    ;reference temperature, one for
> each ;group, in K
>
> ;Pressure coupling is on
> pcoupl          =Parrinello-Rahman      ;Pressure coupling on in NPT
> pcoupltype              =isotropic                      ;uniform scaling of box vectors
> tau_p           =2.0                            ;time constant, in ps
> ref_p                   =1.0                            ;reference pressure, in bar
> compressibility =4.5e-5                 ;isothermal compressibility of h2O, 1/bar
>
> ;Periodic boundary conditions
> pbc                     =xyz                    ; 3-D PBC
>
> ;Dispersion correction
> DispCorr                =EnerPres               ;account for cut-off vdW scheme
>
> ;Velocity generation
> gen_vel         =no             ;Velocity generation is off
>
>
>
> --
> Best regards,
>
> Fabian F. Casteblanco
> Rutgers University --
> Chemical Engineering PhD Student
> C: +908 917 0723
> E:  fabian.casteblanco at gmail.com
>



-- 
Best regards,

Fabian F. Casteblanco
Rutgers University --
Chemical Engineering PhD Student
C: +908 917 0723
E:  fabian.casteblanco at gmail.com