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