From 1bd7709797c7a902321b8aab24d7d1aa49f5842a Mon Sep 17 00:00:00 2001 From: Walter Wissdorf Date: Wed, 18 Mar 2026 11:53:24 +0100 Subject: [PATCH 1/4] Adds gitignore file for Python --- .gitignore | 164 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 164 insertions(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..dda4f8f --- /dev/null +++ b/.gitignore @@ -0,0 +1,164 @@ +### Python template +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# poetry +# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control +#poetry.lock + +# pdm +# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. +#pdm.lock +# pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it +# in version control. +# https://pdm.fming.dev/latest/usage/project/#working-with-version-control +.pdm.toml +.pdm-python +.pdm-build/ + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# PyCharm +# JetBrains specific template is maintained in a separate JetBrains.gitignore that can +# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore +# and can be added to the global gitignore or merged into this file. For a more nuclear +# option (not recommended) you can uncomment the following to ignore the entire idea folder. +#.idea/ + From ebe460499aeb630fe8da399018296efc2133d69a Mon Sep 17 00:00:00 2001 From: Walter Wissdorf Date: Wed, 18 Mar 2026 16:18:48 +0100 Subject: [PATCH 2/4] Updates gitignore file --- .gitignore | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index dda4f8f..b20b59b 100644 --- a/.gitignore +++ b/.gitignore @@ -155,10 +155,14 @@ dmypy.json # Cython debug symbols cython_debug/ +# Results folder +results + # PyCharm # JetBrains specific template is maintained in a separate JetBrains.gitignore that can # be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore # and can be added to the global gitignore or merged into this file. For a more nuclear # option (not recommended) you can uncomment the following to ignore the entire idea folder. -#.idea/ +.idea/ + From 02253191f4dca8f95c28acaee8d7fa0f6ca39e8c Mon Sep 17 00:00:00 2001 From: Walter Wissdorf Date: Wed, 18 Mar 2026 16:20:18 +0100 Subject: [PATCH 3/4] Adds a faster version of the capillary potential script Fast version uses numba jit compilation for accerleration --- capillary_potential_fast.py | 2134 +++++++++++++++++++++++++++++++++++ 1 file changed, 2134 insertions(+) create mode 100644 capillary_potential_fast.py diff --git a/capillary_potential_fast.py b/capillary_potential_fast.py new file mode 100644 index 0000000..c061097 --- /dev/null +++ b/capillary_potential_fast.py @@ -0,0 +1,2134 @@ +#!/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 " +# +# +# 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 +from numba import njit, jit + + + +#------------------------------------------------------------------------------------------ +# 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 +@njit +def calculate_xy_pot_values(n_x, n_y, x_start, y_start, delta_x, delta_y, bead_radius_cap): + xy_pot_values = [] + + 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 + # str( format(x,"12.5f") + format(y,"12.5f") + format(pot,"12.5f") + "\n" ) + xy_pot_values.append( (x,y,pot) ) + + + return(xy_pot_values) + + +#@jit(forceobj=True) +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 + + # 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) ) + + xy_pot_values = calculate_xy_pot_values(n_x, n_y, x_start, y_start, delta_x, delta_y, bead_radius_cap) + xy_pot_content = [f'{v[0]:12.5f} {v[1]:12.5f} {v[2]:12.5f} \n' for v in xy_pot_values] # str( format(x,"12.5f") + format(y,"12.5f") + format(pot,"12.5f") + "\n" )] + + 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). +#@njit +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 + +@njit +def calculate_single_position(i_x, y_pos, n_bead_cap, n_bead_cel, center_pot_field, bead_xyzq, bead_radius_cap): + 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 + + (y_pos-(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 + + +#@jit(forceobj=True) +def calculate_axial_pot_field(rest_of_filename, y_pos=0.0): + + 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 + + calculate_single_position(i_x, y_pos, n_bead_cap, n_bead_cel, center_pot_field, bead_xyzq, bead_radius_cap) + # cutted + + + # 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 + +@njit +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_axial_pot_field("_center_pot_V_field_Vnm_with_img.txt") +if (cel == "no" ): calculate_axial_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_axial_pot_field("_center_pot_V_field_Vnm_no_img_no_cel.txt") + # Calculation of potential in axial direction off axis is possible: + #calculate_axial_pot_field("_center_pot_V_field_Vnm_no_img_no_cel_y3.txt",y_pos=3.0) + + # 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_axial_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 From fa6e0f8be9558e9fd319039656b5b0082b1be2a4 Mon Sep 17 00:00:00 2001 From: Walter Wissdorf Date: Wed, 18 Mar 2026 16:23:11 +0100 Subject: [PATCH 4/4] Adds a simple visualization script for capillary potential data --- visualize_results.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 visualize_results.py diff --git a/visualize_results.py b/visualize_results.py new file mode 100644 index 0000000..b094133 --- /dev/null +++ b/visualize_results.py @@ -0,0 +1,24 @@ +import numpy as np +import matplotlib.pyplot as plt +import sys + +def read_axial_potential(filename): + data = np.genfromtxt(filename, skip_header=1) + return data + +def plot_axial_potential(filename): + data = read_axial_potential(filename) + fig, ax = plt.subplots(2,1,figsize=(10,10)) + + ax[0].plot(data[:,0], data[:,1]) + ax[1].plot(data[:, 0], data[:, 2]) + plt.show() + print(data) + +if __name__ == '__main__': + if len(sys.argv) < 2: + exit('Usage: python visualize_results.py ') + else: + filename = sys.argv[1] + plot_axial_potential(filename) +