# debug
debug = false

# variables
rCut = 8.0
nMol = 3
dt = 10^(-12)
nTime = 10^6

# L-J params
ε = 669.44
σ = 3.205*10^(-9)

# molecular masses
mass = 2.65*10^(-26)*ones(nMol, 1)

# array of molecules' evoulution
r  = zeros(3, nTime + 1, nMol)

# let's define the system's coordinate
# 1st molecule is at (0,0,0)
# 2nd molecule is at (1, 1, 1)
r[1, 1, 2] = 1
r[2, 1, 2] = 1
r[3, 1, 2] = 1
# 3nd molecule is at (-1, -1, 1)
r[1, 1, 3] = -1
r[2, 1, 3] = -1
r[3, 1, 3] = 1

rv = zeros(3, nTime + 1, nMol)
# let's define the system initial velocities
# 1st molecule moving 1754 m/s along the OX axis
rv[1, 1, 1] = 478
# 2nd molecule moving 1754 m/s along the OZ axis
rv[3, 1, 2] = -478
# 3rd molecule moving at 30 degrees slope towards the (0,0,0)
rv[1, 1, 3] = 413.96
rv[3, 1, 3] = -239

rf = zeros(3, nTime + 1, nMol)

# computing forces
function computeForce(x, y, z, fx, fy, fz)
    for i = 1:nMol-1
        for j = i+1:nMol
            dx = abs(x[i] - x[j])
            dy = abs(y[i] - y[j])
            dz = abs(z[i] - z[j])

            r = sqrt(dx^2 + dy^2 + dz^2)

            if debug
                println("Distance between atom no. ", i, " and no. ", j, " is equal to ", r, ".")
            end

            if r < rCut
                fr2 = σ^2 / r^2
                fr6 = fr2^3
                fpr = 48.0 * ε * fr6 * (fr6 - 0.5) / r^2

            if debug
                println("fpr of atoms ", i, " and ", j, " is equal to ", fpr)
            end

                fxi = fpr * dx
                fyi = fpr * dy
                fzi = fpr * dz

                if debug
                    println("Component fxi of ", i, " and ", j, " is equal to ", fxi, ".")
                    println("Component fyi of ", i, " and ", j, " is equal to ", fyi, ".")
                    println("Component fzi of ", i, " and ", j, " is equal to ", fzi, ".")
                end

                fx[i] = fx[i] + fxi
                fx[j] = fx[j] - fxi

                if debug
                    println("Component fx of ", i, " is equal to ", fx[i], ".")
                    println("Component fx of ", j, " is equal to ", fx[j], ".")
                end

                fy[i] = fy[i] + fyi
                fy[j] = fy[j] - fyi

                fz[i] = fz[i] + fzi
                fz[j] = fz[j] - fzi
            end
        end
    end

    if debug
        println("FORCES\nX: ", fx, "\nY: ", fy, "\nZ: ", fz)
    end
end

# let's set initial values of forces according to L-J potential
computeForce(r[1, 1, :], r[2, 1, :], r[3, 1, :], rf[1, 1, :], rf[2, 1, :], rf[3, 1, :])

# computing force
function Verlet(posVec, velVec, forVec, massVec)
    xPos = 1
    yPos = 2
    zPos = 3

    for timeStep = 1:nTime
        if debug
            println("Timestep no. ", timeStep, ".")
        end

        for iAtm = 1:nMol
            posVec[xPos, timeStep + 1, iAtm] = posVec[xPos, timeStep, iAtm] + velVec[xPos, timeStep, iAtm] * dt + forVec[xPos, timeStep, iAtm] / (2 * massVec[iAtm]) * dt^2
            posVec[yPos, timeStep + 1, iAtm] = posVec[yPos, timeStep, iAtm] + velVec[yPos, timeStep, iAtm] * dt + forVec[yPos, timeStep, iAtm] / (2 * massVec[iAtm]) * dt^2
            posVec[zPos, timeStep + 1, iAtm] = posVec[zPos, timeStep, iAtm] + velVec[zPos, timeStep, iAtm] * dt + forVec[zPos, timeStep, iAtm] / (2 * massVec[iAtm]) * dt^2
        end

        for i = 1:nMol-1
            Threads.@threads for j = i+1:nMol
                dx = abs(posVec[xPos,  timeStep + 1, i] - posVec[xPos,  timeStep + 1, j])
                dy = abs(posVec[yPos,  timeStep + 1, i] - posVec[yPos,  timeStep + 1, j])
                dz = abs(posVec[zPos,  timeStep + 1, i] - posVec[zPos,  timeStep + 1, j])

                r = sqrt(dx^2 + dy^2 + dz^2)

                if debug
                    println("Distance between atom no. ", i, " and no. ", j, " is equal to ", r, ".")
                end

                if r < rCut
                    fr2 = σ^2 / r^2
                    fr6 = fr2^3
                    fpr = 48.0 * ε * fr6 * (fr6 - 0.5) / r^2

                    if debug
                        println("fpr of atoms ", i, " and ", j, " is equal to ", fpr)
                    end

                    fxi = fpr * dx
                    fyi = fpr * dy
                    fzi = fpr * dz

                    if debug
                        println("Component fxi of ", i, " and ", j, " is equal to ", fxi, ".")
                        println("Component fyi of ", i, " and ", j, " is equal to ", fyi, ".")
                        println("Component fzi of ",

                         i, " and ", j, " is equal to ", fzi, ".")
                    end

                    forVec[xPos, timeStep, i] = forVec[xPos, timeStep, i] + fxi
                    forVec[xPos, timeStep, j] = forVec[xPos, timeStep, j] - fxi

                    if debug
                        println("Component fx of ", i, " is equal to ", forVec[xPos, timeStep, i], ".")
                        println("Component fx of ", j, " is equal to ", forVec[xPos, timeStep, j], ".")
                    end

                    forVec[yPos, timeStep, i] = forVec[yPos, timeStep, i] + fyi
                    forVec[yPos, timeStep, j] = forVec[yPos, timeStep, j] - fyi

                    forVec[zPos, timeStep, i] = forVec[zPos, timeStep, i] + fzi
                    forVec[zPos, timeStep, j] = forVec[zPos, timeStep, j] - fzi
                end
            end
        end
        #computeForce(posVec[xPos, trv\imeStep + 1, :], posVec[yPos, timeStep + 1, :], posVec[zPos, timeStep + 1, :], forVec[xPos, timeStep + 1, :], forVec[yPos, timeStep + 1, :], forVec[zPos, timeStep + 1, :])

        for iAtm = 1:nMol
            velVec[xPos, timeStep + 1, iAtm] = velVec[xPos, timeStep, iAtm] + (forVec[xPos, timeStep + 1, iAtm] + forVec[xPos, timeStep, iAtm]) / (2 * massVec[iAtm]) * dt
            velVec[yPos, timeStep + 1, iAtm] = velVec[yPos, timeStep, iAtm] + (forVec[yPos, timeStep + 1, iAtm] + forVec[yPos, timeStep, iAtm]) / (2 * massVec[iAtm]) * dt
            velVec[zPos, timeStep + 1, iAtm] = velVec[zPos, timeStep, iAtm] + (forVec[zPos, timeStep + 1, iAtm] + forVec[zPos, timeStep, iAtm]) / (2 * massVec[iAtm]) * dt
        end

        if debug
            println("After step no. ", timeStep + 1, " forces are: 1st atom:\n", forVec[:, timeStep + 1, 1], "\n2nd atom:\n", forVec[:, timeStep + 1, 2], "\n3rd atom:\n", forVec[:, timeStep + 1, 3])
        end
    end
end

@time Verlet(r, rv, rf, mass)