C C $Id: getdir.F,v 1.7 1998/07/16 16:39:54 jjv5 Exp arjan $ C C------------------------------------------------------------------------ SUBROUTINE GETDIR(METHOD,RESET,GRAD1,GNORM1,DIRECT1, . DNORM,DMAX,EHEAT1,IERROR) C C DETERMINES SEARCH DIRECTION FOR OPTIMIZATION. SUPPORTS STEEPEST C DESCENT, CONJUGATE-GRADIENT, AND BFGS. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" LOGICAL RESET CHARACTER*6 METHOD DIMENSION GRAD1(*),DIRECT1(*) C C C LOCAL: C DIMENSION HG(MAXPAR),XYZ0(3,MAXATM) LOGICAL FIRST, newsub DATA FIRST /.TRUE./ DATA DNLAST /100.0D0/ SAVE DNLAST,GNLAST,FIRST LOGICAL STEEP,CNJGRD,BFGS C STEEP = METHOD.EQ.'STEEP ' CNJGRD = METHOD.EQ.'CNJGRD' BFGS = METHOD.EQ.'BFGS ' C C FOR RESETS, STEP SIZE TO TAKE WHEN ESTIMATING HINV: C IF(FIRST)THEN DSTEP = 0.01D0 ELSE DSTEP = 0.002D0 ENDIF FIRST = .FALSE. C IF(STEEP.OR.(RESET.AND.CNJGRD))THEN C C SET DIRECTION TO THE NEGATIVE OF THE GRADIENT. C DO 100 I=1,NPAR DIRECT1(I) = -GRAD1(I) 100 CONTINUE ELSEIF(CNJGRD)THEN C C USE REGULAR CONJUGATE-GRADIENT FORMULA. C G1 = GNLAST*GNLAST G2 = GNORM1*GNORM1 DSCALE = G2/G1 DO 120 I=1,NPAR DIRECT1(I) = -GRAD1(I) + DSCALE*DIRECT1(I) 120 CONTINUE ELSEIF(RESET.AND.BFGS)THEN C C RESET INVERSE HESSIAN MATRIX FOR BFGS OPTIMIZATION AND GET C THE CORRESPONDING SEARCH DIRECTION -HINV*GRAD1. C DO 140 I=1,(NPAR*(NPAR+1))/2 HINV(I) = 0.0D0 140 CONTINUE C C ESTIMATE DIAGONAL ENTRIES BY CHANGING COORDINATES SLIGHTLY AND C USING THE CHANGE IN THE GRADIENT. STORE ORIGINAL INTERNALS C IN DCOORD AND CARTESIANS IN XYZ0. C IJ = 0 DO 160 I=1,NATOMS DO 150 J=1,3 XYZ0(J,I) = XYZ(J,I) IF(IOPT(J,I).EQ.1)THEN IJ = IJ + 1 DCOORD(IJ) = ZMAT(J,I) DIRECT1(IJ) = -SIGN(1.0D0,GRAD1(IJ)) ENDIF 150 CONTINUE 160 CONTINUE CALL GETCRD(DCOORD,DIRECT1,DSTEP,IERROR) IF(IERROR.NE.0) RETURN newsub = .true. C . do setup for pbc's if they have been requested if(PBC)then call setbox(ierror) if(ierror.ne.0) return endif if (newsub) then CALL GENSUB(IERROR) IF(IERROR.NE.0) RETURN endif CALL ENERGY(newsub,EH,IERROR) IF(IERROR.NE.0) RETURN call gcart(dgrad) CALL DOGRAD(DGRAD,DGNORM1,IERROR) IF(IERROR.NE.0) RETURN HMIN = 0.01D0 DGMIN = 0.01D0*DSTEP DXTEST = 0.06D0 II = 0 DO 180 I=1,NPAR II = II + I DX = -SIGN(DSTEP,GRAD1(I)) DG = DGRAD(I) - GRAD1(I) IF(ABS(DG).LT.DGMIN) DG = SIGN(DGMIN,DG) HII = DX/DG DGTEST = ABS(GRAD1(I)) IF(EH.LT.EHEAT1) DGTEST = ABS(DGRAD(I)) IF(HII.LT.0.0D0)THEN IF(DGTEST.LT.DGMIN)THEN HII = HMIN ELSE HII = DXTEST/DGTEST ENDIF ENDIF DGTEST = MAX(DGTEST,DGMIN) HMAX = 0.1D0/DGTEST HINV(II) = MIN(HII,HMAX) DIRECT1(I) = -HINV(II)*GRAD1(I) 180 CONTINUE C C REPLACE ENERGY AND GRADIENT IF THE CHANGE IN THE COORDINATES C LOWERED THE ENERGY. OTHERWISE, RESTORE THE INTERNAL AND C CARTESIAN COORDINATES TO THEIR ORIGINAL STATE. C IF(EH.LT.EHEAT1)THEN EHEAT1 = EH GNORM1 = 0.0D0 DO 210 I=1,NPAR GRAD1(I) = DGRAD(I) GNORM1 = GNORM1 + GRAD1(I)**2 210 CONTINUE GNORM1 = DSQRT(GNORM1) ELSE IJ = 0 DO 230 I=1,NATOMS DO 220 J=1,3 XYZ(J,I) = XYZ0(J,I) IF(IOPT(J,I).EQ.1)THEN IJ = IJ + 1 ZMAT(J,I) = DCOORD(IJ) ENDIF 220 CONTINUE 230 CONTINUE ENDIF ELSEIF(BFGS)THEN C C NORMAL BFGS UPDATE SCHEME FOR THE INVERSE HESSIAN MATRIX. C DG = 0.0D0 GHG = 0.0D0 IJ = 0 II = 0 DO 300 I=1,NPAR DG = DG + DCOORD(I)*DGRAD(I) HGI = 0.0D0 DO 240 J=1,I IJ = IJ + 1 HGI = HGI + HINV(IJ)*DGRAD(J) 240 CONTINUE IF(I.LT.NPAR)THEN II = II + I JI = II + I DO 260 J=I+1,NPAR HGI = HGI + HINV(JI)*DGRAD(J) JI = JI + J 260 CONTINUE ENDIF HG(I) = HGI GHG = GHG + DGRAD(I)*HGI 300 CONTINUE T = 1.0D0 + GHG/DG IJ = 0 DO 400 I=1,NPAR ITYP = IPAR(1,I) IATM = IPAR(2,I) DCRDI = DCOORD(I) HGI = HG(I) DO 380 J=1,I IJ = IJ + 1 JTYP = IPAR(1,J) JATM = IPAR(2,J) HGDIJ = HGI*DCOORD(J) DGHIJ = DCRDI*HG(J) DDIJ = DCRDI*DCOORD(J) HINV(IJ) = HINV(IJ) + (T*DDIJ - HGDIJ - DGHIJ)/DG 380 CONTINUE 400 CONTINUE C C UPDATE COMPLETE. NOW GET SEARCH DIRECTION = -HINV*GRAD1. C IJ = 0 II = 0 DO 500 I=1,NPAR DI = 0.0D0 DO 420 J=1,I IJ = IJ + 1 DI = DI - HINV(IJ)*GRAD1(J) 420 CONTINUE IF(I.LT.NPAR)THEN II = II + I JI = II + I DO 440 J=I+1,NPAR DI = DI - HINV(JI)*GRAD1(J) JI = JI + J 440 CONTINUE ENDIF DIRECT1(I) = DI 500 CONTINUE ENDIF C C GET NORM OF SEARCH DIRECTION VECTOR AND MAXIMUM ENTRY. C DNORM = 0.0D0 DMAX = 0.0D0 DO 600 I=1,NPAR DNORM = DNORM + DIRECT1(I)**2 DMAX = MAX(DMAX,ABS(DIRECT1(I))) 600 CONTINUE DNORM = DSQRT(DNORM) C C IF WE'RE DOING A BFGS, THEN SCALE DOWN DIRECT IF ITS NORM HAS C INCREASED SIGNIFICANTLY FROM THE PREVIOUS CYCLE. C IF(BFGS)THEN IF(DNORM.GT.1.5D0*DNLAST)THEN DFACTR = 1.5D0*DNLAST/DNORM DSCALE = DSCALE*DFACTR DO 840 I=1,NPAR DIRECT1(I) = DIRECT1(I)*DFACTR 840 CONTINUE DNORM = DNORM*DFACTR DMAX = DMAX*DFACTR ENDIF ENDIF DNLAST = DNORM GNLAST = GNORM1 RETURN END CC-------------------------------------------------------------------CC