oqp_resp_charges Subroutine

  subroutine oqp_resp_charges(infos)
    use precision, only: dp
    use io_constants, only: iw
    use oqp_tagarray_driver
    use basis_tools, only: basis_set
    use messages, only: show_message, with_abort
    use types, only: information
    use strings, only: Cstring, fstring
    use mathlib, only: traceprod_sym_packed, triangular_to_full
    use int1, only: electrostatic_potential
    use lebedev, only: lebedev_get_grid
    use elements, only: ELEMENTS_VDW_RADII

    implicit none

    character(len=*), parameter :: subroutine_name = "oqp_resp_charges"

    type(information), target, intent(inout) :: infos

    integer :: nbf, nbf2, ok
    logical :: urohf
    type(basis_set), pointer :: basis
    real(kind=dp), allocatable :: xyz(:,:), wt(:), pot(:)
    real(kind=dp), allocatable :: leb(:,:), lebw(:)
    real(kind=dp), allocatable :: vdwrad(:)
    real(kind=dp), allocatable :: chg(:)
    integer, allocatable :: neigh(:)
    real(kind=dp) :: rms
    real(kind=dp), allocatable :: den(:)
    real(kind=dp), allocatable :: q0(:), alpha(:)

    integer :: nat, npt, nptcur, nadd, nleb
    integer :: i, layer

    integer, parameter :: nlayers = 4
    real(kind=dp), parameter :: &
      layers(nlayers) = [1.4, 1.6, 1.8, 2.0]
    integer, parameter :: npt_layer(nlayers) = [132,152,192,350]
    integer, parameter :: typ_layer(nlayers) = [3,3,3,0]

    logical :: restr

    ! tagarray
    real(kind=dp), contiguous, pointer :: dmat_a(:), dmat_b(:)
    character(len=*), parameter :: tags_alpha(1) = (/ character(len=80) :: &
      OQP_DM_A /)
    character(len=*), parameter :: tags_beta(2) = (/ character(len=80) :: &
      OQP_DM_A, OQP_DM_B /)

    open (unit=IW, file=infos%log_filename, position="append")

    select case(infos%control%esp)
    case(0,1)
      restr = .false.
      write(iw,'(2/)')
      write(iw,'(4x,a)') '======================='
      write(iw,'(4x,a)') 'ESP charges calculation'
      write(iw,'(4x,a)') '======================='
    case (2)
      restr = .true.
      allocate(q0(nat), alpha(nat), source=0.0d0)
      alpha = infos%control%resp_constr
      select case(infos%control%resp_target)
      case(0)
        q0 = 0
      !case(1)
      ! set Mulliken charges
      case default
        error stop 'Unknown RESP charges target'
      end select
      write(iw,'(2/)')
      write(iw,'(4x,a)') '========================'
      write(iw,'(4x,a)') 'RESP charges calculation'
      write(iw,'(4x,a)') '========================'
    case default
      error stop 'Unknown type of ESP charges calculation'
    end select

    call flush(iw)

!   Load basis set
    basis => infos%basis
    basis%atoms => infos%atoms

!   Allocate memory
    nbf = basis%nbf
    nbf2 = nbf*(nbf+1)/2
    nat = ubound(infos%atoms%zn, 1)
    npt = nat * sum(npt_layer)
    urohf = infos%control%scftype == 2 .or. infos%control%scftype == 3

    allocate(xyz(npt,3), &
             wt(npt), &
             pot(npt), &
             leb(maxval(npt_layer),3), &
             lebw(maxval(npt_layer)), &
             vdwrad(nat), &
             chg(nat), &
             neigh(nat), &
             den(nbf2), &
             stat=ok)
    if (ok /= 0) call show_message('Cannot allocate memory', WITH_ABORT)

!   Set up the grid
    nptcur = 0  ! Current number of point in a grid

!   Loop over layers and add spherical grid points on each atom to the molecular grid
    do layer = 1, nlayers

!     Get the atomic radii on which to place new grid layer
      vdwrad = ELEMENTS_VDW_RADII(int(infos%atoms%zn))*layers(layer)

!     Get grid
      nleb = npt_layer(layer)
      leb = 0
      lebw = 0
      call lebedev_get_grid(nleb, leb, lebw, typ_layer(layer))

!     Add new grid layer for each atom, remove inner points
      do i = 1, nat
        call add_atom_grid( &
          x=xyz(nptcur+1:,1), &
          y=xyz(nptcur+1:,2), &
          z=xyz(nptcur+1:,3), &
          wts=wt(nptcur+1:), &
          nadd=nadd, &
          atpts=leb(:nleb,:), &
          atwts=lebw(:nleb), &
          atoms_xyz=infos%atoms%xyz, &
          atoms_rad=vdwrad, &
          cur_atom=i, &
          neighbours=neigh)
        nptcur = nptcur + nadd
      end do
    end do

!   Set all grid weights to be 1 for now
    wt = 1

    deallocate(leb, lebw, vdwrad, neigh)

!   Get density
    if (urohf) then
      call data_has_tags(infos%dat, tags_beta, module_name, subroutine_name, WITH_ABORT)
      call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a)
      call tagarray_get_data(infos%dat, OQP_DM_B, dmat_b)
      den = dmat_a + dmat_b
    else
      call data_has_tags(infos%dat, tags_alpha, module_name, subroutine_name, WITH_ABORT)
      call tagarray_get_data(infos%dat, OQP_DM_A, dmat_a)
      den = dmat_a
    end if

!   Compute electrostatic potential
    pot = 0
    call electrostatic_potential(basis, &
            x=xyz(:nptcur,1), &
            y=xyz(:nptcur,2), &
            z=xyz(:nptcur,3), &
            wt=wt(:nptcur), &
            d=den, &
            pot=pot)

    deallocate(den)

!   Add nuclei contribution to the potential
    call nuc_pot( &
            x=xyz(:nptcur,1), &
            y=xyz(:nptcur,2), &
            z=xyz(:nptcur,3), &
            w=wt(:nptcur), &
            at=infos%atoms%xyz, &
            q=infos%atoms%zn - infos%basis%ecp_zn_num, &
            pot=pot(:nptcur))

!   Fit ESP c8harges
    call chg_fit_mk( &
            x=xyz(:nptcur,1), &
            y=xyz(:nptcur,2), &
            z=xyz(:nptcur,3), &
            w=wt(:nptcur), &
            at=infos%atoms%xyz, &
            pot=pot(:nptcur), &
            chgtot=real(infos%mol_prop%charge, dp), &
            chg=chg, &
            resp=restr, &
            q0=q0, &
            alpha=alpha &
    )

!   Print ESP charges
    call print_charges(infos, chg)

!   Compute RMS error of the ponential induced by ESP charges
    call check_charges( &
            x=xyz(:nptcur,1), &
            y=xyz(:nptcur,2), &
            z=xyz(:nptcur,3), &
            w=wt(:nptcur), &
            at=infos%atoms%xyz, &
            chg=chg, &
            pot=pot(:nptcur), &
            rms=rms)
    write(*,'(x,a,es20.6)') 'rms.err.=', rms

    close(iw)

  end subroutine oqp_resp_charges