All simulations were performed using the Groningen machine for the chemical simulation (GROMACS) program (version 5.1.4). The modeled systems were first energy minimized using the steepest descent algorithm, then equilibrated in an NPT ensemble at constant number of molecules N, temperature T, and pressure P, using a Berendsen weak-coupling thermostat and barostat algorithms. In the production run, a constant temperature of 310K, if not otherwise specified, was maintained using velocity-rescale thermostat with a time constant of 0.1 ps. The Parrinello-Rahman barostat was applied for isotropic pressure coupling with a reference pressure of 1.0 bar. The PME (Particle-Mesh-Ewald) method was used for computing long range electrostatic interactions with a cutoff distance set to 1.4 nm, Fourier spacing of 0.12 nm, and cubic interpolation. The van der Waals interactions were cut-off at 1.4 nm. The Verlet algorithm was employed to list all pairs within the cutoff, beyond which interactions could be ignored. The equations of motion were integrated using the leapfrog scheme except in the case of alchemical calculations where the SD integrator was used instead. Periodic boundary conditions were applied in the three dimensions. A timestep of 2 fs was used.