2.4. Molecular dynamics

Molecular dynamics (MD) simulations were performed with a single-precision version of Gromacs 2020.1.⁹³ The trimeric PSI complex was placed in a tetragonal simulation box with dimensions 26.0 nm × 26.0 nm × 14.5 nm. To minimize steric clashes introduced in the membrane packing process, the total energy was minimized with the steepest descent algorithm until the maximum force was below 1000 kJ mol⁻¹ nm⁻¹. In all following steps, the leap-frog integrator was employed with a time step of 2 fs and bonds to hydrogens were constrained using the LINCS algorithm.¹⁰⁶ Short-range electrostatics were calculated using Verlet lists with a cutoff distance of 1.2 nm. Particle-mesh Ewald summation¹⁰⁷ was conducted for the long-range electrostatics using cubic interpolation and a Fourier grid spacing of 0.16 nm.

Equilibration was conducted in three phases. First, the system was heated from 10 K to 100 K over 50 ps in an NVT ensemble, controlled by the V-rescale thermostat¹⁰⁸ with a time constant of 0.1 ps. The positions of the protein and all cofactors were restrained with a force of 1000 kJ mol⁻¹ nm⁻¹. Two temperature coupling groups were employed, one for the membrane-embedded protein and one for water and ions. The system was propagated for another 50 ps at a constant temperature of 100 K to further minimize clashes in the membrane and solvent. In the second equilibration step, the temperature was raised from 100 K to the production temperature of 300 K over 100 ps in the NPT ensemble, keeping the position restraints. The pressure was controlled by the Berendsen barostat with a semiisotropic reference pressure of 1 bar, a coupling constant of 5 ps and using an isothermal compressibility of 4.5 × 10⁻⁵ bar⁻¹. After annealing, the system was propagated for 900 ps at the target temperature in the NPT ensemble. In the third step, the position restraints were lifted and the system was propagated for 15 ns in an NPT ensemble, controlled by the Nosé–Hoover thermostat^109,110 and the Parrinello–Rahman barostat^111,112 with time constants of 0.5 ps and 5.0 ps, respectively. Proper equilibration was confirmed by the average temperature, pressure and density converging to their target values, as well as the total energy and box vectors remaining stable. The final box dimensions were 26.9 nm × 26.9 nm × 12.6 nm (Fig. 1). From the last nanosecond of this converged NPT ensemble, five production simulations over 15 ns were started. Snapshots were extracted only from the last 10 ns of each production run. During the production trajectories, the RMSDs of the protein backbone as well as that of the cofactors chlorophyll, β-carotene and phylloquinone remained stable and between 1 Å to 2 Å, further indicating that the system was equilibrated. Similarly, the area-per-lipid of the membrane had converged at 0.64 nm², in good agreement with literature values.^113–115