[gmx-users] Re: Re: Inquiry about a completely user defined force field
S. Alireza Bagherzadeh
s.a.bagherzadeh.h at gmail.com
Fri Feb 1 00:17:51 CET 2013
----------------------------------------------------------------------
>
> Message: 1
> Date: Thu, 31 Jan 2013 16:57:56 -0500
> From: Justin Lemkul <jalemkul at vt.edu>
> Subject: Re: [gmx-users] Re: Inquiry about a completely user defined
> force field
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID: <510AE8E4.10304 at vt.edu>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
>
>
> On 1/31/13 1:57 PM, S. Alireza Bagherzadeh wrote:
> > Dear Dr.Abraham,
> >
> > Thank you for taking the time to look at my files.
> > However, I have one question remained in my mind:
> >
> > Where is the "atomtypes.atp" used in the force field?
> >
> > I suspect that this file is only for clarification for the user in
> > order to eliminate/minimize
> > confusion and this file is not read by the program. Am I correct?
> >
>
> No, the .atp file is used by pdb2gmx. See the manual.
>
Dear Justin,
I am not using pdb2gmx in any steps of preparation of my input files. I
have written a short script myself that converts my .xyz files to .gro
files.
So, this (not using the "atomtypes.atp" file) should not be a problem for
me...
> > I would also like to ask you for any insight on why in my system the
> tempera
> > ture keeps increasing linearly with time under NVE ensemble (even up to
> > unrealistic T of 550K)? I have verified this force filed with DL_POLY and
> > my system showed a gradual temperature decrease which is the correct
> > physical behavior of the system.
> >
>
> Sorry, can't comment here. Without an .mdp file, it's hard to say
> anything.
> Apologies if you posted it before and I missed it.
>
Here is my .mdp file. I appreciate it if you could take a look and point
out any possible errors to me:
;-----------------------------
------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; 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 = 1000 ; 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 = 1000 ; frequency to write coordinates to
xtc trajectory [steps]
xtc-precision = 1000 ; precision to write to xtc trajectory
[real]
xtc_grps = HYDW HYDG SOL
energygrps = HYDW HYDG SOL SiO2 SiOH
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; 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.0 ; cut-off distance for the short-range
neighbor list [nm]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Electrostatics
coulombtype = PME
rcoulomb = 1.0 ; distance for the Coulomb cut-off [nm]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; van der Waals
vdw-type = switch ; (the LJ normal out at rvdw_switch to
reach zero at rvdw)
rvdw-switch = 0.8 ; where to strat switching the LJ
potential [nm]
rvdw = 0.9 ; 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 = no ; optimal FFT plan for the grid at
startup (yes / no)
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Temperature coupling
tcoupl = nose-hoover
tc_grps = SiO2 SiOH SOL HYDG HYDW
tau_t = -1.0 -1.0 -1.0 -1.0 -1.0 ; time constant [ps]
ref_t = 300 300 300 300 300 ;
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Pressure coupling
Pcoupl = no
;Pcoupl = Parrinello-Rahman
tau_p = 0.5 ; time constant for coupling
compressibility = 4.5e-05
ref_p = 100.0 ; reference pressure for coupling
[bar]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Velocity generation
gen_vel = no ; generating velocities according to
maxwell distribution
gen_temp = 300 ; [K]
gen_seed = -1 ; initialize random generator based on
the process ID number [integer]
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Energy group exclusions
energygrp_excl = SiOH SiOH SiO2 SiO2 SiOH SiO2
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Bonds
constraints = all-angles
constraint-algorithm = shake
;continuation = yes ; do not apply constraints to the
start configuration
shake-tol = 1e-5
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
; Non-equilibrium MD
freezegrps = SiOH SiO2
freezedim = Y Y Y Y Y Y
;-----------------------------------------------------------------------------------------------------;
;-----------------------------------------------------------------------------------------------------;
Please see the energy trends and T profile for the initial 0.5 NVT
equilibration here (of course, I am using a different mdp file for this
step. All of my residues are thermally coupled to a bath at 300K):
http://s1275.beta.photobucket.com/user/SAlirezaB/media/gromacs%20related%20stuff/05ns_NVTequilibirationwithHYDWfrozen_zps57247dbe.png.html<http://i1275.photobucket.com/albums/y460/SAlirezaB/gromacs%20related%20stuff/05ns_NVTequilibirationwithHYDWfrozen_zps57247dbe.png>
The same trends for NVE production run (above mdp file) can bee seen here:
http://s1275.beta.photobucket.com/user/SAlirezaB/media/gromacs%20related%20stuff/1ns_NVEwithHYDWunfrozen_zps31539555.png.html<http://i1275.photobucket.com/albums/y460/SAlirezaB/gromacs%20related%20stuff/1ns_NVEwithHYDWunfrozen_zps31539555.png>
> > One clue which might be useful is that even though when I hand-calculate
> > the total charge of my system to be zero, the gormacs warns me that my
> > system has the total charge of 0.000108. Could this be the reason? Do you
> > have any suggestions for me to solve this potential problem? (I am not
> > using gromacs in double-precision mode).
> >
>
> Probably just a normal consequence of floating-point math. Check the
> topology
> to be sure.
>
>
> http://www.gromacs.org/Documentation/Errors#System_has_non-zero_total_charge
>
> -Justin
>
>
I believe 0.000108 is close enough to an integer (0) and therefore as
mentioned in that page it is just rounding error and should not be a major
problem.
Thank you,
> >
> > Best regards,
> > Alireza
> >
> >
> >
> >
> >
> >>
>
> >> ------------------------------
> >>
> >> Message: 2
> >> Date: Mon, 28 Jan 2013 11:59:37 +0100
> >> From: Mark Abraham <mark.j.abraham at gmail.com>
> >> Subject: Re: [gmx-users] Inquiry about a completely user defined force
> >> field
> >> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> >> Message-ID:
> >> <
> >> CAMNuMARsFD+WrWqrcPujiPsVcZYPZ8KuJVaxD9h4kO51zLBzrg at mail.gmail.com>
> >> Content-Type: text/plain; charset=ISO-8859-1
> >>
> >> There's nothing that leaps out at me as being wrong - but syntax
> checking
> >> is the job of the parser (grompp) and logic checking can only be done by
> >> you. To do that, I would proceed slowly, demonstrating that small
> chunks of
> >> the force field can be used to generate results you can validate against
> >> something external (e.g. experiment, or a reference implementation of
> the
> >> force field). Don't just grab your simulation system and assume
> everything
> >> works!
> >>
> >> Mark
> >>
> >> On Mon, Jan 28, 2013 at 7:35 AM, S. Alireza Bagherzadeh <
> >> s.a.bagherzadeh.h at gmail.com> wrote:
> >>
> >>> Dear gmx-users,
> >>>
> >>> I am simulating a system in which I have a force field that is not
> >> included
> >>> in gromacs forcefields.
> >>>
> >>> So, I decided to construct my own force field in a sub-directory of my
> >>> actual
> >>> run.
> >>>
> >>> I read chapter 5 of gromacs 4.5.4 user manual (pgs. 107-139) and I
> >> prepared
> >>> the following files.
> >>> Below I am copying the content of a few sample files as well as my
> >> topology
> >>> file.
> >>>
> >>> I appreciate it if you could take a look and advise me whether anything
> >> is
> >>> missing, wrong or incomplete.
> >>>
> >>> Best regards,
> >>>
> >>> *{topol.top}:*
> >>> ; Topology for methane hydrate between silica surfaces in contact with
> >>> water and gas
> >>>
> >>> #include
> >>>
> >>>
> >>
> "/home/alireza/Silica_H2O_Gas_Hyd_Decomp/topol_preparation/myff.ff/forcefield.itp"
> >>>
> >>>
> >>> ;Hydrate
> >>> ;-----------------------------------------------
> >>> ; water topology - hydrate phase
> >>> #include
> >>>
> >>>
> >>
> "/home/alireza/Silica_H2O_Gas_Hyd_Decomp/topol_preparation/myff.ff/spce_hyd.itp"
> >>>
> >>> ; methane topology - large cages of hydrate phase
> >>> #include
> >>>
> >>>
> >>
> "/home/alireza/Silica_H2O_Gas_Hyd_Decomp/topol_preparation/myff.ff/UAmethane_lc.itp"
> >>>
> >>> ; methane topology - small cages of hydrate phase
> >>> #include
> >>>
> >>>
> >>
> "/home/alireza/Silica_H2O_Gas_Hyd_Decomp/topol_preparation/myff.ff/UAmethane_sc.itp"
> >>> ;------------------------------------------------
> >>>
> >>>
> >>> ; water topology - liquid water
> >>> #include
> >>>
> >>>
> >>
> "/home/alireza/Silica_H2O_Gas_Hyd_Decomp/topol_preparation/myff.ff/spce_liq.itp"
> >>>
> >>>
> >>>
> >>> ;Silica
> >>> ;----------------------------
> >>> ; silica topology - bulk Si
> >>> #include
> >>>
> >>>
> >>
> "/home/alireza/Silica_H2O_Gas_Hyd_Decomp/topol_preparation/myff.ff/SilicaLopez_Si.itp"
> >>>
> >>> ; silica topology - bulk O
> >>> #include
> >>>
> >>>
> >>
> "/home/alireza/Silica_H2O_Gas_Hyd_Decomp/topol_preparation/myff.ff/SilicaLopez_Oi.itp"
> >>>
> >>> ; silica topology - surface O
> >>> #include
> >>>
> >>>
> >>
> "/home/alireza/Silica_H2O_Gas_Hyd_Decomp/topol_preparation/myff.ff/SilicaLopez_Os.itp"
> >>>
> >>> ; silica topology - surface H
> >>> #include
> >>>
> >>>
> >>
> "/home/alireza/Silica_H2O_Gas_Hyd_Decomp/topol_preparation/myff.ff/SilicaLopez_Hs.itp"
> >>> ;-----------------------------
> >>>
> >>> [ system ]
> >>> ; Name
> >>> Methane hydrate between 2 silica surfaces in contact with water and gas
> >>>
> >>> [ molecules ]
> >>> ; Compound #mols
> >>> HydrateW 2322
> >>> HydG_sc 81
> >>> HydG_lc 276
> >>> SOL 3413
> >>> Silica_Si 1848
> >>> Silica_Oi 3360
> >>> Silica_Os 672
> >>> Silica_Hs 672
> >>>
> >>>
> >>>
> >>> *SAMPLE FORCE FIELD FILES:*
> >>>
> >>> *{atomtypes.atp}:*
> >>> ;atmnm mass
> >>> Description
> >>> Reference
> >>> O 15.99940 ; water oxygen (hydrate & liquid phase),
> >>> tip4p-ice :
> >>> Ospc 15.99940 ; water oxygen (hydrate & liquid phase),
> >>> spc/e :
> >>> H 1.00800 ; water hydrogen (hydrate & liquid phase),
> >>> tip4p-ice :
> >>> Hspc 1.00800 ; water hydrogen (hydrate & liquid phase),
> >>> spc/e :
> >>> E 0.00000 ; dummy atom
> >>> (tip4p-ice) :
> >>> CH4 16.04300 ; United Atom methane (small/large cages of hydrate
> >> phase &
> >>> gas phase):
> >>> Si 28.08550 ; silica silicon (quartz
> >>> phase) :
> >>> Oi 15.99940 ; silica oxygen, bulk (inside) of the
> >>> crystal :
> >>> Os 15.99940 ; silica oxygen, surface of crystal (hydroxyl
> >>> group) :
> >>> Hs 1.00800 ; silica hydrogen, surface of crystal (hydroxyl
> >>> group) :
> >>>
> >>> *{ffnonbonded.itp}:*
> >>> [ atomtypes ]
> >>> ; name at.num mass charge ptype sigma(nm) epsilon(kj/mol)
> >>> H 1 1.0080 0.5897 A 0.00000e+00 0.00000e+00 ;
> >>> tip4p-ice: JCP 122,234511(2005)
> >>> Hspc 1 1.0080 0.4238 A 0.00000e+00 0.00000e+00 ;
> >> spce:
> >>> Hs 1 1.0080 0.3200 A 0.04000e+00 0.19246e+00 ;
> >> Lopez
> >>> et al: JPCB 110,2782 (2006)
> >>> CH4 6 16.0430 0.0000 A 0.36400e+00 1.36500e+00 ;
> >> United
> >>> Atom: JPC 87,4198(1983)
> >>> O 8 15.9994 0.0000 A 0.31668e+00 0.88217e+00 ;
> >>> tip4p-ice: JCP 122,234511(2005)
> >>> Ospc 8 15.9994 -0.8476 A 0.31656e+00 0.65017e+00 ;
> >> spc/e:
> >>> Oi 8 15.9994 -0.5300 A 0.31538e+00 0.63639e+00 ;
> >> Lopez
> >>> et al: JPCB 110,2782 (2006)
> >>> Os 8 15.9994 -0.6400 A 0.31538e+00 0.63639e+00 ;
> >> Lopez
> >>> et al: JPCB 110,2782 (2006)
> >>> Si 14 28.0855 1.0800 A 0.39200e+00 2.51040e+00 ;
> >> Lopez
> >>> et al: JPCB 110,2782 (2006)
> >>> E 0 0.0000 -1.1794 D 0.00000e+00 0.00000e+00 ;
> >>> tip4p-ice: JCP 122,234511(2005)
> >>>
> >>> *{forcefield.itp}:*
> >>> [ defaults ]
> >>> ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
> >>> 1 2 no 1.0 1.0
> >>> ;MEANING
> >>> ;--------------------------------------------------------------
> >>> ;L-J lorentz-
> >>> ;interaction berthelott
> >>> ;function combination
> >>> ; rule
> >>>
> >>>
> >>> #include "ffnonbonded.itp"
> >>> ;#include "ffbonded.itp"
> >>>
> >>>
> >>> *{SilicaLopez_Hs.itp}:*
> >>> ; methane topology - small cages of hydrate phase
> >>> [ moleculetype ]
> >>> ; Name nrexcl
> >>> Silica_Hs 1
> >>>
> >>> [ atoms ]
> >>> ; nr type resnr residue atom cgnr charge mass
> >>> 1 Hs 1 SiOH Hs 1 0.3200 1.008
> >>>
> >>>
> >>> *{spce_hyd.itp}:*
> >>> ;spce water for Hydrate phase
> >>> ;
> >>> [ moleculetype ]
> >>> ; molname nrexcl
> >>> HydrateW 2
> >>>
> >>> [ atoms ]
> >>> ; nr type resnr residue atom cgnr charge mass qtot
> >>> 1 Ospc 1 HYDW OWc 1 -0.8476 15.994 ;
> -0.8476
> >>> 2 Hspc 1 HYDW HWc 1 0.4238 1.008 ;
> -0.4238
> >>> 3 Hspc 1 HYDW HWc 1 0.4238 1.008 ; 0.0
> >>>
> >>> [ constraints ]
> >>> ; i j func d (nm)
> >>> 1 2 1 0.10000
> >>> 1 3 1 0.10000
> >>> 2 3 1 0.16330
> >>>
> >>> [ exclusions ]
> >>> 1 2 3
> >>> 2 1 3
> >>> 3 1 2
> >>>
> >>>
> >>> *{UAmethane_gas.itp}:*
> >>> ; methane topology - gas phase
> >>> [ moleculetype ]
> >>> ; Name nrexcl
> >>> MethaneG 1
> >>>
> >>> [ atoms ]
> >>> ; nr type resnr residue atom cgnr charge mass
> >>> 1 CH4 1 GAS Cg 1 0.000 16.043
> >>>
> >>>
> >>>
> >>>
> >>>
> >>
> >>> --
> >>> gmx-users mailing list gmx-users at 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-request at gromacs.org.
> >>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >>>
> >>
> >>
> >> ------------------------------
> >>
> >> --
> >> gmx-users mailing list
> >> gmx-users at 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!
> >>
> >> End of gmx-users Digest, Vol 105, Issue 136
> >> *******************************************
> >>
> >
> >
> >
>
> --
> ========================================
>
> Justin A. Lemkul, Ph.D.
> Research Scientist
> 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