      PROGRAM NEB
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAXDIM=50)
      DIMENSION Q(MAXDIM,2),E(MAXDIM),G(MAXDIM,2),F(2)
      DIMENSION TANG(2),DELTAQ(2),Ve(2),Fold(2)
      DIMENSION X(2), gr(2)
C
C  A NEB SCHEME by Euler integrator 
C          W.Quapp 01.03.2010 
c                  (using a program of Antoni Aguilar- Mogas)  
c  for the modified NFK surface
C  nn      DIMENSION OF THE SYSTEM = 2
C  NVER    NUMBER OF IMAGES OF THE STRING (INCLUDING FIRST AND LAST)
C  NITER   MAXIMUM NUMBER OF ITERATIONS over whoule chain
C  DK      MAXIMUM VALUE OF THE SPRING CONSTANT
C  DELTAT  DISPLACEMENT CONSTANT
C  SYSMAS  MASS OF THE SYSTEM = 1 
C  NINT    NUMBER OF INTEGRATION STEPS for one and the same node
c                                         by a direct force step 
C  FMIN    MAXIMUM FORCE NORM AT WHICH THE PROCESS STOPS
c      OPEN(UNIT=8,FILE='force.dat')
      OPEN(UNIT=9,FILE='CHAINnfk03.dat')
      OPEN(UNIT=10,FILE='StartNFK03.dat')
      nn   =    2 
      NVER =   50
c      Write(6,*) 'give NITER now -- number of iterations '
c      Read(5,*)  NITER 
      NITER=  750
      NINT =    2
      DELTAT = 0.0125D0
      FMIN =   1.0D-11
      DK =     0.0D0
      SYSMAS = 1.0D0
C MINIMA
c usual karplus:
c        Q(1,1)=-2.71D0
c        Q(1,2)=0.15D0
c        Q(NVER,1)=2.71D0
c        Q(NVER,2)=-0.15D0
c half karplus with r=0.03:
c behind SP
c Min       
        Q(1,1)= 2.85D0
        Q(1,2)=-0.15D0        
c behind SP
        Q(NVER,1)=-0.61d0
        Q(NVER,2)= 0.63d0
c
C INTERMEDIATE VERTEXES for initial chain
      DO 100 J=1,nn
        DELTAQ(J)=(Q(NVER,J)-Q(1,J))/(NVER-1)
 100  CONTINUE
      DO 107 IVER=2,NVER
       DO 110 J=1,nn
         Q(IVER,J)=Q(1,J)+(IVER-1)*DELTAQ(J)
 110    CONTINUE
 107    CONTINUE
c    Mma input of initial chain
        write(10,*)' { '
        DO 108 IVER=1,NVER
        write(10,*)' { ',Q(IVER,1),'  ,  ',Q(IVER,2), '},'
 108    CONTINUE
        write(10,*)' { ',Q(NVER,1),'  ,  ',Q(NVER,2), '} '
        write(10,*)' } '
 
       DO 10 IVER=1,NVER
        NCOMPU=0
C ----- SLIGHTLY DISPLACE THE CURVE TO THE LEFT(-) OR RIGHT(+) -----
c       Q(IVER,1)=Q(IVER,1)+((IVER-1)*DISINI*(NVER-IVER))
        X(1)=Q(IVER,1)
        X(2)=Q(IVER,2)
        CALL PES(X,Gr,E(IVER))
        G(IVER,1)=gr(1)
        G(IVER,2)=gr(2)
        NCOMPU=NCOMPU+1
 10   CONTINUE
C VARIABLE SPRING CONSTANT --  not used 
c      IF (E(NVER).GE.E(1)) THEN
c          EREF=E(NVER)
c          EMAX=E(1)
c      ELSE
c          EREF=E(1)
c          EMAX=E(NVER)
c      ENDIF
c
C BEGIN LOOP *******************************************************
       DO 999 KITER=1,NITER
        FNMAX=0.0D00
c        DO 20 IVER=2,NVER-1
c          IF (E(IVER).GT.EMAX) EMAX=E(IVER)
c 20     CONTINUE
       write(*,*)
       DO 30 IVER=2,NVER-1
       NINTende=1 
       If(KITER .gt. NITER-4) NINTende=NINT
    
          DO 111 INi=1,NINTende

C SEGMENT TANGENT TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
           if( IVER.eq.2 ) then
             DO  J=1,nn
                TANG(J)=-(Q(IVER-1,J) - Q(IVER,J))
             enddo
             goto 33 
           else                    
            DO  J=1,nn
                TANG(J)= -(Q(IVER-1,J) - Q(IVER+1,J))
            enddo
           endif  
c
 33       Tnorm=DSQRT(SCAPRO(TANG,TANG,nn))
          Do J=1,nn
            TANG(J)= TANG(J)/Tnorm
          enddo
c pes with new point and gradient
          X(1)=Q(IVER,1)
          X(2)=Q(IVER,2)
               CALL PES(X,gr,ener)
          E(IVER)  =ener     
          G(IVER,1)=gr(1)
          G(IVER,2)=gr(2)
          NCOMPU=NCOMPU+1
C SPRING CONSTANT  --  not used
c          IF (E(IVER).GT.E(IVER-1)) THEN
c             EI=E(IVER)
c          ELSE
c             EI=E(IVER-1)
c          ENDIF
c          IF (EI.GT.EREF) THEN
c             DKI=DK-DKMIN*((EMAX-EI)/(EMAX-EREF))
c          ELSE
c             DKI=DK -DKMIN
c          ENDIF
c
C FORCE CALCULATION
cc Projector is 1-( dyadic product of tangent vector ) 
cc Projector (here 2-dimensional)
         DO 49 J=1,nn
          F(J)=-G(IVER,J)+TANG(J)*
     &     ( G(IVER,1)*TANG(1)+G(IVER,2)*TANG(2) )
 49       CONTINUE
cc add  parallel force if DK is not zero:
         XX1=0.0d0
         Do J=1,nn
           XX1=XX1 + (Q(IVER+1,J)-Q(IVER,J))**2
         enddo
         XX2=0.0d0
         Do J=1,nn
           XX2=XX2 + (Q(IVER,J)-Q(IVER-1,J))**2
         enddo
         XX1=DSQRT(XX1)                             
         XX2=DSQRT(XX2)
         DO 50 J=1,nn
           F(J)=F(J) + DK*(XX2-XX1)*TANG(J)
 50      CONTINUE

         FNORM=SCAPRO(F,F,nn)
         IF (FNORM.GT.FNMAX) FNMAX=FNORM
c for start Euler
         If(IVER.eq.2) then
           do J=1,nn
             Ve(J)= 0.0d0
           enddo
         endif                           
          
C PRINT FORCE VECTOR OF ITERATION
c          IF (KITER.EQ.NITER) THEN
c             WRITE(8,9000) E(IVER),(Q(IVER,J),J=1,nn),
c     &                     (F(J),J=1,nn)
c          ENDIF
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c wq:
c Euler integrator - see Sheppard et al. 2008 JCP 128, p.134106
c its use depends from NINT
c
          Fnorm=DSQRT(SCAPRO(F,F,nn))
          If( Fnorm .gt. 1.0D-33) then 
           do J=1,nn
             Fold(J)=F(J)
             F(J)= F(J)/Fnorm
           enddo
          endif                      
          Test=SCAPRO(Ve,F,nn)
          do J=1,nn
            Ve(J)= Test*F(J)
          enddo
          If(Test.lt.0.0d0) then
            do J=1,nn
              Ve(J)=0.0d0
            enddo  
          endif
           DO 70 J=1,nn
            if(NINTende .gt. 1) Then
             Q(IVER,J)= Q(IVER,J) + DELTAT*Fold(J)
            else 
             Q(IVER,J)= Q(IVER,J) + Ve(J)*DELTAT
            endif
             Ve(J)    = Ve(J) + DELTAT*F(J)
 70       CONTINUE
c  Euler integrator - end -----------------------------           
C CONVERGENCE CRITERIA
        If(INi .eq. NINT)  
     &  WRITE(6,*) KITER, IVER, INi,'  FNMAX   ',FNMAX 
        IF (FNMAX.LE.FMIN) goto 777

111     CONTINUE
 30     CONTINUE 
999     CONTINUE

  777  continue
C END OF LOOP. PRINT RESULTS to Mma file 
       write(9,*)' { ' 
      DO 888 IVER=1,NVER
c         WRITE(6,*) E(IVER),(Q(IVER,J),J=1,nn)
          write(9,*)' { ',Q(IVER,1),'  ,  ',Q(IVER,2), '},'
 888  CONTINUE
          write(9,*)' { ',Q(NVER,1),'  ,  ',Q(NVER,2), '} '
          write(9,*)' } '
      WRITE(6,*) ' FMIN = ', FMIN
      WRITE(6,*) '#TOTAL NUMBER OF E & G COMPUTATIONS:',NCOMPU
      WRITE(6,*) '#TOTAL NUMBER OF ITERATIONS:',KITER
c     CLOSE(8)
      CLOSE(9)
      CLOSE(10)
 9000 FORMAT(5(1X,F12.8))
      END
C
      SUBROUTINE PES(COORD,GRAD,ENERGY)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION COORD(2),GRAD(2)
C 
C E. NERIA, S. FISCHER, M. KARPLUS JCP 105(5), 1902(1996)
C     R is modification parameter by M.Hirsch, CPL 2004
C
       R=0.03
       X=COORD(1)
       Y=COORD(2)
       EXPLU =EXP(-(X+3.)**2-Y**2)
       EXMIN =EXP(-(X-3.)**2-Y**2)
       ENERGY=R*(X**2+Y**2)**2+X*Y-9.*EXMIN-9.*EXPLU
       GRAD(1)=4.*R*X*(X**2+Y**2)+Y+18.*(X-3.)*EXMIN+18.*(X+3.)*EXPLU
       GRAD(2)=4.*R*Y*(Y**2+X**2)+X+18.*Y*EXPLU +18*Y*EXMIN
      RETURN
      END
C
      FUNCTION SCAPRO(AVEC,BVEC,N)
C
C COMPUTE THE SCALAR PRODUCT BETWEEN THE VECTORS AVEC AND BVEC
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION AVEC(N), BVEC(N)
      SUM = 0.0D0
      DO 1 I = 1,N
         SUM = SUM + AVEC(I) * BVEC(I)
    1 CONTINUE
      SCAPRO = SUM
      RETURN
      END

