program molx implicit none integer*4 ib,i,j,ii,ix,nat,na,i1,i2,j1,jj,n100, 1N7,ig,LEXCL,ic,IEX,iu,NExc,NSL,iom parameter (n100=100) real*8 p1,p1o,dif real*8,allocatable::rtem(:,:),q(:),xyz(:,:),step(:),r7(:), 1a7(:),t7(:),f0(:),t7n(:) integer*4,allocatable:: back(:),it(:),ifipsi(:),nps(:,:), 1nn(:),i7(:),j7(:),STATE(:),si(:) character*100 s character*80 s80 character*4 s4 character*2 s2 logical lscan,lex,lomit write(6,600) 600 format(' .mol -> .x format conversion',/, 1 ' New mol mcm95 format',/, 1 ' Input : FILE.MOL',/, 1 ' Output: FILE.X') open(9,file='FILE.MOL',status='old') read(9,900)s80 900 format(a80) read(9,*)ib nat=ib+1 allocate(rtem(nat+2,3),q(nat),it(nat),back(nat+2),xyz(3,nat), 1i7(ib),j7(ib),r7(ib),a7(ib),t7(ib),t7n(ib)) read(9,*)it(1),q(1) read(9,*)(rtem(1,ix),ix=1,3) do 8 ix=1,3 8 xyz(ix,1)=rtem(1,ix) do 201 ii=1,ib 201 read(9,*)i7(ii),j7(ii),it(ii+1),r7(ii),a7(ii),t7(ii),q(ii+1) close(9) lscan=.false. lomit=.false. allocate(ifipsi(1),nn(1),step(1),nps(1,1),f0(1)) inquire(file='MOLX.OPT',exist=lex) if(lex)then open(9,file='MOLX.OPT') 102 read(9,444,end=99,err=99)s4 444 format(a4) if(s4.eq.'SCAN')then read(9,*)lscan read(9,*)na deallocate(ifipsi,nn,step,nps,f0) allocate(ifipsi(na),nn(na),step(na),nps(na,nat),f0(na)) if(lscan)write(6,*)na,' classes to scan' do 101 i=1,na read(9,*)ifipsi(i) if(lscan)write(6,*)i,ifipsi(i) do 1011 j=1,ifipsi(i) read(9,*)i1,i2 jj=0 do 1012 ix=1,ib 1012 if(i7(ix).eq.i1.and.j7(ix).eq.i2)jj=ix if(jj.eq.0)then write(6,*)i1,i2 call report('requested atom pair not found in bond list') endif if(lscan)write(6,*)i1,i2,' bond ',jj 1011 nps(i,j)=jj 101 read(9,*)nn(i),f0(i),step(i) endif if(s4.eq.'OMIT')read(9,*)lomit goto 102 99 close(9) endif back(i)=0 back(1)=nat+1 back(nat+1)=nat+2 rtem(nat+1,1)=rtem(1,1) rtem(nat+1,2)=rtem(1,2) rtem(nat+1,3)=rtem(1,3)-1.0d0 rtem(nat+2,1)=rtem(1,1)+1.0d0 rtem(nat+2,2)=rtem(1,2) rtem(nat+2,3)=rtem(1,3)-1.0d0 ig=0 iom=0 if(lscan)then write(6,601) 601 format(' Scan requested') t7n=t7 C Ground state-always/desired angles: do 103 iu=1,na p1=f0(iu) call c360(p1) c loop over all instances of this angle do 103 j1=1,ifipsi(iu) c p1o - original angle, dif - difference: p1o=t7(nps(iu,j1)) dif=p1-p1o t7n(nps(iu,j1))=p1 c look at dependent angles do 103 ii=1,ib if(i7(ii).eq.i7(nps(iu,j1)).and.ii.ne.nps(iu,j1))then t7n(ii)=t7(ii)+dif call c360(t7n(ii)) endif 103 continue call ggx(nat,ib,rtem,xyz,i7,j7,r7,a7,t7n,back,'0.X',s80,q,it, 1 iom,lomit,ig) c c make notation consistent with previous code: LEXCL=0 do 104 iu=1,na 104 LEXCL=LEXCL+nn(iu)-1 N7=1 NSL=na allocate(STATE(na),si(LEXCL)) write(6,607)na,LEXCL,N7 607 format(' na:',i3,' LEXCL:',i3,' N7:',i3) c states Nexc excited do 203 Nexc=1,LEXCL c distribute Nexc excitations upon centers N7..NSL:,angles 1...na c initial indices - all to the first center: do 41 iex=1,Nexc 41 si(iex)=N7 205 STATE=0 do 204 iex=1,Nexc 204 state(si(iex))=state(si(iex))+1 c check that the state is allowed: do 211 I=N7,NSL 211 if(state(I).gt.nn(I)-1)goto 212 c figure out name ic=1 do 110 iu=1,na if(STATE(iu).gt.99)call report('format limit to two characters') write(s2,299)STATE(iu) 299 format(i2) if(STATE(iu).lt.10)then s(ic:ic)=s2(2:2) ic=ic+1 else s(ic:ic+1)=s2(1:2) ic=ic+2 endif if(iu.lt.na)then s(ic:ic)='_' ic=ic+1 endif if(ic+5.gt.n100)call report('name too long') 110 continue s(ic:ic+1)='.X' ic=ic+1 c c set the desired value for each angle and record geometry: do 111 iu=1,na p1=f0(iu)+step(iu)*STATE(iu) call c360(p1) do 111 j1=1,ifipsi(iu) p1o=t7(nps(iu,j1)) dif=p1-p1o t7n(nps(iu,j1))=p1 do 111 ii=1,ib if(i7(ii).eq.i7(nps(iu,j1)).and.ii.ne.nps(iu,j1))then t7n(ii)=t7(ii)+dif call c360(t7n(ii)) endif 111 continue call ggx(nat,ib,rtem,xyz,i7,j7,r7,a7,t7n,back,s(1:ic),s80,q,it, 1 iom,lomit,ig) c find index to be changed 212 do 208 ic=Nexc,1,-1 208 if(si(ic).lt.NSL)goto 209 goto 203 209 do 210 i=ic+1,Nexc 210 si(i)=si(ic)+1 si(ic)=si(ic)+1 c goto 205 c 203 continue write(6,*)ig,' geometries' if(iom.gt.0)write(6,*)iom,' ommited' C else call ggx(nat,ib,rtem,xyz,i7,j7,r7,a7,t7,back,'FILE.X',s80,q,it, 1 iom,lomit,ig) endif end function sp(v1,v2) real*8 v1(*),v2(*),sp sp=v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3) return end subroutine vp(v1,v2,v3) real*8 v1(*),v2(*),v3(*) v3(1)=v1(2)*v2(3)-v1(3)*v2(2) v3(2)=v1(3)*v2(1)-v1(1)*v2(3) v3(3)=v1(1)*v2(2)-v1(2)*v2(1) return end subroutine norm(v) real*8 v(*),a,sp a=1.0d0/dsqrt(sp(v,v)) v(1)=v(1)*a v(2)=v(2)*a v(3)=v(3)*a return end subroutine c360(a) real*8 a if(a.gt. 180.0d0)a=a-360.0d0 if(a.lt.-180.0d0)a=a+360.0d0 return end subroutine report(s) character*(*) s write(6,*)s stop end subroutine ggx(nat,ib,rtem,xyz,i7,j7,r7,a7,t7,back,fn,s80,q,it, 1iom,lomit,ig) implicit none integer*4 ia,ib,l,ii,nat,i7(*),j7(*),back(*),it(*),ix,i,j,ja,iom, 1ig real*8 rtem(nat+2,3),xyz(3,nat),r7(*),a7(*),t7(*),d, 1r3,a,t,r2(3),r1(3),a1,a2,a3,r1p(3),sp,q(*),pi character*(*) fn character*80 s80 logical lomit pi=4.0d0*atan(1.0d0) do 2 ii=1,ib l=ii+1 i=i7(ii) j=j7(ii) r3=r7(ii) a=a7(ii) t=t7(ii) back(j)=i if(back(i).eq.0.or.back(back(i)).eq.0)then write(6,*)'Ill-defined coordinates, bond ',ii stop endif do 3 ia=1,3 r2(ia)=rtem(i ,ia)-rtem(back(i) ,ia) 3 r1(ia)=rtem(back(i),ia)-rtem(back(back(i)),ia) a1=-r3*sin((180.0d0-a)*pi/180.0d0)*cos(t*pi/180.0d0) a2=-r3*cos(a*pi/180.0d0) if(r3**2-a1**2-a2**2.gt.0.0d0)then a3=dsqrt(r3**2-a1**2-a2**2) else a3=0.0d0 endif if(t.lt.0.0d0)a3=-a3 do 4 ia=1,3 4 r1p(ia)=r1(ia)-sp(r1,r2)*r2(ia)/sp(r2,r2) call norm(r1p) call norm(r2) call vp(r1p,r2,r1) do 5 ia=1,3 5 rtem(j,ia)=rtem(i,ia)+a1*r1p(ia)+a2*r2(ia)+a3*r1(ia) do 6 ia=1,3 6 xyz(ia,l)=rtem(j,ia) 2 continue ig=ig+1 c check that atoms are not too close: do 71 ia=1,nat do 71 ja=ia+1,nat d=(xyz(1,ia)-xyz(1,ja))**2+(xyz(2,ia)-xyz(2,ja))**2+ 1(xyz(3,ia)-xyz(3,ja))**2 if(d.lt.0.25d0)then if(lomit)then write(6,609)ia,ja,dsqrt(d) 609 format(2i6,f12.6,' A, atoms too close, geometry ommited') iom=iom+1 return else write(6,619)ig,ia,ja,dsqrt(d) 619 format(' Warning: geom',i6,',',2i6,f12.6,' A, atoms too close') endif endif 71 continue open(10,file=fn) write(10,900)s80 900 format(a80) write(10,*)nat do 7 ia=1,nat 7 write(10,1000)it(ia),(xyz(ix,ia),ix=1,3),q(ia) 1000 format(i6,3f12.6,' 0 0 0 0 0 0 0 ',f10.2) close(10) return end