Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS

CHAPTER 2

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS This chapter has been divided into three parts as follows: Chapter 2.1: Overview of BIOVIA Materials Studio Chapter 2.2: Overview of LAMMPS Chapter 2.3: Overview of GROMACS

Chapter 2.1

Overview of BIOVIA Materials Studio Sumit Sharma, Pramod Kumar, Rakesh Chandra Department of Mechanical Engineering, Dr. B. R. Ambedkar National Institute of Technology, Jalandhar, India

BIOVIA Materials Studio [1] is software for simulating and modeling materials. It is developed and distributed by BIOVIA (formerly Accelrys), a firm specializing in research software for computational chemistry, bioinformatics, cheminformatics, molecular dynamics simulation, and quantum mechanics. This software is used in advanced research of various materials, such as polymers, carbon nanotubes, catalysts, metals, and ceramics, by top universities, research centers, and high-tech companies. BIOVIA Materials Studio is a client-server model software package with Microsoft Windows-based PC clients and Windows- and Linux-based servers running on PCs, Linux IA-64 workstations, and HP XC clusters. BIOVIA Materials Studio [1] offers a rich set of features for material modeling at the atomic level. There are various modules available in BIOVIA Materials Studio for modeling, visualization, and analysis of material systems. Its graphic user interface is very clear and intuitive, which offers a high-quality, Windows-standard environment into which we can plug any BIOVIA Materials Studio product. In the following sections, various modules of BIOVIA Materials Studio have been discussed along with suitable examples wherever required. Molecular Dynamics Simulation of Nanocomposites using BIOVIA Materials Studio, Lammps and Gromacs https://doi.org/10.1016/B978-0-12-816954-4.00002-4 # 2019 Elsevier Inc. All rights reserved.

39

40

Chapter 2

2.1.1 Modules The various modules of BIOVIA Materials Studio have been discussed here. The readers can easily follow these steps for modeling in BIOVIA Materials Studio: (a) Materials Visualizer Materials Visualizer provides fast, interactive tools that enable you to construct graphic models of molecules, crystalline materials, surfaces, interfaces, layers, and polymers. You can manipulate, view, and analyze these models. Materials Visualizer also handles graph, tabular, and textual data and provides the software infrastructure and analysis tools to support the full range of BIOVIA Materials Studio products. Materials Visualizer can be run as a standalone tool for building, visualizing, and editing structures. Interaction with Windows productivity tools allows easy sharing and reporting of results and data. (b) Adsorption Locator Adsorption Locator enables us to simulate a substrate loaded with an adsorbate or an adsorbate mixture of a fixed composition. Adsorption Locator is designed for the study of individual systems, allowing you to find low-energy adsorption sites on both periodic and nonperiodic substrates or to investigate the preferential adsorption of mixtures of adsorbate components, for example. Adsorbates are typically molecular gases or liquids, and substrates are usually porous crystals or surfaces, such as zeolites or CNTs, or amorphous structures, such as silica gel or activated carbon. Adsorption Locator identifies possible adsorption configurations by carrying out Monte Carlo searches of the configurational space of the substrate-adsorbate system as the temperature is slowly decreased. (c) Amorphous Cell Amorphous Cell is a suite of computational tools that allow us to construct representative models of complex amorphous systems and to predict key properties. Among the properties that we can predict and investigate are cohesive energy density, equation-of-state behavior, chain packing, and localized chain motions. The methodology of Amorphous Cell construction is based on an extension of well-established methods for generating bulk-disordered systems containing chain molecules in realistic equilibrium conformations. Other features include provision for construction of arbitrary mixture systems containing any combination of small molecules and polymers, in addition to special capabilities for producing ordered nematic mesophases and slabs of amorphous material suitable for use in creating models of interphases, as would be required to study adhesion or lubrication. (d) COMPASS COMPASS is a powerful force field supporting atomistic simulations of condensed-phase materials. COMPASS stands for condensed-phase optimized molecular potentials for atomistic

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 41 simulation studies. COMPASS is the first ab initio force field that has been parameterized and validated using condensed-phase properties, in addition to various ab initio and empirical data for molecules in isolation. Consequently, this force field enables accurate and simultaneous prediction of structural, conformational, vibrational, and thermophysical properties for a broad range of molecules in isolation and in condensed phases and under a wide range of conditions of temperature and pressure. The COMPASS force field is under continuous review and is updated frequently. The most recent version is denoted as COMPASS on the user interface dialogs. (e) COMPASS II COMPASS II is a significant extension to the COMPASS force field. COMPASS II extends the existing coverage of COMPASS to include a significantly larger number of drugs and compounds of interest to researchers investigating ionic liquids. COMPASS II was developed in collaboration with Prof. Huai Sun of Jiao Tong University, Shanghai. Using tools developed by Prof. Sun the existing COMPASS force field was combined with existing quantum mechanical calculations to obtain new parameters and to resolve a small number of inconsistencies in the existing COMPASS force field. New valance parameters were derived by comparison of equilibrium structures and vibrational frequencies. The training set used to determine these parameters was extended as follows: (i) A polymer database consisting of 430 homo- and copolymers was analyzed using COMPASS, and 454 missing terms were found, represented by 105 fragments. These species were added to the training set. (ii) The NIST ionic liquid molecule database was analyzed in a similar way, and 40 species were added to the training set. (iii) Maybridge database consisting of 59,465 molecules was scanned using COMPASS, and 1257 fragments were added to the training set. (f) Forcite Forcite is a molecular mechanics module for potential energy and geometry optimization calculations of arbitrary molecular and periodic systems using classical mechanics. Forcite offers support for the COMPASS, UFF, and DREIDING force fields. With this wide range of force fields, Forcite can handle essentially any material. The geometry optimization algorithm offers steepest descent, conjugate gradient, and quasi-Newton methods, in addition to the Smart algorithm, which uses these methods sequentially. This allows very accurate energy minimizations to be performed. The Forcite module allows us to perform a wide range of molecular mechanic calculations on both molecular and periodic systems using classical force field-based simulation techniques. Forcite can perform many different tasks such as singlepoint energy calculation, geometry optimization, molecular dynamics, quench dynamics, anneal dynamics, shear, confined shear, cohesive energy density calculation, mechanical property calculation, and solvation free energy calculation.

42

Chapter 2

The Forcite energy task allows us to calculate the total energy of the specified system and report contributions to the energy from different force field terms, such as bond and van der Waals terms. The energy task is useful for assessing whether a force field is applicable to a given system. It can also indicate if the structure has a reasonable geometry and, hence, if a geometry optimization should be performed. The Forcite geometry optimization task allows us to refine the geometry of a structure until it satisfies certain specified criteria. This is done using an iterative process in which the atomic coordinates, and possibly the cell parameters, are adjusted until the total energy of the structure is minimized. In general, therefore, the optimized structure corresponds to a minimum in the potential energy surface. Forcite geometry optimization is based on reducing the magnitude of calculated forces and (where appropriate) stresses until they become smaller than defined convergence tolerances. It is also possible to specify an external model to represent the behavior of the system under tension or compression. The forces on an atom are calculated from the potential energy expression and will, therefore, depend on the force field that is selected. For crystal structures, Forcite geometry optimization honors any symmetry elements defined for the system. The optimization therefore amounts to a constrained optimization with a reduced number of degrees of freedom. The Forcite dynamics task allows you to simulate how the atoms in a structure will move under the influence of computed forces. Before performing a Forcite dynamics calculation a thermodynamic ensemble should be selected, and the associated parameters such as the simulation time step and the simulation temperature should be chosen. It is often difficult to perform experiments for fixed NVE conditions, so making comparisons between experiment and simulation is difficult. Application of classical Legendre transforms allows alternative ensembles to be derived, an example being the canonical, or NVT, ensemble, where a system can exchange heat with the environment and so remain at constant temperature (T). Such an ensemble can be simulated by altering Newton’s equations via the application of a thermostat. Forcite offers a choice of five thermostats: velocity scale, Nos
e, Nos
e-Hoover-Langevin, Andersen, and Berendsen. The time step is an important parameter in an integration algorithm. To make the best use of computation time, a large time step should be used. However, if the time step is too large, it may lead to instability and inaccuracy in the integration process. Typically, this is manifested as a systematic drift in the constant of motion, but it can also lead to the job failing unexpectedly due to a large energy deviation between steps. The Forcite mechanical property task allows us to calculate mechanical properties for a single structure or a trajectory of structures. A Forcite mechanical property calculation may be performed on either a single structure or a series of structures generated, for example, by a dynamics run and stored in a trajectory file (.arc, .his, .trj, and .xtd). The mechanical properties are then calculated in the succeeding text, averaged over all valid configurations, and reported in the output text document.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 43

2.1.2 Simulation Strategy Here, we have explained the general simulation strategy for a CNT-reinforced polymer composite using BIOVIA Materials Studio. The same strategy can be extended for modeling other types of fiber-reinforced composites. The general simulation strategy used in BIOVIA Materials Studio is highlighted in the succeeding text: 1. The CNT is built by using the “Build” nanostructure tool in MS 7.0. Here, we can specify the configuration of the CNT, that is, whether it is armchair (n,n) or zigzag (n,0) or chiral (n,m). 2. Then using “Build symmetry” option under “Build” tool, we increase the length of the CNT to the desired value by entering a numerical value in the options provided. For example, when we click the build symmetry option, we will be asked to enter three values of a, b, and c. Enter the value in “c,” and press enter. 3. A CNT will appear on the screen. 4. To make a polymer, we first model a monomer on the software. This can be done as shown in the succeeding text. 5. Creating and using repeat units. The repeat unit represents the structure of the monomer after it is chemically bonded into a polymer. In BIOVIA Materials Studio the head of the repeat unit is highlighted with a cyan wireframe, and the tail is highlighted with a magenta wireframe. When the polymer is constructed, terminal dangling bonds are replaced with the original atoms (usually hydrogens) or terminal repeat units. 6. Fragments and repeat units. A repeat unit can be created from any fragment that is either sketched or imported into BIOVIA Materials Studio. A repeat unit is different from a fragment in that it has a defined head, tail, and backbone. Initiators and terminators are special types of repeat units. They have a head atom and a backbone defined but no designated tail. The initiators and terminators will bond to the existing polymer chain at the endpoints, removing their head atoms in the process. To create a repeat unit from a fragment: (a) Sketch a fragment or open a 3-D atomistic document containing one. (b) Choose Build j Build Polymers j Repeat Unit from the menu bar to open the Repeat Unit dialog. (c) Select the desired head atom and click the Head Atom button. A cyan wireframe is displayed around the atom to indicate that it is the head atom. (d) Select a different atom and click the Tail Atom button. A magenta wireframe is displayed around the atom to indicate that it is the tail atom. A set of atoms are now selected; this is the automatically defined backbone. Whenever you define a head and a tail, the backbone is automatically defined as the shortest connecting path between the two atoms. To modify

44

Chapter 2

the backbone path, select a set of atoms, and click the Set Backbone Atoms button. The Clear Backbone Atoms button removes all of the selected atoms from the backbone definition. (e) Select a pseudo-chiral backbone atom, with different nonbackbone substituents, and click the Chiral Center button to assign the atom as a chiral center. Chiral inversion and tacticity pertain to repeat units that contain pseudo-chiral centers. 7. Existing and sketched repeat units. Simple polymer chains can be constructed from repeat units stored in repeat unit libraries or from sketched repeat units. To use an existing repeat unit to build a polymer chain: (a) Choose Build j Build Polymers j Homopolymer from the menu bar to open the Homopolymer dialog. (b) Select a repeat unit library from the Library dropdown list. (c) Select a repeat unit from the Repeat unit dropdown list. (d) Set other construction parameters on the Polymerize and Advanced tabs. (e) Click the Build button to create the polymer chain. The new chain will be displayed in a new 3D Atomistic Document. To use a user-defined repeat unit to build a polymer chain: (a) With the desired repeat unit displayed in the 3D Viewer, choose Build j Build Polymers j Homopolymer from the menu bar to open the Homopolymer dialog. (b) Select the Current project entry in the Library dropdown list. This indicates that you want to use a repeat unit that is in your project, rather than one of the repeat units in the set of libraries that is provided. (c) Select the name of the document that contains your repeat unit in the Repeat unit dropdown list. If the document you have selected does not contain a valid repeat unit, you will be notified. (d) Set any construction parameters as earlier. (e) Click the Build button to construct the polymer chain. The constructed chain will be displayed in a new 3D Atomistic Document. 8. After modeling the polymer and CNT, the next step is to pack the polymer around the CNT. This is done using the “Amorphous Cell” module. Under “Amorphous Cell,” click on the “packing” option, keeping the “quality” as “fine.” Specify the required density. Specify the “force field” under “energy” option. Specify the file name of the polymer that is to be packed outside the CNT by selecting from the pull-down menu under “Option” in “Amorphous Cell.” Finally, click “Run.” At the end of this step, we will get a periodic cell in which polymer will be packed around the CNT.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 45 9. Next step is performing the geometry optimization using the “Forcite module.” There is “More” option in the Forcite module. On clicking this option, we can specify the type of algorithm that is to be used for optimization, and we also specify the number of iterations for which the simulation will run. Parameters for optimization of the geometry have been listed in every chapter. Run the optimization till the energy is minimized. 10. After geometry optimization, next step is “Dynamics” run. Again go to the “Forcite” module and select “dynamics” run instead of geometry optimization. Clicking the “More” option will open a new window in which we specify the ensemble, velocities, temperature, time step, and total number of steps. Select the quality as “fine” and click “run.” After the completion of this step, we can view the energy plot and the temperature plot to see whether the system has been stabilized or not. 11. Next step is calculating the mechanical properties. Under “Forcite” module, select “mechanical property” task. Specify the number of steps and maximum strain amplitude. Here, we can again specify whether we want to optimize our structure or not. Click “Run” to begin the task of mechanical property calculation. At the end of this step, we will get a stiffness matrix and the values of Young’s moduli. The whole procedure can be represented as a flow diagram shown in Fig. 2.1.1. In molecular dynamics, mutual atomic interactions are described by force potentials associated with bonding and nonbonding phenomena. The interatomic potential energy is the sum of bonding energy and nonbonding energy: U ¼ Ubonding + Unonbonding

(2.1.1)

For CNTs, the nonbonding term is mainly the energy of van der Waals force, which normally has a weak influence on the overall mechanical behavior among the atomic interactions of the carbon microstructure. The dominant part of the total potential energy, the bonding energy, is a

Fig. 2.1.1 Flowchart of simulation strategy.

46

Chapter 2

sum of three different interactions among atoms: bond stretching, bond bending, and bond torsion: Ubonded ¼ Ubondstretch + Uanglebend + Utorsion

(2.1.2)

Here, Ubondstretch ¼

Xh

k2 ðb  b0 Þ2 + k3 ðb  b0 Þ3 + k4 ðb  b0 Þ4

i (2.1.3)

b

Uanglebend ¼

Xh

k2 ðθ  θ0 Þ2 + k3 ðθ  θ0 Þ3 + k4 ðθ  θ0 Þ4

i (2.1.4)

θ

Utorsion ¼

Xh

i k1 ð1  cos∅Þ + k2 ð1  cos 2∅Þ + k3 ð1  cos 3∅Þ

(2.1.5)

∅

where k1, k2, k3, and k4 are the force constants determined experimentally; b and θ are the bond length and bond angle after stretching and bending, respectively; b0 and θ0 are the equilibrium bond length and equilibrium bond angle, respectively; and φ is the bond torsion angle. The elastic moduli are calculated by directly computing the average mechanical forces developed between carbon atoms in the nanotube. The effective elastic moduli determination using the force approach can be calculated directly from the virial theorem given by Swenson [2] in which the expression of the stress tensor in a macroscopic system is given as the function of atom coordinates and interatomic forces. The method provides a continuum measure of the internal mechanical interactions between atoms. In an atomistic calculation, the internal stress tensor can be obtained using the so-called virial expression given by Swenson [2], which is as shown in the succeeding text: 1 σ¼ V0

"

n X

!# ! X  + mi vi vTi rij fijT 

(2.1.6)

i
i¼1

where index i runs over all particles 1 through N; mi, vi, and fi denote the mass, velocity, and force acting on particle i; and V0 denotes the (undeformed) system volume. The application of stress to a body results in a change in the relative positions of particles within the body, expressed quantitatively via the strain tensor: 2

3 ε11 ε12 ε13 εij ¼ 4 ε21 ε22 ε23 5 ε31 ε32 ε33

(2.1.7)

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 47 The elastic stiffness coefficients, relating the various components of stress and strain, are defined by   ∂σ lm  1 ∂2 A  ¼ (2.1.8) Clmnk ¼ ∂εnk T , znk V0 ∂εlm ∂εnk T , zlm , znk where A denotes the Helmholtz free energy. For small deformations, the relationship between the stresses and strains may be expressed in terms of a generalized Hooke’s law: σ lm ¼ Clmnk εnk

(2.1.9)

To calculate the axial Young’s modulus (E11) the atoms are displaced by u1 ¼ ε011. The average strain and stresses are ε11 ¼ ε011 σ 11 6¼ 0,any other σ ij ¼ 0

(2.1.10)

Through MD simulation the longitudinal elastic modulus, E11, can be calculated as E11 ¼

σ 11 ε011

(2.1.11)

In other simulation runs, the load was applied either in the transverse or shear direction. On the lateral surfaces of the simulation box, stress-free conditions were imposed to satisfy the simple tension condition. This approach is used in MD for predicting the moduli in transverse direction and was also used by Sun and Cho [3]. The internal stress tensor is then obtained from the analytically calculated virial and used to obtain estimates of the six columns of the elastic stiffness coefficient matrix. The engineering constants may be calculated from elastic constants using the following relations given by Christensen [4]:   C212 ðC22 + C23 Þ + C23 C11 C23 + C212 E33 ¼ E22 ¼ C22 + C11 C22  C212 ν12 ¼ ν13 ¼

C12 C22 + C23

1 G23 ¼ ðC22  C23 Þ 2

(2.1.12)

1 K23 ¼ ðC22 + C23 Þ 2 If the strain is only applied in uniaxial direction, for example, in longitudinal direction, the constant E11 is given as E11 ¼

σ 11 ε11

(2.1.13)

48

Chapter 2

This represents the ratio of the longitudinal stress to the corresponding longitudinal strain. Similarly, for transverse loading, the modulus E22 is given as E22 ¼

σ 22 ε22

(2.1.14)

Similarly, the shear moduli can be calculated.

2.1.3 Case Studies For the benefit of the users, in this section, various case studies have been discussed that will be of immense help to the researchers who are new to BIOVIA Materials Studio: (a) Carbon nanotubes: Different types of CNTs can be modeled in BIOVIA Materials Studio such as armchair, zigzag, and chiral. Both single-walled carbon nanotubes (SWCNTs) and multiwalled carbon nanotubes (MWCNTs) can be modeled using the “Build” module of BIOVIA Materials Studio. In this module, there are several options, and one has to select the “Nanostructure” dialog for modeling of SWCNT, MWCNT, nanorope, and nanocluster. Fig. 2.1.2 shows the interface for modeling of CNTs. By changing the values of the indexes “m” and “n,” one can obtain all types of CNTs. Also, the length of the CNT can be changed by removing the periodicity, as shown in Fig. 2.1.2. Fig. 2.1.3 shows all the three types of SWCNTs, which can be obtained from BIOVIA Materials Studio.

Fig. 2.1.2 Interface for modeling different types of CNTs.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 49

Fig. 2.1.3 ˚ and l ¼ 24.60 A ˚ ), (B) (6,0) zigzag SWCNT (d ¼ 4.70 A ˚ (A) (6,6) Armchair SWCNT (d ¼ 8.14 A ˚ ˚ ˚ and l ¼ 42.60 A), and (C) (6,3) chiral SWCNT (d ¼ 6.21 A and l ¼ 112.71 A).

(b) Nanoparticles: Besides CNTs, nanoparticles can also be modeled using the “Build” module of BIOVIA Materials Studio. For example, the silica nanoparticle (SiO2) can be modeled using the steps shown in the succeeding text: (i) (ii) (iii) (iv)

Select “File” option and then go to “structures.” In “glasses” option, select the SiO2_21A_3d file to import in the current project. A structure as shown in Fig. 2.1.4A will be obtained. Using the “Build nanocluster” dialog as shown in Fig. 2.1.4B the final structure of the silica nanoparticle can be obtained, as has been shown in Fig. 2.1.4C.

Using the same approach, any type of nanoparticle can be modeled in BIOVIA Materials Studio. (c) Nanofiber: There is no built-in module in BIOVIA Materials Studio for modeling of nanofibers. Here, we present an approach through which we can model a carbon nanofiber in BIOVIA Materials Studio. Long sheets of plane carbon atoms are drawn in BIOVIA Materials Studio 7.0, using partial double bonds to join each other forming hexagonal rings to the adjacent atoms. Fig. 2.1.5 shows one of the plane sheets. An imaginary wedge of 60

50

Chapter 2

Fig. 2.1.4 (A) SiO2_21A_3d model, (B) BIOVIA Materials Studio interface for converting SiO2_21A_3d model into silica nanoparticle, and (C) final model of silica nanoparticle.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 51 degrees angle is drawn on the modeled plane grapheme sheet, replacing all the carbon atoms falling in it with a single partial double bond. The stable configuration of this molecular model is achieved by performing energy minimization using steepest descent algorithm available in “geometry optimization” menu in “Forcite” module of BIOVIA Materials Studio. After approximately 5000 iterations of geometry optimization, the plane sheet is transformed into a cone as shown in Fig. 2.1.6. The cone thus formed is then taken as the base for modeling of CNFs, creating the second, third, fourth, and other cones by following the same procedure. The cones so obtained are then cut and stacked one over the other having an offset distance ˚ , which is close to the value of 3.4 A ˚ , reported by Zhu et al. [5]. The between them of 3.5 A schematic in Fig. 2.1.7 displays six independent cones stacked over each other having length ˚ and outer and inner diameter of 15.324 and 6 A ˚ , respectively. L ¼ 27 A (d) CNT-polymer composite: Here, we discuss an example of how to model a CNT-reinforced polymer composite using BIOVIA Materials Studio. MD simulations of polymer/CNT composites, composed of (10,10) SWCNT into two different amorphous polymer matrices, polymethyl methacrylate (PMMA) and poly-m-phenylene vinylene (PmPV), have been

Fig. 2.1.5 Plane sheet of carbon atoms.

Fig. 2.1.6 Carbon nanocone after energy minimization.

52

Chapter 2

Fig. 2.1.7 Molecular model of stacked cones.

Fig. 2.1.8 An MD model of methyl methacrylate (MMA) monomer.

performed with different volume fractions, using BIOVIA Materials Studio 7.0. MD simulations are useful to study the time evolution behavior of systems in a variety of states where thermal sampling of configurational space is required. After equilibration at finite temperature, an energy minimization method has been applied to calculate the elastic moduli of the model structures computed from the MD simulations. Periodic boundary conditions have been applied to the models along both the tube axis and transverse directions. Using the “Amorphous” module of BIOVIA Materials Studio 7.0, two chains of PMMA, each with 50 repeat units, have been built. Fig. 2.1.8 shows a methyl methacrylate (MMA) monomer. These units have been inserted in a periodic box with initial density of 0.1 g/cm3 using the DREIDING force field parameters. The model so constructed is then put into an ensemble with constant number of atoms, pressure, and temperature (NPT) simulation with a pressure of 10 atm and a temperature of 500 K for 1 ns with a simulation time step of 1 fs. Through this approach, the structure has been compressed slowly to generate initial amorphous matrix with the correct density and low residual stresses. It has been observed that further increase in pressure is not possible using the DREIDING force field. So, we have used the COMPASS force field under the same conditions for another 500 ps. The simulation temperature was then set to 298 K and pressure of 1 atm for another 500 ps. The density of the final matrix was 1.15 g/cm3, which is very close to the experimental value of

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 53

Fig. 2.1.9 An MD model of phenylene vinylene monomer.

1.19 g/cm3. The total number of atoms in the model is 11,104. MD simulation was performing approximately 50000 steps when the temperature attained a constant value. After this, a constant strain minimization has been carried out with a strain of 0.005. Similar procedure has been repeated to obtain the model of PmPV. A phenylene vinylene monomer is shown in Fig. 2.1.9. Using the procedure described in Section 2.1.2, Young’s modulus of a hypothetical isotropic amorphous PMMA matrix from the present simulations is about 2.8 GPa (compared with an experimental range of between 2.24 and 3.8 GPa). The same process has been applied to PmPV, and Young’s modulus calculated is 1.9 GPa. A (10,10) SWCNT was then placed in the center of a periodic simulation cell. PmPV molecules with different number of repeat units were then placed randomly around the tube in nonoverlapping positions. The total number of atoms in the cell varies from 15,104 to 18,904. Fig. 2.1.10 shows an MD simulation model of a (10,10) SWCNT wrapped helically with PmPV in vacuo after equilibration. Similarly, Fig. 2.1.11 shows an MD simulation model of a (10,10) SWCNT wrapped helically with PMMA in vacuo after equilibration. The mechanical properties can then be calculated using the approach as discussed in Section 2.1.2. (e) CNT metal matrix composite: As a case study, an example of CNT-reinforced titanium (Ti) composite has been discussed here. Ti has hexagonal closed packed (HCP) structure having a density of 4.5 g/cm3 [6]. The MD model of pure Ti with the specified density was created using the “Build” module of BIOVIA Materials Studio with the cell dimensions of ˚ 3. This is shown in Fig. 2.1.12. After modeling the pure Ti, a single38.7  30.2  37.7 A ˚ and having the indexes (6,0) was inserted into the Ti matrix. walled CNT of diameter 4.7 A The length of the CNT was kept equal to the length of the simulation cell so as to have a continuous periodic SWCNT in a Ti matrix. One such CNT-Ti composite is shown in Fig. 2.1.13. The number of atoms in the CNT-Ti composite model having volume fraction (Vf) of 2% was 10,283. A number of cells were created with Vf varying from 0% to 20%.

54

Chapter 2

Fig. 2.1.10 An MD model of (10,10) SWCNT wrapped helically with PmPV.

Fig. 2.1.11 An MD model of (10,10) SWCNT wrapped helically with PMMA.

Fig. 2.1.12 ˚ 3. MD model of pure Ti with the cell dimensions of 38.7  30.2  37.7 A

After constructing the models, geometry optimization was performed using the “smart” algorithm with the energy convergence criterion of 1  104 kcal/mol. The “universal” force field was used in the simulations, and the maximum number of steps for geometry optimization was around 10,000. This was followed by dynamics run during which the composite models

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 55

Fig. 2.1.13 ˚ inserted into the Ti matrix. MD simulation cell with a (6,0) single-walled CNT of diameter 4.7 A

were subjected to NVT ensemble. The time step was taken as 1 fs with the total number of time steps as 50,000. Andersen thermostat was used to stabilize the temperature of the system. The last step was performing the mechanical property calculation. In this step, the MD models were subjected to a constant strain. The strain was varied from 0.1% to 10% for studying the effect of strain on the mechanical properties of CNT-Ti composites. In all the simulations, periodic boundary conditions (PBCs) were used. The elastic moduli were calculated by directly computing the average mechanical forces developed between carbon atoms in the nanotube. The effective elastic moduli could be calculated directly from the virial theorem given by Swenson [2] in which the expression of the stress tensor in a macroscopic system was given as the function of atomic coordinates and interatomic forces. The details of the mechanical property calculation using MD could be found in another work by the authors Sharma et al. [7].

References [1] Dassault Syste`mes BIOVIA, Materials Studio, 7.0, Dassault Syste`mes, San Diego, 2017. [2] R.J. Swenson, Comments on virial theorems for bounded systems, Am. J. Phys. 51 (1983) 940–942. [3] C.T. Sun, J. Cho, A molecular dynamics simulation study of inclusion size effect on polymeric nanocomposites, Comput. Mater. Sci. 41 (2007) 54–62. [4] R. Christensen, Mechanics of Composite Materials, Krieger Publishing Company, Malbar, FL, 1991, p. 74. [5] Y.A. Zhu, Z.J. Sui, T.J. Zhao, Y.C. Dai, Z.M. Cheng, W.K. Yuan, Modeling of fishbone-type carbon nanofibers: a theoretical study, Carbon 43 (8) (2005) 1694–1699. [6] Y. Shimizu, S. Miki, T. Soga, et al., Multi-walled carbon nanotube reinforced magnesium alloy composites, Scr. Mater. 58 (2008) 267–270. [7] S. Sharma, R. Chandra, P. Kumar, N. Kumar, Effect of Stone-Wales and vacancy defects on elastic moduli of carbon nanotubes and their composites using molecular dynamics simulation, Comput. Mater. Sci. 86 (4) (2014) 1–8.

56

Chapter 2

Chapter 2.2

Overview of LAMMPS S.P. Singh*, Apurba Mandal† *

Department of Mechanical Engineering, IIT Delhi, New Delhi, India, †Department of Mechanical Engineering, NIT Uttarakhand, Srinagar, India

2.2.1 Introduction to LAMMPS The acronym LAMMPS stands for Large-scale Atomic/Molecular Massively Parallel Simulator. It is a classical molecular dynamics simulation package especially designed to run efficiently in parallel mode with only a few particles up to millions or billions [1,2]. The package is distributed as open source under the term of GNU General Public License (GPL) by Sandia National Laboratories, which is a lab specializing in defense systems and developed as a department of energy (DOE) facility.

2.2.2 Anatomy of a Nanomechanical System As we know, any molecular system or a nanomechanical system is basically a collection of atoms of different types connected by bonds and having some specific arrangement in space. If the arrangement is regular and periodic, we call it crystal lattice, while if it is nonorganized and nonrepeating, we call it amorphous. Thus sodium chloride is a crystal, steel is crystalline, while plastics and rubber are noncrystalline in general. In addition to these solids, one can also quite easily simulate systems if they are in liquid form.

2.2.3 Internal Working of LAMMPS Calculations In a physical system the bunch of atoms will be held together by forces that may be due to different kinds of interaction potentials, namely, coulomb interactions, van der Waals interactions, and electron cloud and nuclei interactions. LAMMPS normally integrates Newton’s equations of motions by velocity Verlet integrator for collections of atoms, molecules, and macroscopic particles. There are some other integrators available in LAMMPS, for example, Brownian dynamics, rigid body integrator, energy minimization via conjugate gradient or steepest descent relaxation, reversible reference system propagation algorithm (r-RESPA), and hierarchical time stepping. The interaction via short- or long-range forces is

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 57 available in LAMMPS, for example, Lennard-Jones, Buckingham, Morse, embedded atom method (EAM), and consistent valence force field (CVFF). Variety of additional boundary conditions and constraints are needed to be applied depending on the configuration of the system.

2.2.4 Methodology of MD Simulation Using LAMMPS The methodology involves defining the system configuration and boundary conditions, the environmental conditions to be maintained during simulations. The atoms in the system are till now given only geometric conditions based on lattice geometry. These atoms are relaxed to move to those positions where they can exhibit minimum potential energy of the system. After equilibration and defining of loading and boundary conditions the simulation is executed. During the execution as the atomic coordinates evolve with time the positions are recorded in a trajectory file. From the time evolution of the trajectory, the nanocomposite or whatever the system under study can be analyzed for its mechanical properties. The user needs to mention which statistical variables of the ensemble would remain constant during the simulation. Such fixes are introduced in the simulation at relevant stages and remain in force until they are disabled. The various steps have been described in the flowchart shown in Fig. 2.2.1.

2.2.5 Development of the Unit Cell Model of Polymeric Nanocomposite For modeling of a composite, we first need to create a unit cell of the composite. The unit cell consists of a reinforcement (e.g., a CNT) embedded in the matrix material. A very important part of the model involves the nature of boundary between the matrix and the reinforcement. In general, CNT does not distinguish the carbon of polymers any differently from that of its own polymer molecules. Both do exhibit van der Waals interactions. A polymer chain is free to mold itself and near a CNT, and then it will try to wind over and cling around the CNT as shown in Fig. 2.2.2. To ensure that while the van der Waals interaction is encouraged, but the tube is not allowed to disintegrate the CNT (as actually happens in practice), a cutoff distance between the two has to be maintained while doing the equilibration process. The polymer matrix modeling involves broadly two techniques, a united atom model and the full atomistic model. In united atom model, a large pendant atom is considered as one equivalent atom (united atom). The model can be created in commercial software such as BIOVIA Materials Studio [3] or open-source software such as Aten [4] and Avogadro [5]. BIOVIA Materials Studio is a complete modeling and simulation environment. The realistic atomistic model of CNT and polymer chain is easily designed with the help of Build Polymer and Build Nanostructure menus under the Build menu bar in BIOVIA Materials Studio as shown in Fig. 2.2.3.

58

Chapter 2 Start simulation

A

Define the systems (atomic positions, lattice coordinates, connectivity)

Establish the boundary conditions of the system (periodic or shrink wrap or fixed or free)

Record the new positions. As per new positions and time, calculate the new interaction forces and external loads

Restart the iteration loop

Equilibrate at specific temperature and pressure (Fix the environmental variables. These conditions on the ensemble would be checked from time to time during the simulation)

Start the iteration loop

Plot the trajectory. From the evaluted variables, estimate specific material properties of interest

Initial position and velocities of atoms and molecules apply loads, boundary conditions, Use integral of Newton’s law, d2x/dt2 = S F to estimate the new positions and velocities

At time t+dt, Recheck the new velocities for environmental (temperatures and pressures) conditions of the ensemble. If needed, scale them according to predetermined logic (thermostats and barostats)

A

Fig. 2.2.1 Flowchart of molecular dynamics procedure using LAMMPS.

The amorphous cell module provides a comprehensive set of tools to perform atomistic simulations on complex systems containing dense amorphous polymers, liquids, and other noncrystalline materials. After creating the monomer, the polymeric index is to be provided. After making a polymer chain, the number of chains is to be packed in a

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 59

Fig. 2.2.2 Wrapping of polymer around CNT.

Fig. 2.2.3 Single-walled CNT and single polymer chain of polypropylene (PP).

simulation box. It can be smoothly performed with modules of amorphous cell in BIOVIA Materials Studio. The amorphous cell module is utilized as a comprehensive set of tools to perform atomistic simulations on complex systems, for example, dense amorphous polymer, liquid, and other noncrystalline material. The simulation box and the number of chains are so related that after packing, we get the required material density of the polymer as is practically observed as shown in Fig. 2.2.4.

60

Chapter 2

Fig. 2.2.4 Unit cell of polymer nanocomposite with centrally placed CNT.

2.2.6 Setting the Conditions of Simulation To study the deformation mechanisms during uniaxial tensile deformation of an amorphous polymer, a full atomistic model is used. Interatomic force field of polypropylene is modeled with CVFF potential [6]. CVFF force field defines atom types for most hydrocarbons and many others organic molecules. The typical input script of LAMMPS contains four parts: i. Initialization: Parameters need to define before atoms of molecules are created or read from standard data file. ii. Atom definition: Atoms can be defined in three ways: lattice command, from data file, and read_restart command. iii. Settings: The main controlling section in the input file for simulation. Here, one can set the force field coefficients, simulation parameters and output options, etc. iv. Run a simulation: Dynamic simulation is operated by run command, and energy minimization (molecular statics) is performed using minimized command. v. Postsimulation processing: The standard outputs from the complete running of LAMMPS simulation are basically three kinds: restart file, thermodynamic output on the screen and into log file, and finally the trajectories of simulation that are stored in dump file. These are analyzed using specialized analysis tools. Under the postprocessing of LAMMPS simulation, standard thermodynamic output can be processed to quantify the stress-strain behavior, structural properties, etc., as has been discussed in the next section.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 61

2.2.7 Structural Properties Glass transition temperature (Tg) is the gradual and reversible transition in amorphous materials, from a hard and relatively brittle “glassy” state to a viscous or rubbery state as the temperature is increased. The molecular chains do not have enough energy present to allow them to move around. When heat is applied, the polymer molecules gain some energy, and they can start to move around. At some point the heat energy is enough to change the amorphous rigid structure to a flexible structure. The polymer molecules move freely around each other. This transition point is called the glass transition temperature. As the temperature is increased and once it reaches near the glass transition temperature of PP, the polymer chains start moving around the free volume and the conformational rearrangement of chain segment alters. As a result the polymer became less stiff at higher temperature. Also it reduces the yield point. Over the range of temperature around the Tg the volumetric data of polymer are recorded as shown in Fig. 2.2.5. In this plot the change of slope is the indication of glass transition point; as per the date obtained the Tg is about 250 K. The glass transition of polymer nanocomposite (PNC) has also been measured as shown in Fig. 2.2.6 in the same way as mentioned for PP. It is noticed that the Tg of PNC is not a fixed value like pure PP (Fig. 2.2.5). As the volume fraction of CNT is increased, the Tg of PNC is increased. This is due to the fact that the increased nanofiller normally restricted the chain segment mobility as reported in the literature [7,8].

2.2.8 Stress-Strain Behavior Parallel molecular dynamic code, LAMMPS, was used for studying the deformation of the PNC unit cell with different thermostat. This thermostat is used to regulate the temperature of the system under canonical ensemble; it is basically a system coupled with heat bath (thermostat).

Fig. 2.2.5 Volume evolution as a function of temperature for PP.

62

Chapter 2

Fig. 2.2.6 Volume evolution as a function of temperature for PNC.

Fig. 2.2.7 Stress-strain plot of PNC.

The deformation under a uniaxial tensile strain is applied at a constant strain rate with a zeropressure condition for the two lateral simulation cell faces. This deformation condition was implemented in LAMMPS by decoupling the boundary in the loading direction from the NPT equations of motion as proposed by Hossain et al. [9]. The small simulations (lesser degree of polymerization (DOP), e.g., 60, 200, and 1000, and total simulation time) have captured both the static and dynamic behavior of PNC though longer simulation time, and realistic DOP of polymer is expected to move closer to experimental results. Due to this fact, the results are fluctuating in nature. To remove these high-frequency components, a smoothing function of Matlab [10] is operated on the raw data. This “smooth” function is a local regression technique, and it uses a weighted least square and second-degree polynomial model to reduce the fluctuations. A stress-strain plot of PNC with 3.0% vol. fractions is shown in Fig. 2.2.7.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 63

2.2.9 Stress Relaxation Stress relaxation of PNC with different volume fraction of CNT has been performed with temperature below and above the glass transition temperature of pure polypropylene. Before the stress relaxation, every sample was strained within the elastic limit. Thereafter keeping that strain constant, by fixing the end faces, the bulk of the unit cell was allowed to relax for sufficient time. During these relaxation processes, it was observed how the stresses in the PNC have changed as shown in Fig. 2.2.8. In a typical stress relaxation experiment for a macrosample at room temperature the time span is about 105–107 s. For an MD simulation, that much of real-time length is not possible to simulate with currently available resources. MD simulation of stress relaxation is done in two stages. First a small strain (17.0%) is applied to the system, and then keeping that strain constant the system is allowed to relax for sufficient time (1.72 ns). With small time simulation of stress relaxation with different temperatures, the results are fitted with an empirical equation as given by Eq. (2.2.1): t

σ ¼ Ke τ

(2.2.1)

where σ is stress along z-axis, τ is the relaxation time, t is the time for stress relaxation, and K is the constant. Due to viscoelastic nature of PNC, the stress variations show approximately an exponential decay in all the case studies. These isothermal stress relaxation curves help to find the relaxation time and viscoelastic properties of nanocomposite for different vol. fractions of CNT. With exponential fitting of all isothermal curves the relaxation time is calculated as tabulated in Table 2.2.1 along with other parameters.

Fig. 2.2.8 Stress relaxation of PNC nanocomposite with 1.0% volume fraction of SWCNT.

64

Chapter 2 Table 2.2.1 Evaluated relaxation parameters for PNC

Vol. Fraction 1.0

2.5

3.0

Temp., T (K) 200 253 300 350 400 200 253 300 350 400 200 253 300 350 400

Modulus, Y (GPa) 1.21 0.98 2.62 1.10 1.05 1.44 2.18 2.65 3.62 1.33 3.07 1.09 2.27 1.64 1.13

Relax. Time, τ (s) 8

2.30e 4.26e9 1.89e9 1.36e9 1.71e9 5.75e9 2.23e9 9.32e9 2.15e9 1.21e9 6.62e8 8.33e9 6.12e9 2.88e9 7.98e9

Viscosity, η (Pa s) 27.83 4.17 4.95 1.50 1.80 8.28 4.86 24.69 7.78 1.61 203.23 9.08 13.89 4.72 9.02

2.2.10 LAMMPS Input File (a) Program initialization: The input file shown in Fig. 2.2.9 is used for analysis of a polypropylene nanocomposite unit cell when subjected to a uniaxial tension. The results were compared with a research paper by Hossein et al. [9]. As we can see the input file has number of sections. Each denoted by the symbol “#.” In the first section, we declare some system variables by the keyword “variable” and “variable name” and what that variable name will be whatever string. In the initialization routine the units correspond to a specific type of units for mass (grams/mole), length (angstroms), time (femtoseconds), pressure (atmospheres), temperature (Kelvin), and other derived quantities. In the “boundary” keyword, we need to specify the boundaries in three directions. For example, p p p states the system under consideration extends in all the three direction, and periodic boundary conditions are applicable. In the atom_style, we mention the style of atom, which essentially tells the program that attributes of the atom are going to be used during the simulation. One important parameter of MD simulation, which is to be predecided related to how far the interactions are to be considered for this purpose the keyword “neighbor,” is specified that tell how much distance from the atom under consideration the atom interactions are to be considered. After some movements of the atoms, the nearest neighbors might change. So here again the user specifies when this list of neighbor atoms are to be updated. Bond style harmonic uses to compute bond interactions between pairs of atoms. The read data keyword tells the program where to pick the initial structure. For simple structures (such as fcc, bcc, or hcp metals) the initial structure can be generated within LAMMPS with commands like create_atom and lattice and create_box. However for more complex structure the

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 65

Fig. 2.2.9 Input file used for analysis of a polypropylene nanocomposite unit cell when subjected to a uniaxial tension.

structure file (atomic positions and their bond definitions) is generated outside and then called into LAMMPS via read_data and read_restart commands. (b) Equilibration: Fig. 2.2.10 shows the input file used for equilibration of the PNC. The first important step for any simulation is to equilibrate the structure. When the initial configuration is made the atoms are not set as their minimum energy level as per the potential filed. If we start simulation on a nonequilibrated structure, we shall get wrong and absurd results since the relaxation process will take place along with the loading procedure and will affect the results, since most of the polymers are made from long-chain molecules that have rather limited mobility. Thus simple relaxation process may not lead to minimum energy configuration. For this purpose, we generally adopt a four three-step relaxation process in which we take the system first to a large temperature where the chains have high mobility, relax it at that temperature, and then bring down slowly to room temperature, or the whatever temperature, where we want to study the material behavior. The keyword “velocity” initializes the atomic velocity to set values such that the temperature of the ensemble is 500 K. Then at that temperature the system is allowed to equilibrate for 10000 iterations. In the second process the temperature is slowly reduced from 500 K to the working temperature of 300 K as can be seen by the fix command. Fix commands are used for setting the environmental variables. More details on the syntax of the commands can be obtained from the documentation manual in the official website of LAMMPS [11]. Fig. 2.2.11 shows the input file for the final stages of equilibration of the PNC. The thermo_style command allows defining the style and content for printing thermodynamic data to the screen and log file. Write restart is one useful command for creating current

66

Chapter 2

Fig. 2.2.10 A part of the input file used for equilibration of the PNC.

Fig. 2.2.11 Final part of the input file used for equilibration of the PNC.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 67

Fig. 2.2.12 An input file used for performing dynamics on the PNC.

state of simulation in binary form that can be used subsequently, if the simulation failed during the long period of running or if the simulation is to be extended. (c) Dynamics: After the equilibration stage is completed the stage is set for the next dynamic simulation of the system when subject to any user-defined input. For this purpose, we first need to know what kind of properties we are interested in calculating during the simulation run. For example, to calculate the modulus of the composite, we need to know the global average stress and the overall strain in the sample. For any strain energy-related calculations, we need to store the energies of the system. For Poisson ratio calculation, we need to know the lateral and longitudinal strain. The code given in Fig. 2.2.12 first defines the different variables, which need to be output as the simulation grows. (d) Tensile loading: The next part in the code, as shown in Fig. 2.2.13, applies a uniform strain (“erate”) at a prefixed rate. While the straining is done, the temperature of the atomic ensemble might increase, which might affect the results. Thus while the straining process is going on the temperature constant is also maintained. Accordingly, different thermostats are used. For uniaxial tensile simulation an isothermal-isobaric (NPT) ensemble is used.

2.2.11 LAMMPS Output File (a) Log file: Part of a sample log file has been shown in Fig. 2.2.14, where the initial part of the file is basically the detail of simulation condition and default setting. This file also stores the data of thermo_style command.

68

Chapter 2

Fig. 2.2.13 An input file used for tensile loading of the PNC.

Fig. 2.2.14 An output log file obtained for tensile loading of the PNC.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 69

Fig. 2.2.15 An output trajectory file obtained for tensile loading of the PNC.

(b) Trajectory file: A sample trajectory file is shown in Fig. 2.2.15 with normal store unit cell detail, the atomic coordinate at each time step as specified with dump command.

References [1] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Comput. Phys. 117 (1995) 1–19. [2] S. Plimpton, R. Pollock, M. Stevens, Particle-mesh Ewald and rRESPA for parallel molecular dynamics simulations, in: The English SIAM Conference on Parallel Processing for Scientific Computing, SIAM, 1997. pp. 1–13. [3] Dassault Syste`mes BIOVIA, Materials Studio, 7.0, Dassault Syste`mes, San Diego, 2017. [4] T.G.A. Youngs, Aten-An application for the creation, editing, and visualization of coordinates for glasses, liquids, crystals, and molecules, J. Comput. Chem. 31 (2010) 639–648. [5] M.D. Hanwell, D.E. Curtis, D.C. Lonie, T. Vandermeersch, E. Zurek, G.R. Hutchison, Avogadro: an advanced semantic chemical editor, visualization, and analysis platform, J. Cheminform. 4 (2012) 17. [6] P. Dauber-Osguthorpe, V.A. Roberts, D.J. Osguthorpe, J. Wol, M. Genest, A.T. Hagler, Structure and energetics of ligand binding to proteins: escherichia coli dihydrofolate reductase-trimethoprim, a drug-receptor system, Proteins Struct. Funct. Genet. 4 (1988) 31–47. [7] M.K. Seo, J.R. Lee, S.J. Park, Crystallization kinetics and interfacial behaviors of polypropylene composites reinforced with multi-walled carbon nanotubes, Mater. Sci. Eng. A 404 (2005) 79–84. [8] M. Ganβ, B.K. Satapathy, M. Thunga, R. Weidisch, P. Ptschke, A. Janke, Temperature dependence of creep behavior of PP-MWNT nanocomposites, Macromol. Rapid Commun. 28 (2007) 1624–1633. [9] D. Hossain, M.A. Tschopp, D.K. Ward, P. Wang, M.F. Horstemeyer, Molecular dynamics simulations of deformation mechanisms of amorphous polyethylene, Polymer 51 (2010) 6071–6083. ® [10] Matlab , Matlab—The Language of Technical Computing, The Math Works, Karnataka, India, 2010.http:// www.mathworks.in/. [11] Official Web of LAMMPS, https://lammps.sandia.gov/doc/Manual.html.

70

Chapter 2 Chapter 2.3

Overview of GROMACS Raja Sekhar Dondapati School of Mechanical Engineering, Lovely Professional University, Phagwara, India

2.3.1 Introduction GROningen MAchine for Chemical Simulation (GROMACS) is a software package, which is used to create physical systems consisting of particles, from hundreds to millions in number, and to solve Newtonian equations of motion [1]. For relatively small number of particles or molecules involved, manual methods may serve for the purpose; however, when the number of molecules in the system is large, this method is not feasible. Hence, we need to search for computational code, which could solve the governing equations of motion in relatively small timescale. GROMACS is one of the suitable codes, which could serve for enhancing our understanding on the behavior of matter at molecular level. Fig. 2.3.1 shows the utility of GROMACS for different scope in molecular dynamics study. GROMACS is one of the widely accepted open-source software code, primarily built for dynamic simulations of biomolecule, and complex motion of a system of particles [2], which provides rich set of calculation methodology, preparation, and analysis tools [3]. Apart from bonded or nonbonded interaction of molecules, GROMACS can handle rectangular periodic boundary condition, with temperature and pressure scaling [4]. Kuling et al. [5] used the GROMACS package for simulation of lipid molecules, under the effect of external force field. Topologies were provided for four phosphatidylcholines: saturated DPPC, mono-cis-unsaturated POPC and DOPC, and mono-trans-unsaturated PEPC. Dynamics and simulation studies were carried out for human α-amino-β-carboxymuconate-ε-semialdehyde decarboxylase (hACMSD), which is a zinc, containing amidohydrolase, for calculating the binding free energy of the complexes [6].

Fig. 2.3.1 Application of GROMACS in various fields of sciences.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 71 Extension to GROMACS, such as ProtoMD [7], a toolkit written in python and freely available under the GNU GPL, facilitates in the development of algorithm for molecular dynamics simulation. Sampling of the data is of prime importance in molecular dynamics simulation, and subsequently, implementation of various algorithms and techniques is proposed, like distributed sampling, enhanced sampling methods, and the use of coarse-grained (CG) models [8]. Applications in the field of nuclear physics experiment, like neutron scattering experiments, have also utilized the approach of molecular dynamics, for simulation purposes [9]. Ion channels, which are selective membrane proteins, for ions like sodium, potassium, and chloride simulations, which are of universal importance in the field of cellular physiology and pathology, have also shown earlier simulation techniques with the aid of GROMACS [10]. Additional parallel software has also been implemented to enhance the data, for further technical output, where single or distributed atomic clusters can be simulated to decrease the calculation time, by dividing the frames across the processors available [11]. Mutated enzyme models were also simulated using GROMACS, for effective drug-delivery process [12]. Shortrange interatomic interactions, which governs many biomolecular processes, are of key importance for system comprising of solute molecule, which complexly interact with solvent atoms [13]. Computational modeling of polymers, which plays fundamental role in multidiverse fields of oil industry, pharmaceutical industry, and the development of materials, is also achieved using GROMACS, which however was poorly understood at atomistic level [14]. Asymmetrical lipid bilayers were also simulated using GROMACS, whose studies are gaining attention of researchers, to identify the detailed lipid compositions in cell membranes [15]. Molecular dynamics simulations are also useful in simulating variable external fields. Vagadia et al. [16] investigated the stability of soybean trypsin inhibitor (STI) protein, under the effects of temperature and oscillating electric fields. Thermal bioflow problems were also investigated, to analyze heat and mass transfer phenomenon, which is the key requirement for the establishment of correlation between the flow fields and thermal gradient developed in the fluid domain, at the fluid-wall contact boundary layer [17]. GROMACS has also been used to simulate ion/water exchange process, in a double-membrane simulation setup, where steady-state continuous flow of ions through the channels is achieved [18]. Following the success of the simulation setup, simulations for hydrogen degrees of freedom removal, from a lipid bilayers, were also achieved [19]. Macromolecular rotor motions in atomistic simulations were also performed, which are the essential for functioning of many motor proteins [20]. Lung surfactant protein B (SP-B), a complex mixture of proteins and lipids that covers the air-water interface of the alveoli, which reduces surface tension and aids in the breathing process, is also simulated using GROMACS [21]. In this section, we present an overview of the steps involved in the functionality of GROMACS, and an introduction to computational chemistry [22], molecular dynamics [23] simulations, and energy minimization [24] is illustrated.

72

Chapter 2

2.3.2 Working Principle of GROMACS In the procedural steps for GROMACS (see Fig. 2.3.2) the initial step is the selection of a group or cluster of molecules for simulation purposes. The molecule under consideration can be predefined, or a new molecular system can be defined, based on the fact that its stability based on valence electrons must be satisfied. The next step is the creation of simulation cell on which the condition of periodicity is imposed. Therefore, the readers are suggested to take precaution while selecting the type of simulation cell. Then, the initial conditions of the coordinates and velocities of all the atoms must be set. However, if the velocity is unknown, statistical mechanics approach is utilized in initializing the velocities to the atoms. Moreover, geometry optimization of the molecular system is an important step, when predicting the velocities of the atom, to achieve the minimal energy levels of the developed system, for stability purposes. Then the atoms are coupled to one or more group, based on the requirement of the user. For example, when computing a system consisting of two immiscible phases, then the surface tension coupling must be imposed as an interfacial condition. Later, during the simulation, the update algorithm is followed, which writes the atomistic details, such as position, velocity, and acceleration as an output file.

Structure topology

Load molecule

Molecular database

Polymer database

Creation of system Simulation cell

Assigning initial conditions Coordinates and velocities

Computation of forces

Updating algorithm Position of all atom at time t 1 Δt Velocities ν of all atoms at time t – –– 2 Acceleration of all atoms at time t

Fig. 2.3.2 Overview of the fundamental steps involved in the computational code of GROMACS.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 73

2.3.3 Computational Chemistry and Molecular Modeling GROMACS perform molecular dynamic simulations and energy minimization, which are two of the many approaches utilized in the field of computational chemistry and molecular modeling. In computational chemistry, methods are based on quantum mechanics of molecules to the molecular simulation of large aggregates. However, in molecular modeling, the objective is to build realistic atomic model, for predicting macroscopic properties of the developed system. Macroscopic properties can be further classified as (a) static equilibrium properties, such as potential energy of a system [25] or radial distribution of a liquid [26] and (b) dynamic or nonequilibrium properties, such as viscosity, diffusion processes [27], phase change dynamics [28], reaction kinetics, or dynamics of defects in crystals [29]. However, the relativistic timedependent Schr€ odinger equation [30] aids in describing the properties of molecular systems, which implies that the transient simulations yield better results compared with the simulation done at equilibrium state. However, more complexities are involved when transient simulations are involved. Thus approximations are required to reduce the complexities of the developed systems. The longer the time span of the simulation, the more number of approximations is required for the reduction of terms involved in governing equations. For molecular modeling, two consequences follow: (i) It is important to generate a representative ensemble to retrieve the macroscopic properties. However, it is not enough to compute thermodynamic equilibrium properties, which are based on free energies [31], phase equilibrium [25], binding constant [32], etc. (ii) During the computation process, the atomistic details of structure and molecules are often irrelevant for accounting the macroscopic properties. This opens the possibility to simplify the atomic interactions and average the involved processes. The concepts of statistical mechanics [33] approach can be utilized to achieve the required simplification. To generate representative equilibrium ensemble, two methods are available: (a) Monte Carlo simulations [34] and (b) molecular dynamics simulations [21]. When transient case is involved, Monte Carlo simulations may fail to predict the macroscopic properties, and molecular dynamics simulations may produce significantly better results. That is, for nonequilibrium ensembles and analysis of dynamic systems, molecular dynamics approach is a more appropriate methodology. Hence in the present context, we present a discussion on molecular dynamics. However, when the initial configuration of the system is very far from equilibrium, a robust energy minimization is required. Moreover, energy minimization also helps in reducing the thermal noise of the system, by removing the kinetic energy and reducing the potential energy.

2.3.3.1 Molecular Dynamics Simulations Molecular dynamics solves Newton’s equations of motion for a system of N interacting atoms, given by Eq. (2.3.1):

74

Chapter 2 mi

∂2 ri ¼ Fi , i ¼ 1,…,N ∂t2

(2.3.1)

where the forces being the negative derivates of a potential function V (r1, r2, …, rN) defined as Fi ¼ 

∂V ∂ri

(2.3.2)

The equations are solved simultaneously in small time steps, and the temperature and pressure are monitored, to retain its required values. Furthermore, the coordinates are written to an output file at fixed intervals of time. The coordinates of the system, as a function of time, represents the trajectory of the system.

2.3.3.2 Molecular Dynamics Approximation In this section the approximations that are incorporated for the simulation intent are presented, to illustrate the limitations of the software. A constant verification with the experimental data should be undertaken to the check the accuracy of the simulation: i. The simulations are classical Newtonian equations are incorporated to describe the motion of atoms. However, this approach is useful at normal temperature, and deviation from actual results may occur at very low temperature. The quantum mechanical phenomenon, like quantum tunneling of proton through a potential barrier [35], cannot be handled using Newtonian approach. Another example is of helium at low temperature [36], where the classical approach breaks down. ii. Electrons are in ground state In molecular dynamics, the motion of electrons is not considered. The electrons are supposed to adjust their position when a change in atomic position is encountered (Bohr-Oppenheimer approximation) [37] and remain in their respective ground state. However, following the consequences, phenomenon like electron excitation or electron transfer processes cannot be taken into account. Moreover, for the same reason, the chemical reactions cannot be handled properly. iii. Force fields are approximate The force is provided by the force fields. These forces are not a part of simulation method; consequently the parameters can be changed accordingly. The existing version of GROMACS includes pair additive properties, with an exception existing for long-range coulomb forces. Moreover, polarizabilities are not incorporated, and fine-tuning of bonded interactions is not possible. iv. Long-range interactions are cut off For the Lennard-Jones interactions, a cutoff radius is always used by GROMACS. A single image of each particle in the periodic boundary condition is used by GROMACS for the “minimum image convention.” As a result the cutoff radius is never exceeded half of the simulation cell.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 75 v. Boundary conditions are unnatural The intrinsic small nature of system causes a cluster of particles to interact with the environment (here vacuum), due to which the system develops unwanted boundaries. To compensate the stated problem, periodic boundary conditions are imposed, to handle the problem of interphase interaction [38], as shown in Fig. 2.3.3. Moreover, this method is additionally beneficial when dealing with crystalline structures or compounds.

Fig. 2.3.3 Schematic representation of the procedure for MD simulations, which consist of (A) creation of simulation cell, (B) assigning of molecules, and (C) imposing constrains and periodic boundary conditions.

76

Chapter 2

However, problem arises when noncrystalline structure, such as water, is taken into account, where the condition of periodicity is irrelevant.

2.3.3.3 Energy Minimization Energy minimization [39] is a technique that is utilized to achieve lowest possible energy levels of the simulated cell. In this approach, the process is to find the possible arrangement in a space-filling volume, to check which configuration yields the lowest possible energy state. In this process, the configuration has one global minimum (one deepest point), a large number of local minima, where the derivate of the potential energy is zero and the second derivative is nonnegatives. The matrix of second derivative is called Hessian matrix [40]. In between the local minima, saddle points exist, where the Hessian matrix has only one negative eigen value. The dynamics of structural configuration and transition can be completely defined when one has the knowledge of entire local minima, global, and saddle points. However, GROMACS can calculate the nearest local minimum only. If one wishes to find the other minima, eventually the global minimum, then it is advisable to use the temperature-coupled molecular simulations [41]. In this process, the simulated system is run at high temperature and slowly reducing the temperature to the predetermined temperature, a process called simulated annealing [42]. Three possible energy minimization methods are as follows: i. Methods that require only function evaluations. Examples are simplex method [43] and its variant. ii. Use of derivative information: The partial derivatives of the potential energy with respect to spatial coordinates are known in MD simulations. Hence, this method could also be employed for energy minimization process. iii. Methods that also require second derivate information. These methods are superior in convergence criteria; however, extra computational cost is required to obtain the second derivatives [33]. This method could be beyond the available capacity of most systems.

2.3.4 Algorithms In GROMACS software, generalized concepts of periodic boundary conditions and group concept are utilized for molecular dynamics. In the present section the aspects of the algorithms are described, along with pair-list generation, velocity and position update, external pressure, and temperature coupling.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 77

2.3.4.1 Periodic Boundary Conditions The application of periodic boundary condition helps in minimizing the edge effects. Under this scheme, atoms are packed into space-filling box, which is surrounded by translated copies of itself, which imposes no boundaries condition. If this condition is not imposed, then the simulated molecular geometry will develop interphase with the environment (here vacuum). Physically, this situation implies that the phase will be in contact with the environment. However, for simulation criteria, the environment would be vacuum, which would demand additional boundary conditions to be applied and possibly produce erroneous results. This condition is further beneficial when the system to be analyzed is crystalline in nature. However, for nonperiodic systems, such as simulation of liquids, the condition of periodicity seems to be irrelevant in nature, due to the anisotropic distribution of molecules throughout the volume of matter phase. Hence, the condition of periodicity may cause erroneous results. GROMACS supports triclinic boxes of any shape. Considering a simulation box (unit cell), defined by three box vectors a, b, and c, the following conditions must be satisfied: ay ¼ a z ¼ bz ¼ 0

(2.3.3)

ax > 0, by > 0, cz > 0

(2.3.4)

  1 1 1 jbx j  ax , jcx j  ax , cy   by 2 2 2

(2.3.5)

2.3.4.2 The Group Concept The GROMACS MD uses user-defined groups of atoms for simulation purposes. The maximum of 256 groups are allowed; however, each atom can be subjected to six different groups: i. Temperature-coupling group The temperature coupling parameters, such as reference temperature, time constant, and number of degrees of freedom, are considered [44]. This parameter allows to define a state of temperature, where there could exist a differential temperature range. This group is essential in implementing state of the systems, where the constituent structures have different temperature magnitude. For example, temperature coupling group is beneficial in modeling systems, where the surface has to be kept at a lower temperature, than the adsorbing molecule, as shown in Fig. 2.3.4. ii. Freeze group The atoms belonging to freeze group are kept stationary during the simulation process. This methodology is particularly useful to model equilibration [45]. Moreover, to avoid

78

Chapter 2

Fig. 2.3.4 Illustration of heat surface and adsorbing molecule in thermal contact.

unreasonable movements of molecules during MD simulations, this methodology is useful. Also, an atom has six degrees of freedom, and its motion can be constrained such that the atom has the provision to move either in a particular plane or line. In other words, user-defined directional constraint can be imposed to the collections of atoms, to restrict its motion. Hence, a fully frozen atom cannot be moved by any constraint. iii. Accelerate group When external force is to be imposed to the system, “accelerate group” is useful to achieve the required acceleration of the atoms. This methodology is equivalent to the system subjected to external force condition. Each atom in an “accelerate group” is imposed with an acceleration ag condition. Hence, it allows to develop nonequilibrium simulation [46]. Moreover, this methodology is useful in obtaining transport properties, like thermal conductivity. iv. Energy-monitor group Mutual interactions among all the energy groups are compiled during the simulation process. All the nonbonded interactions between the pairs of energy-monitor group can be excluded. v. Center-of-mass group This classification of group is useful when low-viscous medium is to be modeled; an example would be a gas system. Low viscosity would refer to a state when the intermolecular distance is large, such that the collision between the atoms doesn’t occur, which could have given rise to frictional forces. Owing to this frictional forces, the center of mass motion is prevented [47]. vi. Compressed position output group To reduce the size of the trajectory file (.xtc or .tng), a subset of all particles can also be stored.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 79

2.3.5 Molecular Dynamics A global flow scheme of MD is shown in Fig. 2.3.5. Each simulation requires a set of initial coordinates. Additionally, for transient simulations, initial velocities of all the particles are also required as input parameter.

2.3.5.1 Initial Conditions i. Topology and force field

1. Input initial conditions Potential interaction V of atoms Position r of atoms Velocities n of atoms

2. Compute forces The force on atom is defined as, Fi =–

JV Jri

3. Update configuration The movement of atoms is simulated by numerically solving Newton’s equations of motion F d2ri = i mi dt2

4. Output step (if required) Write positions, velocities, energies, temperature, pressure, etc.

Fig. 2.3.5 Global molecular dynamics (MD) algorithm.

80

Chapter 2

The system topology and the description of the force field are required for simulation intent. According to the desired condition, the molecules constituting in a cell might be subjected to external force. Under such circumstances, this methodology is useful. ii. Coordinates and velocities Before the start of the simulation, the size of the cell, coordinates, and the velocities of all the particles are needed to be assigned. The box size is determined by three vectors b1, b2, and b3. If the velocities remain unknown, the program uses the initial atomic velocities vi, where i ¼ 1,…,3 N, at a given absolute temperature T, defined as rffiffiffiffiffiffiffiffiffiffiffi   mi mi v2i exp  (2.3.6) pðvi Þ ¼ 2πkT 2kT where k is Boltzmann’s constant. To retrieve the velocity, normally random distributed numbers as shown in Fig. 2.3.6 are generated by adding 12 random numbers Rk, defined as 0  Rk < 1

(2.3.7)

iii. Center-of-mass motion The center of mass is usually kept stationary at each step to carry out simulation under no external forces, that is, the velocity of the center of mass should remain constant. However, in 3

×10–6

Probability distribution (-)

2.5

2

1.5

1

0.5

0 0

1

2

3

4 5 6 Velocity (m/s)

7

8

9

10

Fig. 2.3.6 Maxwell-Boltzmann velocity distribution curve generated from random numbers.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 81 practice the algorithm introduces a small slow change in the center of mass; as a result the total kinetic energy of the system remains unfixed. This could lead to large deviation in the required simulation results in the long run, especially when the temperature parameter is coupled to the cell.

2.3.5.2 Neighbour Searching Internal forces are either generated from fixed (static) lists or from dynamic lists. The forces from dynamic lists contain nonbonded interactions between any pair of particles. When nonbonded forces are considered, it is convenient to have all particles in a rectangular box. A triclinic box can be transformed into a rectangular box. It may be necessary because output coordinates are always in a rectangular box, even when a dodecahedron or triclinic box is used for the simulation purposes.

2.3.5.3 Pair Lists Generation The nonbonded pair forces are needed to be calculated for those pairs i and j for which the distance rij is less than the given cutoff radius Rc. When the interaction between the particle pairs is fully accounted by bonded interactions, then the particles are excluded. A pair list is employed by GROMACS for which the nonbonded interactions are calculated. The pair list would consist of particles i, a displacement vector for particle i, including the particles j, that are within “rlist” of this particular image of particle i. The list is updated in every “nstlist” steps, where “nstlist” is typically 10. A provision is there to calculate the total nonbonded force on each particle, as all the particles are in a shell around the list cutoff, that is, at a distance between “rlist” and “rlistlong.” In order a make a neighbor list, all particles that are close, that is, within the neighbor list cutoff, to a given particle must be found. This methodology of searching is generally referred as neighbor search (NS) or pair search, which involves period boundary condition and determination of image. The search algorithm is O (N), along with a simpler O (N2) that is also available under some conditions.

2.3.5.4 Cut-Off Schemes: Group Versus Verlet GROMACS supports two different cutoff scheme setups, from its version 4.6 release. The original one is based on particle group and the other on Verlet buffer. Fundamental differences are present between the two schemes, which affect the results, performance, and feature support. The group scheme can almost work like the Verlet scheme; however, a decrease in performance is expected. When the simulation deals with molecules, which are commonly abundant in many simulations, such as water molecules, the group scheme is especially fast. A neighbor list is generated in the group scheme, which consists of pair of groups at one particle. Charge groups were originally present but with proper treatment of long-range order of electrostatics. These simulations were performed in unbuffered status, which was the only

82

Chapter 2

advantage. The center of geometry is within the cutoff distance when groups are put into the neighbor list. MD steps are performed until the interactions between all particles are calculated, until the neighbor list is updated. This setup is efficient, as the neighbor search only check distance between charge-group pair and not particle pairs. This setup is efficient as the neighbor search only checks the distance between the charge-pair group and not particle pairs. Without explicit buffering, this setup leads to the phenomenon of energy drift, as some particle pairs that are within the cutoff do not react and some outside the cutoff do react. Such situations occurred when (i) particles move across the cutoff neighbor search steps and (ii) when more than one particle or particle pairs move in/out of the cutoff, with their geometric center distance outside/inside of the cutoff. The intervention of user in adding a buffer to the neighbor list removes such problems; however, high computational cost is involved. The severity of such problems depends on the system, the properties in which we are interested, and the cutoff setup. A buffered pair list is present in Verlet cutoff scheme by default. However, the cluster of particles in this scheme is not static, as in the case of group scheme. To ensure that the particles move between pair search steps, forces between almost all the particles within the cutoff distance are calculated. However, the size of the applied buffer depends upon various issues, such as the implementation of temperature parameter.

2.3.5.5 Energy Drift and Pair-List Buffering When dealing with canonical (NVT) ensemble, the atomic displacement and the shape of the potential at the cutoff can be used to determine the average error caused by the finite Verlet buffer size. The displacement distribution obtained along one dimension for a particle moving freely, with mass m, in a time period t at temperature T is Gaussian with zero mean and variance σ 2 ¼ tκBT/m. The variance changes to σ 2 ¼ σ 212 ¼ tκ BT(1/m1 + 1/m2), for the distance between two particles. However, in actual practice, particles do interact with other particle over the specified time period, and hence the real displacement distribution is much narrower. When a nonbonded interaction cutoff distance of rc and a pair-list cutoff rl ¼ rc + rb are provided, then the average error after time t for pair interactions between one particle of type 1 surrounded by particles of type 2 with number density ρ2, when the particle distance changes from r0 to rt, can be expressed as ð ðrc ∞ 4πr02 ρ2 V ðrt ÞG

hΔV i ¼ 0 rl

r  r t 0 dr 0 dr t σ

(2.3.8)

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 83 ð ðrc ∞ 4πr02 ρ2

 ∞ rl

ðrc ∞ ð

2

 4π ðrl + σ Þ ρ2 ∞ rl

 r  r 1 t 0 2 0 00 V ðrc Þðrt  rc Þ + V ðrc Þ ðrt  rc Þ G dr 0 dr t σ 2

(2.3.9)

  1 1 rt  r0 dr 0 dr t V 0 ðrc Þðrt  rc Þ + V 00 ðrc Þ ðrt  rc Þ2 + V 000 ðrc Þ ðrt  rc Þ3 G σ 2 6

8 9 h r   rb i 1 0 b 2 2 > > > > E ð r Þ r σG + σ  r V c b b > > > > σ σ 2 > > > > > > > > < = h   i     1 r r 2 b b 00 2 2 2 2 ¼ 4π ðrl + σ Þ ρ2 + V ðrc Þ σ rb + 2σ G  rb rb + 3σ E > > σ σ 6 > > > > > > > > > > h   i > >     1 r r > : + V 000 ðrc Þ rb σ r2 + 5σ 2 G b  r 4 + 6r 2 σ 2 + 3σ 4 E b > ; b b b σ σ 24 (2.3.10) where G is a Gaussian distribution, along with zero mean and unit variance. It is always desired to obtain least energy error, so σ will be small compared with rc and rl. The energy error would be averaged over all particles and weighted with particle counts. However, when dealing with condensed matter phase, the estimated error will be much greater than in actual practice, as the displacement will be much smaller than for freely moving particles. When constraints are present, such as bonded molecules, some particles will have fewer degree of freedom, which will tend to reduce the energy error. Furthermore, the displacement in an arbitrary direction of a particle with two degrees of freedom is not Gaussian, rather it follows the complimentary error function:   pffiffiffi π jr j pffiffiffi erfc pffiffiffi (2.3.11) 2σ 2 2σ where σ 2 is κBT/m. However, this distribution cannot be integrated analytically to obtain the energy error. One important implementation is present that reduces the energy error caused by the finite Verlet buffer list size. The earlier derivation assumes a particle pair list. However, for more efficiency, GROMACS uses a cluster pair list. The list consists of clusters of four particles (in most cases), called a 4  4 list. Slightly beyond the pair-list cutoff, there will still be a large fraction of particle pairs present. This unknown fraction can be determined in a simulation and accurately predicted under reasonable assumptions.

2.3.5.6 Cut-Off Artifacts and Switched Interactions Using the Verlet scheme, the pair potentials are shifted to be zero at the cutoff, making the potential the integral of the force. Such situation is possible in the group scheme, if the shape of the potential corresponds to zero at the cutoff distance. Moreover, energy drift is possible when

84

Chapter 2

the forces are nonzero at the cutoff. Such effects are negligibly small and may be even skipped, under other dominant factors, like integration errors from constants. To completely avoid such cutoff artifacts, the nonbonded forces can be completely switched to zero at distances smaller than the neighbor list cutoff.

2.3.5.7 Grid Search All particles are put on a grid, with the smallest spacing Rc/2 in each of the directions. Along the direction of each box vector, a particle i has three images. Each direction of the image is associated with 1, 0, and 1, associated to a transition of each box vector. The searching method involves the construction of images first and then searching the neighbors corresponding to the image of i. On such occasion, some grid cells may be searched more than once. The grid search methodology is equally fast for rectangular and triclinic boxes.

2.3.5.8 Charge Groups For the reduction of cutoff artifacts of coulomb interactions, charge groups were introduced. When dealing with plain cutoff, significant jumps in the potential and forces arise when atoms with (partial) charges translates in and out of cutoff radius. These jumps can be reduced by moving groups of atoms with net charge zero, called charge group, in and out of the neighbor list. Such condition reduces the cutoff effects from charge-charge level to the dipoledipole level. However, under this circumstance, a slight negative result may be produced, depending on the neighbor list generated and the process of calculating interactions. Some important reasons are present for using the “charge group.” In molecular topology the neighbor searching is carried out on the basis of charge groups. If the image distance between the geometric centers of the atoms is less than the cutoff radius, all atom pairs are included in the pair list. Furthermore, the charge groups are ignored, when the Verlet cutoff scheme is implemented.

2.3.6 Compute Forces In this section, the methodology for calculation of the forces is discussed, along with the various mathematical models involved in this scheme.

2.3.6.1 Potential Energy When the forces are computed, the potential energy resulting from each interaction is computed as well. The summation of all the potential energy arises from various contributors, such as Lennard-Jones, coulomb, and bonded terms.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 85

2.3.6.2 Kinetic Energy and Temperature The temperature is calculated using the total kinetic energy of the N particle: Ekin ¼

N 1X mi v2i 2 i¼1

(2.3.12)

The absolute temperature can be computed using 1 Ndf kT ¼ Ekin 2

(2.3.13)

where k is Boltzmann’s constant and Ndf is the number of degrees of freedom, which is given by Ndf ¼ 3N  Nc  Ncom

(2.3.14)

Here the number of constraints on the system is given by Nc. When more than one temperature coupling group is used the number of degrees of freedom for group i is given as   3N  Nc  Ncom i ¼ 3N i  Nci Ndf 3N  Nc

(2.3.15)

When pressure calculation is required, the kinetic energy can also be written in the form of tensor, defined as Ekin ¼

N 1X mi vi vi 2 i

(2.3.16)

2.3.6.3 Pressure and Virial The difference between the kinetic energy Ekin and the virial gives the pressure tensor P as 2 P ¼ ðEkin  ΞÞ V

(2.3.17)

where V represents volume of the computational box. The scalar pressure P, which can be utilized for pressure coupling for isotropic systems, is given by P ¼ traceðPÞ=3

(2.3.18)

The virial Ξ is defined as Ξ¼

1X rij Fij 2 i
(2.3.19)

86

Chapter 2

2.3.7 The Leap-Frog Integrator The leapfrog algorithm is the default MD integrator used in GROMACS, utilized for the integration of the equations of motion. However, when accurate integration with temperature and/or pressure coupling is required, the velocity Verlet integrators are preferable. Using the forces F(t) determined by the positions at time t, the leapfrog integrator utilizes positions r at time t and velocity at time t  12 Δt, to update the positions and velocities, using the following relations:     1 1 Δt (2.3.20) v t + Δt ¼ v t  Δt + FðtÞ 2 2 m   1 rðt + ΔtÞ ¼ rðtÞ + Δtv t + Δt (2.3.21) 2 It produces trajectories, which are identical to the Verlet algorithm, and the position update relation is given by rðt + ΔtÞ ¼ 2rðtÞ  rðt  ΔtÞ +

  1 FðtÞΔt2 + O Δt4 m

(2.3.22)

2.3.8 The Velocity Verlet Integrator To integrate the positions r and velocity v at time t for equations of motion, GROMACS also employs the velocity Verlet algorithm, given by   1 Δt v t + Δt ¼ vðtÞ + FðtÞ (2.3.23) 2 2m   1 rðt + ΔtÞ ¼ rðtÞ + Δtv t + Δt (2.3.24) 2   1 Δt vðt + ΔtÞ ¼ v t + Δt + Fðt + ΔtÞ (2.3.25) 2 2m

or, equivalently, rðt + ΔtÞ ¼ rðtÞ + Δtv + vðt + ΔtÞ ¼ vðtÞ +

Δt2 FðtÞ 2m

Δt ½FðtÞ + Fðt + ΔtÞ 2m

(2.3.26) (2.3.27)

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 87 When no pressure and temperature coupling is involved, the Verlet and leapfrog will produce the same trajectory. However, with a single starting file along with the same points x(0) and v(0), the velocity Verlet and leapfrog will not give identical trajectories, as a consequence of the fact that the Verlet will solve the equations corresponding to t ¼ 0 and the leapfrog at t ¼  12 Δt.

2.3.9 Reversible Integrators: The Trotter Decomposition To establish the relation between velocity Verlet and leapfrog integration, the reversible Trotter formulation of dynamics is useful. The evolution operator can be applied for a system of coupled, first-order differential equations, given by ΓðtÞ ¼ exp ðiLtÞΓð0Þ

(2.3.28)

iL ¼ Γ_ rΓ

(2.3.29)

where L is the Liouville operator and Γ is the multidimensional vector of independent variables (position and velocities): ΓðtÞ ¼

p Y

exp ðiLΔtÞΓð0Þ

(2.3.30)

i¼1

For a given time step size, the half-step-averaged kinetic energy and temperature are slightly more accurate. The kinetic energy obtained using the half-step-averaged kinetic energies will be closer to the kinetic energy obtained in the limit of small step size than the kinetic energy acquired with the full-step kinetic energy. Employing the NVE simulations, the difference is not a significant figure, since the positions and velocities of the particles remain same; however, the difference lies in the approach of interpreting the temperature. Furthermore, the trajectories produced by the integrators remain identical. An additional advantage of producing more accurate kinetic energy by half-step-averaged method is that the obtained kinetic energy will change less as the time step gets large. Furthermore, for NVT simulations, there will be a difference, due to the fact that the velocities of the particles are modified based on the distribution to be attained, corresponding to the temperature. Under such circumstances, the three methods will not produce same results. Pressure control methods can be accurately controlled by the velocity Verlet integrator as both the velocity and position can be defined at the same time t. Specifically, for large systems where the communication speed is important for parallelization, NVT ensembles are preferable. However, when fine thermodynamic details are important, velocity Verlet is desirable for calculation purpose.

88

Chapter 2

2.3.10 Temperature Coupling When implementing the concepts of molecular dynamics using NVE, that is, constant number, constant volume, and constant energy ensemble, most quantities that we intend to calculate are fundamentally obtained using constant-temperature (NVT) ensemble, which is also known as the canonical ensemble. Various approaches are available to simulate constant-temperature phenomenon, like weak coupling scheme of Berendsen, velocity-rescaling scheme, or stochastic randomization through the Andersen thermostat, with each methodology holding respective advantages. However, with the formulation of phenomenon such as heating due to viscous or frictional forces, the processes to maintain a constant temperature are not relevant, from thermodynamic point of view. Such methodology will implicitly incorporate the assumption of constant temperature, and any temperature deviation will be considered negligible.

2.3.10.1 Berendsen Temperature Coupling The Berendsen coupling algorithm implements weak coupling using first-order kinetics to an external heat bath with the given temperature T0. Using this algorithm, the system temperature is slowly corrected with the implementation of the relation: dT T0  T ¼ dt τ

(2.3.31)

which implies that the deviation in temperature decays exponentially with a time constant τ. The advantage of implementing this scheme is that the strength of the coupling can be based on user intervention. For example, the coupling can be relatively short, for example, 0.01 ps, and for other equilibrium simulation, it can be extended till 0.5 ps, which does not affect the conservation dynamics. Furthermore, the Berendsen thermostat suppresses the variations in the kinetic energy, which results in the improper generation of canonical ensemble and incorrect sampling. Moreover, the fluctuation properties, such as heat capacity, will be affected. A related thermostat that produces a correct ensemble is the velocity-rescaling thermostat, described as 8 931=2 2 > > > > =7 6 nTC Δt < T 0 6  1 7 (2.3.32) λ ¼ 41 + 5 > 1 τT > > > :T t  Δt ; 2 τ ¼ 2CV τT =Ndf k

(2.3.33)

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 89

2.3.10.2 Velocity-Rescaling Temperature Coupling A velocity-rescaling thermostat is essentially a Berendsen thermostat, with the add-on stochastic term, which helps in predicting the correct canonical ensemble, given by sffiffiffiffiffiffiffiffiffi dt KK 0 dW (2.3.34) dK ¼ ðK0  K Þ + 2 pffiffiffiffiffi Nf τ T τT Moreover, this thermostat possesses the advantage of Berendsen thermostat, which is firstorder decay, with no oscillations. Furthermore, when the NVT assemble is employed, the energy quantity that is conserved is written to the energy and log file.

2.3.10.3 Andersen Thermostat Various methods are available that can be utilized for maintaining a constant temperature. One scheme is to take an NVE integrator and then reselecting the particles from a MaxwellBoltzmann distribution periodically. This is achieved by randomizing all the velocity simultaneously, at every τ/Δt steps, or by initializing every particle with small probability at every time step of Δt/τ, where Δt represents the time step. Due to the predefined condition on the way constraints operate, simultaneous randomization should be done of all the particles. However, due to parallelization issues, the Andersen version is not currently utilized. However, the Andersen-massive can be implemented despite of the involved constraint. Furthermore, this thermostat can be utilized only with velocity Verlet algorithms, as at each time step, it operates directly on velocities. Additionally, this algorithm avoids some ergodicity issues of other thermosetting algorithms. It is a consequence of the fact that between energetically decoupled components of a system, energy cannot oscillate back and forth, as in velocity scaling motion. Moreover, the kinetics of system slows down by randomizing correlated motions of the system.

2.3.10.4 Nose-Hoover Temperature Coupling For maintaining a targeted temperature, Berendsen weak coupling is exceptionally efficient for the relaxation of a system. However, once the equilibrium is reached, it is also important to explore correct canonical ensemble. To enable canonical ensemble at such equilibrium condition, GROMACS further supports the extended-ensemble approach, proposed by Nos
e [48], which was modified later by Hoover [49]. The Hamiltonian system is continued by the introduction of thermal reservoir and a friction term, in the equations of motion.

90

Chapter 2

The friction factor is calculated as the product friction parameter and velocity of each particle. The friction parameter, or “heat bath” variable, is defined as a fully dynamic quantity, along with its momentum and equation of motion. The difference between the current kinetic energy and reference temperature yields the time derivative. In this context, the equations of motion of particle are replaced by d 2 ri Fi pξ dr i ¼  dt2 mi Q dt

(2.3.35)

where the equation of motion of friction factor parameter is given by dpξ ¼ ðT  T0 Þ dt

(2.3.36)

where T0 denotes the reference temperature and T is the current instantaneous temperature of the system. The constant Q is usually termed as “mass parameter of the reservoir.” The conserved quantity for the Nose-Hoover equations of motion is defined as H¼

N X p2ξ pi + U ðri , r2 , …, rN Þ + + Nf kTξ 2mi 2Q i¼1

(2.3.37)

where Nf is the total number of degrees of freedom. However, due to the dependence of mass parameter on reference temperature, it is some undesirable for describing the coupling strength. For maintaining the coupling strength, Q must be changed with reference temperature. Hence, it is preferable to work with the period τT of the oscillations of kinetic energy between the system and the reservoir. Its mathematical formulation is given by Eq. (2.3.38): Q¼

τ2T T0 4π 2

(2.3.38)

A strongly damped exponential relaxation can be achieved using weak coupling scheme, whereas an oscillatory relaxation is achieved using Nose-Hoover approach. The realistic time in the relaxation using Nose-Hoover coupling is several times larger than the period of the oscillations selected. Such oscillations also indicate that the time constant should be 4–5 times larger compared with relaxation time with weak coupling. Nose-Hoover dynamics consists of a system with collections of harmonic oscillators, where only a subsection of phase space is ever sampled, even if infinite long simulation process exists. Under such circumstances, the Nose-Hoover chain approach was developed, with each having its own Nose-Hoover thermostat. In the present scenario, the default number of chains is 10,

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 91 which can be modified by user. To include a chain of thermosetting particles, the existing equations can be modified as follows: d 2 ri Fi pξ1 dr i ¼  dt2 mi Q1 dt dpξ1 pξ ¼ ðT  T0 Þ  pξ1 2 dt Q2 ! p2ξi1 dpξi¼2…N pξ  kT  pξi i + 1 ¼ dt Qi1 Qi + 1 ! p2ξN1 dpξN  kT ¼ dt QN1

(2.3.39) (2.3.40)

(2.3.41)

(2.3.42)

The quantity that is conserved for Nose-Hoover chains is H¼

N M p2 M X X X pi ξk + U ðr1 , r2 , …, rN Þ + + N kTξ + kT ξk f 1 2mi 2Q0k i¼1 k¼1 k¼2

(2.3.43)

In the output the velocity data of Nose-Hoover thermostat are not included, due to the considerable consumption of space. In velocity Verlet and leapfrog dynamics, calculations for matching the temperature with the reference temperature are different. Examining the Trotter decomposition for better understanding the difference between constanttemperature integrators, consider the case of Nose-Hoover dynamics, where the splitting of Liouville operator as iL ¼ iL1 + iL2 + iLNHC

(2.3.44)

where iL1 ¼

N  X pi i¼1

iL2 ¼

N X i¼1

iLNHC ¼

mi Fi :



∂ ∂ri

∂ ∂pi

N X pξ pξ ∂ ∂  vi rvi + + ðT  T0 Þ Q Q ∂ξ ∂p ξ i¼1

(2.3.45)

(2.3.46)

(2.3.47)

When using standard velocity Verlet with Nose-Hoover temperature control the expression becomes   exp ðiLΔtÞ ¼ exp ðiLNHC Þ exp ðiL2 Δt=2Þ exp ðiL1 ΔtÞexp ðiLNHC Δt=2Þ + O Δt3 (2.3.48)

92

Chapter 2

With half-step-averaged temperature control using md-vv-avek, this decomposition will not work. This is attributed to the fact that until the second velocity step, we do not have the fullstep temperature. An alternate decomposition can be constructed, which will be reversible, by changing the place of N-heterocyclic carbenes (NHC) and velocity portions of the decomposition:   exp ðiLΔtÞ ¼ exp ðiL2 Δt=2Þ exp ðiLNHC Δt=2Þ exp ðiL1 ΔtÞexp ðiLNHC Δt=2Þ exp ðiL2 Δt=2Þ + O Δt3 (2.3.49) This formulation eases to visualize the difference between the difference genres of velocity Verlet integrator.

2.3.10.5 Group Temperature Coupling Using GROMACS, a temperature coupling can be achieved on groups of atoms, which are typically a protein or solvent. Due to imperfect energy exchange at molecular level due to effects including cutoffs, such algorithms were introduced. In this scheme the whole system is coupled to one heat bath. An additional feature of this methodology is the assignment of temperature partially to the system. In other words a part of the simulation cell can be temperature coupled while other parts not. It is achieved by specifying 1 for the time constant τT for the group, which is not to be maintained to a reference temperature. When a part of the system is assigned to a temperature, the system will still be converging to an NVT system. To minimize the errors in the temperature, arose due to discretized time step, only those molecules should be thermostatic and not the remaining molecules.

2.3.11 Pressure Coupling As the simulation cell can be coupled to a “temperature bath,” similarly, a “pressure bath” can also be assigned to it. GROMACS supports both the extended-ensemble Parrinello-Rahman approach [50,51] and the Berendsen algorithm [52] and, for velocity Verlet, the MartynaTuckerman-Tobias-Klein (MTTK) [53]. Moreover, any coupling method can be utilized for combining Parrinello-Rahman and Berendsen. However, MTTK can only be applied with Nose-Hoover temperature control.

2.3.11.1 Berendsen Pressure Coupling The box vectors and the coordinates are rescaled by the Berendsen algorithm at every step or nPC steps, with a matrix μ, which consists of the influence of first-order kinetic relaxation of the pressure toward a given pressure P0, according to

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 93 dP P0  P ¼ dt τp

(2.3.50)

nPC Δt βij P0ij  Pij ðtÞ 3τp

(2.3.51)

The scaling matrix μ is defined as μij ¼ δij 

where β is the isothermal compressibility of the system. In most likely cases, it would correspond to be a diagonal matrix. When the system is completely anisotropic, it has to be rotated. Such rotation, in the first order in the scaling, is approximated, which is typically less than 10-4. However, the actual scaling matrix μ’ is given by 0 1 μxx μxy + μyx μxz + μzx (2.3.52) μyy μyz + μzy A μ0 ¼ @ 0 0 0 μzz Here the velocities are neither rotated nor scaled. The Berendsen scaling can be implemented isotropically, which implies that a diagonal matrix with elements of size trace (P)/3 is used. Semiisotropic scaling can be used for systems with interfaces. In this particular case, isotropically scaling is achieved for x/y-directions and independent scaling for z-direction. Furthermore, compressibility in the x/y- or z-directions can be initialized to zero, for scaling in the remaining direction(s). However, if full anisotropic deformations are allowed, along with the usage of constraints, then either scaling must be allowed slowly or decrease the time step for avoiding errors from constraint algorithms. Moreover the Berendsen pressure control algorithm does not yield the exact NPT ensemble, even though the algorithm produces correct average pressure.

2.3.11.2 Parrinello-Rahman Pressure Coupling When the cases deal with fluctuation in pressure or volume, for example, calculating the thermodynamic properties preferably for small systems, it poses a problem that the exact ensemble may not be properly defined, and the true NPT ensemble is not simulated. In this context a constant-pressure simulation is also supported by GROMACS, using the ParrinelloRahman approach, which has similarity with Nose-Hoover coupling and theoretically provides true NPT ensemble. Implementing the Parrinello-Rahman barostat the box vectors represented by matrix b obey the matrix equation of motion: db2 1 ¼ VW 1 b0 ðP  Pref Þ 2 dt

(2.3.53)

94

Chapter 2

where V denotes the volume of the box and W represents the matrix parameter, which determines the strength of the coupling. The matrices P and Pref are the current and reference pressures, respectively. The governing equations of motion are also changed. The Parrinello-Rahman modification is d 2 r Fi dr i ¼ M 2 mi dt dt

(2.3.54)

 db0 db 0 0 1 + b b M¼b b dt dt

(2.3.55)

1

The inverse parameter matrix W1 specifies the strength of the coupling and the nature of box deformation. When the elements of W1 are zero, the box restriction will be automatically fulfilled. The coupling strength can be automatically calculated in GROMACS, and the approximate isothermal compressibility β and the pressure time constant τp must be provided in the input file:  1  4π 2 βij W ij ¼ 2 3τp L

(2.3.56)

If the pressure is very far from equilibrium, very large box oscillations may be resulted from Parrinello-Rahman coupling, which could even crash the run. In such situations, either the time constant have to be increased or the weak coupling scheme may be utilized to reach the target pressure, and then the Parrinello-Rahman coupling can be implemented, once the system is in equilibrium. In addition, using the leapfrog algorithm, the pressure at time t is not available, until the time step has completed. Hence the pressure from the previous step must be utilized, which makes the algorithm not directly reversible. Such situations may be inappropriate for high precision thermodynamic calculations.

2.3.11.3 Surface-Tension Coupling Considering a system that consists of more than one distinguishable phase, which are separated by surfaces, parallel to the x-y plane, then the surface tension and the z-component of the pressure cab are coupled to a pressure bath. Till the present version of GROMACS, this scheme only works with Berendsen pressure coupling. The average surface tension γ(t) can be calculated from the difference between the normal and lateral pressure: 1 γ ðtÞ ¼ n

L ðz 

 Pxx ðz, tÞ + Pyy ðz, tÞ Pzz ðz, tÞ  dz 2

0

(2.3.57)

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 95   Pxx ðtÞ + Pyy ðtÞ L ¼ Pzz ðtÞ  2 n

(2.3.58)

where Lz represents the height of the box and n is the number of surfaces. Moreover, by scaling the height of the box with μzz, the pressure in the z-direction can be corrected, given by ΔPzz ¼

Δt fP0zz  Pzz ðtÞg τp

μzz ¼ 1 + βzz ΔPzz

(2.3.59) (2.3.60)

This formulation is similar to the normal pressure coupling, expect the omission of 1/3 factor. The pressure correction in the z-direction is then utilized to get the correct convergence for the surface tension value γ 0. The correction factor for the box length, in the x/y-direction, is given by    Pxx ðtÞ + Pyy ðtÞ Δt nγ 0 (2.3.61) β  Pzz ðtÞ + ΔPzz  μx=y ¼ 1 + 2 2τp x=y μzz Lz In most cases, the τp values will scale, under incorrect compressibility. The convergence of the surface tension is affected by compressibility, when surface tension coupling is applied. Moreover, under constant box length (βzz ¼ 0), the ΔPzz is also set to zero, which is necessary for obtaining correct surface tension.

2.3.12 The Complete Update Algorithm For velocities and coordinate update, the complete algorithm is given using leapfrog. GROMACS provides the provision of preventing the motion of selected particles, which are defined as a “freeze group.” Under such situations, a freeze factor fg is a vector quantity. Hence, this freeze factor and the acceleration ah must be taken into account for updating the algorithm of the velocities, which becomes

     Δt Δt FðtÞ (2.3.62) ¼ fg λ v t + Δt + ah Δt v t+ 2 2 m

2.3.13 Output Step The most important step of the MD run is the creation of trajectory file, which contains the coordinates of the particles and (optionally) the velocities at fixed intervals. The trajectory file contains the following information: positions, velocities, forces, dimensions of the simulation cell, integration step, integration time, etc. When velocity Verlet integrators are used, velocity

96

Chapter 2

labeled at time t is for that time, whereas for other integrators (e.g., stochastic dynamics), velocities at time t are for time t  1/2 Δt. Moreover, the trajectory files are lengthy in nature; therefore it is preferable to omit steps. To retain information, it is adequate to write a frame every 15 steps, since for the highest frequency in the system, 30 steps are made per period.

2.3.14 Advantage and Functional Characteristics The following advantages are offered by the GROMACS software: i. GROMACS provides high performance for computational purpose. ii. GROMACS is user-friendly, along with the clear text formatting of topologies and parameter files. A number of consistency checking are involved, and unambiguous error messages are issued, when something is wrong. iii. With no scripting language involved, all programs use an interface, along with command line options for input and output files. iv. Elapsed simulation time will be continuously available in the interface, along with the expected time for the final completion and logging of data. v. GROMACS is a free software, available under Lesser General Public License (LGPL), version 2.1. The GROMACS software consists of a preprocessor, a parallel MD, and an energy minimization program, which can utilize an arbitrary number of processors, an optional monitor, and several analysis tools. The programming for the language is written in ANSI C. The following are the functional characteristics of the software: i. Nature of physical problem Analysis and prediction of dynamic behavior of macromolecules can be made under external driving forces. ii. Method of solution GROMACS uses classical Newtonian equations of motion and force fields, based on variable or fixed bonded or nonbonded molecular interaction. The developed system of molecules is coupled to external bath of constant temperature and pressure. Periodic conditions (rectangular) can be applied to the simulation, which would provide an ease to establishing repetitive behavior of molecular structure. iii. Restriction on the complexity of the problem The computational file size, limited by the memory and the number of processors, depends on the complexities of the molecular interactions.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 97 iv. Running time A typical small biomolecule (a peptide of 20 residues) in water (800 atoms) runs in 100-time steps (0.2 ps) in 1 min on a 32-i860 processor system. In the past years, GROMACS has evolved and responded to rapid growth, to be a large project, with almost two million lines of code. From the first release of GROMACS the most significant change in the software is that it has moved to C++. While many parts of the code are still coded using C algorithm, it will take a couple of years to have a complete transition to C++. Incorporating the transition to a more versatile language has led to the involvement of modularity in the code. Partitioning of functionality and consequently in the memory has led to better handling of memory and errors. Within the simulation, parallel computer functionality can be enabled, which would lead to the splitting of work, into independent entity of modules, for faster transaction toward less computational time lapse.

2.3.15 Application of GROMACS GROMACS has emerged to be an extremely useful computational code for computationally replicating the system of molecular cluster and aiding in fetching technical data, which is beneficial for gaining deeper understanding at the atomic scale. Working at the molecular scale aids in optimizing of devices, for creating state-of-the-art equipment. From the modeling of medical devices and implants to simulation of nanoscale flows, GROMACS is in continuous process, toward our understanding at small scale, and the verification of physical principles, which are visually prominent at macroscale. In this section, a list on the applications of GROMACS is presented; however, this is a small note on the applications, and the readers are suggested not to be confined with the discussed utility of the software.

2.3.15.1 Biochip Devices GROMACS provides molecular dynamics simulations for the construction and optimization of biochip devices, with miniature passages as shown in Fig. 2.3.7. With the molecular dynamics approach, the flow problem can be accurately measured and resolved. The thermal flow problems is inherent to the gradients of both velocity and thermal profile; however, interlink between both the gradient can be understood on the basis of molecular approach. Physical conditions of frostbite or burns suffered by the patients can also be studied by the simulation approach. The simulations mainly focus on the effect of effective microchannel width and the thermal boundary conditions at the solid-fluid interface on the velocity profile. GROMACS provides useful insights on the influence of global and local effects on bioflow models, important for configuration and the operation of biochip devices.

98

Chapter 2 Temperature boundary (T) Wall atom

Direction of flow

Effective channel width (D) Molecule of fluid

y

x

Fig. 2.3.7 Schematic diagram of a simulation geometry.

2.3.15.2 Molecular Modeling of Biomolecules Calculating the molecular dynamics data and solving the energy equation are an essential feature of GROMACS. Modeling of complex biomolecule like protein, which is the most important nutrient present in the food apart from carbohydrates and fats, can be simulated with molecular dynamics (MD). The functional properties exhibited by the molecules vary at different processing and storage stages. During processing of various molecules, under various external thermal fields, like boiling and roasting. or nonthermal fields, like pulsed electric fields and pressure processing, external stresses are applied, which would affect the shelf life of the product. Therefore efforts are made by the researchers in understanding the complex dynamics of processing protein structures in recent years. Usage of magnetic resonance imaging (MRI) and X-ray diffraction has been for studying the dynamics of molecular structure; however, due to various limitations depending upon the technique adopted and their intrinsic expensive nature, MD simulation techniques are emerging as a viable alternative to overcome the standard issues existing due to the prevailing techniques.

References [1] S. Kemmerer, J.C. Voss, R. Faller, Molecular dynamics simulation of dipalmitoylphosphatidylcholine modified with a MTSL nitroxide spin label in a lipid membrane, Biochim. Biophys. Acta Biomembr. 1828 (2013) 2770–2777. [2] M. Yu, Computational Modeling of Protein Dynamics With GROMACS and Java, Master’s Projects, 2012, p. 267. http://scholarworks.sjsu.edu/etd_projects/267. [3] M.J. Abraham, T. Murtola, R. Schulz, et al., GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX 1–2 (2015) 19–25.

Overview of BIOVIA Materials Studio, LAMMPS, and GROMACS 99 [4] H.J.C. Berendsen, D. Van Der Spoel, R. Van Drunen, GROMACS: a message-passing parallel molecular dynamics implementation, Comput. Phys. Commun. 91 (1995) 43–56. [5] Q.W. Kulig, M. Pasenkiewicz-gierula, T. Ro´g, Topologies, structures and parameter files for lipid simulations in GROMACS with the OPLS-aa force field: DPPC, POPC, DOPC, PEPC, and cholesterol, Data Brief 5 (2015) 1–4. [6] N. Singh, V. Dalal, P. Kumar, Structure based mimicking of Phthalic acid esters (PAEs) and inhibition of hACMSD, an important enzyme of the tryptophan kynurenine metabolism pathway, Int. J. Biol. Macromol. 108 (2018) 214–224. [7] E. Somogyi, A. Abi, P.J. Ortoleva, ProtoMD: a prototyping toolkit for multiscale molecular dynamics, Comput. Phys. Commun. 202 (2016) 337–350. [8] D.H. De Jong, S. Baoukina, H.I. Ingo´lfsson, S.J. Marrink, Martini straight: boosting performance using a shorter cutoff and GPUs, Comput. Phys. Commun. 199 (2016) 1–7. [9] N.P. Walter, A. Jaiswal, Z. Cai, Y. Zhang, LiquidLib: a comprehensive toolbox for analyzing classical and ab initio molecular dynamics simulations of liquids and liquid-like matter with applications to neutron scattering experiments, Comput. Phys. Commun. 228 (2018) 209–218. [10] C. Kutzner, D.A. K€opfer, J. Machtens, B.L. De Groot, C. Song, U. Zachariae, Insights into the function of ion channels by computational electrophysiology simulations, Biochim. Biophys. Acta 1858 (2016) 1741–1752. [11] M. Ferna´ndez-penda´s, E. Akhmatskaya, J.M. Sanz-serna, Adaptive multi-stage integrators for optimal energy conservation in molecular simulations, J. Comput. Phys. 327 (2016) 434–449. [12] J. Vidya, S. Sajitha, M.V. Ushasree, et al., Genetic and metabolic engineering approaches for the production and delivery of L. asparaginases: an overview, Bioresour. Technol. 245 (2017) 1775–1781. [13] C. Blau, H. Grubmuller H., g_contacts: fast contact search in bio-molecular ensemble data, Comput. Phys. Commun. 184 (2013) 2856–2859. [14] M.T. Degiacomi, V. Erastova, M.R. Wilson, Easy creation of polymeric systems for molecular dynamics with Assemble! Comput. Phys. Commun. 202 (2016) 304–309. [15] T. Ro´g, A. Or, A. Llorente, Data including GROMACS input files for atomistic molecular dynamics simulations of mixed, asymmetric bilayers including molecular topologies, equilibrated structures and force field for lipids compatible with OPLS-AA parameters, Data Brief 7 (2016) 1171–1174. [16] B. Harish, S. Kranthi, A. Singh, V. Raghavan, Effects of thermal and electric fields on soybean trypsin inhibitor protein: a molecular modeling study, Innov. Food Sci. Emerg. Technol. 35 (2016) 9–20. [17] D.T.W. Lin, A study of local effect and global effect on the microthermal bio-flows by molecular dynamics, Int. J. Biol. Macromol. 41 (2007) 260–265. [18] C. Kutzner, R.T. Ullmann, B.L. De Groot, U. Zachariae, H. Grubmueller, Ions in action-studying ion channels by computational electrophysiology in GROMACS, Biophys. J. 112 (2017) 139a. [19] B. Loubet, W. Kopec, H. Khandelia, Application of the virtual site technique to lipids in gromacs, hydrogens degrees of freedom removal and performance increase, Biophys. J. 106 (2014) 17a. [20] C. Kutzner, J. Czub, H. Grubmueller, Keep it flexible: driving macromolecular rotary motions in atomistic simulations with gromacs, Biophys. J. 102 (2012) 171a. [21] M. Hassan, I. Saika-voivod, V. Booth, All-atom molecular dynamics simulations of lung surfactant protein B: structural features of SP-B promote lipid reorganization, Biochim. Biophys. Acta Biomembr. 1858 (2016) 3082–3092. [22] F. Jensen, Introduction to Computational Chemistry, thir ed., John Wiley & Sons Ltd, West Sussex, UK, 2015. [23] T. Hansson, C. Oostenbrink, W.F. Van Gunsteren, Molecular dynamics simulations, Curr. Opin. Struct. Biol. 12 (2002) 190–196. [24] C.B. Roth, Polymer Glasses, CRC Press, Taylor & Francis Group, Boca Raton, FL, 2017. [25] A.L. Galbraith, C.K. Hall, Vapor-liquid phase equilibria for mixtures containing diatomic Lennard-Jones molecules, Fluid Phase Equilib. 241 (2006) 175–185. [26] J.G. Kirkwood, E.M. Boggs, The radial distribution function in liquids, J. Chem. Phys. 10 (1942) 394–402. [27] H. Kuhn, Viscosity, sedimentation, and diffusion of long-chain molecules in solution as determined by experiments on large-scale models, J. Colloid Sci. 5 (1950) 331–348. [28] M. Matsumoto, Molecular dynamics of fluid phase change, Fluid Phase Equilib. 144 (1998) 307–314. [29] J. Nord, Molecular dynamics study of defect formation in GaN cascades, Nucl. Instrum. Methods Phys. Res. Sect. B 202 (2003) 93–99.

100 Chapter 2 [30] F.J. Lin, J.T. Muckerman, Solution of the time-dependent Schrodinger equation employing a basis of explicit discrete-coordinate eigen functions: spherical and azimuthal symmetry, adiabaticity, and multiphoton excitation of a rotating Morse oscillator, Comput. Phys. Commun. 63 (1991) 538–568. [31] A. Farhi, B. Singh, Calculation of molecular free energies in classical potentials, New J. Phys. 18 (2016). 023039 (1-13). [32] C.J. van Oss, J. Walker, Concentration dependence of the binding constant of antibodies, Mol. Immunol. 24 (1987) 715–717. [33] S. Kossida, D. Vlachakis, E. Bencurova, N. Papangelopoulos, Current state-of-the-art molecular dynamics methods and applications, Adv. Protein Chem. Struct. Biol. 94 (2014) 269–313. [34] R.B. Pearson, J.L. Richardson, D. Toussain, A fast processor for Monte-Carlo simulation, J. Comput. Phys. 51 (1983) 241–249. [35] M. Neumanna, M. Johnsonb, L. Von Lauea, H.P. Trommsdorff, Proton tunneling in molecular solids, J. Lumin. 67 (1995) 146–151. [36] E.F. Hammel, The low temperature properties of helium three, Progr. Low Temp. Phys. 1 (1955) 78–107 (Chapter V). [37] R.G. Woolley, Molecular structure and the Born-Oppenheimer approximation, Chem. Phys. Lett. 45 (1977) 393–398. [38] L.E. Bilston, K. Mylvaganam, Molecular simulations of the large conductance mechanosensitive (MscL) channel under mechanical loading, FEBS Lett. 512 (2002) 185–190. [39] M. Levitt, Protein folding by restrained energy minimization and molecular dynamics, J. Mol. Biol. 170 (1983) 723–764. [40] R. Lindh, A. Bernhardsson, G. Karlstr€om, P.A. Malmqvist, On the use of a Hessian model function in molecular geometry optimizations, Chem. Phys. Lett. 241 (1995) 423–428. [41] M.B. Kubitzki, B.L. De Groot, Molecular dynamics simulations using temperature-enhanced essential dynamics replica exchange, Biophys. J. 92 (2007) 4262–4270. [42] J.Y. Yi, J. Bernholc, P. Salamon, Simulated annealing strategies for molecular dynamics, Comput. Phys. Commun. 66 (1991) 177–180. [43] H. Wan, F.Y. Hunt, A.J. Kearsley, An optimization approach to multiple sequence alignment, Appl. Math. Lett. 16 (2003) 785–790. [44] S.J. Irudayam, M.L. Berkowitz, Binding and reorientation of melittin in a POPC bilayer: computer simulations, Biochim. Biophys. Acta Biomembr. 1818 (2012) 2975–2981. [45] A.T. Frank, Q. Zhang, H.M. Al-hashimi, I. Andricioaei, Slowdown of interhelical motions induces a glass transition in RNA, Biophys. J. 108 (2015) 2876–2885. [46] L. Sun, J.K. Noel, J.I. Sulkowska, H. Levine, J.N. Onuchic, Connecting thermal and mechanical protein (un) folding landscapes, Biophys. J. 107 (2014) 2950–2961. [47] S.M. Cory, Y. Liu, M.I. Glavinovic, Interfacial interactions of glutamate, water and ions with carbon nanopore evaluated by molecular dynamics simulations, Biochim. Biophys. Acta Biomembr. 1768 (2007) 2319–2341. [48] S. Nos
e, A molecular dynamics method for simulations in the canonical ensemble, Mol. Phys. 52 (1984) 255–268. [49] W.G. Hoover, Canonical dynamics: equilibrium phase-space distributions, Phys. Rev. A 31 (1985) 1695–1697. [50] M. Parrinello, A. Rahman, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys. 52 (1981) 7182–7190. [51] S. Nos
e, M.L. Klein, Constant pressure molecular dynamics for molecular systems, Mol. Phys. 50 (1983) 1055–1076. [52] H.J.C. Berendsen, J.P.M. Postma, W.F. Van Gunsteren, A. Dinola, J.R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys. 81 (1984) 3684–3690. [53] G.J. Martyna, M.E. Tuckerman, D.J. Tobias, M.L. Klein, Explicit reversible integrators for extended systems dynamics, Mol. Phys. 87 (1996) 1117–1157.