program rv implicit none integer*4 ND,n,i,n3, 1j,nm,nat,nd0,ne1,ne2,nv1,nv2,nn,ic,Jb,Kb,je,jv,ke, 2kv,KJ,jq parameter (nd0=21) real*8 CM,ecm,pi,po,RS,DS,CLIGHT,DEBYE,ux,uy,uz,mx,my,mz,pis real*8,allocatable::C(:,:),s(:,:), 1E(:),u(:,:,:),JdK(:),dd(:),ddq(:), 2dl(:,:),dv(:,:),g(:,:),qv(:,:,:),p(:,:,:),a(:,:,:), 3pq(:,:),aq(:,:),ew(:),m(:,:,:),mn(:,:,:),un(:,:,:), 4t(:,:) integer*4 ,allocatable::q(:) CM=219474.63d0 CLIGHT=137.035999173d0 DEBYE=2.54158d0 pi=4.0d0*atan(1.0d0) po=dsqrt(2.0d0*pi) open(9,file='ENBO.SCR',form='unformatted') read(9)ND close(9) allocate(E(ND),C(ND,ND)) call rre('ENBO.SCR',E,ND) call rre2('CNBO.SCR',C,ND) write(6,*)'Energies read from ENBO.SCR' write(6,*)'Eigenvectors read from CNBO.SCR' write(6,601)'CG 1 :',C(1 ,1 ) write(6,601)'C133 133:',C(133,133) write(6,601)'C133 134:',C(134,133) write(6,601)'C134 134:',C(134,134) write(6,601)'C134 133:',C(133,134) 601 format(1x,a9,1x,g11.5) allocate(s(1,1),q(1),ew(1)) call reads(s,1,nat,nm,q,ew,0) deallocate(s,ew,q) n3=3*nat allocate(ew(nm),q(nat),s(n3,n3)) call reads(s,n3,nat,nm,q,ew,1) n=10 ne1=1 ne2=3 nv1=1 nv2=133 if(nv2.eq.0)nv2=nm nn=ne2-ne1+1 allocate(JdK((nn+1)**2*nd0)) call vz(JdK,(nn+1)**2*nd0) allocate(dl(1,1),dv(1,1),g(1,1),qv(1,1,1)) dl(1,1)=0.0d0 dv(1,1)=0.0d0 g(1,1)=0.0d0 qv(1,1,1)=0.0d0 call rwdip(-1,n,dl,dv,g,qv) deallocate(dl,dv,g,qv) allocate(dl(n,3),dv(n,3),g(n,3),qv(n,3,3)) call rwdip(1,n,dl,dv,g,qv) i=0 nn=ne2-ne1+1 do 6 j=1,nn JdK( 1+nd0*i+nd0*(nn+1)*j)=dl(j,1) JdK( 2+nd0*i+nd0*(nn+1)*j)=dl(j,2) JdK( 3+nd0*i+nd0*(nn+1)*j)=dl(j,3) JdK( 4+nd0*i+nd0*(nn+1)*j)=dv(j,1) JdK( 5+nd0*i+nd0*(nn+1)*j)=dv(j,2) JdK( 6+nd0*i+nd0*(nn+1)*j)=dv(j,3) JdK( 7+nd0*i+nd0*(nn+1)*j)= g(j,1) JdK( 8+nd0*i+nd0*(nn+1)*j)= g(j,2) JdK( 9+nd0*i+nd0*(nn+1)*j)= g(j,3) JdK(13+nd0*i+nd0*(nn+1)*j)=qv(j,1,1) JdK(14+nd0*i+nd0*(nn+1)*j)=qv(j,1,2) JdK(15+nd0*i+nd0*(nn+1)*j)=qv(j,1,3) JdK(16+nd0*i+nd0*(nn+1)*j)=qv(j,2,2) JdK(17+nd0*i+nd0*(nn+1)*j)=qv(j,2,3) 6 JdK(18+nd0*i+nd0*(nn+1)*j)=qv(j,3,3) deallocate(dl,dv,g,qv) write(6,*)' <0|u|J> filled' write(6,*)' loading derivative <0|u|J> elements' allocate(dd(1),ddq(1)) dd(1)=0.0d0 ddq(1)=0.0d0 call rwdd(-1,n,nat,nd0,dd) n3=3*nat deallocate(dd,ddq) allocate(dd(nd0*n*n3),ddq(nd0*n*nm)) call rwdd( 1,n,nat,nd0,dd) call tv3(dd,ddq,n,nm,nd0,n3,s) write(6,*)' <0|u|J> derivatives transformed to normal modes' allocate(p(nat,3,3),a(nat,3,3),pq(nm,3),aq(nm,3)) call raptaat(nat,p,a) call taptaat(nat,p,a,pq,aq,s,nm) allocate(u(3,ND,ND),m(3,ND,ND)) call vz(u,3*ND**2) call vz(m,3*ND**2) c dipole elements in the BO basis: do 1 ke=ne1-1,ne2 do 1 kv=nv1-1,nv2 Kb=(nv2-nv1+2)*(ke-ne1+1)+kv-nv1+2 do 1 je=ne1-1,ne2 KJ=nd0*(ke-ne1+1)+nd0*(nn+1)*(je-ne1+1) do 1 jv=nv1-1,nv2 Jb=(nv2-nv1+2)*(je-ne1+1)+jv-nv1+2 if(ke.eq.ne1-1)then c |jv>: c du_KJ du_KJ c |jv> = u_KJ del(kv,jv) + ----- + ----- c dQm dPm if(je.eq.ne1-1)then c <0|0> if(kv.eq.nv1-1)then c <0|<0|u|J>|jv>: if(jv.gt.nv1-1)then c <0|<0|u|0>|jv=1>: if(CM*ew(jv).gt.1.0d0)then c <0|Q|1>=1/sqrt(2w) u(1,Kb,Jb)=u(1,Kb,Jb)+pq(jv,1)/sqrt(2.0d0*ew(jv)) u(2,Kb,Jb)=u(2,Kb,Jb)+pq(jv,2)/sqrt(2.0d0*ew(jv)) u(3,Kb,Jb)=u(3,Kb,Jb)+pq(jv,3)/sqrt(2.0d0*ew(jv)) c <0|P|1>=sqrt(w/2) pis=sqrt(ew(jv)/2.0d0) * 2.0d0/CLIGHT m(1,Kb,Jb)=m(1,Kb,Jb)+aq(jv,1)*pis m(2,Kb,Jb)=m(2,Kb,Jb)+aq(jv,2)*pis m(3,Kb,Jb)=m(3,Kb,Jb)+aq(jv,3)*pis call symd(ND,u,m,Kb,Jb) endif endif endif else c <0|J> if(kv.eq.jv)then u(1,Kb,Jb)=u(1,Kb,Jb)+JdK(1+KJ) u(2,Kb,Jb)=u(2,Kb,Jb)+JdK(2+KJ) u(3,Kb,Jb)=u(3,Kb,Jb)+JdK(3+KJ) m(1,Kb,Jb)=m(1,Kb,Jb)-JdK(7+KJ)/CLIGHT m(2,Kb,Jb)=m(2,Kb,Jb)-JdK(8+KJ)/CLIGHT m(3,Kb,Jb)=m(3,Kb,Jb)-JdK(9+KJ)/CLIGHT call symd(ND,u,m,Kb,Jb) endif if(kv.eq.nv1-1)then c <0|<0|u|J>|jv=1>: if(jv.gt.nv1-1)then if(CM*ew(jv).gt.1.0d0)then jq=nd0*(je-1)+nd0*n*(jv-1) u(1,Kb,Jb)=u(1,Kb,Jb)+ddq(1+jq)/sqrt(2.0d0*ew(jv)) u(2,Kb,Jb)=u(2,Kb,Jb)+ddq(2+jq)/sqrt(2.0d0*ew(jv)) u(3,Kb,Jb)=u(3,Kb,Jb)+ddq(3+jq)/sqrt(2.0d0*ew(jv)) pis=sqrt(ew(jv)/2.0d0) * 2.0d0/CLIGHT m(1,Kb,Jb)=m(1,Kb,Jb)+ddq(7+jq)*pis m(2,Kb,Jb)=m(2,Kb,Jb)+ddq(8+jq)*pis m(3,Kb,Jb)=m(3,Kb,Jb)+ddq(9+jq)*pis endif call symd(ND,u,m,Kb,Jb) endif endif endif else c if(ke.eq.ne1-1)then if(je.eq.ne1-1)then c else if(je.eq.ke)then c else c endif endif endif 1 continue write(6,602)'',(u(i,1,133),i=1,3) write(6,602)'',(m(i,1,133),i=1,3) write(6,601)'ds :', 1u(1,1,133)*u(1,1,133)+u(2,1,133)*u(2,1,133)+u(3,1,133)*u(3,1,133) write(6,601)'rs :', 1u(1,1,133)*m(1,1,133)+u(2,1,133)*m(2,1,133)+u(3,1,133)*m(3,1,133) write(6,602)'',(u(i,1,134),i=1,3) write(6,602)'',(m(i,1,134),i=1,3) write(6,601)'ds :', 1u(1,1,134)*u(1,1,134)+u(2,1,134)*u(2,1,134)+u(3,1,134)*u(3,1,134) write(6,601)'rs :', 1u(1,1,134)*m(1,1,134)+u(2,1,134)*m(2,1,134)+u(3,1,134)*m(3,1,134) 602 format(1x,a14,1x,3g11.5) write(6,*)'BO basis elements assigned' c transform to NBO wavefunction: allocate(un(3,ND,ND),mn(3,ND,ND),t(ND,ND)) do 3 ic=1,3 do 21 i=1,ND do 21 j=1,ND 21 t(i,j)=u(ic,i,j) call tunb(ND,t,C) do 22 i=1,ND do 22 j=1,ND 22 un(ic,i,j)=t(i,j) do 23 i=1,ND do 23 j=1,ND 23 t(i,j)=m(ic,i,j) call tunb(ND,t,C) do 3 i=1,ND do 3 j=1,ND 3 mn(ic,i,j)=t(i,j) end subroutine rre(f,a,n) implicit none integer*4 n,i character*(*) f real*8 a(*) open(9,file=f,form='unformatted') read(9)n read(9)(a(i),i=1,n) close(9) return end subroutine rre2(f,a,n) implicit none integer*4 i,j,n character*(*) f real*8 a(n,n) open(9,file=f,form='unformatted') read(9)n read(9)((a(i,j),i=1,n),j=1,n) close(9) return end subroutine reads(s,n30,nat,nm,q,e,ic) implicit none integer*4 n30,nm,nat,i,q(*),ia,k,ir,ic,n3 real*8 s(n30,n30),SFACR,e(*),CM SFACR=0.02342179d0 CM=219474.0d0 open(9,file='F.INP',status='old') read(9,*)nm,nat,nat if(ic.eq.0)goto 99 n3=3*nat do 8 i=1,nat 8 read(9,*)q(i) read(9,*) do 3 ia=1,nat do 3 i=1,nm ir=nm-i+1 3 read(9,*)(s(3*(ia-1)+k,ir),k=1,2),(s(3*(ia-1)+k,ir),k=1,3) read(9,*) read(9,*)(e(nm-i+1),i=1,nm) do 4 i=1,nm e(i)=e(i)/CM do 4 k=1,n3 4 S(k,i)=S(k,i)*SFACR write(6,*)'F.INP read, changed to au' 99 close(9) return end subroutine rwdip(ic,n,dl,dv,m,qv) implicit none integer*4 n,l,i,j,ic real*8 dl(n,3),dv(n,3),m(n,3),qv(n,3,3) open(9,file='DIPOLES.SCR.TXT') if(ic.eq.0)then write(9,*)n,' n' do 903 l=1,n 903 write(9,9023)l,(dl(l,i),i=1,3),(dv(l,i),i=1,3), 1 (m(l,i),i=1,3),((qv(l,i,j),i=1,3),j=1,3) 9023 format(i4,18f20.10) write(6,*)'DIPOLES.SCR.TXT written' else read(9,*)n if(ic.eq.-1)goto 99 do 904 l=1,n 904 read(9,*)dl(l,1),(dl(l,i),i=1,3),(dv(l,i),i=1,3), 1 (m(l,i),i=1,3),((qv(l,i,j),i=1,3),j=1,3) write(6,*)'DIPOLES.SCR.TXT read' endif 99 close(9) return end subroutine rwdd(ic,n,nat,nd0,dd) c ic=-1:read dimensions only c 0:write c 1:read implicit none integer*4 ic,n,nat,l,ia,ix,j,nd0,i real*8 dd(*) open(9,file='DIPOLEDER1.SCR.TXT') if(ic.eq.0)then write(9,*)'First dipole derivatives:' write(9,900)n,nat,nd0 900 format(3i5,' (n nat nd0)') do 904 l=1,n do 904 ia=1,nat do 904 ix=1,3 i=ix+3*(ia-1) write(9,90241)l,ia,ix 90241 format(2i4,i2,$) do 902 j=1,nd0 902 write(9,90242)dd(j+nd0*(l-1)+nd0*n*(i-1)) 90242 format(f20.10,$) 904 write(9,*) write(6,*)' Dipole derivatives written to DIPOLEDER1.SCR.TXT' else read(9,*) read(9,*)n,nat if(ic.eq.-1)then goto 99 else do 901 l=1,n do 901 ia=1,nat do 901 ix=1,3 i=ix+3*(ia-1) 901 read(9,*)(dd(j+nd0*(l-1)+nd0*n*(i-1)),j=1,3), 1 (dd(j+nd0*(l-1)+nd0*n*(i-1)),j=1,nd0) write(6,*)' Dipole derivatives read from DIPOLEDER1.SCR.TXT' endif endif 99 close(9) return end subroutine tv3(dd,ddq,n,nm,nd0,n3,s) implicit none integer*4 n,nm,nd0,n3,l,iq,jq,jj,i,j real*8 dd(*),ddq(*),sji,s(n3,n3) do 11 j=1,nm*nd0*n 11 ddq(j)=0.0d0 do 1 l=1,n do 1 iq=1,nm jq=nd0*(l-1)+nd0*n*(iq-1) do 1 i=1,n3 j =nd0*(l-1)+nd0*n*( i-1) sji=s(i,iq) do 1 jj=1,nd0 1 ddq(jj+jq)=ddq(jj+jq)+sji*dd(jj+j) return end subroutine raptaat(nat,p,a) implicit none integer*4 nat,L,J,I real*8 p(nat,3,3),a(nat,3,3) OPEN(15,FILE='FILE.TEN',STATUS='OLD') READ (15,*) C C READ DIPOLE DERIVATIVES = TENSOR P(A,B) C DO 10 L=1,nat DO 10 J=1,3 10 READ (15,*) (p(L,J,I),I=1,3) C C READ TENSORS I(A,B) AND J(A,B) FOR THE 0 0 0 ORIGIN: C (=ELECTRONIC & NUCLEAR TERMS) C DO 220 L=1,nat DO 220 J=1,3 220 READ (15,*) (a(L,J,I),I=1,3) C CLOSE(15) WRITE(*,*)' APT AAT read in' return end subroutine symd(ND,u,m,Kb,Jb) implicit none integer*4 ND,Kb,Jb real*8 u(3,ND,ND),m(3,ND,ND) u(1,Jb,Kb)= u(1,Kb,Jb) u(2,Jb,Kb)= u(2,Kb,Jb) u(3,Jb,Kb)= u(3,Kb,Jb) m(1,Jb,Kb)=-m(1,Kb,Jb) m(2,Jb,Kb)=-m(2,Kb,Jb) m(3,Jb,Kb)=-m(3,Kb,Jb) return end subroutine tunb(ND,t,C) implicit none integer*4 ND,i,ii,j,jj real*8 t(ND,ND),C(ND,ND),y real*8,allocatable::w(:,:) allocate(w(ND,ND)) do 1 i=1,ND do 1 jj=1,ND y=0.0d0 do 11 ii=1,ND 11 y=y+C(ii,i)*t(ii,jj) 1 w(i,jj)=y do 2 i=1,ND do 2 j=1,ND y=0.0d0 do 22 jj=1,ND 22 y=y+C(jj,j)*w(i,jj) 2 t(i,j)=y return end subroutine vz(x,n) implicit none integer*4 n,i real*8 x(*) do 1 i=1,n 1 x(i)=0.0d0 return end subroutine taptaat(nat,p,a,pq,aq,s,nm) implicit none integer*4 nat,I,nm,j,B,ii,l real*8 p(nat,3,3),a(nat,3,3),s(3*nat,3*nat), 1pq(nm,3),aq(nm,3) DO 56 I=1,nm DO 56 B=1,3 pq(I,B)=0.0d0 aq(I,B)=0.0d0 DO 56 l=1,nat DO 56 j=1,3 ii=j+3*(l-1) pq(I,B)=pq(I,B)+p(l,j,B)*s(ii,I) 56 aq(I,B)=aq(I,B)+a(l,j,B)*s(ii,I) return end