Exploring the Polymorphism of Drostanolone Propionate

2α-Methyl-4,5α-dihydrotestosterone 17β-propionate, known as drostanolone propionate or masteron, is a synthetic anabolic-androgenic steroid derived from dihydrotestosterone. The crystal structures of two polymorphs of drostanolone propionate have been determined by single crystal X-ray diffraction and both crystallizes in the monoclinic crystal system. One is belonging to the P21 space group, Z = 2, and has one molecule in the asymmetric unit while the second belongs to the I2 space group, Z = 4, and contains two molecules in the asymmetric unit. Another polymorph has been investigated by an X-ray powder diffraction method and solved by Parallel tempering/Monte Carlo technique and refined with the Rietveld method. This polymorph crystallizes in the orthorhombic P212121 space group, Z = 4 having one molecule in the asymmetric unit. The structural configuration analysis shows that the A, B, and C steroid rings exist as chair geometry, while ring D adopts a C13 distorted envelope configuration in all structures. For all polymorphs, the lattice energy has been computed by CLP (Coulomb-London-Pauli), and tight-binding density functional theory methods. Local electron correlation methods were used to estimate the role of electron correlation in the magnitude of the dimer energies. The nature of the intermolecular interactions has been analyzed by the SAPT0 energy decomposition methods as well as by Hirshfeld surfaces.


Introduction
The polymorphism in organic compounds is the ability to exist in several crystalline forms. The polymorphs have similar chemical content but exhibit different crystal packing and arrangement [1]. Androgens represent a class of compounds, which can be synthetic or naturally found in vertebrates that are responsible for development and maintenance of male characteristics. The anabolic component of a certain anabolic-androgenic steroid is related to increased protein synthesis in muscle tissues and bones, while the androgenic component is responsible for the development of secondary male characteristics such as deepening of the voice, facial and body hair growth, and aggressiveness. Synthetic androgens are chemically-modified forms of the primary male hormone, testosterone, with the purpose to lower the androgenic characteristics and increase the anabolic properties [2]. Despite the fact that anabolic-androgenic steroids fulfill certain functions in vertebrates and are very effective to boost sports performance, their misuse and abuse can lead to undesirable, serious negative side effects on health. Heart diseases (hypertension, left ventricular hypertrophy, impaired diastolic filling, polycythemia, and thrombosis) are known to be related to the long-term consumption of anabolic steroids [3]. A decrease in the high density lipoprotein (HDL), concomitant with an increase in a low density lipoprotein (LDL) and total cholesterol, are related to the consumption as well, which may anabolic steroids [3]. A decrease in the high density lipoprotein (HDL), concomitant with an increase in a low density lipoprotein (LDL) and total cholesterol, are related to the consumption as well, which may increase the risk of atherosclerosis in the coronary arteries [4][5][6]. Hypogonadism induced by anabolic steroids in males and the development of male characteristics in women and children are regarded as having milder side effects and they are reversible to some extent [7]. Besides the anabolic properties of the derivatives in this group, they can be useful in certain medical conditions, which cause undesired catabolism and loss of muscle mass. They act as anti-catabolic, which negates the effects of cortisol and the derivatives from glucocorticoids class [7,8].
Drostanolone propionate (2α-Methyl-4,5α-dihydrotestosterone 17β-propionate), known with the trade name as masteron (Figure 1b), is an androstane steroid, derived from dihydrotestosterone ( Figure 1a) [9,10]. The scheme of atoms and rings labeling was made according to the established notations of the compounds in this class [9]. It works in the same manner as any androgenic steroid, being an agonist of the androgen receptor [11]. Medically, it was used in breast cancer treatment [9] and, in addition, is often used in sport, bodybuilding, and powerlifting as a performance enhancer, while providing increased protein synthesis, which reflects in the gain of lean muscle mass tissue and aids recovery [9]. Often anabolic-androgenic steroids are esterified with the intention to increase the duration of action by intramuscular or subcutaneous administration [12,13]. Drostanolone propionate is an injectable steroid, which is modified by esterification of the parent hormone (drostanolone) in the position of the O-H hydroxyl group ( Figure 1). Several esterized forms of drostanolone are available on the market. These include drostanolone propionate that has been investigated for treating breast cancer [14], drostanolone enanthate, which by microbial transformation, led to the synthesis of eight potentially anti-cancer metabolites [15]. Since there are no studies regarding the crystal structure of drostanolone propionate and its polymorphism, this paper will focus on the structural aspects of this anabolic-androgenic steroid. For the starting compound (denoted Drost 1), no single crystals were obtained and the structure was solved from powder diffraction data. Using samples from the starting compound (Drost 1), by recrystallization in ethanol, the polymorph Drost 2 was obtained and, by recrystallization of the starting compound in acetone, was obtained by polymorph Drost 3.

Results and Discussion
The chemical configurations of drostanolone and drostanolone propionate are depicted in Figure 1, while overall conformations of the molecules in the asymmetric unit are shown in Figure 2. Since, in Drost 3, there are two molecules in the asymmetric unit (denoted with molecule A and molecule B), in Figure 2c, both molecules are presented. The most important crystallographic data are collected in Table 1.

Results and Discussion
The chemical configurations of drostanolone and drostanolone propionate are depicted in Figure 1, while overall conformations of the molecules in the asymmetric unit are shown in Figure 2. Since, in Drost 3, there are two molecules in the asymmetric unit (denoted with molecule A and molecule B), in Figure 2c, both molecules are presented. The most important crystallographic data are collected in Table 1.     Figure 3 presents the molecular packing of Drost 1, viewed along the a-axis direction. The molecules of Drost 1 are linked by C6-H...O3 contacts formed between the carbonyl O3 oxygen and C6 carbon belonging to the B steroid ring and are building molecular sinusoidal chains in the direction of the b-axis. Along the c-axis, the molecules are linked by C5-H...H-C19 short contacts formed between C5 carbon, which belongs both to A and B rings toward the C19 methyl group. This pattern is repeating by translations equal to lattice parameters both along the c-axis and b-axis.  The unit cell of Drost 2 contains two molecules related by a 21 screw axis linked by C5-H...H-C12 short contact, formed between C5 carbon, which belongs to both A and B steroid rings and C12 carbon of the C ring. The pattern is further extended by translations and forms an arrangement parallel to the ob direction. At the same time, along the c-axis, an infinite arrangement of molecules is formed, which are linked between O3 carbonyl oxygen and C19 methyl group by C19-H...O3 hydrogen bond, as shown in Figure 4. The unit cell of Drost 2 contains two molecules related by a 2 1 screw axis linked by C5-H...H-C12 short contact, formed between C5 carbon, which belongs to both A and B steroid rings and C12 carbon of the C ring. The pattern is further extended by translations and forms an arrangement parallel to the ob direction. At the same time, along the c-axis, an infinite arrangement of molecules is formed, which are linked between O3 carbonyl oxygen and C19 methyl group by C19-H...O3 hydrogen bond, as shown in Figure 4.  Figure 3 presents the molecular packing of Drost 1, viewed along the a-axis direction. The molecules of Drost 1 are linked by C6-H...O3 contacts formed between the carbonyl O3 oxygen and C6 carbon belonging to the B steroid ring and are building molecular sinusoidal chains in the direction of the b-axis. Along the c-axis, the molecules are linked by C5-H...H-C19 short contacts formed between C5 carbon, which belongs both to A and B rings toward the C19 methyl group. This pattern is repeating by translations equal to lattice parameters both along the c-axis and b-axis.  The unit cell of Drost 2 contains two molecules related by a 21 screw axis linked by C5-H...H-C12 short contact, formed between C5 carbon, which belongs to both A and B steroid rings and C12 carbon of the C ring. The pattern is further extended by translations and forms an arrangement parallel to the ob direction. At the same time, along the c-axis, an infinite arrangement of molecules is formed, which are linked between O3 carbonyl oxygen and C19 methyl group by C19-H...O3 hydrogen bond, as shown in Figure 4. O1B oxygen. Within the structure, the B molecules are held together by a combination of C4B-H . . . O1B between C4B carbon of the ring A with the oxygen O1B and the C21B-H . . . H-C22B short contact interaction extended between C21B and C22B carbons of the propionate ester ( Figure 5). Mercury software [16] was used to generate perspective molecular views and the packing diagrams.

Crystal Structure Analysis
Molecules 2020, 25, x FOR PEER REVIEW 5 of 21 The unit cell of Drost 3 fits eight drostanolone propionate molecules. The two independent molecules in the asymmetric unit of Drost 3 are held by bifurcated C5B-H...O1A and C1B-H...O1A hydrogen bonds between the carbonyl O1A oxygen with C1B and C5B carbon atoms of the A ring in molecule B. The C4A-H…O1B hydrogen bond links the C4A carbon of molecule A with the carbonyl O1B oxygen. Within the structure, the B molecules are held together by a combination of C4B-H…O1B between C4B carbon of the ring A with the oxygen O1B and the C21B-H…H-C22B short contact interaction extended between C21B and C22B carbons of the propionate ester ( Figure 5). Mercury software [16] was used to generate perspective molecular views and the packing diagrams.

Crystal Structure Determination of Drost 1 by X-ray Powder Diffraction (XRPD)
The crystallization attempts to obtain suitable single crystals of Drost 1 have failed and it was needed to undertake the crystal structure determination by the XRPD method. This is a multi-step method and involves the following: the X-ray pattern indexing, Pawley refinement, space group assignment, structural model determination, and Rietveld refinement [17].
The pattern indexing was accomplished in the Reflex module, which is implemented in Materials Studio software [18]. Analyzing the unit cell solutions and based on the figure of merit (F.O.M), an orthorhombic common solution was obtained by using four different programs: DICVOL96 [19], TREOR90 [20], ITO15 [21], and X-cell [22]. Taking into account the first 30 lines, a very high FOM is obtained as follows: FOM = 50 with TREOR90, FOM = 49 with DICVOL96, FOM = 9.7 with X-cell, while with ITO15, the FOM = 73 was obtained by indexing the first 20 diffraction lines. The orthorhombic solution (a = 27.2543 Å, b = 12.078 Å, c = 6.4178 Å) has successfully indexed all selected lines.
Pawley refinement procedure (Rwp = 4.02%) has confirmed that the crystal system is orthorhombic P212121 as being the most likely space group with one molecule of drostanolone propionate in the asymmetric unit.
As a starting point for finding the structural model, the CIF (Crystallographic Information File) file of the Drost 2 polymorph was obtained by single crystal X-ray diffraction in order to have the most reasonable bond distances. The model was optimized by altering torsion angles translations and, respectively, rotations. A single molecule was considered in the asymmetric unit because the calculated density using the adopted spatial group is 1.14 g/cm 3 , which is a reasonable value considering the chemical composition of the compound. The structural solution model was determined by using direct-space Monte-Carlo with parallel tempering methods in Fox software [23]. This is a trial and error method, which works by comparison between experimental and calculated patterns and, when they are in agreement, the solution is subjected to the Rietveld refinement.
In the Rietveld procedure, which was carried out with a Reflex module from Materials Studio software, was included as the refinement of the following parameters: diffraction peaks, which were approximated as Pseudo-Voight function, (U, V, W) parameters from the Caglioti formula [24], profile parameters NA, NB in Bragg-Brentano geometry, zero point shift parameter, peaks

Crystal Structure Determination of Drost 1 by X-ray Powder Diffraction (XRPD)
The crystallization attempts to obtain suitable single crystals of Drost 1 have failed and it was needed to undertake the crystal structure determination by the XRPD method. This is a multi-step method and involves the following: the X-ray pattern indexing, Pawley refinement, space group assignment, structural model determination, and Rietveld refinement [17].
The pattern indexing was accomplished in the Reflex module, which is implemented in Materials Studio software [18]. Analyzing the unit cell solutions and based on the figure of merit (F.O.M), an orthorhombic common solution was obtained by using four different programs: DICVOL96 [19], TREOR90 [20], ITO15 [21], and X-cell [22]. Taking into account the first 30 lines, a very high FOM is obtained as follows: FOM = 50 with TREOR90, FOM = 49 with DICVOL96, FOM = 9.7 with X-cell, while with ITO15, the FOM = 73 was obtained by indexing the first 20 diffraction lines. The orthorhombic solution (a = 27.2543 Å, b = 12.078 Å, c = 6.4178 Å) has successfully indexed all selected lines.
Pawley refinement procedure (R wp = 4.02%) has confirmed that the crystal system is orthorhombic P2 1 2 1 2 1 as being the most likely space group with one molecule of drostanolone propionate in the asymmetric unit.
As a starting point for finding the structural model, the CIF (Crystallographic Information File) file of the Drost 2 polymorph was obtained by single crystal X-ray diffraction in order to have the most reasonable bond distances. The model was optimized by altering torsion angles translations and, respectively, rotations. A single molecule was considered in the asymmetric unit because the calculated density using the adopted spatial group is 1.14 g/cm 3 , which is a reasonable value considering the chemical composition of the compound. The structural solution model was determined by using direct-space Monte-Carlo with parallel tempering methods in Fox software [23]. This is a trial and error method, which works by comparison between experimental and calculated patterns and, when they are in agreement, the solution is subjected to the Rietveld refinement.
In the Rietveld procedure, which was carried out with a Reflex module from Materials Studio software, was included as the refinement of the following parameters: diffraction peaks, which were approximated as Pseudo-Voight function, (U, V, W) parameters from the Caglioti formula [24], profile parameters NA, NB in Bragg-Brentano geometry, zero point shift parameter, peaks asymmetry parameters from Berar-Baldinozzi approximation (P1, P2, P3, P4), the parameters that describe preferred orientation (a * , b * , c * and R0) in March-Dollase correction, and the coefficients describing the Molecules 2020, 25, 1436 6 of 20 background profile. The final fit between calculated and measured patterns shows a good agreement, with R wp = 5.65% (Table 2 and Figure 6). asymmetry parameters from Berar-Baldinozzi approximation (P1, P2, P3, P4), the parameters that describe preferred orientation (a * , b * , c * and R0) in March-Dollase correction, and the coefficients describing the background profile. The final fit between calculated and measured patterns shows a good agreement, with Rwp = 5.65% (Table 2 and Figure 6).

Hirshfeld Surfaces
Front views and back views of the three-dimensional Hirshfeld surfaces (mapped with dnorm) for studied drostanolone propionate structures are given in Figure 7. The contacts mapped in red represent the interactions with the distances between atoms smaller than the sum of van der Waals radii, white areas highlight intermolecular contacts close to the sum of van der Waals radii, and blue is used for longer contacts [25]. The Hirshfeld surfaces were mapped with dnorm between −0.151 (red) and 1.628 (blue). Figure 7 shows the Hirshfeld surfaces' front views, respectively, back views of each polymorph, and intermolecular interactions smaller than the sum of van der Waals radii showed as labelled arrows. Each label illustrated on the Hirshfeld surface has the interaction geometry detailed in Table 3.

Hirshfeld Surfaces
Front views and back views of the three-dimensional Hirshfeld surfaces (mapped with d norm ) for studied drostanolone propionate structures are given in Figure 7. The contacts mapped in red represent the interactions with the distances between atoms smaller than the sum of van der Waals radii, white areas highlight intermolecular contacts close to the sum of van der Waals radii, and blue is used for longer contacts [25]. The Hirshfeld surfaces were mapped with d norm between −0.151 (red) and 1.628 (blue). Figure 7 shows the Hirshfeld surfaces' front views, respectively, back views of each polymorph, and intermolecular interactions smaller than the sum of van der Waals radii showed as labelled arrows. Each label illustrated on the Hirshfeld surface has the interaction geometry detailed in Table 3.  In Table 4, the intermolecular contacts X...H/H...X, (where X can be O or C). The first letter signifies the atom located inside the Hirshfeld surface and the second signifies the atom outside the surface.
The fingerprint plot of Drost 1 presents small spikes (labelled as 1 and 2), which are related to O…H/H…O interactions. Spike 2 is caused by the C6-H6A donor group, which is located inside the Hirshfeld surface and the O3 acceptor of the carboxyl group is situated outside the surface. The contact is taking place at distance di + de~2.6 Å. The complementary spike 1 accounts for the same interaction but the O3 acceptor is inside the surface and the donor outside. The spike labelled as 1 is related to H...H interactions, which extends to the smallest distance of di + de~2.2 Å with the hydrogens being inside and outside the surface as well.
The fingerprint plot of Drost 2 show similar spikes (labelled as 1 and 2) as Drost 1, which are characteristic to O…H/H…O inter-contacts. The spike denoted by 2 indicates the C19-H19B donor, which is situated inside the Hirshfeld surface and connects with O3 carboxyl outside the surface. The intercontact distance extends to di + de~2.57 Å. Spike 1 is complementary with spike 2, but the O3 carboxyl is inside the surface, while the donor is situated outside. The H...H interactions, displayed as label 1 on fingerprint plot, are characterized by the sum of di + de~2.17 Å.  In Table 4    The percentage contributions to the Hirshfeld surfaces areas for studied crystals are represented in Table 4. The conclusions resulting from fingerprint plots analysis include the following. (i) The plots shape and features are different in all three compounds and indicate that the supramolecular assemblies are different for each crystal structure.
(ii) The top end values of de and di in the fingerprint plots of Drost 2 are slightly smaller in comparison with Drost 1 and the two independent molecules of Drost 3, which conclude that Drost 2 has higher packing efficiency [26]. This is already in good agreement with the Kitaigorodskii packing index.  with Drost 1 and the two independent molecules of Drost 3, which conclude that Drost 2 has higher packing efficiency [26]. This is already in good agreement with the Kitaigorodskii packing index.

Lattice Energy Evaluation by the Coulomb-London-Pauli (CLP) Method
The CLP method (See details at Section 3.6), which is based on atom-atom type potentials, has been shown that the formation of two drostanolone propionate polymorphs has led to structures that display similar lattice energies (−156.3 kJ/mol in Drost 1; −159.2 kJ/mol in Drost 2 and −151.6 kJ/mol in Drost 3, respectively).
The partitioned packing energies calculated for the polymorphs driven by slow evaporation in ethanol (Drost 2) and acetone (Drost3) shows similar values as the start compound (Drost 1). The total lattice energy breakdown in individual terms is given in Table 5a. The electrostatic terms, which account for Coulombic and polarization term have a small contribution, while the dispersion energy plays the major role. Specific interatomic interactions that contribute mainly into the total intermolecular interaction energy are given in Table 3. Although the O...H contacts have a lower preponderance than the H...H contacts (Table 4), they have a high energy and contribute in the greatest extent to lattice energy.  [28]. Such compounds were investigated by adsorption and Raman scattering, which is assigned to intermolecular interaction and specific effects of crystal packing by comparing experimental spectra with quantum chemical calculations [29].
The Kitaigorodskii packing index [30] is a measure of packing efficiency and usually has a value of 65%. It was evaluated by PLATON software [31] and the results obtained are as follows: 60.13% for Drost 3, 60.25% for Drost 1, and the highest value of 61.64% for Drost 2. Between the packing index and total CLP lattice energy, there is a correlation. The higher the packing index is, the greater the absolute value of the lattice energy is. From Table 5a, it is observed that this correlation between packing index and lattice energy is noticed.

Lattice Energy Evaluation by a Density-Functional Tight-Binding Model
In order to validate our method based on the atom-atom potential-type CLP model, a higher level theoretical method was also considered. Accordingly, the density-functional tight-binding (DFTB) model in its self-consistent charge corrected the variant (SCC-DFTB) [32] implemented in the DFTB+ code [33] and applied in order to compute lattice energies for the three crystal configurations [34]. Several supercells with different cell sizes ( [35] and their electronic energies were computed, including the dispersion effects via the Slater-Kirkwood (SK) dispersion correction scheme [36]. Based on our past experience [37], compared with results computed using the second order Møller-Plesset perturbation theory, the SK-dispersion model can well reproduce the dispersion effects in the case of large molecular clusters. The largest supercell structure (3 × 3 × 3) for Drost 2 crystal conformation includes 87 monomers (5394 atoms) while those for Drost 1 and Drost 3 (2 × 3 × 3) contain 134 (8308 atoms) and 164 monomers (10168 atoms), respectively. Then, the total electronic energy values of the supercells with different cell sizes defined by the number of monomers in the supercell were extrapolated into the infinite large monomer numbers and the lattice energy, as the size independent parameter in the interpolation scheme (parameter a), was obtained [34]. The best performant fitting function (R 2 ≈ 0.94) was found as: where N is the number of the unit cells, a is the fitting parameter, and E latt means the lattice energy for an infinite number of unit cells. The dispersion component of the lattice energy considering the same interpolation scheme (See Equation (1)) was also computed. Accordingly, the lattice total energies and their dispersion parts computed for the three unit cell configurations with the above described calculation scheme are presented in Table 5b. Comparing the results obtained based on the CLP model [38,39] and by SCC-DFTB theory, relatively good agreement between the two theories can be found for the lattice energy values. This agreement looks effective for the case of Drost 2 crystal configuration (the energy difference is 7.5 kJ/mol), while, for the Drost 1 and Drost 3, the energy deviations are a bit larger: 13.7 and 13.9 kJ/mol, respectively. Furthermore, both theoretical methods give the Drost 2 conformation as the most bounded crystal structure followed by Drost 1, while the Drost 3 crystal has the smallest cohesive energy value. Furthermore, both methods suggest that the main attractive forces, which keep the crystalline structures together are the dispersion effects.
In order to better understand the nature of the intermolecular forces, which govern the formation of the molecular crystals, higher level electronic structure theory was also considered. Accordingly, the molecular dimer energies for different dimer geometry conformations taken from the Drost 2 crystal supercell were computed using the second order Møller-Plesset perturbation theory based on localized molecular orbitals and density-fitting technique (DF-LMP2) [40]. The local correlation treatment also offers the possibility to decompose the intermolecular interaction energy into intramolecular, dispersive, and ionic components of the correlation contribution. Five different dimer geometries, relevant for the unit cell configuration, were selected (see Figure 9). The selection criterium is based on the Drost 2 unit cell configuration, where all close contacts (defined as the sum of the vdW radii + 0.2 Å) of the Drost 2 molecule from the unit cell were generated and, from the resulted oligomer cluster, all the possible dimer configurations were kept.
Molecules 2020, 25, x FOR PEER REVIEW 12 of 21 localized molecular orbitals and density-fitting technique (DF-LMP2) [40]. The local correlation treatment also offers the possibility to decompose the intermolecular interaction energy into intramolecular, dispersive, and ionic components of the correlation contribution. Five different dimer geometries, relevant for the unit cell configuration, were selected (see Figure 9). The selection criterium is based on the Drost 2 unit cell configuration, where all close contacts (defined as the sum of the vdW radii + 0.2 Å) of the Drost 2 molecule from the unit cell were generated and, from the resulted oligomer cluster, all the possible dimer configurations were kept. The intermolecular interaction energies between the drostanolone propionate monomers were computed considering the DF-LMP2 method and using the def2-tzvp [41] basis set, as implemented in the Molpro program package [42]. Binding energies were estimated in the supramolecular approximation (i.e., as a difference between the energy of a given dimer complex and the sum of energies of the isolated molecules constituting it). The Drost 2 crystal structure is built by parallel layers grown in the aob crystal plane, where layers are kept together by the side-chain interactions of the drostanolone molecules along the oc crystal axis. Three dimer configurations (see Figures 9a-c) were identified as relevant pair conformations for the layer structure stability (parallel stacking, Tshape, and antiparallel stacking). Their intermolecular interaction energies obtained at the DF-LMP2/def2-tzvp level of theory are: ΔE a) = -21.67 kJ/mol, ΔE b) = -15.37 kJ/mol and ΔE c) = -18.08 kJ/mol, respectively. On the other hand, the interaction between the layers is mainly built by two characteristic pair configurations (See Figures 9d,e). Their intermolecular interaction energies computed at DF-LMP2/def2-tzvp level of theory are: ΔE d) = -8.25 kJ/mol and ΔE e) = -14.13 kJ/mol, respectively. To estimate the contribution of higher (other than pair correlation) electron correlation effects to the intermolecular interaction energy, the coupled cluster level of theory based on localized molecular orbitals and density-fitting technique (DF-LCCSD(T)) for the a) dimer configuration was considered. The result shows that higher level electron correlation effects do not change the magnitude of the intermolecular interaction energy significantly, which gives only a small contribution (+1.18 kJ/mol) to the final energy value (−21.67 kJ/mol for DF-LMP2 versus -20.49 kJ/mol for DF-LCCSD(T)). Since advanced ab initio theories including high-level electron correlation effects are not suitable to directly compute the lattice energies due to the large amount of computation, they can be used indirectly to estimate the lattice energy as a sum of the energies of intermolecular (noncovalent) pairwise interactions between the considered molecule and its neighbors [43]. The intermolecular interaction energies between the drostanolone propionate monomers were computed considering the DF-LMP2 method and using the def2-tzvp [41] basis set, as implemented in the Molpro program package [42]. Binding energies were estimated in the supramolecular approximation (i.e., as a difference between the energy of a given dimer complex and the sum of energies of the isolated molecules constituting it). The Drost 2 crystal structure is built by parallel layers grown in the aob crystal plane, where layers are kept together by the side-chain interactions of the drostanolone molecules along the oc crystal axis. Three dimer configurations (see Figure 9a)-c)) were identified as relevant pair conformations for the layer structure stability (parallel stacking, T-shape, and antiparallel stacking). Their intermolecular interaction energies obtained at the DF-LMP2/def2-tzvp level of theory are: ∆E a) = −21.67 kJ/mol, ∆E b) = −15.37 kJ/mol and ∆E c) = −18.08 kJ/mol, respectively. On the other hand, the interaction between the layers is mainly built by two characteristic pair configurations (See Figure 9d,e). Their intermolecular interaction energies computed at DF-LMP2/def2-tzvp level of theory are: ∆E d) = −8.25 kJ/mol and ∆E e) = −14.13 kJ/mol, respectively. To estimate the contribution of higher (other than pair correlation) electron correlation effects to the intermolecular interaction energy, the coupled cluster level of theory based on localized molecular orbitals and density-fitting technique (DF-LCCSD(T)) for the a) dimer configuration was considered. The result shows that higher level electron correlation effects do not change the magnitude of the intermolecular interaction energy significantly, which gives only a small contribution (+1.18 kJ/mol) to the final energy value (−21.67 kJ/mol for DF-LMP2 versus −20.49 kJ/mol for DF-LCCSD(T)). Since advanced ab initio theories including high-level electron correlation effects are not suitable to directly compute the lattice energies due to the large amount of computation, they can be used indirectly to estimate the lattice energy as a sum of the energies of intermolecular (noncovalent) pairwise interactions between the considered molecule and its neighbors [43]. Accordingly, in the case of the Drost 2 molecule, the total pairwise interaction energy between a monomer from the crystal and its neighbors (close contacts) can be computed as the sum of the binding energies of the five dimers but consider each of them twice, since in the close contact configuration, each of them occurs two times. In this way, the sum of the pairwise interaction energies is −155.00 kcal/mol, which is consistent with previous lattice energy calculations for the Drost 2 case based on the CLP and SCC-DFTB methods (−159.2 kcal/mol and −151.7 kcal/mol, respectively). On the other hand, the Espinoza's empirical formula [44] defined in the framework of QTAIM (or quantum theory of atoms in molecules) theory [45] can be used as another possible solution for estimating the lattice energy. This solution has been successfully applied in several cases [46][47][48] based on Density Functional Theory calculations, but applying it in the case of the DF-LMP2 theory is not a simple task.
Using the CLP and the SCC-DFTB simple models, only the global contribution of the dispersion effects was taken into account. However, it would be interesting to analyze how the dispersion effects manifest along the different crystal axes, since it is shown in Figure 4 that different relative pair conformations can be observed along the oa and ob or oc axes. More precisely, the binding energies of dimers a)-c) build layers along the oab crystal plane, while interactions found in case of dimers d) and e) enhance the crystal cohesion along the oc direction. Accordingly, for the five dimer cases of the Drost 2 crystal configuration, the symmetry adapted perturbation theory (SAPT) method [49] was applied to decompose their intermolecular interaction energy into physically meaningful energy components (electrostatic, exchange, induction, and dispersion) defined in the framework of the SAPT theory. Using the SAPT energy decomposition method, our goal was to get rather realistic values than to directly compare them with results obtained by other theoretical methods, like DF-LMP2. Due to the large numbers of atoms in the dimer geometry, applying the so-called "gold" and "silver" standards for the SAPT method are not feasible. Only the "bronze" standard is feasible. This standard is defined as the zero-order SAPT expansion combined with the exchange-scaling (sSAPT0) approximation [50] and used together with the jun-cc-pVDZ basis set [51]. The intermolecular energy values and their different energy components computed at sSAPT0 as well as the intermolecular energy values computed at DF-LMP2 levels of theory are presented in Table 6. Table 6. The intermolecular energy values and their different energy components computed at sSAPT0 (E electr -electrostatic, E exch -exchange, E ind -induction and E disp -dispersion) as well as the intermolecular energy values obtained at DF-LMP2 levels of theory for the five dimer configurations taken from the Drost2 crystal conformations (all energy values are given in kJ/mol). In the case of dimer configurations having parallel stacking or a), T-shape or b) and antiparallel stacking or c) conformations, the most dominant energy component is the dispersion part, which is almost equal or even lower than the total sSAPT0 energy, but it is almost half balanced by the exchange repulsion effects. The electrostatic and induction contributions individually have a relatively small contribution to the sSAPT0 energy, but their common effect is no longer a negligible contribution. In the case of layer-layer interaction, which includes the d) and e) dimer configurations, in addition to dispersion and exchange interactions, electrostatic contributions have also become important. All these results show us that the layer structure defined by the a)-c) interaction configurations (oa and ob crystal axes directions) are much stronger bound than the interaction between the layers (along the oc crystal axis).

Conformational Analysis
The A, B and C rings of the steroid skeletons were found as chair conformation, whereas D rings are adopting a C13 envelope conformation in all three structures.
The displacement in the ring's geometry from the ideal chair conformation was described by the asymmetry parameter ∆C s evaluation [52]. The calculated values of the asymmetry parameter ∆C s show that the geometry of the A, B, and C rings for all polymorphs are close to the ideal chair configuration, which would display an asymmetry parameter equal to zero. Each ring has three pseudo mirror planes (Figure 10). The calculated asymmetry parameter ∆C s values for each of them are listed for comparison in Table 7. For a ring, the three values of the asymmetry parameter ∆C s have different values and it makes sense to define the average value of the asymmetry parameter as the sum of the asymmetry parameters divided by three.

Conformational Analysis
The A, B and C rings of the steroid skeletons were found as chair conformation, whereas D rings are adopting a C13 envelope conformation in all three structures.
The displacement in the ring's geometry from the ideal chair conformation was described by the asymmetry parameter ΔCs evaluation [52]. The calculated values of the asymmetry parameter ΔCs show that the geometry of the A, B, and C rings for all polymorphs are close to the ideal chair configuration, which would display an asymmetry parameter equal to zero. Each ring has three pseudo mirror planes ( Figure 10). The calculated asymmetry parameter ΔCs values for each of them are listed for comparison in Table 7. For a ring, the three values of the asymmetry parameter ΔCs have different values and it makes sense to define the average value of the asymmetry parameter as the sum of the asymmetry parameters divided by three. The closest value to the ideal geometry is found in ring B of Drost 3, molecule A (ΔCs = 0.3), and, from the data listed in Table 6, it is concluded that, overall, in all structures, the asymmetry parameter ΔCs has the smallest values in ring A, while greater are present in C rings.
The maximal torsional parameter, τm (Table 8), in the studied compounds, was found to be relatively constant and close to 47°, which is a common value found in all D steroid rings [53] and other androstane derivatives [54]. Overlapping Drost 1, Drost 2, and the two molecules in the corresponding asymmetric unit of Drost 3, it is observed that A, B, C, and D rings overlap quite well, while the largest differences are observed in the propanoic acid terminals ( Figure 11). Table 7. Calculated asymmetry parameter values in studied compounds.

Mirror Planes Drost 1 Drost 2 Drost 3-Mol A Drost 3-Mol B
Ring A ΔCs (C3-C10) ΔCs (C4-C1) ΔCs (C5-C2) Average ΔCs   The closest value to the ideal geometry is found in ring B of Drost 3, molecule A (∆C s = 0.3), and, from the data listed in Table 6, it is concluded that, overall, in all structures, the asymmetry parameter ∆C s has the smallest values in ring A, while greater are present in C rings.
The maximal torsional parameter, τ m (Table 8), in the studied compounds, was found to be relatively constant and close to 47 • , which is a common value found in all D steroid rings [53] and other androstane derivatives [54]. Overlapping Drost 1, Drost 2, and the two molecules in the corresponding asymmetric unit of Drost 3, it is observed that A, B, C, and D rings overlap quite well, while the largest differences are observed in the propanoic acid terminals (Figure 11).

Materials
White crystalline powder of drostanolone propionate was received from the Chinese company Wuhan Shu Mai Technology Co. Ltd. (Wuhan, China).

Crystal Growth
Needle-shaped suitable single crystal for X-ray experiments were obtained for Drost 2 in ethanol solution by a slow evaporation method and plate-like crystals from acetone solution for Drost 3. No suitable single crystals were obtained for the starting compound known as Drost 1.

X-ray Powder Diffraction (XRPD)
X-ray powder diffraction pattern of Drost 1 was recorded with monochromatic radiation (CuKα1 radiation) obtained using a germanium monochromator on a Brucker D8 Advance diffractometer (tube operating at 40 kV, 40 mA), equipped with a LYNXEYE detector. The sample was scanned in the range 2θ = 5.5-40° with a step of 0.005 and 3 s/step.

Single Crystal X-ray Diffraction
Suitable single crystals of Drost 2 and Drost 3 were selected and mounted on a SuperNova diffractometer goniometer. The diffractometer was equipped with dual micro-sources, Eos CCD detector, with the experimental data being collected using CuKα radiation. Data collection and data reduction was performed with CrysAllis PRO software [55]. The crystal structures of both polymorphs were solved with Olex2 software [56] by Intrinsic Phasing method with ShelXT [57] structure solution program for Drost 2, while Drost 3 was solved using Direct Methods with SHELXS [58], which are both refined with the ShelXL [59] refinement package using Least Squares minimisation.
All non-hydrogenoid atoms were localized in the Fourier difference map and refined anisotropically with the displacement isotropic parameter Uiso(H) = 1.2Ueq(C) for all CH, CH2 groups and 1.5Ueq(C) for all CH3 groups. Hydrogen atoms were placed in idealised positions and treated as riding as follows: ternary CH refined with riding coordinates (C-H = 0.98 Å), secondary CH2 refined

Materials
White crystalline powder of drostanolone propionate was received from the Chinese company Wuhan Shu Mai Technology Co. Ltd. (Wuhan, China).

Crystal Growth
Needle-shaped suitable single crystal for X-ray experiments were obtained for Drost 2 in ethanol solution by a slow evaporation method and plate-like crystals from acetone solution for Drost 3. No suitable single crystals were obtained for the starting compound known as Drost 1.

X-ray Powder Diffraction (XRPD)
X-ray powder diffraction pattern of Drost 1 was recorded with monochromatic radiation (CuKα1 radiation) obtained using a germanium monochromator on a Brucker D8 Advance diffractometer (tube operating at 40 kV, 40 mA), equipped with a LYNXEYE detector. The sample was scanned in the range 2θ = 5.5-40 • with a step of 0.005 and 3 s/step.

Single Crystal X-ray Diffraction
Suitable single crystals of Drost 2 and Drost 3 were selected and mounted on a SuperNova diffractometer goniometer. The diffractometer was equipped with dual micro-sources, Eos CCD detector, with the experimental data being collected using CuKα radiation. Data collection and data reduction was performed with CrysAllis PRO software [55]. The crystal structures of both polymorphs were solved with Olex2 software [56] by Intrinsic Phasing method with ShelXT [57] structure solution program for Drost 2, while Drost 3 was solved using Direct Methods with SHELXS [58], which are both refined with the ShelXL [59] refinement package using Least Squares minimisation.
All non-hydrogenoid atoms were localized in the Fourier difference map and refined anisotropically with the displacement isotropic parameter U iso (H) = 1.2U eq (C) for all CH, CH 2 groups and 1.5U eq (C) for all CH 3 groups. Hydrogen atoms were placed in idealised positions and treated as riding as follows: ternary CH refined with riding coordinates (C-H = 0.98 Å), secondary CH 2 refined with riding coordinates (C-H = 0.97 Å), and idealized CH 3 methyl groups refined as a rotating group (C-H = 0.96 Å).

Evaluation of Intermolecular Interactions by Hirshfeld Surfaces and Fingerprint Plots
The analysis of three-dimensional Hirshfeld surfaces mapped with the d norm function [60] offers the possibility to compare intermolecular contacts based on van der Waals radii in an interactive way with a red, white, and blue color mapping on the surface [61]. Hirschfeld surfaces and fingerprint plots are tools used in order to explore intermolecular interactions in polymorphs [62] and were generated in Crystal Explorer17 [63] software. The software uses CIF files as input and, during the computation, the C-H bond lengths are moved to the well-known standard distances determined by neutron diffraction (C-H = 1.083 Å). The d norm function can be expressed as Equation (2), where d i is considered the distance from the surface to the atom inside the surface, while d e represents the distance from the surface toward the exterior atom and the van der Waals radii for both atoms inside and outside the surface.

Lattice Energy Evaluation by the CLP Model
The crystal lattice energy was calculated using the Coulomb-London-Pauli approximation, developed by Gavezotti and implemented in CLP software [64] as a sum of the following terms (Equation (3)), where I and j represent pairs of atoms belonging to two different molecules.
In Equation (3), the first term is the Coulombian energy, which is treated according to Coulomb's law. The second term is the polarization energy and is treated in the approximation of a linear dipole and depends on the inverse fourth power of distance. The third term represents the dispersion term and is approximated by the inverse of the distance to the sixth power. The last term represents the repulsion energy and is considered as a modulation of the overlapping wave function and depends on the distance between atoms at the twelfth-power. The F P , F D , F R coefficients in Equation (3) are empirical scale parameters and the P ij , D ij , and T ij coefficients are dependent on the vicinity of the atom in the molecule. F Q coefficient is included in Equation (4), q i being the rescaled net charge population on atom i, and q 0 i is the charge in each atomic basin.

Evaluation of Intermolecular Interactions by First Principal Methods
The intermolecular interaction energies between the drostanolone monomers were computed considering the DF-LMP2 method and using the Def2-TZVP [40] basis set as implemented in the Molpro program package [41]. To compute lattice energies for the three crystal configurations [34], the density-functional tight-binding (DFTB) model in its self-consistent charge corrected variant (SCC-DFTB) [32] implemented in the DFTB+ code [33] was applied. The sSAPT0 energy values were computed using the PSI4 software [65].

Conformational Analysis
The conformational analysis of the A, B, and C steroid rings has been accomplished based on the asymmetry parameter ∆C s , which is defined in Equation (5) [46], where φ i and φ i represent the symmetry related torsion angles and m represents the pair numbers of torsion angles (Figure 12a). They have alternate opposite signs. For steroids, the torsion angle θ 0 is considered as the angle corresponding to the common side between C and D rings. torsion angle τm defined as Equation (7) and were derived by C. Altona and M. Sundaralingam [66]. In Equation (6), the torsion angles , , , , and related to the five-membered D ring are considered as shown in Figure 12b.

Conclusions
Starting from drostanolone propionate, two more polymorphs were obtained by recrystallization. The crystal structure for the starting compound was determined by X-ray powder diffraction using Parallel Tempering and refined by the Rietveld method, while, for two single crystals, the structures were determined by single crystal X-ray diffraction. The steroid skeletons were found to be very similar in all three polymorphs, with the A, B, and C rings depicting chair geometries and D rings envelope conformations. The most noticeable difference in molecular configurations were found to be in the propanoic acid terminals. Based on the calculation of the lattice energies with the CLP method, it was concluded that polymorphs obtained by recrystallization have almost the same lattice energies with respect to the starting compound. Hirschfeld surfaces analysis proved that the aggregates are held by combinations of C-H…O hydrogen bonds and by weak C-H…H-C contacts. The evaluation of the fingerprint plots showed that the interactions in the  The puckering in the five membered D rings of the steroid skeleton can be characterized by two distinct parameters: phase angle of pseudorotation, P defined in Equation (6) and the maximum torsion angle τ m defined as Equation (7) and were derived by C. Altona and M. Sundaralingam [66]. In Equation (6), the torsion angles θ 0 , θ 1 , θ 2 , θ 3 , and θ 4 related to the five-membered D ring are considered as shown in Figure 12b. τ m can be roughly approximated as the angle between the C13-C14-C17 and C14-C15-C16-C17 planes.

Conclusions
Starting from drostanolone propionate, two more polymorphs were obtained by recrystallization. The crystal structure for the starting compound was determined by X-ray powder diffraction using Parallel Tempering and refined by the Rietveld method, while, for two single crystals, the structures were determined by single crystal X-ray diffraction. The steroid skeletons were found to be very similar in all three polymorphs, with the A, B, and C rings depicting chair geometries and D rings envelope conformations. The most noticeable difference in molecular configurations were found to be in the propanoic acid terminals. Based on the calculation of the lattice energies with the CLP method, it was concluded that polymorphs obtained by recrystallization have almost the same lattice energies with respect to the starting compound. Hirschfeld surfaces analysis proved that the aggregates are held by combinations of C-H . . .