From cdd0a85600bbf0a039ca6cf02ae31ea0b776d38a Mon Sep 17 00:00:00 2001 From: Marvin Friede <51965259+marvinfriede@users.noreply.github.com> Date: Wed, 9 Aug 2023 20:05:11 +0200 Subject: [PATCH] Print atomic dispersion energies (#46) --- app/driver.f90 | 14 +++++++++---- src/dftd3.f90 | 3 ++- src/dftd3/disp.f90 | 49 +++++++++++++++++++++++++++++++++++++++----- src/dftd3/output.f90 | 30 +++++++++++++++++++++++++++ 4 files changed, 86 insertions(+), 10 deletions(-) diff --git a/app/driver.f90 b/app/driver.f90 index 3e19f036..23d5b11f 100644 --- a/app/driver.f90 +++ b/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) diff --git a/src/dftd3.f90 b/src/dftd3.f90 index 3f375cb5..bd3436e0 100644 --- a/src/dftd3.f90 +++ b/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 diff --git a/src/dftd3/disp.f90 b/src/dftd3/disp.f90 index 458d6d5f..90c67c0d 100644 --- a/src/dftd3/disp.f90 +++ b/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 diff --git a/src/dftd3/output.f90 b/src/dftd3/output.f90 index 10222899..b13513d3 100644 --- a/src/dftd3/output.f90 +++ b/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