1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1990 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ############################################################# 9c ## ## 10c ## program xyzedit -- editing of Cartesian coordinates ## 11c ## ## 12c ############################################################# 13c 14c 15c "xyzedit" provides for modification and manipulation 16c of the contents of Cartesian coordinates files 17c 18c 19 program xyzedit 20 use atomid 21 use atoms 22 use bound 23 use boxes 24 use couple 25 use files 26 use inform 27 use iounit 28 use limits 29 use math 30 use molcul 31 use ptable 32 use titles 33 use units 34 use usage 35 implicit none 36 integer i,j,k,it 37 integer ixyz,imod 38 integer init,stop 39 integer nmode,mode 40 integer natom,atmnum 41 integer nlist 42 integer offset,origin 43 integer oldtype,newtype 44 integer freeunit 45 integer trimtext 46 integer, allocatable :: list(:) 47 integer, allocatable :: keep(:) 48 real*8 xi,yi,zi 49 real*8 xr,yr,zr 50 real*8 xcm,ycm,zcm 51 real*8 xnew,ynew,znew 52 real*8 xorig,yorig,zorig 53 real*8 ri,rij,dij 54 real*8 phi,theta,psi 55 real*8 cphi,ctheta,cpsi 56 real*8 sphi,stheta,spsi 57 real*8 dist2,cut2 58 real*8 random,norm,weigh 59 real*8, allocatable :: rad(:) 60 real*8 a(3,3) 61 logical exist,query 62 logical opened,multi 63 logical append 64 character*3 symb 65 character*240 xyzfile 66 character*240 modfile 67 character*240 record 68 character*240 string 69 external random,merge 70c 71c 72c initialize various constants and the output flags 73c 74 call initial 75 opened = .false. 76 multi = .false. 77 nmode = 24 78 offset = 0 79c 80c try to get a filename from the command line arguments 81c 82 call nextarg (xyzfile,exist) 83 if (exist) then 84 call basefile (xyzfile) 85 call suffix (xyzfile,'xyz','old') 86 inquire (file=xyzfile,exist=exist) 87 end if 88c 89c ask for the user specified input structure filename 90c 91 do while (.not. exist) 92 write (iout,10) 93 10 format (/,' Enter Cartesian Coordinate File Name : ',$) 94 read (input,20) xyzfile 95 20 format (a240) 96 call basefile (xyzfile) 97 call suffix (xyzfile,'xyz','old') 98 inquire (file=xyzfile,exist=exist) 99 end do 100c 101c open and then read the Cartesian coordinates file 102c 103 ixyz = freeunit () 104 open (unit=ixyz,file=xyzfile,status='old') 105 rewind (unit=ixyz) 106 call readxyz (ixyz) 107c 108c get the force field definition and assign atom types 109c 110 call attach 111 call field 112 call katom 113c 114c present a list of possible coordinate modifications 115c 116 write (iout,30) 117 30 format (/,' The Tinker XYZ File Editing Utility Can :', 118 & //,4x,'(1) Offset the Numbers of the Current Atoms', 119 & /,4x,'(2) Remove User Specified Individual Atoms', 120 & /,4x,'(3) Remove User Specified Types of Atoms', 121 & /,4x,'(4) Delete Inactive Atoms Beyond Cutoff Range', 122 & /,4x,'(5) Insertion of Individual Specified Atoms', 123 & /,4x,'(6) Replace Old Atom Type with a New Type', 124 & /,4x,'(7) Assign Connectivities for Linear Chain', 125 & /,4x,'(8) Assign Connectivities Based on Distance', 126 & /,4x,'(9) Assign Atom Types for BASIC Force Field', 127 & /,3x,'(10) Convert Units from Bohrs to Angstroms', 128 & /,3x,'(11) Invert thru Origin to Give Mirror Image', 129 & /,3x,'(12) Translate All Atoms by an X,Y,Z-Vector', 130 & /,3x,'(13) Translate Center of Mass to the Origin', 131 & /,3x,'(14) Translate a Specified Atom to the Origin', 132 & /,3x,'(15) Translate and Rotate to Inertial Frame', 133 & /,3x,'(16) Move to Specified Rigid Body Coordinates', 134 & /,3x,'(17) Move Stray Molecules into Periodic Box', 135 & /,3x,'(18) Trim a Periodic Box to a Smaller Size', 136 & /,3x,'(19) Make Truncated Octahedron from Cubic Box', 137 & /,3x,'(20) Make Rhombic Dodecahedron from Cubic Box', 138 & /,3x,'(21) Append a Second XYZ File to Current One', 139 & /,3x,'(22) Create and Fill a Periodic Boundary Box', 140 & /,3x,'(23) Soak Current Molecule in Box of Solvent', 141 & /,3x,'(24) Place Monoatomic Ions around a Solute') 142c 143c get the desired type of coordinate file modification 144c 145 40 continue 146 abort = .false. 147 mode = -1 148 query = .true. 149 call nextarg (string,exist) 150 if (exist) then 151 read (string,*,err=50,end=50) mode 152 if (mode.ge.0 .and. mode.le.nmode) query = .false. 153 end if 154 50 continue 155 if (query) then 156 do while (mode.lt.0 .or. mode.gt.nmode) 157 mode = 0 158 write (iout,60) 159 60 format (/,' Number of the Desired Choice [<Enter>=Exit]', 160 & ' : ',$) 161 read (input,70,err=40,end=80) mode 162 70 format (i10) 163 80 continue 164 end do 165 end if 166c 167c open the file to be used for the output coordinates 168c 169 if (mode.gt.0 .and. .not.opened) then 170 opened = .true. 171 imod = freeunit () 172 modfile = filename(1:leng)//'.xyz' 173 call version (modfile,'new') 174 open (unit=imod,file=modfile,status='new') 175 end if 176c 177c get the offset value to be used in atom renumbering 178c 179 if (mode .eq. 1) then 180 90 continue 181 offset = 0 182 query = .true. 183 call nextarg (string,exist) 184 if (exist) then 185 read (string,*,err=100,end=100) offset 186 query = .false. 187 end if 188 100 continue 189 if (query) then 190 write (iout,110) 191 110 format (/,' Offset used to Renumber the Atoms [0] : ',$) 192 read (input,120,err=90) offset 193 120 format (i10) 194 end if 195 do while (.not. abort) 196 call makeref (1) 197 call readxyz (ixyz) 198 if (.not. abort) multi = .true. 199 if (multi) then 200 call makeref (2) 201 call getref (1) 202 call prtmod (imod,offset) 203 call getref (2) 204 end if 205 end do 206 if (.not. multi) then 207 call getref (1) 208 goto 40 209 end if 210 end if 211c 212c remove a specified list of individual atoms 213c 214 if (mode .eq. 2) then 215 allocate (list(n)) 216 nlist = 0 217 do i = 1, n 218 list(i) = 0 219 end do 220 write (iout,130) 221 130 format (/,' Numbers of the Atoms to be Removed : ',$) 222 read (input,140) record 223 140 format (a240) 224 read (record,*,err=150,end=150) (list(i),i=1,n) 225 150 continue 226 do while (list(nlist+1) .ne. 0) 227 nlist = nlist + 1 228 end do 229 do i = 1, nlist 230 if (list(i) .gt. n) list(i) = n 231 if (list(i) .lt. -n) list(i) = -n 232 end do 233 call sort4 (nlist,list) 234 do while (.not. abort) 235 do i = nlist, 1, -1 236 if (i .gt. 1) then 237 if (list(i-1) .lt. 0) then 238 do j = abs(list(i)), abs(list(i-1)), -1 239 call delete (j) 240 end do 241 else if (list(i) .gt. 0) then 242 call delete (list(i)) 243 end if 244 else if (list(i) .gt. 0) then 245 call delete (list(i)) 246 end if 247 end do 248 call makeref (1) 249 call readxyz (ixyz) 250 if (.not. abort) multi = .true. 251 if (multi) then 252 call makeref (2) 253 call getref (1) 254 call prtmod (imod,offset) 255 call getref (2) 256 end if 257 end do 258 deallocate (list) 259 if (.not. multi) then 260 call getref (1) 261 goto 40 262 end if 263 end if 264c 265c remove atoms with any of a specified list of atom types 266c 267 if (mode .eq. 3) then 268 allocate (list(n)) 269 nlist = 0 270 do i = 1, n 271 list(i) = 0 272 end do 273 write (iout,160) 274 160 format (/,' Atom Types to be Removed : ',$) 275 read (input,170) record 276 170 format (a240) 277 read (record,*,err=180,end=180) (list(i),i=1,n) 278 180 continue 279 do while (list(nlist+1) .ne. 0) 280 nlist = nlist + 1 281 end do 282 natom = n 283 do while (.not. abort) 284 do i = natom, 1, -1 285 it = type(i) 286 do j = 1, nlist 287 if (list(j) .eq. it) then 288 call delete (i) 289 goto 190 290 end if 291 end do 292 190 continue 293 end do 294 call makeref (1) 295 call readxyz (ixyz) 296 if (.not. abort) multi = .true. 297 if (multi) then 298 call makeref (2) 299 call getref (1) 300 call prtmod (imod,offset) 301 call getref (2) 302 end if 303 end do 304 deallocate (list) 305 if (.not. multi) then 306 call getref (1) 307 goto 40 308 end if 309 end if 310c 311c remove atoms that are inactive and lie outside all cutoffs 312c 313 if (mode .eq. 4) then 314 call active 315 call cutoffs 316 cut2 = 0.0d0 317 if (vdwcut .le. 1000.0d0) cut2 = max(vdwcut**2,cut2) 318 if (chgcut .le. 1000.0d0) cut2 = max(chgcut**2,cut2) 319 if (dplcut .le. 1000.0d0) cut2 = max(dplcut**2,cut2) 320 if (mpolecut .le. 1000.0d0) cut2 = max(mpolecut**2,cut2) 321 if (cut2 .eq. 0.0d0) cut2 = 1.0d16 322 allocate (list(n)) 323 allocate (keep(n)) 324 do while (.not. abort) 325 nlist = 0 326 do i = 1, n 327 keep(i) = 0 328 end do 329 do i = 1, n 330 if (.not. use(i)) then 331 do j = 1, n12(i) 332 keep(i12(j,i)) = i 333 end do 334 do j = 1, n13(i) 335 keep(i13(j,i)) = i 336 end do 337 do j = 1, n14(i) 338 keep(i14(j,i)) = i 339 end do 340 xi = x(i) 341 yi = y(i) 342 zi = z(i) 343 do j = 1, n 344 if (use(j)) then 345 if (keep(j) .eq. i) goto 200 346 dist2 = (x(j)-xi)**2+(y(j)-yi)**2+(z(j)-zi)**2 347 if (dist2 .le. cut2) goto 200 348 end if 349 end do 350 nlist = nlist + 1 351 list(nlist) = i 352 200 continue 353 end if 354 end do 355 do i = nlist, 1, -1 356 call delete (list(i)) 357 end do 358 call makeref (1) 359 call readxyz (ixyz) 360 if (.not. abort) multi = .true. 361 if (multi) then 362 call makeref (2) 363 call getref (1) 364 call prtmod (imod,offset) 365 call getref (2) 366 end if 367 end do 368 deallocate (list) 369 deallocate (keep) 370 if (.not. multi) then 371 call getref (1) 372 goto 40 373 end if 374 end if 375c 376c insert a specified list of individual atoms 377c 378 if (mode .eq. 5) then 379 allocate (list(n)) 380 nlist = 0 381 do i = 1, n 382 list(i) = 0 383 end do 384 write (iout,210) 385 210 format (/,' Numbers of the Atoms to be Inserted : ',$) 386 read (input,220) record 387 220 format (a240) 388 read (record,*,err=230,end=230) (list(i),i=1,n) 389 230 continue 390 do while (list(nlist+1) .ne. 0) 391 nlist = nlist + 1 392 end do 393 call sort4 (nlist,list) 394 do while (.not. abort) 395 do i = nlist, 1, -1 396 if (i .gt. 1) then 397 if (list(i-1) .lt. 0) then 398 do j = abs(list(i-1)), abs(list(i)) 399 call insert (j) 400 end do 401 else if (list(i) .gt. 0) then 402 call insert (list(i)) 403 end if 404 else if (list(i) .gt. 0) then 405 call insert (list(i)) 406 end if 407 end do 408 call makeref (1) 409 call readxyz (ixyz) 410 if (.not. abort) multi = .true. 411 if (multi) then 412 call makeref (2) 413 call getref (1) 414 call prtmod (imod,offset) 415 call getref (2) 416 end if 417 end do 418 deallocate (list) 419 if (.not. multi) then 420 call getref (1) 421 goto 40 422 end if 423 end if 424c 425c get an old atom type and new atom type for replacement 426c 427 if (mode .eq. 6) then 428 240 continue 429 oldtype = 0 430 newtype = 0 431 call nextarg (string,exist) 432 if (exist) read (string,*,err=250,end=250) oldtype 433 call nextarg (string,exist) 434 if (exist) read (string,*,err=250,end=250) newtype 435 250 continue 436 if (oldtype.eq.0 .or. newtype.eq.0) then 437 write (iout,260) 438 260 format (/,' Numbers of the Old and New Atom Types : ',$) 439 read (input,270) record 440 270 format (a240) 441 end if 442 read (record,*,err=240,end=240) oldtype,newtype 443 do while (.not. abort) 444 do i = 1, n 445 if (type(i) .eq. oldtype) then 446 type(i) = newtype 447 end if 448 end do 449 call katom 450 call makeref (1) 451 call readxyz (ixyz) 452 if (.not. abort) multi = .true. 453 if (multi) then 454 call makeref (2) 455 call getref (1) 456 call prtmod (imod,offset) 457 call getref (2) 458 end if 459 end do 460 if (.not. multi) then 461 call getref (1) 462 goto 40 463 end if 464 end if 465c 466c assign atom connectivities to produce a linear chain 467c 468 if (mode .eq. 7) then 469 do while (.not. abort) 470 do i = 1, n 471 n12(i) = 0 472 if (i .ne. 1) then 473 n12(i) = n12(i) + 1 474 i12(n12(i),i) = i - 1 475 end if 476 if (i .ne. n) then 477 n12(i) = n12(i) + 1 478 i12(n12(i),i) = i + 1 479 end if 480 end do 481 call makeref (1) 482 call readxyz (ixyz) 483 if (.not. abort) multi = .true. 484 if (multi) then 485 call makeref (2) 486 call getref (1) 487 call prtmod (imod,offset) 488 call getref (2) 489 end if 490 end do 491 if (.not. multi) then 492 call getref (1) 493 goto 40 494 end if 495 end if 496c 497c assign atom connectivities based on interatomic distances 498c 499 if (mode .eq. 8) then 500 allocate (rad(n)) 501 do while (.not. abort) 502 call unitcell 503 call lattice 504 do i = 1, n 505 rad(i) = 0.76d0 506 atmnum = atomic(i) 507 if (atmnum .ne. 0) rad(i) = covrad(atmnum) 508 rad(i) = 1.15d0 * rad(i) 509 end do 510 do i = 1, n 511 n12(i) = 0 512 end do 513 do i = 1, n-1 514 xi = x(i) 515 yi = y(i) 516 zi = z(i) 517 ri = rad(i) 518 do j = i+1, n 519 xr = x(j) - xi 520 yr = y(j) - yi 521 zr = z(j) - zi 522 rij = ri + rad(j) 523 dij = sqrt(xr*xr + yr*yr + zr*zr) 524 if (dij .lt. rij) then 525 n12(i) = n12(i) + 1 526 i12(n12(i),i) = j 527 n12(j) = n12(j) + 1 528 i12(n12(j),j) = i 529 end if 530 end do 531 end do 532 do i = 1, n 533 call sort (n12(i),i12(1,i)) 534 end do 535 call makeref (1) 536 call readxyz (ixyz) 537 if (.not. abort) multi = .true. 538 if (multi) then 539 call makeref (2) 540 call getref (1) 541 call prtmod (imod,offset) 542 call getref (2) 543 end if 544 end do 545 deallocate (rad) 546 if (.not. multi) then 547 call getref (1) 548 goto 40 549 end if 550 end if 551c 552c assign atom types for the Tinker BASIC force field 553c 554 if (mode .eq. 9) then 555 do while (.not. abort) 556 do i = 1, n 557 k = atomic(i) 558 if (k .eq. 0) then 559 symb = name(i) 560 call upcase (symb(1:1)) 561 call lowcase (symb(2:3)) 562 do j = 1, maxele 563 if (symb .eq. elemnt(j)) then 564 k = j 565 goto 280 566 end if 567 end do 568 280 continue 569 end if 570 type(i) = 10*k + n12(i) 571 end do 572 call makeref (1) 573 call readxyz (ixyz) 574 if (.not. abort) multi = .true. 575 if (multi) then 576 call makeref (2) 577 call getref (1) 578 call prtmod (imod,offset) 579 call getref (2) 580 end if 581 end do 582 if (.not. multi) then 583 call getref (1) 584 goto 40 585 end if 586 end if 587c 588c convert the coordinate units from Bohrs to Angstroms 589c 590 if (mode .eq. 10) then 591 do while (.not. abort) 592 do i = 1, n 593 x(i) = x(i) * bohr 594 y(i) = y(i) * bohr 595 z(i) = z(i) * bohr 596 end do 597 call makeref (1) 598 call readxyz (ixyz) 599 if (.not. abort) multi = .true. 600 if (multi) then 601 call makeref (2) 602 call getref (1) 603 call prtmod (imod,offset) 604 call getref (2) 605 end if 606 end do 607 if (.not. multi) then 608 call getref (1) 609 goto 40 610 end if 611 end if 612c 613c get mirror image by inverting coordinates through origin 614c 615 if (mode .eq. 11) then 616 do while (.not. abort) 617 do i = 1, n 618 x(i) = -x(i) 619 y(i) = -y(i) 620 z(i) = -z(i) 621 end do 622 call makeref (1) 623 call readxyz (ixyz) 624 if (.not. abort) multi = .true. 625 if (multi) then 626 call makeref (2) 627 call getref (1) 628 call prtmod (imod,offset) 629 call getref (2) 630 end if 631 end do 632 if (.not. multi) then 633 call getref (1) 634 goto 40 635 end if 636 end if 637c 638c translate the entire system by a specified x,y,z-vector 639c 640 if (mode .eq. 12) then 641 xr = 0.0d0 642 yr = 0.0d0 643 zr = 0.0d0 644 call nextarg (string,exist) 645 if (exist) read (string,*,err=290,end=290) xr 646 call nextarg (string,exist) 647 if (exist) read (string,*,err=290,end=290) yr 648 call nextarg (string,exist) 649 if (exist) read (string,*,err=290,end=290) zr 650 290 continue 651 if (xr.eq.0.0d0 .and. yr.eq.0.0d0 .and. zr.eq.0.0d0) then 652 write (iout,300) 653 300 format (/,' Enter Translation Vector Components : ',$) 654 read (input,310) record 655 310 format (a240) 656 read (record,*,err=320,end=320) xr,yr,zr 657 320 continue 658 end if 659 do while (.not. abort) 660 do i = 1, n 661 x(i) = x(i) + xr 662 y(i) = y(i) + yr 663 z(i) = z(i) + zr 664 end do 665 call makeref (1) 666 call readxyz (ixyz) 667 if (.not. abort) multi = .true. 668 if (multi) then 669 call makeref (2) 670 call getref (1) 671 call prtmod (imod,offset) 672 call getref (2) 673 end if 674 end do 675 if (.not. multi) then 676 call getref (1) 677 goto 40 678 end if 679 end if 680c 681c translate the center of mass to the coordinate origin 682c 683 if (mode .eq. 13) then 684 do while (.not. abort) 685 xcm = 0.0d0 686 ycm = 0.0d0 687 zcm = 0.0d0 688 norm = 0.0d0 689 do i = 1, n 690 weigh = mass(i) 691 xcm = xcm + x(i)*weigh 692 ycm = ycm + y(i)*weigh 693 zcm = zcm + z(i)*weigh 694 norm = norm + weigh 695 end do 696 xcm = xcm / norm 697 ycm = ycm / norm 698 zcm = zcm / norm 699 do i = 1, n 700 x(i) = x(i) - xcm 701 y(i) = y(i) - ycm 702 z(i) = z(i) - zcm 703 end do 704 call makeref (1) 705 call readxyz (ixyz) 706 if (.not. abort) multi = .true. 707 if (multi) then 708 call makeref (2) 709 call getref (1) 710 call prtmod (imod,offset) 711 call getref (2) 712 end if 713 end do 714 if (.not. multi) then 715 call getref (1) 716 goto 40 717 end if 718 end if 719c 720c translate to place a specified atom at the origin 721c 722 if (mode .eq. 14) then 723 origin = 0 724 call nextarg (string,exist) 725 if (exist) read (string,*,err=330,end=330) origin 726 330 continue 727 if (origin .eq. 0) then 728 write (iout,340) 729 340 format (/,' Number of the Atom to Move to the Origin', 730 & ' : ',$) 731 read (input,350) origin 732 350 format (i10) 733 end if 734 do while (.not. abort) 735 xorig = x(origin) 736 yorig = y(origin) 737 zorig = z(origin) 738 do i = 1, n 739 x(i) = x(i) - xorig 740 y(i) = y(i) - yorig 741 z(i) = z(i) - zorig 742 end do 743 call makeref (1) 744 call readxyz (ixyz) 745 if (.not. abort) multi = .true. 746 if (multi) then 747 call makeref (2) 748 call getref (1) 749 call prtmod (imod,offset) 750 call getref (2) 751 end if 752 end do 753 if (.not. multi) then 754 call getref (1) 755 goto 40 756 end if 757 end if 758c 759c translate and rotate into standard orientation 760c 761 if (mode .eq. 15) then 762 do while (.not. abort) 763 call inertia (2) 764 call makeref (1) 765 call readxyz (ixyz) 766 if (.not. abort) multi = .true. 767 if (multi) then 768 call makeref (2) 769 call getref (1) 770 call prtmod (imod,offset) 771 call getref (2) 772 end if 773 end do 774 if (.not. multi) then 775 call getref (1) 776 goto 40 777 end if 778 end if 779c 780c translate and rotate to specified rigid body coordinates 781c 782 if (mode .eq. 16) then 783 xcm = 0.0d0 784 ycm = 0.0d0 785 zcm = 0.0d0 786 phi = 0.0d0 787 theta = 0.0d0 788 psi = 0.0d0 789 call nextarg (string,exist) 790 if (exist) read (string,*,err=360,end=360) xcm 791 call nextarg (string,exist) 792 if (exist) read (string,*,err=360,end=360) ycm 793 call nextarg (string,exist) 794 if (exist) read (string,*,err=360,end=360) zcm 795 call nextarg (string,exist) 796 if (exist) read (string,*,err=360,end=360) phi 797 call nextarg (string,exist) 798 if (exist) read (string,*,err=360,end=360) theta 799 call nextarg (string,exist) 800 if (exist) read (string,*,err=360,end=360) psi 801 360 continue 802 if (min(xcm,ycm,zcm,phi,theta,psi).eq.0.0d0 .and. 803 & max(xcm,ycm,zcm,phi,theta,psi).eq.0.0d0) then 804 write (iout,370) 805 370 format (/,' Enter Rigid Body Coordinates : ',$) 806 read (input,380) record 807 380 format (a240) 808 read (record,*,err=390,end=390) xcm,ycm,zcm,phi,theta,psi 809 390 continue 810 end if 811 call inertia (2) 812 phi = phi / radian 813 theta = theta / radian 814 psi = psi / radian 815 cphi = cos(phi) 816 sphi = sin(phi) 817 ctheta = cos(theta) 818 stheta = sin(theta) 819 cpsi = cos(psi) 820 spsi = sin(psi) 821 a(1,1) = ctheta * cphi 822 a(2,1) = spsi*stheta*cphi - cpsi*sphi 823 a(3,1) = cpsi*stheta*cphi + spsi*sphi 824 a(1,2) = ctheta * sphi 825 a(2,2) = spsi*stheta*sphi + cpsi*cphi 826 a(3,2) = cpsi*stheta*sphi - spsi*cphi 827 a(1,3) = -stheta 828 a(2,3) = ctheta * spsi 829 a(3,3) = ctheta * cpsi 830 do while (.not. abort) 831 do i = 1, n 832 xorig = x(i) 833 yorig = y(i) 834 zorig = z(i) 835 x(i) = a(1,1)*xorig + a(2,1)*yorig + a(3,1)*zorig + xcm 836 y(i) = a(1,2)*xorig + a(2,2)*yorig + a(3,2)*zorig + ycm 837 z(i) = a(1,3)*xorig + a(2,3)*yorig + a(3,3)*zorig + zcm 838 end do 839 call makeref (1) 840 call readxyz (ixyz) 841 if (.not. abort) multi = .true. 842 if (multi) then 843 call makeref (2) 844 call getref (1) 845 call prtmod (imod,offset) 846 call getref (2) 847 end if 848 end do 849 if (.not. multi) then 850 call getref (1) 851 goto 40 852 end if 853 end if 854c 855c move stray molecules back into original periodic box 856c 857 if (mode .eq. 17) then 858 do while (.not. abort) 859 call unitcell 860 if (use_bounds) then 861 call lattice 862 call molecule 863 call bounds 864 end if 865 call makeref (1) 866 call readxyz (ixyz) 867 if (.not. abort) multi = .true. 868 if (multi) then 869 call makeref (2) 870 call getref (1) 871 call prtmod (imod,offset) 872 call getref (2) 873 end if 874 end do 875 if (.not. multi) then 876 call getref (1) 877 goto 40 878 end if 879 end if 880c 881c remove molecules to trim periodic box to smaller size 882c 883 if (mode .eq. 18) then 884 xnew = 0.0d0 885 ynew = 0.0d0 886 znew = 0.0d0 887 call nextarg (string,exist) 888 if (exist) read (string,*,err=400,end=400) xnew 889 call nextarg (string,exist) 890 if (exist) read (string,*,err=400,end=400) ynew 891 call nextarg (string,exist) 892 if (exist) read (string,*,err=400,end=400) znew 893 400 continue 894 do while (xnew .eq. 0.0d0) 895 write (iout,410) 896 410 format (/,' Enter Periodic Box Dimensions (X,Y,Z) : ',$) 897 read (input,420) record 898 420 format (a240) 899 read (record,*,err=430,end=430) xnew,ynew,znew 900 430 continue 901 end do 902 if (ynew .eq. 0.0d0) ynew = xnew 903 if (znew .eq. 0.0d0) znew = xnew 904 xbox = xnew 905 ybox = ynew 906 zbox = znew 907 call lattice 908 call molecule 909 allocate (list(n)) 910 allocate (keep(n)) 911 do while (.not. abort) 912 do i = 1, n 913 list(i) = 1 914 end do 915 do i = 1, nmol 916 init = imol(1,i) 917 stop = imol(2,i) 918 xcm = 0.0d0 919 ycm = 0.0d0 920 zcm = 0.0d0 921 do j = init, stop 922 k = kmol(j) 923 weigh = mass(k) 924 xcm = xcm + x(k)*weigh 925 ycm = ycm + y(k)*weigh 926 zcm = zcm + z(k)*weigh 927 end do 928 weigh = molmass(i) 929 xcm = xcm / weigh 930 ycm = ycm / weigh 931 zcm = zcm / weigh 932 if (abs(xcm).gt.xbox2 .or. abs(ycm).gt.ybox2 933 & .or. abs(zcm).gt.zbox2) then 934 do j = init, stop 935 k = kmol(j) 936 list(k) = 0 937 end do 938 end if 939 end do 940 k = 0 941 do i = 1, n 942 if (list(i) .ne. 0) then 943 k = k + 1 944 keep(k) = i 945 list(i) = k 946 end if 947 end do 948 n = k 949 do k = 1, n 950 i = keep(k) 951 name(k) = name(i) 952 x(k) = x(i) 953 y(k) = y(i) 954 z(k) = z(i) 955 type(k) = type(i) 956 n12(k) = n12(i) 957 do j = 1, n12(k) 958 i12(j,k) = list(i12(j,i)) 959 end do 960 end do 961 call makeref (1) 962 call readxyz (ixyz) 963 if (.not. abort) then 964 multi = .true. 965 xbox = xnew 966 ybox = ynew 967 zbox = znew 968 call lattice 969 call molecule 970 end if 971 if (multi) then 972 call makeref (2) 973 call getref (1) 974 call prtmod (imod,offset) 975 call getref (2) 976 end if 977 end do 978 deallocate (list) 979 deallocate (keep) 980 if (.not. multi) then 981 call getref (1) 982 goto 40 983 end if 984 end if 985c 986c trim cube to truncated octahedron or rhombic dodecahedron 987c 988 if (mode.eq.19 .or. mode.eq.20) then 989 call unitcell 990 dowhile (xbox .eq. 0.0d0) 991 write (iout,440) 992 440 format (/,' Enter Edge Length of Cubic Periodic Box : ',$) 993 read (input,450) record 994 450 format (a240) 995 read (record,*,err=460,end=460) xbox 996 460 continue 997 end do 998 if (mode .eq. 18) octahedron = .true. 999 if (mode .eq. 19) dodecadron = .true. 1000 call molecule 1001 allocate (list(n)) 1002 allocate (keep(n)) 1003 do while (.not. abort) 1004 do i = 1, n 1005 list(i) = 1 1006 end do 1007 do i = 1, nmol 1008 init = imol(1,i) 1009 stop = imol(2,i) 1010 xcm = 0.0d0 1011 ycm = 0.0d0 1012 zcm = 0.0d0 1013 do j = init, stop 1014 k = kmol(j) 1015 weigh = mass(k) 1016 xcm = xcm + x(k)*weigh 1017 ycm = ycm + y(k)*weigh 1018 zcm = zcm + z(k)*weigh 1019 end do 1020 weigh = molmass(i) 1021 xcm = xcm / weigh 1022 ycm = ycm / weigh 1023 zcm = zcm / weigh 1024 if (octahedron) then 1025 xcm = xcm - xbox*nint(xcm/xbox) 1026 ycm = ycm - ybox*nint(ycm/ybox) 1027 zcm = zcm - zbox*nint(zcm/zbox) 1028 if (abs(xcm)+abs(ycm)+abs(zcm) .gt. 0.75*xbox) then 1029 do j = init, stop 1030 k = kmol(j) 1031 list(k) = 0 1032 end do 1033 end if 1034 else if (dodecadron) then 1035 xcm = xcm - xbox*nint(xcm/xbox) 1036 ycm = ycm - ybox*nint(ycm/ybox) 1037 zcm = zcm - root2*zbox*nint(zcm/(zbox*root2)) 1038 if (abs(xcm)+abs(ycm)+abs(root2*zcm) .gt. xbox) then 1039 do j = init, stop 1040 k = kmol(j) 1041 list(k) = 0 1042 end do 1043 end if 1044 end if 1045 end do 1046 k = 0 1047 do i = 1, n 1048 if (list(i) .ne. 0) then 1049 k = k + 1 1050 keep(k) = i 1051 list(i) = k 1052 end if 1053 end do 1054 n = k 1055 do k = 1, n 1056 i = keep(k) 1057 name(k) = name(i) 1058 x(k) = x(i) 1059 y(k) = y(i) 1060 z(k) = z(i) 1061 type(k) = type(i) 1062 n12(k) = n12(i) 1063 do j = 1, n12(k) 1064 i12(j,k) = list(i12(j,i)) 1065 end do 1066 end do 1067 call makeref (1) 1068 call readxyz (ixyz) 1069 if (.not. abort) multi = .true. 1070 if (multi) then 1071 call makeref (2) 1072 call getref (1) 1073 call prtmod (imod,offset) 1074 call getref (2) 1075 end if 1076 end do 1077 deallocate (list) 1078 deallocate (keep) 1079 if (.not. multi) then 1080 call getref (1) 1081 goto 40 1082 end if 1083 end if 1084c 1085c append a second file to the current coordinates file 1086c 1087 if (mode .eq. 21) then 1088 append = .false. 1089 do while (.not. abort) 1090 call makeref (1) 1091 if (append) then 1092 call getref (3) 1093 else 1094 call getxyz 1095 call makeref (3) 1096 append = .true. 1097 end if 1098 call merge (1) 1099 call makeref (1) 1100 call readxyz (ixyz) 1101 if (.not. abort) multi = .true. 1102 if (multi) then 1103 call makeref (2) 1104 call getref (1) 1105 call prtmod (imod,offset) 1106 call getref (2) 1107 end if 1108 end do 1109 if (.not. multi) then 1110 call getref (1) 1111 goto 40 1112 end if 1113 end if 1114c 1115c create random box full of the current coordinates file 1116c 1117 if (mode .eq. 22) then 1118 call makebox 1119 end if 1120c 1121c solvate the current system by insertion into a solvent box 1122c 1123 if (mode .eq. 23) then 1124 call soak 1125 end if 1126c 1127c replace random solvent molecules outside solute with ions 1128c 1129 if (mode .eq. 24) then 1130 call molecule 1131 call addions 1132 end if 1133c 1134c output final coordinates for single frame and print info 1135c 1136 if (opened .and. .not.multi) then 1137 call prtmod (imod,offset) 1138 end if 1139 if (opened) then 1140 close (unit=imod) 1141 write (iout,470) modfile(1:trimtext(modfile)) 1142 470 format (/,' New Coordinates Written To : ',a) 1143 end if 1144 close (unit=ixyz) 1145c 1146c perform any final tasks before program exit 1147c 1148 call final 1149 end 1150c 1151c 1152c ################################################################ 1153c ## ## 1154c ## subroutine prtmod -- Cartesian coordinates with offset ## 1155c ## ## 1156c ################################################################ 1157c 1158c 1159c "prtmod" writes out a set of modified Cartesian coordinates 1160c with an optional atom number offset to an external disk file 1161c 1162c 1163 subroutine prtmod (imod,offset) 1164 use atomid 1165 use atoms 1166 use bound 1167 use boxes 1168 use couple 1169 use inform 1170 use titles 1171 implicit none 1172 integer i,k,imod 1173 integer offset 1174 integer size,crdsiz 1175 real*8 crdmin,crdmax 1176 character*2 atmc 1177 character*2 crdc 1178 character*2 digc 1179 character*25 fstr 1180c 1181c 1182c check for large systems needing extended formatting 1183c 1184 atmc = 'i6' 1185 if (n+offset .ge. 100000) atmc = 'i7' 1186 if (n+offset .ge. 1000000) atmc = 'i8' 1187 crdmin = 0.0d0 1188 crdmax = 0.0d0 1189 do i = 1, n 1190 crdmin = min(crdmin,x(i),y(i),z(i)) 1191 crdmax = max(crdmax,x(i),y(i),z(i)) 1192 end do 1193 crdsiz = 6 1194 if (crdmin .le. -1000.0d0) crdsiz = 7 1195 if (crdmax .ge. 10000.0d0) crdsiz = 7 1196 if (crdmin .le. -10000.0d0) crdsiz = 8 1197 if (crdmax .ge. 100000.0d0) crdsiz = 8 1198 crdsiz = crdsiz + max(6,digits) 1199 size = 0 1200 call numeral (crdsiz,crdc,size) 1201 if (digits .le. 6) then 1202 digc = '6 ' 1203 else if (digits .le. 8) then 1204 digc = '8' 1205 else 1206 digc = '10' 1207 end if 1208c 1209c write out the number of atoms and the title 1210c 1211 if (ltitle .eq. 0) then 1212 fstr = '('//atmc//')' 1213 write (imod,fstr(1:4)) n 1214 else 1215 fstr = '('//atmc//',2x,a)' 1216 write (imod,fstr(1:9)) n,title(1:ltitle) 1217 end if 1218c 1219c write out the periodic cell lengths and angles 1220c 1221 if (use_bounds) then 1222 fstr = '(1x,6f'//crdc//'.'//digc//')' 1223 write (imod,fstr) xbox,ybox,zbox,alpha,beta,gamma 1224 end if 1225c 1226c write out the coordinate line for each atom 1227c 1228 fstr = '('//atmc//',2x,a3,3f'//crdc// 1229 & '.'//digc//',i6,8'//atmc//')' 1230 do i = 1, n 1231 write (imod,fstr) i+offset,name(i),x(i),y(i),z(i),type(i), 1232 & (i12(k,i)+offset,k=1,n12(i)) 1233 end do 1234 return 1235 end 1236c 1237c 1238c ################################################################ 1239c ## ## 1240c ## subroutine makebox -- build periodic box from monomers ## 1241c ## ## 1242c ################################################################ 1243c 1244c 1245c "makebox" builds a periodic box of a desired size by randomly 1246c copying a specified number of monomers into a target box size, 1247c followed by optional excluded volume refinement 1248c 1249c 1250 subroutine makebox 1251 use atomid 1252 use atoms 1253 use boxes 1254 use couple 1255 use iounit 1256 implicit none 1257 integer i,j,k,m 1258 integer ncopy 1259 integer offset 1260 real*8 xcm,ycm,zcm 1261 real*8 phi,theta,psi 1262 real*8 cphi,ctheta,cpsi 1263 real*8 sphi,stheta,spsi 1264 real*8 random,reduce 1265 real*8 norm,weigh 1266 real*8, allocatable :: x0(:) 1267 real*8, allocatable :: y0(:) 1268 real*8, allocatable :: z0(:) 1269 real*8 a(3,3) 1270 logical exist,query 1271 logical refine 1272 character*1 answer 1273 character*240 record 1274 character*240 string 1275c 1276c 1277c get the number of copies of the monomer to be used 1278c 1279 ncopy = 0 1280 call nextarg (string,exist) 1281 if (exist) read (string,*,err=10,end=10) ncopy 1282 10 continue 1283 if (ncopy .eq. 0) then 1284 write (iout,20) 1285 20 format (/,' Enter Number of Copies to Put in Box : ',$) 1286 read (input,30) ncopy 1287 30 format (i10) 1288 end if 1289c 1290c find the size of the periodic box to be constructed 1291c 1292 xbox = 0.0d0 1293 ybox = 0.0d0 1294 zbox = 0.0d0 1295 call nextarg (string,exist) 1296 if (exist) read (string,*,err=40,end=40) xbox 1297 call nextarg (string,exist) 1298 if (exist) read (string,*,err=40,end=40) ybox 1299 call nextarg (string,exist) 1300 if (exist) read (string,*,err=40,end=40) zbox 1301 40 continue 1302 do while (xbox .eq. 0.0d0) 1303 write (iout,50) 1304 50 format (/,' Enter Periodic Box Dimensions (X,Y,Z) : ',$) 1305 read (input,60) record 1306 60 format (a240) 1307 read (record,*,err=70,end=70) xbox,ybox,zbox 1308 70 continue 1309 end do 1310 if (ybox .eq. 0.0d0) ybox = xbox 1311 if (zbox .eq. 0.0d0) zbox = xbox 1312 orthogonal = .true. 1313c 1314c decide whether to use excluded volume refinement 1315c 1316 refine = .true. 1317 answer = 'Y' 1318 query = .true. 1319 call nextarg (string,exist) 1320 if (exist) then 1321 read (string,*,err=80,end=80) answer 1322 query = .false. 1323 end if 1324 80 continue 1325 if (query) then 1326 write (iout,90) 1327 90 format (/,' Refine the Periodic Box Configuration', 1328 & ' [Y] : ',$) 1329 read (input,100) answer 1330 100 format (a1) 1331 end if 1332 call upcase (answer) 1333 if (answer .eq. 'N') refine = .false. 1334c 1335c center the monomer and reduce its size to avoid overlap 1336c 1337 xcm = 0.0d0 1338 ycm = 0.0d0 1339 zcm = 0.0d0 1340 norm = 0.0d0 1341 do i = 1, n 1342 weigh = mass(i) 1343 xcm = xcm + x(i)*weigh 1344 ycm = ycm + y(i)*weigh 1345 zcm = zcm + z(i)*weigh 1346 norm = norm + weigh 1347 end do 1348 xcm = xcm / norm 1349 ycm = ycm / norm 1350 zcm = zcm / norm 1351 allocate (x0(n)) 1352 allocate (y0(n)) 1353 allocate (z0(n)) 1354 reduce = 0.001d0 1355 do i = 1, n 1356 x(i) = x(i) - xcm 1357 y(i) = y(i) - ycm 1358 z(i) = z(i) - zcm 1359 if (refine) then 1360 x(i) = reduce * x(i) 1361 y(i) = reduce * y(i) 1362 z(i) = reduce * z(i) 1363 end if 1364 x0(i) = x(i) 1365 y0(i) = y(i) 1366 z0(i) = z(i) 1367 end do 1368c 1369c randomly place monomer copies in the periodic box 1370c 1371 do k = 1, ncopy 1372 offset = (k-1) * n 1373 xcm = xbox * (random()-0.5d0) 1374 ycm = ybox * (random()-0.5d0) 1375 zcm = zbox * (random()-0.5d0) 1376 phi = 360.0d0 * random () 1377 theta = 360.0d0 * random () 1378 psi = 360.0d0 * random () 1379 cphi = cos(phi) 1380 sphi = sin(phi) 1381 ctheta = cos(theta) 1382 stheta = sin(theta) 1383 cpsi = cos(psi) 1384 spsi = sin(psi) 1385 a(1,1) = ctheta * cphi 1386 a(2,1) = spsi*stheta*cphi - cpsi*sphi 1387 a(3,1) = cpsi*stheta*cphi + spsi*sphi 1388 a(1,2) = ctheta * sphi 1389 a(2,2) = spsi*stheta*sphi + cpsi*cphi 1390 a(3,2) = cpsi*stheta*sphi - spsi*cphi 1391 a(1,3) = -stheta 1392 a(2,3) = ctheta * spsi 1393 a(3,3) = ctheta * cpsi 1394 do i = 1, n 1395 j = i + offset 1396 name(j) = name(i) 1397 type(j) = type(i) 1398 mass(j) = mass(i) 1399 x(j) = a(1,1)*x0(i) + a(2,1)*y0(i) + a(3,1)*z0(i) + xcm 1400 y(j) = a(1,2)*x0(i) + a(2,2)*y0(i) + a(3,2)*z0(i) + ycm 1401 z(j) = a(1,3)*x0(i) + a(2,3)*y0(i) + a(3,3)*z0(i) + zcm 1402 n12(j) = n12(i) 1403 do m = 1, n12(i) 1404 i12(m,j) = i12(m,i) + offset 1405 end do 1406 end do 1407 end do 1408 deallocate (x0) 1409 deallocate (y0) 1410 deallocate (z0) 1411 offset = 0 1412 n = ncopy * n 1413c 1414c optionally perform excluded volume coordinate refinement 1415c 1416 if (refine) then 1417 call boxfix 1418 call bounds 1419 else 1420 call lattice 1421 call molecule 1422 call bounds 1423 end if 1424 return 1425 end 1426c 1427c 1428c ################################################################# 1429c ## ## 1430c ## subroutine boxfix -- expand molecules into periodic box ## 1431c ## ## 1432c ################################################################# 1433c 1434c 1435c "boxfix" uses minimization of valence and vdw potential energy 1436c to expand and refine a collection of solvent molecules in a 1437c periodic box 1438c 1439c 1440 subroutine boxfix 1441 use atomid 1442 use atoms 1443 use boxes 1444 use inform 1445 use limits 1446 use linmin 1447 use minima 1448 use neigh 1449 use output 1450 use potent 1451 use scales 1452 use vdw 1453 use vdwpot 1454 implicit none 1455 integer i,j,k,nvar 1456 integer ii,kk 1457 real*8 minimum 1458 real*8 boxfix1 1459 real*8 grdmin 1460 real*8 boxsiz 1461 real*8, allocatable :: xx(:) 1462 external boxfix1 1463 external optsave 1464c 1465c 1466c setup for minimization with only valence and vdw terms 1467c 1468 call mechanic 1469 call potoff 1470 use_bond = .true. 1471 use_angle = .true. 1472 use_opbend = .true. 1473 use_opdist = .true. 1474 use_improp = .true. 1475 use_imptor = .true. 1476 use_tors = .true. 1477 use_vdw = .true. 1478c 1479c set artificial Lennard-Jones vdw values for the system 1480c 1481 vdwtyp = 'LENNARD-JONES' 1482 nvdw = n 1483 do i = 1, n 1484 ivdw(i) = i 1485 jvdw(i) = class(i) 1486 ired(i) = i 1487 end do 1488 do i = 1, n-1 1489 ii = jvdw(i) 1490 do k = i+1, n 1491 kk = jvdw(k) 1492 if (atomic(i).eq.1 .and. atomic(k).eq.1) then 1493 radmin(ii,kk) = 2.90d0 1494 epsilon(ii,kk) = 0.016d0 1495 radmin4(ii,kk) = 2.90d0 1496 epsilon4(ii,kk) = 0.016d0 1497 else if (atomic(i).eq.1 .or. atomic(k).eq.1) then 1498 radmin(ii,kk) = 3.35d0 1499 epsilon(ii,kk) = 0.040d0 1500 radmin4(ii,kk) = 3.35d0 1501 epsilon4(ii,kk) = 0.040d0 1502 else 1503 radmin(ii,kk) = 3.80d0 1504 epsilon(ii,kk) = 0.100d0 1505 radmin4(ii,kk) = 3.80d0 1506 epsilon4(ii,kk) = 0.100d0 1507 end if 1508 end do 1509 end do 1510c 1511c cutoff values and neighbor lists for vdw interactions 1512c 1513 use_list = .false. 1514 use_vlist = .false. 1515 vdwcut = 5.0d0 1516 vdwtaper = 4.5d0 1517 lbuffer = 1.0d0 1518 boxsiz = min(xbox,ybox,zbox) 1519 if (boxsiz .gt. 2.0d0*(vdwcut+lbuffer)) then 1520 use_list = .true. 1521 use_vlist = .true. 1522 dovlst = .true. 1523 lbuf2 = (0.5d0*lbuffer)**2 1524 vbuf2 = (vdwcut+lbuffer)**2 1525 vbufx = (vdwcut+2.0d0*lbuffer)**2 1526 maxvlst = int(sqrt(vbuf2)**3) + 100 1527 end if 1528c 1529c perform dynamic allocation of some global arrays 1530c 1531 if (use_vlist) then 1532 if (allocated(nvlst)) deallocate (nvlst) 1533 if (allocated(vlst)) deallocate (vlst) 1534 if (allocated(xvold)) deallocate (xvold) 1535 if (allocated(yvold)) deallocate (yvold) 1536 if (allocated(zvold)) deallocate (zvold) 1537 allocate (nvlst(n)) 1538 allocate (vlst(maxvlst,n)) 1539 allocate (xvold(n)) 1540 allocate (yvold(n)) 1541 allocate (zvold(n)) 1542 end if 1543c 1544c perform dynamic allocation of some global arrays 1545c 1546 if (.not. allocated(scale)) allocate (scale(3*n)) 1547c 1548c mark for use of all atoms, and set scale factors 1549c 1550 nvar = 0 1551 do i = 1, n 1552 do j = 1, 3 1553 nvar = nvar + 1 1554 scale(nvar) = 12.0d0 1555 end do 1556 end do 1557c 1558c perform dynamic allocation of some local arrays 1559c 1560 allocate (xx(nvar)) 1561c 1562c scale the coordinates of each active atom 1563c 1564 nvar = 0 1565 do i = 1, n 1566 nvar = nvar + 1 1567 xx(nvar) = x(i) * scale(nvar) 1568 nvar = nvar + 1 1569 xx(nvar) = y(i) * scale(nvar) 1570 nvar = nvar + 1 1571 xx(nvar) = z(i) * scale(nvar) 1572 end do 1573c 1574c make the call to the optimization routine 1575c 1576 iprint = 100 1577 maxiter = 10000 1578 stpmax = 10.0d0 1579 grdmin = 1.0d0 1580 coordtype = 'NONE' 1581 call lbfgs (nvar,xx,minimum,grdmin,boxfix1,optsave) 1582c 1583c unscale the final coordinates for active atoms 1584c 1585 nvar = 0 1586 do i = 1, n 1587 nvar = nvar + 1 1588 x(i) = xx(nvar) / scale(nvar) 1589 nvar = nvar + 1 1590 y(i) = xx(nvar) / scale(nvar) 1591 nvar = nvar + 1 1592 z(i) = xx(nvar) / scale(nvar) 1593 end do 1594c 1595c perform deallocation of some local arrays 1596c 1597 deallocate (xx) 1598 return 1599 end 1600c 1601c 1602c ############################################################ 1603c ## ## 1604c ## function boxfix1 -- energy and gradient for boxfix ## 1605c ## ## 1606c ############################################################ 1607c 1608c 1609c "boxfix1" is a service routine that computes the energy and 1610c gradient during refinement of a periodic box 1611c 1612c 1613 function boxfix1 (xx,g) 1614 use atoms 1615 use energi 1616 use potent 1617 use repel 1618 use scales 1619 implicit none 1620 integer i,nvar 1621 real*8 e,boxfix1 1622 real*8 xx(*) 1623 real*8 g(*) 1624 real*8, allocatable :: derivs(:,:) 1625c 1626c 1627c convert optimization parameters to atomic coordinates 1628c 1629 nvar = 0 1630 do i = 1, n 1631 nvar = nvar + 1 1632 x(i) = xx(nvar) / scale(nvar) 1633 nvar = nvar + 1 1634 y(i) = xx(nvar) / scale(nvar) 1635 nvar = nvar + 1 1636 z(i) = xx(nvar) / scale(nvar) 1637 end do 1638c 1639c perform dynamic allocation of some local arrays 1640c 1641 allocate (derivs(3,n)) 1642c 1643c compute and store the energy and gradient 1644c 1645 call gradient (e,derivs) 1646 boxfix1 = e 1647c 1648c convert gradient components to optimization parameters 1649c 1650 nvar = 0 1651 do i = 1, n 1652 nvar = nvar + 1 1653 g(nvar) = derivs(1,i) / scale(nvar) 1654 nvar = nvar + 1 1655 g(nvar) = derivs(2,i) / scale(nvar) 1656 nvar = nvar + 1 1657 g(nvar) = derivs(3,i) / scale(nvar) 1658 end do 1659c 1660c perform deallocation of some local arrays 1661c 1662 deallocate (derivs) 1663 return 1664 end 1665c 1666c 1667c ############################################################## 1668c ## ## 1669c ## subroutine soak -- place a solute into a solvent box ## 1670c ## ## 1671c ############################################################## 1672c 1673c 1674c "soak" takes a currently defined solute system and places 1675c it into a solvent box, with removal of any solvent molecules 1676c that overlap the solute 1677c 1678c 1679 subroutine soak 1680 use atomid 1681 use atoms 1682 use bound 1683 use boxes 1684 use couple 1685 use iounit 1686 use molcul 1687 use refer 1688 implicit none 1689 integer i,j,k 1690 integer ii,jj 1691 integer n12i,n12k 1692 integer isolv,icount 1693 integer ntot,freeunit 1694 integer, allocatable :: map(:) 1695 real*8 xi,yi,zi 1696 real*8 xr,yr,zr 1697 real*8 rik2,close2 1698 real*8 dxx,dxx2 1699 real*8 dxh,dxh2 1700 real*8 dhh,dhh2 1701 logical exist,header 1702 logical, allocatable :: remove(:) 1703 character*240 solvfile 1704 external merge 1705c 1706c 1707c make a copy of the solute coordinates and connectivities 1708c 1709 call makeref (1) 1710c 1711c get the filename for the solvent box coordinates 1712c 1713 call nextarg (solvfile,exist) 1714 if (exist) then 1715 call basefile (solvfile) 1716 call suffix (solvfile,'xyz','old') 1717 inquire (file=solvfile,exist=exist) 1718 end if 1719 do while (.not. exist) 1720 write (iout,10) 1721 10 format (/,' Enter Name of Solvent Box Coordinates : ',$) 1722 read (input,20) solvfile 1723 20 format (a240) 1724 call basefile (solvfile) 1725 call suffix (solvfile,'xyz','old') 1726 inquire (file=solvfile,exist=exist) 1727 end do 1728c 1729c read the coordinate file containing the solvent atoms 1730c 1731 isolv = freeunit () 1732 open (unit=isolv,file=solvfile,status='old') 1733 rewind (unit=isolv) 1734 call readxyz (isolv) 1735 close (unit=isolv) 1736c 1737c combine solute and solvent into a single coordinate set 1738c 1739 call merge (1) 1740 call basefile (solvfile) 1741 call getkey 1742c 1743c reset the default values for unitcell dimensions 1744c 1745 xbox = 0.0d0 1746 ybox = 0.0d0 1747 zbox = 0.0d0 1748 alpha = 0.0d0 1749 beta = 0.0d0 1750 gamma = 0.0d0 1751c 1752c count number of molecules and set lattice parameters 1753c 1754 call molecule 1755 call unitcell 1756 call lattice 1757c 1758c set distance cutoffs for solute-solvent close contacts 1759c 1760 dxx = 2.40d0 1761 dxh = 2.19d0 1762 dhh = 1.82d0 1763 dxx2 = dxx * dxx 1764 dxh2 = dxh * dxh 1765 dhh2 = dhh * dhh 1766c 1767c perform dynamic allocation of some local arrays 1768c 1769 allocate (map(n)) 1770 allocate (remove(nmol)) 1771c 1772c initialize the list of solvent molecules to be deleted 1773c 1774 do i = 1, nmol 1775 remove(i) = .false. 1776 end do 1777c 1778c print header information when processing large systems 1779c 1780 icount = 0 1781 header = .true. 1782 if (n-nref(1) .ge. 10000) then 1783 write (iout,30) 1784 30 format (/,' Scan for Solvent Molecules to be Removed :') 1785 end if 1786c 1787c OpenMP directives for the major loop structure 1788c 1789!$OMP PARALLEL default(private) 1790!$OMP& shared(nref,n,x,y,z,n12,molcule,dxx2,dxh2,dhh2,remove, 1791!$OMP& header,icount,iout) 1792!$OMP DO schedule(guided) 1793c 1794c search for close contacts between solute and solvent 1795c 1796 do i = nref(1)+1, n 1797 if (.not. remove(molcule(i))) then 1798 xi = x(i) 1799 yi = y(i) 1800 zi = z(i) 1801 n12i = n12(i) 1802 do k = 1, nref(1) 1803 n12k = n12(k) 1804 xr = x(k) - xi 1805 yr = y(k) - yi 1806 zr = z(k) - zi 1807 call imagen (xr,yr,zr) 1808 rik2 = xr*xr + yr*yr + zr*zr 1809 if (n12i.gt.1 .and. n12k.gt.1) then 1810 close2 = dxx2 1811 else if (n12i.gt.1 .or. n12k.gt.1) then 1812 close2 = dxh2 1813 else 1814 close2 = dhh2 1815 end if 1816 if (rik2 .lt. close2) then 1817 remove(molcule(i)) = .true. 1818 goto 40 1819 end if 1820 end do 1821 40 continue 1822 end if 1823 icount = icount + 1 1824 if (mod(icount,10000) .eq. 0) then 1825 if (header) then 1826 header = .false. 1827 write (iout,50) 1828 50 format () 1829 end if 1830 write (iout,60) 10000*(icount/10000) 1831 60 format (' Solvent Atoms Processed',i15) 1832 end if 1833 end do 1834c 1835c OpenMP directives for the major loop structure 1836c 1837!$OMP END DO 1838!$OMP END PARALLEL 1839c 1840c print final status when processing large systems 1841c 1842 icount = n - nref(1) 1843 if (mod(icount,10000).ne.0 .and. icount.gt.10000) then 1844 write (iout,70) icount 1845 70 format (' Solvent Atoms Processed',i15) 1846 end if 1847c 1848c delete solvent molecules that are too close to the solute 1849c 1850 k = nref(1) 1851 ntot = k 1852 do i = nref(1)+1, n 1853 map(i) = 0 1854 if (.not. remove(molcule(i))) then 1855 k = k + 1 1856 map(i) = k 1857 ntot = k 1858 end if 1859 end do 1860 do i = nref(1)+1, n 1861 ii = map(i) 1862 if (ii .ne. 0) then 1863 x(ii) = x(i) 1864 y(ii) = y(i) 1865 z(ii) = z(i) 1866 name(ii) = name(i) 1867 type(ii) = type(i) 1868 k = 0 1869 do j = 1, n12(i) 1870 jj = map(i12(j,i)) 1871 if (jj .ne. 0) then 1872 k = k + 1 1873 i12(k,ii) = jj 1874 end if 1875 end do 1876 n12(ii) = k 1877 end if 1878 end do 1879 n = ntot 1880c 1881c perform deallocation of some local arrays 1882c 1883 deallocate (map) 1884 deallocate (remove) 1885 return 1886 end 1887c 1888c 1889c ############################################################### 1890c ## ## 1891c ## subroutine addions -- placement of ions around solute ## 1892c ## ## 1893c ############################################################### 1894c 1895c 1896c "addions" takes a currently defined solvated system and 1897c places ions, with removal of solvent molecules 1898c 1899c 1900 subroutine addions 1901 use atomid 1902 use atoms 1903 use couple 1904 use iounit 1905 use katoms 1906 use molcul 1907 implicit none 1908 integer i,j,k 1909 integer nsolute,size 1910 integer start,stop 1911 integer icount,iontyp 1912 integer ncopy,ranatm 1913 integer, allocatable :: list(:) 1914 integer, allocatable :: isolute(:) 1915 real*8 xi,yi,zi 1916 real*8 xr,yr,zr,rik2 1917 real*8 close,close2 1918 real*8 rand,random 1919 real*8 weigh,xmid,ymid,zmid 1920 real*8, allocatable :: xion(:) 1921 real*8, allocatable :: yion(:) 1922 real*8, allocatable :: zion(:) 1923 logical exist,header,done 1924 logical, allocatable :: remove(:) 1925 character*240 record 1926 character*240 string 1927 external random 1928c 1929c 1930c perform dynamic allocation of some local arrays 1931c 1932 size = 40 1933 allocate (list(size)) 1934 allocate (isolute(n)) 1935c 1936c get the range atoms numbers constituting the solute 1937c 1938 do i = 1, size 1939 list(i) = 0 1940 end do 1941 i = 0 1942 do while (exist) 1943 call nextarg (string,exist) 1944 if (exist) then 1945 read (string,*,err=10,end=10) list(i+1) 1946 i = i + 1 1947 end if 1948 end do 1949 10 continue 1950 if (i .eq. 0) then 1951 write (iout,20) 1952 20 format (/,' Enter Atom Numbers in Solute Molecules : ',$) 1953 read (input,30) record 1954 30 format (a240) 1955 read (record,*,err=40,end=40) (list(i),i=1,size) 1956 40 continue 1957 end if 1958 i = 1 1959 nsolute = 0 1960 do while (list(i) .ne. 0) 1961 list(i) = max(-n,min(n,list(i))) 1962 if (list(i) .gt. 0) then 1963 k = list(i) 1964 nsolute = nsolute + 1 1965 isolute(nsolute) = k 1966 i = i + 1 1967 else 1968 list(i+1) = max(-n,min(n,list(i+1))) 1969 do k = abs(list(i)), abs(list(i+1)) 1970 nsolute = nsolute + 1 1971 isolute(nsolute) = k 1972 end do 1973 i = i + 2 1974 end if 1975 end do 1976c 1977c get the atom type of ion to be added and number of copies 1978c 1979 50 continue 1980 iontyp = 0 1981 ncopy = 0 1982 call nextarg (string,exist) 1983 if (exist) read (string,*,err=60,end=60) iontyp 1984 call nextarg (string,exist) 1985 if (exist) read (string,*,err=60,end=60) ncopy 1986 60 continue 1987 if (iontyp.eq.0 .or. ncopy.eq.0) then 1988 write (iout,70) 1989 70 format (/,' Enter Ion Atom Type Number & Copies to Add : ',$) 1990 read (input,80) record 1991 80 format (a240) 1992 end if 1993 read (record,*,err=50,end=50) iontyp,ncopy 1994c 1995c set minimum distance cutoff for solute-ion contacts 1996c 1997 close = 6.0d0 1998 close2 = close * close 1999c 2000c perform dynamic allocation of some local arrays 2001c 2002 allocate (remove(nmol)) 2003 allocate (xion(ncopy)) 2004 allocate (yion(ncopy)) 2005 allocate (zion(ncopy)) 2006c 2007c initialize the list of solvent molecules to be deleted 2008c 2009 do i = 1, nmol 2010 remove(i) = .false. 2011 end do 2012c 2013c print header information when processing large systems 2014c 2015 icount = 0 2016 header = .true. 2017 if (n .ge. 10000) then 2018 write (iout,90) 2019 90 format (/,' Scan for Available Locations to Place Ions :') 2020 end if 2021c 2022c OpenMP directives for the major loop structure 2023c 2024!$OMP PARALLEL default(private) 2025!$OMP& shared(n,x,y,z,molcule,close2,remove,header,nsolute, 2026!$OMP& isolute,icount,iout) 2027!$OMP DO schedule(guided) 2028c 2029c search for short distance between solute and solvent 2030c 2031 do i = 1, n 2032 if (.not. remove(molcule(i))) then 2033 xi = x(i) 2034 yi = y(i) 2035 zi = z(i) 2036 do k = 1, nsolute 2037 j = isolute(k) 2038 xr = x(j) - xi 2039 yr = y(j) - yi 2040 zr = z(j) - zi 2041 call imagen (xr,yr,zr) 2042 rik2 = xr*xr + yr*yr + zr*zr 2043 if (rik2 .lt. close2) then 2044 remove(molcule(i)) = .true. 2045 goto 100 2046 end if 2047 end do 2048 100 continue 2049 end if 2050 icount = icount + 1 2051 if (mod(icount,10000) .eq. 0) then 2052 if (header) then 2053 header = .false. 2054 write (iout,110) 2055 110 format () 2056 end if 2057 write (iout,120) 10000*(icount/10000) 2058 120 format (' Solvent Atoms Processed',i15) 2059 end if 2060 end do 2061c 2062c OpenMP directives for the major loop structure 2063c 2064!$OMP END DO 2065!$OMP END PARALLEL 2066c 2067c perform deallocation of some local arrays 2068c 2069 deallocate (list) 2070 deallocate (isolute) 2071c 2072c print final status when processing large systems 2073c 2074 if (mod(n,10000).ne.0 .and. n.gt.10000) then 2075 write (iout,130) n 2076 130 format (' Solvent Atoms Processed',i15) 2077 end if 2078c 2079c randomly replace the solvent molecules with ions 2080c 2081 do i = 1, ncopy 2082 done = .false. 2083 do while (.not. done) 2084 rand = random () 2085 ranatm = int(rand*dble(n)) + 1 2086c 2087c check solute distance, then delete polyatomic molecule 2088c 2089 if (.not. remove(molcule(ranatm))) then 2090 start = imol(1,molcule(ranatm)) 2091 stop = imol(2,molcule(ranatm)) 2092 if (start .eq. stop) then 2093 done = .false. 2094 else 2095 done = .true. 2096 xmid = 0.0d0 2097 ymid = 0.0d0 2098 zmid = 0.0d0 2099 do k = stop, start, -1 2100 weigh = mass(k) 2101 xmid = xmid + x(k)*weigh 2102 ymid = ymid + y(k)*weigh 2103 zmid = zmid + z(k)*weigh 2104 call delete (k) 2105 end do 2106 weigh = molmass(molcule(ranatm)) 2107 xion(i) = xmid / weigh 2108 yion(i) = ymid / weigh 2109 zion(i) = zmid / weigh 2110 end if 2111 end if 2112 end do 2113 end do 2114c 2115c insert new monoatomic ions at saved centers of mass 2116c 2117 do i = 1, ncopy 2118 n = n + 1 2119 name(n) = symbol(iontyp) 2120 x(n) = xion(i) 2121 y(n) = yion(i) 2122 z(n) = zion(i) 2123 type(n) = iontyp 2124 n12(n) = 0 2125 end do 2126c 2127c perform deallocation of some local arrays 2128c 2129 deallocate (remove) 2130 deallocate (xion) 2131 deallocate (yion) 2132 deallocate (zion) 2133 return 2134 end 2135