From 251f633fc6b4b0b98869085afd8257326b4bf03c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 13:58:59 +0200
Subject: [PATCH 01/13] Initialize within a periodic box.

---
 src/initialization.f90 | 16 ++++++++--------
 src/main.f90           |  3 ++-
 2 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/initialization.f90 b/src/initialization.f90
index a5dee9d..82d28a7 100644
--- a/src/initialization.f90
+++ b/src/initialization.f90
@@ -2,19 +2,19 @@ module initialization_module
     implicit none
     
 contains
-    subroutine initialize(pos, vel, n)
+    subroutine initialize(pos, vel, n, box_length)
 
         real(8), dimension(n,3), intent(out) :: pos, vel
         integer, intent(in) :: n
+        real(8), intent(in) :: box_length
         integer :: i, j
 
-        ! Initialize the positions and velocities between -10 and 10
-        do i = 1, n
-            do j = 1, 3
-                pos(i,j) = 20.0 * (rand() - 0.5)
-                vel(i,j) = 20.0 * (rand() - 0.5)
-            end do
-        end do
+        call random_seed() ! Initialize the random number generator
+        call random_number(pos) ! Initialize the positions between 0 and 1
+        pos = pos * box_length
+
+        call random_number(vel)
+        vel = vel - 0.5 ! Shift the velocities to be between -0.5 and 0.5
 
     end subroutine initialize
     
diff --git a/src/main.f90 b/src/main.f90
index 120bcf1..572271a 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -9,6 +9,7 @@ program md_simulation
     integer, parameter :: n_atoms = 100
     integer, parameter :: n_steps = 1000
     real(8), parameter :: dt = 0.001
+    real(8), parameter :: box_length = 10.0
     character(len=100) :: output_file
 
     ! Define position, velocity and forces arrays
@@ -18,7 +19,7 @@ program md_simulation
     integer :: step
 
     ! Initialize the positions and velocities
-    call initialize(positions, velocities, n_atoms)
+    call initialize(positions, velocities, n_atoms, box_length)
 
     ! Open the output file
     output_file = 'trajectories.dat'

From a3382b281479265fe05e8bfeca35173e3212c63a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 14:07:28 +0200
Subject: [PATCH 02/13] Forces using PBC

---
 src/forces.f90 | 12 ++++++++----
 src/main.f90   |  2 +-
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/src/forces.f90 b/src/forces.f90
index 475b458..531149f 100644
--- a/src/forces.f90
+++ b/src/forces.f90
@@ -5,12 +5,14 @@ module forces_module
     
 contains
 
-    subroutine compute_forces(pos, frc, n)
+    subroutine compute_forces(pos, frc, n, box_length)
         real(8), dimension(n, 3), intent(in) :: pos
         real(8), dimension(n, 3), intent(out) :: frc
         integer, intent(in) :: n
+        real(8), intent(in) :: box_length
         integer :: i, j
         real(8) :: r, r2, r6, lj_force
+        real(8), dimension(3) :: dist
 
         ! Initialize forces to zero
         frc = 0.0
@@ -18,15 +20,17 @@ subroutine compute_forces(pos, frc, n)
         ! Compute Leenard-Jones forces
         do i = 1, n-1
             do j = i+1, n
-                r = sqrt(sum((pos(i, :) - pos(j, :))**2))
+                dist = pos(i, :) - pos(j, :)
+                dist = dist - nint(dist / box_length) * box_length
+                r = sqrt(sum(dist**2))
                 r2 = r**2
                 r6 = r2**3
 
                 lj_force = 24.0d0 * epsilon * (2.0d0 * sigma**6 / r6**2 - sigma**12 / r6)
 
                 ! Add forces to both atoms
-                frc(i, :) = frc(i, :) + lj_force * (pos(i, :) - pos(j, :)) / r
-                frc(j, :) = frc(j, :) - lj_force * (pos(i, :) - pos(j, :)) / r
+                frc(i, :) = frc(i, :) + lj_force * dist / r
+                frc(j, :) = frc(j, :) - lj_force * dist / r
             end do
         end do
 
diff --git a/src/main.f90 b/src/main.f90
index 572271a..8fbc988 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -27,7 +27,7 @@ program md_simulation
     ! Perform the molecular dynamics simulation
     do step = 1, n_steps
         ! Compute the forces
-        call compute_forces(positions, forces, n_atoms)
+        call compute_forces(positions, forces, n_atoms, box_length)
 
         ! Integrate positions and velocities
         call integrate(positions, velocities, forces, dt, n_atoms)

From d45b8ef8a88163f42443c7b8155ded96baa4e264 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 14:19:20 +0200
Subject: [PATCH 03/13] PBC for integration.

---
 src/integration.f90 | 21 +++++++++++++++------
 src/main.f90        |  7 +++++--
 2 files changed, 20 insertions(+), 8 deletions(-)

diff --git a/src/integration.f90 b/src/integration.f90
index f6949c9..25853f7 100644
--- a/src/integration.f90
+++ b/src/integration.f90
@@ -2,17 +2,26 @@ module integration_module
     implicit none
     
 contains
-    subroutine integrate(pos, vel, frc, dt, n)
+    subroutine integrate(pos, vel, frc, dt, n, box_length)
         real(8), dimension(n, 3), intent(inout) :: pos, vel
         real(8), dimension(n, 3), intent(in) :: frc
-        real(8), intent(in) :: dt
+        real(8), intent(in) :: dt, box_length
         integer, intent(in) :: n
-        integer :: i
+        integer :: i, j
 
-        ! Integrate using velocity verlet
         do i = 1, n
-            pos(i, :) = pos(i, :) + vel(i, :) * dt + 0.5 * frc(i, :) * dt**2
-            vel(i, :) = vel(i, :) + 0.5 * frc(i, :) * dt
+            do j = 1, 3
+                pos(i, j) = pos(i, j) + vel(i, j) * dt + 0.5 * frc(i, j) * dt**2
+                vel(i, j) = vel(i, j) + 0.5 * frc(i, j) * dt
+                
+                ! Apply periodic boundary conditions
+                if (pos(i, j) < 0.0d0) then
+                    pos(i, j) = pos(i, j) + box_length
+                else if (pos(i, j) >= box_length) then
+                    pos(i, j) = pos(i, j) - box_length
+                end if
+
+            end do
         end do
 
     end subroutine integrate
diff --git a/src/main.f90 b/src/main.f90
index 8fbc988..bb38617 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -26,14 +26,17 @@ program md_simulation
 
     ! Perform the molecular dynamics simulation
     do step = 1, n_steps
+        print *, 'Step: ', step
+
         ! Compute the forces
         call compute_forces(positions, forces, n_atoms, box_length)
 
         ! Integrate positions and velocities
-        call integrate(positions, velocities, forces, dt, n_atoms)
-
+        call integrate(positions, velocities, forces, dt, n_atoms, box_length)
+        
         ! Output the positions
         call output_positions(step, positions, n_atoms, output_file)
+
     end do
 
     print *, 'Simulation finished'

From b1355d67f62cf99579a9a70bba7b7162fc5dd037 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 14:23:52 +0200
Subject: [PATCH 04/13] Ignore output files

---
 .gitignore | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/.gitignore b/.gitignore
index 24184aa..1a00e82 100644
--- a/.gitignore
+++ b/.gitignore
@@ -31,3 +31,9 @@
 *.out
 *.app
 bin/
+
+# Results
+*.png
+*.gif
+*.mp4
+*.dat

From 6913debf19c8620c1fc25281e80e20757107d40d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 14:45:39 +0200
Subject: [PATCH 05/13] Clean output files

---
 src/makefile | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/makefile b/src/makefile
index 21ffbfe..6de57c0 100644
--- a/src/makefile
+++ b/src/makefile
@@ -20,5 +20,5 @@ output_module.o: output_module.f90
 	$(FC) -c output_module.f90
 
 clean:
-	rm -f *.o *.mod md_simulation
+	rm -f *.o *.mod md_simulation trajectories.dat
 

From a1b0c144461a1f66e9caffad3a1345be2289e278 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 14:46:10 +0200
Subject: [PATCH 06/13] Save trajectories in MP4 format.

---
 src/plot_trajectories.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/plot_trajectories.py b/src/plot_trajectories.py
index 5a4e4cc..32bfc46 100644
--- a/src/plot_trajectories.py
+++ b/src/plot_trajectories.py
@@ -28,8 +28,8 @@ def animate(i):
     ax.clear()
     pos = positions[i]
     ax.scatter(pos[:, 0], pos[:, 1], s=10)
-    ax.set_xlim(0, 10)
-    ax.set_ylim(0, 10)
+    ax.set_xlim(-1, 11)
+    ax.set_ylim(-1, 11)
     ax.set_title(f'Time Step: {steps[i]}')
     ax.set_xlabel('X')
     ax.set_ylabel('Y')
@@ -47,6 +47,7 @@ def main():
     # Save the animation
     writer = PillowWriter(fps=24)
     ani.save('particle_trajectories.gif', writer=writer)
+    ani.save('particle_trajectories.mp4', writer='ffmpeg', fps=24)
 
 if __name__ == '__main__':
     main()

From 2277e0006f73d65038eb9956b0e9cbbb2ea5c770 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 15:25:12 +0200
Subject: [PATCH 07/13] Set the x and y limits as dependent of the positions.

---
 src/plot_trajectories.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/plot_trajectories.py b/src/plot_trajectories.py
index 32bfc46..562dd46 100644
--- a/src/plot_trajectories.py
+++ b/src/plot_trajectories.py
@@ -25,22 +25,28 @@ def read_trajectories(filename):
     return steps, positions
 
 def animate(i):
+
     ax.clear()
     pos = positions[i]
     ax.scatter(pos[:, 0], pos[:, 1], s=10)
-    ax.set_xlim(-1, 11)
-    ax.set_ylim(-1, 11)
+    ax.set_xlim(0, x_max)
+    ax.set_ylim(0, y_max)
     ax.set_title(f'Time Step: {steps[i]}')
     ax.set_xlabel('X')
     ax.set_ylabel('Y')
 
 def main():
-    global steps, positions, ax
+    global steps, positions, ax, x_max, y_max
     steps, positions = read_trajectories('trajectories.dat')
 
     steps = steps[::20]
     positions = positions[::20]
 
+    # x and y limits depend on the positions
+    positions = np.array(positions)
+    x_max = np.max(positions[:, :, 0])
+    y_max = np.max(positions[:, :, 1])
+
     fig, ax = plt.subplots(figsize=(8, 8))
     ani = FuncAnimation(fig, animate, frames=len(steps), repeat=False)
 

From 170857f4b94d2c52a6e78ba44c41c2f6ee0a5326 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 15:25:37 +0200
Subject: [PATCH 08/13] Read a input.dat file for parameters

---
 src/main.f90 | 39 ++++++++++++++++++++++++++++++++-------
 1 file changed, 32 insertions(+), 7 deletions(-)

diff --git a/src/main.f90 b/src/main.f90
index bb38617..4d22e44 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -6,23 +6,45 @@ program md_simulation
     implicit none
 
     ! Define parameters
-    integer, parameter :: n_atoms = 100
-    integer, parameter :: n_steps = 1000
-    real(8), parameter :: dt = 0.001
-    real(8), parameter :: box_length = 10.0
+    integer :: n_atoms, n_steps, step
+    real(8) :: dt, box_length
     character(len=100) :: output_file
+    logical :: file_exists
 
     ! Define position, velocity and forces arrays
-    real(8), dimension(n_atoms, 3) :: positions, velocities, forces
+    real(8), allocatable:: positions(:,:), velocities(:,:), forces(:,:)
+
+    ! Read the input parameters
+    open(unit=20, file='input.dat')
+    read(20, *) n_atoms
+    read(20, *) n_steps
+    read(20, *) dt
+    read(20, *) box_length
+    close(20)
+
+    ! Print the input parameters
+    print *, 'Number of atoms: ', n_atoms
+    print *, 'Number of steps: ', n_steps
+    print *, 'Time step: ', dt
+    print *, 'Box length: ', box_length
+
+    ! Allocate the arrays
+    allocate(positions(3, n_atoms))
+    allocate(velocities(3, n_atoms))
+    allocate(forces(3, n_atoms))
 
-    ! Other variables
-    integer :: step
 
     ! Initialize the positions and velocities
     call initialize(positions, velocities, n_atoms, box_length)
 
     ! Open the output file
     output_file = 'trajectories.dat'
+    ! Check if the file already exists
+    inquire(file=output_file, exist=file_exists)
+    if (file_exists) then
+        print *, 'The file ', output_file, ' already exists. Removing it.'
+        call system('rm ' //trim(output_file))
+    end if
 
     ! Perform the molecular dynamics simulation
     do step = 1, n_steps
@@ -40,6 +62,9 @@ program md_simulation
     end do
 
     print *, 'Simulation finished'
+
+    ! Deallocate the arrays
+    deallocate(positions, velocities, forces)
     
 
     

From 16e2733cfabf00a72f92e670b0b0235ec0aceeda Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 15:35:49 +0200
Subject: [PATCH 09/13] Read input parameters from a file

---
 src/main.f90 | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/src/main.f90 b/src/main.f90
index 4d22e44..722aec5 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -6,14 +6,23 @@ program md_simulation
     implicit none
 
     ! Define parameters
-    integer :: n_atoms, n_steps, step
+    integer :: n_atoms, n_steps, step, num_args
     real(8) :: dt, box_length
-    character(len=100) :: output_file
+    character(len=100) :: input_file, output_file
     logical :: file_exists
 
     ! Define position, velocity and forces arrays
     real(8), allocatable:: positions(:,:), velocities(:,:), forces(:,:)
 
+    num_args = command_argument_count()
+    if (num_args /= 1) then
+        print *, 'Usage: ./md_simulation input_file'
+        stop
+    else
+        call get_command_argument(1, input_file)
+        print *, 'Input file: ', input_file
+    end if
+
     ! Read the input parameters
     open(unit=20, file='input.dat')
     read(20, *) n_atoms

From c7df45b6009f7108d14edb38ee607fb33eab2419 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 15:43:18 +0200
Subject: [PATCH 10/13] Updated gitignore

---
 .gitignore | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.gitignore b/.gitignore
index 1a00e82..312a991 100644
--- a/.gitignore
+++ b/.gitignore
@@ -36,4 +36,5 @@ bin/
 *.png
 *.gif
 *.mp4
-*.dat
+trajectories.dat
+input.dat

From 6e8eea810ec861976a328c2d7909932be7137969 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 15:43:36 +0200
Subject: [PATCH 11/13] Documentation for input file

---
 README.md | 27 +++++++++++++++++++++------
 1 file changed, 21 insertions(+), 6 deletions(-)

diff --git a/README.md b/README.md
index 610be02..245e44c 100644
--- a/README.md
+++ b/README.md
@@ -1,23 +1,38 @@
 # F-MD, Fortran Molecular Dynamic
+
 A simple Molecular Dynamic implementation in Fortran.
 
 ## Installation
+
 To install, clone this repository:
 
-```
+```bash
 git clone git@github.com:comecattin/F-MD.git
 cd F-MD
 ```
 
 Then compile the source:
-```
+
+```bash
 cd src
 make
 ```
+
 ## Run an MD simulation
-To run an MD simulation simply run `./md_simulation` in the `src` directory.
 
-Parameters for the MD simulation can be changed inside the `src/main.f90` file. Compilation is needed after making changes.
+To run an MD simulation simply run `./md_simulation input.dat` in the `src` directory.
+
+The file `input.dat` is the input file and contain in the following order:
+
+- The number of particles
+- The number of time steps
+- The time step
+- The box length
+
+An example is given in the `example_input.dat` file.
+
+## Output and plot the trajectories
+
+By default, positions along time steps are outputted in the `trajectories.dat` file.
 
-## Ouput and plot the trajectories
-To plot an animation of the trajectories, run `python plot_trajectories.py` in the `src` directory. 
+To plot an animation of the trajectories, run `python plot_trajectories.py` in the `src` directory.

From 83c0c8e7488e16f04184e083bc78029d94837aa2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 15:43:57 +0200
Subject: [PATCH 12/13] Example of an input file.

---
 src/example_input.dat | 4 ++++
 1 file changed, 4 insertions(+)
 create mode 100644 src/example_input.dat

diff --git a/src/example_input.dat b/src/example_input.dat
new file mode 100644
index 0000000..ca0742d
--- /dev/null
+++ b/src/example_input.dat
@@ -0,0 +1,4 @@
+50       ! Number of particles
+10000    ! Number of steps
+0.001    ! Time step
+10.0     ! Box length

From 2676cda52ff218fc52d541383ab90d9ebd8ba5b4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?C=C3=B4me=20Cattin?= <come.cattin@ens-lyon.fr>
Date: Mon, 27 May 2024 15:44:38 +0200
Subject: [PATCH 13/13] Ignore the main executable.

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index 312a991..65321fa 100644
--- a/.gitignore
+++ b/.gitignore
@@ -38,3 +38,4 @@ bin/
 *.mp4
 trajectories.dat
 input.dat
+md_simulation