A personal tribute to the founder of MD simulation of biological molecules: Prof Herman J.C. Berendsen (1934-2019)

On the 7 October 2019, Prof Dr Herman Johan Christiaan Berendsen passed away just shortly after his 85 birthday. Prof Berendsen is considered the founder of the molecular dynamics simulation of biological system: the area of theoretical research that also shaped my scientific career. He was working at the University of Groningen in the picturesque Northern part of the Netherlands. It was there that I meet him the first time as it allowed me to conduct research in his lab during the last year of my doctorate researches training at the University of Rome “La Sapienza”. After I completed my doctorate, Herman gave me the opportunity to continue working in his group with a postdoc position within the “Protein Folding” EU Training network. This happens just two years before his retirement and therefore I was also one of his last postdocs. After retirement, Herman dedicated himself to write two books that distillate all his experience in the area of molecular simulation [1] and in the education [2]. He stated in a project on the social scientific network Researchgate that “I am retired and work occasionally on methods for multiscale simulations.”

Prof. Berendsen will be remembered as a great scientist and a kind and generous person. I do not remember seen him getting irritated by anyone during group meeting, or in the famous joint BIOMOS meeting in Burg Arras (Germany) with the group of his historical collaborator Prof. Wilfred van Gunsteren from the ETH in Zurich. He could give you sharp criticism without let you feel an idiot but instead trying to help to solve your problem with insightfully suggestions. When I was working with him in Groningen, I remember that having meeting with him it was like going to have a good harvest of ideas: he was always ready to give you the best grapes but it was up to you get good wine out of it.

Prof. van Gunsteren in a essay dedicated to his 85 birthday recall the contribution given by Herman to the establishment of the computer simulation of biological system [3]. All started with the CECAM conference in Paris organized by Herman in 1976, gathering together 20 scientists, two of them awarded in 2013 with the Nobel Prize in Chemistry for their contribution to the field of computer simulations of biological system.

Since then Prof. Berendsen has given many contributions to the development of the molecular dynamics simulation of biological molecules. Here, I am just going to recall and describe some of the most important.

Berendsen Thermostat and Barostat

H. J. C. Berendsen, J. P. M. Postma, A. DiNola and J. R. Haak. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81 (1984) pp. 3684-3690.

According to Web of Science (accessed on Nov 2019), this is the most cited paper of Herman with 18242 citations.

In molecular dynamics simulation, you need to control the temperature and the pressure of the system to reproduce the same experimental conditions. In the 80s of the last century, Berendsen invented and simple and robust algorithm to achieve this purpose. He called this method the weak-coupling method to an external thermal bath but it soon becomes known ad the Berendsen method. So far, it is one of the most widely used  in  MD simulations. One of the advantages of this method is its simplicity. The algorithm simulates the coupling of the system to an external thermal bath with fixed temperature. The equations of the motion are modified in such a way that the current temperature of  the system T has a first order relaxation  towards the reference value T_0

\frac{dT}{dt} =  \frac{1}{\tau} [T_0-T(t)].

The intensity of the interaction between the thermostat bath and the system is determined by the time constant, \tau. Schematically the coupling is obtained by modifying the velocity of the atoms at each step.  This can be done by multiplying the velocities of  each particle by a constant \lambda. The new temperature is given by the relation

T'(t)= \sqrt{\left[1+  \frac{\Delta t}{\tau}\left(\frac{T_0}{T(t)}-1\right) \right]} T(t)

The weak-coupling method is also used in order to maintain the pressure of the system at a constant value. The equations of the motion are modified in such a way that the pressure of  the system has a first order relaxation  towards the reference value: 

\frac{dP(t)}{dt} = \frac{P_0-P(t)}{\tau_p}

the hydrostatic pressure P is defined for an isotropic system by the relation:

P(t) = \frac{2}{3}\frac{E_{K}(t) - \Theta(t)}{V_{s}(t)}

in which V_{s}  represents the volume of the system and the virial of intermoleculat forces \theta(t). The isothermal compressibility (\beta_T) is  the proportionality constant that relates the change of pressure  to the change of volume by the relation:

\Delta P(t) = -\frac{\Delta V(t)}{\beta_T V(t)}.

The change of volume can be obtained by multiplying it by a factor \mu:

\Delta V(t) = (\mu(t)^3 - 1 )V(t).

given by the expression

\mu(t) = \left\{ 1- \beta_T \frac{\delta t}{\tau_P} \left[ P(t)-P_0 \right] \right\}^{1/3}.

Practically this is accomplished by rescaling the coordinate of the atoms in the system, \vec{r}_i(t), using the same value of \mu

\vec{r}_i^N(t)=\mu \vec{r}_i(t)

According to the equation, the relaxation of the pressure $P(t)$ towards the reference value is controlled by the relationship between the isothermal compressibility of the system and the value chosen for relaxation time of the pressure (\tau_p).  It is not necessary to know the value of \beta_t exactly, since \tau_p is a parameter that can always be adjusted to obtain the proper reference pressure. Such a value should be sufficiently small (strong coupling) in order to obtain the equilibrium value of the pressure, but large enough (weak coupling) to avoid perturbing the properties of the system.

Development of the SHAKE algorithm

Ryckaert, J-P; Ciccotti G; Berendsen HJC (1977). “Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes”. Journal of Computational Physics23 (3): 327–341.

According to Web of Science (accessed on Nov 2019), this is the second most cited paper of Herman with 13482 citations.

For flexible molecules the length of the  integration step is limited by the frequency of the faster oscillator present in the molecule. This frequency is related to the type of vibration and to the mass of the oscillator. In fact, in order to integrate these motions correctly (assuring the conservation of the total  energy of the system), it is necessary to choice values of the  integration step that are very small. In the case, we are not interested to study vibrational properties of the system, than it is a good choice to remove these vibrations using  holonomic constraints. The use of these constraints has the advantage to allow the use of larger time steps (till to 10 times), speeding up the simulation. Nowadays as (and even more) at the dawn of MD simulation of biomolecules, extend the simulation time span allow to access the study of dynamic processes a occurring at different time scale. Berendsen and collaborator developed an algorithm, named SHAKE, to constraints bond distances based on the coordinate relaxation method. The SHAKE is an iterative procedure in which after each integration step, the positions \vec{r}_i(t+\Delta t) of atoms is varied until all the constraint conditions are satisfied, and a position \vec{r}_i(t+\Delta t) is found. This procedure is  performed correcting the velocities of the  constrained  atoms, after the application of the SHAKE algorithm, to eliminate velocity  components along the direction of the constraint.  The use of the SHAKE algorithm allows increasing the integration step in MD simulations of $2-4$ times.

Simple Point Charge (SPC) models of water

H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, The missing term in effective pair potentials, Journal of Physical Chemistry 91 (1987) 6269-6271.

According to Web of Science (accessed on Nov 2019), this is the most cited paper of Herman with 8023 citations.

Realistic simulations of biomolecules require the presence of an explicitly solvent model. The model of water parameterized by Berendsen and his group in Groningen is a simple rigid model. Still, capable of reproducing with reasonable accuracy the experimental density, enthalpy of evaporation and the variant SPC/E also the diffusion coefficient of water at 298 K. Given its simplicity, the models are still broadly used for simulation of biological molecules. The parameters of the model are reported in the following figure.

Figure 1: SPC and SPC/E water models

The GROMACS Package for MD Simulations

H. J. C. Berendsen, D. van der Spoel and R. van Drunen. GROMACS: A message-passing parallel molecular dynamics implementation. Comp. Phys. Comm. 91 (1995) pp. 43-56.

Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark A.E, Berendsen H.J.C (2005). “GROMACS: fast, flexible, and free”. J Comput Chem26(16): 1701–18

According to Web of Science (accessed on Nov 2019), these two papers are the sixth and fifth most cited with 5048 and 7048 citations, respectively. 

The name of this program is an acronym for Groningen Machine for Chemical Simulations. At beginning of 1990, Prof Berendsen initiated the development at University of Groningen (Netherlands) of a parallel computer called Groningen Machine for Chemical Simulations or Gromacs ring processor. It was a parallel computer based on the i860 transputers and specifically designed to perform MD simulations with high computational efficiency [4]. The specific parallel program for MD simulations written in the C programming language was also written by mainly David van der Spoel a talented PhD student now full professor at the Uppsala University in Sweden. The rapid evolution of low-cost parallel computing based on PC clusters made disadvantageous to keep continuing developing the specialized Gromacs hardware. Nevertheless, the Gromacs software package was further developed and implemented on different parallel architectures initially based on Silicon Graphics clusters. The first release of the program was the 1.6 and it was released on the 27/11/1997 (see in Figure 1 the original webpage and in Figure 2 an example of header comments in the source codes).

Figure 2: the first web page of the Gromacs 1.6 distribution.
Figure 3: Header comments in the source code file mdrun.c of Gromacs 1.6.

The contemporary developed of the Linux OS on ordinary PC rapidly spread the use of the program in the scientific community of molecular simulation. The introduction of new smart ideas for the calculation of the non-bonded pairs list, the virial of the forces and the use of assembler level optimized routines for the non-bonded interaction calculations on Intel and AMD processors enormously improved the performance of the program on affordable PC clusters and consequently its extensive use. In the spirit of the freeware movement of Herman wanted that the program remains public domain free available to everyone working in academics. Since 2001, the program has been licensed as GNU Public Licence (GPL) freeware and its source code is still freely available and open for development. Currently, David van der Spoel and his collaborators continue to maintain and develop it at Uppsala University with other contributors worldwide. As stated in the current documentation, the name Gromacs has not associated any more with its original acronym, another example of Herman’s scientific generosity and modesty. The official website of Gromacs is http://www.gromacs.org

Herman, with his group, has contributed to many other breakthrough achievements in MD simulation of a biological system that it would be too long to report here. I have co-authored 5 papers with Herman, and I am delighted that one of them is also my first published one [5].

There is no doubt that his legacy will continue to impact the field of MD simulations for many years to come.


  1. H.J.C. Berendsen. Simulating the Physical World: Hierarchical Modeling from Quantum Mechanics to Fluid Dynamics. Cambridge University Press, (12 July 2007).
  2. H.J.C. Berendsen.A Student’s Guide to Data and Error Analysis (Student’s Guides).Cambridge University Press, (7 April 2011).
  3. W. van Gunsteren. The Roots of Bio‐Molecular Simulation: The Eight‐Week CECAM Workshop ‘Models for Protein Dynamics’ of 1976. Helvetica Chimica Acta 102(8), 2019.
  4. Bekker, H., Berendsen, H.J.C., Dijkstra, E.J., Achterop, S., van Drunen, R., Van der Spoel, D., Sijbers, A., Keegstra, H., Reitsma, B. and Renardus, M.K.R., 1993. Gromacs: A parallel computer for molecular dynamics simulations. In Physics computing (Vol. 92, pp. 252-256).
  5. A. Di Nola, D. Roccatano, H. J. C. BerendsenMolecular dynamics simulation of the docking of substrates to proteins. Proteins: Struct., Funct., Genet., 19, 174-182 (1994).

2 thoughts on “A personal tribute to the founder of MD simulation of biological molecules: Prof Herman J.C. Berendsen (1934-2019)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.