******************************************************************************** * * * 2DHF version 1-2003 * * Copyright (C) 1996 Jacek Kobus, Leif Laaksonen, Dage Sundholm * * * * This software may be used and distributed according to the terms * * of the GNU General Public License, see README and COPYING. * * * ******************************************************************************** c ### exchmom ### c c Calculates multipole moments to k_max=7 and Delta_m=4 c subroutine exchmom (iorb1,iorb2,psi,f4,wgt2,d0,d1,d2,d3, & d4,d5,d6,d7,wk1,wk2) implicit integer*4 (i-n) implicit real*8 (a-h,o-z) include 'commons8.inc' dimension psi(*),f4(*),wgt2(*),wk1(*),wk2(*), & d0(nni,mxnmu),d1(nni,mxnmu),d2(nni,mxnmu), & d3(nni,mxnmu),d4(nni,mxnmu), & d5(nni,mxnmu),d6(nni,mxnmu),d7(nni,mxnmu) dimension dome(10) if (itouch(iorb1).eq.0.and.itouch(iorb2).eq.0) return if (iorb1.eq.iorb2.and.mgx(6,iorb2).eq.0) return ido=0 1234 ido=ido+1 c calculate m idel=iabs(mgx(6,iorb2)-mgx(6,iorb1)) if (iorb1.eq.iorb2) idel=2*mgx(6,iorb2) ipc=iorb1+iorb2*(iorb2-1)/2 if ((iorb1.eq.iorb2).and.(ilc(ipc).lt.1)) return c second values of expansion coefficients for the same pair of c orbitals are stored in the second half of the arrays. if (ido.eq.2) then idel=mgx(6,iorb2)+mgx(6,iorb1) ipc=ipc+norb*(norb+1)/2 endif c homonuclear case is treated as a heteronuclear one ibeg1=i1b (iorb1) ibeg2=i1b (iorb2) xr2=r/2.d0 do mu=2,mxnmu do ni=2,nni rr=sqrt (vxisq(mu)+vetasq(ni)-1.d0) xrr=xr2*rr do i=1,8 dome(i)=0.d0 enddo call mulex(ni,mu,idel,dome) xw=xrr d0(ni,mu)=dome(1)*xw xw=xw*xrr d1(ni,mu)=dome(2)*xw xw=xw*xrr d2(ni,mu)=dome(3)*xw xw=xw*xrr d3(ni,mu)=dome(4)*xw xw=xw*xrr d4(ni,mu)=dome(5)*xw xw=xw*xrr d5(ni,mu)=dome(6)*xw xw=xw*xrr d6(ni,mu)=dome(7)*xw xw=xw*xrr d7(ni,mu)=dome(8)*xw c if ((iorb1.eq.1).and.(iorb2.eq.1).and.(ni.eq.nni/2) c & .and.(mu.eq.mxnmu/2)) then c write(*,'(10d16.7)') xrr,xw c write(*,'(10d16.7)') (dome(i),i=1,8) c endif c xw=xw*xrr c d8(ni,mu)=dome(9)*xw c xw=xw*xrr c d9(ni,mu)=dome(10)*xw enddo enddo c Now d0, d1, d2,... are equal to P0, P1*r, P2*r**2, ..., c respectively. P0, P1,.. are Legendre polynomials of a specified c order. c Prepare integrands and calculate moments. Multiplication by f4 is c necessary because of the factor 2/(r*cosh(mu)) incorporated in wgt2. ngrid = min (i1si(iorb1),i1si(iorb2)) call prod2 (ngrid,psi(ibeg1),psi(ibeg2),wk1) call prod (ngrid,f4,wk1) call prod2 (ngrid,d0,wk1,wk2) excdi(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d1,wk1,wk2) excqu(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d2,wk1,wk2) excoc(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d3,wk1,wk2) exche(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d4,wk1,wk2) exc5(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d5,wk1,wk2) exc6(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d6,wk1,wk2) exc7(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d7,wk1,wk2) exc8(ipc) = ddot (ngrid,wgt2,1,wk2,1) if (idbg(102).ne.0) then write(*,1000) iorn(iorb1),bond(iorb1),gut(iorb1), & iorn(iorb2),bond(iorb2),gut(iorb2),idel 1000 format(/,i4,1x,a8,a1,3x,i4,1x,a8,a1,3x,i5) write(*,1010) excdi(ipc),excqu(ipc),excoc(ipc),exche(ipc), & exc5(ipc),exc6(ipc),exc7(ipc),exc8(ipc) 1010 format(4d16.8) endif if (ilc(iorb1+iorb2*(iorb2-1)/2).eq.2.and.ido.eq.1) go to 1234 return end c ### exchmom1 ### c c Calculates multipole moments to l=8 and m=4 (iorb1<=iorb2) c subroutine exchmom1 (iorb1,iorb2,psi,f4,wgt2,d0,d1,d2,d3, & d4,d5,d6,d7,wk1,wk2) implicit integer*4 (i-n) implicit real*8 (a-h,o-z) include 'commons8.inc' dimension psi(*),f4(*),wgt2(*),wk1(*),wk2(*), & d0(nni,mxnmu),d1(nni,mxnmu),d2(nni,mxnmu), & d3(nni,mxnmu),d4(nni,mxnmu), & d5(nni,mxnmu),d6(nni,mxnmu),d7(nni,mxnmu) dimension dome(10) if (itouch(iorb1).eq.0.and.itouch(iorb2).eq.0) return if (iorb1.eq.iorb2.and.mgx(6,iorb2).eq.0) return ido=0 1234 ido=ido+1 c calculate m idel=iabs(mgx(6,iorb2)-mgx(6,iorb1)) if (iorb1.eq.iorb2) idel=2*mgx(6,iorb2) ipc=iorb1+iorb2*(iorb2-1)/2 if ((iorb1.eq.iorb2).and.(ilc(ipc).lt.1)) return c second values of expansion coefficients for the same pair of c orbitals are stored in the second half of the arrays. if (ido.eq.2) then idel=mgx(6,iorb2)+mgx(6,iorb1) ipc=ipc+norb*(norb+1)/2 endif c homonuclear case is treated as a heteronuclear one ibeg1=i1b (iorb1) ibeg2=i1b (iorb2) do mu=2,mxnmu do ni=2,nni xrr=r2*sqrt (vxisq(mu)+vetasq(ni)-1.d0) call mulex1(ni,mu,idel,dome) xw=xrr d0(ni,mu)=dome(1)*xw xw=xw*xrr d1(ni,mu)=dome(2)*xw xw=xw*xrr d2(ni,mu)=dome(3)*xw xw=xw*xrr d3(ni,mu)=dome(4)*xw xw=xw*xrr d4(ni,mu)=dome(5)*xw xw=xw*xrr d5(ni,mu)=dome(6)*xw xw=xw*xrr d6(ni,mu)=dome(7)*xw xw=xw*xrr d7(ni,mu)=dome(8)*xw c if ((iorb1.eq.1).and.(iorb2.eq.1).and.(ni.eq.nni/2) c & .and.(mu.eq.mxnmu/2)) then c write(*,'(10d16.7)') xrr,xw c write(*,'(10d16.7)') (dome(i),i=1,8) c endif c print *,'iorb1,iorb2,ni,mu,nnu,mxnmu', c & iorb1,iorb2,ni,mu,nni,mxnmu c stop "exchmom1" enddo enddo c Now d0, d1,... are equal to N*P(1,q)*r, N*P(2,q)*r**2, ..., c respectively. P(1,q), P(2,q),... are associate Legendre polynomials c of a specified order (k) and degree q (q==idel). N is equal to c [(k-|q|)!/(k+|q|)!]^{1/2}. c Note that N*N is present in eq. (19). One N is taken care of when c evaluating Q_kq (Q_{k\Delta m}^{ab}) and the other when generating c Pkm (P_{k|\Delta m|). c Prepare integrands and calculate moments. Multiplication by f4 is c necessary because of the factor 2/(r*cosh(mu)) incorporated in wgt2. ngrid = min (i1si(iorb1),i1si(iorb2)) call prod2 (ngrid,psi(ibeg1),psi(ibeg2),wk1) call prod (ngrid,f4,wk1) call prod2 (ngrid,d0,wk1,wk2) excdi(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d1,wk1,wk2) excqu(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d2,wk1,wk2) excoc(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d3,wk1,wk2) exche(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d4,wk1,wk2) exc5(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d5,wk1,wk2) exc6(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d6,wk1,wk2) exc7(ipc) = ddot (ngrid,wgt2,1,wk2,1) call prod2 (ngrid,d7,wk1,wk2) exc8(ipc) = ddot (ngrid,wgt2,1,wk2,1) if (idbg(102).ne.0) then write(*,1000) iorn(iorb1),bond(iorb1),gut(iorb1), & iorn(iorb2),bond(iorb2),gut(iorb2),idel 1000 format(/,i4,1x,a8,a1,3x,i4,1x,a8,a1,3x,i5) write(*,1010) excdi(ipc),excqu(ipc),excoc(ipc),exche(ipc), & exc5(ipc),exc6(ipc),exc7(ipc),exc8(ipc) 1010 format(4d16.8) endif if (ilc(iorb1+iorb2*(iorb2-1)/2).eq.2.and.ido.eq.1) go to 1234 return end