1! 2! PCMSolver, an API for the Polarizable Continuum Model 3! Copyright (C) 2020 Roberto Di Remigio, Luca Frediani and contributors. 4! 5! This file is part of PCMSolver. 6! 7! PCMSolver is free software: you can redistribute it and/or modify 8! it under the terms of the GNU Lesser General Public License as published by 9! the Free Software Foundation, either version 3 of the License, or 10! (at your option) any later version. 11! 12! PCMSolver is distributed in the hope that it will be useful, 13! but WITHOUT ANY WARRANTY; without even the implied warranty of 14! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15! GNU Lesser General Public License for more details. 16! 17! You should have received a copy of the GNU Lesser General Public License 18! along with PCMSolver. If not, see <http://www.gnu.org/licenses/>. 19! 20! For information on the complete list of contributors to the 21! PCMSolver API, see: <http://pcmsolver.readthedocs.io/> 22! 23 24 module pedra_cavity_derivatives 25 26 use pedra_precision 27 28 implicit none 29 30 public cavder 31 32 private 33 34 contains 35 36 subroutine cavder(nsj, nsjr, icoord, intsph, newsph) 37 38#include "pcm_mxcent.inc" 39#include "pcm_nuclei.inc" 40#include "pcm_pcmdef.inc" 41#include "pcm_pcm.inc" 42 43 integer(kind=regint_k) :: nsj, nsjr, icoord 44 integer(kind=regint_k) :: intsph(mxts, 10), newsph(mxsp, 2) 45 46 real(kind=dp), parameter :: d0 = 0.0d0 47 integer(kind=regint_k) :: alge(63), casca(10) 48 real(kind=dp) :: ddr, ddx, ddy, ddz, der, dr1, dx, dy, dz 49 real(kind=dp) :: fact 50 integer(kind=regint_k) :: i, ii, index, k, livel, ll, max, icont, min 51 integer(kind=regint_k) :: ns1, ns2, nsa, nsub, number 52 53 54! Le derivate contengono termini dovuti direttamente allo 55! spostamento del centro della sfera NSJ, e termini "mediati" dagli 56! spostamenti del centro e dal cambiamento del raggio delle sfere 57! "aggiunte" (create da PEDRA, oltre a quelle originarie). 58 59! Memorizza in DERRAD(NS,NSJ,ICOORD) la derivata del raggio di 60! NS e in DERCEN(NS,NSJ,ICOORD,3) le derivate delle 61! coordinate del centro di NS rispetto alla coord. ICOORD della 62! sfera NSJ. 63 64! Se NS e' una sfera originaria queste derivate sono 0, tranne 65! DERCEN(NSJR,NSJ,ICOORD,ICOORD)=1: 66 67 68 69 DERCEN(NSJR,NSJ,ICOORD,ICOORD) = 1.0D+00 70 71! 2) Effetti indiretti. 72! Loop sulle sfere aggiunte 73 74 DO 500 NSA = NESFP+1, NESF 75 DO II = 1, 63 76 ALGE(II) = 0 77 end do 78 79 ! Costruiamo l'"albero genealogico" della sfera NSA 80 81 ALGE(1) = NSA 82 ALGE(2) = ABS(NEWSPH(NSA,1)) 83 ALGE(3) = ABS(NEWSPH(NSA,2)) 84 LIVEL = 3 85 NUMBER = 2 86 510 NSUB = 1 87 DO II = LIVEL-NUMBER+1, LIVEL 88 IF(ALGE(II) > NESFP) THEN 89 ALGE(LIVEL+NSUB) = ABS(NEWSPH(ALGE(II),1)) 90 ALGE(LIVEL+NSUB+1) = ABS(NEWSPH(ALGE(II),2)) 91 end if 92 NSUB = NSUB + 2 93 end do 94 NUMBER = NUMBER * 2 95 LIVEL = LIVEL + NUMBER 96 IF(NUMBER < 32) go to 510 97 98 ! Quando un elemento di ALGE e' = NSJR, costruisce la corrispondente 99 ! "cascata" di sfere aggiunte che collega NSJR a NSA 100 101 DO 600 LIVEL = 2, 6 102 MIN = 2**(LIVEL-1) 103 MAX = (2**LIVEL) - 1 104 DO 700 II = MIN, MAX 105 IF(ALGE(II) /= NSJR) go to 700 106 DO K = 1, 10 107 CASCA(K) = 0 108 end do 109 CASCA(1) = NSJR 110 INDEX = II 111 K = 2 112 DO LL = LIVEL, 2, -1 113 FACT = (INDEX - 2**(LL-1)) / 2.0D+00 114 INDEX = INT(2**(LL-2) + FACT) 115 CASCA(K) = ALGE(INDEX) 116 K = K + 1 117 end do 118 ! Contiamo gli elementi diversi da 0 in CASCA 119 ICONT = 0 120 DO K = 1, 10 121 IF(CASCA(K) /= 0) ICONT = ICONT + 1 122 end do 123 124 ! Costruiamo le derivate composte del raggio e delle coordinate di 125 ! NSA (ultimo elemento di CASCA) 126 ! rispetto alla coordinata ICOORD di NSJ (primo elemento di CASCA) 127 128 NS1 = CASCA(1) 129 NS2 = CASCA(2) 130 CALL DRRDCN(NS2,ICOORD,NS1,DR1,NEWSPH) 131 CALL DRCNCN(1,NS2,ICOORD,NS1,DX,NEWSPH) 132 CALL DRCNCN(2,NS2,ICOORD,NS1,DY,NEWSPH) 133 CALL DRCNCN(3,NS2,ICOORD,NS1,DZ,NEWSPH) 134 DO 800 I = 3, ICONT 135 DDR = D0 136 DDX = D0 137 DDY = D0 138 DDZ = D0 139 NS1 = CASCA(I-1) 140 NS2 = CASCA(I) 141 142 ! Derivata del raggio dell'elemento I di CASCA rispetto 143 ! alla coord. ICOORD dell'elemento 1 di CASCA 144 145 CALL DRRDRD(NS2,NS1,DER,NEWSPH) 146 DDR = DER * DR1 147 CALL DRRDCN(NS2,1,NS1,DER,NEWSPH) 148 DDR = DDR + DER * DX 149 CALL DRRDCN(NS2,2,NS1,DER,NEWSPH) 150 DDR = DDR + DER * DY 151 CALL DRRDCN(NS2,3,NS1,DER,NEWSPH) 152 DDR = DDR + DER * DZ 153 154 ! Derivata della coord. X dell'elemento I di CASCA rispetto 155 ! alla coord. ICOORD dell'elemento 1 di CASCA 156 157 CALL DRCNRD(1,NS2,NS1,DER,NEWSPH) 158 DDX = DER * DR1 159 CALL DRCNCN(1,NS2,1,NS1,DER,NEWSPH) 160 DDX = DDX + DER * DX 161 CALL DRCNCN(1,NS2,2,NS1,DER,NEWSPH) 162 DDX = DDX + DER * DY 163 CALL DRCNCN(1,NS2,3,NS1,DER,NEWSPH) 164 DDX = DDX + DER * DZ 165 166 ! Derivata della coord. Y dell'elemento I di CASCA rispetto 167 ! alla coord. ICOORD dell'elemento 1 di CASCA 168 169 CALL DRCNRD(2,NS2,NS1,DER,NEWSPH) 170 DDY = DER * DR1 171 CALL DRCNCN(2,NS2,1,NS1,DER,NEWSPH) 172 DDY = DDY + DER * DX 173 CALL DRCNCN(2,NS2,2,NS1,DER,NEWSPH) 174 DDY = DDY + DER * DY 175 CALL DRCNCN(2,NS2,3,NS1,DER,NEWSPH) 176 DDY = DDY + DER * DZ 177 178 ! Derivata della coord. Z dell'elemento I di CASCA rispetto 179 ! alla coord. ICOORD dell'elemento 1 di CASCA 180 181 CALL DRCNRD(3,NS2,NS1,DER,NEWSPH) 182 DDZ = DER * DR1 183 CALL DRCNCN(3,NS2,1,NS1,DER,NEWSPH) 184 DDZ = DDZ + DER * DX 185 CALL DRCNCN(3,NS2,2,NS1,DER,NEWSPH) 186 DDZ = DDZ + DER * DY 187 CALL DRCNCN(3,NS2,3,NS1,DER,NEWSPH) 188 DDZ = DDZ + DER * DZ 189 190 DR1 = DDR 191 DX = DDX 192 DY = DDY 193 DZ = DDZ 194 800 end do 195 196 ! Se NS e' una sfera aggiunta, memorizza le derivate del raggio 197 ! e delle coordinate del centro: 198 199 DERRAD(NSA,NSJ,ICOORD) = DR1 200 DERCEN(NSA,NSJ,ICOORD,1) = DX 201 DERCEN(NSA,NSJ,ICOORD,2) = DY 202 DERCEN(NSA,NSJ,ICOORD,3) = DZ 203 700 end do 204 600 end do 205 500 end do 206 207 end subroutine cavder 208 209 subroutine drcnrd(jj,nsi,nsj,dc,newsph) 210 211#include "pcm_mxcent.inc" 212#include "pcm_nuclei.inc" 213#include "pcm_pcmdef.inc" 214#include "pcm_pcm.inc" 215 216 integer(kind=regint_k) :: jj, nsi, nsj 217 real(kind=dp) :: coordj(3), coordk(3) 218 integer(kind=regint_k) :: intsph(mxts, 10), newsph(mxsp, 2) 219 220 real(kind=dp), parameter :: d0 = 0.0d0 221 real(kind=dp) :: dc, d, d2 222 integer(kind=regint_k) :: nsk 223 224! Trova la derivata della coordinata JJ del centro della sfera 225! NSI rispetto al raggio dellla sfera NSJ. 226 227! La sfera NSI (che appartiene alle sfere "aggiunte" da PEDRA_) 228! dipende dalle due sfere "precedenti" NSJ e NSK 229 230! Se NSJ o NSK sono negativi, la sfera aggiunta e' di tipo C 231! e la generatrice "principale" corrisponde al label negativo 232! (cfr. JCC 11, 1047 (1990)) 233 234 IF(NEWSPH(NSI,1) < 0 .OR. NEWSPH(NSI,2) < 0) go to 100 235 NSK = NEWSPH(NSI,1) 236 IF(NSK == NSJ) NSK = NEWSPH(NSI,2) 237 COORDJ(1) = XE(NSJ) 238 COORDJ(2) = YE(NSJ) 239 COORDJ(3) = ZE(NSJ) 240 COORDK(1) = XE(NSK) 241 COORDK(2) = YE(NSK) 242 COORDK(3) = ZE(NSK) 243 D2 = (XE(NSJ)-XE(NSK))**2 + (YE(NSJ)-YE(NSK))**2 + & 244 (ZE(NSJ)-ZE(NSK))**2 245 D = SQRT(D2) 246 DC = - (COORDJ(JJ) - COORDK(JJ)) / (2.0D+00 * D) 247 go to 200 248 249 100 CONTINUE 250 NSK = NEWSPH(NSI,1) 251 IF(ABS(NSK) == NSJ) NSK = NEWSPH(NSI,2) 252 DC = D0 253 IF(NSK < D0) go to 200 254 COORDJ(1) = XE(NSJ) 255 COORDJ(2) = YE(NSJ) 256 COORDJ(3) = ZE(NSJ) 257 COORDK(1) = XE(NSK) 258 COORDK(2) = YE(NSK) 259 COORDK(3) = ZE(NSK) 260 D2 = (XE(NSJ)-XE(NSK))**2 + (YE(NSJ)-YE(NSK))**2 + & 261 (ZE(NSJ)-ZE(NSK))**2 262 D = SQRT(D2) 263 DC = - ( COORDJ(JJ) - COORDK(JJ) ) / D 264 265 200 CONTINUE 266 end subroutine drcnrd 267 268 subroutine drcncn(jj,nsi,icoord,nsj,dc,newsph) 269 270#include "pcm_mxcent.inc" 271#include "pcm_nuclei.inc" 272#include "pcm_pcmdef.inc" 273#include "pcm_pcm.inc" 274 275 integer(kind=regint_k) :: jj, nsi, icoord, nsj 276 real(kind=dp) :: dc 277 integer(kind=regint_k) :: newsph(mxsp,2) 278 279 real(kind=dp) :: coordj(3), coordk(3) 280 real(kind=dp), parameter :: d0 = 0.0d0 281 real(kind=dp) :: d, d2 282 integer(kind=regint_k) :: k, nsk 283 284! Trova la derivata della coordinata JJ del centro della sfera 285! NSI rispetto alla coordinata ICOORD di NSJ, che interseca NSI. 286 287! La sfera NSI (che appartiene alle sfere "aggiunte" da PEDRA_) 288! dipende dalle due sfere "precedenti" NSJ e NSK 289 290! Se NSJ o NSK sono negativi, la sfera aggiunta e' di tipo C 291! e la generatrice "principale" corrisponde al label negativo 292! (cfr. JCC 11, 1047 (1990)) 293 294 IF(NEWSPH(NSI,1) < 0 .OR. NEWSPH(NSI,2) < 0) go to 100 295 K = NEWSPH(NSI,1) 296 IF(K == NSJ) K = NEWSPH(NSI,2) 297 COORDJ(1) = XE(NSJ) 298 COORDJ(2) = YE(NSJ) 299 COORDJ(3) = ZE(NSJ) 300 COORDK(1) = XE(K) 301 COORDK(2) = YE(K) 302 COORDK(3) = ZE(K) 303 D2 = (XE(NSJ)-XE(K))**2 + (YE(NSJ)-YE(K))**2 + (ZE(NSJ)-ZE(K))**2 304 D = SQRT(D2) 305 DC = (RE(NSJ)-RE(K)) * (COORDJ(ICOORD)-COORDK(ICOORD)) * & 306 (COORDJ(JJ) - COORDK(JJ)) / (2.0D+00 * D**3) 307 IF(JJ == ICOORD)DC = DC + 0.5D+00 - (RE(NSJ)-RE(K)) / (2.0D+00*D) 308 go to 200 309 310 100 CONTINUE 311 NSK = NEWSPH(NSI,1) 312 IF(ABS(NSK) == NSJ) NSK = NEWSPH(NSI,2) 313 COORDJ(1) = XE(NSJ) 314 COORDJ(2) = YE(NSJ) 315 COORDJ(3) = ZE(NSJ) 316 COORDK(1) = XE(ABS(NSK)) 317 COORDK(2) = YE(ABS(NSK)) 318 COORDK(3) = ZE(ABS(NSK)) 319 D2 = (COORDJ(1)-COORDK(1))**2 + (COORDJ(2)-COORDK(2))**2 + & 320 (COORDJ(3)-COORDK(3))**2 321 D = SQRT(D2) 322 IF(NSK > 0) THEN 323 DC = RE(NSJ) * (COORDJ(JJ)-COORDK(JJ)) * (COORDJ(ICOORD)- & 324 COORDK(ICOORD)) / D**3 325 IF(ICOORD == JJ) DC = DC + 1.0D+00 - RE(NSJ) / D 326 ELSE 327 DC = - RE(ABS(NSK)) * (COORDK(JJ)-COORDJ(JJ)) * (COORDK(ICOORD)- & 328 COORDJ(ICOORD)) / D**3 329 IF(ICOORD == JJ) DC = DC + RE(ABS(NSK)) / D 330 end if 331 332 200 CONTINUE 333 end subroutine drcncn 334 335 subroutine drrdrd(nsi,nsj,dr1,newsph) 336 337#include "pcm_mxcent.inc" 338#include "pcm_nuclei.inc" 339#include "pcm_pcmdef.inc" 340#include "pcm_pcm.inc" 341 342 integer(kind=regint_k) :: nsi, nsj, newsph(mxsp, 2) 343 real(kind=dp) :: dr1 344 345 real(kind=dp), parameter :: d0 = 0.0d0 346 real(kind=dp) :: d, d2, ri, rj, rk, rs 347 integer(kind=regint_k) :: nsk 348 349! Trova la derivata del raggio della sfera NSI rispetto al raggio 350! della sfera NSJ. 351 352! La sfera NSI (che appartiene alle sfere "aggiunte" da PEDRA_) 353! dipende dalle due sfere "precedenti" NSJ e NSK 354! Se NSJ o NSK sono negativi, la sfera aggiunta e' di tipo C 355! e la generatrice "principale" corrisponde al label negativo 356! (cfr. JCC 11, 1047 (1990)) 357 358 IF(NEWSPH(NSI,1) < 0 .OR. NEWSPH(NSI,2) < 0) go to 100 359 NSK = NEWSPH(NSI,1) 360 IF(NSK == NSJ) NSK = NEWSPH(NSI,2) 361 RS = RSOLV 362 RJ = RE(NSJ) + RS 363 RK = RE(NSK) + RS 364 RI = RE(NSI) + RS 365 D2 = (XE(NSJ)-XE(NSK))**2 + (YE(NSJ)-YE(NSK))**2 + & 366 (ZE(NSJ)-ZE(NSK))**2 367 D = SQRT(D2) 368 DR1 = (-3.0D+00*RJ*RJ + RK*RK + 2.0D+00*RJ*RK & 369 + 3.0D+00*D*RJ - D*RK) / (4.0D+00*D*RI) 370 go to 200 371 372 100 CONTINUE 373 NSK = NEWSPH(NSI,1) 374 IF(ABS(NSK) == NSJ) NSK = NEWSPH(NSI,2) 375 376 IF(NSK > 0) THEN 377 RS = RSOLV 378 RJ = RE(NSJ) + RS 379 RK = RE(NSK) + RS 380 RI = RE(NSI) + RS 381 D2 = (XE(NSJ)-XE(NSK))**2 + (YE(NSJ)-YE(NSK))**2 + & 382 (ZE(NSJ)-ZE(NSK))**2 383 D = SQRT(D2) 384 DR1 = ( 2.0D+00*D*RJ + 2.0D+00*D*RE(NSJ) - 2.0D+00*RJ*RE(NSJ) + & 385 D*D - RJ*RJ - RK*RK ) / (2.0D+00*D*RI) 386 ELSE 387 RS = RSOLV 388 RJ = RE(NSJ) + RS 389 RI = RE(NSI) + RS 390 D2 = (XE(NSJ)-XE(ABS(NSK)))**2 + (YE(NSJ)-YE(ABS(NSK)))**2 + & 391 (ZE(NSJ)-ZE(ABS(NSK)))**2 392 D = SQRT(D2) 393 DR1 = ( RE(ABS(NSK)) * RJ ) / ( D*RI) 394 end if 395 200 CONTINUE 396 end subroutine drrdrd 397 398 subroutine drrdcn(nsi, icoord, nsj, dr1, newsph) 399 400#include "pcm_mxcent.inc" 401#include "pcm_nuclei.inc" 402#include "pcm_pcmdef.inc" 403#include "pcm_pcm.inc" 404 405 integer(kind=regint_k) :: nsi, icoord, nsj, newsph(mxsp, 2) 406 real(kind=dp) :: dr1 407 408 real(kind=dp) :: coordj(3), coordk(3) 409 real(kind=dp), parameter :: d0 = 0.0d0 410 real(kind=dp) :: a, b, d, d2, diff, fac, ri, rj, rk, rs 411 integer(kind=regint_k) :: k, nsk 412 413! Trova la derivata del raggio della sfera NSI rispetto alla 414! coordinata ICOORD (1=X, 2=Y, 3=Z) della sfera NSJ, che interseca 415! NSI. 416 417! La sfera NSI (che appartiene alle sfere "aggiunte" da PEDRA_) 418! dipende dalle due sfere "precedenti" NSJ e K 419 420! Se NSJ o NSK sono negativi, la sfera aggiunta e' di tipo C 421! e la generatrice "principale" corrisponde al label negativo 422! (cfr. JCC 11, 1047 (1990)) 423 424 IF(NEWSPH(NSI,1) < 0 .OR. NEWSPH(NSI,2) < 0) go to 100 425 K = NEWSPH(NSI,1) 426 IF(K == NSJ) K = NEWSPH(NSI,2) 427 COORDJ(1) = XE(NSJ) 428 COORDJ(2) = YE(NSJ) 429 COORDJ(3) = ZE(NSJ) 430 COORDK(1) = XE(K) 431 COORDK(2) = YE(K) 432 COORDK(3) = ZE(K) 433 D2 = (XE(NSJ)-XE(K))**2 + (YE(NSJ)-YE(K))**2 + (ZE(NSJ)-ZE(K))**2 434 D = SQRT(D2) 435 B = 0.5D+00 * (D + RE(NSJ) - RE(K)) 436 RS = RSOLV 437 A = ((RE(NSJ)+RS)**2 + D2 - (RE(K)+RS)**2) / D 438 DR1 = (2.0D+00*A*B - 2.0D+00*B*D - A*D) * & 439 (COORDJ(ICOORD)-COORDK(ICOORD)) / (4.0D+00*D2*(RE(NSI)+RS)) 440 go to 200 441 442 100 CONTINUE 443 NSK = NEWSPH(NSI,1) 444 IF(ABS(NSK) == NSJ) NSK = NEWSPH(NSI,2) 445 COORDJ(1) = XE(NSJ) 446 COORDJ(2) = YE(NSJ) 447 COORDJ(3) = ZE(NSJ) 448 COORDK(1) = XE(ABS(NSK)) 449 COORDK(2) = YE(ABS(NSK)) 450 COORDK(3) = ZE(ABS(NSK)) 451 RI = RE(NSI) + RSOLV 452 RJ = RE(NSJ) + RSOLV 453 RK = RE(ABS(NSK)) + RSOLV 454 DIFF = COORDJ(ICOORD) - COORDK(ICOORD) 455 D2 = (COORDJ(1)-COORDK(1))**2 + (COORDJ(2)-COORDK(2))**2 + & 456 (COORDJ(3)-COORDK(3))**2 457 D = SQRT(D2) 458 FAC = RE(NSJ) * ( RJ*RJ - D*D - RK*RK ) 459 IF(NSK < 0) FAC = RE(ABS(NSK)) * (RK*RK - D*D - RJ*RJ ) 460 DR1 = DIFF * FAC / ( 2.0D+00 * D**3 * RI ) 461 462 200 CONTINUE 463 end subroutine drrdcn 464 465 end module pedra_cavity_derivatives 466