program x00 implicit none integer*4 i,j,nat,iccub,n,ice,ii,ii1,iof,ix,jj,jj1,jx,nat1,jof,iu, 1n27,i1,i2,i3,iargc,ifind,ij,k real*8 u(3,3),x,y,z,s1,s2,s3 real*8,allocatable::r(:,:),P(:,:,:),A(:,:,:), 1f(:,:),f1(:,:),P1(:,:,:),A1(:,:,:),r1(:,:),r0(:,:), 1alpha(:,:,:),gp(:,:,:),aa(:,:,:,:), 1alpha1(:,:,:),gp1(:,:,:),aa1(:,:,:,:) integer*4,allocatable::ik(:),ik1(:),ik0(:) character*10 s10 logical lfix write(6,601) 601 format(' Input: FILE.X = CRYSTAL.X',/, 1 ' FILE.FC',/, 1 ' FILE.TEN',/, 1 ' FILE.TTT',/, 1 ' (CELL.TXT,CELL.X)',/, 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',/) 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)) do 1 i=1,nat 1 read(9,*)ik(i),(r(i,j),j=1,3) close(9) n27=27 u=0.0d0 lfix=iargc().eq.3 if(iargc().ge.1)then if(lfix)then call getarg(1,s10) read(s10,*)s1 call getarg(2,s10) read(s10,*)s2 call getarg(3,s10) read(s10,*)s3 write(6,600)s1,s2,s3 600 format(' Shifts:',3f10.2) 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) write(6,*)'CELL.TXT read' open(19,file='CELL.X',status='old') read(19,*) read(19,*)nat1 allocate(ik0(nat1),r0(nat1,3)) do 111 i=1,nat1 111 read(19,*)ik0(i),(r0(i,j),j=1,3) close(19) write(6,*)'CELL.X read,',nat1,' atoms' else if(mod(nat,nat1).ne.0)call report('mod nat nat1 <> 0') n27=nat/nat1 write(6,*)n27,' cells' write(6,*)nat1,' atoms in cell' endif else nat1=nat/n27 if(mod(nat,n27).ne.0)call report('mod nat n27 <> 0') write(6,*)n27,' cells' write(6,*)nat1,' atoms in cell' endif CALL READFF(n,f) CALL READTEN(nat,P,A) CALL READTTT(ALPHA,AA,GP,NAT) allocate(ik1(nat1),r1(nat1,3)) 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 take geometry and tensors from the central cube: f1=0.0d0 if(.not.lfix)then c determine central cube: ice=iccub(nat,r,n27) iof=nat1*(ice-1) do 2 i=1,nat1 ik1(i)=ik(i+iof) call apol(i,i+iof,nat1,nat,P1,A1,alpha1,gp1,aa1,P,A,alpha,gp,aa) do 2 ix=1,3 r1(i,ix)=r(i+iof,ix) ii1=ix+3*(i -1) ii =ix+3*(i+iof-1) c dynamic matrix to other cubes: do 2 j=1,nat1 do 2 jx=1,3 jj1=jx+3*(j -1) do 2 iu=1,n27 jof=nat1*(iu-1) jj =jx+3*(j+jof-1) 2 f1(ii1,jj1)=f1(ii1,jj1)+f(ii,jj) else do 29 i=1,nat1 ik1(i)=ik0(i) do 27 k=1,3 27 r1(i,k)=r0(i,k)+s1*u(1,k)+s2*u(2,k)+s3*u(3,k) ii=ifind(i,r1(i,1),r1(i,2),r1(i,3),r,nat,ik1(i),ik) if(ii.eq.0)then write(6,603)i,ik1(i),r1(i,1),r1(i,2),r1(i,3) 603 format(' Atom',i6,i3,3f10.4,' not found') stop endif 29 call apol(i,ii,nat1,nat,P1,A1,alpha1,gp1,aa1,P,A,alpha,gp,aa) write(6,*)'Diagonal terms done' c dynamic matrix to other cubes: do 3 j=1,nat1 do 3 i1=-1,1 do 3 i2=-1,1 do 3 i3=-1,1 c coordinates of atom j in cube iu=(i1,i2,i3): x=r1(j,1)+u(1,1)*dble(i1)+u(2,1)*dble(i2)+u(3,1)*dble(i3) y=r1(j,2)+u(1,2)*dble(i1)+u(2,2)*dble(i2)+u(3,2)*dble(i3) z=r1(j,3)+u(1,3)*dble(i1)+u(2,3)*dble(i2)+u(3,3)*dble(i3) ij=ifind(j,x,y,z,r,nat,ik1(j),ik) if(ij.eq.0)then write(6,602)j,i1,i2,i3 602 format(i6,3i3,' skipped') else do 31 i=1,nat1 ii=ifind(i,r1(i,1),r1(i,2),r1(i,3),r,nat,ik1(i),ik) do 31 ix=1,3 ii1=ix+3*(i-1) do 31 jx=1,3 jj1=jx+3*(j-1) 31 f1(ii1,jj1)=f1(ii1,jj1)+f(ix+3*(ii-1),jx+3*(ij-1)) endif 3 continue endif open(9,file='FILE0.X') write(9,*)'0 elementary cell' write(9,*)nat1 do 4 i=1,nat1 4 write(9,900)ik1(i),(r1(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 function iccub(nat,r,n27) c determine central cube implicit none integer*4 nat,nat1,ix,iu,ii,iccub,ice,i,n27 real*8 cet(3),r(nat,3),xi,yi,zi,d real*8,allocatable::cei(:,:) allocate(cei(n27,3)) nat1=nat/n27 c center of all and of each one: do 40 ix=1,3 cet(ix)=0.0d0 do 42 iu=1,n27 cei(iu,ix)=0.0d0 ii=nat1*(iu-1) do 41 i=ii+1,ii+nat1 cet( ix)=cet( ix)+r(i,ix) 41 cei(iu,ix)=cei(iu,ix)+r(i,ix) 42 cei(iu,ix)=cei(iu,ix)/dble(nat1) 40 cet(ix)=cet(ix)/dble(nat) write(6,600)cet 600 format(/,' CRYSTAL.X center (A):',7x,3f7.2,/,' Cube centers (A):') write(6,601)(iu,(cei(iu,ix),ix=1,3),iu=1,n27) 601 format(3(i3,':',3f7.2)) d=1.0d10 ice=0 do 43 iu=1,n27 xi=cei(iu,1)-cet(1) yi=cei(iu,2)-cet(2) zi=cei(iu,3)-cet(3) if(dsqrt(xi**2+yi**2+zi**2).lt.d)then d=dsqrt(xi**2+yi**2+zi**2) ice=iu endif 43 continue if(ice.eq.0)then call report('Central cube not found') else write(6,602)ice,d 602 format(' Central cube:',i3,f10.2,' A from the center') iccub=ice endif return end C 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 function ifind(i,x,y,z,r,nat,ik1i,ik) implicit none integer*4 i,nat,ia,ik1i,ik(nat),ifind real*8 x,y,z,r(nat,3),d do 25 ia=1,nat d=(r(ia,1)-x)**2+(r(ia,2)-y)**2+(r(ia,3)-z)**2 if(d.lt.1.0d-5.and.ik(ia).eq.ik1i)then ifind=ia return endif 25 continue ifind=0 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