ddCOSMO Surface‐Area (Gradient) Is Noisy/Discontinuous

[lukaswittmann]

Description:
When scanning the O–H1 distance in H1-OH from 0.8 Å → 1.9 Å (as in J. Comp. Chem. 20, 428–446 (1999), DOI:10.1002/jcc.70099, Figure 4), the surface area A (and the numerical gradient ∂A/∂R) shows larger oscillations and sharp kinks (and possibly discontinuities at $n\approx80$ and $n\approx330$). In contrast, an ISWIG‐style area (using the same quadrature grid) is smooth over the same range.

• Total surface area:
  
• Numerical gradient of the surface area (or here more precisely for simplicity A/nsteps; nsteps are all equidistant):
  

Steps to Reproduce:

I used ddinit and ddrun and used the results for the calculation of the surface area via $a_i = w_i u_i R_I$ (by using ui_cav).
This was done using the bleeding egde version of ddx.

(Implementational details.)

      ! allocate multipoles
      allocate(self%multipoles(1, mol%nat))
      self%multipoles(1, :) = qat(:) / sqrt(4.0_wp * pi)

      ! Initialize ddX
      call ddinit(  &
         1,         &
         mol%nat,   &
         mol%xyz,   &
         rvdw,      &
         1000.0_wp, & ! This does not matter for cpcm
         self%ddx,  &
         error,     & 
         eta=0.1_wp,& 
         ngrid=110)
      call check_error(error)

      ! Allocate state
      call allocate_state(   &
         self%ddx%params,    &
         self%ddx%constants, &
         self%state,         &
         error)
      call check_error(error)

      call multipole_electrostatics(self%ddx%params, self%ddx%constants, &
         & self%ddx%workspace, self%multipoles, 0, self%electrostatics, error)
      call multipole_psi(self%ddx%params, self%multipoles, 0, self%state%psi)
      call check_error(error)

...

      call ddrun(             &
         self%ddx,            &
         self%state,          &  
         self%electrostatics, &
         self%state%psi,      &
         1E-10_wp,            &
         ddx_energy,          &
         error)
      call check_error(error)

!> Surface area per grid point of a dimension (ncav)
allocate(self%a(ddx_model%ddx%constants%ncav))
do k = 1, ddx_model%ddx%constants%ncav
   gidx = self%cav2ngrid(k)
   iat  = self%owner(k)
   self%a(k) = ddx_model%ddx%constants%wgrid(gidx) &
      * ddx_model%ddx%constants%ui_cav(k) &
      * (ddx_model%ddx%params%rsph(iat)**2)
   ! write(*, '(3i5, 4f10.4)') k, gidx, iat, ddx_model%ddx%constants%wgrid(gidx), &
   !    ddx_model%ddx%constants%ui_cav(k), self%a(k)
end do
self%totalArea = sum(self%a)

(Calculation of cav2grid and surface-point owner.)

The mapping from ncav to ngrid was done using cav2grid, similarly with mapping the surface point $i$ to the owner atom $I$ with owner.

!> Build cavity owner and grid point mapping
allocate(self%cav2ngrid(self%ncav))
allocate(self%owner(self%ncav))
! Walk the CSR structure:
!  for sphere i, the cavity points are in k = icav_ia(i) : icav_ia(i+1)-1
do i = 1, self%nsph
   do k = ddx_model%ddx%constants%icav_ia(i), ddx_model%ddx%constants%icav_ia(i+1) - 1
      self%cav2ngrid(k) = ddx_model%ddx%constants%icav_ja(k)
      self%owner(k) = i
   end do
end do

Is there something I am overlooking in the implementation of the surface area (i.e. needed additional terms)?

I also realized that the area is very dependent on the choice of $\eta$, is there a way to circumvent this?

Thank you for your help!