program numforce
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      parameter (nat0=1000)
      character*80 filename,st,standa
      character*1 ok
      character*2 at
      real*8 r0(3,nat0),g0(3,nat0),r(3),g(3,nat0),f(3*nat0,3*nat0)
      dimension iz(nat0)
      logical lex
c
      write(6,*)' FF from gradients'
      write(6,*)' Z-Matrix orientation required !!!'
      write(6,*)
c
      write(6,*)' Geometry ( .x):'
      read(5,5000)filename
5000  format(a80)
      open(15,file=filename)
      read(15,*)
      read(15,*)nat
      if(nat.gt.nat0)then
       write(6,*)'Too many atoms'
       close(15)
       stop
      endif 
      do 1 i=1,nat
1     read(15,*)iz(i),(r0(j,i),j=1,3)
      close(15)
c
      write(6,*)nat,' atoms read'
      write(6,*)' Diff step (A):'
      read(5,*)dd
      write(6,*)'Standard options:'
      standa='#PM3 FORCE NOSYMM iop(2/11=1)'
 2    write(6,5000)standa
      write(6,*)' Change options (y/n)?'
      read(5,'(a)')ok
      if(ok.eq.'y'.or.ok.eq.'Y')then
       write(6,*)'Input options:'
       read(5,5000)standa
       goto 2
      endif
      write(6,*)' G98 input filename:'
      read(5,5000)filename
      open(15,file=filename)
      write(15,*)'%chk=x.chk'
      write(15,*)'%nproc=1'
      write(15,5000)standa
      write(15,505)
505   format(/,'zero',/,/,'0 1')
      do 3 i=1,nat
      at='xx'
      if(iz(i).eq.1)at='h'
      if(iz(i).eq.6)at='c'
      if(iz(i).eq.7)at='n'
      if(iz(i).eq.8)at='o'
3     write(15,5002)at,(r0(j,i),j=1,3)
5002  format(a2,3f15.8)
      write(15,*)
      do 4 ia=1,nat
      do 4 ix=1,3
      write(15,503)
503   format('--link1--')
      write(15,*)'%chk=x.chk'
      write(15,*)'%nproc=1'
      write(15,5000)standa
      write(15,506)ia,ix
506   format(/,'atom',i6,' coord',i3,/,/,'0 1')
      do 31 i=1,nat
      at='xx'
      if(iz(i).eq.1)at='h'
      if(iz(i).eq.6)at='c'
      if(iz(i).eq.7)at='n'
      if(iz(i).eq.8)at='o'
      x=r0(1,i)
      y=r0(2,i)
      z=r0(3,i)
      if(i.eq.ia.and.ix.eq.1)x=x+dd
      if(i.eq.ia.and.ix.eq.2)y=y+dd
      if(i.eq.ia.and.ix.eq.3)z=z+dd
31    write(15,5002)at,x,y,z
4     write(15,*)
      close(15)
      write(6,*)'Input written'
      stop
      END