PROGRAM DODCG2
      IMPLICIT none
      integer*4 iper,nt0,it,i,N,J,nc,ns,jx,MXB,n6,npar0,
     1ig,is,ix,NG,nat,ncs,ierr,ITER,NATB,nwl,nwl0,n4,n40,
     2ivel,imag,ntr,nd,nt
      PARAMETER (nt0=550,MXB=50,npar0=14,n40=6,nwl0=nt0)
c     Limits:
c     nt0 -# of transitions on one chromophore
c     MXB -# of overlaped atoms
c     npar0 -# of parameters to follow
c     n40 -# of potential sites
c     nwl0 -# of fitted transitions
      integer*4 IP(npar0),im,imax,imin,nm,nmol
      REAL*8 rc(3),mag(nt0,3),elel(nt0,3),
     1elev(nt0,3),en(nt0),en1(nt0,3),w,A,XS,XB,RTOL,um(3),ue(3),dco,
     1PG(nt0,npar0),zim(n40),dco0,eau,ecm(nt0),eev(nt0),enm(nt0),
     1buvcd(nwl0,n40),fi(n40),F(npar0),enc(nt0),
     1dcn,dcn0,dhn,dhn0,DPCX,DPCY,DPCZ,elelc(nt0,3),rcs(3)
      REAL*8,allocatable:: R(:),rn(:),wc(:)
      integer*4,allocatable:: ic(:),ics(:),ttype(:),id(:),tl(:)
      character*80 s80,fn
      character*1 prog
      logical lper,lfdp(4),lvel,llen,lene,lmag,lgeo,take
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NATB,ITER
c
c     takes electronic dipoles from G.OUT, spreads them on large peptide
c
c     Input: G.OUT       --   TDDFT G Output
c            or moments  --   TDDFT Turbomole output after transformation byc            program turbomoments
c            or F.INP files for vibrations and more chromophores
c            polymer.x  --- the whole system
c
c     Read in transition energies and dipoles from Gaussian or Turbomole:
      lvel=.false.
      llen=.false.
      lene=.false.
      lmag=.false.
      lgeo=.false.

c     read-in geometry
      OPEN(2,FILE='BIG.X',STATUS='OLD')
      read(2,*)
      read(2,*)N
      allocate (R(3*N),ttype(N),rn(3*N))
      do 1 i=1,N
1     READ(2,*)ttype(i),(R(J+3*(i-1)),J=1,3)
      close(2)
      write(6,*)N, ' atoms in BIG.X'
c
      open(4,file='DODC.CON')
      read(4,401)prog
401   format(a1)
      read(4,*)ivel,imag
      if(imag.eq.1)then
       write(6,*)' magnetic moments added'
      else
       write(6,*)' magnetic moments skipped'
      endif

      if(prog.eq.'g'.or.prog.eq.'G')then
c      gaussian electric:
       open(30,file='G.OUT',status='old')
10     read(30,300,end=949,err=949)s80
300    format(a80)
       if(s80(12:44).eq.'excited state Transition electric')then
        call rd(nt0,it,elel,llen)
        write(6,*)it,' electric dipoles'
        goto 10
       endif
       if(s80(12:44).eq.'excited state transition electric')then
        call rd(nt0,it,elel,llen)
        write(6,*)it,' electric dipoles'
        goto 10
       endif
       if(s80(12:46).eq.'excited state transition velocity d')then
        call rd(nt0,it,elev,lvel)
        write(6,*)it,' electric velocity dipoles'
        goto 10
       endif
       if(s80(12:44).eq.'excited state transition magnetic')then
        call rd(nt0,it,mag,lmag)
        write(6,*)it,' magnetic dipoles'
        goto 10
       endif
       if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s')then
        lene=.true.
        read(s80(15:18),*)nd
        do 203 i=1,79
203     if(s80(i:i+1).eq.'nm')read(s80(i- 9:i-2),*)en(nd)
       endif
       if(s80(26:46).eq.'Standard orientation:')call ro(rn,nat,lgeo)
       goto 10
949    close(30)
       write(6,*)nd,' energy transitions in G.OUT'
       do 481 i=1,it
c      nm to hartrees:
       eau=1.0d7/(8065.54476345045d0*27.211384205943d0)/en(i)
       do 481 ix=1,3
c      correct also gaussian sign:
481    elev(i,ix)=-elev(i,ix)/eau
      else
       if(prog.eq.'t'.or.prog.eq.'T')then
c       turbomole electric:
        open(30,file='moments',status='old')
100     read(30,300,end=9490,err=9490)s80
        if(s80(2:46).eq.'Transition dipole moments in length formalism')
     1  call rdt(nt0,it,elel,llen)
        if(s80(2:46).eq.'Transition dipole moments in velocity formali')
     1  call rdt(nt0,it,elev,lvel)
        if(s80(2:28).eq.'Transition magnetic moments')
     1  call rdt(nt0,it,mag,lmag)
        if(s80(2:24).eq.'Energies of Transitions')
     1  call rdt(nt0,it,en1,lene)
        if(s80(2:19).eq.'Atomic coordinates')call rot(rn,nat,lgeo)
        goto 100
9490    close(30)
        do 9432 i=1,it
9432    en(i)=en1(i,2)
        write(6,*)it,' transitions in moments'
       else
        if(prog.eq.'f'.or.prog.eq.'F')then
c        vibrational from F.INP:
c
c        find total number of molecules and groups:
         NG=0
         rewind 4
         read(4,*)
         read(4,*)
         read(4,*)nmol
         write(6,*)nmol,' source molecules'

         do 41 im=1,nmol
c        read corresponding list of atoms in BIG.X, to find out
c        number of groups:
         do 411 it=1,4
411      read(4,*)
         read(4,400)fn
         open(8,file=fn)
         it=0
42       read(8,*,end=88,err=88)ix
         it=it+1
         goto 42
88       close(8)
         write(6,*)fn
         NG=NG+it
         write(6,*)it,' groups'
41       read(4,*)

c        write DC.SC, reread DODC.CON:
         open(3,file='DC.SC')
         write(3,*)NG,' total groups'
         write(6,*)NG,' total groups'

         IG=0
         rewind 4
         do 412 it=1,3
412      read(4,*)
         do 40 im=1,nmol
         read(4,*)imin,imax
         if(imin.lt.0)then
c         list of transitions instead of interval:
          nt=-imin
          allocate(tl(nt))
          backspace 4
          read(4,*)tl(1),(tl(i),i=1,nt)
         endif
         read(4,400)fn
400      format(a80)
         write(6,*)fn
c        read dipoles and energies from F.INP:
         open(30,file=fn,status='old')
         read(30,*)nm,nat,nat
         allocate(id(nat))
         close(30)
         call rvd(fn,nt0,nat,nm,rn,id,mag,elel,en)
c        chromophore atoms, e.g.:   H N C O
         read(4,*)nc
         allocate (ics(nc),wc(nc),ic(nc))
         backspace 4
         read(4,*)nc,(ics(i),i=1,nc)
c        atomic weights for placing of the dipole:
         read(4,*)(wc(i),i=1,nc)
         if(nc.gt.MXB)call report('nc>MXB')
c        make center of the "small" part:
         call makecenter(rcs,rn,wc,nc,ics)
c        list of atoms in BIG.X:        
         read(4,400)fn
         open(8,file=fn)
43       read(8,*,end=98,err=98)(ic(i),i=1,nc)
         IG=IG+1
         write(3,*)IG,' . group:'
c        make center of the "big" part:
         call makecenter(rc,R,wc,nc,ic)
         write(3,1000)(rc(i)/10.0d0,i=1,3)
c        write group reference atoms into XS and corresponding ones in XB:
         do 49 i=1,nc
         do 49 ix=1,3
         XB(i,ix)= R(ix+3*(ic( i)-1))-rc( ix)
49       XS(i,ix)=rn(ix+3*(ics(i)-1))-rcs(ix)
         NATB=nc
c        overlap XS,XB, make the rotation matrix, store in A:
         call domatrix(ierr)
         if(ierr.ne.0)call report('rotation not successful')
         write(6,*)' Transformation matrix'
         do 67671 i=1,3
67671    write(6,6767)(A(i,j),j=1,3)

         ntr=0
         do 472 i=1,nm
         take=.false.
         if(imin.lt.0)then
c         list of transitions:
          do 62 ix=1,nt
62        if(tl(ix).eq.i)take=.true.
         else
c         interval of transitions:
          take=i.ge.imin.and.i.le.imax
         endif
         if(take)then
          ntr=ntr+1
          ecm(ntr)=en(i)
          enm(ntr)=1.0d7/ecm(ntr)
          eev(ntr)=ecm(ntr)/8065.54093d0
          do 473 ix=1,3
c         borrow elev variable for magnetic:
          if(imag.eq.0)mag(i,ix)=0.0d0
          elev(ntr,ix)=mag(i,ix)
473       elelc(ntr,ix)=elel(i,ix)
         endif
472      continue

         ns=ntr+1
         write(3,*)ns,' states'
         write(3,*)
c        ground state:
         write(3,*)'0 0.0'
         do 48 is=2,ns
48       write(3,1001)is,eev(is-1),ecm(is-1),enm(is-1)
         write(3,*)
c        write rotated dipole moments:
         do 52 is=2,ns
         do 53 ix=1,3
         ue(ix)=0.0d0
         um(ix)=0.0d0
         do 53 jx=1,3
         ue(ix)=ue(ix)+A(ix,jx)*elelc(is-1,jx)
53       um(ix)=um(ix)+A(ix,jx)*elev( is-1,jx)
52       write(3,1002)is,(ue(ix),ix=1,3),(um(ix),ix=1,3)

         goto 43
98       close(8)
c
         deallocate(id,ics,ic,wc)
         if(imin.lt.0)deallocate(tl)
40       read(4,*)

         close(3)
         close(4)
         call report('DC.SC written')
         stop

        else
         write(6,*)prog
         call report('unknown source')
        endif
       endif
      endif

      if(.not.llen)call report(' Electric length dipoles not found')
      if(.not.lvel)call report(' Electric velocity dipoles not found')
      if(.not.lmag)call report(' Magnetic dipoles not found')
      if(.not.lene)call report(' Energies not found')
      if(.not.lgeo)call report(' Geometry not found')
c
      if(ivel.eq.1)then
       write(6,*)' velocity dipoles used'
       do 54 i=1,it
       do 54 ix=1,3
54     elel(i,ix)=elev(i,ix)
      else
       write(6,*)' length dipoles used'
      endif
c
c     transform aus to Bohr magnetons, debyes
      do 56 i=1,it
      do 56 ix=1,3
c     this is not sure, check the magnetic moment units!:
      mag(i,ix)=mag(i,ix)/2.0d0*2.541765d0
56    elel(i,ix)=elel(i,ix)*2.541765d0
      if(imag.eq.0)then
       do 55 i=1,it
       do 55 ix=1,3
55     mag(i,ix)=0.0d0
      endif
      open(3,file='DC.SC')
c     chromophore atoms:   H N C O
      read(4,*)ncs
      allocate (ic(ncs),wc(ncs),ics(ncs))
      backspace 4
      read(4,*)ncs,(ics(i),i=1,ncs)
      dco0=dsqrt((rn(1+3*(ics(4)-1))-rn(1+3*(ics(3)-1)))**2+
     1           (rn(2+3*(ics(4)-1))-rn(2+3*(ics(3)-1)))**2+
     1           (rn(3+3*(ics(4)-1))-rn(3+3*(ics(3)-1)))**2)
      dcn0=dsqrt((rn(1+3*(ics(2)-1))-rn(1+3*(ics(3)-1)))**2+
     1           (rn(2+3*(ics(2)-1))-rn(2+3*(ics(3)-1)))**2+
     1           (rn(3+3*(ics(2)-1))-rn(3+3*(ics(3)-1)))**2)
      dhn0=dsqrt((rn(1+3*(ics(2)-1))-rn(1+3*(ics(1)-1)))**2+
     1           (rn(2+3*(ics(2)-1))-rn(2+3*(ics(1)-1)))**2+
     1           (rn(3+3*(ics(2)-1))-rn(3+3*(ics(1)-1)))**2)


      read(4,*)NG
      write(3,*)NG,' groups'
      write(6,*)NG,' groups'
c
c     geometry and ES field corrections:
      read(4,*)(IP(i),i=1,npar0)
      write(6,6009)(IP(i),i=1,npar0)
6009  format(14i2)
      read(4,*)ntr
      do 6 it=1,ntr
6     read(4,*)(PG(it,i),i=1,npar0)
c
c     determine which reference value to use
      if(IP(1).lt.0)then
       write(6,*)' Ab initio value of d(C=O):'
      else
       dco0=1.231d0
       write(6,*)' MD value of d(C=O):'
      endif
      write(6,6070)dco0
 6070 format(f12.3,' A')
      if(IP(3).lt.0)then
       write(6,*)' Ab initio value of d(C=N):'
      else
       dcn0=1.335d0
       write(6,*)' MD value of d(C=N):'
      endif
      write(6,6070)dcn0
c
c     periodic boundary conditions
      read(4,*)iper,DPCX,DPCY,DPCZ
      lper=iper.eq.1
c
c     potential correction:
      if(IP(npar0).eq.1)then
       open(55,file='BIJ.SCR',status='old',form='unformatted')
       read(55)n4,n6,(zim(i),i=1,n4),(lfdp(i),i=1,4)
       write(6,*)'Zs:'
       write(6,6071)(zim(i),i=1,n4)
 6071  format(6f12.4)
       if(lfdp(1))read(55)
       if(lfdp(2))read(55)
       if(lfdp(3))read(55)
       if(lfdp(4))read(55)nwl,((buvcd(i,j),i=1,nwl),j=1,n6)
       close(55)
      endif

      do 2 IG=1,NG
      write(3,*)IG,' . group:'
      write(6,*)IG,' . group:'
c     Number of atoms determining center and their numbers:
      read(4,*)nc
      if(nc.ne.ncs)call report('nc<>ncs!')
      backspace 4
      read(4,*)nc,(ic(i),i=1,nc)
      write(6,*)nc,(ic(i),i=1,nc)
c     center weights:
      read(4,*)(wc(i),i=1,nc)
      call makecenter(rc ,R ,wc,nc,ic )
      call makecenter(rcs,rn,wc,nc,ics)
      write(3,1000)(rc(i)/10.0d0,i=1,3)
1000  format(3f15.6,' nm')
      
      ns=ntr+1
      write(3,*)ns,' states'
      write(3,*)
c
c     ground state:
      write(3,*)'0 0.0'
c
c     correct energies and dipoles according to
c     geometry and ES field:w=w0*(1+dd/d*p) dip=dip0*(1+dd/d*p)
c     first par-energy, second dipole
c      
      do 58 is=2,ns

      do 57 i=1,npar0
57    F(i)=0.0d0

      if(IP(1).ne.0)then
       dco=dsqrt((R(1+3*(ic(4)-1))-R(1+3*(ic(3)-1)))**2+
     1           (R(2+3*(ic(4)-1))-R(2+3*(ic(3)-1)))**2+
     1           (R(3+3*(ic(4)-1))-R(3+3*(ic(3)-1)))**2)
       F(1)=(dco-dco0)/dco0*PG(is-1,1)
       F(2)=(dco-dco0)/dco0*PG(is-1,2)
       if(is.eq.2)write(6,6001)1,dco,dco0
6001   format(i6,' d(C=O) = ',f10.5,'; D0 = ',f10.5)
      endif

c     C=N bond
      if(IP(3).ne.0)then
       dcn=dsqrt((R(1+3*(ic(2)-1))-R(1+3*(ic(3)-1)))**2+
     1           (R(2+3*(ic(2)-1))-R(2+3*(ic(3)-1)))**2+
     1           (R(3+3*(ic(2)-1))-R(3+3*(ic(3)-1)))**2)
       F(3)=(dcn-dcn0)/dcn0*PG(is-1,3)
       F(4)=(dcn-dcn0)/dcn0*PG(is-1,4)
       if(is.eq.2)write(6,6003)3,dcn,dcn0
6003   format(i6,' d(C-N) = ',f10.5,'; D0 = ',f10.5)
      endif

c     H-N bond
      if(IP(5).ne.0)then
       dhn=dsqrt((R(1+3*(ic(2)-1))-R(1+3*(ic(1)-1)))**2+
     1           (R(2+3*(ic(2)-1))-R(2+3*(ic(1)-1)))**2+
     1           (R(3+3*(ic(2)-1))-R(3+3*(ic(1)-1)))**2)
       F(5)=(dhn-dhn0)/dhn0*PG(is-1,5)
       F(6)=(dhn-dhn0)/dhn0*PG(is-1,6)
       if(is.eq.2)write(6,6004)5,dhn,dhn0
6004   format(i6,' d(N-H) = ',f10.5,'; D0 = ',f10.5)
      endif

      enc(is-1)=en(is-1)*(1.0d0+F(1)+F(3)+F(5))

      do 58 ix=1,3 
58    elelc(is-1,ix)=elel(is-1,ix)*(1.0d0+F(2)+F(4)+F(6))
      if(IP(npar0).ne.1)then
       call fw(n6,fi,N,ttype,R,ic(1),ic(2),ic(3),ic(4),lper,
     1 DPCX,DPCY,DPCZ)
       do 59 is=2,ns
       if(en(is-1).gt.200.0d0)then
        do 60 i=1,4
60      enc(is-1)=enc(is-1)+buvcd(1,i)*fi(i)
       else
        do 61 i=1,4
61      enc(is-1)=enc(is-1)+buvcd(2,i)*fi(i)
       endif
59     continue
      endif

      do 4 is=2,ns
      w=1.0d7/enc(is-1)
4     write(3,1001)is,w/8065.54093D0,w,enc(is-1)
1001  format(i3,f11.6,' eV =',f11.2,' cm-1 =',f11.2,' nm')
      write(3,*)
c
c     write group reference atoms into XS and corresponding ones in XB:
      if(ncs.gt.MXB)call report('ncs>MXB')
      do 5 i=1,ncs
      do 51 ix=1,3
      XB(i,ix)= R(ix+3*(ic(i)-1))-rc(ix)
51    XS(i,ix)=rn(ix+3*(ics(i)-1))-rcs(ix)
5     continue
      NATB=ncs
c
c     overlap XS,XB, make the rotation matrix, store in A:
      call domatrix(ierr)
      if(ierr.ne.0)call report('rotation not successful')
      write(6,*)' Transformation matrix'
      do 67673 i=1,3
67673 write(6,6767)(A(i,j),j=1,3)
6767  format(3f10.5)
c
c     write rotated dipole moments:
      do 2 is=2,ns

      do 8 ix=1,3
      ue(ix)=0.0d0
      um(ix)=0.0d0
      do 8 jx=1,3
      ue(ix)=ue(ix)+A(ix,jx)*elelc(is-1,jx)
8     um(ix)=um(ix)+A(ix,jx)* mag(is-1,jx)
c
2     write(3,1002)is,(ue(ix),ix=1,3),(um(ix),ix=1,3)
1002  format(i3,6g12.4)

      close(3)
      close(4)
      call report('DC.SC written')
      end
c
      SUBROUTINE REPORT(RS)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER*(*) RS
      WRITE(6,3000)
3000  FORMAT(/,80(1H*))
      WRITE(6,*)RS
      WRITE(6,3001)
3001  FORMAT(80(1H*),/,/,'PROGRAM STOPPED')
      STOP
      END

      subroutine rd(nt0,it,e,l)
      implicit none
      integer*4 nt0,it,idumm,i
      character*80 s80
      real*8 e(nt0,3)
      logical nu,l
      l=.true.
      it=0
      read(30,*)
31    read(30,300)s80
300   format(a80)
      if(nu(s80(10:10)))then
       it=it+1
       if(it.gt.nt0)call report(' too many transitions !')
       read(s80,*)idumm,(e(it,i),i=1,3)
       goto 31
      else
       backspace 30
       return
      endif
      end

      subroutine rdt(nt0,it,e,l)
      implicit none
      integer*4 nt0,it,idumm,i
      character*10 s10
      real*8 e(nt0,3)
      logical nu,l
      l=.true.
      it=0
      read(30,*)
31    read(30,300)s10
300   format(a10)
      if(nu(s10(4:4)))then
       it=it+1
       if(it.gt.nt0)then
        write(6,*)' too many transitions !'
        stop
       endif
       backspace 30
       read(30,*)idumm,(e(it,i),i=1,3)
       goto 31
      endif
      backspace 30
      return
      end


      subroutine ro(rn,nat,l)
      implicit none
      integer*4 dum,i,nat
      character*80 s80
      real*8 rn(*)
      logical nu,l

      l=.true.
      nat=0
      read(30,*)
      read(30,*)
      read(30,*)
      read(30,*)
31    read(30,300)s80
300   format(a80)
      if(nu(s80(5:5)).or.nu(s80(7:7)))then
       nat=nat+1
       read(s80,*)dum,dum,dum,(rn(i+3*(nat-1)),i=1,3)
       goto 31
      endif
      write(6,*)nat,' atoms in chromophore'
      return
      end

      subroutine rot(rn,nat,l)
      implicit none
      integer*4 i,nat
      character*80 s80
      character*15 dum
      real*8 rn(*)
      logical l
      l=.true.
      nat=0
      read(30,*)
31    read(30,300)s80
300   format(a80)
      if(s80(1:15).ne.'               ')then
       nat=nat+1
       read(s80,*)dum,(rn(i+3*(nat-1)),i=1,3)
       goto 31
      endif
      write(6,*)nat,' atoms in chromophore'
      return
      end

      function nu(c)
      character*1 c
      logical nu
      nu=c.eq.'1'.or.c.eq.'2'.or.c.eq.'3'.or.c.eq.'4'.or.c.eq.'5'.
     1or.c.eq.'6'.or.c.eq.'7'.or.c.eq.'8'.or.c.eq.'9'.or.c.eq.'0'
      return
      end

      SUBROUTINE DOMATRIX(IERR)
C     This is a Fortran version of Petr Malon's subroutine Amoeba
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MXB=50)
      DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3)
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER
      IERR=0
C     STARTING VALUES FOR THE ITERATION:
      ANG(1)=0.2d0
      ANG(2)=0.0d0
      ANG(3)=0.0d0
      DO 1 I=1,3
1     P(1,I)=ANG(I)
      Y(1)=FU(ANG)
      ANG(2)=ANG(1)
      ANG(1)=0.0d0
      DO 2 I=1,3
2     P(2,I)=ANG(I)
      Y(2)=FU(ANG)
      ANG(3)=ANG(2)
      ANG(2)=0.0d0
      DO 3 I=1,3
3     P(3,I)=ANG(I)
      Y(3)=FU(ANG)
      ANG(1)=ANG(3)
      ANG(2)=ANG(3)
      ANG(3)=0.0d0
      DO 4 I=1,3
4     P(4,I)=ANG(I)
      Y(4)=FU(ANG)
      FTOL=0.0000001d0
      ITMAX=1500
      ITER=0
99999 ILO=1
      IF(Y(1).GT.Y(2))THEN
       IHI=1
       INHI=2
      ELSE
       IHI=2
       INHI=1
      ENDIF
      DO 5 I=1,4
      IF(Y(I).LT.Y(ILO))ILO=I
      IF(Y(I).GT.Y(IHI))THEN
       INHI=IHI
       IHI=I
      ELSE
       IF(Y(I).GT.Y(INHI))THEN
        IF(I.NE.IHI)INHI=I
       ENDIF
      ENDIF
5     CONTINUE
      IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.0.00000000001d0)THEN
       RTOL=0.0d0
      ELSE
       RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0
      ENDIF
      IF(RTOL.LT.FTOL)RETURN
      IF(ITER.EQ.ITMAX)THEN
       WRITE(3,*)' Rotation has not converged !'
       IERR=1
       RETURN
      ENDIF
      ITER=ITER+1
      DO 55 I=1,3
55    PBAR(I)=0.0d0
      DO 6 I=1,4
      IF(I.NE.IHI)THEN
       DO 7 J=1,3
7      PBAR(J)=PBAR(J)+P(I,J)
      ENDIF
6     CONTINUE
      DO 8 J=1,3
      PBAR(J)=PBAR(J)/3.0d0
8     PR(J)=2.0d0*PBAR(J)-P(IHI,J)
      YPR=FU(PR)
      IF(YPR.LE.Y(ILO))THEN
       DO 9 J=1,3
9      PRR(J)=2.0d0*PR(J)-PBAR(J)
       YPRR=FU(PRR)
       IF(YPRR.LT.Y(ILO))THEN
        DO 10 J=1,3
10      P(IHI,J)=PRR(J)
        Y(IHI)=YPRR
       ELSE
        DO 11 J=1,3
11      P(IHI,J)=PR(J)
        Y(IHI)=YPR
       ENDIF
      ELSE
       IF(YPR.GE.Y(INHI))THEN
        IF(YPR.LT.Y(IHI))THEN
         DO 12 J=1,3
12       P(IHI,J)=PR(J)
         Y(IHI)=YPR
        ENDIF
        DO 13 J=1,3
13      PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J)
        YPRR=FU(PRR)
        IF(YPRR.LT.Y(IHI))THEN
         DO 14 J=1,3
14       P(IHI,J)=PRR(J)
         Y(IHI)=YPRR
        ELSE
         DO 15 I=1,4
         IF(I.NE.ILO)THEN
          DO 16 J=1,3
          PR(J)=0.5d0*(P(I,J)+P(ILO,J))
16        P(I,J)=PR(J)
          Y(I)=FU(PR)
         ENDIF
15       CONTINUE
        ENDIF
       ELSE
        DO 17 J=1,3
17      P(IHI,J)=PR(J)
        Y(IHI)=YPR
       ENDIF
      ENDIF
      GOTO 99999
      END

      FUNCTION FU(ANG)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (MXB=50)
      COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER
      DIMENSION ANG(3)
      S1=SIN(ANG(1))
      S2=SIN(ANG(2))
      S3=SIN(ANG(3))
      C1=COS(ANG(1))
      C2=COS(ANG(2))
      C3=COS(ANG(3))
C     ang1 = theta
C     ang2 = phi
C     ang3 = cappa
      A(1,1)=C1*C2*C3-S2*S3
      A(1,2)=C1*S2*C3+C2*S3
      A(1,3)=-S1*C3
      A(2,1)=-C1*C2*S3-S2*C3
      A(2,2)=-C1*S2*S3+C2*C3
      A(2,3)=S1*S3
      A(3,1)=S1*C2
      A(3,2)=S1*S2
      A(3,3)=C1
      R=0.0d0
      DO 1 I=1,NAT
      TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3)
      TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3)
      TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3)
1     R=R+(TX-XB(I,1))*(TX-XB(I,1))+
     1    (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3))
      FU=R
      RETURN
      END
c
      subroutine fw(n6,fi,nat,ity,r,iih,iin,iic,iio,lper,
     1DPCX,DPCY,DPCZ)
c     find waters and their potentials
      implicit none
      integer*4 i,n6,jj,ity(*),ityo,ityh,j1,nh,kk,k1,iih,iin,iic,iio,
     1ii,nat
      logical lper
      real*8 fi(*),xo,yo,zo,r(*),xh,yh,zh,d2,x1,y1,z1,
     1x2,y2,z2,qh,qo,DPCX,DPCY,DPCZ,dx2,dy2,dz2
      data ityo,ityh,qh,qo/649,650,0.41d0,-0.82d0/

      dx2=DPCX/2.0d0
      dy2=DPCY/2.0d0
      dz2=DPCZ/2.0d0
      do 4 i=1,n6
4     fi(i)=0.0d0
      do 6 jj=1,nat
      if(ity(jj).eq.ityo)then
       j1=3*(jj-1)
       xo=r(j1+1)
       yo=r(j1+2)
       zo=r(j1+3)
       nh=0
       do 7 kk=1,nat
       if(ity(kk).eq.ityh)then
        k1=3*(kk-1)
        xh=r(k1+1)
        yh=r(k1+2)
        zh=r(k1+3)
        d2=(xo-xh)**2+(yo-yh)**2+(zo-zh)**2
        if(d2.gt.0.64d0.and.d2.lt.1.44d0)then
         nh=nh+1
         if(nh.gt.2)call report('Too many Hs around one O')
         if(nh.eq.1)then
          x1=xh
          y1=yh
          z1=zh
         else 
          x2=xh
          y2=yh
          z2=zh
         endif
        endif
       endif
7      continue
       if(nh.eq.2)then
        call add(iih ,1,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
     1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
        call add(iin ,2,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
     1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
        call add(iic ,3,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
     1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
        call add(iio ,4,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
     1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
c       call add(iian,5,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
c    1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
c       call add(iiac,6,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
c    1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
       endif 
      endif
6     continue
      write(6,*)' Potentials:'
      write(6,1603)(fi(ii),ii=1,n6)
1603  format(8f9.3)
      return
      end
c
      subroutine add(ii,ic,r,xo,yo,zo,x1,y1,z1,x2,y2,z2,fi,qh,qo,
     1lper,dx2,dy2,dz2,DPCX,DPCY,DPCZ)
      implicit none
      real*8 r(*),xo,yo,zo,x1,y1,z1,x2,y2,z2,fi(*),x,y,z,ro,r1,r2,qh,qo,
     1dx2,dy2,dz2,DPCX,DPCY,DPCZ
      integer*4 ii,ic
      logical lper
c      
      x=r(3*(ii-1)+1)
      y=r(3*(ii-1)+2)
      z=r(3*(ii-1)+3)
      if(lper)then
       if(x-xo.gt.dx2.and.x-x1.gt.dx2.and.x-x2.gt.dx2)x=x-DPCX
       if(y-yo.gt.dy2.and.y-y1.gt.dy2.and.y-y2.gt.dy2)y=y-DPCY
       if(z-zo.gt.dz2.and.z-z1.gt.dz2.and.z-z2.gt.dz2)z=z-DPCZ
       if(xo-x.gt.dx2.and.x1-x.gt.dx2.and.x2-x.gt.dx2)x=x+DPCX
       if(yo-y.gt.dy2.and.y1-y.gt.dy2.and.y2-y.gt.dy2)y=y+DPCY
       if(zo-z.gt.dz2.and.z1-z.gt.dz2.and.z2-z.gt.dz2)z=z+DPCZ
      endif 
      ro=dsqrt((xo-x)**2+(yo-y)**2+(zo-z)**2)
      r1=dsqrt((x1-x)**2+(y1-y)**2+(z1-z)**2)
      r2=dsqrt((x2-x)**2+(y2-y)**2+(z2-z)**2)
c      
      fi(ic)=fi(ic)+qh/r1+qh/r2+qo/ro
      return
      end

      subroutine rvd(fn,nt0,nat,nm,rn,id,m,u,en)
      implicit none
      integer nat,nm,ia,id(*),ix,it,nt0
      real*8 rn(*),u(nt0,3),m(nt0,3),en(*)
      character*(*) fn
      open(30,file=fn,status='old')
      read(30,*)nm,nat,nat
      if(nm.gt.nt0)call report('Too many transitions')
      do 1 ia=1,nat
1     read(30,*)id(ia),(rn(ix+(ia-1)*3),ix=1,3)
      do 2 ia=1,2+nat*nm
2     read(30,*)
      read(30,300)(en(it),it=1,nm)
300   format(6f11.3)
      do 3 it=1,nm
3     read(30,*)(u(it,ix),ix=1,3),(m(it,ix),ix=1,3)
      close(30)
      return
      end

      subroutine makecenter(c,r,w,n,ic)
      implicit none
      integer*4 n,i,ii,ic(*)
      real*8 c(*),r(*),w(*),sum
      sum=0.0d0
      do 45 ii=1,n
45    sum=sum+w(ii)
      do 46 i=1,3
      c(i)=0.0d0
      do 47 ii=1,n
47    c(i)=c(i)+r(i+3*(ic(ii)-1))*w(ii)
46    c(i)=c(i)/sum
      return
      end