1c 2c 3c ################################################################ 4c ## COPYRIGHT (C) 1990 by Craig Kundrot & Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################################ 7c 8c ################################################################ 9c ## ## 10c ## subroutine volume -- excluded volume term via Connolly ## 11c ## ## 12c ################################################################ 13c 14c 15c "volume" calculates the excluded volume via the Connolly 16c analytical volume and surface area algorithm 17c 18c 19 subroutine volume (volume_tot,radius,exclude) 20 implicit none 21 real*8 exclude,probe 22 real*8 volume_tot 23 real*8 area_tot 24 real*8 radius(*) 25c 26c 27c make call to the volume and surface area routine 28c 29 probe = 0.0d0 30 call connolly (volume_tot,area_tot,radius,probe,exclude) 31 return 32 end 33c 34c 35c ################################################################ 36c ## ## 37c ## subroutine volume1 -- Cartesian excluded volume derivs ## 38c ## ## 39c ################################################################ 40c 41c 42c "volume1" calculates first derivatives of the total excluded 43c volume with respect to the Cartesian coordinates of each atom 44c 45c literature reference: 46c 47c C. E. Kundrot, J. W. Ponder and F. M. Richards, "Algorithms for 48c Calculating Excluded Volume and Its Derivatives as a Function 49c of Molecular Conformation and Their Use in Energy Minimization", 50c Journal of Computational Chemistry, 12, 402-409 (1991) 51c 52c 53 subroutine volume1 (radius,probe,dex) 54 use atoms 55 use iounit 56 use math 57 use sizes 58 implicit none 59 integer maxcube,maxarc 60 parameter (maxcube=15) 61 parameter (maxarc=1000) 62 integer i,j,k,m 63 integer io,ir,in 64 integer narc,nx,ny,nz 65 integer istart,istop 66 integer jstart,jstop 67 integer kstart,kstop 68 integer mstart,mstop 69 integer isum,icube,itemp 70 integer inov(maxarc) 71 integer, allocatable :: itab(:) 72 integer cube(2,maxcube,maxcube,maxcube) 73 real*8 xmin,ymin,zmin 74 real*8 xmax,ymax,zmax 75 real*8 aa,bb,temp,phi_term 76 real*8 theta1,theta2,dtheta 77 real*8 seg_dx,seg_dy,seg_dz 78 real*8 pre_dx,pre_dy,pre_dz 79 real*8 rinsq,rdiff 80 real*8 rsecn,rsec2n 81 real*8 cosine,ti,tf 82 real*8 alpha,beta 83 real*8 ztop,zstart 84 real*8 ztopshave 85 real*8 phi1,cos_phi1 86 real*8 phi2,cos_phi2 87 real*8 zgrid,pix2 88 real*8 rsec2r,rsecr 89 real*8 rr,rrx2,rrsq 90 real*8 rmax,edge 91 real*8 xr,yr,zr 92 real*8 dist2,vdwsum 93 real*8 probe,zstep 94 real*8 arci(maxarc) 95 real*8 arcf(maxarc) 96 real*8 dx(maxarc) 97 real*8 dy(maxarc) 98 real*8 dsq(maxarc) 99 real*8 d(maxarc) 100 real*8, allocatable :: volrad(:) 101 real*8 radius(*) 102 real*8 dex(3,*) 103 logical, allocatable :: skip(:) 104c 105c 106c fix the stepsize in the z-direction; this value sets 107c the accuracy of the numerical derivatives; zstep=0.06 108c is a good balance between compute time and accuracy 109c 110 zstep = 0.0601d0 111c 112c initialize minimum and maximum ranges of atoms 113c 114 pix2 = 2.0d0 * pi 115 rmax = 0.0d0 116 xmin = x(1) 117 xmax = x(1) 118 ymin = y(1) 119 ymax = y(1) 120 zmin = z(1) 121 zmax = z(1) 122c 123c perform dynamic allocation of some local arrays 124c 125 allocate (itab(n)) 126 allocate (volrad(n)) 127 allocate (skip(n)) 128c 129c assign van der Waals radii to the atoms; note that 130c the radii are incremented by the size of the probe; 131c then get the maximum and minimum ranges of atoms 132c 133 do i = 1, n 134 volrad(i) = radius(i) 135 if (volrad(i) .eq. 0.0d0) then 136 skip(i) = .true. 137 else 138 skip(i) = .false. 139 volrad(i) = volrad(i) + probe 140 if (volrad(i) .gt. rmax) rmax = volrad(i) 141 if (x(i) .lt. xmin) xmin = x(i) 142 if (x(i) .gt. xmax) xmax = x(i) 143 if (y(i) .lt. ymin) ymin = y(i) 144 if (y(i) .gt. ymax) ymax = y(i) 145 if (z(i) .lt. zmin) zmin = z(i) 146 if (z(i) .gt. zmax) zmax = z(i) 147 end if 148 end do 149c 150c load the cubes based on coarse lattice; first of all 151c set edge length to the maximum diameter of any atom 152c 153 edge = 2.0d0 * rmax 154 nx = int((xmax-xmin)/edge) + 1 155 ny = int((ymax-ymin)/edge) + 1 156 nz = int((zmax-zmin)/edge) + 1 157 if (max(nx,ny,nz) .gt. maxcube) then 158 write (iout,10) 159 10 format (/,' VOLUME1 -- Increase the Value of MAXCUBE') 160 call fatal 161 end if 162c 163c initialize the coarse lattice of cubes 164c 165 do i = 1, nx 166 do j = 1, ny 167 do k = 1, nz 168 cube(1,i,j,k) = 0 169 cube(2,i,j,k) = 0 170 end do 171 end do 172 end do 173c 174c find the number of atoms in each cube 175c 176 do m = 1, n 177 if (.not. skip(m)) then 178 i = int((x(m)-xmin)/edge) + 1 179 j = int((y(m)-ymin)/edge) + 1 180 k = int((z(m)-zmin)/edge) + 1 181 cube(1,i,j,k) = cube(1,i,j,k) + 1 182 end if 183 end do 184c 185c determine the highest index in the array "itab" for the 186c atoms that fall into each cube; the first cube that has 187c atoms defines the first index for "itab"; the final index 188c for the atoms in the present cube is the final index of 189c the last cube plus the number of atoms in the present cube 190c 191 isum = 0 192 do i = 1, nx 193 do j = 1, ny 194 do k = 1, nz 195 icube = cube(1,i,j,k) 196 if (icube .ne. 0) then 197 isum = isum + icube 198 cube(2,i,j,k) = isum 199 end if 200 end do 201 end do 202 end do 203c 204c "cube(2,,,)" now contains a pointer to the array "itab" 205c giving the position of the last entry for the list of 206c atoms in that cube of total number equal to "cube(1,,,)" 207c 208 do m = 1, n 209 if (.not. skip(m)) then 210 i = int((x(m)-xmin)/edge) + 1 211 j = int((y(m)-ymin)/edge) + 1 212 k = int((z(m)-zmin)/edge) + 1 213 icube = cube(2,i,j,k) 214 itab(icube) = m 215 cube(2,i,j,k) = icube - 1 216 end if 217 end do 218c 219c set "cube(2,,,)" to be the starting index in "itab" 220c for atom list of that cube; and "cube(1,,,)" to be 221c the stop index 222c 223 isum = 0 224 do i = 1, nx 225 do j = 1, ny 226 do k = 1, nz 227 icube = cube(1,i,j,k) 228 if (icube .ne. 0) then 229 isum = isum + icube 230 cube(1,i,j,k) = isum 231 cube(2,i,j,k) = cube(2,i,j,k) + 1 232 end if 233 end do 234 end do 235 end do 236c 237c process in turn each atom from the coordinate list; 238c first select the potential intersecting atoms 239c 240 do ir = 1, n 241 pre_dx = 0.0d0 242 pre_dy = 0.0d0 243 pre_dz = 0.0d0 244 if (skip(ir)) goto 50 245 rr = volrad(ir) 246 rrx2 = 2.0d0 * rr 247 rrsq = rr * rr 248 xr = x(ir) 249 yr = y(ir) 250 zr = z(ir) 251c 252c find cubes to search for overlaps of current atom 253c 254 istart = int((xr-xmin)/edge) 255 istop = min(istart+2,nx) 256 istart = max(istart,1) 257 jstart = int((yr-ymin)/edge) 258 jstop = min(jstart+2,ny) 259 jstart = max(jstart,1) 260 kstart = int((zr-zmin)/edge) 261 kstop = min(kstart+2,nz) 262 kstart = max(kstart,1) 263c 264c load all overlapping atoms into "inov" 265c 266 io = 0 267 do i = istart, istop 268 do j = jstart, jstop 269 do k = kstart, kstop 270 mstart = cube(2,i,j,k) 271 if (mstart .ne. 0) then 272 mstop = cube(1,i,j,k) 273 do m = mstart, mstop 274 in = itab(m) 275 if (in .ne. ir) then 276 io = io + 1 277 if (io .gt. maxarc) then 278 write (iout,20) 279 20 format (/,' VOLUME1 -- Increase ', 280 & ' the Value of MAXARC') 281 call fatal 282 end if 283 dx(io) = x(in) - xr 284 dy(io) = y(in) - yr 285 dsq(io) = dx(io)**2 + dy(io)**2 286 dist2 = dsq(io) + (z(in)-zr)**2 287 vdwsum = (rr+volrad(in))**2 288 if (dist2.gt.vdwsum .or. dist2.eq.0.0d0) then 289 io = io - 1 290 else 291 d(io) = sqrt(dsq(io)) 292 inov(io) = in 293 end if 294 end if 295 end do 296 end if 297 end do 298 end do 299 end do 300c 301c determine resolution along the z-axis 302c 303 if (io .ne. 0) then 304 ztop = zr + rr 305 ztopshave = ztop - zstep 306 zgrid = zr - rr 307c 308c half of the part not covered by the planes 309c 310 zgrid = zgrid + 0.5d0*(rrx2-(int(rrx2/zstep)*zstep)) 311 zstart = zgrid 312c 313c section atom spheres perpendicular to the z axis 314c 315 do while (zgrid .le. ztop) 316c 317c "rsecr" is radius of circle of intersection 318c of "ir" sphere on the current sphere 319c 320 rsec2r = rrsq - (zgrid-zr)**2 321 if (rsec2r .lt. 0.0d0) rsec2r = 0.000001d0 322 rsecr = sqrt(rsec2r) 323 if (zgrid .ge. ztopshave) then 324 cos_phi1 = 1.0d0 325 phi1 = 0.0d0 326 else 327 cos_phi1 = (zgrid + 0.5d0*zstep - zr) / rr 328 phi1 = acos(cos_phi1) 329 end if 330 if (zgrid .eq. zstart) then 331 cos_phi2 = -1.0d0 332 phi2 = pi 333 else 334 cos_phi2 = (zgrid - 0.5d0*zstep - zr) / rr 335 phi2 = acos(cos_phi2) 336 end if 337c 338c check intersections of neighbor circles 339c 340 narc = 0 341 do k = 1, io 342 in = inov(k) 343 rinsq = volrad(in)**2 344 rsec2n = rinsq - (zgrid-z(in))**2 345 if (rsec2n .gt. 0.0d0) then 346 rsecn = sqrt(rsec2n) 347 if (d(k) .lt. rsecr+rsecn) then 348 rdiff = rsecr - rsecn 349 if (d(k) .le. abs(rdiff)) then 350 if (rdiff .lt. 0.0d0) then 351 narc = 1 352 arci(narc) = 0.0d0 353 arcf(narc) = pix2 354 end if 355 goto 40 356 end if 357 narc = narc + 1 358 if (narc .gt. maxarc) then 359 write (iout,30) 360 30 format (/,' VOLUME1 -- Increase', 361 & ' the Value of MAXARC') 362 call fatal 363 end if 364c 365c initial and final arc endpoints are found for intersection 366c of "ir" circle with another circle contained in same plane; 367c the initial endpoint of the enclosed arc is stored in "arci", 368c the final endpoint in "arcf"; get "cosine" via law of cosines 369c 370 cosine = (dsq(k)+rsec2r-rsec2n) 371 & / (2.0d0*d(k)*rsecr) 372 cosine = min(1.0d0,max(-1.0d0,cosine)) 373c 374c "alpha" is the angle between a line containing either point 375c of intersection and the reference circle center and the 376c line containing both circle centers; "beta" is the angle 377c between the line containing both circle centers and x-axis 378c 379 alpha = acos(cosine) 380 beta = atan2(dy(k),dx(k)) 381 if (dy(k) .lt. 0.0d0) beta = beta + pix2 382 ti = beta - alpha 383 tf = beta + alpha 384 if (ti .lt. 0.0d0) ti = ti + pix2 385 if (tf .gt. pix2) tf = tf - pix2 386 arci(narc) = ti 387c 388c if the arc crosses zero, then it is broken into two segments; 389c the first ends at two pi and the second begins at zero 390c 391 if (tf .lt. ti) then 392 arcf(narc) = pix2 393 narc = narc + 1 394 arci(narc) = 0.0d0 395 end if 396 arcf(narc) = tf 397 40 continue 398 end if 399 end if 400 end do 401c 402c find the pre-area and pre-forces on this section (band), 403c "pre-" means a multiplicative factor is yet to be applied 404c 405 if (narc .eq. 0) then 406 seg_dz = pix2 * (cos_phi1**2 - cos_phi2**2) 407 pre_dz = pre_dz + seg_dz 408 else 409c 410c sort the arc endpoint arrays, each with "narc" entries, 411c in order of increasing values of the arguments in "arci" 412c 413 k = 1 414 do while (k .lt. narc) 415 aa = arci(k) 416 bb = arcf(k) 417 temp = 1000000.0d0 418 do i = k, narc 419 if (arci(i) .le. temp) then 420 temp = arci(i) 421 itemp = i 422 end if 423 end do 424 arci(k) = arci(itemp) 425 arcf(k) = arcf(itemp) 426 arci(itemp) = aa 427 arcf(itemp) = bb 428 k = k + 1 429 end do 430c 431c consolidate arcs by removing overlapping arc endpoints 432c 433 temp = arcf(1) 434 j = 1 435 do k = 2, narc 436 if (temp .lt. arci(k)) then 437 arcf(j) = temp 438 j = j + 1 439 arci(j) = arci(k) 440 temp = arcf(k) 441 else if (temp .lt. arcf(k)) then 442 temp = arcf(k) 443 end if 444 end do 445 arcf(j) = temp 446 narc = j 447 if (narc .eq. 1) then 448 narc = 2 449 arcf(2) = pix2 450 arci(2) = arcf(1) 451 arcf(1) = arci(1) 452 arci(1) = 0.0d0 453 else 454 temp = arci(1) 455 do k = 1, narc-1 456 arci(k) = arcf(k) 457 arcf(k) = arci(k+1) 458 end do 459 if (temp.eq.0.0d0 .and. arcf(narc).eq.pix2) then 460 narc = narc - 1 461 else 462 arci(narc) = arcf(narc) 463 arcf(narc) = temp 464 end if 465 end if 466c 467c compute the numerical pre-derivative values 468c 469 do k = 1, narc 470 theta1 = arci(k) 471 theta2 = arcf(k) 472 if (theta2 .ge. theta1) then 473 dtheta = theta2 - theta1 474 else 475 dtheta = (theta2+pix2) - theta1 476 end if 477 phi_term = phi2 - phi1 - 0.5d0*(sin(2.0d0*phi2) 478 & -sin(2.0d0*phi1)) 479 seg_dx = (sin(theta2)-sin(theta1)) * phi_term 480 seg_dy = (cos(theta1)-cos(theta2)) * phi_term 481 seg_dz = dtheta * (cos_phi1**2 - cos_phi2**2) 482 pre_dx = pre_dx + seg_dx 483 pre_dy = pre_dy + seg_dy 484 pre_dz = pre_dz + seg_dz 485 end do 486 end if 487 zgrid = zgrid + zstep 488 end do 489 end if 490 50 continue 491 dex(1,ir) = 0.5d0 * rrsq * pre_dx 492 dex(2,ir) = 0.5d0 * rrsq * pre_dy 493 dex(3,ir) = 0.5d0 * rrsq * pre_dz 494 end do 495c 496c perform deallocation of some local arrays 497c 498 deallocate (itab) 499 deallocate (volrad) 500 deallocate (skip) 501 return 502 end 503c 504c 505c ################################################################# 506c ## ## 507c ## subroutine volume2 -- Cartesian excluded volume Hessian ## 508c ## ## 509c ################################################################# 510c 511c 512c "volume2" calculates second derivatives of the total excluded 513c volume with respect to the Cartesian coordinates of the atoms 514c 515c literature reference: 516c 517c C. E. Kundrot, J. W. Ponder and F. M. Richards, "Algorithms for 518c Calculating Excluded Volume and Its Derivatives as a Function 519c of Molecular Conformation and Their Use in Energy Minimization", 520c Journal of Computational Chemistry, 12, 402-409 (1991) 521c 522c 523 subroutine volume2 (iatom,radius,probe,xhess,yhess,zhess) 524 use atoms 525 use iounit 526 use math 527 implicit none 528 integer maxarc 529 parameter (maxarc=1000) 530 integer i,j,k,m 531 integer in,iaa,ibb 532 integer iatom,narc 533 integer iblock,itemp 534 integer idtemp,idfirst 535 integer nnear,id(0:2) 536 integer inear(maxarc) 537 integer arciatom(maxarc) 538 integer arcfatom(maxarc) 539 real*8 probe,zstep,xr,yr,zr 540 real*8 ztop,ztopshave,zstart 541 real*8 aa,bb,temp,tempf 542 real*8 phi1,phi2,phiold 543 real*8 theta1,theta2,firsti 544 real*8 zgrid,rsec2r,rsecr 545 real*8 pix2,dist2,rcut2 546 real*8 rr,rrx2,rrsq 547 real*8 alpha,beta,gamma 548 real*8 ti,tf,ri,s2,b,cosine 549 real*8 rinsq,rsecn,rsec2n 550 real*8 cos1,cos2,sin1,sin2 551 real*8 phi_xy,phi_z 552 real*8 delx(2),dely(2),delz(2) 553 real*8 r_s(2),r_s2(2),u(2) 554 real*8 r(0:2),r_r(0:2) 555 real*8 duds(2),dudr(2) 556 real*8 u_term(2) 557 real*8 dfdtheta(3,2) 558 real*8 dthetadx(2,3,0:2) 559 real*8 dalphdx(2,3,0:2) 560 real*8 dbetadx(2,2,0:2) 561 real*8 dudx(2,3,0:2) 562 real*8 dsdx(2,2,0:2) 563 real*8 drdz(2,0:2) 564 real*8 arci(maxarc) 565 real*8 arcf(maxarc) 566 real*8 dx(maxarc) 567 real*8 dy(maxarc) 568 real*8 dsq(maxarc) 569 real*8 d(maxarc) 570 real*8 radius(*) 571 real*8, allocatable :: volrad(:) 572 real*8 xhess(3,*) 573 real*8 yhess(3,*) 574 real*8 zhess(3,*) 575 logical covered 576c 577c 578c fix the stepsize in the z-direction; this value sets 579c the accuracy of the numerical derivatives; zstep=0.06 580c is a good balance between compute time and accuracy 581c 582 zstep = 0.0601d0 583c 584c zero out the Hessian elements for current atom 585c 586 do i = 1, n 587 do j = 1, 3 588 xhess(j,i) = 0.0d0 589 yhess(j,i) = 0.0d0 590 zhess(j,i) = 0.0d0 591 end do 592 end do 593 if (radius(iatom) .eq. 0.0d0) return 594 pix2 = 2.0d0 * pi 595c 596c perform dynamic allocation of some local arrays 597c 598 allocate (volrad(n)) 599c 600c assign van der Waals radii to the atoms; note that 601c the radii are incremented by the size of the probe 602c 603 do i = 1, n 604 volrad(i) = radius(i) 605 if (volrad(i) .ne. 0.0d0) volrad(i) = volrad(i) + probe 606 end do 607c 608c set the radius and coordinates for current atom 609c 610 rr = volrad(iatom) 611 rrx2 = 2.0d0 * rr 612 rrsq = rr**2 613 xr = x(iatom) 614 yr = y(iatom) 615 zr = z(iatom) 616c 617c select potential intersecting atoms 618c 619 nnear = 1 620 do j = 1, n 621 if (j.ne.iatom .and. volrad(j).ne.0.0d0) then 622 dx(nnear) = x(j) - xr 623 dy(nnear) = y(j) - yr 624 dsq(nnear) = dx(nnear)**2 + dy(nnear)**2 625 dist2 = dsq(nnear) + (z(j)-zr)**2 626 rcut2 = (volrad(j) + rr)**2 627 if (dist2 .lt. rcut2) then 628 d(nnear) = sqrt(dsq(nnear)) 629 inear(nnear) = j 630 nnear = nnear + 1 631 if (nnear .gt. maxarc) then 632 write (iout,10) 633 10 format (/,' VOLUME2 -- Increase', 634 & ' the Value of MAXARC') 635 call fatal 636 end if 637 end if 638 end if 639 end do 640 nnear = nnear - 1 641c 642c determine the z resolution 643c 644 if (nnear .ne. 0) then 645 ztop = zr + rr 646 ztopshave = ztop - zstep 647 zgrid = zr - rr 648c 649c half of the part not covered by the planes 650c 651 zgrid = zgrid + (0.5d0*(rrx2-(int(rrx2/zstep)*zstep))) 652 zstart = zgrid 653c 654c section atom spheres perpendicular to the z axis 655c 656 do while (zgrid .le. ztop) 657c 658c "rsecr" is radius of current atom sphere on the z-plane 659c 660 rsec2r = rrsq - (zgrid-zr)**2 661 if (rsec2r .lt. 0.0d0) then 662 rsec2r = 0.000001d0 663 end if 664 rsecr = sqrt(rsec2r) 665 if (zgrid .ge. ztopshave) then 666 phi1 = 0.0d0 667 else 668 phi1 = acos(((zgrid+0.5d0*zstep)-zr) / rr) 669 end if 670 if (zgrid .eq. zstart) then 671 phi2 = pi 672 else 673 phi2 = phiold 674 end if 675c 676c check intersections of neighbor circles 677c 678 k = 0 679 narc = 0 680 covered = .false. 681 do while (.not.covered .and. k.lt.nnear 682 & .and. narc.lt.maxarc) 683 k = k + 1 684 in = inear(k) 685 rinsq = volrad(in)**2 686 rsec2n = rinsq - (zgrid-z(in))**2 687 if (rsec2n .gt. 0.0d0) then 688 rsecn = sqrt(rsec2n) 689 if (d(k) .lt. rsecr+rsecn) then 690 b = rsecr - rsecn 691 if (d(k) .le. abs(b)) then 692 if (b .lt. 0.0d0) then 693 narc = 1 694 arci(narc) = 0.0d0 695 arcf(narc) = pix2 696 arciatom(narc) = in 697 arcfatom(narc) = in 698 covered = .true. 699 end if 700 else 701 narc = narc + 1 702 if (narc .gt. maxarc) then 703 write (iout,20) 704 20 format (/,' VOLUME2 -- Increase', 705 & ' the Value of MAXARC') 706 call fatal 707 else 708c 709c initial and final arc endpoints are found for intersection 710c of "ir" circle with another circle contained in same plane; 711c the initial endpoint of the enclosed arc is stored in "arci", 712c the final endpoint in "arcf"; get "cosine" via law of cosines 713c 714 cosine = (dsq(k)+rsec2r-rsec2n) / 715 & (2.0d0*d(k)*rsecr) 716 cosine = min(1.0d0,max(-1.0d0,cosine)) 717c 718c "alpha" is the angle between a line containing either point 719c of intersection and the reference circle center and the 720c line containing both circle centers; "beta" is the angle 721c between the line containing both circle centers and x-axis 722c 723 alpha = acos(cosine) 724 if (dx(k) .eq. 0.0d0) then 725 gamma = 0.5d0 * pi 726 else 727 gamma = atan(abs(dy(k)/dx(k))) 728 end if 729 if (dy(k) .gt. 0.0d0) then 730 if (dx(k) .gt. 0.0d0) then 731 beta = gamma 732 else 733 beta = pi - gamma 734 end if 735 else 736 if (dx(k) .gt. 0.0d0) then 737 beta = pix2 - gamma 738 else 739 beta = pi + gamma 740 end if 741 end if 742c 743c finally, the arc endpoints 744c 745 ti = beta - alpha 746 tf = beta + alpha 747 if (ti .lt. 0.0d0) ti = ti + pix2 748 if (tf .gt. pix2) tf = tf - pix2 749 arci(narc) = ti 750 arciatom(narc) = in 751 arcfatom(narc) = in 752 if (tf .lt. ti) then 753 arcf(narc) = pix2 754 narc = narc + 1 755 arci(narc) = 0.0d0 756 arciatom(narc) = in 757 arcfatom(narc) = in 758 end if 759 arcf(narc) = tf 760 end if 761 end if 762 end if 763 end if 764 end do 765c 766c find the pre-area and pre-forces on this section (band) 767c through sphere "ir"; the "pre-" means a multiplicative 768c factor is yet to be applied 769c 770 if (narc .ne. 0) then 771c 772c general case; sort arc endpoints 773c 774 k = 1 775 do while (k .lt. narc) 776 aa = arci(k) 777 bb = arcf(k) 778 iaa = arciatom(k) 779 ibb = arcfatom(k) 780 temp = 10000000.0d0 781 do i = k, narc 782 if (arci(i) .le. temp) then 783 temp = arci(i) 784 itemp = i 785 end if 786 end do 787 arci(k) = arci(itemp) 788 arcf(k) = arcf(itemp) 789 arciatom(k) = arciatom(itemp) 790 arcfatom(k) = arcfatom(itemp) 791 arci(itemp) = aa 792 arcf(itemp) = bb 793 arciatom(itemp) = iaa 794 arcfatom(itemp) = ibb 795 k = k + 1 796 end do 797c 798c eliminate overlapping arc endpoints; 799c first, consolidate the occluded arcs 800c 801 m = 1 802 tempf = arcf(1) 803 idtemp = arcfatom(1) 804 do k = 2, narc 805 if (tempf .lt. arci(k)) then 806 arcf(m) = tempf 807 arcfatom(m) = idtemp 808 m = m + 1 809 arci(m) = arci(k) 810 arciatom(m) = arciatom(k) 811 tempf = arcf(k) 812 idtemp = arcfatom(k) 813 else if (tempf .lt. arcf(k)) then 814 tempf = arcf(k) 815 idtemp = arcfatom(k) 816 end if 817 end do 818 arcf(m) = tempf 819 arcfatom(m) = idtemp 820 narc = m 821c 822c change occluded arcs to accessible arcs 823c 824 if (narc .eq. 1) then 825 if (arci(1).eq.0.0d0 .and. arcf(1).eq.pix2) then 826 narc = 0 827 else 828 firsti = arci(1) 829 idfirst = arciatom(1) 830 arci(1) = arcf(1) 831 arciatom(1) = arcfatom(1) 832 arcf(1) = firsti + pix2 833 arcfatom(1) = idfirst 834 end if 835 else 836 firsti = arci(1) 837 idfirst = arciatom(1) 838 do k = 1, narc-1 839 arci(k) = arcf(k) 840 arciatom(k) = arcfatom(k) 841 arcf(k) = arci(k+1) 842 arcfatom(k) = arciatom(k+1) 843 end do 844c 845c check gap between first and last arcs; if the 846c occluded arc crossed zero, then no accessible arc 847c 848 if (firsti.eq.0.0d0 .and. arcf(narc).eq.pix2) then 849 narc = narc - 1 850 else 851 arci(narc) = arcf(narc) 852 arciatom(narc) = arcfatom(narc) 853 arcf(narc) = firsti 854 arcfatom(narc) = idfirst 855 end if 856 end if 857c 858c setup prior to application of chain rule 859c 860 do k = 1, narc 861 ri = sqrt(rrsq - (zgrid-zr)**2) 862 do i = 1, 2 863 if (i .eq. 1) then 864 id(1) = arciatom(k) 865 else 866 id(2) = arcfatom(k) 867 end if 868 delx(i) = x(id(i)) - xr 869 dely(i) = y(id(i)) - yr 870 delz(i) = zgrid - z(id(i)) 871 s2 = delx(i)**2 + dely(i)**2 872 r_s(i) = 1.0d0 / sqrt(s2) 873 r_s2(i) = r_s(i)**2 874 r(i) = sqrt(volrad(id(i))**2 - delz(i)**2) 875 r_r(i) = 1.0d0 / r(i) 876 u(i) = (ri**2+s2-r(i)**2) * (0.5d0*r_s(i)/ri) 877 end do 878c 879c apply the chain rule repeatedly 880c 881 theta1 = arci(k) 882 theta2 = arcf(k) 883 cos1 = cos(theta1) 884 cos2 = cos(theta2) 885 sin1 = sin(theta1) 886 sin2 = sin(theta2) 887 phi_xy = phi2 - phi1 - 0.5d0*(sin(2.0d0*phi2) 888 & -sin(2.0d0*phi1)) 889 phi_z = sin(phi2)**2 - sin(phi1)**2 890 phi_xy = 0.5d0 * rrsq * phi_xy 891 phi_z = 0.5d0 * rrsq * phi_z 892 dfdtheta(1,1) = -cos1 * phi_xy 893 dfdtheta(2,1) = -sin1 * phi_xy 894 dfdtheta(3,1) = -phi_z 895 dfdtheta(1,2) = cos2 * phi_xy 896 dfdtheta(2,2) = sin2 * phi_xy 897 dfdtheta(3,2) = phi_z 898 do i = 1, 2 899 dbetadx(i,1,0) = dely(i) * r_s2(i) 900 dbetadx(i,2,0) = -delx(i) * r_s2(i) 901 dbetadx(i,1,i) = -dbetadx(i,1,0) 902 dbetadx(i,2,i) = -dbetadx(i,2,0) 903 end do 904 do i = 1, 2 905 duds(i) = (1.0d0/ri) - (u(i)*r_s(i)) 906 dsdx(i,1,i) = delx(i) * r_s(i) 907 dsdx(i,2,i) = dely(i) * r_s(i) 908 dsdx(i,1,0) = -dsdx(i,1,i) 909 dsdx(i,2,0) = -dsdx(i,2,i) 910 dudr(i) = -r(i) * r_s(i) / ri 911 drdz(i,i) = delz(i) * r_r(i) 912 drdz(i,0) = -drdz(i,i) 913 end do 914 do m = 0, 2 915 do i = 1, 2 916 dudx(i,1,m) = duds(i) * dsdx(i,1,m) 917 dudx(i,2,m) = duds(i) * dsdx(i,2,m) 918 dudx(i,3,m) = dudr(i) * drdz(i,m) 919 end do 920 end do 921 do i = 1, 2 922 u_term(i) = -1.0d0 / sqrt(1.0d0-u(i)**2) 923 end do 924 do j = 1, 3 925 do m = 0, 2 926 do i = 1, 2 927 dalphdx(i,j,m) = u_term(i) * dudx(i,j,m) 928 end do 929 end do 930 end do 931 do j = 1, 2 932 do m = 0, 2 933 dthetadx(1,j,m) = dbetadx(1,j,m) 934 & + dalphdx(1,j,m) 935 dthetadx(2,j,m) = dbetadx(2,j,m) 936 & - dalphdx(2,j,m) 937 end do 938 end do 939 do m = 0, 2 940 dthetadx(1,3,m) = dalphdx(1,3,m) 941 dthetadx(2,3,m) = -dalphdx(2,3,m) 942 end do 943c 944c partials with respect to coordinates of serial atom id(m) 945c 946 id(0) = iatom 947 do m = 0, 2 948 iblock = id(m) 949 do j = 1, 3 950 xhess(j,iblock) = xhess(j,iblock) 951 & + dfdtheta(1,1)*dthetadx(1,j,m) 952 & + dfdtheta(1,2)*dthetadx(2,j,m) 953 yhess(j,iblock) = yhess(j,iblock) 954 & + dfdtheta(2,1)*dthetadx(1,j,m) 955 & + dfdtheta(2,2)*dthetadx(2,j,m) 956 zhess(j,iblock) = zhess(j,iblock) 957 & + dfdtheta(3,1)*dthetadx(1,j,m) 958 & + dfdtheta(3,2)*dthetadx(2,j,m) 959 end do 960 end do 961 end do 962 end if 963 zgrid = zgrid + zstep 964 phiold = phi1 965 end do 966 end if 967c 968c perform deallocation of some local arrays 969c 970 deallocate (volrad) 971 return 972 end 973