program x00C implicit none integer*4 i,j,nat,n,ii,ii1,ix,jj,jj1,jx,nat1,ir, 1i1,i2,i3,ic,iargc real*8 u(3,3),x,y,z,cm(3),ac,bc,cc,ap,bp,cp,h, 1xi,yi,zi,dc,d,tol real*8,allocatable::r(:,:),P(:,:,:),A(:,:,:), 1f(:,:),f1(:,:),P1(:,:,:),A1(:,:,:), 1alpha(:,:,:),gp(:,:,:),aa(:,:,:,:), 1alpha1(:,:,:),gp1(:,:,:),aa1(:,:,:,:) integer*4,allocatable::ik(:),id(:) character*10 s10 h=0.5d0 tol=0.1d0 write(6,601) 601 format(' Input: FILE.X = part of a crystal',/, 1 ' FILE.FC',/, 1 ' FILE.TEN',/, 1 ' FILE.TTT',/, 1 ' CELL.TXT',/, 1 ' Output: FILE0.X, el. cell geom.',/, 1 ' FILE0.FC, el. cell FF',/, 1 ' FILE0.TTT, el. cell PD',/, 1 ' FILE0.TEN, el. cell DD',/) nat1=0 if(iargc().gt.0)then call getarg(1,s10) read(s10,*)nat1 endif open(9,file='FILE.X',status='old') read(9,*) read(9,*)nat n=3*nat allocate(r(nat,3),ik(nat),f(n,n),P(nat,3,3),A(nat,3,3), 1alpha(n,3,3),gp(n,3,3),aa(n,3,3,3),id(nat)) cm=0.0d0 do 1 i=1,nat read(9,*)ik(i),(r(i,j),j=1,3) do 1 j=1,3 1 cm(j)=cm(j)+r(i,j) cm=cm/dble(nat) close(9) write(6,605)nat 605 format(' FILE.X,',i5,' atoms') open(12,file='CELL.TXT',status='old') read(12,*)(u(1,i),i=1,3) read(12,*)(u(2,i),i=1,3) read(12,*)(u(3,i),i=1,3) close(12) ac=dsqrt(u(1,1)**2+u(1,2)**2+u(1,3)**2) bc=dsqrt(u(2,1)**2+u(2,2)**2+u(2,3)**2) cc=dsqrt(u(3,1)**2+u(3,2)**2+u(3,3)**2) write(6,604)ac,bc,cc 604 format(' CELL.TXT read, a b c:',3f12.4) if(nat1.eq.0)then nat1=0 do 5 i=1,nat x=r(i,1)-cm(1) y=r(i,2)-cm(2) z=r(i,3)-cm(3) c elementary cell projections cp= ((x*u(1,2)-y*u(1,1))*(u(2,2)*u(1,3)-u(2,3)*u(1,2)) 1 -(y*u(1,3)-z*u(1,2))*(u(2,1)*u(1,2)-u(2,2)*u(1,1)))/ 1 ((u(3,1)*u(1,2)-u(3,2)*u(1,1))*(u(2,2)*u(1,3)-u(2,3)*u(1,2)) 1 -(u(3,2)*u(1,3)-u(3,3)*u(1,2))*(u(2,1)*u(1,2)-u(2,2)*u(1,1))) bp= (x*u(1,2)-y*u(1,1)-cp*(u(3,1)*u(1,2)-u(3,2)-u(1,1)))/ 1 (u(2,1)*u(1,2)-u(2,2)*u(1,1)) ap=(x-bp*u(2,1)-cp*u(3,1))/u(1,1) ap=x/u(1,1) bp=y/u(2,2) cp=z/u(3,3) if(ap.gt.-h.and.ap.le.h.and. 1 bp.gt.-h.and.bp.le.h.and. 1 cp.gt.-h.and.cp.le.h)then nat1=nat1+1 id(nat1)=i endif 5 continue else do 51 i=1,nat1 51 id(i)=i endif write(6,*)nat1,' atoms in the central cell' write(6,766)(id(i),i=1,nat1) 766 format(20i4) CALL READFF(n,f) CALL READTEN(nat,P,A) CALL READTTT(ALPHA,AA,GP,NAT) allocate(f1(nat1*3,nat1*3),P1(nat1,3,3),A1(nat1,3,3), 1alpha1(3*nat1,3,3),gp1(3*nat1,3,3),aa1(3*nat1,3,3,3)) c rewrite tensors and geometry do 7 i=1,nat1 7 call apol(i,id(i),nat1,nat,P1,A1,alpha1,gp1,aa1,P,A,alpha,gp,aa) c dynamic matrix: f1=0.0d0 do 2 i=1,nat1 do 2 j=1,nat c atom j, what is the related one in the cell? ir=0 dc=1.0d9 ic=0 do 22 i1=-3,3 do 22 i2=-3,3 do 22 i3=-3,3 x=r(j,1)+dble(i1)*u(1,1)+dble(i2)*u(2,1)+dble(i3)*u(3,1) y=r(j,2)+dble(i1)*u(1,2)+dble(i2)*u(2,2)+dble(i3)*u(3,2) z=r(j,3)+dble(i1)*u(1,3)+dble(i2)*u(2,3)+dble(i3)*u(3,3) do 22 ii=1,nat1 xi=r(id(ii),1) yi=r(id(ii),2) zi=r(id(ii),3) d=dsqrt((x-xi)**2+(y-yi)**2+(z-zi)**2) if(ik(j).eq.ik(id(ii)))then if(d.lt.dc)then dc=d ic=ii endif if(d.lt.tol)then ir=ii goto 221 endif endif 22 continue 221 if(ir.eq.0)then write(6,664)j,ic,dc 664 format(' Atom ',i4,' not assigned, closest to',i4,f8.4,' A') stop endif do 24 ix=1,3 ii1=ix+3*(i-1) ii =ix+3*(id(i)-1) do 24 jx=1,3 jj1=jx+3*(ir-1) jj =jx+3*(j-1) 24 f1(ii1,jj1)=f1(ii1,jj1)+f(ii,jj) 2 continue open(9,file='FILE0.X') write(9,*)'0 elementary cell' write(9,*)nat1 do 4 i=1,nat1 4 write(9,900)ik(id(i)),(r(id(i),j),j=1,3) 900 format(i5,3f12.6,' 0 0 0 0 0 0 0 0.0') close(9) write(6,*)'FILE0.X written' call WRITEFF(3*nat1,f1) CALL WRITETEN(P1,A1,NAT1) call WRITETTT(NAT1,ALPHA1,AA1,GP1,'FILE0.TTT') end SUBROUTINE WRITETEN(P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(NAT,3,3),A(NAT,3,3) OPEN(15,FILE='FILE0.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 220 L=1,NAT DO 220 J=1,3 220 WRITE(15,1501) (A(L,J,I),I=1,3),L Z0=0.0d0 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) WRITE(6,*)' FILE0.TEN written ...' RETURN END SUBROUTINE READFF(N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) open(20,file='FILE.FC') N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,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+1,N 31 FCAR(I,J)=FCAR(J,I) close(20) write(6,*)'FILE.FC read' RETURN END SUBROUTINE READTEN(nat,P,A) IMPLICIT none integer*4 nat,L,I,J REAL*8 P(nat,3,3),A(nat,3,3) open(15,file='FILE.TEN',status='old') READ (15,*) C 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 DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (A(L,J,I),I=1,3) WRITE(6,*)' VCD and IR parameters read in ...' C close(15) RETURN END C SUBROUTINE report(s) character*(*) s write(6,*)s stop end SUBROUTINE WRITEFF(N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) open(20,file='FILE0.FC') 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) write(6,*)' FILE0.FC written' close(20) RETURN END SUBROUTINE READTTT(ALPHA,A,G,NAT) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3) logical lex inquire(file='FILE.TTT',exist=lex) if(lex)then open(15,file='FILE.TTT',status='old') READ(15,*) READ(15,*)NATC if(NAT.NE.NATC)call report('Wrong number of atoms in .TTT file') 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) write(6,*)'FILE.TTT read' else DO 4 I=1,3 DO 4 L=1,NAT DO 4 IX=1,3 IIND=3*(L-1)+IX DO 4 J=1,3 ALPHA(IIND,I,J)=0.0d0 4 G(IIND,I,J)=0.0d0 DO 5 I=1,3 DO 5 J=1,3 DO 5 L=1,NAT DO 5 IX=1,3 IIND=3*(L-1)+IX DO 5 K=1,3 5 A(IIND,I,J,K)=0.0d0 write(6,*)'FILE.TTT not found' endif RETURN END SUBROUTINE WRITETTT(NAT,ALPHA,A,G,s) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*nat,3,3),G(3*nat,3,3),A(3*nat,3,3,3) character*(*)s OPEN(2,FILE=s) 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,3F18.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) CLOSE(2) WRITE(6,*)s//' written ' RETURN END subroutine apol(i,iu,nat1,nat,P1,A1,ala1,gp1,aa1,P,A,ala,gp,aa) implicit none integer*4 i,iu,ix,ii1,ii,jx,kx,lx,nat1,nat real*8 P1(nat1,3,3),A1(nat1,3,3),ala1(3*nat1,3,3), 1gp1(3*nat1,3,3),aa1(3*nat1,3,3,3),P(nat,3,3),A(nat,3,3), 1ala(3*nat,3,3),gp(3*nat,3,3),aa(3*nat,3,3,3) do 31 ix=1,3 ii1=ix+3*(i -1) ii =ix+3*(iu-1) do 31 jx=1,3 P1(i,ix,jx)=P(iu,ix,jx) A1(i,ix,jx)=A(iu,ix,jx) do 31 kx=1,3 ala1(ii1,jx,kx)=ala(ii,jx,kx) gp1( ii1,jx,kx)=gp( ii,jx,kx) do 31 lx=1,3 31 aa1( ii1,jx,kx,lx)=aa( ii,jx,kx,lx) return end