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