Fixed sign error in Lennard-Jones force and slight modifications to the class verson

2 files changed, 21 insertions, 18 deletions

diff --git a/class/class-langevin-motion.f90 b/class/class-langevin-motion.f90
index b52396c..0b97d2d 100644
--- a/class/class-langevin-motion.f90
+++ b/class/class-langevin-motion.f90
@@ -39,6 +39,7 @@ function get_neighbour_ids(p, N) result (neighbours)
     integer :: max_list(8)=-1
 
     call ONETO2(p,N,x_cell,y_cell)
+    max_list(
 
     !> Start by just getting all of the neighbours
     do i = -1, 1
diff --git a/new/self-langevin-motion.f90 b/new/self-langevin-motion.f90
index 6c048d9..92ac84d 100644
--- a/new/self-langevin-motion.f90
+++ b/new/self-langevin-motion.f90
@@ -11,27 +11,30 @@ program langevin_motion
     real(real64), parameter :: L = 1.0_real64
 
     !> Numerical constants and such
-    real(real64), parameter :: dt = 0.005_real64
-    real(real64), parameter :: t_max = 10.0_real64
+    real(real64), parameter :: dt = 0.00001_real64
+    real(real64), parameter :: t_max = 1.0_real64
     real(real64), parameter :: kT = 10.0_real64
     real(real64), parameter :: g = 1.0_real64
     real(real64), parameter :: m = 1.0_real64
     real(real64), parameter :: eps = 1.0_real64
     real(real64), parameter :: sigma = 1E-03_real64
     real(real64), parameter :: rc = sigma*2**(1.0_real64/6.0_real64)
-    real(real64), parameter :: pref1 = g, pref2 = sqrt(24.0_real64*kT*g/dt)
+    real(real64), parameter :: pref1 = g, pref3 = sqrt(24.0_real64*kT*g/dt)
 
     !> And non-constants
-    real(real64) :: t
+    real(real64) :: t = 0.0_real64
     real(real64), dimension(n) :: dx, dy, dn
     real(real64), dimension(n) :: F
-    real(real64) :: ran1, ran2
+    real(real64), dimension(n) :: ran1, ran2
     integer(int64) :: i
 
     !> Vectors for positon, velocity, acceleration, and tracking 
     real(real64), dimension(n) :: px, py, vx, vy, ax, ay, vhx, vhy
     logical, dimension(n) :: is_tracked = .TRUE., are_interacting
 
+    !> Prepare random numbers
+    call random_seed()
+
     !> Start by initializing the particles with a random pos. and vel.
     call initialize_particles(px, py, vx, vy, ax, ay, kT, L, pi)
 
@@ -53,21 +56,22 @@ program langevin_motion
         ax = 0.0_real64
         ay = 0.0_real64
 
+        !> Random numbers
         call random_number(ran1)
         ran1 = ran1 - 0.5_real64
         call random_number(ran2)
         ran2 = ran2 - 0.5_real64
 
         !> Random force
-        ax = ax - pref1*vhx + ran1
-        ay = ay - pref1*vhy + ran2
+        !ax = ax - pref1*vhx + pref2*ran1
+        !ay = ay - pref1*vhy + pref2*ran2
 
         !> Particle-particle interaction force
         interactions: do i = 1,n
             
             !> Calculate the distance between the particles:
-            dx = px - px(i)
-            dy = py - py(i)
+            dx = px(i) - px
+            dy = py(i) - py
             dn = sqrt(dx**2 + dy**2)
 
             !> Check which particles are interacting
@@ -77,10 +81,8 @@ program langevin_motion
             F = 0.0_real64
 
             where (are_interacting)
-
-                !> Lennardn-Jones force
-                F = 4.0_real64*eps*(-12.0_real64*sigma**12/dn**13 + 6.0_real64*sigma**6/dn**7)
-
+                !> Lennard-Jones force
+                F = 4.0_real64*eps*(12.0_real64*sigma**12/dn**13 - 6.0_real64*sigma**6/dn**7)
             end where
 
             !> Update the accelerations with the new force when a force was calculated
@@ -97,11 +99,13 @@ program langevin_motion
         t = t + dt
 
         !> Print posiitons to stdout
-        print '(e24.12,e24.12)', (px(i), py(i), i=1,n)
+        !print '(e24.12,e24.12)', (px(i), py(i), i=1,n)
+
+        print *, sum(0.5_real64*m*vx**2 + 0.5_real64*m*vy**2, mask=is_tracked)
 
         !> This splits timesteps into data blocks for Gnuplot
-        print *
-        print *
+        !print *
+        !print *
 
     end do
 
@@ -113,8 +117,6 @@ contains
         real(real64) :: ran1, ran2
 
         !> Random number seed
-        call random_seed()
-
         !> And 2 random numbers
         call random_number(ran1)
         call random_number(ran2)