C C $Id: geoopt.F,v 1.9 1999/03/02 20:02:25 arjan Exp arjan $ C C------------------------------------------------------------------------ SUBROUTINE GEOOPT(EHEAT1,IERROR) C C GEOMETRY OPTIMIZATION ROUTINE. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) #include "divcon.dim" #include "divcon.h" C C LOCAL ARRAYS AND VARIABLES: C CHARACTER METHOD*6 LOGICAL STEEP,CNJGRD,BFGS,RESET,DIIS,DIISOK,FORCIT, . THERE,newsub DIMENSION DIRECT1(MAXPAR),OLDGRD(MAXPAR) C CALL ETIMER(TPREV) IERROR = 0 C STEEP = INDEX(KEYWRD,'OPT=STEEP').NE.0 CNJGRD = INDEX(KEYWRD,'OPT=CONJGRAD').NE.0 BFGS = INDEX(KEYWRD,'OPT=BFGS').NE.0 IF(STEEP)THEN METHOD = 'STEEP ' ELSEIF(CNJGRD)THEN METHOD = 'CNJGRD' ELSEIF(BFGS)THEN METHOD = 'BFGS ' ELSE METHOD = 'CNJGRD' ENDIF C C CRITERIA FOR OPTIMIZATION: ETEST --> CHANGE IN ENERGY C XTEST --> CHANGE IN COORDINATES C GTEST --> GRADIENT ENTRY TEST C IETEST = INDEX(KEYWRD,'ETEST=') IF(IETEST.NE.0)THEN CALL RDNUM(KEYWRD,IETEST,ETEST,IERROR) ELSE ETEST = 0.002D0 ENDIF C IXTEST = INDEX(KEYWRD,'XTEST=') IF(IXTEST.NE.0)THEN CALL RDNUM(KEYWRD,IXTEST,XTEST,IERROR) ELSE XTEST = 0.001D0 ENDIF C IGTEST = INDEX(KEYWRD,'GTEST=') IF(IGTEST.NE.0)THEN CALL RDNUM(KEYWRD,IGTEST,GTEST,IERROR) ELSE GTEST = 0.5D0 ENDIF C C SET UP CORRESPONDENCE BETWEEN THE OPTIMIZABLE PARAMETER NUMBER C AND THE Z-MATRIX OR XYZ ARRAY. C CALL SETOPT IF(NPAR.EQ.0) RETURN C C GRADIENT NORM CRITERION: C SQNPAR = DSQRT(DFLOAT(NPAR)) GNTEST = MAX(GTEST*SQNPAR*0.25D0,GTEST) C-RDC WRITE(IOUT,'(/" ETEST =",F10.6, C-RDC . /" XTEST =",F10.6, C-RDC . /" GNTEST =",F10.6, C-RDC . /" GTEST =",F10.6)') ETEST,XTEST,GNTEST,GTEST C C BINARY SEARCH CRITERA: C BXTEST = 0.001D0 BGTEST = 0.5D0 C C C DIRECT INVERSION OF ITERATIVE SUBSPACE FLAGS: C C DIISOK - DETERMINES WHETHER IT'S OKAY TO USE DIIS IN OPTIMIZATION. C SHOULD BE .FALSE. IF THE USER DOESN'T WANT DIIS, OR IF C WE'RE NOT DOING A BFGS, OR IF THERE ARE VERY FEW PARAMETERS. C C DIIS - TELLS LINMIN AND DIISUB WHETHER OR NOT THE DIIS STEP C WILL BE USED. C DIISOK = INDEX(KEYWRD,'NODIIS').EQ.0.AND.BFGS.AND.NPAR.GT.5 DIIS = DIISOK C C GEOMETRY OPTIMIZATION WILL HALT IF THE ENERGY INCREASES ON THREE C SUCCESSIVE CYCLES UNLESS THE USER OVERRIDES IT WITH THE KEYWORD C 'FORCE-IT'. C NRAISE = 0 FORCIT = INDEX(KEYWRD,'FORCE-IT').NE.0 C C DETERMINE MAXIMUM NUMBER OF CYCLES IN OPTIMIZATION. C MAXCYC = 10*NPAR IMAX = INDEX(KEYWRD,'MAXOPT=') IF(IMAX.NE.0)THEN CALL RDNUM(KEYWRD,IMAX,OPTMAX,IERROR) MAXCYC = OPTMAX ENDIF C C FREQUENCY FOR WRITING RESTART FILE: C NDUMP = 100 IDUMP = INDEX(KEYWRD,'DUMP=') IF(IDUMP.NE.0)THEN CALL RDNUM(KEYWRD,IDUMP,DUMP,IERROR) NDUMP = DUMP ENDIF C C GET INITIAL ENERGY AND GRADIENT. C NSCF = 0 newsub = .true. call setbox(ierror) if (ierror.ne.0) return CALL GENSUB(IERROR) IF(IERROR.NE.0) RETURN CALL ENERGY(newsub,EHEAT1,IERROR) IF(IERROR.NE.0) RETURN call gcart(grad) CALL DOGRAD(GRAD,GNORM,IERROR) IF(IERROR.NE.0) RETURN ELAST = 0.0D0 DELTAE = 0.0D0 DEPRED = 100.0D0 DELTAX = 0.0D0 DNORM = 1.0D0 RESET = .TRUE. C C BEGIN OPTIMIZATION. C DO 1000 ICYCLE=0,MAXCYC IF(ICYCLE.NE.0.AND.MOD(ICYCLE,NDUMP).EQ.0) then CALL WRTRST c call wrtdmx(ierror) endif C sld C C ONLY FOR LARGE SYSTEMS: C INQUIRE(FILE='stop.dat',EXIST=THERE) IF(THERE)THEN C-RDC WRITE(IOUT,'(" FILE stop.dat FOUND, RETURNING")') C-RDC WRITE(0,'(" FILE stop.dat FOUND, RETURNING")') RETURN ENDIF C C end sld C C C GET MAXIMUM COMPONENT OF GRADIENT. C GRDMAX = 0.0D0 DO 320 I=1,NPAR GRDMAX = MAX(GRDMAX,ABS(GRAD(I))) 320 CONTINUE C CALL ETIMER(TNOW) TCYCLE = TNOW - TPREV TPREV = TNOW C C OUTPUT LATEST INFORMATION. C C-RDC WRITE(IOUT,'(/" CYCLE =",I6,9X,"TIME =",F11.3, C-RDC . " ENERGY =",F15.6/" DELTAE =",F12.6, C-RDC . " GNORM =",F11.3," GRDMAX =",F10.4, C-RDC . " DELTAX =",F9.6)') ICYCLE,TCYCLE,EHEAT1, C-RDC . DELTAE,GNORM,GRDMAX,DELTAX IF(SCREEN)THEN C-RDC WRITE(ISCR,'(/" CYCLE =",I6,9X,"TIME =",F11.3, C-RDC . " ENERGY =",F15.6/" DELTAE =",F12.6, C-RDC . " GNORM =",F11.3," GRDMAX =",F10.4, C-RDC . " DELTAX =",F9.6)') ICYCLE,TCYCLE,EHEAT1, C-RDC . DELTAE,GNORM,GRDMAX,DELTAX ENDIF C IF(ICYCLE.EQ.0)THEN DELTAE = EHEAT1 DELTAX = 100.0D0 ENDIF C C SEE IF VARIOUS OPTIMIZATION TESTS HAVE BEEN PASSED. C NPASS = 0 IF(ABS(DELTAE).LT.ETEST.AND.DELTAE.LT.0.0D0)THEN NPASS = NPASS + 1 C-RDC WRITE(IOUT,'(" ENERGY TEST PASSED")') C-RDC IF(SCREEN) WRITE(ISCR,'(" ENERGY TEST PASSED")') ENDIF IF(DELTAX.LT.XTEST)THEN NPASS = NPASS + 1 C-RDC WRITE(IOUT,'(" COORDINATE TEST PASSED")') C-RDC IF(SCREEN) WRITE(ISCR,'(" COORDINATE TEST PASSED")') ENDIF IF(GNORM.LT.GNTEST)THEN NPASS = NPASS + 1 C-RDC WRITE(IOUT,'(" GRADIENT NORM TEST PASSED")') C-RDC IF(SCREEN) WRITE(ISCR,'(" GRADIENT NORM TEST PASSED")') ENDIF IF(GRDMAX.LT.GTEST)THEN NPASS = NPASS + 1 C-RDC WRITE(IOUT,'(" GRADIENT COMPONENT TEST PASSED")') C-RDC IF(SCREEN) WRITE(ISCR,'(" GRADIENT COMPONENT TEST PASSED")') ENDIF IF(NPASS.EQ.4)THEN C-RDC WRITE(IOUT,'(/" -- GEOMETRY OPTIMIZED --")') C-RDC IF(SCREEN) WRITE(ISCR,'(/" -- GEOMETRY OPTIMIZED --")') IF(INDEX(KEYWRD,'DUMP=').NE.0) then CALL WRTRST c call wrtdmx(ierror) endif RETURN ENDIF C C NO NEED TO CONTINUE OF THE MAXIMUM NUMBER OF OPTIMIZATION CYCLES C HAS ELAPSED. C IF(ICYCLE.EQ.MAXCYC) GO TO 1000 C C IF ENERGY INCREASES ON SUCCESSIVE CYCLES, KEEP TRACK OF IT. C MAY WANT TO RESET SEARCH DIRECTION OR EXIT. C IF(DELTAE.GT.0.0D0)THEN NRAISE = NRAISE + 1 ELSE NRAISE = 0 ENDIF C IF(NRAISE.GE.3.AND.DELTAX.LT.0.001D0.AND..NOT.FORCIT)THEN C-RDC WRITE(IOUT,'(/" THE ENERGY HAS INCREASED ON THREE OR MORE", C-RDC . " CONSECUTIVE CYCLES, BUT THE", C-RDC . /" GEOMETRY IS CHANGING BY VERY LITTLE. THE", C-RDC . " GRADIENT IS PROBABLY NOT RELIABLE", C-RDC . /" AT THIS POINT, SO THE OPTIMIZATION IS BEING", C-RDC . " HALTED. OVERRIDE WITH THE", C-RDC . /" KEYWORD ''FORCE-IT''")') RETURN ELSEIF(NRAISE.GE.2.AND.DELTAX.LT.0.001D0.AND.FORCIT)THEN C-RDC WRITE(IOUT,'(/" ENERGY IS INCREASING BUT ''FORCE-IT'' WAS", C-RDC . " SPECIFIED . . .", C-RDC . /" IF THE ENERGY CONTINUES TO INCREASE, THE", C-RDC . " USER CAN TERMINATE THE", C-RDC . /" OPTIMIZATION GRACEFULLY BY CREATING THE", C-RDC . " FILE stop.opt")') IF(SCREEN)THEN C-RDC WRITE(ISCR,'(/" ENERGY IS INCREASING BUT ''FORCE-IT'' WAS", C-RDC . " SPECIFIED . . .", C-RDC . /" IF THE ENERGY CONTINUES TO INCREASE, THE", C-RDC . " USER CAN TERMINATE THE", C-RDC . /" OPTIMIZATION GRACEFULLY BY CREATING THE", C-RDC . " FILE stop.opt")') ENDIF INQUIRE(FILE='stop.opt',EXIST=THERE) IF(THERE)THEN C-RDC WRITE(IOUT,'(" TERMINATING OPTIMIZATION BECAUSE THE FILE", C-RDC . " stop.opt WAS FOUND")') RETURN ENDIF ENDIF C C IF WE'RE DOING A BFGS, THEN GET THE CHANGE IN THE COORDINATES C AND GRADIENT TO UPDATE THE INVERSE HESSIAN IN THE CALL TO C GETDIR. C IF(BFGS.AND.ICYCLE.GT.0)THEN DO 350 I=1,NPAR DCOORD(I) = ALPHA*DIRECT1(I) DGRAD(I) = GRAD(I) - OLDGRD(I) 350 CONTINUE ENDIF C C GET LATEST SEARCH DIRECTION. RESET IF ENERGY HAS INCREASED C ON TWO SUCCESSIVE CYCLES. C IF(.NOT.RESET)THEN RESET = NRAISE.GE.2 IF(RESET)THEN C-RDC WRITE(IOUT,'(/" RESETTING SEARCH DIRECTION")') C-RDC IF(SCREEN) WRITE(ISCR,'(/" RESETTING SEARCH DIRECTION")') NRAISE = 0 ENDIF ENDIF DNPREV = DNORM ELAST = EHEAT1 CALL GETDIR(METHOD,RESET,GRAD,GNORM,DIRECT1, . DNORM,DMAX,EHEAT1,IERROR) IF(IERROR.NE.0) RETURN C C SEE IF SEARCH DIRECTION HAS BECOME ORTHOGONAL TO GRADIENT. C IF SO, THEN RESET. WILL NEED GDOTD BELOW TO GET PREDICTED C ENERGY DROP, EVEN IF WE'RE ON A RESET ALREADY. C GDOTD = 0.0D0 DO 400 I=1,NPAR GDOTD = GDOTD + GRAD(I)*DIRECT1(I) 400 CONTINUE IF(.NOT.RESET.AND..NOT.STEEP)THEN COSINE = -GDOTD/(GNORM*DNORM) IF((COSINE.LT.0.15D0.AND.DEPRED.GT.1.0D0).OR. . COSINE.LT.0.05D0)THEN RESET = .TRUE. C-RDC WRITE(IOUT,'(/" RESETTING SEARCH DIRECTION")') C-RDC IF(SCREEN) WRITE(ISCR,'(/" RESETTING SEARCH DIRECTION")') CALL GETDIR(METHOD,RESET,GRAD,GNORM,DIRECT1, . DNORM,DMAX,EHEAT1,IERROR) IF(IERROR.NE.0) RETURN ENDIF ENDIF C C STORE LATEST GRADIENT FOR USE IN NEXT CYCLE OF BFGS. C IF(BFGS)THEN DO 450 I=1,NPAR OLDGRD(I) = GRAD(I) 450 CONTINUE ENDIF C C GET APPROPRIATE STEP SIZE TO TAKE IN LINE SEARCH. C DIIS = DIISOK IF(ICYCLE.GT.0)THEN IF(ALPHA.GT.1.0D-5)THEN C C COUNTERACT FLUCTUATIONS IN THE SIZE OF DIRECT. C STEP = MAX(ALPHA*DNPREV/DNORM,0.001D0/DMAX) ELSE C C WE PROBABLY TOOK TOO LARGE OF A STEP LAST CYCLE. C STEP = MAX(0.1D0*STEP*DNPREV/DNORM,0.001D0) ENDIF ELSE STEP = 0.01D0/DMAX ENDIF STEP = MIN(STEP,0.02D0/DMAX) DEPRED = ABS(GDOTD*STEP) C C DO LINE MINIMIZATION, DIIS, OR BINARY SEARCH. C IF(GRDMAX.GE.BGTEST.OR.DELTAX.GE.BXTEST)THEN C C IF DIIS IS TURNED ON, THEN THE NORMAL PROCEDURE IN LINMIN WOULD C BE TO TAKE ONE STEP AND RETURN. THEN, AFTER COMPUTING GRADIENT, C DIISUB WOULD BE CALLED TO GET INTERPOLATED COORDINATES AND C GRADIENT. IF, HOWEVER, THE FIRST STEP IN LINMIN CAUSES THE C ENERGY TO GO UP, THEN LINMIN WILL BE ALLOWED TO PROCEED WITH C A FULL LINE MINIMIZATION. LINMIN WILL ALSO SET THE DIIS FLAG C TO .FALSE. IN THIS CASE, SO THAT THE UPCOMING CALL TO DIISUB C WILL RESULT ONLY IN THE COORDINATES, GRADIENT, AND ENERGY BEING C STORED. C CALL LINMIN(DIIS,DIRECT1,GRAD,STEP,EHEAT1,ALPHA,IERROR) IF(IERROR.NE.0)THEN C-RDC WRITE(IOUT,'(/" *** LINMIN ERROR ***")') C-RDC IF(SCREEN) WRITE(ISCR,'(/" *** LINMIN ERROR ***")') RETURN ENDIF ELSE CALL BSRCH(DIRECT1,STEP,EHEAT1,ALPHA,IERROR) IF(IERROR.NE.0)THEN C-RDC WRITE(IOUT,'(/" *** BINARY SEARCH ERROR ***")') C-RDC IF(SCREEN) WRITE(ISCR,'(/" *** BINARY SEARCH ERROR ***")') RETURN ENDIF C C SINCE WE'RE DOING A BINARY SEARCH, WE ARE PROBABLY SO C CLOSE TO THE MINIMUM THAT USING DIIS IS NOT ADVISABLE. C DIIS = .FALSE. ENDIF C C MAKE A CALL TO DIISUB IF WE ARE NOT PROHIBITED FROM DOING SO. C IF(DIISOK)THEN C C NEED LATEST GRADIENT FOR DIISUB. C call gcart(grad) CALL DOGRAD(GRAD,GNORM,IERROR) IF(IERROR.NE.0) RETURN C C IF DIIS IS .FALSE. THEN DIISUB WILL MERELY STORE COORDINATES, C GRADIENT, AND ENERGY. THIS WILL BE THE CASE WHENEVER A FULL C LINE MINIMIZATION HAS ALREADY BEEN DONE, OR WHEN WE ARE SO C CLOSE TO THE MINIMUM THAT WE DID A BINARY SEARCH. ALSO, DIIS C CANNOT BE USED ON THE FIRST CALL, AND DIISUB WILL REALIZE C THIS, EVEN IF DIIS IS .TRUE. C CALL DIISUB(DIIS,GRAD,EHEAT1,DIRECT1,DNORM,GNORM,ALPHA) ENDIF C C GET THE GRADIENT, UNLESS IT WAS ALREADY COMPUTED IN DIISUB. C IF(.NOT.DIISOK)THEN call gcart(grad) CALL DOGRAD(GRAD,GNORM,IERROR) IF(IERROR.NE.0) RETURN ENDIF C C DELTAX IS THE RMS COORDINATE CHANGE. C DELTAX = DNORM*ALPHA/SQNPAR DELTAE = EHEAT1 - ELAST RESET = .FALSE. 1000 CONTINUE C-RDC WRITE(IOUT,'(" MAXIMUM NUMBER OF OPTIMIZATION CYCLES", C-RDC . " COMPLETED"/" -- INCREASE BY USING THE", C-RDC . " ''MAXOPT='' KEYWORD")') RETURN END