Efficient Computer Simulation of Polymer Conformation. I. Geometric Properties of the Hard-Sphere Model

A system of efficient computer programs has been developed for simulating the conformations of macromolecules. The conformation of an individual polymer is defined as a point in conformation space, whose mutually orthogonal axes represent the successive dihedral angles of the backbone chain. The statistical-mechanical average of any property is obtained as the usual configuration integral over this space. A Monte Carlo method for estimating averages is used because of the impossibility of direct numerical integration. Monte Carlo corresponds to the execution of a Markoffian random walk of a representative point through the conformation space. Unlike many previous Monte Carlo studies of polymers, which sample conformation space indiscriminately, importance sampling increases efficiency because selection of new polymers is biased to reflect their Boltzmann probabilities in the canonical ensemble, leading to reduction of sampling variance and hence to greater accuracy in given computing time. The simulation is illustrated in detail. Overall running time is proportional to dj4, where n is the chain length. Results are presented for a hard-sphere linear polymer of n atoms, with free dihedral rotation, with tz = 20-298. The fraction of polymers accepted in the importance sampling scheme, fA, is fit to a Fisher-Sykes attrition relation, giving an effective attrition constant of zero. fA is itself an upper bound to the partition function, Q, relative to the unrestricted walk. The mean-squared end-to-end distance and radius of gyration exhibit the expected exponential dependence, but with exponent for the radius of gyration significantly greater than that of the end-to-end distance. The 90% confidence limits calculated for both exponents did not include either 6/15 or 4/3, the lattice and zero-order perturbation values, respectively. A self-correcting scheme for generating coordinates free of roundoff error is given in an Appendix. he Monte Carlo method’ was firmly established by Wall T and his successors2-1O as a valuable tool for investigating * Address correspondence to this author at the Department of Biochemical Sciences, Frick Chemical Laboratory, Princeton University, Princeton, N. J . 08540. (1) J. M. Hammersley and D. C. Handscomb, “Monte Carlo Methods,” Wiley, New York, N. Y. , 1964. (2) (a) F. T. Wall, L. A. Hiller, Jr., and D. J. Wheeler, J . Chem. Phj,s., 22, 1036 (1954); (b) F. T. Wall, L. A. Hiller, Jr., and W. F. Atchison, ibid., 23, 913, 2314 (1955); 26, 1742 (1957). (3) F. T. Wall and J. Mazur, Ann. N. Y . Acad. Sci., 89, 608 (1961). (4) (a) F. T. Wall, S. Windwer, and P. J. Gans, J . Chem. Phys., 38, 2220 (1963); (b) ibid., 38,2228 (1963). the geometric and thermodynamic properties of polymers. The majority of previous Monte Carlo studies have been confined to the exploration of extremely simple models, such as random walks on lattices. In this paper we offer a more complete model of polymeric systems which includes lattice as well as off-lattice polymers on special cases, an ideal frame(5) J. Mazur and F. L. McCrackin, ibid., 49, 648 (1968). (6) E. Loftus and P. J. Gans, ibid., 49,3828 (1968). (7) P. J. Gans, ibid., 42, 4159 (1965), and references cited therein. (8) C. Domb, Adcan. Chem. Ph.l,s., 15,229 (1969). (9) J. Mazur, ibid., 15,261 (1969). (10) I<. I<. Knaell and R . A. Scott, J . Chem. Phys., 54, 566 (1971). Voi. 5, NO. 4 , Jdy-August I972 COMPUTER SIMULATION OF POLYMER CONFORMATION 5 17 work for Monte Carlo calculations and the use of importance sampling for the reduction of sampling variance. The remainder of this paper is devoted to elaboration of the model and method and presentation of detailed geometrical results for the hard-sphere off-lattice model. Distribution function properties of this model, and temperature-dependent properties of soft-core potential models are discussed in future papers of this series. The primary objective is the calculation of statistical-mechanical configuration integrals of the form where P is the value of some conformational property, such as the squared end-to-end distance, and U is the molecular potential energy. (P) denotes the canonical ensemble average of the property P. Numerical polymer theory depends on representation of the energy of the molecule in terms of structural parameters. Using the adiabatic treatment of Scheraga and coworkers,11 the conformational energy is decoupled from high-frequency bending and stretching modes, permitting a classical description in terms of nonbonded and rotational barrier potential energy functions, with semiempirically determined parameters. As shown by G o and Scheraga,12 when bond lengths and bond angles are held constant, the conformational partition function can be written Q = (constant) s. . s exp( U/kT)dr ( 2 ) where T i s the absolute temperature and U is the potential energy of the conformation. The element d r represents the minimum volume element necessary for this representation and defines the subset of phase space correspondiRg to the “conformation space” of the molecule. With fixed bond lengths and bond angles, the only remaining degrees of freedom are the dihedral angles about successive backbone-chain bonds. For n atoms in the backbone, n 3 dihedral angles define the relative position of every atom to every other atom (ignoring multiatomic side chains). These angles can be considered independent variables, and from them the coordinates and potential energy can be calculated using classical potential energy functions.I3 The ( n 3)-dimensional vector o, whose elements are the individual dihedral angles mi, defines a single conformation; all conformations are represented by the set of all possible vectors, { 0 ) . The volume element in eq 2 can be written dr = J(q)dwidwi. . , dWn--a where the range of each wi extends from 0 t o 29 . J(q) is the Jacobian for the transformation from generalized coordinates to the w representation. As a simple numerical factor,14 it cancels out of all equations for averages in which Q I appears, such as eq 1, and hence will be ignored. Lattice calculations restrict the values of the individual w.I)s to a discrete set of suitable values. For example, the tetrahedral lattice (11) H. A. Scheraga, Adran. P h j . ~ . Org. Chem., 6 , 103 (1968); Chem. (12) N. Go and H. A. Scheraga, J . Chem. Phys., 51,4751 (1969). Rec., 71, 195 (1971). (13) G. N. Ramachandran and V. Sasisekheran, “Conformations of Biopolymers,” G. N. Ramachandran, Ed., Vol. 1, Academic Press, London, 1967, p 283. (14) J. E. Mayer and M. G. Mayer, “Statistical Mechanics,” Wiley, New York, N. Y . , 1940, p 230. results from fixing bond angles a t 109.47‘ and taking W , E (60”, 180”, 300”). For any given structural model, the choice of energy function U ( o ) determines the degree to which the properties of the random walk mimic the conformational properties of a real polymer. The simplest energy function occurs when U is identically zero, for which the polymer dimensions correspond exactly with the statistics of the unrestricted random walk, whose mean-squared end-to-end distance is given (for large n) byl5,16