Skip to content

Commit 9096eee

Browse files
MtoLStoNMtoLStoN
authored
Prepare cosmo implementation for use with post scf methods. (#666)
1 parent 3915ff0 commit 9096eee

File tree

7 files changed

+243
-29
lines changed

7 files changed

+243
-29
lines changed

‎include/param_cosmo.fh‎

Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
type(gbsa_parameter),parameter :: gfn_cosmo = gbsa_parameter ( &
2+
76.36000000_wp, &
3+
18.00000000_wp, &
4+
1.00000000_wp, &
5+
1.00000000_wp, &
6+
1.30000000_wp, &
7+
0.00000000_wp, &
8+
0.00000000_wp, &
9+
0.00000000_wp, &
10+
(/ -0.08499967_wp, 0.46780225_wp, -2.87013596_wp, -3.95935069_wp, -0.29783987_wp, &
11+
-0.48323273_wp, 0.00133622_wp, 0.20448945_wp, 0.20150600_wp, 0.36379863_wp, &
12+
-3.47082133_wp, -0.93451053_wp, -1.46342018_wp, -0.32774697_wp, -0.38015204_wp, &
13+
-0.35311116_wp, -0.19972593_wp, -0.12891363_wp, -1.19450558_wp, -1.61289300_wp, &
14+
-4.68022533_wp, 3.22705056_wp, -1.76349679_wp, -1.78246886_wp, -1.41389733_wp, &
15+
-2.03424961_wp, -13.15585203_wp, -0.60232543_wp, -0.51255518_wp, -1.05389890_wp, &
16+
-0.80992362_wp, -0.04534252_wp, -0.48402589_wp, -0.34738945_wp, -0.31467674_wp, &
17+
-0.13545936_wp, -1.19984632_wp, -1.57439157_wp, -7.77905850_wp, -0.25386776_wp, &
18+
-1.36181862_wp, -1.89986750_wp, -1.66264106_wp, 0.63693473_wp, -0.42449931_wp, &
19+
-0.73041159_wp, -0.48383194_wp, -0.60860050_wp, -1.34579404_wp, -0.28314431_wp, &
20+
-0.58997478_wp, -0.50506060_wp, -0.36680822_wp, -0.28828614_wp, -1.50384745_wp, &
21+
-1.72772005_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
22+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
23+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
24+
0.00000000_wp, -0.95945946_wp, -0.97138135_wp, -2.83432331_wp, -2.01934203_wp, &
25+
-1.19296724_wp, -0.73245256_wp, -0.58420150_wp, -0.43386368_wp, -0.70607621_wp, &
26+
-0.73268812_wp, -1.00876694_wp, -0.52518304_wp, 0.00000000_wp, 0.00000000_wp, &
27+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
28+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp /), &
29+
(/ 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, &
30+
0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, &
31+
0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, &
32+
0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, &
33+
0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, &
34+
0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, &
35+
0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, &
36+
0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, &
37+
0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, &
38+
0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, &
39+
0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, 0.30000000_wp, &
40+
0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, &
41+
0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, &
42+
0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, &
43+
0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, &
44+
0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, &
45+
0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, &
46+
0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp, &
47+
0.80000000_wp, 0.80000000_wp, 0.80000000_wp, 0.80000000_wp /), &
48+
(/ 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
49+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
50+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
51+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
52+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
53+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
54+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
55+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
56+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
57+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
58+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
59+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
60+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
61+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
62+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
63+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
64+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
65+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp, &
66+
0.00000000_wp, 0.00000000_wp, 0.00000000_wp, 0.00000000_wp /) )

‎src/param/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@ list(APPEND srcs
2222
"${dir}/paulingen.f90"
2323
"${dir}/sqrtzr4r2.f90"
2424
"${dir}/vdwradd3.f90"
25+
"${dir}/vdwradcosmo.f90"
2526
)
2627

2728
set(srcs ${srcs} PARENT_SCOPE)

‎src/param/meson.build‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,4 +20,5 @@ srcs += files(
2020
'paulingen.f90',
2121
'sqrtzr4r2.f90',
2222
'vdwradd3.f90',
23+
'vdwradcosmo.f90',
2324
)

‎src/param/vdwradcosmo.f90‎

Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,103 @@
1+
! This file is part of xtb.
2+
!
3+
! Copyright (C) 2019-2020 Sebastian Ehlert
4+
!
5+
! xtb is free software: you can redistribute it and/or modify it under
6+
! the terms of the GNU Lesser General Public License as published by
7+
! the Free Software Foundation, either version 3 of the License, or
8+
! (at your option) any later version.
9+
!
10+
! xtb is distributed in the hope that it will be useful,
11+
! but WITHOUT ANY WARRANTY; without even the implied warranty of
12+
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+
! GNU Lesser General Public License for more details.
14+
!
15+
! You should have received a copy of the GNU Lesser General Public License
16+
! along with xtb. If not, see <https://www.gnu.org/licenses/>.
17+
18+
!> Cosmo van-der-Waals radii
19+
20+
module xtb_param_vdwradcosmo
21+
use xtb_mctc_accuracy, only : wp
22+
use xtb_mctc_convert, only : aatoau
23+
use xtb_mctc_symbols, only : toNumber
24+
implicit none
25+
private
26+
27+
public :: getVanDerWaalsRadCosmo, vanDerWaalsRadCosmo
28+
29+
30+
!> Get van-der-Waals Rad for a species
31+
interface getVanDerWaalsRadCosmo
32+
module procedure :: getVanDerWaalsRadCosmoSymbol
33+
module procedure :: getVanDerWaalsRadCosmoNumber
34+
end interface getVanDerWaalsRadCosmo
35+
36+
37+
!> Default value for unoptimized van-der-Waals radii
38+
real(wp), parameter :: cosmostub = 2.223_wp
39+
40+
!> COSMO optimized van-der-Waals radii
41+
real(wp), parameter :: vanDerWaalsRadCosmo(94) = aatoau * [ &
42+
& 1.3000_wp, 1.6380_wp, 1.5700_wp, 1.0530_wp, & ! h-be
43+
& 2.0480_wp, 2.0000_wp, 1.8300_wp, 1.7200_wp, & ! B-O
44+
& 1.7200_wp, 1.8018_wp, 1.8000_wp, 1.6380_wp, & ! F-Mg
45+
& 2.1530_wp, 2.2000_wp, 2.1060_wp, 2.1600_wp, & ! Al-S
46+
& 2.0500_wp, 2.2000_wp, 2.2230_wp, cosmostub, & ! Cl-Ca
47+
& cosmostub, 2.2930_wp, cosmostub, cosmostub, & ! Sc-Cr
48+
& cosmostub, cosmostub, cosmostub, cosmostub, & ! Mn-Ni
49+
& cosmostub, 1.6260_wp, cosmostub, 2.7000_wp, & ! Cu-Ge
50+
& 2.3500_wp, 2.2000_wp, 2.1600_wp, 2.3630_wp, & ! As-Kr
51+
& cosmostub, cosmostub, cosmostub, cosmostub, & ! Rb-Zr
52+
& cosmostub, cosmostub, cosmostub, cosmostub, & ! Nb-Ru
53+
& cosmostub, cosmostub, cosmostub, cosmostub, & ! Rh-Cd
54+
& 2.2580_wp, 2.5500_wp, 2.4100_wp, 2.4100_wp, & ! In-Te
55+
& 2.3200_wp, 2.5270_wp, cosmostub, cosmostub, & ! I-Ba
56+
& cosmostub, cosmostub, cosmostub, cosmostub, & ! La-Nd
57+
& cosmostub, cosmostub, cosmostub, cosmostub, & ! Pm-Gd
58+
& cosmostub, cosmostub, cosmostub, cosmostub, & ! Tb-Er
59+
& cosmostub, cosmostub, cosmostub, cosmostub, & ! Tm-Hf
60+
& cosmostub, cosmostub, cosmostub, cosmostub, & ! Ta-Os
61+
& cosmostub, cosmostub, cosmostub, cosmostub, & ! Ir-Hg
62+
& cosmostub, 2.3600_wp, 2.4220_wp, 2.3050_wp, & ! Tl-Po
63+
& 2.3630_wp, 2.5740_wp, cosmostub, cosmostub, & ! At-Ra
64+
& cosmostub, cosmostub, cosmostub, cosmostub, & ! Ac-U
65+
& cosmostub, cosmostub] ! Np-Pu
66+
67+
68+
contains
69+
70+
71+
!> Get van-der-Waals radius for species with a given symbol
72+
elemental function getVanDerWaalsRadCosmoSymbol(symbol) result(rad)
73+
74+
!> Element symbol
75+
character(len=*), intent(in) :: symbol
76+
77+
!> van-der-Waals radius
78+
real(wp) :: rad
79+
80+
rad = getVanDerWaalsRadCosmo(toNumber(symbol))
81+
82+
end function getVanDerWaalsRadCosmoSymbol
83+
84+
85+
!> Get van-der-Waals radius for species with a given atomic number
86+
elemental function getVanDerWaalsRadCosmoNumber(number) result(rad)
87+
88+
!> Atomic number
89+
integer, intent(in) :: number
90+
91+
!> van-der-Waals radius
92+
real(wp) :: rad
93+
94+
if (number > 0 .and. number <= size(vanDerWaalsRadCosmo, dim=1)) then
95+
rad = vanDerWaalsRadCosmo(number)
96+
else
97+
rad = -1.0_wp
98+
end if
99+
100+
end function getVanDerWaalsRadCosmoNumber
101+
102+
103+
end module xtb_param_vdwradcosmo

‎src/solv/model.f90‎

Lines changed: 46 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
12
! This file is part of xtb.
23
!
34
! Copyright (C) 2019-2020 Sebastian Ehlert
@@ -23,13 +24,15 @@ module xtb_solv_model
2324
use xtb_mctc_strings, only : lowercase
2425
use xtb_mctc_systools, only : rdpath
2526
use xtb_param_vdwradd3, only : vanDerWaalsRadD3
27+
use xtb_param_vdwradcosmo, only : vanDerWaalsRadCosmo
2628
use xtb_solv_gbsa, only : TBorn, init_ => init
2729
use xtb_solv_cosmo, only : TCosmo, init_ => init
2830
use xtb_solv_input, only : TSolvInput
2931
use xtb_solv_kernel, only : gbKernel
3032
use xtb_solv_state, only : solutionState, getStateShift
3133
use xtb_type_environment, only : TEnvironment
3234
use xtb_type_solvation, only : TSolvation
35+
use ieee_arithmetic, only : ieee_value,ieee_positive_inf
3336
implicit none
3437
private
3538

@@ -143,6 +146,8 @@ module xtb_solv_model
143146
include 'param_alpb_methanol.fh'
144147
include 'param_alpb_ethanol.fh'
145148

149+
include 'param_cosmo.fh'
150+
146151
!> Solvent density (g/cm^3) and molar mass (g/mol)
147152
real(wp), parameter :: molcm3toau = 8.92388e-2_wp
148153

@@ -191,7 +196,7 @@ subroutine initSolvModel(self, env, input, level)
191196
call getParamFile(env, solvent, 'cosmo', level, self%paramFile)
192197
end if
193198

194-
if (.not.allocated(self%paramFile)) then
199+
if (.not.allocated(self%paramFile) .and. (.not.input%cosmo)) then
195200
call getParamFile(env, solvent, 'gbsa', level, self%paramFile)
196201
end if
197202

@@ -238,7 +243,12 @@ subroutine loadInternalParam(self, env, solvent, level)
238243
!> Method level for solvation model
239244
integer, intent(in) :: level
240245

241-
type(gbsa_parameter), allocatable :: param
246+
!> I/O handling
247+
248+
integer :: iostat
249+
250+
type(gbsa_parameter), allocatable :: param, stub
251+
242252

243253
select case(self%kernel)
244254
case (gbKernel%still)
@@ -376,6 +386,34 @@ subroutine loadInternalParam(self, env, solvent, level)
376386
param = gfn1_alpb_hexane
377387
end select
378388
end if
389+
390+
if (self%cosmo) then
391+
if (allocated(param)) then
392+
self%paramFile = self%paramFile//"+COSMO"
393+
stub = param
394+
param = gfn_cosmo
395+
param%epsv = stub%epsv
396+
param%smass = stub%smass
397+
param%rhos = stub%rhos
398+
param%gamscale = stub%gamscale
399+
else if ((solvent .eq. "inf") .or. (solvent .eq. "infinity")) then
400+
self%paramFile = "internal GFN-xTB/COSMO"
401+
param = gfn_cosmo
402+
param%epsv = ieee_value(param%epsv,ieee_positive_inf)
403+
else
404+
self%paramFile = "internal GFN-xTB/COSMO"
405+
param = gfn_cosmo
406+
read(solvent,*,iostat=iostat) param%epsv
407+
if (iostat .ne. 0) then
408+
Call env%error(solvent// " is neither a solvent name nor a dielectric constant.",source)
409+
return
410+
end if
411+
if (param%epsv .eq. 0) then
412+
Call env%error("You really chose 0 as your dielectric constant?",source)
413+
return
414+
end if
415+
end if
416+
end if
379417

380418
if (.not.allocated(param)) then
381419
call env%error("solvent: '"//solvent//"' is not parametrized", source)
@@ -459,7 +497,12 @@ subroutine paramToModel(self, param)
459497
self%bornOffset = param%soset * 0.1_wp*aatoau
460498
self%probeRad = param%rprobe * aatoau
461499

462-
self%vdwRad = vanDerWaalsRadD3
500+
if (self%cosmo) then
501+
self%vdwRad=vanDerWaalsRadCosmo
502+
else
503+
self%vdwRad=vanDerWaalsRadD3
504+
end if
505+
463506
self%surfaceTension = param%gamscale * fourpi*surfaceTension
464507
self%descreening = param%sx
465508
allocate(self%hBondStrength(size(param%tmp)))

‎test/unit/test_gfn1.f90‎

Lines changed: 13 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -874,23 +874,23 @@ subroutine test_gfn1_mindless_cosmo(error)
874874
& "h2o", "chcl3", "thf", "acetonitrile", "toluene", &
875875
& "ch2cl2", "ether", "methanol", "cs2", "dmso"]
876876
real(wp), parameter :: ref_energies(10) = &
877-
&[-33.068149119769_wp,-26.879297949029_wp,-25.848339775486_wp, &
878-
& -24.920550779459_wp,-29.057643142396_wp,-20.616717615153_wp, &
879-
& -35.781900213083_wp,-33.108563594503_wp,-21.434739333961_wp, &
880-
& -26.668994477103_wp]
877+
&[-33.086826283006_wp,-26.892421337119_wp,-25.862366080135_wp, &
878+
& -24.929812068243_wp,-29.068535381601_wp,-20.631250827831_wp, &
879+
& -35.780608195840_wp,-33.119209493948_wp,-21.444833411574_wp, &
880+
& -26.674508354162_wp]
881881
real(wp), parameter :: ref_gnorms(10) = &
882-
&[0.048610989970_wp, 0.061037608875_wp, 0.051110430194_wp, &
883-
& 0.093510458216_wp, 0.048462578319_wp, 0.056499605067_wp, &
884-
& 0.062179565785_wp, 0.060672469106_wp, 0.041589454825_wp, &
885-
& 0.042982140079_wp]
882+
&[0.056910263744_wp, 0.059732400683_wp, 0.052197041661_wp, &
883+
& 0.097844384419_wp, 0.047961348424_wp, 0.057230114090_wp, &
884+
& 0.060889237549_wp, 0.068600226343_wp, 0.042102755940_wp, &
885+
& 0.042610278777_wp]
886886
real(wp), parameter :: ref_hlgaps(10) = &
887-
&[3.512016853496_wp, 1.675002766293_wp, 2.500980871958_wp, &
888-
& 1.659094616440_wp, 2.420893146112_wp, 3.802456715570_wp, &
889-
& 2.599886803836_wp, 1.354752587871_wp, 4.107155093819_wp, &
890-
& 1.379801456095_wp]
887+
&[3.560961680882_wp, 1.701178515210_wp, 2.575575698307_wp, &
888+
& 1.642444939544_wp, 2.389954524574_wp, 3.624912973388_wp, &
889+
& 2.770751688023_wp, 1.357191094207_wp, 4.102623151933_wp, &
890+
& 1.391461153653_wp]
891891

892892
call init(env)
893-
do iMol = 1, 2
893+
do iMol = 1,10
894894

895895
call getMolecule(mol, mindless(iMol))
896896

‎test/unit/test_gfn2.f90‎

Lines changed: 13 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -793,23 +793,23 @@ subroutine test_gfn2_mindless_cosmo(error)
793793
& "h2o", "chcl3", "thf", "acetonitrile", "toluene", &
794794
& "ch2cl2", "ether", "methanol", "cs2", "dmso"]
795795
real(wp), parameter :: ref_energies(10) = &
796-
&[-30.376518115875_wp,-24.102578449616_wp,-23.754139493763_wp, &
797-
& -22.766727017808_wp,-27.757889554228_wp,-18.587445483567_wp, &
798-
& -33.445663939387_wp,-30.003461899598_wp,-20.573488191967_wp, &
799-
& -25.670456523922_wp]
796+
&[-30.398257570644_wp,-24.125868305041_wp,-23.787177902522_wp, &
797+
& -22.779539008049_wp,-27.783054701847_wp,-18.610486346192_wp, &
798+
& -33.456896298026_wp,-30.008556230192_wp,-20.588957648989_wp, &
799+
& -25.684670492104_wp]
800800
real(wp), parameter :: ref_gnorms(10) = &
801-
&[0.065115699688_wp, 0.057742250527_wp, 0.041885898146_wp, &
802-
& 0.072134973561_wp, 0.048230827198_wp, 0.052250268149_wp, &
803-
& 0.045597725657_wp, 0.062065139853_wp, 0.054691676440_wp, &
804-
& 0.045829509846_wp]
801+
&[0.070117675966_wp, 0.060272066809_wp, 0.048063628868_wp, &
802+
& 0.076296084561_wp, 0.050041916284_wp, 0.051542825738_wp, &
803+
& 0.046066119467_wp, 0.060436009320_wp, 0.058023381592_wp, &
804+
& 0.050440402702_wp]
805805
real(wp), parameter :: ref_hlgaps(10) = &
806-
&[2.608938030439_wp, 1.390465293950_wp, 0.169294532274_wp, &
807-
& 1.254915268362_wp, 2.137937542332_wp, 2.267712835326_wp, &
808-
& 2.669493177775_wp, 0.840230142925_wp, 3.341412667263_wp, &
809-
& 0.393537743933_wp]
806+
&[3.173028297842_wp, 1.553006040621_wp, 0.324761161079_wp, &
807+
& 1.233091187215_wp, 1.918282889245_wp, 2.067664999523_wp, &
808+
& 2.757505697757_wp, 0.846361565116_wp, 3.336544801399_wp, &
809+
& 0.315118656281_wp]
810810

811811
call init(env)
812-
do iMol = 1, 2
812+
do iMol = 1, 10
813813

814814
call getMolecule(mol, mindless(iMol))
815815

0 commit comments

Comments
 (0)