Computational study of the water-driven graphene wrinkle life-cycle towards applications in flexible electronics

The ubiquitous presence of wrinkles in two-dimensional materials alters their properties significantly. It is observed that during the growth process of graphene, water molecules, sourced from ambient humidity or transferred method used, can get diffused in between graphene and the substrate. The water diffusion causes/assists wrinkle formation in graphene, which influences its properties. The diffused water eventually dries, altering the geometrical parameters and properties of wrinkled graphene nanoribbons. Our study reveals that the initially distributed wrinkles tend to coalesce to form a localized wrinkle whose configuration depends on the initial wrinkle geometry and the quantity of the diffused water. The movement of the localized wrinkle is categorized into three modes—bending, buckling, and sliding. The sliding mode is characterized in terms of velocity as a function of diffused water quantity. Direct bandgap increases linearly with the initial angle except the highest angle considered (21°), which can be attributed to the electron tunneling effect observed in the orbital analysis. The system becomes stable with an increase in the initial angle of wrinkle as observed from the potential energy plots extracted from MD trajectories and confirmed with the DOS plot. The maximum stress generated is less than the plastic limit of the graphene.

. (a) Schematic diagram of the wrinkle formation process through the ice-like water formation and liquid water diffusion at high humidity, and liquid water evaporation in a dry environment. (b) Height changes in a graphene/SiO 2 structure after exposure to high humidity. AFM topographic images of a monolayer graphene mechanically exfoliated on SiO 2 (b.1) before and (b.2) after the high humidity exposure. (c) Averaged height profiles measured for the areas marked with black and red rectangle in (b.1) and (b.2), respectively. G denotes monolayer graphene. © Copyright : Springer Science + Business Media, Lee et al. 40 .
Scientific RepoRtS | (2020) 10:11315 | https://doi.org/10.1038/s41598-020-68080-5 www.nature.com/scientificreports/ in geometry of the wrinkles, respectively, resulting in Wrinkled Graphene Nanoribbon (WGNR). It is essential to study wrinkle transformation lead by water diffusion and evaporation in the system. In this work, we provided a detailed explanation of water diffusion and its impact on a wrinkled bilayer graphene system characterized by geometrical parameters and electronic properties. Our results can explain the phenomenon of the water interaction with the graphene systems 36,37 , and change in electronic properties of different 2D materials caused by wrinkles 38,40 . Our study is composed of three phases. Phase I addresses the wrinkle evolution due to diffused water in between graphene with distributed wrinkles and a graphene substrate. We had varied the amount of diffused water and the initial angle of the distributed wrinkles in graphene and studied the formation of final localized wrinkle configuration and the stress generated in wrinkled graphene. Phase II deals with the evaporation of diffused water, and the effect of evaporation on the collapsing of the wrinkle, i.e. completion of wrinkle life-cycle. We observed that, upon complete evaporation, the resulting localized wrinkle undergoes motion that is composed of three fundamental modes-bending, buckling, and sliding as studied theoretically by Im et al. 41 Lastly, Phase III uses Density Functional Theory (DFT) while considering the free-standing WGNR with four different initial angles and studies its influence on the electronic properties. Only quasi-periodic wrinkles are considered instead of a complex network of wrinkles that exist in the experimental setups. However, we believe that our studies will encourage researchers to study more complicated wrinkle structures.

Molecular dynamics (MD).
This study is divided into three phases based on which part of the wrinkle lifecycle is being studied, i.e., formation, evolution, collapse, and the approach used for analysis, i.e., MD or DFT. In the first phase, we investigated the wrinkle formation and evolution. In this phase, as shown in Fig. 2a, using the in-house MATLAB code, we created a Graphene with Distributed Wrinkles (GDW) supported by a Flat Graphene (FG) with Periodic Boundary Conditions (PBC) (width = 20 Å ). Here, we considered different Initial Angle of Wrinkle (IAW) denoted as θ IAW ( θ IAW = 6°, 11°, 16°, 21°) to generate GDW with four different configurations. Figure 2a shows a representative case of θ IAW = 11 • . While changing θ IAW , we kept the number of carbon (C) atoms in GDW the same, but the size of the FG is reduced. For example, if we take a flat graphene sheet and keep compressing to get wrinkles of different angles, the number of atoms in the sheet remains unchanged, but the overall length of the sheet is reduced. Table S1 of the Supplementary Information (SI) shows the C atoms for different θ IAW . We define a parameter called Carbon Ratio (CR), which is the ratio of the number of C atoms in GDW to that of FG. For θ IAW of 6°, 11°, 16°, 21°; CR values are 1.017, 1.053, 1.11, 1.20, respectively (Table S1).
The Adaptive Intermolecular Reactive Bond Order (AIREBO) potential was used for the C-C interactions 42 . We performed an energy minimization of this system (Fig. 2a) using the conjugate gradient method, as www.nature.com/scientificreports/ implemented in the MD package: LAMMPS 43 . Figure 2b shows the minimized/equilibrium configuration of GDW supported by the FG. At this point, we studied wrinkling evolution for four cases-(a) without water (Fig. 2b), and (b) two-, (c) four-, and (d) six-layer of water in between the GDW and FG (Fig. 2c-e). In-house MATLAB code was used to insert the water layer in the wrinkled system shown in Fig. 2b. The density of the added water molecules, i.e., the number of water molecules per unit volume, is the same for all cases. Table S2 of SI shows the number of water molecules for different cases. The water molecules were added to the system (Fig. 2b) in two different ways-(1) While adding the water, the gap between the graphene layers was increased to maintain a fixed distance between the water layers and the adjacent graphene sheet (Figs. 2c-e). (2) Water layers were added while keeping the distance between GDW and FG to a fixed distance corresponding to the maximum water case considered (i.e., gap shown in Fig. 2e is considered for all other cases). However, in both cases, we qualitatively observed the same results. Therefore, we only discuss the results for the cases shown in Figs. 2b-e. We used TIP4P potential for water (H 2 O) 44 , and Lennard-Jones type of pair potential for C and H 2 O interactions 45 . We used TIP4P potential before for water modeling 3,37 , and our computational results were in good agreement with the experimental findings 3 .
The methodology of Phase I is summarized in the 'Wrinkle Evolution' part of the flow-chart shown in Fig. 3. The temperature of the system was controlled at 300 K using the Berendsen thermostat 46 . We then performed MD simulation using NVE ensemble with the timestep of 1 femtosecond. Sufficient MD steps (at least one million steps) were performed to make sure the system reached its lowest energy configurations (discussed later in detail). After the end of the first phase simulations, we stored the system configuration and reused it as the initial structure for the next phase of the project, i.e., water drying case. We performed the stress analysis of the final stabilized graphene sheet with a localized wrinkle using the in-built stress computation algorithm of LAMMPS 47 . The detail of the stress calculation methods is provided in the supplementary information (section D).
In Phase II, we studied how wrinkle collapses due to the drying of water. The 'Evaporation' part of Fig. 3 summarizes the methodology for studying the drying phase. For drying modeling, we started with the last configuration of the system in the first phase, e.g., schematic in Fig. 3b1. In our LAMMPS script, the water molecules at the outer periphery were identified. The identification process was accomplished by locating and then grouping those water molecules having lesser coordination numbers as compared to their counterparts inside the waterbody. The identified surface water molecules were deleted, followed by running the MD simulation for 2000 steps. The same loop was performed again (water deletion and MD simulation) until all water molecules were deleted. www.nature.com/scientificreports/ Density functional theory (DFT). In Phase III, the DFT based analyses were performed on one frame from each case considered from Phase I, as shown in Figs. 11, 12, 13. In this phase, we tried to understand the impact of geometrical parameters of WGNR on their fundamental electronic properties. We implemented first-principles DFT with plane-wave basis sets and pseudopotentials to describe the electron-ion interactions as implemented in Quantum Espresso (QE) 48 . All calculations were done using the projector augmented wave (PAW) pseudopotentials and the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 49 . The planewave basis sets were used with a plane-wave cutoff energy of 100 eV, and a kinetic energy cutoff for augmentation charges of 400 eV. The convergence threshold of Kohn-Sham equations was set to 1e−06. The gammacentered k-point sampling grids obtained using the Monkhorst-Pack method, were 8 × 8 × 8 with a unit offset for the graphene unit cell. The valence electrons contain s and p orbitals for carbon. Before DFT calculations, all atoms in the cell, as well as the lattice dimensions and angles, were relaxed to the equilibrium configurations by using MD. For consistency, only one case corresponding to two water (H 2 O) layers was studied from each considered IAW, i.e., 6°, 11°, 16°, 21° (Fig. 11). Furthermore, the systems from all four cases were trimmed down to 120 atoms to eliminate the system size dependency, and conform to the available computational resources. For the band structure calculations, the symbols and coordinates of the high-symmetry points in the first Brillouin zone of the crystals were taken from Y Hinuma et al. 50 . MATLAB and VESTA codes were used for the postprocessing of the results.

Results and discussions
We discuss the results for three phases-(I) The formation and evolution of wrinkles for no-water, and water inside GDW and supported FG (without drying), (II) The collapse of the localized wrinkle during water drying, and (III) The electronic structure of free-standing WGNR for four different geometries.
Phase I: The formation and evolution of wrinkles due to diffused water. As discussed in "Introduction" section, during the synthesis 51 and/or transfer process 52 , graphene can inherent distributed wrinkles, and they tend to coalescence to form a localized structures, i.e., from 'Graphene with Distributed Wrinkles (GDW)' to 'Graphene with Localized Wrinkle (GLW)' . For the simplicity of analysis, we considered here GDW only among other available morphologies (Fig. 2b), i.e., ripples, crumples, folds, etc. First, we analyzed how GDW transforms into GLW when there is no diffused water (Fig. 2b). Figure 4a shows the case of distributed wrinkles for the representative case of θ IAW = 21 • and CR = 1.21. As mentioned earlier in "Models and methods" section, Carbon Ratio (CR) is the ratio of C atoms in the upper sheet with wrinkles and the underlying flat sheet (see Table S1 of SI). The video related to this phenomenon is shown in the SI (Movie_Figure4a.mp4). The distributed wrinkles on same-side amalgamate to form a localized wrinkle, i.e., GDW transforms to GLW. In crystalline solids, same-sign dislocations repel according to Frank's rule 53 . However, in van der Waals (vdW) www.nature.com/scientificreports/ layered structures, such as graphene, same-sign (same-side) winkles attract each other to form a localized wrinkle. A similar observation was reported for the MoS 2 bilayer system 54 . As shown in Fig. 4a, for θ IAW = 21 • , first, the initially distributed wrinkles transform into intermediate wrinkle localization and finally single 'localized wrinkle' . During the intermediate steps (step = 10,000, 25,000), initially distributed wrinkles form three followed by two wrinkles and, finally, one 'localized wrinkle' . There is a drop in potential energy during the transformation of three to two and, finally, one wrinkle. In SI, Figures S1, S2, S3 show the evolution of initial wrinkle for three other cases-(1) θ IAW = 6 • , CR = 1.017; (2) θ IAW = 11 • , CR = 1.053; (3) θ IAW = 16 • , CR = 1.11. As shown in Figure S1, when the θ IAW is less (e.g., θ IAW = 6 • ), and CR is close to 1, initially distributed wrinkles do not converge to single localized wrinkle, but tries to transform into a highly strained flat sheet. Figure 4b represents a case-study where water molecules are initially placed on monolayer FG with no graphene sheet on top. The SI contains the related video (Movie_Figure4b.mp4). The MD simulation shows that water molecules coalescence to form a water droplet because of the hydrophobicity of graphene 55 with a contact angle of (≈)125°, which matches with the literature 56 . In Fig. 4c, we considered a representative case of θ IAW = 6 • (two-layer water) to discuss how diffused water molecules inside GDW and FG influence the wrinkle evolution. The SI contains the related video (Movie_Figure4c.mp4). The presence of diffused water causes 'competition' between two phenomena-(i) wrinkle amalgamation, and (ii) droplet formation. The final localized wrinkle configuration is primarily determined by the dominant phenomenon. Figure S1 shows that for GDW with θ IAW = 6 • , initial wrinkles tend to merge together to form a flat sheet. Figure 4b shows the water droplet formation starting from the initial water layer. When these phenomena, i.e., wrinkling formation and droplet formation interact with each other, the droplet formation process dominates. However, the final droplet shape in this case (Fig. 4c) is not the same as in Fig. 4b. Because in this case, the upper graphene sheet tries to compress the water droplet, as observed by McKenzie et al. 57 . However, since the droplet formation phenomenon is dominant here, the final wrinkle configuration is primarily determined by the water droplet. Figures S4, S5, S6, S7 show the final wrinkle configuration of graphene starting from the initial wrinkle angle of θ IAW = 6 • , 11 • , 16 • , and 21 • respectively for 0-, 2-, 4-, and 6-layer water cases. The stress generated during the wrinkle formation is shown here. We considered a representative case of θ IAW = 11 • for the 2-layer water case ( Figure S6b) for analyzing the stress generation during the wrinkle formation. Figure 5 shows the z-component of the Cauchy stress. Considering the stress-field in Fig. 5a, we note that atoms in the flat regions and at the top of the wrinkle have no z-component stress. Only the sidewalls, i.e., the base of the localized wrinkle, are stressed. Here, stresses are averaged upon the whole line of atoms along the periodic y-axis. We note that for all three stress cases ( σ zz , σ xz , σ yz ), the stress is high in the sidewalls (the base) of the wrinkle as expected 58 . The maximum individual atomic stress, observed in the wrinkle, varies from around − 17 GPa (compressive) to 18 GPa (tensile) (Fig. 5a). When averaged over the stripe of atoms along the periodic y-axis, the maximum magnitude of stress is around 4 GPa (Fig. 5b). In summary, stress generation in graphene during the wrinkle formation is much lower than the stress that can cause plastic deformation 59,60 . Hence, water diffusion-induced wrinkle formation in graphene does not cause in-plane lattice-strain high enough to break the bonds and cause plasticity. As shown in Figures S5, S6, S7,S8 irrespective of the number of the diffused water layers (i.e., 2-, 4-, or 6-layer), the pattern of the final configuration is the same for all water layer cases. We define two new parameters, as shown in Fig. 6a-(1) Aspect Ratio (AR) as the ratio of the height of the wrinkle to its base length, and (2) 'Final Angle of Wrinkle, θ FAW ' . If AR = 0, it means the sheet is flat (height = 0). The higher the value of AR, the lesser the radius of the curvature of wrinkle-vertex. Moreover, there is a directly proportional relationship between AR and θ FAW . Figure 6b,c show how AR and θ FAW vary with θ IAW (initial angle of water) for the differing number of the water layers. When there is no diffused water ( Figures S5a, S6a, S7a, S8a), as θ IAW increases from 6 • to 21 • , the final localized wrinkle becomes sharper. Therefore, the blue curve (no water case) in Fig. 6b shows that the AR and θ FAW increase as θ IAW increases. As mentioned before, the presence of diffused water starts the complex interplay between the wrinkle amalgamation and droplet formation 61 . The final shape is the result of the interaction between two phenomena. For 2-, 4-, or 6-layer diffused water cases ( Figures S5b-d, S6b-d, S7b-d, S8b-d), AR and θ FAW linearly increases as θ IAW increases from 6 • to 21 • . As shown in Figure S5a, for θ IAW = 6 • case, the graphene sheet with initially distributed wrinkle flattens (AR = 0, θ FAW = 0 • in Fig. 6b.1) for no water case. With the insertion of 2-, 4-, or 6-layer diffused water ( Figure S5b-d), the droplet and wrinkle size increases. However, the decrease in AR for a different water-layer case for θ IAW = 6 • is not significant (black line in Fig. 6c.1) as the length and height of the wrinkle increase almost proportionately. For the θ IAW = 21 • case ( Figure S8), we have the spiker localized wrinkle for no-water case ( Figure S8a). The AR, in this case, is high (blue line in Fig. 6b.1), since the base length of wrinkle is moderately short as compared to the height of the wrinkle. With the increase in diffused water amount, the base length of the wrinkle keeps increasing while the height is almost the same ( Figure S8b-d). Thus, AR and θ FAW values drop significantly (cyan curve in Fig. 6c).
In Phase I, it is observed that for all cases considered, distributed wrinkles evolve into a localized wrinkle. The way it appears is either the intra-attraction between water molecules to make one droplet, and the vdW attraction between the wrinkle and water molecules that "drag" the graphene along with water droplet 62 . Another possibility is the vdW attraction between the upper graphene sheet and graphene substrate, causing that drag. www.nature.com/scientificreports/ The probability of the vdW attraction between walls of the same wrinkle, during evolution, are slim to none given the wrinkle wavelength. Furthermore, it may be a combination of all these factors, which may be the function of temporal and spatial variables as well. Though the detailed analysis is beyond the scope of this study, in Phase I, we tried to get an insight into the underlying mechanisms into the wrinkling formation and evolution kinetics.
Phase II: The dynamics of the localized wrinkle due to water evaporation/drying. As shown in Fig. 1a, the diffused water evaporates, which changes the structure (AR and θ FAW ) of the localized wrinkle. We consider the final localized wrinkle structure for θ IAW = 11 • case with 2-layer water diffusion ( Figure S6b) and study the variation in wrinkle configurations upon evaporation. In Fig. 7a, we note the sequence of steps of water evaporation and the change in wrinkle configuration. The SI contains the related video (Movie_Figure7a.mp4). Figure 7b shows the corresponding change in potential energy. We performed the comparative analysis in this plot (as any other potential energy plot in this work), with the least interest in absolute values. When water starts evaporating, the potential energy drops drastically as there is a significant perturbation in the system. Upon drying, there is an approximately 25% decrease in the angle of curvature of wrinkle, which is analogous to the experimental observation shown in Fig. 1. Similarly, we studied the water evaporation and wrinkle formation for other initial angle cases ( θ IAW = 6 • , 16 • , 21 • ). Figure S8 shows the wrinkle for 6-layer diffused water case for θ IAW = 21 • . The SI contains the related video (Movie_FigureS8.mp4). As expected, there was an inverse proportionality in the angle of curvature (Table S3). Here, angle of curvature decreases as opposed to what was observed for other cases, i.e. ( θ IAW = 6 • , 11 • , 16 • ). That is attributed to the high slenderness/aspect ratio because of the highest CR among cases considered. We had performed the angle of curvature analysis for all θ IAW considered, as shown in Table S3. After the complete evaporation, the localized wrinkle undergoes structural changes, which results in three kinds of fundamental movements-(1) Buckling Mode (Fig. 8a): The two sides of the wrinkle buckles and reduces the base-length of the wrinkle. However, the height of the wrinkle remains the same (AR increases). (2) Bending Mode (Fig. 8b): The base length and height of the wrinkle remain the same (AR = constant). However, the wrinkle bends sideways. (3) Sliding Mode (Fig. 8c): In this case, after the complete evaporation, the final localized wrinkle does not remain static at a particular location. Instead, it keeps moving left or right. The sliding mode is more or less present in all cases beyond θ IAW = 6 • , i.e., after complete evaporation, wrinkle moves sideways. The bending and buckling modes dominate as the initial angle increases. Figure S9 of SI shows the evaporation for a 6-layer water diffusion case. Soon after the complete drying, we notice the buckling at the side of wrinkle. The buckling mode is accompanied by sliding and bending mode (See the related video in SI-Movie_FigureS8. mp4). The bending mode is less dominant for low θ IAW cases ( 6 • , 11 • ) when the height of the wrinkle is lesser.
A critical study is needed to analyze which mode (sliding, bending, buckling) is more dominant at different stages of evaporation, given various initial angles, to tune the wrinkle according to the needs as targeted by other methods 63 . We attempted to get insights into the domination of sliding mode during drying. As mentioned before, for a lower initial angle (e.g., θ IAW = 11 • ), the sliding mode is more prominent. Figure 9 shows the velocity of the wrinkle as a function of water content. The 'wrinkle-velocity' is computed with reference to the atom present at the crest of the wrinkle. The identity of the atom located at the crest is changing as wrinkle is moving under the sliding mode. We utilized the 'in-built' python-scripting facility in visualization package OVITO 64 to identify the atomic position at the crest for a different time frame and monitor its movement to compute the velocity. In Figure S10 of the SI, we plotted the number of water molecules present and the corresponding 'wrinkle-velocity' . Because of the atomic scale movement, there is a lot of noise in the velocity data. Here, positive data represents velocity when wrinkle moves towards the positive x-axis, while negative represents when wrinkle moves along the negative x-axis. In Fig. 9, the velocity is averaged over every 50,000 steps. We note that as long as water molecules are diffused inside the wrinkle, the speed ranges around zero (until 0.5 × 1e5 MD steps). However, once the water molecules are obliterated, then not only velocity jumps to a higher value of magnitude, but it also shows a www.nature.com/scientificreports/ positive gradient with MD time steps, and the wrinkle prefers to move towards the positive x-axis. This 'wrinkle-velocity' pattern indicates that the water molecules diffused under the wrinkle behaves as 'wrinkleanchor' (speed breaker). After the water molecules are completely evaporated, the sideways movement of the wrinkle dominates because the 'wrinkle-anchor' is removed. We hypothesize that the direction of the wrinkle movement after the complete evaporation is decided by the complex interplay of three different modes at the moment when the "last-set" of water-molecules is evaporated. However, further critical analysis is needed to validate this proposition.
To understand further the 'wrinkle-anchor' phenomenon by the diffused water molecules inside the wrinkle, we considered another set-up, as shown in Fig. 10. We artificially created a scenario where two wrinkles are placed apart, and the different number of water molecules are diffused inside. Different snap-shots are taken until the  www.nature.com/scientificreports/ system has reached the local minima. SI contains the related video (Movie_Figure10a.mp4). Wrinkle-A and Wrinkle-B do not amalgamate until one of them lost all the water molecules diffused underneath. As shown in Fig. 10, until time = 18,000 fs, the vdW interaction between the diffused water and the substrate acts as an 'anchor support' to provide the spatial rigidity to the wrinkle. At time = 19,000 fs, Wrinkle-A has the complete evaporation, while for Wrinkle-B, some diffused water molecules are still present. The diffused water molecules act as 'wrinkle-anchor' (speed breaker) for Wrinkle-B. Moreover, the presence of water molecules inside Wrinkle-B also prevents the sliding mode of Wrinkle-A. In fact, instead of sliding mode, here Wrinkle-A starts collapsing after complete evaporation since Wrinkle-B prevents its sideway motion. After complete evaporation for Wrinkle-B (time = 20,000 fs), Wrinkle-A collapses completely and merges with Wrinkle-B. Therefore, unlike single wrinkle cases (Fig. 7), when two wrinkles are present, the motion of one wrinkle not only depends on the amount of water diffused underneath it but also on the water underneath the adjacent wrinkle. We hypothesize that in a network of wrinkles in a 2D material, we can control the evolution of the wrinkle-network by controlling the  www.nature.com/scientificreports/ amount of water content trapped under fewer wrinkles than the total number of wrinkles participating in the wrinkle-network. It will be valuable observation 1 if we can control the evolution of the whole wrinkle-network by controlling the water content trapped under just one wrinkle, which may be or may not be the member of the original wrinkle-network.

Phase III: DFT study of the electronic properties of winkled graphene. Crystal orbitals (HOMO and LUMO).
We performed the Molecular Orbital (MO) analysis of WGNR for all four angles ( 6 • , 11 • , 16 • , 21 • ). We considered the localized wrinkle structure for four-layer of water cases shown in Figures S5c, S6c, S7c and S8c. Then we performed the complete evaporation. The resulting localized wrinkle after complete evaporation is considered here, excluding the underlying substrate, i.e., free-standing WGNR (Fig. 11). We considered the same number of atoms in all the cases, i.e., 120 atoms. The system is optimized with the MD code before performing electronic calculations with DFT. Figure 11 shows the isometric view and Fig. S12 of SI shows the front view. Crystal orbitals are demonstrated by the iso-surfaces with yellow color. A carbon atom has six electrons, out of which two lowest electrons are considered as core electrons, and the next four electrons are regarded as valence electrons. These valence electrons participate in the interaction with other atoms in the crystal. Consequently, we are considering these four electrons in the outmost shell in the DFT calculations. We observed the Highest Occupied Molecular Orbital (HOMO) at the portion of the WGNR, where the the radius of curvature is high, i.e., the flat portion (Fig. 11). For IAW of 6°, 11°, & 16°, HOMO is observed to be only on the flat part of the GNR. This indicates that lattice distortion, present due to the wrinkle, is increasing the energy level of the molecular orbitals forming up the wrinkled portion of GNR. Hence, the electrons prefer to fill up the molecular orbitals formed by the flat portion of the GNR. The same phenomenon is observed in the sidewalls of the wrinkle for the highest angle considered (21°). In this case, the wrinkle is slender enough so that its sidewalls act like flat graphene and form the molecular orbitals having less energy. Consequently, electrons prefer to stay in the molecular orbitals formed by the sidewalls of the wrinkle instead of the vertex.
The Lowest Unoccupied Molecular Orbital (LUMO) is mostly observed in the portion of the GNR having high curvature (Fig. 11). From IAW of 6° to 16°, LUMO can be found in the wrinkled section of the GNR but definitely not in the flat portion. It indicates that if the GNR gains an electron, then it will be "absorbed" by the atoms forming the wrinkle. The effect of curvature plays a dominant role when we increase the IAW beyond 16°. In the case of 21°, we observe LUMO in the form of tunnel shape 65 . The wrinkle is slender enough to confine the lowest unoccupied molecular orbital parallel to its direction. Further investigation is necessary to understand whether it's just the effect of curvature or the presence of HOMO along the sidewalls of the wrinkle or both of these factors contributing to some extent to the formation of "electron-tunnel"-LUMO. This observation proves that not only the presence of wrinkles but its slenderness also plays a crucial role in the generation of "tunnel" LUMO orbital. It can be utilized to control the direction of electron flow in a (semi)conducting system, as observed by Sulpizio et al. in the form of electron flow following a macro hydrodynamic flow pattern, i.e., Poiseuille flow 66 .
Similar to these 'electron-tunnels' , 'electron-channels' are obtained in the IAW values of 6°, 11°, 16°. The 'electron-channels' is the overlapping of the atomic orbitals of the same energy level forming an "electron-flowchannel" along the direction of wrinkle, which can act as a tunnel to the flow of every extra electron added to the system. Consequently, by controlling the course of these channels, the direction of electron-flow can also be controlled, which has multiple applications in the semiconductor industries. The studies aimed for similar outcomes are available in the literature 66,67 . For IAW of 16°, the asymmetricity is observed for LUMO. As shown in Fig. 11c (SI Figure S12c), the wrinkle is bent on one side (left). That causes the LUMO to only form on the curved sidewall completely and partially on the right sidewall. There is a portion of the right sidewall of the wrinkle, where no LUMO or HOMO is observed (SI Figure S12c). Further investigation is necessary to determine the energy level of the orbitals, possessed by the corresponding set of atoms forming the portion of the right sidewall of the wrinkle, where neither LUMO nor HOMO is obtained.
Band structure. In Fig. 12, band structures are plotted for localized wrinkled graphene structures shown in Fig. 11. Each band structure is based upon the system consisting of 240 bands corresponding to 480 electrons constituting the crystal. The direct bandgap is observed for all of the cases. The observed bandgap values, based upon the difference between HOMO and LUMO, are 0.5492, 0.5453, 0.4801, and 0.5266 eV for the localized wrinkled structure corresponding to the IAW cases of 6°, 11°, 16°, 21°, respectively. Figure S11 in SI shows the band structure of the flat pristine graphene. As expected, no bandgap is observed for this case.
The bandgap is observed to be decreasing linearly with IAW (6°, 11°, 16°), but then it increases suddenly for 21°. This can be attributed to the observation made in Fig. 11d for the 21° case, i.e., distortion of the HOMO-LUMO. As explained before, because of slenderness, the HOMO is observed on the sidewalls of the wrinkle. Moreover, because of high curvature, LUMO shifted from the tip to under the wrinkle to form an "electron tunnel". The delocalization of the HOMO and LUMO orbital for the highest angle case (21°) may be the reason behind non-linearity in the proportional relationship of IAW and bandgap.
Density of states (DOS). In Fig. 13, the density of states (DOS) is plotted for four wrinkled graphene structures mentioned in Fig. 11. An inversely proportional relationship is observed between Fermi energy and the IAW.  www.nature.com/scientificreports/ Consequently, the strained lattice increases the energy levels of the localized crystal orbitals, which results in changing the Fermi energy as a function of IAW. In Fig. 13, we notice two prominent peaks in DOS plots for all cases (IAW = 11°, 16°, 21°) except for the least IAW case (6°). We have the second peak in all the cases at around − 2 eV to 0 eV, and the first peak at around − 5 eV except for the least IAW case (6°). Further work can be conducted to find out the first value of IAW, which will generate the first peak in between IAW of 6° to 11° and results in two peaks.
As can be seen from potential energy plots ( Figure S1−S4), the local minima of potential energy increase (− 7.427 eV < − 7.426 eV < − 7.424 eV < − 7.423 eV) as the initial angle of wrinkle increases (6° < 11° < 16° < 21°). And the corresponding relationship can be observed in band structure plots while observing the Valence Band Maximum (VBM) in Fig. 12. As the initial angle of the wrinkle is increased, so does the VBM from approximately 0-1.5 eV. Synchronously, in total DOS plot in Fig. 13, the energy corresponding to the highest energy state is also increased from ≈0 eV to ≈2 eV.

conclusions
In conclusion, we have performed MD and DFT simulations to geometrically and electronically characterize the wrinkle formation, evolution, and collapse in graphene when water is diffused in between graphene and substrate, and while it is evaporating/had evaporated. Our key findings are summarized below: • The distributed wrinkles in graphene coalesce to form a 'localized wrinkle' , whose configuration significantly depends on the initial pattern of distributed wrinkles and the amount of diffused water present in between the graphene and substrate. Upon evaporation of diffused water, the localized wrinkles, except with the highest initial angle, lose their angle of curvature values in the range from 7 to 63%. • When water is distributed on a graphene sheet, due to hydrophobicity, distributed water molecules condense to form a water droplet. In the presence of diffused water, the final wrinkle configuration is the result of the competition between the water droplet formation process and the localized wrinkle formation process. For the lower initial angle of wrinkle, droplet formation is dominant, resulting in wrinkle with a high angle of curvature values, i.e., non-sharp vertex. For higher initial wrinkle angle, wrinkle formation dominates, and the resulted localized wrinkle tends to possess a lower angle of curvature, i.e., sharp vertex. • The drying of water, diffused under the localized wrinkle, changes its configuration, e.g., curvature/vertex angle. The wrinkle is static until the complete evaporation, and after that, it acquires a motion. In case of diffused water under two localized wrinkles placed side-by-side, even after the complete evaporation under one wrinkle, it remains static until the water diffused under the neighboring wrinkle is completely dried. • The wrinkle movement is the combination of three fundamental modes-bending, buckling, and sliding.
The dominant mode depends on the configuration of the wrinkle., e.g., relatively narrowed or wrinkles with relatively a lower angle of curvature tends to bend. Further investigation is necessary to understand how three different modes influence wrinkle dynamics. • The stress analysis reveals that the maximum stress is at the base of the wrinkle, and always below its plasticity limit. Therefore, wrinkle formation due to water diffusion does not cause any significant bond-breaking. • "Electron-tunnel/channel" is observed along the vertex of the WGNR. By controlling the direction of WGNR, we can control the direction of the electron flow. This concept can be used in designing flexible electronics and optoelectronics. • The direct bandgap is observed for all the wrinkle cases considered. The change in wrinkle geometry also generates lattice strains, which further changes the energy level of the crystal orbitals. The strain in lattice changes the Fermi energies which results in Fermi level modulations.
Our study addresses the fundamental questions regarding the wrinkle formation, evolution, and collapse due to water diffusion/evaporation, and its importance in flexible electronics.

Data availability
The data reported in this paper is available from the corresponding author upon reasonable request.

Code availability
The pre-and post-processing codes used in this paper are available from the corresponding author upon reasonable request. Restrictions apply to the availability of the simulation codes, which were used under license for this study.