program local c c ROA tensor derivatives/transform to local/common ors c implicit real*8(a-h,o-z) implicit integer*4(i-n) parameter (nat0=80,N30=3*nat0) dimension alder(N30,3,3),gder(N30,3,3),aader(N30,3,3,3), 1xv(3,nat0),adv(N30,3,3) dimension al(N30,3,3),g(N30,3,3),a(N30,3,3,3), 1alv(N30,3,3) character*80 fttt,fx c write(6,*)'Tensor file (with .TTT):' read(5,'(a)')fttt OPEN(2,FILE=fttt,STATUS='OLD') call READTEN(N30,NAT,alder,gder,aader,adv) CLOSE(2) WRITE(6,*)NAT,' atoms, tensors read in' write(6,*)'Coordinate file (with .X):' read(5,'(a)')fx OPEN(2,FILE=fx) call READX(xv,NAT,N30) CLOSE(2) WRITE(6,*)NAT,' atoms, geometry read in' write(6,*) write(6,609) 609 format('Choose the kind of transformation',/, 1' (1-into local',/, 1' 2-into common',/, 1' 3-zero out A, G',/, 1' 4-improve with ALPHA.TTT):') read(5,*)ico if(ico.eq.4)then OPEN(2,FILE='ALPHA.TTT',STATUS='OLD') call READTEN(N30,NAT,al,g,a,alv) CLOSE(2) call ttai(N30,al,alder,Aader,Gder,nat,5,xv) open(2,file='RES.TTT') call writettt(3*nat0,NAT,al,aader,gder,alder) close(2) else call ttai(N30,alder,adv,aader,gder,NAT,ico,xv) open(2,file='RES.TTT') call writettt(3*nat0,NAT,alder,aader,gder,adv) close(2) endif write(6,*)'Tensors written to RES.TTT' stop 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 READX(R,NAT,N30) IMPLICIT REAL*8(A-H,O-Z) implicit integer*4 (i-n) DIMENSION R(3,N30/3) BOHR=0.529177D0 READ(2,*) READ(2,*)NAT DO 1 I=1,NAT 1 READ(2,*)IDUM,(R(J,I),J=1,3) DO 2 I=1,NAT DO 2 J=1,3 2 R(J,I)=R(J,I)/BOHR RETURN END C ===================== SUBROUTINE READTEN(N30,NAT,ALPHA,G,A,alphav) IMPLICIT REAL*8 (A-H,O-Z) implicit integer*4 (i-n) DIMENSION ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3),alphav(N30,3,3) READ(2,*) READ(2,*)NAT READ(2,*) READ(2,*) DO 1 I=1,3 READ(2,*) DO 1 L=1,NAT DO 1 IX=1,3 IND=3*(L-1)+IX READ(2,*)LDUM,IXDUM,(ALPHA(IND,I,J),J=1,3) DO 1 J=1,3 1 alphav(IND,I,J)=ALPHA(IND,I,J) READ(2,*) READ(2,*) DO 2 I=1,3 READ(2,*) DO 2 L=1,NAT DO 2 IX=1,3 IND=3*(L-1)+IX c last fastest index magnetic (inner loop) 2 READ(2,*)LDUM,IXDUM,(G(IND,I,J),J=1,3) READ(2,*) READ(2,*) DO 3 I=1,3 DO 3 J=1,3 READ(2,*) DO 3 L=1,NAT DO 3 IX=1,3 IND=3*(L-1)+IX 3 READ(2,*)LDUM,IXDUM,(A(IND,I,J,K),K=1,3) READ(2,*) READ(2,*) DO 4 I=1,3 READ(2,*) DO 4 L=1,NAT DO 4 IX=1,3 IND=3*(L-1)+IX 4 READ(2,*)LDUM,IXDUM,(alphav(IND,I,J),J=1,3) RETURN END C ========================================== SUBROUTINE WRITETTT(N30,NAT,ALPHA,A,G,alphav) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3), 1alphav(N30,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 (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 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) WRITE(2,2009) 2009 FORMAT('The velocity form of alpha',/, 1' Atom/x jx jy jz') DO 11 I=1,3 WRITE(2,2012)I 2012 FORMAT(' AlphaV(',I1,',J):') DO 11 L=1,NAT DO 11 IX=1,3 IIND=3*(L-1)+IX 11 WRITE(2,2001)L,IX,(ALPHAV(IIND,I,J),J=1,3) RETURN END C ========================================== SUBROUTINE ttai(N,ALPHA,AV,A,G,NAT,ICO,R) C C ICO=1 ... transform A and G into local origins C ICO=2 ... transform A and G into common origin C ICO=3 ... zero out A and G C ICO=4 ... distributed origin transformation for G C ICO=5 ... improvement with better alpha, in AV is the worse one C C R supposed to be passed in atomic units, G=G/w C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,NAT),AV(N,3,3) C IF(ICO.EQ.1)SG = 1.0D0 IF(ICO.EQ.1)SG2 = 0.0D0 IF(ICO.EQ.2)SG = -1.0D0 IF(ICO.EQ.2)SG2 = 0.0D0 IF(ICO.EQ.4)SG = 1.0D0 IF(ICO.EQ.4)SG2 = 1.0D0 IF(ICO.EQ.5)SG = 1.0D0 IF(ICO.EQ.5)SG2 = 1.0D0 C DO 1 IA=1,3 DO 1 IB=1,3 DO 1 L=1,NAT DO 1 IE=1,3 K=3*(L-1)+IE if(ICO.EQ.3)then G(K,IA,IB)=0.0d0 else DO 2 IG=1,3 DO 2 ID=1,3 2 G(K,IA,IB)=G(K,IA,IB) 1 +SG *0.5D0*EPS(IB,IG,ID)*R(IG,L)*AV( K,IA,ID) 2 -SG2*0.5D0*EPS(IB,IG,ID)*R(IG,L)*ALPHA(K,IA,ID) c last index of G magnetic endif 1 continue c IF(ICO.EQ.1)WRITE(6,*)' G transformed into atomic origins with AV' IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin with AV' IF(ICO.EQ.3)WRITE(6,*)' G set to zero' IF(ICO.EQ.4)then WRITE(6,*)' DOG for G' return ENDIF IF(ICO.EQ.5)WRITE(6,*)' G inproved with ALPHA.TTT' 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 K=3*(L-1)+IE IF(ICO.EQ.3)then A(K,IA,IB,IC)=0.0D0 else S1=0.0D0 S2=0.0D0 IF(IB.EQ.IC)THEN DO 4 ID=1,3 S2=S2+R(ID,L)* AV(K,IA,ID) 4 S1=S1+R(ID,L)*ALPHA(K,IA,ID) ENDIF if(ICO.EQ.5)then A(K,IA,IB,IC)=A(K,IA,IB,IC) 1 -(S1-1.5D0*(R(IB,L)*ALPHA(K,IA,IC)+R(IC,L)*ALPHA(K,IA,IB))) 1 +(S2-1.5D0*(R(IB,L)*AV( K,IA,IC)+R(IC,L)*AV( K,IA,IB))) else A(K,IA,IB,IC)=A(K,IA,IB,IC) 1 +SG *(S1-1.5D0*(R(IB,L)*ALPHA(K,IA,IC)+R(IC,L)*ALPHA(K,IA,IB))) endif endif 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' IF(ICO.EQ.5)WRITE(6,*)' A inproved with ALPHA.TTT' c RETURN END