PROGRAM stringMB C CCC for mb potential C Vorlage: Stacho / Ban ... DDRP C Computers Chem 17, No1 (1993) 21-25 C ************* C W.Quapp 12.12.2003 for IRC, RGF, GE and all that C CCCCCC REAL*8 D,Dh, Y(250,2),Z(250,2),ETA,R,EPS Integer LY, LZ, L, N , M, J, IGES, LMAX, ITERMAX c Mathematica - file of the result: OPEN(7,FILE='mbChain.dat') c if using a given initial chain c OPEN(8,FILE='mbStart.dat') rewind 7 rewind 8 C Constants c Dimension N=2 C Damping ETA=0.125d0 C Convergence criterion eps=0.04d0 C Chain length L= 33 LZ=33 LY=33 LMAX=133 ITERMAX=100 R=2.75d0 M=1 c SPs of mb: c (-0.822, 0.624) & ( 0.212, 0.293) c Minima of mb: c (-0.558, 1.442) & (-0.050, 0.467) & ( 0.624, 0.028) write(7,*) '{ ' PRINT*,' Program for chain from Min1 to Min3 on MB' PRINT*,' ' c Put a Chain of Points between the Minima Do 444 J=1,L Y(J,1)= (-0.558D0*(J-1) + (L-J)*0.624D0)/(L-1) Y(J,2)= ( 1.442D0*(J-1) + (L-J)*0.028D0)/(L-1) c or use: c Read(8,*) Y(J,1),Y(J,2) 444 continue IGES=1 443 continue c CALL NEWCURVE(N,L,Y,Z,ETA,M,R) LY=L LZ=L call HAUSD(N,Y,LY,Z,LZ,D) WRITE(6, *) 'Iter=',IGES,' Distance of Chains Y , Z = ', D c c Dh is an alternate convergence criterion: call HeighestPoint(N,Z,LZ,ETA,Dh) WRITE(6, *) 'Iter=',IGES,' Proj at SP = ', Dh c c for homogenization only, not used here: c CALL HOMOGEN(N,Y,LY,Z,LZ,EPS,LMAX) c WRITE(6, *) ' LYneu = ', LY c Do 440 J=1,L Do 440 II=1,N Y(J,II)=Z(J,II) 440 continue L=LY If( (IGES .gt. 1) .AND. (D .lt. EPS ) ) GOTO 9 cc If( (IGES .gt. 1) .AND. (Dh .lt. EPS ) ) GOTO 9 If( LY .GT. LMAX ) GOTO 999 IGES = IGES +1 IF( IGES.lt.ITERMAX ) goto 443 CC Iteration-Cycle goto 999 9 WRITE(6, *) ' --- Fini by Convergence --- ' 999 CONTINUE 22 FORMAT(3H { , F13.9, 3H , , F13.9, 3H }, ) 21 FORMAT(3H { , F13.9, 3H , , F13.9, 3H }} ) c write a Mathematica - file for visualization: Do 446, J=1,LY-1 WRITE(7, 22) Y(J,1), Y(J,2) 446 continue WRITE(7, 21) Y(LY,1), Y(LY,2) STOP END C en=energy, g=gradient, t=normed grad, h=Hessian C SUBROUTINE ENERGY(N,XX,en,grad) REAL*8 XX(2),grad(2),en,gnorm,x,y,x1,y1,x5,y5 Integer N ,I C mb - PES: Mueller/Brown x=XX(1) y=XX(2) I=N en= 1./100.*( 1 -170.*DEXP(-6.5*(0.5+x)**2+11.*(0.5+x)*(-1.5+y)-6.5*(-1.5+y)**2) 2 +15. *DEXP(0.7*(1.+x)**2+0.6*(1.+x)*(-1.+y)+0.7*(-1.+y)**2) 3 -100.*DEXP(-1.*x**2 - 10.*(-0.5+y)**2) 4 -200.*DEXP(-1.*(-1.+x)**2 - 10.*y**2) ) c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc x5= 0.5d0+x y5= -1.5d0 +y x1= 1.0d0+x y1= -1.0d0+y grad(1)= 1.0d0/100.0d0*( * -170.*DEXP(-6.5*x5**2+11.*x5*y5-6.5*y5**2)*(-13.*x5+11.*y5) 1 + 15.*DEXP(0.7*x1**2+0.6*x1*y1+0.7*y1**2) * (1.4*x1+0.6*y1) 2 -100.*DEXP(-1.*x**2 -10.*(-0.5+y)**2)* (-2.*x ) 3 -200.*DEXP(-1.*(-1.+x)**2 -10.*y**2)*(-2.*(-1.+x) ) ) grad(2)= 1.0d0/100.0d0*( * -170.*DEXP(-6.5*x5**2+11.*x5*y5-6.5*y5**2)*(11.*x5-13.*y5) 1 +15.*DEXP(0.7*x1**2+0.6*x1*y1+0.7*y1**2)*(0.6*x1+1.4*y1) 2 -100.*DEXP(-1.*x**2 -10.*(-0.5+y)**2)* ( - 20.*(-0.5+y)) 3 -200.*DEXP(-1.*(-1.+x)**2 - 10.*y**2)* ( - 20.*y) ) c gnorm=DSQRT( gx**2+gy**2 ) c tx=gx/gnorm c ty=gy/gnorm Return END C SUBROUTINE HESS(N,XX,HE) REAL*8 XX(2),HE(2,2),x,y,x1,y1,x5,y5 Integer N ,I I=N x=XX(1) y=XX(2) c x5= 0.5d0+x y5= -1.5d0 +y x1= 1.d0+x y1= -1.d0+y HE(1,1)=1.d0/100.*( 2210.*DEXP(-6.5*x5**2+11.*x5*y5-6.5*y5**2)+ 2 21.* DEXP(0.7*x1**2+0.6*x1*y1+0.7*y1**2)+ 3 200.*DEXP(-1.*x**2 - 10.*(-0.5+y)**2)+ 4 400.*DEXP(-1.*(-1.+x)**2 - 10.*y**2) - 5 170.*DEXP(-6.5*x5**2+11.*x5*y5-6.5*y5**2)*(-13.*x5+11.*y5)**2+ 1 15.*DEXP(0.7*x1**2+0.6*x1*y1+0.7*y1**2)*(1.4*x1+0.6*y1)**2 2 -100.*DEXP(-1.*x**2 -10.*(-0.5+y)**2 ) * 4.*x**2 - 3 200.*DEXP( -1.*(-1.+x)**2 - 10.*y**2 ) * 4.*(-1.+x)**2) HE(2,2)=1.d0/100.*( 2210.*DEXP(-6.5*x5**2+11.*x5*y5- 6.5*y5**2)+ 2 21.*DEXP(0.7*x1**2+0.6*x1*y1+0.7*y1**2)+ 3 2000.*DEXP(-1.*x**2 - 10.*(-0.5+y)**2)+ 4 4000.*DEXP(-1.*(-1.+x)**2 - 10.*y**2) - 5 170.*DEXP(-6.5*x5**2+11.*x5*y5- 6.5*y5**2)*(11.*x5-13.*y5)**2+ 1 15.*DEXP(0.7*x1**2+0.6*x1*y1+0.7*y1**2)* (0.6*x1+1.4*y1)**2 2 -100.*DEXP(-1.*x**2 - 10.*(-0.5+y)**2)*(20.*(-0.5+y))**2- 3 200.*DEXP(-1.*(-1.+x)**2 - 10.*y**2)*(20.*y)**2) HE(2,1)=1.d0/100.*(-1870.*DEXP(-6.5*x5**2+11.*x5*y5 -6.5*y5**2)+ 2 9.*DEXP(0.7*x1**2+0.6*x1*y1+0.7*y1**2)- 3 170.*DEXP(-6.5*x5**2+11.*x5*y5-6.5*y5**2)*(11.*x5-13.*y5)* 5 (-13.*x5+11.*y5)+ 6 15.*DEXP(0.7*x1**2+0.6*x1*y1+0.7*y1**2)* 7 (1.4*x1+0.6*y1)*(0.6*x1+1.4*y1) - 8 100.*DEXP(-x**2 - 10.*(-0.5+y)**2)*(-20.*(-0.5+y))*(-2.*x)- 9 200.*DEXP(-(-1.+x)**2 - 10.*y**2)*(-20.*y)*(-2.*(-1.+x)) ) HE(1,2)=HE(2,1) Return END SUBROUTINE HeighestPoint(N,Z,LZ,ETA,Dh) REAL*8 Z(250,2),X0(2),VV(2),ETA,R,Dh,en REAL*8 ProGrad(2),Tang(2),TN,ProMat(2,2) Integer I,J,N , JJ ,LZ c Projector at SP of Chain c R=-100000000.0d0 Dh=1.0d20 I=1 1 Continue If(I.gt.LZ) GOTO 99 Do 2 J=1,N 2 X0(J)= Z(I,J) Call Energy(N,X0,en,VV) If(en.gt.R) Then R=en I=I+1 Else GOTO 3 ENDIF GOTO 1 3 continue Write(6,*) 'maximales I ' , I Do 21 J=1,N 21 Tang(J)= Z(I,J)- Z(I+1,J) TN=0.0d0 Do 22 J=1,N 22 TN= TN+ Tang(J)*Tang(J) TN=DSQRT(TN) Do 23 J=1,N 23 Tang(J)=Tang(J)/TN C Projektionsoperator Do 24 J=1,N Do 24 JJ=1,N ProMat(J,JJ)=0.0d0 24 Continue Do 25 J=1,N 25 ProMat(J,J)=1.0d0 Do 26 J=1,N Do 261 JJ=1,N ProMat(J,JJ)= ProMat(J,JJ) - Tang(J)*Tang(JJ) 261 continue 26 continue C Projektion Gradient-Schritt Do 27 J=1,N ProGrad(J)=0.0D0 Do 271 JJ=1,N 271 ProGrad(J)=ProGrad(J) + VV(JJ)*ProMat(J,JJ) 27 Continue cc kein Schritt zur Neuen Kette: ( Proj ) grad TN=0.0d0 Do 28 J=1,N 28 TN= TN+ ProGrad(J)*ProGrad(J) Dh= ETA*DSQRT(TN) 99 continue Return END SUBROUTINE NEWCURVE(N,L,Y,Z,ETA,M,R) REAL*8 Y(250,2),Z(250,2),X(2),X0(2),VV(2),ETA,R,R2,D REAL*8 ProGrad(2),Tang(2),TN,ProMat(2,2),Rs(2),Pro2(2) REAL*8 HE(2,2),en ,Ad(2,2) Integer N ,I,J,K,M, method, L , JJ c if search direction, r, is input: c OPEN(19,FILE='rsearch.dat') c rewind 19 c Do 1119,J=1,N c READ(19,*) Rs(J) c1119 continue c 2 Minima MB: (-0.558, 1.442) & (0.624, 0.028) Rs(1)= 0.624d0 +0.558d0 Rs(2)= 0.028d0 -1.442d0 TN=0.0d0 Do 32 J=1,N 32 TN= TN+ Rs(J)*Rs(J) TN=DSQRT(TN) Do 33 J=1,N 33 Rs(J)=Rs(J)/TN C Projektionsoperator Do 34 J=1,N Do 34 JJ=1,N ProMat(J,JJ)=0.0d0 34 Continue Do 35 J=1,N 35 ProMat(J,J)=1.0d0 Do 36 J=1,N Do 361 JJ=1,N ProMat(J,JJ)= ProMat(J,JJ) - Rs(J)*Rs(JJ) 361 continue 36 continue C ETA=0.01d0 R2=R*R c Do 1 J=1,N 1 X0(J)= Y(1,J) PRINT*,' Methode: ' PRINT*,' 1:IRC, 2:IRC-Projector, 3:RGF-Projektor ' PRINT*,' 4:GE-Projector, 5:GE-Projektor2 ' C C READ(5,*) method cccccccccccccccccccccccccccccc C GE mit Projektor P_Ag ( g ) c method =5 cccccccccccccccccccccccccccccc C GE mit Projektor P_g ( H g ) c method =4 ccccccccccccccccccccccccccc c c method =3 c c C RGF mit Projektor P_r ccccccccccccccccccccccccccc c method = 2 C IRC mit Projektor P_t ccccccccccccccccccccccccccc c method =1 c IRC by DDRP: P = I ccccccccccccccccccccccccccc PRINT*,' Methode: ',method c Do 7 I=1,L Do 2 J=1,N 2 X(J)= Y(I,J) Do 5 K=1,M Call Energy(N,X,en,VV) c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC goto(351,352,353,354,355) method ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 351 Continue C C Gradient- Schritt fuer IRC Do 12 J=1,N 12 X(J)= X(J) - ETA*VV(J) D=0.0d0 Do 13 J=1,N 13 D=D+X(J)*X(J) If(D.GT.R2)THEN D=DSQRT(D) Do 14 J=1,N 14 X(J)= R* X(J)/D ENDIF GOTO 6 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 352 Continue C Tangent to curve If(I.lt.L) Then Do 21 J=1,N 21 Tang(J)= X(J)- Y(I+1,J) Endif If(I.eq.L) Then Do 211 J=1,N 211 Tang(J)= -X(J)+ Y(I-1,J) Endif c TN=0.0d0 Do 22 J=1,N 22 TN= TN+ Tang(J)*Tang(J) TN=DSQRT(TN) Do 23 J=1,N 23 Tang(J)=Tang(J)/TN C Projection operator Do 24 J=1,N Do 24 JJ=1,N ProMat(J,JJ)=0.0d0 24 Continue Do 25 J=1,N 25 ProMat(J,J)=1.0d0 Do 26 J=1,N Do 261 JJ=1,N ProMat(J,JJ)= ProMat(J,JJ) - Tang(J)*Tang(JJ) 261 continue 26 continue C WRITE(6, *) ' - ProMat - ', ((ProMat(J,JJ),j=1,N),JJ=1,N) C Projection: Gradient-Step Do 27 J=1,N ProGrad(J)=0.0D0 Do 271 JJ=1,N 271 ProGrad(J)=ProGrad(J) + VV(JJ)*ProMat(J,JJ) 27 Continue ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Step to new chain: ( Proj )* grad Do 28 J=1,N 28 X(J)= X(J) - ETA*ProGrad(J) CC D=0.0d0 Do 29 J=1,N 29 D=D+X(J)*X(J) If(D.GT.R2) THEN D=DSQRT(D) Do 30 J=1,N 30 X(J)= R* X(J)/D ENDIF GOTO 6 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc C Newton trajectory, with RGF Projector from the beginning 353 Continue C Projection Gradient Step Do 37 J=1,N ProGrad(J)=0.0D0 Do 371 JJ=1,N 371 ProGrad(J)=ProGrad(J) + VV(JJ)*ProMat(J,JJ) 37 Continue c Step to new chain: Proj * grad Do 38 J=1,N 38 X(J)= X(J) - ETA * ProGrad(J) GOTO 6 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C GE with Projector P_g ( H g ) 354 Continue Call HESS(N,X,HE) c WRITE(6, *) ' HesseMat ', ((HE(J,JJ),j=1,N),JJ=1,N) TN=0.0d0 Do 422 J=1,N 422 TN= TN+ VV(J)*VV(J) TN=DSQRT(TN) Do 423 J=1,N 423 Tang(J)=VV(J)/TN C use Tang for normed Gradient C Projection operator Do 424 J=1,N Do 424 JJ=1,N ProMat(J,JJ)=0.0d0 424 Continue Do 425 J=1,N 425 ProMat(J,J)=1.0d0 Do 46 J=1,N Do 45 JJ=1,N ProMat(J,JJ)= ProMat(J,JJ) - Tang(J)*Tang(JJ) 45 continue 46 continue C form H * g Do 47 J=1,N ProGrad(J)=0.0D0 Do 48 JJ=1,N 48 ProGrad(J)=ProGrad(J) + VV(JJ)*HE(J,JJ) 47 Continue C Projection - Step Do 41 J=1,N Pro2(J)=0.0D0 Do 42 JJ=1,N 42 Pro2(J)=Pro2(J) + ProGrad(JJ)*ProMat(J,JJ) 41 Continue c c Step to new chain: - Proj1 * grad Do 43 J=1,N 43 X(J)= X(J) - ETA * Pro2(J) GOTO 6 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 355 Continue C GE with Projector P_Ag ( g ) Call HESS(N,X,HE) 360 Continue C 2-Dim adjoint Ad is trivial: CCCCCCCCCCCCCCC Ad(1,1)= HE(2,2) Ad(2,1)= -HE(2,1) Ad(1,2)= -HE(2,1) Ad(2,2)= HE(1,1) cCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC c Form Adj * Gr Do 57 J=1,N ProGrad(J)=0.0D0 Do 58 JJ=1,N 58 ProGrad(J)=ProGrad(J) + VV(JJ)*Ad(J,JJ) 57 Continue TN=0.0d0 Do 522 J=1,N 522 TN= TN+ ProGrad(J)*ProGrad(J) IF(TN .lt. 1.1d-11) TN=1.0d0 TN=DSQRT(TN) Do 523 J=1,N 523 Tang(J)=ProGrad(J)/TN c use Tang for normed vector (Adj . Gradient) c Projection operator Do 524 J=1,N Do 524 JJ=1,N ProMat(J,JJ)=0.0d0 524 Continue Do 525 J=1,N 525 ProMat(J,J)=1.0d0 Do 56 J=1,N Do 55 JJ=1,N ProMat(J,JJ)= ProMat(J,JJ) - Tang(J)*Tang(JJ) 55 continue 56 continue c C Projection - Step Do 51 J=1,N ProGrad(J)=0.0D0 Do 52 JJ=1,N 52 ProGrad(J)=ProGrad(J) + VV(JJ)*ProMat(J,JJ) 51 Continue c Step to new chain: - Proj1 * grad c Do 53 J=1,N 53 X(J)= X(J) - ETA * ProGrad(J) GOTO 6 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 6 CONTINUE 5 CONTINUE Do 116 J=1,N 116 Z(I,J)= X(J) 7 CONTINUE c IF(method.eq.4) Then c Do 117 J=1,N c117 Z(1,J)= X0(J) c Endif RETURN END SUBROUTINE HOMOGEN(N,Y,LY,Z,LZ,EPS, LMAX) REAL*8 Y(250,2),Z(250,2),Y1(250,2),D(250),S,DD,EPS Integer M1,L, LZ, LY, I1,K,I,J,N,LY1, amm, LMAX C L Length of Chain L=LZ-1 Do 2 I=1,L S=0.0d0 I1=I+1 Do 1 J=1,N 1 S= S+(Z(I1,J)-Z(I,J))**2 2 D(I)= DSQRT(S) Do 3 J=1,N 3 Y1(1,J)=Z(1,J) LY1=1 S=0.0d0 Do 8 I=1,L I1=I+1 DD=D(I) If((DD.GE.EPS).AND.(S.GT. 0.0d0))THEN LY1=LY1+1 IF(LY1.gt. LMAX) GOTO 8 Do 4 J=1,N 4 Y1(LY1,J)=Z(I,J) ENDIF If(DD.GE.EPS) THEN amm= AINT(DD/EPS) M1=INT(amm) If(DD.GT.M1*EPS) M1=M1+1 Do 5 K=1,M1 LY1=LY1+1 IF(LY1.gt. LMAX) GOTO 8 Do 5 J=1,N 5 Y1(LY1,J)= ((M1-K)*Z(I,J)+ K*Z(I1,J))/M1 S=0.0d0 GOTO 8 ENDIF S=S+DD If(S.GT.EPS)THEN LY1=LY1+1 IF(LY1.gt. LMAX) GOTO 8 Do 6 J=1,N 6 Y1(LY1,J)= Z(I,J) S=DD ENDIF If(I.EQ.L)THEN LY1=LY1+1 Do 7 J=1,N 7 Y1(LY1,J)= Z(I1,J) ENDIF 8 CONTINUE Do 9 I=1,LY1 Do 9 J=1,N LY=I 9 Y(I,J)=Y1(I,J) RETURN END SUBROUTINE HAUSD(N,Y,LY,Z,LZ,D) REAL*8 Y(250,2),Z(250,2),D,D1 Integer N,LY,LZ CALL HAUSD2(N,Y,LY,Z,LZ,D1) CALL HAUSD2(N,Z,LY,Y,LY,D) If(D1.GT.D) D=D1 D=DSQRT(D) RETURN END SUBROUTINE HAUSD2(N,G,LG,H,LH,D) REAL*8 G(250,2),H(250,2),D,D1,S ,P ,S1 Integer N,LG,LH, I,J,K,I1 D=0.0d0 Do 4 J=1,LH Do 3 I=2,LG I1=I-1 P=0.0d0 S1=0.0d0 Do 1 K=1,N P = P +(H(J,K)-G(I,K))* (G(I1,K)-G(I,K)) 1 S1= S1+(G(I,K)-G(I1,K))**2 If((P.LE.0.0d0).OR.(S1.EQ.0.0d0)) P=0.0d0 If(S1.GT.0.0d0) P= P/S1 If(P.GT.1.0d0) P=1.0d0 S=0.0d0 Do 2 K=1,N 2 S= S+(G(I,K) +P*(G(I1,K)-G(I,K)) - H(J,K))**2 If((I.EQ.2).OR.(S.LT.D1)) D1=S 3 If(D1.LE.D) GOTO 4 D=D1 4 CONTINUE RETURN END