import espressomd
import espressomd.magnetostatics
import espressomd.visualization
import numpy as np

np.random.seed(seed=42)
system = espressomd.System(box_l=[16., 16., 16.])
system.time_step = 0.001
system.cell_system.skin = 0.2

system.part.add(pos=np.random.random((100, 3)) * system.box_l,
                dip=np.random.random((100, 3)) - 0.5)

system.non_bonded_inter[0, 0].wca.set_params(epsilon=1., sigma=1.)
system.integrator.set_steepest_descent(f_max=10., gamma=1., max_displacement=0.01)
system.integrator.run(100)
system.integrator.set_vv()
system.thermostat.set_langevin(kT=1., seed=42, gamma=10.)
system.magnetostatics.solver = espressomd.magnetostatics.DipolarP3M(
    prefactor=1., accuracy=1e-2, mesh=16)

visualizer = espressomd.visualization.openGLLive(system, director_arrows=True)
visualizer.run(1)