Polymer glasses are non-equilibrium materials with complex mechanical behaviour. Although significant advances have been achieved in modelling deformation, yield and fracture of polymer glasses at the macroscopic level , connecting these properties to the chemical constitution and to the formation and processing history of a glass is still a challenge. This is because of the extremely broad spectra of characteristic times that govern molecular motion in the glassy state. Conventional atomistic simulation techniques, such as molecular dynamics (MD), employ an integration time step of approximately 10-15 s to track fast vibrational motions and thus can only address time scales up to 100 ns and length scales up to 10 nm with currently available computational means. As a consequence, they face two very serious challenges: (a) It is impossible to obtain a computer glass with a formation history that is both well-defined and realistic: MD vitrification experiments necessarily impose cooling rates of at least 108 K s-1, nine orders of magnitude higher than the rates of typical laboratory experiments. (b) Even if one is able to form molecular configurations that are truly representative of a real glass, an MD deformation experiment has to be performed at a strain rate of at least 106 s-1, which is much higher than those encountered in most applications.
In this work we discuss a strategy for addressing challenge (b) described above leading to a better understanding of yield and strain softening phenomena, as well as physical ageing. It is based on the idea that the local configuration of a glass is trapped in the vicinity of a local minimum of the energy, undergoing infrequent transitions to neighbouring minima across free energy barriers that may vary widely in height. This “energy landscape” picture focusses on the determination of representative energy minima and of the transition paths leading from those to neighbouring minima in the multidimensional configuration space of the glass. Thermodynamic properties and elastic constants in the individual energy minima are estimated by invoking a quasiharmonic approximation for the energy, and the corresponding properties of the glass are obtained through arithmetic (“quenched”) averaging over all minima. The rate constants for transitions from a minimum to neighbouring minima are estimated using the principles of multidimensional transition-state theory and the temporal evolution of the system, in the presence or absence of external stress, is tracked by Kinetic Monte Carlo (KMC) simulation as a succession of transitions between the minima. Such “quasi-dynamics” simulations can deal with arbitrarily slow transition rates and thus overcome the long-time problems of “brute-force” MD. The steps used in our approach could be described as follows:We focus on a small region of a polymer glass containing a few hundreds of atoms. We envision that the configuration of this region fluctuates in the vicinity of a local minimum of the potential energy, or “inherent structure” . Transitions between minima are largely inhibited by the presence of energy barriers that are high relative to kBT. Configuration space is thus partitioned into “basins of attraction”. The measured volumetric properties and elastic constants of the glass are shaped by the restricted probability distributions of configurations within individual basins. Physical ageing brings about a gradual redistribution among the basins through infrequent transitions across the energy barriers, and therefore a gradual change in the properties of the glass.
We invoke a quasi-harmonic approximation (QHA) by assuming that the potential energy of a glassy region of given spatial extent, while it fluctuates in the vicinity of its inherent structure, is well approximated by a Taylor expansion to second order around the inherent structure. At given stress, the glassy region will adopt that strainwhich minimizes G.
For the simulation of time-dependent plastic deformation and physical ageing phenomena, it is essential to determine all relevant transition pathways out of a given minimum. We begin by finding as many as possible saddle points around the minimum in the multidimensional configuration space of the polymer using the “dimer method” of Henkelman and Jónsson , which does not require second derivatives. In this procedure it is essential to recognize already visited minima on the basis of their energy and configuration. For each transition path between energy minima A and B, the rate constant is estimated according to transition-state theory, computed via the QHA as described above. Tracking the temporal evolution of the system thus reduces to a KMC simulation of a sequence of elementary transitions between basins, their rate constants being computed “on the fly”.
 Meijer, H.E.H. and Govaert, L.E. Prog. Polym. Sci. 2005, 30, 915-938.
 Stillinger, F.H. Science 1995, 267, 1935-1939.
 Kopsias, N.P. and Theodorou, D.N. J. Chem. Phys. 1998, 109, 8573-8582.
 Lyulin, A.V. Europhys. Lett. 2005, 71, 618-624.
 Henkelman, G. and Jónsson, H. J. Chem. Phys. 1999, 111, 7010-7022.