Thermodynamic Integration Research Articles

Molecular simulations of understanding the Na+ ion structure, dynamic and thermodynamic behavior in ionic liquids: Butyl ammonium hydrogen bisulfate and tri-butyl ammonium hydrogen bisulfate

This manuscript presents the all-atom molecular dynamics simulations to investigate intermolecular structure and solvation thermodynamics of Na+ ion in two different ammonium-based protic ionic liquids (1) Butyl Ammonium hydrogen bisulfate [BA+][HSO4−], (2) Tri-butyl ammonium hydrogen bisulfate [TBA+][HSO4−]. The ionic liquid [BA+][HSO4−] show a more coordinated behavior when compared to [TBA+][HSO4−], which is observed over the temperature range from 278 K to 348 K. Hydrogens of the cations show a hydrogen bonding interaction with oxygens of anions. The cationic [TBA+] molecules show more solvation behavior with anions when compared to the [BA+]. The Na+ ion show a strong coordination structure with [HSO4−] in [TBA+][HSO4−] when compared to the [BA+][HSO4−]. We further calculate the detailed solvation free energy (ΔG) calculations using thermodynamic integration. We found that the ΔG of Na+ is more favorable in [TBA+][HSO4−] when compared to [BA+][HSO4−] in the temperature range varying from 278 K to 348 K. With the temperature rise, we observe the more favorable solvation of Na+ in both ionic liquids. On the other hand, the solvation of Cl− becomes less favorable. Overall, this manuscript provides detailed molecular level structural and thermodynamic origins of Na+ in protic ionic liquids useful for designing and developing sustainable electrolytes for Na+ battery applications.

Application of referenced thermodynamic integration to Bayesian model selection.

Evaluating normalising constants is important across a range of topics in statistical learning, notably Bayesian model selection. However, in many realistic problems this involves the integration of analytically intractable, high-dimensional distributions, and therefore requires the use of stochastic methods such as thermodynamic integration (TI). In this paper we apply a simple but under-appreciated variation of the TI method, here referred to as referenced TI, which computes a single model's normalising constant in an efficient way by using a judiciously chosen reference density. The advantages of the approach and theoretical considerations are set out, along with pedagogical 1 and 2D examples. The approach is shown to be useful in practice when applied to a real problem -to perform model selection for a semi-mechanistic hierarchical Bayesian model of COVID-19 transmission in South Korea involving the integration of a 200D density.

Assessing the Performance of Non-Equilibrium Thermodynamic Integration in Flavodoxin Redox Potential Estimation.

Flavodoxins are enzymes that contain the redox-active flavin mononucleotide (FMN) cofactor and play a crucial role in numerous biological processes, including energy conversion and electron transfer. Since the redox characteristics of flavodoxins are significantly impacted by the molecular environment of the FMN cofactor, the evaluation of the interplay between the redox properties of the flavin cofactor and its molecular surroundings in flavoproteins is a critical area of investigation for both fundamental research and technological advancements, as the electrochemical tuning of flavoproteins is necessary for optimal interaction with redox acceptor or donor molecules. In order to facilitate the rational design of biomolecular devices, it is imperative to have access to computational tools that can accurately predict the redox potential of both natural and artificial flavoproteins. In this study, we have investigated the feasibility of using non-equilibrium thermodynamic integration protocols to reliably predict the redox potential of flavodoxins. Using as a test set the wild-type flavodoxin from Clostridium Beijerinckii and eight experimentally characterized single-point mutants, we have computed their redox potential. Our results show that 75% (6 out of 8) of the calculated reaction free energies are within 1 kcal/mol of the experimental values, and none exceed an error of 2 kcal/mol, confirming that non-equilibrium thermodynamic integration is a trustworthy tool for the quantitative estimation of the redox potential of this biologically and technologically significant class of enzymes.

Anharmonicity in grain boundary energy for Al: Thermodynamic integration with artificial-neural-network potential

Anharmonicity in grain boundary (GB) energies for Al is elucidated by employing a thermodynamic integration method with a high-accuracy artificial-neural-network potential. Symmetric tilt GBs with the [001] and [11¯0] rotational axes are examined. It is found that anharmonic components tend to be small at low temperatures and increasingly contribute to reduction of GB free energies at elevated temperatures. The magnitude of anharmonicity is roughly correlated with excess volume at GBs. This is presumably because GB atoms tend to have longer bond lengths than the bulk atom, and hence the inclination of the potential energy surface at GBs is more moderate than in the bulk. On the other hand, the Σ13(510) GB takes non-negligible values of the anharmonic component even at intermediate temperatures. This GB contains a distinct atom around which the potential-energy curve significantly deviates from that of the harmonic approximation. Similar trends are also observed for other [001] GBs.

Exact tunneling splittings from symmetrized path integrals.

We develop a new simulation technique based on path-integral molecular dynamics for calculating ground-state tunneling splitting patterns from ratios of symmetrized partition functions. In particular, molecular systems are rigorously projected onto their J = 0 rotational state by an "Eckart spring" that connects two adjacent beads in a ring polymer. Using this procedure, the tunneling splitting can be obtained from thermodynamic integration at just one (sufficiently low) temperature. Converged results are formally identical to the values that would have been obtained by solving the full rovibrational Schrödinger equation on a given Born-Oppenheimer potential energy surface. The new approach is showcased with simulations of hydronium and methanol, which are in good agreement with wavefunction-based calculations and experimental measurements. The method will be of particular use for the study of low-barrier methyl rotations and other floppy modes, where instanton theory is not valid.

Solubility of NaCl under anisotropic stress state.

Salt solubility is generally determined under isotropic stress conditions. Yet, in the context of salt weathering of porous media, mechanical constraints on the in-pore growth of salt crystals are likely to be orientation-dependent, resulting in an anisotropic stress state on the crystal. In this paper, we determine by molecular simulation the solubility of NaCl in water when the crystal is subjected to anisotropic stress. Such anisotropy causes the chemical potential of the crystal to be orientation-dependent, and proper thermodynamic formulation requires describing the chemical potential as a tensor. The solute and crystal chemical potentials are computed from free energy calculations using Hamiltonian thermodynamic integration, and the usual condition of solubility is reformulated to account for the tensorial nature of the crystal chemical potential. We investigate in detail how the uniaxial compression of the crystal affects its solubility. The molecular simulation results led to revisiting the Correns law under anisotropic stress. Regarding the solute, the non-ideal behavior of the liquid phase is captured using Pitzer's ion interaction approach up to high concentrations of interest for in-pore crystallization and beyond the concentrations addressed in the existing literature. Regarding NaCl crystals, the validity of the generalized Gibbs-Duhem equation for a tensorial chemical potential is carefully verified, and it is found that crystallization progresses almost orthogonally to the crystal surface even under high shear stresses. Comparing uniaxial and isotropic compression highlights the major differences in solubility caused by stress anisotropy, and the revisited Correns law offers an appropriate framework to capture this phenomenon.

Thermodynamic stability of cytosine tetramers mediated by partially reduced silver trimers: QM/MM free energy analysis

AbstractPartially reduced silver clusters growing in a polycytosine matrix could simultaneously possess stabilizing ability and remarkable optical features. The mechanisms of this phenomenon are not fully understood so far. We examined the thermodynamic stability of a cytosine–silver associate comprising four cytosines and three silver atoms that can be considered the smallest known luminescent silver marker stabilized on DNA. Quantum and molecular mechanical representations thermodynamic integration toward the dissociation pathway of the cytosine tetramer, time‐dependent density functional theory calculations of the lowest energy excited states and quantum theory of atoms in molecules analysis of bonding showed that the reduction of two ions within the silver trimer mediating cytosine tetramer leads to the appearance of the green‐emitting silver cluster with low stabilization ability of cytosines relative to a tetramer stabilized by two silver ions.

Entropy-driven phases at high coverage adsorption of straight rigid rods on three-dimensional cubic lattices.

Combining Monte Carlo simulations and thermodynamic integration method, we study the configurational entropy per site of straight rigid rods of length k (k-mers) adsorbed on three-dimensional (3D) simple cubic lattices. The process is monitored by following the dependence of the lattice coverage θ on the chemical potential μ (adsorption isotherm). Then, we perform the integration of μ(θ) over θ to calculate the configurational entropy per site of the adsorbed phase s(k,θ) as a function of the coverage. Based on the behavior of the function s(k,θ), different phase diagrams are obtained according to the k values: k≤4, disordered phase; k=5,6, disordered and layered-disordered phases; and k≥7, disordered, nematic and layered-disordered phases. In the limit of θ→1 (full coverage), the configurational entropy per site is determined for values of k ranging between 2 and 8. For k≥6, MC data coincide (within the statistical uncertainty) with recent analytical predictions [D. Dhar and R. Rajesh, Phys. Rev. E 103, 042130 (2021)2470-004510.1103/PhysRevE.103.042130] for very large rods. This finding represents the first numerical validation of the expression obtained by Dhar and Rajesh for d-dimensional lattices with d>2. In addition, for k≥5, the values of s(k,θ→1) for simple cubic lattices are coincident with those values reported in [P. M. Pasinetti et al., Phys. Rev. E 104, 054136 (2021)2470-004510.1103/PhysRevE.104.054136] for two-dimensional (2D) square lattices. This is consistent with the picture that at high densities and k≥5, the layered-disordered phase is formed on the lattice. Under these conditions, the system breaks to 2D layers, and the adsorbed phase becomes essentially 2D. The 2D behavior of the fully covered lattice reinforces the conjecture that the large-k behavior of entropy per site is superuniversal, and holds on d-dimensional hypercubical lattices for all d≥2.

Different metal coordination in sub- and super-critical fluids: Do molybdenum(IV) chloride complexes contribute to mass transfer in magmatic systems?

Molybdenum (Mo) is extracted predominantly from magmatic-hydrothermal ore deposits, and knowledge of the thermodynamic properties of aqueous Mo complexes over a wide range in pressures (P) and temperatures (T) is required to predict the transport and deposition of Mo in ore-forming environments. These properties are usually derived from experimental studies conducted over a limited P-T range, with uncertainties increasing as temperature and pressure increase, because of inherent technical difficulties and decreasing number of studies. There is a discrepancy between numerous experimental studies at T < 400 °C and a handful of studies at magmatic-hydrothermal conditions. The low T (<400 °C) experiments suggested that Mo(VI) chloride complexes are unlikely to play a significant role in nature, whereas some studies in magmatic-hydrothermal conditions showed the predominance of Mo(VI) chloride complexes. In an attempt to resolve these inconsistencies, we investigated Mo(VI)-Cl complexations from 350 to 750 °C (500 to 2000 bar) using ab initio molecular dynamics (MD) simulations. The simulations reveal that in high-temperature magmatic-hydrothermal fluids (750 °C, 2000 bar), a tetrahedral Mo(VI) chloride complex, MoO2(OH)Cl(aq), becomes stable; this coordination differs from the octahedral Mo(VI)-chloride complexes identified in highly acidic sub-critical fluids. The ab initio thermodynamic integration method was employed to calculate the equilibrium constants of the reaction.MoO2(OH)Cl(aq) + H2O = MoO2(OH)2(aq) + HCl(aq)Thermodynamic modelling of the speciation using the calculated equilibrium constants shows that at lower temperatures (350 °C, 500 bar), Mo(VI)-chloro complexes are only important at highly acidic conditions (pH < 1), in accordance with previous studies. However, the MD results show that at higher temperatures, Mo-chloro-complexes may be important in transporting Mo in mildly acidic conditions up to the quartz-muscovite-K-feldspar buffered pH. Hence, Mo(VI) chloro-complexes could play an important role in the formation of magmatic-hydrothermal Mo deposits (e.g., Climax-type Mo deposits). We conclude that the mode of metal transport can change dramatically at high temperatures as a result of coordination changes caused by changes in the solvation properties of water; and that quantitative MD is an important tool to support the interpretation of scant experimental results available for magmatic-hydrothermal systems.

A Consistent Thermodynamic Characterization of the Adsorption Process through the Calculation of Free Energy Contributions.

We apply in this study different methodologies based on thermodynamic integration (TI), free energy perturbation (FEP), and potential of mean force (PMF) to address the challenging issue of the calculation of the free energy of adsorption. A model system composed of a solid substrate, an adsorbate, and solvent particles is specifically designed to reduce the dependence of our free energy results on the sampling of the phase space and the choice of the pathway. The reliability and efficiency of these alchemical free energy simulations are established through the closure of a thermodynamic cycle describing the adsorption process in solution and in a vacuum. We complete this study by the calculation of free energy contributions related to phenomena of desorption of solvent molecules and desolvation of the adsorbate upon adsorption. This calculation relies on the work of adhesion, the interfacial tension of the liquid-vapor of the solvent, and the free energy of solvation of the substrate. The different ways of calculating the free energy of adsorption are in excellent agreement and could complete experiments in the field of adsorption by giving quantitative data on the different energy contributions involved in the process.

An optimized thermodynamics integration protocol for identifying beneficial mutations in antibody design.

Accurate identification of beneficial mutations is central to antibody design. Many knowledge-based (KB) computational approaches have been developed to predict beneficial mutations, but their accuracy leaves room for improvement. Thermodynamic integration (TI) is an alchemical free energy algorithm that offers an alternative technique for identifying beneficial mutations, but its performance has not been evaluated. In this study, we developed an efficient TI protocol with high accuracy for predicting binding free energy changes of antibody mutations. The improved TI method outperforms KB methods at identifying both beneficial and deleterious mutations. We observed that KB methods have higher accuracies in predicting deleterious mutations than beneficial mutations. A pipeline using KB methods to efficiently exclude deleterious mutations and TI to accurately identify beneficial mutations was developed for high-throughput mutation scanning. The pipeline was applied to optimize the binding affinity of a broadly sarbecovirus neutralizing antibody 10-40 against the circulating severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) omicron variant. Three identified beneficial mutations show strong synergy and improve both binding affinity and neutralization potency of antibody 10-40. Molecular dynamics simulation revealed that the three mutations improve the binding affinity of antibody 10-40 through the stabilization of an altered binding mode with increased polar and hydrophobic interactions. Above all, this study presents an accurate and efficient TI-based approach for optimizing antibodies and other biomolecules.

Alchemical Free Energy and Hamiltonian Replica Exchange Molecular Dynamics to Compute Hydrofluorocarbon Isotherms in Imidazolium-Based Ionic Liquids.

Ionic liquids (ILs) have shown promise for applications that leverage differential gas solubility in an IL solvent, e.g., gas separations. Although most available literature provides Henry's law constants, the ability to efficiently estimate full isotherms is important for engineering design calculations. Molecular simulation can be used as a tool to predict full isotherms of gas in ILs. However, particle insertions or deletions in a charge-dense IL medium and the sluggish conformational dynamics of ILs present two sampling challenges for these systems. We therefore devised a method that uses Hamiltonian replica exchange (HREX) molecular dynamics (MD) combined with alchemical free energy calculations to compute full solubility isotherms of two different hydrofluorocarbons (HFCs) in imidazolium-based IL binary mixtures. This workflow is significantly faster than the Gibbs ensemble Monte Carlo (GEMC) simulations which fail to deal with the slow conformational relaxation caused by the sluggish dynamics of ILs. Multiple free energy estimators, including thermodynamic integration, free energy perturbation, and multistate Bennett acceptance ratio method, provided consistent results. Overall, the simulated Henry's law constant, isotherm curvature, and solubility trends match experimental results reasonably well. We close by calculating the full solubility isotherms of two HFCs in IL mixtures that have not been reported in the literature, demonstrating the potential of this method to be used for solubility prediction and setting the stage for future computational screening studies that search for the "best" IL to separate azeotropic HFC mixtures.

Co-segregation of Y and Zr in W-Cr-Y-Zr alloys: First-principles modeling at finite temperature and application to SMART materials

Spinodal phase separation in SMART (Self-passivating Metal Alloys with Reduced Thermo-oxidation) materials based on binary W-Cr with alloying elements Y and Zr is systematically investigated by a combination of Density Functional Theory with Cluster Expansion Hamiltonian and large-scale Monte Carlo simulations with thermodynamic integration. Comparing alloys of Zr with those of Y in ternary compositions with W and Cr, it is shown that there is a significant difference in the average short-range order between W and Zr, which changes from positive to negative from 500 K to higher temperatures in W70Cr29Zr1 alloys compared with the positive SRO (Short-Range Order) between W and Y in W70Cr29Y1 alloys at all temperatures. This change, however, does not affect the segregation behavior between W and Cr in the whole range of temperature between these alloys. Importantly, it is found that co-segregation of Y and Zr due to the negative SRO between them, in W70Cr29Y0.5Zr0.5 and W70Cr28Y1Zr1 considered alloys, increases the W-Cr positive SRO closer to the binary W70Cr30. Our modeling results reveal a significant impact on spinodal phase segregation of W and Cr in designing SMART materials for DEMO devices of future fusion power plants. Findings from recent experimental studies validate the modeling results, particularly the high-temperature behavior of the W-Cr-Zr alloys and the microstructure evolution of the W-Cr-Y-Zr alloys. The results provide a fundamental understanding of the radiation resistance phenomena observed in recent experiments on neutron-irradiated self-passivating alloys.

Efficient Training of the Machine-Learning Interatomic Potential Based on an Artificial Neural Network for Estimating the Helmholtz Free Energy of Alkali Metals

Although the Helmholtz free energy can be accurately calculated using the thermodynamic integration (TI) method based on first-principles molecular dynamics (FPMD) simulations, the computational cost of this method is very high. By training the data obtained by FPMD simulations based on an artificial neural network (ANN), an interatomic potential with a very low computational cost can be constructed. To learn the FPMD data efficiently, it is important to elucidate the relation between the number of training data and accuracy of the free energy obtained by the TI method. We statistically evaluated the minimum number of training data by constructing multiple ANN potentials to accurately estimate the Helmholtz free energy of rubidium. The root mean square errors (RMSEs) of energy, atomic force, and pressure between the ANN and FPMD converged as the number of training data increased. The radial distribution function of the FPMD can be reproduced, even if the RMSEs do not converge. However, the number of training data for which the RMSEs converge is required to calculate the free energy with sufficient accuracy. This study demonstrated that the quality of an ANN potential can be estimated by statistically investigating the dependence of the RMSEs on the number of training data.

Temperature Transferable and Thermodynamically Consistent Coarse-Grained Model for Binary Polymer Systems

For the simulations of polymer systems, coarse-grained (CG) models are often developed to tackle the length and time scale limitations that are not feasible by all-atom molecular dynamic simulations. However, due to the necessary simplification or reduction in atomistic degrees of freedom, CG models usually have poor transferability over various temperatures, compositions, or thermodynamic states. In this work, the structure-based iterative Boltzmann inversion (IBI) method is further developed to obtain temperature transferable and thermodynamically consistent CG potentials for binary polymer systems. In particular, the conventional IBI procedure is first utilized to optimize single-component CG potentials from homopolymer systems while cross-interaction potentials from a random block copolymer system. Afterward, a thermodynamic integration method is applied to refine the cross-interaction potential between two constituent components in binary systems until the Flory–Huggins parameter (χ) between two components is matched with experimental value. We apply the above optimization procedure for the polystyrene-block-poly(methyl methacrylate) (PS-b-PMMA) as an example, which is widely studied experimentally, and the χ value can be easily obtained in the literature. To validate the transferability of the obtained CG potential, systems of both PS and PMMA homopolymers, random block copolymers, diblock copolymers, and PS/PMMA blends with various compositions are simulated. Many properties, including the mass density, radius distribution function, and bonded distribution functions, between CG beads generated from CG simulations match very well with those from atomistic simulations in the simulated temperature range, i.e., from 453 to 500 K. More importantly, phase morphologies of diblock copolymers and phase diagram of PS/PMMA binary blends over wide temperature and compositional ranges generated from CG simulations using our model are in good agreement with experimental results and predictions from self-consistent field theory.

Accelerated Enveloping Distribution Sampling to Probe the Presence of Water Molecules.

Determining the presence of water molecules at protein-ligand interfaces is still a challenging task in free-energy calculations. The inappropriate placement of water molecules results in the stabilization of wrong conformational orientations of the ligand. With classical alchemical perturbation methods, such as thermodynamic integration (TI), it is essential to know the amount of water molecules in the active site of the respective ligands. However, the resolution of the crystal structure and the correct assignment of the electron density do not always lead to a clear placement of water molecules. In this work, we apply the one-step perturbation method named accelerated enveloping distribution sampling (AEDS) to determine the presence of water molecules in the active site by probing them in a fast and straightforward way. Based on these results, we combined the AEDS method with standard TI to calculate accurate binding free energies in the presence of buried water molecules. The main idea is to perturb the water molecules with AEDS such that they are allowed to alternate between regular water molecules and non-interacting dummy particles while treating the ligand with TI over an alchemical pathway. We demonstrate the use of AEDS to probe the presence of water molecules for six different test systems. For one of these, previous calculations showed difficulties to reproduce the experimental binding free energies, and here, we use the combined TI-AEDS approach to tackle these issues.

Modeling the high-pressure solid and liquid phases of tin from deep potentials with ab initio accuracy

Constructing an accurate atomistic model for the high-pressure phases of tin (Sn) is challenging because properties of Sn are sensitive to pressures. We develop machine-learning-based deep potentials for Sn with pressures ranging from 0 to 50 GPa and temperatures ranging from 0 to 2000 K. In particular, we find the deep potential, which is obtained by training the ab initio data from density functional theory calculations with the state-of-the-art SCAN exchange-correlation functional, is suitable to characterize high-pressure phases of Sn. We systematically validate several structural and elastic properties of the {\alpha} (diamond structure), {\beta}, bct, and bcc structures of Sn, as well as the structural and dynamic properties of liquid Sn. The thermodynamics integration method is further utilized to compute the free energies of the {\alpha}, {\beta}, bct, and liquid phases, from which the deep potential successfully predicts the phase diagram of Sn including the existence of the triple-point that qualitatively agrees with the experiment.

Study of SQ109 analogs binding to mycobacterium MmpL3 transporter using MD simulations and alchemical relative binding free energy calculations

N-geranyl-N΄-(2-adamantyl)ethane-1,2-diamine (SQ109) is a tuberculosis drug that has high potency against Mycobacterium tuberculosis (Mtb) and may function by blocking cell wall biosynthesis. After the crystal structure of MmpL3 from Mycobacterium smegmatis in complex with SQ109 became available, it was suggested that SQ109 inhibits Mmpl3 mycolic acid transporter. Here, we showed using molecular dynamics (MD) simulations that the binding profile of nine SQ109 analogs with inhibitory potency against Mtb and alkyl or aryl adducts at C-2 or C-1 adamantyl carbon to MmpL3 was consistent with the X-ray structure of MmpL3 – SQ109 complex. We showed that rotation of SQ109 around carbon–carbon bond in the monoprotonated ethylenediamine unit favors two gauche conformations as minima in water and lipophilic solvent using DFT calculations as well as inside the transporter’s binding area using MD simulations. The binding assays in micelles suggested that the binding affinity of the SQ109 analogs was increased for the larger, more hydrophobic adducts, which was consistent with our results from MD simulations of the SQ109 analogues suggesting that sizeable C-2 adamantyl adducts of SQ109 can fill a lipophilic region between Y257, Y646, F260 and F649 in MmpL3. This was confirmed quantitatively by our calculations of the relative binding free energies using the thermodynamic integration coupled with MD simulations method with a mean assigned error of 0.74 kcal mol−1 compared to the experimental values.Graphical abstract

Computing chemical potentials of adsorbed or confined fluids.

The chemical potential of adsorbed or confined fluids provides insight into their unique thermodynamic properties and determines adsorption isotherms. However, it is often difficult to compute this quantity from atomistic simulations using existing statistical mechanical methods. We introduce a computational framework that utilizes static structure factors, thermodynamic integration, and free energy perturbation for calculating the absolute chemical potential of fluids. For demonstration, we apply the method to compute the adsorption isotherms of carbon dioxide in a metal-organic framework and water in carbon nanotubes.

Target-aware Bayesian inference via generalized thermodynamic integration

In Bayesian inference, we are usually interested in the numerical approximation of integrals that are posterior expectations or marginal likelihoods (a.k.a., Bayesian evidence). In this paper, we focus on the computation of the posterior expectation of a function f(textbf{x}). We consider a target-aware scenario where f(textbf{x}) is known in advance and can be exploited in order to improve the estimation of the posterior expectation. In this scenario, this task can be reduced to perform several independent marginal likelihood estimation tasks. The idea of using a path of tempered posterior distributions has been widely applied in the literature for the computation of marginal likelihoods. Thermodynamic integration, path sampling and annealing importance sampling are well-known examples of algorithms belonging to this family of methods. In this work, we introduce a generalized thermodynamic integration (GTI) scheme which is able to perform a target-aware Bayesian inference, i.e., GTI can approximate the posterior expectation of a given function. Several scenarios of application of GTI are discussed and different numerical simulations are provided.