Method of undetermined parameters

The basic constraint dynamics method employed There are two distinct methods of constraint dynamics the analytical method and the method of undetermined parameters. The reason for the names will become apparent. 

The analytical method deserves a detailed discussion for at least two major reasons. First, if used in conjunction with some constraint correction scheme," it is important in its own right as a practical method of solution of the constrained dynamics problem. Second, the method of undetermined parameters, central to the subject matter of this chapter, is an outgrowth of the analytical method hence a thorough understanding of the analytical method is an essential prerequisite for understanding the method of undetermined parameters. 

As mentioned above, Ryckaert et al. introduced the method of undetermined parameters, which is essentially the analytical method modified to ensure that the constraints are satisfied exactly at each time step. The method of undetermined parameters is discussed in detail, in its most general form, in a later section. Basically, the aforementioned modification requires that the highest derivatives with respect to time of the Lagrangian multipliers (i.e., of order be replaced by a set of undetermined parameters with values to be determined such that the constraints are satisfied exactly at each time step. As will be seen, the algorithm does not introduce into the trajectories additional numerical errors any worse than the error already present in the integration scheme itself. 

Again, the method of undetermined parameters deserves a detailed discussion in its most general form because most researchers are interested in its implementation with different integration algorithms (e.g., basic Verlet,- i velocity Verlet ), with holonomic constraints of various types (e.g., bond-stretch, angle-bend, and torsional constraints), and with particular techniques of solution (e.g., the ma-... 

The numerical procedure used to solve the final equations The analytical method leads to a system of equations linear in the unknowns (i.e., the Lagrangian multipliers and their time derivatives up to order Therefore standard numerical techniques for solving such systems can be employed. The method of undetermined parameters leads to an additional system of equations generally nonlinear in the unknowns (i.e., the derivatives of the Lagrange multipliers of order s ,3x)- The order of nonlinearity depends on the particular... 

This section describes in detail the method of undetermined parameters for any form of holonomic constraints and with an integration algorithm requiring derivatives of the forces up to arbitrary order The treatment is again based mainly on the work of Ryckaert et al., but more detailed mathematical derivations are given. As before, where fundamental equations are reached, corresponding equations in Reference 5 (often with some notational differences) are noted for comparison. 

The most general form of holonomic constraint is nonlinear in the particle positions. Even the simple bond-stretch constraint is nonlinear. Consequently, Eq. [39] is in general a system of / coupled nonlinear equations, to be solved for the / unknowns (7). This nonlinear system of equations must be contrasted with the linear system of equations Eqs. [10] and [11] (which is also in general part of the method of undetermined parameters) used in the analytical method to solve for the Lagrangian multipliers and their derivatives. A solution of Eq. [39] can be achieved in two steps ... 

The (7), obtained by solving Eq. [39], are substituted into Eq. [37] to provide the displacements necessary to satisfy the constraints. Subsequently, the constrained position vectors r(tQ -I- 5f) are obtained from Eq. [38] by adding these constraint corrections to the partially constrained position vectors. The actual derivatives of order of the forces of constraint must be computed (a priori) in the analytical method, whereas in the method of undetermined parameters, the approximate derivatives of order of the forces of constraint can be computed (a posteriori) if desired, by replacing the X (fo)) by the (7). 

The error in the method of undetermined parameters is shown to be of the same order as the error in the analytical method, but with the former approach the constraints are satisfied exactly at every time step. 

Inserting Eq. [42] back into Eq. [37] of the method of undetermined parameters gives ... 

Therefore, the difference between the trajectory computed with the method of undetermined parameters (i.e., using the ["y]) and the trajectory computed with the analytical method (i.e., using the [) (to))) is of 0(8f ). This is the same order of error assumed in the integration algorithm and present in the analytical method before this modification. Hence the modification in the analytical method leading to the method of undetermined parameters does not introduce an error of lower order in 8r than the error already present. The errors in the trajectory computed with the analytical method and the trajectory computed with the method of undetermined parameters are both of 0(8t ), whereas in the latter method the constraints are perfectly satisfied at every time step (compare Eq. [13] with Eq. [36]). 

USING THE METHOD OF UNDETERMINED PARAMETERS WITH THE BASIC VERLET INTEGRATION ALGORITHM... 

Ryckaert et al. incorporated initially the basic Verlet integration algorithm, known also as the Stormer algorithm,into the method of undetermined parameters. In the basic Verlet scheme, the highest time derivative of the coordinates is of second order, and Eq. [37] with = 0 reduces to ... 

For the Verlet scheme with = 0, the are replaced by the (k = 1,. ..,/), and the method of undetermined parameters incorporating the basic Verlet scheme could equally well be termed, and is often referred to in the literature as, the method of undetermined (Lagrangian) multipliers. However, in general where > 0, the Lagrangian multipliers and their derivatives [X °>(tQ),. . ., computed, in addition to the undeter-... 

Inserting Eq. [49] into Eq. [46] shows that the coordinates given by the method of undetermined parameters are accurate to 0(ht ). The error in the coordinates (local error) is of the 0(bt ) present in the basic Verlet scheme. Therefore, consistent with the preceding error analysis, no additional error is introduced in the method of undetermined parameters, and the constraints are exactly satisfied at every time step. If the approximate forces of constraints are desired, they can be computed a posteriori as... 

This section begins by applying the method of undetermined parameters, with the basic Verlet algorithm, to the treatment of bond-stretch constraints. Since they are the most common type of constraint in MD simulations and were... 

Either matrix inversion or SHAKE can be used to solve the set of / linear equations Eq. [130] for the t. As discussed before, solution for the 7 and (t)] by matrix techniques becomes computationally expensive for systems with large numbers of coupled constraints, so the following presentation is confined to the solution by the SHAKE procedure. Andersen discusses why the velocity Verlet algorithm cannot be incorporated into the method of undetermined parameters in as straightforward a manner as can the basic Verlet algorithm. 

---