[gmx-users] Re: Unsteady Density
Justin A. Lemkul
jalemkul at vt.edu
Mon May 2 21:50:25 CEST 2011
Fabian Casteblanco wrote:
> Hello,
>
> Thanks for your response. Here are my responses to your points:
>
> nstlist=1 is overkill. I will change this back to 5 except for EM
> like you said. I originially thought that maybe the density not
> settling was because the neighbor searching frequency had not been
> updated frequently enough causing molecule density to change.
>
Not likely.
> When I performed this specific NPT run, the average was off due to the
> fact that the density started off in the 400s and rised up to ~high
> 700s. The last time I had performed this with nstlist=5 and for the
> actual MD run, I had a density of approximately 0.781 g/cm3 compared
> to literature value of 0.791 g/cm3 at 298K. I was getting similar
> results for 1-propanol with the difference much greater. I was hoping
> that if I could somehow fix the methanol system, I would be able to
> fix the 1-propanol system.
>
If you're taking averages over the whole time, that probably accounts for the
discrepancy. The initial values are clearly not correct, and thus those frames
should be discarded as equilibration. Block averaging will show you if there
are any systematic trends in the data (i.e. drifts that indicate you're not yet
converged).
> I created the topologies using the OPLSAA *.itp file (to keep the
> structure correctly, bonds, angles, etc) and I reassigned the values
> corresponding to CHARMM (attached). At first, I had placed the values
> straight from the literature found on the CHARMM website,
> http://mackerell.umaryland.edu/CHARMM_ff_params.html, but I realized
> they were nearly identical and were infact the same values in the
> CHARMM.ff files in GROMACS (i.e. ffbonded.itp, ffnonbonded.itp) All I
> did was reassign the functions to the functions used in the
> ffbonded.itp and ffnonbonded.itp found in CHARMM files on GROMACS.
> The charges were all from the CHARMM parameter website that had this
> molecule already documented (a copy of the lines is placed at the top
> of the *.itp file attached).
>
Based on the fact that you're taking the parameters straight from CHARMM, I
don't suspect anything wrong here.
> I guess the compressibility factor might make a difference. The
> manual said that the value is fairly close to most liquids so I had
> simply ignored the difference but I will look into it.
>
You're looking at a difference of 1.3% (or less, depending on whether or not you
include unequilibrated frames in your analysis), so even something that may seem
trivial could cause this difference.
> After reading more up on it, I was thinking maybe accounting for
> gradual cutoffs near the ends using rcoulomb_switch and similarly for
> rvdw.
>
Switching functions have their own issues associated with them (discontinuities
at the longest cutoff), so for now, stick with PME unless you have specific
reason to use a switching function (i.e., that's how the parameters were
designed to be used).
Otherwise, this might all be a matter of incorrect analysis, using frames that
are not supposed to be included.
-Justin
> Anything wrong with the way I created the topology? Dr. Vitaly, I
> also attached the *.itp files like you mentioned.
>
> Thanks again for all your help.
>
>
>
> 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