MD simulations were performed with the GROMACS 2016.4 suite, using the GAFF force field obtained from the AmberTools 18 package with the AM1-BCC charge model. Each simulation was first energy minimized using the steepest descents algorithm. FS and PA molecules were then positionally restrained using the LINCS algorithm, while the solvent molecules were allowed to move freely and relax into equilibrium positions. The system was preequilibrated for 1 ns in the NPT ensemble with the velocity-rescale Berendsen thermostat at 120 K, with a temperature coupling constant of 0.1 ps, and an isotropic Berendsen barostat applied with a reference pressure of 1 bar and a pressure coupling constant of 20 ps. The positional restraints on the FS PA were then removed after the solvent molecule relaxation, and the whole system allowed to equilibrate for an additional 1 ns with the same protocol as above, with the temperature set to 120 K, as in the original crystal structure determination. This initial equilibration was followed by an annealing run heating the system from 120 K to 350 K, in steps of 30 or 50 K over 10 ns, followed by 10 ns relaxation to prevent a hysteresis before the next temperature increment. Hheteroatom bonds were constrained using the LINCS algorithm, other parameters as above. Structures were extracted at the temperatures of 150 K, 200 K, 250 K, 300 K and 350 K, and then simulated in the NPT ensemble with the Nosé-Hoover thermostat at the given temperature with a temperature coupling constant of 1 ps, and an anisotropic Parrinello-Rahman barostat employed at 1 bar, with a pressure coupling constant of 20 ps. Production simulations runs were performed for 200 ns without bond constraints. The 150 K simulation was extended to 400 ns for improved sampling statistics.