PROGRAM lj7sc4
C CCC for LJ_7 potential N=21 CCCCCCCCCCCCCCCCCCCCCCCCCC
C
C W.Quapp  30.06.2004  for RGF=Newton Trajectory
C GS nach: Baron Peters et al. JCP 120 (2004) 7877-7886
c script version - modular structure 
c  
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C EPREV=energy,V=gradient
c SUBROUTINE ENERGY(N,X,EPREV,V)                               
c 
       REAL*8 MASS,EPREV,X(21),V(21),gnorm
       INTEGER J1, MAX, N
       OPEN(13,FILE='lj7point.dat') 
       OPEN(14,FILE='lj7ener.dat') 
       OPEN(15,FILE='lj7grad.dat') 
       N=21
       rewind 13 
       READ(13,*)(X(J1),J1=1,N)
c
C********************************************************************
C  This is a program to find the gradient 
C  for a given LJ geometry, assume Argon, from D.Wales
C*********************************************************************
C
C   MASS=atomic mass in amu
C   Input parameters:
C: Befehl nicht gefunden.
C
       MASS=39.95D0
       MAX= N/3 
c
C      PRINT*
C      PRINT*,'    Input coordinates'
C      PRINT*
C      READ(5,*) MAX 
C      MAX=4 
C      DO 30 J1=1,MAX
C       N=3*(J1-1)+1
C       READ(10,*,END=40) X(N),X(N+1),X(N+2)
C       WRITE(6,20) X(N),X(N+1),X(N+2)
C 20    FORMAT(1X,3(F20.10,3X))
C 30    CONTINUE
C 40    N=J1-1
C       PRINT*
C       PRINT '(A30,I6)','Number of points=',N
C       PRINT*
c
      CALL PAIRS(3*MAX, X, EPREV)
c
       rewind 14
       write(14,*) EPREV
c 
       WRITE(6,*)' energy= ',EPREV
C      PRINT*
C      PRINT '(A30)','Gradient of the Energy : '
      CALL NDIFF(MAX, X, V)
C     WRITE(*,50) (V(J1),J1=1,N)
C 50    FORMAT(1X,3F20.10)
      gnorm = 0.D0
      DO 60 J1=1,N
60    gnorm = gnorm + V(J1)*V(J1)
      gnorm = DSQRT(gnorm) 
cc      DO 70 J1=1,N  
cc70    V(J1)=  V(J1)/ gnorm 
c
       rewind 15
       WRITE(15,*)(V(I),I=1,N)
c       Return 
      Return 
      END
C
C*********************************************************************
C
C     THIS SUBROUTINE CALCULATES THE PAIR POTENTIAL
C                                        
C     LET'S CONSIDER A LENNARD-JONES SYSTEM
C      THE FORM OF THE POTENTIAL IS:
C          U(R)=4E[(S/R)**12-(S/R)**6], WHERE E IS IN UNITS OF ENERGY,
C          S IS A CONSTANT LENGTH, AND R IS THE RADIUS.
C                                                          
      SUBROUTINE PAIRS (N, X, POTEL)
      IMPLICIT NONE 
      INTEGER N, I, J, J1, J2, MAX
      DOUBLE PRECISION X(21), SIG, EPS, POTEL, R
      MAX=7
      SIG=3.4D0
      EPS=1.0D0
      POTEL=0.0D0
      DO 20 I=1,N/3-1
         J1=3*(I-1)+1
         DO 10 J=I+1,N/3
            J2=3*(J-1)+1
            R=(X(J1)-X(J2))**2+(X(J1+1)-X(J2+1))**2+(X(J1+2)-X(J2+2))**2
            POTEL=POTEL+4.D0*EPS*(SIG**12/R**6-SIG**6/R**3)
10       CONTINUE
20    CONTINUE
      RETURN
      END
C
C*************************************************************************
C
C  Subroutine NDIFF calculates the cartesian gradient
C
C*************************************************************************
C
      SUBROUTINE NDIFF(N, X, V)
      IMPLICIT NONE
      INTEGER N, J1, J2, J3, J4 
      DOUBLE PRECISION X(21), 
     1                 V(21), DIST, S6, S12,
     2                 R8(7,7), R10(7,7), 
     3                 R14(7,7), R16(7,7)
C 
C  Note that S6 and S12 subsume a factor of 4 epsilon
C
      PARAMETER (S6=6179.217664D0, S12=9545682.734D0)
C
C  Store distance matrices.
C
      DO J1=1,N
         R8(J1,J1)=0.0D0
         R10(J1,J1)=0.0D0
         R14(J1,J1)=0.0D0
         R16(J1,J1)=0.0D0
         DO J2=1,J1-1
            DIST=(X(3*(J1-1)+1)-X(3*(J2-1)+1))**2
     1          +(X(3*(J1-1)+2)-X(3*(J2-1)+2))**2
     2          +(X(3*(J1-1)+3)-X(3*(J2-1)+3))**2
            R8(J1,J2)=DIST**(-4)
            R10(J1,J2)=DIST**(-5)
            R14(J1,J2)=DIST**(-7)
            R16(J1,J2)=DIST**(-8)
            R8(J2,J1)=R8(J1,J2)
            R10(J2,J1)=R10(J1,J2)
            R14(J2,J1)=R14(J1,J2)
            R16(J2,J1)=R16(J1,J2)
         ENDDO
      ENDDO
C
C  First calculate the gradient analytically.
C
      DO J1=1,N
         DO J2=1,3
            J3=3*(J1-1)+J2
            V(J3)=0.0D0
            DO J4=1,N
               V(J3)=V(J3)-(12*S12*R14(J4,J1)-6*S6*R8(J4,J1))
     1              *(X(J3)-X(3*(J4-1)+J2))
            ENDDO
cc grad:           PRINT*,'V(',J3,')=',V(J3)
        ENDDO
      ENDDO
      RETURN
      END