1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1994 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ################################################################# 9c ## ## 10c ## subroutine egauss1 -- Gaussian vdw energy & derivatives ## 11c ## ## 12c ################################################################# 13c 14c 15c "egauss1" calculates the Gaussian expansion van der Waals 16c interaction energy and its first derivatives with respect 17c to Cartesian coordinates 18c 19c 20 subroutine egauss1 21 use limits 22 use warp 23 implicit none 24c 25c 26c choose the method for summing over pairwise interactions 27c 28 if (use_smooth) then 29 call egauss1d 30 else if (use_vlist) then 31 call egauss1c 32 else if (use_lights) then 33 call egauss1b 34 else 35 call egauss1a 36 end if 37 return 38 end 39c 40c 41c ################################################################ 42c ## ## 43c ## subroutine egauss1a -- double loop Gaussian vdw derivs ## 44c ## ## 45c ################################################################ 46c 47c 48c "egauss1a" calculates the Gaussian expansion van der Waals 49c interaction energy and its first derivatives using a pairwise 50c double loop 51c 52c 53 subroutine egauss1a 54 use atomid 55 use atoms 56 use bound 57 use cell 58 use couple 59 use deriv 60 use energi 61 use group 62 use shunt 63 use usage 64 use vdw 65 use vdwpot 66 use virial 67 implicit none 68 integer i,j,k,m 69 integer ii,iv,it 70 integer kk,kv,kt 71 integer, allocatable :: iv14(:) 72 real*8 e,de,rdn 73 real*8 eps,rad2,fgrp 74 real*8 xi,yi,zi 75 real*8 xr,yr,zr 76 real*8 redi,rediv 77 real*8 redk,redkv 78 real*8 dedx,dedy,dedz 79 real*8 expcut,expterm 80 real*8 taper,dtaper 81 real*8 rik,rik2,rik3 82 real*8 rik4,rik5 83 real*8 vxx,vyy,vzz 84 real*8 vyx,vzx,vzy 85 real*8 a(maxgauss) 86 real*8 b(maxgauss) 87 real*8, allocatable :: vscale(:) 88 logical proceed,usei 89 character*6 mode 90c 91c 92c zero out the van der Waals energy and first derivatives 93c 94 ev = 0.0d0 95 do i = 1, n 96 dev(1,i) = 0.0d0 97 dev(2,i) = 0.0d0 98 dev(3,i) = 0.0d0 99 end do 100 if (nvdw .eq. 0) return 101c 102c perform dynamic allocation of some local arrays 103c 104 allocate (iv14(n)) 105 allocate (vscale(n)) 106c 107c set arrays needed to scale connected atom interactions 108c 109 do i = 1, n 110 iv14(i) = 0 111 vscale(i) = 1.0d0 112 end do 113c 114c set the coefficients for the switching function 115c 116 mode = 'VDW' 117 call switch (mode) 118 expcut = -50.0d0 119c 120c apply any reduction factor to the atomic coordinates 121c 122 do k = 1, nvdw 123 i = ivdw(k) 124 iv = ired(i) 125 rdn = kred(i) 126 xred(i) = rdn*(x(i)-x(iv)) + x(iv) 127 yred(i) = rdn*(y(i)-y(iv)) + y(iv) 128 zred(i) = rdn*(z(i)-z(iv)) + z(iv) 129 end do 130c 131c find van der Waals energy and derivatives via double loop 132c 133 do ii = 1, nvdw-1 134 i = ivdw(ii) 135 iv = ired(i) 136 redi = kred(i) 137 rediv = 1.0d0 - redi 138 it = jvdw(i) 139 xi = xred(i) 140 yi = yred(i) 141 zi = zred(i) 142 usei = (use(i) .or. use(iv)) 143c 144c set exclusion coefficients for connected atoms 145c 146 do j = 1, n12(i) 147 vscale(i12(j,i)) = v2scale 148 end do 149 do j = 1, n13(i) 150 vscale(i13(j,i)) = v3scale 151 end do 152 do j = 1, n14(i) 153 vscale(i14(j,i)) = v4scale 154 iv14(i14(j,i)) = i 155 end do 156 do j = 1, n15(i) 157 vscale(i15(j,i)) = v5scale 158 end do 159c 160c decide whether to compute the current interaction 161c 162 do kk = ii+1, nvdw 163 k = ivdw(kk) 164 kv = ired(k) 165 proceed = .true. 166 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 167 if (proceed) proceed = (usei .or. use(k) .or. use(kv)) 168c 169c compute the energy contribution for this interaction 170c 171 if (proceed) then 172 kt = jvdw(k) 173 xr = xi - xred(k) 174 yr = yi - yred(k) 175 zr = zi - zred(k) 176 call image (xr,yr,zr) 177 rik2 = xr*xr + yr*yr + zr*zr 178c 179c check for an interaction distance less than the cutoff 180c 181 if (rik2 .le. off2) then 182 rad2 = radmin(kt,it)**2 183 eps = epsilon(kt,it) 184 if (iv14(k) .eq. i) then 185 rad2 = radmin4(kt,it)**2 186 eps = epsilon4(kt,it) 187 end if 188 eps = eps * vscale(k) 189 do j = 1, ngauss 190 a(j) = igauss(1,j) * eps 191 b(j) = igauss(2,j) / rad2 192 end do 193 e = 0.0d0 194 de = 0.0d0 195 rik = sqrt(rik2) 196 do j = 1, ngauss 197 expterm = -b(j) * rik2 198 if (expterm .gt. expcut) then 199 expterm = a(j)*exp(expterm) 200 e = e + expterm 201 de = de - 2.0d0*b(j)*rik*expterm 202 end if 203 end do 204c 205c use energy switching if near the cutoff distance 206c 207 if (rik2 .gt. cut2) then 208 rik3 = rik2 * rik 209 rik4 = rik2 * rik2 210 rik5 = rik2 * rik3 211 taper = c5*rik5 + c4*rik4 + c3*rik3 212 & + c2*rik2 + c1*rik + c0 213 dtaper = 5.0d0*c5*rik4 + 4.0d0*c4*rik3 214 & + 3.0d0*c3*rik2 + 2.0d0*c2*rik + c1 215 de = e*dtaper + de*taper 216 e = e * taper 217 end if 218c 219c scale the interaction based on its group membership 220c 221 if (use_group) then 222 e = e * fgrp 223 de = de * fgrp 224 end if 225c 226c find the chain rule terms for derivative components 227c 228 de = de / rik 229 dedx = de * xr 230 dedy = de * yr 231 dedz = de * zr 232c 233c increment the total van der Waals energy and derivatives 234c 235 ev = ev + e 236 if (i .eq. iv) then 237 dev(1,i) = dev(1,i) + dedx 238 dev(2,i) = dev(2,i) + dedy 239 dev(3,i) = dev(3,i) + dedz 240 else 241 dev(1,i) = dev(1,i) + dedx*redi 242 dev(2,i) = dev(2,i) + dedy*redi 243 dev(3,i) = dev(3,i) + dedz*redi 244 dev(1,iv) = dev(1,iv) + dedx*rediv 245 dev(2,iv) = dev(2,iv) + dedy*rediv 246 dev(3,iv) = dev(3,iv) + dedz*rediv 247 end if 248 if (k .eq. kv) then 249 dev(1,k) = dev(1,k) - dedx 250 dev(2,k) = dev(2,k) - dedy 251 dev(3,k) = dev(3,k) - dedz 252 else 253 redk = kred(k) 254 redkv = 1.0d0 - redk 255 dev(1,k) = dev(1,k) - dedx*redk 256 dev(2,k) = dev(2,k) - dedy*redk 257 dev(3,k) = dev(3,k) - dedz*redk 258 dev(1,kv) = dev(1,kv) - dedx*redkv 259 dev(2,kv) = dev(2,kv) - dedy*redkv 260 dev(3,kv) = dev(3,kv) - dedz*redkv 261 end if 262c 263c increment the internal virial tensor components 264c 265 vxx = xr * dedx 266 vyx = yr * dedx 267 vzx = zr * dedx 268 vyy = yr * dedy 269 vzy = zr * dedy 270 vzz = zr * dedz 271 vir(1,1) = vir(1,1) + vxx 272 vir(2,1) = vir(2,1) + vyx 273 vir(3,1) = vir(3,1) + vzx 274 vir(1,2) = vir(1,2) + vyx 275 vir(2,2) = vir(2,2) + vyy 276 vir(3,2) = vir(3,2) + vzy 277 vir(1,3) = vir(1,3) + vzx 278 vir(2,3) = vir(2,3) + vzy 279 vir(3,3) = vir(3,3) + vzz 280 end if 281 end if 282 end do 283c 284c reset exclusion coefficients for connected atoms 285c 286 do j = 1, n12(i) 287 vscale(i12(j,i)) = 1.0d0 288 end do 289 do j = 1, n13(i) 290 vscale(i13(j,i)) = 1.0d0 291 end do 292 do j = 1, n14(i) 293 vscale(i14(j,i)) = 1.0d0 294 end do 295 do j = 1, n15(i) 296 vscale(i15(j,i)) = 1.0d0 297 end do 298 end do 299c 300c for periodic boundary conditions with large cutoffs 301c neighbors must be found by the replicates method 302c 303 if (.not. use_replica) return 304c 305c calculate interaction energy with other unit cells 306c 307 do ii = 1, nvdw 308 i = ivdw(ii) 309 iv = ired(i) 310 redi = kred(i) 311 rediv = 1.0d0 - redi 312 it = jvdw(i) 313 xi = xred(i) 314 yi = yred(i) 315 zi = zred(i) 316 usei = (use(i) .or. use(iv)) 317c 318c set exclusion coefficients for connected atoms 319c 320 do j = 1, n12(i) 321 vscale(i12(j,i)) = v2scale 322 end do 323 do j = 1, n13(i) 324 vscale(i13(j,i)) = v3scale 325 end do 326 do j = 1, n14(i) 327 vscale(i14(j,i)) = v4scale 328 iv14(i14(j,i)) = i 329 end do 330 do j = 1, n15(i) 331 vscale(i15(j,i)) = v5scale 332 end do 333c 334c decide whether to compute the current interaction 335c 336 do kk = ii, nvdw 337 k = ivdw(kk) 338 kv = ired(k) 339 proceed = .true. 340 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 341 if (proceed) proceed = (usei .or. use(k) .or. use(kv)) 342c 343c compute the energy contribution for this interaction 344c 345 if (proceed) then 346 kt = jvdw(k) 347 do m = 2, ncell 348 xr = xi - xred(k) 349 yr = yi - yred(k) 350 zr = zi - zred(k) 351 call imager (xr,yr,zr,m) 352 rik2 = xr*xr + yr*yr + zr*zr 353c 354c check for an interaction distance less than the cutoff 355c 356 if (rik2 .le. off2) then 357 rad2 = radmin(kt,it)**2 358 eps = epsilon(kt,it) 359 if (use_polymer) then 360 if (rik2 .le. polycut2) then 361 if (iv14(k) .eq. i) then 362 rad2 = radmin4(kt,it)**2 363 eps = epsilon4(kt,it) 364 end if 365 eps = eps * vscale(k) 366 end if 367 end if 368 do j = 1, ngauss 369 a(j) = igauss(1,j) * eps 370 b(j) = igauss(2,j) / rad2 371 end do 372 e = 0.0d0 373 de = 0.0d0 374 rik = sqrt(rik2) 375 do j = 1, ngauss 376 expterm = -b(j) * rik2 377 if (expterm .gt. expcut) then 378 expterm = a(j)*exp(expterm) 379 e = e + expterm 380 de = de - 2.0d0*b(j)*rik*expterm 381 end if 382 end do 383c 384c use energy switching if near the cutoff distance 385c 386 if (rik2 .gt. cut2) then 387 rik3 = rik2 * rik 388 rik4 = rik2 * rik2 389 rik5 = rik2 * rik3 390 taper = c5*rik5 + c4*rik4 + c3*rik3 391 & + c2*rik2 + c1*rik + c0 392 dtaper = 5.0d0*c5*rik4 + 4.0d0*c4*rik3 393 & + 3.0d0*c3*rik2 + 2.0d0*c2*rik + c1 394 de = e*dtaper + de*taper 395 e = e * taper 396 end if 397c 398c scale the interaction based on its group membership 399c 400 if (use_group) then 401 e = e * fgrp 402 de = de * fgrp 403 end if 404c 405c find the chain rule terms for derivative components 406c 407 de = de / rik 408 dedx = de * xr 409 dedy = de * yr 410 dedz = de * zr 411c 412c increment the total van der Waals energy and derivatives 413c 414 if (i .eq. k) e = 0.5d0 * e 415 ev = ev + e 416 if (i .eq. iv) then 417 dev(1,i) = dev(1,i) + dedx 418 dev(2,i) = dev(2,i) + dedy 419 dev(3,i) = dev(3,i) + dedz 420 else 421 dev(1,i) = dev(1,i) + dedx*redi 422 dev(2,i) = dev(2,i) + dedy*redi 423 dev(3,i) = dev(3,i) + dedz*redi 424 dev(1,iv) = dev(1,iv) + dedx*rediv 425 dev(2,iv) = dev(2,iv) + dedy*rediv 426 dev(3,iv) = dev(3,iv) + dedz*rediv 427 end if 428 if (i .ne. k) then 429 if (k .eq. kv) then 430 dev(1,k) = dev(1,k) - dedx 431 dev(2,k) = dev(2,k) - dedy 432 dev(3,k) = dev(3,k) - dedz 433 else 434 redk = kred(k) 435 redkv = 1.0d0 - redk 436 dev(1,k) = dev(1,k) - dedx*redk 437 dev(2,k) = dev(2,k) - dedy*redk 438 dev(3,k) = dev(3,k) - dedz*redk 439 dev(1,kv) = dev(1,kv) - dedx*redkv 440 dev(2,kv) = dev(2,kv) - dedy*redkv 441 dev(3,kv) = dev(3,kv) - dedz*redkv 442 end if 443 end if 444c 445c increment the internal virial tensor components 446c 447 vxx = xr * dedx 448 vyx = yr * dedx 449 vzx = zr * dedx 450 vyy = yr * dedy 451 vzy = zr * dedy 452 vzz = zr * dedz 453 vir(1,1) = vir(1,1) + vxx 454 vir(2,1) = vir(2,1) + vyx 455 vir(3,1) = vir(3,1) + vzx 456 vir(1,2) = vir(1,2) + vyx 457 vir(2,2) = vir(2,2) + vyy 458 vir(3,2) = vir(3,2) + vzy 459 vir(1,3) = vir(1,3) + vzx 460 vir(2,3) = vir(2,3) + vzy 461 vir(3,3) = vir(3,3) + vzz 462 end if 463 end do 464 end if 465 end do 466c 467c reset exclusion coefficients for connected atoms 468c 469 do j = 1, n12(i) 470 vscale(i12(j,i)) = 1.0d0 471 end do 472 do j = 1, n13(i) 473 vscale(i13(j,i)) = 1.0d0 474 end do 475 do j = 1, n14(i) 476 vscale(i14(j,i)) = 1.0d0 477 end do 478 do j = 1, n15(i) 479 vscale(i15(j,i)) = 1.0d0 480 end do 481 end do 482c 483c perform deallocation of some local arrays 484c 485 deallocate (iv14) 486 deallocate (vscale) 487 return 488 end 489c 490c 491c ############################################################### 492c ## ## 493c ## subroutine egauss1b -- Gaussian vdw derivs via lights ## 494c ## ## 495c ############################################################### 496c 497c 498c "egauss1b" calculates the Gaussian expansion van der Waals 499c energy and its first derivatives with respect to Cartesian 500c coordinates using the method of lights 501c 502c 503 subroutine egauss1b 504 use atomid 505 use atoms 506 use bound 507 use boxes 508 use cell 509 use couple 510 use deriv 511 use energi 512 use group 513 use light 514 use shunt 515 use usage 516 use vdw 517 use vdwpot 518 use virial 519 implicit none 520 integer i,j,k,m 521 integer ii,iv,it 522 integer kk,kv,kt 523 integer kgy,kgz 524 integer start,stop 525 integer, allocatable :: iv14(:) 526 real*8 e,de,rdn 527 real*8 eps,rad2,fgrp 528 real*8 xi,yi,zi 529 real*8 xr,yr,zr 530 real*8 redi,rediv 531 real*8 redk,redkv 532 real*8 dedx,dedy,dedz 533 real*8 expcut,expterm 534 real*8 taper,dtaper 535 real*8 rik,rik2,rik3 536 real*8 rik4,rik5 537 real*8 vxx,vyy,vzz 538 real*8 vyx,vzx,vzy 539 real*8 a(maxgauss) 540 real*8 b(maxgauss) 541 real*8, allocatable :: vscale(:) 542 real*8, allocatable :: xsort(:) 543 real*8, allocatable :: ysort(:) 544 real*8, allocatable :: zsort(:) 545 logical proceed,usei,prime 546 logical unique,repeat 547 character*6 mode 548c 549c 550c zero out the van der Waals energy and first derivatives 551c 552 ev = 0.0d0 553 do i = 1, n 554 dev(1,i) = 0.0d0 555 dev(2,i) = 0.0d0 556 dev(3,i) = 0.0d0 557 end do 558 if (nvdw .eq. 0) return 559c 560c perform dynamic allocation of some local arrays 561c 562 allocate (iv14(n)) 563 allocate (vscale(n)) 564 allocate (xsort(8*n)) 565 allocate (ysort(8*n)) 566 allocate (zsort(8*n)) 567c 568c set arrays needed to scale connected atom interactions 569c 570 do i = 1, n 571 iv14(i) = 0 572 vscale(i) = 1.0d0 573 end do 574c 575c set the coefficients for the switching function 576c 577 mode = 'VDW' 578 call switch (mode) 579 expcut = -50.0d0 580c 581c apply any reduction factor to the atomic coordinates 582c 583 do j = 1, nvdw 584 i = ivdw(j) 585 iv = ired(i) 586 rdn = kred(i) 587 xred(j) = rdn*(x(i)-x(iv)) + x(iv) 588 yred(j) = rdn*(y(i)-y(iv)) + y(iv) 589 zred(j) = rdn*(z(i)-z(iv)) + z(iv) 590 end do 591c 592c transfer the interaction site coordinates to sorting arrays 593c 594 do i = 1, nvdw 595 xsort(i) = xred(i) 596 ysort(i) = yred(i) 597 zsort(i) = zred(i) 598 end do 599c 600c use the method of lights to generate neighbors 601c 602 unique = .true. 603 call lights (off,nvdw,xsort,ysort,zsort,unique) 604c 605c loop over all atoms computing the interactions 606c 607 do ii = 1, nvdw 608 i = ivdw(ii) 609 iv = ired(i) 610 redi = kred(i) 611 rediv = 1.0d0 - redi 612 it = jvdw(i) 613 xi = xsort(rgx(ii)) 614 yi = ysort(rgy(ii)) 615 zi = zsort(rgz(ii)) 616 usei = (use(i) .or. use(iv)) 617c 618c set exclusion coefficients for connected atoms 619c 620 do j = 1, n12(i) 621 vscale(i12(j,i)) = v2scale 622 end do 623 do j = 1, n13(i) 624 vscale(i13(j,i)) = v3scale 625 end do 626 do j = 1, n14(i) 627 vscale(i14(j,i)) = v4scale 628 iv14(i14(j,i)) = i 629 end do 630 do j = 1, n15(i) 631 vscale(i15(j,i)) = v5scale 632 end do 633c 634c loop over method of lights neighbors of current atom 635c 636 if (kbx(ii) .le. kex(ii)) then 637 repeat = .false. 638 start = kbx(ii) + 1 639 stop = kex(ii) 640 else 641 repeat = .true. 642 start = 1 643 stop = kex(ii) 644 end if 645 10 continue 646 do m = start, stop 647 kk = locx(m) 648 kgy = rgy(kk) 649 if (kby(ii) .le. key(ii)) then 650 if (kgy.lt.kby(ii) .or. kgy.gt.key(ii)) goto 20 651 else 652 if (kgy.lt.kby(ii) .and. kgy.gt.key(ii)) goto 20 653 end if 654 kgz = rgz(kk) 655 if (kbz(ii) .le. kez(ii)) then 656 if (kgz.lt.kbz(ii) .or. kgz.gt.kez(ii)) goto 20 657 else 658 if (kgz.lt.kbz(ii) .and. kgz.gt.kez(ii)) goto 20 659 end if 660 k = ivdw(kk-((kk-1)/nvdw)*nvdw) 661 kv = ired(k) 662 prime = (kk .le. nvdw) 663c 664c decide whether to compute the current interaction 665c 666 proceed = .true. 667 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 668 if (proceed) proceed = (usei .or. use(k) .or. use(kv)) 669c 670c compute the energy contribution for this interaction 671c 672 if (proceed) then 673 kt = jvdw(k) 674 xr = xi - xsort(m) 675 yr = yi - ysort(kgy) 676 zr = zi - zsort(kgz) 677 if (use_bounds) then 678 if (abs(xr) .gt. xcell2) xr = xr - sign(xcell,xr) 679 if (abs(yr) .gt. ycell2) yr = yr - sign(ycell,yr) 680 if (abs(zr) .gt. zcell2) zr = zr - sign(zcell,zr) 681 if (monoclinic) then 682 xr = xr + zr*beta_cos 683 zr = zr * beta_sin 684 else if (triclinic) then 685 xr = xr + yr*gamma_cos + zr*beta_cos 686 yr = yr*gamma_sin + zr*beta_term 687 zr = zr * gamma_term 688 end if 689 end if 690 rik2 = xr*xr + yr*yr + zr*zr 691c 692c check for an interaction distance less than the cutoff 693c 694 if (rik2 .le. off2) then 695 rad2 = radmin(kt,it)**2 696 eps = epsilon(kt,it) 697 if (prime) then 698 if (iv14(k) .eq. i) then 699 rad2 = radmin4(kt,it)**2 700 eps = epsilon4(kt,it) 701 end if 702 eps = eps * vscale(k) 703 end if 704 do j = 1, ngauss 705 a(j) = igauss(1,j) * eps 706 b(j) = igauss(2,j) / rad2 707 end do 708 e = 0.0d0 709 de = 0.0d0 710 rik = sqrt(rik2) 711 do j = 1, ngauss 712 expterm = -b(j) * rik2 713 if (expterm .gt. expcut) then 714 expterm = a(j)*exp(expterm) 715 e = e + expterm 716 de = de - 2.0d0*b(j)*rik*expterm 717 end if 718 end do 719c 720c use energy switching if near the cutoff distance 721c 722 if (rik2 .gt. cut2) then 723 rik3 = rik2 * rik 724 rik4 = rik2 * rik2 725 rik5 = rik2 * rik3 726 taper = c5*rik5 + c4*rik4 + c3*rik3 727 & + c2*rik2 + c1*rik + c0 728 dtaper = 5.0d0*c5*rik4 + 4.0d0*c4*rik3 729 & + 3.0d0*c3*rik2 + 2.0d0*c2*rik + c1 730 de = e*dtaper + de*taper 731 e = e * taper 732 end if 733c 734c scale the interaction based on its group membership 735c 736 if (use_group) then 737 e = e * fgrp 738 de = de * fgrp 739 end if 740c 741c find the chain rule terms for derivative components 742c 743 de = de / rik 744 dedx = de * xr 745 dedy = de * yr 746 dedz = de * zr 747c 748c increment the total van der Waals energy and derivatives 749c 750 ev = ev + e 751 if (i .eq. iv) then 752 dev(1,i) = dev(1,i) + dedx 753 dev(2,i) = dev(2,i) + dedy 754 dev(3,i) = dev(3,i) + dedz 755 else 756 dev(1,i) = dev(1,i) + dedx*redi 757 dev(2,i) = dev(2,i) + dedy*redi 758 dev(3,i) = dev(3,i) + dedz*redi 759 dev(1,iv) = dev(1,iv) + dedx*rediv 760 dev(2,iv) = dev(2,iv) + dedy*rediv 761 dev(3,iv) = dev(3,iv) + dedz*rediv 762 end if 763 if (k .eq. kv) then 764 dev(1,k) = dev(1,k) - dedx 765 dev(2,k) = dev(2,k) - dedy 766 dev(3,k) = dev(3,k) - dedz 767 else 768 redk = kred(k) 769 redkv = 1.0d0 - redk 770 dev(1,k) = dev(1,k) - dedx*redk 771 dev(2,k) = dev(2,k) - dedy*redk 772 dev(3,k) = dev(3,k) - dedz*redk 773 dev(1,kv) = dev(1,kv) - dedx*redkv 774 dev(2,kv) = dev(2,kv) - dedy*redkv 775 dev(3,kv) = dev(3,kv) - dedz*redkv 776 end if 777c 778c increment the internal virial tensor components 779c 780 vxx = xr * dedx 781 vyx = yr * dedx 782 vzx = zr * dedx 783 vyy = yr * dedy 784 vzy = zr * dedy 785 vzz = zr * dedz 786 vir(1,1) = vir(1,1) + vxx 787 vir(2,1) = vir(2,1) + vyx 788 vir(3,1) = vir(3,1) + vzx 789 vir(1,2) = vir(1,2) + vyx 790 vir(2,2) = vir(2,2) + vyy 791 vir(3,2) = vir(3,2) + vzy 792 vir(1,3) = vir(1,3) + vzx 793 vir(2,3) = vir(2,3) + vzy 794 vir(3,3) = vir(3,3) + vzz 795 end if 796 end if 797 20 continue 798 end do 799 if (repeat) then 800 repeat = .false. 801 start = kbx(ii) + 1 802 stop = nlight 803 goto 10 804 end if 805c 806c reset exclusion coefficients for connected atoms 807c 808 do j = 1, n12(i) 809 vscale(i12(j,i)) = 1.0d0 810 end do 811 do j = 1, n13(i) 812 vscale(i13(j,i)) = 1.0d0 813 end do 814 do j = 1, n14(i) 815 vscale(i14(j,i)) = 1.0d0 816 end do 817 do j = 1, n15(i) 818 vscale(i15(j,i)) = 1.0d0 819 end do 820 end do 821c 822c perform deallocation of some local arrays 823c 824 deallocate (iv14) 825 deallocate (vscale) 826 deallocate (xsort) 827 deallocate (ysort) 828 deallocate (zsort) 829 return 830 end 831c 832c 833c ############################################################# 834c ## ## 835c ## subroutine egauss1c -- Gaussian vdw derivs via list ## 836c ## ## 837c ############################################################# 838c 839c 840c "egauss1c" calculates the Gaussian expansion van der Waals 841c energy and its first derivatives with respect to Cartesian 842c coordinates using a pairwise neighbor list 843c 844c 845 subroutine egauss1c 846 use atomid 847 use atoms 848 use bound 849 use couple 850 use deriv 851 use energi 852 use group 853 use neigh 854 use shunt 855 use usage 856 use vdw 857 use vdwpot 858 use virial 859 implicit none 860 integer i,j,k 861 integer ii,iv,it 862 integer kk,kv,kt 863 integer, allocatable :: iv14(:) 864 real*8 e,de,rdn 865 real*8 eps,rad2,fgrp 866 real*8 xi,yi,zi 867 real*8 xr,yr,zr 868 real*8 redi,rediv 869 real*8 redk,redkv 870 real*8 dedx,dedy,dedz 871 real*8 expcut,expterm 872 real*8 taper,dtaper 873 real*8 rik,rik2,rik3 874 real*8 rik4,rik5 875 real*8 vxx,vyy,vzz 876 real*8 vyx,vzx,vzy 877 real*8 a(maxgauss) 878 real*8 b(maxgauss) 879 real*8, allocatable :: vscale(:) 880 logical proceed,usei 881 character*6 mode 882c 883c 884c zero out the van der Waals energy and first derivatives 885c 886 ev = 0.0d0 887 do i = 1, n 888 dev(1,i) = 0.0d0 889 dev(2,i) = 0.0d0 890 dev(3,i) = 0.0d0 891 end do 892 if (nvdw .eq. 0) return 893c 894c perform dynamic allocation of some local arrays 895c 896 allocate (iv14(n)) 897 allocate (vscale(n)) 898c 899c set arrays needed to scale connected atom interactions 900c 901 do i = 1, n 902 iv14(i) = 0 903 vscale(i) = 1.0d0 904 end do 905c 906c set the coefficients for the switching function 907c 908 mode = 'VDW' 909 call switch (mode) 910 expcut = -50.0d0 911c 912c apply any reduction factor to the atomic coordinates 913c 914 do k = 1, nvdw 915 i = ivdw(k) 916 iv = ired(i) 917 rdn = kred(i) 918 xred(i) = rdn*(x(i)-x(iv)) + x(iv) 919 yred(i) = rdn*(y(i)-y(iv)) + y(iv) 920 zred(i) = rdn*(z(i)-z(iv)) + z(iv) 921 end do 922c 923c OpenMP directives for the major loop structure 924c 925!$OMP PARALLEL default(private) shared(nvdw,ivdw,jvdw,ired, 926!$OMP& kred,xred,yred,zred,use,nvlst,vlst,n12,n13,n14,n15, 927!$OMP& i12,i13,i14,i15,v2scale,v3scale,v4scale,v5scale,use_group, 928!$OMP& off2,radmin,epsilon,radmin4,epsilon4,ngauss,igauss,expcut, 929!$OMP& cut2,c0,c1,c2,c3,c4,c5) 930!$OMP& firstprivate(vscale,iv14) shared(ev,dev,vir) 931!$OMP DO reduction(+:ev,dev,vir) schedule(guided) 932c 933c find van der Waals energy and derivatives via neighbor list 934c 935 do ii = 1, nvdw 936 i = ivdw(ii) 937 iv = ired(i) 938 redi = kred(i) 939 rediv = 1.0d0 - redi 940 it = jvdw(i) 941 xi = xred(i) 942 yi = yred(i) 943 zi = zred(i) 944 usei = (use(i) .or. use(iv)) 945c 946c set exclusion coefficients for connected atoms 947c 948 do j = 1, n12(i) 949 vscale(i12(j,i)) = v2scale 950 end do 951 do j = 1, n13(i) 952 vscale(i13(j,i)) = v3scale 953 end do 954 do j = 1, n14(i) 955 vscale(i14(j,i)) = v4scale 956 iv14(i14(j,i)) = i 957 end do 958 do j = 1, n15(i) 959 vscale(i15(j,i)) = v5scale 960 end do 961c 962c decide whether to compute the current interaction 963c 964 do kk = 1, nvlst(ii) 965 k = ivdw(vlst(kk,ii)) 966 kv = ired(k) 967 proceed = .true. 968 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 969 if (proceed) proceed = (usei .or. use(k) .or. use(kv)) 970c 971c compute the energy contribution for this interaction 972c 973 if (proceed) then 974 kt = jvdw(k) 975 xr = xi - xred(k) 976 yr = yi - yred(k) 977 zr = zi - zred(k) 978 call image (xr,yr,zr) 979 rik2 = xr*xr + yr*yr + zr*zr 980c 981c check for an interaction distance less than the cutoff 982c 983 if (rik2 .le. off2) then 984 rad2 = radmin(kt,it)**2 985 eps = epsilon(kt,it) 986 if (iv14(k) .eq. i) then 987 rad2 = radmin4(kt,it)**2 988 eps = epsilon4(kt,it) 989 end if 990 eps = eps * vscale(k) 991 do j = 1, ngauss 992 a(j) = igauss(1,j) * eps 993 b(j) = igauss(2,j) / rad2 994 end do 995 e = 0.0d0 996 de = 0.0d0 997 rik = sqrt(rik2) 998 do j = 1, ngauss 999 expterm = -b(j) * rik2 1000 if (expterm .gt. expcut) then 1001 expterm = a(j)*exp(expterm) 1002 e = e + expterm 1003 de = de - 2.0d0*b(j)*rik*expterm 1004 end if 1005 end do 1006c 1007c use energy switching if near the cutoff distance 1008c 1009 if (rik2 .gt. cut2) then 1010 rik3 = rik2 * rik 1011 rik4 = rik2 * rik2 1012 rik5 = rik2 * rik3 1013 taper = c5*rik5 + c4*rik4 + c3*rik3 1014 & + c2*rik2 + c1*rik + c0 1015 dtaper = 5.0d0*c5*rik4 + 4.0d0*c4*rik3 1016 & + 3.0d0*c3*rik2 + 2.0d0*c2*rik + c1 1017 de = e*dtaper + de*taper 1018 e = e * taper 1019 end if 1020c 1021c scale the interaction based on its group membership 1022c 1023 if (use_group) then 1024 e = e * fgrp 1025 de = de * fgrp 1026 end if 1027c 1028c find the chain rule terms for derivative components 1029c 1030 de = de / rik 1031 dedx = de * xr 1032 dedy = de * yr 1033 dedz = de * zr 1034c 1035c increment the total van der Waals energy and derivatives 1036c 1037 ev = ev + e 1038 if (i .eq. iv) then 1039 dev(1,i) = dev(1,i) + dedx 1040 dev(2,i) = dev(2,i) + dedy 1041 dev(3,i) = dev(3,i) + dedz 1042 else 1043 dev(1,i) = dev(1,i) + dedx*redi 1044 dev(2,i) = dev(2,i) + dedy*redi 1045 dev(3,i) = dev(3,i) + dedz*redi 1046 dev(1,iv) = dev(1,iv) + dedx*rediv 1047 dev(2,iv) = dev(2,iv) + dedy*rediv 1048 dev(3,iv) = dev(3,iv) + dedz*rediv 1049 end if 1050 if (k .eq. kv) then 1051 dev(1,k) = dev(1,k) - dedx 1052 dev(2,k) = dev(2,k) - dedy 1053 dev(3,k) = dev(3,k) - dedz 1054 else 1055 redk = kred(k) 1056 redkv = 1.0d0 - redk 1057 dev(1,k) = dev(1,k) - dedx*redk 1058 dev(2,k) = dev(2,k) - dedy*redk 1059 dev(3,k) = dev(3,k) - dedz*redk 1060 dev(1,kv) = dev(1,kv) - dedx*redkv 1061 dev(2,kv) = dev(2,kv) - dedy*redkv 1062 dev(3,kv) = dev(3,kv) - dedz*redkv 1063 end if 1064c 1065c increment the internal virial tensor components 1066c 1067 vxx = xr * dedx 1068 vyx = yr * dedx 1069 vzx = zr * dedx 1070 vyy = yr * dedy 1071 vzy = zr * dedy 1072 vzz = zr * dedz 1073 vir(1,1) = vir(1,1) + vxx 1074 vir(2,1) = vir(2,1) + vyx 1075 vir(3,1) = vir(3,1) + vzx 1076 vir(1,2) = vir(1,2) + vyx 1077 vir(2,2) = vir(2,2) + vyy 1078 vir(3,2) = vir(3,2) + vzy 1079 vir(1,3) = vir(1,3) + vzx 1080 vir(2,3) = vir(2,3) + vzy 1081 vir(3,3) = vir(3,3) + vzz 1082 end if 1083 end if 1084 end do 1085c 1086c reset exclusion coefficients for connected atoms 1087c 1088 do j = 1, n12(i) 1089 vscale(i12(j,i)) = 1.0d0 1090 end do 1091 do j = 1, n13(i) 1092 vscale(i13(j,i)) = 1.0d0 1093 end do 1094 do j = 1, n14(i) 1095 vscale(i14(j,i)) = 1.0d0 1096 end do 1097 do j = 1, n15(i) 1098 vscale(i15(j,i)) = 1.0d0 1099 end do 1100 end do 1101c 1102c OpenMP directives for the major loop structure 1103c 1104!$OMP END DO 1105!$OMP END PARALLEL 1106c 1107c perform deallocation of some local arrays 1108c 1109 deallocate (iv14) 1110 deallocate (vscale) 1111 return 1112 end 1113c 1114c 1115c ################################################################## 1116c ## ## 1117c ## subroutine egauss1d -- Gaussian vdw derivs for smoothing ## 1118c ## ## 1119c ################################################################## 1120c 1121c 1122c "egauss1d" calculates the Gaussian expansion van der Waals 1123c interaction energy and its first derivatives for use with 1124c potential energy smoothing 1125c 1126c 1127 subroutine egauss1d 1128 use atomid 1129 use atoms 1130 use couple 1131 use deriv 1132 use energi 1133 use group 1134 use math 1135 use usage 1136 use vdw 1137 use vdwpot 1138 use virial 1139 use warp 1140 implicit none 1141 integer i,j,k,ii,kk 1142 integer iv,kv,it,kt 1143 integer, allocatable :: iv14(:) 1144 real*8 e,de,rdn 1145 real*8 eps,rad2,fgrp 1146 real*8 xi,yi,zi 1147 real*8 xr,yr,zr 1148 real*8 redi,rediv 1149 real*8 redk,redkv 1150 real*8 dedx,dedy,dedz 1151 real*8 expcut,broot 1152 real*8 erf,term,term2 1153 real*8 expterm,expterm2 1154 real*8 width,wterm 1155 real*8 rik,rik2 1156 real*8 vxx,vyy,vzz 1157 real*8 vyx,vzx,vzy 1158 real*8 t1,t2 1159 real*8 a(maxgauss) 1160 real*8 b(maxgauss) 1161 real*8, allocatable :: vscale(:) 1162 logical proceed,usei 1163 external erf 1164c 1165c 1166c zero out the van der Waals energy and first derivatives 1167c 1168 ev = 0.0d0 1169 do i = 1, n 1170 dev(1,i) = 0.0d0 1171 dev(2,i) = 0.0d0 1172 dev(3,i) = 0.0d0 1173 end do 1174 if (nvdw .eq. 0) return 1175c 1176c perform dynamic allocation of some local arrays 1177c 1178 allocate (iv14(n)) 1179 allocate (vscale(n)) 1180c 1181c set arrays needed to scale connected atom interactions 1182c 1183 do i = 1, n 1184 iv14(i) = 0 1185 vscale(i) = 1.0d0 1186 end do 1187c 1188c set the extent of smoothing to be performed 1189c 1190 expcut = -50.0d0 1191 width = 0.0d0 1192 if (use_dem) then 1193 width = 4.0d0 * diffv * deform 1194 else if (use_gda) then 1195 wterm = (2.0d0/3.0d0) * diffv 1196 else if (use_tophat) then 1197 width = max(diffv*deform,0.0001d0) 1198 end if 1199c 1200c apply any reduction factor to the atomic coordinates 1201c 1202 do k = 1, nvdw 1203 i = ivdw(k) 1204 iv = ired(i) 1205 rdn = kred(i) 1206 xred(i) = rdn*(x(i)-x(iv)) + x(iv) 1207 yred(i) = rdn*(y(i)-y(iv)) + y(iv) 1208 zred(i) = rdn*(z(i)-z(iv)) + z(iv) 1209 end do 1210c 1211c find van der Waals energy and derivatives via double loop 1212c 1213 do ii = 1, nvdw-1 1214 i = ivdw(ii) 1215 iv = ired(i) 1216 redi = kred(i) 1217 rediv = 1.0d0 - redi 1218 it = jvdw(i) 1219 xi = xred(i) 1220 yi = yred(i) 1221 zi = zred(i) 1222 usei = (use(i) .or. use(iv)) 1223c 1224c set exclusion coefficients for connected atoms 1225c 1226 do j = 1, n12(i) 1227 vscale(i12(j,i)) = v2scale 1228 end do 1229 do j = 1, n13(i) 1230 vscale(i13(j,i)) = v3scale 1231 end do 1232 do j = 1, n14(i) 1233 vscale(i14(j,i)) = v4scale 1234 iv14(i14(j,i)) = i 1235 end do 1236 do j = 1, n15(i) 1237 vscale(i15(j,i)) = v5scale 1238 end do 1239c 1240c decide whether to compute the current interaction 1241c 1242 do kk = ii+1, nvdw 1243 k = ivdw(kk) 1244 kv = ired(k) 1245 proceed = .true. 1246 if (use_group) call groups (proceed,fgrp,i,k,0,0,0,0) 1247 if (proceed) proceed = (usei .or. use(k) .or. use(kv)) 1248c 1249c compute the energy contribution for this interaction 1250c 1251 if (proceed) then 1252 kt = jvdw(k) 1253 xr = xi - xred(k) 1254 yr = yi - yred(k) 1255 zr = zi - zred(k) 1256 rik2 = xr*xr + yr*yr + zr*zr 1257c 1258c check for an interaction distance less than the cutoff 1259c 1260 rad2 = radmin(kt,it)**2 1261 eps = epsilon(kt,it) 1262 if (iv14(k) .eq. i) then 1263 rad2 = radmin4(kt,it)**2 1264 eps = epsilon4(kt,it) 1265 end if 1266 eps = eps * vscale(k) 1267 do j = 1, ngauss 1268 a(j) = igauss(1,j) * eps 1269 b(j) = igauss(2,j) / rad2 1270 end do 1271 e = 0.0d0 1272 de = 0.0d0 1273 rik = sqrt(rik2) 1274c 1275c transform the potential function via smoothing 1276c 1277 if (use_tophat) then 1278 rik = sqrt(rik2) 1279 do j = 1, ngauss 1280 broot = sqrt(b(j)) 1281 expterm = -b(j) * (rik+width)**2 1282 if (expterm .gt. expcut) then 1283 expterm = exp(expterm) 1284 else 1285 expterm = 0.0d0 1286 end if 1287 expterm2 = -b(j) * (width-rik)**2 1288 if (expterm2 .gt. expcut) then 1289 expterm2 = exp(expterm2) 1290 else 1291 expterm2 = 0.0d0 1292 end if 1293 term = broot * (expterm-expterm2) 1294 term = term + rootpi*b(j)*rik 1295 & * (erf(broot*(rik+width)) 1296 & +erf(broot*(width-rik))) 1297 e = e + term*a(j)/(b(j)*b(j)*broot) 1298 term = expterm * (2.0d0*rik*b(j)*width+1.0d0) 1299 term2 = expterm2 * (2.0d0*rik*b(j)*width-1.0d0) 1300 de = de + a(j)*(term+term2)/(b(j)*b(j)) 1301 end do 1302 term = 3.0d0 / (8.0d0*rik*width**3) 1303 e = e * term 1304 de = -de * term / rik 1305 else 1306 if (use_gda) width = wterm * (m2(i)+m2(k)) 1307 do j = 1, ngauss 1308 t1 = 1.0d0 + b(j)*width 1309 t2 = sqrt(t1**3) 1310 expterm = -b(j) * rik2 / t1 1311 if (expterm .gt. expcut) then 1312 expterm = (a(j)/t2)*exp(expterm) 1313 e = e + expterm 1314 de = de - (2.0d0*b(j)*rik/t1)*expterm 1315 end if 1316 end do 1317 end if 1318c 1319c scale the interaction based on its group membership 1320c 1321 if (use_group) then 1322 e = e * fgrp 1323 de = de * fgrp 1324 end if 1325c 1326c find the chain rule terms for derivative components 1327c 1328 de = de / rik 1329 dedx = de * xr 1330 dedy = de * yr 1331 dedz = de * zr 1332c 1333c increment the total van der Waals energy and derivatives 1334c 1335 ev = ev + e 1336 if (i .eq. iv) then 1337 dev(1,i) = dev(1,i) + dedx 1338 dev(2,i) = dev(2,i) + dedy 1339 dev(3,i) = dev(3,i) + dedz 1340 else 1341 dev(1,i) = dev(1,i) + dedx*redi 1342 dev(2,i) = dev(2,i) + dedy*redi 1343 dev(3,i) = dev(3,i) + dedz*redi 1344 dev(1,iv) = dev(1,iv) + dedx*rediv 1345 dev(2,iv) = dev(2,iv) + dedy*rediv 1346 dev(3,iv) = dev(3,iv) + dedz*rediv 1347 end if 1348 if (k .eq. kv) then 1349 dev(1,k) = dev(1,k) - dedx 1350 dev(2,k) = dev(2,k) - dedy 1351 dev(3,k) = dev(3,k) - dedz 1352 else 1353 redk = kred(k) 1354 redkv = 1.0d0 - redk 1355 dev(1,k) = dev(1,k) - dedx*redk 1356 dev(2,k) = dev(2,k) - dedy*redk 1357 dev(3,k) = dev(3,k) - dedz*redk 1358 dev(1,kv) = dev(1,kv) - dedx*redkv 1359 dev(2,kv) = dev(2,kv) - dedy*redkv 1360 dev(3,kv) = dev(3,kv) - dedz*redkv 1361 end if 1362c 1363c increment the internal virial tensor components 1364c 1365 vxx = xr * dedx 1366 vyx = yr * dedx 1367 vzx = zr * dedx 1368 vyy = yr * dedy 1369 vzy = zr * dedy 1370 vzz = zr * dedz 1371 vir(1,1) = vir(1,1) + vxx 1372 vir(2,1) = vir(2,1) + vyx 1373 vir(3,1) = vir(3,1) + vzx 1374 vir(1,2) = vir(1,2) + vyx 1375 vir(2,2) = vir(2,2) + vyy 1376 vir(3,2) = vir(3,2) + vzy 1377 vir(1,3) = vir(1,3) + vzx 1378 vir(2,3) = vir(2,3) + vzy 1379 vir(3,3) = vir(3,3) + vzz 1380 end if 1381 end do 1382c 1383c reset exclusion coefficients for connected atoms 1384c 1385 do j = 1, n12(i) 1386 vscale(i12(j,i)) = 1.0d0 1387 end do 1388 do j = 1, n13(i) 1389 vscale(i13(j,i)) = 1.0d0 1390 end do 1391 do j = 1, n14(i) 1392 vscale(i14(j,i)) = 1.0d0 1393 end do 1394 do j = 1, n15(i) 1395 vscale(i15(j,i)) = 1.0d0 1396 end do 1397 end do 1398c 1399c perform deallocation of some local arrays 1400c 1401 deallocate (iv14) 1402 deallocate (vscale) 1403 return 1404 end 1405