Multimer Embedding Approach for Molecular Crystals up to Harmonic Vibrational Properties

Accurate calculations of molecular crystals are crucial for drug design and crystal engineering. However, periodic high-level density functional calculations using hybrid functionals are often prohibitively expensive for the relevant systems. These expensive periodic calculations can be circumvented by the usage of 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 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 evaluate this approach for the X23 benchmark set of molecular crystals by approximating a periodic hybrid density functional (PBE0+MBD) by embedding multimers into less expensive calculations using a generalized-gradient approximation functional (PBE+MBD). We show that trimer interactions are crucial for accurately approximating lattice energies within 1 kJ/mol and might also be needed for further improvement of lattice constants and hence cell volumes. Finally, the vibrational properties are already very well captured at the monomer and dimer level, making it possible to approximate vibrational free energies at room temperature within 1 kJ/mol.


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 harmonic phonon calculations for several different unit-cell volumes.
Given the periodic nature and the often large unit-cell sizes of practically 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 a) johannes.hoja@uni-graz.atb) adrian daniel.boese@uni-graz.atcomputational 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 current blind test just having completed in 2022.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][19][20][21][22][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 memory requirements can often not be satisfied for large unit-cell sizes.
One possible solution to this problem is the usage of embedding schemes, which approximate the periodic hybrid DFT calculation with less expensive calculations.Such embedding approaches typically make use of a molecular many-body expansion, i.e., they involve monomers, molecular dimers, trimers, etc., and are often also referred to as fragment methods [24][25][26][27][28][29][30][31][32][33] .Any periodic high-level method can be approximated by either an additive or a subtractive scheme.In the additive case [34][35][36] , 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 [37][38][39][40][41][42] 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 embedding approaches have been developed for molecular crystals utilizing even up to MP2 or CCSD(T) as high-level method in a subtractive scheme [37][38][39][40]43,44 . Also,Chen and Xu 45 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 42 .
One of us has introduced a subtractive embedding scheme for approximating hybrid density functionals 41 , which consists of a periodic GGA calculation and a monomer and dimer correction utilizing the hybrid functional.This methodology was implemented in Ref. 41  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 [46][47][48] multimers into periodic PBE+MBD 49 calculations and comparing with explicit periodic PBE0+MBD results utilizing the X23 8,18,19,50 benchmark set, which has been extensively used to test and develop methods for molecular crystals 45,[51][52][53][54][55][56][57][58][59][60][61][62] .

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 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 label embedding results MEX(PBE0+MBD:PBE+MBD), where the first term in parenthesis 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 the information in parenthesis but rather use this way to specify the used cutoff distance.
From the input unit cell a supercell of sufficient size is 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 is 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, 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 (with i ̸ = j)

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 or the stress tensor for fragment/multimer methods have for instance been published by Nanda and Beran 39 and by Loboda et al. 41 .
Here, we approximate the stress tensor components of our high-level method with where σ low pq are the obtained stress tensor components 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.The meaning of n is the same as above and accounts for the number of 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 63 to create the necessary finite displacements within sufficiently large supercells.The corresponding high-level atomic forces are approximated as described in section B, but this time applying the correction to the whole phonopy supercell.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 64,65 (version 0.2.0).All electronic structure calculations were performed by utilizing FHI-aims [66][67][68][69][70][71] (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 72 (ASE) for several of its functionalities, including an interface with FHI-aims.
Throughout, we used either the PBE 49 or the PBE0 46 density functional approximation together with the many-body dispersion 47,48 (MBD) method (rsSCS version) for proper accounting of van-der-Waals dispersion interactions.Most 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 for the multimer embedding approach.Performing all 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.However, in order to validate the multimer embedding approach also for tight settings, we have additionally performed single-point energy calculations and lattice relaxations with tight settings (version 2020 defaults plus one additional auxiliary g function to improve the resolution of identity approximation 71 ) but lattice relaxations of the canonical PBE0+MBD/tight method were restricted to a subset of 8 small structures from X23.
Note that these embedding calculations differ from the corresponding results by Loboda et al. 41 .Therein, periodic calculations were performed using plane waves and pseudopotentials in VASP, while the isolated multimers were calculated with a TZVPPD 73,74 Gaussian basis set within TURBOMOLE.This lead to some inconsistency since in this case the low-level description of the periodic system and the multimers is not completely identical.Furthermore, all dispersion corrections were calculated using the PBE0 range-separation parameter for MBD 41 , effectively using the multimer scheme only for the DFT part.Herein, by using FHI-aims, we are able to perform all calculations (periodic and isolated multimers) on a completely equal footing in terms of software and basis sets since we utilize all-electron calculations with numeric atom-centered basis functions.In addition, 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 kgrids 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 MBD energies and forces a tighter k-grid satisfying nx > 25 Å was used.All phonon calculations were performed utilizing finite displacements of 0.005 Å and appropriate large supercells were created so that the length in every direction exceeds 12 Å.The qgrid 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 75 class, as well as the FixSymmetry class to maintain symmetry.For the calculation of lattice energies, we utilized for the isolated monomers the structures provided in Ref. 18 as a starting point and re-optimized them with the respective methods using a convergence criteria of 0.001 eV/ Å.

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 only single point calculations carried out with different approaches on top of the PBE0+MBD/lightoptimized 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 S1 in the Supporting Information (SI).Among the X23 systems n varies between 3 atoms (carbon dioxide) and 26 atoms (adamantane), while Z varies between 2 and 8, yielding unit cells containing between 12 and 72 atoms.
In order to illustrate the complexity of the evaluated embedding approaches, Table S1 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.For a succinic acid crystal, which involves slightly larger and less symmetric molecules with 14 atoms, symmetry can also significantly reduce the costs of a ME3 calculation with 3 Å cutoff from 19 dimers and 22 trimers to only 5 dimers and 4 trimers.We start by first discussing lattice energies obtained with light species default settings since these settings can be used to evaluate the performance of our multimer embedding approach all the way up to harmonic vibrational properties due to the fact that corresponding canonical periodic PBE0+MBD calculations are still computationally feasible.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 I 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.The corresponding lattice energies of all systems are listed in Tables S2 and S3 in the SI.All structures and detailed calculation results including all energies of all isolated multimers are further available in a Zenodo repository 76 as ASE database files.
It can be seen that PBE+MBD leads for almost all systems to 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.2 kJ/mol or 5.9 % in terms of the mean absolute relative error (MARE).In the worst case (anthracene), the largest observed absolute difference between PBE+MBD and PBE0+MBD is 13.3 kJ/mol.
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, which accounts for about 30 % of the difference between PBE0+MBD and PBE+MBD.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 observe 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 fixed geometries.Increasing the multimer cutoff further leads to a small increase in all errors.
In terms of convergence of the lattice energies with the used multimer order, a part of the systems always improve with increasing multimer order, while others basically follow a damped oscillation.The latter can be seen for instance for ammonia.There, the addition of dimers leads to a quite substantial underbinding, resulting in larger errors as at the ME1 level.Adding trimers leads to a smaller correction in the opposite direction, which significantly reduces the error.In general, utilizing trimers with a larger multimer cutoff leads in this case for X23 to a small overbinding.So in order to further reduce the errors at large multimer cutoffs, tetramer energies would probably need to be included.However, there seems to be a quite beneficial error cancellation at small multimer cutoffs, so that 3 or 4 Å are sufficient here.
After having established the convergence behavior of our multimer embedding approach with light species default settings, we discuss now lattice energies obtained with converged tight species default settings, which were calculated on top of PBE0+MBD/light-optimized structures.The obtained individual lattice energies are listed in Table S4 in the SI and the corresponding statistical errors w.r.t.periodic PBE0+MBD are given in Table S5 in the SI.For tight settings the difference between PBE0+MBD and PBE+MBD is for X23 smaller than for the light settings with a MAE of 3.0 kJ/mol.Monomer embedding (ME1) reduces the MAE to 2.2 kJ/mol.Adding now dimer corrections actually increases the MAE to for instance 2.8 kJ/mol at a cutoff of 3 Å.As soon as trimer interactions are corrected, the errors drop significantly -just like in the case of light settings.Here, we even reach already at the ME3(3 Å) level a MAE below 1 kJ/mol (0.7 kJ/mol) and at the ME3(4 Å) level the convergence is very similar to the one with light settings, with now ME, MAE, and MAX values of only -0.5, 0.5, and 1.4 kJ/mol, respectively.

B. Embedding for Forces and Stress
After having evaluated the energies at fixed geometries, we now study forces and stress tensors.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 II.
It can be seen that this resulted in an average absolute force component difference between PBE0+MBD and PBE+MBD of about 0.2 eV/ Å.The monomer em- bedding already significantly reduces the MAE by a factor of ten to 0.023 eV/ Å and with dimer embedding, it is further reduced by another factor of four to 0.006 eV/ Å, which is close to our optimization convergence criterion.This implies that the atomic force contributions are already 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 Å with light settings.The errors of the resulting cell volumes and the corresponding lattice energies compared to the optimized PBE0+MBD/light values are listed in Table III and the individual volumes can be found in Table S6 in the SI.
PBE+MBD/light overestimates the X23 cell volumes compared to PBE0+MBD/light with a mean relative error (MRE) of 3.8 %.All shown embedding approaches also always overestimate the cell volume compared to PBE0+MBD.In case of monomer embedding, the mean absolute relative error (MARE) can be reduced to 2.9 % and the best dimer embedding (4 Å) leads to a MARE of 2.3 %.Hence, the accuracy of the stress tensor is not yet sufficient at the dimer level to accurately approximate the PBE0+MBD lattice constants.
Table III also shows the errors of the resulting lattice energies at the optimized cells when compared to those of canonical PBE0+MBD lattice relaxations.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 probably 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.
In order to illustrate the impact of trimer interactions on lattice constants and unit-cell volumes, we have calculated as example single-point energies of cubic ammonia at varying lattice constants (see Fig. 2).Our dimer embedding significantly overestimates the lattice constant and hence the cell volume compared to PBE0+MBD as well as to PBE+MBD for the ammonia crystal.In fact, ammonia is for all cutoffs at the ME2 level the system with the worst agreement with the periodic PBE0+MBD results in terms of the cell volume.However, we note that the optimal unit-cell volumes for ME2(3 Å) and ME2(4 Å) obtained via a Murnaghan equation-of-state 77 fit from the single-point energies in Fig. 2 respectively agree within 0.3 % and 0.1 % with the corresponding embedding optimizations, which further validates the used stress-tensor approximation.
When trimers are 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 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.Since the number of trimers to be considered can be quite large and is increasing heavily with increasing cutoff, an efficient inclusion of trimer interactions for the calculation of forces and stress tensors requires an indepth study of which trimers are important and which could be omitted.Therefore, we will discuss explicit lattice relaxations with trimer interactions in a follow-up publication investigating also the performance of different dimer/trimer cutoff combinations and exploring other ways to reduce the number of considered trimers.
Next, we have also performed lattice relaxations utilizing converged tight species default settings.Given the massive computation time of the canonical periodic PBE0+MBD/tight calculations, we only compare a subset of X23 containing 8 structures, for which the reference calculations were still computationally feasible.The resulting individual unit-cell volumes are given in Table S7 and the errors compared to PBE0+MBD/tight in Table S8 in the SI.In this case, PBE+MBD/tight has a MARE of 2.5 % for the cell volume and both monomer embedding as well as dimer embedding reduce the error to 1.8 %.The convergence behavior for dimer embedding is very similar between light and tight settings; in both cases the PBE+MBD MARE of the cell volumes is 33 % larger than the corresponding ME2(4 Å) MARE.
So far, we have only evaluated the performance of the multimer 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 evaluate static lattice energies, sublimation enthalpies can be back-corrected for vibrational contributions 18,50 .Here, we utilize our recently introduced X23b reference data 8 , which also includes a back-correction of experimental volumes in terms of the average thermal expansion of three density functionals, so that the results of lattice relaxations can directly be compared.Our PBE0+MBD/light cell volumes have a MARE of 2.4 %, while PBE+MBD/light has a MARE of 5.9 %, and with ME2(4 Å)/light we obtain a MARE of 4.4 %.In comparison, the PBE0:PBE+MBD approach by Loboda et al. 41 reaches an accuracy of 3.6 % compared to the X23b reference.However, the two approaches are not directly comparable since in this work, we are utilizing light settings for numeric atom-centered basis functions, no pseudopotentials, and standard range-separation parameters for MBD in all calculations.
In terms of lattice energies, our ME2(4 Å) approach utilizing light settings reaches a MAE of 4.9 kJ/mol.We note that while geometries are often already well described at the light level, it is by far not sufficient for obtaining converged values for energetics.Utilizing converged tight species default settings the MAE w.r.t.X23b lattice energies amounts for PBE+MBD (op-TABLE IV.Mean errors (ME) and mean absolute errors (MAE) of the gamma-point vibrational/phonon frequencies of the X23 set in cm −1 and vibrational free energies at 300 K normalized per molecule in kJ/mol (converged q-grid) compared to PBE0+MBD results (light settings).In one case the structures are internally optimized utilizing the PBE0+MBD lattice constants and in the other case the structures are fully optimized.timized) to 3.8 kJ/mol and for PBE0+MBD (calculated on top of PBE0+MBD/light structures) to 3.2 kJ/mol.At fully optimized ME2(4 Å)/tight structures this error even decreases to 2.9 kJ/mol, where we seem to benefit from a certain error cancellation.Small differences in the cell volumes seem to virtually have no effect on the overall accuracy of the corresponding lattice energies when comparing to X23b since the ME2(4 Å)/tight MAE when using PBE0+MBD/light-optimized structures amounts to also 2.9 kJ/mol.In comparison, the PBE0:PBE+MBD approach by Loboda et al. 41 (when correcting the isolated monomer energies for oxalic acid) reaches an accuracy of 3.6 kJ/mol.

C. Embedding for Harmonic Vibrational Properties
Next, we evaluate the performance of our embedding approach for vibrations/phonons using the X23 benchmark set.Therefore, we have calculated with several methods using light species default settings the gammapoint 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 IV).It can be seen that there are substantial 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.Monomer embedding already significantly improves the internal vibrational modes and hence produces frequencies with a MAE 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 statistics, 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.
In order to illustrate the differences in vibrational frequencies, we have plotted in Fig. 3 the phonon density of states between 400 and 610 cm −1 for completely optimized structures of uracil.It can be seen that in this range the peaks for PBE+MBD are shifted quite significantly to lower wave numbers compared to PBE0+MBD.Since the main difference in this range comes from intramolecular modes, monomer embedding already corrects for most of the differences, and dimer embedding leads to a small further improvement.The same plot for the frozen PBE0+MBD cells can be found in Fig. S1 in the SI.However, in this frequency range the effect of the slightly different cell parameters is rather small.In addition, we have also plotted the low-frequency phonon density of states up to about 160 cm −1 for optimized and frozen cells in Figs.S2 and S3 in the SI, respectively.It can be seen that in this frequency range the small differences in cell parameters lead to some frequency shifts (Fig. S2) while -when calculated at the PBE0+MBD cells -all embedding approaches already nicely match the PBE0+MBD result.
Finally, we discuss the accuracy for vibrational free energies.In Table IV 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.Note that all thermal properties (in 10 K steps from 0 to 300 K) as well as all data to further post process the results using phonopy are available in ASE databases 76 .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 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 at room temperature are already tightly converged at a multimer cutoff of 3 Å, increasing it does not lead to further changes in terms of the MAE.Since the vibrational free energy consists of the zero-point vibrational energy (ZPVE) and a thermal contribution, we now analyze the accuracy of the ZPVE for the fully optimized cells.For PBE+MBD the MAE amounts to 9.8 kJ/mol (per molecule), suggesting that most of the corresponding vibrational free energy error at room temperature originates in fact from the ZPVE.After monomers are corrected (ME1) the ZPVE MAE drops down to 0.3 kJ/mol, which is a quarter of the corresponding vibrational free energy error at room temperature.Including dimer corrections then further reduces the ZPVE MAE to 0.2 kJ/mol.

D. Timings
After having discussed the accuracy of our multimer embedding approach, we illustrate on two examples how much computation time can be saved at a certain embedding level compared to the canonical periodic Due to recent advances in the implementation of hybrid density functionals in FHI-aims, the two shown light PBE0+MBD calculations are respectively only 3.4 and 5.7 times more expensive than the corresponding PBE+MBD calculations, which make the light setting extremely useful for comparing with the canonical PBE0+MBD calculations.Given the fast implementation, only a small speedup is possible using light settings and for succinic acid the inclusion of trimer interactions leads already to a similar computation time as the canonical method.However, real computational advantage can be achieved using tight settings, which are typically used for accurate energetics, since in this case PBE0+MBD calculations for the two shown small examples are already almost 100 times more expensive than the PBE+MBD calculations.In addition to the CPU time, larger converged periodic PBE0+MBD calculations are also often prohibitively expensive due to massive memory demands, especially when forces and stress tensors are required.In contrast, multimer embedding does not suffer any real memory issues since the largest hybrid calculation at for instance the ME3 level is an isolated trimer.
In the case of ammonia, which is a highly symmetric crystal with 16 atoms per unit cell and 4 atoms per molecule, tight ME2 calculations are about 70 times faster than the canonical PBE0+MBD approach and even ME3(4 Å) is still about 25 times faster.When the involved molecules get larger and less symmetric -like in the case of succinic acid with 14 atoms per moleculethe speedup is not as massive anymore but ME2(3 Å) calculations are still 10 times faster as the canonical methods and even ME3(3 Å) is 4 times faster.

IV. CONCLUSION
We have introduced trimer interactions and harmonic vibrational properties for a subtractive multimer embedding scheme in order to enable larger calculations for molecular crystals utilizing hybrid density functionals, including a new open-source implementation.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 multimer calculations.Herein, we approximated 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.However, we note that this approach can in principle be used for any combination of methodologies but convergence will most likely be slower when using less compatible methods than PBE0+MBD and PBE+MBD.
The performance of the shown approach was evaluated by directly comparing the multimer embedding results for the X23 benchmark set of molecular crystals with canonical periodic hybrid calculations.In order to accurately approximate lattice energies, the newly incorporated trimer energies are crucial, enabling an agreement within 1 kJ/mol.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 further improve the description of the cell volume.Hence, the next crucial step towards improving this multimer embedding methodology is the explicit inclusion of trimer interactions for gradients and stress tensors and to reduce the number of multimers that need to be considered to improve the efficiency of this methodology.
Furthermore, we have also newly introduced the calculation of vibrational properties utilizing multimer embedding.We are able to approximate gamma-point vibrational/phonon frequencies with an accuracy of a few wave numbers using monomer or dimer embedding.This enables a very accurate approximation of room temperature vibrational free energies within 1 kJ/mol on average when normalized per molecule in the unit cell.
This multimer embedding approach at the dimer level can already for single-point energies be up to 70 times faster than the canonical high-level periodic calculation in the case of ammonia when embedding PBE0+MBD into PBE+MBD using converged tight settings within FHI-aims.Since the largest speedup is observed for small monomers, this could potentially be especially relevant for modeling hydrates.gramme 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).

FIG. 2 .
FIG.2.Relative energies of a cubic ammonia unit cell w.r.t the lattice constant for several embedding approaches (light settings).

FIG. 3 .
FIG.3.Phonon density of states of uracil calculated on top of optimized structures for several methods (only frequencies between 400 and 610 cm −1 are shown).

TABLE I .
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 %.

TABLE II .
Mean absolute errors (MAE) and maximal errors (MAX) of the calculated atomic force components and nonzero stress tensor components of the X23 set compared to PBE0+MBD calculated at PBE+MBD-optimized structures.

TABLE V .
Relative timings of single-point energy calculations normalized to PBE+MBD/light or PBE+MBD/tight calculations of ammonia calculated on 4 cores (light) and 24 cores (tight).