[gmx-users] Re: Unsteady Density

Justin A. Lemkul jalemkul at vt.edu
Mon May 2 23:02:15 CEST 2011



Fabian Casteblanco wrote:
> 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?
> 

No, in fact it will probably make it worse.  Use the values prescribed by the 
force field.

-Justin

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

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================



More information about the gromacs.org_gmx-users mailing list