MODULE DEFINITIONS
     !
     ! All constants and variables defined in this module
     ! will be known by all program units that USE this module
     !
     INTEGER, PARAMETER :: LONG=selected_real_kind(9,99)
     INTEGER, PARAMETER :: NX=1000   ! Max number of particles allowed
     REAL(KIND=LONG), PARAMETER :: ZERO = 0.0_LONG
     REAL(KIND=LONG), PARAMETER :: UNO  = 1.0_LONG
     REAL(KIND=LONG), PARAMETER :: PI   = ACOS(-UNO)
     REAL(KIND=LONG), PARAMETER :: SQ3  = SQRT(3.0_LONG)  
     REAL(KIND=LONG), DIMENSION(3), PARAMETER :: &
                      R0=(/ UNO/SQ3,UNO/SQ3,UNO/SQ3 /)
     REAL(KIND=LONG)  Delta_t, dt, dt2, dt22, time, time_samp, dts
     REAL(KIND=LONG)  ENK,ENP
     REAL(KIND=LONG), DIMENSION(2*NX,3) :: U
     REAL(KIND=LONG), DIMENSION(NX,3) :: F
     ! 
     ! U(I,J) contains the three Cartesian components (J=1,2,3) of 
     !        the particles displacements for 1 <= I <= N and 
     !        the three Cartesian components of the particles
     !        velocities for N+1 <= I <= 2*N
     ! F(I,J) contains the three Cartesian components (J=1,2,3) of 
     !        the forces on the I-th particle (1 <= I <= N) 
     !
     INTEGER I,J,K,M,N,Ns
     LOGICAL SAMPLE
END MODULE DEFINITIONS
!=======================================================================
PROGRAM MAIN 
    USE DEFINITIONS
    N  = 1                             ! Number of particles
    Delta_t = 200*PI                     ! Integration time
    dt = 0.001_LONG                    ! Time step
    dt2 = dt / 2.0_LONG
    dt22 = dt*dt / 2.0_LONG            !
    dts = 1000 * dt                     ! Sampling time step
    Ns = INT(Delta_t/dt)               ! Number of integration steps
    CALL INIT                          ! Initial conditions
    PRINT*, '# U(0),V(0) =', U(1,:), U(2,:)
    PRINT*, '# |V(0)|    =', (U(2,1)+U(2,2)+U(2,3))/SQ3
    CALL FORCES                        ! Initial forces
    time = 0.0_LONG
    time_samp = time
    CALL OPEN_FILES
!=======================================================================
!    
    DO II=1,Ns                         ! Loop over number of int. steps
      time = time + dt
      CALL INTEGR                      ! Perform one step of integration
      SAMPLE = time > time_samp        ! Sampling switch
      IF ( SAMPLE ) THEN               ! Sampling time?
         CALL COMPUTE_OBSERVABLES      ! Compute relevant observables
         CALL WRITE_DATA               ! Write required data on file
         time_samp = time_samp + dts   ! Next sampklig time
         SAMPLE = .FALSE.              ! Reset sampling switch
      END IF
    END DO                             ! End loop over int. steps
!=======================================================================
END PROGRAM MAIN 
!=======================================================================
! Procedures externes 
!
SUBROUTINE INIT
    USE DEFINITIONS
    U = ZERO
    DO J=1,3
      U(2,J) = 1.0_LONG * RAND()
    END DO  
END SUBROUTINE INIT
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE INTEGR
    USE DEFINITIONS
    !
    ! This subroutine implements your integration algorithm 
    ! Simple velocity Verlet
    !
    ! Advance displacement and half-time-step velocities
    DO I=1,N
       DO J=1,3
           U(I,J)   = U(I,J)   + dt * U(2*I,J) + dt22 * F(I,J)
           U(2*I,J) = U(2*I,J) + dt2* F(I,J) 
       END DO
    END DO
    CALL FORCES ! New forces 
    ! The remaining half time-step for velocities
    DO I=1,N
       DO J=1,3
           U(2*I,J) = U(2*I,J) + dt2 * F(I,J) 
       END DO 
    END DO
END SUBROUTINE INTEGR
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE FORCES
    USE DEFINITIONS
    REAL(KIND=LONG) r
    ! This subroutine compute the forces
    r = (R0(1) + U(1,1))**2 + &
        (R0(2) + U(1,2))**2 + &
        (R0(3) + U(1,3))**2
    r = SQRT(r) 
    DO J=1,3
       F(1,J) = -(1.0_LONG - 1.0_LONG/r)*(R0(J)+U(1,J))
    END DO
END SUBROUTINE FORCES   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE ENERGY
    USE DEFINITIONS
    ! This subroutine computes the system energy
    REAL(KIND=LONG) r
    r = (R0(1) + U(1,1))**2 + &
        (R0(2) + U(1,2))**2 + &
        (R0(3) + U(1,3))**2
    r = SQRT(r)
    ENK = ZERO
    DO J=1,3
       ENK = ENK + U(2,J)**2
    END DO
    ENK = ENK / 2.0_LONG
    ENP = 0.5_LONG*(r-1.0_LONG)**2 
END SUBROUTINE ENERGY   
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE COMPUTE_OBSERVABLES
    USE DEFINITIONS
    ! This subroutine computes the desired observables
    CALL ENERGY
END SUBROUTINE COMPUTE_OBSERVABLES  
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE OPEN_FILES
    OPEN(1,FILE='XV.DAT')
END SUBROUTINE OPEN_FILES
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE WRITE_DATA
    USE DEFINITIONS
    ! This subroutine writes your data on file
    PRINT*, time, ENK, ENP, ENK + ENP
    WRITE (1,'(10(G12.6,3X))') time,R0+U(1,:), U(2,:)
END SUBROUTINE WRITE_DATA  

