program gnmrpara
      implicit none
      character*80 s80,fn
      integer*4,allocatable::iz(:)
      real*8,allocatable::r(:,:),fc(:)
      real*8 g(3,3),a(3),b(3),cm(3)
      integer*4 l,i,j,ix,it,ie,ia,iu,ih
      real*8 s,u,giso,acm,akt,t,ak,uN,uB,x,y,z,xi,yi,zi,ri,ch,
     1gN,ppm,sp
      character*4 s4

      write(6,*)
     1' Exctracts paramagnetic NMR parameters from Gaussian output.'
      write(6,*)'Filename:'
      read(5,'(a)')fn
      it=0
      ie=0
      ia=0
      iu=0
      s=0.0d0
      u=1.0d0
88    open(9,file=fn)
      l=0
1     l=l+1
      read(9,900,end=99,err=99)s80
900   format(a80)
      if(s80(14:27).eq.'Multiplicity =')then
       iu=iu+1
       read(s80(28:29),*)u
       s=(u-1.0d0)/2.0d0
       if(ia.gt.0)write(6,900)s80
      endif
      if(s80(27:59).eq.'Isotropic Fermi Contact Couplings')then
       write(6,*)
       write(6,900)s80
       read(9,*)
       do 5 i=1,ia
       read(9,900)s80
5      read(s80(63:76),*)fc(i)
       ie=ie+1
       write(6,*)
      endif
      if(s80(1:10).eq.' g tensor ')then
       read(9,901)((g(i,j),i=1,3),j=1,3)
901    format(3(6x,E15.8))
       write(6,*)
       write(6,900)s80
       write(6,901)((g(i,j),i=1,3),j=1,3)
       it=it+1
      endif
      if(s80(27:44).eq.'Input orientation:')then
       do 2 i=1,4
2      read(9,*)
       if(ia.eq.0)then
3       read(9,900)s80
        if(s80(2:3).ne.'--')then
         ia=ia+1
         goto 3
        endif
        write(6,*)ia,' atoms'
        if(ia.eq.0)stop
        allocate(iz(ia),r(ia,3),fc(ia))
        close(9)
        goto 88
       else
        do 4 i=1,ia
4       read(9,*)iz(i),iz(i),r(i,1),(r(i,j),j=1,3)
        write(6,*)ia,' atoms read'
       endif
      endif
      goto 1
99    close(9)
      write(6,*)l,' lines'
      if(ia.gt.0.and.it.gt.0.and.ie.gt.0)then
       ih=1
       do 7 i=2,ia
7      if(iz(i).gt.iz(ih))ih=i
       write(6,*)'Heaviest atom:',ih
       do 8 ix=1,3
       cm(ix)=0.0d0
       do 81 i=1,ia
81     if(i.ne.ih)cm(ix)=cm(ix)+r(i,ix)
8      cm(ix)=cm(ix)/dble(ia-1)

       do 82 ix=1,3
82     a(ix)=cm(ix)-r(ih,ix)

c      Bohr magneton:
       uB=9.27400968d-24
c      Nuclear magneton:
       uN=5.05078353d-27
c      Boltzmann constant:
       ak=1.38064852d-23
c      temperature:
       t=300.0d0
c      kT:
       akt=t*ak
c      kt/cm-1:
       acm=akt*5.034d22
c      isotropic g-tensor:
       giso=(g(1,1)+g(2,2)+g(3,3))/3.0d0
       write(6,607)uB,uN,s,giso,acm
607    format(' uBohr: ',e12.4,' SI',/,
     1        ' uNucl: ',e12.4,' SI',/,
     1        ' S    : ',f12.1,' ',/,
     1        ' g_iso: ',f12.3,' ',/,
     1        ' kT   : ',f12.1,' cm-1 ',/)

       if(dabs(giso).gt.20)then
        write(6,*)
        write(6,*)'g_iso unreasonably big, 2.00231930436182 taken'
        write(6,*)
        giso=2.00231930436182d0
       endif

       write(6,*)'Paramagnetic NMR parameters'
       write(6,*)'ATOM  IZ       shift_P        (3cos(th)^2 - 1)/ri^3 '
       x=r(ih,1)
       y=r(ih,2)
       z=r(ih,3)
       do 6 i=1,ia
       gN=0.0d0
       if(iz(i).eq.1)then
        gN=5.5748d0
        s4=' 1H '
       endif
       if(iz(i).eq.6)then
        gN=1.04044d0
        s4='13C '
       endif
       if(iz(i).eq.7)then
        gN=0.40362d0
        s4='14N '
       endif
       if(iz(i).eq.9)then
        gN=5.2553d0
        s4='19F '
       endif
       if(gN.ne.0.0d0)then
        ppm=uB/uN/gN*s*(s+1.0d0)/3.0d0*giso*(fc(i)*1.0d-4/acm)*1.0d6
        xi=r(i,1)
        yi=r(i,2)
        zi=r(i,3)
        ri=sqrt((x-xi)**2+(y-yi)**2+(z-zi)**2)
        do 83 ix=1,3
83      b(ix)=r(i,ix)-r(ih,ix)
        ch=sp(a,b)/dsqrt(sp(a,a)*sp(b,b))
        write(6,609)i,iz(i),s4,ppm,(3.0d0*ch**2-1.0d0)/ri**3
609     format(i6,i3,1x,a4,1x,f12.1,f12.4)
       endif
6      continue
      endif
      end
      
      function sp(a,b)
      implicit none
      real*8 a(3),b(3),sp
      sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
      return
      end