Adhesion of bis-salphen-based coordination polymers to graphene: insights from free energy perturbation study

Citation: Pyrlin, S.; Lenzi, V.; Silva, A.; Ramos, M.; Marques, L. Adhesion of Bis-Salphen-Based Coordination Polymers to Graphene: Insights from Free Energy Perturbation Study. Polymers 2022,14, 4525. https:// doi.org/10.3390/polym14214525 Academic Editor: Aleksandar Y. Mehandzhiyski Received: 26 July 2022 Accepted: 18 October 2022 Published: 26 October 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). polymers Article Adhesion of Bis-Salphen-Based Coordination Polymers to Graphene: Insights from Free Energy Perturbation Study Sergey Pyrlin * , Veniero Lenzi , Alexandre Silva, Marta Ramos and Luís Marques Physics Center of Minho and Porto Universities (CF-UM-UP), University of Minho, Campus de Gualtar, 4710-057 Braga, Portugal *Correspondence: [email protected] Abstract: Manipulation of nanoscale objects using molecular self-assembly is a potent tool to achieve large scale nanopatterning with small effort. Coordination polymers of bis-salphen compounds based on zinc have demonstrated their ability to align carbon nanotubes into micro-scale networks with an unusual “rings-and-rods” pattern. This paper investigates how the compounds interact with pristine and functionalized graphene using density functional theory calculations and molecular dynamic simulations. Using the free energy perturbation method we will show how the addition of phenyl side groups to the core compound and functionalization of graphene affect the stability, mobility and conformation adopted by a dimer of bis-(Zn)salphen compound adsorbed on graphene surface and what it can reveal about the arrangement of chains of bis-(Zn)salphen polymer around carbon nanotubes during the self-assembly of microscale networks. Keywords: salphen base; metal-organic; molecular dynamics; free energy; self-assembly 1. Introduction The last decades brought significant breakthroughs in nanoscale synthesis and property engineering [ 1 – 3 ]. However, exploiting these advances often requires positioning, alignment of nanoscale objects and assembling them into higher order structures, such as lattices of nanocrystals for thermoelectrics [ 4 , 5 ], networks of nanofillers for optoelectronics and photovoltaic applications [6,7] etc. One particularly interesting way of constructing interconnected networks of nanotubes was reported by Escarcega-Bobadilla et al. [ 8 ], who observed that certain metal organic compounds can self-assemble into networks upon drop-casting from dichloromethane (DCM) upon solvent evaporation (Figure 1a). The networks were composed of microscale rings and rods interconnected in a neuron-like fashion. While the thickness of the observed rings and rods was estimated to be around 100 nm, the networks remained connected across areas of hundreds of micrometers. Even more intriguing, such self-assembling networks have demonstrated an ability to incorporate carbon nanotubes (CNTs) and align them both on surface and in the bulk of polymer, which resulted in significant boost of electrical conductivity. While both the arrangement of carbon nanotubes into circular patterns [ 9 , 10 ] and self-assembly of polymers and metal-organic compounds into circular structures [ 11 , 12 ] has been studied previously, the self-assembly of bis-(Zn)salphen compounds deserves a special attention both from fundamental point of view due to its unusual neuron-like pattern of interconnected “rings-and-rods” as well as from application standpoint as an approach to produce conductive networks from very diluted solutions of conductive fillers. The chemical structure of the compounds, for which self-assembly of networks was achieved, is shown in the Figure 1b: it consists of the two symmetrical salphen bases binding metal (Zn) cations through coordination interaction with oxygen and nitrogen atoms, connected with central biphenyl bond and functionalized with additional phenyl rings. Polymers 2022,14, 4525. https://doi.org/10.3390/polym14214525 https://www.mdpi.com/journal/polymers Polymers 2022,14, 4525 2 of 18 Polymers 2022, 14, x FOR PEER REVIEW 2 of 19 Figure 1. (a) Self-assembled molecular network of bis-(Zn)salphen compound, obtained by dropcasting from DCM (preprinted with a permission from Ref. [8]); (b) Chemical structure of bis- (Zn)salphen compound, functional phenyl rings are highlighted with yellow; (c), 3D rendering of a coordination bonded dimer of the compound and a fragment of a polymer chain of such dimers Ref [13]. The chemical structure of the compounds, for which self-assembly of networks was achieved, is shown in the Figure 1b: it consists of the two symmetrical salphen bases binding metal (Zn) cations through coordination interaction with oxygen and nitrogen atoms, connected with central biphenyl bond and functionalized with additional phenyl rings. A subsequent study [13] has explained this phenomenon by the formation of a coordination polymer, the periodic units of which are dimers of the bis-(Zn)salphen compounds (Figure 1c). Using atomistic simulations it was shown that such chains can form semiflexible fibers via π-π interactions due to phenyl side groups. It is expected that upon solvent evaporation such fibers aggregate and form networks of interconnected loops and toroidal globules in accordance with experimental observations. However, the interaction of bis-(Zn)salphen compounds and carbon nanotubes, more specifically their incorporation in the core of the self-assembled networks, remains unexplored. This aspect is especially crucial for applications as electrical conductivity of CNT-based films and composites depends drastically on the contact between nanotubes [14,15]. Furthermore, molecular rearrangement in the nanoscale gaps between carbon nanotubes can alter mechanical and transport properties of surrounding matrix material[16]. In this paper we will provide an insight to this phenomenon from atomistic calculations. The free energy is a key quantity to understand the physical processes in a solvent environment, such as conformational changes, protein-ligand binding and moleculesubstrate adsorption[17]. A numerical estimate of the free energy difference (Δ𝐺→) between two states i and f, characterized by potential energies 𝑈 and 𝑈, close enough in the configuration space, is given by the Zwanzig formula[18], which is at the basis of many free energy estimation methods: Δ𝐺→ =−𝑘𝑇 𝑙𝑛〈𝑒𝑥𝑝󰇧−𝑈−𝑈 𝑘𝑇󰇨〉, (1) where T is the ambient temperature, 𝑘 is the Boltzmann constant and <> denotes an ensemble average by the microstates of the starting state i. The main factor hindering the application of Equation (1) is high computational cost required for its calculation, as a good sampling over the configuration space is required. As the Zwanzig equation is valid only for states that are close enough in the configuration space, a typical transition needs to be split into many intermediate steps for which Equation (1) holds, and in doing so the free energy along the resulting path can be calculated, which is referred to as the potential of mean force (PMF). For instance, Umbrella Sampling method [19,20] consists of simulations of the system dynamics, constrained by an external elastic force to one of a set of overlapping windows, subdividing the path from initial to final state. Then a weighted Figure 1. ( a ) Self-assembled molecular network of bis-(Zn)salphen compound, obtained by drop-casting from DCM (preprinted with a permission from Ref. [ 8 ]); ( b ) Chemical structure of bis-(Zn)salphen compound, functional phenyl rings are highlighted with yellow; ( c ) 3D rendering of a coordination bonded dimer of the compound and a fragment of a polymer chain of such dimers Ref [13]. A subsequent study [ 13 ] has explained this phenomenon by the formation of a coordination polymer, the periodic units of which are dimers of the bis-(Zn)salphen compounds (Figure 1c). Using atomistic simulations it was shown that such chains can form semiflexible fibers via π - π interactions due to phenyl side groups. It is expected that upon solvent evaporation such fibers aggregate and form networks of interconnected loops and toroidal globules in accordance with experimental observations. However, the interaction of bis-(Zn)salphen compounds and carbon nanotubes, more specifically their incorporation in the core of the self-assembled networks, remains unexplored. This aspect is especially crucial for applications as electrical conductivity of CNTbased films and composites depends drastically on the contact between nanotubes [ 14 , 15 ]. Furthermore, molecular rearrangement in the nanoscale gaps between carbon nanotubes can alter mechanical and transport properties of surrounding matrix material [ 16 ]. In this paper we will provide an insight to this phenomenon from atomistic calculations. The free energy is a key quantity to understand the physical processes in a solvent environment, such as conformational changes, protein-ligand binding and molecule-substrate adsorption [ 17 ]. A numerical estimate of the free energy difference ( ∆Gi→f ) between two states iand f, characterized by potential energies Ui and Uf , close enough in the configuration space, is given by the Zwanzig formula [ 18 ], which is at the basis of many free energy estimation methods: ∆Gi→f=−kBT lnhexp −(Ui−Uf) kBT!i, (1) where Tis the ambient temperature, kB is the Boltzmann constant and <> denotes an ensemble average by the microstates of the starting state i. The main factor hindering the application of Equation (1) is high computational cost required for its calculation, as a good sampling over the configuration space is required. As the Zwanzig equation is valid only for states that are close enough in the configuration space, a typical transition needs to be split into many intermediate steps for which Equation (1) holds, and in doing so the free energy along the resulting path can be calculated, which is referred to as the potential of mean force (PMF). For instance, Umbrella Sampling method [ 19 , 20 ] consists of simulations of the system dynamics, constrained by an external elastic force to one of a set of overlapping windows, subdividing the path from initial to final state. Then a weighted histogram analysis method (WHAM) [ 21 ] is used to reconstruct PMF, taking into account the restraining potentials applied within each window. Another way to estimate ∆Gi→f is to estimate the work needed to force the system from one state to another ( Wi→f ) in an non-equilibrium process, which is related to a free energy change via Jarzynski’s equality ( exp−∆Gi→f/kBT = hexp−Wi→f/kBTi ) [ 22 ]. Then one could perform a simulation including a forcing term to drive the dynamics and track the forces and displacements Polymers 2022,14, 4525 3 of 18 to estimate the work. This approach is called steered molecular dynamics (SMD) [ 23 , 24 ]. On the other hand, the adaptive biasing force method [ 25 ] implies altering the free energy landscape of a process in such a way to nullify any difference in free energy along the path. This way, the PMF is reconstructed from the amount of biasing potential used. Alternatively, the free energy perturbation method (FEP) [ 26 – 28 ] suggests to extend the systems potential energy into a function of a coupling parameter ( λ ) in such a way that it gradually changes between the initial and final states: U(λ) = λUf+(1−λ)Ui . In this case the transition consists of simulations of system dynamics under potential corresponding to a set of intermediate values of λ , rather than forcefully drugging the molecules along a predefined path. The transition free energy can then be evaluated using the Bennett acceptance ratio method [ 29 , 30 ]. FEP was shown to give very precise results with the relative errors of the order of 1–2 kcal/mol [26]. An intriguing way to use FEP to reduce the computational cost of free energy calculations bases on the fact that the Gibbs free energy differences depend only on the initial and final states and not on the path connecting them. Such approach, known as alchemical transformations [ 31 ] substitutes estimating ∆Gi→f of a transition that can be harder to trace directly by introducing an intermediate state, related to the initial and final states by transitions, which may be unphysical, but computationally more convenient. Thus, solvation of a new compound can be studied by placing both solute and solvent molecules into the same simulation box and vary the strength of the solute-solvent interaction potential by using a coupling parameter λ as a scaling factor. Scaling interaction up ( λ=0→1 ) will correspond to immersion of a solute into solvent (solvation). In certain cases, however, it is more convenient to model the opposite process—decoupling (or extraction) of solute from solvent ( λ=1→0 ). Thus, free energy of a reaction in solvent can be estimated using the solvation or decoupling free energies of the reactants and products and free energy of the same reaction in vacuum. This approach has been successfully used in computational drug discovery simulations to determine drug-ligand binding free energies [ 32 , 33 ], as well as to investigate the exfoliation properties of 2D materials in various solvents [34]. In this study we combine density functional theory and free energy perturbation calculations to compare possible mechanisms of adhesion of Zn-salphen compounds to carbon nanotube surface. Given a significant difference in scale between CNT diameter (tens of nm) and size of a bis-(Zn)salphen compound (>2 nm), we replace nanotube surface with a periodic sheet of plane graphene. Focusing on the atomistic details of coordination polymer adhesion to CNT we model the interacting tip of a polymer chain with its single building block—a dimer of bis-(Zn)salphen compounds. In the following sections we will describe the calculation setup used in our simulations and the effect of compound and graphene functionalization on the strength of adhesion and geometry of the adsorbed compounds in dichloromethane solvent. 2. Materials and Methods 2.1. Ab Initio Energy of Adhesion To estimate the energy of adhesion of the (Zn)salphen-based compounds to graphene surface we have performed ab initio calculations for the systems consisting of single salphen base compound of Zn and a periodic sheet of pristine or functionalized graphene. Since the chemical structure of a salphen ligand is rich in aromatic rings, π - π interaction is a primary factor of adhesion. However, compound’s cation center surrounded by negative oxygen and nitrogen atoms would be attracted to charged impurities on graphene surface. For this reason, besides pristine graphene we include sheets modified with: boron and nitrogen substitutions, chemisorbed oxygen, hydroxyl and carboxyl groups, which have been reported to improve graphene adhesive properties [35–37]. System under study (Figure 2) was enclosed in (approximately) 2.48 × 2.58 × 1.29 nm periodic box, which is enough to accommodate a supercell of graphene with 240 atoms (6 by 10 rectangular unit cells of the hexagonal graphene structure) and ensure lack of interaction of the compound with its periodic images. Polymers 2022,14, 4525 4 of 18 Polymers 2022, 14, x FOR PEER REVIEW 4 of 19 negative oxygen and nitrogen atoms would be attracted to charged impurities on graphene surface. For this reason, besides pristine graphene we include sheets modified with: boron and nitrogen substitutions, chemisorbed oxygen, hydroxyl and carboxyl groups, which have been reported to improve graphene adhesive properties [35–37]. System under study (Figure 2) was enclosed in (approximately) 2.48 × 2.58 × 1.29 nm periodic box, which is enough to accommodate a supercell of graphene with 240 atoms (6 by 10 rectangular unit cells of the hexagonal graphene structure) and ensure lack of interaction of the compound with its periodic images. Figure 2. Supercell of graphene (grey tubes) with adsorbed (Zn)salphen compound (ball&stick, color by element) used in DFT simulations: (a) side view, (b) top view. Periodic image of the graphene beyond the simulation region (blue box) lattice is shown with light grey. Graphene atoms shown with ball and stick representation enclosed in a grey rectangle indicate a rectangle unit cell of a graphene lattice. The calculations were performed with the SIESTA DFT package (v. 4.1)[38]. Exchange and correlation functionals were described by generalized gradient approximation in Perdew-Becke-Ernzerhof approximation[39]. Valence electrons orbitals were simulated by a basis set of triple-zeta triply polarized numerical atomic orbitals with energy shift 10 meV [40]. Core electrons were replaced by Troullier-Martins pseudopotentials[41]. Cut-off of 1000 Ry was used to define integration mesh. Van der Waals interaction was included using Grimme potentials [42]. [42] Counterpoise correction to exclude basis superposition error [43]. Energy minimization was performed with respect to atomic positions and lengths of periodic box sides. Geometry of standalone graphene and salphen-Zn compound was optimized until total energy was converged to 0.001 eV, maximum atomic force did not exceed 0.01 eV/Å and maximum atomic displacement was 0.01 Bohr. Correspondingly, the weakly interacting graphene–salphen system was converged to 0.01 eV, 0.1 eV/Å and 0.05 Bohr. The set of calculations included geometry optimizations to obtain relaxed structures followed by single point energy calculations to get ground state energies of the standalone compound (𝐸 ) and graphene sheet (𝐸 ), as well as of the complex of the compound adsorbed on the graphene surface (𝐸 ). Furthermore, additional single point energy calculations were performed for each of compound and graphene in the geometry of the complex to estimate the basis superposition error for each counterpart (𝐸 , ). Then the energy of adhesion for graphene-compound system can be calculated according to: 𝐸 =𝐸 −𝐸 −𝐸 −󰇛𝐸 +𝐸 󰇜. (2) 2.2. Molecular Dynamic Simulations To estimate the free energy of adhesion of bis-(Zn)salphen dimers to graphene surface we have performed all-atom free energy perturbation simulations with LAMMPS package [44]. The solvent (Dichloromethane, DCM) molecules were treated explicitly. Non-bonded (van der Waals and electrostatic) interactions were modelled using 6–12 Figure 2. Supercell of graphene (grey tubes) with adsorbed (Zn)salphen compound (ball&stick, color by element) used in DFT simulations: ( a ) side view, ( b ) top view. Periodic image of the graphene beyond the simulation region (blue box) lattice is shown with light grey. Graphene atoms shown with ball and stick representation enclosed in a grey rectangle indicate a rectangle unit cell of a graphene lattice. The calculations were performed with the SIESTA DFT package (v. 4.1) [ 38 ]. Exchange and correlation functionals were described by generalized gradient approximation in Perdew-Becke-Ernzerhof approximation [ 39 ]. Valence electrons orbitals were simulated by a basis set of triple-zeta triply polarized numerical atomic orbitals with energy shift 10 meV [ 40 ]. Core electrons were replaced by Troullier-Martins pseudopotentials [ 41 ]. Cutoff of 1000 Ry was used to define integration mesh. Van der Waals interaction was included using Grimme potentials [ 42 ]. Counterpoise correction to exclude basis superposition error [ 43 ]. Energy minimization was performed with respect to atomic positions and lengths of periodic box sides. Geometry of standalone graphene and salphen-Zn compound was optimized until total energy was converged to 0.001 eV, maximum atomic force did not exceed 0.01 eV/ ˚ A and maximum atomic displacement was 0.01 Bohr. Correspondingly, the weakly interacting graphene–salphen system was converged to 0.01 eV, 0.1 eV/ ˚ A and 0.05 Bohr. The set of calculations included geometry optimizations to obtain relaxed structures followed by single point energy calculations to get ground state energies of the standalone compound ( EC g ) and graphene sheet ( EG g ), as well as of the complex of the compound adsorbed on the graphene surface ( EC+G g ). Furthermore, additional single point energy calculations were performed for each of compound and graphene in the geometry of the complex to estimate the basis superposition error for each counterpart ( EC,G BSSE ). Then the energy of adhesion for graphene-compound system can be calculated according to: Eadh =EC+G g−EC g−EG g−EC BSSE +EG BSSE. (2) 2.2. Molecular Dynamic Simulations To estimate the free energy of adhesion of bis-(Zn)salphen dimers to graphene surface we have performed all-atom free energy perturbation simulations with LAMMPS package [ 44 ]. The solvent (Dichloromethane, DCM) molecules were treated explicitly. Non-bonded (van der Waals and electrostatic) interactions were modelled using 6–12 Lennard-Jones and Coulomb potentials. A cut-off of 10 Å was used for short range part of both potentials. The long range interaction was treated using particle-particle particle-mesh solver [ 45 ]. Electrostatic interactions were modelled using fixed atomic charges, fitted by RESP method [ 46 ]. Electrostatic potential grid was obtained from ab initio calculations using Gaussian package [ 47 ]. The interactions between Zn 2+ ion and organic ligand were treated within a free cation model [ 48 ] in order to allow the compound to adopt more relaxed conformation while interacting with hydroxyl and carboxyl groups of functionalized graphene [ 46 ]. During FEP simulation “soft” version of van der Waals potential was used Polymers 2022,14, 4525 5 of 18 to gradually scale down the strength of interactions without creating a large deviation from equilibrium [49]. Both covalent and van der Waals interactions were described by OPLS-AA force field [ 50 , 51 ], with the exception for the atoms of the imine group (-C=N-), explicit parametrization of which was introduced in OPLS_2005 force field [ 52 ] and is essential for accurate description of salphen compounds. While many improvements were introduced over the years with more recent parametrizations [ 53 – 56 ], such as more precise torsions for DNA backbone chain and better description of nucleic acid interaction, the OPLS-AA force field is still widely used for molecular simulations in liquid state, such as including recent studies of SARS-CoV-2 inhibitors [ 57 ], chromophores [ 58 ] and electrolytes for novel batteries [ 59 ]. It is also implemented in many open source toolkits for building input files for molecular dynamics simulations, such as Moltemplate [ 60 ], used in this work. Although virtual charge sites, introduced in latest parametrizations [ 55 , 56 ] could further improve the description of hydrogen bonding to functional groups—this will be investigated in future work. The system under study consisted of a dimer of bis-(Zn)salphen compound, put at random orientation on top of 4 layered graphene sheet of 12 by 20 rectangular unit cells of the hexagonal graphene structure (approx. 5 × 5 nm), periodic in 2 in-plane directions (x&y). Due to the symmetry of the bis-salphen dimer the number of its orientations relative to graphene surface which needs to be studied can be reduced. To obtain initial positions the dimer was placed in one of the orientations shown in the Figure 3and rotated randomly around the axis perpendicular to the graphene surface (z). The dimer was placed in each orientation and relaxed 5 times, resulting in 20 independent simulations. Polymers 2022, 14, x FOR PEER REVIEW 5 of 19 Lennard-Jones and Coulomb potentials. A cut-off of 10 Å was used for short range part of both potentials. The long range interaction was treated using particle-particle particlemesh solver [45]. Electrostatic interactions were modelled using fixed atomic charges, fitted by RESP method [46]. Electrostatic potential grid was obtained from ab initio calculations using Gaussian package [47]. The interactions between Zn 2+ ion and organic ligand were treated within a free cation model [48] in order to allow the compound to adopt more relaxed conformation while interacting with hydroxyl and carboxyl groups of functionalized graphene [46]. During FEP simulation “soft” version of van der Waals potential was used to gradually scale down the strength of interactions without creating a large deviation from equilibrium [49]. Both covalent and van der Waals interactions were described by OPLS-AA force field [50,51], with the exception for the atoms of the imine group (-C=N-), explicit parametrization of which was introduced in OPLS_2005 force field [52] and is essential for accurate description of salphen compounds. While many improvements were introduced over the years with more recent parametrizations [53–56], such as more precise torsions for DNA backbone chain and better description of nucleic acid interaction, the OPLS-AA force field is still widely used for molecular simulations in liquid state, such as including recent studies of SARS-CoV-2 inhibitors [57], chromophores [58] and electrolytes for novel batteries [59]. It is also implemented in many open source toolkits for building input files for molecular dynamics simulations, such as Moltemplate [60], used in this work. Although virtual charge sites, introduced in latest parametrizations [55,56] could further improve the description of hydrogen bonding to functional groups— this will be investigated in future work. The system under study consisted of a dimer of bis-(Zn)salphen compound, put at random orientation on top of 4 layered graphene sheet of 12 by 20 rectangular unit cells of the hexagonal graphene structure (approx. 5 × 5 nm), periodic in 2 in-plane directions (x & y). Due to the symmetry of the bis-salphen dimer the number of its orientations relative to graphene surface which needs to be studied can be reduced. To obtain initial positions the dimer was placed in one of the orientations shown in the Figure 3 and rotated randomly around the axis perpendicular to the graphene surface (z). The dimer was placed in each orientation and relaxed 5 times, resulting in 20 independent simulations. Figure 3. Placement of the bis-(Zn)salphen dimer with functional phenyl rings (highlighted with yellow) atop of graphene surface: (a) a compound dimer with axes indicating orientation, (b–e) variants of placement (after relaxation), axes indicate rotation of dimer. Only topmost graphene layer is shown for clarity, also solvent molecules and hydrogen atoms of compound are omitted. For each independent run the movement of atoms of the dimer were first simulated in vacuum at 300K then laid down onto graphene surface and relaxed. After which the rest of the space on top of graphene was filled with DCM molecules up to 5 nm in height. The initial positions of solvent molecules were taken from an independent simulation of dichloromethane squeezed between graphene in the same simulation box to ensure close to equilibrium initial conditions, the solvent molecules intersecting with the compound Figure 3. Placement of the bis-(Zn)salphen dimer with functional phenyl rings (highlighted with yellow) atop of graphene surface: ( a ) a compound dimer with axes indicating orientation, ( b – e ) variants of placement (after relaxation), axes indicate rotation of dimer. Only topmost graphene layer is shown for clarity, also solvent molecules and hydrogen atoms of compound are omitted. For each independent run the movement of atoms of the dimer were first simulated in vacuum at 300K then laid down onto graphene surface and relaxed. After which the rest of the space on top of graphene was filled with DCM molecules up to 5 nm in height. The initial positions of solvent molecules were taken from an independent simulation of dichloromethane squeezed between graphene in the same simulation box to ensure close to equilibrium initial conditions, the solvent molecules intersecting with the compound molecules were excluded. Temperature and pressure of solvent were controlled by NoseHoover thermostat and barostat [ 61 ], with the atmospheric pressure applied in the z direction, perpendicular to the graphene surface, while the simulation box bounds the two in-plane directions were kept fixed. The motion of the atoms of dimers and graphene were controlled separately by Langevin thermostats [62]. Furthermore, during the first 500 fs of solvent equilibration elastic restraints with gradually decreasing spring constants were applied to the heavy atoms of the compound to avoid the initial relaxed state being compromised by excessive stress from unequilibrated solvent. After initial stress relaxation and equilibration during 250 ps with timestep of Polymers 2022,14, 4525 6 of 18 0.5 fs, the movement of hydrogen atoms was simulated using SHAKE algorithm [ 62 ] and the timestep increased to 2 fs for a 2 ns equilibration run. Besides that, to get the reference solvated state the same bis-salphen dimer was placed in a 5 × 5 × 5 nm solvent box without graphene for another 5 independent simulations. Periodic boundary conditions were applied to all surfaces of simulation box. After initial 2 ns equilibration a FEP simulation was conducted for each system. The simulation consisted of the two stages. On the first the electrostatic interaction between the compound and surrounding environment was gradually scaled down during approximately 2.0 ns of simulation time (20 intermediate λ -values each simulated during 0.1 ns). On the second stage van der Waals interactions between the atoms of compound and environment were scaled down. To make the transition smoother on the second stage the compound was decoupled atom by atom starting from exterior hydrogen atoms and finishing with the core Zn, O, N atoms. For each atom 10 intermediate λ -values were used, each simulated during 0.025 ns. Total simulation time was 40–50 ns depending on the variation of the bis-salphen compound (with or without phenyl groups). Every 10 simulation steps (20 fs) backward and forward difference of interaction potential was recorded and subsequently processed using the Bennett acceptance ratio method to estimate accumulated change of Gibbs free energy ∆G . Accumulated change of the Gibbs free energy was estimated from the recorded backward and forwards differences of interaction potential using ParseFEP plugin of VMD [27,63]. The thermodynamic cycle corresponding to this study is shown in the Figure 4. The red arrow indicates the transition of interest between solvated (c) and adsorbed compound (a) as the free energy change between these states contains the free energy of adhesion ( ∆Gadh ). While it is possible to estimate it directly by steered molecular dynamics, due to large capacity of the compound dimer to deform and change conformation a large number of slow simulations would be required to accurately sample the PMF change during slow approximation or detachment of the compound molecules. Instead, we estimate the ∆Gadh indirectly by: ∆Gadh +∆Gsolv G=−∆Gdec S&G+∆Gdec G+∆Gsolv S, (3) where ∆Gsolv G —free energy of solvation of graphene sheets, ∆Gdec S&G —free energy change during decoupling of the adsorbed compound dimer from surrounding solvent and graphene molecules, ∆Gdec G —free energy change during decoupling of the graphene from solvent and ∆Gsolv S —free energy of solvation of the compound. Since within FEP formalism, decoupling of the compound from solvent is the opposite of solvation (immersion into solvent), the solvation free energy is the negative of the decoupling free energy: ∆Gsolv =−∆Gdec . Using this substitution in Equation (3) the desired free energy of adhesion to graphene in solvent can be estimated as follows: ∆Gadh =∆Gdec S−∆Gdec S&G. (4) Using the decoupling transition to estimate the free energy of adhesion in this case is technically simpler than solvation (immersion) as in the latter process it is not guaranteed that during finite simulation time the compound will adopt the desired conformation on top of the graphene surface. On the contrary, in the case of decoupling process both starting and final states are well defined and if the set of initial positions of the adsorbed compound molecules contained the most favorable position, the largest estimated ∆Gdec across the simulations will produce a good estimate of −∆Gadh. Polymers 2022,14, 4525 7 of 18 Polymers 2022, 14, x FOR PEER REVIEW 7 of 19 Figure 4. Thermodynamic cycle used to compare the adhesion free energy of the compound to graphene in solvent: (a) a bis-(Zn)salphen dimer interacting with graphene surface, (b) graphene with the dimer decoupled (shown with transparency), (c) dimer solvated in DCM, (d) dimer decoupled from solvent. The red arrow (c-a) indicates the transformation, the free energy of which contains the target contribution, black arrows (a-b and d-c) indicate transformations which are simulated with FEP method and the gray arrow (b-d)—transformation, the free energy of which is cancelled out. Using the decoupling transition to estimate the free energy of adhesion in this case is technically simpler than solvation (immersion) as in the latter process it is not guaranteed that during finite simulation time the compound will adopt the desired conformation on top of the graphene surface. On the contrary, in the case of decoupling process both starting and final states are well defined and if the set of initial positions of the adsorbed compound molecules contained the most favorable position, the largest estimated ∆𝐺 across the simulations will produce a good estimate of −∆𝐺. A slightly different approach was used in case of the functionalized graphene (Figure 5). Rather than simulating again the lengthy decoupling of the compound dimer from graphene and solvent, we have estimated the free energy change during removing the functional groups (-OH and -COOH) from graphene and adsorbed compound (∆𝐺∗ ). In this case both Coulomb and van der Waals interactions were scaled (up or down) simultaneously for all group atoms through 20 intermediate λ -values, simulated during 0.1 ns each. Figure 5. Thermodynamic cycle used to compare the adhesion free energy of the compound to functionalized graphene in solvent: (a) a bis-(Zn)salphen dimer interacting with functionalized (- Figure 4. Thermodynamic cycle used to compare the adhesion free energy of the compound to graphene in solvent: (a) a bis-(Zn)salphen dimer interacting with graphene surface, (b) graphene with the dimer decoupled (shown with transparency), (c) dimer solvated in DCM, (d) dimer decoupled from solvent. The red arrow (c-a) indicates the transformation, the free energy of which contains the target contribution, black arrows (a-b and d-c) indicate transformations which are simulated with FEP method and the gray arrow (b-d)—transformation, the free energy of which is cancelled out. A slightly different approach was used in case of the functionalized graphene (Figure 5). Rather than simulating again the lengthy decoupling of the compound dimer from graphene and solvent, we have estimated the free energy change during removing the functional groups (-OH and -COOH) from graphene and adsorbed compound ( ∆Gdec F∗ ). In this case both Coulomb and van der Waals interactions were scaled (up or down) simultaneously for all group atoms through 20 intermediate λ -values, simulated during 0.1 ns each. Polymers 2022, 14, x FOR PEER REVIEW 7 of 19 Figure 4. Thermodynamic cycle used to compare the adhesion free energy of the compound to graphene in solvent: (a) a bis-(Zn)salphen dimer interacting with graphene surface, (b) graphene with the dimer decoupled (shown with transparency), (c) dimer solvated in DCM, (d) dimer decoupled from solvent. The red arrow (c-a) indicates the transformation, the free energy of which contains the target contribution, black arrows (a-b and d-c) indicate transformations which are simulated with FEP method and the gray arrow (b-d)—transformation, the free energy of which is cancelled out. Using the decoupling transition to estimate the free energy of adhesion in this case is technically simpler than solvation (immersion) as in the latter process it is not guaranteed that during finite simulation time the compound will adopt the desired conformation on top of the graphene surface. On the contrary, in the case of decoupling process both starting and final states are well defined and if the set of initial positions of the adsorbed compound molecules contained the most favorable position, the largest estimated ∆𝐺 across the simulations will produce a good estimate of −∆𝐺. A slightly different approach was used in case of the functionalized graphene (Figure 5). Rather than simulating again the lengthy decoupling of the compound dimer from graphene and solvent, we have estimated the free energy change during removing the functional groups (-OH and -COOH) from graphene and adsorbed compound (∆𝐺∗ ). In this case both Coulomb and van der Waals interactions were scaled (up or down) simultaneously for all group atoms through 20 intermediate λ -values, simulated during 0.1 ns each. Figure 5. Thermodynamic cycle used to compare the adhesion free energy of the compound to functionalized graphene in solvent: (a) a bis-(Zn)salphen dimer interacting with functionalized (- Figure 5. Thermodynamic cycle used to compare the adhesion free energy of the compound to functionalized graphene in solvent: (a) a bis-(Zn)salphen dimer interacting with functionalized (-COOH) graphene, (b) dimer and graphene with the functional group decoupled (shown with transparency), (c) functionalized graphene surface in DCM solvent, (d) graphene surface with the functional group decoupled from surrounding. The red arrow (c-a) indicates the transformation, the free energy of which is the target, black arrows (a-b and d-c) indicate transformations which are simulated with FEP method and the gray arrow (b-d)—transformation, the free energy of which is cancelled out. Polymers 2022,14, 4525 8 of 18 In this case, approximating the free energy of immersion of the compound dimer on top of functionalized graphene by ∆Gsolv S&G+∆∆Gadh , where ∆∆Gadh is the correction to adhesive energy due to interaction of compound with the functional group, the transition free energies are linked with the following relation: ∆Gsolv S&G+∆∆Gadh =−∆Gdec F∗+∆Gdec S&G+∆Gsolv F. (5) Using the same equality between solvation and decoupling free energies, it is straightforward from (5) that the adhesive energy correction can be estimated as: ∆∆Gadh =−∆Gdec F∗+∆Gsolv F. (6) where ∆Gsolv F —is the free energy change due to immersion of the functional group on top of graphene in solvent. Since the position of the functional group is well defined by a chemical bond to graphene, here we calculate ∆Gsolv F directly, rather than through decoupling process as we did for compound solvation. Furthermore, immersion (in case of ∆Gsolv F ) and decoupling (in case of ∆Gdec F∗ ) of both electrostatic and van der Waals interactions are performed simultaneously during 2 ns simulations. A set of 10 immersion and 20 decoupling simulations were performed for each of the selected functional groups. The initial guesses for the positions of the bis-salphen dimer, interacting with the functional groups were obtained by placing of the separately optimized geometries of the dimers so that the coordinates of the Zn, O and N atoms of one of its non-bonded salphen groups closely approximate the position of the corresponding atoms of the single Zn-salphen base interacting with -OH or -COOH functional group of graphene patch from DFT study. To obtain the relaxed initial positions vacuum simulations were performed in which the hydrogen and oxygen atoms of the functional groups were initially attached to the oxygen and zinc atoms of the compound with gradually relaxing spring force. 3. Results 3.1. Binding of a Single Salphen Complex to Functionalized Graphene To compare the influence of functionalization on the adhesion of the (Zn)salphenbased compounds to graphene surface we have performed ab initio calculations for a single compound adsorbed on the surface of pristine and functionalized graphene as shown in the Figure 2. These calculations served to analyze which type of interactions is mostly responsible for the adhesion and select the systems of interest for molecular dynamics study. The optimized geometries of the compound-graphene system are shown in the Figure 6. The estimated energies of adhesion are summarized in the Table 1. On the surface of a pristine graphene the adsorbed (Zn)salphen compound assumes completely flat shape lying parallel to the graphene surface at the distance of 3.04 Å. This together with the high negative value of the energy of adhesion of a (Zn)salphen compound to pristine graphene, indicates that the π - π interactions between graphene surface and aromatic rings of the salphen base play the major part in attraction of the compound. Doping graphene with B or N impurities provides a moderate improvement of adhesion due to electrostatic attraction between the impurity and the Zn or O atoms of salphen compound without affecting the π - π interactions: the flat shape of the compound and separating distance of 3.04 Å is maintained. At the same time, the presence of a carbonyl oxygen, i.e., connected to the two graphene atoms, has an opposite effect. Despite the ability of Zn to form additional coordination bonds, connecting to out-of-plane O 2− center of the oxidized graphene requires the salphen compound to leave its favorable flat conformation and overcome the potential barrier created by N and O atoms, surrounding Zn in salphen complex. Ass the result, the average distance between salphen aromatic rings and graphene surface increases to 3.23 Å and adhesion energy decreases by 3.3 kcal/mol. Polymers 2022,14, 4525 9 of 18 Polymers 2022, 14, x FOR PEER REVIEW 9 of 19 Figure 6. Found relaxed geometries for a single base Zn-salphen compound adsorbed to a pristine graphene (a), as well as boron (b) and nitrogen (c) substituted, oxidized (d) and functionalized with hydroxyl (e) and carboxyl (f) groups. Annotations in (e) and (f) show hydrogen bond lengths. Table 1. Ab initio energy of adhesion (kcal/mol). Functionalization Type E adh E adh -E adh (Pristine) Pristine −38.8 - B-doped −40.7 −1.9 N-doped −41.1 −2.2 =O −35.5 +3.3 -OH −40.5 −1.6 -COOH −44.6 −5.6 On the surface of a pristine graphene the adsorbed (Zn)salphen compound assumes completely flat shape lying parallel to the graphene surface at the distance of 3.04 Å. This together with the high negative value of the energy of adhesion of a (Zn)salphen compound to pristine graphene, indicates that the π-π interactions between graphene surface and aromatic rings of the salphen base play the major part in attraction of the compound. Doping graphene with B or N impurities provides a moderate improvement of adhesion due to electrostatic attraction between the impurity and the Zn or O atoms of salphen compound without affecting the π-π interactions: the flat shape of the compound and separating distance of 3.04 Å is maintained. At the same time, the presence of a carbonyl oxygen, i.e., connected to the two graphene atoms, has an opposite effect. Despite the ability of Zn to form additional coordination bonds, connecting to out-of-plane O 2center of the oxidized graphene requires the salphen compound to leave its favorable flat conformation and overcome the potential barrier created by N and O atoms, surrounding Zn in salphen complex. Ass the result, the average distance between salphen aromatic rings and graphene surface increases to 3.23 Å and adhesion energy decreases by 3.3 kcal/mol. On the contrary, functionalizing graphene with hydroxyl (-OH) and carboxyl (- COOH) groups makes the compound to bind better to the graphene surface, as these groups form hydrogen bonds with the salphen compound. Hydrogen bond lengths are estimated as 1.75 Å in case of hydroxyl group and 1.50 Å in case of carboxyl group. However, while both groups provide a single hydrogen bond, the benefit is significantly higher in case of carboxyl group. The reason for this is that hydroxyl group is situated Figure 6. Found relaxed geometries for a single base Zn-salphen compound adsorbed to a pristine graphene ( a ), as well as boron ( b ) and nitrogen ( c ) substituted, oxidized ( d ) and functionalized with hydroxyl (e) and carboxyl (f) groups. Annotations in (e) and (f) show hydrogen bond lengths. Table 1. Ab initio energy of adhesion (kcal/mol). Functionalization Type Eadh Eadh-Eadh (Pristine) Pristine −38.8 - B-doped −40.7 −1.9 N-doped −41.1 −2.2 =O −35.5 +3.3 -OH −40.5 −1.6 -COOH −44.6 −5.6 On the contrary, functionalizing graphene with hydroxyl (-OH) and carboxyl (-COOH) groups makes the compound to bind better to the graphene surface, as these groups form hydrogen bonds with the salphen compound. Hydrogen bond lengths are estimated as 1.75 Å in case of hydroxyl group and 1.50 Å in case of carboxyl group. However, while both groups provide a single hydrogen bond, the benefit is significantly higher in case of carboxyl group. The reason for this is that hydroxyl group is situated closer to the graphene surface and attraction of one of salphen’s oxygens to hydroxyl hydrogen simultaneously with repulsion of the other salphen’s oxygen from the oxygen of the functional group creates a distortion of the plane shape of salphen compound, counteracting with the π - π interactions. On the contrary, the higher position of the carboxyl group and flexibility of rotation of its bonds allows the hydrogen bond to be formed without distorting the π - π interactions: salphen compound maintains its flat shape and distance to graphene of 3.06 Å. Although the observed improvement of binding energy is comparatively small in comparison with the impact of aromatic interactions in vacuum, in a good solvent, such as dichloromethane (DCM), the molecules of the compound will experience a competing attraction by the solvent molecules. In such case these additional interactions with functionalized graphene may give a significant boost in favor of adsorption. 3.2. Free Energy of Adhesion of a Bis-Saplhen Dimer A better criterion for the tendency of a molecule to be adsorbed on a surface is the Gibbs free energy of adsorption ( ∆Gadh ). Besides potential energy of attraction it includes temperature and entropy contributions, which are crucial to describe the behavior of polymers in solution. Here using the FEP method we compare the effect of adding functional groups to the chemical structure of bis-salphen compounds and graphene Polymers 2022,14, 4525 16 of 18 3. Xu, Z.; Buehler, M.J. Nanoengineering Heat Transfer Performance at Carbon Nanotube Interfaces. ACS Nano 2009 ,3, 2767–2775. [CrossRef] [PubMed] 4. Musland, L.; Flage-Larsen, E. Thermoelectric Effect in Superlattices; Applicability of Coherent and Incoherent Transport Models. Comput. Mater. Sci. 2018,153, 88–96. [CrossRef] 5. Coropceanu, I.; Boles, M.A.; Talapin, D.V. Systematic Mapping of Binary Nanocrystal Superlattices: The Role of Topology in Phase Selection. J. Am. Chem. Soc. 2019; Just Accepted Manuscript. [CrossRef] 6. Wang, A.; Ye, J.; Humphrey, M.G.; Zhang, C. Graphene and Carbon-Nanotube Nanohybrids Covalently Functionalized by Porphyrins and Phthalocyanines for Optoelectronic Properties. Adv. Mater. 2018,30, 1705704. [CrossRef] 7. Mustonen, K.; Susi, T.; Kaskela, A.; Laiho, P.; Tian, Y.; Nasibulin, A.G.; Kauppinen, E.I. Influence of the Diameter of Single-Walled Carbon Nanotube Bundles on the Optoelectronic Performance of Dry-Deposited Thin Films. Beilstein J. Nanotechnol. 2012 ,3, 692–702. [CrossRef] 8. Escárcega-Bobadilla, M.V.; Zelada-Guillén, G.A.; Pyrlin, S.V.; Wegrzyn, M.; Ramos, M.M.D.; Giménez, E.; Stewart, A.; Maier, G.; Kleij, A.W. Nanorings and Rods Interconnected by Self-Assembly Mimicking an Artificial Network of Neurons. Nat Commun 2013,4, 2648. [CrossRef] 9. Hong, S.W.; Jeong, W.; Ko, H.; Kessler, M.R.; Tsukruk, V.V.; Lin, Z. Directed Self-Assembly of Gradient Concentric Carbon Nanotube Rings. Adv. Funct. Mater. 2008,18, 2114–2122. [CrossRef] 10. Basu, S.; Patra, P.; Sarkar, J. Dewetting Assisted Self-Assembly of Carbon Nanotube into Circular Nanorings. Chem Eng Sci 2022 , 261, 117961. [CrossRef] 11. Koner, K.; Karak, S.; Kandambeth, S.; Karak, S.; Thomas, N.; Leanza, L.; Perego, C.; Pesce, L.; Capelli, R.; Moun, M.; et al. Porous Covalent Organic Nanotubes and Their Assembly in Loops and Toroids. Nat. Chem. 2022,14, 507–514. [CrossRef] [PubMed] 12. Datta, S.; Kato, Y.; Higashiharaguchi, S.; Aratsu, K.; Isobe, A.; Saito, T.; Prabhu, D.D.; Kitamoto, Y.; Hollamby, M.J.; Smith, A.J.; et al. Self-Assembled Poly-Catenanes from Supramolecular Toroidal Building Blocks. 400 Nat. 2020 ,583, 400–405. [CrossRef] [PubMed] 13. Pyrlin, S.V.; Hine, N.D.M.; Kleij, A.W.; Ramos, M.M.D. Self-Assembly of Bis-Salphen Compounds: From Semiflexible Chains to Webs of Nanorings. Soft Matter 2018,14, 1181–1194. [CrossRef] [PubMed] 14. Nirmalraj, P.N.; Lyons, P.E.; De, S.; Coleman, J.N.; Boland, J.J. Electrical Connectivity in Single-Walled Carbon Nanotube Networks. Nano. Lett. 2009,9, 3890–3895. [CrossRef] [PubMed] 15. Buldum, A.; Lu, J.P. Contact Resistance between Carbon Nanotubes. Phys. Rev. B 2001,63, 161403. [CrossRef] 16. Alberti, S.A.N.; Schneider, J.; Müller-Plathe, F. Mobility of Polymer Melts in a Regular Array of Carbon Nanotubes. J. Chem. Theory Comput. 2021,18, 3295. [CrossRef] 17. Skyner, R.E.; McDonagh, J.L.; Groom, C.R.; van Mourik, T.; Mitchell, J.B.O. A Review of Methods for the Calculation of Solution Free Energies and the Modelling of Systems in Solution. Phys. Chem. Chem. Phys. 2015,17, 6174–6191. [CrossRef] 18. Zwanzig, R.W. High-Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. J. Chem. Phys. 1954 ,22, 1420. [CrossRef] 19. Kästner, J. Umbrella Sampling. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011,1, 932–942. [CrossRef] 20. Virnau, P.; Muller, M. Calculation of Free Energy through Successive Umbrella Sampling. J. Chem. Phys. 2004,120, 10925. [CrossRef] 21. Kumar, S.; Rosenberg, J.M.; Bouzida, D.; Swendsen, R.H.; Kollman, P.A. THE Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992,13, 1011–1021. [CrossRef] 22. Jarzynski, C. Nonequilibrium Equality for Free Energy Differences. Phys. Rev. Lett. 1997,78, 2690–2693. [CrossRef] 23. Park, S.; Schulten, K. Calculating Potentials of Mean Force from Steered Molecular Dynamics Simulations. J. Chem. Phys. 2004 , 120, 5946–5961. [CrossRef] 24. Park, S.; Khalili-Araghi, F.; Tajkhorshid, E.; Schulten, K. Free Energy Calculation from Steered Molecular Dynamics Simulations Using Jarzynski’s Equality. J. Chem. Phys. 2003,119, 3559. [CrossRef] 25. Comer, J.; Gumbart, J.C.; Hénin, J.; Lelievre, T.; Pohorille, A.; Chipot, C. The Adaptive Biasing Force Method: Everything You Always Wanted to Know but Were Afraid to Ask. J. Phys. Chem. B 2015,119, 1129–1151. [CrossRef] [PubMed] 26. Shivakumar, D.; Williams, J.; Wu, Y.; Damm, W.; Shelley, J.; Sherman, W. Prediction of Absolute Solvation Free Energies Using Molecular Dynamics Free Energy Perturbation and the Opls Force Field. J. Chem. Theory Comput. 2010,6, 1509–1519. [CrossRef] [PubMed] 27. Pohorille, A.; Jarzynski, C.; Chipot, C. Good Practices in Free-Energy Calculations. J. Phys. Chem. B 2010 ,114, 10235–10253. [CrossRef] [PubMed] 28. Gumbart, J.C.; Roux, B.; Chipot, C. Standard Binding Free Energies from Computer Simulations: What Is the Best Strategy? J. Chem. Theory Comput. 2013,9, 794–802. [CrossRef] 29. Bennett, C.H. Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Comput. Phys. 1976,22, 245–268. [CrossRef] 30. Kim, I.; Allen, T.W. Bennett’s Acceptance Ratio and Histogram Analysis Methods Enhanced by Umbrella Sampling along a Reaction Coordinate in Configurational Space. J. Chem. Phys. 2012,136, 164103. [CrossRef] 31. Abrams, J.B.; Rosso, L.; Tuckerman, M.E. Efficient and Precise Solvation Free Energies via Alchemical Adiabatic Molecular Dynamics. J. Chem. Phys. 2006,125, 074115. [CrossRef] [PubMed] 32. Cournia, Z.; Allen, B.; Sherman, W. Relative Binding Free Energy Calculations in Drug Discovery: Recent Advances and Practical Considerations. J. Chem. Inf. Model 2017,57, 2911–2937. [CrossRef] [PubMed] 33. Williams-Noonan, B.J.; Yuriev, E.; Chalmers, D.K. Free Energy Methods in Drug Design: Prospects of “Alchemical Perturbation” in Medicinal Chemistry. J. Med. Chem. 2018,61, 638–649. [CrossRef] [PubMed] Polymers 2022,14, 4525 17 of 18 34. Mukhopadhyay, T.K.; Datta, A. Disentangling the Liquid Phase Exfoliation of Two-Dimensional Materials: An “in Silico” Perspective. Phys. Chem. Chem. Phys. 2020,22, 22157–22179. [CrossRef] [PubMed] 35. Liu, B.; Salgado, S.; Maheshwari, V.; Liu, J. DNA Adsorbed on Graphene and Graphene Oxide: Fundamental Interactions, Desorption and Applications. Curr. Opin. Colloid Interface Sci. 2016,26, 41–49. [CrossRef] 36. Ghaderi, N.; Peressi, M. First-Principle Study of Hydroxyl Functional Groups on Pristine, Defected Graphene, and Graphene Epoxide. J. Phys. Chem. C 2010,114, 21625–21630. [CrossRef] 37. Mehandzhiyski, A.Y.; Morita, M.; Oya, Y.; Kato, N.; Mori, K.; Koyanagi, J. Effect of Electrostatic Interactions on the Interfacial Energy between Thermoplastic Polymers and Graphene Oxide: A Molecular Dynamics Study. Polymers 2022 ,14, 2579. [CrossRef] 38. Soler, J.M.; Artacho, E.; Gale, J.D.; García, A.; Junquera, J.; Ordejón, P.; Sánchez-Portal, D. The SIESTA Method for Ab Initio OrderN Materials Simulation. J. Phys. Condens. Matter 2002,14, 2745–2779. [CrossRef] 39. Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996 ,77, 3865–3868. [CrossRef] 40. Junquera, J.; Paz, Ó.; Sánchez-Portal, D.; Artacho, E. Numerical Atomic Orbitals for Linear-Scaling Calculations. Phys. Rev. B 2001,64, 235111. [CrossRef] 41. Troullier, N.; Martins, J.L. Efficient Pseudopotentials for Plane-Wave Calculations. Phys. Rev. B 1991 ,43, 1993–2006. [CrossRef] [PubMed] 42. Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006,27, 1787–1799. [CrossRef] [PubMed] 43. Sordo, J.A.; Chin, S.; Sordo, T.L. On the Counterpoise Correction for the Basis Set Superposition Error in Large Systems. Theor. Chim. Acta 1988,74, 101–110. [CrossRef] 44. Thompson, A.P.; Aktulga, H.M.; Berger, R.; Bolintineanu, D.S.; Brown, W.M.; Crozier, P.S.; in’t Veld, P.J.; Kohlmeyer, A.; Moore, S.G.; Nguyen, T.D.; et al. LAMMPS—A Flexible Simulation Tool for Particle-Based Materials Modeling at the Atomic, Meso, and Continuum Scales. Comput. Phys. Commun. 2022,271, 108171. [CrossRef] 45. Luty, B.A.; Davis, M.E.; Tironi, I.G.; van Gunsteren, W.F. A Comparison of Particle-Particle, Particle-Mesh and Ewald Methods for Calculating Electrostatic Interactions in Periodic Molecular Systems. Mol. Simul. 1994,14, 11–20. [CrossRef] 46. Dupradeau, F.-Y.; Pigache, A.; Zaffran, T.; Savineau, C.; Lelong, R.; Grivel, N.; Lelong, D.; Rosanski, W.; Cieplak, P. The R.E.D. Tools: Advances in RESP and ESP Charge Derivation and Force Field Library Building. Phys. Chem. Chem. Phys. 2010 ,12, 7821. [CrossRef] [PubMed] 47. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Montgomery, J.A., Jr.; Vreven, T.; Kudin, K.N.; Burant, J.C.; et al. Gaussian 03, Revision C.02. J. Comput. Chem. 2004,24, 1748–1757. 48. Stote, R.H.; Karplus, M. Zinc Binding in Proteins and Solution: A Simple but Accurate Nonbonded Representation. Proteins Struct. Funct. Genet. 1995,23, 12–31. [CrossRef] 49. Beutler, T.C.; Mark, A.E.; van Schaik, R.C.; Gerber, P.R.; van Gunsteren, W.F. Avoiding Singularities and Numerical Instabilities in Free Energy Calculations Based on Molecular Simulations. Chem. Phys. Lett. 1994,222, 529–539. [CrossRef] 50. Jorgensen, W.L.; Maxwell, D.S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996,118, 11225–11236. [CrossRef] 51. Kaminski, G.A.; Friesner, R.A.; Tirado-Rives, J.; Jorgensen, W.L. Evaluation and Reparametrization of the OPLS-AA Force Field for Proteins via Comparison with Accurate Quantum Chemical Calculations on Peptides. J. Phys. Chem. B 2001 ,105, 6474–6487. [CrossRef] 52. Banks, J.L.; Beard, H.S.; Cao, Y.; Cho, A.E.; Damm, W.; Farid, R.; Felts, A.K.; Halgren, T.A.; Mainz, D.T.; Maple, J.R.; et al. Integrated Modeling Program, Applied Chemical Theory (IMPACT). J. Comput. Chem. 2005 ,26, 1752–1780. [CrossRef] [PubMed] 53. Robertson, M.J.; Skiniotis, G. Development of OPLS-AA/M Parameters for Simulations of G Protein-Coupled Receptors and Other Membrane Proteins. J. Chem. Theory Comput. 2022,18, 4482–4489. [CrossRef] [PubMed] 54. Doherty, B.; Zhong, X.; Gathiaka, S.; Li, B.; Acevedo, O. Revisiting OPLS Force Field Parameters for Ionic Liquid Simulations. J. Chem. Theory Comput. 2017,13, 6131–6135. [CrossRef] 55. Harder, E.; Damm, W.; Maple, J.; Wu, C.; Reboul, M.; Xiang, J.Y.; Wang, L.; Lupyan, D.; Dahlgren, M.K.; Knight, J.L.; et al. OPLS3: A Force Field Providing Broad Coverage of Drug-like Small Molecules and Proteins. J. Chem. Theory Comput. 2016 ,12, 281–296. [CrossRef] [PubMed] 56. Lu, C.; Wu, C.; Ghoreishi, D.; Chen, W.; Wang, L.; Damm, W.; Ross, G.A.; Dahlgren, M.K.; Russell, E.; von Bargen, C.D.; et al. OPLS4: Improving Force Field Accuracy on Challenging Regimes of Chemical Space. J. Chem. Theory Comput. 2021 ,17, 4291–4300. [CrossRef] [PubMed] 57. Aljuhani, A.; Ahmed, H.E.A.; Ihmaid, S.K.; Omar, A.M.; Althagfan, S.S.; Alahmadi, Y.M.; Ahmad, I.; Patel, H.; Ahmed, S.; Almikhlafi, M.A.; et al. In Vitro and Computational Investigations of Novel Synthetic Carboxamide-Linked Pyridopyrrolopyrimidines with Potent Activity as SARS-CoV-2-M Pro Inhibitors. RSC Adv. 2022,12, 26895–26907. [CrossRef] 58. Monti, M.; Stener, M.; Aschi, M. A Computational Approach for Modeling Electronic Circular Dichroism of Solvated Chromophores. J. Comput. Chem. 2022,15, 2023–2036. [CrossRef] 59. Ding, K.; Xu, C.; Peng, Z.; Long, X.; Shi, J.; Li, Z.; Zhang, Y.; Lai, J.; Chen, L.; Cai, Y.-P.; et al. Tuning the Solvent Alkyl Chain to Tailor Electrolyte Solvation for Stable Li-Metal Batteries. ACS Appl. Mater. Interfaces 2022,14, 44470–44478. [CrossRef] Polymers 2022,14, 4525 18 of 18 60. Jewett, A.I.; Stelter, D.; Lambert, J.; Saladi, S.M.; Roscioni, O.M.; Ricci, M.; Autin, L.; Maritan, M.; Bashusqeh, S.M.; Keyes, T.; et al. Moltemplate: A Tool for Coarse-Grained Modeling of Complex Biological Matter and Soft Condensed Matter Physics. J. Mol. Biol. 2021,433, 166841. [CrossRef] 61. Holian, B.L.; Voter, A.F.; Ravelo, R. Thermostatted Molecular Dynamics: How to Avoid the Toda Demon Hidden in Nose-Hoover Dynamics. Phys. Rev. E 1995,52, 2338. [CrossRef] [PubMed] 62. Schneider, T.; Stoll, E. Molecular-Dynamics Study of a Three-Dimensional One-Component Model for Distortive Phase Transitions. Phys. Rev. B 1978,17, 1302. [CrossRef] 63. Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph. 1996,14, 33–38. [CrossRef] 64. Steinbrecher, T.; Mobley, D.L.; Case, D.A. Nonlinear Scaling Schemes for Lennard-Jones Interactions in Free Energy Calculations. J. Chem. Phys. 2007,127, 214108. [CrossRef] [PubMed] 65. Shirts, M.R.; Pande, V.S. Comparison of Efficiency and Bias of Free Energies Computed by Exponential Averaging, the Bennett Acceptance Ratio, and Thermodynamic Integration. J. Chem. Phys. 2005,122, 144107. [CrossRef] [PubMed] 66. Klimovich, P.V.; Shirts, M.R.; Mobley, D.L. Guidelines for the Analysis of Free Energy Calculations. J. Comput. Aided Mol. Des. 2015,29, 397–411. [CrossRef] 67. Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R.A.; et al. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009,180, 1961–1972. [CrossRef] 68. Kharisov, B.I.; Kharissova, O.V.; Dimas, A.V. The Dispersion, Solubilization and Stabilization in “Solution” of Single-Walled Carbon Nanotubes. RSC Adv. 2016,6, 68760–68787. [CrossRef] 69. Cao, X.Z.; Merlitz, H.; Wu, C.X.; Ungar, G.; Sommer, J.U. A Theoretical Study of Dispersion-to-Aggregation of Nanoparticles in Adsorbing Polymers Using Molecular Dynamics Simulations. Nanoscale 2016,8, 6964–6968. [CrossRef] 70. Elkashef, M.; Wang, K.; Abou-Zeid, M.N. Acid-Treated Carbon Nanotubes and Their Effects on Mortar Strength. Front. Struct. Civ. Eng. 2015,10, 180–188. [CrossRef] 71. Liang, S.; Li, G.; Tian, R. Multi-Walled Carbon Nanotubes Functionalized with a Ultrahigh Fraction of Carboxyl and Hydroxyl Groups by Ultrasound-Assisted Oxidation. J. Mater. Sci. 2016,51, 3513–3524. [CrossRef] 72. Datsyuk, V.; Kalyva, M.; Papagelis, K.; Parthenios, J.; Tasis, D.; Siokou, A.; Kallitsis, I.; Galiotis, C. Chemical Oxidation of Multiwalled Carbon Nanotubes. Carbon N. Y. 2008,46, 833–840. [CrossRef]