PROGRAM cvib IMPLICIT none integer*4 nat,n,ic,m,nc real*8 u(3,3) real*8 ,allocatable::r(:),c(:) integer*4 ,allocatable::ty(:),bt(:,:),nt(:),tc(:),bc(:,:),ntc(:) logical lcorner,ledge c WRITE(6,3000) 3000 FORMAT(' Takes "crystal" geometry CRYSTAL.X,',/, 1 ' makes pairs of CELL.X geometries',/ 1 ' to calculate its force field.',/, 1 ' Input: CRYSTAL.X ',/, 1 ' CELL.X ',/, 1 ' CELL.TXT',/, 3 ' Output: cluster geometries, TEM.CCT',/,/) call ropt(ic,lcorner,ledge) call rdim(nat,n,nc,m) allocate(ty(nat),r(n),bt(nat,7),nt(nat)) allocate(tc(nc), c(m),bc(nc ,7),ntc(nc)) call rc(nat,ty,r,'CRYSTAL.X') call mkbonds(nat,ty,r,bt,nt) call rc(nc ,tc,c,'CELL.X' ) call mkbonds(nc,tc,c,bc,ntc) call capc(nat,ty,r,bt,nt ,ic,'CRYSTAL-CAP.X') call capc(nc ,tc,c,bc,ntc,ic,'CELL-CAP.x') call rv(u) call mult(lcorner,ledge,u,ic) call mult1(u) end subroutine mult1(u) IMPLICIT none integer*4 nc,m,i,j,k,ia,io,nu, 1nat,iu,ix,iov real*8 u(3,3) real*8,allocatable:: rt(:),c(:),r(:) integer*4 ,allocatable::tt(:),bt(:,:),nt(:),ind(:), 1bb(:,:),tc(:),ty(:) open(8,file='CRYSTAL-CAP.X') read(8,*) read(8,*)nat close(8) open(8,file='CELL-CAP.x') read(8,*) read(8,*)nc close(8) m=3*nc allocate(rt(m),tt(nc),nt(nc),bt(nc,7),ind(nat),c(m), 1tc(nc),r(3*nat),ty(nat)) call rc(nc ,tc,c,'CELL-CAP.x' ) call rc(nat,ty,r,'CRYSTAL-CAP.X' ) do 3 i=1,nc 3 tt(i )=tc( i) open(88,file='TEM1.CCT') write(88,80)27 80 format('LNEW',/,'t',/,'POLYMER',/,i6) iov=0 do 1 i=-1,1 do 1 j=-1,1 do 1 k=-1,1 iov=iov+1 do 2 ia=1,nc do 2 ix=1,3 2 rt(ix+3*(ia-1))=c(ix+3*(ia-1))+u(1,ix)*dble(i+1) 1 +u(2,ix)*dble(j+1) 1 +u(3,ix)*dble(k+1) call mkbonds(nc,tt,rt,bt,nt) write(88,81)iov,'CELL-CAP' 81 format(i6,/,a8) call setind(nat,ind,nc,m,r,rt,ty,tt) nu=0 do 8 iu=1,nat 8 if(ind(iu).ne.0)nu=nu+1 allocate(bb(2,nu)) nu=0 do 9 iu=1,nat if(ind(iu).ne.0)then nu=nu+1 bb(1,nu)=iu bb(2,nu)=ind(iu) endif 9 continue write(88,821)nu 821 format(i6,$) write(88,82)((bb(io,iu),io=1,2),iu=1,nu) 82 format(10(2i6,' ')) deallocate(bb) 1 continue close(88) write(6,600)27 600 format(i6,' overlaps for CELL-CAP.x') return end subroutine setind(nat,ind,nc,m,r,rt,ty,tt) implicit none integer*4 ind(nat),ia,nc,ii,ja,nat,m,ty(nat),tt(nc) real*8 x,y,z,r(3*nat),rt(m) ind=0 do 6 ia=1,nc ii=0 do 7 ja=1,nat x=r(1+3*(ja-1))-rt(1+3*(ia-1)) y=r(2+3*(ja-1))-rt(2+3*(ia-1)) z=r(3+3*(ja-1))-rt(3+3*(ia-1)) if(x**2+y**2+z**2.lt.1.0d-4)then ii=ii+1 ind(ja)=ia if(ty(ja).ne.tt(ia))then write(6,*)'small',ia,tt(ia),', big ',ja,ty(ja),' differ' stop endif endif 7 continue if(ii.eq.0)then write(6,*)'small',ia,', counterpart not found' endif if(ii.gt.1)then write(6,*)'small',ia,', big ',ja,' not unique' stop endif 6 continue return end subroutine mult(lcorner,ledge,u,ic) IMPLICIT none integer*4 nc,m,i,j,k,ia,isum,nov,io,nu,nz, 1nat,iu,ic,iov,ix,iovu real*8 u(3,3) real*8,allocatable:: rt(:),c(:),r(:) integer*4 ,allocatable::tt(:),bt(:,:),nt(:),ind(:), 1bb(:,:),bbt(:,:),tc(:),ty(:) logical lcorner,ledge,lu character*1 ch(3) character*3 nam data ch/'m','0','p'/ open(8,file='CRYSTAL-CAP.X') read(8,*) read(8,*)nat close(8) if(ic.eq.2)then open(8,file='CELL-CAP.x') else open(8,file='CELL.X') endif read(8,*) read(8,*)nc close(8) m=3*nc allocate(rt(2*m),tt(2*nc),nt(2*nc),bt(2*nc,7),ind(nat),c(m), 1tc(nc),r(3*nat),ty(nat)) if(ic.eq.2)then call rc(nc ,tc,c,'CELL-CAP.x' ) else call rc(nc ,tc,c,'CELL.X' ) endif call rc(nat,ty,r,'CRYSTAL-CAP.X' ) do 3 i=1,nc tt(i )=tc( i) tt(i+nc)=tc( i) do 3 ix=1,3 3 rt(ix+3*(i-1))=c(ix+3*(i-1))+u(1,ix)+u(2,ix)+u(3,ix) nov=6 if(lcorner)nov=nov+ 8 if(ledge )nov=nov+12 open(88,file='TEM.CCT') write(88,80)nov 80 format('LNEW',/,'t',/,'POLYMER',/,i6) iov=0 iovu=0 do 1 i=-1,1 do 1 j=-1,1 do 1 k=-1,1 nz=0 if(i.eq.0)nz=nz+1 if(j.eq.0)nz=nz+1 if(k.eq.0)nz=nz+1 do 2 ia=1,nc do 2 ix=1,3 2 rt(ix+3*(nc+ia-1))=c(ix+3*(ia-1))+u(1,ix)*dble(i+1) 1 +u(2,ix)*dble(j+1) 1 +u(3,ix)*dble(k+1) isum=abs(i)+abs(j)+abs(k) if(isum.eq.0)goto 1 if(isum.eq.3.and..not.lcorner)goto 1 if(isum.eq.2.and..not.ledge)goto 1 iov=iov+1 nam=ch(i+2)//ch(j+2)//ch(k+2) call mkbonds(2*nc,tt,rt,bt,nt) lu=.true. if(nz.eq.0)lu=i.lt.0 if(nz.eq.1)lu=(i.eq.0.and.((j.lt.k).or.(j.eq.k.and.j.lt.0))).or. 1 (j.eq.0.and.((i.lt.k).or.(i.eq.k.and.i.lt.0))).or. 1 (k.eq.0.and.((i.lt.j).or.(i.eq.j.and.i.lt.0))) if(nz.eq.2)lu=(i.ne.0.and.i.lt.0).or. 1 (j.ne.0.and.j.lt.0).or. 1 (k.ne.0.and.k.lt.0) if(lu)then call capc(2*nc,tt,rt,bt,nt,ic,nam//'.u.x') else call capc(2*nc,tt,rt,bt,nt,ic,nam//'.x') nam=ch(-i+2)//ch(-j+2)//ch(-k+2) write(6,*)ch(i+2)//ch(j+2)//ch(k+2)//' same as '//nam endif write(88,81)iov,nam//'.u' 81 format(i6,/,a5) if(lu)iovu=iovu+1 call setind(nat,ind,2*nc,2*m,r,rt,ty,tt) nu=0 do 8 iu=1,nat 8 if(ind(iu).ne.0)nu=nu+1 allocate(bb(2,nu),bbt(2,nu)) nu=0 do 9 iu=1,nat if(ind(iu).ne.0)then nu=nu+1 bb(1,nu)=iu bb(2,nu)=ind(iu) endif 9 continue write(88,821)nu 821 format(i6,$) bbt=bb if(.not.lu)then c switch elementary cells do 10 iu=1,nu if(bb(2,iu).le.nc)bbt(2,iu)=bb(2,iu)+nc 10 if(bb(2,iu).gt.nc)bbt(2,iu)=bb(2,iu)-nc endif write(88,82)((bbt(io,iu),io=1,2),iu=1,nu) 82 format(10(2i6,' ')) deallocate(bb,bbt) 1 continue close(88) write(6,600)iov,iovu 600 format(i6,' overlaps,',i6,' unique') return end subroutine capc(nat,ty,r,bt,nt,ic,f) IMPLICIT none integer*4 nat,ty(nat),bt(nat,7),nt(nat),i,j,k,ix,ic,b7(7), 1ianew,ixx,natn,ini,ica real*8 r(3*nat),vik(3),rcc,r2(3),g,ro(3),r1(3),u(3),rnc character*(*) f integer*4,allocatable::tyn(:),btn(:,:),ntn(:) real*8,allocatable::rn(:) rcc=1.519d0 rnc=1.454d0 if(ic.ne.0.and.ic.ne.1.and.ic.ne.2.and.ic.ne.3)then write(6,*)ic,' unrecognized capping type' stop endif ianew=0 ini=0 ica=0 open(9,file='scr') do 1 i=1,nat c recognize -C.=O: if(ty(i).eq.6.and.nt(i).eq.2)then do 2 ix=1,nt(i) j=bt(i,ix) if(ty(j).eq.8.and.nt(j).eq.1)then ica=ica+1 call setvk(nat,ix,bt,r,i,j,k,vik,u,g,r1,ro) ianew=ianew+1 bt(i,3)=nat+ianew nt(i)=3 c methyl if(ic.eq.0)call am(nat,r,i,g,rcc,vik,u,r1,r2,ro,ianew) c hydrogen if(ic.eq.1)call ah(nat,r,i,g,vik,u) c hydrogen replaces -C.=O if(ic.eq.2)call ahr(nat,r,i,j,k,vik,ty,nt,bt,ianew) c -NH-CH3 if(ic.eq.3)call a3(nat,i,g,vik,u,r1,r2,ianew) endif 2 continue endif c recognize -N.-H: if(ty(i).eq.7.and.nt(i).eq.2)then do 5 ix=1,nt(i) j=bt(i,ix) if(ty(j).eq.1.and.nt(j).eq.1)then ini=ini+1 call setvk(nat,ix,bt,r,i,j,k,vik,u,g,r1,ro) ianew=ianew+1 bt(i,3)=nat+ianew nt(i)=3 c methyl if(ic.eq.0)call am(nat,r,i,g,rnc,vik,u,r1,r2,ro,ianew) c hydrogen if(ic.eq.1)call ah(nat,r,i,g,vik,u) c hydrogen replaces -N.-H if(ic.eq.2)call ahr(nat,r,i,j,k,vik,ty,nt,bt,ianew) c -CO-CH3 if(ic.eq.3)call a4(nat,i,g,vik,u,r1,r2,ianew) endif 5 continue endif 1 continue close(9) allocate(tyn(nat),rn(3*nat),btn(nat,7),ntn(nat)) natn=nat tyn=ty rn=r btn=bt ntn=nt c eliminate left-over atoms: if(ic.eq.2)then 111 do 11 i=1,natn if(tyn(i).eq.0)then do 12 j=1,natn do 12 ix=1,7 if(btn(j,ix).eq.i)then do 13 ixx=ix,6 13 btn(j,ixx)=btn(j,ixx+1) btn(j,7)=0 ntn(j)=ntn(j)-1 endif 12 continue do 14 j=1,natn do 14 ix=1,7 14 if(btn(j,ix).gt.i-1)btn(j,ix)=btn(j,ix)-1 do 15 j=i,natn-1 tyn(j)=tyn(j+1) do 16 ix=1,7 16 btn(j,ix)=btn(j+1,ix) do 15 ix=1,3 15 rn(ix+3*(j-1))=rn(ix+3*(j+1-1)) natn=natn-1 goto 111 endif 11 continue write(6,*)'number of atoms ',nat,' reduced to ',natn endif if(f.eq.'CRYSTAL-CAP.X')then write(6,*)ianew,' atoms added' write(6,*)ica,' -C.O radicals' write(6,*)ini,' -N.H radicals' endif open(8,file='scr') open(9,file=f) write(9,901)natn+ianew 901 format('Capped crystal',/,i6) do 8 i=1,natn 8 write(9,902)tyn(i),(rn(ix+3*(i-1)),ix=1,3),(btn(i,ix),ix=1,7) 902 format(i6,3f12.6,7i5,' 0.0') do 9 i=1,ianew read(8,*)j,ro,b7 9 write(9,902)j,ro,b7 close(8) close(9) write(6,*)f//' written' return end subroutine ahr(nat,r,i,j,k,vik,ty,nt,bt,ianew) c hydrogen replaces i, type of j = 0 implicit none integer*4 nat,i,j,k,ty(nat),nt(nat),bt(nat,7),ianew,ix real*8 rch,r(3*nat),vik(3) ianew=ianew-1 ty(i)=1 ty(j)=0 rch=1.080d0 do 1 ix=1,3 1 r(ix+3*(i-1))=r(ix+3*(k-1))+rch*vik(ix) bt(i,3)=0 do 2 ix=1,2 2 if(bt(i,ix).eq.j)bt(i,ix)=0 nt(i)=1 return end subroutine ah(nat,r,i,g,vik,u) implicit none integer*4 nat,a,i real*8 g,r2(3),r(3*nat),rch,vik(3),u(3) rch=1.080d0 do 4 a=1,3 4 r2(a)=r(a+3*(i-1))+rch*(cos(g)*vik(a)+sin(g)*u(a)) write(9,900)6,r2,i,0,0,0,0,0,0 900 format(i6,3f12.6,7i6) return end subroutine am(nat,r,i,g,rcc,vik,u,r1,r2,ro,ianew) integer*4 i,a,ianew,nat real*8 r(*),g,rcc,vik(3),u(3),r1(3),ro(3),xh1(3),xh2(3),xh3(3), 1r2(3) do 11 a=1,3 11 r2(a)=r(a+3*(i-1))+rcc*(cos(g)*vik(a)+sin(g)*u(a)) write(9,900)6,r2,i,nat+ianew+1,nat+ianew+2,nat+ianew+3,0,0,0 900 format(i6,3f12.6,7i6) call methyl(ro,r1,r2,xh1,xh2,xh3) ianew=ianew+1 write(9,900)1,xh1,nat+ianew-1,0,0,0,0,0,0 ianew=ianew+1 write(9,900)1,xh2,nat+ianew-2,0,0,0,0,0,0 ianew=ianew+1 write(9,900)1,xh3,nat+ianew-3,0,0,0,0,0,0 return end subroutine a3(nat,i,g,vik,u,r1,r2,ianew) c add ..NH-CH3 integer*4 i,a,ianew,nat real*8 g,vik(3),u(3),r1(3),xh1(3),xh2(3),xh3(3), 1r2(3),w(3),up(3),spuw,sp,al,pi,alp,r3(3),b,bp,rnc,rcn,rnh,r4(3) pi=4.0d0*datan(1.0d0) rcn=1.369d0 rnc=1.451d0 rnh=1.008d0 al=119.0d0*pi/180.0d0 b=122.0d0*pi/180.0d0 bp=pi-b alp=pi-al do 11 a=1,3 w(a)=rcn*(cos(g)*vik(a)+sin(g)*u(a)) 11 r2(a)=r1(a)+w(a) write(9,900)7,r2,i,nat+ianew+1,nat+ianew+2,0,0,0,0 900 format(i6,3f12.6,7i6) call norm(w) spuw=sp(u,w) do 1 a=1,3 1 up(a)=u(a)-w(a)*spuw call norm(up) do 2 a=1,3 2 r3(a)=r2(a)+rnh*(w(a)*cos(alp)+up(a)*sin(alp)) ianew=ianew+1 write(9,900)1,r3,nat+ianew-1,0,0,0,0,0,0 do 3 a=1,3 3 r4(a)=r2(a)+rnc*(w(a)*cos(bp)-up(a)*sin(bp)) ianew=ianew+1 write(9,900)6,r4,nat+ianew-2,nat+ianew+1,nat+ianew+2, 1nat+ianew+3,0,0,0 call methyl(r3,r2,r4,xh1,xh2,xh3) ianew=ianew+1 write(9,900)1,xh1,nat+ianew-1,0,0,0,0,0,0 ianew=ianew+1 write(9,900)1,xh2,nat+ianew-2,0,0,0,0,0,0 ianew=ianew+1 write(9,900)1,xh3,nat+ianew-3,0,0,0,0,0,0 return end subroutine a4(nat,i,g,vik,u,r1,r2,ianew) c add ..NH-CH3 integer*4 i,a,ianew,nat real*8 g,rcc,vik(3),u(3),r1(3),xh1(3),xh2(3),xh3(3), 1r2(3),w(3),up(3),spuw,sp,al,pi,alp,r3(3),b,bp,rcn,rco,r4(3) pi=4.0d0*datan(1.0d0) rcn=1.369d0 rcc=1.522d0 rco=1.225d0 al=123.0d0*pi/180.0d0 b=115.0d0*pi/180.0d0 bp=pi-b alp=pi-al do 11 a=1,3 w(a)=rcn*(cos(g)*vik(a)+sin(g)*u(a)) 11 r2(a)=r1(a)+w(a) write(9,900)6,r2,i,nat+ianew+1,nat+ianew+2,0,0,0,0 900 format(i6,3f12.6,7i6) call norm(w) spuw=sp(u,w) do 1 a=1,3 1 up(a)=u(a)-w(a)*spuw call norm(up) do 2 a=1,3 2 r3(a)=r2(a)+rco*(w(a)*cos(alp)+up(a)*sin(alp)) ianew=ianew+1 write(9,900)8,r3,nat+ianew-1,0,0,0,0,0,0 do 3 a=1,3 3 r4(a)=r2(a)+rcc*(w(a)*cos(bp)-up(a)*sin(bp)) ianew=ianew+1 write(9,900)6,r4,nat+ianew-2,nat+ianew+1,nat+ianew+2, 1nat+ianew+3,0,0,0 call methyl(r3,r2,r4,xh1,xh2,xh3) ianew=ianew+1 write(9,900)1,xh1,nat+ianew-1,0,0,0,0,0,0 ianew=ianew+1 write(9,900)1,xh2,nat+ianew-2,0,0,0,0,0,0 ianew=ianew+1 write(9,900)1,xh3,nat+ianew-3,0,0,0,0,0,0 return end subroutine setvk(nat,ix,bt,r,i,j,k,vik,u,g,r1,ro) implicit none integer*4 nat,ix,bt(nat,7),a,i,j,k real*8 r(3*nat),vik(3),vji(3),v(3),u(3),al,b,g,r1(3),ro(3),pi,sp pi=4.0d0*datan(1.0d0) if(ix.eq.1)then k=bt(i,2) else k=bt(i,1) endif do 10 a=1,3 r1(a)=r(a+3*(i-1)) ro(a)=r(a+3*(j-1)) vik(a)=r(a+3*(i-1))-r(a+3*(k-1)) 10 vji(a)=r(a+3*(j-1))-r(a+3*(i-1)) call vp(v,vik,vji) call vp(u,vik,v) call norm(u) call norm(vik) call norm(vji) al=dacos(-sp(vik,vji)) b=(2.0d0*pi-al)/2.0d0 g=pi-b return end subroutine methyl(ro,r1,r2,xh1,xh2,xh3) c add methyl hydrogens to the 1-2 bond implicit none integer*4 ix real*8 v12(3),rch,xh1(3),xh2(3),xh3(3),t(3),cv(3), 1vco(3),ro(3),r1(3),r2(3),w(3),pi,c pi=4.0d0*atan(1.0d0) rch=1.085d0 do 1 ix=1,3 1 v12(ix)=r2(ix)-r1(ix) call norm(v12) c get CO direction do 6 ix=1,3 6 vco(ix)=ro(ix)-r1(ix) call vp(w,v12,vco) call norm(w) call vp(t,v12,w) call norm(vco) call norm(t) c=pi-acos(-1.0d0/3.0d0) do 7 ix=1,3 7 cv(ix)=cos(c)*v12(ix)+sin(c)*t(ix) do 8 ix=1,3 8 xh1(ix)=r2(ix)+rch*cv(ix) do 11 ix=1,3 xh2(ix)=r2(ix)+(v12(ix)/3.0d0-t(ix)*dsqrt(8.0d0)/6.0d0 2-dsqrt(2.0d0/3.0d0)*w(ix))*rch 11 xh3(ix)=r2(ix)+(v12(ix)/3.0d0-t(ix)*dsqrt(8.0d0)/6.0d0 2+dsqrt(2.0d0/3.0d0)*w(ix))*rch return end subroutine rc(nat,ty,r,f) integer*4 nat,ty(nat),ia,ix real*8 r(3*nat) character*(*) f open(2,file=f,status='old') read(2,*) read(2,*) do 1 ia=1,nat 1 read(2,*)ty(ia),(r(ix+3*(ia-1)),ix=1,3) close(2) return end subroutine vp(vy,vz,vx) real*8 vy(3),vz(3),vx(3) vy(1)=vz(2)*vx(3)-vz(3)*vx(2) vy(2)=vz(3)*vx(1)-vz(1)*vx(3) vy(3)=vz(1)*vx(2)-vz(2)*vx(1) return end subroutine rdim(nat,n,nc,m) integer*4 nat,n,nc,m open(2,file='CRYSTAL.X',status='old') read(2,*) read(2,*)nat open(2,file='CELL.X',status='old') read(2,*) read(2,*)nc close(2) close(2) n=3*nat m=3*nc write(6,*)nat,' atoms in CRYSTAL.X',nc,' atoms in CELL.X' return end function sp(a,b) real*8 a(*),b(*),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine norm(v) real*8 v(*),n,sp n=dsqrt(sp(v,v)) v(1)=v(1)/n v(2)=v(2)/n v(3)=v(3)/n return end subroutine ropt(ic,lcorner,ledge) integer*4 ic logical lex,lcorner,ledge character*4 key lcorner=.true. ledge=.true. ic=0 inquire(file='CVIB.OPT',exist=lex) if(lex)then open(9,file='CVIB.OPT') 1 read(9,900,end=99,err=99)key 900 format(a4) if(key.eq.'CTYP')read(9,*)ic c ic=0... methyls adde as caps c ic=1 .. hydrogens added c ic=2 .. hydrogens just replacing troubled group c ic=3 .. rebuild the entire amide groups if(key.eq.'CORN')read(9,*)lcorner c lcorner ... interactions with the corner cells in 3x3x3 supercell if(key.eq.'EDGE')read(9,*)ledge c ledge ... interactions with the edger cells in 3x3x3 supercell goto 1 99 close(9) endif return end subroutine rv(u) real*8 u(3,3) open(9,file='CELL.TXT') read(9,*)u(1,1),u(1,2),u(1,3) read(9,*)u(2,1),u(2,2),u(2,3) read(9,*)u(3,1),u(3,2),u(3,3) end subroutine mkbonds(nat,ty,r,bt,nt) implicit none integer*4 nat,ty(nat),bt(nat,7),nt(nat),i,j real*8 r(3*nat),rb,rt,xi,yi,zi,d real*4 bonding(88) c Du H He Li Be B C N O F c Ne Na Mg Al Si P S Cl Ar K Ca Sc Ti c V Cr Mn Fe Co Ni Cu Zn Ga Ge As 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/ nt=0 bt=0 do 3 i=1,nat if(ty(i)+1.lt.1.or.ty(i)+1.gt.89)call report('type out of range') rb=dble(bonding(ty(i)+1)) xi=r(1+3*(i-1)) yi=r(2+3*(i-1)) zi=r(3+3*(i-1)) do 3 j=i+1,nat if(ty(i)+1.ne.2.or.ty(j)+1.ne.2)then rt=(rb+dble(bonding(ty(j)+1)))**2 d=(xi-r(1+3*(j-1)))**2+(yi-r(2+3*(j-1)))**2+(zi-r(3+3*(j-1)))**2 if(d.lt.rt)then nt(i)=nt(i)+1 if(nt(i).gt.7)call report('too many (>7) bonds') nt(j)=nt(j)+1 if(nt(j).gt.7)call report('too many (>7) bonds') bt(i,nt(i))=j bt(j,nt(j))=i endif endif 3 continue return end subroutine report(s) character*(*) s write(6,*)s stop end