 Welcome back. In the first part of this lecture we have seen how in a QM scheme a molecular system is usually partitioned that is a small region treated at a quantum level, computational expensive, and the rest of the system with a classical resolution, computational less demanding. Obviously special attention has to be paid when coupling the quantum and the classical regions. In particular in the Hamiltonian of this fictitious QMM system fictitious because the real system is not separated into parts with two different physics. The potential energies contain three types of interactions, interactions between particles in the quantum region that is between electrons and nuclei, interactions between atoms in the classical region, and interactions between quantum particle and classical particles. According to how one takes into account these last interactions, the different QMM schemes can be grouped in subtractive and additive schemes. In the subtractive scheme the potential energy of the QMM system is obtained in three steps by performing three independent calculations. That is first the energy of the total system is evaluated at the mm level, for example at force field level. Then the energy of just the isolated quantum subsystem is calculated at quantum level and added to the previous result. Finally the mm energy of the quantum subsystem is subtractive. This to avoid to countwise the interaction within the quantum subsystem. The main advantage of this QMM coupling scheme is that in the code no communication is required between the quantum mechanical and the classical routines. This makes the implementation relatively straightforward. However, there are also several drawbacks. A major disadvantage is that in the first and in the third step above a force field is required for the quantum subsystem, which may not always be available. In addition, this force field needs to be sufficiently flexible to describe, for example, the effect of chemical changes when a reaction occurs. A further drawback of this class of methods is the total absence of the polarization of electron density due to the classical environment. This shortcoming can be particularly problematic when modelling, for example, biological charge transfer processes, since these are usually mediated by the protein environment. Therefore, for a realistic description of such cases, a more consistent treatment of the interactions between the electrons and their surrounding environment is needed. This can be obtained through an additive scheme. In such schemes, there are no calculations at different resolutions taking place separately, like in the previous scheme. But the potential energy for the whole QMM system is a sum of three contributions, quantum energy terms, classical energy terms and QMM classical coupling terms. In contrast to the subtractive schemes, the interactions between particles in the quantum region and the classical atom in the MMM region are treated here explicitly in the last term of this expression. This additive approach is the coupling scheme implemented in the QMM interface of the Cp2K code that will be described in the tutorial. The quantum energy term usually comes from the DFD-Konsham and Miltonian in such an additive scheme. The classical energy terms comes almost always from a classical force field and its choice can be limited by the availability of force field implemented in the code that you want to use it. While the third term describing the interaction between quantum and classical regions is usually decomposed into parts, the bonded and the non-bonded parts. The bonded part is the part that describes the covalent interactions between quantum atoms, that is the atoms in the quantum regions, and classical atoms in the MMM region. Therefore, this bonded part is a present only when the boundary between the quantum and the classical regions cuts a covalent bond connecting a quantum atom to a classical atom, as shown in the picture here if you look at the yellow fragments. Moreover, in this case, care has to be taken when evaluating the quantum wave function, that is when you solve the Schrodinger equation for the quantum part. In fact, a straightforward cut through the QMM boundary will create one or more unpaid electrons in the quantum subsystem. In reality, these electrons are paired in bonding orbitals with electrons belonging to the atoms in the MMM sides. However, now those electrons do not exist in the MMM region due to our artificial partitioning and the lower level of resolution that we decided to use in this region. In literature, a number of approaches have been proposed to remedy the artefact that originates from such open valences, many of which are available in C2K. For example, we can saturate the dangling valences with a monovalent capping atom, usually a hydrogen atom at an appropriate position along the bond vector. Alternatively, it is possible to use the link atom pseudo-potential approach, which consists in introducing in the QMM system a description of the classical atom at the border bounded to the quantum atom through a special pseudo-potential with the required valence charge. It can be shown that this method requires constraining the bond distance appropriately. Another approach is to place a set of hybrid orbitals on the boundary atom between the quantum and the MMM fragments. One of these orbitals is included in the QM region for the self-consistent optimization that we have seen in the first part of this lecture to find the wave function of the quantum system. While the other orbitals are treated as auxiliary orbitals that do not participate in the QMM optimization, but they provide an effective electric field that contributes to the external potential felt by the dynamical electrons. Let's come back to the QMM energy coupling term. Beyond the bond part, there is also a non-bonded part that describes the interactions between the QM and the MMM regions not connected through a covalent bond. This part is usually formed by a steric term that is a term that takes into account the van der Waals interactions and the pure electrostatic term. In most of the QMM approaches, including in CP2K, the interactions in EQMM-bonded and EQMM-steric are handled at force field level. This means that, for example, in the bonded term, you can find chemical bonds between QM and MMM atoms modeled by harmonic potential, as well as the angles defined by one quantum atom and two MMM atoms, but also torsions involving at most two quantum atoms, which are commonly modeled by a periodic potential function. Similarly, the steric term EQMM-steric is usually described by a Lennard-Jones potential, and the parameters come from the force field used in the MMM part. The rest of the non-bonded interaction, that is the interaction between quantum and classical atoms separated by three or more atoms in the topology, are included in the electrostatic term. This is really an interaction between classical partial charges in the MMM region and the quantum charge density. But how is this term calculated in practice? Really, several approaches have been proposed for this term. In the so-called mechanical embedded, the simplest electrostatic coupling scheme, the electronic wave function is computed for an isolated quantum subsystem. In other words, the classical environment cannot induce polarization on the electron density in the quantum region. Within this approach, different ways to get EQMM electrostatic have been proposed. The simplest and very rough approximation consists to completely negate it, so in this case EQMM electrostatic equals to zero. Alternatively, one can assign classical charges to the atoms in the quantum region and then evaluate EQMM electrostatic in the usual way that is as a sum of per vice-coulomb terms between these quantum atoms and the classical atoms. And for the assignment of classical charges to quantum atoms, people have used either a fixed set of charges for the EQMM atoms, for example those given by the force field, or they have computed the charges somehow by, for example, the least squares fitting procedure to optimally reproduce some quantities such as the electrostatic potential at the surface of the quantum subsystem, and then recompute the charges ideally even at each integration step. An improvement of this mechanical embedding is to include polarization effects of the quantum region due to the presence of the classical atom in the EQMM region. In fact, in this so-called electrostatic embedding scheme the electrostatic interactions between the quantum and the EQMM subsystem are handled during the computation of the electronic wave function. That is, the EQMM electrostatic, which I remind you has to be headed in the quantum Hamiltonian, depends on both the classical charges and the quantum charge density. Increasing further the level of sophistication implies to include in the model also the polarizability of the classical atoms. In this electrostatic coupling approach, called polarization embedding scheme, both regions, quantum and classical, can mutually polarize each other. Although this last embedding offers the most realistic electrostatic coupling between the quantum and the classical regions, polarizable forces filled for biomolecular simulations are not so effective yet. Therefore, despite the progress in the development of such force fields, QMM studies with polarizable mm regions are so far not so popular. The Cp2k QMM interface adopts the electrostatic embedding strategy that now I will describe in more details. By construction, the QMM electrostatic energy term that in the electrostatic embedding one should add to the quantum Hamiltonian is the red expression in this slide. This means that the electrons that define rho see the classical atoms as special nuclei with non-integer and possibly even negative charges QI. As it is, this expression immediately arises two problems originating from the peculiar short range and the long range behavior of this term. That is the so-called electron spillout and the very large computational cost required to evaluate this term. Let's start with the issue at the short range. A problem that may arise when using standard partial charges to describe the charge distribution in the mm subsystem is the risk of overpolarization near the boundary. Note that now there are no covalent bonds in between quantum and classical regions. However, also in this case, the point charges on the mm side may attract or repel the quantum electrons too strongly, which could lead to electron density spilling out into the mm region as sketching the picture. In reality, this would not be possible due to the Pauli repartion of the electrons of the atoms in the mm regions. That, however, now are missing in a QM scheme. This phenomenon of electronic spillout at the boundary can become serious if a large flexible basis set or worst a non-local basis set is used in the QM calculation. As happened, for example, in CP2K, because with such basis sets the electrons are fully free to delocalize. This artifact can be avoided by modeling classical charges as smeared out charges instead of the traditional point charges, where the smearing function can be a Gaussian distribution or another suitable function centered on the mm atom. The smearing function used in CP2K to the same is this one. For the ones of you that are very interested in the theoretical details, this is the exact potential energy function generated by a Gaussian charge distribution. In contrast to the point charged model, with the smearing function, the coulomb interaction between the electrons and the smeared charge distribution does not diverge if the electrons approach the mm atoms. As you can see in the plot sketched here. Let's now move to the second issues in the calculation of EQM electrostatic in the electrostatic embedding. The one at long range. In fact, due to the coulomb long range behavior, the computational cost to compute the electrostatic EQM energy term, EQM electrostatic, is surprisingly large. Let's try to evaluate it. Typically, this term is calculated by collocating on the nodes of a grid in the quantum region the contribution coming from the mm potential. Therefore, the number of operation that a direct evaluation of this quantity requires can be estimated as the number of grid points times the number of mm classical atoms. Typical grid points number are of order of millions, while classical atoms are order of 10,000 or even more in system of biological interest. Therefore, it is evident that in a real system a brute force computation of this term is impractical. Many QMM codes face the problem by adopting strategies such as hierarchical methods or multiple techniques, which are based on the observation that contribution to EQM electrostatic of the classical atoms far and far away from the quantum region can be replaced by less expensive expressions. These methods are very effective, but require a fine tuning of parameters to yield optimal performance, which leads also to a loss of accuracy that makes error control sometimes difficult. Instead, Cp2k solves this problem by adopting a different approach based on the representation of the electrostatic potential of the classical atoms that is the electrostatic charge times this meaning function as a sum of function with different widths, sum of gaussians of different widths derived from the so-called gaussian expansion of the electrostatic potential. This in combination with the use of a real space multigrid technique. Even if the details of this approach are a bit beyond what I can introduce here, let's try to provide a few more elements to help you understanding this a bit better. Multigrid approach means that we use different grids, finer and coarser to represent on a computer continuous function in real space, like for example the function in the GIP expression. In the right hand side of the GIP expression, you can see NG gaussian function and another function called R loop. Now, the gaussian function in the reciprocal space have an intrinsic cut-off. Now remember what we have seen in the first part of this lecture. Each grid in the real space correspond to a cut-off in the reciprocal space. Of course, the computation on the finest grid is the most accurate one. If we map a gaussian function on the first grid whose reciprocal cut-off is equal to or bigger than the intrinsic cut-off of that gaussian function, then every gaussian function will be represented on the same number of grid points irrespective of its width. In practice, a grid with 25 points per side is usually sufficient for an optimal gaussian representation. But if each gaussian can be represented on the same number of grid points, this means that the gaussian can be considered a compact domain function, that is, it is zero beyond a certain distance. The smitting function on the left side and side of the GIP expansion instead cannot be correctly approximated as a function with a compact domain. Therefore, in practice, when calculating the EQMM electrostatic, the GIP decomposition allows converting the problem of mapping a non-compact smitting function, left-hand side term, on a fine grid, into the mapping of NG-approximately compact function on grids lower or at least equal to the fine grid, plus an additional non-compact function, RLU, which, however, is very smooth and therefore it can be mapped on the coarsest available grid. Why do we do that? Because the sum of the contribution of all the grids suitably interpolated is approximately equal to the function mapped analytically on the finest grid. That is the result that we want to calculate in practice. But the total cost, in this case, is one or two orders of magnitude smaller because the total number of operation for the grid points evaluation results to be much smaller. The algorithm to evaluate the EQMM electrostatic potential on the finest grid can be outlined as follows. First, each EQMM atom is represented as continuous Gaussian-like charge distribution via a-smitting function. Then, the electrostatic potential generating from this charge is fitted through a GIP expansion. Then, every Gaussian function of the GIP expansion is mapped on the first grid level whose reciprocal cut-off is equal to or bigger than the cut-off of that particular Gaussian function. This way, we can limit each Gaussian to only a finite domain without loss of accuracy, as mentioned before, so that only EQMM atoms embedded in the EQMM box or close to it will contribute to the finest grid as sketch in the picture. Finally, in the last step of the algorithm, the contribution at each grid level is interpolated starting from the co-orchest grid level up to the finest level because the EQMM electrostatic potential we are looking for, that is the one with the higher accuracy, is the one on the finest grid. Unlike other approaches, the lack of tuning parameters here makes this multi-grid implementation a totally free parameter scheme without any significant loss of accuracy. Consequently, very stable simulation can be obtained with also optimal energy conservation properties. Okay, we have arrived at the end of this brief introduction to EQMM approaches. If you have a question or doubts on this part, I'll be waiting for you in the Q&A session in a few minutes.