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