Implement surface dipole correction

.github/workflows/makefile.yml

@@ -142,6 +142,12 @@ jobs:
         mpirun -np ${{matrix.np}} ../../bin/Conquest
         cat Conquest_out
 
+    - name: Run test 008
+      working-directory: ${{github.workspace}}/testsuite/test_008_surface_dipole
+      run: |
+        mpirun -np ${{matrix.np}} ../../bin/Conquest
+        cat Conquest_out
+
     - name: Check test results
       working-directory: ${{github.workspace}}/testsuite
       run: pytest test_check_output.py

docs/groundstate.rst

@@ -410,6 +410,34 @@ beyond the valence electrons.
 
 Go to :ref:`top <groundstate>`.
 
+.. _gs_surf_dip:
+
+Surface dipole correction
+-------------------------
+
+If your simulation involves a slab calculation, then unless the slab is
+perfectly symmetrical there can be a dipole moment which gives an unphysical
+field across the periodic cell even for a neutral system.  A correction can be applied, following the
+method outlined in :cite:`g-Neugebauer1992,g-Bengtsson1999` (where we use the energy from the
+Bengtsson paper, which is correct).  This method *only* works for slab
+calculations, and calculates a dipole correction potential, placing the
+necessary discontinuity in the vacuum.  The user can specify the location
+of the discontinuity (in *fractional* coordinates), or the code will place it at the point where the
+planar averaged density is a minimum (which is a good way to find the
+centre of the vacuum).
+
+::
+
+   SC.SurfaceDipoleCorrection T/F
+   SC.SurfaceNormal x, y, or z
+   SC.DiscontinuityLocation (*real*)
+
+It is also possible to write out the planar-averaged potential and charge density
+(averaged in the plane perpendicular to the surface normal) using the parameter
+``SC.OutputAveragePotential T/F``.  This generates the file ``AveragedPotential.dat``.
+
+Go to :ref:`top <groundstate>`.
+
 .. _gs_spin:
 
 Spin polarisation

docs/references.bib

@@ -679,3 +679,23 @@ @article{Boys:1970aa
   year = {1970},
   doi = {10.1080/00268977000101561},
   pages = {553}}
+@article{Bengtsson1999,
+  title = {Dipole correction for surface supercell calculations},
+  author = {Bengtsson, Lennart},
+  journal = {Phys. Rev. B},
+  volume = {59},
+  issue = {19},
+  pages = {12301--12304},
+  year = {1999},
+  doi = {10.1103/PhysRevB.59.12301}
+}
+@article{Neugebauer1992,
+  title = {Adsorbate-substrate and adsorbate-adsorbate interactions of Na and K adlayers on Al(111)},
+  author = {Neugebauer, J\"org and Scheffler, Matthias},
+  journal = {Phys. Rev. B},
+  volume = {46},
+  issue = {24},
+  pages = {16067--16080},
+  year = {1992},
+  doi = {10.1103/PhysRevB.46.16067}
+}

src/DiagModule.f90

@@ -514,7 +514,6 @@ subroutine FindEvals(electrons)
     use functions_on_grid,           only: atomfns, &
          allocate_temp_fn_on_grid,    &
          free_temp_fn_on_grid
-    use density_module, ONLY: get_band_density
     use io_module, ONLY: write_eigenvalues, write_eigenvalues_format_ase
     use pao_format, ONLY: pao
     use ELPA_module, ONLY: flag_use_elpa, init_ELPA, end_ELPA

src/H_matrix_module.f90

@@ -630,7 +630,9 @@ subroutine get_h_on_atomfns(output_level, fixed_potential, &
     use block_module,                only: n_blocks, n_pts_in_block
     use primary_module,              only: domain
     use set_blipgrid_module,         only: naba_atoms_of_blocks
-    use density_module,              only: density_pcc, density_atom
+    use density_module,              only: density_pcc, density_atom, &
+         flag_surface_dipole_correction, get_surface_dipole, get_average_potential, &
+         flag_output_average_potential
     use GenComms,                    only: gsum, inode, ionode, cq_abort
     use energy,                      only: hartree_energy_total_rho,  &
                                            xc_energy,       &
@@ -799,6 +801,11 @@ subroutine get_h_on_atomfns(output_level, fixed_potential, &
     end if
     call gsum(delta_E_xc)
     delta_E_xc = delta_E_xc * grid_point_volume
+    if(flag_surface_dipole_correction) then
+       call get_surface_dipole(h_potential, rho_tot, size)
+    else if (flag_output_average_potential) then
+       call get_average_potential(h_potential, rho_tot, size)
+    end if
     !
     !
     ! Make total potential

src/control.f90

@@ -158,6 +158,7 @@ subroutine control_run(fixed_potential, vary_mu, total_energy)
     use store_matrix,         only: dump_pos_and_matrices
     use io_module,            only: write_extxyz
     use md_control,           only: flag_write_extxyz
+    use density_module,       only: flag_output_average_potential, write_average_potential
 
     implicit none
 
@@ -247,7 +248,7 @@ subroutine control_run(fixed_potential, vary_mu, total_energy)
 !****lat<$
     call stop_backtrace(t=backtrace_timer,who='control_run',echo=.true.)
 !****lat>$
-
+    if(flag_output_average_potential) call write_average_potential
     ! Added if, otherwise step numbering is broken on MD restart - zamaan
     if (.not. leqi(runtype, 'md')) call dump_pos_and_matrices
     return