C     program IDFF_mi_ctat
C
C     Iduced Dipole Force Field 
C     
C     Plarizable Molecular Block (PMB) model       (Matrix Inversion)
C       
C     idff_intra_2020_lin --> idff_minv_2022
C     
C     May 1994 - May 2022   by S.Nakagawa
C
C     2023  2/26 Cleaned up the source code.
C           2/26 A model in which the positive charge of the induced dipole is
C                placed on the atom is added (ictat=1).
C     ---------------------------------------------------------
C   
C     Intramolecular polarization model (Induced point dipole model)
C
C     irest
C           = 3 intra & intermolecular atom polarizability fit Simple scaling
C           = 4 Thole Exponential model  van Duijnen type a**3/8piexp(-au)   default a=1.9088 (Thole 2.089)
C           = 5 Thole Exponential model  Ponder type      3a/4piexp(-au**3)  default a=0.39 (Thole 0.572)
C           =13 Molecular Block option for charge-induced dipole applied to Simple 1-2 scaling 
C               (intra & intermolecular atom polarizability fit Simple scaling)
C           =14 Molecular Block option for charge-induced dipole applied to Thole molde van Duijnen type
C           =15 Molecular Block option for charge-induced dipole applied to Thole molde Ponder type
C           =143 Molecular Block option for charge-induced dipole applied to Thole molde van Duijnen type
C                with simple scaling
C           =153 Molecular Block option for charge-induced dipole applied to Thole molde Ponder type
C                with simple scaling
C           =163 Molecular Block option for charge-induced dipole applied to Gaussian screening  Masia type
C                with simple scaling  (0.8456)
C           =173 Molecular Block option for charge-induced dipole applied to Thole Linear model 
C                with simple scaling
C           =24 Molecular Block option for es 
C               + block option for charge-induced dipole
C               (intra & intermolecular atom polarizability fit Simple scaling)
C     
C     Input data
C
C       fort.1 (iu1) / topology file
C         1) Residue_name total_charge #_bond_charge irest plz_sfactor (segment)
C            (a4,f10.5,2i4,f10.5 (i3))
C         2) #_atom_plz #_bond_plz #_bond_center_plz #_linkatm 
C         3) atom charge atm_plz 0.0 0.0 0.0 rmin emin (2a2,7f10.5)
C         4) bond_1 bond_2 bond_plz 0.0 0.0 0.0 bond_charge (*)
C         5) Title                                          (*)
C            cfld dfld factor for 1-5... and  es_14 vdw_14_factor    cfld: charge-induced dipole
C            (default -1.0=>1/1.2 for es_14 & -1.0=>0.5 for vdw_14)  dfld: induced dipole-induced dipole
C            1.0>=factor>=0.0                               (*)
C         6) Title                                          (*)
C            cfld dfld                                      (*)
C            link12                                         (50i4)
C         7) Title                                          (*)
C            cfld dfld                                      (*)
C         8) Title                                          (*)
C            cfld dfld                                      (*)
C
C       fort.2 (iu2) / sequence file 
C
C       fort.3 (iu3) / Complex molecule Rest_mp2 coordinate and esp file 
C
C       fort.4       scp data and print option (f6.3,a5)  1.0 or 0.9 etc. 'print' or '     '
C
C       fort.10      / output data
C
C      #
C      foreach scp (09)
C      foreach dis(030)
C      foreach dat (xxx)
C      cp b3_GGGNa3CCCNa3_{$dis}.rest REST_1
C      ln -s topseg.inp fort.1
C      ln -s seq_{$dat}.inp fort.2
C      ln -s REST_1 fort.3
C      ln -s scp_{$scp}.dat fort.4
C      ~/source/EspOpt/idff_mi.exe
C      cp fort.10 {$scp}_{$dat}_{$dis}.out
C      rm -fr fort.*
C      rm -f REST_1
C      end
C      end
C      end
C
      call cpu_time(tt1)
C
      call rddata
C     write(6,*)' % call rddata end'
      call precoo
C     write(6,*)' % call precoo end'
      call predis
C     write(6,*)' % call predis end'
      call elesta
C     write(6,*)' % call elesta end'
      call vdwene
C     write(6,*)' % call vdwene end'
      call cnsdip
C     write(6,*)' % call cnsdip end'
      call fldcrg
C     write(6,*)' % call fldcrg end'
      call plzcyc
C     write(6,*)' % call plzcyc end'
      call popsrf
C     write(6,*)' % call popsrf end'
      write(10,*) '  '
      write(10,*) ' stop idff_minv_ctat '
C
      call cpu_time(tt2)
      write(10,*) '  '
      write(10,*) 'CPU time:',tt2-tt1,' second'
C
      stop
      end
************************************************************************
      subroutine rddata
C     
C     ReaD DATA
C
      iu1=1
      iu2=2
      iu3=3
C
C     open(unit=iu1,form='formatted',status='old',readonly)
C
      call rdtop(iu1)
      call rdseq(iu2)
      call rdesp(iu3)
C
C     close(unit=iu1)
C
      return
      end
************************************************************************
      subroutine rdtop(iu1)
      parameter(k100=30,k50=250)
      parameter(kr50=50,ka200=300,k4=4,k21=21)
C 
C     k100 Number of residues
C      k50 atoms or bonds
C       k4 alpha, beta/2, gamma/6, delta/24 polz and hyperpolz parameters
C
      implicit real*8 (a-h,o-z)
      character*5 prin
      character*4 resi,dummy
      character*2 atname(k100,k50,2),atom1(k21),atom2(3)
      dimension atmmass(k21)
      common/pr/ prin
      common/resid/ numres,resi(kr50)
      common/topol1/ tcrg(kr50),crg(kr50,ka200),
     &  atmplz(kr50,ka200,k4),bonplz(kr50,ka200,k4),natmplz(kr50),
     &  nbonplz(kr50),nbon1(kr50,ka200),nbon2(kr50,ka200),
     &  atmass(kr50,ka200),
     &  xrmin(kr50,ka200),xemin(kr50,ka200),mbon1(kr50,ka200),
     &  mbon2(kr50,ka200),boaplz(kr50,ka200,4),mbonplz(kr50),
     &  mseg(kr50,ka200)
      common/topol2/
     &  nlnk(kr50),nl12(kr50,ka200,ka200),nl13(kr50,ka200,ka200),
     &  nl14(kr50,ka200,ka200),scac14(kr50),scav14(kr50),
     &  facp12(kr50),facp13(kr50),facp14(kr50),facp15(kr50),
     &  facc12(kr50),facc13(kr50),facc14(kr50),facc15(kr50),ibonc(kr50),
     &  bchg(kr50,ka200),irest(kr50),psfac(kr50),factt(kr50)

      dimension rmin(k21),emin(k21)
      CHARACTER*4 title(20)
      data zero/0.0d0/,au/627.509541D0/,bohr/0.52917706D0/,
     &     one/1.0d0/
      data atom1/'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne',
     &                     'Na','Mg','Al','Si','P ','S ','Cl','Ar'
     &                    ,'K ','Ca','X '/
      data atom2/'  ','+ ','- '/
      data atmmass/1.00794d0,4.002602d0,
     & 6.941d0,9.012182d0,10.811d0,12.011d0,14.00674d0,15.9994d0,
     & 18.9984032d0,20.1797d0,
     & 22.989768d0,24.3050d0,26.981539d0,28.0855d0,30.973762d0,
     & 32.066d0,35.4527d0,39.948d0,
     & 39.0983d0,40.0780d0,0.0d0/
C
C     topology    k50 residues are read
C
C
C     read scp (scale factor for plz parm)
C
      read(4,*) scp
      if(scp.lt.0.0d0) scp=1.0d0
      write(10,'(a34,f10.5)') "scale factor for plz parm(scp) is ",scp
      write(10,'(a5)') '    '
C
C    scp : scale factor for plz   read from fort.4
C    psfac(i) : scale factor for plz  read from iu1 parameter&topology file
C
      numres=0
      do i=1,k50
        numres=numres+1
        read(iu1,'(A4,f10.5,2i4,f10.5)') resi(i),tcrg(i),ibonc(i),
     &  irest(i),psfac(i)
        if(psfac(i).le.zero.or.psfac(i).gt.one) psfac(i)=one
        read(iu1,*) na,nb,mb,nlnk(i)
        natmplz(i)=na
        nbonplz(i)=nb
        mbonplz(i)=mb
        totc=zero
        write(10,'(i3,A5,f10.5,i4,a8,i4,f10.5,4i6)') i,resi(i),tcrg(i),
     &   ibonc(i),
     &   ' #irest=',irest(i),psfac(i),na,nb,mb,nlnk(i)
        do j=1,na
          if(irest(i).eq.13.or.irest(i).eq.24.or.irest(i).eq.14.or.
     &      irest(i).eq.15.or.irest(i).eq.143.or.
     &      irest(i).eq.153.or.irest(i).eq.163.or.irest(i).eq.173)then
            read(iu1,'(2a2,7f10.5,i3)') (atname(i,j,k),k=1,2),crg(i,j),
     &        (atmplz(i,j,k),k=1,k4),xrmin(i,J),xemin(i,j),mseg(i,j)
            do k=1,k4
              atmplz(i,j,k)=atmplz(i,j,k)*scp
            enddo
          else
            read(iu1,'(2a2,7f10.5)') (atname(i,j,k),k=1,2),crg(i,j),
     &        (atmplz(i,j,k),k=1,k4),xrmin(i,J),xemin(i,j)
            do k=1,k4
              atmplz(i,j,k)=atmplz(i,j,k)*scp
            enddo
            mseg(i,j)=0
          endif
          xrmin(i,j)=xrmin(i,j)/bohr
          xemin(i,j)=xemin(i,j)/au
          totc=totc+crg(i,j)
          do k=1,k21
            if(atname(i,J,1).eq.atom1(k))then
              atmass(i,j)=atmmass(k)
            endif
          enddo
        enddo
        do j=1,nb
          read(iu1,*) nbon1(i,j),nbon2(i,j),
     &      (bonplz(i,j,k),k=1,k4),bchg(i,j)
          totc=totc+bchg(i,j)
        enddo
        do j=1,mb
          read(iu1,*) mbon1(i,j),mbon2(i,j),
     &      (boaplz(i,j,k),k=1,k4)
        enddo
C
C      check total charge
C
        if(totc.ne.zero)then
          write(10,'(a20,f10.5)') 
     &    '   #Total charge is ',totc 
        endif
C
C       Intramolecular interactin
C         scac14  scale factor for 1-4 electrostatic
C         scav14  scale factor for 1-4 vdw
C
C         facc12  scale factor for 1-2 charge-induced dipole
C         facc13  scale factor for 1-3 charge-induced dipole
C         facc14  scale factor for 1-4 charge-induced dipole
C         facc15  scale factor for 1-5... charge-induced dipole
C
C         facp12  scale factor for 1-2 induced dipole-induced dipole (polarization)
C         facp13  scale factor for 1-3 induced dipole-induced dipole (polarization)
C         facp14  scale factor for 1-4 induced dipole-induced dipole (polarization)
C         facp15  scale factor for 1-5... induced dipole-induced dipole (polarization)
C
C         factt   Thole type a factor
C
        if(nlnk(i).gt.0)then
C         write(6,*) 
C    &    '#topology  Read Intramolecular 1-5 1-2 1-3 1-4 data!'
          read(iu1,'(20a4)') title
          read(iu1,*) facc15(i),facp15(i),scac14(i),scav14(i),factt(i)
          read(iu1,'(20a4)') title
          read(iu1,*) facc12(i),facp12(i)
          do j=1,nlnk(i)
            read(iu1,'(200i4)') (nl12(i,j,k),k=1,nlnk(i))
          enddo
C%SN
C%SN
          write(10,'(5(a10,f10.5))') 
     &      '   facc15=',facc15(i),'   facp15=',facp15(i),
     &      '   scac14=',scac14(i),'   scav1 =',scav14(i),
     &      '   factt =',factt(i)
          write(10,'(a5)') '     '
C%SN
          call lnkplz(i)
C
          read(iu1,'(20a4)') title
          read(iu1,*) facc13(i),facp13(i)
C         do j=1,nlnk(i)
C           read(iu1,'(200i4)') (nl13(i,j,k),k=1,nlnk(i))
C         enddo
          read(iu1,'(20a4)') title
          read(iu1,*) facc14(i),facp14(i)
C         do j=1,nlnk(i)
C           read(iu1,'(200i4)') (nl14(i,j,k),k=1,nlnk(i))
C         enddo
        endif
C
        read(iu1,'(a4)') dummy
        if(dummy.eq.'end '.or.dummy.eq.'END '.or.dummy.eq.'End ')then
          go to 999
        else if(dummy.ne.'    ') then
          write(10,*) 'Read Error in Topology File!! STOP!!'
          STOP
        endif
      enddo
C
C     write(6,*) '#irest ',(irest(i),i=1,numres)  
C      
  999 continue
      return
      end
************************************************************************
      subroutine rdseq(iu2)
      parameter(ks30=30)
C 
C     read sequence of cluster
C
C     ex)   3
C          NA+ METNMETN    (10a4)
C
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      read(iu2,*) nseq
      read(iu2,'(10a4)') (seq(i),i=1,nseq)
      write(10,'(a8,i5)') 'Sequence',nseq
      write(10,'(20(a4,1x))') (seq(i),i=1,nseq)
      return
      end
************************************************************************
      subroutine rdesp(iu3)
      parameter(k300=300,k3000=5000)
C 
C     read coordinates and esp data of cluster
C  
C     (g09 Rest_scf or Rest_mp2)
C
C     cx,cy,cz 100 atomic coordinates of cluster
C     sx,sy,sz 500 coordinates on 1.8 times vdw surface
C     spot     500 Surface esp potential from ab initio calculations
C
      implicit real*8 (a-h,o-z)
      common/espot/ cx(k300),cy(k300),cz(k300),sx(k3000),sy(k3000),
     &  sz(k3000),spot(k3000),gcrg(k300),icoor,isurf,
     &  pop(k3000),gpop(k3000)
      read(iu3,*) icoor
      do i=1,icoor
        read(iu3,'(2x,4(f16.10,2x))') cx(i),cy(i),cz(i),gcrg(i)
C       write(6,'(2x,4(f16.10,2x))') cx(i),cy(i),cz(i),gcrg(i)
      enddo
      read(iu3,*) isurf
      do i=1,isurf
        read(iu3,'(4(f15.9))') sx(i),sy(i),sz(i),spot(i)
      enddo
      return
      end
************************************************************************
      subroutine precoo
      parameter(k300=300,k50=250,k3000=5000,k100=30)
      parameter(ks30=30,k200=300)
      parameter(kr50=50,ka200=300,k4=4,k21=21)
C 
C     preparation of working coordinates and parameters
C  
C                 mol coor
C     wx,wy,wz 100,50 working coordinates of atom 
C
C     wxmass,wymass,wzmass 100 mass center coordinates of residue
C
      implicit real*8 (a-h,o-z)
      character*5 prin
      character*4 resi,dummy
      character*4 seq(ks30)
      common/pr/ prin
      common/resid/ numres,resi(kr50)
      common/topol1/ tcrg(kr50),crg(kr50,ka200),
     &  atmplz(kr50,ka200,k4),bonplz(kr50,ka200,k4),natmplz(kr50),
     &  nbonplz(kr50),nbon1(kr50,ka200),nbon2(kr50,ka200),
     &  atmass(kr50,ka200),
     &  xrmin(kr50,ka200),xemin(kr50,ka200),mbon1(kr50,ka200),
     &  mbon2(kr50,ka200),boaplz(kr50,ka200,4),mbonplz(kr50),
     &  mseg(kr50,ka200)
      common/topol2/
     &  nlnk(kr50),nl12(kr50,ka200,ka200),nl13(kr50,ka200,ka200),
     &  nl14(kr50,ka200,ka200),scac14(kr50),scav14(kr50),
     &  facp12(kr50),facp13(kr50),facp14(kr50),facp15(kr50),
     &  facc12(kr50),facc13(kr50),facc14(kr50),facc15(kr50),ibonc(kr50),
     &  bchg(kr50,ka200),irest(kr50),psfac(kr50),factt(kr50)
      common/espot/ cx(k300),cy(k300),cz(k300),sx(k3000),sy(k3000),
     &  sz(k3000),spot(k3000),gcrg(k300),icoor,isurf,
     &  pop(k3000),gpop(k3000)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/sequ/ nseq,seq
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
      dimension flnk(k50,k50)
      data one/1.0d0/,zero/0.0d0/
      write(10,*) nseq,numres
      k1=0
      do i=1,nseq
        kk=0
        do j=1,numres
          if(seq(i).eq.resi(j))then
            k2=0
            nwaplz(i)=natmplz(j)
            nwbplz(i)=nbonplz(j)
            mwbplz(i)=mbonplz(j)
            nwlnk(i)=nlnk(j)
            iwbonc(i)=ibonc(j)
            iwrest(i)=irest(j)
            wfacp(i)=factt(j)
C%SN
C%SN  2019/8/14
C%SN
Cxxx        facp15(j)=1.0d0
C%SN
            do k=1,natmplz(j)
              k1=k1+1
              k2=k2+1
              wx(i,k2)=cx(k1)
              wy(i,k2)=cy(k1)
              wz(i,k2)=cz(k1)
              wcrg(i,k2)=crg(j,k)
              nseg(i,k2)=mseg(j,k)
              wamass(i,k2)=atmass(j,k)
              wrmin(i,k2)=xrmin(j,k)
              wemin(i,k2)=xemin(j,k)
C
C       !!!! scale factor 
C
              do l=1,k4
                wplz(i,k2,l)=psfac(j)*atmplz(j,k,l)
              enddo
            enddo
            do k=1,nbonplz(j)
              k2=k2+1
              wx(i,k2)=(wx(i,(nbon1(j,k)))+wx(i,(nbon2(j,k))))*0.5d0
              wy(i,k2)=(wy(i,(nbon1(j,k)))+wy(i,(nbon2(j,k))))*0.5d0
              wz(i,k2)=(wz(i,(nbon1(j,k)))+wz(i,(nbon2(j,k))))*0.5d0
              wcrg(i,k2)=bchg(j,k)
              nwbon1(i,k)=nbon1(j,k)
              nwbon2(i,k)=nbon2(j,k)
C
C       !!!! scale factor 
C
              do l=1,k4
                wplz(i,k2,l)=psfac(j)*bonplz(j,k,l)
              enddo
            enddo
            do k=1,mbonplz(j)
              k2=k2+1
              wx(i,k2)=(wx(i,(mbon1(j,k)))+wx(i,(mbon2(j,k))))*0.5d0
              wy(i,k2)=(wy(i,(mbon1(j,k)))+wy(i,(mbon2(j,k))))*0.5d0
              wz(i,k2)=(wz(i,(mbon1(j,k)))+wz(i,(mbon2(j,k))))*0.5d0
              wcrg(i,k2)=0.0d0
              mwbon1(i,k)=mbon1(j,k)
              mwbon2(i,k)=mbon2(j,k)
              do l=1,k4
                wplz(i,k2,l)=boaplz(j,k,l)
              enddo
            enddo
C
C           intramolecular scale factor for es scac & vdw scav 
C              plz facc (charge-induced dipole)
C              plz facp (induced dipole-induced dipole) 
C 
            kk=nlnk(j)
            if(kk.gt.0) then
C
              if(scac14(j).gt.1.0d0.or.scac14(j).lt.0.0d0) then
                scac14(j)=1.0d0/1.2d0
              endif
              call fscale(i,j,kk,zero,zero,zero,one,flnk) 
              do ii=1,kk
                do jj=1,kk
                  scac(i,ii,jj)=flnk(ii,jj)
                  if(iwrest(i).eq.24)then
                    if(nseg(i,ii).eq.nseg(i,jj)) scac(i,ii,jj)=0.0d0
                  endif
                enddo
                if(iwrest(i).eq.24.or.iwrest(i).eq.13.or.
     &            iwrest(i).eq.14.or.iwrest(i).eq.15.or.
     &            iwrest(i).eq.143.or.iwrest(i).eq.153.or.
     &            iwrest(i).eq.163.or.iwrest(i).eq.173.or.
     &            iwrest(i).eq.3)then
                  if(prin.eq.'print')then
                    write(10,'(a6,150f2.0)')
     &                '%scac ',(scac(i,ii,jj),jj=1,kk)
                  endif
                endif
              enddo
              call fscale(i,j,kk,zero,zero,scac14(j),zero,flnk) 
              do ii=1,kk
                do jj=1,kk
                  se14(i,ii,jj)=flnk(ii,jj)
                  if(iwrest(i).eq.24)then
                    if(nseg(i,ii).eq.nseg(i,jj)) se14(i,ii,jj)=0.0d0
                  endif
                enddo
                if(iwrest(i).eq.24.or.iwrest(i).eq.13.or.
     &            iwrest(i).eq.14.or.iwrest(i).eq.15.or.
     &            iwrest(i).eq.143.or.iwrest(i).eq.153.or.
     &            iwrest(i).eq.163.or.iwrest(i).eq.173.or.
     &            iwrest(i).eq.3)then
                  if(prin.eq.'print')then
                    write(10,'(a6,150f2.0)')
     &                '%se14 ',(se14(i,ii,jj),jj=1,kk)
                  endif
                endif
              enddo
C
              if(scav14(j).gt.1.0d0.or.scav14(j).lt.0.0d0) then
                scav14(j)=0.5d0
              endif
              call fscale(i,j,kk,zero,zero,zero,one,flnk) 
              do ii=1,kk
                do jj=1,kk
                  scav(i,ii,jj)=flnk(ii,jj)
                enddo
                if(iwrest(i).eq.24.or.iwrest(i).eq.13.or.
     &            iwrest(i).eq.14.or.iwrest(i).eq.15.or.
     &            iwrest(i).eq.143.or.iwrest(i).eq.153.or.
     &            iwrest(i).eq.163.or.iwrest(i).eq.173.or.
     &            iwrest(i).eq.3)then
                  if(prin.eq.'print')then
                    write(10,'(a6,150f2.0)')
     &                '#scav ',(scav(i,ii,jj),jj=1,kk)
                  endif
                endif
              enddo
              call fscale(i,j,kk,zero,zero,scav14(j),zero,flnk) 
              do ii=1,kk
                do jj=1,kk
                  sv14(i,ii,jj)=flnk(ii,jj)
                enddo
                if(iwrest(i).eq.24.or.iwrest(i).eq.13.or.
     &            iwrest(i).eq.14.or.iwrest(i).eq.15.or.
     &            iwrest(i).eq.143.or.iwrest(i).eq.153.or.
     &            iwrest(i).eq.163.or.iwrest(i).eq.173.or.
     &            iwrest(i).eq.3)then
                  if(prin.eq.'print')then
                    write(10,'(a6,150f2.0)')
     &                '#sv14 ',(sv14(i,ii,jj),jj=1,kk)
                  endif
                endif
              enddo
C
              call fscale
     &          (i,j,kk,facp12(j),facp13(j),facp14(j),facp15(j),flnk) 
              do ii=1,kk
                do jj=1,kk
                  facp(i,ii,jj)=flnk(ii,jj)
                enddo
                if(iwrest(i).eq.24.or.iwrest(i).eq.13.or.
     &            iwrest(i).eq.14.or.iwrest(i).eq.15.or.
     &            iwrest(i).eq.143.or.iwrest(i).eq.153.or.
     &            iwrest(i).eq.163.or.iwrest(i).eq.173.or.
     &            iwrest(i).eq.3)then
                  if(prin.eq.'print')then
                    write(10,'(a6,150f2.0)')
     &                '#facp ',(facp(i,ii,jj),jj=1,kk)
                  endif
                endif
              enddo
C
              call fscale
     &          (i,j,kk,facc12(j),facc13(j),facc14(j),facc15(j),flnk)
              do ii=1,kk
                do jj=1,kk
                  facc(i,ii,jj)=flnk(ii,jj)
                  if(iwrest(i).eq.13.or.iwrest(i).eq.24.or.
     &              iwrest(i).eq.14.or.iwrest(i).eq.15.or.
     &              iwrest(i).eq.143.or.iwrest(i).eq.153.or.
     &              iwrest(i).eq.163.or.iwrest(i).eq.173)then
                    if(nseg(i,ii).eq.nseg(i,jj)) facc(i,ii,jj)=0.0d0
                  endif
                enddo
                if(iwrest(i).eq.24.or.iwrest(i).eq.13.or.
     &            iwrest(i).eq.14.or.iwrest(i).eq.15.or.
     &            iwrest(i).eq.143.or.iwrest(i).eq.153.or.
     &            iwrest(i).eq.163.or.iwrest(i).eq.173.or.
     &            iwrest(i).eq.3)then
                  if(prin.eq.'print')then
                    write(10,'(a6,150f2.0)')
     &                '#facc ',(facc(i,ii,jj),jj=1,kk)
                  endif
                endif
              enddo
C
              write(10,*) 
     &          'scale factor es1-4 vdw1-4 induced dipole-induced dipole
     & 1-2 1-3 1-4 1-5...       charge-induced diple 1-2 1-3 1-4 1-5...'
              write(10,'(a4,2f10.5,4x,4f10.5,4x,4f10.5)') 
     &          resi(j),scac14(j),scav14(j),
     &          facp12(j),facp13(j),facp14(j),facp15(j),
     &          facc12(j),facc13(j),facc14(j),facc15(j)
            endif
C   
            go to 10
          else
            kk=kk+1
            if(kk.eq.numres) then 
              write(10,'(a30,a4)') '!! Can not find the residue!! ',
     &          seq(i)
              STOP
            endif
          endif
        enddo
   10 continue
      enddo
C
C
C     Mass Center of Residue
C
      xmass=zero
      ymass=zero
      zmass=zero
      tmass=zero
      do i=1,nseq
        wxmass(i)=zero
        wymass(i)=zero
        wzmass(i)=zero
        wtmass=zero
        na=nwaplz(i)
        do j=1,na
          x=wx(i,j)*wamass(i,j)
          y=wy(i,j)*wamass(i,j)
          z=wz(i,j)*wamass(i,j)
          wxmass(i)=wxmass(i)+x
          wymass(i)=wymass(i)+y
          wzmass(i)=wzmass(i)+z
          wtmass=wtmass+wamass(i,j)
          xmass=xmass+x
          ymass=ymass+y
          zmass=zmass+z
          tmass=tmass+wamass(i,j)
        enddo
        wxmass(i)=wxmass(i)/wtmass
        wymass(i)=wymass(i)/wtmass
        wzmass(i)=wzmass(i)/wtmass
      enddo
      xmass=xmass/tmass
      ymass=ymass/tmass
      zmass=zmass/tmass
      write(10,'(a17,3f10.5)') 'xmass ymass zmass',xmass,ymass,zmass
C  
C     preperation of intra molecular bond vector 
C
      do i=1,nseq
        do j=1,nwbplz(i)
          bvecx(i,j)=wx(i,(nwbon1(i,j)))-wx(i,(nwbon2(i,j)))
          bvecy(i,j)=wy(i,(nwbon1(i,j)))-wy(i,(nwbon2(i,j)))
          bvecz(i,j)=wz(i,(nwbon1(i,j)))-wz(i,(nwbon2(i,j)))
          bonlen(i,J)=dsqrt(bvecx(i,j)**2+bvecy(i,j)**2+bvecz(i,j)**2)
          bvecx(i,j)=bvecx(i,j)/bonlen(i,J)
          bvecy(i,j)=bvecy(i,j)/bonlen(i,J)
          bvecz(i,j)=bvecz(i,j)/bonlen(i,J)
        enddo
      enddo
      return
      end
************************************************************************
      subroutine predis
C     
C     preparation of inter- & intramolecular atom-atom distances
C
      IMPLICIT REAL*8(A-H,O-Z)
      parameter (k50=250,k300=300,k3000=5000,k100=30)
      parameter(ks30=30,k200=300)
      parameter(kr50=50,ka200=300,k4=4,k21=21)
      character*4 resi,dummy
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/resid/ numres,resi(kr50)
      common/topol1/ tcrg(kr50),crg(kr50,ka200),
     &  atmplz(kr50,ka200,k4),bonplz(kr50,ka200,k4),natmplz(kr50),
     &  nbonplz(kr50),nbon1(kr50,ka200),nbon2(kr50,ka200),
     &  atmass(kr50,ka200),
     &  xrmin(kr50,ka200),xemin(kr50,ka200),mbon1(kr50,ka200),
     &  mbon2(kr50,ka200),boaplz(kr50,ka200,4),mbonplz(kr50),
     &  mseg(kr50,ka200)
      common/topol2/
     &  nlnk(kr50),nl12(kr50,ka200,ka200),nl13(kr50,ka200,ka200),
     &  nl14(kr50,ka200,ka200),scac14(kr50),scav14(kr50),
     &  facp12(kr50),facp13(kr50),facp14(kr50),facp15(kr50),
     &  facc12(kr50),facc13(kr50),facc14(kr50),facc15(kr50),ibonc(kr50),
     &  bchg(kr50,ka200),irest(kr50),psfac(kr50),factt(kr50)
      common/espot/ cx(k300),cy(k300),cz(k300),sx(k3000),sy(k3000),
     &  sz(k3000),spot(k3000),gcrg(k300),icoor,isurf,
     &  pop(k3000),gpop(k3000)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
C
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      data au/627.509541D0/,bohr/0.52917706D0/
      do i=1,nseq
        do j=1,nseq
          if(i.le.j)then
            n1=nwaplz(i)+nwbplz(i)+mwbplz(i)
            n2=nwaplz(j)+nwbplz(j)+mwbplz(j)
            do k=1,n1
              do l=1,n2
                rx(i,j,k,l)=wx(i,k)-wx(j,l) 
                ry(i,j,k,l)=wy(i,k)-wy(j,l) 
                rz(i,j,k,l)=wz(i,k)-wz(j,l) 
                rdis(i,j,k,l)=dsqrt(rx(i,j,k,l)**2
     &              +ry(i,j,k,l)**2+rz(i,j,k,l)**2)
                rx(j,i,l,k)=-rx(i,j,k,l) 
                ry(j,i,l,k)=-ry(i,j,k,l) 
                rz(j,i,l,k)=-rz(i,j,k,l) 
                rdis(j,i,l,k)=rdis(i,j,k,l)
              enddo
            enddo
          endif
        enddo
      enddo
      write(10,*) '------ input data ---------------------------------'
      write(10,*) '         x         y         z     p charge  plz alpha
     &    beta     gamma     delta     rmin       emin      segment'
      do i=1,nseq
        write(10,*) '%% nwaplz nwbplz mwbplz %%',
     &    i,nwaplz(i),nwbplz(i),mwbplz(i)
        nn=nwaplz(i)+nwbplz(i)+mwbplz(i)
        do j=1,nn
          write(10,'(2i3,3f10.5,5f10.3,2f10.5,i3)') 
     &      i,j,wx(i,j),wy(i,j),wz(i,j),
     &      wcrg(i,j),wplz(i,j,1),wplz(i,j,2),wplz(i,j,3),wplz(i,j,4),
     &      wrmin(i,J)*bohr,wemin(i,j)*au,nseg(i,j)
        enddo
      enddo
      write(10,*) '---------------------------------------------------'
      RETURN
      END
************************************************************************
      subroutine vdwene
C     
C     inter molecular vdw energy
C
      parameter (k50=250,k300=300,k3000=5000,k100=30)
      parameter(ks30=30,k200=300)
      parameter(kr50=50,ka200=300,k4=4,k12=12)
      implicit real*8 (a-h,o-z)
      character*4 resi,dummy
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/resid/ numres,resi(kr50)
      common/topol1/ tcrg(kr50),crg(kr50,ka200),
     &  atmplz(kr50,ka200,k4),bonplz(kr50,ka200,k4),natmplz(kr50),
     &  nbonplz(kr50),nbon1(kr50,ka200),nbon2(kr50,ka200),
     &  atmass(kr50,ka200),
     &  xrmin(kr50,ka200),xemin(kr50,ka200),mbon1(kr50,ka200),
     &  mbon2(kr50,ka200),boaplz(kr50,ka200,4),mbonplz(kr50),
     &  mseg(kr50,ka200)
      common/topol2/
     &  nlnk(kr50),nl12(kr50,ka200,ka200),nl13(kr50,ka200,ka200),
     &  nl14(kr50,ka200,ka200),scac14(kr50),scav14(kr50),
     &  facp12(kr50),facp13(kr50),facp14(kr50),facp15(kr50),
     &  facc12(kr50),facc13(kr50),facc14(kr50),facc15(kr50),ibonc(kr50),
     &  bchg(kr50,ka200),irest(kr50),psfac(kr50),factt(kr50)
      common/espot/ cx(k300),cy(k300),cz(k300),sx(k3000),sy(k3000),
     &  sz(k3000),spot(k3000),gcrg(k300),icoor,isurf,
     &  pop(k3000),gpop(k3000)
C
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/energy/ dees,es(k100,k100),devdw,vdw(k100,k100),
     &  plztot,plz(k100,k100),totene,tot(k100,k100) 
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
      data zero/0.0d0/,au/627.509541D0/,bohr/0.52917706D0/
*     write(10,*) 'enter vdwene'
*     write(10,'(10f10.5)') wrmin
*     write(10,'(10f10.5)') wemin
      vdwtot=zero
      do i=1,k100
        do j=1,k100
          vdw(i,j)=zero
        enddo
      enddo
      do i=1,nseq
        do j=i+1,nseq
          do k=1,nwaplz(i)
            do l=1,nwaplz(j)
              r=rdis(i,j,k,l)
              rmin12=wrmin(i,k)+wrmin(j,l)
              emin12=dsqrt(wemin(i,k)*wemin(j,l))
              r=rmin12/r
              r3=r*r*r
              r6=r3*r3
              r12=r6*r6
              v=emin12*(r12-2.0d0*r6)
              vdw(i,j)=vdw(i,j)+v
            enddo
          enddo
          vdwtot=vdwtot+vdw(i,j)
        enddo
      enddo
      write(10,*) 'inter vdwtot   kcal/mol'
      write(10,'(f10.5)') vdwtot*au
C
C     intramolecular interaction
C
      if(nwlnk(1).ne.0)then
        vdwin15=zero
        vdwin14=zero
        do i=1,nseq
          do j=1,nseq
            if(i.eq.j)then
              do k=1,nwaplz(i)
                do l=k+1,nwaplz(j)
                  r=rdis(i,j,k,l)
                  rmin12=wrmin(i,k)+wrmin(j,l)
                  emin12=dsqrt(wemin(i,k)*wemin(j,l))
                  r=rmin12/r
                  r3=r*r*r
                  r6=r3*r3
                  r12=r6*r6
                  v=emin12*(r12-2.0d0*r6)*scav(i,k,l)
                  vdwin15=vdwin15+v
                  v=emin12*(r12-2.0d0*r6)*sv14(i,k,l)
                  vdwin14=vdwin14+v
                enddo
              enddo
            endif
          enddo
        enddo
        write(10,*) 'inter vdw, intra 1-5..., intra 1-4 vdw 
     & and sum vdw  kcal/mol'
        devdw=vdwtot+vdwin15+vdwin14
        write(10,'(a10,4f15.5)') ' VDWcompo ', 
     &  vdwtot*au,vdwin15*au,vdwin14*au,devdw*au
      endif
*     write(10,*) 'end of subr.vdwene'
      return
      end
************************************************************************
      subroutine elesta
C     
C     inter & intramolecular ELEctroSTAtic energy
C
      parameter (k50=250,k300=300,k3000=5000,k100=30)
      parameter(kr50=50,ka200=300,k21=21,k4=4)
      parameter(ks30=30,k200=300)
      implicit real*8 (a-h,o-z)
      character*4 resi,dummy
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/resid/ numres,resi(kr50)
      common/topol1/ tcrg(kr50),crg(kr50,ka200),
     &  atmplz(kr50,ka200,k4),bonplz(kr50,ka200,k4),natmplz(kr50),
     &  nbonplz(kr50),nbon1(kr50,ka200),nbon2(kr50,ka200),
     &  atmass(kr50,ka200),
     &  xrmin(kr50,ka200),xemin(kr50,ka200),mbon1(kr50,ka200),
     &  mbon2(kr50,ka200),boaplz(kr50,ka200,4),mbonplz(kr50),
     &  mseg(kr50,ka200)
      common/topol2/
     &  nlnk(kr50),nl12(kr50,ka200,ka200),nl13(kr50,ka200,ka200),
     &  nl14(kr50,ka200,ka200),scac14(kr50),scav14(kr50),
     &  facp12(kr50),facp13(kr50),facp14(kr50),facp15(kr50),
     &  facc12(kr50),facc13(kr50),facc14(kr50),facc15(kr50),ibonc(kr50),
     &  bchg(kr50,ka200),irest(kr50),psfac(kr50),factt(kr50)
      common/espot/ cx(k300),cy(k300),cz(k300),sx(k3000),sy(k3000),
     &  sz(k3000),spot(k3000),gcrg(k300),icoor,isurf,
     &  pop(k3000),gpop(k3000)
C
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/energy/ dees,es(k100,k100),devdw,vdw(k100,k100),
     &  plztot,plz(k100,k100),totene,tot(k100,k100) 
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
      data au/627.509541D0/,zero/0.0d0/
      estot=zero
      do i=1,k100
        do j=1,k100
          es(i,j)=zero
        enddo
      enddo
C
C     intermolecular interaction
C
      do i=1,nseq
        do j=i+1,nseq
          do k=1,nwaplz(i)+iwbonc(i)
*         do k=1,nwaplz(i)
            do l=1,nwaplz(j)+iwbonc(j)
*           do l=1,nwaplz(j)
              r=rdis(i,j,k,l)
              es(i,j)=es(i,j)+wcrg(i,k)*wcrg(j,l)/r
            enddo
          enddo
          estot=estot+es(i,j)
        enddo
      enddo
      write(10,*) 'estot   kcal/mol'
      write(10,'(f14.5)') estot*au
C
C     intramolecular interaction 
C
      if(nwlnk(1).ne.0)then
        esin15=zero
        esin14=zero
        do i=1,nseq
          do j=1,nseq
            if(i.eq.j)then
            do k=1,nwaplz(i)+iwbonc(i)
*           do k=1,nwaplz(i)
              do l=k+1,nwaplz(j)+iwbonc(j)
*             do l=k+1,nwaplz(j)
                r=rdis(i,j,k,l)
                esin15=esin15+scac(i,k,l)*wcrg(i,k)*wcrg(j,l)/r
                esin14=esin14+se14(i,k,l)*wcrg(i,k)*wcrg(j,l)/r
              enddo
            enddo
            endif
          enddo
        enddo
C
        write(10,*) 'inter es, intra 1-5..., intra 1-4 es 
     & and sum ES  kcal/mol'
       dees=estot+esin15+esin14
        write(10,'(a10,4f15.5)') ' EScompo  ', 
     &    estot*au,esin15*au,esin14*au,dees*au
      endif
*     write(10,*) 'end of subr.elesta'
      return
      end
************************************************************************
      subroutine cnsdip
C     
C     DIPole moments for constant charges tcdipx y z 
C
C     DIPole moments for gaussian charges gdipx y z 
C
      parameter (k50=250,k300=300,k3000=5000,k100=30)
      parameter(kr50=50,ka200=300,k21=21,k4=4)
      parameter(ks30=30,k200=300)
      implicit real*8 (a-h,o-z)
      character*4 resi,dummy
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/resid/ numres,resi(kr50)
      common/topol1/ tcrg(kr50),crg(kr50,ka200),
     &  atmplz(kr50,ka200,k4),bonplz(kr50,ka200,k4),natmplz(kr50),
     &  nbonplz(kr50),nbon1(kr50,ka200),nbon2(kr50,ka200),
     &  atmass(kr50,ka200),
     &  xrmin(kr50,ka200),xemin(kr50,ka200),mbon1(kr50,ka200),
     &  mbon2(kr50,ka200),boaplz(kr50,ka200,4),mbonplz(kr50),
     &  mseg(kr50,ka200)
      common/topol2/
     &  nlnk(kr50),nl12(kr50,ka200,ka200),nl13(kr50,ka200,ka200),
     &  nl14(kr50,ka200,ka200),scac14(kr50),scav14(kr50),
     &  facp12(kr50),facp13(kr50),facp14(kr50),facp15(kr50),
     &  facc12(kr50),facc13(kr50),facc14(kr50),facc15(kr50),ibonc(kr50),
     &  bchg(kr50,ka200),irest(kr50),psfac(kr50),factt(kr50)
      common/espot/ cx(k300),cy(k300),cz(k300),sx(k3000),sy(k3000),
     &  sz(k3000),spot(k3000),gcrg(k300),icoor,isurf,
     &  pop(k3000),gpop(k3000)
C
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/energy/ dees,es(k100,k100),devdw,vdw(k100,k100),
     &  plztot,plz(k100,k100),totene,tot(k100,k100) 
      common/dipole/tcdipx,tcdipy,tcdipz,tmdipx,tmdipy,tmdipz,
     &  gdipx,gdipy,gdipz,gdicx,gdicy,gdicz
      dimension dipx(k50),dipy(k50),dipz(k50)
      dimension dimx(k50),dimy(k50),dimz(k50)
      dimension xxtcdipx(k300),xxtcdipy(k300),xxtcdipz(k300)
      data au/627.509541D0/
      data debye/2.54158059d0/
      tcdipx=0.0d0
      tcdipy=0.0d0
      tcdipz=0.0d0
      tmdipx=0.0d0
      tmdipy=0.0d0
      tmdipz=0.0d0
      ix=0
      do i=1,nseq
        dipx(i)=0.0d0
        dipy(i)=0.0d0
        dipz(i)=0.0d0
        dimx(i)=0.0d0
        dimy(i)=0.0d0
        dimz(i)=0.0d0
        na=nwaplz(i)
        do j=1,na+iwbonc(i)
          ix=ix+1
          dipx(i)=dipx(i)+(wx(i,j)-wxmass(i))*wcrg(i,j)
          dipy(i)=dipy(i)+(wy(i,j)-wymass(i))*wcrg(i,j)
          dipz(i)=dipz(i)+(wz(i,j)-wzmass(i))*wcrg(i,j)
          dimx(i)=dimx(i)+wx(i,j)*wcrg(i,j)
          dimy(i)=dimy(i)+wy(i,j)*wcrg(i,j)
          dimz(i)=dimz(i)+wz(i,j)*wcrg(i,j)
          tcdipx=tcdipx+(wx(i,j)-xmass)*wcrg(i,j)
          tcdipy=tcdipy+(wy(i,j)-ymass)*wcrg(i,j)
          tcdipz=tcdipz+(wz(i,j)-zmass)*wcrg(i,j)
          xxtcdipx(ix)=tcdipx
          xxtcdipy(ix)=tcdipy
          xxtcdipz(ix)=tcdipz
          tmdipx=tmdipx+wx(i,j)*wcrg(i,j)
          tmdipy=tmdipy+wy(i,j)*wcrg(i,j)
          tmdipz=tmdipz+wz(i,j)*wcrg(i,j)
C         write(6,'(a6,3i4,2f6.2,a2,f6.2)')'&cc&gc',i,j,ix,
C    &    wcrg(i,j),gcrg(ix),'! ',wcrg(i,j)-gcrg(ix)
C         write(6,'(a7,6f8.3)')'&codxyz',
C    &    wx(i,j),wy(i,j),wz(i,j),cx(ix),cy(ix),cz(ix)
        enddo
      enddo
      gdipx=0.0d0
      gdipy=0.0d0
      gdipz=0.0d0
      gdicx=0.0d0
      gdicy=0.0d0
      gdicz=0.0d0
      do i=1,icoor
        gdipx=gdipx+(cx(i)-xmass)*gcrg(i)
        gdipy=gdipy+(cy(i)-ymass)*gcrg(i)
        gdipz=gdipz+(cz(i)-zmass)*gcrg(i)
        gdicx=gdicx+cx(i)*gcrg(i)
        gdicy=gdicy+cy(i)*gcrg(i)
        gdicz=gdicz+cz(i)*gcrg(i)
C       write(6,'(a6,i4,6f6.2,a2,3f6.2)') '&cd&gd',i,
C    &  xxtcdipx(i),xxtcdipy(i),xxtcdipz(i),
C    &  gdipx,gdipy,gdipz,'! ',
C    &  xxtcdipx(i)-gdipx,xxtcdipy(i)-gdipy,xxtcdipz(i)-gdipz
      enddo
      return
      end
************************************************************************
      subroutine fldcrg
C
C     Field From Point Fractional Charge
C
      parameter (k50=250,k300=300,k3000=5000,k100=30)
      parameter(kr50=50,ka200=300,k4=4,k21=21)
      parameter(ks30=30,k200=300)
      implicit real*8 (a-h,o-z)
      character*4 resi,dummy
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/resid/ numres,resi(kr50)
      common/topol1/ tcrg(kr50),crg(kr50,ka200),
     &  atmplz(kr50,ka200,k4),bonplz(kr50,ka200,k4),natmplz(kr50),
     &  nbonplz(kr50),nbon1(kr50,ka200),nbon2(kr50,ka200),
     &  atmass(kr50,ka200),
     &  xrmin(kr50,ka200),xemin(kr50,ka200),mbon1(kr50,ka200),
     &  mbon2(kr50,ka200),boaplz(kr50,ka200,4),mbonplz(kr50),
     &  mseg(kr50,ka200)
      common/topol2/
     &  nlnk(kr50),nl12(kr50,ka200,ka200),nl13(kr50,ka200,ka200),
     &  nl14(kr50,ka200,ka200),scac14(kr50),scav14(kr50),
     &  facp12(kr50),facp13(kr50),facp14(kr50),facp15(kr50),
     &  facc12(kr50),facc13(kr50),facc14(kr50),facc15(kr50),ibonc(kr50),
     &  bchg(kr50,ka200),irest(kr50),psfac(kr50),factt(kr50)
      common/espot/ cx(k300),cy(k300),cz(k300),sx(k3000),sy(k3000),
     &  sz(k3000),spot(k3000),gcrg(k300),icoor,isurf,
     &  pop(k3000),gpop(k3000)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
C
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/energy/ dees,es(k100,k100),devdw,vdw(k100,k100),
     &  plztot,plz(k100,k100),totene,tot(k100,k100) 
      common/fldcharg/
     &    cafldx(k100,k50),cafldy(k100,k50),
     &    cafldz(k100,k50)
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
C
C     cafldx,y,z : field on atom  
C
      do i=1,nseq
        do j=1,nwaplz(i)+nwbplz(i)+mwbplz(i)
          cafldx(i,j)=0.0d0
          cafldy(i,j)=0.0d0
          cafldz(i,j)=0.0d0
        enddo
      enddo
      do i=1,nseq
        do j=1,nwaplz(i)+nwbplz(i)+mwbplz(i)
          do k=1,nseq
            if(i.ne.k)then
              do l=1,nwaplz(k)+nwbplz(k)+mwbplz(k)
                xl=rx(i,k,j,l) 
                yl=ry(i,k,j,l) 
                zl=rz(i,k,j,l) 
                r=rdis(i,k,j,l)
                r=r**3
                cafldx(i,j)=cafldx(i,j)+wcrg(k,l)*xl/r
                cafldy(i,j)=cafldy(i,j)+wcrg(k,l)*yl/r
                cafldz(i,j)=cafldz(i,j)+wcrg(k,l)*zl/r
              enddo
            endif
          enddo
        enddo
      enddo
C
C     intramolecular charge field only for atomic polz
C
      if(nwlnk(1).ne.0)then
        do i=1,nseq
          do j=1,nwaplz(i)
            do k=1,nseq
              if(i.eq.k)then
                do l=1,nwaplz(k)
                  if(j.ne.l)then
                    xl=rx(i,k,j,l) 
                    yl=ry(i,k,j,l) 
                    zl=rz(i,k,j,l) 
                    r=rdis(i,k,j,l)
                    r=r**3
                    cafldx(i,j)=cafldx(i,j)+facc(i,j,l)*wcrg(k,l)*xl/r
                    cafldy(i,j)=cafldy(i,j)+facc(i,j,l)*wcrg(k,l)*yl/r
                    cafldz(i,j)=cafldz(i,j)+facc(i,j,l)*wcrg(k,l)*zl/r
                  endif
                enddo
              endif
            enddo
          enddo
        enddo
      endif
      return
      end
************************************************************************
      subroutine idfield(i,j,k,l,dx,dy,dz,fx,fy,fz)
C
C     calculation of induced dipole field
C     Induced delta charge model
C
C     dx,dy,dz: induced dipole j,l
C     fx,fy,fz: field i,k from induced dipole j,l
C
      parameter(ks30=30,k200=300)
      parameter (k50=250,k100=30,k500=500,k4=4)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      data d05/0.5d0/,one/1.0d0/
      r=dsqrt(dx**2+dy**2+dz**2)
      na=nwaplz(j)
      nb=nwbplz(j)
      if(l.le.na.or.l.gt.na+nb)then
        q=r/one
      else
        q=r/bonlen(j,l-na)
      endif
      xh=d05*dx/q
      yh=d05*dy/q
      zh=d05*dz/q
      xp=wx(j,l)+xh
      yp=wy(j,l)+yh
      zp=wz(j,l)+zh
      rp=dsqrt((wx(i,k)-xp)**2+(wy(i,k)-yp)**2+(wz(i,k)-zp)**2)
      rp3=rp*rp*rp
      fx=q*(wx(i,k)-xp)/rp3
      fy=q*(wy(i,k)-yp)/rp3
      fz=q*(wz(i,k)-zp)/rp3
      xm=wx(j,l)-xh
      ym=wy(j,l)-yh
      zm=wz(j,l)-zh
      rm=dsqrt((wx(i,k)-xm)**2+(wy(i,k)-ym)**2+(wz(i,k)-zm)**2)
      rm3=rm*rm*rm
      fx=fx-q*(wx(i,k)-xm)/rm3
      fy=fy-q*(wy(i,k)-ym)/rm3
      fz=fz-q*(wz(i,k)-zm)/rm3
      return
      end
************************************************************************
      subroutine idfielp(ii,jj,kk,ll)
C
C     calculation of intermolecular induced dipole field
C     induced point dipole
C
C     i1.ne.jj
C
C     i block molecule     dx,dy,dz: induced dipole j,l
C     j block molecule     fx,fy,fz: field i,k from induced dipole j,l
C
      parameter (k50=250,k100=30,k500=500,k4=4)
      parameter(ks30=30,k200=300)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
cc 
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
      common/tens/
     &  txx(k50,k50),txy(k50,k50),
     &  txz(k50,k50),
     &  tyx(k50,k50),tyy(k50,k50),
     &  tyz(k50,k50),
     &  tzx(k50,k50),tzy(k50,k50),
     &  tzz(k50,k50)
cc
      data one/1.0d0/,three/3.0d0/
      ni=0
      do i=1,ii-1
        ni=ni+nwaplz(i)
      enddo
      nj=0
      do j=1,jj-1
        nj=nj+nwaplz(j)
      enddo
*     write(6,*)'#fieldp ii jj kk ll ni nj',ii,jj,kk,ll,ni,nj
      do k=1,kk
        do l=1,ll
          ax=wx(ii,k)-wx(jj,l)
          ay=wy(ii,k)-wy(jj,l)
          az=wz(ii,k)-wz(jj,l)
          r1=dsqrt(ax**2+ay**2+az**2)
          r2=r1*r1
          r3=r1*r2
          r5=r2*r3
          s1=one/r3
          s2=three/r5
cc
          txx(k+ni,l+nj)=-(s2*ax*ax-s1)*facp(ii,k,l)
          txy(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
          txz(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
          tyx(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
          tyy(k+ni,l+nj)=-(s2*ay*ay-s1)*facp(ii,k,l)
          tyz(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
          tzx(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
          tzy(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
          tzz(k+ni,l+nj)=-(s2*az*az-s1)*facp(ii,k,l)
        enddo
      enddo
cc
*     write(6,*) '#fieldp txx ...nii njj',nii,njj
*     do k=1,kk
*       do l=1,ll
*          write(6,'(2i5,9f10.5)') k+ni,l+nj,
*    &      txx(k+ni,l+nj),txy(k+ni,l+nj),txz(k+ni,l+nj),
*    &      tyx(k+ni,l+nj),tyy(k+ni,l+nj),tyz(k+ni,l+nj),
*    &      tzx(k+ni,l+nj),tzy(k+ni,l+nj),tzz(k+ni,l+nj)
*       enddo
*     enddo
cc
cc      
C     fx=((s2*ax*ax-s1)*dx+
C    &        s2*ax*ay*dy+
C    &        s2*ax*az*dz)
C     fy=(s2*ax*ay*dx+
C    &       (s2*ay*ay-s1)*dy+
C    &        s2*ay*az*dz)
C     fz=(s2*ax*az*dx+
C    &        s2*ay*az*dy+
C    &       (s2*az*az-s1)*dz)
      return
      end
************************************************************************
      subroutine imdflds(ii,jj,kk,ll)
C
C     calculation of intramolecular induced dipole field
C     Simple Scaling model
C
C     irest 13
C
C     ii.eq.jj
C
C     dx,dy,dz: induced dipole j,l
C     fx,fy,fz: field i,k from induced dipole j,l
C
      parameter(ks30=30,k200=300)
      parameter (k50=250,k100=30,k500=500,k4=4)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
cc
      common/tens/
     &  txx(k50,k50),txy(k50,k50),txz(k50,k50),
     &  tyx(k50,k50),tyy(k50,k50),tyz(k50,k50),
     &  tzx(k50,k50),tzy(k50,k50),tzz(k50,k50)
cc
      data zero/0.0d0/,one/1.0d0/,three/3.0d0/
      ni=0
      do i=1,ii-1
        ni=ni+nwaplz(i)
      enddo
      nj=0
      do j=1,jj-1
        nj=nj+nwaplz(j)
      enddo
*     write(6,*) '#imdflds ii jj kk ll ni nj ', ii,jj,kk,ll,ni,nj
      do k=1,kk
        do l=1,ll
          if(k.ne.l.or.ii.ne.jj)then
            ax=wx(ii,k)-wx(jj,l)
            ay=wy(ii,k)-wy(jj,l)
            az=wz(ii,k)-wz(jj,l)
            r1=dsqrt(ax**2+ay**2+az**2)
            r2=r1*r1
            r3=r1*r2
            r5=r2*r3
            s1=one/r3
            s2=three/r5
cc
            txx(k+ni,l+nj)=-(s2*ax*ax-s1)*facp(ii,k,l)
            txy(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
            txz(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
            tyx(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
            tyy(k+ni,l+nj)=-(s2*ay*ay-s1)*facp(ii,k,l)
            tyz(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
            tzx(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
            tzy(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
            tzz(k+ni,l+nj)=-(s2*az*az-s1)*facp(ii,k,l)
          else
            alp=wplz(jj,l,1)
CSN 2022
            if(alp.eq.zero)alp=0.00000001d0
CSN 2022
            txx(k+ni,l+nj)=one/alp
            txy(k+ni,l+nj)=zero
            txz(k+ni,l+nj)=zero
            tyx(k+ni,l+nj)=zero
            tyy(k+ni,l+nj)=one/alp
            tyz(k+ni,l+nj)=zero
            tzx(k+ni,l+nj)=zero
            tzy(k+ni,l+nj)=zero
            tzz(k+ni,l+nj)=one/alp
          endif 
        enddo
      enddo
cc
      return
      end
************************************************************************
      subroutine imdfldli(ii,jj,kk,ll)
C
C     calculation of intramolecular induced dipole field
C     
C     Thole Linear model  van Duijnen a=1.7823  original Thole a=1.662
C     with simple scaling
C 
C     irest 173          2020/12/18
C
C     dx,dy,dz: induced dipole j,l
C     fx,fy,fz: field i,k from induced dipole j,l
C
      parameter (ks30=30,k200=300)
      parameter (k50=250,k100=30,k500=500,k4=4)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
cc
      common/tens/
     &  txx(k50,k50),txy(k50,k50),txz(k50,k50),
     &  tyx(k50,k50),tyy(k50,k50),tyz(k50,k50),
     &  tzx(k50,k50),tzy(k50,k50),tzz(k50,k50)
cc
      data d05/0.5d0/,one/1.0d0/,three/3.0d0/,bohr/0.52917706D0/,
     &     six/6.0d0/,two/2.0d0/,four/4.0d0/,zero/0.0d0/
      ni=0
      do i=1,ii-1
        ni=ni+nwaplz(i)
      enddo
      nj=0
      do j=1,jj-1
        nj=nj+nwaplz(j)
      enddo
*     write(6,*) '#imdflds ii jj kk ll ni nj ', ii,jj,kk,ll,ni,nj
      do k=1,kk
        do l=1,ll
          if(k.ne.l.or.ii.ne.jj)then
            ax=wx(ii,k)-wx(jj,l)
            ay=wy(ii,k)-wy(jj,l)
            az=wz(ii,k)-wz(jj,l)
            r1=dsqrt(ax**2+ay**2+az**2)
            r2=r1*r1
            r3=r1*r2
            r5=r2*r3
C   
C           avoid /zero      
C   
            aa=(wplz(ii,k,1)*wplz(jj,l,1))
            if(aa.gt.0.0001d0)then
              aa=aa**(one/six)
              s=wfacp(ii)*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
cc
            txx(k+ni,l+nj)=-(s2*ax*ax-s1)*facp(ii,k,l)
            txy(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
            txz(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
            tyx(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
            tyy(k+ni,l+nj)=-(s2*ay*ay-s1)*facp(ii,k,l)
            tyz(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
            tzx(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
            tzy(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
            tzz(k+ni,l+nj)=-(s2*az*az-s1)*facp(ii,k,l)

          else
            alp=wplz(jj,l,1)
            txx(k+ni,l+nj)=one/alp
            txy(k+ni,l+nj)=zero
            txz(k+ni,l+nj)=zero
            tyx(k+ni,l+nj)=zero
            tyy(k+ni,l+nj)=one/alp
            tyz(k+ni,l+nj)=zero
            tzx(k+ni,l+nj)=zero
            tzy(k+ni,l+nj)=zero
            tzz(k+ni,l+nj)=one/alp
          endif
        enddo
      enddo
      return
      end
************************************************************************
      subroutine imdfldt(ii,jj,kk,ll)
C
C     calculation of intramolecular induced dipole field
C     Thole's model  van Duijnen Exp type   1.9088  original Thole 2.089
C
C     irest 14
C
C     dx,dy,dz: induced dipole j,l
C     fx,fy,fz: field i,k from induced dipole j,l
C
      parameter (ks30=30,k200=300)
      parameter (k50=250,k100=30,k500=500,k4=4)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
cc
      common/tens/
     &  txx(k50,k50),txy(k50,k50),txz(k50,k50),
     &  tyx(k50,k50),tyy(k50,k50),tyz(k50,k50),
     &  tzx(k50,k50),tzy(k50,k50),tzz(k50,k50)
cc
      data d05/0.5d0/,one/1.0d0/,three/3.0d0/,bohr/0.52917706D0/,
     &     six/6.0d0/,two/2.0d0/,zero/0.0d0/
      ni=0
      do i=1,ii-1
        ni=ni+nwaplz(i)
      enddo
      nj=0
      do j=1,jj-1
        nj=nj+nwaplz(j)
      enddo
*     write(6,*) '#imdflds ii jj kk ll ni nj ', ii,jj,kk,ll,ni,nj
      do k=1,kk
        do l=1,ll
          if(k.ne.l.or.ii.ne.jj)then
            ax=wx(ii,k)-wx(jj,l)
            ay=wy(ii,k)-wy(jj,l)
            az=wz(ii,k)-wz(jj,l)
            r1=dsqrt(ax**2+ay**2+az**2)
            r2=r1*r1
            r3=r1*r2
            r5=r2*r3
C   
C           avoid /zero      check aa bohr  *bohr can be remove
C   
            aa=(wplz(ii,k,1)*wplz(jj,l,1))
            if(aa.gt.0.0001d0)then
              aa=aa**(one/six)
              rr=r1*bohr
              u=rr/(aa*bohr)
              ar1=wfacp(ii)*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
cc
            txx(k+ni,l+nj)=-(s2*ax*ax-s1)
            txy(k+ni,l+nj)=-(s2*ax*ay)
            txz(k+ni,l+nj)=-(s2*ax*az)
            tyx(k+ni,l+nj)=-(s2*ax*ay)
            tyy(k+ni,l+nj)=-(s2*ay*ay-s1)
            tyz(k+ni,l+nj)=-(s2*ay*az)
            tzx(k+ni,l+nj)=-(s2*ax*az)
            tzy(k+ni,l+nj)=-(s2*ay*az)
            tzz(k+ni,l+nj)=-(s2*az*az-s1)
          else
            alp=wplz(jj,l,1)
            txx(k+ni,l+nj)=one/alp
            txy(k+ni,l+nj)=zero
            txz(k+ni,l+nj)=zero
            tyx(k+ni,l+nj)=zero
            tyy(k+ni,l+nj)=one/alp
            tyz(k+ni,l+nj)=zero
            tzx(k+ni,l+nj)=zero
            tzy(k+ni,l+nj)=zero
            tzz(k+ni,l+nj)=one/alp
          endif
        enddo
      enddo
C   
      return
      end
************************************************************************
      subroutine imdfldp(ii,jj,kk,ll)
C
C     calculation of intramolecular induced dipole field
C     Thole's model  Ponder Exp type   0.39  original Thole 0.572
C
C     irest 15
C
C     dx,dy,dz: induced dipole j,l
C     fx,fy,fz: field i,k from induced dipole j,l
C
      parameter(ks30=30,k200=300)
      parameter (k50=250,k100=30,k500=500,k4=4)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
cc
      common/tens/
     &  txx(k50,k50),txy(k50,k50),txz(k50,k50),
     &  tyx(k50,k50),tyy(k50,k50),tyz(k50,k50),
     &  tzx(k50,k50),tzy(k50,k50),tzz(k50,k50)
cc
      data d05/0.5d0/,one/1.0d0/,three/3.0d0/,bohr/0.52917706D0/,
     &     six/6.0d0/,two/2.0d0/,zero/0.0d0/
      ni=0
      do i=1,ii-1
        ni=ni+nwaplz(i)
      enddo
      nj=0
      do j=1,jj-1
        nj=nj+nwaplz(j)
      enddo
*     write(6,*) '#imdflds ii jj kk ll ni nj ', ii,jj,kk,ll,ni,nj
      do k=1,kk
        do l=1,ll
          if(k.ne.l.or.ii.ne.jj)then
            ax=wx(ii,k)-wx(jj,l)
            ay=wy(ii,k)-wy(jj,l)
            az=wz(ii,k)-wz(jj,l)
            r1=dsqrt(ax**2+ay**2+az**2)
            r2=r1*r1
            r3=r1*r2
            r5=r2*r3
C   
C           avoid /zero       *bohr can be remove
C   
            aa=(wplz(ii,k,1)*wplz(jj,l,1))
            if(aa.gt.0.0001d0)then
              aa=aa**(one/six)
              rr=r1*bohr
              u=rr/(aa*bohr)
              au3=wfacp(ii)*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
cc
            txx(k+ni,l+nj)=-(s2*ax*ax-s1)
            txy(k+ni,l+nj)=-(s2*ax*ay)
            txz(k+ni,l+nj)=-(s2*ax*az)
            tyx(k+ni,l+nj)=-(s2*ax*ay)
            tyy(k+ni,l+nj)=-(s2*ay*ay-s1)
            tyz(k+ni,l+nj)=-(s2*ay*az)
            tzx(k+ni,l+nj)=-(s2*ax*az)
            tzy(k+ni,l+nj)=-(s2*ay*az)
            tzz(k+ni,l+nj)=-(s2*az*az-s1)
          else
            alp=wplz(jj,l,1)
            txx(k+ni,l+nj)=one/alp
            txy(k+ni,l+nj)=zero
            txz(k+ni,l+nj)=zero
            tyx(k+ni,l+nj)=zero
            tyy(k+ni,l+nj)=one/alp
            tyz(k+ni,l+nj)=zero
            tzx(k+ni,l+nj)=zero
            tzy(k+ni,l+nj)=zero
            tzz(k+ni,l+nj)=one/alp
          endif
        enddo
      enddo
C
      return
      end
************************************************************************
      subroutine imdfldts(ii,jj,kk,ll)
C
C     calculation of intramolecular induced dipole field
C     Thole's model  van Duijnen Exp type   1.9088  original Thole 2.089
C     with simple scaling
C 
C     irest 143 
C
C     dx,dy,dz: induced dipole j,l
C     fx,fy,fz: field i,k from induced dipole j,l
C
      parameter (ks30=30,k200=300)
      parameter (k50=250,k100=30,k500=500,k4=4)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
cc
      common/tens/
     &  txx(k50,k50),txy(k50,k50),txz(k50,k50),
     &  tyx(k50,k50),tyy(k50,k50),tyz(k50,k50),
     &  tzx(k50,k50),tzy(k50,k50),tzz(k50,k50)
cc
      data d05/0.5d0/,one/1.0d0/,three/3.0d0/,bohr/0.52917706D0/,
     &     six/6.0d0/,two/2.0d0/,zero/0.0d0/
      ni=0
      do i=1,ii-1
        ni=ni+nwaplz(i)
      enddo
      nj=0
      do j=1,jj-1
        nj=nj+nwaplz(j)
      enddo
*     write(6,*) '#imdflds ii jj kk ll ni nj ', ii,jj,kk,ll,ni,nj
      do k=1,kk
        do l=1,ll
          if(k.ne.l.or.ii.ne.jj)then
            ax=wx(ii,k)-wx(jj,l)
            ay=wy(ii,k)-wy(jj,l)
            az=wz(ii,k)-wz(jj,l)
            r1=dsqrt(ax**2+ay**2+az**2)
            r2=r1*r1
            r3=r1*r2
            r5=r2*r3
C   
C           avoid /zero      2019.4.4       2020.8.5 *bohr can be remove
C   
            aa=(wplz(ii,k,1)*wplz(jj,l,1))
            if(aa.gt.0.0001d0)then
              aa=aa**(one/six)
              rr=r1*bohr
              u=rr/(aa*bohr)
              ar1=wfacp(i)*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
cc
            txx(k+ni,l+nj)=-(s2*ax*ax-s1)*facp(ii,k,l)
            txy(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
            txz(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
            tyx(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
            tyy(k+ni,l+nj)=-(s2*ay*ay-s1)*facp(ii,k,l)
            tyz(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
            tzx(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
            tzy(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
            tzz(k+ni,l+nj)=-(s2*az*az-s1)*facp(ii,k,l)
          else
            alp=wplz(jj,l,1)
            txx(k+ni,l+nj)=one/alp
            txy(k+ni,l+nj)=zero
            txz(k+ni,l+nj)=zero
            tyx(k+ni,l+nj)=zero
            tyy(k+ni,l+nj)=one/alp
            tyz(k+ni,l+nj)=zero
            tzx(k+ni,l+nj)=zero
            tzy(k+ni,l+nj)=zero
            tzz(k+ni,l+nj)=one/alp
          endif
        enddo
      enddo
C
      return
      end
************************************************************************
      subroutine imdfldms(ii,jj,kk,ll)
C
C     calculation of intramolecular induced dipole field
C     Masia's model  Gaussian density screening  Atom Type Beta=0.8456
C     Beta from J.Wang JTCT 2019 15 1146
C     with simple scaling
C 
C     irest 163      2019.8.30  2020.8.5
C
C     dx,dy,dz: induced dipole j,l
C     fx,fy,fz: field i,k from induced dipole j,l
C
      parameter (ks30=30,k200=300)
      parameter (k50=250,k100=30,k500=500,k4=4)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
cc
      common/tens/
     &  txx(k50,k50),txy(k50,k50),txz(k50,k50),
     &  tyx(k50,k50),tyy(k50,k50),tyz(k50,k50),
     &  tzx(k50,k50),tzy(k50,k50),tzz(k50,k50)
cc
      data d05/0.5d0/,one/1.0d0/,three/3.0d0/,bohr/0.52917706D0/,
     &     six/6.0d0/,two/2.0d0/,pa/3.1415926d0/,zero/0.0d0/
C 
      ni=0
      do i=1,ii-1
        ni=ni+nwaplz(i)
      enddo
      nj=0
      do j=1,jj-1
        nj=nj+nwaplz(j)
      enddo
*     write(6,*) '#imdflds ii jj kk ll ni nj ', ii,jj,kk,ll,ni,nj
      do k=1,kk
        do l=1,ll
          if(k.ne.l.or.ii.ne.jj)then
            co0=two/(three*dsqrt(two*pa))
            othr=one/three
            co1=co0**(-othr)
            co2=two/(dsqrt(pa))
            osix=one/six
C
            ax=wx(ii,k)-wx(jj,l)
            ay=wy(ii,k)-wy(jj,l)
            az=wz(ii,k)-wz(jj,l)
            r1=dsqrt(ax**2+ay**2+az**2)
            r2=r1*r1
            r3=r1*r2
            r5=r2*r3
C
C     2020.8.5 Gaussian model
C
            bi=wfacp(ii)*co1*(wplz(ii,k,1)**(-othr))
            bj=wfacp(jj)*co1*(wplz(jj,l,1)**(-othr))
            bb=dsqrt(bi**2+bj**2)
            if(bb.gt.0.0000001d0)then
              bij=bi*bj/bb
              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
            else
              s1=one/r3
              s2=three/r5
            endif
cc
            txx(k+ni,l+nj)=-(s2*ax*ax-s1)*facp(ii,k,l)
            txy(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
            txz(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
            tyx(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
            tyy(k+ni,l+nj)=-(s2*ay*ay-s1)*facp(ii,k,l)
            tyz(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
            tzx(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
            tzy(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
            tzz(k+ni,l+nj)=-(s2*az*az-s1)*facp(ii,k,l)
          else
            alp=wplz(jj,l,1)
            txx(k+ni,l+nj)=one/alp
            txy(k+ni,l+nj)=zero
            txz(k+ni,l+nj)=zero
            tyx(k+ni,l+nj)=zero
            tyy(k+ni,l+nj)=one/alp
            tyz(k+ni,l+nj)=zero
            tzx(k+ni,l+nj)=zero
            tzy(k+ni,l+nj)=zero
            tzz(k+ni,l+nj)=one/alp
          endif
        enddo
      enddo
C
      return
      end
************************************************************************
      subroutine imdfldps(ii,jj,kk,ll)
C
C     calculation of intramolecular induced dipole field
C     Thole's model  Ponder Exp type   0.39  original Thole 0.572
C     with simple scaling
C
C     irest 153
C
C     dx,dy,dz: induced dipole j,l
C     fx,fy,fz: field i,k from induced dipole j,l
C
      parameter(ks30=30,k200=300)
      parameter (k50=250,k100=30,k500=500,k4=4)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
cc
      common/tens/
     &  txx(k50,k50),txy(k50,k50),txz(k50,k50),
     &  tyx(k50,k50),tyy(k50,k50),tyz(k50,k50),
     &  tzx(k50,k50),tzy(k50,k50),tzz(k50,k50)
cc
      data d05/0.5d0/,one/1.0d0/,three/3.0d0/,bohr/0.52917706D0/,
     &     six/6.0d0/,two/2.0d0/,zero/0.0d0/
      ni=0
      do i=1,ii-1
        ni=ni+nwaplz(i)
      enddo
      nj=0
      do j=1,jj-1
        nj=nj+nwaplz(j)
      enddo
*     write(6,*) '#imdflds ii jj kk ll ni nj ', ii,jj,kk,ll,ni,nj
      do k=1,kk
        do l=1,ll
          if(k.ne.l.or.ii.ne.jj)then
            ax=wx(ii,k)-wx(jj,l)
            ay=wy(ii,k)-wy(jj,l)
            az=wz(ii,k)-wz(jj,l)
            r1=dsqrt(ax**2+ay**2+az**2)
            r2=r1*r1
            r3=r1*r2
            r5=r2*r3
C   
C           avoid /zero      2019.4.4
C   
            aa=(wplz(ii,k,1)*wplz(jj,l,1))
            if(aa.gt.0.0001d0)then
              aa=aa**(one/six)
              rr=r1*bohr
              u=rr/(aa*bohr)
              au3=wfacp(i)*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
cc
            txx(k+ni,l+nj)=-(s2*ax*ax-s1)*facp(ii,k,l)
            txy(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
            txz(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
            tyx(k+ni,l+nj)=-(s2*ax*ay)*facp(ii,k,l)
            tyy(k+ni,l+nj)=-(s2*ay*ay-s1)*facp(ii,k,l)
            tyz(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
            tzx(k+ni,l+nj)=-(s2*ax*az)*facp(ii,k,l)
            tzy(k+ni,l+nj)=-(s2*ay*az)*facp(ii,k,l)
            tzz(k+ni,l+nj)=-(s2*az*az-s1)*facp(ii,k,l)
          else
            alp=wplz(jj,l,1)
            txx(k+ni,l+nj)=one/alp
            txy(k+ni,l+nj)=zero
            txz(k+ni,l+nj)=zero
            tyx(k+ni,l+nj)=zero
            tyy(k+ni,l+nj)=one/alp
            tyz(k+ni,l+nj)=zero
            tzx(k+ni,l+nj)=zero
            tzy(k+ni,l+nj)=zero
            tzz(k+ni,l+nj)=one/alp
          endif
        enddo
      enddo
C
      return
      end
************************************************************************
      subroutine plzcyc
C
C     Iterative solution   max cycle is i50
C
C     Matrix inversion for induced dipole model (irest 3,4,5,13,14,15,143,153,163,173)
C
      parameter(ks30=30,k200=300)
      parameter(k50=250,k100=30,k500=500,k4=4,i50=500)
      parameter(i200=250)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
C
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/fldcharg/
     &    cafldx(k100,k50),cafldy(k100,k50),
     &    cafldz(k100,k50)
      common/dipole/tcdipx,tcdipy,tcdipz,tmdipx,tmdipy,tmdipz,
     &  gdipx,gdipy,gdipz,gdicx,gdicy,gdicz
      common/energy/ dees,es(k100,k100),devdw,vdw(k100,k100),
     &  plztot,plz(k100,k100),totene,tot(k100,k100) 
      common/inddip/atdipx(k100,k50),atdipy(k100,k50),atdipz(k100,k50),
     &  nnot
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
cc
      common/tens/
     &  txx(k50,k50),txy(k50,k50),
     &  txz(k50,k50),
     &  tyx(k50,k50),tyy(k50,k50),
     &  tyz(k50,k50),
     &  tzx(k50,k50),tzy(k50,k50),
     &  tzz(k50,k50)
cc
      dimension dafldx(k50),dafldy(k50),dafldz(k50),
     &          atdipox(k100,k50),atdipoy(k100,k50),atdipoz(k100,k50),
     &          xatdip(k50),yatdip(k50),zatdip(k50),
     &          vpolatm(k100)
      dimension xdip(i200),ydip(i200),zdip(i200),
     &          xe(i200),ye(i200),ze(i200)
      data zero/0.0d0/,two/2.0d0/
      data au/627.509541D0/
      data debye/2.54158059d0/
C
      numab=0
      do i=1,nseq
        do j=1,nwaplz(i)+nwbplz(i)+mwbplz(i)
          numab=numab+1
          atdipx(i,j)=zero
          atdipy(i,j)=zero
          atdipz(i,j)=zero
        enddo
      enddo
cc
C
C     Field From Induced Dipole
C
C     dafldx,y,z : field on atom  
C
      if(iwrest(1).gt.2)then
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        jd=0
        do i=1,nseq
          do j=1,nwaplz(i)
            jd=jd+1
            dafldx(jd)=cafldx(i,j)
            dafldy(jd)=cafldy(i,j)
            dafldz(jd)=cafldz(i,j)
          enddo
        enddo
        do i=1,nseq
C
C     intermolecular induced dipole field for irest=3,4,5,13,14,15,143,153,163,173
C   
C     Block Molecule model       imdfldt & imdfldp added
C   
          do j=1,nseq
              nni=nwaplz(i)
              nnj=nwaplz(j)
              if(iwrest(i).eq.24.or.iwrest(i).eq.13.or.
     &          iwrest(i).eq.3)then
                call imdflds(i,j,nni,nnj)
              else if(iwrest(i).eq.14)then
                call imdfldt(i,j,nni,nnj)
              else if(iwrest(i).eq.15)then
                call imdfldp(i,j,nni,nnj)
              else if(iwrest(i).eq.143)then
                call imdfldts(i,j,nni,nnj)
              else if(iwrest(i).eq.153)then
                call imdfldps(i,j,nni,nnj)
              else if(iwrest(i).eq.163)then
                call imdfldms(i,j,nni,nnj)
              else if(iwrest(i).eq.173)then
                call imdfldli(i,j,nni,nnj)
              else if(iwrest(i).eq.4)then
                call imdfldt(i,j,nni,nnj)
              else if(iwrest(i).eq.5)then
                call imdfldp(i,j,nni,nnj)
              endif
          enddo
        enddo 
        call matinv(dafldx,dafldy,dafldz,xdip,ydip,zdip,xe,ye,ze)
        ka=0 
        do i=1,nseq
           vpolatm(i)=zero
          do j=1,nwaplz(i)
            ka=ka+1
            atdipx(i,j)=xdip(ka)
            atdipy(i,j)=ydip(ka)
            atdipz(i,j)=zdip(ka)
            vpolatm(i)=vpolatm(i)+xe(ka)+ye(ka)+ze(ka)
          enddo
          vpolatm(i)=-vpolatm(i)/two
*         write(10,*)' # plzcyc ',i,vpolatm(i)*au
          enddo
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      else
        nnot=0
        do n=1,i50
          do i=1,nseq
            do j=1,nwaplz(i)+nwbplz(i)+mwbplz(i)
              dafldx(j)=cafldx(i,j)
              dafldy(j)=cafldy(i,j)
              dafldz(j)=cafldz(i,j)
            enddo
            do j=1,nseq
              if(i.ne.j)then
                do k=1,nwaplz(i)+nwbplz(i)+mwbplz(i)
                  do l=1,nwaplz(j)+nwbplz(j)+mwbplz(j)
                    dx=atdipx(j,l)
                    dy=atdipy(j,l)
                    dz=atdipz(j,l)
                    if(dx.eq.zero.and.dy.eq.zero.and.dz.eq.zero)then
                      fx=zero
                      fy=zero
                      fz=zero
                    else
                      call idfield(i,j,k,l,dx,dy,dz,fx,fy,fz)
                    endif
                    dafldx(k)=dafldx(k)+fx
                    dafldy(k)=dafldy(k)+fy
                    dafldz(k)=dafldz(k)+fz
                  enddo
                enddo
              endif
            enddo
            call indatm(i,dafldx,dafldy,dafldz,
     &                xatdip,yatdip,zatdip,vpolat)
            vpolatm(i)=vpolat
            do k=1,nwaplz(i)+nwbplz(i)+mwbplz(i)
              atdipox(i,k)=atdipx(i,k)
              atdipoy(i,k)=atdipy(i,k)
              atdipoz(i,k)=atdipz(i,k)
              atdipx(i,k)=xatdip(k)
              atdipy(i,k)=yatdip(k)
              atdipz(i,k)=zatdip(k)
            enddo
          enddo
cccc seq
          sumdip=zero
          do i=1,nseq
            do k=1,nwaplz(i)+nwbplz(i)+mwbplz(i)
              sumdip=sumdip+dsqrt((atdipx(i,k)-atdipox(i,k))**2
     &                         +(atdipy(i,k)-atdipoy(i,k))**2
     &                         +(atdipz(i,k)-atdipoz(i,k))**2)
            enddo
          enddo
          sumdip=sumdip/dfloat(numab)
          write(10,*) 'sumdip=',sumdip
C         if(sumdip.lt.0.004d0) then
          if(sumdip.lt.0.000001d0) then
            write(10,'(a20,i5)') ' converged !  cycle=',n
            go to 999
          else if(sumdip.gt.2.0d0) then
            write(10,'(a47,i5)')
     &      '!!!! Started to diverge! Stop! > 2.0 !!! cycle=',n
            nnot=400
            go to 999
          endif
        enddo
        write(10,*) ' !!!!!!! NOT converged !!!!!!!',n
        nnot=500
      endif
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  999 continue
      plztot=zero
      totdpx=zero
      totdpy=zero
      totdpz=zero
      write(10,*) ' plztot vpolatm in kcal/mol'
      do i=1,nseq
        plztot=plztot+vpolatm(i)
        do k=1,nwaplz(i)+nwbplz(i)+mwbplz(i)
          totdpx=totdpx+atdipx(i,k)
          totdpy=totdpy+atdipy(i,k)
          totdpz=totdpz+atdipz(i,k)
          atdipt=dsqrt(atdipx(i,k)**2+atdipy(i,k)**2+atdipz(i,k)**2)
        enddo
        write(10,'(i5,2f10.5)') i,vpolatm(i)*au
      enddo 
C
C     final induced dipole moment
C
      alldipx=tcdipx+totdpx
      alldipy=tcdipy+totdpy
      alldipz=tcdipz+totdpz
      alldimx=tmdipx+totdpx
      alldimy=tmdipy+totdpy
      alldimz=tmdipz+totdpz
      tcdip= dsqrt(tcdipx**2+tcdipy**2+tcdipz**2)
      gdip=  dsqrt(gdipx**2+gdipy**2+gdipz**2)
      alldip=dsqrt(alldipx**2+alldipy**2+alldipz**2)
      tmdip= dsqrt(tmdipx**2+tmdipy**2+tmdipz**2)
      gdic=  dsqrt(gdicx**2+gdicy**2+gdicz**2)
      alldim=dsqrt(alldimx**2+alldimy**2+alldimz**2)
      totdp= dsqrt(totdpx**2+totdpy**2+totdpz**2)
C
C     write 10 
C
      write(10,*) ' '
      write(10,*) 'delEes delEvdw plztot sum E kcal/mol'
      write(10,'(a23,4f14.5)')
     &  ' dEcompo es vdw plz tot',dees*au,devdw*au,plztot*au,
     &  (dees+devdw+plztot)*au
      write(10,*) 
     & '------------------------------------------------------------'
      write(10,*) 
     & '  Dipole moments of Mass Center (debye) x y z |mu|'
      write(10,*) 
     & '------------------------------------------------------------'
      write(10,*) '       Gaussian/ monomer charges'
      write(10,'(a10,10x,4f10.5)') 'GMC MC   ',
     &  tcdipx*debye,tcdipy*debye,tcdipz*debye,tcdip*debye
      write(10,*) '       Gaussian/ cluster charges'
      write(10,'(a10,10x,4f10.5)') ' GCC MC   ', 
     &  gdipx*debye,gdipy*debye,gdipz*debye,gdip*debye
      write(10,*) '       Estimated total'  
      write(10,'(a10,10x,4f10.5)') ' ET  MC   ',alldipx*debye,
     &   alldipy*debye,alldipz*debye,alldip*debye
      write(10,*) 
     & '------------------------------------------------------------'
      write(10,*) 
     & '  Coordinate Dependent Dipole moments (debye) x y z |mu|'
      write(10,*)
     &  '-----------------------------------------------------------'
      write(10,*) '       Gaussian/ monomer charges'
      write(10,'(a10,10x,4f10.5)') ' GMC CD   ', 
     &  tmdipx*debye,tmdipy*debye,tmdipz*debye,tmdip*debye
      write(10,*) '       Gaussian/ cluster charges'
      write(10,'(a10,10x,4f10.5)') ' GCC CD   ',
     &  gdicx*debye,gdicy*debye,gdicz*debye, gdic*debye
      write(10,*) '       Estimated total'   
      write(10,'(a10,10x,4f10.5)') ' ET  CD   ',alldimx*debye,
     &   alldimy*debye,alldimz*debye,alldim*debye
      write(10,*)
     &  '-----------------------------------------------------------'
      write(10,*) '       Estimated induced dipole moment (debye)'
      write(10,*) 
     &  '-----------------------------------------------------------'
      write(10,'(a10,10x,4f10.5)') ' EID      ',totdpx*debye,
     &   totdpy*debye,totdpz*debye,totdp*debye
      write(10,*) 
     &  '-----------------------------------------------------------'
*     write(10,*) 'end of subr.plzcyc'
      return
      end 
************************************************************************
      subroutine indatm(n,dafldx,dafldy,dafldz,
     &                xatdip,yatdip,zatdip,vpol)
      parameter(ks30=30,k200=300)
      parameter (k50=250,k100=30,k500=500,k4=4,i50=500)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
C
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/fldcharg/
     &    cafldx(k100,k50),cafldy(k100,k50),
     &    cafldz(k100,k50)
      dimension dafldx(k50),dafldy(k50),dafldz(k50),
     &          xatdip(k50),yatdip(k50),zatdip(k50) 
      data one/1.0d0/,two/2.0d0/
      data zero/0.0d0/,bohr/0.52917706D0/,aukcal/627.509541d0/
C
      druk=1000.0d0/(bohr*bohr*aukcal)
      na=nwaplz(n)
      nb=nwbplz(n)
      mb=mwbplz(n)
      vpol=zero
      do i=1,nwaplz(n)
        fldc=dsqrt(cafldx(n,i)**2+cafldy(n,i)**2+cafldz(n,i)**2)
        fld1=dsqrt(dafldx(i)**2+dafldy(i)**2+dafldz(i)**2)
        fld2=fld1*fld1
        fld3=fld2*fld1
        alp=wplz(n,i,1)+wplz(n,i,2)*fld1+wplz(n,i,3)*fld2+
     &      wplz(n,i,4)*fld3
        xatdip(i)=alp*dafldx(i)
        yatdip(i)=alp*dafldy(i)
        zatdip(i)=alp*dafldz(i)
        if(iwrest(n).ne.7)then
          vpol=vpol-(wplz(n,i,1)/2.0d0 
     &              +wplz(n,i,2)*fld1/3.0d0
     &              +wplz(n,i,3)*fld2/4.0d0
     &              +wplz(n,i,4)*fld3/5.0d0)
     &    *(dafldx(i)*cafldx(n,i)+dafldy(i)*cafldy(n,i)
     &     +dafldz(i)*cafldz(n,i))
        else
          q=dsqrt(alp*druk)
          disx=xatdip(i)/q
          disy=yatdip(i)/q
          disz=zatdip(i)/q
          vpol=vpol-druk*(disx**2+disy**2+disz**2)/two
        endif
      enddo
      do i=1,nwbplz(n)
        fldc=dsqrt(cafldx(n,i+na)**2+cafldy(n,i+na)**2+
     &      cafldz(n,i+na)**2)
        cafx=-cafldx(n,i+na)/fldc
        cafy=-cafldy(n,i+na)/fldc
        cafz=-cafldz(n,i+na)/fldc
        canifac=(bvecx(n,i)*cafx+bvecy(n,i)*cafy+bvecz(n,i)*cafz)
        cafx=cafldx(n,i+na)*dabs(canifac)
        cafy=cafldy(n,i+na)*dabs(canifac)
        cafz=cafldz(n,i+na)*dabs(canifac)
C 
        fld0=dsqrt(dafldx(i+na)**2+dafldy(i+na)**2+dafldz(i+na)**2)
        dafx=-dafldx(i+na)/fld0
        dafy=-dafldy(i+na)/fld0
        dafz=-dafldz(i+na)/fld0
        danifac=(bvecx(n,i)*dafx+bvecy(n,i)*dafy+bvecz(n,i)*dafz)
        if(danifac.ge.zero)then
          sigd=1.0d0
        else
          sigd=-1.0d0
        endif
        bx=bvecx(n,i)*sigd
        by=bvecy(n,i)*sigd
        bz=bvecz(n,i)*sigd
        dafx=dafldx(i+na)*dabs(danifac)
        dafy=dafldy(i+na)*dabs(danifac)
        dafz=dafldz(i+na)*dabs(danifac)
        fld1=dsqrt(dafx**2+dafy**2+dafz**2)
        fld2=fld1*fld1
        fld3=fld2*fld1
        alp=wplz(n,i+na,1)+wplz(n,i+na,2)*fld1+wplz(n,i+na,3)*fld2+
     &      wplz(n,i+na,4)*fld3
        xatdip(i+na)=-fld1*alp*bx
        yatdip(i+na)=-fld1*alp*by
        zatdip(i+na)=-fld1*alp*bz
        vpol=vpol-(wplz(n,i+na,1)/2.0d0 
     &            +wplz(n,i+na,2)*fld1/3.0d0
     &            +wplz(n,i+na,3)*fld2/4.0d0
     &            +wplz(n,i+na,4)*fld3/5.0d0)
     &  *(dafx*cafx+dafy*cafy+dafz*cafz)
      enddo
C
      do i=1,mb
        j=i+na+nb
        fldc=dsqrt(cafldx(n,j)**2+cafldy(n,j)**2+cafldz(n,j)**2)
        fld1=dsqrt(dafldx(j)**2+dafldy(j)**2+dafldz(j)**2)
        fld2=fld1*fld1
        fld3=fld2*fld1
        alp=wplz(n,j,1)+wplz(n,j,2)*fld1+wplz(n,j,3)*fld2+
     &      wplz(n,j,4)*fld3
        xatdip(j)=alp*dafldx(j)
        yatdip(j)=alp*dafldy(j)
        zatdip(j)=alp*dafldz(j)
        vpol=vpol-(wplz(n,j,1)/2.0d0 
     &            +wplz(n,j,2)*fld1/3.0d0
     &            +wplz(n,j,3)*fld2/4.0d0
     &            +wplz(n,j,4)*fld3/5.0d0)
     &  *(dafldx(j)*cafldx(n,j)+dafldy(j)*cafldy(n,j)
     &      +dafldz(j)*cafldz(n,j))
      enddo
      return
      end
************************************************************************
      subroutine popsrf
      parameter(ks30=30,k200=300)
      parameter (k50=250,k300=300,k3000=5000,k4=4,i50=500,k100=30)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/espot/ cx(k300),cy(k300),cz(k300),sx(k3000),sy(k3000),
     &  sz(k3000),spot(k3000),gcrg(k300),icoor,isurf,
     &  pop(k3000),gpop(k3000)
      common/inddip/atdipx(k100,k50),atdipy(k100,k50),atdipz(k100,k50),
     &  nnot
      dimension cpop(k3000)
      data zero/0.0d0/,au/627.509541D0/,oneh/100.0d0/,
     &  bohr/0.52917706D0/
C
      rmsp=zero
      rmsg=zero
      rmsc=zero
      srms=zero
      drmsp=zero
      drmsg=zero
      drmsc=zero
      do i=1,isurf
        pop(i)=zero
        gpop(i)=zero
        cpop(i)=zero
        kk=0
        do j=1,nseq
          na=nwaplz(j)
          nb=nwbplz(j)
          mb=mwbplz(j)
          do k=1,na+iwbonc(j)
            r=dsqrt((wx(j,k)-sx(i))**2+(wy(j,k)-sy(i))**2
     &             +(wz(j,k)-sz(i))**2)
            e=wcrg(j,k)/r
            pop(i)=pop(i)+e
            if(k.le.na)then
              kk=kk+1
              g=gcrg(kk)/r
              gpop(i)=gpop(i)+g
            endif
            cpop(i)=cpop(i)+e
          enddo
C
C         intermolecular delta charge   model a
C
          if(iwrest(1).lt.2)then
            do k=1,na
              q=dsqrt(atdipx(j,k)**2+atdipy(j,k)**2+atdipz(j,k)**2)
              q=q/1.0d0
              if(q.ne.zero)then
                x=0.5d0*atdipx(j,k)/q
                y=0.5d0*atdipy(j,k)/q
                z=0.5d0*atdipz(j,k)/q
              else
                x=zero
                y=zero
                z=zero
              endif
              xp=wx(j,k)+x
              yp=wy(j,k)+y
              zp=wz(j,k)+z
              xm=wx(j,k)-x
              ym=wy(j,k)-y
              zm=wz(j,k)-z
              r1=dsqrt((xp-sx(i))**2+(yp-sy(i))**2+(zp-sz(i))**2)
              r2=dsqrt((xm-sx(i))**2+(ym-sy(i))**2+(zm-sz(i))**2)
              p=q/r1-q/r2
              pop(i)=pop(i)+p
            enddo
C
C         intra- and intermolecular polz model
C
          else if((iwrest(1).ge.3.and.iwrest(1).lt.7).or.
     &            (iwrest(1).eq.13).or.
     &            (iwrest(1).eq.14).or.(iwrest(1).eq.15).or.
     &            (iwrest(1).eq.143).or.(iwrest(1).eq.153).or.
     &            (iwrest(1).eq.163).or.(iwrest(1).eq.173).or.
     &            (iwrest(1).eq.24))then
            do k=1,na
              dx=atdipx(j,k)
              dy=atdipy(j,k)
              dz=atdipz(j,k)
              ax=sx(i)-wx(j,k)
              ay=sy(i)-wy(j,k)
              az=sz(i)-wz(j,k)
              a=dsqrt(ax**2+ay**2+az**2)
              d=dsqrt(dx**2+dy**2+dz**2)
C
              p=(ax*dx+ay*dy+az*dz)/a**3
              pop(i)=pop(i)+p
            enddo
          else if(iwrest(1).eq.7)then
            druk=1000.0d0/(bohr*bohr*au)
            do k=1,na
              q=dsqrt(wplz(j,k,1)*druk)
              if(q.ne.zero)then
                x=atdipx(j,k)/q
                y=atdipy(j,k)/q
                z=atdipz(j,k)/q
              else
                x=zero
                y=zero
                z=zero
              endif
              xp=wx(j,k)
              yp=wy(j,k)
              zp=wz(j,k)
              xm=wx(j,k)-x
              ym=wy(j,k)-y
              zm=wz(j,k)-z
              r1=dsqrt((xp-sx(i))**2+(yp-sy(i))**2+(zp-sz(i))**2)
              r2=dsqrt((xm-sx(i))**2+(ym-sy(i))**2+(zm-sz(i))**2)
              p=q/r1-q/r2
              pop(i)=pop(i)+p
            enddo
          endif
C
C         model b
C
          do k=1,nb
            q=dsqrt(atdipx(j,k+na)**2+atdipy(j,k+na)**2+
     &        atdipz(j,k+na)**2)/bonlen(j,k)
            x=0.5d0*atdipx(j,k+na)/q
            y=0.5d0*atdipy(j,k+na)/q
            z=0.5d0*atdipz(j,k+na)/q
            xp=wx(j,k+na)+x
            yp=wy(j,k+na)+y
            zp=wz(j,k+na)+z
            xm=wx(j,k+na)-x
            ym=wy(j,k+na)-y
            zm=wz(j,k+na)-z
            r1=dsqrt((xp-sx(i))**2+(yp-sy(i))**2+(zp-sz(i))**2)
            r2=dsqrt((xm-sx(i))**2+(ym-sy(i))**2+(zm-sz(i))**2)
            p=q/r1-q/r2
            pop(i)=pop(i)+p
          enddo
C
C         model c
C
          do k=1,mb
            l=k+na+nb
            q=dsqrt(atdipx(j,l)**2+atdipy(j,l)**2+
     &        atdipz(j,l)**2)
            x=0.5d0*atdipx(j,l)/q
            y=0.5d0*atdipy(j,l)/q
            z=0.5d0*atdipz(j,l)/q
            xp=wx(j,l)+x
            yp=wy(j,l)+y
            zp=wz(j,l)+z
            xm=wx(j,l)-x
            ym=wy(j,l)-y
            zm=wz(j,l)-z
            r1=dsqrt((xp-sx(i))**2+(yp-sy(i))**2+(zp-sz(i))**2)
            r2=dsqrt((xm-sx(i))**2+(ym-sy(i))**2+(zm-sz(i))**2)
            p=q/r1-q/r2
            pop(i)=pop(i)+p
          enddo
        enddo
        difpot=pop(i)+spot(i)
        gdifpot=gpop(i)+spot(i)
        cdifpot=cpop(i)+spot(i)
C       write(6,'(i3,7f8.3)') 
C    & i,spot(i),pop(i),gpop(i),cpop(i),difpot,gdifpot,cdifpot
        drmsp=drmsp+difpot**2
        drmsg=drmsg+gdifpot**2
        drmsc=drmsc+cdifpot**2
        srms=srms+spot(i)**2
      enddo
      drmsp=dsqrt(drmsp/isurf)
      drmsg=dsqrt(drmsg/isurf)
      drmsc=dsqrt(drmsc/isurf)
      srms=dsqrt(srms/isurf)
      srmsp=oneh*drmsp/srms
      srmsg=oneh*drmsg/srms
      srmsc=oneh*drmsc/srms
      write(10,*) '                POP       cluster   monomer' 
      write(10,'(a15,3f10.5)') 'rms dev au     ',
     &     drmsp,drmsg,drmsc
      drmsp=drmsp*au
      drmsg=drmsg*au
      drmsc=drmsc*au
      write(10,'(a15,3f10.5)') 'rms dev kca/mol',
     &     drmsp,drmsg,drmsc
      if(nnot.eq.500)then
        write(10,'(a15,3f10.5,a19)') 'rms dev % spot ',
     &  srmsp,srmsg,srmsc,'    NOT converged! '
      else if(nnot.eq.400)then
        write(10,'(a15,3f10.5,a19)') 'rms dev % spot ',
     &  srmsp,srmsg,srmsc,'    STOP diverged! '
      else
        write(10,'(a15,3f10.5)') 'rms dev % spot ',
     &     srmsp,srmsg,srmsc
      endif
      return
      end
************************************************************************
      subroutine spoppd(i,j,k,dx,dy,dz,p)
C
C     point dipole model 
C
      parameter(ks30=30,k200=300)
      parameter (k50=250,k300=300,k3000=5000,k4=4,i50=500,k100=30)
      implicit real*8 (a-h,o-z)
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
      common/espot/ cx(k300),cy(k300),cz(k300),sx(k3000),sy(k3000),
     &  sz(k3000),spot(k3000),gcrg(k300),icoor,isurf,
     &  pop(k3000),gpop(k3000)
      data one/1.0d0/,three/3.0d0/
      ax=sx(i)-wx(j,l)
      ay=sy(i)-wy(j,l)
      az=sz(i)-wz(j,l)
      a=dsqrt(ax**2+ay**2+az**2)
      d=dsqrt(dx**2+dy**2+dz**2)
C
      p=(ax*dx+ay*dy+az*dz)/a**3
      return
      end
***********************************************************************
      subroutine fscale(i,j,k,f12,f13,f14,f15,flnk)
C
C     Preparation for 1-2,1-3,1-4,1-5... link data 
C         factor of 1-2 is f12
C         factor of 1-3 is f13
C         factor of 1-4 is f14 
C         factor of 1-5... is f15 
C
      IMPLICIT REAL*8(A-H,O-Z)
      parameter(k100=30,k50=250)
      parameter(kr50=50,ka200=300,k4=4,k21=21)
      common/topol1/ tcrg(kr50),crg(kr50,ka200),
     &  atmplz(kr50,ka200,k4),bonplz(kr50,ka200,k4),natmplz(kr50),
     &  nbonplz(kr50),nbon1(kr50,ka200),nbon2(kr50,ka200),
     &  atmass(kr50,ka200),
     &  xrmin(kr50,ka200),xemin(kr50,ka200),mbon1(kr50,ka200),
     &  mbon2(kr50,ka200),boaplz(kr50,ka200,4),mbonplz(kr50),
     &  mseg(kr50,ka200)
      common/topol2/
     &  nlnk(kr50),nl12(kr50,ka200,ka200),nl13(kr50,ka200,ka200),
     &  nl14(kr50,ka200,ka200),scac14(kr50),scav14(kr50),
     &  facp12(kr50),facp13(kr50),facp14(kr50),facp15(kr50),
     &  facc12(kr50),facc13(kr50),facc14(kr50),facc15(kr50),ibonc(kr50),
     &  bchg(kr50,ka200),irest(kr50),psfac(kr50),factt(kr50)
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
      dimension flnk(k50,k50)
      data zero/0.0d0/,one/1.0d0/
C
C     clear matrix 0.0 or f15
C
C     write(6,*)' % subr fscale '
      do i1=1,k
        do j1=i1,k
          flnk(i1,j1)=f15
          if(i1.eq.j1) flnk(i1,j1)=zero
        enddo
      enddo
C
C     scale f12 for 1-2 intramolecular interaction
C
      do i1=1,k
        do j1=2,k
          if(nl12(j,i1,j1).ne.0)then
            flnk(nl12(j,i1,1),nl12(j,i1,j1))=f12
          endif
        enddo
      enddo
C
C     scale f13 for 1-3 intramolecular interaction
C
      do i1=1,k
        do j1=2,k
          if(nl13(j,i1,j1).ne.0)then
            flnk(nl13(j,i1,1),nl13(j,i1,j1))=f13
          endif
        enddo
      enddo
C
C     scale f14 for 1-4 intramolecular interaction
C
      if(f14.gt.1.0d0.or.f14.lt.0.0d0) f14=1.0d0/1.2d0
      do i1=1,k
        do j1=2,k
          if(nl14(j,i1,j1).ne.0)then
            flnk(nl14(j,i1,1),nl14(j,i1,j1))=f14
          endif
        enddo
      enddo
C
      do i1=1,k
        do j1=i1,k
          flnk(j1,i1)=flnk(i1,j1)
        enddo
      enddo
C
C     do i1=1,k
C       write(6,'(50f5.2)') (flnk(i1,j1),j1=1,k)
C     enddo
      return
      end
***********************************************************************
      subroutine lnkplz(ii)
c
c     preparation for connectivity and esp optimized parmanent charges
c
      implicit real*8(a-h,o-z)
      parameter(k100=30,k50=250)
      parameter(kr50=50,ka200=300,k4=4,k21=21)
      common/topol1/ tcrg(kr50),crg(kr50,ka200),
     &  atmplz(kr50,ka200,k4),bonplz(kr50,ka200,k4),natmplz(kr50),
     &  nbonplz(kr50),nbon1(kr50,ka200),nbon2(kr50,ka200),
     &  atmass(kr50,ka200),
     &  xrmin(kr50,ka200),xemin(kr50,ka200),mbon1(kr50,ka200),
     &  mbon2(kr50,ka200),boaplz(kr50,ka200,4),mbonplz(kr50),
     &  mseg(kr50,ka200)
      common/topol2/
     &  nlnk(kr50),nl12(kr50,ka200,ka200),nl13(kr50,ka200,ka200),
     &  nl14(kr50,ka200,ka200),scac14(kr50),scav14(kr50),
     &  facp12(kr50),facp13(kr50),facp14(kr50),facp15(kr50),
     &  facc12(kr50),facc13(kr50),facc14(kr50),facc15(kr50),ibonc(kr50),
     &  bchg(kr50,ka200),irest(kr50),psfac(kr50),factt(kr50)
      dimension lnk1234(k50,k50,3),ln(k50),kk(k50),k3(k50),kk2(k50)
      dimension lnk12bf(k50,k50,3)
      data zero/0.0d0/,one/1.0d0/
C  
C     initial clear or set & distance check
C
C
C     read 1-2 bond data
C
      iuniq=nlnk(ii)
      do i=1,iuniq
        do j=1,iuniq
          lnk1234(i,j,1)=nl12(ii,i,j)
        enddo
      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
C       write(6,'(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
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
      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
C!                if(lnk12bf(k,l,1).ne.lnk12bf(i,1,1))then
                  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
          if(lnk1234(i,j,3).ne.0)then
            kka=kka+1
            ln(kka)=lnk1234(i,j,3)
          endif
        enddo
        do j=1,iuniq
          lnk1234(i,j,3)=0
        enddo
        do j=1,kka
          lnk1234(i,j,3)=ln(j)
        enddo 
C       write(6,'(i4,a5,90i4)') kka,'%1-4 ',(lnk1234(i,j,3),j=1,kka)
      enddo
C
      do i=1,iuniq
        do j=1,iuniq
          nl13(ii,i,j)=lnk1234(i,j,2)
          nl14(ii,i,j)=lnk1234(i,j,3)
        enddo
      enddo
C
      return
      end
***********************************************************************
      subroutine matinv(dafldx,dafldy,dafldz,xdip,ydip,zdip,xe,ye,ze)
      parameter(ks30=30,k200=300)
      parameter(k50=250,k100=30,k500=500,k4=4,i50=500)
      parameter(i200=250)
      parameter(k300=300,k3000=5000)
      implicit real*8 (a-h,o-z)
      character*4 seq(ks30)
      common/sequ/ nseq,seq
      common/wcoor/ wx(ks30,k200),wy(ks30,k200),wz(ks30,k200),
     &  wcrg(ks30,k200),wplz(ks30,k200,k4),nwaplz(ks30),nwbplz(ks30),
     &  wxmass(ks30),wymass(ks30),wzmass(ks30),xmass,ymass,zmass,
     &  wamass(ks30,k200),bvecx(ks30,k200),bvecy(ks30,k200),
     &  bvecz(ks30,k200)
     & ,bonlen(ks30,k200),nwbon1(ks30,k200),nwbon2(ks30,k200),
     &  wrmin(ks30,k200),wemin(ks30,k200),mwbplz(ks30),
     &  mwbon1(ks30,k200),
     &  mwbon2(ks30,k200),nwlnk(ks30),iwbonc(ks30),iwrest(ks30),
     &  wfacp(ks30),nseg(ks30,k200)
C
      common/dist/ rx(k100,k100,k50,k50),ry(k100,k100,k50,k50),
     &             rz(k100,k100,k50,k50),rdis(k100,k100,k50,k50)
      common/fldcharg/
     &    cafldx(k100,k50),cafldy(k100,k50),
     &    cafldz(k100,k50)
      common/dipole/tcdipx,tcdipy,tcdipz,tmdipx,tmdipy,tmdipz,
     &  gdipx,gdipy,gdipz,gdicx,gdicy,gdicz
      common/energy/ dees,es(k100,k100),devdw,vdw(k100,k100),
     &  plztot,plz(k100,k100),totene,tot(k100,k100) 
      common/inddip/atdipx(k100,k50),atdipy(k100,k50),atdipz(k100,k50),
     &  nnot
      common/wintra/scac(k100,k50,k50),scav(k100,k50,k50),
     &  facp(k100,k50,k50),facc(k100,k50,k50),
     &  se14(k100,k50,k50),sv14(k100,k50,k50)
cc
      common/espot/ cx(k300),cy(k300),cz(k300),sx(k3000),sy(k3000),
     &  sz(k3000),spot(k3000),gcrg(k300),icoor,isurf,
     &  pop(k3000),gpop(k3000)
      common/tens/
     &  txx(k50,k50),txy(k50,k50),
     &  txz(k50,k50),
     &  tyx(k50,k50),tyy(k50,k50),
     &  tyz(k50,k50),
     &  tzx(k50,k50),tzy(k50,k50),
     &  tzz(k50,k50)
!*--------------------------------------------------------------------------
      dimension A(i200*3,i200*2*3),DA(i200*3),E(i200*3)
!*--------------------------------------------------------------------------
cc
      dimension dafldx(k50),dafldy(k50),dafldz(k50),
     &          atdipox(k100,k50),atdipoy(k100,k50),atdipoz(k100,k50),
     &          xatdip(k50),yatdip(k50),zatdip(k50),
     &          vpolatm(k100)
      dimension xdip(i200),ydip(i200),zdip(i200),
     &          xe(i200),ye(i200),ze(i200)
      data zero/0.0d0/,two/2.0d0/
      data au/627.509541D0/
      data debye/2.54158059d0/
C

!*----------------------------------------------------------
!*    Matrix A
!*----------------------------------------------------------
      k=1
      do i=1,icoor
        l=1
        do j=1,icoor
          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,icoor
          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,icoor
          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(dafldx,dafldy,dafldz,icoor,A,DA,E)
!*   
!*----------------------------------------------------------
c
c     calculation for intramolecular induced dipoles
c
c
      j=0
      do i=1,icoor*3,3
       j=j+1
       xdip(j)=DA(i)
       ydip(j)=DA(i+1)
       zdip(j)=DA(i+2)
       xe(j)=E(i)
       ye(j)=E(i+1)
       ze(j)=E(i+2)
      enddo
C
      return
      end
***********************************************************************
      Subroutine houhol0(dafldx,dafldy,dafldz,nn,C,DI,E)
      implicit real*8(a-h,o-z)
      parameter(ks30=30,k200=300)
      parameter(k50=250,k100=30,k500=500,k4=4,i50=500)
      parameter(i200=250)
C     integer n   !size of matrix A
cc
      common/fldcharg/
     &    cafldx(k100,k50),cafldy(k100,k50),
     &    cafldz(k100,k50)
cc
      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),E(i200*3)
      dimension dafldx(k50),dafldy(k50),dafldz(k50) 
*     data zero/0.0d0/
      data au/627.509541D0/,zero/0.0d0/,two/2.0d0/
!*
!*    A Matrix to be inverted
!*
      n=nn*3
      do i=1,n
        do j=1,n
          A(i,j)=C(i,j)
        enddo
      enddo
c
!*     
!*    dafldx y z  F Field from a fld charge
!* 
      k=1
      do i=1,n,3
        F(i)=dafldx(k)
        F(i+1)=dafldy(k)
        F(i+2)=dafldz(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     ' Inverted matrix:'
cc
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)=zero
      enddo
      do i=1,n
        do j=1,n
          DI(i)=DI(i)+AI(i,j)*F(j)        
        enddo
      enddo
cc
      call molpol(n,AI)
      call dipf(n,DI)
      call upol(n,DI,F,E)
C
      return
      End 
***********************************************************************
      Subroutine S1000(n,A,V)
      implicit real*8(a-h,o-z)
      parameter(i200=250)

C integer n
      dimension A(i200*3,i200*2*3),V(i200*3)
      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
c
      end do !k loop
      return
      End
***********************************************************************
      Subroutine S2000(n,A,B,X)
      implicit real*8(a-h,o-z)
      parameter(i200=250)
C  integer n
      dimension A(i200*3,i200*2*3),B(i200*3),X(i200*3)
      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=250)
C integer n
      dimension A0(i200*3,i200*3),AI(i200*3,i200*3),AP(i200*3,i200*3)
      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=250)

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
!end of file householder.f90
***********************************************************************
      Subroutine molpol(n,AI)
C
C     Molecular Polarizability
C
      implicit real*8(a-h,o-z)
      parameter(i200=250)

C integer n
      dimension AI(i200*3,i200*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)+AI(i,j) 
          A1(i,2)=A1(i,2)+AI(i,j+1) 
          A1(i,3)=A1(i,3)+AI(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
cc
      write(10,*) '  '
      write(10,*) 'Molecular Polarizability  lower triangle'
      write(10,'(a10,6f12.3)')' Mol Pol  ',
     &   AM(1,1),AM(2,1),AM(2,2),AM(3,1),AM(3,2),AM(3,3)
cc
      return
      End
***********************************************************************
      Subroutine dipf(n,DI)
C     
C     Induced Dipole Moment
C
      implicit real*8(a-h,o-z)
      parameter(i200=250)
      dimension DI(i200*3),DP(3)
      data debye/2.54158059d0/
cc
      do i=1,3
        DP(i)=zero
      enddo
      do i=1,n,3
        DP(1)=DP(1)+DI(i)
        DP(2)=DP(2)+DI(i+1)
        DP(3)=DP(3)+DI(i+2)
      end do
      DT=dsqrt(DP(1)**2+DP(2)**2+DP(3)**2)
cc
      write(10,*) '  '
      write(10,*) 'Induced Dipole Moments from Matrix Inversion (debye)'
      write(10,'(a10,4f10.5)')' x y z tot', 
     &  DP(1)*debye,DP(2)*debye,DP(3)*debye,DT*debye
cc
      return
      End
***********************************************************************
      Subroutine upol(n,DI,F,E)
C
C     Polarization Energy 
C
      implicit real*8(a-h,o-z)
      parameter(i200=250)
      dimension F(i200*3),DI(i200*3),E(i200*3)
      data au/627.509541D0/,zero/0.0d0/,two/2.0d0/
cc
      vpol=zero
      do i=1,n
        E(i)=DI(i)*F(i)
        vpol=vpol-E(i)
      end do
      vpol=vpol*au/two
cc
      write(10,*) '  '
      write(10,*) 'Polarization Energy (kcal/mol)'
      write(10,'(a10,f10.5)')' Pol Ene  ',vpol
      write(10,*) '  '
cc
      return
      End
***********************************************************************