Abstract:
A method for dynamics simulation involves maintaining quantities according to a floating point binary format quantized to a first precision lower than the precision supported by the floating point format. For example, although an IEEE floating point number can represent numbers with a precision of one part in 2ˆ24, the quantities are quantized to a lower precision, such as one part in 2ˆ22. Operations are applied to sets of the quantities by quantizing the intermediate results of the operations to the lower precision than the precision supported by the floating point format.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS  
       [0001]     This application claims the benefit of U.S. Provisional Application No. 60/827,547, filed Sep. 29, 2006, titled “NUMERICAL TECHNIQUES FOR DYNAMICS SIMULATION REVERSIBILITY,” which is incorporated herein by reference.  
         [0002]     This application is related to U.S. Provisional Application No. 60/827,539, titled “CONSTRAINT STABILIZATION,” filed Sep. 29, 2006, and to U.S. application Ser. No. ______ (attorney docket 30023-007001) titled “CONSTRAINT STABILIZATION,” filed concurrently with the present application. These applications are incorporated herein by reference.  
         [0003]     This application is also related to PCT Application No. PCT/US2006/032498, titled “ARCHITECTURES FOR COMPUTATION OF PARTICLE INTERACTIONS,” filed on Aug. 18, 2006, which is incorporated herein by reference. 
     
    
     BACKGROUND  
       [0004]     This invention relates to numerical techniques for dynamics simulation.  
         [0005]     In principle, a dynamics simulation, such as molecular dynamics (MD) simulation can be time reversible because the classical equations of motion are time reversible. For example, if one runs an MD simulation for a thousand time steps, negates the instantaneous velocities for all particles, and runs for another thousand time steps, one should exactly recover the initial particle positions. However, many simulation software implementations fail to achieve reversibility for two reasons. First, roundoff error during integration leads to a loss of state information; particle positions and velocities at one time step cannot be exactly reconstructed from those at the next time step. Second, lack of associativity of floating point summation leads to a situation where computed forces can depend on various factors that affect order of summation, such that these forces are not uniquely defined by particle positions and force field parameters in an MD simulation.  
         [0006]     Large floating point code, scalar or parallel, that ensures exact reversibility for non-trivial dynamics simulation runs is not known to the inventors.  
       SUMMARY  
       [0007]     In one aspect, in general, a method for dynamics simulation involves maintaining quantities according to a floating point binary format quantized to a first precision lower than the precision supported by the floating point format. An operation is applied to a set of the quantities by quantizing the intermediate results of the operation to the lower precision than the precision supported by the floating point format.  
         [0008]     In some examples, although an IEEE floating point number can represent numbers with a precision of one part in 2ˆ24, the quantities are quantized to a lower precision, such as one part in 2ˆ22.  
         [0009]     Aspects may include one or more of the following features.  
         [0010]     A precision is selected to be lower than the precision supported by the floating point format.  
         [0011]     Selecting the precision includes dynamically varying that precision during the course of a dynamic:; simulation.  
         [0012]     Varying the precision includes varying that precision according to a characteristic of a system being simulated in the dynamics simulation. For example, the characteristics may include a temperature at which the simulation is being performed, or a dimension of the simulation volume. For example, at a high temperature or at a large dimension, a relatively greater precision may be selected.  
         [0013]     Selecting the precision includes selecting different precisions for different of the quantities.  
         [0014]     At least some of the quantities are maintained such that during the dynamics simulation the values can exceed a determined range and are then represented at a lower precision than when quantized. For example, if the values do not exceed the determined range the simulation remains reversible, but when the values exceed the range, reversibility may be lost but the accuracy of the simulation is largely maintained thereby providing a graceful degradation of the characteristics of the simulation.  
         [0015]     The maintained quantities include quantities that represent physical features of bodies in the dynamics simulation.  
         [0016]     The maintained quantities comprise integer quantities.  
         [0017]     The maintained quantities comprise quantities with fractional components.  
         [0018]     The maintained quantities comprise both integer quantities and quantities with fractional components.  
         [0019]     Quantizing the intermediate results comprises performing multiple addition and/or subtraction floating point operations.  
         [0020]     Quantizing an intermediate result comprises sequentially adding and subtracting a floating point constant number.  
         [0021]     Applying the operation comprises computing a quantity that exceeds a range in which quantities are exactly represented  
         [0022]     Supporting bitwise exact time reversed simulation using the maintained quantities.  
         [0023]     Applying an operation to quantities maintained at the first precision and to quantities maintained at the second precision, including quantizing intermediate results of the operation.  
         [0024]     Fixed point numbers are stored in a floating point representation transparently. For example, a fixed point number can be treated as a seemingly normal floating point number without need for special treatment or resealing.  
         [0025]     Fixed point number stored in a floating point representation are used as a hybrid fixed point/floating point number. For example, using fixed point arithmetic in the non-exceptional case while keeping the ability to accommodate exceptional out-of-fixed-point-range results gracefully and without loss of performance.  
         [0026]     “Magic number” shifts perform fixed point rounding efficiently.  
         [0027]     The above techniques are used to make a simulation time-reversible.  
         [0028]     The above techniques are used to reconfigure the effective precision used in a simulation without loss of performance. Note that these techniques allow the precision of various quantities can be set independently (e.g., positions have 40-bit resolution and velocities are 3 1-bit resolution).  
         [0029]     The numerical techniques described in the context of dynamics simulation can be applied to other computations, for example, where exact, repeatable, or reversible calculations are required.  
         [0030]     In another aspect, a system includes computational modules that implement the steps of any of the methods set forth above.  
         [0031]     In another aspect, software, which may be stored on a computer readable medium, includes instruction: for causing a data processing system to perform the steps of any of the methods set forth above.  
         [0032]     Aspects can have one or more of the following advantages.  
         [0033]     Reversibility can help ensure that an implementation will accurately model thermodynamic relationships that depend on detailed balance and can result in very low energy drift.  
         [0034]     The numerical techniques described that can be used to provide bitwise reversibility are also applicable in situations in which bitwise reversibility is not strictly achievable for other reasons. MD simulations that employ constraints and that perform pair list updates infrequently or other such techniques may not be able to be rendered reversible by these techniques. Nevertheless, these techniques may be beneficial in those cases, resulting, for (example, in significantly lower energy drift. As another example, simulations using these numerical techniques can be made deterministic in a parallel implementation and/or not dependent on the order or other arrangement of the calculations involved.  
         [0035]     As a technique for controlling energy drift, the numerical techniques described herein can be combined with constrain stabilization techniques described in the copending utility application “CONSTRAINT STABILIZATION.” For example, the procedures outlined in the copending application can be implemented using the numerical techniques described herein to control precision of the computations.  
         [0036]     Other features and advantages of the invention are apparent from the following description, and from, the claims. 
     
    
     DESCRIPTION  
       [0037]     A new numerical technique for dynamics simulation can be used to preserve exact bitwise reversibility for simulations that do not use constraints and in which pairlist updating is performed at every time step. This approach is applicable to at least some MD integration schemes. In one embodiment, loss of information is avoided during integration by performing updates to particle coordinates using fixed point arithmetic. Problems due to non-associativity are avoided by maintaining a consistent ordering of the particles and computations such that computed forces are unique functions of the particle positions and force field parameters.  
         [0038]     Consider a simple unconstrained NVE (constant particle number, constant volume, constant energy) simulation integrated with the velocity Verlet method (in general, in equations below the notation a_b denotes a subscript b and the notation aˆb denotes a superscript b). Let:  
         [0039]     F(X) be a function that produces the vector of forces given the vector of positions.  
         [0040]     X — 0 and V — 0 be vectors that give all the particle positions and velocities at the beginning of a timestep.  
         [0041]     DT be the timestep.  
         [0042]     M be the particle mass matrix (and in the mass of particle i).  
         [0043]     X — 1 and V — 1 be the positions and velocities at the end of a timestep.  
         [0044]     Vectors like V_h anti so forth be computational intermediates.  
         [0045]     The NVE velocity verlet algorithm to advance (X — 0,V — 0) to (X — 1,V — 1) can be expressed as: 
 
 V   —   h=V   — 0 +Mˆ− 1 F ( X   — 0) DT/ 2 . . . half advance velocities 
 
 X   — 1 =X   — 0 +V   —   hDT . . .  advance positions 
 
 V   — 1 =V   —   h+Mˆ− 1 F ( X   — 1) DT/ 2 . . . half advance velocities 
 
         [0046]     Suppose you nun this algorithm starting from (X — 1,−V — 1). It can be proved that algorithm advances this starting point to (X — 0,−V — 0) if done in exact arithmetic. That is, the above algorithm is time reversible in exact arithmetic. It is straightforward to show that the above is still reversible if the following properties hold for inexact arithmetic: 
 
( X+ ( VDT ))+((− V ) DT )= X  
 
(−( V+ ( Mˆ− 1 F ( X ) DT/ 2)))+( Mˆ− 1 F ( X ) DT/ 2)=− V  
 
         [0047]     Above, all operations are done in the order implied by the parenthesis with the inexact operations. No assumptions have been made about the behavior of the inexact operations.  
         [0048]     For standard IEEE floating point arithmetic, the above properties do not hold due to floating point non-associativity. That is (a+b)+c does not necessarily equal a+(b+c) in IEEE floating point, dramatically in situations in which b and c are much smaller than a.  
         [0049]     Suppose the individual components of X and V can be bounded to be less than XMAX and VMAX respectively. This is relatively straightforward for most MD simulations. For example, periodic simulations already have a natural bound for particle positions (e.g., all particles reside within or near the simulation unit cell). Likewise, velocities can be strictly bounded in a constant energy simulation (e.g. for systems with a minimal potential energy configuration, a strict but overly conservative bound can be obtained by concentrating all the kinetic energy and excess potential energy into a single velocity component of the lightest particle). Much tighter, though often probabilistic, bounds can be found. quickly using thermodynamic arguments (e.g. VMAX is such that one should expect to wait several billion years before seeing a particle with such a large velocity).  
         [0050]     With X&#39;s and V&#39;s components bounded, it is natural to represent X and V via vectors XI and VI whose components can take a finite number of integral values. For example, let XI and VI be vectors of conventional signed twos complement integers of bit width N_BX and N_BV respectively and X and V obtained from XI and VI from: 
 
X=XIDX V=VIDV 
 
         [0051]     where: 
 
 DX=XMAX/ 2ˆ( N   —   BX− 1) DV=VMAX/ 2ˆ( N   —   BV− 1) 
 
         [0052]     It can be proved that in such a representation, the following NVE equivalent integer math integrator is exactly reversible: 
 
 XI   —   h=VI   — 0+ROUND( Mˆ− 1 F ( X   — 0) DT/ (2 DV )) 
 
 XI   — 1= XI   — 0+ROUND( VI   —   hDVDT/DX ) 
 
 VI   — 1 =VI   —   h +ROUND( Mˆ− 1 F ( X   — 1) DT/ 2 DV )) 
 
         [0053]     The operation ROUND( D ) above numerically rounds D (whose representation is arbitrary) to the nearest integer (with ties broken toward the even integer).  
         [0054]     Implementation of simulation software using integer arithmetic can have one or more disadvantages. For example, most quantities of interest in MD are naturally represented as a floating point number. Furthermore, modern processors are typically heavily optimized floating point arithmetic. For example, support is provided for vectorizing floating point operations but little support is provided for vectorizing integer operations. Finally, mixing integer and floating point arithmetic is typically very slow, and operations such as ROUND and converting between floating point and integer numbers are very inefficient.  
         [0055]     In another embodiment of a time-reversible simulation, deficiencies of an all-integer embodiment can be alleviated. The basic steps for constructing such a reversible unconstrained NVE integrator are detailed below.  
         [0056]     First we store XI and VI as IEEE floating point numbers. Note that in a single precision IEEE floating point number, 1 bit is used for the sign, 8 bits are used to the exponent, e, (biased by adding 127), and 23 bits for the mantissa, f, (the fractional part of the significant, m). That is, m=1+f, and the floating point value is 2 e (1+f).  
         [0057]     Using floating point representations for XI and VI can be done without loss of information if it is noted that IEEE single precision can exactly represent all integers whose magnitude is less than or equal to 2ˆ24 because there are 23 bits of mantissa (with the most significant bit of the mantissa implied for normalized values) in the 32-bit floating point format. (This limit can be understood by considering the floating point representation of 2ˆ24+1 as a single precision IEEE floating point number, which results in an exponent of 24, and a mantissa of 2ˆ−24—but the mantissa has only a precision of 2ˆ−23 so the fractional part is truncated to 0, which yields the same representation as 2ˆ24.) Thus, this requires N_BX and N_BV to be less than or equal to 24. For double precision, the corresponding restriction is less than or equal to 53.  
         [0058]     Next we note that IEEE floating point specifies that basic binary operations (+,−,*,/) between floating point numbers produce the floating point result closest to the exact result with ties broken to the “even” floating point number in non-exceptional cases. Since the additions in the integer based Verlet algorithm are between bounded representable integers and produce bounded representable integer results, IEEE rounding requires these operations be performed exactly. Thus, floating point non-associativity will not affect the addition operations in the NVE integer math integrator.  
         [0059]     Next we note that between 2ˆ23 and 2ˆ24 in single precision (and 2ˆ52 and 2ˆ53 in double precision), the only numbers that can be represented are integers and that all integers in that range can be represented. If we restrict N_BX and N_BV to less than or equal to 23 for single precision (or 52 for double precision) and restrict the input to the ROUND operation lo less than 2ˆ22 for single precision (or 2ˆ51 for double precision) we can implement the ROUND operation with great efficiency. Namely, in single precision, if |D|&lt;2ˆ22 
 
ROUND( D )=( D+ (2ˆ23+2ˆ22))−(2ˆ23+2ˆ22) 
 
         [0060]     To explain, with the given restriction on D, the result of D+(2ˆ23+2ˆ22) is in interval [2ˆ23,2ˆ24]. All integers can be represented and only integers can be represented in this interval. Thus, the result of this operation must be ROUND(D)+(2ˆ23+2ˆ22) (exactly). Since, ROUND(D)+(2ˆ23+2ˆ22), (2ˆ23+2ˆ22) and ROUND(D) can all be represented exactly, the subsequent subtraction must be done exactly.  
         [0061]     In short, this exploits floating point non-associativity to do integer rounding for us in just two normal IEEE adds or subtracts. Note that ROUND can be implemented almost as efficiently if the above additional restrictions are not met.  
         [0062]     In another embodiment, if XMAX and VMAX are restricted to be powers of 2, then DX and DV are also powers of 2. So instead of multiplying and dividing the inputs to the ROUND operation by various factors of DX and DV, the algorithm can be written directly in terms of X and V with scale factors absorbed into a magic shift factors of the form VELOCITY_MAGIC=DV (2ˆ23+2ˆ22 ) and POSITION_MAGIC=DX (2ˆ23+2ˆ22).  
         [0063]     Thus, if X — 0 and V — 0 are position of the form XI — 0 DX and VI — 0 DV and the various restrictions on ranges for positions, velocities and what not are met, the following algorithm produces correctly rounded X — 1 and V — 1 of the same form and is exactly reversible if implemented in IEEE floating point arithmetic: 
 
 X   —   h=V   — 0+(( Mˆ −1 F ( X   — 0) DT +VELOCITY_MAGIC)−VELOCITY_MAGIC) 
 
 X   — 1 =X   — 0+(( V   —   hDT +POSITION_MAGIC)−POSITION_MAGIC) 
 
 V   — 1 =V   —   h +(( Mˆ− 1 F ( X   — 1) DT +VELOCITY_MAGIC)−VELOCITY_MAGIC) 
 
         [0064]     That is at a cost of two adds per phase space component update, the conventional NVE algorithm is rendered exactly reversible and no other modifications are necessary algorithms that use this integrator.  
         [0065]     A side benefit of this approach is that if any of the magnitude restrictions above are violated, the result is not catastrophic. The simulation will just not be fully reversible. Thus, exceptional results (abnormally large velocities for example) will be handled gracefully. (In practice, it is easy to insure that exceptional results occur so rarely that reversibility is maintained for any conceivable simulation.)  
         [0066]     As an implementation issue, it is important that an optimizing compiler not recognize that an expression of the form y=(x+MAGIC)−MAGIC can be algebraically optimized at compile time to be equal to y=x. A number of engineering techniques can be used, including (1) using function calls MAGIC( ) that return the constant value, (2) introducing a procedure call in the following sequence  
         [0067]     y=x+MAGIC;  
         [0068]     null_proc(&amp;y);  
         [0069]     y=y−MAGIC;  
         [0070]     where null_proc( ) is a procedure that receives the address of y but takes no action before returning, but nevertheless preventing the compiler from combining the assignment statements to y, and (;3) using an explicit compiler directive that prevents the optimization, such as  
         [0071]     volatile float tmp[1];  
         [0072]     tmp[0]=x+MAGIC;  
         [0073]     y=tmp[0]−MAGIC;  
         [0000]     which prevents the compiler from assuming that the value of tmp[0] doesn&#39;t change between execution of the two assignment statements.  
         [0074]     In using these ideas in a parallel MD code, the evaluation of F(X) should give a result that is strictly a function of X. Because various calculations needed to compute F(X) are continually being reassigned and reordered among the nodes, floating point non-associativity may cause the calculation of F to depend on other factors. This can be avoided by maintaining a consistent ordering of the particles and computations such that computed forces are unique functions of the particle positions and force field parameters. For example, in a distributed implementation, a node may receive elements of the operand X, such as quantities associated with various particles, in an order that can vary from run to run of the simulation, or can differ in a forward time versus a reverse time timestep. By defining a canonical ordering of the operations for implementing F(X) based, for example, on the positions of the particles, the same computation is performed regardless of the order in which the elements of the operand are received.  
         [0075]     Another approach to maintaining a consistent evaluation of F(X) is to make the relevant operations in the computation of F(X) resilient against such reorderings via the above rounding techniques (e.g. round the forces output by various computations to fixed point numbers as well).  
         [0076]     Another reason for maintaining consistent evaluation of a function, such as F(X), is in a distributed implementation in which it may be preferable to perform various computations redundantly at two locations rather than exchange a result computed at one location. Although high numerical accuracy may be possible with a numerically small error, any difference in the values of F(X) computed at different location can cause significant errors. For example, if two processors compute F(X) and compare the value to a fixed threshold to determine whether to perform a certain computation, with the desired result being that one and only one of the processors performs that computation, even a small numerical difference near the threshold could result in no processor or two processors performing the computation. However, by using a technique to provide consistent evaluation of F(X), the consistency is maintained across space in the distributed system as well as across time.  
         [0077]     In some examples, the precision of computation is not the same for all types of quantities. Similarly, the precision of a simulation to be reconfigured dynamically (i.e., with different precision at different time intervals, in different spatial regions, or for different types of quantities) without loss of efficiency. For instance, precision of computations can be adjusted according to characteristics of the computation system, or according to a state or other characteristics of the system being simulated. These ideas can also be used to investigate or assess the impact of precision on position, velocity, force and other MD variables on MD simulation quality without significantly degrading performance.  
         [0078]     In some examples, the variables are allowed to exceed the “maximum” values in a hybrid approach. As long as the values are within the maximum, the precision is determined by the MAGIC values. If the values exceed the maximum value, the bitwise reversibility or other desirable properties may be lost or may not be guaranteed, but the computations can continue with only some degradation in numerical properties. Therefore, such a hybrid approach can provided desired numerical properties in typical situations while providing robustness if values exceed their typical ranges.  
         [0079]     Other embodiments are within the scope of the appended claims.  
       Appendix  
       [0080]     Below is software source code for an embodiment of the techniques described above.  
                                                                                                                                                                                                                                                                                           /**********************************************************************       *       **        * POSITION PHASE SPACE FLOATING POINT REVERSIBILITY AND CACHE       COHERENCE       *       ***********************************************************************       *       */            /*    Do not touch this unless you really understand the floating point                cache coherence voodoo magic and are well versed in IEEE-754           standard, ANSI-C89 standard, POSIX-1 standard and common compiler           and floating point compliance issues. */            /*   FLOAT32_INT_MAX is the float32 number that is equal to the largest                exactly representable positive integer such that all positive           integers smaller than FLOAT32_INT_MAX are also exactly           representable. (Floating point adds and multiples performed on           integers less than FLOAT32_INT_MAX in magnitude are performed           _exactly_.) Note: In IEEE floating point, 2/FLOAT32_INT_MAX is           equivalent to the ANSI-C “FLOAT32”_EPSILON. */            #ifndef FLOAT32_INT_MAX       # if FLOAT32_SIZE==4            #   define FLOAT32_INT_MAX (16777216.)            # elif FLOAT32_SIZE==8            #   define FLOAT32_INT_MAX (9007199254740992.)       # else       #   error “Unknown floating point format. Define FLOAT32_INT_MAX.”       # endif       #endif            /*    POSITION_RESOLUTION is the number of representable positions across                a homebox. Note in normalized coordinates, the homebox has unit           length and the edges of the homebox are at +/− 0.5. For the           midpoint method math and coordinate translation math to be           reliable, this number must be a power 2 less than but not equal to           FLOAT32_INT_MAX. This insures that any given coordinate of an owned           particle can be exactly represented on its neighbors and that any           borderline midpoints can be unambiguously resolved. */            #ifndef POSITION_RESOLUTION       #define POSITION_RESOLUTION (0.5*FLOAT32_INT_MAX)       #endif            /*   The interval [POSITION_WHITE_MAGIC,2*POSITION_WHITE_MAGIC] is the                interval containing all positive floating point numbers such that           the spacing between all exactly representable adjacent numbers in           the interval is exactly 1/POSITION_RESOLUTION. POSITION_BLACK_MAGIC           is the center of this interval. Given the restrictions on           POSITION_RESOLUTION above, the smallest interval that might occur           is [1,2] and the distance from the center to the edge of the           interval is at least 0.5. */            #define POSITION_WHITE_MAGIC ( 0.5 *FLOAT32_INT_MAX/POSITION_RESOLUTION       )       #define POSITION_BLACK_MAGIC ( 0.75*FLOAT32_INT_MAX/POSITION_RESOLUTION       )            /*   INCREMENT_POSITION increment a rationalized position by an                _arbitrary_ dx to the nearest rationalized position assuming:           1. |dx| &lt; 1/2 (no increments larger than half a homebox)           2. X_old is a rationalized coordinate.           3. |X_new| &lt; 2 (the output position is within two homeboxes of the           center of the current homebox)           If these are true the operation is exactly (bit-level) reversible,           the rounding is exactly symmetric and the output is a rationalized           coordinate. If these are not true, the operation will still produce           reasonable results but may not be reversible any more. However, if           any of these are not true, the simulation is likely           self-destructing due to numerical instability.           How this works:           Since |dx|&lt;1/2 and the distance from POSITION_BLACK_MAGIC to the           edges of the interval [POSITION_WHITE_MAGIC,2*POSITION_WHITE_MAGIC]           is at least 1/2, the result forced assigment                magic[0] = dx + POSITION_BLACK_MAGIC                is in the interval           [POSITION_WHITE_MAGIC,2*POSITION_WHITE_MAGIC]. In this interval,           only numbers with the desired adjacent spacing can be represented.           Since POSITION_BLACK_MAGIC is exactly representable in this range           then, magic[0] is exactly set to:                fixp_round_even(dx) + POSITION_BLACK_MAGIC.                fixp_round_even(dx) can then be found as POSITION_BLACK_MAGIC can           be exactly subtracted off this quantity to yield an increment to a           rationalized position that will result in another rationalized           position. This increment is also exactly reversible as           fixp_round_even is symmetric under dx --&gt; −dx. */            #define INCREMENT_POSITION(x,dx) BEGIN_PRIMITIVE {       \                 volatile float32 _magic[1];            \                 _magic[0] = (dx) + POSITION_BLACK_MAGIC;            \                 (x) += _magic[0] − POSITION_BLACK_MAGIC;            \                } END_PRIMITIVE