A Multimer Embedding Approach for Molecular Crystals up to Harmonic Vibrational Properties

Accurate calculations of molecular crystals would be crucial for drug design and crystal engineering. However, periodic high-level density functional calculations using hybrid functionals are often prohibitively expensive for relevant systems. Therefore, these expensive periodic calculations are sometimes circumvented by the usage of quantum-mechanical embedding methods, in which for instance the periodic calculation is only performed at a lower-cost level and then monomer energies and dimer interactions are replaced by those of the higher-level method. Herein, we extend upon such a multimer embedding approach to enable energy corrections for trimer interactions and the calculation of harmonic vibrational properties up to the dimer level. We show that trimer interactions are crucial for accurately approximating lattice energies and might also be needed for a better description of lattice constants. Finally, vibrational properties are already very well captured at the dimer level.


I. INTRODUCTION
The capability to accurately but still efficiently model molecular crystals would be invaluable for crystal engineering 1 and drug design 2 . However, the individual molecules within molecular crystals are only weakly held together by non-covalent interactions and for many molecules different crystal-packing arrangements are possible. Such different polymorphs can have very similar lattice energies 3 , which often differ by only a few kJ/mol. Therefore, it is vital to accurately capture the subtle interplay of intermolecular interactions. Furthermore, single-point energy calculations or simple lattice relaxations are also often insufficient since many properties of molecular crystals can highly depend on temperature and pressure 4 . For instance, the actual relative stability of polymorphs can often not be determined by static lattice energies alone, but rather free energies have to be considered [5][6][7] , which means that computationally expensive vibrational free energies have to be calculated as well. In addition, it might often also be necessary to explicitly account for the thermal expansion of the crystal 4, [8][9][10][11] , which requires at the minimum several harmonic phonon calculations.
Given the periodic nature and the often large unitcell sizes of relevant molecular crystals, highly accurate wave-function methods cannot routinely be used and we have to rely on more approximate methods. Currently, the main workhorse for high-level calculations of molecular crystals is periodic density functional theory (DFT). One important way of assessing the quality of computational methods under realistic conditions are the regular crystal structure prediction blind tests organized by the Cambridge Crystallographic Data Centre [12][13][14][15][16][17] , with the a) johannes.hoja@uni-graz.at b) adrian daniel.boese@uni-graz.at current blind test just being completed in June this year. Therein, van-der-Waals dispersion inclusive density functional approximations are often very successfully used in the final steps of such crystal structure prediction procedures.
Among the density functional approximations, hybrid functionals are generally more accurate than functionals solely based on the generalized gradient approximation (GGA) -but also significantly more expensive. Several examples indicate that for instance the hybrid PBE0+MBD approach can yield more accurate results for molecular crystals and improve upon the PBE+MBD description at the GGA level 7,18-23 . However, fully converged periodic hybrid calculations can easily become prohibitively expensive for practically relevant systems or in cases when a huge number of calculations are required, for instance during crystal structure predictions. In addition to the immense increase in CPU time compared to GGAs, also the memory requirements can often not be satisfied for large unit-cell sizes.
One possible solution to this problem is the usage of quantum-mechanical embedding schemes, which approximate the periodic hybrid calculation with less expensive calculations. These approaches typically make use of a molecular many-body expansion, i.e., they involve monomers, molecular dimers, trimers, etc., and are often referred to as fragment methods 24-28 . Any periodic high-level method can now be approximated by either an additive or a subtractive scheme. In the additive case 29-31 , monomer energies, dimer interactions, trimer interactions, etc. are summed up, eventually converging to the periodic result and thereby completely circumventing an explicitly periodic calculation. In contrast, a subtractive scheme [32][33][34][35] involves an explicit periodic calculation at a lower-level method followed by replacing monomer energies, dimer interactions, trimer interactions, etc. with the values from the high-level method.
Recently, several such quantum-mechanical embedding approaches have been developed for molecular crystals utilizing even up to MP2 or CCSD(T) as high-level method 31,33 . Also, Chen and Xu 36 have introduced a different fragmentation scheme involving only parts of a molecule. Instead of including the most expensive methods, embedding is also utilized to approximate for instance GGA calculations by GGA fragments and periodic density-functional tight binding calculations in order to enable very large calculations 35 .
One of us has introduced a subtractive embedding scheme for approximating hybrid density functionals 34 , which consists of a periodic GGA calculation and a monomer and dimer correction utilizing the hybrid functional. This methodology was implemented in Ref. 34 for energies and lattice relaxations and the following hybrid:GGA combinations were tested: PBE0:PBE+D3, PBE0:PBE+MBD, and B3LYP:BLYP+D3.
In this paper, we present an extension and a new opensource implementation of the above-mentioned embedding approach in order to speed up or enable hybrid calculations for larger molecular crystals. Specifically, we extend the energy calculation up to trimers and enable harmonic phonon calculations, which can now be performed utilizing up to dimers. We test the performance of the resulting multimer embedding approach by embedding PBE0+MBD 37-39 multimers into periodic PBE+MBD 40 calculations and comparing with explicit periodic PBE0+MBD results utilizing the X23 8,18,41 benchmark set.

A. Energy
Within our subtractive multimer embedding scheme, the periodic high-level energy E high per is approximated according to The first term E low per refers to the fully periodic calculation utilizing a computationally more efficient lower-level method, while the following terms then replace certain energies with results from a high-level method and are explained in detail below.
Note that we prefer to utilize here the term multimer embedding in order to indicate that our fragments are complete molecules and in order to refrain from the term many-body expansion for specifying multimers, since this may cause confusion with the used many-body dispersion (MBD) method, wherein a body refers to an atom. We suggest to label results with this scheme as for instance MEX(PBE0+MBD:PBE+MBD), where the first term in paranthesis is the high-level method, the second term the low-level method, and the X refers to the utilized multimer order, i.e., 1 if only monomers are included, 2 when up to dimers are considered, and 3 when up to trimers are used. Dimers and trimers are considered if their shortest intermolecular distance is smaller than a defined multimer cutoff distance. Currently, for trimers to be considered, the distance between all pairs of molecules within the trimer must have an intermolecular distance smaller than the multimer cutoff value. Since we will only utilize PBE0+MBD embedded into PBE+MBD in this paper, we will omit this information in parenthesis further on but rather use this way here to specify the used cutoff distance.
The input unit cell is unwrapped and then a supercell of sufficient size created based on the given multimer cutoff distance, from which then all necessary multimers are extracted. The sums in eq. 1 always run over the monomers within the created supercell and the respective n amounts to the number of monomers of a given multimer that belong to the central unit cell. For the first sum, the value of n i is simply 1 in case monomer i belongs to the central unit cell and 0 otherwise. For the second term, n ij is 2 if both monomer i and monomer j belong to the central unit cell, 1 if only monomer i or j belongs to the central unit cell, and 0 otherwise. The same is true for the third sum containing trimers. Any ∆E term always refers to the difference between the high-level and the low-level quantity The term E int ij corresponds to the dimer interaction energy where ij is the dimer consisting of monomers i and j. Similarly, the term E int ijk refers to the trimer interaction energy where ijk refers to the trimer consisting of monomers i, j, and k. The lattice energy E high latt ist then calculated according to where Z is the number of molecules within the unit cell and E high mon is the high-level energy of one monomer in its most stable gas-phase conformation.

B. Forces
For any given atom a belonging to monomer i within the central unit cell (i = j), the approximated high-level force acting on it is given by with f low per (a) being the force on atom a from the lowlevel periodic calculation; ∆f i (a) refers to the difference in the forces on atom a between the high-level and lowlevel calculation of monomer i, and f int ij (a) is the dimer interaction force at atom a C. Stress The stress tensor σ is a second-order symmetrical tensor consisting of nine components six of which are unique. Expressions for the derivative of lattice parameters and the stress tensor for fragment/multimer methods have for instance been published by Nanda and Beran 32 and Loboda et al. 34 . Here, we approximate the stress tensor of our high-level method with where σ low pq is the obtained stress tensor from the low-level periodic calculations and V is the unit cell volume. The first summation in both terms sums up over all monomers and dimers, respectively. The second summation runs over all atoms a in the respective multimer, with atom a being within monomer i. The meaning of n is the same as above and accounts for monomers within the central unit cell. The indices p and q range from 1 to 3 and indicate weather the x, y, or z component of the position r and force f of atom a is to be used.

D. Harmonic Vibrational Properties
For the calculation of vibrational properties, we utilize phonopy 42 to create the necessary finite displacements within sufficiently large supercells. The corresponding high-level forces for each atom in these supercells are approximated as described in section B, treating each unit cell as if it were the central unit cell. So for every displaced atom, all corresponding displaced multimers are calculated and the periodic force constants corrected accordingly. With the approximated high-level force sets, the vibrational properties are then calculated within phonopy.

E. Computational Details
All shown multimer embedding calculations were performed by using our new open-source code MEmbed 43 (version 0.2.0). All electronic structure calculations were performed by utilizing FHI-aims [44][45][46][47][48][49] (version 210716 2), which enables the calculation of isolated and periodic systems on an equal footing due to the employed numeric atom-centered basis functions. MEmbed makes use of the Atomic Simulation Environment 50 (ASE) for several of its functionalities, including the interface with FHIaims.
Throughout, we used either the PBE 40 or the PBE0 37 density functional approximation together with the many-body dispersion 38,39 (MBD) method for proper accounting of van-der-Waals dispersion interactions. The calculations were performed by utilizing the light species default settings within FHI-aims for integration grids and basis functions (version 2020) in order to allow for the calculation of the canonical periodic PBE0+MBD method as reference. Performing all shown calculations in a fully periodic way with PBE0+MBD using converged tight settings would not be feasible due to the massive amount of required CPU time and memory, especially for the supercells needed for the phonon calculations.
Note that these embedding calculations differ from the corresponding results by Loboda et al. 34 . Therein, periodic calculations were performed by utilizing plane waves in VASP, while the multimers were calculated with a large TZVPPD 51,52 Gaussian basis set within TURBO-MOLE. Furthermore, all dispersion corrections were calculated using the PBE0 range-separation parameter for MBD 34 , effectively using the multimer scheme only for the DFT part. Herein, we perform all calculations on a completely equal footing in terms of software and basis sets, and we use for all PBE and PBE0 calculations MBD with the respective default value for the range-separation parameter.
In all FHI-aims single-point calculations the total energy, the forces, the charge density, and the sum of eigenvalues were converged to 10 −6 eV, 10 −4 eV/Å, 10 −5 electrons/Å 3 , and 10 −5 eV, respectively. The k-grids for the periodic DFT calculations were set to always satisfy nx > 18Å, with x being the cell length in the respective direction and n the number of k-points in that direction. For the MBD energy a tighter k-grid satisfying nx > 25 A was used. For phonon calculations appropriate supercells with displacements of 0.005Å were created so that the length in every direction exceeds 12Å and the q-grid used for the evaluation of the vibrational free energy was set to satisfy nx > 50Å. This resulted in all acoustic modes being smaller than 0.8 cm −1 in magnitude at the gamma point. All lattice relaxations were performed until a convergence of 0.005 eV/Å was reached using ASE with the BFGS algorithm and the ExpCellFilter 53 class, as well as the FixSymmetry class to maintain symmetry.

A. Embedding for Lattice Energies
First, we assess the performance of our multimer embedding approach for lattice energies. In order to eliminate any geometry effects within this comparison, we compare here now only single point calculations carried out with different approaches on top of the PBE0+MBD/light-optimized structures. The systems of the X23 benchmark set are listed together with their number of molecules within the unit cell (Z) and their number of atoms per molecule (n) in Table I.
In order to illustrate the complexity of the evaluated embedding approaches, Table I also contains for each multimer cutoff distance the number of identified unique dimers and trimers. Note that, whenever possible, symmetry constraints are utilized in order to limit the number of unique multimers to calculate. For instance, without symmetry a 3Å cutoff ME3 calculation for the ammonia crystal would already require 42 dimer and 76 trimer calculations. In contrast, with symmetry only 2 dimers and 7 trimers are necessary in this case. Fig. 1 shows relative lattice energies of our low-level method (PBE+MBD/light) and a selection of embedding approaches for all X23 systems compared with the used high-level method (PBE0+MBD/light) and Table II lists several error statistics for all considered approaches. Since the lattice energies have a negative sign, positive values in the plot and the reported mean errors (ME) indicate a smaller interaction magnitude than the reference.
It can be seen that PBE+MBD yields in almost all cases an underbinding with a ME of 4.9 kJ/mol and a mean relative error (MRE) of -5.2 % compared to PBE0+MBD; only for ammonia and cyanamide PBE+MBD yields a larger lattice energy in magnitude. Hence, the mean absolute error (MAE) is slightly larger than the ME, amounting to 5.3 kJ/mol or 5.9 % in terms of the mean absolute relative error (MARE). In the worst case (anthracene), the largest observed difference between PBE+MBD and PBE0+MBD is 13.1 kJ/mol or 11.9 %.
The monomer embedding (ME1) shows the same qualitative trend as PBE+MBD. However, the inclusion of only monomers at the PBE0+MBD level already reduces the ME to 3.3 kJ/mol as well as the MAE to 3.7 kJ/mol. When dimers are included, all lattice energies are now smaller in magnitude than the reference and at a multimer cutoff of 3Å the MAE is reduced to 3.3 kJ/mol. When utilizing larger cutoff values the errors remain quite similar but actually increase slightly.
As soon as trimers are included, we see a significant improvement in terms of all errors and now some lattice energies become again larger in magnitude than the reference value. By only considering dimers and trimers up to 3Å the ME and MAE can be reduced to 1.0 and 1.1 kJ/mol, respectively. When moving on to 4Å, we ob-serve the best agreement with PBE0+MBD with the ME, MAE, and the maximal error being only -0.2, 0.4, and 1.2 kJ/mol, respectively. At this level the high-level lattice energy is extremely well approximated when considering frozen geometries. Increasing the multimer cutoff further leads then to a small increase in all errors.
It can be clearly seen that the lattice energies often follow the common convergence behavior of a perturbational approach. So in order to reduce again the errors at large multimer cutoff, tetramer energies would probably need to be included. However, there seems to be a quite beneficial error cancellation at small supercell cutoffs, so that 3 or 4Å seem to be sufficient here.

B. Embedding for Forces and Stress
After having evaluated the energies at fixed geometries, we now study how well the forces and the stress tensor of the high level can be approximated. First, we have calculated the atomic forces and the stress tensors at PBE+MBD-optimized structures, so that the components are non-zero for our PBE0+MBD reference method, while they are virtually zero for our utilized lowlevel method. The resulting errors are shown in Table III.
It can be seen that this resulted in an average absolute force component difference between PBE0+MBD and PBE+MBD of 0.2 eV/Å. The monomer embedding already significantly reduces the MAE to 0.023 eV/Å and with dimer embedding we can reduce it with 0.006 eV/Å almost to our optimization convergence criterion. This implies that the atomic force contributions can be very well approximated by dimer embedding.
In terms of the non-zero components of the stress tensor, we are able to reproduce them with a MAE of 7 × 10 −4 eV/Å in the case of ME2(3Å) and ME2(4Å). Note that the dimers have a small effect on stress tensor components within the currently used approximation. While these errors might look very promising, they are still large enough to lead to quite different cell volumes for systems with flat potential energy surfaces like molecular crystals, as we will discuss below.
Next, we have performed lattice relaxations utilizing monomer and dimer embedding up to 5Å and the errors of the resulting cell volumes and the corresponding lattice energies compared to the optimized PBE0+MBD values are listed in Table IV. PBE+MBD/light overestimates in all cases the PBE0+MBD/light cell volume with a mean error of 3.8 %. All shown embedding approaches also always overestimate here the cell volume compared to PBE0+MBD. In the case of monomer embedding the mean error can be reduced to 2.9 % and the best dimer embedding (4Å) leads to a MAE of 2.3 %. Hence, we see that the accuracy of the stress tensor is not sufficient at the dimer level to satisfactory approximate the PBE0+MBD lattice constants.
Table IV also shows the errors of the resulting lattice   II. Errors of the calculated lattice energies of the X23 set compared to PBE0+MBD. All calculations were done with light settings on top of the PBE0+MBD-optimized structures. The mean error (ME), the mean absolute error (MAE), and the maximal error (MAX) are given in kJ/mol while the mean relative error (MRE), the mean absolute relative error (MARE), and the maximal relative error (RMAX) are given in %. energies at the optimized cells when compared to those of canonical PBE0+MBD lattice relaxations. It can be seen that these errors are very similar to the results for the frozen PBE0+MBD structures. This indicates that the potential energy surfaces of these systems are indeed very flat and that a very high accuracy of the stress tensor is likely needed to actually reproduce the PBE0+MBD cell volumes. Interestingly, the PBE+MBD mean errors are now much smaller than for the frozen PBE0+MBD structures, which is due to the fact that now at the actual PBE+MBD equilibrium structures the resulting lattice energies have increased in magnitude. So far we have only evaluated the performance of the embedding by comparison with the canonical PBE0+MBD calculations. To put the results a bit more into perspective, we briefly also mention the performance when comparing to reference values derived from experimental sublimation enthalpies. In order to directly evalu- These results suggest that forces and stress for trimers are necessary to further improve the description of the lattice constants. Therefore, we perform for cubic ammonia a preliminary test how including trimers would modify the lattice parameters (see Fig. 2). It can be seen that our dimer embedding significantly overestimates the lattice constant and hence the cell volume compared to PBE0+MBD as well as to PBE+MBD. When trimers are now included, the lattice constant gets significantly smaller. At the ME3(3Å) level there is an excellent agreement with the PBE0+MBD value; moving to larger multimer cutoffs seems to lead then to a slight underestimation of the lattice constant. This illustrates that utilizing trimers within lattice relaxation could indeed significantly improve our description of lattice constants.  Next, we evaluate the performance of our embedding approach for vibrations/phonons using the X23 benchmark set. Therefore, we have calculated with several methods the gamma-point frequencies for the respective optimized lattice constants and for an internally relaxed structure with the PBE0+MBD lattice constants in order to determine if there are significant changes due to the different lattice constants (see Table V). It can be seen that there are significant differences between PBE+MBD and PBE0+MBD in terms of vibrational frequencies with a MAE of almost 50 cm −1 . These large differences originate mainly from the higher-frequency modes; when comparing only the first 300 cm −1 the resulting MAEs for the PBE0+MBD cell and for the optimized structure are only 3.3 and 6.6 cm −1 , respectively.

Method
Monomer embedding already significantly improves the internal vibrational modes and hence produces frequencies with an accuracy of about 4 cm −1 when evaluated for the entire frequency range. When dimers are considered, we can reach an overall MAE of about 3 cm −1 and a ME of about 1 cm −1 for the optimized structures already at a multimer cutoff of 3Å. For the overall statistic the difference in geometries does not seem to have a large impact, but that is simply due to the large number of internal modes. When considering only frequencies up to 300 cm −1 we observe slightly larger errors with ME2(3Å) and ME2(4Å) having a MAE of 6 and 4 cm −1 at the optimized structures, respectively.
Finally, we discuss the accuracy for vibrational free energies. In Table V we compare vibrational free energies evaluated at a converged q-grid at a temperature of 300 K and normalized per molecule with the respective PBE0+MBD result. It can be seen that PBE+MBD deviates by about 10 kJ/mol from PBE0+MBD and that the monomer embedding already provides a very accurate approximation of the PBE0+MBD vibrational free energy. For the optimized structures the error is about 1 kJ/mol and for the PBE0+MBD lattice constants the MAE amounts to only 0.3 kJ/mol. Including dimers slightly decreases this error further, leading to a MAE of only 0.8 and 0.2 kJ/mol for the optimized structure and the PBE0+MBD lattice constants, respectively. The vibrational free energies are already tightly converged at a multimer cutoff of 3Å, increasing it does not lead to further changes in terms of the MAE.

IV. CONCLUSION
We have introduced here an extension and new implementation of a subtractive multimer embedding scheme based on the work of Loboda et al. 34 , in order to enable larger calculations for molecular crystals utilizing hybrid functionals. Due to the fact that only up to trimers have to be calculated with the high-level method (hybrid functional), this approach is very memory efficient and can also be easily parallelized over the multimer calculations. Herein, we approximated the periodic PBE0+MBD results by performing periodic calculations using only the more efficient PBE+MBD approach and then introducing the effects of PBE0+MBD by improving monomer energies, dimer interaction energies, and trimer interaction energies.
The performance of the shown approach was then evaluated by directly comparing the embedding results with canonical periodic hybrid calculations. In order to accurately approximate lattice energies, the newly incorporated trimer energies are crucial. For lattice relaxations, the dimer embedding yields an error of about 2 % in terms of the cell volume. A numerical test on the ammonia crystal illustrated that trimer interactions can significantly improve the description of the cell volume. Hence, the next crucial step towards improving this multimer embedding methodology is the development of gradients and stress including trimer interactions.
Finally, we have also newly introduced here the calculation of vibrational properties utilizing multimer embedding. Herein, we are able to approximate gamma-point vibrational/phonon frequencies with an accuracy of a few wave numbers using dimer embedding. This enables an accurate approximation of room temperature vibrational free energies within 1 kJ/mol on average when normalized per molecule in the unit cell.

ACKNOWLEDGMENTS
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 890300. The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC).