      program test_potential
      implicit none
      integer iout, iprint, iprint1, np, i, j,  k
      parameter(iout=6)
      real*8 coords(7), ang(6), r, pi,e
c
      iprint=100
      iprint1=2
      pi=dacos(-1d0)
      do k=1,4
         if (k.eq.1) then
            call sapt5s_init
            call ss_print(iout,iprint)
         elseif (k.eq.2) then
            call sapt5st_init
            call ss_print(iout,iprint)
         elseif (k.eq.3) then
            call sapt5s_init
            call ss_dshift
            call ss_print(iout,iprint1)
         else
            call sapt5st_init
            call ss_dshift
	    call ss_print(iout,iprint1)
         endif

         read(*,*) np    ! the number of points
         write(*,*)
         write(*,*)'Number of points = ',np
         write(iout,'(a11,7a13)') 'alpha1(rad)', 'beta_1','gamma_1',
     +        'alpha_2','beta_2','gamma_2','R (bohr)','V (hartree)'
         do i=1,np
            read(*,*) r,(ang(j),j=1,6)
            do j=1,6
               coords(j)=ang(j)*pi/180d0    ! convert from degrees to radians
            end do
            coords(7)=r
            call ss_eval(coords,e)
            write(iout,'(7(e12.5,1x),e15.8)')coords,e
         end do
         rewind(5)
         write(iout,*)
      end do
      end
c
c***********************************************************************
c
      subroutine sapt5s_init
      implicit none
c
c  Initialize SAPT5-s potential. Use SS_EVAL to evaluate the potential.
c
c  Monomer geometry (in x-z plane, O on positive z-axis)
c
c  rOH  =  1.836106337 bohr
c  aHOH =  104.69      degrees
c
c  Sites
c  1:  O
c  2:  H1
c  3:  H2
c  4:  Bunny2 1 ("lone pair repulsion")
c  5:  Bunny2 2
c  6:  M        Center of mass for mass_O=16, mass_H=1
c  7:  Bunny1 1 (with charge)
c  8:  Bunny1 2
c
c     Integer parameters
c
      integer nip
      parameter(nip=12+2*( 5 + 6 + 3 )) ! (12+2*(nq+nr+nl))
      integer ip(nip)
      data ip/
     + 1,          ! Version
     + 8,          ! Number of sites (nsite)
     + 5,          ! Number of sites with charges (nq)   
     + 3,          ! Different types of sites with charges (mq)
     + 6,          ! Number of sites with repulsion (nr)
     + 4,          ! Different types of sites with repulsion (mr)
     + 4,          ! Order of polynomial for repulsive part (mrorder)
     + 3,          ! Number of sites with long range (nl)
     + 2,          ! Different types of sites with long range (ml)
     + 3,6,2,      ! nc,nlow,nstep:   n=nlow+(0:nc-1)*nstep, R^-n

     + 1,2,3,7,8,  ! listq(nq), charges on O, H1, H2, B1_1, B1_2
     + 1,2,2,3,3,  ! itypeq(nq), 1=O, 2=H, 3=B1

     + 1,2,3,6,4,5,! listr(nr), repulsion on O, H1, H2, M, B2_1, B2_1
     + 1,2,2,3,4,4,! ityper(nr), 1=O, 2=H, 3=M, 4=B2

     + 1,2,3,      ! listl(nl), O, H1, H2
     + 1,2,2/      ! itypel(nl), 1=O, 2=H
c
c     Real*8 parameters
c
      integer np
      parameter(np=3*8+3+6+10+10*4+3*3+3*3)
      real*8 p(np)
      data p/
c
c     cartesian coordinates of the sites
c
     +     0.00000000000000d+00,     0.00000000000000d+00,  ! O
     +     1.24631924913843d-01,    -1.45365196228170d+00,  !   H1
     +     0.00000000000000d+00,    -9.97055399310745d-01,
     +     1.45365196228170d+00,     0.00000000000000d+00,  ! H2
     +    -9.97055399310745d-01,     0.00000000000000d+00,  !   B2
     +    -2.00000000000000d-01,     2.50000000000000d-01,
     +     0.00000000000000d+00,     2.00000000000000d-01,  ! B2
     +     2.50000000000000d-01,     0.00000000000000d+00,  !   M
     +     0.00000000000000d+00,     0.00000000000000d+00,  !
     +     0.00000000000000d+00,    -0.2067213d0         ,  ! B1
     +    -0.247160511d0,            0.00000000000000d+00,  !   B1
     +     0.2067213d0,             -0.247160511d0,
     +     0.258504607653328389d0,   0.564050362569283648d0,! charge
     +    -0.693302666395947953d0,
     +    5.30188399705912605d0,  ! betaq
     +    1.40264620725287381d0,
     +    1.37777021367524322d0,
     +    4.88257857151977070d0,
     +    1.45060834365989932d0,
     +    1.6247997153704614d0,
     +    3.96533329909905d0,     ! betar
     +    3.31344984120762d0,
     +    2.41429207550821d0,
     +    3.97675261600551d0,
     +    3.17256531346957d0,
     +    3.98715016325772d0,
     +    3.96341198097565d0,
     +    3.42171736051266d0,
     +    3.96806674259008d0,
     +    3.96112259788312d0,
     +    2.29382210802996d+06,   ! arep(10,4)
     +    1.83391990344517d+03,
     +    1.55448423438490d+03,
     +    7.07598200679908d+02,
     +    1.46465352548225d+02,
     +    4.00984358322281d+01,
     +    7.35714039293257d+04,
     +    1.74134617587565d+04,
     +    7.32996751673306d+02,
     +    2.61027567540975d+02,
c Raw coefficients:
c a1
     +    4.24349609210067d+01,
     +    5.53348022483608d+01,
     +   -3.33594721363884d-01,
     +   -7.58030380141890d+04,
     +   -2.20441451218667d+02,
     +    6.98033826871118d+05,
     +   -3.45744186471392d+02,
     +   -3.31423873332427d+00,
     +    1.75746806682120d+04,
     +    2.38407041000669d+04,
c a2
     +   -5.79143175436583d+01,
     +   -4.90389806476621d+01,
     +   -2.08430904465023d-02,
     +    1.00152026628681d+05,
     +    1.52610558646430d+02,
     +   -9.21648631994710d+05,
     +    4.58827725477179d+02,
     +    2.62590835009244d+00,
     +   -2.33976987314305d+04,
     +   -3.19575522797722d+04,
c a3
     +    2.18222940530393d+01,
     +    2.89395321747882d+01,
     +    6.74078416427915d-02,
     +   -3.79108764417642d+04,
     +   -1.14235392332530d+02,
     +    3.54509635717971d+05,
     +   -1.70081179834320d+02,
     +   -1.16337257478537d+00,
     +    8.85278685268623d+03,
     +    1.18066845162708d+04,
     +  -418.942352042298694d0, ! cn(3,3)
     +   410.613776430728421d0,
     +  -41.9833633271704727d0,
     +   38830.4051652831113d0,
     +  -12799.7555426009494d0,
     +   1069.18148179751506d0,
     +  -462206.934091776377d0,
     +   145609.693716387614d0,
     +  -4294.7493548165002d0,
     +    1.28174004865438040d0, ! betan(3,3)
     +    9.37167755081878084d0,
     +    5.82920320846327833d0,
     +    2.11674870899066869d0,
     +    1.03598365915185830d0,
     +    4.95061807875526849d0,
     +    6.24058015066356209d0,
     +    3.04097343017986255d0,
     +    3.54381750322741196d0/
c
      logical first
      data first/.true./
      real*8 bohr2a
      data bohr2a/0.529177249d0/
      integer i,j,ii
      real*8 f
      save
c
      write(*,'(3h** ,1a)')'Initializing SAPT5-s potential',
     + 'G. C. Groenenboom, E. M. Mas, R. Bukowski, K. Szalewicz,',
     + 'P. E. S. Wormer, and A. van der Avoird, Phys. Rev. Lett. (2000)'
c
      if (first) then
         first=.false.
c        convert parameters to atomic units
         do i=28,33
            p(i)=p(i)*bohr2a               ! betaq(1..6)
         end do
         do i=34,43
            p(i)=p(i)*bohr2a               ! betar(1..10)
         end do
         ii=43
         do i=1,10
            ii=ii+1
            p(ii)=p(ii)/627.510d0          ! arep(:,1)
         end do
         f=bohr2a
         do j=2,4
            do i=1,10
               ii=ii+1
c turn arep(:,1:4) into linear parameters again
               p(ii)=p(43+i)*p(ii)*f               ! arep(:,2:4)
            end do
            f=f*bohr2a
         end do
         f=627.51d0*bohr2a**6
         do j=1,3
            do i=1,3
               ii=ii+1
               p(ii)=p(ii)/f               ! cn(i,j)
            end do
            f=f*bohr2a*bohr2a
         end do
         do i=1,9
            ii=ii+1
            p(ii)=p(ii)*bohr2a             ! betan
         end do
      end if
      call ss_init(nip,ip,np,p)
      end
c
c***********************************************************************
c
      subroutine sapt5st_init
      implicit none
c
c  Initialize SAPT5-st potential. Use SS_EVAL to evaluate the potential.
c
c  Use SS_EVAL to evaluate the potential.
c
c  Monomer geometry (in x-z plane, O on positive z-axis)
c
c  rOH  =  1.836106337 bohr
c  aHOH =  104.69      degrees
c
c  Sites
c  1:  O
c  2:  H1
c  3:  H2
c  4:  Bunny2 1 ("lone pair repulsion")
c  5:  Bunny2 2
c  6:  M        Center of mass for mass_O=16, mass_H=1
c  7:  Bunny1 1 (with charge)
c  8:  Bunny1 2
c
c     Integer parameters
c
      integer nip
      parameter(nip=12+2*( 5 + 6 + 3 )) ! (12+2*(nq+nr+nl))
      integer ip(nip)
      data ip/
     + 1,          ! Version
     + 8,          ! Number of sites (nsite)
     + 5,          ! Number of sites with charges (nq)   
     + 3,          ! Different types of sites with charges (mq)
     + 6,          ! Number of sites with repulsion (nr)
     + 4,          ! Different types of sites with repulsion (mr)
     + 4,          ! Order of polynomial for repulsive part (mrorder)
     + 3,          ! Number of sites with long range (nl)
     + 2,          ! Different types of sites with long range (ml)
     + 3,6,2,      ! nc,nlow,nstep:   n=nlow+(0:nc-1)*nstep, R^-n

     + 1,2,3,7,8,  ! listq(nq), charges on O, H1, H2, B1_1, B1_2
     + 1,2,2,3,3,  ! itypeq(nq), 1=O, 2=H, 3=B1

     + 1,2,3,6,4,5,! listr(nr), repulsion on O, H1, H2, M, B2_1, B2_1
     + 1,2,2,3,4,4,! ityper(nr), 1=O, 2=H, 3=M, 4=B2

     + 1,2,3,      ! listl(nl), O, H1, H2
     + 1,2,2/      ! itypel(nl), 1=O, 2=H
c
c     Real*8 parameters
c
      integer np
      parameter(np=3*8+3+6+10+10*4+3*3+3*3)
      real*8 p(np)
      data p/
c
c     cartesian coordinates of the sites
c
     +     0.00000000000000d+00,     0.00000000000000d+00,  ! O
     +     1.24631924913843d-01,    -1.45365196228170d+00,  !   H1
     +     0.00000000000000d+00,    -9.97055399310745d-01,
     +     1.45365196228170d+00,     0.00000000000000d+00,  ! H2
     +    -9.97055399310745d-01,     0.00000000000000d+00,  !   B2
     +    -2.00000000000000d-01,     2.50000000000000d-01,
     +     0.00000000000000d+00,     2.00000000000000d-01,  ! B2
     +     2.50000000000000d-01,     0.00000000000000d+00,  !   M
     +     0.00000000000000d+00,     0.00000000000000d+00,
     +     0.00000000000000d+00,    -0.2067213d0         ,  ! B1
     +    -0.247160511d0,            0.00000000000000d+00,  !   B1
     +     0.2067213d0,             -0.247160511d0,
     +     0.258504607653328389d0,   0.564050362569283648d0,! charge
     +    -0.693302666395947953d0,
     +     5.30188399705912605d0, ! betaq
     +     1.40264620725287381d0,
     +     1.37777021367524322d0,
     +     4.88257857151977070d0,
     +     1.45060834365989932d0,
     +     1.6247997153704614d0,
     +     3.96533329909905d0,    ! betar
     +     3.31344984120762d0,
     +     2.41429207550821d0,
     +     3.97675261600551d0,
     +     3.17256531346957d0,
     +     3.98715016325772d0,
     +     3.96341198097565d0,
     +     3.42171736051266d0,
     +     3.96806674259008d0,
     +     3.96112259788312d0,
     +     0.473804466702517D+08, ! arep(10,4)
     +     0.232304456812340D+05,
     +     0.157036053427091D+04,
     +     0.212327780148652D+08,
     +    -0.894572572154939D+04,
     +    -0.335855976729718D+08,
     +    -0.296420327686655D+08,
     +     0.109476792484018D+05,
     +     0.515786593402049D+07,
     +     0.107616042114306D+08,
     +    -0.170254540191501D+00,
     +     0.195092272355747D+01,
     +    -0.363612765499975D+00,
     +    -0.331177830910078D+01,
     +     0.104171475021330D+01,
     +    -0.218302327014551D+01,
     +    -0.103137210054502D+01,
     +    -0.368179002721265D+01,
     +     0.776827754184826D-01,
     +    -0.123761562940866D+01,
     +    -0.107953750104456D+01,
     +    -0.166060145096589D+01,
     +     0.241016384146484D-01,
     +     0.316570351062773D+01,
     +    -0.211412510719483D+00,
     +     0.164634960500608D+01,
     +     0.962834404645738D-01,
     +     0.269715281697085D+01,
     +    -0.140429859686551D+01,
     +     0.364493206178939D+00,
     +     0.561987349572409D+00,
     +     0.154719094083673D+01,
     +     0.516343451845443D-01,
     +    -0.103077437966963D+01,
     +     0.113965727605120D+01,
     +    -0.461360908362144D+00,
     +     0.122550135968013D+00,
     +    -0.134191579566847D+01,
     +     0.696071765603755D+00,
     +     0.249105358806801D-01,
     +    -418.942352042298694d0, ! cn(3,3)
     +     410.613776430728421d0,
     +    -41.9833633271704727d0,
     +     38830.4051652831113d0,
     +    -12799.7555426009494d0,
     +     1069.18148179751506d0,
     +    -462206.934091776377d0,
     +     145609.693716387614d0,
     +    -4294.7493548165002d0,
     +     1.28174004865438040d0, ! betan(3,3)
     +     9.37167755081878084d0,
     +     5.82920320846327833d0,
     +     2.11674870899066869d0,
     +     1.03598365915185830d0,
     +     4.95061807875526849d0,
     +     6.24058015066356209d0,
     +     3.04097343017986255d0,
     +     3.54381750322741196d0/
c
      logical first
      data first/.true./
      real*8 bohr2a
      data bohr2a/0.529177249d0/
      integer i,j,ii
      real*8 f
      save
c
      write(*,'(3h** ,1a)')'Initializing SAPT5-st (tuned) potential',
     + 'G. C. Groenenboom, E. M. Mas, R. Bukowski, K. Szalewicz,',
     + 'P. E. S. Wormer, and A. van der Avoird, Phys. Rev. Lett. (2000)'

      if (first) then
         first=.false.
c        convert parameters to atomic units

         do i=28,33
            p(i)=p(i)*bohr2a               ! betaq(1..6)
         end do
         do i=34,43
            p(i)=p(i)*bohr2a               ! betar(1..10)
         end do
         ii=43
         do i=1,10
            ii=ii+1
            p(ii)=p(ii)/627.510d0          ! arep(:,1)
         end do
         f=bohr2a
         do j=2,4
            do i=1,10
               ii=ii+1
c turn arep(:,1:4) into linear parameters again
               p(ii)=p(43+i)*p(ii)*f               ! arep(:,2:4)
            end do
            f=f*bohr2a
         end do
         f=627.51d0*bohr2a**6
         do j=1,3
            do i=1,3
               ii=ii+1
               p(ii)=p(ii)/f               ! cn(i,j)
            end do
            f=f*bohr2a*bohr2a
         end do
         do i=1,9
            ii=ii+1
            p(ii)=p(ii)*bohr2a             ! betan
         end do
      end if
      call ss_init(nip,ip,np,p)
      end
c
c**********************************************************************
c
      subroutine ss_dshift
      implicit none
c
c     shift for D2O-D2O
c
c     First call SS_INIT(NIP,IP,NP,P) to initialize (H2O)2
c     surface, then call SS_DSHIFT to shift sites for (D2O)2.
c
      real*8 O1zp, H1zp, dist, mH, mD, mO, shift
      parameter(mH=1.007825d0, mD=2.0140d0, mO=15.99491d0)
c
      O1zp=0.1246319249d0
      H1zp=-0.9970553993d0
c
      dist=O1zp-H1zp
      shift=dist*2d0*mO*(mD-mH)/( (2*mH+mO)*(2*mD+mO))
      write(*,*)'SS_DSHIFT, D2O-D2O, shifting by ',shift
      call ss_shift(shift)
      end
c
c***********************************************************************
c
      subroutine ss_eval(coords,etotal)
      implicit none
      real*8 coords(7), etotal
c
c  <coords(7)>
c     input: coords(1..3), zyz euler angles monomer A (rad).
c            coords(4..6), zyz euler angles monomer B (rad).
c            coords(7),    distance (bohr).
c  <etotal>
c     output: total energy in Hartree
c
c  First, call entry point SS_INIT(nip,ip,np,p) (see SAPT5S_INIT
c  or SAPT5ST_INIT)
c
c  Gerrit C. Groenenboom, July 1999.
c
      save
c
c     name convention:
c      q for charge
c      r for repulsion
c      l for long range
c
      integer maxsite, maxq, maxr, maxl, maxc, maxrorder
      parameter(maxsite=8, maxq=3, maxr=4, maxl=3, maxc=3,
     +   maxrorder=4)
c
c     derived dimensions
c
      integer maxq2, maxr2, maxl2
      parameter(maxq2=maxq*(maxq+1)/2,
     +          maxr2=maxr*(maxr+1)/2,
     +          maxl2=maxl*(maxl+1)/2)
c
      real*8 geom(3,maxsite)
      integer itab(maxsite,maxsite), nsite
c
c  Coulomb stuff
c
      integer listq(maxsite), itypeq(maxsite), nq, mq, mq2
      real*8 charge(maxq)
      real*8 betaq(maxq2)
c
c  Repulsion
c
      integer listr(maxsite), nr, mr2, mr
      integer ityper(maxsite)
      real*8  betar(maxr2), arep(maxr2,maxrorder)
c
c  Long range stuff
c
      integer listl(maxsite), itypel(maxsite), nl, ml2, ml, nc, nlow,
     +        nstep,mrorder
      real*8  cn(maxl2,maxc), betan(maxl2,maxc)
c
      real*8 rot(3,3), coor1(3,maxsite), coor2(3,maxsite),
     +       rpt(maxsite,maxsite)
c
      real*8 sum, f_damp, r, x, r2, rn, tmp
      integer i,j,i1,k,ind,n,j1
c
c     entry ss_init(nip,ip,np,p)
c
      integer nip, ip(nip), np
      real*8  p(np)
      integer ii, iversion
      logical init
      data init/.false./
c
c     entry ss_shift(shift)
c
      real*8 shift
c
c     entry ss_print(iout,iprint)
c
      integer iout, iprint

c     rotate monomers
      if (.not.init) then
         write(*,*)'** ERROR in SS_EVAL: first call ss_init'
         stop 300
      endif

      call g_euler(coords(1),coords(2),coords(3),rot)
      do i=1,3
         do k=1,nsite
            sum=0d0
            do j=1,3
               sum=sum+rot(i,j)*geom(j,k)
            end do
            coor1(i,k)=sum
         end do
      end do
      call g_euler(coords(4),coords(5),coords(6),rot)
      do i=1,3
         do k=1,nsite
            sum=0d0
            do j=1,3
               sum=sum+rot(i,j)*geom(j,k)
            end do
            coor2(i,k)=sum
         end do
      end do
      do k=1,nsite
         coor2(3,k)=coor2(3,k)+coords(7)
      end do

c     compute distance matrix

      do i=1,nsite
         do j=1,nsite
            sum=0d0
            do k=1,3
               sum=sum+(coor1(k,i)-coor2(k,j))**2
            end do
            rpt(i,j)=dsqrt(sum)
         end do
      end do

      etotal=0d0

c     coulomb contribution

      do i=1,nq
         i1=itypeq(i)
         sum=0d0
         do j=1,nq
            j1=itypeq(j)
            r=rpt(listq(i),listq(j))
            ind=itab(i1,j1)
            x=betaq(ind)*r
            sum=sum+charge(j1)*(1d0-(1d0+x)*dexp(-x))/r
         end do
         etotal=etotal+charge(i1)*sum
      end do

c     repulsion

      do i=1,nr
         i1=ityper(i)
         do j=1,nr
            r=rpt(listr(i),listr(j))
            ind=itab(i1,ityper(j))
c            etotal=etotal+dexp(-betar(ind)*r)*
c     +      (arep(ind,1)+arep(ind,2)*r+arep(ind,3)*r2+arep(ind,4)*r3)
            tmp=arep(ind,mrorder)
            do k=mrorder-1,1,-1
               tmp=arep(ind,k)+r*tmp
            end do
            etotal=etotal+dexp(-betar(ind)*r)*tmp
         end do
      end do

c     long range

      do i=1,nl
         i1=itypel(i)
         do j=1,nl
            r=rpt(listl(i),listl(j))
            ind=itab(i1,itypel(j))
            r2=r**nstep
            rn=r**nlow
            n=nlow
            do k=1,nc
               etotal=etotal-f_damp(r*betan(ind,k),n)*cn(ind,k)/rn
               n=n+nstep ! n=6,8,10
               rn=rn*r2  ! rn=n**6, n**8, n**10
            end do
         end do
      end do


      return

c......................................................................

      entry ss_init(nip,ip,np,p)
c
c      integer nip, ip(nip), np
c      real*8  p(np)
c
c   <ip(nip)>
c      input: integer parameters
c   <p(np)>
c      input: real*8 parameters (all converted to atomic units)
c      See SAPT5S_INIT and SAPT5ST_INIT.
c
      if (nip.lt.12) then
         write(*,*)'** ERROR in SS_INIT, nip = ',nip,' < 12'
         stop 301
      endif
c
      ii=0
      iversion= ip(1)
      if (iversion.ne.1) then
         write(*,*)'** ERROR: iversion = ',iversion,' must be 1'
         write(*,*)'** ERROR in SS_EVAL entry point SS_INIT'
         stop 302
      end if
      nsite   = ip(2)
      nq      = ip(3)
      mq      = ip(4)
      nr      = ip(5)
      mr      = ip(6)
      mrorder = ip(7)
      nl      = ip(8)
      ml      = ip(9)
      nc      = ip(10)
      nlow    = ip(11)
      nstep   = ip(12)

      mq2     = mq*(mq+1)/2
      mr2     = mr*(mr+1)/2
      ml2     = ml*(ml+1)/2

      if (nsite.gt.maxsite) then
         write(*,*)'** ERROR: nsite=',nsite,'>maxsite=',maxsite
         write(*,*)'** ERROR in SS_EVAL, entry point SS_INIT'
         stop 303
      endif
      if (nq.gt.maxsite) then
         write(*,*)'** ERROR: nq=',nq,'>maxsite=',maxsite
         write(*,*)'** ERROR in SS_EVAL, entry point SS_INIT'
         stop 304
      end if
      if (nr.gt.maxsite) then
         write(*,*)'** ERROR: nr=',nr,'>maxsite=',maxsite
         write(*,*)'** ERROR in SS_EVAL, entry point SS_INIT'
         stop 305
      end if
      if (nl.gt.maxsite) then
         write(*,*)'** ERROR: nl=',nl,'>maxsite=',maxsite
         write(*,*)'** ERROR in SS_EVAL, entry point SS_INIT'
         stop 306
      end if
      if (mq.gt.maxq .or. mr.gt.maxr .or. ml.gt.maxl ) then
        write(*,*)'** mq  =',mq,'mr  =',mr,  'ml  =',ml
        write(*,*)'** maxq=',mq,'maxr=',maxr,'maxl=',maxl
        write(*,*)'** ERROR in SS_EVAL: increase dimensions'
        stop 307
      endif
      if (nc.gt.maxc) then
         write(*,*)'** ERROR: nc=',nc,'>maxc=',maxc
         write(*,*)'** ERROR in SS_EVAL, entry point SS_INIT'
         stop 308
      end if
      if (nip.lt.12+2*(nq+nr+nl)) then
         write(*,*)'** ERROR: nip = ',nip,' < ',12+2*(nq+nr+nl)
         write(*,*)'** ERROR in SS_EVAL, entry point SS_INIT'
         stop 309
      end if

      ii=12
      do i=1,nq
         listq(i)=ip(ii+i)
         itypeq(i)=ip(ii+i+nq)
      end do
      ii=ii+2*nq
      do i=1,nr
         listr(i)=ip(ii+i)
         ityper(i)=ip(ii+i+nr)
      end do
      ii=ii+2*nr
      do i=1,nl
         listl(i)=ip(ii+i)
         itypel(i)=ip(ii+i+nl)
      end do
c
      do i=1,maxsite
         do j=1,i
            ind=i*(i-1)/2+j
            itab(i,j)=ind
            itab(j,i)=ind
         end do
      end do
c
      ii=0
      do j=1,nsite
         do i=1,3
            ii=ii+1
            geom(i,j)=p(ii)
         end do
      end do
      do i=1,mq
         ii=ii+1
         charge(i)=p(ii)
      end do
      do i=1,mq2
         ii=ii+1
         betaq(i)=p(ii)
      end do
      do i=1,mr2
         ii=ii+1
         betar(i)=p(ii)
      end do
      do j=1,mrorder
         do i=1,mr2
            ii=ii+1
            arep(i,j)=p(ii)
         end do
      end do
      do j=1,nc
         do i=1,ml2
            ii=ii+1
            cn(i,j)=p(ii)
         end do
      end do
      do j=1,nc
         do i=1,ml2
            ii=ii+1
            betan(i,j)=p(ii)
         end do
      end do
      init=.true.
      return
c
c......................................................................
c
      entry ss_shift(shift)
c
c     isotope shift

      if (.not.init) then
         write(*,*)'** ERROR: first call SS_INIT'
         write(*,*)'** ERROR in SS_EVAL entry point SS_SHIFT'
         stop 310
      endif
      do j=1,nsite
         geom(3,j)=geom(3,j)+shift
      end do
      return
c
c......................................................................
c
      entry ss_print(iout,iprint)
c
c     print parameters of potential
c
c  <iout>
c     input: unit number of standard output (usually iout=6)
c  <iprint>
c     input: print option: iprint=0 print nothing
c                          iprint=1 only scalar information
c                          iprint>1 also the geometry
c                          iprint>2 also other vectors
c                          iprint>9 also matrices
c
      if (iprint.le.0) return
      write(iout,*)'Site-site potential parameters (in atomic units)'
      write(iout,*)'Geometry'
      write(iout,*)'========'
      write(iout,*)'Number of sites: ',nsite
      if (iprint.gt.1) then
        write(iout,'(1a3,1a14,2a23)')'i','x(i)','y(i)','z(i)'
        do i=1,nsite
          write(iout,'(i3,3(1x,e22.15))')i,geom(1,i),geom(2,i),geom(3,i)
        end do
      end if
      write(iout,*)'Coulomb part'
      write(iout,*)'============'
      write(iout,*)'Number of sites with charges:          nq = ',nq
      write(iout,*)'Different types of sites with charges: mq = ',mq
      if (iprint.gt.2) then
         write(iout,*)'Charges:'
         write(iout,'(1x,2a6,a14)')'site','type','charge'
         do i=1,nq
            write(iout,'(2(i6),2x,e22.15)')listq(i),itypeq(i),
     +          charge(itypeq(i))
         end do
      end if
      if (iprint.gt.9) then
         write(iout,*)'Damping parameters:'
         write(iout,*)'type_a type_b ind     betaq(ind)'
         do i=1,mq
            do j=1,i
               ind=itab(i,j)
               write(iout,'(i6,i6,i5,2x,e22.15)')i,j,ind,betaq(ind)
            end do
         end do
      end if
      write(iout,*)'Repulsion'
      write(iout,*)'========='
      write(iout,*)'Number of sites with repulsion:           nr = ',nr
      write(iout,*)'Different types of sites with repulsion:  mr = ',mr
      if (iprint.gt.2) then
        write(iout,'(1x,2a6)')'site','type'
        do i=1,nr
           write(iout,'(i6,i6))')listr(i),ityper(i)
        end do
      end if
      if (iprint.gt.9) then
        write(iout,*)'arep:'
        write(iout,*)'type_a type_b ind   (arep(ind,k),k=1,mrorder)'
        do i=1,mr
           do j=1,i
              ind=itab(i,j)
              write(iout,'(i6,i6,i5,20(2x,e22.15))')i,j,ind,
     +             (arep(ind,k),k=1,mrorder)
           end do
        end do
        write(iout,*)'Damping parameters:'
        write(iout,*)'type_a type_b ind  betar(ind)'
        do i=1,mr
           do j=1,i
              ind=itab(i,j)
              write(iout,'(i6,i6,i5,2x,e22.15)')i,j,ind,betar(ind)
           end do
        end do
      end if
      write(iout,*)'Long range'
      write(iout,*)'=========='
      write(iout,*)'Number of sites with long range:          nl = ',nl
      write(iout,*)'Different types of sites with long range: ml = ',ml
      if (iprint.gt.2) then
        write(iout,'(1x,2a6)')'site','type'
        do i=1,nl
           write(iout,'(2i6)')listl(i),itypel(i)
        end do
        write(iout,'(2a6)')'k','n'
        do i=1,nc
           write(iout,'(2i6)') i,nlow+(i-1)*nstep
        end do
      end if
      if (iprint.gt.9) then
        write(iout,*)'Long range coefficients:'
        write(iout,*)'type_a type_b ind (cn(ind,k),k=1,nc)'
        do i=1,ml
           do j=1,i
              ind=itab(i,j)
              write(iout,'(i6,i6,i5,20(2x,e22.15))')i,j,ind,
     +         (cn(ind,k),k=1,nc)
           end do
        end do
        write(iout,*)'Damping parameters:'
        write(iout,*)'type_a type_b ind  (betan(ind,k),k=1,nc)'
        do i=1,ml
           do j=1,i
              ind=itab(i,j)
              write(iout,'(i6,i6,i5,20(2x,e22.15))')i,j,ind,
     +          (betan(ind,k),k=1,nc)
           end do
        end do
      end if
      end
c
c**********************************************************************
c
      real*8 function f_damp(br,n)
      implicit none
      real*8 br
      integer n
c
c  Tang-Toennies damping function, careful implementation for small R.
c
      real*8 d, sum, term
      integer i, ncn
c
      if(br.eq.0d0) then
         d=0d0
         f_damp=0d0
         return
      endif
      sum=1d0
      term=1d0
      ncn=n
      do i=1,ncn
         term=term*br/i
         sum=sum+term
      enddo
      d=1d0-dexp(-br)*sum
c     in case of d --> 0 use
c     d=1d0 - dexp(-br)*sum = sum_m=ncn+1^\infty br^m/m!
      if(dabs(d).lt.1d-8) then
         d=0d0
         do i=ncn+1,1000
            term=term*br/i
            d=d+term
            if(term/d.lt.1d-8) go to 111
         enddo
         write(*,*)'** BR = ',BR,' N=',N
         write(*,*)'** ERROR: No convergence in d'
         write(*,*)'** ERROR in F_DAMP'
         stop 400
c
  111    continue
         d=d*dexp(-br)
      endif
      f_damp=d
      end
c
c**********************************************************************
c
      subroutine g_euler(alpha,beta,gamma,R)
      implicit none
      real*8 alpha, beta, gamma, R(3,3)
c
c  Compute Euler rotation matrix in the zyz parameterization:
c  R=R_z(alpha)*R_y(beta)*R_z(gamma)
c
c  <alpha>, <beta>
c     input: polar coordinates of rotation axis
c  <gamma>
c     input: "spin" angle
c  <R(3,3)>
c     output: Euler rotation matrix.
c
      real*8 sa, ca, sb, cb, sc, cc
c
      ca=cos(alpha)
      cb=cos(beta)
      cc=cos(gamma)
      sa=sin(alpha)
      sb=sin(beta)
      sc=sin(gamma)
      R(1,1)=cc*cb*ca-sc*sa
      R(2,1)=cc*cb*sa+sc*ca
      R(3,1)=-cc*sb
      R(1,2)=-sc*cb*ca-cc*sa
      R(2,2)=-sc*cb*sa+cc*ca
      R(3,2)=sc*sb
      R(1,3)=sb*ca
      R(2,3)=sb*sa
      R(3,3)=cb
      end
c
c***********************************************************************
c
