ParAMS: Parameter Optimization for Atomistic and Molecular Simulations
Leonid Komissarov, Robert Rüger, Matti Hellström, Toon Verstraelen
PParAMS: Parameter Optimization forAtomistic and Molecular Simulations
Leonid Komissarov , Robert R¨uger , Matti Hellstr¨om , and ToonVerstraelen ∗ Center for Molecular Modeling (CMM), Ghent University,Technologiepark-Zwijnaarde 46, B-9052, Ghent, Belgium Software for Chemistry & Materials (SCM) B.V., De Boelelaan 1083, 1081HV Amsterdam, The Netherlands
Abstract
This work introduces ParAMS – a versatile Python package that aimsto make parameterization workflows in computational chemistry and physicsmore accessible, transparent and reproducible. We demonstrate how ParAMSfacilitates the parameter optimization for potential energy surface (PES)models, which can otherwise be a tedious specialist task. Because of thepackage’s modular structure, various functionality can be easily combinedto implement a diversity of parameter optimization protocols. For example,the choice of PES model and the parameter optimization algorithm can beselected independently. An illustration of ParAMS’ strengths is provided intwo case studies: i) a density functional-based tight binding (DFTB) repul-sive potential for the inorganic ionic crystal ZnO, and ii) a ReaxFF forcefield for the simulation of organic disulfides. ∗ [email protected] a r X i v : . [ phy s i c s . c h e m - ph ] F e b hroughout the years, the use of predictive computational models has be-come standard practice in many professional fields such as logistics, economics orR&D , with computational chemistry being no exception. Predictive modelsin this field often approximate the potential energy surface (PES) and derivedproperties for a given chemical system, and can be broadly categorized basedon their level of theory: Quantum mechanical (QM) approaches such as wavefunction or density functional theory explicitly model the electronic structure,generally resulting in a high accuracy and broad applicability. However, explicitelectronic treatments come with a high computational price tag, making calcula-tions unfeasible for larger systems. This limitation can be overcome by introduc-ing more empiricism to the model on the electronic ( e.g. tight binding), atomic( e.g. molecular mechanics) or molecular ( e.g. coarse-graining) level. Such empir-ical models attempt to strike the balance between speed and prediction accuracyby approximating the PES description through the introduction of parameters.Prominent examples include approximate functionals in density functional theory(DFT) , density functional tight binding (DFTB) , machine learning (ML)potentials , and force fields (FF). While some of the best empirical models can closely approximate higher-level QM theories at only a fraction of the computational time , their qualitystrongly varies with different sets of parameters . Additionally, many empir-ical models can only deliver accurate results for a comparably limited chemicalspace, either due to their functional form or parameters being specific to certain(combinations of) elements. These limitations can lead to the existence of mul-tiple parameter sets based on the desired application, e.g.
ReaxFF has distinctfamilies of combustion or condensed phase parameterizations. Such lack of gen-2ral parameters gives rise to the research field of parameter fitting , wherethe task is to find an optimal parameter set, given training data constructed fromchemical systems and their properties of interest. Although the fitting processcan be an appealing solution to the above shortcomings, its practical implemen-tation remains hardly accessible to the broader audience and instead is almostexclusively carried out by specialized research groups. In our experience, the ma-jority of researchers, although being interested in individual parameter fitting,are discouraged by the high barrier that comes with it. The main reason for thisbeing a lack of generalization and transparency: Training data often comes in avariety of formats, optimizers expect a different input all together and the formatin which parameters are stored is specific to each method.
The combinationof these oftentimes results in works that can be hardly comprehended and repro-duced by third parties. In an effort to address the above issues, we introduce theParAMS scripting package for Python. We demonstrate how the package can beused to i) generate density functional-based tight binding (DFTB) two-body re-pulsive potentials for an ionic material, and ii) reparameterize a ReaxFF reactiveforce field for organic disulfides. Step-by-step Jupyter notebooks for the twocase studies are provided in the supplementary material, as well as on GitHub .Additional application examples are available in the package’s documentation .The general optimization problem in the scope of ParAMS, including definitionsof jobs, reference data, training set, and loss function is described in the Methodssection. 3 Results
Here, we illustrate how the ParAMS package can be used to train a two-bodySCC-DFTB repulsive potential . For simplicity, we use ZnO, for which severalprevious parametrizations already exist in the literature . We will approximatelyfollow the approach from Ref. 41, in which the authors reused the electronicparameters from the znorg-0-1 DFTB parameter set and reparametrized thetwo-body Zn-O repulsive potential to reference data calculated for the wurtziteand rocksalt polymorphs of ZnO. With the znorg-0-1 parameters, the rocksaltpolymorph is predicted to be more stable than the wurtzite polymorph, but inexperiments and DFT calculations, the opposite is true, which motivates thereparametrization.The entire code needed for the parametrization is given in the supplementarymaterial . It fits the Zn-O pairwise repulsive potential V rep as a tapered doubleexponential function of the form V rep ( r ) = [ A exp( − A r ) + A exp( − A r )] f cut ( r ) (1)where A , A , A , A are the parameters and f cut ( r ) is a tapering function of theform f cut ( r ) = 12 (cid:18) cos( πrr cut ) + 1 (cid:19) (2)with the cutoff distance r cut = 5 .
67 bohr.With ParAMS, it is possible to either directly define the reference values forthe training set (for example, from literature values) or to automatically calcu-4ate them if no reference values have been given. Here, we illustrate the secondapproach, and perform the reference calculations using the periodic DFT codeBAND in the Amsterdam Modeling Suite (AMS).The training set is chosen to be the a and c lattice parameters of wurtziteZnO, the bulk modulus B of wurtzite ZnO, as well as the relative energies of thewurtzite and rocksalt polymorphs of ZnO, ∆ E = E wurtzite − E rocksalt (per ZnOformula unit). We do not need to specify the reference values themselves, as theywill automatically be calculated by ParAMS.The job collection contains two jobs: lattice optimizations of the wurtzite androcksalt polymorphs. From these two jobs, the a , c , B , and ∆ E quantities canbe extracted. For wurtzite, B can be extracted by requesting that the elastictensor be calculated at the end of the lattice optimization.The reference DFT calculations were run with the PBE exchange-correlationfunctional, a triple- ζ (TZP) basis set, and “Good“ numerical quality (dense k-space and integration grids). For the parametrized DFTB engine, a “Good“(dense) k-space grid was also used, since the results of lattice optimizations canbe quite sensitive to the k-space grid.The optimization was done with the Nelder-Mead algorithm from scipy,with a sum-of-squared-errors loss function. The smallest loss function value wasobtained for A = 0 . A = 1 . A = 0 . A = 0 .
40. The resulting repulsivepotential is shown in Figure 1. Typical Zn-O distances in wurtzite (“w”) androcksalt (“rs”) are indicated with gray lines at 3.8 bohr and 4.1 bohr, respectively.With znorg-0-1 (red line), the repulsive potential decays very rapidly between thetypical wurtzite and rocksalt distances. This decrease is not as pronounced with5 (˚A) c (˚A) B (GPa) ∆ E (eV)DFT (this work) 3.30 5.32 126 − − − − − a and c , bulk modulus B ,and energy relative to the rocksalt polymorph ∆ E = E wurtzite − E rocsksalt (perZnO formula unit)either the potential in this work (black line) or znopt (blue line). This affects therelative stability of the wurtzite and rocksalt ZnO polymorphs.Table 1 compares the resulting ZnO properties from the parametrization in thiswork to the DFT reference data, as well as previous DFTB ZnO parametrizations(note: znorg-0-1 was not parametrized to the DFT data in Table 1, and znoptwas parametrized to also reproduce some adsorption energies). The repulsivepotential in this work closely reproduces the wurtzite lattice parameters a and c and the relative energy ∆ E , and provides a good estimate of the bulk modulus,compared to the DFT reference to which it was trained.With ParAMS it is additionally possible to evaluate the loss function, and in-dividual training set entries, using any of the engines in the Amsterdam ModelingSuite. For comparison, Table 1 also gives the corresponding quantities for theUFF (Universal Force Field) engine, which performs significantly worse than theDFTB parametrizations for these quantities.6 a) ZnO wurtzite (b) ZnO rocksalt(c) Twobody Zn-O repulsive potentials V r ep ( H a ) Zn-O distance (bohr)
Figure 1:
ZnO parameterization data. (a-b) pictures of the unit cells of ZnOwurtzite and rocksalt polymorphs. Zn is shown in blue and O in red. (c) DFTBZn-O pairwise repulsive potentials for the znorg-0-1 (red), znopt (blue), andrefitted repulsive potential from this work (black). The gray lines mark typicalZn-O distances in the wurtzite (w) and rocksalt (rs) polymorphs, respectively.7 .2 ReaxFF parameterization ReaxFF is another contemporary example of an empirical model. . This for-malism has been applied to a wide range of chemical problems, and consequentlyhas seen a lot of new parameter development (for a general overview ofReaxFF and its development, see Senftle et al. ). In this section, we demon-strate the parameterization of ReaxFF with a training set previously publishedby M¨uller and Hartke (MH) . The optimized parameter vector x ∗ , as found byMH, is called Mue2016. An overview of the training set can be found in Table S1.It features a total of 231 geometries needed for the computation of 4875 chemicalproperties. Additionally, MH included a validation set to check for overfitting.Examples of three structures included in the data are depicted in Fig. S1, showingcyclopentathione, diphenyl disulfide and dimethyl disulfide. Prior to the param-eter optimization, we evaluate the training and validation sets with the Mue2016ReaxFF parameters and report sum-of-squared-errors (SSE) losses of 14441 and14451 respectively. Note that in the original publication, MH report a trainingset loss of 12400 , while in a more recent work, Shchygol et al. calculate a loss of16300 . Such differences are expected, because numerical instabilities inherentto ReaxFF and software improvements (mostly related to geometry optimization)may result in different optimized geometries. In our setup, we use the covariance matrix adaptation evolution strategy(CMA-ES) as the optimization algorithm with Mue2016 as the initial point.CMA-ES is gradient-free, and relies on a population to sample new parametervectors from an adapted, n-dimensional Gaussian. It does not require additionalhyperparameters other than a population size and an initial width of the Gaus-8ian distribution, σ . Here, we use a population size of 36 and an initial σ of 0.3.Furthermore, we limit the optimization to 24 hours and set up an early stoppingmechanism based on the validation set. The optimization is set up to stop earlyonly if there has been no improvement in the validation set error for the last 6000evaluations.Rather than optimizing the same 87 parameters as MH, we perform a one-dimensional scan on all parameters and select the 35 most sensitive with respectto the training set: Although Mue2016 is a set of 701 parameters in total, onlya subset of these significantly affects the overall cost function value. This is forexample the case when a model includes parameters for each chemical element( e.g. C, H, O), but the total training set of systems R can be constructed fromfewer elements ( e.g. C, H). In such cases, the dimensionality of the problem canbe reduced by scanning for a relevant parameter subset which yields the biggestchange in the cost function value. The simplest setting, which we used in thiscase study, only modifies one parameter at a time to determine its influence onthe objective function. It is also possible to scan all parameter combinations, todiscover coupling between parameters, albeit at a highly increased computationalcost. Out of the 35 parameters selected this way, 16 have also been optimized byMH. We list all optimized parameters in the provided Python notebooks .Parameter bounds are set to be relative to the initial values such that x ± = x ± . | x | . In addition to box constraints, ParAMS enables a definition ofinequality constraints. As the ReaxFF formalism works with bond orders, welimit the parameters responsible for the covalent radii of σ , π and ππ bonds to r σ ≥ r π ≥ r ππ for every atom and atom pair defined in the force field. Thisapproach effectively limits the search space and is available in combination with9ll optimizers.A summary of all settings is provided in Table S2. To compensate for therandomness of CMA-ES, we repeat the optimization set-up nine times. For thebest solution, we report improved training and validation set losses of 11877 and5377 respectively. We make this work’s optimized parameter set available throughthe supplementary information under the title MueParAMS. Correlation plotsbetween reference and predicted values for the new parameters are presented inFigure 2, showing very good agreement to the reference data. Moreover, Figure3 compares the S-S dissociation curve of diphenyl disulfide, as computed withMue2016 and the new MueParAMS parameters, showing an improved agreementto the reference data for this case. With ParAMS, we have presented a modern Python package that allows the cre-ation of versatile parameterization workflows with minimal effort. Its integrationwith the Amsterdam Modeling Suite adds a high amount of flexibility throughthe number of properties that can be fitted alongside the support for multiplecodes when it comes to the model, optimization algorithm and reference data se-lection. Features such as highly customizable optimizations, support for multiplevalidation sets, or the intuitive processing of data aim to make ParAMS accessi-ble to both, advanced and less experienced users. At the same time, developerscan easily extend existing functionality. In Sec. 1.1, we showed how an SCC-DFTB repulsive potential could easily be parameterized for the inorganic crystalZnO. The reference data was calculated automatically using a DFT engine withinAMS. Table 1 also demonstrates how ParAMS can be used to compare the accu-10igure 2:
Training data correlation plots.
Showing (a,b) energy differences,(c,d) atomic forces, (e,f) atomic distances, (g,h) interatomic angles and (i,k) di-hedral angles as calculated with the Mue2016 (left) and MueParAMS (right) pa-rameters. X and Y axes depict reference and predicted values respectively.11igure 3:
S-S dissociation curve of diphenyl disulfide.
Computed withthe Mue2016 and MueParAMS ReaxFF parametrizations. For details about thereference data, see Ref. 48.racy of different chemical simulation packages given the same training set. Sec.1.2 demonstrated how the package can be used to easily process, set up and starta fitting procedure for ReaxFF. Using previously published data by M¨uller andHartke , we were able to find parameters that produce a considerably lower errorfor the validation set while maintaining a similar accuracy in the training data.In future, we hope to extend the number of empirical models that can be fittedwith ParAMS and further improve the ease of use through the introduction ofadditional shortcut functions for training set building. We also expect additionsin other functionality such as optimization algorithms or extractors based on userfeedback and wishes as the project matures. The package is included in all AMSreleases since 2020. 12 Methods
This section provides a mathematical framework for the parameterization prob-lem. We assume that the training data can be defined as a set of physico-chemicalproperties for a number of isolated or periodic systems. Examples for relevantproperties are energy differences, nuclear gradients or system geometries. In thecontext of ParAMS, we define an arbitrary property P that can be expressed asthe output of a computational job. When working with multiple jobs as part ofa training set, a job function can be defined as J ( R j , S j , M ) = (cid:0) R j , P nj (cid:1) ∀ j ∈ { . . . N job } , n ∈ { . . . N prop ( j ) } , (3)calculating the output geometry R j and all properties P nj of a job j . Theinput for every job in J consists of the input geometry R j , the job settings S j ( e.g. geometry optimization and frequencies) and the computational model M .Note that a parametric model is additionally a function of the parameter vector x , in which case the outputs of the above equation can be denoted with the hatoperator ( i.e. ˆ R j , ˆ P nj ), as to distinguish between reference properties and proper-ties predicted by the parametric model. Training set entries can be constructed,for example, from a linear combination of multiple properties y i = N lc ( i ) X k =1 c i,k P n ( i,k ) j ( i,k ) ∀ i ∈ { . . . N data } , (4)where c i,k is the coefficient for term k of training set entry i and N lc ( i ) is the totalnumber of terms per entry. Non-linear combinations of properties to construct13 i are also possible. Such a formulation offers a high degree of flexibility for theconstruction of a training set. One example is the combination of multiple systemenergies into one reaction energy. It should be noted that a training set entry,as defined in Eq. 4, does not have to originate from the results of computationaljobs. The reference value can instead be provided directly, making it easy to workwith experimental or external data.While training set entries y have to be defined only once, their predictedcounterpart ˆ y has to be re-calculated every time the model parameters change.For this purpose, we introduce a Data Set function DS ( x | y ) = y − ˆ y , (5)which extracts all properties needed for the calculation of ˆ y based on a parameterset x and returns the respective vector of residuals. A metric in the form of a lossfunction L ( y − ˆ y , w ) is then applied to the residuals for a qualitative measure ofhow close reference and predicted values are. The additional weights vector vector w can be used to balance possibly different orders of magnitude in the data setor make certain entries more relevant for the fitting process than others.Finally, the optimization algorithm can be defined as a function that minimizes L with respect to the parameters O ( x , L ) = arg min x L = x ∗ , (6)finding an optimal solution x ∗ from an initial point x .14 .2 Implementation ParAMS follows a modular package structure with well-defined application pro-gramming interfaces (APIs), which allows components to be treated indepen-dently. This is essential for future development, as individual sub-modules canbe easily worked on and extended. We describe the main components and theirfunctionality below.
Job Collection and Data Set classes are responsible for the input / output(IO) of relevant data. In the context of ParAMS, these will be collections ofchemical systems R , properties P and settings S alongside optional metadata.A Job Collection clearly defines job entries by combining systems and settings.The Data Set defines which properties of a job are relevant to the optimizationand stores the reference values of all entries in y . Reference values can be addedfrom any source: Experimental / external results, or a high-level calculation. Toensure reproducibility and ease-of-use, ParAMS makes use of the YAML data-serialization format as the default for all IO operations. Extractors tell ParAMS how to extract a property of interest P from a cal-culated job. Technically, they are small standalone Python modules that readthe native ( e.g. stream, text or binary file) output of the Amsterdam Model-ing Suite into Python variables. Examples for implemented extractors are:Interatomic distances, valence angles, dihedral angles, atomic charges, reactionenergies, linear transit energy profiles, lattice vectors, bulk moduli, atomic forces,stress tensors, Hessian matrices and vibrational frequencies. New extractors areeasily written by the user, effectively allowing any property that can be calculatedwith AMS to be fitted within the scope of ParAMS. Extractors also support more15laborate cases that require additional processing before a comparison. This isfor example needed when computing the minimal root-mean-square deviation ofatomic positions. Loss Functions implement various metrics that describe the distance be-tween two vectors. In the context of parameter fitting, a vector consisting of allreference values y has to be compared to the predictions vector ˆ y , as generatedgiven a specific set of parameters, in order to measure the quality of the fit. Im-plemented metrics are least absolute error (LAE), mean absolute error (MAE),root-mean-square error (RMSE) and residual sum of squares (RSS). Additionally,user-defined metrics are supported. Optimizers provide a unified interface to a variety of optimization algorithms.Currently, the following are supported: Covariance Matrix Adaptation EvolutionStrategy (CMA-ES) , Adaptive Rate Monte Carlo (ARMC) and optimizersavailable through the Nevergrad and SciPy packages. To guarantee parame-terization support for a wider range of models, the current version of ParAMS isdesigned to work with gradient-free optimization algorithms only. Parameter interfaces translate the parameter vector x into the native for-mat of the empirical model ( e.g. a file on disk). Any existing parameter interfacecan be parameterized. At the time of writing, ParAMS supports interfaces toReaxFF , SCC-DFTB repulsive potentials, GFN1-xTB , and Lennard-Jonespotentials. Callbacks allow the interaction with a parameterization at runtime. Suchinteractions can be progress loggers, timeouts, early stopping criteria or plottingfunctions . 16igure 4: A schematic representation of the ParAMS parameterizationloop.
Highlighting the interplay between the three main components: Optimizer,Job Collection and Data Set. Every parameter set suggested by the optimizerproduces different job results, which are evaluated by the Data Set based on aunique jobID . The results are then compared with the reference values, weightedand combined into the loss function value. A more detailed description is availablein the package documentation.The main script is a command line interface to ParAMS for users who donot wish to spend much time writing their own parameterization scripts. It allowsthe setup and execution of the most common tasks through a configuration file.Figure 4 shows the general parameterization loop and highlights the maininput-output relationships. The implementation allows users already experiencedin other data science packages to prepare inputs and process results in a familiarway. Techniques like training- and validation set splitting, cross-validation, earlystopping, batching or outlier detection are supported out of the box. Parameterconstraints can be further used to limit the search space of a problem. In addi-17ion to regular box constraints, users can express inequality constraints involvingmultiple parameters ( e.g. p ≤ p ), as was demonstrated in Section 1.2. ParAMSimplements two levels of parallelism. Multiple parameter vectors and multiplejobs per vector can be evaluated at the same time, resulting in workloads thatcan be distributed effectively based on the training set size and optimization algo-rithm. Moreover, since the signatures of all classes in a submodule are the same,different models and optimizers can be effortlessly compared and deployed ( e.g. comparison of different levels of theory in Tab. 1). Daily regression tests areperformed to guarantee an error-free functionality of the package. Acknowledgements
We thank Micha l Handzlik for the initial design of the ParAMS package and Dr.Tom´aˇs Trnka for the implementation of the AMSWorker interface, resulting ina considerable speed-up of computation times. This project has received fund-ing from the European Union’s Horizon 2020 research and innovation programmeunder grant agreement No 814143 (L.K. & T.V.), No 798129 (M.H). T.V. alsoacknowledges funding of the research board of Ghent University. R.R. and M.H.have received funding from the Netherlands Enterprise Agency (RVO) and Stimu-lus under the MIT R&D Collaboration programme, project number PROJ-02612.The computational resources (Stevin Supercomputer Infrastructure) and servicesused in this work were provided by the VSC (Flemish Supercomputer Center),funded by Ghent University, FWO and the Flemish Government – departmentEWI. 18 uthor Contributions
L.K. and R.R. developed the ParAMS Python package. L.K. and M.H. performedthe case studies. L.K., M.H. and T.V. wrote the paper. T.V. oversaw the project.All authors read and approved the final manuscript.
Competing Interests
Authors L.K., R.R. and M.H. were employed by the company Software for Chem-istry and Materials (SCM). SCM develops and commercializes the AmsterdamModeling Suite, of which ParAMS is a new module.
Data Availability
All data is available alongside the supporting information at github.com/oiao/params si.The package’s documentation is available at scm.com/doc.trunk/params.
Code Availability
The ParAMS package is available upon free registration for a Amsterdam Mod-eling Suite trial at scm.com.
References
1. Bruzda, J. Multistep quantile forecasts for supply chain and logistics op-erations: bootstrapping, the GARCH model and quantile regression basedapproaches.
Cent. Eur. J. Oper. Res. , 309–336 (2018).19. Wolfers, J. & Zitzewitz, E. Prediction markets. J. Econ. Perspect. , 107–126(2004).3. Lewis, B. P., hung Shih, I., Jones-Rhoades, M. W., Bartel, D. P. & Burge,C. B. Prediction of mammalian MicroRNA targets. Cell , 787–798 (2003).4. Wilson, R. L. & Sharda, R. Bankruptcy prediction using neural networks.
Decis. Support Syst. , 545–557 (1994).5. Geller, R. J. Earthquake prediction: a critical review. Geophys. J. Int. ,425–450 (1997).6. Becke, A. D. A new mixing of hartree–fock and local density-functional the-ories.
J. Chem. Phys. , 1372–1377 (1993).7. Lee, C., Yang, W. & Parr, R. G. Development of the colle-salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B ,785–789 (1988).8. Stephens, P. J., Devlin, F. J., Chabalowski, C. F. & Frisch, M. J. Ab ini-tio calculation of vibrational absorption and circular dichroism spectra usingdensity functional force fields. J. Phys. Chem. , 11623–11627 (1994).9. Chai, J.-D. & Head-Gordon, M. Long-range corrected hybrid density func-tionals with damped atom–atom dispersion corrections. Phys. Chem. Chem.Phys. , 6615 (2008).10. Elstner, M. et al. Self-consistent-charge density-functional tight-bindingmethod for simulations of complex materials properties.
Phys. Rev. B ,7260–7268 (1998). 201. Yang, Yu, H., York, D., Cui, Q. & Elstner, M. Extension of the self-consistent-charge density-functional tight-binding method: Third-order expansion of thedensity functional theory total energy and introduction of a modified effectivecoulomb interaction. J. Phys. Chem. A , 10861–10873 (2007).12. Gaus, M., Cui, Q. & Elstner, M. Dftb3: Extension of the self-consistent-charge density-functional tight-binding method (scc-dftb).
J. Chem. TheoryComput. , 931–948 (2011).13. Grimme, S., Bannwarth, C. & Shushkov, P. A robust and accurate tight-binding quantum chemical method for structures, vibrational frequencies, andnoncovalent interactions of large molecular systems parametrized for all spd-block elements (z = 1–86). J. Chem. Theory Comput. , 1989–2009 (2017).14. Bannwarth, C., Ehlert, S. & Grimme, S. Gfn2-xtb—an accurate and broadlyparametrized self-consistent tight-binding quantum chemical method withmultipole electrostatics and density-dependent dispersion contributions. J.Chem. Theory Comput. , 1652–1671 (2019).15. Behler, J. & Parrinello, M. Generalized neural-network representation of high-dimensional potential-energy surfaces. Phys. Rev. Lett. , 146401 (2007).16. Sch¨utt, K. T., Sauceda, H. E., Kindermans, P.-J., Tkatchenko, A. & M¨uller,K.-R. Schnet – a deep learning architecture for molecules and materials. J.Chem. Phys. , 241722 (2018).17. Smith, J. S., Nebgen, B., Lubbers, N., Isayev, O. & Roitberg, A. E. Less ismore: Sampling chemical space with active learning.
J. Chem. Phys. ,241733 (2018). 218. Shao, Y., Hellstr¨om, M., Mitev, P. D., Knijff, L. & Zhang, C. PiNN: A pythonlibrary for building atomic neural networks of molecules and materials.
J.Chem. Inf. Model. , 1184–1193 (2020).19. Cornell, W. D. et al. A second generation force field for the simulation ofproteins, nucleic acids, and organic molecules.
J. Am. Chem. Soc. , 5179–5197 (1995).20. MacKerell, A. D. et al.
All-atom empirical potential for molecular modelingand dynamics studies of proteins.
J. Phys. Chem. B , 3586–3616 (1998).21. van Duin, A. C. T., Dasgupta, S., Lorant, F. & Goddard, W. A. Reaxff:A reactive force field for hydrocarbons.
J. Phys. Chem. A , 9396–9409(2001).22. Chenoweth, K., van Duin, A. C. T. & Goddard, W. A. Reaxff reactive forcefield for molecular dynamics simulations of hydrocarbon oxidation.
J. Phys.Chem. A , 1040–1053 (2008).23. Smith, J. S. et al.
Approaching coupled cluster accuracy witha general-purpose neural network potential through transfer learn-ing. Preprint at https://chemrxiv.org/articles/Outsmarting_Quantum_Chemistry_Through_Transfer_Learning/6744440 (2019).24. Bannwarth, C., Ehlert, S. & Grimme, S. GFN2-xTB—an accurate andbroadly parametrized self-consistent tight-binding quantum chemical methodwith multipole electrostatics and density-dependent dispersion contributions.
J. Chem. Theory Comput. , 1652–1671 (2019).225. Shchygol, G., Yakovlev, A., Trnka, T., van Duin, A. C. T. & Verstraelen, T.Reaxff parameter optimization with monte-carlo and evolutionary algorithms:Guidelines and insights. J. Chem. Theory Comput. , 6799–6812 (2019).26. Dieterich, J. M. & Hartke, B. Ogolem: Global cluster structure optimisationfor arbitrary mixtures of flexible molecules. a multiscaling, object-orientedapproach. Mol. Phys. , 279–291 (2010).27. Wang, L.-P., Chen, J. & Voorhis, T. V. Systematic parametrization of polar-izable force fields from quantum chemistry data.
J. Chem. Theory Comput. , 452–460 (2012).28. Brommer, P. et al. Classical interaction potentials for diverse materials fromabinitiodata: a review ofpotfit.
Model. Simul. Mater. Sc. , 074002 (2015).29. Jaramillo-Botero, A., Naserifar, S. & Goddard, W. A. General multiobjectiveforce field optimization framework, with application to reactive force fieldsfor silicon carbide. J. Chem. Theory Comput. , 1426–1439 (2014).30. Senftle, T. P. et al. The reaxff reactive force-field: development, applicationsand future directions. npj Comput. Mat. , 15011 (2016).31. Guvench, O., Hatcher, E., Venable, R. M., Pastor, R. W. & MacKerell, A. D.Charmm additive all-atom force field for glycosidic linkages between hexopy-ranoses. J. Chem. Theory Comput. , 2353–2370 (2009).32. Gaus, M., Chou, C.-P., Witek, H. & Elstner, M. Automatized parametrizationof scc-dftb repulsive potentials: Application to hydrocarbons. J. Phys. Chem.A , 11866–11881 (2009). 233. van Beest, B. W. H., Kramer, G. J. & van Santen, R. A. Force fields for silicasand aluminophosphates based on ab initio calculations.
Phys. Rev. Lett. ,1955–1958 (1990).34. Ashraf, C. & van Duin, A. C. Extension of the reaxff combustion force fieldtoward syngas combustion and initial oxidation kinetics. J. Phys. Chem. A , 1051–1068 (2017).35. van Duin, A. C. T. ReaxFF User Manual. (2002).36. Wang, L.-P. Forcebalance: Main page. http://leeping.github.io/forcebalance/doc/html/index.html .37. Dieterich, J. M. & Hartke, B. Ogolem.org homepage. .38. Martinez, J. A. et al.
Potential optimization software for materials (POSMat).
Comput. Phys. Commun. , 201–211 (2016).39. Komissarov, L. & R¨uger, R. Params documentation. .40. Komissarov, L., R¨uger, R., Hellstr¨om, M. & Verstraelen, T. Params support-ing information. https://doi.org/10.5281/zenodo.4543520 (2020).41. Hellstr¨om, M. et al.
An scc-dftb repulsive potential for various zno polymorphsand the zno–water system.
J. Phys. Chem. C , 17004–17015 (2013). https://doi.org/10.1021/jp404095x .242. Moreira, N. H., Dolgonos, G., Aradi, B., da Rosa, A. L. & Frauenheim,T. Toward an accurate density-functional tight-binding description of zinc-containing compounds.
J. Chem. Theory Comput. , 605–614 (2009).43. te Velde, G. & Baerends, E. J. Precise density-functional method for periodicstructures. Phys. Rev. B , 7888–7903 (1991).44. R¨uger et al. Amsterdam modeling suite. https://scm.com (2019).45. Nelder, J. A. & Mead, R. A simplex method for function minimization.
Comput. J. , 308–313 (1965).46. Strachan, A., van Duin, A. C. T., Chakraborty, D., Dasgupta, S. & Goddard,W. A. Shock waves in high-energy materials: The initial chemical events innitramine rdx. Phys. Rev. Lett. , 098301 (2003).47. Fogarty, J. C., Aktulga, H. M., Grama, A. Y., van Duin, A. C. T. & Pandit,S. A. A reactive molecular dynamics simulation of the silica-water interface. J. Chem. Phys. , 174704 (2010).48. M¨uller, J. & Hartke, B. reaxff reactive force field for disulfide mechanochem-istry, fitted to multireference ab initio data.
J. Chem. Theory Comput. ,3913–3925 (2016).49. Furman, D. & Wales, D. J. Transforming the accuracy and numerical stabilityof reaxff reactive force fields. J. Phys. Chem. Lett. , 7215–7223 (2019).50. Hansen, N. & Ostermeier, A. Completely derandomized self-adaptation inevolution strategies. Evol. Comput. , 159–195 (2001).251. Hansen, N. & Kern, S. Evaluating the CMA evolution strategy on multimodaltest functions. In Yao, X. et al. (eds.) Parallel Problem Solving from NaturePPSN VIII , vol. 3242 of
LNCS , 282–291 (Springer, 2004).52. Hansen, N., Akimoto, Y. & Baudis, P. CMA-ES/pycma on Github. https://doi.org/10.5281/zenodo.2559634 (2019).53. d¨ot Net, I., Evans, C. & Ben-Kiki, O. The official yaml web site. https://yaml.org/ (2019).54. van Duin et al.
Reaxff 2019.4. (2019).55. Cosseddu, S. & Infante, I. Force field parametrization of colloidal CdSenanocrystals using an adaptive rate monte carlo optimization algorithm.
J.Chem. Theory Comput. , 297–308 (2016).56. Rapin, J. & Teytaud, O. Nevergrad - A gradient-free optimization platform. https://GitHub.com/FacebookResearch/Nevergrad (2018).57. Virtanen, P. et al. SciPy 1.0: Fundamental Algorithms for Scientific Com-puting in Python.
Nat. Methods , 261–272 (2020).58. Kamat, A. M., van Duin, A. C. T. & Yakovlev, A. Molecular dynamicssimulations of laser-induced incandescence of soot using an extended reaxffreactive force field. J. Phys. Chem. A , 12561–12572 (2010).59. Hunter, J. D. Matplotlib: A 2d graphics environment.
Comput. Sci. Eng. ,90–95 (2007). 26 upporting Information for ‘ParAMS:Parameter Optimization for Atomistic andMolecular Simulations’ Leonid Komissarov , Robert R¨uger , Matti Hellstr¨om , and ToonVerstraelen ∗ Center for Molecular Modeling (CMM), Ghent University,Technologiepark-Zwijnaarde 46, B-9052, Ghent, Belgium Software for Chemistry & Materials (SCM) B.V., De Boelelaan 1083, 1081HV Amsterdam, The Netherlands ∗ [email protected] S1 a r X i v : . [ phy s i c s . c h e m - ph ] F e b able S1: Composition of the reference data published by M¨uller and Hartke ,split by the computational tasks Single Point (SP) and Geometry Optimization(GO). For each of the two sets, the upper part describes the chemical systems,while the lower breaks down the individual entries in the training and validationsets. Note that some entries might be a function of multiple chemical systems,meaning that the sum of SP+GO is not necessarily equal to the total number ofentries for that row ( cf. Sec. 3.1 in the main text).
Training Set SP GO Total
Number of systems 222 9 231Mean system size (atoms) 6.6 11.4 6.8Std. dev. (atoms) 2.9 7.7 3.3Total number of entries 4620 317 4875Energies 219 62 219Forces 4401 0 4401Atomic distances 0 94 94Angles 0 85 85Dihedrals 0 76 76
Validation Set
Number of systems 200 24 224Mean system size (atoms) 24.0 12.7 22.8Std. dev. (atoms) 0.0 5.9 4.0Total number of entries 199 771 970Energies 199 0 199Forces 0Atomic distances 0 281 281Angles 0 257 257Dihedrals 0 233 233S2able S2: Summary of relevant ParAMS settings used for the re-parameterizationof Mue2016.
Setting Value
Number of optimizations 9Number of parameters to optimize 35Lower / upper parameter bounds x ± . | x | Optimization timeout 24 hoursCMA-ES population size 36CMA-ES sigma 0.3Loss function sum of squared errorsEarly stopping patience 6000 evaluationsConstraints r σ ≥ r π and r π ≥ r ππ Figure S1: From top to bottom: Example structures of cyclopentathione, diphenyldisulfide and dimethyl disulfide, containing S (yellow), C (black), and H (white),included in the data provided by MH. The fitted properties include bond distances,angles, relative energies and atomic forces.S3 eferences
1. M¨uller, J. & Hartke, B. reaxff reactive force field for disulfide mechanochem-istry, fitted to multireference ab initio data.
J. Chem. Theory Comput.12