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
jalem...@outerbanks.umaryland.edu | (410) 706-7441
==================================================
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
* Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
