c c search for a molecule, using a list of seed atoms c subroutine psmol(nseed,iseed,ipun) implicit real*8(a-h,o-z) include 'crystal.sizes' common /xbasat/exx(lim042),c1(lim042),c2(lim042),c3(lim042), * xa(3,lim016),aznuc(lim016), * che(lim015),xl(3,lim015),exad(lim015+1), * nat(lim016),nsh(lim016),nshpri(lim016+1),naopri(lim016+1), * ipseud(lim016), * laa(lim015+1),lan(lim015),lat(lim015),latao(lim015), * nds(lim015+1),jan(lim015),latoat(lim015) common /gvect/xg(3,lim006),paret(3,3),w1r(3,3),gmodu(lim007), * gmodus(lim007),lg(3,lim006),nm(lim007),nn1(lim006) common /parinf/ *par(50),lprint(200),iunit(50),itol(50),inf(200),itit(20),iin,iout logical omore, omol(lim016) dimension t(3), de(3,3), xv(3,3) common/molecu/xt(3,lim016),nmol,omol logical opointc common/pointc/opointc common/charges/ucha(lim016) dimension rcov(92) data rcov/ 0.65d0,1.40d0,1.57d0,1.12d0,0.89d0,0.77d0,0.74d0,0.74d0 *,0.72d0,1.60d0, 1.91d0,1.50d0,1.43d0,1.17d0,1.10d0,1.04d *0,0.99d0,1.88d0,2.20d0,1.97d0, 1.60d0,1.40d0,1.35d0,1.40 *d0,1.40d0,1.40d0,1.35d0,1.35d0,1.35d0,1.35d0, 1.30d0,1.2 *5d0,1.15d0,1.15d0,1.14d0,2.00d0,2.20d0,2.00d0,1.85d0,1.55d0, * 1.45d0,1.45d0,1.35d0,1.30d0,1.35d0,1.40d0,1.60d0,1.55d0,1.55d0,1. *45d0,1.45d0,1.40d0,1.40d0,2.53d0,2.60d0,2.00d0,1.75d0,1.55d0,1.55d *0,1.55d0,1.55d0,1.55d0,1.55d0,1.55d0,1.55d0,1.55d0,1.55d0,1.55d0,1 *.55d0,1.55d0,1.55d0,1.55d0,1.45d0,1.35d0,1.35d0,1.30d0,1.35d0,1.35 *d0,1.35d0,1.50d0,1.90d0,1.80d0,1.60d0,1.55d0,1.55d0,1.55d0,2.80d0, c1.60d0,1.95d0,1.55d0,1.55d0,1.55d0/ character*2 symbat(0:92) data symbat/ 'xx','h ','he','li','be','b ','c ','n ', 1'o ','f ','ne','na','mg','al','si','p ','s ','cl','ar', 2'k ','ca','sc','ti','v ','cr','mn','fe','co','ni','cu', 3'zn','ga','ge','as','se','br','kr','rb','sr','y ','zr', 4'nb','mo','tc','ru','rh','pd','ag','cd','in','sn','sb', 5'te','i ','xe','cs','ba','la','ce','pr','nd','pm','sm', 6'eu','gd','tb','dy','ho','er','tm','yb','lu','hf','ta', 7'w ','re','os','ir','pt','au','hg','tl','pb','bi','po', 8'at','rn','fr','ra','ac','th','pa','u '/ character*80 title common/ztitl/title dimension iseed(*) tol=0.5d0 do j=1,inf(24) do k=1,3 xt(k,j)=xa(k,j) enddo omol(j)=.false. enddo write(6,*)'seeds to be used :' do i = 1,nseed j = iseed(i) omol(j) = .true. write(6,100)j,nat(j),xt(1,j),xt(2,j),xt(3,j) enddo write(6,*)'other atoms in molecule :' nmol=nseed ipass=0 10 continue omore=.false. ipass=ipass+1 if(ipass.gt.20)then write(6,*)'too many passes' call exit(1) endif do j=1,inf(24) if(.not.omol(j))then do i = 1,inf(24) if(omol(i).and.i.ne.j)then rr=aadist(xt(1,i),xt(1,j)) cov=(rcov(nat(i))+rcov(nat(j)))/0.57+tol r2=cov*cov rtest=rr+cov ifin=0 do 45 istar=1,inf(5) rstar = gmodu(istar) if(rstar.gt.rtest)goto 20 istart=ifin+1 ifin=istart+nm(istar)-1 do 40 ine=istart,ifin if(ine.gt.999.or.ine.lt.0)then write(6,*)'crap vector in psmol' call exit(1) endif t(1)=xt(1,j)+xg(1,ine) t(2)=xt(2,j)+xg(2,ine) t(3)=xt(3,j)+xg(3,ine) if(dist2(t,xt(1,i)).lt.r2)then xt(1,j)=t(1) xt(2,j)=t(2) xt(3,j)=t(3) write(6,100)j,nat(j),xt(1,j),xt(2,j),xt(3,j) 100 format(1x,i4,i3,2x,3f12.8) omol(j)=.true. omore=.true. nmol=nmol+1 goto 50 endif 40 continue 45 continue 20 continue endif enddo endif 50 continue enddo if(omore)goto 10 if(ipun.ne.0)then write(7,99)title(1:40),nmol 99 format(1x,'block=fragment records=0 index=1',/, & 1x,'block=title records= 1 index=1',/, & 1x,'home molecule ',a40,/, & 1x,'block=coordinates records=',i4,' index=1') do j=1,inf(24) if(omol(j))write(7,101)symbat(nat(j)),xt(1,j),xt(2,j),xt(3,j) enddo write(7,199)nmol 199 format(1x,'block= potential_derived_charges records=',i4) tcha = 0.0 do j=1,inf(24) if(omol(j)) then write(7,102)ucha(j) 102 format(f10.5) tcha = tcha + ucha(j) endif enddo write(7,198) tcha 198 format(1x,'block= charge records= 1',/,f10.5) c c primitive unit cell vectors c ndim = inf(10) write(7,130)ndim write(7,140)(paret(1,j),j=1,3) if(ndim.gt.1)write(7,140)(paret(2,j),j=1,3) if(ndim.gt.2)write(7,140)(paret(3,j),j=1,3) endif if(ipun.eq.2 .and. opointc )then write(7,299)nmol 299 format(1x,'block= point_charges records=',i4) do j=1,inf(24) if(omol(j))write(7,103)ucha(j),xt(1,j),xt(2,j),xt(3,j) 103 format(f10.5,3(1x,f20.10)) enddo endif 101 format(1x,a2,2x,3f22.14) 130 format(1x,'block=cell_vectors records=',i4,' index=1') 140 format(1x,3f20.12) return end function aadist(c1,c2) implicit real*8(a-h,o-z) dimension c1(3),c2(3) aadist=dsqrt(dist2(c1,c2)) return end function dist2(c1,c2) implicit real*8(a-h,o-z) dimension c1(3),c2(3) dist2 = (c1(1)-c2(1))**2 + (c1(2)-c2(2))**2 + & (c1(3)-c2(3))**2 return end