# Free energy calculations of defects in 12-fold quasicrystals
###### tags: `Post-doc projects`
---
### Owners (the only one with the permission to edit the main text, others can comment)
Giovanni, Laura, Frank, Alptug
---
### TODO:
- [ ] Plug in the vibrational free energy to the lattice sim to see the equib. defect density. (There is a slight $\Delta F$ between shields and eggs.)
- [ ] Can we transform $\epsilon_D \rightarrow n_D$ rigorously?
- [ ] Write.
- [ ] Publish.
- [ ] Profit.
---
## Notation and definitions
##### Square shoulder potential
\begin{equation}
V_\text{SS}(\vec{r}_{ij}) =
\begin{cases}
\infty & \quad , \quad |\vec{r}_j-\vec{r}_i| \leq \sigma \\
\varepsilon & \quad , \quad \sigma < |\vec{r}_j-\vec{r}_i| \leq \delta \\
0 & \quad , \quad |\vec{r}_j-\vec{r}_i| > \delta
\end{cases}
\end{equation}
$\sigma$: core's diameter
$\delta$: potential well's diameter
##### Simulation units and parameters
$\sigma$: length
$\varepsilon$: energy
$\mu$: mass
$\tau = \sqrt{\frac{\varepsilon}{\mu\sigma^2}}$: time
$\delta = 1.4\sigma$
$k_BT = 0.1\varepsilon$
$F_M(N,V,T) \equiv$ free energy of a system of $N<M$ particles on a lattice of $M$ sites in a volume $V$ at temperature $T$
$N=$ # of occupied sites in the lattice
$M=$ total # of lattice sites (# particles + # vacancies)
## Free energy of a lattice with vacancies
Generally, the free energy of a crystalline lattice of $M$ sites with $M-N$ vacancies depends both on the vibrational properties and their mutual location. We write the free energy gain for the creation of a single vacancy at a specific lattice site $i$ as:
$$
-f^v_M(V,T) = F^v_M(N-1,V,T;i) - F_M(N,V,T)
$$
where I keep the sign convention as in Sec.~10.4 of Frenkel-Smit's book.
We consider the situation $N \sim M$ (non-interacting defects), for which the $f^v_M$ does not depend on $N$ (how many vacancies) or $i$ (their location), so that we can write:
$$
F^v_M(N,V,T;\{i\}) = M f_M(M,V,T) - (M-N) f^v_M(V,T)
$$
as the free energy of the system with $M-N$ vacancies at sites $\{i\}$, where $f_M$ is the defect-free free energy per particle.
In such an approximation, it is possible to partition the phase space in many different basins of equal hypervolume (same vibrational free energy) associated to many distinct configurations (configurational free energy).
Then we are allowed to write:
$$
F_M(N,V,T) = M f_M(M,V,T) - (M-N) f^v_M(V,T) - T S_M^\text{conf}(N)
$$
where $S_M^\text{conf}(N)=k_B^{-1}\log\Omega_{M,N}$ contains the counting of the all distinct configurations of (M-N) vacancies on a lattice of $M$ sites.
In the case of a 12-fold quasi-crystalline lattice further considerations are needed. The subtraction of a particle generates two half vacancies in the quasi-lattice in the form of eggs (e) or shields (s), that in principle should have different vibrational properties. Then, for the case $N=M-1$ there will be in principle three different possible values for the free energy of a vacancy: $f^v_\text{ss},f^v_\text{se},f^v_\text{ee}$, that in the non-interacting approximation between tiles should be possible to decompose into:
\begin{eqnarray}
f^v_\text{ss} &=& 2 f^{v/2}_s \\
f^v_\text{se} &=& f^{v/2}_s + f^{v/2}_e \\
f^v_\text{ee} &=& 2 f^{v/2}_e
\end{eqnarray}
where $f^{v/2}$ is the free energy associated to the creation of a defective egg or shield tile.
In the end we should be able to correctly write the total free energy for a 1-vacancy system as:
$$
e^{-\beta F_M(N)} = e^{-\beta M f_M(M)} \left[
e^{\frac{S^\text{conf}_\text{ss}(N)}{k_B}+\beta f^v_\text{ss}} +
e^{\frac{S^\text{conf}_\text{es}(N)}{k_B}+\beta f^v_\text{es}} +
e^{\frac{S^\text{conf}_\text{ee}(N)}{k_B}+\beta f^v_\text{ee}}\right] \equiv e^{-\beta M f_M(M)} Z^v(V,T)
$$
with the relative weight of the exponentials on the r.h.s. determining the appearance frequency of each defects configuration pair.
[Probably, according to the factorization properties of $e^{S^\text{conf}_{xx}(N)/k_B}$ it is possible to express it as a function of single eggs and shields? In that way it should be possible to analogously express the free energy of a system with $M-N>1$ vacancies, accounting for all the possible mixing of $2(M-N)$ eggs and shields, isn't it?]
Lastly, to be able to calculate the equilibrium fraction of defects, we need to calculate the Gibbs free energy:
\begin{eqnarray}
G_M(N,P,T) &=& \inf_V\left\{F_M(N,V,T)+P~V\right\} \\
&=& F_M(N,V_M(N),T)+P~V_M(N)
\end{eqnarray}
Then, we can calculate the Gibbs free energy for the creation of a vacancy at constant pressure and density as:
\begin{eqnarray}
g^v(P) &=& G_M(M-1,P)-G_{M-1}(M-1,P) =\\
& & \\
&=& F_M(M-1,V_M(M-1)) - F_{M-1}(M-1,V_{M-1}(M-1)) +\\
& & + P\left[ V_M(M-1) - V_{M-1}(M-1) \right] =\\
& & \\
&=& F_M(M-1,V_M(M-1)) - F_M(M-1,V_M(M)) +\\
& & + F_M(M-1,V_M(M)) - F_M(M,V_M(M)) +\\
& & + F_M(M,V_M(M)) - F_{M-1}(M-1,V_{M-1}(M-1)) +\\
& & + P\left[ V_M(M-1) - V_{M-1}(M-1) \right]
\end{eqnarray}
where $V_M(M)$ and $V_M(M-1)$ are the equilibrium volumes of the one-vacancy and defect-free system at same $N,P,T$; in the last equality I added and subtracted terms relative to a system in which a vacancy is created at constant volume in a lattice of $M$ sites, and to its defect-free analogous, following the same procedure as in Daan Frenkel's textbook.
In such a way we obtain:
\begin{eqnarray}
g^v(P) &=& P\left[ V_M(M) - V_M(M-1) \right] +\\
& & - k_BT\log Z^v(V,T) +\\
& & + f_M(\rho^*,T) +\\
& & + P\left[ V_M(M-1) - V_{M-1}(M-1) \right] \\
& & \\
&=& - k_BT\log Z^v(V,T) +\\
& & + f_M(\rho^*,T) + P/\rho^*
\end{eqnarray}
where $\mu^*=f_M(\rho^*,T) + P/\rho^* = G_M(M,P,T)/M$ is the chemical potential of the defect-free lattice. Finally, we can maybe write:
$$
e^{-\beta g^v(P,T)} =
e^{\frac{S^\text{conf}_\text{ss}(N)}{k_B}-\beta (\mu^*-f^v_\text{ss})} +
e^{\frac{S^\text{conf}_\text{es}(N)}{k_B}-\beta (\mu^*-f^v_\text{es})} +
e^{\frac{S^\text{conf}_\text{ee}(N)}{k_B}-\beta (\mu^*-f^v_\text{ee})}
$$
The same way we can probably split also the different $g^v_{xx}$:
\begin{eqnarray}
g^v_\text{ss} &=& \mu^* - f^v_\text{ss} &=& 2 g^{v/2}_s \\
g^v_\text{se} &=& \mu^* - f^v_\text{se} &=& g^{v/2}_s + g^{v/2}_e \\
g^v_\text{ee} &=& \mu^* - f^v_\text{ee} &=& 2 g^{v/2}_e
\end{eqnarray}
So that we can probably write in the general case:
$$
e^{-\beta G_M(N,P,T)} = e^{-\beta G_N(N,P,T)} \sum_{n_e,n_s}^{n_e+n_s=2(M-N)} \left[ e^{S^\text{conf}(M,n_s,n_e)/k_B - \beta (n_sg^{v/2}_s + n_eg^{v/2}_e)} \right]
$$
and probably:
$$
G_M(N,P,T) \approx G_N(N,P,T) + \overline{n}_s g^{v/2}_s + \overline{n}_e g^{v/2}_e - T\Delta S^\text{conf}(M,\overline{n}_s,\overline{n}_e)
$$
where $\Delta S^\text{conf}(M,n_s,n_e)$ is the configurational entropy of $M-N$ vacancies found in a lattice of $M$ sites as $n_s$ shields and $n_e$ eggs.
### Can we transform $\epsilon_D \rightarrow n_D$ rigorously?
Let's say that our configurational entropy in the open boundary system is:
\begin{equation}
S(N,\epsilon_D,\gamma) = \log\sum_{\mathcal{C}_{i,N}} e^{-n_D(i)\epsilon_D}
\end{equation}
where $n_D$ is the number of defects, $\epsilon_D$ is the energy cost divided by $k_BT$ per single defect, $\gamma$ the surface tension, and $\mathcal{C}_{i,N}$ is the $i$-th configuration of N vertices with $n_D(i)$ defects. (Maybe unit-wise a $k_B$ is missing?)
Can we say that since
\begin{equation}
\left\langle n_D\right\rangle = -\frac{\partial S(N,\epsilon_D,\gamma)}{\partial \epsilon_D}
\end{equation}
$n_D$ and $\epsilon_D$ are conjugated variables, then allowing to make a Legendre transform? To obtain:
\begin{equation}
S(N,n_D,\gamma) = S(N,\epsilon_D,\gamma) + n_D\epsilon_D
\end{equation}
This way $S(N,n_D)$ would be more or less of the same kind of $\Delta S^\text{conf}(M,\overline{n}_s,\overline{n}_e)$, that we can recast equivalently as $\Delta S^\text{conf}(N,n_d,n_s)$, where $n_D=2(M-N)=n_s+n_e$.
Then, we know that for same $\epsilon_D$ (that should be $\beta g^{v/2}$?) $n_s=1/3$. If we have the fluctuations of $n_s$ in Alptug's simulations, can we calculate the shift in the average value of $n_s$ without repeating the simulations? (like in the spirit of a sort of histogram reweighting?)
That way we would be only left with the minimization along $n_D$ to find the average defect concentration, isn't it?
### Calculation of the single-vacancy free energy cost
To calculate the free energy cost to produce a vacancy at constant volume at specific positions in a quasi-lattice, as a first step I calculate the free energy of a defect-free configuration of $M=N=2911$ sites, as previously done with Alptug, with the Frenkel-Ladd method.
I used the usual biasing potential:
$$
V_\text{Ein} = \lambda\sum_{i}^{N}\frac{(\vec{r}_i-\vec{R}_i)^2}{\sigma^2}
$$
and calculated the integral in $\lambda$ with the Gauss-Legendre quadrature method with $n=50$ nodes and $c=0.1$ (see the note on [patchy particles](https://hackmd.io/Bqxq7vZ6Q1SDurPZnG7-LA) for more details).
The final formula I used to obtain the free energy per particle is:
``` python //
loglmin = np.log( c )
loglmax = np.log( lambdamax + c )
integral = np.add.reduce( lambda_weights * 0.5 * (loglmax-loglmin) *
* (lambda_vals+c) * lambda_unit_nodes , axis=0 )
Vlambda = /* time serie of the values of the potential energy at lambdamax*/
avgV = np.mean( Vlambda , axis=0 )
deltaF_ideal = ( avgV/kT - np.log(
np.mean( np.exp( -(Vlambda-avgV)/kT ) , axis=0 ) ) ) / M
deltaF = deltaF_ideal - integral/kT +
+ 0.5*D*(1./M-1.)*np.log(np.pi*kT/lambdamax) +
+ np.log(M/volume)/M - 0.5*D*np.log(M)/M
```
I calculated the values of $F_M(M,V,T)$ for 4 values of density $\rho$, averaged over 10 different realizations.
>[name=Frank] For safety: have you checked whether you can integrate between these using the Equation of State?
>[name=Giovanni] Yes, the discrepancy is minimal (order of the error on the integrated values)
#### Free energy of a configuration with 1 vacancy
I calculated the free energy of a configuration by starting from the 10 configurations used in the previous section, by generating 2 rhombi defects and letting them evolve until they were far apart of $~L/2$ half simulation box. Finally, I popped out two rhombi vertices to suitably create shield and egg defects, resulting in a lattice of N=M-1=2910 sites.
I considered 4 cases: configurations with 2 shields, an egg and a shield, two eggs of same chirality and two eggs with opposite chirality.
To calculate the Free energy $F^v_M(M-1,V,T)$ I used the same formula of previous section, by changing $M$ with $M-1$.
###### Simulations.
I performed simulations with the constrained center of mass and biasing potential $V_\text{Ein}$ by further imposing the constraint that each particle should not go out from its Wigner-Seitz (WS) cell, and close to defects they should not invade the WS cell of every possible free rhombus vertex (in that case it would individuate a different configuration).
To do that I slightly changed the code I used to simulate the defect-free system. Here I highlight some changes in the initialization:
``` C++ //
lambda = /* V_Ein spring constant */ ;
cutoff = /* cutoff to check cells occupancy */ ;
reference_density = 1.0773486881311538 ; /* max QC12 SS density */
cutoff *= sqrt( reference_density / N * L^2 ) ;
// non - generation of the verlet list, used to map neighbors to check cell's escaping
verlet_list = new verlet<particle>( conf_ptr , parts_number, cutoff , 0.0 , "NO_DISPLACEMENT_VECTOR" ) ;
/* ... */
// initialization of the array of starting positions
starting_positions = new particle [N] ;
for( int i=0; i<conf_ptr->numberOfPart; i++ )
starting_positions[i] = particles[i]->unwrapped( box_sides ) ;
// This is a vacancies neighbour list for every particle:
neigh_ghost_verts = new particle** [N] ;
/* ... */
// This is a list of all the Nvac possible rhombi vertices within defects
empty_sites = new particle [Nvac] ;
/* Then I load the positions of the empty rhombi sites from file ... */
int vac_Nneighatoms = 0 , *neigh_array = NULL , neigh_array_length = 0 , neigh_atom = 0 ;
int *neigh_Nvac = new int [N] ;
for( int i=0; i<Nvac; i++ ) {
/* Here I find the particles close to the empty rhombus site */
vac_Nneighatoms = verlet_list->get_neighbour_indices( empty_sites+i , &neigh_array_length , &neigh_array ) ;
/* and for every neighbour I insert the empty site i in the vacancies neighbour list : */
for( int j=0; j<vac_Nneighatoms; j++ ) {
neigh_atom = neigh_array[j] ;
neigh_ghost_verts[ neigh_atom ][ neigh_Nvac[neigh_atom] ] = empty_sites+i ;
neigh_Nvac[ neigh_atom ]++ ;
neigh_ghost_verts[ neigh_atom ] = (particle **)realloc( neigh_ghost_verts[ neigh_atom ] , (neigh_Nvac[neigh_atom]+1)*sizeof(particle *) ) ;
neigh_ghost_verts[ neigh_atom ][ neigh_Nvac[neigh_atom] ] = NULL ;
}
}
```
In the following, instead, the body of the function computing the interactions between the particles after each MC move (all the rest of the algorithm is left unchanged, which is reasonably correct, after comparison with Alptug some months ago).
``` C++ //
idx = /* index of the particle that has been moved */ ;
particle upart , image_part ;
upart = particles[ idx ]->unwrapped( box_sidxes ) ;
energy = lambda_potential( &upart , &(equilibrium_positions[idx]) ) ;
// check if the particle is still in its own Wigner-Seitz cell
// CoM_displacement is the displacement of the center of mass since the beginning
upart -= ( equilibrium_positions[idx] + CoM_displacement ) ;
walker = verlet_list->atom[idx].neighbours ;
while( walker != NULL ) {
verlet_list->nearest_image( &(equilibrium_positions[walker->neigh]) , &(equilibrium_positions[idx]) , &image_part ) ;
image_part -= equilibrium_positions[idx] ;
if( ( (upart * image_part) / image_part.square_norm() ) > 0.5 ) return INFINITY ;
walker = walker->next ;
}
// check if the particle does not invade defects' Wigner-Seitz cells
neigh_ghost = neigh_ghost_verts[ idx ] ;
while( *neigh_ghost ) {
verlet_list->nearest_image( *neigh_ghost , &(equilibrium_positions[idx]) , &image_part ) ;
image_part -= equilibrium_positions[idx] ;
if( ( (upart * image_part) / image_part.square_norm() ) > 0.5 ) return INFINITY ;
neigh_ghost++ ;
}
```
##### CHECKS
To check that the results obtained are correct, I changed different things and verified that everything stays the same.
###### Different cutoff
In the first simulations I used a cutoff of $1.6$ lattice spacings, which should be enough to grasp all the neighbours determining the Wigner-Seitz cell. I also tried to put it to $2.6$ lattice spacings (with the actual value in length units depending on the density), and the simulations do not change.
###### Different WS cells check for vacancies
In another batch of simulations, I substituted the check on the WS cell with ghost rhombi vertices with the procedure suggested by Frank, based on the projection of the 4D WS cell around the half-vacancy.
I defined in the code the vectors:
$$
\vec{z}_i = \sqrt{\frac{\rho_0}{\rho}} \frac{ \vec{e}_i - \vec{e}_{(i-1)\%12} }{2}
$$
which are half the minor diagonal of rhombi at density $\rho$ ($\rho_0$ is the density at which $\vec{e}_i$ have unit norm). In the following the part that differ in the calculation with respect to the previous case:
``` C++ //
zeta_sqnorm = zeta[0].square_norm() ;
neigh_ghost = neigh_ghost_verts[idx] ;
upart += equilibrium_positions[idx] ;
while( *neigh_ghost ) {
verlet_list->nearest_image( &upart , *neigh_ghost , &image_part ) ;
image_part -= (*neigh_ghost) ;
flag = 1 ;
for( int i=0; i<12; i++ ) flag = flag & ( ( (image_part * zeta[i]) / zeta_sqnorm ) < 1.0 ) ;
if( flag ) return INFINITY ;
neigh_ghost++ ;
}
```
###### Assessment of the statistical error
The error bars on the calculation of $f^v_{xx}$ are quite large, in comparison with the differences between $f^{v/2}_s$ and $f^{v/2}_e$, then I tried to separately estimate the statistical variance on calculations over different trajectories of the same configuration. It comes out that the standard deviation on the free energy of shield-shield defects at the lowest density over 10 different configurations, is exactly the same as over 10 trajectories from the same configuration ($30.73833 \pm 0.00020$ versus $30.73827\pm 0.00020$).
###### Defects statistics at $k_BT=0.1\varepsilon$
I am running $NVT$ simulations at the four different densities to gather statistics on the defects counting, to check if the predictions from $f^v$ and configurational entropy will be consistent.
I tried with EDMD, but simulations were to slow, since the temperature is too low, then I moved to MonteCarlo with an additional jumping move by one of the twelve vectors $2\vec{z}_i$.
### Results for $f^v$ and $g^v$ at $\beta\epsilon = 10$
Finally, with these calculations, by applying the first formula for $f^v_{xx}$ I obtained these values, for the free energies of vacancies of different kind ($s-s$, $s-e_x$, $e_x-e_x$, $e_L-e_R$).
<figure>
**Figure** Free energy cost to make a vacancy $f^v_{xx}$ at constant $(V,T)$
</figure>
<figure>
**Figure** Free energy cost to make a vacancy $g^v_{xx}$ at constant $(\rho,P,T)$
</figure>
Things are quite messy at densities around $\rho\sim1.05$. I am checking that no mistakes are present in the simulations. In that case we would need more statistics (longer simulations? More replies?), or to use a larger system size.
As we can see in the figure below showing the free energy densities of the systems with vacancies with respect to the one with two shields as a reference, especially in the equilibrium region the differences are around or less than the error on the free energies ($\sim0.0002$), and this strongly affects the calculation of $f^v$, since:
\begin{equation}
f^v = f_M +(M-1)(f_M - f_{M-1})
\end{equation}
where $f_M$ and $f_{M-1}$ are the per-particle free energies of the defect free system and the 1-vacancy system, respectively.
In this case is easy to see that the error on $f^v$ scales as $\delta f^v\sim M\delta f$. This could not necessarily constitute a big problem, but the fundamental thing would be that we need to have as accurate as possible estimates of the various $f$, I guess$\ldots$
### Results for $g^{v/2}$ of defective tiles at different state points
| $\beta\epsilon$ | $\rho^\text{df}$ | $g^v_\text{egg}$ | error | $g^v_\text{shield}$ | error |
|--| ----------------- | ----- | ----- | ----- | ----- |
| $10$ | $1.0485$ | $16.11$ | $0.31$ | $15.65$ | $0.40$ |
| $5$ | $1.0186$ | $6.73$ | $0.33$ | $6.34$ | $0.41$ |
| $10/3$ | $0.9906$ | $4.00$ | $0.34$ | $3.65$ | $0.42$ |
How we can see, the statistical error is quite large on $g^{v/2}$, and this could depend on the error on $f^v$.
### Calculation of $f^v$ and $g^v$ at $\beta\epsilon = 5$
To stay within the equilibrium region of the quasicrystal, I calculated the free energy of the systems with and without vacancies at density $\rho=1.0485$ (the density refers to the defect-free system, vacancies are created by keeping V constant):
| | $\beta f(\rho,T)$ | error |
|--| ----------------- | ----- |
| defect free | $32.45728$ | $0.00016$ |
| eggs (same chirality) | $32.44347$ | $0.00020$ |
| eggs (opposite chirality) | $32.44376$ | $0.00030$ |
| egg and shield | $32.44360$ | $0.00018$ |
| shields | $32.44330$ | $0.00022$ |
and integrated over a linear path $(\beta\epsilon = 10, \rho = 1.0485)_i \longrightarrow (5,1.0186)_f$ of 21 points (considering initial and final point).
Along isochores we would have:
\begin{equation}
\beta_f f_f = \beta_i f_i + \int_{\beta_i}^{\beta_f}\left\langle U\right\rangle_\beta d\beta
\end{equation}
Since in our system the partition function is a function of $\beta\epsilon$, with a potential energy $U=\epsilon n_b$, where $n_b$ is the number of bonds ($|\vec{r}_j-\vec{r}_i|\leq \sigma+\delta$), we can rewrite:
\begin{equation}
\beta_f f_f = \beta_i f_i + \int_{(\beta\epsilon)_i}^{(\beta\epsilon)_f}\left\langle n_b\right\rangle_\beta d(\beta\epsilon)
\end{equation}
Then the complete relation I used over the considered integration path is:
\begin{eqnarray}
\beta f(\rho_f,(\beta\epsilon)_f) &=&
\beta f(\rho_i,(\beta\epsilon)_i) +
\int_{\rho_i}^{\rho_f} \frac{\beta P}{\rho^2}d\rho +
\int_{(\beta\epsilon)_i}^{(\beta\epsilon)_f}\frac{\beta\left\langle U\right\rangle_{\beta\epsilon}}{\beta\epsilon} d(\beta\epsilon) \\
&=& \beta f_i + \beta\Delta f_\rho + \beta\Delta f_T
\end{eqnarray}
where the quantities I directly get from the EDMD simulations are $\beta P$ and $\beta U$.
>[name=Giovanni] memo: have to check in the code how the energy is updated, since apparently some times it happens that it looses a bond (see the anomalously large error bars)$\ldots$
>[name=Giovanni] I checked, and at every collision at the shoulder's edge I update both the energy and the momentum exchanged for the calculation of P$\ldots$ Ther should be some lost collision, or I have to check if something weird appears in the system, but it is improbable$\ldots$
Both the contributions are considered negative, since the direction of the integration path.
##### Defect free system
For the defect free system I obtain a final vibrational free energy per particle of $\beta f(\rho=1.0186,\beta\epsilon=5) = 18.350714\pm0.000083$, which seems to be consistent with the value calculated directly at the same state point by thermodynamic integration of $\beta f(\rho=1.0186,\beta\epsilon=5) = 18.35077\pm0.00044$.
##### 1-vacancy system
For the system with 1 vacancy I had to add two further contributions, since I integrated it over the same path, but we have to consider that the free energy values that we need to calculate $g^v$ needs to refer to systems with and without vacancy at the same volume, then the starting and ending densities are $\rho_{0,i/f}=\frac{M-1}{M}\rho_{i/f}$. Then I added an initial compression and final dilatation step to the integration, at constant $\beta\epsilon$, to avoid redoing all the calculations:
\begin{equation}
\beta f_f = \beta f_i + \beta\Delta f_c + \beta\Delta f_\rho + \beta\Delta f_T + \beta\Delta f_e
\end{equation}
where on top of the previous contributions, we also have the free energy of the initial compression:
\begin{equation}
\beta\Delta f_c (\rho_{0,i}\rightarrow\rho_i) = \frac{1}{2}\left(\frac{\beta_i P_i}{\rho_i^2}+\frac{\beta_{0,i} P_{0,i}}{\rho_{0,i}^2}\right)(\rho_{i}-\rho_{0,i})
\end{equation}
which is positive, and the free energy difference due to the final expansion (negative contribution):
\begin{equation}
\beta\Delta f_e (\rho_f\rightarrow\rho_{0,f}) = \frac{1}{2}\left(\frac{\beta_{0,f} P_{0,f}}{\rho_{0,f}^2} + \frac{\beta_f P_f}{\rho_f^2}\right)(\rho_{0,f}-\rho_f)
\end{equation}
To calculate the contribution for the eggs I averaged the initial free energies over the two different chirality couplings.
**Warning:** Along the integration path the defects start to diffuse, then in these calculations we assume that the pressure of the system with one vacancy is insensitive to the species of defects (eggs or shields), since the integrand functions are averages over configurations containing both. Is that correct?
### Free energies at $\beta\epsilon=5$
| | Defect free system ($\rho_f=1.0186\sigma^{-2}$) | error |
|--| ----------------- | ----- |
| $\beta f_i$ | $32.45728$ | $0.00016$ |
| $\beta\Delta f_\rho$ | $-1.42576$ | $0.00005$ |
| $\beta\Delta f_T$ | $-12.68015$ | $0.00002$ |
| $\beta f_f$ | $18.35137$ | $0.00017$ |
| | Eggs ($\rho_f=1.01825\sigma^{-2}$) | error |
|--| ----------------- | ----- |
| $\beta f_i$ | $32.44361$ | $0.00014$ |
| $\beta\Delta f_c$ | $0.024888$ | $0.000002$ |
| $\beta\Delta f_\rho$ | $-1.43827$ | $0.00006$ |
| $\beta\Delta f_T$ | $-12.67420$ | $0.00002$ |
| $\beta\Delta f_e$ | $-0.012197$ | $0.000001$ |
| $\beta f_f$ | $18.34383$ | $0.00016$ |
| | Shields ($\rho_f=1.01825\sigma^{-2}$) | error |
|--| ----------------- | ----- |
| $\beta f_i$ | $32.44330$ | $0.00022$ |
| $\beta\Delta f_c$ | $0.024885$ | $0.000003$ |
| $\beta\Delta f_\rho$ | $-1.43820$ | $0.00006$ |
| $\beta\Delta f_T$ | $-12.67422$ | $0.00001$ |
| $\beta\Delta f_e$ | $-0.012197$ | $0.000001$ |
| $\beta f_f$ | $18.34357$ | $0.00023$ |
### Free energies at $\beta\epsilon=10/3$
| | Defect free system ($\rho_f=0.9906\sigma^{-2}$) | error |
|--| ----------------- | ----- |
| $\beta f_i$ | $32.45728$ | $0.00016$ |
| $\beta\Delta f_\rho^{(1)}$ | $-1.42576$ | $0.00005$ |
| $\beta\Delta f_T^{(1)}$ | $-12.68015$ | $0.00002$ |
| $\beta\Delta f_\rho^{(2)}$ | $-0.80847$ | $0.00003$ |
| $\beta\Delta f_T^{(2)}$ | $-4.23047$ | $0.00002$ |
| $\beta f_f$ | $13.31243$ | $0.00017$ |
| | Eggs ($\rho_f=0.9902597\sigma^{-2}$) | error |
|--| ----------------- | ----- |
| $\beta f_i$ | $32.44361$ | $0.00014$ |
| $\beta\Delta f_c$ | $0.024888$ | $0.000002$ |
| $\beta\Delta f_\rho^{(1)}$ | $-1.43827$ | $0.00006$ |
| $\beta\Delta f_T^{(1)}$ | $-12.67420$ | $0.00002$ |
| $\beta\Delta f_\rho^{(2)}$ | $-0.81232$ | $0.00004$ |
| $\beta\Delta f_T^{(2)}$ | $-4.22852$ | $0.00002$ |
| $\beta\Delta f_e$ | $-0.0083063$ | $0.0000007$ |
| $\beta f_f$ | $13.30688$ | $0.00016$ |
| | Shields ($\rho_f=0.9902597\sigma^{-2}$) | error |
|--| ----------------- | ----- |
| $\beta f_i$ | $32.44330$ | $0.00022$ |
| $\beta\Delta f_c$ | $0.024885$ | $0.000003$ |
| $\beta\Delta f_\rho^{(1)}$ | $-1.43820$ | $0.00006$ |
| $\beta\Delta f_T^{(1)}$ | $-12.67422$ | $0.00001$ |
| $\beta\Delta f_\rho^{(2)}$ | $-0.81230$ | $0.00004$ |
| $\beta\Delta f_T^{(2)}$ | $-4.22851$ | $0.00002$ |
| $\beta\Delta f_e$ | $-0.0083066$ | $0.0000007$ |
| $\beta f_f$ | $13.30665$ | $0.00023$ |
Sign in with Wallet
Connect another wallet