diff --git a/Readme.md b/Readme.md
index e2ebb54..0f493ae 100644
--- a/Readme.md
+++ b/Readme.md
@@ -4,7 +4,7 @@ This is a python wrapper for [MECSim](http://www.garethkennedy.net/MECSim.html)
It works completely in python in a Linux environment. I wrote this while working on [GPCV](https://github.com/kiranvad/gpcv) related work.
If you use this software in your work please cite the original MECSim software along with this repository:
-```
+```bibtex
@misc{pymecsim,
author = {Kiran Vaddi},
title = {{pyMECSim: A Python wrapper for MECSim}},
@@ -15,17 +15,17 @@ If you use this software in your work please cite the original MECSim software a
url = {https://github.com/kiranvad/pyMECSim}
}
```
-
-This repository is arranged as follows:
-1. src : Contains the original MECSim software distributed under the same license as MECSim
-2. notebooks : An example usage of pyMECSim is shown.
-3. mechanisms : Folder where you can host a original mechanism as an input file
+To install as a package, run
+```bash
+pip install git+https://github.com/kiranvad/pyMECSim#egg=pyMECSIM.`
+```
+Dependencies will be checked and installed from the setup.py file.
A sample usage is as follows:
Import `pymecsim` using the following:
```python
-from src.pymecsim import MECSIM, pysed, plotcv
+from pymecsim import MECSIM, pysed
```
We can perform a simulation on a one electron transfer mechanism and visualize the effect of changing the formal potential using the following code:
@@ -42,8 +42,8 @@ dirname = os.getcwd()
for i,e0 in enumerate(E0):
outfile = dirname + '/outfile.sk'
pysed('$E0', str(e0), configfile, outfile)
- out = MECSIM(outfile)
- ax = plotcv(out['current'],out['voltage'], ax = ax)
+ model = MECSIM(outfile)
+ ax = model.plot(ax = ax)
ax.set_label("E0 = "+str(e0))
plt.legend([r'$E_0=0.5$',r'$E_0=0.1$',r'$E_0=1e-2$'],loc='lower right')
#plt.savefig('cvexample.png',dpi=500,bbox_inches='tight')
@@ -51,17 +51,23 @@ plt.show()
```
This will plot the following:
-
+
Naturally, you would want to be able to run simulations for different mechanisms and confgurations. You can do that by just defining a mechanism that MECSim can model(see examples for some possible reaction mechanisms [here](http://www.garethkennedy.net/MECSimScripts.html)).
Once you have the mechanism file in say `/path/to/folder/mechanism.sk` format, turn it in as an input to pyMECSim using the following:
```python
-out = MECSIM('/path/to/folder/mechanism.sk')
-plotcv(out['current'],out['voltage'])
+model = MECSIM('/path/to/folder/mechanism.sk')
+model.plot()
+plt.show()
```
+One can also get concentration profiles by first indicating `MECSIM` to return concentration profiles in the configuration file by setting `1 ! show debug output files as well as MECSimOutput.txt (1=yes; 0=no)`. `pymecsim` will then be able to return concentration profiles as numpy arrays. see `notebooks/Cyclic Voltammetry Simulation Example for Single Electron Transfer Mechanism.ipynb` for an example use case.
+
+
+## Notes
+Please free to contribute to this repository both interms of code and documetation or simple example use cases in jupyter notebook. Submit a pull request and I would be happy to integrate into this repository.
diff --git a/__init__.py b/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/mechanisms/cvexamples.sk b/mechanisms/cvexamples.sk
index 4472c4d..dbffdac 100644
--- a/mechanisms/cvexamples.sk
+++ b/mechanisms/cvexamples.sk
@@ -15,7 +15,7 @@
10.0 ! Dstar_min
0.005e0 ! max voltage step
25.6e0 ! time resolution experimentally to correct vscan/f (us)
-0 ! show debug output files as well as MECSimOutput.txt (1=yes; 0=no)
+1 ! show debug output files as well as MECSimOutput.txt (1=yes; 0=no)
0 ! use advanced voltage ramp (0 = E_start=E_end, 1 = use advanced ramp below, 2=From file "EInput.txt")
2 ! number of E_rev lines for advanced ramp (if enter 0 then first E_rev value is the final time
0.0 ! E_start (V)
diff --git a/notebooks/.ipynb_checkpoints/Cyclic Voltammetry Simulation Example for Single Electron Transfer Mechanism-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/Cyclic Voltammetry Simulation Example for Single Electron Transfer Mechanism-checkpoint.ipynb
index 6583564..027fca5 100644
--- a/notebooks/.ipynb_checkpoints/Cyclic Voltammetry Simulation Example for Single Electron Transfer Mechanism-checkpoint.ipynb
+++ b/notebooks/.ipynb_checkpoints/Cyclic Voltammetry Simulation Example for Single Electron Transfer Mechanism-checkpoint.ipynb
@@ -5,7 +5,7 @@
"metadata": {},
"source": [
"## Using MECSim\n",
- "To be finished soon"
+ "In this example we look at a simple example case of single electron trasnfer mechanism. We simulate Cyclic Voltametry response of this particualr mechanism with a planar electrode assumption."
]
},
{
@@ -14,27 +14,12 @@
"metadata": {},
"outputs": [],
"source": [
- "%load_ext autoreload\n",
- "%autoreload 2"
+ "from pymecsim import pysed, MECSIM"
]
},
{
"cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "import sys\n",
- "if '../' not in sys.path:\n",
- " sys.path.insert(0,'../')\n",
- "import os\n",
- "\n",
- "from src.pymecsim import pysed, MECSIM, plotcv"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 3,
+ "execution_count": 2,
"metadata": {},
"outputs": [
{
@@ -74,14 +59,14 @@
" \n",
" \n",
+ "\" id=\"mb92f4f7d36\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n",
" \n",
" \n",
- " \n",
" \n",
" \n",
- " \n",
+ " \n",
" \n",
" \n",
" \n",
" \n",
- " \n",
" \n",
" \n",
- " \n",
+ " \n",
" \n",
" \n",
" \n",
" \n",
- " \n",
" \n",
" \n",
- " \n",
+ " \n",
" \n",
" \n",
" \n",
- " \n",
" \n",
" \n",
" \n",
" \n",
" \n",
- " \n",
+ " \n",
" \n",
" \n",
" \n",
@@ -1192,6 +1155,9 @@
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"%config InlineBackend.figure_formats = ['svg']\n",
+ "import os\n",
+ "import warnings\n",
+ "warnings.filterwarnings(\"ignore\")\n",
"\n",
"configfile = '../mechanisms/cvexamples.sk'\n",
"E0 = [-0.25,0.0,0.25]\n",
@@ -1201,16 +1167,2369 @@
"for i,e0 in enumerate(E0):\n",
" outfile = dirname + '/outfile.sk'\n",
" pysed('$E0', str(e0), configfile, outfile)\n",
- " out = MECSIM(outfile)\n",
- " T = out['time']\n",
+ " model = MECSIM(outfile)\n",
+ " T = model.T\n",
" forward_sweep = np.arange(len(T)/2,len(T)).astype(int)\n",
- " ax = plotcv(out['current'],out['voltage'], ax = ax)\n",
+ " model.plot(ax=ax)\n",
" ax.set_label(\"E0 = \"+str(e0))\n",
"plt.legend([r'$E_0=0.5$',r'$E_0=0.1$',r'$E_0=1e-2$'],loc='lower right')\n",
"#plt.savefig('cvexample.png',dpi=500,bbox_inches='tight')\n",
"plt.show()"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can visualize the concentration profiles of the species in the bulk and at the surface. In this example, we simulated `A +1e = B` thus we have two species. Since this is a simple one step reaction mechanism and diffusin driven, we should see that increase in concentration of species `A` should lead to decrease in concentration of species `B` at the surface where the electron transfer takes place. Similarily in the bulk, one should observe similar behavior upuntil a `diffusion layer` length scale."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "species = [r'$A$',r'$B$'] # in the same order as the mechanism\n",
+ "N = len(species)\n",
+ "C_bulk = model.get_bulk_concentrations()\n",
+ "\n",
+ "fig, ax = plt.subplots()\n",
+ "for i in range(N):\n",
+ " ax.plot(C_bulk[:,0], C_bulk[:,i+1], label=species[i])\n",
+ "ax.set_xlabel('Normalized length')\n",
+ "ax.set_ylabel('Normalized concentration')\n",
+ "ax.legend()\n",
+ "ax.set_title('Species concentration in the bulk (Final)')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "species = [r'$A$',r'$B$'] # in the same order as the mechanism\n",
+ "N = len(species)\n",
+ "C_surface = model.get_surface_concentrations()\n",
+ "\n",
+ "fig, ax = plt.subplots()\n",
+ "for i in range(N):\n",
+ " ax.plot(model.T, C_surface[1:,i], label=species[i])\n",
+ "ax.set_xlabel('time (s)')\n",
+ "ax.set_ylabel('Normalized concentration')\n",
+ "ax.legend()\n",
+ "ax.set_title('Species concentration at the surface')\n",
+ "plt.show()"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
@@ -1235,7 +3554,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.3"
+ "version": "3.7.7"
}
},
"nbformat": 4,
diff --git a/notebooks/Cyclic Voltammetry Simulation Example for Single Electron Transfer Mechanism.ipynb b/notebooks/Cyclic Voltammetry Simulation Example for Single Electron Transfer Mechanism.ipynb
index 6583564..027fca5 100644
--- a/notebooks/Cyclic Voltammetry Simulation Example for Single Electron Transfer Mechanism.ipynb
+++ b/notebooks/Cyclic Voltammetry Simulation Example for Single Electron Transfer Mechanism.ipynb
@@ -5,7 +5,7 @@
"metadata": {},
"source": [
"## Using MECSim\n",
- "To be finished soon"
+ "In this example we look at a simple example case of single electron trasnfer mechanism. We simulate Cyclic Voltametry response of this particualr mechanism with a planar electrode assumption."
]
},
{
@@ -14,27 +14,12 @@
"metadata": {},
"outputs": [],
"source": [
- "%load_ext autoreload\n",
- "%autoreload 2"
+ "from pymecsim import pysed, MECSIM"
]
},
{
"cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "import sys\n",
- "if '../' not in sys.path:\n",
- " sys.path.insert(0,'../')\n",
- "import os\n",
- "\n",
- "from src.pymecsim import pysed, MECSIM, plotcv"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 3,
+ "execution_count": 2,
"metadata": {},
"outputs": [
{
@@ -74,14 +59,14 @@
" \n",
" \n",
+ "\" id=\"mb92f4f7d36\" style=\"stroke:#000000;stroke-width:0.8;\"/>\n",
" \n",
" \n",
- " \n",
" \n",
" \n",
- " \n",
+ " \n",
" \n",
" \n",
" \n",
" \n",
- " \n",
" \n",
" \n",
- " \n",
+ " \n",
" \n",
" \n",
" \n",
" \n",
- " \n",
" \n",
" \n",
- " \n",
+ " \n",
" \n",
" \n",
" \n",
- " \n",
" \n",
" \n",
" \n",
" \n",
" \n",
- " \n",
+ " \n",
" \n",
" \n",
" \n",
@@ -1192,6 +1155,9 @@
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"%config InlineBackend.figure_formats = ['svg']\n",
+ "import os\n",
+ "import warnings\n",
+ "warnings.filterwarnings(\"ignore\")\n",
"\n",
"configfile = '../mechanisms/cvexamples.sk'\n",
"E0 = [-0.25,0.0,0.25]\n",
@@ -1201,16 +1167,2369 @@
"for i,e0 in enumerate(E0):\n",
" outfile = dirname + '/outfile.sk'\n",
" pysed('$E0', str(e0), configfile, outfile)\n",
- " out = MECSIM(outfile)\n",
- " T = out['time']\n",
+ " model = MECSIM(outfile)\n",
+ " T = model.T\n",
" forward_sweep = np.arange(len(T)/2,len(T)).astype(int)\n",
- " ax = plotcv(out['current'],out['voltage'], ax = ax)\n",
+ " model.plot(ax=ax)\n",
" ax.set_label(\"E0 = \"+str(e0))\n",
"plt.legend([r'$E_0=0.5$',r'$E_0=0.1$',r'$E_0=1e-2$'],loc='lower right')\n",
"#plt.savefig('cvexample.png',dpi=500,bbox_inches='tight')\n",
"plt.show()"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can visualize the concentration profiles of the species in the bulk and at the surface. In this example, we simulated `A +1e = B` thus we have two species. Since this is a simple one step reaction mechanism and diffusin driven, we should see that increase in concentration of species `A` should lead to decrease in concentration of species `B` at the surface where the electron transfer takes place. Similarily in the bulk, one should observe similar behavior upuntil a `diffusion layer` length scale."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "species = [r'$A$',r'$B$'] # in the same order as the mechanism\n",
+ "N = len(species)\n",
+ "C_bulk = model.get_bulk_concentrations()\n",
+ "\n",
+ "fig, ax = plt.subplots()\n",
+ "for i in range(N):\n",
+ " ax.plot(C_bulk[:,0], C_bulk[:,i+1], label=species[i])\n",
+ "ax.set_xlabel('Normalized length')\n",
+ "ax.set_ylabel('Normalized concentration')\n",
+ "ax.legend()\n",
+ "ax.set_title('Species concentration in the bulk (Final)')\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "species = [r'$A$',r'$B$'] # in the same order as the mechanism\n",
+ "N = len(species)\n",
+ "C_surface = model.get_surface_concentrations()\n",
+ "\n",
+ "fig, ax = plt.subplots()\n",
+ "for i in range(N):\n",
+ " ax.plot(model.T, C_surface[1:,i], label=species[i])\n",
+ "ax.set_xlabel('time (s)')\n",
+ "ax.set_ylabel('Normalized concentration')\n",
+ "ax.legend()\n",
+ "ax.set_title('Species concentration at the surface')\n",
+ "plt.show()"
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
@@ -1235,7 +3554,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.7.3"
+ "version": "3.7.7"
}
},
"nbformat": 4,
diff --git a/notebooks/log.txt b/notebooks/log.txt
new file mode 100644
index 0000000..e5d33c9
--- /dev/null
+++ b/notebooks/log.txt
@@ -0,0 +1,116 @@
+
+ *****************************************************
+ MECSim - GK 7/08/2016
+ *****************************************************
+
+
+
+ General parameters:
+ Temperature = 298.2000
+ Uncompensated R = 0.0
+
+ Voltage ramp (V):
+ E_start = 0.5000
+ E_rev = -0.5000
+ E_end = 0.5000
+ Cycles = 1
+ Total range = 2.0000
+
+ Scan rate (V/s) = 1.00000000E+01
+ Final time (s) = 2.00000000E-01
+
+Found 1 AC signal(s). Max frequency = 0.0000 Hz
+ 0.0000 mV at 180.0000 Hz
+
+ No capacitor requested
+
+
+ Found a total of 2 species
+
+ No pre-equilibriation; use entered concentrations
+
+ Summary of solution reactions:
+ Charge transfer = 1
+ 1st Order Chem = 0
+ 2nd Order Chem = 0
+ Summary of surface confined reactions:
+ Charge transfer = 0
+ 1st Order Chem = 0
+ 2nd Order Chem = 0
+
+ Charge transfer reaction(s):
+ A + 1e = B ; E0 = 0.250 , ks = 1.0000E+04 cm/s , alpha = 0.50
+ ; K_eqm = 1.6794E+04
+
+ Electrode geometry:
+ Planar with area = 1.00000000E+00 cm^2
+
+ Technical details of simulation
+ Using Butler-Volmer Theory
+ Points in time = 2^14 = 16384
+ delt = 1.2207E-05 and delE = 1.2207E-04
+ Exponential spatial grid with 57 points and delx(0) = 3.4939E-06
+ x = 1.7033E-06 5.3763E-06 9.4355E-06 1.3922E-05 | 9.4112E-03
+
+ Solution phase:
+ Initial scaled concentrations and diffusion coeff (cm2/s)
+ [A] = 1.00000000E-06 ; D = 1.00000000E-05
+ (1.00000000)
+ [B] = 0.00000000E+00 ; D = 1.00000000E-05
+ (0.00000000)
+ Concentrations scaled by 1.00000000E-06 mol/cm3
+ 1.00000000E+00 mM
+
+
+ ************** Starting Time Loop ******************
+
+ Done 10%; t,Eapp,i = 2.00E-02 3.00E-01 -7.11E-04, C_err= -2.87E-12
+ 0.8752 0.1248 0.0000 0.0000 0.0000 0.0000
+ Done 20%; t,Eapp,i = 4.00E-02 1.00E-01 -1.47E-03, C_err= -1.78E-12
+ 0.0029 0.9971 0.0000 0.0000 0.0000 0.0000
+ Done 30%; t,Eapp,i = 6.00E-02 -1.00E-01 -9.27E-04, C_err= -1.26E-12
+ 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000
+ Done 40%; t,Eapp,i = 8.00E-02 -3.00E-01 -7.36E-04, C_err= -1.42E-12
+ 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000
+ Done 50%; t,Eapp,i = 1.00E-01 -5.00E-01 -6.29E-04, C_err= -1.62E-12
+ 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000
+ Done 60%; t,Eapp,i = 1.20E-01 -3.00E-01 -5.59E-04, C_err= -1.83E-12
+ 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000
+ Done 70%; t,Eapp,i = 1.40E-01 -1.00E-01 -5.08E-04, C_err= -2.06E-12
+ 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000
+ Done 80%; t,Eapp,i = 1.60E-01 1.00E-01 -4.51E-04, C_err= -3.06E-12
+ 0.0029 0.9971 0.0000 0.0000 0.0000 0.0000
+ Done 90%; t,Eapp,i = 1.80E-01 3.00E-01 2.10E-03, C_err= -6.10E-13
+ 0.8747 0.1253 0.0000 0.0000 0.0000 0.0000
+ Done 100%; t,Eapp,i = 2.00E-01 5.00E-01 6.93E-04, C_err= -2.34E-12
+ 0.9999 0.0001 0.0000 0.0000 0.0000 0.0000
+
+ Concentration min/max values:
+ 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
+ 1.0000 1.0000 0.0000 0.0000 0.0000 0.0000
+
+ ****************************************************
+
+ Current minimum = -2.68622429E-03 and maximum = 2.24577576E-03
+ Voltage minimum = 2.21557617E-01 and maximum = 2.78686523E-01
+
+ Average number of iterations at electrode surface = 1.000
+ Average number of iterations for chemical reactions = 1.000
+ Average ODE calculations for iterating surface conf = 0.000
+ Average ODE calculations for bounding surface conf = 0.000
+
+ Final scaled concentrations on spatial grid
+ [A] = 9.99940454E-01 9.98716251E-01 9.96076410E-01 | 9.99998347E-01
+ [B] = 5.95460837E-05 1.28374882E-03 3.92359046E-03 | 1.65337954E-06
+
+ Error in concentration sum (=0 if D_A = D_B etc)
+ -2.34245956E-12 -2.34312569E-12 -2.34601227E-12 | 6.50146603E-13
+
+
+ Output file written to MECSimOutput_Pot.txt
+
+
+ Additional files written to
+ EC_Model.tvc - time, Eapp and concentrations for all time steps (first number is #species)
+ EC_Model.fin - final concentrations against distance from electrode
+
diff --git a/notebooks/outfile.sk b/notebooks/outfile.sk
index 7eb35de..a060cdc 100644
--- a/notebooks/outfile.sk
+++ b/notebooks/outfile.sk
@@ -15,7 +15,7 @@
10.0 ! Dstar_min
0.005e0 ! max voltage step
25.6e0 ! time resolution experimentally to correct vscan/f (us)
-0 ! show debug output files as well as MECSimOutput.txt (1=yes; 0=no)
+1 ! show debug output files as well as MECSimOutput.txt (1=yes; 0=no)
0 ! use advanced voltage ramp (0 = E_start=E_end, 1 = use advanced ramp below, 2=From file "EInput.txt")
2 ! number of E_rev lines for advanced ramp (if enter 0 then first E_rev value is the final time
0.0 ! E_start (V)
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..6bc0d93
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,15 @@
+from setuptools import setup,find_packages
+import sys, os
+
+setup(name="pymecsim",
+ description="Python Cyclic Voltammetry Simulator",
+ version='1.0',
+ author='Kiran Vaddi',
+ author_email='kiranvad@buffalo.edu',
+ license='MIT',
+ python_requires='>=3.6',
+ install_requires=['numpy','scipy', 'pandas'],
+ extras_require = {},
+ packages=find_packages(),
+ long_description=open('Readme.md').read()
+)
\ No newline at end of file
diff --git a/src/__init__.py b/src/__init__.py
index 6e45a10..ebea8f6 100644
--- a/src/__init__.py
+++ b/src/__init__.py
@@ -1 +1,3 @@
-from .pymecsim import *
\ No newline at end of file
+#from .pymecsim import *
+from .utils import pysed
+from .core import MECSIM
\ No newline at end of file
diff --git a/src/__pycache__/pymecsim.cpython-37.pyc b/src/__pycache__/pymecsim.cpython-37.pyc
deleted file mode 100644
index fd65b3a..0000000
Binary files a/src/__pycache__/pymecsim.cpython-37.pyc and /dev/null differ
diff --git a/src/core.py b/src/core.py
new file mode 100644
index 0000000..65ac5cf
--- /dev/null
+++ b/src/core.py
@@ -0,0 +1,206 @@
+import subprocess, shlex
+import pdb, os, sys
+import numpy as np
+import pandas as pd
+import os
+dirname = os.path.dirname(__file__)
+
+from shutil import copyfile
+import warnings
+
+import seaborn as sns
+import matplotlib.pyplot as plt
+from matplotlib import rc
+rc('text', usetex=True)
+
+class MECSIM:
+ def __init__(self, configfile):
+ """
+ Simulate a voltammetry response for a given configuration.
+
+ Input:
+ -----
+ configfile : An MECSIM configuration file usually in the .inp form
+
+ Attributes:
+ -----------
+ T : Time profile
+ I : Current profile
+ V : Voltage profile
+
+ Note
+ -----
+ This class is a python wrapper for the MECSIM software by Gareth F Kennerdy (http://www.garethkennedy.net/MECSim.html).
+ It simulates a mechanism with the configurations specified and throws errors or warnings based on the results of simulations.
+ While the python wrapper provides access to many functions of MECSIM, it sometimes fails to understand the errors from MECSIM.
+ Best practice is to always look at the log.txt file generated in your working directory and figure out what has gone wrong.
+ """
+ self.configfile = configfile
+ self.dirname = os.path.dirname(__file__)
+ self.solve()
+ self._get_errors()
+ self.T,self.V,self.I = self._read_mecsim_out(self._get_filename('./'+self.outfile))
+
+ def solve(self):
+ main_configfile = os.path.join(self.dirname, 'Master.inp')
+ copyfile(self.configfile, main_configfile)
+
+ args = shlex.split('chmod u+x ./MECSim')
+ process = subprocess.Popen(args,stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=dirname)
+ args = shlex.split('./MECSim')
+ process = subprocess.Popen(args,stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=dirname)
+ self.stdout, self.stderr = process.communicate()
+
+ logfile = os.path.join(os.getcwd(),'log.txt')
+ copyfile(os.path.join(dirname, 'log.txt'), logfile)
+ file_path = self._get_filename('log.txt')
+ f = open(file_path, 'r')
+ self.logs = list(f)
+ self.flag_concs = False
+ for line in reversed(self.logs):
+ if 'Output file' in line:
+ self.outfile = line.split(' ')[-1].split('\n')[0]
+
+ elif 'Additional files' in line:
+ self.flag_concs= True
+
+ def _get_errors(self):
+ for line in reversed(self.logs):
+ if "Error" in line:
+ error = line.replace('Error: ','')
+ break
+
+ if not len(self.stderr)==0:
+ raise RuntimeError('MECSIM threw the following error that pyMECSIM could not understand \n'+
+ self.stderr.decode("utf-8"))
+
+ if self.stderr.decode("utf-8") == 'STOP 4\n':
+ raise RuntimeError('pyMECSIM could not understand the error message from MECSIM.')
+
+ if len(self.logs)<30:
+ if 'error' in locals():
+ raise RuntimeError(error)
+ else:
+ raise RuntimeError('Could not identify a clear error from MECSIM but here is what is avaiable: \n'
+ + logs[-1])
+ else:
+ if 'error' in locals():
+ warnings.warn(error)
+
+ def get_current_profile(self):
+ """ return current values over simulated time scale """
+ return self.I
+ def get_voltage_profile(self):
+ """
+ return voltage applied over simulated time scale.
+ Read MECSIM paper for more details on what is applied voltage.
+ https://doi.org/10.1016/j.coelec.2016.12.001
+
+ """
+ return self.V
+
+ def get_time_profile(self):
+ """
+ return time profile of the simulation.
+ """
+ return self.T
+
+ def _read_mecsim_out(self, filename):
+ """ internal function to read the text output """
+ f = open(filename, 'r')
+ # search for last line of header that is made by MECSim (always this line)
+ for line in f:
+ if line.strip() == "Post(ms): 0.000000E+00": break
+ time = []
+ eapp = []
+ current = []
+ for line in f:
+ columns = line.split()
+ time.append(float(columns[2]))
+ eapp.append(float(columns[0]))
+ current.append(float(columns[1]))
+
+ return np.asfarray(time), np.asfarray(eapp), np.asfarray(current)
+
+ def get_bulk_concentrations(self):
+ """
+ returns bulk concentrations C(x,t=t_end).
+
+ output is a numpy array of shape (N,P) where N is number of grid points in x-direction and P is number of species
+ """
+ if not self.flag_concs:
+ raise KeyError('MECSIM did not output any concentration profiles.'+
+ 'Did you specifiy in the configuration file to output the concentration profiles?')
+
+ fin_file = self._get_filename('EC_Model.fin')
+ f = open(fin_file, 'r')
+ concs = []
+ for i, line in enumerate(f):
+ columns = line.split()
+ concs.append(columns)
+ concs = np.asarray(concs)
+ concs = self._get_float_matrix(concs)
+
+ return concs
+
+ def get_surface_concentrations(self):
+ """
+ returns surface concentrations profiles over time C(x=0,t).
+
+ output is a numpy array of shape (N,P) where N is number of time points and P is number of species
+ """
+ if not self.flag_concs:
+ raise KeyError('MECSIM did not output any concentration profiles.'+
+ 'Did you specifiy in the configuration file to output the concentration profiles?')
+ tcv_file = self._get_filename('EC_Model.tvc')
+ f = open(tcv_file, 'r')
+ concs = []
+ for i, line in enumerate(f):
+ if i>0:
+ columns = line.split()
+ concs.append(columns[3:])
+ concs = np.asarray(concs)
+ concs = self._get_float_matrix(concs)
+
+ return concs
+
+ def plot(self,ax = None):
+ """ simple visualization tool for CV curve """
+ if ax is None:
+ fig = plt.figure(figsize = (4,4))
+ ax = fig.add_subplot(111)
+ ax.plot(self.V,self.I)
+
+ ax.set_xlabel('Voltage(V)')
+ ax.set_ylabel('Current(A)')
+
+ plt.tight_layout()
+ sns.despine()
+
+ return ax
+
+ def _get_filename(self,filename):
+ file_path = os.path.join(dirname, './', filename)
+
+ if os.path.exists(file_path):
+ return file_path
+ else:
+ return False
+
+ def _get_float_matrix(self, matrix):
+ """
+ Given a numpy matrix of floats, returns a float matrix
+ This is useful because sometimes you get a string that can not be made to a float
+ """
+ def to_float(string):
+ try:
+ f = float(string)
+ except ValueError:
+ f = string.replace("-", 'E-')
+
+ return f
+
+ df = pd.DataFrame(matrix)
+ df = df.applymap(to_float)
+
+ return df.to_numpy()
\ No newline at end of file
diff --git a/src/pymecsim.py b/src/pymecsim.py
deleted file mode 100644
index 1f4efbf..0000000
--- a/src/pymecsim.py
+++ /dev/null
@@ -1,117 +0,0 @@
-import subprocess, shlex
-import pdb
-import os
-import sys
-if '../' not in sys.path:
- sys.path.insert(0,'../')
-
-import os
-dirname = os.path.dirname(__file__)
-
-from shutil import copyfile
-
-def MECSIM(configfile='Master.sk'):
- """
- This is a python wrapper for the original MECSIM software re-distrubuted along with this packages as a .exe file
- It takes an MECSim configuration file in the configfile as input and returns a dictonarry.
-
- INPUT
- -----
- configfile : an MECSim configuration file (see ../mechanisms/cvexamples.sk for an example)
- For advanced usage, refer to the original MECSim user manual
-
- OUTPUT
- ------
- output : a dictonary with the 'current', 'voltage', 'time' and 'info':output from MECSim software
-
-
- NOTES:
- ------
- A log.txt file will be created in the `src` folder with detailed output from MECSim
- (in future this willbe moved to be as an output file)
-
- """
- main_configfile = os.path.join(dirname, 'Master.inp')
- copyfile(configfile, main_configfile)
-
- args = shlex.split('chmod u+x ./MECSim')
- process = subprocess.Popen(args,stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=dirname)
- args = shlex.split('./MECSim')
- process = subprocess.Popen(args,stdout=subprocess.PIPE, stderr=subprocess.PIPE, cwd=dirname)
- stdout, stderr = process.communicate()
-
- T,V,I = ReadMECSimOut(os.path.join(dirname, './MECSimOutput_Pot.txt'))
-
- output = {'current':I,'voltage':V,'time':T,'info':[stdout,stderr]}
-
- return output
-
-
-import re
-
-def pysed(oldstr, newstr, infile, outfile):
- """
- Sed-like Replace function..
- Usage: pysed.replace(, , )
- Example: pysed.replace('xyz', 'XYZ', '/path/to/file.txt')
- Example 'DRYRUN': pysed.replace('xyz', 'XYZ', '/path/to/file.txt', dryrun=True)
- #This will dump the output to STDOUT instead of changing the input file.
-
- A Python equivalent of sed from bash:
- source : https://github.com/mahmoudadel2/pysed/blob/master/pysed.py
-"""
-
- linelist = []
- with open(infile) as f:
- for item in f:
- newitem = item.replace(oldstr, newstr)
- linelist.append(newitem)
- with open(outfile, "w") as f:
- f.truncate()
- for line in linelist: f.writelines(line)
- f.writelines('\n')
-
-
-import numpy as np
-
-def ReadMECSimOut(filename):
- """
- Internal function used to read output from MECSim to simple numpy arrays
- """
- f = open(filename, 'r')
- # search for last line of header that is made by MECSim (always this line)
- for line in f:
- if line.strip() == "Post(ms): 0.000000E+00": break
- time = []
- eapp = []
- current = []
- for line in f:
- columns = line.split()
- time.append(float(columns[2]))
- eapp.append(float(columns[0]))
- current.append(float(columns[1]))
-
- return np.asfarray(time), np.asfarray(eapp), np.asfarray(current)
-
-import seaborn as sns
-import matplotlib.pyplot as plt
-from matplotlib import rc
-rc('text', usetex=True)
-
-def plotcv(I,V,ax = None):
- """
- A utility function to plot CV curve. I is output['current'] and V is output['voltage']
- """
- if ax is None:
- fig = plt.figure(figsize = (4,4))
- ax = fig.add_subplot(111)
- ax.plot(V,I)
-
- ax.set_xlabel('Voltage(V)')
- ax.set_ylabel('Current(mA)')
-
- plt.tight_layout()
- sns.despine()
-
- return ax
-
\ No newline at end of file
diff --git a/src/utils.py b/src/utils.py
new file mode 100644
index 0000000..f60e4ab
--- /dev/null
+++ b/src/utils.py
@@ -0,0 +1,39 @@
+import re
+
+def pysed(oldstr, newstr, infile, outfile):
+ """
+ A utility function to modify a skeleton configuration file iteratively.
+ See examples for a sample usage.
+
+ Inputs:
+ --------
+ oldstr : string uses as a variable in .sk file (usually starts with a $ sign)
+ newstr : string used to replace the `oldstr` variable (although its float one should pass a float)
+ infile : file with variables defined
+ outfile : file with variables replaced by the floats we want.
+
+
+ --------------------------------------------------------------------
+ Sed-like Replace function..
+ Usage: pysed.replace(, , )
+ Example: pysed.replace('xyz', 'XYZ', '/path/to/file.txt')
+ Example 'DRYRUN': pysed.replace('xyz', 'XYZ', '/path/to/file.txt', dryrun=True)
+ #This will dump the output to STDOUT instead of changing the input file.
+
+ A Python equivalent of sed from bash:
+ source : https://github.com/mahmoudadel2/pysed/blob/master/pysed.py
+"""
+
+ linelist = []
+ with open(infile) as f:
+ for line in f:
+ if oldstr in line:
+ newline = line.replace(oldstr, newstr)
+ linelist.append(newline)
+ else:
+ linelist.append(line)
+
+ with open(outfile, "w") as f:
+ f.truncate()
+ for line in linelist: f.writelines(line)
+ f.writelines('\n')
\ No newline at end of file
+