**3. Software that implements molecular dynamics simulation**

AMBER(Case, 2005; 2010; Pearlman, 1995) designates two different things: a force-field and a package for MD simulation. AMBER the package uses AMBER the force-field for its calculations, but it is entirely possible to use AMBER the force-field in a non-AMBER package, such as GROMACS or CHARMM. Beware, though, that using the same force-field in two different packages will not necessarily get identical results. AMBER the package is currently in version 11. Its learning curve is steep. It is possible to simulate small species, with at least one tutorial showing how to do it.

CHARMMBrooks (2009) also shares the situation of AMBER, in that the name designates both a force-field and a computational engine for MD simulation. It is also possible to use CHARMM the force-field in a non-CHARMM engine. Similar caveats apply, although these authors could not find any information on using CHARMM for small molecules.

Gabedit(Allouche, 2011) makes quantum chemistry software accessible to the novice modeller. It presents to the user a rather limited array of options, making for a less confusing experience. On the other hand, this means that to take full advantage of the capabilities of the quantum chemistry software, the user needs to be well versed in the respective manuals. Gabedit can perform MD simulations by itself, using the AMBER99 force-field. It can also setup MD simulations using semiempirical quantum mechanical energy evaluations, and submit them to a variety of computational engies, such as MOPAC2009(Stewart, 2009), ORCA(Radoul, 2010) or FireFly(Granovsky, 2011). When using quantum-mechanical energy evaluations, the user should keep in mind that these methods allow for bond breaking and formation, so it is entirely possible to end up with an isomerized structure.

Ghemical(Hassinen, 2001) can perform MD simulations, and defines a graphical user interface for this. It has a convenient facility to generate water solvation boxes. The user has to make sure that the force-field contains appropriate parameters. New parameters can be added by editing some configuration files. Its graphical user interface makes Ghemical very accessible to the beginning modeller.

GROMACS(Berendsen, 1995; Hess, 2008; Lindahl, 2001; van der Spoel, 2005) is a molecular dynamics software tailored to simulations with hundreds of millions of particles. It is certainly not made with the novice user in mind, and its learning curve is steep. However, it is an extremely fast computational engine. It is not recommended for dynamics of small species because then the advantages of the fast computation are lost, and because it is command-line based. Building the required topology before running the simulation can be daunting to a novice.

AMBER, CHARMM and GROMACS are tailored towards simulation of bio-macromolecules, and the force-fields included reflect this. They contain force-fields highly optimized for aminoacids or nucleotides or carbohydrates.

HyperChem(Hyperchem, 2011) is a commercial software product has a module for MD simulation, in addition to ab initio, density functional and semiempirical capabilities. It can 6 Will-be-set-by-IN-TECH

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

AMBER(Case, 2005; 2010; Pearlman, 1995) designates two different things: a force-field and a package for MD simulation. AMBER the package uses AMBER the force-field for its calculations, but it is entirely possible to use AMBER the force-field in a non-AMBER package, such as GROMACS or CHARMM. Beware, though, that using the same force-field in two different packages will not necessarily get identical results. AMBER the package is currently in version 11. Its learning curve is steep. It is possible to simulate small species, with at least

CHARMMBrooks (2009) also shares the situation of AMBER, in that the name designates both a force-field and a computational engine for MD simulation. It is also possible to use CHARMM the force-field in a non-CHARMM engine. Similar caveats apply, although these

Gabedit(Allouche, 2011) makes quantum chemistry software accessible to the novice modeller. It presents to the user a rather limited array of options, making for a less confusing experience. On the other hand, this means that to take full advantage of the capabilities of the quantum chemistry software, the user needs to be well versed in the respective manuals. Gabedit can perform MD simulations by itself, using the AMBER99 force-field. It can also setup MD simulations using semiempirical quantum mechanical energy evaluations, and submit them to a variety of computational engies, such as MOPAC2009(Stewart, 2009), ORCA(Radoul, 2010) or FireFly(Granovsky, 2011). When using quantum-mechanical energy evaluations, the user should keep in mind that these methods allow for bond breaking and formation, so it is

Ghemical(Hassinen, 2001) can perform MD simulations, and defines a graphical user interface for this. It has a convenient facility to generate water solvation boxes. The user has to make sure that the force-field contains appropriate parameters. New parameters can be added by editing some configuration files. Its graphical user interface makes Ghemical very accessible

GROMACS(Berendsen, 1995; Hess, 2008; Lindahl, 2001; van der Spoel, 2005) is a molecular dynamics software tailored to simulations with hundreds of millions of particles. It is certainly not made with the novice user in mind, and its learning curve is steep. However, it is an extremely fast computational engine. It is not recommended for dynamics of small species because then the advantages of the fast computation are lost, and because it is command-line based. Building the required topology before running the simulation can be daunting to a

AMBER, CHARMM and GROMACS are tailored towards simulation of bio-macromolecules, and the force-fields included reflect this. They contain force-fields highly optimized for

HyperChem(Hyperchem, 2011) is a commercial software product has a module for MD simulation, in addition to ab initio, density functional and semiempirical capabilities. It can

authors could not find any information on using CHARMM for small molecules.

entirely possible to end up with an isomerized structure.

temperature is reduced. This adds a cooling step to the simulation.

**3. Software that implements molecular dynamics simulation**

one tutorial showing how to do it.

to the beginning modeller.

aminoacids or nucleotides or carbohydrates.

novice.

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

MacroModel(Mohamadi, 1990) is a well-developed molecular dynamics package for biomolecules, and it includes a polished interface known as Maestro.

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 interface. Its learning curve is somewhat steep.

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.

#### **4. Conformational analysis of heterocycles**

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 favored by convenience, given that many different software programs include it.

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 to the relative configuration of the stereocenters in the molecule.

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.

**5. Configurational search of clusters by molecular dynamics simulation**

the configurational space of clusters.

A lot of work on clusters is based on atomic clusters. In this case it may suffice to develop a random-placement algorithm to generate structural variety. Many researchers have taken this route. Chen et al.(Chen, 2011) employed VASP to study clusters of metal carbides. They first generated a variety of Ca8 clusters and used step-wise addition of carbon atoms and geometry optimization after each step. MD simulations were used to evaluate the thermodynamic stability of several cage structures. The presence of only small distortions in the cage structures at 400 K was taken as evidence of their thermodynamic stability. Fujima and Oda employed VASP to study titanium clusters adsorbed on a single wall nano-capsule. There is no description of the parameters used in the MD calculations (temperature, time step, total simulation time or any others). The only configurational searching took place by putting Ti atoms on different adsorption sites on the carbon wall and doing geometry optimization. Given that they attempted to maximize the contact area between the cluster and the nano-wall, such approach seems justified(Fujima, 2009). Jian-Song and Li(Jiang-Song, 2010) did a configurational search for Ga7As7 clusters by randomly choosing points in space from a tridimensional box, cage or sphere, applying distance constraints to keep the atoms at chemically reasonable distances. This is not MD, and they had to generate thousands of structures, although the original paper is sketchy on the details of how many structures were generated. It could well be that the imposed distance constraints biased the resulting structures. This method does have the advantage of not requiring force-field parameters for generating the structures. Jiménez-Sáez(Jiménez-Sáez, 2006) studied the equilibrium structures of copper clusters as a function of the kinetic energy of deposition on a gold surface. For the starting geometries, no configurational search was done. The configurations of the deposited clusters were analyzed in terms of the deformation produced as the kinetic energy of deposition varied. Kuzmin et al.(Kuzmin, 2008) used software developed in-house to study configurations of silver clusters using MD simulations with the embedded atom model on clusters between 13 and 2057 atoms. They used temperatures between 0 and 1300 K. Li et al.(Li, 2007) used full-potential linear-muffin-tin-orbital MD (FP-LMTO-MD) calculations to study the effect of Al impurities on Si clusters. The method is suitable for semiconductor and metal clusters. As for the generation of the initial structures, the authors took the reported ground states of the silicon clusters and added or substituted Al atoms in all possible positions of each cluster. Yang and Xiong(Yang, 2008) used a similar method to generate the initial geometries for FeB*<sup>n</sup>* clusters. So there appears to be a need for more systematically searching

Application of Molecular Dynamics Simulation to Small Systems 65

In our literature review, we found only one recent example of a configurational search on a cluster using MD simulations, that of Chandrachud et al.(Chandrachud, 2009). These researchers employed VASP to study gold cages using Born-Oppenheimer MD to generate initial configurations. They did MD runs at four different constant temperatures, for 60 ps each, obtaining 600 structures (150 for each constant temperature run). Geometry

Clusters are not limited to groups of separate individual particles like in the case of metal clusters. When we have polyatomic species involved in cluster formation, it becomes important to maintain chemical bonds intact. Shiroishi et al.(Shiroishi, 2005) used Car-Parrinello Molecular Dynamics to study iron oxide clusters. Doll et al.(Doll, 2010) used Ab Initio Molecular Dynamics calculations to study clusters of lithium fluoride. They used a two-part protocol: first generation of candidate structures by means of simulated annealing at

optimization of these structures yielded 50 distinct isomeric clusters.

Fig. 2. Some heterocycles and fused cycles studied by molecular dynamics simulations.

8 Will-be-set-by-IN-TECH

Fig. 2. Some heterocycles and fused cycles studied by molecular dynamics simulations.

#### **5. Configurational search of clusters by molecular dynamics simulation**

A lot of work on clusters is based on atomic clusters. In this case it may suffice to develop a random-placement algorithm to generate structural variety. Many researchers have taken this route. Chen et al.(Chen, 2011) employed VASP to study clusters of metal carbides. They first generated a variety of Ca8 clusters and used step-wise addition of carbon atoms and geometry optimization after each step. MD simulations were used to evaluate the thermodynamic stability of several cage structures. The presence of only small distortions in the cage structures at 400 K was taken as evidence of their thermodynamic stability. Fujima and Oda employed VASP to study titanium clusters adsorbed on a single wall nano-capsule. There is no description of the parameters used in the MD calculations (temperature, time step, total simulation time or any others). The only configurational searching took place by putting Ti atoms on different adsorption sites on the carbon wall and doing geometry optimization. Given that they attempted to maximize the contact area between the cluster and the nano-wall, such approach seems justified(Fujima, 2009). Jian-Song and Li(Jiang-Song, 2010) did a configurational search for Ga7As7 clusters by randomly choosing points in space from a tridimensional box, cage or sphere, applying distance constraints to keep the atoms at chemically reasonable distances. This is not MD, and they had to generate thousands of structures, although the original paper is sketchy on the details of how many structures were generated. It could well be that the imposed distance constraints biased the resulting structures. This method does have the advantage of not requiring force-field parameters for generating the structures. Jiménez-Sáez(Jiménez-Sáez, 2006) studied the equilibrium structures of copper clusters as a function of the kinetic energy of deposition on a gold surface. For the starting geometries, no configurational search was done. The configurations of the deposited clusters were analyzed in terms of the deformation produced as the kinetic energy of deposition varied. Kuzmin et al.(Kuzmin, 2008) used software developed in-house to study configurations of silver clusters using MD simulations with the embedded atom model on clusters between 13 and 2057 atoms. They used temperatures between 0 and 1300 K. Li et al.(Li, 2007) used full-potential linear-muffin-tin-orbital MD (FP-LMTO-MD) calculations to study the effect of Al impurities on Si clusters. The method is suitable for semiconductor and metal clusters. As for the generation of the initial structures, the authors took the reported ground states of the silicon clusters and added or substituted Al atoms in all possible positions of each cluster. Yang and Xiong(Yang, 2008) used a similar method to generate the initial geometries for FeB*<sup>n</sup>* clusters. So there appears to be a need for more systematically searching the configurational space of clusters.

In our literature review, we found only one recent example of a configurational search on a cluster using MD simulations, that of Chandrachud et al.(Chandrachud, 2009). These researchers employed VASP to study gold cages using Born-Oppenheimer MD to generate initial configurations. They did MD runs at four different constant temperatures, for 60 ps each, obtaining 600 structures (150 for each constant temperature run). Geometry optimization of these structures yielded 50 distinct isomeric clusters.

Clusters are not limited to groups of separate individual particles like in the case of metal clusters. When we have polyatomic species involved in cluster formation, it becomes important to maintain chemical bonds intact. Shiroishi et al.(Shiroishi, 2005) used Car-Parrinello Molecular Dynamics to study iron oxide clusters. Doll et al.(Doll, 2010) used Ab Initio Molecular Dynamics calculations to study clusters of lithium fluoride. They used a two-part protocol: first generation of candidate structures by means of simulated annealing at

of interest and how many of them are to be added. Ghemical is less flexible in this regard, because it has tools only for building solvation shells, and the number of water molecules added is not specified directly by the user, but by the volume specified for the water box. Adding a precise number of water molecules in Ghemical can be a hit-and-miss experience.

Application of Molecular Dynamics Simulation to Small Systems 67

The set of 'physical' information contains: temperature, pressure, relaxation times, compresibilities, whether the simulation will be at constant temperature or constant pressure and–probably most important of all–the force-field used to evaluate the interactions within a molecule and between molecules. Given that the trajectories should have physical meaning, how do we know what kind of experimental conditions are we simulating? This corresponds to the choice of the ensemble. We should be familiar with three ensembles: the microcanonical

The microcanonical ensemble maintains constant number of particles (N), constant volume (V) and constant total energy (E), so it is also known as the NVE ensemble. The canonical ensemble keeps a constant number of particles, constant volume and constant temperature

The isothermal-isobaric ensemble keeps constant number of particles, constant pressure and constant temperature, so it is also known as the NPT ensemble. Physically, the NPT ensemble is the most important in chemistry, because many chemical processes are performed under

For our particular situation, when we deal with so few molecules that even the concepts of pressure and temperature are not well defined, it suffices to say that these ensembles are different ways to give energy to the system and any one of them can accomplish the task of

In our group we typically choose the program defaults, as we are not interested in the physical meaning of the trajectories, but only in the energetic minima resulting from the dynamic

Choosing a force-field can be daunting to a novice, because of all the options available. In terms of the specific strengths and weaknesses of each force-field, the reader is referred to the literature. However, the main roadblock in using molecular mechanics force-fields is that, sooner or later, one wants to study a molecule lacking adequate parameters in any force-field. Here, the user of MD software needs to know that some software packages, particularly the most friendly to the user, sometimes allow a dynamics calculation to run substituting default values for the missing parameters. Such calculations have practically no value at all. Ghemical uses the TRIPOS force-field, but the user should be aware of the error messages because usually many parameters are missing, and Ghemical will substitute default values. Gabedit will not run a dynamics unless all the parameters are defined or one decides to use semiempirical methods. The lack of adequate parameters usually requires doing ab initio calculations on a model compound, so the parameters can be generated. This route is reasonable when the molecules of interest are large compared to the model molecule, but what

are we supposed to do if the molecule and the model compound are the same?

ensemble, the canonical ensemble and the isothermal-isobaric ensemble.

taking the system out of an energy well and into another one.

**6.2 Choice of physical conditions**

(T), so it is also called the NVT ensemble.

constant pressure and temperature.

search.

**6.3 Choice of force-field**

a low level of theory (Hartree-Fock) and, second, optimization of the obtained structures using the Local Density Approximation. Takayanagi(Takayanagi, 2008) studied clusters of solvated glycine using the PM6 Hamiltonian. Their semiempirical MD simulations were performed at 300 K, and the initial geometries were taken from previously reported higher-level results and reoptimized using PM6. They observed dissociation of the proton from a carboxylate group, although could not observe formation of the zwitterion. In our group we have studied calcium carbonate clusters(Rosas-Garcia, 2011), and we have explored the configurational space by means of semiempirical MD simulations, using the PM6 Hamiltonian in MOPAC2009. Pang et al.(Pang, 1994) studied inclusion compounds in cycloalkanes by simulated annealing using HyperChem using the MM+ force-field (a variant of MM2. We cannot recommend the use of MM+ due to the lack of a published description of the modifications) and varying the temperature from 300 to 1000 K in 100 K increments. The dynamics revealed that there was orientational flexibility within the cycle and that the interconversion barriers were as low as 1 kcal/mole.
