[AMBER-Developers] a note regarding qm/mm shake

From: Timothy J Giese <timothyjgiese.gmail.com>
Date: Thu, 25 Feb 2016 16:06:01 -0500

I'm using gcc 5.3.1
I haven't run any of the tests.
I've been working on an out-of-main-stream branch that is based on a
commit made several months ago. So, maybe this has been fixed in the
master branch.

Having said that, if anyone else happens to notice funny behavior with
shake and qm/mm when using something like

   ntf=1,ntc=2, noshakemask=".1-72"

then this message may offer you a clue.

sander/shake.F90 loops through bonds and looks up info about the current
bond from the ib, jb, ifstwt, and noshake arrays.
When sander.F90 calls setbon (see set.F90), it may delete qm bonds
from ib and jb, but it doesn't look like those changes are propagated to
the ifstwt and noshake arrays, so the information gets offset during the
shake. (Actually, there's something more subtle here because the offset
seems to apply to the X-H bonds, in particular, and not the heavy-atom
bonds, but I haven't been able to perform a complete pathology)

The sloppy work-around that *appears* to work for me (but not well
tested) is to add the following to sander.F90, below the two calls to
setbon, which should recompute the ifstwt and noshake arrays (...and
possibly screw other stuff up instead?):

       call timer_start(TIME_FASTWT)
       call fastwat(ih(m04),nres,ix(i02),ih(m02), &
            nbonh,nbona,ix(iibh),ix(ijbh),ibelly,ix(ibellygp), &
            iwtnm,iowtnm,ihwtnm,jfastw,ix(iifstwt), &
            ix(iifstwr),ibgwat,ienwat,ibgion,ienion,iorwat, &
       call timer_stop(TIME_FASTWT)

       if( len_trim(noshakemask) > 0 ) then
          allocate( noshakegp(natom) )
          call atommask( natom, nres, 0, ih(m04), ih(m06), &
               ix(i02), ih(m02), x(lcrd), noshakemask, noshakegp )
          natnos = sum(noshakegp(1:natom))
          write(6,'(a,a,a,i5,a)') 'Noshake mask ', &
               noshakemask(1:len_trim(noshakemask)), ' matches
',natnos,' atoms'
          call setnoshake(ix,noshakegp,ntc,num_noshake)
          if( ntf > 1 ) then
             write(6,'(a)') ' Setting ntf to 1'
             ntf = 1
          end if
       end if

Near the top of subroutine sander, you'd need to add
   use findmask
    integer,allocatable :: noshakegp(:)
    integer :: natnos

(Obviously those print statements would break the tests, but whatever)

I believe the above described issue would normally be benign when all
H's in the MM region are shaked because the array offsets should simply
cause some of the water-bonds to (redundantly) shake.
Sander appears to correctly count the number of degrees of freedom in
the system with or without making the above change.


AMBER-Developers mailing list
Received on Thu Feb 25 2016 - 13:30:03 PST
Custom Search