1c 2c 3c ############################################################ 4c ## COPYRIGHT (C) 2018 by Joshua Rackers & Jay W. Ponder ## 5c ## All Rights Reserved ## 6c ############################################################ 7c 8c ################################################################ 9c ## ## 10c ## subroutine echgtrn1 -- charge transfer energy & derivs ## 11c ## ## 12c ################################################################ 13c 14c 15c "echgtrn1" calculates the charge transfer energy and first 16c derivatives with respect to Cartesian coordinates 17c 18c 19 subroutine echgtrn1 20 use limits 21 implicit none 22c 23c 24c choose method for summing over charge transfer interactions 25c 26 if (use_mlist) then 27 call echgtrn1b 28 else 29 call echgtrn1a 30 end if 31 return 32 end 33c 34c 35c ################################################################# 36c ## ## 37c ## subroutine echgtrn1a -- charge transfer derivs via loop ## 38c ## ## 39c ################################################################# 40c 41c 42c "echgtrn1a" calculates the charge transfer interaction energy 43c and first derivatives using a double loop 44c 45c 46 subroutine echgtrn1a 47 use atoms 48 use bound 49 use chgpot 50 use chgtrn 51 use cell 52 use couple 53 use ctrpot 54 use deriv 55 use energi 56 use group 57 use mplpot 58 use mpole 59 use shunt 60 use usage 61 use virial 62 implicit none 63 integer i,j,k 64 integer ii,kk 65 integer jcell 66 real*8 e,de,f,fgrp 67 real*8 rr1,r,r2 68 real*8 r3,r4,r5 69 real*8 xi,yi,zi 70 real*8 xr,yr,zr 71 real*8 chgi,chgk 72 real*8 chgik 73 real*8 alphai,alphak 74 real*8 alphaik 75 real*8 expi,expk 76 real*8 expik 77 real*8 frcx,frcy,frcz 78 real*8 vxx,vyy,vzz 79 real*8 vxy,vxz,vyz 80 real*8 taper,dtaper 81 real*8, allocatable :: mscale(:) 82 logical proceed,usei 83 character*6 mode 84c 85c 86c zero out the charge transfer energy and first derivatives 87c 88 ect = 0.0d0 89 do i = 1, n 90 dect(1,i) = 0.0d0 91 dect(2,i) = 0.0d0 92 dect(3,i) = 0.0d0 93 end do 94 if (nct .eq. 0) return 95c 96c perform dynamic allocation of some local arrays 97c 98 allocate (mscale(n)) 99c 100c initialize connected atom exclusion coefficients 101c 102 do i = 1, n 103 mscale(i) = 1.0d0 104 end do 105c 106c set conversion factor, cutoff and switching coefficients 107c 108 f = electric / dielec 109 mode = 'CHGTRN' 110 call switch (mode) 111c 112c calculate the charge transfer energy and derivatives 113c 114 do ii = 1, npole-1 115 i = ipole(ii) 116 xi = x(i) 117 yi = y(i) 118 zi = z(i) 119 chgi = chgct(ii) 120 alphai = dmpct(ii) 121 if (alphai .eq. 0.0d0) alphai = 1000.0d0 122 usei = use(i) 123c 124c set exclusion coefficients for connected atoms 125c 126 do j = 1, n12(i) 127 mscale(i12(j,i)) = m2scale 128 end do 129 do j = 1, n13(i) 130 mscale(i13(j,i)) = m3scale 131 end do 132 do j = 1, n14(i) 133 mscale(i14(j,i)) = m4scale 134 end do 135 do j = 1, n15(i) 136 mscale(i15(j,i)) = m5scale 137 end do 138c 139c evaluate all sites within the cutoff distance 140c 141 do kk = ii+1, npole 142 k = ipole(kk) 143 proceed = .true. 144 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 145 if (.not. use_intra) proceed = .true. 146 if (proceed) proceed = (usei .or. use(k)) 147 if (proceed) then 148 xr = x(k) - xi 149 yr = y(k) - yi 150 zr = z(k) - zi 151 if (use_bounds) call image (xr,yr,zr) 152 r2 = xr*xr + yr* yr + zr*zr 153 if (r2 .le. off2) then 154 r = sqrt(r2) 155 rr1 = 1.0d0 / r 156 chgk = chgct(kk) 157 alphak = dmpct(kk) 158 if (alphak .eq. 0.0d0) alphak = 1000.0d0 159 if (ctrntyp .eq. 'SEPARATE') then 160 expi = exp(-alphai*r) 161 expk = exp(-alphak*r) 162 e = -chgi*expk - chgk*expi 163 de = chgi*expk*alphak + chgk*expi*alphai 164 else 165 chgik = sqrt(abs(chgi*chgk)) 166 alphaik = 0.5d0 * (alphai+alphak) 167 expik = exp(-alphaik*r) 168 e = -chgik * expik 169 de = -e * alphaik 170 end if 171 e = f * e * mscale(k) 172 de = f * de * mscale(k) 173c 174c use energy switching if near the cutoff distance 175c 176 if (r2 .gt. cut2) then 177 r3 = r2 * r 178 r4 = r2 * r2 179 r5 = r2 * r3 180 taper = c5*r5 + c4*r4 + c3*r3 181 & + c2*r2 + c1*r + c0 182 dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3 183 & + 3.0d0*c3*r2 + 2.0d0*c2*r + c1 184 de = e*dtaper + de*taper 185 e = e * taper 186 end if 187c 188c scale the interaction based on its group membership 189c 190 if (use_group) then 191 e = e * fgrp 192 de = de * fgrp 193 end if 194c 195c compute the force components for this interaction 196c 197 frcx = de * xr * rr1 198 frcy = de * yr * rr1 199 frcz = de * zr * rr1 200c 201c increment the total charge transfer energy and derivatives 202c 203 ect = ect + e 204 dect(1,i) = dect(1,i) - frcx 205 dect(2,i) = dect(2,i) - frcy 206 dect(3,i) = dect(3,i) - frcz 207 dect(1,k) = dect(1,k) + frcx 208 dect(2,k) = dect(2,k) + frcy 209 dect(3,k) = dect(3,k) + frcz 210c 211c increment the internal virial tensor components 212c 213 vxx = xr * frcx 214 vxy = yr * frcx 215 vxz = zr * frcx 216 vyy = yr * frcy 217 vyz = zr * frcy 218 vzz = zr * frcz 219 vir(1,1) = vir(1,1) + vxx 220 vir(2,1) = vir(2,1) + vxy 221 vir(3,1) = vir(3,1) + vxz 222 vir(1,2) = vir(1,2) + vxy 223 vir(2,2) = vir(2,2) + vyy 224 vir(3,2) = vir(3,2) + vyz 225 vir(1,3) = vir(1,3) + vxz 226 vir(2,3) = vir(2,3) + vyz 227 vir(3,3) = vir(3,3) + vzz 228 end if 229 end if 230 end do 231c 232c reset exclusion coefficients for connected atoms 233c 234 do j = 1, n12(i) 235 mscale(i12(j,i)) = 1.0d0 236 end do 237 do j = 1, n13(i) 238 mscale(i13(j,i)) = 1.0d0 239 end do 240 do j = 1, n14(i) 241 mscale(i14(j,i)) = 1.0d0 242 end do 243 do j = 1, n15(i) 244 mscale(i15(j,i)) = 1.0d0 245 end do 246 end do 247c 248c for periodic boundary conditions with large cutoffs 249c neighbors must be found by the replicates method 250c 251 if (use_replica) then 252c 253c calculate interaction components with other unit cells 254c 255 do ii = 1, npole 256 i = ipole(ii) 257 xi = x(i) 258 yi = y(i) 259 zi = z(i) 260 chgi = chgct(ii) 261 alphai = dmpct(ii) 262 if (alphai .eq. 0.0d0) alphai = 1000.0d0 263 usei = use(i) 264c 265c set exclusion coefficients for connected atoms 266c 267 do j = 1, n12(i) 268 mscale(i12(j,i)) = m2scale 269 end do 270 do j = 1, n13(i) 271 mscale(i13(j,i)) = m3scale 272 end do 273 do j = 1, n14(i) 274 mscale(i14(j,i)) = m4scale 275 end do 276 do j = 1, n15(i) 277 mscale(i15(j,i)) = m5scale 278 end do 279c 280c evaluate all sites within the cutoff distance 281c 282 do kk = ii, npole 283 k = ipole(kk) 284 proceed = .true. 285 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 286 if (.not. use_intra) proceed = .true. 287 if (proceed) proceed = (usei .or. use(k)) 288 if (proceed) then 289 do jcell = 2, ncell 290 xr = x(k) - xi 291 yr = y(k) - yi 292 zr = z(k) - zi 293 call imager (xr,yr,zr,jcell) 294 r2 = xr*xr + yr* yr + zr*zr 295 if (r2 .le. off2) then 296 r = sqrt(r2) 297 rr1 = 1.0d0 / r 298 chgk = chgct(kk) 299 alphak = dmpct(kk) 300 if (alphak .eq. 0.0d0) alphak = 1000.0d0 301 if (ctrntyp .eq. 'SEPARATE') then 302 expi = exp(-alphai*r) 303 expk = exp(-alphak*r) 304 e = -chgi*expk - chgk*expi 305 de = chgi*expk*alphak + chgk*expi*alphai 306 else 307 chgik = sqrt(abs(chgi*chgk)) 308 alphaik = 0.5d0 * (alphai+alphak) 309 expik = exp(-alphaik*r) 310 e = -chgik * expik 311 de = -e * alphaik 312 end if 313 if (use_group) then 314 e = e * fgrp 315 de = de * fgrp 316 end if 317 e = f * e * mscale(k) 318 de = f * de * mscale(k) 319 if (i .eq. k) then 320 e = 0.5d0 * e 321 de = 0.5d0 * de 322 end if 323c 324c use energy switching if near the cutoff distance 325c 326 if (r2 .gt. cut2) then 327 r3 = r2 * r 328 r4 = r2 * r2 329 r5 = r2 * r3 330 taper = c5*r5 + c4*r4 + c3*r3 331 & + c2*r2 + c1*r + c0 332 dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3 333 & + 3.0d0*c3*r2 + 2.0d0*c2*r + c1 334 de = e*dtaper + de*taper 335 e = e * taper 336 end if 337c 338c scale the interaction based on its group membership 339c 340 if (use_group) then 341 e = e * fgrp 342 de = de * fgrp 343 end if 344c 345c compute the force components for this interaction 346c 347 frcx = de * xr * rr1 348 frcy = de * yr * rr1 349 frcz = de * zr * rr1 350c 351c increment the total charge transfer energy and derivatives 352c 353 ect = ect + e 354 dect(1,i) = dect(1,i) - frcx 355 dect(2,i) = dect(2,i) - frcy 356 dect(3,i) = dect(3,i) - frcz 357 dect(1,k) = dect(1,k) + frcx 358 dect(2,k) = dect(2,k) + frcy 359 dect(3,k) = dect(3,k) + frcz 360c 361c increment the internal virial tensor components 362c 363 vxx = xr * frcx 364 vxy = yr * frcx 365 vxz = zr * frcx 366 vyy = yr * frcy 367 vyz = zr * frcy 368 vzz = zr * frcz 369 vir(1,1) = vir(1,1) + vxx 370 vir(2,1) = vir(2,1) + vxy 371 vir(3,1) = vir(3,1) + vxz 372 vir(1,2) = vir(1,2) + vxy 373 vir(2,2) = vir(2,2) + vyy 374 vir(3,2) = vir(3,2) + vyz 375 vir(1,3) = vir(1,3) + vxz 376 vir(2,3) = vir(2,3) + vyz 377 vir(3,3) = vir(3,3) + vzz 378 end if 379 end do 380 end if 381 end do 382c 383c reset exclusion coefficients for connected atoms 384c 385 do j = 1, n12(i) 386 mscale(i12(j,i)) = 1.0d0 387 end do 388 do j = 1, n13(i) 389 mscale(i13(j,i)) = 1.0d0 390 end do 391 do j = 1, n14(i) 392 mscale(i14(j,i)) = 1.0d0 393 end do 394 do j = 1, n15(i) 395 mscale(i15(j,i)) = 1.0d0 396 end do 397 end do 398 end if 399c 400c perform deallocation of some local arrays 401c 402 deallocate (mscale) 403 return 404 end 405c 406c 407c ################################################################# 408c ## ## 409c ## subroutine echgtrn1b -- charge transfer derivs via list ## 410c ## ## 411c ################################################################# 412c 413c 414c "echgtrn1b" calculates the charge transfer energy and first 415c derivatives using a pairwise neighbor list 416c 417c 418 subroutine echgtrn1b 419 use atoms 420 use bound 421 use chgpot 422 use chgtrn 423 use cell 424 use couple 425 use ctrpot 426 use deriv 427 use energi 428 use group 429 use mplpot 430 use mpole 431 use neigh 432 use shunt 433 use usage 434 use virial 435 implicit none 436 integer i,j,k 437 integer ii,kk,kkk 438 real*8 e,de,f,fgrp 439 real*8 rr1,r,r2 440 real*8 r3,r4,r5 441 real*8 xi,yi,zi 442 real*8 xr,yr,zr 443 real*8 chgi,chgk 444 real*8 chgik 445 real*8 alphai,alphak 446 real*8 alphaik 447 real*8 expi,expk 448 real*8 expik 449 real*8 frcx,frcy,frcz 450 real*8 vxx,vyy,vzz 451 real*8 vxy,vxz,vyz 452 real*8 taper,dtaper 453 real*8, allocatable :: mscale(:) 454 logical proceed,usei 455 character*6 mode 456c 457c 458c zero out the charge transfer energy and first derivatives 459c 460 ect = 0.0d0 461 do i = 1, n 462 dect(1,i) = 0.0d0 463 dect(2,i) = 0.0d0 464 dect(3,i) = 0.0d0 465 end do 466 if (nct .eq. 0) return 467c 468c perform dynamic allocation of some local arrays 469c 470 allocate (mscale(n)) 471c 472c initialize connected atom exclusion coefficients 473c 474 do i = 1, n 475 mscale(i) = 1.0d0 476 end do 477c 478c set conversion factor, cutoff and switching coefficients 479c 480 f = electric / dielec 481 mode = 'CHGTRN' 482 call switch (mode) 483c 484c OpenMP directives for the major loop structure 485c 486!$OMP PARALLEL default(private) 487!$OMP& shared(npole,ipole,x,y,z,chgct,dmpct,n12,i12,n13,i13, 488!$OMP& n14,i14,n15,i15,m2scale,m3scale,m4scale,m5scale,nelst, 489!$OMP& elst,use,use_group,use_intra,use_bounds,ctrntyp,f,off2, 490!$OMP& cut2,c0,c1,c2,c3,c4,c5) 491!$OMP& firstprivate(mscale) shared(ect,dect,vir) 492!$OMP DO reduction(+:ect,dect,vir) schedule(guided) 493c 494c compute the charge transfer energy and derivatives 495c 496 do ii = 1, npole 497 i = ipole(ii) 498 xi = x(i) 499 yi = y(i) 500 zi = z(i) 501 chgi = chgct(ii) 502 alphai = dmpct(ii) 503 if (alphai .eq. 0.0d0) alphai = 1000.0d0 504 usei = use(i) 505c 506c set exclusion coefficients for connected atoms 507c 508 do j = 1, n12(i) 509 mscale(i12(j,i)) = m2scale 510 end do 511 do j = 1, n13(i) 512 mscale(i13(j,i)) = m3scale 513 end do 514 do j = 1, n14(i) 515 mscale(i14(j,i)) = m4scale 516 end do 517 do j = 1, n15(i) 518 mscale(i15(j,i)) = m5scale 519 end do 520c 521c evaluate all sites within the cutoff distance 522c 523 do kkk = 1, nelst(ii) 524 kk = elst(kkk,ii) 525 k = ipole(kk) 526 proceed = .true. 527 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 528 if (.not. use_intra) proceed = .true. 529 if (proceed) proceed = (usei .or. use(k)) 530 if (proceed) then 531 xr = x(k) - xi 532 yr = y(k) - yi 533 zr = z(k) - zi 534 if (use_bounds) call image (xr,yr,zr) 535 r2 = xr*xr + yr* yr + zr*zr 536 if (r2 .le. off2) then 537 r = sqrt(r2) 538 rr1 = 1.0d0 / r 539 chgk = chgct(kk) 540 alphak = dmpct(kk) 541 if (alphak .eq. 0.0d0) alphak = 1000.0d0 542 if (ctrntyp .eq. 'SEPARATE') then 543 expi = exp(-alphai*r) 544 expk = exp(-alphak*r) 545 e = -chgi*expk - chgk*expi 546 de = chgi*expk*alphak + chgk*expi*alphai 547 else 548 chgik = sqrt(abs(chgi*chgk)) 549 alphaik = 0.5d0 * (alphai+alphak) 550 expik = exp(-alphaik*r) 551 e = -chgik * expik 552 de = -e * alphaik 553 end if 554 e = f * e * mscale(k) 555 de = f * de * mscale(k) 556c 557c use energy switching if near the cutoff distance 558c 559 if (r2 .gt. cut2) then 560 r3 = r2 * r 561 r4 = r2 * r2 562 r5 = r2 * r3 563 taper = c5*r5 + c4*r4 + c3*r3 564 & + c2*r2 + c1*r + c0 565 dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3 566 & + 3.0d0*c3*r2 + 2.0d0*c2*r + c1 567 de = e*dtaper + de*taper 568 e = e * taper 569 end if 570c 571c scale the interaction based on its group membership 572c 573 if (use_group) then 574 e = e * fgrp 575 de = de * fgrp 576 end if 577c 578c compute the force components for this interaction 579c 580 frcx = de * xr * rr1 581 frcy = de * yr * rr1 582 frcz = de * zr * rr1 583c 584c increment the total charge transfer energy and derivatives 585c 586 ect = ect + e 587 dect(1,i) = dect(1,i) - frcx 588 dect(2,i) = dect(2,i) - frcy 589 dect(3,i) = dect(3,i) - frcz 590 dect(1,k) = dect(1,k) + frcx 591 dect(2,k) = dect(2,k) + frcy 592 dect(3,k) = dect(3,k) + frcz 593c 594c increment the internal virial tensor components 595c 596 vxx = xr * frcx 597 vxy = yr * frcx 598 vxz = zr * frcx 599 vyy = yr * frcy 600 vyz = zr * frcy 601 vzz = zr * frcz 602 vir(1,1) = vir(1,1) + vxx 603 vir(2,1) = vir(2,1) + vxy 604 vir(3,1) = vir(3,1) + vxz 605 vir(1,2) = vir(1,2) + vxy 606 vir(2,2) = vir(2,2) + vyy 607 vir(3,2) = vir(3,2) + vyz 608 vir(1,3) = vir(1,3) + vxz 609 vir(2,3) = vir(2,3) + vyz 610 vir(3,3) = vir(3,3) + vzz 611 end if 612 end if 613 end do 614c 615c reset exclusion coefficients for connected atoms 616c 617 do j = 1, n12(i) 618 mscale(i12(j,i)) = 1.0d0 619 end do 620 do j = 1, n13(i) 621 mscale(i13(j,i)) = 1.0d0 622 end do 623 do j = 1, n14(i) 624 mscale(i14(j,i)) = 1.0d0 625 end do 626 do j = 1, n15(i) 627 mscale(i15(j,i)) = 1.0d0 628 end do 629 end do 630c 631c OpenMP directives for the major loop structure 632c 633!$OMP END DO 634!$OMP END PARALLEL 635c 636c perform deallocation of some local arrays 637c 638 deallocate (mscale) 639 return 640 end 641