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 <[email protected]>
Subject: Re: [gmx-users] Re: Inquiry about a completely user defined
force field
To: Discussion list for GROMACS users <[email protected]>
Message-ID: <[email protected]>
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 <[email protected]>
Subject: Re: [gmx-users] Inquiry about a completely user defined force
field
To: Discussion list for GROMACS users <[email protected]>
Message-ID:
<
camnumarsfd+wrwqrcpujipsvczypz8kujvaxd9h4ko51zlb...@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 <
[email protected]> 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 [email protected]
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 [email protected].
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
------------------------------
--
gmx-users mailing list
[email protected]
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
========================================
--
gmx-users mailing list [email protected]
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 [email protected].
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
