Development of a new force field for the family of primary aliphatic amines using the three steps systematic parameterization procedure

The applicability of the three steps systematic parametrization procedure (3SSPP) to develop a force field for primary amines was evaluated in the present work. Previous simulations of primary amines show that current force fields (FF) can underestimate some experimental values under room conditions. Therefore, we propose a new set of parameters, for an united atom (UA) model, that can be used for short and long amines which predict correctly thermodynamic and dynamical properties. Following the 3SSPP methodology, the partial charges are chosen to match the experimental dielectric constant whereas the Lennard-Jones (LJ) parameters, $\epsilon$ and $\sigma$, are fitted to reproduce the surface tension at the vapor-liquid interface and the liquid density, respectively. Simulations were initially conducted for the propylamine molecule by introducing three different types of carbon atoms, C$_\alpha$ and C$_\beta$, with electric charges, and C$_n$, without charge. Then, modifying the charges of the carbons and using the transferable LJ parameters, the new set of constants for long amines were found. The results show good agreement for the experimental dielectric constant and mass density with a percentage error less than 1%, whereas for the surface tension the error is up to 4%. For the short amines, methylamine and ethylamine, the new charges were obtained from a fitting function calculated from the long amines results. For these molecules, the values of the dielectric constant and the surface tension present errors of the order of 10% with the experimental data. Miscibility of the amines was also tested with the new parameters and the results show reasonable agreement with experiments.


Introduction
Amines are molecules derived from ammonia with one or more alkyl or aryl groups with great interest in areas of chemical engineering. They are used in several industrial applications, such as carbon dioxide retainers [1] and to remove hydrogen sulfides and carbon dioxide from natural gas [2,3]. Due to their biological activity they have also been used in drugs and medicines [4,5]. All of these attributes make amines very interesting molecules to study not only from the experimental but also from the theoretical and computational points of view.
Nowadays, computer simulations have become an important tool to study complex systems. In particular, molecular dynamics simulations are a good alternative to understanding and obtaining complementary information that can be difficult to collect from a laboratory. However, in order to have reliable computational results it is necessary to have good force fields that can reproduce several thermodynamic, dynamic and structural properties.
For the case of amines, several properties were predicted using different force fields. For instance, Rizzo et al. used an all atom (AA) OPLS model to report the liquid density and enthalpy of vaporization at one temperature showing reasonable agreement with experiments [6]. Wick et al. used a different AA-model, TraPPe-EH, to conduct simulations of amines at different temperatures above the boiling point [7].
A different approach to the AA model is to consider the united atom (UA) model where each chemical group is reduced to one single site with appropriate parameters. In particular, a few years ago the UA model with shifted Lennard-Jones centers (AUA) [8], and reparametrized time after (AUA4) [9], was used to study linear and branched hydrocarbons with good results to predict fluid densities and pressures [10][11][12][13]. Those models were also tested for alkylamides and alkanols giving partial agreement with actual experiments [14]. Brad and Patel, studied several thermodynamic and dynamic properties of methylamine and ethylamine using charge equilibration force fields [15] and they report discrepancies above 15% in some properties, e.g., enthalpy of vaporization.
On the other hand, Orozco et al. [16] proposed an anisotropic united atom force field (AUA4) with four charges to study phase equilibrium of several primary amines. They used displaced force centers for the Lennard-Jones potential and predicted thermodynamic and transport properties with acceptable results compared to experimental data. Later, the same authors using the same model, studied equilibrium and transport properties of primary, secondary and tertiary amines and they obtained good agreement with real experiments [17].
In the present study, we propose a simple united atom force field for primary amines using the 3SSPP method reported a few years ago [18]. The methodology was tested in several systems with very good results [19][20][21], and then we used the same methodology for the primary amine molecules.

Computational model
Primary amines were constructed using an amino group, NH 2 , with the nitrogen and hydrogen atoms explicitly modelled, attached to a hydrocarbon tail of united atoms for the CH 2 and CH 3 alkyl groups. In the united atom model, the first CH 2 group attached to the NH 2 group is named C , in the case of long amines the next CH 2 is named C , the rest CH 2 are C 2 and the last CH 3 group is called C 3 , see figure 1. In the present work, the amine force field considers inter-molecular and intra-molecular interactions, (2.1) The first three elements, bond , ang , dih , correspond to the intra-molecular interactions (bond, angle and torsional potentials, respectively) and the last two terms, LJ and coul , represent the inter-molecular interactions (Lennard-Jones and coulombic potentials, respectively). Several years ago, Rizzo and Jorgensen showed good hydrogen-bond strengths and hydration results for amines simulations using the OPLS force field [6]. Therefore, for the intra-molecular interactions we chose the OPLS model as initial force field. For the inter-molecular interactions we decided to use the TraPPe model used by Wick et al. [7] since they obtained good densities and critical temperatures, i.e., bulk properties, in their simulations of amines. On the other hand, the initial charges were obtained with the RESP (restricted electrostatics potential) method [22]. All simulations were conducted using the molecular dynamics method in the isothermal-isobaric (NPT) and canonical (NVT) ensembles. The dielectric constants, densities and diffusion coefficients were calculated in the NPT simulations with 500 amine molecules at constant temperature, = 298.15 K, and pressure, = 1 bar. The internal temperature was coupled to a Nosé-Hoover thermostat [23] with a relaxation time of = 1 ps, while the pressure was coupled to a Parrinello-Rahman barostat [24] with a time parameter of = 1.0 ps. The short range interactions were cutoff at 20 Å and simulations were performed up to 50 ns after 10 ns of equilibration.
The NVT simulations were used to calculate surface tensions in a box of dimensions = = 5.0 nm and = 15.0 nm with 1000 amine molecules in a liquid phase in the middle of the box, i.e., two liquidvapor interfaces were constructed. Simulations were conducted at constant temperature, = 298.15 K, using the Nosé-Hoover thermostat [23] with a relaxation time constant of = 1 ps. In this case, the short range interactions were cutoff at 25 Å to avoid any size effects on this property, as recommended in previous works [25]. Then, simulations run for 30 ns after 10 ns of equilibration.
The simulations to measure miscibility of amines in water were carried out in the NPT ensemble with 4000 water molecules and 100 amines, placed initially in an homogeneous mixture, in a box of initial dimensions of 5 × 5 × 5 nm. The pressure was = 1 bar and temperature = 298 K, with relaxation times of = 1 ps and = 1 ps, respectively. The short range interactions were cutoff at 20 Å and simulations were performed up to 50 ns after 10 ns of equilibration. For these simulations, we chose a water model which correctly reproduced several experimental data, including our target properties, the dielectric constant and the surface tension, TIP4P/ [26].
All simulations were carried out with GROMACS-2021 [27] software applying periodic boundary conditions in all directions using the leap-frog algorithm [28] with a time step of = 2 fs to solve the equations of motion with compressibility value of 4.5 × 10 −5 bar −1 . The unlike interactions between distinct atoms were calculated using the Lorentz-Berthelot combination rules [28] and the electrostatic interactions were handled with the particle mesh Ewald method [29] (fourth order with a Fourier spacing 0.16) whereas bond lengths were constrained using the Lincs algorithm [30].

Optimization procedure
The new force field was constructed following the three steps systematic parameterization procedure developed a few years ago [18]. In that method three experimental properties are chosen as target quantities to be fitted, the dielectric constant, the surface tension and the liquid density by scaling the charges and Lennard-Jones parameters of all atoms in the molecules.
In the first step of the 3SSPP scheme all the partial charges are scaled from the original ones to match the experimental dielectric constant. Then, using the NPT ensemble, simulations were conducted with a set of charges, and the static dielectric constant, , was estimated by the time average of the fluctuations of the dipolar moment, M, of the whole system, where B is Boltzmann's constant, is the absolute temperature, is the volume of the simulation cell, and M is the summation of the dipole moment vectors of all the atoms, where and r are the charge and position of atom . The angled brackets in equation (2.2) indicate a time average. For those calculations, partial charges were imposed in the nitrogen and hydrogens (of the NH 2 group), C and C whereas the C sites had zero charge.
The second step of the 3SSPP consists in parameterization of the LJ Lennard-Jones parameters to fit the experimental surface tension ( ) using the mechanical expression, The ⟨ ⟩ are the components of the stress pressure and is the length of the simulation box. The factor 1/2 is for the two liquid-vapor interfaces in the simulation box.

23603-3
In the third step, all the LJ Lennard-Jones parameters are scaled to obtain the experimental mass density.

Dielectric constant
Simulations started for the propylamine molecule. Firstly, the original partial charges for the nitrogen, hydrogen (of the NH 2 group) and C sites, (the C and C 3 charges were zero) were scaled until the experimental dielectric constant was fitted with an error less than 5%. In table 1 the new and original charges are shown. As the hydrocarbon chain increases (for long amines), the NH 2 group modifies its polar activity, i.e., a negative charge arises in the C carbon and consequently the charge of the C changes to keep the neutrality of the system. Then, the C and C charges of the long amines should be modified to fit the experimental dielectric constants.
In table 2, the values of C and C charges are shown for the -amines from three (propylamine) to ten (decylamine) carbons in the alkyl chains. The results for the dielectric constants are given in table 3 and figure 2. In the same figure 2 comparisons with other force field are included where it is observed that the values with the new parameters have an error less than 1% with the experiments. Table 2. Charges of the different carbons (C 3 , C 2 , C and C ) in the amine molecules for all the amines from 1 to 10 carbons in the hydrocarbon tail. Charges of the system with 2 carbons were obtained from equation (3.2).
Amines, num. of carbons In figure 3, a plot of the chain length (number of carbons) as function of the C and C carbon charges are shown along with the best fitting curves to the data. Several functions were tested to fit the data and the best one was a third order polynomial.
For the C , the function that best fits the data is, whereas for the C the fitting curves is, where , ( ) represent the number of carbons in the molecule.

Parameterization of Lennard-Jones parameters
Once the dielectric constant of the propylamine is obtained, the next step in the SSPP3 method is the evaluation of the surface tension. For that calculation, the charges are fixed and the simulations are carried out by scaling the LJ values of all the atoms in the molecule. However, the surface tension did not change significantly with any variation of the LJ , i.e., the best results were obtained with the original parameter.

23603-5
With the LJ parameters and charges, the LJ were now modified, i.e., LJ of all atoms in the molecule were scaled, in increments of 10%, and mass density was calculated until the experimental data were obtained, within an error less than 5%. The results for all the amines are given in table 3. Transferability of the new LJ and LJ was evaluated by simulating the rest of the long amines and considering the charges described in the previous section. Then, the surface tensions and densities were calculated using the LJ and LJ parameters of the carbons in the long amines. The surface tension results are plotted in figure 4 where it is possible to see that the new reparameterization gives slightly better agreement with the experiments [32] (ID: 6089, 6101, 7564, 7716, 7769, 7811, 7835, 7851, 15387, 8576), although in some cases the error is about 5% as in the butylamine (see table 3).  In figure 5, the mass density data are plotted and compared with actual experiments [33]. It can be observed that the new values are in much better agreement than the results with the old parameters and other force fields. In our case, the OPLS-AA data for the propylamine and butylamine differ by 4% from those reported in references [6,34].
With the above results, transferability of the force field was also tested for the shortest amines, the methylamine and ethylamine. The methylamine molecule has only one C 3 site and its charge was chosen to have a neutral charged molecule. In the case of the ethylamine, it has carbon C , where its charge was estimated from equation (3.2), and carbon C 3 , with an electric charge determined to have a molecule with zero total charge. In table 2, the values of the charges for the methylamine and ethylamine are shown. Using those charges and the LJ and LJ Lennard-Jones parameters, the target properties were calculated, i.e., the dielectric constant, the surface tension and the density. In the case of the dielectric constant we found errors of 5% and 0.86% with the experiments for the methylamine and ethylamine, respectively. For the densities, the errors were ≈ 1% and 6%, for the methylamine and ethylamine, respectively. For the surface tension, the errors were about 10% in both molecules, see table 3. Once the target properties were evaluated, the new parameters were also tested with other thermodynamic and dynamical quantities. In figure 6 and table 4, the heat capacities and enthalpies of vaporization were calculated where we observed that the results were well compared to experiments.  The heat capacity at constant pressure ( ) was calculated using the energy file obtained from the simulations in the NPT ensemble described in section 2. Then, the was evaluated from GROMACS utilities with the defaults settings and without the quantum corrections. For the enthalpy of vaporization, the following equation was used [35], is the potential energy, is the gas constant and is the temperature. Here, the enthalpy was evaluated with simulations of two boxes in the liquid and gas states. The liquid state was simulated as described in section 2, for the calculation of the dielectric constant. For the gas state, it was assumed that the gas is ideal and simulations were performed using a stochastic dynamics integrator in the canonical ensemble as implemented in the GROMACS software. The system contains one molecule and the simulations were run for 10 ns after an equilibration period of 1 ns. Then, the potential energy was obtained and then used to calculate the heat of vaporization.
The diffusion coefficients were also calculated and plotted in figure 7 and in table 4. They were evaluated with the Einstein relation, i.e., with the mean square displacements, considering the linear region of the plots [28].
The calculations of the diffusion coefficients were obtained over the entire production time using the trajectory file and the defaults settings of the GROMACS utilities for the separation of time origins. From table 4 and figure 7 we can see that the diffusion coefficients calculated in this work have significant errors compared to the experiments [36]. Nevertheless, they have better agreement with them than the other force fields reported in the literature [36].

Miscibility
From a previous work we know that miscibility of a solute in water can be used as another target property to evaluate the force field [40]. For the primary amines, we found that the short ones, up to the pentylamine, are miscible in water whereas the long ones are not [31]. Here, in the present paper, miscibility was studied in terms of density profiles, ( ), calculated from the NPT simulations in boxes with homogeneous water/amine mixtures as initial configurations, as described in section 2, where ( ) is the number of amine molecules in a volume of area and thickness Δ . is the factor to convert the number density to mass density. The density profiles were measured along the -axis using the GROMACS utility with 200 slices. Then, by analyzing the profiles of both components, water and amines, it was determined whether or not they were still mixing.
The miscibility was evaluated from the amount of mass of the solute, , dissolved in water, where is the mass of the solvent, ( ) and ( ) are the density profiles of the solute and water, respectively. Since the units are given in g/l then, for a litter of water, is equal to 1000 g, i.e, The results suggest that methylamine and ethylamine are miscible in water as can be seen in the snapshots of figures 8a and 8b, i.e., water and amines are mixed in the whole simulation box. By calculating the density profiles it is observed that water and the amines present uniform distributions along the simulation box indicating that both systems are well mixed (figures 9a and 9b) in agreement with the experiments. From equation (3.6), the values for the methylamine and ethylamine are 360 g/l and 471 g/l, respectively, which suggest miscibility as reported in the literature for those amines [31]. On the other hand, propylamine appears to be partial miscible in water, only a few amine molecules are located in the bulk water phase (figure 8c) whereas butylamine looks immiscible since two regions are well defined in the simulation box (figure 8d). The same conclusions can be stated from the density profiles. For the water/propylamine system, two regions of high and small densities are formed, suggesting that the two components are not completely mixed (figure 9c), i.e., there are regions rich in water or rich in amines.
In the case of the water/butylamine mixture, the system separates into two well defined bulk phases as observed in the density profiles of figure 9d, there are observed regions of nearly zero density for water and butylamine, i.e., the system forms a liquid/liquid interface. For long amines, the same trends are obtained as butylamine (not shown here), i.e., those amines are not miscible in water.

Conclusions
In the present work molecular dynamics simulations were conducted to construct a transferable force field for primary aliphatic amines. The new force field was built up following the 3SSPP method where charges and Lennard-Jones parameters were scaled to fit three target experimental properties, dielectric constant, surface tension and density.
Different force fields have been reported in the literature that calculate only a few thermodynamic or dynamical properties of amines. In some cases they are not good enough when those quantities are compared to actual experiments. Then, in this work we intend to generalize a model that predicts several experimental properties with good accuracy. The procedure to obtain the new parameters of the amine molecules is systematic and initially requires a reliable force field for the subsequent reparameterization.
The results show that the new force field adequately reproduces the target and other properties, although the shortest amines present slightly larger errors than the long ones compared to the experimental data. When miscibility is evaluated, the results are good for the shortest and largest amines, i.e., they show miscibility and no miscibility respectively, although for the intermediate amines the results are not good enough, since they show partial miscibility when they should be completely miscible in water. Finally, since the present force field shows good agreement with experiments it could be used as the starting one to characterize another type of amines, such as secondary or tertiary amines, as well as branched and/or aromatic amines, since they share an amino group with fixed partial charges.