From 2a97e2fc004788cdf0a0a6239e438e5d660776d1 Mon Sep 17 00:00:00 2001
From: Tjerk Straatsma <straatsmatp@ornl.gov>
Date: Mon, 18 Sep 2006 15:45:00 +0000
Subject: [PATCH] New Q-HOP implementation

---
 src/cafe/GNUmakefile  |   3 +-
 src/cafe/cf_common.fh |   9 +-
 src/cafe/cf_f.fh      | 284 ++++++++++++++----------------------------
 src/cafe/cf_finish.F  |   7 +-
 src/cafe/cf_forces.fh |   4 +-
 src/cafe/cf_hopping.F | 134 +++++++++++++-------
 src/cafe/cf_init.F    |  30 ++++-
 src/cafe/cf_list.F    |  65 +++++++---
 src/cafe/cf_pairs.F   |   5 +-
 src/cafe/cf_rdtop.F   |  15 +--
 src/cafe/cf_shake.F   |  16 ++-
 src/nwmd/md_main.F    |  20 ++-
 src/qhop/qhop_asite.F | 250 ++++++++++++++++++++++++++++++++-----
 src/qhop/qhop_prob.F  |  19 +--
 src/space/GNUmakefile |  13 +-
 src/space/sp_main.F   |  61 ++++++++-
 16 files changed, 596 insertions(+), 339 deletions(-)

diff --git a/src/cafe/GNUmakefile b/src/cafe/GNUmakefile
index b44a073918..6210ed5ac9 100644
--- a/src/cafe/GNUmakefile
+++ b/src/cafe/GNUmakefile
@@ -1,5 +1,5 @@
 #
-# $Id: GNUmakefile,v 1.8 2003-07-17 19:44:23 d3j191 Exp $
+# $Id: GNUmakefile,v 1.9 2006-09-18 15:42:14 d3j191 Exp $
 #
 OBJ_OPTIMIZE=cf_d.o cf_fp.o cf_f.o\
              cf_finish.o\
@@ -15,6 +15,7 @@ OBJ_OPTIMIZE=cf_d.o cf_fp.o cf_f.o\
              cf_shake.o\
              cf_grid.o\
              cf_rdtop.o\
+             cf_uhop.o\
              pme_init.o\
              pme_fft.o\
              pme_grid.o
diff --git a/src/cafe/cf_common.fh b/src/cafe/cf_common.fh
index 4e1bfa9ed5..a45c7cdf5e 100644
--- a/src/cafe/cf_common.fh
+++ b/src/cafe/cf_common.fh
@@ -1,4 +1,4 @@
-c $Id: cf_common.fh,v 1.52 2005-09-14 22:15:56 d3j191 Exp $
+c $Id: cf_common.fh,v 1.53 2006-09-18 15:42:14 d3j191 Exp $
 c
       real*8 zero,tiny,small
       real*8 fourth,third,half,one,two,three,four,six,twelve
@@ -18,6 +18,9 @@ c
       real*8 avoga,rgas
       parameter(avoga=6.022045d-07)
       parameter(rgas=8.31441d-03)
+c
+      integer mis2
+      parameter(mis2=15)
 c
       integer mbs,mpe
       parameter(mbs=2)
@@ -72,7 +75,7 @@ c
      + i_s3r1,l_s3r1,i_s3r2,l_s3r2,i_smr1,l_smr1,
      + i_smr2,l_smr2,i_smr3,l_smr3,i_smr4,l_smr4,
      + i_smr5,l_smr5,i_smr6,l_smr6,
-     + i_list,l_list,i_snam,l_snam,i_sti,l_sti,
+     + i_list,l_list,i_snam,l_snam,i_sti,l_sti,i_mprot,l_mprot,
      + i_pmew,l_pmew,i_pmes,l_pmes,
      + i_theta,l_theta,i_dtheta,l_dtheta,
      + i_ipmf,l_ipmf,i_jpmf,l_jpmf,i_rpmf,l_rpmf,i_xpmf,l_xpmf,
@@ -96,7 +99,7 @@ c
      + i_s3r1,l_s3r1,i_s3r2,l_s3r2,i_smr1,l_smr1,
      + i_smr2,l_smr2,i_smr3,l_smr3,i_smr4,l_smr4,
      + i_smr5,l_smr5,i_smr6,l_smr6,
-     + i_list,l_list,i_snam,l_snam,i_sti,l_sti,
+     + i_list,l_list,i_snam,l_snam,i_sti,l_sti,i_mprot,l_mprot,
      + i_pmew,l_pmew,i_pmes,l_pmes,
      + i_theta,l_theta,i_dtheta,l_dtheta,
      + i_ipmf,l_ipmf,i_jpmf,l_jpmf,i_rpmf,l_rpmf,i_xpmf,l_xpmf,
diff --git a/src/cafe/cf_f.fh b/src/cafe/cf_f.fh
index bafe055e8f..94c8833da4 100644
--- a/src/cafe/cf_f.fh
+++ b/src/cafe/cf_f.fh
@@ -13,7 +13,7 @@ c,rdf,ngt,iagc,jagc,igrc,u,uwmw)
 c error
 #endif
 c
-c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
+c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
 c
       implicit none
 c
@@ -759,7 +759,7 @@ c
 c error
 #endif
 c
-c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
+c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
 c
       implicit none
 c
@@ -786,7 +786,7 @@ c
       integer isal(mscr),isrx(mscr)
 c
       integer lswj(*)
-      integer nums
+      integer nums,i
       integer lswndx(0:msa,2),lswjpt(nums,2),lswin(nums,2)
       integer list(0:msa)
 c
@@ -821,6 +821,7 @@ c
       integer is,inum,jnum,jwm,iptr,jwa,iss,isg
       logical lid,ljd
 #endif
+      real*8 rtmp
 c
 #include "cf_funcs_sfn.fh"
 #include "bitops_funcs.fh"
@@ -833,6 +834,7 @@ c
       boxi(1)=one/box(1)
       boxi(2)=one/box(2)
       boxi(3)=one/box(3)
+      print*,'icntrl=',icntrl
       if(npbtyp.eq.1.and.nbxtyp.eq.0.and.icntrl.eq.2) then
       do 101 ipsw=1,npsw
       do 102 is=1,nums
@@ -843,6 +845,7 @@ c
       inum=lswjpt(is,ipsw)-1
       if(iand(isss(isa),6).eq.2) iss=-1
       if(iand(isss(isa),6).eq.4) iss=-1
+c      write(*,'(a,4i5)') 'iss ',isa,isss(isa),iand(isss(isa),6),iss
       isg=isga(isa)
       lid=iand(isdt(isa),mdynam).eq.ldynam
       qi=chg(isq1(isa),1,iset)
@@ -1174,6 +1177,9 @@ c      isrx(nax+ismn)=isrst
 c   10 continue
 c      endif
 c
+c      write(*,'(4i5,2f12.6)') 
+c     + lssscl,isga(isa),isss(isfr+isa),iand(isss(isfr+isa),6),
+c     + shift0(1),shift1(1)
       if(lssscl) then
       do 10 ismn=1,lswin(isa,ipsw)
 c      isrx(nax+ismn)=isss(isfr+isa)
@@ -1285,8 +1291,11 @@ c      endif
       rwx(iax,2)=xi(iax,2)-xj(iax,2,iwa)
       rwx(iax,3)=xi(iax,3)-xj(iax,3,iwa)
       rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
+      rtmp=rwi2(iax)
       if(isrx(iax).gt.0) rwi2(iax)=one/(one/rwi2(iax)+shift0(1))
       if(isrx(iax).lt.0) rwi2(iax)=one/(one/rwi2(iax)+shift1(1))
+c      write(*,'(3i5,2f12.6)')
+c     + isga(isal(iax)),isal(iax),isrx(iax),rtmp,rwi2(iax)
    27 continue
 c
 c     Lennard-Jones interactions
@@ -1889,13 +1898,14 @@ c
 c error
 #endif
 c
-c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
+c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
 c
       implicit none
 c
 #include "cf_common.fh"
 #include "mafdecls.fh"
 c
+      real*8 rtmp
       real*8 xs(msa,3),xsm(msm,3),fs(msa,3,2)
       real*8 zs(msf,3,3,2),ess(msf,msf,mpe,2)
       real*8 fss(msf,msf,3,2)
@@ -1920,8 +1930,8 @@ c
 c
 #if defined(CAFE_FORCES)
       real*8 dera(6,nsatot)
-      integer lda(10,*)
-      real*8 rda(5,*),uda(4,*)
+      integer lda(16,*)
+      real*8 rda(11,*),uda(4,*)
 #endif
 #if defined(CAFE_POLARIZATION)
       real*8 ps(msa,3,2),psp(msa,3,2,2)
@@ -1933,11 +1943,7 @@ c
       integer ipss,number,isslen,nax,jsaptr
       integer jnum,lssptr,iax
       real*8 dercon
-#if !defined(CAFE_POLARIZATION)
-      integer iss,jss,ilda,ndxqi,ndxqj
-      integer idon,iacc,idong,iaccg,idons,iaccs,idonp,iaccp
-      real*8 uhop(4),xacc(3),racci1
-#endif
+c
       real*8 c6,c12,cf6,cf12
       real*8 c64,c124
       real*8 q14,sumen1,sumen2,sumen3
@@ -2095,6 +2101,9 @@ c
       jstt=iand(isss(lssptr),384)
       if(istt.eq.128.and.jstt.eq.256) isrx(nax+jnum)=999
       if(istt.eq.256.and.jstt.eq.128) isrx(nax+jnum)=999
+c
+c      write(*,'(5i5)')
+c     + isga(isfr+isa),isga(lssptr),istt,jstt,isrx(nax+jnum)
 c
       endif
 c
@@ -2236,6 +2245,7 @@ c
       rwi2(iax)=zero
       isrx(iax)=0
       endif
+      rtmp=rwi2(iax)
       if(isrx(iax).gt.0) then
       rwi2(iax)=one/(one/rwi2(iax)+shift0(1))
       elseif(isrx(iax).lt.0) then
@@ -2244,6 +2254,8 @@ c
       rwi1(iax)=sqrt(rwi2(iax))
       isa=isal(iax)
       jsa=jsal(iax)
+c      write(*,'(3i5,4f12.6)')
+c     + isga(isa),isga(jsa),isrx(iax),shift0(1),shift1(1),rtmp,rwi2(iax)
       if(jfal(iax).ne.isf) then
       q=qsa2(isa)*qsa2(jsa)
       else
@@ -2574,177 +2586,6 @@ cx      endif
       nax=nax+lssin(isa,ipss)
    36 continue
 c
-#if defined(CAFE_FORCES)
-c
-c     qhopping
-c
-c     lda(1,i)  : donor global atom number
-c     lda(2,i)  : donor local atom number
-c     lda(3,i)  : donor segment number
-c     lda(4,i)  : donor protonation state after hopping
-c     lda(5,i)  : donor charge mutation
-c     lda(6,i)  : acceptor global atom number
-c     lda(7,i)  : acceptor local atom number
-c     lda(8,i)  : acceptor segment number
-c     lda(9,i)  : acceptor protonation state after hopping
-c     lda(10,i) : acceptor charge mutation
-c
-      if(lhop) then
-      do 2001 ilda=1,nldat
-c
-      uhop(1)=0.0d0
-      uhop(2)=0.0d0
-      uhop(3)=0.0d0
-      uhop(4)=0.0d0
-c
-      idon=lda(2,ilda)
-      idong=lda(1,ilda)
-      idons=lda(3,ilda)
-      idonp=lda(4,ilda)
-      iacc=lda(7,ilda)
-      iaccg=lda(6,ilda)
-      iaccs=lda(8,ilda)
-      iaccp=lda(9,ilda)
-c
-c      write(*,'(a,8i5)') 'DA ',idon,idong,idons,idonp,
-c     + iacc,iaccg,iaccs,iaccp
-c
-c     Energy of the donor and acceptor segment before hopping
-c
-      do 2002 iax=1,nax
-      isa=isal(iax)
-      jsa=jsal(iax)
-      iss=isgm(isa)
-      jss=isgm(jsa)
-c      write(*,'(4i4)') isa,jsa,iss,jss
-c
-c     skip all intrasegment interactions
-c
-      if(iss.ne.jss.and.(iss.eq.idons.or.jss.eq.idons.or.
-     + iss.eq.iaccs.or.jss.eq.iaccs)) then
-c
-      rwi2(iax)=one/(rwx(iax,1)**2+rwx(iax,2)**2+rwx(iax,3)**2)
-      rwi1(iax)=sqrt(rwi2(iax))
-c      write(*,'(4i4,f12.6)') isa,jsa,iss,jss,rwi1(iax)
-c
-c     donor energy
-c
-      if((iss.eq.idons.and.jss.ne.iaccs).or.
-     +   (jss.eq.idons.and.iss.ne.iaccs)) then
-      if(jfal(iax).ne.isf) then
-      q=chg(isq2(isa),1,lseq(iss))*
-     + chg(isq2(jsa),1,lseq(jss))
-      else
-      q=chg(isq3(isa),1,lseq(iss))*
-     + chg(isq3(jsa),1,lseq(jss))
-      endif
-      uhop(1)=uhop(1)+facu(iax)*q*rwi1(iax)
-c      write(*,'(8i5,2f12.3)') 
-c     + isa,lseq(isgm(isa)),isq2(isa),isq3(isa),
-c     + jsa,lseq(isgm(jsa)),isq2(jsa),isq3(jsa),q,uhop(1)
-      endif
-c
-c     acceptor energy
-c
-      if((iss.eq.iaccs.and.jss.ne.idons).or.
-     +   (jss.eq.iaccs.and.iss.ne.idons)) then
-      if(jfal(iax).ne.isf) then
-      q=chg(isq2(isa),1,lseq(isgm(isa)))*
-     + chg(isq2(jsa),1,lseq(isgm(jsa)))
-      else
-      q=chg(isq3(isa),1,lseq(isgm(isa)))*
-     + chg(isq3(jsa),1,lseq(isgm(jsa)))
-      endif
-      uhop(2)=uhop(2)+facu(iax)*q*rwi1(iax)
-      endif
-c
-c     donor energy after hopping
-c
-      if(iss.eq.idons.and.jss.ne.iaccs) then
-      ndxqi=isa
-      if(isga(isa).eq.idong) ndxqi=isa+lda(5,ilda)
-      if(isga(isa).eq.idong+lda(5,ilda)) ndxqi=isa-lda(5,ilda)
-      if(jfal(iax).ne.isf) then
-      q=chg(isq2(ndxqi),1,idonp)*chg(isq2(jsa),1,lseq(isgm(jsa)))
-      else
-      q=chg(isq3(ndxqi),1,idonp)*chg(isq3(jsa),1,lseq(isgm(jsa)))
-      endif
-      uhop(3)=uhop(3)+facu(iax)*q*rwi1(iax)
-      endif
-      if(jss.eq.idons.and.iss.ne.iaccs) then
-      ndxqj=jsa
-      if(isga(jsa).eq.idong) ndxqj=jsa+lda(5,ilda)
-      if(isga(jsa).eq.idong+lda(5,ilda)) ndxqj=jsa-lda(5,ilda)
-      if(jfal(iax).ne.isf) then
-      q=chg(isq2(isa),1,lseq(isgm(isa)))*chg(isq2(ndxqj),1,idonp)
-      else
-      q=chg(isq3(isa),1,lseq(isgm(isa)))*chg(isq3(ndxqj),1,idonp)
-      endif
-      uhop(3)=uhop(3)+facu(iax)*q*rwi1(iax)
-      endif
-c
-c     acceptor energy after hopping
-c
-      if(iss.eq.iaccs.and.jss.ne.idons) then
-      if(isga(isa).eq.iaccg) then
-      xacc(1)=rda(1,ilda)-xj(iax,1)
-      xacc(2)=rda(2,ilda)-xj(iax,2)
-      xacc(3)=rda(3,ilda)-xj(iax,3)
-      call cf_pbc(1,xacc,1,xacc,1,0,1,1)
-      else
-      xacc(1)=xi(iax,1)-xj(iax,1)
-      xacc(2)=xi(iax,2)-xj(iax,2)
-      xacc(3)=xi(iax,3)-xj(iax,3)
-      endif
-      racci1=one/sqrt(xacc(1)**2+xacc(2)**2+xacc(3)**2)
-      ndxqi=isa
-      if(isga(isa).eq.iaccg) ndxqi=isa+lda(10,ilda)
-      if(isga(isa).eq.iaccg+lda(10,ilda)) ndxqi=isa-lda(10,ilda)
-      if(jfal(iax).ne.isf) then
-      q=chg(isq2(ndxqi),1,iaccp)*chg(isq2(jsa),1,lseq(isgm(jsa)))
-      else
-      q=chg(isq3(ndxqi),1,iaccp)*chg(isq3(jsa),1,lseq(isgm(jsa)))
-      endif
-      uhop(4)=uhop(4)+facu(iax)*q*racci1
-      endif
-      if(jss.eq.iaccs.and.iss.ne.idons) then
-      if(isga(jsa).eq.iaccg) then
-      xacc(1)=xi(iax,1)-rda(1,ilda)
-      xacc(2)=xi(iax,2)-rda(2,ilda)
-      xacc(3)=xi(iax,3)-rda(3,ilda)
-      call cf_pbc(1,xacc,1,xacc,1,0,1,1)
-      else
-      xacc(1)=xi(iax,1)-xj(iax,1)
-      xacc(2)=xi(iax,2)-xj(iax,2)
-      xacc(3)=xi(iax,3)-xj(iax,3)
-      endif
-      racci1=one/sqrt(xacc(1)**2+xacc(2)**2+xacc(3)**2)
-      ndxqj=jsa
-      if(isga(jsa).eq.iaccg) ndxqj=jsa+lda(5,ilda)
-      if(isga(jsa).eq.iaccg+lda(10,ilda)) ndxqj=jsa-lda(10,ilda)
-      if(jfal(iax).ne.isf) then
-      q=chg(isq2(isa),1,lseq(isgm(isa)))*chg(isq2(ndxqj),1,iaccp)
-      else
-      q=chg(isq3(isa),1,lseq(isgm(isa)))*chg(isq3(ndxqj),1,iaccp)
-      endif
-      uhop(4)=uhop(4)+facu(iax)*q*racci1
-      endif
-c
-      endif
- 2002 continue
-c
-      uda(1,ilda)=uda(1,ilda)+uhop(1)
-      uda(2,ilda)=uda(2,ilda)+uhop(2)
-      uda(3,ilda)=uda(3,ilda)+uhop(3)
-      uda(4,ilda)=uda(4,ilda)+uhop(4)
-c
-c      write(*,'(a,8f12.3)') 'uhop ',(uda(iax,ilda),iax=1,4),uhop
-c
- 2001 continue
-      endif
-c
-#endif
-c
 c     thermodynamic integration
 c
       if(ithint) then
@@ -2783,6 +2624,9 @@ c
      + rwi2(nax+jnum)*rwi6(nax+jnum)*(c64-c124*rwi6(nax+jnum))
       endif
       deriv(15,ipss)=deriv(15,ipss)+dercon
+c      write(*,'(a,3i5,4f12.6)') 'gv ',
+c     + isga(isfr+isa),isga(jsa),isrx(nax+jnum),shift0(4),shift1(4),
+c     + dercon,deriv(15,ipss)
 #if defined(CAFE_FORCES)
       if(npgdec.gt.1) then
       dera(3,isga(isfr+isa))=dera(3,isga(isfr+isa))+half*dercon
@@ -2816,6 +2660,9 @@ c
       if(ireact.ne.0) then
       deriv(17,ipss)=deriv(17,ipss)+dercon*rffss/rwi2(nax+jnum)
       endif
+c      write(*,'(a,3i5,4f12.6)') 'gq ',
+c     + isga(isfr+isa),isga(jsa),isrx(nax+jnum),shift0(4),shift1(4),
+c     + dercon,deriv(17,ipss)
 #if defined(CAFE_FORCES)
       if(npgdec.gt.1) then
       dera(4,isga(isfr+isa))=dera(4,isga(isfr+isa))+
@@ -3222,7 +3069,7 @@ c
      + natoms,ndim,igan,isgm,imol,idyn,ichg,chg,xs,fs,ess,lpbc,lpbcs,
      + lupden,lupdti,dera,lseq)
 c
-c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
+c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
 c
 c     cf_fsb returns forces and energies for solute bonds
 c
@@ -3486,7 +3333,7 @@ c
       endif
 c
       if(ipme.ne.0) then
-      qij=chg(ichg(isa),1,1)*chg(ichg(jsa),1,1)
+      qij=chg(ichg(isa),1,ibset)*chg(ichg(jsa),1,ibset)
       ferfc=one-erfc(ealpha*rss)
       fderfc=-(ealpha*derfc(ealpha*rss))
       epmecs=epmecs-ferfc*qij*rssi
@@ -3532,7 +3379,7 @@ c
      + natoms,ndim,igan,isgm,imol,idyn,ichg,chg,xs,fs,ess,lpbc,lpbcs,
      + lupden,lupdti,dera,lseq)
 c
-c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
+c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
 c
 c     cf_fsh returns forces and energies for solute angles
 c
@@ -3955,7 +3802,7 @@ c
      + natoms,ndim,igan,isgm,imol,idyn,xs,fs,ess,lpbc,lpbcs,lupden,
      + lupdti,dera,lseq)
 c
-c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
+c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
 c
 c     cf_fsd returns forces and energies for solute angles
 c
@@ -4301,7 +4148,7 @@ c
      + natoms,ndim,igan,isgm,imol,idyn,xs,fs,ess,lpbc,lpbcs,lupden,
      + lupdti,dera,lseq)
 c
-c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
+c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
 c
 c     cf_fso returns forces and energies for solute improper dihedral angles
 c
@@ -4628,9 +4475,10 @@ c
       return
       end
       subroutine cf_fst(nthrds,indexl,mst,idtn,vdw,chg,natoms,ndim,
-     + iatt,igan,isgm,imol,idyn,ichg,xs,fs,ess,lpbc,lpbcs,dera,esa,lseq)
+     + iatt,igan,isgm,imol,idyn,ichg,isss,xs,fs,ess,lpbc,lpbcs,dera,
+     + esa,lseq)
 c
-c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
+c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
 c
 c     cf_fst returns forces and energies for solute third neighbors
 c
@@ -4682,6 +4530,7 @@ c
       integer natoms,ndim
       real*8 vdw(mat,mat,map,mset),chg(mqt,mqp,mset)
       integer iatt(ndim),igan(ndim),imol(ndim),idyn(ndim),ichg(ndim)
+      integer isss(msa)
       real*8 xs(ndim,3),fs(ndim,3)
       real*8 ess(msf,msf,mpe,2),esa(nsa)
 c
@@ -4694,7 +4543,7 @@ c
 c     declaration of local variables
 c     ------------------------------
 c
-      integer i,j,ist,isa,jsa,isf,jsf,isad
+      integer i,j,ist,isa,jsa,isf,jsf,isad,istt,jstt,isrx
       real*8 factu,etermq,eterml
       real*8 rxx,rxy,rxz,r2,r2i,r1i,dfs,dfsu,r6i,c6,c12,cf6,cf12,q
       real*8 ferfc,fderfc,qij,dfsp,c64,c124,dercon,qi,qi4,qj,qj4,q14
@@ -4746,6 +4595,30 @@ c     find solute molecule numbers involved
 c
       isf=imol(isa)
       jsf=imol(jsa)
+c
+      if(lssscl) then
+c
+      istt=iand(isss(isa),48)
+      jstt=iand(isss(jsa),48)
+      if(isf.ne.jsf) then
+      if(istt.eq.16.or.jstt.eq.16) isrx=-1
+      if(istt.eq.32.or.jstt.eq.32) isrx=1
+      endif
+c
+      istt=iand(isss(isa),384)
+      jstt=iand(isss(jsa),384)
+      if(istt.eq.128.or.jstt.eq.128) isrx=-2
+      if(istt.eq.256.or.jstt.eq.256) isrx=2
+c
+      istt=iand(isss(isa),384)
+      jstt=iand(isss(jsa),384)
+      if(istt.eq.128.and.jstt.eq.256) isrx=999
+      if(istt.eq.256.and.jstt.eq.128) isrx=999
+c
+c      write(*,'(5i5)')
+c     + isga(isfr+isa),isga(lssptr),istt,jstt,isrx(nax+jnum)
+c
+      endif
 c
 c     determine energy factor depending on atoms being dynamic
 c     or fixed
@@ -4779,6 +4652,19 @@ c
 c
       r2=rxx*rxx+rxy*rxy+rxz*rxz
       r2i=one/r2
+c
+      if(lssscl) then
+      if(isrx.eq.999) then
+      r2i=zero
+      isrx=0
+      endif
+      if(isrx.gt.0) then
+      r2i=one/(one/r2i+shift0(1))
+      elseif(isrx.lt.0) then
+      r2i=one/(one/r2i+shift1(1))
+      endif
+      endif
+c
       r1i=sqrt(r2i)
       dfs=zero
       dfsu=zero
@@ -4894,6 +4780,15 @@ c      c124=cb12(iatt(isa),iatt(jsa),4)
       c64=vdw(iatt(isa),iatt(jsa),2,4)
       c124=vdw(iatt(isa),iatt(jsa),4,4)
       dercon=(c124*r6i-c64)*r6i
+      if(isrx.gt.0) then
+      c64=three*vdw(iatt(isa),iatt(jsa),1,iset)
+      c124=six*vdw(iatt(isa),iatt(jsa),3,iset)
+      dercon=dercon+shift0(4)*r2i*r6i*(c64-c124*r6i)
+      elseif(isrx.lt.0) then
+      c64=three*vdw(iatt(isa),iatt(jsa),1,iset)
+      c124=six*vdw(iatt(isa),iatt(jsa),3,iset)
+      dercon=dercon+shift1(4)*r2i*r6i*(c64-c124*r6i)
+      endif
       deriv(15,1)=deriv(15,1)+dercon
 #if defined(CAFE_FORCES)
       if(npgdec.gt.1) then
@@ -4915,6 +4810,13 @@ c      qj4=qs(jsa,4)
       dercon=q14fac*(qi*qj4+qj*qi4)*r1i
       else
       dercon=q14fac*(qi*qj4+qj*qi4)*r1i
+      endif
+      if(isrx.gt.1) then
+      dercon=dercon-half*shift0(4)*
+     + chg(ichg(isa),1,iset)*chg(ichg(jsa),1,iset)*r2i
+      elseif(isrx.lt.-1) then
+      dercon=dercon-half*shift1(4)*
+     + chg(ichg(isa),1,iset)*chg(ichg(jsa),1,iset)*r2i
       endif
       deriv(17,1)=deriv(17,1)+dercon
       if(ireact.ne.0) then
@@ -4993,7 +4895,7 @@ c
      + mwd,nwd,ndp,idih,dih,rdih,mwo,nwo,nop,iimp,dimp,rimp,
      + mwt,nwt,idwt,mwn,nwn,idwn)
 c
-c $Id: cf_f.fh,v 1.63 2006-02-07 20:17:04 d3j191 Exp $
+c $Id: cf_f.fh,v 1.64 2006-09-18 15:42:14 d3j191 Exp $
 c
       implicit none
 c
diff --git a/src/cafe/cf_finish.F b/src/cafe/cf_finish.F
index be42713998..3be7174163 100644
--- a/src/cafe/cf_finish.F
+++ b/src/cafe/cf_finish.F
@@ -1,5 +1,5 @@
       subroutine cf_finish()
-c $Id: cf_finish.F,v 1.22 2003-11-25 16:40:12 d3j191 Exp $
+c $Id: cf_finish.F,v 1.23 2006-09-18 15:42:14 d3j191 Exp $
       implicit none
 c
 #include "cf_common.fh"
@@ -142,6 +142,11 @@ c
       if(.not.ma_pop_stack(l_sti))
      + call md_abort('Failed to deallocate sti',0)
 c
+      if(lqhop) then
+      if(.not.ma_pop_stack(l_mprot))
+     + call md_abort('Failed to deallocate mprot',0)
+      endif
+c 
       if(.not.ma_pop_stack(l_snam))
      + call md_abort('Failed to deallocate snam',0)
 c
diff --git a/src/cafe/cf_forces.fh b/src/cafe/cf_forces.fh
index 2009605f30..cd042d5160 100644
--- a/src/cafe/cf_forces.fh
+++ b/src/cafe/cf_forces.fh
@@ -1,4 +1,4 @@
-c $Id: cf_forces.fh,v 1.23 2005-01-08 19:17:20 d3j191 Exp $
+c $Id: cf_forces.fh,v 1.24 2006-09-18 15:42:14 d3j191 Exp $
 #if defined(CAFE_POLARIZATION)
       subroutine forcep(lself,lpbcs,
      + xw,xwm,fw,pw,pwp,zwi,rtos,iwdt,iwz,
@@ -349,7 +349,7 @@ c
       if(npairs.gt.0) then
       call cf_fst(npairs,int_mb(lptr+1),mtt(2),int_mb(i_itrd(2)),
      + dbl_mb(i_vdw),dbl_mb(i_chg),
-     + max(isto,jsto),msa,isat,isga,isgm,ismf,isdt,isq3,
+     + max(isto,jsto),msa,isat,isga,isgm,ismf,isdt,isq3,isss,
      + xs,fs,ess,lpbc,lpbcs,dbl_mb(i_sti),esa,int_mb(i_lseq))
       lptr=lptr+1+npairs
       else
diff --git a/src/cafe/cf_hopping.F b/src/cafe/cf_hopping.F
index 06b1d94b72..3f43dcae83 100644
--- a/src/cafe/cf_hopping.F
+++ b/src/cafe/cf_hopping.F
@@ -1,7 +1,7 @@
-      logical function cf_hopping(stimei,issgm,isgan,ishop,
-     + xs,nsaloc)
+      logical function cf_hopping(lpbc,lpbcs,stimei,
+     + isl,issgm,isgan,isq3,ishop,xs,nsaloc)
 c
-c $Id: cf_hopping.F,v 1.16 2005-10-03 21:03:12 edo Exp $
+c $Id: cf_hopping.F,v 1.17 2006-09-18 15:42:14 d3j191 Exp $
 c
       implicit none
 c
@@ -13,13 +13,16 @@ c
       logical cf_hop
       external cf_hop
 c
-      integer issgm(msa),nsaloc,isgan(msa),ishop(msa)
+      logical lpbc,lpbcs
+      integer isl(msa,mis2)
+      integer issgm(msa),nsaloc,isgan(msa),ishop(msa),isq3(msa)
       real*8 xs(msa,3)
       real*8 stimei
 c
       integer i_itmp,l_itmp,i_dtmp,l_dtmp
 c
       stime=stimei
+      cf_hopping=.false.
 c
       if(lpair.and.nhop.gt.0) then
       lhop=.true.
@@ -27,13 +30,20 @@ c
       if(.not.ma_push_get(mt_int,np,'itmp',l_itmp,i_itmp))
      + call md_abort('Failed to allocate itmp',me)
 c
-      call cf_hoplist(int_mb(i_itmp),int_mb(i_lda),dbl_mb(i_rda))
+      call cf_hoplist(issgm,nsaloc,int_mb(i_itmp),int_mb(i_lda),
+     + dbl_mb(i_rda))
 c
       if(.not.ma_pop_stack(l_itmp))
      + call md_abort('Failed to deallocate itmp',me)
 c
       cf_hopping=.true.
       elseif(lhop) then
+c
+c     evaluate the E12
+c
+      call cf_uhop(lpbc,lpbcs,int_mb(i_lseq),isl,issgm,isgan,isq3,
+     + xs,nsaloc,
+     + dbl_mb(i_chg),int_mb(i_lda),dbl_mb(i_rda),dbl_mb(i_uda))
 c
       if(.not.ma_push_get(mt_dbl,nldat,'dtmp',l_dtmp,i_dtmp))
      + call md_abort('Failed to allocate dtmp',me)
@@ -54,7 +64,7 @@ c
 c
       return
       end
-      subroutine cf_hoplist(nhopl,lda,rda)
+      subroutine cf_hoplist(issgm,nsaloc,nhopl,lda,rda)
 c
       implicit none
 c
@@ -62,10 +72,11 @@ c
 #include "global.fh"
 #include "mafdecls.fh"
 #include "msgids.fh"
+#include "util.fh"
 c
-      integer nhopl(np)
-      integer lda(10,*)
-      real*8 rda(5,*)
+      integer issgm(msa),nsaloc,nhopl(np)
+      integer lda(16,*)
+      real*8 rda(11,*)
 c
       integer i,j,ioff
 c
@@ -73,52 +84,82 @@ c
       nhopl(i)=0
     1 continue
       nhopl(me+1)=nlda
+c      write(*,'(i5,a,i5)') me,' nlda= ',nlda
 c
       call ga_igop(mcf_76,nhopl,np,'+')
 c
       do 2 i=2,np
       nhopl(i)=nhopl(i)+nhopl(i-1)
     2 continue
-c      if(me.eq.0) write(*,'(a,10i5)') 'nhopl ',(nhopl(i),i=1,np)
+c      write(*,'(i5,a,10i5)') me,' nhopl ',(nhopl(i),i=1,np)
 c
-c      do 222 i=1,nlda
-c      write(*,'(12i5)') me,i,(lda(j,i),j=1,10)
-c  222 continue
 c
       ioff=0
       if(me.gt.0) then
       ioff=nhopl(me)
       do 3 i=1,nlda
-      do 4 j=1,10
+      do 4 j=1,16
       lda(j,i+ioff)=lda(j,i)
     4 continue
-      do 5 j=1,5
+      do 5 j=1,8
       rda(j,i+ioff)=rda(j,i)
     5 continue
     3 continue
       do 6 i=1,nhopl(me)
-      do 7 j=1,10
+      do 7 j=1,16
       lda(j,i)=0
     7 continue
-      do 8 j=1,5
+      do 8 j=1,8
       rda(j,i)=0.0d0
     8 continue
     6 continue
       endif
       if(me.lt.np-1) then
       do 9 i=nhopl(me+1)+1,nhopl(np)
-      do 10 j=1,10
+      do 10 j=1,16
       lda(j,i)=0
    10 continue
-      do 11 j=1,5
+      do 11 j=1,8
       rda(j,i)=0.0d0
    11 continue
     9 continue
       endif
 c
       nldat=nhopl(np)
-      call ga_igop(mcf_77,lda,10*nldat,'+')
-      call ga_dgop(mcf_78,rda,5*nldat,'+')
+      call ga_igop(mcf_77,lda,16*nldat,'+')
+c
+c      do 322 i=1,nldat
+c      write(*,'(18i5)') me,i,(lda(j,i),j=1,16)
+c  322 continue
+c
+      do 201 i=1,nldat
+      do 202 j=1,nsaloc
+      if(lda(3,i).eq.issgm(j)) then
+      if(lda(12,i).eq.0) lda(12,i)=j
+      lda(13,i)=j
+      endif
+      if(lda(8,i).eq.issgm(j)) then
+      if(lda(15,i).eq.0) lda(15,i)=j
+      lda(16,i)=j
+      endif
+  202 continue
+  201 continue
+c
+      if(me.gt.0) then
+      do 221 j=1,nldat
+      do 222 i=1,11
+      lda(i,j)=0
+  222 continue
+      lda(14,j)=0
+  221 continue
+      endif
+c
+      call ga_igop(mcf_77,lda,16*nldat,'+')
+      call ga_dgop(mcf_78,rda,11*nldat,'+')
+c
+c      do 2222 i=1,nlda
+c      write(*,'(18i5)') me,i,(lda(j,i),j=1,16)
+c 2222 continue
 c
       nlda=0
       do 12 i=1,nldat
@@ -126,24 +167,31 @@ c
       if(lda(1,i).eq.lda(1,j).and.lda(6,i).eq.lda(6,j)) goto 12
    13 continue
       nlda=nlda+1
-      do 14 j=1,10
+      do 14 j=1,16
       lda(j,nlda)=lda(j,i)
    14 continue
-      do 15 j=1,5
+      do 15 j=1,8
       rda(j,nlda)=rda(j,i)
    15 continue
    12 continue
       nldat=nlda
 c
-ctest
-c      if(me.eq.0) then
-c      print*,'Number of donor-acceptor pairs is ',nldat
-c      do 12 i=1,nldat
-c      write(*,111) (lda(j,i),j=1,10),(rda(j,i),j=1,3)
-c  111 format(10x,10i5,3f12.6)
-c   12 continue
-c      endif
-ctest
+c
+      if(me.eq.0) then
+      if(util_print('qhop',print_high)) then
+      write(lfnhop,110) nldat
+  110 format(/,'Number of donor-acceptor pairs is ',i5,/,
+     + '---------- Donor --------- --------- Acceptor ----',
+     + ' --- Donor ---  ---Acceptor---',/,
+     + ' glob  loc  sgm prot  off glob  loc  sgm prot  off',
+     + '  own from   to  own from   to',/)
+      do 212 i=1,nldat
+      write(lfnhop,111) (lda(j,i),j=1,16),(rda(j,i),j=1,3)
+  111 format(16i5,3f12.6)
+  212 continue
+      endif
+      endif
+c
       return
       end
       logical function cf_hop(lseq,issgm,lda,rda,uda,ndx,pda,
@@ -160,40 +208,40 @@ c
 #include "bitops_funcs.fh"
 c
       integer lseq(mseq),issgm(msa)
-      integer lda(10,*),nsaloc,isgan(msa),ishop(msa)
-      real*8 rda(5,*),uda(4,*),pda(*),xs(msa,3),x(3)
+      integer lda(16,*),nsaloc,isgan(msa),ishop(msa)
+      real*8 rda(11,*),uda(4,*),pda(*),xs(msa,3),x(3)
       integer ndx(*)
       integer lsthop(2,*)
       real*8 timhop(*)
 c
       integer i,j,k,ndxrel
-      real*8 dprob,drand,e12
+      real*8 dprob,drand,e12,e120
 c
       call ga_dgop(mcf_79,uda,4*nldat,'+')
 c
       if(me.eq.0) then
       if(util_print('qhop',print_high)) then
-      write(lfnhop,110) 
-  110 format(/,
+      write(lfnhop,110) stime 
+  110 format(/,' Time ',f12.6,/,
      + ' QHOP  ----donor---   --acceptor-- ',
      + '     pre-hop energies        post-hop energies  ',/,
      + ' numb atom  res prot atom  res prot',
      + '    donor      acceptor     donor      acceptor ',
-     + '     e12        distance     angle   probability',/)
+     + '     e12        e120        distance     angle   probability',/)
       endif
       do 1 i=1,nldat
-c      write(*,'(i5,5f12.6)') i,(rda(j,i),j=1,5)
+c      write(*,'(i5,8f12.6)') i,(rda(j,i),j=1,8)
       call qhop_prob(lda(1,i),lda(3,i),lseq(lda(3,i)),
-     + lda(6,i),lda(8,i),lseq(lda(8,i)),uda(1,i),rda(4,i),pda(i))
-c      write(*,'(i5,5f12.6)') i,(rda(j,i),j=1,5)
+     + lda(6,i),lda(8,i),lseq(lda(8,i)),uda(1,i),rda(4,i),pda(i),e120)
+c      write(*,'(i5,8f12.6)') i,(rda(j,i),j=1,8)
 c      pda(i)=0.0d0
 c
       if(util_print('qhop',print_high)) then
       e12=uda(3,i)+uda(4,i)-uda(1,i)-uda(2,i)
       write(lfnhop,111) i,lda(1,i),lda(3,i),lseq(lda(3,i)),
      + lda(6,i),lda(8,i),lseq(lda(8,i)),
-     + (uda(j,i),j=1,4),e12,rda(4,i),rda(5,i),pda(i)
-  111 format(7i5,7f12.6,f12.9)
+     + (uda(j,i),j=1,4),e12,e120*4.184,rda(4,i),rda(5,i),pda(i)
+  111 format(7i5,8f12.6,f12.9)
       endif
 c
     1 continue
diff --git a/src/cafe/cf_init.F b/src/cafe/cf_init.F
index 269a59e666..0f4705ed98 100644
--- a/src/cafe/cf_init.F
+++ b/src/cafe/cf_init.F
@@ -7,7 +7,7 @@
      + ipset2,issscl,delta,nfanal,lpbc,npgi,fldi,fvec,ffrq,npenrg,ictrl,
      + nbiasi,mropti,incl,ltwn,nseqi,i_lseqi,nfhopi,rhopi,thopi,ndumsi,
      + ipbtpi,lfnhopi,iradgi,nbgeti,npreci)
-c $Id: cf_init.F,v 1.108 2005-09-14 22:15:56 d3j191 Exp $
+c $Id: cf_init.F,v 1.109 2006-09-18 15:42:14 d3j191 Exp $
       implicit none
 c
 #include "cf_common.fh"
@@ -228,6 +228,11 @@ c
 c
       if(.not.ma_push_get(mt_byte,16*nsatot,'snam',l_snam,i_snam))
      + call md_abort('Failed to allocate snam',me)
+c
+      if(lqhop) then
+      if(.not.ma_push_get(mt_int,mseq,'mprot',l_mprot,i_mprot))
+     + call md_abort('Failed to allocate mprot',me)
+      endif
 c
       if(npgdec.gt.1) then
       if(.not.ma_push_get(mt_dbl,6*nsatot,'sti',l_sti,i_sti))
@@ -278,9 +283,9 @@ c     proton hopping: donor-acceptor pair list allocation
 c     ---------------------------------------------------
 c
       if(nhop.gt.0) then
-      if(.not.ma_push_get(mt_int,10*nhop*3,'lda',l_lda,i_lda))
+      if(.not.ma_push_get(mt_int,16*nhop*3,'lda',l_lda,i_lda))
      + call md_abort('Failed to allocate lda',0)
-      if(.not.ma_push_get(mt_dbl,5*nhop*3,'rda',l_rda,i_rda))
+      if(.not.ma_push_get(mt_dbl,11*nhop*3,'rda',l_rda,i_rda))
      + call md_abort('Failed to allocate rda',0)
       if(.not.ma_push_get(mt_dbl,4*nhop*3,'uda',l_uda,i_uda))
      + call md_abort('Failed to allocate uda',0)
@@ -704,7 +709,7 @@ c
 #include "mafdecls.fh"
 c
       integer ia
-      real*8 chgi(24,2)
+      real*8 chgi(24,3)
 c
       integer i
 c
@@ -730,7 +735,7 @@ c
 #include "mafdecls.fh"
 c
       integer ia
-      real*8 chgi(24,2),chg(mqt,mqp,mset)
+      real*8 chgi(24,3),chg(mqt,mqp,mset)
 c
       integer i
 c
@@ -738,11 +743,13 @@ c
       do 1 i=1,nparms
       chg(ia,1,i)=qfac*chgi(i,1)
       chg(ia,2,i)=chgi(i,2)
+      chg(ia,3,i)=qfac*chgi(i,3)
     1 continue
       else
       do 2 i=1,3
       chg(ia,1,i)=qfac*chgi(i,1)
       chg(ia,2,i)=chgi(i,2)
+      chg(ia,3,i)=qfac*chgi(i,3)
     2 continue
       do 3 i=1,2
       chg(ia,i,4)=chg(ia,i,3)-chg(ia,i,2)
@@ -1082,6 +1089,19 @@ c
       call cf_pcopy(idp,ip,p,6,mdp(id),int_mb(i_idih(id)),
      + dbl_mb(i_dih(id)),mdt(id),
      + ip2(9+(id-1)*12),ip3(9+(id-1)*12),ith(9+(id-1)*12))
+c
+      return
+      end
+      subroutine cf_parseq(id,idp)
+c
+c     cf_parseq
+c
+#include "cf_common.fh"
+#include "mafdecls.fh"
+c
+      integer id,idp
+c
+      int_mb(i_mprot+id-1)=idp
 c
       return
       end
diff --git a/src/cafe/cf_list.F b/src/cafe/cf_list.F
index 5fb4932277..6c5d528261 100644
--- a/src/cafe/cf_list.F
+++ b/src/cafe/cf_list.F
@@ -1,6 +1,6 @@
       subroutine cf_lww(lself,lpbc,xwm,idt,iwfr,iwto,jwfr,jwto,
      + nlocw,mpairs,npairs,lwwjpt,lwwin,lwwj,list,rwx,rw)
-c $Id: cf_list.F,v 1.32 2005-01-08 19:17:20 d3j191 Exp $
+c $Id: cf_list.F,v 1.33 2006-09-18 15:42:14 d3j191 Exp $
       implicit none
 c
 #include "cf_common.fh"
@@ -351,7 +351,7 @@ c
       subroutine cf_lss(lself,lpbc,lpbcs,xs,iga,idt,igr,isgm,iss,ish,
      + isz,
      + isfr,isto,jsfr,jsto,nlocs,mpairs,npairs,lssjpt,lssin,lssj,list,
-     + rwx,rw,ndxt,nt,ndxx,nx,lda,rda,lseq,lsthop)
+     + rwx,rw,ndxt,nt,ndxx,nx,lda,rda,lseq,lsthop,mprot)
 c
       implicit none
 c
@@ -368,19 +368,19 @@ c
       integer lssj(*),lssjpt(nlocs,2),lssin(nlocs,2)
       integer list(msa,3),iga(msa),igr(msa),idt(msa),iss(msa),isz(msa)
       integer ish(msa),isgm(msa)
-      integer lseq(mseq),lsthop(2,*)
-      logical lself,lpbc,lpbcs
+      integer lseq(mseq),lsthop(2,*),mprot(mseq)
+      logical lself,lpbc,lpbcs,lsame
 c
       integer isafr,isfr,isato,isto,jsafr,jsfr,jsato,jsto
       integer mpairs,npairs,isa,jsa,nlist1,nlist2,idoff,iaoff
       integer ifrsx,itosx,ifrst,itost,jfrsx,jtosx,jfrst,jtost
       integer iskip1,iskip2,istsum,ist,jst,isxsum,isx,jsx,ilist
       integer nlist0,isa0,iisa,isglob,ipss,jsglob,nlist,leq,nssl0
-      integer lda(10,*)
-      real*8 rda(5,*)
+      integer lda(16,*)
+      real*8 rda(11,*)
 c
-      real*8 xacc(5)
-      integer idon,iacc,i
+      real*8 xacc(11)
+      integer idon,iacc,i,idown,idfr,idto,iaown,iafr,iato
 c
 c
 c     this subroutine evaluates the solute-solute pairlist
@@ -391,6 +391,8 @@ c
       jsato=jsto
       if(lself) isato=isato-1
 c
+c      write(*,'(i5,a,i5)') me,' cf_lss nseq=',nseq
+c      write(*,'(20i5)') (lseq(i),i=1,nseq)
 c
       npairs=0
 c
@@ -570,11 +572,15 @@ c
       if(nhop.gt.0) then
       idon=0
       iacc=0
+      idown=0
+      iaown=0
       if(ish(isa).gt.0) then
       if(iand(ish(isa),1).eq.0) then
       idon=isa
+      idown=me
       else
       iacc=isa
+      iaown=me
       endif
       endif
       if(idon.gt.0.or.iacc.gt.0) then
@@ -585,29 +591,39 @@ c
       if(ish(jsa).gt.0) then
       if(iacc.gt.0.and.iand(ish(jsa),1).eq.0) then
       idon=jsa
+      if(lself) idown=me
       elseif(idon.gt.0.and.iand(ish(jsa),1).eq.1) then
       iacc=jsa
+      if(lself) iaown=me
       endif
       endif
       if(idon.gt.0.and.iacc.gt.0.and.rw(jsa).lt.rhop2) then
-ctmp
-c     temporary change to exclude same to same protonation state
 c
-      if(lseq(isgm(idon)).ne.lseq(isgm(iacc))) then
-ctmp
+      if(lseq(isgm(idon))+ish(idon)/2.le.mprot(idon)
+     + .or.lseq(isgm(iacc))-ish(iacc)/2.ge.1) then
 c
       do 226 i=1,nhops
       if(lsthop(1,i).eq.isgm(iacc).and.lsthop(2,i).eq.isgm(idon)) then
       goto 227
       endif
   226 continue
+c
+c
 c
       call qhop_asite(msa,xs,idon,iga(idon),isgm(idon),lseq(isgm(idon)),
-     + iacc,iga(iacc),isgm(iacc),lseq(isgm(iacc)),xacc,idoff,iaoff)
+     + iacc,iga(iacc),isgm(iacc),lseq(isgm(iacc)),xacc,idoff,iaoff,
+     + lsame)
+c
+c     do not allow hops to the same protonation state if the type of
+c     residue/molecule is identical (e.q. water dissociating into OH and H3O)
+c
+      if(.not.lsame.or.lseq(isgm(idon)).ne.lseq(isgm(iacc))) then
+c
       nlda=nlda+1
-c      write(*,2222) nlda,idon,iga(idon),isgm(idon),lseq(isgm(idon)),
-c     + iacc,iga(iacc),isgm(iacc),lseq(isgm(iacc)),sqrt(rw(jsa))
- 2222 format(9i5,f12.6)
+c      write(*,2222) me,nlda,idon,iga(idon),isgm(idon),lseq(isgm(idon)),
+c     + ish(idon),iacc,iga(iacc),isgm(iacc),lseq(isgm(iacc)),ish(iacc),
+c     + sqrt(rw(jsa))
+c 2222 format(12i5,f12.6)
       if(nlda.gt.3*nhop) call md_abort('Increase nhop',nhop)
       lda(1,nlda)=iga(idon)
       lda(2,nlda)=idon
@@ -619,19 +635,32 @@ c     + iacc,iga(iacc),isgm(iacc),lseq(isgm(iacc)),sqrt(rw(jsa))
       lda(8,nlda)=isgm(iacc)
       lda(9,nlda)=lseq(isgm(iacc))-ish(iacc)/2
       lda(10,nlda)=iaoff
+      lda(11,nlda)=idown
+      lda(12,nlda)=0
+      lda(13,nlda)=0
+      lda(14,nlda)=iaown
+      lda(15,nlda)=0
+      lda(16,nlda)=0
       rda(1,nlda)=xacc(1)
       rda(2,nlda)=xacc(2)
       rda(3,nlda)=xacc(3)
       rda(4,nlda)=xacc(4)
       rda(5,nlda)=xacc(5)
+      rda(6,nlda)=xacc(6)
+      rda(7,nlda)=xacc(7)
+      rda(8,nlda)=xacc(8)
+      rda(9,nlda)=xacc(9)
+      rda(10,nlda)=xacc(10)
+      rda(11,nlda)=xacc(11)
+c
+      endif
 c
   227 continue
 c
       endif
       endif
-ctmp
       endif
-ctmp
+c
   225 continue
       endif
       endif
diff --git a/src/cafe/cf_pairs.F b/src/cafe/cf_pairs.F
index 847f6f6f53..d8f87b48b4 100644
--- a/src/cafe/cf_pairs.F
+++ b/src/cafe/cf_pairs.F
@@ -30,7 +30,7 @@ c     in/out int : lstptr          : list pointer
 c
 c     dimensions nwm,nwa and nsa need to have been given by a call to cf_initx
 c
-c $Id: cf_pairs.F,v 1.14 2003-09-22 22:49:12 d3j191 Exp $
+c $Id: cf_pairs.F,v 1.15 2006-09-18 15:42:14 d3j191 Exp $
       implicit none
 c
 #include "cf_common.fh"
@@ -142,7 +142,8 @@ c
      + int_mb(lptr+1+4*nsloc),
      + int_mb(i_s2i1),dbl_mb(i_s3r1),dbl_mb(i_s1r1),
      + int_mb(i_itrd(2)),mtt(2),int_mb(i_ixcl(2)),mxt(2),
-     + int_mb(i_lda),dbl_mb(i_rda),lseq,int_mb(i_lsthop))
+     + int_mb(i_lda),dbl_mb(i_rda),lseq,int_mb(i_lsthop),
+     + int_mb(i_mprot))
       endif
       int_mb(lptr)=npairs
       int_mb(lptrn+3)=lptr
diff --git a/src/cafe/cf_rdtop.F b/src/cafe/cf_rdtop.F
index 0ba964d23d..b62b7fb8ce 100644
--- a/src/cafe/cf_rdtop.F
+++ b/src/cafe/cf_rdtop.F
@@ -1,5 +1,5 @@
       subroutine cf_rdtop(lfntop,filtop,snam)
-c $Id: cf_rdtop.F,v 1.35 2005-01-08 20:27:06 d3j191 Exp $
+c $Id: cf_rdtop.F,v 1.36 2006-09-18 15:42:14 d3j191 Exp $
       implicit none
 c
 #include "cf_common.fh"
@@ -15,14 +15,14 @@ c
       real*8 releps
       character*16 stemp,snam(nsatot)
       integer natyps,nqtyps
-      real*8 rdata(24),vdata(24,4),qdata(24,2),pdata(4,24)
+      real*8 rdata(24),vdata(24,4),qdata(24,3),pdata(4,24)
       integer idata(24)
       character*6 cdata(24)
       integer i,ia,iq,j,k,l,m,m2,m3,m4,itemp(24),i_i,i_j,l_i,l_j
       integer naw,nbw,nhw,ndw,now,ntw,nnw
       integer nhp
       integer nas,nbt,nhs,nds,nos,nts,nxs
-      integer numbl,numbi,numbd,iqfr,iqto,nsad
+      integer numbl,numbi,numbd,iqfr,iqto,nsad,maxpro
       character*10 string
 c
       real*8 cf_charge
@@ -89,16 +89,17 @@ c
       endif
       do 4 i=1,nqtyps
       read(lfntop,1007,end=9997,err=9998)
-     + ((qdata(k,l),l=1,2),k=1,nparms)
- 1007 format(5x,f12.6,e12.5)
+     + ((qdata(k,l),l=1,3),k=1,nparms)
+ 1007 format(5x,f12.6,e12.5,f12.6)
       call cf_parq(i,qdata)
     4 continue
       if(.not.ma_verify_allocator_stuff()) then
       call md_abort('ERROR IN MEMORY IN RDTOP CHARGE TYPES',0)
       endif
       do 4004 i=1,nseq
-      read(lfntop,4008,end=9997,err=9998) string
- 4008 format(a10)
+      read(lfntop,4008,end=9997,err=9998) string,maxpro
+ 4008 format(a10,33x,i5)
+      call cf_parseq(i,maxpro)
  4004 continue
       read(lfntop,1008,end=9997,err=9998) naw,nbw,nhw,ndw,now,ntw,nnw
  1008 format(5i7,2i10)
diff --git a/src/cafe/cf_shake.F b/src/cafe/cf_shake.F
index 759a2de166..f03fa270df 100644
--- a/src/cafe/cf_shake.F
+++ b/src/cafe/cf_shake.F
@@ -1,5 +1,5 @@
       subroutine cf_shakew(xw,yw,iwgm,iwdt,numw)
-c $Id: cf_shake.F,v 1.26 2003-10-19 03:30:56 d3j191 Exp $
+c $Id: cf_shake.F,v 1.27 2006-09-18 15:42:15 d3j191 Exp $
       implicit none
 c
 #include "cf_common.fh"
@@ -514,7 +514,21 @@ c
       endif
       endif
       endif
+c
+c     place dummy hydrogens (acceptor sites) onto bound heavy atoms
+c
+      if(iand(ishp(isa),1).eq.1.and.iand(ishp(jsa),1).ne.1) then
+      xs(isa,1)=xs(jsa,1)
+      xs(isa,2)=xs(jsa,2)
+      xs(isa,3)=xs(jsa,3)
       endif
+      if(iand(ishp(isa),1).ne.1.and.iand(ishp(jsa),1).eq.1) then
+      xs(jsa,1)=xs(isa,1)
+      xs(jsa,2)=xs(isa,2)
+      xs(jsa,3)=xs(isa,3)
+      endif
+      endif
+c
     2 continue
 c
       niters=niters+1
diff --git a/src/nwmd/md_main.F b/src/nwmd/md_main.F
index 3b06bcb639..ce605c5d51 100644
--- a/src/nwmd/md_main.F
+++ b/src/nwmd/md_main.F
@@ -1,6 +1,6 @@
       subroutine md_main()
 c
-c $Id: md_main.F,v 1.211 2006-06-05 00:12:29 marat Exp $
+c $Id: md_main.F,v 1.212 2006-09-18 15:42:03 d3j191 Exp $
 c
       implicit none
 c
@@ -97,19 +97,17 @@ c
      + epot,epotw,epotsw,epots,volume,dwr,dbl_mb(i_dsr),ekin,etot,
      + npolit,dbl_mb(i_gsm),dbl_mb(i_esa),box,dbl_mb(i_xsm))
 c
-      if(.not.lqmmm) 
-     >       call prp_print()
+      if(.not.lqmmm) call prp_print()
 c
       call rtdb_put(irtdb,'md:energy',mt_dbl,1,epot)
 c
 c     print energies
 c
-      if(.not.lqmmm) 
-     >  call cf_print_energy(lfnout)
-c
-      if(.not.lqmmm) 
-     >       call sp_printf(filtop,lfntop,
+      if(.not.lqmmm) then
+      call cf_print_energy(lfnout)
+      call sp_printf(filtop,lfntop,
      + int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_fs),npener,dbl_mb(i_esa))
+      endif
 c
       if(ifidi.ne.0) then
       call md_fd(int_mb(i_is),dbl_mb(i_xs),dbl_mb(i_vs),dbl_mb(i_fs),
@@ -2481,8 +2479,8 @@ c
       endif
 c
       if(.not.ldone) then
-      if(cf_hopping(stime,isl(1,lssgm),isl(1,lsgan),isl(1,lshop),
-     + xs,nsaloc,lfnhop)) then
+      if(cf_hopping(lpbc,lpbcs,stime,isl,isl(1,lssgm),isl(1,lsgan),
+     + isl(1,lsct3),isl(1,lshop),xs,nsaloc,lfnhop)) then
       if(lpair) then
       lpair=.false.
       lforce=.true.
@@ -2902,7 +2900,7 @@ c
       end
       subroutine md_test
 c
-c $Id: md_main.F,v 1.211 2006-06-05 00:12:29 marat Exp $
+c $Id: md_main.F,v 1.212 2006-09-18 15:42:03 d3j191 Exp $
 c
       implicit none
 c
diff --git a/src/qhop/qhop_asite.F b/src/qhop/qhop_asite.F
index 0207ee575f..ee521b8633 100644
--- a/src/qhop/qhop_asite.F
+++ b/src/qhop/qhop_asite.F
@@ -1,7 +1,7 @@
       subroutine qhop_asite(msa,xs,idon,idong,resd,istatd,iacc,
-     + iaccg,resa,istata,xacc,doff,aoff)
+     + iaccg,resa,istata,xacc,doff,aoff,flada)
 c
-c $Id: qhop_asite.F,v 1.3 2003-08-11 19:43:18 d3j191 Exp $
+c $Id: qhop_asite.F,v 1.4 2006-09-18 15:45:00 d3j191 Exp $
 c
       implicit none
 c
@@ -11,11 +11,11 @@ c
 c     declaration of arguments
 c
       integer msa,idon,idong,resd,iacc,iaccg,resa,istatd,istata,doff,
-     + aoff
-      real*8 xacc(5),xs(msa,3)
+     + aoff,flada
+      real*8 xacc(11),xs(msa,3)
 c
       call qhop_asite2(msa,xs,idon,idong,resd,istatd,iacc,
-     + iaccg,resa,istata,xacc,doff,aoff,
+     + iaccg,resa,istata,xacc,doff,aoff,flada,
      + int_mb(i_ptseq),int_mb(i_ptarat),int_mb(i_iarat),
      + dbl_mb(i_racs),dbl_mb(i_racs+mxar),int_mb(i_ptpar),
      + dbl_mb(i_par),dbl_mb(i_par+7*maxpar),dbl_mb(i_par+14*maxpar),
@@ -25,11 +25,11 @@ c
       return
       end
       subroutine qhop_asite2(msa,xs,idon,idong,resd,istatd,iacc,
-     + iaccg,resa,istata,xacc,doff,aoff,
+     + iaccg,resa,istata,xacc,doff,aoff,flada,
      + ptseq,ptarat,arat,deq,aneq,ptpar,tunnel,defe12,zpef,
      + tdsgl,tstval,e12fxy)
 c
-c $Id: qhop_asite.F,v 1.3 2003-08-11 19:43:18 d3j191 Exp $
+c $Id: qhop_asite.F,v 1.4 2006-09-18 15:45:00 d3j191 Exp $
 c
       implicit none
 c
@@ -38,8 +38,8 @@ c
 c     declaration of arguments
 c
       integer msa,idon,idong,resd,iacc,iaccg,resa,istatd,istata,doff,
-     + aoff
-      real*8 xacc(5),xs(msa,3)
+     + aoff,flada
+      real*8 xacc(11),xs(msa,3)
 c
       real*8 tunnel(maxpar,7),defe12(maxpar,7),zpef(maxpar,3),
      + tdsgl(maxpar,6),tstval(maxpar,5),e12fxy(maxpar,3)
@@ -49,15 +49,15 @@ c
 c     declaration of variables
 c
       integer typ,i,idr(2),idhv(2),inh,inhv,numh,nhv(2)
-      real*8 dist,ang
-      real*8 cost,sint,cosp,sinp,th,ph
+      real*8 dist,ang,d2
+      real*8 cost,sint,cosp,sinp,th,ph,alpha
       real*8 trsl,rds,fct
       real*8 xh2(3),xacc2(2),xdr(3),xr12(3),xar(3),xdrp(2),norm,r1(3),
-     + sid1(3),sid2(3),nsid(2),xrh(3),nrh
+     + sid1(3),sid2(3),nsid(2),xrh(3),nrh,r1pbc(3)
       real*8 tol,pi
 c
-      parameter(tol=1E-10)
-      parameter(pi=3.14159)      
+      parameter(tol=1E-15)
+      parameter(pi=3.14159265)      
 c
       typ=arat(ptarat(ptseq(resa,1))+iaccg-
      + ptseq(resa,2),2)  
@@ -72,7 +72,10 @@ c
      + ptseq(resa,2))
       dist=deq(ptarat(ptseq(resa,1))+iaccg-
      + ptseq(resa,2))
+c      write (*,98),'angdi ',ang,dist
+c 98   format(a6,2f10.6)
 c
+      ang=1.910610
       idhv(1)=arat(ptarat(ptseq(resd,1))+idong-
      + ptseq(resd,2),1)+idon
       idhv(2)=arat(ptarat(ptseq(resa,1))+iaccg-
@@ -102,21 +105,13 @@ c
       rds=dist*sin(pi-ang)
 c  
       if(xdr(1).lt.tol.and.xdr(2).lt.tol)then
-c   
-c     no rotation needed
-c        
-      fct=sqrt(xh2(1)**2+xh2(2)**2)
-      if(fct.gt.tol)then
-      xacc(1)=xs(iacc,1)*(rds/fct)
-      xacc(2)=xs(iacc,2)*(rds/fct)
-      else
 c
 c     H projection on the circle plane coincides with the center of 
 c     the circle, all points are equally distant
 c     
-      xacc(1)=rds
-      xacc(2)=0.
-      endif
+      xacc(1)=rds+xs(idhv(2),1)
+      xacc(2)=xs(idhv(2),2)
+      xacc(3)=xs(idhv(2),3)+sign(trsl,xs(idhv(2),3)-xs(idr(1),3))
       return
       else
       ph=atan(-(xdr(2)/(xdr(1))))
@@ -196,8 +191,126 @@ c
 c
 c     H2O: idr(1) and idr(2) are the indexes of the two non dummy H
 c
-      else       
-c 
+      else if(typ.eq.4)then      
+c
+c     if acceptor is water then fix exactly the 2 angles
+c
+      xdr(1)=xs(idhv(2)+1,1)-xs(idhv(2),1)
+      xdr(2)=xs(idhv(2)+1,2)-xs(idhv(2),2)
+      xdr(3)=xs(idhv(2)+1,3)-xs(idhv(2),3)     
+      call cf_pbc(1,xdr,1,xdr,1,0,1,1)
+c      write (*,91),'xdr1 ',xdr
+      
+      if(abs(xdr(1)).gt.tol)then                                  
+         if(abs(xdr(2)).gt.tol)then
+            th=atan(-xdr(1)/xdr(2))
+c            write (*,92),'theta ',th,xdr(1),xdr(2)
+ 92         format(a6,3f10.6)
+         else
+            th=-sign(1.,xdr(1))*(pi/2)
+c            write *,'thetanotol ',xdr(1)-tol,xdr(2)-tol
+         endif
+      else
+         th=0                                              
+      endif 
+      xdrp(1)=-xdr(1)*sin(th)+xdr(2)*cos(th)                     
+                                                              
+      if(abs(xdrp(1)).gt.tol)then                               
+         if(abs(xdr(3)).gt.tol)then                              
+            ph=atan(-xdrp(1)/xdr(3)) 
+c            write (*,92),'phi ',ph,xdrp(1),xdr(3)
+         else                                                 
+            ph=-sign(1.,xdrp(1))                      
+c            write *,'phinotol ',xdrp(1)-tol,xdr(3)-tol
+         endif                                                
+      else                                                    
+         ph=0                                                
+      endif                                                   
+      xdrp(2)=-xdrp(1)*sin(ph)+xdr(3)*cos(ph)
+      if(xdrp(2).le.0)ph=ph+pi
+
+      xdr(1)=xs(idhv(2)+2,1)-xs(idhv(2),1)
+      xdr(2)=xs(idhv(2)+2,2)-xs(idhv(2),2)
+      xdr(3)=xs(idhv(2)+2,3)-xs(idhv(2),3)
+      call cf_pbc(1,xdr,1,xdr,1,0,1,1)
+
+      call rot(th,xdr(1),xdr(2),xdr(3),3)
+      call rot(ph,xdr(1),xdr(2),xdr(3),1)
+c      write (*,91),'xdr2 ',xdr
+      if(abs(xdr(2)).gt.tol)then
+         if(abs(xdr(1)).gt.tol)then
+            alpha=atan(xdr(2)/xdr(1))
+c            write (*,92),'alpha ',alpha,xdr(2),xdr(1)
+         else
+            alpha=sign(1.,xdr(2))
+         endif
+      else
+         alpha=0
+      endif
+      xdrp(1)=xdr(1)*cos(alpha)+xdr(2)*sin(alpha)
+      if(xdrp(1).le.0)then
+         if(alpha.le.0)then
+            alpha=alpha+pi
+         else
+            alpha=alpha-pi
+         endif
+      endif
+c
+c     Found the angles, perform rotation
+c
+      xdr(1)=xs(idhv(1),1)-xs(idhv(2),1)
+      xdr(2)=xs(idhv(1),2)-xs(idhv(2),2)
+      xdr(3)=xs(idhv(1),3)-xs(idhv(2),3)
+      call cf_pbc(1,xdr,1,xdr,1,0,1,1)
+
+      call rot(th,xdr(1),xdr(2),xdr(3),3)
+      call rot(ph,xdr(1),xdr(2),xdr(3),1)
+      call rot(alpha,xdr(1),xdr(2),xdr(3),3)
+      xr12(1)=xdr(1)
+      xr12(2)=xdr(2)
+      xr12(3)=xdr(3)
+      xdr(1)=xs(idhv(2)+2,1)-xs(idhv(2),1)
+      xdr(2)=xs(idhv(2)+2,2)-xs(idhv(2),2)
+      xdr(3)=xs(idhv(2)+2,3)-xs(idhv(2),3)
+      call cf_pbc(1,xdr,1,xdr,1,0,1,1)
+
+      call rot(th,xdr(1),xdr(2),xdr(3),3)
+      call rot(ph,xdr(1),xdr(2),xdr(3),1)
+      call rot(alpha,xdr(1),xdr(2),xdr(3),3)
+c      write (*,91),'xdr3 ',xdr
+      d2=dist**2
+      xacc(3)=dist*cos(ang)
+      xacc(1)=(1/xdr(1))*xacc(3)*(0.1-xdr(3))
+      xacc(2)=sqrt(d2-xacc(3)**2-xacc(1)**2)
+c      write (*,93),'xacc1 ',xacc(1),xacc(2),xacc(3),d2,cos(ang)
+c      write(*,94),'norm ',xacc(1)**2+xacc(2)**2+xacc(3)**2
+      d2=abs(xr12(2)-xacc(2))-abs(xr12(2)+xacc(2))
+      if(d2.ge.0)xacc(2)=-xacc(2)
+      xdr(1)=xacc(1)
+      xdr(2)=xacc(2)
+      xdr(3)=xacc(3)
+      call rot(-alpha,xdr(1),xdr(2),xdr(3),3)
+      call rot(-ph,xdr(1),xdr(2),xdr(3),1)
+      call rot(-th,xdr(1),xdr(2),xdr(3),3)
+c      write (*,94),'norm ',xacc(1)**2+xacc(2)**2+xacc(3)**2
+      xacc(1)=xdr(1)+xs(idhv(2),1)
+      xacc(2)=xdr(2)+xs(idhv(2),2)
+      xacc(3)=xdr(3)+xs(idhv(2),3)
+c      write (*,93),'xacc2 ',xacc(1),xacc(2),xacc(3),d2,cos(ang)
+ 91   format(a5,3f10.6)
+ 93   format(a6,5f10.6)
+ 94   format(a5,f10.6)
+
+      else if(typ.eq.5)then
+         
+         xar(1)=xs(idon,1)-xs(idhv(2),1)
+         xar(2)=xs(idon,2)-xs(idhv(2),2)
+         xar(3)=xs(idon,3)-xs(idhv(2),3)
+         call cf_pbc(1,xar,1,xar,1,0,1,1)
+         norm=sqrt(xar(1)**2+xar(2)**2+xar(3)**2)
+         xacc(1)=xs(idhv(2),1)+dist*(xar(1)/norm)
+         xacc(2)=xs(idhv(2),2)+dist*(xar(2)/norm)
+         xacc(3)=xs(idhv(2),3)+dist*(xar(3)/norm)
       endif
 c
 c     Calculating angle and donor-acceptor dist
@@ -205,8 +318,37 @@ c
       r1(1)=xs(idhv(1),1)-xs(idhv(2),1)
       r1(2)=xs(idhv(1),2)-xs(idhv(2),2)
       r1(3)=xs(idhv(1),3)-xs(idhv(2),3)
+      r1pbc(1)=r1(1)
+      r1pbc(2)=r1(2)
+      r1pbc(3)=r1(3)
       call cf_pbc(1,r1,1,r1,1,0,1,1)
       xacc(4)=sqrt(r1(1)**2+r1(2)**2+r1(3)**2)
+c  
+c     calculating midpoint of donor-acceptor
+c
+      if(abs(r1pbc(1)-r1(1)).lt.0.001) then
+      xacc(6)=(xs(idhv(1),1)+xs(idhv(2),1))/2      
+      xacc(9)=xacc(6)
+      else
+      xacc(6)=(xs(idhv(1),1)+(r1pbc(1)-r1(1))+xs(idhv(2),1))/2
+      xacc(9)=xacc(6)-(r1pbc(1)-r1(1))
+      endif
+c
+      if(abs(r1pbc(2)-r1(2)).lt.0.001) then
+      xacc(7)=(xs(idhv(1),2)+xs(idhv(2),2))/2      
+      xacc(10)=xacc(7)
+      else
+      xacc(7)=(xs(idhv(1),2)+(r1pbc(2)-r1(2))+xs(idhv(2),2))/2
+      xacc(10)=xacc(7)-(r1pbc(2)-r1(2))
+      endif
+c
+      if(abs(r1pbc(3)-r1(3)).lt.0.001) then
+      xacc(8)=(xs(idhv(1),3)+xs(idhv(2),3))/2      
+      xacc(11)=xacc(8)
+      else
+      xacc(8)=(xs(idhv(1),3)+(r1pbc(3)-r1(3))+xs(idhv(2),3))/2
+      xacc(11)=xacc(8)-(r1pbc(3)-r1(3))
+      endif
 c
       sid1(1)=xs(idhv(1),1)-xs(idon,1)
       sid1(2)=xs(idhv(1),2)-xs(idon,2)
@@ -224,6 +366,15 @@ c
      +         (1/(nsid(1)*nsid(2)))
       xacc(5)=acos(xacc(5))
 c
+c
+c     debug: print shift of proton
+c
+      r1(1)=xs(idon,1)-xacc(1)
+      r1(2)=xs(idon,2)-xacc(2)
+      r1(3)=xs(idon,3)-xacc(3)
+      call cf_pbc(1,r1,1,r1,1,0,1,1)
+c      dshop=sqrt(r1(1)**2+r1(2)**2+r1(3)**2)
+c      print*,"mov ",resd,idong,resa,norm,xacc(4),xacc(5)
 c     Shift for permutation of donor H
 c
       numh=arat(ptarat(ptseq(resd,1))+idong-
@@ -232,8 +383,15 @@ c
      + ptseq(resd,2),1)
       inhv=arat(ptarat(ptseq(resd,1))+idong-
      + ptseq(resd,2),3)
-      doff=numh-inh-iand(((istatd-1)/(2**(inhv-1))+1),1)
-      if(typ.eq.3)doff=numh-inh-(istatd-1)
+      if(numh.eq.1)then
+      doff=0
+      else if(numh.eq.2)then
+      doff=numh-inh
+c      doff=numh-inh-iand(((istatd-1)/(2**(inhv-1))+1),1)
+      else
+      doff=numh-inh-(istatd-1)
+      endif
+c      print*,'don ',resd,inh,doff,istatd
 c
 c     Shift for permutation of acceptor dummy site
 c     Currently if nhv>1,only one dummy is allowed
@@ -246,11 +404,39 @@ c
      + ptseq(resa,2),3)
       aoff=inh-numh-1+iand(((istata-1)/(2**(inhv-1))+1),1)
       aoff=0
-      if(typ.eq.3)aoff=numh-inh-(istata-2)
-c      print*,numh,inh,istata,aoff
+c      if(numh>2)aoff=numh-inh-(istata-2)
+c     print*,'acc ',resa,inh,doff,istata
 c
 c      print 1021,typ,doff,aoff,xacc,numh,inh,istatd
- 1021 format(3i5,5f8.3,3i5)
+c 1021 format(3i5,5f8.3,3i5)
+c
+c      flag to see if res are same type
+c
+      flada=0
+      if(ptseq(resa,1).eq.ptseq(resd,1))flada=1
       return
       end
 
+      subroutine rot(angle,x,y,z,axis)
+
+      implicit none
+
+      real*8 angle,x,y,z,x1,y1,z1
+      integer axis
+
+      x1=x
+      y1=y
+      z1=z
+
+      if(axis.eq.3)then
+         x=x1*cos(angle)+y1*sin(angle)
+         y=-x1*sin(angle)+y1*cos(angle)
+      else if(axis.eq.1) then
+         y=y1*cos(angle)+z1*sin(angle)
+         z=-y1*sin(angle)+z1*cos(angle)
+      else
+
+      endif
+      end
+
+
diff --git a/src/qhop/qhop_prob.F b/src/qhop/qhop_prob.F
index 9892b51e71..f590cefa2e 100644
--- a/src/qhop/qhop_prob.F
+++ b/src/qhop/qhop_prob.F
@@ -1,4 +1,4 @@
-      subroutine qhop_prob(gd,rsd,prd,ga,rsa,pra,uda,dang,prob)
+      subroutine qhop_prob(gd,rsd,prd,ga,rsa,pra,uda,dang,prob,e120)
 c     
       implicit none
 c
@@ -7,7 +7,7 @@ c
 c
 c     declaration of arguments
 c
-      real*8 uda(4),dang(2),prob
+      real*8 uda(4),dang(2),prob,e120
       integer gd,rsd,prd,ga,rsa,pra
 c
       call qhop_prob2(gd,rsd,prd,ga,rsa,pra,uda,dang,prob,
@@ -15,15 +15,17 @@ c
      + dbl_mb(i_racs),dbl_mb(i_racs+mxar),int_mb(i_ptpar),
      + dbl_mb(i_par),dbl_mb(i_par+7*maxpar),dbl_mb(i_par+14*maxpar),
      + dbl_mb(i_par+17*maxpar),dbl_mb(i_par+23*maxpar),
-     + dbl_mb(i_par+28*maxpar))
+     + dbl_mb(i_par+28*maxpar),e120)
+c
+      e120=e120*4.184
 c
       return
       end
       subroutine qhop_prob2(gd,rsd,prd,ga,rsa,pra,uda,dang,prob,
      + ptseq,ptarat,arat,deq,aneq,ptpar,tunnel,defe12,zpef,
-     + tdsgl,tstval,e12fxy)
+     + tdsgl,tstval,e12fxy,e120)
 c
-c $Id: qhop_prob.F,v 1.6 2003-09-22 15:20:10 d3j191 Exp $
+c $Id: qhop_prob.F,v 1.7 2006-09-18 15:45:00 d3j191 Exp $
 c
 c     
 c     The hopping probability prob is calculated given the difference in 
@@ -46,7 +48,7 @@ c
 c
 c     declaration of arguments
 c
-      real*8 uda(4),dang(2),prob
+      real*8 uda(4),dang(2),prob,e120
       integer gd,rsd,prd,ga,rsa,pra
 c
       real*8 tunnel(maxpar,7),defe12(maxpar,7),zpef(maxpar,3),
@@ -100,9 +102,8 @@ c 1026 format('prob ',9i5)
 c
 c     Calculating e12 
 c
-      e12=(uda(3)+uda(4)-uda(1)-uda(2))/4.184
-      e12=e12+e12fxy(ipar,1)+e12fxy(ipar,2)*dist+
-     + e12fxy(ipar,3)*(dist**2)
+      e120=e12fxy(ipar,1)+e12fxy(ipar,2)*dist+e12fxy(ipar,3)*(dist**2)
+      e12=e12+e120+(uda(3)+uda(4)-uda(1)-uda(2))/4.184
 c      
 c      lhv(1)=ld+arat(ptarat(ptseq(iseq(1),1)+gd-ptseq(iseq(1),2),1)
 c      lhv(2)=la+arat(ptarat(ptseq(iseq(2),1)+ga-ptseq(iseq(2),2),1)
diff --git a/src/space/GNUmakefile b/src/space/GNUmakefile
index e5637ef38a..b9fd0863dd 100644
--- a/src/space/GNUmakefile
+++ b/src/space/GNUmakefile
@@ -1,5 +1,5 @@
 #
-# $Id: GNUmakefile,v 1.3 2006-08-25 20:24:45 d3p307 Exp $
+# $Id: GNUmakefile,v 1.4 2006-09-18 15:44:09 d3j191 Exp $
 #
 #
 OBJ_OPTIMIZE=sp_init.o sp_main.o
@@ -16,15 +16,4 @@ ifeq ($(TARGET),LINUX)
 DEFINES += -DNEED_LOC
 endif
 
-ifeq ($(FC),ifc)
-FOPTIMIZE += -Ob0
-endif
-
-ifeq ($(FC),ifort)
- ifeq ($(_CPU),x86_64)
-  FOPTIMIZE += -O0
- endif
-endif
-
-
 include ../config/makelib.h
diff --git a/src/space/sp_main.F b/src/space/sp_main.F
index 9bc2fefc23..8527176f98 100644
--- a/src/space/sp_main.F
+++ b/src/space/sp_main.F
@@ -1,7 +1,7 @@
       subroutine sp_setup(node,iwl,xw,xwcr,vw,gw,numwm,
      + isl,xs,xscr,vs,gs,numsa,lpack)
 c
-c $Id: sp_main.F,v 1.97 2006-03-07 19:44:28 d3j191 Exp $
+c $Id: sp_main.F,v 1.98 2006-09-18 15:44:09 d3j191 Exp $
 c
       implicit none
 c
@@ -2866,6 +2866,65 @@ c
      + int_mb(i_bb),int_mb(i_ipl),int_mb(i_jpl),lself,.false.)
 c
       lpbcs=lpbc9
+c
+      return
+      end
+      subroutine sp_getda(jnode,jsfr,jsto,isl,isga,isq3,xs,ndxda)
+c
+      implicit none
+c
+#include "sp_common.fh"
+#include "mafdecls.fh"
+c
+      integer jnode,jsfr,jsto,ndxda
+      integer isga(msa),isq3(msa),isl(msa,mis2)
+      real*8 xs(msa,3)
+c
+      integer il,ih,jl,jh
+      integer l_islt,i_islt
+c
+      if(.not.ma_push_get(mt_int,msa*mis2,'islt',l_islt,i_islt))
+     + call md_abort('Failed to allocate memory for islt',0)
+      call sp_gtda(jnode,jsfr,jsto,int_mb(i_islt),int_mb(i_pack),
+     + isga,isq3,xs,ndxda)
+      if(.not.ma_pop_stack(l_islt))
+     + call md_abort('Failed to deallocate memory for islt',0)
+c
+      return
+      end
+      subroutine sp_gtda(jnode,jsfr,jsto,isl,islp,isga,isq3,xs,ndxda)
+c
+      implicit none
+c
+#include "sp_common.fh"
+#include "global.fh"
+c
+      integer jnode,jsfr,jsto,ndxda
+      integer isl(msa,mis2),islp(msa,npack)
+      integer isga(msa),isq3(msa)
+      real*8 xs(msa,3)
+c
+      integer i,il,ih,jl,jh
+c
+      if(npack.eq.0) then
+      call ga_distribution(ga_is,jnode,il,ih,jl,jh)
+      call ga_get(ga_is,il+jsfr-1,il+jsto-1,jl,jl,
+     + isga(ndxda+1),msa)
+      call ga_get(ga_is,il+jsfr-1,il+jsto-1,jl+lsct3-1,jl+lsct3-1,
+     + isq3(ndxda+1),msa)
+      else
+      call ga_distribution(ga_is,jnode,il,ih,jl,jh)
+      call ga_get(ga_is,il+jsfr-1,il+jsto-1,jl,jl+npack-1,
+     + islp(nsaloc+1,1),msa)
+      call sp_unpack(jsto-jsfr+1,isl(nsaloc+1,1),islp(nsaloc+1,1))
+      do 1 i=1,jsto-jsfr
+      isga(ndxda+i)=isl(nsaloc+i,lsgan)
+      isq3(ndxda+i)=isl(nsaloc+i,lsct3)
+    1 continue
+      endif
+      call ga_distribution(ga_s,jnode,il,ih,jl,jh)
+      call ga_get(ga_s,il+jsfr-1,il+jsto-1,jl,jl+2,
+     + xs(ndxda+1,1),msa)
 c
       return
       end