Generalized Potential Energy Finite Elements for Modeling Molecular

Article pubs.acs.org/jcim

Generalized Potential Energy Finite Elements for Modeling Molecular Nanostructures Stavros Chatzieleftheriou,† Matthew R. Adendorff,‡ and Nikos D. Lagaros*,† †

Institute of Structural Analysis & Antiseismic Research, Department of Structural Engineering, School of Civil Engineering, National Technical University of Athens, 9 Heroon Polytechniou Street, Zografou Campus, GR-15780 Athens, Greece ‡ Laboratory of Computational Biology & Biophysics, Department of Biological Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139-4307, United States S Supporting Information *

ABSTRACT: The potential energy of molecules and nanostructures is commonly calculated in the molecular mechanics formalism by superimposing bonded and nonbonded atomic energy terms, i.e. bonds between two atoms, bond angles involving three atoms, dihedral angles involving four atoms, nonbonded terms expressing the Coulomb and Lennard-Jones interactions, etc. In this work a new, generalized numerical simulation is presented for studying the mechanical behavior of three-dimensional nanostructures at the atomic scale. The energy gradient and Hessian matrix of such assemblies are usually computed numerically; a potential energy finite element model is proposed herein where these two components are expressed analytically. In particular, generalized finite elements are developed that express the interactions among atoms in a manner equivalent to that invoked in simulations performed based on the molecular dynamics method. Thus, the global tangent stiffness matrix for any nanostructure is formed as an assembly of the generalized finite elements and is directly equivalent to the Hessian matrix of the potential energy. The advantages of the proposed model are identified in terms of both accuracy and computational efficiency. In the case of popular force fields (e.g., CHARMM), the computation of the Hessian matrix by implementing the proposed method is of the same order as that of the gradient. This analysis can be used to minimize the potential energy of molecular systems under nodal loads in order to derive constitutive laws for molecular systems where the entropy and solvent effects are neglected and can be approximated as solids, such as double stranded DNA nanostructures. In this context, the sequence dependent stretch modulus for some typical base pairs step is calculated.

1. INTRODUCTION During the last 20 years, significant progress in the field of nanotechnology has been established. Despite the fact that nanotchnology covers a very broad range of scientific topics and applications, the principal characteristic of extremely small dimension(s) (1−100 nm) unifies this diverse field. The applications are numerous, steaming from various fields, such as surface science, organic chemistry, molecular biology, semiconductor physics, microfabrication, and many others. This inherent diversity and the goal to manipulate objects on the nanoscale often demands a multidisciplinary approach dissolving the boundaries between the various fields of traditional science. In the field of structural engineering, as nanostructures find numerous applications in forming new categories of composites (i.e., nanocomposites), methods developed by chemists, chemical engineers, and biophysicists to study chemical and biological systems (i.e., molecular mechanics and dynamics) have been incorporated in the traditional engineering analyses leading to new hybrid methods, such as the multiscale FEM (finite element method), MSM (molecular structural mechanics), and others. One important milestone for this development © XXXX American Chemical Society

was the discovery of CNT (carbon nanotubes), a famous carbon allotrope cylindrical nanostructure with extraordinary properties which finds applications in nanotechnology, electronics, optics, materials science, and others. Efforts to simulate the material properties at the nanoscale stem from various fields,1−11 such as molecular dynamics, molecular mechanics, molecular structural mechanics equivalent continuum formulations, and others. In biological space, nucleic acid nanotechnology has developed into a mature synthetic technique for the design of hierarchical molecular assemblies,12−18 and these constructs have found application to a number of contemporary research endeavors, including, but not limited to, biomimicry,19 macromolecular sensing,20 energy transfer/photonics,21,22 drug delivery,23,24 and actuated nanoreactor design.25 Continuum model representations of the nucleic acid building blocks that comprise these assemblies, based on Kirchoff’s theory of elastic rods, have been successfully applied to studying DNA molecules for many years.26−34 For large DNA Received: June 16, 2016

A

DOI: 10.1021/acs.jcim.6b00356 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling nanostructures, the finite element-based35 computational modeling framework CanDo has been the most successful.36−38 In its original implementation36 DNA origami nanostructures were modeled as bundles of isotropic elastic rods that were rigidly constrained to their nearest neighbors at specific crossover positions. Later modifications included the introduction of nicked-duplex and single-stranded DNA modeling,38 and extension to lattice-free designs via parametrization of junction alignment elements.37 The introduction of alignment elements required parametrization of the mechanical junction models by all-atom simulations, an example of the possible synergy between the finite element and all-atom molecular dynamics methods. Very recently, a Brownian Dynamics (BD) model has been added to the CanDo framework to simulate the equilibrium dynamics of large nanostructures in solution,39 providing a powerful predictive tool for understanding the dynamic behavior of these large assemblies. The success of the aforementioned continuum models is strongly dependent on the assigned material properties of the coarse-grained elements used to represent a particular nanostructure, and efficient finegrained models are needed that can accurately parametrize them. In light of these new developments, nanostructure analysis and design can benefit from an intedisciplinary approach incorporating traditional structural and mechanical engineering principles for buildnig structures and machines in the nanoscale. Energy minimization is a very important procedure in the analysis of molecular systems. In order to perform molecular dynamics simulations, the system has to be at a point not far from a local minimum of the potential energy,40 and for specific applications, e.g. normal-mode analysis,41 the system has to be at a minimum. From a structural mechanics point of view, the minimization of energy suggests equilibrium between nodal loads and internal forces of the structure. A fundamental principle of structural mechanics and elasticity theory is subdomain equilibrium under general equilibrium.42 Considering equilibrium of a structure suggests equilibrium for every subpart of that structure, and even if the structure is not in equilibrium under zero loads, fictitious nodal forces can be considered, that put the structure in equilibrium. This property allows partitioning of the structure into an appropriate number of elements (one element for each energy term). The analysis of the internal forces is equivalent to the computation of the energy gradient, and the global tangent stiffness matrix assembly is equivalent to the Hessian matrix calculation. In order to verify that argument, butane and DNA molecules were analyzed, and the numerically computed Hessian was compared to the structure global stiffness matrix by the introduction of appropriate metrics and found to be equivalent. Furthermore, as was indicated from the numerical examples examined herein but also from a combinatorial analysis, the analytical protocol proposed here was several orders of magnitude faster than the conventional numerical computation. For molecular systems that can be approximated as solids revolving in a potential energy minimum, where the entropic effects can be neglected, this analysis can incorporate nodal loads and, via minimization, derive constitutive laws. In this context, the sequence dependent stretch modulus for indicating two base pair double stranded DNA molecules was derived.

2. METHODOLOGY 2.1. Generalized Formulation of Potential Energy Finite Elements. 2.1.1. Force Fields as a Reference. In the context of molecular modeling, a force field refers to the functional form and parameter sets used to calculate the potential energy of a system of atoms or coarse-grained particles in molecular mechanics and molecular dynamics simulations. Potential energy force fields provide parameters for every type of atom in a system, including hydrogen. The basic functional form of potential energy in molecular mechanics includes bonded terms for interactions of atoms that are linked by covalent bonds (intramolecular bonds, i.e. internal terms) and nonbonded (intermolecular interactions, i.e. external terms also called “noncovalent”) terms. The intramolecular portion of the potential energy function includes terms for the bonds, valence (in plane) angles, dihedral (torsion) angles, improper dihedral (out of plane bending) angles, and the Urey−Bradley 1,3-term, while the intermolecular terms include the long-range electrostatic and van der Waals forces. The general form for the total potential energy can be written as the sum of bonded and nonbonded energy terms, where the components of the covalent energy term are the bond, angle, dihedral, and Urey−Bradley contributions while those of the noncovalent energy term are the electrostatic and van der Waals ones: V = Vbonded + Vnonbonded = Vbonds + Vangles + Vdihedrals + Vimpropers + Vnonbonded

(1)

The bond and angle terms are usually modeled by quadratic energy functions that do not allow bond breaking, while the functional form for dihedral energy is highly variable. Additional, improper torsional terms can also be added to enforce planarity and “cross-terms” (Urey−Bradley 1,3-term) that describe coupling of different internal variables. Some force fields also include explicit terms for hydrogen bonds. Nonbonded energy terms are most computationally intensive, and a popular choice is to limit interactions to pairwise energies. The van der Waals (vdW) term is usually computed with a Lennard-Jones (LJ) potential and the electrostatic term with the Coulombic (CL) potential law. CHARMM (Chemistry at HARvard Macromolecular Mechanics)43,44 is a widely used set of potential energy force fields for molecular dynamics and is the one that was adopted herein in order to implement the proposed formulation of the generalized potential energy finite element. 2.1.2. Generalized Internal and Nodal Forces for Energy Terms Assuming Equilibrium. Based on the proposed generalized formulation, for each component of the energy terms (both bonded and nonbonded ones), a special finite element35 is formulated composed of nodes equal to the number of atoms (a node is located at each atom) expressing the contribution of the specific energy term. Every energy term implies a generalized internal force. At an arbitrary configuration, fictitious nodal forces can be considered that put the structure in equilibrium.42 These nodal forces satisfy equilibrium with the generalized internal force. Through the equilibrium equation implemented for the ith generalized finite element (corresponding to any energy term), the reaction forces f i are computed according to the following expression:

fi = ni ·LTi ·cTi B

(2) DOI: 10.1021/acs.jcim.6b00356 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling where f i is the vector of element forces for element i, ni is the internal generalized force (force Ni or moment Mi) for this element, LTi is the matrix of transformation to the global coordinate system for element i, and cTi is the equilibrium vector. The proposed formulation is general, since it is applicable to any force field. Bond Stretches and Nonbonded Interaction. The first term in the energy function (eq 1) accounts for the bond stretches, while the last term represents nonbonded interactions between pairs of atoms (i, j). For the case of the bonded and nonbonded components of the energy terms, as shown in Figure 1 the

2 the interaction between three atoms (represented by nodes I, J, and K) is described. While the expressions corresponding to

Figure 2. In plane angle bond interaction.

the energy function along with the internal force (Min) and stiffness (k̅θ) are given below, together with the transformation matrix L and the corresponding equilibrium vector c:

Figure 1. Bond stretches and nonbonded interaction.

interaction between two atoms (represented by nodes I and J) is described. The expressions corresponding to the energy function along with the internal force (Nin) and tangent stiffness (k̅b) of the interaction are given below, together with the transformation matrix L and the corresponding equilibrium vector c: ⎧1 2 ⎪ kb·(b‐bο) 2 ⎪ 12 ⎪ q ·q ⎡⎛ R ⎛ R min, ij ⎞6 ⎤ i j min, ij ⎞ ⎪ ⎢ + εij · ⎜ ⎟ − 2⎜ ⎟⎥ V (b) = ⎨  ⎢⎣⎝ b ⎠ 4·π ·D·b ⎝ b ⎠ ⎥⎦ ⎪  ⎪Coulombic Lennard ‐Jones ⎪ ⎪ ⎪ D ·(1 − e−a(b − bο))2 ⎩ e ∂Vb , ∂b

kb̅ =

⎡ λT 0 ⎤ 1 1×3⎥ L=⎢ , T ⎥ ⎢0 ⎣ 1 × 3 λ1 ⎦

01 × 3 ⎤ ⎥ 01 × 3 ⎥ ⎥, λ1T ⎥⎥ ⎥ λ 2T ⎦

(4)

where kθ is the angle force constant and while θ0 and θ are the reference and current in-plane angles, respectively. λ1 and λ2 define directional vectors perpendicular to the vectors along axes defined between nodes I-J and J-K (unit vectors, see Figure 2), respectively; and la, lβ are the distances between nodes J-K and I-J, respectively. Dihedral Angle Bond Interaction (Proper and Improper). The third term is used for the dihedrals (i.e., torsion angles), while the next one accounts for the improper dihedrals that represent out of plane bending. For the case of the dihedral angle bond interaction component of the energy terms, as shown in Figure 3, the interaction between four atoms (represented by nodes I, J, K, and L) is described. While the expressions corresponding to the energy function along with the internal force (Min) and stiffness (k̅φ) are given below:

(3)

Nin =

1 kθ ·(θ − θο)2 2

⎡ λT 0 1×3 ⎢ 1 ⎢ dVangle 01 × 3 λ 2T = kθ ·(θ − θο), L = ⎢ Min = ⎢0 dθ ⎢ 1 × 3 01 × 3 2 ∂ Vangle ⎢ ⎣ 01 × 3 01 × 3 kθ̅ = 2 ∂θ ⎡ 1/la ⎤ ⎢ ⎥ ⎢−1/lβ ⎥ ⎥ cT = ⎢ ⎢ −1/la ⎥ ⎢ ⎥ ⎢⎣ 1/lβ ⎥⎦ Vangle(θ ) =

∂ 2Vb ∂b2

⎡−1⎤ cT = ⎢ ⎥ ⎣1⎦

where kb is the bond force constant, De is the well depth, a controls the width of the potential, while b0 and b are the reference and current bond lengths, respectively. In expression for the energy function, the top term is for the bonded term, the middle term is for the nonbonded term, and the bottom term is the Morse potential for hydrogen bonds. By definition, the nonbonded forces are only applied to atom pairs separated by at least three bonds. The van der Waals energy is calculated with a standard 12−6 Lennard-Jones potential (VLJ) and the electrostatic energy with a Coulombic potential (VCL). In the nonbonded term of the energy function, qi and qj are the partial atomic charges of atoms i and j, respectively, εij is the well depth, Rmin,ij is the radius in the Lennard-Jones 6−12 term used to treat the van der Waals interactions, and b (also denoted with rij in the literature) is the distance between atoms i and j. λ1 is the unit vector of the axis defined between nodes I and J (see Figure 1). In Plane Angle Bond Interaction. The second term in eq 1 accounts for the bond angles. For the case of the bonded angle interaction component of the energy terms, as shown in Figure

⎧ kφ·[1 + cos(n·φ − δ)] for the proper term ⎪ Vdihedral(φ) = ⎨ ⎪ kφ·(φ − φ )2 for the improper term ⎩ 0 Min =

∂Vdihedral , ∂φ

kφ̅ =

∂ 2Vdihedral ∂φ 2 (5a)

where kφ is the dihedral force constant for the proper term while for the case of the improper one it represents the improper force constant. In the proper term, n and δ are constants (n is the multiplicity of the energy function and δ is the phase shift) while φ is the dihedral angle. In the improper term, φ0 and φ are the reference and current out of plane angles, respectively; and C


Article

Journal of Chemical Information and Modeling

Figure 3. Dihedral angle bond interaction (proper and improper).

⎡ λT ⎢ 1 ⎢0 ⎢ 1×3 ⎢0 1×3 L=⎢ ⎢ ⎢ 01 × 3 ⎢ ⎢ 01 × 3 ⎢ ⎣ 01 × 3

01 × 3 01 × 3 01 × 3 ⎤ ⎥ λ1T 01 × 3 01 × 3 ⎥ ⎥ 01 × 3 λ1T 01 × 3 ⎥⎥ , ⎥ 01 × 3 01 × 3 λ4T ⎥ ⎥ λ4T 01 × 3 01 × 3 ⎥ ⎥ 01 × 3 λ4T 01 × 3 ⎦

⎤ ⎡ 1 ⎥ ⎢ la·sin θa ⎥ ⎢ ⎢ 1 1 ⎥ − ⎥ ⎢− lb·tan θa ⎥ ⎢ la·sin θα ⎥ ⎢ 1 ⎥ ⎢ lb·tan θa ⎥ ⎢ T c =⎢ ⎥ 1 ⎥ ⎢ ⎥ ⎢ lc·sin θc ⎥ ⎢ 1 ⎥ ⎢ ⎥ ⎢ lb·tan θc ⎥ ⎢ ⎢ 1 1 ⎥ − ⎥ ⎢− lb tan θc ⎦ ⎣ lc·sin θc

order to restrain the motions of the bonds involved in the angle. This potential function is often referred to as the Urey− Bradley potential. The Urey−Bradley component (cross-term accounting for angle bending using I−K nonbonded interactions) comprises the fifth term in the energy function (eq 1). The Urey−Bradley energy term describes the interaction between two atoms (represented by nodes I and K). The expressions corresponding to the energy function along with the internal force (Nin) and stiffness (k̅UB) of the interaction are given below, together with the transformation matrix L and the corresponding equilibrium vector c: VUB(b) =

1 kUB· (b − bο)2 2 2

Nin =

∂Vb ∂ VUB , kUB ̅ = ∂b ∂b2

⎡ λT 0 ⎤ 1 1×3⎥ L=⎢ , ⎢0 T ⎥ ⎣ 1 × 3 λ1 ⎦

⎡− 1⎤ cT = ⎢ ⎥ ⎣1⎦ (6)

where kUB is the respective force constant, while b0 and b are the reference and current bond lengths between atoms I and K in the harmonic potential (of the in plane angle bond interaction configuration of the three atoms I, J, and K). λ1 is the unit vector of the axis defined between nodes I and K. Superposition of Potential Terms through Forces. With appropriate summation of the nodal forces at the corresponding degrees of freedom for each element, the energy gradient can be obtained:

(5b)

nele

are the transformation matrix L and the corresponding equilibrium vector c, where λ1 and λ4 represent directional vectors perpendicular to planes defined by nodes J-K-L and I-JK, respectively (unit vectors, see Figure 3), and la, lβ, lc, θa, θβ, and θc are the distances between the four nodes (I, J, K, and L) and the angles shaped by the axes defined between the four atoms as depicted in Figure 3. Urey−Bradley Bond Interaction. In the CHARMM force field an additional potential function is added to the angle, in

Fglob =

∑ fi i=1

=

∂V ∂x

(7)

where V is the total potential energy of the system and nele is the number of all the potential energy finite elements. 2.2. Tangent Stiffness Matrix of the Generalized Finite Elements. If vector X represents the initial configuration and vector x denotes the final configuration, then the displacement field can be considered as D


Article

Journal of Chemical Information and Modeling U=x−X dU = d x

projections of the infinitesimal displacements at the global coordinate system of nodes 1 and 2 (i.e., J and I) along the axis that is defined with the unit vector λ1, du1 and du2 are the infinitesimal displacement vectors at the global coordinate system of nodes 1 and 2 (i.e., J and I), and dU is the infinitesimal displacement vector at the global coordinate system. For the second part of eq 9, it can be written that

(8) 45

In order to derive the tangent stiffness matrix for every energy term (i.e., for every potential energy finite element described in the previous section), the force vector fi is differentiated as follows: Ki =

r

d fi ∂f ∂n = i i + dU ∂ni ∂U

∑ k=1

∂fi ∂ck + ∂ck ∂U

m

∑ j=1

∂fi ∂λj ∂λj ∂U

r

∑

(9)

k=1

∂fi ∂ck =0 ∂ck ∂U

(11)

th

where Ki is the global tangent stiffness matrix for the i element, ni is the generalized internal force (moment Mi or axial force Ni) for the ith element, ck is the corresponding term of the r × 1 vector c, and r is the number of DOFs of the element’s local system. Depending on the potential energy finite element (two atom, angle, dihedral, etc.) the number of DOFs varies (r = 2 for the case of two atom interaction element, r = 4 for the in plane angle element, and r = 6 for the dihedral angle element); these terms express the equilibrium between internal and nodal forces for each element. λj is the corresponding element axis vector, m is the number of element axis vectors (m = 1 for the inline element, m = 2 for the inplane angle element, and m = 2 for the dihedral angle element), and U is the element’s displacement vector at the global system. The tangent stiffness matrix for the ith element is comprised of three parts; the first one is equivalent to the structural mechanics stiffness matrix when considering small displacements. Here the projection of the global displacements to the conjugate forces comes into play. The other two terms come into play when there are significant internal forces present (which is the case for molecular assemblies). The second one expresses changes in equilibrium due to variation of the displacement field, and the third one changes to geometry. 2.2.1. Bond, Nonbonded, and Urey−Bradley Potential Energy Finite Elements. For the case of interaction between two atoms (bond, nonbonded and Urey−Bradley terms) the differentiation of the corresponding internal force is expressed as follows:

For the third part of eq 9, two infinitesimal rotations are considered. As it is well established, the unit vector derivative expresses change in orientation and is perpendicular to the unit vector. An orthonormal basis matrix [λ1 λ2 λ3] is constructed and the projection of the global displacements at the appropriate axes λ3 and λ2 is considered to cause rotations dθ2 and dθ3, respectively. Therefore, it can be written that dλ1 = (λ 2 × λ1) ·dθ2 + (λ3 × λ1) ·dθ3 =

= B ·dU m

∑ j=1

3

⎡− B ⎤ K i = LT ·(kb̅ ·c T ·c) ·L + Ni·⎢ ⎥ ⎣B⎦

3

∑ λ1i ·du1i)

i=1

i=1

(10a)

⎡ du11 ⎤ ⎢ ⎥ ⎢ du12 ⎥ ⎢ ⎥ 0 ⎤ ⎢ du13 ⎥ ⎥·⎢ ⎥ λ13 ⎥⎦ ⎢ du 21 ⎥ ⎢ du ⎥ ⎢ 22 ⎥ ⎢⎣ du 23 ⎥⎦

Τ

Τ

= L ·(kb̅ ·c ·c)·L = L ·KLoc·L

(13)

dMi = kθ̅ ·dθ = kθ̅ ·c ·L ·d U dMi = kθ̅ ·c ·L dU

(14)

where dU is the infinitesimal displacement vector at the global coordinate system. The first part of the stiffness matrix is given by

= kb̅ ·c ·L ·d U ∂f ∂Ni dNi dN = kb̅ ·c ·L and i = LΤ·c Τ · i ∂Ni ∂U dU dU Τ

(12)

2.2.2. In-Plane Angle Potential Energy Finite Element. For the angle bond interaction (three atoms) there are two element axis, namely λ1 and λ2. For the case of interaction between three atoms (in plane angle term), the differentiation of the corresponding internal force is expressed as follows:

while dNi can also be written as

⎡ λ11 λ12 λ13 0 0 dNi = kb̅ ·[− 1 1]·⎢ ⎢⎣ 0 0 0 λ11 λ12

⎡−I3 × 3 ⎤ ∂λj ⎡− B ⎤ ∂fi ∂λj ⎥· = Ni·⎢ ⎥ = Ni·⎢ ⎣B⎦ ∂λj ∂U ⎣ I3 × 3 ⎦ ∂U

where ( × ) denotes the vector cross product; dλ1 is the infinitesimal change of the λ1 axis, considered to be a superposition of the two infinitesimal rotations dθ2 and dθ3 around axes λ2 and λ3, respectively; and B is the 3 × 6 matrix expressing the element axis derivative (λ1) with respect to the global displacements of nodes I and J. Therefore, based on the three parts of eq 9 provided by eqs 10b, 11, and 12, the global tangent stiffness matrix for the two-atom interaction (bonded or nonbonded) is given by the following expression:

dNi = kb̅ ·(dub2 − dub1) = kb̅ ·(∑ λ1i ·du 2i −

⎛1 T T ⎜ (λ × λ ) · [(λ × λ ) − (λ × λ ) ] 1 1 2 1 2 ⎝b 2 ⎞ 1 + (λ3 × λ1) ·[(λ1 × λ3)T − (λ1 × λ3)T ]⎟ ·dU ⎠ b

∂fi ∂ni dM = LΤ·c Τ · = LΤ·(kθ̅ ·c Τ ·c) ·L = LΤ·KLoc·L dU ∂ni ∂U

(10b)

(15)

The first part of eq 9 is derived in eqs 10, where dNi is the infinitesimal change of the internal force, dub1 and dub2 are the

The second term can be derived by differentiating c: E


Article


Figure 4. Part of the stiffness matrix for the in-plane angle element due to rotations of the plane. r

∑ k=1

∂Fi ∂ck = ∂ck ∂U

r

∑ k=1

∂Fi ⎛ ∂ck ∂lα ∂c ∂lβ ⎞ ⎜⎜ ⎟⎟ + k ∂ck ⎝ ∂lα ∂U ∂lβ ∂U ⎠

m

∑ j=1

∂Fi ∂λj = ∂λj ∂U

m

∑ j=1

2 ∂λ ∂λj ∂φ ⎞ ∂λ ∂ω ∂Fi ⎛ ∂θ ⎜⎜∑ j l + j ⎟ + ∂φ ∂U ⎟⎠ ∂ω ∂U ∂λj ⎝ l = 1 ∂θl ∂U

⎡ c1·I3 × 3 ⎤ ⎡ 03 × 3 ⎤ ⎢ ⎥ ∂Fi ⎢ ⎥ ∂Fi = ⎢ c 2·I3 × 3 ⎥ , = ⎢ 03 × 3 ⎥ , ∂λ1 ⎢ ⎥ ∂λ 2 ⎢ ⎥ ⎣ c3·I3 × 3 ⎦ ⎣ c4·I3 × 3 ⎦

∂Fi ∂l = M ·LT , a = [ λαT 01 × 3 −λαT ], ∂c ∂U ∂lβ = [01 × 3 λβΤ −λβΤ ] ∂U ∂c ∂c = [−1/lα2 0 1/lα2 0 ], = [0 1/lβ2 0 −1/lβ2 ] ∂la ∂lβ

∂λ ∂λ ∂λ1 ∂λ = (R1 × λ1), 1 = 0, 2 = (R1 × λ 2), 2 = 0 ∂θ2 ∂θ1 ∂θ1 ∂θ2

(16a)

(18a)

and in a more compact form: r

∑ k=1

and ⎡ λΤ ∂θ1 1 = [−1/la 1/la ]·⎢ ⎢0 ∂U ⎣

∂Fi ∂ck = M ·LT ·KELOC ·L 2 ∂ck ∂U

1×3

KELOC = ⎡ −1/l 2 0 1/lα2 α ⎢ ⎢ 0 1/lβ2 0 ⎢ ⎢ 2 −1/lα2 0 ⎢ 1/lα ⎢ −1/lβ2 0 ⎣ 0 ⎡ λT 0 ⎤ 1 × 3 01 × 3 ⎥ ⎢ α ⎢0 λ T 01 × 3 ⎥ ⎢ 1×3 β ⎥ L2 = ⎢ T ⎥ ⎢ 01 × 3 01 × 3 λα ⎥ ⎢ T ⎥ ⎣ 01 × 3 01 × 3 λβ ⎦

Τ ⎡0 01 × 3 ⎤ ∂θ2 1×3 λ2 ⎢ ⎥ = [−1/lβ 1/lβ ]· ⎢0 Τ ⎥ ∂U ⎣ 1 × 3 01 × 3 λ 2 ⎦

⎤ ⎥ 2 ⎥ −1/lβ ⎥, ⎥ 0 ⎥ ⎥ 1/lβ2 ⎦ 0

B1 = m

∑ j=1

(16b)

⎡ ⎤ c1·B1 2 ∂λ ⎞ ⎢ ⎥ ∂Fi ⎛ j ∂θl ⎜⎜∑ ⎟⎟ = M ·⎢ c 2· B2 ⎥ ∂λj ⎝ l = 1 ∂θl ∂U ⎠ ⎢ ⎥ ⎣ c3·B1 + c4·B2 ⎦

R3 =

λβ × λα |λβ × λα|

∂λ1 ∂θ1 ∂λ ∂θ , B2 = 2 2 ∂θ1 ∂U ∂θ2 ∂U

(18b)

For the λ1 and λ2 axis derivatives the definition of infinitesimal rotations can be considered; the infinitesimal rotations dθ1 and dθ2 are caused by the global displacement projections on the λ1 and λ2 axes at the start and end nodes of vectors JI and JK, respectively (see Figure 2). For the second and third parts, the rotation of axes λ1 and λ2 due to rotations of the plane is considered; thus, the displacements perpendicular to the plane uKproj, uIproj, and uJproj are of importance. These displacements can be decomposed into three parts, one translation and two rotations along two perpendicular axes, i.e. dφ and dω along R2 (unit vector KI; see Figures 2 and 4) and R3, respectively, where

where lα and lβ are the length of vectors λα and λβ, respectively. For the third term, the change of the orientation of axis λ1 and λ2 is of importance. As this is an in-plane angle element, at an imposed displacement field, the plane will have a new position. One way to express the changes of the element axis is relative to the plane. The change is considered in two terms; the first term is due to in plane displacements. These displacement cause λ1 and λ2 to rotate in plane along the plane vector,

R1 =

01 × 3 01 × 3 ⎤ ⎥ Τ ⎥ 01 × 3 λ1 ⎦

R1 × R 2 |R1 × R 2|

(19)

which are considered to have as origin the center C0 of the parallelogram of the triangle’s ABC coordinates (see Figure 4). The second and third parts are calculated by the following relations:

(17)

by dθ1 and dθ2, respectively. The second term is due to rotations ω and φ of the plane. F


Article


⎡ ⎤T 1 1 ∂Fi ∂C 0 0 0 0⎥ = M · LT , = ⎢− 2 2 ⎣ la sin θa la sin θa ⎦ ∂la ∂C

∂λ1 ∂λ = (R 2 × λ1), 1 = (R3 × λ1), ∂φ ∂ω ∂λ 2 ∂λ = (R 2 × λ 2), 2 = (R3 × λ 2) ∂φ ∂ω

T ⎡ 1 1 ⎤ ∂C 0 2 = ⎢0 0 0 − 2 ⎥ lc sin θc lc sin θc ⎦ ⎣ ∂lc

∂λ1 ∂φ ∂λ ∂ω ∂λ ∂φ ∂λ ∂ω + 1 = Bl1 and 2 + 2 = Bl 2 ∂φ ∂U ∂ω ∂U ∂φ ∂U ∂ω ∂U

⎡ ⎤T 1 1 1 1 ∂C 0 − 2 − 2 = ⎢0 2 ⎥ lb tan θa lb tan θc lb 2 tan θc ⎦ ⎣ lb tan θa ∂lβ

(20a)

thus m

∑ j=1

∂la = [ λαΤ − λαΤ 01 × 3 01 × 3 ], ∂U ∂lc = [ 01 × 3 01 × 3 − λcΤ λcΤ ] ∂U

⎤ ⎡ c1·Bl1 ⎥ ⎢ ∂λj ∂φ ⎞ ∂Fi ⎛ ∂λj ∂ω c 2·Bl 2 ⎥ + ⎜ ⎟ = M·⎢ ∂φ ∂U ⎠ ∂λj ⎝ ∂ω ∂U ⎥ ⎢ ⎣ c3·Bl1 + c4·Bl 2 ⎦

∂lβ ∂U

(20b)

= [01 × 3 λβΤ − λβΤ 01 × 3 ] (23b)

For an analytical calculation of matrices Bl1 and Bl2, see the Supporting Information. The global tangent stiffness matrix for the angle bond interaction is given by the following expression: K i = LΤ·(K theta·CΤ·C) ·L + M ·LT ·KELOC ·L 2 ⎡ ⎤ c1·(B1 + Bl1) ⎢ ⎥ ⎥ c 2·(B2 + Bl 2) + M·⎢ ⎢ ⎥ ⎢⎣ c3·(B1 + Bl1) + c4·(B2 + Bl 2)⎥⎦

and

∂C ∂θa ⎡ cos θa ⎤T cos θa 1 1 0 0 0⎥ + − = ⎢− 2 2 2 2 lβ sin θa lβ sin θa ⎢⎣ la sin θa la sin θa ⎥⎦

(21)

2.2.3. Dihedral Angle Potential Energy Finite Element. For the case of interaction between four atoms (proper and improper dihedral terms), the differentiation of the corresponding internal force is expressed as follows: There are two element axes normal to each plane, namely λ1 and λ2. The interatomic potential, the corresponding internal force (torque), and the first term of the stiffness matrix are dM =

∂ 2Vdihedral ∂chi 2

∂C ∂θc ⎡ ⎤T cos θc cos θc 1 1 0 0 0 − − + ⎥ ⎢ = lc sin 2 θc lβ sin 2 θc lc sin 2 θc lβ sin 2 θc ⎥⎦ ⎢⎣

dchi = Kchi·C·L ·dU

∂M = Kchi·C·L ∂U ∂Fi ∂ni = LΤ·(Kchi·CΤ·C) ·L = LΤ·KLoc ·L ∂ni ∂U

∂θa Τ Τ = [(1/la)λΝΑ 1 − (1/ la)λΝΑ1 01 × 3 01 × 3 ] ∂UKL ∂θc Τ Τ = [ 01 × 3 01 × 3 − (1/lc)λΝΑ 4 (1/ lc)λΝΑ4 ] ∂UIJ

(22)

where dU is the element infinitesimal global displacement. The second term can be derived by differentiating C with respect to la, lβ, lc, θa, and θc or by realizing that, for displacements along each plane (see Figure 3), the internal torque should not change and should appropriately define the changes dla, dlβ, dlc, dθa, and dθc. The second term is divided into two terms: one for changes of la, lβ, and lc and one for changes of θa and θc (see Figure 3). Due to the fact that the change of the level arm of the forces at nodes L and I (see Figure 3) is what is of interest herein, the angle derivatives should be expressed regarding the displacements of the appropriate arm of each angle. r

∑ k=1

∂Fi ∂ck = ∂ck ∂U

+

r

∑ k=1

(23c)

and in a more compact form:

r

∑

∂Fi ⎛ ∂ck ∂la ∂c ∂lβ ∂c ∂l ⎜⎜ + k + k c ∂ck ⎝ ∂la ∂U ∂lβ ∂U ∂lc ∂U

∂ck ∂θa ∂c ∂θβ ⎞ ⎟ + k ∂θa ∂UKL ∂θβ ∂UIJ ⎟⎠

k=1

∂Fi ∂ck = M ·LT ·Θ ∂ck ∂U

Θ = KELOC1·LA + KELOC 2·LB

(23d)

(23a)

where

where G


Article

Journal of Chemical Information and Modeling ⎡−KE1 KE1 0 0 0 0 ⎤ ⎢ ⎥ 0 0 ⎥ ⎢ KE1 −KE1 KE2 −KE2 ⎢ ⎥ −KE2 KE2 0 0 0 ⎥ ⎢ 0 KELOC1 = ⎢ , KE3 −KE3 ⎥ 0 0 0 0 ⎢ ⎥ ⎢ 0 −KE4 KE4 0 0 0 ⎥ ⎢ ⎥ ⎢⎣ 0 KE4 −KE4 −KE3 KE3 ⎥⎦ 0 ⎡ λT ⎢ a ⎢0 ⎢ 1×3 ⎢0 1×3 LA = ⎢ ⎢ ⎢ 01 × 3 ⎢ ⎢ 01 × 3 ⎢ ⎣ 01 × 3

01 × 3 01 × 3 01 × 3 ⎤ ⎥ T λa 01 × 3 01 × 3 ⎥ ⎥ T λb 01 × 3 01 × 3 ⎥⎥ ⎥ 01 × 3 λbT 01 × 3 ⎥ ⎥ 01 × 3 λcT 01 × 3 ⎥ ⎥ 01 × 3 01 × 3 λcT ⎦

For the third term the unit vectors define the two planes of the dihedral angle. The unit vector derivatives are derived from the dihedral unit vector (element’s axis) definition. The unit vector derivative expresses orientation change, and the rotations of each of the two planes are of importance; thus, displacements perpendicular to each plane are considered. ⎡ C1·I3 × 3 ⎤ ⎡ 03 × 3 ⎤ ⎢ ⎥ ⎢ ⎥ ⎢C2·I3 × 3 ⎥ ∂Fi ⎢C5·I3 × 3 ⎥ ∂Fi = M·⎢ = M⎢ ⎥, ⎥ ∂λ1 ⎢ C3·I3 × 3 ⎥ ∂λ4 ⎢C6·I3 × 3 ⎥ ⎢ 0 ⎥ ⎢C ·I ⎥ ⎣ 3×3 ⎦ ⎣ 4 3×3 ⎦

also dλ1 = −

1 ((b1 + u11 − u 21) × (b2 + u31 − u 21)) |b1′ × b2′| 1 b1 × b2 + |b1 × b2|

(23e)

1 1 1 , KE2 = 2 , KE3 = 2 , la2 sin(θα) lb tan(θα) lc sin(θc) 1 KE4 = 2 lb tan(θc)

KE1 =

≅− =

KEB3 =

cos(θc) lc2 sin 2(θc)

⎡ λT ⎢ NA1 ⎢0 1×3 LB = ⎢ ⎢0 ⎢ 1×3 ⎢ ⎣ 01 × 3

, KEB2 =

, KE4 =

(25b)

and dλ4 ≅ =

1 , lalb sin 2(θα)

(25c)

1 lclb sin 2(θc)

thus

01 × 3 01 × 3 01 × 3 ⎤ ⎥ T λNA1 01 × 3 01 × 3 ⎥ ⎥ T 01 × 3 λNA 4 01 × 3 ⎥⎥ T ⎥ 01 × 3 01 × 3 λNA 4⎦

m

∑ j=1

T ∂Fi ∂λj ∂λj ∂U

⎡ ⎤ c1·Bl1 ⎢ ⎥ ⎢ c 2·Bl1 + c5·Bl4 ⎥ = M·⎢ ⎥ ⎢ c3·Bl1 + c6·Bl4 ⎥ ⎢ ⎥ c4·Bl4 ⎣ ⎦

(25d)

where b1, b2, and b3 are the vectors defined by the 4 atoms (b1 = −λα, b2 = −λβ, and b3 = λc) in the initial configuration and M1, M2, and M3 are the matrices expressing the cross product of each of these vectors (skew-symmetric matrices), b1′ , b2′ , and b3′ are the vectors defined by the 4 atoms in the final configuration, and uXY are the projections of the global displacements of node X (1 for node L, 2 for K, 3 for J, and 4 for I) at vector Y (1 for λ1, 4 for λ4), while the matrices Bl1 and Bl4 express the unit vector derivatives λ1 and λ4, respectively. Finally, as dU is small, it is considered that |b′1 × b′2| ≅ |b1 × b2| and |b′2 × b′3| ≅ |b2 × b3|. The global tangent stiffness matrix for the dihedral bond interaction is given by

λ1 × λα |λ1 × λα| λ4 × λc |λ4 × λc|

1 (M 2·λ4 ·([01 × 9 λ4T ]·dU − [01 × 6 λ4T 01 × 3 ]·dU ) |b2 × b3|

= Bl4 ·dU

and λNA 4 =

1 [b2 × (u44 − u34) + (u34 − u 24) × b3] |b2 × b3|

− M3·λ4 ·([01 × 6 λ4T 01 × 3 ]·dU − [01 × 3 λ4T 01 × 6 ]·dU ))

where λα is the vector KL, λβ is the vector JK, and λc is the vector JI (see Figure 3): λNA1 =

1 (M1·λ1·( −[ 01 × 6 λ1T 01 × 3 ]·dU |b1 × b2| + [ 01 × 3 λ1T 01 × 6 ]·dU ) −M 2 ·λ1·

= Bl1·dU

(23f)

cos(θα)

1 [b1 × (u31 − u 21) + (u11 − u 21) × b2] |b1 × b2|

( −[ λ1T 01 × 9 ]·dU − [ 01 × 3 λ1T 01 × 6 ]·dU ))

KELOC2 = ⎤ ⎡ − KEB1 KEB1 0 0 ⎥ ⎢ ⎥ ⎢(KEB1+KEB2 ) − (KEB1+KEB2 ) 0 0 ⎥ ⎢ KEB2 0 0 ⎥ ⎢ − KEB2 ⎥ ⎢ 0 0 KEB KEB − 3 3 ⎥ ⎢ ⎢ 0 0 KEB4 − KEB4 ⎥ ⎥ ⎢ ⎢⎣ 0 0 − (KEB3+KEB4 ) (KEB3+KEB4 )⎥⎦

la2 sin 2(θα)

1 1 b1′ × b2′ + b1 × b2 ′ ′ |b1 × b2| |b1 × b2|

=−

and

KEB1 =

(25a)

(24) H


Article


Figure 5. Schematic representation for the element assembly of the global stiffness matrix for the nonbonded and Urey−Bradley interactions.

⎡ ⎤ c1·Bl1 ⎢ ⎥ ⎢ c 2·Bl1 + c5·Bl4 ⎥ T Τ Τ K i = L ·(Κchi·C ·C) ·L + M ·L ·Θ + M ·⎢ ⎥ ⎢ c3·Bl1 + c6·Bl4 ⎥ ⎢ ⎥ c4·Bl4 ⎣ ⎦

2.3. On the Calculation of the Hessian Matrix as a Finite Element Assembly. Based on the proposed generalized formulation, for each component of the energy terms (both bonded and nonbonded ones), a specialized finite element is formulated composed by nodes equal to the number of atoms (a node is located at the position of each atom) expressing the specific component of the energy terms. The total potential energy can be expressed as a summation of the energy terms:

(26)

With appropriate summation at the corresponding degrees of freedom for each element,35 the energy Hessian matrix is obtained: nele

K glob =

∑ Ki = i=1

∂ 2V ∂x 2

nele

V=

∑ Vi i=1

(27)

Vi = Vi (zi) = Vi (Xi) ∂V = ∂X

The element assembly for the formation of the global stiffness matrix is not restricted only to the bonded atoms. The same summation procedure is followed for the nonbonded terms as well as for the Urey−Bradley terms involved in the angle bond interaction among three atoms. A schematic representation of the global element assembly for such cases is presented in Figure 5.

nele

∑ i=1

∂Vi ∂X

⎧ ∂Vi ∂zi , for x = xi ⎪ ∂Vi = ⎨ ∂zi ∂Xi ∂X ⎪ ⎩ 0, otherwise I

(28) DOI: 10.1021/acs.jcim.6b00356 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling where nele is the number of energy terms (or energy elements) and zi is the corresponding internal degree of freedom for the ith energy term, while Xi are the corresponding Cartesian coordinates for the ith energy term (zi is a function of Xi). As can be recognized for the expressions of eq 28, the gradient can be obtained as the proper summation for all energy elements. Likewise, for the Hessian matrix it can be written that ∂ 2V = ∂X2

nele

∑ i=1

burden of a naive numerical computation, the connectivity of each degree of freedom was taken into account when computing energy pertrubations. In order to validate the proposed modeling approach and asses its efficiency in the case of energy minimization, the following molecular systems were studied: 2.4.1. Single Butane Molecule. A single butane molecule (Figure 6) is used as a proof-of-concept system, where the

∂ 2Vi ∂X2

⎧ ∂Vi ∂zi ∂ ⎪ ⎪ ∂zi ∂Xi ∂ Vi , for x = xi ⎨ = ∂X2 ⎪ ∂X ⎪ 0, otherwise ⎩

(

2

(

∂

∂Vi ∂zi ∂zi ∂X i

)

) = ⎛⎜ ∂ V ∂z ⎞⎟

T ∂zi ∂V i i + i 2 X X ∂X ∂ ∂ ∂zi z ∂ ⎝ i ⎠ i 2 ⎛ 2 T ⎞ ∂ Vi ∂zi ∂zi ∂V ∂ zi ⎟+ i ⎜ = 2 ∂zi ∂Xi2 ∂zi ⎝ ∂Xi ∂Xi ⎠ 2

( )

∂

∂zi ∂X i

∂X

thus ⎧ ∂ 2V ⎛ ∂z T ∂z ⎞ ∂V ∂ 2z i i ⎪ i⎜ i ⎟+ i , for 2 2 ⎪ X X z ∂ ∂ ∂ z X ∂ ∂ ⎠ ⎝ i i i i i 2 ⎪ ∂ Vi 2 = ⎨ ⎛ ∂ Vi ⎞ 2 ⎜ ⎟ ∈ Xi ⊗ Xi ∂X ⎪ 2 X ∂ ⎝ ⎠mn ⎪ ⎪ 0, otherwise ⎩

Figure 6. Butane molecule.

number of nodes totalled 14 (4 carbon atoms and 10 hydrogen atoms) and the corresponding degrees of freedom of the proposed FE model are 42. The number of bonds are 13, those of angle and Urey−Bradley bonds are 24, those of dihedral angle bonds 27, and those of nonbonded terms 54. The single butane molecule case was considered in order to validate the analytical expressions presented herein in order to calculate the energy gradient and Hessian matrix. Atomic parameters from the CHARMM2743,44 field were used. 2.4.2. Double Stranded B-DNA. The second system investigated was a canonical B-form DNA duplex. The number of nodes (or atoms) for AT (adenine-thymine) and TA bases are 64, and those for CG (cytosine-guanine) and GC bases are 63 (Figure 7 denotes graphically the nucleobases of DNA). However, in the model used here, we used 41 for AT and CG as the hydrogen atoms where excluded and the hydrogen bonds among the bases described explicitly using a Morse potential. The number of bonds are 61 for CG and 62 for AT, the number of angle bonds are 58 for CG and AT, and the number of dihedral angle bonds are 59 for CG and 56 for AT; finally, the number of improper dihedral angle bonds are 4 for CG and AT. Regarding the nonbonded terms that express Coulombic and Lenard-Jones interactions, a cut off radius of 8 Å was used. DNA cases of 2, 10, and 30 base pairs were used to compare the performance of the method when incorporated in energy minimization. The solution of equilibrium equations was used in order to derive the sequence dependent stretch modulus for indicative base pair steps. Atomic parameters from the CHARMM27 field43,44 were used. In the case of DNA, one of the authors has seen successful use of the CHARMM field to study DNA nanostructure conformations.37,39,46

(29)

A combinatorial comparison to a conventional numerical comutation for full nonbonded interaction cases is presented below. The main argument is to demonstrate that, in this case, the computation cost is of the same order of that of the energy gradient. In the full connectivity case (upper limit), each atom is connected to at least (n − 3)/3 other atoms via nonbonded interactions, where n is the total number of configurational degrees of freedom. In this case, the number of energy terms (i.e., number of elements) is nele ∼ n2 (where ∼ denotes that nele is proportional to n2). Let O(NumHess) denote the number of operations for the conventional numerical computation of the Hessian Matrix. If k operations are performed for each component of the Hessian matrix, hence the total number of operations: O(NumHess) ∼ kn2. The number of operations k (energy perturbations) are not fixed when considering the full nonbonded interaction. The operations for the calculation of the energy difference of the connected atoms is proportional to n; thus, O(NumHess) ∼ n3. If O(PEEM) are the number of operations for the proposed method where the Hessian matrix is computed as an element assembly of nele ∼ n2 energy terms, it is proved that according to the proposed modeling approach, Hessian matrix computing requirements O(PEEM) ∼ n2 are of the same order to that of gradient and energy computations for a full nonbonded connections case. In this context the use of an analytically computed Hessian matrix should be reconsidered when dealing with energy minimization problems, since fast and efficient algorithms can be developed. 2.4. Computational Details and Test Systems. All the computations were performed in Matlab version 15b.45 The Hessian numerical computation was also done in Matlab by implementing central differences, and in order to reduce the

3. RESULTS AND DISCUSSION 3.1. Validation. For the butane molecule case, the energy gradient and Hessian matrix were computed both numerically and analytically, implementing the proposed modeling approach. For the numerical differentiation, the centered J


Article


Figure 7. DNA bases (a) adenine, (b) thymine, (c) cytosine and (d) guanine.

differencing technique was used, aiming to achieve increased accuracy. The difference between the gradient calculated numerically and that achieved by means of the assembly of internal force vector is practically equal to zero (the larger difference is equal to 5 × 10−07%), proving the equivalence of the internal force analysis with the energy differentiation. The energy Hessian matrix computed numerically and the tangent global stiffness matrix are identical (the differences are of the order of E−06% while the large elements of the matrices are of the order of E+03; see Figures 8a and 8b for the graphical presentation of the form of the Hessian matrix and that of the differences). In order to verify the equivalence of the two matrices (calculated numerically and analytically), the following energy comparison was performed. The energy is approximated locally as a quadratic function of displacement. Thus, in an arbitrary direction, the energy can be computed by 1 V1 = V + g T ·U + UT ·H ·U (30) 2 where V and V1 are the potential energy at the initial and final configurations, respectively, U is an arbitrary displacement vector while g and H are the energy gradient vector and Hessian matrix, respectively. The measure for comparing the Hessian matrices considered herein is based on the maximum difference (ρ) calculated along every possible direction that is defined according to the following expression: ρ = max

UT ·H ·U − UT ·K·U ,∀U UT ·H ·U

Figure 8. Butane molecule-performance of the proposed modeling approach: (a) graphical representation of the Hessian Matrix and (b) element by element difference of Hessian matrices (analytical calculation obtained based on the new modeling approach and numerically computed).

effort required for computing the Hessian matrix. The CPU time required for computing the Hessian matrix based on the proposed analytical calculation and that of the numerical computation for the case of double stranded DNA molecules (varying on the number of base pairs considered) is depicted in Figure 9. As was mentioned previously for the numerical computation, the central differencing technique for each DOF was applied while the energy differences were computed by involving only the relevant energy terms. As is already mentioned in section 4 and depicted in Figure 9, the proposed method, as well as any method that assembles the Hessian matrix as a summation from the matrices of second derivatives for every energy term, is significantly faster than the numerical computation. The method was integrated in the energy minimization algorithm in order to assess its performance when compared with the BFGS minimization scheme, which is one of the standards for this kind of system. The algorithm

(31)

where K is the tangent global stiffness matrix. The value of ρ for the butane molecule case is equal to 1.26 × 10−08%, proving the equivalence of the proposed method. As already mentioned, the proposed modeling approach resulting in analytical calculation of the energy gradient vector and Hessian matrix was significantly faster than the case of the numerical computation. 3.2. Performance in Terms of Speed and Energy Minimization. Initially, the performance of the proposed modeling approach was assessed in terms of computational K


Article


BFGS minimization does not use the entire Hessian matrix but rather approximates it by the results obtained through the successive iterations. In the BFGS scheme used for this comparative study, the analytical gradients (residual forces) were provided, so the method could be comparable in terms of speed (numerical gradients demand an appropriate energy computation in each DOF involving only the relevant energy terms; otherwise, they are much slower). The results of energy minimization implementing both schemes in three double stranded DNA molecules (composed by 2, 10, and 30 base pairs) are presented in Figure 10. Both methods were approximately equivalent in terms of the speed of each iteration. However, the BFGS minimization scheme did not manage to reach a satisfactory minimization level (the minimized configuration contains large remaining forces). For both minimization schemes, the maximum number of iterations allowed was 1000. In all cases, BFGS terminated when the number of maximum iterations was reached. The trust-region algorithm, integrating the Hessian matrix and gradient vectors calculated in accordance with the proposed method, performed significantly better in terms of minimization level (the remaining forces were almost equal to zero in all three cases) as well as in the rate of convergence (number of iterations needed to reach the minimum). In all three cases examined here, the Hessian matrix was storable; however, in cases where such an option is not feasible, a sparse version of the Hessian matrix was also examined. For the three DNA cases, a banded version of the tangent stiffness matrix (Hessian) was used. The 11 main diagonals were stored and used in the minimization process. This amount of memory is, for the 2 base pair case, 4.8% of the matrix, in the 10 base pair case, 0.97% of the matrix, and in the 30 base pair case, 0.32% of the matrix. In all the cases, the scheme implemented using the banded form of the Hessian matrix performed better in terms of residual forces than the one using the BFGS method, as can be seen in Tables 1−3. 3.3. Equivalent Finite Element Simulation and Modulus of Elasticity of Double Stranded DNA. The proposed analysis is in essence a structural molecular model.

Figure 9. DNA molecules-performance of the proposed modeling approach (i.e., analytical calculation based on the proposed FE type approach) vs numerical calculation of the Hessian matrix in terms of computation effort.

integrating the proposed method is of the trust-region family.47 In general, the problem is stated as

{ 12 s Hs + s g , such that Ds ≤ Δ}

min

T

T

(32)

where Δ is the trust-region dimension. In the particular variant48 used, the trust-region subproblem is restricted into a two-dimensional subspace S {s1, s2} where s1 is the direction of the gradient g and s2 is defined by

Hs2 = −g or s2T Hs2 < 0

(33)

and obtained by solving iteratively via the preconditioned conjugate gradient algorithm.

Figure 10. Thirty base pair double stranded DNA: (a) minimized; (b) minimized under load on the top in y−y′ direction and fixed at the base. L


Article

Journal of Chemical Information and Modeling Table 1. Minimization Results50 for Double Stranded DNA of DNA Molecule of 2 Base Pairs (CT)

No. of Iterations Potential energy value Max. Residual Force (First order optimality)

BFGS

Trust-Region (Full Hessian)

Trust-Region (Banded Hessian)

1001 2.86 × 1004 13.425

237 2.85 × 1004 9.21 × 10−06

1001 2.86 × 1004 6.689

(gradient equals 0). Then external forces (F) are applied at the top two edges of the DNA molecule backbone (see Figure 11 and Table 4) while the DOFs of the bottom base pair are fixed along the direction of F along with all the DOFs (x,y,z) of the connecting to the backbone atoms of the bottom base pair. In order to calculate the mean distance (l ̅) between the base pairs and the mean displacement (u)̅ of the top base pair, two schemes are implemented. In the first one, l ̅ is the difference of the mean values of the z-coordinate of the base pairs while u̅ is the mean value of the z-direction displacement of the top base pair (with the bottom base pair being fixed on this direction). In the second, l ̅ is the mean value of the z-coordinate of the distance between the atoms connecting the backbone for top and bottom base pairs; for example, the atoms which are fixed and the atoms where the external forces are applied and as displacement the mean value of the z-direction displacement of the two atoms where the external forces are applied. As the stretch modulus is of interest here, there is no need for nonlinear minimization; thus, solving the linear system in eq 36 is straightforward.

Table 2. Minimization Results50 for Double Stranded DNA of DNA Molecule of 10 Base Pairs (AGTCTAGTCT)


BFGS



1001 1.466 × 1005 68.71

285 1.462 × 1005 7.43 × 10−06

1001 1.466 × 1005 16.10

K mU = F

Table 3. Minimization Results50 for Double Stranded DNA of DNA Molecule of 30 Base Pairs (AAGGGCCTACGTACGATTCGAGTCTAGTCT)


BFGS



1001 4.518 × 1005 162.01

491 4.500 × 1005 2.53 × 10−05

1001 4.518 × 1005 35.10

where Km is the reordered stiffness matrix. The stretch modulus EA is computed by the following expression:

Fl ̅ (37) u̅ The values for the stretch modulus for four different base pair steps are depicted in Table 2. The exact stretch modulus is not independent of the profile of the external forces; for the four cases examined herein, the external forces were applied at the edges of the backbone of the top base pair. EA =

4. CONCLUSIONS In this work a new modeling approach for simulating the potential energy interactions of molecular nanostructures is proposed. Finite element method based approaches are usually employed in solving differential equations in continuum media. The proposed numerical modeling approach devises the finite element method in order to develop a discrete model that accurately captures the interactions at the atomic level. The new modeling approach relies on the development of generalized potential energy finite elements and allows taking into account all types of force fields and parameter sets used to calculate the potential energy of a system of atoms or coarsegrained particles. The analytical expression of both gradient vector and Hessian matrix is the main advantage of the proposed modeling approach, leading to drastic reduction of the computational effort required. Two molecular systems are considered in order to assess the performance of the proposed modeling approach: a single butane molecule and double stranded B-DNA sequences (composed of 2, 10, and 30 basepairs). The fast and efficient performance of the method can find applications in directly improving the energy minimization algorithms, especially in cases of rigorous minimization. Furthermore, the sequence dependent stretching modulus for four DNA base pair steps was calculated. Finally, one other application of the proposed formulation when coupled with an appropriate coarse graining orientation scheme (orientation and displacement of groups of atoms) is the development of a systematic definition of coarse grained models. This could be possibly achieved by coupling already developed methodologies

The force and stiffness matrix analysis was used in the double stranded DNA model using the CHARMM force field. The deformed configuration for a double stranded DNA chain under external loads can be derived (neglecting entropic effects) by implementing the following analysis: As the tangent stiffness matrix is a structural model, arbitrary DOFs can be “fixed” by eliminating the corresponding lines and columns and reordering. Nodal loads are considered, and the potential energy for minimization is E = V − FT ·U

(34)

where E is the potential energy, V is the internal potential energy, F is the nodal loads vector, and U is the displacement vector. At equilibrium, as the potential energy exists as a function of displacement and as Castigliano’s first theorem implies: ∇E = ∇V − F = 0

(36)

(35)

The double stranded DNA sequence was minimized at extensive loads of increasing magnitude. The minimization procedure is nonlinear and equivalent to energy minimization usually performed before molecular dynamics analysis, though much more rigorous, as the energy must be at an absolute minimum (see Figures 9a−9b). 3.4. Sequence Dependent DNA Stretch Modulus. For the derivation of the modulus of elasticity in vacuum and low temperature (neglecting entropic effects), a set of 4 base pair steps is examined. Initially, energy minimization is performed in order to render the global stiffness matrix positive-definite (eigenvalues greater or equal to 0) and vanish external forces M


Article


Figure 11. DNA base pair steps along with deformed configurations (a) AA, (b) AT, (c) GC, and (d) GG.

■

Table 4. Stretch Modulus for Four Different DNA Base Pair Steps along with Deformed Configurations51 DNA Base Pair Step

Stretch modulus (1st method) (pN)

Stretch modulus (2nd method) (pN)

AA AT GC GG

1182 2785 2605 3525

533 1108 1235 429

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The first and third authors would like to thank the people of the Laboratory for Computational Biology and Biophysics of the Department of Biological Engineering at MIT for sharing important aspects of the fascinating field of DNA nanotechnology and for their help in the development of the presented work.

of static reduction in the field of structural mechanics with existing methods of absolute orientation49 that can help determine the “stiff” parts of the nanostructure to be condensed.

■

AUTHOR INFORMATION

Corresponding Author

■

ASSOCIATED CONTENT

S Supporting Information *

REFERENCES

(1) Agrawal, P. M.; Sudalayandi, B. S.; Raff, L. M.; Komanduri, R. Molecular Dynamics (MD) Simulations of the Dependence of C-C Bond Lengths and Bond Angles on the Tensile Strain in Single-Wall Carbon Nanotubes (SWCNT). Comput. Mater. Sci. 2008, 41, 450− 456. (2) Chen, W. H.; Cheng, H. C.; Hsu, Y. C. Mechanical Properties of Carbon Nanotubes Using Molecular Dynamics Simulations with the

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.6b00356. Analytical calculation of matrices Bl1 and Bl2 used for deriving the in-plane ange potential energy finite element (PDF) N


Article

Journal of Chemical Information and Modeling Inlayer Van Der Waals Interactions. Comput. Model. Eng. Sci. 2007, 20, 123−145. (3) Cheng, H.-C.; Liu, Y.-L.; Hsu, Y.-C.; Chen, W.-H. AtomisticContinuum Modeling for Mechanical Properties of Single-Walled Carbon Nanotubes. Int. J. Solids Struct. 2009, 46, 1695−1704. (4) Gao, G.; Ç aǧin, T.; Goddard, W. A., III Energetics, Structure, Mechanical and Vibrational Properties of Single-Walled Carbon Nanotubes. Nanotechnology 1998, 9, 184−191. (5) Hu, N.; Nunoya, K.; Pan, D.; Okabe, T.; Fukunaga, H. Prediction of Buckling Characteristics of Carbon Nanotubes. Int. J. Solids Struct. 2007, 44, 6535−6550. (6) Krasheninnikov, A. V.; Banhart, F.; Li, J. X.; Foster, A. S.; Nieminen, R. M. Stability of Carbon Nanotubes Under Electron Irradiation: Role of Tube Diameter and Chirality. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 125428. (7) Nasdala, L.; Ernst, G. Development of a 4-Node Finite Element for the Computation of Nano-Structured Materials. Comput. Mater. Sci. 2005, 33, 443−458. (8) Odegard, G. M.; Gates, T. S.; Nicholson, L. M.; Wise, K. E. Equivalent-Continuum Modeling of Nano-Structured Materials. Compos. Sci. Technol. 2002, 62, 1869−1880. (9) Wang, Y.; Sun, C.; Sun, X.; Hinkley, J.; Odegard, G. M.; Gates, T. S. Nano-Scale Finite Element Analysis of Polymer Networks. In 43rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Denver, Colorado, USA, April 22−25, 2002. (10) Yakobson, B. I.; Brabec, C. J.; Bernholc, J. Nanomechanics of Carbon Tubes: Instabilities Beyond Linear Response. Phys. Rev. Lett. 1996, 76, 2511−2514. (11) Merli, R.; Lázaro, C.; Monleón, S.; Domingo, A. A Molecular Structural Mechanics Model Applied to the Static Behavior of SingleWalled Carbon Nanotubes: New General Formulation. Comput. Struct. 2013, 127, 68−87. (12) Seeman, N. C.; Kallenbach, N. R. Design of Immobile Nucleic Acid Junctions. Biophys. J. 1983, 44, 201−209. (13) Rothemund, P. W. K. Folding DNA to Create Nanoscale Shapes and Patterns. Nature 2006, 440, 297−302. (14) Douglas, S. M.; Dietz, H.; Liedl, T.; Hoegberg, B.; Graf, F.; Shih, W. M. Self-Assembly of DNA into Nanoscale Three-Dimensional Shapes. Nature 2009, 459, 414−418. (15) Dietz, H.; Douglas, S. M.; Shih, W. M. Folding DNA into Twisted and Curved Nanoscale Shapes. Science 2009, 325, 725−730. (16) Han, D.; Pal, S.; Nangreave, J.; Deng, Z.; Liu, Y.; Yan, H. DNA Origami with Complex Curvatures in Three-Dimensional Space. Science 2011, 332, 342−346. (17) Zhang, F.; Jiang, S.; Wu, S.; Li, Y.; Mao, C.; Liu, Y.; Yan, H. Complex Wireframe DNA Origami Nanostructures with Multi-Arm Junction Vertices. Nat. Nanotechnol. 2015, 10, 779−784. (18) Geary, C.; Rothemund, P. W. K.; Andersen, E. S. A SingleStranded Architecture for Cotranscriptional Folding of RNA Nanostructures. Science 2014, 345, 799−804. (19) Delebecque, C. J.; Lindner, A. B.; Silver, P. A.; Aldaye, F. A. Organization of Intracellular Reactions with Rationally Designed RNA Assemblies. Science 2011, 333, 470−474. (20) Wei, R.; Martin, T. G.; Rant, U.; Dietz, H. DNA Origami Gatekeepers for Solid-State Nanopores. Angew. Chem., Int. Ed. 2012, 51, 4864−4867. (21) Dutta, P. K.; Varghese, R.; Nangreave, J.; Lin, S.; Yan, H.; Liu, Y. DNA-Directed Artificial Light-Harvesting Antenna. J. Am. Chem. Soc. 2011, 133, 11985−11993. (22) Pan, K.; Boulais, E.; Yang, L.; Bathe, M. Structure-Based Model for Light-Harvesting Properties of Nucleic Acid Nanostructures. Nucleic Acids Res. 2014, 42, 2159−2170. (23) Mohri, K.; Nishikawa, M.; Takahashi, N.; Shiomi, T.; Matsuoka, N.; Ogawa, K.; Endo, M.; Hidaka, K.; Sugiyama, H.; Takahashi, Y. Design and Development of Nanosized DNA Assemblies in PolypodLike Structures as Efficient Vehicles for Immunostimulatory CpG Motifs to Immune Cells. ACS Nano 2012, 6, 5931−5940.

(24) Davis, M. E.; Chen, Z.; Shin, D. M. Nanoparticle Therapeutics: An Emerging Treatment Modality for Cancer. Nat. Rev. Drug Discovery 2008, 7, 771−782. (25) Liu, M.; Fu, J.; Hejesen, C.; Yang, Y.; Woodbury, N. W.; Gothelf, K.; Liu, Y.; Yan, H. A DNA Tweezer-Actuated Enzyme Nanoreactor. Nat. Commun. 2013, 4, 2127. (26) Wadati, M.; Tsuru, H. Elastic Model of Looped DNA. Phys. D 1986, 21, 213−226. (27) Westcott, T. P.; Tobias, I.; Olson, W. K. Elasticity Theory and Numerical Analysis of DNA Supercoiling: An Application to DNA Looping. J. Phys. Chem. 1995, 99, 17926−17935. (28) Purohit, P. K.; Nelson, P. C. Effect of Supercoiling on Formation of Protein-Mediated DNA Loops. Phys. Rev. E 2006, 74, 061907. (29) Balaeff, A.; Mahadevan, L.; Schulten, K. Elastic Rod Model of a DNA Loop in the Lac Operon. Phys. Rev. Lett. 1999, 83, 4900−4903. (30) Balaeff, A.; Koudella, C. R.; Mahadevan, L.; Schulten, K. Modelling DNA Loops Using Continuum and Statistical Mechanics. Philos. Trans. R. Soc., A 2004, 362, 1355−1371. (31) Balaeff, A.; Mahadevan, L.; Schulten, K. Modelling DNA Loops Using the Theory of Elasticity. Phys. Rev. E 2006, 73, 031919. (32) Villa, E.; Balaeff, A.; Mahadevan, L.; Schulten, K. Multiscale Method for Simulating Protein-DNA Complexes. Multiscale Model. Simul. 2004, 2, 527−553. (33) Starostin, E. L. Symmetric Equilibria of a Thin Elastic Rod with Self-Contacts. Philos. Trans. R. Soc., A 2004, 362, 1317−1334. (34) Coleman, B. D.; Swigon, D. Theory of Self-Contact in Kirchhoff Rods with Applications to Supercoiling of Knotted and Unknotted DNA Plasmids. Philos. Trans. R. Soc., A 2004, 362, 1281−1299. (35) Finite Element Procedures; Bathe, K.-J., Ed.; Prentice Hall Inc.: Upper Saddle River, New Jersey, NJ, 1996. (36) Castro, C. E.; Kilchherr, F.; Kim, D.-N.; Shiao, E. L.; Wauer, T.; Wortmann, P.; Bathe, M.; Dietz, H. A Primer to Scaffolded DNA Origami. Nat. Methods 2011, 8, 221−229. (37) Pan, K.; Kim, D.-N.; Zhang, F.; Adendorff, M. R.; Yan, H.; Bathe, M. Lattice-Free Prediction of Three-Dimensional Structure of Programmed DNA Assemblies. Nat. Commun. 2014, 5, 5578. (38) Kim, D.-N.; Kilchherr, F.; Dietz, H.; Bathe, M. Quantitative Prediction of 3D Solution Shape and Flexibility of Nucleic Acid Nanostructures. Nucleic Acids Res. 2012, 40, 2862−2868. (39) Sedeh, R. S.; Pan, K.; Adendorff, M. R.; Hallatschek, O.; Bathe, K.-J.; Bathe, M. Computing Nonequilibrium Conformational Dynamics of Structured Nucleic Acid Assemblies. J. Chem. Theory Comput. 2016, 12, 261−273. (40) Understanding Molecular Simulation: From Algorithms to Applications; Frenkel, D., Smit, B., Eds.; Academic Press: San Diego, CA, 1996. (41) Hayward, S.; de Groot, B. L. Normal Modes and Essential Dynamics. Methods Mol. Biol. 2008, 443, 89−106. (42) Elasticity, Theory, Applications and Numerics; Sadd, M. H., Ed.; Elsevier: Amsterdam, 2009. (43) Foloppe, N.; MacKerell, A. D. All-Atom Empirical Force Field for Nucleic Acids: I. Parameter Optimization Based on Small Molecule and Condensed Phase Macromolecular Target Data. J. Comput. Chem. 2000, 21, 86−104. (44) MacKerell, A. D.; Banavali, N. K. All-Atom Empirical Force Field for Nucleic Acids: II. Application to Molecular Dynamics Simulations of DNA and RNA in Solution. J. Comput. Chem. 2000, 21, 105−120. (45) MATLAB Primer, R2015b; The Math Works Inc.: 3 Apple Hill Drive Natick, MA 01760-2098, 2015. (46) Dhakal, S.; Adendorff, M. R.; Liu, M.; Yan, H.; Bathe, M.; Walter, N. G. Rational Design of DNA-actuated Enzyme Nanoreactors Guided by Single Molecule Analysis. Nanoscale 2016, 8, 3125−3137. (47) Moré, J. J.; Sorensen, D. C. Computing a Trust Region Step. SIAM J. Sci. Stat. Comp. 1983, 4 (3), 553−572. (48) Byrd, R. H.; Schnabel, R. B.; Shultz, G. A. Approximate Solution of the Trust Region Problem by Minimization Over Two-Dimensional Subspaces. Math. Program. 1988, 40, 247−263. O


Article

Journal of Chemical Information and Modeling (49) Horn, B. K. P. Closed-Form Solution of Absolute Orientation Using Unit Quaternions. J. Opt. Soc. Am. A 1987, 4 (4), 629−642. (50) The results are obtained implementing the BFGS method, the trust region method where the Hessian matrix is computed using the proposed finite element formulation, and the trust region method where the Hessian matrix stored is a banded form (11 main diagonals) of that computed using the proposed finite element formulation (the potential energy values are in kcal/mol, while the forces are in kcal/(A mol) = 4.1868E25 pN/mol). (51) Stretch moduli are derived from static equilibrium by assuming external forces at the DNA backbone of the top base pair and fixing the z-direction DOFs of all the bottom base pair atoms. The x,ydirections for the atoms connected to the backbone of the bottom base pairs were also fixed. The first number of the stretch modulus is derived by assuming as initial length the difference of the mean values of the z-coordinate of the base pairs and as displacement the mean value of the z-direction displacement of the top base pair (with the bottom base pair being fixed in this direction). The second number in parentheses is the stretch modulus derived assuming as initial length the mean value of the z-coordinate of the distance between the atoms connecting the backbone for top and bottom base pairs, e.g. the atoms which are fixed and the atoms where the external forces are applied and as displacement the mean value of the z-direction displacement of the two atoms where the external forces are applied.

P


Generalized Potential Energy Finite Elements for Modeling Molecular

Recommend Documents