1 subroutine reddis(vr) 2 implicit real (a-h,o-z), integer (i-n) 3 logical box,cell,fast 4 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 5 dimension vr(3) 6 7 8 9 do i=1,3 10 vr(i) = mod(vr(i),abc(i)) 11 if (vr(i).gt.abc2(i)) then 12 vr(i) = vr(i) - abc(i) 13 endif 14 if (vr(i).lt.-abc2(i)) then 15 vr(i) = vr(i) + abc(i) 16 endif 17 end do 18 19 return 20 end 21 22 subroutine rddisr(vr) 23 implicit real (a-h,o-z), integer (i-n) 24 real abc,abc2,angles 25 logical box,cell,fast 26 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 27 dimension vr(3) 28 29 do i=1,3 30 do while (vr(i).gt.real(abc2(i))) 31 vr(i) = vr(i) - real(abc(i)) 32 end do 33 do while (vr(i).lt.-real(abc2(i))) 34 vr(i) = vr(i) + real(abc(i)) 35 end do 36 end do 37 38 return 39 end 40 41 subroutine appbnd(coo,ityp) 42 implicit real (a-h,o-z), integer (i-n) 43 integer*2 ityp 44 logical box,cell,fast,outbox 45 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 46 common /athlp/ iatoms, mxnat 47 dimension coo(3,*),cw(3,3),vec(3),ityp(*) 48 49c put water molecules which are entirely out of the box 50c back in the box 51 52 do i=1,iatoms 53 54 if (ityp(i) .eq.649.and. 55 & ityp(i+1).eq.650.and. 56 & ityp(i+2).eq.650) then 57 58 outbox = .true. 59 60 do j=1,3 61 vec(j) = 0.0e0 62 end do 63 64 do k=1,3 65 66 l = 0 67 68 do j=1,3 69 cw(j,k) = coo(j,i+k-1) 70 if (cw(j,k).gt.abc2(j)) then 71 l = l + 1 72 if (k.eq.1) vec(j) = -abc(j) 73 endif 74 if (cw(j,k).lt.-abc2(j)) then 75 l = l + 1 76 if (k.eq.1) vec(j) = abc(j) 77 endif 78 end do 79 80 if (l.eq.0) outbox = .false. 81 82 end do 83 84c water out of the box, put it back 85 86 if (outbox) then 87 do k=1,3 88 do j=1,3 89 coo(j,i+k-1) = cw(j,k) + vec(j) 90 end do 91 end do 92 endif 93 94 endif 95 end do 96 97 return 98 end 99 100 subroutine makbod(coo) 101 implicit real (a-h,o-z), integer (i-n) 102 logical box,cell,fast 103 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 104 dimension coo(3,*),vec(3) 105 106c get largest diameter protein and add a default 107c set abc and abc2 108c set protein in center of the box 109 110 offs = 14.0e0 111 112 call docnt(vec,coo) 113 114 do i=1,3 115 abc(i) = 2.0e0*vec(i) + offs 116 abc2(i) = 0.5e0*abc(i) 117 end do 118 119 return 120 end 121 122 subroutine allbox 123 implicit real (a-h,o-z), integer (i-n) 124 parameter (mxion=2000) 125 logical box,cell,fast 126 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 127 common /h2oer/ numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat 128 common /athlp/ iatoms, mxnat 129 dimension nbox(3) 130 131 watd = 18.620e0 132 nwater = 648/3 133 134 do i=1,3 135 t = abc(i)/watd 136 nbox(i) = int(t) + 1 137 end do 138 139 newat = nbox(1)*nbox(2)*nbox(3)*nwater*3 140 141 call allcoo(newat) 142 if (nion.gt.0) call allq 143 144 return 145 end 146 147 subroutine filbod(water,coo,iconn,iresid,ityp, 148 & nac,iac,nad, 149 & nbnd,ibnd,bl,bk, 150 & nang,iang,ango,ak,q,iopt,iwtpr) 151 implicit real (a-h,o-z), integer (i-n) 152 parameter (mxcon=10) 153 parameter (mxac=3*mxcon) 154 parameter (numres=50000) 155 parameter (mxion=2000) 156 logical box,cell,fast,chkwat,chkbox,chkion 157 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 158 common /athlp/ iatoms, mxnat 159 common /residu/ qtot,ihsres, 160 & ires(numres),ibeg(numres),iend(numres) 161 common /h2oer/ numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat 162 integer*2 ityp 163 dimension water(3,*), coo(3,*), nbox(3), doff(3) 164 dimension ityp(*),iconn(mxcon+1,*),iresid(*) 165 dimension nac(*),nad(*),iac(mxac,*) 166 dimension ibnd(2,*), bl(*), bk(*) 167 dimension iang(3,*),ango(*),ak(*) 168 dimension q(*),iopt(*),iwtpr(*),cwat(3,3),ibox(3) 169 170 watd = 18.620e0 171 nwater = 648/3 172 173c 18.620 ang box 174c take prefilled water box and tile the real box with this 175c sub box 176c remove water that overlap with protein 177 178 do i=1,3 179 t = abc(i)/watd 180 nbox(i) = int(t) + 1 181 end do 182 183 ishoh = ires(ihsres) 184 if (ishoh.gt.0) then 185 ishoh = -4 186 else 187 ishoh = ishoh - 1 188 endif 189 190 natnow = iatoms 191 ioff = iatoms 192 numwat = 0 193 nwat = 0 194 195 ibox(1) = 0 196 197 do i=1,nbox(1) 198 if (i.eq.nbox(1)) ibox(1) = 1 199 ibox(2) = 0 200 201 do j=1,nbox(2) 202 if (j.eq.nbox(2)) ibox(2) = 1 203 ibox(3) = 0 204 205 do k=1,nbox(3) 206 if (k.eq.nbox(3)) ibox(3) = 1 207 208 doff(1) = -abc2(1) + 0.5e0*watd + watd*dble(i-1) 209 doff(2) = -abc2(2) + 0.5e0*watd + watd*dble(j-1) 210 doff(3) = -abc2(3) + 0.5e0*watd + watd*dble(k-1) 211 212 chkbox = .false. 213 if (i.eq.nbox(1).or.j.eq.nbox(2).or.k.eq.nbox(3)) 214 & chkbox = .true. 215 216 do n=1,nwater 217 218 in = (n-1)*3 219 220 do m=1,3 221 cwat(m,1) = water(m,in+1) + doff(m) 222 cwat(m,2) = water(m,in+2) + doff(m) 223 cwat(m,3) = water(m,in+3) + doff(m) 224 end do 225 226 if (chkwat(cwat,coo,ityp,iopt,ioff,ibox, 227 & iwt,chkbox)) then 228 229 230 nwat = nwat + 1 231 232 if (.not.chkion(nwat)) then 233 234c ADD WATER 235 236 numwat = numwat + 1 237 iwtpr(numwat) = iwt 238 239 if (ihsres.lt.numres) then 240 ihsres = ihsres + 1 241 endif 242 243 ires(ihsres) = ishoh 244 ibeg(ihsres) = ioff+1 245 iend(ihsres) = ioff+3 246 iresid(ioff+1) = ishoh 247 iresid(ioff+2) = ishoh 248 iresid(ioff+3) = ishoh 249 ishoh = ishoh - 1 250 251 do m=1,3 252 coo(m,ioff+1) = cwat(m,1) 253 coo(m,ioff+2) = cwat(m,2) 254 coo(m,ioff+3) = cwat(m,3) 255 end do 256 257 ityp(ioff+1) = 649 258 ityp(ioff+2) = 650 259 ityp(ioff+3) = 650 260 261 iwtpr(numwat) = 0 262 263 iconn(1,ioff+1) = 2 264 iconn(2,ioff+1) = ioff+2 265 iconn(3,ioff+1) = ioff+3 266 267 iconn(1,ioff+2) = 1 268 iconn(2,ioff+2) = ioff+1 269 270 iconn(1,ioff+3) = 1 271 iconn(2,ioff+3) = ioff+1 272 273 ibnd(1,nbnd+1) = ioff+1 274 ibnd(2,nbnd+1) = ioff+2 275 bl(nbnd+1) = 0.9572e0 276 bk(nbnd+1) = 553.0e0 277 278 ibnd(1,nbnd+2) = ioff+1 279 ibnd(2,nbnd+2) = ioff+3 280 bl(nbnd+2) = 0.9572e0 281 bk(nbnd+2) = 553.0e0 282 nbnd = nbnd + 2 283 284 nang = nang + 1 285 iang(1,nang) = ioff+2 286 iang(2,nang) = ioff+1 287 iang(3,nang) = ioff+3 288 ango(nang) = 104.52e0 289 ak(nang) = 100.00e0 290 291 q(ioff+1) = -0.8340e0 292 q(ioff+2) = 0.4170e0 293 q(ioff+3) = 0.4170e0 294 295 if (i.eq.1.or.j.eq.1.or.k.eq.1) then 296 iopt(ioff+1) = 1 297 iopt(ioff+2) = 1 298 iopt(ioff+3) = 1 299 else 300 iopt(ioff+1) = 0 301 iopt(ioff+2) = 0 302 iopt(ioff+3) = 0 303 endif 304 305 nac(ioff+1) = 0 306 nac(ioff+2) = 1 307 nac(ioff+3) = 1 308 iac(1,ioff+2) = ioff+3 309 iac(1,ioff+3) = ioff+2 310 nad(ioff+1) = 0 311 nad(ioff+2) = 0 312 nad(ioff+3) = 0 313 314 ioff = ioff + 3 315 316 endif 317 318 endif 319 320 end do 321 322 end do 323 end do 324 end do 325 326 do i=iatoms,ioff 327 iopt(i) = 1 328 end do 329 330 iatoms = ioff 331 332 if (nion.ne.0) 333 & call addions(coo,iconn,iresid,ityp,q,iopt) 334 335 return 336 end 337 338 logical function chkion(iwat) 339 implicit real (a-h,o-z), integer (i-n) 340 parameter (mxion=2000) 341 common /h2oer/ numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat 342 343 chkion = .false. 344 if (nion.eq.0) return 345 346 do i=1,nion 347 if (iwat.eq.iw(i)) then 348 chkion = .true. 349 return 350 endif 351 end do 352 353 return 354 end 355 356 subroutine addions(coo,iconn,iresid,ityp,q,iopt) 357 implicit real (a-h,o-z), integer (i-n) 358 parameter (mxcon=10) 359 parameter (numres=50000) 360 parameter (mxion=2000) 361 common /athlp/ iatoms, mxnat 362 common /residu/ qtot,ihsres, 363 & ires(numres),ibeg(numres),iend(numres) 364 common /h2oer/ numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat 365 integer*2 ityp 366 dimension coo(3,*), ityp(*),iconn(mxcon+1,*),iresid(*) 367 dimension q(*),iopt(*) 368 369c ADD ION 370 371 ioff = iatoms 372 ishoh = iresid(ioff) - 1 373 374 do i=1,nion 375 ihsres = ihsres + 1 376 ires(ihsres) = ishoh 377 ibeg(ihsres) = ioff+1 378 iend(ihsres) = ioff+1 379 iresid(ioff+1) = ishoh 380 381 ishoh = ishoh - 1 382 383 iatptr = natnow + (iw(i)-1)*3 + 1 384 385 do j=1,3 386 coo(j,ioff+1) = coo(j,iatptr) 387 end do 388 389 if (iontyp.eq.1) then 390 ityp(ioff+1) = 659 391 q(ioff+1) = -1.0e0 392 else 393 ityp(ioff+1) = 652 394 q(ioff+1) = +1.0e0 395 endif 396 397 iconn(1,ioff+1) = 0 398 iopt(ioff+1) = 1 399 400 ioff = ioff + 1 401 402 end do 403 404 iatoms = ioff 405 406 return 407 end 408 409 subroutine cntwad(water,coo,iconn,iresid,ityp,q,iopt,qpot) 410 implicit real (a-h,o-z), integer (i-n) 411 parameter (mxcon=10) 412 parameter (mxac=3*mxcon) 413 parameter (numres=50000) 414 parameter (mxion=2000) 415 logical doind,exclu,box,cell,fast,chkwat,chkbox 416 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 417 common /athlp/ iatoms, mxnat 418 common /residu/ qtot,ihsres, 419 & ires(numres),ibeg(numres),iend(numres) 420 common /h2oer/ numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat 421 integer*2 ityp 422 dimension water(3,*), coo(3,*), nbox(3), doff(3) 423 dimension ityp(*),iconn(mxcon+1,*),iresid(*) 424 dimension q(*),iopt(*),cwat(3,3),c1(3),c2(3),ibox(3),qpot(*) 425 426 watd = 18.620e0 427 nwater = 648/3 428 429c 18.620 ang box 430c take prefilled water box and tile the real box with this 431c sub box 432c remove water that overlap with protein 433 434 do i=1,3 435 t = abc(i)/watd 436 nbox(i) = int(t) + 1 437 end do 438 439 ioff = iatoms 440 nitmp = 0 441 442 ibox(1) = 0 443 444 do i=1,nbox(1) 445 if (i.eq.nbox(1)) ibox(1) = 1 446 ibox(2) = 0 447 448 do j=1,nbox(2) 449 if (j.eq.nbox(2)) ibox(2) = 1 450 ibox(3) = 0 451 452 do k=1,nbox(3) 453 if (k.eq.nbox(3)) ibox(3) = 1 454 455 doff(1) = -abc2(1) + 0.5e0*watd + watd*dble(i-1) 456 doff(2) = -abc2(2) + 0.5e0*watd + watd*dble(j-1) 457 doff(3) = -abc2(3) + 0.5e0*watd + watd*dble(k-1) 458 459 chkbox = .false. 460 if (i.eq.nbox(1).or.j.eq.nbox(2).or.k.eq.nbox(3)) 461 & chkbox = .true. 462 463 do n=1,nwater 464 465 in = (n-1)*3 466 467 do m=1,3 468 cwat(m,1) = water(m,in+1) + doff(m) 469 cwat(m,2) = water(m,in+2) + doff(m) 470 cwat(m,3) = water(m,in+3) + doff(m) 471 end do 472 473 if (chkwat(cwat,coo,ityp,iopt,ioff,ibox, 474 & iwt,chkbox)) then 475 476 numwat = numwat + 1 477 478 if (i.eq.1.or.j.eq.1.or.k.eq.1) then 479 iopt(ioff+1) = 1 480 iopt(ioff+2) = 1 481 iopt(ioff+3) = 1 482 else 483 iopt(ioff+1) = 0 484 iopt(ioff+2) = 0 485 iopt(ioff+3) = 0 486 endif 487 488 ityp(ioff+1) = 649 489 ityp(ioff+2) = 650 490 ityp(ioff+3) = 650 491 492 do m=1,3 493 coo(m,ioff+1) = cwat(m,1) 494 coo(m,ioff+2) = cwat(m,2) 495 coo(m,ioff+3) = cwat(m,3) 496 end do 497 498 if (ionpl.eq.1) then 499 call clmond(cwat(1,1),pot,coo,q) 500 qpot(numwat) = pot 501 endif 502 503 ioff = ioff + 3 504 505 endif 506 507 end do 508 509 end do 510 end do 511 end do 512 513 if (ionpl.eq.1) then 514 515 nitmp = 0 516 517 do while (nitmp.lt.nion) 518 519 doind = .false. 520 potl = 1.0e10 521 iptr = 0 522 523 do i=1,numwat 524 525 pot = qpot(i) 526 527 exclu = .false. 528 do j=1,nitmp 529 if (iw(j).eq.i) exclu = .true. 530 end do 531 532 if (pot.lt.potl.and..not.exclu) then 533 doind = .true. 534 potl = pot 535 iptr = i 536 endif 537 538 end do 539 540 if (doind) then 541 nitmp = nitmp + 1 542 iw(nitmp) = iptr 543 544 do i=1,numwat 545 iatptr = iatoms + (i-1)*3 + 1 546 do j=1,3 547 c1(j) = coo(j,iatptr) 548 end do 549 550 jatptr = iatoms + (iptr-1)*3 + 1 551 do j=1,3 552 c2(j) = coo(j,jatptr) 553 end do 554 r2 = dist2(c1,c2) 555 r = sqrt(r2) 556 qpot(i) = qpot(i) + 1.0e0/r 557 558 end do 559 560 endif 561 562 end do 563 564 else 565 566 call crionp 567 568 endif 569 570 return 571 end 572 573 logical function chkvec(vec,cwat,coo,iopt,ityp,ioff) 574 implicit real (a-h,p-w),integer (i-n),logical (o) 575 parameter (mxamb=1590) 576 parameter (mxgff=72) 577 integer ambcls 578 common /typpar/ ambchg(mxamb),ambcls(mxamb) 579 parameter (mxcls=49) 580 common /clspar/ ambvdwr(mxcls),ambvdwe(mxcls),mapagf(mxcls) 581 integer amb2gf 582 common /gftyp/ gfvdw(2,mxgff),amb2gf(mxamb) 583 common /athlp/ iatoms, mxnat 584 integer*2 ityp 585 dimension v(3),vec(3),cwat(3,3),coo(3,*),ityp(*) 586 dimension iopt(*) 587 588 chkvec = .true. 589 590 vdwat = 2.7683e0 591 592 do i=iatoms,ioff 593 594 if (iopt(i).eq.1) then 595 596 i1 = int(ityp(i)) 597 if (i1.gt.0) then 598 il = ambcls(i1) 599 vdwr = ambvdwr(il) 600 else 601 il = iabs(i1) 602 vdwr = gfvdw(1,il) 603 endif 604 605 606 do k=1,3 607 v(k) = coo(k,i) + vec(k) - cwat(k,1) 608 end do 609 610 rab2 = v(1)*v(1) + v(2)*v(2) + v(3)*v(3) 611 dmaxsq = vdwr + vdwat 612 dmaxsq = dmaxsq * dmaxsq 613 614 if (rab2.lt.dmaxsq) then 615 chkvec = .false. 616 return 617 endif 618 619 620 endif 621 622 end do 623 624 return 625 end 626 627 logical function chkwat(cwat,coo,ityp,iopt,ioff,ibox,iwrpr,chkbox) 628 implicit real (a-h,p-w),integer (i-n),logical (o) 629 parameter (mxamb=1590) 630 integer ambcls 631 common /typpar/ ambchg(mxamb),ambcls(mxamb) 632 parameter (mxcls=49) 633 common /clspar/ ambvdwr(mxcls),ambvdwe(mxcls),mapagf(mxcls) 634 parameter (mxgff=72) 635 integer amb2gf 636 common /gftyp/ gfvdw(2,mxgff),amb2gf(mxamb) 637 integer*2 ityp 638 logical chkbox,chkvec 639 common /athlp/ iatoms, mxnat 640 logical box,cell,fast,owat 641 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 642 dimension coo(3,*),ityp(*),cwat(3,3),v(3),vdwat(3),ibox(3) 643 dimension vec(3),iopt(*) 644 645 chkwat = .false. 646 647 648 owat = .true. 649 iwrpr = 0 650 651 if (chkbox) then 652 653 do i=1,3 654 do j=1,3 655 if (cwat(j,i).gt.abc2(j)) owat = .false. 656 end do 657 end do 658 659 if (.not.owat) return 660 661 do i=1,3 662 663 do j=1,3 664 vec(j) = 0.0e0 665 end do 666 667 if (ibox(i).eq.1) then 668 vec(i) = abc(i) 669 if (.not.chkvec(vec,cwat,coo,iopt,ityp,ioff)) return 670 endif 671 672 end do 673 674c chkwat = .true. 675c return 676 677 endif 678 679c check overlap with protein atoms 680 681 owat = .true. 682 683 vdwat(1) = 2.7683e0 684 685 do i=1,iatoms 686 687 i1 = int(ityp(i)) 688 if (i1.gt.0) then 689 il = ambcls(i1) 690 vdwr = ambvdwr(il) 691 elseif (i1.le.0) then 692 i1 = iabs(i1) 693 vdwr = gfvdw(1,i1) 694 endif 695 696 do k=1,3 697 v(k) = coo(k,i) - cwat(k,1) 698 end do 699 rab2 = v(1)*v(1) + v(2)*v(2) + v(3)*v(3) 700 dmaxsq = vdwr + vdwat(1) 701 dmaxsq = dmaxsq * dmaxsq 702 if (rab2.lt.dmaxsq) owat = .false. 703 if (rab2.lt.81.0e0) iwrpr = 1 704 705 end do 706 707 if (owat) chkwat = .true. 708 709 return 710 end 711 712 subroutine docnt(vecm,coo) 713 implicit real (a-h,p-w),integer (i-n),logical (o) 714 common /athlp/ iatoms, mxnat 715 dimension vec(3),vecm(3) 716 dimension coo(3,*) 717 718 do i=1,3 719 vecm(i) = 0.0e0 720 end do 721 722 call cntvec(vec,coo,iatoms) 723 724 do i=1,iatoms 725 do j=1,3 726 coo(j,i) = coo(j,i) - vec(j) 727 da = abs(coo(j,i)) 728 if (da.gt.vecm(j)) vecm(j) = da 729 end do 730 end do 731 732 return 733 end 734 735 subroutine cntvec(vec,coo,iatoms) 736 implicit real (a-h,p-w),integer (i-n),logical (o) 737 dimension vec(3), coo(3,*) 738 739 do i=1,3 740 vec(i) = 0.0e0 741 end do 742 743 if (iatoms.le.0) return 744 745 do i=1,iatoms 746 do j=1,3 747 vec(j) = vec(j) + coo(j,i) 748 end do 749 end do 750 751 do j=1,3 752 vec(j) = vec(j) / dble(iatoms) 753 end do 754 755 return 756 end 757 758 subroutine setop(a,b,c,alpha,beta,gamma) 759 implicit real (a-h,p-z),integer (i-n),logical (o) 760 common /cell/ xa,ya,yb,za,zb,zc 761 real rxa,rya,ryb,rza,rzb,rzc 762 common /cellr/ rxa,rya,ryb,rza,rzb,rzc 763 764 torad = atan(1.0e0) / 45.0e0 765 alpha = alpha*torad 766 beta = beta*torad 767 gamma = gamma*torad 768 769 ca = cos(alpha) 770 cb = cos(beta) 771 cc = cos(gamma) 772 sc = sin(gamma) 773 sa = sin(alpha) 774 cvol = a*b*c* 775 & sqrt(1-ca**2-cb**2-cc**2+2.0e0*ca*cb*cc) 776 xa = cvol / (b*c*sa) 777 ya = (a*(cc-(ca*cb)))/sa 778 yb = (b*sa) 779 za = (a*cb) 780 zb = (b*ca) 781 zc = c 782 783 rxa = real(xa) 784 rya = real(ya) 785 ryb = real(yb) 786 rza = real(za) 787 rzb = real(zb) 788 rzc = real(zc) 789 790c print*,' ' 791c write(*,'(a,f9.3)') ' Cell Volume = ',cvol 792 793 return 794 end 795 796 subroutine getcel(a,b,c,alpha,beta,gamma) 797 implicit real (a-h,p-z),integer (i-n),logical (o) 798 common /cell/ xa,ya,yb,za,zb,zc 799 800 torad = atan(1.0e0) / 45.0e0 801 802 a = sqrt(xa**2+ya**2+za**2) 803 b = sqrt(yb**2+zb**2) 804 c = zc 805 806 sa = yb/b 807 ca = zb/b 808 cb = za/a 809 cc = ca*cb + ya*sa/a 810 811 alpha = acos(ca) / torad 812 beta = acos(cb) / torad 813 gamma = acos(cc) / torad 814 815 return 816 end 817 818 subroutine fr2crt(vec) 819 implicit real (a-h,p-z),integer (i-n),logical (o) 820 common /cell/ xa,ya,yb,za,zb,zc 821 dimension vec(3) 822 823 vec(3) = za*vec(1) + zb*vec(2) + zc*vec(3) 824 vec(2) = ya*vec(1) + yb*vec(2) 825 vec(1) = xa*vec(1) 826 827 return 828 end 829 830 subroutine crt2fr(veci,veco) 831 implicit real (a-h,p-z),integer (i-n),logical (o) 832 common /cell/ xa,ya,yb,za,zb,zc 833 dimension veci(3),veco(3) 834 835 veco(1) = veci(1) / xa 836 veco(2) = (veci(2) - ya*veco(1)) / yb 837 veco(3) = (veci(3) - za*veco(1) - zb*veco(2)) / zc 838 839 return 840 end 841 842 subroutine fr2crr(vec) 843 implicit real (a-h,p-z),integer (i-n),logical (o) 844 common /cellr/ xa,ya,yb,za,zb,zc 845 dimension vec(3) 846 847 vec(3) = za*vec(1) + zb*vec(2) + zc*vec(3) 848 vec(2) = ya*vec(1) + yb*vec(2) 849 vec(1) = xa*vec(1) 850 851 return 852 end 853 854 subroutine crr2fr(veci,veco) 855 implicit real (a-h,p-z),integer (i-n),logical (o) 856 common /cellr/ xa,ya,yb,za,zb,zc 857 dimension veci(3),veco(3) 858 859 veco(1) = veci(1) / xa 860 veco(2) = (veci(2) - ya*veco(1)) / yb 861 veco(3) = (veci(3) - za*veco(1) - zb*veco(2)) / zc 862 863 return 864 end 865 866 subroutine redfr(xyz) 867 implicit real (a-h,p-z),integer (i-n),logical (o) 868 logical box,cell,fast 869 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 870 dimension xyz(3),fr(3) 871 872 call crt2fr(xyz,fr) 873 874 do i=1,3 875 frt = 1.0e0 876 if (fr(i).lt.0.0e0) then 877 fr(i) = abs(fr(i)) 878 frt = -1.0e0 879 endif 880 881 fr(i) = mod(fr(i),1.0e0) 882 if (fr(i) .gt. 0.5e0) fr(i) = fr(i) - 1.0e0 883 884 xyz(i) = frt*fr(i) 885 end do 886 887 call fr2crt(xyz) 888 889 return 890 end 891 892 subroutine rddfrr(xyz) 893 implicit real (a-h,p-z),integer (i-n),logical (o) 894 real abc,abc2,angles 895 logical box,cell,fast 896 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 897 common /debug/ ideb 898 dimension xyz(3),fr(3) 899 900 901 call crr2fr(xyz,fr) 902 903 904 do i=1,3 905 frt = 1.0e0 906 if (fr(i).lt.0.0e0) then 907 fr(i) = abs(fr(i)) 908 frt = -1.0e0 909 endif 910 911 fr(i) = amod(fr(i),1.0e0) 912 if (fr(i) .gt. 0.5e0) fr(i) = fr(i) - 1.0e0 913 914 xyz(i) = frt*fr(i) 915 end do 916 917 918 call fr2crr(xyz) 919 920 return 921 end 922 923 subroutine adcldr(cellder,ded,de,fct,lex) 924 implicit real (a-h,p-z),integer (i-n),logical (o) 925 common /cell/ xa,ya,yb,za,zb,zc 926 common /expnd/ addfr(3,125),ie(5),izero 927 dimension ded(3),cellder(6),addf(3) 928 929 do j=1,3 930 addf(j) = addfr(j,lex) 931 end do 932 933 cellder(1) = cellder(1) + addf(1)*ded(1)*de*fct 934 935 cellder(2) = cellder(2) + addf(1)*ded(2)*de*fct 936 cellder(3) = cellder(3) + addf(2)*ded(2)*de*fct 937 938 cellder(4) = cellder(4) + addf(1)*ded(3)*de*fct 939 cellder(5) = cellder(5) + addf(2)*ded(3)*de*fct 940 cellder(6) = cellder(6) + addf(3)*ded(3)*de*fct 941 942 return 943 end 944 945 subroutine prcldr(cellder,par) 946 implicit real (a-h,p-z),integer (i-n),logical (o) 947 common /athlp/ iatoms, mxnat 948 dimension cellder(6),par(3,*) 949 950c do i=1,6 951c print*,'cellder (',i,')=',cellder(i) 952c end do 953 954 print*,'val ',par(1,iatoms+1),' cellder (1)=',cellder(1) 955 print*,'val ',par(2,iatoms+1),' cellder (2)=',cellder(2) 956 print*,'val ',par(3,iatoms+1),' cellder (3)=',cellder(3) 957 print*,'val ',par(1,iatoms+2),' cellder (4)=',cellder(4) 958 print*,'val ',par(2,iatoms+2),' cellder (5)=',cellder(5) 959 print*,'val ',par(3,iatoms+2),' cellder (6)=',cellder(6) 960 961 return 962 end 963 964 subroutine expfr(in,cxyz,exyz,lex,nexpnd,oje) 965 implicit real (a-h,p-z),integer (i-n),logical (o) 966 common /expnd/ addfr(3,125),ie(5),izero 967 logical oje(125) 968 dimension lex(125),cxyz(3,125),exyz(3,125,*) 969 970 l = 0 971 972 do i=1,125 973 974 if (i.ne.izero) then 975 l = l + 1 976 do j=1,3 977 cxyz(j,l) = exyz(j,i,in) - exyz(j,izero,in) 978 end do 979 lex(l) = i 980 oje(l) = .false. 981 endif 982 983 end do 984 985 return 986 end 987 988 subroutine setfrac(coo,frac,exyz) 989 implicit real (a-h,p-z),integer (i-n),logical (o) 990 common /athlp/ iatoms, mxnat 991 common /expnd/ addfr(3,125),ie(5),izero 992 dimension coo(3,*),frac(3,*),exyz(3,125,*),cxyz(3,125),dfr(3) 993 994 do j=1,3 995 dfr(j) = 0.0e0 996 end do 997 998 do i=1,iatoms 999 1000 do j=1,3 1001 if (abs(coo(j,i)).lt.0.001e0) coo(j,i) = 0.0e0 1002 end do 1003 1004 call crt2fr(coo(1,i),frac(1,i)) 1005 1006 do j=1,3 1007 if (frac(j,i).lt.dfr(j)) dfr(j) = frac(j,i) 1008 end do 1009 1010 end do 1011 1012c huh ??? 1013 1014 do i=1,5 1015 do j=1,3 1016 if ((dfr(j).ge.dble(-1*i)).and.(dfr(j).lt.dble(-1*(i-1)))) 1017 & dfr(j) = dble(-1*i) 1018 end do 1019 end do 1020 1021 do i=1,iatoms 1022 1023 do j=1,3 1024 frac(j,i) = frac(j,i) - dfr(j) 1025 end do 1026 1027 do k=1,125 1028 do j=1,3 1029 exyz(j,k,i) = frac(j,i) + addfr(j,k) 1030 end do 1031 call fr2crt(exyz(1,k,i)) 1032 end do 1033 1034 end do 1035 1036 return 1037 end 1038 1039 subroutine xpfr(in,kn,cxyz,exyz,lex,nexpnd,oje) 1040 implicit real (a-h,p-z),integer (i-n),logical (o) 1041 common /expnd/ addfr(3,125),ie(5),izero 1042 logical oje(125) 1043 dimension cxyz(3,125),exyz(3,125,*),lex(125) 1044 1045 do i=1,nexpnd+1 1046 oje(i) = .false. 1047 if (i.eq.izero) oje(i) = .true. 1048 do j=1,3 1049 cxyz(j,i) = exyz(j,i,in) - exyz(j,izero,kn) 1050 end do 1051 lex(i) = i 1052 end do 1053 1054 return 1055 end 1056 1057 subroutine initxp 1058 implicit real (a-h,p-z),integer (i-n),logical (o) 1059 common /expnd/ addfr(3,125),ie(5),izero 1060 data (ie(j),j=1,5) / -2,-1,0,1,2/ 1061 1062 izero = 0 1063 1064 l = 0 1065 do i=1,5 1066 do j=1,5 1067 do k=1,5 1068 l = l + 1 1069 addfr(1,l) = dble(ie(i)) 1070 addfr(2,l) = dble(ie(j)) 1071 addfr(3,l) = dble(ie(k)) 1072 if (ie(i).eq.0.and.ie(j).eq.0.and.ie(k).eq.0) izero = l 1073 end do 1074 end do 1075 end do 1076 1077 return 1078 end 1079 1080 subroutine var2cl(x,n) 1081c copy cell variables to cell common 1082 implicit real (a-h,p-z),integer (i-n),logical (o) 1083 logical box,cell,fast 1084 common /pbc/ abc(3),abc2(3),angles(3),box,cell,fast 1085 common /cell/ xa,ya,yb,za,zb,zc 1086 real rxa,rya,ryb,rza,rzb,rzc 1087 common /cellr/ rxa,rya,ryb,rza,rzb,rzc 1088 dimension x(*) 1089 1090 n3 = n*3 1091 1092 xa = x(n3+1) 1093 ya = x(n3+2) 1094 yb = x(n3+3) 1095 za = x(n3+4) 1096 zb = x(n3+5) 1097 zc = x(n3+6) 1098 1099 abc(1) = sqrt(xa**2+ya**2+za**2) 1100 abc(2) = sqrt(yb**2+zb**2) 1101 abc(3) = zc 1102 1103 do i=1,3 1104 abc2(i) = abc(i)*0.5e0 1105 end do 1106 1107 rxa = real(x(n3+1)) 1108 rya = real(x(n3+2)) 1109 ryb = real(x(n3+3)) 1110 rza = real(x(n3+4)) 1111 rzb = real(x(n3+5)) 1112 rzc = real(x(n3+6)) 1113 1114 return 1115 end 1116 1117 subroutine grd2var(cellder,f,n) 1118c copy cell variables to cell common 1119 implicit real (a-h,p-z),integer (i-n),logical (o) 1120 common /cell/ xa,ya,yb,za,zb,zc 1121 dimension cellder(6),f(*) 1122 1123c print*,'cellval ',xa,ya,yb,za,zb,zc 1124 1125 n3 = n*3 1126 1127 do i=1,6 1128 f(n3+i) = -cellder(i) 1129 end do 1130 1131 return 1132 end 1133 1134 subroutine uinner(uin,cellder,coo,f,q,n) 1135 implicit real (a-h,p-z),integer (i-n),logical (o) 1136 common /cell/ xa,ya,yb,za,zb,zc 1137 dimension cellder(6),pole(3),coo(3,*),f(3,*),q(*) 1138 1139 pi = 3.141592654e0 1140 pi43 = 4.0e0*pi/3.0e0 1141 econv = 332.05382e0 1142 fin = pi43*econv/(xa*yb*zc) 1143 1144 do j=1,3 1145 pole(j) = 0.0e0 1146 end do 1147 1148 do i=1,n 1149 do j=1,3 1150 pole(j) = pole(j) + coo(j,i)*q(i) 1151 end do 1152 end do 1153 1154 poll = 0.0e0 1155 do j=1,3 1156 poll = poll + pole(j)*pole(j) 1157 end do 1158 1159 uin = -0.5e0*fin*poll 1160 1161 cellder(1) = -uin/xa 1162 cellder(3) = -uin/yb 1163 cellder(6) = -uin/zc 1164 1165 do i=1,n 1166 do j=1,3 1167 f(j,i) = f(j,i) - fin*pole(j)*q(i) 1168 end do 1169 end do 1170 1171 return 1172 end 1173 1174 subroutine crionp 1175 parameter (mxion=2000) 1176 implicit real (a-h,p-z),integer (i-n),logical (o) 1177 common /h2oer/ numwat,natnow,nion,iontyp,ionpl,iw(mxion),niwat 1178 logical doit 1179 1180 nitmp = 0 1181 1182 do while (nitmp.lt.nion) 1183 1184 doit = .true. 1185 irnd = numwat*random() 1186 1187 do i=1,nitmp 1188 if (iw(i).eq.irnd) doit = .false. 1189 end do 1190 1191 if (doit) then 1192 nitmp = nitmp + 1 1193 iw(nitmp) = irnd 1194 endif 1195 1196 end do 1197 1198 return 1199 end 1200 1201 real function random() 1202 parameter (mplier=16807,modlus=2147483647,mobymp=127773, 1203 & momdmp=2836) 1204 common /seed/jseed,ifrst,nextn 1205c mseed comes from alloc.c 1206 integer hvlue, lvlue, testv, nextn, mseed 1207 1208 if (ifrst.eq.0) then 1209 jseed = mseed() 1210 nextn = jseed 1211 ifrst = 1 1212 endif 1213 1214 hvlue = nextn / mobymp 1215 lvlue = mod(nextn, mobymp) 1216 testv = mplier*lvlue - momdmp*hvlue 1217 1218 if (testv.gt.0) then 1219 nextn = testv 1220 else 1221 nextn = testv + modlus 1222 endif 1223 1224 random = dble(nextn)/dble(modlus) 1225 1226 return 1227 end 1228 1229