Access the full text.
Sign up today, get DeepDyve free for 14 days.
Loading next page...
/lp/wiley/continuum-modeling-with-functional-lennard-jones-parameters-for-dna-mwRIcgbvSB
- Publisher
- Wiley
- Copyright
- © 2023 Wiley‐VCH GmbH
- eISSN
- 2513-0390
- DOI
- 10.1002/adts.202200896
- Publisher site
- See Article on Publisher Site
Abstract
IntroductionDue to their high surface area, carbon materials have been explored as platforms for molecular adsorption for uses in a variety of applications. Gas storage,[1–4] pollution capture,[5–7] biomedical uses,[8–10] gas sensing,[11–13] and drug encapsulation[14–16] are some of the applications that make use of the non‐covalent interactions between carbon nanostructures and stored materials. Since carbon compounds are typically nonpolar, the non‐covalent interactions between carbon nanostructures and other molecules are dominated by van der Waals forces. To model these interactions and obtain the potential energy for the system, the 6‐12 Lennard‐Jones potential is commonly used due to its simple form and computational efficiency.[17]The 6‐12 Lennard‐Jones potential for two atoms at a distance ρ apart is given by1Φ(ρ)=−Aρ6+Bρ12$$\begin{equation} \hspace*{80pt}\mathrm{\Phi}(\rho )=-\frac{A}{{\rho}^{6}}+\frac{B}{{\rho}^{12}} \end{equation}$$where A and B are the attractive and repulsive constants, respectively. As shown in Equation (1), the Lennard‐Jones potential is an interatomic potential, so for intermolecular interactions, one method to obtain the total potential energy is to sum over the pairwise interaction between atoms i and j of the two interacting molecules, namely2E=∑i∑jΦ(ρij)$$\begin{equation} \hspace*{80pt}E=\sum _{i}\sum _{j}\mathrm{\Phi}({\rho}_{\textit{ij}}) \end{equation}$$This method is referred to as a discrete approach. Another method to obtain the total interaction energy, which we refer to as a continuum approach, is to replace the above double summations with double surface integrals,3E=η1η2∫S1∫S2Φ(ρ)dS1dS2$$\begin{equation} E={\eta}_{1}{\eta}_{2}\int_{{S}_{1}}\int_{{S}_{2}}\mathrm{\Phi}(\rho )\text{d}{S}_{1}\text{d}{S}_{2} \end{equation}$$where ρ here is a distance between typical surface elements dS1$\text{d}{S}_{1}$ and dS2$\text{d}{S}_{2}$, η1 and η2 are the surface densities of atoms on the two molecules. This approach was first proposed by Girifalco[18] to determine the potential energy between interacting fullerenes. Over the past decade, extensive development of this approach has been made to determine interaction energies for a variety of molecular structures with various shapes, sizes, and chemical compositions.[19–22] The key advantage of a continuum model lies in the resultant analytical expression of the potential energy, which can be readily used to determine characteristics of an interaction, such as critical distance, energy minima, and forces at an arbitrary distance and configuration. However, the continuum model is not straightforwardly used for irregular shaped structures.In the standard continuum approach mentioned above, the coefficients (A and B) in the Lennard‐Jones potential are constants, which is suitable for modeling interactions between mononuclear molecules since they can be thought of as homogeneous surfaces with uniform atomic distribution.[23,24] Examples include interactions between two carbon nanostructures, such as fullerene–fullerene, fullerene‐nanotube, fullerene‐nanocone, nanotube–nanotube, and graphene–graphene.[25] For interactions involving heterogeneous structures or heteronuclear molecules, the standard continuum approach can be adapted by either averaging contributions of the constituent atoms as a homogeneous approximation,[26,27] or dividing the molecules into multiple surfaces of the same atomic type.[28,29] Even though these extensions are relatively simple to implement, the homogeneous approximation produces results that poorly agree with molecular dynamics simulations, and the multi‐surface method requires multiple integral computations which can be cumbersome for large molecules. Recently, a new development for the continuum approach has been introduced to capture heterogeneous structure where the constants A and B are replaced by functional coefficients known as interaction functions, A(ρ)$A(\rho )$ and B(ρ)$B(\rho )$.[30] These functions represent the change in the interaction profile depending on the location of different types of atoms across the heterogeneous molecule. Although excellent agreements to molecular dynamics simulations are obtained, the integration involved becomes more complicated than having A and B as constants. Importantly, care must be taken to determine the interaction functions so that the integration is tractable.We comment that the method of averaging contribution from all constituent atoms is the simplest to upscale for studying large heterogeneous molecules, whereas the multi‐surface method is not as straightforward due to the summation of all interactions from all the surfaces. As such, the multi‐surface method can get computationally intensive for large heterogeneous structures. The method with the interaction functions is relatively easy to upscale compared to the multi‐surface method as there is only one surface to handle.In refs. [30–32], we use the continuum approach with the interaction functions to model the interactions of carbon structures with coronene and methane, which are small heterogeneous structures comprising only two types of atoms. In this paper, we investigate potential use of interaction functions to model larger heterogeneous structures comprising more than two types of atoms.In particular, we consider a heterogeneous molecule with clear patterns or repeated structures due to the organization of their atoms, such as a periodic structure in polymers. One example of such structures is a single‐stranded DNA (ssDNA) molecule, and so we model its interaction with a graphene sheet. Based on the patterns within a ssDNA molecule, we propose a technique where we combine the method of averaging over regions of homogeneous surfaces and using the interaction functions to describe the interaction profile across the regions within a heterogeneous molecule.We note that we particularly choose to study the interaction between a ssDNA molecule and a graphene sheet in order to compare our results to that of Alshehri et al.[26] who adopt a continuum approach with homogenous approximation, assuming atoms on the ssDNA molecule are uniformly smeared across DNA's helical surface.We also note a number of proposed applications involving DNA‐graphene interactions, including molecular logic gates,[33] material synthesis,[34] and electronic biosensing.[35] Thus, understanding adsorption mechanisms of DNA onto graphene sheet is crucial for the development of such applications.In the following section, we describe the newly combined approach taken to obtain analytical expressions for the interaction energy. Note that full details for evaluating integral expressions are given in Appendices A and B, while detail for molecular dynamics simulations used to benchmark analytical solution is given in Appendix C. In Section 3, we investigate the effects of the parameters on the interaction energy, as well as discuss improvements that can be made when implementing this modeling technique to other molecules in the future.Modeling ApproachIn this section, we model the interaction between a single‐stranded B‐DNA molecule and a graphene sheet (see Figure 1). We assume that a graphene sheet is lying flat on the xy-$\textit{xy}\text{-}$plane and that the DNA is a helix of length c and radius b. Following Alshehri et al.,[26] a general point on the surface of the helix has coordinates4cν2π,btsinν,btcosν+δ$$\begin{align} \left(\frac{c\nu}{2\pi},\textit{bt}\sin\nu ,\textit{bt}\cos\nu +\delta \right) \end{align}$$where ν∈[0,2π]$\nu \in [0,2\pi ]$, t∈[0,1]$t\in [0,1]$ are surface parameters and δ is the perpendicular distance between the xy-$\textit{xy}\text{-}$plane and the helix axis. While Alshehri et al.[26] average the contribution of atoms on the DNA to obtain constants A and B, this paper considers the atomic structure of the DNA and thus, we use the Lennard‐Jones potential with functional interaction coefficients to model the DNA. Here, the potential energy between a ssDNA molecule and a graphene sheet is evaluated by the expression5E=ηg−K3+K6$$\begin{equation} E={\eta}_{\text{g}}\left(-{K}_{3}+{K}_{6}\right) \end{equation}$$where ηg${\eta}_{\text{g}}$ is the atomic density of the graphene sheet, and Kn${K}_{n}$ is determined by the computation of the integral6Kn=η1∫SDNA∫Sgfn(ρ)ρ−2ndSDNAdSg$$\begin{equation} {K}_{n}={\eta}_{1}\int_{{S}_{\text{DNA}}}\int_{{S}_{\text{g}}}{f}_{n}(\rho ){\rho}^{-2n}\text{d}{S}_{\text{DNA}}\text{d}{S}_{\text{g}} \end{equation}$$where η1 is the atomic density of the ssDNA molecule, ρ is the distance between a point on the surface of the ssDNA and the graphene sheet, dSDNA$\text{d}{S}_{\text{DNA}}$ and dSg$\text{d}{S}_{\text{g}}$ are the surface elements of the DNA helix and graphene sheet respectively, and fn${f}_{n}$ for n=3,6$n=3,6$ are the interaction functions replacing the standard constants A and B. The physical constants of DNA and graphene used in this model are found in Table 1.1TableNumerical values of parameters used in calculationsParameterValuessDNA radius, b10.2Å$10.2\, \AA $region split point, τ0.6ssDNA length, c34Å$34\, \AA $base pair frequency, m111backbone frequency, m211ssDNA density, η10.7069Å−2$0.7069\hspace*{3.33333pt}{\mathrm{\AA}}^{-2}$base pair region density, η1,i${\eta}_{1,i}$0.6285Å−2$0.6285\hspace*{3.33333pt}{\mathrm{\AA}}^{-2}$backbone region density, η1,o${\eta}_{1,o}$0.8262Å−2$0.8262\hspace*{3.33333pt}{\mathrm{\AA}}^{-2}$graphene sheet density, ηg${\eta}_{\text{g}}$0.3812Å−2$0.3812\hspace*{3.33333pt}{\mathrm{\AA}}^{-2}$1FigureAtomic representation of a ssDNA molecule above a graphene sheet and its corresponding diagram representing the ssDNA as a helical surface and the graphene sheet as a flat plane.Since a graphene sheet is a homogeneous surface of carbon atoms, the choice of interaction function only depends on the structure of the ssDNA molecule. The ssDNA molecule consists of two main regimes along the length of the strand, namely the phosphate‐deoxyribose backbone, comprising alternating molecules of phosphate and deoxyribose in a periodic fashion, and the inner region comprising the four (sometimes five) types of nucleobases (see Figure 2a). A single turn of ssDNA typically has either 10 or 11 nucleobases, often written as 10.5.[36] This allows us to assume the interaction function that is periodic in ν with a frequency of 11. While there is a variability in the distribution of the atoms throughout the DNA as shown in Figure 2a, the regularity in the distribution of the constituent molecular groups is clear, as shown in Figure 2b. Thus, it seems intuitive to model each molecular group on the DNA as a homogeneous region. So, for example, rather than attempting to capture the heterogeneity within a phosphate group, we treat the entire phosphate as a homogeneous surface with average contribution of phosphorous and oxygen. Additionally, since we do not consider a particular DNA, specific sequence of nucleobases is not relevant in the model. As such, we assume that all four nucleobases have an equal likelihood to be on the DNA and thus we approximate an averaged nucleobase as a homogeneous surface over which its constituent atoms are smeared and the contributions of all four bases are averaged.2FigureVisual representations of a ssDNA molecule. a) A stick and ball model ( nitrogen; oxygen; carbon; hydrogen; phosphorous). b) A point plot with groups of molecular structures represented by colors ( deoxyribose; phosphate; nucleobase).Based on these assumptions, there are three possible candidates for the interaction function: a) a piecewise interaction function that varies in ν (i.e., fn,1(ν)${f}_{n,1}(\nu )$ for t∈[0,τ]$t\in [0,\tau ]$ and fn,2(ν)${f}_{n,2}(\nu )$ for t∈[τ,1]$t\in [\tau ,1]$); b) a piecewise function that varies in t (i.e., fn,1(t)${f}_{n,1}(t)$ for ν∈[(2i)2π11,(2i+1)2π11]$\nu \in [\frac{(2i)2\pi}{11},\frac{(2i+1)2\pi}{11}]$ and fn,2(t)${f}_{n,2}(t)$ for ν∈[(2i+1)2π11,(2i+2)2π11]$\nu \in [\frac{(2i+1)2\pi}{11},\frac{(2i+2)2\pi}{11}]$ for i=0,1,⋯,5$i=0,1,\text{\ensuremath{\cdots}},5$, not including the final interval as there is an odd number of strips); and finally c) an interaction function that varies continuously in both variables fn(ν,t)${f}_{n}(\nu ,t)$. Visually, (a) and (b) correspond to the images in Figure 3, where the first is two continuous surfaces separated along t=τ$t=\tau $, the second image has 11 continuous surfaces stacked vertically along the helix, and (c) would be a fully continuous combination of (a) and (b). This division into regions also naturally leads to using individual densities (ηi${\eta}_{i}$) for each region rather than a single density for the whole DNA strand.3FigurePossible configurations of how a ssDNA molecule is divided into regions a) a series of stacked horizontal strips and b) two vertically aligned strips.According to (a), we consider interaction functions of the form αcos(mν)+β$\alpha \cos(m\nu )+\beta $ in order to capture the alternating nature of the two regimes, phosphate and deoxyribose in the outer band and nucleobases and empty space in the inner band. The resulting integral under this choice of the interaction function is given by7Kn=η1,i∫0τ∫02π∫Sgα1cosm1ν+β1ρ2ndSDNAdSg︸inner band+η1,o∫τ1∫02π∫Sgα2cosm2ν+β2ρ2ndSDNAdSg︸outer band,$$\begin{eqnarray} {K}_{n} &=& \underset{\text{inner band}}{\underbrace{{\eta}_{1,i}\int_{0}^{\tau}\int_{0}^{2\pi}\int_{{S}_{\text{g}}}\frac{{\alpha}_{1}\cos\left({m}_{1}\nu \right)+{\beta}_{1}}{{\rho}^{2n}}\text{d}{S}_{\text{DNA}}\text{d}{S}_{\text{g}}}} \nonumber\\ &&+\,\underset{\text{outer band}}{\underbrace{{\eta}_{1,o}\int_{\tau}^{1}\int_{0}^{2\pi}\int_{{S}_{\text{g}}}\frac{{\alpha}_{2}\cos\left({m}_{2}\nu \right)+{\beta}_{2}}{{\rho}^{2n}}\text{d}{S}_{\text{DNA}}\text{d}{S}_{\text{g}}}}, \end{eqnarray}$$where the values of the coefficients αi${\alpha}_{i}$ and βi${\beta}_{i}$ (i=1,2$i=1,2$) are determined from the interaction coefficients (A and B) of the constituent molecules interacting with graphene, and mi${m}_{i}$ is determined by the number of base pairs in one turn.According to (b), we consider interaction functions that are sigmoidal in shape since we need to capture the transition from nucleobase to deoxyribose and empty space to phosphate across the strips. Interaction functions that are sigmoidal in shape have previously been used to capture transition between two types of atoms[30] so we assume similar behavior to hold when transitioning between two regions within a molecule. In particular we use functions of the form αarctan(μ(t−τ))+β$\alpha \arctan(\mu (t-\tau ))+\beta $, which implies that the resulting integral for evaluation is of the form8Kn=η1,1∫01∫2i2π112i+12π11∫Sgα1arctanμ1t−τ+β1ρ2ndSDNAdSg+η1,2∫01∫2i+12π112i+22π11∫Sgα2arctanμ2t−τ+β2ρ2ndSDNAdSg$$\begin{eqnarray} {K}_{n}&=& {\eta}_{1,1}\int_{0}^{1}\!\!\int_{\frac{\left(2i\right)2\pi}{11}}^{\frac{\left(2i+1\right)2\pi}{11}}\!\!\!\int_{{S}_{\text{g}}}\frac{{\alpha}_{1}\arctan\left({\mu}_{1}\left(t-\tau \right)\right)+{\beta}_{1}}{{\rho}^{2n}}\text{d}{S}_{\text{DNA}}\text{d}{S}_{\text{g}} \nonumber\\ && +\,{\eta}_{1,2}\!\int_{0}^{1}\!\!\int_{\frac{\left(2i+1\right)2\pi}{11}}^{\frac{\left(2i+2\right)2\pi}{11}}\!\!\!\int_{{S}_{\text{g}}}\!\frac{{\alpha}_{2}\arctan\left({\mu}_{2}\left(t-\tau \right)\right)+{\beta}_{2}}{{\rho}^{2n}}\text{d}{S}_{\text{DNA}}\text{d}{S}_{\text{g}} \end{eqnarray}$$where the coefficients αi${\alpha}_{i}$ and βi${\beta}_{i}$ are determined in a similar manner to the first choice above, η1, 1 and η1, 2 are atomic densities for each of the horizontal strips, and μi${\mu}_{i}$ is determined by the slope of the transition from the inner to outer region.According to (c), we consider an interaction function that has combined properties from both (a) and (b), that is, we expect a sigmoidal profile when fixing ν and a sinusoidal shape when fixing t. This behavior can be achieved by taking the cosine function above and replacing the coefficients α and β with arctan functions, namely (α1arctan(μ1(t−τ))+β1)cos(mν)+(α2arctan(μ2(t−τ))+β2)$({\alpha}_{1}\arctan({\mu}_{1}(t-\tau ))+{\beta}_{1})\cos(m\nu )+({\alpha}_{2}\arctan({\mu}_{2}(t-\tau ))+{\beta}_{2})$, where the coefficients of the arctan functions are determined from the coefficients of the original cosine functions from the first choice. Note that this only works if both cosine functions have matching frequencies m, otherwise the frequency will need to be modified by some function of t that smoothly transitions from m1 to m2. With this choice of interaction functions, the integral Kn${K}_{n}$ of the form9Kn=η1∫01∫02π∫Sgα1arctanμ1t−τ+β1cosmν+α2arctanμ2t−τ+β2ρ2ndSDNAdSg$$\begin{equation} {K}_{n}={\eta}_{1}\int_{0}^{1}\int_{0}^{2\pi}\int_{{S}_{\text{g}}}\frac{\left({\alpha}_{1}\arctan\left({\mu}_{1}\left(t-\tau \right)\right)+{\beta}_{1}\right)\cos\left(m\nu \right)+\left({\alpha}_{2}\arctan\left({\mu}_{2}\left(t-\tau \right)\right)+{\beta}_{2}\right)}{{\rho}^{2n}}\text{d}{S}_{\text{DNA}}\text{d}{S}_{\text{g}} \end{equation}$$Of all these three choices of interaction functions and the corresponding integrals, we find that only Equation (7) can be evaluated analytically, namely10Kn=bc2n−22mαη1,i∑i=0m∑j=0∞A1J1+4πβη1,iδ2n−2∑i=0∞A2J2+2mαη1,o∑i=0m∑j=0∞A1J3+4πβη1,oδ2n−2∑i=0∞A2J4$$\begin{align} {K}_{n}& =\frac{\textit{bc}}{2n-2}\left\{2m\alpha {\eta}_{1,i}\sum _{i=0}^{m}\sum _{j=0}^{\infty}{A}_{1}{J}_{1}+\frac{4\pi \beta {\eta}_{1,i}}{{\delta}^{2n-2}}\sum _{i=0}^{\infty}{A}_{2}{J}_{2}+2m\alpha {\eta}_{1,o}\sum _{i=0}^{m}\sum _{j=0}^{\infty}{A}_{1}{J}_{3}+\frac{4\pi \beta {\eta}_{1,o}}{{\delta}^{2n-2}}\sum _{i=0}^{\infty}{A}_{2}{J}_{4}\right\} \end{align}$$where11A1=−1i22m−i2m−i−1!B12,m−i+122n−2j12j2bji!j!2m−2i!m−i+1jA2=n−1in−12ib2ii!2δ2i$$\begin{align} {A}_{1}& =\frac{{\left(-1\right)}^{i}{2}^{2\left(m-i\right)}\left(2m-i-1\right)!B\left(\frac{1}{2},m-i+\frac{1}{2}\right){\left(2n-2\right)}_{j}{\left(\frac{1}{2}\right)}_{j}{\left(2b\right)}^{j}}{i!j!\left(2m-2i\right)!{\left(m-i+1\right)}_{j}} \nonumber\\ {A}_{2}& =\frac{{\left(n-1\right)}_{i}{\left(n-\frac{1}{2}\right)}_{i}{b}^{2i}}{{\left(i!\right)}^{2}{\delta}^{2i}} \end{align}$$and12J1=δ−j+2n−2∑k=0∞∑l=0k−12k−4b2π2c2k1+4b2π2c2t012−k−t0k−ll!k−l!τj+2l+1j+2l+12F1j+2n−2,j+2l+1;j+2l+2;−bτδJ2=τ2i+12i+12F1−12,2i+12;2i+32;−4b2π2c2τ2J3=1−τδ+bτj+2n−2∑k=0∞∑l=0k−12k−4b2π2c2k1+4b2π2c2t012−k−t0k−lτj+2ll!k−l!F11,−j−2l,j+2n−2,2;τ−1τ,bτ−1δ+bτJ4=∑k=0∞∑l=0k−12k−4b2π2c2k1+4b2π2c2t012−k−t0k−ll!k−l!1−τ2i+2l+12i+2l+1$$\begin{eqnarray} \def\eqcellsep{&}\begin{array}{l} \displaystyle {J}_{1} ={\delta}^{-\left(j+2n-2\right)}\sum _{k=0}^{\infty}\sum _{l=0}^{k}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}\right)}^{k}{\left(1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}_{0}\right)}^{\frac{1}{2}-k}{\left(-{t}_{0}\right)}^{k-l}}{l!\left(k-l\right)!}\frac{{\tau}^{j+2l+1}}{j+2l+1}{}_{2}{F}_{1}\left(j+2n-2,j+2l+1;j+2l+2;-\frac{b\tau}{\delta}\right)\\[12pt] \displaystyle {J}_{2} =\frac{{\tau}^{2i+1}}{2i+1}{}_{2}{F}_{1}\left(-\frac{1}{2},\frac{2i+1}{2};\frac{2i+3}{2};-\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{\tau}^{2}\right)\\[12pt] \displaystyle {J}_{3} =\frac{\left(1-\tau \right)}{{\left(\delta +b\tau \right)}^{j+2n-2}}\sum _{k=0}^{\infty}\sum _{l=0}^{k}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}\right)}^{k}{\left(1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}_{0}\right)}^{\frac{1}{2}-k}{\left(-{t}_{0}\right)}^{k-l}{\tau}^{j+2l}}{l!\left(k-l\right)!}{F}_{1}\left(1,-j-2l,j+2n-2,2;\frac{\tau -1}{\tau},\frac{b\left(\tau -1\right)}{\delta +b\tau}\right)\\[12pt] \displaystyle {J}_{4} =\sum _{k=0}^{\infty}\sum _{l=0}^{k}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}\right)}^{k}{\left(1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}_{0}\right)}^{\frac{1}{2}-k}{\left(-{t}_{0}\right)}^{k-l}}{l!\left(k-l\right)!}\frac{\left(1-{\tau}^{2i+2l+1}\right)}{2i+2l+1} \end{array} \end{eqnarray}$$where 2F1 denotes the standard hypergeometric function and F1 denotes an Appell's hypergeometric function. Detailed derivation of these expressions is given in Appendix A. For Equation (8), we obtain a semi‐analytical solution where one integration is remained to be evaluated numerically (see Appendix B). Due to the high computational time to compute the semi‐analytical solution numerically, we do not adopt this case here. We also do not progress with Equation (9) since we are only interested in analytical solution for Kn${K}_{n}$. High computational time is also expected for Equation (9).Determination of Coefficients α1, 2 and β1, 2 Used in Equation (10)In order to adopt (10) to determine the potential energy between a DNA and a graphene sheet, we first need to determine the coefficients αi${\alpha}_{i}$ and βi${\beta}_{i}$ (i=1,2$i=1,2$), which are the interaction coefficients of the constituent molecular groups interacting with a graphene sheet. We note that two sets of the interaction coefficients are considered. The first set is derived from the carbon–carbon values in Table 2, while the second set uses the graphene–graphene values.2TableValues of the van der Waals radius (σ) and the energy well depth (ε) for each interaction used to determine the Lennard‐Jones's attractive and repulsive constants: A=4εσ6$A=4\varepsilon {\sigma}^{6}$ and B=4εσ12$B=4\varepsilon {\sigma}^{12}$Interactionσ [Å]ε [kcal mol−1]A [kcal mol−1${}^{-1}\hspace*{3.33333pt}$Å6]B [kcal mol−1${}^{-1}\hspace*{3.33333pt}$Å12]C−Ca)3.4310.105684.9521117047$1\hspace*{3.33333pt}117\hspace*{3.33333pt}047$P−Pa)3.6950.3053102.6597890568$7\hspace*{3.33333pt}890\hspace*{3.33333pt}568$N−Na)3.2610.069331.7154398677.9$398\hspace*{3.33333pt}677.9$O−Oa)3.1180.060220.5919202753.2$202\hspace*{3.33333pt}753.2$H−Ha)2.5710.04450.8464514689.55$14\hspace*{3.33333pt}689.55$Graphene−Grapheneb)3.4150.055350.23555299.54$555\hspace*{3.33333pt}299.54$a)From ref. [37];b)From ref. [38].In refs. [30, 31] methane‐nanotube and coronene‐graphene interactions were considered. In these studies, the heterogeneous molecule (methane/coronene) can be modeled as two regions (inner region of carbon atoms and outer region of hydrogen atoms). The interaction functions are then used to describe the interaction profile of a carbon surface (nanotube/graphene) interacting with the two regions. In refs. [30–32] the sigmoidal functions were used to represent a smooth transition from the inner carbon core to the outer hydrogen region. As shown in refs. [30–32] the coefficients for the interaction functions are relatively simple to determine since they can be restricted to match the Lennard‐Jones constants AC--C${A}_{\text{C--C}}$ (or BC--C${B}_{\text{C--C}}$) at the carbon region and AC--H${A}_{\text{C--H}}$ (or BC--H${B}_{\text{C--H}}$) at the hydrogen region.For our choice of interaction functions presented in Equation (7) for modeling a ssDNA molecule, we determine αi${\alpha}_{i}$ and βi${\beta}_{i}$ (i=1,2$i=1,2$) such that when n=3$n=3$, α1+β1=Anucleobase-x${\alpha}_{1}+{\beta}_{1}={A}_{\text{nucleobase-x}}$ and α2+β2=Adeoxyribose-x${\alpha}_{2}+{\beta}_{2}={A}_{\text{deoxyribose-x}}$ and when n=6$n=6$, α1+β1=Bnucleobase-x${\alpha}_{1}+{\beta}_{1}={B}_{\text{nucleobase-x}}$ and α2+β2=Bdeoxyribose-x${\alpha}_{2}+{\beta}_{2}={B}_{\text{deoxyribose-x}}$. Note that x represents interaction involving either carbon or graphene. Similarly, we prescribe that β1−α1=0${\beta}_{1}-{\alpha}_{1}=0$ (when n=3$n=3$) or 0 (when n=6$n=6$) and β2−α2=Aphosphate-x${\beta}_{2}-{\alpha}_{2}={A}_{\text{phosphate-x}}$ (when n=3$n=3$) or Bphosphate-x${B}_{\text{phosphate-x}}$ (when n=6$n=6$).We derive the Lennard‐Jones constants of each region by employing the homogeneous smearing method using the parameters given in Table 2. Using the standard combination rules for Lennard‐Jones parameters, we generate two sets of interatomic interaction constants Ai-C,Bi-C${A}_{\text{i-C}},{B}_{\text{i-C}}$, and Ai-graphene,Bi-graphene${A}_{\text{i-graphene}},{B}_{\text{i-graphene}}$, where i is any other type of molecular groups. For example, the attractive constant for carbon and phosphate (PO4$\text{P}{\text{O}}_{4}$) is calculated from Aphosphate-C=15AP-C+45AO-C${A}_{\text{phosphate-C}}=\frac{1}{5}{A}_{\text{P-C}}+\frac{4}{5}{A}_{\text{O-C}}$. Using the values from Table 2, we derive the interaction coefficients for the constituent molecules of DNA interacting with either carbon or graphene as given in Table 3.3TableInteraction coefficients for the constituent molecules of DNA interacting with either carbon or grapheneMoleculeA★−C${A}_{\star -C}$B★−C${B}_{\star -C}$A★−G${A}_{\star -G}$B★−G${B}_{\star -G}$Deoxyribose400.49544857.74$544\hspace*{3.33333pt}857.74$360.70499546.24$499\hspace*{3.33333pt}546.24$Phosphate605.86984648.81$984\hspace*{3.33333pt}648.81$439.97673728.00$673\hspace*{3.33333pt}728.00$Nucleobase456.18647484.77$647\hspace*{3.33333pt}484.77$390.01565623.69$565\hspace*{3.33333pt}623.69$Using the values shown in Table 3, we obtain α1, 2 and β1, 2 for both sets of coefficients as presented in Table 4. In the following section, we plot the energy obtained from Equation (10) using these two set of coefficients. The results are compared with molecular dynamics simulations and the homogeneous model where atoms in the DNA are smeared over the helical surface.4TableNumerical values for α1, 2 and β1, 2 for the interaction function (7) derived from the values in Table 3. Note the set of coefficients (i) and (ii) are based on the interaction involving carbon–carbon and graphene–graphene, respectivelySetnα1β1α2β2(i)3228.09228.09−102.69503.186323742.39$323\hspace*{3.33333pt}742.39$323742.39$323\hspace*{3.33333pt}742.39$−219895.54$-219\hspace*{3.33333pt}895.54$764753.28$764\hspace*{3.33333pt}753.28$(ii)3195.01195.01−39.64400.33448996282811.85$282\hspace*{3.33333pt}811.85$282811.85$282\hspace*{3.33333pt}811.85$−87090.88$-87\hspace*{3.33333pt}090.88$586637.12$586\hspace*{3.33333pt}637.12$Results and DiscussionIn this section, the interaction between a ssDNA and a flat graphene sheet is considered. Particular attention is made to the DNA structure in order to use a continuum approach to study such interaction. In Alshehri et al.,[26] the DNA is modeled as a homogeneous helical surface where the effect of constituent atoms is averaged throughout the entire structure. As shown in Figure 4, the potential energy obtained from the homogeneous model does not agree well with the results from molecular dynamics (MD) simulations (see Appendix C for details of the MD study). We compare our new results from the expression derived from Equation (10) using the two sets of coefficients found in Table 4, denoted as models (i) and (ii) respectively, with the results from the homogeneous model described in Alshehri et al.[26] and the numerical results from the MD study. For the interaction coefficients used in the homogeneous model, Alshehri et al.[26] considered a specific DNA sequence, 5′‐CCACTAGTGG‐3′, so the attractive and repulsive constants are obtained by taking the arithmetic mean of all the individual atomic interaction coefficients (e.g., ADNA−C=1320(97AC−C+10AP−C+114AH−C+61AO−C+38AN−C)${A}_{\text{DNA}-\text{C}}=\frac{1}{320}(97{A}_{\text{C}-\text{C}}+10{A}_{\text{P}-\text{C}}+114{A}_{\text{H}-\text{C}}+61{A}_{\text{O}-\text{C}}+38{A}_{\text{N}-\text{C}})$). We note that the total number of atoms in a DNA is determined from the number of atoms in each molecular group in its bonded state, meaning that there are less hydrogen atoms than the standalone molecules. In general, for every base pair there is a corresponding deoxyribose molecule (C5OH7${\text{C}}_{5}{\text{OH}}_{7}$) and a phosphate molecule (PO4${\text{PO}}_{4}$). Of the base pairs, there are guanine (C5H4N5O${\text{C}}_{5}{\text{H}}_{4}{\text{N}}_{5}\text{O}$), adenine (C5H4N5${\text{C}}_{5}{\text{H}}_{4}{\text{N}}_{5}$), cytosine (C4H4N3O${\text{C}}_{4}{\text{H}}_{4}{\text{N}}_{3}\text{O}$), and thymine (C5H5N2O2${\text{C}}_{5}{\text{H}}_{5}{\text{N}}_{2}{\text{O}}_{2}$). Last there are the lone hydrogen and hydroxide at the 5′ and 3′ end, respectively. The interaction coefficients we use for the homogeneous model in this paper are derived in the same manner, though rather than considering a specific sequence of nucleobases, eleven averaged nucleobases are used instead. Specifically, we use C194H174N154O${\text{C}}_{\frac{19}{4}}{\text{H}}_{\frac{17}{4}}{\text{N}}_{\frac{15}{4}}\text{O}$ as the chemical composition of an averaged nucleobase when determining the interaction coefficients.4FigureEnergy profiles comparing models (i) and (ii) with the homogeneous model and the MD results ( molecular dynamics simulations; model (i); model (ii); and − homogeneous model).It can be seen that model (i) performs even more poorly compared to the homogeneous model. For model (ii), the result is closer to the MD simulations. This indicates that using coefficients based on graphene–graphene interaction is more appropriate than using coefficients based on C–C interaction for modeling interactions involving graphene sheet.Next, we explore the impact of other parameters in the model on the interaction energy profile. We first investigate the two shape parameters of the helix, namely the radius b and the length c. From Figure 4, we find the location of the minimum energy (δmin${\delta}_{\text{min}}$) and the corresponding minimum energy (Emin${E}_{\text{min}}$) as shown in Table 5. The values of δmin${\delta}_{\text{min}}$ from the three models are found to be in good agreement with the MD results, despite varying widely in the values of Emin${E}_{\text{min}}$.5TableEnergy minima and critical distance of the models and numerical resultsModelδmin${\delta}_{\text{min}}$ [Å]Emin${E}_{\text{min}}$ [kcal mol−1]MD12.939−17.600Model (i)12.881−26.495Model (ii)12.864−21.496Homogeneous12.846−23.544We show in Figure 5 that the location δmin${\delta}_{\text{min}}$ is strongly affected by the radius of the ssDNA (b). In Figure 5, we demonstrate that varying the values of b by 0.2 Å results in the same amount of displacement in δmin${\delta}_{\text{min}}$, as well as a slight change in Emin${E}_{\text{min}}$ due to helical surface being further away from the graphene sheet.5FigureEnergy profiles of model (ii) for varying the radius of the helix, b, where 10.2Å$10.2\, \AA $ is the value used in other calculations ( b=10$b= 10$ Å; b=10.2$b=10.2$ Å; −b=10.4$b= 10.4$ Å).On the other hand, it is shown in Figure 6 that an increase in the length of the helix, c, over 10 Å is required to produce a decrease in Emin${E}_{\text{min}}$ of a similar magnitude to that produced by increasing the radius, b by 0.2 Å.6FigureEnergy profiles of model (ii) with varied lengths of the helix, c. A value of 34 Å is used for c in the other calculations ( c=32$c= 32$ Å; c=34$c= 34$ Å; −c=50$c= 50$ Å).These results demonstrate the significance of DNA's radius over its length in the model. However, these findings indicate that the size of the helix largely impacts the location of the minimum energy (δmin${\delta}_{\text{min}}$), but only minimally affects the value of minimum energy (Emin${E}_{\text{min}}$).To further investigate the cause of the discrepancies in the energy profiles shown in Figure 4, we consider interaction coefficients, A and B, used for interactions involving phosphate, deoxyribose, and nucleobase. We comment that the coefficients A and B calculated from averaging the contribution of constituent atoms within a molecule do not provide a result that is in good agreement with the MD simulations. Due to the non‐linear effect that the distance has on the Lennard‐Jones potential, the configuration of the atoms within the molecule needs to be accounted for. In Figure 7, two homogeneous sphere models for phosphate interacting with graphene are compared against simulation results. One has coefficients determined from averaging the contribution of each constituent atoms (i.e., A=15APC+45AOC$A=\frac{1}{5}{A}_{\text{PC}}+\frac{4}{5}{A}_{\text{OC}}$) and the other has coefficients determined by fitting a standard sphere‐plane continuum model[25] with the externally obtained well depth (ε) and the van der Waals radius (σ) (i.e., Emin=−ε${E}_{\text{min}}=-\varepsilon $ and E′(216σ)=0${E}^{\prime}({2}^{\frac{1}{6}}\sigma )=0$). From the energy profiles, it can be seen that the average atomic contribution approach leads to interaction coefficients Aphosphate--carbon${A}_{\text{phosphate--carbon}}$ and Bphosphate--carbon${B}_{\text{phosphate--carbon}}$ that do not produce result that is in good agreement to simulation results. This is contrasted with the coefficients obtained from fitting a sphere‐plane model, which matches the simulation results more closely. Even though there is a big difference between the two approaches, the fitted coefficients approach comes with some drawbacks. These include the needs for a means of accounting for the molecule's shape for a homogeneous model and for information about the well depth and the location of the well depth for the interaction considered. Thus, this method would only be suitable for determining coefficients of smaller molecular units as part of a larger molecule. As a result, we do not pursue the fitted coefficients approach for DNA molecule studied here as it requires further assumptions for the shapes of nucleobases and deoxyribose and the determination of ε and σ for these molecular groups and graphene.7FigureEnergy profiles for phosphate group interacting with graphene (, molecular dynamics simulations; model with coefficients deriving from averaging contribution of constituent atoms; model with coefficients fitted from a sphere‐plane model).Another possibility for the discrepancy between the continuum model and the simulation results is the assumption made to the structure of the DNA molecule. Figure 8a shows how the helical‐surface model, with the right parameters, captures the overall shape of the ssDNA molecule but not the volume. As the nucleobases are orthogonal to the backbone and the backbone itself is not flat, a model comprising two helical cylinders capturing nucleobases and backbone atoms may be more suitable for modeling a ssDNA strand. An example of this modeling approach is shown in Figure 8b.8FigureOne turn of a ssDNA molecule's atomic coordinates plotted along with a) a helical‐surface and b) two helical cylinders.Concluding RemarksThe Lennard‐Jones potential is widely used for determining interaction energy between non‐bonded molecules. This potential together with a continuum approach are adopted when dealing with regular shaped molecules since the resultant energy can be found explicitly as a function of distance between two interacting molecules. Challenges arise when the interaction involves irregular shaped structures or heterogeneous molecules, which comprise more than one types of atoms. For heterogeneous molecules, the assumptions underlying traditional continuum approach that atoms are evenly smeared throughout the molecule and the effects of atoms are the same regardless of its location on the molecule seem invalid. To address this issue, we propose a new continuum approach for a heterogeneous structure, where the Lennard‐Jones coefficients are replaced by interaction functions.[30–32] While this approach has been shown to work well with methane and coronene, we note their relatively simple structures, which can be modeled as a sphere and a disk, respectively. Furthermore, their atomic arrangement of carbon atoms at the central region and hydrogen atoms at the perimeter enables the interaction functions in the form of a sigmoidal profile, which makes the integrals in the interaction energy traceable yielding analytical outcomes. For larger and more complex heterogeneous molecules, there has yet to be a continuum method that takes into account both atomic arrangement and structure of a molecule. As such, this paper can be viewed as an introduction to a continuum approach using the Lennard‐Jones potential with interaction functions to model complicated heterogeneous molecules, such as polymers and DNAs.A continuum model for heterogeneous molecules is proposed here for the modeling of single‐stranded DNA interacting with a graphene sheet. By recognizing the patterns of atoms on the DNA molecule, we are able to develop a technique combining the method of averaging atomic contributions over regions of surfaces and using the interaction functions to describe the interaction profile across the regions within a heterogeneous molecule. This new approach appears to improve the accuracy of the model when benchmarking with molecular dynamics simulations. Despite some improvement shown here, the model can be further benefited by a systematic method for determining the interaction coefficients between the constituent molecules (namely phosphate, deoxyribose, and the nucleobases) and the graphene sheet. Different interaction functions may also be explored to better capture the heterogeneity of the DNA molecule. Additionally, the model can also be further improved by using a volume‐based model to better capture all atoms of the DNA.These suggested improvements to the continuum modeling for DNAs provide areas for future research directions. Finally, we note that the approach presented in this paper can be extended to model polymers or other heterogeneous molecules with repeated patterns and structures.AAppendixComputation of the Integral Kn${K}_{n}$ in Equation (7)A ssDNA helix has coordinates (cν2π,btsinν,btcosν+δ)$(\frac{c\nu}{2\pi},\textit{bt}\sin\nu ,\textit{bt}\cos\nu +\delta )$, where c is the length of the DNA strand, b is the radius of the DNA strand, δ is the distance from the helix axis to graphene sheet. Here, we assume that δ>b$\delta >b$ since the DNA molecule cannot intersect the graphene sheet. For a graphene sheet lying flat on the xy-$\textit{xy}\text{-}$plane, its coordinates are given by (x,y,0)$(x,y,0)$. Typical surface elements for graphene and DNA are given respectively by dSg=dxdy$\text{d}{S}_{\text{g}}=\text{d}x\text{d}y$, dSDNA=bb2t2+c24π2dνdt$\text{d}{S}_{\text{DNA}}=b\sqrt{{b}^{2}{t}^{2}+\frac{{c}^{2}}{4{\pi}^{2}}}\text{d}\nu \text{d}t$, where x,y∈(−∞,∞)$x,y\in (-\infty ,\infty )$, t∈[0,1]$t\in [0,1]$ and ν∈[0,2π]$\nu \in [0,2\pi ]$.From the coordinates of graphene and the DNA, distance between two typical surface element is given by ρ2=(cν2π−x)2+(btsinν−y)2+(btcosν+δ)2${\rho}^{2}={(\frac{c\nu}{2\pi}-x)}^{2}+{(\textit{bt}\sin\nu -y)}^{2}+{(\textit{bt}\cos\nu +\delta )}^{2}$. Thus, the integral Kn${K}_{n}$ over both surfaces can be written asA1Kn=η1b∫01∫02π∫−∞∞∫−∞∞fnν,tρ2nb2t2+c24π2dxdydνdt=η1bc2π∫01∫02πfnν,t1+4b2π2c2t2∫−∞∞∫−∞∞ρ−2ndxdydνdt$$\begin{eqnarray} {K}_{n} &=& {\eta}_{1}b\int_{0}^{1}\int_{0}^{2\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{{f}_{n}\left(\nu ,t\right)}{{\rho}^{2n}}\sqrt{{b}^{2}{t}^{2}+\frac{{c}^{2}}{4{\pi}^{2}}}\text{d}x\text{d}y\text{d}\nu \text{d}t \nonumber\\ &=& {\eta}_{1}\frac{\textit{bc}}{2\pi}\int_{0}^{1}\int_{0}^{2\pi}{f}_{n}\left(\nu ,t\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}{\rho}^{-2n}\text{d}x\text{d}y\text{d}\nu \text{d}t \end{eqnarray}$$Next, we write ρ2=Λ12+(cν2π−x)2${\rho}^{2}={\mathrm{\Lambda}}_{1}^{2}+{(\frac{c\nu}{2\pi}-x)}^{2}$, and make the substitution x−cν2π=Λ1tanϕ$x-\frac{c\nu}{2\pi}={\mathrm{\Lambda}}_{1}\tan\phi $. This changes the limits of the integration from (−∞,∞)$(-\infty ,\infty )$ to [−π2,π2]$[-\frac{\pi}{2},\frac{\pi}{2}]$ and the differential dx=Λ1sec2ϕdϕ$\text{d}x={\mathrm{\Lambda}}_{1}{\operatorname{sec}}^{2}\phi \text{d}\phi $. With the substitution the integral Kn${K}_{n}$ becomesA2Kn=η1bc2π∫01∫02πfnν,t1+4b2π2c2t2∫−∞∞∫−π2π2Λ1sec2ϕΛ12nsec2nϕdϕdydνdt=η1bc2π∫01∫02πfnν,t1+4b2π2c2t2×∫−∞∞1Λ12n−1∫−π2π2cos2n−2ϕdϕdydνdt=η1bc22n−3!!2n−2!!∫01∫02πfnν,t1+4b2π2c2t2∫−∞∞1Λ12n−1dydνdt$$\begin{eqnarray} {K}_{n} &=& \frac{{\eta}_{1}\textit{bc}}{2\pi}\!\int_{0}^{1}\!\!\int_{0}^{2\pi}{f}_{n}\left(\nu ,t\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\int_{-\infty}^{\infty}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\frac{{\mathrm{\Lambda}}_{1}{\operatorname{sec}}^{2}\phi}{{\mathrm{\Lambda}}_{1}^{2n}{\operatorname{sec}}^{2n}\phi}\text{d}\phi \text{d}y\text{d}\nu \text{d}t \nonumber\\ &=& \frac{{\eta}_{1}\textit{bc}}{2\pi}\int_{0}^{1}\int_{0}^{2\pi}{f}_{n}\left(\nu ,t\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}} \nonumber\\ &&\times\,\int_{-\infty}^{\infty}\frac{1}{{\mathrm{\Lambda}}_{1}^{2n-1}}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}{\cos}^{2n-2}\phi \text{d}\phi \text{d}y\text{d}\nu \text{d}t \nonumber\\ &=& \frac{{\eta}_{1}\textit{bc}}{2}\frac{\left(2n-3\right)!!}{\left(2n-2\right)!!}\int_{0}^{1}\int_{0}^{2\pi}{f}_{n}\left(\nu ,t\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\int_{-\infty}^{\infty}\frac{1}{{\mathrm{\Lambda}}_{1}^{2n-1}}\text{d}y\text{d}\nu \text{d}t \nonumber\\ \end{eqnarray}$$Note that ∫0π2cos2nxdx=(2n−1)!!(2m)!!π2$\int_{0}^{\frac{\pi}{2}}{\cos}^{2n}x\text{d}x=\frac{(2n-1)!!}{(2m)!!}\frac{\pi}{2}$ and n!! is the double factorial. Now, we again write Λ12=Λ22+(y−btsinν)2${\mathrm{\Lambda}}_{1}^{2}={\mathrm{\Lambda}}_{2}^{2}+{(y-\textit{bt}\sin\nu )}^{2}$ and substitute y−btsinν=Λ2tanϕ$y-\textit{bt}\sin\nu ={\mathrm{\Lambda}}_{2}\tan\phi $, soA3Kn=η1bc22n−3!!2n−2!!∫01∫02πfnν,t1+4b2π2c2t2×∫−π2π2Λ2sec2ϕΛ22sec2ϕn−12dϕdνdt=η1bc22n−3!!2n−2!!∫01∫02πfnν,tΛ22n−21+4b2π2c2t2∫−π2π2cos2n−3ϕdϕdνdt=η1bc2n−2∫01∫02πfnν,tbtcosν+δ2n−21+4b2π2c2t2dνdt$$\begin{eqnarray} {K}_{n} &=& \frac{{\eta}_{1}\textit{bc}}{2}\frac{\left(2n-3\right)!!}{\left(2n-2\right)!!}\int_{0}^{1}\int_{0}^{2\pi}{f}_{n}\left(\nu ,t\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}} \nonumber\\ &&\times\,\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\frac{{\mathrm{\Lambda}}_{2}{\operatorname{sec}}^{2}\phi}{{\left({\mathrm{\Lambda}}_{2}^{2}{\operatorname{sec}}^{2}\phi \right)}^{n-\frac{1}{2}}}\text{d}\phi \text{d}\nu \text{d}t \nonumber\\ &=& \frac{{\eta}_{1}\textit{bc}}{2}\frac{\left(2n\,{-}\,3\right)!!}{\left(2n\,{-}\,2\right)!!}\int_{0}^{1}\!\!\int_{0}^{2\pi}\!\frac{{f}_{n}\left(\nu ,t\right)}{{\mathrm{\Lambda}}_{2}^{2n-2}}\sqrt{1{+}\,\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\!\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\!{\cos}^{2n-3}\phi \text{d}\phi \text{d}\nu \text{d}t \nonumber\\ &=& \frac{{\eta}_{1}\textit{bc}}{2n-2}\int_{0}^{1}\int_{0}^{2\pi}\frac{{f}_{n}\left(\nu ,t\right)}{{\left(\textit{bt}\cos\nu +\delta \right)}^{2n-2}}\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\text{d}\nu\text{d}t \end{eqnarray}$$since ∫0π2cos2n+1xdx=(2n)!!(2m+1)!!$\int_{0}^{\frac{\pi}{2}}{\cos}^{2n+1}x\text{d}x=\frac{(2n)!!}{(2m+1)!!}$The rest of the analytical calculation depends on the choice of fn(ν,t)${f}_{n}(\nu ,t)$. Dividing the DNA strand into two regimes along the t-$t\text{-}$axis, we have the nucleobases in the interval [0, τ] and the phosphate‐deoxyribose backbone in the interval [τ, 1]. Accordingly, fn${f}_{n}$ is set to be two different functions of ν in these intervals, namely fn,1=α1cos(m1ν)+β1${f}_{n,1}={\alpha}_{1}\cos({m}_{1}\nu )+{\beta}_{1}$ in [0, τ] and fn,2=α2cos(m2ν)+β2${f}_{n,2}={\alpha}_{2}\cos({m}_{2}\nu )+{\beta}_{2}$ in [τ, 1]. This also splits Kn${K}_{n}$ into two double integrals,A4Kn=bc2n−2η1,i∫0τ1+4b2π2c2t2∫02πα1cosm1ν+β1btcosν+δ2n−2dνdt+η1,o∫τ11+4b2π2c2t2∫02πα2cosm2ν+β2btcosν+δ2n−2dνdt$$\begin{eqnarray} {K}_{n} &=& \frac{\textit{bc}}{2n-2}\left\{\eta_{1,i}\int_{0}^{\tau}\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\int_{0}^{2\pi}\frac{{\alpha}_{1}\cos\left({m}_{1}\nu \right)+{\beta}_{1}}{{\left(\textit{bt}\cos\nu +\delta \right)}^{2n-2}}\text{d}\nu \text{d}t \right. \nonumber\\[6pt] &&+\left.{\eta}_{1,o}\int_{\tau}^{1}\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\int_{0}^{2\pi}\frac{{\alpha}_{2}\cos\left({m}_{2}\nu \right)+{\beta}_{2}}{{\left(\textit{bt}\cos\nu +\delta \right)}^{2n-2}}\text{d}\nu \text{d}t\right\} \end{eqnarray}$$The integral in ν is effectively the same for both regions, so we only need to solve the integral of the formA5I=∫02παcosmν+βbtcosν+δ2n−2dν=α∫02πcosmνbtcosν+δ2n−2dν+β∫02π1btcosν+δ2n−2dν=αI1+βI2$$\begin{eqnarray} I &=&\int_{0}^{2\pi}\frac{\alpha \cos\left(m\nu \right)+\beta}{{\left(\textit{bt}\cos\nu +\delta \right)}^{2n-2}}\text{d}\nu \nonumber\\ &=&\alpha \int_{0}^{2\pi}\frac{\cos\left(m\nu \right)}{{\left(\textit{bt}\cos\nu +\delta \right)}^{2n-2}}\text{d}\nu +\beta \int_{0}^{2\pi}\frac{1}{{\left(\textit{bt}\cos\nu +\delta \right)}^{2n-2}}\text{d}\nu \nonumber\\ &=& \alpha {I}_{1}+\beta {I}_{2} \end{eqnarray}$$After applying a double‐angle trigonometric identity, the integral I2 can be readily computed as a hypergeometric function using a special case shown in equation 3.681.1 in ref. [39],A6∫0π2dx1−k2sin2xϱ=π22F1ϱ,12;1;k2$$\begin{equation} \int_{0}^{\frac{\pi}{2}}\frac{\text{d}x}{{\left(1-{k}^{2}{\sin}^{2}x\right)}^{\varrho}}=\frac{\pi}{2}{}_{2}{F}_{1}\left(\varrho ,\frac{1}{2};1;{k}^{2}\right) \end{equation}$$which we then expand into an infinite sum in order to facilitate the integration with respect to t:A7I2=2∫0π1δ+bt−2btsin2ν22n−2dν=4δ+bt2n−2∫0π211−2btδ+btsin2ν2n−2dν=2πδ+bt2n−22F12n−2,12;1;2btδ+bt$$\begin{eqnarray} {I}_{2}&=&2\int_{0}^{\pi}\frac{1}{{\left(\delta +\textit{bt}-2\textit{bt}{\sin}^{2}\frac{\nu}{2}\right)}^{2n-2}}\text{d}\nu \nonumber\\ &=&\frac{4}{{\left(\delta +\textit{bt}\right)}^{2n-2}}\int_{0}^{\frac{\pi}{2}}\frac{1}{{\left(1-\frac{2\textit{bt}}{\delta +\textit{bt}}{\sin}^{2}\nu \right)}^{2n-2}}\text{d}\nu \nonumber\\ &=& \frac{2\pi}{{\left(\delta +\textit{bt}\right)}^{2n-2}}{}_{2}{F}_{1}\left(2n-2,\frac{1}{2};1;\frac{2\textit{bt}}{\delta +\textit{bt}}\right) \end{eqnarray}$$Using a quadratic transformation of the 2F1 hypergeometric function,A82F1a,b;2b;z=1−z2−a2F1a2,a+12;b+12;z2−z2$$\begin{equation} {}_{2}{F}_{1}\left(a,b;2b;z\right)={\left(1-\frac{z}{2}\right)}^{-a}{}_{2}{F}_{1}\left(\frac{a}{2},\frac{a+1}{2};b+\frac{1}{2};{\left(\frac{z}{2-z}\right)}^{2}\right) \end{equation}$$both the coefficient and the argument of the hypergeometric function in I2 are simplified toA9I2=2πδ+bt2n−21−btδ+bt2n−22F12n−22,2n−12;1;2btδ+bt2−2btδ+bt2=2πδ2n−22F1n−1,n−12;1;btδ2=2πδ2n−2∑i=0∞n−1in−12ii!2btδ2i=2πδ2n−2∑i=0∞A2t2i$$\begin{eqnarray} {I}_{2}& =&\frac{2\pi}{{\left(\delta +\textit{bt}\right)}^{2n-2}}{\left(1-\frac{\textit{bt}}{\delta +\textit{bt}}\right)}^{2n-2}{}_{2}{F}_{1}\left(\frac{2n-2}{2},\frac{2n-1}{2};1;{\left(\frac{\frac{2\textit{bt}}{\delta +\textit{bt}}}{2-\frac{2\textit{bt}}{\delta +\textit{bt}}}\right)}^{2}\right) \nonumber\\ & =&\frac{2\pi}{{\delta}^{2n-2}}{}_{2}{F}_{1}\left(n-1,n-\frac{1}{2};1;{\left(\frac{\textit{bt}}{\delta}\right)}^{2}\right) \nonumber\\ & =&\frac{2\pi}{{\delta}^{2n-2}}\sum _{i=0}^{\infty}\frac{{\left(n-1\right)}_{i}{\left(n-\frac{1}{2}\right)}_{i}}{{\left(i!\right)}^{2}}{\left(\frac{\textit{bt}}{\delta}\right)}^{2i} \nonumber\\ & =&\frac{2\pi}{{\delta}^{2n-2}}\sum _{i=0}^{\infty}{A}_{2}{t}^{2i} \end{eqnarray}$$The integral I1 is more interesting due to the cosine term. By making the assumption that m is an integer renders it tractable as we can make use of De Moivre's formula:A10cos2mx=Rcos2mx+isin2mx=Rcosx+isinx2m=R∑i=02m2micos2m−ixiisinix=∑i=0m2m2icos2n−2ixi2isin2ix=∑i=0m−1i2m2icos2m−ixsin2ix$$\begin{align} \cos\left(2\textit{mx}\right)& =\mathfrak{R}\left(\cos\left(2\textit{mx}\right)+\mathrm{i}\sin\left(2\textit{mx}\right)\right)\nonumber\\ & =\mathfrak{R}\left({\left(\cos x+\mathrm{i}\sin x\right)}^{2m}\right)\nonumber\\ & =\mathfrak{R}\left(\sum _{i=0}^{2m}\left(\genfrac{}{}{0pt}{}{2m}{i}\right){\cos}^{2m-i}\left(x\right){\mathrm{i}}^{i}{\sin}^{i}\left(x\right)\right)\nonumber\\ & =\sum _{i=0}^{m}\left(\genfrac{}{}{0pt}{}{2m}{2i}\right){\cos}^{2n-2i}\left(x\right){\mathrm{i}}^{2i}{\sin}^{2i}\left(x\right)\nonumber\\ & =\sum _{i=0}^{m}{\left(-1\right)}^{i}\left(\genfrac{}{}{0pt}{}{2m}{2i}\right){\cos}^{2\left(m-i\right)}x{\sin}^{2i}x \end{align}$$where i is the imaginary unit. Thus, we haveA11I1=2∫0πcosmνbtcosν+δ2n−2dν=4∫0π2cos2mνδ+bt−2btsin2ν2n−2dν=4δ+bt2n−2∫0π2∑i=0m−1i2m2i×cos2m−iνsin2iν11−2btδ+btsin2ν2n−2dν=4δ+bt2n−2∑i=0m−1i2m2i∫0π2cos2m−iνsin2iν1−2btδ+btsin2ν2n−2dν=42m!δ+bt2n−2∑i=0m−1i2i!2m−2i!∫0π2cos2m−iνsin2iν1−2btδ+btsin2ν2n−2dν=42m!δ+bt2n−2∑i=0m−1i2i!2m−2i!12Bi+12,m−i+12×2F12n−2,i+12;m+1;2btδ+bt=22m!∑i=0m∑j=0∞−1iBi+12,m−i+122n−2ji+12j2bj2i!j!2m−2i!m+1j×tjδ+btj+2n−2=22m!∑i=0m∑j=0∞A1tjδ+btj+2n−2$$\begin{align} {I}_{1}& =2\int_{0}^{\pi}\frac{\cos\left(m\nu \right)}{{\left(\textit{bt}\cos\nu +\delta \right)}^{2n-2}}\text{d}\nu \nonumber\\ & =4\int_{0}^{\frac{\pi}{2}}\frac{\cos\left(2m\nu \right)}{{\left(\delta +\textit{bt}-2\textit{bt}{\sin}^{2}\nu \right)}^{2n-2}}\text{d}\nu \nonumber\\ & =\frac{4}{{\left(\delta +\textit{bt}\right)}^{2n-2}}\int_{0}^{\frac{\pi}{2}}\sum _{i=0}^{m}{\left(-1\right)}^{i}\left(\genfrac{}{}{0pt}{}{2m}{2i}\right) \nonumber\\ &\ \ \ \times \,{\cos}^{2\left(m-i\right)}\nu {\sin}^{2i}\nu \frac{1}{{\left(1-\frac{2\textit{bt}}{\delta +\textit{bt}}{\sin}^{2}\nu \right)}^{2n-2}}\text{d}\nu \nonumber\\ & =\frac{4}{{\left(\delta +\textit{bt}\right)}^{2n-2}}\sum _{i=0}^{m}{\left(-1\right)}^{i}\left(\genfrac{}{}{0pt}{}{2m}{2i}\right)\int_{0}^{\frac{\pi}{2}}\frac{{\cos}^{2\left(m-i\right)}\nu {\sin}^{2i}\nu}{{\left(1-\frac{2\textit{bt}}{\delta +\textit{bt}}{\sin}^{2}\nu \right)}^{2n-2}}\text{d}\nu \nonumber\\ & =\frac{4\left(2m\right)!}{{\left(\delta +\textit{bt}\right)}^{2n-2}}\sum _{i=0}^{m}\frac{{\left(-1\right)}^{i}}{\left(2i\right)!\left(2m-2i\right)!}\int_{0}^{\frac{\pi}{2}}\frac{{\cos}^{2\left(m-i\right)}\nu {\sin}^{2i}\nu}{{\left(1-\frac{2\textit{bt}}{\delta +\textit{bt}}{\sin}^{2}\nu \right)}^{2n-2}}\text{d}\nu \nonumber\\ & =\frac{4\left(2m\right)!}{{\left(\delta +\textit{bt}\right)}^{2n-2}}\sum _{i=0}^{m}\frac{{\left(-1\right)}^{i}}{\left(2i\right)!\left(2m-2i\right)!}\frac{1}{2}B\left(i+\frac{1}{2},m-i+\frac{1}{2}\right) \nonumber\\ &\ \ \ \times \,{}_{2}{F}_{1}\left(2n-2,i+\frac{1}{2};m+1;\frac{2\textit{bt}}{\delta +\textit{bt}}\right)\nonumber\\ & =2\left(2m\right)!\sum _{i=0}^{m}\sum _{j=0}^{\infty}\frac{{\left(-1\right)}^{i}B\left(i+\frac{1}{2},m-i+\frac{1}{2}\right){\left(2n-2\right)}_{j}{\left(i+\frac{1}{2}\right)}_{j}{\left(2b\right)}^{j}}{\left(2i\right)!j!\left(2m-2i\right)!{\left(m+1\right)}_{j}} \nonumber\\ &\ \ \ \times \,\frac{{t}^{j}}{{\left(\delta +\textit{bt}\right)}^{j+2n-2}}\nonumber\\ & =2\left(2m\right)!\sum _{i=0}^{m}\sum _{j=0}^{\infty}{A}_{1}\frac{{t}^{j}}{{\left(\delta +\textit{bt}\right)}^{j+2n-2}} \end{align}$$where B(x,y)$B(x,y)$ is the beta function that comes from equation 3.681.1 in ref. [39]A12∫0π2sin2μ−1xcos2ν−1x1−k2sin2xϱdx=12Bμ,ν2F1ϱ,μ;μ+ν;k2$$\begin{equation} \int_{0}^{\frac{\pi}{2}}\frac{{\sin}^{2\mu -1}x{\cos}^{2\nu -1}x}{{\left(1-{k}^{2}{\sin}^{2}x\right)}^{\varrho}}\text{d}x=\frac{1}{2}B\left(\mu ,\nu \right){}_{2}{F}_{1}\left(\varrho ,\mu ;\mu +\nu ;{k}^{2}\right) \end{equation}$$Combining the results for I1 and I2, the integral I can be written as a function of t with parameters α, β, and m, namelyA13It;α,β,m=2mα∑i=0m∑j=0∞A1tjδ+btj+2n−2+4πβδ2n−2∑i=0∞A2t2i$$\begin{align} I\left(t;\alpha ,\beta ,m\right)& =2m\alpha \sum _{i=0}^{m}\sum _{j=0}^{\infty}{A}_{1}\frac{{t}^{j}}{{\left(\delta +\textit{bt}\right)}^{j+2n-2}}+\frac{4\pi \beta}{{\delta}^{2n-2}}\sum _{i=0}^{\infty}{A}_{2}{t}^{2i} \end{align}$$Thus, the expression for Kn${K}_{n}$ is given byA14Kn=bc2n−2η1,i∫0τIt;α1,β1,m11+4b2π2c2t2dt+η1,o∫τ1It;α2,β2,m21+4b2π2c2t2dt$$\begin{eqnarray} {K}_{n}&=&\frac{\textit{bc}}{2n-2}\left\{{\eta}_{1,i}\int_{0}^{\tau}I\left(t;{\alpha}_{1},{\beta}_{1},{m}_{1}\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\text{d}t \right. \nonumber\\ &&+\left.{\eta}_{1,o}\int_{\tau}^{1}I\left(t;{\alpha}_{2},{\beta}_{2},{m}_{2}\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\text{d}t\right\} \end{eqnarray}$$which means only four integrals remain, namelyA15J1=∫0τtjδ+btj+2n−21+4b2π2c2t2dt,J2=∫0τt2i1+4b2π2c2t2dt,J3=∫τ11δ+btj+2n−21+4b2π2c2t2dtJ4=∫τ1t2i1+4b2π2c2t2dt$$\begin{align} {J}_{1}& =\int_{0}^{\tau}\frac{{t}^{j}}{{\left(\delta +\textit{bt}\right)}^{j+2n-2}}\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\text{d}t, \nonumber\\ {J}_{2}& =\int_{0}^{\tau}{t}^{2i}\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\text{d}t,\nonumber\\ {J}_{3}& =\int_{\tau}^{1}\frac{1}{{\left(\delta +\textit{bt}\right)}^{j+2n-2}}\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\text{d}t\nonumber\\ {J}_{4}& =\int_{\tau}^{1}{t}^{2i}\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\text{d}t \end{align}$$In some of the above integrals, namely J1 and J3, we approximate 1+4b2π2c2t2$\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}$ by its Taylor series about a point t0 to render them more tractable. The general series is given byA161±γtr=∑i=0∞−ri∓γi1±γt0r−ii!t−t0i$$\begin{equation} {\left(1\pm \gamma t\right)}^{r}=\sum _{i=0}^{\infty}\frac{{\left(-r\right)}_{i}{\left(\mp \gamma \right)}^{i}{\left(1\pm \gamma {t}_{0}\right)}^{r-i}}{i!}{\left(t-{t}_{0}\right)}^{i} \end{equation}$$so, for γ=4b2π2c2$\gamma =\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}$ and r=12$r=\frac{1}{2}$, we haveA171+4b2π2c2t2=∑k=0∞−12k−4b2π2c2k1+4b2π2c2t012−kk!t2−t0k$$\begin{eqnarray} \sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}=\sum _{k=0}^{\infty}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}\right)}^{k}{\left(1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}_{0}\right)}^{\frac{1}{2}-k}}{k!}{\left({t}^{2}-{t}_{0}\right)}^{k} \nonumber\\ \end{eqnarray}$$As a general rule, we select t0 to be the midpoint of the interval over which we are integrating. This means t0=τ2${t}_{0}=\frac{\tau}{2}$ for J1 and t0=1−τ2${t}_{0}=\frac{1-\tau}{2}$ for J3.We first evaluate J2 since it is the most straightforward integral compared to the others. From equation 3.254.1 in ref. [39] we have ∫0uxλ−1(u−x)μ−1(x2+β2)νdx=β2νuλ+μ−1B(λ,μ)3F2(−ν,λ2,λ+12;λ+μ2,λ+μ+12;−u2β2)$\int_{0}^{u}{x}^{\lambda -1}{(u-x)}^{\mu -1}{({x}^{2}+{\beta}^{2})}^{\nu}\text{d}x={\beta}^{2\nu}{u}^{\lambda +\mu -1}B(\lambda ,\mu ){}_{3}{F}_{2}(-\nu ,\frac{\lambda}{2},\frac{\lambda +1}{2};\frac{\lambda +\mu}{2},\frac{\lambda +\mu +1}{2};\frac{-{u}^{2}}{{\beta}^{2}})$, where pFq${}_{p}{F}_{q}$ is a generalized hypergeometric function given by pFq(a1,⋯,ap;b1,⋯,bq;z)=∑i=0∞(a1)i⋯(ap)i(b1)i⋯(bq)izii!${}_{p}{F}_{q}({a}_{1},\text{\ensuremath{\cdots}},{a}_{p};{b}_{1},\text{\ensuremath{\cdots}},{b}_{q};z)={\sum}_{i=0}^{\infty}\frac{{({a}_{1})}_{i}\text{\ensuremath{\cdots}}{({a}_{p})}_{i}}{{({b}_{1})}_{i}\text{\ensuremath{\cdots}}{({b}_{q})}_{i}}\frac{{z}^{i}}{i!}$.Thus,A18J2=∫0τtλ1+γ2t2dt=γ∫0τtλ1γ2+t2dt=γ1γτλ+1Bλ+1,13F2−12,λ+12,λ+22;λ+22,λ+32;−γ2τ2=τλ+1λ+12F1−12,λ+12;λ+32;−γ2τ2$$\begin{align} {J}_{2}& =\int_{0}^{\tau}{t}^{\lambda}\sqrt{1+{\gamma}^{2}{t}^{2}}\text{d}t\nonumber\\ & =\gamma \int_{0}^{\tau}{t}^{\lambda}\sqrt{\frac{1}{{\gamma}^{2}}+{t}^{2}}\text{d}t\nonumber\\ & =\gamma \frac{1}{\gamma}{\tau}^{\lambda +1}B\left(\lambda +1,1\right){}_{3}{F}_{2}\left(-\frac{1}{2},\frac{\lambda +1}{2},\frac{\lambda +2}{2};\frac{\lambda +2}{2},\frac{\lambda +3}{2};-{\gamma}^{2}{\tau}^{2}\right)\nonumber\\ & =\frac{{\tau}^{\lambda +1}}{\lambda +1}{}_{2}{F}_{1}\left(-\frac{1}{2},\frac{\lambda +1}{2};\frac{\lambda +3}{2};-{\gamma}^{2}{\tau}^{2}\right) \end{align}$$where the reduction in order of the 3F2 to a 2F1 occurs as the λ+22$\frac{\lambda +2}{2}$ terms cancel, as implied by the series expansion of the general pFq${}_{p}\textit{Fq}$.For the evaluation of J4, we use the substitution γt=sinhθ$\gamma t=\operatorname{\mathrm{\sinh}}\theta $, the recurrence relation ∫sinhpθcoshqθdθ=sinhp−1θcoshq−1θp+q+q−1p+q∫sinhpθcoshq−2θdθ$\int {\operatorname{\mathrm{\sinh}}}^{p}\theta {\cosh}^{q}\theta \text{d}\theta =\frac{{\operatorname{\mathrm{\sinh}}}^{p-1}\theta {\cosh}^{q-1}\theta}{p+q}+\frac{q-1}{p+q}\int {\operatorname{\mathrm{\sinh}}}^{p}\theta {\cosh}^{q-2}\theta \text{d}\theta $ and the antiderivative ∫sinh2nθdθ=(−1)n(2nn)θ22n+122n−1∑i=0n−1(−1)n(2ni)sinh((2n−2i)θ)2n−2i$\int {\operatorname{\mathrm{\sinh}}}^{2n}\theta \text{d}\theta ={(-1)}^{n}(\genfrac{}{}{0pt}{}{2n}{n})\frac{\theta}{{2}^{2n}}+\frac{1}{{2}^{2n-1}}{\sum}_{i=0}^{n-1}{(-1)}^{n}(\genfrac{}{}{0pt}{}{2n}{i})\frac{\operatorname{\mathrm{\sinh}}((2n-2i)\theta )}{2n-2i}$. As such,A19J4=∫τ1tλ1+γ2t2dt=1γλ+1∫t=τt=1sinhλθcosh2θdθ=1γλ+1sinhλ−1θcoshθλ+2+1λ+2∫sinhλθdθt=τt=1=1γλ+1sinhλ−1θcoshθλ+2+1λ+2−1n2nnθ22n∑i=0n−1+122n−1∑i=0n−1−1n2nisinh2n−2iθ2n−2iθ=arcsinhγτθ=arcsinhγ$$\begin{align} {J}_{4}& =\int_{\tau}^{1}{t}^{\lambda}\sqrt{1+{\gamma}^{2}{t}^{2}}\text{d}t\nonumber\\ & =\frac{1}{{\gamma}^{\lambda +1}}\int_{t=\tau}^{t=1}{\operatorname{\mathrm{\sinh}}}^{\lambda}\theta {\cosh}^{2}\theta \text{d}\theta \nonumber\\ & =\frac{1}{{\gamma}^{\lambda +1}}{\left.\left(\frac{{\operatorname{\mathrm{\sinh}}}^{\lambda -1}\theta \cosh\theta}{\lambda +2}+\frac{1}{\lambda +2}\int {\operatorname{\mathrm{\sinh}}}^{\lambda}\theta \text{d}\theta \right)\right|}_{t=\tau}^{t=1}\nonumber\\ & =\frac{1}{{\gamma}^{\lambda +1}}\left.\left(\frac{{\operatorname{\mathrm{\sinh}}}^{\lambda -1}\theta \cosh\theta}{\lambda +2}+\frac{1}{\lambda +2}\left({\left(-1\right)}^{n}\left(\genfrac{}{}{0pt}{}{2n}{n}\right)\frac{\theta}{{2}^{2n}}\vphantom{\sum _{i=0}^{n-1}} \right.\right.\right.\nonumber\\ &\ \ \left.\left.\left.+\,\frac{1}{{2}^{2n-1}}\sum _{i=0}^{n-1}{\left(-1\right)}^{n}\left(\genfrac{}{}{0pt}{}{2n}{i}\right)\frac{\operatorname{\mathrm{\sinh}}\left(\left(2n-2i\right)\theta \right)}{2n-2i}\right)\right)\right|_{\theta =\text{arcsinh}\gamma \tau}^{\theta =\text{arcsinh}\gamma} \end{align}$$Next, we consider J1. From equation 3.194.1 in ref. [39] we have ∫0uxa−1(1+bx)νdx=uaa2F1(ν,a;1+a;−bu)$\int_{0}^{u}\frac{{x}^{a-1}}{{(1+\textit{bx})}^{\nu}}\text{d}x=\frac{{u}^{a}}{a}{}_{2}{F}_{1}(\nu ,a;1+a;-\textit{bu})$, so thatA20J1=∫0τtλδ+btσ1+γ2t2dt=∫0τtλδ+btσ∑k=0∞−12k−γ2k1+γ2t012−kk!t2−t0kdt=∑k=0∞−12k−γ2k1+γ2t012−kk!×∫0τtλδ+btσ∑l=0kk!l!k−l!−t0k−lt2ldt=δ−σ∑k=0∞∑l=0k−12k−γ2k1+γ2t012−k−t0k−ll!k−l!∫0τtλ+2l1+bδtσdt=δ−σ∑k=0∞∑l=0k−12k−γ2k1+γ2t012−k−t0k−ll!k−l!×τλ+2l+1λ+2l+12F1σ,λ+2l+1;λ+2l+2;−bτδ$$\begin{align} {J}_{1}& =\int_{0}^{\tau}\frac{{t}^{\lambda}}{{\left(\delta +\textit{bt}\right)}^{\sigma}}\sqrt{1+{\gamma}^{2}{t}^{2}}\text{d}t\nonumber\\ & =\int_{0}^{\tau}\frac{{t}^{\lambda}}{{\left(\delta +\textit{bt}\right)}^{\sigma}}\sum _{k=0}^{\infty}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-{\gamma}^{2}\right)}^{k}{\left(1+{\gamma}^{2}{t}_{0}\right)}^{\frac{1}{2}-k}}{k!}{\left({t}^{2}-{t}_{0}\right)}^{k}\text{d}t\nonumber\\ & =\sum _{k=0}^{\infty}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-{\gamma}^{2}\right)}^{k}{\left(1+{\gamma}^{2}{t}_{0}\right)}^{\frac{1}{2}-k}}{k!} \nonumber\\ &\ \ \ \times \,\int_{0}^{\tau}\frac{{t}^{\lambda}}{{\left(\delta +\textit{bt}\right)}^{\sigma}}\sum _{l=0}^{k}\frac{k!}{l!\left(k-l\right)!}{\left(-{t}_{0}\right)}^{k-l}{t}^{2l}\text{d}t\nonumber\\ & ={\delta}^{-\sigma}\sum _{k=0}^{\infty}\sum _{l=0}^{k}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-{\gamma}^{2}\right)}^{k}{\left(1+{\gamma}^{2}{t}_{0}\right)}^{\frac{1}{2}-k}{\left(-{t}_{0}\right)}^{k-l}}{l!\left(k-l\right)!}\int_{0}^{\tau}\frac{{t}^{\lambda +2l}}{{\left(1+\frac{b}{\delta}t\right)}^{\sigma}}\text{d}t\nonumber\\ & ={\delta}^{-\sigma}\sum _{k=0}^{\infty}\sum _{l=0}^{k}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-{\gamma}^{2}\right)}^{k}{\left(1+{\gamma}^{2}{t}_{0}\right)}^{\frac{1}{2}-k}{\left(-{t}_{0}\right)}^{k-l}}{l!\left(k-l\right)!} \nonumber\\ &\ \ \ \times \,\frac{{\tau}^{\lambda +2l+1}}{\lambda +2l+1}{}_{2}{F}_{1}\left(\sigma ,\lambda +2l+1;\lambda +2l+2;-\frac{b\tau}{\delta}\right) \end{align}$$Last, we evaluate J3 by using a linear substitution, t=(1−τ)u+τ$t=(1-\tau )u+\tau $, and rewrite the integral in the form of an Appell's F1 function, that is, ∫01xλ−1(1−x)μ−1(1−ux)−ϱ(1−vx)−σdx=B(μ,λ)F1(λ,ϱ,σ,λ+μ;u,v)$\int_{0}^{1}{x}^{\lambda -1}{(1-x)}^{\mu -1}{(1-\textit{ux})}^{-\varrho}{(1-\textit{vx})}^{-\sigma}\textit{dx}=B(\mu ,\lambda ){F}_{1}(\lambda ,\varrho ,\sigma ,\lambda +\mu ;u,v)$. Appell's F1 function is a generalized 2F1 function with two variables instead of one. Thus,A21J3=∫τ1tλδ+btσ1+γ2t2dt=∑k=0∞∑l=0k−12k−γ2k1+γ2t012−k−t0k−ll!k−l!∫τ1tλ+2lδ+btσdt=1−τ∑k=0∞∑l=0k−12k−γ2k1+γ2t012−k−t0k−ll!k−l!×∫011−τu+τλ+2lδ+bτ+b1−τuσdu=1−τ∑k=0∞∑l=0k−12k−γ2k1+γ2t012−k−t0k−ll!k−l!τλ+2lδ+bτσ×∫011+1−ττuλ+2l1+b1−τδ+bτuσdu=1−τδ+bτσ∑k=0∞∑l=0k−12k−γ2k1+γ2t012−k−t0k−lτλ+2ll!k−l!×F11,−λ−2l,σ,2;τ−1τ,bτ−1δ+bτ.$$\begin{align} {J}_{3}& =\int_{\tau}^{1}\frac{{t}^{\lambda}}{{\left(\delta +\textit{bt}\right)}^{\sigma}}\sqrt{1+{\gamma}^{2}{t}^{2}}\text{d}t\nonumber\\ & =\sum _{k=0}^{\infty}\sum _{l=0}^{k}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-{\gamma}^{2}\right)}^{k}{\left(1+{\gamma}^{2}{t}_{0}\right)}^{\frac{1}{2}-k}{\left(-{t}_{0}\right)}^{k-l}}{l!\left(k-l\right)!}\int_{\tau}^{1}\frac{{t}^{\lambda +2l}}{{\left(\delta +\textit{bt}\right)}^{\sigma}}\text{d}t\nonumber\\ & =\left(1-\tau \right)\sum _{k=0}^{\infty}\sum _{l=0}^{k}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-{\gamma}^{2}\right)}^{k}{\left(1+{\gamma}^{2}{t}_{0}\right)}^{\frac{1}{2}-k}{\left(-{t}_{0}\right)}^{k-l}}{l!\left(k-l\right)!} \nonumber\\ &\ \ \ \times \,\int_{0}^{1}\frac{{\left(\left(1-\tau \right)u+\tau \right)}^{\lambda +2l}}{{\left(\delta +b\tau +b\left(1-\tau \right)u\right)}^{\sigma}}\text{d}u\nonumber\\ & =\left(1-\tau \right)\sum _{k=0}^{\infty}\sum _{l=0}^{k}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-{\gamma}^{2}\right)}^{k}{\left(1+{\gamma}^{2}{t}_{0}\right)}^{\frac{1}{2}-k}{\left(-{t}_{0}\right)}^{k-l}}{l!\left(k-l\right)!}\frac{{\tau}^{\lambda +2l}}{{\left(\delta +b\tau \right)}^{\sigma}} \nonumber\\ &\ \ \ \times \,\int_{0}^{1}\frac{{\left(1+\frac{\left(1-\tau \right)}{\tau}u\right)}^{\lambda +2l}}{{\left(1+\frac{b\left(1-\tau \right)}{\delta +b\tau}u\right)}^{\sigma}}\text{d}u\nonumber\\ & =\frac{\left(1-\tau \right)}{{\left(\delta +b\tau \right)}^{\sigma}}\sum _{k=0}^{\infty}\sum _{l=0}^{k}\frac{{\left(-\frac{1}{2}\right)}_{k}{\left(-{\gamma}^{2}\right)}^{k}{\left(1+{\gamma}^{2}{t}_{0}\right)}^{\frac{1}{2}-k}{\left(-{t}_{0}\right)}^{k-l}{\tau}^{\lambda +2l}}{l!\left(k-l\right)!} \nonumber\\ &\ \ \ \times \,{F}_{1}\left(1,-\lambda -2l,\sigma ,2;\frac{\tau -1}{\tau},\frac{b\left(\tau -1\right)}{\delta +b\tau}\right). \end{align}$$Using the results of J1, J2, J3, and J4, we obtain analytical expression for Kn${K}_{n}$ as presented in Equation (10).BAppendixEvaluation of integral Kn${K}_{n}$ Given in Equation (8)For a semi‐analytical solution to the integral Kn${K}_{n}$ in Equation (8), we start from Equation (A3) in the vertical strip derivation, then we break the integral in ν into m intervals, where m is the number of nucleobases in the turn of DNA being modeled. For these horizontal strips, two interaction functions, Tn,1${T}_{n,1}$ and Tn,2${T}_{n,2}$, are used to represent the transition from empty space to phosphate and from nucleobase to deoxyribose, respectively. As the interaction alternates between these regions along the ν interval, we capture this alternating behavior by using Tn,1${T}_{n,1}$ on [(2i)2πm,(2i+1)2πm]$[\frac{(2i)2\pi}{m},\frac{(2i+1)2\pi}{m}]$ and Tn,2${T}_{n,2}$ on [(2j+1)2πm,(2j+2)2πm]$[\frac{(2j+1)2\pi}{m},\frac{(2j+2)2\pi}{m}]$ where i∈{0,1,⋯,⌊m−12⌋}$i\in \{0,1,\text{\ensuremath{\cdots}},\lfloor \frac{m-1}{2}\rfloor \}$ and j∈{0,1,⋯,⌊m−22⌋}$j\in \{0,1,\text{\ensuremath{\cdots}},\lfloor \frac{m-2}{2}\rfloor \}$. Both Tn${T}_{n}$ functions have the same form of αarctan(m(t−t0))+β$\alpha \arctan(m(t-{t}_{0}))+\beta $ to approximate the smooth transition from inner strip to outer strip of the ssDNA molecule.As the integral is no longer invariant to rotations, we must also rotate the helical surface about its own axis by some angle ϕ, thus a strip of the helical surface is parametrized by (cν2π,bt(cosϕsinν+sinϕcosν),bt(−sinϕsinν+cosϕcosν)+δ)$(\frac{c\nu}{2\pi},\textit{bt}(\cos\phi \sin\nu +\sin\phi \cos\nu ),\textit{bt}(-\sin\phi \sin\nu +\cos\phi \cos\nu )+\delta )$. In the x-$x\text{-}$coordinate, we add either a 2icm$\frac{2\textit{ic}}{m}$ or a (2j+1)cm$\frac{(2j+1)c}{m}$ term, depending on which interval the strip segment is in.Combining the above terms, the integral Kn${K}_{n}$ becomesB1Kn=bc2n−2η1,1∑i=0⌊m−12⌋∫01∫2i2πm2i+12πmTn,1t1+4b2π2c2t2Λ22n−2dνdt+η1,2∑i=0⌊m−22⌋∫01∫2i+12πm2i+22πmTn,2t1+4b2π2c2t2Λ22n−2dνdt,$$\begin{align} {K}_{n}& =\frac{\textit{bc}}{2n-2}\left({\eta}_{1,1}\sum _{i=0}^{\lfloor \frac{m-1}{2}\rfloor}\int_{0}^{1}\int_{\frac{\left(2i\right)2\pi}{m}}^{\frac{\left(2i+1\right)2\pi}{m}}\frac{{T}_{n,1}\left(t\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}}{{\mathrm{\Lambda}}_{2}^{2n-2}}\text{d}\nu \text{d}t \right. \nonumber\\ &\ \ \ \ +\left.{\eta}_{1,2}\sum _{i=0}^{\lfloor \frac{m-2}{2}\rfloor}\int_{0}^{1}\int_{\frac{\left(2i+1\right)2\pi}{m}}^{\frac{\left(2i+2\right)2\pi}{m}}\frac{{T}_{n,2}\left(t\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}}{{\mathrm{\Lambda}}_{2}^{2n-2}}\text{d}\nu \text{d}t\right), \end{align}$$where Λ22=(bt(−sinϕsinν+cosϕcosν)+δ)2${\mathrm{\Lambda}}_{2}^{2}={(\textit{bt}(-\sin\phi \sin\nu +\cos\phi \cos\nu )+\delta )}^{2}$. More generally, we compute the integral Kn${K}_{n}$ of the formB2Kn=∫01Tnt1+4b2π2c2t2×∫abbt−sinϕsinν+cosϕcosν+δ2−2ndνdt$$\begin{align} {K}_{n}& =\int_{0}^{1}{T}_{n}\left(t\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}} \nonumber\\ &\ \ \ \times \int_{a}^{b}{\left(\textit{bt}\left(-\sin\phi \sin\nu +\cos\phi \cos\nu \right)+\delta \right)}^{2-2n}\text{d}\nu \text{d}t \end{align}$$To evaluate this integral, we first combine the sinν$\sin\nu $ and cosν$\cos\nu $ terms with phasor addition,B3Kn=∫01Tnt1+4b2π2c2t2∫abδ+btcosν+θ2−2ndνdt=∫01Tnt1+4b2π2c2t2∫a+θb+θδ+btcosν2−2ndνdt$$\begin{align} {K}_{n}& =\int_{0}^{1}{T}_{n}\left(t\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\int_{a}^{b}{\left(\delta +\textit{bt}\cos\left(\nu +\theta \right)\right)}^{2-2n}\text{d}\nu \text{d}t\nonumber\\ & =\int_{0}^{1}{T}_{n}\left(t\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}\int_{a+\theta}^{b+\theta}{\left(\delta +\textit{bt}\cos\left(\nu \right)\right)}^{2-2n}\text{d}\nu \text{d}t \end{align}$$The above simplification arises asB4−sinϕsinν+cosϕcosν=−sinϕcos−π2+cosϕcos02+−sinϕsin−π2+cosϕcossin02cosν+θ$$\begin{eqnarray} &&{-}\!\sin\phi \sin\nu +\cos\phi \cos\nu \nonumber\\ &&\quad=\sqrt{ \def\eqcellsep{&}\begin{array}{l} {\left({-}\!\sin\phi \cos\left({-}\frac{\pi}{2}\right)+\cos\phi \cos 0\right)}^{2} \\ +{\left({-}\!\sin\phi \sin\left({-}\frac{\pi}{2}\right)+\cos\phi \cos\sin0\right)}^{2} \end{array} }\,\cos\left(\nu +\theta \right) \end{eqnarray}$$where θ=arctan(tan(ϕ))$\theta =\arctan(\tan(\phi ))$ when cosϕ>0$\cos\phi >0$, θ=π+arctan(tan(ϕ))$\theta =\pi +\arctan(\tan(\phi ))$ when cosϕ<0$\cos\phi <0$ and θ=sgn(sinϕ)π2$\theta =\text{sgn}(\sin\phi )\frac{\pi}{2}$ when cosϕ=0$\cos\phi =0$. Next, to evaluate Kn${K}_{n}$, we consider the integral I given byB5I=∫a+θb+θdνδ+btcosν2n−2$$\begin{align} I& =\int_{a+\theta}^{b+\theta}\frac{\text{d}\nu}{{\left(\delta +\textit{bt}\cos\nu \right)}^{2n-2}} \end{align}$$We note that since δ>bt$\delta >\textit{bt}$, the integrand do not have singularities for all values of ν. Using the standard tangent half‐angle substitution, we obtainB6I=∫z1z221+z2−1δ+bt1−z21+z22n−2dz=2∫z1z21+z22n−3δ1+z2+bt1−z22n−2dz=2∫z1z21+z22n−3δ+bt+δ−btz22n−2dz=2∫u1u21+u2n−3δ+bt+δ−btu2n−2du2u$$\begin{align} I& =\int_{{z}_{1}}^{{z}_{2}}\frac{2{\left(1+{z}^{2}\right)}^{-1}}{{\left(\delta +\textit{bt}\frac{1-{z}^{2}}{1+{z}^{2}}\right)}^{2n-2}}\text{d}z\nonumber\\ & =2\int_{{z}_{1}}^{{z}_{2}}\frac{{\left(1+{z}^{2}\right)}^{2n-3}}{{\left(\delta \left(1+{z}^{2}\right)+\textit{bt}\left(1-{z}^{2}\right)\right)}^{2n-2}}\text{d}z\nonumber\\ & =2\int_{{z}_{1}}^{{z}_{2}}\frac{{\left(1+{z}^{2}\right)}^{2n-3}}{{\left(\delta +\textit{bt}+\left(\delta -\textit{bt}\right){z}^{2}\right)}^{2n-2}}\text{d}z\nonumber\\ & =2\int_{{u}_{1}}^{{u}_{2}}\frac{{\left(1+u\right)}^{2n-3}}{{\left(\delta +\textit{bt}+\left(\delta -\textit{bt}\right)u\right)}^{2n-2}}\frac{\text{d}u}{2\sqrt{u}} \end{align}$$We now apply a linear substitution u=(u2−u1)ω+u1$u=({u}_{2}-{u}_{1})\omega +{u}_{1}$ to integrate over [0, 1], after which we binomially expand the term with power 2n−3$2n-3$ to put the integral in the form of an Appell's F1 function,B7I=u2−u1∫01u2−u1ω+u1−121+u2−u1ω+u12n−3δ+bt+δ−btu2−u1ω+u12n−2dω=u2−u1∑i=02n−32n−3!1+u12n−3−iu2−u1ii!2n−3−i!δ1+u1+bt1−u12n−2×∫01ωi1−u1−u2u1ω−121−δ−btu1−u2δ1+u1+bt1−u1ω2n−2dω=u2−u1∑i=02n−32n−3!1+u12n−3−iu2−u1ii!2n−3−i!δ1+u1+bt1−u12n−2i+2×F1i+1,12,2n−2,i+2;u1−u2u1,δ−btu1−u2δ1+u1+bt1−u1$$\begin{align} I& =\left({u}_{2}-{u}_{1}\right)\int_{0}^{1}\frac{{\left(\left({u}_{2}-{u}_{1}\right)\omega +{u}_{1}\right)}^{-\frac{1}{2}}{\left(1+\left({u}_{2}-{u}_{1}\right)\omega +{u}_{1}\right)}^{2n-3}}{{\left(\delta +\textit{bt}+\left(\delta -\textit{bt}\right)\left(\left({u}_{2}-{u}_{1}\right)\omega +{u}_{1}\right)\right)}^{2n-2}}\text{d}\omega \nonumber\\ & =\left({u}_{2}-{u}_{1}\right)\sum _{i=0}^{2n-3}\frac{\left(2n-3\right)!{\left(1+{u}_{1}\right)}^{2n-3-i}{\left({u}_{2}-{u}_{1}\right)}^{i}}{i!\left(2n-3-i\right)!{\left(\delta \left(1+{u}_{1}\right)+\textit{bt}\left(1-{u}_{1}\right)\right)}^{2n-2}} \nonumber\\ &\ \ \ \times \,\int_{0}^{1}\frac{{\omega}^{i}{\left(1-\frac{{u}_{1}-{u}_{2}}{{u}_{1}}\omega \right)}^{-\frac{1}{2}}}{{\left(1-\frac{\left(\delta -\textit{bt}\right)\left({u}_{1}-{u}_{2}\right)}{\delta \left(1+{u}_{1}\right)+\textit{bt}\left(1-{u}_{1}\right)}\omega \right)}^{2n-2}}\text{d}\omega \nonumber\\ & =\left({u}_{2}-{u}_{1}\right)\sum _{i=0}^{2n-3}\frac{\left(2n-3\right)!{\left(1+{u}_{1}\right)}^{2n-3-i}{\left({u}_{2}-{u}_{1}\right)}^{i}}{i!\left(2n-3-i\right)!{\left(\delta \left(1+{u}_{1}\right)+\textit{bt}\left(1-{u}_{1}\right)\right)}^{2n-2}\left(i+2\right)} \nonumber\\ &\ \ \ \times \,{F}_{1}\left(i+1,\frac{1}{2},2n-2,i+2;\frac{{u}_{1}-{u}_{2}}{{u}_{1}},\frac{\left(\delta -\textit{bt}\right)\left({u}_{1}-{u}_{2}\right)}{\delta \left(1+{u}_{1}\right)+\textit{bt}\left(1-{u}_{1}\right)}\right) \end{align}$$where u1=tan2(a+θ2)${u}_{1}={\tan}^{2}(\frac{a+\theta}{2})$ and u2=tan2(b+θ2)${u}_{2}={\tan}^{2}(\frac{b+\theta}{2})$.Now Kn${K}_{n}$ only has a single integral in t remained, namelyB8Kn=∫01Tnt1+4b2π2c2t2Itdt$$\begin{align} {K}_{n}& =\int_{0}^{1}{T}_{n}\left(t\right)\sqrt{1+\frac{4{b}^{2}{\pi}^{2}}{{c}^{2}}{t}^{2}}I\left(t\right)\text{d}t \end{align}$$We note that this integration is unlikely to yield analytical results so to progress further, numerical integration may be required.CAppendixMolecular Dynamics SimulationsThe large‐scale atomic/molecular massively parallel simulator (LAMMPS) software package [40] was used to perform the simulations for the above described ssDNA‐graphene interaction. Results of these simulations are reported in this paper. The system is simulated in a domain of size 100 Å ×$\ensuremath{\times{}}$ 100 Å ×$\ensuremath{\times{}}$ 100 Å. The Lennard‐Jones pair potential is used with a cut‐off distance of 14 Å. To compare simulation results with continuum models, the molecules are forced to move along prescribed trajectories as opposed to allowing the program to determine the movement from a set of initial conditions. This method of simulation is identical to a numeric solution of the discrete case of each interaction. To capture the invariance of the model to rotations about the ssDNA axis, the simulation results are the average potentials from ten random ssDNA strands simulated over ten different rotational configurations about the ssDNA axis (100 simulations in total).AcknowledgementsThe authors are grateful for the support of the Australian Research Council through funding of a Discovery Project DP170102705.Open access publishing facilitated by The University of Newcastle, as part of the Wiley ‐ The University of Newcastle agreement via the Council of Australian University Librarians.Conflict of InterestThe authors declare no conflict of interest.Data Availability StatementThe data that support the findings of this study are available from the corresponding author upon reasonable request.M. K. Antoniou, E. K. Diamanti, A. Enotiadis, A. Policicchio, K. Dimos, F. Ciuchi, E. Maccallini, D. Gournis, R. G. Agostino, Microporous Mesoporous Mater. 2014, 188, 16.P. Bénard, R. Chahine, Scr. Mater. 2007, 56, 803.M. E. Casco, M. Martínez‐Escandell, E. Gadea‐Ramos, K. Kaneko, J. Silvestre‐Albero, F. Rodríguez‐Reinoso, Chem. Mater. 2015, 27, 959.M. Feroldi, A. Neves, C. Borba, H. Alves, J. Cleaner Prod. 2018, 172, 921.P.‐S. Choi, J.‐M. Jeong, Y.‐K. Choi, M.‐S. Kim, G.‐J. Shin, S.‐J. Park, Carbon Lett. 2016, 17, 18.F. Su, C. Lu, W. Cnen, H. Bai, J. F. Hwang, Sci. Total Environ. 2009, 407, 3017.M. González Barriuso, C. Pesquera González, F. González Martínez, Á. Yedra Martínez, C. Blanco Delgado, Chem. Eng. Trans. 2019, 75, 637.J. Simon, E. Flahaut, M. Golzio, Materials 2019, 12, 624.N. Anzar, R. Hasan, M. Tyagi, N. Yadav, J. Narang, Sens. Int. 2020, 1, 100003.C. Cha, S. R. Shin, N. Annabi, M. R. Dokmeci, A. Khademhosseini, ACS Nano 2013, 7, 2891.E. Llobet, Sens. Actuators, B 2013, 179, 32.S. W. Lee, W. Lee, Y. Hong, G. Lee, D. S. Yoon, Sens. Actuators, B 2018, 255, 1788.C. Cantalini, L. Valentini, I. Armentano, J. Kenny, L. Lozzi, S. Santucci, J. Eur. Ceram. Soc. 2004, 24, 1405.A. Z. Wilczewska, K. Niemirowicz, K. H. Markiewicz, H. Car, Pharmacol. Rep. 2012, 64, 1020.M. Martincic, G. Tobias, Expert Opin. Drug Delivery 2015, 12, 563.C. Tripisciano, K. Kraemer, A. Taylor, E. Borowiak‐Palen, Chem. Phys. Lett. 2009, 478, 200.J. A. Keith, V. Vassilev‐Galindo, B. Cheng, S. Chmiela, M. Gastegger, K.‐R. Muller, A. Tkatchenko, Chem. Rev. 2021, 121, 9816.L. A. Girifalco, J. Phys. Chem. 1992, 96, 858.D. A. Barlow, F. E. Colson, Chem. Phys. Lett. 2022, 805, 139959.P. Sarapat, D. Baowan, J. M. Hill, Appl. Nanosci. 2018, 8, 537.A. Farrokhabadi, A. Koochi, M. Abadyan, Microsyst. Technol. 2014, 20, 291.E. Alvira, Chem. Phys. Lett. 2007, 439, 252.B. J. Cox, N. Thamwattana, J. M. Hill, Proc. R. Soc. London, Ser. A 2007, 463, 461.B. J. Cox, N. Thamwattana, J. M. Hill, J. Phys. A: Math. Theor. 2008, 41, 235209.D. Baowan, B. J. Cox, T. A. Hilder, J. M. Hill, N. Thamwattana, Modelling and Mechanics of Carbon‐Based Nanostructured Materials, William Andrew Publishing, Norwich, NY 2017.M. H. Alshehri, B. J. Cox, J. M. Hill, Eur. Phys. J. D 2013, 67, 226.H. Al Garalleh, N. Thamwattana, B. J. Cox, J. M. Hill, Bull. Math. Biol. 2015, 77, 184.T. Tran‐Duc, N. Thamwattana, B. J. Cox, J. M. Hill, Comput. Mater. Sci. 2010, 49, S307.O. O. Adisa, B. J. Cox, J. M. Hill, Micro Nano Lett. 2010, 5, 291.K. Stevens, N. Thamwattana, T. Tran‐Duc, Jpn. J. Appl. Phys. 2020, 128, 204301.K. Stevens, T. Tran‐Duc, N. Thamwattana, J. M. Hill, Computation 2020, 8, 107.K. Stevens, N. Thamwattana, T. Tran‐Duc, ACS Omega 2022, 7, 29773.F. Pu, J. Ren, X. Qu, Adv. Mater. 2014, 26, 5742.K. Qu, L. Wu, J. Ren, X. Qu, ACS Appl. Mater. Interfaces 2012, 4, 5001.K. Xu, X. Meshik, B. M. Nichols, E. Zakar, M. Dutta, M. A. Stroscio, Nanotechnology 2014, 25, 205501.J. C. Wang, Proc. Natl. Acad. Sci. U. S. A. 1979, 76, 200.A. K. Rappé, C. J. Casewit, K. Colwell, W. A. Goddard III, W. M. Skiff, J. Am. Chem. Soc. 1992, 114, 10024.L. A. Girifalco, M. Hodak, R. S. Lee, Phys. Rev. B 2000, 62, 13104.I. S. Gradshteyn, I. M. Ryzhik, Table of Integrals, Series, and Products, Academic Press, San Diego, CA 2014.A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, S. J. Plimpton, Comp. Phys. Comm. 2022, 271, 108171.
Journal
Advanced Theory and Simulations – Wiley
Published: May 1, 2023
Keywords: continuum modeling; DNA; graphene; Lennard‐Jones potential