[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