Improving 1-propanol force field: a new methodology

A new force field for 1-propanol, in the united and all atom models, has been obtained by combining two different empirical methodologies. The first was developed by scaling atom charges and Lennard-Jones parameters to fit the dielectric constant, surface tension, and density; this methodology is named three steps systematic parameterization procedure (3SSPP), as reported by Pérez de la Luz et al. (J Chem Theory Comput 14:5949–5958, 2018). The second methodology consists of moving these parameters and together with the bond distance to obtain the liquid-vapor phase diagram of the CO2 molecule as discussed by Harris and Yung (J Phys Chem 99:12021–12024, 1995). The last methodology is used to obtain the self-diffusion coefficient, which was not consider in the 3SSPP. The 3SSPP/bond methodology is the 3SSPP plus the bond distance scaling. With this new methodology, the experimental density, dielectric constant, surface tension, and self-diffusion coefficient at ambient temperature could be achieved. Furthermore, we show the temperature dependence of the aforementioned properties. The static structure factors are in accordance with the experimental spectrum. Solubility is increased to the experimental value for the united atom (UA) model after applying this methodology and for all atom (AA) scheme, the experimental solubility value is maintained. Graphical abstract The reduction in bond distance of the 1-propanol molecule does not modify the structure factor. The reduction in bond distance of the 1-propanol molecule does not modify the structure factor.


Introduction
Alcohols are the most used solvents in science and engineering, with 1-propanol being the water-miscible compound with the largest aliphatic chain. Reproducing most of the thermodynamic, structural, and dynamic properties of any molecule is one of the main goals in molecular simulation. Many force fields reproduce some of these properties, mainly conformational energies, density, and the enthalpy of vaporization ( H v ) are reproduced using the all-atom optimized potentials for liquid simulations (OPLS/AA) force field [1][2][3]; H v , Gibbs free energy of solvation ( G solv ), and density are the target properties for the Groningen molecular simulation (GROMOS) force field [4]; free energies and many conformational, bulk, and interface properties in biological macromolecules are estimated from chemistry at Harvard macromolecular mechanics (CHARMM) [5]; and conformational energies and root-mean-square displacement to build molecular structures of proteins and nucleic acids for organic molecules are predicted by the general Amber force field (GAFF) [6].
The parameterization methods have been applied to reproduce the properties that current force fields do not consider. Cole et al. [7] used quantum mechanics to obtain partial charges and Lennard-Jones potential parameters (LJPP) to model organic and biomolecular systems, including proteins. One hundred forty-six organic liquids are benchmarked in the OPLS/AA, CHARMM general force field (CGenFF), and GAFF force fields by Caleman et al. [8]. The methodology 3SSPP was proposed to develop force fields [9], where intermolecular parameters were used to fit some target properties, such as the charge molecule distribution, the energy, and distance LJPP for obtaining dielectric constant, surface tension, and density, respectively. Recently, a methodology has been developed to improve the 3SSPP to force fields by assigning the Hirshfeld electronic density charge distribution in some polar liquids to obtain the experimental self-diffusion coefficient [10]. New force fields for water model, the simple point charge (SPC) [12] and a four-site rigid one, the TIP4P/2005 [13], were developed by considering that there is a dipole moment of minimum density at a fixed temperature to match the dielectric constant at room temperature and the temperature of maximum density.
On the other hand, Ríos-López et al. [14,15] parameterized tensoactives by scaling the distance and energy LJPP to obtain the micelle radius and surface tension, respectively. Different force field models have been proposed by modifying the force field of semi-flexible CO 2 molecules, such as scaling the charge and LJPP [16], and by changing the angle force constant and bond length [17], to obtain thermodynamics properties and the liquid-vapor coexistence curve.
With alcohols, a four-site methanol force field was developed by Martínez-Jiménez et al. [18]. This UA model has an extra site to match the dipole moment for methanol so the temperature dependence of dielectric constant and surface tension fit experimental data. It was proposed for mixtures, a linear scaling for the geometrical combination rule, for both mixtures, methanol oxygen with water oxygen and propanol methyl group with water oxygen, to improve the density, dielectric constant, excess volume ( V ), and excess enthalpy ( H ) for different molar methanol/water concentrations.
In the case of 1-propanol in the AA scheme, a force field focuses on modifying only the partial charges of the hydroxyl functional group, thereby achieving self-diffusion, which is close to the experimental value [19]. However, the obtained dielectric constant is lower than the experimental value after parameterizing (Table 1). Meanwhile, for UA model, there was a parameterization for 1-propanol force field [20] using the 3SSPP methodology [9]. The obtained results are listed in Table 1.
Recently, Zangi [21] combined the OPLS/AA force field to model hydroxyl groups for alcohols and the optimized OPLS/AA for long hydrocarbons (L-OPLS) for the hydrocarbon tails to overcome the high attraction and prevent liquid-solid phase transition at ambient temperature. Although bulk density and enthalpy of vaporization were calculated, this force field did not improve self-diffusion coefficients; they were over or under predicted values for most of the studied alcohol molecules.
Finally, a new parameterization method, four steps systematic parametrization procedure (4SSPP), was developed by García-Melgarejo et al. [11]. The 4SSPP fitted the solubility, self-diffusion coefficient, surface tension, and dielectric constant for 1-propanol and other compounds in the TraPPE-UA force field. This methodology improves the 3SSPP force field [9] where solubility and self-diffusion coefficient are close to experimental values.
The rest of the work is summarized as follows: the computational details are described in "Computational details", and the 3SSPP/bond methodology to obtain the 1-propanol force fields is described in "Methodology". Then, the results obtained and discussions are given in "Results and discussions". Finally, conclusions are made in "Conclusions".

Computational details
We developed molecular dynamics simulations [22] in GROMACS software version 5 [23]. The isothermalisobaric (NP T ) ensemble at T = 298.15 K and P = 1 atm was used to calculate the density, dielectric constant, and self-diffusion coefficient. The canonical (NV T ) ensemble was used to obtain surface tension. Nosé-Hoover thermostat (τ T = 0.5 ps) and Parrinello-Rahman barostat (τ P = 1 ps) were used to keep temperature and pressure constant, respectively. In addition, a time step of t = 0.002 ps and periodic boundary conditions were considered. A cutoff radius of 1.2 nm and 2.5 nm for the bulk and liquid-vapor was used, respectively. The linear constraint solver (LINCS) algorithm and the particle-mesh Ewald (PME) method were established.
A system with 500 1-propanol molecules in the NP T ensemble for the OPLS/AA [1][2][3] and UA [20] force fields was built to obtain the density, dielectric constant, and selfdiffusion coefficient. LigParGen software [1][2][3] was used to get the force field and molecule configuration in the AA scheme. An NV T system with 2000 molecules in a box of volume, V = L x × L y × L z = 5.3 × 5.3 × 15.000 nm 3 , was created to obtain the surface tension, and a reciprocal grid of 0.12 nm in the xand y-directions and 0.294 nm in z-direction was set in NV T simulations.
Solubility is obtained in the NP N AT ensemble [24], where N is the particle number, P N is the normal pressure (z-direction), A is the interfacial area (x − y plane), and T the temperature. To keep the P N constant, the L z was varing. This liquid-liquid simulations were carried out at constant molecules (4000 water molecules and 1000 alcohols), where temperature and pressure were set at 298 K and 1 bar, respectively, and the surface area was equal to A = L x × L y = 5.1×5.1 nm 2 of dimension. For both, NV T and NP N AT ensembles, the truncation distance was 2.5 nm, and long-range corrections for energy and pressure were not considered.
All simulations in the NP T ensemble were running up to 200 ns, and in the last 5 ns, the density, self-diffusion constant, and dielectric constant were obtained. In addition, surface tension and solubility were calculated with 5 ns after 10 ns of equilibration.

Methodology
A molecular dynamic simulation was executed to obtain the aforementioned target properties, and the obtained results are listed in Table 1. The original OPLS/AA force field (column 3) gave a high percentage error (> 5%) for the dielectric constant and self-diffusion coefficient, and low percentage error (< 5%) for surface tension. In addition, the percentage error in density was < 2% because this force field was obtained by adjusting this property. The OPLS/AA force field was compared with a modified OPLS/AA force field [19], where the selfdiffusion coefficient has a percentage error of 11.8-1.6%, with respect to three experimental data, and the dielectric constant has a percentage error of 28%. This model showed a dielectric constant with a high percentage error with no surface tension value. The parameterized UA 1-propanol is a modified force field that improves the original force field [20] by applying the 3SSPP methodology [9], as observed from the last column of Table 1. This methodology involves the following procedure: 1. All partial charges of the molecule are scaled to obtain the dielectric constant in the NP T ensemble. 2. The energy LJPP, , is modified to obtain the surface tension in the NV T ensemble. 3. In the NP T ensemble, the distance LJPP, σ , is scaled to obtain the density of the system. 4. Finally, we repeat the 1-3 steps, to obtain < 5% deviations from experimental values.
We observed in the 3SSPP methodology that the dielectric constant, density, and surface tension have percentage errors < 5%, but the self-diffusion coefficient has a high percentage error, see [9] and Table 1.
On the contrary, Harris and Yung [17], modified the LJPP and reduced the bond distance by 1%, to improve the CO 2 force field. In this work, we used the same scale to modify all bonds at the same time to be consistent with the 3SSPP methodology where all sigmas, all epsilons, and all charges are modified or scaled in the same way. In the case of 1-propanol bonds, reduction was used, we also used the 3SSPP methodology to formulate a new methodology called 3SSPP/bond which is summarized below: 1. The 3SSPP methodology is applied to a force field. 2. After, bond distance is modified to increase the selfdiffusion. 3. Finally, the 3SSPP methodology is applied again to obtain experimental values with 5% or less percentage error.
The step 1 of the 3SSPP/bond methodology applied to 1-propanol force fields to reproduce experimental surface tension, dielectric constant, and density but self-diffusion, to improving this property, scaling is applied to all bonds.
Step 2, of the proposed methodology, increases the selfdiffusion coefficient, see Table 2; the second column has the goal properties obtained by applying the usual 3SSPP methodology to the 1-propanol original AA force field (step 1 of the new methodology). The reduction of bond distance for the molecule monotonically increases the density and self-diffusion coefficient. On the contrary, the dielectric constant and surface tension reach their maximum value for a bond length of −2%.
Step 2 arises as an empirical procedure where we note that self-diffusion increases when The structure factor (S(Q)) and potential energy surface were calculated to measure the effect of bond reduction of the molecule. The molecular structure was calculated from the X-ray structure factor S(Q) in the reciprocal space Q. Figure 1A shows the original and modified force field spectrum. The modified structure factors were not shifted; however, they were slightly modified from the original ones in parts of the spectrum. The S(Q) is obtained by considering the combinations of the Fourier transform,ĝ i j (Q), for the radial distribution functions, g ij (r), where i and j run for all atom types in the system. In this way, we have where n is the number of atom species and c i are the molar fraction for each atom. The X-ray scattering form factors f i (Q) are defined by where the parameters a i , b i , and c are taken from [32].
On the other side, the energetic effect of the modified bond distance was measured by scanning the potential energy surface using quantum mechanics, which indicates that the potential energy was increasing, as seen in Fig. 1B,   Fig. 1 A Structure factor for 3SSPP AA 1-propanol for 0%, 1%, 2%, and 3% bond distance reduction. B Distance scan for AA 1-propanol molecule: comparing bonding distance and energy between original and parameterized force fields where it formed a parabola-like shape at distances close to the original bond (1.529 nm). This original bond distance was the lowest point of the curve, and the distance where the bond was reduced by −2% had a difference of 1.05 kJ/mol in the potential energy with an error of 0.0002% with respect to the original distance. This reduced distance was taken where the dielectric constant and surface tension have maximum values. Gaussian 09 package was used to obtain the potential energy with the theory b3lyp/6-311++g**.
In the UA force field, the modification of the force field was carried out by taking the parameterized force field [20] and applying a bond distance reduction to the molecule, a similar way to the AA force field; the maximum in dielectric constant and surface tension were obtained corresponding to −2% reduction in the bond distance using this process. Moreover, Table 3 suggests that density and self-diffusion coefficient increased as the bond distance decreased.

Results and discussions
The effect of applying this new methodology in the AA and UA 1-propanol force field models can be deduced from  Table 4, in which the estimated properties by reducing 2% of the bond distance for both force fields are listed.
With the parameterized OPLS/AA force field, the dielectric constant, density, and surface tension exhibited percentage errors < 3% of all experimental data. For the self-diffusion coefficient, the simulational data was close to the lowest experimental value, with < 1% error. Using the same molecular dynamic parameters and a force field for UA, [20], we obtained the results tabulated in Table 4, in which percentage errors < 3% were obtained using 3SSPP/bond UA (third column) for almost all experimental data through the 3SSPP methodology [20]. The self-diffusion coefficient was the only exception, being lower than the experimental data and slightly higher than the modified force field. Hence, the proposed procedure to parameterize the 1propanol force field was applied to AA to reproduce the surface tension, density, dielectric constant, and selfdiffusion coefficient. The original and new force field for AA and UA schemes are shown in the supplementary information.

Structure factor and bond distance
In the AA model, the structure factors for pure substances in the original and the modified force fields were compared with the experimental spectrum obtained by X-ray diffraction [33] (Fig. 2A). The difference between both simulations was imperceptible; the molecular structure was not affected by reducing the bond distance. An insignificant difference between the first maximum of both simulations and the experimental spectrum was the most meaningful change. For the UA model, depicted in Fig. 2B, we noted a similar behavior with respect to the original force field. The structure factor of the parameterized force field had fewer deviations from the experimental spectrum.

Temperature dependence
The temperature dependence of self-diffusion coefficient, density, and dielectric constant is shown in Fig. 3A   Temperature dependence of A self-diffusion coefficient and B density obtained using AA and UA models in the 1-propanol force field. The black line is the experimental data, and the circles and triangles correspond to the data obtained from the UA and AA models, respectively. The red and blue lines correspond to the values obtained from the original and 3SSP/bond force field, respectively were close to the experimental self-diffusion coefficient and their behavior at low temperatures was similar to that of the original one. The diffusion values obtained from AA were closer to the experimental value than those obtained from the UA-modified force field. The UA (3SSPP/bond) obtained self-diffusion coefficients greater than the 3SSPP UA force field (Tables 1 and 4). The values of self-diffusion coefficient obtained of the UA 3SSPP/bond force field are not better than the original UA value, but they are better than the 3SSPP methodology in 24%. Also, for density (Fig. 3B), dielectric constant (Fig. 4), and solubility (Fig. 6) shown that the 3SSPP/bond is better and close to the experimental values for all the temperatures than the original force field. That was the improvement that was achieved, but apparently, the methodology is better applicable to AA than UA models. The error bars for density and self-diffusion coefficient (Fig. 3) are so small, and they were removed from the graphs. The density was obtained by adjusting the force fields to the experimental values, with low errors (around 2% or less). Figure 3B shows that the results obtained from all Fig. 4 Dependence of dielectric constant on temperature obtained using AA and UA model in the 1-propanol force field. The black line reflects the experimental data, the circles represent the UA model, and the triangles correspond to the AA model. The red and blue lines represent the original and 3SSP/bond force field, respectively the force fields are close to the experimental values; the obtained results are more accurate than the original ones for the AA and UA parameterized force fields. The data display a similar constant slope after parameterizing the force fields.
Low dielectric constant and/or super self-diffusion coefficients are the common characteristics across all force fields. In Fig. 4, the dielectric constant in the new force field is greater than that in the original force field, whereas the self-diffusion coefficients were higher for both the AA and UA models, particularly the AA force field.

Solubility and heat capacity
The solubility was used to parameterize force fields and was related with charge distribution [11]. This property was calculated for the original and parameterized force fields and the water models used are SPCE [34] and TIP4P/ [13]; both are the best models that fit basic experimental properties (Supplementary Information). The solubility of a solute i in solvent j was obtained from the density profiles, where M i is the solute molar mass in kg/mol, and ρ i (z) and ρ j (z) are the density profiles for the solute and solvent, respectively. The solubility is calculated by considering the density rate of the solute plateau divided by the solvent, then the result is multiplied by 1000 g/L because 1propanol is miscible in water so the maximum solubility of 1-propanol in water would be 1000 g/L. The initial configuration to calculated solubility was a system with the same volume of 1-propanol and water molecules separated by 0.5 nm between them, then a molecular simulation was implemented. Figure 5 shows the density profiles of parameterized and original AA 1-propanol force fields in water. The mixture of 1-propanol with SPCE is shown in Fig. 5A and for TIP4P/ in Fig. 5B. Figure 5A shows complete miscibility before modifying the force field and it is maintained after applying the 3SSPP/bond methodology. On the other hand, Fig. 5B shows a similar behavior that concludes that AA model is soluble for both water models.
On the contrary, Fig. 6A shows density profiles of parameterized and original UA 1-propanol force fields and SPCE water. The functions oscillated indicating partial miscibility. For the original 1-propanol, the solubility was 194.2 g/L; after parameterizing the force field, this value increased to 322.4 g/L. The solubility was improved by considering the TIP4P/ water model (Fig. 6B) where the value for solubility before applying the 3SSPP/bond methodology was 280 g/L, and it is greater than that with the SPCE water model. However, in the parameterized scheme, the 1-propanol is miscible in TIP4P/ water model.
The solubility for UA was modified after changing the force field, as observed in Fig. 6 that behavior is due to the change in the charge distribution as a result of  Density profiles were obtained from the UA force field in water, for A SPCE and B TIP4P/ . The solid and dashed lines indicated the density profile for water and the 1-propanol, respectively. The red and the blue lines correspond to the data obtained using the original and 3SSPP/bond force fields, respectively the reparameterization. A similar study was developed by García-Melgarejo V et al. [11], who observed that the different charge distribution modified the molecule distribution. In the case of AA, the solubility remains constant and close to experimental value after and before reparameterization.
The criteria to consider one component to be miscible in another are that the density of one component is lower that the other. In the case of Fig. 5 and reparmaterized models of Fig. 6, we considered they miscible because not a density profile greater than, or there are not regions where there is more concentration of one component than other as the original force fields shows in Fig. 6, in this case the miscibility is partial.
However, this methodology has not been tested for the AA scheme. The OPLS/AA 1-propanol force field is watersoluble before applying parameterization, see Fig. 5, so the 4SSPP methodology cannot be applied to this force field.
There is not a study of miscibility in AA, in the case of UA was shown in [11] the 4SSPP methodology to improve force fields that do not reproduce experimental solubility value, they proved different charge distribution to reach this value after that they applied the 3SSPP methodology. In that work, they conclude that solubility is function of charge distribution. However, this methodology has not been tested for the AA scheme. The OPLS/AA 1-propanol force field is water-soluble before applying parameterization, see Fig. 5, so the 4SSPP methodology cannot be applied to this force field.
The heat capacity at constant pressure was calculated for 1-propanol at 298.15 K, we observed that the UA values for both original and reparametrized were closed to the experimental ones; while for AA models, both were to high with respect to the experiment [35]. Table 5 below shows the experimental heat capacity of 144.96 J· K −1 · mol −1 at 299.25 K compared with simulations.
This methodology was being proved to other force fields in the AA and UA molecular models, and we obtained better results with AA than with UA scheme, mainly with the self-diffusion coefficient.

Conclusions
In this work, we studied and proposed a new force field for 1-propanol AA and UA schemes by applying a new methodology called 3SSPP/bond. We noted that the proposed parameterization method provided consistent simulation results with experimental findings for AA with minimal error pertaining to the dielectric constant, density, surface tension, and self-diffusion coefficient. For the UA model, the self-diffusion coefficient was lower than the experimental value but increased with the 3SSPP methodology.
The structure factor for the new force field used in both molecular models did not change significantly with reduction in bond distance and the difference in bond potential energy was not modified significantly; these configurations indicated that this parametrized model was suitable to parameterize the 1-propanol force field.
The temperature dependence of density and dielectric constant fitted the experimental data using the AA and UA models. The temperature dependence of the self-diffusion coefficient exhibited the right behavior for the temperature range obtained for both atomistic 1-propanol models.
In the case of solubility, the AA force field was watersoluble before applying the 3SSPP/bond methodology. This solubility remained constant after parameterizing by using two different water models, TIP4P/ and SPCE. For the UA model, the obtained solubility was low for the original force field; the parameterization increased the solubility and for TIP4P/ water model, this property reached the experimental value.