1C> 2C> \ingroup cosmo 3C> @{ 4C> 5C> \file hnd_cosmo_lib.F 6C> Library of HONDO routines used to implement COSMO 7C> 8C> \brief Setup the COSMO cavity surface 9C> 10 subroutine hnd_cosset(rtdb,nat,c,radius) 11 implicit none 12c 13c ----- starting from -icosahedron- ----- 14c 15c pass, napex, nface, error = 0 12 20 20 16c pass, napex, nface, error = 1 42 80 100 0.4982 17c pass napex, nface, error = 2 162 320 420 0.1848 18c pass napex, nface, error = 3 642 1280 1700 0.0523 19c pass napex, nface, error = 4 2562 5120 6820 0.0135 20c pass napex, nface, error = 5 10242 20480 27300 0.0034 21c 22c ----- starting from -octahedron- ----- 23c 24c pass, napex, nface, error = 0 6 8 8 25c pass, napex, nface, error = 1 18 32 40 0.8075 26c pass napex, nface, error = 2 66 128 168 0.4557 27c pass napex, nface, error = 3 258 512 680 0.1619 28c pass napex, nface, error = 4 1026 2048 2728 0.0451 29c pass napex, nface, error = 5 4098 8192 10920 0.0116 30c pass napex, nface, error = 6 16386 32768 43688 0.0029 31c 32#include "cosmo_params.fh" 33#include "errquit.fh" 34#include "global.fh" 35#include "rtdb.fh" 36#include "mafdecls.fh" 37#include "nwc_const.fh" 38#include "stdio.fh" 39c 40 integer rtdb !< [Input] The RTDB handle 41 integer nat !< [Input] The number of atoms 42 double precision c(3,nat) !< [Input] The atomic coordinates 43 double precision radius(nat) !< [Input] The atomic radii 44c 45c 46 integer mxface, mxapex 47 parameter (mxface=43688) 48 parameter (mxapex=16386) 49 logical dbug, stat 50 integer l_i10, i10 51 integer l_i20, i20 52 integer l_i30, i30 53 integer l_i40, i40 54 integer l_i50, i50 55 integer l_i60, i60 56 integer l_i70, i70 57 integer l_i80, i80 58 integer l_i90, i90 59 integer l_i100, i100 60 integer l_i110, i110 61 integer l_i120, i120 62 integer l_i130, i130 63 integer need 64c 65 dbug=.false. 66 if(dbug.and.ga_nodeid().eq.0) then 67 write(luout,9999) 68 endif 69c 70 if(ificos.eq.0.and.maxbem.gt.6) then 71 write(luout,*) '-maxbem- too large for parameters in -cosset-' 72 call errquit('hnd_cosset, -maxbem- too large',911,0) 73 elseif(ificos.ne.0.and.maxbem.gt.7) then 74 write(luout,*) '-maxbem- too large for parameters in -cosset-' 75 call errquit('hnd_cosset, -maxbem- too large',911,0) 76 endif 77c 78c ----- partition memory ----- 79c 80 need = 6*nat + 7*mxface + 7*mxface*nat + 3*mxapex 81c 82c ----- allocate memory block ----- 83c 84c if(.not.ma_push_get(mt_dbl,need,'mem init:cosmo:hnd_cosset:1', 85c & i_init,init)) 86c & call errquit('hnd_cosset, malloc of init failed',911,MA_ERR) 87c 88 if(.not.ma_push_get(mt_dbl,3*nat,"xyzatm",l_i10,i10)) 89 c call errquit('hndcosset: not enuf mem',0,MA_ERR) 90 if(.not.ma_push_get(mt_dbl, nat,"ratm",l_i20,i20)) 91 c call errquit('hndcosset: not enuf mem',1,MA_ERR) 92 if(.not.ma_push_get(mt_int, nat,"nspa",l_i30,i30)) 93 c call errquit('hndcosset: not enuf mem',2,MA_ERR) 94 if(.not.ma_push_get(mt_int, nat,"nppa",l_i40,i40)) 95 c call errquit('hndcosset: not enuf mem',3,MA_ERR) 96 if(.not.ma_push_get(mt_int,3*mxface,"ijkfac",l_i50,i50)) 97 c call errquit('hndcosset: not enuf mem',4,MA_ERR) 98 if(.not.ma_push_get(mt_dbl,3*mxface,"xyzseg",l_i60,i60)) 99 c call errquit('hndcosset: not enuf mem',5,MA_ERR) 100 if(.not.ma_push_get(mt_int, mxface,"ijkseg",l_i70,i70)) 101 c call errquit('hndcosset: not enuf mem',6,MA_ERR) 102 if(.not.ma_push_get(mt_log, mxface*nat,"insseg",l_i80,i80)) 103 c call errquit('hndcosset: not enuf mem',7,MA_ERR) 104 if(.not.ma_push_get(mt_dbl,3*mxface*nat,"xyzspa",l_i90,i90)) 105 c call errquit('hndcosset: not enuf mem',8,MA_ERR) 106 if(.not.ma_push_get(mt_int, mxface*nat,"ijkspa",l_i100,i100)) 107 c call errquit('hndcosset: not enuf mem',9,MA_ERR) 108 if(.not.ma_push_get(mt_int, mxface*nat,"numpps",l_i110,i110)) 109 c call errquit('hndcosset: not enuf mem',10,MA_ERR) 110 if(.not.ma_push_get(mt_dbl,3*mxapex ,"apex",l_i120,i120)) 111 c call errquit('hndcosset: not enuf mem',11,MA_ERR) 112 if(.not.ma_push_get(mt_dbl, mxface*nat,"xyzff",l_i130,i130)) 113 c call errquit('hndcosset: not enuf mem',12,MA_ERR) 114c i10 =init ! xyzatm(3,nat) 115c i20 =i10 +3*nat ! ratm( nat) 116c i30 =i20 + nat ! nspa( nat) 117c i40 =i30 + nat ! nppa( nat) 118c i50 =i40 + nat ! ijkfac(3,mxface) 119c i60 =i50 +3*mxface ! xyzseg(3,mxface) 120c i70 =i60 +3*mxface ! ijkseg( mxface) 121c i80 =i70 + mxface ! insseg( mxface,nat) 122c i90 =i80 + mxface*nat ! xyzspa(3,mxface,nat) 123c i100=i90 +3*mxface*nat ! ijkspa( mxface,nat) 124c i110=i100+ mxface*nat ! numpps( mxface,nat) 125c i120=i110+ mxface*nat ! apex(3,mxapex) 126c 127c ----- get -cosmo- surface ----- 128c 129 call hnd_cossrf(nat,c,radius,nat,mxface,mxapex, 130 1 dbl_mb(i10),dbl_mb(i20),int_mb(i30),int_mb(i40), 131 2 int_mb(i50),dbl_mb(i60),int_mb(i70),log_mb(i80), 132 3 dbl_mb(i90),int_mb(i100),int_mb(i110), 133 4 dbl_mb(i120),dbl_mb(i130),rtdb) 134 135c 136c ----- release memory block ----- 137c 138 if(.not.ma_chop_stack(l_i10)) call 139 & errquit('hnd_cosset, ma_pop_stack of init failed',911,MA_ERR) 140c 141 return 142 9999 format(/,10x,15(1h-), 143 1 /,10x,'-cosmo- surface', 144 2 /,10x,15(1h-)) 145 end 146c 147C> \brief Generate the COSMO cavity surface 148C> 149 subroutine hnd_cossrf(nat,c,radius,mxatm,mxfac,mxapx, 150 1 xyzatm,ratm,nspa,nppa, 151 2 ijkfac,xyzseg,ijkseg,insseg, 152 3 xyzspa,ijkspa,numpps,apex,xyzff,rtdb) 153 implicit none 154c 155#include "nwc_const.fh" 156#include "cosmo_params.fh" 157#include "rtdb.fh" 158#include "global.fh" 159#include "stdio.fh" 160#include "cosmoP.fh" 161#include "mafdecls.fh" 162#include "util_params.fh" 163 integer rtdb, nat 164 integer mxatm 165 integer mxfac 166 integer mxapx 167 double precision c(3,nat ) 168 double precision radius( nat) 169 double precision xyzatm(3,mxatm) 170 double precision ratm( mxatm) 171 integer nspa( mxatm) 172 integer nppa( mxatm) 173 integer ijkfac(3,mxfac) 174 double precision xyzseg(3,mxfac) 175 integer ijkseg( mxfac) 176 logical insseg( mxfac,mxatm) 177 double precision xyzspa(3,mxfac,mxatm) 178 integer ijkspa( mxfac,mxatm) 179 integer numpps( mxfac,mxatm) 180 double precision apex(3,mxapx) 181 double precision xyzff( mxfac,mxatm) 182c 183 logical some 184 logical out 185 logical more 186 logical dbug 187c 188 integer i, iat, lfac, lseg, ndiv, nfac, nseg 189 integer mfac 190c 191 dbug=.false. 192 more=.false. 193 more=more.or.dbug 194 out =.false. 195 out =out.or.more 196 some=.false. 197 some=some.or.out 198c 199c ----- approximate sphere with segments and points ----- 200c 201 do iat = 1, mxatm 202 nspa(iat) = 0 203 nppa(iat) = 0 204 enddo 205 nseg = 0 206 nfac = 0 207 ndiv = 0 208 call hnd_cossph(nseg,nfac,ndiv, 209 1 ijkfac,xyzseg,ijkseg,mxfac,apex,mxapx, 210 2 dsurf,dvol,adiag) 211 ptspatm = dble(nseg) 212c 213c ----- debug printing ----- 214c 215 if(out.and.ga_nodeid().eq.0) then 216 write(luout,9999) nseg,nfac,ndiv,dsurf,dvol 217 write(luout,9995) adiag 218 if(more) then 219 write(luout,9998) 220 do lseg=1,nseg 221 write(luout,9997) lseg,xyzseg(1,lseg),xyzseg(2,lseg), 222 1 xyzseg(3,lseg),ijkseg( lseg) 223 enddo 224 endif 225 if(dbug) then 226 write(luout,9996) 227 do lfac=1,nfac 228 mfac=lfac+nseg 229 write(luout,9997) mfac,xyzseg(1,mfac),xyzseg(2,mfac), 230 1 xyzseg(3,mfac),ijkseg( mfac) 231 enddo 232 endif 233 endif 234c 235c ----- set molecule ----- 236c 237 do iat=1,nat 238 do i=1,3 239 xyzatm(i,iat)=c(i,iat) 240 enddo 241 enddo 242 do iat=1,nat 243 if(radius(iat).eq.0.0d0) then 244 ratm(iat)=0.0d0 245 else 246 if (do_cosmo_model.eq.DO_COSMO_KS) then 247 ratm(iat)=(radius(iat)+rsolv)/cau2ang 248 else 249 ratm(iat)=radius(iat)/cau2ang 250 endif 251 endif 252 enddo 253c 254c ----- create -solvent accessible surface- of the molecule ----- 255c 256 257 call hnd_cossas(nat,xyzatm,ratm,mxatm, 258 1 nspa,nppa,xyzspa,ijkspa, 259 2 nseg,nfac,xyzseg,ijkseg,insseg, 260 3 numpps,xyzff,mxfac,rtdb) 261c 262 return 263 9999 format(' nseg,nfac,ndiv=nfac/nseg,dsurf,dvol = ',3i7,2f10.6) 264 9998 format(' seg ',' x ',' y ',' z ', 265 1 ' seg ',/,1x,47(1h-)) 266 9997 format(i7,3f12.8,i5,f12.8) 267 9996 format(' fac ',' x ',' y ',' z ', 268 1 ' seg ',/,1x,47(1h-)) 269 9995 format(' adiag = ',f12.6) 270 end 271C> 272C> \brief Construct the Solvent Accessible Surface (SAS) from the 273C> triangulated spheres 274C> 275C> ## The legacy of Klamt and Schüürmann ## 276C> 277C> This subroutine was originally written to implement the algorithm 278C> to construct the Solvent Accessible Surface as proposed by 279C> Klamt and Schüürmann [1]. This algorithm worked as follows: 280C> 281C> If two spheres partially overlap then parts of the surface need 282C> to be eliminated. Segments that are contained entirely within the 283C> sphere of another atom will be eliminated completely. Segments 284C> that straddle the boundary between two spheres will have their 285C> surface reduced proportional to the fraction that resides within the 286C> sphere of the other atom. This fraction is established by counting 287C> the number of faces that fall within the sphere of the other atom. 288C> In addition the location of the charge representing a segment should 289C> be calculated as the center of the remaining points representing the 290C> faces (see [1] page 802, 2nd column, 2nd paragraph), but currently 291C> that is not done. 292C> 293C> To understand the approach suggested above it is important to know 294C> the concepts "segments" and "faces". 295C> - "Segments" are parts of the surface of the sphere that will be 296C> represented by a single COSMO charge. 297C> - "Faces" are further refinements of segments. I.e. segments have 298C> been partitioned into a number of faces. The faces are used to 299C> eliminate parts of a segment that are within the sphere of another 300C> atom. By counting the remaining faces the surface area of a segment 301C> can be adjusted. 302C> The trick with segments and faces is needed to create a smoother 303C> boundary between neighboring spheres without introducing large 304C> numbers of COSMO charges. The smooth boundary is needed to keep 305C> discretization errors small, whereas "small" numbers of COSMO charges 306C> are needed to keep the cost of calculating the COSMO charges low. 307C> 308C> The segments have been generated in `hnd_cossph` by partitioning the 309C> triangles of the original polyhedron `minbem` times. The faces have 310C> generated by partitioning the segments an additional `maxbem-minbem` 311C> times. 312C> 313C> ## The new smooth approach of York and Karplus ## 314C> 315C> The approach by Klamt and Schüürmann [1] led to problems 316C> because the corresponding potential energy surface was not 317C> continuous. York and Karplus [2] proposed a method that provides 318C> a smooth potential energy surface and this subroutine was changed 319C> to implement this new approach. This meant that some things stayed 320C> the same as before (for example minbem still works the same way 321C> to generate the surface charge positions), other things changed 322C> significantly (maxbem and rsolv are not used anymore, also the 323C> elimination of point charges is no longer based on reducing the 324C> segment surface area until it vanishes, instead the surface charge 325C> of a segment is quenched with a switching function to eliminate the 326C> contribution of a surface point). 327C> 328C> ## Popular demand ## 329C> 330C> Due to popular demand this routine can do either the original 331C> Klamt-Schuurmann approach or the York-Karplus approach. The approach 332C> used is dictated by the `do_cosmo_model` variable. 333C> 334C> ### References ### 335C> 336C> [1] A. Klamt, G. Schüürmann, 337C> "COSMO: a new approach to dielectric screening in solvents with 338C> explicit expressions for the screening energy and its gradient", 339C> <i>J. Chem. Soc., Perkin Trans. 2</i>, 1993, pp 799-805, DOI: 340C> <a href="https://doi.org/10.1039/P29930000799"> 341C> 10.1039/P29930000799</a>. 342C> 343C> [2] D.M. York, M. Karplus, 344C> "A smooth solvation potential based on the conductor-like 345C> screening model", <i>J. Phys. Chem. A</i> (1999) <b>103</b>, 346C> pp 11060-11079, DOI: 347C> <a href="https://doi.org/10.1021/jp992097l"> 348C> 10.1021/jp992097l</a>. 349C> 350 subroutine hnd_cossas(nat,xyzatm,ratm,mxatom, 351 1 nspa,nppa,xyzspa,ijkspa, 352 2 nseg,nfac,xyzseg,ijkseg,insseg, 353 3 numpps,xyzff,mxface,rtdb) 354 implicit none 355#include "cosmo_params.fh" 356#include "errquit.fh" 357#include "rtdb.fh" 358#include "mafdecls.fh" 359#include "global.fh" 360#include "stdio.fh" 361#include "bq.fh" 362#include "prop.fh" 363cnew 364#include "cosmoP.fh" 365#include "util_params.fh" 366c 367 integer rtdb !< [Input] The RTDB handle 368 integer nat !< [Input] The actual number of atoms 369 integer mxface !< [Input] The maximum number of faces 370 integer mxatom !< [Input] The maximum number of atoms 371 integer nseg !< [Input] The actual number of segments 372 integer nfac !< [Input] The actual number of faces 373c 374 logical dbug 375 double precision xyzatm(3,mxatom) !< [Input] The atom positions 376 double precision ratm( mxatom) !< [Input] The atom radii 377 integer nspa( mxatom) !< [Output] The number of 378 !< segments remaining for 379 !< each atom 380 integer nppa( mxatom) !< [Output] The number of faces 381 !< remaining for each atom 382 double precision xyzseg(3,mxface) !< [Input] The coordinates of 383 !< the segment and face points 384 !< on the unit sphere 385 integer ijkseg( mxface) !< [Input] List for every 386 !< face what the corresponding 387 !< segment is, if ijkseg(ii) is 388 !< 0 then face ii should be 389 !< ignored (has been eliminated) 390 logical insseg( mxface,mxatom) !< [Output] If .false. 391 !< keep the segment or 392 !< face, discard it 393 !< otherwise 394 double precision xyzspa(3,mxface,mxatom) 395 integer ijkspa( mxface,mxatom) 396 integer numpps( mxface,mxatom) 397 double precision xyzff( mxface,mxatom) 398 double precision zero, one 399 data one /1.0d+00/ 400 integer l_efcc, k_efcc, l_efcs, k_efcs, l_efcz, k_efcz 401 integer l_efclb, k_efclb, k_efciat, l_efciat 402 double precision ratm_real,dij,dum,cavdsp,pi,zetai,zetaii 403 integer m,mfac,mseg 404 integer nefc,iat,jat,npp,i,iseg,ifac,ief,ipp 405c 406 double precision cosff 407 external cosff 408c 409c MN solvation models --> 410c 411 logical do_cosmo_smd 412c 413c <-- MN solvation models 414c 415 double precision dist, xi, yi, zi, xj, yj, zj, rin, rout, alphai 416 parameter (alphai = 0.5d0) 417 dist(xi,yi,zi,xj,yj,zj)=sqrt((xj-xi)**2+(yj-yi)**2+(zj-zi)**2) 418 rin(iat)=ratm(iat)*(1.0d0-alphai*gammas*sqrt(0.25d0**minbem)) 419 rout(iat)=ratm(iat)*(1.0d0+(1.0d0-alphai)*gammas* 420 & sqrt(0.25d0**minbem)) 421c 422 dbug=.false. 423 pi = acos(-1.0d0) 424c 425 if(ga_nodeid().eq.0) then 426 write(luout,9999) 427 endif 428c 429c ----- print atomic centers ----- 430c 431 if(ga_nodeid().eq.0) then 432 write(luout,9998) 433 do iat=1,nat 434 if (do_cosmo_model.eq.DO_COSMO_KS) then 435 write(luout,9997) iat,xyzatm(1,iat),xyzatm(2,iat), 436 1 xyzatm(3,iat), 437 2 (ratm(iat)*cau2ang-rsolv) 438 else if (do_cosmo_model.eq.DO_COSMO_YK) then 439 write(luout,9997) iat,xyzatm(1,iat),xyzatm(2,iat), 440 1 xyzatm(3,iat), 441 2 (ratm(iat)*cau2ang) 442 endif 443 enddo 444 endif 445c 446c ----- clear arrays ..... ----- 447c 448 do iat=1,nat 449 do i=1,mxface 450 ijkspa(i,iat)=0 451 numpps(i,iat)=0 452 xyzff(i,iat)=0d0 453 enddo 454 enddo 455c 456c ----- sift through atomic centers and decide if a face ----- 457c belongs to the -sas- or is inside the molecule. 458c 459 do iat=1,nat 460c 461 if(ratm(iat).ne.0d0) then 462 do iseg=1,nseg 463 ijkspa(iseg,iat)=ijkseg(iseg) 464 xyzff(iseg,iat)=one 465 do m=1,3 466 xyzspa(m,iseg,iat)=xyzseg(m,iseg)*ratm(iat) 467 1 +xyzatm(m,iat) 468 enddo 469 enddo 470 do ifac=1,nfac 471 ijkspa(ifac+nseg,iat)=ijkseg(ifac+nseg) 472 do m=1,3 473 xyzspa(m,ifac+nseg,iat)=xyzseg(m,ifac+nseg)*ratm(iat) 474 1 +xyzatm(m,iat) 475 enddo 476 enddo 477 if(dbug.and.ga_nodeid().eq.0) then 478 write(luout,9996) iat 479 write(luout,9995) (ijkspa(ifac+nseg,iat),ifac=1,nfac) 480 endif 481 do jat=1,nat 482 dij=dist(xyzatm(1,iat), 483 1 xyzatm(2,iat), 484 2 xyzatm(3,iat), 485 3 xyzatm(1,jat), 486 4 xyzatm(2,jat), 487 5 xyzatm(3,jat)) 488 if (do_cosmo_model.eq.DO_COSMO_KS) then 489 if(jat.ne.iat.and.(dij.lt.(ratm(iat)+ratm(jat)))) then 490 do ifac=1,nfac 491 dum=dist(xyzspa(1,ifac+nseg,iat), 492 1 xyzspa(2,ifac+nseg,iat), 493 2 xyzspa(3,ifac+nseg,iat), 494 3 xyzatm(1,jat), 495 4 xyzatm(2,jat), 496 5 xyzatm(3,jat)) 497 if(dum.lt.ratm(jat)) then 498 ijkspa(ifac+nseg,iat)=0 499 endif 500 enddo 501 endif 502 else if (do_cosmo_model.eq.DO_COSMO_YK) then 503 if((jat.ne.iat).and.(ratm(jat).ne.0d0) 504 1 .and.(dij.lt.(ratm(iat)+rout(jat)))) then 505 do iseg=1,nseg 506 dum=dist(xyzspa(1,iseg,iat), 507 1 xyzspa(2,iseg,iat), 508 2 xyzspa(3,iseg,iat), 509 3 xyzatm(1,jat), 510 4 xyzatm(2,jat), 511 5 xyzatm(3,jat)) 512 xyzff(iseg,iat) = xyzff(iseg,iat) * 513 1 cosff((dum-rin(jat))/(rout(jat)-rin(jat))) 514 enddo 515 endif 516 endif 517 enddo 518 if(dbug.and.ga_nodeid().eq.0) then 519 write(luout,9996) iat 520 write(luout,9995) (ijkspa(ifac+nseg,iat),ifac=1,nfac) 521 endif 522c 523c ----- check segments belonging to -sas- ----- 524c 525 if (do_cosmo_model.eq.DO_COSMO_KS) then 526 do ifac=1,nseg+nfac 527 insseg(ifac,iat)=.true. 528 enddo 529 do ifac=1,nfac 530 iseg=ijkspa(ifac+nseg,iat) 531 if(iseg.ne.0) then 532 insseg(ifac+nseg,iat)=.false. 533 insseg( iseg,iat)=.false. 534 endif 535 enddo 536 else if (do_cosmo_model.eq.DO_COSMO_YK) then 537 do iseg=1,nseg 538 insseg(iseg,iat)=.not.(xyzff(iseg,iat).ge.swtol) 539 enddo 540 endif 541 if(dbug.and.ga_nodeid().eq.0) then 542 write(luout,9994) iat 543 if (do_cosmo_model.eq.DO_COSMO_KS) then 544 write(luout,9993) (insseg(ifac,iat),ifac=1,nseg+nfac) 545 else if (do_cosmo_model.eq.DO_COSMO_YK) then 546 write(luout,9993) (insseg(iseg,iat),iseg=1,nseg) 547 endif 548 endif 549 mseg=0 550 do iseg=1,nseg 551 if(.not.insseg(iseg,iat)) mseg=mseg+1 552 enddo 553 mfac=0 554 if (do_cosmo_model.eq.DO_COSMO_KS) then 555 do ifac=1,nfac 556 if(.not.insseg(ifac+nseg,iat)) mfac=mfac+1 557 enddo 558 endif 559 nspa(iat)=mseg 560 nppa(iat)=mfac 561c 562c ----- surface area of segments ----- 563c 564 if (do_cosmo_model.eq.DO_COSMO_KS) then 565 do iseg=1,nseg 566 numpps(iseg,iat)=0 567 enddo 568 do ifac=1,nfac 569 iseg=ijkspa(ifac+nseg,iat) 570 if(iseg.ne.0) numpps(iseg,iat)=numpps(iseg,iat)+1 571 enddo 572 endif 573c 574 endif 575c 576 enddo 577c 578 if(ga_nodeid().eq.0) then 579 write(luout,9985) nseg,nfac 580 write(luout,9992) 581 do iat=1,nat 582 npp=0 583 do iseg=1,nseg 584 npp=npp+numpps(iseg,iat) 585 enddo 586 write(luout,9991) iat,nspa(iat),nppa(iat),npp 587 enddo 588 endif 589 if(dbug.and.ga_nodeid().eq.0) then 590 write(luout,9987) 591 do iat=1,nat 592 do iseg=1,nseg 593 write(luout,9986) iat,iseg,numpps(iseg,iat) 594 enddo 595 enddo 596 endif 597c 598c Count the number of surface points, i.e. number of point charges 599c and generate memory to store them 600c 601 nefc = 0 602 do iat=1,nat 603 if(ratm(iat).ne.0d0) then 604 do iseg=1,nseg 605 if(.not.insseg(iseg,iat)) nefc = nefc+1 606 enddo 607 endif 608 enddo 609c 610c Allocate memory for point charges 611c 612 if(.not.ma_push_get(mt_dbl,nefc*3,'cosmo efcc',l_efcc,k_efcc)) 613 & call errquit('cosmo_cossas malloc k_efcc failed',911,MA_ERR) 614 if(.not.ma_push_get(mt_dbl,nefc,'cosmo efcs',l_efcs,k_efcs)) 615 & call errquit('cosmo_cossas malloc k_efcs failed',911,MA_ERR) 616 if(.not.ma_push_get(mt_int,nefc,'cosmo efciat',l_efciat,k_efciat)) 617 & call errquit('cosmo_cossas malloc k_efciat failed',911,MA_ERR) 618 if(.not.ma_push_get(mt_dbl,nefc,'cosmo efcz',l_efcz,k_efcz)) 619 & call errquit('cosmo_cossas malloc k_efcz failed',911,MA_ERR) 620 if(.not.ma_push_get(mt_byte,nefc*8,'cosmo tags',l_efclb,k_efclb)) 621 & call errquit('cosmo_cossas malloc k_tag failed',911,MA_ERR) 622c 623c ----- save coordinates of surface points ----- 624c save segment surfaces 625c save segment to atom mapping 626c 627 srfmol=0d0 628 volmol=0d0 629 ief =0 630 do iat=1,nat 631 if(ratm(iat).ne.0d0) then 632 ratm_real=-1d99 633 if (do_cosmo_model.eq.DO_COSMO_KS) then 634 ratm_real=ratm(iat)-rsolv/cau2ang 635 else if (do_cosmo_model.eq.DO_COSMO_YK) then 636 ratm_real=ratm(iat) 637 endif 638 do iseg=1,nseg 639 if(.not.insseg(iseg,iat)) then 640 ief=ief+1 641 do i=1,3 642 dbl_mb(k_efcc+3*(ief-1)+i-1)=xyzatm(i,iat) 643 1 +xyzseg(i,iseg)*ratm_real 644 enddo 645 ipp=numpps(iseg,iat) 646 if (do_cosmo_model.eq.DO_COSMO_KS) then 647 dbl_mb(k_efcs+ief-1) = dble(ipp)*dsurf*ratm_real**2 648 srfmol = srfmol + dble(ipp)*dsurf*ratm_real**2 649 volmol = volmol + dble(ipp)*dvol *ratm_real**3 650 else if (do_cosmo_model.eq.DO_COSMO_YK) then 651c 652c --- eval eq.(67) from [2] --- 653c 654 dum=4.00d0**(maxbem-minbem) 655 dum=16.0d0 ! MAXBEM is obsolete in York and Karplus approach 656 zetai=zeta*sqrt(nseg*dum)/(ratm_real*sqrt(2.0d0*pi)) 657 zetaii=zetai/sqrt(2.0d0) 658 dbl_mb(k_efcs+ief-1) = (1.0d0/xyzff(iseg,iat)) 659 1 * 2.0d0*zetaii/sqrt(pi) 660 srfmol = srfmol + xyzff(iseg,iat)*dsurf*ratm_real**2 661 volmol = volmol + xyzff(iseg,iat)*dvol *ratm_real**3 662 endif 663 int_mb(k_efciat+ief-1)=iat 664 endif 665 enddo 666 endif 667 enddo 668 srfmol=srfmol*(cau2ang**2) 669 volmol=volmol*(cau2ang**3) 670c 671 if(ga_nodeid().eq.0) then 672 write(luout,9990) nefc 673 write(luout,9984) srfmol 674 write(luout,9983) volmol 675 endif 676c 677c ----- Cavity/Dispersion free energy --- 678c Sitkoff, Sharp, and Honig, 679c J.Phys.Chem. 98, 1978 (1994) 680c 681 cavdsp=0.860+0.005*srfmol 682c 683c MN solvation models --> 684c 685c if(ga_nodeid().eq.0) then 686c write(luout,9981) cavdsp 687c endif 688 if (.not. 689 $ rtdb_get(rtdb,'cosmo:do_cosmo_smd',mt_log,1,do_cosmo_smd)) 690 $ call errquit('hnd_cossas: cannot get do_cosmo_smd from rtdb', 691 $ 0,rtdb_err) 692 if(ga_nodeid().eq.0) then 693 if (.not.do_cosmo_smd) write(luout,9981) cavdsp 694 endif 695c 696c <-- MN solvation models 697c 698c ----- print segment surfaces ----- 699c 700 if(dbug.and.ga_nodeid().eq.0) then 701 write(luout,9989) 702 do ief=1,nefc 703 write(luout,9988) ief,dbl_mb(k_efcs+ief-1), 704 & int_mb(k_efciat+ief-1) 705 enddo 706 endif 707c 708 do ief=1,nefc 709 dbl_mb(k_efcz+ief-1)=0d0 710 enddo 711 do ief=1,nefc 712 byte_mb(k_efclb+(ief-1)*8)=' ' 713 enddo 714c 715c ----- write out to -rtdb- ----- 716c 717 if(.not.rtdb_put(rtdb,'cosmo:nefc',mt_int,1 ,nefc)) 718 $ call errquit('hnd_cossas: rtdb put failed for nefc ',911, 719 & rtdb_err) 720 if(.not.rtdb_put(rtdb,'cosmo:efcc',mt_dbl,3*nefc,dbl_mb(k_efcc))) 721 $ call errquit('hnd_cossas: rtdb put failed for efcc ',912, 722 & rtdb_err) 723 if(.not.rtdb_put(rtdb,'cosmo:efcz',mt_dbl, nefc,dbl_mb(k_efcz))) 724 $ call errquit('hnd_cossas: rtdb put failed for efcz ',913, 725 & rtdb_err) 726 if(.not.rtdb_put(rtdb,'cosmo:efcs',mt_dbl, nefc,dbl_mb(k_efcs))) 727 $ call errquit('hnd_cossas: rtdb put failed for efcs ',914, 728 & rtdb_err) 729c 730c ----- reset cosmo:rawt to avoid trouble in cosmo charge 731c calculation ----- 732c 733 if(.not.rtdb_put(rtdb,'cosmo:rawt',mt_dbl, nefc,dbl_mb(k_efcz))) 734 $ call errquit('hnd_cossas: rtdb put failed for rawt ',915, 735 & rtdb_err) 736c 737c We will need the next bit of information to calculate the analytic 738c COSMO gradients. This table describes the relationship between 739c the COSMO charges and the associated atoms. So we better save this 740c info. 741c 742 if(.not.rtdb_put(rtdb,'cosmo:efciat',mt_int,nefc, 743 $ int_mb(k_efciat))) 744 $ call errquit('hnd_cossas: rtdb put failed for iatefc',916, 745 & rtdb_err) 746c if(.not.rtdb_cput(rtdb,'char variable',nefc,byte_mb(k_efclb))) 747c $ call errquit('hnd_cossas: rtdb put failed for efclab',917, 748c & rtdb_err) 749c 750 if(.not.ma_pop_stack(l_efclb)) call 751 & errquit('cosmo_cossas chop stack k_efclb failed',911,MA_ERR) 752 if(.not.ma_pop_stack(l_efcz)) call 753 & errquit('cosmo_cossas chop stack k_efcz failed',911,MA_ERR) 754 if(.not.ma_pop_stack(l_efciat)) call 755 & errquit('cosmo_cossas chop stack k_efciat failed',911,MA_ERR) 756 if(.not.ma_pop_stack(l_efcs)) call 757 & errquit('cosmo_cossas chop stack k_efcs failed',911,MA_ERR) 758 if(.not.ma_pop_stack(l_efcc)) call 759 & errquit('cosmo_cossas chop stack k_efcc failed',911,MA_ERR) 760c 761 return 762 9999 format(/,1x,'solvent accessible surface',/,1x,26(1h-)) 763 9998 format(/,1x,'---------- ATOMIC COORDINATES (A.U.) ----------', 764 1 '-- VDWR(ANG.) --') 765 9997 format( 1x,i5,3f14.8,f10.3) 766 9996 format(/,1x,'---------- SEGMENTS FOR -IAT- = ',i5) 767 9995 format(16i4) 768 9994 format(/,1x,'-INSSEG- FACES FOR IAT = ',i5) 769 9993 format(16l4) 770 9992 format( 1x,'atom',' ( ',' nspa',',',' nppa',' )',/,1x,22(1h-)) 771 9991 format( 1x, i4 ,' ( ', i6 ,',', i6 ,' )',i8) 772 9990 format( 1x,'number of -cosmo- surface points = ',i8) 773 9989 format( 1x,'SEGMENT SURFACES =',/,1x,18(1h-)) 774 9988 format(i8,f10.5,i5) 775 9987 format( 1x,'NUMBER OF FACES / SEGMENT =',/,1x,27(1h-)) 776 9986 format(3i5) 777 9985 format(' number of segments per atom = ',i10,/, 778 1 ' number of points per atom = ',i10) 779 9984 format(' molecular surface = ',f10.3,' angstrom**2') 780 9983 format(' molecular volume = ',f10.3,' angstrom**3') 781 9981 format(' G(cav/disp) = ',f10.3,' kcal/mol') 782 end 783c 784C> \brief Triangulate a sphere using the Boundary Element Method (BEM) 785C> 786C> This routine approximates a sphere starting from either an 787C> octahedron or an icosahedron and partitioning the triangles that 788C> make up these polyhedra. Each triangle is partitioned into four 789C> triangles at each level in the recursion. The procedure is starting 790C> from an equal sided triangle, select the midpoints of all three 791C> sides, and draw a triangle through the three midpoints. This way four 792C> triangles of equal size are obtained. However, the midpoints of the 793C> original sides do not lie on the surface of the sphere and therefore 794C> they need to be projected outwards. This outwards projection changes 795C> the size of the central triangle more than that of the other three. 796C> So in the final triangulation the triangles are not all of the same 797C> size, but this is ignored in the COSMO formalism. 798C> 799C> Ultimately we are interested only in the triangles at 2 levels of 800C> granularity: 801C> 802C> - minbem: these triangles are referred to as "segments" and 803C> represent the sphere and their centers become the positions for 804C> the COSMO charges. 805C> 806C> - maxbem: these triangles are referred to as "faces" and they are 807C> used to adjust the surface of the segments in regions where two 808C> atomic spheres meet and the segments straddle the boundary between 809C> both spheres. 810C> 811C> All other triangles are reduced to mere artefacts of the triangle 812C> generation algorithm. The array `ijkseg` administrates what the 813C> status of a triangle is. It lists for each face which segment it is 814C> part of. 815C> 816C> In addition this routine computes `adiag` of [1] Eq.(B1). The 817C> expression in this routine can be obtained from Eq.(B1) as 818C> \f{eqnarray*}{ 819C> \frac{1}{2R} 820C> &=&\frac{M}{2}\sum_{\nu=2}^M\frac{M^{-2}}{||t_1-t_\nu||} 821C> + MM^{-2}\frac{a_{\mathrm{diag}}}{2} \\\\ 822C> \frac{a_{\mathrm{diag}}}{M} 823C> &=&\frac{1}{R}-\sum_{\nu=2}^M\frac{M^{-1}}{||t_1-t_\nu||} \\\\ 824C> a_{\mathrm{diag}} 825C> &=&\frac{M}{R}-\sum_{\nu=2}^M\frac{1}{||t_1-t_\nu||} 826C> \f} 827C> The expression implemented in this routine can be mapped onto this 828C> by 829C> \f{eqnarray*}{ 830C> a_{\mathrm{diag}}' &=& \left(\frac{4\pi}{M}\right)^{1/2} 831C> \left(M-\sum_{\nu=2}\frac{1}{||t_1'-t_\nu'||}\right) 832C> \f} 833C> where \f$a_{\mathrm{diag}}'\f$ is `adiag` as calculated in this 834C> routine, the \f$t'\f$ are points in the unit sphere (as opposed to 835C> \f$t\f$ which are points on the sphere with radius \f$R\f$). 836C> In the routines `hnd_coschg`, `hnd_cosaxd` and `hnd_cosxad` this is 837C> multiplied with 838C> \f$|S_\mu|^{-1/2} = \left(\frac{4\pi R^2}{M}\right)^{-1/2}\f$ to give 839C> the proper \f$a_{\mathrm{diag}}\f$ 840C> \f{eqnarray*}{ 841C> a_{\mathrm{diag}} &=& \left(\frac{4\pi}{M}\right)^{1/2} 842C> \left(M-\sum_{\nu=2}\frac{1}{||t_1'-t_\nu'||}\right) 843C> \left(\frac{M}{4\pi R^2}\right)^{1/2} \\\\ 844C> &=& \frac{1}{R} 845C> \left(M-\sum_{\nu=2}\frac{1}{||t_1'-t_\nu'||}\right) \\\\ 846C> &=& \left(\frac{M}{R}-\sum_{\nu=2}\frac{1}{||t_1-t_\nu||}\right) 847C> \f} 848C> In ref.[1] Eq.(B2) is wrong because of a spurious scale factor 849C> \f$4\pi R^2/M\f$. 850C> 851C> ### References ### 852C> 853C> [1] A. Klamt, G. Schüürmann, 854C> "COSMO: a new approach to dielectric screening in solvents with 855C> explicit expressions for the screening energy and its gradient", 856C> <i>J. Chem. Soc., Perkin Trans. 2</i>, 1993, pp 799-805, DOI: 857C> <a href="https://doi.org/10.1039/P29930000799"> 858C> 10.1039/P29930000799</a>. 859C> 860 subroutine hnd_cossph(nseg,nfac,ndiv, 861 1 ijkfac,xyzseg,ijkseg,mxface,apex,mxapex, 862 2 dsurf_in,dvol_in,adiag_in) 863 implicit double precision (a-h,o-z) 864#include "cosmo_params.fh" 865#include "global.fh" 866#include "stdio.fh" 867cnew 868#include "cosmoP.fh" 869c 870c ----- starting from -icosahedron- ----- 871c 872c pass, napex, nface, error = 0 12 20 20 873c pass, napex, nface, error = 1 42 80 100 0.4982 874c pass napex, nface, error = 2 162 320 420 0.1848 875c pass napex, nface, error = 3 642 1280 1700 0.0523 876c pass napex, nface, error = 4 2562 5120 6820 0.0135 877c pass napex, nface, error = 5 10242 20480 27300 0.0034 878c 879c ----- starting from -octahedron- ----- 880c 881c pass, napex, nface, error = 0 6 8 8 882c pass, napex, nface, error = 1 18 32 40 0.8075 883c pass napex, nface, error = 2 66 128 168 0.4557 884c pass napex, nface, error = 3 258 512 680 0.1619 885c pass napex, nface, error = 4 1026 2048 2728 0.0451 886c pass napex, nface, error = 5 4098 8192 10920 0.0116 887c pass napex, nface, error = 6 16386 32768 43688 0.0029 888c 889 dimension apex(3,*) 890 dimension ijkfac(3,*) 891 dimension ijkseg( *) 892 dimension xyzseg(3,*) 893 parameter (mxpass= 7) 894 dimension minfac(mxpass) 895 dimension maxfac(mxpass) 896 dimension minico(mxpass) 897 dimension maxico(mxpass) 898 dimension minoct(mxpass) 899 dimension maxoct(mxpass) 900 dimension ijknew(3) 901 dimension ijkold(3) 902 equivalence (ijkold(1),iold),(ijkold(2),jold),(ijkold(3),kold) 903 equivalence (ijknew(1),inew),(ijknew(2),jnew),(ijknew(3),knew) 904 logical icos 905 logical octa 906 logical some,out,dbug 907 data minico / 1, 21, 101, 421, 1701, 6821, 0/ 908 data maxico / 20, 100, 420, 1700, 6820,27300, 0/ 909 data minoct / 1, 9, 41, 169, 681, 2729,10921/ 910 data maxoct / 8, 40, 168, 680, 2728,10920,43688/ 911 data zero /0.0d+00/ 912 data one /1.0d+00/ 913 data two /2.0d+00/ 914 data three /3.0d+00/ 915 data four /4.0d+00/ 916c 917 dist(xi,yi,zi,xj,yj,zj)=sqrt((xj-xi)**2+(yj-yi)**2+(zj-zi)**2) 918c 919 dbug=.false. 920 out =.false. 921 out =out.or.dbug 922 some=.false. 923 some=some.or.out 924c 925 pi=four*atan(one) 926 rad=one 927 cir= two*pi*rad 928 srf=four*pi*rad**2 929 vol=four*pi*rad**3/three 930c 931 npass=maxbem 932c 933c ----- define hedron ----- 934c define -minfac- 935c define -maxfac- 936c 937 icos=ificos.ne.0 938 octa=.not.icos 939 if(icos) then 940 call hnd_sphico(apex,napex,ijkfac,ijkseg,nface) 941 do ipass=1,mxpass 942 minfac(ipass)=minico(ipass) 943 maxfac(ipass)=maxico(ipass) 944 enddo 945 endif 946 if(octa) then 947 call hnd_sphoct(apex,napex,ijkfac,ijkseg,nface) 948 do ipass=1,mxpass 949 minfac(ipass)=minoct(ipass) 950 maxfac(ipass)=maxoct(ipass) 951 enddo 952 endif 953 if(some.or.out.or.dbug.and.ga_nodeid().eq.0) then 954 if(icos) then 955 write(luout,9994) 956 endif 957 if(octa) then 958 write(luout,9982) 959 endif 960 if(out) then 961 write(luout,*) '-minbem- = ',minbem 962 write(luout,*) '-maxbem- = ',maxbem 963 write(luout,*) '-minfac- = ',minfac 964 write(luout,*) '-maxfac- = ',maxfac 965 write(luout,*) '-npass - = ',npass 966 write(luout,9999) 967 do iapex=1,napex 968 write(luout,9998) iapex,apex(1,iapex), 969 1 apex(2,iapex), 970 2 apex(3,iapex) 971 enddo 972 endif 973 endif 974c 975c ----- loop over divisions to create sphere ----- 976c 977 mxfac=0 978 ipass=1 979 100 ipass=ipass+1 980 mnfac=mxfac+1 981 mxfac=nface 982 if(out.and.ga_nodeid().eq.0) then 983 write(luout,9996) ipass,napex,nface,mnfac,mxfac 984 endif 985c 986 dmin =one 987 mapex=napex 988 mface=nface 989 do lface=mnfac,mxfac 990 iold=ijkfac(1,lface) 991 jold=ijkfac(2,lface) 992 kold=ijkfac(3,lface) 993 call hnd_sphapx(apex,mapex,ijkfac,ijkseg,mface,lface, 994 1 ijkold,ijknew,dijk) 995 dmin=min(dmin,dijk) 996 enddo 997 napex=mapex 998 nface=mface 999 if(out.and.ga_nodeid().eq.0) then 1000 write(luout,9995) napex,nface 1001 endif 1002c 1003c ----- print out new apeces ----- 1004c 1005 if(dbug.and.ga_nodeid().eq.0) then 1006 do iapex=1,napex 1007 write(luout,9998) iapex,apex(1,iapex),apex(2,iapex), 1008 1 apex(3,iapex) 1009 enddo 1010 endif 1011c 1012c ----- print approximate volume ----- 1013c 1014 radapp= dmin 1015 radrat= dmin 1016 raderr=one-radrat 1017 srfapp=srf*dmin**2 1018 srfrat= dmin**2 1019 srferr=one-srfrat 1020 volapp=vol*dmin**3 1021 volrat= dmin**3 1022 volerr=one-volrat 1023 if(out.and.ga_nodeid().eq.0) then 1024 write(luout,9997) vol,volapp,volrat,volerr 1025 write(luout,9992) srf,srfapp,srfrat,srferr 1026 write(luout,9991) rad,radapp,radrat,raderr 1027 endif 1028c 1029c ----- assign early segment to each face ----- 1030c 1031 if(ipass.gt.(minbem+1)) then 1032 if(dbug.and.ga_nodeid().eq.0) then 1033 write(luout,9981) ipass 1034 write(luout,9980) (minfac(i),i=1,ipass) 1035 write(luout,9979) (maxfac(i),i=1,ipass) 1036 endif 1037 maxseg=maxfac(minbem) 1038 lfacmn=minfac(ipass) 1039 lfacmx=maxfac(ipass) 1040 if(dbug.and.ga_nodeid().eq.0) then 1041 write(luout,9990) ipass 1042 write(luout,9988) (ijkseg(lface),lface=lfacmn,lfacmx) 1043 endif 1044 do lface=lfacmn,lfacmx 1045 ijkseg(lface)=ijkseg(ijkseg(lface)) 1046 if(ijkseg(lface).gt.maxseg.and.ga_nodeid().eq.0) then 1047 write(luout,9987) lface,ijkseg(lface) 1048 endif 1049 enddo 1050 if(dbug.and.ga_nodeid().eq.0) then 1051 write(luout,9989) ipass 1052 write(luout,9988) (ijkseg(lface),lface=lfacmn,lfacmx) 1053 endif 1054 endif 1055c 1056 if(ipass.lt.npass) go to 100 1057c 1058c ----- end of loop over tessalating passes ----- 1059c 1060 if(dbug.and.ga_nodeid().eq.0) then 1061 do ipass=1,npass 1062 lfacmn=minfac(ipass) 1063 lfacmx=maxfac(ipass) 1064 write(luout,9989) ipass 1065 write(luout,*) '-lfacmn- = ',lfacmn 1066 write(luout,*) '-lfacmx- = ',lfacmx 1067 write(luout,9988) (ijkseg(lface),lface=lfacmn,lfacmx) 1068 enddo 1069 endif 1070 if(some.or.out.or.dbug.and.ga_nodeid().eq.0) then 1071 write(luout,9993) npass,napex,minfac(npass),maxfac(npass), 1072 1 radapp,raderr,srfapp,srferr,volapp,volerr 1073 endif 1074c 1075c ----- at this point each of the faces is assigned to one ----- 1076c segment. now define centers of segments ... 1077c 1078 third =one/three 1079 lfacmn= minfac(minbem) 1080 lfacmx= maxfac(minbem) 1081 do lface=lfacmn,lfacmx 1082 mface=lface-lfacmn+1 1083 ijkseg(mface)=mface 1084 i=ijkfac(1,lface) 1085 j=ijkfac(2,lface) 1086 k=ijkfac(3,lface) 1087 do m=1,3 1088 xyzseg(m,mface)=(apex(m,i)+apex(m,j)+apex(m,k))*third 1089 enddo 1090 dseg=one/dist(xyzseg(1,mface),xyzseg(2,mface),xyzseg(3,mface), 1091 1 zero,zero,zero) 1092 do m=1,3 1093 xyzseg(m,mface)=xyzseg(m,mface)*dseg 1094 enddo 1095 enddo 1096 nseg=(lfacmx-lfacmn+1) 1097c 1098 if(dbug.and.ga_nodeid().eq.0) then 1099 lfacmn=1 1100 lfacmx=nseg 1101 write(luout,*) 'segment to segment mapping = ' 1102 write(luout,9988) (ijkseg(lface),lface=lfacmn,lfacmx) 1103 endif 1104c 1105c ----- now the faces ... ----- 1106c 1107 if(npass.gt.minbem) then 1108 lfacmn=minfac(minbem+1) 1109 lfacmx=maxfac(npass ) 1110 do lface=lfacmn,lfacmx 1111 mface=lface-lfacmn+1 1112 1 +(maxfac(minbem)-minfac(minbem)+1) 1113 ijkseg(mface)=ijkseg(lface) 1114 1 -( minfac(minbem)-1) 1115 i=ijkfac(1,lface) 1116 j=ijkfac(2,lface) 1117 k=ijkfac(3,lface) 1118 do m=1,3 1119 xyzseg(m,mface)=(apex(m,i)+apex(m,j)+apex(m,k))*third 1120 enddo 1121 dseg=one/dist(xyzseg(1,mface), 1122 1 xyzseg(2,mface), 1123 2 xyzseg(3,mface),zero,zero,zero) 1124 do m=1,3 1125 xyzseg(m,mface)=xyzseg(m,mface)*dseg 1126 enddo 1127 enddo 1128 nfac=(lfacmx-lfacmn+1) 1129c 1130c ----- only keep the faces at granularity maxbem ----- 1131c ----- discard all other faces ----- 1132c 1133c do lface=1,maxfac(maxbem-1)-minfac(minbem+1)+1 1134c ijkseg(nseg+lface) = 0 1135c enddo 1136 else 1137 do iseg=1,nseg 1138 ifac=iseg+nseg 1139 ijkseg(ifac)=ijkseg(iseg) 1140 do m=1,3 1141 xyzseg(m,ifac)=xyzseg(m,iseg) 1142 enddo 1143 enddo 1144 nfac=nseg 1145 endif 1146c 1147 if(dbug.and.ga_nodeid().eq.0) then 1148 lfacmn=nseg+1 1149 lfacmx=nseg+nfac 1150 write(luout,*) ' face to segment mapping = ' 1151 write(luout,9988) (ijkseg(lface),lface=lfacmn,lfacmx) 1152 endif 1153c 1154c ----- calculate -dsurf dvol- for the -cosmo- theory ----- 1155c 1156 if (do_cosmo_model.eq.DO_COSMO_YK) nfac =nseg 1157 ndiv =nfac/nseg 1158 dsurf_in=srf/dble(nfac) 1159 dvol_in =vol/dble(nfac) 1160 if(some.or.out.or.dbug.and.ga_nodeid().eq.0) then 1161 write(luout,9986) nseg,nfac,ndiv,dsurf_in,dvol_in 1162 endif 1163 if(out.and.ga_nodeid().eq.0) then 1164 write(luout,9985) 1165 do i=1,nseg 1166 done=dist(xyzseg(1,i),xyzseg(2,i),xyzseg(3,i), 1167 1 zero,zero,zero) 1168 write(luout,9984) i, 1169 1 xyzseg(1,i),xyzseg(2,i),xyzseg(3,i), 1170 2 ijkseg(i),done 1171 enddo 1172 endif 1173 if(dbug.and.ga_nodeid().eq.0) then 1174 write(luout,9985) 1175 do i=nseg+1,nseg+nfac 1176 done=dist(xyzseg(1,i),xyzseg(2,i),xyzseg(3,i), 1177 1 zero,zero,zero) 1178 write(luout,9984) (i-nseg), 1179 1 xyzseg(1,i),xyzseg(2,i),xyzseg(3,i), 1180 2 ijkseg(i),done 1181 enddo 1182 endif 1183c 1184c ----- calculate -adiag- of the -cosmo- theory ----- 1185c 1186 avgdia=zero 1187 avgfac=zero 1188 do mseg=1,nseg 1189 sum=zero 1190 do lseg=1,nseg 1191 if(lseg.ne.mseg) then 1192 l1=mseg 1193 l2=lseg 1194 sum=sum+rad/dist(xyzseg(1,l2),xyzseg(2,l2),xyzseg(3,l2), 1195 1 xyzseg(1,l1),xyzseg(2,l1),xyzseg(3,l1)) 1196 endif 1197 enddo 1198 fac=(dble(nseg)-sum)/sqrt(dble(nseg)) 1199 adiag_in=sqrt(four*pi)*fac 1200 if(some.or.out.or.dbug.and.ga_nodeid().eq.0) then 1201 write(luout,9983) mseg,adiag_in,fac,dble(nseg),sum 1202 endif 1203 avgdia=avgdia+adiag_in 1204 avgfac=avgfac+fac 1205 enddo 1206 adiag_in=avgdia/dble(nseg) 1207 fac =avgfac/dble(nseg) 1208 if(some.or.out.or.dbug.and.ga_nodeid().eq.0) then 1209 write(luout,9978) adiag_in,fac 1210 endif 1211c 1212 return 1213 9999 format(/,1x,' apex',5x,'x',6x,5x,'y',6x,5x,'z',6x,/,1x,42(1h-)) 1214 9998 format(1x,i6,f12.8,f12.8,f12.8) 1215 9997 format(' vol, approx., ratio, error = ',2f12.8,2 f8.4) 1216 9996 format(' pass, napex, nface, mnfac, mxfac = ',i3,4i8) 1217 9995 format(' napex, nface = ',3x,2i8) 1218 9994 format(1x,'sphere from -icosahedron-',/,1x,25(1h-)) 1219 9993 format(' npass = ',i2,' napex = ',i8, 1220 1 ' minfac = ',i8,' maxfac = ',i8,/, 1221 2 ' rad = ',f10.6,' error = ',f8.4,/, 1222 3 ' srf = ',f10.6,' error = ',f8.4,/, 1223 4 ' vol = ',f10.6,' error = ',f8.4) 1224 9992 format(' srf, approx., ratio, error = ',2f12.8,2 f8.4) 1225 9991 format(' rad, approx., ratio, error = ',2f12.8,2 f8.4) 1226 9990 format(' absolute -ijkseg- , for -ipass- = ',i5) 1227 9989 format(' relative -ijkseg- , for -ipass- = ',i5) 1228 9988 format(12i6) 1229 9987 format(' assigned segment for -lface- = ',i7, 1230 1 ' is = ',i7,' ( greater than -maxseg- = ',i4,' )') 1231 9986 format(' nseg,nfac,ndiv=nfac/nseg,dsurf,dvol = ',3i7,2f10.6) 1232 9985 format(' pt ',' x ',' y ',' z ', 1233 1 ' seg ',' norm ',/,1x,59(1h-)) 1234 9984 format(i7,3f12.8,i5,f12.8) 1235 9983 format(' mseg,adiag,fac,m,sum = ',i7,4f12.6) 1236 9982 format(1x,'sphere from -octahedron-',/,1x,24(1h-)) 1237 9981 format(' pass # = ',i5) 1238 9980 format(' minfac = ',10i5) 1239 9979 format(' maxfac = ',10i5) 1240 9978 format(' adiag,fac = ', 2f12.6) 1241 end 1242C> 1243C> \brief Setup the data for an Octahedron 1244C> 1245C> This routine initiates the segments of an octahedron. The output 1246C> is split over a few arrays. One array `apex` holds the coordinates 1247C> of the corners of the triangles, another array `ijkfac` lists 1248C> for each triangle which points in the `apex` array hold the 1249C> corresponding corner points, and finally `ijkseg` record the 1250C> mapping between faces and segments. 1251C> 1252 subroutine hnd_sphoct(apex,napex,ijkfac,ijkseg,nface) 1253 implicit none 1254#include "global.fh" 1255#include "stdio.fh" 1256c 1257 logical out 1258 integer napex !< [Output] The number of apeces 1259 integer nface !< [Output] The number of faces 1260 double precision xyz(3,6) 1261 integer ijk(3,8) 1262 double precision apex(3,*) !< [Output] The corners of the 1263 !< triangles 1264 integer ijkfac(3,*) !< [Output] For each triangle which 1265 !< points make up its corners 1266 integer ijkseg( *) !< [Output] For each face the number 1267 !< of the segment it belongs to 1268 integer iapex, lface 1269 data xyz / 1.0d+00, 0.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 1270 1 -1.0d+00, 0.0d+00, 0.0d+00, 0.0d+00,-1.0d+00, 0.0d+00, 1271 2 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00,-1.0d+00/ 1272 data ijk / 5, 1, 2, 5, 2, 3, 5, 3, 4, 5, 4, 1, 1273 1 6, 1, 2, 6, 2, 3, 6, 3, 4, 6, 4, 1/ 1274c 1275 out=.false. 1276 out=out.and.ga_nodeid().eq.0 1277c 1278 if(out) then 1279 write(luout,9997) 1280 endif 1281c 1282c ----- set the 6 apeces of an octahedron ----- 1283c 1284c 1 1. 0. 0. 1285c 2 0. 1. 0. 1286c 3 -1. 0. 0. 1287c 4 0. -1. 0. 1288c 5 0. 0. 1. 1289c 6 0. 0. -1. 1290c 1291 napex=6 1292 do iapex=1,napex 1293 apex(1,iapex)=xyz(1,iapex) 1294 apex(2,iapex)=xyz(2,iapex) 1295 apex(3,iapex)=xyz(3,iapex) 1296 enddo 1297 if(out) then 1298 write(luout,9999) 1299 do iapex=1,napex 1300 write(luout,9998) iapex,apex(1,iapex),apex(2,iapex), 1301 1 apex(3,iapex) 1302 enddo 1303 endif 1304c 1305 nface=8 1306 do lface=1,nface 1307 ijkfac(1,lface)=ijk(1,lface) 1308 ijkfac(2,lface)=ijk(2,lface) 1309 ijkfac(3,lface)=ijk(3,lface) 1310 ijkseg( lface)= lface 1311 enddo 1312c 1313 if(out) then 1314 write(luout,*) '...... end of -sphoct- ......' 1315 endif 1316 return 1317 9999 format(/,1x,' apex',5x,'x',6x,5x,'y',6x,5x,'z',6x,/,1x,42(1h-)) 1318 9998 format(1x,i6,f12.8,f12.8,f12.8) 1319 9997 format(/,1x,'octahedron',/,1x,10(1h-)) 1320 end 1321c 1322 subroutine hnd_sphico(apex,napex,ijkfac,ijkseg,nface) 1323 implicit double precision (a-h,o-z) 1324#include "global.fh" 1325#include "stdio.fh" 1326c 1327 logical out 1328 dimension c(3,12) 1329 dimension s(3,12) 1330 dimension ijk(3,20) 1331 dimension apex(3,*) 1332 dimension ijkfac(3,*) 1333 dimension ijkseg( *) 1334 data c / 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 1335 1 0.0d+00,-1.0d+00, 0.0d+00, 0.0d+00,-1.0d+00, 0.0d+00, 1336 2 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 1337 3 0.0d+00, 0.0d+00,-1.0d+00, 0.0d+00, 0.0d+00,-1.0d+00, 1338 4 1.0d+00, 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00, 1339 5 -1.0d+00, 0.0d+00, 0.0d+00,-1.0d+00, 0.0d+00, 0.0d+00/ 1340 data s / 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00,-1.0d+00, 1341 1 0.0d+00, 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00,-1.0d+00, 1342 2 1.0d+00, 0.0d+00, 0.0d+00,-1.0d+00, 0.0d+00, 0.0d+00, 1343 3 1.0d+00, 0.0d+00, 0.0d+00,-1.0d+00, 0.0d+00, 0.0d+00, 1344 4 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00,-1.0d+00, 0.0d+00, 1345 5 0.0d+00, 1.0d+00, 0.0d+00, 0.0d+00,-1.0d+00, 0.0d+00/ 1346 data ijk / 1, 2, 9, 1, 9, 5, 1, 5, 6, 1, 6,11, 1,11, 2, 1347 1 2, 9, 7, 2, 7, 8, 2, 8,11, 1348 2 3, 4,10, 3,10, 5, 3, 5, 6, 3, 6,12, 3,12, 4, 1349 3 4,10, 7, 4, 7, 8, 4, 8,12, 1350 4 9,10, 7, 9, 5,10,11,12, 8,11, 6,12/ 1351 data one /1.0d+00/ 1352 data two /2.0d+00/ 1353 data five /5.0d+00/ 1354c 1355 out=.false. 1356 out=out.and.ga_nodeid().eq.0 1357c 1358 if(out) then 1359 write(luout,9997) 1360 endif 1361c 1362c ----- set the 12 apeces of an icosahedron ----- 1363c 1364c 1 0. cosa sina 1365c 2 0. cosa -sina 1366c 3 0. -cosa sina 1367c 4 0. -cosa -sina 1368c 5 sina 0. cosa 1369c 6 -sina 0. cosa 1370c 7 sina 0. -cosa 1371c 8 -sina 0. -cosa 1372c 9 cosa sina 0. 1373c 10 cosa -sina 0. 1374c 11 -cosa sina 0. 1375c 12 -cosa -sina 0. 1376c 1377 ang=acos(one/sqrt(five))/two 1378 cosa=cos(ang) 1379 sina=sin(ang) 1380 napex=12 1381 do iapex=1,napex 1382 apex(1,iapex)=cosa*c(1,iapex)+sina*s(1,iapex) 1383 apex(2,iapex)=cosa*c(2,iapex)+sina*s(2,iapex) 1384 apex(3,iapex)=cosa*c(3,iapex)+sina*s(3,iapex) 1385 enddo 1386 if(out) then 1387 write(luout,9999) 1388 do iapex=1,napex 1389 write(luout,9998) iapex,apex(1,iapex),apex(2,iapex), 1390 1 apex(3,iapex) 1391 enddo 1392 endif 1393c 1394 nface=20 1395 do lface=1,nface 1396 ijkfac(1,lface)=ijk(1,lface) 1397 ijkfac(2,lface)=ijk(2,lface) 1398 ijkfac(3,lface)=ijk(3,lface) 1399 ijkseg( lface)= lface 1400 enddo 1401c 1402 if(out) then 1403 write(luout,*) '...... end of -sphico- ......' 1404 endif 1405 return 1406 9999 format(/,1x,' apex',5x,'x',6x,5x,'y',6x,5x,'z',6x,/,1x,42(1h-)) 1407 9998 format(1x,i6,f12.8,f12.8,f12.8) 1408 9997 format(/,1x,'icosahedron',/,1x,11(1h-)) 1409 end 1410C> 1411C> \brief Partition a given triangle into four triangles and project 1412C> them outward onto the unit sphere 1413C> 1414 subroutine hnd_sphapx(apex,mapex,ijkfac,ijkseg,mface,lface, 1415 1 ijkold,ijknew,dmin) 1416 implicit double precision (a-h,o-z) 1417#include "global.fh" 1418#include "stdio.fh" 1419c 1420 logical out 1421 logical duplic 1422 dimension apex(3,*) 1423 dimension ijkfac(3,*) 1424 dimension ijkseg( *) 1425 dimension ijkold(3) 1426 dimension ijknew(3) 1427 dimension xyz(3,3) 1428 dimension d(3) 1429 dimension xyzijk(3) 1430 equivalence (xyz(1,1),xij),(xyz(2,1),yij),(xyz(3,1),zij), 1431 1 (xyz(1,2),xjk),(xyz(2,2),yjk),(xyz(3,2),zjk), 1432 2 (xyz(1,3),xki),(xyz(2,3),yki),(xyz(3,3),zki) 1433 data zero /0.0d+00/ 1434 data pt5 /0.5d+00/ 1435 data one /1.0d+00/ 1436 data three /3.0d+00/ 1437 data tol /1.0d-04/ 1438c 1439 dist(x1,y1,z1,x2,y2,z2)=sqrt((x2-x1)**2+(y2-y1)**2+(z2-z1)**2) 1440c 1441 out=.false. 1442 out=out.and.ga_nodeid().eq.0 1443c 1444c ----- create mid-point of the 3 edges ----- 1445c 1446 iold=ijkold(1) 1447 jold=ijkold(2) 1448 kold=ijkold(3) 1449 do m=1,3 1450 xyz(m,1)=(apex(m,iold)+apex(m,jold))*pt5 1451 xyz(m,2)=(apex(m,jold)+apex(m,kold))*pt5 1452 xyz(m,3)=(apex(m,kold)+apex(m,iold))*pt5 1453 enddo 1454c 1455c ----- project onto sphere ----- 1456c 1457 d(1)=dist(xij,yij,zij,zero,zero,zero) 1458 d(2)=dist(xjk,yjk,zjk,zero,zero,zero) 1459 d(3)=dist(xki,yki,zki,zero,zero,zero) 1460 d(1)=one/d(1) 1461 d(2)=one/d(2) 1462 d(3)=one/d(3) 1463 do l=1,3 1464 do m=1,3 1465 xyz(m,l)=xyz(m,l)*d(l) 1466 enddo 1467 enddo 1468c 1469c ----- check for duplicate apeces ----- 1470c 1471 newapx=0 1472 do iapx=1,3 1473 duplic=.false. 1474 lduplc=0 1475 do lapex=1,mapex 1476 dd=dist(xyz(1, iapx),xyz(2, iapx),xyz(3, iapx), 1477 1 apex(1,lapex),apex(2,lapex),apex(3,lapex)) 1478 if(abs(dd).lt.tol) then 1479 duplic=.true. 1480 lduplc=lapex 1481 endif 1482 enddo 1483 if(duplic) then 1484 ijknew(iapx)=lduplc 1485 if(out) then 1486 write(luout,9999) iapx,ijkold,lduplc 1487 endif 1488 else 1489 newapx=newapx+1 1490 japx=mapex+newapx 1491 ijknew(iapx)=japx 1492 do m=1,3 1493 apex(m,japx)=xyz(m,iapx) 1494 enddo 1495 if(out) then 1496 write(luout,9998) iapx,ijkold,japx, 1497 1 apex(1,japx),apex(2,japx),apex(3,japx) 1498 endif 1499 endif 1500 enddo 1501 mapex=mapex+newapx 1502c 1503c ----- make up new faces and their centers ----- 1504c 1505 third=one/three 1506 dmin =one 1507c 1508 mface=mface+1 1509 ijkseg( mface)=lface 1510 ijkfac(1,mface)=ijkold(1) 1511 ijkfac(2,mface)=ijknew(1) 1512 ijkfac(3,mface)=ijknew(3) 1513 do m=1,3 1514 xyzijk(m)=(apex(m,iold)+apex(m,jold)+apex(m,kold))*third 1515 enddo 1516 dijk=dist(xyzijk(1),xyzijk(2),xyzijk(3),zero,zero,zero) 1517 dmin=min(dmin,dijk) 1518c 1519 mface=mface+1 1520 ijkseg( mface)=lface 1521 ijkfac(1,mface)=ijkold(2) 1522 ijkfac(2,mface)=ijknew(1) 1523 ijkfac(3,mface)=ijknew(2) 1524 do m=1,3 1525 xyzijk(m)=(apex(m,iold)+apex(m,jold)+apex(m,kold))*third 1526 enddo 1527 dijk=dist(xyzijk(1),xyzijk(2),xyzijk(3),zero,zero,zero) 1528 dmin=min(dmin,dijk) 1529c 1530 mface=mface+1 1531 ijkseg( mface)=lface 1532 ijkfac(1,mface)=ijkold(3) 1533 ijkfac(2,mface)=ijknew(2) 1534 ijkfac(3,mface)=ijknew(3) 1535 do m=1,3 1536 xyzijk(m)=(apex(m,iold)+apex(m,jold)+apex(m,kold))*third 1537 enddo 1538 dijk=dist(xyzijk(1),xyzijk(2),xyzijk(3),zero,zero,zero) 1539 dmin=min(dmin,dijk) 1540c 1541 mface=mface+1 1542 ijkseg( mface)=lface 1543 ijkfac(1,mface)=ijknew(1) 1544 ijkfac(2,mface)=ijknew(2) 1545 ijkfac(3,mface)=ijknew(3) 1546 do m=1,3 1547 xyzijk(m)=(apex(m,iold)+apex(m,jold)+apex(m,kold))*third 1548 enddo 1549 dijk=dist(xyzijk(1),xyzijk(2),xyzijk(3),zero,zero,zero) 1550 dmin=min(dmin,dijk) 1551c 1552 if(out) then 1553 write(luout,9997) dmin,mface 1554 endif 1555c 1556 return 1557 9999 format(' duplicated apex =',i2,' for face ',3i5,'. same as = ',i5) 1558 9998 format(' new apex =',i2,' for face ',3i5,'. newapx = ',i5, 1559 1 /,7x,3f12.8) 1560 9997 format(' --- dmin = ',f12.8,' --- mface = ',i10) 1561 end 1562C> 1563C> \brief Evaluate COSMO related energy terms 1564C> 1565C> Based on the COSMO charges a number of energy terms can be evaluated, 1566C> including: 1567C> 1568C> - the nuclear - COSMO charge interaction energy 1569C> 1570C> - the electron - COSMO charge interaction energy 1571C> 1572C> - the total solute charge density - COSMO charge interaction energy 1573C> 1574C> - the COSMO charge - COSMO charge interaction energy 1575C> 1576C> These terms are inherent in Ref.[1] Eq.(8) and in Ref.[2] Eq.(39). 1577C> 1578C> ### References ### 1579C> 1580C> [1] A. Klamt, G. Schüürmann, 1581C> "COSMO: a new approach to dielectric screening in solvents with 1582C> explicit expressions for the screening energy and its gradient", 1583C> <i>J. Chem. Soc., Perkin Trans. 2</i>, 1993, pp 799-805, DOI: 1584C> <a href="https://doi.org/10.1039/P29930000799"> 1585C> 10.1039/P29930000799</a>. 1586C> 1587C> [2] D.M. York, M. Karplus, 1588C> "A smooth solvation potential based on the conductor-like 1589C> screening model", <i>J. Phys. Chem. A</i> (1999) <b>103</b>, 1590C> pp 11060-11079, DOI: 1591C> <a href="https://doi.org/10.1021/jp992097l"> 1592C> 10.1021/jp992097l</a>. 1593C> 1594 subroutine hnd_cos_energy(nat,nefc,chgscr,efcc,efcs,efcz,efciat, 1595 & ratm,catm,zatm,pot,allefc,atmefc,elcefc,efcefc) 1596 implicit none 1597c 1598#include "cosmoP.fh" 1599#include "cosmo_params.fh" 1600c 1601 integer nat !< [Input] The number of electrons 1602 integer nefc !< [Input] The number of COSMO charges 1603 integer efciat(nefc) !< [Input] Mapping of COSMO charges to the 1604 !< atom that generated the corresponding 1605 !< solvent accessible surface area segment 1606c 1607 double precision chgscr !< [Input] The dielectric screening 1608 !< factor \f$f(\epsilon)\f$ (see Ref.[2] 1609 !< Eq.(46) with \f$\epsilon_1=1\f$ and 1610 !< \f$\epsilon_2=\epsilon\f$) 1611 double precision efcc(3,nefc) !< [Input] The COSMO charge 1612 !< coordinates 1613 double precision efcs(nefc) !< [Input] The COSMO charge surface 1614 !< area \f$|S_\mu|\f$ in Ref.[1] 1615 !< (see e.g. Eq.(7b)) or Ref.[2] 1616 !< Eq.(67). 1617 double precision efcz(nefc) !< [Input] The COSMO charges 1618 double precision ratm(nat) !< [Input] The atomic radii 1619 double precision catm(3,nat) !< [Input] The atomic coordinates 1620 double precision zatm(nat) !< [Input] The nuclear charges 1621 double precision pot(nefc) !< [Input] The molecular potential 1622 !< at the COSMO charge positions 1623c 1624 double precision allefc !< [Output] The total solute charge 1625 !< density - COSMO charge interaction 1626 !< energy 1627 double precision atmefc !< [Output] The nuclear - COSMO charge 1628 !< interaction energy 1629 double precision elcefc !< [Output] The electron - COSMO charge 1630 !< interaction energy 1631 double precision efcefc !< [Output] The COSMO charge - COSMO 1632 !< charge interaction energy 1633c 1634c 1635 integer ief, jef !< Counters over COSMO charges 1636 integer iat !< Counter over atoms 1637c 1638 double precision xi, yi, zi, qi 1639 double precision xj, yj, zj, qj 1640 double precision aii, aij, bij, dij 1641c 1642 double precision zetai, zetaj, zetaij 1643c 1644 double precision zero, one, two 1645 parameter (zero=0.0d0, one=1.0d0, two=2.0d0) 1646 double precision pi 1647 double precision derf 1648c 1649 pi =acos(-1.0d0) 1650 allefc=zero 1651 atmefc=zero 1652 efcefc=zero 1653 do jef=1,nefc 1654 xj=efcc(1,jef) 1655 yj=efcc(2,jef) 1656 zj=efcc(3,jef) 1657 qj=efcz( jef) 1658c 1659 allefc=allefc+qj*pot(jef) 1660c 1661 do iat=1,nat 1662 xi=catm(1,iat) 1663 yi=catm(2,iat) 1664 zi=catm(3,iat) 1665 qi=zatm( iat) 1666 dij=sqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2) 1667 bij=one/dij 1668 atmefc=atmefc+qi*bij*qj 1669 enddo 1670 enddo 1671c 1672 if (do_cosmo_model.eq.DO_COSMO_KS) then 1673 do jef=1,nefc 1674 xj=efcc(1,jef) 1675 yj=efcc(2,jef) 1676 zj=efcc(3,jef) 1677 qj=efcz( jef) 1678 efcefc=efcefc+qj*adiag*qj/sqrt(efcs(jef)) 1679 do ief=jef+1,nefc 1680 qi=efcz( ief) 1681 xi=efcc(1,ief) 1682 yi=efcc(2,ief) 1683 zi=efcc(3,ief) 1684 dij=sqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2) 1685 aij=one/dij 1686 efcefc=efcefc+2*qi*aij*qj 1687 enddo 1688 enddo 1689 else if (do_cosmo_model.eq.DO_COSMO_YK) then 1690 do jef=1,nefc 1691 xj=efcc(1,jef) 1692 yj=efcc(2,jef) 1693 zj=efcc(3,jef) 1694 qj=efcz( jef) 1695 zetaj=zeta*sqrt(ptspatm)/(ratm(efciat(jef))*sqrt(2.0d0*pi)) 1696 efcefc=efcefc+qj*efcs(jef)*qj 1697c 1698 do ief=jef+1,nefc 1699 zetai=zeta*sqrt(ptspatm) 1700 & /(ratm(efciat(ief))*sqrt(2.0d0*pi)) 1701 xi=efcc(1,ief) 1702 yi=efcc(2,ief) 1703 zi=efcc(3,ief) 1704 qi=efcz( ief) 1705 dij=sqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2) 1706 zetaij=zetai*zetaj/sqrt(zetai**2+zetaj**2) 1707 if (dij.lt.1.0d-8) then 1708 aij=2.0d0*zetaij/sqrt(pi) 1709 else 1710 aij=derf(zetaij*dij)/dij 1711 endif 1712 efcefc=efcefc+2*qi*aij*qj 1713 enddo 1714 enddo 1715 endif 1716 efcefc= efcefc/(two*chgscr) 1717 elcefc= allefc-atmefc 1718c 1719 end 1720C> 1721C> \brief Compute matrix A 1722C> 1723C> Compute matrix `A` as needed for the in-memory COSMO 1724C> charge fitting. 1725C> 1726 subroutine hnd_cosmata(nat,nefc,efcc,efcs,efciat,ratm,a) 1727 implicit none 1728#include "cosmoP.fh" 1729#include "cosmo_params.fh" 1730 integer nat !< [Input] The number of atoms 1731 integer nefc !< [Input] The number of surface charges 1732 double precision efcc(3,nefc) !< [Input] The surface charge coords 1733 double precision efcs(nefc) !< [Input] The surface areas 1734 integer efciat(nat) !< [Input] The atom of a surface 1735 !< charge 1736 double precision ratm(nat) !< [Input] The atom radii 1737 double precision a(nefc,nefc) !< [Output] Matrix `A` 1738c 1739 integer jef, ief 1740 double precision aii, xi, xj, yi, yj, zi, zj, dij, one 1741 double precision zetai, zetaj, zetaij 1742 parameter (one = 1.0d0) 1743c 1744c 1745 double precision derf 1746 double precision factor, factor2 1747 1748 double precision pi 1749 pi = acos(-1.0d0) 1750c 1751 if (do_cosmo_model.eq.DO_COSMO_KS) then 1752 do jef=1,nefc 1753 xj=efcc(1,jef) 1754 yj=efcc(2,jef) 1755 zj=efcc(3,jef) 1756 do ief=1,nefc 1757 if(ief.eq.jef) then 1758 aii=adiag/sqrt(efcs(ief)) 1759 a(ief,jef)=aii 1760 else 1761 xi=efcc(1,ief) 1762 yi=efcc(2,ief) 1763 zi=efcc(3,ief) 1764 dij=sqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2) 1765 a(ief,jef)=one/dij 1766 endif 1767 enddo 1768 enddo 1769 else if (do_cosmo_model.eq.DO_COSMO_YK) then 1770 factor = zeta*dsqrt(ptspatm/(2d0*pi)) 1771 factor2 = 2d0/sqrt(pi) 1772 do jef=1,nefc 1773 zetaj=factor/ratm(efciat(jef)) 1774 xj=efcc(1,jef) 1775 yj=efcc(2,jef) 1776 zj=efcc(3,jef) 1777 do ief=1,nefc 1778 if(ief.eq.jef) then 1779 aii=efcs(ief) 1780 a(ief,jef)=aii 1781 else 1782 zetai=factor/ratm(efciat(ief)) 1783 zetaij=zetai*zetaj/sqrt(zetai**2+zetaj**2) 1784 xi=efcc(1,ief) 1785 yi=efcc(2,ief) 1786 zi=efcc(3,ief) 1787 dij=sqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2) 1788 if (dij.lt.1.0d-5) then 1789 a(ief,jef)=factor2*zetaij*(1d0 - (zetaij*dij)**2/3d0) 1790 else 1791 a(ief,jef)=derf(zetaij*dij)/dij 1792 endif 1793 endif 1794 enddo 1795 enddo 1796 endif 1797c 1798 end 1799C> 1800C> \brief On the fly matrix-vector multiplication 1801C> 1802C> This routine multiplies the COSMO matrix with the current guess 1803C> for the charge vector using a dot-product based algorithm. The matrix 1804C> is generated on the fly. For performance reasons the routine is 1805C> replicated data parallel. 1806C> 1807 subroutine hnd_cosaxd(nat,x,ax,nefc,efcc,efcs,efciat,ratm) 1808 implicit none 1809#include "global.fh" 1810#include "msgids.fh" 1811#include "cosmoP.fh" 1812#include "cosmo_params.fh" 1813c 1814 integer nat !< [Input] The number of atoms 1815 integer nefc !< [Input] The number of surface charges 1816 double precision efcc(3,nefc) !< [Input] The surface charge coords 1817 double precision efcs(nefc) !< [Input] The surface areas 1818 integer efciat(nefc) !< [Input] The atom of a surface 1819 !< charge 1820 double precision ratm(nat) !< [Input] The atom radii 1821 double precision x(nefc) !< [Input] Vector `x` 1822 double precision ax(nefc) !< [Output] Matrix-vector product 1823 !< `Ax` 1824c 1825 double precision zetai, zetaj, zetaij, pi, aij, dij, dum, xj 1826 double precision zero, one 1827 parameter (zero=0.0d+00, one=1.0d+00) 1828c 1829c 1830 double precision derf 1831c 1832c Introduced a trivial replicated data parallelization of this 1833c matrix-vector multiplication 1834c 1835 double precision r, d, factor, factor2 1836 integer i, j 1837 r (i,j)=sqrt((efcc(1,i)-efcc(1,j))**2+ 1838 1 (efcc(2,i)-efcc(2,j))**2+ 1839 2 (efcc(3,i)-efcc(3,j))**2) 1840 d (i )=adiag/sqrt(efcs(i)) 1841c 1842 pi=acos(-1.0d0) 1843 call dfill(nefc,0.0d0,ax,1) 1844c 1845 if (do_cosmo_model.eq.DO_COSMO_KS) then 1846 do i=ga_nodeid()+1,nefc,ga_nnodes() 1847 ax(i) = ax(i) + d(i)*x(i) 1848 do j=i+1,nefc 1849 aij = one/r(i,j) 1850 ax(i) = ax(i) + x(j)*aij 1851 ax(j) = ax(j) + x(i)*aij 1852 enddo 1853 enddo 1854 else if (do_cosmo_model.eq.DO_COSMO_YK) then 1855 factor = zeta*dsqrt(ptspatm/(2d0*pi)) 1856 factor2 = 2d0/sqrt(pi) 1857 do i=ga_nodeid()+1,nefc,ga_nnodes() 1858 zetai = ratm(efciat(i))**2 1859 ax(i) = ax(i) + efcs(i)*x(i) 1860 do j=i+1,nefc 1861 zetaj = ratm(efciat(j))**2 1862 zetaij= factor/sqrt(zetai+zetaj) 1863 dij=r(i,j) 1864 if (dij.lt.1.0d-5) then 1865 aij=factor2*zetaij*(1d0 - (zetaij*dij)**2/3d0) 1866 else 1867 aij=derf(zetaij*dij)/dij 1868 endif 1869 ax(i)=ax(i)+aij*x(j) 1870 ax(j)=ax(j)+aij*x(i) 1871 enddo 1872 enddo 1873 endif 1874c 1875 call ga_dgop(msg_cosaxd,ax,nefc,'+') 1876c 1877 return 1878 end 1879C> 1880C> \brief On the fly vector-matrix multiplication 1881C> 1882C> This routine multiplies the current guess for the COSMO charges 1883C> with the matrix. The routine is replicated data parallel. 1884C> The matrix `A` is symmetric so we simply call the matrix-vector 1885C> product. 1886C> 1887 subroutine hnd_cosxad(nat,x,xa,nefc,efcc,efcs,efciat,ratm) 1888 implicit none 1889 integer nat !< [Input] The number of atoms 1890 integer nefc !< [Input] The number of surface charges 1891 double precision efcc(3,nefc) !< [Input] The surface charge coords 1892 double precision efcs(nefc) !< [Input] The surface areas 1893 integer efciat(nefc) !< [Input] The atom of a surface 1894 !< charge 1895 double precision ratm(nat) !< [Input] The atom radii 1896 double precision x(nefc) !< [Input] Vector `x` 1897 double precision xa(nefc) !< [Output] Vector-matrix product 1898 !< `xA` 1899c 1900 call hnd_cosaxd(nat,x,xa,nefc,efcc,efcs,efciat,ratm) 1901c 1902 return 1903 end 1904C> 1905C> \brief Solve a linear system of equations using an iterative 1906C> procedure 1907C> 1908C> This routine solves a linear system of equations \f$A\cdot x = b\f$ 1909C> using an iterative procedure based on [1] page 70. 1910C> 1911C> ### References ### 1912C> 1913C> [1] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, 1914C> "Numerical Recipes: in Fortran 77", 2nd edition, Volume 1, 1915C> Cambridge University Press, 1992, ISBN: 0-521-43064-X, 1916C> <a href="http://apps.nrbook.com/fortran/index.html">nrbook</a>. 1917C> 1918C> [2] A. Klamt, G. Schüürmann, 1919C> "COSMO: a new approach to dielectric screening in solvents with 1920C> explicit expressions for the screening energy and its gradient", 1921C> <i>J. Chem. Soc., Perkin Trans. 2</i>, 1993, pp 799-805, DOI: 1922C> <a href="https://doi.org/10.1039/P29930000799"> 1923C> 10.1039/P29930000799</a>. 1924C> 1925 subroutine hnd_cosequ(nat,b,x,nefc,g,h,xi,xj,efcc,efcs,efciat, 1926 & ratm) 1927 implicit none 1928c 1929c ----- solve A * x = b , using an iterative procedure ----- 1930c 1931c ----- numerical recipes (p.70), cambridge university press ----- 1932c w.h.press, b.p.flannery, s.a.teukolsky, w.t.vetterling 1933c 1934#include "errquit.fh" 1935#include "stdio.fh" 1936#include "global.fh" 1937c 1938 logical dbug 1939 integer nat !< [Input] The number of atoms 1940 integer nefc !< [Input] The number of surface charges 1941 double precision efcc(3,nefc) !< [Input] The surface charge coords 1942 double precision efcs(nefc) !< [Input] The surface areas 1943 integer efciat(nefc) !< [Input] The atom of a surface 1944 !< charge 1945 double precision ratm(nat) !< [Input] The atom radii 1946 double precision b(nefc) !< [Input] The RHS of Ax=b 1947 double precision x(nefc) !< [Output] The solution of Ax=b 1948c 1949 double precision g(nefc) !< [Scratch] 1950 double precision h(nefc) !< [Scratch] 1951 double precision xi(nefc) !< [Scratch] 1952 double precision xj(nefc) !< [Scratch] 1953c 1954 double precision rp, bsq, anum, aden, rsq, gg, dgg, gam 1955 double precision zero, eps, eps2 1956 data zero /0.0d+00/ 1957 data eps /1.0d-07/ 1958c 1959 integer i, j, irst, iter 1960c 1961 dbug=.false. 1962 if(dbug) then 1963 write(luout,*) 'in -cosequ-' 1964 do i=1,nefc 1965 write(luout,9999) i,b(i) 1966 enddo 1967 endif 1968c 1969 eps2=nefc*eps*eps 1970 irst=0 1971 10 irst=irst+1 1972 call hnd_cosaxd(nat,x,xi,nefc,efcc,efcs,efciat,ratm) 1973 rp=zero 1974 bsq=zero 1975 do j=1,nefc 1976 bsq=bsq+b(j)*b(j) 1977 xi(j)=xi(j)-b(j) 1978 rp=rp+xi(j)*xi(j) 1979 enddo ! j 1980 call hnd_cosxad(nat,xi,g,nefc,efcc,efcs,efciat,ratm) 1981 do j=1,nefc 1982 g(j)=-g(j) 1983 h(j)= g(j) 1984 enddo ! j 1985 do iter=1,10*nefc 1986 call hnd_cosaxd(nat,h,xi,nefc,efcc,efcs,efciat,ratm) 1987 anum=zero 1988 aden=zero 1989 do j=1,nefc 1990 anum=anum+g(j)*h(j) 1991 aden=aden+xi(j)*xi(j) 1992 enddo ! j 1993 if(aden.eq.zero) then 1994 write(luout,*) 'very singular matrix' 1995 call errquit('hnd_cosequ: singular matrix',911,UERR) 1996 endif 1997 anum=anum/aden 1998 do j=1,nefc 1999 xi(j)=x(j) 2000 x(j)=x(j)+anum*h(j) 2001 enddo ! j 2002 call hnd_cosaxd(nat,x,xj,nefc,efcc,efcs,efciat,ratm) 2003 rsq=zero 2004 do j=1,nefc 2005 xj(j)=xj(j)-b(j) 2006 rsq=rsq+xj(j)*xj(j) 2007 enddo ! j 2008 if (iter.gt.10*nefc-min(10,nefc/10)) then 2009 if (ga_nodeid().eq.0) then 2010 write(luout,'(" hnd_cosequ: it,residue,thresh = ", 2011 & i5,2f12.5)')iter,rsq,bsq*eps2 2012 endif 2013 endif 2014 if(rsq.eq.rp.or.rsq.le.bsq*eps2) return 2015 if(rsq.gt.rp) then 2016 do j=1,nefc 2017 x(j)=xi(j) 2018 enddo ! j 2019 if(irst.ge.3) return 2020 go to 10 2021 endif 2022 rp=rsq 2023 call hnd_cosxad(nat,xj,xi,nefc,efcc,efcs,efciat,ratm) 2024 gg=zero 2025 dgg=zero 2026 do j=1,nefc 2027 gg=gg+g(j)*g(j) 2028 dgg=dgg+(xi(j)+g(j))*xi(j) 2029 enddo ! j 2030 if(gg.eq.zero) return 2031 gam=dgg/gg 2032 do j=1,nefc 2033 g(j)=-xi(j) 2034 h(j)=g(j)+gam*h(j) 2035 enddo ! j 2036 enddo ! iter 2037 write(luout,*) 'too many iterations' 2038 call errquit('hnd_cosequ: too many iters',911,UERR) 2039 return 2040 9999 format(i8,f16.8) 2041 end 2042C> 2043C> \brief Direct linear system solver 2044C> 2045C> This is direct linear system solver, solving \f$Ax=b\f$ for \f$x\f$. 2046C> On return the matrix \f$A\f$ has been destroyed, and \f$b\f$ has been 2047C> overwritten with the answer \f$x\f$. 2048C> 2049 subroutine hnd_linequ(a,lda,b,n,ib,t,ierr,nodcmp) 2050#ifdef USE_OPENMP 2051 use omp_lib 2052#endif 2053 implicit none 2054#include "errquit.fh" 2055#include "mafdecls.fh" 2056#include "global.fh" 2057 integer lda !< [Input] The leading dimension of \f$A\f$ 2058 integer n !< [Input] The dimension of matrix \f$A\f$ 2059 double precision a(lda,n) !< [In/Output] On input the matrix 2060 !< \f$A\f$; On output the matrix has been destroyed 2061 double precision b(n,*) !< [In/Output] On input the right-hand-side 2062 !< \f$b\f$; on output the solution \f$x\f$ 2063 double precision t(n) 2064 integer ib(n) 2065 integer ierr !< [Output] Error code 2066 integer nodcmp !< [Input] Flag, if \f$\mathrm{nodcmp}\le 1\f$ 2067 !< then skip the decomposition 2068 integer j !< Counter 2069 integer lwork, iwork, kwork 2070 double precision t0 2071 integer iMaxThreads 2072 2073#ifdef USE_OPENMP 2074 iMaxThreads = omp_get_max_threads() 2075 call util_blas_set_num_threads(iMaxThreads) 2076#endif 2077c 2078c ----- solve a * x = b , with x returned in b ----- 2079c 2080 ierr = 0 2081 2082 ! Save diagonal in case dposv fails 2083 do j=1,n 2084 t(j) = a(j,j) 2085 enddo 2086 2087 call dposv('u', n, 2, a, lda, b, n, ierr) 2088 2089 if (ierr.ne.0) then 2090 2091 lwork = 64*n 2092 if (.not.ma_push_get(mt_dbl,lwork,"dsysv work",iwork,kwork)) 2093 $ call errquit("hnd_linequ: could not allocate work",lwork,MA_ERR) 2094 2095 ! Restore diagonal 2096 do j=1,n 2097 a(j,j) = t(j) 2098 enddo 2099 2100 call dsysv('l',n,2,a,lda,ib,b,n,dbl_mb(kwork),lwork,ierr) 2101 2102 if (.not.ma_chop_stack(iwork)) 2103 $ call errquit("hnd_linequ: could not chop stack",0,MA_ERR) 2104 2105 endif 2106c 2107#ifdef USE_OPENMP 2108 call util_blas_set_num_threads(1) 2109#endif 2110 return 2111 end 2112C> 2113C> \brief Compute the function \f$f(r)\f$ 2114C> 2115C> Computes the function \f$f(r)\f$ as discussed with function 2116C> `hnd_coschg`. This function is Eq.(64) in [1]. 2117C> 2118C> ### References ### 2119C> 2120C> [1] D.M. York, M. Karplus, 2121C> "A smooth solvation potential based on the conductor-like 2122C> screening model", <i>J. Phys. Chem. A</i> (1999) <b>103</b>, 2123C> pp 11060-11079, DOI: 2124C> <a href="https://doi.org/10.1021/jp992097l"> 2125C> 10.1021/jp992097l</a>. 2126C> 2127 double precision function cosff(r) 2128 implicit none 2129#include "errquit.fh" 2130 double precision r !< [Input] penetration fraction 2131c 2132 if (r.lt.0.0d0) then 2133 cosff = 0.0d0 2134 else if (r.gt.1.0d0) then 2135 cosff = 1.0d0 2136 else 2137 cosff = r**3*(10.0d0-15.0d0*r+6.0d0*r**2) 2138 endif 2139c 2140 return 2141 end 2142C> 2143C> \brief Compute the function \f$\frac{\partial f(r)}{\partial r}\f$ 2144C> 2145C> Computes the function \f$\frac{\partial f(r)}{\partial r}\f$ as 2146C> discussed with function `hnd_coschg`. This function is the 2147C> derivative of Eq.(64) in [1]. 2148C> 2149C> ### References ### 2150C> 2151C> [1] D.M. York, M. Karplus, 2152C> "A smooth solvation potential based on the conductor-like 2153C> screening model", <i>J. Phys. Chem. A</i> (1999) <b>103</b>, 2154C> pp 11060-11079, DOI: 2155C> <a href="https://doi.org/10.1021/jp992097l"> 2156C> 10.1021/jp992097l</a>. 2157C> 2158 double precision function cosdff(r) 2159 implicit none 2160#include "errquit.fh" 2161 double precision r !< [Input] penetration fraction 2162c 2163 if (r.lt.0.0d0) then 2164 cosdff = 0.0d0 2165 else if (r.gt.1.0d0) then 2166 cosdff = 0.0d0 2167 else 2168 cosdff = 30.0d0*(r**2)*((1.0d0-r)**2) 2169 endif 2170c 2171 return 2172 end 2173 2174 subroutine hnd_cg(nat,b,x,nefc,y,r,p,ap,efcc,efcs,efciat,ratm) 2175#ifdef USE_OPENMP 2176 use omp_lib 2177#endif 2178 implicit none 2179#include "cosmoP.fh" 2180#include "cosmo_params.fh" 2181#include "stdio.fh" 2182#include "global.fh" 2183 logical dbug 2184 integer nat !< [Input] The number of atoms 2185 integer nefc !< [Input] The number of surface charges 2186 double precision efcc(3,nefc) !< [Input] The surface charge coords 2187 double precision efcs(nefc) !< [Input] The surface areas 2188 integer efciat(nefc) !< [Input] The atom of a surface 2189 !< charge 2190 double precision ratm(nat) !< [Input] The atom radii 2191 double precision b(nefc) !< [Input] The RHS of Ax=b 2192 double precision x(nefc) !< [Output] The solution of Ax=b 2193c 2194 double precision y(nefc) !< [Scratch] 2195 double precision ap(nefc) !< [Scratch] 2196 double precision p(nefc) !< [Scratch] 2197 double precision r(nefc) !< [Scratch] 2198 2199 double precision rtol2, rsqnew, alpha, beta 2200 double precision xnorm, facold, facnew 2201c 2202 integer iter, maxiter 2203c 2204 data rtol2 / 1d-13 / 2205 data maxiter / 100 / 2206 integer iMaxThreads, i 2207 2208#ifdef USE_OPENMP 2209 iMaxThreads = omp_get_max_threads() 2210 call util_blas_set_num_threads(iMaxThreads) 2211#endif 2212 2213 iter = 0 2214 2215 xnorm = dot_product(x,x) 2216 if (xnorm.eq.0d0) then 2217 r(:) = b(:) 2218 else 2219 call hnd_cosaxd(nat,x,r,nefc,efcc,efcs,efciat,ratm) 2220 r(:) = b(:) - r(:) 2221 endif 2222 2223 rsqnew = dot_product(r,r) 2224 if (rsqnew.lt.rtol2) return 2225 2226 2227 if (do_cosmo_model.eq.DO_COSMO_YK) then 2228 p(:) = r(:)/efcs(:) 2229 elseif (do_cosmo_model.eq.DO_COSMO_KS) then 2230 p(:) = dsqrt(efcs(:))*r(:)/adiag 2231 endif 2232 2233 facnew = dot_product(r,p) 2234 2235 do while(iter.le.nefc) 2236 2237 iter = iter + 1 2238 call hnd_cosaxd(nat,p,ap,nefc,efcc,efcs,efciat,ratm) 2239 2240 alpha = facnew/dot_product(p,ap) 2241 call daxpy(nefc,alpha,p,1,x,1) 2242 call daxpy(nefc,-alpha,ap,1,r,1) 2243 2244 rsqnew = dot_product(r,r) 2245 if (rsqnew.lt.rtol2) goto 2000 2246 2247 if (do_cosmo_model.eq.DO_COSMO_YK) then 2248 y(:) = r(:)/efcs(:) 2249 elseif (do_cosmo_model.eq.DO_COSMO_KS) then 2250 y(:) = dsqrt(efcs(:))*r(:)/adiag 2251 endif 2252 2253 facold = facnew 2254 facnew = dot_product(r,y) 2255 beta = facnew/facold 2256 2257 p(:) = y(:) + beta*p(:) 2258 2259 enddo 2260 2261 if (ga_nodeid().eq.0) then 2262 write(luout,1000) dsqrt(rsqnew) 2263 endif 2264 1000 format(2x,"Warning! hnd_cg residual: ",G8.3) 2265 2266 2000 continue 2267 2268#ifdef USE_OPENMP 2269 call util_blas_set_num_threads(1) 2270#endif 2271 return 2272 end 2273 2274C> 2275C> @} 2276c $Id$ 2277