MDBox module

Simple 2D Lennard-Jones liquid simulation

box -- boundary conditions in (x,y) direction: ( x > 0 , _ ) = closed boundary at x and -x ( x = 0 , _ ) = no boundary in x-direction ( x < 0 , _ ) = periodic boundary at x and -x n_bubbles -- number of bubbles (int, default 100) masses -- array containing particle masses (default 1) pos -- array of shape (3, n_bubbles) containing positions vel -- array of scalars with maximum random velocity radius -- relax -- minimize potential energy on initialization using simulated annealing grid -- initial spacing in a grid (default True)

Example usage (in a notebook):

import hylleraas.bubblebox as bb

%matplotlib notebook

system = mdbox(n_bubbles = 100, size = (10,10)) #initialize 10 by 10 closed box containing 100 bubbles

system.run() #run simulation interactively

Source code in bubblebox/mdbox.py

class mdbox():
    """
    Simple 2D Lennard-Jones liquid simulation

    Keyword arguments:

    box              -- boundary conditions in (x,y) direction:
                        ( x > 0 , _ ) = closed boundary at x and -x
                        ( x = 0 , _ ) = no boundary in x-direction
                        ( x < 0 , _ ) = periodic boundary at x and -x
    n_bubbles        -- number of bubbles (int, default 100)
    masses           -- array containing particle masses (default 1)
    pos              -- array of shape (3, n_bubbles) containing positions
    vel              -- array of scalars with maximum random velocity
    radius           -- 
    relax            -- minimize potential energy on initialization using simulated annealing
    grid             -- initial spacing in a grid (default True)

    Example usage (in a notebook):

    import hylleraas.bubblebox as bb

    %matplotlib notebook

    system = mdbox(n_bubbles = 100, size = (10,10)) #initialize 10 by 10 closed box containing 100 bubbles

    system.run() #run simulation interactively 


    """

    def __init__(self, n_bubbles = 100, masses = None, vel = 0.0, size = (0,0), grid = True, pair_list = True, fields = False, sphere_collisions = False):
        # Initialize system
        self.sphere_collisions = sphere_collisions


        # Boundary conditions

        self.size = np.array(size)
        self.size2 = np.array(size)*2
        self.ndim = len(size)

        # obsolete parameter for hard-sphere collisions
        self.radius = .1
        self.n_bubbles = n_bubbles

        # list to hold force algorithms 
        self.force_stack = (no_force, lj_force)

        # array to keep track of forces and interaction parameters
        self.interactions = np.ones((self.n_bubbles, self.n_bubbles, 3), dtype = float)

        # array to keep track of the positions of the bubbles
        self.pos = np.zeros((self.ndim,self.n_bubbles), dtype = float)

        # arrange bubbles in grid
        n_bubbles_axis = int(np.ceil(n_bubbles**(1/self.ndim)))
        grid_axes = []
        for i in range(self.size.shape[0]):
            grid_axes.append(np.linspace(-self.size[i], self.size[i], n_bubbles_axis+1)[:-1])

        self.pos = np.array(np.meshgrid(*grid_axes)).reshape(self.ndim, int(n_bubbles_axis**self.ndim))[:,:self.n_bubbles]

        # move to center
        self.pos = self.pos - np.mean(self.pos, axis = 1)[:, None]




        if masses is None:
            self.masses = np.ones(self.n_bubbles, dtype = int)

            self.masses_inv = np.array(self.masses, dtype = float)**-1
            #self.n_bubbles = n_bubbles

        else:
            self.masses = masses
            self.n_bubbles = len(masses)
            self.set_interactions(self.masses)
            self.masses_inv = np.array(self.masses, dtype = float)**-1







        self.pos_old = self.pos*1 # retain previous timesteps to compute velocities

        # all bubbles active by default
        self.active = np.ones(self.pos.shape[1], dtype = bool)


        # Set velocities (and half-step velocity for v-verlet iterations)
        self.vel = np.random.multivariate_normal(np.zeros(self.ndim), vel*np.eye(self.ndim), self.n_bubbles).T
        self.vel[:] -= np.mean(self.vel, axis = 1)[:, None]
        self.vel_ = self.vel # verlet integrator velocity at previous timestep

        # Integrator 
        self.advance = self.advance_vverlet

        # Algorithm for force calculation
        self.forces = forces
        self.force = lj_force
        self.r2_cut = 9.0 #distance cutoff squared for force-calculation

        # Time and timestep
        self.t = 0
        self.dt = 0.001

        # Collision counter
        self.col = 0

        # Prime for special first iteration
        self.first_iteration = True
        self.iteration = 0

        # Pair list for efficient force-calculation
        self.pair_list = None
        if pair_list:
            self.pair_list_update_frequency = 10
            self.pair_list_buffer = 1.2
            self.pair_list = compute_pair_list(self.pos, self.r2_cut*self.pair_list_buffer, self.size2)

        # fields
        self.fields = fields
        if self.fields:
            self.nbins, self.Z = 20, np.zeros((20,20,2), dtype = float)



    def resize_box(self, size):
        # New boundary conditions
        #self.Lx = size[0]
        #self.Ly = size[1]

        #self.L2x = 2*size[0]
        #self.L2y = 2*size[1]

        self.size = np.array(size)
        self.size2 = np.array(size)*2







    def advance_vverlet(self):
        """
        Advance one step in time according to the Velocity-Verlet algorithm
        """
        if self.first_iteration:
            self.Fn = self.forces(self.pos, self.size2, self.interactions, r2_cut = self.r2_cut, pair_list = self.pair_list)
            self.first_iteration = False
        if self.pair_list is not None:
            if self.iteration % self.pair_list_update_frequency == 0:
                self.pair_list = compute_pair_list(self.pos, self.r2_cut*self.pair_list_buffer, self.size2)


        Fn = self.Fn

        # field
        if self.fields:
            self.Z, dv = gen_field(self, self.nbins, self.Z)
            #self.vel_ = .99*self.vel_ + .01*dv.T
            Fn += .1*dv.T

        self.d_pos = self.vel_*self.dt + .5*Fn*self.dt**2*self.masses_inv

        pos_new = self.pos + self.d_pos

        forces_new = self.forces(pos_new, self.size2, self.interactions, r2_cut = self.r2_cut, pair_list = self.pair_list)

        self.vel_ = self.vel_ + .5*(forces_new + Fn)*self.dt*self.masses_inv

        self.Fn = forces_new


        # impose PBC
        for i in range(self.ndim):
            if self.size[i]<0:
                pos_new[i, :]  = (pos_new[i,:] + self.size[i]) % (self.size2[i]) - self.size[i]


        # impose wall and collision boundary conditions
        self.vel_, self.col = collisions(pos_new, self.vel_, screen = 10.0, radius = self.radius, size2 = self.size2, masses = self.masses, pair_list = self.pair_list, sphere_collisions = self.sphere_collisions)

        #update arrays (in order to retain velocity)
        self.vel = (pos_new - self.pos_old)/(2*self.dt)
        self.pos_old[:] = self.pos
        self.pos[:, self.active] = pos_new[:, self.active]

        # Track time
        self.t += self.dt
        self.iteration += 1


    def advance_euler(self):
        """
        Advance one step in time according to the explicit Euler algorithm
        """
        self.vel_ += self.forces(self.pos, self.size2, self.interactions, r2_cut = self.r2_cut, pair_list = self.pair_list)*self.dt*self.masses_inv
        #self.vel += self.forces(self.pos, self.size2, self.interactions, self.r2_cut, self.force)*self.dt*self.masses_inv
        self.pos[:,self.active] += self.vel_[:,self.active]*self.dt


        # impose PBC
        for i in range(self.ndim):
            if self.size[i]<0:
                self.pos[i, :]  = (self.pos[i,:] + self.size[i]) % (self.size2[i]) - self.size[i]

        # impose wall and collision bounary conditions
        self.vel_, self.col = collisions(self.pos, self.vel_, screen = 10.0, radius = self.radius, size2 = self.size2, masses = self.masses, pair_list = self.pair_list, sphere_collisions = self.sphere_collisions)

        #update arrays (in order to retain velocity)
        #self.vel = (pos_new - self.pos_old)/(2*self.dt)
        #self.pos_old[:] = self.pos
        #self.pos[:, self.active] = pos_new[:, self.active]

        # Track time
        self.t += self.dt
        self.iteration += 1

    def compute_energy(self):
        """
        Compute total energy of system
        """
        return self.compute_potential_energy() + self.compute_kinetic_energy()


    def compute_potential_energy(self):
        """
        Compute total potential energy in system
        """
        #return lj_potential(self.pos, self.interactions, L2x = self.L2x, L2y = self.L2y, r2_cut = self.r2_cut)
        return lj_potential(self.pos, self.size2, interactions = self.interactions, r2_cut = self.r2_cut, force = self.force, pair_list = self.pair_list)

    def compute_kinetic_energy(self):
        """
        Compute total kinetic energy in system
        """
        # Vektorisert funksjon med hensyn på ytelse
        return .5*np.sum(self.masses*np.sum(self.vel_**2, axis = 0))

    def kinetic_energies(self):
        """
        Compute kinetic energy of each bubble
        """
        # Vektorisert funksjon med hensyn på ytelse
        return .5*self.masses*np.sum(self.vel_**2, axis = 0)




    def evolve(self, t = 1.0):
        """
        Let system evolve in time for t seconds
        """
        t1 = self.t+t
        while self.t<t1:
            self.advance()

    """
    Setters
    """

    def set_size(self, size):
        """
        Set new size of box / alias for resize
        """
        self.resize_box(size)

    def set_charges(self, charges):
        """
        Set charges for Coulomb interactions

        Arguments
        ===
        - charges: a numpy.ndarray (or list) containing float or integer charges for all particles

        """
        assert(len(charges) == self.n_bubbles), "Number of charges must equal number of bubbles (%i)." %self.n_bubbles
        for i in range(self.n_bubbles):
            for j in range(self.n_bubbles):
                self.interactions[i,j] = np.array([charges[i],charges[j]])

    def set_masses(self, masses, bubbles = None):
        """
        Set masses 

        Arguments
        ===
        - masses: a numpy.ndarray (or list) containing float or integer charges for all particles

        """
        if bubbles is None:
            self.masses[:] = masses
        else:
            self.masses[bubbles] = masses

        self.masses_inv = np.array(self.masses, dtype = float)**-1


    def set_forces(self, force,  bubbles_a = None, bubbles_b = None, force_params = np.array([0,0])):
        """
        Set the force acting between bubbles_a and bubbles_b (or all)
        """


        if bubbles_a is None:
            bubbles_a_ =  np.array(np.arange(self.n_bubbles), dtype = int)
        else:
            if type(bubbles_a) is np.ndarray and bubbles_a.dtype is np.dtype('bool'):
                bubbles_a_ = np.array(np.arange(self.n_bubbles)[bubbles_a], dtype = int)
            else:
                bubbles_a_ = bubbles_a

        if bubbles_b is None:
            bubbles_b_ =  np.array(np.arange(self.n_bubbles), dtype = int)
        else:
            if type(bubbles_b) is np.ndarray and bubbles_b.dtype is np.dtype('bool'):
                    bubbles_b_ = np.array(np.arange(self.n_bubbles)[bubbles_b], dtype = int)
            else:
                bubbles_b_ = bubbles_b





        nx, ny = np.array(np.meshgrid(bubbles_a_,bubbles_b_)).reshape(2,-1)
        # set pointer to the force function in the force stack
        self.interactions[nx,ny,:] = force()

        # set parameters to use for the specific pairs of particles
        self.interactions[nx,ny,1] = force_params[0]
        self.interactions[nx,ny,2] = force_params[1]
    def set_vel(self, vel):
        """
        set manually the velocities of the system
        """
        assert(np.all(vel.shape==self.vel.shape)), "incorrect shape of velocities"
        self.vel = vel
        self.vel_ = self.vel # verlet integrator velocity at previous timestep

    def set_pos(self, pos):
        assert(pos.shape[0] == len(self.size)), "incorrect shape of pos"
        assert(pos.shape[1] == self.n_bubbles), "incorrect shape of pos"
        self.pos = pos







    def set_interactions(self, masses):
        """
        Placeholder for proper parametrization of LJ-interactions
        """
        #self.interactions = np.ones((self.n_bubbles, self.n_bubbles, 3), dtype = float)


        epsilons = np.linspace(.5,10,100)
        sigmas   = np.linspace(1,np.sqrt(2), 100)



        for i in range(self.n_bubbles):
            mi = int(masses[i])
            for j in range(i+1, self.n_bubbles):

                mj = int(masses[j])

                eps = np.sqrt(epsilons[mi]*epsilons[mj])
                sig = sigmas[mi] + sigmas[mj]


                self.interactions[i,j] = [1, eps, sig]



    """
    Visualization tools (some obsolete to be deleted)
    """
    def visualize_state(self, axis = False, figsize = None):
        """
        Show an image of the current state with positions, velocities (as arrows) and boundaries of the box.
        """
        if figsize is None:
            figsize = (6,6)
            if self.L2x != 0 and self.L2y != 0:
                figsize = (4, 4*np.abs(self.L2y/self.L2x))

            plt.rcParams["figure.figsize"] = figsize



        col = colorscheme()

        plt.figure(figsize = figsize)
        plt.plot([-self.Lx, self.Lx, self.Lx, -self.Lx, -self.Lx],[-self.Ly, -self.Ly, self.Ly, self.Ly, -self.Ly], color = (0,0,0), linewidth = 2)
        plt.plot(self.pos[0], self.pos[1], 'o', alpha = .4, markersize = 8*1.8, color = col.getcol(.5))
        plt.plot(self.pos[0], self.pos[1], '.', alpha = 1, markersize = 10, color = (0,0,0))


        for i in range(len(self.vel[0])):
            plt.plot([self.pos[0,i], self.pos[0,i] + self.vel[0,i]],[self.pos[1,i], self.pos[1,i] + self.vel[1,i]], "-", color = (0,0,0))

            th = np.arctan2(self.vel[1,i],self.vel[0,i])
            plt.text(self.pos[0,i] + self.vel[0,i],self.pos[1,i] + self.vel[1,i], "▲", rotation = -90+360*th/(2*np.pi),ha = "center", va = "center") #, color = (0,0,0), fontsize = 20, rotation=0, ha = "center", va = "center")

        plt.xlim(-self.Lx-1, self.Lx+1)
        if self.Lx == 0:
            plt.xlim(-11, 11)
        plt.ylim(-self.Ly-1, self.Ly+1)
        if self.Ly == 0:
            plt.ylim(-11, 11)

        if not axis:
            plt.axis("off")
        plt.show()

    #def run(self, n_steps_per_vis = 5, interval = 1):
    #    run_system = animated_system(system = self, n_steps_per_vis=n_steps_per_vis, interval = interval)
    #    plt.show()

    def run(self, nsteps, n_iterations_per_step = 1):
        for i in range(nsteps):
            for j in range(n_iterations_per_step):
                self.advance()
            self.update_view()

    def view(self, viewer = ev.MDView):
        self.mview = viewer(self)
        return self.mview

    def evince(self, realism = True, dof = True, sao = True, focus = 10, aperture = 0.0001, max_blur = 0.001):
        """
        Create high-quality 3D view using Evince
        """
        # extract bonds (harmonic oscillator interactions)
        bonds = ev.spotlight.extract_bonds(self)

        self.mview = ev.SpotlightView(self, bonds = bonds, realism = realism, dof = dof, sao = sao, focus = focus, aperture=aperture, max_blur=max_blur)
        return self.mview


    def update_view(self):
        self.mview.pos = self.pos.T.tolist()


    def save_view(self, filename, title =""):


        embed_minimal_html(filename, [self.mview], title)


    # "algebra"

    def __add__(self, other):
        if type(other) in [float, int, np.ndarray]:
            ret = copy.deepcopy(self)
            if type(other) is np.ndarray:
                ret.pos = self.pos + other[:, None]
            else:
                ret.pos = self.pos + other
            return ret


        else: 


            #np.concatenate([b.interactions, b.interactions], axis = 3).shape


            interactions = np.zeros((self.n_bubbles+other.n_bubbles, self.n_bubbles+other.n_bubbles, 3), dtype = float )
            interactions[:self.n_bubbles,:self.n_bubbles,:] = self.interactions
            interactions[self.n_bubbles:,self.n_bubbles:,:] = other.interactions



            # determine unique interactions
            masses = np.concatenate([self.masses, other.masses])
            ui = np.unique(masses )

            UI = np.zeros((len(ui), len(ui), 3), dtype = float)
            #UI = {}


            # locate force definitions
            for i in range(len(ui)):
                for j in range(1,len(ui)):
                    for l in range(len(masses)):
                        if l==i:
                            for k in range(len(masses)):
                                if k==j:
                                    UI[i,j] = interactions[l,k]
                                    UI[j,i] = interactions[l,k]
                                    break
                            break


            # set interactions accordingly
            for i in range(len(self.masses), len(self.masses)+len(other.masses)):
                for j in range(len(self.masses)):
                    interactions[i,j] = UI[ui==masses[i], ui==masses[j]]

            # determine size
            size = []
            for i in range(len(self.size)):
                size.append(np.sign(np.sign(self.size[i])+np.sign(other.size[i]))*max(abs(self.size[i]), abs(other.size[i])))




            ret = mdbox(self.n_bubbles+other.n_bubbles, size = size)

            ret.interactions = interactions
            ret.set_masses(masses)
            ret.vel_ = np.concatenate([self.vel_, other.vel_], axis = 1)
            ret.vel = np.concatenate([self.vel, other.vel], axis = 1)
            ret.pos = np.concatenate([self.pos, other.pos], axis = 1)

            return ret

    def __radd__(self, other):
        return self.__add__(other)

    def __mul__(self, other):
        if type(other) is int:
            return extend_system(self, other)

        if type(other) is float:
            ret = copy.deepcopy(self)
            ret.pos *= other
        return ret

    def __rmul__(self, other):
        return self.__mul__(other)

    def __matmul__(self, other):
        ret = copy.deepcopy(self)
        ret.pos = self.pos.T.dot(other).T
        return ret

    def __rmatmul__(self, other):
        ret = copy.deepcopy(self)
        ret.pos = other.dot(self.pos)
        return ret

advance_euler()

Advance one step in time according to the explicit Euler algorithm

Source code in bubblebox/mdbox.py

def advance_euler(self):
    """
    Advance one step in time according to the explicit Euler algorithm
    """
    self.vel_ += self.forces(self.pos, self.size2, self.interactions, r2_cut = self.r2_cut, pair_list = self.pair_list)*self.dt*self.masses_inv
    #self.vel += self.forces(self.pos, self.size2, self.interactions, self.r2_cut, self.force)*self.dt*self.masses_inv
    self.pos[:,self.active] += self.vel_[:,self.active]*self.dt


    # impose PBC
    for i in range(self.ndim):
        if self.size[i]<0:
            self.pos[i, :]  = (self.pos[i,:] + self.size[i]) % (self.size2[i]) - self.size[i]

    # impose wall and collision bounary conditions
    self.vel_, self.col = collisions(self.pos, self.vel_, screen = 10.0, radius = self.radius, size2 = self.size2, masses = self.masses, pair_list = self.pair_list, sphere_collisions = self.sphere_collisions)

    #update arrays (in order to retain velocity)
    #self.vel = (pos_new - self.pos_old)/(2*self.dt)
    #self.pos_old[:] = self.pos
    #self.pos[:, self.active] = pos_new[:, self.active]

    # Track time
    self.t += self.dt
    self.iteration += 1

advance_vverlet()

Advance one step in time according to the Velocity-Verlet algorithm

Source code in bubblebox/mdbox.py

def advance_vverlet(self):
    """
    Advance one step in time according to the Velocity-Verlet algorithm
    """
    if self.first_iteration:
        self.Fn = self.forces(self.pos, self.size2, self.interactions, r2_cut = self.r2_cut, pair_list = self.pair_list)
        self.first_iteration = False
    if self.pair_list is not None:
        if self.iteration % self.pair_list_update_frequency == 0:
            self.pair_list = compute_pair_list(self.pos, self.r2_cut*self.pair_list_buffer, self.size2)


    Fn = self.Fn

    # field
    if self.fields:
        self.Z, dv = gen_field(self, self.nbins, self.Z)
        #self.vel_ = .99*self.vel_ + .01*dv.T
        Fn += .1*dv.T

    self.d_pos = self.vel_*self.dt + .5*Fn*self.dt**2*self.masses_inv

    pos_new = self.pos + self.d_pos

    forces_new = self.forces(pos_new, self.size2, self.interactions, r2_cut = self.r2_cut, pair_list = self.pair_list)

    self.vel_ = self.vel_ + .5*(forces_new + Fn)*self.dt*self.masses_inv

    self.Fn = forces_new


    # impose PBC
    for i in range(self.ndim):
        if self.size[i]<0:
            pos_new[i, :]  = (pos_new[i,:] + self.size[i]) % (self.size2[i]) - self.size[i]


    # impose wall and collision boundary conditions
    self.vel_, self.col = collisions(pos_new, self.vel_, screen = 10.0, radius = self.radius, size2 = self.size2, masses = self.masses, pair_list = self.pair_list, sphere_collisions = self.sphere_collisions)

    #update arrays (in order to retain velocity)
    self.vel = (pos_new - self.pos_old)/(2*self.dt)
    self.pos_old[:] = self.pos
    self.pos[:, self.active] = pos_new[:, self.active]

    # Track time
    self.t += self.dt
    self.iteration += 1

compute_energy()

Compute total energy of system

Source code in bubblebox/mdbox.py

def compute_energy(self):
    """
    Compute total energy of system
    """
    return self.compute_potential_energy() + self.compute_kinetic_energy()

compute_kinetic_energy()

Compute total kinetic energy in system

Source code in bubblebox/mdbox.py

def compute_kinetic_energy(self):
    """
    Compute total kinetic energy in system
    """
    # Vektorisert funksjon med hensyn på ytelse
    return .5*np.sum(self.masses*np.sum(self.vel_**2, axis = 0))

compute_potential_energy()

Compute total potential energy in system

Source code in bubblebox/mdbox.py

def compute_potential_energy(self):
    """
    Compute total potential energy in system
    """
    #return lj_potential(self.pos, self.interactions, L2x = self.L2x, L2y = self.L2y, r2_cut = self.r2_cut)
    return lj_potential(self.pos, self.size2, interactions = self.interactions, r2_cut = self.r2_cut, force = self.force, pair_list = self.pair_list)

evince(realism=True, dof=True, sao=True, focus=10, aperture=0.0001, max_blur=0.001)

Create high-quality 3D view using Evince

Source code in bubblebox/mdbox.py

def evince(self, realism = True, dof = True, sao = True, focus = 10, aperture = 0.0001, max_blur = 0.001):
    """
    Create high-quality 3D view using Evince
    """
    # extract bonds (harmonic oscillator interactions)
    bonds = ev.spotlight.extract_bonds(self)

    self.mview = ev.SpotlightView(self, bonds = bonds, realism = realism, dof = dof, sao = sao, focus = focus, aperture=aperture, max_blur=max_blur)
    return self.mview

evolve(t=1.0)

Let system evolve in time for t seconds

Source code in bubblebox/mdbox.py

def evolve(self, t = 1.0):
    """
    Let system evolve in time for t seconds
    """
    t1 = self.t+t
    while self.t<t1:
        self.advance()

kinetic_energies()

Compute kinetic energy of each bubble

Source code in bubblebox/mdbox.py

def kinetic_energies(self):
    """
    Compute kinetic energy of each bubble
    """
    # Vektorisert funksjon med hensyn på ytelse
    return .5*self.masses*np.sum(self.vel_**2, axis = 0)

set_charges(charges)

Set charges for Coulomb interactions

Arguments

• charges: a numpy.ndarray (or list) containing float or integer charges for all particles

Source code in bubblebox/mdbox.py

def set_charges(self, charges):
    """
    Set charges for Coulomb interactions

    Arguments
    ===
    - charges: a numpy.ndarray (or list) containing float or integer charges for all particles

    """
    assert(len(charges) == self.n_bubbles), "Number of charges must equal number of bubbles (%i)." %self.n_bubbles
    for i in range(self.n_bubbles):
        for j in range(self.n_bubbles):
            self.interactions[i,j] = np.array([charges[i],charges[j]])

set_forces(force, bubbles_a=None, bubbles_b=None, force_params=np.array([0, 0]))

Set the force acting between bubbles_a and bubbles_b (or all)

Source code in bubblebox/mdbox.py

def set_forces(self, force,  bubbles_a = None, bubbles_b = None, force_params = np.array([0,0])):
    """
    Set the force acting between bubbles_a and bubbles_b (or all)
    """


    if bubbles_a is None:
        bubbles_a_ =  np.array(np.arange(self.n_bubbles), dtype = int)
    else:
        if type(bubbles_a) is np.ndarray and bubbles_a.dtype is np.dtype('bool'):
            bubbles_a_ = np.array(np.arange(self.n_bubbles)[bubbles_a], dtype = int)
        else:
            bubbles_a_ = bubbles_a

    if bubbles_b is None:
        bubbles_b_ =  np.array(np.arange(self.n_bubbles), dtype = int)
    else:
        if type(bubbles_b) is np.ndarray and bubbles_b.dtype is np.dtype('bool'):
                bubbles_b_ = np.array(np.arange(self.n_bubbles)[bubbles_b], dtype = int)
        else:
            bubbles_b_ = bubbles_b





    nx, ny = np.array(np.meshgrid(bubbles_a_,bubbles_b_)).reshape(2,-1)
    # set pointer to the force function in the force stack
    self.interactions[nx,ny,:] = force()

    # set parameters to use for the specific pairs of particles
    self.interactions[nx,ny,1] = force_params[0]
    self.interactions[nx,ny,2] = force_params[1]

set_interactions(masses)

Placeholder for proper parametrization of LJ-interactions

Source code in bubblebox/mdbox.py

def set_interactions(self, masses):
    """
    Placeholder for proper parametrization of LJ-interactions
    """
    #self.interactions = np.ones((self.n_bubbles, self.n_bubbles, 3), dtype = float)


    epsilons = np.linspace(.5,10,100)
    sigmas   = np.linspace(1,np.sqrt(2), 100)



    for i in range(self.n_bubbles):
        mi = int(masses[i])
        for j in range(i+1, self.n_bubbles):

            mj = int(masses[j])

            eps = np.sqrt(epsilons[mi]*epsilons[mj])
            sig = sigmas[mi] + sigmas[mj]


            self.interactions[i,j] = [1, eps, sig]

set_masses(masses, bubbles=None)

Set masses

Arguments

• masses: a numpy.ndarray (or list) containing float or integer charges for all particles

Source code in bubblebox/mdbox.py

def set_masses(self, masses, bubbles = None):
    """
    Set masses 

    Arguments
    ===
    - masses: a numpy.ndarray (or list) containing float or integer charges for all particles

    """
    if bubbles is None:
        self.masses[:] = masses
    else:
        self.masses[bubbles] = masses

    self.masses_inv = np.array(self.masses, dtype = float)**-1

set_size(size)

Set new size of box / alias for resize

Source code in bubblebox/mdbox.py

def set_size(self, size):
    """
    Set new size of box / alias for resize
    """
    self.resize_box(size)

set_vel(vel)

set manually the velocities of the system

Source code in bubblebox/mdbox.py

def set_vel(self, vel):
    """
    set manually the velocities of the system
    """
    assert(np.all(vel.shape==self.vel.shape)), "incorrect shape of velocities"
    self.vel = vel
    self.vel_ = self.vel # verlet integrator velocity at previous timestep

visualize_state(axis=False, figsize=None)

Show an image of the current state with positions, velocities (as arrows) and boundaries of the box.

Source code in bubblebox/mdbox.py

def visualize_state(self, axis = False, figsize = None):
    """
    Show an image of the current state with positions, velocities (as arrows) and boundaries of the box.
    """
    if figsize is None:
        figsize = (6,6)
        if self.L2x != 0 and self.L2y != 0:
            figsize = (4, 4*np.abs(self.L2y/self.L2x))

        plt.rcParams["figure.figsize"] = figsize



    col = colorscheme()

    plt.figure(figsize = figsize)
    plt.plot([-self.Lx, self.Lx, self.Lx, -self.Lx, -self.Lx],[-self.Ly, -self.Ly, self.Ly, self.Ly, -self.Ly], color = (0,0,0), linewidth = 2)
    plt.plot(self.pos[0], self.pos[1], 'o', alpha = .4, markersize = 8*1.8, color = col.getcol(.5))
    plt.plot(self.pos[0], self.pos[1], '.', alpha = 1, markersize = 10, color = (0,0,0))


    for i in range(len(self.vel[0])):
        plt.plot([self.pos[0,i], self.pos[0,i] + self.vel[0,i]],[self.pos[1,i], self.pos[1,i] + self.vel[1,i]], "-", color = (0,0,0))

        th = np.arctan2(self.vel[1,i],self.vel[0,i])
        plt.text(self.pos[0,i] + self.vel[0,i],self.pos[1,i] + self.vel[1,i], "▲", rotation = -90+360*th/(2*np.pi),ha = "center", va = "center") #, color = (0,0,0), fontsize = 20, rotation=0, ha = "center", va = "center")

    plt.xlim(-self.Lx-1, self.Lx+1)
    if self.Lx == 0:
        plt.xlim(-11, 11)
    plt.ylim(-self.Ly-1, self.Ly+1)
    if self.Ly == 0:
        plt.ylim(-11, 11)

    if not axis:
        plt.axis("off")
    plt.show()