C 30 APR 07 - SPW,TDI - NEW MODULE FOR NEO GEOMETRY SEARCHES
C
C*MODULE NEOOPT  *DECK NEOSPT
      SUBROUTINE NEOSPT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /OPTGRD/ X(3*MXATM),ES,FE(20),
     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGH
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA OPTMIZ /8HOPTIMIZE/
C
      NSERCH = 0
      CALL UNFREZ(.FALSE.)
      RUNSV = RUNTYP
      RUNTYP = OPTMIZ
      IF(RUNTYP.EQ.OPTMIZ) CALL SIGX(.FALSE.)
      IF(RUNTYP.NE.OPTMIZ) THEN
        WRITE(*,*) 'PROBLEMS CHANGIN THE RUNTYP IN NEOSPT'
        CALL ABRT
      END IF
      RUNTYP = RUNSV
C
      RETURN
      END
C*MODULE NEOOPT  *DECK NEOGOP
      SUBROUTINE NEOGOP
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL NEOLOOP,NUCOPT,CLSTEP,NEOTS,CVGDQM,CVGDCL
C
      PARAMETER (MXATM=2000)
      PARAMETER (ONEPT5=1.5D+00, THREE=3.00D+00)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NEOOPT/ NUCOPT,CLSTEP,NEOTS
      COMMON /OPTGRD/ X(3*MXATM),ES,FE(20),
     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGH
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /TMVALS/ TI,TX,TIM
C
      NSERCH_CL = 0
      NSERCH_QM = 0
      NEOLOOP = .TRUE.
C
      DO WHILE (NEOLOOP)
C
C        --- OPTIMIZE THE QM. NUCLEAR CENTERS ---
C
         NSERCH = 0
         CLSTEP = .FALSE.
         CALL UNFREZ(.FALSE.)
         CALL SIGX(.FALSE.)
         CVGDQM = (FRMS .LT. CONVF/THREE) .AND. (FMAX .LT. CONVF)
         NSERCH_QM = NSERCH
C
C        --- OPTIMIZE THE CLASSICAL NUCLEI  ---
C
         NSERCH = 0
         CLSTEP = .TRUE.
C
         CALL UNFREZ(.TRUE.)
         CALL SIGX(.FALSE.)
C
         NSERCH_CL = NSERCH_CL + 1
         CVGDCL= (FRMS .LT. CONVF/THREE) .AND. (FMAX .LT. CONVF)
         IF(CVGDQM .AND. CVGDCL .AND. NSERCH_QM.EQ.0) THEN
            NEOLOOP = .FALSE.
            WRITE(IW,9010)
         ELSE
            WRITE(IW,9020) NSERCH_CL
         END IF
C
C        --- CHECK TIME ---
C
         TNEED = ONEPT5 * (TIM/NSERCH_CL)
         IF((TIMLIM-TIM) .LT. TNEED) THEN
            WRITE(IW,9030) NSERCH_CL,TIM,TIMLIM,TNEED
            GO TO 100
         END IF
C
      END DO
C
  100 CONTINUE
C
      RETURN
 9010 FORMAT(/1X,'... NEO GEOMETRY OPTIMIZATION IS NOW COMPLETE ...'/)
 9020 FORMAT(/5X,57(1H-)/,5X,'CLASSICAL STEP ',I4,' DONE - NOW ',
     *       'OPTIMIZE QM NUCLEI CENTERS'/5X,57(1H-),/)
 9030 FORMAT(6X,'***** FAILURE TO LOCATE STATIONARY POINT,',
     *           ' NSERCH_CL=',I5,' *****'/
     *           1X,'TOO LITTLE TIME TO DO ANOTHER POINT,',
     *           1X,'TIME USED=',F15.1/' TIMLIM GIVEN=',F15.1,
     *              ' NEXT POINT NEEDS=',F15.1,' SECONDS')
      END
C
C*MODULE NEOOPT  *DECK NEOSTP
      SUBROUTINE NEOSTP
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL NUCOPT,CLSTEP,NEOTS
C
      PARAMETER (MXATM=2000)
C
      COMMON /NEOOPT/ NUCOPT,CLSTEP,NEOTS
      COMMON /OPTGRD/ X(3*MXATM),ES,FE(20),
     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGH
C
C    --- CLSTEP IS THE CLASS. NUC. OPTIMIZATION STEP ---
C
      NSTSV = NSTEP
      IF (CLSTEP) THEN
         NSTSV = NSTEP
         NSTEP = 0
      ELSE
         NSTEP = NSTSV
      END IF
C
      RETURN
      END
C
CspwC*MODULE NEOOPT  *DECK NEOSAD
Cspw      SUBROUTINE NEOSAD
CspwC
Cspw      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
CspwC
Cspw      LOGICAL NEOLOOP,NUCOPT,CLSTEP,NEOTS,NEOHSS,CVGDQM,CVGDCL
CspwC
Cspw      PARAMETER (MXATM=2000)
Cspw      PARAMETER (ONEPT5=1.5D+00, THREE=3.00D+00)
CspwC
CspwC     DOUBLE PRECISION CHECK,RUNTYP,EXETYP
CspwC
Cspw      COMMON /FRZCRT/ IFCART(3*MXATM),NFCART
Cspw      COMMON /FROZNC/ IFREEZ(3*MXATM),NFRZ
Cspw      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
Cspw     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
Cspw      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
Cspw      COMMON /NEOOPT/ NUCOPT,CLSTEP,NEOTS
Cspw      COMMON /NUCHSS/ HSSINI,HSSUPD,NEOHSS
Cspw      COMMON /OPTGRD/ X(3*MXATM),ES,FE(20),
Cspw     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
Cspw     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
Cspw     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
Cspw     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGH
Cspw      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
Cspw      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
Cspw      COMMON /TMVALS/ TI,TX,TIM
CspwC
Cspw      DATA SADPT/8HSADPOINT/,OPT     /8HOPTIMIZE/,CALCUP/8HCALCUP  /
Cspw      DATA CALCH /8HCALCH   /
CspwC
CspwC     --- NEOSAD IS NOT CURRENTLY USED - IT MAYBE HOOKED UP AT SOME POINT
CspwC         SPW ---
CspwC
Cspw      NCOORD = NAT*3
Cspw      NSERCH_NU = 0
CspwC
CspwC     --- CALCULATE THE INITIAL HESSIAN IF REQUESTED ---
CspwC
Cspw      IF(HSSINI.EQ.CALCH) THEN
CspwC
CspwC        --- OPTIMIZE THE QM. NUC. CENTERS---
CspwC
CspwC        WRITE(*,*) '** RUNNING NEOSAD OPTIMIZE THE QM.N.B.F.C. **'
Cspw         WRITE(IW,9140)
Cspw         NEOTS  = .FALSE.
Cspw         CLSTEP = .FALSE.
Cspw         RUNTYP = OPT
CspwC
Cspw         IF(RUNTYP.EQ.OPT) THEN
Cspw           CALL UNFREZ(.FALSE.)
Cspw           CALL SIGX(.FALSE.)
Cspw         ELSE
CspwC          WRITE(*,*) 'PROBLEMS CHANGIN THE RUNTYP IN NEOTS'
Cspw           WRITE(IW,9190)
Cspw           CALL ABRT
Cspw         END IF
Cspw         CVGDQM = (FRMS .LT. CONVF/THREE) .AND. (FMAX .LT. CONVF)
CspwC
CspwC        --- CALCULATE THE HESSIAN ---
CspwC
Cspw         NEOTS  = .FALSE.
Cspw         CLSTEP = .FALSE.
Cspw         CALL NEOSTP
CspwC        WRITE(*,*) '**RUNNING NEOSAD CALCULATING THE HESSIAN**'
Cspw         WRITE(IW,9150)
Cspw         CALL VCLR(IFCART,1,NCOORD)
Cspw         CALL VCLR(IFREEZ,1,NCOORD)
CspwC
Cspw         CALL HESSX(.FALSE.,.FALSE.)
Cspw      END IF
CspwC
CspwC     --- START NEO TS SEARCH ---
CspwC
Cspw      NSERCH = 0
Cspw      NEOLOOP = .TRUE.
CspwC
Cspw      DO WHILE (NEOLOOP)
CspwC
CspwC     --- OPTIMIZE THE QM. NUCLEI CENTERS ---  !!SKIP THIS STEP IF HESS=CALC
CspwC
CspwC        WRITE(*,*) '**RUNNING NEOSAD OPTIMIZE THE QM NUCLEAR CENTERS**'
Cspw         WRITE(IW,9160)
Cspw         NEOTS  = .FALSE.
Cspw         CLSTEP = .FALSE.
Cspw         RUNTYP = OPT
CspwC
Cspw         IF(RUNTYP.EQ.OPT) THEN
Cspw           CALL UNFREZ(.FALSE.)
Cspw           CALL SIGX(.FALSE.)
Cspw         ELSE
CspwC          WRITE(*,*) 'PROBLEMS CHANGIN THE RUNTYP IN NEOTS'
Cspw           WRITE(IW,9190)
Cspw           CALL ABRT
Cspw         END IF
CspwC
Cspw         CVGDQM = (FRMS .LT. CONVF/THREE) .AND. (FMAX .LT. CONVF)
Cspw         NSERCH_QM = NSERCH
Cspw         WRITE(*,*) 'NSERCH_QM=',NSERCH_QM
CspwC
CspwC        --- CALCULATE THE HESSIAN ---
CspwC
Cspw         IF(HSSUPD.EQ.CALCUP .AND. NSERCH_NU.GT.0) THEN
CspwC
Cspw            NEOTS  = .FALSE.
Cspw            CLSTEP = .FALSE.
Cspw            CALL NEOSTP
CspwC           WRITE(*,*) '**RUNNING NEOSAD CALCULATING THE HESSIAN**'
Cspw            WRITE(IW,9170)
Cspw            CALL VCLR(IFCART,1,NCOORD)
Cspw            CALL VCLR(IFREEZ,1,NCOORD)
CspwC
Cspw            CALL HESSX(.FALSE.,.FALSE.)
CspwC
CspwC           ---  SADDLE POINT SEARCH FOR THE CLASSICAL NUCLEI  ---
CspwC
CspwC           ---  GET THE VIB0 GEOMETRY --- !! THAT MUST BE DONE IN HESS.SRC
CspwC
Cspw            CALL VIBZRO(.FALSE.)
CspwC
Cspw         END IF
CspwC
CspwC        --- CORRECT THE OPTIMIZATION STEP NSERCH
CspwC
Cspw         NSERCH = NSERCH_NU
Cspw         NEOTS  = .TRUE.
Cspw         CLSTEP = .TRUE.
CspwC
Cspw         CALL NEOSTP
CspwC
CspwC        WRITE(*,*) '**RUNNING NEOSAD SEARCH FOR THE CLASSICAL NUCLEI**'
Cspw         WRITE(IW,9180)
Cspw         RUNTYP = SADPT
CspwC
Cspw         IF(RUNTYP.EQ.SADPT) THEN
Cspw            CLSTEP = .TRUE.
Cspw            NEOTS = .TRUE.
Cspw            CALL UNFREZ(.TRUE.)
Cspw            CALL SIGX(.FALSE.)
Cspw            NSERCH_NU = NSERCH_NU + 1
Cspw            NSERCH = 0
Cspw         ELSE
CspwC           WRITE(*,*) 'PROBLEMS CHANGIN THE RUNTYP IN NEOTS'
Cspw            WRITE(IW,9190)
Cspw            CALL ABRT
Cspw         END IF
CspwC
Cspw         CVGDCL= (FRMS .LT. CONVF/THREE) .AND. (FMAX .LT. CONVF)
CspwC
Cspw         IF(CVGDQM .AND. CVGDCL.AND. NSERCH_QM.EQ.0) THEN
Cspw            NEOLOOP=.FALSE.
CspwC           WRITE(*,*) '****NEO SADDLE POINT SEARCH IS COMPLETED****'
Cspw            WRITE(IW,9210)
Cspw         ELSE
CspwC           WRITE(*,*) '****NEO SADDLE POINT SEARCH WILL CONTINUE****'
Cspw            WRITE(IW,9200)
Cspw         END IF
CspwC
CspwC        --- CHECK TIME LIMIT ---
CspwC
Cspw         TNEED = ONEPT5 * (TIM/NSERCH_NU)
Cspw         WRITE(*,*) 'ONEPT5=',ONEPT5
Cspw         WRITE(*,*) 'TIM=',TIM
Cspw         WRITE(*,*) 'NSERCH_NU=',NSERCH_NU
Cspw         WRITE(*,*) 'TIMLIM-TIM=',TIMLIM-TIM
Cspw         WRITE(*,*) 'TNEED =',TNEED
CspwC
Cspw         IF((TIMLIM-TIM) .LT. TNEED) THEN
Cspw            WRITE(IW,9220) NSERCH_NU,TIM,TIMLIM,TNEED
Cspw            GO TO 100
Cspw         END IF
Cspw      END DO
CspwC
Cspw  100 CONTINUE
CspwC
Cspw      RETURN
Cspw 9140 FORMAT(6X,'RUNNING NEOSAD INITIAL OPTIMIZATION OF QM.N.B.F.C')
Cspw 9150 FORMAT(6X,'RUNNING NEOSAD CALCULATING INITIAL HESSIAN       ')
Cspw 9160 FORMAT(6X,'RUNNING NEOSAD OPTIMIZE THE QM NUCLEAR CENTERS   ')
Cspw 9170 FORMAT(6X,'RUNNING NEOSAD CALCULATING THE NEO-HESSIAN       ')
Cspw 9180 FORMAT(6X,'RUNNING NEOSAD SEARCH FOR THE CLASSICAL NUCLEI   ')
Cspw 9190 FORMAT(6X,'***** PROBLEMS CHANGIN THE RUNTYP IN NEOSAD *****')
Cspw 9200 FORMAT(6X,'***** NEO SADDLE POINT SEARCH WILL CONTINUE *****')
Cspw 9210 FORMAT(6X,'***** NEO SADDLE POINT SEARCH IS COMPLETED  *****')
Cspw 9220 FORMAT(6X,'***** FAILURE TO LOCATE STATIONARY POINT,',
Cspw     *           ' NSERCH_NU=',I5,' *****'/
Cspw     *           1X,'TOO LITTLE TIME TO DO ANOTHER POINT,',
Cspw     *           1X,'TIME USED=',F15.1/' TIMLIM GIVEN=',F15.1,
Cspw     *              ' NEXT POINT NEEDS=',F15.1,' SECONDS')
Cspw      END
CspwC
CspwC*MODULE NEOOPT  *DECK BKRUPM_NEO
Cspw      SUBROUTINE BKRUPM_NEO(NPUN,OUT)
CspwC
CspwCSPW  THIS SUBROUTINE IS NOT CALLED AT ALL.
CspwCSPW  I LEFT IT HERE AS IT MAYBE NEEDED IF NUCOPT=.T. SADPOINT
CspwCSPW  IS IMPLEMENTED IN THE FUTURE.
CspwC
Cspw      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CspwC
Cspw      LOGICAL OUT
Cspw      LOGICAL NUCOPT,CLSTEP,NEOTS
CspwC
Cspw      PARAMETER (MXATM=2000)
Cspw      PARAMETER (MXFRG=1050,
Cspw     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
CspwC
Cspw      COMMON /FMCOM / XX(1)
Cspw      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
Cspw     *                NPPTS(MXFRG),NPTTPT,IEFP,
Cspw     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
Cspw     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
Cspw      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
Cspw     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
Cspw      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
Cspw     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
Cspw     *                NNA,NNB,NTAUXB
Cspw      COMMON /NEOOPT/ NUCOPT,CLSTEP,NEOTS
CspwC
CspwC     ----- ALLOCATE MEMORY -----
CspwC          -FC-      AT I10 (NEO CLFCM/MIGHT BE QMFCM)
CspwC          -FCCART-  AT I20
CspwC          -SVEC-    AT I30
CspwC          -TVEC-    AT I40
CspwC
CspwC--TDI
Cspw      WRITE(*,*) 'IN BKRUPM_NEO NEOTS,CLSTEP',NEOTS,CLSTEP
CspwC     IF(NEOTS .AND. CLSTEP) THEN
Cspw        NCL3=(NAT-NQMNUC-NTAUXB)*3
Cspw        NC1 = NCL3
CspwC     ELSE IF (NEOTS .AND. .NOT.CLSTEP)
CspwC       NQM3=(NQMNUC+NTAUXB)*3
CspwC       NC1 = NQM3
CspwC     END IF
CspwC--TDI
Cspw      NC3 = NC1*NC1
Cspw      NCOORD = 3*NAT
CspwC
Cspw      CALL VALFM(LOADFM)
Cspw      I10 = 1 + LOADFM
Cspw      I20 = I10 + NC1
Cspw      I30 = I20 + NC3
Cspw      I40 = I30 + (NCOORD+6*NFRG)**2
Cspw      I50 = I40 + NC1
Cspw      LAST= I50 + NC1
Cspw      NEED = LAST - LOADFM - 1
Cspw      CALL GETFM(NEED)
CspwC
Cspw      CALL BKRUPD_NEO(NC1,NCOORD,NPUN,OUT,
Cspw     *            XX(I10),XX(I20),XX(I30),XX(I40),XX(I50))
CspwC
Cspw      CALL RETFM(NEED)
CspwC
Cspw      RETURN
Cspw      END
CspwC
CspwC*MODULE NEOOPT  *DECK BKRUPD_NEO
Cspw      SUBROUTINE BKRUPD_NEO(NCVAL,NCOORD,NPUN,OUT,
Cspw     *                  F,FC,FCCART,SVEC,TVEC)
CspwC
Cspw      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CspwC
Cspw      LOGICAL LINEAR,OUT,GOPARR,MASWRK,DSKWRK
Cspw      LOGICAL NUCOPT,CLSTEP,NEOTS
CspwC
CspwC
Cspw      PARAMETER (MXATM=2000)
Cspw      PARAMETER (MXFRG=1050, MXFGPT=12000)
Cspw      PARAMETER (MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
Cspw      PARAMETER (ZERO=0.0D+00)
CspwC
Cspw      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
Cspw     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
Cspw     *                ATORQ(3,MXFRG)
Cspw      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
Cspw     *                NPPTS(MXFRG),NPTTPT,IEFP,
Cspw     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
Cspw     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
Cspw      COMMON /FUNCT / ENERGY,EG(3*MXATM)
Cspw      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
Cspw      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
Cspw     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
Cspw      COMMON /NEOOPT/ NUCOPT,CLSTEP,NEOTS
Cspw      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
Cspw     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
Cspw     *                NNA,NNB,NTAUXB
Cspw      COMMON /OPTEF / D(3*MXATM+6*MXFRG),OLDF(3*MXATM+6*MXFRG),
Cspw     *                VMODE(3*MXATM+6*MXFRG),RADIUS
Cspw      COMMON /OPTGRD/ X(3*MXATM),ES,FE(20),
Cspw     *                CONVF,FMAXT,DXMAXT,RMAX,RMIN,RLIM,
Cspw     *                EIGMAX,EIGMIN,GRDERR,FRMS,FMAX,TRMAX,TRMIN,
Cspw     *                IC(20),NSTEP,NSERCH,NPMAX,NP,IFOLOW,
Cspw     *                NNEG,IUPHSS,IEXIT,ITRUPD,IPAD,KDIAGH
Cspw      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
Cspw      COMMON /SIMDAT/ NACC,NREJ,IGOMIN,NRPA,IBWM,NACCT,NREJT,NRPAT,
Cspw     *                NPRTGO,IDPUNC,IGOFLG
Cspw      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
CspwC
Cspw      DIMENSION F(*),FC(NCVAL+6*NFRG,NCVAL+6*NFRG),
Cspw     *          SVEC(NCVAL+6*NFRG),TVEC(NCVAL+6*NFRG),
Cspw     *          FCCART(NCOORD+6*NFRG,NCOORD+6*NFRG)
CspwC
CspwCSPW  THIS SUBROUTINE IS NOT CALLED AT ALL.
CspwCSPW  I LEFT IT HERE AS IT MAYBE NEEDED IF NUCOPT=.T. SADPOINT
CspwCSPW  IS IMPLEMENTED IN THE FUTURE.
CspwC
CspwC     ----- UPDATE THE HESSIAN MATRIX -----
CspwC     2 UPDATING PROCEDURES ARE AVAILABLE, DEPENDING ON UPHESS.
CspwC
CspwC     UPDATE USES CURRENT GRADIENT, GRADIENT FROM PREVIOUS POINT
CspwC     AND CORRECTION VECTOR USED ON LAST CYCLE
CspwC     SVEC AND TVEC ARE FOR TEMPORARY STORAGE
CspwC     D    -- CORRECTION VECTOR (H)
CspwC     SVEC -- V = G'-G
CspwC     TVEC -- H(OLD)H
CspwC
CspwC     NC1 = NCVAL+6*NFRG
CspwC     NC3 = NC1*NC1
CspwC--TDI
Cspw      NCL3=(NAT-NQMNUC-NTAUXB)*3
Cspw      NC1 = NCL3
Cspw      NC3 = NC1*NC1
CspwC--TDI
CspwC
CspwC     ----- READ IN FC MATRIX ------
CspwC
CspwC--TDI
Cspw      IF(NZVAR.EQ.0 .AND. CLSTEP .AND. NEOTS) THEN
CspwCSPW CHECK THIS
Cspw         CALL DAREAD(IDAF,IODA,FC,NC3, 43,0)
Cspw      ELSE
Cspw         CALL DAREAD(IDAF,IODA,FC,NC3,43,0)
Cspw      END IF
Cspw      IF(OUT) WRITE(IW,9020)
Cspw      IF(OUT) CALL PRSQ(FC,NC1,NC1,NC1)
CspwC--TDI
Cspw      WRITE(*,*) 'THE CLFCM :-) IN BKRUPD_NEO'
Cspw      CALL FCMOUT(FC,NCL3)
CspwC
Cspw      DO 120 I=1,NC1
Cspw         TVEC(I)=ZERO
Cspw         DO 110 J=1,NC1
Cspw            TVEC(I)=TVEC(I) + FC(I,J)*D(J)
Cspw  110    CONTINUE
Cspw  120 CONTINUE
CspwC
Cspw         CALL EGMOVE(F,EG,DEFT,TORQ,NCVAL,NFRG)
CspwC
CspwC     ----- POWELL UPDATE (USUAL FOR RUNTYP=SADPOINT) -----
CspwC     ----- MIXED MS-POWELL UPDATE (SEE BOFILL) -----
CspwC
Cspw      IF(IUPHSS.EQ.1  .OR.  IUPHSS.EQ.5) THEN
Cspw         IF (MASWRK  .AND.  NPRTGO.GT.0) WRITE(IW,9000)
Cspw         DO 200 I=1,NC1
Cspw            TVEC(I)=F(I)-OLDF(I)-TVEC(I)
Cspw  200    CONTINUE
Cspw         DDS=RADIUS*RADIUS
Cspw         DDT=DDOT(NC1,TVEC,1,TVEC,1)
Cspw         DDTD=DDOT(NC1,TVEC,1,D,1)
Cspw         IF(IUPHSS.EQ.1) THEN
Cspw            DDTD=DDTD/DDS
Cspw         ELSE
Cspw            ONE=1.0D+00
Cspw            COEF1=ONE-DDTD*DDTD/(DDS*DDT)
Cspw            COEF2=(ONE-COEF1)*DDS/DDTD
Cspw            COEF3=COEF1*DDTD/DDS
Cspw         END IF
Cspw         DO 220 I=1,NC1
Cspw            DO 210 J=1,I
Cspw               IF(IUPHSS.EQ.1) THEN
Cspw                  TEMP=TVEC(I)*D(J) + D(I)*TVEC(J) - D(I)*DDTD*D(J)
Cspw               ELSE
Cspw                  TEMP=(TVEC(I)*D(J) + D(I)*TVEC(J))*COEF1
Cspw     *                 -D(I)*D(J)*COEF3+TVEC(I)*TVEC(J)*COEF2
Cspw               END IF
Cspw               FC(I,J)=FC(I,J)+TEMP/DDS
Cspw               FC(J,I)=FC(I,J)
Cspw  210       CONTINUE
Cspw  220    CONTINUE
Cspw      END IF
CspwC
CspwC     ----- BFGS UPDATE.  NORMAL UPDATE FOR RUNTYP=OPTIMIZE -----
CspwC
Cspw      IF(IUPHSS.EQ.2) THEN
Cspw         IF (MASWRK  .AND.  NPRTGO.GT.0) WRITE(IW,9010)
Cspw         DO 300 I=1,NC1
Cspw             SVEC(I)=F(I)-OLDF(I)
Cspw  300    CONTINUE
Cspw         DDS=DDOT(NC1,SVEC,1,D,1)
CspwC
CspwC           IF DDS IS NEGATIVE, RETENSION OF POSITIVE
CspwC           DEFINITENESS IS NOT GUARANTEED.
CspwC
Cspw         IF(DDS.LT.ZERO .AND. MASWRK .AND. NPRTGO.GT.-1) WRITE(IW,9050)
Cspw         DDTD=DDOT(NC1,TVEC,1,D,1)
CspwC
Cspw         DO 320 I=1,NC1
Cspw            DO 310 J=1,I
Cspw               TEMP= (SVEC(I)*SVEC(J))/DDS - (TVEC(I)*TVEC(J))/DDTD
Cspw               FC(I,J)=FC(I,J)+TEMP
Cspw               FC(J,I)=FC(I,J)
Cspw  310       CONTINUE
Cspw  320    CONTINUE
Cspw      END IF
CspwC
CspwC     ----- SAVE UPDATED FORCE CONSTANT MATRIX -----
CspwC        (TRANSFORM FC TO CARTESIAN SPACE BEFORE PUNCHING)
CspwC
Cspw      IF(NZVAR.EQ.0 .AND. CLSTEP .AND. NEOTS) THEN
Cspw         WRITE(*,*) 'THE UPDATED CLFCM IN BKRUPD_NEO'
Cspw         CALL FCMOUT(FC,NCL3)
Cspw         CALL DAWRIT(IDAF,IODA,FC,NC3, 410,0)
Cspw         IF(NPUN.EQ.1  .OR.  NPUN.EQ.3) THEN
Cspw            IF(MASWRK) WRITE(IP,*) 'THE APPROXIMATE HESSIAN IS'
Cspw            CALL FCMPUN(FC,NC1)
Cspw         END IF
Cspw      ELSE
Cspw         CALL DAWRIT(IDAF,IODA,FC,NC3,43,0)
CspwC
Cspw         IF(NPUN.EQ.1  .OR.  NPUN.EQ.3) THEN
Cspw            CALL TFHB(FC,FCCART,NC1,NCOORD,NC1,NCOORD,49)
Cspw            IF(MASWRK) WRITE(IP,*) 'THE APPROXIMATE HESSIAN IS'
Cspw            CALL FCMPUN(FCCART,NCOORD)
Cspw         END IF
Cspw      END IF
Cspw      IF(OUT) WRITE (IW,9030)
Cspw      IF(OUT) CALL PRSQ(FC,NC1,NC1,NC1)
Cspw      RETURN
CspwC
Cspw 9000 FORMAT(10X,'HESSIAN UPDATED USING THE POWELL FORMULA')
Cspw 9050 FORMAT(10X,'WARNING! HEREDITARY POSITIVE DEFINITENESS ENDANGERED')
Cspw 9010 FORMAT(10X,'HESSIAN UPDATED USING THE BFGS FORMULA')
Cspw 9020 FORMAT(1X,'IN BKRUPD, FC =')
Cspw 9030 FORMAT(1X,'IN BKRUPD, THE UPDATED FC=')
Cspw      END
C
C*MODULE NEOOPT  *DECK UNFREZ
      SUBROUTINE UNFREZ(QMNFRZ)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL QMNFRZ
C
      PARAMETER (MXATM=2000)
C
      COMMON /FRZCRT/ IFCART(3*MXATM),NFCART
      COMMON /FROZNC/ IFREEZ(3*MXATM),NFRZ
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
C
C     --- QMNFRZ = .TRUE.
C         FREEZE TNE CARTESIAN COORDINATES OF THE QUANTUM NUCLEI ---
C
C     --- QMNFRZ = .FALSE.
C         FREEZE TNE CARTESIAN COORDINATES OF THE CLASSICAL NUCLEI ---
C
C     --- ZERO OUT THE IFCART AND IFREEZ VECTORS ---
C     (NEED CLEAR IFREEZ IN ORDER TO AVOID THE PARTIAL HESSIAN)
C
      NCOORD=3*NAT
C
      CALL VCLR(IFCART,1,NCOORD)
      CALL VCLR(IFREEZ,1,NCOORD)
C
      IFRZZ  = 0
      ICONUM = 0
CTDI  NTOTCN = NUNIQN + NAUXNB
      NTOTCN = NQMNUC + NTAUXB
      DO 150 IAT=1, NAT
         IF (QMNFRZ) THEN
            DO 110 I=1, NTOTCN
               IQM = IQMNUC(I)
               IQMC= IAUXNB(I)
               IF (IAT.EQ.IQM .OR. IAT.EQ.IQMC) THEN
                  ICONUM=IAT*3-2
                  DO 100 J=1,3
                    IFRZZ  = IFRZZ+1
                    IFCART(IFRZZ)=ICONUM
CTDI                IFREEZ(IFRZZ)=ICONUM
                    ICONUM = ICONUM +1
  100             CONTINUE
               END IF
  110       CONTINUE
         ELSE
            DO 120 I=1,NTOTCN
               IQM=IQMNUC(I)
               IQMC= IAUXNB(I)
               IF (IAT.EQ.IQM .OR. IAT.EQ.IQMC) GO TO 140
  120       CONTINUE
            ICONUM=IAT*3-2
            DO 130 J=1,3
               IFRZZ  = IFRZZ+1
               IFCART(IFRZZ)=ICONUM
CTDI           IFREEZ(IFRZZ)=ICONUM
               ICONUM = ICONUM +1
  130       CONTINUE
         END IF
  140 CONTINUE
  150 CONTINUE
      INDEX=IFRZZ/3
      WRITE(IW,9010)
      DO 160 I=1,INDEX
         JMIN = I*3-2
         JMAX = JMIN+2
         WRITE(IW,9000) (IFCART(J),J=JMIN,JMAX)
  160 CONTINUE
      WRITE(IW,*)
      RETURN
 9000 FORMAT(8X,I5,1X,I5,1X,I5)
 9010 FORMAT(/10X,17(1H-)/10X,'FROZEN CARTESIANS'/10X,17(1H-)/)
C              WRITE(*,*) 'IFCART(',IFRZZ,')',IFCART(IFRZZ)
      END
C
C*MODULE NEOOPT  *DECK NEOIZ
      SUBROUTINE NEOIZ(ZANI,ZANJ,IZI,IZJ,I,J)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
C
      IZI = INT(ZANI)
      IZJ = INT(ZANJ)
C
      IF(IZI.EQ.0) THEN
         IZI = 1
      ELSE
         IZI = IZI + IZCORE(I)
      END IF
C
      IF(IZJ.EQ.0) THEN
         IZJ = 1
      ELSE
         IZJ = IZJ + IZCORE(J)
      END IF
C
      RETURN
      END
