PROGRAM MAKECRYSTAL IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) parameter (nt0=7000) character*80 name dimension r(nt0,3),ik(nt0),PS(nt0,3,3),AS(nt0,3,3),rt(3), 1at(3),bt(3),ct(3), 1ALPHAT(3*nt0,3,3),G(3*nt0,3,3),ATT(3*nt0,3,3,3), 1ALPHAbig(3*nt0,3,3),Gbig(3*nt0,3,3),Abig(3*nt0,3,3,3), 1PB(nt0,3,3),ABT(nt0,3,3),j33(nt0*3),a33(nt0*3) allocatable ::fcar(:,:),fbig(:,:) logical lex c WRITE(6,3000) 3000 FORMAT(' Takes geometry of a crystal cell, and propagates it',/, 1 ' Input: CELL.X(.FC,,.TEN,.TTT)',/, 2 ' CELL.TXT elementary cell vectors',/, 3 ' XYZ.TXT number of cell in x,y, and z',/, 3 ' Output: CRYSTAL.X(.FC,.TTT,.TEN,ACRYSTAL.FC)',/,/) open(2,file='CELL.TXT') read(2,*)(at(j),j=1,3) read(2,*)(bt(j),j=1,3) read(2,*)(ct(j),j=1,3) close(2) write(6,601)at,bt,ct 601 format(' Translational vectors: ',/,3(3f12.4,/)) pi=4.0d0*atan(1.0d0) a=dsqrt(sp(at,at)) b=dsqrt(sp(bt,bt)) c=dsqrt(sp(ct,ct)) ab=sp(at,bt) ac=sp(at,ct) bc=sp(bt,ct) alpha=180.0d0/pi*acos(bc/(b*c)) beta =180.0d0/pi*acos(ac/(a*c)) gamma=180.0d0/pi*acos(ab/(b*a)) write(6,600)alpha,a,beta,b,gamma,c 600 format(/,' Crystall cell parameters:',/,/, 1 ' alpha: ',f10.2,' a: ',f10.4,/, 1 ' beta: ',f10.2,' b: ',f10.4,/, 1 ' gamma: ',f10.2,' c: ',f10.4,/,/) open(2,file='XYZ.TXT') read(2,*)nx,ny,nz write(6,500)nx,ny,nz 500 format(/,10x,i2,'a X ',i2,'b X ',i2,'c',/) close(2) open(2,file='CELL.X') read(2,2000)name 2000 format(a80) read(2,*)nat if(nat.gt.nt0)then close(2) write(6,*)'too many atoms' stop endif do 3 i=1,nat 3 read(2,*)ik(i),(r(i,j),j=1,3) close(2) open(2,file='CRYSTAL.X') write(2,2000)name write(2,*)nat*nx*ny*nz do 2 ix=0,nx-1 do 2 iy=0,ny-1 do 2 iz=0,nz-1 do 2 i=1,nat rt(1)=r(i,1)+ix*at(1)+iy*bt(1)+iz*ct(1) rt(2)=r(i,2)+ix*at(2)+iy*bt(2)+iz*ct(2) rt(3)=r(i,3)+ix*at(3)+iy*bt(3)+iz*ct(3) 2 write(2,20)ik(i),(rt(j),j=1,3) 20 format(i2,3f15.8,' 0 0 0 0 0 0 0 0.0') close(2) write(6,*)' CRYSTAL.X was written' open(2,file='CCT.INP') write(2,2001)nx*ny*nz 2001 format('LDIA',/,'f',/,'LOFF',/,'f',/,'LNAMES',/, q 'f',/,'LRAM',/,'t',/,'OLDFORMAT',/,'f',/, 1 'LROA',/,'t',/,'LABS',/,'t',/,'LVCD',/,'t',/,'POLYMER',/, 1 i4) io=0 do 21 ix=0,nx-1 do 21 iy=0,ny-1 do 21 iz=0,nz-1 io=io+1 write(2,*)io,' overlap:' write(2,2002)(I,I=1,nat*nx*ny*nz) 2002 format(20I6) do 21 i=1,nat*nx*ny*nz if(i.gt.(io-1)*nat.and.i.le.io*nat)then iw=i-(io-1)*nat else iw=0 endif write(2,2003)iw 2003 format(i6,' ',$) if(mod(i,20).eq.0.or.i.eq.nat*nx*ny*nz)write(2,*) 21 continue close(2) write(6,*)'CCT.INP written' inquire(file='CELL.FC',exist=lex) if (lex)then if(nat*nx*ny*nz.gt.nt0)then write(6,*)'too many atoms' stop endif c ic=0 allocate(FCAR(1:3*nat,1:3*nat),stat=ic) if(ic.ne.0)call report('fcar cannot be allocated') allocate(FBIG(1:3*nat*nx*ny*nz,1:3*nat*nx*ny*nz),stat=ic) if(ic.ne.0)call report('fbig cannot be allocated') do 41 ii=1,3*nat*nx*ny*nz do 41 jj=1,3*nat*nx*ny*nz 41 fbig(jj,ii)=0.0d0 open(20,file='CELL.FC') call READFF(3*nat,3*nat,FCAR) close(20) do 4 ic=1,nx*ny*nz do 4 i=1,nat ia=(ic-1)*nat+i do 4 j=1,nat ja=(ic-1)*nat+j do 4 ix=1,3 ii=ix+3*(ia-1) do 4 jx=1,3 jj=jx+3*(ja-1) 4 fbig(ii,jj)=fcar(ix+3*(i-1),jx+3*(j-1)) open(20,file='CRYSTAL.FC') call WRITEFF(3*nat*nx*ny*nz,3*nat*nx*ny*nz,fbig) close(20) write(6,*)' CRYSTAL.FC written' open(33,file='ACRYSTAL.FC',form='unformatted') flim=2.0d-7 do 8 i=1,3*nat*nx*ny*nz nn=0 do 81 j=i,3*nat*nx*ny*nz ft=fbig(i,j) if(abs(ft).gt.flim)then nn=nn+1 j33(nn)=j if(i.eq.j)ft=ft/2.0d0 a33(nn)=ft endif 81 continue 8 write(33)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) close(33) write(6,*)' ACRYSTAL.FC written' deallocate(FCAR,fbig) endif inquire(file='CELL.TTT',exist=lex) if (lex)then call READTTT(3*nt0,NAT,ALPHAT,G,ATT) ia=0 do 6 i1=0,nx-1 do 6 i2=0,ny-1 do 6 i3=0,nz-1 do 6 i=1,nat ia=ia+1 do 6 ix=1,3 ii=ix+3*(ia-1) do 6 jx=1,3 do 6 kx=1,3 alphabig(ii,jx,kx)=ALPHAT(ix+3*(i-1),jx,kx) Gbig(ii,jx,kx)=G(ix+3*(i-1),jx,kx) do 6 lx=1,3 6 Abig(ii,jx,kx,lx)=ATT(ix+3*(i-1),jx,kx,lx) call WRITETTT(nt0,NAT*nx*ny*nz,ALPHABIG,Abig,Gbig) write(6,*)' CRYSTAL.TTT written' endif inquire(file='CELL.TEN',exist=lex) if (lex)then CALL READTEN(nt0,PS,AS,NAT) ia=0 do 7 i1=0,nx-1 do 7 i2=0,ny-1 do 7 i3=0,nz-1 do 7 i=1,nat ia=ia+1 do 7 ix=1,3 do 7 jx=1,3 PB(ia,ix,jx)=PS(i,ix,jx) 7 ABT(ia,ix,jx)=AS(i,ix,jx) call WRITETEN(nt0,PB,ABT,NAT*nx*ny*nz) endif stop 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) open(15,file='CRYSTAL.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,*)' CRYSTAL.TEN written ...' RETURN END SUBROUTINE READFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) 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) WRITE(*,*)' Cartesian FF read in ... ' RETURN END SUBROUTINE WRITEFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) 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 READTTT(N30,NAT,ALPHA,G,A) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3) OPEN(2,FILE='CELL.TTT',STATUS='OLD') 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 1 READ(2,*)LDUM,IXDUM,(ALPHA(IND,I,J),J=1,3) 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 index magnetic, inner 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) CLOSE(2) RETURN END 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) OPEN(2,FILE='CRYSTAL.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(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) RETURN END SUBROUTINE READTEN(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) open(15,file='CELL.TEN') READ (15,*) NOAT NAT=NOAT 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,NOAT DO 220 J=1,3 220 READ (15,*) (A(L,J,I),I=1,3) close(15) RETURN END subroutine report(s) character*(*) s write(6,*)s stop end function sp(a,b) real*8 sp,a(*),b(*) sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end