capillary_potential.py

#!/usr/bin/python
#
# capillary_potential.py
# Written by: Lars Konermann, June 2025
#
# -------------------------------------------------------------------------
# MODELING THE ELECTROSTATICS OF AN ESI CAPILLARY AND ITS COUNTER ELECTRODE
# -------------------------------------------------------------------------
# PLEASE CITE:
# 
# "A Simple Method for Determining Charge Distributions, Potentials, and Electric Fields
#  Associated with an Electrospray Capillary"
#  Lars Konermann
#  J. Am. Soc. Mass Spectrom. 2025
#
#
# See below for
# - what this program does
# - input file info
# - output file info
#
# run this program by typing "python3 ./capillary_potential.py   <input_file.txt>"
#
# 
# WHAT THIS PROGRAM DOES:
# =======================
#
# This program builds a single layer ESI capillary consisting of concentric rings (cap_rings) of beads.
# Because of symmetry, all beads in a given cap_ring have the same charge.
# These charges are optimized, such that the surface potential is equal everywhere,
# as expected for an electrically conducting capillary (e.g. iron) connected to
# a high voltage power supply.
# The surface potential at each cap_probe_out site j is phi = SUM(qi/rij)
# 
# The program can run without or with counter electrode (cel). The cel (if present) is modeled as a grounded
# planar metal disk
#
# Every capillary bead has the same size, and the spherical capillary diameter at every z-value assumes that
# beads are in direct contact, in a hexagonally packed pattern.
# Arbitrary capillary shapes can be generated by specifying the number of beads in each cap_ring.
#
# The program can operate in two modes:
#
# Mode 1: NO COUNTER ELECTRODE
# cel = "no"  means that we are dealing with an isolated capillary,
# where the potential becomes zero at infinite distance.
#
# capillary layout and cap_ring numbering (*s represent cap_probe_out sites):
# (cap_probe_ins sites are not shown, they are just on the INSide of the capillary)
#                                          
#   CAPILLARY (pos charge)                 
#   *   *   *   *                                    .
#   9+  8+  7+  6+  *                                .
#   9+  8+  7+  6+  5+  *   *   *   *   *            .
#   9+  8+  7+  6+  5+  4+  3+  2+  1+  0+           .
#   9+  8+  7+  6+  5+  4+  3+  2+  1+  0+           .
#   9+  8+  7+  6+  5+  4+  3+  2+  1+  0+           .
#   9+  8+  7+  6+  5+                               . 
#   9+  8+  7+                                       .
#                                                    .
# <---------negative---------------------------------0--------------------positive---->  x-axis  
#                                       |                        
#                                 -(x_cel_dist)            
#
#  (when used in this mode, x_cel_dist specifies the distance from ring0 to yz plane at x=0)
#                                  
# Mode 2: WITH GROUNDED COUNTER ELECTRODE
# cel = "yes"  means that the capillary is facing a grounded flat counter electrode
# The potential at the counter electrode is zero.
# The distance from capillary tip (cap_ring 0) to the counter electrode is x_cel_dist.
# We first deal with this situation using the image charge method, where every bead has an oppositely
# charged mirror image.
#
# capillary layout and cap_ring numbering (dots represent cap_probe_out sites):                  
#
#
#                                                     
#   ACTUAL CAPILLARY (pos charge)                    | grounded                 IMAGE CHARGE CAPILLARY
#   *   *   *   *                                    | counter
#   9+  8+  7+  6+  *                                | electrode has pot = 0             6-  7-  8-  9-
#   9+  8+  7+  6+  5+  *   *   *   *   *            |                               5-  6-  7-  8-  9-
#   9+  8+  7+  6+  5+  4+  3+  2+  1+  0+           |           0-  1-  2-  3-  4-  5-  6-  7-  8-  9-
#   9+  8+  7+  6+  5+  4+  3+  2+  1+  0+           |           0-  1-  2-  3-  4-  5-  6-  7-  8-  9-
#   9+  8+  7+  6+  5+  4+  3+  2+  1+  0+           |           0-  1-  2-  3-  4-  5-  6-  7-  8-  9-
#   9+  8+  7+  6+  5+                               |                               5-  6-  7-  8-  9-
#   9+  8+  7+                                       |                                   6-  7-  8-  9-
#
# <---------negative---------------------------------0--------------------positive---->  x-axis  
#                                       |                        |
#                                 -(x_cel_dist)            +(x_cel_dist)
#
# Once charges on the actual capillary are optimized in the presence of the image charge capillary,
# actual capillary charges are retained and an explicit counter electrode is added at x = 0.
# The (negative) charges on the explicit counter electrode are adjusted to ensure that the counter
# electrode has a zero potential.
#
#
#                                                     
#   ACTUAL CAPILLARY (pos charge)                   -| grounded                 
#   *   *   *   *                                   -| counter
#   9+  8+  7+  6+  *                               -| electrode has pot = 0 and carries negative charges           
#   9+  8+  7+  6+  5+  *   *   *   *   *           -|                              
#   9+  8+  7+  6+  5+  4+  3+  2+  1+  0+          -|           (image capillary is no longer present)
#   9+  8+  7+  6+  5+  4+  3+  2+  1+  0+          -|           
#   9+  8+  7+  6+  5+  4+  3+  2+  1+  0+          -|           
#   9+  8+  7+  6+  5+                              -|                               
#   9+  8+  7+                                      -|                                 
#
# <---------negative---------------------------------0--------------------positive---->  x-axis  
#                                       |                        |
#                                 -(x_cel_dist) +(x_cel_dist)#
#                                     
#
#
# CHANGING SYSTEM DIMENSIONS:
# 1. By default,
#    potentials, distances, fields, charges are reported in output files assuming that 1 a.u. length = 1 nm:
#    [V]         [nm]       [V/nm]  [e] 
#
# 2. Changing the ESI voltage: Could be done by running program with a specific cap_voltage value.
#                      BETTER: Always run with  cap_voltage = 1 (ESI voltage = 1 V).
#                      For scaling to different ESI voltage = X V: multiply potentials by X
#                                                                  multiply fields     by X
#                                                                  multiply charges    by X
#                                                                  leave distances unchanged               
#
# 3. Changing length units:
#    * Potentials [V] - no change needed (or multiply by X as explained below)
#
#    * Field = -delta_potential / delta_distance; switch from length unit 1 a.u. =    1 nm to another unit
#                                                            for example  1 a.u. = 0.01 mm = 1E4 nm  
#      This means that fields will retain their numerical values, but units switch from V/nm to V/(1E4 nm) = V/(0.01 mm).
#      In the example used here, 1 V/nm becomes V/(0.01 mm) = 100 V/mm
#      BOTTOM LINE: SIZE INCREASE BY 1E4 LOWERS FIELD BY 1E4
#
#    * Charges: Start from charges in e (generated for length in nm)
#                                                 switch from length unit 1 a.u. =    1 nm to another unit
#                                                            for example  1 a.u. = 0.01 mm
# 
#      potential = SUM (qj/rj); so for potential to stay constant, rj and qj must be scaled by the same factor 
#      In the example used here, 0.01 mm / 1 nm = 1E-5 m / 1E-9 m = 1E4, so charges have to be multiplied by 1E4 to retain the same potential 
#      The rescaled charges obtained in this way are still in e. To convert to C, multiply by 1.6E-19.
#      BOTTOM LINE: SIZE INCREASE BY 1E4 REQUIRES CHARGE INCREASE BY 1E4, IF POTENTIAL IS TO BE CONSTANT
# 
#
#
#
# INPUT FILE:
# ===========
#
# All the parameters required to perform a run are summarized in an input file
# Example (remove all # symbols when using this text as imput file):
#
#=== input file for capillary_potential.py === DO NOT CHANGE ORDER OF INPUT PARAMETERS ===
#===  output prefix; used for all files produced by the program                                  >output_prefix
#test_with_cel
#=== mage capillary and/or counter electrode (cel) present (yes or no)                           >cel
#yes
#===  dist of cap tip from cel (or from x = 0 if cel = "no"; same units as bead_radius_cap)      >x_cel_dist
#100.
#===  number determines counter electrode disk size                                              >n_cel_grid_size
#10
#===  distance between back wall and cap exit (0. = no back wall; same units as bead_radius_cap) >bwa_dist
#10
#===  capillary voltage set by power supply (V)                                                  >cap_voltage
#1.
#===  # of capillary fitting iterations                                                          >n_capfit_iter
#20
#===  # of cel fit iterations (0 = no charge transfer from img cap to cel)                       >n_celfit_iter
#30
#===  radius of capillary beads (a.u. or nm)                                                     >bead_radius_cap
#1.
#===  radius of counter electrode beads (same units as bead_radius_cap)                          >bead_radius_cel
#6.
#=== radius of back wall beads (same units as bead_radius_cap)                                   >bead_radius_bwa
#0.2
#===  number of beads in each cap_ring (comma delimited, start & end [] are optional, NO LINE BREAK ALLOWED)
#[10,10,10,10,10,10,10,10,10,10,10,10,10,10,12,14,16,18,20,22,24,26,28,30,30,30,30,30,30,30,30,30,30]
#
#
#
# OUTPUT FILES:
# =============
#
# The following output files can be expected (names preceded by output_prefix)
#
# ..._center_pot_V_field_Vnm_no_img_no_cel.txt
# potential and field for optimized charges, contribution from capillary only
#
# t..._center_pot_V_field_Vnm_with_cel.txt
# potential and field for optimized charges, contribution from capillary and cel
#
# ..._center_pot_V_field_Vnm_with_img.txt 
# potential and field for optimized charges, contribution from capillary and img capillary
#
# ..._fit_summary_cap.txt
# summary of all fitting iterations for capillary charges: energies, pot on each ring, charges on each ring
#
# test_with_cel_fit_summary_cel.txt
# same as above, but for cel
#
# ..._inp_outp_summary.txt
# summary of all input and output parameters. Also: capillary outline coordinates
#
# ..._rad_cel_pot_V.txt
# radial positions of cel rings vs optimized potentials
# 
# ..._rad_cel_pot_V_no_img_no_cel.txt
# as above, but only considering the contributions from the capillary
#
# ..._system_with_charges.gro
# xyz coordinates, charges in a.u., charges in e for all atoms
# Use this file to visualize the system, e.g. using Pymol or VMD
# 
# ..._xy_cap_charge_e.txt
# 2D data of capillary charges
#
# ..._xy_elfield_Vnm_dist_1_with_cel.txt
# electric field at capillary surface, in the presence of cel
# 
# ..._xy_elfield_Vnm_dist_1_with_img.txt
# same as above, but in the presence of image capillary
#
# ..._xy_pot_V_no_img_no_cel.txt
# 2D data of potential in the plance of the capillary centerline, img cap and cel contributions omitted
#
# ..._xy_pot_V_with_cel.txt
# as above, with cel included (no img cap)
# 
# ..._xy_pot_V_with_img.txt
# as above, with img cap included (no cel)



import sys
import numpy as np
import os.path
from os import path
import random


#------------------------------------------------------------------------------------------
# Start of function read_input():
# CHANGE PARAMETERS HERE

def read_input(input_file_name):

  global gro_out_name
  global n_capfit_iter
  global n_celfit_iter
  global x_cel_dist
  global n_bead_in_ring_cap
  global cel
  global n_cel_grid_size
  global bwa_dist
  global cap_voltage
  global output_prefix
  global bead_radius_cap
  global bead_radius_cel
  global bead_radius_bwa

  # read all lines of input_file_name
  input         = open(input_file_name,"r")         
  input_content = input.readlines()  
  input.close()

  # output prefix is used for all files produced by the program
  output_prefix   = input_content[2].rstrip()
  print("  output_prefix                                           : ",output_prefix)

  # image capillary and/or counter electrode present ("yes") or not ("no")
  cel             = input_content[4].rstrip()
  print("  counter electr / img cap present                        : ",cel)

  # distance of capillary tip from counter-electrode (or from x = 0 if cel = "no")
  x_cel_dist      = float(input_content[6])
  print("  cap tip dist from counter electr (or from x=0 if cel=no): ",x_cel_dist)

  # number determines counter electrode disk size
  n_cel_grid_size = int(input_content[8])
  print("  n_cel_grid_size (determines  size of the counter electr): ",n_cel_grid_size)

  # distance between back wall and cap exit (0. = no back wall)
  bwa_dist        = float(input_content[10])
  print("  dist between back wall and cap exit (0. = no back wall) : ",bwa_dist)

  # capillary voltage set by power supply (V)
  cap_voltage     = float(input_content[12])
  print("  capillary voltage (V)                                   : ",cap_voltage) 

  # number of capillary fit iterations
  n_capfit_iter   = int(input_content[14])
  print("  number of capillary fitting iterations                  : ",n_capfit_iter) 

  # number of counter electr fit iterations (0 = no charge transfer from img cap to counter electr)
  n_celfit_iter   = int(input_content[16])
  print("  number of counter electrode fitting iterations          : ",n_celfit_iter) 

 # radius of capillary beads; these beads are the building blocks of the capillary.
 # The units (a.u.) assumed for this parameter determine the dimensions of the capillary (it is ok to assume 1 a.u. = 1 nm).
 # It may be best to leave this at 1 or 0.1 (this facilitates subsequent changes to other length units, see program header)
  bead_radius_cap   = float(input_content[18])
  print("  radius of capillary beads                               : ",bead_radius_cap) 

  # radius of counter electrode beads (same units as bead_radius_cap)
  bead_radius_cel   = float(input_content[20])
  print("  radius of counter electrode beads                       : ",bead_radius_cel) 

  # radius of back wall beads (same units as bead_radius_cap; bead_radius_cap/5 might be a good value)
  bead_radius_bwa   = float(input_content[22])
  print("  radius of back wall beads                               : ",bead_radius_bwa)

  # number of beads in each cap_ring (comma delimited, square brackets mark beginning and end)
  n_bead_in_ring_cap_string = input_content[24]
  n_bead_in_ring_cap_string = n_bead_in_ring_cap_string.replace("[","")
  n_bead_in_ring_cap_string = n_bead_in_ring_cap_string.replace("]","")
#  n_bead_in_ring_cap = [int(number) for number in n_bead_in_ring_cap_string.split(",") if number.strip().isdigit()]
  n_bead_in_ring_cap = [int(number) for number in n_bead_in_ring_cap_string.split(",")]
  print("  number of beads in each ring                            : ",n_bead_in_ring_cap)

#------------------------------------------------------------------------------------------
# Calculating an xy - potential map at z = 0 (in the capillary center)
# for the optimized charge distribution.
# Also: writing data to output file

def calculate_xy_pot_map(rest_of_filename): 

  print("  xy mapping of potential  ...")

  xy_pot_name = output_prefix + rest_of_filename

  n_x             = 51                                # number of x data points
  n_y             = 51                                # number of y data points
  xy_pot_content  = []

  # scanning from negative x (left of capillary inlet) to zero (counter electrode plane)
  x_start              = bead_xyzq[n_bead_cap-1,0] *1.1    
  delta_x              = abs( x_start/ (n_x - 1) )
  
  #scanning from negative to positive y
  y_start              = 0.5 * bead_xyzq[n_bead_cap-1,0] *1.1    
  delta_y              = abs( 2 * y_start/ (n_y - 1) )

  for i_x in range(0,n_x):
    x = x_start + i_x * delta_x  

    for i_y in range(0,n_y):
      y = y_start + i_y * delta_y
    
      pot = 0.

      for j_bead in range(0,2*n_bead_cap + n_bead_cel):

        d_ij = np.sqrt(   (x - bead_xyzq[j_bead,0])**2
                        + (y - bead_xyzq[j_bead,1])**2
                        + (    bead_xyzq[j_bead,2])**2  )

        if (d_ij < 0.5*bead_radius_cap): d_ij = 0.5*bead_radius_cap 

        # pot = pot + bead_xyzq[j_bead,3] / d_ij  # this is in a.u.

        pot = pot + bead_xyzq[j_bead,3] / d_ij * volt_conv_factor   # Potential in Volts


      line_out =str( format(x,"12.5f") + format(y,"12.5f") + format(pot,"12.5f") + "\n" )
      xy_pot_content.append(line_out)

  xy_pot_out = open(xy_pot_name,"w+")
  xy_pot_out.writelines(xy_pot_content)
  xy_pot_out.close()
  print("  ... saved as :",xy_pot_name)
  print("")

#------------------------------------------------------------------------------------------
# Calculating xy - electric field map at z = 0 (in the capillary center plane)
# for the optimized charge distribution.
# One could use
# E_x = -(pot[x+dx,y] - pot[x,y])/dx
# E_y = -(pot[x,y+dy] - pot[x,y])/dy ... but this would yield the potential slightly next to point x/y.
# Better:
# E_x = -(pot[x+dx,y] - pot[x-dx,y])/(2.*dx)
# E_y = -(pot[x,y+dy] - pot[x,y-dy])/(2.*dy)  ... to truly get the derivative AT point x/y.
#
# Also: writing data to output file: x  y  E_x [V/nm]  E_y [V/nm]  angle [degrees rel. to x-axis]  field_strength [V/nm]
# 
# dist_factor determines how far from the capillary the scan is performed

def calculate_xy_elfield_map(rest_of_filename,dist_factor): 

  print("  xy mapping of electric field  ...")

  xy_elfield_name = output_prefix + rest_of_filename

  xy_elfield_content  = []

  # generating xy points that trace the capillary outline
  # 
  #      -->  -->  -->  --> 
  #      4   3   2   1   0    *
  #      4   3   2   1   0  
  #      4   3   2   1   0    *  
  #      4   3   2   1   0
  #      4   3   2   1   0    *
  #      <--  <--  <--  <--                   

  i_ring_cap = n_ring_cap

  for i_xy_point in range(0,2*n_ring_cap+5):
    if(i_xy_point < n_ring_cap  ): i_ring_cap = i_ring_cap - 1
    if(i_xy_point > n_ring_cap+5): i_ring_cap = i_ring_cap + 1
    x = cap_ring_xcoord[i_ring_cap]
    if(i_xy_point < n_ring_cap+5): y =     cap_ring_radius[i_ring_cap] + dist_factor*bead_radius_cap
    if(i_xy_point > n_ring_cap+5): y = -1.*cap_ring_radius[i_ring_cap] - dist_factor*bead_radius_cap

    if(i_xy_point == n_ring_cap):
      x = cap_ring_xcoord[0]       + dist_factor*bead_radius_cap
      y = cap_ring_radius[0]       + dist_factor*bead_radius_cap

    if(i_xy_point == n_ring_cap+1):
      x = cap_ring_xcoord[0]       + dist_factor*bead_radius_cap
      y = 0.5*(cap_ring_radius[0]  + dist_factor*bead_radius_cap)

    if(i_xy_point == n_ring_cap+2):
      x = cap_ring_xcoord[0]       + dist_factor*bead_radius_cap
      y = 0.

    if(i_xy_point == n_ring_cap+3):
      x =       cap_ring_xcoord[0] + dist_factor*bead_radius_cap
      y = -0.5*(cap_ring_radius[0] + dist_factor*bead_radius_cap)

    if(i_xy_point == n_ring_cap+4):
      x = cap_ring_xcoord[0]      + dist_factor*bead_radius_cap
      y = -1.*cap_ring_radius[0]  - dist_factor*bead_radius_cap

    dx = 0.05 * bead_radius_cap
    dy = dx

    
    pot_xmdx_y = 0.   #potential at point x-dx,y
    pot_xpdx_y = 0.   #potential at point x+dx,y
    pot_x_ymdy = 0.   #potential at point x,y-dy
    pot_x_ypdy = 0.   #potential at point x,y+dy

    for j_bead in range(0,2*n_bead_cap + n_bead_cel):

      d_ij_xmdx_y = np.sqrt(   (x-dx - bead_xyzq[j_bead,0])**2  #for calculating potential at point x-dx,y
                          +    (y    - bead_xyzq[j_bead,1])**2
                          +    (       bead_xyzq[j_bead,2])**2  )

      d_ij_xpdx_y = np.sqrt(   (x+dx - bead_xyzq[j_bead,0])**2  #for calculating potential at point x+dx,y
                             + (y    - bead_xyzq[j_bead,1])**2
                             + (       bead_xyzq[j_bead,2])**2  )


      d_ij_x_ymdy = np.sqrt(   (x    - bead_xyzq[j_bead,0])**2  #for calculating potential at point x,y-dy
                             + (y-dy - bead_xyzq[j_bead,1])**2
                             + (       bead_xyzq[j_bead,2])**2  )

      d_ij_x_ypdy = np.sqrt(   (x    - bead_xyzq[j_bead,0])**2  #for calculating potential at point x,y+dy
                             + (y+dy - bead_xyzq[j_bead,1])**2
                             + (       bead_xyzq[j_bead,2])**2  )

      if (d_ij_xmdx_y < 0.5*bead_radius_cap): d_ij_xmdx_y = 0.5*bead_radius_cap 
      if (d_ij_xpdx_y < 0.5*bead_radius_cap): d_ij_xpdx_y = 0.5*bead_radius_cap
      if (d_ij_x_ymdy < 0.5*bead_radius_cap): d_ij_x_ymdy = 0.5*bead_radius_cap  
      if (d_ij_x_ypdy < 0.5*bead_radius_cap): d_ij_x_ypdy = 0.5*bead_radius_cap 
        
      # these potentials are in a.u.
      pot_xmdx_y = pot_xmdx_y + bead_xyzq[j_bead,3] / d_ij_xmdx_y
      pot_xpdx_y = pot_xpdx_y + bead_xyzq[j_bead,3] / d_ij_xpdx_y
      pot_x_ymdy = pot_x_ymdy + bead_xyzq[j_bead,3] / d_ij_x_ymdy   
      pot_x_ypdy = pot_x_ypdy + bead_xyzq[j_bead,3] / d_ij_x_ypdy

    # potentials are converted to V (assuming distances are in nm), so that field is in V/nm
    elfield_x     = -1.*( pot_xpdx_y - pot_xmdx_y ) * volt_conv_factor / (2.*dx)  #electric field component in x-direction
    elfield_y     = -1.*( pot_x_ypdy - pot_x_ymdy ) * volt_conv_factor / (2.*dy)  #electric field component in y-direction
        

    elfield_norm  = np.sqrt(elfield_x**2 + elfield_y**2)   #length of elfield vector

    # angle of  elfield vector relative to x-axis (converted to degrees)
    # (parallel = 0 deg; up = 90 deg; antiparallel = -180)
      
    if (elfield_y >= 0.) :  # angles 0 to 180 deg
      elfield_angle =        np.arccos(elfield_x/elfield_norm)*360./2/3.14159

    if (elfield_y <  0.) :  # angles 0 to -180 deg
      elfield_angle =    -1.*np.arccos(elfield_x/elfield_norm)*360./2/3.14159


    line_out =str( format(x,"12.5f")             + format(y,"12.5f") 
                 + format(elfield_x,"12.5f")     + format(elfield_y,"12.5f") 
                 + format(elfield_angle,"12.5f") + format(elfield_norm,"12.5f") + "\n" )
    xy_elfield_content.append(line_out)

  xy_elfield_out = open(xy_elfield_name,"w+")
  xy_elfield_out.writelines(xy_elfield_content)
  xy_elfield_out.close()
  print("  ... saved as :",xy_elfield_name)
  print("")

#------------------------------------------------------------------------------------------
# writing overall input and some results to output file

def write_input_output_summary(): 

  print("  writing summary of input and output parameters ...")

  input_output_summary_name = output_prefix + "_inp_outp_summary.txt"

  input_output_summary_content = []

  output_line = "  === capillary_potential.py parameters === " + "\n" 
  input_output_summary_content.append(output_line)

  output_line = "  " + "\n" 
  input_output_summary_content.append(output_line)

  output_line = "  INPUT  ****************" + "\n" 
  input_output_summary_content.append(output_line)

  output_line = "  output_prefix     = " + output_prefix + "\n"
  input_output_summary_content.append(output_line)    

  output_line = "  n_bead_in_ring_cap= " + str(n_bead_in_ring_cap) + "\n" 
  input_output_summary_content.append(output_line)
                     
  output_line = "  n_capfit_iter     = " + str(n_capfit_iter) + "\n"
  input_output_summary_content.append(output_line)  

  output_line = "  n_celfit_iter     = " + str(n_celfit_iter) + "\n"
  input_output_summary_content.append(output_line) 

  output_line = "  cel               = " + cel + "\n"
  input_output_summary_content.append(output_line) 

  output_line = "  bwa_dist          = " + str(bwa_dist) + "\n"
  input_output_summary_content.append(output_line) 

  output_line = "  x_cel_dist        = " + str(x_cel_dist) + "\n"
  input_output_summary_content.append(output_line) 

  output_line = "  n_cel_grid_size   = " + str(n_cel_grid_size) + "\n"
  input_output_summary_content.append(output_line) 

  output_line = "  cap_voltage       = " + str(cap_voltage) + "\n"
  input_output_summary_content.append(output_line) 

  output_line = "  bead_radius_cap   = " + str(bead_radius_cap) + "\n"
  input_output_summary_content.append(output_line)

  output_line = "  bead_radius_cel   = " + str(bead_radius_cel) + "\n"
  input_output_summary_content.append(output_line)

  output_line = "  bead_radius_bwa   = " + str(bead_radius_bwa) + "\n"
  input_output_summary_content.append(output_line)

  output_line = "  " + "\n" 
  input_output_summary_content.append(output_line)

  output_line = "  OUTPUT **************** " + "\n" 
  input_output_summary_content.append(output_line)

  output_line = "  capillary inlet x-position [units]    = " + format(bead_xyzq[n_bead_cap-1,0],"12.3f") + "\n"
  input_output_summary_content.append(output_line)

  cap_length = abs(bead_xyzq[n_bead_cap-1,0] + x_cel_dist)
  output_line = "  capillary length [units]              = " + format(cap_length,"12.3f") + "\n"
  input_output_summary_content.append(output_line)

  output_line = "  capillary outlet i.d. [units]         = " + format(2*cap_ring_radius[0],"12.3f") + "\n"
  input_output_summary_content.append(output_line)      

  output_line = "  capillary  inlet i.d. [units]         = " + format(2*cap_ring_radius[n_ring_cap-1],"12.3f") + "\n"
  input_output_summary_content.append(output_line) 

  # Determining maximum system dimensions (largest bead-bead distance, excl. image beads).
  # Also determining counter electrode diameter.

  max_sys_dist = 0.
  cel_diam      = 0.
  for i_bead in range(0,2*n_bead_cap + n_bead_cel):
    for j_bead in range(i_bead+1,2*n_bead_cap + n_bead_cel):

      dist = np.sqrt(  (bead_xyzq[i_bead,0] - bead_xyzq[j_bead,0])**2
                     + (bead_xyzq[i_bead,1] - bead_xyzq[j_bead,1])**2
                     + (bead_xyzq[i_bead,2] - bead_xyzq[j_bead,2])**2  )

      if (dist > max_sys_dist):
        if (cel == "yes"):
          # by excluding positive x-values, image beads are excluded
          if (bead_xyzq[i_bead,0] <= 0.) and (bead_xyzq[j_bead,0] <= 0.):max_sys_dist = dist 
        if (cel == "no"):
          # by excluding positive x-values and x = 0, image and cel beads are excluded
          # (this works properly only for x_cel_dist > 0, otherwise the first ring is not considered)
          if (bead_xyzq[i_bead,0] < 0.) and (bead_xyzq[j_bead,0] < 0.):max_sys_dist = dist 

      if (dist > cel_diam):
        if (cel == "yes"):
          # by focusing on x = 0 values, only counter electrode beads are considered
          if (bead_xyzq[i_bead,0] == 0.) and (bead_xyzq[j_bead,0] == 0.):cel_diam = dist 


  output_line = "  maxdist. (cap and cel beads) [units]  = " + format(max_sys_dist,"12.3f") + "\n"
  input_output_summary_content.append(output_line) 

  output_line = "  n_ring_cap                            =           " + str(n_ring_cap) + "\n"
  input_output_summary_content.append(output_line) 

  output_line = "  n_ring_cel                            =           " + str(n_ring_cel) + "\n"
  input_output_summary_content.append(output_line) 

  output_line = "  n_bead_cap                            =           " + str(n_bead_cap) + "\n"
  input_output_summary_content.append(output_line) 

  output_line = "  n_bead_cel                            =           " + str(n_bead_cel) + "\n"
  input_output_summary_content.append(output_line)

  output_line = "  counter electrode diameter [units]    = " + format(cel_diam,"12.3f") + "\n"
  input_output_summary_content.append(output_line) 

  output_line = "  charge recovery on counter electrode  = " + format(charge_recov_on_cel,"12.3f") + "\n"
  input_output_summary_content.append(output_line)

  #save coordinates of capillary outline

  output_line = "\n"
  input_output_summary_content.append(output_line)

  output_line = "  capillary x/y outline [units]" + "\n"
  input_output_summary_content.append(output_line)

  for i_ring_cap in range (n_ring_cap-1,-1,-1):
    output_line = format(cap_ring_xcoord[i_ring_cap],"12.3f") + format(cap_ring_radius[i_ring_cap],"12.3f") + "\n"
    input_output_summary_content.append(output_line)

  for i_ring_cap in range (0,n_ring_cap):
    output_line = format(cap_ring_xcoord[i_ring_cap],"12.3f") + format(-1*cap_ring_radius[i_ring_cap],"12.3f") + "\n"
    input_output_summary_content.append(output_line)

  output_line = format(cap_ring_xcoord[n_ring_cap-1],"12.3f") + format(cap_ring_radius[n_ring_cap-1],"12.3f") + "\n"
  input_output_summary_content.append(output_line)



  line_out = "\n" 
  input_output_summary_content.append(line_out)
  line_out = "CAPILLARY CHARGES"  + "\n" 
  input_output_summary_content.append(line_out)
  line_out = " ring#     #beads  x-coord (units)   q per bead (a.u.)   q per bead (e)"  + "\n"
  input_output_summary_content.append(line_out)

  for i_ring_cap in range(0,n_ring_cap):

    # charge values prior to unit conversion:
    # These values were used internally during fitting, assuming that every capillary
    # bead initially had a charge of 1 a.u.
    charge_in_au   = charge_per_bead_in_ring_cap[i_ring_cap]
    
    #charge values after unit conversion
    charge_in_e = charge_per_bead_in_ring_cap[i_ring_cap] / 1.44 * volt_conv_factor
    line_out =str( format(i_ring_cap                    ,"6") 
                 + format(n_bead_in_ring_cap[i_ring_cap],"6") 
                 + format(cap_ring_xcoord[i_ring_cap]   ,"16.7f")
                 + format(charge_in_au                  ,"16.7f")
                 + format(charge_in_e                   ,"20.7f") + "\n" )
    input_output_summary_content.append(line_out)

  line_out = "                         "  + "\n" 
  input_output_summary_content.append(line_out)

  line_out = "COUNTER ELECTRODE CHARGES (data are listed twice (first with neg. radial_dist) for plotting)"  + "\n" 
  input_output_summary_content.append(line_out)
  line_out = " ring#     #beads  rad dist (units)  q per bead (a.u.)   q per bead (e)"  + "\n"
  input_output_summary_content.append(line_out)

  for i_ring_cel in range(n_ring_cel-1,0,-1):                # printing in reverse order, from n_ring_cel downward

    charge_in_au   = charge_per_bead_in_ring_cel[i_ring_cel] # nalogous to cap charges in previous section

    charge_in_e = charge_per_bead_in_ring_cel[i_ring_cel] / 1.44 * volt_conv_factor
    line_out =str( format(i_ring_cel                            ,"6") 
                 + format(n_bead_in_ring_cel[i_ring_cel]        ,"6")
                 + format(-1.*cel_ring_radius[i_ring_cel]       ,"16.7f")
                 + format(charge_in_au                          ,"16.7f") 
                 + format(charge_in_e                           ,"20.7f") + "\n" )
    input_output_summary_content.append(line_out)


  for i_ring_cel in range(0,n_ring_cel):                     # printing again, in proper order from 0 upward

    charge_in_au   = charge_per_bead_in_ring_cel[i_ring_cel] # analogous to cap charges in previous section

    charge_in_e = charge_per_bead_in_ring_cel[i_ring_cel] / 1.44 * volt_conv_factor
    line_out =str( format(i_ring_cel                        ,"6") 
                 + format(n_bead_in_ring_cel[i_ring_cel]    ,"6")
                 + format(cel_ring_radius[i_ring_cel],       "16.7f")
                 + format(charge_in_au                      ,"16.7f") 
                 + format(charge_in_e                       ,"20.7f") + "\n" )
    input_output_summary_content.append(line_out)



  input_output_summary_out = open(input_output_summary_name,"w+")
  input_output_summary_out.writelines(input_output_summary_content)
  input_output_summary_out.close()
  print("  ... saved as :",input_output_summary_name)
  print("")

#------------------------------------------------------------------------------------------
# writing capillary fitting statistics to output file

def write_fit_summary_cap(rest_of_filename): 

  print("  writing capillary fitting summary ...")

  fit_summary_cap_name = output_prefix + rest_of_filename

  fit_summary_cap_content = []

  for i_ring_cap in range(0,n_ring_cap+1):
    output_line = str(format(i_ring_cap-1,"6"))
    if (i_ring_cap == 0): output_line = " ring#"
    for j_capfit_iter in range(0,n_capfit_iter+1):
      output_line = output_line + str( format(fit_summary_cap[i_ring_cap,2*j_capfit_iter  ],"11.2f") 
                                     + format(fit_summary_cap[i_ring_cap,2*j_capfit_iter+1],"11.2f")) 
    output_line = output_line+ "\n" 
    fit_summary_cap_content.append(output_line)


  fit_summary_cap_out = open(fit_summary_cap_name,"w+")
  fit_summary_cap_out.writelines(fit_summary_cap_content)
  fit_summary_cap_out.close()
  print("  ... saved as :",fit_summary_cap_name)
  print("")


#------------------------------------------------------------------------------------------
# writing counter electrode fitting summary to output file

def write_fit_summary_cel(rest_of_filename): 

  print("  writing counter electrode fitting summary ...")

  fit_summary_cel_name = output_prefix + rest_of_filename

  fit_summary_cel_content = []

  for i_ring_cel in range(0,n_ring_cel+1):
    output_line = str(format(i_ring_cel-1,"6"))
    if (i_ring_cel == 0): output_line = " ring#"
    for j_celfit_iter in range(0,n_celfit_iter+1):
      output_line = output_line + str( format(fit_summary_cel[i_ring_cel,2*j_celfit_iter  ],"11.2f") 
                                     + format(fit_summary_cel[i_ring_cel,2*j_celfit_iter+1],"11.2f")) 
    output_line = output_line+ "\n" 
    fit_summary_cel_content.append(output_line)


  fit_summary_cel_out = open(fit_summary_cel_name,"w+")
  fit_summary_cel_out.writelines(fit_summary_cel_content)
  fit_summary_cel_out.close()
  print("  ... saved as :",fit_summary_cel_name)
  print("")

#------------------------------------------------------------------------------------------
# Calculating an xy - charge map at z = 0 (in the capillary center)
# for the optimized capillary charges (no other charges are considered here).

def calculate_xy_cap_charge_map(rest_of_filename): 

  print("  xy mapping of capillary charge ...")

  xy_cap_charge_name = output_prefix + rest_of_filename

  n_x                    = n_ring_cap+10      # number of x data points
  n_y                    = n_ring_cap+5       # number of y data points
  xy_cap_charge_content  = []

  # finding minimum y value for grid sizing
  y_min = 0.
  for j_bead in range(0,n_bead_cap):
    if (bead_xyzq[j_bead,1] < y_min): y_min = bead_xyzq[j_bead,1]
  y_min = y_min * 1.1

  # scanning from negative x (left of cap inlet) to just past the cap outlet
  x_start          = bead_xyzq[n_bead_cap-1,0] - 10. * 0.866 * bead_radius_cap  
  delta_x          = 0.866 * 2. * bead_radius_cap
  
  #scanning from negative to positive y
  y_start          = y_min - 2. * 0.866 * bead_radius_cap   
  delta_y          = 2.* (abs(y_min - 2. * 0.866 * bead_radius_cap) / (n_y-1) )

  for i_x in range(0,n_x):
    x = x_start + i_x * delta_x  

    for i_y in range(0,n_y):
      y = y_start + i_y * delta_y
    
      charge = 0.

      for j_bead in range(0,n_bead_cap):
         
        # distances in xy projected capillary  
        d_ij = np.sqrt(   (x - bead_xyzq[j_bead,0])**2
                        + (y - bead_xyzq[j_bead,1])**2  )
        # charges are converted to e, assuming distances are in nm
        if (d_ij < bead_radius_cap): charge = bead_xyzq[j_bead,3] / 1.44 * volt_conv_factor


      line_out =str( format(x,"12.5f") + format(y,"12.5f") + format(charge,"12.5f") + "\n" )
      xy_cap_charge_content.append(line_out)

  xy_cap_charge_out = open(xy_cap_charge_name,"w+")
  xy_cap_charge_out.writelines(xy_cap_charge_content)
  xy_cap_charge_out.close()
  print("  ... saved as :",xy_cap_charge_name)
  print("")

#------------------------------------------------------------------------------------------
# Calculating a radial potential map of the counter electrode at x = -0.5 * bead_radius_cel 
#                                                                y = 0
# (all charges are considered here). 

def calculate_rad_cel_pot_map(rest_of_filename): 

  print("  mapping radial counter electrode potential ...")

  rad_cel_pot_name = output_prefix + rest_of_filename

  n_z                  = 101     # number of z data points (radial scan is performed along z-axis)
  rad_cel_pot_content  = []

  # finding minimum y value for grid sizing
  z_min = 0.
  for j_bead in range(2*n_bead_cap,2*n_bead_cap + n_bead_cel):
    if (bead_xyzq[j_bead,2] < z_min): z_min = bead_xyzq[j_bead,2]
  z_min = z_min * 1.1

   
  delta_z = 2.* abs(z_min / (n_z-1) )

  for i_z in range(0,n_z):
    z = z_min + i_z * delta_z  

    pot = 0.

    for j_bead in range(0,2*n_bead_cap + n_bead_cel):
       
      # xyz distances in counter electrode
      # Scanning is performed not IN the xy plane (x = 0), but with a slight x-offset to avoid singularities.  
      #                  <--probe coordinates--> <--bead coordinates-->
      d_ij = np.sqrt(   (-0.5 * bead_radius_cel  - bead_xyzq[j_bead,0])**2
                      + (0.                      - bead_xyzq[j_bead,1])**2
                      + (z                       - bead_xyzq[j_bead,2])**2  )

      if (d_ij < bead_radius_cel): d_ij = bead_radius_cel  # for smoothing the scan 

      pot = pot + bead_xyzq[j_bead,3]/d_ij      # pot = SUM(q_j_d_ij)

    pot = pot * volt_conv_factor   # Potential is now in Volts

    line_out = str( format(z,"12.5f") + format(pot,"12.5f") + "\n" )
    rad_cel_pot_content.append(line_out)

  rad_cel_pot_out = open(rad_cel_pot_name,"w+")
  rad_cel_pot_out.writelines(rad_cel_pot_content)
  rad_cel_pot_out.close()
  print("  ... saved as :",rad_cel_pot_name)
  print("")


#------------------------------------------------------------------------------------------
# Calculating potential and field along the center of the capillary (= x axis)
# for the optimized charge distribution, and writting data to output file

def calculate_center_pot_field(rest_of_filename):

  global volt_conv_factor # This is also needed for generate_charge_output()
                          # and                     calculate_rad_cel_pot_map()
                          # and                     calculate_xy_cap_charge_map


  print("  calculating potential along center (x-axis) ...")

  output_name = output_prefix + rest_of_filename                                    
                                            

  x_start              = bead_xyzq[n_bead_cap-1,0]*1.1   # we scan from here  
  delta_x              = abs( bead_radius_cap/5. )       # scanning increment
  n_x                  = 2* int(abs(x_start/delta_x)+1)  # number of data points
  center_pot_field     = np.zeros((n_x,3))               # center_pot_field[i_x,0] = x position
  center_pot_av        = 0.                              # center_pot_field[i_x,1] = potential
  n_points_within_cap  = 0                               # center_pot_field[i_x,2] = electric field
  center_out_content   = []

  center_out_content.append("   x(a.u.)    pot(V)       field(V/[a.u.])  along x-axis"  + "\n" )

  # this is the x-position in the middle of the capillary
  x_middle_cap = bead_xyzq[n_bead_cap-1,0] + 0.5 * (bead_xyzq[0,0]-bead_xyzq[n_bead_cap-1,0])
  #print("  x_middle_cap = ",x_middle_cap)
  min_x_dist_from_x_middle_cap = 1E40
  x_middle_pot                 = 0.

  for i_x in range(0,n_x):
    # defining x values for scan line
    center_pot_field[i_x,0] = x_start + i_x * delta_x

    # Calculate the electric potential at each point along the x-axis
    # as     potential = SUM(all j_bead) q_j/d_ij    where q_j [elemtary charges] and d_ij [nm]
    # and the sum includes all beads and image beads

    for j_bead in range(0,2*n_bead_cap + n_bead_cel):

      d_ij = np.sqrt(   (center_pot_field[i_x,0] - bead_xyzq[j_bead,0])**2
                      + (                          bead_xyzq[j_bead,1])**2
                      + (                          bead_xyzq[j_bead,2])**2  )

      if (d_ij < 0.5 * bead_radius_cap):d_ij = 0.5 * bead_radius_cap
      center_pot_field[i_x,1] = center_pot_field[i_x,1] + bead_xyzq[j_bead,3] / d_ij

    # Finding the potential value (in a.u.) in the middle of the capillary
    test_x_dist_from_x_middle_cap = abs(center_pot_field[i_x,0] - x_middle_cap)
    if(test_x_dist_from_x_middle_cap <= min_x_dist_from_x_middle_cap):
      min_x_dist_from_x_middle_cap = test_x_dist_from_x_middle_cap
      x_middle_pot = center_pot_field[i_x,1]

  # Converting potentials to Volts, such that x_middle_pot is equal to the user-defined cap_voltage
  volt_conv_factor = cap_voltage / x_middle_pot

  # converting centerline potentials to Volts
  for i_x in range(0,n_x):
    center_pot_field[i_x,1] = center_pot_field[i_x,1]*volt_conv_factor   # Potential is now in Volts

  # calculating centerline electric field
  for i_x in range(0,n_x):
    if (   (i_x > 0) and ( i_x < (n_x-1) )   ):
      # Electric field = E_x = -d_potential/d_x
      # First and last E_x values are skipped in this loop
      center_pot_field[i_x,2] =  -1. * (center_pot_field[i_x+1,1] - center_pot_field[i_x-1,1] ) / (2 * delta_x)
                                                                 
  center_pot_field[0    ,2] = center_pot_field[1    ,2]          # First and last E_x values are 
  center_pot_field[n_x-1,2] = center_pot_field[n_x-2,2]          # assigned their neighbor values.


  # write center output 

  for i_x in range(0,n_x):
    line_out =str( format(center_pot_field[i_x,0],"12.5f")  
                 + format(center_pot_field[i_x,1],"12.5f") 
                 + format(center_pot_field[i_x,2],"12.5f") + "\n" )
    center_out_content.append(line_out)

  center_out = open(output_name,"w+")
  center_out.writelines(center_out_content)
  center_out.close()
  print("  ... saved as : ",output_name)
  print("") 


#------------------------------------------------------------------------------------------
# Generating .gro output file:  x/y/z coordinates, charge in a.u.; charge in e
#                                                                  (applies if length is in nm)
# Species included        1CAP capillary beads
#                         2IMG image capillary beads
#                         3CEL counter electrode beads
#                         4BWA back wall beads
#                         5CPO probe locations on outside    for measuring capillary potential
#                         5CPI probe locations on inside     for measuring capillary potential
#                         6EPP probe locations at positive x for measuring cel potential
#                         6EPP probe locations at negative x for measuring cel potential

def generate_gro_output(rest_of_filename):

  print("  saving coordinates in .gro format ...")

  gro_out_name = output_prefix + rest_of_filename

  gro_out_content    = [] # will be written to                .gro output
  bead_xyzq_text     = [] # capillary bead    coordinates for .gro output
  cap_probe_xyz_text = [] # cap_probe     coordinates for .gro output
  cel_probe_xyz_text = [] # cel_probe         coordinates for .gro output
  cel_xyz_text       = [] # counter electrode coordinates for .gro output
  bwa_xyz_text       = [] # back wall         coordinates for .gro output   


  # transferring capillary coordinates
  i_bead          = 0                
  for i_ring_cap in range(0,n_ring_cap):                       
    for i_bead_in_ring_cap in range(0,n_bead_in_ring_cap[i_ring_cap]):

      bead_type     = "C" + str(i_ring_cap).zfill(3)        # cap_ring number is coded in bead_type
                                                            # C for "Capillary"
      # converting charge to e for the chosen cap_voltage (assuming length units are nm)
      charge_in_e = bead_xyzq[i_bead,3] / 1.44 * volt_conv_factor
      bead_xyzq_text.append(str( "    1CAP   " + bead_type + format(i_bead,"5") 
                                                           + format(bead_xyzq[i_bead,0],"8.3f")  
                                                           + format(bead_xyzq[i_bead,1],"8.3f") 
                                                           + format(bead_xyzq[i_bead,2],"8.3f")
                                                           + format(bead_xyzq[i_bead,3],"8.3f") 
                                                           + format(charge_in_e        ,"8.3f") 
                                                           + "\n") )
      i_bead = i_bead + 1

  # transferring capillary probe coordinates 
  for i_cap_probe in range(0,n_ring_cap):                      
                                                            # P for capillary "P"robe, O for "outside"
    cap_probe_xyz_text.append(str( "    5CPO      P"  + format(i_bead,"5") 
                                                      + format( cap_probe_out_xyz[i_cap_probe,0],"8.3f")  
                                                      + format( cap_probe_out_xyz[i_cap_probe,1],"8.3f") 
                                                      + format( cap_probe_out_xyz[i_cap_probe,2],"8.3f")
                                                      + format(0.                               ,"8.3f") 
                                                      + format(0.                               ,"8.3f") 
                                                      + "\n") )

  for i_cap_probe in range(0,n_ring_cap):                      
                                                            # P for capillary "P"robe, I for "inside"
    cap_probe_xyz_text.append(str( "    5CPI      P"  + format(i_bead,"5") 
                                                      + format( cap_probe_ins_xyz[i_cap_probe,0],"8.3f")  
                                                      + format( cap_probe_ins_xyz[i_cap_probe,1],"8.3f") 
                                                      + format( cap_probe_ins_xyz[i_cap_probe,2],"8.3f")
                                                      + format(0.                               ,"8.3f") 
                                                      + format(0.                               ,"8.3f") 
                                                      + "\n") )
  
  # transferring image capillary coordinates
  i_bead     = n_bead_cap
  for i_ring_cap in range(0,n_ring_cap):                       
    for i_bead_in_ring_cap in range(0,n_bead_in_ring_cap[i_ring_cap]):

      bead_type     = "I" + str(i_ring_cap).zfill(3)        # cap_ring number is coded in bead_type
                                                            # I for "Image"
      bead_xyzq_text.append(str( "    2IMG   " + bead_type + format(i_bead,"5") 
                                                           + format(bead_xyzq[i_bead,0],"8.3f")  
                                                           + format(bead_xyzq[i_bead,1],"8.3f") 
                                                           + format(bead_xyzq[i_bead,2],"8.3f")
                                                           + format(0.                 ,"8.3f") 
                                                           + format(0.                 ,"8.3f") 
                                                           + "\n") )
      i_bead = i_bead + 1


  # transferring counter electrode coordinates
  i_bead_cel = 0
  for i_ring_cel in range(0,n_ring_cel)              :         
    for i_bead_in_ring_cel in range(0,n_bead_in_ring_cel[i_ring_cel]):
      bead_type     = "N" + str(i_ring_cel).zfill(3)        # i_ring_cel is coded in bead_type
                                                            # N for "couNter electrode"
      charge_in_e = bead_xyzq[i_bead,3] / 1.44 * volt_conv_factor
      cel_xyz_text.append(str( "    3CEL   " 
                                              + bead_type + format(i_bead,"5") 
                                                          + format(bead_xyzq[i_bead,0],"8.3f")  
                                                          + format(bead_xyzq[i_bead,1],"8.3f") 
                                                          + format(bead_xyzq[i_bead,2],"8.3f")
                                                          + format(bead_xyzq[i_bead,3],"8.3f")
                                                          + format(charge_in_e        ,"8.3f")  
                                                          + "\n") )
      i_bead = i_bead + 1

  # transferring back wall coordinates.
  if (bwa_dist  != 0.):
    for i_bead_bwa in range(0,n_bead_bwa):         
      bwa_xyz_text.append(str( "    4BWA   " 
                                              + bead_type + format(i_bead,"5") 
                                                          + format(bwa_xyz[i_bead_bwa,0],"8.3f")  
                                                          + format(bwa_xyz[i_bead_bwa,1],"8.3f") 
                                                          + format(bwa_xyz[i_bead_bwa,2],"8.3f")
                                                          + format(0.                   ,"8.3f") 
                                                          + format(0.                   ,"8.3f") 
                                                          + "\n") )
      i_bead = i_bead + 1


  # transferring counter electrode probe coordinates 
  for i_cel_probe in range(0,n_ring_cel):                      
                                                            # B for counter electrode Brobe
    cel_probe_xyz_text.append(str( "    6EPP      B"  + format(i_bead,"5") 
                                                      + format( cel_probe_beh_xyz[i_cel_probe,0],"8.3f")  
                                                      + format( cel_probe_beh_xyz[i_cel_probe,1],"8.3f") 
                                                      + format( cel_probe_beh_xyz[i_cel_probe,2],"8.3f")
                                                      + format(0.                               ,"8.3f")
                                                      + format(0.                               ,"8.3f")  
                                                      + "\n") )

    cel_probe_xyz_text.append(str( "    6EPN      B"  + format(i_bead,"5") 
                                                      + format( cel_probe_fro_xyz[i_cel_probe,0],"8.3f")  
                                                      + format( cel_probe_fro_xyz[i_cel_probe,1],"8.3f") 
                                                      + format( cel_probe_fro_xyz[i_cel_probe,2],"8.3f")
                                                      + format(0.                               ,"8.3f")
                                                      + format(0.                               ,"8.3f")  
                                                      + "\n") )

  # Building .gro output
  gro_out_content.append("coordinates and charges: x/y/z/q. [charge] in a.u. / e (for nm) " + "\n")

  n_bead_in_gro = n_bead_cap + 2*n_ring_cap
  if(cel == "yes"):         n_bead_in_gro = n_bead_in_gro + n_bead_cap
  if(cel == "yes"):         n_bead_in_gro = n_bead_in_gro + n_bead_cel + 2*n_ring_cel
  if(bwa_dist != 0.): n_bead_in_gro = n_bead_in_gro + n_bead_bwa

  gro_out_content.append(str(n_bead_in_gro) + "\n")

  for i_bead in range(0,n_bead_cap):                          # writing capillary bead coordinates
    gro_out_content.append(bead_xyzq_text[i_bead])


  if(cel == "yes"):
    for i_bead in range(n_bead_cap,2*n_bead_cap):             # writing image bead coordinates
      gro_out_content.append(bead_xyzq_text[i_bead])

  if(cel == "yes"):
    for i_bead_cel in range(0,n_bead_cel):                    # writing counter electrode coordinates
      gro_out_content.append(cel_xyz_text[i_bead_cel])

  if(bwa_dist != 0.):
    for i_bead_bwa in range(0,n_bead_bwa):                    # writing back wall coordinates
      gro_out_content.append(bwa_xyz_text[i_bead_bwa])
 
  for i_ring_cap in range(0,2*n_ring_cap):                    # writing cap_probe coordinates
    gro_out_content.append(cap_probe_xyz_text[i_ring_cap])

  if(cel == "yes"):
    for i_ring_cel in range(0,2*n_ring_cel):                  # writing cel_probe coordinates
      gro_out_content.append(cel_probe_xyz_text[i_ring_cel])

  gro_out_content.append(' 999.9 999.9 999.9' + "\n")

  # write gro output 
  gro_out = open(gro_out_name,"w+")
  gro_out.writelines(gro_out_content)
  gro_out.close()
  print("    ... saved as :",gro_out_name)
  print("")            

#-------------------------------------------------------------------------------------------
# This function generates coordinates of all the beads constituting 
# the ESI capillary, the image capillary, and the counter electrode
#
# Format is as follows:
#   ACTUAL CAPILLARY beadS
#   i_bead=0                         j_coord=0: x  j_coord=1: y  j_coord=2: z  j_coord=3: q
#   i_bead=1                         j_coord=0: x  j_coord=1: y  j_coord=2: z  j_coord=3: q 
#   ...                                                                          
#   i_bead=n_bead_cap-1              j_coord=0: x  j_coord=1: y  j_coord=2: z  j_coord=3: q    
#  ---------------------------------------------------------------------------------
#  IMAGE CAPILLARY beadS
#  i_bead=n_bead_cap                 j_coord=0: x  j_coord=1: y  j_coord=2: z  j_coord=3: q
#  i_bead=n_bead_cap+1               j_coord=0: x  j_coord=1: y  j_coord=2: z  j_coord=3: q
#  ...                                                                        
#  i_bead=2*n_bead_cap-1             j_coord=0: x  j_coord=1: y  j_coord=2: z  j_coord=3: q   
#  ---------------------------------------------------------------------------------
#  COUNTER ELECTRODE beadS
#  i_bead=2*n_bead_cap               j_coord=0: x  j_coord=1: y  j_coord=2: z  j_coord=3: q 
#  i_bead=2*n_bead_cap+1             j_coord=0: x  j_coord=1: y  j_coord=2: z  j_coord=3: q   
#   ...                                                                          
#  i_bead=2*n_bead_cap+n_bead_cel-1  j_coord=0: x  j_coord=1: y  j_coord=2: z   _coord=3: q        
#                                                                   
# bead numbering (i_bead) is as follows (schematically for n_bead_cap = 9):            
#
#   capillary   6   3   0    |    9   12  15    image
#               7   4   1    |    10  13  16
#               8   5   2    |    11  14  17
# 
#                         symmetry
#                          plane   = location of counter electrode
#
#   counter electrode numbering starts in the center and continues outward
#                        23  24
#                      22  18  19              
#                        21  20 
#
# Counter electrode beads are arranged in rings, numbered using i_ring_cel = 0.. (n_ring_cel-1)
# Most rings (defined by a unique distance from midpoint = 0/0/0) have 6 beads. Some have 12, or 18, or 24
#
# In addition: capillary probe points (cap_probe_out_xyz/ap_probe_ins_xyz) are generated, where the potential of
# the capillary is measured. These points are placed just outside/inside of bead "0" in each cap_ring
#
# In addition: counter electrode probe points (cel_probe_beh_xyz/cel_probe_fro_xyz) are generated, where the potential of
# the cel is measured. These points are placed just behind (or in front of)bead "0" in each cel_ring. 
# This allows calculating the potential AT the the cel bead, by taking the average of cel_probe_beh and cel_probe_fro values. 
#

def generate_xyz_coordinates():


  global bead_xyzq                         # capillary/image capillary/counter elctr bead coord and charges 
  global cap_probe_out_xyz                 # capillary probe coordinates on outside
  global cap_probe_ins_xyz                 # capillary probe coordinates on  inside
  global cap_ring_radius                   # radial distance of capillary rings from x-axis
  global cel_probe_beh_xyz                 # counter electrode probe coordinates
  global cel_probe_fro_xyz                 # counter electrode probe coordinates
  global cel_ring_radius                   # radial distance of counter eletrode rings from x-axis
  global n_bead_cel                        # number of counter electrode beads
  global n_bead_bwa                        # number of back wall beads
  global cel_bead_x,cel_bead_y,cel_bead_z  # counter electrode coordinate beads
  global n_ring_cel                        # number of rings in counter electrode
  global n_bead_in_ring_cel                # number of beads in each counter electrode ring
  global bwa_xyz                           # back wall coordinates
  global cap_ring_xcoord                   # x position of each ring


      
  # initiating array  bead_xyzq  for bead coordinates
  n_bead_cel_max    = 5000                                      # max number of counter electr beads
  n_bead_bwa_max    = 5000                                      # max number of back wall beads
  bead_xyzq         = np.zeros((2*n_bead_cap+n_bead_cel_max,4)) # capillary bead coordinates and charges
  cap_probe_out_xyz = np.zeros((n_ring_cap,3))                  # cap_probe_out coordinate points
  cap_probe_ins_xyz = np.zeros((n_ring_cap,3))                  # cap_probe_ins coordinate points
  cel_probe_beh_xyz = np.zeros((n_ring_cel_max,3))              # cel_probe_beh coordinate points 
  cel_probe_fro_xyz = np.zeros((n_ring_cel_max,3))              # cel_probe_fro coordinate points 
  bwa_xyz           = np.zeros((n_bead_bwa_max,3))              # back wall bead coordinates 

  i_bead           = 0                                         # counts beads in capillary/image /counter electr
  i_bead_cap_probe = 0                                         # counts cap_probe beads (for both _out and _ins)

  delta_x          = 0.866 * 2. * bead_radius_cap

  # regardless whether there is a counter electrode, the position of the capillary tip 
  #(cap_ring[0]) is determined by the variable x_cel_dist: x = -x_cel_dist
  x_coord_start = -1. * x_cel_dist
#  # with    counter electrode, the capillary tip (cap_ring[0]) is at x = -x_cel_dist
#  # without counter electrode, the capillary tip (cap_ring[0]) is at x = 0
#  if(cel == "yes"): x_coord_start = -1. * x_cel_dist  
#  if(cel == "no"):  x_coord_start =  0.

  x_coord   = x_coord_start
  phistart  = 0.
  cap_ring_xcoord = []
  cap_ring_radius = []


  # Generating capillary coordinates (actual beads, placed at negative x values)

  for i_ring_cap in range(0,n_ring_cap): # capillary beads are placed in stacked cap_rings
    # current capillary radius is determined by the number of beads in current cap_ring
    cap_radius      = bead_radius_cap / ( np.sin(2*np.pi/2./n_bead_in_ring_cap[i_ring_cap]) )
    cap_ring_radius.append(cap_radius)

    for i_bead_in_ring_cap in range(0,n_bead_in_ring_cap[i_ring_cap]):

      delta_phi     = 2 * np.pi / n_bead_in_ring_cap[i_ring_cap]  # angle increment 
      phi           = phistart + delta_phi  * i_bead_in_ring_cap  # angle for polar coordinates

      bead_xyzq[i_bead,0] = x_coord
      bead_xyzq[i_bead,1] = cap_radius * np.sin(phi)
      bead_xyzq[i_bead,2] = cap_radius * np.cos(phi)

      # this allows the use of ""rings"" with single beads
      if (n_bead_in_ring_cap[i_ring_cap] == 1):
        bead_xyzq[i_bead,1] = 0.
        bead_xyzq[i_bead,2] = 0.

      i_bead = i_bead + 1

      if (i_bead_in_ring_cap == 0):
 
       # cap_probe beads are placed outside and inside of "0" bead in each cap_ring
        cap_probe_out_xyz[i_bead_cap_probe,0] = x_coord
        cap_probe_out_xyz[i_bead_cap_probe,1] = (cap_radius +  0.5 * bead_radius_cap) * np.sin(phi)
        cap_probe_out_xyz[i_bead_cap_probe,2] = (cap_radius +  0.5 * bead_radius_cap) * np.cos(phi)
        cap_probe_ins_xyz[i_bead_cap_probe,0] = x_coord
        cap_probe_ins_xyz[i_bead_cap_probe,1] = (cap_radius -  0.5 * bead_radius_cap) * np.sin(phi)
        cap_probe_ins_xyz[i_bead_cap_probe,2] = (cap_radius -  0.5 * bead_radius_cap) * np.cos(phi)

        # this allows the use of "rings" with single beads
        if (n_bead_in_ring_cap[i_ring_cap] == 1):
          cap_probe_out_xyz[i_bead_cap_probe,1] = 0.
          cap_probe_out_xyz[i_bead_cap_probe,2] = 0.5 * bead_radius_cap
          cap_probe_ins_xyz[i_bead_cap_probe,1] = 0.
          cap_probe_ins_xyz[i_bead_cap_probe,2] =-0.5 * bead_radius_cap

        i_bead_cap_probe = i_bead_cap_probe + 1

    cap_ring_xcoord.append(x_coord)
    x_coord  = x_coord - delta_x
    phistart = phistart + delta_phi/2


  # Generating capillary coordinates for image beads (positive x values)

  i_bead     = n_bead_cap
  for i_ring_cap in range(0,n_ring_cap): # capillary beads are placed in stacked cap_rings
    for i_bead_in_ring_cap in range(0,n_bead_in_ring_cap[i_ring_cap]):

      image_bead_x = abs(bead_xyzq[i_bead-n_bead_cap,0])

      # Image bead coordinates are always present. If the program runs without counter electrode,
      # these image beads are shifted  to the right, so that they don't clash with actual beads
      # during distance calculations.  
      # [note that, for cel == "no", all image charges are set to zero]
      if(cel == "no"): image_bead_x =  image_bead_x + 10.
      bead_xyzq[i_bead,0] = image_bead_x
      bead_xyzq[i_bead,1] = bead_xyzq[i_bead-n_bead_cap,1]
      bead_xyzq[i_bead,2] = bead_xyzq[i_bead-n_bead_cap,2]

      i_bead = i_bead + 1


  n_bead_cel = 0

  # Generating counter electrode beads as a flat disk at x = 0.
  # These beads are arranged in n_ring_cel concentric rings
  # Also: generating counter electrode probe sites (cel_probe_beh and cel_probe_fro)
  # Counter electrode coordinates are always generated (even for cel = "no").


  # First generating a rectangular grid of data points in the yz plane
  # with hexagonal packing. Rings of counter electrode beads (cel_ring) will then be selected from
  # these grid points.


  yz_grid_coord_cel = np.zeros(( (4*n_cel_grid_size+1)**2, 2)) 

  i_grid_point = 0

  for i_z in range(-n_cel_grid_size,n_cel_grid_size+1):
    for i_y in range(-3*n_cel_grid_size,3*n_cel_grid_size+1):
      yz_grid_coord_cel[i_grid_point,0] =  i_y * bead_radius_cel
      if ( (i_y % 2) == 0): factor = 0.   #i_y is even
      if ( (i_y % 2) != 0): factor = 0.5  #i_y is odd
      yz_grid_coord_cel[i_grid_point,1] =  (i_z + factor) * bead_radius_cel * (np.sqrt(12)) 
      i_grid_point =  i_grid_point + 1

  n_grid_point = i_grid_point - 1

  # The grid midpoint (located at x = 0, y = 0, z = 0) has i_grid_point = int(n_grid_point/2)
  # The following is just an internal check, probably not necessary
  grid_center_y = yz_grid_coord_cel[int(n_grid_point/2),0]
  grid_center_z = yz_grid_coord_cel[int(n_grid_point/2),1]
  if ( (grid_center_y != 0.) or (grid_center_z != 0.) ):
    print("  ERROR: counter electrode is off-center")
    sys.exit()


  # Identifying concentric cel_rings in the grid of points in the x = 0 plane

  cel_radial_dist = []
  for i_grid_point in range(0,n_grid_point):
    cel_radial_dist.append( np.sqrt( yz_grid_coord_cel[i_grid_point,0]**2 
                                        + yz_grid_coord_cel[i_grid_point,1]**2) )
      
  cel_radial_dist_sorted = sorted(cel_radial_dist)
  

  # Identifying unique distance values (generating a list where each distance
  # is only listed once)
    
  cel_ring_radius = []
  concentric_distance_test    = 1E99

  for i_grid_point in range(0,n_grid_point):
    if (concentric_distance_test != 0.):
      dist_ratio = cel_radial_dist_sorted[i_grid_point] / concentric_distance_test
    if ( (dist_ratio < 0.999999999) or (dist_ratio > 1.000000001) ):  #allowing for some wiggle room bc of numerical issues
      cel_ring_radius.append(cel_radial_dist_sorted[i_grid_point])
      concentric_distance_test         = cel_radial_dist_sorted[i_grid_point]

  # Now that the unique distances are known, cel_rings can be identified
  n_bead_in_ring_cel = []

    
  # the first cel_ring has just one bead (located at y=0, z=0)
  n_bead_in_ring_cel.append(1)

  bead_xyzq[2*n_bead_cap,0] = 0.
  bead_xyzq[2*n_bead_cap,1] = 0.
  bead_xyzq[2*n_bead_cap,2] = 0.
  
  #cel_probe_beh sites are located at positive x, just behind cel beads
  cel_probe_beh_xyz[0,0] = 0.5 * bead_radius_cap
  cel_probe_beh_xyz[0,1] = 0.
  cel_probe_beh_xyz[0,2] = 0.
  #cel_probe_fro sites are located at negative x, just behind cel beads
  cel_probe_fro_xyz[0,0] =-0.5 * bead_radius_cap
  cel_probe_fro_xyz[0,1] = 0.
  cel_probe_fro_xyz[0,2] = 0.

  n_bead_cel = 1
  i_bead     = 2*n_bead_cap
  n_ring_cel = 1

  for i_ring_cel in range(1,len(cel_ring_radius)):
    i_bead_in_ring_cel = 0
    for i_grid_point in range(0,n_grid_point):
      dist_ratio = cel_radial_dist[i_grid_point] / cel_ring_radius[i_ring_cel]
      if ( (dist_ratio > 0.999999) and (dist_ratio < 1.000001) ):  #allowing for some wiggle room bc of numerical issues
        # only transfering data points that form complete circles
        if(cel_radial_dist[i_grid_point] < 2*n_cel_grid_size*bead_radius_cel): 
          i_bead              = i_bead + 1
          i_bead_in_ring_cel  = i_bead_in_ring_cel + 1
          
          # assigning counter electrode coordinates
          bead_xyzq[i_bead,0] = 0.
          bead_xyzq[i_bead,1] = yz_grid_coord_cel[i_grid_point,0]
          bead_xyzq[i_bead,2] = yz_grid_coord_cel[i_grid_point,1]

          n_ring_cel          = i_ring_cel+1

          if(n_ring_cel > n_ring_cel_max-1):
            print("")
            print("  ERROR: n_ring_cel too large")
            sys.exit()

    n_bead_in_ring_cel.append(i_bead_in_ring_cel)
    n_bead_cel = n_bead_cel + i_bead_in_ring_cel


  i_bead = 2*n_bead_cap
  for i_ring_cel in range(1,n_ring_cel):    
    # assigning counter electrode probe coordinates
    cel_probe_beh_xyz[i_ring_cel,0] = 0.5 * bead_radius_cap
    cel_probe_beh_xyz[i_ring_cel,1] = bead_xyzq[i_bead+1,1]
    cel_probe_beh_xyz[i_ring_cel,2] = bead_xyzq[i_bead+1,2]
    cel_probe_fro_xyz[i_ring_cel,0] =-0.5 * bead_radius_cap
    cel_probe_fro_xyz[i_ring_cel,1] = bead_xyzq[i_bead+1,1]
    cel_probe_fro_xyz[i_ring_cel,2] = bead_xyzq[i_bead+1,2]
    i_bead = i_bead + n_bead_in_ring_cel[i_ring_cel]

  #discarding outermost cel ring, because it might be incomplete
  n_bead_cel = n_bead_cel - n_bead_in_ring_cel[n_ring_cel-1]
  n_ring_cel = n_ring_cel-1   

  if(cel == "yes"):
    print("  number of rings in counter electr: ",n_ring_cel)
    print("")



  # Generating back wall (bwa) coordinates from yz sheet of hexagonally packed beads

  yz_grid_coord_bwa = np.zeros((n_bead_bwa_max, 2))
  x_bwa                  = bead_xyzq[0,0] - abs(bwa_dist)
  
  # Finding x_bwa_true (which will be slightly different from x_bwa), to ensure that 
  # the back wall is exactly aligned with the x-position of one of the rings
  x_diff_min  = 1E20
  for i_bead_cap in range(0,n_bead_cap):
    x_diff = abs(x_bwa - bead_xyzq[i_bead_cap,0])
    if (x_diff < x_diff_min):
      x_diff_min             = x_diff
      x_bwa_true             = bead_xyzq[i_bead_cap,0]
      #this is the capillary radius at the location of the back wall
      cap_radius_at_bwa_xpos = np.sqrt(bead_xyzq[i_bead_cap,1]**2 + bead_xyzq[i_bead_cap,2]**2)


  i_grid_point = 0
  i_bead_bwa   = 0
  for i_z in range(-20,20+1):
    for i_y in range(-3*20,3*20+1):
      yz_grid_coord_bwa[i_grid_point,0] =  i_y * bead_radius_bwa
      if ( (i_y % 2) == 0): factor = 0.   #i_y is even
      if ( (i_y % 2) != 0): factor = 0.5  #i_y is odd
      yz_grid_coord_bwa[i_grid_point,1] =  (i_z + factor) * bead_radius_bwa * (np.sqrt(12))
      # selecting grid points that are inside the capillary ...
      if( np.sqrt(yz_grid_coord_bwa[i_grid_point,0]**2 + yz_grid_coord_bwa[i_grid_point,1]**2) 
          <= cap_radius_at_bwa_xpos ):
        # ... and that do not clash with any capillary beads
        dist_min = 1E20
        for i_bead_cap in range(0,n_bead_cap):
          dist = np.sqrt(  (x_bwa_true                        - bead_xyzq[i_bead_cap,0])**2
                         + (yz_grid_coord_bwa[i_grid_point,0] - bead_xyzq[i_bead_cap,1])**2
                         + (yz_grid_coord_bwa[i_grid_point,1] - bead_xyzq[i_bead_cap,2])**2 )
          if (dist < dist_min):
            dist_min = dist
        if ( dist_min > (bead_radius_bwa + bead_radius_bwa) ):
          bwa_xyz[i_bead_bwa,0] = x_bwa_true
          bwa_xyz[i_bead_bwa,1] = yz_grid_coord_bwa[i_grid_point,0]
          bwa_xyz[i_bead_bwa,2] = yz_grid_coord_bwa[i_grid_point,1]
          i_bead_bwa            = i_bead_bwa + 1

      i_grid_point = i_grid_point + 1

  n_bead_bwa   = i_bead_bwa


#------------------------------------------------------------------------------------------
# This function calculates distances between capillary surface probes and all beads, to save
# time for future calculations (these distances will be needed many times)
#
# cap_probe_out_dist contains r_ij distances (nm) for each cap_probe_out bead (i_ring_cap) 
# with all capillary beads AND all image beads AND all counter electrode beads (j_bead)
#     j_bead-->  0   1   2   3   4   5  ...
#  i_ring_cap  0    00  01  02  03  04  05 ...
#  i_ring_cap  1    10  11  12  13  14  15 ...
#  i_ring_cap  2    20  21  22  23  24  25 ...
#  i_ring_cap  3    30  31  32  33  34  35
#  ...
# The same steps are performed for each cap_probe_ins bead
#
def calculate_cap_probe_dist():

  global cap_probe_out_dist
  global cap_probe_ins_dist

  print("  calculating distances for capillary probes ...")
  
  #initiating matrix with zeroes
  cap_probe_out_dist = np.zeros((n_ring_cap,2*n_bead_cap + n_bead_cel))
  cap_probe_ins_dist = np.zeros((n_ring_cap,2*n_bead_cap + n_bead_cel))


  for i_ring_cap in range(0,n_ring_cap):
    for j_bead in range(0,2*n_bead_cap + n_bead_cel):

      cap_probe_out_dist[i_ring_cap,j_bead] = np.sqrt(   (cap_probe_out_xyz[i_ring_cap,0] - bead_xyzq[j_bead,0])**2
                                                       + (cap_probe_out_xyz[i_ring_cap,1] - bead_xyzq[j_bead,1])**2
                                                       + (cap_probe_out_xyz[i_ring_cap,2] - bead_xyzq[j_bead,2])**2  )

      cap_probe_ins_dist[i_ring_cap,j_bead] = np.sqrt(   (cap_probe_ins_xyz[i_ring_cap,0] - bead_xyzq[j_bead,0])**2
                                                       + (cap_probe_ins_xyz[i_ring_cap,1] - bead_xyzq[j_bead,1])**2
                                                       + (cap_probe_ins_xyz[i_ring_cap,2] - bead_xyzq[j_bead,2])**2  )
  print("  ... done")
  print("  ")

#------------------------------------------------------------------------------------------
# This function calculates distances between counter electrode probes and all beads, to save
# time for future calculations (these distances will be needed many times)
#
# cel_probe_beh_dist contains r_ij distances (nm) for each cel_probe_beh bead (i_ring_cel) 
# with all capillary beads AND all image beads AND all counter electrode beads (j_bead)
#     j_bead-->      0   1   2   3   4   5  ...
#  i_ring_cel  0    00  01  02  03  04  05 ...
#  i_ring_cel  1    10  11  12  13  14  15 ...
#  i_ring_cel  2    20  21  22  23  24  25 ...
#  i_ring_cel  3    30  31  32  33  34  35
#  ...
#
#
def calculate_cel_probe_dist():

  global cel_probe_beh_dist
  global cel_probe_fro_dist

  print("  calculating distances for counter electrode probes ...")
  
  #initiating matrix with zeroes
  cel_probe_beh_dist = np.zeros((n_ring_cel,2*n_bead_cap + n_bead_cel))
  cel_probe_fro_dist = np.zeros((n_ring_cel,2*n_bead_cap + n_bead_cel))


  for i_ring_cel in range(0,n_ring_cel):
    for j_bead in range(0,2*n_bead_cap + n_bead_cel):

      cel_probe_beh_dist[i_ring_cel,j_bead] = np.sqrt(   (cel_probe_beh_xyz[i_ring_cel,0] - bead_xyzq[j_bead,0])**2
                                                       + (cel_probe_beh_xyz[i_ring_cel,1] - bead_xyzq[j_bead,1])**2
                                                       + (cel_probe_beh_xyz[i_ring_cel,2] - bead_xyzq[j_bead,2])**2  )

      if (cel_probe_beh_dist[i_ring_cel,j_bead] < bead_radius_cel):
        cel_probe_beh_dist[i_ring_cel,j_bead] = bead_radius_cel  # for smoothing the scan 

      cel_probe_fro_dist[i_ring_cel,j_bead] = np.sqrt(   (cel_probe_fro_xyz[i_ring_cel,0] - bead_xyzq[j_bead,0])**2
                                                       + (cel_probe_fro_xyz[i_ring_cel,1] - bead_xyzq[j_bead,1])**2
                                                       + (cel_probe_fro_xyz[i_ring_cel,2] - bead_xyzq[j_bead,2])**2  )

      if (cel_probe_fro_dist[i_ring_cel,j_bead] < bead_radius_cel):
        cel_probe_fro_dist[i_ring_cel,j_bead] = bead_radius_cel  # for smoothing the scan 

  print("  ...done")
  print("")

#------------------------------------------------------------------------------------------
# This function calculates the electric potential at the site of each cap_probe_out:
#     cap_probe_pot[i_ring_cap] = SUM(qj/dij)    where dij are distances.
# This sum extends over all capillary beads AND image beads AND counter electrode beads.
# In this version, each ap_probe_pot[i_ring_cap] value is the AVERAGE of the potential measured by
# cap_probe_out_pot and cap_probe_ins_pot
#
# The following matrix elements are relevant:
#    j_bead-->  0       1       2     ...
#  i_ring_cap 0    q0/d00  q1/d01  q2/d02 ...    SUM = cap_probe_out_pot[0]
#  i_ring_cap 1    q1/d10  q1/d11  q2/d12 ...    SUM = cap_probe_out_pot[1]
#  i_ring_cap 2    q2/d20  q1/d21  q2/d22 ...    SUM = cap_probe_out_pot[2]
#  ...
#
#
def calculate_cap_probe_pot():

  global cap_probe_pot

  cap_probe_out_pot =  np.zeros(n_ring_cap) # initiating cap_probe_out potentials
  cap_probe_ins_pot =  np.zeros(n_ring_cap) # initiating cap_probe_ins potentials
  cap_probe_pot     =  np.zeros(n_ring_cap) # initiating cap_probe average potentials

  for i_ring_cap in range(0,n_ring_cap):
    for j_bead in range(0,2*n_bead_cap + n_bead_cel):

      if(bead_xyzq[j_bead,3] != 0.):  #this if statement speeds things up a bit
 
        q_over_r          = bead_xyzq[j_bead,3] / cap_probe_out_dist[i_ring_cap,j_bead]
        cap_probe_out_pot[i_ring_cap] = cap_probe_out_pot[i_ring_cap] + q_over_r

        q_over_r          = bead_xyzq[j_bead,3] / cap_probe_ins_dist[i_ring_cap,j_bead]
        cap_probe_ins_pot[i_ring_cap] = cap_probe_ins_pot[i_ring_cap] + q_over_r 

  # Now calculating average potentials
  for i_ring_cap in range(0,n_ring_cap):
    cap_probe_pot[i_ring_cap] = 0.5*( cap_probe_out_pot[i_ring_cap] + cap_probe_ins_pot[i_ring_cap] )
    #cap_probe_pot[i_ring_cap] = cap_probe_out_pot[i_ring_cap]

#------------------------------------------------------------------------------------------
# This function calculates the electric potential at the site of each cel_probe_beh:
#     cel_probe_beh_pot[i_ring_cel] = SUM(qj/dij)    where dij are distances.
# This sum extends over all capillary beads AND image beads AND counter electrode beads
#
# The following matrix elements are relevant:
#    j_bead-->  0       1       2     ...
#  i_ring_cel 0    q0/d00  q1/d01  q2/d02 ...    SUM = cel_probe_beh_pot[0]
#  i_ring_cel 1    q1/d10  q1/d11  q2/d12 ...    SUM = cel_probe_beh_pot[1]
#  i_ring_cel 2    q2/d20  q1/d21  q2/d22 ...    SUM = cel_probe_beh_pot[2]
#  ...
#
# Then the same procedure is performed for the cel_probe_fro sites.
# Then the potential of each cel site is calculated as the average of the pos and neg values
# (keeping in mind that the pos/neg pairs sandwich the cel points of interest)
def calculate_cel_probe_pot():

  global cel_probe_pot

  cel_probe_pot     = np.zeros(n_ring_cel) # initiating cel_probe     potentials

  cel_probe_beh_pot = np.zeros(n_ring_cel) # initiating cel_probe_beh potentials
  for i_ring_cel in range(0,n_ring_cel):
    for j_bead in range(0,2*n_bead_cap + n_bead_cel):
      if(bead_xyzq[j_bead,3] != 0.):  #this if statement speeds things up a bit
        q_over_r                  = bead_xyzq[j_bead,3] / cel_probe_beh_dist[i_ring_cel,j_bead]
        cel_probe_beh_pot[i_ring_cel] = cel_probe_beh_pot[i_ring_cel] + q_over_r

  cel_probe_fro_pot = np.zeros(n_ring_cel) # initiating cel_probe_fro potentials
  for i_ring_cel in range(0,n_ring_cel):
    for j_bead in range(0,2*n_bead_cap + n_bead_cel):
      if(bead_xyzq[j_bead,3] != 0.):  #this if statement speeds things up a bit
        q_over_r                  = bead_xyzq[j_bead,3] / cel_probe_fro_dist[i_ring_cel,j_bead]
        cel_probe_fro_pot[i_ring_cel] = cel_probe_fro_pot[i_ring_cel] + q_over_r

  #Calculating average value for each cel probe point
  for i_ring_cel in range(0,n_ring_cel):
    cel_probe_pot[i_ring_cel] = 0.5*(cel_probe_beh_pot[i_ring_cel] + cel_probe_fro_pot[i_ring_cel])
    #cel_probe_pot[i_ring_cel] = cel_probe_beh_pot[i_ring_cel]
#------------------------------------------------------------------------------------------
# This function calculates the total potential energy as SUM( Vij = q_i * q_j / r_ij )
# All terms above the i=j diagonal are included.
#
#       j=0  j=1  j=2  j=3  j=4  j=5
#  i=0  x    V01  V02  V03  V04  V05
#  i=1  x    x    V12  V13  V14  V15 
#  i=2  x    x    x    V23  V24  V25
#  i=3  x    x    x    x    V34  V35    
#  i=4  x    x    x    x    x    V45
#  i=5  x    x    x    x    x    x
#
# This includes capillary beads, image beads, and counter electrode beads,
# subject to the "if" stattements below
# 
# Note that there are either image charges OR counter electrode charges, but never both.
# We use the same function regardless, because 
# if there are image             charges all counter electrode charges are zero
# if there are counter electrode charges all             image charges are zero

def calculate_E_pot():

  result = 0.

  for i_bead in range(0,2*n_bead_cap + n_bead_cel):
    for j_bead in range(i_bead+1,2*n_bead_cap + n_bead_cel):

      r_ij = np.sqrt(  (bead_xyzq[i_bead,0] - bead_xyzq[j_bead,0])**2
                     + (bead_xyzq[i_bead,1] - bead_xyzq[j_bead,1])**2
                     + (bead_xyzq[i_bead,2] - bead_xyzq[j_bead,2])**2  )

      V_ij = bead_xyzq[i_bead,3] * bead_xyzq[j_bead,3] / r_ij

      # When calculating the energy of the system, interactions among image charges or counter electrode
      # charges do not count.
      # Reason: Bringing two image charges from inifinity into a certain distance with one another does  
      # not cost energy, as they move in a field-free region. See Griffiths/Electrodynamics
      if( i_bead > (n_bead_cap-1)) and (j_bead > (n_bead_cap-1) ): V_ij = 0.


      # The interaction of an actual charge on the capillary with an image charge is only half of what 
      # it would be between two actual charges (same reason as above, see Griffiths/Electrodynamics)
      # The interaction of an actual charge on the capillary with an induced charge on the counter electrode
      # is also divided by a factor of two.

      # selecting capillary beads    
      if( i_bead < n_bead_cap):
        # selecting image beads or induced charges 
        if( j_bead > (n_bead_cap-1) ):
          V_ij = V_ij / 2

      result = result + V_ij

#      # selecting capillary beads    
#      if( i_bead < n_bead_cap):
#        # selecting image beads
#        if( (j_bead > (n_bead_cap-1)) and (j_bead < 2*n_bead_cap) ):
#          V_ij = V_ij / 2
#
#      result = result + V_ij

  return result


#------------------------------------------------------------------------------------------
# start of main program

global n_bead_cap
global n_ring_cap
global fit_summary_cap
global fit_summary_cel
global n_ring_cel_max                    # maximum nunber of rings in counter electrode
global charge_recov_on_cel



print("                              ")
print("  ############################")
print("  #                          #")
print("  #  capillary_potential.py  #")
print("  #                          #")
print("  ############################")
print("                              ")


n_command_line = len(sys.argv)

if(n_command_line != 2):
  print("  ")
  print("  ### ERROR: Command_line not ok.")
  print("             Expected 2 arguments, found ",n_command_line)
  print("  ")
  print(sys.argv)
  sys.exit()

input_file_name        = sys.argv[1]   

print("  input file       : ",input_file_name)
print("")


path.exists(input_file_name)
if (os.path.isfile(input_file_name) == False):
  print("  ")
  print("  ### ERROR: file not found ",input_file_name)
  print("  ")
  sys.exit()

#reading input parameters
read_input(input_file_name)

x_cel_dist = abs(x_cel_dist) # Making sure that this is a positive number,
                             # so that capillary tip is at negative x-value


n_ring_cel_max    = 300                   # max number of counter electr rings

n_ring_cap      = len(n_bead_in_ring_cap) # number of cap_rings (excluding image beads)

n_ring_cap_max  = 350                     # raise this number if necessary 
if (n_ring_cap >= n_ring_cap_max-2):
  print("  ")
  print("  ERROR: n_ring_cap greater than n_ring_cap_max")
  sys.exit()

# Initiate array fit_summary_cap (contains E_pot,
#                              cap_probe_pot[i_ring_cap],
#                              charge_per_bead_in_ring_cap[i_ring_cap])
#                              for fitting of capillary charge distribution
fit_summary_cap = np.zeros((n_ring_cap_max+1,2*n_capfit_iter+2))

# Initiate array fit_summary_cel. Same as above, but for
#                              for fitting of counter electrode charge distribution
fit_summary_cel = np.zeros((n_ring_cel_max+1,2*n_celfit_iter+2))


n_bead_cap  = sum(n_bead_in_ring_cap)     # number of beads in actual capillary (excluding image beads)
print("")
print("  number of beads in capillary     : ",n_bead_cap )
print("  number of rings in capillary     : ",n_ring_cap)

#generating capillary coordinates
generate_xyz_coordinates()

#calculate distances for capillary probes
calculate_cap_probe_dist()

#calculate distances for counter electrode probes
if (cel == "yes"):
  calculate_cel_probe_dist()   

# initiate charge_per_bead_in_ring_cap (similar information as above,
# keeping in mind that all beads in a given cap_ring have the same charge)
charge_per_bead_in_ring_cap = np.zeros(n_ring_cap)
# ... and the same for charge_per_bead_in_ring_cel
charge_per_bead_in_ring_cel = np.zeros(n_ring_cel_max)


# Initially, each capillary bead has the same charge ( =1).
# This is in dimensionless units, but may be interpreted as elementary charges (e) if length is in nm

charge_cap = 0
for i_bead in range(0,n_bead_cap):
  bead_xyzq[i_bead,3] = 1. 
  charge_cap          = charge_cap + bead_xyzq[i_bead,3]       # total charge of capillary (in dimensionless units or e)
                                                               # excluding image charges
for i_ring_cap in range(0,n_ring_cap):
  charge_per_bead_in_ring_cap[i_ring_cap] = 1. 

if(cel == "yes"):                                              # image charges are assigned 
  for i_bead in range(n_bead_cap,2*n_bead_cap):                # for i_bead n_bead_cap to (2*n_bead_cap-1).      
    bead_xyzq[i_bead,3] = -1. * bead_xyzq[i_bead-n_bead_cap,3] # If there is no counter electrode, all image charges stay at zero
                                                               # (but neutral image beads are still present)


# calculate initial potential energy of entire capillary
E_pot = calculate_E_pot()

# calculate initial potential at the surface of each cap_ring
calculate_cap_probe_pot()

# keep track of initial E_tot, cap_probe_pot, charge_per_bead_in_ring_cap
fit_summary_cap[0,0] = -1.    # current round of fitting (-1 means unmodified input parameters) 
fit_summary_cap[0,1] = E_pot  # current value of E_pot
for i_ring_cap in range(0,n_ring_cap):
  fit_summary_cap[i_ring_cap+1,0] = cap_probe_pot              [i_ring_cap]
  fit_summary_cap[i_ring_cap+1,1] = charge_per_bead_in_ring_cap[i_ring_cap]
 

print("")
print("                             initial E_pot = ",E_pot)

##############################################################################
# Now fitting the charges on the capillary, guided by the principle that all #
# cap_probe_pot[i_ring_cap] values should be as similar as possible (meaning #
# that the potential on the metal capillary is constant).                    #
#                                                                            #
#       START OF FITTING LOOP TO DETERMINE CHARGES ON CAP_RINGS              #
#                   (without or with image capillary)                        #
#                                                                            #
##############################################################################

print("  ")
print("  Fitting charge distribution on capillary ...")
print("  ============================================")
print("  ")
for i_capfit_iter in range(0,n_capfit_iter):

  # calculate average potential of all cap_rings, based on charge
  # distribution from previous iteration (or initial charge distribution, for i_capfit_iter = 0)
  cap_probe_pot_av = 0.
  for i_ring_cap in range(0,n_ring_cap):
    cap_probe_pot_av = cap_probe_pot_av + cap_probe_pot[i_ring_cap]

  cap_probe_pot_av = cap_probe_pot_av/n_ring_cap
 

  # now attempting to make the potential of each cap_ring equal to the average potential

  j_bead_end      = 0
  charge_cap_temp = 0.

  for i_ring_cap in range(0,n_ring_cap):

    j_bead_start = j_bead_end
    j_bead_end   = j_bead_end   + n_bead_in_ring_cap[i_ring_cap]

    cap_probe_pot_fac = cap_probe_pot[i_ring_cap] / cap_probe_pot_av  # cap_probe_pot_fac > 1: potential on cap_ring is too high
                                                                      # cap_probe_pot_fac < 1: potential on cap_ring is too low

    if(i_ring_cap == 0           ): cap_probe_pot_fac = cap_probe_pot_fac**2 # First and last ring need a bit of extra help during
    if(i_ring_cap == n_ring_cap-1): cap_probe_pot_fac = cap_probe_pot_fac**2 # rescaling in next step. Without this modification,
                                                                             # first and last ring tend to retain slightly below-average
                                                                             # potentials.

    # each bead charge is scaled up or down; all beads within a cap_ring are treated equally
    for j_bead in range(j_bead_start,j_bead_end):
      bead_xyzq[j_bead,3] = bead_xyzq[j_bead,3] / cap_probe_pot_fac
      charge_cap_temp     = charge_cap_temp + bead_xyzq[j_bead,3]   # charge_cap_temp is the total capillary charge after rescaling
      if(j_bead == j_bead_start):
        charge_per_bead_in_ring_cap[i_ring_cap] = bead_xyzq[j_bead,3]
    
    charge_corr_fac = charge_cap_temp/charge_cap          # charge_corr_fac > 1: capillary charge is too large after rescaling
                                                          # charge_corr_fac < 1: capillary charge is too small after rescaling

  # charges are rescaled, so that the proper (initial) total capillary charge is restored
  for j_bead in range(0,n_bead_cap):
    bead_xyzq[j_bead,3] = bead_xyzq[j_bead,3] / charge_corr_fac # capillary bead charges are assigned 
  for i_ring_cap in range(0,n_ring_cap):                        # for j_bead 0 to n_bead_cap-1
    charge_per_bead_in_ring_cap[i_ring_cap] = charge_per_bead_in_ring_cap[i_ring_cap] / charge_corr_fac

  if(cel == "yes"):                                              # Image bead charges are assigned
    for j_bead in range(n_bead_cap,2*n_bead_cap):                # for j_bead n_bead_cap to (2*n_bead_cap-1)  
      bead_xyzq[j_bead,3] = -1. * bead_xyzq[j_bead-n_bead_cap,3] # If there is no counter electrode, all image charges stay at zero
                                                                 # (but neutral image beads are still present)

  # calculate new total potential energy
  E_pot = calculate_E_pot()
  # calculate new potential on each cap_ring
  calculate_cap_probe_pot()

  # keep track of current E_tot, cap_probe_pot, charge_per_bead_in_ring_cap
  fit_summary_cap[0,2*i_capfit_iter+2] = i_capfit_iter  # current round of fitting 
  fit_summary_cap[0,2*i_capfit_iter+3] = E_pot          # current value of E_pot
  for i_ring_cap in range(0,n_ring_cap):
    fit_summary_cap[i_ring_cap+1,2*i_capfit_iter+2] = cap_probe_pot              [i_ring_cap]
    fit_summary_cap[i_ring_cap+1,2*i_capfit_iter+3] = charge_per_bead_in_ring_cap[i_ring_cap]


  print("    capillay fitting iteration = ",i_capfit_iter," / ",n_capfit_iter,"   E_pot = ",E_pot)

####### END OF FITTING LOOP (TO DETERMINE CHARGES ON CAP_RINGS)

print("")
print("  RESULTS AFTER CAPILLARY FITTING (dimensionless units):")
print("  E_pot = ",E_pot)
print("")
print("n_bead_in_ring_cap  cap_probe_pot       charge_per_bead_in_ring_cap")
for i_ring_cap in range(0,n_ring_cap):
  line =str( "  " + format(n_bead_in_ring_cap[i_ring_cap],            "10") 
                  + format(cap_probe_pot[i_ring_cap],              "19.8f")  
                  + format(charge_per_bead_in_ring_cap[i_ring_cap],"19.8f") )
  print(line)


# Calculate centerline potential (V) and field (V/[a.u.] or V/nm)
# in the presence of image capillary (note that image charges are zero for cel = "no")
if (cel == "yes"): calculate_center_pot_field("_center_pot_V_field_Vnm_with_img.txt")
if (cel == "no" ): calculate_center_pot_field("_center_pot_V_field_Vnm.txt")

# perform charge scan of capillary
calculate_xy_cap_charge_map("_xy_cap_charge_e.txt")

# perform potential scan in xy plane
#print("mapping deactivated")
if (cel == "yes"): calculate_xy_pot_map("_xy_pot_V_with_img.txt")
if (cel == "no" ): calculate_xy_pot_map("_xy_pot_V.txt")

# perform electric field scans in xy plane
#print("mapping deactivated")
if (cel == "yes"): calculate_xy_elfield_map("_xy_elfield_Vnm_dist_1_with_img.txt",1)
if (cel == "no" ): calculate_xy_elfield_map("_xy_elfield_Vnm_dist_1.txt",1)


charge_recov_on_cel = 0.

# write cap fitting stats to output file
write_fit_summary_cap("_fit_summary_cap.txt")





###############################################################################
# The optimal capillary charges are now known, from                           #
# the preceding image charge-based fitting routine. Now we retain these       #
# optimized capillary charges, and we delete all image charges.               #
# Instead, the program now determines the induced charges on the grounded     #
# planar counter electrode. At the end of the procedure                       #
# 1. actual counter electrode charges will be known                           #
# 2. the potential anywhere on the counter electrode will be zero             #
#  and, for a very large counter electrode:                                   #
#  3. the overall energy of the system will be the same as with image charges #
#  4. the total counter electrode charge equals the negative capillary charge #
#                                                                             #
#       START OF FITTING LOOP (TO DETERMINE CHARGES ON CEL_RINGS)             #
#                                                                             #
###############################################################################

if(cel == "yes") and (n_celfit_iter > 0):
  print("  ")
  print("  Fitting charge distribution on counter electrode ...")
  print("  ====================================================")
  print("  ")


  for i_bead in range(n_bead_cap,2*n_bead_cap):                     
    bead_xyzq[i_bead,3] = 0.     # setting all image charges to zero


  # calculate initial potential energy after removing image charges
  E_pot = calculate_E_pot()

  print("   E_pot (without image or counter electrode charges) = ",E_pot)

  # calculate potential at the surface of each cel_ring after removing image charges
  calculate_cel_probe_pot()


  #for i_ring_cel in range(0,n_ring_cel):
  #  line =str( "    " + format(n_bead_in_ring_cel[i_ring_cel],            "10") 
  #                    + format(cel_probe_pot[i_ring_cel],              "19.8f")  
  #                    + format(charge_per_bead_in_ring_cel[i_ring_cel],"19.8f") )
  #  print(line)

  # keep track of initial E_tot, cap_probe_pot, charge_per_bead_in_ring_cel
  fit_summary_cel[0,0] = -1.    # current round of fitting (-1 means unmodified input parameters) 
  fit_summary_cel[0,1] = E_pot  # current value of E_pot
  for i_ring_cel in range(0,n_ring_cel):
    fit_summary_cel[i_ring_cel+1,0] = cel_probe_pot              [i_ring_cel]
    fit_summary_cel[i_ring_cel+1,1] = charge_per_bead_in_ring_cel[i_ring_cel]


  # Calculate centerline potential and field in proper units
  # in the absence of induced counter electrode charges and without image charges)
  calculate_center_pot_field("_center_pot_V_field_Vnm_no_img_no_cel.txt")

  # perform potential scan in xy plane 
  calculate_xy_pot_map("_xy_pot_V_no_img_no_cel.txt") 

  # perform potential scan in yz counter electrode plane 
  calculate_rad_cel_pot_map("_rad_cel_pot_V_no_img_no_cel.txt")


  for i_celfit_iter in range(0,n_celfit_iter): # Start of the actual cel fitting loop

    # calculate average potential of all cel_rings, based on charge
    # distribution from previous iteration (or initial charge distribution, for i_celfit_iter = 0)
    
    cel_probe_pot_av = 0.
    for i_ring_cel in range(0,n_ring_cel):
      cel_probe_pot_av = cel_probe_pot_av + cel_probe_pot[i_ring_cel]
    cel_probe_pot_av = cel_probe_pot_av/n_ring_cel


    # now bringing the potential of each cel_ring down to zero
    j_bead_end      = 2*n_bead_cap
    charge_cel_temp = 0.

    for i_ring_cel in range(0,n_ring_cel):
      j_bead_start = j_bead_end
      j_bead_end   = j_bead_end   + n_bead_in_ring_cel[i_ring_cel]

      # each cel charge is scaled up or down; all beads within a cel_ring are treated equally
      for j_bead in range(j_bead_start,j_bead_end):

# used this for GM 007__cel
#        bead_xyzq[j_bead,3] = ( bead_xyzq[j_bead,3]         #0.5 (for smaller cel)      0.9 (for larger cel)
#                               - cel_probe_pot[i_ring_cel] * 0.9/n_cel_grid_size * bead_radius_cap 
#                               * (1. + cel_ring_radius[i_ring_cel] / cel_ring_radius[n_ring_cel])    )
        # used this for 009_cel
        #bead_xyzq[j_bead,3] = bead_xyzq[j_bead,3] - 0.6 * cel_probe_pot[i_ring_cel] * (0.5*bead_radius_cap)


        #used this for 009a_cel
        if(i_ring_cel/n_ring_cel <= 0.9):
          bead_xyzq[j_bead,3] = bead_xyzq[j_bead,3] - 0.6 * cel_probe_pot[i_ring_cel] * (0.5*bead_radius_cap)
        if(i_ring_cel/n_ring_cel >  0.9):
          bead_xyzq[j_bead,3] = bead_xyzq[j_bead,3] - 0.7 * cel_probe_pot[i_ring_cel] * (0.5*bead_radius_cap)


        if(j_bead == j_bead_start):
          charge_per_bead_in_ring_cel[i_ring_cel] = bead_xyzq[j_bead,3]


    # calculate new total potential energy
    E_pot = calculate_E_pot()
    # calculate new potential on each cel_ring
    calculate_cel_probe_pot()

    # keep track of current E_tot, cel_probe_pot, charge_per_bead_in_ring_cel
    fit_summary_cel[0,2*i_celfit_iter+2] = i_celfit_iter  # current round of fitting 
    fit_summary_cel[0,2*i_celfit_iter+3] = E_pot          # current value of E_pot
    for i_ring_cel in range(0,n_ring_cel):
      fit_summary_cel[i_ring_cel+1,2*i_celfit_iter+2] = cel_probe_pot              [i_ring_cel]
      fit_summary_cel[i_ring_cel+1,2*i_celfit_iter+3] = charge_per_bead_in_ring_cel[i_ring_cel]


    cel_q_sum = 0.
    for j_bead in range(2*n_bead_cap,2*n_bead_cap+n_bead_cel): 
      cel_q_sum = cel_q_sum + bead_xyzq[j_bead,3]

    print("    counter electrode fitting iteration = ",i_celfit_iter," / ",n_celfit_iter,"   E_pot = ",E_pot)


####### END OF FITTING LOOP (TO DETERMINE CHARGES ON CEL_RINGS)


  print("")
  charge_recov_on_cel = abs(cel_q_sum/n_bead_cap)
  print(" charge recovery on counter electrode = ",charge_recov_on_cel)
  print("n_bead_in_ring_cel   cel_probe_pot       charge_per_bead_in_ring_cel" )
  for i_ring_cel in range(0,n_ring_cel):
    line =str( "    " + format(n_bead_in_ring_cel[i_ring_cel],            "10") 
                      + format(cel_probe_pot[i_ring_cel],              "19.8f")  
                      + format(charge_per_bead_in_ring_cel[i_ring_cel],"19.8f") )
    print(line)


  # Calculate centerline potential (V) and field (V/[length unit])
  # in the presence of image capillary (note that image charges are zero for cel = "no")
  calculate_center_pot_field("_center_pot_V_field_Vnm_with_cel.txt")

  # perform potential scan in xy plane
  calculate_xy_pot_map("_xy_pot_V_with_cel.txt")

  # write cap fitting stats to output file
  write_fit_summary_cel("_fit_summary_cel.txt")

  # perform potential scan in yz counter electrode plane 
  calculate_rad_cel_pot_map("_rad_cel_pot_V.txt")

  # scan electric field
  calculate_xy_elfield_map("_xy_elfield_Vnm_dist_1_with_cel.txt",1)

# writing overall summary including charges
write_input_output_summary()

# Generate output gro file that contains all beads,
# with charge values for capillary and counter electrode
generate_gro_output("_system_with_charges.gro")
  
print("")
print("  PROGRAM EXECUTED SUCCESSFULLY")
print("")
# End of Main Program