C     program POPopt_mi_2023
C
C     Polarized one-electron potential (POP) fitting    (Matrix Inversion)
C
C     June 1992 - Feb 2021  by S. Nakagawa
C
C     2023  2/22 Cleaned up the source code.
C           2/23 A model in which the positive charge of the induced dipole is 
C                placed on the atom is added (ictat=1).
c     ---------------------------------------------------------
c     irest = 3 intra & intermolecular atom polarizability fit simple scaling 
c           = 4 thole exponential model  van duijnen type a**3/8piexp(-au)
c           = 5 thole exponential model  ponder type      3a/4piexp(-au**3)
C           = 6 Masia's model  Gaussian density screening  Atom Type a=0.926
C           =43 Thole exponential model  van duijnen type a**3/8piexp(-au)
C               with simple scaling
C           =53 Thole exponential model  ponder type      3a/4piexp(-au**3)
C               with simple scaling
C           =63 Masia's model  Gaussian density screening  Atom Type a=0.926
C               with simple scaling
C           =73 Thole Linear model     fe=4v**3-3v**4   ft=v**4   (rij<s)  
C               with simple scaling
c     ---------------------------------------------------------
c     ivary = 0 the variable of the center is varying
c           = 1 the variable of the center is frozen
c           = 2 the variable of the center is equal to one of 
c               the other centre    
c     ---------------------------------------------------------
c
c     files 
c
c     unit 1  file name is rest_1  
c             output pop data by fitted polarizabilitiy,
c             induction energies and induced dipole moments 
c
c     unit 2  file name is rest_2  
c             input  esp & pop data from ab initio calculations(Gaussian 09)
c
c     unit 3  file name is rest_3 : input data    read (3,*)
c   
c    1. ishl   : number of vdw surface shell
c       ifit   : type of fitting / e: field potential
c                0 alpha=a
c                1 alpha=a+be
c       irest  : 2 molecular polarization model 
c                  (no intramolecular induced dipoles)
c                3 inter- and intramolecular polarization model
c                  (scale 1-2,1-3,1-4)
c                4 thole's exponential type intramolecular a**3/8piexp(-au)
c                  a=1.9088(default)(van Duijnen) a=2.089(Thole)
c                5 thole's exponential type intramolecular 3a/4piexp(-au**3)
c                  a=0.39(default)(Ponder) a=0.572(Thole)
c                  (irest= 3 or 4 or 5 require link list (unit 4))
C                6 Masia's model  Gaussian density screening  Atom Type a=0.926
C               43 Thole's exponential type intramolecular a**3/8piexp(-au)
C                  with simple scaling
C               53 Thole's exponential type intramolecular 3a/4piexp(-au**3)
C                  with simple scaling
C               63 Masia's model  Gaussian density screening  Atom Type a=0.926
C                  with simple scaling
C        ictat : 0 A model in which the center of the induced dipole is placed on the atom 
C                1 A model in which the positive charge of the induced dipole is placed on the atom 
c
c    2. fpchrg : electric field point charge for the shell  (A test charge 0.1e is recommended)
c   
c    3. iptyp  : type of polarizability
c                0 transvers + anisotropy (bond polz) 
c                1 transvers  
c                2 anisotropy 
c                3 atom + anisotropy (atomic polz)
c                4 atom
c   
c    4. ibpair : # of bond pair for bond polarizability
c       iatid  : # of atom for atom polarizability 
c   
c    4-1. ibpair & iatid data
c   
c    5. ifld   : number of field point charge on shell
c                999 positions of fpc equal esp)
c                888 symmetry is used for positions of fpc)
c                777 real symmetry is used for positions of fpc/
c                    faster than 888!)
c   
c    5-1. fpcxyz or 999 or 888 or 777  data
c   
c      =999 no additional data
c      =888 or =777
c           number of symmetry data (free format)
c           symmetry data (5x,2i5)
c      =number fpcx   : xyz coordinates of fpc (a.u.)
c              fpcy     
c              fpcz
c   
c     unit 4  input  atom link data
c
c     unit 7  output data
c
c     unit 8  output final polz parameters
c
      implicit real*8(a-h,o-z)
      parameter (i345=345,i346=i345+1,i347=i345+2)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/espdat/esfact,iestat,irest,ictat
      common/comp12/cxnp(i340),cynp(i340),cznp(i340),
     & fpcx(i340),fpcy(i340),fpcz(i340),
     & qupot(i340,i340,i20)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
C   
      common/cominv/v(i100*i100*i100)
C   
      common/comp34/fcal(i345*i346*i20)
      dimension vv(i340*i340*i20)
C   
      dimension w(i340*i340*i100)
      dimension var(i200*4)
C
      call cpu_time(tt1)
c
c       rest_1 : write restart data to   rest_1
c       rest_2 : read  restart data from rest_2
c
      open(1,file='rest_1',access='sequential',form='formatted')
      open(2,file='rest_2',access='sequential',form='formatted')
      open(3,file='rest_3',access='sequential',form='formatted')
c
      call esread
      close(unit=2)
      close(unit=3)
      call setup(var)
      i0= 1+nvar*(nvar+1)/2
      i1=i0+nvar*(nvar+1)/2
      i2=i1+nvar
      i3=i2+nvar
      i4=i3+nvar
      i5=i4+nvar
      j1=i5+nvar
      j2=j1+iestat
      ibig=i340*i340*i20
      if(nvar.gt.70)
     &   stop 'memory over  nvar'
      if(iestat.gt.ibig)
     &stop 'memory over  iestat'
      call setbck(var)
      call sumsq (var,fcal)
      call optssq(var,fcal,w,
     &   v,v(i0),v(i1),v(i2),v(i3),v(i4),v(i5),v(j1),vv)
      call potout(var,fcal)
c
      if(irest.eq.3.or.irest.eq.4.or.irest.eq.5.or.irest.eq.6
     &  .or.irest.eq.43.or.irest.eq.53
     &  .or.irest.eq.63.or.irest.eq.73) then
        call matout
      endif
C
      call cpu_time(tt2)
      write(7,*) '  '
      write(7,*) 'CPU time:',tt2-tt1,' second'
C
      stop
      end
c**********************************************************************
      subroutine esread
      implicit real*8(a-h,o-z)
      parameter (i345=345,i346=i345+1,i347=i345+2)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/consts/natm,ndet,nshl,noca,nocb,nocc,mdet,ngto,natt,nover
      common/espdat/esfact,iestat,irest,ictat
      common/comp12/cxnp(i340),cynp(i340),cznp(i340),
     & fpcx(i340),fpcy(i340),fpcz(i340),
     & qupot(i340,i340,i20)
      common/cvlist/cx(3200),cy(3200),cz(3200),cv(3200),npoint
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt(i200)
      common/sekibn/diag(i345),fchrg(3)
c
      common/pdpol1/ibon(i200),jbon(i200),ieqbon(i200),ibpair,iptyp,
     & nvarp,ivarp(i200),iestah,iatid,
     & iatm(i200),
     & nvara,nvart,nvaran,ishl,nvarh,
     & ifit,ifld,jfld,isym2(i340)
      common/pdpol2/
     & fpchrg(i20),pkpt(i200*2,4),
     & eqatmi(i200),varai(i200),
     & alptrn(i200),alpani(i200),alpatm(i200)
c
      dimension xyzkpt(i200,4),isym1(i340)
      equivalence (xyzkpt(1,1),xkpt(1))
c
      iatid=0
      ibpair=0
      iptyp=0
      read (3,*) ishl,ifit,irest,ictat
      read (3,*) (fpchrg(i),i=1,ishl)
      write(7,*) '  '
      write(7,'(4(a7,i3))') 
     &  '  ishl=',ishl,'  ifit=',ifit,' irest=',irest,
     &  ' ictat=',ictat
      write(7,'(a15,10f8.3)') '  field charge:',(fpchrg(i),i=1,ishl)
      read (3,*) iptyp
      write(7,*) '  '
      write(7,'(a15,i3,f8.3)') '  iptyp',
     & iptyp
      read (3,*) ibpair,iatid
      write(7,*) ' polarizability fitting'
      if(ibpair.gt.0)then
        read (3,*) (ibon(i),jbon(i),ieqbon(i),i=1,ibpair)
      endif
      if(iatid.gt.0)then
        read (3,*) (iatm(i),eqatmi(i),i=1,iatid)
      endif
      read (3,*) ifld
      if(ifld.le.700)then
        do i=1,ifld
          read (3,*) fpcx(i),fpcy(i),fpcz(i)
        enddo
      else if(ifld.eq.888.or.ifld.eq.777)then
        read(3,*) ksym
        do i=1,ksym
          read(3,'(5x,2i5)') isym1(i),isym2(i) 
        enddo
      endif
      write(7,*) '  '
      write(7,'(a7,i5)') '  ifld=',ifld
      if(ifld.le.700)then
      do i=1,ifld
        write (7,'(a6,3f10.5)') 'fpcxyz',fpcx(i),fpcy(i),fpcz(i)
      enddo
      else if(ifld.eq.888)then
        write(7,*) 'symmetrized fpcxyz'
      else if(ifld.eq.777)then
        write(7,*) 'real symmetrized fpcxyz '
      else
        write(7,*) 'fpcxyz equal cxnpyz'
      endif
      write(7,*) '  '
      write(7,*) ' bond/atom polarizability fit '
      if(iptyp.eq.1)then
        write(7,*) '  type  transvers'
      else if(iptyp.eq.2) then
        write(7,*) '  type  anisotropy'
      else if(iptyp.eq.3) then
        write(7,*) '  type  atom and anisotropy'
      else if(iptyp.eq.4) then
        write(7,*) '  type  atom'
      else 
        write(7,*) '  type  transvers and anisotropy'
      endif
      if(ibpair.ne.0)then
        write(7,*) '  '
        write(7,*) '  no. bond  i -- j  (similarity)'
        do i=1,ibpair,2
          ii=i+1
          write(7,'(2(i5,2x,i5,a3,i5,5x,a1,i2,a1,8x))') 
     &     i,ibon(i),' --',jbon(i),'(',ieqbon(i),')',
     &     ii,ibon(ii),' --',jbon(ii),'(',ieqbon(ii),')' 
        enddo
      endif
      if(iatid.ne.0)then
        write(7,*) '  '
        write(7,*) '  no. atom  (similarity) ' 
        do i=1,iatid
          write(7,'(i5,2x,i5,5x,a1,f7.3,a1)') 
     &     i,iatm(i),'(',eqatmi(i),')'
        enddo
      endif
      write(7,*) '  '
      jfld=ifld
c
      do l=1,ishl
c   
c       read first data 
c   
        read (2,'(8x,i5,e24.12)') iuniq,fchrg(2)
        read (2,'(i2,f16.10,i2,f16.10,i2,f16.10,i2,f16.10,i2)') 
     &   (iatom(i),xkpt(i),ivary(i,1),ykpt(i),
     &   ivary(i,2),zkpt(i),ivary(i,3),qkpt(i),
     &   ivary(i,4),i=1,iuniq)
        read (2,'(8x,i5)') iestah
        do i=1,iestah
          read (2,'(4f15.10)') 
     &     cxnp(i),cynp(i),cznp(i),qupot(i,1,l)
        enddo
        if(ifld.eq.999) then
          ifld=iestah
          do i=1,iestah
            fpcx(i)=cxnp(i)
            fpcy(i)=cynp(i)
            fpcz(i)=cznp(i)
          enddo
        else if(ifld.eq.888)then
          ifld=iestah
          ks=0
          write(7,*) 'symmetry fpcx'
          do i=1,ksym
            do j=1,isym2(i)
              ks=ks+1
              fpcx(ks)=cxnp(isym1(i))
              fpcy(ks)=cynp(isym1(i))
              fpcz(ks)=cznp(isym1(i))
            enddo
          enddo
        else if(ifld.eq.777)then
          ifld=ksym
          write(7,*) 'real symmetry fpcx'
          do i=1,ksym
              fpcx(i)=cxnp(isym1(i))
              fpcy(i)=cynp(isym1(i))
              fpcz(i)=cznp(isym1(i))
          enddo
        endif
c   
c       read data 2-ifld
c   
        if(jfld.eq.888)then
          ks=1
          do j=2,isym2(1)
            ks=ks+1
            do i=1,iestah
              qupot(i,ks,l)=qupot(i,1,l)
            enddo
          enddo
          do k=2,ksym
            ks=ks+1
            do i=1,iuniq+2
              read (2,'(a4)') dummy
            enddo 
            do i=1,iestah
              read (2,'(45x,f15.10)') qupot(i,ks,l)
            enddo
            kso=ks
            do j=2,isym2(k)
              ks=ks+1
              do i=1,iestah
                qupot(i,ks,l)=qupot(i,kso,l)
              enddo
            enddo
          enddo
        else
          do i=2,ifld
            do j=1,iuniq+2
              read (2,'(a4)') dummy
            enddo 
            do j=1,iestah
              read (2,'(45x,f15.10)') qupot(j,i,l)
            enddo
          enddo
        endif
      enddo
c
c     irest = 3 or 4 or 5  or 6 or 43 or 53 or 63 or 73
c             read scaling param or conectivity 1-2 1-3 1-4 link list
c
      if(irest.eq.3.or.irest.eq.4.or.irest.eq.5.or.irest.eq.6
     &  .or.irest.eq.43.or.irest.eq.53
     &  .or.irest.eq.63.or.irest.eq.73) then
        call lnkplz
      endif
      call prepol
      iestat=iestah*ifld*ishl
c   
  999 do 300 k=1,4
      neq(k)=0
      do 301 i=1,iuniq
  301 if(ivary(i,k).ne.1.and.abs(ivary(i,k)).lt.8) ivary(i,k)=0
      do 310 i=2,iuniq
      if(ivary(i,k).ne.0                       ) go to 310
c     if(dabs(xyzkpt(i,k)            ).lt.1.d-5) go to 310
      do 311 j=1,i-1
      if(ivary(j,k).ne.0                       ) go to 311
      if(dabs(xyzkpt(j,k)            ).lt.1.d-5) go to 311
      if(dabs(xyzkpt(i,k)-xyzkpt(j,k)).gt.1.d-5) go to 312
      neq(k)=neq(k)+1
      iateq(neq(k),k)=i
      jateq(neq(k),k)=j
      ivary(i,k)=2
      ivary(j,k)=0
                                                 go to 311
  312 if(dabs(xyzkpt(i,k)+xyzkpt(j,k)).gt.1.d-5) go to 311
      neq(k)=neq(k)+1
      iateq(neq(k),k)=i
      jateq(neq(k),k)=j
      ivary(i,k)=-2
      ivary(j,k)= 0
  311 continue
  310 continue
  300 continue
c
c    find fixed charge
c   
      iall=0
      ifix=0
      do i=1,iuniq
       if(ivary(i,4).eq.0) then
        ifix=i
        iall=iall+1
       endif
      enddo
      if(iall.eq.0)then
       write(7,*) '**all the charges are fixed**'
      else if(iall.eq.1) then
       ivary(ifix,4)=0
      else
       ivary(ifix,4)=3
       do i=1,neq(4)
        if(jateq(i,4).eq.ifix) then
         if(ivary(iateq(i,4),4).eq.2) ivary(iateq(i,4),4)=3
         if(ivary(iateq(i,4),4).eq.-2) ivary(iateq(i,4),4)=-3
        endif
       enddo
      endif
  220 if(fchrg(2).lt.0.d0) ichrg=fchrg(2)-0.01d0
      if(fchrg(2).ge.0.d0) ichrg=fchrg(2)+0.01d0
                  fchrg(2)=ichrg
c   
      charge=0.d0
      ntst=0
      do i=1,iuniq
       if(ivary(i,4).ne.3.and.ivary(i,4).ne.-3) then
        charge=charge+qkpt(i)
       else
        if(ivary(i,4).eq.3)ntst=ntst+1
        if(ivary(i,4).eq.-3)ntst=ntst-1
       endif
      enddo
      if(ntst.ne.0)then
       qfix=(fchrg(2)-charge)/ntst
       do i=1,iuniq
        if(ivary(i,4).eq.3) qkpt(i)=qfix
        if(ivary(i,4).eq.-3) qkpt(i)=-1.0d0*qfix
       enddo
      endif
      write(7,1)
      write(7,'(9x,2i7,2x,f15.7,i2,f15.7,i2,f15.7,i2,f15.7,i2))')
     &(i,iatom(i),xkpt(i),ivary(i,1),ykpt(i),ivary(i,2),
     &            zkpt(i),ivary(i,3),qkpt(i),ivary(i,4),i=1,iuniq)
      write(7,'(13x,''total'',51x,f22.7)') fchrg(2)
      write(7,'(12x,42(''* ''))')
    1 format(/12x,42(2h* )/13x,'atom   nuclear',18x,
     &     'coordinates (au)',21x,'initial'/20x,'charge',
     &  10x,'x',14x,'y',14x,'z',15x,'charge'/12x,42(2h* ))
      return
      end
c**********************************************************************
      subroutine setup(var)
      implicit real*8(a-h,o-z)
      parameter (i345=345,i346=i345+1,i347=i345+2)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/consts/natm,ndet,nshl,noca,nocb,nocc,mdet,ngto,natt,nover
      common/cvlist/cx(3200),cy(3200),cz(3200),cv(3200),npoint
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/sekibn/diag(i345),fchrg(3)
c
      common/pdpol1/ibon(i200),jbon(i200),ieqbon(i200),ibpair,iptyp,
     & nvarp,ivarp(i200),iestah,iatid,
     & iatm(i200),
     & nvara,nvart,nvaran,ishl,nvarh,
     & ifit,ifld,jfld,isym2(i340)
      common/pdpol2/
     & fpchrg(i20),pkpt(i200*2,4),
     & eqatmi(i200),varai(i200),
     & alptrn(i200),alpani(i200),alpatm(i200)
c
      dimension xyzkpt(i200,4)
      dimension jeqbon(i200),eqatmj(i200)
      equivalence (xyzkpt(1,1),xkpt(1))
      dimension var(1)
      nvarp =0
      nvarq =0
      nvar  =0
      do 1000 i=1,iuniq
      if(ivary(i,4).ne.0.and.ivary(i,4).lt.8) go to 1000
      nvarq=nvarq+1
      nvar =nvar +1
      var(nvar) = qkpt(i)
 1000 continue
      do 1100 k=1,3
      do 1100 i=1,iuniq
      if(ivary(i,k).ne.0.and.ivary(i,k).lt.8) go to 1100
      nvar=nvar+1
      var(nvar) = xyzkpt(i,k)
 1100 continue
      nvarc=nvar-nvarq
c   
      write(7,1) nvarq,nvarc
    1 format(/10x,'no. of varying charges   =',i5
     &       /10x,'no of varying coordinates=',i5)
c      
c     for bond polarizability
c     
      k=0
      if(ibpair.gt.0)then
        do i=1,ibpair
          jeqbon(i)=ieqbon(i)
        enddo
        do i=1,ibpair
          if(jeqbon(i).ne.-100) then
            do j=i+1,ibpair
              if(jeqbon(i).eq.jeqbon(j)) jeqbon(j)=-100 
            enddo
            k=k+1
            ivarp(k)=jeqbon(i)
          endif
        enddo
        nvart=0
        nvaran=0
        if(iptyp.eq.0.or.iptyp.eq.1)then
          nvart=(ifit+1)*k
        write(7,5) nvart
    5   format(/10x,'no of varing transverse polarizabirity=',i5)
        endif
        if(iptyp.eq.0.or.iptyp.eq.2.or.iptyp.eq.3)then
          nvaran=(ifit+1)*k
        write(7,2) nvaran
    2   format(/10x,'no of varing anisotropy polarizabirity=',i5)
        endif
        nvarp=nvart+nvaran
        do i=1,nvarp
          nvar=nvar+1
          var(nvar) = 0.000d0  
        enddo
      endif  
c      
c     for atom polarizability
c     
      k=0
      nvara=0
      if(iatid.gt.0)then
        do i=1,iatid
          eqatmj(i)=eqatmi(i)
        enddo
        do i=1,iatid
          if(eqatmj(i).ne.-100.0d0) then
            do j=i+1,iatid
              if(eqatmj(i).eq.eqatmj(j)) eqatmj(j)=-100.0d0 
            enddo
            k=k+1
            varai(k)=eqatmj(i)
          endif
        enddo
        nvara=(ifit+1)*k
        write(7,4) nvara
    4   format(/10x,'no of varing atom polarizabirity=',i5)
        do i=1,nvara
          nvar=nvar+1
          var(nvar) = varai(i)  
        enddo
      endif  
      write(7,3) nvar
    3 format(/10x,'total no of variables (not including frozen variables!) =
     &',i5)
c   
      return
c
      entry setbck(var)
c
      jfit=ifit+1
      nrst  =0
      do 2000 i=1,iuniq
      if(ivary(i,4).ne.0.and.ivary(i,4).lt.8) go to 2000
      nrstq=nrstq+1
      nrst=nrst+1
      qkpt(i)=var(nrst)
 2000 continue
      do 2100 k=1,3
      do 2100 i=1,iuniq
      if(ivary(i,k).ne.0.and.ivary(i,k).lt.8) go to 2100
      nrst=nrst+1
      xyzkpt(i,k)=var(nrst)
 2100 continue
c
      do 2500 k=1,4
      if(neq(k).eq.0) go to 2500
      do 2510 i=1,neq(k)
      j=iateq(i,k)
                          xyzkpt(j,k)= xyzkpt(jateq(i,k),k)
 2510 if(ivary(j,k).lt.0) xyzkpt(j,k)=-xyzkpt(j         ,k)
 2500 continue
      do i=1,iuniq
       if(ivary(i,4).eq.-9.or.ivary(i,4).eq.-8)then
        qkpt(i)=qkpt(i+1)*(-1.0d0) 
       endif
      enddo
c   
      charge=0.d0
      ntst=0
      do i=1,iuniq
       if(ivary(i,4).ne.3.and.ivary(i,4).ne.-3) then
        charge=charge+qkpt(i)
       else
        if(ivary(i,4).eq.3)ntst=ntst+1
        if(ivary(i,4).eq.-3)ntst=ntst-1
       endif
      enddo
      if(ntst.ne.0)then
       qfix=(fchrg(2)-charge)/ntst
       do i=1,iuniq
        if(ivary(i,4).eq.3) qkpt(i)=qfix
        if(ivary(i,4).eq.-3) qkpt(i)=-1.0d0*qfix
       enddo
      endif
c  
      ka=0
      if(ibpair.gt.0)then
        if(iptyp.eq.1.or.iptyp.eq.2.or.iptyp.eq.3)then
          k=nvarp
          ka=ibpair
        else
          k=nvarp/2
          ka=ibpair*2
        endif
        do i=1,k,jfit
          do j=1,ibpair
            if(ieqbon(j).eq.ivarp((i-1)/jfit+1))then
              if(ifit.ge.0) pkpt(j,1)=var(nrst+1)
              if(ifit.ge.1) pkpt(j,2)=var(nrst+2)
              if(ifit.ge.2) pkpt(j,3)=var(nrst+3)
              if(ifit.ge.3) pkpt(j,4)=var(nrst+4)
            endif
          enddo
          nrst=nrst+jfit
        enddo
        if(iptyp.eq.0)then
          do i=1,k,jfit
            do j=1,ibpair
              if(ieqbon(j).eq.ivarp((i-1)/jfit+1))then
                if(ifit.ge.0) pkpt(ibpair+j,1)=var(nrst+1)
                if(ifit.ge.1) pkpt(ibpair+j,2)=var(nrst+2)
                if(ifit.ge.2) pkpt(ibpair+j,3)=var(nrst+3)
                if(ifit.ge.3) pkpt(ibpair+j,4)=var(nrst+4)
              endif
            enddo
            nrst=nrst+jfit
          enddo
        endif
      endif
c   
      if(iatid.gt.0)then
        do i=1,nvara,jfit
          do j=1,iatid
            if(eqatmi(j).eq.varai((i-1)/jfit+1))then
              if(ifit.ge.0) pkpt(ka+j,1)=var(nrst+1)
              if(ifit.ge.1) pkpt(ka+j,2)=var(nrst+2)
              if(ifit.ge.2) pkpt(ka+j,3)=var(nrst+3)
              if(ifit.ge.3) pkpt(ka+j,4)=var(nrst+4)
            endif
          enddo
          nrst=nrst+jfit
        enddo
      endif
c   
      return
      end
c**********************************************************************
      subroutine potout(var,fcal)
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/espdat/esfact,iestat,irest,ictat
      common/comp12/cxnp(i340),cynp(i340),cznp(i340),
     & fpcx(i340),fpcy(i340),fpcz(i340),
     & qupot(i340,i340,i20)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/constr/dipcal(3),qucal(3),cmas(3),dipsum
c
      common/pdpol1/ibon(i200),jbon(i200),ieqbon(i200),ibpair,iptyp,
     & nvarp,ivarp(i200),iestah,iatid,
     & iatm(i200),
     & nvara,nvart,nvaran,ishl,nvarh,
     & ifit,ifld,jfld,isym2(i340)
      common/pdpol2/
     & fpchrg(i20),pkpt(i200*2,4),
     & eqatmi(i200),varai(i200),
     & alptrn(i200),alpani(i200),alpatm(i200)
c
      common/prpol/ xmid(i200),ymid(i200),zmid(i200),
     &          trnx1(i340,i200),trny1(i340,i200),trnz1(i340,i200),
     &          trnx2(i340,i200),trny2(i340,i200),trnz2(i340,i200),
     &          bonx(i200),bony(i200),bonz(i200),bonlen(i200),
     &          fldlen(i340,i200),anifac(i340,i200),ealen(i340,i200),
     &          falen(i340,i200),
     &          signa(i340,i200),ftmp(i340),
     &          atmx1(i340,i200),atmy1(i340,i200),atmz1(i340,i200),
     &          atmx2(i340,i200),atmy2(i340,i200),atmz2(i340,i200) 
      common/enedip/ upol(i340,i20),
     &  cdpx(i340,i20),cdpy(i340,i20),cdpz(i340,i20)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/odipl/oxdip(i20,i340,i200),oydip(i20,i340,i200),
     &  ozdip(i20,i340,i200)
      dimension var(1),fcal(1)
      dimension alph(i200)
      dimension oupol(i20,i340),txdip(i20,i340),tydip(i20,i340),
     &  tzdip(i20,i340)
      character*4 dlabl(3)/ 'x  =', 'y  =', 'z  ='/
      character*5 qlabl(3)/' aa =',' bb =',' cc ='/
      data aukcal /627.509541d0/,zero/0.0d0/,debye/2.54158059d0/
c
      charge=0.d0
      do 20 i=1,iuniq
   20 charge=charge+qkpt(i)
c   
      j=0
      if(nvart.gt.0)then
        write(7,*) ' '
        write(7,*) '      parameter of transverse polarizability'
        write(7,*) '          no. bond  i -- j  (parm. trans.)'
        do i=1,ibpair
          j=j+1
          write(7,'(8x,i5,2x,i5,a3,i5,5x,a1,4f14.4,a1)') 
     &     i,ibon(i),' --',jbon(i),'(',pkpt(j,1),pkpt(j,2),
     &     pkpt(j,3),pkpt(j,4),')' 
        enddo
      endif
      if(nvaran.gt.0)then
        write(7,*) ' '
        write(7,*) '      parameter of anisotropy polarizability'
        write(7,*) '          no. bond  i -- j  (parm. aniso.)'
        do i=1,ibpair
          j=j+1
          write(7,'(8x,i5,2x,i5,a3,i5,5x,a1,4f14.4,a1,8x)') 
     &     i,ibon(i),' --',jbon(i),'(',pkpt(j,1),pkpt(j,2),
     &     pkpt(j,3),pkpt(j,4),')' 
        enddo
      endif
      if(nvara.gt.0)then
        write(7,*) ' '
        write(7,*) '      parameter of atom polarizability'
        write(7,*) '          no.     atom  simil  (parm. atom) ' 
        sump1=0.0d0
        do i=1,iatid
          j=j+1
          write(7,'(8x,i5,2x,i5,2x,f7.3,5x,a1,4f14.4,a1)') 
     &     iatm(i),iatom(i),eqatmi(i),'(',pkpt(j,1),pkpt(j,2),
     &     pkpt(j,3),pkpt(j,4),')' 
          write(8,'(8x,i5,2x,i5,2x,f7.3,5x,a1,4f14.4,a1)') 
     &     iatm(i),iatom(i),eqatmi(i),'(',pkpt(j,1),pkpt(j,2),
     &     pkpt(j,3),pkpt(j,4),')' 
           sump1=sump1+pkpt(j,1)
        enddo
        write(7,'(10x,a8,17x,f14.4)') 'sum plz1',sump1
      endif
      write(7,*) ' '
      ssq = zero
      vsq = zero
      smin=99999.d0
      smax=zero
      smaxn=0.d-5
      kmax=0
      kmaxn=0
      if(ibpair.eq.0.and.iatid.eq.0)then
        do i=1,iestat
          ssq = ssq + fcal(i)*fcal(i)
          vsq = vsq + qupot(i,1,l)*qupot(i,1,l)
          if(dabs(fcal(i)).lt.smin) smin=dabs(fcal(i))
          if(dabs(fcal(i)).gt.smax) then
            smax=dabs(fcal(i))
            kmax=i
          endif
          if(dabs(fcal(i)).lt.smax.and.dabs(fcal(i)).gt.smaxn) then
            smaxn=dabs(fcal(i))
            kmaxn=i
          endif
        enddo
      else
        kk=0
        write(7,*) '   #      sigma(kcal/mol)    s(%)   max.#    maximum 
     & (kcal/mol)'
c
        do l=1,ishl
          ssqs=zero
          vsqs=zero
          do i=1,ifld
            ssqe=zero
            vsqe=zero
            smaxe=zero
            kmaxe=0
            do j=1,iestah
              kk=kk+1
              ff =  fcal(kk)*fcal(kk)
              qq =  qupot(j,i,l)*qupot(j,i,l)
              ssqe= ssqe+ ff
              vsqe= vsqe+ qq
              if(jfld.ne.777)then
                ssqs= ssqs+ ff
                vsqs= vsqs+ qq
                ssq = ssq + ff
                vsq = vsq + qq
              else
                ssqs= ssqs+ ff*dfloat(isym2(i))
                vsqs= vsqs+ qq*dfloat(isym2(i))
                ssq = ssq + ff*dfloat(isym2(i))
                vsq = vsq + qq*dfloat(isym2(i))
              endif
              if(dabs(fcal(kk)).lt.smin) smin=dabs(fcal(kk))
              if(dabs(fcal(kk)).gt.smax) then
                smax=dabs(fcal(kk))
                kmax=kk
              endif
              if(dabs(fcal(kk)).lt.smax.and.dabs(fcal(kk)).gt.smaxn)then
                smaxn=dabs(fcal(kk))
                kmaxn=kk
              endif
              if(dabs(fcal(kk)).gt.smaxe) then
                smaxe=dabs(fcal(kk))
                kmaxe=j
              endif
            enddo
            ss = dsqrt(ssqe/dfloat(iestah))
            rr   = ss/dsqrt(vsqe/dfloat(iestah))/aukcal
            write(7,'(i5,2f15.4,i5,f15.4)') i,ss,rr*100,kmaxe,smaxe
          enddo
          if(jfld.ne.777)then
            sss = dsqrt(ssqs/dfloat(iestah*ifld))
            rrr = sss/dsqrt(vsqs/dfloat(iestah*ifld))/aukcal
          else
            sss = dsqrt(ssqs/dfloat(iestah**2))
            rrr = sss/dsqrt(vsqs/dfloat(iestah**2))/
     &            aukcal
          endif
        enddo
c
      endif
      if(jfld.ne.777)then
        sigma = dsqrt(ssq/dfloat(iestat))
        rms   = sigma/dsqrt(vsq/dfloat(iestat))/aukcal
      else
        sigma = dsqrt(ssq/dfloat(ishl*iestah**2))
        rms   = sigma/dsqrt(vsq/dfloat(ishl*iestah**2))/aukcal
      endif
c                    
      write(7,*) '  '
      write(7,1080) sigma,rms*100,smin,kmax,smax,kmaxn,smaxn
 1080 format( 10x,'standard deviation(sigma) =',f15.4,'  kcal/mol'
     &       /10x,'percent  deviation   (s)  =',f15.4,'   %      '  
     &       /10x,'minimum  deviation        =',f15.4,'  kcal/mol'  
     &       /10x,'maximum  deviation ',i6,' =',f15.4,'  kcal/mol'  
     &       /10x,'next max deviation ',i6,' =',f15.4,'  kcal/mol')  
c
      if(irest.eq.0.or.irest.eq.1)then
        write(1,*) iuniq,charge
        write(1,'(i2,f16.10,i2,f16.10,i2,f16.10,i2,f16.10,i2)') 
     &   (iatom(i),xkpt(i),ivary(i,1),ykpt(i),ivary(i,2),
     &        zkpt(i),ivary(i,3),qkpt(i),ivary(i,4),i=1,iuniq)
        write(1,*) iestat
        write(1,'(4f15.10)') 
     &   (cxnp(i),cynp(i),cznp(i),qupot(i,1,1),i=1,iestat)
c           rewind 1
      else
        kk=0
c
        do l=1,ishl 
          do i=1,ifld
            write(1,*) iuniq,charge
            write(1,'(i2,f16.10,i2,f16.10,i2,f16.10,i2,f16.10,i2)') 
     &       (iatom(k),xkpt(k),ivary(k,1),ykpt(k),ivary(k,2),
     &       zkpt(k),ivary(k,3),qkpt(k),ivary(k,4),k=1,iuniq)
            write(1,*) iestah
            do j=1,iestah
              kk=kk+1
              qq=qupot(j,i,l)*aukcal
              write(1,'(5f15.9)') 
     &         cxnp(j),cynp(j),cznp(j),qq,fcal(kk)
            enddo
          enddo
        enddo
      endif
c
c     write induction energy and induced dipole moment to rest_1
c
      do l=1,ishl
        do i=1,iestah
          upol(i,l)=zero
          cdpx(i,l)=zero
          cdpy(i,l)=zero
          cdpz(i,l)=zero
        enddo
      enddo
c   
c      bond/atom polarizability
c        transverse, induced diople length is 1.0 
c   
      jfirs=0
      if(ibpair.ne.0)then
        kk=0
        if(iptyp.eq.0.or.iptyp.eq.1)then
          jfirs=ibpair
          do l=1,ishl
            do k=1,ifld
              do j=1,ibpair
                fld1=fpchrg(l)/(fldlen(k,j)*fldlen(k,j))
                fld2=fld1*fld1
                fld3=fld2*fld1
                fld4=fld3*fld1
                fld5=fld4*fld1
                alp=pkpt(j,1)+pkpt(j,2)*fld1+pkpt(j,3)*fld2
     &             +pkpt(j,4)*fld3
                dltaq=-fld1*alp
                upol(k,l)=upol(k,l)-(0.5d0*alp*fld2)
                cdpx(k,l)=cdpx(k,l)-dltaq*trnx1(k,j)
     &                   +dltaq*trnx2(k,j)
                cdpy(k,l)=cdpy(k,l)-dltaq*trny1(k,j)
     &                   +dltaq*trny2(k,j)
                cdpz(k,l)=cdpz(k,l)-dltaq*trnz1(k,j)
     &                   +dltaq*trnz2(k,j)
              enddo
            enddo
          enddo
        endif
c   
c        anisotoropy, induced dipole length .eq. bond length
c   
        if(iptyp.eq.0.or.iptyp.eq.2.or.iptyp.eq.3)then
          kk=0
          kfirs=jfirs
          jfirs=jfirs+ibpair
          do l=1,ishl
            do k=1,ifld
              do j=1,ibpair
                jj=kfirs+j
                fld1=fpchrg(l)*anifac(k,j)/(fldlen(k,j)*fldlen(k,j))
                fld2=fld1*fld1
                fld3=fld2*fld1
                fld4=fld3*fld1
                fld5=fld4*fld1
                alp=pkpt(jj,1)+pkpt(jj,2)*fld1+pkpt(jj,3)*fld2
     &             +pkpt(jj,4)*fld3
                dltaq=-fld1*signa(k,j)*alp/bonlen(j)
                upol(k,l)=upol(k,l)-(0.5d0*fld2*alp)
                cdpx(k,l)=cdpx(k,l)-dltaq*xkpt(ibon(j))
     &                 +dltaq*xkpt(jbon(j))
                cdpy(k,l)=cdpy(k,l)-dltaq*ykpt(ibon(j))
     &                 +dltaq*ykpt(jbon(j))
                cdpz(k,l)=cdpz(k,l)-dltaq*zkpt(ibon(j))
     &                 +dltaq*zkpt(jbon(j))
              enddo
            enddo
          enddo
        endif
      endif
c   
c        atom, induced diople 
c   
      if(iatid.ne.0)then
        if(irest.eq.3.or.irest.eq.4.or.irest.eq.5.or.irest.eq.6
     &  .or.irest.eq.43.or.irest.eq.53
     &  .or.irest.eq.63.or.irest.eq.73) then

c
c         intramolecular polarization case  (irest=3 4 5 6 43 53 63 73)
c
          call intralst(jfirs,f)
        else
c   
c        atom, induced diople  
c   
          kk=0
          do l=1,ishl
            do k=1,ifld
              do j=1,iatid
                jj=jfirs+j
                fld1=fpchrg(l)/(falen(k,j)*falen(k,j))
                fld2=fld1*fld1
                fld3=fld2*fld1
                fld4=fld3*fld1
                fld5=fld4*fld1
                alp=pkpt(jj,1)+pkpt(jj,2)*fld1+pkpt(jj,3)*fld2
     &             +pkpt(jj,4)*fld3
                dltaq=-fld1*alp
                upol(k,l)=upol(k,l)-(0.5d0*alp*fld2)
                cdpx(k,l)=cdpx(k,l)-dltaq*atmx1(k,j)
     &                   +dltaq*atmx2(k,j)
                cdpy(k,l)=cdpy(k,l)-dltaq*atmy1(k,j)
     &                   +dltaq*atmy2(k,j)
                cdpz(k,l)=cdpz(k,l)-dltaq*atmz1(k,j)
     &                   +dltaq*atmz2(k,j)
              enddo
            enddo
          enddo
        endif
      endif
c
c     write induction energy and induced dipole moment
c
      if(irest.eq.3.or.irest.eq.4.or.irest.eq.5.or.irest.eq.6
     &  .or.irest.eq.43.or.irest.eq.53
     &  .or.irest.eq.63.or.irest.eq.73) then
        do l=1,ishl
          do k=1,ifld
             oupol(l,k)=zero
             txdip(l,k)=zero
             tydip(l,k)=zero
             tzdip(l,k)=zero
          enddo
        enddo
        do l=1,ishl
          do k=1,ifld
            do j=1,iatid
              oupol(l,k)=oupol(l,k)-
     &          0.5d0*(oxdip(l,k,j)*pfldx(l,k,j)+
     &                 oydip(l,k,j)*pfldy(l,k,j)+
     &                 ozdip(l,k,j)*pfldz(l,k,j))
              txdip(l,k)=txdip(l,k)+oxdip(l,k,j)
              tydip(l,k)=tydip(l,k)+oydip(l,k,j)
              tzdip(l,k)=tzdip(l,k)+ozdip(l,k,j)
c             write(7,*) l,k,j,oxdip(l,k,j),oxdip(l,k,j),oxdip(l,k,j)
            enddo
          enddo
        enddo

        write(7,*) '  '
        write(7,*) 'Induction Energy (kcal/mol) Induced dipole moment 
     &    x y z tot (debye)'
        write(7,'(2i5)') ishl,ifld
        do l=1,ishl
          do k=1,ifld
            dt=dsqrt(txdip(l,k)**2+tydip(l,k)**2+tzdip(l,k)**2)
            write(7,'(i5,f10.5,a3,4f10.5)') l*k,oupol(l,k)*aukcal,'   ',
     &        txdip(l,k)*debye,tydip(l,k)*debye,tzdip(l,k)*debye,
     &        dt*debye
          enddo
        enddo
      endif
      return
      end
!********1*********2*********3*********4*********5*********6*********7**
      subroutine matout
C
C     Write Matrix A for householder Matrix inversion 
C          
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/tens/
     &  txx(i200,i200),txy(i200,i200),txz(i200,i200),
     &  tyx(i200,i200),tyy(i200,i200),tyz(i200,i200),
     &  tzx(i200,i200),tzy(i200,i200),tzz(i200,i200)
      common/pdpol1/ibon(i200),jbon(i200),ieqbon(i200),ibpair,iptyp,
     & nvarp,ivarp(i200),iestah,iatid,
     & iatm(i200),
     & nvara,nvart,nvaran,ishl,nvarh,
     & ifit,ifld,jfld,isym2(i340)
       common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
C
      dimension A(i200*3,i200*3),F(i20*i340,i200*3),P(3,3),D(i20*i340,3)
C
C     Matrix A
C
      k=1
      do i=1,iuniq
        l=1
        do j=1,iuniq
          A(k,l)=txx(i,j)
          A(k,l+1)=txy(i,j)
          A(k,l+2)=txz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+1,l)=tyx(i,j)
          A(k+1,l+1)=tyy(i,j)
          A(k+1,l+2)=tyz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+2,l)=tzx(i,j)
          A(k+2,l+1)=tzy(i,j)
          A(k+2,l+2)=tzz(i,j)
          l=l+3
        enddo
        k=k+3
      enddo
C
C      write field from a fld charge on molecular surfaces 
C          
      m=0
      do i=1,ishl
        do j=1,ifld
c         write(9,'(180f10.5)')
c    &    (-pfldx(i,j,k),-pfldy(i,j,k),-pfldz(i,j,k),k=1,iuniq)  
          l=1
          m=m+1
          do k=1,iuniq
            F(m,l)=-pfldx(i,j,k)
            F(m,l+1)=-pfldy(i,j,k)
            F(m,l+2)=-pfldz(i,j,k)
            l=l+3
          enddo
        enddo
      enddo
      call houhol(ishl,ifld,iuniq,A,F,P,D)
      return
      end
c**********************************************************************
      subroutine sumsq(var,f)
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/espdat/esfact,iestat,irest,ictat
      common/comp12/cxnp(i340),cynp(i340),cznp(i340),
     & fpcx(i340),fpcy(i340),fpcz(i340),
     & qupot(i340,i340,i20)
c
      common/pdpol1/ibon(i200),jbon(i200),ieqbon(i200),ibpair,iptyp,
     & nvarp,ivarp(i200),iestah,iatid,
     & iatm(i200),
     & nvara,nvart,nvaran,ishl,nvarh,
     & ifit,ifld,jfld,isym2(i340)
      common/pdpol2/
     & fpchrg(i20),pkpt(i200*2,4),
     & eqatmi(i200),varai(i200),
     & alptrn(i200),alpani(i200),alpatm(i200)
c
      common/prpol/ xmid(i200),ymid(i200),zmid(i200),
     &          trnx1(i340,i200),trny1(i340,i200),trnz1(i340,i200),
     &          trnx2(i340,i200),trny2(i340,i200),trnz2(i340,i200),
     &          bonx(i200),bony(i200),bonz(i200),bonlen(i200),
     &          fldlen(i340,i200),anifac(i340,i200),ealen(i340,i200),
     &          falen(i340,i200),
     &          signa(i340,i200),ftmp(i340),
     &          atmx1(i340,i200),atmy1(i340,i200),atmz1(i340,i200),
     &          atmx2(i340,i200),atmy2(i340,i200),atmz2(i340,i200) 
      dimension var(1),f(1)
      dimension dltaq(i200)
      data aukcal /627.509541d0/
      if(ibpair.eq.0.and.iatid.eq.0)then
        do 100 i=1,iestat
        sump =0.d0
        do 110 j=1,iuniq
        xc     =fpcx(i)-xkpt(j)
        yc     =fpcy(i)-ykpt(j)
        zc     =fpcz(i)-zkpt(j)
        ra     =dsqrt(xc*xc+yc*yc+zc*zc)
  110   sump   =sump-qkpt(j)/ra
  100   f(i)   =(qupot(i,1,1) -sump)*aukcal
c                 
      else 
c   
c      calculation of electrostatic potential for induced dipoles
c   
c      bond/atom polarizability
c   
        kk=0
c
        do l=1,ishl
          do i=1,ifld
            do j=1,iestah
              kk=kk+1  
              f(kk)=qupot(j,i,l)+ftmp(j)
            enddo
          enddo
        enddo
c
c   
c        transverse, induced diople  
c   
        il=0
        jfirs=0
        if(ibpair.ne.0)then
          kk=0
          if(iptyp.eq.0.or.iptyp.eq.1)then
            jfirs=ibpair
            do j=1,ibpair
            enddo
c
            do l=1,ishl
              do k=1,ifld
                do j=1,ibpair
                  fld1=fpchrg(l)/(fldlen(k,j)*fldlen(k,j))
                  fld2=fld1*fld1
                  fld3=fld2*fld1
                  if(ifit.eq.0) then
                    alp=pkpt(j,1)
                  else if(ifit.eq.1) then
                    alp=pkpt(j,1)+pkpt(j,2)*fld1
                  else if(ifit.eq.2) then
                    alp=pkpt(j,1)+pkpt(j,2)*fld1
     &                 +pkpt(j,3)*fld2
                  else if(ifit.eq.3) then
                    alp=pkpt(j,1)+pkpt(j,2)*fld1
     &                 +pkpt(j,3)*fld2+pkpt(j,4)*fld3
                  endif
                  dltaq(j)=-fld1*alp
                enddo
                do i=1,iestah
                  kk=kk+1
                  do j=1,ibpair
                    flen1=dsqrt((cxnp(i)-trnx1(k,j))**2+
     &                          (cynp(i)-trny1(k,j))**2+
     &                          (cznp(i)-trnz1(k,j))**2)
                    flen2=dsqrt((cxnp(i)-trnx2(k,j))**2+
     &                          (cynp(i)-trny2(k,j))**2+
     &                          (cznp(i)-trnz2(k,j))**2)
                    f(kk)=f(kk)-dltaq(j)/flen1+dltaq(j)/flen2
                  enddo
                enddo
              enddo
            enddo
c
          endif
c   
c        anisotoropy, induced dipole length .eq. bond length
c   
          if(iptyp.eq.0.or.iptyp.eq.2.or.iptyp.eq.3)then
            kk=0
            kfirs=jfirs
            jfirs=jfirs+ibpair
            do j=1,ibpair
              jj=kfirs+j
            enddo
c
            do l=1,ishl
              do k=1,ifld
                do j=1,ibpair
                  jj=kfirs+j
                  fld1=fpchrg(l)*anifac(k,j)/(fldlen(k,j)*fldlen(k,j))
                  fld2=fld1*fld1
                  fld3=fld2*fld1
                  if(ifit.eq.0) then
                    alp=pkpt(jj,1)
                  else if(ifit.eq.1) then
                    alp=pkpt(jj,1)+pkpt(jj,2)*fld1
                  else if(ifit.eq.2) then
                    alp=pkpt(jj,1)+pkpt(jj,2)*fld1
     &                 +pkpt(jj,3)*fld2
                  else if(ifit.eq.3) then
                    alp=pkpt(jj,1)+pkpt(jj,2)*fld1
     &                 +pkpt(jj,3)*fld2+pkpt(jj,4)*fld3
                  endif
                  dltaq(j)=-fld1*signa(k,j)*alp/bonlen(j)
                enddo
                do i=1,iestah
                  kk=kk+1
                  do j=1,ibpair
                    flen1=ealen(i,ibon(j))
                    flen2=ealen(i,jbon(j))
                    f(kk)=f(kk)-dltaq(j)/flen1+dltaq(j)/flen2
                  enddo
                enddo
              enddo
            enddo
c
          endif
        endif
c   
c        atom, induced diople  
c   
        if(iatid.ne.0)then
c
          if(irest.eq.3.or.irest.eq.4.or.irest.eq.5.or.irest.eq.6
     &    .or.irest.eq.43.or.irest.eq.53
     &    .or.irest.eq.63.or.irest.eq.73) then
c
c         intramolecular polarization case  (irest=3 or 4 or 5 or 6 or 43 or 53 or 63 or 73)
c
            call intrapop(jfirs,f)
          else
c   
c        atom, induced diople  
c   
            kk=0
            do l=1,ishl
              do k=1,ifld
                do j=1,iatid
                  jj=jfirs+j
                  fld1=fpchrg(l)/(falen(k,j)*falen(k,j))
                  fld2=fld1*fld1
                  fld3=fld2*fld1
                  if(ifit.eq.0) then
                    alp=pkpt(jj,1)
                  else if(ifit.eq.1) then
                    alp=pkpt(jj,1)+pkpt(jj,2)*fld1
                  else if(ifit.eq.2) then
                    alp=pkpt(jj,1)+pkpt(jj,2)*fld1
     &                 +pkpt(jj,3)*fld2
                  else if(ifit.eq.3) then
                    alp=pkpt(jj,1)+pkpt(jj,2)*fld1
     &                 +pkpt(jj,3)*fld2+pkpt(jj,4)*fld3
                  endif
                  dltaq(j)=-fld1*alp
                enddo
                do i=1,iestah
                  kk=kk+1
                  do j=1,iatid
                    flen1=dsqrt((cxnp(i)-atmx1(k,j))**2+
     &                          (cynp(i)-atmy1(k,j))**2+
     &                          (cznp(i)-atmz1(k,j))**2)
                    flen2=dsqrt((cxnp(i)-atmx2(k,j))**2+
     &                          (cynp(i)-atmy2(k,j))**2+
     &                          (cznp(i)-atmz2(k,j))**2)
                    f(kk)=f(kk)-dltaq(j)/flen1+dltaq(j)/flen2
                  enddo
                enddo
              enddo
            enddo
          endif
c   
        endif
        kk=0
c
        do l=1,ishl
          do k=1,ifld
            do i=1,iestah
              kk=kk+1
              f(kk)=f(kk)*aukcal
*             write(7,*) '### sumsq ##',kk,f(kk)
            enddo
          enddo
        enddo
c
      endif
      return
      end
c**********************************************************************
      subroutine optssq(x,f,xjac,xjtj,work,w1,w2,w3,w4,w5,f1,f2)
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/espdat/esfact,iestat,irest,ictat
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt(i200)
c
      common/pdpol1/ibon(i200),jbon(i200),ieqbon(i200),ibpair,iptyp,
     & nvarp,ivarp(i200),iestah,iatid,
     & iatm(i200),
     & nvara,nvart,nvaran,ishl,nvarh,
     & ifit,ifld,jfld,isym2(i340)
      common/pdpol2/
     & fpchrg(i20),pkpt(i200*2,4),
     & eqatmi(i200),varai(i200),
     & alptrn(i200),alpani(i200),alpatm(i200)
c
      real*8 parm(4)/1.2d0,2.d0,120.d0,1.d-4/
      data ax/0.1d0/,maxfn/2000/
      dimension x(1),f(1),xjac(iestat,1),xjtj(1),work(1)
      dimension w1(1),w2(1),w3(1),w4(1),w5(1),f1(1),f2(1)
c
c     ----- machine dependent constants -----
c
      data ieps/9/,idelta/5/,nsig/6/,ksig/16/
      eps    = 10.d0**(-ieps)
      delta  = 10.d0**(-idelta)
      prec   = 10.d0**(-ksig-1)
      rel    = 10.d0**(-ksig/2)
      relcon = 10.d0**(-nsig)
c
      if(nvar.le.0.or.parm(2).le.1.d0.or.parm(1).le.0.d0) go to 9305
c
      ier = 0
      al     = parm(1)
      erl2 = 1.d10
      ibad = -99
      isw = 1
      infer = 0
      ier = 0
      call setbck(x)
      call sumsq (x,f1)
      ieval = 1
      ssqold = 0.d0
      if(jfld.ne.777)then
        do i=1,iestat
          ssqold = ssqold+f1(i)*f1(i)
        enddo
      else
        ix=0
        do i=1,ishl
          do j=1,ifld
            xtmp=0.0d0
            do k=1,iestah
              ix=ix+1
              xtmp = xtmp+f1(ix)*f1(ix)
            enddo
            xtmp=xtmp*dfloat(isym2(j))
            ssqold = ssqold+xtmp
          enddo
        enddo
      endif
      write(7,3) ieval,ssqold,(x(i),i=1,min0(nvar,12))
    3 format(//'  ieval',7x,'sum.sq',8x,'norm of grad.',6x,'variables'
     &       /i5,f17.7,22x,12f11.5)
                                                     go to 1100
c
c********** main iteration loop **********
c
 1000 ssqold = ssq
      if(infer.gt.0.or.ijac.ge.nvar.or.icount.gt.0) go to 1100
c
c     ----- rank one update to jacobian -----
c
      ijac = ijac+1
      dsq = 0.d0
      do 35 j=1,nvar
   35 dsq = dsq+w2(j)*w2(j)
      if(dsq.le.0.d0) go to 1100
      if(jfld.ne.777)then
        do i=1,iestat
          g = f(i)-f2(i)
          do j=1,nvar
            g = g+xjac(i,j)*w2(j)
          enddo
          g = g/dsq
          do j=1,nvar
            xjac(i,j) = xjac(i,j)-g*w2(j)
          enddo
        enddo
      else
        ix=0
        do i=1,ishl
          do jj=1,ifld
            do k=1,iestah
              ix=ix+1
              g = f(ix)*dfloat(isym2(jj))-f2(ix)*dfloat(isym2(jj))
              do j=1,nvar
                g = g+xjac(ix,j)*w2(j)
              enddo
              g = g/dsq
              do j=1,nvar
                xjac(ix,j) = xjac(ix,j)-g*w2(j)
              enddo
            enddo
          enddo
        enddo
      endif
                                                     go to 1200
c
c     ----- jacobian by incrementing x -----
c
 1100 ijac = 0
      do 75 j=1,nvar
      xdabs = dabs(x(j))
      hh = rel*(dmax1(xdabs,ax))
      xhold = x(j)
      x(j) = x(j)+hh
      call setbck(x)
      call sumsq (x,f1)
      ieval = ieval+1
      x(j) = xhold
      if(isw.eq.1) go to 65
c
c     ----- central differences -----
c
      x(j) = xhold-hh
      call setbck(x)
      call sumsq (x,f2)
      ieval = ieval+1
      x(j) = xhold
      rhh = 0.5d0/hh
      if(jfld.ne.777)then
        do i=1,iestat
          xjac(i,j) = (f1(i)-f2(i))*rhh
        enddo
      else
        ix=0
        do i=1,ishl
          do jj=1,ifld
            do k=1,iestah
              ix=ix+1
              xjac(ix,j) = (f1(ix)*dfloat(isym2(jj))-
     &                      f2(ix)*dfloat(isym2(jj)))*rhh
            enddo
          enddo
        enddo
      endif
      go to 75
c
c     ----- forward differences -----
c
   65 rhh = 1.d0/hh
      if(jfld.ne.777)then
        do i=1,iestat
          xjac(i,j) = (f1(i)-f(i))*rhh
        enddo
      else
        ix=0
        do i=1,ishl
          do jj=1,ifld
            do k=1,iestah
              ix=ix+1
              xjac(ix,j) = (f1(ix)*dfloat(isym2(jj))
     &                     -f(ix)*dfloat(isym2(jj)))*rhh
            enddo
          enddo
        enddo
      endif
   75 continue
c
c     ----- calculate gradient -----
c
 1200 erl2x = erl2
      erl2 = 0.d0
      sumx=0.0d0
      sumy=0.0d0
      do j=1,nvar
        if(jfld.ne.777)then
          do i=1,iestat
            sumx=sumx+xjac(i,j)
            sumy=sumy+f(i)*xjac(i,j)
          enddo
        else
          ix=0
          do i=1,ishl
            do jj=1,ifld
              do k=1,iestah
                ix=ix+1
                sumx = sumx+xjac(ix,j)
                sumy = sumy+f(ix)*xjac(ix,j)
              enddo
            enddo
          enddo
        endif
      enddo
      sumx=0.0d0
      if(jfld.ne.777)then
        do i=1,iestat
          sumx=sumx+f(i)
        enddo
      else
        ix=0
        do i=1,ishl
          do jj=1,ifld
            do k=1,iestah
              ix=ix+1
              sumx = sumx+f(ix)*dfloat(isym2(jj))
            enddo
          enddo
        enddo
      endif
      do 90 j=1,nvar
      sum = 0.d0
      if(jfld.ne.777)then
        do i=1,iestat
          sum = sum+xjac(i,j)*f(i)
        enddo
      else
        ix=0
        do i=1,ishl
          do jj=1,ifld
            do k=1,iestah
              ix=ix+1
              sum = sum+xjac(ix,j)*f(ix)
            enddo
          enddo
        enddo
      endif
      w1(j) = sum
   90 erl2 = erl2+sum*sum
      erl2 = dsqrt(erl2)
c
c    -----  convergence test for norm of gradient -----
c
      if(ijac  .gt.  0    ) go to 95
      if(erl2*2.le.delta  ) infer = infer+4
      if(erl2*2.le.parm(4)) isw = 2
c
c     ----- calculate the lower super triange of
c           jacobian (transposed) * jacobian -----
c
   95 l=0
      do 100 i=1,nvar
      do 100 j=1,i
      l=l+1
      sum = 0.d0
      if(jfld.ne.777)then
        do k=1,iestat
          sum = sum+xjac(k,i)*xjac(k,j)
        enddo
      else
        ix=0
        do ii=1,ishl
          do jj=1,ifld
            do k=1,iestah
              ix=ix+1
              sum = sum+xjac(ix,i)*xjac(ix,j)/dfloat(isym2(jj))
            enddo
          enddo
        enddo
      endif
  100 xjtj(l) = sum
c
c     ----- convergence checks -----
c
      if(infer.gt.  0  )                                  go to 9000
      if(ieval.ge.maxfn)                                  go to 9290
c
c     ----- compute scaling vector -----
c
      k=0
      do 115 j=1,nvar
      k=k+j
  115 w3(j) = xjtj(k)
c
c     ----- add l-m factor to diagonal -----
c
      icount = 0
 1300 k = 0
      do 150 i=1,nvar
      do 145 j=1,i
      k = k+1
  145 work(k) = xjtj(k)
      work(k) = work(k)+w3(i)*al
  150 w2(i) = w1(i)
c
c     ----- cholesky decomposition -----
c
c           (1) decompose a -----
c
      call ludecp (work,work,nvar,g,xhold,ier)
      if(ier .eq.0) go to 160
      ier = 0
      if(ijac.gt.0)                                        go to 1100
      if(ibad.le.0)                                        go to 1900
      if(ibad.ge.2)                                        go to 9310
                                                           go to 1400
c
c           (2) elimination -----
c
  160 call luelmp (work,w2,nvar,w2)
      if(ibad.ne.-99) ibad = 0
c
c     ----- calculate sum of squares -----
c
      do 170 j=1,nvar
  170 w4(j) = x(j)-w2(j)
      call setbck(w4)
      call sumsq (w4,f1)
      ieval = ieval+1
      ssq = 0.d0
      if(jfld.ne.777)then
      do i=1,iestat
        ssq = ssq+f1(i)*f1(i)
      enddo
      else
        ix=0
        do i=1,ishl
          do j=1,ifld
            xtmp=0.0d0
            do k=1,iestah
              ix=ix+1
              xtmp = xtmp+f1(ix)*f1(ix)
            enddo
            xtmp=xtmp*dfloat(isym2(j))
            ssq = ssq+xtmp
          enddo
        enddo
      endif
      if(ssq.le.ssqold*1.000001d0) go to 205
c
c     ----- increase parameter and try again -----
c
 1400 icount = icount+1
      al = al*parm(2)**2
      if(ijac  .eq.0                 ) go to  195
      if(icount.ge.4.or.al.gt.parm(3)) go to  200
  195 if(               al.le.parm(3))                     go to 1300
      if(ibad  .eq.1                 )                     go to 9310
                                                           go to 9300
  200 al = al/parm(2)**8
                                                           go to 1100
c
c     ----- adjust marquardt parameter -----
c
  205 if(icount.eq.0) al = al/parm(2)
      if(erl2x.le.0.d0) go to 210
      g = erl2/erl2x
      if(erl2.lt.erl2x) al = al*dmax1(1.d0/parm(2),g)
      if(erl2.gt.erl2x) al = al*dmin1(parm(2),g)
  210 al = dmax1(al,prec)
c
c************ one iteration cycle completed ************
c
      do 220 j=1,nvar
  220 x(j) = w4(j)
      if(jfld.ne.777)then
        do i=1,iestat
          f2(i) = f (i)
          f (i) = f1(i)
        enddo
      else
        ix=0
        do i=1,ishl
          do jj=1,ifld
            do k=1,iestah
              ix=ix+1
              f2(ix) = f (ix)
              f (ix) = f1(ix)
            enddo
          enddo
        enddo
      endif
c
      write(7,1) ieval,ssq,erl2*2,(x(i),i=1,min0(nvar,12))
    1 format(i5,f17.7,f19.6,3x,12f11.5)
c
c     ----- relative convergence test for x -----
c
      if(icount.gt.0.or.ijac.gt.0)                        go to 1000
      xdifmx=0.d0
      do 230 j=1,nvar
      xdif = dabs(w2(j))/dmax1(dabs(x(j)),ax)
      xdifmx=dmax1(xdifmx,xdif)
      if(xdif.gt.relcon) go to 235
  230 continue
      infer = 1
c
c     ----- relative convergence test for ssq -----
c
  235 sqdif = dabs(ssq-ssqold)/dmax1(ssqold,ax)
      if(sqdif.le.eps) infer = infer+2
                                                           go to 1000
c
c     ----- singular decomposition -----
c
 1900 if(ibad) 1910,1920,1930
c
c     ----- check to see if current iterate has cycled back to
c           the last singular point -----
c
 1920 do 250 j=1,nvar
      xhold = w5(j)
      if(dabs(x(j)-xhold).gt.relcon*dmax1(ax,dabs(xhold)))go to 1910
  250 continue
                                                           go to 9295
c
c     ----- update the bad x values -----
c
 1910 do 260 j=1,nvar
  260 w5(j) = x(j)
      ibad = 1
c
c     ----- replace zeroes on hessian diagonal -----
c
 1930 izro = 0
      do 285 j=1,nvar
      if(w3(j).gt.0.d0) go to 285
      izro = izro+1
      w3(j) = 1.d0
  285 continue
      if(izro.lt.nvar)                                     go to 1300
                                                           ier = 38
                                                           go to 9000
c
c     ----- error termination -----
c
 9290 ier = ier+1
 9295 ier = ier+1
 9300 ier = ier+1
 9305 ier = ier+1
 9310 ier = ier+129
 9000 continue
           if(ier.eq.129) then
                                 write(7, 8)
      else if(ier.eq.130) then
                                 write(7,18)
                                 return
      else if(ier.eq.131) then
                                 write(7,28)
C%SN
                  write(7,'(a21,f10.5)') '          al > 120.0 ',al
C%SN
      else if(ier.eq.132) then
                                 write(7,38)
      else if(ier.ne. 0 ) then
                                 write(7,48) ier
                          endif
      write(7,2) nsig,xdifmx,ieps,sqdif,idelta,erl2*2
c
    2 format(//28x,'*** results of levenberg-marquardt algorithm ***',
     *//22x,'convergence limit for the variables =',' 1.d-0',i1,e15.3
     * /22x,'convergence limit for the ssq       =',' 1.d-0',i1,e15.3
     * /22x,'convergence limit for the norm of g =',' 1.d-0',i1,e15.3)
    8 format(/10x,'failed to recover from singular jacobian')
   18 format(/10x,'one of m,n,parm,is incorrect')
   28 format(/10x,'marquardt parameter exceeded')
   38 format(/10x,'the solution cycled backto singularity')
   48 format(/10x,'error indicating parameter      =',i5)
      return
      end
c**********************************************************************
      subroutine luelmp (a,b,n,x)
      implicit real*8(a-h,o-z)
c
      dimension a(1),b(1),x(1)
c
c     ----- solution of ly = b -----
c
      ip=1
      iw = 0
      do 15 i=1,n
      t=b(i)
      im1 = i-1
      if(iw .eq. 0) go to 9
      ip=ip+iw-1
      do 5 k=iw,im1
      t = t-a(ip)*x(k)
      ip=ip+1
    5 continue
      go to 10
    9 if(t .ne. 0.d0) iw = i
      ip = ip+im1
   10 x(i)=t*a(ip)
      ip=ip+1
   15 continue
c
c     ----- solution of ux = y -----
c
      n1 = n+1
      do 30 i = 1,n
      ii = n1-i
      ip=ip-1
      is=ip
      iq=ii+1
      t=x(ii)
      if(n.lt.iq) go to 25
      kk = n
      do 20 k=iq,n
      t = t - a(is) * x(kk)
      kk = kk-1
      is = is-kk
   20 continue
   25 x(ii)=t*a(is)
   30 continue
      return
      end
c**********************************************************************
      subroutine ludecp (a,ul,n,d1,d2,ier)
      implicit real*8(a-h,o-z)
c
      dimension a(1),ul(1)
c
      kount=0
      d1=1.d0
      d2=0.d0
      rn = 1.d0/(float(n)*16.d0)
      ip = 1
      ier=0
      do 45 i = 1,n
      iq = ip
      ir = 1
      do 40 j = 1,i
      x = a(ip)
      if(j .eq. 1) go to 10
      do 5  k=iq,ip1
      x = x - ul(k) * ul(ir)
      ir = ir+1
    5 continue
   10 if(i.ne.j) go to 30
      d1 = d1*x
      if(a(ip) + x*rn .le. a(ip)) go to 50
   15 if(dabs(d1).le.1.d0) go to 20
      d1 = d1 * 0.0625d0
      d2 = d2 + 4.d0
      go to 15
   20 if(dabs(d1) .ge. 0.0625d0) go to 25
      d1 = d1 * 16.d0
      d2 = d2 - 4.d0
      go to 20
   25 ul(ip) = 1.d0/dsqrt(x)
      go to 35
   30 ul(ip) = x * ul(ir)
   35 ip1 = ip
      ip = ip+1
      ir = ir+1
   40 continue
   45 continue
      go to 60
   50 ier = 129
   60 continue
      return
      end
c**********************************************************************
      subroutine prepol
c   
c     preparation for bond polarizability fit
c   
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/espdat/esfact,iestat,irest,ictat
      common/comp12/cxnp(i340),cynp(i340),cznp(i340),
     & fpcx(i340),fpcy(i340),fpcz(i340),
     & qupot(i340,i340,i20)
c
      common/pdpol1/ibon(i200),jbon(i200),ieqbon(i200),ibpair,iptyp,
     & nvarp,ivarp(i200),iestah,iatid,
     & iatm(i200),
     & nvara,nvart,nvaran,ishl,nvarh,
     & ifit,ifld,jfld,isym2(i340)
      common/pdpol2/
     & fpchrg(i20),pkpt(i200*2,4),
     & eqatmi(i200),varai(i200),
     & alptrn(i200),alpani(i200),alpatm(i200)
c
      common/prpol/ xmid(i200),ymid(i200),zmid(i200),
     &          trnx1(i340,i200),trny1(i340,i200),trnz1(i340,i200),
     &          trnx2(i340,i200),trny2(i340,i200),trnz2(i340,i200),
     &          bonx(i200),bony(i200),bonz(i200),bonlen(i200),
     &          fldlen(i340,i200),anifac(i340,i200),ealen(i340,i200),
     &          falen(i340,i200),
     &          signa(i340,i200),ftmp(i340),
     &          atmx1(i340,i200),atmy1(i340,i200),atmz1(i340,i200),
     &          atmx2(i340,i200),atmy2(i340,i200),atmz2(i340,i200) 
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      data zero/0.0d0/
c           
c     midpoint of bond   xmid,ymid,zmid
c  
      do i=1,ibpair
        xmid(i)=0.5d0*(xkpt(ibon(i))+xkpt(jbon(i)))
        ymid(i)=0.5d0*(ykpt(ibon(i))+ykpt(jbon(i)))
        zmid(i)=0.5d0*(zkpt(ibon(i))+zkpt(jbon(i)))
      enddo
c  
c     electrostatic potential for fixed parameters
c  
      do i=1,iestah
        sump=zero
        do j=1,iuniq
          xc     =cxnp(i)-xkpt(j)
          yc     =cynp(i)-ykpt(j)
          zc     =cznp(i)-zkpt(j)
          ealen(i,j)=dsqrt(xc*xc+yc*yc+zc*zc)
          sump   =sump-qkpt(j)/ealen(i,j)
        enddo
        ftmp(i)   = -sump
      enddo
c  
c     atom atom distances and intramolecular charge field
c  
      if(irest.eq.3.or.irest.eq.4.or.irest.eq.5.or.irest.eq.6
     &  .or.irest.eq.43.or.irest.eq.53
     &  .or.irest.eq.63.or.irest.eq.73) then
        do i=1,iuniq
          do j=1,iuniq
            ax(i,j)=xkpt(i)-xkpt(j)
            ay(i,j)=ykpt(i)-ykpt(j)
            az(i,j)=zkpt(i)-zkpt(j)
            ar(i,j)=sqrt(ax(i,j)**2+ay(i,j)**2+az(i,j)**2)
          enddo
        enddo
        call imcfld
        do l=1,ishl
          do k=1,ifld
            do j=1,iuniq
                fldx0(l,k,j)=zero
                fldy0(l,k,j)=zero
                fldz0(l,k,j)=zero
                pfldx(l,k,j)=cfldx(j)
                pfldy(l,k,j)=cfldy(j)
                pfldz(l,k,j)=cfldz(j)
            enddo
          enddo
        enddo
      endif
      if(irest.eq.3.or.irest.eq.4.or.irest.eq.5.or.irest.eq.6
     &  .or.irest.eq.43.or.irest.eq.53
     &  .or.irest.eq.63.or.irest.eq.73) then
        do l=1,ishl
          do k=1,ifld
            do j=1,iuniq
              xc     =fpcx(k)-xkpt(j)
              yc     =fpcy(k)-ykpt(j)
              zc     =fpcz(k)-zkpt(j)
              falen(k,j)=dsqrt(xc*xc+yc*yc+zc*zc)
              r3=falen(k,j)*falen(k,j)*falen(k,j)
              fldx0(l,k,j)=fldx0(l,k,j)+fpchrg(l)*xc/r3
              fldy0(l,k,j)=fldy0(l,k,j)+fpchrg(l)*yc/r3
              fldz0(l,k,j)=fldz0(l,k,j)+fpchrg(l)*zc/r3
              pfldx(l,k,j)=pfldx(l,k,j)+fpchrg(l)*xc/r3
              pfldy(l,k,j)=pfldy(l,k,j)+fpchrg(l)*yc/r3
              pfldz(l,k,j)=pfldz(l,k,j)+fpchrg(l)*zc/r3
            enddo
          enddo
        enddo
c
      else
        do i=1,ifld
        sump=0.d0
          do j=1,iuniq
            xc     =fpcx(i)-xkpt(j)
            yc     =fpcy(i)-ykpt(j)
            zc     =fpcz(i)-zkpt(j)
            falen(i,j)=dsqrt(xc*xc+yc*yc+zc*zc)
            fldl2=falen(i,j)*2.d0
            atmx1(i,j)=xkpt(j)-xc/fldl2
            atmx2(i,j)=xkpt(j)+xc/fldl2
            atmy1(i,j)=ykpt(j)-yc/fldl2
            atmy2(i,j)=ykpt(j)+yc/fldl2
            atmz1(i,j)=zkpt(j)-zc/fldl2
            atmz2(i,j)=zkpt(j)+zc/fldl2
          enddo
        enddo
      endif
c  
c     bond length 
c  
      do i=1,ibpair
        bonx(i)=xkpt(ibon(i))-xkpt(jbon(i))
        bony(i)=ykpt(ibon(i))-ykpt(jbon(i))
        bonz(i)=zkpt(ibon(i))-zkpt(jbon(i))
        bonlen(i)=dsqrt(bonx(i)**2+bony(i)**2+bonz(i)**2)
      enddo
c   
c       electric field length & anisotoropy factor & 
c       x y z coordinates of transverse induced dipoles
c   
      do i=1,ifld
        do j=1,ibpair
          fldx=fpcx(i)-xmid(j)
          fldy=fpcy(i)-ymid(j)
          fldz=fpcz(i)-zmid(j)
          fldlen(i,j)=dsqrt(fldx**2+fldy**2+fldz**2)
          if(falen(i,ibon(j)).ge.falen(i,jbon(j)))
     &                                               then
            anifac(i,j)=(-bonx(j)*fldx-bony(j)*fldy
     &               -bonz(j)*fldz)/(fldlen(i,j)*bonlen(j))
            signa(i,j)=1.0d0
          else
            anifac(i,j)=(bonx(j)*fldx+bony(j)*fldy
     &               +bonz(j)*fldz)/(fldlen(i,j)*bonlen(j))
            signa(i,j)=-1.0d0
          endif
          fldl2=fldlen(i,j)*2.d0
          trnx1(i,j)=xmid(j)-fldx/fldl2
          trnx2(i,j)=xmid(j)+fldx/fldl2
          trny1(i,j)=ymid(j)-fldy/fldl2
          trny2(i,j)=ymid(j)+fldy/fldl2
          trnz1(i,j)=zmid(j)-fldz/fldl2
          trnz2(i,j)=zmid(j)+fldz/fldl2
        enddo
      enddo
      return
      end
***********************************************************************
      subroutine lnkplz
c
c     preparation for connectivity and optimized esp charges
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/espdat/esfact,iestat,irest,ictat
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      character*4 title(20)
      dimension lnk1234(i200,i200,3),ln(90),kk(i200),k3(i200),kk2(i200)
      dimension lnk12bf(i200,i200,3)
      data zero/0.0d0/,one/1.0d0/
C  
C     initial clear or set & distance check
C
      read(4,'(20a4)') title
      write(6,*) title
      read(4,*) allcc,alldd,cutc,cutd,factt
      write(6,*) allcc,alldd,cutc,cutd,factt
      numc=0
      numd=0
      do i=1,iuniq
        do j=1,iuniq
          flnkc(i,j)=allcc
          flnkd(i,j)=alldd
          if(i.eq.j)then
            flnkc(i,j)=zero
            flnkd(i,j)=zero
          endif
          dis=dsqrt((xkpt(i)-xkpt(j))**2+(ykpt(i)-ykpt(j))**2+
     &        (zkpt(i)-zkpt(j))**2)
          if(dis.gt.cutc) then
            flnkc(i,j)=one
            numc=numc+1
          endif
          if(dis.gt.cutd) then
            flnkd(i,j)=zero
            numd=numd+1
          endif
        enddo
      enddo
      write(7,*) '   '
      write(7,'(a31,f6.2,a34,f7.2,a12,i6)')
     &  '  initial clear charge         ',
     &  allcc,' distance > cut distance  set 1.0 ',cutc,
     &  ' cut number ',numc
      write(7,'(a31,f6.2,a34,f7.2,a12,i6)')
     &  '  initial clear induced dipole ',
     &  alldd,' distance > cut distance  set 0.0 ',cutd,
     &  ' cut number ',numd
      write(7,*) '   '
C
C     read 1-2 bond data
C
      read(4,'(20a4)') title
      read(4,*) facc(1),facd(1)
      write(7,*) '   '
      write(7,'(a34,2f7.4)') '  s-factor charge & induced dipole',
     &  facc(1),facd(1)
      write(7,*) '   '
      do i=1,iuniq
        read(4,'(200i4)') (lnk1234(i,j,1),j=1,iuniq)
      enddo
      do i=1,iuniq
        kk(i)=0
        do j=1,iuniq
          if(lnk1234(i,j,1).ne.0)then
            kk(i)=kk(i)+1
            ln(kk(i))=lnk1234(i,j,1)
          endif
        enddo
        write(7,'(i4,a5,90i4)') kk(i),'%1-2 ',(ln(k),k=1,kk(i))
      enddo
C
C     make for & back 1-2 bond data
C
      do i=1,iuniq
        kk2(i)=kk(i)
        do j=1,iuniq
          lnk12bf(i,j,1)=lnk1234(i,j,1)
        enddo
      enddo
      do i=1,iuniq
        do j=2,iuniq
          if(lnk1234(i,j,1).ne.0)then
            k2=lnk1234(i,j,1)
            kk2(k2)=kk2(k2)+1
            lnk12bf(k2,kk2(k2),1)=lnk1234(i,1,1)
          endif
        enddo
      enddo
C     do i=1,iuniq
C       write(7,'(i4,a5,90i4)') kk2(i),'%1-2b',
C    &     (lnk12bf(i,k,1),k=1,kk2(i))
C     enddo
C
C     make 1-3 data from 1-2 
C
      do i=1,iuniq
        lnk12bf(i,1,2)=i
        lnk1234(i,1,2)=i
        do j=2,iuniq
          lnk12bf(i,j,2)=0
          lnk1234(i,j,2)=0
        enddo
      enddo
C
      do i=1,iuniq
        k3(i)=1
      enddo
      do i=1,iuniq
        do j=2,iuniq
          if(lnk12bf(i,j,1).ne.0)then
            do k=1,iuniq
              if(lnk12bf(i,j,1).eq.lnk12bf(k,1,1))then
                do l=2,kk2(k)
                  if(lnk12bf(k,l,1).ne.lnk12bf(i,1,1))then
                    k3(i)=k3(i)+1
                    lnk12bf(i,k3(i),2)=lnk12bf(k,l,1)
                  endif
                enddo
              endif
            enddo
          endif
        enddo
      enddo
C
      do i=1,iuniq
        kka=1
        do j=2,iuniq
           if(lnk12bf(i,j,2).gt.lnk12bf(i,1,2))then
             kka=kka+1
             lnk1234(i,kka,2)=lnk12bf(i,j,2)
           endif
        enddo
      enddo
C
      do i=1,iuniq
        kka=0
        do j=1,iuniq
          if(lnk1234(i,j,2).ne.0)then
            kka=kka+1
            ln(kka)=lnk1234(i,j,2)
          endif
        enddo
      write(7,'(i4,a5,90i4)') kka,'%1-3 ',(lnk1234(i,j,2),j=1,kka)
      enddo
C
C     make 1-4 data from 1-3 
C
      do i=1,iuniq
        lnk12bf(i,1,3)=i
        lnk1234(i,1,3)=i
        do j=2,iuniq
          lnk12bf(i,j,3)=0
          lnk1234(i,j,3)=0
        enddo
      enddo
C
      do i=1,iuniq
        k3(i)=1
      enddo
      do i=1,iuniq
        do j=2,iuniq
          if(lnk12bf(i,j,2).ne.0)then
            do k=1,iuniq
              if(lnk12bf(i,j,2).eq.lnk12bf(k,1,1))then
                do l=2,kk2(k)
                  mk=0
                  do m=2,iuniq
                    if(lnk12bf(k,l,1).eq.lnk12bf(i,m,1))then
                      mk=1
                    endif
                  enddo
                  if(mk.eq.0)then
                    k3(i)=k3(i)+1
                    lnk12bf(i,k3(i),3)=lnk12bf(k,l,1)
                  endif
                enddo
              endif
            enddo
          endif
        enddo
      enddo
C
      do i=1,iuniq
        kka=1
        do j=2,iuniq
           if(lnk12bf(i,j,3).gt.lnk12bf(i,1,3))then
             kka=kka+1
             lnk1234(i,kka,3)=lnk12bf(i,j,3)
           endif
        enddo
      enddo
C
      do i=1,iuniq
        kka=0
        do j=1,iuniq
C
C         for 5 menber ring
C
          if(j.ne.1)then
            do k=2,iuniq
              if(lnk1234(i,j,3).eq.lnk1234(i,k,2)) lnk1234(i,j,3)=0
            enddo
          endif
C
C         for 6 menber ring
C
          do k=j+1,iuniq
            if(lnk1234(i,j,3).eq.lnk1234(i,k,3)) lnk1234(i,k,3)=0
          enddo
C%
          if(lnk1234(i,j,3).ne.0)then
            kka=kka+1
            ln(kka)=lnk1234(i,j,3)
          endif
        enddo
C%
        do j=1,iuniq
          lnk1234(i,j,3)=0
        enddo
        do j=1,kka
          lnk1234(i,j,3)=ln(j)
        enddo
C%
      write(7,'(i4,a5,90i4)') kka,'%1-4 ',(lnk1234(i,j,3),j=1,kka)
      enddo
C
C     set scaling 1   1-2
C                 2   1-3
C                 3   1-4
C
      do k=2,3
        read(4,'(20a4)') title
        read(4,*) facc(k),facd(k)
        write(7,*) '   '
        write(7,'(a34,2f7.4)') '  s-factor charge & induced dipole',
     &    facc(k),facd(k)
        write(7,*) '   '
      enddo
      do k=3,1,-1
        do i=1,iuniq
          do j=2,iuniq
            if(lnk1234(i,j,k).ne.0)then
              flnkc(lnk1234(i,1,k),lnk1234(i,j,k))=facc(k)
              flnkd(lnk1234(i,1,k),lnk1234(i,j,k))=facd(k)
            endif
          enddo
        enddo
        do i=1,iuniq
          do j=i,iuniq
            flnkc(j,i)=flnkc(i,j)
            flnkd(j,i)=flnkd(i,j)
          enddo
        enddo
      enddo
C
C     atomic charge read
C
      read(4,'(20a4)') title
      pchtot=zero
      do i=1,iuniq
        read(4,'(i6,3x,f12.6)') lnk,pch
        pchin(lnk)=pch
        pchtot=pchtot+pch
      enddo
      if(pchtot.ne.zero)then
        write(7,*) 'total of esp charges is not zero! ',pchtot
      endif
      if(allcc.ne.zero.or.facc(1).ne.zero.or.
     & facc(2).ne.zero.or.facc(3).ne.zero)then
        write(7,*) '  '
        write(7,*) 
     &  '  scale factor for intramolecular charge interaction '
        do i=1,iuniq
          write(7,'(a5,200f5.2)') '  %c ',(flnkc(i,j),j=1,iuniq)
        enddo
      endif
      write(7,*) '  '
c
c     set default dumping value a
c
      if(irest.eq.4.or.irest.eq.5.or.irest.eq.6 
     &   .or.irest.eq.43.or.irest.eq.53
     &   .or.irest.eq.63.or.irest.eq.73) then
        if(factt.lt.zero)then
          if(irest.eq.4.or.irest.eq.43) factt=1.9088d0
          if(irest.eq.5.or.irest.eq.53) factt=0.39d0
          if(irest.eq.6)  factt=0.8456d0
          if(irest.eq.63) factt=0.8456d0
          if(irest.eq.73) factt=1.7823d0
        endif
        write(7,*) '  '
        write(7,*) 'irest=',irest,' % Dumping Parameter a=',
     & factt
      else
        write(7,*) '  '
        write(7,*) 
     & '  scale factor for intramolecular induced dipole interaction'
        iiu1=1
        iiu30=30
        numr=iuniq/iiu30
        write(7,*) numr
        numx=iuniq-numr*30
        do ii=1,numr
          write(7,'(i5,a5,i5)') iiu1,' --- ',iiu30
        do i=1,iuniq
          write(7,'(a5,i4,90f4.1)') '  %d ',i,(flnkd(i,j),j=iiu1,iiu30)
        enddo
          iiu1=iiu1+30
          iiu30=iiu30+30
        enddo
          write(7,'(i5,a5,i5)') iiu1,' --- ',iiu1+numx-1
        do i=1,iuniq
          write(7,'(a5,i4,200f4.1)') '  %d ',
     & i,(flnkd(i,j),j=iiu1,iiu1+numx-1)
        enddo

      endif
      if(allcc.ne.zero.or.facc(1).ne.zero.or.
     & facc(2).ne.zero.or.facc(3).ne.zero)then
        write(7,*) '  '
        write(7,*)' scaled charge '
        do i=1,iuniq
          write(7,'(i4,f10.5)') i,pchin(i)
        enddo
      endif
      write(7,*) '  '
      return
      end
***********************************************************************
      subroutine imcfld
c
c     calculation for intramolecular atomic charge field
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension x(i200,i200),y(i200,i200),z(i200,i200),r3(i200,i200)
      data zero /0.0d0/
      do i=1,iuniq
        cfldx(i)=zero
        cfldy(i)=zero
        cfldz(i)=zero
      enddo
      do i=1,iuniq
        do j=i,iuniq
          if(i.ne.j)then
            x(i,j)=xkpt(i)-xkpt(j)
            x(j,i)=-x(i,j)
            y(i,j)=ykpt(i)-ykpt(j)
            y(j,i)=-y(i,j)
            z(i,j)=zkpt(i)-zkpt(j)
            z(j,i)=-z(i,j)
            r3(i,j)=(dsqrt(x(i,j)**2+y(i,j)**2+z(j,i)**2))**3
            r3(j,i)=r3(i,j)
          else
            x(i,j)=zero
            y(i,j)=zero
            z(i,j)=zero
            r3(i,j)=zero
          endif
        enddo
      enddo
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            cfldx(i)=cfldx(i)+flnkc(i,j)*pchin(j)*x(i,j)/r3(i,j)
            cfldy(i)=cfldy(i)+flnkc(i,j)*pchin(j)*y(i,j)/r3(i,j)
            cfldz(i)=cfldz(i)+flnkc(i,j)*pchin(j)*z(i,j)/r3(i,j)
          endif
        enddo
      enddo
      return
      end
***********************************************************************
      subroutine implst(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
c
c       intramolecular polarization last 
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/espdat/esfact,iestat,irest,ictat
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          odipx(i200),odipy(i200),odipz(i200),dipmom(i200)
      data zero /0.0d0/
c
c     adipx adipy adipz    intramolecular induced dipole moments
c
      if(irest.eq.63)then
        call lstfldms(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
      else if(irest.eq.73)then
        call lstfldli(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
      else if(irest.eq.4)then
        call lstfldt(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
      else if(irest.eq.5)then
        call lstfldp(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
      else if(irest.eq.6)then
        call lstfldm(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
      else if(irest.eq.43)then
        call lstfldts(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
      else if(irest.eq.53)then
        call lstfldps(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
      else if(irest.eq.3)then
        call lstflds(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
      endif
      return
      end
***********************************************************************
      subroutine impcyc(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
c
c       intramolecular polarization cycle
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/espdat/esfact,iestat,irest,ictat
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          odipx(i200),odipy(i200),odipz(i200),dipmom(i200)
      data zero /0.0d0/
c
c     adipx adipy adipz    intramolecular induced dipole moments
c
        if(irest.eq.63)then
          call imdfldms(ll,kk,alpi,xdip,ydip,zdip)
        else if(irest.eq.73)then
          call imdfldli(ll,kk,alpi,xdip,ydip,zdip)
        else if(irest.eq.4)then
          call imdfldt(ll,kk,alpi,xdip,ydip,zdip)
        else if(irest.eq.5)then
          call imdfldp(ll,kk,alpi,xdip,ydip,zdip)
        else if(irest.eq.6)then
          call imdfldm(ll,kk,alpi,xdip,ydip,zdip)
        else if(irest.eq.43)then
          call imdfldts(ll,kk,alpi,xdip,ydip,zdip)
        else if(irest.eq.53)then
          call imdfldps(ll,kk,alpi,xdip,ydip,zdip)
        else if(irest.eq.3)then
          call imdflds(ll,kk,alpi,xdip,ydip,zdip)
        endif
C
      return
      end
***********************************************************************
      subroutine imdfldt(ll,kk,alpi,xdip,ydip,zdip)
c
c     irest=4  Matrix Inversion
c
c     calculations for intramolecular dipole field
c     thole's dipole field tensor scaling
c
C     dumping r3: 1-(a2u2/2+au+1)e(-au)
C             r5: 1-(a3u3/6+a2u2/2+au+1)e(-au)
C
C     a=2.089 A   Thole chem. phys. 59 (1981) 341-350.
C     a=1.9088 A  Duijnen & Swart J. Phys. Chem. A 102 (1998) 2399-2407.
C
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/tens/
     &  txx(i200,i200),txy(i200,i200),
     &  txz(i200,i200),
     &  tyx(i200,i200),tyy(i200,i200),
     &  tyz(i200,i200),
     &  tzx(i200,i200),tzy(i200,i200),
     &  tzz(i200,i200)
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
!*--------------------------------------------------------------------------
      dimension A(i200*3,i200*3),DA(i200*3)
!*--------------------------------------------------------------------------
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/
c
c     thole's exponential scaling parameter is 2.089.
c     van duijinen                             1.9088 
c
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            rx=ar(i,j)*bohr
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            aa=(alpi(i)*alpi(j))
C%SN
C%SN        avoid /zero      2019.3.31
C%SN
            if(aa.gt.0.00001d0)then
              aa=aa**(one/six)
              u=rx/(aa*bohr)
              ar1=factt*u
              ar2=ar1*ar1
              ar3=ar1*ar2
              ex=exp(-ar1)
              s0=ar2/two+ar1+one
              s1=(one-s0*ex)/r3
              s2=three*(one-(ar3/six+s0)*ex)/r5
            else
              s1=one/r3
              s2=three/r5
            endif
C
            txx(i,j)=-(s2*ax(i,j)*ax(i,j)-s1)
            txy(i,j)=-(s2*ax(i,j)*ay(i,j))
            txz(i,j)=-(s2*ax(i,j)*az(i,j))
            tyx(i,j)=-(s2*ay(i,j)*ax(i,j))
            tyy(i,j)=-(s2*ay(i,j)*ay(i,j)-s1)
            tyz(i,j)=-(s2*ay(i,j)*az(i,j))
            tzx(i,j)=-(s2*az(i,j)*ax(i,j))
            tzy(i,j)=-(s2*az(i,j)*ay(i,j))
            tzz(i,j)=-(s2*az(i,j)*az(i,j)-s1)

C
          else
            txx(i,j)=one/alpi(i)
            txy(i,j)=zero
            txz(i,j)=zero
            tyx(i,j)=zero
            tyy(i,j)=one/alpi(i)
            tyz(i,j)=zero
            tzx(i,j)=zero
            tzy(i,j)=zero
            tzz(i,j)=one/alpi(i)
          endif
        enddo
      enddo
!*----------------------------------------------------------
!*    Matrix A
!*----------------------------------------------------------
      k=1
      do i=1,iuniq
        l=1
        do j=1,iuniq
          A(k,l)=txx(i,j)
          A(k,l+1)=txy(i,j)
          A(k,l+2)=txz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+1,l)=tyx(i,j)
          A(k+1,l+1)=tyy(i,j)
          A(k+1,l+2)=tyz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+2,l)=tzx(i,j)
          A(k+2,l+1)=tzy(i,j)
          A(k+2,l+2)=tzz(i,j)
          l=l+3
        enddo
        k=k+3
      enddo
!*
      call houhol0(ll,kk,iuniq,A,DA)
!*  
!*----------------------------------------------------------
      j=0
      do i=1,iuniq*3,3
       j=j+1
       xdip(j)=DA(i)
       ydip(j)=DA(i+1)
       zdip(j)=DA(i+2)
      enddo
!*----------------------------------------------------------
      return
      end
***********************************************************************
      subroutine imdfldp(ll,kk,alpi,xdip,ydip,zdip)
c
c     irest=5    Matrix Inversion
c
c     calculations for intramolecular dipole field
c     thole's dipole field tensor scaling
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/tens/
     &  txx(i200,i200),txy(i200,i200),
     &  txz(i200,i200),
     &  tyx(i200,i200),tyy(i200,i200),
     &  tyz(i200,i200),
     &  tzx(i200,i200),tzy(i200,i200),
     &  tzz(i200,i200)
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
!*--------------------------------------------------------------------------
      dimension A(i200*3,i200*3),DA(i200*3)
!*--------------------------------------------------------------------------
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/,four/4.0d0/
c
c     Thole's expo -au3 scaling      parameter is 0.572
c     Ren & Ponder                                0.39
c
C     r3: 1-e(-au3)
C     r5: 1-(1+au3)e(-au3)
C
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            rx=ar(i,j)*bohr
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            aa=alpi(i)*alpi(j)
            if(aa.gt.0.00001d0)then
              aa=aa**(one/six)
              u=rx/(aa*bohr)
              au3=factt*u*u*u
              ex=exp(-au3)
              s1=(one-ex)/r3
              s2=three*(one-(one+au3)*ex)/r5
            else
              s1=one/r3
              s2=three/r5
            endif
C
            txx(i,j)=-(s2*ax(i,j)*ax(i,j)-s1)
            txy(i,j)=-(s2*ax(i,j)*ay(i,j))
            txz(i,j)=-(s2*ax(i,j)*az(i,j))
            tyx(i,j)=-(s2*ax(i,j)*ay(i,j))
            tyy(i,j)=-(s2*ay(i,j)*ay(i,j)-s1)
            tyz(i,j)=-(s2*ay(i,j)*az(i,j))
            tzx(i,j)=-(s2*ax(i,j)*az(i,j))
            tzy(i,j)=-(s2*ay(i,j)*az(i,j))
            tzz(i,j)=-(s2*az(i,j)*az(i,j)-s1)
C
          else
            txx(i,j)=one/alpi(i)
            txy(i,j)=zero
            txz(i,j)=zero
            tyx(i,j)=zero
            tyy(i,j)=one/alpi(i)
            tyz(i,j)=zero
            tzx(i,j)=zero
            tzy(i,j)=zero
            tzz(i,j)=one/alpi(i)
          endif
        enddo
      enddo
!*----------------------------------------------------------
!*    Matrix A
!*----------------------------------------------------------
      k=1
      do i=1,iuniq
        l=1
        do j=1,iuniq
          A(k,l)=txx(i,j)
          A(k,l+1)=txy(i,j)
          A(k,l+2)=txz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+1,l)=tyx(i,j)
          A(k+1,l+1)=tyy(i,j)
          A(k+1,l+2)=tyz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+2,l)=tzx(i,j)
          A(k+2,l+1)=tzy(i,j)
          A(k+2,l+2)=tzz(i,j)
          l=l+3
        enddo
        k=k+3
      enddo
!*
      call houhol0(ll,kk,iuniq,A,DA)
!*  
!*----------------------------------------------------------
c     calculation for intramolecular induced dipoles
!*----------------------------------------------------------
      j=0
      do i=1,iuniq*3,3
       j=j+1
       xdip(j)=DA(i)
       ydip(j)=DA(i+1)
       zdip(j)=DA(i+2)
      enddo
!*----------------------------------------------------------
      return
      end
***********************************************************************
      subroutine imdfldli(ll,kk,alpi,xdip,ydip,zdip)
c
c     irest=73     Matrix Inversion        
c
c     calculations for intramolecular dipole field
c     thole's dipole field tensor Linear scaling
c
C     dumping r3: 4v**3-3v**4
C             r5: v**4
C
C     a=1.662    Thole chem. phys. 59 (1981) 341-350.
C     a=1.7823   Duijnen & Swart J. Phys. Chem. A 102 (1998) 2399-2407.
C
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/tens/
     &  txx(i200,i200),txy(i200,i200),
     &  txz(i200,i200),
     &  tyx(i200,i200),tyy(i200,i200),
     &  tyz(i200,i200),
     &  tzx(i200,i200),tzy(i200,i200),
     &  tzz(i200,i200)
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
!*--------------------------------------------------------------------------
      dimension A(i200*3,i200*3),DA(i200*3)
!*--------------------------------------------------------------------------
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/,four/4.0d0/
c
c     thole's exponential scaling parameter is 1.662
c     van duijinen                             1.7823 
c
c     do i=1,iuniq
c         xf(i)=pfldx(ll,kk,i)
c         yf(i)=pfldy(ll,kk,i)
c         zf(i)=pfldz(ll,kk,i)
c     enddo
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            rx=ar(i,j)*bohr
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            aa=(alpi(i)*alpi(j))
C%SN
C%SN        avoid /zero      2019.3.31
C%SN
            if(aa.gt.0.0001d0)then
              aa=aa**(one/six)
              s=factt*aa
              v=r1/s
              if(r1.lt.s) then
                s1=(four*(v**3)-three*(v**4))/r3
                s2=three*(v**4)/r5
              else
                s1=one/r3
                s2=three/r5
              endif
            else
              s1=one/r3
              s2=three/r5
            endif
C
            txx(i,j)=-(s2*ax(i,j)*ax(i,j)-s1)*flnkd(i,j)
            txy(i,j)=-(s2*ax(i,j)*ay(i,j))*flnkd(i,j)
            txz(i,j)=-(s2*ax(i,j)*az(i,j))*flnkd(i,j)
            tyx(i,j)=-(s2*ax(i,j)*ay(i,j))*flnkd(i,j)
            tyy(i,j)=-(s2*ay(i,j)*ay(i,j)-s1)*flnkd(i,j)
            tyz(i,j)=-(s2*ay(i,j)*az(i,j))*flnkd(i,j)
            tzx(i,j)=-(s2*ax(i,j)*az(i,j))*flnkd(i,j)
            tzy(i,j)=-(s2*ay(i,j)*az(i,j))*flnkd(i,j)
            tzz(i,j)=-(s2*az(i,j)*az(i,j)-s1)*flnkd(i,j)

C
          else
            txx(i,j)=one/alpi(i)
            txy(i,j)=zero
            txz(i,j)=zero
            tyx(i,j)=zero
            tyy(i,j)=one/alpi(i)
            tyz(i,j)=zero
            tzx(i,j)=zero
            tzy(i,j)=zero
            tzz(i,j)=one/alpi(i)
          endif
        enddo
      enddo
!*----------------------------------------------------------
!*    Matrix A
!*----------------------------------------------------------
      k=1
      do i=1,iuniq
        l=1
        do j=1,iuniq
          A(k,l)=txx(i,j)
          A(k,l+1)=txy(i,j)
          A(k,l+2)=txz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+1,l)=tyx(i,j)
          A(k+1,l+1)=tyy(i,j)
          A(k+1,l+2)=tyz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+2,l)=tzx(i,j)
          A(k+2,l+1)=tzy(i,j)
          A(k+2,l+2)=tzz(i,j)
          l=l+3
        enddo
        k=k+3
      enddo
!*
      call houhol0(ll,kk,iuniq,A,DA)
!*
!*----------------------------------------------------------
c     calculation for intramolecular induced dipoles
!*----------------------------------------------------------
      j=0
      do i=1,iuniq*3,3
       j=j+1
       xdip(j)=DA(i)
       ydip(j)=DA(i+1)
       zdip(j)=DA(i+2)
      enddo
!*----------------------------------------------------------
      return
      end
***********************************************************************
      subroutine imdfldts(ll,kk,alpi,xdip,ydip,zdip)
c
c     irest=43  Matrix Inversion     
c
c     calculations for intramolecular dipole field
c     thole's dipole field tensor scaling
c
C     dumping r3: 1-(a2u2/2+au+1)e(-au)
C             r5: 1-(a3u3/6+a2u2/2+au+1)e(-au)
C
C     a=2.089 A   Thole chem. phys. 59 (1981) 341-350.
C     a=1.9088 A  Duijnen & Swart J. Phys. Chem. A 102 (1998) 2399-2407.
C
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/tens/
     &  txx(i200,i200),txy(i200,i200),
     &  txz(i200,i200),
     &  tyx(i200,i200),tyy(i200,i200),
     &  tyz(i200,i200),
     &  tzx(i200,i200),tzy(i200,i200),
     &  tzz(i200,i200)
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
!*--------------------------------------------------------------------------
      dimension A(i200*3,i200*3),DA(i200*3)
!*--------------------------------------------------------------------------
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/
c
c     thole's exponential scaling parameter is 2.089.
c     van duijinen                             1.9088 
c
c     do i=1,iuniq
c         xf(i)=pfldx(ll,kk,i)
c         yf(i)=pfldy(ll,kk,i)
c         zf(i)=pfldz(ll,kk,i)
c     enddo
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            rx=ar(i,j)*bohr
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            aa=(alpi(i)*alpi(j))
C   
C           avoid /zero      2019.3.31
C   
            if(aa.gt.0.0001d0)then
              aa=aa**(one/six)
              u=rx/(aa*bohr)
              ar1=factt*u
              ar2=ar1*ar1
              ar3=ar1*ar2
              ex=exp(-ar1)
              s0=ar2/two+ar1+one
              s1=(one-s0*ex)/r3
              s2=three*(one-(ar3/six+s0)*ex)/r5
            else
              s1=one/r3
              s2=three/r5
            endif
C
            txx(i,j)=-(s2*ax(i,j)*ax(i,j)-s1)*flnkd(i,j)
            txy(i,j)=-(s2*ax(i,j)*ay(i,j))*flnkd(i,j)
            txz(i,j)=-(s2*ax(i,j)*az(i,j))*flnkd(i,j)
            tyx(i,j)=-(s2*ax(i,j)*ay(i,j))*flnkd(i,j)
            tyy(i,j)=-(s2*ay(i,j)*ay(i,j)-s1)*flnkd(i,j)
            tyz(i,j)=-(s2*ay(i,j)*az(i,j))*flnkd(i,j)
            tzx(i,j)=-(s2*ax(i,j)*az(i,j))*flnkd(i,j)
            tzy(i,j)=-(s2*ay(i,j)*az(i,j))*flnkd(i,j)
            tzz(i,j)=-(s2*az(i,j)*az(i,j)-s1)*flnkd(i,j)
C  
          else
            txx(i,j)=one/alpi(i)
            txy(i,j)=zero
            txz(i,j)=zero
            tyx(i,j)=zero
            tyy(i,j)=one/alpi(i)
            tyz(i,j)=zero
            tzx(i,j)=zero
            tzy(i,j)=zero
            tzz(i,j)=one/alpi(i)
          endif
        enddo
      enddo
!*----------------------------------------------------------
!*    Matrix A
!*----------------------------------------------------------
      k=1
      do i=1,iuniq
        l=1
        do j=1,iuniq
          A(k,l)=txx(i,j)
          A(k,l+1)=txy(i,j)
          A(k,l+2)=txz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+1,l)=tyx(i,j)
          A(k+1,l+1)=tyy(i,j)
          A(k+1,l+2)=tyz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+2,l)=tzx(i,j)
          A(k+2,l+1)=tzy(i,j)
          A(k+2,l+2)=tzz(i,j)
          l=l+3
        enddo
        k=k+3
      enddo
!*
      call houhol0(ll,kk,iuniq,A,DA)
!*
!*----------------------------------------------------------
c     calculation for intramolecular induced dipoles
!*----------------------------------------------------------
      j=0
      do i=1,iuniq*3,3
       j=j+1
       xdip(j)=DA(i)
       ydip(j)=DA(i+1)
       zdip(j)=DA(i+2)
      enddo
!*----------------------------------------------------------
      return
      end
***********************************************************************
      subroutine imdfldps(ll,kk,alpi,xdip,ydip,zdip)
c
c     irest=53  Matrix Inversion   
c
c     calculations for intramolecular dipole field
c     thole's dipole field tensor scaling
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/tens/
     &  txx(i200,i200),txy(i200,i200),
     &  txz(i200,i200),
     &  tyx(i200,i200),tyy(i200,i200),
     &  tyz(i200,i200),
     &  tzx(i200,i200),tzy(i200,i200),
     &  tzz(i200,i200)
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
!*--------------------------------------------------------------------------
      dimension A(i200*3,i200*3),DA(i200*3)
!*--------------------------------------------------------------------------
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/,four/4.0d0/
c
c     Thole's expo -au3 scaling      parameter is 0.572
c     Ren & Ponder                                0.39
c
C     r3: 1-e(-au3)
C     r5: 1-(1+au3)e(-au3)
C
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            rx=ar(i,j)*bohr
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            aa=alpi(i)*alpi(j)
            if(aa.gt.0.0001d0)then
              aa=aa**(one/six)
              u=rx/(aa*bohr)
              au3=factt*u*u*u
              ex=exp(-au3)
              s1=(one-ex)/r3
              s2=three*(one-(one+au3)*ex)/r5
            else
              s1=one/r3
              s2=three/r5
            endif
C
            txx(i,j)=-(s2*ax(i,j)*ax(i,j)-s1)*flnkd(i,j)
            txy(i,j)=-(s2*ax(i,j)*ay(i,j))*flnkd(i,j)
            txz(i,j)=-(s2*ax(i,j)*az(i,j))*flnkd(i,j)
            tyx(i,j)=-(s2*ax(i,j)*ay(i,j))*flnkd(i,j)
            tyy(i,j)=-(s2*ay(i,j)*ay(i,j)-s1)*flnkd(i,j)
            tyz(i,j)=-(s2*ay(i,j)*az(i,j))*flnkd(i,j)
            tzx(i,j)=-(s2*ax(i,j)*az(i,j))*flnkd(i,j)
            tzy(i,j)=-(s2*ay(i,j)*az(i,j))*flnkd(i,j)
            tzz(i,j)=-(s2*az(i,j)*az(i,j)-s1)*flnkd(i,j)

C
          else
            txx(i,j)=one/alpi(i)
            txy(i,j)=zero
            txz(i,j)=zero
            tyx(i,j)=zero
            tyy(i,j)=one/alpi(i)
            tyz(i,j)=zero
            tzx(i,j)=zero
            tzy(i,j)=zero
            tzz(i,j)=one/alpi(i)
          endif
        enddo
      enddo
!*----------------------------------------------------------
!*    Matrix A
!*----------------------------------------------------------
      k=1
      do i=1,iuniq
        l=1
        do j=1,iuniq
          A(k,l)=txx(i,j)
          A(k,l+1)=txy(i,j)
          A(k,l+2)=txz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+1,l)=tyx(i,j)
          A(k+1,l+1)=tyy(i,j)
          A(k+1,l+2)=tyz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+2,l)=tzx(i,j)
          A(k+2,l+1)=tzy(i,j)
          A(k+2,l+2)=tzz(i,j)
          l=l+3
        enddo
        k=k+3
      enddo
!*
      call houhol0(ll,kk,iuniq,A,DA)
!*
!*----------------------------------------------------------
c     calculation for intramolecular induced dipoles
!*----------------------------------------------------------
      j=0
      do i=1,iuniq*3,3
       j=j+1
       xdip(j)=DA(i)
       ydip(j)=DA(i+1)
       zdip(j)=DA(i+2)
      enddo
!*----------------------------------------------------------
*     write(6,'(a27,10f13.6)') '1 xyzdip alpi pfldxyz xyzf ',
*    &  xdip(1),ydip(1),zdip(1),alpi(1),
*    &  pfldx(ll,kk,1),pfldy(ll,kk,1),pfldz(ll,kk,1),
*    &  xf(1),yf(1),zf(1)
      return
      end
***********************************************************************
      subroutine imdfldm(ll,kk,alpi,xdip,ydip,zdip)
c
c     irest=6     Matrix Inversion
c
c     calculations for intramolecular dipole field
C     Masia's model  Gaussian density screening  Atom Type a=0.8456
C     0.8456 from J. Wang st al. JCTC 2019 15 1146
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/tens/
     &  txx(i200,i200),txy(i200,i200),
     &  txz(i200,i200),
     &  tyx(i200,i200),tyy(i200,i200),
     &  tyz(i200,i200),
     &  tzx(i200,i200),tzy(i200,i200),
     &  tzz(i200,i200)
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
!*--------------------------------------------------------------------------
      dimension A(i200*3,i200*3),DA(i200*3)
!*--------------------------------------------------------------------------
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/,pa/3.1415926d0/
C
      co0=two/(three*dsqrt(two*pa))
      othr=one/three
      co1=co0**(-othr)
      co2=two/(dsqrt(pa))
C
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
C 
            bi=factt*co1*(alpi(i)**(-othr))
            bj=factt*co1*(alpi(j)**(-othr))
            bibj=dsqrt(bi**2+bj**2)
            bij=bi*bj/bibj
            sij=bij*r1
            ex=co2*sij*exp(-(sij**2))
            eef=derf(sij)
            epa=one+(two*(sij**2)/three)
            if(bibj.gt.0.0001d0)then
              s1=(eef-ex)/r3
              s2=three*(eef-ex*epa)/r5
            else
              s1=one/r3
              s2=three/r5
            endif
C
            txx(i,j)=-(s2*ax(i,j)*ax(i,j)-s1)
            txy(i,j)=-(s2*ax(i,j)*ay(i,j))
            txz(i,j)=-(s2*ax(i,j)*az(i,j))
            tyx(i,j)=-(s2*ax(i,j)*ay(i,j))
            tyy(i,j)=-(s2*ay(i,j)*ay(i,j)-s1)
            tyz(i,j)=-(s2*ay(i,j)*az(i,j))
            tzx(i,j)=-(s2*ax(i,j)*az(i,j))
            tzy(i,j)=-(s2*ay(i,j)*az(i,j))
            tzz(i,j)=-(s2*az(i,j)*az(i,j)-s1)
C
          else
            txx(i,j)=one/alpi(i)
            txy(i,j)=zero
            txz(i,j)=zero
            tyx(i,j)=zero
            tyy(i,j)=one/alpi(i)
            tyz(i,j)=zero
            tzx(i,j)=zero
            tzy(i,j)=zero
            tzz(i,j)=one/alpi(i)
          endif
        enddo
      enddo
!*----------------------------------------------------------
!*    Matrix A
!*----------------------------------------------------------
      k=1
      do i=1,iuniq
        l=1
        do j=1,iuniq
          A(k,l)=txx(i,j)
          A(k,l+1)=txy(i,j)
          A(k,l+2)=txz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+1,l)=tyx(i,j)
          A(k+1,l+1)=tyy(i,j)
          A(k+1,l+2)=tyz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+2,l)=tzx(i,j)
          A(k+2,l+1)=tzy(i,j)
          A(k+2,l+2)=tzz(i,j)
          l=l+3
        enddo
        k=k+3
      enddo
!*
      call houhol0(ll,kk,iuniq,A,DA)
!*
!*----------------------------------------------------------
c     calculation for intramolecular induced dipoles
!*----------------------------------------------------------
      j=0
      do i=1,iuniq*3,3
       j=j+1
       xdip(j)=DA(i)
       ydip(j)=DA(i+1)
       zdip(j)=DA(i+2)
      enddo
!*----------------------------------------------------------
      return
      end
***********************************************************************
      subroutine imdfldms(ll,kk,alpi,xdip,ydip,zdip)
c
c     irest=63   Matrix Inversion  
c
c     calculations for intramolecular dipole field
C     Masia's model  Gaussian density screening  Atom Type a=0.8456
C     0.8456 from J. Wang st al. JCTC 2019 15 1146
c     with simple scaling 
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/tens/
     &  txx(i200,i200),txy(i200,i200),
     &  txz(i200,i200),
     &  tyx(i200,i200),tyy(i200,i200),
     &  tyz(i200,i200),
     &  tzx(i200,i200),tzy(i200,i200),
     &  tzz(i200,i200)
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
!*--------------------------------------------------------------------------
      dimension A(i200*3,i200*3),DA(i200*3)
!*--------------------------------------------------------------------------
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/,pa/3.1415926d0/
C
      co0=two/(three*dsqrt(two*pa))
      othr=one/three
      co1=co0**(-othr)
      co2=two/(dsqrt(pa))
C
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
C 
            bi=factt*co1*(alpi(i)**(-othr))
            bj=factt*co1*(alpi(j)**(-othr))
            bibj=dsqrt(bi**2+bj**2)
            bij=bi*bj/bibj
            sij=bij*r1
            ex=co2*sij*exp(-(sij**2))
            eef=derf(sij)
            epa=one+(two*(sij**2)/three)
            if(bibj.gt.0.0001d0)then
              s1=(eef-ex)/r3
              s2=three*(eef-ex*epa)/r5
            else
              s1=one/r3
              s2=three/r5
            endif
C
            txx(i,j)=-(s2*ax(i,j)*ax(i,j)-s1)*flnkd(i,j)
            txy(i,j)=-(s2*ax(i,j)*ay(i,j))*flnkd(i,j)
            txz(i,j)=-(s2*ax(i,j)*az(i,j))*flnkd(i,j)
            tyx(i,j)=-(s2*ax(i,j)*ay(i,j))*flnkd(i,j)
            tyy(i,j)=-(s2*ay(i,j)*ay(i,j)-s1)*flnkd(i,j)
            tyz(i,j)=-(s2*ay(i,j)*az(i,j))*flnkd(i,j)
            tzx(i,j)=-(s2*ax(i,j)*az(i,j))*flnkd(i,j)
            tzy(i,j)=-(s2*ay(i,j)*az(i,j))*flnkd(i,j)
            tzz(i,j)=-(s2*az(i,j)*az(i,j)-s1)*flnkd(i,j)
C
          else
            txx(i,j)=one/alpi(i)
            txy(i,j)=zero
            txz(i,j)=zero
            tyx(i,j)=zero
            tyy(i,j)=one/alpi(i)
            tyz(i,j)=zero
            tzx(i,j)=zero
            tzy(i,j)=zero
            tzz(i,j)=one/alpi(i)

          endif
        enddo
      enddo
!*----------------------------------------------------------
!*    Matrix A
!*----------------------------------------------------------
      k=1
      do i=1,iuniq
        l=1
        do j=1,iuniq
          A(k,l)=txx(i,j)
          A(k,l+1)=txy(i,j)
          A(k,l+2)=txz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+1,l)=tyx(i,j)
          A(k+1,l+1)=tyy(i,j)
          A(k+1,l+2)=tyz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+2,l)=tzx(i,j)
          A(k+2,l+1)=tzy(i,j)
          A(k+2,l+2)=tzz(i,j)
          l=l+3
        enddo
        k=k+3
      enddo
!*
      call houhol0(ll,kk,iuniq,A,DA)
!*
!*----------------------------------------------------------
c     calculation for intramolecular induced dipoles
!*----------------------------------------------------------
      j=0
      do i=1,iuniq*3,3
       j=j+1
       xdip(j)=DA(i)
       ydip(j)=DA(i+1)
       zdip(j)=DA(i+2)
      enddo
!*----------------------------------------------------------
      return
      end
***********************************************************************
      subroutine indchg(adipx,adipy,adipz,xp,yp,zp,xm,ym,zm,q)
c
c     Calculations for induced charges and positions    ictat=0
c     induced dipole length is 1.0  
c     The center of the induced dipole is placed on the atom. (- 0.5 atm 0.5 +)
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      dimension adipx(i200),adipy(i200),adipz(i200),
     & xp(i200),yp(i200),zp(i200),xm(i200),ym(i200),zm(i200),
     & q(i200)
      data d05/0.5d0/,zero/0.0d0/
        do i=1,iuniq
        q(i)=dsqrt(adipx(i)**2+adipy(i)**2+adipz(i)**2)
        if(q(i).eq.zero)then
          xh=zero
          yh=zero
          zh=zero
        else
          xh=d05*adipx(i)/q(i)
          yh=d05*adipy(i)/q(i)
          zh=d05*adipz(i)/q(i)
        endif
        xp(i)=xkpt(i)+xh
        yp(i)=ykpt(i)+yh
        zp(i)=zkpt(i)+zh
        xm(i)=xkpt(i)-xh
        ym(i)=ykpt(i)-yh
        zm(i)=zkpt(i)-zh
      enddo
      return
      end
***********************************************************************
      subroutine indchga(adipx,adipy,adipz,xp,yp,zp,xm,ym,zm,q)
c
c     Calculations for induced charges and positions   ictat=1
c     induced dipole length is 1.0
c     The center of the induced dipole is placed on the atom.
C     A + delta charge is placed on the atom.
C     A - delta charge is located 1 Bohr from the atom.
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      dimension adipx(i200),adipy(i200),adipz(i200),
     & xp(i200),yp(i200),zp(i200),xm(i200),ym(i200),zm(i200),
     & q(i200)
      data d05/0.5d0/,zero/0.0d0/
        do i=1,iuniq
        q(i)=dsqrt(adipx(i)**2+adipy(i)**2+adipz(i)**2)
        if(q(i).eq.zero)then
          xh=zero
          yh=zero
          zh=zero
        else
          xh=adipx(i)/q(i)
          yh=adipy(i)/q(i)
          zh=adipz(i)/q(i)
        endif
        xp(i)=xkpt(i)+xh
        yp(i)=ykpt(i)+yh
        zp(i)=zkpt(i)+zh
        xm(i)=xkpt(i)
        ym(i)=ykpt(i)
        zm(i)=zkpt(i)
      enddo
      return
      end
***********************************************************************
      subroutine intrapop(jfirs,f)
c
c     In addition to a field test charge, intramolecular polarization 
c     from esp charges & induced dipoles is included.
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/espdat/esfact,iestat,irest,ictat
      common/comp12/cxnp(i340),cynp(i340),cznp(i340),
     & fpcx(i340),fpcy(i340),fpcz(i340),
     & qupot(i340,i340,i20)
c
      common/pdpol1/ibon(i200),jbon(i200),ieqbon(i200),ibpair,iptyp,
     & nvarp,ivarp(i200),iestah,iatid,
     & iatm(i200),
     & nvara,nvart,nvaran,ishl,nvarh,
     & ifit,ifld,jfld,isym2(i340)
      common/pdpol2/
     & fpchrg(i20),pkpt(i200*2,4),
     & eqatmi(i200),varai(i200),
     & alptrn(i200),alpani(i200),alpatm(i200)
c
      common/prpol/ xmid(i200),ymid(i200),zmid(i200),
     &          trnx1(i340,i200),trny1(i340,i200),trnz1(i340,i200),
     &          trnx2(i340,i200),trny2(i340,i200),trnz2(i340,i200),
     &          bonx(i200),bony(i200),bonz(i200),bonlen(i200),
     &          fldlen(i340,i200),anifac(i340,i200),ealen(i340,i200),
     &          falen(i340,i200),
     &          signa(i340,i200),ftmp(i340),
     &          atmx1(i340,i200),atmy1(i340,i200),atmz1(i340,i200),
     &          atmx2(i340,i200),atmy2(i340,i200),atmz2(i340,i200) 
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/odipl/oxdip(i20,i340,i200),oydip(i20,i340,i200),
     &  ozdip(i20,i340,i200)
      dimension var(1),f(1)
      dimension dltaq(i200)
      dimension dx(i200),dy(i200),dz(i200),xp(i200),yp(i200),zp(i200),
     &  xm(i200),ym(i200),zm(i200),q(i200),alpi(i200),
     &  ddx(i200),ddy(i200),ddz(i200)
      data aukcal /627.509541d0/
c
      kk=0
      do l=1,ishl
        do k=1,ifld
          do j=1,iatid
            jj=jfirs+j
            alpi(j)=pkpt(jj,1)
            dx(j)=alpi(j)*pfldx(l,k,j)
            dy(j)=alpi(j)*pfldy(l,k,j)
            dz(j)=alpi(j)*pfldz(l,k,j)
          enddo
          call impcyc(l,k,dx,dy,dz,alpi,ddx,ddy,ddz)
C
C       dipole moment --> common odipl --> subr potout
C
          do  j=1,iatid
            oxdip(l,k,j)=ddx(j)
            oydip(l,k,j)=ddy(j)
            ozdip(l,k,j)=ddz(j)
          enddo
C   
          if(ictat.eq.1)then
             call indchga(ddx,ddy,ddz,xp,yp,zp,xm,ym,zm,q)
          else
             call indchg(ddx,ddy,ddz,xp,yp,zp,xm,ym,zm,q)
          endif
C
          do i=1,iestah
            kk=kk+1
            do j=1,iatid
              xa     =cxnp(i)-xp(j)
              ya     =cynp(i)-yp(j)
              za     =cznp(i)-zp(j)
              ra     =dsqrt(xa*xa+ya*ya+za*za)
              xb     =cxnp(i)-xm(j)
              yb     =cynp(i)-ym(j)
              zb     =cznp(i)-zm(j)
              rb     =dsqrt(xb*xb+yb*yb+zb*zb)
              f(kk)=f(kk)-q(j)/ra+q(j)/rb
            enddo
          enddo
        enddo
      enddo
      return
      end
***********************************************************************
      subroutine intralst(jfirs,f)
c
c     in addition to a field test charge, intramolecular polarization 
c     from esp charges & induced dipoles is included.
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/espdat/esfact,iestat,irest,ictat
      common/comp12/cxnp(i340),cynp(i340),cznp(i340),
     & fpcx(i340),fpcy(i340),fpcz(i340),
     & qupot(i340,i340,i20)
c
      common/pdpol1/ibon(i200),jbon(i200),ieqbon(i200),ibpair,iptyp,
     & nvarp,ivarp(i200),iestah,iatid,
     & iatm(i200),
     & nvara,nvart,nvaran,ishl,nvarh,
     & ifit,ifld,jfld,isym2(i340)
      common/pdpol2/
     & fpchrg(i20),pkpt(i200*2,4),
     & eqatmi(i200),varai(i200),
     & alptrn(i200),alpani(i200),alpatm(i200)
c
      common/prpol/ xmid(i200),ymid(i200),zmid(i200),
     &          trnx1(i340,i200),trny1(i340,i200),trnz1(i340,i200),
     &          trnx2(i340,i200),trny2(i340,i200),trnz2(i340,i200),
     &          bonx(i200),bony(i200),bonz(i200),bonlen(i200),
     &          fldlen(i340,i200),anifac(i340,i200),ealen(i340,i200),
     &          falen(i340,i200),
     &          signa(i340,i200),ftmp(i340),
     &          atmx1(i340,i200),atmy1(i340,i200),atmz1(i340,i200),
     &          atmx2(i340,i200),atmy2(i340,i200),atmz2(i340,i200) 
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/enedip/ upol(i340,i20),
     &  cdpx(i340,i20),cdpy(i340,i20),cdpz(i340,i20)
      dimension dltaq(i200)
      dimension dx(i200),dy(i200),dz(i200),xp(i200),yp(i200),zp(i200),
     &  xm(i200),ym(i200),zm(i200),q(i200),alpi(i200),
     &  ddx(i200),ddy(i200),ddz(i200)
      data aukcal /627.509541d0/
c
      kk=0
      do l=1,ishl
        do k=1,ifld
          do j=1,iatid
            jj=jfirs+j
            fldx=pfldx(l,k,j)
            fldy=pfldy(l,k,j)
            fldz=pfldz(l,k,j)
            alpi(j)=pkpt(jj,1)
            dx(j)=alpi(j)*fldx
            dy(j)=alpi(j)*fldy
            dz(j)=alpi(j)*fldz
          enddo
          call implst(l,k,dx,dy,dz,alpi,ddx,ddy,ddz)
          do j=1,iatid
            cdpx(k,l)=cdpx(k,l)-ddx(j)          
            cdpy(k,l)=cdpy(k,l)-ddy(j)          
            cdpz(k,l)=cdpz(k,l)-ddz(j)          
            ene=ddx(j)*fldx0(l,k,j)+ddy(j)*fldy0(l,k,j)
     &       +ddz(j)*fldz0(l,k,j)
            upol(k,l)=upol(k,l)-0.5d0*ene
          enddo
        enddo
      enddo
      return
      end
***********************************************************************
      subroutine lstflds(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
c
c     irest=3  last
c
c     calculations for intramolecular dipole field
c     simple scaling model
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/tens/
     &  txx(i200,i200),txy(i200,i200),
     &  txz(i200,i200),
     &  tyx(i200,i200),tyy(i200,i200),
     &  tyz(i200,i200),
     &  tzx(i200,i200),tzy(i200,i200),
     &  tzz(i200,i200)
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/
c
      do i=1,iuniq
          xf(i)=pfldx(ll,kk,i)
          yf(i)=pfldy(ll,kk,i)
          zf(i)=pfldz(ll,kk,i)
      enddo
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            s1=one/r3
            s2=three/r5
C
            txx(i,j)=-(s2*ax(i,j)*ax(i,j)-s1)*flnkd(i,j)
            txy(i,j)=-(s2*ax(i,j)*ay(i,j))*flnkd(i,j)
            txz(i,j)=-(s2*ax(i,j)*az(i,j))*flnkd(i,j)
            tyx(i,j)=-(s2*ax(i,j)*ay(i,j))*flnkd(i,j)
            tyy(i,j)=-(s2*ay(i,j)*ay(i,j)-s1)*flnkd(i,j)
            tyz(i,j)=-(s2*ay(i,j)*az(i,j))*flnkd(i,j)
            tzx(i,j)=-(s2*ax(i,j)*az(i,j))*flnkd(i,j)
            tzy(i,j)=-(s2*ay(i,j)*az(i,j))*flnkd(i,j)
            tzz(i,j)=-(s2*az(i,j)*az(i,j)-s1)*flnkd(i,j)
            
C
            xf(i)=xf(i)-(txx(i,j)*adipx(j)+
     &                   txy(i,j)*adipy(j)+
     &                   txz(i,j)*adipz(j))
            yf(i)=yf(i)-(tyx(i,j)*adipx(j)+
     &                   tyy(i,j)*adipy(j)+
     &                   tyz(i,j)*adipz(j))
            zf(i)=zf(i)-(tzx(i,j)*adipx(j)+
     &                   tzy(i,j)*adipy(j)+
     &                   tzz(i,j)*adipz(j))
          else
            txx(i,j)=one/alpi(i)
            txy(i,j)=zero
            txz(i,j)=zero
            tyx(i,j)=zero
            tyy(i,j)=one/alpi(i)
            tyz(i,j)=zero
            tzx(i,j)=zero
            tzy(i,j)=zero
            tzz(i,j)=one/alpi(i)
          endif
        enddo
      enddo
c
c     calculation for intramolecular induced dipoles
c     
      do i=1,iuniq
       xdip(i)=alpi(i)*xf(i)
       ydip(i)=alpi(i)*yf(i)
       zdip(i)=alpi(i)*zf(i)
      enddo
      return
      end
***********************************************************************
      subroutine lstfldt(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
c
c     irest=4
c
c     calculations for intramolecular dipole field
c     thole's dipole field tensor scaling
c
C     dumping r3: 1-(a2u2/2+au+1)e(-au)
C             r5: 1-(a3u3/6+a2u2/2+au+1)e(-au)
C
C     a=2.089 A   Thole chem. phys. 59 (1981) 341-350.
C     a=1.9088 A  Duijnen & Swart J. Phys. Chem. A 102 (1998) 2399-2407.
C
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/
c
c     thole's exponential scaling parameter is 2.089.
c     van duijinen                             1.9088 
c
      do i=1,iuniq
          xf(i)=pfldx(ll,kk,i)
          yf(i)=pfldy(ll,kk,i)
          zf(i)=pfldz(ll,kk,i)
      enddo
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            rx=ar(i,j)*bohr
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            aa=(alpi(i)*alpi(j))
C%SN
C%SN        avoid /zero      2019.3.31
C%SN
            if(aa.gt.0.00001d0)then
              aa=aa**(one/six)
              u=rx/(aa*bohr)
              ar1=factt*u
              ar2=ar1*ar1
              ar3=ar1*ar2
              ex=exp(-ar1)
              s0=ar2/two+ar1+one
              s1=(one-s0*ex)/r3
              s2=three*(one-(ar3/six+s0)*ex)/r5
            else
              s1=one/r3
              s2=three/r5
            endif
C%SN
            xf(i)=xf(i)+((s2*ax(i,j)*ax(i,j)-s1)*adipx(j)+
     &                    s2*ax(i,j)*ay(i,j)*adipy(j)+
     &                    s2*ax(i,j)*az(i,j)*adipz(j))
            yf(i)=yf(i)+ (s2*ax(i,j)*ay(i,j)*adipx(j)+
     &                   (s2*ay(i,j)*ay(i,j)-s1)*adipy(j)+
     &                    s2*ay(i,j)*az(i,j)*adipz(j))
            zf(i)=zf(i)+ (s2*ax(i,j)*az(i,j)*adipx(j)+
     &                    s2*ay(i,j)*az(i,j)*adipy(j)+
     &                   (s2*az(i,j)*az(i,j)-s1)*adipz(j))
          endif
        enddo
      enddo
c
c     calculation for intramolecular induced dipoles
c
      do i=1,iuniq
       xdip(i)=alpi(i)*xf(i)
       ydip(i)=alpi(i)*yf(i)
       zdip(i)=alpi(i)*zf(i)
      enddo
      return
      end
***********************************************************************
      subroutine lastfldp(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
c
c     irest=5
c
c     calculations for intramolecular dipole field
c     thole's dipole field tensor scaling
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/,four/4.0d0/
c
c     Thole's expo -au3 scaling      parameter is 0.572
c     Ren & Ponder                                0.39
c
C     r3: 1-e(-au3)
C     r5: 1-(1+au3)e(-au3)
C
      do i=1,iuniq
          xf(i)=pfldx(ll,kk,i)
          yf(i)=pfldy(ll,kk,i)
          zf(i)=pfldz(ll,kk,i)
      enddo
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            rx=ar(i,j)*bohr
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            aa=alpi(i)*alpi(j)
            if(aa.gt.0.00001d0)then
              aa=aa**(one/six)
              u=rx/(aa*bohr)
              au3=factt*u*u*u
              ex=exp(-au3)
              s1=(one-ex)/r3
              s2=three*(one-(one+au3)*ex)/r5
            else
              s1=one/r3
              s2=three/r5
            endif
            xf(i)=xf(i)+((s2*ax(i,j)*ax(i,j)-s1)*adipx(j)+
     &                    s2*ax(i,j)*ay(i,j)*adipy(j)+
     &                    s2*ax(i,j)*az(i,j)*adipz(j))
            yf(i)=yf(i)+ (s2*ax(i,j)*ay(i,j)*adipx(j)+
     &                   (s2*ay(i,j)*ay(i,j)-s1)*adipy(j)+
     &                    s2*ay(i,j)*az(i,j)*adipz(j))
            zf(i)=zf(i)+ (s2*ax(i,j)*az(i,j)*adipx(j)+
     &                    s2*ay(i,j)*az(i,j)*adipy(j)+
     &                   (s2*az(i,j)*az(i,j)-s1)*adipz(j))
          endif
        enddo
      enddo
c
c     calculation for intramolecular induced dipoles
c
      do i=1,iuniq
       xdip(i)=alpi(i)*xf(i)
       ydip(i)=alpi(i)*yf(i)
       zdip(i)=alpi(i)*zf(i)
      enddo
      return
      end
***********************************************************************
      subroutine lstfldli(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
c
c     irest=73          2020.12.18  
c
c     calculations for intramolecular dipole field
c     thole's dipole field tensor Linear scaling
c
C     dumping r3: 4v**3-3v**4
C             r5: v**4
C
C     a=1.662    Thole chem. phys. 59 (1981) 341-350.
C     a=1.7823   Duijnen & Swart J. Phys. Chem. A 102 (1998) 2399-2407.
C
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/,four/4.0d0/
c
c     thole's exponential scaling parameter is 1.662
c     van duijinen                             1.7823 
c
      do i=1,iuniq
          xf(i)=pfldx(ll,kk,i)
          yf(i)=pfldy(ll,kk,i)
          zf(i)=pfldz(ll,kk,i)
      enddo
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            rx=ar(i,j)*bohr
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            aa=(alpi(i)*alpi(j))
C%SN
C%SN        avoid /zero      2019.3.31
C%SN
            if(aa.gt.0.0001d0)then
              aa=aa**(one/six)
              s=factt*aa
              v=r1/s
              if(r1.lt.s) then
                s1=(four*(v**3)-three*(v**4))/r3
                s2=three*(v**4)/r5
              else
                s1=one/r3
                s2=three/r5
              endif
            else
              s1=one/r3
              s2=three/r5
            endif
C%SN
            xf(i)=xf(i)+((s2*ax(i,j)*ax(i,j)-s1)*adipx(j)+
     &                  s2*ax(i,j)*ay(i,j)*adipy(j)+
     &                  s2*ax(i,j)*az(i,j)*adipz(j))*flnkd(i,j)
            yf(i)=yf(i)+ (s2*ax(i,j)*ay(i,j)*adipx(j)+
     &                  (s2*ay(i,j)*ay(i,j)-s1)*adipy(j)+
     &                  s2*ay(i,j)*az(i,j)*adipz(j))*flnkd(i,j)
            zf(i)=zf(i)+ (s2*ax(i,j)*az(i,j)*adipx(j)+
     &                  s2*ay(i,j)*az(i,j)*adipy(j)+
     &                 (s2*az(i,j)*az(i,j)-s1)*adipz(j))*flnkd(i,j)
          endif
        enddo
      enddo
c
c     calculation for intramolecular induced dipoles
c
      do i=1,iuniq
       xdip(i)=alpi(i)*xf(i)
       ydip(i)=alpi(i)*yf(i)
       zdip(i)=alpi(i)*zf(i)
      enddo
      return
      end
***********************************************************************
      subroutine lstfldts(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
c
c     irest=43          2019.8.15  9.21 dabs
c
c     calculations for intramolecular dipole field
c     thole's dipole field tensor scaling
c
C     dumping r3: 1-(a2u2/2+au+1)e(-au)
C             r5: 1-(a3u3/6+a2u2/2+au+1)e(-au)
C
C     a=2.089 A   Thole chem. phys. 59 (1981) 341-350.
C     a=1.9088 A  Duijnen & Swart J. Phys. Chem. A 102 (1998) 2399-2407.
C
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/
c
c     thole's exponential scaling parameter is 2.089.
c     van duijinen                             1.9088 
c
      do i=1,iuniq
          xf(i)=pfldx(ll,kk,i)
          yf(i)=pfldy(ll,kk,i)
          zf(i)=pfldz(ll,kk,i)
      enddo
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            rx=ar(i,j)*bohr
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            aa=(alpi(i)*alpi(j))
C%SN
C%SN        avoid /zero      2019.3.31
C%SN
            if(aa.gt.0.0001d0)then
              aa=aa**(one/six)
              u=rx/(aa*bohr)
              ar1=factt*u
              ar2=ar1*ar1
              ar3=ar1*ar2
              ex=exp(-ar1)
              s0=ar2/two+ar1+one
              s1=(one-s0*ex)/r3
              s2=three*(one-(ar3/six+s0)*ex)/r5
            else
              s1=one/r3
              s2=three/r5
            endif
C%SN
            xf(i)=xf(i)+((s2*ax(i,j)*ax(i,j)-s1)*adipx(j)+
     &                  s2*ax(i,j)*ay(i,j)*adipy(j)+
     &                  s2*ax(i,j)*az(i,j)*adipz(j))*flnkd(i,j)
            yf(i)=yf(i)+ (s2*ax(i,j)*ay(i,j)*adipx(j)+
     &                  (s2*ay(i,j)*ay(i,j)-s1)*adipy(j)+
     &                  s2*ay(i,j)*az(i,j)*adipz(j))*flnkd(i,j)
            zf(i)=zf(i)+ (s2*ax(i,j)*az(i,j)*adipx(j)+
     &                  s2*ay(i,j)*az(i,j)*adipy(j)+
     &                 (s2*az(i,j)*az(i,j)-s1)*adipz(j))*flnkd(i,j)
          endif
        enddo
      enddo
c
c     calculation for intramolecular induced dipoles
c
      do i=1,iuniq
       xdip(i)=alpi(i)*xf(i)
       ydip(i)=alpi(i)*yf(i)
       zdip(i)=alpi(i)*zf(i)
      enddo
      return
      end
***********************************************************************
      subroutine lstfldps(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
c
c     irest=53               2019.8.15  9.21 dabs
c
c     calculations for intramolecular dipole field
c     thole's dipole field tensor scaling
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/,four/4.0d0/
c
c     Thole's expo -au3 scaling      parameter is 0.572
c     Ren & Ponder                                0.39
c
C     r3: 1-e(-au3)
C     r5: 1-(1+au3)e(-au3)
C
      do i=1,iuniq
          xf(i)=pfldx(ll,kk,i)
          yf(i)=pfldy(ll,kk,i)
          zf(i)=pfldz(ll,kk,i)
      enddo
*     write(6,'(16f10.5)') (alpi(i),i=1,16)
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            rx=ar(i,j)*bohr
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            aa=alpi(i)*alpi(j)
            if(aa.gt.0.0001d0)then
              aa=aa**(one/six)
              u=rx/(aa*bohr)
              au3=factt*u*u*u
              ex=exp(-au3)
              s1=(one-ex)/r3
              s2=three*(one-(one+au3)*ex)/r5
            else
              s1=one/r3
              s2=three/r5
            endif
            xf(i)=xf(i)+((s2*ax(i,j)*ax(i,j)-s1)*adipx(j)+
     &                    s2*ax(i,j)*ay(i,j)*adipy(j)+
     &                  s2*ax(i,j)*az(i,j)*adipz(j))*flnkd(i,j)
            yf(i)=yf(i)+ (s2*ax(i,j)*ay(i,j)*adipx(j)+
     &                   (s2*ay(i,j)*ay(i,j)-s1)*adipy(j)+
     &                  s2*ay(i,j)*az(i,j)*adipz(j))*flnkd(i,j)
            zf(i)=zf(i)+ (s2*ax(i,j)*az(i,j)*adipx(j)+
     &                  s2*ay(i,j)*az(i,j)*adipy(j)+
     &               (s2*az(i,j)*az(i,j)-s1)*adipz(j))*flnkd(i,j)
          endif
        enddo
      enddo
c
c     calculation for intramolecular induced dipoles
c
      do i=1,iuniq
       xdip(i)=alpi(i)*xf(i)
       ydip(i)=alpi(i)*yf(i)
       zdip(i)=alpi(i)*zf(i)
      enddo
      return
      end
***********************************************************************
      subroutine lstfldm(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
c
c     irest=6      2019.9.17 ---> 2020.8.4   Sij=Betajj*rij
c
c     calculations for intramolecular dipole field
C     Masia's model  Gaussian density screening  Atom Type a=0.8456
C     0.8456 from J. Wang st al. JCTC 2019 15 1146
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/,pa/3.1415926d0/
c
      do i=1,iuniq
          xf(i)=pfldx(ll,kk,i)
          yf(i)=pfldy(ll,kk,i)
          zf(i)=pfldz(ll,kk,i)
      enddo
C
      co0=two/(three*dsqrt(two*pa))
      othr=one/three
      co1=co0**(-othr)
      co2=two/(dsqrt(pa))
c     osix=one/six
C
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
C 
            bi=factt*co1*(alpi(i)**(-othr))
            bj=factt*co1*(alpi(j)**(-othr))
            bij=bi*bj/(dsqrt(bi**2+bj**2))
            sij=bij*r1
            ex=co2*sij*exp(-(sij**2))
            eef=derf(sij)
            epa=one+(two*(sij**2)/three)
            s1=(eef-ex)/r3
            s2=three*(eef-ex*epa)/r5
C
            xf(i)=xf(i)+((s2*ax(i,j)*ax(i,j)-s1)*adipx(j)+
     &                    s2*ax(i,j)*ay(i,j)*adipy(j)+
     &                    s2*ax(i,j)*az(i,j)*adipz(j))
            yf(i)=yf(i)+ (s2*ax(i,j)*ay(i,j)*adipx(j)+
     &                   (s2*ay(i,j)*ay(i,j)-s1)*adipy(j)+
     &                    s2*ay(i,j)*az(i,j)*adipz(j))
            zf(i)=zf(i)+ (s2*ax(i,j)*az(i,j)*adipx(j)+
     &                    s2*ay(i,j)*az(i,j)*adipy(j)+
     &                   (s2*az(i,j)*az(i,j)-s1)*adipz(j))
          endif
        enddo
      enddo
c
c     calculation for intramolecular induced dipoles
c
      do i=1,iuniq
       xdip(i)=alpi(i)*xf(i)
       ydip(i)=alpi(i)*yf(i)
       zdip(i)=alpi(i)*zf(i)
      enddo
      return
      end
***********************************************************************
      subroutine lstfldms(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
c
c     irest=63      2019.9.2
C                   2020.8.4 
c
c     calculations for intramolecular dipole field
C     Masia's model  Gaussian density screening  Atom Type a=0.8456
C     0.8456 from J. Wang st al. JCTC 2019 15 1146
c     with simple scaling 
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/,pa/3.1415926d0/
c
      do i=1,iuniq
          xf(i)=pfldx(ll,kk,i)
          yf(i)=pfldy(ll,kk,i)
          zf(i)=pfldz(ll,kk,i)
      enddo
C
      co0=two/(three*dsqrt(two*pa))
      othr=one/three
      co1=co0**(-othr)
      co2=two/(dsqrt(pa))
C
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
C 
            bi=factt*co1*(alpi(i)**(-othr))
            bj=factt*co1*(alpi(j)**(-othr))
            bij=bi*bj/(dsqrt(bi**2+bj**2))
            sij=bij*r1
            ex=co2*sij*exp(-(sij**2))
            eef=derf(sij)
            epa=one+(two*(sij**2)/three)
            s1=(eef-ex)/r3
            s2=three*(eef-ex*epa)/r5
C
            xf(i)=xf(i)+((s2*ax(i,j)*ax(i,j)-s1)*adipx(j)+
     &                    s2*ax(i,j)*ay(i,j)*adipy(j)+
     &                    s2*ax(i,j)*az(i,j)*adipz(j))*flnkd(i,j)
            yf(i)=yf(i)+ (s2*ax(i,j)*ay(i,j)*adipx(j)+
     &                   (s2*ay(i,j)*ay(i,j)-s1)*adipy(j)+
     &                    s2*ay(i,j)*az(i,j)*adipz(j))*flnkd(i,j)
            zf(i)=zf(i)+ (s2*ax(i,j)*az(i,j)*adipx(j)+
     &                    s2*ay(i,j)*az(i,j)*adipy(j)+
     &                   (s2*az(i,j)*az(i,j)-s1)*adipz(j))*flnkd(i,j)
          endif
        enddo
      enddo
c
c     calculation for intramolecular induced dipoles
c
      do i=1,iuniq
       xdip(i)=alpi(i)*xf(i)
       ydip(i)=alpi(i)*yf(i)
       zdip(i)=alpi(i)*zf(i)
      enddo
      return
      end
***********************************************************************
      subroutine lstfldp(ll,kk,adipx,adipy,adipz,alpi,xdip,ydip,zdip)
c
c     irest=5
c
c     calculations for intramolecular dipole field
c     thole's dipole field tensor scaling
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/,bohr/0.52917706d0/,four/4.0d0/
c
c     Thole's expo -au3 scaling      parameter is 0.572
c     Ren & Ponder                                0.39
c
C     r3: 1-e(-au3)
C     r5: 1-(1+au3)e(-au3)
C
      do i=1,iuniq
          xf(i)=pfldx(ll,kk,i)
          yf(i)=pfldy(ll,kk,i)
          zf(i)=pfldz(ll,kk,i)
      enddo
*     write(6,'(16f10.5)') (alpi(i),i=1,16)
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            rx=ar(i,j)*bohr
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            aa=alpi(i)*alpi(j)
            if(aa.gt.0.00001d0)then
              aa=aa**(one/six)
              u=rx/(aa*bohr)
              au3=factt*u*u*u
              ex=exp(-au3)
              s1=(one-ex)/r3
              s2=three*(one-(one+au3)*ex)/r5
            else
              s1=one/r3
              s2=three/r5
            endif
            xf(i)=xf(i)+((s2*ax(i,j)*ax(i,j)-s1)*adipx(j)+
     &                    s2*ax(i,j)*ay(i,j)*adipy(j)+
     &                    s2*ax(i,j)*az(i,j)*adipz(j))
            yf(i)=yf(i)+ (s2*ax(i,j)*ay(i,j)*adipx(j)+
     &                   (s2*ay(i,j)*ay(i,j)-s1)*adipy(j)+
     &                    s2*ay(i,j)*az(i,j)*adipz(j))
            zf(i)=zf(i)+ (s2*ax(i,j)*az(i,j)*adipx(j)+
     &                    s2*ay(i,j)*az(i,j)*adipy(j)+
     &                   (s2*az(i,j)*az(i,j)-s1)*adipz(j))
          endif
        enddo
      enddo
c
c     calculation for intramolecular induced dipoles
c
      do i=1,iuniq
       xdip(i)=alpi(i)*xf(i)
       ydip(i)=alpi(i)*yf(i)
       zdip(i)=alpi(i)*zf(i)
      enddo
      return
      end
***********************************************************************
!* Ref.: "Méthodes de calcul numérique - Tome 2 - By Claude    *
!*        Nowakowski, PSI Editions, 1984".                     *
!*                                                             *
!*                  Fortran 90 Version By J-P Moreau, Paris.   *
!*                            (www.jpmoreau.fr)                *
!***************************************************************
CProgram Householder
!*  1. dimension up
!*     dimension*8 A(50, 100), A0(50, 50), V(50) -> dimension*8 A(50, 100), A0(50, 50), V(50)  
!*  2. inverse matrix -> Alpha_mol AM
!*
*********1*********2*********3*********4*********5*********6*********7**
      Subroutine houhol(is,if,nn,C,E,AM,DI)
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
C     integer n   !size of matrix A
      dimension A(i200*3, i200*2*3), A0(i200*3, i200*3), V(i200*3)
      dimension B(i200*3), AI(i200*3, i200*3), 
     &  AP(i200*3, i200*3), X(i200*3)
      dimension AM(3, 3), F(i340*i20, i200*3),  
     &  DD(i340*i20, i200*3), DI(i340*i20,3)
      dimension C(i200*3,i200*3),E(i340*i20, i200*3)
      data aukcal /627.509541d0/,zero/0.0d0/,debye/2.54158059d0/
C    !read matrix to be inverted
C    !Open input/output text files
!*
      n=nn*3
      do i=1,n
        do j=1,n
          A(i,j)=C(i,j)
        enddo
      enddo
!* 
      l=is*if
      do i=1,l
        do j=1,n
          F(i,j)=E(i,j)
        enddo
      enddo
!* 
      do i = 1, n
        do j = 1, n
          A(i, j + n) = 0.d0
          A0(i, j) = A(i, j)
        end do
        A(i, i + n) = 1.d0
      end do
C   !Transform A into triangular matrix
      call S1000(n,A,V)
C   !N linear systems to solve
      do k = 1, n
        do i = 1, n
          B(i) = A(i, K + n)
        end do
C     !Solve triangular system
        call S2000(n,A,B,X)
        do i = 1, n
          AI(i, K) = X(i)
        end do
      end do !k loop
C   !Calculate determinant
      call S4000(n,A,D)
C   !Check AP=A*A^-1=Unity Matrix
      call S3000(n,A0,AI,AP)
!
!   Alpha_mol
!
C   ! TENSOR AM   lower matrix elements
      write(7,*) ' '
      call molpol(n,AI,AM)
      TAM=(AM(1,1)+AM(2,2)+AM(3,3))/3.0d0
      write(7,30) AM(1,1),AM(2,1),AM(2,2),AM(3,1),AM(3,2),AM(3,3),TAM
!
!  Dipole moment 
!
      write(7,*) ' '
      call dipf(n,l,AI,F,DD,DI)
!
!
   10 format(' ',f10.5)
   20 format('Determinant: ',f10.5)
   30 format(' Molecular Polarizability (Bohr**3):',6f10.4,
     &   '  isotropic',f10.4)
      return
      End !of main program

      Subroutine S1000(n,A,V)
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)

C integer n
      dimension A(i200*3,i200*2*3),V(i200*3)
C     dimension S,XA,XB,XG,XL
      do k = 1, n - 1
        S = 0.d0
        do i = k, n
          S = S + A(i, K) * A(i, K)
        end do

        XL = -sign(1.d0,A(K, K)) * dsqrt(S)
        XA = XL * (XL - A(K, K))
        V(K) = A(K, K) - XL
        do i = k + 1, n
          V(i) = A(i, K)
        end do
        do j = k + 1, 2*n
          XB = 0
          do i = k, n
            XB = XB + V(i) * A(i, j)
          end do
          XG = XB / XA
          do i = k, n
            A(i, j) = A(i, j) - XG * V(i)
          end do
        end do
        A(K, K) = XL
        do i = k + 1, n
          A(i, K) = 0.d0
        end do
      end do !k loop
      return
      End

      Subroutine S2000(n,A,B,X)
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
C  integer n
      dimension A(i200*3,i200*2*3),B(i200*3),X(i200*3)
C     integer i, j
C     dimension S
      do i = n, 1, -1
        j = n
        S = 0.d0
   20   If (i == j) goto 40
         S = S - A(i, j) * X(j)
         j = j - 1
         goto 20
   40   X(i) = (B(i) + S) / A(i, i)
      end do
      return
      End

      Subroutine S3000(n,A0,AI,AP)
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
C integer n
      dimension A0(i200*3,i200*3),AI(i200*3,i200*3),AP(i200*3,i200*3)
C     dimension S
      do i = 1, n
        do j = 1, n
          S = 0.d0
          do k = 1, n
            S = S + A0(i, k) * AI(k, j)
          end do
          AP(i, j) = S
        end do
      end do
      return
      End

      Subroutine S4000(n,A,D)
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)

C integer n
      dimension A(i200*3,i200*2*3)
C , D
      D = 1.d0
      do i = 1, n
        D = D * A(i, i)
      end do
      return
      End
 
      Subroutine molpol(n,A,AM)
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)

C integer n
      dimension A(i200*3,i200*2*3), AM(3,3), A1(i200*3,3)
      k=1
      do i = 1, n
        A1(i,1)=0.0d0
        A1(i,2)=0.0d0
        A1(i,3)=0.0d0
        do j = 1,n,3
          A1(i,1)=A1(i,1)+A(i,j) 
          A1(i,2)=A1(i,2)+A(i,j+1) 
          A1(i,3)=A1(i,3)+A(i,j+2) 
        end do
      end do
      do i=1,3
        do j=1,3
          AM(i,j)=0.0d0
        end do
      end do
      do i = 1, n, 3
        do j = 1,3
          AM(1,j)=AM(1,j)+A1(i,j)
          AM(2,j)=AM(2,j)+A1(i+1,j)
          AM(3,j)=AM(3,j)+A1(i+2,j)
        end do
      end do
      return
      End

      Subroutine dipf(n,m,A,F,D,DI)
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)

C integer n, m
      dimension A(i200*3,i200*3),F(i340*i20,i200*3),
     &  D(i340*i20,i200*3),DI(i340*i20,3)
      k=1
      do i = 1, m
        do j = 1, n
          D(i,j)=0.0d0
        end do
      end do
      do i = 1,n
        do j=1,m
          do k=1,n
            D(j,k)=D(j,k)+A(i,k)*F(j,k)
          end do
        end do
      end do
!
      do i = 1, m
        do j = 1, 3
          DI(i,j)=0.0d0
        end do
      end do
      do i=1,m
        do j=1,n,3
          DI(i,1)=DI(i,1)+D(i,j)
          DI(i,2)=DI(i,2)+D(i,j+1)
          DI(i,3)=DI(i,3)+D(i,j+2)
        end do
      end do
      return
      End

!end of file householder.f90
***********************************************************************
      subroutine imdflds(ll,kk,alpi,xdip,ydip,zdip)
c                             
c     irest=3  Matrix Inversion
c
c     calculations for intramolecular dipole field
c     simple scaling model
c
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
      common/chcor/iateq(i200,4),jateq(i200,4),neq(5),
     & iuniq,nvar,ifix
     & ,iatom(i200),ivary(i200,4),xkpt(i200),ykpt(i200),zkpt(i200)
     & ,qkpt (i200)
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      common/tens/
     &  txx(i200,i200),txy(i200,i200),
     &  txz(i200,i200),
     &  tyx(i200,i200),tyy(i200,i200),
     &  tyz(i200,i200),
     &  tzx(i200,i200),tzy(i200,i200),
     &  tzz(i200,i200)
      dimension adipx(i200),adipy(i200),adipz(i200),alpi(i200),
     &          xdip(i200),ydip(i200),zdip(i200),
     &          xf(i200),yf(i200),zf(i200)
!*--------------------------------------------------------------------------
      dimension A(i200*3,i200*3),DA(i200*3)
!*--------------------------------------------------------------------------
      data zero/0.0d0/,tho/2.089d0/,one/1.0d0/,two/2.0d0/,six/6.0d0/,
     & three/3.0d0/
c
      do i=1,iuniq
        do j=1,iuniq
          if(i.ne.j)then
            r1=ar(i,j)
            r2=r1*r1
            r3=r1*r2
            r5=r1*r2*r2
            s1=one/r3
            s2=three/r5
C
            txx(i,j)=-(s2*ax(i,j)*ax(i,j)-s1)*flnkd(i,j)
            txy(i,j)=-(s2*ax(i,j)*ay(i,j))*flnkd(i,j)
            txz(i,j)=-(s2*ax(i,j)*az(i,j))*flnkd(i,j)
            tyx(i,j)=-(s2*ay(i,j)*ax(i,j))*flnkd(i,j)
            tyy(i,j)=-(s2*ay(i,j)*ay(i,j)-s1)*flnkd(i,j)
            tyz(i,j)=-(s2*ay(i,j)*az(i,j))*flnkd(i,j)
            tzx(i,j)=-(s2*az(i,j)*ax(i,j))*flnkd(i,j)
            tzy(i,j)=-(s2*az(i,j)*ay(i,j))*flnkd(i,j)
            tzz(i,j)=-(s2*az(i,j)*az(i,j)-s1)*flnkd(i,j)
            
C
          else
            txx(i,j)=one/alpi(i)
            txy(i,j)=zero
            txz(i,j)=zero
            tyx(i,j)=zero
            tyy(i,j)=one/alpi(i)
            tyz(i,j)=zero
            tzx(i,j)=zero
            tzy(i,j)=zero
            tzz(i,j)=one/alpi(i)
          endif
        enddo
      enddo
!*----------------------------------------------------------
!*    Matrix A
!*----------------------------------------------------------
      k=1
      do i=1,iuniq
        l=1
        do j=1,iuniq
          A(k,l)=txx(i,j)
          A(k,l+1)=txy(i,j)
          A(k,l+2)=txz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+1,l)=tyx(i,j)
          A(k+1,l+1)=tyy(i,j)
          A(k+1,l+2)=tyz(i,j)
          l=l+3
        enddo
        l=1
        do j=1,iuniq
          A(k+2,l)=tzx(i,j)
          A(k+2,l+1)=tzy(i,j)
          A(k+2,l+2)=tzz(i,j)
          l=l+3
        enddo
        k=k+3
      enddo
!*
      call houhol0(ll,kk,iuniq,A,DA)
!*   
!*----------------------------------------------------------
c     calculation for intramolecular induced dipoles
c 
      j=0
      do i=1,iuniq*3,3
       j=j+1
       xdip(j)=DA(i)
       ydip(j)=DA(i+1)
       zdip(j)=DA(i+2)
      enddo 
      return
      end
***********************************************************************
      Subroutine houhol0(ll,kk,nn,C,DI)
      implicit real*8(a-h,o-z)
      parameter (i200=60,i340=240,i20=30,i100=98)
C     integer n   !size of matrix A
      common/poplnk/flnkc(i200,i200),flnkd(i200,i200),facc(3),facd(3),
     &  pchin(i200),cfldx(i200),cfldy(i200),cfldz(i200),
     &  pfldx(i20,i340,i200),pfldy(i20,i340,i200),pfldz(i20,i340,i200),
     &  fldx0(i20,i340,i200),fldy0(i20,i340,i200),fldz0(i20,i340,i200),
     &  ax(i200,i200),ay(i200,i200),az(i200,i200),ar(i200,i200),factt
      dimension A(i200*3, i200*2*3), A0(i200*3, i200*3), V(i200*3)
      dimension B(i200*3), AI(i200*3, i200*3), 
     &  AP(i200*3, i200*3), X(i200*3)
      dimension  F(i200*3), DI(i200*3),C(i200*3,i200*3)
!*
!*    A Matrix to be inverted
!*
      n=nn*3
      do i=1,n
        do j=1,n
          A(i,j)=C(i,j)
        enddo
      enddo
!*     
!*    is if F Field from a fld charge
      k=1
      do i=1,n,3
        F(i)=pfldx(ll,kk,k)
        F(i+1)=pfldy(ll,kk,k)
        F(i+2)=pfldz(ll,kk,k)
        k=k+1
      enddo
!*
!*    A Matrix 
!* 
      do i = 1, n
        do j = 1, n
          A(i, j + n) = 0.d0
          A0(i, j) = A(i, j)
        end do
        A(i, i + n) = 1.d0
      end do
c
C     !Transform A into triangular matrix
      call S1000(n,A,V)
C     !N linear systems to solve
      do k = 1, n
        do i = 1, n
          B(i) = A(i, K + n)
        end do
C       !Solve triangular system
        call S2000(n,A,B,X)
        do i = 1, n
          AI(i, K) = X(i)
        end do
      end do !k loop
C     !Calculate determinant
      call S4000(n,A,D)
C     !Check AP=A*A^-1=Unity Matrix
      call S3000(n,A0,AI,AP)
!*
!*    AI  Inverted Matrix
!*    DI  Induced Dipole Moment/atom  a.u. 
!*    F   Field from a fld charge 
      do i=1,n
        DI(i)=0.0D0
      enddo
      do i=1,n
        do j=1,n
          DI(i)=DI(i)+AI(i,j)*F(j)        
        enddo
      enddo
      return
      End