program environ implicit none integer*4 nat,i,ice,iccub,nat1,nmol,nc,ic,imol,j,jj, 1natc,tt(88) real*4 xi,yi,zi,xj,yj,zj,rb,btol,bonding(88),d,rt, 1dij2,cei(27,3),d2 integer*4 ,allocatable::ik(:),iout(:),natm(:),alist(:,:), 1it(:,:),bt(:,:) real*4, allocatable::r(:,:),cm(:,:) logical bb,ltab character*8 s8 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/ WRITE(6,3000) 3000 FORMAT(' Takes "crystal" geometry CRYSTAL.X, selects',/, 1 ' central cube and close molecules up to a distance',/,/, 1 ' Input: CRYSTAL.X',/, 3 ' Output: EN.X',/,/) if(iargc().lt.1)then write(6,*)'Usage: environ <> <>' stop endif call getarg(1,s8) read(s8,*)d d2=d**2 btol=1. if(iargc().gt.1)then call getarg(2,s8) read(s8,*)btol endif ltab=.false. if(iargc().gt.2)then call getarg(3,s8) read(s8,*)ltab endif open(2,file='CRYSTAL.X',status='old') read(2,*) read(2,*)nat allocate(ik(nat),r(nat,3),iout(nat),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' cei=0. ice=iccub(nat,r,cei) nat1=nat/27 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 dij2=xj**2+yj**2+zj**2 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=dij2.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, atom lists:') 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 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) c natm .. number of atoms in a molecule c alist .. list of atoms in a molecule write(6,*) c centers of masses in the first cube: allocate(cm(nmol,3)) cm=0. do 12 imol=1,nmol do 12 i=1,3 do 5 jj=1,natm(imol) 5 cm(imol,i)=cm(imol,i)+r(alist(imol,jj),i) 12 cm(imol,i)=cm(imol,i)/real(natm(imol)) write(6,*)'Centers in the first cube:' do 15 imol=1,nmol 15 write(6,612)imol,(cm(imol,i),i=1,3) 612 format(i4,3f10.3) allocate(it(27,nmol)) it=0 c number of molecules closer than d to central: nc=0 c number of atoms in extended system: natc=nat1 c loop over all cubes do 6 ic=1,27 if(ic.ne.ice)then c loop over all molecules: do 4 imol=1,nmol c center of mass, absolute: xi=cm(imol,1)+cei(ic,1)-cei(1,1) yi=cm(imol,2)+cei(ic,2)-cei(1,2) zi=cm(imol,3)+cei(ic,3)-cei(1,3) c loop over atoms in the central cube: do 7 j=nat1*(ice-1)+1,nat1*ice xj=r(j,1)-xi yj=r(j,2)-yi zj=r(j,3)-zi dij2=xj**2+yj**2+zj**2 if(dij2.lt.d2)then nc=nc+1 natc=natc+natm(imol) it(ic,imol)=1 goto 4 endif 7 continue 4 continue endif 6 continue tt=0 open(9,file='EN.X') write(9,900)natc 900 format('Central cube and environment',/,i6) do 11 i=nat1*(ice-1)+1,nat1*ice tt(ik(i))=tt(ik(i))+1 11 write(9,901)ik(i),(r(i,j),j=1,3) 901 format(i6,3f12.6,' 0 0 0 0 0 0 0 0.0') do 14 ic=1,27 do 14 imol=1,nmol if(it(ic,imol).eq.1)then do 10 jj=1,natm(imol) xi=r(alist(imol,jj),1)+cei(ic,1)-cei(1,1) yi=r(alist(imol,jj),2)+cei(ic,2)-cei(1,2) zi=r(alist(imol,jj),3)+cei(ic,3)-cei(1,3) tt(ik(alist(imol,jj)))=tt(ik(alist(imol,jj)))+1 10 write(9,901)ik(alist(imol,jj)),xi,yi,zi endif 14 continue close(9) write(6,*)nc,' molecules of ',natc,' atoms selected to EN.X' do 16 i=1,88 16 if(tt(i).ne.0)write(6,616)tt(i),i 616 format(i4,' x ',i2) end function iccub(nat,r,cei) implicit none integer*4 nat,nat1,ix,iu,ii,iccub,ice,i real*4 cet(3),cei(27,3),r(nat,3),xi,yi,zi 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) 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,27) 601 format(3(i3,':',3f7.2)) 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') iccub=ice endif return end subroutine report(s) character*(*) s write(6,*)s stop end