EchemDID
Testing with variable cut-offs:
/home/kristian/sshfs/rocket/echemtest/validation/forces_capa/test_cutoff
The default cut-off for the Coulombic interaction for the Qeq scheme used in EchemDID is 15 Å. Did tests with 20, 30,40.... 100 Å with the distance between plates being 18 Å. Found out that at higher cut-off values, the scheme breaks down entirely. The force and charge distributions become increasingly unphysical. Probably better to use the pre-determined cut-off value parametrized for the force field, even though this might cause the known problems.
Potential landscapes for different electric fields
Did tests with different electric fields by moving a test atom on 100 surface by one lattice constant. We can see, that greater the electric field, lower the potential energy of an adatom.
16.10 Minimization and dynamics of adatom; general review QeQ+Efield papers.
Ran the minimization of an adatom on the tiltsurface. It worked, but does not give anything that interesting. The forces go to 0 if the tolerance is low enough, as should be. At least it does not crash, meaning it is possbile to simulate such a system. Single atom temperature has no meaning - which should be well understood. Also thats the reason it is not possible or worthwhile to run dynamics with only the adatom unfrozen. Tried to run the whole tiltplane system at 300 K, which gave bogus results. Furthermore, more sensitive tests should be run only on minimized systems. Read https://doi.org/10.1016/j.ijsolstr.2007.09.010. It seems, that the electric field gradient is the driving force, although they look from contiinum-analytical perspective. Also it seems from their arguments, that the electric field gradient runs the diffusion anyways, without the necessary parallel components(dipoles)
Found several other methods which include the effect of the electric field. For example: https://doi.org/10.1063/1.5029877 Says, that they got rid of the cut-off problem. Hard to understand how from the paper, seems just a bit modified $\chi_{effective} -V/2$. They also have a fully filled system(no vacuum)
Also, the Van Duin paper for Ni/H2O interface supports the $\chi_i\rightarrow \chi_i + \mathbf{r_i}\cdot \mathbf{E}$ formulation. Although they have a tiny slab and fully filled, so the cut-off problem might not arise as heavily. They introduced additional force q*E also. This also makes sense (outside the cut-off the forces from the distant plane/point are not taken into account)
The main question: in the case of p-plate cap, is the effect of the electric field to change: $$\chi_i \rightarrow \chi_i \pm \frac{V}{2}$$
Or $$\chi_i \rightarrow \chi_i \pm \mathbf{r_i} \cdot \mathbf{E}$$
Where $V$ is the external potential difference, Both cannot be right at the same time. Or do they correspond to different physical models? If so, to which models they correspond? How do they interact with the finite cut-off? After all, the cut-off is there to make problems computationally feasible..
Comparisons of EchemDID and FEM
Did 3 comparisons with different geometries, some with adatoms, some without.
Triangular prism
Triangular prism in MD Constructed a right triangular prism with a side length of 70 angstrom. The width in the periodic direction was 35 angstroms. The surface on the diagonal was the [110] surface. The counter-electrode was constructed as a single point and set outside the cutoff of the QeQ , to mimic a bias voltage on the system with respect to infinity. An adatom was placed in the middle of the prism to allow to study the effects on the adatom within the same system, as it did not change the charge distribution significantly. voltage of 20V was applied to the prism.
Triangular prism in comsol Constructed a prism in comsol to match the dimensions of the MD system. To mimic the potential with respect to infinity, a large cylinder(100x the size of the prism) was created. Voltage of 20V was applied between the prism and the cylinder.
Comparison
The electric field values in MD were calculated in two ways: first way: [math]E=F/q[/math], where F is the electrostatic force acting on a surface atom and q is its charge. The other method is to look as the electrostatic force as a primary output, and the electric field as calculated: [math]E=\sqrt(2f/(s \varepsilon_0))[/math], where s is the mean atomic surface area, f is the electrostatic force acting on a particular atom. The electric field values in FEM in the middle of the prism, far from the edges were around [math]0.9 GV/m [/math] . The values in MD, calcuated from the first method were around [math]3.6 GV/m[/math]. The value calculated from the second method was display="inline">1.06 GV/m</math>, which is much closer to the field calculated by FEM. Conclusion: The electric field values calculated from the forces and mean areas give predictions much closer to the FEM values. The downside of this method is that in the case of deformed geometries, it might be hard to define a sensible mean surface area.