1 subroutine eem(iop,iasel,istat) 2c 3c program for the calculation of eem parameters 4c patrick bultinck, jrf, 2001 5c 6c J. Phys. Chem. A (2002) 7c 8 implicit double precision(a-h,o-z) 9 parameter (mxeat=300) 10 dimension var(mxeat,2),ipntr(mxeat) 11 12 istat = 0 13 numat = 0 14 call valdis(var,ipntr,numat,iop,iasel,istat) 15 if (istat.eq.0) call eemcalc(var,ipntr,numat) 16 17 18 return 19 20 end 21 22 subroutine valdid(var,ipntr,numat,iop,iasel,istat, 23 & ianz,iresid,qat) 24 25 implicit double precision(a-h,o-z), integer (i-n) 26 parameter (mxeat=300) 27 parameter (maxtyp=19) 28 parameter (mxtyp4=4*maxtyp) 29 parameter (mxel=100) 30 parameter (mexcl=10) 31 32 common /eempar/param(maxtyp,2),pparam(maxtyp,2), 33 & peparm(10,2),esffprm(19,2),ipt(mxel,4) 34 common /athlp/ iatoms, mxnat 35 character*2 elemnt 36 common /elem/elemnt(mxel) 37 common /charge/ dipo(3),ihasq,ihsdp,iqon,idipon 38 common /metexc/ qexcl(mexcl),ianexc(mexcl) 39 dimension ianz(*),iresid(*),qat(*),var(mxeat,2),ipntr(mxeat) 40 41c eem NPA 42 43 data ((param(i,k),k=1,2),i=1,maxtyp) / 44 & -0.00298,0.89729,0.54541,0.00000,0.93882,0.00000, 45 & 0.00000,0.10570,0.00000,0.61634,0.35970,0.33639, 46 & 0.53938,0.37574,1.04219,0.72138,1.44195,1.57956, 47 & 0.28351,0.00000,0.00000,0.00000,0.00000,0.00000, 48 & 0.00000,0.00000,0.00000,0.00000,0.00000,0.00000, 49 & 0.00000,0.00000,0.00000,0.00000,0.00000,0.00000, 50 & 0.00000,0.00000/ 51c eem Mull + P en CL estimated 52 data ((pparam(i,k),k=1,2),i=1,maxtyp) / 53 & 0.20606,0.65971,0.00000,0.00000,0.00000,0.75331, 54 & 0.00000,0.00000,0.00000,0.00000,0.36237,0.32966, 55 & 0.49279,0.34519,0.73013,0.54428,0.72052,0.72664, 56 & 0.78058,0.75467,0.00000,0.00000,0.00000,0.00000, 57 & 0.00000,0.00000,0.00000,0.00000,0.36200,0.32000, 58 & 0.00000,0.00000,0.41200,0.36600,0.00000,0.00000, 59 & 0.00000,0.00000/ 60c 61c Pesp, eem parameters 62c khi eta atom MMtype Definition 63c 0.68903 0.39888 1 H common hydrogen 64c 0.71998 0.34395 6 C common carbon 65c 0.97104 0.46763 7 N common nitrogen 66c 0.95419 0.43623 8 O common oxigen 67c 0.78155 0.24617 16 S common sulphur 68c 0.67920 0.17446 12 Mg common magnesium 69c 0.98387 0.13229 15 P common phosporous 70c 0.86408 0.40514 9 F common fluorine 71c 0.76957 0.17213 17 Cl common chlorine 72c 0.77709 0.15474 35 Br common bromine 73c 0.31627 0.90661 1 HAc acidic proton 74c 0.58475 0.54282 1 Hpo H bound to electrophile 75c 0.61286 0.82238 1 Hpm H bound to weak electrophile 76c 0.67301 0.72481 1 HC H bound to alkane carbon 77c 0.68654 0.67974 1 HC3 H in -CH3 78c 0.54822 0.44969 6 CO carbonil carbon 79c 0.79906 0.33085 6 CAr aromatic carbon 80c 0.75858 0.28429 6 CArp polarized aromatic carbon 81c 0.72259 0.28054 6 CT sp3 carbon 82c 0.84720 0.70094 6 CT3 carbon in -CH3 83c 0.78794 0.30423 6 CX sp2 carbon 84c 0.85476 0.38483 6 CZ sp carbon 85c 0.91123 0.40572 7 NT sp3 nitrogen 86c 0.48744 0.43578 7 Np quaterner positive nitrogen 87c 0.80779 0.28151 7 NX sp2 nitrogen 88c 0.98950 0.48874 7 NZ sp nitrogen 89c 1.15452 0.50429 7 NAr aromatic nitrogen 90c 0.77898 0.35097 7 NArp aromatic-positive nitrogen 91c 0.71947 0.31954 7 NCO amide nitrogen 92c 0.96737 0.43334 8 OH OH oxigen 93c 0.65563 0.29674 8 Op positive oxigen 94c 1.10775 0.41545 8 On oxidanione 95c 0.81599 0.35269 8 OAr aromatic oxigen 96c 0.82099 0.26674 8 Oox oxo-oxigen 97c 0.84823 0.33950 8 OS aetheric oxigen 98c 0.69661 0.25867 8 OAcH acidic OH oxigen 99c 0.77413 0.25308 16 SAc sulphate sulphur 100c 0.88299 0.41097 16 ST sp3 sulphur 101c 0.84700 0.26207 16 SX sp2 sulphur 102c 103c currently only common types implemented HCNOF Mg P S Cl Br 104c 105 data ((peparm(i,k),k=1,2),i=1,10) / 106 & 0.68903,0.39888, 107 & 0.71998,0.34395, 108 & 0.97104,0.46763, 109 & 0.95419,0.43623, 110 & 0.86408,0.40514, 111 & 0.67920,0.17446, 112 & 0.98387,0.13229, 113 & 0.78155,0.24617, 114 & 0.76957,0.17213, 115 & 0.77709,0.15474/ 116 117c from ESFF parameterisation 118c 119c 1,2 H,H+ 120c 3-5 Csp3, Csp2, Csp 121c 6-8 Nsp3, Nsp2, Nsp 122c 9-11 Osp3, Osp2, Osp 123c 12 F 124c 13-14 Psp3, Psp2 125c 15-16 Ssp3, Ssp2 126c 17 Cl 127c 18 Br 128c 19 I 129c 130c hybridisation specific currently not used, only first atom 131 132 data ((esffprm(i,k),k=1,2),i=1,19) / 133 & 0.233,0.440,0.252,0.440, 134 & 0.2904,0.382,0.316,0.384,0.366,0.389, 135 & 0.392,0.454,0.419,0.454,0.490,0.464, 136 & 0.480,0.518,0.527,0.519,0.625,0.529, 137 & 0.431,0.582, 138 & 0.306,0.309,0.323,0.311, 139 & 0.365,0.346,0.399,0.349, 140 & 0.328,0.374, 141 & 0.300,0.343, 142 & 0.358,0.304/ 143 144c eem NPA supported elements HCNOF 145 146 data (ipt(i,1),i=1,mxel) /1,0, 147 1 0,0,0,6,7,8,9,0, 148 2 0,0,0,0,0,0,0,0, 149 3 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 150 4 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 151 5 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 152 3 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 153 4 0,0,0,0,0/ 154 155c eem Mull supported elements HCNOF + est. P Cl 156 157 data (ipt(i,2),i=1,mxel) /1,0, 158 1 0,0,0,6,7,8,9,0, 159 2 0,0,0,0,15,0,17,0, 160 3 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 161 4 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 162 5 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 163 3 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 164 4 0,0,0,0,0/ 165 166c eem PESP supported elements HCNOF Mg P S Cl Br 167 168 data (ipt(i,3),i=1,mxel) /1,0, 169 1 0,0,0,2,3,4,5,0, 170 2 0,6,0,0,7,8,9,0, 171 3 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10,0, 172 4 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 173 5 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 174 3 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 175 4 0,0,0,0,0/ 176 177 178c eem ESFF supported elements HCNOF P S Cl Br I 179c plus hybridisation specific, currently not used 180 181 data (ipt(i,4),i=1,mxel) /1,0, 182 1 0,0,0,3,6,9,12,0, 183 2 0,0,0,0,13,15,17,0, 184 3 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,18,0, 185 4 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,19,0, 186 5 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 187 3 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 188 4 0,0,0,0,0/ 189 190 istat = 0 191 192 icnt = 0 193 ioverf = 0 194 195 do i=1,iatoms 196 if (iasel.gt.0.or.(iasel.lt.-3.and.iasel.eq.iresid(i)) ) 197 & then 198 199 iexcl = 0 200 do k=1,mexcl 201 if (ianz(i).eq.ianexc(k)) then 202 iexcl = 1 203 qat(i) = qexcl(k) 204 print*,' ' 205 print*,'Excluded metal center ',ianz(i) 206 print*,'Set formal charge: ',qat(i) 207 print*,' ' 208 endif 209 end do 210 211 if (ianz(i).eq.100) iexcl = 1 212 213 if (iexcl.eq.0) then 214 icnt = icnt + 1 215 if (icnt.lt.mxeat) then 216 ipntr(icnt) = i 217 do j=1,2 218 ipatm = ipt(ianz(i),iop) 219 if (ipatm.eq.0) then 220 call inferr( 221 & 'no parameters for element '//elemnt(ianz(i)),0) 222 ihasq = 0 223 istat = 1 224 return 225 endif 226 if (iop.eq.1) then 227c NPA 228 var(icnt,j) = param(ipatm,j) 229 elseif (iop.eq.2) then 230c Mull 231 var(icnt,j) = pparam(ipatm,j) 232 elseif (iop.eq.3) then 233c PESP 234 var(icnt,j) = peparm(ipatm,j) 235 elseif (iop.eq.4) then 236c ESFF 237 var(icnt,j) = esffprm(ipatm,j) 238 endif 239 end do 240c print*,i,' ',ianz(i),(var(icnt,j),j=1,2) 241 else 242 icnt = mxeat 243 ioverf = 1 244 endif 245 246 endif 247 endif 248 end do 249 250 if (ioverf.eq.1) then 251 call inferr( 252 & 'increase parameter mxeat in subroutine valdis',0) 253 endif 254 255 numat = icnt 256 257 return 258 end 259 260 subroutine eemcald(var,ipntr,numat,ianz,coo,qat) 261c 262 implicit double precision(a-h,o-z), integer (i-n) 263 parameter (mxeat=300) 264 common /charge/ dipo(3),ihasq,ihsdp,iqon,idipon 265 common /totchg/ itot 266 common /athlp/ iatoms, mxnat 267 dimension det(2) 268 dimension coo(3,*),qat(*),ianz(*) 269 dimension var(mxeat,2),ipntr(mxeat),work(mxeat+1,mxeat+1) 270 dimension x(mxeat+1,mxeat+1),y(mxeat+1), ipvt(mxeat+1),q(mxeat+1) 271 272 do i=1,(numat-1) 273 x(i,i) = 2*var(i,2) 274 x(i,numat+1) = -1 275 x(numat+1,i) = 1 276 277 do j=i+1,numat 278 x(i,j) = ( coo(1,ipntr(i)) - coo(1,ipntr(j)) )**2 + 279 & ( coo(2,ipntr(i)) - coo(2,ipntr(j)) )**2 + 280 & ( coo(3,ipntr(i)) - coo(3,ipntr(j)) )**2 281 x(i,j) = 1.0d0/dsqrt(x(i,j)) 282 x(j,i) = x(i,j) 283 end do 284 285 end do 286 x(numat,numat) = 2*var(numat,2) 287 x(numat,numat+1) = -1 288 x(numat+1,numat) = 1 289 x(numat+1,numat+1) = 0 290 291 y(numat+1) = dble(itot) 292 293 do i=1,numat 294 y(i) = -var(i,1) 295 end do 296 297 call dgefa(x,mxeat+1,numat+1,ipvt,info) 298 call dgedi(x,mxeat+1,(numat+1),ipvt,det,work,01) 299 300 call mtmul(x,y,q,numat) 301 302 write (*,'(a)') ' ' 303 write (*,'(a)') 'EEM charges' 304 write (*,'(a)') ' ' 305 qtot = 0.0d0 306 do i=1,(numat) 307 j = ipntr(i) 308 qat(j) = q(i) 309 qtot = qtot + q(i) 310 write (*,'(i5,1x,i3,1x,f10.5)') j,ianz(j),qat(j) 311 end do 312 write (*,'(a)') ' ' 313 write (*,'(a,f10.3)') 'Sum of EEM charges = ',qtot 314 ihasq = 1 315 316 return 317 end 318 319c 320 subroutine mtmul(a,y,b,numat) 321c 322 parameter (mxeat=300) 323c 324 integer i,j,numat 325 double precision a(mxeat+1,mxeat+1),y(mxeat+1),b(mxeat+1) 326c 327 do i=1,numat+1 328 329 b(i) = 0.0 330 331 do j=1,numat+1 332 b(i) = b(i) + a(i,j)*y(j) 333 end do 334 335 end do 336c 337 return 338 end 339c 340 subroutine dgedi(a,lda,n,ipvt,det,work,job) 341 implicit double precision(a-h,o-z) 342 dimension a(lda,*),det(2),work(*),ipvt(*) 343c 344c dgedi computes the determinant and inverse of a matrix 345c using the factors computed by dgeco or dgefa. 346c 347c on entry 348c 349c a double precision(lda, n) 350c the output from dgeco or dgefa. 351c 352c lda integer 353c the leading dimension of the array a . 354c 355c n integer 356c the order of the matrix a . 357c 358c ipvt integer(n) 359c the pivot vector from dgeco or dgefa. 360c 361c work double precision(n) 362c work vector. contents destroyed. 363c 364c job integer 365c = 11 both determinant and inverse. 366c = 01 inverse only. 367c = 10 determinant only. 368c 369c on return 370c 371c a inverse of original matrix if requested. 372c otherwise unchanged. 373c 374c det double precision(2) 375c determinant of original matrix if requested. 376c otherwise not referenced. 377c determinant = det(1) * 10.0**det(2) 378c with 1.0 .le. abs(det(1)) .lt. 10.0 379c or det(1) .eq. 0.0 . 380c 381c error condition 382c 383c a division by zero will occur if the input factor contains 384c a zero on the diagonal and the inverse is requested. 385c it will not occur if the subroutines are called correctly 386c and if dgeco has set rcond .gt. 0.0 or dgefa has set 387c info .eq. 0 . 388c 389c linpack. this version dated 08/14/78 . 390c cleve moler, university of new mexico, argonne national lab. 391c 392c subroutines and functions 393c 394c blas daxpy,dscal,dswap 395c fortran abs,mod 396c 397c compute determinant 398c 399 if (job/10 .eq. 0) go to 70 400 det(1) = 1.0d+00 401 det(2) = 0.0d+00 402 ten = 10.0d+00 403 do 50 i = 1, n 404 if (ipvt(i) .ne. i) det(1) = -det(1) 405 det(1) = a(i,i)*det(1) 406c ...exit 407 if (det(1) .eq. 0.0d+00) go to 60 408 10 if (abs(det(1)) .ge. 1.0d+00) go to 20 409 det(1) = ten*det(1) 410 det(2) = det(2) - 1.0d+00 411 go to 10 412 20 continue 413 30 if (abs(det(1)) .lt. ten) go to 40 414 det(1) = det(1)/ten 415 det(2) = det(2) + 1.0d+00 416 go to 30 417 40 continue 418 50 continue 419 60 continue 420 70 continue 421c 422c compute inverse(u) 423c 424 if (mod(job,10) .eq. 0) return 425 do k = 1, n 426 a(k,k) = 1.0d+00/a(k,k) 427 t = -a(k,k) 428 call dscal(k-1,t,a(1,k),1) 429 kp1 = k + 1 430 if (n .lt. kp1) go to 90 431 do j = kp1, n 432 t = a(k,j) 433 a(k,j) = 0.0d+00 434 call daxpy(k,t,a(1,k),1,a(1,j),1) 435 end do 436 90 continue 437 end do 438c 439c form inverse(u)*inverse(l) 440c 441 nm1 = n - 1 442 if (nm1 .lt. 1) return 443 do kb = 1, nm1 444 k = n - kb 445 kp1 = k + 1 446 do i = kp1, n 447 work(i) = a(i,k) 448 a(i,k) = 0.0d+00 449 end do 450 do j = kp1, n 451 t = work(j) 452 call daxpy(n,t,a(1,j),1,a(1,k),1) 453 end do 454 l = ipvt(k) 455 if (l .ne. k) call dswap(n,a(1,k),1,a(1,l),1) 456 end do 457 458 return 459 end 460 461 subroutine dgefa(a,lda,n,ipvt,info) 462 implicit double precision(a-h,o-z) 463 dimension a(lda,*),ipvt(*) 464c 465c dgefa factors a double precision matrix by gaussian elimination. 466c 467c dgefa is usually called by dgeco, but it can be called 468c directly with a saving in time if rcond is not needed. 469c (time for dgeco) = (1 + 9/n)*(time for dgefa) . 470c 471c on entry 472c 473c a double precision(lda, n) 474c the matrix to be factored. 475c 476c lda integer 477c the leading dimension of the array a . 478c 479c n integer 480c the order of the matrix a . 481c 482c on return 483c 484c a an upper triangular matrix and the multipliers 485c which were used to obtain it. 486c the factorization can be written a = l*u where 487c l is a product of permutation and unit lower 488c triangular matrices and u is upper triangular. 489c 490c ipvt integer(n) 491c an integer vector of pivot indices. 492c 493c info integer 494c = 0 normal value. 495c = k if u(k,k) .eq. 0.0 . this is not an error 496c condition for this subroutine, but it does 497c indicate that dgesl or dgedi will divide by zero 498c if called. use rcond in dgeco for a reliable 499c indication of singularity. 500c 501c linpack. this version dated 08/14/78 . 502c cleve moler, university of new mexico, argonne national lab. 503c 504c subroutines and functions 505c 506c blas daxpy,dscal,idamax 507c 508c gaussian elimination with partial pivoting 509c 510 info = 0 511 nm1 = n - 1 512 513 if (nm1.lt.1) goto 70 514 515 do k=1,nm1 516 kp1 = k + 1 517 518c find l = pivot index 519 520 l = idamax(n-k+1,a(k,k),1) + k - 1 521 ipvt(k) = l 522 523c zero pivot implies this column already triangularized 524 525 if (a(l,k).eq.0.0d0) then 526 info = k 527 else 528 529c interchange if necessary 530 531 if (l.ne.k) then 532 t = a(l,k) 533 a(l,k) = a(k,k) 534 a(k,k) = t 535 endif 536 537c compute multipliers 538 539 t = -1.0d+00/a(k,k) 540 call dscal(n-k,t,a(k+1,k),1) 541 542c row elimination with column indexing 543 544 do j = kp1, n 545 546 t = a(l,j) 547 548 if (l.ne.k) then 549 a(l,j) = a(k,j) 550 a(k,j) = t 551 endif 552 553 call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 554 end do 555 556 endif 557 558 end do 559 560 70 continue 561 562 ipvt(n) = n 563 if (a(n,n) .eq. 0.0d+00) info = n 564 565 return 566 end 567 568 double precision function dasum(n,dx,incx) 569 570c takes the sum of the absolute values. 571 572 double precision dx(1),dtemp 573 integer i,incx,m,mp1,n,nincx 574 575 dasum = 0.0d+00 576 dtemp = 0.0d+00 577 578 if (n.le.0) return 579 580 if (incx.ne.1) then 581 582 nincx = n*incx 583 584 do i=1,nincx,incx 585 dtemp = dtemp + abs(dx(i)) 586 end do 587 588 dasum = dtemp 589 590 else 591 592 m = mod(n,6) 593 594 if (m.ne.0) then 595 596 do i = 1,m 597 dtemp = dtemp + abs(dx(i)) 598 end do 599 600 if (n.lt.6) then 601 dasum = dtemp 602 return 603 endif 604 605 endif 606 607 mp1 = m + 1 608 609 do i = mp1,n,6 610 dtemp = dtemp + abs(dx(i)) + abs(dx(i + 1)) + abs(dx(i + 2)) 611 & + abs(dx(i + 3)) + abs(dx(i + 4)) + abs(dx(i + 5)) 612 end do 613 614 dasum = dtemp 615 616 endif 617 618 return 619 end 620 621 subroutine daxpy(n,da,dx,incx,dy,incy) 622 implicit double precision(a-h,o-z) 623 dimension dx(1),dy(1) 624c 625c constant times a vector plus a vector. 626c dy(i) = dy(i) + da * dx(i) 627c uses unrolled loops for increments equal to one. 628c jack dongarra, linpack, 3/11/78. 629c 630 if (n.le.0) return 631 if (da .eq. 0.0d0) return 632 633 if (incx.eq.1.and.incy.eq.1) then 634 m = mod(n,4) 635 if (m.ne.0) then 636 637 do i = 1,m 638 dy(i) = dy(i) + da*dx(i) 639 end do 640 641 if (n.lt.4) return 642 endif 643 644 mp1 = m + 1 645 646 do i = mp1,n,4 647 dy(i) = dy(i) + da*dx(i) 648 dy(i + 1) = dy(i + 1) + da*dx(i + 1) 649 dy(i + 2) = dy(i + 2) + da*dx(i + 2) 650 dy(i + 3) = dy(i + 3) + da*dx(i + 3) 651 end do 652 653 else 654 655c code for unequal increments or equal increments 656c not equal to 1 657 658 ix = 1 659 iy = 1 660 661 if (incx.lt.0) ix = (-n+1)*incx + 1 662 if (incy.lt.0) iy = (-n+1)*incy + 1 663 664 do i = 1,n 665 dy(iy) = dy(iy) + da*dx(ix) 666 ix = ix + incx 667 iy = iy + incy 668 end do 669 670 endif 671 672 return 673 end 674 675 subroutine dcopy(n,dx,incx,dy,incy) 676 implicit double precision(a-h,o-z) 677 dimension dx(*),dy(*) 678c 679c copies a vector. 680c dy(i) <== dx(i) 681c uses unrolled loops for increments equal to one. 682c jack dongarra, linpack, 3/11/78. 683c 684 if(n.le.0)return 685 if(incx.eq.1.and.incy.eq.1)go to 20 686c 687c code for unequal increments or equal increments 688c not equal to 1 689c 690 ix = 1 691 iy = 1 692 if(incx.lt.0)ix = (-n+1)*incx + 1 693 if(incy.lt.0)iy = (-n+1)*incy + 1 694 do 10 i = 1,n 695 dy(iy) = dx(ix) 696 ix = ix + incx 697 iy = iy + incy 698 10 continue 699 return 700c 701c code for both increments equal to 1 702c 703c 704c clean-up loop 705c 706 20 m = mod(n,7) 707 if( m .eq. 0 ) go to 40 708 do 30 i = 1,m 709 dy(i) = dx(i) 710 30 continue 711 if( n .lt. 7 ) return 712 40 mp1 = m + 1 713 do 50 i = mp1,n,7 714 dy(i) = dx(i) 715 dy(i + 1) = dx(i + 1) 716 dy(i + 2) = dx(i + 2) 717 dy(i + 3) = dx(i + 3) 718 dy(i + 4) = dx(i + 4) 719 dy(i + 5) = dx(i + 5) 720 dy(i + 6) = dx(i + 6) 721 50 continue 722 return 723 end 724 725 double precision function ddot(n,dx,incx,dy,incy) 726 implicit double precision(a-h,o-z) 727 dimension dx(1),dy(1) 728c 729c forms the dot product of two vectors. 730c dot = dx(i) * dy(i) 731c uses unrolled loops for increments equal to one. 732c jack dongarra, linpack, 3/11/78. 733c 734 ddot = 0.0d+00 735 dtemp = 0.0d+00 736 if(n.le.0)return 737 if(incx.eq.1.and.incy.eq.1)go to 20 738c 739c code for unequal increments or equal increments 740c not equal to 1 741c 742 ix = 1 743 iy = 1 744 if(incx.lt.0)ix = (-n+1)*incx + 1 745 if(incy.lt.0)iy = (-n+1)*incy + 1 746 do 10 i = 1,n 747 dtemp = dtemp + dx(ix)*dy(iy) 748 ix = ix + incx 749 iy = iy + incy 750 10 continue 751 ddot = dtemp 752 return 753c 754c code for both increments equal to 1 755c 756c 757c clean-up loop 758c 759 20 m = mod(n,5) 760 if( m .eq. 0 ) go to 40 761 do 30 i = 1,m 762 dtemp = dtemp + dx(i)*dy(i) 763 30 continue 764 if( n .lt. 5 ) go to 60 765 40 mp1 = m + 1 766 do 50 i = mp1,n,5 767 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + 768 * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 769 50 continue 770 60 ddot = dtemp 771 return 772 end 773 774 double precision function dnrm2 ( n, dx, incx) 775 integer next 776 double precision dx(1), cutlo, cuthi, hitest, sum, xmax,zero,one 777 data zero, one /0.0d+00, 1.0d+00/ 778c 779c euclidean norm of the n-vector stored in dx() with storage 780c increment incx . 781c if n .le. 0 return with result = 0. 782c if n .ge. 1 then incx must be .ge. 1 783c 784c c.l.lawson, 1978 jan 08 785c 786c four phase method using two built-in constants that are 787c hopefully applicable to all machines. 788c cutlo = maximum of sqrt(u/eps) over all known machines. 789c cuthi = minimum of sqrt(v) over all known machines. 790c where 791c eps = smallest no. such that eps + 1. .gt. 1. 792c u = smallest positive no. (underflow limit) 793c v = largest no. (overflow limit) 794c 795c brief outline of algorithm.. 796c 797c phase 1 scans zero components. 798c move to phase 2 when a component is nonzero and .le. cutlo 799c move to phase 3 when a component is .gt. cutlo 800c move to phase 4 when a component is .ge. cuthi/m 801c where m = n for x() real and m = 2*n for complex. 802c 803c values for cutlo and cuthi.. 804c from the environmental parameters listed in the imsl converter 805c document the limiting values are as follows.. 806c cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are 807c univac and dec at 2**(-103) 808c thus cutlo = 2**(-51) = 4.44089e-16 809c cuthi, s.p. v = 2**127 for univac, honeywell, and dec. 810c thus cuthi = 2**(63.5) = 1.30438e19 811c cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. 812c thus cutlo = 2**(-33.5) = 8.23181d-11 813c cuthi, d.p. same as s.p. cuthi = 1.30438d+19 814c data cutlo, cuthi / 8.232d-11, 1.304d+19 / 815c data cutlo, cuthi / 4.441e-16, 1.304e19 / 816 data cutlo, cuthi / 8.232d-11, 1.304d+19 / 817c 818 j = 0 819 820 if (n.le.0) then 821 dnrm2 = zero 822 return 823 endif 824c 825 assign 30 to next 826 sum = zero 827 nn = n * incx 828c begin main loop 829 i = 1 830 20 go to next,(30, 50, 70, 110) 831 30 if( abs(dx(i)) .gt. cutlo) go to 85 832 assign 50 to next 833 xmax = zero 834c 835c phase 1. sum is zero 836c 837 50 if( dx(i) .eq. zero) go to 200 838 if( abs(dx(i)) .gt. cutlo) go to 85 839c 840c prepare for phase 2. 841 assign 70 to next 842 go to 105 843c 844c prepare for phase 4. 845c 846 100 i = j 847 assign 110 to next 848 sum = (sum / dx(i)) / dx(i) 849 105 xmax = abs(dx(i)) 850 go to 115 851c 852c phase 2. sum is small. 853c scale to avoid destructive underflow. 854c 855 70 if( abs(dx(i)) .gt. cutlo ) go to 75 856c 857c common code for phases 2 and 4. 858c in phase 4 sum is large. scale to avoid overflow. 859c 860 110 if (abs(dx(i)).le.xmax) go to 115 861 sum = one + sum * (xmax / dx(i))**2 862 xmax = abs(dx(i)) 863 go to 200 864c 865 115 sum = sum + (dx(i)/xmax)**2 866 go to 200 867c 868c 869c prepare for phase 3. 870c 871 75 sum = (sum * xmax) * xmax 872c 873c 874c for real or d.p. set hitest = cuthi/n 875c for complex set hitest = cuthi/(2*n) 876c 877 85 hitest = cuthi/n 878c 879c phase 3. sum is mid-range. no scaling. 880c 881 do 95 j =i,nn,incx 882 if(abs(dx(j)) .ge. hitest) go to 100 883 95 sum = sum + dx(j)**2 884 dnrm2 = sqrt( sum ) 885 return 886c 887 200 continue 888 i = i + incx 889 if ( i .le. nn ) go to 20 890c 891c end of main loop. 892c 893c compute square root and adjust for scaling. 894c 895 dnrm2 = xmax * sqrt(sum) 896 300 continue 897 return 898 end 899 900 subroutine drot(n,dx,incx,dy,incy,c,s) 901 implicit double precision(a-h,o-z) 902 dimension dx(1),dy(1) 903c 904c applies a plane rotation. 905c dx(i) = c*dx(i) + s*dy(i) 906c dy(i) = -s*dx(i) + c*dy(i) 907c jack dongarra, linpack, 3/11/78. 908c 909 if(n.le.0)return 910 if(incx.eq.1.and.incy.eq.1)go to 20 911c 912c code for unequal increments or equal increments not equal 913c to 1 914c 915 ix = 1 916 iy = 1 917 if(incx.lt.0)ix = (-n+1)*incx + 1 918 if(incy.lt.0)iy = (-n+1)*incy + 1 919 do i = 1,n 920 dtemp = c*dx(ix) + s*dy(iy) 921 dy(iy) = c*dy(iy) - s*dx(ix) 922 dx(ix) = dtemp 923 ix = ix + incx 924 iy = iy + incy 925 end do 926 927 return 928c 929c code for both increments equal to 1 930c 931 20 do i = 1,n 932 dtemp = c*dx(i) + s*dy(i) 933 dy(i) = c*dy(i) - s*dx(i) 934 dx(i) = dtemp 935 end do 936 937 return 938 end 939 940 subroutine drotg(da,db,c,s) 941c 942c construct givens plane rotation. 943c jack dongarra, linpack, 3/11/78. 944c 945 double precision da,db,c,s,roe,scale,r,z 946 double precision zero, one 947c 948 parameter (zero=0.0d+00, one=1.0d+00) 949c 950c----------------------------------------------------------------------- 951c 952c 953 roe = db 954 if( abs(da) .gt. abs(db) ) roe = da 955 scale = abs(da) + abs(db) 956 if ( scale .eq. zero ) then 957 c = one 958 s = zero 959 r = zero 960 else 961 r = scale*sqrt((da/scale)**2 + (db/scale)**2) 962 r = sign(one,roe)*r 963 c = da/r 964 s = db/r 965 endif 966 z = one 967 if( abs(da) .gt. abs(db) ) z = s 968 if( abs(db) .ge. abs(da) .and. c .ne. zero ) z = one/c 969 da = r 970 db = z 971 972 return 973 end 974 975 subroutine dscal(n,da,dx,incx) 976 implicit double precision(a-h,o-z) 977 dimension dx(1) 978c 979c scales a vector by a constant. 980c dx(i) = da * dx(i) 981c uses unrolled loops for increment equal to one. 982c jack dongarra, linpack, 3/11/78. 983c 984 if(n.le.0)return 985 if(incx.eq.1)go to 20 986c 987c code for increment not equal to 1 988c 989 nincx = n*incx 990 991 do i = 1,nincx,incx 992 dx(i) = da*dx(i) 993 end do 994 995 return 996c 997c code for increment equal to 1 998c 999c 1000c clean-up loop 1001c 1002 20 m = mod(n,5) 1003 if ( m .eq. 0 ) go to 40 1004 do i = 1,m 1005 dx(i) = da*dx(i) 1006 end do 1007 if( n .lt. 5 ) return 1008 1009 40 mp1 = m + 1 1010 do i = mp1,n,5 1011 dx(i) = da*dx(i) 1012 dx(i + 1) = da*dx(i + 1) 1013 dx(i + 2) = da*dx(i + 2) 1014 dx(i + 3) = da*dx(i + 3) 1015 dx(i + 4) = da*dx(i + 4) 1016 end do 1017 1018 return 1019 end 1020 1021 subroutine dswap(n,dx,incx,dy,incy) 1022 implicit double precision(a-h,o-z) 1023 dimension dx(1),dy(1) 1024c 1025c interchanges two vectors. 1026c dx(i) <==> dy(i) 1027c uses unrolled loops for increments equal one. 1028c jack dongarra, linpack, 3/11/78. 1029c 1030 if(n.le.0)return 1031 if(incx.eq.1.and.incy.eq.1)go to 20 1032c 1033c code for unequal increments or equal increments not equal 1034c to 1 1035c 1036 ix = 1 1037 iy = 1 1038 if(incx.lt.0)ix = (-n+1)*incx + 1 1039 if(incy.lt.0)iy = (-n+1)*incy + 1 1040 do i = 1,n 1041 dtemp = dx(ix) 1042 dx(ix) = dy(iy) 1043 dy(iy) = dtemp 1044 ix = ix + incx 1045 iy = iy + incy 1046 end do 1047 1048 return 1049c 1050c code for both increments equal to 1 1051c 1052c 1053c clean-up loop 1054c 1055 20 m = mod(n,3) 1056 if( m .eq. 0 ) go to 40 1057 do 30 i = 1,m 1058 dtemp = dx(i) 1059 dx(i) = dy(i) 1060 dy(i) = dtemp 1061 30 continue 1062 if( n .lt. 3 ) return 1063 40 mp1 = m + 1 1064 do 50 i = mp1,n,3 1065 dtemp = dx(i) 1066 dx(i) = dy(i) 1067 dy(i) = dtemp 1068 dtemp = dx(i + 1) 1069 dx(i + 1) = dy(i + 1) 1070 dy(i + 1) = dtemp 1071 dtemp = dx(i + 2) 1072 dx(i + 2) = dy(i + 2) 1073 dy(i + 2) = dtemp 1074 50 continue 1075 return 1076 end 1077 1078 integer function idamax(n,dx,incx) 1079 implicit double precision(a-h,o-z) 1080 dimension dx(1) 1081c 1082c finds the index of element having max. absolute value. 1083c jack dongarra, linpack, 3/11/78. 1084c 1085 idamax = 0 1086 if( n .lt. 1 ) return 1087 idamax = 1 1088 if(n.eq.1)return 1089 if(incx.eq.1)go to 20 1090c 1091c code for increment not equal to 1 1092c 1093 ix = 1 1094 rmax = abs(dx(1)) 1095 ix = ix + incx 1096 1097 do i = 2,n 1098 if (abs(dx(ix)).gt.rmax) then 1099 idamax = i 1100 rmax = abs(dx(ix)) 1101 endif 1102 ix = ix + incx 1103 end do 1104 1105 return 1106c 1107c code for increment equal to 1 1108c 1109 20 rmax = abs(dx(1)) 1110 do i = 2,n 1111 if(abs(dx(i)).gt.rmax) then 1112 idamax = i 1113 rmax = abs(dx(i)) 1114 endif 1115 end do 1116 1117 return 1118 end 1119 1120 subroutine dgemv(forma,m,n,alpha,a,lda,x,incx,beta,y,incy) 1121 implicit double precision(a-h,o-z) 1122 character*1 forma 1123 dimension a(lda,*),x(*),y(*) 1124 parameter (zero=0.0d+00, one=1.0d+00) 1125c 1126c clone of -dgemv- written by mike schmidt 1127c 1128 locy = 1 1129 if(forma.eq.'t') go to 200 1130c 1131c y = alpha * a * x + beta * y 1132c 1133 if (alpha.eq.one .and. beta.eq.zero) then 1134 do i=1,m 1135 y(locy) = ddot(n,a(i,1),lda,x,incx) 1136 locy = locy+incy 1137 end do 1138 else 1139 do i=1,m 1140 y(locy) = alpha*ddot(n,a(i,1),lda,x,incx) + beta*y(locy) 1141 locy = locy+incy 1142 end do 1143 end if 1144 return 1145c 1146c y = alpha * a-transpose * x + beta * y 1147c 1148 200 continue 1149 1150 if(alpha.eq.one .and. beta.eq.zero) then 1151 do i=1,n 1152 y(locy) = ddot(m,a(1,i),1,x,incx) 1153 locy = locy+incy 1154 end do 1155 else 1156 do i=1,n 1157 y(locy) = alpha*ddot(m,a(1,i),1,x,incx) + beta*y(locy) 1158 locy = locy+incy 1159 end do 1160 end if 1161 1162 return 1163 end 1164 1165 subroutine filnam(filn,nnn) 1166 character name*70,filn*(*) 1167 1168 i1 = 0 1169 do i = 1,nnn 1170 if (filn(i:i).ne.' ') then 1171 i1 = i1 + 1 1172 name(i1:i1) = filn(i:i) 1173 endif 1174 end do 1175 1176 do i = 1,nnn 1177 filn(i:i) = ' ' 1178 end do 1179 1180 do i = 1,i1 1181 filn(i:i) = name(i:i) 1182 end do 1183 1184 return 1185 end 1186 1187 subroutine distchd(coo) 1188c 1189c 1190 implicit double precision(a-h,o-z) 1191 common /athlp/ iatoms, mxnat 1192 dimension coo(3,*) 1193 1194 tol = 0.1d0 1195 do i=1,iatoms 1196 do j=i+1,iatoms 1197 dd = 0.0d0 1198 do k=1,3 1199 d1 = coo(k,i) - coo(k,j) 1200 d2 = d1*d1 1201 dd = dd + d2 1202 end do 1203 if (dd.lt.tol) print*,'close ',i,' ',j 1204 end do 1205 end do 1206 1207 return 1208 end 1209 1210 subroutine sigini 1211 implicit double precision(a-h,o-z) 1212 parameter (mxsigm=16) 1213 parameter (mxmol2=41) 1214 common /sigma/ siga(mxsigm),sigb(mxsigm),sigc(mxsigm), 1215 & sigd(mxsigm),impmol2(mxmol2) 1216c data sigs /"H ","C3","C2","C1","N3","N2","N1","O3","O2","F ", 1217c & "Cl","Br","I ","S3","P ","DU"/ 1218 data siga /7.17,7.98,8.79,10.39,11.54,12.87,15.68,14.18, 1219 & 17.07,14.66, 11.00, 10.08,9.90,10.14,8.9,0.0/ 1220 data sigb /6.24,9.18,9.32,9.45,10.82,11.15,11.70,12.92, 1221 & 13.79,13.85,9.69,8.47,7.96,9.13,8.24,0.0/ 1222 data sigc /-0.56,1.88,1.51,0.73,1.36,0.85,-0.27,1.39, 1223 & 0.47,2.31,1.35,1.16,0.96,1.38,0.96,0.0/ 1224 data impmol2 /16,16,16,16,2,3,4,3,3,5,6,7,6,6,6,5, 1225 & 8,9,9,8,8,14,14,14,14,15,1,1,1,10,11,12,13, 1226 & 16,16,16,16,16,16,16,16/ 1227 1228c sybyl atom types to gasteiger types 1229c1 any 1230c2 hal 1231c3 het 1232c4 hev 1233c5 C.3 2 1234c6 C.2 3 1235c7 C.1 4 1236c8 C.ar 3 1237c9 C.cat ? 1238c10 N.3 5 1239c11 N.2 6 1240c12 N.1 7 1241c13 N.ar 6 1242c14 N.am 6 1243c15 N.pl3 6 1244c16 N.4 5 1245c17 O.3 8 1246c18 O.2 9 1247c19 O.co2 9 1248c20 O.spc 8 1249c21 O.t3p 8 1250c22 S.3 14 1251c23 S.2 14 1252c24 S.O 14 1253c25 S.O2 14 1254c26 P.3 15 1255c27 H 1 1256c28 H.spc 1 1257c29 H.t3p 1 1258c30 F 10 1259c31 Cl 11 1260c32 Br 12 1261c33 I 13 1262c34 Si 1263c35 Lp 1264c36 Du 1265c37 Na 1266c38 K 1267c39 Ca 1268c40 Li 1269c41 Al 1270 1271c unknown is DU 16 1272 1273 do i=1,mxsigm 1274 sigd(i) = siga(i) + sigb(i) + sigc(i) 1275 end do 1276 sigd(15) = 1.0d0 1277 1278 return 1279 end 1280 1281 subroutine inigad(qtot,qat,ianz,iconn,ityp) 1282 implicit double precision(a-h,o-z) 1283 parameter (numat1=20000) 1284 parameter (mxcon=10) 1285 parameter (mxsigm=16) 1286 parameter (mxmol2=41) 1287 common /sigma/ siga(mxsigm),sigb(mxsigm),sigc(mxsigm), 1288 & sigd(mxsigm),impmol2(mxmol2) 1289 common /athlp/ iatoms, mxnat 1290 common /charge/ dipo(3),ihasq,ihsdp,iqon,idipon 1291 integer*2 ityp 1292 common /types/ iff 1293 dimension qat(*),ianz(*),iconn(mxcon+1,*),ityp(*) 1294 1295 ihasq = 0 1296 iff = 5 1297 call dotyp(0) 1298 qtot = 0.0d0 1299 1300 do i=1,iatoms 1301 qat(i) = 0.0d0 1302 ia = ianz(i) 1303 1304 if (ia.eq.6) then 1305 if (ityp(i).eq.9) qat(i) = 1.0d0 1306 elseif (ia.eq.7) then 1307 ibnd = 0 1308 do j=1,iconn(1,i) 1309 if (iconn(1+j,i).gt.0) ibnd = ibnd + 1 1310 end do 1311 if (ibnd.eq.4) qat(i) = 1.0d0 1312 elseif (ia.eq.8) then 1313c carboxyl 1314 if (ityp(i).eq.19) qat(i) = -0.5d0 1315 endif 1316 qtot = qtot + qat(i) 1317 end do 1318 1319 return 1320 end 1321 1322 subroutine clqgas(ishoh) 1323 implicit double precision(a-h,o-z) 1324 logical dozme 1325 common /getpnt/irtype,ipdbon,ipdbgro,ifav,ioxyz, 1326 & iconv,ircus,dozme 1327 parameter (mxheta=150) 1328 character*3 hetz 1329 common /clfstr/ ihashz,ihetq(mxheta),ihqset(mxheta),ihhadd(mxheta) 1330 & ,labhet(mxheta),ilcset,ligcat(mxheta),hetz(mxheta) 1331 istat = 0 1332 1333 if (ipdbon.eq.1) then 1334 1335 call numhet(nhmol) 1336 do i=4,nhmol 1337 if (i.ne.ishoh) then 1338 print*,' ' 1339 if (ihashz.eq.1) then 1340 print*,'HETATM residue ',hetz(i+1) 1341 else 1342 print*,'HETATM residue ',(i+1-4) 1343 endif 1344 print*,' ' 1345 call calgas(-i,istat) 1346 endif 1347 end do 1348 1349 print*,' ' 1350 print*,'WARNING: Total charges of the different HETATM residues' 1351 print*,' assumed zero. If this is NOT correct, assign ' 1352 print*,' HETATM charges individually, by clicking with' 1353 print*,' 2nd mouse button on the residue and select ' 1354 print*,' Calculate Charges' 1355 print*,' ' 1356 1357 else 1358 call calgas(1,istat) 1359 endif 1360 1361 if (istat.eq.1) call messg(14) 1362 1363 return 1364 end 1365 1366 subroutine calgad(iasel,istat,qat,ianz,iconn,iresid,ityp) 1367 implicit double precision(a-h,o-z) 1368 parameter (numat1=20000) 1369 parameter (mxcon=10) 1370 parameter (mxsigm=16) 1371 parameter (mxmol2=41) 1372 common /sigma/ siga(mxsigm),sigb(mxsigm),sigc(mxsigm), 1373 & sigd(mxsigm),impmol2(mxmol2) 1374 common /athlp/ iatoms, mxnat 1375 common /charge/ dipo(3),ihasq,ihsdp,iqon,idipon 1376 integer*2 ityp 1377 common /types/ iff 1378 dimension zz(numat1) 1379 dimension qat(*),ianz(*),iconn(mxcon+1,*),ityp(*),iresid(*) 1380 1381 call sigini 1382 1383 ihasq = 0 1384 istat = 0 1385 iff = 5 1386 call dotyp(0) 1387 1388 do i=1,iatoms 1389 if (iasel.gt.0.or.(iasel.lt.-3.and.iasel.eq.iresid(i)) ) 1390 & then 1391 1392 it = ityp(i) 1393 1394 if (it.gt.0.and.it.le.mxmol2) then 1395 l = impmol2(it) 1396 else 1397 print*,'Gasteiger: Untyped mol2 atom type.' 1398 istat = 1 1399 return 1400 endif 1401 1402 if (i.le.numat1) then 1403 zz(i) = siga(l) 1404 else 1405 print*,'array for gasteiger calculation to small' 1406 istat = 1 1407 return 1408 endif 1409 1410 qat(i) = 0.0d0 1411 ia = ianz(i) 1412 1413 if (ia.eq.6) then 1414 1415 if (ityp(i).eq.9) qat(i) = 1.0d0 1416 1417 elseif (ia.eq.7) then 1418 1419 ibnd = 0 1420 do j=1,iconn(1,i) 1421 if (iconn(1+j,i).gt.0) ibnd = ibnd + 1 1422 end do 1423 if (ibnd.eq.4) qat(i) = 1.0d0 1424 1425c if (N3+) qat(i) = 1.0d0 1426 1427 elseif (ia.eq.8) then 1428c carboxyl 1429 if (ityp(i).eq.19) qat(i) = -0.5d0 1430 1431 endif 1432 endif 1433 end do 1434 1435 fac = 1.0d0 1436 icnt = 0 1437 1438 do while (.true.) 1439 1440 fac = fac*0.5d0 1441 sd1 = 0.0d0 1442 1443 do i=1,iatoms 1444 1445 if (iasel.gt.0.or.(iasel.lt.-3.and.iasel.eq.iresid(i)) ) 1446 & then 1447 1448 l = impmol2(ityp(i)) 1449 1450 if (sigd(l).ne.1.0d0) then 1451 1452 qt = qat(i) 1453 do j=1,iconn(1,i) 1454 if (iconn(1+j,i).gt.0) then 1455 jj = iabs(iconn(1+j,i)) 1456 ll = impmol2(ityp(jj)) 1457 if (sigd(ll).ne.1.0d0) then 1458 sd2 = sigd(ll) 1459 if (i.le.numat1.and.jj.le.numat1) then 1460 if (zz(jj).gt.zz(i)) sd2 = sigd(l) 1461 endif 1462 if (ianz(i).eq.1.or.ianz(jj).eq.1) 1463 & sd2 = 20.02d0 1464 if (i.le.numat1.and.jj.le.numat1) then 1465 qat(i) = qat(i) + (zz(jj) - zz(i))*fac/sd2 1466 endif 1467 endif 1468 endif 1469 end do 1470 1471 qt = dabs(qat(i) - qt) 1472 if (qt.gt.sd1) sd1 = qt 1473 1474 endif 1475 1476 endif 1477 end do 1478 1479 if (sd1.ge.1.0d-3) then 1480 1481 do i=1,iatoms 1482 1483 if (iasel.gt.0.or.(iasel.lt.-3.and.iasel.eq.iresid(i)) ) 1484 & then 1485 l = impmol2(ityp(i)) 1486 zz(i) = siga(l) + sigb(l)*qat(i) + sigc(l)*qat(i)*qat(i) 1487 endif 1488 1489 end do 1490 1491 endif 1492 1493 icnt = icnt + 1 1494 1495 if (.not.(sd1.gt.1.0d-3.and.icnt.le.5)) goto 100 1496 end do 1497 1498100 continue 1499 1500 write (*,'(a)') ' ' 1501 write (*,'(a)') 'Gasteiger charges' 1502 write (*,'(a)') ' ' 1503 qtot = 0.0d0 1504 do i=1,iatoms 1505 if (iasel.gt.0.or.(iasel.lt.-3.and.iasel.eq.iresid(i)) ) 1506 & then 1507 qtot = qtot + qat(i) 1508 write (*,'(i5,1x,i3,1x,f10.5)') i,ianz(i),qat(i) 1509 endif 1510 end do 1511 write (*,'(a)') ' ' 1512 write (*,'(a,f10.3)') 'Sum of Gasteiger charges = ',qtot 1513 1514 ihasq = 1 1515 return 1516 end 1517 1518 subroutine clqeem(iop,ishoh) 1519 implicit double precision(a-h,o-z) 1520 logical dozme 1521 common /getpnt/irtype,ipdbon,ipdbgro,ifav,ioxyz, 1522 & iconv,ircus,dozme 1523 parameter (mxheta=150) 1524 character*3 hetz 1525 common /clfstr/ ihashz,ihetq(mxheta),ihqset(mxheta),ihhadd(mxheta) 1526 & ,labhet(mxheta),ilcset,ligcat(mxheta),hetz(mxheta) 1527 1528 istat = 0 1529 1530 if (ipdbon.eq.1) then 1531 1532 call numhet(nhmol) 1533 do i=4,nhmol 1534 if (i.ne.ishoh) then 1535 print*,' ' 1536 if (ihashz.eq.1) then 1537 print*,'HETATM residue ',hetz(i+1) 1538 else 1539 print*,'HETATM residue ',(i+1-4) 1540 endif 1541 print*,' ' 1542 call eem(iop,-i,istat) 1543 endif 1544 end do 1545 1546 print*,' ' 1547 print*,'WARNING: Total charges of the different HETATM residues' 1548 print*,' assumed zero. If this is NOT correct, assign ' 1549 print*,' HETATM charges individually, by clicking with' 1550 print*,' 2nd mouse button on the residue and select ' 1551 print*,' Calculate Charges' 1552 print*,' ' 1553 1554 else 1555 call eem(iop,1,istat) 1556 endif 1557 1558 if (istat.eq.1) call messg(14) 1559 1560 return 1561 end 1562 1563 subroutine clceem(ishoh) 1564 implicit double precision(a-h,o-z) 1565 logical dozme 1566 common /getpnt/irtype,ipdbon,ipdbgro,ifav,ioxyz, 1567 & iconv,ircus,dozme 1568 common /types/ iff 1569 istat = 0 1570 1571 if (iff.eq.7) then 1572 if (ipdbon.eq.1) then 1573 call numhet(nhmol) 1574 do i=4,nhmol 1575 if (i.ne.ishoh) then 1576 call eem(3,-i,istat) 1577 endif 1578 end do 1579 else 1580 call eem(3,1,istat) 1581 endif 1582 endif 1583 1584 print*,' ' 1585 print*,'EEM: Electronegativity Equalization Method:' 1586 print*,' ' 1587 print*,' Patrick Bultinck et al' 1588 print*,' J. Phys. Chem. A (2002)' 1589 print*,' ' 1590 1591 if (istat.eq.1) call messg(14) 1592 1593 return 1594 end 1595 1596