demo2 morse chain

%23%20%2F%2F%2F%20script%0A%23%20requires-python%20%3D%20%22%3E%3D3.11%22%0A%23%20dependencies%20%3D%20%5B%0A%23%20%20%20%20%20%22marimo%22%2C%0A%23%20%20%20%20%20%22numpy%22%2C%0A%23%20%20%20%20%20%22matplotlib%22%2C%0A%23%20%20%20%20%20%22marimo-precompute%3E%3D0.2.3%22%2C%0A%23%20%5D%0A%23%20%2F%2F%2F%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.21.1%22%0Aapp%20%3D%20marimo.App()%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Demo%202%3A%20Morse%20chain%20%E2%80%94%20which%20bond%20snaps%3F%0A%0A%20%20%20%20A%20chain%20of%206%20atoms%20connected%20by%20Morse%20bonds%20is%20pulled%20apart%20by%20a%0A%20%20%20%20force%20applied%20to%20the%20right%20end.%20Each%20nearest-neighbour%20bond%20has%20a%0A%20%20%20%20different%20well%20depth%20%24D%24%2C%20giving%20each%20a%20different%20maximum%20sustainable%0A%20%20%20%20tension%20%24F_%7B%5Cmax%7D%20%3D%20D%5Calpha%2F2%24.%20The%20**weakest%20bond**%20(smallest%20%24D%24)%0A%20%20%20%20reaches%20its%20fold%20first%20and%20snaps.%0A%0A%20%20%20%20In%20this%20notebook%20we%20use%20LACT%20to%20trace%20the%20force%E2%80%93extension%20S-curve%0A%20%20%20%20through%20the%20fold%2C%20revealing%20the%20full%20unstable%20branch%20that%0A%20%20%20%20quasi%20static%20loading%20cannot%20reach.%20The%20bond-length%20plot%20confirms%20that%20only%0A%20%20%20%20the%20weakest%20bond%20snaps%20%E2%80%94%20the%20others%20remain%20near%20equilibrium.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20We%20first%20load%20the%20necessary%20dependencies%20and%20define%20some%20helper%20functions%2C%20which%20include%20creating%20a%20LAMMPS%20instance%2C%20passing%20on%20the%20simulation%20setup%20and%20the%20update%20command%20associated%20with%20the%20continuation%20parameter%20(magnitude%20of%20the%20force%20applied%20to%20the%20atom%20on%20the%20right)%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20from%20marimo_precompute%20import%20persistent_cache%2C%20prefetch_all%0A%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20from%20LACT.precomputed%20import%20PrecomputedSystem%0A%20%20%20%20except%20ImportError%3A%0A%20%20%20%20%20%20%20%20class%20PrecomputedSystem%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20def%20__init__(self%2C%20d)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20self.natoms%20%3D%20int(d%5B%22natoms%22%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20self.U_0%20%3D%20np.asarray(d%5B%22U_0%22%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20self.data%20%3D%20%7B%22Y_s%22%3A%20%5Bnp.asarray(y)%20for%20y%20in%20d%5B%22Y_s%22%5D%5D%2C%20%22ds_s%22%3A%20list(d.get(%22ds_s%22%2C%20%5B%5D))%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20%22energies%22%20in%20d%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20self.data%5B%22energies%22%5D%20%3D%20list(d%5B%22energies%22%5D)%0A%20%20%20%20return%20PrecomputedSystem%2C%20mo%2C%20np%2C%20persistent_cache%2C%20plt%2C%20prefetch_all%0A%0A%0A%40app.cell%0Aasync%20def%20_(prefetch_all)%3A%0A%20%20%20%20await%20prefetch_all()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20n_atoms%20%3D%206%0A%20%20%20%20alpha%20%3D%206.0%20%20%20%20%20%20%20%23%20Morse%20width%20parameter%0A%20%20%20%20r0%20%3D%201.12%20%20%20%20%20%20%20%20%20%23%20Morse%20equilibrium%20distance%0A%20%20%20%20%23%20Well%20depths%20per%20bond%20%E2%80%94%20ordered%20weakest%20to%20strongest%0A%20%20%20%20bond_info%20%3D%20%5B%0A%20%20%20%20%20%20%20%20((1%2C%202)%2C%200.85)%2C%20%20%23%20weakest%20%E2%80%94%20snaps%20first%0A%20%20%20%20%20%20%20%20((2%2C%203)%2C%200.90)%2C%0A%20%20%20%20%20%20%20%20((3%2C%204)%2C%200.95)%2C%0A%20%20%20%20%20%20%20%20((4%2C%205)%2C%201.00)%2C%0A%20%20%20%20%20%20%20%20((5%2C%206)%2C%201.05)%2C%20%20%23%20strongest%0A%20%20%20%20%5D%0A%0A%20%20%20%20make_chain%20%3D%20None%0A%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20from%20lammps%20import%20lammps%0A%20%20%20%20%20%20%20%20from%20LACT%20import%20atom_cont_system%0A%0A%20%20%20%20%20%20%20%20def%20make_chain()%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_lmp%20%3D%20lammps(cmdargs%3D%5B%22-screen%22%2C%20%22none%22%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20_lmp.commands_string(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22%22%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20units%20%20%20%20%20%20%20%20%20lj%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20dimension%20%20%20%20%203%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20boundary%20%20%20%20%20%20f%20f%20f%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20atom_style%20%20%20%20atomic%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20atom_modify%20%20%20map%20yes%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20region%20box%20block%20-50.0%2050.0%20-50.0%2050.0%20-50.0%2050.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20create_box%20%7Bn_atoms%7D%20box%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20mass%20*%201.0%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20create_atoms%201%20single%20%7B0%20*%20r0%3A.6f%7D%200.0%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20create_atoms%202%20single%20%7B1%20*%20r0%3A.6f%7D%200.0%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20create_atoms%203%20single%20%7B2%20*%20r0%3A.6f%7D%200.0%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20create_atoms%204%20single%20%7B3%20*%20r0%3A.6f%7D%200.0%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20create_atoms%205%20single%20%7B4%20*%20r0%3A.6f%7D%200.0%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20create_atoms%206%20single%20%7B5%20*%20r0%3A.6f%7D%200.0%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20reset_atom_ids%20sort%20yes%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pair_style%20morse%205.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pair_coeff%20*%20*%200.001%20%7Balpha%7D%20%7Br0%7D%200.001%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pair_coeff%201%202%20%7Bbond_info%5B0%5D%5B1%5D%7D%20%7Balpha%7D%20%7Br0%7D%205.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pair_coeff%202%203%20%7Bbond_info%5B1%5D%5B1%5D%7D%20%7Balpha%7D%20%7Br0%7D%205.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pair_coeff%203%204%20%7Bbond_info%5B2%5D%5B1%5D%7D%20%7Balpha%7D%20%7Br0%7D%205.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pair_coeff%204%205%20%7Bbond_info%5B3%5D%5B1%5D%7D%20%7Balpha%7D%20%7Br0%7D%205.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pair_coeff%205%206%20%7Bbond_info%5B4%5D%5B1%5D%7D%20%7Balpha%7D%20%7Br0%7D%205.0%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20group%20left_end%20id%201%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20group%20right_end%20id%206%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fix%20freeze%20left_end%20setforce%200.0%200.0%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fix%20pull%20right_end%20addforce%200.0%200.0%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fix_modify%20pull%20energy%20yes%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20compute%20forces%20all%20property%2Fatom%20fx%20fy%20fz%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20compute%20ids%20all%20property%2Fatom%20id%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20def%20update_command(force)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20return%20f%22%22%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20unfix%20pull%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fix%20pull%20right_end%20addforce%20%7Bforce%7D%200.0%200.0%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20fix_modify%20pull%20energy%20yes%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20atom_cont_system(_lmp%2C%20update_command)%0A%20%20%20%20except%20ImportError%3A%0A%20%20%20%20%20%20%20%20pass%0A%20%20%20%20return%20alpha%2C%20bond_info%2C%20make_chain%2C%20n_atoms%2C%20r0%0A%0A%0A%40app.cell%0Adef%20_(alpha%2C%20n_atoms%2C%20np%2C%20r0)%3A%0A%20%20%20%20def%20get_chain_data(sys_obj)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Extract%20extensions%2C%20forces%2C%20and%20bond%20lengths%20from%20a%20chain%20system.%22%22%22%0A%20%20%20%20%20%20%20%20_U0%20%3D%20sys_obj.U_0%0A%20%20%20%20%20%20%20%20_extensions%2C%20_forces%2C%20_bonds%20%3D%20%5B%5D%2C%20%5B%5D%2C%20%5B%5D%0A%20%20%20%20%20%20%20%20for%20_Y%20in%20sys_obj.data%5B%22Y_s%22%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_pos%20%3D%20_U0%20%2B%20_Y%5B%3A-1%5D.reshape(n_atoms%2C%203)%0A%20%20%20%20%20%20%20%20%20%20%20%20_extensions.append(_pos%5B-1%2C%200%5D%20-%20_U0%5B-1%2C%200%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20_forces.append(_Y%5B-1%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20_bonds.append(%5B_pos%5Bi%20%2B%201%2C%200%5D%20-%20_pos%5Bi%2C%200%5D%20for%20i%20in%20range(n_atoms%20-%201)%5D)%0A%20%20%20%20%20%20%20%20return%20np.array(_extensions)%2C%20np.array(_forces)%2C%20np.array(_bonds)%0A%0A%20%20%20%20%23%20Analytical%20fold%20distance%20and%20fold%20forces%20per%20bond%0A%20%20%20%20r_fold%20%3D%20r0%20%2B%20np.log(2)%20%2F%20alpha%0A%20%20%20%20return%20get_chain_data%2C%20r_fold%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Quasi-static%20loading%20(for%20comparison)%0A%0A%20%20%20%20First%20we%20try%20the%20naive%20approach%3A%20ramp%20the%20applied%20force%20in%20small%20steps%0A%20%20%20%20from%200%20to%204.0%20using%20LAMMPS%20energy%20minimisation%20at%20each%20step.%0A%20%20%20%20This%20works%20up%20to%20the%20fold%20(%24F_%7B%5Cmax%7D%20%3D%200.85%20%5Ctimes%206%2F2%20%3D%202.55%24)%2C%0A%20%20%20%20where%20the%20bond%20snaps%20and%20the%20atom%20flies%20to%20the%20pair-potential%20cutoff.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(PrecomputedSystem%2C%20get_chain_data%2C%20make_chain%2C%20persistent_cache)%3A%0A%20%20%20%20def%20_run_qs()%3A%0A%20%20%20%20%20%20%20%20_s%20%3D%20make_chain()%0A%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_s.quasi_static_run(0.0%2C%200.02%2C%20200%2C%20verbose%3DFalse)%0A%20%20%20%20%20%20%20%20except%20Exception%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20pass%20%20%23%20LAMMPS%20crashes%20past%20the%20fold%20%E2%80%94%20keep%20the%20points%20we%20got%0A%20%20%20%20%20%20%20%20return%20%7B%22natoms%22%3A%20_s.natoms%2C%20%22U_0%22%3A%20_s.U_0%2C%20%22Y_s%22%3A%20_s.data%5B%22Y_s%22%5D%7D%0A%20%20%20%20with%20persistent_cache(name%3D%22demo2_qs%22)%3A%0A%20%20%20%20%20%20%20%20qs_data%20%3D%20_run_qs()%0A%20%20%20%20qs_sys%20%3D%20PrecomputedSystem(qs_data)%0A%20%20%20%20qs_extensions%2C%20qs_forces%2C%20qs_bonds%20%3D%20get_chain_data(qs_sys)%0A%20%20%20%20return%20qs_bonds%2C%20qs_extensions%2C%20qs_forces%2C%20qs_sys%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Arclength%20continuation%0A%0A%20%20%20%20Now%20we%20use%20the%20continuation%20routine%20in%20LACT.%20A%20short%20quasi-static%0A%20%20%20%20ramp%20seeds%20the%20path%20with%20a%20few%20points%2C%20then%20%60continuation_run%60%0A%20%20%20%20traces%20smoothly%20through%20the%20fold%20and%20onto%20the%20unstable%20branch.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(PrecomputedSystem%2C%20get_chain_data%2C%20make_chain%2C%20persistent_cache)%3A%0A%20%20%20%20def%20_run_cont()%3A%0A%20%20%20%20%20%20%20%20_s%20%3D%20make_chain()%0A%20%20%20%20%20%20%20%20_s.quasi_static_run(0.0%2C%200.5%2C%205%2C%20verbose%3DFalse)%0A%20%20%20%20%20%20%20%20_n_qs%20%3D%20len(_s.data%5B%22Y_s%22%5D)%0A%20%20%20%20%20%20%20%20_s.continuation_run(%0A%20%20%20%20%20%20%20%20%20%20%20%20n_iter%3D200%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20ds_default%3D0.01%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20ds_smallest%3D0.001%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20ds_largest%3D0.05%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20verbose%3DFalse%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20checkpoint_freq%3D0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20cont_target%3D0.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20target_tol%3D0.1%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20return%20%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%22natoms%22%3A%20_s.natoms%2C%20%22U_0%22%3A%20_s.U_0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22Y_s%22%3A%20_s.data%5B%22Y_s%22%5D%2C%20%22n_qs%22%3A%20_n_qs%2C%0A%20%20%20%20%20%20%20%20%7D%0A%20%20%20%20with%20persistent_cache(name%3D%22demo2_cont%22)%3A%0A%20%20%20%20%20%20%20%20cont_data%20%3D%20_run_cont()%0A%20%20%20%20cont_sys%20%3D%20PrecomputedSystem(cont_data)%0A%20%20%20%20n_qs%20%3D%20cont_data%5B%22n_qs%22%5D%0A%20%20%20%20cont_extensions%2C%20cont_forces%2C%20cont_bonds%20%3D%20get_chain_data(cont_sys)%0A%20%20%20%20return%20cont_bonds%2C%20cont_extensions%2C%20cont_forces%2C%20cont_sys%2C%20n_qs%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20We%20can%20compare%20the%20two%20approaches%20visually%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(cont_extensions%2C%20cont_forces%2C%20n_qs%2C%20plt%2C%20qs_extensions%2C%20qs_forces)%3A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(7%2C%205))%0A%0A%20%20%20%20_ax.plot(qs_forces%2C%20qs_extensions%2C%20%22s%22%2C%20ms%3D5%2C%20color%3D%22C2%22%2C%20alpha%3D0.7%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22Quasi-static%20(load-stepping)%22)%0A%20%20%20%20_ax.plot(cont_forces%5Bn_qs%3A%5D%2C%20cont_extensions%5Bn_qs%3A%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%22-%22%2C%20lw%3D1.5%2C%20color%3D%22C0%22%2C%20label%3D%22Arclength%20continuation%22)%0A%0A%20%20%20%20_ax.set_xlabel(%22Applied%20force%20%24F%24%22)%0A%20%20%20%20_ax.set_ylabel(%22Extension%20(displacement%20of%20atom%206)%22)%0A%20%20%20%20_ax.set_title(%22Force%E2%80%93extension%20through%20the%20first%20fold%22)%0A%20%20%20%20_ax.legend()%0A%20%20%20%20_ax.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(%0A%20%20%20%20alpha%2C%0A%20%20%20%20bond_info%2C%0A%20%20%20%20cont_bonds%2C%0A%20%20%20%20cont_forces%2C%0A%20%20%20%20n_qs%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20qs_bonds%2C%0A%20%20%20%20qs_forces%2C%0A%20%20%20%20r0%2C%0A%20%20%20%20r_fold%2C%0A)%3A%0A%20%20%20%20_fig%2C%20(_ax_top%2C%20_ax_bot)%20%3D%20plt.subplots(2%2C%201%2C%20figsize%3D(8%2C%208)%2C%20sharex%3DTrue)%0A%0A%20%20%20%20for%20_idx%2C%20((_i%2C%20_j)%2C%20_D%2C%20)%20in%20enumerate(bond_info)%3A%0A%20%20%20%20%20%20%20%20_c%20%3D%20f%22C%7B_idx%7D%22%0A%20%20%20%20%20%20%20%20_Fmax%20%3D%20_D%20*%20alpha%20%2F%202%0A%20%20%20%20%20%20%20%20_ax_top.plot(qs_forces%2C%20qs_bonds%5B%3A%2C%20_idx%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22s%22%2C%20ms%3D4%2C%20color%3D_c%2C%20alpha%3D0.5)%0A%20%20%20%20%20%20%20%20_ax_top.plot(cont_forces%5Bn_qs%3A%5D%2C%20cont_bonds%5Bn_qs%3A%2C%20_idx%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22-%22%2C%20lw%3D1.5%2C%20color%3D_c%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22Bond%20%7B_i%7D%E2%80%93%7B_j%7D%20(%24D%24%3D%7B_D%3A.2f%7D%2C%20%24F_%7B%7Bmax%7D%7D%24%3D%7B_Fmax%3A.2f%7D)%22)%0A%0A%20%20%20%20_ax_top.axhline(r_fold%2C%20color%3D%22k%22%2C%20ls%3D%22--%22%2C%20lw%3D0.8%2C%20alpha%3D0.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22Fold%20at%20%24r%5E*%24%3D%7Br_fold%3A.3f%7D%22)%0A%20%20%20%20_ax_top.set_ylabel(r%22Bond%20length%20%24r%24%22)%0A%20%20%20%20_ax_top.set_title(%22Bond%20lengths%20vs.%20applied%20force%22)%0A%20%20%20%20_ax_top.legend(fontsize%3D7%2C%20loc%3D%22upper%20left%22)%0A%20%20%20%20_ax_top.grid(True%2C%20alpha%3D0.3)%0A%0A%20%20%20%20for%20_idx%2C%20((_i%2C%20_j)%2C%20_D)%20in%20enumerate(bond_info)%3A%0A%20%20%20%20%20%20%20%20_c%20%3D%20f%22C%7B_idx%7D%22%0A%20%20%20%20%20%20%20%20_ax_bot.plot(qs_forces%2C%20qs_bonds%5B%3A%2C%20_idx%5D%20-%20r0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22s%22%2C%20ms%3D4%2C%20color%3D_c%2C%20alpha%3D0.5)%0A%20%20%20%20%20%20%20%20_ax_bot.plot(cont_forces%5Bn_qs%3A%5D%2C%20cont_bonds%5Bn_qs%3A%2C%20_idx%5D%20-%20r0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22-%22%2C%20lw%3D1.5%2C%20color%3D_c%2C%20label%3Df%22Bond%20%7B_i%7D%E2%80%93%7B_j%7D%22)%0A%20%20%20%20_ax_bot.set_xlabel(%22Applied%20force%20%24F%24%22)%0A%20%20%20%20_ax_bot.set_ylabel(r%22Bond%20stretch%20%24r%20-%20r_0%24%22)%0A%20%20%20%20_ax_bot.set_title(%22Bond%20stretch%3A%20weakest%20bond%20absorbs%20most%20deformation%22)%0A%20%20%20%20_ax_bot.legend(fontsize%3D7)%0A%20%20%20%20_ax_bot.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Interpreting%20the%20results%0A%0A%20%20%20%20**Force%E2%80%93extension%20curve%3A**%20The%20green%20squares%20(quasi-static%0A%20%20%20%20loading)%20follow%20the%20stable%20branch%20up%20to%20the%20fold%20at%0A%20%20%20%20%24F_%7B%5Cmax%7D%20%3D%202.55%24%2C%20where%20the%20LAMMPS%20minimiser%20crashes%20%E2%80%94%20an%20atom%0A%20%20%20%20is%20lost.%20The%20blue%20line%20(arclength%20continuation)%20traces%20smoothly%0A%20%20%20%20through%20the%20fold%20onto%20the%20unstable%20branch%20(decreasing%20force).%0A%0A%20%20%20%20**Bond%20lengths%3A**%20Under%20increasing%20applied%20force%2C%20all%20bonds%20stretch%2C%0A%20%20%20%20but%20the%20weakest%20bond%20(%24D%3D0.85%24)%20stretches%20the%20most.%20At%20the%20fold%2C%20it%0A%20%20%20%20passes%20the%20inflection%20point%20%24r%5E*%24%20and%20begins%20to%20snap%20%E2%80%94%20its%20length%0A%20%20%20%20increases%20rapidly%20while%20the%20other%20bonds%20contract%20(since%20the%20total%0A%20%20%20%20force%20is%20now%20decreasing).%20This%20confirms%20that%20**disorder%20determines%0A%20%20%20%20which%20bond%20breaks**%3A%20the%20weakest%20link%20in%20the%20chain.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Representative%20configurations%0A%0A%20%20%20%20The%20chain%20at%20key%20points%20along%20the%20continuation%20curve%3A%20near%20equilibrium%2C%0A%20%20%20%20at%20the%20fold%2C%20and%20on%20the%20unstable%20branch.%20The%20arrow%20shows%20the%20applied%0A%20%20%20%20force%20on%20the%20rightmost%20atom.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(cont_extensions%2C%20cont_forces%2C%20cont_sys%2C%20n_atoms%2C%20n_qs%2C%20np%2C%20plt)%3A%0A%20%20%20%20_n_cont%20%3D%20len(cont_sys.data%5B%22Y_s%22%5D)%0A%20%20%20%20_i_fold%20%3D%20n_qs%20%2B%20np.argmax(cont_forces%5Bn_qs%3A%5D)%0A%20%20%20%20_snapshots%20%3D%20%5B%0A%20%20%20%20%20%20%20%20(0%2C%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Equilibrium%20(F%3D0)%22)%2C%0A%20%20%20%20%20%20%20%20(_i_fold%2C%20%20%20%20%20%20%20%20%22At%20the%20fold%22)%2C%0A%20%20%20%20%20%20%20%20(_n_cont%20-%201%2C%20%20%20%20%22Unstable%20branch%22)%2C%0A%20%20%20%20%5D%0A%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(1%2C%203%2C%20figsize%3D(15%2C%203.5))%0A%20%20%20%20_U0%20%3D%20cont_sys.U_0%0A%20%20%20%20_max_force%20%3D%20np.max(np.abs(cont_forces))%0A%0A%20%20%20%20for%20_ax%2C%20(_idx%2C%20_title)%20in%20zip(_axes%2C%20_snapshots)%3A%0A%20%20%20%20%20%20%20%20_Y%20%3D%20cont_sys.data%5B%22Y_s%22%5D%5B_idx%5D%0A%20%20%20%20%20%20%20%20_dev%20%3D%20_Y%5B%3A-1%5D.reshape(n_atoms%2C%203)%0A%20%20%20%20%20%20%20%20_pos%20%3D%20_U0%20%2B%20_dev%0A%20%20%20%20%20%20%20%20_F%20%3D%20_Y%5B-1%5D%0A%20%20%20%20%20%20%20%20_ext%20%3D%20cont_extensions%5B_idx%5D%0A%0A%20%20%20%20%20%20%20%20%23%20Atoms%3A%20fixed%20%3D%20grey%2C%20others%20%3D%20C0%0A%20%20%20%20%20%20%20%20_colors%20%3D%20%5B%22grey%22%5D%20%2B%20%5B%22C0%22%5D%20*%20(n_atoms%20-%201)%0A%20%20%20%20%20%20%20%20_ax.scatter(_pos%5B%3A%2C%200%5D%2C%20%5B0%5D%20*%20n_atoms%2C%20s%3D250%2C%20c%3D_colors%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20edgecolors%3D%22k%22%2C%20lw%3D0.8%2C%20zorder%3D5)%0A%20%20%20%20%20%20%20%20_ax.annotate(%22fixed%22%2C%20(_pos%5B0%2C%200%5D%2C%20-0.2)%2C%20ha%3D%22center%22%2C%20fontsize%3D7%2C%20color%3D%22grey%22)%0A%20%20%20%20%20%20%20%20_ax.annotate(%22pulled%22%2C%20(_pos%5B-1%2C%200%5D%2C%20-0.2)%2C%20ha%3D%22center%22%2C%20fontsize%3D7%2C%20color%3D%22C0%22)%0A%0A%20%20%20%20%20%20%20%20%23%20Bond%20lines%0A%20%20%20%20%20%20%20%20for%20_i%20in%20range(n_atoms%20-%201)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.plot(%5B_pos%5B_i%2C%200%5D%2C%20_pos%5B_i%20%2B%201%2C%200%5D%5D%2C%20%5B0%2C%200%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22k--%22%2C%20lw%3D0.6%2C%20alpha%3D0.4)%0A%0A%20%20%20%20%20%20%20%20%23%20Force%20arrow%0A%20%20%20%20%20%20%20%20if%20abs(_F)%20%3E%200.01%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_arrow_len%20%3D%200.4%20*%20_F%20%2F%20_max_force%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax.annotate(%22%22%2C%20xy%3D(_pos%5B-1%2C%200%5D%20%2B%20_arrow_len%2C%200)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20xytext%3D(_pos%5B-1%2C%200%5D%2C%200)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20arrowprops%3Ddict(arrowstyle%3D%22-%3E%22%2C%20color%3D%22red%22%2C%20lw%3D2))%0A%0A%20%20%20%20%20%20%20%20_ax.set_title(f%22%7B_title%7D%5Cnext%20%3D%20%7B_ext%3A.3f%7D%2C%20%24F%20%3D%20%7B_F%3A.2f%7D%24%22)%0A%20%20%20%20%20%20%20%20_ax.set_xlim(_pos%5B0%2C%200%5D%20-%200.4%2C%20_pos%5B-1%2C%200%5D%20%2B%200.8)%0A%20%20%20%20%20%20%20%20_ax.set_ylim(-0.4%2C%200.4)%0A%20%20%20%20%20%20%20%20_ax.set_aspect(%22equal%22)%0A%20%20%20%20%20%20%20%20_ax.set_yticks(%5B%5D)%0A%20%20%20%20%20%20%20%20_ax.grid(True%2C%20alpha%3D0.2)%0A%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Interactive%20explorer%0A%0A%20%20%20%20Pick%20a%20system%20and%20slide%20along%20the%20solution%20curve.%20The%20left%20panel%0A%20%20%20%20shows%20the%20chain%20configuration%3B%20the%20right%20panel%20shows%20the%0A%20%20%20%20force%E2%80%93extension%20curve%20with%20a%20red%20dot%20marking%20the%20current%20point.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20sys_dropdown%20%3D%20mo.ui.dropdown(%0A%20%20%20%20%20%20%20%20options%3D%7B%22Quasi-static%22%3A%20%22qs%22%2C%20%22Continuation%22%3A%20%22cont%22%7D%2C%0A%20%20%20%20%20%20%20%20value%3D%22Continuation%22%2C%0A%20%20%20%20%20%20%20%20label%3D%22System%22%2C%0A%20%20%20%20)%0A%20%20%20%20sys_dropdown%0A%20%20%20%20return%20(sys_dropdown%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(cont_sys%2C%20mo%2C%20qs_sys%2C%20sys_dropdown)%3A%0A%20%20%20%20_systems%20%3D%20%7B%22qs%22%3A%20qs_sys%2C%20%22cont%22%3A%20cont_sys%7D%0A%20%20%20%20_n_pts%20%3D%20len(_systems%5Bsys_dropdown.value%5D.data%5B%22Y_s%22%5D)%0A%20%20%20%20idx_slider%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%20start%3D0%2C%20stop%3D_n_pts%20-%201%2C%20step%3D1%2C%20value%3D0%2C%0A%20%20%20%20%20%20%20%20label%3D%22Index%20along%20curve%22%2C%0A%20%20%20%20)%0A%20%20%20%20idx_slider%0A%20%20%20%20return%20(idx_slider%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(%0A%20%20%20%20cont_extensions%2C%0A%20%20%20%20cont_forces%2C%0A%20%20%20%20cont_sys%2C%0A%20%20%20%20get_chain_data%2C%0A%20%20%20%20idx_slider%2C%0A%20%20%20%20n_atoms%2C%0A%20%20%20%20n_qs%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20qs_extensions%2C%0A%20%20%20%20qs_forces%2C%0A%20%20%20%20qs_sys%2C%0A%20%20%20%20sys_dropdown%2C%0A)%3A%0A%20%20%20%20_systems%20%3D%20%7B%22qs%22%3A%20qs_sys%2C%20%22cont%22%3A%20cont_sys%7D%0A%20%20%20%20_active%20%3D%20_systems%5Bsys_dropdown.value%5D%0A%20%20%20%20_idx%20%3D%20idx_slider.value%0A%20%20%20%20_Y%20%3D%20_active.data%5B%22Y_s%22%5D%5B_idx%5D%0A%20%20%20%20_U0%20%3D%20_active.U_0%0A%20%20%20%20_dev%20%3D%20_Y%5B%3A-1%5D.reshape(n_atoms%2C%203)%0A%20%20%20%20_pos%20%3D%20_U0%20%2B%20_dev%0A%20%20%20%20_F%20%3D%20_Y%5B-1%5D%0A%20%20%20%20_act_ext%2C%20_%2C%20_%20%3D%20get_chain_data(_active)%0A%20%20%20%20_ext%20%3D%20_act_ext%5B_idx%5D%0A%20%20%20%20_max_force%20%3D%20max(np.max(np.abs(cont_forces))%2C%20np.max(np.abs(qs_forces)))%0A%0A%20%20%20%20_fig%2C%20(_ax_l%2C%20_ax_r)%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(14%2C%205)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20gridspec_kw%3D%7B%22width_ratios%22%3A%20%5B1%2C%201.5%5D%7D)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Left%3A%20chain%20configuration%20%E2%94%80%E2%94%80%0A%20%20%20%20_colors%20%3D%20%5B%22grey%22%5D%20%2B%20%5B%22C0%22%5D%20*%20(n_atoms%20-%201)%0A%20%20%20%20_ax_l.scatter(_pos%5B%3A%2C%200%5D%2C%20%5B0%5D%20*%20n_atoms%2C%20s%3D350%2C%20c%3D_colors%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20edgecolors%3D%22k%22%2C%20lw%3D0.8%2C%20zorder%3D5)%0A%20%20%20%20_ax_l.annotate(%22fixed%22%2C%20(_pos%5B0%2C%200%5D%2C%20-0.25)%2C%20ha%3D%22center%22%2C%20fontsize%3D8%2C%20color%3D%22grey%22)%0A%20%20%20%20_ax_l.annotate(%22pulled%22%2C%20(_pos%5B-1%2C%200%5D%2C%20-0.25)%2C%20ha%3D%22center%22%2C%20fontsize%3D8%2C%20color%3D%22C0%22)%0A%0A%20%20%20%20for%20_i%20in%20range(n_atoms%20-%201)%3A%0A%20%20%20%20%20%20%20%20_ax_l.plot(%5B_pos%5B_i%2C%200%5D%2C%20_pos%5B_i%20%2B%201%2C%200%5D%5D%2C%20%5B0%2C%200%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22k--%22%2C%20lw%3D0.6%2C%20alpha%3D0.4)%0A%0A%20%20%20%20if%20abs(_F)%20%3E%200.01%3A%0A%20%20%20%20%20%20%20%20_arrow_len%20%3D%200.4%20*%20_F%20%2F%20_max_force%0A%20%20%20%20%20%20%20%20_ax_l.annotate(%22%22%2C%20xy%3D(_pos%5B-1%2C%200%5D%20%2B%20_arrow_len%2C%200)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20xytext%3D(_pos%5B-1%2C%200%5D%2C%200)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20arrowprops%3Ddict(arrowstyle%3D%22-%3E%22%2C%20color%3D%22red%22%2C%20lw%3D2.5))%0A%0A%20%20%20%20_ax_l.set_title(f%22ext%20%3D%20%7B_ext%3A.3f%7D%2C%20%24F%20%3D%20%7B_F%3A.2f%7D%24%22)%0A%20%20%20%20_ax_l.set_xlim(_pos%5B0%2C%200%5D%20-%200.4%2C%20_pos%5B-1%2C%200%5D%20%2B%200.8)%0A%20%20%20%20_ax_l.set_ylim(-0.5%2C%200.5)%0A%20%20%20%20_ax_l.set_aspect(%22equal%22)%0A%20%20%20%20_ax_l.set_yticks(%5B%5D)%0A%20%20%20%20_ax_l.grid(True%2C%20alpha%3D0.2)%0A%0A%20%20%20%20%23%20%E2%94%80%E2%94%80%20Right%3A%20force-extension%20with%20current%20point%20%E2%94%80%E2%94%80%0A%20%20%20%20_ax_r.plot(qs_forces%2C%20qs_extensions%2C%20%22s%22%2C%20ms%3D4%2C%20color%3D%22C2%22%2C%20alpha%3D0.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22Quasi-static%22)%0A%20%20%20%20_ax_r.plot(cont_forces%5Bn_qs%3A%5D%2C%20cont_extensions%5Bn_qs%3A%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22-%22%2C%20lw%3D1.2%2C%20color%3D%22C0%22%2C%20alpha%3D0.5%2C%20label%3D%22Continuation%22)%0A%0A%20%20%20%20_ax_r.plot(_F%2C%20_ext%2C%20%22o%22%2C%20ms%3D12%2C%20color%3D%22red%22%2C%20zorder%3D10%2C%20label%3D%22Current%20point%22)%0A%0A%20%20%20%20_ax_r.set_xlabel(%22Applied%20force%20%24F%24%22)%0A%20%20%20%20_ax_r.set_ylabel(%22Extension%20(displacement%20of%20atom%206)%22)%0A%20%20%20%20_ax_r.set_title(%22Force%E2%80%93extension%20curve%22)%0A%20%20%20%20_ax_r.legend(fontsize%3D8)%0A%20%20%20%20_ax_r.grid(True%2C%20alpha%3D0.3)%0A%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A