**Application of Molecular Dynamics Simulation to Small Systems**

Víctor M. Rosas-García and Isabel Sáenz-Tavera

*Facultad de Ciencias Químicas, Universidad Autónoma de Nuevo León, San Nicolás de los Garza, N. L. México*

#### **1. Introduction**

14 Will-be-set-by-IN-TECH

56 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

[30] G. R. W. Quispel & C. P. Dyt (1998). Volume-preserving integrators have linear error

[31] R. D. Ruth (1983). A canonical integration technique, *IEEE Trans. Nucl. Sci.* NS-30:

[32] J. M. Sanz-Serna (1992). Symplectic integrators for Hamiltonian problems: an overview,

[34] C. Tsallis (1988). Possible generalization of Boltzmann-Gibbs statistics, *J. Stat. Phys.* 52:

[35] M. Tuckerman, B. J. Berne, & G. J. Martyna (1992). Reversible multiple time scale

[36] H. Yoshida (1990). Construction of higher order symplectic integrators, *Phys. Lett.* A 150:

[37] S. Zai-jiu (1994). Construction of volume-preserving difference-schemes for source-free

[29] G. R. W. Quispel (1995). Volume-preserving integrators, *Phys. Lett. A* 206: 26-30.

[33] T. Schlick (2006). *Molecular Modeling and Simulation*, Springer-Verlag, Berlin.

molecular-dynamics, *J. Chem. Phys.* 97: 1990-2001.

systems via generating-functions, *J. Comp. Math.* 12: 265-272.

growth, *Phys. Lett.* A 242: 25-30.

*Acta Numerica* 1: 243-286.

2669-2671.

479-487.

262-268.

The study of chemical behavior includes answering questions as 'which isomer is the most stable?', 'which relative orientation is the most favorable for such-and-such interaction?', 'which conformer is the global minimum?', 'what are the lowest energy configurations and their relative energies?'. The answer to these questions–and many others–depends on the ability to find and study a variety of configurations of the system of interest. Recently, (Atilgan, 2007) briefly reviewed the use of molecular dynamics simulation for conformational search in the process of drug design, concluding that its use could reduce the errors in estimating binding affinities and finding more viable conformations. In addition, (Corbeil, 2009) considered the need to include ring flexibility in the conformational searches used in flexible docking. Most of the flexible docking algorithms skip searching for conformations in rings, even though a protein may stabilize a conformation other than the most stable one.

The need for a tool to examine the diverse configurations of the constituent particles of a system becomes obvious even as we consider relatively small systems (10-100 atoms). Finding by hand all the conformers of cyclohexane is feasible and maybe even instructive; this is somewhat more complex when doing morpholine and even something as small as an eight-membered heterocycle can be prohibitively complex to analyze by manipulating molecular models by hand. One of the methods available to the researcher to tackle this kind of problem is Molecular Dynamics (MD) simulation. MD allows an exploration of the configurational space of a system, respecting chemical constraints. Chemical constraints–such as atomic connectivities–are needed in cases such as conformations of molecular rings and configurations of molecular clusters, e.g., solvation shells. In both cases, atomic connectivities must be kept intact, otherwise we risk breaking up the ring or the molecular constituents of the cluster. In our research group we routinely employ MD simulation, sometimes in its semiempirical variety, to study small systems (systems with fewer than 100 atoms) whether they be solvation shells, inorganic clusters or heterocycles.

This review is narrowly focused on current software and methods appropriate for doing MD simulations of small heterocycles and clusters composed of 10-100 molecules. In particular, we try to systematize the tools available to tackle the problem of searching for minima in heterocycles and in molecular clusters. We want to use MD simulations as a tool to explore the energy landscape of a small system, so we can locate the global minimum. We do not include the vast literature simulating water solvation by MD. Even though the aggregates of water molecules around a solute can be considered clusters, we focus our attention on non-water clusters and we only borrow some tools for our purposes.

system of interest, such as boiling point or density, to name a few. The L-J potential illustrates one of the main stumbling blocks faced by the simulator, that is, the need to determine values

Application of Molecular Dynamics Simulation to Small Systems 59

The L-J potential as such is of little use in chemical systems because it does not consider chemically important bonds, such as covalent or hydrogen bonds. To take into account the *intra*molecular interactions, as well as the *inter*molecular interactions, the potential has to be defined in terms of the bond lengths, bond angles and torsions present in the molecule (its internal coordinates). When the potential energy functions are defined in terms of the molecular internal coordinates, the potential function is called a 'force-field'(Engler, 1973). There are many force-fields available for simulation of chemical systems: AMBER(Case, 2005; 2010; Pearlman, 1995), CHARMM(Brooks, 2009), MM3(Allinger, 1989), OPLS(Jorgensen, 1988) and TRIPOS(Clark, 1989), among others. In a MD simulation we use the force-field together with the equations of motion to obtain the dynamic behavior of the system (the trajectories of

The MD algorithm is a way to compute such trajectories and analyze them to extract physical information about a system. How do we describe the particle trajectories? At any given moment, each particle can be described in tridimensional space by 3 numbers, the x-coordinate, the y-coordinate, the z-coordinate. Since we are interested in the trajectories, we need to know the location of a particle and where it is headed at each point in time. So, we also need the momenta along each axis: *px*, *py* and *pz*. The first three numbers (XYZ coordinates) are known as *configuration space*; while the three momenta define the *momentum space*. Taken together they form the *phase space*. In the same way that a set of values for XYZ coordinates defines a point in configuration space, any set of values for these 6 coordinates (three XYZ and three momenta) defines a point in phase space. Inasmuch as both the positions and momenta depend on the time, as time changes the phase space coordinates also change, thus defining a trajectory in phase space. So a MD simulation is about computing the trajectories of particles

Fig. 1. Phase space coordinates of a particle with position in x,y,z and momenta *px*, *py*, *pz*

of empirical parameters, so the simulation has physical meaning.

the particles).

in phase space.

**2.3 Description of the trajectories**

When MD is based on molecular mechanics force-fields, it cannot model bond breaking-forming processes. For these cases, there are mixed methods such as Quantum Mechanics/Molecular Mechanics (QM/MM), where the bond breaking and formation is taken care of by the quantum mechanical part, while everything else is handled by the molecular mechanics/dynamics part. Another method that considers bond breaking-formation is Car-Parrinello Molecular Dynamics, where the electrons are also considered particles in the dynamics. Such methods are beyond the scope of this chapter.

#### **2. Basic concepts of molecular dynamics simulation**

Let us imagine a container (an imaginary box) with a finite number of particles. These particles move about the box at varying velocities (both speed and direction can vary), continuously colliding and bouncing off each other. The trajectories of the particles, taken as a set, contain valuable physical information. For example, if the particles interact with each other, making the particles move more slowly (equivalent to lowering the temperature) the interaction may allow the particles to associate and form a liquid. In a condensed phase, collisions happen more often, and the distance a particle can travel before colliding is much shorter than in the gas phase. It follows that, the movements of the particles correspond to the state of the system. It is by analyzing trajectories that we can compute properties. One can also hope that the results of such analysis will yield *physical insight* about the behavior of a system. For a full introduction to MD simulation see the book by Haile(Haile, 1992). It is of the foremost importance to have a correct description of the interactions between the particles.

#### **2.1 Is MD ready for general consumption?**

For a long time, MD simulation has been the province of specialists, and the literature is still packed with obscure vocabulary and long descriptions of complicated algorithms. All these are necessary when the purpose is the calculation of dynamical macroscopic properties of the system, such as viscosities, surface tensions or rheological properties. Unfortunately, finding such a scenario can put off the non-specialist that does not want to become an expert in the theory of molecular dynamics before attempting some configurational search.

When narrowly restricted to the task of configurational search, we can prescind with many of the details like periodic boundary conditions, equilibration, thermostats and the like. In addition, the field has produced software packages increasingly friendly to the user so, doing a MD simulation becomes just a little harder than using a molecular visualization program.

#### **2.2 Description of the interactions**

Every MD simulation depends on specifying an interaction model, that is, some analytical function that calculates the energy of interaction between the particles of a system (a potential). The Lennard-Jones (L-J) potential is probably the most popular potential in MD studies, (eq 1)

$$V\_{Lf}(r) = 4\varepsilon \left(\frac{\sigma}{r^{12}}\right) - \left(\frac{\sigma}{r^6}\right) \tag{1}$$

where *σ* is the particle diameter, *ε* is the depth of the energy well and *r* is the distance between the centers of the particles. *σ* and *ε* are determined empirically, to fit observed properties of the system of interest, such as boiling point or density, to name a few. The L-J potential illustrates one of the main stumbling blocks faced by the simulator, that is, the need to determine values of empirical parameters, so the simulation has physical meaning.

The L-J potential as such is of little use in chemical systems because it does not consider chemically important bonds, such as covalent or hydrogen bonds. To take into account the *intra*molecular interactions, as well as the *inter*molecular interactions, the potential has to be defined in terms of the bond lengths, bond angles and torsions present in the molecule (its internal coordinates). When the potential energy functions are defined in terms of the molecular internal coordinates, the potential function is called a 'force-field'(Engler, 1973). There are many force-fields available for simulation of chemical systems: AMBER(Case, 2005; 2010; Pearlman, 1995), CHARMM(Brooks, 2009), MM3(Allinger, 1989), OPLS(Jorgensen, 1988) and TRIPOS(Clark, 1989), among others. In a MD simulation we use the force-field together with the equations of motion to obtain the dynamic behavior of the system (the trajectories of the particles).

#### **2.3 Description of the trajectories**

2 Will-be-set-by-IN-TECH

molecules around a solute can be considered clusters, we focus our attention on non-water

When MD is based on molecular mechanics force-fields, it cannot model bond breaking-forming processes. For these cases, there are mixed methods such as Quantum Mechanics/Molecular Mechanics (QM/MM), where the bond breaking and formation is taken care of by the quantum mechanical part, while everything else is handled by the molecular mechanics/dynamics part. Another method that considers bond breaking-formation is Car-Parrinello Molecular Dynamics, where the electrons are also considered particles in the dynamics. Such methods are beyond the scope of this chapter.

Let us imagine a container (an imaginary box) with a finite number of particles. These particles move about the box at varying velocities (both speed and direction can vary), continuously colliding and bouncing off each other. The trajectories of the particles, taken as a set, contain valuable physical information. For example, if the particles interact with each other, making the particles move more slowly (equivalent to lowering the temperature) the interaction may allow the particles to associate and form a liquid. In a condensed phase, collisions happen more often, and the distance a particle can travel before colliding is much shorter than in the gas phase. It follows that, the movements of the particles correspond to the state of the system. It is by analyzing trajectories that we can compute properties. One can also hope that the results of such analysis will yield *physical insight* about the behavior of a system. For a full introduction to MD simulation see the book by Haile(Haile, 1992). It is of the foremost

importance to have a correct description of the interactions between the particles.

theory of molecular dynamics before attempting some configurational search.

*VL J*(*r*) = 4*ε*

For a long time, MD simulation has been the province of specialists, and the literature is still packed with obscure vocabulary and long descriptions of complicated algorithms. All these are necessary when the purpose is the calculation of dynamical macroscopic properties of the system, such as viscosities, surface tensions or rheological properties. Unfortunately, finding such a scenario can put off the non-specialist that does not want to become an expert in the

When narrowly restricted to the task of configurational search, we can prescind with many of the details like periodic boundary conditions, equilibration, thermostats and the like. In addition, the field has produced software packages increasingly friendly to the user so, doing a MD simulation becomes just a little harder than using a molecular visualization program.

Every MD simulation depends on specifying an interaction model, that is, some analytical function that calculates the energy of interaction between the particles of a system (a potential). The Lennard-Jones (L-J) potential is probably the most popular potential in MD

> *σ r*<sup>12</sup> − *σ r*6

where *σ* is the particle diameter, *ε* is the depth of the energy well and *r* is the distance between the centers of the particles. *σ* and *ε* are determined empirically, to fit observed properties of the

(1)

clusters and we only borrow some tools for our purposes.

**2. Basic concepts of molecular dynamics simulation**

**2.1 Is MD ready for general consumption?**

**2.2 Description of the interactions**

studies, (eq 1)

The MD algorithm is a way to compute such trajectories and analyze them to extract physical information about a system. How do we describe the particle trajectories? At any given moment, each particle can be described in tridimensional space by 3 numbers, the x-coordinate, the y-coordinate, the z-coordinate. Since we are interested in the trajectories, we need to know the location of a particle and where it is headed at each point in time. So, we also need the momenta along each axis: *px*, *py* and *pz*. The first three numbers (XYZ coordinates) are known as *configuration space*; while the three momenta define the *momentum space*. Taken together they form the *phase space*. In the same way that a set of values for XYZ coordinates defines a point in configuration space, any set of values for these 6 coordinates (three XYZ and three momenta) defines a point in phase space. Inasmuch as both the positions and momenta depend on the time, as time changes the phase space coordinates also change, thus defining a trajectory in phase space. So a MD simulation is about computing the trajectories of particles in phase space.

#### **2.4 Calculation of the trajectories**

The forces acting on that particle dictate where a particle is located in space and its direction of movement at each point in time. The classical way to deal with forces uses Newton's laws of motion. The first law states that "every object persists in its state of rest or uniform motion in a straight line unless it is compelled to change that state by forces impressed on it". So, if **r** is a vector that contains the particle coordinates at a given moment, and its first derivative with respect to time (its velocity) is symbolized by **˙r**, then the first law of motion can be mathematically expressed as

$$
\dot{\mathbf{r}} = \mathbf{c} constant \tag{2}
$$

where the new variable *rij* denotes the center-to-center distance between each pair of particles. This formidable-looking equation tells us that the total energy is calculated by adding the contribution of each particle and of each pair of particles to the kinetic and potential energy,

Application of Molecular Dynamics Simulation to Small Systems 61

Once an initial configuration of atoms is specified and values for location and momenta have been assigned to each atom in the system, the system is allowed to evolve as time progresses. This evolution causes a redistribution of energy, and allows the formation of an energy distribution characteristic of the temperature. This step is known as equilibration. The key step in the calculation of an equilibrated distribution is the determination of the time between collisions and the pairs of colliding particles, because the collisions are the ones responsible for the energy redistribution. After the system has achieved equilibration, we can register the trajectories of the particles. This is the simulation step, and it is the only stage when the trajectories have physical meaning. Once the trajectories have been calculated, properties can be estimated, as long as they can be formulated as averages over time. In dealing with small systems, let us say, macrocycles, we do not expect ever to achieve equilibration, because all the atoms that form the system have restrictions on them that preclude an accurate calculation

All the calculations are performed using finite-difference methods, of which Runge-Kutta is probably the best known, although the Runge-Kutta family of methods finds little use in MD simulations because of the large computational demands. One of the most widely used finite-difference methods is Verlet's algorithm, a third-order Störmer algorithm. It is not as

Each atom has its own velocity, which could take it in a direction very different from that of the other atoms so, what happens when a particle, atom or molecule, moves far away from the others? The usual way to deal with this problem is to employ *periodic boundary conditions* (PBC). In PBC we formally consider the system as made up by multiple copies of itself along all three X, Y, and Z axes. With this setup, if a particle wanders far enough from the others in one direction as to be located outside the box that contains the particles, another, identical, particle comes into the system from the opposite direction, bearing the same velocity. In general, when dealing with a single molecule–within the molecular mechanics formalism–we do not have to worry about losing atoms, because all the atoms are connected by chemical bonds, and molecular mechanics does not allow for bond breaking. In the case of clusters, it is conceivable that a single group (either a neutral molecule or an ion) might wander off the box limits, but that could give us information about the intensity of the interaction and about

Simulated annealing is a technique able to locate the global minimum of a system of particles. The concept is obtained by analogy with the process of annealing a metal, where the metal is heated to high temperature and then suddenly cooled down by submersion in water. By raising the temperature of the system, it leaves the local minimum where we happened to find it (or build it), and is able to sample the configuration space so it can find another energy minimum when lowering the temperature. Hopefully the new energetic minimum will be

stable as a Runge-Kutta, but its computational demands are much lower.

respectively.

of the parameters that indicate equilibrium.

**2.5 Keeping the system in one piece**

the optimum equilibrium geometry.

**2.6 Simulated annealing**

keeping in mind that quantities written in bold are vectors.

The second law of motion states that "force is equal to the change in momentum (mV) per change in time. For a constant mass, force equals mass times acceleration". If **F** is the force, *m* is the particle mass, and **¨r** indicates a second derivative of position with respect to time (or the first derivative of velocity with respect to time), then **¨r** corresponds to an acceleration. The second law can be mathematically expressed as follows,

$$\mathbf{F} = m\ddot{\mathbf{r}}\tag{3}$$

The third law states that "For every action, there is an equal and opposite re-action". Assuming an isolated system of identical particles, where the total net force, **F***total*, is zero, **F***total* = 0, then any force exerted by particle 1 on particle 2, **F**1, is compensated by an equal and opposite force, **F**2, exerted by particle 2 on particle 1. Using the same notation, this law can be expressed as,

$$\mathbf{F}\_1 = -\mathbf{F}\_2 \tag{4}$$

To calculate the trajectory in phase space, a MD simulation relies on solving Newton's equations of motion. We just need to use a slightly modified notation. Since we need to keep track of each particle, we use subindices, just as we did for the two particles used in the third law

$$\mathbf{F}\_{\mathbf{i}} = m\_{\mathbf{i}} \ddot{\mathbf{r}}\_{\mathbf{i}} \tag{5}$$

We now turn our attention to an interesting fact: the trajectories depend on time, but the mathematical form of the second law (see equation 3 is time-independent, that is, at any moment the relationship between forces, masses and accelerations is expressed by the same formula. So we expect to find a quantity that remains constant with time for the whole system of particles. In an isolated system, the total energy is constant with time, so this means that the sum of kinetic and potential energies for all the particles in the system is constant. This invariant quantity is also known as the Hamiltonian. The kinetic energy for each particle can be expressed by

$$E\_{ki} = \frac{1}{2} m\_i \dot{\mathbf{r}}\_i \tag{6}$$

while the potential energy is calculated according the the model for description of the interactions, such as the Lennard-Jones potential (see equation 1), although it could be any of the force-fields available in the literature. So the form of the equation for the total energy might be

$$E\_{\text{total}} = E\_{\text{ki}} + E\_{\text{pot}} = \Sigma\_i \frac{1}{2} m\_i \dot{\mathbf{r}}\_i + \Sigma\_i \Sigma\_j 4\varepsilon \left(\frac{\sigma}{r\_{ij}^{12}}\right) - \left(\frac{\sigma}{r\_{ij}^6}\right) \tag{7}$$

4 Will-be-set-by-IN-TECH

The forces acting on that particle dictate where a particle is located in space and its direction of movement at each point in time. The classical way to deal with forces uses Newton's laws of motion. The first law states that "every object persists in its state of rest or uniform motion in a straight line unless it is compelled to change that state by forces impressed on it". So, if **r** is a vector that contains the particle coordinates at a given moment, and its first derivative with respect to time (its velocity) is symbolized by **˙r**, then the first law of motion can be

The second law of motion states that "force is equal to the change in momentum (mV) per change in time. For a constant mass, force equals mass times acceleration". If **F** is the force, *m* is the particle mass, and **¨r** indicates a second derivative of position with respect to time (or the first derivative of velocity with respect to time), then **¨r** corresponds to an acceleration. The

The third law states that "For every action, there is an equal and opposite re-action". Assuming an isolated system of identical particles, where the total net force, **F***total*, is zero, **F***total* = 0, then any force exerted by particle 1 on particle 2, **F**1, is compensated by an equal and opposite force, **F**2, exerted by particle 2 on particle 1. Using the same notation, this law can be expressed as,

To calculate the trajectory in phase space, a MD simulation relies on solving Newton's equations of motion. We just need to use a slightly modified notation. Since we need to keep track of each particle, we use subindices, just as we did for the two particles used in the

We now turn our attention to an interesting fact: the trajectories depend on time, but the mathematical form of the second law (see equation 3 is time-independent, that is, at any moment the relationship between forces, masses and accelerations is expressed by the same formula. So we expect to find a quantity that remains constant with time for the whole system of particles. In an isolated system, the total energy is constant with time, so this means that the sum of kinetic and potential energies for all the particles in the system is constant. This invariant quantity is also known as the Hamiltonian. The kinetic energy for each particle can

> *Eki* <sup>=</sup> <sup>1</sup> 2

> > 1 2

*Etotal* = *Eki* + *Epoti* = Σ*<sup>i</sup>*

while the potential energy is calculated according the the model for description of the interactions, such as the Lennard-Jones potential (see equation 1), although it could be any of the force-fields available in the literature. So the form of the equation for the total energy

*mi***˙r***<sup>i</sup>* + Σ*i*Σ*j*4*ε*

 *σ r*12 *ij*

 − *σ r*6 *ij*

**˙r** = *constant* (2)

**F** = *m***¨r** (3)

**F**<sup>1</sup> = −**F**<sup>2</sup> (4)

**F***<sup>i</sup>* = *mi***¨r***<sup>i</sup>* (5)

*mi***˙r***<sup>i</sup>* (6)

(7)

**2.4 Calculation of the trajectories**

mathematically expressed as

third law

be expressed by

might be

keeping in mind that quantities written in bold are vectors.

second law can be mathematically expressed as follows,

where the new variable *rij* denotes the center-to-center distance between each pair of particles. This formidable-looking equation tells us that the total energy is calculated by adding the contribution of each particle and of each pair of particles to the kinetic and potential energy, respectively.

Once an initial configuration of atoms is specified and values for location and momenta have been assigned to each atom in the system, the system is allowed to evolve as time progresses. This evolution causes a redistribution of energy, and allows the formation of an energy distribution characteristic of the temperature. This step is known as equilibration. The key step in the calculation of an equilibrated distribution is the determination of the time between collisions and the pairs of colliding particles, because the collisions are the ones responsible for the energy redistribution. After the system has achieved equilibration, we can register the trajectories of the particles. This is the simulation step, and it is the only stage when the trajectories have physical meaning. Once the trajectories have been calculated, properties can be estimated, as long as they can be formulated as averages over time. In dealing with small systems, let us say, macrocycles, we do not expect ever to achieve equilibration, because all the atoms that form the system have restrictions on them that preclude an accurate calculation of the parameters that indicate equilibrium.

All the calculations are performed using finite-difference methods, of which Runge-Kutta is probably the best known, although the Runge-Kutta family of methods finds little use in MD simulations because of the large computational demands. One of the most widely used finite-difference methods is Verlet's algorithm, a third-order Störmer algorithm. It is not as stable as a Runge-Kutta, but its computational demands are much lower.

#### **2.5 Keeping the system in one piece**

Each atom has its own velocity, which could take it in a direction very different from that of the other atoms so, what happens when a particle, atom or molecule, moves far away from the others? The usual way to deal with this problem is to employ *periodic boundary conditions* (PBC). In PBC we formally consider the system as made up by multiple copies of itself along all three X, Y, and Z axes. With this setup, if a particle wanders far enough from the others in one direction as to be located outside the box that contains the particles, another, identical, particle comes into the system from the opposite direction, bearing the same velocity. In general, when dealing with a single molecule–within the molecular mechanics formalism–we do not have to worry about losing atoms, because all the atoms are connected by chemical bonds, and molecular mechanics does not allow for bond breaking. In the case of clusters, it is conceivable that a single group (either a neutral molecule or an ion) might wander off the box limits, but that could give us information about the intensity of the interaction and about the optimum equilibrium geometry.

#### **2.6 Simulated annealing**

Simulated annealing is a technique able to locate the global minimum of a system of particles. The concept is obtained by analogy with the process of annealing a metal, where the metal is heated to high temperature and then suddenly cooled down by submersion in water. By raising the temperature of the system, it leaves the local minimum where we happened to find it (or build it), and is able to sample the configuration space so it can find another energy minimum when lowering the temperature. Hopefully the new energetic minimum will be lower in energy than the previous one. In a MD simulation, we can maintain the system at high temperature (even unrealistically high temperatures, like 2000 K) and then the system temperature is reduced. This adds a cooling step to the simulation.

simulate chemical reactions by molecular dynamics, because it is not limited to molecular

Application of Molecular Dynamics Simulation to Small Systems 63

MacroModel(Mohamadi, 1990) is a well-developed molecular dynamics package for

TINKER(Ponder, 2011) is a molecular dynamics package created and maintained by the group of Jay W. Ponder. It employs molecular mechanics and currently lacks a graphical user

VASP(Kresse, 1993; 1994; 1996a;b) seems to be very popular in the metal clusters community. It can do MD simulations using density functional theory (considered as part of Ab Initio Molecular Dynamics) and its use does require the skills of an expert computational chemist.

It is well known that MD is a very inefficient way to search for minima in large rings(Saunders, 1990). However, it should not be underestimated for searches in medium and small rings (rings < 10 atoms). The use of MD simulations for conformational search is in large part

Isayev et al.(Isayev, 2007) claim that the pyrimidine ring in nucleic acid bases has a range of effective non-planar conformations under ambient conditions. Saiz et al.(Saiz, 1996) demonstrated that MD simulation with the Tripos force-field was good at reproducing the conformational behavior of a dioxo ring in aqueous solution. Rosas-García et al.(Rosas-Garcia, 2010) studied conformations of fosforinanes using MD simulations in Ghemical, although some parameters had to be determined at the Hartree-Fock level (basis set 6-31G\*). An extensive search by MD found all the conformers for two diastereomers, and a comparison of the global minima allowed to explain why the axial preference of the phenyl ring was linked

Sometimes a side chain can modify the conformational behavior of small rings, like the case of Tosco et al. (Tosco, 2008) who used MD simulations in CHARMM to do conformational search on a series of cyclic oxadiazolol, and thiadiazolol isosteres of carboxylate, where most of the conformational freedom came from the side chain attached to the ring. A similar case is that of Bombasaro et al.(Bombasaro, 2008) who used GROMACS in combination with a systematic grid conformational search to study bullacin. Bullacin contains three five-membered rings, with one 11-carbon side chain, and two rings joined by a 12-carbon chain (see Figure 2).

Regarding structures with fused rings, Aleksandrov and Simonson(Aleksandrov, 2006; 2009) reported development of CHARMM22 parameters (employing the TIP3P(Jorgensen, 1983) model for water) for several tetracycline derivatives and a tetracycline/Mg2<sup>+</sup> complex (see the structure in Figure 2). Tetracyclines exist in two conformations, twisted and extended. In this case the interest was not so much the conformational variety but the possibility of several protonation sites, and the uncertainty of the binding sites for metals such as Mg2<sup>+</sup>. They employed TINKER and the MM3 force-field for MD/simulated annealing. On the other hand, Kiliç et al.(Kiliç, 2000) found out by simulations using quantum MD that cubane and

their group 14 analogs can convert to eight-membered rings at high temperature.

favored by convenience, given that many different software programs include it.

to the relative configuration of the stereocenters in the molecule.

biomolecules, and it includes a polished interface known as Maestro.

interface. Its learning curve is somewhat steep.

**4. Conformational analysis of heterocycles**

mechanics parameters.
