c42b1

Parameter Fitting Procedure By Victor Anisimov (victor@outerbanks.umaryland.edu) and Alex MacKerell Jr. (alex@outerbanks.umaryland.edu); December 2007 FITPARAM is a parameter fitting procedure that is primarily designed to fit partial atomic charges and atomic polarziabilities based on the Drude oscillator model to interaction energy data and dipole moments, though it may be applied to fitting of other parameters (see below). It supports optimization of multiple parameters for a series of model compounds sharing common parameter sets. Different weights can be assigned to different target data. Optimized parameters can be restrained to their corresponding initial values by using a parabolic penalty function. FITPARAM performs non-linear least square fitting using the Levenberg-Marquardt algorithm. * Introduction / Overview of functionality * Syntax / Syntax of commands * Keywords / Description of keywords * Format / File format * Example / Input examples * Limitations / Known limitations

Top Overview of functionality The primary purpose of FITPARAM is charge derivation using various types of interaction energy as a target data. In this goal it complements the functionality of the FITCHARGE module which derives charges solely from fitting to electrostatic potentials (see fitcharge.doc). However, FITPARAM is not limited to charge derivation and it can be used to optimize any sort of parameters using any sort of relevant target data. Although implemented in CHARMM the FITPARAM module is in fact a stand-alone procedure. The user is not required to load topology and generate a molecule and if loaded this information will be ignored. The command FITPARAM operates with generic parameters without the awareness whether these are point charges, atomic polarizabilities or dihedral force constants. Thus, any parameter may be optimized using FITPARAM based on basic least square fitting. Correspondingly, calculation of the target property and its gradient has to be performed separately by external routines. Preparing the necessary scripts to accomplish this task becomes the user assignment. Despite this obvious inconvenience, such architectural design is necessary to allow for flexibility in the target data used, which otherwise would not be accessible for the parameter optimization. It also makes possible parameter fitting for a variety of molecules to be performed simultaneously. In order to start parameter optimization FITPARAM needs an initial guess for parameters. It reads the parameter values from the input file, performs a cycle of least square optimization and saves the updated parameters in the output file. Then FITPARAM executes the external job, which performs target property and gradient calculations using the current parameter values and returns the control back to FITPARAM. These iterations are continued until convergence criterion is satisfied. Using FITPARAM to optimize charges requires the user to preserve the total charge of the optimized molecule. This can be accomplished by excluding one charge from the fitting and computing the value of the excluded charge in the external procedure as the difference in total charge of the molecule and the sum of optimized charges. In general it is not important which atomic charge to exclude from the fitting. The general suggestion is to exclude that charge which appears least important in reproducing the value of the target data. Dipole moment information can be provided to FITPARAM in two different ways. A dipole may be presented in scalar form as the total dipole moment or in vector form representing three dipole moment components. These two ways of dipole moment specification are handled independently by FITPARAM. When dealing with multiple molecules each molecule may be represented by its own total dipole moment and dipole components. The total dipole moment value may be defined as the gas-phase experimental value, or QM computed value, or manually defined value the user is targeting as the result of the fitting. Providing manually increased (or HF/6-31G* computed) dipole moments is a standard practice when deriving charges for a non-polarizable model. Providing the dipole moment in vector form has the purpose of minimizing the angle between the model and target dipole moment vectors.

Top Syntax of commands [SYNTAX FITPARAM - parameter optimization] FITPARAM { [NITEr int] [TOLErance real] [COUNter int] NEXPeriments int NPARameters int GUESs int EXPData int PARM int [EXPWeight int] [RESTraints int] { [ANTOINE] [ [NDIPoles int] [FVALues int REXternal str [MULTiplicity int] ] } }

Top Description of keywords NITEr - (optional) maximum number of Levenberg-Marquardt (least square) iterations. Default value NITE=100. TOLErance - (optional) relative tolerance on the convergence of the minimized function. Default value TOLE=1.0E-4. COUNter - (optional) is number of iterations under the tolerance it needs to converge. Default value COUN=2. Greater values can be used to assure the convergence. NEXPeriments - (mandatory) number of target data. Care must be taken when handling dipole moments which can be provided in scalar (single value for total dipole) or vector form (three dipole component values in a row). Three components of a vector are considered as one target value for NEXP keyword. Total dipole moment value counts independently as one experimental value. This is also discussed in the examples section. NPARameters - (mandatory) number of parameters to optimize. GUESs - (mandatory) file unit number for initial parameter guess. EXPData - (mandatory) file unit number for target ("experimental") data. PARM - (mandatory) file unit number to save optimized parameter values. EXPWeight - (optional) file unit number for individual weights assigned to each target value. Default value is 1.0 for each target value. RESTraints int - (optional) file unit number for individual restraint factors assigned to each optimized parameter. Default value is 0.0, which means free fitting without restraints. Following keywords are exclusive. ANTOINE - (mandatory if Antoine Constant calculation is expected) switches on Antoine parameter fitting to find analytical dependence in experimental vapor pressure - temperature dependence. This information is useful for calculation of experimental heat of vaporization. or NDIPoles - (mandatory, if dipoles are provided in vector form) number of dipole moments specified in vector form (X Y Z form). Default value is 0. Note, NDIP keyword counts only the number of dipoles provided in vector form (X Y Z components). Correspondingly, NDIP does not count the number of dipoles provided in scalar form (total dipole moment). FVALues - file unit number. This file provides function values and gradients to FITPARAM. The data in this file have to be computed externally using updated parameter values (which are saved in PARM unit by FITPARAM). FITPARAM needs current target function values and gradients to perform next parameter optimization iteration. REXternal - external procedure provided in the form of a string enclosed in double quotes. This string will be executed by CHARMM. This external procedure is in charge of calculating current function values and gradients. MULTiplicity - (optional) file unit number. The data in this file specify th multiplicity of the optimized parameters. Default value of multiplicity for each parameter value is 1. FITPARAM computes the total charge using the information about charge multiplicity. This data makes no influence on the progress of parameter optimization and is implemented for debugging purpose only. For example, three hydrogen atoms in methyl group carrying the same charge value imply a multiplicity of 3. Correct accounting for charge multiplicity will help FITPARAM to provide meaningful information about total charge of the optimized molecule. Non-charge parameters should be assigned multiplicity 0. This is particularly helpful to when FITPARAM optimizes charges simultaneously with non-charge parameters.

Top File format GUESs unit Includes NPAR number of lines. Format: string (up to 20 characters before real number is encountered), real (recognized by decimal point). Example: "A 17.81671". The string portion of the data will be printed by FITPARAM as the parameter name. One may use the string portion to do some useful work for external job, e.g. "set A 17.81671" would turn the parameter file into a functional CHARMM script file. EXPD unit Includes NEXP number of lines. Format: one real value per string. (Exception is Antoine input file which contains two real numbers in free format. First value is temperature; the second one is vapor pressure.) Total dipole moment is defined by one value. Dipole moment components are specified by three numbers in free format. Example: 1.63 ! dipoleTotal MP2 value -1.5155 0.1868 -0.5583 ! dipoleXYZ Comments can be placed after exclamation mark. This information will be ignored by FITPARAM. PARM unit Includes NPAR number of lines. Format: A20,F16.8 This is an output file which will be created by FITPARAM. FVAL unit Includes NEXP + NEXP * NPAR number of lines. Format: The format is free. This file is created by external job and read by FITPARAM. The external job creating this file has to take care of the following requirements. One real number per string is expected for scalar values. Three numbers are provided for a vector value. First NEXP data are the function values which are computed for NEXP target data using the current parameter values. The order of computed values must be the same as in the EXPD file. Next NPAR * NEXP data are partial derivatives computed in the order of experimental data (first running index) and parameter data (second running index). Example: E1 E2 E3 dE1/dp1 dE2/dp1 dE3/dp1 dE1/dp2 dE2/dp2 dE3/dp2 Where E1, E2, E3 are target (experimental) data; p1 and p2 are optimized parameters. Derivatives of vector properties are provided by three real numbers in a row (see example below). EXPW unit Includes NEXP number of lines. Format: one real number in a row. MULT unit Includes NPAR number of lines. Format: one integer number in a row. REST unit Includes NPAR number of lines. Format: one real number in a row.

Top Input Examples Two examples given in this section illustrate basic functionality of FITPARAM. First example covers optimization of Antoine function parameters. This is an example where the function value and gradient computations are implemented inside the FITPARAM module so there is no need to call an external procedure. Therefore this is an exception to the standard use of FITPARAM. This example is discussed here because of its simplicity and because it is included in The Antoine function is a simple parametric analytical function which is used to describe experimental vapor pressure - temperature dependence. It has the following form: lnP = A + [B / (T + C)], where A,B, and C are fitted parameters, P - pressure, T- temperature. Having this function gives the opportunity to compute derivative of pressure over temperature which is necessary to compute the heat of vaporization of the pure liquid. Following is the content of the antoine.inp script: open unit 11 read form name antoine.ini ! initial guess for parameters open unit 12 read form name antoine.exp ! experimental (target data) open unit 13 write form name antoine.prm ! storage for optimized parameters FITPARAM - NITE 50 - ! maximum number of iterations TOLE 0.001 - ! chi^2 convergence threshold COUN 2 - ! number of consecutive successful steps before convergence NEXP 8 - ! number of experimental data NPAR 3 - ! number of parameters to fit (2 or 3 Antoine coefficients) ANTOINE - ! Antoine coefficient fitting GUES 11 - ! input: initial guess for parameters EXPD 12 - ! input: data to fit to PARM 13 ! output: file to store optimized parameters stop The computation starts from opening two input files antoine.ini and antoine.exp, which contain initial guess for parameters and target experimental data, respectively. Following is the content of antoine.ini file: A 17.81671 B 4705.03330 C -60.75000 Three parameters, NEXP=3, will be optimized starting form the above values. The antoine.exp file contains the following data: 393.15 3.649359 398.15 3.877432 403.15 4.076690 408.15 4.264087 413.15 4.461877 418.15 4.651099 423.15 4.825109 428.15 5.018603 Here we have 8 experimental data, NEXP=8, representing temperature - vapor pressure data. Following parameter values are obtained after the execution of the above script: A 17.84653417 B 4706.72855119 C -61.45937115 The second example illustrates how to set up FITPARAM calculation for a charge parameter optimization calculation. The example is based on hydroxyl charge optimization in ethanol targeting interactions with water, which is a standard charge derivation procedure in the additive CHARMM force field. In the example the alkane charges in methyl group are constrained to their standard values. The methyl group is also kept electro-neutral. Two lone pairs (OLP) are assigned to oxygen atom in this example; therefore the oxygen atom is represented by two point charges (qOLP). Note, the central O atom carries zero charge. The charges to be determined are qC, two qHC, two qOLP, and qHO. From these variables one should be excluded. This charge should be assigned manually to maintain a total charge of zero. In the present example we exclude the methylene group hydrogen atom charge: qH = -1/2 * (2*qOLP + qC + qHO). Remaining variables define NPARM=3. The lone pairs are equivalent therefore they contribute as one parameter qOLP. Correspondingly, the initial guess parameter file "parameters.ini" contains the following data: set qOLP -0.23 set qHO 0.36 set qC -0.06 Because the qOLP parameter has the multiplicity of 2 we declare this in the "multiplicity" file: group Here the parameter multiplicity is preceded with the "group" keyword which indicates the program that an electro-neutral group is being treated. If the methyl group charges were also included in the fitting, then the "multiplicity" file would contain the following data: group group Here we optimize the charge on methyl hydrogens (which has multiplicity 3), whereas the charge on methyl carbon has to be computed externally to keep the CH3 group electro-neutral. The charges specified in the "parameters.ini" file are taken from the CHARMM22 force field. To restrain the charges to their initial values we specify "restraints" file using the following restraint factors: 0.1 0.1 0.1 The smaller the number the weaker the restraint. To choose a particular value one needs to do some experimenting. Next we define the target data in "expdata" file: -4.89 ! C-O-H angle bisector -2.51 ! along C-O line -4.80 ! lone-pair position -4.36 ! along O-H line 1.8 ! dipoleTotal MP2 = 1.6259 -1.5155 0.1868 -0.5583 ! dipoleXYZ Presented is ethanol - water QM interaction energies and the ethanol dipole. First four values in the expdata file are interaction energies for the corresponding four water orientations. The fifth value is the ethanol dipole moment set to 1.8 Debye. Note, that the MP2 dipole is 1.63. The dipole moment is purposefully increased to provide an illustrative example that we are free to define any target data we want our model to reproduce. While we are experimenting with the magnitude of the total dipole moment we certainly want to minimize the angle between the target and empirical dipole moment vectors. Therefore the sixth line of data is the dipole moment specified in vector form. Overall, we have six target data, NEXP=6, and one dipole vector, NDIP=1. Weighting of data is performed via the "expweight" file that instructs FITPARAM of the level of priority assigned to the target data: 1. 1. 1. 1. 100. 100000. Each line in this file applies to corresponding target data in the "expdata" file. The value 1.0 means a standard weight. The larger the value we specify the stronger FITPARAM will try to match the target data during the parameter optimization. One needs to make some empirical adjustments in order to define the optimal value for the each weighting factor. The FITPARAM script to run the optimization has the following form: open unit 11 read form name parameters.ini open unit 12 read form name expdata open unit 13 write form name parameters open unit 14 write form name fvalues open unit 15 read form name expweight open unit 16 read form name multiplicity open unit 17 read form name restraints FITPARAM - NITE 100 - ! max number of iterations TOLE 0.0001 - ! convergence tolerance COUN 4 - ! steps to test the convergence NEXP 6 - ! number of target data NDIP 1 - ! number of dipole moments among the target data NPAR 3 - ! number of parameters to optimize GUES 11 - ! initial guess for parameters EXPD 12 - ! target data PARM 13 - ! place to store optimized parameters between loops FVAL 14 - ! function values and derivatives REXT "./run.pl ./parameters ./fvalues" - ! script to compute derivatives EXPW 15 - ! weight of individual target data, optional MULT 16 - ! parameter multiplicity, optional REST 17 ! weight of parameter restraint, optional stop After reading initial guesses for parameters specified in the "parameters.ini" file FITPARAM will write updated parameter values in the "parameters" file (see unit PARM 13). The "parameters" file can be streamed into a CHARMM script to be run externally to compute energies and gradients using the current parameter values. This is done with the help of the REXT keyword (which stands for Run EXTernal process), which tells FITPARAM how to invoke the external process to compute the function values and derivatives. In this particular example we invoke the string "./run.pl ./parameters ./fvalues". This says that computation will be managed by perl script "run.pl" which we assume will invoke individual computations. The perl script takes two arguments. The "parameter" file contains current parameter values computed by FITPARAM. The "fvalues" file instructs "run.pl" script to save the function values and gradients in the "fvalues" file, because FITPARAM will be reading it (see FVAL 14 unit) after "run.pl" script finishes its computation. The perl script or other user selected external routine has to be written individually for each optimization job. Due to the nature of individual computations and because of the multitude of possibilities it is difficult to provide one set of rules on how to do this in the best way. Preparing such script the user should take care about writing the "fvalues" file in the right format (see above). The data in this file are rewritten with each optimization cycle. For the example shown above the "fvalues" file had the following intermediate values: -5.34805 -2.72626 -5.07762 -4.62209 1.86600 -1.72300 0.03500 -0.71400 30.85000 27.06500 30.33000 -15.22500 -14.00000 0.89123 -6.56414 -2.47438 5.55000 9.83500 5.61500 -25.45000 -4.00000 0.55868 -4.51683 -1.57192 2.59500 1.12500 2.67500 -1.37500 -3.00000 -0.09208 -0.23845 0.21026 Relating these data to the data in "expdata" file helps to reinforce the understanding of the format of the "fvalues" file.

Top Known limitations The present implementation does not preserve the total charge. Therefore, the user is responsible to manage this problem in designing the external procedure. Implementing the charge constraint mechanism is in future development plans.