Add "--tmcosmo" mode for writing .cosmo files with TM convention. (#864)

* Add "--tmcosmo" mode for writing .cosmo files with TM convention.
* Add ddvolume module for Volume output (and handling improvement)

---------

Signed-off-by: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
Co-authored-by: albert <92109627+Albkat@users.noreply.github.com>

src/prog/main.F90

@@ -1572,11 +1572,11 @@ subroutine parseArguments(env, args, inputFile, paramFile, lgrad, &
             call env%error("No solvent name provided for ALPB", source)
          end if
 
-      case('--cosmo')
+      case('--cosmo','--tmcosmo')
          call args%nextArg(sec)
          if (allocated(sec)) then
             call set_gbsa(env, 'solvent', sec)
-            call set_gbsa(env, 'cosmo', 'true')
+            call set_gbsa(env, flag(3:), 'true')
             call args%nextArg(sec)
             if (allocated(sec)) then
                if (sec == 'reference') then

src/set_module.f90

@@ -2021,6 +2021,7 @@ subroutine set_gbsa(env,key,val)
    logical,save :: set6 = .true.
    logical,save :: set7 = .true.
    logical,save :: set8 = .true.
+   logical,save :: set9 = .true.
    select case(key)
    case default ! do nothing
       call env%warning("the key '"//key//"' is not recognized by gbsa",source)
@@ -2074,8 +2075,15 @@ subroutine set_gbsa(env,key,val)
    case('cosmo')
       if (getValue(env,val,ldum).and.set7) set%solvInput%cosmo = ldum
       set7 = .false.
+   case('tmcosmo')
+      if (getValue(env,val,ldum).and.set8) then
+         set%solvInput%cosmo = ldum
+         set%solvInput%tmcosmo = .true.
+      end if
+      set8 = .false.
    case('cpcmx')
-      if (set8) set%solvInput%cpxsolvent = val
+      if (set9) set%solvInput%cpxsolvent = val
+      set9 = .false.
    end select
 end subroutine set_gbsa
 

src/solv/CMakeLists.txt

@@ -31,6 +31,7 @@ list(APPEND srcs
   "${dir}/sasa.f90"
   "${dir}/state.f90"
   "${dir}/cpx.F90"
+  "${dir}/ddvolume.f90"
 )
 
 set(srcs ${srcs} PARENT_SCOPE)

src/solv/cosmo.f90

@@ -95,6 +95,12 @@ module xtb_solv_cosmo
       real(wp), allocatable :: dsdr(:, :)
       real(wp), allocatable :: dsdrt(:, :, :)
 
+      !> TM convention? 
+      logical :: tmcosmo = .false.
+      !> Write volume in COSMO file? Volume resolution in Angstrom.
+      logical :: wrvolume = .false.
+      real(wp) :: volumeres = 0.1_wp
+
    contains
 
       !> Update coordinates and internal state
@@ -827,6 +833,7 @@ end subroutine update_nnlist_sasa
 
 !> Write a COSMO file output
 subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy)
+   use xtb_solv_ddvolume, only: ddvolume
 
    !> COSMO container
    class(TCosmo), intent(in) :: self
@@ -851,7 +858,7 @@ subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy)
 
    integer :: ii, ig, iat
    real(wp) :: dielEnergy, keps
-   real(wp), allocatable :: phi(:), zeta(:), area(:)
+   real(wp), allocatable :: phi(:), zeta(:), area(:), volume(:)
 
    allocate(phi(self%ddCosmo%ncav), zeta(self%ddCosmo%ncav), area(self%ddCosmo%ncav))
    ! Reset potential on the cavity, note that the potential is expected in e/Å
@@ -870,6 +877,10 @@ subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy)
       end do
    end do
 
+   !! Switch convention for TM mode
+   if (self%tmcosmo) zeta=-zeta
+   if (self%wrvolume) call ddvolume(xyz,self%rvdw,(self%volumeres)/autoaa,volume)
+
 
     ! Dielectric energy is the energy on the dielectric continuum
    keps = 0.5_wp * (1.0_wp - 1.0_wp/self%dielectricConst)
@@ -886,9 +897,16 @@ subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy)
 
    write(unit, '(a)') &
       & "$cosmo_data"
-   write(unit, '(2x, a:, "=", g0)') &
-      & "fepsi", keps, &
-      & "area", sum(area)
+   if (self%wrvolume) then
+      write(unit, '(2x, a:, "=", g0)') &
+         & "fepsi", keps, &
+         & "area", sum(area), &
+         & "volume", sum(volume)
+   else
+      write(unit, '(2x, a:, "=", g0)') &
+         & "fepsi", keps, &
+         & "area", sum(area)
+   end if
 
    write(unit, '(a)') &
       & "$coord_rad", &
@@ -940,6 +958,7 @@ subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy)
       & "#", &
       & "#"
 
+
    ii = 0
    do iat = 1, self%ddCosmo%nat
       do ig = 1, self%ddCosmo%ngrid

src/solv/ddvolume.f90

@@ -0,0 +1,200 @@
+! This file is part of xtb.
+! SPDX-Identifier: LGPL-3.0-or-later
+!
+! xtb is free software: you can redistribute it and/or modify it under
+! the terms of the GNU Lesser General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! xtb is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU Lesser General Public License for more details.
+!
+! You should have received a copy of the GNU Lesser General Public License
+! along with xtb.  If not, see <https://www.gnu.org/licenses/>.
+!
+! MS, 09/2023: Initial implementation
+
+
+!> domain-decomposed volume calculation for molecule
+module xtb_solv_ddvolume
+    use xtb_mctc_accuracy, only: wp
+    implicit none
+    private
+
+    public :: ddvolume
+
+    type :: sphere_type
+        real(wp) :: x,y,z,r
+        integer, allocatable :: neighbours(:)
+        integer :: num_neighbours
+        real(wp) :: volume
+    contains
+        procedure :: init => sphere_init
+    end type sphere_type
+
+    interface distance2
+        module procedure :: distance2_spheres
+        module procedure :: distance2_spherepoint
+    end interface
+
+contains 
+
+    !> Calculates the volumes of an array of xyz-coordinates with a given radii and resolution based on a custom dd algorithm
+    subroutine ddvolume(xyz,r,resolution,volume)
+        !> xyz-coordinates (3,nat)
+        real(wp), dimension(:,:), intent(in) :: xyz
+        !> Radius (nat)
+        real(wp), dimension(:), intent(in) :: r
+        !> Resolution
+        real(wp), intent(in) :: resolution
+        !> Volume (nat)
+        real(wp), dimension(:), allocatable, intent(out) :: volume
+
+        !> Transfer input to sphere type
+        type(sphere_type), allocatable :: spheres(:)
+        integer :: i,s,x,y,z
+        real(wp) :: x_min,x_max,y_min,y_max,z_min,z_max
+
+        allocate(spheres(size(xyz,2)))
+        allocate(volume(size(xyz,2)))
+        volume = 0.0_wp
+        !! Initialize spheres
+        do i=1,size(r)
+            call spheres(i)%init(xyz(1,i),xyz(2,i),xyz(3,i),r(i),size(xyz,2))
+        end do
+
+        call find_neighbours(spheres)
+        !$OMP PARALLEL DO DEFAULT(NONE) SHARED(spheres, volume,resolution) PRIVATE(s,x,y,z,x_min,x_max,y_min,y_max,z_min,z_max)
+        do s=1,size(spheres)
+            x_min=(spheres(s)%x-spheres(s)%r)
+            x_max=(spheres(s)%x+spheres(s)%r)
+            y_min=(spheres(s)%y-spheres(s)%r)
+            y_max=(spheres(s)%y+spheres(s)%r)
+            z_min=(spheres(s)%z-spheres(s)%r)
+            z_max=(spheres(s)%z+spheres(s)%r)
+            do x = floor(x_min/resolution), ceiling(x_max/resolution)
+                do y = floor(y_min/resolution), ceiling(y_max/resolution)
+                    do z = floor(z_min/resolution), ceiling(z_max/resolution)
+                        if (is_inside(x*resolution,y*resolution,z*resolution,spheres,s)) then
+                            spheres(s)%volume = spheres(s)%volume + resolution**3
+                        end if
+                    end do
+                end do
+            end do
+            volume(s) = spheres(s)%volume
+      end do
+      !$OMP END PARALLEL DO
+    end subroutine ddvolume
+
+    !> Initializes a sphere
+    subroutine sphere_init(self,x,y,z,r,count)
+        class(sphere_type), intent(inout) :: self
+        !> Coordinates and radii
+        real(wp), intent(in) :: x,y,z,r
+        !> Number of maximum possible neighbours (amount of spheres)
+        integer, intent(in) :: count
+
+        self%x = x
+        self%y = y
+        self%z = z
+        self%r = r
+        self%num_neighbours = 0
+        self%volume = 0.0_wp
+        allocate(self%neighbours(count))
+        self%neighbours = -1
+
+    end subroutine sphere_init
+
+    !> Set up neighbour list for a given array of spheres
+    subroutine find_neighbours(spheres)
+        !> Array of Sphere Type
+        type(sphere_type), intent(inout) :: spheres(:)
+
+        integer :: i,j,k
+        real(wp) :: dist2
+
+        do i=1,size(spheres)-1
+            do j=i+1,size(spheres)
+                dist2 = distance2(spheres(i),spheres(j))
+                if (dist2 < (spheres(i)%r+spheres(j)%r)**2) then
+                    do k=1,size(spheres(i)%neighbours)
+                        if (spheres(i)%neighbours(k) == -1) then
+                            spheres(i)%neighbours(k) = j
+                            spheres(i)%num_neighbours = spheres(i)%num_neighbours + 1
+                            exit
+                        end if
+                    end do
+                    do k=1,size(spheres(j)%neighbours)
+                        if (spheres(j)%neighbours(k) == -1) then
+                            spheres(j)%neighbours(k) = i
+                            spheres(j)%num_neighbours = spheres(j)%num_neighbours + 1
+                            exit
+                        end if
+                    end do
+                end if
+            end do
+        end do
+    end subroutine
+
+    !> Quadratic distance between two spheres
+    function distance2_spheres(sphere1,sphere2) result(distance2)
+        type(sphere_type), intent(in) :: sphere1,sphere2
+        real(wp) :: distance2
+
+        distance2 = (sphere1%x-sphere2%x)**2 + (sphere1%y-sphere2%y)**2 + (sphere1%z-sphere2%z)**2
+    end function distance2_spheres
+
+    !> Quadratic distance between a sphere and a point
+    function distance2_spherepoint(sphere1,x,y,z) result(distance2)
+        type(sphere_type), intent(in) :: sphere1
+        real(wp), intent(in) :: x,y,z
+        real(wp) :: distance2
+
+        distance2 = (sphere1%x-x)**2 + (sphere1%y-y)**2 + (sphere1%z-z)**2
+    end function distance2_spherepoint
+
+    !> Checks if a point is inside a sphere or not
+    function is_inside(x,y,z,spheres,this) result(inside)
+        implicit none
+        !> Array of Sphere Type (all spheres)
+        type(sphere_type), intent(in) :: spheres(:)
+        !> Index of the sphere to check
+        integer, intent(in) :: this
+        !> Coordinates of the point to check
+        real(wp), intent(in) :: x,y,z
+        !> Is it inside?
+        logical :: inside
+
+
+        integer :: i,n
+        real(wp) :: dist2, distn2
+
+
+        !! Check if point is even inside the sphere
+        dist2 = distance2(spheres(this),x,y,z)
+        if (dist2 <= spheres(this)%r**2) then
+            inside = .true.
+        else
+            inside = .false.
+            return
+        end if
+
+        !! Check if point is nearer to one of the neighbours than to the sphere itself
+        do i=1,spheres(this)%num_neighbours
+            n = spheres(this)%neighbours(i)
+            distn2 = distance2(spheres(n),x,y,z)
+            if (distn2 < dist2 .and. distn2 <= spheres(n)%r**2) then
+                inside = .false.
+                return
+            end if
+        end do
+      end function is_inside
+
+
+end module xtb_solv_ddvolume
+
+
+
+

src/solv/input.F90

@@ -42,6 +42,9 @@ module xtb_solv_input
       !> Use COSMO solvation model instead of Born based one
       logical :: cosmo = .false.
 
+      !> Use TM convention for COSMO solvation model
+      logical :: tmcosmo = .false.
+
       !> Additional CPCM-X mode
       character(len=:),allocatable :: cpxsolvent
 

src/solv/meson.build

@@ -29,4 +29,5 @@ srcs += files(
   'model.f90',
   'sasa.f90',
   'state.f90',
+  'ddvolume.f90',
 )

src/solv/model.f90

@@ -693,6 +693,10 @@ subroutine newSolvationModel(self, env, model, num)
       allocate(cosmo)
       call init_(cosmo, env, num, self%dielectricConst, self%nAng, self%bornScale, &
          & self%vdwRad, self%surfaceTension, self%probeRad, srcut)
+      if (self%tmcosmo) then
+         cosmo%tmcosmo=.true.
+         cosmo%wrvolume=.true.
+      end if
       call move_alloc(cosmo, model)
    else
       allocate(born)

src/xhelp.f90

@@ -165,6 +165,9 @@ subroutine help(iunit)
    "    Additionally, the dielectric constant can be set manually or an ideal conductor", &
    "    can be chosen by setting epsilon to infinity.",&
    "",&
+   "--tmcosmo SOLVENT/EPSILON",&
+   "    same as --cosmo, but uses TM convention for writing the .cosmo files.",&
+   "",&
    "--cpcmx SOLVENT",&
    "    extended conduction-like polarizable continuum solvation model (CPCM-X),",&
    "    available solvents are all solvents included in the Minnesota Solvation Database.",&