program bder c from computations on displaced geometries make c Cartesian derivatives of polarizability implicit none integer*4 nat,i,n,ic,ne,ix,ia,ii,jj,is,ie real*8 step,f,bohr real*8,allocatable::a(:),al(:),gp(:),aa(:),ab(:) character*10 s10 bohr=0.529177d0 write(6,600) 600 format(' Number of atoms, step (A): ',$) read(5,*)nat,step step=step/bohr n=3*nat allocate(al(6*n),gp(9*n),aa(18*n),ab(27*n)) call vz(al, 6*n) call vz(gp, 9*n) call vz(aa,18*n) call vz(ab,27*n) do 200 ic=1,4 if(ic.eq.1)then ne= 6 write(6,*)' Reading electric dipole polarizabilities (alpha)' endif if(ic.eq.2)then ne= 9 write(6,*)' Reading the magnetic polarizability (Gp)' endif if(ic.eq.3)then write(6,*)' Reading quadrupole polarizability (A)' ne=18 endif if(ic.eq.4)then write(6,*)' Reading magnetic tensor (AB)' ne=27 endif allocate(a(ne)) ix=0 ia=1 f=0.5d0/step c i=1 .. zero geometry, skip here do 1 i=2,2*n+1 ix=ix+1 if(ix.eq.4)then ia=ia+1 ix=1 endif if(ia.gt.nat)then ia=1 ix=1 f=-f endif jj=ix+3*(ia-1) write(6,609)i,ia,ix,f 609 format(' Point ',i3,' atom ',i3,' coord ',i1,' f:',f4.0) write(s10,500)i 500 format(i10) do 2 is=1,len(s10) 2 if(s10(is:is).ne.' ')goto 3 3 if(i.gt.0)then if(ic.eq.1)then call readpol(a,s10(is:len(s10))//'/POL.TTT','POL',ie) if(ie.ne.0)then write(6,*)'skip this' goto 200 endif do 41 ii=1,ne 41 al(ii+ne*(jj-1))=al(ii+ne*(jj-1))+f*a(ii) endif if(ic.eq.2)then call readpol(a,s10(is:len(s10))//'/GP.TTT','GP',ie) if(ie.ne.0)then write(6,*)'skip this' goto 200 endif do 42 ii=1,ne 42 gp(ii+ne*(jj-1))=gp(ii+ne*(jj-1))+f*a(ii) endif if(ic.eq.3)then call readpol(a,s10(is:len(s10))//'/A.TTT','A',ie) if(ie.ne.0)then write(6,*)'skip this' goto 200 endif do 43 ii=1,ne 43 aa(ii+ne*(jj-1))=aa(ii+ne*(jj-1))+f*a(ii) endif if(ic.eq.4)then call readpol(a,s10(is:len(s10))//'/AB.TTT','AB',ie) if(ie.ne.0)then write(6,*)'skip this' goto 200 endif do 44 ii=1,ne 44 ab(ii+ne*(jj-1))=ab(ii+ne*(jj-1))+f*a(ii) endif endif 1 continue 200 deallocate(a) CALL WRITETTT09(NAT,al,gp,aa) CALL WRITETAB(NAT,ab) end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE vz(v,n) real*8 v(*) integer*4 n,i do 1 i=1,n 1 v(i)=0.0d0 return end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE WRITETTT09(NAT,Al,Gp,aa) IMPLICIT none INTEGER*4 NAT,I,ia,IX,J,K,jj real*8 al(*),Gp(*),aa(*) real*8,allocatable::l(:,:,:),g(:,:,:),a(:,:,:,:) OPEN(2,FILE='FILE.TTT') WRITE(2,2000)NAT,1,0.0d0 2000 FORMAT(' ROA tensors, cartesian derivatives',/, 1I4,' atoms, freq. ',i2,f11.6/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') allocate(l(3,3,3*NAT)) DO 11 ia=1,NAT DO 11 ix=1,3 jj=ix+3*(ia-1) l(1,1,jj)=al(1+6*(jj-1)) l(1,2,jj)=al(2+6*(jj-1)) l(2,1,jj)=al(2+6*(jj-1)) l(2,2,jj)=al(3+6*(jj-1)) l(1,3,jj)=al(4+6*(jj-1)) l(3,1,jj)=al(4+6*(jj-1)) l(2,3,jj)=al(5+6*(jj-1)) l(3,2,jj)=al(5+6*(jj-1)) 11 l(3,3,jj)=al(6+6*(jj-1)) DO 1 I=1,3 WRITE(2,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 ia=1,NAT DO 1 ix=1,3 jj=ix+3*(ia-1) 1 WRITE(2,2001)ia,ix,(L(I,J,jj),J=1,3) 2001 FORMAT(I5,1H ,I1,3g15.7) WRITE(2,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx(Bx) jy(By) jz(Bz)') allocate(g(3,3,3*NAT)) DO 21 ia=1,NAT DO 21 ix=1,3 jj=ix+3*(ia-1) do 21 i=1,3 do 21 j=1,3 21 g(i,j,jj)=gp(i+3*(j-1)+9*(jj-1)) DO 2 I=1,3 WRITE(2,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 ia=1,NAT DO 2 ix=1,3 jj=ix+3*(ia-1) 2 WRITE(2,2001)ia,ix,(G(i,j,jj),j=1,3) WRITE(2,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' Atom/x kx ky kz') allocate(a(3,3,3,3*nat)) DO 31 ia=1,NAT DO 31 ix=1,3 jj=ix+3*(ia-1) do 31 i=1,3 a(i,1,1,jj)=aa(i+3*(1-1)+6*(jj-1)) a(i,1,2,jj)=aa(i+3*(2-1)+6*(jj-1)) a(i,2,1,jj)=aa(i+3*(2-1)+6*(jj-1)) a(i,2,2,jj)=aa(i+3*(3-1)+6*(jj-1)) a(i,1,3,jj)=aa(i+3*(4-1)+6*(jj-1)) a(i,3,1,jj)=aa(i+3*(4-1)+6*(jj-1)) a(i,2,3,jj)=aa(i+3*(5-1)+6*(jj-1)) a(i,3,2,jj)=aa(i+3*(5-1)+6*(jj-1)) 31 a(i,3,3,jj)=aa(i+3*(6-1)+6*(jj-1)) DO 3 I=1,3 DO 3 J=1,3 WRITE(2,2006)I,J 2006 FORMAT(' A(',I1,',',I1,',K):') DO 3 ia=1,NAT DO 3 ix=1,3 jj=ix+3*(ia-1) 3 WRITE(2,2007)ia,ix,(A(K,J,I,jj)*3.d0/2.d0,K=1,3),ia,ix,I,J 2007 FORMAT(I5,I2,3g15.7,' ',4i3) write(2,*) write(2,*)'dummy alpha v:' DO 4 I=1,3 WRITE(2,2002)I DO 4 ia=1,NAT DO 4 IX=1,3 jj=3*(ia-1)+IX 4 WRITE(2,2001)ia,IX,(L(I,J,jj),J=1,3) CLOSE(2) RETURN END c ============================================================ subroutine readpol(a,s,t,ie) implicit none real*8 a(*) character*(*) s,t integer*4 IX,IY,ie logical lex ie=0 inquire(file=s,exist=lex) if(.not.lex)then write(6,*)s,' does not exist' ie=1 return endif OPEN(90,FILE=s) if(t.eq.'A')then read(90,*,end=99,err=99) read(90,*,end=99,err=99) do 1 IX=1,3 1 read(90,*,end=99,err=99)(a(IX+3*(IY-1)),IY=1,6) else if(t.eq.'POL')then read(90,*,end=99,err=99) read(90,*,end=99,err=99) read(90,*,end=99,err=99)(a(IX),IX=1,6) else if(t.eq.'GP')then read(90,*,end=99,err=99) read(90,*,end=99,err=99) read(90,*,end=99,err=99)(a(IX),IX=1,9) else if(t.eq.'AB')then read(90,*,end=99,err=99) read(90,*,end=99,err=99) do 2 IX=1,3 read(90,*,end=99,err=99) 2 read(90,*,end=99,err=99)(a(IX+3*(IY-1)),IY=1,9) else close(90) call report('Unknown polarizability requested') endif endif endif endif close(90) return 99 ie=1 write(6,*)'error in reading' close(90) return end c ============================================================ subroutine report(s) character*(*) s write(6,*)s stop end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE WRITETAB(NAT,ab) IMPLICIT none INTEGER*4 NAT,I,ia,IX,J,K,jj real*8 ab(*) real*8,allocatable::l(:,:,:,:) OPEN(2,FILE='FILE.TTB') WRITE(2,2000)NAT,1,0.0d0 2000 FORMAT( 1' Magnetically perturbed polarizability Cartesian derivatives',/, 1I4,' atoms, freq. ',i2,f11.6/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') allocate(l(3,3,3,3*NAT)) DO 11 ia=1,NAT DO 11 ix=1,3 jj=ix+3*(ia-1) DO 11 I=1,3 DO 11 J=1,3 DO 11 K=1,3 11 l(I,J,K,jj)=ab(I+3*(J-1)+9*(K-1)+27*(jj-1)) DO 1 I=1,3 WRITE(2,2002)I 2002 FORMAT(' Field ',I1,':') DO 1 ia=1,NAT DO 1 ix=1,3 jj=ix+3*(ia-1) DO 1 J=1,3 1 WRITE(2,2001)ia,ix,(L(I,J,K,jj),K=1,3) 2001 FORMAT(I5,1H ,I1,3g15.7) CLOSE(2) RETURN END c ============================================================