Hi Tim,
I don't think noshakemask was ever tested in the approach you are trying here. That is that noshakemask overlaps the QM region. Out of interest if you were to change your input to be:
ntf=1,ntc=2,noshakemask=".25-72"
and
qmmask=".1-24"
qmshake=0
would you still see an issue or does this give you the behavior you'd expect - that is that atoms 1 to 72 don't get shaken.
I never really understood the purpose of noshakemask since you have to run with a 1fs timestep if you use it so you might as well shake nothing at all. Except if you still want to use a triangulated water model rather than a flexible water model. But even then one would sure just want shake on the waters and everything else unshaken no?
If the above suggestion works we can just add a check for overlap of noshakemask and qmmask - if it doesn't then I can go through what you suggest below in more detail.
I put this in bugzilla to track - see: http://bugzilla.ambermd.org/show_bug.cgi?id=329 (login for the page itself is amber / amberbugs) - you can create your own account once in if you don't have one.
All the best
Ross
> On Feb 25, 2016, at 1:06 PM, Timothy J Giese <timothyjgiese.gmail.com> wrote:
>
> 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"
> and
> qmmask=".1-24"
> qmshake=0
>
> 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, &
> 6,natom)
> 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,*)
> write(6,'(a,a,a,i5,a)') 'Noshake mask ', &
> noshakemask(1:len_trim(noshakemask)), ' matches
> ',natnos,' atoms'
> call setnoshake(ix,noshakegp,ntc,num_noshake)
> deallocate(noshakegp)
> 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
> and
> !!!!!!!!!!!!
> 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.
>
> -Tim
>
>
> _______________________________________________
> AMBER-Developers mailing list
> AMBER-Developers.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber-developers
_______________________________________________
AMBER-Developers mailing list
AMBER-Developers.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber-developers
Received on Thu Feb 25 2016 - 17:30:03 PST
