[gmx-users] Re: energy conservation / frozen atoms
S. Alireza Bagherzadeh
s.a.bagherzadeh.h at gmail.com
Sun Aug 4 22:08:35 CEST 2013
> 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.
--
Alireza
More information about the gromacs.org_gmx-users
mailing list