Merge pull request #18 from comecattin/yukawa

New water model

README.md

@@ -6,7 +6,38 @@
 
 A simple Molecular Dynamic implementation in Fortran.
 
-Molecules are randomly initialized in a simulation box with Periodic Bondary Condition using minimal image. Lennard-Jones force are then propagated along the time using a Velocity Verlet integrator.
+Water molecules are randomly initialized in a simulation box with Periodic Bondary Condition using minimal image. This work relies on the SPC/Fw water model, where interaction potential is described as:
+
+$$V = V^\text{intra} + V^\text{inter}$$
+
+with
+$$
+  V^\text{intra} = \cfrac{k_b}{2}\left[
+    (r_{\text{OH}_1} - r_\text{OH}^\text{eq})^2
+    + (r_{\text{OH}_2} - r_\text{OH}^\text{eq})^2
+    \right]
+    + \cfrac{k_a}{2} (\theta_\text{HOH} - \theta_\text{HOH}^\text{eq})^2
+$$
+
+and
+$$
+  V^\text{inter} = \sum_{ij}^\text{all pairs} \left\{
+    4\epsilon_{ij} \left[
+      \left(\cfrac{\sigma_{ij}}{R_{ij}} \right)^{12} -
+      \left(\cfrac{\sigma_{ij}}{R_{ij}} \right)^{6}
+      \right] -
+      q_i q_j \cfrac{e^{-\alpha m R_{ij}}}{R_{ij}}
+  \right\}
+$$
+
+Parameters are given in the following table:
+|$k_b$|$r_\text{OH}^\text{eq}$ (Å) | $k_a$ | $\theta_\text{HOH}^\text{eq}$ (deg) | $\sigma_\text{OO}$ (Å) | $\epsilon_\text{OO}$ (kcal.mol-1) | $q(\text{O})$ (e) | $q(\text{H})$ (e) |
+|---------|---------|--------|---------|--------|--------|--------|--------|
+| 1059.162 | 1.012 | 75.90 | 113.24 | 3.165492 | 0.1554253 | -0.82 | 0.41 |
+
+- Bonded and angle interactions are modeled by harmonic potentials
+- van der Waals interactions are modeled by a Lennard-Jones potential
+- Coulomb interactions are modeled by a Yukawa potential, *ie.* screened Coulomb potential
 
 ## Installation
 
@@ -34,8 +65,7 @@ The file `input.in` is the input file and contain:
 - The number of time steps (`n_steps`, default `1000`)
 - The time step (`dt`, default `0.001`)
 - The box length (`box_length`, default `10.0`)
-- The tolerance of the SHAKE algorithm for angle and bonds length constraints (`tolerance`, default `1e-6`)
-- The maximum number of iteration for the SHAKE algorithm (`max_iter`, default `100`)
+- Stride, write position, energies and print them every nth step (`stride`, default `1`)
 
 An example is given in the `examples/example_input.in` file.
 

examples/example_input.in

@@ -1,6 +1,5 @@
-n_atoms : 30          # Number of atoms in the system
-n_steps : 10000       # Number of time steps
-dt : 0.001            # Time step
-box_length : 10.0     # Length of the PBC box
-tolerance : 1.0e-6    # SHAKE tolerance
-max_iter : 100        # Maximum iteration for SHAKE algorithm
+n_atoms : 60           # Number of atoms in the system
+n_steps : 100000       # Number of time steps
+dt : 1e-5              # Time step
+box_length : 10.0      # Length of the PBC box
+stride : 100           # Output information every nth step

src/energies.f90

@@ -1,71 +1,136 @@
 module energies_module
     implicit none
-    real(8), parameter :: epsilon = 1.0 ! LJ potential well depth
-    real(8), parameter :: sigma = 1.0 ! LJ distance where potential is zero
+    real(8), parameter :: pi = 3.14159265358979323846
+    real(8), parameter :: epsilon = 0.1554253          ! kcal/mol
+    real(8), parameter :: sigma = 3.165492             ! Angstrom 
+    real(8), parameter :: ka_HOH = 75.90               ! kcal.mol-1.rad-2
+    real(8), parameter :: kb_OH = 1059.162             ! kcal.mol-1.angstrom-2
+    real(8), parameter :: alpha = 1.0
+    real(8), parameter :: r_eq_OH = 1.012              ! Angstrom
+    real(8), parameter :: angle_eq = 113.24            ! deg
 
 contains 
-    subroutine compute_energies(pos, vel, charges, n, ke, pe, te)
+    subroutine compute_energies(pos, vel, mass, charges, n, box_length, ke, pe, te)
         real(8), dimension(n, 3), intent(in) :: pos, vel
-        real(8), dimension(n), intent(in) :: charges
+        real(8), dimension(n), intent(in) :: charges, mass
         integer, intent(in) :: n
+        real(8), intent(in) :: box_length
         real(8), intent(out) :: ke, pe, te
 
-        call compute_kinetic_energy(vel, n, ke)
-        ! call compute_potential_energy(pos, n, pe)
-        call compute_potential_energy_water(pos, n, charges, pe)
+        call compute_kinetic_energy(vel, mass, n, ke)
+        call compute_potential_energy(pos, charges, n, box_length, pe)
         te = ke + pe
 
-        print *, "Kinetic energy: ", ke
-        print *, "Potential energy: ", pe
-        print *, "Total energy: ", te
     end subroutine compute_energies
 
-    subroutine compute_kinetic_energy(vel, n, ke)
+    subroutine compute_kinetic_energy(vel, mass, n, ke)
         real(8), dimension(n, 3), intent(in) :: vel
         integer, intent(in) :: n
+        real(8), dimension(n), intent(in) :: mass
         real(8), intent(out) :: ke
         integer :: i
 
         ke = 0.0
         do i = 1, n
-            ke = ke + 0.5 * sum(vel(i, :) ** 2)
+            ke = ke + 0.5 * mass(i) * sum(vel(i, :) ** 2)
         end do
     end subroutine compute_kinetic_energy
 
-    subroutine compute_potential_energy(pos, n, pe)
+    subroutine compute_potential_energy(pos, charges, n, box_length, pe)
         real(8), dimension(n, 3), intent(in) :: pos
+        real(8), dimension(n), intent(in) :: charges
+        real(8), intent(in) :: box_length
         integer, intent(in) :: n
         real(8), intent(out) :: pe
+        real(8) :: harmonic_p_bond, harmonic_p_angle, lennard_jones_p, yukawa_p
         integer :: i, j
         real(8) :: r
+        real(8), dimension(3) :: vec1, vec2, dist
+        real(8) :: theta, cos_theta, angle_eq_rad
+
+        angle_eq_rad = angle_eq * pi / 180.0
 
         pe = 0.0
-        do i = 1, n
-            do j = i + 1, n
-                r = sqrt(sum((pos(i, :) - pos(j, :)) ** 2))
-                pe = pe + 4.0 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)
+
+        ! Intra molecular potential
+        do i = 1, n-1, 3
+            ! Harmonic potential on HOH angle
+            vec1 = pos(i+1, :) - pos(i, :)
+            vec2 = pos(i+2, :) - pos(i, :)
+            ! Minimum image convention
+            vec1 = vec1 - nint(vec1 / box_length) * box_length
+            vec2 = vec2 - nint(vec2 / box_length) * box_length
+            cos_theta = dot_product(vec1, vec2) / (norm2(vec1) * norm2(vec2))
+            theta = acos(cos_theta)
+            call harmonic_potential(theta, angle_eq_rad, ka_HOH, harmonic_p_angle)
+            pe = pe + harmonic_p_angle
+
+            ! Harmonic potential on O-H bond
+            !    O-H1 bond
+            call harmonic_potential(norm2(vec1), r_eq_OH, kb_OH, harmonic_p_bond)
+            pe = pe + harmonic_p_bond
+            !    O-H2 bond
+            call harmonic_potential(norm2(vec2), r_eq_OH, kb_OH, harmonic_p_bond)
+            pe = pe + harmonic_p_bond
+        end do
+
+        ! Lennard-Jones potential
+        do i = 1, n-4, 3
+            do j = i+3, n, 3
+                dist = pos(i, :) - pos(j, :)
+                ! Minimum image convention
+                dist = dist - nint(dist / box_length) * box_length
+                r = norm2(dist)
+                call lennard_jones_potential(r, lennard_jones_p)
+                pe = pe + lennard_jones_p
             end do
         end do
-    end subroutine compute_potential_energy
-
-    subroutine compute_potential_energy_water(pos, n, charges, pe)
-        real(8), dimension(n, 3), intent(in) :: pos
-        integer, intent(in) :: n
-        real(8), dimension(n), intent(in) :: charges
-        real(8), intent(out) :: pe
-        integer :: i, j
-        real(8) :: r
 
-        pe = 0.0
-        do i = 1, n
-            do j = i + 1, n
-                r = sqrt(sum((pos(i, :) - pos(j, :)) ** 2))
-                ! Lennard Jones potential
-                pe = pe + 4.0 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)
-                ! Coulomb potential
-                pe = pe + charges(i) * charges(j) / r
+        ! Coulomb-Yukawa potential
+        do i = 1, n-1
+            do j = i+1, n
+                dist = pos(i, :) - pos(j, :)
+                ! Minimum image convention
+                dist = dist - nint(dist / box_length) * box_length
+                r = norm2(dist)
+                call yukawa_potential(r, 1.0_8, charges(i), charges(j), alpha, yukawa_p)
+                pe = pe + yukawa_p
             end do
         end do
-    end subroutine compute_potential_energy_water
+
+
+
+    end subroutine compute_potential_energy
+
+    subroutine yukawa_potential(r, mass, q1, q2, alpha_param, yukawa_p)
+
+        real(8), intent(in) :: r
+        real(8), intent(in) :: q1, q2
+        real(8), intent(in) :: mass, alpha_param
+        real(8), intent(out) :: yukawa_p
+
+        yukawa_p = - q1*q2 * exp(-alpha_param * mass * r) / r
+
+    end subroutine yukawa_potential
+
+    subroutine lennard_jones_potential(r, lj_p)
+
+        real(8), intent(in) :: r
+        real(8), intent(out) :: lj_p
+
+        lj_p = 4 * epsilon * ((sigma / r)**12 - (sigma / r)**6)
+
+    end subroutine lennard_jones_potential
+
+    subroutine harmonic_potential(r, r_eq, kb, harmonic_p)
+
+        real(8), intent(in) :: r
+        real(8), intent(in) :: r_eq, kb
+        real(8), intent(out) :: harmonic_p
+
+        harmonic_p = 0.5 * kb * (r - r_eq)**2
+
+    end subroutine harmonic_potential
+
 
 end module energies_module