qfqwf

MODULE CONSTANTS
    IMPLICIT NONE
    INTEGER, PARAMETER :: dp                       = SELECTED_REAL_KIND(16)

    REAL(KIND = dp), PARAMETER :: PI               = 3.1415926536_dp
    REAL(KIND = dp), PARAMETER :: BOLTZMANN        = 1.38064852E-23_dp ! J/K

    REAL(KIND = dp), PARAMETER :: EPS              = 1.65E-21_dp ! J
    REAL(KIND = dp), PARAMETER :: M                = 6.6335209E-26_dp ! kg
    REAL(KIND = dp), PARAMETER :: SIGMA            = 3.4E-10_dp ! m
    REAL(KIND = dp), PARAMETER :: LIQUID_DENSITY   = 1401 ! kg / m^3

    ! Units - not yet used
    REAL(KIND = dp), PARAMETER :: UNIT_LENGTH      = SIGMA ! m
    REAL(KIND = dp), PARAMETER :: UNIT_TEMPERATURE = EPS / BOLTZMANN ! K
    REAL(KIND = dp), PARAMETER :: UNIT_ENERGY      = EPS ! J
    REAL(KIND = dp), PARAMETER :: UNIT_TIME        = SIGMA * SQRT(M / EPS) ! 1/s=Hz
    REAL(KIND = dp), PARAMETER :: UNIT_FORCE       = EPS / SIGMA ! N
    REAL(KIND = dp), PARAMETER :: UNIT_PRESSURE    = EPS / (SIGMA ** 3) ! N * m**-2
    REAL(KIND = dp), PARAMETER :: UNIT_SPEED       = SQRT(EPS / M) ! m/s
END MODULE

MODULE SIMULATION_BOX
    USE CONSTANTS
    INTEGER :: N

    REAL(KIND = dp) :: BOX ! m
    REAL(KIND = dp) :: DENSITY ! 1 / (m ** 3)
    REAL(KIND = dp) :: VOLUME ! m ** 3

    REAL(KIND = dp), DIMENSION(:), ALLOCATABLE :: RX, RY, RZ ! m
    REAL(KIND = dp), DIMENSION(:), ALLOCATABLE :: VX, VY, VZ ! m / s
    REAL(KIND = dp), DIMENSION(:), ALLOCATABLE :: AX, AY, AZ ! m / (s ** 2)
    REAL(KIND = dp), DIMENSION(:), ALLOCATABLE :: FX, FY, FZ ! eV / (m ** 2)

    REAL(KIND = dp) :: RCUT
CONTAINS

    !
    ! Prints out pair distances
    !
    SUBROUTINE PRINT_PAIR_DISTANCES()
        IMPLICIT NONE

        REAL(KIND = dp) :: SCLX, SCLY, SCLZ
        INTEGER         :: I, J

        WRITE (*, *)
        WRITE (*, '(''DISTANCE MATRIX'')')
        WRITE (*, '(A)', ADVANCE = 'NO') "                   1"
        DO I = 2, N
            WRITE (*, '(I10)', ADVANCE = 'NO') I
        END DO

        WRITE (*, *)

        DO I = N, 1, -1
            WRITE (*, '(I10)', ADVANCE = 'NO') I

            DO J = 1, I - 1
                SCLX = (RX(J) - RX(I)) ** 2.0
                SCLY = (RY(J) - RY(I)) ** 2.0
                SCLZ = (RZ(J) - RZ(I)) ** 2.0

                WRITE (*, '(F10.5)', ADVANCE = 'NO') SQRT(SCLX + SCLY + SCLZ)
            END DO

            WRITE (*, *)
        END DO
    END SUBROUTINE

    !
    ! Returns true if the box needs position correction.
    ! This function is used when generating an initial configuration.
    !
    LOGICAL FUNCTION NEEDS_CORRECTION(THRESHOLD) RESULT(CORR)
        IMPLICIT NONE

        INTEGER         :: I, J
        REAL(KIND = dp) :: DIST
        REAL(KIND = dp) :: THRESHOLD

        DO I = 1, N - 1
            DO J = I + 1, N
                DIST = SQRT((RX(J) - RX(I)) ** 2 + (RY(J) - RY(I)) ** 2 + (RZ(J) - RZ(I)) ** 2)

                IF (DIST .LT. THRESHOLD) THEN
                    CORR = .TRUE.
                    return
                END IF
            END DO
        END DO

        CORR = .FALSE.
        return
    END FUNCTION

    !
    ! Allocates runtime arrays
    !
    SUBROUTINE INIT_BOX()
        INTEGER :: ALLOC_STATUS

        ALLOCATE ( VX(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP
        ALLOCATE ( VY(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP
        ALLOCATE ( VZ(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP

        ALLOCATE ( RX(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP
        ALLOCATE ( RY(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP
        ALLOCATE ( RZ(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP

        ALLOCATE ( AX(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP
        ALLOCATE ( AY(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP
        ALLOCATE ( AZ(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP

        ALLOCATE ( FX(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP
        ALLOCATE ( FY(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP
        ALLOCATE ( FZ(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP
    END SUBROUTINE

    !
    ! Cleanup routine
    !
    SUBROUTINE TERMINATE_BOX()
        DEALLOCATE (VX)
        DEALLOCATE (VY)
        DEALLOCATE (VZ)

        DEALLOCATE (RX)
        DEALLOCATE (RY)
        DEALLOCATE (RZ)

        DEALLOCATE (AX)
        DEALLOCATE (AY)
        DEALLOCATE (AZ)

        DEALLOCATE (FX)
        DEALLOCATE (FY)
        DEALLOCATE (FZ)
    END SUBROUTINE

    !
    ! Applies periodic boundary conditions
    !
    SUBROUTINE BOUNDARY_CONDITIONS()
        IMPLICIT NONE
        INTEGER :: I

        DO I = 1, N
            RX(I) = RX(I) - ANINT( RX(I) / BOX ) * BOX
            RY(I) = RY(I) - ANINT( RY(I) / BOX ) * BOX
            RZ(I) = RZ(I) - ANINT( RZ(I) / BOX ) * BOX
        END DO
    END SUBROUTINE
END MODULE

MODULE LENNARD_JONES
    USE CONSTANTS
    IMPLICIT NONE

    REAL(KIND = dp) :: VLRC, WLRC
    REAL(KIND = dp) :: V, VC, W ! eV

    REAL(KIND = dp), DIMENSION(:), ALLOCATABLE :: DFX, DFY, DFZ
CONTAINS

    !
    ! Initializes long-range lennard-jones correction terms
    !
    SUBROUTINE INIT_LJ()
        USE SIMULATION_BOX
        USE CONSTANTS
        IMPLICIT NONE

        REAL(KIND = dp) :: SR3, SR9
        INTEGER :: ALLOC_STATUS, I

        SR3  = ( 1.0_dp / RCUT ) ** 3
        SR9  = SR3 ** 3
        VLRC = (  8.0_dp / 9.0_dp ) * PI * DENSITY * REAL( N ) * ( SR9 - 3.0_dp * SR3 )
        WLRC = ( 16.0_dp / 9.0_dp ) * PI * DENSITY * REAL( N ) * ( 2.0_dp * SR9 - 3.0_dp * SR3 )

        ALLOCATE ( DFX(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP
        ALLOCATE ( DFY(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP
        ALLOCATE ( DFZ(N), STAT = ALLOC_STATUS )
        IF (ALLOC_STATUS /= 0) STOP

        DO I = 1, N
            DFX(I) = 0.0_dp
            DFY(I) = 0.0_dp
            DFZ(I) = 0.0_dp
        END DO
    END SUBROUTINE

    SUBROUTINE TERMINATE_LJ()
        DEALLOCATE(DFX)
        DEALLOCATE(DFY)
        DEALLOCATE(DFZ)
    END SUBROUTINE

    !
    ! Applies long-range corrections
    !
    SUBROUTINE CORRECT_LJ()
        V = V + VLRC
        W = W + WLRC
    END SUBROUTINE

    !
    ! Calculates lennard-jones interactions between particle pairs,
    ! computing the total value of potential, shifted potential and virial,
    ! as well as updating forces on each particle
    !
    SUBROUTINE CALCULATE_LJ()
        USE SIMULATION_BOX

        IMPLICIT NONE

        REAL(KIND = dp) :: RIJX, RIJY, RIJZ, RIJSQ
        REAL(KIND = dp) :: RCUTSQ
        REAL(KIND = dp) :: SR2, SR6, SR12
        REAL(KIND = dp) :: VIJ, WIJ
        REAL(KIND = dp) :: DFXIJ, DFYIJ, DFZIJ

        INTEGER :: NCUT, I, J

        RCUTSQ  = RCUT ** 2

        DO I = 1, N
            FX(I) = 0.0_dp
            FY(I) = 0.0_dp
            FZ(I) = 0.0_dp

            DFX(I) = 0.0_dp
            DFY(I) = 0.0_dp
            DFZ(I) = 0.0_dp
        END DO

        V = 0.0
        VC = 0.0
        W = 0.0
        NCUT = 0

        DO I = 1, N - 1
            DO J = I + 1, N
                RIJX = RX(J) - RX(I)
                RIJY = RY(J) - RY(I)
                RIJZ = RZ(J) - RZ(I)

                RIJX = RIJX - ANINT(RIJX / BOX) * BOX
                RIJY = RIJY - ANINT(RIJY / BOX) * BOX
                RIJZ = RIJZ - ANINT(RIJZ / BOX) * BOX

                RIJSQ = RIJX ** 2 + RIJY ** 2 + RIJZ ** 2

                IF (RIJSQ .LT. RCUTSQ) THEN
                    SR2 = 1.0_dp / RIJSQ
                    SR6 = SR2 * SR2 * SR2
                    SR12 = SR6 * SR6

                    VIJ = SR12 - SR6
                    WIJ = 2.0_dp * SR12 - SR6

                    ! Forces
                    FX(I) = FX(I) + WIJ * (RIJX / RIJSQ)
                    FY(I) = FY(I) + WIJ * (RIJY / RIJSQ)
                    FZ(I) = FZ(I) + WIJ * (RIJZ / RIJSQ)

                    FX(J) = FX(J) - WIJ * (RIJX / RIJSQ)
                    FY(J) = FY(J) - WIJ * (RIJY / RIJSQ)
                    FZ(J) = FZ(J) - WIJ * (RIJZ / RIJSQ)

                    ! Configurational temperature
                    DFXIJ = 48.0_dp * (RIJSQ - 14.0_dp * (RIJX ** 2)) / (RIJSQ ** 8) - 24.0_dp * (RIJSQ - 8.0_dp * (RIJX ** 2)) / (RIJSQ ** 5)
                    DFYIJ = 48.0_dp * (RIJSQ - 14.0_dp * (RIJY ** 2)) / (RIJSQ ** 8) - 24.0_dp * (RIJSQ - 8.0_dp * (RIJY ** 2)) / (RIJSQ ** 5)
                    DFZIJ = 48.0_dp * (RIJSQ - 14.0_dp * (RIJZ ** 2)) / (RIJSQ ** 8) - 24.0_dp * (RIJSQ - 8.0_dp * (RIJZ ** 2)) / (RIJSQ ** 5)

                    DFX(I) = DFX(I) + DFXIJ
                    DFY(I) = DFY(I) + DFYIJ
                    DFZ(I) = DFZ(I) + DFZIJ

                    DFX(J) = DFX(J) - DFXIJ
                    DFY(J) = DFY(J) - DFYIJ
                    DFZ(J) = DFZ(J) - DFZIJ

                    V = V + VIJ
                    W = W + WIJ

                    NCUT = NCUT + 1
                END IF
            END DO
        END DO

        ! Shifted Potential

        SR2 = 1.0 / RCUTSQ
        SR6 = SR2 * SR2 * SR2
        SR12 = SR6 ** 2
        VC = V - REAL(NCUT) * (SR12 - SR6)

        DO I = 1, N
            FX(I) = FX(I) * 24.0_dp
            FY(I) = FY(I) * 24.0_dp
            FZ(I) = FZ(I) * 24.0_dp
        END DO

        V  = V * 4.0_dp
        W  = W * 24.0_dp / 3.0_dp
        VC = VC * 4.0_dp

    END SUBROUTINE

END MODULE

MODULE GEAR
    USE CONSTANTS
    REAL(KIND = dp), PARAMETER :: GEAR0 = 19.0_dp / 120.0_dp
    REAL(KIND = dp), PARAMETER :: GEAR1 =  3.0_dp /   4.0_dp
    REAL(KIND = dp), PARAMETER :: GEAR3 =  1.0_dp /   2.0_dp
    REAL(KIND = dp), PARAMETER :: GEAR4 =  1.0_dp /  12.0_dp

    REAL(KIND = dp), DIMENSION(:), ALLOCATABLE :: BX, BY, BZ
    REAL(KIND = dp), DIMENSION(:), ALLOCATABLE :: CX, CY, CZ
CONTAINS

    !
    ! Allocation
    !
    SUBROUTINE INIT_GEAR()
        USE SIMULATION_BOX
        IMPLICIT NONE
        INTEGER :: STAT_ALLOC

        ALLOCATE (BX(N), STAT = STAT_ALLOC)
        IF (STAT_ALLOC /= 0) STOP
        ALLOCATE (BY(N), STAT = STAT_ALLOC)
        IF (STAT_ALLOC /= 0) STOP
        ALLOCATE (BZ(N), STAT = STAT_ALLOC)
        IF (STAT_ALLOC /= 0) STOP

        ALLOCATE (CX(N), STAT = STAT_ALLOC)
        IF (STAT_ALLOC /= 0) STOP
        ALLOCATE (CY(N), STAT = STAT_ALLOC)
        IF (STAT_ALLOC /= 0) STOP
        ALLOCATE (CZ(N), STAT = STAT_ALLOC)
        IF (STAT_ALLOC /= 0) STOP
    END SUBROUTINE

    !
    ! Cleanup routine
    !
    SUBROUTINE TERMINATE_GEAR()
        DEALLOCATE (BX)
        DEALLOCATE (BY)
        DEALLOCATE (BZ)

        DEALLOCATE (CX)
        DEALLOCATE (CY)
        DEALLOCATE (CZ)
    END SUBROUTINE

    !
    ! Implements the prediction part of Gear's 5-value algorithm
    !
    SUBROUTINE PREDICT(DT)
        USE SIMULATION_BOX
        IMPLICIT NONE

        REAL(KIND = dp), INTENT(IN) :: DT
        REAL(KIND = dp) :: C2, C3, C4

        INTEGER :: I

        C2 = DT * DT / 2.0_dp
        C3 = C2 * DT / 3.0_dp
        C4 = C3 * DT / 4.0_dp

        DO I = 1, N
            RX(I) = RX(I) + DT * VX(I) + C2 * AX(I) + C3 * BX(I) + C4 * CX(I)
            RY(I) = RY(I) + DT * VY(I) + C2 * AY(I) + C3 * BY(I) + C4 * CY(I)
            RZ(I) = RZ(I) + DT * VZ(I) + C2 * AZ(I) + C3 * BZ(I) + C4 * CZ(I)

            VX(I) = VX(I) + DT * AX(I) + C2 * BX(I) + C3 * CX(I)
            VY(I) = VY(I) + DT * AY(I) + C2 * BY(I) + C3 * CY(I)
            VZ(I) = VZ(I) + DT * AZ(I) + C2 * BZ(I) + C3 * CZ(I)

            AX(I) = AX(I) + DT * BX(I) + C2 * CX(I)
            AY(I) = AY(I) + DT * BY(I) + C2 * CY(I)
            AZ(I) = AZ(I) + DT * BZ(I) + C2 * CZ(I)

            BX(I) = BX(I) + DT * CX(I)
            BY(I) = BY(I) + DT * CY(I)
            BZ(I) = BZ(I) + DT * CZ(I)
        END DO

    END SUBROUTINE

    !
    ! Implements the corrector part of Gear's 5-value algorithm
    !
    SUBROUTINE CORRECT(DT)
        USE SIMULATION_BOX
        IMPLICIT NONE

        REAL(KIND = dp), INTENT(IN) :: DT

        REAL(KIND = dp) :: C2, C3, C4
        REAL(KIND = dp) :: CR, CV, CB, CC
        REAL(KIND = dp) :: AXI, AYI, AZI
        REAL(KIND = dp) :: CORRX, CORRY, CORRZ

        INTEGER :: I

        C2 = DT * DT / 2.0_dp
        C3 = C2 * DT / 3.0_dp
        C4 = C3 * DT / 4.0_dp

        CR = GEAR0 * C2
        CV = GEAR1 * C2 / DT
        CB = GEAR3 * C2 / C3
        CC = GEAR4 * C2 / C4

        DO I = 1, N
            AXI = FX(I) / M
            AYI = FY(I) / M
            AZI = FZ(I) / M

            CORRX = AXI - AX(I)
            CORRY = AYI - AY(I)
            CORRZ = AZI - AZ(I)

            RX(I) = RX(I) + CR * CORRX
            RY(I) = RY(I) + CR * CORRY
            RZ(I) = RZ(I) + CR * CORRZ

            VX(I) = VX(I) + CV * CORRX
            VY(I) = VY(I) + CV * CORRY
            VZ(I) = VZ(I) + CV * CORRZ

            BX(I) = BX(I) + CB * CORRX
            BY(I) = BY(I) + CB * CORRY
            BZ(I) = BZ(I) + CB * CORRZ

            CX(I) = CX(I) + CC * CORRX
            CY(I) = CY(I) + CC * CORRY
            CZ(I) = CZ(I) + CC * CORRZ

            AX(I) = AXI
            AY(I) = AYI
            AZ(I) = AZI
        END DO
    END SUBROUTINE
END MODULE

PROGRAM MAIN
    USE GEAR
    USE LENNARD_JONES
    USE SIMULATION_BOX
    USE CONSTANTS

    IMPLICIT NONE

    INTEGER   :: CNUNIT               ! Unite de fichier (FORTRAN Specific)
    CHARACTER :: CNFILE * 30
    PARAMETER ( CNUNIT = 10 )

    INTEGER :: N_STEP, I, J, WRITEFREQ, FSTAT

    REAL(KIND = dp) :: K, DT, TC, TK
    REAL(KIND = dp) :: SF, SGF

    WRITE (*, '(''MOLECULAR RIJYNAMICS CONSTANT NVE'')')
    WRITE (*, '(''PARTICLE COUNT '')', ADVANCE = 'NO')
    READ  (*, *) N
    WRITE (*, '(''STEP COUNT '')', ADVANCE = 'NO')
    READ  (*, *) N_STEP
    WRITE (*, '(''TIME STEP '')', ADVANCE = 'NO')
    READ  (*, *) DT
    WRITE (*, '(''CUTOFF DISTANCE (ANGSTROM) '')', ADVANCE = 'NO')
    READ  (*, *) RCUT
    WRITE (*, '(''OUTPUT FREQUENCY '')', ADVANCE = 'NO')
    READ  (*, *) WRITEFREQ

    ! Initializes everything
    RCUT = RCUT * 1E-10_dp / UNIT_LENGTH
    DENSITY = LIQUID_DENSITY * (SIGMA ** 3) / M
    BOX = (REAL(N) / DENSITY) ** (1.0_dp / 3.0_dp)

    CALL INIT_BOX()
    CALL INIT_GEAR()
    CALL GENERATE_INIT_CONF()
    CALL INIT_LJ()

    VOLUME = BOX ** 3

    DT = DT * 1E-15_dp / UNIT_TIME

    WRITE (*, '(''******************* INITIAL PARAMETERS *******************'')')
    WRITE (*, '(''NUMBER OF ATOMS            = '', I5    )') N
    WRITE (*, '(''NUMBER OF STEPS            = '', I10   )') N_STEP
    WRITE (*, '(''DENSITY                    = '', F10.5 )') DENSITY * 1E-30_dp / (SIGMA ** 3)
    WRITE (*, '(''TIME STEP                  = '', F10.5 )') DT
    WRITE (*, '(''BOX SIZE                   = '', 2F10.5)') BOX, BOX * SIGMA / 1E-10_dp
    WRITE (*, '(''VOLUME                     = '', 2F10.5)') VOLUME, VOLUME * SIGMA ** 3 / 1E-30_dp
    WRITE (*, '(''CUTOFF RADIUS              = '', F10.5 )') RCUT * SIGMA / 1E-10_dp
    WRITE (*, '(''**********************************************************'')')
    WRITE (*, *)
    WRITE (*, '(''INITIAL CONFIGURATION'')')
    DO I = 1, N
        WRITE (*, *) RX(I), RY(I), RZ(I)
    END DO

    CALL PRINT_PAIR_DISTANCES()

    WRITE (*, *)

    TC = 0.0_dp
    TK = 0.0_dp
    SF = 0.0_dp
    SGF = 0.0_dp

    OPEN(UNIT = 50, FILE = 'output.txt', IOSTAT = FSTAT, STATUS = 'old')
    IF (FSTAT .EQ. 0) THEN
        CLOSE(UNIT = 50, STATUS = 'delete')
    END IF

    OPEN(UNIT = 50, FILE = 'output.txt', STATUS = 'new', ACTION = 'write')

    DO J = 1, N_STEP
        CALL PREDICT(DT)
        CALL CALCULATE_LJ()
        CALL CORRECT(DT)
        CALL BOUNDARY_CONDITIONS()
        CALL CORRECT_LJ()

        K = 0.0
        DO I = 1, N
            K = K + VX(I) ** 2 + VY(I) ** 2 + VZ(I) ** 2

            SGF = SGF + DFX(I) + DFY(I) + DFZ(I)
            SF = SF + (FX(I) ** 2 + FY(I) ** 2 + FZ(I) ** 2)
        END DO

        K = 1.0_dp / 2.0_dp * K

        TC = (1.0_dp / J) * (- SF / SGF)
        TK = (2.0_dp / 3.0_dp) * K / N

        IF (MOD(J, WRITEFREQ) .EQ. 0) THEN
            WRITE (*, *) J, V, K, (V + K), TC, TK
            WRITE (50, *) J, V, K, (V + K), TC, TK
        END IF
    END DO

    WRITE (*, *)
    WRITE (*, *) "Final configuration results"
    WRITE (*, *) "Energy", V, K, VC, V + K
    WRITE (*, *) "Temperature", TC, TK
    DO I = 1, N
        WRITE (*, *) "P", I, RX(I), RY(I), RZ(I)
        WRITE (*, *) "V", I, VX(I), VY(I), VZ(I)
        WRITE (*, *) "F", I, FX(I), FY(I), FZ(I)
        WRITE (*, *) "A", I, AX(I), AY(I), AZ(I)
        WRITE (*, *) "B", I, BX(I), BY(I), BZ(I)
        WRITE (*, *) "C", I, CX(I), CY(I), CZ(I)
    END DO

    WRITE (*, *) "Press any key to exit..."
    READ (*,*)

    CLOSE(UNIT = 50)

    CALL TERMINATE_BOX()
    CALL TERMINATE_GEAR()
    CALL TERMINATE_LJ()
END

SUBROUTINE GENERATE_INIT_CONF()
    USE SIMULATION_BOX
    USE GEAR
    USE LENNARD_JONES
    USE CONSTANTS
    IMPLICIT NONE

    INTEGER :: I, J
    REAL(KIND = dp) :: D, SCLX, SCLY, SCLZ

    CALL RANDOM_NUMBER(RX)
    CALL RANDOM_NUMBER(RY)
    CALL RANDOM_NUMBER(RZ)

    DO I = 1, N
        RX(I) = -BOX + RX(I) * BOX / 2.0_dp
        RY(I) = -BOX + RY(I) * BOX / 2.0_dp
        RZ(I) = -BOX + RZ(I) * BOX / 2.0_dp

        VX(I) = 0.0_dp
        VY(I) = 0.0_dp
        VZ(I) = 0.0_dp
        AX(I) = 0.0_dp
        AY(I) = 0.0_dp
        AZ(I) = 0.0_dp
        BX(I) = 0.0_dp
        BY(I) = 0.0_dp
        BZ(I) = 0.0_dp
        CX(I) = 0.0_dp
        CY(I) = 0.0_dp
        CZ(I) = 0.0_dp
        FX(I) = 0.0_dp
        FY(I) = 0.0_dp
        FZ(I) = 0.0_dp
    END DO

    DO WHILE (NEEDS_CORRECTION(1.0_dp))
        DO I = 1, N - 1
            DO J = I + 1, N
                D = (RX(I) - RX(J)) ** 2 + (RY(I) - RY(J)) ** 2 + (RZ(I) - RZ(J)) ** 2
                IF (D .LT. 1.0_dp) THEN
                    SCLX = 1.05_dp / ABS(RX(I) - RX(J))
                    SCLY = 1.05_dp / ABS(RY(I) - RY(J))
                    SCLZ = 1.05_dp / ABS(RZ(I) - RZ(J))

                    RX(J) = RX(J) * SCLX
                    RY(J) = RY(J) * SCLY
                    RZ(J) = RZ(J) * SCLZ
                END IF
            END DO
        END DO
        CALL BOUNDARY_CONDITIONS()
    END DO

END SUBROUTINE