Rotation by deformation and time-crystalline dynamics of cyclopropane molecule

A deformable body can rotate even with no angular momentum simply by changing its shape. Here the first all-atom level molecular dynamics example of this phenomenon is presented. For this the thermal vibrations of individual atoms in an isolated cyclopropane molecule are simulated in vacuum and at ultra-low internal temperature values. When the molecule is observed stroboscopically, at discrete equidistant time steps, the random thermal vibrations of the individual atoms become self-organized into a collective oscillatory motion of the entire molecule. The period of oscillation is emergent and intrinsic to the molecule so that this self-organization bears resemblance to a driven time crystal. The oscillation period increases in a self-similar manner when the length of the stroboscopic time step is increased. In the limit of very long stroboscopic time steps the entire molecule can then rotate in an apparent uniform fashion, but with no angular momentum. It is proposed that the observed behavior is universal in the case of triangular molecules. Moreover, it is shown that the emergent uniform rotation without any angular momentum, can be described in an effective theory approach as an autonomous Hamiltonian time crystal. The emergent oscillatory motion appears to be highly sensitive to temperature. This proposes that potential applications could be found from the development of molecular level detector to sensor and control technologies.

Geometric mechanics states that in the case of a deformable body, the vibrational and rotational motions are not separable [1][2][3][4][5][6][7][8].Even with no angular momentum, local vibrations in parts of the body can become self-organized into an emergent, global rotation of the entire body.A good example is a falling cat, how it maneuvers to rotate in air and lands on its feet.The effect is already in use for a variety of control purposes, for example the attitude of satellites can be controlled by periodic motions of parts of a satellite such as spinning rotors [6].For this kind of an effective rotational motion to emerge, the deformable body needs to have at least three movable components [1].
The original article [1] is entitled On rotation and vibration motions of molecules but it appears that until now, no first principles atomic level molecular example of the phenomenon has been described.In the present article the first atomic level example of the phenomenon is constructed, using the methods of all-atom molecular dynamics (MD) [9].All-atom MD is a methodology that is designed for a realistic and accurate description of a molecular system, at a semiclassical level.Here MD is employed as a tool, to scrutinize the collective response of an entire molecule to random thermal (or quantum) fluctuations of its individual atoms.The example that is considered is a single cyclopropane (C 3 H 6 ) molecule, shown in figure 1. Cyclopropane is a ring molecule that is made of three methylenes (CH 2 ) at the corners of an equilateral triangle.As such, it is the simplest organic ring molecule.In biochemistry cyclopropane appears as a building block of organic synthesis [10,11], and in medicine it is used as an anesthetics [12].
In order to observe the collective effects of the individual atom random thermal (quantum) fluctuations, the cyclopropane is monitored stroboscopically, at different equidistant time steps, and also by varying (successively increasing) the length of the time step.First, when the stroboscopic time scale becomes clearly longer than the characteristic period of a covalent bond oscillation, the random thermal (or quantum) fluctuations of the individual atoms cause the molecule to oscillate collectively.When the stroboscopic time scale increases further, these oscillations start self-organizing into apparent small amplitude back-and-forth rotational motions around the normal of the triangular plane.Finally, in the limit of very large stroboscopic time scales these rotational oscillations become a uniform rotation, veritably with no angular momentum: the apparent uniform rotation around the triangular symmetry axis that is observed in the limit of large stroboscopic time steps, is entirely due to rapid small amplitude shape deformations.This rotational motion can be modeled using an effective theory description that engages the three methyles as point-like interaction centers, at the corners of an equilateral triangle.Moreover, the effective theory description coincides with the Hamiltonian time crystal model that was previously introduced in [13]: a Hamiltonian time crystal describes autonomous energy conserving Hamiltonian dynamics, such that the minimum of the energy is in a state in periodic motion.Both autonomous and driven time crystals are presently under vigorous investigations, see for example [14][15][16][17][18][19][20][21][22][23][24].
Additional examples of simple triangular molecules that can be analyzed similarly and where effective theory timecrystalline dynamics can be expected to appear, include aziridine, oxirane and 1,2-dimethylcyclopropane [25].
The section 2 sets the stage by summarizing known results on the relation between vibrational and rotational motions [1][2][3][4][5][6][7][8].It is shown how, in the case of a deformable body with three identical point-like particles, the effective rotation angle relates to a Dirac monopole connection in the space of shapes; the Dirac monopole connection describes nontrivial holonomy which is the cause of rotation by shape deformation.The section then continues with a presentation of the simulation tools and methodology.The section 3 is the main body of the article.It presents all the simulations results, in the case of a cyclopropane molecule.The section 4 recalls the concept of a time crystal and in particular the effective theory description of Hamiltonian time crystals [8,13,24].The large stroboscopic time scale simulation results are interpreted in terms of this effective theory.The supplementary material (https://stacks.iop.org/NJP/23/073024/mmedia)includes three movies that show how the cyclopropane rotates, at different internal temperature values.

Theoretical background
The original articles on deformation induced rotation are [1][2][3][4] and reviews can be found in [5][6][7][8], and the following presentation follows closely [8]: consider the (time) t-evolution of three equal mass point-like objects, at the corners r i (i = 1, 2, 3) of a triangle.There are no external forces, thus the center of mass is stationary at all times r 1 (t) + r 2 (t) + r 3 (t) = 0. (1) The total angular momentum is assumed to vanish, Nevertheless, as first shown in [1,2] the triangle can rotate, by shape changes.To describe this rotation, the triangle is oriented to always lie on the z = 0 plane.Two triangles have the same shape when they only differ by a rigid rotation on the plane, around the z-axis.The shape is described by assigning internal shape coordinates s i (t) to the three vertices.These shape coordinates describe all possible triangular shapes, in an unambiguous fashion and in particular with no extrinsic rotational motion, when one demands that the vertex s 1 (t) always coincides with the positive x-axis so that s 1x (t) > 0 and s 1y (t) = 0.The vertex s 2 (t) is then restricted so that s 2y (t) > 0 and the third vertex s 3 (t) is fully determined as in (1), by demanding that the center of mass is at origin: Now let the shape of the triangle change arbitrarily but in a T-periodic fashion.Thus the s i (t) evolve also in a T-periodic fashion, The triangular shape then traces a closed loop Γ in the space of all possible triangular shapes.
At each time the shape coordinates s i (t) are related to the space coordinates r i (t), at most by a spatial rotation on the z = 0 plane, Initially θ(0) = 0.But if there is any net rotation due to shape changes, it is found that θ(T) = 0 so that the triangle has rotated during the period, by an angle θ(T), with θ(T)/T the stroboscopic angular velocity.The rotation angle θ(T) is obtained by substituting (3) into (2): Here the value of θ(T) is expressed as a line integral of a connection one-form [1][2][3][4][5][6][7][8].
in the space of all triangular shapes.Stokes theorem is used to convert the line integral into an integral over a surface Σ bounded by the loop Γ, Whenever (4) does not vanish the triangle returns to its initial shape, but with a net rotation by θ(T) = 0 around the axis through its center of mass.Since (2) vanishes there is no angular momentum and the rotation is emergent, it is entirely due to shape changes.
To demonstrate that the angle θ(T) is in general non-vanishing, the connection one-form A is evaluated as follows: start from standard Jacobi coordinates and change variables as follows: This is recognized as the connection one-form of a Dirac monopole at the origin r = 0 of R 3 with the Dirac string placed along the negative z-axis.In particular, since does not vanish, the rotation angle ( 4) is generically non-vanishing as it coincides with the flux of the Dirac monopole through the surface Σ with boundary Γ in the space of all triangular shapes.Notably, the Dirac string concurs with a shape where two corners overlap, and at the location of the monopole all three corners of the triangle coincide.

Simulation methods
The all-atom MD simulations are performed using the GROMACS package [26].The force field is CHARMM36m [27] and parameters are generated using CGenFF [28,29].A single cyclopropane molecule is simulated in vacuum under NVE conditions, with no ambient temperature or pressure coupling in the simulations.The boundary conditions are free and both Coulomb and van der Waals interactions extend over all atom pairs in the molecule with no cut-off approximation.The initial cyclopropane structure is taken from the PubChem web-site [30] where it is specified with 10 −14 m precision in both the C and H atomic coordinates.
The simulation process starts with a search of the minimum energy configuration of the potential energy contribution in CHARMM36m free energy.For this the initial PubChem configuration is first subjected to dissipation using 100 consecutive, repeated GROMACS simulations.Each of the simulations uses the final coordinates of the atoms obtained in the previous simulation as its initial configuration, but for a dissipation effect all the velocities of all atoms are set to zero; this exhausts the kinetic energy step-wise.Each trajectory has a length of 1.0 picoseconds (ps), and the time step is 0.01 femtoseconds (fs).Double precision floating point format is used in all these simulations.At the end of the 100th simulation the atoms are static, the structure is fully D 3h symmetric and in particular there is no kinetic energy left, with double precision accuracy.The configuration is then even further refined, using the L-BFGS optimization algorithm [31] until the maximal force becomes less than 10 −12 kJ mol −1 nm −1 .Finally, it is confirmed that the structure is indeed a locally stable potential energy minimum, by performing a 10 μs double precision simulation with 1.0 fs time step and vanishing initial velocities.During this entire simulation the total energy is conserved, there is no movement of the atoms observed.This confirms that the configuration is a static D 3h symmetric local minimum energy of the potential energy, with double precision accuracy.
The final configuration specifies all the atomic positions with better than Δx = 10 −24 m double precision accuracy.Similarly, all the atomic velocities vanish with better than Δv = 10 −11 m s −1 accuracy.Thus, both in the case of hydrogen and carbon atoms ΔxΔp and as such the structure is specified much more precisely than should be reasonable, by quantum mechanical considerations.To alleviate this, and in order to introduce minuscule thermal (or semiclassical quantum) fluctuations, the next step of preparatory simulations proceeds with single precision floating point format: even with single precision accuracy the atomic positions are determined with a precision that is better than e.g. the proton radius.
A single precision simulation is initiated from the double precision static D 3h symmetric local minimum energy of the potential energy, with all velocities set to zero.The single precision simulation uses a 1.0 fs time step.It is observed that during the first few time steps, the single precision potential energy slightly decreases from its initial double precision value.At the same time, all the covalent bond lengths start to oscillate, but with extremely tiny amplitudes.It is concluded that this single precision instability of the initial double precision minimum energy configuration is due to accumulation of round-up errors.The instability and the ensuing tiny atomic oscillations is caused by an incompatibility of the single precision and double precision floating point formats in combination with CHARMM36m parameter values.
After the first few initial steps during which the potential energy decreases, the total energy starts increasing but at a very slow rate: the single precision simulation is subject to a very slow energy drift, a phenomenon that is quite common in all-atom MD simulations.See GROMACS manual for detailed description [26].
When the covalent bonds of the atoms start oscillating, due to the energy drift, equipartition theorem can be used to assign the entire molecule an internal temperature value T. A GROMACS routine is used to estimate the temperature.The single precision simulation with energy drift is then continued until the cyclopropane molecule reaches a desired value of internal temperature T, at which point the configuration is saved.This saving is done at 50 different values of T, all with a T value below 1.0 K as determined by the equipartition theorem.These 50 configurations are the initial configurations for the 50 different final double precision production runs: each of the final double precision production run extends a full 10 μs of in silico time, with 1.0 fs time step.During these double precision production run simulations only insignificant energy drift is observed, the internal temperature values is very stable, and fluctuations around the average temperature value is small; despite the small number of atoms, it is concluded that the molecule is near a thermal equilibrium state during the entire production run.Tests at different (shorter) time steps have been performed, to conclude that results obtained using 1.0 fs time step remain numerically stable, over the entire simulation period; for single precision this is not the case, even at substantially shorter time step.

Results
In the case of a cyclopropane molecule, the time evolution θ(t) of the angle of rotation ( 4) is evaluated numerically (for t T, the characteristic time period for an individual atom oscillation) using all-atom MD.The instantaneous, stroboscopic values θ(t) are then reported, at increasing time scales.
The initial configurations for the production runs that evaluate this angle, are described in the methods section: they are a set of 50 'random' initial configurations each with a different internal temperature value T < 1.0 K with temperature determined by the equipartition theorem.As explained in the methods section, each initial configuration is designed to minimize the mechanical CHARMM36m free energy at the given internal temperature value T, with single precision floating point format.In particular, the angular momentum of each initial configuration vanishes with comparable numerical accuracy.
In the production runs double precision floating point format is used in combination with Δτ = 1.0 fs time steps, which is quite standard in high precision all-atom simulations.This is sufficient to stabilize the internal molecular temperature so that T remains essentially constant during each of the production runs, with vanishingly small energy drift.Each of the 50 production run trajectories that simulate 10 μs in silico, take about 90 h of wall clock time with the available processors.
To verify that the angular momentum indeed remains vanishingly small during all of the production runs, a methodology is introduced using a reference structure that builds on the concept of Eckart framing [32]; it should be noted that Eckart frames assume that the connection one-form (5), ( 7) is flat and removable which it is not the case [3].For the reference structure, a rotational motion is expressed as an accumulation of infinitesimal rotations, with each infinitesimal rotation describing the rotation of a 'virtual' rigid body that is defined by the instantaneous conformation of the molecule, at the given time.Thus, the instantaneous positions and velocities of all the individual C and H atoms are recorded at every femtosecond time step n of all production run trajectories.The instantaneous values L z (n) of angular momentum component that is normal to the plane of the three carbon atoms is then evaluated for the entire molecule, together with the corresponding instantaneous moment of inertia values I(n).This yields an estimate of the instantaneous apparent angular velocity values These instantaneous values ω(n) are then summed up, to compute the accumulated 'rigid body' rotation angle, that is due to the instantaneous angular momentum: In full compliance with the condition (2), in all the production run simulations the accumulated values (8) always remain less than ∼ 10 −6 radians, for all n and figure 2 shows a typical example: in the production simulations there is no observable rotation of the cyclopropane molecule, due to angular momentum.Accordingly any systematic rotational motion that exceeds ∼ 10 −6 radians during a production run must be emergent, and entirely due to shape deformations.
To observe such a systematic rotational motion the orientation of the cyclopropane structure is monitored stroboscopically, using time steps that are large in comparison to the frequency of a typical covalent bond oscillation.It is found that the cyclopropane molecule always rotates, and the rotation is always around the center of mass axis that is normal to the plane of the three carbon atoms.Moreover, the  8) in the production runs: there is no observable net rotation due to angular momentum.In this example the internal temperature is T = 0.065 K; the same trajectory is analyzed in figure 4.
angular velocity of the effective rotational motion is very sensitive to the value T of internal molecular temperature.
In the panels of figure 3 three different, generic examples of production run trajectories are presented; movies of these trajectories are available as supplementary material.Each trajectory describes how the molecule is observed stroboscopically, with 100 ps time steps.
In figure 3(a) the internal molecular temperature has the value T = 0.066 K.The molecule rotates in a uniform fashion clockwise i.e. with increasing θ(t) at a constant rate of ∼10 rad μs −1 around the axis that is normal to the plane of the carbon; see also the sketch in figure 5.In figure 3(b) T = 0.0087 K. Now the rotational motion is ratcheting, with a drift toward decreasing values of θ(t) at around ∼0.07 rad μs −1 .Such a combination of ratcheting and drifting motion is quite commonly observed, at low T-values and with those relatively short stroboscopic time steps that are attainable in the simulations.The drift implies that with larger stroboscopic time steps, in excess of 10 μs which is the length of a production trajectory, the rotation becomes increasingly uniform; this will be exemplified in figures 4(a)-(d).Finally, in figure 3(c) T = 0.94 K which is close to the upper bound of T-values accessible in the present simulations, while maintaining a good numerical stability.Now the molecule rotates in the counterclockwise direction at a rate around ∼15 rad μs −1 during the 10 μs simulation.The rotation remains slightly erratic at the stroboscopic time scale.Presumably, an increase in the stroboscopic time scale in combination with a longer MD trajectory produces a uniform rotation also in this case.
Unfortunately, with the available computers it is not feasible to try and produce trajectories that are much longer than 10 μs, with a high level of numerical stability.Such trajectories would be needed, to analyze in detail higher T-values.The present simulation results propose that there is a critical value T c probably around a few Kelvin, above which uniform stroboscopic rotation can no longer be observed.This could be a consequence of simulation artifacts that are due to issues such as energy drift [26].
In figures 4(a)-(d) the trajectory shown in figure 2 is analyzed in more detail; the internal temperature value T = 0.065 K in these panels is very close to the T-value of the trajectory in panel (a) of figure 3.In figure 4(a) θ(t) is sampled with very short stroboscopic time steps Δt = 200 fs.Oscillatory rotational motion is observed, with an amplitude that also oscillates but at much lower frequency.In figures 4(b) and (c) the time steps are increased to Δt = 2.0ps and Δt = 20ps, respectively.Now an increasingly regular ratcheting motion is observed, with apparent clockwise rotation, and with a decreasing relative amplitude; see also figure 3(b).Finally, in figure 4(d) where Δt = 100 ns the cyclopropane structure apparently rotates and in a uniform manner, akin the rotational motion shown in figure 3(a).See also figure 5. Despite a very small difference in T-values, the angular velocity is now clearly slower.Indeed, a peak of angular velocity that is centered at T ≈ 0.066 K is very narrow.Thus the rotational motion is apparently a very sensitive observable, in temperature.

Analysis
The individual atom covalent bond oscillations have a period that is near the femtosecond scale.The figures 4(a)-(d) reveal that there are also emergent (quasi) periodic oscillations that are detected at various very different stroboscopic time scales; those that are displayed in the panels extend from 10 −6 μs to 1.0 μs.This emergence of apparently periodic oscillations with an increase of the stroboscopic time scale is remarkably similar to the observations reported in [8].There, the rotation angle θ(t) in (4) was evaluated for three equal mass point particles that are subject to an external periodic driving force.Notably, there was no intrinsic dynamics and in particular no Hamiltonian was present in [8].Instead the shape coordinates (3) evolved as follows with variable amplitude |a| < 1 and frequencies ω 1,2 .The t-evolution of the angle θ(t) was then evaluated in [8] by a direct substitution of ( 9) in (4).
In figures 4(e)-(h) the t-evolution of the angle θ(t) as computed from ( 9) is shown for a = 0.1 and with ω 2 = 2ω 1 = 3; in [8] the values ω 2 = 2ω 1 = 2 were used.A comparison between the figures 4(a)-(d) and (e)-(h) shows remarkable qualitative similarity.This emergence of qualitatively similar periodic oscillations at different stroboscopic time scales, in two very different scenarios, is suggestive of universality.Moreover, in both cases the emergent large time scale periodicity is reminiscent of the phenomenon of discrete time crystals [20,21,23].The present results demonstrate that cyclopropane can realize emergent discrete timecrystalline dynamics, with different periodicity at different stroboscopic time scales.It is notable that this occurs for very specific parameter values: in the case of cyclopropane T ≈ 0.066 K, and in the case of driven triangle ω 2 = 2ω 1 .
In the final figures 4(d) and (h) the rotational motion is akin that of an equilateral triangle, that rotates uniformly at constant angular velocity around the symmetry axis that is normal to the plane of the triangle.Thus in the large time scale limit a cyclopropane molecule becomes a dynamical system in the shape of an equilateral triangle, with point-like interaction centers at the vertices, that rotates at constant angular velocity around the triangular symmetry axis.The limit is governed by an effective theory Hamiltonian model, in terms of a reduced set of variables; the effective theory coincides with the timecrystalline Hamiltonian model analyzed in [8,13,24].The model describes the time evolution of three point-like interaction centers x 1 , x 2 and x 3 in terms of the three link vectors The link vectors are subject to the Lie-Poisson brackets where abc is the standard permutation symbol and the two signs are related by parity.This bracket is designed to generate any kind of dynamics of the vertices x i except for stretching and shrinking of the links, {n k , |x i+1 − x i |} = 0 for all k, i.
This holds independently of the Hamiltonian.This approximation is appropriate in the case of a cyclopropane at very long time scales, where the oscillations of various covalent bond lengths can be replaced by their time averaged values.For clarity all |x i+1 − x i | are then chosen equal, so that the structure is an equilateral triangle.Hamilton's equation is With Hamiltonian the equation ( 10) does indeed describe uniform rotation of the triangle, with constant angular velocity around the symmetry axis that is normal to the triangle [13], as shown in figure 5. See supplementary material movie 1.
The two rotation directions around the equilateral symmetry axis of the cyclopropane are related by parity.In the present simulations the direction becomes selected randomly, essentially by floating point round-up errors, while in the effective Hamiltonian theory the direction changes with the sign in the Lie-Poisson bracket.In the case of actual cyclopropane, the molecule is quantum mechanical and it displays anomalous energetic and magnetic behaviors that are suggestive of a ring current [33][34][35].More generally, the ability to sustain a diamagnetic (diatropic) ring current in the presence of an external magnetic field is often considered as a defining characteristics of aromatic ring molecules [36,37].Ring current is a quantum manifestation of an effective electron current flow that appears in a semiclassical treatment of the molecule, and all-atom MD captures this current flow as an effective counter-rotational motion of the nuclear backbone.In the absence of a magnetic field, the two directions of electron current flow specify two degenerate quantum mechanical states.The ground state is a positive superposition of the two, and it is separated by a small energy gap from the negative superposition.An applied external magnetic field selects an energetically preferred direction for the current to flow.This is the essence of spontaneous symmetry breakdown.

Discussion and conclusions
In summary, the article describes results from state-of-art all-atom MD simulations that model the ground state properties of a single cyclopropane molecule in vacuum, at ultra-low internal molecular temperature values.Even when there is no angular momentum, an emergent oscillatory, ratcheting rotational motion is commonly observed in the simulations.This effective rotational motion is due to a collective self-organization of individual atom thermal vibrations, it takes place at characteristic time scales that are large in comparison to the frequency scale of covalent bond oscillations; alternatively, such individual atom vibrations could also have a quantum mechanical origin.When the molecule is followed with sufficiently long stroboscopic time steps the emergent rotation appears increasingly uniform, in a manner that can be reproduced by an effective theory timecrystalline Hamiltonian dynamics.The angular velocity of the rotation is found to be very sensitive to the internal molecular temperature.In the case of a single cyclopropane molecule studied here, the ability to produce uniform rotational motion without angular momentum appears to be limited to ultra-low internal molecular temperature values.But in the case of a larger molecular system the pertinent temperature scale for rotation should increase with increase in the molecular size.The effect could be possibly observed, e.g. using NMR spectroscopy.
The emergent timecrystalline rotation could be commonplace in the case of aromatic ring molecules.As a manifestation of aromatic ring current it could even serve as their defining characteristic.Its apparent high temperature sensitivity proposes that the phenomenon could be employed for a variety of detection, sensor and control purposes, in situations where extremely high precision is desirable.The phenomenon could even govern the in vivo working of biological macromolecules.Finally, for a proper analysis a full all-atom quantum MD should be performed.Unfortunately, such a simulation is out of reach for presently available computers.

Figure 1 .
Figure 1.A cyclopropane molecule C 3 H 6 has three covalently connected carbon atoms at the corners of an equilateral triangle, and each carbon atom is covalently bonded to two hydrogen atoms.The molecule has D 3h symmetry.

Figure 2 .
Figure 2. Generic time evolution of the angle ϑ(n) in equation (8) in the production runs: there is no observable net rotation due to angular momentum.In this example the internal temperature is T = 0.065 K; the same trajectory is analyzed in figure4.

Figure 3 .
Figure 3. Evolution of rotation angle (4) in cyclopropane, at different internal T-values and recorded at a stroboscopic time step Δt = 10 −10 s (with simulation time step 10 −15 s).In (a) T = 0.066 K, in (b) T = 0.0087 K, in (c) T = 0.94 K. See supplementary material for movies.

Figure 4 .
Figure 4. Panels (a)-(d) are for cyclopropane, and show the time evolution of rotation angle θ(t) along a T = 0.065 K trajectory.Sampling is done at different stroboscopic time steps Δt (with simulation time step 10 −15 s).In (a) Δt = 10 −13 s, in (b) Δt = 10 −12 s, in (c) Δt = 10 −11 s and in (d) Δt = 10 −7 s.Panels (e)-(h) are for a driven triangle of [8], as described by equation (9) with a = 0.1 and ω 2 = 2ω 1 = 3.The panels show the evolution of rotation angle θ(t) at increasing time scales; the time scale increases by a factor of ten between each panel.

Figure 5 .
Figure 5.Time evolution of the solution to (10) is akin the rotating cyclopropane shown in figures 3(a) and in 4(d); see supplementary material movie 1.