******************************************************************************** * * * 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 ### propet #### c c Calculates dipole, quadrupole, octopole and hexadecapole moments c relative to the geometrical centre and the centre of mass. c c Last modification: 27/02/01 c subroutine propet (iscf,cw_orb,cw_coul,cw_exch, & f4,wgt2,wk0,wk1,wk2,wk3,wk4,wk5) implicit integer*4 (i-n) implicit real*8 (a-h,o-z) include 'commons8.inc' dimension cw_orb(*),cw_coul(*),cw_exch(*),f4(*), & wgt2(*),wk0(*),wk1(*),wk2(*),wk3(*),wk4(*),wk5(*) dimension elect(10),electDA(10),total(10),totalDA(10) c 1 bohr electron = 2.541765 Debye -- Gaussian94 User's Reference c data au2Debye /2.541765d0/ c data bohr2ang /0.529177249d0/ mpole=4 cm1=0.d0 cm2=0.d0 cm3=0.d0 cm4=0.d0 zcm=0.d0 qtot=0.d0 do iorb=1,norb qq=-occ(iorb) qtot=qtot+qq cm1=cm1+qq*cmulti(iorb) cm2=cm2+qq*cmulti(iorb+norb) cm3=cm3+qq*cmulti(iorb+2*norb) cm4=cm4+qq*cmulti(iorb+3*norb) enddo elect(1)=cm1 elect(2)=cm2 elect(3)=cm3 elect(4)=cm4 do k=1,mpole qkz1=z1*abs(r/2.d0)**k*plegendg(k,0,-1.d0) qkz2=z2*abs(r/2.d0)**k*plegendg(k,0, 1.d0) qktot1=cm1+qkz1+qkz2 total(k)=qktot1 totalDA(k)=total(k)*au2Debye*bohr2ang**(k-1) electDA(k)=elect(k)*au2Debye*bohr2ang**(k-1) enddo write(*,'(/," multipole moments relative to geometrical centre:")') write(*,'(35x,"dipole",14x,"quadrupole",11x,"octopole",12x, & "hexadecapole")') write(*,'(35x,"(Mu_z)",14x,"(Theta_zz)",10x,"(Omega_zzz)",11x, & "(Phi_zzzz)")') write(*,'(2x,"electronic (au) ",4e21.12)') (elect(k),k=1,mpole) write(*,'(2x,"electronic (Debye-Ang^k) ",4e21.12)') (electDA(k),k=1,mpole) write(*,'(2x,"total (au) ",4e21.12)') (total(k),k=1,mpole) write(*,'(2x,"total (Debye-Ang^k) ",4e21.12)') (totalDA(k),k=1,mpole) write(*,'(//," multipole moments relative to centre of mass:")') izz1=nint(z1) izz2=nint(z2) atw1=atweight(izz1) atw2=atweight(izz2) zcm=(-atw1+atw2)/(atw1+atw2)*r/2.d0 write(iout6,1000) element(izz1),atweight(izz1),-r/2.d0, & element(izz2), atweight(izz2),r/2.d0,zcm 01000 format(/,6x,'centre atomic weight z',/ & 8x, a2,4x,f14.8,3x,f10.4,/, & 8x, a2,4x,f14.8,3x,f10.4,/, & 2x,'centre-of-mass',15x,f10.4/) c multipole moments relative to the centre of mass do k=1,mpole c loop over grid points do imu=1,mxnmu inioff=(imu-1)*nni do in=1,nni igp=inioff+in c rr=(r/2.d0)*sqrt(vxisq(imu)+vetasq(in)-1.d0) c x=(r/2.d0)*vxi1(imu)*veta1(in) z=(r/2.d0)*vxi(imu)*veta(in)-zcm xxplusyy=((r/2.d0)*vxi1(imu)*veta1(in))**2 rr=sqrt(xxplusyy+z*z) if (abs(rr).lt.precis) then costh=0.d0 else costh=z/rr endif wk0(igp)=rr**k*plegendg(k,0,costh) enddo enddo qtot=0.d0 cm2zz=0.d0 do iorb=1,norb i1beg=i1b(iorb) ngorb=i1si(iorb) qq=-occ(iorb) qtot=qtot+qq call prod2 (ngorb,cw_orb(i1beg),cw_orb(i1beg),wk3) call prod (ngorb,wk0,wk3) call prod (ngorb,f4,wk3) w=ddot(ngorb,wgt2,1,wk3,1) cm2zz=cm2zz+qq*w enddo elect(k)=cm2zz electDA(k)=elect(k)*au2Debye*bohr2ang**(k-1) qkz1=z1*abs(r/2.d0+zcm)**k*plegendg(k,0,-1.d0) qkz2=z2*abs(r/2.d0-zcm)**k*plegendg(k,0, 1.d0) qktot=cm2zz+qkz1+qkz2 total(k)=qktot totalDA(k)=total(k)*au2Debye*bohr2ang**(k-1) enddo write(*,'(35x,"dipole",14x,"quadrupole",11x,"octopole",12x,"hexadecapole")') write(*,'(35x,"(Mu_z)",14x,"(Theta_zz)",10x,"(Omega_zzz)",11x,"(Phi_zzzz)")') write(*,'(2x,"electronic (au) ",4e21.12)') (elect(k),k=1,mpole) write(*,'(2x,"electronic (Debye-Ang^k) ",4e21.12)') (electDA(k),k=1,mpole) write(*,'(2x,"total (au) ",4e21.12)') (total(k),k=1,mpole) write(*,'(2x,"total (Debye-Ang^k) ",4e21.12)') (totalDA(k),k=1,mpole) if (idbg(221).ne.0) then write(*,*) write(*,*) write(*,*) 'Expectation values of r_1^k' do iorb=1,norb i1beg=i1b(iorb) ngorb=i1si(iorb) write(*,*) write(*,*) 'orbital #',iorb do mu=1,mxnmu imu=(mu-1)*nni do ni=1,nni r1=(r/2.d0)*(vxi(mu)+veta(ni)) if (r1.lt.precis) then wk0(imu+ni)=0.d0 else wk0(imu+ni)=1.d0/r1 endif wk1(imu+ni)=r1 wk2(imu+ni)=r1*r1 enddo enddo call prod2 (ngorb,cw_orb(i1beg),wk0,wk3) call prod (ngorb,cw_orb(i1beg),wk3) call prod (ngorb,f4,wk3) w=ddot(ngorb,wgt2,1,wk3,1) write(*,*) ' <1/r> =',w call prod2 (ngorb,cw_orb(i1beg),wk1,wk3) call prod (ngorb,cw_orb(i1beg),wk3) call prod (ngorb,f4,wk3) w=ddot(ngorb,wgt2,1,wk3,1) write(*,*) ' < r > =',w call prod2 (ngorb,cw_orb(i1beg),wk2,wk3) call prod (ngorb,cw_orb(i1beg),wk3) call prod (ngorb,f4,wk3) w=ddot(ngorb,wgt2,1,wk3,1) write(*,*) ' =',w enddo write(*,*) write(*,*) write(*,*) 'Expectation values of r^k' do iorb=1,norb i1beg=i1b(iorb) ngorb=i1si(iorb) write(*,*) write(*,*) 'orbital #',iorb do mu=1,mxnmu imu=(mu-1)*nni do ni=1,nni rr=sqrt(vxisq(mu)+vetasq(ni)-1.d0)*r/2.d0 if (rr.lt.precis) then wk0(imu+ni)=0.d0 else wk0(imu+ni)=1.d0/rr endif wk1(imu+ni)=rr wk2(imu+ni)=rr*rr enddo enddo call prod2 (ngorb,cw_orb(i1beg),wk0,wk3) call prod (ngorb,cw_orb(i1beg),wk3) call prod (ngorb,f4,wk3) w=ddot(ngorb,wgt2,1,wk3,1) write(*,*) ' <1/r> =',w call prod2 (ngorb,cw_orb(i1beg),wk1,wk3) call prod (ngorb,cw_orb(i1beg),wk3) call prod (ngorb,f4,wk3) w=ddot(ngorb,wgt2,1,wk3,1) write(*,*) ' < r > =',w call prod2 (ngorb,cw_orb(i1beg),wk2,wk3) call prod (ngorb,cw_orb(i1beg),wk3) call prod (ngorb,f4,wk3) w=ddot(ngorb,wgt2,1,wk3,1) write(*,*) ' =',w call prod2 (ngorb,cw_orb(i1beg),cw_orb(i1beg),wk3) call prod (ngorb,f4,wk3) w=ddot(ngorb,wgt2,1,wk3,1) write(*,*) ' < > =',w enddo write(*,*) write(*,*) write(*,*) 'Expectation values of r_2^k' do iorb=1,norb i1beg=i1b(iorb) ngorb=i1si(iorb) write(*,*) write(*,*) 'orbital #',iorb do mu=1,mxnmu imu=(mu-1)*nni do ni=1,nni r1=(r/2.d0)*(vxi(mu)-veta(ni)) if (r1.lt.precis) then wk0(imu+ni)=0.d0 else wk0(imu+ni)=1.d0/r1 endif wk1(imu+ni)=r1 wk2(imu+ni)=r1*r1 enddo enddo call prod2 (ngorb,cw_orb(i1beg),wk0,wk3) call prod (ngorb,cw_orb(i1beg),wk3) call prod (ngorb,f4,wk3) w=ddot(ngorb,wgt2,1,wk3,1) write(*,*) ' <1/r> =',w call prod2 (ngorb,cw_orb(i1beg),wk1,wk3) call prod (ngorb,cw_orb(i1beg),wk3) call prod (ngorb,f4,wk3) w=ddot(ngorb,wgt2,1,wk3,1) write(*,*) ' < r > =',w call prod2 (ngorb,cw_orb(i1beg),wk2,wk3) call prod (ngorb,cw_orb(i1beg),wk3) call prod (ngorb,f4,wk3) w=ddot(ngorb,wgt2,1,wk3,1) write(*,*) ' =',w enddo endif return end