2a. MD computational details. As done in our previous works , all MD simulations were performed with an MD code parallelized to fit the embarrassingly parallel nature of the crystallization trajectories . The details of the code can be found in the literature , and can either be obtained on github or through contacting the authors. Short-range non-bonded forces were truncated at 1.2 nm. Reciprocal space electrostatic forces were computed using a particle mesh Ewald (PME) solver with a 1 x 10 -6 error tolerance. 100 grid points were used in each direction for the PME solver. The electric field augmented equations of motion, equation (2), were solved using a velocity Verlet algorithm with a 2 fs timestep to generate atomic trajectories. The temperature was controlled using a Langevin thermostat with a relaxation time of 2000 fs. The isotropic Berendsen barostat was used with a 1000 fs relaxation time and a set point of 1 atm. A value of 21454.85 atm was used for the bulk modulus of water to determine the relaxation time constant in the Berendsen barostat. The shake algorithm was employed to constrain the motion of all bonds containing hydrogen. A shake tolerance of 1 x 10 -4 was employed for all simulations, and a maximum number of 25 shake iterations were allowed. Periodic boundary conditions were employed in all directions. All production simulations were performed at 330 K. 10 trajectories were simulated for each polymorph and electric field intensity to provide average growth and dissolution system dynamics.