The workflow of a MiMiC-based QM/MM BO-MD simulation using CPMD and GROMACS is illustrated in Figure (an analogous workflow would apply for CP-MD). Both CPMD and GROMACS are executed and run simultaneously. After initialization, MiMiC collects system data from GROMACS (number of atoms, atom types, coordinates, etc.) and propagates it to CPMD. Once in the MD loop, the first two steps of the velocity Verlet integrator are performed as usual. After the coordinates have been updated, they are distributed to GROMACS through MiMiC using the communication library (CommLib) which is part of the MiMiC framework (see below for further details). Then, both CPMD and GROMACS proceed computing energies and forces simultaneously. In the case of BO-MD, as in the example shown in Figure , this includes a self-consistent field (SCF) procedure to solve the KS equations in order to determine the electronic density. This requires also the electrostatic potential from the MM part (eq 3), which is computed by MiMiC and added to the KS-DFT external potential in CPMD. After the SCF, CPMD computes the QM forces and receives the QM/MM part of the forces (eqs 4-7) from MiMiC. Finally, MiMiC collects the energy and forces from GROMACS and propagates them to CPMD, which then performs the second velocity update and goes back to the beginning of the MD loop. In the present work, the equations of motion are integrated by CPMD, including those of the MM atoms. This implies that bond, angle, and other constraints among MM atoms must be enforced by CPMD. The (serial) RATTLE algorithm 5 implemented in CPMD as the standard constraint solver is suitable for systems of up to a few hundred atoms which is typical of pure QM systems but may become a computational bottleneck for larger systems (hundreds of thousands of atoms are not uncommon in MM). To overcome this problem, we implemented an improved version of RATTLE proposed by Weinbach and Elber. This improved version takes advantage of sparse linear algebra to impose the constraints, which greatly reduces the associated computational cost. Even though the method is in principle fully parallel, here we implement a version relying only on OpenMP parallel linear algebra. MPI and hybrid MPI/OpenMP versions of the method will be implemented in the future to further enhance the scalability.