Aurélien Doitrand
email: aurelien.doitrand@grenoble-inp.fr

Rafael, Estevez

Dominique Leguillon
email: dominique.leguillon@upmc.fr

Comparison between cohesive zone and coupled criterion modeling of crack initiation in rhombus hole specimens under quasi-static compression

Keywords: Cohesive zone, Coupled criterion, Crack initiation

Comparison between cohesive zone and coupled criterion modeling of crack initiation in rhombus hole specimens under quasi-static compression

Introduction

In a industrial design process, crack nucleation and growth determination based on predictive modeling tools appears to be a main challenge so as to reduce costly and time consuming experiments. Classical linear elastic fracture mechanics allows the growth of a pre-existing crack to be studied but fails to predict crack nucleation. To this aim, approaches based on Finite Fracture Mechanics (FFM) were developed so as to study finite crack increments [START_REF] Hashin | Finite thermoelastic fracture criterion with application to laminate cracking analysis[END_REF][START_REF] Nairn | Exact and variational theorems for fracture mechanics of composites with residual stresses, traction-loaded cracks, and imperfect interfaces[END_REF]. In this context, the abrupt nucleation a crack over a finite length can be modeled by the so-called Coupled Criterion (CC) initially proposed by Leguillon [20]. This approach is applicable for crack initiation prediction in various configurations (V-notch [START_REF] Cornetti | Short cracks and V-notches: Finite Fracture Mechanics vs[END_REF][START_REF] Leguillon | Size effects for crack initiation at blunt notches or cavities[END_REF][START_REF] Leguillon | Failure initiation at V-notch tips in quasi-brittle materials[END_REF], interfaces [START_REF] Mantič | Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion[END_REF][START_REF] Mantič | Crack onset and growth at the fibre-matrix interface under a remote biaxial transverse load. Application of a coupled stress and energy criterion[END_REF][START_REF] Martin | Competition between deflection and penetration at an interface in the vicinity of a main crack[END_REF][START_REF] Stein | A model for brittle failure in adhesive lap joints of arbitrary joint configuration[END_REF], etc) and bulk materials (composites [START_REF] Doitrand | Damage onset modeling in woven composites based on a coupled stress and energy criterion[END_REF][START_REF] Doitrand | Mesoscale analysis of damage growth in woven composites Compos Part A[END_REF][START_REF] García | A numerical study of transverse cracking in cross-ply laminates by 3D finite fracture mechanics[END_REF][START_REF] Martin | A coupled strength and toughness criterion for the prediction of the open hole tensile strength of a composite plate[END_REF], ceramics [START_REF] Leguillon | Application of the coupled stress-energy criterion to predict the fracture behaviour of layered ceramics designed with internal compressive stresses[END_REF][START_REF] Leguillon | What is the tensile strength of a ceramic to be used in numerical models for predicting crack initiation?[END_REF], etc). An extensive review of the CC applications was recently proposed by Weißgraeber et al. [START_REF] Weißgraeber | A review of Finite Fracture Mechanics: crack initiation at singular and non-singular stress raisers[END_REF].

Another approach that is usually employed in order to describe material or structure failure is Cohesive Zone Model (CZM) [START_REF] Alfano | Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues[END_REF][START_REF] Barenblatt | The formation of equilibrium cracks during brittle fracture: general ideas and hypotheses[END_REF][START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF][START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF][START_REF] Tvergaard | The relation between crack growth resistance and fracture process parameters in elastic-plastic solids[END_REF][START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF] that consists in defining a traction-opening relationship that mimicks mechanically the failure process. CZM describes the relation between the cohesive stress and the displacement jump between both surfaces. When the traction reaches a critical value, progressive debonding takes place up to a critical opening corresponding to the local nucleation of a crack. A main difference between CZM and CC is the definition of a process zone between undamaged and completely failed zones. Nevertheless, both approaches basically require the same input data, namely the use of (at least) two failure parameters: a critical traction σ c and a critical energy release rate G c . CC and CZM can be employed in order to model the same fracture problem, which allows comparing both approaches.

Several authors compared crack initiation prediction using CC and CZM. For instance, Cornetti et al. [START_REF] Cornetti | Short cracks and V-notches: Finite Fracture Mechanics vs[END_REF] predicted similar initiation loading levels for both models in the case of short cracks and V-notch. Both methods were also employed for modeling fiber/matrix interface cracking in composite materials [START_REF] García | Fiber-size effects on the onset of fibermatrix debonding under transverse tension: A comparison between cohesive zone and finite fracture mechanics models[END_REF][START_REF] Poitou | Amorçage de fissure au voisinage des interfaces fibre/matrice : comparaison dun critère mixte et des modèles de zone cohésive. 18 eme Congrè Français de Mécanique[END_REF].

Poitou et al. [START_REF] Poitou | Amorçage de fissure au voisinage des interfaces fibre/matrice : comparaison dun critère mixte et des modèles de zone cohésive. 18 eme Congrè Français de Mécanique[END_REF] obtained similar initiation strains using either Dugdale model [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] or CC. García et al. [START_REF] García | Fiber-size effects on the onset of fibermatrix debonding under transverse tension: A comparison between cohesive zone and finite fracture mechanics models[END_REF] showed different degrees of agreement between both models, depending on some geometric features such as the fiber radius size. They also highlighted the fact that CC predicts an abrupt initiation on a finite length whereas CZM is based on a progressive debonding. However, Henninger et al. [START_REF] Henninger | Crack initiation at a V-notch-comparison between a brittle fracture criterion and the Dugdale cohesive model[END_REF] demonstrated that either CC or CZM result in crack initiation being forthwith followed by an unstable crack propagation that should lead the crack to directly jump to a finite length, but no details were given about this unstable crack propagation following nucleation. The spontaneity of crack nucleation was also discussed by Martin et al. [START_REF] Martin | Initiation of edge debonding: coupled criterion versus cohesive zone model[END_REF] in the case of edge debonding. Similar predictions of crack initiation level and load-displacement curve between CC and CZM were obtained provided that the characteristic fracture length of the material remained small enough compared to a characteristic length of the studied configuration. Otherwise, noticeable differences were observed since using CZM led to a large process zone and thus postponed crack initiation compared to CC. Cornetti et al. [START_REF] Cornetti | Finite Fracture Mechanics and cohesive crack model: weight functions vs. cohesive laws[END_REF] recently showed that the loading level predicted by FFM and CZM could even match in the case of large process zones, which was achieved by introducing weight functions in the stress condition of the CC. Most of the previously cited works focused on comparing the initiation loading level obtained using CC or CZM but include few analysis about the initiation crack length.

The objective of the present work is to provide additional insight about the comparison between CC and CZM approach for crack initiation prediction, especially concerning the spontaneous nature of crack nucleation, as well as the resulting crack length. The studied configuration is rhombus hole specimens under compression, which are chosen because they lead to a crack arrest after initiation, followed by a stable crack propagation. The CC is described in Section 2 and the CZM in Section 3. Numerical considerations for both methods are discussed in section 4 and the numerical prediction of crack initiation using both methods are compared in Section 5 .

Coupled criterion

Crack initiation determination

The initiation of a crack under quasi-static loading requires two separate conditions to be simultaneously fulfilled [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF]. The first requirement is a high enough stress state over the whole expected crack surface prior to crack nucleation. For an isotropic homogeneous material, this condition consists in comparing the stress normal to the crack plane to the material strength (Eq. 1) σ nn σ c over the whole expected crack surface.

(1)

The second requirement is derived from a balance, between an elastic state before and after crack initiation, of the potential energy (δW p ), the kinetic energy (δW k ) and the energy dissipated by crack nucleation (G c S, G c being the material toughness and S the newly created crack surface). Under quasi-static initial conditions, there is a production of kinetic energy (δW k > 0), which leads to the energy condition that compares the incremental energy release rate, i.e. the ratio between the energy dissipated by the crack opening (which corresponds to the change in potential energy ∆W (S) = W (0) -W (S) ) and the surface S of the crack, to the material critical energy release rate G c (Eq. 2),

G inc (S) = - ∆W (S) S = W (0) -W (S) S ≥ G c . (2)
Under the assumption of linear elasticity and small deformation, the stress is proportional to the applied load P whereas the potential energy of both the damaged and the undamaged material, and hence the incremental energy release rate, are proportional to the square of applied load (Eq. 3). In practice, only one calculation is required in order to compute k(S) and hence the stress criterion. Several calculations with varying crack surface allow A(S) and thus the energy criterion to be computed,

   σ nn (S) = k(S)P = σ c G inc = A(S)P 2 = G c (3)
The crack configuration at damage initiation (i.e. S * ) can be determined by solving Eq. 4, as the configuration for which both the energy and the stress conditions are satisfied for a minimum imposed loading P * [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum specimens under four point bending[END_REF][START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF],

as

A(S * ) k(S * ) 2 = G c σ 2 c . (4)

Stability of the crack at initiation

The coupled criterion allows the determination of finite crack nucleation configuration, i.e. a crack surface S * at a prescribed loading P * . The growth of the nucleated crack can then be studied considering the differential energy release rate G = -dW dS , which can be expressed as a function of the incremental energy release rate, its derivative and the crack surface (Eq. 5).

G(S) = - dW (S) dS = dG inc (S) dS S + G inc (S) (5)
It is worth noting that since G inc has already been computed to determine the nucleation configuration, no extra calculations are needed to compute G which is obtained using Eq. In the case of stable propagation, an increase in the loading is required so that the crack can grow, whereas an unstable crack propagation is achieved without any increase in the loading. At crack initiation, G inc = G c , therefore the stability of the propagation depends on the sign of dG inc (S)

dS

, i.e. on the monotony of the incremental energy release rate. A strictly increasing monotonic variation of G inc leads to an unstable propagation of the crack just after initiation which may even trigger the specimen failure [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum specimens under four point bending[END_REF][START_REF] Doitrand | Comparison between 2D and 3D applications of the coupled criterion to crack initiation prediction in scarf adhesive joints[END_REF]. Fig. 1 shows the case of a non monotonic variation of G inc , i.e. G inc attains a maximum for a crack surface S m [START_REF] Doitrand | Damage onset modeling in woven composites based on a coupled stress and energy criterion[END_REF][START_REF] Martin | Energetic conditions for interfacial failure in the vicinity of a matrix crack in brittle matrix composites[END_REF][START_REF] Weissgraeber | Crack nucleation in negative geometries[END_REF]. In most cases, crack initiation surface verifies S * < S m (Fig. 1) so that dG inc dS > 0 and

G > G inc = G c
, which leads to an unstable propagation of the crack. The crack will propagate at least up to a surface S min satisfying G(S min ) = G c . However, during crack propagation from S * to S min , the released energy is higher than the energy required for the crack growth (G > G c ). This energy excess may be consumed for further crack propagation up to a surface S max verifying Eq. 6. Besides, it can be shown that Eq. 6 is equivalent to

G inc (S max ) = G c . Smin S * (G(S) -G c )dS = Smax Smin (G c -G(S))dS (6)
In practice, only a part of the energy excess is consumed for further crack propagation while the other part may be dissipated [START_REF] Weissgraeber | Crack nucleation in negative geometries[END_REF]. Since it is really difficult to quantify the ratio between the energy excess that is dissipated or used for further crack propagation, the crack surface observed experimentally may be compared Crack initiation (S = S * ) is followed by unstable crack propagation up to an arrest surface comprised between S min and Smax.

the crack does not propagate without any increase in the loading, which means that the initiation crack surface is also an arrest surface (S * = S min = S max ) . In this case, the crack initiation surface predicted with the CC can directly be compared to that observed experimentally.

Cohesive zone model

A cohesive zone model defines a relation between the cohesive traction T and the opening displacement δ of a process zone. It is usually composed of a first part corresponding to an elastic undamaged behavior, for which the traction increases up to a critical value σ c corresponding to a characteristic opening displacement δ 0 . When the traction has reached σ c , debonding occurs with a decrease in the cohesive traction until a critical opening displacement δ c for which the traction becomes null, hence corresponding to the nucleation of a crack locally. The area under the traction-opening profile correspond to the work for separation and creation of free surfaces, equal to G c for an elastic material with no dissipation effects. Fig. 2a shows an example of traction-opening displacement curve in the case of a bilinear cohesive zone model. In addition to G c , σ c and δ c , the traction-opening profile can be different depending on the underlying failure mechanism [START_REF] Alfano | Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues[END_REF][START_REF] Alfano | On the influence of the shape of the interface law on the application of cohesive-zone models[END_REF][START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF][START_REF] Park | A unified potential-based cohesive model of mixed-mode fracture[END_REF][START_REF] Tvergaard | The relation between crack growth resistance and fracture process parameters in elastic-plastic solids[END_REF][START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF]. Alfano [START_REF] Alfano | On the influence of the shape of the interface law on the application of cohesive-zone models[END_REF] compared several cohesive models and concluded that the traction-opening profile can influence the response of the studied configuration, depending on the boundary value problem.

Moreover, Alfano showed that this influence is only noticeable in the vicinity of the maximum applied load (i.e., at crack initiation in the studied configuration). This conclusion is in agreement with that of Acary and Monerie [START_REF] Acary | Nonsmooth fracture dynamics using a cohesive zone approach[END_REF] who showed that all the cohesive zone models are equivalent for studying stable rectilinear crack propagation but not for crack initiation and branching. The influence of the traction-opening displacement profile on crack initiation is studied in Section 4.3 by employing the model of Park et al. [START_REF] Park | A unified potential-based cohesive model of mixed-mode fracture[END_REF][START_REF] Park | Computational implementation of the PPR potential-based cohesive model in ABAQUS: Educational perspective[END_REF] which enables different traction-opening displacement profiles to be obtained by varying a single shape parameter (for fixed fracture parameters G c and σ c ) . For the configuration under investigation in this study, cracking occurs in mode I. Therefore, mode mixity, which can also be taken into account using CC [START_REF] Carrère | Comparison between models based on a coupled criterion for the prediction of the failure of adhesively bonded joint[END_REF][START_REF] García | A model for the prediction of debond onset in sphericalparticle-reinforced composites under tension. Application of a coupled stress and energy criterion[END_REF][START_REF] Martin | Initiation of edge debonding: coupled criterion versus cohesive zone model[END_REF] or CZM [START_REF] Alfano | Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues[END_REF][START_REF] Alfano | On the influence of the shape of the interface law on the application of cohesive-zone models[END_REF][START_REF] Park | A unified potential-based cohesive model of mixed-mode fracture[END_REF], will not be considered in this work. Contrary to the CC presented in Section 2 for which crack process zone undamaged zone the definition of a crack is straightforward, the definition of a cohesive zone model introduces the concept of process zone (cf. Fig. 2b). In the following, the crack length for CZM is defined as the traction-free zone which does not include the process zone (cf. Fig. 2b). The initiation loading is therefore defined as the minimum loading for which a crack nucleates (i.e. δ = δ c and T = 0 in the CZM).

T=0 δ<δ 0 δ 0 <δ<δ c T>0 (a) (b)

Numerical implementation of CC and CZM

Finite Element models

The configuration chosen for the comparison of CC and CZM is quasi-static compression of rhombus hole specimen, whose geometry is depicted in Fig. 3a. A uniaxial compression in the y-direction is generated by imposing a vertical displacement on the specimen top edge. Under such a loading, a mode I vertical crack initiates at the rhombus hole upper and lower V-notch due to the lateral expansion. The symmetric configuration allows to model only one specimen quarter, for which a mesh example is presented in Fig. 3b.

The specimen length and width respecively are l = 60mm and w = 40mm. The mechanical properties chosen for both CC and CZM computation are E = 1600GPa, ν = 0.37, σ c = 80MPa and G c = 250J/m 2 , which are representative of PMMA. All the calculations presented in this work have been performed in 2D using Abaqus standard v6.13. Since only 2D models are considered in this work, in the following, the crack is described by its length L rather than its surface S. The initiation crack length (solution of Eq. 4) is noted L * , and the crack arrest length L arrest .

CC implementation

The implementation of the CC for crack initiation prediction in rhombus hole specimens requires both the stress and the energy conditions described in Section 2 to be computed. This is done by FE calculations on undamaged and damaged configurations of rhombus hole specimens, whose dimensions are depicted in Fig. 3a. The stress criterion is computed from a calculation on undamaged specimen meshed with 4 node elements (Fig. 3b). The energy criterion is obtained by computing the potential energy of both the undamaged and damaged configuration for several crack lengths, which requires the generation of the corresponding crack meshes. In 2D, the crack is only defined by one parameter, namely its length L along (Oy)-axis so that 2D cracked mesh generation is straightforward by doubling the nodes on the presupposed crack path. Once the stress and the energy conditions have been computed, the method presented in Section 2 is applied in order to determine the crack initiation and arrest lengths, as well as the force at initiation. Fig. 4 and Table 1 show the incremental energy release rate and the initiation force and crack arrest lengths predicted by the CC using different mesh refinement levels along the crack path in the case of a hole with β=90 deg. angle and s × cos( β 2 )=4.95mm half-diagonal. The number of calculations corresponds to the number of nodes that are successively doubled in order to compute the energy release rate as a function of the crack length. maximum mesh size (mm) 2. 1. 0.5 0.5 0.5 0.5 minimum mesh size (mm) 0.5 ×10 -3 0.5 ×10 -3 0. It is clear that a too coarse mesh does not allow correctly capturing the incremental energy release rate maximum (e.g. for a 0.5mm minimum mesh size). It can however be noted that the part of the curve corresponding to crack lengths larger than 1mm is well captured whatever the mesh size. Another paramount feature concerns the prediction of the initiation crack length L * (solution of Eq. 4), which adds a second level of mesh refinement. Indeed, it can be shown either analytically [START_REF] Martin | A coupled strength and toughness criterion for the prediction of the open hole tensile strength of a composite plate[END_REF] or numerically [START_REF] Martin | Understanding the tensile strength of ceramics in the presence of small critical flaws[END_REF] that the initiation length is a fraction of the material characteristic length

L mat = EG c σ 2 c
. Therefore, fine enough mesh size is required to catch this length when solving Eq. 4. The material under investigation exhibits a characteristic length around 60 µm, hence leading to initiation lengths as small as 20 µm. A too large minimum mesh size leads to an erroneous estimate of the crack initiation force and arrest lengths (Table 1). Since the mesh size must not be as small to capture the post-maximum part of the curve as it must be to catch the initiation length, it is possible to unrefine progressively the mesh along the crack path, which makes the computation of the CC faster. Indeed, it allows decreasing both the number of degrees of freedom and the number of calculations to solve the CC. For instance, using a mesh including a 0.5 × 10 -3 mm minimum size, similar predicted initiation forces and crack arrest lengths are obtained imposing maximum mesh sizes varying from 0.5mm to 2mm. However, increasing the maximum mesh size from 0.5mm to 2mm allows reducing both the number of calculations needed to solve the CC and the number of degrees of freedom in each calculations (cf. Table 1). In the following, all the calculations for the CC solution have been performed using a 0.5µm minimum mesh size as well as a progressive mesh unrefinement along the crack path to a 2mm maximum mesh size, so that around 100 linear elastic calculations are needed to solve the CC.

Fig. 5 shows the crack length as a function of the force predicted using the CC for a hole with β = 90deg.

angle and s × cos β 2 = 4.95mm half-diagonal. The crack arrest lengths obtained taking into account or not the energy excess available for further propagation of the nucleated crack respectively correspond to the crack arrest upper and lower bounds. Crack propagation (i.e. the variation of L as a function of F after initiation) has been computed using the same incremental algorithm as proposed by Doitrand et al. [START_REF] Doitrand | Mesoscale analysis of damage growth in woven composites Compos Part A[END_REF] to study crack propagation in woven composites. At crack initiation, if the energy surplus available for further crack propagation is not used (cf. Fig. 1), the crack reaches its arrest length (i.e the lower bound predicted by the CC) so that G = G c and dG dL ≤ 0, corresponding to a stable crack propagation. However, if further crack propagation is allowed by taking into account the energy surplus, the crack arrest length (i.e. the upper bound predicted by the CC) corresponds to a configuration for which G < G c . Therefore, a sufficient increase in the loading is necessary for the crack to propagate again (i.e. so that G = G c ). Since except just after initiation (in the case the energy excess is taken into account), no unstable crack propagation phase is observed, similar variations of the crack length as a function of force is obtained whether the external dissipation is taken into account or not [START_REF] Doitrand | Mesoscale analysis of damage growth in woven composites Compos Part A[END_REF].

CZM implementation

Contrary to the CC approach, which requires several linear elastic calculations to be computed, only one non-linear calculation is needed in the case of CZM. The use of CZM requires sufficiently fine mesh discretization so as to accurately represent the process zone ahead of the crack. The order of magnitude of the process zone length is generally a fraction of the material characteristic length L mat [START_REF] Acary | Nonsmooth fracture dynamics using a cohesive zone approach[END_REF][START_REF] Barenblatt | The formation of equilibrium cracks during brittle fracture: general ideas and hypotheses[END_REF][START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF][START_REF] Martin | Initiation of edge debonding: coupled criterion versus cohesive zone model[END_REF][START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF]. In the present case, the process zone length is 50 µm, therefore a mesh size of 5µm has been imposed along the crack path so that the process zone contains around 10 elements. The other numerical details of the computations are given in the rest of this section. Fig. 6 shows the crack length as a function of force obtained in the case of a bilinear CZM. The elastic slope of the CZM, whose influence is studied in the following, is set to K CZM = 10 8 MPa/mm. Similarly to the CC approach, a crack jump at initiation is obtained. Indeed, just after nucleation, the crack propagates in an unstable manner and directly reaches a length of 0.6mm in a really short period (corresponding to an imposed displacement increment lower than 10 -5 mm). This result supports one of the main assumptions made when using the CC, namely the spontaneous crack initiation over a given length. The initiation configuration is described by the force level together with the crack arrest length after the jump. These quantities will be compared to the CC results in next section.

The softening behavior of the CZM may lead to convergence issue in quasi-static numerical simulations [START_REF] Acary | Nonsmooth fracture dynamics using a cohesive zone approach[END_REF][START_REF] Chaboche | Interface debonding model: a viscous regularization with a limited rate dependency[END_REF][START_REF] Gao | A simple technique for avoiding convergence problems in finite element simulations of crack nucleation and growth on cohesive interfaces[END_REF]. In particular, this is observed in the case of unstable crack propagation, for which solution jumps occur. Several methods exist in order to overcome this problem such as introducing a fictitious viscous regularization [START_REF] Chaboche | Interface debonding model: a viscous regularization with a limited rate dependency[END_REF][START_REF] Gao | A simple technique for avoiding convergence problems in finite element simulations of crack nucleation and growth on cohesive interfaces[END_REF]. An other way consists in performing dynamic rather than quasi-static numerical simulations [START_REF] Acary | Nonsmooth fracture dynamics using a cohesive zone approach[END_REF], which naturally introduces a relaxation that allows avoiding convergence issues. Such a choice can be justified by the fact that an unstable crack propagation is itself a dynamic process, even under quasi-static loading. All the computations involving CZM have been performed via dynamic numerical simulations, using the PMMA density ρ = 1.18 × 10 -3 g/mm 3 .

Since we assume a quasi-static loading state in the CC solution (even if crack initiation is considered as a spontaneous and dynamic phenomenon), it is paramount to study the influence of the imposed displacement rate on crack initiation in CZM simulation in order to provide a consistent comparison between both methods.

In order to check the overall quasi-static state, several calculations with varying imposed displacement rate have been performed. The initiation force and crack arrest length variations as a function of the imposed displacement rate are shown in Fig. 7. The displacement increment during which the crack jump occurs are summarized in Table 2. The initiation force and the crack arrest length decrease with decreasing displacement rate and reaches a plateau for displacement rates lower than 1mm/s, corresponding to an overall quasistatic state. The displacement increment corresponding to crack initiation also decreases with decreasing displacement rate. Therefore, whatever the applied displacement rate, the crack jump occurs within a time increment of 5+1×10 -5 s. In the following, a 1×10 -1 mm/s displacement rate has been imposed in simulations involving CZM so that the predicted initiation force and crack arrest length do not depend on the displacement rate.

In CZM models, the elastic part of the traction-opening displacement curve (i.e for T < σ c and δ < δ 0 ) introduces an artificial cohesive element stiffness K CZM = σc δ0 . However, the purpose of the cohesive zone elements is to simulate fracture which means that the stiffness of the cohesive element must be high enough so that it does not affect the overall rigidity of the model [START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF]. Fig. 8 shows the initiation force and crack arrest length obtained varying the CZM stiffness from 10 5 MPa/mm to 10 10 MPa/mm. Too small values of the CZM stiffness delay crack initiation and lead to a slight underestimate of the crack arrest length. However, Displacement rate (mm/s) 1×10 -2 1×10 -1 1 10 100 1000

Crack jump displacement increment (mm) 4×10 -7 5×10 -6 6×10 -5 4×10 -4 5×10 -3 5×10 -2 similar initiation forces and crack arrest lengths are obtained using a stiffness larger than 10 7 MPa/mm. In the following, all the calculations have been performed using a CZM stiffness K CZM = 1 × 10 8 MPa/mm.

It can be noted that increasing the CZM stiffness also increases the number of iterations and decreases the minimum required time increment (cf. Table 3), as already observed by Turon et al. [START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF] Previous analyses concern a bilinear CZM, corresponding to a linear softening part of the traction-opening displacement curve. However, for given fracture parameters G c and σ c , it is possible to vary the CZM profile, which may influence crack initiation [START_REF] Acary | Nonsmooth fracture dynamics using a cohesive zone approach[END_REF][START_REF] Alfano | On the influence of the shape of the interface law on the application of cohesive-zone models[END_REF]. Park et al. [START_REF] Park | A unified potential-based cohesive model of mixed-mode fracture[END_REF][START_REF] Park | Computational implementation of the PPR potential-based cohesive model in ABAQUS: Educational perspective[END_REF] proposed a model that is able to represent numerous CZM profiles by varying a single parameter α. A Dugdale-like model is obtained when α tends towards 1 (α must be strictly larger than 1 for numerical reasons). If α = 2, a bilinear model is obtained and α > 2 leads to trilinear-like models. Graphical representations of the traction-opening displacement curve obtained for α=1.01, 2 and 5 are shown in Fig. 9a. Whatever the CZM profile, the elastic slope for opening displacement lower than the characteristic opening displacement δ 0 is also set to K CZM = 1 × 10 8 MPa/mm so that the overall rigidity of the model is not altered. Keeping the values of K CZM , G c and σ c constant, the only parameter varying when changing the CZM profile is the critical opening displacement δ c , whose value is given as a function of the shape parameter α in Table 4.

Fig. 9 shows the force and the initiation crack arrest length obtained for different CZM profiles as a function of the shape parameter. Both crack initiation force and arrest length increases as the CZM profile tends to a Dugdale-like profile. This might be explained by the fact that when the critical opening displacement decreases (i.e. when α decreases), the length on which the material strength σ c has to be attained before failure increases (in order to keep a constant fracture energy), which required a larger imposed loading and thus leads to a larger crack arrest length.

K CZM 1×10 5 1×10 6 1×10 7 1×10 8 1×10 9 1×10

Comparison between CC and CZM

Crack initiation prediction

The initiation forces and crack arrest lengths predicted using CC and CZM are compared in this section for several rhombus hole specimen configuration. The hole angle β varies from 70 deg . to 110 deg . while its half diagonal s × cos β 2 is kept constant. Fig. 10 depicts the initiation force (Fig. 10a) and crack arrest lengths (Fig. 10b) obtained using CC or CZM with various profiles as a function of the hole angle β. The initiation force decreases with increasing angles. This trends is predicted using both CC or CZM with different profiles. The initiation force predicted using the CC is comprised between the initiation force obtained with a Dugdale-like (α = 1.01) or a trilinear-like CZM (α = 5). The differences between the bilinear model and CC initiation force predictions are around 3%. The crack arrest lengths increase with increasing hole angle. This trend is, once again, predicted using both CZM or CC. The difference between the upper and lower

Discussion

Previous works about the comparison between CC and CZM for crack initiation prediction were mainly focused about confronting the initiation loading predicted by both methods [START_REF] Cornetti | Short cracks and V-notches: Finite Fracture Mechanics vs[END_REF][START_REF] García | Fiber-size effects on the onset of fibermatrix debonding under transverse tension: A comparison between cohesive zone and finite fracture mechanics models[END_REF][START_REF] Henninger | Crack initiation at a V-notch-comparison between a brittle fracture criterion and the Dugdale cohesive model[END_REF][START_REF] Martin | Initiation of edge debonding: coupled criterion versus cohesive zone model[END_REF][START_REF] Poitou | Amorçage de fissure au voisinage des interfaces fibre/matrice : comparaison dun critère mixte et des modèles de zone cohésive. 18 eme Congrè Français de Mécanique[END_REF]. In all of these studies, only one CZM was employed so that varying the CZM profile was not accounted for. In the present study, we highlight that the choice of the CZM has a non-negligible influence on the predicted initiation force.

For instance, differences up to 15% (cf. Fig. 10a) are obtained between the initiation force predicted using a Dudgale-like or trilinear-like CZM.

An important aspect that can also be compared is the computational cost of both approaches. As regards CC, following the mesh analysis presented in section 4.2, the whole CC analysis (generation of the model, calculations and post-treatment) requires around 15 min to be completed. Since a sufficiently refined mesh is required over the whole crack path when using CZM (cf. section 4.3), meshes with around 500000 degrees of freedom are obtained, corresponding to calculations duration longer than 10H. Such a refinement along the crack path is required if the crack arrest length is not known, since large crack arrest lengths (e.g., for β = 110deg.) can be obtained. The calculation time can be reduced to around 1H in the case of smaller crack arrest lengths (e.g., for β = 70deg.) by unrefining the mesh along the crack path at a certain distance of the crack arrest length and only considering the initiation phase, which is only possible if the crack arrest length is known a priori.

One of the main advantages of the CC compared to CZM is that the toughness and strength values are only implemented during the post-processing step. Therefore, provided that the mesh size is small enough to correctly represent G inc and capture the initiation length L * (cf. section 4.2), the FE calculations with different crack lengths can be computed once to obtain initiation force and crack arrest length for different values of G c and σ c . This is particularly convenient in the case where these fracture parameters are unknown and have to be determined by inverse confrontation to experimental results [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum specimens under four point bending[END_REF][START_REF] Doitrand | Comparison between 2D and 3D applications of the coupled criterion to crack initiation prediction in scarf adhesive joints[END_REF]. This is not the case employing CZM, which requires several calculations with varying G c and σ c to be computed. Conversely, in the case where several crack configurations have to be studied for given fracture parameters, the number of calculations to perform in order to solve the CC exponentially increases with the number of parameters describing the crack [START_REF] Doitrand | Damage onset modeling in woven composites based on a coupled stress and energy criterion[END_REF]. In this case, the CZM can be a powerfull approach since it allows taking into account, for instance, different crack paths in a single computation, obviously leading to more costly computations due to an increase in the number of dof and more restrictive convergence conditions.

Conclusion

Crack nucleation in rhombus hole specimens under compression modeled using either CZM or CC leads to the conclusion that at initiation, there is a crack jump to a given length at a certain loading. Even under overall quasi-static loading conditions, the crack jump is a dynamic and instantaneous phenomenon. The CC allows determining the initiation loading as well as a lower and upper bounds for crack arrest by considering the unstable crack propagation just after initiation, taking into account or not the energy surplus for further crack propagation. The initiation loading and the crack arrest length can also be determined using CZM.

For fixed fracture parameters G c and σ c , both quantities depends on the CZM profile.

CC and CZM predict similar trends of the initiation force and crack arrest length as a function of the rhombus hole angle. With increasing hole angle, the initiation force decreases whereas the crack arrest length increases. For various geometrical configurations, the initiation force predicted using the CC lies within the range of those obtained with different CZM profiles. Similarly, the range of crack arrest lengths obtained with different CZM profiles comprises the crack arrest length lower bound obtained with the CC. For the configurations under investigation, similar results are obtained with the CC and bilinear CZM. Confronting CZM and CC appears to be a way to identify the CZM profile for given material parameters.

CZM requires a sufficiently refined mesh over the whole crack path in order to correctly account for the process zone ahead of the crack tip. For CC computation, the mesh must be sufficiently fine to catch the initiation length and correctly describe the incremental energy release rate but it is possible to unrefine it along the crack path to decrease both the number of calculations and the number of degrees of freedom in each calculations (and thus to save computational time). In both cases, the minimum mesh size can be related to the material characteristic length EGc . In the present study, the CC appears to be more computationally efficient that CZM. A main advantage of CC approach compared to CZM is that all the calculations can be done without assuming any values of G c and σ c (except in order to ensure a sufficiently fine minimum mesh size), that are only implemented while in the calculation post-treatment. Therefore, this approach can be really convenient in order to determine unknown fracture parameters by inverse confrontation to experimental results [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum specimens under four point bending[END_REF][START_REF] Doitrand | Comparison between 2D and 3D applications of the coupled criterion to crack initiation prediction in scarf adhesive joints[END_REF][START_REF] Leguillon | An attempt to extend the 2d coupled criterion for crack nucleation in brittle materials to the 3d case[END_REF]. On the contrary, if many possible crack configurations have to be studied, the number of calculations rapidly increases using CC, therefore CZM might be preferred in this case.

The present analysis exclusively concerns the comparison of two numerical methods. Future work will cover an experimental analysis of crack initiation in rhombus hole specimens so as to evaluate the validity of the models, especially regarding the predicted spontaneous crack initiation.

5 .

 5 After crack nucleation, several propagation scenarii may happen depending on the differential energy release rate: o No propagation : G < G c , o Stable propagation: G = G c and dG dS 0, o Unstable propagation: G > G c or (G = G c and dG dS 0).

  to the lower (S min verifying G(S min ) = G c ) and upper (S max verifying G inc (S max ) = G c ) bound estimates of the crack arrest surface. A particular case concerns crack initiation surface corresponding to the maximum of the incremental energy release rate (S * = S m ) (Fig. 1), leading to dG inc dS = 0 and G = G inc = G c . Therefore,

Figure 1 :

 1 Figure 1: Stress to strength, incremental and differential energy release rate to toughness ratios as a function of the crack surface.

Figure 2 :

 2 Figure 2: (a) Traction-opening displacement curve (example of bilinear model) and (b) schematic representation of the process zone.

Figure 3 :

 3 Figure 3: (a) Geometry of the rhombus hole specimen and (b) mesh of a quarter of specimen.

Table 1 :

 1 Crack initiation force and arrest length lower and upper bounds predicted using CC as a function of the minimum mesh size along the crack path in the case of a rhombus hole with β=90 deg. angle and s × cos( β 2 )=4.95mm half-diagonal.

Figure 4 :

 4 Figure 4: Incremental energy release rate as a function of the crack length obtained for different mesh refinement levels (The size indicated in the legend is the minimum mesh size).

Figure 5 :

 5 Figure 5: Crack length as a function of force obtained employing the CC in the case of a hole with β=90 deg. angle and s × cos( β 2 )=4.95mm half-diagonal. The initiation crack arrest upper and lower bounds are obtained if the energy surplus at initiation is used (no external dissipation) or not (external dissipation) for further crack propagation.

Figure 6 :

 6 Figure 6: Crack length as a function of force obtained with bilinear CZM in the case of a hole with β=90 deg. angle and s × cos( β 2 )=4.95mm half-diagonal.

Figure 7 :

 7 Figure 7: (a) Initiation force and (b) crack arrest length as a function of the imposed displacement rate in the case of a hole with β=90 deg. angle and s × cos( β 2 )=4.95mm half-diagonal using a bilinear CZM model.

Figure 8 :

 8 Figure 8: (a) Initiation force and (b) crack arrest length as a function of the bilinear CZM stiffness in the case of a hole with β=90 deg. angle and s × cos( β 2 )=4.95mm half-diagonal.

Figure 9 :

 9 Figure 9: Initiation (a) force and (b) crack length as a function of the shape parameter α of the CZM.

Figure 10 :

 10 Figure 10: Initiation (a) force and (b) crack arrest length obtained using CC or CZM with different profiles as a function of rhombus hole angle β for a rhombus half diagonal s × cos( β 2 ) = 4.95mm.

σ 2 c

 2

Table 2 :

 2 Displacement and time increment during which occurs the initiation crack jump for several imposed displacement rates.

  10

	Number of iterations	155	180	187	206	219	328
	Minimum time increment (s) 1.2×10 -7 1.0×10 -7 1.4×10 -8 1.5×10 -9 4.1×10 -10 2.4×10 -10

Table 3 :

 3 Critical opening displacement of the cohesive zone as a function of the shape parameter α.

	α	1.01	1.2	2	3.5	5
	δ c (mm) 0.003 0.004 0.006 0.01 0.015

Table 4 :

 4 Critical opening displacement of the cohesive zone as a function of the shape parameter α.


R Niranjana

David Stephen
email: davidstephen_dav@yahoo.co.in

P Justin

K Saravanan

Piero Macchi

R Niranjana Devi

A David Stephen

Christian Jelsch

Topological and electrostatic properties of diclofenac molecule as a non-steroidal anti-inflammatory drug: An experimental and theoretical study

Keywords: Experimental Charge density, Hirshfeld surface analysis, topological properties, electrostatic interaction, Molecular Dynamics, Quantum Crystallography

Introduction

Diclofenac is the most prominent member of the Non-Steroidal Anti-Inflammatory Drugs (NSAIDs) [START_REF] Yamazaki | Diclofenac, a non-steroidal anti-inflammatory drug, suppresses apoptosis induced by endoplasmic reticulum stresses by inhibiting caspase signaling[END_REF], as it inhibits the enzyme cyclo-oxygenase COX-2 [START_REF] Hinz | Aceclofenac spares cyclooxygenase 1 as a result of limited but sustained biotransformation to diclofenac[END_REF]. Diclofenac does also inhibit the Transthyretin (TTR) [START_REF] Connelly | Structure-based design of kinetic stabilizers that ameliorate the transthyretin amyloidoses[END_REF] protein that carries the thyroid hormone thyroxine (T4) and retinol-binding protein bound to retinol. This is how transthyretin gained its name: transports thyroxine and retinol. The TTR protein exhibits a β-sheet homo-tetramericthree-dimensional structure and has two channels at the dimer-dimer interface. This protein may be the main reason for familial amyloid polyneuropathy (FAP), an amyloid diseases caused by the TTR fibril formation over the heart and peripheral nerves. Diclofenac have been evaluated as a strong in vitro TTR amyloid inhibitor in the reverse mode which makes the rotation of dihalogen-substituted ring by 90° with respect to phenylacetic moiety [START_REF] Lima | Identification of a novel ligand binding motif in the transthyretin channel[END_REF]. Diclofenac is the best alternative for NSAID drug flufenamic (FLU) [START_REF] López-Mejías | Nonamorphism in flufenamic acid and a new record for a polymorphic compound with solved structures[END_REF], which is also a biphenylamine derivative with carboxylate substituent, however FLU binds with TTR in "forward mode". The in vitro analysis of this inhibition results with 17-56% reduction of TTR fibril formation [6].

To the best of our knowledge, the experimental electron density study for this compound has not been carried out till date. In order to shed light on the electronic structure of molecules and on the molecular recognition processes, high resolution X-ray diffraction measurements on crystals and charge density analysis were carried out. These kinds of studies enable the determination of the electron density distribution in molecular solids and, from this, a careful analysis of the electrostatic interactions. The charge density derived from experimental models can be compared with the theoretically derived, calculated in vacuo or under periodic boundary conditions, typically using density functional theory (DFT). The analysis of the static density, however, may not be enough for a full comprehension of the molecular aggregation and the intermolecular interactions. For this reason, molecular dynamics (MD) simulations may be useful to refine the binding conformation of the complex structure and to investigate the structural ability of a protein in the presence of ligands. Having this strategy in mind, we have undertaken a charge density analysis of diclofenac in order to obtain, from experiment and theory, an overview of its structural, topological, reactivity and electrostatic properties. In addition, molecular dynamic simulations were carried out to depict a sensible model of the actual binding mechanism. This study will provide pathway for exploration of new compounds with fewer side effects and a more sophisticated profile than already existing drug molecules, which is the focus of chemists and pharmacologists.

Experimental details

Crystallization, X-ray data collection and structure solution

The crystals of diclofenac were grown by slow evaporation method using ethanol solution.

Among the harvested crystals, a good quality single crystal (ca. 0.12 x 0.14 x 0.19 mm) was chosen for high resolution X-ray diffraction experiment. All measurements were made on an Oxford Diffraction Super Nova area-detector diffractometer using mirror optics monochromated Mo Kα radiation (λ= 0.71073 Å) and Al filtered [START_REF] Macchi | Low-energy contamination of Mo microsource X-ray radiation: analysis and solution of the problem[END_REF]. A total of 3600 frames were collected in 31 runs with data exposure of 10+10 and 45+45 seconds up to θ<48°, with a crystal-detector distance of 65.0 mm, at T = 110(2) K. Data reduction was performed using the CrysAlisPro program [START_REF] Macchi | Low-energy contamination of Mo microsource X-ray radiation: analysis and solution of the problem[END_REF]. The intensities were corrected for Lorentz and polarization effects, and an analytical absorption correction was applied, using SCALE3 ABSPACK in CrysAlisPro [8]. The smart-background data correction was adopted.

Experimental multipole refinement

The model was refined with MoPro software [START_REF] Jelsch | Advances in protein and smallmolecule charge-density refinement methods usingMoPro[END_REF] using Hansen and Coppens [START_REF] Hansen | Testing aspherical atom refinements on small-molecule data sets[END_REF] multipole model for pseudo-atom electron density. The C, N and O atoms were refined up to octapole level, hydrogen atoms up to dipole level and chlorine atom was refined up to hexadecapole.

The core and valence spherical scattering factors were calculated from Su & Coppens [START_REF] Su | Relativistic X-ray Elastic Scattering Factors for Neutral AtomsZ= 1-54 from Multiconfiguration Dirac-Fock Wavefunctions in the 0-12Å-1sinθ /λ Range, and Six-Gaussian Analytical Expressions in the 0-6Å-1Range[END_REF] and the anomalous dispersion coefficients were taken from the International Tables for Crystallography [START_REF] Koch | Twinning[END_REF]. The parameters of the Slater functions used to describe the multipoles are available in the deposited CIF file. The weighting scheme applied was Whkl = 1/ 2 (I) and all reflections were used in all refinements. At first, after structural refinement, the Plm multipole populations were refined. In the subsequent steps, the structural and charge density parameters were refined successively, in several cycles. In the final stages, all parameters were refined simultaneously till convergence.

According to the values found by the SHADE3 server [START_REF] Madsen | SHADE3server: a streamlined approach to estimate H-atom anisotropic displacement parameters using periodicab initiocalculations or experimental information[END_REF], the Uij thermal displacements parameters were fixed for H atoms. The H-X distances of H atoms were restrained to the standard values obtained from neutron diffraction studies [START_REF] Allen | Tables of bond lengths determined by x-ray and neutron diffraction. Part 1. Bond lengths in organic compounds[END_REF] with a restraint sigma of 0.01Å.

The crystallographic details of the experimental refinements are given in Table 1. The chemically equivalent atoms in the structure were restrained to have similar charge density parameters (,',Pval, Plm) with a restraint sigma r=0.01 [START_REF] Zarychta | Charge-density analysis of 1-nitroindoline: refinement quality using freeRfactors and restraints[END_REF]. The electron density of atoms was restrained to follow the local symmetry (r=0.01). For hydrogen, the  and ' parameters of an atom were restrained to have similar values (k=0.01) and the  values were restrained to be equal to 1.16 (k=0.03).

The most important properties such as the total electron density ρbcp(r), the Laplacian , the ellipticity ε at the bond critical points and the eigenvalues (λ1,λ2,λ3) were calculated and the experimental deformation maps were estimated and plotted using VMoPro module of the MoPro software [START_REF] Jelsch | Advances in protein and smallmolecule charge-density refinement methods usingMoPro[END_REF]. The residual indexes and the other parameters of multipole refinement are presented in Table 1. The normal probability plot of Iobs-Icalc and the evolution of Iobs/Icalc as a function of reciprocal resolution sin/ is shown in Fig. Sup1 and Sup2 [START_REF] Zhurov | Optimization and evaluation of data quality for charge density studies[END_REF].

The fractal dimension plot of the residual electron density is shown in Fig. Sup3 [START_REF] Meindl | Foundations of residual-density analysis[END_REF].

Theoretical calculations

A density functional calculation (DFT) was performed in order to understand the electronic and reactivity properties for the drug diclofenac using GAUSSIAN09 software [START_REF] Frisch | Gaussian09 Revision D.01[END_REF]. Singlepoint calculations were performed once the convergence was reached. The absence of imaginary frequencies showed that the minimum energy structure had been achieved.

Similarly, after the MD simulation, the diclofenac molecule with active site amino acids was taken for the single point energy minimization to generate the wave function. The obtained wave function for the gas phase of diclofenac molecule and active site region were given as an input to the AIM2000 package [START_REF] König | AIM2000-a program to analyze and visualize atoms in molecules[END_REF] for calculating the topological properties.

Molecular Dynamics Simulation

To perform MD simulation, three-dimensional coordinates of the human transthyretin complexed with diclofenac (1DVX) complex were retrieved from the Protein Databank. The input parameters for the MD simulations were prepared using the default protein force fields and the ligand parameters were generated using GAFF force field. The complex was placed at the centre of the box and TIP3P water molecules [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF] from the AMBER solvent library were added into it. Necessary sodium ions were added to neutralize the system. Periodic boundary conditions were used to treat the boundary of the system to minimize the edge effects. Particle Mesh Ewald (PME) was used to treat electrostatic interactions, the PME grid size was 73×69×82Ǻ. The resulting system was minimized for 1000 steps and the temperature was maintained at 300K throughout the simulation and pressure (1 bar) using the same Langevin thermostat and Berendsen barostat used in the heating process [START_REF] Jones | Adaptive stochastic methods for sampling driven molecular systems[END_REF][START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF]. 20 ns MD simulation was carried out to find out the stability of diclofenac at the active site of human transthyretin using the molecular dynamics software NAMD [START_REF] Phillips | Scalable molecular dynamics with NAMD[END_REF] with the time step of one femto-second.

VMD [START_REF] Humphrey | VMD: Visual molecular dynamics[END_REF] and CPPTRAJ [START_REF] Roe | PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data[END_REF] software were used to generate the MD trajectory.

Results and discussion

Crystal structure description

The crystal volume of diclofenac at 110 K is 3.9% smaller than at ambient conditions, as reported by Muangsinet al. [START_REF] Muangsin | Structure determination of diclofenac in a diclofenaccontaining chitosan matrix using conventional X-ray powder diffraction data[END_REF]. The molecular structure is shown in Figure 1 as an ORTEP plot. The 2,6-dichlorophenyl ring and the acid chain lie almost perpendicular to the central ring. Despite forming a strong hydrogen bonded dimer, the single and double bond characters of the C-O bonds are quite evident: C1-O2 is shorter [1.2314(6)Å] than C1-O1 [1.3080(8)Å] (see also Table Sup 1). Compared with the standard C-N bond lengths [START_REF] Fábián | Cocrystals of fenamic acids with nicotinamide[END_REF], it is evident that some delocalization occurs in the C-N bonds [C8-N1=1.414(5)Å, C9-N1=1.3939 [START_REF] López-Mejías | Nonamorphism in flufenamic acid and a new record for a polymorphic compound with solved structures[END_REF]Å]. This is further confirmed by C8-N1-C9 angle [122.73(3) o ] being larger than that of a tetrahedral angle which is expected for N1 atom. The electron lone pair of N1 is weakly visible in the electron density map and is involved in the resonance. Another reason for the large angle might be the C9 atom is attached with the 2,6-dichlorophenyl ring. In order to get details regarding the distortion in bond angles, the deformation of C-C-C valence angles has been calculated for the 2,6-dichlorophenyl ring and the other phenyl ring separately. The root means square deviations from the ideal 120° value of the hexagonal angles is computed as: It is found that the distortion in angle for the 2,6-dichlorophenyl ring is higher (σhexagon=2.2°) than that of the other phenyl ring of the molecule (σhexagon=0.65°) and the reason might be presence of the two Cl atoms and their electron withdrawing effect on the C6 ring.

Crystal packing and intermolecular interactions

The main feature of the supramolecular packing is the coupling of two carboxylic groups forming a strongly hydrogen bonded dimer (Fig. 2a surface analysis [START_REF] Mckinnon | Novel tools for visualizing and exploring intermolecular interactions in molecular crystals[END_REF][START_REF] Spackman | A novel definition of a molecule in a crystal[END_REF] summarises all the contributions and it certainly highlights major role played by O1-H1B•••O2 interaction (Fig. 3). The hydrogen bonding details are given in Table 2 There is limited non-parallel cycle stacking present between 2,6-dichlorophenyl rings. A dimer formed via an inversion centre results in a Cl•••Cl distance of 3.576 Å which is shorter than the sum of Van der Waals radii (Fig. 3b). The experimental Cl•••Cl interaction which is traced in our molecule is comparable with the type I trans Cl•••Cl geometry as reported by Hathwar and Guru Row [START_REF] Hathwar | Nature of ClꞏꞏꞏCl Intermolecular Interactions via Experimental and Theoretical Charge Density Analysis: Correlation of Polar Flattening Effects with Geometry[END_REF].

The fingerprint analysis (Fig. 4) [START_REF] Spackman | From charge densities to crystal engineering[END_REF] shows that H•••H contacts (Fig. 4b) give the largest contribution, which is about 32.5%. The Cl 4d) also render a significant contribution of 18.7% and 13.6% respectively to the packing.

•••H/H•••Cl (Fig.4c) and O•••H/H•••O hydrogen bonds (Fig.

Contact Enrichment Ratio Analysis

The contacts enrichments [START_REF] Jelsch | The enrichment ratio of atomic contacts in crystals, an indicator derived from the Hirshfeld surface analysis[END_REF]have been analysed with software MoPro Viewer [START_REF] Guillot | MoProViewer: a molecule viewer for theMoProcharge-density analysis program[END_REF] According to the recent study by Eskandariet al. [START_REF] Eskandari | Halogen bonding: A lump-hole interaction[END_REF], halogen bonding can be described as an interaction between a depleted region also known as a -hole on the halogen atom and a charge concentrated region also known as a lump on the other molecule. One Cl•••Cl halogen bond occurs in the crystal (Fig. 3b) but this type of contact is under-represented, similarly to the Besides that, the oxygen atoms are mostly involved in weak C-H•••O hydrogen bonds which are slightly favoured at EO.Hc= 1.06. The nitrogen atom, despite its hydrogen bond acceptor potential, is mostly surrounded by carbon atoms. On the other hand, the H-N nitrogen atoms which forms an intramolecular hydrogen bonds with the C1=O2 carbonyl is only surrounded by Hc atoms. Globally, the molecular surface contains 86 % of hydrophobic atoms (C, Hc, Cl) and 14% of hydrogen bond donors/acceptors (Hp, O). When this partition of the molecule is considered, it can be derived from Table 3, that the contacts between hydrophilic contacts appear quite enriched at E=2. 6. The hydrophobic ones are slightly overrepresented at E=1.04 while cross contacts between polar/non-polar atoms are significantly disfavoured at E=0.76.

Molecular electron density

The quality of the high resolution data and the multipole model is attested by the residual density map for the phenyl ring in the molecule and is shown in Fig. 6.

The topological properties of the electron density are listed in Table 4 and in Table S2 of the Supporting information. The most relevant results concern the heteronuclear Cl-C and the C-N bonds. The C-Cl bonds show lower magnitudes of ρ(r) and  2 ρ(r) at the bond critical point. This reflects the weak shared interactions occurring between these atoms. Moreover, the bonding density peaks are closer to the carbon atom and their positions do not correspond to that of the CPs, resulting in weaker ρ(rcp) and  2 ρ(rcp) values. The electron density was found to be polarized towards the chlorine atom, bearing the negative charge. The corresponding experimental/theoretical electron density and Laplacian of electron density for the C14-Cl1 and C10-Cl2 bond critical points are 1.29/1.35eÅ -3 ; 1.31/1.38 eÅ -3 and -1.8/-7.2 eÅ -5 ;-2.2/-7.6 eÅ - 5 respectively, and this was well agreed with the earlier reported values [START_REF] Hathwar | Halogen bonding in 2,5-dichloro-1,4-benzoquinone: Insights from experimental and theoretical charge density analysis[END_REF]. We note that theoretical calculations predict the critical point closer to the region of negative Laplacian of the Chlorine atoms, at variance from the multipolar model adopted for the experimental electron density.

As anticipated in section 4.1, the C-N bonds display values which are larger than for typical single bonds and effects of the electron delocalization are visible. Anyway, the two bonds are not identical: C8-N1 [1.414( 6)Å] is comparatively longer than C9-N1 [1.394( 5)Å] which shows a weaker π-donating effect of the amine nitrogen atom towards the phenyl ring compared to the 2,6-dichlorophenylring.

As for the two C-O bonds, C1-O2 is stronger and it clearly shows the localized nature of the double bond, compared with C1-O1, a localized single bond. In fact, at the bond critical points, the electron density is 2.83e/Å 3 (Expt.) or 2.70 e/Å 3 (Theor.) for C1-O2 and 2.33 e/Å 3 (Expt.) or 2.17 e/Å 3 (Theor.) for C1-O1.

The deformation maps in the 2,6-dichlorophenyl (Fig. 7a) and phenyl ring (Fig. 7b) provide a pictorial view of the electron density depicting the bonding regions and electron lone pairs of the molecule. This information is in keeping with that provided by the Laplacian maps(Fig. 8a-8c).

Molecular electrostatics

The molecule features some regions of strongly negative electrostatic potentials and others dominated by a positive potential. Of course, this mirrors the distribution of charges in the molecule. The carboxylic group consists of carbon C1, which bears the highest positive QTAIM charge (+1.536e, see Table 5) and of the two electronegative O atoms O1 (-1.071e) and O2 (-1.245e). Comparatively, the O2 atom (O2=C) has slightly more negative charge than the O1 atom (C1-O1-H1B). The aminic N1 atom bears the second most negative charge and the electron rich region around the atom is clearly seen in the Fig. 9b. The group charges on the 2,6-dichlorophenyl ring, the other phenyl ring and the acid chain are calculated as +0.285e, +0.367e and +0.030e, respectively. The phenyl ring has an electropositive contribution on the molecule while the 2,6-dichlorophenyl ring is more electronegative. The Cl2 atom has slightly larger volume (109.6Å 3 ) than the Cl1 atom (102.1Å 3 ) (as listed in Table 5).

The electrostatic potential maps of the intermolecular interactions are shown in Fig. 9c in which, the high negative potential is observed around the vicinity of carboxyl group of the diclofenac molecule in both the crystal and TTR protein environment. From this, it was confirmed that the carboxyl group of diclofenac molecule plays an important role in the drugreceptor interactions.

Topological properties for Intermolecular interactions in crystal and active site phase

Topological properties of the strong intermolecular interactions in the crystal and protein environment were determined by using charge density analysis to understand the strength of intermolecular interactions at electronic level, which plays vital role in the molecular recognition of bio-systems [START_REF] Saravanan | Understanding the conformational flexibility and electrostatic properties of curcumin in the active site of rhAChE via molecular docking, molecular dynamics, and charge density analysis[END_REF][START_REF] Saravanan | Investigation of intermolecular interactions and stability of verubecestat in the active site of BACE1: Development of first model from QM/MM-based charge density and MD analysis[END_REF]. Numerous reports clearly explained about the topological properties of intermolecular interactions [START_REF] Arputharaj | Topological Electron Density Analysis and Electrostatic Properties of Aspirin: An Experimental and Theoretical Study[END_REF][START_REF] Rajalakshmi | Topological analysis of electron density and the electrostatic properties of isoniazid: an experimental and theoretical study[END_REF][START_REF] Zhurova | Charge Density and Electrostatic Potential Study of 16α,17β-Estriol and the Binding of Estrogen Molecules to the Estrogen Receptors ERα and ERβ[END_REF][START_REF] Wu | A charge density study of π-delocalization and intermolecular interactions[END_REF][START_REF] Hibbs | X-N Charge density analysis of the hydrogen bonding motif in 1-(2-hydroxy-5-nitrophenyl)ethanoneElectronic supplementary information (ESI) available: multipole population coefficients and pseudoatom parameterization[END_REF]. However, the same properties were calculated and compared for diclofenac molecule in the solid state and human transthyretin enzyme environment with the help of (3,-1) type of critical point search. In the crystal phase, the electron density and Laplacian of electron density of H(1B)•••O(2) interaction is 0.30 eÅ -3 and 5.19 eÅ -5 respectively. The Laplacian map of the interaction is shown in Fig. 10a,b; a lone pair of the oxygen atom is directed towards the hydrogen atom. The carboxyl group of the diclofenac molecule forms strong dimer interactions in the solid state and the same group forms strong interaction in the protein environment. In the human transthyretin-diclofenac complex, the electron density and Laplacian of electron density of O(2)•••H-O Ser117 interaction is 0.149 eÅ -3 and 2.113 eÅ -5 respectively and the interaction is shown in Fig. 10c. On comparing both the electron density and Laplacian of electron density of the particular interaction in the solid state and protein phase, the ρ(rbcp) and  2 ρ(rbcp) values are decreased when the molecule enters into the active site of human transthyretin enzyme. This is strongly attributed to the effect of the protein environment. Furthermore, the total energy densities of the intermolecular interactions were calculated from the sum of kinetic energy and potential energy as stated by viral theorem [START_REF] Bader | The characterization of atomic interactions[END_REF][START_REF] Abramov | On the Possibility of Kinetic Energy Density Evaluation from the Experimental Electron-Density Distribution[END_REF][START_REF] Espinosa | Hydrogen bond strengths revealed by topological analyses of experimentally observed electron densities[END_REF]. The topological properties [ρ(rbcp) and 2 ρ(rbcp)], the total energy density (Hr), potential energy density (Vr), kinetic energy density (Gr), bond dissociation energy (D), bond path length (Rij) and bond critical point distances (d1 and d2) of intermolecular interactions are presented in Table 6. The negative kinetic energy dominates the potential energy and, as a result, the total energies are negative.

Reactivity properties

The analysis on HOMO-LUMO plays a vital role in explaining the process of charge transfer within the molecule. The HOMO-LUMO energy gap characterizes the kinetic stability and chemical reactivity of the molecule. Some global reactivity descriptors such as electron affinity [START_REF] Sarkar | Toxicity analysis of polychlorinated dibenzofurans through global and local electrophilicities[END_REF], ionization potential, hardness [START_REF] Priyakumar | Theoretical Study of Silabenzene and Its Valence Isomers[END_REF][START_REF] Ghiasi | A computational study of the arsabenzenes: Structure, properties and aromaticity[END_REF], softness, electronegativity [START_REF]The nature of the chemical bond and the structure of molecules and crystals[END_REF] and electrophilicity index [START_REF] Parr | Electrophilicity Index[END_REF][START_REF] Hazarika | Molecular structure and reactivity of antituberculosis drug molecules isoniazid, pyrazinamide, and 2methylheptylisonicotinate: a density functional approach[END_REF] help us to characterize the nature and interaction of the molecule. The reactivity descriptors have been listed in Table 7. The low energy gap between the donor and acceptor shows that the transfer of energy takes place well between them. The very low electron affinity (0.031a.u.) and low ionization potential (0.213a.u.) says that the molecule has low kinetic stability as well as high chemical reactivity. The very low global hardness value (0.091a.u.) and very high softness value (5.495a.u.) pronounces that the molecule is ready to share electrons and it is termed as very soft molecule. The very low electronegativity (0.122a.u.)

predicts that the molecule has very low attraction for electrons. Moreover, the low electrophilicity index (0.086a.u.) tells that the molecule can act as nucleophile which could donate electrons while interacting with the bio-molecule.

Intermolecular interactions and stability of the complex

MD simulation is a trustful tool to validate the conformational changes, stability of human transthyretin-diclofenac complex and intermolecular interactions. The stable conformation of diclofenac was positioned initially in the active site, as in the PDB coordinates file.

There are few statistical parameters to evaluate the conformational stability and strength of intermolecular interactions such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg) and number of hydrogen bonding in the proteinligand complex. The RMSD of the transthyretin-diclofenac complex was calculated with respect to initial coordinates to understand the deviation of the backbone of the protein. The maximum RMSD of the complex is 1.5 Å, which indicates that the complex is stable during the 20 ns of MD simulation (Fig. 11 and Fig. 12). The RMSF of the transthyretin-diclofenac complex clearly indicates that larger fluctuation was noticed in the loop regions of the complex than in the α-helix and β-sheet domain. The maximum fluctuation of the complex during the simulation is ~3 Å in the loop region which is not stabilized by intermolecular interactions.

Particularly, the active site amino acids are not fluctuating highly; the RMSF of these amino acids is ~2 Å. The radius of gyration of the transthyretin-diclofenac complex is determined to understand the compactness of the protein when the ligand molecule enters into the active site of the protein. The average Rggyration radius value of the transthyretin-diclofenac complex is ~17.7663 Å (with standard deviation of 0.10848); which indicates that the complex is highly stable during the MD simulation. The number of hydrogen bonding interactions in the transthyretin-diclofenac complex comprise ~120 (within the protein) and 3 (between the protein and ligand) (Fig. 11d). Notably, the carboxyl group of the diclofenac molecule is forming strong interactions with active site residues which is clearly shown in the Fig. 12.

Apart from the direct hydrogen bonds, water mediated hydrogen bonds and hydrophobic clusters stabilize also the orientation. The H-C hydrogen atoms involved in the cluster formation belong to amino acids SER 108, THR 109, THR 110 and SER 221. The electrophilic and nucleophilic regions of the ligands interact with the nucleophilic and electrophilic regions of the protein active site, respectively. The Connolly surface view of the transthyretindiclofenac complex clearly shows the binding cavity. The superimposed view of the transthyretin-diclofenac complex also confirms its stability (Fig. 12). On comparing the crystal and simulated transthyretin-diclofenac complex, the relative orientation of the dichlorobenzene and phenyl ring of the diclofenac molecule is different; this can clearly be attributed to the effect of interactions and to the dynamical behaviour.

Conclusion

The experimental electron density analysis of the diclofenac molecule was performed at 110K using high-resolution X-ray diffraction data. The aminic N1 atom, which has three bonds, could be expected to have tetrahedral geometry (sp 3 ). The structure analysis shows that it is, instead, closer to a trigonal geometry (sp 2 ) but it still has a weak electron lone pair. This can be attributed to the involvement of N1 atom in resonance primarily with the aromatic 2,6dichlorophenyl ring. The C8-N1 bond is comparatively longer than the bond C9-N1 which shows a weaker π-donating effect of the amine nitrogen atom towards the phenyl ring than towards the 2,6-dichlorophenyl ring.

The interaction between two parallel aromatic rings results in a Cl…Cl distance shorter than the sum of van der Waals radii. The O2=C1 oxygen atom is considered to be the most nucleophilic and it is the preferred site of electrophilic attack. The drug proves as a potent drug and physiologically active according to the Lipinski rule of five parameters. The reactivity descriptors describes that diclofenac has high chemical reactivity, is very soft and is a fast interacting drug with the biomolecules. The molecular dynamics simulation, started from crystal structure of the complex, shows a clear picture of the interaction between the diclofenac molecule and the amino acid residues present in the human transthyretin protein.

Bonds D-H(Å) H•••A(Å) D•••A(Å) D-H•••A( o )

O1-H1B

  ), through O1-H1B•••O2 interactions. The packing of the molecule along b axis is given in Fig.2b.There are anyway other contacts, such as C-H•••Cl and C-H•••O, that also give contributions to the packing of the crystal. The Hirshfeld

  . The enrichment ratio, which compares the actual contacts with equi-probable theoretical contacts (derived from the chemical content on the surface), enables to see if interactions are over or under-represented. Hydrophobic hydrogen atoms Hc bound to C were distinguished from polar Hp atoms bound to O or N atoms. The four major contacts are hydrophobic, of type C•••Hc, C•••Cl, Cl•••Hc and Hc•••Hc (Table 3, Fig. 5). Among these four contacts, three are close to a balanced representation (E1) while Cl•••C contacts are significantly over-represented with EC.Cl=1.54 due to numerous C-Cl•••C interactions. Weak hydrogen bonds of type Cl•••Hc are also well represented.

  Cl•••O halogen bonds. The C•••C contacts are impoverished at E=0.68 due to the occurrence of very limited non-parallel cycle stacking. C•••Hc contacts are more represented as significant C-H•••π interactions are seen, notably between phenyl rings forming an angle of 85.1° and related by two-fold helical axis of direction b (Fig. 2b). The most enriched contact at E= 5.0 is between polar Hp (H-O and H-N) hydrogen atoms and oxygen, although it represents less than 5% of the contact surface. The O•••H-O strong double hydrogen bond constitutes the most attractive electrostatic interaction in the crystal packing.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .

 123456 Fig. 1. ORTEP view of the diclofenac molecule, showing the atom-numbering scheme.Displacement ellipsoids are drawn at the 50% probability level.

Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .Fig. 11 .

 7891011 Fig. 7. Deformation electron density map for a) 2, 6-dichlorophenyl and b) Phenyl ring. Contour intervals are 0.05 eÅ -3 ; blue lines are positive contours, red lines are negative contours and green lines are zero contours.

Fig. 12 .

 12 Fig.12. The diclofenac molecule at the binding site of human transthyretin enzyme in different views.

Figure Sup1 .

 Sup1 Figure Sup1.Normal probability plot of (F 2 o-F 2 c) using program DRKplot (Zhurov, Zhurova& Pinkerton, J. Appl. Cryst. (2008) 41, 340-349). Weighting scheme used: W= 1/ [σ(F 2 )] 2 .

Figure Sup2 .

 Sup2 Figure Sup2. Variation of scale factor, Σ (Fo 2 )/ Σ (Fc 2 ) with respect to resolution, obtained using program DRKplot (Zhurov, Zhurova& Pinkerton, J. Appl. Cryst. (2008) 41, 340-349).

Figure Sup3 .

 Sup3 Figure Sup3. Fractal dimension (d f ) vs. residual electron density ρ(Fo-Fc) using all reflections. (Meindl, K & Henn, J., Acta Cryst. (2008). A64, 404-418). sin max / is 1.04 Å -1 .

Table 1

 1 Crystallographic and refinement details.

	Crystal data
	Chemical formula	C14H11Cl2NO2
	Mr	296.14
	Crystal system, space group	Monoclinic, C2/c
	Temperature (K)	110
	a, b, c (Å)	20.0817 (2), 6.8873(4), 19.8429 (2)
	β (°) V (Å 3 )	109.749 (1) 2583.02(4)
	Z	8
	Radiation type µ (mm -1 )	Mo Kα 0.50
	Crystal size (mm)	0.12 × 0.14 × 0.19
	Data collection
	Diffractometer	Oxford Diffraction SuperNova
	Absorption correction	Multiscan(Abspack, CrysAlisPro 1.171.38.43)
	Tmin, Tmax	0.705, 1.00
	No. measured;unique;	109947; 12271;9581
	observed [I> 2σ(I)] reflections (sin θ/λ)max (Å -1 )	1.043
	Refinement R[F 2 > 2σ(F 2 )], wR(F 2 ), S	0.034, 0.054, 0.982
	No. of reflections	11436
	No. of parameters	609
	H-atom treatment	H atoms treated by a mixture of independent and
	Δρmax, Δρmin (e Å -3 )	constrained refinement 0.41, -0.41

Table 2

 2 Hydrogen bonding details.

Table 3

 3 Nature of intermolecular contacts on the Hirshfeld surface by chemical types. The second row contains the contribution S X of each chemical type X on the Hirshfeld. The third part of the Table show the %C xy of the contact types on the surface. The lower part of the Table shows the E xy enrichment ratios of contact types. The major C xy contact types and the E xy ratios larger than unity (corresponding to the significantly enriched contacts) are highlighted in bold characters. The hydrophobic Hc atoms bound to carbon were distinguished from the more polar Hp hydrogen atom from the COOH and NH groups. Chemical types have been regrouped in hydrophobic (C, Hc, Cl) and charged (Ho, O) atoms. N chemical type is omitted as it represents only 0.5% on the contact surface.

	chem.	C	Hc	Cl	Hp	O
	surf. %	29.5	35.1	20.8	5.6	8.3
	C	5.7
	Hc	21.8	12.6		% contacts
	Cl	18.8	14.7	2.3
	Hp	1.3	2.6	2.5	0.2
	O	3.6	6.2	1.9	4.8	0.0
	C	0.68
	Hc	1.07	1.02		enrichment
	Cl	1.54	0.98	0.50
	Hp	0.37	0.63	1.03	0.61
	O	0.75	1.06	0.53	5.0	0.00

Table 4

 4 Topological properties of the electron density at the bond critical points of Diclofenac.

	Bonds	ρbcp (e/Å 3 )	 2 ρbcp (e/Å 5 )	ε	d1 (Å)	d2 (Å)
	(C-C) aro	2.1(2)	-18.4(6)	0.198(2) 0.697 0.617
	C2-C3	1.65(2)	-9.4(5)	0.0(2) 0.728 0.785
	C14-Cl1	1.29(1)	-1.8(4)	0.11(5)	0.96 0.777
	C10-Cl2	1.31(2)	-2.2(5)	0.02(5) 0.956 0.769
	C1-O1	2.33(4)	-27.7(2)	0.09(3) 0.813 0.495
	C1-O2	2.83(8)	-22.2(7)	0.18(4) 0.814 0.417
	C8-N1	2.00(2)	-16.6(1)	0.12(2) 0.822 0.593
	C9-N1	2.13(3)	-15.7(1)	0.14(4)	0.78 0.614
	(C-H) aro	1.82(2)	-16.8(7)	0.04(2) 0.723 0.359
	(C-H) sp3	1.74(0)	-13.9(0)	0.01(1) 0.716 0.375
	O1-H1B	2.10(4)	-23.7(3)	0.003(4) 0.754 0.226
	N1-H1A	2.03(4)	-21.2(1)	0.05(2) 0.756 0.259

Table 5

 5 Experimental integrated topological atomic charges and volumes.

	Atom	Q(e)	Atom	Q(e)
	Cl1	-0.157	C11	-0.022
	Cl2	-0.177	C12	-0.124
	O1	-1.071	C13	-0.024
	O2	-1.245	C14	-0.009
	N1	-1.111	H2A	0.132
	C1	1.536	H2B	0.132
	C2	-0.048	H4	0.147
	C3	-0.146	H5	0.129
	C4	-0.218	H6	0.129
	C5	-0.009	H7	0.189
	C6	-0.121	H11	0.148
	C7	0.013	H12	0.14
	C8	0.253	H13	0.16
	C9	0.28	H1A	0.444
	C10	0.053	H1B	0.597

Table 6 .

 6 Topological properties of electron density for the main intermolecular interactions in the crystal phase and active site region. G and V are the kinetic, potential energy respectively.

	Bond	cp(rcp)	 2 (rcp)	Ellip.	G(rcp)	V(rcp)	d1	d2
		e/Å 3	e/Å 5		au	au	(Å)	(Å)
	O1-H(1B)•••O2 #1	0.301(3)	5.19(5)	0.03(2) 0.052 -0.050 1.649 1.114
	C2-H(2A)•••O2 #2	0.033(2)	0.64(2)	0.11(1) 0.005 -0.003 1.554 1.023
	C6-H6•••Cl1 #3	0.038(2)	0.66(1)	4.33(5) 0.005 -0.003 1.812 1.034
	O2•••H-O Ser117	0.149	2.113	0.07	0.03	-0.04	0.756 1.248
	#1. -X+1/2 ; -Y+5/2 ; -Z+1 #2.-X+1/2 ; -Y+3/2 ; -Z+1, #3. -X+1 ; Y ; -Z+3/2

Table 7 .

 7 Reactivity descriptors of the Diclofenac molecule.

	Molecular descriptor Energy(a.u.)
	Electron affinity
	A=[-E LUMO ]	0.031
	Ionization potential
	I=[-E HOMO ]	0.213
	Global hardness
	Η=(I-A)/2	0.091
	Softness
	S= ½ η	5.495
	Electronegativity
	Χ=(I+A)/2	0.122
	Electrophilicity index
	Ω=μ 2 /2η	0.086
	HOMO energy	-0.213
	LUMO energy	-0.031

  The first line indicates the experimental values; the second line indicates the values obtained from the gas phase theoretical calculation. ρ BCP (r), ∇ 𝜌 𝐫 , electron density and the Laplacian of electron density at the BCPs; ε, bond ellipticity; d 1 , d 2 , the distance between BCP and each bonded atom.

	Bonds C4-H4	ρBCP (e/Å 3 ) 1.85(3) 2.63	𝜵 𝟐 𝝆 ( e/Å 5 ) -18.1(7) -42.6	d1 0.05(1) 0.723 ε (Å) 0.001 0.545	d2 0.360 (Å) 0.367
	C2-C3 C5-H5	1.65(2) 1.81(2)	-9.4(8) -15.8(7)	0.00(2) 0.728 0.02(1) 0.708	0.785 0.374
		1.66 2.63	-13.4 -42.5	0.04 0.004	0.758 0.546	0.755 0.367
	C3-C4 C6-H6	2.18(2) 1.85(2)	-18.9(7) -17.6(8)	0.26(2) 0.686 0.05(2) 0.718	0.707 0.365
		2.09 2.64	-20.7 -42.9	0.22 0.003	0.685 0.366	0.707 0.547
	C3-C8 C7-H7	2.11(2) 1.78(2)	-17.7(7) -16.5(7)	0.29(2) 0.693 0.04(2) 0.740	0.714 0.343
		2.06 2.63	-20.2 -42.7	0.24 0.008	0.717 0.364	0.69 0.549
	C4-C5 C11-H11	2.17(2) 1.82(1)	-18.8(5) -16.1(6)	0.19(4) 0.708 0.03(2) 0.727	0.685 0.356
		2.09 2.64	-20.9 -43.2	0.19 0.008	0.689 0.552	0.702 0.36
	C5-C6 C13-H13	2.17(2) 1.81(2)	-19.3(7) -16.1(6)	0.17(2) 0.685 0.01(1) 0.729	0.360 0.354
		2.11 2.64	-21.2 -43.2	0.19 0.008	0.697 0.359	0.691 0.553
	C6-C7 C12-H12	1.81(2) 1.83(2)	-15.8(8) -17.7(7)	0.02(2) 0.708 0.06(2) 0.719	0.374 0.364
		2.08 2.63	-20.7 -42.8	0.20 0.003	0.700 0.363	0.692 0.549
	C7-C8 C2-H2A	1.85(3) 1.72(2)	-17.6(9) -13.7(5)	0.05(2) 0.718 0.01(1) 0.716	0.365 0.376
		2.08 2.42	-20.6 -36.4	0.22 0.01	0.673 0.358	0.724 0.596
	C9-C10 C2-H2B	2.10(1) 1.76(3)	-17.5(4) -14.1(6)	0.27(2) 0.698 0.01(1) 0.716	0.709 0.375
		2.05 2.42	-19.8 -36.5	0.27 0.012	0.700 0.591	0.707 0.363
	C9-C14 O1-H1B	2.14(2) 2.10(4)	-17.9(5) -23.7(3)	0.27(3) 0.697 0.003(4) 0.754	0.709 0.226
		2.06 2.68	-19.9 -70.7	0.27 0.02	0.707 0.151	0.698 0.755
	C10-C11 N1-H1A	2.18(3) 2.03(4)	-19.7(8) -21.2(1)	0.23(2) 0.709 0.05(2) 0.756	0.681 0.259
		2.09 3.95	-20.8 -114.2	0.23 0.04	0.663 0.148	0.727 0.626
	C11-C12	2.16(2)	-19.397) 0.21(1) 0.680	0.710
		2.10	-21.1	0.19	0.688	0.700
	C12-C13	2.15(2)	-19.1(7)	0.19(1) 0.711	0.679
		2.10	-21.1	0.20	0.700	0.688
	C13-C14	2.18(2)	-19.4(7)	0.23(2) 0.674	0.716
		2.10	-20.8	0.23	0.663	0.728
	C14-Cl1	1.29(1)	-1.8(4)	0.11(5) 0.960	0.777
		1.35	-7.2	0.07	0.995	0.743
	C10-Cl2	1.31(2)	-2.2(5)	0.02(5) 0.956	0.769
		1.38	-7.6	0.07	0.737	0.989
	C1-O1	2.33(4)	-27.7(2)	0.09(3) 0.813	0.495
		2.17	-7.6	0.04	0.444	0.867
	C1-O2	2.83(8)	-22.2(7)	0.18(4) 0.814	0.417
		2.70	-9.0	0.08	0.807	0.423
	C8-N1	2.00(2)	-16.6(1)	0.12(2) 0.822	0.593
		1.92	-19.3	0.08	0.551	0.865
	C9-N1	2.13(3)	-15.7(1)	0.14(4) 0.780	0.614
		2.02	-21.1	0.11	0.548	0.848

Supporting information