********>Bugfix 9:
Author: Dave Case, based on a report by Jan Fredin
Date: 08/10/06
                                                                                
Programs: sander, pbsa
                                                                                
Description: There is some IPS code in the pbsa routines that should never
             be executed.  But if it is, a divide by zero error can occur.
             The patch removes the offending code.
                                                                                
Fix:  apply the following patch to amber9/src/sander/pb_direct.f  *and* to
      amber9/src/pbsa/pb_direct.f

------------------------------------------------------------------------------
*** pb_direct.f	3 Apr 2006 23:35:55 -0000	9.0
--- pb_direct.f	11 Aug 2006 06:04:14 -0000
***************
*** 490,499 ****
      
     implicit none
  
-    ! Common variables
-    !WXW  sgld and ips head file
- #  include "sgld.h"  
-     
     ! Passed variables
      
     integer natom, iprshrt(*), iar1pb(6,0:natom)
--- 490,495 ----
***************
*** 510,521 ****
     _REAL_ dx, dy, dz, d2inv, r6
     _REAL_ df2, f2, f1, df, fw1, fw2, fw3
  
- !WXW
-    _REAL_ uips,uips2,uips6,twou,twou2,twou6,twou12,rinv
-    _REAL_ pipse,dpipse,eipse,dipse,PVC,PVCU,DVCU,PVA,PVAU,DVAU
- !WXW
- 
-     
     ! initialization
      
     eel     = ZERO; enb     = ZERO
--- 506,511 ----
***************
*** 528,536 ****
        jfirst1 = iar1pb(1, i) + 1; jlast1 = iar1pb(2, i)
        jfirst2 = iar1pb(2, i) + 1; jlast2 = iar1pb(4, i)
  
- !WXW  Get IPS NB count
-       IF(TEIPS.OR.TVIPS)NNBIPS=NNBIPS+(jlast2-jfirst1+1)*2
-        
        ! loop over direct coulombic pairs
         
        do jp = jfirst1, jlast1
--- 518,523 ----
***************
*** 539,582 ****
           d2inv = ONE/(dx**2+dy**2+dz**2)
           df2 = cn3pb(jp)*sqrt(d2inv)
           eel = eel+df2
- !WXW     IPS for electrostatic
-          IF(TEIPS)THEN
-             uips=1.0d0/rips/rinv
-             uips2=uips*uips
-             twou2=2.0d0-uips2
-             twou=sqrt(twou2)
-             pipse=bipse0+uips2*(bipse1+uips2*(bipse2+uips2*bipse3))
-             dpipse=2.0d0*bipse1+uips2*(4.0d0*bipse2+6.0d0*uips2*bipse3)
-             dipse=df2*uips*(dpipse+pipse/twou2)/twou*uips2
-             eipse=df2*uips*(pipse/twou-pipsec)
-             eel=eel+eipse
-             df2=df2-dipse
-          ENDIF
- !WXW
           r6 = d2inv**3
           f2 = cn2pb(jp)*r6
           f1 = cn1pb(jp)*(r6*r6)
           enb = enb + (f2-f1)
           df = (df2+SIX*((f2-f1)-f1))*d2inv
- !WXW     IPS for L-J
-          IF(TVIPS)THEN
-             ! L-J r6 term
-             UIPS2=1.0D0/(D2INV*RIPS2)
-             UIPS6=UIPS2*UIPS2*UIPS2
-             TWOU2=2.0D0-UIPS2
-             TWOU6=TWOU2*TWOU2*TWOU2
-             TWOU12=TWOU6*TWOU6
-             PVC=BIPSVC0+UIPS2*(BIPSVC1+UIPS2*(BIPSVC2+UIPS2*BIPSVC3))
-             PVCU=2.0D0*BIPSVC1+UIPS2*(4.0D0*BIPSVC2+6.0D0*BIPSVC3*UIPS2)
-             DVCU=(PVCU+6.0D0*PVC/TWOU2)/TWOU6
-             ! L-J r12 term
-             PVA=BIPSVA0+UIPS2*(BIPSVA1+UIPS2*(BIPSVA2+UIPS2*BIPSVA3))
-             PVAU=2.0D0*BIPSVA1+UIPS2*(4.0D0*BIPSVA2+6.0D0*BIPSVA3*UIPS2)
-             DVAU=(PVAU+12.0D0*PVA/TWOU2)/TWOU12
-             enb=enb-(f1*(PVA/TWOU12-PIPSVAC)*UIPS6-f2*(PVC/TWOU6-PIPSVCC))*UIPS6
-             df=df-(f1*DVAU*uips6-f2*DVCU)*uips6/RIPS2
-          ENDIF
- !WXW
           fw1 = dx*df; fw2 = dy*df; fw3 = dz*df
           dumx = dumx + fw1
           dumy = dumy + fw2
--- 526,536 ----
***************
*** 597,621 ****
           f1 = cn1pb(jp)*(r6*r6)
           enb = enb + (f2-f1)
           df = SIX*((f2-f1)-f1)*d2inv
- !WXW     IPS for L-J
-          IF(TVIPS)THEN
-             ! L-J r6 term
-             UIPS2=1.0D0/(D2INV*RIPS2)
-             UIPS6=UIPS2*UIPS2*UIPS2
-             TWOU2=2.0D0-UIPS2
-             TWOU6=TWOU2*TWOU2*TWOU2
-             TWOU12=TWOU6*TWOU6
-             PVC=BIPSVC0+UIPS2*(BIPSVC1+UIPS2*(BIPSVC2+UIPS2*BIPSVC3))
-             PVCU=2.0D0*BIPSVC1+UIPS2*(4.0D0*BIPSVC2+6.0D0*BIPSVC3*UIPS2)
-             DVCU=(PVCU+6.0D0*PVC/TWOU2)/TWOU6
-             ! L-J r12 term
-             PVA=BIPSVA0+UIPS2*(BIPSVA1+UIPS2*(BIPSVA2+UIPS2*BIPSVA3))
-             PVAU=2.0D0*BIPSVA1+UIPS2*(4.0D0*BIPSVA2+6.0D0*BIPSVA3*UIPS2)
-             DVAU=(PVAU+12.0D0*PVA/TWOU2)/TWOU12
-             enb=enb-(f1*(PVA/TWOU12-PIPSVAC)*UIPS6-f2*(PVC/TWOU6-PIPSVCC))*UIPS6
-             df=df-(f1*DVAU*uips6-f2*DVCU)*uips6/RIPS2
-          ENDIF
- !WXW
           fw1 = dx*df; fw2 = dy*df; fw3 = dz*df
           dumx = dumx + fw1
           dumy = dumy + fw2
--- 551,556 ----

------------------------------------------------------------------------------
                                                                                
Temporary workarounds: none