[gmx-users] Re: energy conservation / frozen atoms
Justin Lemkul
jalemkul at vt.edu
Mon Aug 5 13:45:23 CEST 2013
On 8/4/13 4:08 PM, S. Alireza Bagherzadeh wrote:
>> On 8/2/13 3:19 PM, S. Alireza Bagherzadeh wrote:
>>> Thanks for your notes.
>>>
>>>
>>> I did a diagnosis test which could be of relevance here.
>>>
>>> I set up the following system:
>>> [ gas | liquid water (solid water) liquid water | gas ]
>>>
>>> gas is united atom methane.
>>> liquid water is tip4p-ice model and solid water is a cage-like
>> crystalline
>>> structure of water and methane called gas hydrate.
>>>
>>> Now, in order to test the effect of freezing and position restraining on
>>> the performance of nve I did two tests at 370 K.
>>>
>>> Test 1 (freezing):
>>> Solid water was kept frozen in all 3 dimensions (Y Y Y).
>>> First I ran a nvt for 250 ps for equilibration (potential and total
>> energy
>>> both converged after 250 ps, Pressure equilibrated at ~ 3950 bar). Then I
>>> started a 1ns nve.
>>> Similar to my other simulation, the total energy linearly decreased
>> (0.84%
>>> per ns) as well as potential energy. Pressure remained around 3950 bar;
>>> however, the temperature decreased from 370 to 364 K (physically, this
>>> should not happen).
>>>
>>>
>>> Test 2 (position restraining):
>>> Oxygen of solid water was strongly restrained to a point (fc of 100000).
>>> Similar to the previous test, first I ran a nvt for 250 ps for
>>> equilibration (potential and total energy both converged after 250 ps,
>>> Pressure equilibrated at about 0 bar with fluctuations of ~ 2000 bar).
>> Then
>>> I started a 1ns nve.
>>> Again, similar to test 1, the total energy linearly decreased (1.33% per
>>> ns) as well as potential energy. Pressure remain around 0 bar; however,
>> the
>>> temperature initially dropped from 370 K to 355K within 1 ps, then
>>> increased to 358 K during the next 50 ps and thereafter kept linearly
>>> decreasing to 353 K until the end of 1 ns run (physically and
>> intuitively,
>>> this should not happen).
>>>
>>>
>>>
>>> (In both of the tests, I kept the methane inside the cages of solid
>> water
>>> position-restrained to a point by fc = 1000).
>>>
>>> If needed I can post the .mdp and .top files too.
>>>
>>
>> An .mdp file would be useful, otherwise a demonstration that these
>> parameters
>> actually produce an energy-conserving NVE ensemble for a simple system.
>>
>> -Justin
>>
>> --
>>
>
> Here is the .mdp file for freeze test:
>
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; Run control
> integrator = md ; Leap-frog algorithm
> tinit = 0 ; starting time [ps]
> dt = 0.001 ; time step [ps]
> nsteps = 1000000 ; number of steps
> nstcomm = 100 ; frequency for center of mass motion
> removal [steps]
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; Output control
> nstxout = 0 ; frequency to write coordinates to
> output trajectory file [steps]
> nstvout = 0 ; frequency to write velocities to
> output trajectory file [steps]
> nstfout = 0 ; frequency to write forces to output
> trajectory file [steps]
> nstlog = 500 ; frequency to write energies to log
> file [steps]
> nstenergy = 500 ; frequency to write energies to energy
> file [steps]
> nstxtcout = 0 ; frequency to write coordinates to xtc
> trajectory [steps]
> xtc-precision = 1000 ; precision to write to xtc trajectory
> [real]
> xtc_grps = HYDW HYDG SOL GAS
> energygrps = HYDW HYDG SOL GAS
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; Neighborsearching and short-range nonbonded interactions
> nstlist = 10 ; frequency to update the neighbor list
> [steps]
> ns_type = grid ; (grid / simple) search for
> neighboring list
> pbc = xyz ; priodic boundary conditions (xyz / no
> / xy)
> rlist = 1.7 ; cut-off distance for the short-range
> neighbor list [nm]
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; Electrostatics
> coulombtype = PME-Switch
> rcoulomb_switch = 1.3 ; where to switch the Coulomb potential
> [nm]
> rcoulomb = 1.5 ; distance for the Coulomb cut-off [nm]
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; van der Waals
> vdw-type = shift ; (the LJ normal out at rvdw_switch to
> reach zero at rvdw)
> rvdw-switch = 1.3 ; where to strat switching the LJ
> potential [nm]
> rvdw = 1.5 ; cut-off distance for vdw potenrial
> [nm]
> DispCorr = EnerPres ; (Apply long range dispersion
> corrections for Energy and Pressure)
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; EWALD/PME/PPPM parameters
> fourierspacing = 0.12 ; maximum grid spacing for the FFT when
> using PPPM or PME [nm]
> pme_order = 6 ; interpolation order for PME
> ewald_rtol = 1e-06 ; relative strength of the
> Ewald-shifted direct potential at rcoulomb
> epsilon_surface = 0 ; dipole correction to the Ewald
> summation in 3d
> optimize_fft = yes ; optimal FFT plan for the grid at
> startup (yes / no)
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; Temperature coupling
> tcoupl = no
> ;tcoupl = nose-hoover
> nh-chain-length = 1
> tc_grps = system
> tau_t = 0.2 ; time constatn for coupling [ps], 1
> for each group
> ref_t = 370 ; refernce temperature [K], 1 for each
> group
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; Pressure coupling
> Pcoupl = no
> ;Pcoupl = Parrinello-Rahman
> tau_p = 0.5 ; time constant for coupling
> compressibility = 4.5e-05
> ref_p = 40.0 ; reference pressure for coupling [bar]
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; Velocity generation
> gen_vel = no ; generating velocities according to
> maxwell distribution?
> gen_temp = 370 ; [K]
> gen_seed = -1 ; initialize random generator based on
> the process ID number [integer]
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; Bonds
> constraints = none
> constraint-algorithm = shake
> continuation = no ; DO apply constraints to the start
> configuration
> shake-tol = 1e-10
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; Non-equilibrium MD
> freezegrps = HYDW
> freezedim = Y Y Y
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
>
>
>
> Here is the .top file for freeze test:
>
> ; Topology for methane hydrate between silica surfaces in contact with
> water and gas
>
> #include "/home/alireza/gromacs/myff.ff/forcefield.itp"
>
> ;Hydrate
> ;-----------------------------------------------
> ; water topology - hydrate phase
> #include "/home/alireza/gromacs/myff.ff/tip4p-ice_hyd.itp"
>
>
> ; methane topology - hydrate phase
> #include "/home/alireza/gromacs/myff.ff/UAmethane_hyd.itp"
>
> [position_restraints]
> ; ai funct fc
> ; 1 1 1000 1000 1000 ; Restrain Carbon to a point
>
>
> ;SOL
> ;-----------------------------------------------
> ; water topology - liquid phase
> #include "/home/alireza/gromacs/myff.ff/tip4p-ice_sol.itp"
>
> ;------------------------------------------------
>
>
> ;GAS
> ;-----------------------------------------------
> ; methane topology - gas phase
> #include "/home/alireza/gromacs/myff.ff/UAmethane_gas.itp"
>
> ;------------------------------------------------
>
>
>
> [ system ]
> ; Name
> Methane hydrate in contact with water and gas
>
> [ molecules ]
> ; Compound #mols
> HydrateW 2322
> HydrateG 357
> SOL 8659
> MethaneG 133
>
>
> I am also running an nve in which I unfroze the solid water phase after
> equilibration in nvt and it seems that the total energy is very well
> conserved.
>
This may be buggy behavior, so you're welcome to file an issue on
redmine.gromacs.org that provides a description of the problem and necessary
input files to reproduce the issue. It's still not clear to me whether frozen
groups (which are considered a non-equilibrium perturbation) should even truly
be compatible with NVE. In such cases, you have inelastic collisions with the
frozen groups, which should likely lead to a decrease in kinetic energy. I know
you said DL_POLY does the simulations just fine, but I'm not familiar with that
code, so I don't know if there are any tricks that are being played under such
circumstances.
-Justin
--
==================================================
Justin A. Lemkul, Ph.D.
Postdoctoral Fellow
Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201
jalemkul at outerbanks.umaryland.edu | (410) 706-7441
==================================================
More information about the gromacs.org_gmx-users
mailing list