The interest in the structure and thermodynamic properties of interfaces stems from the fact that they are the place where a myriad of natural phenomena and technologically important processes take place. Knowledge of interfaces behaviour is also essential in the design of the new generation of optical materials. These materials have interesting optical properties [1] stemming from a complex arrangement in space of elementary components at the colloid scales (i.e., from tens of nanometers to micrometers). Studies have focused on different spatial arrangements of such structures, and one of the most promising is the Cubic Diamonds crystal (CDc) [1, 2]. CDc stability and kinetic attainability  with respect to liquid and other competing ordered lattices depend critically on its Surface Free Energy (SFE).

One of the tools available to describe the properties of the interfaces is represented by Molecular Dynamics (MD) calculations, where the time evolution of each component of the system (being it a molecule composed of atoms or a larger object described in a coarse-grained fashion) can be followed in detail. However, their applications to large molecules available in the literature mainly focus on the description of the molecule's spatial structure, whereas most of the MD models dealing with the description of surface properties have been developed for  atomic systems.

The standard methodologies to obtain surface properties do not handle the additional difficulty imposed by the study of (macro-)molecular systems, because of their molecular size and complex interactions. In 2022, Di~Pasquale and Davidchack published an extension of the cleaving technique (sketched in Fig. 1) which can be used to calculate the SFE between two different phases for molecular crystals [3]. 

Figure 1: Sketch of the cleaving model 

In its most general form, we start (Step1) from two different phases (labelled ɑ and β) without any interface. We cleave the two systems (Step 2) by inserting a ‘wall’ (as some external potential). The four parts thus created are then rearranged in Step 3 to put the two phases in contact. The walls are removed (Step 4), and we obtain two interfaces between the phases ɑ and β. If these operations are carried out reversibly, the SFE is then the total work needed to perform such operations.

The generality of such framework has two advantages: i) while the first incarnation of the cleaving methodology [3] uses the Thermodynamic Integration (TI) methods, it can be easily extended to Non-Equilibrium MD simulations which are preferable in systems where the large relaxation times can make TI methods inefficient; ii) it can be easily applied to different systems, and therefore can be of importance to different communities which need to characterise the interface properties of the system they are studying.

Following these considerations, the addition of a package for the calculation of surface properties to a widely used MD code such as LAMMPS is essential. Indeed, such a package could be used by different communities to investigate systems where interfaces are important, regardless of the level of detail at which they are described (atoms, molecules, colloids, etc.).