PROGRAM CONTACTCC IMPLICIT none integer*4 ik,iout,i,j,ic,nat,nc,alist,natm,ia,iks,ixx,iyy,izz, 1imol,jj,nmol,istart,ind,nat1,iu,ice,icx,icy,icz,no,ii,iku, 1nx,ny,nz,lmo,mmoin,nain,ns,ju,js,iup,iue,nss,nua,bt real*4 bonding(88),r,xi,xj,yi,yj,zi,zj,rt,rb,ce(3),rs, 1tol,px,py,pz,btol,cet(3),am(3),cmin(3),cmax(3),aind, 1tx,ty,tz,d2,ru allocatable ::ik(:),r(:,:),iout(:),alist(:,:),natm(:), 1rs(:,:),ind(:,:),lmo(:),iks(:),aind(:,:),iku(:),ru(:,:), 1bt(:,:) 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/ logical lwr,ltab,bb character*10 s10 tol=0.001 c WRITE(6,3000) 3000 FORMAT(' Takes "crystal" geometry CRYSTAL.X, make',/, 1 ' clusters by shifts of the central cube',/ 1 ' Input: CRYSTAL.X',/, 1 ' CELL.TXT',/, 1 ' CONTACTCC.OPT',/, 3 ' Output: cluster geometries, CCT.TEM',/,/) call readopt(px,py,pz,lwr,btol,am,nx,ny,nz,no,ltab) open(2,file='CRYSTAL.X') read(2,*) read(2,*)nat allocate(ik(nat),r(nat,3),iout(nat),ind(nat,2),lmo(nat),iks(nat), 1rs(nat,3),bt(nat,7)) do 3 i=1,nat iout(i)=0 if(ltab)then read(2,*)ik(i),(r(i,j),j=1,3),(bt(i,j),j=1,7) else read(2,*)ik(i),(r(i,j),j=1,3) endif 3 if(ik(i)+1.lt.1.or.ik(i)+1.gt.88)call report('type out of range') close(2) write(6,*)nat,' atoms in CRYSTAL.X' call iccub(ice,nat,r,cet,lwr) nat1=nat/27 write(6,602) 602 format(' Finding molecules') c find free atom, and connected atoms, in the first cube actually: nmol=1 iout(1)=1 991 nc=1 99 ic=0 do 13 i=1,nat1 c is atom i close to some in the molecule?: if(iout(i).eq.0)then xi=r(i,1) yi=r(i,2) zi=r(i,3) rb=bonding(ik(i)+1)*btol do 2 j=1,nat1 if(iout(j).eq.nmol)then rt=(rb+btol*bonding(ik(j)+1))**2 xj=r(j,1)-xi yj=r(j,2)-yi zj=r(j,3)-zi if(ltab)then bb=bt(i,1).eq.j.or.bt(i,2).eq.j.or.bt(i,3).eq.j.or. 1 bt(i,4).eq.j.or.bt(i,5).eq.j.or.bt(i,6).eq.j.or.bt(i,7).eq.j else bb=xj**2+yj**2+zj**2.lt.rt endif if(bb)then c atom i is bound to j, which is in, take i, too: iout(i)=nmol ic=ic+1 nc=nc+1 goto 13 endif endif 2 continue endif 13 continue c if still adding to nmol, continue: if(ic.ne.0)goto 99 do 131 i=1,nat1 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,606)nmol 606 format(i5,' molecules found:') c alist .. atom list related to the first cube!!! allocate(alist(nmol,nat),natm(nmol)) do 8 imol=1,nmol 8 natm(imol)=0 do 1 i=1,nat1 natm( iout(i))=natm(iout(i))+1 1 alist(iout(i),natm(iout(i)))=i if(lwr)then do 9 imol=1,nmol write(6,604)imol 604 Format(i4,':',$) 9 write(6,603)(alist(imol,jj),jj=1,natm(imol)) 603 format(15i5) endif c natm .. number of atoms in a molecule c alist .. list of atoms in a molecule c offset to the central cube: alist=alist+(ice-1)*nat1 if(lwr)then write(6,*)'Offset to central cube:' do 91 imol=1,nmol write(6,604)imol 91 write(6,603)(alist(imol,jj),jj=1,natm(imol)) endif c write final atomic map: open(19,file='CCT.TEM') write(19,901)no 901 format('LNAMES',/,'T',/,'POLYMER',/,I9) iu=0 ju=0 js=0 allocate(aind(no,3)) c loop over cube shifts do 101 ixx=-nx,nx cmin(1)=cet(1)+px*am(1)*real(ixx)-px/2.0 cmax(1)=cet(1)+px*am(1)*real(ixx)+px/2.0 do 101 iyy=-ny,ny cmin(2)=cet(2)+py*am(2)*real(iyy)-py/2.0 cmax(2)=cet(2)+py*am(2)*real(iyy)+py/2.0 do 101 izz=-nz,nz cmin(3)=cet(3)+pz*am(3)*real(izz)-pz/2.0 cmax(3)=cet(3)+pz*am(3)*real(izz)+pz/2.0 iu=iu+1 aind(iu,1)=am(1)*real(ixx) aind(iu,2)=am(2)*real(iyy) aind(iu,3)=am(3)*real(izz) c does it already exist? iue=0 do 1013 iup=1,iu-1 tx=aind(iup,1)-aind(iu,1) ty=aind(iup,2)-aind(iu,2) tz=aind(iup,3)-aind(iu,3) if(abs(real(nint(tx))-tx).lt.tol.and. 1 abs(real(nint(ty))-ty).lt.tol.and. 1 abs(real(nint(tz))-tz).lt.tol)then iue=iup goto 1014 endif 1013 continue c number of molecules inside: 1014 mmoin=0 c number of atoms inside: nain=0 ia=0 c loop over molecules do 1011 imol=1,nmol c loop over cubes do 1011 icx=-1,1 do 1011 icy=-1,1 do 1011 icz=-1,1 c center of this molecule in its cube: ce=0. do 102 i=1,natm(imol) ce(1)=ce(1)+r(alist(imol,i),1)+px*real(icx) ce(2)=ce(2)+r(alist(imol,i),2)+py*real(icy) 102 ce(3)=ce(3)+r(alist(imol,i),3)+pz*real(icz) ce=ce/real(natm(imol)) c is it in? if(ce(1).gt.cmin(1).and.ce(1).le.cmax(1).and. 1 ce(2).gt.cmin(2).and.ce(2).le.cmax(2).and. 1 ce(3).gt.cmin(3).and.ce(3).le.cmax(3))then mmoin=mmoin+1 lmo(mmoin)=imol nain=nain+natm(imol) do 1021 i=1,natm(imol) ia=ia+1 iks(ia)=ik(alist(imol,i)) rs(ia,1)=r(alist(imol,i),1)+px*real(icx) rs(ia,2)=r(alist(imol,i),2)+py*real(icy) 1021 rs(ia,3)=r(alist(imol,i),3)+pz*real(icz) endif 1011 continue ns=ia if(lwr)then write(6,644)iu,ixx,iyy,izz,ns,mmoin 644 format(/,' Cube',i6,3i3,',',i5,' atoms,',i5,' molecules:',$) do 104 jj=1,mmoin 104 write(6,605)lmo(jj) 605 format(i5,$) write(6,*) write(6,645)cmin(1),cmax(1),cmin(2),cmax(2), 1 cmin(3),cmax(3) 645 format('Min-max xyz:',3(f6.2,' - ',f6.2,1x)) endif call sets(s10,iu,istart) write(19,400)iu if(iue.eq.0)then if(lwr)write(6,*)'unique' open(40,file='u.'//s10(istart:len(s10))//'.x') write(19,'(a)')'u.'//s10(istart:len(s10)) write(40,'(a)')'u.'//s10(istart:len(s10)) ju=ju+1 else if(lwr)write(6,655)iue 655 format(' same as',i5) open(40,file='s.'//s10(istart:len(s10))//'.x') write(40,'(a)')'s.'//s10(istart:len(s10)) call sets(s10,iue,istart) write(19,'(a)')'u.'//s10(istart:len(s10)) js=js+1 nss=ns c load in the unique geometry again: open(41,file='u.'//s10(istart:len(s10))//'.x') read(41,*) read(41,*)nua allocate(iku(nua),ru(nua,3)) do 1017 i=1,nua 1017 read(41,*)iku(i),(ru(i,j),j=1,3) close(41) c check that all atoms are also in the unique one: 1022 do 1018 i=1,ns xi=rs(i,1) yi=rs(i,2) zi=rs(i,3) do 1019 j=1,nua xj=xi-ru(j,1)+tx*px yj=yi-ru(j,2)+ty*py zj=zi-ru(j,3)+tz*pz 1019 if(xj**2+yj**2+zj**2.lt.tol.and.iks(i).eq.iku(j))goto 1018 c atom i is missing, delete it: do 1020 j=i,ns-1 iks(j)=iks( j+1) rs( j,1)=rs(j+1,1) rs( j,2)=rs(j+1,2) 1020 rs( j,3)=rs(j+1,3) ns=ns-1 if(ns.gt.0)goto 1022 1018 continue deallocate(iku,ru) if(lwr)then if(ns.eq.nss)then write(6,607) 607 format(' All atoms kept') else write(6,608)ns,nss 608 format(i6,' atoms kept of',i6) endif endif endif write(40,400)ns c atomic index ii=0 do 1012 i=1,ns xi=rs(i,1) yi=rs(i,2) zi=rs(i,3) do 1015 j=1,nat d2=(r(j,1)-xi)**2+(r(j,2)-yi)**2+(r(j,3)-zi)**2 if(d2.lt.tol.and.iks(i).eq.ik(j))then ii=ii+1 ind(ii,1)=j ind(ii,2)=i goto 1012 endif 1015 continue 1012 write(40,400)iks(i),xi,yi,zi 400 format(i6,3f12.6,' 0 0 0 0 0 0 0 0.0') write(19,191)ii 191 format(i5,$) do 1016 i=1,ii 1016 write(19,192)ind(i,1),ind(i,2) 192 format(1x,2i5,$) write(19,*) 101 close(40) close(19) write(6,3434)ju,js 3434 format(/,i6,' unique,',i5,' shifted,',/, 1'CCT.TEM written') end subroutine report(s) character*(*) s write(6,*)s stop end subroutine sets(s10,i,istart) implicit none character*10 s10 integer*4 i,istart write(s10,5002)i 5002 format(i10) do 504 istart=1,len(s10) 504 if(s10(istart:istart).ne.' ')goto 505 505 return end subroutine readopt(px,py,pz,lwr,btol,am,nx,ny,nz,no,ltab) implicit none real*4 px,py,pz,btol,am(3) logical lwr,lex,ltab integer*4 nx,ny,nz,no character*4 key c bonding cutoff multiplied by btol, to accomodate for MD distorted c bonds btol=1.1 c use bond table from CRYSTAL.X: ltab=.false. lwr=.true. px=0. py=0. pz=0. am=0.5 nx=2 ny=2 nz=2 inquire(file='CELL.TXT',exist=lex) if(lex)then open(2,file='CELL.TXT') read(2,*)px read(2,*)py,py read(2,*)pz,pz,pz close(2) endif inquire(file='CONTACTCC.OPT',exist=lex) if(lex)then open(2,file='CONTACTCC.OPT') 100 read(2,2000,end=2200,err=2200)key 2000 format(a4) if(key.eq.'CELL')read(2,*)px,py,pz if(key.eq.'NXYZ')read(2,*)nx,ny,nz if(key.eq.'SHIF')read(2,*)am if(key.eq.'WRIT')read(2,*)lwr if(key.eq.'BTOL')read(2,*)btol if(key.eq.'TABL')read(2,*)ltab goto 100 2200 close(2) endif if(px.lt.0.1)call report('PX must be set') if(py.lt.0.1)call report('PY must be set') if(pz.lt.0.1)call report('PZ must be set') no=(2*nx+1)*(2*ny+1)*(2*nz+1) write(6,600)no,px,py,pz,nx,ny,nz,am,lwr,btol 600 format(i6,' overlaps',/, 1' CELL :',3f9.2,5x,' NXYZ:',3i3,/, 1' SHIFTS:',3f9.2,5x,' WRIT:',l3,/, 1' BTOL :',f10.4,/) return end subroutine iccub(ice,nat,r,cet,lwr) implicit none integer*4 nat,nat1,ix,iu,ii,ice,i real*4 cet(3),cei(27,3),r(nat,3),xi,yi,zi logical lwr if(mod(nat,27).ne.0)call report('mod nat 27 <> 0') nat1=nat/27 c determine central cube do 40 ix=1,3 cet(ix)=0.0 do 42 iu=1,27 cei(iu,ix)=0.0 ii=nat1*(iu-1) do 41 i=ii+1,ii+nat1 cet( ix)=cet( ix)+r(i,ix) 41 cei(iu,ix)=cei(iu,ix)+r(i,ix) 42 cei(iu,ix)=cei(iu,ix)/real(nat1) 40 cet(ix)=cet(ix)/real(nat) if(lwr)then write(6,600)cet 600 format(/,' CRYSTAL.X center (A):',7x,3f7.2,/, 1 ' Cube centers (A):') write(6,601)(iu,(cei(iu,ix),ix=1,3),iu=1,27) 601 format(3(i3,':',3f7.2)) endif ice=0 do 43 iu=1,27 xi=cei(iu,1)-cet(1) yi=cei(iu,2)-cet(2) zi=cei(iu,3)-cet(3) 43 if(sqrt(xi**2+yi**2+zi**2).lt.0.1)ice=iu if(ice.eq.0)then call report('Central cube not found') else write(6,602)ice 602 format(' Central cube:',i3,', finding molecules') endif return end