31:T736,lj_plot = { let marks = [ Plot.ruleX([1], { stroke: "black" }), Plot.ruleY([0], { stroke: "black" }) ]; if (lj_plot_inputs.includes("show schematics")) { marks = [ ...marks, [ Plot.image(spring_svgs, { x: "x", y: "y", src: "src", width: 200 }) ] ]; } if (lj_plot_inputs.includes("plot energy")) { marks = [ ...marks, [ Plot.line(lj_data, { x: "rho", y: "energy", stroke: "#4e79a7" //tip: true }), Plot.crosshairX(lj_data, { x: "rho", y: "energy", ruleStroke: "#4e79a7", ruleStrokeWidth: 2, ruleStrokeOpacity: 0.25 }), Plot.text(lj_data.slice(lj_data.length - 1), { x: "rho", y: "energy", text: ["Energy"], fill: "#4e79a7", dy: 10 }) ] ]; } if (lj_plot_inputs.includes("plot force")) { marks = [ ...marks, [ Plot.line(lj_data, { x: "rho", y: "force", stroke: "#e15759" //tip: true }), Plot.crosshairX(lj_data, { x: "rho", y: "force", ruleStroke: "#e15759", ruleStrokeWidth: 2, ruleStrokeOpacity: 0.25 }), Plot.text(lj_data.slice(lj_data.length - 1), { x: "rho", y: "force", text: ["Force"], fill: "#e15759", dy: 15 }) ] ]; } return Plot.plot({ width: 500, aspectRatio: 12.5, x: { domain: [0.75, 1.5], label: "interatomic distance" }, y: { domain: [-2.75, 2.75], grid: true, label: "potential" }, style: { fontSize: "12px" }, marks: marks }); }32:T78e,lj_harmonic_plot = { let marks = [ Plot.ruleX([1], { stroke: "black" }), Plot.ruleY([0], { stroke: "black" }) ]; if (lj_harmonic_plot_inputs.includes("show schematics")) { marks = [ ...marks, [ Plot.image(spring_svgs, { x: "x", y: "y", src: "src", width: 200 }) ] ]; } if (lj_harmonic_plot_inputs.includes("plot energy")) { marks = [ ...marks, [ Plot.line(lj_data, { x: "rho", y: "energy", stroke: "#4e79a7" //tip: true }), Plot.crosshairX(lj_data, { x: "rho", y: "energy", ruleStroke: "#4e79a7", ruleStrokeWidth: 2, ruleStrokeOpacity: 0.25 }), Plot.text(lj_data.slice(lj_data.length - 1), { x: "rho", y: "energy", text: ["Energy"], fill: "#4e79a7", dy: 10 }) ] ]; } if (lj_harmonic_plot_inputs.includes("plot harmonic approximation")) { marks = [ ...marks, [ Plot.line(lj_data, { x: "rho", y: "harmonic", stroke: "#e15759" //tip: true }), Plot.crosshairX(lj_data, { x: "rho", y: "harmonic", ruleStroke: "#e15759", ruleStrokeWidth: 2, ruleStrokeOpacity: 0.25 }), Plot.text(lj_data.slice(lj_data.length - 175, lj_data.length - 174), { x: "rho", y: "harmonic", text: ["Harmonic"], fill: "#e15759", dx: 35 }) ] ]; } return Plot.plot({ width: 500, aspectRatio: 12.5, style: { fontSize: "12px" }, x: { domain: [0.75, 1.5], label: "interatomic distance" }, y: { domain: [-2.75, 2.75], grid: true, label: "potential" }, marks: marks }); }33:T576,harmonic_animation={ // init plot + get relevant data and keyframe const svg = d3.select(DOM.svg(chart_param.width, chart_param.width/8)), nodes = one_dimensional_lattice.map(d => Object.create(d)), edges = one_dimensional_lattice_links.map(d => Object.create(d)), kf = restart_button_1d; let counter = 0; svg.attr('viewBox',[-2.5, -2.5, num_of_atoms + 5, 5]); // create simulation using: // linking force with user-defined distance and collide + charge forces to seperate nodes const sim = d3.forceSimulation(nodes) .force("link", d3.forceLink(edges).distance(1).iterations(0)) .alphaDecay(0) .velocityDecay(0) .dt(0.5); // draw nodes (with uniform radius and no edges) const [node, edge] = draw_nodes(svg, nodes, null, false, false); node .attr('opacity',d=>d.fx == undefined ? 1/2 : 1) .attr('fill',d=>d.fx == undefined ? 'gray' : 'black'); // update each node upon simulation tick sim.on("tick", () => { if (counter++ > 750 || kf ==0) sim.stop(); if (oscillation_direction == "horizontal"){ node .attr("cx", d => (d.x+(d.x-d.index)/oscillation_magnitude_static)) .attr("cy", d => d.y); } else { node .attr("cx", d => d.index) .attr("cy", d => (d.x-d.index)/oscillation_magnitude_static); } }); invalidation.then(() => sim.stop()); return svg.node(); }34:T4ea,function compute_single_period(num_of_atoms,sub_harmonic,oscillation_magnitude,dt,max_iterations) { if (sub_harmonic == 0 || Math.abs(sub_harmonic) == num_of_atoms) return Infinity; // zone-center //if (Math.abs(sub_harmonic) == num_of_atoms/2) return Infinity; // zone-edge if (Math.abs(sub_harmonic) == num_of_atoms/2) sub_harmonic+=1; // zone-edge if (oscillation_magnitude == undefined) oscillation_magnitude = oscillation_magnitude_static; if (dt == undefined) dt = 0.0075; if (max_iterations == undefined) max_iterations = 5000; let steps = [oscillation_magnitude],i=0,nodes_template,edges_template; [nodes_template,edges_template] = make_one_dimensional_lattice(num_of_atoms,sub_harmonic,oscillation_magnitude); const nodes = nodes_template.map(d => Object.create(d)), edges = edges_template.map(d => Object.create(d)); const simulation = d3.forceSimulation(nodes) .force("link", d3.forceLink(edges).distance(1).iterations(0)) .alphaDecay(0) .velocityDecay(0) .dt(dt); do{ i++; simulation.tick(); const copies = nodes.map(n => n.x-n.index); steps.push(copies.reduce((sum,n)=>sum+=n*n)); } while (i< max_iterations && (steps[i]-steps[i-1])<0) simulation.stop(); return (i-1)*dt }35:T535,dispersion_animation = { // init plot + get relevant data and keyframe const svg = d3.select(DOM.svg(chart_param.width, chart_param.width/8)); const [nodes, edges] = make_one_dimensional_lattice( num_of_atoms,dispersion_subharmonic,oscillation_magnitude_static ); let counter = 0; svg.attr('viewBox',[-2.5, -2.5, num_of_atoms + 5, 5]); // create simulation using: // linking force with user-defined distance and collide + charge forces to seperate nodes const sim = d3.forceSimulation(nodes) .force("link", d3.forceLink(edges).distance(1).iterations(0)) .alphaDecay(0) .velocityDecay(0) .dt(0.25); // draw nodes (with uniform radius and no edges) const [node, edge] = draw_nodes(svg, nodes, null, false, false); node .attr('opacity',d=>d.fx == undefined ? 1/2 : 1) .attr('fill',d=>d.fx == undefined ? 'gray' : 'black'); // update each node upon simulation tick sim.on("tick", () => { if (counter++ > 750) sim.stop(); if (oscillation_direction == "horizontal"){ node .attr("cx", d => (d.x+(d.x-d.index)/oscillation_magnitude)) .attr("cy", d => d.y); } else { node .attr("cx", d => d.index) .attr("cy", d => (d.x-d.index)/oscillation_magnitude); } }); invalidation.then(() => sim.stop()); return svg.node(); }36:T559,lennard_jones_animation = { // init plot + get relevant data and keyframe const svg = d3.select(DOM.svg(chart_param.width, chart_param.height/3)), kf=restart_button_1d_lj; const [nodes,edges] = make_one_dimensional_lattice(num_of_atoms,sub_harmonic,oscillation_magnitude) let counter = 0; svg.attr('viewBox',[-2.5,-2.5, num_of_atoms + 5, 5]); // create simulation using: // linking force with user-defined distance and collide + charge forces to seperate nodes const sim = d3.forceSimulation(nodes) .force("lennard-jones-BH",d3.forceLennardJonesPotential().theta(0.25)) .alphaDecay(0) .velocityDecay(0) .dt(0.025); // draw nodes (with uniform radius and no edges) const [node, edge] = draw_nodes(svg, nodes,null, false, false, colorScale1D); // update each node upon simulation tick sim.on("tick", () => { if (counter++ > 750 || kf == 0) sim.stop(); if (oscillation_direction == "horizontal"){ node .attr("cx", d => (d.x+(d.x-d.index)/oscillation_magnitude)) .attr("cy", d => d.y) .attr("fill",d=>colorScale1D(d.energy)); } else { node .attr("cx", d => d.index) .attr("cy", d => (d.x-d.index)/oscillation_magnitude) .attr("fill",d=>colorScale1D(d.energy)); } }); invalidation.then(() => sim.stop()); return svg.node(); }37:T4f7,surface_energy_slide_graphic = { // init plot + get relevant data and keyframe const svg = d3.select(DOM.svg(chart_param.width, chart_param.height)); const theta = theta_in_deg/180*Math.PI; const [nodes, nodesRight] = slice_particle(square_lattice,theta) svg.attr('viewBox',[-10, -10, 20, 20]); d3.forceSimulation(nodes) .force("lennard-jones-BH",d3.forceLennardJonesPotential().theta(0)) .alphaDecay(0) .velocityDecay(0) .dt(0) .tick() .stop(); d3.forceSimulation(nodesRight) .force("lennard-jones-BH",d3.forceLennardJonesPotential().theta(0)) .alphaDecay(0) .velocityDecay(0) .dt(0) .tick() .stop(); // draw nodes (with uniform radius and no edges) nodes.forEach(d => {d.x-=Math.sin(theta),d.y+=Math.cos(theta)}); nodesRight.forEach(d => {d.x+=Math.sin(theta),d.y-=Math.cos(theta)}); const [node, edge] = draw_nodes(svg, [...nodes,...nodesRight], null, false, false,colorScale); svg.append('line') .style("stroke", "black") .style("stroke-width", 0.05) .style("stroke-dasharray",0.25) .attr("x1", 12.5*Math.cos(theta)) .attr("y1", 12.5*Math.sin(theta)) .attr("x2", -12.5*Math.cos(theta)) .attr("y2", -12.5*Math.sin(theta)); return svg.node(); }38:T456,generalized_lj_animation = { // init plot + get relevant data and keyframe const svg = d3.select(DOM.svg(chart_param.width, chart_param.height)), nodes = hexagonal_lattice.map(d => Object.create(d)), kf = restart_button_hexagonal_generalized; svg.attr('viewBox',[-12.5, -12.5, 25, 25]); const sim = d3.forceSimulation(nodes) .force("lennard-jones-BH", d3.forceLennardJonesPotential() .theta(0.125) .distanceMax(1.75) .attractivePower(attractivePower) .repulsivePower(repulsivePower) ) .alphaDecay(0) .velocityDecay(0) .dt(0.05); sim.counter = 0; // draw nodes (with uniform radius and no edges) const [node, edge] = draw_nodes(svg, nodes, null, false, false,colorScale); node.call(drag(sim)); // update each node upon simulation tick sim.on("tick", () => { if (sim.counter++ > 250 || kf == 0) sim.stop(); node .attr("cx", d => d.x) .attr("cy", d => d.y) .attr("fill", d => colorScale(d.energy)); }); invalidation.then(() => sim.stop()); return svg.node(); }39:T443,morse_animation = { // init plot + get relevant data and keyframe const svg = d3.select(DOM.svg(chart_param.width, chart_param.height)), nodes = hexagonal_lattice.map(d => Object.create(d)), kf = restart_button_hexagonal_morse; svg.attr('viewBox',[-12.5, -12.5, 25, 25]); const sim = d3.forceSimulation(nodes) .force("morse-potential", d3.forceCentralPotential() .theta(0.125) .distanceMax(1.75) .potentialFunction(morsePotential) .gradientFunction(morseGradient) ) .alphaDecay(0) .velocityDecay(0) .dt(0.05); sim.counter = 0; // draw nodes (with uniform radius and no edges) const [node, edge] = draw_nodes(svg, nodes, null, false, false,colorScale); node.call(drag(sim)); // update each node upon simulation tick sim.on("tick", () => { if (sim.counter++ > 250 || kf == 0) sim.stop(); node .attr("cx", d => d.x) .attr("cy", d => d.y) .attr("fill", d => colorScale(d.energy)); }); invalidation.then(() => sim.stop()); return svg.node(); }3a:T7dc,playground_canvas = { restart; let polygon = []; mutable polygon_mutable = []; const svg = d3.select(DOM.svg(chart_param.width/2, chart_param.height/2)); const aspect_ratio = chart_param.height/chart_param.width; svg.attr('viewBox',[-10.5, -10.5, 21, 21]); const rect = svg.append("rect") .attr("x",-10.5/aspect_ratio) .attr("y",-10.5) .attr("width", 21/aspect_ratio) .attr("height", 21) .on("mousemove", function(event) { const tracker = svg.selectAll(".tracker") .data([d3.pointer(event)]); tracker.enter().append("circle") .attr("class", "tracker") .attr("r", 0.25) .attr("stroke-width",0.05) .merge(tracker) .attr("transform", d => `translate(${d})`); if (polygon.length) { const polygonVirtual = polygon.slice(); polygonVirtual.push(d3.pointer(event)); polygonVirtual.push(polygon[0]); const polyLine = svg.selectAll(".polyline") .data([polygonVirtual]); polyLine.enter().append("polygon") .attr("class", "polyline") .attr("stroke-width",0.1) .merge(polyLine) .attr("points", d => d.join(" ")); } }) .on("click", function(event){ polygon.push(d3.pointer(event)); mutable polygon_mutable = [...polygon,d3.pointer(event)]; const polyVertices = svg.selectAll(".vertex") .data(polygon); polyVertices.enter().append("circle") .attr("class", "vertex") .attr("r", 0.25) .attr("stroke-width",0.05) .merge(polyVertices) .attr("transform", d => `translate(${d})`); const polyShape = svg.selectAll(".polygon") .data([polygon]); polyShape.enter().append("polygon") .attr("class", "polygon") .merge(polyShape) .attr("stroke-width",0.1) .attr("points", d => d.join(" ")); }) return svg.node(); }3b:T424,playground_animation = { // init plot + get relevant data and keyframe const svg = d3.select(DOM.svg(chart_param.width, chart_param.height)), nodes = large_hexagonal_lattice .map(d => Object.create(d)) .filter(d=>geometric.pointInPolygon([d.x,d.y], polygon_mutable)); const kf = restart_button_polygon; svg.attr('viewBox',[-10.5, -10.5, 21, 21]); const sim = d3.forceSimulation(nodes) .force("lennard-jones-BH",d3.forceLennardJonesPotential().theta(0.125).distanceMax(1.5)) .alphaDecay(0) .velocityDecay(0) .dt(0.05) .tick(); sim.counter = 0; // draw nodes (with uniform radius and no edges) const [node, edge] = draw_nodes(svg, nodes, null, false, false,colorScale); node.call(drag(sim)); // update each node upon simulation tick sim.on("tick", () => { if (sim.counter++ > 500) sim.stop(); node .attr("cx", d => d.x) .attr("cy", d => d.y) .attr("fill", d => colorScale(d.energy)); }); invalidation.then(() => sim.stop()); return svg.node(); }3c:[{"id":"9eb118bf23ac092df9a421a8d9c916f5d61293edfead564ff112dfd7efbbea3911a4beff311352eb82f84204d446d7f8f01cf79b8559607d0562c5de9bc3a265","href":"https://static.observableusercontent.com/files/9eb118bf23ac092df9a421a8d9c916f5d61293edfead564ff112dfd7efbbea3911a4beff311352eb82f84204d446d7f8f01cf79b8559607d0562c5de9bc3a265","name":"spring_compressed@1.svg","type":"image/svg+xml","lastModified":1708981372629.801,"size":7690},{"id":"115407433e203d685394f8455161fcc5b2823344ef100e809f1ff68b126f1821875db46b7d6d69dabb8c5f7935986d8ac39ee699041300278439618f869b7c89","href":"https://static.observableusercontent.com/files/115407433e203d685394f8455161fcc5b2823344ef100e809f1ff68b126f1821875db46b7d6d69dabb8c5f7935986d8ac39ee699041300278439618f869b7c89","name":"spring_equilibrium@1.svg","type":"image/svg+xml","lastModified":1708981395863.615,"size":5944},{"id":"57b4d0b587021c3c4d9b5e438a6d1bdb36a5a5a9fcb3a4ca018a6c2447fda605df1aa210536410ab27db71f276688bd88c396afaa589092ef4d78ff1a940559c","href":"https://static.observableusercontent.com/files/57b4d0b587021c3c4d9b5e438a6d1bdb36a5a5a9fcb3a4ca018a6c2447fda605df1aa210536410ab27db71f276688bd88c396afaa589092ef4d78ff1a940559c","name":"spring_expanded@1.svg","type":"image/svg+xml","lastModified":1708981408246.922,"size":7252}] 2f:{"title":"Teaching Materials Science Concepts using the (d3-)force","cells":[{"id":0,"value":"md`# Teaching Materials Science Concepts using the (d3-)force\n#### Georgios Varnavides* \n---\n*Content initially created in collaboration with Craig Carter and Amina Matt using the Wolfram Language\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":8,"value":"md`\n## Introduction\n>If, in some cataclysm, all of scientific knowledge were to be destroyed, and only one sentence passed on to the next generations of creatures, what statement would contain the most information in the fewest words? I believe it is the atomic hypothesis (or the atomic fact, or whatever you wish to call it) that all things are made of atoms - little particles that move around in perpetual motion, attracting each other when they are a little distance apart, but repelling upon being squeezed into one another. In that one sentence, you will see, there is an enormous amount of information about the world, if just a little imagination and thinking are applied.\n\nRichard Feynman, The Feynman Lectures on Physics, Chapter 1, Volume I (Feynman, Leighton, and Sands)\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1532,"value":"md`\n## Reference\n\n- [Lennard-Jones Potential](#lennardJonesPotentialDoc)\n- [1D Lattice Vibrations](#latticeVibrationsDoc)\n- [Surface Energy](#surfaceEnergyDoc)\n- [Compression Waves And Symmetry Breaking](#compressionWavesDoc)\n- [Dislocations](#dislocationsDoc)\n- [Beyond Lennard-Jones Potentials](#beyondLJDoc)\n- [Playground](#playgroundDoc)\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":30,"value":"lennardJonesPotentialDoc = md`## [#](#lennardJonesPotentialDoc) Lennard-Jones Potential [<>](https://observablehq.com/@gvarnavi/teaching-materials-science-concepts-using-the-d3-force#lennardJonesPotentialDoc)\n\nThe [Lennard-Jones potential](https://en.wikipedia.org/wiki/Lennard-Jones_potential) is a model for a particle's potential energy when it is separated from another particle. \n\nThe model is simple - particles attract each other when they are far apart and repel when they are very near. Nevertheless, the information that we will compute and visualize is enormous - especially when there are many interacting particles.\n\nAll of the physics and its consequences derive from the shape of the potential:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1751,"value":"viewof lj_plot_inputs = Inputs.checkbox(\n [\"show schematics\", \"plot energy\", \"plot force\"],\n {\n value: [\"show schematics\", \"plot energy\"]\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":38,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":79,"value":"md`When the distance is small, the potential **energy** decreases with increasing distance. The **force** is positive and repulsive because force is minus the gradient of potential energy.\n\nWhen the distance is large, the potential **energy** slope is positive and so the **force** is negative: the particles attract. \n\nWe could construct manyfunctions that behave like this; the Lennard-Jones potential is a simple choice. \n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":124,"value":"md`\n### Non-Dimensional Form\nMathematically, the Lennard-Jones potential is given by a linear combination of two terms:\n- An attractive term that goes to zero as ${tex `-\\frac{1}{r^6}`} as ${tex `r \\rightarrow \\infty`}\n- A repulsive term that goes to infinity as ${tex `\\frac{1}{r^{12}}`} as ${tex `r \\rightarrow 0`}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":126,"value":"tex.block`U_{ljp} = \\frac{a}{r^{12}} + \\frac{b}{r^6}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":132,"value":"md`\nwith ${tex`a`} and ${tex`b`} being material-properties we could extract using experimental observations.\n\nWhile this form is useful in investigating material-specific behaviour, we wish to understand the dynamics of materials with this **family** of interatomic potentials. To do this, we can instead derive a general form by non-dimensionalizing our potential by:\n1. introducing a non-dimensional length ${tex `\\rho = \\frac{r}{r_{min}}`}\n2. dividing the potential by an energy minimum ${tex `e_{min}`}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":137,"value":"tex.block`\\tilde{U}_{ljp} = \\frac{1}{\\rho^{12}}-\\frac{2}{\\rho^6}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":179,"value":"md`\n### Harmonic Approximation\nThe plotting range above is slightly misleading. Most of the time, atoms like to be close the energy minimum and fluctuate around that position.\n\nTherefore, before we tackle the full problem - it's instructive to look for an approximate solution near the energy minimum. We can do so by expanding near the bottom of the energy well to second order and obtain the following parabolic (or harmonic) approximation:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":184,"value":"tex.block`\\tilde{U}_{\\textrm{harmonic}} = -1 + 36\\left(\\rho -1\\right)^2 + \\mathcal{O}\\left(\\rho-1\\right)^3`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1781,"value":"viewof lj_harmonic_plot_inputs = Inputs.checkbox(\n [\"show schematics\", \"plot energy\", \"plot harmonic approximation\"],\n {\n value: [\"show schematics\", \"plot energy\", \"plot harmonic approximation\"]\n }\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":141,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":187,"value":"latticeVibrationsDoc = md`---\n\n## [#](#latticeVibrationsDoc) 1D Lattice Vibrations [<>](https://observablehq.com/@gvarnavi/teaching-materials-science-concepts-using-the-d3-force#latticeVibrationsDoc)\n\nLet's start by looking at a (finite) one-dimensional lattice with fixed ends. We will consider two treatments:\n1. Small amplitude vibrations\n - atoms interact only with their nearest neighbors\n - as if they're masses attached with springs (Harmonic approximation)\n2. Larger amplitude vibrations (atoms interact via the Lennard-Jones potential)\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":192,"value":"md`\n### 1. Harmonic Approximation\nWe make a one-dimensional lattice with ${num_of_atoms+1} atoms and start the simulation by displacing the atoms using an initial sine-wave pattern.\n\n__Note:__ We're exaggerating the displacements by a factor of ${Math.round(1/oscillation_magnitude_static)} for visual clarity.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":240,"value":"viewof restart_button_1d = Inputs.button(\"(re)Start Simulation\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":206,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":248,"value":"md`\nWe chose an initial sine-wave displacement in anticipation of the behaviour being periodic in space (and time). In particular, we moved each atom from its equilibrium position by the function:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":259,"value":"sine_wave_plot = Plot.plot({\n width:chart_param.width,\n height:chart_param.width/8,\n marginRight:40,\n x: {grid: true, label:\"Normalized position along chain →\"},\n y: {grid: true,label: \"↑ Normalized displacement\"},\n marks: [\n Plot.line(sine_wave, {x: \"x\", y: \"y\"})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":303,"value":"md`\nWhat if we started with fewer/more 'nodes' along our initial displacement?\nChange the value fo the sub-harmonic mode below to see how this influences the vibrations!\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":214,"value":"viewof sub_harmonic = Inputs.range([1, 6], {step: 1,label: \"Sub-harmonic\",value:3})","pinned":false,"hidden":false,"mode":"ojs"},{"id":309,"value":"md`\nThe visualization using horizontal displacements is hard to see. Let's instead plot the displacements of each atom away from its equilibrium position vertically. \n\n__Note:__ The physics is still in the horizontal direction!`","pinned":false,"hidden":false,"mode":"ojs"},{"id":207,"value":"viewof oscillation_direction = Inputs.radio([\"horizontal\", \"vertical\"], {label: \"Oscillation direction\",value: \"horizontal\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":316,"value":"md`\n### Dispersion Relations\nWe were running the above simulation for a fixed time, but notice how the behaviour is periodic in time. Let's write a helper function to extract that period and stop the simulation.\n\nThe function __compute_single_period()__ calculates one period of oscillation given a particular sub-harmonic mode number:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":319,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":325,"value":"md`\nFor example it takes ${compute_single_period(num_of_atoms,sub_harmonic)} (dimensionless) time-units to do a single periodic oscillation for the ${sub_harmonic}th sub-harmonic mode for ${num_of_atoms} atoms!\n\n\\`\\`\\`javascript\ncompute_single_period(num_of_atoms,sub_harmonic)\n\\`\\`\\`\n${compute_single_period(num_of_atoms,sub_harmonic)}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":334,"value":"md`\nLet's do this for all sub-harmonic modes in the range [${-num_of_atoms/2-5},${num_of_atoms+5}] \nand plot the frequency ${tex`\\left(\\frac{2\\pi}{\\mathrm{period}}\\right)`} against the wave-vector ${tex`\\left( \\frac{2 \\pi \\times \\mathrm{sub-harmonic\\;mode}}{\\mathrm{number\\;of\\;atoms}}\\right)`} \nfor each sub-harmonic mode:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":401,"value":"viewof dispersion_plot = {\n let selection = mouseover_selection;\n return vl.markCircle({size:50})\n .data(dispersion_data)\n .params(selection)\n .encode(\n vl.x().fieldQ('wavevector'),\n vl.y().fieldQ('frequency'),\n vl.color().if(selection, vl.value('steelblue')).value('lightgrey')\n )\n .width(chart_param.width/2)\n .render()\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":512,"value":"md`\nWhat we just plotted is called a dispersion relation. \n\nLet's get a sense of what it means for our lattice vibrations. Try hovering over one of the points in the dispersion plot and notice how the wavelength (inversely related to the wavevector) and frequency (inversely related to the period) change below!\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":462,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":488,"value":"md`\nSub-harmonic mode:${dispersion_subharmonic}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":520,"value":"md`\nThe dispersion behaviour is periodic in the domain ${tex`\\left(-\\frac{\\pi}{a},\\frac{\\pi}{a}\\right)`}, where _a_ is the lattice spacing (${tex`a=1`} in our case).\n\nThis periodic domain is called the __first Brillouin zone__ and it turns out to be a very useful property in materials science.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":526,"value":"md`\n### 2. Lennard Jones Vibrations\nLet's run the same 1D-lattice simulation using the Lennard-Jones potential instead.\nThis means our particles interact according to second, third, etc. nearest neighbor interactions - making the behaviour nonlinear!\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":533,"value":"viewof restart_button_1d_lj = Inputs.button(\"(re)Start Simulation\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":536,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":197,"value":"viewof oscillation_magnitude = Inputs.range([0.0125, 0.05], {step: 0.0125,label: \"Oscillation Magnitude\",value:0.0125})","pinned":false,"hidden":false,"mode":"ojs"},{"id":530,"value":"md`\nThe larger the initial oscillation magnitude - the further away from the harmonic approximation we sample, and thus the more nonlinear behaviour. Try running the simulation above using larger and larger oscillation magnitudes. \n\n__Note:__ The actual oscillation is normalized to be the same in all cases (to allow easier comparison), but notice how the vibrations quickly become nonlinear for larger values!`","pinned":false,"hidden":false,"mode":"ojs"},{"id":659,"value":"legend({color:colorScale1D, title: \"Lennard-Jones energy 1D\",width:chart_param.width/2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":569,"value":"md`\n__Note:__ Finally, notice we've colored each atom according to their Lennard-Jones potential energy.\nThe colors fluctuate as atoms move from high (red) to low (blue) energies, but our fixed edge atoms appear to be mostly in a state of higher energy.\n\nThis because they only have one nearest neighbor (as opposed to two like the 'bulk' atoms), and give us a hint as to the physical origin of surface tension!\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":678,"value":"surfaceEnergyDoc = md`---\n\n## [#](#surfaceEnergyDoc) Surface Energy [<>](https://observablehq.com/@gvarnavi/teaching-materials-science-concepts-using-the-d3-force#surfaceEnergyDoc)\n\nLet's investigate this surface-energy effect by looking at larger interacting systems.\nFor example, some atoms on a square lattice.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":694,"value":"square_surface_energy = {\n // init plot + get relevant data and keyframe\n const svg = d3.select(DOM.svg(chart_param.width, chart_param.height)),\n nodes = square_lattice.map(d => Object.create(d));\n\n svg.attr('viewBox',[-10, -10, 20, 20]);\n\n d3.forceSimulation(nodes)\n .force(\"lennard-jones-BH\",d3.forceLennardJonesPotential().theta(0))\n .alphaDecay(0)\n .velocityDecay(0)\n .dt(0.0125)\n .tick()\n .stop();\n\n // draw nodes (with uniform radius and no edges)\n const [node, edge] = draw_nodes(svg, nodes, null, false, false,colorScale);\n\n mutable static_bulk_energy = d3.sum(nodes.map(d=>d.energy));\n \n return svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":981,"value":"legend({color:colorScale, title: \"Lennard-Jones energy 2D\",width:chart_param.width/2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":706,"value":"md`\nWe've again colored each atom by their total energy, with red being 'higher' energy and dark blue being 'lower' energy. It's evident that the atoms on the surface have a higher energy! \n\nWe can think about this by counting 'satisfied' and 'dangling' bonds. The atoms inside the 'bulk' of the particle all have 4 nearest neighbors (blue), the atoms on the edges of the particle have 3 nearest neighbors (green), while the atoms in the corners of the particle only have two nearest neighbors (orange).\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":709,"value":"md`\n### Surface energy anisotropy\nWe started with a square-particle above, and found that energy is related to its geometry and nearest-neighbour configurations. \n\nWhat can this tell us about the directional anisotropy of the particle's facets? E.g. are 90 degree facets higher or lower energy than 45 degree facets? Let's find out!\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":755,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":720,"value":"viewof theta_in_deg = Inputs.range([-90, 90], {label: \"Slice angle\",value:-30})","pinned":false,"hidden":false,"mode":"ojs"},{"id":827,"value":"md`\nTry moving the slice angle slider above, which splits our square particle down the middle. Can you guess which angles will lead to higher/lower energy facets?\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":946,"value":"md`\nThe surface energy is given by the additional energy the sliced particle has, as compared to the bulk energy of the particle before slicing:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":956,"value":"tex.block`\\frac{E_{\\mathrm{particle}} - N E_{\\mathrm{bulk}}}{2 A}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":964,"value":"md`\nwhere ${tex`N`}, ${tex`A`} are the number of particles and area of the cleaved surface respectively. \nWe can write a function to compute the surface energy for a particular orientation:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":832,"value":"compute_average_energy = (particle, theta_in_deg) => {\n const theta = (theta_in_deg / 180) * Math.PI;\n let length = 16;\n if (Math.abs(theta_in_deg) < 45) {\n length = Math.sqrt(256 + 256 * Math.abs(Math.tan(theta)));\n } else {\n length = Math.sqrt(256 + 256 * Math.abs(1 / Math.tan(theta)));\n }\n\n const [nodes, nodesRight] = slice_particle(square_lattice, theta);\n\n d3.forceSimulation(nodes)\n .force(\"lennard-jones-BH\", d3.forceLennardJonesPotential().theta(0))\n .alphaDecay(0)\n .velocityDecay(0)\n .dt(0)\n .tick()\n .stop();\n\n d3.forceSimulation(nodesRight)\n .force(\"lennard-jones-BH\", d3.forceLennardJonesPotential().theta(0))\n .alphaDecay(0)\n .velocityDecay(0)\n .dt(0)\n .tick()\n .stop();\n\n return (\n (d3.sum([...nodes, ...nodesRight].map((d) => d.energy)) -\n static_bulk_energy) /\n 2\n );\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":971,"value":"md`\nFor example cleaving along ${theta_in_deg}° leads to a (dimensionless) surface energy cost of ${compute_average_energy(square_lattice,theta_in_deg)}\n\n\\`\\`\\`javascript\ncompute_average_energy(square_lattice,theta_in_deg)\n\\`\\`\\`\n${compute_average_energy(square_lattice,theta_in_deg)}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":988,"value":"md`\nLet's do this for all angles in the range [-180°,180°] and plot the surface energy as a polar plot!\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":850,"value":"surface_energy_plot =\nPlot.plot({\n x: {axis:null},\n y: {axis:null},\n marks: [\n Plot.ruleX([0], {stroke: \"black\"}),\n Plot.ruleY([0], {stroke: \"black\"}),\n Plot.dot(surface_energy_data, {x: \"x\", y: \"y\",fill: \"steelblue\",r: 7.5})\n ],\n width:400,height:400,marginLeft:10,marginTop:10,marginRight:10,marginBottom:10\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":997,"value":"md`\nNote that:\n1. Our surface energy plot has (as expected) square symmetry\n2. 90° facets are indeed the lowest energy surfaces\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1000,"value":"md`\n### Lattice Instabilities\nThe square lattice we've constructed above however, cannot be stabilized purely by central forces (like the Lennard-Jones) and is mechanically unstable in 2D.\n\nIt turns out, any slight perturbation (here due to our surface-energy terms, but the same would also happen due to thermal fluctuations in large ('infinite') lattices), will lead to the atoms re-arranging in a 2D closed-packed arrangement!\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1007,"value":"viewof restart_button_square = Inputs.button(\"(re)Start Simulation\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1004,"value":"square_instability_animation = {\n // init plot + get relevant data and keyframe\n const svg = d3.select(DOM.svg(chart_param.width, chart_param.height)),\n nodes = square_lattice.map(d => Object.create(d)),\n kf = restart_button_square;\n \n let counter = 0;\n\n svg.attr('viewBox',[-10, -10, 20, 20]);\n\n const sim = d3.forceSimulation(nodes)\n .force(\"lennard-jones-BH\",\n d3.forceLennardJonesPotential()\n .theta(0.25)\n .distanceMax(1.5)\n )\n .alphaDecay(0)\n .velocityDecay(0)\n .dt(0.0125);\n\n // draw nodes (with uniform radius and no edges)\n const [node, edge] = draw_nodes(svg, nodes, null, false, false,colorScale);\n \n // update each node upon simulation tick\n sim.on(\"tick\", () => {\n \n if (counter++ > 250 || kf == 0) sim.stop();\n \n node\n .attr(\"cx\", d => d.x)\n .attr(\"cy\", d => d.y)\n .attr(\"fill\", d => colorScale(d.energy)); \n\n });\n\n invalidation.then(() => sim.stop());\n return svg.node();\n \n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1560,"value":"legend({color:colorScale, title: \"Lennard-Jones energy 2D\",width:chart_param.width/2})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1015,"value":"md`\n__Note:__ Most atoms (especially in the corners) have re-arranged themselves to have ~6 nearest neighbors!`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1039,"value":"compressionWavesDoc = md`---\n\n## [#](#compressionWavesDoc) Compression Waves and Symmetry Breaking [<>](https://observablehq.com/@gvarnavi/teaching-materials-science-concepts-using-the-d3-force#compressionWavesDoc)\n\nWe'll continue our exploration by looking at stable hexagonal close-packed lattices.\nIn particular, we'll focus on a rotating nanoparticle.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1045,"value":"md`### Compression Waves\nLet's start our nanoparticle spinning at a slow* rate.\nTry running the animation below, and notice the compression waves which form as the particle spins.\n\n*Under normal conditions, most particles would not spin this fast - but we're exploring the physics here, and it helps to exaggerate.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1069,"value":"viewof restart_button_hexagonal = Inputs.button(\"(re)Start Simulation\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1065,"value":"compression_waves_animation = {\n // init plot + get relevant data and keyframe\n const svg = d3.select(DOM.svg(chart_param.width, chart_param.height)),\n nodes = hexagonal_lattice.map((d) => Object.create(d)),\n kf = restart_button_hexagonal;\n\n svg.attr(\"viewBox\", [-12.5, -12.5, 25, 25]);\n\n const sim = d3\n .forceSimulation(nodes)\n .force(\n \"lennard-jones-BH\",\n d3.forceLennardJonesPotential().theta(0.125).distanceMax(1.75)\n )\n .alphaDecay(0)\n .velocityDecay(0)\n .dt(0.05);\n\n sim.counter = 0;\n\n // draw nodes (with uniform radius and no edges)\n const [node, edge] = draw_nodes(svg, nodes, null, false, false, colorScale);\n\n node.call(drag(sim));\n\n // update each node upon simulation tick\n sim.on(\"tick\", () => {\n if (sim.counter++ > 500 || kf == 0) sim.stop();\n\n node\n .attr(\"cx\", (d) => d.x)\n .attr(\"cy\", (d) => d.y)\n .attr(\"fill\", (d) => colorScale(d.energy));\n });\n\n invalidation.then(() => sim.stop());\n return svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1413,"value":"compression_legend = legend({\n color: colorScale,\n title: \"Lennard-Jones energy 2D\",\n width: chart_param.width / 2\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1059,"value":"viewof nanoparticle_initial_velocity_magnitude = Inputs.range([0.225, 0.325], {\n step: 0.0125,\n label: \"Nanoparticle Initial Velocity\",\n value: 0.25\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1091,"value":"md`Let's crank up the initial rotation speed to see how the forces 'rip' the particle apart.\n\nYou might notice the end result no-longer has six-fold symmetry.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1109,"value":"md`__Note:__ The simulation above is interactive! \nTry running the simulation again and dragging one of the atoms (ever so slightly, closed-packed particles experience high forces!)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1243,"value":"dislocationsDoc = md`---\n\n## [#](#dislocationsDoc) Dislocations [<>](https://observablehq.com/@gvarnavi/teaching-materials-science-concepts-using-the-d3-force#dislocationsDoc)\n\nLet's push our simulation to the limit, by _literally_ pushing particles into one another - e.g. in an indentation experiment.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1350,"value":"viewof animation_frame = sim_scrubber(75)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1348,"value":"dislocations_animation = {\n // init plot + get relevant data and keyframe\n const svg = d3.select(DOM.svg(chart_param.width, chart_param.height)),\n nodes = precomputed_animation[animation_frame];\n\n svg.attr('viewBox',[-15, -7.5, 30, 22.5]);\n\n // draw nodes (with uniform radius and no edges)\n const [node, edge] = draw_nodes(svg, nodes, null, false, false,colorScale);\n\n return svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1416,"value":"dislocation_legend = legend({\n color: colorScale,\n title: \"Lennard-Jones energy 2D\",\n width: chart_param.width / 2\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1411,"value":"md`\n__Note:__ Run the slider again slowly, and notice the lines of high energy starting to form near the indenter.\nAs particles get pushed out of their equilibrium positions, the indenter creates extra 'half-planes' of atoms which need to be accommodated in the lattice.\n\nThese are dislocations starting to form, and it's how materials deform plastically!\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1447,"value":"beyondLJDoc = md`---\n\n## [#](#beyondLJDoc) Beyond Lennard-Jones Potentials [<>](https://observablehq.com/@gvarnavi/teaching-materials-science-concepts-using-the-d3-force#beyondLJDoc)\n\nThe Lennard-Jones potential is just a model to guide our exploration of simple physics. There is no (real) justification for its precise form. \n\nIn fact, the choice of the function won't affect the results very much - any potential function with the repulsive/attractive shape above will produce similar consequences.\n\nLet's test that!\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1454,"value":"md`### Generalized Lennard-Jones Potential\nSometimes the Lennard-Jones Potential is called a 6-12 potential because of the exponents used.\nThis generalizes nicely to _nm_-potentials, where _n_ and _m_ refer to the exponents in the interatomic potential:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1458,"value":"tex.block`\\tilde{U}_{g-ljp}(\\rho) = \\frac{n \\rho^{-m}-m \\rho^{-n}}{m-n}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1462,"value":"lj_generalized_plot = Plot.plot({\n x: {\n domain: [0.75, 1.5],\n label: \"Normalized Interatomic Distance, ρ →\"\n },\n y: {\n domain: [-2.75, 2.75],\n grid: true,\n label: null\n },\n marks: [\n Plot.line(generalized_lj_data, {x: \"rho\", y: \"lj_energy\",stroke: \"#4e79a7\"}),\n Plot.line(generalized_lj_data, {x: \"rho\", y: \"generalized_energy\",stroke: \"#e15759\"}),\n Plot.ruleX([1], {stroke: \"black\"}),\n Plot.text(generalized_lj_data.slice(generalized_lj_data.length-1), \n {x: \"rho\",y: \"lj_energy\",text: [\"LJ Energy\"],fill: \"#4e79a7\",dy: \"-1.5em\"}),\n Plot.text(generalized_lj_data.slice(generalized_lj_data.length-50,generalized_lj_data.length-49), \n {x: \"rho\",y: \"generalized_energy\",text: [\"Generalized Energy\"],fill: \"#e15759\",dy: \"-5em\"})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1466,"value":"viewof repulsivePower = Inputs.range([6, 16], {value: 14, step:1, label: \"Repulsive Power\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1471,"value":"viewof attractivePower = Inputs.range([3, 8], {value: 4, step:1, label: \"Attractive Power\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1493,"value":"viewof restart_button_hexagonal_generalized = Inputs.button(\"(re)Start Simulation\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1490,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":1499,"value":"md`### Morse Potential\nWe're also not restricted to LJ-type potentials. The extended \\`d3-force-md\\` module we're using can also handle general central-body potentials given a general function and its (2D) gradient. \nE.g. we can look at the closely related Morse-potential:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1503,"value":"tex.block`\n\\begin{aligned}\n\\tilde{U}_{Morse}(\\rho) &= e^{-2 a (\\rho -1)}-2 e^{a (1-\\rho )} \\\\\n\\frac{\\tilde{U}_{Morse}'(\\rho)}{\\rho} &= \\frac{2 a e^{a (1-\\rho )}-2 a e^{-2 a (\\rho -1)}}{\\rho }\n\\end{aligned}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1515,"value":"morse_plot = Plot.plot({\n x: {\n domain: [0.75, 1.5],\n label: \"Normalized Interatomic Distance, ρ →\"\n },\n y: {\n domain: [-2.75, 2.75],\n grid: true,\n label: null\n },\n marks: [\n Plot.line(morse_data, {x: \"rho\", y: \"lj_energy\",stroke: \"#4e79a7\"}),\n Plot.line(morse_data, {x: \"rho\", y: \"morse_energy\",stroke: \"#e15759\"}),\n Plot.ruleX([1], {stroke: \"black\"}),\n Plot.text(morse_data.slice(morse_data.length-1), \n {x: \"rho\",y: \"lj_energy\",text: [\"LJ Energy\"],fill: \"#4e79a7\",dy: \"-1.5em\"}),\n Plot.text(morse_data.slice(morse_data.length-50,morse_data.length-49), \n {x: \"rho\",y: \"morse_energy\",text: [\"Morse Energy\"],fill: \"#e15759\",dy: \"-5em\"})\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1509,"value":"morsePotential= rho => Math.exp(20-20*rho)-2*Math.exp(10-10*rho)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1511,"value":"morseGradient = rho => 20*Math.exp(10-20*rho)*(Math.exp(10*rho)-Math.exp(10))/rho","pinned":true,"hidden":false,"mode":"ojs"},{"id":1529,"value":"viewof restart_button_hexagonal_morse = Inputs.button(\"(re)Start Simulation\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1525,"value":"$39","pinned":false,"hidden":false,"mode":"ojs"},{"id":1372,"value":"playgroundDoc = md`---\n\n## [#](#playgroundDoc) Playground [<>](https://observablehq.com/@gvarnavi/teaching-materials-science-concepts-using-the-d3-force#playgroundDoc)\n\nYour turn! \n\nUse the canvas to below to draw a particle shape by clicking to add vertices.\nOnce you're satisfied with your shape, scroll down to play around with your particle by dragging (remember, ever so slightly) some of the atoms. \n\nHave fun!\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1119,"value":"viewof restart = Inputs.button(\"Clear Drawing\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1172,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs"},{"id":1375,"value":"viewof restart_button_polygon = Inputs.button(\"(re)Start Simulation\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1186,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs"},{"id":1436,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1433,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":44,"value":"md`\n---\n## Utility Functions`","pinned":false,"hidden":false,"mode":"ojs"},{"id":218,"value":"md`### Imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":41,"value":"d3=require(\"d3@6\",\"d3-force-md@2.2.2\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1556,"value":"md`[\\`d3-force-md\\` Documentation](https://observablehq.com/@gvarnavi/d3-force-md-documentation?collection=@gvarnavi/d3-force-md) `","pinned":false,"hidden":false,"mode":"ojs"},{"id":656,"value":"import {legend} from \"@d3/color-legend\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1114,"value":"geometric = require(\"geometric@2\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1352,"value":"import {Scrubber} from '@mbostock/scrubber'","pinned":false,"hidden":false,"mode":"ojs"},{"id":221,"value":"md`### Plots`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1791,"value":"spring_svgs = [\n {\n x: 0.875,\n y: -2,\n src: await FileAttachment(\"spring_compressed@1.svg\").url()\n },\n {\n x: 1.1,\n y: 2,\n src: await FileAttachment(\"spring_equilibrium@1.svg\").url()\n },\n {\n x: 1.375,\n y: 1,\n src: await FileAttachment(\"spring_expanded@1.svg\").url()\n }\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":48,"value":"lj_data=d3.range(0.75,1.5,0.00125)\n .map(d=> ({\n rho:d,\n energy:Math.pow(d,-12)-2*Math.pow(d,-6),\n force:12*Math.pow(d,-13)-12*Math.pow(d,-7),\n harmonic: 36*Math.pow(d-1,2)-1}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1464,"value":"generalized_lj_data=d3.range(0.75,1.5,0.00125)\n.map(d=> ({\n rho:d,\n lj_energy:Math.pow(d,-12)-2*Math.pow(d,-6),\n generalized_energy:(\n repulsivePower*Math.pow(d,-attractivePower)-attractivePower*Math.pow(d,-repulsivePower))/\n (attractivePower-repulsivePower)}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1517,"value":"morse_data=d3.range(0.75,1.5,0.00125)\n.map(d=> ({\n rho:d,\n lj_energy:Math.pow(d,-12)-2*Math.pow(d,-6),\n morse_energy:morsePotential(d)}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":257,"value":"sine_wave = [...Array(201)]\n.map((d,i) => ({\n x: i/200,\n y: Math.sin(sub_harmonic*(2*Math.PI*i)/200)\n}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":335,"value":"dispersion_data=d3.range(-num_of_atoms/2-5,num_of_atoms+5)\n .map(d => ({wavevector:2*Math.PI*d/num_of_atoms,\n frequency:2*Math.PI/compute_single_period(num_of_atoms,d)}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":595,"value":"mouseover_selection = vl.selectPoint().on('mouseover').fields('wavevector').nearest('true')","pinned":false,"hidden":false,"mode":"ojs"},{"id":593,"value":"dispersion_mouseover = Generators.observe(next => {\n next([]);\n const signal = mouseover_selection.name();\n const change = (name, value) => next(value?.vlPoint?.or ?? []);\n dispersion_plot.addSignalListener(signal, change);\n return () => dispersion_plot.removeSignalListener(signal, change)[0];\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":224,"value":"md`### 1D Lattice`","pinned":false,"hidden":false,"mode":"ojs"},{"id":194,"value":"num_of_atoms = 50","pinned":false,"hidden":false,"mode":"ojs"},{"id":552,"value":"oscillation_magnitude_static = 0.0125","pinned":false,"hidden":false,"mode":"ojs"},{"id":211,"value":"one_dimensional_lattice = [...Array(num_of_atoms+1)]\n.map((d,i) => ({\n x: i + oscillation_magnitude_static*Math.sin(sub_harmonic*(2*Math.PI*i)/num_of_atoms),\n y: 0, fy: 0\n}))\n// fix end nodes\n.map(d => (d.x == 0 || d.x == num_of_atoms) ? {x:d.x,fx:d.x,y:d.y,fy:d.y} : d)","pinned":false,"hidden":false,"mode":"ojs"},{"id":230,"value":"one_dimensional_lattice_links = [...Array(num_of_atoms)]\n.map((d,i) => ({source:i,target:i+1}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":321,"value":"function make_one_dimensional_lattice(num_of_atoms,sub_harmonic,oscillation_magnitude) {\n \n const nodes = [...Array(num_of_atoms+1)]\n .map((d,i) => ({\n x: i + oscillation_magnitude*Math.sin(sub_harmonic*(2*Math.PI*i)/num_of_atoms),\n y: 0, fy: 0\n }))\n .map(d => (d.x == 0 || d.x == num_of_atoms) ? {x:d.x,fx:d.x,y:d.y,fy:d.y} : d);\n\n const edges = [...Array(num_of_atoms)]\n .map((d,i) => ({source:i,target:i+1}));\n \n return [nodes,edges]\n \n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":443,"value":"dispersion_subharmonic = {\n if (dispersion_mouseover.length == 0) {\n return 0 \n }\n else { \n if ('wavevector' in dispersion_mouseover[0]) {\n return Math.round(dispersion_mouseover[0].wavevector*num_of_atoms/(2*Math.PI))\n } \n else {\n return 0\n }\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":540,"value":"colorScale1D = d3.scaleSequential(d3.interpolateTurbo)\n .domain([-3.075,-3.025])","pinned":false,"hidden":false,"mode":"ojs"},{"id":681,"value":"md`### Square Lattice`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1073,"value":"colorScale = d3.scaleSequential(d3.interpolateTurbo)\n .domain([-7,-2])","pinned":false,"hidden":false,"mode":"ojs"},{"id":684,"value":"square_lattice = [...Array((8*2 + 1)**2)]\n.map((d,i) => ({\n x: 0.975*(i%17 - 8),\n y: 0.975*(Math.floor(i/17)-8)\n}))","pinned":true,"hidden":false,"mode":"ojs"},{"id":726,"value":"slice_particle = (particle, theta) => {\n const particle_copy = particle.map(d => Object.create(d));\n let left,right;\nleft = particle_copy.filter(d => Math.tan(theta%(Math.PI))*(d.x+1e-6) < d.y);\n right = particle_copy.filter(x => !left.includes(x));\n return [left, right]\n}\n ","pinned":false,"hidden":false,"mode":"ojs"},{"id":854,"value":"mutable static_bulk_energy = null","pinned":false,"hidden":false,"mode":"ojs"},{"id":922,"value":"surface_energy_data_quarter=d3.range(0,90,7.5)\n .map(d => ({theta:d,energy:compute_average_energy(square_lattice,d)}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":928,"value":"surface_energy_data = d3.range(0,360,7.5)\n .map( (d,i) => ({\n x:surface_energy_data_quarter[i%surface_energy_data_quarter.length].energy*Math.cos(d*Math.PI/180),\n y:surface_energy_data_quarter[i%surface_energy_data_quarter.length].energy*Math.sin(d*Math.PI/180)\n }))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1048,"value":"md`### Hexagonal Lattice`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1055,"value":"calc_squared_distance = (x1, x2) => x1.reduce(\n (sum, x1_i, i) => sum + Math.pow(x1_i - x2[i],2),\n 0\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1051,"value":"hexagonal_lattice_vectors = [[1,0],[1/2,Math.sqrt(3)/2]]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1054,"value":"hexagonal_lattice = [...Array((12*2 + 1)**2)]\n// Create square lattice sites\n.map((d,i) => ({\n x: (i%25 - 12),\n y: (Math.floor(i/25)-12)\n}))\n// Transform to hexagonal lattice sites\n.map(d => ({\n x: d.x*hexagonal_lattice_vectors[0][0]+d.y*hexagonal_lattice_vectors[1][0],\n y: d.x*hexagonal_lattice_vectors[0][1]+d.y*hexagonal_lattice_vectors[1][1]\n}))\n// Filter according to norm\n.filter(d => calc_squared_distance([d.x,d.y],[0,0])<10**2)\n// Start it spinning\n.map(d=>({...d,\n vx:nanoparticle_initial_velocity_magnitude*Math.sqrt(calc_squared_distance(\n [d.x,d.y],[0,0]))*Math.cos(Math.atan2(d.y, d.x)+Math.PI/2),\n vy:nanoparticle_initial_velocity_magnitude*Math.sqrt(calc_squared_distance(\n [d.x,d.y],[0,0]))*Math.sin(Math.atan2(d.y, d.x)+Math.PI/2)\n }))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1615,"value":"inside_range = (dx, r = 0.0175) => (Math.abs(dx) < r ? dx : Math.sign(dx) * r)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1095,"value":"drag = (simulation, r = 0.0175) => {\n function dragstarted(event) {\n if (!event.active) simulation.restart(), (simulation.counter = 100);\n event.subject.fx = event.subject.x;\n event.subject.fy = event.subject.y;\n }\n\n function dragged(event) {\n event.subject.fx += inside_range(event.x - event.subject.fx, r);\n event.subject.fy += inside_range(event.y - event.subject.fy, r);\n }\n\n function dragended(event) {\n event.subject.fx = null;\n event.subject.fy = null;\n }\n\n return d3\n .drag()\n .on(\"start\", dragstarted)\n .on(\"drag\", dragged)\n .on(\"end\", dragended);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1245,"value":"md`### Beam`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1318,"value":"large_hex_particles = [...Array((30*2 + 1)**2)]\n// Create square lattice sites\n.map((d,i) => ({\n x: (i%61 - 30),\n y: (Math.floor(i/61)-30)\n}))\n// Transform to hexagonal lattice sites\n.map(d => ({\n x: d.x*hexagonal_lattice_vectors[0][0]+d.y*hexagonal_lattice_vectors[1][0],\n y: d.x*hexagonal_lattice_vectors[0][1]+d.y*hexagonal_lattice_vectors[1][1]\n}))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1248,"value":"beam_particles = large_hex_particles\n// Carve out beam\n.filter(d => d.y>-0.5 && d.y<12.5 && d.x<17.5 && d.x>-17.5)\n// Fix beam edges\n.map(d => ((d.x <-15) || (d.x >15)) ? {x:d.x,y:d.y,fx:d.x,fy:d.y} : d)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1314,"value":"indenter_particles = large_hex_particles\n// Carve out indenter\n.filter(d => \n (15+6*d.x >= Math.sqrt(3)*(25+2*d.y)) \n && (25+2*d.y >= 0) \n && (6*d.x + Math.sqrt(3)*(25+2*d.y)<= 15)\n )\n// Start indenter moving\n.map(d => (d.y <-5) ? {x:d.x,y:(d.y+7.5),vx:0,vy:0.75} : d)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1342,"value":"precomputed_animation = {\n let steps = [];\n const nodes = beam_particles.map((d) => Object.create(d)),\n indenter_nodes = indenter_particles.map((d) => Object.create(d));\n\n const sim = d3\n .forceSimulation([...nodes, ...indenter_nodes])\n .force(\n \"lennard-jones-BH\",\n d3.forceLennardJonesPotential().theta(0.125).distanceMax(1.75)\n )\n .force(\"zero\", () => {\n indenter_nodes.forEach((n) => {\n n.force_x = 0;\n n.force_y = 0;\n });\n })\n .alphaDecay(0)\n .velocityDecay(0)\n .dt(0.0375);\n\n for (let i = 0; i < 75; i++) {\n sim.tick();\n const copies = new Map(\n [...nodes, ...indenter_nodes].map((n) => [\n n,\n { x: n.x, y: n.y, energy: n.energy }\n ])\n );\n steps.push(Array.from(copies.values()));\n sim.stop();\n }\n\n return steps;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1212,"value":"md`### Make Your Own`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1193,"value":"large_hexagonal_lattice = [...Array((20*2 + 1)**2)]\n// Create square lattice sites\n.map((d,i) => ({\n x: (i%41 - 20),\n y: (Math.floor(i/41)-20)\n}))\n// Transform to hexagonal lattice sites\n.map(d => ({\n x: d.x*hexagonal_lattice_vectors[0][0]+d.y*hexagonal_lattice_vectors[1][0],\n y: d.x*hexagonal_lattice_vectors[0][1]+d.y*hexagonal_lattice_vectors[1][1]\n}))\n.filter(d => (d.x >-10 && d.x < 10 && d.y>-10 && d.y < 10))","pinned":false,"hidden":false,"mode":"ojs"},{"id":1220,"value":"mutable polygon_mutable = []","pinned":false,"hidden":false,"mode":"ojs"},{"id":227,"value":"md`### General Utilities`","pinned":false,"hidden":false,"mode":"ojs"},{"id":233,"value":"draw_nodes = (svg, node_data, edge_data, w_r, w_l, colorScale) => {\n \n if (colorScale == null) {\n colorScale = e => null\n }\n \n let edge = null;\n if (edge_data) {\n edge = svg.selectAll(\".edge\")\n .data(edge_data).enter()\n .append(\"line\")\n .classed(\"edge\", true)\n .attr(\"x1\", d => node_data[d.source].x)\n .attr(\"y1\", d => node_data[d.source].y)\n .attr(\"x2\", d => node_data[d.target].x)\n .attr(\"y2\", d => node_data[d.target].y)\n .style(\"stroke\", w_l ? \"#bbb\" : \"none\");\n }\n \n const node = svg.selectAll(\".node\")\n .data(node_data).enter()\n .append(\"circle\")\n .classed(\"node\", true)\n .attr(\"cx\", d => d.x)\n .attr(\"cy\", d => d.y)\n .attr(\"fill\",d=>colorScale(d.energy))\n .attr(\"r\", d => w_r ? d.r : 0.5);\n \n return [node, edge];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":237,"value":"chart_param = ({\n width: width,\n height: 500,\n margin: {\n top: 10,\n right: 10,\n bottom: 10,\n left: 10,\n center: 150\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1355,"value":"sim_scrubber = (n_frames) =>\n Scrubber(\n Array.from({ length: n_frames }, (_, i) => i),\n {\n format: (i) => \"frame of simulation: \" + (i + 1),\n autoplay: false,\n loop: true,\n alternate: true\n }\n )","pinned":false,"hidden":false,"mode":"ojs"},{"id":1127,"value":"html`\n`","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W3c","stdlib":"1","titling":"auto"}