Solution to project 6

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

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)

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

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()

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)

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)

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

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

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()