      MODULE argonne
       USE precision
       USE gauss
       USE spline

       IMPLICIT NONE
       PRIVATE
       PUBLIC potav18lsj

C     printlevel 
       INTEGER,PARAMETER :: printlevel=0,igreen=0,ispin=1

C     variables to keep track of last parameters used in calls 
       LOGICAL :: first=.true.     ! first ever called 
       INTEGER :: oldl=-1,olds=-1,oldj=-1,oldchan=-1
       INTEGER :: imodel=-1,empotflag=-1

C     r-space grids used 
C     to be defined by input file 
       INTEGER :: nr1pot=20,nr2pot=10,nrpot=30
       REAL(dpreal) :: r1pot=1.0,r2pot=5.0,r3pot=20.0

       INTEGER :: nr1bess=200,nr2bess=100,nrbess=300
       REAL(dpreal) :: r1bess=1.0,r2bess=5.0,r3bess=20.0

C     arrays for r-space grids for potential and besselfunction 

       REAL(dpreal),ALLOCATABLE :: rpot(:),wrpot(:) 
       REAL(dpreal),ALLOCATABLE :: rbess(:),wrbess(:)

C     arrays for splines from rpot to rbess 
       REAL(dpreal),ALLOCATABLE :: rspl(:,:)

C     array to store the operator valued r-space potentials 
       REAL(dpreal),ALLOCATABLE :: v_c(:,:,:),v_t(:,:,:),v_so(:,:,:),
     $                             v_s(:,:,:),v_qso(:,:,:),v_l2(:,:,:),
     $                             v_l2s(:,:,:),v_ls2(:,:,:)
       REAL(dpreal),ALLOCATABLE :: fi_c(:,:,:),fi_t(:,:,:),fi_so(:,:,:),
     $                             fi_s(:,:,:),fi_qso(:,:,:)
       
C     arrays to store the current partial wave projected rspace potential 
       REAL(dpreal),ALLOCATABLE :: v_lsj(:,:,:) 
       REAL(dpreal),ALLOCATABLE :: fi_lsj(:,:,:) 

       REAL(dpreal) :: massfakt   ! to be set in initav18 and initchan 
       REAL(dpreal) :: hbarc=197.327053 ! fm to MeV 

C     arrays to store bessel functions 
       REAL(dpreal),ALLOCATABLE :: besf(:,:),besfp(:,:),besspot(:)

       CONTAINS
C     routine for the momentum space potential 
C     this version is preliminary and only prepares part of the potential
C     pp,p outgoing and incoming momentum in fm-1 
C     lfix,sfix,jfix  orbital, spin and total angular momentum of the NN system
C     V00,V20,V02,V22, coupled channel matrix elements of V in order 
C        V(1,1),V(2,1),V(1,2) and V(2,2) 
C     iforce = 1,2,3 = pp,np and nn system

       SUBROUTINE potav18lsj(pp,p,lfix,sfix,jfix,V00,V02,V20,V22,koup,
     $                        iforce) 
        IMPLICIT NONE 
        INTEGER  :: lfix,sfix,jfix,iforce
        REAL(dpreal) :: pp,p,V00,V20,V02,V22
        INTEGER irbess,irpot
        REAL(dpreal) :: momfakt
        LOGICAL koup

        IF(first) THEN
         CALL initav18
         first=.false.
        END IF

        IF(oldl.NE.lfix.OR.olds.NE.sfix
     $    .OR.oldj.NE.jfix.OR.oldchan.NE.iforce) THEN
         oldl=lfix
         olds=sfix
         oldj=jfix
         oldchan=iforce
         IF(oldl.NE.oldj) THEN  ! force l on lowest value for the calculation
          oldl=abs(oldj-olds)
         END IF
         CALL initav18chan
        END IF

        IF(oldj.NE.0.AND.oldl.NE.oldj) THEN ! coupled channels ?

C     besfp
         CALL SBF1(oldl,pp,rbess,nrbess,besfp(:,1))
         CALL SBF1(oldl+2,pp,rbess,nrbess,besfp(:,2))

         IF(printlevel.GT.2) THEN
          WRITE(*,*) 'Bessel function  for pp'
          WRITE(*,'(A15,3X,A15)')
     $      'R in fm','BESFP'
          DO irbess=1,nrbess
           WRITE(*,'(I6,E15.6,3X,6E15.6)')
     $          irbess,rbess(irbess),besfp(irbess,1),besfp(irbess,2)
          END DO
         END IF

         CALL SBF1(oldl,p,rbess,nrbess,besf(:,1))
         CALL SBF1(oldl+2,p,rbess,nrbess,besf(:,2))

         IF(printlevel.GT.2) THEN
          WRITE(*,*) 'Bessel function  for p'
          WRITE(*,'(A15,3X,A15)')
     $      'R in fm','BESF'
          DO irbess=1,nrbess
           WRITE(*,'(I6,E15.6,3X,6E15.6)')
     $          irbess,rbess(irbess),besf(irbess,1),besf(irbess,2)
          END DO
         END IF

C     V00
C     perform integration with splines 
         besspot=0.0
         DO irpot=1,nrpot
          DO irbess=1,nrbess
           besspot(irpot)=besspot(irpot)+
     $        besfp(irbess,1)*besf(irbess,1)
     $          *rspl(irbess,irpot)*wrbess(irbess)
     $          *rbess(irbess)**2
          END DO
         END DO

C do the interpolation for V
         V00=0.0
         DO irpot=1,nrpot
          V00=V00+besspot(irpot)*v_lsj(irpot,0,0)
         END DO

C     and for phi 
         momfakt=0.5**(pp**2+p**2)*hbarc**2
         DO irpot=1,nrpot
          V00=V00+momfakt*besspot(irpot)*fi_lsj(irpot,0,0)
         END DO

C     V02
C     perform integration with splines 
         besspot=0.0
         DO irpot=1,nrpot
          DO irbess=1,nrbess
           besspot(irpot)=besspot(irpot)+
     $        besfp(irbess,1)*besf(irbess,2)
     $          *rspl(irbess,irpot)*wrbess(irbess)
     $          *rbess(irbess)**2
          END DO
         END DO

C do the interpolation for V
         V02=0.0
         DO irpot=1,nrpot
          V02=V02+besspot(irpot)*v_lsj(irpot,0,1)
         END DO

C     and for phi 
         momfakt=0.5**(pp**2+p**2)*hbarc**2
         DO irpot=1,nrpot
          V02=V02+momfakt*besspot(irpot)*fi_lsj(irpot,0,1)
         END DO

C     V20
C     perform integration with splines 
         besspot=0.0
         DO irpot=1,nrpot
          DO irbess=1,nrbess
           besspot(irpot)=besspot(irpot)+
     $        besfp(irbess,2)*besf(irbess,1)
     $          *rspl(irbess,irpot)*wrbess(irbess)
     $          *rbess(irbess)**2
          END DO
         END DO

C do the interpolation for V
         V20=0.0
         DO irpot=1,nrpot
          V20=V20+besspot(irpot)*v_lsj(irpot,1,0)
         END DO

C     and for phi 
         momfakt=0.5**(pp**2+p**2)*hbarc**2
         DO irpot=1,nrpot
          V20=V20+momfakt*besspot(irpot)*fi_lsj(irpot,1,0)
         END DO

C     V22
C     perform integration with splines 
         besspot=0.0
         DO irpot=1,nrpot
          DO irbess=1,nrbess
           besspot(irpot)=besspot(irpot)+
     $        besfp(irbess,2)*besf(irbess,2)
     $          *rspl(irbess,irpot)*wrbess(irbess)
     $          *rbess(irbess)**2
          END DO
         END DO

C do the interpolation for V
         V22=0.0
         DO irpot=1,nrpot
          V22=V22+besspot(irpot)*v_lsj(irpot,1,1)
         END DO

C     and for phi 
         momfakt=0.5**(pp**2+p**2)*hbarc**2
         DO irpot=1,nrpot
          V22=V22+momfakt*besspot(irpot)*fi_lsj(irpot,1,1)
         END DO

         V00=V00*2.0/pi*p*pp*massfakt
         V02=-V02*2.0/pi*p*pp*massfakt
         V20=-V20*2.0/pi*p*pp*massfakt
         V22=V22*2.0/pi*p*pp*massfakt


        ELSE ! uncoupled 
C     prepare the besselfunctions for the Fourier transformation              

         CALL SBF1(oldl,pp,rbess,nrbess,besfp(:,1))

         IF(printlevel.GT.2) THEN
          WRITE(*,*) 'Bessel function  for pp'
          WRITE(*,'(A15,3X,A15)')
     $      'R in fm','BESFP'
          DO irbess=1,nrbess
           WRITE(*,'(I6,E15.6,3X,5E15.6)')
     $          irbess,rbess(irbess),besfp(irbess,1)
          END DO
         END IF

         CALL SBF1(oldl,p,rbess,nrbess,besf(:,1))

         IF(printlevel.GT.2) THEN
          WRITE(*,*) 'Bessel function  for p'
          WRITE(*,'(A15,3X,A15)')
     $      'R in fm','BESF'
          DO irbess=1,nrbess
           WRITE(*,'(I6,E15.6,3X,5E15.6)')
     $          irbess,rbess(irbess),besf(irbess,1)
          END DO
         END IF

C     perform integration with splines 
         besspot=0.0
         DO irpot=1,nrpot
          DO irbess=1,nrbess
           besspot(irpot)=besspot(irpot)+
     $        besfp(irbess,1)*besf(irbess,1)
     $          *rspl(irbess,irpot)*wrbess(irbess)
     $          *rbess(irbess)**2
          END DO
         END DO

C do the interpolation for V
         V00=0.0
         DO irpot=1,nrpot
          V00=V00+besspot(irpot)*v_lsj(irpot,0,0)
         END DO

C     and for phi 
         momfakt=0.5**(pp**2+p**2)*hbarc**2
         DO irpot=1,nrpot
          V00=V00+momfakt*besspot(irpot)*fi_lsj(irpot,0,0)
         END DO


         V00=V00*2.0/pi*p*pp*massfakt
         V20=0.0
         V02=0.0
         V22=0.0
        END IF ! coupled/uncoupled 


       END SUBROUTINE 


C     routine that prepare the lsj matrix elements from the operators 
C     call each time the partial wave channel changes 

      SUBROUTINE initav18chan
       IMPLICIT NONE
       REAL(dpreal) :: lsfakt,qsofakt,tfakt
       INTEGER irpot
       INTEGER t

       t=(1+(-1)**(olds+oldl))/2
C     set mass faktor for current channel
       
       SELECT CASE(oldchan)
        CASE(1) 
           massfakt=hbarc/41.4996
        CASE(2) 
           massfakt=hbarc/41.47105
        CASE(3)
           massfakt=hbarc/41.4425
        CASE DEFAULT
           WRITE(*,*) 'channel in potescnnlsj'
           CALL abort
       END SELECT

C     need to distinguish coupled channels 
C     from uncoupled 

       v_lsj(1:nrpot,0:1,0:1)=0.0 ! set partial wave potential to zero 
       fi_lsj(1:nrpot,0:1,0:1)=0.0

       IF(oldj.NE.0.AND.oldl.NE.oldj) THEN   ! coupled channels 

        v_lsj(1:nrpot,0,0)=v_c(1:nrpot,oldchan,t) ! add central part for current charge channel  l=j-1
        fi_lsj(1:nrpot,0,0)=fi_c(1:nrpot,oldchan,t)

        v_lsj(1:nrpot,1,1)=v_c(1:nrpot,oldchan,t) ! add central part for current charge channel  l=j+1 
        fi_lsj(1:nrpot,1,1)=fi_c(1:nrpot,oldchan,t)
 
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)        ! and spin-spin part l=j-1
     $     + (2*olds*(olds+1)-3)*v_s(1:nrpot,oldchan,t)
        fi_lsj(1:nrpot,0,0)=fi_lsj(1:nrpot,0,0)
     $     + (2*olds*(olds+1)-3)*fi_s(1:nrpot,oldchan,t)

        v_lsj(1:nrpot,1,1)=v_lsj(1:nrpot,1,1)        ! and spin-spin part l=j+1
     $     + (2*olds*(olds+1)-3)*v_s(1:nrpot,oldchan,t)
        fi_lsj(1:nrpot,1,1)=fi_lsj(1:nrpot,1,1)
     $     + (2*olds*(olds+1)-3)*fi_s(1:nrpot,oldchan,t)

        IF(olds.EQ.1) THEN
         tfakt=-2.*(oldj-1.)/(2.*oldj+1.)
         v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)        ! and tensor part if S=1 
     $      + tfakt*v_t(1:nrpot,oldchan,t)                  ! lp=j-1  l=j-1
         fi_lsj(1:nrpot,0,0)=fi_lsj(1:nrpot,0,0)
     $      + tfakt*fi_t(1:nrpot,oldchan,t)

         tfakt=6.*sqrt(oldj*(oldj+1.))/(2.*oldj+1)
         v_lsj(1:nrpot,1,0)=v_lsj(1:nrpot,1,0)        ! and tensor part if S=1 
     $      + tfakt*v_t(1:nrpot,oldchan,t)                  ! lp=j+1  l=j-1
         fi_lsj(1:nrpot,1,0)=fi_lsj(1:nrpot,1,0)
     $      + tfakt*fi_t(1:nrpot,oldchan,t)

         tfakt=6.*sqrt(oldj*(oldj+1.))/(2.*oldj+1)
         v_lsj(1:nrpot,0,1)=v_lsj(1:nrpot,0,1)        ! and tensor part if S=1 
     $      + tfakt*v_t(1:nrpot,oldchan,t)                  ! lp=j-1  l=j+1
         fi_lsj(1:nrpot,0,1)=fi_lsj(1:nrpot,0,1)
     $      + tfakt*fi_t(1:nrpot,oldchan,t)

         tfakt=-2.*(oldj+2.)/(2.*oldj+1.)
         v_lsj(1:nrpot,1,1)=v_lsj(1:nrpot,1,1)        ! and tensor part if S=1 
     $      + tfakt*v_t(1:nrpot,oldchan,t)                  ! lp=j+1  l=j+1
         fi_lsj(1:nrpot,1,1)=fi_lsj(1:nrpot,1,1)
     $      + tfakt*fi_t(1:nrpot,oldchan,t)
        END IF

        lsfakt=0.5*(oldj*(oldj+1)-oldl*(oldl+1)-olds*(olds+1))
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)         ! and spin-orbit part l=j-1
     $     + lsfakt*v_so(1:nrpot,oldchan,t)
        fi_lsj(1:nrpot,0,0)=fi_lsj(1:nrpot,0,0)
     $     + lsfakt*fi_so(1:nrpot,oldchan,t)


        qsofakt=2*lsfakt**2-oldl*(oldl+1)+lsfakt
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)         ! and nijm quadratic spin-orbit part
     $     + qsofakt*v_qso(1:nrpot,oldchan,t)           !    l=j-1
        fi_lsj(1:nrpot,0,0)=fi_lsj(1:nrpot,0,0)
     $     + qsofakt*fi_qso(1:nrpot,oldchan,t)

        qsofakt=(oldl*(oldl+1))
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)         ! and av18 quadratic l term 
     $     + qsofakt*v_l2(1:nrpot,oldchan,t)           !    l=j-1

        qsofakt=(oldl*(oldl+1))*(2*olds*(olds+1)-3)
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)         ! and av18 quadratic l spin-spin term 
     $     + qsofakt*v_l2s(1:nrpot,oldchan,t)           !    l=j-1

        qsofakt=lsfakt**2
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)         ! and av18 quadratic spin-orbit term 
     $     + qsofakt*v_ls2(1:nrpot,oldchan,t)           !    l=j-1
        
        lsfakt=0.5*(oldj*(oldj+1)-(oldl+2)*(oldl+3)-olds*(olds+1))
        v_lsj(1:nrpot,1,1)=v_lsj(1:nrpot,1,1)         ! and spin-orbit part l=j+1
     $     + lsfakt*v_so(1:nrpot,oldchan,t)
        fi_lsj(1:nrpot,1,1)=fi_lsj(1:nrpot,1,1)
     $     + lsfakt*fi_so(1:nrpot,oldchan,t)

        qsofakt=2*lsfakt**2-(oldl+2)*(oldl+3)+lsfakt
        v_lsj(1:nrpot,1,1)=v_lsj(1:nrpot,1,1)         ! and quadratic spin-orbit part
     $     + qsofakt*v_qso(1:nrpot,oldchan,t)           !    l=j+1
        fi_lsj(1:nrpot,1,1)=fi_lsj(1:nrpot,1,1)
     $     + qsofakt*fi_qso(1:nrpot,oldchan,t)

        qsofakt=((oldl+2)*(oldl+3))
        v_lsj(1:nrpot,1,1)=v_lsj(1:nrpot,1,1)         ! and av18 quadratic l term 
     $     + qsofakt*v_l2(1:nrpot,oldchan,t)           !    l=j+1

        qsofakt=((oldl+2)*(oldl+3))*(2*olds*(olds+1)-3)
        v_lsj(1:nrpot,1,1)=v_lsj(1:nrpot,1,1)         ! and av18 quadratic l spin-spin term 
     $     + qsofakt*v_l2s(1:nrpot,oldchan,t)           !    l=j+1

        qsofakt=lsfakt**2
        v_lsj(1:nrpot,1,1)=v_lsj(1:nrpot,1,1)         ! and av18 quadratic spin-orbit term 
     $     + qsofakt*v_ls2(1:nrpot,oldchan,t)           !    l=j+1

       ELSE                                  ! uncoupled channels 

        v_lsj(1:nrpot,0,0)=v_c(1:nrpot,oldchan,t) ! add central part for current charge channel
        fi_lsj(1:nrpot,0,0)=fi_c(1:nrpot,oldchan,t)
 
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)        ! and spin-spin part 
     $     + (2*olds*(olds+1)-3)*v_s(1:nrpot,oldchan,t)
        fi_lsj(1:nrpot,0,0)=fi_lsj(1:nrpot,0,0)
     $     + (2*olds*(olds+1)-3)*fi_s(1:nrpot,oldchan,t)

        IF(olds.EQ.1) THEN
         IF(oldl.ne.oldj) THEN
          tfakt=-2.*(oldj+2.)/(2.*oldj+1.)
         ELSE
          tfakt=2.
         END IF
         v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)        ! and tensor part if S=1 
     $      + tfakt*v_t(1:nrpot,oldchan,t)
         fi_lsj(1:nrpot,0,0)=fi_lsj(1:nrpot,0,0)
     $      + tfakt*fi_t(1:nrpot,oldchan,t)
        END IF

        lsfakt=0.5*(oldj*(oldj+1)-oldl*(oldl+1)-olds*(olds+1))
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)         ! and spin-orbit part
     $     + lsfakt*v_so(1:nrpot,oldchan,t)
        fi_lsj(1:nrpot,0,0)=fi_lsj(1:nrpot,0,0)
     $     + lsfakt*fi_so(1:nrpot,oldchan,t)

        qsofakt=2*lsfakt**2-oldl*(oldl+1)+lsfakt
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)         ! and quadratic s=n-orbit part
     $     + qsofakt*v_qso(1:nrpot,oldchan,t)
        fi_lsj(1:nrpot,0,0)=fi_lsj(1:nrpot,0,0)
     $     + qsofakt*fi_qso(1:nrpot,oldchan,t)
        
        qsofakt=(oldl*(oldl+1))
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)         ! and av18 quadratic l term 
     $     + qsofakt*v_l2(1:nrpot,oldchan,t)            

        qsofakt=(oldl*(oldl+1))*(2*olds*(olds+1)-3)
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)         ! and av18 quadratic l spin-spin term 
     $     + qsofakt*v_l2s(1:nrpot,oldchan,t)           

        qsofakt=lsfakt**2
        v_lsj(1:nrpot,0,0)=v_lsj(1:nrpot,0,0)         ! and av18 quadratic spin-orbit term 
     $     + qsofakt*v_ls2(1:nrpot,oldchan,t)           

       END IF                   ! coupled/uncoupled 

       IF(printlevel.GT.1) THEN
        WRITE(*,*) 'PW form of the potential (regular part)'
        WRITE(*,'(A15,3X,7A15)')
     $      'R in fm','V00','V01','V10','V00'
        DO irpot=1,nrpot
         WRITE(*,'(E15.6,3X,5E15.6)')
     $      rpot(irpot),v_lsj(irpot,0,0),
     $      v_lsj(irpot,0,1),v_lsj(irpot,1,0),v_lsj(irpot,1,1)
        END DO
        WRITE(*,*) 'PW form of the potential (derivative part)'
        WRITE(*,'(A15,3X,7A15)')
     $      'R in fm','V00','V01','V10','V00'
        DO irpot=1,nrpot
         WRITE(*,'(E15.6,3X,5E15.6)')
     $      rpot(irpot),fi_lsj(irpot,0,0),
     $      fi_lsj(irpot,0,1),fi_lsj(irpot,1,0),fi_lsj(irpot,1,1)
        END DO
       END IF


      END SUBROUTINE

C     routine that initializes the grids and the potential at the grid points 
C     and the splines 

       SUBROUTINE initav18
        IMPLICIT NONE
        INTEGER irbess,is,irpot,ichan,t
        REAL(dpreal),ALLOCATABLE :: spl(:,:)
        INTEGER,ALLOCATABLE :: indx(:,:)
C     interface to configuration space potential
        REAL(dpreal) VC,VSIG,VTEN,VSO,VSO2,FI,DFI,DDFI,                  
     .     FISS,DFISS,DDFISS,FIT,DFIT,DDFIT,                
     .     FISO,DFISO,DDFISO
        REAL(dpreal) VV(18),VEM(14)
        REAL(dpreal) :: isofakt(4,3)


C     read input file 'av18-para.dat'

        OPEN(UNIT=75,FILE='av18-para.dat',
     $      FORM='FORMATTED',STATUS='OLD')
        READ(75,*,END=100,ERR=100)
        READ(75,*,END=100,ERR=100) imodel,empotflag                  
        READ(75,*,END=100,ERR=100) r1pot,r2pot,r3pot,nr1pot,nr2pot 
        READ(75,*,END=100,ERR=100) r1bess,r2bess,r3bess,nr1bess,nr2bess 
        CLOSE(75)

        nrpot=nr1pot+nr2pot 
        nrbess=nr1bess+nr2bess 
        
        WRITE(*,*) 'Parameters of the Fouriertransform of'
        WRITE(*,*) 'AV18 imodel-empot: ',imodel,empotflag
        WRITE(*,'(A,3E15.6,3I5)') 'RPOT:  ',r1pot,r2pot,r3pot,
     $                             nr1pot,nr2pot,nrpot 
        WRITE(*,'(A,3E15.6,3I5)') 'RBESS: ',r1bess,r2bess,r3bess,
     $                             nr1bess,nr2bess,nrbess 

C     now allocate storage for mesh and potential
C     meshes 
        IF(allocated(rpot)) DEALLOCATE(rpot)
        ALLOCATE(rpot(nrpot))
        IF(allocated(wrpot)) DEALLOCATE(wrpot)
        ALLOCATE(wrpot(nrpot))

        IF(allocated(rbess)) DEALLOCATE(rbess)
        ALLOCATE(rbess(nrbess))
        IF(allocated(wrbess)) DEALLOCATE(wrbess)
        ALLOCATE(wrbess(nrbess))
C     splines
        IF(allocated(rspl)) DEALLOCATE(rspl)
        ALLOCATE(rspl(nrbess,nrpot))
C     besselfunctions etc 
        IF(allocated(besf)) DEALLOCATE(besf)
        ALLOCATE(besf(nrbess,2))
        IF(allocated(besfp)) DEALLOCATE(besfp)
        ALLOCATE(besfp(nrbess,2))
        IF(allocated(besspot)) DEALLOCATE(besspot)
        ALLOCATE(besspot(nrpot))

C     potential operators 
        IF(allocated(v_c)) DEALLOCATE(v_c)
        ALLOCATE(v_c(nrpot,3,0:1))
        IF(allocated(v_t)) DEALLOCATE(v_t)
        ALLOCATE(v_t(nrpot,3,0:1))
        IF(allocated(v_so)) DEALLOCATE(v_so)
        ALLOCATE(v_so(nrpot,3,0:1))
        IF(allocated(v_s)) DEALLOCATE(v_s)
        ALLOCATE(v_s(nrpot,3,0:1))
        IF(allocated(v_qso)) DEALLOCATE(v_qso)
        ALLOCATE(v_qso(nrpot,3,0:1))
        IF(allocated(v_l2)) DEALLOCATE(v_l2)
        ALLOCATE(v_l2(nrpot,3,0:1))
        IF(allocated(v_l2s)) DEALLOCATE(v_l2s)
        ALLOCATE(v_l2s(nrpot,3,0:1))
        IF(allocated(v_ls2)) DEALLOCATE(v_ls2)
        ALLOCATE(v_ls2(nrpot,3,0:1))

        IF(allocated(fi_c)) DEALLOCATE(fi_c)
        ALLOCATE(fi_c(nrpot,3,0:1))
        IF(allocated(fi_t)) DEALLOCATE(fi_t)
        ALLOCATE(fi_t(nrpot,3,0:1))
        IF(allocated(fi_so)) DEALLOCATE(fi_so)
        ALLOCATE(fi_so(nrpot,3,0:1))
        IF(allocated(fi_s)) DEALLOCATE(fi_s)
        ALLOCATE(fi_s(nrpot,3,0:1))
        IF(allocated(fi_qso)) DEALLOCATE(fi_qso)
        ALLOCATE(fi_qso(nrpot,3,0:1))        
C     potential in partial wave basis
        IF(allocated(v_lsj)) DEALLOCATE(v_lsj)
        ALLOCATE(v_lsj(nrpot,0:1,0:1))        
        IF(allocated(fi_lsj)) DEALLOCATE(fi_lsj)
        ALLOCATE(fi_lsj(nrpot,0:1,0:1))        
        
C     prepare meshed according to read parameters 
        CALL TRNS(nr1pot,nr2pot,nrpot,
     $            r1pot,r2pot,r3pot,rpot,wrpot)
        CALL TRNS(nr1bess,nr2bess,nrbess,
     $            r1bess,r2bess,r3bess,rbess,wrbess)

C     print mesh if printlevel large enough 
        IF(printlevel.GT.0) THEN
         WRITE(*,*) 'R-Grid for potential'
         DO irpot=1,nrpot
          WRITE(*,'(A,I5,A,E15.6,A,I5,A,E15.6)')
     $                   'RPOT(',irpot,')= ',rpot(irpot),
     $              '    WRPOT(',irpot,')= ',wrpot(irpot)
         END DO

         WRITE(*,*) 'R-Grid for bessel functions'
         DO irbess=1,nrbess
          WRITE(*,'(A,I5,A,E15.6,A,I5,A,E15.6)')
     $                   'RBESS(',irbess,')= ',rbess(irbess),
     $              '    WRBESS(',irbess,')= ',wrbess(irbess)
         END DO
        END IF 

C     prepare splines using hermitian splines but storing like basis splines 
        
        ALLOCATE(spl(nrbess,4),indx(nrbess,4))
        CALL cubherm_dp(rpot,nrpot,rbess,nrbess,spl,indx)
C     prepare the rspl array that include zeros, where the index is not defined 
        rspl=0.0
        DO irbess=1,nrbess
         DO is=1,4
          rspl(irbess,indx(irbess,is))
     $       =rspl(irbess,indx(irbess,is))+spl(irbess,is)
         END DO
        END DO
        DEALLOCATE(spl,indx)

C     now prepare the isospin factors assuming that I do not need to symmetrize particle space 


C     now prepare the configuration space potential 
        DO irpot=1,nrpot
         CALL av18op(rpot(irpot),VV)
         IF(empotflag.NE.0) THEN
          CALL empot(1,rpot(irpot),VEM)
         ELSE 
          VEM=0.0_dpreal
         END IF

C     convert from to fm-1 
         VV=VV/hbarc
         VEM=VEM/hbarc
         
         DO t=0,1               !! isospin

          isofakt(1,1:3)=1.0_dpreal !! scalar 

          IF(t.EQ.0) THEN
           isofakt(2,1:3)=-3.0_dpreal !! tau1*tau2 
          ELSE
           isofakt(2,1:3)=1.0_dpreal !! tau1*tau2 
          END IF

          isofakt(3,1)=3.0_dpreal-isofakt(2,1)  !! T12
          isofakt(3,2)=-3.0_dpreal-isofakt(2,2)  !! T12
          isofakt(3,3)=3.0_dpreal-isofakt(2,3)  !! T12

          isofakt(4,1)=2.0_dpreal !! T1Z+T2Z
          isofakt(4,2)=0.0_dpreal !! T1Z+T2Z
          isofakt(4,3)=-2.0_dpreal !! T1Z+T2Z

          DO ichan=1,3          !! pp,np and nn 
          
          

c order of operators l in vnn(l):
c l:    1=1                              2=t1.t2
c       3=s1.s2                          4=(s1.s2)(t1.t2)
c       5=S12 [=3(s1.r)(s2.r)-s1.s2]     6=S12(t1.t2)
c       7=L.S                            8=L.S(t1.t2)
c       9=L**2                          10=L**2(t1.t2)
c      11=L**2(s1.s2)                   12=L**2(s1.s2)(t1.t2)
c      13=(L.S)**2                      14=(L.S)**2(t1.t2)
c      15=T12 [=3*t1z*t2z-t1.t2]        16=(s1.s2)T12
c      17=S12*T12                       18=t1z+t2z

c order of operators in vem(l)
c l:    1=C1    (pp)          2=DF    (pp)          3=C2      (pp)
c       4=VP    (pp)                                5=C1      (np)
c       6=s1.s2 (pp)          7=s1.s2 (nn)          8=s1.s2   (np)
c       9=S12   (pp)         10=S12   (nn)         11=S12     (np)
c      12=L.S   (pp)         13=L.S   (nn)         14=L.S     (np)
c C1 = one-photon-exchange Coulomb with form factor
c C2 = two-photon-exchange Coulomb
c DF = Darwin-Foldy
c VP = vacuum polarization (short-range approximation)
c all other terms from magnetic moment (MM) interactions
c ----------------------------------------------------------------------

C     set all central parts 
           v_c(irpot,ichan,t)=isofakt(1,ichan)*VV(1)
     $                       +isofakt(2,ichan)*VV(2)
     $                       +isofakt(3,ichan)*VV(15)
     $                       +isofakt(4,ichan)*VV(18)
           SELECT CASE(ichan)
            CASE(1) ! pp 
             v_c(irpot,ichan,t)=v_c(irpot,ichan,t)
     $                         +VEM(1)+VEM(2)+VEM(3)+VEM(4)           
            CASE(2) ! np 
             v_c(irpot,ichan,t)=v_c(irpot,ichan,t)
     $                         +VEM(5)           
            CASE(3) ! nn 
             v_c(irpot,ichan,t)=v_c(irpot,ichan,t)          
            CASE DEFAULT
             WRITE(*,*) 'channel in potav18lsj'
             CALL abort
           END SELECT

C     set all spin-spin parts 

           v_s(irpot,ichan,t)=isofakt(1,ichan)*VV(3)
     $                       +isofakt(2,ichan)*VV(4)
     $                       +isofakt(3,ichan)*VV(16)

           SELECT CASE(ichan)
            CASE(1) ! pp 
             v_s(irpot,ichan,t)=v_s(irpot,ichan,t)+VEM(6)
            CASE(2) ! np 
             v_s(irpot,ichan,t)=v_s(irpot,ichan,t)+VEM(8)           
            CASE(3) ! nn 
             v_s(irpot,ichan,t)=v_s(irpot,ichan,t)+VEM(7)           
            CASE DEFAULT
             WRITE(*,*) 'channel in potav18lsj'
             CALL abort
           END SELECT


C     set all tensor parts 

           v_t(irpot,ichan,t)=isofakt(1,ichan)*VV(5)
     $                     +isofakt(2,ichan)*VV(6)
     $                     +isofakt(3,ichan)*VV(17)

           SELECT CASE(ichan)
            CASE(1) ! pp 
             v_t(irpot,ichan,t)=v_t(irpot,ichan,t)+VEM(9)
            CASE(2) ! np 
             v_t(irpot,ichan,t)=v_t(irpot,ichan,t)+VEM(11)           
            CASE(3) ! nn 
             v_t(irpot,ichan,t)=v_t(irpot,ichan,t)+VEM(10)           
            CASE DEFAULT
             WRITE(*,*) 'channel in potav18lsj'
             CALL abort
           END SELECT


C     set all spin-orbit parts 

           v_so(irpot,ichan,t)=isofakt(1,ichan)*VV(7)
     $                      +isofakt(2,ichan)*VV(8)

           SELECT CASE(ichan)
            CASE(1) ! pp 
             v_so(irpot,ichan,t)=v_so(irpot,ichan,t)+VEM(12)
            CASE(2) ! np 
             v_so(irpot,ichan,t)=v_so(irpot,ichan,t)+VEM(14)           
            CASE(3) ! nn 
             v_so(irpot,ichan,t)=v_so(irpot,ichan,t)+VEM(13)           
            CASE DEFAULT
             WRITE(*,*) 'channel in potav18lsj'
             CALL abort
           END SELECT


C     set no nijm qso  parts 

           v_qso(irpot,ichan,t)=0.0


C     set all av18 l2 parts 

           v_l2(irpot,ichan,t)=isofakt(1,ichan)*VV(9)
     $                     +isofakt(2,ichan)*VV(10)


C     set all av18 l2 spin-spin parts 

           v_l2s(irpot,ichan,t)=isofakt(1,ichan)*VV(11)
     $                      +isofakt(2,ichan)*VV(12)

C     set all av18 ls2 parts 

           v_ls2(irpot,ichan,t)=isofakt(1,ichan)*VV(13)
     $                      +isofakt(2,ichan)*VV(14)


           fi_c(irpot,ichan,t)=0.0_dpreal
           fi_s(irpot,ichan,t)=0.0_dpreal
           fi_t(irpot,ichan,t)=0.0_dpreal
           fi_so(irpot,ichan,t)=0.0_dpreal
           fi_qso(irpot,ichan,t)=0.0_dpreal
          END DO                 ! ichan
         END DO                  ! t
        END DO                  ! irpot

        IF(printlevel.GT.1) THEN
         WRITE(*,*) 'Operator form of the potential (regular part)'
         WRITE(*,'(2A5,A15,3X,8A15)')
     $      'ICH','T','R in fm','VC','VS','VT','VSO','VQSO','VL2','VL2S','VLS2'
         DO t=0,1               !! isospin 
          DO ichan=1,3          !! pp,np and nn 
           DO irpot=1,nrpot
            WRITE(*,'(2I5,E15.6,3X,8E15.6)')
     $               ichan,t,rpot(irpot),v_c(irpot,ichan,t),
     $               v_s(irpot,ichan,t),v_t(irpot,ichan,t),
     $               v_so(irpot,ichan,t),v_qso(irpot,ichan,t),
     $               v_l2(irpot,ichan,t),v_l2s(irpot,ichan,t),
     $               v_ls2(irpot,ichan,t)
           END DO
          END DO
         END DO
         WRITE(*,*) 'Operator form of the potential (derivative part)'
         WRITE(*,'(2A5,A15,3X,5A15)')
     $      'ICH','T','R in fm','FIC','FIS','FIT','FISO','FIQSO'
         DO t=0,1               !! isospin 
          DO ichan=1,3          !! pp,np and nn 
           DO irpot=1,nrpot
            WRITE(*,'(2I5,E15.6,3X,5E15.6)')
     $         ichan,t,rpot(irpot),fi_c(irpot,ichan,t),
     $         fi_s(irpot,ichan,t),fi_t(irpot,ichan,t),
     $         fi_so(irpot,ichan,t),fi_qso(irpot,ichan,t)
           END DO
          END DO
         END DO
        END IF

        RETURN
 100    CONTINUE
        STOP 'problem with av18.para.dat'

       END SUBROUTINE initav18

      
       SUBROUTINE SBF1(LM,P,R,NR,BL)
        IMPLICIT NONE
        INTEGER LM,NR
        REAL(dpreal) P,X
        REAL(dpreal) BL(NR)
        REAL(dpreal) R(NR)
        REAL(dpreal),SAVE :: DF(8) 
        DATA DF/1.,3.,15.,105.,945.,10395.,135135. ,2027025. /
        INTEGER MRM,IR,III
C
        III=1
        IF(P*R(1).EQ.0. .AND. LM.EQ.0 ) THEN
         III=2
         BL(1)=1.0
        END IF
        IF(P*R(1).EQ.0. .AND. LM.NE.0 ) THEN
         III=2
         BL(1)=0.0
        END IF
C
C      MRM=NR+1
        MRM=NR
        DO IR=1,NR
         IF(P*R(IR).GT.LM*0.2) THEN
          MRM=IR
          GOTO 1000
         END IF
        END DO
 1000   CONTINUE
        DO IR=III,MRM
         X=P*R(IR)
         BL(IR)
     &      =(X**LM/DF(LM+1))*(1.-0.5*X**2/(2.*LM+3.)+
     1      (0.5*X**2)**2/(2.*(2.*LM+3.)*(2.*LM+5.))-
     2      (0.5*X**2)**3/(6.*(2.*LM+3.)*(2.*LM+5.)*(2.*LM+7.))+
     3      (0.5*X**2)**4/(24.*(2.*LM+3.)*(2.*LM+5.)*(2.*LM+7.)*(2.*
     4      LM+9.)))
        END DO
 10     GOTO(1,2,3,4,5,6,7,8)LM+1
 1      DO IR=MRM,NR
         X=P*R(IR)
         BL(IR)=SIN(X)/X
        END DO
        GOTO 2000
 2      DO IR=MRM,NR
         X=P*R(IR)
         BL(IR)=SIN(X)/X**2-COS(X)/X
        END DO
        GOTO 2000
 3      DO IR=MRM,NR
        X=P*R(IR)
        BL(IR)=(3./X**3-1./X)*SIN(X)-3.*COS(X)/X**2
       END DO
       GOTO 2000
 4     DO IR=MRM,NR
        X=P*R(IR)
        BL(IR)=(15./X**4-6./X**2)*SIN(X)-(15./X**3-1./X)*COS(X)
       END DO
       GOTO 2000
 5     DO IR=MRM,NR
        X=P*R(IR)
        BL(IR)=(105./X**5-45./X**3+1./X)*SIN(X)+(10./X**2-
     1     105./X**4)*COS(X)
       END DO
       GOTO 2000
 6     DO IR=MRM,NR
        X=P*R(IR)
        BL(IR)=(15./X**2-420./X**4+945./X**6)*SIN(X)+(-1./X+
     1     105./X**3-945./X**5)*COS(X)
       END DO
       GOTO 2000
 7     DO IR=MRM,NR
        X=P*R(IR)
        BL(IR)=
     -     ((-10395./x**5 + 1260./x**3 - 21./x)*COS(X) +(- 1. + 
     -     10395./x**6 - 4725./x**4 + 210./x**2)
     -     *Sin(x) )/x
       END DO
       GOTO 2000
 8     DO IR=MRM,NR
        X=P*R(IR)
        BL(IR)=
     -     ((1. - 135135./x**6 + 17325./x**4 - 378./x**2)*Cos(x) + 
     -     (  135135./x**7 - 62370./x**5 + 3150./x**3 - 
     -     28./x)*Sin(x))/x
       END DO

 2000  CONTINUE
       RETURN
      END SUBROUTINE

C====================================================================C
C     AV18 ORIGINAL AREA IS STARTING.
C     1994.2.11 : neuteron-neuteron parts are added.
C     Change point is indicated by *HK*
C====================================================================C
C
C
c *id* av18op **********************************************************
c subroutine for strong interaction part of argonne v18 potential 
c in operator format
c calls subroutine consts
c ----------------------------------------------------------------------
c arguments for av18pot
c r:   separation in fm
c vnn: output potential in MeV (18 component array)
c ----------------------------------------------------------------------
c order of operators l in vnn(l):
c l:    1=1                              2=t1.t2
c       3=s1.s2                          4=(s1.s2)(t1.t2)
c       5=S12 [=3(s1.r)(s2.r)-s1.s2]     6=S12(t1.t2)
c       7=L.S                            8=L.S(t1.t2)
c       9=L**2                          10=L**2(t1.t2)
c      11=L**2(s1.s2)                   12=L**2(s1.s2)(t1.t2)
c      13=(L.S)**2                      14=(L.S)**2(t1.t2)
c      15=T12 [=3*t1z*t2z-t1.t2]        16=(s1.s2)T12
c      17=S12*T12                       18=t1z+t2z
c where s1=sigma_1, t1=tau_1, t1z=tau_1(z), etc.
c ----------------------------------------------------------------------
      subroutine av18op(r,vnn)
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      dimension vnn(18)
      real*8 mpi0,mpic,mp,mn,mup,mun
      real*8 mpi,mu0,muc,mu
      data small/1e-4/
      do 5 l=1,18
        vnn(l)=0
    5 continue
      call consts(hc,mpi0,mpic,mp,mn,alpha,mup,mun)
      mpi=(mpi0+2.*mpic)/3.
      mu0=mpi0/hc
      muc=mpic/hc
      mu=mpi/hc
      fsq=.075
      cpi=2.1
      rws=.5
      aiws=5.
      x=mu*r
      x0=mu0*r
      xc=muc*r
      if (r.le.small) then
        tpi=3*cpi**2*r/mu**3
        ypi0=(mpi0/mpic)**2*(mpi0/3)*cpi*r/mu0
        tpi0=3*cpi*ypi0/mu0**2
        ypic=(mpic/3)*cpi*r/muc
        tpic=3*cpi*ypic/muc**2
      else
        rcut=1-exp(-cpi*r*r)
        ypi=exp(-x)*rcut/x
        tpi=(1+(3+3/x)/x)*ypi*rcut
        ypi0=(mpi0/mpic)**2*(mpi0/3)*exp(-x0)*rcut/x0
        tpi0=(1+(3+3/x0)/x0)*ypi0*rcut
        ypic=(mpic/3)*exp(-xc)*rcut/xc
        tpic=(1+(3+3/xc)/xc)*ypic*rcut
      end if
      ypi0=fsq*ypi0
      ypic=fsq*ypic
      tpi0=fsq*tpi0
      tpic=fsq*tpic
      tpi2=tpi*tpi
      ws=1/(1+exp((r-rws)*aiws))
      ws0=1/(1+exp(-rws*aiws))
      wsp=ws*(1+aiws*exp(-rws*aiws)*ws0*r)
      wsx=ws*x
      wsx2=wsx*x
      dypi00=(mpi0/mpic)**2*(mpi0/3)*cpi/mu0
      dypic0=(mpic/3)*cpi/muc
      ypi0p=ypi0-fsq*dypi00*ws*r/ws0
      ypicp=ypic-fsq*dypic0*ws*r/ws0
      ypi=(ypi0+2*ypic)/3
      tpi=(tpi0+2*tpic)/3
      p11pp=  -7.62701*tpi2+1815.4920*wsp+1847.8059*wsx2+ypi0p
      p11np=  -7.62701*tpi2+1813.5315*wsp+1847.8059*wsx2-ypi0p+2*ypicp
      p11nn=  -7.62701*tpi2+1811.5710*wsp+1847.8059*wsx2+ypi0p
      pt1pp=   1.07985*tpi2 -190.0949*wsx -811.2040*wsx2+tpi0
      pt1np=   1.07985*tpi2 -190.0949*wsx -811.2040*wsx2-tpi0+2*tpic
      pt1nn=   1.07985*tpi2 -190.0949*wsx -811.2040*wsx2+tpi0
      pls1=    -.62697*tpi2 -570.5571*wsp +819.1222*wsx2
      pl211=    .06709*tpi2 +342.0669*wsp -615.2339*wsx2
      pls21=    .74129*tpi2   +9.3418*wsp -376.4384*wsx2
      p10=    -8.62770*tpi2+2605.2682*wsp +441.9733*wsx2-ypi0p-2*ypicp
      pt0=    1.485601*tpi2-1126.8359*wsx +370.1324*wsx2-tpi0-2*tpic
      pls0=     .10180*tpi2  +86.0658*wsp -356.5175*wsx2
      pl210=   -.13201*tpi2 +253.4350*wsp   -1.0076*wsx2
      pls20=    .07357*tpi2 -217.5791*wsp  +18.3935*wsx2
      p01pp= -11.27028*tpi2+3346.6874*wsp-3*ypi0p
      p01np= -10.66788*tpi2+3126.5542*wsp-3*(-ypi0p+2*ypicp)
      p01nn= -11.27028*tpi2+3342.7664*wsp-3*ypi0p
      pl201=    .12472*tpi2  +16.7780*wsp
      p00=    -2.09971*tpi2+1204.4301*wsp-3*(-ypi0p-2*ypicp)
      pl200=   -.31452*tpi2 +217.4559*wsp
      p11=(p11pp+p11nn+p11np)/3
      p11cd=(.5*(p11pp+p11nn)-p11np)/6
      p11cs=(p11pp-p11nn)/4
      pt1=(pt1pp+pt1nn+pt1np)/3
      pt1cd=(.5*(pt1pp+pt1nn)-pt1np)/6
      pt1cs=(pt1pp-pt1nn)/4
      p01=(p01pp+p01nn+p01np)/3
      p01cd=(.5*(p01pp+p01nn)-p01np)/6
      p01cs=(p01pp-p01nn)/4
      vnn(1)=.0625*(9*p11+3*p10+3*p01+p00)
      vnn(2)=.0625*(3*p11-3*p10  +p01-p00)
      vnn(3)=.0625*(3*p11  +p10-3*p01-p00)
      vnn(4)=.0625*(  p11  -p10  -p01+p00)
      vnn(5)=.25*(3*pt1+pt0)
      vnn(6)=.25*(  pt1-pt0)
      vnn(7)=.25*(3*pls1+pls0)
      vnn(8)=.25*(  pls1-pls0)
      vnn(9)= .0625*(9*pl211+3*pl210+3*pl201+pl200)
      vnn(10)=.0625*(3*pl211-3*pl210+  pl201-pl200)
      vnn(11)=.0625*(3*pl211+  pl210-3*pl201-pl200)
      vnn(12)=.0625*(  pl211-  pl210-  pl201+pl200)
      vnn(13)=.25*(3*pls21+pls20)
      vnn(14)=.25*(  pls21-pls20)
      vnn(15)=.25*(3*p11cd+p01cd)
      vnn(16)=.25*(  p11cd-p01cd)
      vnn(17)=pt1cd
      vnn(18)=p01cs
      return
      end subroutine
c *id* consts **********************************************************
c subroutine for constants in av18 potential
c ----------------------------------------------------------------------
c arguments for consts
c hc:    output value for hbar*c (MeV-fm)
c mpi0:    "      "    "  neutral pion mass (MeV)
c mpic:    "      "    "  charged pion mass (MeV)
c mp:      "      "    "  proton mass (MeV)
c mn:      "      "    "  neutron mass (MeV)
c alpha:   "      "    "  electromagnetic constant alpha
c mup:     "      "    "  proton magnetic moment (nm)
c mun:     "      "    "  neutron magnetic moment (nm)
c ----------------------------------------------------------------------
      subroutine consts(hc,mpi0,mpic,mp,mn,alpha,mup,mun)
      implicit real*8 (a-h,o-z)
      real*8 hc,mpi0,mpic,mp,mn,alpha,mup,mun
      hc=197.327053
      mpi0=134.9739
      mpic=139.5675
      mp=938.27231
      mn=939.56563
      alpha=1./137.035989
      mup= 2.7928474
      mun=-1.9130427
      return
      end subroutine
C===================================================================C
C   ORIGINAL CODE END
C
C===================================================================C
C
c *id* empot ***********************************************************
c subroutine for electromagnetic potential
c arguments for empot
c lemp: 0=Coulomb w/ff only
c       1=full electromagnetic potential
c rr:   position r in fm
c vem:  returned potential in MeV (14 component array)
c ----------------------------------------------------------------------
c order of operators in vem(l)
c l:    1=C1    (pp)          2=DF    (pp)          3=C2      (pp)
c       4=VP    (pp)                                5=C1      (np)
c       6=s1.s2 (pp)          7=s1.s2 (nn)          8=s1.s2   (np)
c       9=S12   (pp)         10=S12   (nn)         11=S12     (np)
c      12=L.S   (pp)         13=L.S   (nn)         14=L.S     (np)
c C1 = one-photon-exchange Coulomb with form factor
c C2 = two-photon-exchange Coulomb
c DF = Darwin-Foldy
c VP = vacuum polarization (short-range approximation)
c all other terms from magnetic moment (MM) interactions
c ----------------------------------------------------------------------
      subroutine empot(lemp,rr,vem)
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension vem(14)
      real*8 kr,me,mp,mn,mr,mun,mup
      data small/1e-4/,vsmall/1e-10/
      alpha=1./137.03599
      hc=197.327053
      b=4.27
      rrsave=rr
      if (rr.le.small) rr=small
      br=b*rr
      fcoul=1-(1+11*br/16+3*br**2/16+br**3/48)*exp(-br)
      vem(1)=alpha*hc*fcoul/rr
      if (lemp.ge.1) then
CAN        pi=acos(-1.) defined in precision
        me=0.510999
        mp=938.27231
        mn=939.56563
        mr=mp*mn/(mp+mn)
        mup=2.79285
        mun=-1.91304
        gamma=0.577216
        beta=.0189
        kr=me*rr/hc
        ftensor=1-(1+br+br**2/2+br**3/6+br**4/24+br**5/144)*exp(-br)
        fspinor=1-(1+br+br**2/2+7*br**3/48+br**4/48)*exp(-br)
        fdelta=b**3*(1+br+br**2/3)*exp(-br)/16
        fnp=b**2*(15*br+15*br**2+6*br**3+br**4)*exp(-br)/384
        fivp=-gamma-5./6.+abs(log(kr))+6*pi*kr/8
        vem(2)=-alpha*hc**3*fdelta/(4*mp**2)
        vem(3)=-vem(1)**2/mp
        vem(4)=2*alpha*vem(1)*fivp/(3*pi)
        vem(5)=alpha*hc*beta*fnp/rr
        vem(6)=-alpha*hc**3*mup**2*fdelta/(6*mp**2)
        vem(7)=-alpha*hc**3*mun**2*fdelta/(6*mn**2)
        vem(8)=-alpha*hc**3*mup*mun*fdelta/(6*mn*mp)
        vem(9)=-alpha*hc**3*mup**2*ftensor/(4*mp**2*rr**3)
        vem(10)=-alpha*hc**3*mun**2*ftensor/(4*mn**2*rr**3)
        vem(11)=-alpha*hc**3*mup*mun*ftensor/(4*mp*mn*rr**3)
        vem(12)=-alpha*hc**3*(4*mup-1)*fspinor/(2*mp**2*rr**3)
        vem(13)=0
        vem(14)=-alpha*hc**3*mun*fspinor/(2*mn*mr*rr**3)
      else
        do 10 i=2,14
          vem(i)=0
   10   continue
      end if
      rr=rrsave
      return
      end subroutine

      END MODULE argonne
