Merge pull request #5 from comecattin/read_file

Read input file and PBC conditions

.gitignore

@@ -31,3 +31,11 @@
 *.out
 *.app
 bin/
+
+# Results
+*.png
+*.gif
+*.mp4
+trajectories.dat
+input.dat
+md_simulation

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.

src/example_input.dat

@@ -0,0 +1,4 @@
+50       ! Number of particles
+10000    ! Number of steps
+0.001    ! Time step
+10.0     ! Box length

src/forces.f90

@@ -5,28 +5,32 @@ 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
 
         ! 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
 

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
 

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

src/main.f90

@@ -6,36 +6,74 @@ program md_simulation
     implicit none
 
     ! Define parameters
-    integer, parameter :: n_atoms = 100
-    integer, parameter :: n_steps = 1000
-    real(8), parameter :: dt = 0.001
-    character(len=100) :: output_file
+    integer :: n_atoms, n_steps, step, num_args
+    real(8) :: dt, box_length
+    character(len=100) :: input_file, 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(:,:)
+
+    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
+    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)
+    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
+        print *, 'Step: ', step
+
         ! 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)
-
+        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'
+
+    ! Deallocate the arrays
+    deallocate(positions, velocities, forces)
 
 
 

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
 

src/plot_trajectories.py

@@ -25,28 +25,35 @@ 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(0, 10)
-    ax.set_ylim(0, 10)
+    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)
 
     # 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()