1subroutine test_eeq_model_gbsa
2   use xtb_mctc_accuracy, only : wp
3   use xtb_mctc_io, only : stdout
4   use assertion
5   use xtb_type_environment
6   use xtb_type_molecule
7   use xtb_type_param
8   use xtb_solv_gbsa
9   use xtb_solv_input
10   use xtb_solv_model
11   use xtb_solv_state
12   use xtb_disp_coordinationnumber, only : getCoordinationNumber, cnType
13   use xtb_eeq
14   use xtb_chargemodel
15   use xtb_solv_kernel
16   implicit none
17
18   real(wp),parameter :: thr = 1.0e-10_wp
19   integer, parameter :: nat = 8
20   integer, parameter :: at(nat) = [82,1,1,1,1,52,1,1]
21   real(wp),parameter :: xyz(3,nat) = reshape(&
22      &[0.00607452450000_wp,-4.06569893412500_wp, 0.00000000000000_wp, &
23      &-0.10304014450000_wp,-0.81095583012500_wp, 0.00000000000000_wp, &
24      &-1.48984957750000_wp,-5.20124851012500_wp, 2.65857423800000_wp, &
25      &-1.48984957750000_wp,-5.20124851012500_wp,-2.65857423800000_wp, &
26      & 3.11171824650000_wp,-5.04151186112500_wp, 0.00000000000000_wp, &
27      & 0.00607452450000_wp, 5.31621802187500_wp, 0.00000000000000_wp, &
28      & 2.16097631450000_wp, 7.52891897087500_wp, 0.00000000000000_wp, &
29      &-2.20210431050000_wp, 7.47552665287500_wp, 0.00000000000000_wp], shape(xyz))
30   integer, parameter :: iunit = stdout
31   real(wp),parameter :: temp = 298.15_wp
32
33   type(TMolecule) :: mol
34   type(chrg_parameter) :: chrgeq
35   type(TBorn) :: gbsa
36   type(TSolvInput) :: input
37   type(TSolvModel) :: model
38   type(TEnvironment) :: env
39   real(wp), parameter :: trans(3, 1) = 0.0_wp
40   real(wp) :: es,gsolv,sigma(3,3)
41   real(wp),allocatable :: cn(:),dcndr(:,:,:),dcndL(:,:,:)
42   real(wp),allocatable :: q(:),dqdr(:,:,:),dqdL(:,:,:),ges(:,:)
43
44   allocate(cn(nat),dcndr(3,nat,nat),dcndL(3,3,nat), &
45      & q(nat),dqdr(3,nat,nat),ges(3,nat))
46   es  = 0.0_wp
47   ges = 0.0_wp
48
49   call init(env)
50   call init(mol, at, xyz)
51
52   call getCoordinationNumber(mol, trans, 40.0_wp, cnType%erf, cn, dcndr, dcndL)
53   call assert_close(cn(1),3.9364408722599_wp,thr)
54   call assert_close(cn(6),1.9734817313511_wp,thr)
55
56   ! test CN derivative, correct summation of diagonals
57   call assert_close(dcndr(2,1,5),-0.30966754386801E-01_wp,thr)
58   call assert_close(dcndr(1,3,2),-0.00000000000000E+00_wp,thr)
59   call assert_close(dcndr(2,1,6), 0.24290126624329E-05_wp,thr)
60
61   call new_charge_model_2019(chrgeq,mol%n,mol%at)
62   call eeq_chrgeq(mol,env,chrgeq,cn,dcndr,dcndL,q,dqdr,dqdL,es,ges,sigma, &
63      &            .false.,.true.,.true.)
64
65   ! test electrostatic energy
66   call assert_close(es,-0.10559262715710E-01_wp,thr)
67
68   ! test molecular gradient of ES, also check symmetry
69   call assert_close(ges(2,5), 0.45485199335534E-04_wp,thr)
70   call assert_close(ges(3,3),-0.62757274919919E-04_wp,thr)
71   call assert_close(ges(1,2), 0.70062774412852E-05_wp,thr)
72
73   ! test for charge constraint
74   call assert_close(sum(q),0.0_wp,            thr)
75   call assert_close(q(1),0.21873929670215_wp,thr)
76
77   ! test derivative of partial charges
78   call assert_close(dqdr(2,1,6),-0.32458122507411E-02_wp,thr)
79   call assert_close(dqdr(3,2,1), 0.00000000000000E+00_wp,thr)
80   call assert_close(dqdr(1,3,2), 0.27025784902288E-03_wp,thr)
81   call assert_close(dqdr(2,8,5), 0.36651303109315E-03_wp,thr)
82
83   es = 0.0_wp
84   ges = 0.0_wp
85   q = 0.0_wp
86   dqdr = 0.0_wp
87
88   input = TSolvInput(solvent='ch2cl2', alpb=.false., kernel=gbKernel%still, &
89      & state=solutionState%gsolv)
90   call init(model, env, input, 2)
91   call newBornModel(model, env, gbsa, mol%at)
92   call gbsa%update(env, mol%at, mol%xyz)
93
94   es = gbsa%gsasa + gbsa%gshift
95   ges = gbsa%dsdr
96   call eeq_chrgeq(mol,env,chrgeq,gbsa,cn,dcndr,q,dqdr,es,gsolv,ges, &
97      &            .false.,.true.,.true.)
98
99   ! test electrostatic energy
100   call assert_close(es,-0.18723100944150E-01_wp,thr)
101
102   ! test molecular gradient of ES, also check symmetry
103   call assert_close(ges(2,5), 0.15098175186158E-03_wp,thr)
104   call assert_close(ges(3,3),-0.42131892534669E-03_wp,thr)
105   call assert_close(ges(1,2), 0.28893943115425E-04_wp,thr)
106
107   ! test for charge constraint
108   call assert_close(sum(q),0.0_wp,thr)
109   call assert_close(q(1), 0.22653931001924E+00_wp,thr)
110   call assert_close(q(2),-0.20322834228298E-01_wp,thr)
111   call assert_close(q(3),-0.26032042423717E-01_wp,thr)
112   call assert_close(q(4),-0.26032042423717E-01_wp,thr)
113   call assert_close(q(5),-0.25976606121313E-01_wp,thr)
114   call assert_close(q(6),-0.18580233359855E+00_wp,thr)
115   call assert_close(q(7), 0.28815815576667E-01_wp,thr)
116   call assert_close(q(8), 0.28810733199693E-01_wp,thr)
117
118   !call assert_close(dqdr(2,1,6),-0.21098409244566E-02_wp,thr)
119   !call assert_close(dqdr(3,2,1), 0.00000000000000E+00_wp,thr)
120   !call assert_close(dqdr(1,3,2), 0.37653484040321E-03_wp,thr)
121   !call assert_close(dqdr(2,8,5), 0.82628766939384E-03_wp,thr)
122
123   call mol%deallocate
124
125   ! done: everythings fine
126   call terminate(afail)
127
128end subroutine test_eeq_model_gbsa
129
130subroutine test_eeq_model_hbond
131   use xtb_mctc_accuracy, only : wp
132   use xtb_mctc_io, only : stdout
133   use assertion
134   stop 77
135   call terminate(afail)
136end subroutine test_eeq_model_hbond
137
138subroutine test_eeq_model_salt
139   use xtb_mctc_accuracy, only : wp
140   use xtb_mctc_io, only : stdout
141   use xtb_mctc_convert, only : aatoau
142   use assertion
143   use xtb_type_environment
144   use xtb_type_molecule
145   use xtb_type_param
146   use xtb_solv_gbsa
147   use xtb_solv_input
148   use xtb_solv_model
149   use xtb_solv_state
150   use xtb_disp_coordinationnumber, only : getCoordinationNumber, cnType
151   use xtb_eeq
152   use xtb_chargemodel
153   use xtb_solv_kernel
154   implicit none
155
156   real(wp),parameter :: thr = 1.0e-10_wp
157   integer, parameter :: nat = 12
158   integer, parameter :: at(nat) = [9,1,6,6,6,6,6,6,9,1,1,9]
159   real(wp),parameter :: xyz(3,nat) = reshape([&
160      &-3.61707571665846_wp,-2.21326241566733_wp,-0.20914111039455_wp, &
161      &-4.12085912923274_wp, 2.64361027892242_wp,-0.77360733295900_wp, &
162      & 0.70571179547784_wp,-1.46804568562192_wp, 0.40260319330016_wp, &
163      &-1.70284271199247_wp,-0.55881754958407_wp,-0.07245619686109_wp, &
164      &-2.22721391221737_wp, 1.98624548249273_wp,-0.40689082126896_wp, &
165      &-0.21016646823535_wp, 3.64828289401216_wp,-0.24136575031326_wp, &
166      & 2.24393351740138_wp, 2.86731529508407_wp, 0.23334807696843_wp, &
167      & 2.63467443489279_wp, 0.29590222227677_wp, 0.54871426109542_wp, &
168      &-0.65510226134586_wp, 6.12711283170949_wp,-0.54503487549757_wp, &
169      & 3.78023043231610_wp, 4.20058933929303_wp, 0.35763959771486_wp, &
170      & 1.05737533649814_wp,-3.45504988937755_wp, 0.68084171455529_wp, &
171      & 4.98699896397970_wp,-0.51782531551223_wp, 1.02289296328971_wp], shape(xyz))
172   integer, parameter :: iunit = stdout
173   real(wp),parameter :: temp = 298.15_wp
174
175   type(TMolecule) :: mol
176   type(chrg_parameter) :: chrgeq
177   type(TBorn) :: gbsa
178   type(TSolvInput) :: input
179   type(TSolvModel) :: model
180   type(TEnvironment) :: env
181   real(wp), parameter :: trans(3, 1) = 0.0_wp
182   real(wp) :: es,gsolv,sigma(3,3)
183   real(wp),allocatable :: cn(:),dcndr(:,:,:),dcndL(:,:,:)
184   real(wp),allocatable :: q(:),dqdr(:,:,:),dqdL(:,:,:),ges(:,:)
185   real(wp),allocatable :: fgb(:,:),fhb(:)
186
187   allocate( cn(nat),dcndr(3,nat,nat),dcndL(3,3,nat), &
188      & q(nat),dqdr(3,nat,nat),ges(3,nat))
189   es  = 0.0_wp
190   ges = 0.0_wp
191
192   call init(env)
193   call init(mol, at, xyz)
194
195   call getCoordinationNumber(mol, trans, 40.0_wp, cnType%erf, cn, dcndr, dcndL)
196
197   call new_charge_model_2019(chrgeq,mol%n,mol%at)
198   call eeq_chrgeq(mol,env,chrgeq,cn,dcndr,dcndL,q,dqdr,dqdL,es,ges,sigma, &
199      &            .false.,.true.,.true.)
200
201   ! test electrostatic energy
202   call assert_close(es,-0.64149462344256E-01_wp,thr)
203
204   ! test molecular gradient of ES, also check symmetry
205   call assert_close(ges(2,7),-0.52237685114757E-03_wp,thr)
206   call assert_close(ges(3,3),-0.10347464904462E-03_wp,thr)
207   call assert_close(ges(1,9),-0.19202182317373E-03_wp,thr)
208
209   ! test for charge constraint
210   call assert_close(sum(q),0.0_wp,            thr)
211   call assert_close(q(1),-0.20644568417986E+00_wp,thr)
212   call assert_close(q(8), 0.83875312119020E-01_wp,thr)
213
214   ! test derivative of partial charges
215   call assert_close(dqdr(2,2,6),-0.93834552520241E-02_wp,thr)
216   call assert_close(dqdr(3,3,1),-0.24876158910325E-03_wp,thr)
217   call assert_close(dqdr(1,5,2),-0.15306250365591E-01_wp,thr)
218   call assert_close(dqdr(2,8,7),-0.92052318691156E-02_wp,thr)
219
220   es = 0.0_wp
221   ges = 0.0_wp
222   q = 0.0_wp
223   dqdr = 0.0_wp
224
225   input = TSolvInput(solvent='ch2cl2', alpb=.false., kernel=gbKernel%still, &
226      & state=solutionState%gsolv, ionStrength=0.001_wp, ionRad=1.0_wp*aatoau)
227   call init(model, env, input, 2)
228   call newBornModel(model, env, gbsa, mol%at)
229   call gbsa%update(env, mol%at, mol%xyz)
230
231   es = gbsa%gsasa + gbsa%gshift
232   ges = gbsa%dsdr
233   call eeq_chrgeq(mol,env,chrgeq,gbsa,cn,dcndr,q,dqdr,es,gsolv,ges, &
234      &            .false.,.true.,.true.)
235
236   ! test electrostatic energy
237   call assert_close(es,-0.73041696453862E-01_wp,thr)
238
239   ! test molecular gradient of ES, also check symmetry
240   call assert_close(ges(2,7),-0.54144355306156E-03_wp,thr)
241   call assert_close(ges(3,3),-0.13994375398115E-03_wp,thr)
242   call assert_close(ges(1,9),-0.72309374268259E-04_wp,thr)
243
244   ! test for charge constraint
245   call assert_close(sum(q),0.0_wp,thr)
246   call assert_close(q( 1),-0.22628341725089E+00_wp,thr)
247   call assert_close(q( 2), 0.17181969742365E+00_wp,thr)
248   call assert_close(q( 3),-0.31688852764536E-01_wp,thr)
249   call assert_close(q( 4), 0.86131431150190E-01_wp,thr)
250   call assert_close(q( 5),-0.31623215923196E-01_wp,thr)
251   call assert_close(q( 6), 0.86034414210273E-01_wp,thr)
252   call assert_close(q( 7),-0.31614086140679E-01_wp,thr)
253   call assert_close(q( 8), 0.86132083253082E-01_wp,thr)
254   call assert_close(q( 9),-0.22624788057207E+00_wp,thr)
255   call assert_close(q(10), 0.17181804147276E+00_wp,thr)
256   call assert_close(q(11), 0.17180510824691E+00_wp,thr)
257   call assert_close(q(12),-0.22628332310549E+00_wp,thr)
258
259   ! test derivative of partial charges (still incorrect)
260   !call assert_close(dqdr(2,2,6),-0.88602533712134E-02_wp,thr)
261   !call assert_close(dqdr(3,3,1),-0.38467542758235E-03_wp,thr)
262   !call assert_close(dqdr(1,5,2),-0.13407102470387E-01_wp,thr)
263   !call assert_close(dqdr(2,8,7),-0.84239153450647E-02_wp,thr)
264
265   call mol%deallocate
266
267   call terminate(afail)
268end subroutine test_eeq_model_salt
269