      subroutine redundant(icycle,natoms,nmodes,nbonds,nangles,
     &                     ndihedrals,nimpropers,ifrozen,natomtype,
     &                     i_bond,i_angle,i_dihedral,i_improper,ifreeze,
     &                     dmaxstep,grms,gmax,crms,cmax,coord,q_0,q_1,
     &                     deltaq,deltax,bmatrix,gmatrix,ginvmatrix,
     &                     gcart,gint,gintp,pmatrix,amatrix,hessian,
     &                     wkmat1,wkmat2,coordsave,wkvec)

c  this subroutine written by alain st-amant of the
c  department of chemistry, university of ottawa, ottawa, canada.
c  all rights reserved.  this is part of the DeFT project.

      implicit real*8(a-h,o-z)

      dimension natomtype(*),i_bond(2,*),i_angle(3,*),i_dihedral(4,*),
     &          i_improper(4,*),ifreeze(*)

      dimension coord(3,*),q_0(*),q_1(*),deltaq(*),deltax(*),
     &          bmatrix(nmodes,*),gmatrix(nmodes,*),
     &          ginvmatrix(nmodes,*),gcart(*),gint(*),gintp(*),
     &          pmatrix(nmodes,*),amatrix(3*natoms,*),hessian(nmodes,*),
     &          wkmat1(nmodes,*),wkmat2(nmodes,*),coordsave(3,*),
     &          wkvec(*)

      data zero,one,two,pi,bohr/0.0d0,1.0d0,
     &                          2.0d0,3.1415926535898d0,0.529177d0/
      data tolerance,big/0.000001d0,1000.0d0/
      data ang_cutoff,one80/170.0d0,180.0d0/

      fac=(one80/pi)

      call get_bmatrix(natoms,nmodes,nbonds,nangles,ndihedrals,
     &                 nimpropers,i_bond,i_angle,i_dihedral,i_improper,
     &                 coord,bmatrix)

      do 1001 i=1,nbonds
 1001 q_0(i)=dist(coord(1,i_bond(1,i)),coord(1,i_bond(2,i)))

      do 1002 i=1,nangles
      test=ang(coord(1,i_angle(1,i)),
     &         coord(1,i_angle(2,i)),
     &         coord(1,i_angle(3,i)))
      if(test.gt.ang_cutoff) icycle=0
      if(test.gt.ang_cutoff) return
      if(test.lt.(one80-ang_cutoff)) icycle=0
      if(test.lt.(one80-ang_cutoff)) return
 1002 q_0(nbonds+i)=dang(coord(1,i_angle(1,i)),
     &                   coord(1,i_angle(2,i)),
     &                   coord(1,i_angle(3,i)))

      do 1003 i=1,ndihedrals
      test1=ang(coord(1,i_dihedral(1,i)),
     &          coord(1,i_dihedral(2,i)),
     &          coord(1,i_dihedral(3,i)))
      test2=ang(coord(1,i_dihedral(2,i)),
     &          coord(1,i_dihedral(3,i)),
     &          coord(1,i_dihedral(4,i)))
      test=max(test1,test2)
      if(test.gt.ang_cutoff) icycle=0
      if(test.gt.ang_cutoff) return
      test=min(test1,test2)
      if(test.lt.(one80-ang_cutoff)) icycle=0
      if(test.lt.(one80-ang_cutoff)) return
 1003 q_0(nbonds+nangles+i)=ddihed(coord(1,i_dihedral(1,i)),
     &                             coord(1,i_dihedral(2,i)),
     &                             coord(1,i_dihedral(3,i)),
     &                             coord(1,i_dihedral(4,i)))

      do 1004 i=1,nimpropers
      test1=ang(coord(1,i_improper(1,i)),
     &          coord(1,i_improper(4,i)),
     &          coord(1,i_improper(2,i)))
      test2=ang(coord(1,i_improper(1,i)),
     &          coord(1,i_improper(4,i)),
     &          coord(1,i_improper(3,i)))
      test=max(test1,test2)
      if(test.gt.ang_cutoff) icycle=0
      if(test.gt.ang_cutoff) return
      test=min(test1,test2)
      if(test.lt.(one80-ang_cutoff)) icycle=0
      if(test.lt.(one80-ang_cutoff)) return
 1004 q_0(nbonds+nangles+ndihedrals+i)=dimp(coord(1,i_improper(1,i)),
     &                                      coord(1,i_improper(2,i)),
     &                                      coord(1,i_improper(3,i)),
     &                                      coord(1,i_improper(4,i)))

      do 1005 j=1,nmodes
      do 1005 i=1,nmodes
 1005 gmatrix(i,j)=zero

      do 1006 k=1,3*natoms
      do 1006 j=1,nmodes
      do 1006 i=1,nmodes
 1006 gmatrix(i,j)=gmatrix(i,j)+bmatrix(i,k)*bmatrix(j,k)

      ja=1
      jb=ja+nmodes*nmodes
      jc=jb+nmodes*nmodes
      jd=jc+nmodes

      call ginverse(nmodes,gmatrix,ginvmatrix,wkvec(ja),wkvec(jb),
     &              wkvec(jc),wkvec(jd))

      do 1007 i=1,nmodes
 1007 wkvec(i)=zero

      do 1008 j=1,3*natoms
      do 1008 i=1,nmodes
 1008 wkvec(i)=wkvec(i)+bmatrix(i,j)*gcart(j)

      do 1009 i=1,nmodes
 1009 gint(i)=zero

      do 1010 j=1,nmodes
      do 1010 i=1,nmodes
 1010 gint(i)=gint(i)+ginvmatrix(i,j)*wkvec(j)

      do 1011 j=1,nmodes
      do 1012 i=1,nmodes
 1012 pmatrix(i,j)=zero
      do 1011 k=1,nmodes
      do 1011 i=1,nmodes
 1011 pmatrix(i,j)=pmatrix(i,j)+gmatrix(i,k)*ginvmatrix(k,j)

      do 1013 i=1,nmodes
      wkvec(i)=zero
      do 1013 j=1,ifrozen
 1013 if(i.eq.ifreeze(j)) wkvec(i)=one

      do 1014 j=1,nmodes
      do 1015 i=1,nmodes
 1015 wkmat1(i,j)=zero
      do 1014 i=1,nmodes
 1014 wkmat1(i,j)=wkmat1(i,j)+wkvec(i)*pmatrix(i,j)

      do 1016 j=1,nmodes
      do 1017 i=1,nmodes
 1017 gmatrix(i,j)=zero
      do 1016 i=1,nmodes
 1016 gmatrix(i,j)=gmatrix(i,j)+pmatrix(i,j)*wkvec(j)

      do 1018 j=1,nmodes
      do 1019 i=1,nmodes
 1019 wkmat2(i,j)=zero
      do 1018 i=1,nmodes
 1018 wkmat2(i,j)=wkmat2(i,j)+wkvec(i)*gmatrix(i,j)

      ja=nmodes+1
      jb=ja+nmodes*nmodes
      jc=jb+nmodes*nmodes
      jd=jc+nmodes

      call ginverse(nmodes,wkmat2,gmatrix,wkvec(ja),wkvec(jb),
     &              wkvec(jc),wkvec(jd))

      do 1020 j=1,nmodes
      do 1021 i=1,nmodes
 1021 wkmat2(i,j)=zero
      do 1020 k=1,nmodes
      do 1020 i=1,nmodes
 1020 wkmat2(i,j)=wkmat2(i,j)+gmatrix(i,k)*wkmat1(k,j)

      do 1022 j=1,nmodes
      do 1023 i=1,nmodes
 1023 gmatrix(i,j)=zero
      do 1022 i=1,nmodes
 1022 gmatrix(i,j)=gmatrix(i,j)+wkvec(i)*wkmat2(i,j)

      do 1024 j=1,nmodes
      do 1025 i=1,nmodes
 1025 wkmat1(i,j)=zero
      do 1024 k=1,nmodes
      do 1024 i=1,nmodes
 1024 wkmat1(i,j)=wkmat1(i,j)+pmatrix(i,k)*gmatrix(k,j)

      do 1026 j=1,nmodes
      do 1026 i=1,nmodes
 1026 pmatrix(i,j)=pmatrix(i,j)-wkmat1(i,j)

      do 1027 i=1,nmodes
 1027 gintp(i)=zero

      do 1028 j=1,nmodes
      do 1028 i=1,nmodes
 1028 gintp(i)=gintp(i)+pmatrix(i,j)*gint(j)

      do 1029 j=1,nmodes
      do 1029 i=1,3*natoms
      amatrix(i,j)=zero
      do 1029 k=1,nmodes
 1029 amatrix(i,j)=amatrix(i,j)+bmatrix(k,i)*ginvmatrix(k,j)

      if(icycle.eq.1) then
                            call hguess(nmodes,nbonds,nangles,
     &                                  ndihedrals,nimpropers,natomtype,
     &                                  i_bond,i_angle,i_dihedral,
     &                                  i_improper,coord,hessian)
                      else
                            goto 1036
                      endif

      do 1030 j=1,nmodes
      do 1031 i=1,nmodes
 1031 wkmat1(i,j)=zero
      do 1030 k=1,nmodes
      do 1030 i=1,nmodes
 1030 wkmat1(i,j)=wkmat1(i,j)+hessian(i,k)*pmatrix(k,j)

      do 1032 j=1,nmodes
      do 1033 i=1,nmodes
 1033 hessian(i,j)=zero
      do 1032 k=1,nmodes
      do 1032 i=1,nmodes
 1032 hessian(i,j)=hessian(i,j)+pmatrix(i,k)*wkmat1(k,j)

      do 1034 j=1,nmodes
      do 1034 i=1,nmodes
      factor=zero
      if(i.eq.j) factor=one
 1034 hessian(i,j)=hessian(i,j)+big*(factor-pmatrix(i,j))

      ja=1
      jb=ja+nmodes*nmodes
      jc=jb+nmodes*nmodes
      jd=jc+nmodes

      call ginverse(nmodes,hessian,wkmat1,wkvec(ja),wkvec(jb),wkvec(jc),
     &              wkvec(jd))

      do 1035 j=1,nmodes
      do 1035 i=1,nmodes
 1035 hessian(i,j)=wkmat1(i,j)

 1036 continue

      ja=1
      jb=ja+nmodes
      jc=jb+nmodes
      jd=jc+nmodes
      je=jd+nmodes
      jf=je+nmodes*nmodes

      if(icycle.gt.1) then
                            call bfgs_r(icycle,nmodes,hessian,gintp,q_0,
     &                                  wkvec(ja),wkvec(jb),wkvec(jc),
     &                                  wkvec(jd),wkvec(je),wkvec(jf))
                      endif

      do 1037 i=1,nmodes
 1037 deltaq(i)=zero

      do 1038 j=1,nmodes
      do 1038 i=1,nmodes
 1038 deltaq(i)=deltaq(i)-hessian(i,j)*gintp(j)

      do 1039 i=1,nmodes
 1039 wkvec(i)=deltaq(i)

      do 1040 i=1,natoms
      do 1040 j=1,3
 1040 coordsave(j,i)=coord(j,i)

      check_save=one/tolerance

 1041 continue

      do 1042 i=1,3*natoms
 1042 deltax(i)=zero

      do 1043 j=1,nmodes
      do 1043 i=1,3*natoms
 1043 deltax(i)=deltax(i)+amatrix(i,j)*wkvec(j)

      check=zero

      do 1044 i=1,3*natoms
 1044 check=check+deltax(i)**2

      check=sqrt(check/real64(3*natoms))

      if(check.gt.check_save) write(6,1045)
 1045 format(/,' divergence in redundant internals --- resetting ',/)
      if(check.gt.check_save) goto 1054

      check_save=check

      do 1046 i=1,natoms
      do 1046 j=1,3
 1046 coord(j,i)=coord(j,i)+deltax(3*(i-1)+j)

      do 1047 i=1,nbonds
 1047 q_1(i)=dist(coord(1,i_bond(1,i)),coord(1,i_bond(2,i)))

      do 1048 i=1,nangles
 1048 q_1(nbonds+i)=dang(coord(1,i_angle(1,i)),
     &                   coord(1,i_angle(2,i)),
     &                   coord(1,i_angle(3,i)))

      do 1049 i=1,ndihedrals
 1049 q_1(nbonds+nangles+i)=ddihed(coord(1,i_dihedral(1,i)),
     &                             coord(1,i_dihedral(2,i)),
     &                             coord(1,i_dihedral(3,i)),
     &                             coord(1,i_dihedral(4,i)))

      do 1050 i=1,nimpropers
 1050 q_1(nbonds+nangles+ndihedrals+i)=dimp(coord(1,i_improper(1,i)),
     &                                      coord(1,i_improper(2,i)),
     &                                      coord(1,i_improper(3,i)),
     &                                      coord(1,i_improper(4,i)))

      do 1051 i=1,nmodes
 1051 wkvec(i)=q_1(i)-q_0(i)

      do 1052 i=1,nmodes
      if(wkvec(i).lt.(-pi)) wkvec(i)=wkvec(i)+two*pi
 1052 if(wkvec(i).gt.(+pi)) wkvec(i)=wkvec(i)-two*pi

      do 1053 i=1,nmodes
 1053 wkvec(i)=deltaq(i)-wkvec(i)

      if(check.gt.tolerance) goto 1041

 1054 continue

      cmax=zero

      do 1055 i=1,natoms
      difference=sqrt((coordsave(1,i)-coord(1,i))**2
     &               +(coordsave(2,i)-coord(2,i))**2
     &               +(coordsave(3,i)-coord(3,i))**2)
 1055 if(difference.gt.cmax) cmax=difference

      if(cmax.gt.dmaxstep) then

              do 1056 i=1,natoms
              do 1056 j=1,3
 1056         coord(j,i)=coordsave(j,i)
     &                  +(dmaxstep/cmax)*(coord(j,i)-coordsave(j,i))

                           endif

      do 1057 i=1,nmodes
 1057 wkvec(i)=q_1(i)-q_0(i)

      do 1058 i=1,nmodes
      if(wkvec(i).lt.(-pi)) wkvec(i)=wkvec(i)+two*pi
 1058 if(wkvec(i).gt.(+pi)) wkvec(i)=wkvec(i)-two*pi

      do 1059 i=1,nmodes
 1059 wkvec(i)=zero-wkvec(i)

 1060 continue

      do 1061 i=1,3*natoms
 1061 deltax(i)=zero

      do 1062 j=1,nmodes
      factor=zero
      do 1063 k=1,ifrozen
 1063 if(j.eq.ifreeze(k)) factor=one
      do 1062 i=1,3*natoms
 1062 deltax(i)=deltax(i)+amatrix(i,j)*wkvec(j)*factor

      do 1064 i=1,natoms
      do 1064 j=1,3
 1064 coord(j,i)=coord(j,i)+deltax(3*(i-1)+j)

      do 1065 i=1,nbonds
 1065 q_1(i)=dist(coord(1,i_bond(1,i)),coord(1,i_bond(2,i)))

      do 1066 i=1,nangles
 1066 q_1(nbonds+i)=dang(coord(1,i_angle(1,i)),
     &                   coord(1,i_angle(2,i)),
     &                   coord(1,i_angle(3,i)))

      do 1067 i=1,ndihedrals
 1067 q_1(nbonds+nangles+i)=ddihed(coord(1,i_dihedral(1,i)),
     &                             coord(1,i_dihedral(2,i)),
     &                             coord(1,i_dihedral(3,i)),
     &                             coord(1,i_dihedral(4,i)))

      do 1068 i=1,nimpropers
 1068 q_1(nbonds+nangles+ndihedrals+i)=dimp(coord(1,i_improper(1,i)),
     &                                      coord(1,i_improper(2,i)),
     &                                      coord(1,i_improper(3,i)),
     &                                      coord(1,i_improper(4,i)))

      do 1069 i=1,nmodes
 1069 wkvec(i)=q_1(i)-q_0(i)

      do 1070 i=1,nmodes
      if(wkvec(i).lt.(-pi)) wkvec(i)=wkvec(i)+two*pi
 1070 if(wkvec(i).gt.(+pi)) wkvec(i)=wkvec(i)-two*pi

      do 1071 i=1,nmodes
 1071 wkvec(i)=zero-wkvec(i)

      check=zero

      do 1072 i=1,3*natoms
 1072 check=check+deltax(i)**2

      check=sqrt(check/real64(3*natoms))

      if(check.gt.tolerance**2) goto 1060

      rewind 21

      call fastwr(21,hessian,nmodes*nmodes)
      call fastwr(21,gintp,nmodes)
      call fastwr(21,q_0,nmodes)

      grms=zero

      do 1073 i=1,nmodes
 1073 grms=grms+gintp(i)*gintp(i)

      grms=sqrt(grms/real64(nmodes))

      gmax=zero

      do 1074 i=1,nmodes
 1074 if(abs(gintp(i)).gt.gmax) gmax=abs(gintp(i))

      crms=zero

      do 1075 i=1,natoms
 1075 crms=crms+(coordsave(1,i)-coord(1,i))**2
     &         +(coordsave(2,i)-coord(2,i))**2
     &         +(coordsave(3,i)-coord(3,i))**2

      crms=sqrt(crms/real64(natoms))

      cmax=zero

      do 1076 i=1,natoms
      difference=sqrt((coordsave(1,i)-coord(1,i))**2
     &               +(coordsave(2,i)-coord(2,i))**2
     &               +(coordsave(3,i)-coord(3,i))**2)
 1076 if(difference.gt.cmax) cmax=difference

      write(6,1077)
 1077 format(//,' set of redundant internal coordinates:')

      write(6,1078)
 1078 format(/,'      bonds:',/)

      do 1079 i=1,nbonds
 1079 write(6,1080) i,i_bond(1,i),i_bond(2,i),q_1(i)*bohr
 1080 format('            bond #',i4,' : ',2i5,f15.5)

      write(6,1081)
 1081 format(/,'      angles:',/)
 
      do 1082 i=1,nangles
      ii=nbonds+i
 1082 write(6,1083) i,i_angle(1,i),i_angle(2,i),i_angle(3,i),q_1(ii)*fac
 1083 format('            angle #',i4,' : ',3i5,f15.3)

      write(6,1084)
 1084 format(/,'      dihedrals:',/)
 
      do 1085 i=1,ndihedrals
      ii=nbonds+nangles+i
 1085 write(6,1086) i,i_dihedral(1,i),i_dihedral(2,i),
     &                i_dihedral(3,i),i_dihedral(4,i),q_1(ii)*fac
 1086 format('            dihedral #',i4,' : ',4i5,f15.3)

      write(6,1087)
 1087 format(/,'      improper torsions:',/)
 
      do 1088 i=1,nimpropers
      ii=nbonds+nangles+ndihedrals+i
 1088 write(6,1089) i,i_improper(1,i),i_improper(2,i),
     &                i_improper(3,i),i_improper(4,i),q_1(ii)*fac
 1089 format('            improper torsion #',i4,' : ',4i5,f15.3)

      write(6,1090)
 1090 format(//)

      return
      end
