        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
