mopac606_nbo/src/linmin.f

   1       SUBROUTINE LINMIN(XPARAM,ALPHA,PVECT,NVAR,FUNCT,OKF,IC, DOTT)
   2       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
   3       INCLUDE 'SIZES'
   4       DIMENSION XPARAM(NVAR),PVECT(NVAR)
   5       COMMON /GRAVEC/ COSINE
   6       COMMON /NUMCAL/ NUMCAL
   7 C*********************************************************************
   8 C
   9 C  LINMIN DOES A LINE MINIMISATION.
  10 C
  11 C  ON INPUT:  XPARAM = STARTING COORDINATE OF SEARCH.
  12 C             ALPHA  = STEP SIZE FOR INITIATING SEARCH.
  13 C             PVECT  = DIRECTION OF SEARCH.
  14 C             NVAR   = NUMBER OF VARIABLES IN XPARAM.
  15 C             FUNCT  = INITIAL VALUE OF THE FUNCTION TO BE MINIMIZED.
  16 C             ISOK   = NOT IMPORTANT.
  17 C             COSINE = COSINE OF ANGLE OF CURRENT AND PREVIOUS GRADIENT.
  18 C
  19 C  ON OUTPUT: XPARAM = COORDINATE OF MINIMUM OF FUNCTI0N.
  20 C             ALPHA  = NEW STEP SIZE, USED IN NEXT CALL OF LINMIN.
  21 C             FUNCT  = FINAL, MINIMUM VALUE OF THE FUNCTION.
  22 C             OKF    = TRUE IF LINMIN IMPROVED FUNCT, FALSE OTHERWISE.
  23 C
  24 C**********************************************************************
  25       COMMON /KEYWRD/ KEYWRD
  26 C
  27 C  THE FOLLOWING COMMON IS USED TO FIND OUT IF A NON-VARIATIONALLY
  28 C  OPTIMIZED WAVE-FUNCTION IS BEING USED.
  29 C
  30       COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
  31      1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
  32      2                NCLOSE,NOPEN,NDUMY,FRACT
  33       CHARACTER KEYWRD*241
  34       DIMENSION PHI(3), VT(4)
  35       DIMENSION XSTOR(MAXPAR), XPAREF(MAXPAR)
  36       INTEGER LEFT,RIGHT,CENTER
  37       LOGICAL PRINT,OKF, HALFE, DIIS
  38       SAVE ICALCN, PRINT, HALFE, XMAXM, I
  39       SAVE MAXLIN,  YMAXST
  40       DATA ICALCN /0/
  41       IF (ICALCN.NE.NUMCAL) THEN
  42          HALFE =(INDEX(KEYWRD,'C.I.') .NE. 0 .OR. NCLOSE.NE.NOPEN)
  43          IF(INDEX(KEYWRD,'GNORM') .NE. 0)
  44      1DROP=DROP*MIN(READA(KEYWRD,INDEX(KEYWRD,'GNORM')),1.D0)
  45          XMAXM  = 0.4D0
  46          DELTA2 = 0.001D0
  47          IF(INDEX(KEYWRD,'NOTH') .EQ. 0) THEN
  48             DELTA1 = 0.5D0
  49          ELSE
  50             DELTA1 = 0.1D0
  51          ENDIF
  52          ALPHA  = 1.D0
  53          MAXLIN = 15
  54          IF(NVAR.EQ.1)THEN
  55             PVECT(1)=0.01D0
  56             DROP=0.01D0
  57             ALPHA=1.D0
  58             DELTA1 = 0.00005D0
  59             DELTA2 = 0.00001D0
  60             IF(INDEX(KEYWRD,'PREC') .NE. 0) DELTA1=0.0000005
  61             MAXLIN=30
  62          ENDIF
  63          COSINE=99.99D0
  64 C
  65          YMAXST  = 0.4D0
  66          PRINT=(INDEX(KEYWRD,'LINMIN') .NE. 0)
  67          ICALCN=NUMCAL
  68       ENDIF
  69       DO 10 I=1,NVAR
  70          XPAREF(I)=XPARAM(I)
  71    10 CONTINUE
  72       XMAXM=0.D0
  73       DO 20 I=1,NVAR
  74          PABS=ABS(PVECT(I))
  75    20 XMAXM=MAX(XMAXM,PABS)
  76       XMAXM=YMAXST/XMAXM
  77       IF(NVAR.EQ.1)
  78      1CALL COMPFG(XPARAM, .TRUE., FUNCT,.TRUE.,GRAD,.FALSE.)
  79       FIN=FUNCT
  80       SSQLST=FUNCT
  81       DIIS=IC.EQ.1.AND.NVAR.GT.1
  82       PHI(1)=FUNCT
  83       ALPHA=1.D0
  84       VT(1)=0.0D00
  85       VT(2)=ALPHA
  86       IF (VT(2).GT.XMAXM) VT(2)=XMAXM
  87       FMAX=FUNCT
  88       FMIN=FUNCT
  89       ALPHA=VT(2)
  90       DO 30 I=1,NVAR
  91    30 XPARAM(I)=XPAREF(I)+ALPHA*PVECT(I)
  92       CALL COMPFG(XPARAM, .TRUE., PHI(2),.TRUE.,GRAD,.FALSE.)
  93       IF(PHI(2).GT.FMAX) FMAX=PHI(2)
  94       IF(PHI(2).LT.FMIN) FMIN=PHI(2)
  95       CALL EXCHNG (PHI(2),SQSTOR,ENERGY,ESTOR,XPARAM,XSTOR,
  96      1ALPHA,ALFS,NVAR)
  97       IF(DIIS)GOTO 190
  98       IF(NVAR.GT.1)THEN
  99 C
 100 C   CALCULATE A NEW ALPHA BASED ON THIEL'S FORMULA
 101 C
 102          ALPHA=-ALPHA**2*DOTT/(2.D0*(PHI(2)-SSQLST-ALPHA*DOTT))
 103          IF(ALPHA.GT.2.D0)ALPHA=2.D0
 104       ELSE
 105          IF(PHI(2).LT.PHI(1))THEN
 106             ALPHA=2*ALPHA
 107          ELSE
 108             ALPHA=-ALPHA
 109          ENDIF
 110       ENDIF
 111 C#      IF(PRINT)WRITE(6,'(3(A,F12.6))')' ESTIMATED DROP:',DOTT*0.5D0,
 112 C#     1'  ACTUAL: ',PHI(2)-SSQLST, '  PREDICTED ALPHA',ALPHA
 113       OKF=OKF.OR.PHI(2).LT.SSQLST
 114       IF(DELTA1.GT.0.3D0)THEN
 115 C
 116 C  THIEL'S TESTS # 18 AND 19
 117 C
 118          IF(OKF.AND.ALPHA.LT.2.D0)GOTO 190
 119       ENDIF
 120       VT(3)=ALPHA
 121       IF (VT(3).LE.1.D0) THEN
 122          LEFT=3
 123          CENTER=1
 124          RIGHT=2
 125       ELSE
 126          LEFT=1
 127          CENTER=2
 128          RIGHT=3
 129       ENDIF
 130       DO 40 I=1,NVAR
 131    40 XPARAM(I)=XPAREF(I)+ALPHA*PVECT(I)
 132       CALL COMPFG (XPARAM, .TRUE., FUNCT,.TRUE.,GRAD,.FALSE.)
 133       IF(FUNCT.GT.FMAX) FMAX=FUNCT
 134       IF(FUNCT.LT.FMIN) FMIN=FUNCT
 135       IF (FUNCT.LT.SQSTOR) CALL EXCHNG (FUNCT,SQSTOR,ENERGY,
 136      1ESTOR,XPARAM,XSTOR,ALPHA,ALFS,NVAR)
 137       OKF=(OKF.OR.FUNCT.LT.FIN)
 138       PHI(3)=FUNCT
 139       IF (PRINT)WRITE (6,50) VT(1),PHI(1),PHI(1)-FIN,
 140      1                        VT(2),PHI(2),PHI(2)-FIN,
 141      2                        VT(3),PHI(3),PHI(3)-FIN
 142    50 FORMAT ( ' ---QLINMN ',/5X, 'LEFT   ...',F17.8,2F17.11/5X,
 143      1 'CENTER ...',F17.8,2F17.11,/5X, 'RIGHT  ...',F17.8,2F17.11,/)
 144       DO 180 ICTR=3,MAXLIN
 145          ALPHA=VT(2)-VT(3)
 146          BETA=VT(3)-VT(1)
 147          GAMMA=VT(1)-VT(2)
 148          IF(ABS(ALPHA*BETA*GAMMA) .GT. 1.D-4)THEN
 149             ALPHA=-(PHI(1)*ALPHA+PHI(2)*BETA+PHI(3)*GAMMA)/(ALPHA*BETA*G
 150      1AMM   A)
 151          ELSE
 152 C
 153 C   FINISH BECAUSE TWO POINTS CALCULATED ARE VERY CLOSE TOGETHER
 154 C
 155             GOTO 190
 156          ENDIF
 157          BETA=((PHI(1)-PHI(2))/GAMMA)-ALPHA*(VT(1)+VT(2))
 158          IF (ALPHA) 60,60,90
 159    60    IF (PHI(RIGHT).GT.PHI(LEFT)) GO TO 70
 160          ALPHA=3.0D00*VT(RIGHT)-2.0D00*VT(CENTER)
 161          GO TO 80
 162    70    ALPHA=3.0D00*VT(LEFT)-2.0D00*VT(CENTER)
 163    80    S=ALPHA-ALPOLD
 164          IF (ABS(S).GT.XMAXM) S=SIGN(XMAXM,S)*(1+0.01*(XMAXM/S))
 165          ALPHA=S+ALPOLD
 166          GO TO 100
 167    90    ALPHA=-BETA/(2.0D00*ALPHA)
 168          S=ALPHA-ALPOLD
 169          XXM=2.0D00*XMAXM
 170          IF (ABS(S).GT.XXM) S=SIGN(XXM,S)*(1+0.01*(XXM/S))
 171          ALPHA=S+ALPOLD
 172   100    CONTINUE
 173 C
 174 C   FINISH IF CALCULATED POINT IS NEAR TO POINT ALREADY CALCULATED
 175 C
 176          DO 110 I=1,3
 177   110    IF (ABS(ALPHA-VT(I)).LT.DELTA1*(1.D0+VT(I)).AND.OKF) GOTO 190
 178          DO 120 I=1,NVAR
 179   120    XPARAM(I)=XPAREF(I)+ALPHA*PVECT(I)
 180          FUNOLD=FUNCT
 181          CALL COMPFG (XPARAM, .TRUE., FUNCT,.TRUE.,GRAD,.FALSE.)
 182          IF(FUNCT.GT.FMAX) FMAX=FUNCT
 183          IF(FUNCT.LT.FMIN) FMIN=FUNCT
 184          IF (FUNCT.LT.SQSTOR) CALL EXCHNG (FUNCT,SQSTOR,ENERGY,ESTOR,
 185      1   XPARAM,XSTOR,ALPHA,ALFS,NVAR)
 186          OKF=OKF .OR. (FUNCT.LT.FIN)
 187          IF (PRINT) WRITE(6,130) VT(LEFT),PHI(LEFT), PHI(LEFT)-FIN,
 188      1                           VT(CENTER),PHI(CENTER),PHI(CENTER)-FIN,
 189      2                           VT(RIGHT),PHI(RIGHT),PHI(RIGHT)-FIN,
 190      3                           ALPHA,FUNCT,FUNCT-FIN
 191   130    FORMAT (5X,'LEFT    ...',F17.8,2F17.11,/5X,'CENTER  ...',
 192      1F17.8,2F17.11,/5X,'RIGHT   ...',F17.8,2F17.11,/5X,
 193      2 'NEW     ...',F17.8,2F17.11,/)
 194 C
 195 C TEST TO EXIT FROM LINMIN IF NOT DROPPING IN VALUE OF FUNCTION FAST.
 196 C
 197          IF(ABS(FUNOLD-FUNCT) .LT. DELTA2 .AND. OKF) GOTO 190
 198          ALPOLD=ALPHA
 199          IF ((ALPHA.GT.VT(RIGHT)).OR.(ALPHA.GT.VT(CENTER)
 200      1        .AND.FUNCT.LT.PHI(CENTER)).OR.(ALPHA.GT.VT(LEFT)
 201      2        .AND.ALPHA.LT.VT(CENTER).AND.FUNCT.GT.PHI(CENTER)))
 202      3         GOTO 140
 203          VT(RIGHT)=ALPHA
 204          PHI(RIGHT)=FUNCT
 205          GO TO 150
 206   140    VT(LEFT)=ALPHA
 207          PHI(LEFT)=FUNCT
 208   150    IF (VT(CENTER).LT.VT(RIGHT)) GO TO 160
 209          I=CENTER
 210          CENTER=RIGHT
 211          RIGHT=I
 212   160    IF (VT(LEFT).LT.VT(CENTER)) GO TO 170
 213          I=LEFT
 214          LEFT=CENTER
 215          CENTER=I
 216   170    IF (VT(CENTER).LT.VT(RIGHT)) GO TO 180
 217          I=CENTER
 218          CENTER=RIGHT
 219          RIGHT=I
 220   180 CONTINUE
 221   190 CONTINUE
 222 C
 223 C  IC=1 IF THE LAST POINT CALCULATED WAS THE BEST POINT, IC=2 OTHERWISE
 224 C
 225       IC=2
 226       IF(ABS(ESTOR-ENERGY).LT.1.D-12)IC=1
 227       CALL EXCHNG (SQSTOR,FUNCT,ESTOR,ENERGY,XSTOR,XPARAM,
 228      1             ALFS,ALPHA,NVAR)
 229       OKF = (FUNCT.LT.SSQLST.OR.DIIS)
 230       IF (FUNCT.GE.SSQLST) RETURN
 231       IF (ALPHA) 200,220,220
 232   200 ALPHA=-ALPHA
 233       DO 210 I=1,NVAR
 234   210 PVECT(I)=-PVECT(I)
 235   220 CONTINUE
 236       RETURN
 237 C
 238 C
 239       END