************************************************************************ * Reid93 soft core phenomenological potential, updated version, * including one-pion-exchange with neutral-pion and charged-pion * masses; coupling f^2=0.075. Tensor potential is regularized * to equal zero at r=0 * Reference: V.G.J. Stoks et al., Phys. Rev. C 49, 2950 (1994). *----------------------------------------------------------------------- * Fortran 77 version 1.0: July 1994 * version 1.1: March 1995 * version 1.2: August 1995 * version 2.0: April 2000 * * E-mail: info@nn-online.org * WWW: http://nn-online.org *----------------------------------------------------------------------- * IN: real*8 r = distance in fermi * character*2 type = reaction, 'PP', 'NP', 'PN', or 'NN' * character*3 name = name of the partial wave (see below), * maximum total angular momentum j = 9 * OUT: real*8 vpot(2,2) = potential matrix in MeV on LSJ-basis *----------------------------------------------------------------------- * The variable 'name' contains the name of the partial wave in spectral * notation: singlets: 1S0 1P1 1D2 1F3 1G4 ... * triplets uncoupled: 3P0 3P1 3D2 3F3 3G4 ... * triplets coupled: 3C1 3C2 3C3 3C4 ... * where 3C1 denotes 3S1 - epsilon_1 - 3D1 channel * 3C2 denotes 3P2 - epsilon_2 - 3F2 channel * ... ************************************************************************ subroutine rreid93(r,name,type,vpot) implicit none real*8 r,vpot(2,2) character*2 type character*3 name character*12 large,small integer l,spin,j,nchan,iso logical first real*8 ri,rj,vc,vt,vls,vtpi,x,x0,xc,mpi,vspi,vspis real*8 r1,r2,r3,r4,r6 real*8 f0pi,fcpi,hbc,mpi0,mpic,mpis real*8 a(5,5),b(5,5),parspp(5,5),parsnp(5,5) real*8 r93_y,r93_yp,r93_w,r93_z save a,b data f0pi /0.075d0/ data fcpi /0.075d0/ data hbc /197.327053d0/ data mpi0 /134.9739d0/ data mpic /139.5675d0/ data mpis /139.5675d0/ data r1 /1.0d0/ data r2 /2.0d0/ data r3 /3.0d0/ data r4 /4.0d0/ data r6 /6.0d0/ data first /.true./ data large /'CSPDFGHIKLMN'/ data small /'cspdfghiklmn'/ data parspp ./ 0.1756084d0,-0.1414234d2, 0.1518489d3,-0.6868230d3, 0.1104157d4, . -0.4224976d2, 0.2072246d3,-0.3354364d3,-0.1989250d1,-0.6178469d2, . 0.2912845d2, 0.1511690d3, 0.8151964d1, 0.5832103d2,-0.2074743d2, . -0.5840566d0,-0.1029310d2, 0.2263391d2, 0.2316915d2,-0.1959172d1, . -0.2608488d1, 0.1090858d2,-0.4374212d0,-0.2148862d2,-0.6584788d0/ data parsnp ./-0.2234989d2, 0.2551761d3,-0.1063549d4, 0.1609196d4,-0.3505968d1, . -0.4248612d1,-0.5352001d1, 0.1827642d3,-0.3927086d3, 0.5812273d2, . -0.2904577d1, 0.3802497d2, 0.3395927d0, 0.8318097d0, 0.1923895d1, . 0.0913746d0,-0.1274773d2, 0.1458600d3,-0.6432461d3, 0.1022217d4, . -0.0461640d0, 0.7950192d1,-0.1925573d1, 0.5066234d2, 0.8359896d1/ * The parameters if (first) then do j=1,5 do l=1,5 a(j,l) = parspp(l,j) b(j,l) = parsnp(l,j) enddo enddo first = .false. endif * Determine the quantumnumbers do l=1,12 if (name(2:2).eq.small(l:l)) name(2:2) = large(l:l) enddo nchan = 1 if (name(2:2).eq.'C') nchan = 2 if (name(1:1).eq.'1') spin = 0 if (name(1:1).eq.'3') spin = 1 read(name,'(2x,i1)') j l = j if (name.eq.'3P0') l = 1 if (nchan.eq.2) l = j - 1 iso = mod(spin+l+1,2) ri = real(iso) rj = real(j) * OPE potential x0 = mpi0/hbc * r vspis = f0pi*(mpi0/mpis)**2*mpi0/r3*r93_yp(1,8,x0) vspi = f0pi*(mpi0/mpis)**2*mpi0/r3*r93_y(1,8,x0) vtpi = f0pi*(mpi0/mpis)**2*mpi0/r3*r93_z(1,8,x0) if (type.eq.'NP' .or. type.eq.'PN') then xc = mpic/hbc * r vspis = (r4*ri-r2)*fcpi*(mpic/mpis)**2*mpic/r3*r93_yp(1,8,xc) . - vspis vspi = (r4*ri-r2)*fcpi*(mpic/mpis)**2*mpic/r3*r93_y(1,8,xc) . - vspi vtpi = (r4*ri-r2)*fcpi*(mpic/mpis)**2*mpic/r3*r93_z(1,8,xc) . - vtpi endif * Start clean vpot(1,1) = 0.0d0 vpot(1,2) = 0.0d0 vpot(2,1) = 0.0d0 vpot(2,2) = 0.0d0 * Non-OPE contribution to the potential mpi = (mpi0 + r2*mpic)/r3 x = mpi/hbc * r if (name.eq.'1S0') then if (type.eq.'PP' .or. type.eq.'NN') then vpot(1,1) = a(1,1)*r93_y(2,8,x) . + a(1,2)*r93_y(3,8,x) . + a(1,3)*r93_y(4,8,x) . + a(1,4)*r93_y(5,8,x) . + a(1,5)*r93_y(6,8,x) elseif (type.eq.'NP' .or. type.eq.'PN') then vpot(1,1) = b(1,1)*r93_y(3,8,x) . + b(1,2)*r93_y(4,8,x) . + b(1,3)*r93_y(5,8,x) . + b(1,4)*r93_y(6,8,x) endif elseif (name.eq.'1D2') then vpot(1,1) = a(2,1)*r93_y(4,8,x) . + a(2,2)*r93_y(5,8,x) . + a(2,3)*r93_y(6,8,x) elseif (name.eq.'1G4') then vpot(1,1) = a(2,4)*r93_y(3,8,x) elseif (name.eq.'3P0') then vpot(1,1) = a(3,1)*r93_y(3,8,x) . + a(3,2)*r93_y(5,8,x) . + a(2,5)*r93_z(3,8,x)/r3 elseif (name.eq.'3P1') then vpot(1,1) = a(3,3)*r93_y(3,8,x) . + a(3,4)*r93_y(5,8,x) . + a(3,5)*r93_z(3,8,x)/r3 elseif (name.eq.'3F3') then vpot(1,1) = a(4,5)*r93_y(3,8,x) elseif (name.eq.'3C2' .or. name.eq.'3C4') then vc = a(4,1)*r93_y(3,8,x) . + a(4,2)*r93_y(4,8,x) . + a(4,3)*r93_y(5,8,x) . + a(4,4)*r93_y(6,8,x) vt = ( a(5,1)*r93_z(4,8,x) + a(5,2)*r93_z(6,8,x) )/r3 if (name.eq.'3C2') vls = a(5,3)*r93_w(3,8,x) + . a(5,4)*r93_w(5,8,x) if (name.eq.'3C4') vls = a(5,5)*r93_w(3,8,x) vpot(1,1) = vc + (rj-r1)*vls - r2*(rj-r1)/(r2*rj+r1)*vt vpot(2,2) = vc - (rj+r2)*vls - r2*(rj+r2)/(r2*rj+r1)*vt vpot(1,2) = r6*sqrt(rj*(rj+r1))/(r2*rj+r1)*vt vpot(2,1) = vpot(1,2) elseif (name.eq.'1P1') then vpot(1,1) = b(2,1)*r93_y(3,8,x) . + b(2,2)*r93_y(4,8,x) . + b(2,3)*r93_y(5,8,x) . + b(2,4)*r93_y(6,8,x) elseif (name.eq.'1F3') then vpot(1,1) = b(1,5)*r93_y(3,8,x) . + b(2,5)*r93_y(5,8,x) elseif (name.eq.'3D2') then vpot(1,1) = b(3,1)*r93_y(3,8,x) . + b(3,2)*r93_y(5,8,x) . + b(3,3)*r93_z(3,8,x)/r3 elseif (name.eq.'3G4') then vpot(1,1) = b(3,4)*r93_y(3,8,x) elseif (name.eq.'3C1' .or. name.eq.'3C3') then vc = b(4,1)*r93_y(2,8,x) . + b(4,2)*r93_y(3,8,x) . + b(4,3)*r93_y(4,8,x) . + b(4,4)*r93_y(5,8,x) . + b(4,5)*r93_y(6,8,x) vt = ( b(3,5)*r93_z(4,8,x) + b(5,5)*r93_z(6,8,x) )/r3 if (name.eq.'3C1') vls = b(5,1)*r93_w(3,8,x) + . b(5,2)*r93_w(5,8,x) if (name.eq.'3C3') vls = b(5,3)*r93_w(3,8,x) + . b(5,4)*r93_w(5,8,x) vpot(1,1) = vc + (rj-r1)*vls - r2*(rj-r1)/(r2*rj+r1)*vt vpot(2,2) = vc - (rj+r2)*vls - r2*(rj+r2)/(r2*rj+r1)*vt vpot(1,2) = r6*sqrt(rj*(rj+r1))/(r2*rj+r1)*vt vpot(2,1) = vpot(1,2) elseif (spin.eq.0) then * j >= 5, spin = 0 if (iso.eq.1) then vpot(1,1) = a(1,1)*r93_y(2,8,x) . + a(1,2)*r93_y(3,8,x) . + a(1,3)*r93_y(4,8,x) . + a(1,4)*r93_y(5,8,x) . + a(1,5)*r93_y(6,8,x) elseif (iso.eq.0) then vpot(1,1) = b(2,1)*r93_y(3,8,x) . + b(2,2)*r93_y(4,8,x) . + b(2,3)*r93_y(5,8,x) . + b(2,4)*r93_y(6,8,x) endif elseif (spin.eq.1) then * j >= 5, spin = 1 if (iso.eq.1) then vc = a(4,1)*r93_y(3,8,x) . + a(4,2)*r93_y(4,8,x) . + a(4,3)*r93_y(5,8,x) . + a(4,4)*r93_y(6,8,x) vt = ( a(5,1)*r93_z(4,8,x) + a(5,2)*r93_z(6,8,x) )/r3 elseif (iso.eq.0) then vc = b(4,1)*r93_y(2,8,x) . + b(4,2)*r93_y(3,8,x) . + b(4,3)*r93_y(4,8,x) . + b(4,4)*r93_y(5,8,x) . + b(4,5)*r93_y(6,8,x) vt = ( b(3,5)*r93_z(4,8,x) + b(5,5)*r93_z(6,8,x) )/r3 endif if (nchan.eq.1) then if (l.eq.(j-1)) vpot(1,1) = vc - r2*(rj-r1)/(r2*rj+r1)*vt if (l.eq. j) vpot(1,1) = vc + r2*vt if (l.eq.(j+1)) vpot(1,1) = vc - r2*(rj+r2)/(r2*rj+r1)*vt elseif (nchan.eq.2) then vpot(1,1) = vc - r2*(rj-r1)/(r2*rj+r1)*vt vpot(2,2) = vc - r2*(rj+r2)/(r2*rj+r1)*vt vpot(1,2) = r6*sqrt(rj*(rj+r1))/(r2*rj+r1)*vt vpot(2,1) = vpot(1,2) endif endif vpot(1,1) = mpi * vpot(1,1) vpot(1,2) = mpi * vpot(1,2) vpot(2,1) = mpi * vpot(2,1) vpot(2,2) = mpi * vpot(2,2) * All together if (nchan.eq.1) then if (spin.eq.0) then if (l.eq.0) then vpot(1,1) = vpot(1,1) - r3*vspis else vpot(1,1) = vpot(1,1) - r3*vspi endif elseif (l.eq.j) then vpot(1,1) = vpot(1,1) + vspi + r2*vtpi elseif (name.eq.'3P0') then vpot(1,1) = vpot(1,1) + vspi - r2*(rj+r2)/(r2*rj+r1)*vtpi endif elseif (nchan.eq.2) then if (l.eq.0) then vpot(1,1) = vpot(1,1) + vspis - r2*(rj-r1)/(r2*rj+r1)*vtpi else vpot(1,1) = vpot(1,1) + vspi - r2*(rj-r1)/(r2*rj+r1)*vtpi endif vpot(2,2) = vpot(2,2) + vspi - r2*(rj+r2)/(r2*rj+r1)*vtpi vpot(1,2) = vpot(1,2) + r6*sqrt(rj*(rj+r1))/(r2*rj+r1)*vtpi vpot(2,1) = vpot(1,2) endif end ************************************************************************ function r93_y(in,im,x) implicit none integer in,im real*8 x,r93_y real*8 n,m,r1,r2 data r1 /1.0d0/ data r2 /2.0d0/ n = real(in) m = real(im) if (x < 1.0d-4) then r93_y = - n + m/r2 + n*n/(r2*m) else r93_y = exp(-n*x)/x - exp(-m*x)/x*(r1+(m*m-n*n)*x/(r2*m)) endif end ************************************************************************ function r93_yp(in,im,x) implicit none integer in,im real*8 x,r93_yp real*8 n,m,n2,m2,d,r1,r2 data r1 /1.0d0/ data r2 /2.0d0/ n = real(in) m = real(im) n2 = n*n m2 = m*m if (x < 1.0d-4) then d = m*(m2-n2)/(r2*n2) r93_yp = - n + m - d + x*(n2/r2 + m2/r2 + m*d) else r93_yp = exp(-n*x)/x - exp(-m*x)/x*(r1+(m2-n2)*m*x/(r2*n2)) endif end ************************************************************************ function r93_w(in,im,x) implicit none integer in,im real*8 x,r93_w real*8 n,m,n2,m2,x2,r1,r2,r3,r6,r8 data r1 /1.0d0/ data r2 /2.0d0/ data r3 /3.0d0/ data r6 /6.0d0/ data r8 /8.0d0/ n = real(in) m = real(im) n2 = n*n m2 = m*m if (x < 1.0d-4) then r93_w = (r2*n - r3*m + m2*m/n2)/r6 + . x*(r2*m2 - n2 - m2*m2/n2)/r8 else x2 = x*x r93_w = exp(-n*x)/x*(r1/(n*x)+r1/(n2*x2)) . - exp(-m*x)/x*(r1/(m*x)+r1/(m2*x2))*m2/n2 . - exp(-m*x)/x*(m2-n2)/(r2*n2) endif end ************************************************************************ function r93_z(in,im,x) implicit none integer in,im real*8 r93_z,x real*8 n,m,n2,m2,x2,r1,r2,r3,r4,r8 data r1 /1.0d0/ data r2 /2.0d0/ data r3 /3.0d0/ data r4 /4.0d0/ data r8 /8.0d0/ n = real(in) m = real(im) n2 = n*n m2 = m*m if (x < 1.0d-4) then r93_z = x*(n2 + r3*m2*m2/n2 - r4*m2)/r8 else x2 = x*x r93_z = exp(-n*x)/x*(r1+r3/(n*x)+r3/(n2*x2)) . - exp(-m*x)/x*(r1+r3/(m*x)+r3/(m2*x2))*m2/n2 . - exp(-m*x)*(r1+r1/(m*x))*m*(m2-n2)/(r2*n2) endif end ************************************************************************