proppot.f

c $Id: proppot.f,v 1.10 1999/01/18 09:57:12 ernst Exp $
c Setting of global paramters
c
      subroutine params
      implicit none
      real*8 A0, chi
      include 'coms.f'

c     gw = 0.25 fm^-2 width of the gaussian

      logSky = .true.
      logYuk = .true.
      logCb  = .true.
      logPau = .false.

      gw     = 0.25
      sgw    = sqrt(gw)
      Cb0    = 1.44
      Yuk0   = 0.0 !-85.0
      gamYuk = 1.4
      drPau  = 9.0
      dpPau  = 0.0144
      Pau0   = 0.0 !99.5*(hqc/sqrt(drPau*dpPau))**3

C
C hard Skyrme EOS (usual stuff)
C
c      Sky30  = 70.5
c      gamSky = 2.0
c      A0     = -124.2 * 0.5
C
C hard Skyrme EOS (JK parametrisation corrected for Gausswidth gw=0.25)
C
       Sky30  = 125.93
       gamSky = 1.676
       A0     = -87.67
       chi    = 0.93

c      Sky30  = 303.0 !188.18
c      gamSky = 7.0/6.0 !1.457
c      A0 = -356.0 * 0.5


      Sky20 = chi*2.0*A0
      Yuk0 = (1.0-chi)/(2.0*pi*gamYuk**2)*A0

c      Sky20 = 0.0d0
c      Sky30 = 0.0d0
c      Yuk0  = 0.0d0
c      Cb0   = 0.0d0

      delr = 0.2
      fdel = delr*delr/6.0
      da = -(1.0/delr)
      db = -da

      cutdww = 20.0
      cutPau = 20.0
      cutYuk = 20.0
      cutCb  = 20.0

      dtimestep=0.2
      dt  = 0.02
c      dt2 = 0.5*dt
c      dt6 = dt/6.0

      return
      end


c Reset of all indexed variables
c
      subroutine set0
      implicit none
      integer i, j
      include 'coms.f'

      do 10 i=1,nspl
        spPauy(i) = 0.0
        outPau(i) = 0.0
        spCby(i)  = 0.0
        outCb(i)  = 0.0
        spYuky(i) = 0.0
        outYuk(i) = 0.0
        spSkyy(i) = 0.0
        outSky(i) = 0.0
        spdwwy(i) = 0.0
        outdww(i) = 0.0
  10  continue

      do 20 j=1,nmax
        spin(j)  = 0
        iso3(j)   = 0
        ncoll(j) = 0
        rx(j)    = 0.0
        ry(j)    = 0.0
        rz(j)    = 0.0
        px(j)    = 0.0
        py(j)    = 0.0
        pz(j)    = 0.0
        fmass(j) = 0.0
  20  continue
      return
      end



      subroutine derivs(row)
      implicit none
      integer j, k, index, row
      real*8 spu, spo, outu, outo, tmp, a, b, dy, dp, drdp, dpj
      real*8 rxjku, ryjku, rzjku, rjku, pxjku, pyjku, pzjku, pjku
      logical iPau
      include 'coms.f'

C aopx(j,?) = -dH/dx_j
C aorx(j,?) = dH/dpx_j

      do 10 j=1,nbar
        rww(j) = 0.0
   10 continue
      if(logSky) then
      do 20 j=1,nbar
        do 30 k=j+1,nbar
          rxjku = (airx(j)-airx(k))
          ryjku = (airy(j)-airy(k))
          rzjku = (airz(j)-airz(k))
          rjku = sqrt(rxjku**2+ryjku**2+rzjku**2)
          if(rjku.lt.cutdww) then
            index = int(rjku/delr)+1
            a = dble(index) - rjku/delr
            b = 1.0 - a
            tmp = a*spdwwy(index)+b*spdwwy(index+1)+
     +       ((a**3-a)*outdww(index)+(b**3-b)*outdww(index+1))*fdel
            rww(j) = rww(j) + tmp
            rww(k) = rww(k) + tmp
          end if
  30    continue
  20  continue
      end if

      do 40 j=1,nbar
        aopx(j,row) = 0.0
        aopy(j,row) = 0.0
        aopz(j,row) = 0.0
        dpj = 1.0/sqrt(aipx(j)*aipx(j)+aipy(j)*aipy(j)+aipz(j)*aipz(j)+
     +             fmass(j)*fmass(j))
        aorx(j,row) = aipx(j)*dpj
        aory(j,row) = aipy(j)*dpj
        aorz(j,row) = aipz(j)*dpj
  40  continue

      do 50 j=1,nbar
        do 60 k=j+1,nbar
          rxjku = (airx(j)-airx(k))
          ryjku = (airy(j)-airy(k))
          rzjku = (airz(j)-airz(k))
          rjku = sqrt(rxjku**2+ryjku**2+rzjku**2)
          if (rjku.ge.1.0E-8) then
            rxjku = rxjku/rjku
            ryjku = ryjku/rjku
            rzjku = rzjku/rjku
          else
            rxjku = 0.0
            ryjku = 0.0
            rzjku = 0.0
          end if
          spu  = 0.0
          spo  = 0.0
          outu = 0.0
          outo = 0.0
          dy = 0.0
          index = int(rjku/delr)+1
          a = dble(index)-rjku/delr
          b = 1.0-a
          if(logYuk.and.rjku.lt.cutYuk) then
            spu  = spu  + spYuky(index)
            spo  = spo  + spYuky(index+1)
            outu = outu + outYuk(index)
            outo = outo + outYuk(index+1)
          end if
          if(logSky.and.rjku.lt.cutdww) then
            tmp = Sky20 + Sky30*gamSky/(gamSky+1.0)*
     *            (rww(j)**(gamSky-1.0)+rww(k)**(gamSky-1.0))
            spu  = spu  + spdwwy(index)*tmp
            spo  = spo  + spdwwy(index+1)*tmp
            outu = outu + outdww(index)*tmp
            outo = outo + outdww(index+1)*tmp
          end if
          dy = da*spu+db*spo+
     +      ((3.0*a**2-1.0)*da*outu+(3.0*b**2-1.0)*db*outo)*fdel
          if(logCb) then
          if(rjku.lt.cutCb) then
            dy = dy + (da*spCby(index)+db*spCby(index+1)+
     +                ((3.0*a**2-1.0)*da*outCb(index)+
     +                (3.0*b**2-1.0)*db*outCb(index+1))*fdel)*
     *                dble(charge(j)*charge(k))
          else
            dy = dy - Cb0/rjku/rjku*dble(charge(j)*charge(k))
          end if
          end if
          if(logPau.and.iPau(j,k)) then
            pxjku = (aipx(j)-aipx(k))
            pyjku = (aipy(j)-aipy(k))
            pzjku = (aipz(j)-aipz(k))
            pjku = sqrt(pxjku**2+pyjku**2+pzjku**2)
            if (pjku.ge.1.0E-8) then
              pxjku = pxjku/pjku
              pyjku = pyjku/pjku
              pzjku = pzjku/pjku
            else
              pxjku = 0.0
              pyjku = 0.0
              pzjku = 0.0
            end if
            drdp = 0.5*(pjku*pjku/dpPau+rjku*rjku/drPau)
            if(drdp.lt.cutPau) then
              index = int(drdp/delr)+1
              a = dble(index)-drdp/delr
              b = 1.0-a
              tmp = da*spPauy(index)+db*spPauy(index+1)+
     +              ((3.0*a**2-1.0)*da*outPau(index)+
     +               (3.0*b**2-1.0)*db*outPau(index+1))*fdel
              dy = dy+tmp*rjku/drPau
              dp = tmp*pjku/dpPau*0.001
              aorx(j,row) = aorx(j,row)+dp*pxjku
              aory(j,row) = aory(j,row)+dp*pyjku
              aorz(j,row) = aorz(j,row)+dp*pzjku
              aorx(k,row) = aorx(k,row)-dp*pxjku
              aory(k,row) = aory(k,row)-dp*pyjku
              aorz(k,row) = aorz(k,row)-dp*pzjku
            end if
          end if
          dy = -(0.001*dy)
          aopx(j,row) = aopx(j,row)+dy*rxjku
          aopy(j,row) = aopy(j,row)+dy*ryjku
          aopz(j,row) = aopz(j,row)+dy*rzjku
          aopx(k,row) = aopx(k,row)-dy*rxjku
          aopy(k,row) = aopy(k,row)-dy*ryjku
          aopz(k,row) = aopz(k,row)-dy*rzjku
  60    continue
  50  continue
      return
      end

      real*8 function Ekintot()
      implicit none
      integer j
      real*8 Ekin
      include 'coms.f'

      Ekintot = 0.0
      do 3 j=1,npart
         Ekintot= Ekintot+Ekin(j)
 3    continue
      return
      end

      real*8 function EtotJK()
      implicit none
      real*8 Etot
      integer j
      include 'coms.f'

      EtotJK = Etot()
      do 3 j=1,npart
         EtotJK= EtotJK-fmass(j)
 3    continue
      EtotJK = EtotJK/npart
      return
      end


      real*8 function Etot()
      implicit none
      integer j, k, index
      real*8 a, b, y, drdp, tp, tr, tmp, Ekintot
      real*8 Ekinbar, Ekinmes, ESky2, ESky3, EYuk, ECb, EPau, Ekin
      real*8 rxjku, ryjku, rzjku, rjku, pxjku, pyjku, pzjku, pjku
      logical iPau
      include 'coms.f'
      common /energies/ Ekinbar, Ekinmes, ESky2, ESky3, EYuk, ECb, EPau

      Etot = 0.0
      Ekinbar = 0.0
      Ekinmes = 0.0
      ESky2 = 0.0
      ESky3 = 0.0
      EYuk = 0.0
      ECb = 0.0
      EPau = 0.0

      if(EoS.eq.0) then

c CASCADE mode
         Etot=Ekintot()
         return
      else
c with potentials
c kinetic energies of mesons first
         do 4 j=nbar+1,npart
            Etot= Etot+Ekin(j)
            Ekinmes = Ekinmes+Ekin(j)
 4       continue

      do 10 j=1,nbar
        rww(j) = 0.0
   10 continue
      if(logSky) then
      do 20 j=1,nbar
        do 30 k=j+1,nbar
          rxjku = (rx(j)-rx(k))
          ryjku = (ry(j)-ry(k))
          rzjku = (rz(j)-rz(k))
          rjku = sqrt(rxjku**2+ryjku**2+rzjku**2)
          if(rjku.lt.cutdww) then
            index = int(rjku/delr)+1
            a = dble(index) - rjku/delr
            b = 1.0 - a
            tmp = a*spdwwy(index)+b*spdwwy(index+1)+
     +       ((a**3-a)*outdww(index)+(b**3-b)*outdww(index+1))*fdel
            rww(j) = rww(j) + tmp
            rww(k) = rww(k) + tmp
          end if
  30    continue
  20  continue
      end if

      do 40 j=1,nbar
        Etot = Etot + Ekin(j) + 0.0005*Sky20*rww(j) +
     +         0.001*Sky30/(gamSky+1.0)*rww(j)**gamSky
        Ekinbar = Ekinbar + Ekin(j)
        ESky2 = ESky2 + 0.0005*Sky20*rww(j)
        ESky3 = ESky3 + 0.001*Sky30/(gamSky+1.0)*rww(j)**gamSky
        do 50 k=j+1,nbar
          rxjku = (rx(j)-rx(k))
          ryjku = (ry(j)-ry(k))
          rzjku = (rz(j)-rz(k))
          rjku = sqrt(rxjku**2+ryjku**2+rzjku**2)
          index = int(rjku/delr)+1
          a = dble(index)-rjku/delr
          b = 1.0-a
          if(logYuk.and.rjku.lt.cutYuk) then
            y = a*spYuky(index)+b*spYuky(index+1)+
     +          ((a**3-a)*outYuk(index)+(b**3-b)*outYuk(index+1))*fdel
            Etot = Etot + 0.001*y
            EYuk = EYuk + 0.001*y
          end if
          if(logCb) then
          if(rjku.lt.cutCb) then
            y = (a*spCby(index)+b*spCby(index+1)+
     +          ((a**3-a)*outCb(index)+(b**3-b)*outCb(index+1))*fdel)*
     *          dble(charge(j)*charge(k))
          else
            y = Cb0/rjku*dble(charge(j)*charge(k))
          end if
          Etot = Etot + 0.001*y
          ECb = ECb + 0.001*y
          end if
          if(logPau.and.iPau(j,k)) then
            pxjku = (px(j)-px(k))
            pyjku = (py(j)-py(k))
            pzjku = (pz(j)-pz(k))
            pjku = sqrt(pxjku**2+pyjku**2+pzjku**2)
            tp = pjku
            tr = rjku
            drdp = 0.5*(pjku*pjku/dpPau+rjku*rjku/drPau)
            if(drdp.lt.cutPau) then
              index = int(drdp/delr)+1
              a = dble(index)-drdp/delr
              b = 1.0-a
              y = a*spPauy(index)+b*spPauy(index+1)+
     +            ((a**3-a)*outPau(index)+(b**3-b)*outPau(index+1))*fdel
              Etot = Etot + 0.001*y
              EPau = EPau + 0.001*y
            end if
          end if
  50    continue
  40  continue
      Ekinbar = Ekinbar/dble(nbar)
      Ekinmes = Ekinmes/dble(max(1,npart-nbar))
      ESky2 = ESky2/dble(nbar)
      ESky3 = ESky3/dble(nbar)
      EYuk = EYuk/dble(nbar)
      ECb = ECb/dble(nbar)
      EPau = EPau/dble(nbar)
      end if
      return
      end


      subroutine cascstep(tim,dtime)
      implicit none
      real*8 tim,dtime,energ
      integer j
      include 'coms.f'
      include 'boxinc.f'
      include 'options.f'
ctp060202 to avoid warnings with gfortran compilation
      logical ctp060202
      ctp060202=.false.
      if(ctp060202)write(*,*)tim
ctp060202 end

      do 1 j=1,npart
         energ = sqrt(px(j)**2+py(j)**2+pz(j)**2+fmass(j)**2)
         r0(j) = r0(j) + dtime
         rx(j) = rx(j) + px(j)/energ*dtime
         rz(j) = rz(j) + pz(j)/energ*dtime
         ry(j) = ry(j) + py(j)/energ*dtime
1     continue
      return
      end

      subroutine proprk(tim,dtime)

      implicit none
      real*8 tim,dtime,energ,dt2, dt6
      integer j

      include 'coms.f'
      include 'boxinc.f'
      include 'options.f'
ctp060202 to avoid warnings with gfortran compilation
      logical ctp060202
      ctp060202=.false.
      if(ctp060202)write(*,*)tim
ctp060202 end

      if (EoS.eq.0) then
c  cascade mode
         do 1 j=1,npart
               energ = p0(j)    ! sqrt(px(j)**2+py(j)**2+pz(j)**2+fmass(j)**2)
               r0(j) = r0(j) + dtime
               rx(j) = rx(j) + px(j)/energ*dtime
               ry(j) = ry(j) + py(j)/energ*dtime
               rz(j) = rz(j) + pz(j)/energ*dtime
1        continue
         return
      else
c propagation with potentials
c propagate mesons on straight lines
         do 2 j=nbar+1,npart
            energ = p0(j) ! sqrt(px(j)**2+py(j)**2+pz(j)**2+fmass(j)**2)
            r0(j) = r0(j) + dtime
            rx(j) = rx(j) + px(j)/energ*dtime
            ry(j) = ry(j) + py(j)/energ*dtime
            rz(j) = rz(j) + pz(j)/energ*dtime
2        continue

c propagate baryons
c adjust time-step parameters
      dt = dtime
      dt2 = dtime/2.0d0
      dt6 = dtime/6.0d0

      do 10 j=1,nbar
        airx(j) = rx(j)
        airy(j) = ry(j)
        airz(j) = rz(j)
        aipx(j) = px(j)
        aipy(j) = py(j)
        aipz(j) = pz(j)
  10  continue

      call derivs(1)
      do 20 j=1,nbar
        airx(j) = rx(j) + dt2*aorx(j,1)
        airy(j) = ry(j) + dt2*aory(j,1)
        airz(j) = rz(j) + dt2*aorz(j,1)
        aipx(j) = px(j) + dt2*aopx(j,1)
        aipy(j) = py(j) + dt2*aopy(j,1)
        aipz(j) = pz(j) + dt2*aopz(j,1)
  20  continue

      call derivs(2)

      do 30 j=1,nbar
        airx(j) = rx(j) + dt2*aorx(j,2)
        airy(j) = ry(j) + dt2*aory(j,2)
        airz(j) = rz(j) + dt2*aorz(j,2)
        aipx(j) = px(j) + dt2*aopx(j,2)
        aipy(j) = py(j) + dt2*aopy(j,2)
        aipz(j) = pz(j) + dt2*aopz(j,2)
  30  continue

      call derivs(3)

      do 40 j=1,nbar
        airx(j) = rx(j) + dt*aorx(j,3)
        airy(j) = ry(j) + dt*aory(j,3)
        airz(j) = rz(j) + dt*aorz(j,3)
        aipx(j) = px(j) + dt*aopx(j,3)
        aipy(j) = py(j) + dt*aopy(j,3)
        aipz(j) = pz(j) + dt*aopz(j,3)
  40  continue

      call derivs(4)

      do 50 j=1,nbar
         r0(j) = r0(j) + dtime
        rx(j)=rx(j)+dt6*(aorx(j,1)+2.0*(aorx(j,2)+aorx(j,3))+aorx(j,4))
        ry(j)=ry(j)+dt6*(aory(j,1)+2.0*(aory(j,2)+aory(j,3))+aory(j,4))
        rz(j)=rz(j)+dt6*(aorz(j,1)+2.0*(aorz(j,2)+aorz(j,3))+aorz(j,4))
        px(j)=px(j)+dt6*(aopx(j,1)+2.0*(aopx(j,2)+aopx(j,3))+aopx(j,4))
        py(j)=py(j)+dt6*(aopy(j,1)+2.0*(aopy(j,2)+aopy(j,3))+aopy(j,4))
        pz(j)=pz(j)+dt6*(aopz(j,1)+2.0*(aopz(j,2)+aopz(j,3))+aopz(j,4))
        p0(j)=sqrt(px(j)**2+py(j)**2+pz(j)**2+fmass(j)**2)
 50   continue
      end if

      return
      end


      subroutine potPau
      implicit none
      integer i, ncut, index
      real*8 Ecut, dr, abl0, abln, a, b, y, dy, Pau
      include 'coms.f'

      rx(1) = 0.0d0
      ry(1) = 0.0d0
      rz(1) = 0.0d0
      ry(2) = 0.0d0
      rz(2) = 0.0d0
      px(1) = 0.0d0
      py(1) = 0.0d0
      pz(1) = 0.0d0
      px(1) = 0.0d0
      py(2) = 0.0d0
      pz(2) = 0.0d0
      Ecut = 1.0E-5
      i = 0
  99  i = i+1
      dr = delr*dble(i-1)
      rx(2) = sqrt(2.0*dr*drPau)
      spx(i) = dr
      spPauy(i) = Pau(1,2)
      if(spPauy(i).lt.Ecut) then
        spPauy(i) = 0.0
        cutPau = dr
        abl0 = -Pau0
        abln = 0.0
        ncut = i
      else
        goto 99
      end if
      call spline(spx,spPauy,ncut,abl0,abln,outPau)

      write(6,'(''Pauli-Potential    '',e10.3,i5,f7.1)')
     +      Ecut, ncut, cutPau

      do 10 i=0,20
        dr = 0.323*dble(i)
        if(dr.gt.cutPau) then
          y = 0.0
          dy = 0.0
        else
        rx(2) = dr
        dr = 0.5*dr*dr/drPau
        index = int(dr/delr)+1
        a = dble(index) - dr/delr
        b = 1.0 - a
        y = a*spPauy(index)+b*spPauy(index+1)+
     +     ((a**3-a)*outPau(index)+
     +      (b**3-b)*outPau(index+1))*fdel
        dy = da*spPauy(index)+db*spPauy(index+1)+
     +        ((3.0*a**2-1.0)*da*outPau(index)+
     +         (3.0*b**2-1.0)*db*outPau(index+1))*fdel
        dy = dy*sqrt(2.0*dr*drPau)/drPau
        end if
  10  continue
      return
      end

      subroutine potCb
      implicit none
      integer i, ncut, index
      real*8 Ecut, dr, abl0, abln, a, b, y, dy, dCb, Cb
      include 'coms.f'

      rx(1) = 0.0d0
      ry(1) = 0.0d0
      rz(1) = 0.0d0
      ry(2) = 0.0d0
      rz(2) = 0.0d0
      Ecut = 1.0E-5
      iso3(1) = 1
      iso3(2) = 1
      i = 0
  99  i = i+1
      dr = delr*dble(i-1)
      rx(2) = dr
      spx(i) = dr
      spCby(i) = Cb(1,2)
      if(abs(spCby(i)*dr-Cb0)/max(dr,1.0d-5).lt.Ecut) then
        spCby(i) = Cb0/dr
        cutCb  = dr
        abln = dCb(1,2)
        abl0 = 0.0
        ncut = i
      else
        goto 99
      end if
      call spline(spx,spCby,ncut,abl0,abln,outCb)

      write(6,'(''Coulomb-Potential  '',e10.3,i5,f7.1)')
     +      Ecut, ncut, cutCb

      do 10 i=0,20
        dr = 0.2*dble(i)+0.01212
        rx(2) = dr
        if(dr.ge.cutCb) then
          y = Cb0/dr
          dy = -(Cb0/dr/dr)
        else
          index = int(dr/delr)+1
          a = dble(index) - dr/delr
          b = 1.0 - a
          y = a*spCby(index)+b*spCby(index+1)+
     +       ((a**3-a)*outCb(index)+(b**3-b)*outCb(index+1))*fdel
          dy = da*spCby(index)+db*spCby(index+1)+
     +         ((3.0*a**2-1.0)*da*outCb(index)+
     +           (3.0*b**2-1.0)*db*outCb(index+1))*fdel
        end if
  10  continue
      return
      end

      subroutine potYuk
      implicit none
      integer i, ncut, index
      real*8 Ecut, dr, abl0, abln, a, b, y, dy
      real*8 Yuk
      include 'coms.f'

      rx(1) = 0.0d0
      ry(1) = 0.0d0
      rz(1) = 0.0d0
      ry(2) = 0.0d0
      rz(2) = 0.0d0
      Ecut = 1.0E-5
      i = 0
  99  i = i+1
      dr = delr*dble(i-1)
      rx(2) = dr
      spx(i) = dr
      spYuky(i) = Yuk(1,2)
      if(abs(spYuky(i)).lt.Ecut) then
        spYuky(i) = 0.0
        cutYuk = dr
        abl0 = 0.0
        abln = 0.0
        ncut = i
      else
        goto 99
      end if
      call spline(spx,spYuky,ncut,abl0,abln,outYuk)

      write(6,'(''Yukawa-Potential   '',e10.3,i5,f7.1)')
     +      Ecut, ncut, cutYuk

      do 10 i=0,40
        dr = 0.2*dble(i)
        rx(2) = dr
        if(dr.gt.cutYuk) then
          y = 0.0
          dy = 0.0
        else
        index = int(dr/delr)+1
        a = dble(index) - dr/delr
        b = 1.0 - a
        y = a*spYuky(index)+b*spYuky(index+1)+
     +     ((a**3-a)*outYuk(index)+(b**3-b)*outYuk(index+1))*fdel
        dy = da*spYuky(index)+db*spYuky(index+1)+
     +        ((3.0*a**2-1.0)*da*outYuk(index)+
     +         (3.0*b**2-1.0)*db*outYuk(index+1))*fdel
        end if
  10  continue
      return
      end


      subroutine potdww
      implicit none
      integer i, ncut, index
      real*8 Ecut, dr, abl0, abln, a, b, y, dy
      real*8 dww
      include 'coms.f'

      rx(1) = 0.0d0
      ry(1) = 0.0d0
      rz(1) = 0.0d0
      ry(2) = 0.0d0
      rz(2) = 0.0d0
      Ecut = 1.0E-8
      i = 0
  99  i = i+1
      dr = delr*dble(i-1)
      rx(2) = dr
      spx(i) = dr
      spdwwy(i) = dww(1,2)
      if(abs(spdwwy(i)).lt.Ecut) then
        spdwwy(i) = 0.0
        cutdww = dr
        abl0 = 0.0
        abln = 0.0
        ncut = i
      else
        goto 99
      end if
      call spline(spx,spdwwy,ncut,abl0,abln,outdww)

      write(6,'(''Interaction-Density'',e10.3,i5,f7.1)')
     +      Ecut, ncut, cutdww

      do 10 i=0,20
        dr = 0.295*dble(i)
        rx(2) = dr
        if(dr.gt.cutdww) then
          y = 0.0
          dy = 0.0
        else
        index = int(dr/delr)+1
        a = dble(index) - dr/delr
        b = 1.0 - a
        y = a*spdwwy(index)+b*spdwwy(index+1)+
     +     ((a**3-a)*outdww(index)+(b**3-b)*outdww(index+1))*fdel
        dy = da*spdwwy(index)+db*spdwwy(index+1)+
     +        ((3.0*a**2-1.0)*da*outdww(index)+
     +         (3.0*b**2-1.0)*db*outdww(index+1))*fdel
        end if
  10  continue
      return
      end

c Kinetic Energy
c
      function Ekin(j)
      implicit none
      integer j
      real*8 Ekin
      include 'coms.f'

      Ekin = sqrt((px(j)+ffermpx(j))*(px(j)+ffermpx(j))+
     +            (py(j)+ffermpy(j))*(py(j)+ffermpy(j))+
     +            (pz(j)+ffermpz(j))*(pz(j)+ffermpz(j))+
     +            fmass(j)*fmass(j))

      return
      end

c Derivative for Kinetic Energy
c
      function dEkin(j)
      implicit none
      integer j
      real*8 dEkin
      include 'coms.f'

      dEkin = 1.0/sqrt(px(j)*px(j)+py(j)*py(j)+pz(j)*pz(j)+
     +                 fmass(j)*fmass(j))
      return
      end

c Skyrme Potential (3-body) rwwterm
c
      function dww(j,k)
      implicit none
      integer j, k
      real*8 dww, rjk
      include 'coms.f'

      dww = gw/pi*sqrt(gw/pi)*exp(-(gw*rjk(j,k)*rjk(j,k)))/
     /      rho0
      return
      end

cc Skyrme Potential
cc
      function Sky(j,k)
      implicit none
      integer j, k
      real*8 Sky, rjk
      include 'coms.f'

      Sky = Sky20*gw/pi*sqrt(gw/pi)*exp(-(gw*rjk(j,k)*rjk(j,k)))
      return
      end

c Coulomb Potential
c
      function Cb(j,k)
      implicit none
      integer j, k
      real*8 Cb, rjk
      real*8 erf
      include 'coms.f'

      if (iso3(j).eq.1.and.iso3(k).eq.1) then
        if (rjk(j,k).lt.eps) then
          Cb = Cb0*er0*sgw
        else
          Cb = Cb0/rjk(j,k)*erf(sgw*rjk(j,k))
        end if
      else
        Cb = 0.0
      end if
      return
      end

c Derivative for Coulomb Potential
c
      function dCb(j,k)
      implicit none
      integer j, k
      real*8 dCb, rjk
      real*8 erf
      include 'coms.f'

      if (iso3(j).eq.1.and.iso3(k).eq.1) then
        if (rjk(j,k).lt.eps) then
          dCb = 0.0
        else
          dCb = Cb0*(er0*exp(-(gw*rjk(j,k)*rjk(j,k)))*sgw*rjk(j,k)-
     +               erf(sgw*rjk(j,k)))/rjk(j,k)/rjk(j,k)
        end if
      else
        dCb = 0.0
      end if
      return
      end

c Yukawa Potential
c
      function Yuk(j,k)
      implicit none
      integer j, k
      real*8 Yuk, rjk
      real*8 erf
      include 'coms.f'

      if(rjk(j,k).lt.eps) then
        Yuk = Yuk0*(er0*sgw-exp(0.25/gamYuk/gamYuk/gw)/gamYuk*
     *              (1.0-erf(0.5/gamYuk/sgw)))
      else
        Yuk = Yuk0*0.5/rjk(j,k)*exp(0.25/gamYuk/gamYuk/gw)*
     *           (exp(-(rjk(j,k)/gamYuk))*
     +            (1.0-erf(0.5/gamYuk/sgw-sgw*rjk(j,k)))-
     -            exp(rjk(j,k)/gamYuk)*
     +            (1.0-erf(0.5/gamYuk/sgw+sgw*rjk(j,k))))
      end if
      return
      end

c Derivative for Yukawa Potential
c
      function dYuk(j,k)
      implicit none
      integer j, k
      real*8 dYuk, rjk
      real*8 erf
      include 'coms.f'

      if(rjk(j,k).lt.eps) then
        dYuk = 0.0
      else
        dYuk = 0.5*Yuk0/rjk(j,k)*( exp(0.25/gamYuk/gamYuk/gw)*(
     *          (-(1.0/rjk(j,k))-1.0/gamYuk)*exp(-(rjk(j,k)/gamYuk))*
     *             (1.0-erf(0.5/gamYuk/sgw-sgw*rjk(j,k))) +
     *          (1.0/rjk(j,k)-1.0/gamYuk)*exp(rjk(j,k)/gamYuk)*
     *             (1.0-erf(0.5/gamYuk/sgw+sgw*rjk(j,k))) ) +
     +          sgw*er0*2.0*exp(-(gw*rjk(j,k)*rjk(j,k))) )
      end if
      return
      end


c Pauli Potential
c
      function Pau(j,k)
      implicit none
      integer j, k
      real*8 Pau, pjk, rjk
      include 'coms.f'

      Pau = Pau0*exp(-(0.5*rjk(j,k)*rjk(j,k)/drPau))*
     *           exp(-(0.5*pjk(j,k)*pjk(j,k)/dpPau))
      return
      end

c Derivative (p) for Pauli Potential
c
      function dPaup(j,k)
      implicit none
      integer j, k
      real*8 dPaup, pjk, rjk
      include 'coms.f'

      dPaup = -(Pau0/dpPau*pjk(j,k)*
     *                    exp(-(0.5*rjk(j,k)*rjk(j,k)/drPau))*
     *                    exp(-(0.5*pjk(j,k)*pjk(j,k)/dpPau)))
      return
      end

c Derivative (r) for Pauli Potential
c
      function dPaur(j,k)
      implicit none
      integer j, k
      real*8 dPaur, pjk, rjk
      include 'coms.f'

      dPaur = -(Pau0/drPau*rjk(j,k)*
     *                   exp(-(0.5*rjk(j,k)*rjk(j,k)/drPau))*
     *                   exp(-(0.5*pjk(j,k)*pjk(j,k)/dpPau)))
      return
      end

      function rjk(j,k)
      implicit none
      integer j, k
      real*8 rjk
      include 'coms.f'

      rjk = sqrt((rx(j)-rx(k))**2+(ry(j)-ry(k))**2+(rz(j)-rz(k))**2)
      return
      end

      function pjk(j,k)
      implicit none
      integer j, k
      real*8 pjk
      include 'coms.f'

      pjk = sqrt((px(j)-px(k))**2+(py(j)-py(k))**2+(pz(j)-pz(k))**2)
      return
      end

      function iPau(j,k)
      implicit none
      integer j, k
      logical iPau
      include 'coms.f'

      iPau = .false.
      if (iso3(j).eq.iso3(k).and.ityp(j).eq.ityp(k)) iPau = .true.
      return
      end