program cfcred implicit none integer*4 i,nat,natn,n,ix,iargc,ii,nn,nmol,nmoll,jj,iin,imol integer*4,allocatable::q(:),ind(:),iout(:),alist(:,:),natm(:) real*8,allocatable:: r(:),f(:,:),p(:,:,:),aa(:,:,:) real*8 e(3),a,b,c logical lq character*10 s10 character*1 am write(6,600) 600 format(' Reduces size of a vibrating crystal to ellipsoid',/, 1 ' Input: FILE.FC',/, 1 ' FILE.X',/, 1 ' FILE.TEN',/, 1 ' Output: FILEN.FC',/, 1 ' FILEN.X',/, 1 ' FILEN.TEN',/,/) if(iargc().ne.4)then write(6,*)'Usage: cfcred A B C am' write(6,*)' A, B, C ellipsoid axes in A' write(6,*)' am = a: delete outer atoms' write(6,*)' am = m: delete whole molecules' stop endif call getarg(1,s10) read(s10,*)A call getarg(2,s10) read(s10,*)B call getarg(3,s10) read(s10,*)C call getarg(4,am) if(am.eq.'a')am='A' open(9,file='FILE.X',status='old') read(9,*) read(9,*)nat n=3*nat e(1)=0.0d0 e(2)=0.0d0 e(3)=0.0d0 allocate(r(n),q(nat),f(n,n),ind(nat),p(nat,3,3),aa(nat,3,3)) do 1 i=1,nat read(9,*)q(i),(r(ix+3*(i-1)),ix=1,3) do 1 ix=1,3 1 e(ix)=e(ix)+r(ix+3*(i-1)) do 11 ix=1,3 11 e(ix)=e(ix)/dble(nat) close(9) write(6,*)nat,' atoms in FILE.X' natn=0 if(am.ne.'A')then allocate(iout(nat)) call dom(nmol,iout,nat,r,q) allocate(alist(nmol,nat),natm(nmol)) call dolist(nmol,iout,nat,alist,natm) nmoll=0 do 21 imol=1,nmol iin=1 do 211 jj=1,natm(imol) i=alist(imol,jj) 211 if(.not.lq(r(1+3*(i-1))-e(1),r(2+3*(i-1))-e(2),r(3+3*(i-1))-e(3), 1 A,B,C))iin=0 if(iin.eq.1)then nmoll=nmoll+1 do 212 jj=1,natm(imol) natn=natn+1 212 ind(natn)=alist(imol,jj) endif 21 continue write(6,*)nmoll,' molecules left' else do 2 i=1,nat if(lq(r(1+3*(i-1))-e(1),r(2+3*(i-1))-e(2),r(3+3*(i-1))-e(3), 1 A,B,C))then natn=natn+1 ind(natn)=i endif 2 continue endif nn=3*natn open(9,file='FILEN.X') write(9,*)'FILEN' write(9,*)natn do 3 ii=1,natn i=ind(ii) 3 write(9,900)q(i),(r(ix+3*(i-1)),ix=1,3) 900 format(i5,3f12.6,' 0 0 0 0 0 0 0 0.0') close(9) write(6,*)natn,' atoms left in FILEN.X' call READFF(n,f) call WRITEFF(N,f,ind,nn) CALL READTEN(nat,P,AA,.false.,r) call WRITETEN(P,AA,r,nat,natn,IND) end subroutine dolist(nmol,iout,nat,alist,natm) implicit none integer*4 nmol,nat,alist(nmol,nat),natm(nmol),i,iout(*),imol natm=0 do 1 i=1,nat natm(iout(i))=natm(iout(i))+1 1 alist(iout(i),natm(iout(i)))=i do 9 imol=1,nmol write(6,*)' molecule ',imol,': ',natm(imol),' atoms' 9 write(6,603)(alist(imol,i),i=1,natm(imol)) 603 format(15i5) c natm .. number of atoms in a molecule c alist .. list of atoms in a molecule write(6,*) return end subroutine dom(nmol,iout,nat,r,ik) implicit none integer*4 iout(nat),nmol,nat,nc,ic,i,ik(*),j real*8 r(*) real*4 bonding(88),rfac,rb,rt,xi,yi,zi,xj,yj,zj data bonding/ 10.50,0.32,0.98,1.28,0.95,0.87,0.82,0.80,0.78,0.77, 10.76,1.59,1.41,1.23,1.16,1.11,1.20,1.04,1.03,2.08,1.79,1.49,1.37, 11.27,1.23,1.22,1.22,1.21,1.20,1.22,1.30,1.31,1.27,1.25,1.21,1.19, 11.17,2.21,1.96,1.67,1.50,1.39,1.35,1.32,1.30,1.30,1.33,1.39,1.53, 11.49,1.46,1.45,1.41,1.38,1.36,2.40,2.03,1.74,1.70,1.70,1.69,1.68, 11.67,1.90,1.66,1.64,1.64,1.63,1.62,1.61,1.75,1.61,1.49,1.39,1.35, 11.33,1.31,1.32,1.35,1.39,1.54,1.53,1.52,1.51,1.51,1.50,1.49,0.10/ rfac=1.1 iout=0 c find free atom, and connected atoms nmol=1 iout(1)=1 991 nc=1 99 ic=0 do 13 i=1,nat c is atom i close to some in the molecule?: if(iout(i).eq.0)then xi=real(r(3*(i-1)+1)) yi=real(r(3*(i-1)+2)) zi=real(r(3*(i-1)+3)) rb=bonding(ik(i)+1)*rfac do 2 j=1,nat if(iout(j).eq.nmol)then rt=(rb+bonding(ik(j)+1)*rfac)**2 xj=real(r(3*(j-1)+1)) yj=real(r(3*(j-1)+2)) zj=real(r(3*(j-1)+3)) if((xi-xj)**2+(yi-yj)**2+(zi-zj)**2.lt.rt)then c atom j is bound to i, which is close=take i, too: iout(i)=nmol ic=ic+1 nc=nc+1 goto 13 endif endif 2 continue endif 13 continue if(ic.ne.0)goto 99 do 131 i=1,nat if(iout(i).eq.0)then c some atoms is left, start new molecule: nmol=nmol+1 iout(i)=nmol goto 991 endif 131 continue write(6,*)nmol,' molecules' return end SUBROUTINE WRITETEN(P,A,X,NAT,NATN,IND) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(NAT,3,3),A(NAT,3,3),X(*),IND(*) real*8 ,allocatable::R(:,:),AB(:,:,:) allocate(R(3,NATN)) BOHR=0.52917705993d0 DO 3 Ln=1,NATN L=IND(Ln) DO 3 I=1,3 3 R(I,Ln)=X(I+3*(L-1))/BOHR allocate(AB(NATN,3,3)) DO 2201 Ln=1,NATN L=IND(Ln) DO 2201 J=1,3 DO 2201 I=1,3 AB(Ln,J,I)=A(L,J,I) DO 2201 ID=1,3 DO 2201 IG=1,3 2201 AB(Ln,J,I)=AB(Ln,J,I)+0.25d0*EPS(I,IG,ID)*R(IG,Ln)*P(L,J,ID) WRITE(6,*)' AAT shifted back to laboratory system' OPEN(15,FILE='FILEN.TEN') WRITE(15,1500) NATN,NATN-6,0 1500 FORMAT(3I5) DO 10 Ln=1,NATN L=IND(Ln) 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,NATN DO 220 J=1,3 220 WRITE(15,1501) (AB(L,J,I),I=1,3),L Z0=0.0d0 DO 230 L=1,NATN DO 230 J=1,3 230 WRITE(15,1501) (Z0,I=1,3),L DO 100 Ln=1,NATN L=IND(Ln) DO 100 J=1,3 100 WRITE(15,1501) (P(L,J,I),I=1,3),L CLOSE(15) WRITE(6,*)' FILEN.TEN written ...' 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), 1C(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 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 C SUBROUTINE READTEN(NAT,P,A,LAPT,X) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(NAT,3,3),A(NAT,3,3),X(*) LOGICAL LAPT real*8,allocatable::AI(:,:,:),AJ(:,:,:),PV(:,:,:),R(:,:) 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 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) WRITE(6,*)' Dipole derivatives read - in ' C C If no axial tensor, then construct its polar part only: IF(LAPT)THEN DO 2203 L=1,NAT DO 2203 J=1,3 DO 2203 I=1,3 SUM=0.0d0 DO 2204 ID=1,3 DO 2204 IG=1,3 2204 SUM=SUM+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)=SUM 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) WRITE(6,*)' VCD parameters read in ...' 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 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 to atomic origin' CLOSE(15) WRITE(6,*)' FILE.TEN read ...' RETURN END C SUBROUTINE WRITEFF(N,FCAR,ind,nn) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N),ind(*) OPEN(20,FILE='FILEN.FC') N1=1 1 N3=N1+4 IF(N3.GT.nn)N3=nn DO 130 LN=N1,nn 130 WRITE(20,17)LN,(FCAR(3*(ind((LN+2)/3)-1)+LN-3*((LN-1)/3), 13*(ind((J+2)/3)-1)+J-3*((J-1)/3)),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.nn)GOTO 1 17 FORMAT(I4,5D14.6) close(20) write(6,*)'FILEN.FC written' RETURN END function lq(x,y,z,a,b,c) implicit none real*8 x,y,z,a,b,c logical lq lq=(x/a)**2+(y/b)**2+(z/c)**2.lt.1.0d0 return end SUBROUTINE READFF(N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) logical lex inquire(FILE='FILE.FC',exist=lex) if(.not.lex)then write(6,*)'FILE.FC not found, stopped' stop endif 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