get_charges Subroutine

public subroutine get_charges(mol, qvec, dqdr, dqdL)

Obtain charges from electronegativity equilibration model

Arguments

Type IntentOptional Attributes Name
type(structure_type), intent(in) :: mol

Molecular structure data

real(kind=wp), intent(out), contiguous :: qvec(:)

Atomic partial charges

real(kind=wp), intent(out), optional, contiguous :: dqdr(:,:,:)

Derivative of the partial charges w.r.t. the Cartesian coordinates

real(kind=wp), intent(out), optional, contiguous :: dqdL(:,:,:)

Derivative of the partial charges w.r.t. strain deformations


Source Code

subroutine get_charges(mol, qvec, dqdr, dqdL)
   !DEC$ ATTRIBUTES DLLEXPORT :: get_charges

   !> Molecular structure data
   type(structure_type), intent(in) :: mol

   !> Atomic partial charges
   real(wp), intent(out), contiguous :: qvec(:)

   !> Derivative of the partial charges w.r.t. the Cartesian coordinates
   real(wp), intent(out), contiguous, optional :: dqdr(:, :, :)

   !> Derivative of the partial charges w.r.t. strain deformations
   real(wp), intent(out), contiguous, optional :: dqdL(:, :, :)

   logical :: grad
   type(mchrg_model_type) :: model
   real(wp), parameter :: cn_max = 8.0_wp, cutoff = 25.0_wp
   real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :)
   real(wp), allocatable :: rcov(:), trans(:, :)

   grad = present(dqdr) .and. present(dqdL)

   call new_eeq2019_model(mol, model)
   call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans)

   allocate(cn(mol%nat))
   if (grad) then
      allocate(dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat))
   end if

   rcov = get_covalent_rad(mol%num)
   call get_coordination_number(mol, trans, cutoff, rcov, cn, dcndr, dcndL, cut=cn_max)

   call model%solve(mol, cn, dcndr, dcndL, qvec=qvec, dqdr=dqdr, dqdL=dqdL)

end subroutine get_charges