PROGRAM lj7sc3
C CCC for LJ_7 potential N=21 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 
C W.Quapp  30.06.2004  for RGF=Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
c test for script version - modular structure 
c GS cycle for RP: NewPoint calculation by projected gradient
c    goes along reduced gradient to chain point
Ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       REAL*8 X(21),XN(21),grad(21),PrGrad(21),ProMat(21,21),eta,eps,D1
       Integer N,J,JJ,istatus,ITALL,itGS
ccc
       OPEN(11,FILE='lj7param.dat') 
       OPEN(12,FILE='lj7proj.dat') 
       OPEN(13,FILE='lj7point.dat') 
ccc      OPEN(14,FILE='lj7ener.dat') 
       OPEN(15,FILE='lj7grad.dat') 
       OPEN(44,FILE='protocol.txt',status='unknown',access='append')
c
       N=21
       istatus=0
       rewind 11
       read(11,*) L, LA, eta, eps ,ITALL, itGS
       itGS=itGS+1
       rewind 11
       write(11,*)L, LA, eta, eps ,ITALL, itGS
       write(44,*)L, LA, eta, eps ,ITALL, itGS

       rewind 12
       Do 30 J=1,N
30     READ(12,*)(ProMat(J,I),I=1,N)
       rewind 13
       READ(13,*)(X(I),I=1,N)
c       rewind 14
c       READ(14,*) energy
       rewind 15
       READ(15,*)(grad(I),I=1,N)
c
cccccccccccccccccccccccccccccccccccccccc
C   Newton trajectory with RGF Projector
C   Projektion of Gradient-Step
cccccccccccccccccccccccccccccccccccccccc
c
      Do 37 J=1,N
      PrGrad(J)=0.0D0 
      Do 36 JJ=1,N
36    PrGrad(J)=PrGrad(J) + grad(JJ)*ProMat(J,JJ)
37    Continue
c   Step:  eta  Proj * grad 
      Do 38 J=1,N
38    XN(J)= X(J) - ETA * PrGrad(J)  
      rewind 13       
      WRITE(13,*)(XN(I),I=1,N)

      WRITE(44,*) ' new point '
      WRITE(44,*) (XN(I),I=1,N)
      close(44) 

c   Simple test of difference to former point
c   should be zero if PrGrad is zero      
      D1=0.0d0
      Do 40, I=1,N
40    D1=D1 + (XN(I)-X(I))**2
      D1=DSqrt(D1) 
      WRITE(6,*) 'threshold  ', D1

c   Istatus controls the script loops 
      IF(D1.lt.eps) Istatus=1 
      If(LA.gt.L)   Istatus=10 
      call exit(istatus)
      STOP
      END