Creating Circle Traps using external forces.

import os
import sys

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

sys.path.insert(0, '../../../magcolloids')

import magcolloids as mgc

from IPython.display import HTML

%reload_ext autoreload
%autoreload 2

ureg = mgc.ureg

# parameters
N_particles = 20
radius_circle = 50*ureg.um
radius_particles = 2*ureg.um

#initial_conditions
phi_p = np.linspace(0,2*np.pi,N_particles+1)[:-1]
x = np.cos(phi_p)*radius_circle.to(ureg.um).magnitude
y = np.sin(phi_p)*radius_circle.to(ureg.um).magnitude
z = 0*phi_p

initial_conditions = np.array([x,y,z]).transpose()

region = 3*radius_circle.to(ureg.um).magnitude*np.array([1,1,1])
region[2] = (2*radius_particles+1*ureg.um).to(ureg.um).magnitude

particles = mgc.particles(
    initial_conditions*ureg.um,
    radius = 2*ureg.um,
    susceptibility = 0,
    diffusion=0.07*ureg.um**2/ureg.s,
    density = 1000*ureg.kg/ureg.m**3,
    temperature=300*ureg.K)

k = 1*ureg.fN/ureg.nm
"%2.2e"%k.to("pg/us**2").magnitude+str(k.to("pg/us**2").units)

'1.00e-03picogram / microsecond ** 2'

Defining an external force.

I’ve added the ext_force object to the magcolloids package. It allows us to add a force to every particle in every timestep. Below I show how to calculate a force. In the end we need to create three variables, Fx, Fy and Fz which represent the three components of the force that will be added to each particle. The name of the variable can be changed using the variable argument.

The details of how the forces can be calculated can be found in the lammps document about variables. Roughly, we can define a variable (for example k), and then call it in further variable definitions as v_k. There are a few internal variables for each atom, for example the coordinates x, y and z. These don’t need to be prepended by v. The variable comand includes the specification atom which tells lammps to calculate a value for each atom. There are other styles of defining a variable (for example equal which produces a single number), but I think they don’t work with the addforce fix.

Below I show an example in which I create a ring harmonic trap. This is defined as:

\[F = -k \left(\left|r\right|-R\right) \mathbf{e}_r\]

where \(R\) is the radius of the ring, and \(\mathbf{e}_r\) is the unit vector in the radial direction. There are no vector operations, so this function needs to be calculated for each component.

force_calculation = """
# Remember that time is in microsecond
# Force is in pg*um/us^2
# k should be in pg/us^2
# r, cx and cy are in um

variable k atom 1e-3

variable R atom 50

variable rho atom sqrt(x*x+y*y)
variable ex atom x/v_rho
variable ey atom y/v_rho

variable dr atom v_rho-v_R

variable Fx atom -v_k*v_dr*v_ex
variable Fy atom -v_k*v_dr*v_ey
variable Fz atom 0

"""

traps = mgc.ext_force(calculation = force_calculation, variable = "v_F")

using the variables defined above, we create the ext_force object, and then we give it as input to the world object.

field = mgc.field(magnitude = 0*ureg.mT, frequency = 10*ureg.Hz, angle = 15*ureg.degrees)

world = mgc.world(particles, ext_force = traps, temperature = 300*ureg.K,
                  region=region*ureg.um, boundaries = ['f','f','f'], walls = [False,False,True],
                  dipole_cutoff = 20*ureg.um)

sim = mgc.sim(dir_name = ".", file_name = "ring_trap_test_ring_0",
       timestep = 1e-4*ureg.s, framerate = 30*ureg.Hz, total_time = 60*ureg.s,
       particles = particles, world = world, field = field)

sim.generate_scripts()
sim.run()

sim.load()
trj = sim.lazy_read[::10]

HTML(mgc.display_animation_direct(sim,trj,speedup=1))

Now let’s try to shift the ring in space

The above definition can be extended to shifth the ring to a different center

force_calculation = """
# Remember that time is in microsecond
# Force is in pg*um/us^2
# k should be in pg/us^2
# r, cx and cy are in um

variable k atom 1e-3

variable R atom 50

variable cx atom 5
variable cy atom 0

variable rho atom sqrt((x-v_cx)^2+(y-v_cy)^2)
variable ex atom (x-v_cx)/v_rho
variable ey atom (y-v_cy)/v_rho

variable dr atom v_rho-v_R

variable Fx atom -v_k*v_dr*v_ex
variable Fy atom -v_k*v_dr*v_ey
variable Fz atom 0

"""

traps = mgc.ext_force(calculation = force_calculation, variable = "v_F")

field = mgc.field(magnitude = 0*ureg.mT, frequency = 10*ureg.Hz, angle = 15*ureg.degrees)

world = mgc.world(particles, ext_force = traps, temperature = 300*ureg.K,
                  region=region*ureg.um, boundaries = ['f','f','f'], walls = [False,False,True],
                  dipole_cutoff = 20*ureg.um)

sim = mgc.sim(dir_name = ".", file_name = "ring_trap_test_ring_1",
       timestep = 1e-4*ureg.s, framerate = 30*ureg.Hz, total_time = 6*ureg.s,
       particles = particles, world = world, field = field)

sim.generate_scripts()
sim.run()

sim.load()
trj = sim.lazy_read[::1]

HTML(mgc.display_animation_direct(sim,trj,speedup=1))

Now let’s try to add a second ring

I explored a little bit a possibility to add a second ring, which is tangent to the first at a single point. In this setup, rings will need a cutoff. I haven’t been able to get the particles to diffuse in the second ring. This needs to be thought more carefuly.

force_calculation = """
# Remember that time is in microsecond
# Force is in pg*um/us^2
# k should be in pg/us^2
# r, cx and cy are in um

variable k atom 1e-3
variable cutoff atom 1

variable R1 atom 50
variable R2 atom 2

variable cx1 atom 0
variable cy1 atom 0

variable cx2 atom 54
variable cy2 atom 0

variable rho1 atom sqrt((x-v_cx1)^2+(y-v_cy1)^2)
variable ex1 atom (x-v_cx1)/v_rho1
variable ey1 atom (y-v_cy1)/v_rho1

variable dr1 atom v_rho1-v_R1

variable rho2 atom sqrt((x-v_cx2)^2+(y-v_cy2)^2)
variable ex2 atom (x-v_cx2)/v_rho2
variable ey2 atom (y-v_cy2)/v_rho2

variable dr2 atom v_rho2-v_R2

variable Fx atom -v_k*v_dr1*v_ex1*(v_dr1<v_cutoff)-v_k*v_dr2*v_ex2*(v_dr2<v_cutoff)
variable Fy atom -v_k*v_dr1*v_ey1*(v_dr1<v_cutoff)-v_k*v_dr2*v_ey2*(v_dr2<v_cutoff)
variable Fz atom 0

"""

traps = mgc.ext_force(calculation = force_calculation, variable = "v_F")

field = mgc.field(magnitude = 0*ureg.mT, frequency = 10*ureg.Hz, angle = 15*ureg.degrees)

world = mgc.world(particles, ext_force = traps, temperature = 300*ureg.K,
                  region=region*ureg.um, boundaries = ['f','f','f'], walls = [False,False,True],
                  dipole_cutoff = 20*ureg.um)

sim = mgc.sim(dir_name = ".", file_name = "ring_trap_test_ring_2",
       timestep = 1e-4*ureg.s, framerate = 30*ureg.Hz, total_time = 600*ureg.s,
       particles = particles, world = world, field = field)

sim.generate_scripts()
sim.run()

sim.load()
trj = sim.lazy_read[::100]

HTML(mgc.display_animation_direct(sim,trj,speedup=1))