PROGRAM MAKECRYSTALFC implicit none integer*4 ncon,nat,nxstep,nystep,nzstep, 1ix,iy,iz,iqx,iqy,iqz,xstart,ystart,zstart,ia,iac,ic,iind,iindc, 1jind,jindc,jacn,nacn,iix,jjx,j,ja,jx,kx,lx real*8 qxstep,qystep,qzstep,qxmin,qymin,qzmin,av(3),bv(3),cv(3), 1pi,a,b,c,ab,ac,bc,alpha,beta,gamma,ca,cb,cg,sg,at(3),bt(3),ct(3), 1sp,qx,qy,qz,x,y,z,spr,ssp,csp,qqx,qqy,qqz real*8,allocatable ::FCAR(:,:),DI(:,:),DR(:,:), 1alphabig(:,:,:),abig(:,:,:,:),gbig(:,:,:),PB(:,:,:),ABT(:,:,:), 2alphat(:,:,:),aat(:,:,:,:),g(:,:,:),PS(:,:,:),AS(:,:,:) character*20 sx,sy,sz logical lex c WRITE(6,3000) 3000 FORMAT(' Takes geometry and force field of a crystal cell',/, 1 ' and makes the dynamic matrices',/,/, 1 ' Input: CELL.X elementary cell geometry',/, 2 ' CELL.TXT elementary cell vectors',/, 2 ' CELL.FC elementary cell force field',/, 3 ' XYZ.TXT number of cells in x,y, and z', 3 ' (must be 3x3x3)',/, 3 ' Q.TXT definition of q-vectors', 3 ' (qxmin qxstep nxsteps',/, 3 ' qymin qystep nysteps',/, 3 ' qzmin qzstep nzsteps)',/, 3 ' Output: .FC FFs for individual q-vectors',/,/) open(2,file='CELL.TXT') read(2,*)(av(j),j=1,3) read(2,*)(bv(j),j=1,3) read(2,*)(cv(j),j=1,3) close(2) pi=4.0d0*atan(1.0d0) a=dsqrt(sp(av,av)) b=dsqrt(sp(bv,bv)) c=dsqrt(sp(cv,cv)) ab=sp(av,bv) ac=sp(av,cv) bc=sp(bv,cv) 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,/,/) ca=cos(pi*alpha/180.0d0) cb=cos(pi*beta /180.0d0) cg=cos(pi*gamma/180.0d0) sg=sin(pi*gamma/180.0d0) at(1)=a at(2)=0.0d0 at(3)=0.0d0 bt(1)=b*cg bt(2)=b*sg bt(3)=0.0d0 ct(1)=c*cb ct(2)=c*(ca-cg*cb)/sg ct(3)=dsqrt(c**2-ct(1)**2-ct(2)**2) write(6,601)at,bt,ct 601 format(' Translational vectors: ',/,3(3f12.4,/)) open(2,file='CELL.X') read(2,*) read(2,*)nat close(2) write(6,*)nat,' atoms in CELL.X' inquire(file='Q.TXT',exist=lex) if(lex)then open(2,file='Q.TXT') read(2,*)qxmin,qxstep,nxstep read(2,*)qymin,qystep,nystep read(2,*)qzmin,qzstep,nzstep close(2) else qxmin=0.0d0 qymin=0.0d0 qzmin=0.0d0 qxstep=0.0d0 qystep=0.0d0 qzstep=0.0d0 nxstep=1 nystep=1 nzstep=1 endif write(6,6005)qxmin,qxstep,nxstep,qymin,qystep,nystep, 1qzmin,qzstep,nzstep,nxstep*nystep*nzstep 6005 format(/,' q-vectors:',/, 1' min/a del q nsteps',/, 1'X: ',f10.3,f10.3,i7,/, 1'Y: ',f10.3,f10.3,i7,/, 1'Z: ',f10.3,f10.3,i7,/,/, 1i6,' in total',/) ic=0 allocate(FCAR(1:3*nat,1:3*nat),stat=ic) if(ic.ne.0)call report('fcar cannot be allocated') open(20,file='CELL.FC') call READFF(3*nat,3*nat,FCAR) close(20) write(6,*)'CELL.FC read in' allocate(DI(1:3*nat,1:3*nat),stat=ic) if(ic.ne.0)call report('DI cannot be allocated') allocate(DR(1:3*nat,1:3*nat),stat=ic) if(ic.ne.0)call report('DR cannot be allocated') qx=qxmin-qxstep do 1 iqx=1,nxstep qx=qx+qxstep qy=qymin-qystep do 1 iqy=1,nystep qy=qy+qystep qz=qzmin-qzstep do 1 iqz=1,nzstep qz=qz+qzstep write(6,6006)qx,qy,qz 6006 format('q:',3f10.3) qqx=(qx*at(1)/a**2+qy*bt(1)/b**2+qz*ct(1)/c**2)/2.0d0/pi qqy=(qx*at(2)/a**2+qy*bt(2)/b**2+qz*ct(2)/c**2)/2.0d0/pi qqz=(qx*at(3)/a**2+qy*bt(3)/b**2+qz*ct(3)/c**2)/2.0d0/pi do 3 ix=1,3*nat do 3 iy=1,3*nat DI(ix,iy)=0.0d0 3 DR(ix,iy)=0.0d0 nacn=-nat do 2 ix=1,3 do 2 iy=1,3 do 2 iz=1,3 nacn=nacn+nat x=dble(ix-2)*at(1)+dble(iy-2)*bt(1)+dble(iz-2)*ct(1) y=dble(ix-2)*at(2)+dble(iy-2)*bt(2)+dble(iz-2)*ct(2) z=dble(ix-2)*at(3)+dble(iy-2)*bt(3)+dble(iz-2)*ct(3) spr=x*qqx+y*qqy+z*qqz ssp=sin(spr) csp=cos(spr) do 2 ia=1,nat c atoms in central - 14th -cell, off-set by 13: iac=ia+13*nat do 2 iix=1,3 iind=iix+(ia-1)*3 iindc=iix+(iac-1)*3 do 2 ja=1,nat jacn=nacn+ja do 2 jjx=1,3 jind=jjx+(ja-1)*3 jindc=jjx+(jacn-1)*3 DI(iind,jind)=DI(iind,jind)+FCAR(iindc,jindc)*ssp 2 DR(iind,jind)=DR(iind,jind)+FCAR(iindc,jindc)*csp call sset(sx,iqx,xstart) call sset(sy,iqy,ystart) call sset(sz,iqz,zstart) if(lex)then open(20,file= 1 sx(xstart:20)//'.'//sy(ystart:20)//'.'//sz(zstart:20)//'FI.FC') call WRITEFF(3*nat,3*nat,DI) close(20) open(20,file= 1 sx(xstart:20)//'.'//sy(ystart:20)//'.'//sz(zstart:20)//'FR.FC') call WRITEFF(3*nat,3*nat,DR) close(20) else open(20,file='CELL.FC') call WRITEFF(3*nat,3*nat,DR) close(20) write(6,*)' CELL.FC extracted from CRYSTAL.FC' endif 1 continue stop end c ============================================================ SUBROUTINE WRITETEN(fn,NAT,P,A) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(NAT,3,3),A(NAT,3,3) character*(*) fn open(15,file=fn) 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) RETURN END c ============================================================ SUBROUTINE READTTT(fn,NAT,ALPHA,G,A) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ALPHA(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3) character*(*) fn OPEN(2,FILE=fn) 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 c ============================================================ SUBROUTINE WRITETTT(fn,NAT,ALPHA,A,G) 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*(*) fn OPEN(2,FILE=fn) 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 c ============================================================ SUBROUTINE READTEN(fn,nat,P,A) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(nat,3,3),A(nat,3,3) character*(*) fn open(15,file=fn) 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 c ============================================================ 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 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 subroutine sset(s,i,i1) implicit none character*20 s integer*4 i,i1 write(s,10)i 10 format(i20) do 1 i1=1,len(s) 1 if(s(i1:i1).ne.' ')return return end