amber-developers: RE: mixed 1-4 scaling in Gromacs

From: Ross Walker <ross.rosswalker.co.uk>
Date: Tue, 9 Dec 2008 21:19:09 -0800

Hi Peter,

 

The way this is implemented right now, in AMBER (Sander and PMEMD) (the
version 11 CVS development copy) is that we have an additional definition in
the prmtop file. E.g.:

 

%FLAG SCEE_SCALE_FACTOR

%FORMAT(5E16.8)


  1.20000000E+00 1.30000000E+00 1.00000000E+00 1.70000000E+00
1.20000000E+00

  1.20000000E+00 1.30000000E+00

%FLAG SCNB_SCALE_FACTOR

%FORMAT(5E16.8)


  2.00000000E+00 2.10000000E+00 1.00000000E+00 1.00000000E+00
2.10000000E+00

  2.30000000E+00 2.30000000E+00

 

These two consist of a list of 1-4 SCEE and SCNB scale factors. The list is
ndihedral types long (since 1-4's are done over the dihedral list) -
however, the dihedral type list contains duplicates (for entries with
different PN's etc) in which case only the last value of SCEE/SCNB for that
list of dihedrals is used and the earlier ones are ignored. The same goes
for any improper dihedrals etc. If these FLAG's are not found in the prmtop
file then whatever value the user specified in the mdin file (or the default
of 1.2 / 2.0) is used for all 1-4 interactions. In the code this actually
works in the same way. There is an array (ndihedral types long) containing
the SCEE and SCNB scale factors. It is either filled with the values read
from the prmtop file or with the values of scee and scnb from the mdin file.

 

Here is some sample code from within the 1-4's that uses this:

 

  do n = 1, nb14_cnt

    i = nb14(1, n)

    j = nb14(2, n)

    parm_idx = nb14(3, n)

    scee0 = gbl_one_scee(parm_idx)

    scnb0 = gbl_one_scnb(parm_idx)

    dx = crd(1, j) - crd(1, i)

...

...

    r6 = r2inv * r2inv * r2inv

    g = charge(i) * charge(j) * rinv * scee0

    ic = ico(ntypes_lcl * (iac(i) - 1) + iac(j))

...

...

    df = (g + scnb0 * (12.d0 * f12 - 6.d0 * f6)) * r2inv

...

...

    ee14 = ee14 + g

    enb14 = enb14 + (f12 - f6) * scnb0

    frc(1, i) = frc(1, i) - df * dx

...

   end do

 

At the moment the framework for actually creating the SCEE and SCNB flags
inside the prmtop files does not exist - one has to create these arrays by
hand. My ultimate intention is that this will be created by leap. An frcmod
file would then have for example two extra sections labelled SCEE and SCNB:

 

SCEE
CT-CT-N -C 1.2
CJ-CK-CK-CG 1.0
 
SCNB
CT-CT-N -C 2.0
CJ-CK-CK-CG 2.0

 

Similar sections would exist within the parmxx.dat files. If a dihedral term
was not specified in this section then it would default to the value for the
force field (i.e. 1.2/2.0 for FF99SB etc and 1.0/1.0 for Glycam06 - with a
suitable error message or warning printed if incompatible force fields were
loaded without these additional terms). In this way backwards compatibility
with earlier force field definition files is preserved and one doesn't need
to specify the default for every force field - just for things that differ
like Glycam 06.

 

I hope this makes sense. I can send you some test cases with example outputs
if you'd like to try it. If you have any questions please don't hesitate to
contact me.

 

All the best

Ross

 

 

From: Robert J. Woods [mailto:rwoods.ccrc.uga.edu]
Sent: Tuesday, December 09, 2008 9:51 AM
To: kasson.stanford.edu; Ross Walker
Subject: Re: mixed 1-4 scaling in Gromacs

 

Hi Peter,
I'm going to ask you to direct this question to Ross Walker
(ross.rosswalker.co.uk), who implemented it in AMBER. He can no doubt give
you very precise advice.
Great chatting with you again.
Cheers!
Rob

Peter Kasson wrote:

Hi Rob,
  It was great to see you yesterday at the glycan array workshop. I'm
here with Erik Lindahl, and we're chatting about implementing mixed
1-4 scaling in gromacs. It's easy to enable on a per-molecule basis.
Enabling on a per-residue and per-atom basis is certainly possible,
but there are a bunch of ways to do it. Do you have a sense of how
you'd like to see this done? One idea would be to enable two
different 1-4 pair types, one with default 1-4 scaling and one with
user-specified scaling. And how should the interface be handled? I
assume that it goes to default at that point (we want the protein 1-4
scaling at the NLN).
 
Best wishes,
--Peter
 
 
  

 

-- 
Robert J. Woods, Ph.D.
 
Professor of Biochemistry                                 Voice: (706)
542-4454
  and Molecular Biology                                    FAX: (706)
542-4412
 
SFI Research Professor
 
University of Georgia
http://glycam.ccrc.uga.edu <http://glycam.ccrc.uga.edu/> 
Complex Carbohydrate Research Center            
315 Riverbend Road                                        "One small step
for Man,
Athens, GA 30602                                                       one
giant leap for Man-9"
Received on Wed Dec 10 2008 - 01:25:13 PST
Custom Search