diff --git a/projects/project_4/poisson_2.py b/projects/project_4/poisson_2.py
index 23e14e5..f6d448f 100644
--- a/projects/project_4/poisson_2.py
+++ b/projects/project_4/poisson_2.py
@@ -1,6 +1,7 @@
 import numpy as np
 
 
+
 def create_laplacian_2d(nx, ny, lx, ly, pbc=True):
     """ Computes discrete Laplacian for a 2d
         charge density matrix, ordered row-wise
diff --git a/projects/project_4/test_poisson_2.py b/projects/project_4/test_poisson_2.py
index b616683..f66e826 100644
--- a/projects/project_4/test_poisson_2.py
+++ b/projects/project_4/test_poisson_2.py
@@ -37,10 +37,12 @@ def test_laplacian_2d():
     (3, 3, 1, -1, True, ValueError),
     (3, 3, -1, 1, True, ValueError),
     ])
+
 def test_laplacian_2d_exceptions(nx, ny, lx, ly, pbc, exception):
     with pytest.raises(exception):
         create_laplacian_2d(nx, ny, lx, ly, pbc=pbc)
 
+
 ### test TypeError ###
 list1 = [(i, 3, 1, 1, True, TypeError)
      for i in [[0], 'string', {'dict': 0}, True]]
@@ -52,20 +54,25 @@ def test_laplacian_2d_exceptions(nx, ny, lx, ly, pbc, exception):
     for i in [[0], 'string', {'dict': 0}, True]])
 list1.extend([(3, 3, 1, 1, i, TypeError)
     for i in [[0], 'string', {'dict': 0}]])
+
 @pytest.mark.parametrize('nx, ny, lx, ly,pbc,exception', list1)
 def test_laplacian_2d_exceptions(nx, ny, lx, ly, pbc, exception):
     with pytest.raises(exception):
         create_laplacian_2d(nx, ny, lx, ly, pbc=pbc)
 
+
 ### test periodc grid ###
 @pytest.mark.parametrize('nx, ny', [
     (100, 100),
+    (100, 90),
+    (90, 100),
     ])
+
 def test_laplacian_2d(nx, ny):
     x = np.linspace(-np.pi, np.pi, nx, endpoint=False)
     y = np.linspace(-np.pi, np.pi, ny, endpoint=False)
-    xx, yy = np.meshgrid(x,y, indexing='ij')
+    yy, xx = np.meshgrid(x,y, indexing='xy')
     grid = np.sin(xx) + np.cos(yy)
     laplacian = create_laplacian_2d(nx, ny, 2 * np.pi, 2 * np.pi)
     grid_lap = np.dot(laplacian, grid.reshape(-1))
-    np.testing.assert_almost_equal(grid, -grid_lap.reshape(nx, ny),decimal=3)
+    np.testing.assert_almost_equal(grid, -grid_lap.reshape(ny, nx),decimal=3)
diff --git a/projects/project_5/gauss_seidel.py b/projects/project_5/gauss_seidel.py
new file mode 100644
index 0000000..230dc89
--- /dev/null
+++ b/projects/project_5/gauss_seidel.py
@@ -0,0 +1,62 @@
+import numpy as np
+
+
+def compute_phi_2d(rho, l, max_it=1E4, f_tol=1E-7, epsilon_0=1):
+    """
+        Computes potential from charge distribution using
+        Gauss Seidel Method. Does not check for correct input.
+        Arguments:
+            rho (np.array): Charge distribution. axis=0 is y, axis=1 is x
+                Ensure that sum(rho) = 0
+            l (list): list containing lx, ly. min(l) > 0
+            f_tol(float): Error tolerance
+                Iteration is terminated if change of norm of potenital < f_tol
+            max_it(int): maximum number of iterations
+            epsilon_0(float): vacuum permittivity
+    """
+    ### define relevant variables ###
+    nx = rho.shape[1]
+    ny = rho.shape[0]
+    n_min = min([nx, ny])
+    hx = (nx / l[0]) ** 2
+    hy = (ny / l[1]) ** 2
+    hxhy = 1 / (2 * hx + 2 * hy)
+    phi = np.zeros((ny, nx))
+    norm = 1
+    n_it = 0
+
+    ### define indexing lists for loop ###
+    index_y  = [ny - i if ny - i > 0 else 0 for i in range(1, nx + ny)]
+    index_x  = [i - ny if i - ny > 0 else 0 for i in range(1, nx + ny)]
+    indices  = [(a,b) for a,b in zip(index_y,index_x)]
+    diag_len = [n_min for i in np.ones(nx + ny - 1, dtype=np.int8)]
+    for i in range(n_min):
+        diag_len[i] = i + 1
+        diag_len[-i-1] = i + 1
+
+    ### run main loop ###
+    while norm > f_tol and n_it < max_it:
+        n_it += 1
+        norm_old = np.linalg.norm(phi)
+        for k, l in enumerate(indices):
+            for m in range(diag_len[k]):
+                i = l[0] + m
+                j = l[1] + m
+                phi_up = phi[i-1,j]
+                phi_left = phi[i,j-1]
+                if i < ny - 1:
+                    phi_down = phi[i+1,j]
+                else:
+                    phi_down = phi[0,j]
+                if j < nx - 1:
+                    phi_right = phi[i,j+1]
+                else:
+                    phi_right = phi[i,0]
+                phi[i,j] = (hxhy * ( hx * (phi_right + phi_left)
+                                   + hy * (phi_up + phi_down)
+                                   + 1 / epsilon_0 * rho[i,j]
+                                   ))
+        norm_new = np.linalg.norm(phi)
+        norm = abs(norm_old - norm_new)
+
+    return phi
diff --git a/projects/project_5/poisson_scipy.py b/projects/project_5/poisson_scipy.py
new file mode 100644
index 0000000..a132c48
--- /dev/null
+++ b/projects/project_5/poisson_scipy.py
@@ -0,0 +1,46 @@
+import numpy as np
+from scipy.sparse import diags
+
+
+def create_laplacian_2d(nx, ny, lx, ly, pbc=True):
+    """ Computes discrete Laplacian for a 2d
+        charge density matrix, ordered row-wise
+        Args:
+            nx(int  >= 2): number of grid points along x axis
+            ny(int  >= 2): number of grid points along y axis
+            lx(float > 0): length of grid along x axis
+            ly(float > 0): length of grid along y axis
+            pbc(boolean): periodic boundry conditions
+        output:
+            Laplacian as nx * ny by nx * ny np.array
+    """
+    if type(nx) != int or type(ny) != int:
+        raise TypeError('We need an integer')
+    if type(lx) != int and type(lx) != float:
+        raise TypeError('We need a number')
+    if type(ly) != int and type(ly) != float:
+        raise TypeError('We need a number')
+    if nx < 2 or ny < 2:
+        raise ValueError('We need at least two grid points')
+    if lx <= 0 or ly <= 0:
+        raise ValueError('We need a positive length')
+    if type(pbc) != bool:
+        raise TypeError('We need a boolean as pbc')
+
+    hx = (nx / lx) ** 2
+    hy = (ny / ly) ** 2
+    diag0 = (-2 * hx - 2 * hy) * np.ones(nx * ny)
+    diag1 = hx * np.ones(nx * ny - 1)
+    diag1[nx-1::nx] = 0
+    diag2 = hy * np.ones(nx * ny - nx)
+    diagonals = [diag0, diag1, diag1, diag2, diag2]
+    offsets   = [0, 1, -1, nx, -nx]
+
+    if pbc:
+        diag4 = hy * np.ones(nx)
+        diag5 = hx * np.zeros(nx * ny - nx + 1)
+        diag5[::nx] = hx
+        diagonals.extend([diag4, diag4, diag5, diag5])
+        offsets.extend([nx * ny - nx, nx - nx * ny, nx - 1, 1 - nx])
+
+    return diags(diagonals, offsets).todense()
diff --git a/projects/project_5/test_gauss_seidel.py b/projects/project_5/test_gauss_seidel.py
new file mode 100644
index 0000000..69d72d9
--- /dev/null
+++ b/projects/project_5/test_gauss_seidel.py
@@ -0,0 +1,21 @@
+import numpy as np
+import pytest
+from poisson_scipy import create_laplacian_2d
+from gauss_seidel import compute_phi_2d
+
+
+### test gauss seidel vs laplacian ###
+@pytest.mark.parametrize('nx, ny', [
+    (5, 5),
+    (10, 5),
+    (5, 10),
+    (30, 50)
+    ])
+
+def test_compute_phi_2d(nx, ny):
+    n = (nx * ny - 1) / 2
+    rho = np.linspace(-n , n, nx * ny).reshape(ny , nx)
+    laplacian = create_laplacian_2d(nx, ny, 1, 1)
+    phi = compute_phi_2d(rho, [1, 1])
+    rho_calc = -np.dot(laplacian, phi.reshape(-1)).reshape(ny, nx)
+    np.testing.assert_almost_equal(rho, rho_calc, decimal=4)
diff --git a/projects/project_6/Shipra_Yaoyi_Louis_Lucas/LJ.py b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/LJ.py
new file mode 100644
index 0000000..760d012
--- /dev/null
+++ b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/LJ.py
@@ -0,0 +1,67 @@
+import numpy as np
+
+# temporarily define all constants here
+epsilon_lj = 1.
+sigma_lj = 1.
+cutoff_lj = 3. * sigma_lj
+
+def LJ_potential(q):
+    """Calculate the system Lennard-Jones potential based on given particle coordinates.
+    
+    Input:
+        q (iterable of (#dimension,)-arrays): list or array of particle coordinates.
+    
+    Output:
+        pot (float): sum of all pairwise LJ potentials.    
+    """
+    
+    # TODO: input check
+    
+    n_particle = np.shape(q)[0]
+    pot = 0.
+    for i in range(n_particle - 1):
+        for j in range(1, n_particle):
+            pot += _LJ_potential_pair(q[i], q[j])
+    return pot
+
+def _LJ_potential_pair(qi, qj):
+    """Calculate pairwise Lennard-Jones potential.
+        (require global epsilon_lj, sigma_lj, cutoff_lj)
+    """
+    rij = np.linalg.norm(qi - qj)
+    if rij <= 0 or rij > cutoff_lj:
+        return 0.
+    return 4 * epsilon_lj * (sigma_lj ** 12 / rij ** 12 - sigma_lj ** 6 / rij ** 6)
+
+def LJ_force(q):
+    """Calculate the system Lennard-Jones forces based on given particle coordinates.
+    
+    Input:
+        q (iterable of (#dimension,)-arrays): list or array of particle coordinates.
+        k (int): index of the particle for force calc.
+    
+    Output:
+        force ((n, #dimension)-array): system Lennard-Jones forces on particle *0~(n-1)*.    
+    """
+    
+    # TODO: input check
+    
+    n_particle = np.shape(q)[0]
+    ndim = np.shape(q)[1]
+    force = np.zeros((n_particle, ndim))
+    for i in range(n_particle):
+        for j in range(n_particle):
+            # pairwise force calc will return 0 for identical two particles.
+            force[i] += _LJ_force_pair(q[i], q[j])
+    return force
+
+def _LJ_force_pair(qk, qj):
+    """Calculate pairwise Lennard-Jones force on paricle *k*.
+        (require global epsilon_lj, sigma_lj, cutoff_lj)
+    """
+    rjk = np.linalg.norm(qj - qk)
+    ndim = np.shape(qk)[0]
+    if rjk <= 0 or rjk > cutoff_lj:
+        return np.zeros(ndim)
+    # return 4 * epsilon_lj * (12 * sigma_lj ** 12 / rjk ** 14 - 6 * sigma_lj ** 6 / rjk ** 8) * (qj - qk)
+    return 4 * epsilon_lj * (12 * sigma_lj ** 12 / rjk ** 14 - 6 * sigma_lj ** 6 / rjk ** 8) * (qk - qj)
diff --git a/projects/project_6/Shipra_Yaoyi_Louis_Lucas/langevin.py b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/langevin.py
new file mode 100644
index 0000000..3601a3a
--- /dev/null
+++ b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/langevin.py
@@ -0,0 +1,44 @@
+import numpy as np
+
+
+def langevin(
+        force, n_steps, x_init, v_init, mass,
+        time_step=0.001, damping=0.1, beta=1.0):
+    """Langevin integrator for initial value problems
+
+    This function implements the BAOAB algorithm of Benedict Leimkuhler
+    and Charles Matthews. See J. Chem. Phys. 138, 174102 (2013) for
+    further details.
+
+    Arguments:
+        force (function): computes the forces of a single configuration
+        n_steps (int): number of integration steps
+        x_init (numpy.ndarray(n, d)): initial configuration
+        v_init (numpy.ndarray(n, d)): initial velocities
+        mass (numpy.ndarray(n)): particle masses
+        time_step (float): time step for the integration
+        damping (float): damping term, use zero if not coupled
+        beta (float): inverse temperature
+
+    Returns:
+        x (numpy.ndarray(n_steps + 1, n, d)): configuraiton trajectory
+        v (numpy.ndarray(n_steps + 1, n, d)): velocity trajectory
+    """
+    shape = list(x_init.shape)
+    th = 0.5 * time_step
+    thm = 0.5 * time_step / mass[:, None]
+    edt = np.exp(-damping * time_step)
+    sqf = np.sqrt((1.0 - edt ** 2) / (beta * mass))[:, None]
+    x = np.zeros([n_steps + 1] + shape)
+    v = np.zeros_like(x)
+    x[0, :, :] = x_init
+    v[0, :, :] = v_init
+    f = force(x[0])
+    for i in range(n_steps):
+        v[i + 1, :, :] = v[i] + thm * f
+        x[i + 1, :, :] = x[i] + th * v[i + 1]
+        v[i + 1, :, :] = edt * v[i + 1] + sqf * np.random.randn(*shape)
+        x[i + 1, :, :] = x[i + 1] + th * v[i + 1]
+        f[:, :] = force(x[i + 1])
+        v[i + 1, :, :] = v[i + 1] + thm * f
+    return x, v
diff --git a/projects/project_6/Shipra_Yaoyi_Louis_Lucas/mmc_LJ_trap.py b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/mmc_LJ_trap.py
new file mode 100644
index 0000000..653abd9
--- /dev/null
+++ b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/mmc_LJ_trap.py
@@ -0,0 +1,58 @@
+import numpy as np
+##################### declare constants ########################
+N_particles= 4
+mc_steps = 1000
+particle_radius = 1
+potential_depth = 1
+spring_constant = 1
+beta = 1
+dimension = 2
+inital_density = 0.5
+initial_length = (N_particles/inital_density)**(1/dimension)
+scaling = 0.1                                       
+##################### distances #################################
+tril_vector = np.tril_indices(N_particles,-1)
+def get_distances(coordinates):
+    dist_matrix = np.linalg.norm(
+        coordinates[:, None, :] - coordinates[None, :, :],
+        axis=-1)             
+    return dist_matrix[tril_vector]
+################### Lennard-Jones ###############################
+def LJ_potential(coordinates): 
+	u_LJ = 0
+	d = get_distances(coordinates)
+	for _ in range(len(d)):
+		if d[_] == 0:
+			u_LJ = 1e500
+			print('divided by zero')
+		elif d[_] < 2.5 * particle_radius:
+			attractive_part_of_potential = (particle_radius/d[_])**(6.)
+			u_LJ += (attractive_part_of_potential) ** 2. - attractive_part_of_potential
+	return u_LJ * 4 * potential_depth
+################### harmonic_Potential ###########################	
+def harmonic_Potential(coordinates):
+	u_harmonic = 0
+	for _ in range(N_particles):
+		u_harmonic += spring_constant * 0.5 * np.linalg.norm(coordinates[_])**2. 
+	return u_harmonic
+################## Metropolis Monte Carlo ########################
+def mmc():
+	coordinates = np.zeros([mc_steps, N_particles, dimension])
+	u = np.zeros(mc_steps)
+	initial_coordinates = np.random.rand(N_particles, dimension) * initial_length - 0.5 * initial_length * np.ones([N_particles, dimension])
+	coordinates[0] = initial_coordinates
+	for t in range(1,mc_steps):
+		choose_particle = np.random.randint(0, high=N_particles, size=None, dtype='l')
+		random_vector = np.zeros([N_particles,dimension])	
+		random_vector[choose_particle] = np.random.rand(dimension)-0.5
+		coordinates_ = coordinates[t-1] + scaling * random_vector / np.linalg.norm(random_vector)
+		u_  = LJ_potential(coordinates_)     + harmonic_Potential(coordinates_)
+		u[t-1] = LJ_potential(coordinates[t-1])	   + harmonic_Potential(coordinates[t-1])
+		if u_ <= u[t-1] \
+		or np.random.rand() < np.exp(beta*(u[t-1]-u_)):
+			coordinates[t] = coordinates_
+			u[t] = u_  
+		else:
+			coordinates[t] = coordinates[t-1]
+			u[t] = u[t-1]
+	return u, coordinates
diff --git a/projects/project_6/Shipra_Yaoyi_Louis_Lucas/movie_maker.py b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/movie_maker.py
new file mode 100644
index 0000000..23efd6a
--- /dev/null
+++ b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/movie_maker.py
@@ -0,0 +1,33 @@
+################### create Movie of configs #####################
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib
+from mmc_LJ_trap import *
+matplotlib.use("Agg")
+u, coordinates = mmc()
+import matplotlib.animation as manimation
+
+FFMpegWriter = manimation.writers['ffmpeg']
+metadata = dict(title='Movie Test', artist='Matplotlib',
+                comment='Movie support!')
+writer = FFMpegWriter(fps=20, metadata=metadata)
+
+fig = plt.figure()
+l, = plt.plot([], [], 'ro')
+plt.setp(l, markersize=30)
+plt.setp(l, markerfacecolor='C0') 
+
+
+plt.xlim(-5, 5)
+plt.ylim(-5, 5)
+
+x, y = np.zeros(N_particles), np.zeros(N_particles)
+
+with writer.saving(fig, "movie_of_configs.mp4", 100):
+	for _ in range(mc_steps):
+		coord = coordinates[_]
+		for i in range(N_particles):
+			x[i] = coord[i, -1]
+			y[i] = coord[i, 0]
+		l.set_data(x, y)
+		writer.grab_frame()
diff --git a/projects/project_6/Shipra_Yaoyi_Louis_Lucas/movie_of_configs.mp4 b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/movie_of_configs.mp4
new file mode 100644
index 0000000..3c85421
Binary files /dev/null and b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/movie_of_configs.mp4 differ
diff --git a/projects/project_6/Shipra_Yaoyi_Louis_Lucas/plotting.ipynb b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/plotting.ipynb
new file mode 100644
index 0000000..a6442a0
--- /dev/null
+++ b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/plotting.ipynb
@@ -0,0 +1,192 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from pot_force import all_force\n",
+    "from langevin import langevin\n",
+    "\n",
+    "\n",
+    "\n",
+    "dimensions = 2\n",
+    "time_step = 0.001\n",
+    "damping = 0\n",
+    "beta = 0.01"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def x_init_arange(n_particles):\n",
+    "    return np.arange(n_particles * dimensions).reshape(n_particles, dimensions)\n",
+    "\n",
+    "def x_init_circular(n_particles):\n",
+    "    angle = 0\n",
+    "    x_init = np.zeros((n_particles, dimensions))\n",
+    "    for i in range(n_particles):\n",
+    "        x_init[i, :] = [radius * np.cos(angle), radius * np.sin(angle)]\n",
+    "        angle += 2 * np.pi / n_particles\n",
+    "    return x_init\n",
+    "\n",
+    "def plot(n_particles, x_init_function):\n",
+    "    x_init = x_init_function(n_particles)\n",
+    "    v_init = np.arange(n_particles * dimensions).reshape(n_particles, dimensions)\n",
+    "    mass = np.ones(n_particles)\n",
+    "    \n",
+    "    x_sol, v_sol = langevin(all_force, n_steps, x_init, v_init,\n",
+    "                            mass, time_step * 10, damping, beta)\n",
+    "\n",
+    "    r_sol = np.linalg.norm(x_sol, axis=2)\n",
+    "    r_sol = r_sol.reshape(-1)\n",
+    "    r_sol_mean = np.mean(r_sol)\n",
+    "    r_sol_max = np.ndarray.max(r_sol)\n",
+    "    r_sol_min = np.ndarray.min(r_sol)\n",
+    "\n",
+    "    print('mean distance from centre:', r_sol_mean)\n",
+    "    print('minimum distance from centre:', r_sol_min)\n",
+    "    print('maximum distance from centre:', r_sol_max)\n",
+    "    \n",
+    "\n",
+    "    plt.figure(1, figsize=(12, 5))\n",
+    "    plt.subplot(122)\n",
+    "    plt.title('trajectories')\n",
+    "\n",
+    "    for i in range(n_particles):\n",
+    "        plt.plot(x_sol[:,i,0], x_sol[:,i,1])\n",
+    "    \n",
+    "    histogram = []\n",
+    "    iter_list = np.linspace(r_sol_min, r_sol_max, 40)\n",
+    "    for i in range(39):\n",
+    "        n = 0\n",
+    "        for j in range(len(r_sol)):\n",
+    "            if iter_list[i] < r_sol[j]:\n",
+    "                if r_sol[j] < iter_list[i+1]:\n",
+    "                    n += 1\n",
+    "        histogram.append(n)\n",
+    "    \n",
+    "    y_pos = np.arange(len(histogram))\n",
+    "    plt.subplot(121)\n",
+    "    plt.title('radial distribution')\n",
+    "    plt.bar(iter_list[1:], histogram)\n",
+    "    \n",
+    "    plt.show()\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#plots trajectory for n = 2, 3, 4 partiles with initial position [n, n + 1] for n_steps \n",
+    "\n",
+    "n_steps = int(1E3)\n",
+    "\n",
+    "for i in [2, 3, 4]:\n",
+    "    plot(i, x_init_arange)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#plots trajectory for n = 2, 3, 4 partiles with initial position [n, n + 1] for n_steps \n",
+    "\n",
+    "n_steps = int(3E3)\n",
+    "\n",
+    "for i in [2, 3, 4]:\n",
+    "    plot(i, x_init_arange)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#plots trajectory for n = 2, 3, 4 particles evenly distributed on circle of given radius\n",
+    "\n",
+    "n_steps = int(1E3)\n",
+    "\n",
+    "radius = 5\n",
+    "for i in [2,3,4]:\n",
+    "    plot(i, x_init_circular)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#plots trajectory for n = 2, 3 particles evenly distributed on circle of given radius\n",
+    "\n",
+    "n_steps = int(1E3)\n",
+    "\n",
+    "radius = 10\n",
+    "for i in [2,3]:\n",
+    "    plot(i, x_init_circular)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/projects/project_6/Shipra_Yaoyi_Louis_Lucas/pot_force.py b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/pot_force.py
new file mode 100644
index 0000000..55162ff
--- /dev/null
+++ b/projects/project_6/Shipra_Yaoyi_Louis_Lucas/pot_force.py
@@ -0,0 +1,38 @@
+import numpy as np
+from LJ import LJ_force, LJ_potential
+
+# temporarily define all constants here
+k_harmonic = 1
+
+def harmonic_potential(q):
+    """Calculate a well-like harmonic potential based on given particle coordinates.
+    Centered at the origin.
+    
+    Input:
+        q (iterable of (#dimension,)-arrays): list or array of particle coordinates.
+    
+    Output:
+        pot (float): sum of all harmonic potentials.    
+    """
+    
+    # TODO: input check
+    
+    return k_harmonic * np.square(q).sum()
+
+def harmonic_force(q):
+    """Calculate a well-like harmonic force based on given particle coordinates.
+    
+    Input:
+        q (iterable of (#dimension,)-arrays): list or array of particle coordinates.
+    
+    Output:
+        force (float): force from the harmonic well on each particle.    
+    """
+    
+    return -2 * k_harmonic * q
+
+def all_potential(q):
+    return LJ_potential(q) + harmonic_potential(q)
+
+def all_force(q):
+    return LJ_force(q) + harmonic_force(q)