[gmx-users] Charge groups & Coulomb switch - artifacts in radial distribution function. Atom-based or group-based truncation?
David van der Spoel
spoel at xray.bmc.uu.se
Sun Nov 29 09:22:34 CET 2009
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.
spoel at xray.bmc.uu.se spoel at gromacs.org http://folding.bmc.uu.se
More information about the gromacs.org_gmx-users
mailing list