import molpy as mp

# SMILES → 3D
mol = mp.io.smiles.SmilesIr("CCO").to_atomistic()
mol, _ = mp.conformer.Conformer(seed=42).generate(mol)

# OPLS-AA typing
typifier = mp.ff.typifier.OplsAaTypifier()
frame = typifier.typify(mol).to_frame()
ff = typifier.forcefield()

# minimize
potentials = mp.ff.compile.PotentialCompiler(ff).compile(frame)
frame, report = mp.optimize.Lbfgs(potentials, fmax=0.05).minimize(frame)

# a unit with two ports
eo = mp.conformer.Conformer(seed=42).generate(
    mp.io.smiles.SmilesIr.from_fragment("[<]OCC[>]").to_template()
)[0]
assembler = mp.builder.Assembler({"EO": eo}, mp.builder.GrowthPlacer())

def build(cgsmiles):
    sites = mp.io.cgsmiles.CgSmilesIr(cgsmiles).to_coarsegrain()
    return assembler.assemble(sites, mp.Atomistic)

# shape from CGsmiles
chain = build("{[#EO]|10}")
ring = build("{[#EO]1[#EO][#EO][#EO][#EO][#EO]1}")

from molpy.builder import CarbonTubeBuilder, GrapheneBuilder

# (10,10) nanotube
tube = CarbonTubeBuilder(10, 10, length=50.0, periodic=True)
cnt = tube.build()
cnt.box = tube.cell(vacuum=12.0)

# graphene sheet
sheet = GrapheneBuilder(8, 8, periodic_xy=True).build()

frame.box = mp.Box.cube(30.0)
# molecule IDs
frame["atoms"]["mol_id"] = (
    mp.Topology.from_frame(frame).connected_components() + 1
)
# cutoffs
ff.get_style("pair", "lj/cut")["cutoff"] = 10.0
ff.get_style("pair", "coul/cut")["cutoff"] = 10.0

lammps = mp.engine.LammpsEngine(check_executable=False)
lammps.generate_inputs(frame, ff, "./lammps")

gromacs = mp.engine.GromacsEngine(check_executable=False)
gromacs.generate_inputs(frame, ff, "./gromacs")

from molpy.compute import Rdf

# 8 Å neighbor list
nl = mp.core.NeighborList(8.0)
nl.build(frame.coords, frame.box)

# g(r)
rdf = Rdf(n_bins=50, r_max=8.0, r_min=0.0).compute(pairs=[nl.neighbors()])