[gmx-users] Re: Re: Inquiry about a completely user defined force field
Justin Lemkul
jalemkul at vt.edu
Fri Feb 1 00:57:33 CET 2013
On 1/31/13 6:17 PM, S. Alireza Bagherzadeh wrote:
> ----------------------------------------------------------------------
>>
>> 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...
>
Nor would I expect it. But since you asked... ;)
>
>
>>> 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:
General points:
http://www.gromacs.org/Documentation/Terminology/NVE
Specific comments below.
> ;-----------------------------
> ------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; 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
nstlist may be fine, consider using -1 as suggested by the link above.
> [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
Again, per the link above, rlist should be larger than max(rcoulomb,rvdw) so you
may need to increase this value.
> 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 ;
I wouldn't do this. Setting tau_t = -1 is a trick one can use to selectively
turn off coupling of a group, but it strikes me as a bit unpredictable. Be safe
and set "tcoupl = no" and guarantee that you should be getting the ensemble you
desire.
-Justin
> ;-----------------------------------------------------------------------------------------------------;
> ;-----------------------------------------------------------------------------------------------------;
> ; 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
>>
>> ========================================
>>
>>
>> ------------------------------
>>
--
========================================
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