Will Glover wrote:
Hi,
I have constructed a gromacs MD simulation of liquid tetrahydrofuran (at RTP)
in which I'm treating each molecule as a charge group. When using a Switch for
both the Coulomb and VdW interactions I see artifacts in the radial
distribution function of the Oxygen--alpha-carbon (and other pair) distances
around the switching region.
I understand that these artifacts occur in atom-based truncation schemes and of course PME fixes the problem (see e.g. http://doi.wiley.com/10.1002/jcc.10117); however, I would ultimately like to interface Gromacs with my own mixed quantum/classical code for which the Ewald method is highly nontrivial and a simpler group-based truncation would suffice.
Does Gromacs 4.05 only use atom-based truncation, even when charge groups are
defined? I was unable to clearly determine this from the manual. If so, is
there a work-around to achieve group-based switching that does not involve
hacking the code?
Gromacs uses charge-group based neighborsearching: if the centers of the
charge group are within the cutoff all the interactions are computed.
However since you use a shift function the energies beyond 1.3 nm will
be zero, the shift function is atom based. Would you have selected a
real cut-off then all the atom pairs would be taken into account.
What is to be preferred does however depend on the size of your
molecule. In this case the molecule is roughly 0.4 nm. My general advice
is to make charge groups small, such that you have a homegeneous system.
In your case the oxygen with to flanking CH2 groups could be one, and
the other two CH2 groups could be in a group each. This way you could
also reduce rlist slightly.
By the way, if you use nstlist = 1, neighborsearching is done each step,
and hence rlist > rcoulomb does not make sense. If you set nstlist = 5
or so you do need rlist > rcoulomb for proper energy conservation.
I have pasted condensed configuration files below for those that are
interested...
Many thanks in advance!
*** mdout.mdp ***
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.004
nsteps = 250000
; For exact run continuation or redoing part of a run
; Part index is updated automatically on checkpointing (keeps files separate)
simulation_part = 1
init_step = 0
; mode for center of mass motion removal
comm_mode = None
; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist = 1
; ns algorithm (simple or grid)
ns_type = grid
; Periodic boundary conditions: xyz, no, xy
pbc = xyz
periodic_molecules = no
; nblist cut-off
rlist = 1.625
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = SWITCH
rcoulomb_switch = 1.1
rcoulomb = 1.3
; Relative dielectric constant for the medium and the reaction field
epsilon_r = 1
epsilon_rf = 1
; Method for doing Van der Waals
vdwtype = Switch
; cut-off lengths
rvdw_switch = 1.1
rvdw = 1.3
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = no
; Extension of the potential lookup tables beyond the cut-off
table-extension = 4
; Seperate tables between energy group pairs
energygrp_table =
; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl = Nose-Hoover
; Groups to couple separately
tc_grps = System
; Time constant (ps) and reference temperature (K)
tau_t = 0.5
ref_t = 298
; Pressure coupling
Pcoupl = no
; OPTIONS FOR BONDS
constraints = all-bonds
; Type of constraint algorithm
constraint_algorithm = LINCS
; Do not constrain the start configuration
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs_order = 6
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs_iter = 2
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs_warnangle = 30
; Convert harmonic bonds to morse potentials
morse = no
*** thf.itp (note, my THF is rigid and planar, so I use a 3-membered ring and
virtual sites to achieve this) ***
[ defaults ]
; nbft cr gp fLJ fQQ
1 2 1 1 1
[ atomtypes ]
; at m q pt V W
O 26.51056525660896 -0.500 A 0.3000 0.7112584
T2 22.79841737169552 0.000 A 0.0000 0.0000000
CH2A 0.000000000000000 0.250 V 0.3800 0.4937137
CH2B 0.000000000000000 0.000 V 0.3905 0.4937137
[ moleculetype ]
; Name nrexcl
THF 3
[ atoms ]
; nr type resnr resid atom cgnr charge mass
1 O 1 THF O1 1 -0.500
2 CH2A 1 THF CA1 1 0.250
3 CH2A 1 THF CA2 1 0.250
4 CH2B 1 THF CB1 1 0.000
5 CH2B 1 THF CB2 1 0.000
6 T2 1 THF T21 1 0.000
7 T2 1 THF T22 1 0.000
[ constraints]
; ai aj funct c0 c1
1 6 1 0.2184048538518238
1 7 1 0.2184048538518238
6 7 1 0.2183193593425659
[ virtual_sites3 ]
; Site from funct a b
2 1 6 7 1 0.7438730952158359 -0.3213938637031660
3 1 6 7 1 -0.3213938637031660 0.7438730952158359
4 1 6 7 1 0.9516115829366237 0.2512330161413899
5 1 6 7 1 0.2512330161413898 0.9516115829366238
[ exclusions ]
1 2 3 4 5 6 7
2 1 3 4 5 6 7
3 1 2 4 5 6 7
4 1 2 3 5 6 7
5 1 2 3 4 6 7
6 1 2 3 4 5 7
7 1 2 3 4 5 6
--
David van der Spoel, Ph.D., Professor of Biology
Molec. Biophys. group, Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone: +46184714205. Fax: +4618511755.
sp...@xray.bmc.uu.se sp...@gromacs.org http://folding.bmc.uu.se
--
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/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/mailing_lists/users.php
