new_d4_model Subroutine

subroutine new_d4_model(error, d4, mol, ga, gc, wf, qmod)
   !DEC$ ATTRIBUTES DLLEXPORT :: new_d4_model

   !> Instance of the dispersion model
   type(d4_model), intent(out) :: d4

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

   !> Error handling
   type(error_type), allocatable, intent(out) :: error

   !> Charge scaling height
   real(wp), intent(in), optional :: ga

   !> Charge scaling steepness
   real(wp), intent(in), optional :: gc

   !> Weighting factor for coordination number interpolation
   real(wp), intent(in), optional :: wf

   !> Charge model selection
   integer, intent(in), optional :: qmod

   integer :: isp, izp, iref, jsp, jzp, jref
   integer :: mref, tmp_qmod
   real(wp) :: aiw(23), c6
   real(wp), parameter :: thopi = 3.0_wp/pi

   ! check for unsupported elements (104 (Rf) - 111 (Rg))
   do isp = 1, mol%nid
      if (mol%num(isp) > 103 .and. mol%num(isp) < 112) then
         call fatal_error(error, "Structure contains unsupported element '"//trim(mol%sym(isp))//"'")
         return
      end if
   end do

   d4%ncoup = 1

   if (present(ga)) then
      d4%ga = ga
   else
      d4%ga = ga_default
   end if

   if (present(gc)) then
      d4%gc = gc
   else
      d4%gc = gc_default
   end if

   if (present(wf)) then
      d4%wf = wf
   else
      d4%wf = wf_default
   end if

   allocate(d4%rcov(mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      d4%rcov(isp) = get_covalent_rad(izp)
   end do

   allocate(d4%en(mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      d4%en(isp) = get_electronegativity(izp)
   end do

   allocate(d4%zeff(mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      d4%zeff(isp) = get_effective_charge(izp)
   end do

   allocate(d4%eta(mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      d4%eta(isp) = get_hardness(izp)
   end do

   allocate(d4%r4r2(mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      d4%r4r2(isp) = get_r4r2_val(izp)
   end do

   allocate(d4%ref(mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      d4%ref(isp) = get_nref(izp)
   end do

   mref = maxval(d4%ref)
   allocate(d4%cn(mref, mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      call set_refcn(d4%cn(:, isp), izp)
   end do

   if (present(qmod)) then
      tmp_qmod = qmod
   else
      tmp_qmod = d4_qmod%eeq
   end if

   allocate(d4%q(mref, mol%nid))
   allocate(d4%aiw(23, mref, mol%nid))
   select case(tmp_qmod)
   case default
      call fatal_error(error, "Unsupported option for charge model.")
      return
   case(d4_qmod%eeq)
      do isp = 1, mol%nid
         izp = mol%num(isp)
         call set_refq_eeq(d4%q(:, isp), izp)
         call set_refalpha_eeq(d4%aiw(:, :, isp), d4%ga, d4%gc, izp)
      end do
      ! Setup EEQ model
      call new_eeq2019_model(mol, d4%mchrg, error)
      if(allocated(error)) return
   case(d4_qmod%eeqbc)
      do isp = 1, mol%nid
         izp = mol%num(isp)
         call set_refq_eeqbc(d4%q(:, isp), izp)
         call set_refalpha_eeqbc(d4%aiw(:, :, isp), d4%ga, d4%gc, izp)
      end do
      ! Setup EEQBC model
      call new_eeqbc2025_model(mol, d4%mchrg, error)  
      if(allocated(error)) return
   case(d4_qmod%gfn2)
      do isp = 1, mol%nid
         izp = mol%num(isp)
         call set_refq_gfn2(d4%q(:, isp), izp)
         call set_refalpha_gfn2(d4%aiw(:, :, isp), d4%ga, d4%gc, izp)
      end do
   end select

   allocate(d4%ngw(mref, mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      call set_refgw(d4%ngw(:, isp), izp)
   end do

   allocate(d4%c6(mref, mref, mol%nid, mol%nid))
   do isp = 1, mol%nid
      izp = mol%num(isp)
      do jsp = 1, isp
         jzp = mol%num(jsp)
         do iref = 1, d4%ref(isp)
            do jref = 1, d4%ref(jsp)
               aiw(:) = d4%aiw(:, iref, isp) * d4%aiw(:, jref, jsp)
               c6 = thopi * trapzd(aiw)
               d4%c6(jref, iref, jsp, isp) = c6
               d4%c6(iref, jref, isp, jsp) = c6
            end do
         end do
      end do
   end do

end subroutine new_d4_model