tag:blogger.com,1999:blog-47877516311168623222024-03-14T03:31:35.942-07:00Crows and CatsMathematical Modelling in Nature, Life, and ScienceAndrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.comBlogger36125tag:blogger.com,1999:blog-4787751631116862322.post-25380598613775169942014-06-19T00:12:00.001-07:002014-06-19T00:12:38.907-07:00My blog has moved!My blog has moved to <a href="http://andrew-white.com/blog">here</a>Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-75595784452936450172014-03-03T06:11:00.000-08:002014-04-16T12:34:57.833-07:00Replaying a Lammps Trajectory<b>EDIT: </b>Turns out I just didn't google the correct words, LAMMPS does have one. It's the command rerun. Oops.<br />
<br />
<a href="http://lammps.sandia.gov/">LAMMPS</a>, the molecular dynamics engine, doesn't have a replay function allowing one to go back and calculate properties on a trajectory, especially forces and energies. This came up in my work and I got started on writing some new LAMMPS code to accomplish this. I wanted to be able to calculate forces with a different potential on a trajectory. Anyway, I realized you can accomplish this just by using the LAMMPS looping variables. No new code necessary. It can be done by using the <code>read_dump</code> command to load particular trajectory frame. A sample input script is below <br />
<div>
<br /></div>
<div>
<br /></div>
<a name='more'></a><script src="https://gist.github.com/whitead/9325454.js"></script>Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com1tag:blogger.com,1999:blog-4787751631116862322.post-21965143408093663322014-02-06T04:56:00.002-08:002014-02-06T04:56:28.583-08:00Predicting the Future of SimulationsIn my day job, molecular modeling, I simulate the position and velocity of every atom in biological molecules. As you can imagine, doing this yields nearly complete information about whatever biological system I'm studying and allows me to predict the behavior of proteins. Well, that's the theory. It's more complex than that and often fails spectacularly in practice. Ignore that for now. Let's focus on the systems that can be simulated. In 1982 CHARMM, one the first molecular dynamics tool capable of simulating tiny proteins, was released. It was possible to simulate a small protein for one tenth of one billionth of one second. Now we can simulate a protein for one thousandth of a second. That's an increase of 10^7 in the speed of molecular simulations in about 30 years. In 1990, we were able to simulate lysozyme, a small protein, for few billionths of a seconds. Now we can simulate a virus for a few billionths of a second. That's an increase in about 10^4 mass of molecular simulations in 30 years. Assuming continued exponential growth (see plots below), here are my upcoming interesting milestones in computer simulation<br />
<br />
<br />
<style>
#chemistry {
border-collapse:collapse;
}
#chemistry tr.alt td
{
color:#000000;
}
#chemistry th
{
font-size:1.1em;
text-align:left;
padding-top:5px;
padding-bottom:4px;
background-color:#21bb53;
color:#ffffff;
}
#chemistry td, #chemistry th
{
font-size:1em;
border:1px solid #21bb53;
padding:3px 7px 2px 7px;
}
</style>
<br />
<table id="chemistry">
<tbody>
<tr>
<th>Year </th>
<th>Event</th>
</tr>
<tr>
<td>2023</td>
<td>The first simulation of every atom in a ribosome for 10 milliseconds</td>
</tr>
<tr>
<td>2038</td>
<td>The first simulation of every atom in a virus for 1 second</td>
</tr>
<tr>
<td>2131</td>
<td>The first simulation of every atom of a cell for 1 hour</td>
</tr>
<tr>
<td>2191</td>
<td>We can simulate every atom of every cell in a heart for 1 heartbeat</td>
</tr>
<tr>
<td>2227</td>
<td>Every atom of every cell in a human brain can be simulated as fast they move in a human brain</td>
</tr>
</tbody></table>
<a name='more'></a><br />
<br />
Here's the data used to make these predictions:
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgtgLmiZkn6_qaGz7Dg_gnpERR-eWsWK-_gm9d9IoBJUWvenOqbSCGO8OIgmLjWmAm6CQug78wBigqqen8HFu9Zuen98Y6P0Ed7l8AmJQIQ72lWOQbeIrIQrgQ_p3cW9ZtRzug3W6Gh9yNQ/s1600/plot_2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgtgLmiZkn6_qaGz7Dg_gnpERR-eWsWK-_gm9d9IoBJUWvenOqbSCGO8OIgmLjWmAm6CQug78wBigqqen8HFu9Zuen98Y6P0Ed7l8AmJQIQ72lWOQbeIrIQrgQ_p3cW9ZtRzug3W6Gh9yNQ/s1600/plot_2.png" height="300" width="400" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEipLjA7wtESR5w8xTRMDRIcIKZZTyfYVGDoLqHcD-yTHuOJCWS2Mnl3g-IQijsPDCbKoxiqKyamCgYQrjrnw0XKL2bWrvBEAZ-P-hUy55o1Wc3mTS7Ly9UQGhX5VR5ITEMRpkre4qEZseXR/s1600/plot_3.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEipLjA7wtESR5w8xTRMDRIcIKZZTyfYVGDoLqHcD-yTHuOJCWS2Mnl3g-IQijsPDCbKoxiqKyamCgYQrjrnw0XKL2bWrvBEAZ-P-hUy55o1Wc3mTS7Ly9UQGhX5VR5ITEMRpkre4qEZseXR/s1600/plot_3.png" height="300" width="400" /></a><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjaXkh5YDohiOFVVH6_q_PUZmpxgj7vytnCy9h4TQebQkmdhSDnfoE4StY7vXRbpWRlFDokPKd9nuwV0Ol_6eDgj0eAHSnQrmKQsqRvzpoXyzm_4l2ShabgxPhZXtAzGFKrMBrLqETc8y4R/s1600/plot_4.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjaXkh5YDohiOFVVH6_q_PUZmpxgj7vytnCy9h4TQebQkmdhSDnfoE4StY7vXRbpWRlFDokPKd9nuwV0Ol_6eDgj0eAHSnQrmKQsqRvzpoXyzm_4l2ShabgxPhZXtAzGFKrMBrLqETc8y4R/s1600/plot_4.png" height="300" width="400" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
This post is a work in progress and I'll throw the links to the sources up when I have a chance. The key numbers are doubling time in timescale of simulations every 14 months and atom number/mass every 26 months.Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-62009487385319314232013-12-12T18:36:00.003-08:002013-12-12T18:36:45.167-08:00Radial Distribution Functions of Water from ExperimentSite-site radial distribution functions, sometimes called g(r) for short, are one of the most important measured quantities in molecular simulations. As usual, wikipedia has a <a href="http://en.wikipedia.org/wiki/Radial_distribution_function">nice overview of them.</a> Probably the most famous non-trivial g(r) is the water oxygen-oxygen radial distribution function. It's often taught in statistical mechanics courses that it's possible to compare a simulation g(r) with X-ray or neutron-diffraction experiments with a simple Fourier transform. In fact, this is almost never the case due to sources of error in experiments. Check out the Introduction to <a href="http://www.sciencedirect.com/science/article/pii/S0301010400001798">this article by AK Soper</a> for more info on the kinds of error which prevent the simple Fourier transform from converting experimental structure factors to a g(r).<br />
<br />
In any case, the point of this post is to put the raw data for an experimental g(r) for water. The <a href="http://www.sciencedirect.com/science/article/pii/S0301010400001798">AK Soper</a> article uses the empirical potential structure refinement method for converting experimental neutron diffraction data at 298 K 1 bar into the three radial distribution functions for water (O-O, O-H, H-H). This is the best method I know of for obtaining these radial distribution functions and this is the data used to test the <a href="http://scitation.aip.org.proxy.uchicago.edu/content/aip/journal/jcp/123/23/10.1063/1.2121687">TIP4P/2005 water model.</a> I've extracted the data from their paper and made it <a href="https://gist.github.com/whitead/7938965">available here</a>. Make sure to cite their paper if you use the data in a publication. Hopefully these g(r)s will help you in creating your own new special water models! See the figure below the break:<br />
<br />
<a name='more'></a><br /><br />
I don't think I can legally reproduce the image from the paper on my blog, but compare Figure 6 from the AK Soper article with:<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgKOmWtpxwvniR8tEpiM8wlTprwrfF_MHh0CIwVpjRt_hXyqOknEqcu9vQ8lurh_ijFHwe1WrwpFsgrw0wgc9DaB4Lkcs8tiULCvqo3qhEUkSTXIfr67n-q3p6qLX6zr1yVnohSB8ALwBKm/s1600/reproduced_annotated.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="640" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgKOmWtpxwvniR8tEpiM8wlTprwrfF_MHh0CIwVpjRt_hXyqOknEqcu9vQ8lurh_ijFHwe1WrwpFsgrw0wgc9DaB4Lkcs8tiULCvqo3qhEUkSTXIfr67n-q3p6qLX6zr1yVnohSB8ALwBKm/s640/reproduced_annotated.png" width="492" /></a></div>
<br />Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com1tag:blogger.com,1999:blog-4787751631116862322.post-58166281291823113722013-10-01T09:30:00.001-07:002013-10-01T09:30:37.330-07:00Indexing an Emacs MacroThis post is about how to increment and paste a counter while executing an emacs macro. For some reason, this is not easy to find on the Internet. So I'll write it down on my blog since I keep forgetting.<br />
<br />
<br />
<ol>
<li>Set the start point of the counter <br/> <code> M-: (setq i 1)</code></li>
<li> Start recording macro with <br/> <code> C-x (</code></li>
<li> Insert index as needed with <br/> <code> C-u M-: i</code></li>
<li> Increment index when finished with macro <br/> <code> M-: (setq i (+ i 1)) </code></li>
<li> Complete macro with <br/> <code> C-x )</code></li>
<li> Execute Macro with <br/> <code> C-x e </code></li>
</ol>
To execute the macro multiple times, use <code> ALT-(# of times) C-x e</code>. Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-74381858791518831592013-05-16T23:03:00.001-07:002013-05-17T10:05:14.507-07:00Extending Radial Distributions<a href="http://en.wikipedia.org/wiki/Radial_distribution_function">Radial distribution functions</a> are one of the most important quantities to calculate when running a molecular dynamics simulation. They give good insight into the phase and pair potentials between molecules. Practically, they are calculated according to the following equation:<br /><br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBY2IDtHwJdo4Oj3sJsyDC2QyAmomp7t_XJFIDFW3q4NurDzwLEYX9qSUxx2BZyWMfBRHNXZ-0E3Rl8HoK53QqUIUTtv4lboH4KmHYEGN7sOZURCKNjLpbOUAMWRyAiyEsqLwF_yQ__OBF/s1600/1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBY2IDtHwJdo4Oj3sJsyDC2QyAmomp7t_XJFIDFW3q4NurDzwLEYX9qSUxx2BZyWMfBRHNXZ-0E3Rl8HoK53QqUIUTtv4lboH4KmHYEGN7sOZURCKNjLpbOUAMWRyAiyEsqLwF_yQ__OBF/s320/1.png" /></a></div>
<a name='more'></a><br />
<br />
where <code> v_i </code> is a small volume increment that corresponds to the radius increment <code>r_i</code>, <code> n_i </code> is the average number of particles found in the volume <code>v_i</code>, and <code> rho_bulk</code> is the bulk density of the entire simulation. Essentially, at each incremental volume we're calculating the ratio of particles found to what we would find if the particles in that volume behaved like the bulk density. The incremental volume in which we search for particles may be seen in the image below:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj2G4yy2DeRIb1fEmNRMkoXcT4Nz1AqdWbZQAMQt8eoxho48oP8yF1Kfy9gy6K4HIN6XZjOaL3tsAMTtR5K9iIAtI7Dp3slvNZs85v8ml2uf5k4YghUX3H3lPqsY0YUN_nfHUYtKLJS78JU/s1600/drawing.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj2G4yy2DeRIb1fEmNRMkoXcT4Nz1AqdWbZQAMQt8eoxho48oP8yF1Kfy9gy6K4HIN6XZjOaL3tsAMTtR5K9iIAtI7Dp3slvNZs85v8ml2uf5k4YghUX3H3lPqsY0YUN_nfHUYtKLJS78JU/s320/drawing.png" /></a></div>
However, we notice there is a problem as <code> r_i </code> grows: the concentric circles (or spheres in 3D) begin to intersect with the simulation box that contains the particles. This is is shown below.<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBYJrCTPgSO5bGraPxC5yBOyXbDuQ13TH4lH8cgAKSfnDW0wM1_CHUw88pURfc38zulHvvT5eQuEXSDR8suSNYUnSn3sX3Aboasq_1HLxxIRpAHaCqr_g7JUt2rqdJlaHwrR-7rx_JHpMF/s1600/drawing_2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiBYJrCTPgSO5bGraPxC5yBOyXbDuQ13TH4lH8cgAKSfnDW0wM1_CHUw88pURfc38zulHvvT5eQuEXSDR8suSNYUnSn3sX3Aboasq_1HLxxIRpAHaCqr_g7JUt2rqdJlaHwrR-7rx_JHpMF/s320/drawing_2.png" /></a></div>
Generally this problem is approached by stopping a radial distribution calculation when <code> r_i </code> is equal to half the side length of the simulation box. That actually means we lose a significant amount of statistics, pi / 8 in 3 dimensions to be exact. This in fact may fixed by using the results of my <a href="http://crowsandcats.blogspot.com/2013/04/cube-sphere-intersection-volume.html"> previous post about cube-sphere intersections</a>. In fact, I wrote it specifically about this issue.<br /><br />
The way to allow radial distribution functions to be calculated all the way out until all particles are counted is to calculate the true <code> v_i</code>, which is the intersection volume between the simulation box and the concentric spheres that make up <code> v_i </code>. Using the <a href="https://gist.github.com/whitead/a1dd345a1e9f96805513"> code I wrote</a>, <code> v_i </code> may be calculated according to:<br /><br />
<script src="https://gist.github.com/whitead/67343ed96e5e981eaf28.js"></script>
The relevant radial distribution function calculation now looks like this:<br />
<script src="https://gist.github.com/whitead/766ffca0739fabad6113.js"></script>
Finally, here are the results on Lennard-Jones fluid calculation:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjQDYNwAN-lHU-KTMvEkLhnE-VGcwhJ6r2EIV3TCT7QkudWFv8S3SVqYlBfIWpUQJwe9MMRRvgsCb_ahsrIL97OuAVCaRw7sY1nxinmclOdEpEKZNybyp0pYR1zW6yVupBPLExEl_7dt4gH/s1600/gofr.png" imageanchor="1"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjQDYNwAN-lHU-KTMvEkLhnE-VGcwhJ6r2EIV3TCT7QkudWFv8S3SVqYlBfIWpUQJwe9MMRRvgsCb_ahsrIL97OuAVCaRw7sY1nxinmclOdEpEKZNybyp0pYR1zW6yVupBPLExEl_7dt4gH/s320/gofr.png" /></a>
</div>
<br />
The simulation box ends at the end of the line. As you can see, there is a huge amount of space lost in traditional radial distribution function calculations. Another cool thing about correcting the volume is that the radial distribution does not have to be normalized to converge to 1. Often in these calculations people struggle with double counting and incorrect volumes so that the value is just forced to converge to 1 at the edge of the plot. That is the value which it should converge to. Here's the full radial distribution calculator, which is part of <a href="https://github.com/whitead/Simple-MD"> a simulation engine </a> I'm writing.
<br />
<br />
<script src="https://gist.github.com/whitead/e381e2e6091548bae6a3.js"></script>Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-11296284424130978162013-04-30T18:48:00.001-07:002013-05-13T19:16:40.985-07:00Cube-Sphere Intersection VolumeThis week's math problem was particularly difficult. I wanted to know the intersection volume between a cube and a sphere. Imagine a sphere contained in a cube. The sphere will first touch the cube when its radius is equal to half the width of the cube. Until that point, its volume is just the formula for the volume of a sphere. Once the sphere's radius is <code>2 sqrt(3)</code> times the width of the cube, then it completely engulfs the cube and the intersection of the 2 objects is just the volume of a cube. In between, I've found it looks like this:
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEifyNfwfNGKCsfroY8_vFPVipF8W7N8BDxQIv4vKDcPBOy92pr-7_NXi7P5L4v6yYkCx_Q1EajVhMaC7C3FZ1VdevZU_nfobtpYFh9ed0QQ72cqv-hlr3_-JObP8u4EJFSdu_3SNE-uuu9c/s1600/plot.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEifyNfwfNGKCsfroY8_vFPVipF8W7N8BDxQIv4vKDcPBOy92pr-7_NXi7P5L4v6yYkCx_Q1EajVhMaC7C3FZ1VdevZU_nfobtpYFh9ed0QQ72cqv-hlr3_-JObP8u4EJFSdu_3SNE-uuu9c/s1600/plot.png" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<br />
<a name='more'></a>The gray vertical lines represent when to use which equation. There are in fact 2 other regimes besides the two I discussed. I'm going to discuss these regimes as related to the ratio of the <code> rho </code> and <code> d </code> where the first is the radius of the sphere and the second is half the cube width. Let's see them visually.<br />
<br />
The first regime is the spherical, where <code> rho / d </code> < 1. It looks like a sphere:
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEipghhhqmjwWwAXJc0OIcv_rDR0DKynRpaXz7mxOm8kXz_9l1-epjxWhnKReBlM595wFoTRvWNQZ4NBd-xOlxJB7f0adds1eAVchuURb03bgupX16l_Jo6SyL1lkJ_Gg9excJtdS6yry8fw/s1600/sphere.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="180" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEipghhhqmjwWwAXJc0OIcv_rDR0DKynRpaXz7mxOm8kXz_9l1-epjxWhnKReBlM595wFoTRvWNQZ4NBd-xOlxJB7f0adds1eAVchuURb03bgupX16l_Jo6SyL1lkJ_Gg9excJtdS6yry8fw/s320/sphere.png" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
The second regime has circular "lenses" cut-off from the sphere and occurs while <code> rho / d </code> < <code>sqrt(2):</code><br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiY8M1fwOEb1ndSECaIx6ZpAHNX6ZlOsmxOjmB_I3EgP88Lc2oQq0IpxwF1M1iK2AIcU7qnItX0tDEiR8U6xtMsk7wPMsATjyI72u1K54jWiWPhHyUrbK8i1elCxhyphenhyphenbK2ybjwyuzGMqlycK/s1600/lens.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="180" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiY8M1fwOEb1ndSECaIx6ZpAHNX6ZlOsmxOjmB_I3EgP88Lc2oQq0IpxwF1M1iK2AIcU7qnItX0tDEiR8U6xtMsk7wPMsATjyI72u1K54jWiWPhHyUrbK8i1elCxhyphenhyphenbK2ybjwyuzGMqlycK/s320/lens.png" width="320" /></a></div>
<br />
The third regime is characterized by filling the 8 corners of the cube:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEggoKVH6x8ypVJGMdnimyHZlqELlt8y3yROtOx8AfIXptvsCLT8WEl978IEMBWxitm3TX3yuIpyjM7ECxGWwoVAPWWWqKULISLnnEJRETGVuxt1Be8qm6rG9WkA5jSHK2_-90EKzs49yuDG/s1600/trans.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="180" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEggoKVH6x8ypVJGMdnimyHZlqELlt8y3yROtOx8AfIXptvsCLT8WEl978IEMBWxitm3TX3yuIpyjM7ECxGWwoVAPWWWqKULISLnnEJRETGVuxt1Be8qm6rG9WkA5jSHK2_-90EKzs49yuDG/s320/trans.png" width="320" /></a></div>
<br />
Finally, the last regime is just a cube. This occurs at <code> rho / d </code> > <code> sqrt(3):</code><br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjAuPoAs4c8URcYhXd-UzONwtrL8inIz08FMO36ir4ks3W_QD3aXXmQeXmoyiYYy0gk-gQERhzQvz1JuHKmZUfdVOw16Tgl9RUB9-KmDHg4qd2tE6WZmseXipGoiODZqpoJ56hP8AUa3sXw/s1600/cube.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="180" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjAuPoAs4c8URcYhXd-UzONwtrL8inIz08FMO36ir4ks3W_QD3aXXmQeXmoyiYYy0gk-gQERhzQvz1JuHKmZUfdVOw16Tgl9RUB9-KmDHg4qd2tE6WZmseXipGoiODZqpoJ56hP8AUa3sXw/s320/cube.png" width="320" /></a></div>
<code><br /></code>
The equation for the sphere regime is:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi6ZPeZGpbCOG4oNOyNSZyzGcAyXdvcJH3wCUnxTalnD1eo1ELKWNxyqX9kUw5Bfys07d7vzWvRIwCm2wM1SLOTYHCyzKPXQEoeBJYMNrC7XvABFT6E0h62CkV8juHGaSmiR9Tne8-Vq6pj/s1600/12.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi6ZPeZGpbCOG4oNOyNSZyzGcAyXdvcJH3wCUnxTalnD1eo1ELKWNxyqX9kUw5Bfys07d7vzWvRIwCm2wM1SLOTYHCyzKPXQEoeBJYMNrC7XvABFT6E0h62CkV8juHGaSmiR9Tne8-Vq6pj/s1600/12.png" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
The equation for the "lens" regime may be found by calculating the volume of the removed lenses and then subtracting that from the volume of the sphere.
First, look at the cross section of the lens:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjbYu2VJy-JhmulAcazdXp_4iNs2QrqXJOcMCAXFOqqMZLvV5AzOr_fOJAB26R5sz2FzOVQpPZppaUEuhaUnmfD_mLfkwDz3N4Lpg3cIx7Pbvl4zuB06ihvYtXdgXCqv86P-GbvDlNGvfky/s1600/rect3033.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="114" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjbYu2VJy-JhmulAcazdXp_4iNs2QrqXJOcMCAXFOqqMZLvV5AzOr_fOJAB26R5sz2FzOVQpPZppaUEuhaUnmfD_mLfkwDz3N4Lpg3cIx7Pbvl4zuB06ihvYtXdgXCqv86P-GbvDlNGvfky/s320/rect3033.png" width="320" /></a></div>
<br />
<br />
We can use those rotation transformation we learned in calculus. That is described by the equation:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjr9n-4B0zlwTqb7tbkLOcQ9AcaXASI-1pj_RHndy-gNh57coj8CD8McBfXZyjBkXcSSMFD8cFZMfjRIbi0uD8ceTPhY5tWrrd8I-0qFQdTvgN6NyztHwyxBjB4iwqMTeQiV9A-Auoq9KPn/s1600/13.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjr9n-4B0zlwTqb7tbkLOcQ9AcaXASI-1pj_RHndy-gNh57coj8CD8McBfXZyjBkXcSSMFD8cFZMfjRIbi0uD8ceTPhY5tWrrd8I-0qFQdTvgN6NyztHwyxBjB4iwqMTeQiV9A-Auoq9KPn/s1600/13.png" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg1nQxHw-IW2RZGxhz1AuJfxq7BZU8ngR4Tn0pMPfgtGcWBR0K2jCti9b9bVnDqaJjF9cxeHFvj1bh8vU7mIS84arzCgz5qE2Dm1isQkkfh90bZBzr-wiuHHdMap6C33sqbA5wTdvzft34c/s1600/14.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg1nQxHw-IW2RZGxhz1AuJfxq7BZU8ngR4Tn0pMPfgtGcWBR0K2jCti9b9bVnDqaJjF9cxeHFvj1bh8vU7mIS84arzCgz5qE2Dm1isQkkfh90bZBzr-wiuHHdMap6C33sqbA5wTdvzft34c/s1600/14.png" /></a></div>
<br />
The last regime is much more difficult to deal with. The sphere is first broken into 16ths and then each piece of that is divided into 5 different regions. The 16th piece may be viewed as an animation <a href="http://i.imgur.com/2ToIPEv.gif">here</a> to visualize it. I'll break down each region with an animation showing how it changes as a function of the ratio of the radius to the cube width and the resulting volume of the area.<br />
The first region is defined by:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjOusJZoh-YvDEdhas0D0PDu5nIDosXPFeBjwWunXINkawTsbOg_SzsJ8AmhLAmuqd_6hJQun3V4cRwNoc6ECeFG3djuhSacLJ6yBhDGnRSJGnGf_bm1eVIIs6fvAbVXvLCvauFkAcZSowB/s1600/2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjOusJZoh-YvDEdhas0D0PDu5nIDosXPFeBjwWunXINkawTsbOg_SzsJ8AmhLAmuqd_6hJQun3V4cRwNoc6ECeFG3djuhSacLJ6yBhDGnRSJGnGf_bm1eVIIs6fvAbVXvLCvauFkAcZSowB/s1600/2.png" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhJhrdjrKeK1tnbkpbD4Of58Wyqv1kqLYuNasBUnxYAS3MbQrKP2eQWlBu5-4ALMcGGvH851Mg97I06YOfuIJtwTUFLG96T_VU1E4fT1AG-UauVdLi5yOzFAZs68tesyebqZC8do6wDoGAT/s1600/3.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhJhrdjrKeK1tnbkpbD4Of58Wyqv1kqLYuNasBUnxYAS3MbQrKP2eQWlBu5-4ALMcGGvH851Mg97I06YOfuIJtwTUFLG96T_VU1E4fT1AG-UauVdLi5yOzFAZs68tesyebqZC8do6wDoGAT/s1600/3.png" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgYp8gcN-KRt7dv5dK3XqC3HuZXqc8s3Sxg4QAWTwPPJIk7mTqx2AO9ctJWQukJtl7csWHQNCv7G0avTM7h3PCAHDoI0Cx7hvYNFSL_ml1sD3-9JWZrmiqm_NaB1EyF_8gpqfDo3k1-dQEv/s1600/1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgYp8gcN-KRt7dv5dK3XqC3HuZXqc8s3Sxg4QAWTwPPJIk7mTqx2AO9ctJWQukJtl7csWHQNCv7G0avTM7h3PCAHDoI0Cx7hvYNFSL_ml1sD3-9JWZrmiqm_NaB1EyF_8gpqfDo3k1-dQEv/s1600/1.png" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
Here is a gif of the region:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgbpXLh30XVChbs5d-5SINRX1ctFnTPLy_HcRr_3NRAWd-U3bNBlxVpAwOtf5OqYdmbU66DneZ9AC8cdMZELbPZNZnH3Ze87-LVWJJqUyIeF8uLzNKDcD4w3_NlEBJuzwSFuKGL6S2Rpa57/s1600/region1.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgbpXLh30XVChbs5d-5SINRX1ctFnTPLy_HcRr_3NRAWd-U3bNBlxVpAwOtf5OqYdmbU66DneZ9AC8cdMZELbPZNZnH3Ze87-LVWJJqUyIeF8uLzNKDcD4w3_NlEBJuzwSFuKGL6S2Rpa57/s1600/region1.gif" /></a></div>
<br />
The second region is defined by:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgwJGeMQhXnItGNuYgRfvx0tzxPsf-UO_iFWyKG4fTjXH1kweEe2altSNBnm7zxhyphenhyphencIZPmNqw4QJyAiRzohBszYugooHMlvb8HEQQ4IA-c6Onf3ET1h1Ko3cuZ31naGyr9TvGK-abind8WX/s1600/4.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgwJGeMQhXnItGNuYgRfvx0tzxPsf-UO_iFWyKG4fTjXH1kweEe2altSNBnm7zxhyphenhyphencIZPmNqw4QJyAiRzohBszYugooHMlvb8HEQQ4IA-c6Onf3ET1h1Ko3cuZ31naGyr9TvGK-abind8WX/s1600/4.png" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEijD12a_QYgQXNHQbdDcjPD4NJNNSo6nyybHOje3y8Rbd5qCTOAR7IiUcyQUYCYDtTxk50kj1yfRlw1z5uNcBTZ2bO85Y7uI8J3nZ0ckLp6stmp7ef9ZjBsbeIhrkKmnw60UPStLkjNYP6I/s1600/5.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEijD12a_QYgQXNHQbdDcjPD4NJNNSo6nyybHOje3y8Rbd5qCTOAR7IiUcyQUYCYDtTxk50kj1yfRlw1z5uNcBTZ2bO85Y7uI8J3nZ0ckLp6stmp7ef9ZjBsbeIhrkKmnw60UPStLkjNYP6I/s1600/5.png" /></a></div>
<br />
Here is a gif of the region:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiLWqqR_-xkMq4qklEyUb-qZl-pbNP_RHN0WpS1tocMYy3XCD_-lSkfONYvv1nIxmI_oQhqGdyaO-fWi1FxcLZ1Ma9H0XLf_E9TmmRF6xbjOMPWHY3DEmNtW3NWySe4gYfHRD0sqk2gkWVG/s1600/region2.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiLWqqR_-xkMq4qklEyUb-qZl-pbNP_RHN0WpS1tocMYy3XCD_-lSkfONYvv1nIxmI_oQhqGdyaO-fWi1FxcLZ1Ma9H0XLf_E9TmmRF6xbjOMPWHY3DEmNtW3NWySe4gYfHRD0sqk2gkWVG/s1600/region2.gif" /></a></div>
<br />
The third region is defined by:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi0h0ntKrIEzqMJVGf8uHYlR7H8jIjxDvGnMg0xwTYTp7q3NomKryT-Pp1q7Sp-dCvJuc57u7siH-orgcGfqe8xXrVZalT1HPM8qRK3oVoTBqhHgG_LefYgpZRHh24gEh0hau0dMTbxA0hz/s1600/6.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi0h0ntKrIEzqMJVGf8uHYlR7H8jIjxDvGnMg0xwTYTp7q3NomKryT-Pp1q7Sp-dCvJuc57u7siH-orgcGfqe8xXrVZalT1HPM8qRK3oVoTBqhHgG_LefYgpZRHh24gEh0hau0dMTbxA0hz/s1600/6.png" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjiTx79tt3y7dG90XMq3_VowmB-TbEIib_1JoVAD40aS0cpMGDE7bPfW2aH3ii8Sn5DOR0C0K3QH8Lr7qb4HbOVKuV7GreEWFnkVspjZSpCbPz1l34dhL0twNfs63jOPs3hmUtmBISmCN70/s1600/7.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjiTx79tt3y7dG90XMq3_VowmB-TbEIib_1JoVAD40aS0cpMGDE7bPfW2aH3ii8Sn5DOR0C0K3QH8Lr7qb4HbOVKuV7GreEWFnkVspjZSpCbPz1l34dhL0twNfs63jOPs3hmUtmBISmCN70/s1600/7.png" /></a></div>
<br />
Here is a gif of the region:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg6QcfLaRcPuv0AkJNNNBeW1j-8KRzKyme4J-GIewDuxEMP0b7_Z93x8U3wt8DIxkc0826Tb_8ACq7kxAy3zu6dwvaF_9xSVMgVqy8IZ4RL_i7b8P-RgYZre8_WbpyDKqDmgjbqPZZzyfsD/s1600/region3.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg6QcfLaRcPuv0AkJNNNBeW1j-8KRzKyme4J-GIewDuxEMP0b7_Z93x8U3wt8DIxkc0826Tb_8ACq7kxAy3zu6dwvaF_9xSVMgVqy8IZ4RL_i7b8P-RgYZre8_WbpyDKqDmgjbqPZZzyfsD/s1600/region3.gif" /></a></div>
<br />
The fourth region is defined by:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj2FyxkbRfn7hHDk2iU9IfOlV2eLd8oxpPeOJ_a4vsdz4eupIhnVx7xfS3mhKEAL1Nc8uVArV-0z9V6M6SYVoiyP0AN688avtesTkXqYfGfgvJ4CdZT-EqUFQHNrao18WeAYDUQn9n3nQEr/s1600/8.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj2FyxkbRfn7hHDk2iU9IfOlV2eLd8oxpPeOJ_a4vsdz4eupIhnVx7xfS3mhKEAL1Nc8uVArV-0z9V6M6SYVoiyP0AN688avtesTkXqYfGfgvJ4CdZT-EqUFQHNrao18WeAYDUQn9n3nQEr/s1600/8.png" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjHQSH8SJV4WuBjpOEvXAm_gyyzDeCIzdfevanJ9AFN0FtuGE_YxgIBvi19flZxC8j8SmBwBvkuIAQjHVZtQJVzC8t4AOiMxA4p2R8EOvF5xfPem4J_4CMzHLs3DO4A7RlZ5eGA88UUrPwz/s1600/9.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjHQSH8SJV4WuBjpOEvXAm_gyyzDeCIzdfevanJ9AFN0FtuGE_YxgIBvi19flZxC8j8SmBwBvkuIAQjHVZtQJVzC8t4AOiMxA4p2R8EOvF5xfPem4J_4CMzHLs3DO4A7RlZ5eGA88UUrPwz/s1600/9.png" /></a></div>
<br />
Here is a gif of the region:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgRm64f_j1vVi_yRaMIRrZyxf1sIaLVtzexH9XfqA-ZuIVwHccBrYrTwBdJogk_JvlZh0UyTBPV9HqX_UVakeR6_U9gipIM0zGWNCKFuhdrvqFb6-ITTLV_xSs4gIupICkdAMmvsvs6pMRF/s1600/region4.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgRm64f_j1vVi_yRaMIRrZyxf1sIaLVtzexH9XfqA-ZuIVwHccBrYrTwBdJogk_JvlZh0UyTBPV9HqX_UVakeR6_U9gipIM0zGWNCKFuhdrvqFb6-ITTLV_xSs4gIupICkdAMmvsvs6pMRF/s1600/region4.gif" /></a></div>
<br />
The fifth region is defined by:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhb57RKhMdgVfakrGGi_R6s_TsEEkbxwd5GpEb_624i0NaC0MmdPrFiaG18ClSqIcMpa2orABdbWYbGbJDNlZYsNhrTmL4yjDerSXS-x9vblmqaQSORrjsosqx5FVncJGkOenAWrSnWq7Hq/s1600/10.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhb57RKhMdgVfakrGGi_R6s_TsEEkbxwd5GpEb_624i0NaC0MmdPrFiaG18ClSqIcMpa2orABdbWYbGbJDNlZYsNhrTmL4yjDerSXS-x9vblmqaQSORrjsosqx5FVncJGkOenAWrSnWq7Hq/s1600/10.png" /></a><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhXSFrkS74o15IDEhfgDC9Y6f6DsTj9HR9-v11tXpFrMGa27sspGku179lZWv8cxNZTShNvsUE4eajRvnEwxiRyF8x19QzFzQsbDOaWTkiVFWjdef_JgveTDjiL5Xt97d0_l8Rcw2jW3p9b/s1600/11.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhXSFrkS74o15IDEhfgDC9Y6f6DsTj9HR9-v11tXpFrMGa27sspGku179lZWv8cxNZTShNvsUE4eajRvnEwxiRyF8x19QzFzQsbDOaWTkiVFWjdef_JgveTDjiL5Xt97d0_l8Rcw2jW3p9b/s1600/11.png" /></a></div>
<br />
Here is a gif of the region:<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjVFlcyMqyguKGvFyb1jAOAY80fgvRDdw0W2-Wt4P37ztIUwnO3QvJqL2s-JGE2oBx-_pIM8p58SgizndDikpb-yo-jLS2hpjOSHUU7C8PvG_OIwgW2wt1UfaakY8TBJOsgRaYZu9XuRaOb/s1600/region6.gif" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjVFlcyMqyguKGvFyb1jAOAY80fgvRDdw0W2-Wt4P37ztIUwnO3QvJqL2s-JGE2oBx-_pIM8p58SgizndDikpb-yo-jLS2hpjOSHUU7C8PvG_OIwgW2wt1UfaakY8TBJOsgRaYZu9XuRaOb/s1600/region6.gif" /></a></div>
<br />
Finally, all the regimes/regions may be combined into a python script which calculates the cube-sphere intersection volume. I should mention that I'm not that good at calculus and I think the two unevaluated integrals have analytic forms which I just can't find yet. That's one of the problems of doing integration by hand; you have to be better than me. I'm considering investing in a Mathematica license after doing these integrals. Here's the python code:
<br />
<script src="https://gist.github.com/whitead/a1dd345a1e9f96805513.js"></script>Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com4tag:blogger.com,1999:blog-4787751631116862322.post-14192636662391648722013-04-25T08:05:00.001-07:002013-04-25T08:05:01.535-07:00Bacteria in BlenderI made this while drinking my morning coffee. I just stretched a UV Sphere, stuck two layers of hair on it (flagella, cilia) and added a few modifiers to get it bumpy. The texture is just some glass, Voronoi textures, and diffuse shaders. Meh, it's OK.
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgot9KT8RTZTdJEWVi7y9XOhr6gfsiJvzRi6rjkEshVt_cxgnkrRTGP7uxv06OiLojl96UuwvMnHcvPAJmEdfLZs0fBwC48lNT10diDOVoo7ze4XFwPYSbhuOZPTVMoFMeToz5OPGsz-95B/s1600/bacteri.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="225" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgot9KT8RTZTdJEWVi7y9XOhr6gfsiJvzRi6rjkEshVt_cxgnkrRTGP7uxv06OiLojl96UuwvMnHcvPAJmEdfLZs0fBwC48lNT10diDOVoo7ze4XFwPYSbhuOZPTVMoFMeToz5OPGsz-95B/s400/bacteri.png" width="400" /></a></div>
<br />Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-49705736943829606572013-04-24T08:27:00.000-07:002013-04-24T08:27:11.302-07:00Generating n-Dimensional Lattice CoordinatesIt comes up sometimes in my research that I need to generate coordinates along a lattice. For example, when generating an initial structure for a molecular dynamics simulation. I thought I share my code for this. The code works by recursively enumerating each dimension in a lattice and calling a given function each time it reaches a point. So, it could be used to execute any function along a lattice. Here's what it looks like for a 2D and 3D example. In the 3D example I didn't choose a cubic number of points on the lattice, so that it sort of stops in the middle.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgP_0ANbZx4171nlCQGyefsdyWsvnOy9jlaQQoBSGvAQkvNU3Npio5PfoTk6iIyCD7CYU6VFv0ZV7Vdqt5duJnYvU925pAaJHS-O_mLrsPcsIpkCfENMYB62vGbc2x-1-qxtOuhMbyH5vZC/s1600/2dlattice.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgP_0ANbZx4171nlCQGyefsdyWsvnOy9jlaQQoBSGvAQkvNU3Npio5PfoTk6iIyCD7CYU6VFv0ZV7Vdqt5duJnYvU925pAaJHS-O_mLrsPcsIpkCfENMYB62vGbc2x-1-qxtOuhMbyH5vZC/s320/2dlattice.png" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjEjZWy7ITLkInFaGDqSxcN38_XO_rEMpaUGj-AOxWHRb3nkfuZQgwdUP9eqe4kFNM6JWvMc6flyUtyeVlYBtFoPCvaVkcSn4q6nWaEPNzkoK3xtyOVW5Ofc_fuiGk1wmApP6UYLDOxf6-f/s1600/3dlattice.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="200" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjEjZWy7ITLkInFaGDqSxcN38_XO_rEMpaUGj-AOxWHRb3nkfuZQgwdUP9eqe4kFNM6JWvMc6flyUtyeVlYBtFoPCvaVkcSn4q6nWaEPNzkoK3xtyOVW5Ofc_fuiGk1wmApP6UYLDOxf6-f/s320/3dlattice.png" width="320" /></a></div>
<br />
<a name='more'></a><br />
Here's the code:
<pre>
#!/usr/bin/python
import sys, copy
from math import *
def print_usage():
print "Usage: build_lattice.py [dims] [box_1_size] [box_..._size] [box_dims_size] [number]"
if(len(sys.argv) < 2):
print_usage()
exit()
n_dims = int(sys.argv[1])
if(len(sys.argv) < 2 + n_dims + 1):
print_usage()
exit()
box_size = [float(sys.argv[x]) for x in range(2,2+n_dims)]
number = int(sys.argv[2 + n_dims])
volume = 1
for x in box_size:
volume = volume * x
increment = (volume / number) ** (1. / n_dims)
def prepend_emit(array, element):
array_copy = copy.copy(array)
array_copy.insert(0, element)
return array_copy
def enumerate_grid(fxn, dim, sizes, indices):
if(dim > 0):
for i in range(sizes[dim]):
enumerate_grid(fxn,
dim - 1,
sizes,
prepend_emit(indices, i))
else:
for i in range(sizes[dim]):
fxn(prepend_emit(indices, i))
def print_point(point):
print point,
def print_grid(indices):
global number
#check if to make sure we still need to make particles
if(number > 0):
number -= 1
global increment
[print_point(x * increment) for x in indices]
print ""
enumerate_grid(print_grid, n_dims - 1,
[int(ceil(x / increment)) for x in box_size],
[])
</pre>Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-17697176960816761232013-04-23T20:24:00.000-07:002013-05-13T19:15:40.883-07:00Category Colors in RI thought I'd write down the color palette I use for categories. I've seen many discussions about gradient color choices, but category colors are less often discussed. I have a function given below which generates colors where there is not meant to be an ordering to the color. It is important when data points are near, that they are easy to distinguish with the eye. It makes more sense to see it:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjAxsecuTE1KS6yTt0Sr3n5PYuK8Mo3B2Wp_F04_MqcSxoZ1CPPjyajm5ZhF5saH2xuSFpa-zOQIY5JqvaLFc6TAe2pUa3_Fr9L_zjpSNquf4Nmy6O8KfgfA4Khfpjek0Wf0vagPOORr2m4/s1600/color.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="400" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjAxsecuTE1KS6yTt0Sr3n5PYuK8Mo3B2Wp_F04_MqcSxoZ1CPPjyajm5ZhF5saH2xuSFpa-zOQIY5JqvaLFc6TAe2pUa3_Fr9L_zjpSNquf4Nmy6O8KfgfA4Khfpjek0Wf0vagPOORr2m4/s400/color.png" width="200" /></a></div>
<br />
<a name='more'></a>Here's the function:<br />
<br />
<script src="https://gist.github.com/whitead/9149f62bbbab7a3f8fc9.js"></script>
<br />
Here's a comparison of the function (top) and the built-in rainbow function (bottom) which serves the same purpose.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgypbjZ9C0lyWnqKM6ZwQbhp0iRCRkj-QXy5XDIDLV2nT4m1qYqMVmhdrIXQvSyZHaioyfU5snDm_Pnuj641BQEofd6cnlEyScKp_Rv-z7aXxhczrAVhkU7evj9dVJY7GtcBOPC3DHf0JVm/s1600/example.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgypbjZ9C0lyWnqKM6ZwQbhp0iRCRkj-QXy5XDIDLV2nT4m1qYqMVmhdrIXQvSyZHaioyfU5snDm_Pnuj641BQEofd6cnlEyScKp_Rv-z7aXxhczrAVhkU7evj9dVJY7GtcBOPC3DHf0JVm/s1600/example.png" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhmG-Xh0fmxKyA3AbraxIRdPxVXMhkjuwFoA702Cit_9jbPrmqBhzOo_0DXWO7LMZ6dbSvBiHBkcN9-dj7jGldyAE5XbaYzg77OqANkzbVcyE79tABgXi5a6fvoM_nkF1PWIttt00Kz95Wq/s1600/example_2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhmG-Xh0fmxKyA3AbraxIRdPxVXMhkjuwFoA702Cit_9jbPrmqBhzOo_0DXWO7LMZ6dbSvBiHBkcN9-dj7jGldyAE5XbaYzg77OqANkzbVcyE79tABgXi5a6fvoM_nkF1PWIttt00Kz95Wq/s1600/example_2.png" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com1tag:blogger.com,1999:blog-4787751631116862322.post-20169334469849520252013-04-22T20:40:00.001-07:002013-05-13T19:14:56.451-07:00ggplot Style in Classic R Plots<a href="http://docs.ggplot2.org/current/index.html">ggplot</a> is a highly acclaimed R package for plotting. I have had very little experience with the library because I've mostly memorized all the quirks of normal R plots. I like the ggplot default style though, so I thought write down how to replicate it in R plots. The distinctive feature of ggplot is the gray rectangle and white grid lines. This may be replicated like so:<br />
<br />
<br />
<script src="https://gist.github.com/whitead/5eda03e94ba41205bfd2.js"></script>
<br />
Here's an example<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjeV7Iejex1mvXgqmWQ64QWmGG8CBtTx_3B4SJKaPQani-nk0YRndOyqCprVZuGHfQl5LLRbgrp-mpM16Jumy2Q8x5qIPi7hJUGoBNpMgETM5cHUouXYc_z0VS4MwWMHEc0XJnIk5T0LHYs/s1600/ggplot.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="180" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjeV7Iejex1mvXgqmWQ64QWmGG8CBtTx_3B4SJKaPQani-nk0YRndOyqCprVZuGHfQl5LLRbgrp-mpM16Jumy2Q8x5qIPi7hJUGoBNpMgETM5cHUouXYc_z0VS4MwWMHEc0XJnIk5T0LHYs/s320/ggplot.png" width="320" /></a></div>
<br />
<a name='more'></a><br /><br />
Here's the code to generate that plot:
<br />
<script src="https://gist.github.com/whitead/7a4d4359337b8dc40782.js"></script>
<br />
If you ever dislike that extra space that extends the plot region, it may be removed by adding <code>xaxs="i"</code> to <code>par</code> or <code>plot</code>.Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-40611938445055932812013-04-20T15:37:00.003-07:002013-04-20T15:37:52.349-07:00DNA-Polymer NanoparticleDNA-Polymer nanoparticle are used to deliver DNA to cell nuclei for curing diseases at the genetic level. In my research group, we work with <a href="http://www.sciencedirect.com/science/article/pii/S0142961210001420">carboxy-betaine acrylamide polymers</a>. I modeled an exceptionally large DNA-polymer nanoparticle below. The polymer is 80 kDa and the DNA is a single strand of 65 base pairs. After condensing it (in <a href="http://www.blender.org/">blender</a>, not a true physical simulation), the nanoparticle comes out to about 15nm in diameter. <div>
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi4M2dWQYC6C4zIkbnrLSx5AxphYaR0b2alfmczPogMufXyhCb7qYcxWpIhVWG0GIHfnLhAluPxTsyXhzCRb3XFKlVm2Y519JtVCoSDBwtUghB2JTFKVp6q8CuVZv2-Q2e33fP0P0lNZ-bg/s1600/dna_polymer.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="225" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi4M2dWQYC6C4zIkbnrLSx5AxphYaR0b2alfmczPogMufXyhCb7qYcxWpIhVWG0GIHfnLhAluPxTsyXhzCRb3XFKlVm2Y519JtVCoSDBwtUghB2JTFKVp6q8CuVZv2-Q2e33fP0P0lNZ-bg/s400/dna_polymer.png" width="400" /></a></div>
<div>
<br /></div>
<div>
<br /></div>
Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-74046378561088377082013-04-19T15:31:00.000-07:002013-04-19T15:31:54.955-07:00The Accuracy of Ghose-Crippen ALogP for PeptidesI was using the <a href="http://onlinelibrary.wiley.com/doi/10.1002/jcc.540070419/abstract">Ghose-Crippen ALogP predictor</a> implemented in <a href="http://sourceforge.net/projects/cdk/">CDK</a> as a QSAR for looking at peptide solubilities. Unfortunately, this partition coefficient predictor is fit on small organic molecules and does a pretty bad job at predicting peptide partition coefficients. Here's what it looks like on single amino acids compared with<a href="http://pubs.rsc.org/en/Content/ArticleLanding/1992/P2/p29920000079"> experimental data</a>.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEheipKwh3NnI9SfAhe-_fiDgIfHGfA1LLlkJqu3JzoSBlk14NJ2nAa1BLUnFp6yJjL5tPVmP7uBEvivnRsSJNpFYlcgJHrMUClUMbzyhQNjOkv22STH1XD59wL31nmLrvESzodY1yRFkIOV/s1600/comparison.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEheipKwh3NnI9SfAhe-_fiDgIfHGfA1LLlkJqu3JzoSBlk14NJ2nAa1BLUnFp6yJjL5tPVmP7uBEvivnRsSJNpFYlcgJHrMUClUMbzyhQNjOkv22STH1XD59wL31nmLrvESzodY1yRFkIOV/s400/comparison.png" width="400" /></a></div>
<br />
<a name='more'></a><br />
That's not too bad, all things considered. But, it gets much worse as the peptides get larger. Here's a histogram of the partition coefficients of a collection of 5,000 length 10 peptides.<br />
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhOnotxJaBzk4BXee6H6oS7japv-lZApIQD7V3-p-ETO8utgAz4MP0Q8zK0cSsEjDkKBG0LfEIhe8epM9CgRkeNkZRshUyGsgEpjdszuDpMxLupFuKThw-t6vkBjDz9LZpvrU3khYa2N401/s1600/hist.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhOnotxJaBzk4BXee6H6oS7japv-lZApIQD7V3-p-ETO8utgAz4MP0Q8zK0cSsEjDkKBG0LfEIhe8epM9CgRkeNkZRshUyGsgEpjdszuDpMxLupFuKThw-t6vkBjDz9LZpvrU3khYa2N401/s1600/hist.png" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
As you can see, the partition coefficients are way too low. It turns out the partition coefficient predictions are correlated with the peptide length, but as someone who works with peptides, I know that their solubility is not a really a function of their length. Part of the problem for this is that the surface of these peptides is a linear function of weight because I'm not attempting to find the correct 3D geometry of the peptides.<br />
<br />
How do I deal with this problem? Well, it has to do with my post on<a href="http://crowsandcats.blogspot.com/2013/04/approximate-restricted-integer.html"> approximate integer partitions</a>. I'm not sure if I'll have time to write it up; it'll definitely be in my next publication though! The basic idea is not to look at the absolute partition coefficient, but to compare them with similarly sized peptides. I can do that, because I'm using the predicted partition coefficient as a QSAR and am not actually trying to predict it.Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-5642597418628056272013-04-18T22:03:00.002-07:002013-04-18T22:03:53.722-07:00Visualizing Size in Biochemical Systems - UpdateAfter some helpful suggestions, I've revised the <a href="http://crowsandcats.blogspot.com/2013/04/visualizing-size-in-biochemical-systems.html">previous image I made</a>. This time I've used the atom representation for the proteins and peptides, so atoms are used everywhere for visualization and they have the correct radii. Doesn't look as cool, but it is more consistent.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjh0cnhCH1vAGZslrRZ_Uz31jh1CG4sYRkp95o6EOG4-3hMdM0ODz_vA4iG9a6Db-G1KMRpOh7zGcKnOhhM7RucHaY7e1y0UJN_5qY5QOaVXAIEpoMGo1JBandBZ6CMud_OgFYkoAApr8Tl/s1600/full_annotation_2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgBdZq49JhuB8qpqeqKIv7XJiy9ACPPRN37l_apY6y2ReC2ECQ1hHtOHdOXvTffcTM6iyNtzIlhRa6Zwds4KH1xfurBom0rXnCxRcUCuXHYnTmq5jmKrqz2_4fxHDpy6q5EhRt5q8nkKvfj/s1600/full_annotation_2_thumb.png" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi5AJwIC-YYfhWp6vBiLAE30OJoncF3M2pARN-uE-fJMSFveNOCfVFXZSE4MNdYnBJ4LcGHcQaLavu1A2o_6R7eEpymKbJkwIoyxirgYO3NlLQjsBKvVByCyyVKhgveRkyP0ItO8xfbD8T_/s1600/zoom_2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="225" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi5AJwIC-YYfhWp6vBiLAE30OJoncF3M2pARN-uE-fJMSFveNOCfVFXZSE4MNdYnBJ4LcGHcQaLavu1A2o_6R7eEpymKbJkwIoyxirgYO3NlLQjsBKvVByCyyVKhgveRkyP0ItO8xfbD8T_/s400/zoom_2.png" width="400" /></a></div>
<br />Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com2tag:blogger.com,1999:blog-4787751631116862322.post-15278648840668690742013-04-15T20:58:00.002-07:002013-04-16T13:48:21.787-07:00Visualizing Size in Biochemical SystemsI decided to waste a few hours and make a <a href="http://www.blender.org/">Blender</a> image based on my post on the <a href="http://crowsandcats.blogspot.com/2013/04/relative-size-of-proteins-nanoparticles.html">relative size of common elements in biochemical systems</a>. Using 1 blender unit = 1 nm, I modeled a gold nanoparticle, an antibody, albumin, an iron oxide nanoparticle, a 10 amino acid length peptide, and salicylic acid.<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh3zAPGfo6UkqcQTbyi5pG0ApEXWYWzR5pH1MAyANECuXPBWWp_SmtxcEePSGLJH69YW0WyOJtXd4pFQ51lS96FcgjjsM5YmZ2gG7C3ebfuVAkTnAIzCVBfVKYEgRlXjDcAskH715WkgUId/s1600/full_annotation.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh3zAPGfo6UkqcQTbyi5pG0ApEXWYWzR5pH1MAyANECuXPBWWp_SmtxcEePSGLJH69YW0WyOJtXd4pFQ51lS96FcgjjsM5YmZ2gG7C3ebfuVAkTnAIzCVBfVKYEgRlXjDcAskH715WkgUId/s400/full_annotation_thumb.png"/></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<br />
<br />
<a name='more'></a><br />
<br />
The nanoparticles were created using dupliverts on Icospheres in Blender. Albumin and the Antibody were created by rendering them as STL files in <a href="http://www.ks.uiuc.edu/Research/vmd/">VMD</a> using the surf representation. Once imported into Blender, a decimate modifier followed by remesh (octree depth of 7, smooth) fixed the STL meshes nicely. A 5 iteration smooth modifier fixed any artifacts after the remeshing. Make sure to turn on smooth shading. The peptide and salicylic acid were made the same way, except using the VDW representation instead of the surf. No remeshing is necessary, other than enabling smooth shading. Rendering was done at 150% with 75 samples and the image was scaled down in GIMP to 100% using cubic interpolation (a ghetto denoising technique that works well). Finally, the annotations were added in Inkscape. Here's a zoom of the smaller elements.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEinu1VXJi3dNQvY-yv2b1bMB9YQN6ih5rn09R1YVMmhYVtStRDAx6dYvWFUyWMf6qM6_IGTEot6dy24KlKggWkoc0l1QIwKx5bShI0F4Bx-G8CNX3aO0oLpkaol2C2JCg-eOH3ZX1ZEkwSC/s1600/zoom.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="241" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEinu1VXJi3dNQvY-yv2b1bMB9YQN6ih5rn09R1YVMmhYVtStRDAx6dYvWFUyWMf6qM6_IGTEot6dy24KlKggWkoc0l1QIwKx5bShI0F4Bx-G8CNX3aO0oLpkaol2C2JCg-eOH3ZX1ZEkwSC/s400/zoom.png" width="400" /></a></div>
<br />Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com4tag:blogger.com,1999:blog-4787751631116862322.post-29686350938377702282013-04-15T13:23:00.003-07:002013-04-15T16:14:19.183-07:00Relative Size of Proteins, Nanoparticles, Antibodies, Peptides, and DrugsNeeding to know the relative size of these elements in biochemical systems often comes up in my research, so I've decided to write them.<br />
<br />
<br />
<style>
#chemistry {
border-collapse:collapse;
}
#chemistry tr.alt td
{
color:#000000;
}
#chemistry th
{
font-size:1.1em;
text-align:left;
padding-top:5px;
padding-bottom:4px;
background-color:#21bb53;
color:#ffffff;
}
#chemistry td, #chemistry th
{
font-size:1em;
border:1px solid #21bb53;
padding:3px 7px 2px 7px;
}
</style>
<table id="chemistry">
<tr>
<th> Component </th>
<th> Size [nm] </th>
<th> Reference</th>
</tr>
<tr>
<td> Gold Nanoparticle</td>
<td> ~20</td>
<td>Can vary, most common is 20. <a href="http://www.nature.com/nature-physci/journal/v241/n105/abs/physci241020a0.html"> Frens (1973)</a> </td>
</tr>
<tr>
<td> Iron Oxide Nanoparticle</td>
<td> 5-15</td>
<td>Can be much larger with emulsion methods. This size is most common when magentic properties are desired. Called <a href="http://www.sciencedirect.com/science/article/pii/S0927775704003206">partial reduction coprecipitation method.</a></td>
</tr>
<tr>
<td> Antibody</td>
<td> 15.5 x 11 </td>
<td> <a href="http://pdb.org/pdb/explore/explore.do?structureId=1igt"> PDB: 1igt</a> measured in <a href="http://www.ks.uiuc.edu/Research/vmd/">VMD</a></td>
</tr>
<tr>
<td> Human Serum Albumin</td>
<td> 7.6 x 7 </td>
<td> <a href="http://pdb.org/pdb/explore/explore.do?structureId=1e7i"> PDB: 1e7i</a> measured in <a href="http://www.ks.uiuc.edu/Research/vmd/">VMD</a></td>
</tr>
<tr>
<td> 10 Amino Acid Peptide</td>
<td> ~1.5 (coiled globule sphere)</td>
<td> <a href="http://pubs.acs.org/doi/abs/10.1021/ja3006868"> Paper I wrote on Peptide Conformations</a></td>
</tr>
<tr>
<td> Salicylic Acid (drug) </td>
<td> 0.6 x 0.5</td>
<td> <a href="http://pubchem.ncbi.nlm.nih.gov/summary/summary.cgi?cid=338"> PubChem structure</a> measured in <a href="http://www.ks.uiuc.edu/Research/vmd/">VMD</a></td>
</tr>
</table>Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com1tag:blogger.com,1999:blog-4787751631116862322.post-7824466089946538092013-04-12T08:55:00.002-07:002013-04-12T09:01:22.223-07:00Approximate Restricted Integer Partition With Exact Part NumberI'm interested in a certain type of mathematical structure this week: integer partitions. An integer partition, <code> p(n)</code>, is the set of sets of integers that add up to <code>n</code>. For example,<br />
<br />
<ul>
<li><code>p(3) = 1 + 1 + 1, 1 + 2, 3</code></li>
<li><code>p(4) = 1 + 1 + 1 + 1, 1 + 1 + 2, 2 + 2, ...</code></li>
</ul>
<br />
<a name='more'></a><br />
<br />
Each grouping of integers that add up to <code>n</code> is a partition. I'm interested in partitions where the number of integers is fixed, the order is fixed, and the integers are bounded. This makes more sense if you view it as a problem of putting numbers into boxes. Let's say you have four boxes: <u> </u> <u> </u> <u> </u> <u> </u>. You want to partition 8. An example would be then: <u>2</u> <u>3</u> <u>1</u> <u>2</u>. Or perhaps <u>2</u> <u>2</u> <u>2</u> <u>2</u>. <br />
<br />
The next restriction is on the bounds of each number. With the previous example, let's say I now restrict each integer to be less than 4. So <u>3</u> <u>3</u> <u>0</u> <u>2</u> would work. I'll call this function <code>b(n, m, k)</code> where <code>n</code> is the sum to partition, <code>m</code> is the number of parts, and <code>k</code> is the maximum at each position. The range of <code>b</code> will 0 through infinity because I'm going to have it count the number of partitions, not be the set of partitions. <br />
<br />
An algorithm to calculate <code>b(n, m, k)</code> is given below in python. It works exactly like you would expect: enumerating the exponential number of partitions. The only trick it does is know when a partition it is trying will fail. For example, if you're trying to do the previous example and you have <u>1</u> <u>0</u> <u> </u> <u> </u>, there are no two other numbers that will cause the partition to sum to 8 (<code>1 + 0 + 3 + 3 = 7</code>). That helps reduce the number of partitions it attempts. Here's the code:<br />
<br />
<pre>#Recursive restricted block partition. Each call reduces the sum, n,
# and m, the number of parts.
def blockpart(n, m, k):
#Is the number of parts only 1?
if(m == 1):
#Trivial: to partition n into one part, we need to be able to use n
if(n < k):
return 1
else:
return 0
#what is the highest possible sum we may have
#that will still be a valid partition
smax = m * (k - 1)
#This is the lowest term we can add to the partition (subtract from sum)
start = n - k + 1
#If the lowest is too high, then we aren't able to get a valid partition
if( start > smax):
return 0
paths = 0
#add up all possible other partitions after subtracting choices
#to place in the current part (m)
for i in range(max(0,start), 1 + min(smax, n)):
paths += blockpart(i, m - 1, k)
return paths
</pre>
<br />
<br />
Anyway, it turns out that as <code>n</code> grows, it may be approximated by the equation:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiJ8npo09OhH0kwQRc-hIUdYfqKuZtGcpltWdNcinKnEU57XG82n5D0F4HGH4VKstvT_tRZbyhjhto4pJfqQQ-sZVgBCVdHvN3VW7L78a7VhyphenhyphenDj1VAuwmePRnIpuFg_p85Qt8-yTXySo_LH/s1600/1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiJ8npo09OhH0kwQRc-hIUdYfqKuZtGcpltWdNcinKnEU57XG82n5D0F4HGH4VKstvT_tRZbyhjhto4pJfqQQ-sZVgBCVdHvN3VW7L78a7VhyphenhyphenDj1VAuwmePRnIpuFg_p85Qt8-yTXySo_LH/s1600/1.png" /></a></div>
<br />
Who discovered that equation? Me. And I'm sure many other people. Here's the fit of the equation, where the dots are the results of the python program and the red lines are the approximation.<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjmasWl-spx10qy_H7C2udDbtRvSzk14cu_p1AzSvhySWHr_elSbiXaaBVmP0ykT6XW67h6OtRnW8rcmAfjxMaCg81ORBxsxCyNLVcU-OHkgRNaF-1bOeJFfTprB6a3naHoPqK9lVAHt3hi/s1600/plot_2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjmasWl-spx10qy_H7C2udDbtRvSzk14cu_p1AzSvhySWHr_elSbiXaaBVmP0ykT6XW67h6OtRnW8rcmAfjxMaCg81ORBxsxCyNLVcU-OHkgRNaF-1bOeJFfTprB6a3naHoPqK9lVAHt3hi/s1600/plot_2.png" /></a></div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjifphr_B3VZR-Wa0BXkvd2qZHxnn3_lG_yDn3YRBAEBRKbuSjGjAGRD0bVVxwEItSETRAW2sb3SGW2lEhdeS0Sl9oP49AaaS0o4SHhzWLD5wX6hOewt0b6dXyHg6Hkrh4NLjNwX34HLPz6/s1600/plot_20.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjifphr_B3VZR-Wa0BXkvd2qZHxnn3_lG_yDn3YRBAEBRKbuSjGjAGRD0bVVxwEItSETRAW2sb3SGW2lEhdeS0Sl9oP49AaaS0o4SHhzWLD5wX6hOewt0b6dXyHg6Hkrh4NLjNwX34HLPz6/s1600/plot_20.png" /></a></div>
<br />
<br />
Find out next time on Crows and Cats how I derived that equation!Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com2tag:blogger.com,1999:blog-4787751631116862322.post-48444205875086772872013-04-11T14:58:00.000-07:002013-05-13T19:13:09.280-07:00Speeding up RI've been programming in R for a while now. Whenever I've had performance problems, it almost always is due to <code>data.frame</code> usage. To check what's slowing down your R code, just use the <code>Rprof</code> command like so:<br/>
<br/>
<script src="https://gist.github.com/whitead/963c4408fd1b77611f58.js"></script><a name='more'></a><br/>
I wrote this code to be particularly slow due to its <code>data.frame</code> usage. You may view the results of the profiling by running <code>R CMD Rprof summ.prof</code>:
<pre>
Each sample represents 0.02 seconds.
Total run time: 2.98 seconds.
Total seconds: time spent in function and callees.
Self seconds: time spent in function alone.
% total % self
total seconds self seconds name
81.2 2.42 1.3 0.04 "[<-"
79.9 2.38 65.8 1.96 "[<-.data.frame"
18.1 0.54 10.1 0.30 "[.data.frame"
18.1 0.54 0.0 0.00 "["
12.1 0.36 1.3 0.04 "%in%"
11.4 0.34 9.4 0.28 "match"
4.7 0.14 4.0 0.12 "anyDuplicated"
2.0 0.06 2.0 0.06 "names"
2.0 0.06 2.0 0.06 "sys.call"
1.3 0.04 1.3 0.04 "=="
0.7 0.02 0.7 0.02 ".row_names_info"
0.7 0.02 0.7 0.02 "NROW"
0.7 0.02 0.7 0.02 "anyDuplicated.default"
0.7 0.02 0.7 0.02 "cos"
% self % total
self seconds total seconds name
65.8 1.96 79.9 2.38 "[<-.data.frame"
10.1 0.30 18.1 0.54 "[.data.frame"
9.4 0.28 11.4 0.34 "match"
4.0 0.12 4.7 0.14 "anyDuplicated"
2.0 0.06 2.0 0.06 "names"
2.0 0.06 2.0 0.06 "sys.call"
1.3 0.04 81.2 2.42 "[<-"
1.3 0.04 12.1 0.36 "%in%"
1.3 0.04 1.3 0.04 "=="
0.7 0.02 0.7 0.02 ".row_names_info"
0.7 0.02 0.7 0.02 "NROW"
0.7 0.02 0.7 0.02 "anyDuplicated.default"
0.7 0.02 0.7 0.02 "cos"
</pre>
As you can see, it's very slow. After profiling your own code, if you find that the top calls are <code>[.data.frame</code> or <code> [<-.data.frame</code>, then you have a <code>data.frame</code> problem. Here's how I solve this, in order of things I try:<br/>
<br />
<ol>
<li> avoid loops, use vectorized code (no <code>for</code> loops, no <code>apply</code>, no <code>sapply</code>). In the example, use <code>d[,1]<-</code>, for assigning an entire column</li>
<li> Use numeric indices when possible. In our example, that means using <code>d[i, 1]</code> instead of <code>d[i, "x"]</code></li>
<li> Get rid of the <code>data.frame</code> for heavy calculations by using the <code>data.matrix</code> command. In the example above, just use <code> d.matrix <- data.matrix(d) </code></li>
</ol>
<br /><br/>
Rewriting the silly example above, <br/><br/>
<script src="https://gist.github.com/whitead/8be72442211db64f5ff5.js"></script>
<br/> This time, it runs too quickly for any events to be recorded with the profiling. The execution time is 0.192 seconds, whereas the first version is 3.48 seconds. A pretty good speed-up.Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-73095251135918628822013-04-05T15:15:00.000-07:002013-04-11T21:34:24.507-07:00Changing Image Size in Tachyon in VMDI always forget how to change the image size for rending molecular images in <a href="http://www.ks.uiuc.edu/Research/vmd/">VMD</a> using the (external) Tachyon rendering engine, so I thought I'd write down the steps here:<br />
<br />
In the TCL console, get the aspect ratio correct for the image. So, for example to make a nice <code>1920x1080</code> image, we'll make the view at half of that:<br />
<br />
<pre class="brush: plain; toolbar: false;">
display resize 960 540
</pre>
<br />
Now, add simply add <code>-res 1920 1080</code> to the end of the render command in the file render controls window. Here's an example of one of my render commands:<br />
<pre class="brush: plain; toolbar: false;">
"$HOME/VMD/tachyon_LINUXAMD64" -aasamples 12 %s -format TARGA -o %s.tga -res 1920 1080
</pre>
<br />
Here's the result of one my renders:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjjKmucn-EYOAy41uvi5GcW5q0jlie3AtYMSR0suvGXZ4oqLnJMszciGU-u33dyfeaunH2kYLNAk8W_FjdV-Mv10npCj6H6pqkkq1iUGuSxvSLWcXKzesPzoKMSdtyGOth1KWR1h9046PAP/s1600/quat_little.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="163" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjjKmucn-EYOAy41uvi5GcW5q0jlie3AtYMSR0suvGXZ4oqLnJMszciGU-u33dyfeaunH2kYLNAk8W_FjdV-Mv10npCj6H6pqkkq1iUGuSxvSLWcXKzesPzoKMSdtyGOth1KWR1h9046PAP/s400/quat_little.png" width="400" /></a></div>
Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com4tag:blogger.com,1999:blog-4787751631116862322.post-44849397203409279492013-04-03T15:56:00.002-07:002013-04-05T15:16:09.325-07:00Making nanoparticles in BlenderToday I made some nanoparticles in Blender. I used the node compositor to create the fluorophore effect. Making science in Blender is fun!<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg7deRwG3q6znDSNXXKtds2SMQg1spU3m0AcWNYBWyyAJlQBqiefRY3NtUVBB_vPTnACIggnC3uphxMt2eJGQAJOzIyfxtoTpyYZnlAjthDbB7cmZn3cuaLOwDyuMtQ7WBH7Sz6ybRQLLHj/s1600/nanocomposite_2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="241" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg7deRwG3q6znDSNXXKtds2SMQg1spU3m0AcWNYBWyyAJlQBqiefRY3NtUVBB_vPTnACIggnC3uphxMt2eJGQAJOzIyfxtoTpyYZnlAjthDbB7cmZn3cuaLOwDyuMtQ7WBH7Sz6ybRQLLHj/s400/nanocomposite_2.png" width="400" /></a></div>
<br />Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com9tag:blogger.com,1999:blog-4787751631116862322.post-64918216159861445062013-03-23T22:34:00.001-07:002013-03-23T22:34:28.422-07:00Making Proteins and Polymer in BlenderI was messing around in <a href="http://www.blender.org/">Blender</a> and decided to create some artwork related to a <a href="http://pubs.acs.org/doi/abs/10.1021/bm301335r">recent paper</a> published in my current <a href="http://depts.washington.edu/jgroup/">research group</a>. The idea is to create a two layers of polymer. One layer is to make a biosensor resistant to noise from non-target molecules adsorbing on the sensor surface. A second longer layer of polymer is used to attach anti-body to the sensor. The anti-body allows the sensor to target very specific molecules. Here's the image I made using Blender's cycles render engine:<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgnXz7MU47n7QxqnNQG1Vs3zDHGbv0meAcyAVVMU0H9-lyVzXyJo1FVwkm59GOoUSVXGbxE-Q-KdVKDE8TbNGxlrVk64YXp6Q_kd12MUD9gX0d2KB8CLU7rQMg0haZ2i25cz3xqZMOzkE0V/s1600/bidisperse.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="225" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgnXz7MU47n7QxqnNQG1Vs3zDHGbv0meAcyAVVMU0H9-lyVzXyJo1FVwkm59GOoUSVXGbxE-Q-KdVKDE8TbNGxlrVk64YXp6Q_kd12MUD9gX0d2KB8CLU7rQMg0haZ2i25cz3xqZMOzkE0V/s400/bidisperse.png" width="400" /></a></div>
<br />Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-69552158883922981942013-03-20T16:40:00.001-07:002013-05-13T19:10:54.389-07:00CDK QSARs of PeptidesI was doing some working in <a href="http://sourceforge.net/projects/cdk/">CDK</a> and was pleasantly surprised at how easy it is to calculate descriptors on short peptides. The code to do this is actually only three lines! Adding in a bug fix though, it becomes four. Here's the code to calculate every QSAR descriptor on a peptide
<br />
<br />
<script src="https://gist.github.com/whitead/fa8f5f54e469cdc60b7b.js"></script>
<br />
<a name='more'></a>The complete code with imports and outputting the results is shown below:
<br />
<script src="https://gist.github.com/whitead/fffd4dfcf19b76fc1e53.js"></script>Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-39834389609853669072013-03-11T10:13:00.001-07:002013-04-09T13:52:20.606-07:00Reading Command Line args in RI'm writing this down in my blog because I keep forgetting my standard way for adding command line arguments to R. I like to use executable R scripts that don't require <code>--vanilla --args</code> etc. The syntax in the R script is:<br />
<br />
<a name='more'></a><br />
<br />
<code>script.R:</code>
<br />
<pre class="brush: plain; toolbar: false;">
#!/usr/bin/env Rscript
argv <- commandArgs(trailingOnly=TRUE)
</pre>
<br/>
Make it executable:
<pre class="brush: shell; toolbar: false;">chmod +x script.R
</pre>
Execute with
<pre class="brush: shell; toolbar: false;">./script.R extra arguments go here
</pre><br/>
I prefer this method much more than others because the R script acts like any other command line tool now.
Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-28114969535167177982013-02-28T06:57:00.000-08:002013-08-21T18:02:46.749-07:00Markov Chains in JS — Part 1I'm in the process of writing a new post or two on working with Markov Chains in javascrtipt. I created the javascript files already and I've decided to post them below for now. I'll write up more about it and hopefully use the code for something interesting in a future post. See below for the javascript Markov Chain Simulator! See the code <a href="http://jsbin.com/azurir/22/edit">here.</a><br />
<br />
<a name='more'></a><br />
Well, it turns out that I can't find hosted versions of the libraries I used in creating this. So I guess that's the end of this project.
<br />
<br />
<div style="overflow: hidden;">
<object data="http://jsbin.com/useyav/81/quiet" style="height: 850px; overflow: hidden; width: 520px;" type="text/html"> </object>
</div>
Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0tag:blogger.com,1999:blog-4787751631116862322.post-31599996575331729742012-12-29T08:21:00.000-08:002013-05-13T19:05:26.198-07:00Recovering Android Photos with LinuxAs noticed by others on the Internet, it is possible to recover photos from an Android phone that were deleted. It is impossible, as far as I know, to recover the full resolution photos. You may, however, recover the thumbnails and those are often large. This also allows you to recover photos which weren't saved, like those from text messages. The problem though is that methods I found while googling require you to edit the image files in a hex editor. For example, as demonstrated in <a href="http://www.youtube.com/watch?v=0PH54OmYlqY">this video</a>. I've instead written a python script which automates this greatly simplifying the fact that you need to recover all the pictures to find one in particular.<br />
<br />
<br />
<a name='more'></a><br /><br />
<br />
The approach is to write a short python script which automates the hex editing necessary and then execute the script on all the cache files. Pretty straight forward. As usual, this only works on Linux and I have no idea about other OSes.<br />
<br />
The first step is to find your tec cache files. Plug the phone into a computer and turn on USB storage. Navigate around until you find any tec files, which are cached thumbnails. For me, the path is /sdcard/DCIM/Camera/cache<br />
<br />
Once you have the files on your computer in a folder (do not edit the originals on your phone!), you need to open a text editor to create the python script. This script in python is pretty simple:<br />
<br />
<script src="https://gist.github.com/whitead/b85d0cb344a1b39fe74e.js"></script>
Once you have that python script written in a text file, for example "tec2jpg.py", make it executable with
<br />
<pre class="brush: shell; toolbar: false;">chmod +x tec2jpg.py
</pre>
The final ingredient is to execute the python script on all the tec files in a directory. This can be done using this little bash for loop in the command line when you are in the directory containing the tec files:
<br />
<pre class="brush: shell; toolbar: false;">for i in `ls *tec`;
do NAME=`basename $i .tec`;
path/to/python-script/tec2jpg.py $NAME.tec $NAME.jpg;
done;
</pre>
That's it! You should at this point have jpegs for all the original tec cached thumbnails. I have only tested this on my Samsung galaxy, so I have no idea if it works on other phones. Best of luck!Andrew Whitehttp://www.blogger.com/profile/12233134401389651425noreply@blogger.com0