add_solvation_to_calc Subroutine

   subroutine add_solvation_to_calc(calc, mol, solvent, solvation_model, method, use_cds, use_shift, &
                                    dielectric, cpcm_nang, cpcm_rscale, error)
      !! Add implicit solvation model to XTB calculator
      !!
      !! Adds ALPB, GBSA, or CPCM solvation. For ALPB/GBSA, optionally adds CDS and shift corrections.
      !! CPCM does not support CDS or shift corrections.
      type(xtb_calculator), intent(inout) :: calc
      type(structure_type), intent(in) :: mol
      character(len=*), intent(in) :: solvent           !! Solvent name (can be empty if dielectric > 0)
      character(len=*), intent(in) :: solvation_model   !! "alpb", "gbsa", or "cpcm"
      character(len=*), intent(in) :: method            !! "gfn1" or "gfn2"
      logical, intent(in) :: use_cds, use_shift
      real(wp), intent(in) :: dielectric                !! Direct dielectric constant (-1 = use solvent lookup)
      integer, intent(in) :: cpcm_nang                  !! Angular grid points for CPCM
      real(wp), intent(in) :: cpcm_rscale               !! Radii scaling for CPCM
      type(error_type), allocatable, intent(out) :: error

      type(solvation_input) :: solv_input
      class(container_type), allocatable :: cont
      class(solvation_type), allocatable :: solv
      logical :: use_alpb
      real(wp) :: eps

      ! Handle CPCM model separately
      if (trim(solvation_model) == 'cpcm') then
         ! CPCM does not support CDS or shift - silently skip them
         ! (use_cds and use_shift are ignored for CPCM)

         ! Get dielectric constant (direct value or lookup from solvent name)
         if (dielectric > 0.0_wp) then
            eps = dielectric
         else if (len_trim(solvent) > 0) then
            eps = get_solvent_dielectric(solvent)
            if (eps < 0.0_wp) then
               allocate (error)
               error%message = "Unknown solvent for CPCM dielectric lookup: "//trim(solvent)
               return
            end if
         else
            allocate (error)
            error%message = "CPCM requires either solvent name or dielectric constant"
            return
         end if

         ! Create CPCM solvation
         allocate (solv_input%cpcm)
         solv_input%cpcm%dielectric_const = eps
         solv_input%cpcm%nang = cpcm_nang
         solv_input%cpcm%rscale = cpcm_rscale

         call new_solvation(solv, mol, solv_input, error)
         if (allocated(error)) return
         call move_alloc(solv, cont)
         call calc%push_back(cont)

         return
      end if

      ! For ALPB/GBSA, we need a solvent name
      if (len_trim(solvent) == 0) then
         allocate (error)
         error%message = "ALPB/GBSA solvation requires a solvent name"
         return
      end if

      ! Determine if using ALPB or GBSA (GBSA = ALPB with alpb flag false)
      use_alpb = .true.
      if (trim(solvation_model) == 'gbsa') then
         use_alpb = .false.
      end if

      ! 1. Add ALPB/GBSA (polar electrostatic solvation)
      allocate (solv_input%alpb)
      solv_input%alpb%solvent = solvent
      solv_input%alpb%alpb = use_alpb

      call new_solvation(solv, mol, solv_input, error, method)
      if (allocated(error)) return
      call move_alloc(solv, cont)
      call calc%push_back(cont)

      deallocate (solv_input%alpb)

      ! 2. Add CDS (non-polar: cavity, dispersion, surface) if requested
      if (use_cds) then
         allocate (solv_input%cds)
         solv_input%cds%solvent = solvent

         call new_solvation_cds(solv, mol, solv_input, error, method)
         if (allocated(error)) return
         call move_alloc(solv, cont)
         call calc%push_back(cont)

         deallocate (solv_input%cds)
      end if

      ! 3. Add shift (solution state correction) if requested
      if (use_shift) then
         allocate (solv_input%shift)
         solv_input%shift%solvent = solvent

         call new_solvation_shift(solv, solv_input, error, method)
         if (allocated(error)) return
         call move_alloc(solv, cont)
         call calc%push_back(cont)
      end if
   end subroutine add_solvation_to_calc