program unifyc implicit none integer*4 nat,n,i,ix,iy,n27,nats,ia,ja,ic,ice,icep,icp, 1ie,iep,iofc,iofep,iofe,jc,jcp,nx,ny,nz,nxp,nyp,nzp,iofcp,iz, 1iie,iic real*8 a(3),b(3),c(3) real*8 ,allocatable::r(:),f(:,:),P(:,:,:),X(:,:,:),ALPHA(:,:,:), 1G(:,:,:),AA(:,:,:,:) integer*4 ,allocatable::ip(:,:) logical lff,lten,lttt write(6,600) 600 format(' In a crystal piece distributes force field,',/, 1 ' and tensors from the central cell to the others,',/, 1 ' from the central cell to the others,',/, 1 ' also central-other cell interactions',/,/, 1 ' Input: CELL.TXT cell vectors',/, 1 ' CELL.X cell geometry',/, 1 ' FILE.X geometry',/, 1 ' FILE.FC force field',/, 1 ' FILE.TEN dipole dervs',/, 1 ' FILE.TTT polarizability dervs',/,/, 1 ' Output: FILEU.FC,etc unified files',/,/) open(9,file='CELL.TXT') read(9,*)a read(9,*)b read(9,*)c close(9) open(9,file='CELL.X') read(9,*) read(9,*)nats close(9) open(9,file='FILE.X') read(9,*) read(9,*)nat n=3*nat allocate(r(n),f(n,n),P(nat,3,3),X(nat,3,3),ALPHA(n,3,3), 1G(n,3,3),AA(n,3,3,3)) do 1 i=1,nat 1 read(9,*)(r(ix+3*(i-1)),ix=1,1),(r(ix+3*(i-1)),ix=1,3) close(9) f=0.0d0 inquire(file='FILE.FC',exist=lff) if(lff)CALL READFF(n,f) c load tensors and put them to atomic origins: inquire(file='FILE.TEN',exist=lten) if(lten)CALL READTEN(nat,P,X,.false.,r) inquire(file='FILE.TTT',exist=lttt) if(lttt)then CALL READTTT(NAT,ALPHA,AA,G) CALL TTTT(n,ALPHA,AA,G,1,r) endif n27=nat/nats allocate(ip(n27,3)) ip(1,1)=0 c determine central cube: call iccub(nat,r,n27,ip,a,b,c,ice) iofe=nats*(ice-1) c transfer central FF to all others write(6,608) 608 format(' cells: on:)') do 2 ic=1,n27 iofc=nats*(ic-1) nx=ip(ic,1)-ip(ice,1) ny=ip(ic,2)-ip(ice,2) nz=ip(ic,3)-ip(ice,3) if(ic.ne.ice)then c take tensors in ic from the central: do 9 ia=1,nats ie =ia+iofe jc =ia+iofc if(lten)then do 91 ix=1,3 do 91 iy=1,3 X(jc,ix,iy)=X(ie,ix,iy) 91 P(jc,ix,iy)=P(ie,ix,iy) endif if(lttt)then do 92 ix=1,3 iie=ix+3*(ie-1) iic=ix+3*(jc-1) do 92 iy=1,3 do 92 iz=1,3 ALPHA(iic,iy,iz )=ALPHA(iie,iy,iz) AA( iic,iy,iz,1)=AA( iie,iy,iz,1) AA( iic,iy,iz,2)=AA( iie,iy,iz,2) AA( iic,iy,iz,3)=AA( iie,iy,iz,3) 92 G( iic,iy,iz )=G( iie,iy,iz) endif 9 continue endif if(lff)then c ice...ic to all other pairs: do 21 icep=1,n27 iofep=nats*(icep-1) if(icep.ne.ice)then do 4 icp=1,n27 iofcp=nats*(icp-1) nxp=ip(icp,1)-ip(icep,1) nyp=ip(icp,2)-ip(icep,2) nzp=ip(icp,3)-ip(icep,3) if(nx.eq.nxp.and.ny.eq.nyp.and.nz.eq.nzp)then write(6,607)ice,ic,icep,icp 607 format(2i4,' --> ',2i4) do 8 ia=1,nats ie =ia+iofe iep=ia+iofep do 8 ja=1,nats jc =ja+iofc jcp=ja+iofcp do 8 ix=1,3 do 8 iy=1,3 8 f(ix+3*(iep-1),iy+3*(jcp-1))=f(ix+3*(ie-1),iy+3*(jc-1)) endif 4 continue endif 21 continue endif 2 continue if(lff)CALL WRITEFF(n,f) if(lten)CALL WRITETEN(NAT,P,X,R) if(lttt)then CALL TTTT(n,ALPHA,AA,G,2,r) CALL WRITETTT(n,ALPHA,AA,G,'FILEU.TTT') endif 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',STATUS='OLD') 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) DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) WRITE(*,*)' Cartesian FF read in ... ' CLOSE(20) RETURN END SUBROUTINE WRITEFF(N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) OPEN(20,FILE='FILEU.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) CLOSE(20) write(6,*)'FILEU.FC written.' RETURN END subroutine iccub(nat,r,n27,ip,wi,wj,wk,ice) c determine central cube implicit none integer*4 nat,nat1,ix,iu,ii,ice,i,n27,ip(n27,3) real*8 cet(3),r(*),xi,yi,zi,d,wi(3),wj(3),wk(3),down, 1vi,vj,vk,vii,vij,vik,vjj,vjk,vkk,sp,vu(3),a,b,c 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(3*(i-1)+ix) 41 cei(iu,ix)=cei(iu,ix)+r(3*(i-1)+ix) 42 cei(iu,ix)=cei(iu,ix)/dble(nat1) 40 cet(ix)=cet(ix)/dble(nat) write(6,600)nat,cet 600 format(/,' CRYSTAL.X:',i7,' atoms, center:', 13f7.2,' A',/,' 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,nat1*(ice-1)+1,nat1*ice,d 602 format(' Central cube:',i3,', atoms',i5,' -',i5,',', 1 f8.2,' A from the center') endif c relative positions to the central cell do 1 iu=1,n27 do 2 ix=1,3 2 vu(ix)=cei(iu,ix)-cei(ice,ix) vi=sp(vu,wi) vj=sp(vu,wj) vk=sp(vu,wk) vii=sp(wi,wi) vij=sp(wi,wj) vik=sp(wi,wk) vjj=sp(wj,wj) vjk=sp(wj,wk) vkk=sp(wk,wk) down=(vkk*vii-vik**2)*(vjj*vii-vij**2) 1-(vjk*vii-vij*vik)*(vjk*vii-vik*vij) c write(6,*)vi,vj,vk,vii,vij,vik,vjj,vjk,vkk,down c=(vi*(vjk*vij-vjj*vik)+vj*(vij*vik-vjk*vii)+vk*(vjj*vii-vij**2))* 1vii/down b=(vj*vii-vi*vij-c*(vjk*vii-vik*vij))/(vjj*vii-vij**2) a=(vi-b*vij-c*vik)/vii write(6,603)iu,a,b,c 603 format(i3,3f5.1,$) if(mod(iu,4).eq.0)write(6,*) ip(iu,1)=nint(a) ip(iu,2)=nint(b) 1 ip(iu,3)=nint(c) if(mod(n27,4).ne.0)write(6,*) return end SUBROUTINE report(s) character*(*) s write(6,*)s stop end function sp(v1,v2) real*8 v1(*),v2(*),sp sp=v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3) return end SUBROUTINE READTEN(NAT,P,A,LAPT,X) IMPLICIT none INTEGER*4 NAT,L,I,J,ID,IG REAL*8 P(NAT,3,3),A(NAT,3,3),X(*),BOHR,EPS,SU real*8,allocatable::AI(:,:,:),AJ(:,:,:),PV(:,:,:),R(:,:) LOGICAL LAPT OPEN(15,FILE='FILE.TEN') allocate(AI(NAT,3,3),AJ(NAT,3,3),PV(NAT,3,3),R(3,NAT)) BOHR=0.52917705993d0 DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(I+3*(L-1))/BOHR READ (15,*) DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) IF(LAPT)THEN DO 2203 L=1,NAT DO 2203 J=1,3 DO 2203 I=1,3 SU=0.0d0 DO 2204 ID=1,3 DO 2204 IG=1,3 2204 SU=SU+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) AJ(L,J,I)=0.0d0 PV(L,J,I)=P(L,J,I) 2203 AI(L,J,I)=SU WRITE(6,*)' APT part of AAT constructed' GOTO 1 ENDIF DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (AI(L,J,I),I=1,3) DO 230 L=1,NAT DO 230 J=1,3 230 READ (15,*) (AJ(L,J,I),I=1,3) DO 100 L=1,NAT DO 100 J=1,3 100 READ (15,*) (PV(L,J,I),I=1,3) C Distributed-origin gauge: 1 CALL VCDD0(AI,PV,P,R,NAT) WRITE(6,*)'DOG used' C C Total axial tensor = electronic + nuclear: DO 4 L=1,NAT DO 4 I=1,3 DO 4 J=1,3 4 A(L,I,J)=AI(L,I,J)+AJ(L,I,J) C C Put axial tensor in the origin on the atom L: DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SU=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SU=SU-0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SU WRITE(6,*)' AAT shifted to atomic origin' CLOSE(15) RETURN END 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 SUBROUTINE VCDD0(VCD,VEL,DIP,C,NAT) C STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) INTEGER*4 ALPHA,BETA,GAMMA,DELTA,LAMBDA DIMENSION VCD(NAT,3,3),VEL(NAT,3,3),DIP(NAT,3,3),C(3,NAT) DO 90 BETA =1,3 DO 90 GAMMA=1,3 DO 90 DELTA=1,3 SKEW=0.25D0*EPS(BETA,GAMMA,DELTA) IF (DABS(SKEW).GT.1.0D-10)THEN DO 80 ALPHA=1,3 DO 80 LAMBDA=1,NAT 80 VCD(LAMBDA,ALPHA,BETA)=VCD(LAMBDA,ALPHA,BETA) 1 -SKEW*C(GAMMA,LAMBDA) 2 *(VEL(LAMBDA,ALPHA,DELTA)-DIP(LAMBDA,ALPHA,DELTA)) ENDIF 90 CONTINUE RETURN END C SUBROUTINE READTTT(NAT,ALPHA,A,G) IMPLICIT none INTEGER*4 NAT,I,L,IX,IIND,J,K REAL*8 ALPHA(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3) OPEN(15,FILE='FILE.TTT') READ(15,*) READ(15,*) 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) RETURN END c ============================================================ SUBROUTINE TTTT(N,ALPHA,A,G,ICO,R) c Tensors derivatives C Transforms A and G to local origins (ICO=1) or c common origin (ICO=2); C sets A and G to zero for ICO=3. C R supposed to be passed in in A C IMPLICIT none INTEGER*4 N,NAT,I,J,L,ICO,IA,IB,IC,ID,IE,IG,IIND REAL*8 ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(*),W,EPS,SU,si real*8,allocatable::X(:,:) NAT=N/3 allocate(X(3,NAT)) DO 5 I=1,3 DO 5 J=1,NAT 5 X(I,J)=R(I+3*(J-1))/0.529177D0 C IF(ICO.EQ.1)si=1.0d0 IF(ICO.EQ.2)si=-1.0d0 IF(ICO.EQ.3)si=0.0d0 W=1.0d0 C supposedly the static limit G/w is calculated C DO 2 IA=1,3 DO 2 IB=1,3 DO 2 L=1,NAT DO 2 IE=1,3 IIND=3*(L-1)+IE SU=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SU=SU+si*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(IIND,IA,ID) G(IIND,IA,IB)=G(IIND,IA,IB)+SU IF(ICO.EQ.3)G(IIND,IA,IB)=0.0d0 2 CONTINUE IF(ICO.EQ.1)WRITE(6,*)' G transformed into atomic origins' IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' G set to zero' C c A IA index is electric 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 IIND=3*(L-1)+IE SU=0.0d0 IF(IB.EQ.IC)THEN DO 4 ID=1,3 4 SU=SU+X(ID,L)*ALPHA(IIND,IA,ID)*si ENDIF A(IIND,IA,IB,IC)=A(IIND,IA,IB,IC)+SU 1-1.5d0*si*(X(IB,L)*ALPHA(IIND,IA,IC)+X(IC,L)*ALPHA(IIND,IA,IB)) IF(ICO.EQ.3)A(IIND,IA,IB,IC)=0.0d0 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' RETURN END SUBROUTINE WRITETEN(NAT,P,A,X) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(NAT,3,3),A(NAT,3,3),X(*) real*8,allocatable::R(:,:) allocate(R(3,NAT)) BOHR=0.52917705993d0 DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(I+3*(L-1))/BOHR DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM WRITE(6,*)' AAT shifted back to laboratory system' OPEN(15,FILE='FILEU.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,*)' FILEU.TEN written ...' RETURN END SUBROUTINE WRITETTT(n,ALPHA,A,G,s) 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) character*(*)s OPEN(2,FILE=s) NAT=n/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,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