<jupyter_start><jupyter_text>LJ Tutorial  The target trajectory files and the target rdfs in this example were generated with the code below (hidden, click to reveal):<jupyter_code># This cell takes 40 minutes. Don't uncomment and rerun it unless you have a good reason.
# trajectory{kT}.gsd files were created here

#import freud
#import gsd
#import gsd.hoomd
#import hoomd
#import hoomd.md
#import numpy as np
#import matplotlib.pyplot as plt


#kTs = [0.5, 1.5, 2.0]
#n = 12
#n_particles = n**3

#for i,kT in enumerate(kTs):
#    hoomd.context.initialize("")
#    system = hoomd.init.create_lattice(
#        unitcell = (hoomd.lattice.sc(a=1.58, type_name="A")), 
#        n = n,
#    )
#    
#    nl = hoomd.md.nlist.cell()
#    lj = hoomd.md.pair.lj(r_cut=2.5, nlist=nl)
#    lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0)
#    hoomd.md.integrate.mode_standard(dt=0.001)
#    _all = hoomd.group.all()
#    nvt = hoomd.md.integrate.nvt(group=_all, kT=kT, tau=1)
#    nvt.randomize_velocities(seed=23)
#    hoomd.analyze.log(a
#        filename=f'LJ_kT{kT}.log',
#        quantities=["time", "temperature", "potential_energy"],
#        period=100,
#        overwrite=True
#    )
#    hoomd.dump.gsd(f'trajectory{kT}.gsd', period=5e3, group=_all, overwrite=True)
#    hoomd.run(1e6)

# Check that the system is equilibrated.

#fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True)
#fig.suptitle("Potential Energy vs Timestep")
#for i, kT in enumerate(kTs):   
#    log = np.loadtxt(f'state{i}/LJ_kT{kT}.log', skiprows=1)
#    ax1.plot(log[:,0],log[:,3], label=f"kT={kT}")
#    ax2.plot(log[-50:,0],log[-50:,3], label=f"kT={kT}")
#    ax1.set_xlabel("Timestep")
#    ax2.set_xlabel("Timestep")
#    ax1.set_ylabel("Potential Energy (AU)")
#    ax1.legend()
#    ax2.legend()

#for i,kT in enumerate(kTs):
#    gsdfile = f'trajectory{kT}.gsd'
#    with gsd.hoomd.open(gsdfile) as t:
#        rdf = freud.density.RDF(bins=101, r_max=5)
#        for snap in t[-50:]:
#            rdf.compute(system=snap, reset=False)
#        data = np.stack((rdf.bin_centers, rdf.rdf)).T
#        np.savetxt(f"rdfs/rdftarget{kT}", data)
#        plt.plot(rdf.bin_centers, rdf.rdf)
#        plt.title(f"State {i} kT = {kT}")
#        plt.xlabel("r")
#        plt.ylabel("g(r)")
#        plt.show()
import itertools

import matplotlib.pyplot as plt
import numpy as np

from msibi import MSIBI, State, Pair, mie<jupyter_output><empty_output><jupyter_text>Remove any previous runs.<jupyter_code>%%bash
rm -rf state* potentials<jupyter_output><empty_output><jupyter_text>Initialize the states.<jupyter_code>kTs = [0.5, 1.5, 2.0]<jupyter_output><empty_output><jupyter_text>Create MSIBI instancePass in values that will be used.<jupyter_code>opt = MSIBI(
    rdf_cutoff=5.0,
    n_rdf_points=101,
    max_frames=10,
    pot_cutoff=5.0,
    smooth_rdfs=True,
    verbose=True
)<jupyter_output><empty_output><jupyter_text>Specify statesEach state should have it's own corresponding trajectory (gsd) file. Create a state by passing in a name, trajectory file, and alpha value. The `MSIBI` instance created earlier needs to be passed in as well<jupyter_code>opt.add_state(
    State(name="A", kT=0.5, traj_file="trajectory0.5.gsd", alpha=0.25)
)
opt.add_state(
    State(name="B", kT=1.5, traj_file="trajectory1.5.gsd", alpha=0.50)
)
opt.add_state(
    State(name="C", kT=2.0, traj_file="trajectory2.0.gsd", alpha=0.25)
)<jupyter_output><empty_output><jupyter_text>Specify pairsEach type of possible pair-wise interaction of the system should have its own `Pair` created. Each `Pair` object needs to have the states created above added.<jupyter_code>initial_guess = mie(opt.pot_r, 1.0, 1.0)
opt.add_pair(Pair(type1="A", type2="A", potential=initial_guess))<jupyter_output>/Users/jenny/Projects/msibi/msibi/potentials.py:13: RuntimeWarning: divide by zero encountered in true_divide
  return prefactor * eps * ((sig / r) ** m - (sig / r) ** n)
/Users/jenny/Projects/msibi/msibi/potentials.py:13: RuntimeWarning: invalid value encountered in subtract
  return prefactor * eps * ((sig / r) ** m - (sig / r) ** n)<jupyter_text>Implement MSIBI algorithmRun the optimization steps. Pass in lists of any state and pair objects created above.<jupyter_code>opt.optimize(n_iterations=5, engine='hoomd')<jupyter_output>Using HOOMD version 2.
-------- Iteration 0 --------
Launching 4 CPU threads...
Running state A on CPU
Running state C on CPU
Running state B on CPU<jupyter_text>Plot results<jupyter_code>for state in opt.states:
    plt.title(f'{state.name} kt={state.kT}')
    for step in range(0, opt.n_iterations):
        step_rdf = np.loadtxt(f'{state.dir}/pair_A-A-state_{state.name}-step{step}.txt')
        plt.plot(step_rdf[:,0],step_rdf[:,1], label=f'step {step}')
      
    target = np.loadtxt(f"rdfs/rdftarget{state.kT}")
    plt.plot(target[:,0], target[:,1], label='target')
    plt.legend()
    plt.show()
for step in range(0, opt.n_iterations):
    plt.ylim([-1.5,0.5])
    plt.xlim([0.5,2])
    plt.title("MSIBI potentials by iteration")

    step_pot = np.loadtxt(f'potentials/step{step}.pot.A-A.txt')
    plt.plot(step_pot[:,0],step_pot[:,1], label=f'step {step}')

plt.legend()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1585014409221902, 'fraction_numerical': 0.024015369836695485, 'mean_word_length': 4.139190523198421, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '460565', 'license': 'MIT', 'n_tokens_mistral': 2149, 'n_tokens_neox': 2087, 'n_words': 420, 'path': 'msibi/tutorials/lj/LJ_Tutorial.ipynb', 'repo_name': 'jennyfothergill/msibi', 'size': '5205'}