Molecular dynamics simulations on the effect of solvent and supersaturation on the aggregation behaviour in carbamazepine

Abstract

Polymorphism in carbamazepine is a well-known phenomenon that can occur due to different factors, such as solvent, supersaturation, temperature and external surface. We report molecular dynamics simulation studies on the effect of supersaturation and solvent on the formation of different types of carbamazepine (CBZ) dimers as a function of supersaturation in different solvents. The understanding of the factors that affect the formation of different types of CBZ dimers in solution could help to gain insight into factors that control the nucleation of different polymorphs of carbamazepine in solvent-dependent crystallization. We propose that the presence of different types of CBZ dimers (FI-type, FII-type and FIII-type) with different angles between the 6-membered aromatic rings of different CBZ molecules act as growth synthons for different CBZ polymorphs in solution. The simulation results show that there are four types of species depending on hydrogen bonding in the CBZ clusters – single hydrogen bond (HB), HB dimer, HB trimer and HB tetramer, based on the strength of the solvent–solvent and solvent–CBZ interactions. The HB dimers in the CBZ molecules lead to overlap between the 6-membered rings of CBZ and thus, the formation of FIII-type CBZ dimers in all the supersaturations in acetone and in the high supersaturation system in ethyl lactate. In the presence of a single HB and depending on the strength of the solvent–CBZ interactions there are two possibilities – FII-type CBZ dimers in solvents like ethanol and anisole that have similar CBZ–CBZ and solvent–CBZ interactions or FI-type CBZ dimers in N,N-dimethylformamide if the solvent–CBZ interactions are strong. Solvent–solvent interactions mainly govern the dynamics of the CBZ molecules and CBZ–solvent interactions play an important role in the size of the CBZ clusters.

---

Introduction

Polymorphism has an essential role in the pharmaceutical industry as the bioavailability of pharmaceutically relevant polymorphs could be affected due to different physicochemical properties of the polymorphs.^1–3 The solution crystallization technique is often employed in crystallization of organic molecules, where the solvent plays a substantial role in the polymorphism of the organic compound.^4–8 Carbamazepine (CBZ) is an anticonvulsant drug that is used to treat epilepsy, bipolar disorder and trigeminal neuralgia and shows polymorphism. Carbamazepine has five anhydrous (FI, FII, FIII, FIV and FV) polymorphs and among these polymorphs, the FIII (monoclinic) form is the only bioavailable and thermodynamically stable API form. On the other hand, FI (triclinic), FII (trigonal), FIV (C-centered monoclinic) and FV (orthorhombic) are metastable polymorphs of carbamazepine.^9,10 The FIV and FV forms of CBZ are formed via template growth from a polymeric substrate and single crystal of dihydrocarbamazepine, respectively.^11,12 The stability of the CBZ polymorphs observed in solution decreases as FIII > FI > FII.^13,14 The FI form of CBZ is enantiotropically related to the FIII form with a transition temperature of 78 °C.¹⁴ Polymorphism in CBZ occurs due to several factors, such as supersaturation, solvent, temperature, functional groups of the surface and type of surface.^15,16 Ethyl acetate and 2-butanone exhibit nucleation of only the FII form of CBZ, whereas ethanol and methanol lead to concomitant nucleation of both the FII and FIII forms of CBZ in a medium with supersaturation of 2.0.¹⁵ Carbamazepine shows supersaturation-dependent and solvent-independent polymorphism in cumene where different polymorphs were observed with a change in supersaturation, S and temperature, T; FIII (S < 2.0), FII (S = 2.0 and T = 60 °C), mixture of different polymorphs (60 °C < T < 80 °C) and FI (T > 80 °C, independent of S).¹⁷ Similarly, saturation-dependent polymorphism in carbamazepine was observed in ethanol where concomitant crystallization of both FI and FII forms occurred when S = 1.63, the FIII form when S < 2.0 and concomitant crystallization of the FII and FIII forms when S = 2.0 and polymorphic transformation from FIII to FII was observed when 2.0 < S < 4.0.¹⁸ Interestingly, only the FIII form of CBZ was observed in acetone over a range of supersaturations; 1.6 < S < 4.0.¹⁹

An atomic level understanding of the factors that lead to polymorphism in organic compounds could be beneficial to design experiments for selective crystallization of the desired polymorphs. However, experimental insight into nucleation of small organic molecules is a challenge as nucleation in these systems occurs at smaller length and time scales that are beyond the scope of current state-of-the-art experimental techniques.²⁰ The classical nucleation theory (CNT) assumes that nucleation occurs via a single-step due to simultaneous fluctuations in density and structure that lead to a crystalline nucleus.²⁰ However, recent studies have shown that nucleation in solution occurs via two- or multi-nucleation steps with initial formation of loose aggregates due to density fluctuations and later formation of a crystalline nucleus due to structural fluctuations.^21–24 The nucleation rate equation does not include solute–solvent interactions, and thus cannot capture the role of the solvent in the nucleation in solvent-dependent crystallization of organic molecules.¹⁵ Supramolecular synthons are kinetically formed units, also called growth synthons formed in solution prior to nucleation, and could be related due to the type of interactions with the repeating units in a crystal, known as the structural synthon.^25–28 These growth synthons act as primary units and could participate in the formation of pre-nucleation clusters that would lead to nucleation. Thus, understanding the formation of supramolecular synthons as a function of supersaturation and solvent could be a key to understanding the nucleation of different polymorphs. The molecular dynamics (MD) simulation techniques can be employed to study the formation of supramolecular synthons of API in solution as a function of supersaturation, solvent and temperature. These studies could provide molecular-level insight into factors that influence the formation of supramolecular synthons in solution.^29–31 We report molecular dynamics simulation studies to elucidate the role of supersaturation and solvent in the formation of pre-nucleation aggregates of carbamazepine in a homogeneous medium. The present work explores the atomic factors that affect the formation of FI-, FII- and FIII-type CBZ dimers in solution. FIV- and FV-type polymorphs of CBZ do not crystallize in homogeneous medium and thus, are not studied in the present work.

Simulation details

The simulation systems consist of a supersaturated solution of CBZ in nine different solvents – methanol (METH), ethanol (ETH), tetrahydrofuran (THF), anisole (ANI), acetone (ACE), 2-butanone (2-BUT), N,N-dimethylformamide (DMF), ethyl acetate (EA) and ethyl lactate (EL). The total number of CBZ molecules and solvent molecules for different supersaturations and the final simulation box dimensions of all the systems are reported in Table 1.^15,32 The solubility of CBZ in ethyl lactate is not reported in the literature, and thus three model systems were chosen and labelled as M3, M4 and M5. The polymer consistence force field (PCFF) parameter was employed to model carbamazepine (CBZ) and all the solvent molecules.³³ PCFF is a second-generation forcefield that includes cross-term potentials for bond, angle and dihedral terms and models van der Waals interactions as a 9-6 Lennard Jones potential. All the simulations were performed using the LAMMPS package.³⁴ The timestep for integration was chosen as 1.0 fs. The cut-off distance for both van der Waals and electrostatic interactions was chosen as 14 Å respectively and pppm style was chosen for the calculation of long-range Coulomb interactions. The relaxation times for the Nośe-Hoover barostat and thermostat were chosen as 2 ps and 0.04 ps, respectively. All the systems were initially simulated at 300 K using the NPT ensemble for 4 ns followed by NVT simulations for 4 ns. Furthermore, all the systems were simulated at 265 K using the NPT ensemble for 2 ns followed by NVT simulations for 60 ns. The NVT simulations in a few of the solvents (ethanol, acetone, anisole and ethyl lactate) were further continued for another 10 ns at 265 K. The data were stored at every 1 ps during the NVT simulations at 265 K and were used to calculate the radial distribution and correlation function plots. The final configurations from 265 K NVT runs were used as initial configurations to perform NPT simulations for 2 ns followed by NVT simulations for 20 ns being performed at 298 K in all the systems. The final configurations from 60 ns of NVT simulations were used as the initial configurations to run NVT simulations for 50 ps with a storage interval of 1 fs to store velocities of the lowest supersaturation systems to calculate velocity autocorrelation functions. The final configurations from NVT simulations at 265 K and 298 K, respectively, were used to perform NVT simulations at 50 K for 1 ps to calculate the angle distributions between the 6-membered aromatic rings of different CBZ molecules. The energy analysis was performed using group/group keyword in LAMMPS by simulating the final configurations of 265 K data at 50 K followed by 5 K for 0.5 ns using NVT simulations. The different types of intermolecular interaction energies were also calculated for crystalline forms of FI, FII and FIII of CBZ using 576, 576 and 600 molecules, respectively.

Table 1 Number of carbamazepine (N_CBZ) and solvent (N_Solv) molecules, and final box dimensions (L_x, L_y and L_z) after NPT simulations in systems with different supersaturations (S) and solvents

Solvent	S	N CBZ	N Solv	L x = Ly = Lz (Å^3)
Methanol (METH)	1.5	41	1861	52.3098
Ethanol (ETH)	2.0	17	1289	48.7586
	2.8	24	1289	49.0308
Acetone (ACE)	1.5	16	2032	58.9094
	2.0	21	2032	59.0572
	3.0	31	2032	59.3598
	4.0	41	2032	59.6194
2-Butanone (2-BUT)	1.5	29	1300	55.5670
N,N-Dimethylformamide (DMF)	2.0	30	1861	60.9804
	2.6	40	1861	61.2874
Anisole (ANI)	2.5	20	1000	58.1752
	3.7	30	1000	58.3264
	5.0	41	1000	58.5319
Tetrahydrofuran (THF)	1.5	25	928	52.0416
Ethyl acetate (EA)	3.5	26	1542	62.8562
	5.0	38	1542	63.1506
Ethyl lactate (EL)	M3	27	1542	64.3334
	M4	32	1542	66.6382
	M5	48	1542	64.8540

---

Results and discussion

Interactions between carbamazepine molecules in anhydrous polymorphs of CBZ could be analysed in terms of hydrogen bonding between atoms of the carboxamide group of CBZ and pi–pi interactions between 6-membered rings (6-M) of CBZ molecules as shown in Fig. 1(a–c). There are both symmetric (FII and FIII forms) and asymmetric (FI) hydrogen bonding between CBZ dimers in different polymorphs where the hydrogen bonding distances, d_HB are 1.98 Å and 2.09 Å in FI, 1.94 Å in FII and 2.11 Å in the FIII form of CBZ, respectively. The pi–pi interactions between six-membered (6-M) rings of CBZ molecules occur at a minimum distance, d_PP of 5.05 Å, 5.25 Å and 3.93 Å in the FI, FII and FIII forms of CBZ, respectively. Table 2 shows different intermolecular interaction energies per CBZ molecule where hydrogen bond energy, E_HB is the most favourable in the FI form (−34.09 kcal mol⁻¹) followed by FII (−28.22 kcal mol⁻¹) and less favourable in the FIII form (≈−25 kcal mol⁻¹) of CBZ. These results are consistent with the earlier reported simulation studies where hydrogen bonding between CBZ molecules is more favorable in the FII than the FIII CBZ polymorph.³⁵ On the other hand, pi–pi energy, E_PP is the most favourable in the FIII form (−37.29 kcal mol⁻¹), less favourable in other polymorphs of CBZ between −22.0 kcal mol⁻¹ < E_PP < −26.0 kcal mol⁻¹. These results suggest that the orientation of the 6-M rings in the CBZ dimers is unique for FIII as compared to other anhydrous forms of CBZ. The total energy, E_tot between CBZ molecules also shows a trend similar to E_PP and is the most favourable in the FIII form of CBZ among the four anhydrous polymorphs of CBZ. The hydrogen bond energy is more favourable than the pi–pi energy in the FI and FII forms of CBZ and the opposite is observed in the FIII form of CBZ. The hydrogen bond (HB) distance between the carboxamide atoms of the CBZ molecules shows a small difference among the crystalline forms of CBZ and thus, the HB distance would not be a better choice as a quantitative parameter to differentiate CBZ dimers of different forms in solution. We chose angle, θ between the 6-M rings of the CBZ molecules as a quantitative parameter to differentiate the FI, FII and FIII polymorphs of CBZ. The angle, θ was calculated as a scalar-product between two unit vectors (v_i) that are perpendicular to the planes of the 6-M rings that are at a minimum distance in a CBZ dimer of a given polymorph. The unit vector (v_i) is calculated as a cross product between two vectors, l₁ and l₂ between the chosen carbon atoms of the 6-M ring as shown in Fig. 1d. Fig. 1e reports the probability of angle distribution, P(θ) between different polymorphs of CBZ where θ is between 90°–110° in FI, 60° in FII and 150° in the FIII form of CBZ and we assume these dimers could be one of the possible growth synthons involved in the polymorphism in carbamazepine.

Fig. 1 Hydrogen bonding and pi–pi overlap between dimers in (a) FI, (b) FII and (c) FIII of carbamazepine. Here pi–pi overlap is shown between two CBZ molecules with minimum distance between centroids of 6-membered rings between CBZ molecules as 5.05 Å, 5.25 Å, 3.93 Å and 3.81 Å in the FI, FII and FIII forms of CBZ, respectively, (d) schematic showing the distance between two 6-membered rings (d_PP) of a carbamazepine dimer and vectors l₁ and l₂ used to calculate the unit vector v_k (k = i, j) from the plane of the 6-membered phenyl ring of a carbamazepine molecule and (e) angle distribution between the 6-membered rings at a particular minimum distance, d_PP between different carbamazepine molecules in the F1, FII and FIII forms of carbamazepine, probability of angle distribution, P(θ).

Table 2 Different types of intermolecular interaction energies per CBZ molecule, E_HB, E_PP and E_tot (kcal mol⁻¹) between carbamazepine molecules in different polymorphs (FI, FII and FIII) of carbamazepine

CBZ polymorph	E HB (kcal mol^−1)	E tot (kcal mol^−1)	E PP (kcal mol^−1)
FI	−34.09	−24.24	−8.47
FII	−28.23	−22.36	−8.51
FIII	−24.75	−37.30	−11.80

---

The dynamics of carbamazepine molecules in solution is one of the important factors that could affect the aggregation behaviour of CBZ molecules in solution. Fig. 2a reports the velocity autocorrelation function (VACF) plots of CBZ molecules in low supersaturation systems in different solvents as the probability of large-sized CBZ clusters that could affect diffusivity is low at low supersaturation. The velocity autocorrelation function was calculated using eqn (1), where v_i(0) and v_i(t) are the center-of-mass velocity of the ith carbamazepine molecule at time 0 and t and the diffusion coefficient was calculated using the Green–Kubo relation (eqn (2)).

C_v(t) = 〈v_i(0)·v_i(t)〉 (1)

Fig. 2 (a) Velocity autocorrelation function plot of carbamazepine (CBZ) molecules, (b) diffusion coefficient of CBZ molecules at the lowest supersaturation in different solvents, (c) E_SS (kcal mol⁻¹) solvent–solvent interaction energy per solvent molecule in low supersaturation systems in different solvents and (d) largest size (n^max_cls) of the CBZ cluster as a function of supersaturation in different solvents.

The mobility of CBZ molecules in ethyl lactate is restricted due to large cage effect and the opposite was observed for CBZ molecules in acetone, methanol and ethanol. There is a weak cage effect for CBZ molecules in 2-butanone, ethyl acetate, anisole and N,N-dimethylformamide. Fig. 2b reports self-diffusivity, D of CBZ molecules obtained by integration of VACF plots where the largest and the lowest D values of 13.92 × 10⁻⁷ m² s⁻¹ and 0.71 × 10⁻⁷ m² s⁻¹ were observed in acetone and ethyl lactate, respectively. The range of D values in all other solvents varied from 7.00 to 8.00 × 10⁻⁷ m² s⁻¹ except for 4.9 × 10⁻⁷ m² s⁻¹ in DMF and 6.3 × 10⁻⁷ m² s⁻¹ in ethyl acetate. The diffusion coefficient is not reported for CBZ in tetrahydrofuran due to the presence of large fluctuations in the VACF plot in this system. Fig. 2c reports the solvent–solvent interaction energy per solvent molecule, E_SS, which varies from −67.56 kcal mol⁻¹ in ethyl lactate to −11.4 kcal mol⁻¹ in ethanol and E_SS decreases as EL ≫ EA > DMF > METH > 2-BUT ≈ ACE ≈ ANI > THF > ETH. This is consistent with an increase in the diffusivity of CBZ molecules with a decrease in E_SS and suggests that solvent–solvent interactions are one of the important factors that control the dynamics of CBZ molecules in homogeneous medium. Fig. 2d shows the largest size of CBZ cluster, n^max_cls observed in different systems and the probability distribution plots of cluster size in different systems are reported in Fig. S1. The size of a CBZ cluster was calculated by choosing a minimum cut-off distance, d_cut ≤ 9.0 Å between the center-of-masses of CBZ molecules. The radial distribution function plots between the center-of-mass of CBZ molecules in different systems show no RDF peak at r > 9.0 Å as shown in Fig. S2 and thus, this was chosen as the cut-off distance for cluster size calculations. The size of the largest CBZ cluster, n^max_cls is the same for different supersaturations in a solvent with the exception of ethanol and anisole where aggregate size increases with supersaturation. The largest aggregate with n^max_cls = 10 was observed in the 2-butanone (S = 2.0) and anisole (S = 5.0) systems. On the other hand, the size of the largest CBZ cluster varies between 3 < n^max_cls < 5 in all other solvent systems.

The role of hydrogen bonding between CBZ aggregates is also an important factor that could affect polymorphism. Fig. 3 reports the distribution of hydrogen bond (HB) distances, P(d_HB) and angles, P(θ_HB) between CBZ dimers in different systems. The distributions of HB distances and angles were similar for different supersaturations in a solvent. Hence, we report the distributions of HB distances and angles between CBZ molecules in low supersaturation systems. The minimum value of the hydrogen bond distance, d_HB was observed at 1.5 Å in most of the solvents with the exception of 1.75 Å in alcohol-based solvents. The largest value of d_HB of 2.25 Å was observed only in DMF-based systems as compared to 2.0 Å in other solvent systems. The hydrogen bond angle distribution shows a peak at θ = 165° in acetone-based and ethyl lactate-based systems though the angle distribution is narrower in the former than the latter system. The angle distribution in P(θ_HB) is broader in all other solvents with θ < 150°. The carbamazepine molecules show four different types of hydrogen bondings in different solvents as reported in Fig. 4: (i) a single hydrogen bond between oxygen and one of the hydrogens of carboxamide groups, (ii) hydrogen bond, HB dimer, (iii) HB trimer and (iv) HB tetramer. The hydrogen bond dimers were observed only in S = 2.0 and 3.0 supersaturations in acetone; however, both HB dimers and single hydrogen bonds between CBZ molecules were observed in the S = 1.5 and 4.0 systems in acetone. The observed hydrogen bond angle of θ = 165° in acetone-based systems is expected due to the presence of HB dimers. Hydrogen bond (HB) trimers were observed in ethyl acetate (S = 5.0), anisole, tetrahydrofuran and 2-butanone and only S = 4.0 supersaturation in acetone. The hydrogen bonding between CBZ molecules in a HB trimer occurs due to both the carboxamide hydrogens of one of the CBZ molecules and one of the carboxamide hydrogens of the other two CBZ molecules. On the other hand, only one of the carboxamide hydrogens of each of the four CBZ molecules participates in hydrogen bonding in the HB tetramer and this type of hydrogen bonding was observed only in the S = 1.5 system in acetone. The presence of the HB tetramer of CBZ is in agreement with experimental reports of the observation of a CBZ tetramer from STM images of monolayer deposition of the CBZ solution on the Au(111) surface.³⁶ The alcohol-based solvents show only single hydrogen bonds (HB) between CBZ molecules. The CBZ molecules in tetrahydrofuran show single HBs and HB dimers. Ethyl lactate-based systems show both single hydrogen bonds and HB dimers whereas only one CBZ dimer was observed in the S = 2.6 system in DMF. The presence of hydrogen bond dimers leads to pi–pi overlap between the 6-M rings of the CBZ molecules as shown in acetone in Fig. 4c and d. On the other hand, the presence of a single HB and HB trimers leads to different orientations among the 6-M rings of the CBZ molecules.

Fig. 3 Probability distribution plots for hydrogen bonding (a) distance and (b) angle between carbamazepine molecules in low supersaturation systems in different solvents.

Fig. 4 VMD snapshots of different types of hydrogen bonds (HB); (a) single HB, (b) HB dimer, (c) HB trimers and (d) HB tetramer between carbamazepine molecules in different solvents (METH, ETH, EL, ANI, THF, EA, 2-BUT, DMF and ACE represent methanol, ethanol, ethyl lactate, anisole, tetrahydrofuran, ethyl acetate, 2-butanone, N,N-dimethylformamide and acetone respectively).

Fig. 5 reports the radial distribution function plots between carboxamide oxygen atoms of CBZ molecules, g_O–O(r). The first RDF peak is broad with a peak at r > 4.0 Å in alcohol-based systems due to single hydrogen bonds as shown in Fig. 5a. The intensity of the g_O–O(r) plot is low in methane as in the S = 1.5 system in methanol and DMF-based systems (Fig. 5d) due to a smaller number of hydrogen bonds between CBZ molecules in these systems. Fig. 5b shows that the first RDF peak in the g_O–O(r) plot occurs at 3.73 Å in acetone-based systems due to hydrogen bond dimers and a small shoulder peak at r > 5.0 Å due to the carboxamide oxygen of neighboring CBZ molecules in a cluster. Similarly, Fig. 5c shows the presence of dual peaks in the S = 5.0 system in anisole due to neighboring CBZ molecules, but the intensity of the second peak is low in anisole (S = 2.5 and 3.7) and tetrahydrofuran systems. The g_O–O(r) plot is broad in 2-butanone (Fig. 5b) and ethyl acetate systems (Fig. 5d) due to the presence of HB trimers in these systems. There are two sharp peaks at 3.7 Å and 5.0 Å in the S = 2.5 system in the ethyl lactate systems (Fig. 5e) due to the presence of both single HBs and a HB dimer, respectively.

Fig. 5 Radial distribution function plots, g_O–O(r) between carboxamide oxygens of carbamazepine molecules in different systems; (a) methanol (METH) and ethanol (ETH), (b) acetone (ACE) and 2-butanone (2-BUT), (c) anisole (ANI) and tetrahydrofuran (THF), (d) N,N-dimethylformamide (DMF) and ethyl acetate (EA) and (e) ethyl lactate (EL).

Solvent–CBZ interactions were elucidated by calculating the radial distribution function plots between the carboxamide oxygen (O) and hydrogens (H₁/H₂) of CBZ and polar atomic sites (oxygen, O_S and hydrogen, H_S) of the solvent as shown in Fig. S3. Table 3 reports the distance, r_1st at which the maximum in the first RDF peak was observed in the CBZ–solvent radial distribution function plots. The first RDF peak maximum in the RDF plot of H_S–O occurs at r_1st < 2.0 Å in alcohol-based systems, and M4 and M5 systems in ethyl lactate. Similarly, r_1st < 2.0 Å occurs for O_S–H₁ and O_S–H₂ RDF plots in N,N-dimethylformamide and ethyl acetate-based systems, and thus both CBZ–CBZ and CBZ–solvent hydrogen bonding occurs in these solvents due to which the probability of formation of hydrogen bond dimers is less in these solvents. The solvent-drug RDF plots in acetone, tetrahydrofuran, anisole, and the M3 system in EL show the first RDF peak at r_1st > 2.0 Å, due to which only CBZ–CBZ hydrogen bonding is observed in these systems that enhances the probability of the formation of hydrogen bond dimers. Fig. S4a shows that drug-solvent interaction energy (E_DS) decreases as EA > DMF > EL > METH ≈ ETH > ACE > ANI > 2-BUT > THF, which is consistent with the trend in which r_1st changes in these solvents. The interaction energy between CBZ and the solvent molecules is similar to the change in supersaturation in a solvent and increases with time with the exception of the EA-based systems where E_DS is more favourable for the S = 5.0 than S = 3.7 system with time. Fig. S4b reports that the CBZ–CBZ interaction energy per CBZ molecule, E_DD < −110 kcal mol⁻¹ in most of the systems. The most favorable CBZ–CBZ interactions occur in tetrahydrofuran (E_DD = −183.86 kcal mol⁻¹), followed by ethyl acetate-based (E_DD = −116.28 kcal mol⁻¹) systems and are least favourable in the N,N-dimethylformamide (S = 2.6, E_DD = −71.84 kcal mol⁻¹) system. Thus, CBZ–CBZ interactions are the most favourable in systems that have the least favourable CBZ–solvent interaction energies.

Table 3 Distance, r_1st at which the maximum in the first peak occurs in the radial distribution function plot between atomic sites of the solvent (H_S or O_S) and carboxamide hydrogens (H₁ and H₂) in different systems. Also reported are relaxation times (τ₁ and τ₂) for second order reorientational correlation function plots for O–H₁ and O–H₂ hydrogen bonds in different systems

System	r 1st, Å			τ 1, ns		τ 2, ns
	HS–O	OS–H1	OS–H2	O–H1	O–H2	O–H1	O–H2
METH, S = 1.5	1.88	2.07	2.15	0.0267	0.0239	0.0397	0.0499
ETH, S = 2.0	1.97	2.07	2.07	0.0492	0.0595	0.1113	0.1018
ETH, S = 2.8	1.93	2.17	2.10	0.0558	0.0558	0.1102	0.1102
ACE, S = 1.5	—	2.07	2.05	0.0508	0.0499	0.0181	0.0224
ACE, S = 2.0	—	2.12	2.12	0.0623	0.0646	0.0193	0.0249
ACE, S = 3.0	—	2.17	2.10	0.0646	0.0814	0.0249	0.0222
ACE, S = 4.0	—	2.10	2.12	0.1176	0.1162	0.0222	0.0306
2-BUT, S = 2.0	—	1.90	1.90	0.0832	0.0960	0.4490	0.3830
ANI, S = 2.5	—	2.51	—	0.0755	0.0949	0.2436	0.2206
ANI, S = 3.7	—	2.49	—	0.0915	0.1176	0.2052	0.1876
ANI, S = 5.0	—	2.36	—	0.0497	0.0630	0.2067	0.2494
THF, S = 1.25	—	2.36	2.36	0.0567	0.1291	0.2103	0.1310
DMF, S = 2.0	—	1.95	1.95	1.1237	1.0846	0.1459	0.1785
DMF, S = 2.6	—	1.95	1.95	0.0317	0.0516	0.2902	0.3754
EA, S = 3.5	—	1.90	1.90	0.1247	0.0940	0.0750	0.0952
EA, S = 5.0	—	1.93	1.93	0.0733	0.0905	0.1804	0.1466
EL, M3	2.10	2.10	2.10	2.3455	1.7436	1.4290	1.7268
EL, M4	1.93	2.15	2.17	1.6292	1.7552	1.6289	1.7548
EL, M5	1.93	2.13	2.03	1.6473	1.5627	1.6448	1.5949

---

Fig. 6 shows VMD snapshots of hydrogen bonds between CBZ molecules and solvent–CBZ in different systems at low supersaturations. The hydrogen bonds between the solvent and CBZ molecules exist at d_HB < 2.0 Å in ethanol (Fig. 6a), 2-butanone (Fig. 6c), N,N-dimethylformamide (Fig. 6f) and ethyl acetate (Fig. 6g) and thus, solvent polar atomic sites in these solvents compete with polar sites of carboxamide groups of CBZ, which avoids participation of both oxygen and hydrogen of the carboxamide group in forming a hydrogen bond dimer between CBZ molecules in these solvents. The hydrogen bonds between CBZ and solvents like acetone, tetrahydrofuran and ethyl lactate occur at d_HB > 2.2 Å; as a result, the probability of formation of CBZ hydrogen bond dimers is higher in these systems (Fig. 6b, d and h). There are no hydrogen bonds between the anisole and CBZ molecules as reported in Fig. 6e; however, the presence of interactions between aromatic groups of CBZ and anisole is possible, but is not explored in the present work as systems are modelled using classical force fields and similarly no hydrogen bonds were observed between the tetrahydrofuran and CBZ molecules.

Fig. 6 VMD snapshots of hydrogen bonding distances between carbamazepine molecules in (a) ethanol (S = 2.0), (b) acetone (S = 1.5), (c) 2-butanone (S = 2.0), (d) tetrahydrofuran (S = 1.25), (e) anisole (S = 2.5), (f) N,N-dimethylformamide (S = 2.0), (g) ethyl acetate (S = 3.5) and (h) ethyl lactate (M3).

The lifetime of hydrogen bonds between CBZ molecules could be an important factor that governs the transformation between different types of CBZ aggregates. The hydrogen bond distances between CBZ dimers vary with time, where hydrogen bonds, d_HB > 2.25 Å after 4–5 ps in most of the systems. However, the CBZ molecules in a dimer have a short-memory effect and hydrogen bonds are again reformed (d_HB < 2.25 Å) leading to CBZ dimers. As a result, the lifetime of a hydrogen bonded CBZ dimer in acetone-based systems is ≥2 ns; on the other hand, single hydrogen bonds have lifetime of ≥0.5 ns as shown for the S = 1.5 system in acetone in Fig. S5a. There are a smaller number of single hydrogen bonds and only one dimer in ethyl lactate for the S = 5.0 system, but the lifetime of these hydrogen bonds is longer in ethyl lactate than in acetone as shown in Fig. S5b for the S = 5.0 system in ethyl lactate, where both dimers and single hydrogen bonds have lifetimes of 4 ns. The larger CBZ clusters are mainly formed by CBZ dimers and monomers where the hydrogen bonds between dimers and monomers are dynamic in nature. These dimers and monomers also show a short-memory effect where a cluster is formed, broken, again reformed and is stable for MD steps > 10 ns as shown for aggregates of size, n = 5 in the S = 1.5 system in acetone in Fig. S6. This breaking and reformation of hydrogen bonds and thus, short-memory effect between hydrogen bonded CBZ molecules could be related to the reorientation of the CBZ molecules. Table 3 reports relaxation times, τ₁ and τ₂ of hydrogen bonds, O⋯H₁ and O⋯H₂ between carboxamide groups of CBZ molecules obtained using bi-exponential fit to the second-order reorientational correlation functions in different systems (details are reported in Table S1 and Fig. S7). The relaxation times, τ₁ or τ₂ are similar for both O⋯H₁ and O⋯H₂ hydrogen bonds in different solvents except for THF, DMF and EL. The relaxation times for the O⋯H₁/O⋯H₂ bond between CBZ molecules decreases as EL > DMF > EA > 2-BUT ≈ ANI ≈ THF > ACE > ETH > METH in different solvents, which suggests that the stability of the hydrogen bond between CBZ molecules depends on both solvent–solvent and solvent–CBZ interactions rather than CBZ–CBZ interactions.

Angle distribution, P(θ) between 6-membered rings of a CBZ dimer was calculated based on minimum distance, d_PP observed between 6-membered rings (d_PP = 4.25 Å) in different systems with an exception of d_PP = 5.0 Å in EA-based systems and d_PP = 4.75 Å in M3 and M4 systems of EL-based systems after t = 70 ns as there were no CBZ dimers with d_PP = 4.25 Å in these systems. Fig. 6 reports different types of CBZ dimers observed in different solvents as a function of supersaturation and time (all the angle distribution plots are reported in Fig. S8–S14). There was only one type of CBZ dimer in most of the solvents after 60 ns of NVT simulation – FIII-type CBZ dimers in METH (S = 1.5), THF (S = 1.25), acetone (S = 1.5, 2.0, 3.0 and 4.0), 2-BUT (S = 2.0), ANI (S = 2.5 and 5.0), EA (S = 3.5) and EL (M3 and M5); FII-type CBZ dimers in ETH (S = 2.0), ANI (S = 3.7) and EL (M4) and FI-type CBZ dimers in the DMF (S = 2.6) system. There were mixtures of both FII- and FIII-type CBZ dimers in a few of the systems; DMF (S = 2.0), ETH (S = 2.8) and EA (S = 5.0). FIII-type CBZ dimers were observed in both the low supersaturation systems (THF, S = 1.25 and METH, S = 1.5), which is in agreement with the reported experiments. The type of CBZ dimer observed in a solvent changes with supersaturation in different solvents except for acetone, which is in agreement with the reported experiments.¹⁹ FIII-type CBZ dimers are present at both low and high supersaturations in anisole (S = 2.5 and 5.0) and ethyl lactate (M3 and M5) but the FII-type CBZ dimers exist for the intermediate supersaturations (ANI, S = 3.7 and EL, S = M5) in these systems.

Transformation of CBZ dimers to other forms was further evaluated by simulating different systems (ethanol, 2-butanone, acetone, anisole, ethyl acetate and ethyl lactate) for another 10 ns using the NVT ensemble as reported in Fig. 7. The FIII-type dimers were present only in the S = 2.0 and 3.0 systems in acetone and M5 system in ethyl lactate with an increase in time. There are 20% FII-type CBZ dimers and 80% FIII-type dimers in both the S = 1.5 and 4.0 systems in acetone after 70 ns as shown in Fig. S8. There were only FII-type CBZ dimers observed in the ETH (S = 2.0), which is consistent with reported experiments,¹⁵ ANI (S = 3.7) and EL (M4) systems with an increase in time. There is transformation of the CBZ dimers to other types with time in a few of the systems; FIII- to FI-type CBZ dimers in the S = 3.5 system in EA, a mixture of FIII-type CBZ dimers in the S = 2.8 system in ETH, FIII- to FII-type CBZ dimers at S = 3.5 in EA and S = 2.5 in ANI systems and FIII-type to a mixture of FIII-type and FI-type CBZ dimers in 2-butanone.¹⁵

Fig. 7 CBZ dimers of FI-, FII-, and FIII-type as a function of supersaturation (S) in different solvents after t = 60 ns and 70 ns of NVT simulations; also shown are mixtures of different forms as M; lines are shown as a guide to show different supersaturation systems with the same solvent system.

Hydrogen bonding between CBZ molecules and CBZ–solvent molecules is one of the main factors that affects aggregation of CBZ molecules. The presence of weak CBZ–acetone interactions and HB dimers between CBZ molecules in acetone causes overlap between the 6-M rings of the CBZ molecules that lead to formation of FIII-type CBZ dimers. However, at low (S = 1.5) and high (S = 4.0) supersaturations in acetone, other types of hydrogen bonding were observed – single HB, HB dimers and HB trimers that lead to a mixture of FII-type and FIII-type CBZ dimers. On the other hand, intermediate supersaturations (S = 2.0 and 3.0) in acetone show only HB dimers and thus, FII-type CBZ dimers were observed in these systems. The presence of favourable hydrogen bonding between carboxamide atoms of CBZ with oxygen and polar hydrogen of solvent molecules reduces the probability of HB dimers and thus, single hydrogen bonds were observed in ethanol (Fig. S15a and b). This provides orientation flexibility among CBZ molecules due to which a mixture of CBZ dimers or FII-type CBZ dimers were observed for different supersaturation (S = 2.0 and 2.8) systems in ethanol. There were a smaller number of HBs between CBZ molecules (one hydrogen bond) observed at very low supersaturation in methanol (Fig. S15c) as both oxygen and hydrogens of the carboxamide group of CBZ form hydrogen bonds with methanol, due to which FIII-type CBZ dimers were observed in the S = 1.5 system in methanol. The CBZ–CBZ hydrogen bonding is more favourable in in low supersaturation (S = 1.25) systems in tetrahydrofuran and mainly HB dimers and trimers (Fig. S15d) are present among the CBZ molecules. The tetrahydrofuran–CBZ interactions are very weak, and thus only FIII-type CBZ dimers were observed in tetrahydrofuran. Fig. S16 shows that there are single HBs, HB dimers and trimers in solvents like 2-butanone, anisole and ethyl acetate and solvent–CBZ interactions are also favourable; as a result, a mixture of types of CBZ dimers was observed in these systems both as a function of supersaturation and time. There are no CBZ–CBZ hydrogen bonds in N,N-dimethylformamide but a HB dimer was observed in the S = 2.6 system (Fig. S16a and b). The solvent–CBZ interactions are highly favourable in N,N-dimethylformamide and 2-butanone unlike methanol, due to which a mixture of CBZ dimers were observed in the low supersaturation, S = 2.0 system, and FI-type CBZ dimers were observed in the S = 2.6 system in N,N-dimethylformamide and S = 2.0 in 2-butanone. Interestingly, only HB dimers were observed in low (S = 2.5) and high (S = 5.0) supersaturation systems in ethyl lactate but both single HBs and HB dimers were observed in the intermediate (S = 3.7) system in this solvent as shown in Fig. S16d–f due to which FIII-type was observed in the S = 2.5 and 5.0 systems but FII-type CBZ dimers were observed in the S = 3.7 system. However, single HBs were observed in many systems with an increase in time due to which a mixture of CBZ dimers was observed in these systems. The angle distribution analysis was further performed at 298 K as reported in Table S2, Fig. S17 and S18 and the results were similar to the distributions observed at 265 K except for FI-type dimers observed in the S = 2.0 system in DMF and FII-type CBZ dimers in the S = 5.0 system in ethyl lactate.

Conclusions

Aggregation behaviour of carbamazepine molecules in different solvents was studied using molecular dynamics simulation studies at 265 K and 298 K. The dynamics of CBZ molecules is governed mainly by the solvent–solvent interactions whereas the size of the CBZ aggregates is governed mainly by CBZ–solvent interactions. The FI/FII/FIII-type CBZ dimers could be correlated as growth synthons with the similar structural synthons in different polymorphs of CBZ. The formation of these growth synthons is mainly governed by the type of hydrogen bonding between CBZ molecules and between CBZ–solvent molecules. Future studies will investigate the role of these growth synthons as possible collective variables in governing the formation of different polymorphs of CBZ by calculating the free energy of different CBZ polymorphs using rare event simulation techniques. There are more hydrogen bonded CBZ dimers when solvent–CBZ interactions are weak as observed in acetone, which leads to the formation of FIII-type CBZ dimers both at 265 K and 298 K. Similarly, hydrogen bonded CBZ dimers with a longer lifetime exist in ethyl lactate (high supersaturation system) that also favour the formation of FIII-type CBZ dimers. The single hydrogen bonds and presence of favourable solvent–CBZ interactions favour FII-type CBZ dimers as was observed in ethanol, whereas systems where solvent–CBZ interactions dominate favour FI-type CBZ dimers as observed in N,N-dimethylformamide. The formation of FIII-type CBZ dimers occurs in solvents that are weak hydrogen bond acceptors, such as acetone. However, if the solvent is both a hydrogen bond acceptor and donor but has CBZ hydrogen bonded dimers, it could also facilitate nucleation of FIII-type CBZ dimers at high supersaturation as observed in ethyl lactate. Ethyl lactate is a green solvent and pharmaceutically relevant solvent and further experimental studies are needed to assess if ethyl lactate could be a potential candidate for selective polymorphism of FIII of carbamazepine at high supersaturations.

Conflicts of interest

The authors declare no conflict of interest.

Data availability

Input files for LAMMPS package will be shared upon request to the authors.

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp03059e.

Acknowledgements

MS thanks DST-SERB for funding (Grant No. ECR/2017/002937) and CMSD, University of Hyderabad for computational facilities. SS thanks CSIR-JRF for the fellowship for the PhD Program.

References

1. N. Rodríguez-Hornedo and D. Murphy, J. Pharm. Sci., 1999, 88, 651–660.
2. S. R. Byrn, Annu. Rep. Med. Chem., 1985, 20, 287–294.
3. J. M. Miller, B. M. Collman, L. R. Greene, D. J. W. Grant and A. C. Blackburn, Pharm. Dev. Technol., 2005, 10, 291–297.
4. M. A. O'Mahony, A. Maher, D. M. Croker, Å. C. Rasmuson and B. K. Hodnett, Cryst. Growth Des., 2012, 12, 1925–1932.
5. R. J. Davey, N. Allen, N. Blagden, W. I. Cross, H. F. Lieberman, M. J. Quayle, S. Righini, L. Seton and G. J. T. Tiddy, CrystEngComm, 2002, 4, 257–264.
6. B. Y. Shekunov and P. York, J. Cryst. Growth, 2000, 211, 122–136.
7. L. Yu, J. Am. Chem. Soc., 2003, 125, 6380–6381.
8. C.-H. Gu, V. Young Jr. and D. J. W. Grant, J. Pharm. Sci., 2001, 90, 1878–1890.
9. S. L. Childs, P. A. Woods, N. R. Hornedo, L. S. Reddy and K. I. Hardcastle, Crys. Growth Des., 2009, 9, 1869–1888.
10. A. J. Cruz Cabeza, G. M. Day, W. D. S. Mortherwell and W. Jones, Chem. Commun., 2007, 1600–1602.
11. M. Lang, J. W. Kampf and A. J. Matzger, J Pharm Sci, 2002, 91, 1186–1190.
12. J.-B. Arlin, L. S. Price, S. L. Price and A. J. Florence, Colloids Surf., 2008, 63, 153–158.
13. A. L. Grzesiak, M. Lang, K. Kim and A. J. Matzger, J. Pharm. Sci., 2003, 92, 2260–2271.
14. R. J. Behme and D. J. Brooke, J. Pharm. Sci., 1991, 80, 986–990.
15. R. C. Kelly and N. R. Hornedo, Org. Process Res. Dev., 2009, 13, 1291–1300.
16. E. A. Losev, D. Ya Zheltikova, V. A. Drebushchak and E. V. Boldyreva, CrystEngComm, 2024, 26, 6796–6804.
17. A. Getsoian, R. M. Lodaya and A. C. Blackburn, Int. J. Pharm., 2008, 348, 3–9.
18. D. Ying, W. Kristen, H. Matthew, W. Mahlet, D. Patrick, M. Allan, H. T. Alan and T. Bernhardt, J. Am. Chem. Soc., 2012, 134, 673–684.
19. J. Ouyang, J. Chen, I. Rosbottom, W. Chen, M. Guo and J. Y. Y. Heng, CrystEngComm, 2021, 23, 813–823.
20. A. S. Meyerson and B. L. Trout, Science, 2013, 341, 855–856.
21. P. G. Vekilov, Nanoscale, 2010, 2, 2346–2357.
22. D. Kashchiev, J. Cryst. Growth, 2020, 530, 125300.
23. F. Ito, Y. Suzuki, J. Ichi Fujimori, T. Sagawa, M. Hara, T. Seki, R. Yasukuni and M. L. de la Chapelle, Sci. Rep., 2016, 6, 22918.
24. S. Kakkar, S. Bhattacharya, P. Andre Cazade, D. Thompson and A. Rasmuson, Cryst. Growth Des., 2024, 24, 5740–5753.
25. D. S. Reddy, D. C. Craig and G. R. Desiraju, J. Am. Chem. Soc., 1996, 118, 4090–4093.
26. G. R. Desiraju, J. Am. Chem. Soc., 2013, 135, 9952–9967.
27. J. Rezac, K. E. Riley and P. Hobza, J. Chem. Theory Comput., 2011, 7, 2427–2438.
28. J. D. Dunitz and A. Gavezzotti, Cryst. Growth Des., 2012, 12, 5873–5877.
29. N. A. Gracia, R. I. Malini, C. L. Freeman, R. Demichelis, P. Raiteri, N. A. J. M. Sommerdijk, J. H. Harding and J. D. Gale, Cryst. Growth Des., 2019, 19, 6422–6430.
30. M. Salvalaglio, C. Perego, F. Giberti, M. Mazzotti and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2014, 112, E6–E14.
31. M. Sharma and B. L. Trout, J. Phys. Chem. B, 2015, 119, 8135–8145.
32. W. Liu, L. Dang, S. Black and H. Wei, J. Chem. Eng. Data, 2008, 53, 2204–2206.
33. H. Sun, J. Phys. Chem. B, 1998, 102, 7338–7364.
34. S. Plimpton, J. Comput. Phys., 1995, 117, 1–42.
35. I. Rosbottom, T. N. H. Cheng and J. Y. Y. Heng, Chem. Eng. Technol., 2020, 43, 1152–1159.
36. A. M. S.-Devlin, J. P. Petersen, J. Liu, G. A. Turner, J. C. Poutsma and S. A. Kandel, J. Phys. Chem. C, 2020, 124, 5213–5219.

---