1c 2C> \brief calculate the gradient terms due to the interaction with the 3C> COSMO charges 4C> 5C> Evaluate the gradient contributions from the COSMO embedding. The 6C> original part is from Klamt and Schüürmann [1] 7C> (see Eqs.(13-16)). The derivatives of matrix \f$A\f$ have been 8C> modified by York and Karplus [2] (see Eqs.(73-76)) to obtain smooth 9C> potential energy surfaces. York and Karplus also modified matrix 10C> \f$B\f$ which is easy to do in their classical force field code. 11C> In an ab-initio code this not so easy to do and as it is not 12C> required to eliminate singularities the original expression from [1] 13C> for \f$B\f$ is used here. 14C> 15C> ### References ### 16C> 17C> [1] A. Klamt, G. Schüürmann, 18C> "COSMO: a new approach to dielectric screening in solvents with 19C> explicit expressions for the screening energy and its gradient", 20C> <i>J. Chem. Soc., Perkin Trans. 2</i>, 1993, pp 799-805, DOI: 21C> <a href="https://doi.org/10.1039/P29930000799"> 22C> 10.1039/P29930000799</a>. 23C> 24C> [2] D.M. York, M. Karplus, 25C> "A smooth solvation potential based on the conductor-like 26C> screening model", <i>J. Phys. Chem. A</i> (1999) <b>103</b>, 27C> pp 11060-11079, DOI: 28C> <a href="https://doi.org/10.1021/jp992097l"> 29C> 10.1021/jp992097l</a>. 30C> 31 subroutine grad_hnd_cos ( H, lbuf, scr, lscr, 32 $ dens, frc_cos_nucq, frc_cos_elq, 33 $ frc_cos_qq, 34 $ g_dens, basis, geom, nproc, nat, 35 $ max_at_bf, rtdb, oskel ) 36c$Id: grad1.F 27488 2015-09-10 08:44:11Z alogsdail $ 37 38C COSMO one electron contribution to RHF, ROHF and UHF gradients 39C now also UMP2 ??? unlikely as that requires solutions to the 40c CPHF equation??? 41c 42c Terms included in this subroutine are: 43c 1. Electron - COSMO charge interactions 44c 2. Nuclear - COSMO charge interactions 45c 3. COSMO charge - COSMO charge interactions 46c 47c Terms NOT included are: 48c 1. All regular QM derivatives 49 50 implicit none 51 52#include "nwc_const.fh" 53#include "mafdecls.fh" 54#include "errquit.fh" 55#include "global.fh" 56#include "cosmoP.fh" 57#include "cosmo_params.fh" 58#include "geomP.fh" 59#include "geom.fh" 60#include "bq.fh" 61#include "bas.fh" 62#include "rtdb.fh" 63#include "sym.fh" 64#include "stdio.fh" 65#include "prop.fh" 66 67C-------------------------parameters-------------------------------- 68 integer g_dens !< [Input] the total electron density matrix GA 69 integer basis !< [Input] the basis set handle 70 integer geom !< [Input] the geometry handle 71 integer rtdb !< [Input] the RTDB handle 72 integer lbuf !< [Input] the length of the integral buffer 73 integer lscr, !< [Input] the length of the scratch space 74 $ nproc, nat, max_at_bf 75 76 double precision frc_cos_nucq(3,nat) !< [Output] the forces due 77 !< nuclear-COSMO charge 78 !< interaction 79 double precision frc_cos_elq(3,nat) !< [Output] the forces due 80 !< electron-COSMO charge 81 !< interaction 82 double precision frc_cos_qq(3,nat) !< [Output] the forces due 83 !< COSMO charge-COSMO charge 84 !< interaction 85 double precision H(lbuf) !< [Scratch] the derivative integrals 86 double precision scr(lscr) !< [Scratch] scratch space 87 double precision dens(max_at_bf,max_at_bf) !< [Scratch] local 88 !< density block 89 90 logical oskel ! symmetry? 91c 92 93C-------------------------local variables-------------------------- 94 95 integer ijatom, next, iat1, iat2, iat3, ish1, ish2, 96 $ iab1f, iab1l, iab2f, iab2l, iac1f, iac1l, iac2f, iac2l, 97 $ if1, il1, if2, il2, 98 $ icart, ic, nint, ip1, ip2 99 integer im1, im2, nprim, ngen, sphcart, ityp1, ityp2 100 integer ich1, ich2 101 102 integer nefc ! the number of COSMO charges 103 integer l_efciat ! the handle of the COSMO charge-atom map 104 integer k_efciat ! the index of the COSMO charge-atom map 105 integer l_rad ! the handle of the atom radii 106 integer k_rad ! the index of the atom radii 107 integer nefcl ! the number of COSMO charge for a given atom 108 integer iefc ! counter over COSMO charges 109 integer iefc_c ! memory index for COSMO charge coordinates 110 integer iefc_q ! memory index for COSMO charges 111 112 double precision dE, qfac, fact, dx, dy, dz, rr 113 double precision invscreen, bohr, pi 114 double precision zeta1, zeta2, zeta12 115 parameter (bohr=0.529177249d0) 116 117 logical status, pointforce 118 119 double precision util_erf, cosff, cosdff 120 external util_erf, cosff, cosdff 121 122 integer nxtask, task_size 123 external nxtask 124c 125 double precision rin, rout, alphai, xyzff 126 integer iat 127 parameter (alphai = 0.5d0) 128 rin(iat)=dbl_mb(k_rad-1+iat) 129 & *(1.0d0-alphai*gammas*sqrt(0.25d0**minbem)) 130 rout(iat)=dbl_mb(k_rad-1+iat) 131 & *(1.0d0+(1.0d0-alphai)*gammas*sqrt(0.25d0**minbem)) 132 133c ---- -cosmo- gradient term ----- 134 logical odbug 135 136 pi = acos(-1.0d0) 137 odbug=.false. 138 if(odbug) then 139 write(Luout,*) 'in -grad1_hnd_cos- ...' 140 endif 141c 142 task_size = 1 143 status = rtdb_parallel(.true.) ! Broadcast reads to all processes 144c 145 pointforce = geom_include_bqbq(geom) 146 if (.not.bq_create('cosmo efc bq',cosmo_bq_efc)) 147 $ call errquit("grad_hnd_cos: bq_create on cosmo failed", 148 $ 0,GEOM_ERR) 149 if (.not.bq_rtdb_load(rtdb,cosmo_bq_efc)) 150 $ call errquit('grad_hnd_cos: rtdb load failed for Bq',916, 151 $ rtdb_err) 152 if (.not.bq_ncenter(cosmo_bq_efc,nefc)) 153 $ call errquit('grad_hnd_cos: could not retrieve nefc',917, 154 $ GEOM_ERR) 155 if (.not.bq_index_coord(cosmo_bq_efc,iefc_c)) 156 $ call errquit('grad_hnd_cos: could not get coordinate index Bq', 157 $ cosmo_bq_efc,MA_ERR) 158 if (.not.bq_index_charge(cosmo_bq_efc,iefc_q)) 159 $ call errquit('grad_hnd_cos: could not get charge index Bq', 160 $ cosmo_bq_efc,MA_ERR) 161c 162 if (.not.ma_push_get(MT_DBL,nat,"rad",l_rad,k_rad)) 163 $ call errquit("grad_hnd_cos: could not allocate rad", 164 $ ma_sizeof(MT_BYTE,nat,MT_DBL),MA_ERR) 165 call cosmo_def_radii(rtdb,geom,nat,dbl_mb(k_rad)) 166 status = rtdb_get(rtdb,'cosmo:radius',mt_dbl,nat, 167 $ dbl_mb(k_rad)) 168 do iat1=0,nat-1 169 dbl_mb(k_rad+iat1) = dbl_mb(k_rad+iat1)/bohr 170 enddo 171c 172 if (.not.ma_push_get(MT_INT,nefc,"efciat",l_efciat,k_efciat)) 173 $ call errquit("grad_hnd_cos: could not allocate efciat", 174 $ ma_sizeof(MT_BYTE,nefc,MT_INT),MA_ERR) 175 if(.not.rtdb_get(rtdb,'cosmo:efciat',mt_int,nefc, 176 $ int_mb(k_efciat))) 177 $ call errquit('grad_hnd_cos: rtdb get failed for iatefc',915, 178 $ rtdb_err) 179c 180 call hf_print_set(1) 181 182 ijatom = -1 183 next = nxtask(nproc,task_size) 184 do 90, iat1 = 1, nat 185 do 80, iat2 = 1, iat1 186 187 ijatom = ijatom + 1 188 if ( ijatom .eq. next ) then 189 190 status = bas_ce2bfr(basis,iat1,iab1f,iab1l) 191 status = bas_ce2bfr(basis,iat2,iab2f,iab2l) 192 193 if (iab1f.le.0 .or. iab2f.le.0) then 194c 195c At least one center has no functions on it ... next atom 196c 197 goto 1010 198 endif 199 200 if (oskel) then 201 if (.not. sym_atom_pair(geom, iat1, iat2, qfac)) 202 $ goto 1010 203 else 204 qfac = 1.0d0 205 endif 206 qfac = -qfac 207 208 status = bas_ce2cnr(basis,iat1,iac1f,iac1l) 209 status = bas_ce2cnr(basis,iat2,iac2f,iac2l) 210 211 call ga_get (g_dens, iab1f,iab1l,iab2f,iab2l,dens,max_at_bf) 212 213 do 70, ish1 = iac1f, iac1l 214 if ( iat1.eq.iat2 ) iac2l = ish1 215 do 60, ish2 = iac2f, iac2l 216 217C shell block in atomic (D/Dw)-matrix block 218 status = bas_cn2bfr(basis,ish1,if1,il1) 219 if1 = if1 - iab1f + 1 220 il1 = il1 - iab1f + 1 221c 222c Work out the number of Cartesian basis functions 223c The integrals are evaluated in the Cartesian basis set 224c and then transformed to spherical harmonics. So the 225c buffer size depends on the number of Cartesian functions 226c 227 status = bas_continfo(basis,ish1,ityp1,nprim,ngen, 228 + sphcart) 229 if (sphcart.eq.1.and.ityp1.ge.2) then 230 im1 = if1 + (ityp1+1)*(ityp1+2)/2 - 1 231 else 232 im1 = il1 233 endif 234 status = bas_cn2bfr(basis,ish2,if2,il2) 235 if2 = if2 - iab2f + 1 236 il2 = il2 - iab2f + 1 237c 238c Same Cartesian vs spherical harmonic catastrophy as 239c for ish1. 240c 241 status = bas_continfo(basis,ish2,ityp2,nprim,ngen, 242 + sphcart) 243 if (sphcart.eq.1.and.ityp2.ge.2) then 244 im2 = if2 + (ityp2+1)*(ityp2+2)/2 - 1 245 else 246 im2 = il2 247 endif 248 249 nint = ( im1 - if1 + 1 ) * ( im2 - if2 + 1 ) 250 251 do iefc = 1, nefc 252 253 ic=1 254 do iat3 = 1, 3 ! centers A, B, and C 255 do icart = 1, 3 256 do ip1 = if1, im1 257 do ip2 = if2, im2 258 H(ic)=0.0D0 259 ic = ic + 1 260 enddo 261 enddo 262 enddo 263 enddo 264 265C 1el. -cosmo- derivatives 266c Currently calculated on for every COSMO charge 267c separately. 268 call intd_1epot_cosmo(basis,ish1,basis,ish2,lscr,scr, 269 & lbuf,H,dbl_mb(iefc_c+3*(iefc-1)), 270 & dbl_mb(iefc_q+iefc-1),1) 271 272C D x H 273c 274c Do center A (associated with ish1) 275c 276 ic=1 277 do icart = 1, 3 278 dE = 0.D0 279 do ip1 = if1, il1 280 do ip2 = if2, il2 281 dE = dE + dens(ip1,ip2) * H(ic) 282 ic = ic + 1 283 enddo 284 enddo 285 if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE 286 dE = dE * qfac 287 frc_cos_elq(icart,iat1) 288 & = frc_cos_elq(icart,iat1) - dE 289 enddo 290c 291c Do center B (associated with ish2) 292c 293 do icart = 1, 3 294 dE = 0.D0 295 do ip1 = if1, il1 296 do ip2 = if2, il2 297 dE = dE + dens(ip1,ip2) * H(ic) 298 ic = ic + 1 299 enddo 300 enddo 301 if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) then 302 dE = dE + dE 303 else 304 dE = 0.0d0 305 endif 306 dE = dE * qfac 307 frc_cos_elq(icart,iat2) 308 & = frc_cos_elq(icart,iat2) - dE 309 enddo 310c 311c Do center C, i.e. the Cosmo charge (associated with 312c the atom stored in int_mb(k_efciat)) 313c 314 do icart = 1, 3 315 dE = 0.D0 316 do ip1 = if1, il1 317 do ip2 = if2, il2 318 dE = dE + dens(ip1,ip2) * H(ic) 319 ic = ic + 1 320 enddo 321 enddo 322 if ( iat1.ne.iat2 .or. ish1.ne.ish2 ) dE = dE + dE 323 dE = dE * qfac 324 frc_cos_elq(icart,int_mb(k_efciat+iefc-1)) 325 & = frc_cos_elq(icart,int_mb(k_efciat+iefc-1)) - dE 326 enddo 327 328 enddo 329 330 60 continue 331 70 continue 332 333 1010 continue 334 335 next = nxtask(nproc,task_size) 336 endif 337 338 80 continue 339 90 continue 340 next = nxtask(-nproc,task_size) 341c 342 if (ga_nodeid().eq.0) then 343c 344c Do the Nuclear - Cosmo charge part 345c - 1. The derivative of matrix B (i.e. the Coulomb interaction 346c between the nuclear charge and the surface charge). 347c - 2. The derivative due to the change in the switching 348c function (see [2] Eq.(74)). This only applies for the 349c York-Karplus model. 350c 351 invscreen = 1.0d0/(1.0d0*screen) 352 do ich1 = 1, nefc 353 if (do_cosmo_model.eq.DO_COSMO_YK) then 354 zeta1 = zeta*sqrt(ptspatm) 355 $ / (dbl_mb(k_rad+int_mb(k_efciat+ich1-1)-1) 356 $ * sqrt(2.0d0*pi)) 357 xyzff = 1.0d0 358 do iat1 = 1, nat 359 if (iat1.ne.int_mb(k_efciat+ich1-1)) then 360 dx = coords(1,iat1,geom) - dbl_mb(0+3*(ich1-1)+iefc_c) 361 dy = coords(2,iat1,geom) - dbl_mb(1+3*(ich1-1)+iefc_c) 362 dz = coords(3,iat1,geom) - dbl_mb(2+3*(ich1-1)+iefc_c) 363 rr = sqrt(dx*dx+dy*dy+dz*dz) 364 if (rr.lt.rout(iat1)) then 365 xyzff = xyzff 366 $ * cosff((rr-rin(iat1))/(rout(iat1)-rin(iat1))) 367 endif 368 endif 369 enddo 370 endif 371 do iat1 = 1, nat 372 dx = coords(1,iat1,geom) - dbl_mb(0+3*(ich1-1)+iefc_c) 373 dy = coords(2,iat1,geom) - dbl_mb(1+3*(ich1-1)+iefc_c) 374 dz = coords(3,iat1,geom) - dbl_mb(2+3*(ich1-1)+iefc_c) 375 rr = sqrt(dx*dx+dy*dy+dz*dz) 376c 377c - term 1. 378c 379 fact = -charge(iat1,geom)*dbl_mb(iefc_q+ich1-1) / 380 $ rr**3 381c 382c - term 2. 383c 384 if (do_cosmo_model.eq.DO_COSMO_YK) then 385 if (iat1.ne.int_mb(k_efciat+ich1-1)) then 386 fact = fact - 2.0d0*zeta1*sqrt(2.0d0/pi)*invscreen 387 $ * dbl_mb(iefc_q+ich1-1)**2 388 $ * cosdff((rr-rin(iat1))/(rout(iat1)-rin(iat1))) 389 $ / (rr*xyzff*(rout(iat1)-rin(iat1)) 390 $ *cosff((rr-rin(iat1))/(rout(iat1)-rin(iat1)))) 391 endif 392 endif 393c 394 dE = dx * fact 395 frc_cos_nucq(1,iat1) = frc_cos_nucq(1,iat1) + dE 396 frc_cos_nucq(1,int_mb(k_efciat+ich1-1)) 397 $ = frc_cos_nucq(1,int_mb(k_efciat+ich1-1)) - dE 398 dE = dy * fact 399 frc_cos_nucq(2,iat1) = frc_cos_nucq(2,iat1) + dE 400 frc_cos_nucq(2,int_mb(k_efciat+ich1-1)) 401 $ = frc_cos_nucq(2,int_mb(k_efciat+ich1-1)) - dE 402 dE = dz * fact 403 frc_cos_nucq(3,iat1) = frc_cos_nucq(3,iat1) + dE 404 frc_cos_nucq(3,int_mb(k_efciat+ich1-1)) 405 $ = frc_cos_nucq(3,int_mb(k_efciat+ich1-1)) - dE 406 enddo 407 enddo 408c 409c Do cosmo charge - cosmo charge interaction 410c 411 invscreen = 1.0d0/(1.0d0*screen) 412 do ich1 = 1, nefc 413 if (do_cosmo_model.eq.DO_COSMO_YK) then 414 zeta1 = zeta*sqrt(ptspatm) 415 $ / (dbl_mb(k_rad+int_mb(k_efciat+ich1-1)-1) 416 $ * sqrt(2.0d0*pi)) 417 endif 418 do ich2 = 1, ich1-1 419 if (do_cosmo_model.eq.DO_COSMO_YK) then 420 zeta2 = zeta*sqrt(ptspatm) 421 $ / (dbl_mb(k_rad+int_mb(k_efciat+ich2-1)-1) 422 $ * sqrt(2.0d0*pi)) 423 zeta12 = zeta1*zeta2/sqrt(zeta1**2+zeta2**2) 424 endif 425 if (ich1.ne.ich2) then 426 dx = dbl_mb(0+3*(ich1-1)+iefc_c) 427 $ - dbl_mb(0+3*(ich2-1)+iefc_c) 428 dy = dbl_mb(1+3*(ich1-1)+iefc_c) 429 $ - dbl_mb(1+3*(ich2-1)+iefc_c) 430 dz = dbl_mb(2+3*(ich1-1)+iefc_c) 431 $ - dbl_mb(2+3*(ich2-1)+iefc_c) 432 rr = sqrt(dx*dx+dy*dy+dz*dz) 433 if (rr.lt.1.0d-6) then 434 fact = 0.0d0 435 else 436 if (do_cosmo_model.eq.DO_COSMO_KS) then 437 fact = -invscreen*dbl_mb(iefc_q+ich1-1) 438 $ * dbl_mb(iefc_q+ich2-1)/(rr**3) 439 else if (do_cosmo_model.eq.DO_COSMO_YK) then 440 fact = +invscreen*dbl_mb(iefc_q+ich1-1) 441 $ * dbl_mb(iefc_q+ich2-1) 442 $ * (2.0d0*zeta12/sqrt(pi)/(rr**2) 443 $ * exp(-(zeta12*rr)**2) 444 $ - util_erf(zeta12*rr)/(rr**3)) 445 else 446 call errquit("grad_hnd_cos: panic", 447 $ do_cosmo_model,UERR) 448 endif 449 endif 450 dE = dx * fact 451 frc_cos_qq(1,int_mb(k_efciat+ich1-1)) 452 $ = frc_cos_qq(1,int_mb(k_efciat+ich1-1)) + dE 453 frc_cos_qq(1,int_mb(k_efciat+ich2-1)) 454 $ = frc_cos_qq(1,int_mb(k_efciat+ich2-1)) - dE 455 dE = dy * fact 456 frc_cos_qq(2,int_mb(k_efciat+ich1-1)) 457 $ = frc_cos_qq(2,int_mb(k_efciat+ich1-1)) + dE 458 frc_cos_qq(2,int_mb(k_efciat+ich2-1)) 459 $ = frc_cos_qq(2,int_mb(k_efciat+ich2-1)) - dE 460 dE = dz * fact 461 frc_cos_qq(3,int_mb(k_efciat+ich1-1)) 462 $ = frc_cos_qq(3,int_mb(k_efciat+ich1-1)) + dE 463 frc_cos_qq(3,int_mb(k_efciat+ich2-1)) 464 $ = frc_cos_qq(3,int_mb(k_efciat+ich2-1)) - dE 465 endif 466 enddo 467 enddo 468 endif 469c 470 if (.not.ma_pop_stack(l_efciat)) 471 $ call errquit("grad_hnd_cos: could not deallocate l_efciat", 472 $ 0,MA_ERR) 473 if (.not.ma_pop_stack(l_rad)) 474 $ call errquit("grad_hnd_cos: could not deallocate l_rad", 475 $ 0,MA_ERR) 476 if (.not.bq_destroy(cosmo_bq_efc)) 477 $ call errquit("grad_hnd_cos: bq_destroy on cosmo failed", 478 $ 0,GEOM_ERR) 479 480 return 481 end 482