Print atomic dispersion energies (#46)

app/driver.f90

@@ -51,9 +51,10 @@ subroutine run_driver(config, error)
    class(damping_param), allocatable :: param
    type(d3_param) :: inp
    type(d3_model) :: d3
-   real(wp), allocatable :: energy, gradient(:, :), sigma(:, :)
+   real(wp), allocatable :: energies(:), gradient(:, :), sigma(:, :)
    real(wp), allocatable :: pair_disp2(:, :), pair_disp3(:, :)
    real(wp), allocatable :: s9
+   real(wp) :: energy
    integer :: stat, unit
    logical :: exist
 
@@ -163,7 +164,7 @@ subroutine run_driver(config, error)
    end if
 
    if (allocated(param)) then
-      energy = 0.0_wp
+      allocate(energies(mol%nat))
       if (config%grad) then
          allocate(gradient(3, mol%nat), sigma(3, 3))
       end if
@@ -176,14 +177,19 @@ subroutine run_driver(config, error)
    end if
 
    if (allocated(param)) then
-      call get_dispersion(mol, d3, param, realspace_cutoff(), energy, gradient, &
-         & sigma)
+      call get_dispersion(mol, d3, param, realspace_cutoff(), energies, &
+         & gradient, sigma)
+      energy = sum(energies)
+
       if (config%pair_resolved) then
          allocate(pair_disp2(mol%nat, mol%nat), pair_disp3(mol%nat, mol%nat))
          call get_pairwise_dispersion(mol, d3, param, realspace_cutoff(), pair_disp2, &
             & pair_disp3)
       end if
       if (config%verbosity > 0) then
+         if (config%verbosity > 2) then
+            call ascii_energy_atom(output_unit, mol, energies)
+         end if
          call ascii_results(output_unit, mol, energy, gradient, sigma)
          if (config%pair_resolved) then
             call ascii_pairwise(output_unit, mol, pair_disp2, pair_disp3)

src/dftd3.f90

@@ -31,7 +31,8 @@ module dftd3
    implicit none
    private
 
-   public :: get_dispersion, get_pairwise_dispersion, get_coordination_number
+   public :: get_dispersion, get_pairwise_dispersion
+   public :: get_coordination_number
    public :: realspace_cutoff, get_lattice_points
    public :: damping_param, d3_param
    public :: get_rational_damping, get_zero_damping

src/dftd3/disp.f90

@@ -30,10 +30,17 @@ module dftd3_disp
    public :: get_dispersion, get_pairwise_dispersion
 
 
+   !> Calculate dispersion energy
+   interface get_dispersion
+      module procedure :: get_dispersion_atomic
+      module procedure :: get_dispersion_scalar
+   end interface get_dispersion
+
 contains
 
 
-subroutine get_dispersion(mol, disp, param, cutoff, energy, gradient, sigma)
+!> Calculate atom-resolved dispersion energies.
+subroutine get_dispersion_atomic(mol, disp, param, cutoff, energies, gradient, sigma)
 
    !> Molecular structure data
    class(structure_type), intent(in) :: mol
@@ -48,7 +55,7 @@ subroutine get_dispersion(mol, disp, param, cutoff, energy, gradient, sigma)
    type(realspace_cutoff), intent(in) :: cutoff
 
    !> Dispersion energy
-   real(wp), intent(out) :: energy
+   real(wp), intent(out) :: energies(:)
 
    !> Dispersion gradient
    real(wp), intent(out), contiguous, optional :: gradient(:, :)
@@ -61,7 +68,7 @@ subroutine get_dispersion(mol, disp, param, cutoff, energy, gradient, sigma)
    real(wp), allocatable :: cn(:)
    real(wp), allocatable :: gwvec(:, :), gwdcn(:, :)
    real(wp), allocatable :: c6(:, :), dc6dcn(:, :)
-   real(wp), allocatable :: dEdcn(:), energies(:)
+   real(wp), allocatable :: dEdcn(:)
    real(wp), allocatable :: lattr(:, :)
 
    mref = maxval(disp%ref)
@@ -79,7 +86,6 @@ subroutine get_dispersion(mol, disp, param, cutoff, energy, gradient, sigma)
    if (grad) allocate(dc6dcn(mol%nat, mol%nat))
    call disp%get_atomic_c6(mol, gwvec, gwdcn, c6, dc6dcn)
 
-   allocate(energies(mol%nat))
    energies(:) = 0.0_wp
    if (grad) then
       allocate(dEdcn(mol%nat))
@@ -98,9 +104,42 @@ subroutine get_dispersion(mol, disp, param, cutoff, energy, gradient, sigma)
          & gradient, sigma)
    end if
 
+end subroutine get_dispersion_atomic
+
+
+!> Calculate scalar dispersion energy.
+subroutine get_dispersion_scalar(mol, disp, param, cutoff, energy, gradient, sigma)
+
+   !> Molecular structure data
+   class(structure_type), intent(in) :: mol
+
+   !> Dispersion model
+   class(d3_model), intent(in) :: disp
+
+   !> Damping parameters
+   class(damping_param), intent(in) :: param
+
+   !> Realspace cutoffs
+   type(realspace_cutoff), intent(in) :: cutoff
+
+   !> Dispersion energy
+   real(wp), intent(out) :: energy
+
+   !> Dispersion gradient
+   real(wp), intent(out), contiguous, optional :: gradient(:, :)
+
+   !> Dispersion virial
+   real(wp), intent(out), contiguous, optional :: sigma(:, :)
+
+   real(wp), allocatable :: energies(:)
+
+   allocate(energies(mol%nat))
+
+   call get_dispersion_atomic(mol, disp, param, cutoff, energies, gradient, sigma)
+
    energy = sum(energies)
 
-end subroutine get_dispersion
+end subroutine get_dispersion_scalar
 
 
 !> Wrapper to handle the evaluation of pairwise representation of the dispersion energy

src/dftd3/output.f90

@@ -30,6 +30,7 @@ module dftd3_output
    private
 
    public :: ascii_atomic_radii, ascii_atomic_references, ascii_system_properties
+   public :: ascii_energy_atom
    public :: ascii_results, ascii_damping_param, ascii_pairwise
    public :: turbomole_gradient, turbomole_gradlatt
    public :: json_results, tagged_result
@@ -139,6 +140,35 @@ subroutine ascii_system_properties(unit, mol, disp, cn, c6)
 end subroutine ascii_system_properties
 
 
+!> Print atom-resolved dispersion energies
+subroutine ascii_energy_atom(unit, mol, energies)
+
+   !> Unit for output
+   integer, intent(in) :: unit
+
+   !> Molecular structure data
+   class(structure_type), intent(in) :: mol
+
+   !> Atom-resolved dispersion energies
+   real(wp), allocatable, intent(in) :: energies(:)
+
+   integer :: iat, isp
+
+   write(unit, '(a,":")') "Atom-resolved dispersion energies"
+   write(unit, '(50("-"))')
+   write(unit, '(a6,1x,a4,1x,4x,a15,1x,a15)') "#", "Z", "[Hartree]", "[kcal/mol]"
+   write(unit, '(50("-"))')
+   do iat = 1, mol%nat
+      isp = mol%id(iat)
+      write(unit, '(i6,1x,i4,1x,a4,e15.8,1x,f15.8)') &
+         & iat, mol%num(isp), mol%sym(isp), energies(iat), energies(iat)*autokcal
+   end do
+   write(unit, '(50("-"))')
+   write(unit, '(a)')
+
+end subroutine ascii_energy_atom
+
+
 subroutine ascii_results(unit, mol, energy, gradient, sigma)
 
    !> Unit for output