#!/usr/bin/env python
# coding: utf-8

import hoomd
import hoomd.md
import numpy as np
import sys

# random seeds
seed1 = np.random.randint(0, 2147483647)
seed2 = np.random.randint(0, 2147483647)
seed3 = np.random.randint(0, 2147483647)

# some parameters
d_BH = 1.10688
box = 512
act_vel = 125
rot_dif = 3/(d_BH**2)
dens = 0.3
dist = np.sqrt((np.pi*(0.5*d_BH)**2)/dens)

# hoomd init
hoomd.context.initialize("");
system = hoomd.init.create_lattice(unitcell=hoomd.lattice.sq(a=dist), n=[box,box]);
all = hoomd.group.all();
N = len(all);
nl = hoomd.md.nlist.cell();
lj = hoomd.md.pair.lj(r_cut=2**(1/6), nlist=nl);
lj.set_params(mode='shift');
lj.pair_coeff.set('A', 'A', epsilon=100.0, sigma=1.0);

# set activity with random direction and fixed magnitude act_vel
activity = []
for i in range(N):
    ra = np.random.rand(1)[0] * 2.0 * np.pi
    activity.append((act_vel * np.cos(ra), act_vel * np.sin(ra), 0))
    
hoomd.md.force.active(group=all,
                      seed=seed1,
                      f_lst=activity,
                      rotation_diff=rot_dif,
                      orientation_link=False,
                      orientation_reverse_link=True);

hoomd.md.integrate.mode_standard(dt=1E-6);
hoomd.md.integrate.brownian(group=all, kT=1.0, seed=seed2);

# dumping data
fname = "traj.gsd"
hoomd.dump.gsd(fname, period=4E5, group=all, overwrite=True, dynamic=['attribute', 'momentum', 'topology']);

hoomd.run(1E9, quiet=True);