[gmx-users] Charge groups & Coulomb switch - artifacts in radial distribution function. Atom-based or group-based truncation?

[Will Glover]

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?

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

-- 
Will Glover
UCLA Department of Chemistry & Biochemistry