{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Creating Circle Traps using external forces."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"import os\n",
"import sys\n",
"\n",
"import pandas as pd\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib as mpl\n",
"\n",
"sys.path.insert(0, '../../../magcolloids')\n",
"\n",
"import magcolloids as mgc\n",
"\n",
"from IPython.display import HTML"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%reload_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"ureg = mgc.ureg"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# parameters\n",
"N_particles = 20\n",
"radius_circle = 50*ureg.um\n",
"radius_particles = 2*ureg.um"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"#initial_conditions\n",
"phi_p = np.linspace(0,2*np.pi,N_particles+1)[:-1]\n",
"x = np.cos(phi_p)*radius_circle.to(ureg.um).magnitude\n",
"y = np.sin(phi_p)*radius_circle.to(ureg.um).magnitude\n",
"z = 0*phi_p\n",
"\n",
"initial_conditions = np.array([x,y,z]).transpose()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"region = 3*radius_circle.to(ureg.um).magnitude*np.array([1,1,1])\n",
"region[2] = (2*radius_particles+1*ureg.um).to(ureg.um).magnitude"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"particles = mgc.particles(\n",
" initial_conditions*ureg.um,\n",
" radius = 2*ureg.um,\n",
" susceptibility = 0,\n",
" diffusion=0.07*ureg.um**2/ureg.s,\n",
" density = 1000*ureg.kg/ureg.m**3,\n",
" temperature=300*ureg.K)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'1.00e-03picogram / microsecond ** 2'"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k = 1*ureg.fN/ureg.nm\n",
"\"%2.2e\"%k.to(\"pg/us**2\").magnitude+str(k.to(\"pg/us**2\").units)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Defining an external force.\n",
"\n",
"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. \n",
"\n",
"The details of how the forces can be calculated can be found in the lammps document [about variables](https://lammps.sandia.gov/doc/variable.html). 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. \n",
"\n",
"Below I show an example in which I create a ring harmonic trap. This is defined as:\n",
"$$F = -k \\left(\\left|r\\right|-R\\right) \\mathbf{e}_r$$\n",
"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. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"force_calculation = \"\"\"\n",
"# Remember that time is in microsecond\n",
"# Force is in pg*um/us^2\n",
"# k should be in pg/us^2\n",
"# r, cx and cy are in um\n",
"\n",
"variable k atom 1e-3\n",
"\n",
"variable R atom 50\n",
"\n",
"variable rho atom sqrt(x*x+y*y)\n",
"variable ex atom x/v_rho\n",
"variable ey atom y/v_rho\n",
"\n",
"variable dr atom v_rho-v_R\n",
"\n",
"variable Fx atom -v_k*v_dr*v_ex\n",
"variable Fy atom -v_k*v_dr*v_ey\n",
"variable Fz atom 0\n",
"\n",
"\"\"\""
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"traps = mgc.ext_force(calculation = force_calculation, variable = \"v_F\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"using the variables defined above, we create the `ext_force` object, and then we give it as input to the `world` object. "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"field = mgc.field(magnitude = 0*ureg.mT, frequency = 10*ureg.Hz, angle = 15*ureg.degrees)\n",
"\n",
"world = mgc.world(particles, ext_force = traps, temperature = 300*ureg.K,\n",
" region=region*ureg.um, boundaries = ['f','f','f'], walls = [False,False,True],\n",
" dipole_cutoff = 20*ureg.um)\n",
"\n",
"sim = mgc.sim(dir_name = \".\", file_name = \"ring_trap_test_ring_0\",\n",
" timestep = 1e-4*ureg.s, framerate = 30*ureg.Hz, total_time = 60*ureg.s,\n",
" particles = particles, world = world, field = field)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"sim.generate_scripts()\n",
"sim.run()"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"sim.load()\n",
"trj = sim.lazy_read[::10]"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"HTML(mgc.display_animation_direct(sim,trj,speedup=1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Now let's try to shift the ring in space"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The above definition can be extended to shifth the ring to a different center"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"force_calculation = \"\"\"\n",
"# Remember that time is in microsecond\n",
"# Force is in pg*um/us^2\n",
"# k should be in pg/us^2\n",
"# r, cx and cy are in um\n",
"\n",
"variable k atom 1e-3\n",
"\n",
"variable R atom 50\n",
"\n",
"variable cx atom 5\n",
"variable cy atom 0\n",
"\n",
"variable rho atom sqrt((x-v_cx)^2+(y-v_cy)^2)\n",
"variable ex atom (x-v_cx)/v_rho\n",
"variable ey atom (y-v_cy)/v_rho\n",
"\n",
"variable dr atom v_rho-v_R\n",
"\n",
"variable Fx atom -v_k*v_dr*v_ex\n",
"variable Fy atom -v_k*v_dr*v_ey\n",
"variable Fz atom 0\n",
"\n",
"\"\"\""
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"traps = mgc.ext_force(calculation = force_calculation, variable = \"v_F\")"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"field = mgc.field(magnitude = 0*ureg.mT, frequency = 10*ureg.Hz, angle = 15*ureg.degrees)\n",
"\n",
"world = mgc.world(particles, ext_force = traps, temperature = 300*ureg.K,\n",
" region=region*ureg.um, boundaries = ['f','f','f'], walls = [False,False,True],\n",
" dipole_cutoff = 20*ureg.um)\n",
"\n",
"sim = mgc.sim(dir_name = \".\", file_name = \"ring_trap_test_ring_1\",\n",
" timestep = 1e-4*ureg.s, framerate = 30*ureg.Hz, total_time = 6*ureg.s,\n",
" particles = particles, world = world, field = field)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"sim.generate_scripts()\n",
"sim.run()"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"sim.load()\n",
"trj = sim.lazy_read[::1]"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"HTML(mgc.display_animation_direct(sim,trj,speedup=1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Now let's try to add a second ring"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"force_calculation = \"\"\"\n",
"# Remember that time is in microsecond\n",
"# Force is in pg*um/us^2\n",
"# k should be in pg/us^2\n",
"# r, cx and cy are in um\n",
"\n",
"variable k atom 1e-3\n",
"variable cutoff atom 1\n",
"\n",
"variable R1 atom 50\n",
"variable R2 atom 2\n",
"\n",
"variable cx1 atom 0\n",
"variable cy1 atom 0\n",
"\n",
"variable cx2 atom 54\n",
"variable cy2 atom 0\n",
"\n",
"variable rho1 atom sqrt((x-v_cx1)^2+(y-v_cy1)^2)\n",
"variable ex1 atom (x-v_cx1)/v_rho1\n",
"variable ey1 atom (y-v_cy1)/v_rho1\n",
"\n",
"variable dr1 atom v_rho1-v_R1\n",
"\n",
"variable rho2 atom sqrt((x-v_cx2)^2+(y-v_cy2)^2)\n",
"variable ex2 atom (x-v_cx2)/v_rho2\n",
"variable ey2 atom (y-v_cy2)/v_rho2\n",
"\n",
"variable dr2 atom v_rho2-v_R2\n",
"\n",
"variable Fx atom -v_k*v_dr1*v_ex1*(v_dr1\n",
" \n",
" Your browser does not support the video tag.\n",
""
],
"text/plain": [
""
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"HTML(mgc.display_animation_direct(sim,trj,speedup=1))"
]
},
{
"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.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}