import numpy as np
import espressomd
from espressomd import electrostatics, shapes, electrostatic_extensions


# Parameters 

box_l_xy = 10
box_l_z  = 10
gap_l = 3. 

box_l = np.array((box_l_xy, box_l_xy, box_l_z))
gap = np.array((0., 0., gap_l))

time_step = 1e-10    # Arbitraty small value
lam_bjerrum = 1e-10  # Small, such that particle interactions are negligible

# Set Up System

sys = espressomd.System(box_l = box_l + gap)
sys.periodicity = [1, 1, 1]
sys.time_step = time_step

# Add one positive particle

sys.part.add(
  id = 0,
  type = 0,
  pos = (0, 0, 1),
  q = 1
  )

# Add one negative particle

sys.part.add(
  id = 1,
  type = 1,
  pos = (0, 0, 9),
  q = -1
  )

# Electrostatic Interactions

solver = electrostatics.P3M(
        prefactor = lam_bjerrum,
        accuracy = 1e-3,
        )
sys.actors.add(solver)

# Add metallic electrodes with potential of +-0.5V

solver_elc = electrostatic_extensions.ELC(
        gap_size = gap_l,
        maxPWerror = 1e-3,
        delta_mid_top = -1,
        delta_mid_bot = -1,
        const_pot = 1,
        pot_diff = 1,
        )
sys.actors.add(solver_elc)

# Calculate coulomb energy

E = sys.analysis.energy()

# Here, I get "0.4000000000303713" but would expect 0.8
print(E["coulomb"])