From d5dd8b9bcdc1629a40a415b09cc2ddc84e1405be Mon Sep 17 00:00:00 2001
From: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
Date: Wed, 13 Sep 2023 12:13:52 +0200
Subject: [PATCH 1/6] Add "--tmcosmo" mode for writing .cosmo files with TM
 convention.

Signed-off-by: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
---
 src/prog/main.F90  | 19 +++++++++++++++++++
 src/scf_module.F90 |  2 +-
 src/set_module.f90 | 10 +++++++++-
 src/solv/cosmo.f90 | 18 +++++++++++++++++-
 src/solv/input.F90 |  3 +++
 src/solv/model.f90 |  1 +
 src/xhelp.f90      |  3 +++
 7 files changed, 53 insertions(+), 3 deletions(-)

diff --git a/src/prog/main.F90 b/src/prog/main.F90
index cba5b8e91..aef4b3c2e 100644
--- a/src/prog/main.F90
+++ b/src/prog/main.F90
@@ -1595,6 +1595,25 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
             call env%error("No solvent name provided for COSMO", source)
          end if
       
+      case('--tmcosmo')
+         call args%nextArg(sec)
+         if (allocated(sec)) then
+            call set_gbsa(env, 'solvent', sec)
+            call set_gbsa(env, 'tmcosmo', 'true')
+            call args%nextArg(sec)
+            if (allocated(sec)) then
+               if (sec == 'reference') then
+                  gsolvstate = 1
+               else if (sec == 'bar1M') then
+                  gsolvstate = 2
+               else
+                  call env%warning("Unknown reference state '"//sec//"'", source)
+               end if
+            end if
+         else
+            call env%error("No solvent name provided for COSMO", source)
+         end if
+
       case('--cpcmx')
          if (get_xtb_feature('cpcmx')) then
             call args%nextArg(sec)
diff --git a/src/scf_module.F90 b/src/scf_module.F90
index 94bde4742..02adbad1a 100644
--- a/src/scf_module.F90
+++ b/src/scf_module.F90
@@ -854,7 +854,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, &
          type is (TCosmo)
             call open_file(ich, "xtb.cosmo", 'w')
             call solvation%writeCosmoFile(ich, mol%at, mol%sym, mol%xyz, &
-               & wfn%q, eel + ep + exb + merge(embd, ed + embd, allocated(scD4)))
+               & wfn%q, eel + ep + exb + merge(embd, ed + embd, allocated(scD4)),solvation%tmcosmo)
             call close_file(ich)
          end select
       end if
diff --git a/src/set_module.f90 b/src/set_module.f90
index d9e19ce60..99aaf4c1e 100644
--- a/src/set_module.f90
+++ b/src/set_module.f90
@@ -2020,6 +2020,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)
@@ -2073,8 +2074,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
 
diff --git a/src/solv/cosmo.f90 b/src/solv/cosmo.f90
index fd4e99761..a412c8487 100644
--- a/src/solv/cosmo.f90
+++ b/src/solv/cosmo.f90
@@ -95,6 +95,9 @@ module xtb_solv_cosmo
       real(wp), allocatable :: dsdr(:, :)
       real(wp), allocatable :: dsdrt(:, :, :)
 
+      !> TM convention?
+      logical :: tmcosmo = .false.
+
    contains
 
       !> Update coordinates and internal state
@@ -826,7 +829,7 @@ subroutine update_nnlist_sasa(nat, xyz, srcut, nnsas, nnlists)
 end subroutine update_nnlist_sasa
 
 !> Write a COSMO file output
-subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy)
+subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy, tmcosmo)
 
    !> COSMO container
    class(TCosmo), intent(in) :: self
@@ -849,11 +852,20 @@ subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy)
    !> Total energy
    real(wp), intent(in) :: energy
 
+   !> Switch to TM convention for cosmo file output
+   logical, intent(in), optional :: tmcosmo
+
    integer :: ii, ig, iat
    real(wp) :: dielEnergy, keps
    real(wp), allocatable :: phi(:), zeta(:), area(:)
+   logical :: tm
 
    allocate(phi(self%ddCosmo%ncav), zeta(self%ddCosmo%ncav), area(self%ddCosmo%ncav))
+   if (present(tmcosmo)) then
+      tm = tmcosmo
+   else
+      tm = .false.
+   end if
    ! Reset potential on the cavity, note that the potential is expected in e/Å
    call getPhi(qat, self%jmat, phi)
    ii = 0
@@ -870,6 +882,9 @@ subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy)
       end do
    end do
 
+   !! Switch convention for TM mode
+   if (tm) zeta=-zeta
+
 
     ! Dielectric energy is the energy on the dielectric continuum
    keps = 0.5_wp * (1.0_wp - 1.0_wp/self%dielectricConst)
@@ -940,6 +955,7 @@ subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy)
       & "#", &
       & "#"
 
+
    ii = 0
    do iat = 1, self%ddCosmo%nat
       do ig = 1, self%ddCosmo%ngrid
diff --git a/src/solv/input.F90 b/src/solv/input.F90
index 711f56624..b11918c38 100644
--- a/src/solv/input.F90
+++ b/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
 
diff --git a/src/solv/model.f90 b/src/solv/model.f90
index 322d1bee2..5e8e90734 100644
--- a/src/solv/model.f90
+++ b/src/solv/model.f90
@@ -693,6 +693,7 @@ 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)
+      cosmo%tmcosmo=self%tmcosmo
       call move_alloc(cosmo, model)
    else
       allocate(born)
diff --git a/src/xhelp.f90 b/src/xhelp.f90
index 024aad3d4..4bc70eba0 100644
--- a/src/xhelp.f90
+++ b/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.",&

From 9996906ad706c4fe2faa5c89b29c36e644b41df1 Mon Sep 17 00:00:00 2001
From: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
Date: Mon, 18 Sep 2023 14:51:05 +0200
Subject: [PATCH 2/6] Add ddvolume module for Volume output (and handling
 improvement)

Signed-off-by: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
---
 src/scf_module.F90      |   2 +-
 src/solv/CMakeLists.txt |   1 +
 src/solv/cosmo.f90      |  35 +++----
 src/solv/ddvolume.f90   | 207 ++++++++++++++++++++++++++++++++++++++++
 src/solv/meson.build    |   1 +
 src/solv/model.f90      |   5 +-
 6 files changed, 233 insertions(+), 18 deletions(-)
 create mode 100644 src/solv/ddvolume.f90

diff --git a/src/scf_module.F90 b/src/scf_module.F90
index 02adbad1a..94bde4742 100644
--- a/src/scf_module.F90
+++ b/src/scf_module.F90
@@ -854,7 +854,7 @@ subroutine scf(env, mol, wfn, basis, pcem, xtbData, solvation, &
          type is (TCosmo)
             call open_file(ich, "xtb.cosmo", 'w')
             call solvation%writeCosmoFile(ich, mol%at, mol%sym, mol%xyz, &
-               & wfn%q, eel + ep + exb + merge(embd, ed + embd, allocated(scD4)),solvation%tmcosmo)
+               & wfn%q, eel + ep + exb + merge(embd, ed + embd, allocated(scD4)))
             call close_file(ich)
          end select
       end if
diff --git a/src/solv/CMakeLists.txt b/src/solv/CMakeLists.txt
index d747696a3..98768e731 100644
--- a/src/solv/CMakeLists.txt
+++ b/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)
diff --git a/src/solv/cosmo.f90 b/src/solv/cosmo.f90
index a412c8487..2c0f07428 100644
--- a/src/solv/cosmo.f90
+++ b/src/solv/cosmo.f90
@@ -95,8 +95,11 @@ module xtb_solv_cosmo
       real(wp), allocatable :: dsdr(:, :)
       real(wp), allocatable :: dsdrt(:, :, :)
 
-      !> TM convention?
+      !> TM convention? 
       logical :: tmcosmo = .false.
+      !> Write volume in COSMO file? Volume resolution in Angstrom.
+      logical :: wrvolume = .false.
+      real(wp) :: volumeres = 0.1_wp
 
    contains
 
@@ -829,7 +832,8 @@ subroutine update_nnlist_sasa(nat, xyz, srcut, nnsas, nnlists)
 end subroutine update_nnlist_sasa
 
 !> Write a COSMO file output
-subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy, tmcosmo)
+subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy)
+   use xtb_solv_ddvolume, only: ddvolume
 
    !> COSMO container
    class(TCosmo), intent(in) :: self
@@ -852,20 +856,11 @@ subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy, tmcosmo)
    !> Total energy
    real(wp), intent(in) :: energy
 
-   !> Switch to TM convention for cosmo file output
-   logical, intent(in), optional :: tmcosmo
-
    integer :: ii, ig, iat
    real(wp) :: dielEnergy, keps
-   real(wp), allocatable :: phi(:), zeta(:), area(:)
-   logical :: tm
+   real(wp), allocatable :: phi(:), zeta(:), area(:), volume(:)
 
    allocate(phi(self%ddCosmo%ncav), zeta(self%ddCosmo%ncav), area(self%ddCosmo%ncav))
-   if (present(tmcosmo)) then
-      tm = tmcosmo
-   else
-      tm = .false.
-   end if
    ! Reset potential on the cavity, note that the potential is expected in e/Å
    call getPhi(qat, self%jmat, phi)
    ii = 0
@@ -883,7 +878,8 @@ subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy, tmcosmo)
    end do
 
    !! Switch convention for TM mode
-   if (tm) zeta=-zeta
+   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
@@ -901,9 +897,16 @@ subroutine writeCosmoFile(self, unit, num, sym, xyz, qat, energy, tmcosmo)
 
    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", &
diff --git a/src/solv/ddvolume.f90 b/src/solv/ddvolume.f90
new file mode 100644
index 000000000..754ef0e0c
--- /dev/null
+++ b/src/solv/ddvolume.f90
@@ -0,0 +1,207 @@
+! 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-decomposited 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
+    
+        !! Check if xyz and r have the same size
+        if (size(xyz,2) /= size(r)) then
+            write(*,*) "Error: xyz and r must have the same size"
+            write(*,*) "size(xyz,2) = ", size(xyz,2)
+            write(*,*) "size(r) = ", size(r)
+            stop
+        end if
+
+        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)
+
+        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
+    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
+
+
+
+
diff --git a/src/solv/meson.build b/src/solv/meson.build
index ef980fa44..aab3988db 100644
--- a/src/solv/meson.build
+++ b/src/solv/meson.build
@@ -29,4 +29,5 @@ srcs += files(
   'model.f90',
   'sasa.f90',
   'state.f90',
+  'ddvolume.f90',
 )
diff --git a/src/solv/model.f90 b/src/solv/model.f90
index 5e8e90734..51708c4b4 100644
--- a/src/solv/model.f90
+++ b/src/solv/model.f90
@@ -693,7 +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)
-      cosmo%tmcosmo=self%tmcosmo
+      if (self%tmcosmo) then
+         cosmo%tmcosmo=.true.
+         cosmo%wrvolume=.true.
+      end if
       call move_alloc(cosmo, model)
    else
       allocate(born)

From 63deaf27e3eed49beedd51a9ba6aab63bce4424d Mon Sep 17 00:00:00 2001
From: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
Date: Mon, 18 Sep 2023 14:56:30 +0200
Subject: [PATCH 3/6] Remove debugging code snippet.

Signed-off-by: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
---
 src/solv/ddvolume.f90 | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/src/solv/ddvolume.f90 b/src/solv/ddvolume.f90
index 754ef0e0c..d09018cba 100644
--- a/src/solv/ddvolume.f90
+++ b/src/solv/ddvolume.f90
@@ -57,14 +57,6 @@ subroutine ddvolume(xyz,r,resolution,volume)
         integer :: i,s,x,y,z
         real(wp) :: x_min,x_max,y_min,y_max,z_min,z_max
     
-        !! Check if xyz and r have the same size
-        if (size(xyz,2) /= size(r)) then
-            write(*,*) "Error: xyz and r must have the same size"
-            write(*,*) "size(xyz,2) = ", size(xyz,2)
-            write(*,*) "size(r) = ", size(r)
-            stop
-        end if
-
         allocate(spheres(size(xyz,2)))
         allocate(volume(size(xyz,2)))
         volume = 0.0_wp

From 44fc01118d9c55c3003a57918e6c95f3bde36044 Mon Sep 17 00:00:00 2001
From: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
Date: Mon, 18 Sep 2023 17:01:12 +0200
Subject: [PATCH 4/6] Include parallelization into ddvolume.

Signed-off-by: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
---
 src/solv/ddvolume.f90 | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/solv/ddvolume.f90 b/src/solv/ddvolume.f90
index d09018cba..9d407b57b 100644
--- a/src/solv/ddvolume.f90
+++ b/src/solv/ddvolume.f90
@@ -66,7 +66,7 @@ subroutine ddvolume(xyz,r,resolution,volume)
         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)
@@ -85,6 +85,7 @@ subroutine ddvolume(xyz,r,resolution,volume)
             end do
             volume(s) = spheres(s)%volume
       end do
+      !$OMP END PARALLEL DO
     end subroutine ddvolume
 
     !> Initializes a sphere

From 1dd1e9ed18c11b1bda567ab166dd66cd33d27251 Mon Sep 17 00:00:00 2001
From: albert <92109627+Albkat@users.noreply.github.com>
Date: Wed, 20 Sep 2023 11:43:23 +0200
Subject: [PATCH 5/6] misspelling

---
 src/solv/ddvolume.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/solv/ddvolume.f90 b/src/solv/ddvolume.f90
index 9d407b57b..3f7c14850 100644
--- a/src/solv/ddvolume.f90
+++ b/src/solv/ddvolume.f90
@@ -17,7 +17,7 @@
 ! MS, 09/2023: Initial implementation
 
 
-!> domain-decomposited volume calculation for molecule
+!> domain-decomposed volume calculation for molecule
 module xtb_solv_ddvolume
     use xtb_mctc_accuracy, only: wp
     implicit none

From 1f7a6cb63711d479a12f4f67f0f4de939acf3a01 Mon Sep 17 00:00:00 2001
From: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
Date: Wed, 20 Sep 2023 13:22:31 +0200
Subject: [PATCH 6/6] Remove unneccesary flag.

Signed-off-by: MtoLStoN <70513124+MtoLStoN@users.noreply.github.com>
---
 src/prog/main.F90 | 23 ++---------------------
 1 file changed, 2 insertions(+), 21 deletions(-)

diff --git a/src/prog/main.F90 b/src/prog/main.F90
index aef4b3c2e..122ce99cb 100644
--- a/src/prog/main.F90
+++ b/src/prog/main.F90
@@ -1576,11 +1576,11 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, 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
@@ -1595,25 +1595,6 @@ subroutine parseArguments(env, args, inputFile, paramFile, accuracy, lgrad, &
             call env%error("No solvent name provided for COSMO", source)
          end if
       
-      case('--tmcosmo')
-         call args%nextArg(sec)
-         if (allocated(sec)) then
-            call set_gbsa(env, 'solvent', sec)
-            call set_gbsa(env, 'tmcosmo', 'true')
-            call args%nextArg(sec)
-            if (allocated(sec)) then
-               if (sec == 'reference') then
-                  gsolvstate = 1
-               else if (sec == 'bar1M') then
-                  gsolvstate = 2
-               else
-                  call env%warning("Unknown reference state '"//sec//"'", source)
-               end if
-            end if
-         else
-            call env%error("No solvent name provided for COSMO", source)
-         end if
-
       case('--cpcmx')
          if (get_xtb_feature('cpcmx')) then
             call args%nextArg(sec)