program rtt implicit integer*4 (i-n) implicit real*8 (a-h,o-z) parameter (N0=500,MX3=3*N0) DIMENSION ALPHA(3,3,MX3),GTENS(3,3,MX3),ATENS(3,3,3,MX3) character*1 ok character*2 at character*80 s80 logical num,lex dimension r(5000,3),iz(5000),p(N0,3,3),a(N0,3,3),f(3*N0,3*N0) c write(6,*)'Reads Turbomole tensors, makes FILE.X' write(6,*)' FILE.FC' write(6,*)' FILE.TEN' write(6,*) c bohr=0.529177d0 debye=2.541765d0 ia=0 c c coordinates XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX open(2,file='coord',status='old') c 1 read(2,222,end=1000,err=1000)ok 222 format(a1) if(ok.ne.'$')goto 1 c 2 read(2,*,err=1000,end=1000)x,y,z backspace 2 read(2,28)s80 28 format(a80) is=0 10 is=is+1 if(num(s80(is:is)))goto 10 at(1:1)=' ' at(2:2)=s80(is:is) ia=ia+1 r(ia,1)=x*bohr r(ia,2)=y*bohr r(ia,3)=z*bohr iq=0 if(at.eq.' h')iq=1 if(at.eq.' c')iq=6 if(at.eq.' n')iq=7 if(at.eq.' o')iq=8 if(at.eq.' k')iq=19 if(iq.eq.0)write(6,*)at,' unknown atom' iz(ia)=iq goto 2 1000 nat=ia close(2) write(6,*)nat,' atoms' if(nat.eq.0)stop open(3,file='FILE.X') write(3,*)' Turbomole coordinates' write(3,*)nat do 3 i=1,nat 3 write(3,3000)iz(i),(r(i,j),j=1,3) 3000 format(I3,3f15.8,' 0 0 0 0 0 0 0 0.0') close(3) write(6,*)'coord written into FILE.X' c c c TENSORS XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX inquire(file='aoforce.out', exist=lex) if(.not.lex)write(6,*)'aoforce.out does not exist !' if(lex)then open(3,file='aoforce.out') 100 read(3,28,end=999)s80 if(s80(13:50).eq.'cartesian derivatives of electronic po')then write(6,*)'Alpha derivatives found' do 6 i=1,nat*3 do 6 ix=1,3 do 6 jx=1,3 GTENS(ix,jx,i)=0.0d0 ALPHA(ix,jx,i)=0.0d0 do 6 kx=1,3 6 ATENS(ix,jx,kx,i)=0.0d0 read(3,*) do 22 i=1,6 if(i.eq.1)then ix=1 jx=1 endif if(i.eq.2)then ix=2 jx=2 endif if(i.eq.3)then ix=3 jx=3 endif if(i.eq.4)then ix=1 jx=2 endif if(i.eq.5)then ix=1 jx=3 endif if(i.eq.6)then ix=2 jx=3 endif do 9 j=1,4 9 read(3,*) n1=1 5 n2=min(n1+3,nat) read(3,*) read(3,*) do 4 iax=1,3 4 read(3,3001)(ALPHA(ix,jx,3*(ia-1)+iax),ia=n1,n2) 3001 format(13x,4g15.7) if(n2.lt.nat)then n1=n2+1 goto 5 endif 22 continue c write(6,*)ALPHA(3,3,3*(5-1)+3) do 8 i=1,nat*3 do 8 ix=1,3 do 8 jx=ix+1,3 8 ALPHA(jx,ix,i)=ALPHA(ix,jx,i) call WRITETTT(NAT,ALPHA,GTENS,ATENS,MX3) endif if(s80(9:44).eq.'cartesian dipole derivative matrix (')then write(6,*)'APT found' read(3,*) do 108 i=1,nat/2 read(3,*) read(3,*) read(3,*) read(3,*) do 101 k=1,3 101 read(3,3007)(p((i-1)*2+1,k,j),j=1,3),(p(i*2,k,j),j=1,3) 108 continue 3007 format(6x,6f11.6) if((nat/2)*2.ne.nat)then read(3,*) read(3,*) read(3,*) read(3,*) do 102 k=1,3 102 read(3,3007)(p(nat,k,j),j=1,3) endif write(6,*)'APT read' do 103 i=1,nat do 103 j=1,3 do 103 k=1,3 p(i,j,k)=p(i,j,k)/debye*bohr 103 a(i,j,k)=0.0d0 write(6,*)'Correct X-coordinate of APT (y/n)? ' read(5,'(a)')ok if(ok.eq.'y'.or.ok.eq.'Y')then do 403 i=1,nat do 403 j=1,3 403 p(i,j,1)=-p(i,j,1) write(6,667) 667 format(80(1h*),/,' Sign of Ux has been changed !!!',/,80(1h*)) endif call WRITETEN(N0,p,a,nat) WRITE(6,*)' APTs written into FILE.TEN' goto 999 endif goto 100 c 999 close(3) endif c c hessian XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX inquire(file='hessian', exist=lex) if(.not.lex)write(6,*)'hessian does not exist!' if(lex)then open(3,file='hessian') 105 read(3,28)s80 if(s80(1:5).eq.'level')then read(3,*)((((f(3*(ia-1)+ix,3*(ja-1)+jx),ix=1,3),jx=1,3), 1 ia=1,ja),ja=1,nat) do 88 j=1,3*nat do 88 i=1,j 88 f(j,i)=f(i,j) goto 888 endif if(s80(1:20).ne.'$hessian (projected)')goto 105 write(6,*)'projected hessian exists' do 104 i=1,3*nat js=1 106 je=min(js+4,3*nat) read(3,3008)(f(i,j),j=js,je) 3008 format(5x,5f15.10) if(je.lt.3*nat)then js=min(js+5,3*nat) goto 106 endif 104 continue 888 close(3) open(20,file='FILE.FC') call WRITEFF(3*N0,3*nat,f) close(20) write(6,*)'FF written into FILE.FC' endif c stop end c function num(l) logical num character*1 l num=.false. if(l.eq.'1')num=.true. if(l.eq.'2')num=.true. if(l.eq.'3')num=.true. if(l.eq.'4')num=.true. if(l.eq.'5')num=.true. if(l.eq.'6')num=.true. if(l.eq.'7')num=.true. if(l.eq.'8')num=.true. if(l.eq.'9')num=.true. if(l.eq.'0')num=.true. if(l.eq.'+')num=.true. if(l.eq.'-')num=.true. if(l.eq.'.')num=.true. if(l.eq.' ')num=.true. return end SUBROUTINE WRITEFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION FCAR(MX3,N) C CONST=4.359828/0.5291772**2 CONST=1.0d0 DO 6 I=1,N DO 6 J=1,N 6 FCAR(I,J)=FCAR(I,J)/CONST N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) RETURN END SUBROUTINE WRITETEN(N0,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION P(N0,3,3),A(N0,3,3) Z0=0.0d0 OPEN(15,FILE='FILE.TEN') WRITE(15,1500) NAT,NAT-6,0 1500 FORMAT(3I5) DO 10 L=1,NAT DO 10 J=1,3 10 WRITE(15,1501) (P(L,J,I),I=1,3),L 1501 FORMAT(3F14.8,I5) DO 221 L=1,NAT DO 221 J=1,3 221 WRITE(15,1501) (A(L,J,I),I=1,3),L DO 230 L=1,NAT DO 230 J=1,3 230 WRITE(15,1501) (Z0,I=1,3),L DO 100 L=1,NAT DO 100 J=1,3 100 WRITE(15,1501) (P(L,J,I),I=1,3),L CLOSE(15) RETURN END C SUBROUTINE WRITETTT(NAT,ALPHA,GTENS,ATENS,MX3) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3,3,MX3),GTENS(3,3,MX3),ATENS(3,3,3,MX3) OPEN(2,FILE='FILE.TTT') WRITE(2,2000)NAT 2000 FORMAT(' ROA tensors, cartesian derivatives',/,I4,' atoms',/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') DO 1 I=1,3 WRITE(2,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 WRITE(2,2001)L,IX,(ALPHA(I,J,IIND),J=1,3) 2001 FORMAT(I5,1H ,I1,3F15.7) WRITE(2,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx(Bx) jy(By) jz(Bz)') DO 2 I=1,3 WRITE(2,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 L=1,NAT DO 2 K=1,3 2 WRITE(2,2001)L,K,(GTENS(I,J,3*(L-1)+K),J=1,3) WRITE(2,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' Atom/x kx ky kz') DO 3 I=1,3 DO 3 J=1,3 WRITE(2,2006)I,J 2006 FORMAT(' A(',I1,',',I1,',K):') DO 3 L=1,NAT DO 3 M=1,3 3 WRITE(2,2007)L,M,(ATENS(K,J,I,3*(L-1)+M)*1.5d0,K=1,3),L,M,I,J 2007 FORMAT(I5,1H ,I1,3F15.7,' ',4i3) write(2,*) write(2,*)'dummy alpha v:' DO 4 I=1,3 WRITE(2,2002)I DO 4 L=1,NAT DO 4 IX=1,3 IIND=3*(L-1)+IX 4 WRITE(2,2001)L,IX,(ALPHA(I,J,IIND),J=1,3) CLOSE(2) write(6,*)'FILE.TTT written' RETURN END