SHAKE and the exact constraint satisfaction of the dynamics of semi-rigid molecules in Cartesian coordinates, 1973–1977

This essay traces the history of early molecular dynamics simulations, specifically exploring the development of SHAKE, a constraint-based technique devised in 1976 by Jean-Paul Ryckaert, Giovanni Ciccotti and the late Herman Berendsen at CECAM (Centre Européen de Calcul Atomique et Moléculaire). The work of the three scientists proved to be instrumental in giving impetus to the MD simulation of complex polymer systems and it currently underpins the work of thousands of researchers worldwide who are engaged in computational physics, chemistry and biology. Despite its impact and its role in bringing different scientific fields together, accurate historical studies on the birth of SHAKE are virtually absent. By collecting and elaborating on the accounts of Ryckaert and Ciccotti, this essay aims to fill this gap, while also commenting on the conceptual and computational difficulties faced by its developers.


Introduction
SHAKE is a computational technique for simulating, in Cartesian coordinates, the dynamics of complex atomic and molecular aggregates by introducing, and respecting over time, a set of holonomic constraints. 1 Although its beginnings date back to Atomique et Moléculaire), reaching maturity during a workshop on models for protein dynamics organized in 1976 by that institution. Their results appeared in the Journal of Computational Physics the following year (Ryckaert et al. 1977). 2 By enabling more streamlined computer simulations, SHAKE has given a major boost to fundamental research in computational physics, biology, and chemistry, as well as medicine and drug design. The continuing relevance of SHAKE to a wide range of disciplines lies in the need for multiple researchers to exploit realistic modeling of biopolymers or biomolecules, to study their internal soft modes (e.g., coupled torsional modes). Such slow modes govern the stable or metastable internal configurations of the system under investigation and they are the result of a subtle balance between internal and external forces, such as solvent effects. Slow modes can be approached within the framework of classical mechanics with molecular dynamics (MD) using atomic field models. An inherent problem, however, is that the primary structure of molecules up to biomolecules and biopolymers is maintained by stiff forces (e.g., bonding potentials, bending potentials) that impose a very short time step, whereas the dynamics of soft modes occur on much longer time scales. A natural idea is to freeze all fast motions due to stiff forces by replacing the pseudo-harmonic potential that appears in the primary structure model of the system with holonomic constraints while maintaining the force field governing the softer modes and other long-range internal forces. Yet, such a semi-rigid model leads to equations of motion that are highly complex and difficult to deal with. More specifically, as we shall see, the use of generalized coordinates for systems with many degrees of freedom leads to dynamic equations that are almost impossible to compute in practice, whereas the use of a Cartesian coordinate approach requires explicit constraint forces that depend intricately on both atomic positions and velocities. SHAKE solves the dynamics of semi-rigid molecules in Cartesian coordinates using a general and easy-to-implement numerical approach that simply requires providing as input the list of holonomic constraints that must be incorporated into the molecular model. This result represented a major accomplishment in the MD simulation of complex real systems. Indeed, the publication of SHAKE changed the direction of MD, contributing to the founding of some of the largest research projects in biochemistry and materials science (e.g., Car-Parrinello MD), while also reconfiguring approaches to rare events and attracting much interest from applied mathematicians. 3 "Given the importance of SHAKE for introducing any holonomic constraints on molecular dynamics simulations," commented the major architect of CHARMM, chemist Martin Karplus (born 1930), "knowing its history would be of great interest." 4 The purpose of this article is to outline an overview of SHAKE's development through a series of interviews with two of its creators. 5 To date, there is still a dearth of historical analysis devoted to SHAKE and MD simulations of biomolecules-a circumstance that calls for more attention from historians of science and technology, designating an area of study in need of immediate intervention. 6 The study presented here aims to stimulate new debates on the history and current developments of MD while providing an opportunity to reflect on how statistical mechanics has impinged fields such as chemistry and biophysics through the powerful means offered by computer simulations.

The backdrop and limitations of MD in the mid-1970s
Molecular simulation (MS) represents a set of computational techniques used for the mathematical modeling of polyatomic systems and the simulation of their collective behavior. The MS approach can be divided into two main classes: the MD and Monte Carlo (MC) methods. 7 The latter in particular constitutes a large set of computational techniques, of which those used in MS commonly belong to a variant of MC known as Metropolis MC. 8 While MD methods can be used to study the properties of both 5 As will be clarified throughout this article, attention must be paid to a specific terminological distinction between several concepts, all of which exploit the term "SHAKE." With reference to Ryckaert et al. (1977), the fundamental equation (2.7) expresses that the numerical trajectory must satisfy the constraints at the next time step of the MD algorithm (whatever that may be). This equation is adapted to a specific algorithm in equation (3.4) and is made more specific in equation (3.7), when all holonomic constraints are restricted to binding constraints only. A general numerical method for solving the multivariable quadratic equation (3.7) is given at the end of section 3 of the 1977 article. Section 5, on the other hand, contains an algorithm for solving equation (3.7) in such a convenient and easy-to-use way that today all MD packages dealing with (bio)polymers contain the routine based on this algorithm. This routine was given the name SHAKE, hence "SHAKE routine." Over the years, different names have been introduced to denote the different algorithms used to solve the dynamics of the general equation (2.7) that appears in the original article. For this reason, in July 1997, during an international school of physics on computer simulations of rare events and the dynamics of classical and quantum condensed-phase systems, Ciccotti and physicist Mauro Ferrario (born 1955) recommended using the expression "SHAKE optimization" to refer to any numerical method that solves equation (2.7) and "SHAKE routine" or "SHAKE algorithm" for the specific method described in section 5 of the 1977 article (Ciccotti and Ferrario 1998, p. 165;Ciccotti and Ryckaert 1986). In this paper, we will use the terminology "SHAKE equation," "SHAKE technique," or simply "SHAKE" to refer to the fundamental equation itself, regardless of the specific technique used to solve it, while we will retain the terminology "SHAKE routine" or "SHAKE algorithm" when referring specifically to the algorithm presented in section 5 of the original article. 6 In talking with Ryckaert and Ciccotti, it became apparent that a history of SHAKE's genesis had never been recorded and commented on in detail. The two also confirmed that although they had been interviewed by numerous scientists, no historian of science had yet taken the initiative to survey them in depth, and this circumstance contributed, with some urgency, to the motivation for writing this article. A useful overview, however, appears in . 7 Useful resources for a thorough introduction to MS are Frenkel and Smit (2002); Schlick (2010); Allen and Tildesley (2017); Tuckerman (2023). 8 The original MC method samples uniform numbers while the Metropolis MC method, with all its variants, employs instead a procedure called "Markov chains" to sample a generic probability distribution defined a priori. In computational statistical mechanics (specifically in MS), Metropolis MC, with all its variants, is generally used. Its variants have been developed to accommodate particular problems, generally associated with the fact that, without some additional refinements, the method, despite its power, would not work. equilibrium and dynamical systems, Metropolis MC techniques are generally adopted for the study of thermodynamic equilibrium, for which they can be at times superior to the former. That said, MD simulations represent powerful methods to "compute theory," i.e., to simulate and predict the observable dynamic behavior (mechanical and statistical) of real molecular systems, starting from the fundamental laws of physics. In principle, MD allows one to sample any molecular system at a fixed total energy and density, providing thermodynamic, configurational, and structural data as well as the dynamics of fluctuations useful to compute transport coefficients (e.g., diffusion, viscosity) and spectroscopic data (e.g., Raman spectroscopy, infrared spectroscopy).
It is not the purpose of this article to outline the history of MD, which, however, would deserve more historiographical attention. 9 We can say, however, that the method had its birthplace at the California Institute of Technology (Caltech), where a research group led by chemist Berni Alder (1925Alder ( -2020 focused on the computer simulation of an idealized system of hard spheres. Their results were co-published in 1957, coming about through the collaboration of physicist Thomas Wainwright  and computer engineer Mary Ann Mansigh Karlsen (born 1932) when Alder had already left Caltech and moved to work full time at the University of California Radiation Laboratory at Livermore. 10 The article he co-published with Wainwright, entitled Phase Transition for a Hard Sphere System, set out to solve "exactly (to the number of significant figures carried) the simultaneous classical equations of motion of several hundred particles by means of fast electronic computors [sic]" Wainwright 1957, p. 1208). After Alder and colleagues defined the main tracks, while working on idealized systems, the studies took off. Scientists began to approach real systems, in all their complexity. This trend was successfully developed by physicist Aneesur Rahman (1927-87), who distinguished himself by his simulations of real systems using empirically-derived Lennard-Jones (LJ) potentials. In his study, which appeared in 1964, he operated a CDC 3600 computer and noted the dependence on space and time of the self part of the density correlation function, which determines how slow neutrons are inelastically scattered from real, liquid argon (Rahman 1964). 11 Rahman was then a physicist at Argonne National Laboratory, and at this time he began collaborating with chemist Frank Stillinger (born 1934), a member of the technical staff at Bell Laboratories. 12 With the latter, Rahman would develop a rigid molecule model for 9 In addition to Battimelli et al. (2020) useful resources are Mareschal (2018); Battimelli and Ciccotti (2018); Lévesque and Hansen (2018); Mareschal (2019); Mareschal (2021aMareschal ( , 2021b. 10 In 1953 Alder became a consultant for what was then the University of California Radiation Laboratory at Livermore (now the Lawrence Livermore National Laboratory). In 1955 he became a member of the research staff. The first results of Alder and Wainwright's research were presented in 1956 in Brussels during a symposium on transport processes in statistical mechanics organized by physical chemist Ilya Prigogine  and appeared in Alder and Wainwright (1958). See also Battimelli and Ciccotti (2018); Ceperley (2020); Ceperley and Libby (2021). 11 Some reasons to help explain the time gap, from 1957 to 1964, in MD simulations research are discussed in Verlet (1987, pp. 6-7); Battimelli et al. (2020, pp. 60-61). Verlet's writing mentioned in this footnote appears in a CECAM pamphlet containing testimonies dedicated to Rahman, to access which it is recommended to contact CECAM or search the web for "In Memoriam: Aneesur Rahman 1927Rahman -1987 "Technical staff member" is the exact expression Stillinger uses to describe his title, as can be seen on his personal website. The title, of course, indicates a research position, and Stillinger held it from 1959 to water, to be treated in MD using six generalized coordinates. "Both he [Rahman] and I," remarks Stillinger, happened to attend the Liquids Gordon Conference in August 1969, and as a result of conversations there, we agreed that the idea of a joint simulation project had substantial merit. Bell Labs was convinced to provide financial support to Argonne to cover most of the computer costs. The first of several joint papers that was co-authored with Rahman, covering a wide range of aqueous molecular attributes, appeared in 1971 (Stillinger 2004, p. 19572). 13 In those years, an algorithm developed by physicist Loup Verlet (1931Verlet ( -2019, one of the pioneers of MD in Europe, was gaining momentum (Verlet 1967). Using Rahman's work on liquid argon as a starting point, and also availing himself of the guidance offered by mathematical physicist Joel Lebowitz (born 1930), Verlet compiled an algorithm that allowed him to integrate the Newtonian equations of motion of a system of as many as 864 particles in a relatively simple way. 14 His approach reduced the computation time by a factor of the order of 10, allowing for the accurate calculation of properties such as pressure, internal energy, and high-frequency elastic moduli for liquid argon at different temperatures and densities.
With this background in place, further studies appeared very quickly, initially focusing on water, but then also considering other systems (Ryckaert et al. 1977, p. 327, p. 341). 15 In polyatomic structures, it is necessary to account not only for the stretching of the interatomic bonds, but also for bending motions, which change the angle between the bonds, and for torsional motions, which alter the torsional angles (Allen and Tildesley 2017, p. 113). To avoid the use of very small time steps, to perform MD it is convenient to freeze at least the higher frequency vibrations using holonomic constraints. It is always possible, in principle, to construct a set of generalized coordinates that obey equations of motion in which the constraints do not appear anymore. 16 The pioneering simulations, in generalized coordinates, of a semi-rigid molecule treated as a polymeric structure were those of Ryckaert and physicist André Bellemans (born 1931) (Ryckaert and Bellemans 1975). In their investigation, the two presented: […] a preliminary report on an attempt to simulate liquid n-butane near its boiling point: this molecule may by considered as the simplest example of a flexible chain, with a time-varying configuration and it is of special interest to see how the existence of this internal degree of freedom affects the behavior  13 Consider also Rahman and Stillinger (1971). For more information regarding the 1969 Liquids Gordon Conference see Cruickshank (1969). 14 The number of particles processed by Verlet was the same as that studied by Rahman in 1964. of time-dependent correlation functions and transport properties. (Ryckaert and Bellemans 1975, p. 123).
Inspired by the MD simulation of water by Rahman and Stillinger, Bellemans and Ryckaert treated the CH 2 and CH 3 groups of linear alkanes as LJ units; the C-C bond and C-C-C bending angle were frozen by constraints treated in generalized coordinates, and a torsional potential (now known as the Ryckaert-Bellemans potential) was incorporated into the model to treat the single torsional degree of freedom of the system. 17 This work saw its origins around mid-1973 when Bellemans and his student Ryckaert applied themselves to the development of a semi-rigid model for n-butane that consisted of four point particles (i.e., CH 2 and CH 3 groups) and a single degree of internal rotation (i.e., a dihedral angle between two planes defined by three successive points) in addition to the usual three degrees of global rotation and the three translational degrees. 18 For this semi-rigid model with seven degrees of freedom, in 1973 Ryckaert developed an MD program in generalized coordinates that he applied to a liquid sample of sixty-four molecules. However, moving to longer alkanes with the same kind of approach-which was what Bellemans and Ryckaert hoped to do-was nearly impossible, since the need to explicitly derive Lagrangian equations of the second kind was an impractical task for a complex molecule. In other words, the generalized coordinates employed in the model placed serious limitations on the ability to expand the approach developed by the two. It was this kind of problem that lit the flame from which SHAKE was born, at CECAM, where two of the three researchers mentioned in the introduction-Ryckaert and Ciccotti-were operating and had the opportunity to meet and create a team. 19

The emergence of a European spotlight on MD studies
CECAM is a research center dedicated to the fundamental studies of atomistic and molecular simulations and their application to frontier problems in condensed matter physics, biology and chemistry. 20 Headquartered in Lausanne since 2009, the center was established in October 1969 in Orsay-at the Université Paris-Sud (Paris-XI), Bâtiment 506-through the efforts of quantum chemist Carl Moser , its founding director. 21 The idea of creating what would become CECAM took shape between 29 and 31 March 1967 during a meeting held in Blaricum, the Netherlands, 17 All this appears in Ryckaert and Bellemans (1978), which can be considered the manifesto of their contributions. 18 This semi-rigid n-butane model could be viewed equivalently as four particles connected by three rigid bonding constraints and two rigid bending angle constraints. 19 We will return to this point in more detail in one of the following sections. We can emphasize, however, that Ryckaert and Ciccotti were in contact with Berendsen, who was working at the University of Groningen. In mid-1977, the collaboration between Berendsen and Ryckaert became institutionalized, with the latter becoming a postdoctoral researcher of the former. 20 See www.cecam.org. 21 CECAM's premises were relocated twice. The center remained in Orsay until 1994, when it was transferred to Lyon. "After some turbulence that undermined its survival in Orsay, there were discussions about a new location that would guarantee CECAM's future and independence. Chemist Jean-Pierre Hansen (born 1942), who was at the time the director of the physics department at the newly-founded École Normale organized to discuss the establishment of a European center that could become an international reference in computational atomic and molecular studies (Berendsen 1990, p. 2). 22 The gathering was organized jointly by Moser andchemist Enrico Clementi (1931-2021), along with IBM World Trade. 23 "For obvious reasons," recalls chemist Wim Nieuwpoort (born 1931), who attended the meeting, "IBM Europe was interested in this initiative and proposed to hold a conference at their premises in Blaricum, offering administrative and scientific backing." 24 Also among the participants were chemists Paul Bagus (born 1937) and Geerd Diercksen (born 1936), as well as physicist Robert Nesbet (born 1930). "Herman Berendsen did not join the meeting," recalls Nieuwpoort, as he was an early user of NMR spectroscopy in proteins and not yet involved with MD. Around 1970 Carl Moser visited me in Groningen, and I introduced him to Berendsen who had become interested in MD as a new tool for studying biomolecules, soon to be joined by Wilfred van Gunsteren. From that encounter came a successful series of MD meetings in Orsay attended by an enthusiastic group of international scientists. 25 A typical CECAM workshop consisted of scientists who all shared a well-defined purpose, researchers with different expertise working together to achieve that purpose. 26 "In the beginning," explained Ciccotti, "this included spending a month or two crunching numbers, and establishing potential collaborations that would last over time." Footnote 21 continued Supérieure de Lyon (ENS-L) and director of research at ENS-L, went to great lengths to get CECAM an excellent invitation to Lyon, where the latter in fact ended up settling." Giovanni Ciccotti in conversation with the author. See also Battimelli (2019, pp. 21-24). 22 These accounts have been collated by CECAM in the form of a booklet available in the CECAM archives and also on the website by searching for "Recollections of CECAM for Carl." Consider also the testimonies of Jan J.C. Mulder and Wim Nieuwpoort, paying particular attention to p. 21 and p. 23, respectively. 23 Clementi's involvement as organizer of the meeting with Moser is the reason why the meeting was held on IBM's premises. IBM's official representative was Gerry Short, a not fully identifiable person who we know worked at the time for IBM World Trade in London. Physicist Per-Olov Löwdin (1916-2000 was also supposed to attend the meeting, but was actually represented by his colleague Jean-Louis Calais . This information was obtained in an e-mail exchange with Paul Bagus, who was the secretary of the event and who edited, along with Moser, the meeting report (unpublished). The author became aware of this document during a conversation with Bagus, but was unable to consult it. See also Nieuwpoort (1990, p. 23). 24 Wim Nieuwpoort in conversation with the author. 25 Wim Nieuwpoort in conversation with the author. 26 CECAM was originally located on the highest floor of the same building that housed the Centre Inter-Régional de Calcul Électronique (Interregional Center for Electronic Computing or CIRCÉ). This allowed CECAM to make use of the latter's computing resources (Connes 2022, p. 88), at an affordable price, subsidized by the Centre National de la Recherche Scientifique (the French National Centre for Scientific Research or CNRS). In addition to its workshops, CECAM boasted other activities, such as preparatory meetings and preparatory workshops, as explained in Berendsen (1990). More historical information can be found on the CECAM website, as well as in the pamphlets produced to commemorate Aneesur Rahman's passing and Carl Moser's retirement, respectively In Memoriam: Aneesur Rahman 1927Rahman -1987 and Recollections of CECAM for Carl. A well-grounded contribution is Battimelli (2019) that appeared as a brochure distributed to participants during CECAM's 50th anniversary celebrations held in Lausanne in September 2019. To obtain a copy of this paper, interested readers are invited to contact Giovanni Battimelli or the CECAM archives staff.
Discussion groups, on the other hand, were short meetings to examine a problem to be solved and plan a way to approach it collaboratively. Some discussion groups were very successful; others less so. They never lasted more than two or three days. Often a good discussion group would generate one or more workshops: one of their purposes was precisely to define possible topics for workshops. 27 From Ciccotti's recollections, it appears that during the early years of CECAM there was a lack of focus; furthermore, there was no follow-up on possibly good ideas developed during workshops or after visits by international researchers. Things changed in 1972 with the arrival at CECAM of physicist Gianni Jacucci (born 1943), who was presented to Moser by physicist Franco Bassani , then a professor at the Università degli Studi di Roma "La Sapienza." It is primarily to Jacucci that CECAM's research agenda aligned with MD simulations. 28 Credit must also be given to Jacucci's exchanges with physicist Dominique Lévesque (born 1938), who was then working at Orsay in Verlet's laboratory, a center distinct from CECAM. Through this connection, Jacucci obtained the code to run MD simulations by the classical Verlet method and began his full-scale research. 29 The Italian researcher also managed to create a positive working environment with the staff operating at the Centre Inter-Régional de Calcul Électronique (Interregional Center for Electronic Computing or CIRCÉ) and to establish fruitful collaborations with scholars from all over the world.
The first of the CECAM workshops devoted to MS, entitled "Molecular Dynamics and Monte Carlo Calculations on Water," was organized in 1972 at Moser's behest, following a discussion he had with Berendsen and Nieuwpoort (Battimelli et al. 2020, p. 92). 30 Jacucci was not involved in organizing this workshop in an official way, but his energy helped set the stage for the workshop to happen in a productive way. Jacucci, then a CECAM resident, helped make CECAM not only a gateway for many scientists but also a cohesive community characterized by stable and deep collaborations. His role in this regard, by the way, created friction with Verlet who, on the same campus, had intended to do the same thing without really succeeding. 31 Further providing a strong impetus for the 1972 workshop, upstream of Berendsen's leadership and Jacucci's contribution, was the article on MD simulation of water just published by Rahman and Stillinger (1971). This publication was a key element that contributed to a climate of enthusiasm, allowing Moser to launch at full force the CECAM activities toward the horizons of MS. The 1972 workshop began on June 19 and ended on August 11; the organizer of the event was Berendsen (1990). Rahman also participated, making his work available to the attendees and helping them to understand what it was all about. He, however, was in some ways very unconventional, concocting codes that were not always immediately understandable to his fellow scientists. "His programs and output listings," reported Berendsen,27 Giovanni Ciccotti in conversation with the author. 28 Giovanni Ciccotti in conversation with the author. 29 Giovanni Ciccotti in conversation with the author. 30 The exact title of the workshop is confirmed in Stillinger andRahman (1974, p. 1557). 31 For an account of the problematic relationship between Verlet and CECAM consider Hansen and Lévesque (1990, p. 31); Battimelli (2019, p. 11). were quite unintelligible for others: they did not contain a single comment or heading; indeed only numbers appeared on the output listing. He was extremely pragmatic and would not spend one moment on what he considered irrelevant.
[…] In addition to this inherent complexity, his programs were not precisely the best examples of good structured programming, but they worked and were error-free. (Berendsen 1987, p. 10).
Rahman's direct contact with the workshop participants-mostly younger researchers -were ultimately of great benefit to many.
The very first CECAM workshops were intended for very small groups of about twenty people. Among these researchers, those who had obtained useful results, or had promising problems to suggest, usually offered a lecture to propose a working theme to share. Working groups were then created, comprised of young researchers who had come to CECAM to work on a specific theme, usually accompanied by more senior researchers. Indeed, CECAM workshops were primarily open to senior scientists already deeply embedded in a certain research field. In the early days, Rahman was certainly the figurehead; also worth mentioning is chemist Konrad Singer (1917-2013, who together with physicist Ian McDonald (1938McDonald ( -2020, his researcher, would contribute to the early development of MS in the UK (Battimelli et al. 2020, pp. 77-83). Also in attendance were physicist and mathematician John Barker (1925Barker ( -1995 and chemist Douglas Henderson , both associated with the San José IBM Research Laboratory (later renamed the IBM Research Almaden). 32 Prominent scientists, in addition to Berendsen, included Bellemans and Orban, who took on the task of developing the topic of proton polarization in ice. 33 By carefully analyzing the history of CECAM (e.g., Battimelli 2019; Battimelli et al. 2020), reading the reports and scientific output, one is able to grasp the effervescent atmosphere of those years, animated by highly motivated young researchers who were facing an unexplored field of research. Thanks to the contributions of Rahman and Stillinger, and because it was clear by then that MD simulation of water could be successfully performed, higher aspirations began to emerge, aimed at approaching more complex systems (Berendsen 1987, p. 10). After the successes of 1972, in particular, two new workshops emerged, both held at CECAM in 1974. The first, organized by Singer, was entitled "Study of Ionic Liquids by Computer Simulation." It took place between May 20 and July 13 and focused on theoretical studies and computer simulations of molten salts and electrolyte solutions. 34 The second workshop, organized by Berendsen, was titled "Methods in Molecular Dynamics" and was held between August 19 and September 13. 35 The latter would build the intellectual 32 Barker and Henderson were well-known and respected specialists in the fluid mechanics and classical and quantum MS. Two articles they co-published in 1967 deserve special mention: Henderson (1967a, 1967b). By 1967, they were affiliated with the Division of Applied Chemistry of the Commonwealth Scientific and Industrial Research Organisation (CSIRO) Chemical Research Laboratories in Melbourne, Australia. 33 Jean-Paul Ryckaert in conversation with the author. 34 A short report of the workshop appears in CECAM's 1974 annual scientific report. Interested readers should contact CECAM to obtain the original document. 35 Special thanks to Giovanni Battimelli and Sara Bonella, who are credited with tracking down this information in the CECAM archives. In conversation with the author, Battimelli also noted that a third framework leading to the development of SHAKE. In this first phase, from 1972 to 1974, the protagonist was Ryckaert, a young doctoral student who had just graduated with a master's degree in chemistry from the Université Libre de Bruxelles (ULB). In the second and final phase, between 1974 and 1976, the team was composed mainly of Ryckaert and Ciccotti, while Berendsen, as we shall see, was at first more marginally involved-although his contribution ended up being substantial.
Ryckaert began his Ph.D. in chemistry at ULB in September 1972 under the supervision of Bellemans and spent the first three months defining a doctoral project on alkanes and securing three years of grants for this project. After learning that he had won his fellowship in late December 1972, he began his research activities in January 1973. Bellemans was then the head of a group on statistical mechanics of liquids, belonging to the ULB Physics Department, and Orban was his right-hand man. 36 It was the latter, however, who helped Ryckaert take his very first steps. In early 1973, the two began focusing on an MD program for an n-alkane chain with intramolecular and intermolecular forces and dynamics written in Cartesian coordinates. In their approach, stiff harmonic forces were taken to represent hard degrees of freedom, and thus there were no constraints. In addition, since the forces depended only on positions (and not velocities), Verlet's method was used to produce the MD trajectory. Realizing that the time step was very short, Ryckaert tried to move forward with a semi-rigid model in which C-C bonds and C-C-C bending angles would be frozen by constraints, while explicitly accounting for soft torsional degrees of freedom. "Having no idea how to proceed for a general alkane," Ryckaert commented in a conversation with the author, Bellemans and I decided to focus on the case of butane since it is the simplest semi-rigid alkane, i.e., a four-carbon chain with only one internal degree of freedom in addition to three translational and three rotational degrees of freedom. In other words, to find a way out, we temporarily reduced our ambitions related to the study of the whole series of n-alkanes to focus on the particular case of n-butane, which was nevertheless a remarkable case to treat. Ours, in fact, would have been the first known application of MD to a molecule capable of changing its internal conformation in relation to dihedral angle dynamics. 37 Footnote 35 continued workshop was held at CECAM in 1974, but it is not clear who was the organizer; it was entitled "Direct Methods Applied to the Determination of Protein Structure" and was held between July 15 and September 13. 36 The workshop report states that Orban and Ryckaert were then associated with the "Pool de Physique, Free University, Brussels." At the time of Ryckaert's doctorate, Bellemans was head of the group on Mécanique Statistique et Relaxation Diélectrique (the Statistical Mechanics and Dielectric Relaxation group) at ULB but also a member of a university association that grouped together research groups employed in physics: the Pool de Physique (the Pool of Physics). Members of the association included scientists such as François Englert (born 1932) and Jean Jeener (born 1931); it was created in 1961, primarily at the behest of Prigogine as well as physicist Paul Glansdorff (1904-99), the latter of whom was invited by the ULB science faculty to direct the Pool until 1975, when he retired. That same year, Bellemans became the Pool's coordinator, and his major task was to organize the teaching of physics at ULB. The Pool de Physique was dismantled in 1996, when Bellemans retired. 37 Jean-Paul Ryckaert in conversation with the author.
The idea at this point was to use generalized coordinates and the so-called Lagrangian equations of the second kind. Ryckaert then grew inspired by Rahman-Stillinger's work on water, in which a rigid model with three translational and three rotational degrees was used. The first step was to write the Lagrangian of the semirigid butane molecule in generalized coordinates. "I also made use," Ryckaert pointed out, "of what would become known as the Ryckaert-Bellemans torsional potential to ensure that the difference between the trans and gauche minima as well as the energy barriers were approximately realistic." 38 At this stage of his research, he made use of the Gear predictor-corrector method, also employed by Rahman and Stillinger in their 1971 paper. This choice allowed for high accuracy on the numeric trajectory; unfortunately, however, it caused a systematic drift in the total energy that could be empirically handled by keeping the time step sufficiently small. 39 "I was helped in the derivation of these equations by Jean van Craen," Ryckaert points out. 40 "I struggled a lot to get the MD program debugged, and then I began comparing both models on butane, the one with stiff springs and the one with constrained bonds and bending angles. 41 At ULB, Ryckaert was using a CDC 6500 computer that was somewhat inefficient at calculating long trajectories on the system of 64 butane molecules he was dealing with. 42 I remember that my runs were short (probably 100-500 time steps) and had to be done in sequence to get long enough trajectories and thus good statistics. There were very few trans-gauche isomerization events, and I was quite disappointed with all the divergences, e.g., energy drift, no equilibrium, and other difficulties related to using Lagrangian equations in generalized coordinates. 43 38 The Ryckaert-Bellemans potential was fixed by fitting the minima and maxima of the potential to data published in Scott and Scheraga (1966); Lipnick and Garbisch (1973). Although developed in 1973, the Ryckaert-Bellemans potential made its first appearance in Ryckaert and Bellemans (1975). 39 In dealing with constraints in Cartesian coordinates, Ryckaert used the Gear predictor-corrector method only in the first version of the program, the one produced at the 1974 workshop. The Gear method was later replaced by the Runge-Kutta fourth-order method, which allowed easy propagation of the system dynamics when it was necessary, as in Ryckaert's case, to reset the atomic positions to fully satisfy the constraint conditions every time the bond lengths deviated by more than 0.01% from their set value (in practice every 1000-10,000 steps). Beginning with the collaboration with Ciccotti, the Runge-Kutta method was abandoned in favor of the Verlet method, since this was more directly suitable for dealing with the SHAKE treatment of Lagrange multipliers.  Ryckaert and Bellemans (1975, p. 125). 41 Jean-Paul Ryckaert in conversation with the author. 42 "In 1973, when I started my PhD work," Ryckaert noted in conversation with the author, "the ULB/VUB computer was a CDC 6500." We should also mention that, prior to 1970, there was a single Université Libre de Bruxelles where most of the teaching was done in French, while some courses were taught in Flemish. In 1970, the University split into a purely French university, ULB, and a purely Flemish one, the Vrije Universiteit Brussel (VUB). The computer used by Ryckaert was located in a ULB/VUB Computing Centre (now more commonly known as ULB-HPC or simply HPC, where HPC stands for High-Performance Computing) that still exists today and has always been located on the ULB Solbosch campus. 43 Jean-Paul Ryckaert in conversation with the author.
The young student, in other words, was at a standstill, but soon an insight, born during a discussion with Bellemans and Orban, set him on a different path: I vaguely remember a discussion, probably in late 1973, or early 1974, between Bellemans, Orban and myself about the use of constraint forces along constrained links so as to avoid the use of generalized coordinates: this idea was reformulated more precisely during the 1974 CECAM workshop on long-time-scale events, which I had the opportunity to attend. 44 By the time he went to Orsay, Ryckaert was already planning to work on a new theoretical idea about how to treat constraint forces in Cartesian coordinates. At the 1974 CECAM workshop, he removed the rigid spring forces to include the constraint forces, while continuing to use the Gear predictor-corrector method, since these forces were velocity dependent. At that point, Ryckaert could count on the help provided by Orban.

The Bellemans-Orban-Ryckaert method, 1973-1974
The way for Ryckaert's participation in the 1974 CECAM workshop was prepared by Bellemans and Orban in 1972, thanks to the role the two played in the workshop on Molecular Dynamics and Monte Carlo calculations on water. In 1974 Bellemans and Orban were invited by Berendsen back to CECAM to join the workshop on the study of ionic crystals through computer simulations. The two, in fact, had formed a internationally regarded, and in the late 1960s they had also had a research experience in the United States. When it turned out that Bellemans was not free to go to the workshop because of too many commitments in Brussels, he proposed that Ryckaert replace him, since Orban was really eager to attend. On this occasion, Bellemans suggested that the two collaborate on incorporating constraint forces to replace stiff springs in the MD program (which they had developed earlier in Brussels) for n-alkanes of arbitrary length. As Ryckaert observed, André [Bellemans] was teaching very intensively in Brussels at that time and could not spend four weeks in Orsay. He then asked Herman Berendsen and Carl Moser if his presence could be replaced by a doctoral student who had a promising problem to solve related to long-time-scale events. Bellemans also made it clear that he would supervise such a student remotely and that the latter would work with John Orban, who was more flexible and could spend four weeks at CECAM. Incidentally, that student was me. 45 Ryckaert was then encouraged by his mentor to go to CECAM, where he would also benefit from stimulating scientific contacts and also the opportunity to work with a much more powerful computer than the CDC 6500. The work that Orban and Ryckaert did together during the 1974 CECAM workshop resulted in a paper, written primarily by Orban, that has never been published and contains descriptions of three case studies (Orban and Ryckaert 1974). 47 What the authors called "case A" is a realistic model for a chain of n-alkanes, in which all carbon atoms were considered to be connected by harmonic forces of realistic frequency. Then there were "case B1," restricted to n-butane, and "case B2," written directly for n-alkanes of arbitrary size, both treated with all bonds and bond angles considered constrained. The B-cases, in particular, corresponded to a model that was faster to integrate and therefore more suitable for computer simulations: in case B1, Orban and Ryckaert adapted and exploited the butane program on the IBM computer in Orsay, using generalized coordinates developed earlier in Brussels; case B2, on the other hand, was completely new, putting forth a bold attempt to use Cartesian coordinates. Regarding case B1, which required the use of appropriate equations of motion, it quickly became clear that the approach of using generalized coordinates would be essentially impossible in the case of more complex molecules; the approach seemed to work up to butane, which is characterized by four pseudo-carbon atoms along a chain (representing CH 2 or CH 3 ), but was practically inapplicable for longer chains. In case B2, the authors found a way to introduce constraining forces in addition to the canonical forces in the system. Their goal was to prove that case B1 and case B2 had comparable efficiency, i.e., to demonstrate the consistency of the MD program in generalized coordinates with a new approach of MD in Cartesian coordinates with explicit constraint forces. 48 Although this insight, in retrospect, turned out to be correct, in 1974 they could not convincingly demonstrate this equivalence, and this was the main reason why the two decided not to publish their results. This was also compounded by the fact that Orban was not particularly interested in publishing on this topic. 49 The idea then arose to pursue and refine case B2, as this was reasonably extensible to longer chains. "My concern," Ryckaert points out, "was that, although the programs worked correctly, the energy systematically drifted and that, in the new approach represented by the B2 case, the bond lengths and bond angles systematically shifted, albeit slowly." 50 In other words, the path shown by case B2 was the only viable one for approaching chains of any length, but the constraints necessarily had to be strictly reset about once every thousand steps (never at each time step) since the positions had to be reset to satisfy the constraints. This was because there was a slow but systematic drift in the satisfaction of the constraints. Having found a good topic to give input on, an interesting avenue had opened up.
Footnote 46 continued molecules at liquid density) required 3.25 CPU seconds per MD time step. The time step itself corresponds to 0.36 × 10 -14 s in real time. My runs were about 4000 steps, so the total corresponding time for one state point was 15 picoseconds. This required 3.6 h of computation time on the IBM machine." See also Connes (2022, p. 88). 47 For copies of the manuscript by Orban and Ryckaert, which the author obtained from Ryckaert, please contact the staff of the CECAM archives. 48 The direct equivalence check of case B1 and case B2 could only be performed on the n-butane case, since B1 was restricted to n-butane. 49 Jean-Paul Ryckaert in conversation with the author. 50 Jean-Paul Ryckaert in conversation with the author.
During this workshop, the bulk of the work done by Ryckaert consisted of adapting his MD codes (n-butane in generalized coordinates and n-alkanes in Cartesian coordinates) to the IBM machine in Orsay and to replace, in the program in Cartesian coordinates, the stiff spring forces with constraint forces. By the end of the workshop, the young researcher had managed to produce a FORTRAN program that worked, "more or less, for alkane chains of arbitrary length in which the constraint forces were computed in terms of coordinates and velocities" (Ryckaert 2013, p. 3339). This first phase-the one lasting from the beginning of Ryckaert's doctorate until the end of the 1974 CECAM workshop-led to a result that is occasionally referred to by specialists as the "Bellemans-Orban-Ryckaert method" (abbreviated here as BOR method). Ryckaert's stay at Orsay was productive, and the young Belgian researcher was invited by Moser to return to CECAM as a long-term visitor. After a short interim period at the ULB, Ryckaert returned to Orsay on October 1, 1974, and remained there for a period of nine months-until June 30, 1975. The original purpose of this new stay at CECAM was to work on producing MD data related to the physico-chemical properties of n-butane and n-decane. It was during this period, however, that he met Ciccotti, a circumstance which gave a new twist to his research.
In a relatively short time after Jacucci's arrival at CECAM, Moser realized that the former's approach to research-not only in terms of the direction of his interests but also his flair and ability to make connections-would provide a major impetus for CECAM's growth and gave him carte blanche to bring promising young researchers to Orsay. So it was that in 1973 Jacucci went looking for Ciccotti, then professore incaricato of history of physics at the Università degli Studi di Lecce (now Università del Salento). 51 The two were friends, as well as having been classmates at the Università degli Studi di Roma "La Sapienza." Jacucci handed Ciccotti a handout from a course taught by Verlet and intrigued him with discussions about MD, inviting him to join the CECAM research staff. On March 1, 1974, Ciccotti arrived in Orsay, where he remained until August 31, 1975. 52 During the 1974 CECAM workshop organized by Singer, the two worked together, with McDonald then joining the team. Their attention focused on the application of the so-called "method of small disturbances," developed by Ciccotti and Jacucci shortly before the workshop, to the study of electrical conductivity of some liquid alkali halides (Singer 1974, p. 9;Ciccotti et al. 1976). 53 Their results, computed by averaging the subtraction of the values of certain observables along two trajectories produced, respectively, in the presence of a very 51 The title professore incaricato no longer exists in the Italian academic system; it roughly corresponds to a contract professorship, a non-tenured associate professor. 52 Ciccotti graduated in physics from the Università degli Studi di Roma "La Sapienza" in 1967. After a period, which lasted from 1968 to the end of 1971, as a researcher with a CNR (Consiglio Nazionale delle Ricerche) fellowship in the high energy theoretical group at his alma mater, he joined the physics department at the Università degli Studi di Lecce to teach history of physics, holding this position until 1973. Ciccotti went to CECAM with the explicit intention of collaborating with Jacucci. By the time he arrived at Orsay, he had already moved from the Università degli Studi di Lecce to become professore incaricato of mathematical statistics at the Università degli Studi di Camerino, a position he held from 1973 to 1978. 53 To consult this document, it is necessary to contact CECAM, from which the author received a digital copy. We should note an inaccuracy reported by Singer, who writes that "conductance and shear viscosity of some liquid alkali halides" were addressed by Ciccotti and Jacucci during the 1974 workshop (Singer 1974, p. 9). This information was partially corrected by Ciccotti in conversation with the author. "Our small external field and in its absence, were considered extremely promising, as stated in documents from that time (Singer 1974, p. 9). The subtraction method, as well as the dynamical nonequilibrium molecular dynamics (Ciccotti and Jacucci 1975), later called D-NEMD (Mugnai et al. 2009;Ciccotti and Ferrario 2016), helped raise the bar for the young Italian researcher who had just arrived at Orsay and had previously tried other avenues: those of high energy phenomenology, many-body condensed matter, mathematical statistical mechanics, and the history of physics. It was in this way that Ciccotti and Ryckaert, both long-term CECAM residents, ended up becoming friends. This was late in 1974 or early in 1975, as Ciccotti recalls, although the two had first met in Paris at the closing dinner of the workshop organized by Berendsen in 1974, a reception Ciccotti had been invited to attend.
When Ryckaert began collaborating with Ciccotti in early 1975, he already had a working program on alkanes that could handle constraint forces. In particular, the BOR method allowed for the effective use of Lagrangian equations of the first kind, and in computing a priori the Lagrange multipliers in terms of positions and velocities, requiring a reset of positions due to the accumulation of algorithmic errors. 54 However, the method was far from satisfactory. 55 The main conceptual problem was related to the relaxation times required to achieve configurational equilibrium, thus demanding considerable computation time, which in the 1970s was a real stumbling block even for the powerful computer available to CECAM. Secondly, there was the problem of how to compare the results obtained from the simulations with the experimental data (e.g., how to deal with the link between the various relaxation processes obtained directly from the MD with various spectroscopy techniques). Among the technical weaknesses of the program was the drift of the total energy due to the use of a non-reversible integrator in time and the fact that the constraints were not perfectly conserved along the MD trajectory, due to the accumulation of algorithmic errors. The first technical Footnote 53 continued viscosity studies," he explained, "occurred a little after the end of the workshop (or perhaps even just before its end), when I found out that special perturbations can be used to excite the thermal or viscous response of a system. What Singer reported makes sense, but somehow places that specific development at an earlier time. Certainly in 1974 we still had no results in hand for this generalized approach. Our first article on the subject did not come out until 1978, while a second, more comprehensive study was published in 1979." See Ciccotti et al. (1978); Ciccotti et al. (1979). 54 The algorithmic error is the error of the approximation used to compute the trajectory of the particles under consideration and is normally expressed as the power of the integration step. 55 This point gives us the opportunity to clarify an issue that appears in the secondary literature, which found Ciccotti and Ryckaert in a standoff, with the former saying that Ryckaert was abandoned at CECAM by his supervisor and was working alone (see Battimelli et al. 2020, p. 106). "Bellemans never abandoned me at CECAM," stressed Ryckaert in conversation with the author. Although the verb "abandon" may perhaps be a bit too crude, one cannot help but notice, as Ciccotti argues, that Bellemans was not a fundamental presence in the development of SHAKE. Although Ciccotti and Ryckaert have differing opinions on this matter, it should be noted that Bellemans (as well as Orban) did appear in the acknowledgements of the 1977 article that sanctioned the birth of SHAKE. As Ryckaert explained in a conversation with the author, "At the [1974] workshop, I collaborated with John Orban, who was interested in the challenge we had before us. Very quickly, however-namely, as soon as we returned to Brussels-he no longer showed much interest in my doctoral research. In 1975 and1976, I worked only with Bellemans on the interpretation of the results on n-alkanes: there was no more contact with van Craen, who had moved to VUB and had no more collaboration with Orban, who, in turn, decreased all his research activity from then on. Therefore, in my opinion, Bellemans and Orban deserved to be recognized in the 1977 article for their contributions to the earlier BOR model, which was never published as such.". problem was solved by switching from the Gear predictor-corrector method to the Verlet method (having to iterate over velocities since the constraint forces depend on velocities); the second problem, however, was more thorny and was only partially addressed. As explained earlier, the issue was tackled by resetting the atomic positions in an ad hoc way so that the constraints were exactly met. 56 "From time to time," Ryckaert explained to the author, "I had to shift (reset) the particles slightly in space to re-establish perfect constraint satisfaction. The frequency of these resets was controlled by a tolerance of relative discrepancies of 0.0001 or 0.00001 on the bond length."

SHAKE and the MD of constrained polymer systems, 1975-1977
"In early 1975," Ryckaert recalls, "besides me, there were many scientific visitors at CECAM on temporary contracts, including many Italian chemists and physicists, such as Jacucci, Ciccotti but also Wanda Andreoni, Giuseppe Suffritti, and Mario Raimondi." 57 The specialties of these researchers were quite varied, all related to scientific problems requiring heavy numerical computation, ranging from quantum chemistry to liquid and solid-state physics.
We had many opportunities to meet in the cafeteria on the lower floor of the building where CECAM was located, especially at lunchtime. I was the only Belgian at the time, and it was quite natural for me to integrate myself into the Italian community working there. Most of us researchers were very busy with our own research programs and already channeled into various collaborations. In all this, Giovanni Ciccotti was quite different, very open and flexible, and this was because of his irrepressible curiosity: he was eager to know more about what I was doing, in detail. Personally, and despite the many visitors, I was a bit isolated in the scientific life of CECAM, and so I was happy with this unexpected and new scientific contact. This opportunity allowed me to explain to Giovanni what I was doing, including the various pitfalls I was encountering at that crucial moment in my doctoral research. The hope, of course, was that my connection with him would lead to some key indications on how to solve them. 58 An understanding immediately arose between the two. "Our first scientific discussion," Ryckaert recalled in conversation with the author, "lasted perhaps two hours, but I remember only two points very well. One of the first questions Giovanni asked me was: 'What is an alkane?' This question, so open-ended and out of the blue, puzzled me and confirmed my impression that he was indeed a pure statistical mechanics theorist." The second point is even more curious, revealing Ciccotti's very colorful and impulsive character. "When I told him that there was a drift in constraint satisfaction in my MD trajectories," notes Ryckaert, "and that I reset the coordinates from time to time to 56 As explained before (see footnote 39), the transition to the Verlet method took place at the time of Ryckaert's collaboration with Ciccotti and after Ryckaert had tried, without much success, to experiment with the Runge-Kutta fourth-order method. 57 Jean-Paul Ryckaert in conversation with the author. 58 Jean-Paul Ryckaert in conversation with the author. restore constraint satisfaction, he literally yelled at me: 'This is outrageous! It is a disgrace to treat the work of Newton and Lagrange in this way!'". 59 This comment, it must be understood, did not leave Ryckaert entirely surprised. First, his mathematical training was not comparable to that of a theoretical physicist; moreover, he was well aware of the limitations of his approach. Yet, he was eager to put his work on a solid theoretical basis. "Jean-Paul systematically reset distances," Ciccotti remarked. "I found that unacceptable. I wanted to understand how, in general, one could satisfy constraints: it's this need for the rigor that put us on the right track." 60 The collaboration between the two then found lifeblood in their interest in studying the mechanical laws underlying the MD algorithm, with a marked theoretical approach. "Needless to say, all aspects related to the properties of n-alkanes," Ryckaert points out, "as well as my main concern to interact with Bellemans were considered by Giovanni as totally uninteresting, being chemical." 61 This is a significant testimony that shows us how two young researchers, essentially novices and very different in terms of background and interests, eventually managed to form an efficient team. Indeed, the starting point of their collaboration created a rupture in Ryckaert's work. On the one hand, he continued to work on his doctoral dissertation with Bellemans: this was purely chemical work, as illustrated by their studies on n-butane and n-decane, which had a great impact at the time Bellemans 1975, 1978). The two of them also attended a number of international conferences in the UK and US, thereby creating the conditions for flourishing research in chemistry. "I remember very positive comments on our work on alkanes from John Rowlinson, David Chandler, Keith Gubbins, Doros Theodorou, and others," Ryckaert recalls, citing some of the most prominent names in the discipline. 62 The other avenue that originated in this bifurcation was represented by the interaction with Ciccotti on the theoretical roots of the MD algorithm with constraint forces. The collaboration with Ciccotti, moreover, would allow him to complete his doctoral dissertation: in fact, although the original MD algorithm had been envisioned by Bellemans in his 1973 discussions with Orban and Ryckaert, in early 1975 Bellemans considered it a useful tool for obtaining valid results on a semi-rigid n-alkane model, so that he could close this phase of his career and move on. 63 At this point, Ciccotti and Ryckaert realized that the mathematical background the latter was using was based on Lagrangian equations of the first kind, a topic they 59 Jean-Paul Ryckaert in conversation with the author. See also Ryckaert (2013, p. 3339). 60 Giovanni Ciccotti in conversation with the author. Ciccotti also adds: "So true is this [his penchant for researching principles] that I have used the idea of constraints to study rare events or, more recently, to see how to invent virtual zero-mass dynamics to deal with polarization or motions subjected to extra conditions." 61 Jean-Paul Ryckaert in conversation with the author. 62 Jean-Paul Ryckaert in conversation with the author. 63 "When Bellemans learned that I was beginning to collaborate at CECAM with a scientist from outside our group on the roots of the MD constraint algorithm," Ryckaert recalled in conversation with the author, "he was very pleased, but not to the point of devoting himself to the subject, being at that time about three hundred kilometers from us and very busy with new teaching and administrative tasks at ULB." Contacts between Ciccotti and Bellemans were always good: although the two did not collaborate directly during the writing of Ryckaert's dissertation, they did so later. Consider for example Ryckaert et al. (1981); Ryckaert et al. (1988); Ryckaert et al. (1989). both decided to refresh since they vaguely remembered it from their mechanics course when they had been students. Ryckaert had brought with him the textbook Theoretical Mechanics by physicist Ted Bradbury (1932Bradbury ( -1986, which he had used extensively in the second year of his chemistry degree (Bradbury 1968). 64 Reading chapter eleven of this text together, the two realized that the constraint forces are always orthogonal to the holonomic constraint configuration space surface, and that the Lagrange multipliers are actually regular and continuous-time functions, which Ryckaert was calculating precisely at each time step of the MD trajectory. This approach could be used for any algorithm, including those that would require third-order or higher-order timederived coordinates. "Giovanni and I were discussing mathematical formalism for hours when a striking development occurred as a result of a discussion we had with Singer at CECAM during one of his visits," Ryckaert recalls. On that occasion, Singer informed the two that he was performing MD simulations of rigid diatomic molecules of a given length l, avoiding the use of angular variables. He was, in other words, solving the rotational dynamics of the three-dimensional, constrained vector joining two atoms by first eliminating the external force components acting along the bond and then taking into account the centripetal force indirectly, through its effect; that is, he was forcing the relative vector to rotate on a circle of radius l, determining the amplitude of tangential motion with the external torque. In so doing, there was no longer any explicit dependence of accelerations on velocity, and the Verlet method could be used.
"Singer co-published his results (Singer et al. 1977) in parallel with our paper, so there are no cross-references," Ryckaert explains, "but he had shown his notes in 1975 not only to us but also, importantly, to Berendsen, who was in turn inspired by the work he was conducting." Singer's idea, applied to diatomic molecules, prompted Ryckaert and Ciccotti to impose the constraint conditions on the still unknown positions predicted by the algorithm at the time t + h. The result was a set of equations in the Lagrange multipliers that, when solved, provided the positions and velocities at the time t + h, while rigorously respecting the constraints. The Lagrange multipliers were no longer computed a priori (in terms of positions and velocities at the time t) but treated as adjustable parameters fixed by the satisfaction of the constraint at the next time step. "We had only bond constraints (2nd order in the positions)," remarks Ryckaert. Explicitly, on an n-alkane with N carbon units, the two imposed a set of 2N-3 constraints, thus having a set of 2N-3 gamma unknowns (i.e., the unknown Lagrange multipliers to be determined) and a set of 2N-3 order algebraic equations in the unknowns. "We had for n-decane a 17 × 17 matrix to invert," Ryckaert continues, "and we were solving the gamma with a rapidly converging interaction. This is what we later called the matrix method." That is to say, "Singer's idea was not directly generalizable," Ciccotti explains, "but it gave us the spark to figure out what we needed to do with our formalism. Jean-Paul and I jumped on it, generalized Singer's idea, wrote the SHAKE equation, and gave an iterative way to solve it." 65 This method was very efficient and robust and the two immediately thought of writing an article about it.
Significant developments continued in the period following the end of Ryckaert's long-term residency at CECAM through brief visits the latter made to Orsay, prompted also by his interest in continuing calculations on n-decane. 66 In spring 1976, both Ryckaert and Ciccotti were invited to participate to the CECAM workshop "Models for Protein Dynamics" organized by Berendsen, which was held in Orsay between May 24 and July 17, 1976. Around the beginning of the workshop, Berendsen approached the two, holding up a small bunch of punched cards bound by a rubber band, and said, "I have a routine here that solves the gamma for an arbitrary large molecule with an arbitrary set of binding constraints, avoiding the need to write and invert large matrices!" 67 In Berendsen's routine, one could simply solve a particular constraint as if it were isolated, correcting the positions at time t + h so that that constraint was satisfied. One could then do the same for the subsequent constraint, using the updated positions, and pass all the constraints in succession. One could then iterate the entire procedure over all the constraints until convergence was reached (see also Ciccotti and Ryckaert 1986). Given the many successive corrections of particle positions (i.e., one would take a cluster of positions and make a continuous readjustment, that is, a continuous change of positions), Berendsen called this routine "SHAKE." That same day he said it would be interesting if Ryckaert tested the SHAKE routine on his ongoing n-decane calculations and then checked whether his routine was equivalent to the matrix method that Ciccotti and Ryckaert had previously developed together. "I did the test in a few hours," Ryckaert notes, "and the results were 100% convincing that both iterative methods were leading to the same final positions." 68 "This particular procedure for solving the SHAKE equation," Ciccotti adds, "turned out to be excellent for polymers (and globular proteins are an example) and quite robust even for more complex cases: it was a different way than ours of solving the SHAKE equation." 69 "We were all very excited," continues Ryckaert, "Giovanni and I immediately proposed to Herman to join us as co-author of the draft we were working on, including in it the method of his routine, and he gladly accepted." 70 65 Giovanni Ciccotti in conversation with the author. The SHAKE equation is the equation obtained by inserting the coordinate predictions at time t + h into the constraint equation. This procedure completely cancels the error propagation. See footnote 5. 66 "Most of the results of my doctoral dissertation," stressed Ryckaert in conversation with the author, "were later compiled into the paper I co-published with Bellemans in 1978" (Ryckaert and Bellemans 1978). 67 Ryckaert in conversation with the author. This routine, as a matter of fact, had been written in FORTRAN by physicist Wilfred van Gunsteren (born 1947) in Groningen, applying an idea coming from Berendsen. Van Gunsteren had been until then a nuclear physicist and had just started a postdoc with Berendsen planning to work on the MD of proteins (holding this position from 1976 to 1978). In the 1976 CECAM workshop, which van Gunsteren attended, the two successfully incorporated the routine into an MD code that allowed a small protein (the bovine pancreatic trypsin inhibitor or BPTI) to be processed in vacuum, thus neglecting solvent effects (van Gunsteren and Berendsen 1977). 68 Jean-Paul Ryckaert in conversation with the author. 69 Giovanni Ciccotti in conversation with the author. 70 Jean-Paul Ryckaert in conversation with the author. After the protein dynamics workshop, Ryckaert completed in Brussels his doctoral dissertation (Ryckaert 1976). A few months later, he was contacted by chemical engineer James Haile (born 1946) of Clemson University, who suggested translating the Conventionally, we can place the benefits offered by SHAKE within two categories, which appeared immediately clear to its creators as well as to the scientific community engaged in MD simulations. Firstly, the Lagrange multipliers, calculated with the SHAKE technique, change the prediction of positions and velocities at time t + h to such a small extent that they remain within the algorithmic error. This process then is no longer an alteration, but rather a legitimate optimization. Moreover, by satisfying the constraints at all times, SHAKE does not propagate the error of the algorithm as far as the constraints are concerned (as was the case in the first attempts set by Ryckaert). This means that the SHAKE algorithm does not destroy the model being integrated over time. These two advantages represented, as is easily understood, a landmark result in the history of MD simulations. It was later shown that using the SHAKE technique in conjunction with a symplectic algorithm (such as velocity Verlet's) produces a resulting algorithm that is symplectic. 71 Moreover, all the variants of SHAKE that appeared in later years (e.g., RATTLE (Andersen 1983)) did not in fact represent a substantially new development, but rather new and often very refined applications of SHAKE, i.e., they were all based on the fundamental SHAKE equation and offered different methods for solving it. 72 Here, then, is the originality of the SHAKE technique, which is even today very topical to a large community of scientists.
The publication of SHAKE met with immediate success and resonance: Ryckaert and Ciccotti's work, complemented by Berendsen's contribution, brought to an end the era in which the main chemical bonds (i.e., those that directly build molecules up to polymers) had to be computationally treated in MD with hard pseudo-harmonic potentials. 73 This contribution allowed for enormous gains in computing time and expanded the research possibilities of laboratories all over the world. In fact, in 1985, physicists Roberto Car (born 1947 and Michele Parrinello (born 1945) broke new Footnote 70 continued thesis into English and took care of the translation. This work was completed by 1979, and the English text was reproduced on microfilm, although there appear to be no precise bibliographic references today. In September 1977, after nine months of military service, Ryckaert began a postdoctoral position under Berendsen at the Rijksuniversiteit Groningen (University of Groningen), which lasted until December 1978. In 1979 he joined the Faculty of Motor Sciences and the Faculty of Science at ULB, where he remained until his retirement in 2012. From 1974 to 2014, on the other hand, Ciccotti held various positions at the Università degli Studi di Roma "La Sapienza" becoming associate professor of molecular physics in 1980 and then full professor of structure of matter from 1990 to 2014. 71 A symplectic algorithm is defined as one that preserves the general explicit and exact properties of dynamics (e.g., the conservation of phase space measure, time reversal invariance). The solution of the SHAKE equation, no matter which method is applied to solve it, has been demonstrated to be symplectic when applied together with the velocity Verlet algorithm. For a demonstration of the symplectic nature of the velocity Verlet and the SHAKE techniques see Leimkuhler and Matthews (2015, pp. 161-164). 72 See footnote 5. An interesting anecdote, told by Giovanni Ciccotti to the author, is that the name RATTLE is not accidental, but refers to the song "Shake, Rattle and Roll" originally written by musician and songwriter Jesse Stone (also known as Charles Calhoun, 1901Calhoun, -1999 in 1954. The author tried unsuccessfully to contact physical chemist Hans Andersen (born 1941), creator of RATTLE, to ask for confirmation, but this anecdote could not be ascertained. 73 The 1973-1977 is the period of collaboration between Ryckaert, Ciccotti and Berendsen, which is the core of the evolution of SHAKE. We will now open up to an overview of a later period, which will only be sketched, spanning the entire period of substantial collaboration between Ciccotti and Ryckaert, which lasted until 1986, and then touch on some more recent developments. Note that Ciccotti and Ryckaert worked and published together, albeit more sporadically, even after 1986 (e.g., Ryckaert et al. 1988;Ryckaert et al. 1989;Depaepe et al. 1993;Pierleoni et al. 2015;Perilli et al. 2018). ground for SHAKE by developing classical dynamics to solve the problem in which electrons remain in the ground state and nuclei advance into the field of this state (Car and Parrinello 1985). In this approach, there were several holonomic constraints to satisfy, and the authors, who knew SHAKE, made use of it to compute the evolution of the dynamics within the constraints satisfied. SHAKE thus proved instrumental to the development of Car-Parrinello MD (Car andParrinello 1985, p. 2473; see also Marx and Hutter 2009, pp. 28-31, p. 34, p. 126). 74 Up to 1986, Ciccotti and Ryckaert, also along with other researchers, used SHAKE to study polymers, polyatomic fluids, and other complex systems (e.g., Ciccotti et al. 1982;Ciccotti 1983, 1986;Ferrario and Ryckaert 1985;Ryckaert 1985;Ciccotti and Ryckaert 1986). In 1989, Ciccotti co-published a seminal article in which the authors demonstrated that it was possible to use a holonomic constraint (the reaction coordinate, considered as a function of the configuration space) to force MD to explore regions that are naturally too rare to investigate (Carter et al. 1989). 75 This allowed for the calculation of activation energies of activated processes and rate constants of classical chemical reactions. This approach, called "constrained reaction coordinate dynamics" (CRCD) or "Blue Moon" (Carter et al. 1989, p. 474), was highly successful and led to the development of the theoretical procedure previously formulated by chemical physicists Charles Bennett (born 1943) and David Chandler (1944-2017 opening the door to a flourishing period of studies of rare events (i.e., activated processes) with simulation techniques (Bennet 1975(Bennet , 1977Chandler 1978;Rosenberg et al. 1980). Blue Moon was not based on SHAKE alone; indeed, the theoretical tool needed was the expression of an ensemble in Cartesian coordinates when the system is subject to holonomic constraints. The solution to this problem was presented by Ryckaert and Ciccotti (1983).
Further investigations conducted by physicist Wilfred van Gunsteren (born 1947) and Martin Karplus showed that constraining the molecule too much could alter the physics of the model; consequently, some pseudo-harmonic potentials, less hard than others, were saved (van Gunsteren and Karplus 1982). In 1992, physicist Mark Tuckerman (born 1964), chemical physicist Bruce Berne (born 1940), and chemist Glenn Martyna (born 1963) introduced the reversible multiple-time-step approach that allowed different integration steps to be used for degrees of freedom of different stiffness (Tuckerman et al. 1992). It seemed that this would spell the end of SHAKE, but the new method had inherent limitations and SHAKE eventually survived. In 1981, Ryckaert, Bellemans and Ciccotti had realized that one could use zero-mass material points bound by constraints to real interacting atoms to associate a three-dimensional reference with the diatomic molecules to study how rotations were related to translations (Ryckaert et al. 1981). The idea went unused for a while, having no further 74 In their article the authors wrote "[..] the values of the Lagrange multipliers were determined by the method of Berendsen" (Car andParrinello 1985, p. 2347). They, therefore, did not even use a variant, but remained bound to the original approach of Ryckaert, Ciccotti and Berendsen. 75 Ciccotti, continuing his research in the wake of SHAKE, had meanwhile opened up to a decidedly international perspective, participating in several CECAM workshops and spending extended stays at the University of Colorado Boulder (from October 1, 1987to November 30, 1987and then from October 1, 1988to December 31, 1988 and the University of Toronto (from December 1, 1987to September 30, 1988 applications. In 2013, working with physicists Sara Bonella (born 1971) and Alessandro Coretti (born 1992), Ciccotti began to realize that zero masses could serve to study the evolution of adiabatically separated degrees of freedom in a rigorous way (the Car-Parrinello method allowed this to be done in an approximate way). He then suggested applying SHAKE to polarizable systems and ab initio MD (Bonella et al. 2020). However, there are still works in progress.

Concluding remarks
From previously reported accounts, there are several salient phases that characterize the development of SHAKE. Early on, the Brussels team of Bellemans, Orban, and Ryckaert sensed the possibility of freezing high-frequency vibrations. They tried an initial approach that led to the so-called BOR method (indeed, just the explicit calculation of the constraint forces at the time t), which was for numerical purposes "sufficiently good, but largely empirical." 76 The missing ingredient, in hindsight, was an approach that could frame the problem in a more theoretical and general way. "The ad hoc resetting of the coordinates was one of the other weaknesses in my work," stressed Ryckaert in conversation with the author. By agreeing to tell Ciccotti about his research and the difficulties he was facing, he began a collaboration that took the problem on a more markedly theoretical course. Ryckaert and Ciccotti, in fact, first started down the path of analytical mechanics: they understood the importance of the Lagrange multipliers and that they can be calculated analytically together with their time derivatives. The two scientists then realized that, given an algorithm, it is possible to optimize it by staying within the algorithmic error, so that ultimately one can use any integration algorithm in the presence of holonomic constraints. This marks the transition, enshrined with the publication of SHAKE in 1977, from the initial empirical and heuristic approach to a general theoretical one. "Bellemans, Orban and Ryckaert," stressed Ciccotti in conversation with the author, "realized that instead of generalized coordinates it was possible to adopt Cartesian coordinates." Their contribution led to "rediscover" that the constraining force is computable. However, this was not enough. With SHAKE one does not calculate the Lagrangian multipliers of the constraint forces at time t, but at the time t+h. This circumstance actually leads to the constraint being satisfied exactly without introducing an error in the dynamics larger than that implied by the algorithm itself. This results in the numerical error not propagating, which prevents the unsatisfied constraints from ruining the dynamics. 77 However, there remained the problem, which needed an immediate solution, of formulating a particular ensemble (e.g., canonical, microcanonical), in Cartesian coordinates, if there were constraints to consider. In 1982 the two solved this problem, publishing it the following year, and showed how to perform isobaric MD generalized to molecular systems whose components are represented with constraints instead 76 Jean-Paul Ryckaert in conversation with the author. 77 Giovanni Ciccotti in conversation with the author. of vibrations (Ryckaert and Ciccotti 1983). 78 In 1985, Car and Parrinello, in need of satisfying holonomic constraints, published their landmark paper which emerged as a fundamental contribution to MS since it enabled an ab initio MD approach in which forces are calculated directly from the electronic contribution and not fitted on semi-empirical models (Car and Parrinello 1985).
The collaboration between Ryckaert and Ciccotti in the study of constraints lasted until 1986. Up to that year, Ryckaert recalls, "Giovanni and I, often together with other researchers such as Mauro Ferrario, have worked consistently on the SHAKE methodology for dealing with MD at constant pressure and/or temperature and on the formal statistical mechanical formulation of generalized ensembles." 79 At that point, Ryckaert lost interest in the purely theoretical question of constraints, while Ciccotti maintained the idea that constraints, again holonomic but more general than molecular, might be useful in further developing MD. This possibility came up in 1986 when the latter realized that constraints could bring a system to sample special rare events, setting the stage for the demonstration of the Blue Moon approach that appeared three years later (Carter et al. 1989). Meanwhile, the use of SHAKE in biological simulations began to find widespread use especially by researchers who used MD simulation packages dealing with biomolecules, usually proteins with water/ligand. 80 Nearly thirty years later, it was discovered that constraints can serve to rigorously implement adiabatically separated dynamics (Bonella et al. 2020). This circumstance found immediate application in handling the complexity of polarizable models, providing a tool for a revival of Car-Parrinello MD, which is a heuristic method based on an approximate use of the properties of adiabatic separation.

Conflict of interest
The author declares that there is no conflict of interest regarding the study herein presented.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.