I have a NVT simulation with a constant external electric field using the option “electric-field-z = 1.5 0 0 0”
When calling gmx energy -f (…).edr to calculate the potential energy of the system, is the energy due to interaction with the electric field included or not?
Particularly, I have a lot of molecular dipoles in the system which can align wrt to the electric field. I want to know if the dipole - field interaction energy is contained in the potential from gmx energy.
In general it is not possible to define an energy, so this is not computed by mdrun. If you only have neutral molecules, you can compute this energy afterwards from the net dipole of the system. I forgot if mdrun always stores the dipole in the energy file. If not, you can use an analysis tool.
Note that you need to think about the choice of boundary condition at infinity (surface_epsilon mdp parameter) when using an electric field.
By default PME has conducting boundary conditions. This is fine with normal stimulations. But when adding an external electric field, it starts to matter a lot. With conducting boundary conditions dipoles will align in the field and periodically interact without penalty. With epsilon-surface=1 you get insulating boundary conditions and there is a large penalty to building up a net dipole. In the simple case of a dipole in a box of water, the electric field is amplified by a factor epsilon_r (about 80) when using conducting boundary conditions.
How can I determine the best value of epsilon for my situation? I have a kind of homogeneous crystalline solid. So should it just be the relative permittivity of the substance?
In general insulating boundary conditions gives the correct average polarization and conducting boundary conditions gives the correct fluctuation of the polarization.
But in most cases you can only use conducting boundary conditions as mdrun only supports epsilon-surface>0 when each molecule consists of a single update group (so it doesn’t get broken over PBC). You can get the correction factor for the electric field for use with conducting boundary conditions by computing the dielectric permittivity of your system using gmx dipoles.
Just to clarify: you mean that if I apply a field E with conducting boundary conditions, and the relative dielectric permitivitty (epsilon.xvg from gmx dipoles) is e_r, the “actual” field felt by the system is e_r * E?
I would like to understand a bit better how the electric field has been implemented in GROMACS.
Could you please refer the relevant files that I have to check to proceed with?