program crystaltdc implicit none integer*4 nat,N,i,ix,jx,ii,j,jj,natb,n27,iargc,ic,iu,jc,kc, 1NI,NJ,NK,ice,iccub,ib,jb,Nb,jstart real*8 d1,d2,d5,uij,uir,ujr,epsr,vi(3),vj(3),vk(3), 1AI,AJ,AK,xi,yi,zi,xij,yij,zij real*8,allocatable::r(:),p(:,:),f(:,:),fr(:,:),fn(:,:),rb(:), 1fc(:,:) integer*4,allocatable::ijb(:,:) logical lref character*10 s10 write(6,800) 800 format(/,' Corrects crystal cell FILE.FC by TDC interactions',/, 1 ' based on FILE.TEN for atomic pairs not in CCT.INP',/,/, 1' Input: ',/, 1' Cell variables:',/, 1' ---------------',/, 1' FILE0.FC',/, 1' FILE.TEN',/, 1' FILE.X ',/, 1' CELL.TXT',/, 1' Source FF from:',/, 1' ---------------',/, 1' CRYSTAL.X ',/, 1' CCT.INP',/, 1' Output: FILE.FC',/, 1' REF.TXT (if REF.FC present)',/) if(iargc().ne.4)then write(6,*)'Usage: crystaltdc eps NI NJ NK' stop endif call getarg(1,s10) read(s10,*)epsr call getarg(2,s10) read(s10,*)NI call getarg(3,s10) read(s10,*)NJ call getarg(4,s10) read(s10,*)NK write(6,609)epsr,NI,NJ,NK 609 format(/,' Epsr = ',f6.2,' cell ',i2,' x ',i2,' x ',i2,/) allocate(r(1)) call rdx(0,nat,r,'FILE.X') deallocate(r) N=3*nat allocate(r(N),p(N,3),F(N,N),fr(N,N),fn(N,N), 1fc(N,N)) call rdx(1,nat,r,'FILE.X') allocate(rb(1)) call rdx(0,natb,rb,'CRYSTAL.X') deallocate(rb) Nb=3*natb allocate(rb(Nb),ijb(Nb,Nb)) call rdx(1,natb,rb,'CRYSTAL.X') n27=natb/nat write(6,*)n27,' cells in CRYSTAL.X' if(mod(natb,nat).ne.0)call report('non-integer ratio') call vijk(vi,vj,vk) open(9,file='FILE.TEN',status='old') read(9,*) do 1 i=1,N 1 read(9,*)(p(i,j),j=1,3) close(9) write(6,*)'ATP read' CALL READFF(N,F,'FILE0.FC') inquire(file='REF.FC',exist=lref) if(lref)then CALL READFF(N,fr,'REF.FC') open(20,file='REF.TXT') endif call vz2(fc,N,N) c TDC FF for cell 000 from all the others: do 9 IC=-NI,NI write(6,*)IC,' of ',NI AI=dble(IC) do 9 JC=-NJ,NJ AJ=dble(JC) do 9 KC=-NK,NK AK=dble(KC) if(IC.ne.0.or.JC.ne.0.or.KC.ne.0)then do 10 i=1,nat xi=r(1+3*(i-1)) yi=r(2+3*(i-1)) zi=r(3+3*(i-1)) do 10 j=i,nat xij=xi-r(1+3*(j-1))-AI*vi(1)-AJ*vj(1)-AK*vk(1) yij=yi-r(2+3*(j-1))-AI*vi(2)-AJ*vj(2)-AK*vk(2) zij=zi-r(3+3*(j-1))-AI*vi(3)-AJ*vj(3)-AK*vk(3) d2=xij**2+yij**2+zij**2 d1=dsqrt(d2) d5=d2*d2*d1 jstart=1 if(i.eq.j)jstart=ix do 10 ix=1,3 do 10 jx=jstart,3 ii=ix+3*(i-1) jj=jx+3*(j-1) uij=p(ii,1)*p(jj,1)+p(ii,2)*p(jj,2)+p(ii,3)*p(jj,3) uir=p(ii,1)* xij +p(ii,2)* yij +p(ii,3)* zij ujr=p(jj,1)* xij +p(jj,2)* yij +p(jj,3)* zij 10 fc(ii,jj)=fc(ii,jj)+(uij*d2-3.0d0*uir*ujr)/d5 endif 9 continue fc=fc/epsr do 11 ii=1,N do 11 jj=ii+1,N 11 fc(jj,ii)=fc(ii,jj) write(6,*)'TDC FF made' c delete already made interactions: c 1. get interactions call getint(ijb,natb) c 2. determine central cube: ice=iccub(natb,rb,n27) c 3. delete ff elements: c atom i in central cube do 2 i=1,nat ib=i+nat*(ice-1) xi=rb(1+3*(ib-1)) yi=rb(2+3*(ib-1)) zi=rb(3+3*(ib-1)) do 2 j=1,nat do 2 iu=1,n27 c atom j in cube iu if(iu.ne.ice)then jb=j+nat*(iu-1) if(ijb(ib,jb).ne.0)then c take this interaction out of TDC: xij=xi-rb(1+3*(jb-1)) yij=yi-rb(2+3*(jb-1)) zij=zi-rb(3+3*(jb-1)) d2=xij**2+yij**2+zij**2 d1=dsqrt(d2) d5=d2*d2*d1 do 12 ix=1,3 ii=ix+3*(i-1) do 12 jx=1,3 jj=jx+3*(j-1) uij=p(ii,1)*p(jj,1)+p(ii,2)*p(jj,2)+p(ii,3)*p(jj,3) uir=p(ii,1)* xij +p(ii,2)* yij +p(ii,3)* zij ujr=p(jj,1)* xij +p(jj,2)* yij +p(jj,3)* zij 12 fc(ii,jj)=fc(ii,jj)-(uij*d2-3.0d0*uir*ujr)/d5/epsr endif endif 2 continue write(6,*)'Established elements deleted' c add TDC contribution to orifinal FF: fn=F+fc if(lref)then do 3 i=1,nat do 3 ix=1,3 do 3 j=i,nat do 3 jx=1,3 ii=ix+3*(i-1) jj=jx+3*(j-1) 3 write(20,200)i,j,ix,jx,f(ii,jj),fn(ii,jj),fr(ii,jj) 200 format(4i4,3g12.4) endif if(lref)then close(20) write(6,*)'REF.TXT written' endif OPEN(20,FILE='FILE.FC') CALL WRITEFF(N,N,fn) CLOSE(20) WRITE(*,*)'FF written into FILE.FC' end C SUBROUTINE READFF(N,FCAR,s) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) character*(*) s OPEN(20,FILE=s,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) CLOSE(20) write(6,*)s//' read ' RETURN END C C SUBROUTINE WRITEFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION FCAR(MX3,N) C CONST=4.359828/0.5291772**2 CONST=1.0d0 DO 6 I=1,N DO 6 J=1,N 6 FCAR(I,J)=FCAR(I,J)/CONST 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 rdx(ic,nat,r,f) implicit none integer*4 ic,i,nat,N,ix real*8 r(*),bohr character*(*) f bohr=0.529177d0 open(9,file=f,status='old') read(9,*) read(9,*)nat if(ic.eq.0)goto 99 N=3*nat do 1 i=1,nat 1 read(9,*)r(1+3*(i-1)),(r(ix+3*(i-1)),ix=1,3) do 3 i=1,N 3 r(i)=r(i)/bohr write(6,*)'Geometry read,',nat,' atoms' 99 close(9) return end subroutine vz2(v,n,m) implicit none real*8 v(n,m) integer*4 n,m,i,j do 1 i=1,n do 1 j=1,m 1 v(i,j)=0.0d0 return end subroutine getint(ij,nat) implicit none integer*4 i,j,nat,no,io,n1,n2,ij(nat,nat) integer*4,allocatable::ind(:) logical lnames,old character*10 s10 allocate(ind(nat)) do 8 i=1,nat do 8 j=1,nat 8 ij(i,j)=0 no=0 lnames=.false. old=.true. write(6,*)'reading CCT.INP' open(9,file='CCT.INP',status='old') 4 read(9,90,end=99,err=99)s10 90 format(a10) if(s10(1:6).eq.'LNAMES')read(9,*)lnames if(s10(1:6).eq.'OLDFOR')read(9,*)old if(s10(1:7).eq.'POLYMER')then read(9,*)no do 5 io=1,no read(9,*) if(lnames)read(9,*) if(old)then read(9,91)(ind(i),i=1,nat) read(9,91)(ind(i),i=1,nat) 91 format(20i3) else n1=1 6 n2=min(n1+19,nat) read(9,*)(ind(i),i=n1,n2) if(n2.lt.nat)then n1=n2+1 goto 6 endif n1=1 7 n2=min(n1+19,nat) read(9,*)(ind(i),i=n1,n2) if(n2.lt.nat)then n1=n2+1 goto 7 endif endif do 5 i=1,nat do 5 j=1,nat 5 if(ind(i).ne.0.and.ind(j).ne.0)ij(i,j)=ij(i,j)+1 goto 99 endif goto 4 99 close(9) write(6,*)'CCT.INP ',no,' overlaps' 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(*),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(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)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 subroutine report(s) character*(*) s write(6,*)s stop end subroutine vijk(vi,vj,vk) implicit none integer*4 i real*8 vi(3),vj(3),vk(3),bohr bohr=0.529177d0 open(9,file='CELL.TXT',status='old') read(9,*)vi read(9,*)vj read(9,*)vk close(9) do 1 i=1,3 vi(i)=vi(i)/bohr vj(i)=vj(i)/bohr 1 vk(i)=vk(i)/bohr return end