Structural Dynamics of Carbon Dots in Water and N, N-Dimethylformamide Probed by All-Atom Molecular Dynamics Simulations.

Carbon dots (CDs), one of the youngest members of the carbon nanostructure family, are now widely experimentally studied for their tunable fluorescence properties, bleaching resistance, and biocompatibility. Their interaction with biomolecular systems has also been explored experimentally. However, many atomistic details still remain unresolved. Molecular dynamics (MD) simulations enabling atomistic and femtosecond resolutions simultaneously are a well-established tool of computational chemistry which can provide useful insights into investigated systems. Here we present a full procedure for performing MD simulations of CDs. We developed a builder for generating CDs of a desired size and with various oxygen-containing surface functional groups. Further, we analyzed the behavior of various CDs differing in size, surface functional groups, and degrees of functionalization by MD simulations. These simulations showed that surface functionalized CDs are stable in a water environment through the formation of an extensive hydrogen bonding network. We also analyzed the internal dynamics of individual layers of CDs and evaluated the role of surface functional groups on CD stability. We observed that carboxyl groups interconnected the neighboring layers and decreased the rate of internal rotations. Further, we monitored changes in the CD shape caused by an excess of charged carboxyl groups or carbonyl groups. In addition to simulations in water, we analyzed the behavior of CDs in the organic solvent DMF, which decreased the stability of pure CDs but increased the level of interlayer hydrogen bonding. We believe that the developed protocol, builder, and parameters will facilitate future studies addressing various aspects of structural features of CDs and nanocomposites containing CDs.

Introduction

Although carbon dots[1] (CDs) were discovered rather recently, they have become one of the most widely studied classes of nanomaterials. They represent a perspective material for diagnostic and therapeutic applications (theranostics), especially in optical and photoacoustic imaging, drug delivery, and photothermal therapy.[2,3] CDs have a very low toxicity compared to traditional quantum dots based on metal chalcogenides but possess similar fluorescence advantages, such as low photobleaching.[2−4] Therefore, the field of CD studies has been growing rapidly over the last four years, boosted by expectations of breakthrough applications even in technological fields, e.g., light emitting diodes, water splitting, etc. CDs are quasi-spherical particles with a multilayer graphene structure and sizes below 10 nm in all dimensions (typically 2–3 nm).[5,6] A similar system, graphene dots (GDs), consists of few-layer graphene with sizes up to 100 nm.[5,7,8] Both can be modified either in the core structure by partial substitution of carbon with other elements (i.e., doping by nitrogen, sulfur, or boron) or by surface functionalization. The surface shell contains various functional groups (carbonyl, carboxylic, amine, amide) and can bear significant net charge that is recognizable by enzymes[4] and modifies significantly the fate of CDs in organisms. However, the exact mechanism of CD behavior in solution or biosystems is still not well explained and remains to be clarified.

Although other carbon nanostructures have been extensively studied both theoretically and experimentally, no full model of CDs is available for all atom molecular dynamics (MD) simulations. MD simulations may provide useful and very detailed information (reaching femtosecond time and atomic space resolutions simultaneously) on the interaction of CDs or GDs with other systems, such as bioenvironments or other materials. Theoretical calculations in this field were initially a byproduct of searching for a smaller graphene model[9] to explain the fluorescent properties. Most of the simulations of either CDs or GDs used single-layer graphene flakes (containing ca. 30–50 carbons) studied by quantum-chemical tools.[4,10−12] However, evaluation of CD and GD behavior in complex systems and on larger scales requires more approximate approaches based on, e.g., molecular mechanics. In several MD studies, the behavior of graphene flakes was investigated by focusing on a few graphene layers[13−15] or a graphene flake in water.[16] Theoretical studies were also extended to interactions of graphene flakes with biological systems, such as DNA[17] or lipid membranes.[18−22] Recently, an enantioselective pore in a set of hexagonal graphene sheets was investigated by MD.[23] Nevertheless, to the best of our knowledge, no MD studies of a spherical CD have yet been published.

Here, we present a full procedure for performing CD simulations. We prepared a GUI of a builder provided as a plug-in for the widely used VMD software.[24] We derived parameters for several surface functional groups, i.e., hydroxyls, carbonyls, and protonated and unprotonated carboxyls. Using various models differing in size and surface functionalization, we analyzed the behavior of CDs in water and N,N-dimethylformamide (DMF). We analyzed the geometry and stability of the resulting CDs, identified internal motions of individual layers, and studied differences in the network of interlayer hydrogen bonds in the surrounding solvent. We believe that the provided builder and parameters will aid future MD simulations and atomistic understanding in this novel and rapidly developing area.

Methods

VMD Builder

Based on experimental observations of CDs,[2] the builder uses a hexagonal graphene-like sheet as a basic shape, for which the user can set the edge size in units of number of benzene rings. The size of the CD layers gradually decreases to generate a spherical shape. The user can either choose a level of edge coverage by a chosen functional group and the builder places them randomly or the positions of functional groups can be assigned manually. These approaches can be combined, and an automatically generated CD can then be manually edited and groups can be added or deleted. In this way, a CD with a mixture of functional groups can be prepared with any composition. Future development may focus on enabling doping of the CD core. Using the functionalities of Topotools,[25] it is also possible to save both the structure and the GROMACS topology, i.e., a “fake” one that requires atom type and charge adjustment (see later). These utilities are implemented in the CD VMD builder (http://cd-builder.upol.cz), which was built upon the graphene and nanotube builder[26] already implemented in VMD 1.9.3.[24] A detailed description of the builder algorithm can be found in the Supporting Information.

Charges and Force Field

For assigning the partial charges, we used circumcoronene models as they were small enough to allow quantum-chemical calculation and simultaneously large enough to act as polycyclic aromatic systems. First, all (functionalized, see later) circumcoronene molecules were fully optimized using the Becke three-parameter hybrid density functional B3LYP and 6-31++G(d,p) basis set as we used earlier.[27] Assignment of the partial charges on individual atoms was based on fitting the electrostatic potential calculated on different circumcoronene models using the CHELPG[28] approach at the HF/cc-pVDZ level of theory in a vacuum. For compatibility with available force fields, we calculated the partial charges by HF/6-31G* (for AMBER99SB) and B3LYP/cc-pVTZ at ε = 4 (for AMBER03) and present the results in the Supporting Information.

For each functional group, we prepared several circumcoronene molecules with various positions and number of functional groups. These were further sorted according to their position (in zigzag or armchair conformation), and the partial charges assigned to the atoms of various types and on atoms close to the functional groups (see, e.g., Figure S5) were averaged. However, this level of simplification is not generally relevant because, in a conjugated system, a local change of chemistry can cause a very distant change in the electrostatic potential and calculated partial charges. Thus, this issue needs to be considered by the user for each specific case. Indeed, for simulation of a single layer of coronene or circumcoronene size, a separate partial charge calculation would be appropriate. We aimed to make the parametrization usable in large systems that cannot currently be investigated at an adequate level of theory. We encourage the readers to study the Supporting Information, which contains detailed information about the partial charge assignment procedure, differences between the smaller (coronene) and larger (circumcoronene) CD layer models, the nomenclature used for the carbons in the layers and the individual molecules with calculated partial charges used for the final charge estimations (Figure S7).

Multiple force field parameters are available for benzene, graphite, or graphene simulations.[29−32] Here, we utilized the OPLS all-atom force field[33] with refinements on carbon nonbonded Lennard-Jones parameters proposed by Cheng and Steele,[32] which have been successfully used to study adsorption of small molecules on graphene.[34] Other atom types in functional groups were chosen according to the local chemistry from a regular OPLS all-atom[33] force field (Table ). The output from VMD provided by TopoTools[25] was modified in order to run the simulations—the initial pdb file was sorted to separate each layer and keep the topology ordering. We included proper atom types and charges into the topologies, but dihedral angles were modified to maintain a layer plane (the bash scripts for assigning proper partial charges etc. can be found in the Supporting Information). After such postprocessing, the system consisted of several residues, each of them representing a single CD layer. Our parametrization strategy allows future simulations of complex hybrid systems as we provide the model also in AMBER99SB,[35] can be used also in later force fields, and can allow simulations of CDs in bioenvironments.

Table 1

Atomic Properties Showing for Each Functional Group the Atom Names Assigned by the VMD Builder, Atom Types Used in the Topology,a and the Charges on the Atoms in Armchair and Zigzag Conformations

group/atom	VMD name	atom type	atom charge—armchair	atom charge—zigzag	nearby carbon	next-ring carbon
edge C	CA	Cheng and Steele[32]	–0.210
edge H	HA	OPLS 146	0.179
pure CD	CA	Cheng and Steele[32]	–0.180	–0.400
H	HA	OPLS 146	0.120	0.180
hydroxyl	CA	OPLS 166	0.374	0.153	–0.452	–0.292
O	OH	OPLS 167	–0.605	–0.536
H	HO	OPLS 168	0.400	0.400
carbonyl	CA	OPLS 320	0.705	0.560	–0.371	–0.27
O	ON	OPLS 340	–0.580	–0.540
carboxyl	CA	Cheng and Steele[32]	–0.106	–0.308	–0.223	–0.214
C	CX	OPLS 267	0.766	0.789
O	OX	OPLS 269	–0.610	–0.595
O	OC	OPLS 268	–0.610	–0.595
H	HX	OPLS 270	0.427	0.418
carboxyl—charged	CA	Cheng and Steele[32]	–0.150	–0.222	–0.283	–0.231
C	CR	OPLS 271	1.017	1.096
O	OR	OPLS 272	–0.901	–0.934
O	OK	OPLS 272	–0.901	–0.934

CA atoms in each groups are edge carbon atoms connected to the functional group.

MD Simulations

Using the CD builder, we prepared dots with 3–10 benzene rings (and one additional CD with 18 benzene rings) on the edges for MD simulations with Gromacs 5.1.[36] The CDs were energy minimized in vacuum using the following setup: steepest descent method, cutoff of 1.0 nm for coulomb and van der Waals interactions, Coulombic interactions above the cutoff calculated by the particle mesh Ewald (PME) method,[37] and energy minimization until the total energy difference was less than 10 kJ/mol. Afterward, the dots were solvated (for shorter 50 ns simulations, 1.2 nm of water from the CD was added in each direction; we recommend, as used for longer 1 μs simulations, 2.0 nm of water in each direction; for an impact, see the Results and Discussion section) with the TIP3P water model.[38] When we simulated CDs with charged carboxyls, we also added Na+ cations[39] to neutralize the system and used a physiological concentration of 0.15 M of Na+ cations[39] and Cl– anions.[40] We minimized the system with the same setup as in vacuum, and a MD production run was performed (2 fs time step, coulomb and van der Waals interactions calculated explicitly up to 1.0 nm, and further Coulombic interactions calculated by PME;[37] for van der Waals interactions we used a cutoff scheme; the temperature was maintained by the V-rescale thermostat[41] to 300 K (CD and solvent coupled separately); the Berendsen barostat[42] was applied isotropically to 1 bar; and bonds including hydrogens were constrained by LINCS[43]). The used water and ion model was compatible with the OPLS all-atom force field. We initially chose CDs with 6 benzene rings along the middle edge (2.1 nm gyration diameter) and 30% coverage of individual surface functional groups or a pure dot and simulated these for 1 μs. The structural features and hydration of the CDs converged after 20 ns (Figure S17), indicating that important trends in the size, shape, number of hydrogen bonds, etc., could be evaluated on a time scale of several tens of nanoseconds. Therefore, we simulated all other systems for only 50 ns each. However, long-term internal motion (layer rotation and, consequently, intraparticle hydrogen bond formation) converged on >100 ns time scale, and though some trends were noted, these need to be taken with care. In addition to simulations in water, we performed simulations of CDs in DMF with the same setup as for water. The list of performed simulations can be found in Table S2, and in the Supporting Information we also provide the mdp file with the simulation protocol.

Results and Discussion

Behavior of 2.1 nm CDs in Water

First, we investigated the shape of a CD with ∼2.1 nm diameter of gyration without any surface modification, i.e., terminated only by hydrogen atoms. During the 1 μs long simulation in water, the CD remained stable in a sphere-like shape (Figure ), which was slightly smaller in height than in width (Table ) with radius of gyration of 1.05 ± 0.01 nm. We used the gyration diameter calculated from the radius of gyration in Table instead of the middle layer diameter as a measure of the CD size as it is a more general descriptor of CD size. The mean interlayer distance calculated as 0.34 nm was in accord with experimentally observed values in graphite[44] or CDs.[6,45,46] The individual layers stayed generally flat, and their undulations (calculated as widths of CD density peaks) were lower than 0.15 nm (Figures S5 and S6) and decreased with increasing layer size. We observed that the horizontal positions of individual layers fluctuated and shifted from the middle of the sphere more with decreasing layer size (and increasing distance from the CD center), but all layers remained on average below 0.2 nm from the middle position (Figure S9). We did not observe a strong hydrophobic cage of water molecules around the CD layers (the distance between surface hydrogens and the closest water oxygens was around 0.30 nm without any significant water density increase). On the other hand, we observed a hydrophobic gap and increased water density above the edge layers at approximately the same distance as the layer spacing, i.e., 0.34 nm from the carbon plane (Figure ).

Figure 1

Density plot (left) and side view (right) of a pure CD with diameter of gyration of 2.1 nm (black) in water (red). The layer undulations are reflected in the widths of the density peaks of individual layers. Water molecules surround the CD with a hydrophobic gap above and below the graphene-like planes but do not penetrate inside the CD structure. The structure of the CD is a snapshot from MD simulation displayed as sticks with cyan carbons, white hydrogens, and red oxygens (for clarity only water molecules within 0.5 nm are displayed).

Table 2

Geometry of Built CDsa and Number of Atoms in the Respective Hydrated Systemsb

	no. of layers
edge size (no. of benzene rings)	above middle layer	total	diameter of the central layer (nm)	dot height (nm)	Rgyr	no. of atoms, 103
3	2	5	1.6	1.74	0.63	6
4	2	5	2.0	1.74	0.74	7
5	3	7	2.4	2.34	0.87	9
6	3	7	2.9	2.34	1.05	12
7	4	9	3.3	3.04	1.18	14
8	4	9	3.8	3.04	1.35	20
9	5	11	4.2	3.74	1.48	23
10	5	11	4.7	3.74	1.65	29
18	9	19	8.3	6.44	2.87	99

Dot height was calculated from the density plot; therefore, we added 0.34 nm as twice the carbon van der Waals radius, and Rgyr stands for the radius of gyration.

Some of the CDs are displayed above the table as licorice structures labeled with their edge sizes.

Synthesized CDs bear surface functional groups that typically contain oxygen atoms. Therefore, we analyzed the behavior of CDs containing hydroxyl, carbonyl, and protonated (neutral) and unprotonated (negatively charged) carboxyl groups. The presence of surface functional groups affected the shape of CDs and individual layers. Layers of the 2.1 nm wide CDs with 30% edge coverage stayed generally flat with slight undulations. The smallest layer undulations were observed for the hydroxyl modified CD (0.07 nm in the middle layer), followed by the uncharged carboxyl and carbonyl modified CDs (0.09 and 0.10 nm, respectively). The largest undulation was observed when the CD contained charged carboxyls (0.14). It should be noted that, in all cases, the undulations of the middle layer were lower than 0.2 nm (Figure S8). Similarly to the CD without surface modifications, the undulations decreased with increasing layer size in most of the cases. However, in the case of charged carboxyls, the undulations increased while the layer size decreased on one side of the CD. The presence of charged carboxyls and carbonyls also affected the horizontal shifts of individual layers, which fluctuated significantly around their initial position (Figures S8 and S9), increasing their radii of gyration to 1.27 and 1.14 nm, respectively (vs 1.05 nm for the pure CD).

The surface groups also influenced the intramolecular kinetics. The mutual orientation of neighboring layers was able to be evaluated with a periodicity of 60° due to their hexagonal shape. We observed that the middle layers of a pure CD was oriented by 30° ± k × 60° to each other (Figure ) and rotated with mean time spent in each individual orientation k lower than 2 ns. When surface groups were added, the rotation rates significantly decreased. The mean time spent in each orientation increased to 6 ns for CDs with hydroxyls and carbonyls, whereas in the case of both charged and uncharged carboxyl groups, we were not able to calculate the rotation rate and the layers stayed mostly ±15° from the initial conformation. It should be noted that the smaller, more distant layers rotated more frequently (Figure S10). Because of these rotations, we could not observe stable AB stacking of the layers. However, the distribution of carbon atoms relative to each other (excluding bonded interactions) showed a probable distance between two carbons of 0.37 nm (Figure S11), which corresponds to AB stacking.

Figure 2

Rotation of the nearest middle layers of CDs (left) and distribution of the positions (right). The inset of two CD layers shows the monitored angle.

Relaxation of the structure, i.e., rotation of the layers, also increased the number of interlayer hydrogen bonds, particularly when carboxyl groups covered the edge (Figure ). In the case of hydroxyls, their ability to form hydrogen bonds was limited owing to the short length of hydroxyl groups, which did not span the 0.34 nm interlayer distance. The lifetime of interlayer hydrogen bonds differed significantly. For hydroxyl groups, the lifetime was ∼2 ps compared with ∼50 ns for carboxyls. Whenever possible, CDs formed hydrogen bonds with water molecules (Figure S12). After relaxation, the average number of hydrogen bonds with water generally remained constant. The lifetimes of hydrogen bonds with water were significantly shorter than the interlayer ones in the case of uncharged carboxyls (173 ps vs 349 ps for charged carboxyls). In the case of hydroxyls and carbonyls, the hydrogen bond lifetimes were 96 and 66 ps, respectively.

Figure 3

Number of interlayer hydrogen bonds in CDs covered with protonated (uncharged) carboxyls in water (magenta) and dimethylformamide (DMF, orange). The CD structure relaxation lasted for at least 100 ns. In the right panel, we show a snapshot of part of the CD surface with highlighted hydrogen bonds between water molecules (red) and carboxyl groups (blue). The CD is shown as gray licorice, and highlighted atoms are depicted with colors: carbon in cyan, hydrogen in white, and oxygen in red.

Simulations in DMF

Carbon nanoparticles can be prepared from larger carbon structures by, e.g., sonication in DMF.[47] Therefore, we expected degradation or a generally lower stability of CDs in this and other organic solvents. To analyze the behavior of CDs in such an organic solvent, we performed simulations in DMF representing a polar organic solvent (with relative permittivity of 36.7 under ambient conditions). CDs (2.1 nm, 30% surface coverage of functional groups, and also pure CDs with edge size varying from 3 to 18 benzene rings on the edge, see Table S2) in DMF exhibited a constant interlayer distance of 0.34 nm and adopted a spherical shape, resembling the behavior in water. On the other hand, the smallest dots (1.3 nm in diameter, nonfunctionalized surface) were not stable because we observed significant sliding of individual layers of the CDs on each other, which finally led to dislocation of the outer layer of coronene size. This layer became solvated by DMF and stayed dislocated from its original position on the CD, finally becoming dissolved in DMF (Figure S13). This could explain why CDs prepared in DMF are generally larger than in water.[48] The larger CDs (>1.5 nm in diameter) were stable, but we observed that the rate of rotations was slightly quicker in DMF than in water (mean times in conformation for 2.1 nm CD were 880 and 344 ps in water and DMF, respectively; see Figure S14). In functionalized CDs, we observed a higher amount of interlayer hydrogen bonds compared to the corresponding simulations in water, especially in the case of uncharged carboxyls (Figure ). The lack of hydrogen bonds with water in the organic aprotic solvent therefore seemed to be compensated by the larger number of interlayer hydrogen bonds stabilizing the CD. Therefore, the functionalized CDs with interlayer hydrogen bonds were more stable in DMF than in water, but the pure CDs were less stable in DMF. Another consequence of the absence of hydrogen bonds (with water) in aprotic solvent DMF may be the aggregation of graphene dots observed in DMF experimentally.[49] These observations are in accord with experiments[48,49] and show the reliability of the presented model for CDs.

Simulations of CDs with Various Sizes and Surface Coverage with Oxygen-Containing Groups

To evaluate the role of the amount of oxygen-containing groups on the CD surface, we generated CDs with different levels of surface coverage and diameters of gyration from 1.3 to 5.7 nm and analyzed their behavior during 50 ns simulations. Taking into account the edge functionalization, the dot size limited the possible fraction of oxygen in the CDs. Coverage of ∼45% of possible surface groups resulted in up to 27 wt % of oxygen in the smallest CD covered with carboxyl groups, but the possible oxygen fraction decreased rapidly with increasing dot size up to ∼7% in the case of the largest simulated dot with 5.7 nm diameter (Figure S15). Carboxyl groups provided the highest fraction of oxygen in the CDs as they possess two oxygens instead of the single oxygen atoms in carbonyl or hydroxyl groups. Experimentally, a significant fraction of oxygen (up to 30 wt % in a set of CDs with 3–15 nm in diameter[50]) has been observed. Therefore, we expected a high fraction of the surface to be covered by functional groups in these dots.

The dot size and edge coverage affected the shape of the resulting nanoparticle. The undulations of individual layers increased with increasing size of the dots. In the case of carbonyl and carboxyls (both charged and uncharged), the undulations also increased with increasing coverage of these functional groups (for detailed data, see Table S3), whereas coverage of the dot with hydroxyl groups decreased the undulations slightly. In the case of carbonyls and charged carboxyls, increasing the coverage also increased the horizontal shift of the layers—in the smallest simulation boxes, carbonyls even formed periodic particles (Figure S16).

Whereas the neutral CDs were stable, the behavior of the negatively charged CDs depended on the charge. CDs bearing 15% of negatively charged carboxyl groups were stable, whereas more charged smaller CDs decomposed in the water phase in some cases. From our simulations, it was shown that CDs can have up to ∼35% of edges covered with charged carboxyls and remain stable (except for the smallest 1.3 nm wide CD, Table S3), but a higher degree of surface coverage leads to instability. Smaller CDs (up to 2.4 nm wide) with higher coverage (∼50% of edges) appeared to be unstable and exfoliated into individual layers or parts. In some cases, the CDs dissolved to two particles separated by ions and water (Figure ). In the largest CDs, the layers stayed stable even with a high negative charge but had a slightly deformed shape. Experimentally, a significant negative charge and proportion of carboxyl groups has been observed.[51] Therefore, their presence on CDs is indubitable. The deformation and destabilization of CDs by carboxyls and carbonyls need to be balanced by the hydrophobic effect of the carbon layers or the effects of other functional groups. However, because of the small range of surface coverage used here, we cannot determine the exact ratio between charge and dot size needed to keep the CD stable. However, we can conclude that negatively charged CDs are stabilized by groups forming hydrogen bonds or other interlayer connections, e.g., uncharged carboxyls, hydroxyls on alkyl chains, positively charged amine groups, etc. The network of interlayer covalent or noncovalent bonds affects not only the stability of the CD but also its shape and intraparticle kinetics.

Figure 4

CDs with deprotonated carboxyl groups (sodium counterions shown as dots): (a) large CDs can stay as a deformed dot, (b) layers may move further from each other with ions connecting the edge carboxyl (−1) groups between the layers, (c) dots may become separated into parts that are either fully separated or reconnect, and (d) individual layers (in the case of small layers) may leave the CDs fully and stay solvated in solution.

The layer rotation rate decreased with increasing size of the dot and depended on the type of functional group. In very small CDs (1.3 nm diameter), the middle layers rotated usually with a mean time spent in the position (separated by 60°) of the order of hundreds of picoseconds (notice the difference from the pure CD, which rotated with a mean time of the order of tens of picoseconds, Table S3). With increasing dot size, this time increased steeply (Figure S14). In the case of carbonyl covered CDs, pure CDs, or CDs with hydroxyls, we were able to calculate the rotation time up to a diameter of 2.7 or 3.0 nm. In other cases, the rotation rate was too slow to estimate the rotation time. Generally, functionalization slowed down the rotation, but without clear dependence on the degree of functionalization. It should be noted that, for a mean time between rotations of the order of tens of nanoseconds, the results should be interpreted with care owing to the errors and limited total time of simulations of 50 ns. Apart from the middle layers, the other layers in the CDs (further from the middle) also rotated relative to each other at a higher rate (data not shown)—the rotation rate clearly depended on the layer size.

The edge groups of the CDs interacted significantly with the surrounding solvent through a network of hydrogen bonds. The mean lifetime of hydrogen bonds of water with hydroxyl or carbonyl groups fluctuated between 1 and 4 ps. In contrast, uncharged carboxyls formed hydrogen bonds with lifetimes of up to 10 ps and charged carboxyls formed hydrogen bonds with lifetimes of up to 100 ps. Uncharged carboxyls also formed interlayer hydrogen bonds, and during the simulations, the CDs gradually increased their amount by rotating the layers. The lifetime of interlayer hydrogen bonds was significantly longer in DMF than in water (of the order of hundreds of picoseconds up to nanoseconds), which may have slowed down significantly the rotation rate of the CD layers. Addition of protonated carboxyls or insertion of, e.g., a methanoyl group that could reach the nearby layer would decrease the rotation rate significantly. Changing the position of differently composed domains or layer dipole orientation may affect other CD properties, such as fluorescence. The effect of interlayer hydrogen bonds may even be increased in environments unable to support hydrogen bonds between the CD and the solvent, such as DMF, where hydrogen bonds are formed between the layers only. These hydrogen bonds may be crucial for the particle stability.

Conclusion

Here, we presented a full procedure for all-atom MD simulations of CDs from structure preparation up to MD simulation. We derived partial charges for carbons and various oxygen-containing groups which are compatible with AMBER or OPLS force fields, allowing future simulations of complex molecular systems containing CDs. We developed an intuitive VMD interface for building CD models and generating inputs for MD simulations. Using MD simulations, we analyzed the behavior of CDs differing in size (from 1.3 to 5.7 nm in diameter), functional groups, and surface coverage. The results demonstrated that the suggested parameters yielded behavior comparable to experiments leading to stable spherical CDs in water, but we also observed a destabilizing effect in the presence of excess carbonyl groups or negatively charged carboxyls. The conducted MD simulations provided a detailed insight into the intraparticle dynamics, whereby individual CD layers rotate with respect to each other. The layer rotation rate decreased with increasing size of the layers and number of functional groups, especially those forming hydrogen bonds. Finally, we examined the stabilization effects of interlayer hydrogen bonds in the aprotic solvent DMF. We found that, in this solvent, the lack of hydrogen bonds formed with the solvent compared to simulations in water was balanced by a surplus of interlayer hydrogen bonds, which also explains the preference for aggregation in DMF. Generally, such interactions can affect other CD properties, such as geometry or fluorescence. Knowledge about CDs is growing, and their applicability as a biocompatible marker or carrier is being intensively studied. We believe that the atomic insight provided by our MD simulations may shed light on some of the fascinating phenomena in the novel field of CD studies.