PROGRAM intt implicit none INTEGER*4 MX,NAT,L,I,ip,jp,kp,k,ix,jx,ixp,jxp,i1,i2,j1,j2, 1ia,IIND,j,ja,IINDP,ixx PARAMETER (MX=1000) INTEGER*4 IBOND(4,MX),NB(MX),iii(4,MX),ND(MX) REAL*8 F(3*MX,3*MX),FI(3*MX,3*MX),R(3,MX),dxx,bohr, 1A(MX,3,3),P(MX,3,3),ALPHA(3*MX,3,3),AA(3*MX,3,3,3),G(3*MX,3,3), 1AI(MX,3,3),PI(MX,3,3),ALPHAI(3*MX,3,3),AAI(3*MX,3,3,3), 1GI(3*MX,3,3),v1(3),v2(3),v3(3),u(3,3,MX),pol(3,3),POLI(3,3), 1apol(3,3,3),gpol(3,3),coo(3),GPOLI(3,3),APOLI(3,3,3) logical lxfile,lff,lten,lttt,lpol,lpola,lpolg,lt character*80 fg,fa bohr=0.529177d0 c WRITE(6,3000) 3000 FORMAT(/, 1 ' Rewrites tensors into internal bond-based coordinates',/,/, 9 ' INPUT FILES: FILE.X coordinates',/, 1 ' Cartesian Tensors: FILE.FC force field',/, 1 ' FILE.TEN APT and AAT tensors',/, 3 ' FILE.TTT ROA tensors',/, 3 ' POL.TTT polarizability',/, 3 ' GP.TTT G polarizability',/, 3 ' A.TTT A polarizability',/,/, 1 ' OUTPUT FILES:',/, 1 ' Internal Tensors: INT.FC force field',/, 1 ' INT.TEN APT and AAT tensors',/, 3 ' INT.TTT ROA tensors',/, 3 ' INT.POL polarizability',/) c inquire(FILE='FILE.X',EXIST=lxfile) if(lxfile)then OPEN(2,FILE='FILE.X',STATUS='OLD') CALL READRX(MX,NAT,R,ND,IBOND,NB) CLOSE(2) else call report('FILE.X does not exist') endif WRITE(6,*)NAT, ' atoms in FILE.X' c INQUIRE(FILE='FILE.FC',EXIST=lff) IF(lff)THEN OPEN(2,FILE='FILE.FC',STATUS='OLD') CALL READFF(MX,NAT,F) CLOSE(2) WRITE(6,3001)NAT 3001 FORMAT(' Force field found,',I10,' atoms') ENDIF c INQUIRE(FILE='FILE.TEN',EXIST=lten) IF(lten)THEN write(6,*)'Tensors FILE.TEN found' OPEN(15,FILE='FILE.TEN') CALL READTEN(MX,P,A,.false.,R,NAT) CLOSE(15) ENDIF C INQUIRE(FILE='FILE.TTT',EXIST=lttt) IF(lttt)THEN write(6,*)'Tensors FILE.TTT found' OPEN(15,FILE='FILE.TTT') CALL READTTT(MX,ALPHA,AA,G,NAT) CLOSE(15) CALL TTTT(3*MX,ALPHA,AA,G,NAT,1,R) ENDIF C pol=0.0d0 apol=0.0d0 gpol=0.0d0 INQUIRE(FILE='POL.TTT',EXIST=lpol) INQUIRE(FILE='GP.TTT',EXIST=lpolg) if(lpolg)fg='GP.TTT' INQUIRE(FILE='POL.G.TTT',EXIST=lt) if(lt)fg='POL.G.TTT' lpolg=lpolg.or.lt INQUIRE(FILE='A.TTT',EXIST=lpola) if(lpola)fa='A.TTT' INQUIRE(FILE='POL.A.TTT',EXIST=lt) if(lt)fa='POL.A.TTT' lpola=lpola.or.lt IF(lpol)THEN OPEN(15,FILE='POL.TTT') call readpol(pol) CLOSE(15) write(6,*)'POL.TTT found' ENDIF IF(lpolg)call readpolg(gpol,fg) IF(lpola)call readpola(apol,fa) c transform G, A to molecule center if(lpola.or.lpolg)then do 301 ix=1,3 coo(ix)=0.0d0 do 3011 ia=1,NAT 3011 coo(ix)=coo(ix)+R(ix,ia) 301 coo(ix)=coo(ix)/bohr/dble(NAT) call TTT1(pol,apol,gpol,1,coo,1) endif if(.not.lttt.and..not.lten.and..not.lff) 1call report('Nothing to do') c c Define two coordinate vectors for each atom if(NAT.EQ.2)then do 101 ia=1,NAT if(NB(ia).lt.1)call report('Unbounded atom') i1=ia i2=IBOND(1,IA) do 102 ix=1,3 v2(ix)=0.0d0 102 v1(ix)=R(ix,i2)-R(ix,i1) dxx=dabs(v1(1)) ixx=1 if(dabs(v1(2)).lt.dxx)then ixx=2 dxx=dabs(v1(2)) endif if(dabs(v1(3)).lt.dxx)ixx=3 v2(ixx)=1.0d0 call make3v(v1,v2,v3) iii(1,ia)=i1 iii(2,ia)=i2 iii(3,ia)=0 iii(4,ia)=0 do 101 ix=1,3 u(1,ix,ia)=v1(ix) u(2,ix,ia)=v2(ix) 101 u(3,ix,ia)=v3(ix) else if(NAT.eq.1)then iii(1,ia)=0 iii(2,ia)=0 iii(3,ia)=0 iii(4,ia)=0 do 201 ix=1,3 u(1,ix,ia)=0.0d0 u(2,ix,ia)=0.0d0 201 u(3,ix,ia)=0.0d0 u(1,1,ia)=1.0d0 u(2,2,ia)=1.0d0 u(3,3,ia)=1.0d0 else c v1=i1->i2, v2=j1->j2 do 1 ia=1,NAT if(NB(ia).lt.1)call report('Unbounded atom') if(NB(ia).eq.1)then i1=ia i2=IBOND(1,IA) j1=i2 if(NB(i2).le.1)call report('Not enough bonds') if(IBOND(1,i2).EQ.ia)then j2=IBOND(2,i2) else j2=IBOND(1,i2) endif else i1=ia j1=ia i2=IBOND(1,IA) j2=IBOND(2,IA) endif do 2 ix=1,3 v1(ix)=R(ix,i2)-R(ix,i1) 2 v2(ix)=R(ix,j2)-R(ix,j1) call make3v(v1,v2,v3) iii(1,ia)=i1 iii(2,ia)=i2 iii(3,ia)=j1 iii(4,ia)=j2 do 1 ix=1,3 u(1,ix,ia)=v1(ix) u(2,ix,ia)=v2(ix) 1 u(3,ix,ia)=v3(ix) endif endif C IF(lff)THEN WRITE(6,*)' Transforming FF' do 3 ia=1,NAT do 3 ja=1,NAT do 3 ix=1,3 do 3 jx=1,3 FI(3*(ia-1)+ix,3*(ja-1)+jx)=0.0d0 do 3 ixp=1,3 do 3 jxp=1,3 3 FI(3*(ia-1)+ix,3*(ja-1)+jx)=FI(3*(ia-1)+ix,3*(ja-1)+jx) 1 +F(3*(ia-1)+ixp,3*(ja-1)+jxp)*u(ix,ixp,ia)*u(jx,jxp,ia) OPEN(2,FILE='INT.FC') CALL WRITEFF(MX,NAT,FI) call writeatms(NAT,iii,MX) CLOSE(2) WRITE(6,*)' File INT.FC written ' ENDIF IF(lten)then WRITE(6,*)' Transforming APT,AAT' DO 10 L=1,NAT DO 10 J=1,3 DO 10 I=1,3 PI(L,J,I)=0.0d0 AI(L,J,I)=0.0d0 DO 10 JP=1,3 DO 10 IP=1,3 AI(L,J,I)=AI(L,J,I)+A(L,JP,IP)*u(J,JP,L)*u(I,IP,L) 10 PI(L,J,I)=PI(L,J,I)+P(L,JP,IP)*u(J,JP,L)*u(I,IP,L) OPEN(2,FILE='INT.TEN') CALL WRITETEN(MX,PI,AI,NAT) call writeatms(NAT,iii,MX) CLOSE(2) WRITE(6,*)' File INT.TEN written ' endif IF(lttt)then WRITE(6,*)' Transforming ALPHA, GP, A' DO 11 L=1,NAT DO 11 IX=1,3 IIND=3*(L-1)+IX DO 11 I=1,3 DO 11 J=1,3 ALPHAI(IIND,I,J)=0.0d0 GI(IIND,I,J)=0.0d0 DO 11 IXP=1,3 IINDP=3*(L-1)+IXP DO 11 JP=1,3 DO 11 IP=1,3 ALPHAI(IIND,I,J)=ALPHAI(IIND,I,J) 1 +ALPHA(IINDP,IP,JP)*u(IX,IXP,L)*u(I,IP,L)*u(J,JP,L) 11 GI(IIND,I,J)=GI(IIND,I,J) 1 +G(IINDP,IP,JP)*u(IX,IXP,L)*u(I,IP,L)*u(J,JP,L) DO 13 L=1,NAT DO 13 IX=1,3 IIND=3*(L-1)+IX DO 13 I=1,3 DO 13 J=1,3 DO 13 K=1,3 AAI(IIND,I,J,K)=0.0d0 DO 13 IXP=1,3 IINDP=3*(L-1)+IXP DO 13 IP=1,3 DO 13 JP=1,3 DO 13 KP=1,3 13 AAI(IIND,I,J,K)=AAI(IIND,I,J,K) 1 +AA(IINDP,IP,JP,KP)*u(IX,IXP,L)*u(I,IP,L)*u(J,JP,L)*u(K,KP,L) OPEN(2,FILE='INT.TTT') CALL WRITETTT(MX,NAT,ALPHAI,AAI,GI) call writeatms(NAT,iii,MX) close(2) WRITE(6,*)' File INT.TTT written ' ENDIF IF(lpol)then WRITE(6,*)' Transforming polarizability with the first atom' call ptp(POLI,POL,u,1,MX) call ptp(GPOLI,GPOL,u,1,MX) call ptA(APOLI,APOL,u,1,MX) CALL writepol(POLI,'INT.POL',1,iii,MX) if(lpolg)CALL writepolg(GPOLI,'INTG.POL',1,iii,MX) if(lpola)CALL writepola(APOLI,'INTA.POL',1,iii,MX) ENDIF call report(' Program finished OK.') END subroutine ptp(POLI,POL,u,ia,MX) implicit none integer*4 MX,ia,i,j,ip,jp real*8 POLI(3,3),POL(3,3),u(3,3,MX),ut(3,3),t(3,3) do 1 i=1,3 do 1 j=1,3 1 ut(i,j)=u(i,j,ia) do 2 i=1,3 do 2 jp=1,3 t(i,jp)=0.0d0 do 2 ip=1,3 2 t(i,jp)=t(i,jp)+ut(i,ip)*POL(ip,jp) do 3 i=1,3 do 3 j=1,3 POLI(i,j)=0.0d0 do 3 jp=1,3 3 POLI(i,j)=POLI(i,j)+ut(j,jp)*t(i,jp) return end subroutine ptA(AI,A,u,ia,MX) implicit none integer*4 MX,ia,i,j,ip,jp,k,kp real*8 AI(3,3,3),A(3,3,3),u(3,3,MX),ut(3,3),t(3,3,3) do 1 i=1,3 do 1 j=1,3 1 ut(i,j)=u(i,j,ia) do 2 i=1,3 do 2 jp=1,3 do 2 kp=1,3 t(i,jp,kp)=0.0d0 do 2 ip=1,3 2 t(i,jp,kp)=t(i,jp,kp)+ut(i,ip)*A(ip,jp,kp) do 3 i=1,3 do 3 j=1,3 do 3 kp=1,3 AI(i,j,kp)=0.0d0 do 3 jp=1,3 3 AI(i,j,kp)=AI(i,j,kp)+ut(j,jp)*t(i,jp,kp) do 4 i=1,3 do 4 j=1,3 do 4 k=1,3 t(i,j,k)=0.0d0 do 4 kp=1,3 4 t(i,j,k)=t(i,j,k)+ut(k,kp)*t(i,j,kp) AI=t return end c SUBROUTINE READRX(N0,N,R,ND,IBONDS,NB) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,N0),ND(N0),IBONDS(4,N0),NB(N0) READ(2,*) read(2,*)N IF(N.GT.N0)CALL REPORT(' N > MX !') do 1 i=1,N READ(2,*)ND(i),R(1,i),R(2,i),R(3,i),(IBONDS(IB,i),IB=1,4) NB(i)=0 DO 1 IB=1,4 1 IF(IBONDS(IB,i).NE.0)NB(i)=IB RETURN END SUBROUTINE READFF(N0,NA,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(3*N0,3*N0) NAT3=3*NA N=NAT3 C Read in the lower triangle of FF, C written in parts as n1,n1 C . . C ln,n1 . ln,ln C . . . C . . . n3,n3 C . . . . C n,n1 . . n,n3 N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(2,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) C DO 31 I=1,N DO 31 J=I,N 31 FCAR(I,J)=FCAR(J,I) RETURN END SUBROUTINE WRITEFF(MX,NA,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(3*MX,3*MX) N=NA*3 N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 WRITE(2,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 C SUBROUTINE REPORT(RS) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) RS WRITE(6,3000) 3000 FORMAT(/,80(1H*)) WRITE(6,*)RS WRITE(6,3001) 3001 FORMAT(80(1H*),/,/,'PROGRAM STOPPED') STOP END c SUBROUTINE READTEN(N0,P,A,LAPT,X,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=1000) DIMENSION P(N0,3,3),A(N0,3,3),AI(MX,3,3),AJ(MX,3,3),PV(MX,3,3), 1R(3,MX),X(3,NAT) LOGICAL LAPT BOHR=0.52917705993d0 DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(I,L)/BOHR c READ (15,*) NOAT if(NOAT.ne.NAT)call report('NAT<>NOAT') C Read dipole derivatives: DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) C C If no axial tensor, then construct its polar part only: IF(LAPT)THEN DO 2203 L=1,NAT DO 2203 J=1,3 DO 2203 I=1,3 SUM=0.0d0 DO 2204 ID=1,3 DO 2204 IG=1,3 2204 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) AJ(L,J,I)=0.0d0 PV(L,J,I)=P(L,J,I) 2203 AI(L,J,I)=SUM WRITE(6,*)' APT part of AAT constructed' GOTO 1 ENDIF DO 220 L=1,NOAT DO 220 J=1,3 220 READ (15,*) (AI(L,J,I),I=1,3) DO 230 L=1,NOAT DO 230 J=1,3 230 READ (15,*) (AJ(L,J,I),I=1,3) DO 100 L=1,NOAT DO 100 J=1,3 100 READ (15,*) (PV(L,J,I),I=1,3) C Distributed-origin gauge: 1 CALL VCDD0(AI,PV,P,R,NAT) WRITE(6,700) 700 format(' DOG used.') C C Total axial tensor = electronic + nuclear: DO 4 L=1,NAT DO 4 I=1,3 DO 4 J=1,3 4 A(L,I,J)=AI(L,I,J)+AJ(L,I,J) C C Put axial tensor in the origin on the atom L: DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM-0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM WRITE(6,*)' AAT shifted to atomic origin' RETURN END C 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) WRITE(2,1500) NAT,NAT-6,0 1500 FORMAT(3I5) DO 10 L=1,NAT DO 10 J=1,3 10 WRITE(2,1501) (P(L,J,I),I=1,3),L 1501 FORMAT(3F14.8,I5) DO 220 L=1,NAT DO 220 J=1,3 220 WRITE(2,1501) (A(L,J,I),I=1,3),L Z0=0.0d0 DO 230 L=1,NAT DO 230 J=1,3 230 WRITE(2,1501) (Z0,I=1,3),L DO 100 L=1,NAT DO 100 J=1,3 100 WRITE(2,1501) (P(L,J,I),I=1,3),L RETURN END C FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K EPS=0.0d0 IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0d0 IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0d0 IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0d0 IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0d0 IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0d0 IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0d0 RETURN END C SUBROUTINE VCDD0(VCD,VEL,DIP,C,NAT) C STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) INTEGER*4 ALPHA,BETA,GAMMA,DELTA,LAMBDA PARAMETER (MX=1000) DIMENSION VCD(MX,3,3),VEL(MX,3,3),DIP(MX,3,3), 1C(3,NAT) DO 90 BETA =1,3 DO 90 GAMMA=1,3 DO 90 DELTA=1,3 SKEW=0.25D0*EPS(BETA,GAMMA,DELTA) IF (DABS(SKEW).GT.1.0D-10)THEN DO 80 ALPHA=1,3 DO 80 LAMBDA=1,NAT 80 VCD(LAMBDA,ALPHA,BETA)=VCD(LAMBDA,ALPHA,BETA) 1 -SKEW*C(GAMMA,LAMBDA) 2 *(VEL(LAMBDA,ALPHA,DELTA)-DIP(LAMBDA,ALPHA,DELTA)) ENDIF 90 CONTINUE RETURN END C SUBROUTINE READTTT(MX,ALPHA,A,G,NAT) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*MX,3,3),G(3*MX,3,3),A(3*MX,3,3,3) READ(15,*) READ(15,*)NAT READ(15,*) READ(15,*) DO 1 I=1,3 READ(15,*) DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 READ(15,*)(ALPHA(IIND,I,J),J=1,2),(ALPHA(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 2 I=1,3 READ(15,*) DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 READ(15,*)(G(IIND,I,J),J=1,2),(G(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 3 I=1,3 DO 3 J=1,3 READ(15,*) DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX 3 READ(15,*)(A(IIND,I,J,K),K=1,2),(A(IIND,I,J,K),K=1,3) CLOSE(15) RETURN END C SUBROUTINE WRITETTT(MX,NAT,ALPHA,A,G) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*MX,3,3),G(3*MX,3,3),A(3*MX,3,3,3) 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(IIND,I,J),J=1,3) 2001 FORMAT(I5,1H ,I1,3F13.7) WRITE(2,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx jy jz') DO 2 I=1,3 WRITE(2,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 WRITE(2,2001)L,IX,(G(IIND,I,J),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 IX=1,3 IIND=3*(L-1)+IX 3 WRITE(2,2001)L,IX,(A(IIND,I,J,K),K=1,3) RETURN END SUBROUTINE TTTT(N,ALPHA,A,G,NAT,ICO,R) C Transforms A and G to local origins (ICO=1) c or common origin (ICO=2); C sets A and G to zero for ICO=3. C R supposed to be passed in in A C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=1000) DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,NAT),X(3,MX) C C transfer into atomic coordinates: DO 5 I=1,3 DO 5 J=1,NAT 5 X(I,J)=R(I,J)/0.529177D0 C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 W=1.0d0 C supposedly the static limit G/w is calculated C DO 2 IA=1,3 DO 2 IB=1,3 DO 2 L=1,NAT DO 2 IE=1,3 IIND=3*(L-1)+IE SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(IIND,IA,ID) G(IIND,IA,IB)=G(IIND,IA,IB)+SUM IF(ICO.EQ.3)G(IIND,IA,IB)=0.0d0 2 CONTINUE IF(ICO.EQ.1)WRITE(6,*)' G transformed into atomic origins' IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' G set to zero' C DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 DO 3 L=1,NAT DO 3 IE=1,3 IIND=3*(L-1)+IE SUM=0.0d0 IF(IB.EQ.IC)THEN DO 4 ID=1,3 4 SUM=SUM+X(ID,L)*ALPHA(IIND,IA,ID)*SIGN ENDIF A(IIND,IA,IB,IC)=A(IIND,IA,IB,IC)+SUM 1-1.5d0*SIGN*(X(IB,L)*ALPHA(IIND,IA,IC)+X(IC,L)*ALPHA(IIND,IA,IB)) IF(ICO.EQ.3)A(IIND,IA,IB,IC)=0.0d0 3 CONTINUE IF(ICO.EQ.1)WRITE(6,*)' A transformed into atomic origins' IF(ICO.EQ.2)WRITE(6,*)' A transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' A set to zero' RETURN END c subroutine norm(v) implicit none real*8 v(*),n n=dsqrt(v(1)**2+v(2)**2+v(3)**2) v(1)=v(1)/n v(2)=v(2)/n v(3)=v(3)/n return end subroutine writeatms(NAT,iii,MX) integer*4 MX,NAT,iii(4,MX),ia do 4 ia=1,NAT 4 write(2,200)ia,(iii(ii,ia),ii=1,4) 200 format(5i5) return end subroutine readpol(a) real*8 a(3,3) read(15,*) read(15,*) read(15,*)((a(i,j),i=1,j),j=1,3) do 1 j=1,3 do 1 i=1,j 1 a(j,i)=a(i,j) return end subroutine writepol(a,fn,NAT,iii,MX) implicit none integer*4 MX,i,j,iii(4,MX),NAT character*(*) fn real*8 a(3,3) OPEN(2,FILE=fn) write(2,16) 16 format('polarizability in internal coordinates',/, 1' XX XY YY XZ YZ ZZ') write(2,17)((a(i,j),i=1,j),j=1,3) 17 format(6f9.3) call writeatms(NAT,iii,MX) close(2) WRITE(6,*)fn//' written' return end subroutine writepolg(a,fn,NAT,iii,MX) implicit none integer*4 MX,i,j,iii(4,MX),NAT character*(*) fn real*8 a(3,3) OPEN(2,FILE=fn) write(2,16) 16 format(' G-tens in internal coordinates',/, 1' XX YX ZX', 2' XY YY ZY', 2' XZ YZ ZZ') write(2,17)((a(j,i),i=1,3),j=1,3) 17 format(9f9.3) call writeatms(NAT,iii,MX) close(2) WRITE(6,*)fn//' written' return end subroutine writepola(a,fn,NAT,iii,MX) implicit none integer*4 MX,i,j,k,iii(4,MX),NAT character*(*) fn real*8 a(3,3,3) OPEN(2,FILE=fn) write(2,16) 16 format('polarizability in internal coordinates',/, 1' XX XY YY XZ YZ ZZ') do 1 k=1,3 1 write(2,17)((a(i,j,k),i=1,j),j=1,3) 17 format(6f9.3) call writeatms(NAT,iii,MX) close(2) WRITE(6,*)fn//' written' return end subroutine make3v(v1,v2,v3) c makes system of orthonormal vectors from v1 and v2 implicit none real*8 v1(*),v2(*),v3(*) v3(1)=v1(2)*v2(3)-v1(3)*v2(2) v3(2)=v1(3)*v2(1)-v1(1)*v2(3) v3(3)=v1(1)*v2(2)-v1(2)*v2(1) call norm(v1) call norm(v3) v2(1)=v3(2)*v1(3)-v3(3)*v1(2) v2(2)=v3(3)*v1(1)-v3(1)*v1(3) v2(3)=v3(1)*v1(2)-v3(2)*v1(1) call norm(v2) return end SUBROUTINE readpolg(a,fn) character*(*)fn real*8 a(3,3) integer*4 i,j open(90,file=fn) read(90,*) read(90,*) read(90,*)((a(j,i),i=1,3),j=1,3) close(90) write(6,*)fn//' read' return end c SUBROUTINE readpola(a,fn) character*(*)fn integer*4 i,j,k real*8 a(3,3,3) open(90,file=fn) read(90,*) read(90,*) do 1 k=1,3 read(90,*)((a(i,j,k),i=1,j),j=1,3) do 1 j=1,3 do 1 i=1,j 1 a(j,i,k)=a(i,j,k) close(90) write(6,*)fn//' read' return end SUBROUTINE TTT1(ALPHA,A,G,ICO,coo,iwr) c Tensors, not derivatives C Transforms A and G to local origins (ICO=1), c common origin (ICO=2), C sets them to zero (ICO=3). C coo supposed to be passed in aus C the static limit G/w is in G C implicit none INTEGER*4 ICO,IA,IB,IG,ID,IC,iwr real*8 ALPHA(3,3),G(3,3),A(3,3,3),coo(3),SIGN,SUM,EPS C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 C DO 2 IA=1,3 DO 2 IB=1,3 SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*0.5d0*EPS(IB,IG,ID)*coo(IG)*ALPHA(IA,ID) G(IA,IB)=G(IA,IB)+SUM IF(ICO.EQ.3)G(IA,IB)=0.0d0 2 CONTINUE if(iwr.eq.1)then IF(ICO.EQ.1)WRITE(6,*)' G transformed into overlap origin' IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' G set to zero' endif c A .. last index (IA) is electric DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 SUM=0.0d0 IF(IB.EQ.IC)SUM= 1coo(1)*ALPHA(IA,1)+coo(2)*ALPHA(IA,2)+coo(3)*ALPHA(IA,3) A(IC,IB,IA)=A(IC,IB,IA)+ 1SIGN*(SUM-1.5d0*(coo(IB)*ALPHA(IA,IC)+coo(IC)*ALPHA(IA,IB))) IF(ICO.EQ.3)A(IC,IB,IA)=0.0d0 3 CONTINUE if(iwr.eq.1)then IF(ICO.EQ.1)WRITE(6,*)' A transformed into overlap origin' IF(ICO.EQ.2)WRITE(6,*)' A transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' A set to zero' endif RETURN END c