1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1992 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################## 9c ## ## 10c ## subroutine embed -- structures via distance geometry ## 11c ## ## 12c ############################################################## 13c 14c 15c "embed" is a distance geometry routine patterned after the 16c ideas of Gordon Crippen, Irwin Kuntz and Tim Havel; it takes 17c as input a set of upper and lower bounds on the interpoint 18c distances, chirality restraints and torsional restraints, 19c and attempts to generate a set of coordinates that satisfy 20c the input bounds and restraints 21c 22c literature references: 23c 24c G. M. Crippen and T. F. Havel, "Distance Geometry and Molecular 25c Conformation", Research Studies Press, Letchworth U.K., 1988, 26c John Wiley and Sons, U.S. distributor 27c 28c T. F. Havel, "An Evaluation of Computational Strategies for 29c Use in the Determination of Protein Structure from Distance 30c Constraints obtained by Nuclear Magnetic Resonance", Progress 31c in Biophysics and Molecular Biology, 56, 43-78 (1991) 32c 33c 34 subroutine embed 35 use atoms 36 use disgeo 37 use files 38 use inform 39 use iounit 40 use minima 41 use refer 42 implicit none 43 integer maxeigen 44 parameter (maxeigen=5) 45 integer i,j,nvar,nstep 46 integer lext,igeo,freeunit 47 integer maxneg,nneg 48 integer maxinner,ninner 49 integer maxouter,nouter 50 real*8 fctval,grdmin 51 real*8 dt,wall,cpu 52 real*8 temp_start 53 real*8 temp_stop 54 real*8 rg,rmsorig 55 real*8 rmsflip,mass 56 real*8 bounds,contact 57 real*8 chiral,torsion 58 real*8 local,locerr 59 real*8 bnderr,vdwerr 60 real*8 chirer,torser 61 real*8 evl(maxeigen) 62 real*8, allocatable :: v(:) 63 real*8, allocatable :: a(:) 64 real*8, allocatable :: evc(:,:) 65 real*8, allocatable :: matrix(:,:) 66 real*8, allocatable :: derivs(:,:) 67 logical done,valid 68 logical exist,info 69 character*7 errtyp,ext 70 character*240 title 71 character*240 geofile 72c 73c 74c perform dynamic allocation of some local arrays 75c 76 allocate (evc(n,maxeigen)) 77 allocate (matrix(n,n)) 78c 79c initialize any chirality restraints, then smooth the 80c bounds via triangle and inverse triangle inequalities; 81c currently these functions are performed by "distgeom" 82c 83c call kchiral 84c if (verbose .and. n.le.130) then 85c title = 'Input Distance Bounds :' 86c call grafic (n,bnd,title) 87c end if 88c call geodesic 89c if (verbose .and. n.le.130)) then 90c title = 'Triangle Smoothed Bounds :' 91c call grafic (n,bnd,title) 92c end if 93c 94c generate a distance matrix between the upper and 95c lower bounds, then convert to a metric matrix 96c 97 maxinner = 3 98 maxouter = 3 99 maxneg = 2 100 nouter = 0 101 valid = .false. 102 do while (.not. valid) 103 ninner = 0 104 done = .false. 105 do while (.not. done) 106 ninner = ninner + 1 107 call dstmat (matrix) 108 call metric (matrix,nneg) 109 if (nneg.le.maxneg .or. ninner.eq.maxinner) done = .true. 110 compact = 0.0d0 111 end do 112 if (verbose .and. nneg.gt.maxneg) then 113 write (iout,10) nneg 114 10 format (/,' EMBED -- Warning, Using Metric Matrix', 115 & ' with',i4,' Negative Distances') 116 end if 117c 118c find the principle components of metric matrix, then 119c generate the trial Cartesian coordinates 120c 121 nouter = nouter + 1 122 call eigen (evl,evc,matrix,valid) 123 call coords (evl,evc) 124 if (nouter.eq.maxouter .and. .not.valid) then 125 valid = .true. 126 if (verbose) then 127 write (iout,20) 128 20 format (/,' EMBED -- Warning, Using Poor Initial', 129 & ' Coordinates') 130 end if 131 end if 132 end do 133c 134c superimpose embedded structure and enantiomer on reference 135c 136 info = verbose 137 verbose = .false. 138 call impose (nref(1),xref,yref,zref,n,x,y,z,rmsorig) 139 if (use_invert) then 140 do i = 1, n 141 x(i) = -x(i) 142 end do 143 call impose (nref(1),xref,yref,zref,n,x,y,z,rmsflip) 144 if (rmsorig .lt. rmsflip) then 145 do i = 1, n 146 x(i) = -x(i) 147 end do 148 call impose (nref(1),xref,yref,zref,n,x,y,z,rmsorig) 149 end if 150 write (iout,30) rmsorig,rmsflip 151 30 format (/,' RMS Superposition for Original and', 152 & ' Enantiomer : ',2f12.4) 153 end if 154 verbose = info 155c 156c compute an index of compaction for the embedded structure 157c 158 call chksize 159c 160c write out the unrefined embedded atomic coordinate set 161c 162 if (debug) then 163 i = 0 164 exist = .true. 165 do while (exist) 166 i = i + 1 167 lext = 3 168 call numeral (i,ext,lext) 169 geofile = filename(1:leng)//'-embed'//'.'//ext(1:lext) 170 inquire (file=geofile,exist=exist) 171 end do 172 igeo = freeunit () 173 open (unit=igeo,file=geofile,status='new') 174 call prtxyz (igeo) 175 close (unit=igeo) 176 title = 'after Embedding :' 177 call fracdist (title) 178 end if 179c 180c use majorization to improve initial embedded coordinates 181c 182 do i = 1, n 183 matrix(i,i) = 0.0d0 184 do j = i+1, n 185 matrix(j,i) = matrix(i,j) 186 end do 187 end do 188 call majorize (matrix) 189c 190c square the bounds for use during structure refinement 191c 192 do i = 1, n 193 do j = 1, n 194 dbnd(j,i) = dbnd(j,i)**2 195 end do 196 end do 197c 198c perform dynamic allocation of some local arrays 199c 200 if (use_anneal) then 201 nvar = 3 * n 202 allocate (v(nvar)) 203 allocate (a(nvar)) 204 end if 205c 206c minimize the error function via simulated annealing 207c 208 if (verbose) call settime 209 if (use_anneal) then 210 iprint = 0 211 if (verbose) iprint = 10 212 grdmin = 1.0d0 213 mass = 10000.0d0 214 do i = 1, nvar 215 v(i) = 0.0d0 216 a(i) = 0.0d0 217 end do 218 errtyp = 'FINAL' 219 call refine (errtyp,fctval,grdmin) 220 nstep = 1000 221 dt = 0.04d0 222 temp_start = 200.0d0 223 temp_stop = 200.0d0 224 call explore (errtyp,nstep,dt,mass,temp_start,temp_stop,v,a) 225 nstep = 10000 226 dt = 0.2d0 227 temp_start = 200.0d0 228 temp_stop = 0.0d0 229 call explore (errtyp,nstep,dt,mass,temp_start,temp_stop,v,a) 230 grdmin = 0.01d0 231 call refine (errtyp,fctval,grdmin) 232c 233c minimize the error function via nonlinear optimization 234c 235 else 236 iprint = 0 237 if (verbose) iprint = 10 238 grdmin = 0.01d0 239 errtyp = 'INITIAL' 240 call refine (errtyp,fctval,grdmin) 241 errtyp = 'MIDDLE' 242 call refine (errtyp,fctval,grdmin) 243 errtyp = 'FINAL' 244 call refine (errtyp,fctval,grdmin) 245 end if 246 if (verbose) then 247 call gettime (wall,cpu) 248 write (iout,40) wall 249 40 format (/,' Time Required for Refinement :',10x,f12.2, 250 & ' seconds') 251 end if 252c 253c perform deallocation of some local arrays 254c 255 if (use_anneal) then 256 deallocate (v) 257 deallocate (a) 258 end if 259c 260c perform dynamic allocation of some local arrays 261c 262 allocate (derivs(3,n)) 263c 264c print the final error function and its components 265c 266 bounds = bnderr (derivs) 267 contact = vdwerr (derivs) 268 local = locerr (derivs) 269 chiral = chirer (derivs) 270 torsion = torser (derivs) 271 write (iout,50) fctval,bounds,contact,local,chiral,torsion 272 50 format (/,' Results of Distance Geometry Protocol :', 273 & //,' Final Error Function Value :',10x,f16.4, 274 & //,' Distance Restraint Error :',12x,f16.4, 275 & /,' Hard Sphere Contact Error :',11x,f16.4, 276 & /,' Local Geometry Error :',16x,f16.4, 277 & /,' Chirality-Planarity Error :',11x,f16.4, 278 & /,' Torsional Restraint Error :',11x,f16.4) 279c 280c take the root of the currently squared distance bounds 281c 282 do i = 1, n 283 do j = 1, n 284 dbnd(j,i) = sqrt(dbnd(j,i)) 285 end do 286 end do 287c 288c print the final rms deviations and radius of gyration 289c 290 title = 'after Refinement :' 291 call rmserror (title) 292 call gyrate (rg) 293 write (iout,60) rg 294 60 format (/,' Radius of Gyration after Refinement :',6x,f16.4) 295 if (verbose .and. n.le.130) call dmdump (matrix) 296c 297c print the normalized fractional distance distribution 298c 299 if (debug) then 300 title = 'after Refinement :' 301 call fracdist (title) 302 end if 303c 304c perform deallocation of some local arrays 305c 306 deallocate (evc) 307 deallocate (matrix) 308 deallocate (derivs) 309 return 310 end 311c 312c 313c ############################################################## 314c ## ## 315c ## subroutine kchiral -- chirality restraint assignment ## 316c ## ## 317c ############################################################## 318c 319c 320c "kchiral" determines the target value for each chirality 321c and planarity restraint as the signed volume of the 322c parallelpiped spanned by vectors from a common atom to 323c each of three other atoms 324c 325c 326 subroutine kchiral 327 use atoms 328 use inform 329 use iounit 330 use restrn 331 implicit none 332 integer i,j,ia,ib,ic,id 333 real*8 xad,yad,zad 334 real*8 xbd,ybd,zbd 335 real*8 xcd,ycd,zcd 336 real*8 c1,c2,c3 337c 338c 339c compute the signed volume of each parallelpiped; 340c if the defining atoms almost lie in a plane, then 341c set the signed volume to exactly zero 342c 343 do i = 1, nchir 344 ia = ichir(1,i) 345 ib = ichir(2,i) 346 ic = ichir(3,i) 347 id = ichir(4,i) 348 xad = x(ia) - x(id) 349 yad = y(ia) - y(id) 350 zad = z(ia) - z(id) 351 xbd = x(ib) - x(id) 352 ybd = y(ib) - y(id) 353 zbd = z(ib) - z(id) 354 xcd = x(ic) - x(id) 355 ycd = y(ic) - y(id) 356 zcd = z(ic) - z(id) 357 c1 = ybd*zcd - zbd*ycd 358 c2 = ycd*zad - zcd*yad 359 c3 = yad*zbd - zad*ybd 360 chir(1,i) = 0.1d0 361 chir(2,i) = xad*c1 + xbd*c2 + xcd*c3 362 if (abs(chir(2,i)) .lt. 1.0d0) chir(2,i) = 0.0d0 363 chir(3,i) = chir(2,i) 364 end do 365c 366c print out the results for each restraint 367c 368 if (verbose) then 369 if (nchir .ne. 0) then 370 write (iout,10) 371 10 format (/,' Chirality and Planarity Constraints :') 372 write (iout,20) 373 20 format (/,18x,'Atom Numbers',12x,'Signed Volume',/) 374 end if 375 do i = 1, nchir 376 write (iout,30) i,(ichir(j,i),j=1,4),chir(2,i) 377 30 format (i6,5x,4i6,5x,f12.4) 378 end do 379 end if 380 return 381 end 382c 383c 384c ############################################################## 385c ## ## 386c ## subroutine triangle -- triangle inequality smoothing ## 387c ## ## 388c ############################################################## 389c 390c 391c "triangle" smooths the upper and lower distance bounds via 392c the triangle inequality using a full-matrix variant of the 393c Floyd-Warshall shortest path algorithm; this routine is 394c usually much slower than the sparse matrix shortest path 395c methods in "geodesic" and "trifix", and should be used only 396c for comparison with answers generated by those routines 397c 398c literature reference: 399c 400c A. W. M. Dress and T. F. Havel, "Shortest-Path Problems and 401c Molecular Conformation", Discrete Applied Mathematics, 19, 402c 129-144 (1988) 403c 404c 405 subroutine triangle 406 use atoms 407 use disgeo 408 use iounit 409 implicit none 410 integer i,j,k 411 integer ik1,ik2 412 integer jk1,jk2 413 real*8 eps 414 real*8 lij,lik,ljk 415 real*8 uij,uik,ujk 416c 417c 418c use full-matrix algorithm to smooth upper and lower bounds 419c 420 eps = 1.0d-10 421 do k = 1, n 422 do i = 1, n-1 423 ik1 = min(i,k) 424 ik2 = max(i,k) 425 lik = dbnd(ik2,ik1) 426 uik = dbnd(ik1,ik2) 427 do j = i+1, n 428 lij = dbnd(j,i) 429 uij = dbnd(i,j) 430 jk1 = min(j,k) 431 jk2 = max(j,k) 432 ljk = dbnd(jk2,jk1) 433 ujk = dbnd(jk1,jk2) 434 lij = max(lij,lik-ujk,ljk-uik) 435 uij = min(uij,uik+ujk) 436 if (lij-uij .gt. eps) then 437 write (iout,10) 438 10 format (/,' TRIANGLE -- Inconsistent Bounds;', 439 & ' Geometrically Impossible') 440 write (iout,20) i,j,lij,uij 441 20 format (/,' Error at :',6x,2i6,3x,2f9.4) 442 write (iout,30) i,k,lik,uik,j,k,ljk,ujk 443 30 format (/,' Traced to :',5x,2i6,3x,2f9.4, 444 & /,17x,2i6,3x,2f9.4) 445 call fatal 446 end if 447 if (lij-dbnd(j,i) .gt. eps) then 448 write (iout,40) i,j,dbnd(j,i),lij 449 40 format (' TRIANGLE -- Altered Lower Bound at', 450 & 2x,2i6,3x,f9.4,' -->',f9.4) 451 end if 452 if (dbnd(i,j)-uij .gt. eps) then 453 write (iout,50) i,j,dbnd(i,j),uij 454 50 format (' TRIANGLE -- Altered Upper Bound at', 455 & 2x,2i6,3x,f9.4,' -->',f9.4) 456 end if 457 dbnd(j,i) = lij 458 dbnd(i,j) = uij 459 end do 460 end do 461 end do 462 return 463 end 464c 465c 466c ################################################################# 467c ## ## 468c ## subroutine geodesic -- sparse matrix triangle smoothing ## 469c ## ## 470c ################################################################# 471c 472c 473c "geodesic" smooths the upper and lower distance bounds via 474c the triangle inequality using a sparse matrix version of a 475c shortest path algorithm 476c 477c literature reference: 478c 479c G. M. Crippen and T. F. Havel, "Distance Geometry and Molecular 480c Conformation", Research Studies Press, Letchworth U.K., 1988, 481c John Wiley and Sons, U.S. distributor, see section 6-2 482c 483c 484 subroutine geodesic 485 use atoms 486 use disgeo 487 use restrn 488 implicit none 489 integer i,j,k,nlist 490 integer, allocatable :: list(:) 491 integer, allocatable :: key(:) 492 integer, allocatable :: start(:) 493 integer, allocatable :: stop(:) 494 real*8, allocatable :: upper(:) 495 real*8, allocatable :: lower(:) 496c 497c 498c perform dynamic allocation of some local arrays 499c 500 nlist = 2 * ndfix 501 allocate (list(nlist)) 502 allocate (key(nlist)) 503 allocate (start(n)) 504 allocate (stop(n)) 505 allocate (upper(n)) 506 allocate (lower(n)) 507c 508c build an indexed list of atoms in distance restraints 509c 510 do i = 1, n 511 start(i) = 0 512 stop(i) = -1 513 end do 514 do i = 1, ndfix 515 list(i) = idfix(1,i) 516 list(i+ndfix) = idfix(2,i) 517 end do 518 call sort3 (nlist,list,key) 519 j = -1 520 do i = 1, nlist 521 k = list(i) 522 if (k .ne. j) then 523 start(k) = i 524 j = k 525 end if 526 end do 527 j = -1 528 do i = nlist, 1, -1 529 k = list(i) 530 if (k .ne. j) then 531 stop(k) = i 532 j = k 533 end if 534 end do 535 do i = 1, nlist 536 k = key(i) 537 if (k .le. ndfix) then 538 list(i) = idfix(2,k) 539 else 540 list(i) = idfix(1,k-ndfix) 541 end if 542 end do 543c 544c triangle smooth bounds via sparse shortest path method 545c 546 do i = 1, n 547 call minpath (i,upper,lower,start,stop,list) 548 do j = i+1, n 549 dbnd(i,j) = upper(j) 550 dbnd(j,i) = max(lower(j),dbnd(j,i)) 551 end do 552 end do 553c 554c perform deallocation of some local arrays 555c 556 deallocate (list) 557 deallocate (key) 558 deallocate (start) 559 deallocate (stop) 560 deallocate (upper) 561 deallocate (lower) 562 return 563 end 564c 565c 566c ################################################################ 567c ## ## 568c ## subroutine minpath -- triangle smoothed bounds to atom ## 569c ## ## 570c ################################################################ 571c 572c 573c "minpath" is a routine for finding the triangle smoothed upper 574c and lower bounds of each atom to a specified root atom using a 575c sparse variant of the Bellman-Ford shortest path algorithm 576c 577c literature reference: 578c 579c D. P. Bertsekas, "A Simple and Fast Label Correcting Algorithm 580c for Shortest Paths", Networks, 23, 703-709 (1993) 581c 582c 583 subroutine minpath (root,upper,lower,start,stop,list) 584 use atoms 585 use couple 586 use disgeo 587 implicit none 588 integer i,j,k 589 integer narc,root 590 integer head,tail 591 integer, allocatable :: iarc(:) 592 integer, allocatable :: queue(:) 593 integer start(*) 594 integer stop(*) 595 integer list(*) 596 real*8 big,small 597 real*8 upper(*) 598 real*8 lower(*) 599 logical enter 600 logical, allocatable :: queued(:) 601c 602c 603c perform dynamic allocation of some local arrays 604c 605 allocate (iarc(n)) 606 allocate (queue(n)) 607 allocate (queued(n)) 608c 609c initialize candidate atom queue and the path lengths 610c 611 do i = 1, n 612 queued(i) = .false. 613 upper(i) = 1000000.0d0 614 lower(i) = 0.0d0 615 end do 616c 617c put the root atom into the queue of candidate atoms 618c 619 head = root 620 tail = root 621 queue(root) = 0 622 queued(root) = .true. 623 upper(root) = 0.0d0 624c 625c get the next candidate atom from head of queue 626c 627 do while (head .ne. 0) 628 j = head 629 queued(j) = .false. 630 head = queue(head) 631c 632c make a list of arcs to the current candidate atom 633c 634 narc = 0 635 do i = 1, n12(j) 636 k = i12(i,j) 637 if (k .ne. root) then 638 narc = narc + 1 639 iarc(narc) = k 640 end if 641 end do 642 do i = 1, n13(j) 643 k = i13(i,j) 644 if (k .ne. root) then 645 narc = narc + 1 646 iarc(narc) = k 647 end if 648 end do 649 do i = 1, n14(j) 650 k = i14(i,j) 651 if (k .ne. root) then 652 narc = narc + 1 653 iarc(narc) = k 654 end if 655 end do 656 do i = start(j), stop(j) 657 k = list(i) 658 if (k .ne. root) then 659 narc = narc + 1 660 iarc(narc) = k 661 end if 662 end do 663c 664c check each arc for alteration of the path length bounds 665c 666 do i = 1, narc 667 k = iarc(i) 668 if (k .lt. j) then 669 big = upper(j) + dbnd(k,j) 670 small = max(dbnd(j,k)-upper(j),lower(j)-dbnd(k,j)) 671 else 672 big = upper(j) + dbnd(j,k) 673 small = max(dbnd(k,j)-upper(j),lower(j)-dbnd(j,k)) 674 end if 675 enter = .false. 676 if (upper(k) .gt. big) then 677 upper(k) = big 678 if (.not. queued(k)) enter = .true. 679 end if 680 if (lower(k) .lt. small) then 681 lower(k) = small 682 if (.not. queued(k)) enter = .true. 683 end if 684c 685c enter a new candidate atom at the tail of the queue 686c 687 if (enter) then 688 queued(k) = .true. 689 if (head .eq. 0) then 690 head = k 691 tail = k 692 queue(k) = 0 693 else 694 queue(tail) = k 695 queue(k) = 0 696 tail = k 697 end if 698 end if 699 end do 700 end do 701c 702c perform deallocation of some local arrays 703c 704 deallocate (iarc) 705 deallocate (queue) 706 deallocate (queued) 707 return 708 end 709c 710c 711c ################################################################ 712c ## ## 713c ## subroutine trifix -- update triangle inequality bounds ## 714c ## ## 715c ################################################################ 716c 717c 718c "trifix" rebuilds both the upper and lower distance bound 719c matrices following tightening of one or both of the bounds 720c between a specified pair of atoms, "p" and "q", using a 721c modification of Murchland's shortest path update algorithm 722c 723c literature references: 724c 725c P. A. Steenbrink, "Optimization of Transport Networks", John 726c Wiley and Sons, Bristol, 1974; see section 7.7 727c 728c R. Dionne, "Etude et Extension d'un Algorithme de Murchland", 729c Infor, 16, 132-146 (1978) 730c 731c 732 subroutine trifix (p,q) 733 use atoms 734 use disgeo 735 use inform 736 use iounit 737 implicit none 738 integer i,k,p,q 739 integer ip,iq,np,nq 740 integer, allocatable :: pt(:) 741 integer, allocatable :: qt(:) 742 real*8 eps,ipmin,ipmax 743 real*8 iqmin,iqmax 744 real*8, allocatable :: pmin(:) 745 real*8, allocatable :: pmax(:) 746 real*8, allocatable :: qmin(:) 747 real*8, allocatable :: qmax(:) 748 logical, allocatable :: pun(:) 749 logical, allocatable :: qun(:) 750c 751c 752c perform dynamic allocation of some local arrays 753c 754 allocate (pt(n)) 755 allocate (qt(n)) 756 allocate (pmin(n)) 757 allocate (pmax(n)) 758 allocate (qmin(n)) 759 allocate (qmax(n)) 760 allocate (pun(n)) 761 allocate (qun(n)) 762c 763c initialize the set of nodes that may have changed bounds 764c 765 eps = 1.0d-10 766 np = 0 767 nq = 0 768 do i = 1, n 769 pun(i) = .true. 770 qun(i) = .true. 771 end do 772c 773c store the upper and lower bounds to "p" and "q" 774c 775 do i = 1, p 776 pmin(i) = dbnd(p,i) 777 pmax(i) = dbnd(i,p) 778 end do 779 do i = p+1, n 780 pmin(i) = dbnd(i,p) 781 pmax(i) = dbnd(p,i) 782 end do 783 do i = 1, q 784 qmin(i) = dbnd(q,i) 785 qmax(i) = dbnd(i,q) 786 end do 787 do i = q+1, n 788 qmin(i) = dbnd(i,q) 789 qmax(i) = dbnd(q,i) 790 end do 791c 792c check for changes in the upper bounds to "p" and "q" 793c 794 do i = 1, n 795 ipmax = qmax(p) + qmax(i) 796 if (pmax(i) .gt. ipmax+eps) then 797 np = np + 1 798 pt(np) = i 799 pmax(i) = ipmax 800 pun(i) = .false. 801 end if 802 iqmax = pmax(q) + pmax(i) 803 if (qmax(i) .gt. iqmax+eps) then 804 nq = nq + 1 805 qt(nq) = i 806 qmax(i) = iqmax 807 qun(i) = .false. 808 end if 809 end do 810c 811c for node pairs whose bounds to "p" and "q" have changed, 812c make any needed changes to upper bound of the pair 813c 814 do ip = 1, np 815 i = pt(ip) 816 ipmax = pmax(i) 817 do iq = 1, nq 818 k = qt(iq) 819 if (i .lt. k) then 820 dbnd(i,k) = min(dbnd(i,k),ipmax+pmax(k)) 821 else 822 dbnd(k,i) = min(dbnd(k,i),ipmax+pmax(k)) 823 end if 824 end do 825 end do 826c 827c check for changes in the lower bounds to "p" and "q" 828c 829 do i = 1, n 830 ipmin = max(qmin(p)-qmax(i),qmin(i)-qmax(p)) 831 if (pmin(i) .lt. ipmin-eps) then 832 if (pun(i)) then 833 np = np + 1 834 pt(np) = i 835 end if 836 pmin(i) = ipmin 837 end if 838 iqmin = max(pmin(q)-pmax(i),pmin(i)-pmax(q)) 839 if (qmin(i) .lt. iqmin-eps) then 840 if (qun(i)) then 841 nq = nq + 1 842 qt(nq) = i 843 end if 844 qmin(i) = iqmin 845 end if 846 end do 847c 848c for node pairs whose bounds to "p" and "q" have changed, 849c make any needed changes to lower bound of the pair 850c 851 do ip = 1, np 852 i = pt(ip) 853 ipmin = pmin(i) 854 ipmax = pmax(i) 855 do iq = 1, nq 856 k = qt(iq) 857 if (i .lt. k) then 858 dbnd(k,i) = max(dbnd(k,i),ipmin-pmax(k),pmin(k)-ipmax) 859 else 860 dbnd(i,k) = max(dbnd(i,k),ipmin-pmax(k),pmin(k)-ipmax) 861 end if 862 end do 863 end do 864c 865c update the upper and lower bounds to "p" and "q" 866c 867 do i = 1, p 868 dbnd(p,i) = pmin(i) 869 dbnd(i,p) = pmax(i) 870 end do 871 do i = p+1, n 872 dbnd(i,p) = pmin(i) 873 dbnd(p,i) = pmax(i) 874 end do 875 do i = 1, q 876 dbnd(q,i) = qmin(i) 877 dbnd(i,q) = qmax(i) 878 end do 879 do i = q+1, n 880 dbnd(i,q) = qmin(i) 881 dbnd(q,i) = qmax(i) 882 end do 883c 884c output the atoms updated and amount of work required 885c 886c if (debug) then 887c write (iout,10) p,q,np*nq 888c 10 format (' TRIFIX -- Bounds Update for Atoms',2i6, 889c & ' with',i8,' Searches') 890c end if 891c 892c perform deallocation of some local arrays 893c 894 deallocate (pt) 895 deallocate (qt) 896 deallocate (pmin) 897 deallocate (pmax) 898 deallocate (qmin) 899 deallocate (qmax) 900 deallocate (pun) 901 deallocate (qun) 902 return 903 end 904c 905c 906c ################################################################ 907c ## ## 908c ## subroutine grafic -- schematic graphical matrix output ## 909c ## ## 910c ################################################################ 911c 912c 913c "grafic" outputs the upper & lower triangles and diagonal 914c of a square matrix in a schematic form for visual inspection 915c 916c 917 subroutine grafic (n,a,title) 918 use iounit 919 implicit none 920 integer i,j,k,m,n 921 integer maxj,nrow,ndash 922 integer minrow,maxrow 923 integer trimtext 924 real*8 big,v 925 real*8 amin,dmin,bmin 926 real*8 amax,dmax,bmax 927 real*8 rcl,scl,tcl 928 real*8 ca,cb,cc,cd 929 real*8 cw,cx,cy,cz 930 real*8 a(n,*) 931 character*1 dash 932 character*1 ta,tb,tc,td,te 933 character*1 digit(0:9) 934 character*1 symbol(130) 935 character*240 title 936 data dash / '-' / 937 data ta,tb,tc,td,te / ' ','.','+','X','#' / 938 data digit / '0','1','2','3','4','5','6','7','8','9' / 939c 940c 941c set bounds of length of print row and write the header 942c 943 minrow = 54 944 maxrow = 130 945 ndash = min(max(n,minrow),maxrow) 946 write (iout,10) (dash,i=1,ndash) 947 10 format (/,1x,130a1) 948 write (iout,20) title(1:trimtext(title)) 949 20 format (/,1x,a) 950c 951c find the maximum and minimum elements of the matrix 952c 953 big = 1.0d6 954 dmax = -big 955 dmin = big 956 amax = -big 957 amin = big 958 bmax = -big 959 bmin = big 960 do i = 1, n 961 if (a(i,i) .gt. dmax) dmax = a(i,i) 962 if (a(i,i) .lt. dmin) dmin = a(i,i) 963 do j = 1, i-1 964 if (a(j,i) .gt. amax) amax = a(j,i) 965 if (a(j,i) .lt. amin) amin = a(j,i) 966 if (a(i,j) .gt. bmax) bmax = a(i,j) 967 if (a(i,j) .lt. bmin) bmin = a(i,j) 968 end do 969 end do 970 write (iout,30) amin,amax,dmin,dmax,bmin,bmax 971 30 format (/,' Range of Above Diag Elements : ',f13.4,' to',f13.4, 972 & /,' Range of Diagonal Elements : ',f13.4,' to',f13.4, 973 & /,' Range of Below Diag Elements : ',f13.4,' to',f13.4) 974c 975c now, print out the graphical representation 976c 977 write (iout,40) 978 40 format (/,' Symbol Magnitude Ordering :',14x, 979 & '# > X > + > . > '' ''',/) 980 rcl = (bmax-bmin) / 5.0d0 981 scl = (amax-amin) / 5.0d0 982 tcl = (dmax-dmin) / 9.0d0 983 if (rcl .eq. 0.0d0) rcl = 1.0d0 984 if (scl .eq. 0.0d0) scl = 1.0d0 985 if (tcl .eq. 0.0d0) tcl = 1.0d0 986 ca = amin + scl 987 cb = ca + scl 988 cc = cb + scl 989 cd = cc + scl 990 cw = bmin + rcl 991 cx = cw + rcl 992 cy = cx + rcl 993 cz = cy + rcl 994 do j = 1, n, maxrow 995 maxj = j + maxrow - 1 996 if (maxj .gt. n) maxj = n 997 nrow = maxj - j + 1 998 do i = 1, n 999 do k = j, maxj 1000 m = k - j + 1 1001 if (k .lt. i) then 1002 v = abs(a(i,k)) 1003 if (v .le. cw) then 1004 symbol(m) = ta 1005 else if (v .le. cx) then 1006 symbol(m) = tb 1007 else if (v .le. cy) then 1008 symbol(m) = tc 1009 else if (v .le. cz) then 1010 symbol(m) = td 1011 else 1012 symbol(m) = te 1013 end if 1014 else if (k .eq. i) then 1015 symbol(m) = digit(nint((a(i,i)-dmin)/tcl)) 1016 else if (k .gt. i) then 1017 v = abs(a(i,k)) 1018 if (v .le. ca) then 1019 symbol(m) = ta 1020 else if (v .le. cb) then 1021 symbol(m) = tb 1022 else if (v .le. cc) then 1023 symbol(m) = tc 1024 else if (v .le. cd) then 1025 symbol(m) = td 1026 else 1027 symbol(m) = te 1028 end if 1029 end if 1030 end do 1031 write (iout,50) (symbol(k),k=1,nrow) 1032 50 format (1x,130a1) 1033 end do 1034 write (iout,60) (dash,i=1,ndash) 1035 60 format (/,1x,130a1) 1036 if (maxj .lt. n) then 1037 write (iout,70) 1038 70 format () 1039 end if 1040 end do 1041 return 1042 end 1043c 1044c 1045c ################################################################ 1046c ## ## 1047c ## subroutine dstmat -- choose values for distance matrix ## 1048c ## ## 1049c ################################################################ 1050c 1051c 1052c "dstmat" selects a distance matrix containing values between 1053c the previously smoothed upper and lower bounds; the distance 1054c values are chosen from uniform distributions, in a triangle 1055c correlated fashion, or using random partial metrization 1056c 1057c 1058 subroutine dstmat (dmx) 1059 use atoms 1060 use disgeo 1061 use inform 1062 use iounit 1063 use keys 1064 implicit none 1065 integer i,j,k,m 1066 integer index,next 1067 integer npart,npair 1068 integer nmetrize 1069 integer mik,mjk,nik,njk 1070 integer, allocatable :: list(:) 1071 real*8 random,fraction 1072 real*8 invbeta,alpha,beta 1073 real*8 corr,mean,stdev 1074 real*8 denom,swap,delta 1075 real*8 percent,eps,gap 1076 real*8 wall,cpu 1077 real*8, allocatable :: value(:) 1078 real*8 dmx(n,*) 1079 logical first,uniform 1080 logical update 1081 character*8 method 1082 character*20 keyword 1083 character*240 record 1084 character*240 string 1085 external random,invbeta 1086 save first,method,update 1087 save npart,percent 1088 save mean,stdev 1089 save alpha,beta 1090 data first / .true. / 1091c 1092c 1093c initialize the method for distance element selection 1094c 1095 if (first) then 1096 first = .false. 1097 method = 'PAIRWISE' 1098 uniform = .false. 1099 update = .true. 1100 npart = 0 1101 percent = 0.0d0 1102 mean = 0.0d0 1103 compact = 0.0d0 1104 beta = 4.0d0 1105c 1106c search each line of the keyword file for options 1107c 1108 do i = 1, nkey 1109 next = 1 1110 record = keyline(i) 1111 call gettext (record,keyword,next) 1112 call upcase (keyword) 1113c 1114c get a distance selection method and extent of metrization 1115c 1116 if (keyword(1:15) .eq. 'TRIAL-DISTANCE ') then 1117 call gettext (record,method,next) 1118 call upcase (method) 1119 if (method .eq. 'HAVEL') then 1120 call getnumb (record,npart,next) 1121 else if (method .eq. 'PARTIAL') then 1122 call getnumb (record,npart,next) 1123 else if (method .eq. 'PAIRWISE') then 1124 string = record(next:240) 1125 read (string,*,err=10,end=10) percent 1126 10 continue 1127 end if 1128c 1129c get a choice of initial mean for the trial distribution 1130c 1131 else if (keyword(1:19) .eq. 'TRIAL-DISTRIBUTION ') then 1132 string = record(next:240) 1133 read (string,*,err=20,end=20) mean 1134 20 continue 1135 update = .false. 1136 end if 1137 end do 1138c 1139c set extent of partial metrization during distance selection 1140c 1141 if (method .eq. 'HAVEL') then 1142 if (npart.le.0 .or. npart.ge.n-1) npart = n 1143 else if (method .eq. 'PARTIAL') then 1144 if (npart.le.0 .or. npart.ge.n-1) npart = 4 1145 else if (method .eq. 'PAIRWISE') then 1146 if (percent.le.0.0d0 .or. percent.ge.100.0d0) 1147 & percent = min(100.0d0,2000.0d0/dble(n)) 1148 end if 1149c 1150c set the initial distribution for selection of trial distances 1151c 1152 if (method .eq. 'CLASSIC') uniform = .true. 1153 if (method .eq. 'TRICOR') uniform = .true. 1154 if (method .eq. 'HAVEL') uniform = .true. 1155 if (uniform) update = .false. 1156 if (update) then 1157c mean = 2.35d0 / log(pathmax) 1158 mean = 1.65d0 / (pathmax)**0.25d0 1159c mean = 1.30d0 / (pathmax)**0.20d0 1160 end if 1161 alpha = beta*mean / (1.0d0-mean) 1162 stdev = sqrt(alpha*beta/(alpha+beta+1.0d0)) / (alpha+beta) 1163 end if 1164c 1165c write out the final choice for distance matrix generation 1166c 1167 if (verbose) then 1168 call settime 1169 if (method .eq. 'CLASSIC') then 1170 write (iout,30) 1171 30 format (/,' Distance Matrix via Uniform Random', 1172 & ' Fractions without Metrization :') 1173 else if (method .eq. 'RANDOM') then 1174 write (iout,40) 1175 40 format (/,' Distance Matrix Generated via Normal', 1176 & ' Fractions without Metrization :') 1177 else if (method .eq. 'TRICOR') then 1178 write (iout,50) 1179 50 format (/,' Distance Matrix Generated via Triangle', 1180 & ' Correlated Fractions :') 1181 else if (method.eq.'HAVEL' .and. npart.lt.n) then 1182 write (iout,60) npart 1183 60 format (/,' Distance Matrix Generated via',i4,'-Atom', 1184 & ' Partial Metrization :') 1185 else if (method .eq. 'HAVEL') then 1186 write (iout,70) 1187 70 format (/,' Distance Matrix Generated via Randomized', 1188 & ' Atom-Based Metrization :') 1189 else if (method .eq. 'PARTIAL') then 1190 write (iout,80) npart 1191 80 format (/,' Distance Matrix Generated via',i4,'-Atom', 1192 & ' Partial Metrization :') 1193 else if (method.eq.'PAIRWISE' .and. percent.lt.100.0d0) then 1194 write (iout,90) percent 1195 90 format (/,' Distance Matrix Generated via',f6.2,'%', 1196 & ' Random Pairwise Metrization :') 1197 else 1198 write (iout,100) 1199 100 format (/,' Distance Matrix Generated via Randomized', 1200 & ' Pairwise Metrization :') 1201 end if 1202 end if 1203c 1204c adjust the distribution for selection of trial distances 1205c 1206 if (uniform) then 1207 write (iout,110) 1208 110 format (/,' Trial Distances Selected at Random from', 1209 & ' Uniform Distribution') 1210 else 1211 if (update) then 1212 alpha = alpha - 0.2d0*sign(sqrt(abs(compact)),compact) 1213 mean = alpha / (alpha+beta) 1214 stdev = sqrt(alpha*beta/(alpha+beta+1.0d0)) / (alpha+beta) 1215 end if 1216 write (iout,120) mean,stdev,alpha,beta 1217 120 format (/,' Trial Distance Beta Distribution :', 1218 & 4x,f5.2,' +/-',f5.2,3x,'Alpha-Beta',2f6.2) 1219 end if 1220c 1221c perform dynamic allocation of some local arrays 1222c 1223 npair = n*(n-1) / 2 1224 allocate (list(npair)) 1225 allocate (value(npair)) 1226c 1227c uniform or Gaussian distributed distances without metrization 1228c 1229 if (method.eq.'CLASSIC' .or. method.eq.'RANDOM') then 1230 do i = 1, n 1231 dmx(i,i) = 0.0d0 1232 end do 1233 do i = 1, n-1 1234 do j = i+1, n 1235 fraction = random () 1236 if (method .eq. 'RANDOM') then 1237 fraction = invbeta (alpha,beta,fraction) 1238 end if 1239 delta = dbnd(i,j) - dbnd(j,i) 1240 dmx(j,i) = dbnd(j,i) + delta*fraction 1241 dmx(i,j) = dmx(j,i) 1242 end do 1243 end do 1244c 1245c Crippen's triangle correlated distance selection 1246c 1247 else if (method .eq. 'TRICOR') then 1248 do i = 1, n 1249 dmx(i,i) = 0.0d0 1250 end do 1251 do i = 1, n-1 1252 do j = i+1, n 1253 dmx(j,i) = random () 1254 dmx(i,j) = dmx(j,i) 1255 end do 1256 end do 1257 do i = 1, n-1 1258 do j = i+1, n 1259 denom = 0.0d0 1260 dmx(i,j) = 0.0d0 1261 do k = 1, n 1262 if (k .ne. i) then 1263 mik = max(i,k) 1264 mjk = max(j,k) 1265 nik = min(i,k) 1266 njk = min(j,k) 1267 if (k .eq. j) then 1268 dmx(i,j) = dmx(i,j) + dmx(j,i) 1269 denom = denom + 1.0d0 1270 else if (dbnd(njk,mjk) .le. 1271 & 0.2d0*dbnd(nik,mik)) then 1272 if (i .gt. k) corr = 0.9d0 * dmx(i,k) 1273 if (k .gt. i) corr = 0.9d0 * dmx(k,i) 1274 dmx(i,j) = dmx(i,j) + corr 1275 denom = denom + 0.9d0 1276 else if (dbnd(nik,mik) .le. 1277 & 0.2d0*dbnd(njk,mjk)) then 1278 if (j .gt. k) corr = 0.9d0 * dmx(j,k) 1279 if (k .gt. j) corr = 0.9d0 * dmx(k,j) 1280 dmx(i,j) = dmx(i,j) + corr 1281 denom = denom + 0.9d0 1282 else if (dbnd(mik,nik) .ge. 1283 & 0.9d0*dbnd(njk,mjk)) then 1284 if (j .gt. k) corr = 0.5d0 * (1.0d0-dmx(j,k)) 1285 if (k .gt. j) corr = 0.5d0 * (1.0d0-dmx(k,j)) 1286 dmx(i,j) = dmx(i,j) + corr 1287 denom = denom + 0.5d0 1288 else if (dbnd(mjk,njk) .ge. 1289 & 0.9d0*dbnd(nik,mik)) then 1290 if (i .gt. k) corr = 0.5d0 * (1.0d0-dmx(i,k)) 1291 if (k .gt. i) corr = 0.5d0 * (1.0d0-dmx(k,i)) 1292 dmx(i,j) = dmx(i,j) + corr 1293 denom = denom + 0.5d0 1294 end if 1295 end if 1296 end do 1297 dmx(i,j) = dmx(i,j) / denom 1298 end do 1299 end do 1300 do i = 1, n-1 1301 do j = i+1, n 1302 delta = dbnd(i,j) - dbnd(j,i) 1303 dmx(i,j) = dbnd(j,i) + delta*dmx(i,j) 1304 dmx(j,i) = dmx(i,j) 1305 end do 1306 end do 1307c 1308c Havel/XPLOR atom-based metrization over various distributions 1309c 1310 else if (method.eq.'HAVEL' .or. method.eq.'PARTIAL') then 1311 do i = 1, n 1312 do j = 1, n 1313 dmx(j,i) = dbnd(j,i) 1314 end do 1315 end do 1316 do i = 1, n 1317 value(i) = random () 1318 end do 1319 call sort2 (n,value,list) 1320 gap = 0.0d0 1321 do i = 1, n-1 1322 k = list(i) 1323 do j = i+1, n 1324 m = list(j) 1325 fraction = random () 1326 if (method .eq. 'PARTIAL') then 1327 fraction = invbeta (alpha,beta,fraction) 1328 end if 1329 delta = abs(dbnd(k,m) - dbnd(m,k)) 1330 if (k .lt. m) then 1331 dbnd(k,m) = dbnd(m,k) + delta*fraction 1332 dbnd(m,k) = dbnd(k,m) 1333 else 1334 dbnd(k,m) = dbnd(k,m) + delta*fraction 1335 dbnd(m,k) = dbnd(k,m) 1336 end if 1337 if (i .le. npart) call trifix (k,m) 1338 if (i .gt. npart) gap = gap + delta 1339 end do 1340 end do 1341 do i = 1, n 1342 do j = 1, n 1343 swap = dmx(j,i) 1344 dmx(j,i) = dbnd(j,i) 1345 dbnd(j,i) = swap 1346 end do 1347 end do 1348 if (verbose .and. npart.lt.n-1) then 1349 write (iout,130) gap/dble((n-npart)*(n-npart-1)/2) 1350 130 format (/,' Average Bound Gap after Partial Metrization :', 1351 & 3x,f12.4) 1352 end if 1353c 1354c use partial randomized pairwise distance-based metrization 1355c 1356 else if (method.eq.'PAIRWISE' .and. percent.le.10.0d0) then 1357 npair = n*(n-1) / 2 1358 nmetrize = nint(0.01d0*percent*dble(npair)) 1359 do i = 1, n 1360 do j = 1, n 1361 dmx(j,i) = dbnd(j,i) 1362 end do 1363 end do 1364 do i = 1, nmetrize 1365 140 continue 1366 k = int(dble(n)*random()) + 1 1367 m = int(dble(n)*random()) + 1 1368 if (dbnd(k,m) .eq. dbnd(m,k)) goto 140 1369 if (k .gt. m) then 1370 j = k 1371 k = m 1372 m = j 1373 end if 1374 fraction = random () 1375 fraction = invbeta (alpha,beta,fraction) 1376 delta = dbnd(k,m) - dbnd(m,k) 1377 dbnd(k,m) = dbnd(m,k) + delta*fraction 1378 dbnd(m,k) = dbnd(k,m) 1379 call trifix (k,m) 1380 end do 1381 gap = 0.0d0 1382 do i = 1, n-1 1383 do j = i, n 1384 delta = dbnd(i,j) - dbnd(j,i) 1385 if (delta .ne. 0.0d0) then 1386 gap = gap + delta 1387 fraction = random () 1388 fraction = invbeta (alpha,beta,fraction) 1389 dbnd(i,j) = dbnd(j,i) + delta*fraction 1390 dbnd(j,i) = dbnd(i,j) 1391 end if 1392 end do 1393 end do 1394 do i = 1, n 1395 do j = 1, n 1396 swap = dmx(j,i) 1397 dmx(j,i) = dbnd(j,i) 1398 dbnd(j,i) = swap 1399 end do 1400 end do 1401 if (verbose .and. nmetrize.lt.npair) then 1402 write (iout,150) gap/dble(npair-nmetrize) 1403 150 format (/,' Average Bound Gap after Partial Metrization :', 1404 & 3x,f12.4) 1405 end if 1406c 1407c use randomized pairwise distance-based metrization 1408c 1409 else if (method .eq. 'PAIRWISE') then 1410 npair = n*(n-1) / 2 1411 nmetrize = nint(0.01d0*percent*dble(npair)) 1412 do i = 1, n 1413 do j = 1, n 1414 dmx(j,i) = dbnd(j,i) 1415 end do 1416 end do 1417 do i = 1, npair 1418 value(i) = random () 1419 end do 1420 call sort2 (npair,value,list) 1421 eps = 1.0d-10 1422 gap = 0.0d0 1423 do i = 1, npair 1424 index = list(i) 1425 k = int(0.5d0 * (dble(2*n+1) 1426 & - sqrt(dble(4*n*(n-1)-8*index+9))) + eps) 1427 m = n*(1-k) + k*(k+1)/2 + index 1428 fraction = random () 1429 fraction = invbeta (alpha,beta,fraction) 1430 delta = dbnd(k,m) - dbnd(m,k) 1431 dbnd(k,m) = dbnd(m,k) + delta*fraction 1432 dbnd(m,k) = dbnd(k,m) 1433 if (i .le. nmetrize) call trifix (k,m) 1434 if (i .gt. nmetrize) gap = gap + delta 1435 end do 1436 do i = 1, n 1437 do j = 1, n 1438 swap = dmx(j,i) 1439 dmx(j,i) = dbnd(j,i) 1440 dbnd(j,i) = swap 1441 end do 1442 end do 1443 if (verbose .and. nmetrize.lt.npair) then 1444 write (iout,160) gap/dble(npair-nmetrize) 1445 160 format (/,' Average Bound Gap after Partial Metrization :', 1446 & 3x,f12.4) 1447 end if 1448 end if 1449c 1450c perform deallocation of some local arrays 1451c 1452 deallocate (list) 1453 deallocate (value) 1454c 1455c get the time required for distance matrix generation 1456c 1457 if (verbose) then 1458 call gettime (wall,cpu) 1459 write (iout,170) wall 1460 170 format (/,' Time Required for Distance Matrix :',5x,f12.2, 1461 & ' seconds') 1462 end if 1463 return 1464 end 1465c 1466c 1467c ############################################################### 1468c ## ## 1469c ## subroutine metric -- computation of the metric matrix ## 1470c ## ## 1471c ############################################################### 1472c 1473c 1474c "metric" takes as input the trial distance matrix and computes 1475c the metric matrix of all possible dot products between the atomic 1476c vectors and the center of mass using the law of cosines and the 1477c following formula for the distances to the center of mass: 1478c 1479c dcm(i)**2 = (1/n) * sum(j=1,n)(dist(i,j)**2) 1480c - (1/n**2) * sum(j<k)(dist(j,k)**2) 1481c 1482c upon output, the metric matrix is stored in the lower triangle 1483c plus diagonal of the input trial distance matrix, the upper 1484c triangle of the input matrix is unchanged 1485c 1486c literature reference: 1487c 1488c G. M. Crippen and T. F. Havel, "Stable Calculation of Coordinates 1489c from Distance Information", Acta Cryst., A34, 282-284 (1978) 1490c 1491c 1492 subroutine metric (gmx,nneg) 1493 use atoms 1494 use disgeo 1495 use inform 1496 use iounit 1497 implicit none 1498 integer i,j,nneg 1499 real*8 total,sum,rg 1500 real*8 gmx(n,*) 1501 real*8, allocatable :: dsq(:) 1502 real*8, allocatable :: dcm(:) 1503c 1504c 1505c square and sum trial distances to get radius of gyration 1506c 1507 total = 0.0d0 1508 do i = 1, n 1509 do j = i, n 1510 gmx(j,i) = gmx(j,i)**2 1511 total = total + gmx(j,i) 1512 end do 1513 end do 1514 total = total / dble(n**2) 1515 rg = sqrt(total) 1516 if (verbose) then 1517 write (iout,10) rg 1518 10 format (/,' Radius of Gyration before Embedding :',7x,f16.4) 1519 end if 1520c 1521c perform dynamic allocation of some local arrays 1522c 1523 allocate (dsq(n)) 1524 allocate (dcm(n)) 1525c 1526c sum squared distances from each atom; the center 1527c of mass is derived using the formula shown above 1528c 1529 nneg = 0 1530 do i = 1, n 1531 sum = 0.0d0 1532 do j = 1, i-1 1533 sum = sum + gmx(i,j) 1534 end do 1535 do j = i, n 1536 sum = sum + gmx(j,i) 1537 end do 1538 dsq(i) = sum/dble(n) - total 1539 dcm(i) = sqrt(abs(dsq(i))) 1540 if (dsq(i) .lt. 0.0d0) then 1541 nneg = nneg + 1 1542 dcm(i) = -dcm(i) 1543 end if 1544 end do 1545 if (verbose .and. n.le.130) then 1546 write (iout,20) 1547 20 format (/,' Atomic Distances to the Center of Mass :',/) 1548 write (iout,30) (dcm(i),i=1,n) 1549 30 format (6f13.4) 1550 end if 1551c 1552c calculate the metric matrix using the law of cosines, and 1553c place into the lower triangle of the input distance matrix 1554c 1555 do i = 1, n 1556 do j = i, n 1557 gmx(j,i) = 0.5d0 * (dsq(i)+dsq(j)-gmx(j,i)) 1558 end do 1559 end do 1560c 1561c perform deallocation of some local arrays 1562c 1563 deallocate (dsq) 1564 deallocate (dcm) 1565 return 1566 end 1567c 1568c 1569c ################################################################## 1570c ## ## 1571c ## subroutine eigen -- largest eigenvalues of metric metrix ## 1572c ## ## 1573c ################################################################## 1574c 1575c 1576c "eigen" uses the power method to compute the largest eigenvalues 1577c and eigenvectors of the metric matrix, "valid" is set true if the 1578c first three eigenvalues are positive 1579c 1580c 1581 subroutine eigen (evl,evc,gmx,valid) 1582 use atoms 1583 use inform 1584 use iounit 1585 implicit none 1586 integer i,j,neigen 1587 real*8 wall,cpu 1588 real*8 evl(*) 1589 real*8 evc(n,*) 1590 real*8 gmx(n,*) 1591 logical valid 1592c 1593c 1594c initialize number of eigenvalues and convergence criteria 1595c 1596 if (verbose) call settime 1597 neigen = 3 1598c 1599c compute largest eigenvalues via power method with deflation 1600c 1601 call deflate (n,neigen,gmx,evl,evc) 1602c 1603c check to see if the first three eigenvalues are positive 1604c 1605 valid = .true. 1606 do i = 1, 3 1607 if (evl(i) .lt. 0.0d0) valid = .false. 1608 end do 1609c 1610c print out the eigenvalues and their eigenvectors 1611c 1612 if (verbose) then 1613 write (iout,10) 1614 10 format (/,' Eigenvalues from Metric Matrix :',/) 1615 write (iout,20) (evl(i),i=1,neigen) 1616 20 format (5f15.4) 1617 end if 1618 if (debug) then 1619 write (iout,30) 1620 30 format (/,' Eigenvectors from Metric Matrix :',/) 1621 do i = 1, n 1622 write (iout,40) (evc(i,j),j=1,neigen) 1623 40 format (5f15.4) 1624 end do 1625 end if 1626c 1627c get the time required for partial matrix diagonalization 1628c 1629 if (verbose) then 1630 call gettime (wall,cpu) 1631 write (iout,50) wall 1632 50 format (/,' Time Required for Eigenvalues :',9x,f12.2, 1633 & ' seconds') 1634 end if 1635 return 1636 end 1637c 1638c 1639c ################################################################## 1640c ## ## 1641c ## subroutine coords -- converts eigenvalues to coordinates ## 1642c ## ## 1643c ################################################################## 1644c 1645c 1646c "coords" converts the three principal eigenvalues/vectors from 1647c the metric matrix into atomic coordinates, and calls a routine 1648c to compute the rms deviation from the bounds 1649c 1650c 1651 subroutine coords (evl,evc) 1652 use atoms 1653 use disgeo 1654 use inform 1655 use iounit 1656 implicit none 1657 integer i,j,neigen 1658 real*8 rg 1659 real*8 evl(*) 1660 real*8 evc(n,*) 1661 character*240 title 1662c 1663c 1664c compute coordinates from the largest eigenvalues and vectors 1665c 1666 neigen = 3 1667 do j = 1, neigen 1668 evl(j) = sqrt(abs(evl(j))) 1669 end do 1670 do j = 1, neigen 1671 do i = 1, n 1672 evc(i,j) = evl(j) * evc(i,j) 1673 end do 1674 end do 1675c 1676c transfer the final coordinates back to atomic vectors 1677c 1678 do i = 1, n 1679 x(i) = evc(i,1) 1680 y(i) = evc(i,2) 1681 z(i) = evc(i,3) 1682 end do 1683c 1684c find the rms bounds deviations and radius of gyration 1685c 1686 if (verbose) then 1687 title = 'after Embedding :' 1688 call rmserror (title) 1689 call gyrate (rg) 1690 write (iout,10) rg 1691 10 format (/,' Radius of Gyration after Embedding :',8x,f16.4) 1692 end if 1693 return 1694 end 1695c 1696c 1697c ################################################################ 1698c ## ## 1699c ## subroutine chksize -- estimate compaction or expansion ## 1700c ## ## 1701c ################################################################ 1702c 1703c 1704c "chksize" computes a measure of overall global structural 1705c expansion or compaction from the number of excess upper 1706c or lower bounds matrix violations 1707c 1708c 1709 subroutine chksize 1710 use atoms 1711 use couple 1712 use disgeo 1713 use inform 1714 use iounit 1715 implicit none 1716 integer i,k,npair,nskip 1717 integer nlarge,nsmall 1718 integer, allocatable :: skip(:) 1719 real*8 xi,yi,zi 1720 real*8 dstsq,bupsq,blosq 1721c 1722c 1723c perform dynamic allocation of some local arrays 1724c 1725 allocate (skip(n)) 1726c 1727c zero out the list of atoms locally connected to each atom 1728c 1729 nskip = 0 1730 do i = 1, n 1731 skip(i) = 0 1732 end do 1733c 1734c initialize counters, total pair number, and cutoff distance 1735c 1736 nlarge = 0 1737 nsmall = 0 1738 npair = n*(n-1) / 2 1739c 1740c count the number of excess upper or lower bound violations 1741c 1742 do i = 1, n-1 1743 xi = x(i) 1744 yi = y(i) 1745 zi = z(i) 1746c do k = 1, n12(i) 1747c skip(i12(k,i)) = i 1748c end do 1749c do k = 1, n13(i) 1750c skip(i13(k,i)) = i 1751c end do 1752c do k = 1, n14(i) 1753c skip(i14(k,i)) = i 1754c end do 1755 do k = i+1, n 1756 if (skip(k) .eq. i) then 1757 nskip = nskip + 1 1758 else 1759 dstsq = (x(k)-xi)**2 + (y(k)-yi)**2 + (z(k)-zi)**2 1760 bupsq = dbnd(i,k)**2 1761 blosq = dbnd(k,i)**2 1762 if (dstsq .gt. bupsq) then 1763 nlarge = nlarge + 1 1764 else if (blosq .gt. dstsq) then 1765 nsmall = nsmall + 1 1766 end if 1767 end if 1768 end do 1769 end do 1770c 1771c perform deallocation of some local arrays 1772c 1773 deallocate (skip) 1774c 1775c set the value for the overall index of compaction 1776c 1777 compact = 100.0d0 * dble(nlarge-nsmall)/dble(npair-nskip) 1778 if (verbose) then 1779 write (iout,10) compact 1780 10 format (/,' Index of Structure Expansion/Compaction :', 1781 & 7x,f12.4) 1782 end if 1783 return 1784 end 1785c 1786c 1787c ############################################################### 1788c ## ## 1789c ## subroutine majorize -- Guttman transform majorization ## 1790c ## ## 1791c ############################################################### 1792c 1793c 1794c "majorize" refines the projected coordinates by attempting to 1795c minimize the least square residual between the trial distance 1796c matrix and the distances computed from the coordinates 1797c 1798c literature reference: 1799c 1800c T. F. Havel, "An Evaluation of Computational Strategies for 1801c Use in the Determination of Protein Structure from Distance 1802c Constraints obtained by Nuclear Magnetic Resonance", Progress 1803c in Biophysics and Molecular Biology, 56, 43-78 (1991) 1804c 1805c 1806 subroutine majorize (dmx) 1807 use atoms 1808 use inform 1809 use iounit 1810 implicit none 1811 integer i,k,iter 1812 integer niter,period 1813 real*8 pairs,rg 1814 real*8 dn1,dn2 1815 real*8 wall,cpu 1816 real*8 target,dist,error 1817 real*8 rmserr,average 1818 real*8 xi,yi,zi 1819 real*8, allocatable :: b(:) 1820 real*8, allocatable :: xx(:) 1821 real*8, allocatable :: yy(:) 1822 real*8, allocatable :: zz(:) 1823 real*8 dmx(n,*) 1824 character*240 title 1825c 1826c 1827c set number of iterations and some other needed values 1828c 1829 if (verbose) call settime 1830 niter = 20 1831 period = 5 1832 pairs = dble(n*(n-1)/2) 1833 dn1 = dble(n-1) 1834 dn2 = dble(n*n) 1835c 1836c find the average and rms error from trial distances 1837c 1838 iter = 0 1839 rmserr = 0.0d0 1840 average = 0.0d0 1841 do i = 1, n-1 1842 xi = x(i) 1843 yi = y(i) 1844 zi = z(i) 1845 do k = i+1, n 1846 target = dmx(k,i) 1847 dist = sqrt((x(k)-xi)**2+(y(k)-yi)**2+(z(k)-zi)**2) 1848 error = dist - target 1849 rmserr = rmserr + error**2 1850 average = average + error/target 1851 end do 1852 end do 1853 rmserr = sqrt(rmserr/pairs) 1854 average = 100.0d0 * average / pairs 1855c 1856c write a header with the initial error values 1857c 1858 if (verbose) then 1859 write (iout,10) 1860 10 format (/,' Majorization to Trial Distances using', 1861 & ' Constant Weights :', 1862 & //,4x,'Iteration',6x,'RMS Error',5x,'Ave % Error',/) 1863 write (iout,20) iter,rmserr,average 1864 20 format (5x,i5,2f16.4) 1865 end if 1866c 1867c perform dynamic allocation of some local arrays 1868c 1869 allocate (b(n)) 1870 allocate (xx(n)) 1871 allocate (yy(n)) 1872 allocate (zz(n)) 1873c 1874c initialize the transformed coordinates for each atom 1875c 1876 do iter = 1, niter 1877 do i = 1, n 1878 xi = x(i) 1879 yi = y(i) 1880 zi = z(i) 1881 b(i) = 0.0d0 1882 xx(i) = 0.0d0 1883 yy(i) = 0.0d0 1884 zz(i) = 0.0d0 1885c 1886c form a single row of the B matrix assuming unity weights 1887c 1888 do k = 1, i-1 1889 dist = sqrt((x(k)-xi)**2+(y(k)-yi)**2+(z(k)-zi)**2) 1890 b(k) = -dmx(k,i) / dist 1891 b(i) = b(i) - b(k) 1892 end do 1893 do k = i+1, n 1894 dist = sqrt((x(k)-xi)**2+(y(k)-yi)**2+(z(k)-zi)**2) 1895 b(k) = -dmx(k,i) / dist 1896 b(i) = b(i) - b(k) 1897 end do 1898c 1899c multiply the row of the B matrix by the atomic coordinates 1900c 1901 do k = 1, n 1902 xx(i) = xx(i) + b(k)*x(k) 1903 yy(i) = yy(i) + b(k)*y(k) 1904 zz(i) = zz(i) + b(k)*z(k) 1905 end do 1906 end do 1907c 1908c move the intermediate values into the coordinate arrays 1909c 1910 do i = 1, n 1911 x(i) = xx(i) 1912 y(i) = yy(i) 1913 z(i) = zz(i) 1914 end do 1915c 1916c multiply the inverse weight matrix S+ by the coordinates 1917c 1918 do i = 1, n 1919 xx(i) = (dn1/dn2) * x(i) 1920 yy(i) = (dn1/dn2) * y(i) 1921 zz(i) = (dn1/dn2) * z(i) 1922 do k = 1, i-1 1923 xx(i) = xx(i) - x(k)/dn2 1924 yy(i) = yy(i) - y(k)/dn2 1925 zz(i) = zz(i) - z(k)/dn2 1926 end do 1927 do k = i+1, n 1928 xx(i) = xx(i) - x(k)/dn2 1929 yy(i) = yy(i) - y(k)/dn2 1930 zz(i) = zz(i) - z(k)/dn2 1931 end do 1932 end do 1933c 1934c copy the new coordinates into their permanent arrays 1935c 1936 do i = 1, n 1937 x(i) = xx(i) 1938 y(i) = yy(i) 1939 z(i) = zz(i) 1940 end do 1941c 1942c find the average and rms error from trial distances 1943c 1944 rmserr = 0.0d0 1945 average = 0.0d0 1946 do i = 1, n-1 1947 xi = x(i) 1948 yi = y(i) 1949 zi = z(i) 1950 do k = i+1, n 1951 target = dmx(k,i) 1952 dist = sqrt((x(k)-xi)**2+(y(k)-yi)**2+(z(k)-zi)**2) 1953 error = dist - target 1954 rmserr = rmserr + error**2 1955 average = average + error/target 1956 end do 1957 end do 1958 rmserr = sqrt(rmserr/pairs) 1959 average = 100.0d0 * average / pairs 1960 if (verbose .and. mod(iter,period).eq.0) then 1961 write (iout,30) iter,rmserr,average 1962 30 format (5x,i5,2f16.4) 1963 end if 1964 end do 1965c 1966c perform deallocation of some local arrays 1967c 1968 deallocate (b) 1969 deallocate (xx) 1970 deallocate (yy) 1971 deallocate (zz) 1972c 1973c find the rms bounds deviations and radius of gyration 1974c 1975 if (verbose) then 1976 title = 'after Majorization :' 1977 call rmserror (title) 1978 call gyrate (rg) 1979 write (iout,40) rg 1980 40 format (/,' Radius of Gyration after Majorization :',5x,f16.4) 1981c 1982c get the time required for the majorization procedure 1983c 1984 call gettime (wall,cpu) 1985 write (iout,50) wall 1986 50 format (/,' Time Required for Majorization :',8x,f12.2, 1987 & ' seconds') 1988 end if 1989 return 1990 end 1991c 1992c 1993c ################################################################ 1994c ## ## 1995c ## subroutine refine -- minimization of initial embedding ## 1996c ## ## 1997c ################################################################ 1998c 1999c 2000c "refine" performs minimization of the atomic coordinates 2001c of an initial crude embedded distance geometry structure versus 2002c the bound, chirality, planarity and torsional error functions 2003c 2004c 2005 subroutine refine (mode,fctval,grdmin) 2006 use atoms 2007 use disgeo 2008 use inform 2009 use iounit 2010 use minima 2011 use output 2012 implicit none 2013 integer i,nvar 2014 real*8 initerr,miderr,toterr 2015 real*8 fctval,grdmin 2016 real*8, allocatable :: xx(:) 2017 character*7 mode 2018 external initerr,miderr 2019 external toterr,optsave 2020c 2021c 2022c perform dynamic allocation of some local arrays 2023c 2024 nvar = 3 * n 2025 allocate (xx(nvar)) 2026c 2027c convert atomic coordinates to optimization parameters 2028c 2029 nvar = 0 2030 do i = 1, n 2031 nvar = nvar + 1 2032 xx(nvar) = x(i) 2033 nvar = nvar + 1 2034 xx(nvar) = y(i) 2035 nvar = nvar + 1 2036 xx(nvar) = z(i) 2037 end do 2038c 2039c set values of parameters needed for optimization 2040c 2041 coordtype = 'NONE' 2042 cyclesave = .true. 2043c grdmin = 0.01d0 2044 maxiter = 2 * nvar 2045 iwrite = 0 2046c iprint = 0 2047c if (verbose) iprint = 10 2048c 2049c minimize initially only on the local geometry and torsions, 2050c then on local geometry and chirality, torsions, and finally 2051c minimize on all distance bounds, chirality and torsions 2052c 2053 if (mode .eq. 'INITIAL') then 2054 call lbfgs (nvar,xx,fctval,grdmin,initerr,optsave) 2055 else if (mode .eq. 'MIDDLE') then 2056 call lbfgs (nvar,xx,fctval,grdmin,miderr,optsave) 2057 else if (mode .eq. 'FINAL') then 2058 call lbfgs (nvar,xx,fctval,grdmin,toterr,optsave) 2059 end if 2060c 2061c convert optimization parameters to atomic coordinates 2062c 2063 nvar = 0 2064 do i = 1, n 2065 nvar = nvar + 1 2066 x(i) = xx(nvar) 2067 nvar = nvar + 1 2068 y(i) = xx(nvar) 2069 nvar = nvar + 1 2070 z(i) = xx(nvar) 2071 end do 2072c 2073c perform deallocation of some local arrays 2074c 2075 deallocate (xx) 2076 return 2077 end 2078c 2079c 2080c ############################################################## 2081c ## ## 2082c ## subroutine explore -- simulated annealing refinement ## 2083c ## ## 2084c ############################################################## 2085c 2086c 2087c "explore" uses simulated annealing on an initial crude 2088c embedded distance geoemtry structure to refine versus the 2089c bound, chirality, planarity and torsional error functions 2090c 2091c 2092 subroutine explore (mode,nstep,dt,mass,temp_start,temp_stop,v,a) 2093 use atoms 2094 use inform 2095 use iounit 2096 use math 2097 use units 2098 implicit none 2099 integer i,istep,nstep 2100 integer nvar,period 2101 real*8 error,total 2102 real*8 prior,change 2103 real*8 dt,dt2,dt_2,dt2_2 2104 real*8 xbig,xrms,mass,kinetic 2105 real*8 temp_start,temp_stop 2106 real*8 tau_start,tau_stop 2107 real*8 ratio,sigmoid,scale 2108 real*8 target,temp,tautemp 2109 real*8 initerr,miderr,toterr 2110 real*8 v(*) 2111 real*8 a(*) 2112 real*8, allocatable :: xx(:) 2113 real*8, allocatable :: xmove(:) 2114 real*8, allocatable :: g(:) 2115 character*7 mode 2116c 2117c 2118c set values of the basic simulated annealing parameters 2119c 2120c nstep = 5000 2121c dt = 0.1d0 2122c temp_start = 200.0d0 2123c temp_stop = 0.0d0 2124c mass = 1000.0d0 2125c 2126c perform dynamic allocation of some local arrays 2127c 2128 nvar = 3 * n 2129 allocate (xx(nvar)) 2130 allocate (xmove(nvar)) 2131 allocate (g(nvar)) 2132c 2133c convert atomic coordinates to annealing variables 2134c 2135 nvar = 0 2136 do i = 1, n 2137 nvar = nvar + 1 2138 xx(nvar) = x(i) 2139 nvar = nvar + 1 2140 xx(nvar) = y(i) 2141 nvar = nvar + 1 2142 xx(nvar) = z(i) 2143 end do 2144c 2145c initialize the velocities, accelerations and other parameters 2146c 2147 dt2 = dt * dt 2148 dt_2 = dt / 2.0d0 2149 dt2_2 = dt2 / 2.0d0 2150 period = 100 2151 tau_start = 100.0d0 * dt 2152 tau_stop = 10.0d0 * dt 2153 tautemp = tau_start 2154c 2155c print a header for the simulated annealing protocol 2156c 2157 write (iout,10) 2158 10 format (/,' Molecular Dynamics Simulated Annealing Refinement :') 2159 write (iout,20) nstep,dt,log(mass)/logten,temp_start,temp_stop 2160 20 format (/,' Steps:',i6,3x,'Time/Step:',f6.3,' ps',3x, 2161 & 'LogMass:',f5.2,3x,'Temp:',f6.1,' to',f6.1) 2162c 2163c get the total error and temperature at start of dynamics 2164c 2165 if (mode .eq. 'INITIAL') then 2166 error = initerr (xx,g) 2167 else if (mode .eq. 'MIDDLE') then 2168 error = miderr (xx,g) 2169 else if (mode .eq. 'FINAL') then 2170 error = toterr (xx,g) 2171 end if 2172 kinetic = 0.0d0 2173 do i = 1, nvar 2174 kinetic = kinetic + mass*v(i)**2 2175 end do 2176 kinetic = 0.5d0 * kinetic / ekcal 2177 temp = 2.0d0 * kinetic / (dble(nvar) * gasconst) 2178 total = error + kinetic 2179 prior = total 2180 if (verbose) then 2181 write (iout,30) 2182 30 format (/,' MD Step E Total E Potential E Kinetic', 2183 & ' Temp MaxMove RMS Move',/) 2184 write (iout,40) 0,total,error,kinetic,temp 2185 40 format (i6,2f13.4,f12.4,f11.2) 2186 end if 2187c 2188c find new positions and half-step velocities via Verlet 2189c 2190 do istep = 1, nstep 2191 xbig = 0.0d0 2192 xrms = 0.0d0 2193 do i = 1, nvar 2194 xmove(i) = v(i)*dt + a(i)*dt2_2 2195 xx(i) = xx(i) + xmove(i) 2196 v(i) = v(i) + a(i)*dt_2 2197 if (abs(xmove(i)) .gt. xbig) xbig = abs(xmove(i)) 2198 xrms = xrms + xmove(i)**2 2199 end do 2200 xrms = sqrt(xrms/dble(nvar)) 2201c 2202c get the error function value and gradient 2203c 2204 if (mode .eq. 'INITIAL') then 2205 error = initerr (xx,g) 2206 else if (mode .eq. 'MIDDLE') then 2207 error = miderr (xx,g) 2208 else if (mode .eq. 'FINAL') then 2209 error = toterr (xx,g) 2210 end if 2211c 2212c use Newton's second law to get the next accelerations; 2213c find the full-step velocities using the Verlet recursion 2214c 2215 do i = 1, nvar 2216 a(i) = -ekcal * g(i) / mass 2217 v(i) = v(i) + a(i)*dt_2 2218 end do 2219c 2220c find the total kinetic energy and system temperature 2221c 2222 kinetic = 0.0d0 2223 do i = 1, nvar 2224 kinetic = kinetic + mass*v(i)**2 2225 end do 2226 kinetic = 0.5d0 * kinetic / ekcal 2227 temp = 2.0d0 * kinetic / (dble(nvar) * gasconst) 2228 if (temp .eq. 0.0d0) temp = 0.1d0 2229c 2230c set target temperature and coupling via a sigmoidal cooling 2231c 2232 ratio = dble(istep) / dble(nstep) 2233 ratio = sigmoid (3.5d0,ratio) 2234 target = temp_start*(1.0d0-ratio) + temp_stop*ratio 2235 tautemp = tau_start*(1.0d0-ratio) + tau_stop*ratio 2236c 2237c couple to external temperature bath via velocity scaling 2238c 2239 scale = sqrt(1.0d0 + (dt/tautemp)*(target/temp-1.0d0)) 2240 do i = 1, nvar 2241 v(i) = scale * v(i) 2242 end do 2243c 2244c write results for the current annealing step 2245c 2246 total = error + kinetic 2247 if (verbose .and. mod(istep,period).eq.0) then 2248 write (iout,50) istep,total,error,kinetic,temp,xbig,xrms 2249 50 format (i6,2f13.4,f12.4,f11.2,2f10.4) 2250 end if 2251c 2252c check the energy change for instability in the dynamics 2253c 2254 change = total - prior 2255 if (change .gt. dble(n)) then 2256 do i = 1, nvar 2257 xx(i) = xx(i) - xmove(i) 2258 end do 2259 if (verbose .and. mod(istep,period).ne.0) then 2260 write (iout,60) istep,total,error,kinetic,temp,xbig,xrms 2261 60 format (i6,2f13.4,f12.4,f11.2,2f10.4) 2262 end if 2263 write (iout,70) 2264 70 format (/,' EXPLORE -- Simulated Annealing Unstable;', 2265 & ' Switching to Minimization') 2266 goto 80 2267 end if 2268 end do 2269c 2270c convert annealing variables to atomic coordinates 2271c 2272 80 continue 2273 nvar = 0 2274 do i = 1, n 2275 nvar = nvar + 1 2276 x(i) = xx(nvar) 2277 nvar = nvar + 1 2278 y(i) = xx(nvar) 2279 nvar = nvar + 1 2280 z(i) = xx(nvar) 2281 end do 2282c 2283c perform deallocation of some local arrays 2284c 2285 deallocate (xx) 2286 deallocate (xmove) 2287 deallocate (g) 2288 return 2289 end 2290c 2291c 2292c ################################################################# 2293c ## ## 2294c ## subroutine fracdist -- fractional distance distribution ## 2295c ## ## 2296c ################################################################# 2297c 2298c 2299c "fracdist" computes a normalized distribution of the pairwise 2300c fractional distances between the smoothed upper and lower bounds 2301c 2302c literature reference: 2303c 2304c C. M. Oshiro, J. Thomason and I. D. Kuntz, "Effects of Limited 2305c Input Distance Constraints Upon the Distance Geometry Algorithm", 2306c Biopolymers, 31, 1049-1064 (1991) 2307c 2308c 2309 subroutine fracdist (title) 2310 use atoms 2311 use disgeo 2312 use iounit 2313 implicit none 2314 integer start,stop 2315 parameter (start=-20) 2316 parameter (stop=120) 2317 integer i,j,k,sum 2318 integer trimtext 2319 integer bin(start:stop) 2320 integer bin2(start:stop) 2321 real*8 xi,yi,zi,size 2322 real*8 dist,range,fraction 2323 real*8 fdist(start:stop) 2324 real*8 fdist2(start:stop) 2325 character*240 title 2326c 2327c 2328c set the bin size and zero out the individual bins 2329c 2330 size = 0.01d0 2331 do i = start, stop 2332 bin(i) = 0 2333 bin2(i) = 0 2334 end do 2335c 2336c get distribution of fractional distances between bounds 2337c 2338 do i = 1, n-1 2339 xi = x(i) 2340 yi = y(i) 2341 zi = z(i) 2342 do j = i+1, n 2343 dist = sqrt((x(j)-xi)**2 + (y(j)-yi)**2 + (z(j)-zi)**2) 2344 range = dbnd(i,j) - dbnd(j,i) 2345 if (range .ge. 1.0d0) then 2346 fraction = (dist-dbnd(j,i)) / range 2347 k = nint(fraction / size) 2348 k = max(start,min(stop,k)) 2349 bin(k) = bin(k) + 1 2350 if (range .ge. 0.8d0*pathmax) bin2(k) = bin2(k) + 1 2351 end if 2352 end do 2353 end do 2354c 2355c normalize the fractional distance frequency distribution 2356c 2357 sum = 0 2358 do i = start, stop 2359 sum = sum + bin(i) 2360 end do 2361 do i = start, stop 2362 fdist(i) = dble(bin(i)) / (size*dble(sum)) 2363 end do 2364 sum = 0 2365 do i = start, stop 2366 sum = sum + bin2(i) 2367 end do 2368 do i = start, stop 2369 fdist2(i) = dble(bin2(i)) / (size*dble(sum)) 2370 end do 2371c 2372c print the normalized fractional distance probability 2373c 2374 write (iout,10) title(1:trimtext(title)) 2375 10 format (/,' Fractional Distance Distribution ',a,/) 2376 do i = start, stop 2377 write (iout,20) size*dble(i),fdist(i),fdist2(i) 2378 20 format (8x,f8.4,8x,f8.4,8x,f8.4) 2379 end do 2380 return 2381 end 2382c 2383c 2384c ############################################################## 2385c ## ## 2386c ## subroutine rmserror -- rms bound and restraint error ## 2387c ## ## 2388c ############################################################## 2389c 2390c 2391c "rmserror" computes the maximum absolute deviation and the 2392c rms deviation from the distance bounds, and the number and 2393c rms value of the distance restraint violations 2394c 2395c 2396 subroutine rmserror (title) 2397 use atoms 2398 use disgeo 2399 use iounit 2400 use restrn 2401 implicit none 2402 integer i,j,k,npair 2403 integer nhierr,nloerr 2404 integer ihi,jhi,ilo,jlo 2405 integer trimtext 2406 real*8 rms,himax,lomax 2407 real*8 dist,hierr,loerr 2408 character*240 title 2409c 2410c 2411c search all atom pairs for maximal bounds deviations 2412c 2413 npair = n*(n-1) / 2 2414 nloerr = 0 2415 nhierr = 0 2416 ilo = 0 2417 jlo = 0 2418 ihi = 0 2419 jhi = 0 2420 rms = 0.0d0 2421 lomax = 0.0d0 2422 himax = 0.0d0 2423 do i = 1, n-1 2424 do j = i+1, n 2425 dist = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2 2426 dist = sqrt(dist) 2427 hierr = dist - dbnd(i,j) 2428 if (hierr .gt. 0.0d0) then 2429 nhierr = nhierr + 1 2430 rms = rms + hierr**2 2431 if (hierr .gt. himax) then 2432 himax = hierr 2433 ihi = i 2434 jhi = j 2435 end if 2436 end if 2437 loerr = dbnd(j,i) - dist 2438 if (loerr .gt. 0.0d0) then 2439 nloerr = nloerr + 1 2440 rms = rms + loerr**2 2441 if (loerr .gt. lomax) then 2442 lomax = loerr 2443 ilo = i 2444 jlo = j 2445 end if 2446 end if 2447 end do 2448 end do 2449 rms = sqrt(rms/dble(n*(n-1)/2)) 2450c 2451c print the maximal and rms bound deviations 2452c 2453 write (iout,10) title(1:trimtext(title)) 2454 10 format (/,' Fit to Bounds ',a) 2455 write (iout,20) nhierr,npair,nloerr,npair,himax, 2456 & ihi,jhi,lomax,ilo,jlo,rms 2457 20 format (/,' Num Upper Bound Violations :',4x,i11,' of ',i12, 2458 & /,' Num Lower Bound Violations :',4x,i11,' of ',i12, 2459 & /,' Max Upper Bound Violation :',4x,f12.4,' at ',2i6, 2460 & /,' Max Lower Bound Violation :',4x,f12.4,' at ',2i6, 2461 & /,' RMS Deviation from Bounds :',4x,f12.4) 2462c 2463c search the list of distance restraints for violations 2464c 2465 if (ndfix .gt. 0) then 2466 nloerr = 0 2467 nhierr = 0 2468 ilo = 0 2469 jlo = 0 2470 ihi = 0 2471 jhi = 0 2472 rms = 0.0d0 2473 himax = 0.0d0 2474 lomax = 0.0d0 2475 do k = 1, ndfix 2476 i = idfix(1,k) 2477 j = idfix(2,k) 2478 dist = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2 2479 dist = sqrt(dist) 2480 if (dist .lt. dfix(2,k)) then 2481 nloerr = nloerr + 1 2482 loerr = dfix(2,k) - dist 2483 rms = rms + loerr**2 2484 if (loerr .gt. lomax) then 2485 lomax = loerr 2486 ilo = i 2487 jlo = j 2488 end if 2489 else if (dist .gt. dfix(3,k)) then 2490 nhierr = nhierr + 1 2491 hierr = dist - dfix(3,k) 2492 rms = rms + hierr**2 2493 if (hierr .gt. himax) then 2494 himax = hierr 2495 ihi = i 2496 jhi = j 2497 end if 2498 end if 2499 end do 2500 rms = sqrt(rms/dble(ndfix)) 2501c 2502c print total number and rms value of restraint violations 2503c 2504 write (iout,30) nhierr,ndfix,nloerr,ndfix,himax, 2505 & ihi,jhi,lomax,ilo,jlo,rms 2506 30 format (/,' Num Upper Restraint Violations :',i11,' of ',i12, 2507 & /,' Num Lower Restraint Violations :',i11,' of ',i12, 2508 & /,' Max Upper Restraint Violation :',f12.4,' at ',2i6, 2509 & /,' Max Lower Restraint Violation :',f12.4,' at ',2i6, 2510 & /,' RMS Restraint Dist Violation : ',f12.4) 2511 end if 2512 return 2513 end 2514c 2515c 2516c ############################################################## 2517c ## ## 2518c ## subroutine dmdump -- final distance and error matrix ## 2519c ## ## 2520c ############################################################## 2521c 2522c 2523c "dmdump" puts the distance matrix of the final structure 2524c into the upper half of a matrix, the distance of each atom 2525c to the centroid on the diagonal, and the individual terms 2526c of the bounds errors into the lower half of the matrix 2527c 2528c 2529 subroutine dmdump (dmd) 2530 use atoms 2531 use disgeo 2532 use iounit 2533 implicit none 2534 integer i,j 2535 real*8 sum,rgsq 2536 real*8 dist,dist2 2537 real*8 dmd(n,*) 2538 character*240 title 2539c 2540c 2541c store the final distance matrix and bound violations 2542c 2543 do i = 1, n 2544 dmd(i,i) = 0.0d0 2545 end do 2546 sum = 0.0d0 2547 do i = 1, n-1 2548 do j = i+1, n 2549 dist2 = (x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2 2550 sum = sum + dist2 2551 dmd(i,i) = dmd(i,i) + dist2 2552 dmd(j,j) = dmd(j,j) + dist2 2553 dist = sqrt(dist2) 2554 dmd(i,j) = dist 2555 if (dist .gt. dbnd(i,j)) then 2556 dmd(j,i) = dist - dbnd(i,j) 2557 else if (dist .lt. dbnd(j,i)) then 2558 dmd(j,i) = dbnd(j,i) - dist 2559 else 2560 dmd(j,i) = 0.0d0 2561 end if 2562 end do 2563 end do 2564c 2565c put the distance to the centroid on the diagonal 2566c 2567 rgsq = sum / dble(n**2) 2568 do i = 1, n 2569 dmd(i,i) = sqrt(dmd(i,i)/dble(n) - rgsq) 2570 end do 2571c 2572c write out the interatomic distance and error matrices 2573c 2574 title = 'Final Dist Matrix Above; DCM on Diag; Error Below :' 2575 call grafic (n,dmd,title) 2576 return 2577 end 2578c 2579c 2580c ################################################################# 2581c ## ## 2582c ## function initerr -- initial error function and gradient ## 2583c ## ## 2584c ################################################################# 2585c 2586c 2587c "initerr" is the initial error function and derivatives for 2588c a distance geometry embedding; it includes components from 2589c the local geometry and torsional restraint errors 2590c 2591c 2592 function initerr (xx,g) 2593 use atoms 2594 implicit none 2595 integer i,j,nvar 2596 real*8 initerr 2597 real*8 local,locerr 2598 real*8 torsion,torser 2599 real*8 xx(*) 2600 real*8 g(*) 2601 real*8, allocatable :: derivs(:,:) 2602c 2603c 2604c convert optimization parameters to atomic coordinates 2605c 2606 nvar = 0 2607 do i = 1, n 2608 nvar = nvar + 1 2609 x(i) = xx(nvar) 2610 nvar = nvar + 1 2611 y(i) = xx(nvar) 2612 nvar = nvar + 1 2613 z(i) = xx(nvar) 2614 end do 2615c 2616c perform dynamic allocation of some local arrays 2617c 2618 allocate (derivs(3,n)) 2619c 2620c zero out the values of the atomic gradient components 2621c 2622 do i = 1, n 2623 do j = 1, 3 2624 derivs(j,i) = 0.0d0 2625 end do 2626 end do 2627c 2628c compute the local geometry and the torsional 2629c components of the error function and its gradient 2630c 2631 local = locerr (derivs) 2632 torsion = torser (derivs) 2633 initerr = local + torsion 2634c 2635c store the atomic gradients as the optimization gradient 2636c 2637 nvar = 0 2638 do i = 1, n 2639 nvar = nvar + 1 2640 g(nvar) = derivs(1,i) 2641 nvar = nvar + 1 2642 g(nvar) = derivs(2,i) 2643 nvar = nvar + 1 2644 g(nvar) = derivs(3,i) 2645 end do 2646c 2647c perform deallocation of some local arrays 2648c 2649 deallocate (derivs) 2650 return 2651 end 2652c 2653c 2654c ############################################################### 2655c ## ## 2656c ## function miderr -- second error function and gradient ## 2657c ## ## 2658c ############################################################### 2659c 2660c 2661c "miderr" is the secondary error function and derivatives 2662c for a distance geometry embedding; it includes components 2663c from the distance bounds, local geometry, chirality and 2664c torsional restraint errors 2665c 2666c 2667 function miderr (xx,g) 2668 use atoms 2669 implicit none 2670 integer i,j,nvar 2671 real*8 miderr 2672 real*8 bounds,bnderr 2673 real*8 local,locerr 2674 real*8 chiral,chirer 2675 real*8 torsion,torser 2676 real*8 xx(*) 2677 real*8 g(*) 2678 real*8, allocatable :: derivs(:,:) 2679c 2680c 2681c convert optimization parameters to atomic coordinates 2682c 2683 nvar = 0 2684 do i = 1, n 2685 nvar = nvar + 1 2686 x(i) = xx(nvar) 2687 nvar = nvar + 1 2688 y(i) = xx(nvar) 2689 nvar = nvar + 1 2690 z(i) = xx(nvar) 2691 end do 2692c 2693c perform dynamic allocation of some local arrays 2694c 2695 allocate (derivs(3,n)) 2696c 2697c zero out the values of the atomic gradient components 2698c 2699 do i = 1, n 2700 do j = 1, 3 2701 derivs(j,i) = 0.0d0 2702 end do 2703 end do 2704c 2705c compute the local geometry, chirality and torsional 2706c components of the error function and its gradient 2707c 2708 bounds = bnderr (derivs) 2709 local = locerr (derivs) 2710 chiral = chirer (derivs) 2711 torsion = torser (derivs) 2712 miderr = bounds + local + chiral + torsion 2713c 2714c store the atomic gradients as the optimization gradient 2715c 2716 nvar = 0 2717 do i = 1, n 2718 nvar = nvar + 1 2719 g(nvar) = derivs(1,i) 2720 nvar = nvar + 1 2721 g(nvar) = derivs(2,i) 2722 nvar = nvar + 1 2723 g(nvar) = derivs(3,i) 2724 end do 2725c 2726c perform deallocation of some local arrays 2727c 2728 deallocate (derivs) 2729 return 2730 end 2731c 2732c 2733c ############################################################## 2734c ## ## 2735c ## function toterr -- total error function and gradient ## 2736c ## ## 2737c ############################################################## 2738c 2739c 2740c "toterr" is the error function and derivatives for a distance 2741c geometry embedding; it includes components from the distance 2742c bounds, hard sphere contacts, local geometry, chirality and 2743c torsional restraint errors 2744c 2745c 2746 function toterr (xx,g) 2747 use atoms 2748 implicit none 2749 integer i,j,nvar 2750 real*8 toterr 2751 real*8 bounds,bnderr 2752 real*8 local,locerr 2753 real*8 chiral,chirer 2754 real*8 torsion,torser 2755 real*8 contact,vdwerr 2756 real*8 xx(*) 2757 real*8 g(*) 2758 real*8, allocatable :: derivs(:,:) 2759c 2760c 2761c convert optimization parameters to atomic coordinates 2762c 2763 nvar = 0 2764 do i = 1, n 2765 nvar = nvar + 1 2766 x(i) = xx(nvar) 2767 nvar = nvar + 1 2768 y(i) = xx(nvar) 2769 nvar = nvar + 1 2770 z(i) = xx(nvar) 2771 end do 2772c 2773c perform dynamic allocation of some local arrays 2774c 2775 allocate (derivs(3,n)) 2776c 2777c zero out the values of the atomic gradient components 2778c 2779 do i = 1, n 2780 do j = 1, 3 2781 derivs(j,i) = 0.0d0 2782 end do 2783 end do 2784c 2785c compute the distance bound, vdw, chirality and torsional 2786c components of the error function and its gradient 2787c 2788 bounds = bnderr (derivs) 2789 contact = vdwerr (derivs) 2790 local = locerr (derivs) 2791 chiral = chirer (derivs) 2792 torsion = torser (derivs) 2793 toterr = bounds + contact + local + chiral + torsion 2794c 2795c store the atomic gradients as the optimization gradient 2796c 2797 nvar = 0 2798 do i = 1, n 2799 nvar = nvar + 1 2800 g(nvar) = derivs(1,i) 2801 nvar = nvar + 1 2802 g(nvar) = derivs(2,i) 2803 nvar = nvar + 1 2804 g(nvar) = derivs(3,i) 2805 end do 2806c 2807c perform deallocation of some local arrays 2808c 2809 deallocate (derivs) 2810 return 2811 end 2812c 2813c 2814c ################################################################ 2815c ## ## 2816c ## function bnderr -- computes total distance bound error ## 2817c ## ## 2818c ################################################################ 2819c 2820c 2821c "bnderr" is the distance bound error function and derivatives; 2822c this version implements the original and Havel's normalized 2823c lower bound penalty, the normalized version is preferred when 2824c lower bounds are small (as with NMR NOE restraints), the 2825c original penalty is needed if large lower bounds are present 2826c 2827c 2828 function bnderr (derivs) 2829 use atoms 2830 use restrn 2831 implicit none 2832 integer i,j,k 2833 real*8 bnderr,error,cutsq 2834 real*8 scale,chain,term 2835 real*8 gap,buffer,weigh 2836 real*8 dx,dy,dz,gx,gy,gz 2837 real*8 dstsq,bupsq,blosq 2838 real*8 derivs(3,*) 2839c 2840c 2841c zero out the distance bounds error function 2842c 2843 bnderr = 0.0d0 2844 scale = 10.0d0 2845 cutsq = 40.0d0 2846c 2847c calculate the pairwise distances between atoms 2848c 2849 do k = 1, ndfix 2850 i = idfix(1,k) 2851 j = idfix(2,k) 2852 dx = x(i) - x(j) 2853 dy = y(i) - y(j) 2854 dz = z(i) - z(j) 2855c 2856c calculate squared actual distance and bound distances; 2857c use of a small "buffer" cleans up the final error count 2858c 2859 dstsq = dx*dx + dy*dy + dz*dz 2860 gap = dfix(3,k) - dfix(2,k) 2861 buffer = 0.05d0 * min(1.0d0,gap) 2862 blosq = (dfix(2,k) + buffer)**2 2863 bupsq = (dfix(3,k) - buffer)**2 2864c 2865c error and derivatives for upper bound violation 2866c 2867 if (dstsq .gt. bupsq) then 2868 weigh = scale * dfix(1,k) 2869 term = (dstsq-bupsq) / bupsq 2870 chain = 4.0d0 * weigh * term / bupsq 2871 error = weigh * term**2 2872 gx = dx * chain 2873 gy = dy * chain 2874 gz = dz * chain 2875 bnderr = bnderr + error 2876 derivs(1,i) = derivs(1,i) + gx 2877 derivs(2,i) = derivs(2,i) + gy 2878 derivs(3,i) = derivs(3,i) + gz 2879 derivs(1,j) = derivs(1,j) - gx 2880 derivs(2,j) = derivs(2,j) - gy 2881 derivs(3,j) = derivs(3,j) - gz 2882c 2883c error and derivatives for lower bound violation 2884c 2885 else if (dstsq .lt. blosq) then 2886 weigh = scale * dfix(1,k) 2887 if (blosq .gt. cutsq) then 2888 term = (blosq-dstsq) / dstsq 2889 chain = -4.0d0 * weigh * term * (blosq/dstsq**2) 2890 else 2891 term = (blosq-dstsq) / (blosq+dstsq) 2892 chain = -8.0d0 * weigh * term * (blosq/(blosq+dstsq)**2) 2893 end if 2894 error = weigh * term**2 2895 gx = dx * chain 2896 gy = dy * chain 2897 gz = dz * chain 2898 bnderr = bnderr + error 2899 derivs(1,i) = derivs(1,i) + gx 2900 derivs(2,i) = derivs(2,i) + gy 2901 derivs(3,i) = derivs(3,i) + gz 2902 derivs(1,j) = derivs(1,j) - gx 2903 derivs(2,j) = derivs(2,j) - gy 2904 derivs(3,j) = derivs(3,j) - gz 2905 end if 2906 end do 2907 return 2908 end 2909c 2910c 2911c ############################################################### 2912c ## ## 2913c ## function vdwerr -- computes van der Waals bound error ## 2914c ## ## 2915c ############################################################### 2916c 2917c 2918c "vdwerr" is the hard sphere van der Waals bound error function 2919c and derivatives that penalizes close nonbonded contacts, 2920c pairwise neighbors are generated via the method of lights 2921c 2922c 2923 function vdwerr (derivs) 2924 use atoms 2925 use couple 2926 use disgeo 2927 use light 2928 implicit none 2929 integer i,j,k,kgy,kgz 2930 integer, allocatable :: skip(:) 2931 real*8 vdwerr,error 2932 real*8 scale,chain,term 2933 real*8 xi,yi,zi 2934 real*8 dx,dy,dz,gx,gy,gz 2935 real*8 dstsq,blosq 2936 real*8 radi,radsq 2937 real*8, allocatable :: xsort(:) 2938 real*8, allocatable :: ysort(:) 2939 real*8, allocatable :: zsort(:) 2940 real*8 derivs(3,*) 2941 logical unique 2942c 2943c 2944c zero out the distance van der Waals error function 2945c 2946 vdwerr = 0.0d0 2947 scale = 1.0d0 2948c 2949c perform dynamic allocation of some local arrays 2950c 2951 allocate (skip(n)) 2952 allocate (xsort(n)) 2953 allocate (ysort(n)) 2954 allocate (zsort(n)) 2955c 2956c transfer coordinates and zero out atoms to be skipped 2957c 2958 do i = 1, n 2959 xsort(i) = x(i) 2960 ysort(i) = y(i) 2961 zsort(i) = z(i) 2962 skip(i) = 0 2963 end do 2964c 2965c use the method of lights to generate neighbors 2966c 2967 unique = .true. 2968 call lights (vdwmax,n,xsort,ysort,zsort,unique) 2969c 2970c now, loop over all atoms computing the interactions 2971c 2972 do i = 1, n 2973 radi = georad(i) 2974 xi = xsort(rgx(i)) 2975 yi = ysort(rgy(i)) 2976 zi = zsort(rgz(i)) 2977 do j = 1, n12(i) 2978 skip(i12(j,i)) = i 2979 end do 2980 do j = 1, n13(i) 2981 skip(i13(j,i)) = i 2982 end do 2983 do j = 1, n14(i) 2984 skip(i14(j,i)) = i 2985 end do 2986 do j = kbx(i)+1, kex(i) 2987 k = locx(j) 2988 if (skip(k) .eq. i) goto 10 2989 kgy = rgy(k) 2990 if (kgy.lt.kby(i) .or. kgy.gt.key(i)) goto 10 2991 kgz = rgz(k) 2992 if (kgz.lt.kbz(i) .or. kgz.gt.kez(i)) goto 10 2993c 2994c calculate squared distances and bounds 2995c 2996 dx = xi - xsort(j) 2997 dy = yi - ysort(kgy) 2998 dz = zi - zsort(kgz) 2999 dstsq = dx*dx + dy*dy + dz*dz 3000 radsq = (radi + georad(k))**2 3001 blosq = min(dbnd(k,i),dbnd(i,k),radsq) 3002c 3003c error and derivatives for lower bound violation 3004c 3005 if (dstsq .lt. blosq) then 3006 term = (blosq-dstsq) / (blosq+dstsq) 3007 chain = -8.0d0 * scale * term * (blosq/(blosq+dstsq)**2) 3008 error = scale * term**2 3009 gx = dx * chain 3010 gy = dy * chain 3011 gz = dz * chain 3012 vdwerr = vdwerr + error 3013 derivs(1,i) = derivs(1,i) + gx 3014 derivs(2,i) = derivs(2,i) + gy 3015 derivs(3,i) = derivs(3,i) + gz 3016 derivs(1,k) = derivs(1,k) - gx 3017 derivs(2,k) = derivs(2,k) - gy 3018 derivs(3,k) = derivs(3,k) - gz 3019 end if 3020 10 continue 3021 end do 3022 end do 3023c 3024c perform deallocation of some local arrays 3025c 3026 deallocate (skip) 3027 deallocate (xsort) 3028 deallocate (ysort) 3029 deallocate (zsort) 3030 return 3031 end 3032c 3033c 3034c ################################################################ 3035c ## ## 3036c ## function locerr -- computes local geometry error value ## 3037c ## ## 3038c ################################################################ 3039c 3040c 3041c "locerr" is the local geometry error function and derivatives 3042c including the 1-2, 1-3 and 1-4 distance bound restraints 3043c 3044c 3045 function locerr (derivs) 3046 use angbnd 3047 use atoms 3048 use bndstr 3049 use disgeo 3050 use tors 3051 implicit none 3052 integer i,ia,ib,ic,id 3053 real*8 locerr,error 3054 real*8 scale,chain,term 3055 real*8 dx,dy,dz,gx,gy,gz 3056 real*8 dstsq,bupsq,blosq 3057 real*8 derivs(3,*) 3058c 3059c 3060c zero out the local geometry error function 3061c 3062 locerr = 0.0d0 3063 scale = 10.0d0 3064c 3065c calculate the bounds error for bond length distances 3066c 3067 do i = 1, nbond 3068 ia = min(ibnd(1,i),ibnd(2,i)) 3069 ib = max(ibnd(1,i),ibnd(2,i)) 3070 dx = x(ia) - x(ib) 3071 dy = y(ia) - y(ib) 3072 dz = z(ia) - z(ib) 3073 dstsq = dx*dx + dy*dy + dz*dz 3074 bupsq = dbnd(ia,ib) 3075 blosq = dbnd(ib,ia) 3076 if (dstsq .gt. bupsq) then 3077 term = (dstsq-bupsq) / bupsq 3078 chain = 4.0d0 * scale * term / bupsq 3079 error = scale * term**2 3080 gx = dx * chain 3081 gy = dy * chain 3082 gz = dz * chain 3083 locerr = locerr + error 3084 derivs(1,ia) = derivs(1,ia) + gx 3085 derivs(2,ia) = derivs(2,ia) + gy 3086 derivs(3,ia) = derivs(3,ia) + gz 3087 derivs(1,ib) = derivs(1,ib) - gx 3088 derivs(2,ib) = derivs(2,ib) - gy 3089 derivs(3,ib) = derivs(3,ib) - gz 3090 else if (dstsq .lt. blosq) then 3091 term = (blosq-dstsq) / (blosq+dstsq) 3092 chain = -8.0d0 * scale * term * (blosq/(blosq+dstsq)**2) 3093 error = scale * term**2 3094 gx = dx * chain 3095 gy = dy * chain 3096 gz = dz * chain 3097 locerr = locerr + error 3098 derivs(1,ia) = derivs(1,ia) + gx 3099 derivs(2,ia) = derivs(2,ia) + gy 3100 derivs(3,ia) = derivs(3,ia) + gz 3101 derivs(1,ib) = derivs(1,ib) - gx 3102 derivs(2,ib) = derivs(2,ib) - gy 3103 derivs(3,ib) = derivs(3,ib) - gz 3104 end if 3105 end do 3106c 3107c calculate the bounds error for the bond angle distances 3108c 3109 do i = 1, nangle 3110 ia = min(iang(1,i),iang(3,i)) 3111 ic = max(iang(1,i),iang(3,i)) 3112 dx = x(ia) - x(ic) 3113 dy = y(ia) - y(ic) 3114 dz = z(ia) - z(ic) 3115 dstsq = dx*dx + dy*dy + dz*dz 3116 bupsq = dbnd(ia,ic) 3117 blosq = dbnd(ic,ia) 3118 if (dstsq .gt. bupsq) then 3119 term = (dstsq-bupsq) / bupsq 3120 chain = 4.0d0 * scale * term / bupsq 3121 error = scale * term**2 3122 gx = dx * chain 3123 gy = dy * chain 3124 gz = dz * chain 3125 locerr = locerr + error 3126 derivs(1,ia) = derivs(1,ia) + gx 3127 derivs(2,ia) = derivs(2,ia) + gy 3128 derivs(3,ia) = derivs(3,ia) + gz 3129 derivs(1,ic) = derivs(1,ic) - gx 3130 derivs(2,ic) = derivs(2,ic) - gy 3131 derivs(3,ic) = derivs(3,ic) - gz 3132 else if (dstsq .lt. blosq) then 3133 term = (blosq-dstsq) / (blosq+dstsq) 3134 chain = -8.0d0 * scale * term * (blosq/(blosq+dstsq)**2) 3135 error = scale * term**2 3136 gx = dx * chain 3137 gy = dy * chain 3138 gz = dz * chain 3139 locerr = locerr + error 3140 derivs(1,ia) = derivs(1,ia) + gx 3141 derivs(2,ia) = derivs(2,ia) + gy 3142 derivs(3,ia) = derivs(3,ia) + gz 3143 derivs(1,ic) = derivs(1,ic) - gx 3144 derivs(2,ic) = derivs(2,ic) - gy 3145 derivs(3,ic) = derivs(3,ic) - gz 3146 end if 3147 end do 3148c 3149c calculate the bounds error for the torsion angle distances 3150c 3151 do i = 1, ntors 3152 ia = min(itors(1,i),itors(4,i)) 3153 id = max(itors(1,i),itors(4,i)) 3154 dx = x(ia) - x(id) 3155 dy = y(ia) - y(id) 3156 dz = z(ia) - z(id) 3157 dstsq = dx*dx + dy*dy + dz*dz 3158 bupsq = dbnd(ia,id) 3159 blosq = dbnd(id,ia) 3160 if (dstsq .gt. bupsq) then 3161 term = (dstsq-bupsq) / bupsq 3162 chain = 4.0d0 * scale * term / bupsq 3163 error = scale * term**2 3164 gx = dx * chain 3165 gy = dy * chain 3166 gz = dz * chain 3167 locerr = locerr + error 3168 derivs(1,ia) = derivs(1,ia) + gx 3169 derivs(2,ia) = derivs(2,ia) + gy 3170 derivs(3,ia) = derivs(3,ia) + gz 3171 derivs(1,id) = derivs(1,id) - gx 3172 derivs(2,id) = derivs(2,id) - gy 3173 derivs(3,id) = derivs(3,id) - gz 3174 else if (dstsq .lt. blosq) then 3175 term = (blosq-dstsq) / (blosq+dstsq) 3176 chain = -8.0d0 * scale * term * (blosq/(blosq+dstsq)**2) 3177 error = scale * term**2 3178 gx = dx * chain 3179 gy = dy * chain 3180 gz = dz * chain 3181 locerr = locerr + error 3182 derivs(1,ia) = derivs(1,ia) + gx 3183 derivs(2,ia) = derivs(2,ia) + gy 3184 derivs(3,ia) = derivs(3,ia) + gz 3185 derivs(1,id) = derivs(1,id) - gx 3186 derivs(2,id) = derivs(2,id) - gy 3187 derivs(3,id) = derivs(3,id) - gz 3188 end if 3189 end do 3190 return 3191 end 3192c 3193c 3194c ############################################################## 3195c ## ## 3196c ## function chirer -- computes chirality error function ## 3197c ## ## 3198c ############################################################## 3199c 3200c 3201c "chirer" computes the chirality error and its derivatives 3202c with respect to atomic Cartesian coordinates as a sum the 3203c squares of deviations of chiral volumes from target values 3204c 3205c 3206 function chirer (derivs) 3207 use atoms 3208 use restrn 3209 implicit none 3210 integer i,ia,ib,ic,id 3211 real*8 chirer,error,scale 3212 real*8 vol,dt,dedt 3213 real*8 xad,yad,zad 3214 real*8 xbd,ybd,zbd 3215 real*8 xcd,ycd,zcd 3216 real*8 c1,c2,c3 3217 real*8 dedxia,dedyia,dedzia 3218 real*8 dedxib,dedyib,dedzib 3219 real*8 dedxic,dedyic,dedzic 3220 real*8 dedxid,dedyid,dedzid 3221 real*8 derivs(3,*) 3222c 3223c 3224c zero the chirality restraint error function 3225c 3226 chirer = 0.0d0 3227 scale = 0.1d0 3228c 3229c find signed volume value and compute chirality error 3230c 3231 do i = 1, nchir 3232 ia = ichir(1,i) 3233 ib = ichir(2,i) 3234 ic = ichir(3,i) 3235 id = ichir(4,i) 3236 xad = x(ia) - x(id) 3237 yad = y(ia) - y(id) 3238 zad = z(ia) - z(id) 3239 xbd = x(ib) - x(id) 3240 ybd = y(ib) - y(id) 3241 zbd = z(ib) - z(id) 3242 xcd = x(ic) - x(id) 3243 ycd = y(ic) - y(id) 3244 zcd = z(ic) - z(id) 3245 c1 = ybd*zcd - zbd*ycd 3246 c2 = ycd*zad - zcd*yad 3247 c3 = yad*zbd - zad*ybd 3248 vol = xad*c1 + xbd*c2 + xcd*c3 3249 dt = vol - chir(2,i) 3250 error = scale * dt**2 3251 dedt = 2.0d0 * scale * dt 3252c 3253c compute derivative components for this interaction 3254c 3255 dedxia = dedt * (ybd*zcd - zbd*ycd) 3256 dedyia = dedt * (zbd*xcd - xbd*zcd) 3257 dedzia = dedt * (xbd*ycd - ybd*xcd) 3258 dedxib = dedt * (zad*ycd - yad*zcd) 3259 dedyib = dedt * (xad*zcd - zad*xcd) 3260 dedzib = dedt * (yad*xcd - xad*ycd) 3261 dedxic = dedt * (yad*zbd - zad*ybd) 3262 dedyic = dedt * (zad*xbd - xad*zbd) 3263 dedzic = dedt * (xad*ybd - yad*xbd) 3264 dedxid = -dedxia - dedxib - dedxic 3265 dedyid = -dedyia - dedyib - dedyic 3266 dedzid = -dedzia - dedzib - dedzic 3267c 3268c increment the chirality restraint error and derivatives 3269c 3270 chirer = chirer + error 3271 derivs(1,ia) = derivs(1,ia) + dedxia 3272 derivs(2,ia) = derivs(2,ia) + dedyia 3273 derivs(3,ia) = derivs(3,ia) + dedzia 3274 derivs(1,ib) = derivs(1,ib) + dedxib 3275 derivs(2,ib) = derivs(2,ib) + dedyib 3276 derivs(3,ib) = derivs(3,ib) + dedzib 3277 derivs(1,ic) = derivs(1,ic) + dedxic 3278 derivs(2,ic) = derivs(2,ic) + dedyic 3279 derivs(3,ic) = derivs(3,ic) + dedzic 3280 derivs(1,id) = derivs(1,id) + dedxid 3281 derivs(2,id) = derivs(2,id) + dedyid 3282 derivs(3,id) = derivs(3,id) + dedzid 3283 end do 3284 return 3285 end 3286c 3287c 3288c ############################################################## 3289c ## ## 3290c ## function torser -- computes torsional error function ## 3291c ## ## 3292c ############################################################## 3293c 3294c 3295c "torser" computes the torsional error function and its first 3296c derivatives with respect to the atomic Cartesian coordinates 3297c based on the deviation of specified torsional angles from 3298c desired values, the contained bond angles are also restrained 3299c to avoid a numerical instability 3300c 3301c 3302 function torser (derivs) 3303 use atoms 3304 use couple 3305 use disgeo 3306 use math 3307 use refer 3308 use restrn 3309 implicit none 3310 integer i,j,k,iref 3311 integer ia,ib,ic,id 3312 real*8 torser,error,force 3313 real*8 dt,deddt,dedphi 3314 real*8 angle,target,scale 3315 real*8 xia,yia,zia 3316 real*8 xib,yib,zib 3317 real*8 xic,yic,zic 3318 real*8 xid,yid,zid 3319 real*8 xria,yria,zria 3320 real*8 xrib,yrib,zrib 3321 real*8 xric,yric,zric 3322 real*8 xrid,yrid,zrid 3323 real*8 rba,rcb,rdc 3324 real*8 dot,cosine,sine 3325 real*8 rrba,rrcb,rrdc 3326 real*8 rrca,rrdb 3327 real*8 bndfac,angfac 3328 real*8 xp,yp,zp,rp 3329 real*8 terma,termb 3330 real*8 termc,termd 3331 real*8 angmax,angmin 3332 real*8 cosmax,cosmin 3333 real*8 bamax,bamin 3334 real*8 cbmax,cbmin 3335 real*8 dcmax,dcmin 3336 real*8 camax,camin 3337 real*8 dbmax,dbmin 3338 real*8 xba,yba,zba 3339 real*8 xcb,ycb,zcb 3340 real*8 xdc,ydc,zdc 3341 real*8 xca,yca,zca 3342 real*8 xdb,ydb,zdb 3343 real*8 xt,yt,zt,rt2 3344 real*8 xu,yu,zu,ru2 3345 real*8 xtu,ytu,ztu,rtru 3346 real*8 tf1,tf2,t1,t2 3347 real*8 dedxt,dedyt,dedzt 3348 real*8 dedxu,dedyu,dedzu 3349 real*8 dedxia,dedyia,dedzia 3350 real*8 dedxib,dedyib,dedzib 3351 real*8 dedxic,dedyic,dedzic 3352 real*8 dedxid,dedyid,dedzid 3353 real*8 derivs(3,*) 3354 logical bonded 3355c 3356c 3357c zero the torsional restraint error function 3358c 3359 torser = 0.0d0 3360 scale = 0.01d0 3361c 3362c compute error value and derivs for torsional restraints 3363c 3364 do i = 1, ntfix 3365 ia = itfix(1,i) 3366 ib = itfix(2,i) 3367 ic = itfix(3,i) 3368 id = itfix(4,i) 3369 xia = x(ia) 3370 yia = y(ia) 3371 zia = z(ia) 3372 xib = x(ib) 3373 yib = y(ib) 3374 zib = z(ib) 3375 xic = x(ic) 3376 yic = y(ic) 3377 zic = z(ic) 3378 xid = x(id) 3379 yid = y(id) 3380 zid = z(id) 3381 xba = xib - xia 3382 yba = yib - yia 3383 zba = zib - zia 3384 xcb = xic - xib 3385 ycb = yic - yib 3386 zcb = zic - zib 3387 xdc = xid - xic 3388 ydc = yid - yic 3389 zdc = zid - zic 3390c 3391c find the actual distances between the four atoms 3392c 3393 rba = sqrt(xba*xba + yba*yba + zba*zba) 3394 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb) 3395 rdc = sqrt(xdc*xdc + ydc*ydc + zdc*zdc) 3396c 3397c see if the torsion involves four directly bonded atoms 3398c 3399 k = 0 3400 do j = 1, n12(ia) 3401 if (i12(j,ia) .eq. ib) k = k + 1 3402 end do 3403 do j = 1, n12(ib) 3404 if (i12(j,ib) .eq. ic) k = k + 1 3405 end do 3406 do j = 1, n12(ic) 3407 if (i12(j,ic) .eq. id) k = k + 1 3408 end do 3409 if (k .eq. 3) then 3410 bonded = .true. 3411 else 3412 bonded = .false. 3413 end if 3414c 3415c get maximum and minimum distances from distance matrix 3416c 3417 if (bonded) then 3418 bamax = sqrt(dbnd(min(ib,ia),max(ib,ia))) 3419 bamin = sqrt(dbnd(max(ib,ia),min(ib,ia))) 3420 cbmax = sqrt(dbnd(min(ic,ib),max(ic,ib))) 3421 cbmin = sqrt(dbnd(max(ic,ib),min(ic,ib))) 3422 dcmax = sqrt(dbnd(min(id,ic),max(id,ic))) 3423 dcmin = sqrt(dbnd(max(id,ic),min(id,ic))) 3424 camax = sqrt(dbnd(min(ic,ia),max(ic,ia))) 3425 camin = sqrt(dbnd(max(ic,ia),min(ic,ia))) 3426 dbmax = sqrt(dbnd(min(id,ib),max(id,ib))) 3427 dbmin = sqrt(dbnd(max(id,ib),min(id,ib))) 3428c 3429c get maximum and minimum distances from input structure 3430c 3431 else 3432 iref = 1 3433 xria = xref(ia,iref) 3434 yria = yref(ia,iref) 3435 zria = zref(ia,iref) 3436 xrib = xref(ib,iref) 3437 yrib = yref(ib,iref) 3438 zrib = zref(ib,iref) 3439 xric = xref(ic,iref) 3440 yric = yref(ic,iref) 3441 zric = zref(ic,iref) 3442 xrid = xref(id,iref) 3443 yrid = yref(id,iref) 3444 zrid = zref(id,iref) 3445 rrba = sqrt((xrib-xria)**2+(yrib-yria)**2+(zrib-zria)**2) 3446 rrcb = sqrt((xric-xrib)**2+(yric-yrib)**2+(zric-zrib)**2) 3447 rrdc = sqrt((xrid-xric)**2+(yrid-yric)**2+(zrid-zric)**2) 3448 rrca = sqrt((xric-xria)**2+(yric-yria)**2+(zric-zria)**2) 3449 rrdb = sqrt((xrid-xrib)**2+(yrid-yrib)**2+(zrid-zrib)**2) 3450 bndfac = 0.05d0 3451 angfac = 0.05d0 3452 bamax = (1.0d0 + bndfac) * rrba 3453 bamin = (1.0d0 - bndfac) * rrba 3454 cbmax = (1.0d0 + bndfac) * rrcb 3455 cbmin = (1.0d0 - bndfac) * rrcb 3456 dcmax = (1.0d0 + bndfac) * rrdc 3457 dcmin = (1.0d0 - bndfac) * rrdc 3458 camax = (1.0d0 + angfac) * rrca 3459 camin = (1.0d0 - angfac) * rrca 3460 dbmax = (1.0d0 + angfac) * rrdb 3461 dbmin = (1.0d0 - angfac) * rrdb 3462 end if 3463c 3464c compute the ia-ib-ic bond angle and any error 3465c 3466 dot = xba*xcb + yba*ycb + zba*zcb 3467 cosine = -dot / (rba*rcb) 3468 cosine = min(1.0d0,max(-1.0d0,cosine)) 3469 angle = radian * acos(cosine) 3470 cosmax = (bamin**2+cbmin**2-camax**2) / (2.0d0*bamin*cbmin) 3471 cosmax = min(1.0d0,max(-1.0d0,cosmax)) 3472 angmax = radian * acos(cosmax) 3473 cosmin = (bamax**2+cbmax**2-camin**2) / (2.0d0*bamax*cbmax) 3474 cosmin = min(1.0d0,max(-1.0d0,cosmin)) 3475 angmin = radian * acos(cosmin) 3476 if (angle .gt. angmax) then 3477 dt = angle - angmax 3478 else if (angle .lt. angmin) then 3479 dt = angle - angmin 3480 else 3481 dt = 0.0d0 3482 end if 3483 error = scale * dt**2 3484 deddt = 2.0d0 * radian * scale * dt 3485c 3486c compute derivative components for this interaction 3487c 3488 xp = zcb*yba - ycb*zba 3489 yp = xcb*zba - zcb*xba 3490 zp = ycb*xba - xcb*yba 3491 rp = sqrt(xp*xp + yp*yp + zp*zp) 3492 if (rp .ne. 0.0d0) then 3493 terma = -deddt / (rba*rba*rp) 3494 termc = deddt / (rcb*rcb*rp) 3495 dedxia = terma * (zba*yp-yba*zp) 3496 dedyia = terma * (xba*zp-zba*xp) 3497 dedzia = terma * (yba*xp-xba*yp) 3498 dedxic = termc * (ycb*zp-zcb*yp) 3499 dedyic = termc * (zcb*xp-xcb*zp) 3500 dedzic = termc * (xcb*yp-ycb*xp) 3501 dedxib = -dedxia - dedxic 3502 dedyib = -dedyia - dedyic 3503 dedzib = -dedzia - dedzic 3504c 3505c increment the bond angle restraint error and derivatives 3506c 3507 torser = torser + error 3508 derivs(1,ia) = derivs(1,ia) + dedxia 3509 derivs(2,ia) = derivs(2,ia) + dedyia 3510 derivs(3,ia) = derivs(3,ia) + dedzia 3511 derivs(1,ib) = derivs(1,ib) + dedxib 3512 derivs(2,ib) = derivs(2,ib) + dedyib 3513 derivs(3,ib) = derivs(3,ib) + dedzib 3514 derivs(1,ic) = derivs(1,ic) + dedxic 3515 derivs(2,ic) = derivs(2,ic) + dedyic 3516 derivs(3,ic) = derivs(3,ic) + dedzic 3517 end if 3518c 3519c compute the ib-ic-id bond angle and any error 3520c 3521 dot = xdc*xcb + ydc*ycb + zdc*zcb 3522 cosine = -dot / (rdc*rcb) 3523 cosine = min(1.0d0,max(-1.0d0,cosine)) 3524 angle = radian * acos(cosine) 3525 cosmax = (dcmin**2+cbmin**2-dbmax**2) / (2.0d0*dcmin*cbmin) 3526 cosmax = min(1.0d0,max(-1.0d0,cosmax)) 3527 angmax = radian * acos(cosmax) 3528 cosmin = (dcmax**2+cbmax**2-dbmin**2) / (2.0d0*dcmax*cbmax) 3529 cosmax = min(1.0d0,max(-1.0d0,cosmin)) 3530 angmin = radian * acos(cosmin) 3531 if (angle .gt. angmax) then 3532 dt = angle - angmax 3533 else if (angle .lt. angmin) then 3534 dt = angle - angmin 3535 else 3536 dt = 0.0d0 3537 end if 3538 error = scale * dt**2 3539 deddt = 2.0d0 * radian * scale * dt 3540c 3541c compute derivative components for this interaction 3542c 3543 xp = zdc*ycb - ydc*zcb 3544 yp = xdc*zcb - zdc*xcb 3545 zp = ydc*xcb - xdc*ycb 3546 rp = sqrt(xp*xp + yp*yp + zp*zp) 3547 if (rp .ne. 0.0d0) then 3548 termb = -deddt / (rcb*rcb*rp) 3549 termd = deddt / (rdc*rdc*rp) 3550 dedxib = termb * (zcb*yp-ycb*zp) 3551 dedyib = termb * (xcb*zp-zcb*xp) 3552 dedzib = termb * (ycb*xp-xcb*yp) 3553 dedxid = termd * (ydc*zp-zdc*yp) 3554 dedyid = termd * (zdc*xp-xdc*zp) 3555 dedzid = termd * (xdc*yp-ydc*xp) 3556 dedxic = -dedxib - dedxid 3557 dedyic = -dedyib - dedyid 3558 dedzic = -dedzib - dedzid 3559c 3560c increment the bond angle restraint error and derivatives 3561c 3562 torser = torser + error 3563 derivs(1,ib) = derivs(1,ib) + dedxib 3564 derivs(2,ib) = derivs(2,ib) + dedyib 3565 derivs(3,ib) = derivs(3,ib) + dedzib 3566 derivs(1,ic) = derivs(1,ic) + dedxic 3567 derivs(2,ic) = derivs(2,ic) + dedyic 3568 derivs(3,ic) = derivs(3,ic) + dedzic 3569 derivs(1,id) = derivs(1,id) + dedxid 3570 derivs(2,id) = derivs(2,id) + dedyid 3571 derivs(3,id) = derivs(3,id) + dedzid 3572 end if 3573c 3574c compute the value of the ia-ib-ic-id torsional angle 3575c 3576 xt = yba*zcb - ycb*zba 3577 yt = zba*xcb - zcb*xba 3578 zt = xba*ycb - xcb*yba 3579 xu = ycb*zdc - ydc*zcb 3580 yu = zcb*xdc - zdc*xcb 3581 zu = xcb*ydc - xdc*ycb 3582 xtu = yt*zu - yu*zt 3583 ytu = zt*xu - zu*xt 3584 ztu = xt*yu - xu*yt 3585 rt2 = xt*xt + yt*yt + zt*zt 3586 ru2 = xu*xu + yu*yu + zu*zu 3587 rtru = sqrt(rt2 * ru2) 3588 if (rtru .ne. 0.0d0) then 3589 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb) 3590 cosine = (xt*xu + yt*yu + zt*zu) / rtru 3591 sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru) 3592 cosine = min(1.0d0,max(-1.0d0,cosine)) 3593 angle = radian * acos(cosine) 3594 if (sine .lt. 0.0d0) angle = -angle 3595c 3596c calculate the torsional restraint error for this angle 3597c 3598 force = tfix(1,i) 3599 tf1 = tfix(2,i) 3600 tf2 = tfix(3,i) 3601 if (angle.gt.tf1 .and. angle.lt.tf2) then 3602 target = angle 3603 else if (angle.gt.tf1 .and. tf1.gt.tf2) then 3604 target = angle 3605 else if (angle.lt.tf2 .and. tf1.gt.tf2) then 3606 target = angle 3607 else 3608 t1 = angle - tf1 3609 t2 = angle - tf2 3610 if (t1 .gt. 180.0d0) then 3611 t1 = t1 - 360.0d0 3612 else if (t1 .lt. -180.0d0) then 3613 t1 = t1 + 360.0d0 3614 end if 3615 if (t2 .gt. 180.0d0) then 3616 t2 = t2 - 360.0d0 3617 else if (t2 .lt. -180.0d0) then 3618 t2 = t2 + 360.0d0 3619 end if 3620 if (abs(t1) .lt. abs(t2)) then 3621 target = tf1 3622 else 3623 target = tf2 3624 end if 3625 end if 3626 dt = angle - target 3627 if (dt .gt. 180.0d0) then 3628 dt = dt - 360.0d0 3629 else if (dt .lt. -180.0d0) then 3630 dt = dt + 360.0d0 3631 end if 3632 error = scale * force * dt**2 3633 dedphi = 2.0d0 * radian * scale * force * dt 3634c 3635c chain rule terms for first derivative components 3636c 3637 xca = xic - xia 3638 yca = yic - yia 3639 zca = zic - zia 3640 xdb = xid - xib 3641 ydb = yid - yib 3642 zdb = zid - zib 3643 dedxt = dedphi * (yt*zcb - ycb*zt) / (rt2*rcb) 3644 dedyt = dedphi * (zt*xcb - zcb*xt) / (rt2*rcb) 3645 dedzt = dedphi * (xt*ycb - xcb*yt) / (rt2*rcb) 3646 dedxu = -dedphi * (yu*zcb - ycb*zu) / (ru2*rcb) 3647 dedyu = -dedphi * (zu*xcb - zcb*xu) / (ru2*rcb) 3648 dedzu = -dedphi * (xu*ycb - xcb*yu) / (ru2*rcb) 3649c 3650c compute first derivative components for torsion angle 3651c 3652 dedxia = zcb*dedyt - ycb*dedzt 3653 dedyia = xcb*dedzt - zcb*dedxt 3654 dedzia = ycb*dedxt - xcb*dedyt 3655 dedxib = yca*dedzt - zca*dedyt + zdc*dedyu - ydc*dedzu 3656 dedyib = zca*dedxt - xca*dedzt + xdc*dedzu - zdc*dedxu 3657 dedzib = xca*dedyt - yca*dedxt + ydc*dedxu - xdc*dedyu 3658 dedxic = zba*dedyt - yba*dedzt + ydb*dedzu - zdb*dedyu 3659 dedyic = xba*dedzt - zba*dedxt + zdb*dedxu - xdb*dedzu 3660 dedzic = yba*dedxt - xba*dedyt + xdb*dedyu - ydb*dedxu 3661 dedxid = zcb*dedyu - ycb*dedzu 3662 dedyid = xcb*dedzu - zcb*dedxu 3663 dedzid = ycb*dedxu - xcb*dedyu 3664c 3665c increment the torsional restraint error and derivatives 3666c 3667 torser = torser + error 3668 derivs(1,ia) = derivs(1,ia) + dedxia 3669 derivs(2,ia) = derivs(2,ia) + dedyia 3670 derivs(3,ia) = derivs(3,ia) + dedzia 3671 derivs(1,ib) = derivs(1,ib) + dedxib 3672 derivs(2,ib) = derivs(2,ib) + dedyib 3673 derivs(3,ib) = derivs(3,ib) + dedzib 3674 derivs(1,ic) = derivs(1,ic) + dedxic 3675 derivs(2,ic) = derivs(2,ic) + dedyic 3676 derivs(3,ic) = derivs(3,ic) + dedzic 3677 derivs(1,id) = derivs(1,id) + dedxid 3678 derivs(2,id) = derivs(2,id) + dedyid 3679 derivs(3,id) = derivs(3,id) + dedzid 3680 end if 3681 end do 3682 return 3683 end 3684