#! /usr/bin/env hoomd

from hoomd_script import *
import math

# parameters
phi_P = 0.20
n_poly = 2371
T = 1.2
polymer1 = dict(bond_len=1.2, type=['A']*10 + ['B']*7 + ['A']*10, bond="linear", count=n_poly)

# perform some simple math to find the length of the box
N = len(polymer1['type']) * polymer1['count'];
L = math.pow(math.pi * N / (6.0 * phi_P), 1.0/3.0);

# generate the polymer system
init.create_random_polymers(box=hoomd.BoxDim(L), polymers=[polymer1], separation=dict(A=0.35, B=0.35));

# force field setup
harmonic = bond.harmonic()
harmonic.set_coeff('polymer', k=330.0, r0=0.84)
lj = pair.lj(r_cut=3.0)
lj.pair_coeff.set('A', 'A', epsilon=1.0, sigma=1.0, alpha=0.0)
lj.pair_coeff.set('A', 'B', epsilon=1.0, sigma=1.0, alpha=0.0)
lj.pair_coeff.set('B', 'B', epsilon=1.0, sigma=1.0, alpha=1.0)

# NVT integration
all = group.all()
integrate.mode_standard(dt=0.005)
integrate.nvt(group=all, T=T, tau=0.5)

# warm up the system
run(4000)

nlist.set_params(r_buff=0.4, check_period=5)
# these paramters ^^^         ^^^            are determined automatically by running with the next line uncommented
#tune.r_buff(set_max_check_period=True)

# quick benchmark
run(2000, profile=True)

# run the full benchmark
run(50000, limit_hours = 5.0 / 60.0)