#####################################################
#Trying to reproduce seg error
#####################################################

import sys
import espressomd
import espressomd.magnetostatics as magnetostatics
import numpy as np
import random

#Check relevant features
print(espressomd.features())
required_features = ['EXCLUSIONS', 'FFTW', 'GSL', 'P3M','GAY_BERNE']

espressomd.assert_features(required_features)

#####################################################
#Gay Berne Parameters and Density
#####################################################
n_part = 50

#Zanoni Gay-Berne parameters
sig_s = 1
sige_to_sig_s = 0.345
sige = 0.345
#print(sig_s)

#Mu to nu
mu = 1
nu = 3

#Epse and epss
epse_to_epss = 5
epss = 1
epse = 5
r_cut_gb = 1.4

###############################
#System initialization
##############################

system = espressomd.System(box_l=20*np.ones(3))
system.periodicity = [True, True, True]
system.time_step = 0.001
system.cell_system.skin = 0.4

#Place particles
for p in range(int(n_part)):
        pos1 = [random.random()*20,random.random()*20,random.random()*20]
        system.part.add(type=0, pos=pos1, rotation=[True,True,True], dip=[0,0,1])

print("Set particles")

############################################
#Warm-up & simulation
############################################

#Reset integrator
system.integrator.set_vv()
system.thermostat.set_langevin(kT=1, gamma=1, seed=42)

#Define interactions
system.non_bonded_inter[0,0].gay_berne.set_params(eps=epss,sig=sige,cut=r_cut_gb,k1=sige_to_sig_s,k2=epse_to_epss,mu=mu,nu=nu)

p3m = magnetostatics.DipolarP3M(prefactor=1, mesh=32, accuracy=1E-3)
system.integrator.run(0, recalc_forces=True)
system.actors.add(p3m)

#Warm-up
system.force_cap = 100

for w in range(1000):
    system.integrator.run(1000)
    print(system.analysis.energy())

system.force_cap = 0

print(system.analysis.energy())

print("Survived!")

exit()