import espressomd
import espressomd.interactions

lj_eps = 1.0
lj_sig = 1.0
lj_cut = 1.12246

system = espressomd.System(box_l=[4. * lj_cut, 10.0, 10.0])
system.cell_system.skin = 0.1
system.time_step = .1
max_interaction_cutoff = min(system.box_l / system.cell_system.node_grid) - system.cell_system.skin * (system.cell_system.get_state()['n_nodes'] > 1)
print(f'maximal interaction cutoff = {max_interaction_cutoff}')

system.non_bonded_inter[0, 0].lennard_jones.set_params(
    epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")

harm = espressomd.interactions.HarmonicBond(k=1.0, r_0=0.1, r_cut=lj_cut)
angle = espressomd.interactions.AngleCosine(bend=1., phi0=3.14)
dihe = espressomd.interactions.Dihedral(bend=1., mult=1, phase=0.)
system.bonded_inter.add(harm)
system.bonded_inter.add(angle)
system.bonded_inter.add(dihe)

p1, p2, p3, p4 = system.part.add(
    pos=[[0., 0., 0.], [0.1, 0., 0.], [0.1, 0.1, 0.], [0.1, 0.1, 0.05]])
p1.add_bond((harm, p2))
p2.add_bond((angle, p1, p3))
p2.add_bond((dihe, p1, p3, p4))