1c 2c 3c ############################################################## 4c ## COPYRIGHT (C) 2008 by C. Wu, Zhifeng Jing & Jay Ponder ## 5c ## All Rights Reserved ## 6c ############################################################## 7c 8c ############################################################## 9c ## ## 10c ## program potential -- electrostatic potential utility ## 11c ## ## 12c ############################################################## 13c 14c 15c "potential" calculates the electrostatic potential for a 16c molecule at a set of grid points; optionally compares to a 17c target potential or optimizes electrostatic parameters 18c 19c 20 program potential 21 use atoms 22 use charge 23 use files 24 use inform 25 use iounit 26 use keys 27 use mpole 28 use neigh 29 use potent 30 use potfit 31 use titles 32 use units 33 implicit none 34 integer i,j,k 35 integer ixyz,ipot 36 integer igrd,icub 37 integer next,mode 38 integer nvar,nmodel 39 integer nresid 40 integer maxpgrd 41 integer nglist,nflist 42 integer freeunit 43 integer trimtext 44 integer, allocatable :: glist(:) 45 integer, allocatable :: flist(:) 46 real*8 xi,yi,zi,pot 47 real*8 x0,y0,z0 48 real*8 xx0,xy0,xz0 49 real*8 yy0,yz0,zz0 50 real*8 grdmin 51 real*8, allocatable :: xx(:) 52 real*8, allocatable :: xlo(:) 53 real*8, allocatable :: xhi(:) 54 real*8, allocatable :: gc(:) 55 real*8, allocatable :: presid(:) 56 real*8, allocatable :: pjac(:,:) 57 logical exist,query 58 logical dogrid,docube 59 logical domodel,dopair 60 logical dotarget,dofit 61 logical dofull 62 logical, allocatable :: tmpchg(:) 63 logical, allocatable :: tmppol(:) 64 character*1 answer 65 character*20 keyword 66 character*240 record 67 character*240 string 68 character*240 xyzfile 69 character*240 xyz2file 70 character*240 potfile 71 character*240 gridfile 72 character*240 cubefile 73 external fitrsd 74 external potwrt 75c 76c 77c setup the computation and assign some default values 78c 79 call initial 80 nmodel = 1 81 dogrid = .false. 82 docube = .false. 83 domodel = .false. 84 dopair = .false. 85 dotarget = .false. 86 dofit = .false. 87 resptyp = 'ORIG' 88 wresp = 1.0d0 89c 90c perform dynamic allocation of some global arrays 91c 92 maxpgrd = 100000 93 allocate (ipgrid(maxpgrd,maxref)) 94 allocate (pgrid(3,maxpgrd,maxref)) 95 allocate (epot(2,maxpgrd,maxref)) 96c 97c initialize target molecular dipole and quadrupole values 98c 99 use_dpl = .false. 100 use_qpl = .false. 101 fit_mpl = .true. 102 fit_dpl = .true. 103 fit_qpl = .true. 104 do i = 1, maxref 105 xdpl0(i) = 0.0d0 106 ydpl0(i) = 0.0d0 107 zdpl0(i) = 0.0d0 108 xxqpl0(i) = 0.0d0 109 xyqpl0(i) = 0.0d0 110 xzqpl0(i) = 0.0d0 111 yyqpl0(i) = 0.0d0 112 yzqpl0(i) = 0.0d0 113 zzqpl0(i) = 0.0d0 114 end do 115c 116c find electrostatic potential manipulation to perform 117c 118 mode = 0 119 query = .true. 120 call nextarg (string,exist) 121 if (exist) then 122 read (string,*,err=10,end=10) mode 123 query = .false. 124 end if 125 10 continue 126 if (query) then 127 write (iout,20) 128 20 format (/,' The Tinker Electrostatic Potential Utility Can :', 129 & //,4x,'(1) Create an Input File for Gaussian CUBEGEN', 130 & /,4x,'(2) Get QM Potential from a Gaussian CUBE File', 131 & /,4x,'(3) Calculate the Model Potential for a System', 132 & /,4x,'(4) Compare Two Model Potentials for a System', 133 & /,4x,'(5) Compare a Model Potential to a Target Grid', 134 & /,4x,'(6) Fit Electrostatic Parameters to Target Grid') 135 do while (mode.lt.1 .or. mode.gt.6) 136 mode = 0 137 write (iout,30) 138 30 format (/,' Enter the Number of the Desired Choice : ',$) 139 read (input,40,err=50,end=50) mode 140 40 format (i10) 141 50 continue 142 end do 143 end if 144 if (mode .eq. 1) then 145 dogrid = .true. 146 else if (mode .eq. 2) then 147 docube = .true. 148 else if (mode .eq. 3) then 149 domodel = .true. 150 else if (mode .eq. 4) then 151 nmodel = 2 152 dopair = .true. 153 else if (mode .eq. 5) then 154 dotarget = .true. 155 else if (mode .eq. 6) then 156 dotarget = .true. 157 dofit = .true. 158 end if 159c 160c read electrostatic potential from a Gaussian CUBE file 161c 162 if (docube) then 163 call nextarg (cubefile,exist) 164 if (exist) then 165 call basefile (cubefile) 166 call suffix (cubefile,'cube','old') 167 inquire (file=cubefile,exist=exist) 168 end if 169 do while (.not. exist) 170 write (iout,60) 171 60 format (/,' Enter the Gaussian CUBE File Name : ',$) 172 read (input,70) cubefile 173 70 format (a240) 174 call basefile (cubefile) 175 call suffix (cubefile,'cube','old') 176 inquire (file=cubefile,exist=exist) 177 end do 178 icub = freeunit () 179 open (unit=icub,file=cubefile,status ='old') 180 rewind (unit=icub) 181 read (icub,80) title 182 80 format (1x,a240) 183 ltitle = trimtext (title) 184 read (icub,90) 185 90 format () 186 read (icub,100) n 187 100 format (i5) 188 read (icub,110) npgrid(1) 189 110 format (i5) 190 do i = 1, n+2 191 read (icub,120) 192 120 format () 193 end do 194 do i = 1, npgrid(1) 195 read (icub,130) record 196 130 format (a240) 197 read (record,*) xi,yi,zi,pot 198 pgrid(1,i,1) = xi 199 pgrid(2,i,1) = yi 200 pgrid(3,i,1) = zi 201 epot(1,i,1) = hartree * pot 202 end do 203 close (unit=icub) 204c 205c write the electrostatic potential to a Tinker POT file 206c 207 potfile = filename(1:leng) 208 call suffix (potfile,'pot','new') 209 ipot = freeunit () 210 open (unit=ipot,file=potfile,status ='new') 211 rewind (unit=ipot) 212 write (ipot,140) npgrid(1),title(1:ltitle) 213 140 format (i8,2x,a) 214 do i = 1, npgrid(1) 215 xi = pgrid(1,i,1) 216 yi = pgrid(2,i,1) 217 zi = pgrid(3,i,1) 218 pot = epot(1,i,1) 219 write (ipot,150) i,xi,yi,zi,pot 220 150 format (i8,3x,3f12.6,2x,f12.4) 221 end do 222 close (unit=ipot) 223 write (iout,160) potfile(1:trimtext(potfile)) 224 160 format (/,' Electrostatic Potential Written To : ',a) 225 goto 400 226 end if 227c 228c read the first structure and setup atom definitions 229c 230 call getxyz 231 call field 232 call katom 233c 234c reopen the structure file and get all the structures 235c 236 ixyz = freeunit () 237 xyzfile = filename 238 call suffix (xyzfile,'xyz','old') 239 open (unit=ixyz,file=xyzfile,status ='old') 240 rewind (unit=ixyz) 241 call readxyz (ixyz) 242 nconf = 0 243 namax = n 244 do while (.not. abort) 245 nconf = nconf + 1 246 call makeref (nconf) 247 call readxyz (ixyz) 248 namax = max(namax,n) 249 end do 250 close (unit=ixyz) 251 if (nconf .gt. 1) then 252 write (iout,170) nconf 253 170 format (/,' Structures Used for Potential Analysis :',3x,i16) 254 end if 255c 256c perform dynamic allocation of some global arrays 257c 258 allocate (gatm(namax)) 259 allocate (fatm(namax)) 260c 261c perform dynamic allocation of some local arrays 262c 263 allocate (glist(namax)) 264 allocate (flist(namax)) 265c 266c set defaults for the active grid atoms and fit atoms 267c 268 nglist = 0 269 nflist = 0 270 ngatm = namax 271 nfatm = namax 272 do i = 1, namax 273 glist(i) = 0 274 flist(i) = 0 275 gatm(i) = .true. 276 fatm(i) = .true. 277 end do 278c 279c get control parameters and target values from keyfile 280c 281 do i = 1, nkey 282 next = 1 283 record = keyline(i) 284 call gettext (record,keyword,next) 285 call upcase (keyword) 286 string = record(next:240) 287 if (keyword(1:16) .eq. 'POTENTIAL-ATOMS ') then 288 read (string,*,err=180,end=180) (glist(k),k=nglist+1,namax) 289 180 continue 290 do while (glist(nglist+1) .ne. 0) 291 nglist = nglist + 1 292 glist(nglist) = max(-namax,min(namax,glist(nglist))) 293 end do 294 else if (keyword(1:14) .eq. 'POTENTIAL-FIT ') then 295 read (string,*,err=190,end=190) (flist(k),k=nflist+1,namax) 296 190 continue 297 do while (flist(nflist+1) .ne. 0) 298 nflist = nflist + 1 299 flist(nflist) = max(-namax,min(namax,flist(nflist))) 300 end do 301 else if (keyword(1:9) .eq. 'RESPTYPE ') then 302 call getword (record,resptyp,next) 303 call upcase (resptyp) 304 else if (keyword(1:12) .eq. 'RESP-WEIGHT ') then 305 read (string,*,err=200,end=200) wresp 306 200 continue 307 else if (keyword(1:13) .eq. 'FIX-MONOPOLE ') then 308 fit_mpl = .false. 309 else if (keyword(1:11) .eq. 'FIX-DIPOLE ') then 310 fit_dpl = .false. 311 else if (keyword(1:15) .eq. 'FIX-QUADRUPOLE ') then 312 fit_qpl = .false. 313 else if (keyword(1:14) .eq. 'TARGET-DIPOLE ') then 314 use_dpl = .true. 315 k = 1 316 read (string,*,err=210,end=210) x0,y0,z0,k 317 210 continue 318 xdpl0(k) = x0 319 ydpl0(k) = y0 320 zdpl0(k) = z0 321 else if (keyword(1:18) .eq. 'TARGET-QUADRUPOLE ') then 322 use_qpl = .true. 323 k = 1 324 read (string,*,err=220,end=220) xx0,xy0,xz0,yy0,yz0,zz0,k 325 220 continue 326 xxqpl0(k) = xx0 327 xyqpl0(k) = xy0 328 xzqpl0(k) = xz0 329 yyqpl0(k) = yy0 330 yzqpl0(k) = yz0 331 zzqpl0(k) = zz0 332 end if 333 end do 334c 335c set active grid atoms to only those marked for use 336c 337 i = 1 338 do while (glist(i) .ne. 0) 339 if (i .eq. 1) then 340 ngatm = 0 341 do k = 1, namax 342 gatm(k) = .false. 343 end do 344 end if 345 if (glist(i) .gt. 0) then 346 k = glist(i) 347 if (.not. gatm(k)) then 348 gatm(k) = .true. 349 ngatm = ngatm + 1 350 end if 351 i = i + 1 352 else 353 do k = abs(glist(i)), abs(glist(i+1)) 354 if (.not. gatm(k)) then 355 gatm(k) = .true. 356 ngatm = ngatm + 1 357 end if 358 end do 359 i = i + 2 360 end if 361 end do 362c 363c set active fitting atoms to only those marked for use 364c 365 i = 1 366 do while (flist(i) .ne. 0) 367 if (i .eq. 1) then 368 nfatm = 0 369 do k = 1, namax 370 fatm(k) = .false. 371 end do 372 end if 373 if (flist(i) .gt. 0) then 374 k = flist(i) 375 if (.not. fatm(k)) then 376 fatm(k) = .true. 377 nfatm = nfatm + 1 378 end if 379 i = i + 1 380 else 381 do k = abs(flist(i)), abs(flist(i+1)) 382 if (.not. fatm(k)) then 383 fatm(k) = .true. 384 nfatm = nfatm + 1 385 end if 386 end do 387 i = i + 2 388 end if 389 end do 390c 391c perform deallocation of some local arrays 392c 393 deallocate (glist) 394 deallocate (flist) 395c 396c generate potential grid based on the molecular surface 397c 398 if (.not. dotarget) then 399 do i = 1, nconf 400 call getref (i) 401 call potgrid (i) 402 end do 403 end if 404c 405c get name of optional second structure for comparison 406c 407 if (dopair) then 408 call nextarg (xyz2file,exist) 409 if (exist) then 410 call basefile (xyz2file) 411 call suffix (xyz2file,'xyz','old') 412 inquire (file=xyz2file,exist=exist) 413 end if 414 do while (.not. exist) 415 write (iout,230) 416 230 format (/,' Enter Name of Second Coordinate File : ',$) 417 read (input,240) xyz2file 418 240 format (a240) 419 call basefile (xyz2file) 420 call suffix (xyz2file,'xyz','old') 421 inquire (file=xyz2file,exist=exist) 422 end do 423 end if 424c 425c get optional file with grid points and target potential 426c 427 if (dotarget) then 428 call nextarg (potfile,exist) 429 if (exist) then 430 call basefile (potfile) 431 call suffix (potfile,'pot','old') 432 inquire (file=potfile,exist=exist) 433 end if 434 do while (.not. exist) 435 write (iout,250) 436 250 format (/,' Enter Target Grid/Potential File Name : ',$) 437 read (input,260) potfile 438 260 format (a240) 439 call basefile (potfile) 440 call suffix (potfile,'pot','old') 441 inquire (file=potfile,exist=exist) 442 end do 443 end if 444c 445c decide whether to output potential at each grid point 446c 447 dofull = .false. 448 if (domodel .or. dopair .or. dotarget) then 449 call nextarg (answer,exist) 450 if (.not. exist) then 451 write (iout,270) 452 270 format (/,' Output Potential Value at Each Grid Point', 453 & ' [N] : ',$) 454 read (input,280) record 455 280 format (a240) 456 next = 1 457 call gettext (record,answer,next) 458 end if 459 call upcase (answer) 460 if (answer .eq. 'Y') dofull = .true. 461 end if 462c 463c read grid points where potential will be computed 464c 465 if (dotarget) then 466 ipot = freeunit () 467 open (unit=ipot,file=potfile,status='old') 468 rewind (unit=ipot) 469 do i = 1, nconf 470 call getref (i) 471 call readpot (ipot,i) 472 end do 473 close (unit=ipot) 474 end if 475c 476c output the number of potential grid points to be used 477c 478 do i = 1, nconf 479 if (i .eq. 1) then 480 write (iout,290) 481 290 format () 482 end if 483 if (npgrid(i) .gt. maxpgrd) then 484 write (iout,300) 485 300 format (' POTENTIAL -- Too many Grid Points;', 486 & ' Increase MAXGRID') 487 call fatal 488 else if (nconf .eq. 1) then 489 write (iout,310) npgrid(1) 490 310 format (' Electrostatic Potential Grid Points :',6x,i16) 491 else 492 write (iout,320) i,npgrid(i) 493 320 format (' Potential Grid Points for Structure',i4,' :', 494 & 2x,i16) 495 end if 496 end do 497c 498c output grid points at which to compute QM potential 499c 500 if (dogrid) then 501 igrd = freeunit () 502 gridfile = filename(1:leng) 503 call suffix (gridfile,'grid','new') 504 open (unit=igrd,file=gridfile,status='new') 505 do j = 1, nconf 506 do i = 1, npgrid(j) 507 xi = pgrid(1,i,j) 508 yi = pgrid(2,i,j) 509 zi = pgrid(3,i,j) 510 write (igrd,330) xi,yi,zi 511 330 format (3f15.8) 512 end do 513 end do 514 close (unit=igrd) 515 write (iout,340) gridfile(1:trimtext(gridfile)) 516 340 format (/,' Gaussian CUBEGEN Input Written To : ',a) 517 write (iout,350) 518 350 format (/,' Next, run the Gaussian CUBEGEN program; for', 519 & ' example:', 520 & //,' cubegen 0 potential=MP2 FILE.fchk FILE.cube', 521 & ' -5 h < FILE.grid', 522 & //,' Replace FILE with base file name and MP2 with', 523 & ' density label;', 524 & /,' After CUBEGEN, rerun Tinker POTENTIAL program', 525 & ' using Option 2') 526 end if 527c 528c get termination criterion for fitting as RMS gradient 529c 530 if (dofit) then 531 grdmin = -1.0d0 532 call nextarg (string,exist) 533 if (exist) read (string,*,err=360,end=360) grdmin 534 360 continue 535 if (grdmin .le. 0.0d0) then 536 write (iout,370) 537 370 format (/,' Enter RMS Gradient Termination Criterion', 538 & ' [0.01] : ',$) 539 read (input,380) grdmin 540 380 format (f20.0) 541 end if 542 if (grdmin .le. 0.0d0) grdmin = 0.01d0 543 end if 544c 545c print the parameter restraint value to be used in fitting 546c 547 if (dofit) then 548 write (iout,390) wresp 549 390 format (/,' Electrostatic Parameter Restraint Value :',f18.4) 550 end if 551c 552c setup the potential computation for alternative models 553c 554 if (.not. dogrid) then 555 do k = 1, nmodel 556 ixyz = freeunit () 557 if (k .eq. 1) then 558 call basefile (xyzfile) 559 open (unit=ixyz,file=xyzfile,status='old') 560 else 561 call basefile (xyz2file) 562 open (unit=ixyz,file=xyz2file,status='old') 563 end if 564 rewind (unit=ixyz) 565 do j = 1, nconf 566 call readxyz (ixyz) 567 call makeref (j) 568 end do 569 close (unit=ixyz) 570c 571c get potential for each structure and print statistics 572c 573 do j = 1, nconf 574 call getref (j) 575 call field 576 call setelect 577 if (use_mpole) then 578 call chkpole 579 call rotpole 580 end if 581 if (use_polar) then 582 domlst = .true. 583 doulst = .true. 584 call nblist 585 call induce 586 end if 587!$OMP PARALLEL default(private) shared(j,k,npgrid,pgrid,epot) 588!$OMP DO 589 do i = 1, npgrid(j) 590 xi = pgrid(1,i,j) 591 yi = pgrid(2,i,j) 592 zi = pgrid(3,i,j) 593 call potpoint (xi,yi,zi,pot) 594 epot(k,i,j) = pot 595 end do 596!$OMP END DO 597!$OMP END PARALLEL 598 end do 599 end do 600 call potstat (dofull,domodel,dopair,dotarget) 601 end if 602c 603c perform dynamic allocation of some global arrays 604c 605 if (dofit) then 606 allocate (fit0(12*nconf*namax)) 607 allocate (fchg(maxtyp)) 608 allocate (fpol(13,maxtyp)) 609 allocate (fitchg(maxtyp)) 610 allocate (fitpol(maxtyp)) 611c 612c perform dynamic allocation of some local arrays 613c 614 allocate (xx(12*nconf*namax)) 615 allocate (xlo(12*nconf*namax)) 616 allocate (xhi(12*nconf*namax)) 617 allocate (tmpchg(maxtyp)) 618 allocate (tmppol(maxtyp)) 619c 620c zero the keyfile length to avoid parameter reprocessing 621c 622 nkey = 0 623c 624c set residual count and optimization parameters with bounds 625c 626 do j = 1, maxtyp 627 fitchg(j) = .false. 628 fitpol(j) = .false. 629 end do 630 nvar = 0 631 nresid = 0 632 do j = 1, nconf 633 call getref (j) 634 call setelect 635 call prmvar (nvar,xx) 636 nresid = nresid + npgrid(j) 637 if (fit_mpl) nresid = nresid + 1 638 if (use_dpl) nresid = nresid + 3 639 if (use_qpl) nresid = nresid + 5 640 end do 641 nresid = nresid + nvar 642 do j = 1, nvar 643 fit0(j) = xx(j) 644 xlo(j) = xx(j) - 1000.0d0 645 xhi(j) = xx(j) + 1000.0d0 646 end do 647c 648c perform dynamic allocation of some local arrays 649c 650 allocate (gc(nvar)) 651 allocate (presid(nresid)) 652 allocate (pjac(nresid,nvar)) 653c 654c perform potential fit via least squares optimization 655c 656 call square (nvar,nresid,xlo,xhi,xx,presid,gc,pjac, 657 & grdmin,fitrsd,potwrt) 658c 659c set the final electrostatic parameter values 660c 661 nvar = 0 662 do j = 1, maxtyp 663 fitchg(j) = .false. 664 fitpol(j) = .false. 665 end do 666 do j = 1, nconf 667 call getref (j) 668 call setelect 669 next = nvar 670 do k = 1, maxtyp 671 tmpchg(k) = fitchg(k) 672 tmppol(k) = fitpol(k) 673 end do 674 call varprm (nvar,xx) 675 nvar = next 676 do k = 1, maxtyp 677 fitchg(k) = tmpchg(k) 678 fitpol(k) = tmppol(k) 679 end do 680 call prmvar (nvar,xx) 681 end do 682c 683c get potential for each structure and print statistics 684c 685 nvar = 0 686 do j = 1, maxtyp 687 fitchg(j) = .false. 688 fitpol(j) = .false. 689 end do 690 do j = 1, nconf 691 call getref (j) 692 call setelect 693 call varprm (nvar,xx) 694 call prmvar (nvar,xx) 695 if (use_mpole) then 696 call chkpole 697 call rotpole 698 end if 699 if (use_polar) then 700 domlst = .true. 701 doulst = .true. 702 call nblist 703 call induce 704 end if 705!$OMP PARALLEL default(private) shared(j,npgrid,pgrid,epot) 706!$OMP DO 707 do i = 1, npgrid(j) 708 xi = pgrid(1,i,j) 709 yi = pgrid(2,i,j) 710 zi = pgrid(3,i,j) 711 call potpoint (xi,yi,zi,pot) 712 epot(1,i,j) = pot 713 end do 714!$OMP END DO 715!$OMP END PARALLEL 716 end do 717 call prtfit 718 call potstat (dofull,domodel,dopair,dotarget) 719c 720c perform deallocation of some local arrays 721c 722 deallocate (xx) 723 deallocate (xlo) 724 deallocate (xhi) 725 deallocate (presid) 726 deallocate (gc) 727 deallocate (pjac) 728 deallocate (tmpchg) 729 deallocate (tmppol) 730 end if 731c 732c perform any final tasks before program exit 733c 734 400 continue 735 call final 736 end 737c 738c 739c ############################################################# 740c ## ## 741c ## subroutine readpot -- get and assign potential grid ## 742c ## ## 743c ############################################################# 744c 745c 746c "readpot" gets a set of grid points and target electrostatic 747c potential values from an external disk file 748c 749c 750 subroutine readpot (ipot,iconf) 751 use atoms 752 use katoms 753 use potfit 754 use ptable 755 implicit none 756 integer i,j,k 757 integer ipot,iconf 758 integer npoint,atn 759 real*8 xi,yi,zi 760 real*8 big,small 761 real*8 r2,dist 762 real*8, allocatable :: rad(:) 763 character*240 record 764c 765c 766c read the grid points and target potential from a file 767c 768 npoint = 0 769 read (ipot,10,err=20,end=20) record 770 10 format (a240) 771 read (record,*,err=20,end=20) npoint 772 20 continue 773 do i = 1, npoint 774 pgrid(1,i,iconf) = 0.0d0 775 pgrid(2,i,iconf) = 0.0d0 776 pgrid(3,i,iconf) = 0.0d0 777 epot(2,i,iconf) = 0.0d0 778 read (ipot,30,err=40,end=40) record 779 30 format (a240) 780 read (record,*,err=40,end=40) k,(pgrid(j,i,iconf),j=1,3), 781 & epot(2,i,iconf) 782 40 continue 783 end do 784c 785c perform dynamic allocation of some local arrays 786c 787 allocate (rad(n)) 788c 789c set base atomic radii from consensus vdw values 790c 791 do i = 1, n 792 rad(i) = 0.0d0 793 atn = atmnum(type(i)) 794 if (atn .ne. 0) rad(i) = vdwrad(atn) 795 if (rad(i) .eq. 0.0d0) rad(i) = 1.7d0 796 end do 797c 798c assign each grid point to atom on molecular surface 799c 800 big = 1000.0d0 801 do i = 1, npoint 802 small = big 803 xi = pgrid(1,i,iconf) 804 yi = pgrid(2,i,iconf) 805 zi = pgrid(3,i,iconf) 806 do k = 1, n 807 r2 = (xi-x(k))**2 + (yi-y(k))**2 + (zi-z(k))**2 808 dist = sqrt(r2) - rad(k) 809 if (dist .lt. small) then 810 small = dist 811 ipgrid(i,iconf) = k 812 end if 813 end do 814 end do 815c 816c perform deallocation of some local arrays 817c 818 deallocate (rad) 819c 820c use potential grid points only for active grid atoms 821c 822 k = 0 823 do i = 1, npoint 824 if (gatm(ipgrid(i,iconf))) then 825 k = k + 1 826 ipgrid(k,iconf) = ipgrid(i,iconf) 827 pgrid(1,k,iconf) = pgrid(1,i,iconf) 828 pgrid(2,k,iconf) = pgrid(2,i,iconf) 829 pgrid(3,k,iconf) = pgrid(3,i,iconf) 830 epot(2,k,iconf) = epot(2,i,iconf) 831 end if 832 end do 833 npgrid(iconf) = k 834 return 835 end 836c 837c 838c ############################################################## 839c ## ## 840c ## subroutine potgrid -- generate shells of grid points ## 841c ## ## 842c ############################################################## 843c 844c 845c "potgrid" generates electrostatic potential grid points in 846c radially distributed shells based on the molecular surface 847c 848c 849 subroutine potgrid (iconf) 850 use atoms 851 use iounit 852 use katoms 853 use keys 854 use math 855 use potfit 856 use ptable 857 implicit none 858 integer i,j,k,m 859 integer iconf,next 860 integer npoint,nshell 861 integer maxdot 862 integer ndot,atn 863 real*8 r2,rfactor 864 real*8 roffset 865 real*8 spacing 866 real*8 density 867 real*8 round 868 real*8 xi,yi,zi 869 real*8 xj,yj,zj 870 real*8 xr,yr,zr 871 real*8, allocatable :: rad(:) 872 real*8, allocatable :: rad2(:) 873 real*8, allocatable :: dot(:,:) 874 character*20 keyword 875 character*240 record 876 character*240 string 877c 878c 879c set default values for grid point generation parameters 880c 881 npoint = 0 882 nshell = 4 883 maxdot = 50000 884 spacing = 0.35d0 885 density = 4.0d0 * pi / spacing**2 886 rfactor = 1.0d0 887 roffset = 1.0d0 888 round = 0.000001d0 889c 890c check for keywords containing any altered parameters 891c 892 do i = 1, nkey 893 next = 1 894 record = keyline(i) 895 call gettext (record,keyword,next) 896 call upcase (keyword) 897 string = record(next:240) 898 if (keyword(1:17) .eq. 'POTENTIAL-SHELLS ') then 899 read (string,*,err=10,end=10) nshell 900 else if (keyword(1:18) .eq. 'POTENTIAL-SPACING ') then 901 read (string,*,err=10,end=10) spacing 902 density = 4.0d0 * pi / spacing**2 903 else if (keyword(1:17) .eq. 'POTENTIAL-FACTOR ') then 904 read (string,*,err=10,end=10) rfactor 905 else if (keyword(1:17) .eq. 'POTENTIAL-OFFSET ') then 906 read (string,*,err=10,end=10) roffset 907 end if 908 10 continue 909 end do 910c 911c perform dynamic allocation of some local arrays 912c 913 allocate (rad(n)) 914 allocate (rad2(n)) 915c 916c get modified atomic radii from consensus vdw values 917c 918 do i = 1, n 919 atn = atmnum(type(i)) 920 rad(i) = vdwrad(atn) 921 if (rad(i) .eq. 0.0d0) rad(i) = 1.7d0 922 rad(i) = rfactor*rad(i) + roffset 923 rad2(i) = rad(i) * rad(i) 924 end do 925c 926c perform dynamic allocation of some local arrays 927c 928 allocate (dot(3,maxdot)) 929c 930c find points on each of the molecular surface shells 931c 932 do m = 1, nshell 933 if (m .ne. 1) then 934 do i = 1, n 935 rad(i) = rad(i) + spacing 936 rad2(i) = rad(i) * rad(i) 937 end do 938 end if 939 do i = 1, n 940 xi = x(i) 941 yi = y(i) 942 zi = z(i) 943 ndot = int(density*rad2(i)) 944 if (ndot .gt. maxdot) then 945 write (iout,20) 946 20 format (/,' POTGRID -- Too many Surface Grid Points;', 947 & ' Increase MAXDOT') 948 call fatal 949 end if 950 call sphere (ndot,dot) 951 do j = 1, ndot 952 xj = xi + rad(i)*dot(1,j) 953 yj = yi + rad(i)*dot(2,j) 954 zj = zi + rad(i)*dot(3,j) 955 xj = dble(nint(xj/round)) * round 956 yj = dble(nint(yj/round)) * round 957 zj = dble(nint(zj/round)) * round 958 do k = 1, i-1 959 xr = xj - x(k) 960 yr = yj - y(k) 961 zr = zj - z(k) 962 r2 = xr*xr + yr*yr + zr*zr 963 if (r2 .lt. rad2(k)) goto 30 964 end do 965 do k = i+1, n 966 xr = xj - x(k) 967 yr = yj - y(k) 968 zr = zj - z(k) 969 r2 = xr*xr + yr*yr + zr*zr 970 if (r2 .lt. rad2(k)) goto 30 971 end do 972 npoint = npoint + 1 973 ipgrid(npoint,iconf) = i 974 pgrid(1,npoint,iconf) = xj 975 pgrid(2,npoint,iconf) = yj 976 pgrid(3,npoint,iconf) = zj 977 30 continue 978 end do 979 end do 980 end do 981c 982c perform deallocation of some local arrays 983c 984 deallocate (rad) 985 deallocate (rad2) 986 deallocate (dot) 987c 988c use potential grid points only for active grid atoms 989c 990 k = npoint 991 npoint = 0 992 do i = 1, k 993 if (gatm(ipgrid(i,iconf))) then 994 npoint = npoint + 1 995 ipgrid(npoint,iconf) = ipgrid(i,iconf) 996 pgrid(1,npoint,iconf) = pgrid(1,i,iconf) 997 pgrid(2,npoint,iconf) = pgrid(2,i,iconf) 998 pgrid(3,npoint,iconf) = pgrid(3,i,iconf) 999 end if 1000 end do 1001 npgrid(iconf) = npoint 1002 return 1003 end 1004c 1005c 1006c ################################################################ 1007c ## ## 1008c ## subroutine setelect -- assign electrostatic parameters ## 1009c ## ## 1010c ################################################################ 1011c 1012c 1013c "setelect" assigns partial charge, bond dipole and atomic 1014c multipole parameters for the current structure, as needed 1015c for computation of the electrostatic potential 1016c 1017c 1018 subroutine setelect 1019 implicit none 1020c 1021c 1022c get connectivity info and make parameter assignments 1023c 1024 call attach 1025 call active 1026 call bonds 1027 call angles 1028 call torsions 1029 call bitors 1030 call rings 1031 call cutoffs 1032 call katom 1033 call kcharge 1034 call kdipole 1035 call kmpole 1036 call kpolar 1037 call kchgtrn 1038 return 1039 end 1040c 1041c 1042c ################################################################# 1043c ## ## 1044c ## subroutine potpoint -- electrostatic potential at point ## 1045c ## ## 1046c ################################################################# 1047c 1048c 1049c "potpoint" calculates the electrostatic potential at a grid 1050c point "i" as the total electrostatic interaction energy of 1051c the system with a positive charge located at the grid point 1052c 1053c 1054 subroutine potpoint (xi,yi,zi,pot) 1055 use atoms 1056 use charge 1057 use chgpen 1058 use chgpot 1059 use dipole 1060 use mplpot 1061 use mpole 1062 use polar 1063 use potent 1064 use units 1065 implicit none 1066 integer k,kk,k1,k2 1067 real*8 e,ei,pot 1068 real*8 ec,ed,em,ep 1069 real*8 xi,yi,zi 1070 real*8 xk,yk,zk 1071 real*8 xr,yr,zr 1072 real*8 r,r2,dotk 1073 real*8 rk2,rkr3 1074 real*8 rr1,rr3,rr5 1075 real*8 rr1k,rr3k,rr5k 1076 real*8 corek,valk 1077 real*8 alphak 1078 real*8 f,fi,ci,ck 1079 real*8 dkx,dky,dkz 1080 real*8 ukx,uky,ukz 1081 real*8 qkxx,qkxy,qkxz 1082 real*8 qkyy,qkyz,qkzz 1083 real*8 qkx,qky,qkz 1084 real*8 dkr,qkr,ukr 1085 real*8 dmpk(5) 1086c 1087c 1088c zero out charge, dipole and multipole potential terms 1089c 1090 ec = 0.0d0 1091 ed = 0.0d0 1092 em = 0.0d0 1093 ep = 0.0d0 1094c 1095c set charge of probe site and electrostatic constants 1096c 1097 f = electric / dielec 1098 ci = 1.0d0 1099 fi = f * ci 1100c 1101c calculate the charge contribution to the potential 1102c 1103 do k = 1, nion 1104 kk = iion(k) 1105 xr = x(kk) - xi 1106 yr = y(kk) - yi 1107 zr = z(kk) - zi 1108 r2 = xr*xr + yr* yr + zr*zr 1109 r = sqrt(r2) 1110 e = fi * pchg(k) / r 1111 ec = ec + e 1112 end do 1113c 1114c calculate the bond dipole contribution to the potential 1115c 1116 do k = 1, ndipole 1117 k1 = idpl(1,k) 1118 k2 = idpl(2,k) 1119 xk = x(k2) - x(k1) 1120 yk = y(k2) - y(k1) 1121 zk = z(k2) - z(k1) 1122 xr = x(k1) + xk*sdpl(k) - xi 1123 yr = y(k1) + yk*sdpl(k) - yi 1124 zr = z(k1) + zk*sdpl(k) - zi 1125 r2 = xr*xr + yr* yr + zr*zr 1126 rk2 = xk*xk + yk*yk + zk*zk 1127 rkr3 = sqrt(rk2*r2) * r2 1128 dotk = xk*xr + yk*yr + zk*zr 1129 e = (fi/debye) * bdpl(k) * dotk / rkr3 1130 ed = ed + e 1131 end do 1132c 1133c calculate the multipole contribution to the potential 1134c 1135 do k = 1, npole 1136 kk = ipole(k) 1137 xr = x(kk) - xi 1138 yr = y(kk) - yi 1139 zr = z(kk) - zi 1140 r2 = xr*xr + yr* yr + zr*zr 1141 r = sqrt(r2) 1142 ck = rpole(1,k) 1143 dkx = rpole(2,k) 1144 dky = rpole(3,k) 1145 dkz = rpole(4,k) 1146 qkxx = rpole(5,k) 1147 qkxy = rpole(6,k) 1148 qkxz = rpole(7,k) 1149 qkyy = rpole(9,k) 1150 qkyz = rpole(10,k) 1151 qkzz = rpole(13,k) 1152 if (use_polar) then 1153 ukx = uind(1,k) 1154 uky = uind(2,k) 1155 ukz = uind(3,k) 1156 else 1157 ukx = 0.0d0 1158 uky = 0.0d0 1159 ukz = 0.0d0 1160 end if 1161c 1162c construct some common multipole and distance values 1163c 1164 qkx = qkxx*xr + qkxy*yr + qkxz*zr 1165 qky = qkxy*xr + qkyy*yr + qkyz*zr 1166 qkz = qkxz*xr + qkyz*yr + qkzz*zr 1167 dkr = dkx*xr + dky*yr + dkz*zr 1168 qkr = qkx*xr + qky*yr + qkz*zr 1169 ukr = ukx*xr + uky*yr + ukz*zr 1170 rr1 = 1.0d0 / r 1171 rr3 = rr1 / r2 1172 rr5 = 3.0d0 * rr3 / r2 1173c 1174c compute the potential contributions for this site 1175c 1176 if (use_chgpen) then 1177 corek = pcore(kk) 1178 valk = pval(kk) 1179 alphak = palpha(kk) 1180 call damppot (r,alphak,dmpk) 1181 rr1k = dmpk(1) * rr1 1182 rr3k = dmpk(3) * rr3 1183 rr5k = dmpk(5) * rr5 1184 e = corek*rr1 + valk*rr1k - dkr*rr3k + qkr*rr5k 1185 else 1186 e = ck*rr1 - dkr*rr3 + qkr*rr5 1187 end if 1188 ei = -ukr * rr3 1189c 1190c increment the overall multipole and polarization terms 1191c 1192 e = fi * e 1193 ei = fi * ei 1194 em = em + e 1195 ep = ep + ei 1196 end do 1197c 1198c potential is sum of all interactions with probe site 1199c 1200 pot = ec + ed + em + ep 1201 return 1202 end 1203c 1204c 1205c ################################################################ 1206c ## ## 1207c ## subroutine fitrsd -- residual values for potential fit ## 1208c ## ## 1209c ################################################################ 1210c 1211c 1212c "fitrsd" computes residuals for electrostatic potential fitting 1213c including total charge restraints, dipole and quadrupole moment 1214c targets, and restraints on initial parameter values 1215c 1216c 1217 subroutine fitrsd (nvar,nresid,xx,resid) 1218 use atoms 1219 use moment 1220 use neigh 1221 use potent 1222 use potfit 1223 implicit none 1224 integer i,j,nvar 1225 integer npoint 1226 integer nresid 1227 integer iresid 1228 real*8 xi,yi,zi,pot 1229 real*8 tscale,cscale 1230 real*8 rconf,ratio 1231 real*8 rscale 1232 real*8 xx(*) 1233 real*8 resid(*) 1234c 1235c 1236c initialize least squares residuals and scaling factors 1237c 1238 npoint = 0 1239 do j = 1, nconf 1240 npoint = npoint + npgrid(j) 1241 end do 1242 do j = 1, nresid 1243 resid(j) = 0.0d0 1244 end do 1245 tscale = 300.0d0 1246 cscale = 10000.0d0 1247c 1248c set weight of electrostatic potential vs. parameter restraints 1249c 1250 rconf = dble(nconf) 1251 ratio = dble(npoint) / dble(nvar*nconf) 1252c rscale = 1.0d0 * sqrt(wresp) * rconf * sqrt(ratio) 1253c rscale = 0.1d0 * sqrt(wresp) * rconf * ratio 1254 rscale = 0.006d0 * sqrt(wresp) * rconf * sqrt(ratio**3) 1255c 1256c initialize counters for parameters and residual components 1257c 1258 nvar = 0 1259 iresid = 0 1260 do j = 1, maxtyp 1261 fitchg(j) = .false. 1262 fitpol(j) = .false. 1263 end do 1264c 1265c find least squares residuals via loop over conformations 1266c 1267 do j = 1, nconf 1268 call getref (j) 1269 call setelect 1270 call varprm (nvar,xx) 1271 if (use_mpole) call rotpole 1272 if (use_polar) then 1273 domlst = .true. 1274 doulst = .true. 1275 call nblist 1276 call induce 1277 end if 1278c 1279c get residuals due to potential error over grid points 1280c 1281!$OMP PARALLEL default(private) 1282!$OMP& shared(j,npgrid,pgrid,epot,iresid,resid,rconf) 1283!$OMP DO 1284 do i = 1, npgrid(j) 1285 xi = pgrid(1,i,j) 1286 yi = pgrid(2,i,j) 1287 zi = pgrid(3,i,j) 1288 call potpoint (xi,yi,zi,pot) 1289 epot(1,i,j) = pot 1290 resid(iresid+i) = epot(1,i,j) - epot(2,i,j) 1291 end do 1292!$OMP END DO 1293!$OMP END PARALLEL 1294 iresid = iresid + npgrid(j) 1295c 1296c find moments if they contribute to the overall residue 1297c 1298 if (fit_mpl .or. use_dpl .or. use_qpl) call momfull 1299c 1300c get residual due to total molecular charge restraint 1301c 1302 if (fit_mpl) then 1303 iresid = iresid + 1 1304 resid(iresid) = (netchg-dble(nint(netchg))) * cscale 1305 end if 1306c 1307c get residuals from dipole and quadrupole target violations 1308c 1309 if (use_dpl) then 1310 resid(iresid+1) = (xdpl-xdpl0(j)) * tscale 1311 resid(iresid+2) = (ydpl-ydpl0(j)) * tscale 1312 resid(iresid+3) = (zdpl-zdpl0(j)) * tscale 1313 iresid = iresid + 3 1314 end if 1315 if (use_qpl) then 1316 resid(iresid+1) = (xxqpl-xxqpl0(j)) * tscale 1317 resid(iresid+2) = (xyqpl-xyqpl0(j)) * tscale 1318 resid(iresid+3) = (xzqpl-xzqpl0(j)) * tscale 1319 resid(iresid+4) = (yyqpl-yyqpl0(j)) * tscale 1320 resid(iresid+5) = (yzqpl-yzqpl0(j)) * tscale 1321 resid(iresid+6) = (zzqpl-zzqpl0(j)) * tscale 1322 iresid = iresid + 6 1323 end if 1324 end do 1325c 1326c get residuals due to deviation of initial parameters 1327c 1328 if (resptyp .eq. 'ORIG') then 1329 do i = 1, nvar 1330 iresid = iresid + 1 1331 resid(iresid) = (xx(i)-fit0(i)) * rscale 1332 end do 1333 else if (resptyp .eq. 'ZERO') then 1334 do i = 1, nvar 1335 iresid = iresid + 1 1336 resid(iresid) = xx(i) * rscale 1337 end do 1338 end if 1339 return 1340 end 1341c 1342c 1343c ############################################################# 1344c ## ## 1345c ## subroutine varprm -- optimization to electrostatics ## 1346c ## ## 1347c ############################################################# 1348c 1349c 1350c "varprm" copies the current optimization values into the 1351c corresponding electrostatic potential energy parameters 1352c 1353c 1354 subroutine varprm (nvar,xx) 1355 use atoms 1356 use charge 1357 use mpole 1358 use potent 1359 use potfit 1360 use units 1361 implicit none 1362 integer i,j 1363 integer ii,it 1364 integer nvar 1365 real*8 dterm,qterm 1366 real*8 xx(*) 1367 logical done 1368c 1369c 1370c translate optimization values back to partial charges 1371c 1372 do i = 1, nion 1373 done = .true. 1374 ii = iion(i) 1375 it = type(ii) 1376 if (fatm(ii)) done = .false. 1377 if (.not. done) then 1378 if (fitchg(it)) then 1379 done = .true. 1380 pchg(i) = fchg(it) 1381 end if 1382 end if 1383 if (.not. done) then 1384 if (pchg(i) .ne. 0.0d0) then 1385 nvar = nvar + 1 1386 pchg(i) = xx(nvar) 1387 end if 1388 fitchg(it) = .true. 1389 fchg(it) = pchg(i) 1390 end if 1391 end do 1392c 1393c conversion factors for dipole and quadrupole moments 1394c 1395 dterm = bohr 1396 qterm = bohr**2 / 3.0d0 1397c 1398c translate optimization values back to atomic multipoles 1399c 1400 do i = 1, npole 1401 done = .true. 1402 ii = ipole(i) 1403 it = type(ii) 1404 if (fatm(ii)) done = .false. 1405 if (.not. done) then 1406 if (fitpol(it)) then 1407 done = .true. 1408 do j = 1, 13 1409 pole(j,i) = fpol(j,it) 1410 end do 1411 end if 1412 end if 1413 if (.not. done) then 1414 if (fit_mpl .and. pole(1,i).ne.0.0d0) then 1415 nvar = nvar + 1 1416 pole(1,i) = xx(nvar) 1417 end if 1418 if (fit_dpl .and. pole(2,i).ne.0.0d0) then 1419 nvar = nvar + 1 1420 pole(2,i) = dterm * xx(nvar) 1421 end if 1422 if (fit_dpl .and. pole(3,i).ne.0.0d0) then 1423 nvar = nvar + 1 1424 pole(3,i) = dterm * xx(nvar) 1425 end if 1426 if (fit_dpl .and. pole(4,i).ne.0.0d0) then 1427 nvar = nvar + 1 1428 pole(4,i) = dterm * xx(nvar) 1429 end if 1430 if (fit_qpl .and. pole(5,i).ne.0.0d0) then 1431 if (polaxe(i) .ne. 'Z-Only') then 1432 nvar = nvar + 1 1433 pole(5,i) = qterm * xx(nvar) 1434 end if 1435 end if 1436 if (fit_qpl .and. pole(6,i).ne.0.0d0) then 1437 nvar = nvar + 1 1438 pole(6,i) = qterm * xx(nvar) 1439 pole(8,i) = qterm * xx(nvar) 1440 end if 1441 if (fit_qpl .and. pole(7,i).ne.0.0d0) then 1442 nvar = nvar + 1 1443 pole(7,i) = qterm * xx(nvar) 1444 pole(11,i) = qterm * xx(nvar) 1445 end if 1446 if (fit_qpl .and. pole(9,i).ne.0.0d0) then 1447 if (polaxe(i) .ne. 'Z-Only') then 1448 nvar = nvar + 1 1449 pole(9,i) = qterm * xx(nvar) 1450 end if 1451 end if 1452 if (fit_qpl .and. pole(10,i).ne.0.0d0) then 1453 nvar = nvar + 1 1454 pole(10,i) = qterm * xx(nvar) 1455 pole(12,i) = qterm * xx(nvar) 1456 end if 1457 if (fit_qpl .and. pole(13,i).ne.0.0d0) then 1458 if (polaxe(i) .eq. 'Z-Only') then 1459 nvar = nvar + 1 1460 pole(13,i) = qterm * xx(nvar) 1461 pole(5,i) = -0.5d0 * pole(13,i) 1462 pole(9,i) = pole(5,i) 1463 else 1464 pole(13,i) = -pole(5,i) - pole(9,i) 1465 end if 1466 end if 1467 fitpol(it) = .true. 1468 do j = 1, 13 1469 fpol(j,it) = pole(j,i) 1470 end do 1471 end if 1472 end do 1473c 1474c check chiral multipoles and rotate into global frame 1475c 1476 if (use_mpole) then 1477 call chkpole 1478 call rotpole 1479 end if 1480 return 1481 end 1482c 1483c 1484c ############################################################# 1485c ## ## 1486c ## subroutine prmvar -- electrostatics to optimization ## 1487c ## ## 1488c ############################################################# 1489c 1490c 1491c "prmvar" determines the optimization values from the 1492c corresponding electrostatic potential energy parameters 1493c 1494c 1495 subroutine prmvar (nvar,xx) 1496 use atomid 1497 use atoms 1498 use charge 1499 use iounit 1500 use mpole 1501 use potfit 1502 use units 1503 implicit none 1504 integer i,j,k,m 1505 integer ii,it,kt 1506 integer ktype 1507 integer nvar 1508 integer, allocatable :: equiv(:) 1509 real*8 dterm,qterm 1510 real*8 eps,sum,big 1511 real*8 ival,kval 1512 real*8 xx(*) 1513 logical done 1514 character*17 prmtyp 1515c 1516c 1517c conversion factors for dipole and quadrupole moments 1518c 1519 dterm = 1.0d0 / bohr 1520 qterm = 3.0d0 / bohr**2 1521c 1522c regularize charges, monopoles and diagonal quadrupoles 1523c 1524 eps = 0.00001d0 1525 do i = 1, nion 1526 pchg(i) = dble(nint(pchg(i)/eps)) * eps 1527 end do 1528 do i = 1, npole 1529 pole(1,i) = dble(nint(pole(1,i)/eps)) * eps 1530 pole(5,i) = dble(nint(qterm*pole(5,i)/eps)) * eps/qterm 1531 pole(9,i) = dble(nint(qterm*pole(9,i)/eps)) * eps/qterm 1532 pole(13,i) = dble(nint(qterm*pole(13,i)/eps)) * eps/qterm 1533 end do 1534c 1535c perform dynamic allocation of some local arrays 1536c 1537 allocate (equiv(maxtyp)) 1538c 1539c enforce integer net charge over partial charges 1540c 1541 do i = 1, maxtyp 1542 equiv(i) = 0 1543 end do 1544 ktype = 0 1545 sum = 0.0d0 1546 do i = 1, nion 1547 it = type(iion(i)) 1548 equiv(it) = equiv(it) + 1 1549 sum = sum + pchg(i) 1550 end do 1551 sum = sum - dble(nint(sum)) 1552 k = nint(abs(sum)/eps) 1553 do j = 1, k 1554 m = k / j 1555 if (k .eq. m*j) then 1556 do i = 1, nion 1557 it = type(iion(i)) 1558 if (equiv(it) .eq. m) then 1559 ival = abs(pchg(i)) 1560 if (ktype .eq. 0) then 1561 ktype = it 1562 kval = ival 1563 else if (ival .gt. kval) then 1564 ktype = it 1565 kval = ival 1566 end if 1567 end if 1568 end do 1569 end if 1570 if (ktype .ne. 0) goto 10 1571 end do 1572 10 continue 1573 if (ktype .ne. 0) then 1574 sum = sum / dble(m) 1575 do i = 1, nion 1576 it = type(iion(i)) 1577 if (it .eq. ktype) then 1578 pchg(i) = pchg(i) - sum 1579 fchg(it) = pchg(i) 1580 end if 1581 end do 1582 end if 1583c 1584c enforce integer net charge over atomic multipoles 1585c 1586 do i = 1, maxtyp 1587 equiv(i) = 0 1588 end do 1589 ktype = 0 1590 sum = 0.0d0 1591 do i = 1, npole 1592 it = type(ipole(i)) 1593 equiv(it) = equiv(it) + 1 1594 sum = sum + pole(1,i) 1595 end do 1596 sum = sum - dble(nint(sum)) 1597 k = nint(abs(sum)/eps) 1598 do j = 1, k 1599 m = k / j 1600 if (k .eq. m*j) then 1601 do i = 1, npole 1602 it = type(ipole(i)) 1603 if (equiv(it) .eq. m) then 1604 ival = abs(pole(1,ipole(i))) 1605 if (ktype .eq. 0) then 1606 ktype = it 1607 kval = ival 1608 else if (ival .gt. kval) then 1609 ktype = it 1610 kval = ival 1611 end if 1612 end if 1613 end do 1614 end if 1615 if (ktype .ne. 0) goto 20 1616 end do 1617 20 continue 1618 if (ktype .ne. 0) then 1619 sum = sum / dble(m) 1620 do i = 1, npole 1621 it = type(ipole(i)) 1622 if (it .eq. ktype) then 1623 pole(1,i) = pole(1,i) - sum 1624 fpol(1,it) = pole(1,i) 1625 end if 1626 end do 1627 end if 1628c 1629c perform deallocation of some local arrays 1630c 1631 deallocate (equiv) 1632c 1633c enforce traceless quadrupole at each multipole site 1634c 1635 do i = 1, npole 1636 sum = pole(5,i) + pole(9,i) + pole(13,i) 1637 big = max(abs(pole(5,i)),abs(pole(9,i)),abs(pole(13,i))) 1638 k = 0 1639 if (big .eq. abs(pole(5,i))) k = 5 1640 if (big .eq. abs(pole(9,i))) k = 9 1641 if (big .eq. abs(pole(13,i))) k = 13 1642 if (k .ne. 0) then 1643 it = type(ipole(i)) 1644 pole(k,i) = pole(k,i) - sum 1645 fpol(k,it) = pole(k,i) 1646 end if 1647 end do 1648c 1649c list active atoms when not all are used in optimization 1650c 1651 if (nconf.eq.1 .and. nfatm.ne.n) then 1652 write (iout,30) 1653 30 format (/,' Atomic Parameters Included in Potential Fitting :', 1654 & //,3x,'Atom',10x,'Atom Name',9x,'Atom Type', 1655 & 9x,'Parameters',/) 1656 do i = 1, nion 1657 ii = iion(i) 1658 if (fatm(ii)) then 1659 it = type(ii) 1660 prmtyp = 'Partial Charge' 1661 write (iout,40) ii,name(ii),it,prmtyp 1662 40 format (i6,15x,a3,10x,i6,13x,a) 1663 end if 1664 end do 1665 do i = 1, npole 1666 ii = ipole(i) 1667 if (fatm(ii)) then 1668 it = type(ii) 1669 prmtyp = 'Atomic Multipoles' 1670 write (iout,50) ii,name(ii),it,prmtyp 1671 50 format (i6,15x,a3,10x,i6,13x,a) 1672 end if 1673 end do 1674 end if 1675c 1676c print header information for electrostatic parameters 1677c 1678 if (nvar .eq. 0) then 1679 write (iout,60) 1680 60 format (/,' Potential Fitting of Electrostatic Parameters :', 1681 & //,1x,'Parameter',6x,'Atom Type',9x,'Category', 1682 & 12x,'Value',9x,'Fixed',/) 1683 end if 1684c 1685c get optimization parameters from partial charge values 1686c 1687 do i = 1, nion 1688 done = .true. 1689 ii = iion(i) 1690 it = type(ii) 1691 if (fatm(ii)) done = .false. 1692 if (.not. done) then 1693 if (fitchg(it)) done = .true. 1694 fitchg(it) = .true. 1695 end if 1696 if (.not. done) then 1697 if (pchg(i) .ne. 0.0d0) then 1698 nvar = nvar + 1 1699 xx(nvar) = pchg(i) 1700 write (iout,70) nvar,it,'Charge ',xx(nvar) 1701 70 format (i6,7x,i8,13x,a8,6x,f12.5) 1702 else 1703 write (iout,80) it,'Charge ',pchg(i) 1704 80 format (4x,'--',7x,i8,13x,a8,6x,f12.5,10x,'X') 1705 end if 1706 end if 1707 end do 1708c 1709c get optimization parameters from atomic multipole values 1710c 1711 do i = 1, npole 1712 done = .true. 1713 ii = ipole(i) 1714 it = type(ii) 1715 if (fatm(ii)) done = .false. 1716 if (.not. done) then 1717 if (fitpol(it)) done = .true. 1718 fitpol(it) = .true. 1719 end if 1720 if (.not. done) then 1721 if (fit_mpl .and. pole(1,i).ne.0.0d0) then 1722 nvar = nvar + 1 1723 xx(nvar) = pole(1,i) 1724 write (iout,90) nvar,it,'Monopole',xx(nvar) 1725 90 format (i6,7x,i8,13x,a8,6x,f12.5) 1726 else 1727 write (iout,100) it,'Monopole',pole(1,i) 1728 100 format (4x,'--',7x,i8,13x,a8,6x,f12.5,10x,'X') 1729 end if 1730 if (fit_dpl .and. pole(2,i).ne.0.0d0) then 1731 nvar = nvar + 1 1732 xx(nvar) = dterm * pole(2,i) 1733 write (iout,110) nvar,it,'X-Dipole',xx(nvar) 1734 110 format (i6,7x,i8,13x,a8,6x,f12.5) 1735 else 1736 write (iout,120) it,'X-Dipole',dterm*pole(2,i) 1737 120 format (4x,'--',7x,i8,13x,a8,6x,f12.5,10x,'X') 1738 end if 1739 if (fit_dpl .and. pole(3,i).ne.0.0d0) then 1740 nvar = nvar + 1 1741 xx(nvar) = dterm * pole(3,i) 1742 write (iout,130) nvar,it,'Y-Dipole',xx(nvar) 1743 130 format (i6,7x,i8,13x,a8,6x,f12.5) 1744 else 1745 write (iout,140) it,'Y-Dipole',dterm*pole(3,i) 1746 140 format (4x,'--',7x,i8,13x,a8,6x,f12.5,10x,'X') 1747 end if 1748 if (fit_dpl .and. pole(4,i).ne.0.0d0) then 1749 nvar = nvar + 1 1750 xx(nvar) = dterm * pole(4,i) 1751 write (iout,150) nvar,it,'Z-Dipole',xx(nvar) 1752 150 format (i6,7x,i8,13x,a8,6x,f12.5) 1753 else 1754 write (iout,160) it,'Z-Dipole',dterm*pole(4,i) 1755 160 format (4x,'--',7x,i8,13x,a8,6x,f12.5,10x,'X') 1756 end if 1757 if (fit_qpl .and. pole(5,i).ne.0.0d0) then 1758 if (polaxe(i) .ne. 'Z-Only') then 1759 nvar = nvar + 1 1760 xx(nvar) = qterm * pole(5,i) 1761 write (iout,170) nvar,it,'XX-Quad ',xx(nvar) 1762 170 format (i6,7x,i8,13x,a8,6x,f12.5) 1763 else 1764 write (iout,180) it,'XX-Quad ',qterm*pole(5,i) 1765 180 format (4x,'--',7x,i8,13x,a8,6x,f12.5) 1766 end if 1767 else 1768 write (iout,190) it,'XX-Quad ',qterm*pole(5,i) 1769 190 format (4x,'--',7x,i8,13x,a8,6x,f12.5,10x,'X') 1770 end if 1771 if (fit_qpl .and. pole(6,i).ne.0.0d0) then 1772 nvar = nvar + 1 1773 xx(nvar) = qterm * pole(6,i) 1774 write (iout,200) nvar,it,'XY-Quad ',xx(nvar) 1775 200 format (i6,7x,i8,13x,a8,6x,f12.5) 1776 else 1777 write (iout,210) it,'XY-Quad ',qterm*pole(6,i) 1778 210 format (4x,'--',7x,i8,13x,a8,6x,f12.5,10x,'X') 1779 end if 1780 if (fit_qpl .and. pole(7,i).ne.0.0d0) then 1781 nvar = nvar + 1 1782 xx(nvar) = qterm * pole(7,i) 1783 write (iout,220) nvar,it,'XZ-Quad ',xx(nvar) 1784 220 format (i6,7x,i8,13x,a8,6x,f12.5) 1785 else 1786 write (iout,230) it,'XZ-Quad ',qterm*pole(7,i) 1787 230 format (4x,'--',7x,i8,13x,a8,6x,f12.5,10x,'X') 1788 end if 1789 if (fit_qpl .and. pole(9,i).ne.0.0d0) then 1790 if (polaxe(i) .ne. 'Z-Only') then 1791 nvar = nvar + 1 1792 xx(nvar) = qterm * pole(9,i) 1793 write (iout,240) nvar,it,'YY-Quad ',xx(nvar) 1794 240 format (i6,7x,i8,13x,a8,6x,f12.5) 1795 else 1796 write (iout,250) it,'YY-Quad ',qterm*pole(9,i) 1797 250 format (4x,'--',7x,i8,13x,a8,6x,f12.5) 1798 end if 1799 else 1800 write (iout,260) it,'YY-Quad ',qterm*pole(9,i) 1801 260 format (4x,'--',7x,i8,13x,a8,6x,f12.5,10x,'X') 1802 end if 1803 if (fit_qpl .and. pole(10,i).ne.0.0d0) then 1804 nvar = nvar + 1 1805 xx(nvar) = qterm * pole(10,i) 1806 write (iout,270) nvar,it,'YZ-Quad ',xx(nvar) 1807 270 format (i6,7x,i8,13x,a8,6x,f12.5) 1808 else 1809 write (iout,280) it,'YZ-Quad ',qterm*pole(10,i) 1810 280 format (4x,'--',7x,i8,13x,a8,6x,f12.5,10x,'X') 1811 end if 1812 if (fit_qpl .and. pole(13,i).ne.0.0d0) then 1813 if (polaxe(i) .eq. 'Z-Only') then 1814 nvar = nvar + 1 1815 xx(nvar) = qterm * pole(13,i) 1816 write (iout,290) nvar,it,'ZZ-Quad ',xx(nvar) 1817 290 format (i6,7x,i8,13x,a8,6x,f12.5) 1818 else 1819 write (iout,300) it,'ZZ-Quad ',qterm*pole(13,i) 1820 300 format (4x,'--',7x,i8,13x,a8,6x,f12.5) 1821 end if 1822 else 1823 write (iout,310) it,'ZZ-Quad ',qterm*pole(13,i) 1824 310 format (4x,'--',7x,i8,13x,a8,6x,f12.5,10x,'X') 1825 end if 1826 end if 1827 end do 1828 return 1829 end 1830c 1831c 1832c ################################################################## 1833c ## ## 1834c ## subroutine potstat -- electrostatic potential statistics ## 1835c ## ## 1836c ################################################################## 1837c 1838c 1839c "potstat" computes and prints statistics for the electrostatic 1840c potential over a set of grid points 1841c 1842c 1843 subroutine potstat (dofull,domodel,dopair,dotarget) 1844 use atoms 1845 use files 1846 use iounit 1847 use potfit 1848 use refer 1849 use titles 1850 implicit none 1851 integer i,j,k 1852 integer ipot,npoint 1853 integer freeunit 1854 integer trimtext 1855 integer, allocatable :: natm(:) 1856 real*8 xi,yi,zi 1857 real*8 pave1,pave2 1858 real*8 mave1,mave2 1859 real*8 tave,uave,rmsd 1860 real*8, allocatable :: patm1(:) 1861 real*8, allocatable :: patm2(:) 1862 real*8, allocatable :: rmsa(:) 1863 logical dofull,domodel 1864 logical dopair,dotarget 1865 character*240 potfile 1866c 1867c 1868c output potential values for each model at each point 1869c 1870 if (dofull) then 1871 if (domodel) then 1872 ipot = freeunit () 1873 potfile = filename(1:leng)//'.pot' 1874 call version (potfile,'new') 1875 open (unit=ipot,file=potfile,status='new') 1876 end if 1877 do j = 1, nconf 1878 if (nconf .eq. 1) then 1879 write (iout,10) 1880 10 format (/,' Electrostatic Potential at Each Grid', 1881 & ' Point :', 1882 & /,8x,'(Kcal/mole per unit charge)') 1883 else 1884 write (iout,20) j 1885 20 format (/,' Electrostatic Potential at Grid Points', 1886 & ' for Structure',i4,' :', 1887 & /,12x,'(Kcal/mole per unit charge)') 1888 end if 1889 if (dotarget) then 1890 write (iout,30) 1891 30 format (/,3x,'Point',15x,'XYZ-Coordinates',15x, 1892 & 'Potential',5x,'Target',/) 1893 else if (dopair) then 1894 write (iout,40) 1895 40 format (/,3x,'Point',15x,'XYZ-Coordinates',13x, 1896 & 'Potential 1',3x,'Potential 2',/) 1897 else if (domodel) then 1898 write (iout,50) 1899 50 format (/,3x,'Point',15x,'XYZ-Coordinates',14x, 1900 & 'Potential',/) 1901 write (ipot,60) npgrid(j),title(1:ltitle) 1902 60 format (i8,2x,a) 1903 end if 1904 do i = 1, npgrid(j) 1905 xi = pgrid(1,i,j) 1906 yi = pgrid(2,i,j) 1907 zi = pgrid(3,i,j) 1908 if (dotarget .or. dopair) then 1909 write (iout,70) i,xi,yi,zi,epot(1,i,j),epot(2,i,j) 1910 70 format (i8,3x,3f12.6,2x,2f12.4) 1911 else if (domodel) then 1912 write (iout,80) i,xi,yi,zi,epot(1,i,j) 1913 80 format (i8,3x,3f12.6,2x,f12.4) 1914 write (ipot,90) i,xi,yi,zi,epot(1,i,j) 1915 90 format (i8,3x,3f12.6,2x,f12.4) 1916 end if 1917 end do 1918 end do 1919 if (domodel) then 1920 close (unit=ipot) 1921 write (iout,100) potfile(1:trimtext(potfile)) 1922 100 format (/,' Electrostatic Potential Written To : ',a) 1923 end if 1924 end if 1925c 1926c perform dynamic allocation of some local arrays 1927c 1928 allocate (natm(namax)) 1929 allocate (patm1(namax)) 1930 allocate (patm2(namax)) 1931 allocate (rmsa(namax)) 1932c 1933c find average electrostatic potential around each atom 1934c 1935 write (iout,110) 1936 110 format (/,' Average Electrostatic Potential over Atoms :', 1937 & /,6x,'(Kcal/mole per unit charge)') 1938 if (dotarget) then 1939 write (iout,120) 1940 120 format (/,3x,'Structure',3x,'Atom',6x,'Points', 1941 & 6x,'Potential',8x,'Target',8x,'RMS Diff',/) 1942 else if (dopair) then 1943 write (iout,130) 1944 130 format (/,3x,'Structure',3x,'Atom',6x,'Points', 1945 & 5x,'Potential 1',4x,'Potential 2',6x,'RMS Diff',/) 1946 else if (domodel) then 1947 write (iout,140) 1948 140 format (/,3x,'Structure',3x,'Atom',5x,'Points', 1949 & 6x,'Potential',/) 1950 end if 1951 do j = 1, nconf 1952 call getref (j) 1953 do i = 1, n 1954 natm(i) = 0 1955 patm1(i) = 0.0d0 1956 patm2(i) = 0.0d0 1957 rmsa(i) = 0.0d0 1958 end do 1959 do i = 1, npgrid(j) 1960 k = ipgrid(i,j) 1961 natm(k) = natm(k) + 1 1962 patm1(k) = patm1(k) + epot(1,i,j) 1963 patm2(k) = patm2(k) + epot(2,i,j) 1964 rmsa(k) = rmsa(k) + (epot(1,i,j)-epot(2,i,j))**2 1965 end do 1966 do i = 1, n 1967 if (natm(i) .ne. 0) then 1968 patm1(i) = patm1(i) / dble(natm(i)) 1969 patm2(i) = patm2(i) / dble(natm(i)) 1970 rmsa(i) = sqrt(rmsa(i)/dble(natm(i))) 1971 end if 1972 if (gatm(i)) then 1973 if (dotarget .or. dopair) then 1974 write (iout,150) j,i,natm(i),patm1(i), 1975 & patm2(i),rmsa(i) 1976 150 format (2i9,3x,i9,3x,f12.4,3x,f12.4,3x,f12.4) 1977 else if (domodel) then 1978 write (iout,160) j,i,natm(i),patm1(i) 1979 160 format (2i9,3x,i9,3x,f12.4) 1980 end if 1981 end if 1982 end do 1983 end do 1984c 1985c perform deallocation of some local arrays 1986c 1987 deallocate (natm) 1988 deallocate (patm1) 1989 deallocate (patm2) 1990 deallocate (rmsa) 1991c 1992c overall averages for the sets of electrostatic potentials 1993c 1994 npoint = 0 1995 pave1 = 0.0d0 1996 pave2 = 0.0d0 1997 mave1 = 0.0d0 1998 mave2 = 0.0d0 1999 tave = 0.0d0 2000 uave = 0.0d0 2001 rmsd = 0.0d0 2002 do j = 1, nconf 2003 npoint = npoint + npgrid(j) 2004 do i = 1, npgrid(j) 2005 pave1 = pave1 + epot(1,i,j) 2006 pave2 = pave2 + epot(2,i,j) 2007 mave1 = mave1 + abs(epot(1,i,j)) 2008 mave2 = mave2 + abs(epot(2,i,j)) 2009 tave = tave + epot(1,i,j) - epot(2,i,j) 2010 uave = uave + abs(epot(1,i,j)-epot(2,i,j)) 2011 rmsd = rmsd + (epot(1,i,j)-epot(2,i,j))**2 2012 end do 2013 end do 2014 pave1 = pave1 / dble(npoint) 2015 pave2 = pave2 / dble(npoint) 2016 mave1 = mave1 / dble(npoint) 2017 mave2 = mave2 / dble(npoint) 2018 tave = tave / dble(npoint) 2019 uave = uave / dble(npoint) 2020 rmsd = sqrt(rmsd/dble(npoint)) 2021 if (dopair) then 2022 write (iout,170) pave1,mave1 2023 170 format (/,' Electrostatic Potential over all Grid Points :', 2024 & //,' Average Potential Value for Model 1 :',10x,f12.4, 2025 & /,' Average Potential Magnitude for Model 1 :', 2026 & 6x,f12.4) 2027 else 2028 write (iout,180) pave1,mave1 2029 180 format (/,' Electrostatic Potential over all Grid Points :', 2030 & //,' Average Potential Value for Model :',12x,f12.4, 2031 & /,' Average Potential Magnitude for Model :',8x,f12.4) 2032 end if 2033 if (dotarget) then 2034 write (iout,190) pave2,mave2,tave,uave,rmsd 2035 190 format (' Average Potential Value for Target :',11x,f12.4, 2036 & /,' Average Potential Magnitude for Target :',7x,f12.4, 2037 & //,' Average Signed Potential Difference :',10x,f12.4, 2038 & /,' Average Unsigned Potential Difference :',8x,f12.4, 2039 & /,' Root Mean Square Potential Difference :',8x,f12.4) 2040 else if (dopair) then 2041 write (iout,200) pave2,mave2,tave,uave,rmsd 2042 200 format (' Average Potential Value for Model 2 :',10x,f12.4, 2043 & /,' Average Potential Magnitude for Model 2 :', 2044 & 6x,f12.4, 2045 & //,' Average Signed Potential Difference :',10x,f12.4, 2046 & /,' Average Unsigned Potential Difference :',8x,f12.4, 2047 & /,' Root Mean Square Potential Difference :',8x,f12.4) 2048 end if 2049 return 2050 end 2051c 2052c 2053c ################################################################## 2054c ## ## 2055c ## subroutine prtfit -- create file with optimal parameters ## 2056c ## ## 2057c ################################################################## 2058c 2059c 2060c "prtfit" makes a key file containing results from fitting a 2061c charge or multipole model to an electrostatic potential grid 2062c 2063c 2064 subroutine prtfit 2065 use atoms 2066 use charge 2067 use files 2068 use keys 2069 use mpole 2070 use potfit 2071 use units 2072 implicit none 2073 integer i,j,k,m 2074 integer ii,it 2075 integer ix,iy,iz 2076 integer ikey,size 2077 integer ntot 2078 integer freeunit 2079 integer trimtext 2080 real*8 dterm,qterm 2081 logical done,header 2082 character*4 pa,pb,pc,pd 2083 character*16, allocatable :: pt(:) 2084 character*240 keyfile 2085 character*240 record 2086c 2087c 2088c reread the contents of any previously existing keyfile 2089c 2090 call getkey 2091c 2092c open a new keyfile to contain the optimized parameters 2093c 2094 ikey = freeunit () 2095 keyfile = filename(1:leng)//'.key' 2096 call version (keyfile,'new') 2097 open (unit=ikey,file=keyfile,status='new') 2098c 2099c copy the contents of any previously existing keyfile 2100c 2101 do i = 1, nkey 2102 record = keyline(i) 2103 size = trimtext (record) 2104 write (ikey,10) record(1:size) 2105 10 format (a) 2106 end do 2107c 2108c zero the keyfile length to avoid parameter reprocessing 2109c 2110 nkey = 0 2111c 2112c output optimized partial charge values to the keyfile 2113c 2114 header = .true. 2115 do i = 1, maxtyp 2116 fitchg(i) = .false. 2117 end do 2118 do k = 1, nconf 2119 call getref (k) 2120 call setelect 2121 do i = 1, nion 2122 done = .true. 2123 ii = iion(i) 2124 it = type(ii) 2125 if (fatm(ii)) done = .false. 2126 if (.not. done) then 2127 if (fitchg(it)) done = .true. 2128 fitchg(it) = .true. 2129 end if 2130 if (.not. done) then 2131 pchg(i) = fchg(it) 2132 if (header) then 2133 header = .false. 2134 write (ikey,20) 2135 20 format (/,'#',/,'# Charges from Electrostatic', 2136 & ' Potential Fitting',/,'#',/) 2137 end if 2138 write (ikey,30) it,pchg(i) 2139 30 format ('charge',4x,i5,10x,f11.4) 2140 end if 2141 end do 2142 end do 2143c 2144c conversion factors for dipole and quadrupole moments 2145c 2146 dterm = 1.0d0 / bohr 2147 qterm = 3.0d0 / bohr**2 2148c 2149c get total atoms in all structures used in the fitting 2150c 2151 ntot = 0 2152 do k = 1, nconf 2153 call getref(k) 2154 ntot = ntot + n 2155 end do 2156c 2157c perform dynamic allocation of some local arrays 2158c 2159 allocate (pt(ntot)) 2160c 2161c initialize atom type and local frame defining strings 2162c 2163 do i = 1, ntot 2164 pt(i) = ' ' 2165 end do 2166c 2167c output optimized atomic multipole values to the keyfile 2168c 2169 header = .true. 2170 m = 0 2171 do k = 1, nconf 2172 call getref (k) 2173 call setelect 2174 do i = 1, npole 2175 done = .true. 2176 ii = ipole(i) 2177 it = type(ii) 2178 if (fatm(ii)) done = .false. 2179 if (.not. done) then 2180 iz = zaxis(i) 2181 ix = xaxis(i) 2182 iy = yaxis(i) 2183 if (iz .ne. 0) iz = type(iz) 2184 if (ix .ne. 0) ix = type(ix) 2185 if (iy .ne. 0) iy = type(abs(iy)) 2186 size = 4 2187 call numeral (it,pa,size) 2188 call numeral (iz,pb,size) 2189 call numeral (ix,pc,size) 2190 call numeral (iy,pd,size) 2191 m = m + 1 2192 pt(m) = pa//pb//pc//pd 2193 do j = 1, m-1 2194 if (pt(m) .eq. pt(j)) then 2195 done = .true. 2196 goto 40 2197 end if 2198 end do 2199 40 continue 2200 end if 2201 if (.not. done) then 2202 if (header) then 2203 header = .false. 2204 write (ikey,50) 2205 50 format (/,'#',/,'# Multipoles from Electrostatic', 2206 & ' Potential Fitting',/,'#',/) 2207 end if 2208 pole(1,i) = fpol(1,it) 2209 do j = 2, 4 2210 pole(j,i) = dterm * fpol(j,it) 2211 end do 2212 do j = 5, 13 2213 pole(j,i) = qterm * fpol(j,it) 2214 end do 2215 if (polaxe(i) .eq. 'None') then 2216 write (ikey,60) it,pole(1,i) 2217 60 format ('multipole',1x,i5,21x,f11.5) 2218 else if (polaxe(i) .eq. 'Z-Only') then 2219 write (ikey,70) it,iz,pole(1,i) 2220 70 format ('multipole',1x,2i5,16x,f11.5) 2221 else if (polaxe(i) .eq. 'Z-then-X') then 2222 if (yaxis(i) .eq. 0) then 2223 write (ikey,80) it,iz,ix,pole(1,i) 2224 80 format ('multipole',1x,3i5,11x,f11.5) 2225 else 2226 if (yaxis(i) .lt. 0) then 2227 pole(3,i) = -pole(3,i) 2228 pole(6,i) = -pole(6,i) 2229 pole(8,i) = -pole(8,i) 2230 pole(10,i) = -pole(10,i) 2231 pole(12,i) = -pole(12,i) 2232 end if 2233 write (ikey,90) it,iz,ix,iy,pole(1,i) 2234 90 format ('multipole',1x,4i5,6x,f11.5) 2235 end if 2236 else if (polaxe(i) .eq. 'Bisector') then 2237 if (yaxis(i) .eq. 0) then 2238 write (ikey,100) it,-iz,-ix,pole(1,i) 2239 100 format ('multipole',1x,3i5,11x,f11.5) 2240 else 2241 write (ikey,110) it,-iz,-ix,iy,pole(1,i) 2242 110 format ('multipole',1x,4i5,6x,f11.5) 2243 end if 2244 else if (polaxe(i) .eq. 'Z-Bisect') then 2245 write (ikey,120) it,iz,-ix,-iy,pole(1,i) 2246 120 format ('multipole',1x,4i5,6x,f11.5) 2247 else if (polaxe(i) .eq. '3-Fold') then 2248 write (ikey,130) it,-iz,-ix,-iy,pole(1,i) 2249 130 format ('multipole',1x,4i5,6x,f11.5) 2250 end if 2251 write (ikey,140) pole(2,i),pole(3,i),pole(4,i) 2252 140 format (36x,3f11.5) 2253 write (ikey,150) pole(5,i) 2254 150 format (36x,f11.5) 2255 write (ikey,160) pole(8,i),pole(9,i) 2256 160 format (36x,2f11.5) 2257 write (ikey,170) pole(11,i),pole(12,i),pole(13,i) 2258 170 format (36x,3f11.5) 2259 end if 2260 end do 2261 end do 2262 close (unit=ikey) 2263c 2264c perform deallocation of some local arrays 2265c 2266 deallocate (pt) 2267 return 2268 end 2269c 2270c 2271c ########################################################### 2272c ## ## 2273c ## subroutine potwrt -- least squares output routine ## 2274c ## ## 2275c ########################################################### 2276c 2277c 2278 subroutine potwrt (niter,nresid,xx,gs,resid) 2279 implicit none 2280 integer niter 2281 integer nresid 2282 real*8 xx(*) 2283 real*8 gs(*) 2284 real*8 resid(*) 2285c 2286c 2287c information to be printed at each least squares iteration 2288c 2289 return 2290 end 2291