# Charge regulation Project

###### tags: `charge regulation`

---

---

### Owners (the only one with the permission to edit the main test)

Marielle, Etienne, Frank, Rudy, Giuseppe

---

---

## Background

In this project we are going to build a coarse grained model of a particle system with a LJ interactions plus charge regulation term. We will implement a MC algorithm to solve the problem. The potential energy of the system is given by:

$$
U(\{\mathbf{r}_i, q_i\}) = U_{ LJ}(\{\mathbf{r}_i\}) + U_{ fe} (\{\mathbf{r}_i, q_i\})
$$
the first term is a standard LJ term:
$$
U_{ LJ}(\{\mathbf{r}_i\})=\sum_{i<j} u(|\mathbf{r}_i-\mathbf{r}_j|)
$$

with:

$$
u(r)=4 \epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]
$$

The second one is the main on for the charge fluctuations is given by:

$$
U_{fe} (\{\mathbf{r}_i, q_i\}) = {\textstyle\frac12} \sum_{i=1}^{N} {{\cal C}_i}^{-1} \left( q_i - \overline{q}_i\right)^2 + {\textstyle\frac12} \sum_{i\neq j}^{N} q_i q_j V({\bf r}_i, {\bf r}_j)
$$

with the electrostatic expressed in the Debye-Huckel form:

$$
V({\bf r}_i, {\bf r}_j) = \frac{1}{4\pi\varepsilon\varepsilon_0} \frac{e^{- \kappa \vert {\bf r}_i - {\bf r}_j\vert}}{\vert {\bf r}_i - {\bf r}_j\vert}.
$$

The main paramenters are ${\cal C}_i$ and $\overline{q}_i$ the capacitance and the mean charge.

## Plans

### NVT Monte Carlo Simulation

We will implement a MC method to have the NVT equilibrium of the model. As a first test mode we will investigate the case of zero average charge, i.e. $\overline{q}_i = 0$ and screening length $\kappa^{-1} =0.3$.

Units:

LJ: $\sigma = 1$, $\epsilon = 1$
fe: $\frac{1}{4\pi\varepsilon\varepsilon_0}=1$

## To Do

:::success
green is for finished tasks
:::
### NVT code for simple LJ
:::success
* Implement a working code with cell list
* Test (P and $\langle U\rangle$) the code for different density and temperature
* Test the performance
* Radial distribution function (Marielle you will need to learn about this )
:::

### NPT code for simple LJ
:::success
* Implement a working code with cell list now for NPT
* Verify that you obtain the same results as above
:::

### Gibbs-Ensemble (GE) code for simple LJ
:::success
* Implement a working code with cell list now for GE
* Test the LLPS in LJ
:::

### NVT code for LJ+fe

:::success
* Add the *charge move* into the NVT algorithm
* Measure the $\langle U\rangle$ and P for different ${\cal C}_i$
* Compare to GF results
:::
* Produce data set for $T=1.25$ and $T=2.0$ and $\rho=0.2, 0.4, 0.6, 0.8$ using a fine grid in ${\cal C}_i$.
* Density structure: RDF
* Charge structure: charge distribution

### OPEN QUESTIONS

* Can we look at LLPS for LJ+CHARGE?
* Equation of state (Rudy's note)?
* Polymer

## RESULTS

### simple LJ, NVT
* Lennard Jones potential truncated at $r_c = 2.5$ and shifted by its value at $r_c$
* simulations performed at $T = 2.0$, $T=1.5$, $T=1.2$ (above $T_c$) to validate the code regarding P and U of the system



### simple LJ, NPT


### Gibbs ensemble
system of $N = N_1+N_2= 400$ particles
simulation at 10 different temperatures $T \in [0.64, 1.0]$

* densities for gas and liquid phase predicted correctly (even though larger deviations for the liquid phase)


* pressures underly large fluctuations
* equal pressure values in both phases within the errors
*

# charge fluctuations

* N = 343 particles
* simulations at $\rho = 0.4$; $T= 2.0$ and $T= 1.0$
* truncate the screened Coulomb potential at $r = 2.5$
* fluctuating charges $q \in [-10,10]$


* I will try to assemble some more points within 0.01 and 10


* As far as I remebmer, T=1.0 was the temperature for the plot you showed during the last meeting;

* I'm not too sure whether to trust the simulation at high $C^{-1}$ as this should be the simple-Lennard-Jones limit but for that, $T=1.0$ is below the critical temperature

* I did sth stupid by accidentally only picking even charges during the MC moves; I will correct that and do it again with all charge values allowed

> [name=gfoffi] It seems that we get the same result!
>


* **Energies**

| | |
| -------- | -------- |
|  |  |

* **Separated Energy contributions at $\rho = 0.2$**

| | |
| -------- | -------- |
|  |  |

* **Capacitance Term at $Ţ = 1.25$**

| | |
| --- | --- |
|  |  |

The bump at the left looks a bit weird; it gets bigger with inreasing density; is in the range where Coulomb starts to take over

### Prediction of the shape of the "bump" in the energy at intermediate C

>[name=frank] Quick note on this from me, feel free to move it or reintegrate it somewhere else if that's better.

In the plots above, and in Giuseppe's preliminary data, we see that at intermediate $C^{-1}$, there is an increase in energy due to the charging of the particles. This contribution is likely mainly due to the part of the Hamiltonian that scales with $C^{-1}$. To see the expected shape of it, consider the expected behavior of this part of the energy for an isolated particle. In that case, the energy of the particle is just $U = \frac12{\cal C}^{-1} q^2$ (where we have set $\overline{q}$ to 0, which does not affect the result).

The partition function is:
$$Z = \sum_{q=-q_\mathrm{min}}^{q_\mathrm{max}} \exp(-\frac12\beta {\cal C}^{-1} q^2)$$
This is not a nice sum, it seems (Mathematica does not provide me with any nice solutions...), but we can still also write down the expected energy:
$$\langle U \rangle = \frac{{\cal C}^{-1}}{2Z}\sum_{q=-q_\mathrm{min}}^{q_\mathrm{max}} q^2 \exp(-\beta {\cal C}^{-1} q^2)$$
We can plot this numerically:
(with $\beta = 1/k_B T$ and the charge range [$-n_\mathrm{max}, n_\mathrm{max}$])
It might be interesting to plot this together with the total energy (as an addition to the dashed line), or to just compare it to the energy component that comes from this term in the Hamiltonian.

It seems this might be slightly shifted to the right compared to your results... I might have made a typo somewhere.

---

## A definition of charge-RDF

$$\tilde{q}_i = (q_i - \overline{q}) / \sigma_{q}$$

$\sigma_q$ is the standard deviation of the charge at this state point

$\overline{q}$ here also depends on $T$ and $\rho$.

$$g_q(r) = \frac{\sum_{ij}\tilde{q}_i \tilde{q}_j\delta(r_{ij} - r) }{\sum_{ij}\delta(r_{ij} - r)}$$

## Radial distribution functions

* slight upward shift of the distribution compared to simple LJ; this is more pronounced the lower the temperature, the lower the density and the lower $C^{-1}$

| | |
| ------------------------------------ | ------------------------------------ |
|  |  |
| | |
| | |
| | |

## Charge distribution functions

| | |
| -------------------------- | ------------------------------------ |
|  | |
| | |
| | |

* what looks reasonable is that directly around one particle the opposite charge is assembled, further away the same charge

* the behaviour of the CDF at $C^{-1} = 1$ here (and I checked, it seems the same for other intermediate capacitances 0.6, 2,4,6) is not so clear yet for me

| | |
| -------------------------- | ------------------------------------ |
|  | |
| | |
| ||

* histogram of charges: the spread of charges gets wider the higher the capacitance C to carry a charge; at fixed C the higher the density the more you get of the higher charges

| | |
| ------------------------------------ | --- |
|  | |
| | |
|  | |

## Structural analysis

all values for $k_B T / \epsilon = 1.25$


* total charge within shell increases with decreasing $C^{-1}$ as particles can carry more charge

* the radial distance of the shell boundaries stay roughly the same



upper table 4.2: integration up to the first max of the CDF
lower table 4.2: integration up to the first root of CDF (more reasonable for definition of charge shell)
* charge in the first shell increases
* number of particles in the first charge shell decreases which makes sense as one particle can carry more charge
* the positions of the max/min of the charge distribution function get shifted to the left with decreasing $C^{-1}$ meaning slightly narrower first shell

| | |
| -------- | -------- |
| | |

## Phase Diagram

| | | |
| -------- | -------- | -------- |
|  |  |  |
|phase diagram| shifted phase diagram to match the LJ curve| phase diagram with estimate of the critical temperatures |

* for high $C^{-1} \sigma$ seems to work with only a shift; for lower values the coexistence line is distorted, only the left branch suits th LJ one
* shift values

| $C^{-1} \sigma$ | 0.1 | 0.06 | 0.02|
| --------------- | ---- | --- | ---- |
| $\Delta T$ | 0.019 | 0.05 | 0.15 |

* critical temperatures in third pictures estimated with the assumption that the only change btw the curves was the upward shift (of course that is not really accurate)

**capacitance and coulomb contributions**

| | |
| -------- | -------- |
|  |  |

* at high $C^{-1} \sigma$ no charges, at $C^{-1} \sigma = 0.1$ small charges, $q^2/2 \sim 5$ but both phases almost equally charged
* at low $C^{-1} \sigma$ liquid phase carries more charge than gas phase also leading to stronger Coulomb interactions

* can also be seen in charge histograms

| | |
| ------------------------------------ | ------------------------------------ |
|no charges in both phases| |
|||
| $C^{-1} \sigma = 0.1$ gas phase favours low charges, the more the lower the temperatrue | spread to larges charges is a bit higher than for gas phase but still not too pronounced |
|  |  |
|$C^{-1} \sigma = 0.06$ gas phase, charges get higher but still concentrated around 0 | liquid phase goes to higher charges; still more centered around 0 for $k_B T / \epsilon$ high but more or less equally for $k_B T / \epsilon = 0.8$|
|  | |
| highly charged particles at $C^{-1} \sigma = 0.02$, roughly equal distributed in low density phase | and clearly favouring the maximal charge for liquid phase |
| ||

## "real" DLVO model

potential energy

$\beta U(r) = \beta U_{LJ} + \frac{\beta}{2}\sum_{i = 1}^{N} C_i^{-1} (q_i - \bar{q}_i)^2 + \sum_{i, j > i} \lambda_B q_i q_j \frac{e^{-\kappa r_{ij}}}{r_{ij}} \left( \frac{e^{\kappa a}}{1 +\kappa a } \right)^2$

* think now in terms of $k_B T$
* $\lambda_B / \sigma = 0.1$
* particle extension $a = 0.5 \sigma$
* larger system of $N=2000$ particles

### phase diagram at fixed screening lenghts


I cleaned this up a little bit with some more data for the high screening lengths

* the upward shift with increasing $(\kappa \sigma)^{-1}$ seems to be monotoneous after all

* (comments on the big error bars: the cases where boxes changed identities are less now, but it sometimes still happens; I fitted gaussians to the peaks in the histograms and therefore, depending on the peak-quality, errors from the fit are quite large)

* (comment on the pressure sampling: I still have the case that the pressure values differ quite a lot in the phase separated sytem but I thought that might actually be not surprising as we calculate it from the virial in which we have a term depending on the charges; and if, as we have seen, the phases show different charge behaviour, pressures can hardly be the same; the fact that the differences in pressure get more pronounced the further we phase separate might support that assumption;
but it is still a little confusing in terms of two phases at equilibrium not having the same pressure...)

## phase diagram at fixed capacitance


* (some points are still to be updated as the corresponding simulations did not converge that nicely)

* phase separation towards higher screening lengths
* separates earlier the lower the $C^{-1}$

| | |
| -------- | -------- |
|  |  |

* separate energy contribution; similar to what we observed before: higher charge in the liquid phase and more charge at lower $C^{-1}$

| | |
| -------- | -------- |
| | |
| ||

**average q and standard deviation during the simulation**
| | |
| -------- | -------- |
| ||
|||
|||

* standard deviation and $q^2$ in the liquid phase higher as one might expect from larger charge fluctuations
* but interestingly the mean scatters more around zero for the gas phase

### density histograms
$\beta C^{-1} = 0.04$
| | | | |
| ---- | --- | ---- | ---- |
| $(\kappa \sigma)^{-1}$ = 0.5 | $(\kappa \sigma)^{-1}$ = 0.7 | $(\kappa \sigma)^{-1}$ = 0.9 | $(\kappa \sigma)^{-1}$ = 1.1 |
|  |  | |  |

$(\kappa \sigma)^{-1} = 0.6$
| | | | |
| ---- | --- | ---- | ---- |
| $\beta C^{-1} = 0.07$ | $\beta C^{-1} = 0.05$ | $\beta C^{-1} = 0.03$| $\beta C^{-1} = 0.01$ |
|  |  | |  |

$(\kappa \sigma)^{-1} = 1.4$
| | | | |
| ---- | --- | ---- | ---- |
| $\beta C^{-1} = 0.15$ | $\beta C^{-1} = 0.1$ | $\beta C^{-1} = 0.08$| $\beta C^{-1} = 0.05$ |
|  |  | | |

* less density fluctuations when the phase separation is stronger;
* the peak of the fluid phase is always broader;
* I think this behaviour is more a feature of the Gibbs ensemble simulations than a coupling of charge and density fluctuations