# Eline's Project
###### tags: `masters projects` ,`machine learning` ,`quasicrystal`

### Owners (the only one with the permission to edit the main text, others can comment)

Eline, Alptug, Laura

## Week 0 (15 & 16 September 2022)

:::spoiler To do (click to expand!)
::: info
- [x] try out hackmd
- [x] read paper square shoulder potential quasicrystals (https://aip.scitation.org/doi/pdf/10.1063/1.4977934 )
- [ ] write a project description
- [x] write code for "brute force" MC hard sphere in 2D, after that: cell lists!
- [x] read the part on cell lists in Frenkels book
:::

### Project Description (draft)
**Title: Free energy calculations for quasicrystals using machine learning**

**Project description:**
Quasicrystals are among the most fascinating phases in soft matter, since they have rotational symmetries that are impossible for normal periodic crystals. To be precise, quasicrystals have a long-range positional order without a translational periodicity. Moreover, a relatively simple model of spherical particles with a square shoulder potential, can be used to describe the complex system of quasicrystals[^first].

[^first]: https://aip.scitation.org/doi/pdf/10.1063/1.4977934

Despite the fact of their important role in soft matter, they are one of the least understood phases. Especially the free energy calculation remains to be an important problem. The free energy of quasicrystals consists of the configurational entropy, which is associated with the different realizations of the tiling, and the vibrational free energy, which is associated with the fluctuations of the particles around their lattice sites. It should be noted that these two are not independent of each other. For a specific quasicrystal realisation the vibrational free energy can be calculated using thermodynamic integration. However, the full free energy requires sampling the different quasicrystal realizations, according to their likelyhood. In principle, a Monte Carlo simulation could be created to sample these differnt configurations but it would be unrealistically slow as sampling configurations would require a thermodynamic intergration measurement to weight the configuration. Because of this, it is very difficult to efficiently compute the total free energy of quasicrystals. However, the new developments in fitting free energies using machine learning can provide us a solution to this problem: graph neural networks. Specifically, different configurations of quasicrystals and their vibrational free energies (using the thermodynamic integrations) can be generated to use for training the neural networks. Using this trained neural networks, the free energies of quasicrystals can be calculated more efficiently.

==Feel free to add comments, or change the description :smile:==

## Week 1 (19-23 September 2022)

:::spoiler To do (click to expand!)
::: info
- [x] implement cell lists in NVT hard spheres code
- [ ] catch up on doubly linked lists
- [ ] make the cell lists
- [ ] make a "move particle function"
- [x] check this function with the "brute force method"
- [x] instead of hard spheres, use square shoulder potential!
- [ ] adminstrative stuff (complex systems/extra supervisors)

:::

This week, I worked mainly on making the cell lists in an ordinary hard sphere NVT-ensemble (2D), using doubly linked lists.

:::spoiler Some code snippets:

void make_cell_lists(double r_cutoff)
{
    /*this function makes cell lists, using doubly linked lists */

    ncell = floor(box[0]/r_cutoff); //number of cells needed
    ncell_2 = ncell * ncell; //it is in 2d
    assert(ncell < ncell_max);

    //the box is always a square, so:
    size_cell = box[0] / floor(box[0]/r_cutoff);
    
    for (int x =0 ; x<ncell ; x++)
    {
        for (int y =0; y<ncell; y++)
        {

        head_of_chain[x][y] = NULL; //the cell lists are empty at first
        
        }    
    }
    
    /*loop over the particles to put them into the lookup table and the list */
    for (int i=0; i<n_particles; i++)
    {
        int cell_number_x =floor(r_particles[i][0]/size_cell);
        int cell_number_y = floor(r_particles[i][1]/size_cell);
        struct Node* old_hoc = head_of_chain[cell_number_x][cell_number_y]; //old head of chain


        /* 1. allocate node */
        struct Node* new_node = (struct Node*)malloc(sizeof(struct Node));
    
        /* 2. put in the data  */
        for (int d =0; d<NDIM; d++)
        {
            new_node->r[d] = r_particles[i][d];
        }

        /* 3. Make next of new node as head and previous as NULL*/
            new_node->next = (old_hoc);
            new_node->prev = NULL;
        
        /* 4. change prev of head node to new node */
        if ((old_hoc) != NULL)
            (old_hoc)->prev = new_node;
    
        /* 5. move the head to point to the new node */
        head_of_chain[cell_number_x][cell_number_y] = new_node;

        look_up_table[i] = new_node;
    }

}

int move_particle_cell_lists(void)
{
    double rand_particle[NDIM];
    random_particle = n_particles * dsfmt_genrand(); //the random particle
    struct Node *pntr = look_up_table[random_particle]; //points to the random particle

    for (int d =0 ; d<NDIM; d++) {rand_particle[d] = pntr -> r[d];} 

    //which cell?    
    int cell_number_x =floor(rand_particle[0]/size_cell);
    int cell_number_y =floor(rand_particle[1]/size_cell);


    for (int d =0; d<NDIM; d++)
    {
        new_pos[d]= rand_particle[d] + delta * (2.0 * dsfmt_genrand() - 1.0);
        //periodic boundaries :
        if(new_pos[d]>box[d]) {
            new_pos[d]-=box[d];
        }
        if(new_pos[d]<0.0) {
            new_pos[d]+=box[d];
        }

    }

    int new_cell_number_x = floor(new_pos[0]/size_cell);
    int new_cell_number_y = floor(new_pos[1]/size_cell);

    for(int x_cell =-1 ; x_cell<2 ; x_cell++)
    {
        for (int y_cell = -1; y_cell<2; y_cell++)
        {
            int neighbour_x = new_cell_number_x + x_cell;
            int neighbour_y = new_cell_number_y + y_cell;

            //periodic boundaries:
            if (neighbour_x == -1) {neighbour_x += ncell;}
            if (neighbour_y == -1) {neighbour_y += ncell;}
            if (neighbour_x == ncell) {neighbour_x = 0;}
            if (neighbour_y == ncell){neighbour_y = 0;}
            
            //check for overlap with all the particles in that cell:
            struct Node *current_node = head_of_chain[neighbour_x][neighbour_y];

            while (current_node != NULL)
            {
                if (current_node != pntr)
                {    
                    double dist = 0.0;
                    for(int d = 0; d < NDIM; ++d)
                    {
                        double min_d = new_pos[d] - current_node ->r[d];
                        /*periodic boundaries */
                        if(min_d>0.5*box[d]) { min_d -= box[d];}
                        if(min_d<-0.5*box[d]) {min_d += box[d];}
                        dist += min_d * min_d;
                    }
                if (dist <= diameter_2) {return 0;} //it cannot move
                }
                current_node = current_node -> next;


            }
        }
    }
    /*move the particle! */
    for (int d =0; d<NDIM; d++) {pntr->r[d] = new_pos[d];}

    if ( new_cell_number_x != cell_number_x || new_cell_number_y!= cell_number_y)
    { 
    //that means we need to delete the node in the old chain, and add the node into the new chain....
    /* step 1: delete the node */
 
    /* If node to be deleted is head node */
    if (head_of_chain[cell_number_x][cell_number_y] == pntr)
        head_of_chain[cell_number_x][cell_number_y] = pntr->next;
 
    /* Change next only if node to be deleted is NOT the last node */
    if (pntr->next != NULL)
        pntr->next->prev = pntr->prev;
 
    /* Change prev only if node to be deleted is NOT the first node */
    if (pntr->prev != NULL)
        pntr->prev->next = pntr->next;
 
    /* step 2: add the node to the other cell */
   struct Node* old_hoc = head_of_chain[new_cell_number_x][new_cell_number_y];

    pntr->next = (old_hoc);
    pntr->prev = NULL;
        
    /*change prev of head node to new node */
     if ((old_hoc) != NULL)
         (old_hoc)->prev = pntr;
    
    /* move the head to point to the new node */
    head_of_chain[new_cell_number_x][new_cell_number_y] = pntr;

    }
    return 1;
    }

:::

I have tried the cell list method with 1225 particles, on a square 2D lattice, using 33X33 cells, and compared it with the "brute force". According to my calculations it was approximatly 20 - 30 times as fast (of course - when using more particles, this will improve)!

## Week 2 (26-02 September/Oktober 2022)

:::spoiler To do (click to expand!)
::: info
- [x] fill in the form for thesis
- [x] Structure factor calculation (and a bit of reading)
- [x] implement square shoulder
- [x] improve cell lists code a bit

:::
Some of the plots of the structure factors:
Square lattice:


Hexogonal lattice:

Moreover, I did some small improvements in my code to speed up a little (calculating the neighbours beforehand, make the code work for rectangular lattices)

## Week 3 (3-9 Oktober 2022)

:::spoiler To do (click to expand!)
::: info
- [x] check if the code works using an initial configurations of quasicrystals
- [ ] Einstein/thermodynamic integration

:::

I tried my code with an initial condition of a 12-fold quasicrystal (1 with 836 particles, and one with 2911 particles), and I let it run for 100000 MC cycles. After all these cycles, the system is still a quasicrystal (see Structure factors:)
With 836 particles:


With 2911 particles:

And a GIF of the system:

# Einstein Integration mathematical derivation of the free energy difference
I have worked out the theory in the paper, and I just want to keep this somewhere (I closely follow the derivation in the paper, using the same arguments)

So, we first take the difference between the constrained and the unconstrained crystal:

$$
\beta F - \beta F^{cm} = - \ln(Q/Q^{cm})
$$

Where $Q$ is the partition function of an unconstrained crystal, and $Q^{cm}$ is the partition function of the constrained crystal (center of mass constraint). $Q$ is given by:

$$
Q = Z P
$$
Where $Z$ is the configurational contribution, and $P$ the kinetic contribution. We can calculate the kinetic contribution (in 3 dimensions):
$$
P = \int d^{3N}\mathbf{p} \prod^N_{i=1} e^{- \left({\beta} / {2m_i}\right) p_i^2} = \\
\int d^Np_xd^Np_yd^Np_z \prod^N_{i=1} e^{- \left({\beta} / {2m_i}\right) (p_{x,i}^2 + p_{y,i}^2+ p_{z,i}^2)} = \\
\left(\int d^{N}p \prod^N_{i=1} e^{- \left({\beta} / {2m_i}\right) p_i^2} \right)^3 = \\
\prod^N_{i=1} \left( \frac{2\pi m_i }{\beta} \right) ^\frac{3}{2}
$$

(for the first step: I just wrote out the three dimensions, and since they are independent, you can just put one dimension to the power of3 , and then you are left with an ordinary gaussian integral!)

$Q^{cm}$ is computed in the paper (it is basically the same integration, but using a constraint):
$$
Q^{cm} = Z^{cm} \left(\frac{\beta}{2\pi M}\right)^{3/2} \prod^N_{i=1} \left( \frac{2\pi m_i }{\beta} \right) ^\frac{3}{2}
$$

So lets get back to:

$$
\beta F = - \ln(Q/Q^{cm}) + \beta F^{cm} = \\
\ln(Z^{cm} /Z) - \ln(\frac{2\pi M}{\beta})^{3/2} +\beta F^{cm}
$$

Then, we use that $Z^{cm}/Z = N/V = \rho$ (argumentation in the paper)

so we get:
$$
\beta F = \ln(\rho) - \ln(\frac{2\pi M}{\beta})^{3/2} +\beta F^{cm}
$$

Then, we take a closer look at the free energy difference between the constained crystal and the reference system (Frenkel-Ladd integration):

$$
\beta F^{cm} = \beta F^{cm}_{Ein} - \beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm} = \\
-\ln(Q^{cm}_{ein})- \beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm}
$$

Then, we take a closer look at $\ln(Q^{cm}_{ein})$

$$
Q^{cm}_{ein} = Z^{cm}_{ein}P^{cm}_{ein} = \\
\left(\frac{\alpha \beta}{2 \pi \sum_i \mu_i^2}\right)^{3/2} \left(\frac{ \beta}{2 \pi M}\right)^{3/2} Z_{ein}P_{ein}
$$
(Again, derivation is almost the same as for $Q^{cm}$, and is explained in the paper).

Lets put everything together:

$$
\beta F = \ln(\rho) - \ln(\frac{2\pi M}{\beta})^{3/2} +\beta F^{cm} =\\
\ln(\rho) - \ln(\frac{2\pi M}{\beta})^{3/2} -\ln(Q^{cm}_{ein})- \beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm} = \\
\ln(\rho) - \ln(\frac{2\pi M}{\beta})^{3/2} -\ln(\left(\frac{\alpha \beta}{2 \pi \sum_i \mu_i^2}\right)^{3/2} \left(\frac{ \beta}{2 \pi M}\right)^{3/2} Z_{ein}P_{ein})- \beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm} =\\
\ln(\rho) - \ln(\frac{2\pi M}{\beta})^{3/2} -\ln(\left(\frac{\alpha \beta}{2 \pi \sum_i \mu_i^2}\right)^{3/2}) -\ln( \left(\frac{ \beta}{2 \pi M}\right)^{3/2}) - \ln(Z_{ein}P_{ein})- \beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm} =\\
$$
Until here, I just filled in the expressions we already obtained and split the log
$$
\ln(\rho) - \ln(\frac{2\pi M}{\beta})^{3/2} + \ln( \left(\frac{2 \pi M}{ \beta}\right)^{3/2}) -\ln(\left(\frac{\alpha \beta}{2 \pi \sum_i \mu_i^2}\right)^{3/2}) - \ln(Z_{ein}P_{ein})- \beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm} =\\
\ln(\rho) -\ln(\left(\frac{\alpha \beta}{2 \pi \sum_i \mu_i^2}\right)^{3/2}) - \ln(Z_{ein}P_{ein})- \beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm}
$$
Two terms cancel! Now we want to look at $Z_{ein}$ and $P_{ein}$

$$
Z_{ein} = \int d^{3N}\mathbf{r} \prod^N_{i=1} e^{- \left({\alpha\beta} / {2}\right) r_i^2} = \\
\left(\int d^{N}r \prod^N_{i=1} e^{- \left({\alpha\beta} / {2}\right) r_i^2} \right)^3 = \\
\left( \frac{2\pi }{\alpha\beta} \right) ^\frac{3N}{2}
$$
(note that again, do the same trick as in $P$)
and again:

$$
P_{ein} = \int d^{3N}\mathbf{p} \prod^N_{i=1} e^{- \left({\beta} / {2m_i}\right) p_i^2} = \\
\\
\left(\int d^{N}p \prod^N_{i=1} e^{- \left({\beta} / {2m_i}\right) p_i^2} \right)^3 = \\
\prod^N_{i=1} \left( \frac{2\pi m_i }{\beta} \right) ^\frac{3}{2} =\\
\left( \frac{2\pi m }{\beta} \right) ^\frac{3N}{2}
$$
since, in this case, $m_i = m$.

Lets use this:

$$
\beta F =\ln(\rho) -\ln(\left(\frac{\alpha \beta}{2 \pi \sum_i \mu_i^2}\right)^{3/2}) - \ln(Z_{ein}P_{ein})- \beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm} =\\
\ln(\rho) -\ln(\left(\frac{\alpha \beta}{2 \pi \sum_i \mu_i^2}\right)^{3/2}) - \ln(\left( \frac{2\pi }{\alpha\beta} \right) ^\frac{3N}{2}) - \ln(\left( \frac{2\pi m }{\beta} \right) ^\frac{3N}{2})- \beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm}
$$

Moreover, $\sum_i \mu_i^2 = 1/N$ (in our case, since $m_i = m$)
$$
\beta F =\ln(\rho) -\ln(\left(\frac{\alpha \beta N}{2 \pi }\right)^{3/2}) - \ln(\left( \frac{2\pi }{\alpha\beta} \right) ^\frac{3N}{2}) - \ln(\left( \frac{2\pi m }{\beta} \right) ^\frac{3N}{2})- \beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm}
$$
This is equation (16) from the paper :)

To get to equation (17), we need to calculate $F_{ex} = F - F_{id}$, where:

$$
\beta F_{id} = -\ln(V^N(2 \pi m/\beta)^{3N/2}) + \ln(N!) \approx \\
-N\ln(V) - \ln((2 \pi m/\beta)^{3N/2}) +N\ln(N) - N + (\ln(2\pi N))/2 =\\
N\ln(\rho) - \ln((2 \pi m/\beta)^{3N/2}) - N + (\ln(2\pi N))/2
$$

so, we obtain:
$$
\beta F_{ex} =\ln(\rho) -\ln(\left(\frac{\alpha \beta N}{2 \pi }\right)^{3/2}) - \ln(\left( \frac{2\pi }{\alpha\beta} \right) ^\frac{3N}{2}) - \ln(\left( \frac{2\pi m }{\beta} \right) ^\frac{3N}{2})- \\
\beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm} - N\ln(\rho) + \ln(\left( \frac{2\pi m }{\beta} \right) ^\frac{3N}{2}) + N - (\ln(2\pi N))/2
$$
We can see that two terms cancel, so:

$$
\beta F_{ex} =\ln(\rho) -\ln(\left(\frac{\alpha \beta N}{2 \pi }\right)^{3/2}) - \ln(\left( \frac{2\pi }{\alpha\beta} \right) ^\frac{3N}{2}) - \\
\beta \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm} - N\ln(\rho) + N - (\ln(2\pi N))/2
$$

So:
$$
\frac{\beta F_{ex}}{N} = \frac{\ln(\rho)}{N} -\frac{3}{2N}\ln\left(\frac{\alpha \beta N}{2 \pi }\right) - \frac{3}{2}\ln\left( \frac{2\pi }{\alpha\beta} \right) - \\
\frac{\beta}{N} \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm} - \ln(\rho) + 1 - \frac{ \ln(2\pi N)}{2N}
$$

which we can rewrite such that we got exactly (17) from the paper:

$$
\frac{\beta F_{ex}}{N} = \frac{\ln(\rho)}{N} -\frac{3}{2N}\ln\left(\frac{\alpha \beta }{2 \pi }\right) - \frac{3}{2N}\ln(N) - \frac{3}{2}\ln\left( \frac{2\pi }{\alpha\beta} \right) - \\
\frac{\beta}{N} \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm} - \ln(\rho) + 1 - \frac{ \ln(2\pi )}{2N} - \frac{1}{2N}\ln(N)
$$
so:

$$
\frac{\beta F_{ex}}{N} = \frac{\ln(\rho)}{N} -\frac{3}{2N}\ln\left(\frac{\alpha \beta }{2 \pi }\right) - \frac{4}{2N}\ln(N) - \frac{3}{2}\ln\left( \frac{2\pi }{\alpha\beta} \right) - \\
\frac{\beta}{N} \int_0^1d\lambda \langle \Delta U \rangle _{\lambda}^{cm} - \ln(\rho) + 1 - \frac{ \ln(2\pi )}{2N}
$$

# Free energy calculations

From statistiscal physics we know we can compute the Free energy $F(N,V,T)$ of a system by taking the logarithm of the partition function $Z(N,V,T)$:

$$
F(N,V,T) = -k_b T \ln Z(N,V,T) = -k_bT \ln \left( \frac{\int d \mathbf{r}^N \exp [-\beta U (\mathbf{r}^N)] }{N! \Lambda^{3N}} \right).
$$
Where $k_b$ is the Boltzmann constant, $T$ the temperature, $\beta = \frac{1}{k_bT} - DEFINE U ETC? . However, this partition function cannot be computed directly for most systems, because the number of degrees of freedom is too large. However, we can use thermodynamic integration to compute the differences bbetween the free energy. For this, we can use the potential energy defined as $U(\lambda) = U_A + \lambda (U_B - U_A)$, so when $\lambda = 0$, we are in system A, and when $\lambda =1$ we are in system B. Now we can write for the difference between the free energies:

$$
\Delta F(A \rightarrow B) = \int_0^1 \frac{\partial F(\lambda)}{\partial \lambda}d \lambda
$$

And we can write:

$$
\left(\frac{\partial F(\lambda)}{\partial \lambda} \right)_{N,V,T} = - \frac{1}{\beta} \frac{\partial }{\partial \lambda} \ln Z(N,V,T,\lambda) \\
= -\frac{1}{\beta Z(N,V,T,\lambda)} \frac{\partial Z(N,V,T,\lambda) }{\partial \lambda} \\
= \frac {\int d \mathbf{r}^N \exp [-\beta U (\mathbf{r}^N)] } {}
$$
# Einstein Integration of hard spheres