1      subroutine espod(x,y,z,epot,idebug,
2     &                 p)
3c THIS IS REALLY espot
4      implicit double precision(a-h,o-z)
5      parameter (numatm=2000)
6      parameter (numprm=1600)
7      parameter (numcex=numprm*3)
8      parameter (numtmp=4000)
9      logical oeerst
10      integer shella,shelln,shellt,shellc,shladf,aos
11      common/b/exx(numcex),
12     &         c1(numcex),c2(numcex),c3(numprm),c4(numprm),c5(numprm),
13     &         shladf(numprm),gx(numprm),gy(numprm),gz(numprm),
14     &         jan(numprm),shella(numprm),shelln(numprm),shellt(numprm),
15     &         shellc(numprm),aos(numprm),nshell,maxtyp
16      common /bounds/ ilow(5,5),iupp(5)
17      common /indxe/  jx(35),jy(35),jz(35),ix(20),iy(20),iz(20)
18      common /strt/   oeerst
19      common /moldat/ natoms, norbs, nelecs,nat(numatm)
20      common /pseudo / ipseud,ivale(numatm)
21      common /coord / xyz(3,numatm)
22      common /orbhlp/ mxorb,iuhf,ispd
23      common /extchg/ exchg(3,3),iexchg(3),nexchg
24      dimension tt(4),zz(4),pre(45),ca(35),cb(35),
25     &          coeffx(192),coeffy(192),coeffz(192)
26      dimension h(numtmp),f(numtmp),h2x(9),h2y(9),h2z(9),
27     &          h3x(16),h3y(16),h3z(16)
28      dimension p(*)
29      data jx/1,2,1,1,3,1,1,2,2,1,4,1,1,2,3,3,2,1,1,2,
30     &           5,1,1,4,4,2,1,2,1,3,3,1,3,2,2/
31      data jy/1,1,2,1,1,3,1,2,1,2,1,4,1,3,2,1,1,2,3,2,
32     &           1,5,1,2,1,4,4,1,2,3,1,3,2,3,2/
33      data jz/1,1,1,2,1,1,3,1,2,2,1,1,4,1,1,2,3,3,2,2,
34     &           1,1,5,1,2,1,2,4,4,1,3,3,2,2,3/
35      data ix/0,4,0,0,8,0,0,4,4,0,12,0,0,4,8,8,4,0,0,4/
36      data iy/0,0,4,0,0,8,0,4,0,4,0,12,0,8,4,0,0,4,8,4/
37      data iz/0,0,0,4,0,0,8,0,4,4,0,0,12,0,0,4,8,8,4,4/
38c iupp 1  4  10 20 35
39c      s  3p 6d 10f 15g
40c
41c ilow 0    1       2        3      4    <- shellt
42c---------------------------------------
43c 0    1(s) 1 (sp)  1 (spd)  1      1
44c 1    1    2 (p)   5       11     21
45c 2    1    1       5 (d)   11     21
46c 3    1    1       5       11 (f) 21
47c 4    1    1       5       11     21 (g)
48c---------------------------------------
49c ^
50c | shellc
51c
52      data iupp/1,4,10,20,35/
53      data ilow/5*1,1,2,5,11,21,1,1,5,11,21,1,1,5,11,21,
54     &          1,1,5,11,21/
55
56      if (oeerst) then
57         call epint
58         call setcon
59         oeerst = .false.
60      endif
61
62      epot = 0.0d0
63      pi = 4.d0*datan(1.d0)
64c
65c     loop over shell pairs.
66c
67      do ishell = 1, nshell
68        xa = gx(ishell)
69        ya = gy(ishell)
70        za = gz(ishell)
71        ifrst = shella(ishell)
72        ilast = ifrst + shelln(ishell) - 1
73        itype = shellt(ishell)
74        itype1 = itype + 1
75        istart = ilow(itype1,shellc(ishell)+1)
76        iend = iupp(itype1)
77        iendt = iend
78
79        do jshell = 1,ishell
80          xb = gx(jshell)
81          yb = gy(jshell)
82          zb = gz(jshell)
83          jfsrt = shella(jshell)
84          jlast = jfsrt + shelln(jshell) - 1
85          jtype = shellt(jshell)
86          jtype1 = jtype + 1
87          jstart = ilow(jtype1,shellc(jshell)+1)
88          jend = iupp(jtype1)
89          ijtyp = itype1 + jtype1 - 1
90          ndim = (iend-istart+1)*(jend-jstart+1)
91          iminj = iabs(ishell - jshell)
92          n = (itype + jtype) / 2 + 1
93          xt = xb - xa
94          yt = yb - ya
95          zt = zb - za
96          rsq = xt*xt + yt*yt + zt*zt
97          if (ndim.gt.numtmp)
98     &       print*,'WARNING: exceding array limit in espot'
99          do ii=1,ndim
100             h(ii) = 0.0d0
101          end do
102          jendt = jend
103c
104c         loop over primitive gaussians.
105c
106          do igauss = ifrst, ilast
107            ei = exx(igauss)
108            call fcij(itype,ifrst,igauss,shladf(ishell),ca)
109            do jgauss = jfsrt,jlast
110              iend = iendt
111              jend = jendt
112              ej = exx(jgauss)
113              call fcij(jtype,jfsrt,jgauss,shladf(jshell),cb)
114              ep = ei + ej
115              exptmp = ej*ei*rsq / ep
116              if (exptmp.ge.600.0d0) goto 100
117              expp = 2.0d0*dexp(-exptmp)
118              tx = (ei*xa + ej*xb) / ep
119              ty = (ei*ya + ej*yb) / ep
120              tz = (ei*za + ej*zb) / ep
121              xta = tx - xa
122              yta = ty - ya
123              zta = tz - za
124              xtb = tx - xb
125              ytb = ty - yb
126              ztb = tz - zb
127              call coeffs(coeffx,xta,xtb,itype1,jtype1)
128              call coeffs(coeffy,yta,ytb,itype1,jtype1)
129              call coeffs(coeffz,zta,ztb,itype1,jtype1)
130              if (ndim.gt.numtmp)
131     &           print*,'WARNING: exceding array limit in espot'
132              do ii=1,ndim
133                  f(ii) = 0.0d0
134              end do
135              xct = x - tx
136              yct = y - ty
137              zct = z - tz
138              expont = ep * (xct * xct + yct * yct + zct * zct)
139              call rys(n,expont,tt,zz)
140              epp = 0.5d0 / ep
141              call calct(pre,epp,ijtyp)
142              do ii = 1, n
143                fac1 = (ep + ep)*tt(ii)
144                zf = pi * expp * zz(ii) / ep
145                call twocen(h2x,xct,1.d0,pre,fac1,ijtyp)
146                call twocen(h2y,yct,1.d0,pre,fac1,ijtyp)
147                call twocen(h2z,zct,zf,pre,fac1,ijtyp)
148                call thrcen(h3x,h2x,coeffx,itype1,jtype1)
149                call thrcen(h3y,h2y,coeffy,itype1,jtype1)
150                call thrcen(h3z,h2z,coeffz,itype1,jtype1)
151
152                k = 0
153                do i = istart, iend
154                  do j = jstart, jend
155                    k = k + 1
156                    if (k.gt.numtmp)
157     &                print*,'WARNING: exceding array limit in espot'
158                    f(k) = f(k)-
159     &           (h3x(jx(j)+ix(i))*h3y(jy(j)+iy(i))*h3z(jz(j)+iz(i)))
160                  end do
161                end do
162              end do
163
164              k = 0
165              do i = istart, iend
166                do j = jstart,jend
167                  k = k + 1
168                  if (k.gt.numtmp)
169     &              print*,'WARNING: exceding array limit in espot'
170                  h(k) = h(k) + f(k)*ca(i)*cb(j)
171                end do
172              end do
173              kend = k
174
175100           continue
176              end do
177            end do
178
179            call purdf(itype,jtype,istart,jstart,iend,jend,h)
180c
181c Density matrix contribution
182c
183
184            ii = 0
185            ist = aos(ishell) - 1
186            jst = aos(jshell) - 1
187            do i = istart,iend
188               do j = jstart,jend
189                  ii = ii + 1
190                  if (iminj.eq.0) then
191                     a1or2 = 1.0d0
192                     if (i-j.lt.0) then
193                        pt = p((ist+i-1)*mxorb + (jst+j))
194                     else
195                        pt = p((jst+j-1)*mxorb + (ist+i))
196                     endif
197                  else
198                     a1or2 = 2.0d0
199                     pt = p((jst+j-1)*mxorb + (ist+i))
200                  endif
201                  phelp = pt*a1or2
202                  epot = epot + phelp*h(ii)
203            if (idebug.eq.1) print*,'v(',ist+i,',',jst+j,')=',h(ii)
204               end do
205            end do
206
207            jend = jendt
208            iend = iendt
209        end do
210      end do
211
212c sum in nuclear contribution
213
214c
215      do ii = 1 ,natoms
216          xt = xyz(1,ii) - x
217          yt = xyz(2,ii) - y
218          zt = xyz(3,ii) - z
219          rsq = xt*xt + yt*yt + zt*zt
220          if (rsq.ge.1.0d-8) then
221             if (ipseud.eq.1) then
222                epot = epot + dfloat(ivale(ii))*dsqrt(1.0d0 / rsq)
223             else
224                epot = epot + dfloat(nat(ii))*dsqrt(1.0d0 / rsq)
225             endif
226          endif
227      end do
228
229      if (nexchg.ne.0) then
230          do ii=1,nexchg
231             xt = exchg(1,ii) - x
232             yt = exchg(2,ii) - y
233             zt = exchg(3,ii) - z
234             rsq = xt*xt + yt*yt + zt*zt
235             if (rsq.ge.1.0d-8) then
236                if (iexchg(ii).eq.1) then
237                   epot = epot + dsqrt(1.0d0 / rsq)
238                else
239                   epot = epot - dsqrt(1.0d0 / rsq)
240                endif
241             endif
242          end do
243      endif
244
245      return
246      end
247
248      subroutine espgrdd(npts1,npts2,npts3,idebug,
249     &                   denn,p,xden,yden,zden)
250c THIS IS REALLY espgrd
251      implicit double precision(a-h,o-z)
252      parameter (numatm=2000)
253      parameter (numprm=1600)
254      parameter (numcex=numprm*3)
255      parameter (nmcex2=numcex*numcex)
256      parameter (numtmp=4000)
257      character monstr*100
258      character*2 gstr
259      logical oeerst
260      integer shella,shelln,shellt,shellc,shladf,aos
261      common/b/exx(numcex),
262     &         c1(numcex),c2(numcex),c3(numprm),c4(numprm),c5(numprm),
263     &         shladf(numprm),gx(numprm),gy(numprm),gz(numprm),
264     &         jan(numprm),shella(numprm),shelln(numprm),shellt(numprm),
265     &         shellc(numprm),aos(numprm),nshell,maxtyp
266      common /bounds/ ilow(5,5),iupp(5)
267      common /indxe/  jx(35),jy(35),jz(35),ix(20),iy(20),iz(20)
268      common /strt/   oeerst
269      common /moldat/ natoms, norbs, nelecs,nat(numatm)
270      common /pseudo / ipseud,ivale(numatm)
271      common /coord / xyz(3,numatm)
272      common /orbhlp/ mxorb,iuhf,ispd
273      common /extchg/ exchg(3,3),iexchg(3),nexchg
274      common /plane/  px, py, pz, cx, cy, cz, r(3),v1(3),v2(3),iplat
275      common /grdhlp/ mx3d,mx3d2
276      dimension tt(4),zz(4),pre(45),ca(35),cb(35),
277     &          coeffx(192),coeffy(192),coeffz(192)
278      dimension h(numtmp),f(numtmp),h2x(9),h2y(9),h2z(9),
279     &          h3x(16),h3y(16),h3z(16),v3(3)
280      dimension p(*),denn(*),xden(*),yden(*),zden(*)
281
282      if (oeerst) then
283         call epint
284         call setcon
285         oeerst = .false.
286         do i=1,192
287            coeffx(i) = 0.d0
288         end do
289      endif
290
291
292      call pregrd(npts1,npts2,npts3,xden,yden,zden)
293
294      do kc=1,npts3
295          do ic=1,npts1
296             do jc=1,npts2
297                iadrs = (kc-1)*mx3d2 + ij
298                denn(iadrs) = 0.0d0
299             end do
300          end do
301      end do
302
303      pi = 4.d0*datan(1.d0)
304c
305c     loop over shell pairs.
306c
307      nn = 0
308      iup = 0
309      do ishell = 1, nshell
310
311        xa = gx(ishell)
312        ya = gy(ishell)
313        za = gz(ishell)
314
315        ifrst = shella(ishell)
316        ilast = ifrst + shelln(ishell) - 1
317        itype = shellt(ishell)
318        itype1 = itype + 1
319        istart = ilow(itype1,shellc(ishell)+1)
320        iend = iupp(itype1)
321        iendt = iend
322
323        do jshell = 1,ishell
324
325          xb = gx(jshell)
326          yb = gy(jshell)
327          zb = gz(jshell)
328
329          jfsrt = shella(jshell)
330          jlast = jfsrt + shelln(jshell) - 1
331          jtype = shellt(jshell)
332          jtype1 = jtype + 1
333          jstart = ilow(jtype1,shellc(jshell)+1)
334          jend = iupp(jtype1)
335          ijtyp = itype1 + jtype1 - 1
336          ndim = (iend-istart+1)*(jend-jstart+1)
337          iminj = iabs(ishell - jshell)
338
339          n = (itype + jtype) / 2 + 1
340
341          xt = xb - xa
342          yt = yb - ya
343          zt = zb - za
344
345          rsq = xt*xt + yt*yt + zt*zt
346
347          if (ndim.gt.numtmp)
348     &       print*,'WARNING: exceding array limit in espot'
349
350c deze jendt moet niet in binnenste loop ??
351          jendt = jend
352
353          do kc=1,npts3
354
355             monstr = 'Progress Monitor [1-' // gstr(npts3)//'] ' //
356     &                 gstr(kc)
357             ij = 0
358
359             do ic=1,npts1
360
361                do jc=1,npts2
362
363
364                ij =  ij + 1
365                iadrs = (kc-1)*mx3d2 + ij
366
367                x  = xden(iadrs)
368                y  = yden(iadrs)
369                z  = zden(iadrs)
370
371
372                do ii=1,ndim
373                   h(ii) = 0.0d0
374                end do
375
376
377c         loop over primitive gaussians.
378
379                do igauss = ifrst, ilast
380
381                  ei = exx(igauss)
382                  call fcij(itype,ifrst,igauss,shladf(ishell),ca)
383
384                  do jgauss = jfsrt,jlast
385
386                    iend = iendt
387                    jend = jendt
388
389                    ej = exx(jgauss)
390                    call fcij(jtype,jfsrt,jgauss,shladf(jshell),cb)
391
392                    ep = ei + ej
393                    epr = 1.0d0 / ep
394
395                    exptmp = ej*ei*rsq*epr
396
397                    if (exptmp.ge.600.0d0) goto 100
398
399                    expp = 2.0d0*dexp(-exptmp)
400                    tx = (ei*xa + ej*xb) * epr
401                    ty = (ei*ya + ej*yb) * epr
402                    tz = (ei*za + ej*zb) * epr
403
404                    xta = tx - xa
405                    yta = ty - ya
406                    zta = tz - za
407
408                    xtb = tx - xb
409                    ytb = ty - yb
410                    ztb = tz - zb
411
412                    call coeffs(coeffx,xta,xtb,itype1,jtype1)
413                    call coeffs(coeffy,yta,ytb,itype1,jtype1)
414                    call coeffs(coeffz,zta,ztb,itype1,jtype1)
415
416                    if (ndim.gt.numtmp)
417     &                 print*,'WARNING: exceding array limit in espot'
418
419                    do ii=1,ndim
420                        f(ii) = 0.0d0
421                    end do
422
423                    xct = x - tx
424                    yct = y - ty
425                    zct = z - tz
426
427                    expont = ep * (xct * xct + yct * yct + zct * zct)
428
429                    call rys(n,expont,tt,zz)
430                    epp = 0.5d0 * epr
431                    call calct(pre,epp,ijtyp)
432                    ep2 = ep + ep
433                    pee = pi*expp*epr
434
435                    do ii = 1, n
436                      fac1 = ep2 * tt(ii)
437                      zf = pee * zz(ii)
438                      call twocen(h2x,xct,1.d0,pre,fac1,ijtyp)
439                      call twocen(h2y,yct,1.d0,pre,fac1,ijtyp)
440                      call twocen(h2z,zct,zf,pre,fac1,ijtyp)
441                      call thrcen(h3x,h2x,coeffx,itype1,jtype1)
442                      call thrcen(h3y,h2y,coeffy,itype1,jtype1)
443                      call thrcen(h3z,h2z,coeffz,itype1,jtype1)
444
445                      k = 0
446                      do i = istart, iend
447                        do j = jstart, jend
448                          k = k + 1
449                          if (k.gt.numtmp) print*,
450     &                      'WARNING: exceding array limit in espot'
451                          f(k) = f(k) - (h3x(jx(j)+ix(i))*
452     &                                   h3y(jy(j)+iy(i))*
453     &                                   h3z(jz(j)+iz(i)))
454
455                        end do
456                      end do
457
458                    end do
459
460                    k = 0
461                    do i = istart, iend
462                      do j = jstart,jend
463                        k = k + 1
464
465                        if (k.gt.numtmp) print*,
466     &                    'WARNING: exceding array limit in espot'
467
468                        h(k) = h(k) + f(k)*ca(i)*cb(j)
469
470                      end do
471                    end do
472
473
474100                 continue
475
476                  end do
477                end do
478
479                call purdf(itype,jtype,istart,jstart,iend,jend,h)
480c
481c Density matrix contribution
482c
483
484                ii = 0
485                ist = aos(ishell) - 1
486                jst = aos(jshell) - 1
487
488                do i = istart,iend
489                   do j = jstart,jend
490                      ii = ii + 1
491                      if (iminj.eq.0) then
492                         a1or2 = 1.0d0
493                         if (i-j.lt.0) then
494                            pt = p((ist+i-1)*mxorb + (jst+j))
495                         else
496                            pt = p((jst+j-1)*mxorb + (ist+i))
497                         endif
498                      else
499                         a1or2 = 2.0d0
500                         pt = p((jst+j-1)*mxorb + (ist+i))
501                      endif
502                      phelp = pt*a1or2
503                      denn(iadrs) =  denn(iadrs) + phelp*h(ii)
504
505                      if (idebug.eq.1)
506     &                   print*,'v(',ist+i,',',jst+j,')=',h(ii)
507
508                   end do
509                end do
510
511                jend = jendt
512                iend = iendt
513
514c  end do's kc,ic,jc
515
516                end do
517              end do
518            end do
519
520c  end do's ishell, jshell
521
522        end do
523      end do
524
525c sum in nuclear contribution
526
527      do kc=1,npts3
528          ij = 0
529          do ic=1,npts1
530             do jc=1,npts2
531
532                ij =  ij + 1
533                iadrs = (kc-1)*mx3d2 + ij
534
535                x  = xden(iadrs)
536                y  = yden(iadrs)
537                z  = zden(iadrs)
538
539                do ii = 1,natoms
540
541                    xt = xyz(1,ii) - x
542                    yt = xyz(2,ii) - y
543                    zt = xyz(3,ii) - z
544
545                    rsq = xt*xt + yt*yt + zt*zt
546
547                    if (rsq.ge.1.0d-8) then
548                       if (ipseud.eq.1) then
549                          denn(iadrs) = denn(iadrs) +
550     &                        dfloat(ivale(ii))*dsqrt(1.0d0 / rsq)
551                       else
552                          denn(iadrs) = denn(iadrs) +
553     &                        dfloat(nat(ii))*dsqrt(1.0d0 / rsq)
554                       endif
555                    endif
556                end do
557
558                if (nexchg.ne.0) then
559                    do ii=1,nexchg
560                       xt = exchg(1,ii) - x
561                       yt = exchg(2,ii) - y
562                       zt = exchg(3,ii) - z
563                       rsq = xt*xt + yt*yt + zt*zt
564                       if (rsq.ge.1.0d-8) then
565                          if (iexchg(ii).eq.1) then
566                             denn(iadrs) = denn(iadrs) +
567     &                                       dsqrt(1.0d0 / rsq)
568                          else
569                             denn(iadrs) = denn(iadrs) -
570     &                                       dsqrt(1.0d0 / rsq)
571                          endif
572                       endif
573                    end do
574                endif
575
576             end do
577          end do
578      end do
579
580      return
581      end
582
583      subroutine purdf(itype,jtype,istart,jstart,iend,jend,h)
584c
585c only iend and jend are really nescessary to return
586c
587c iupp 1  4  10 20 35
588c      s  3p 6d 10f 15g
589c
590c ilow 0    1       2        3      4    <- shellt
591c---------------------------------------
592c 0    1(s) 1 (sp)  1 (spd)  1      1
593c 1    1    2 (p)   5       11     21
594c 2    1    1       5 (d)   11     21
595c 3    1    1       5       11 (f) 21
596c 4    1    1       5       11     21 (g)
597c---------------------------------------
598c ^
599c | shellc
600c
601c question does h come filled in shell form ?, eg d is really spd
602c with start at 5 (instead of 1) (we treat d such but not f?)
603c
604c xxxx,yyyy,zzzz *1
605c xxxy ...       *d7
606c xxyy ...       *d5*d7/d3
607c xxyz ...       *d5*d7
608c
609c via fcij ca and cb have these, and therfor h() has them
610c
611      implicit double precision (a-h,o-z)
612      common /slagau/ ihasd,isgau,ido5d,ido7f,ido9g,ihasg
613      common /intcon/ d3,d5,d7,r1,r2,r3,r4,r3ov2,z1,z2,z3,
614     &                g1,g2,g3,g4,g5,g6
615      logical debug
616      dimension h(*),inc(14)
617      debug = .false.
618c
619c     convert  gaussians to pure angular functions.
620c     (6D -> 5D, 10F -> 7F, 15G -> 9G)
621c
622
623      b1 = 1.0d0/d7
624      b2 = d3/(d5*d7)
625      b3 = 1.0d0/(d5*d7)
626
627      idim = iend - istart + 1
628      jdim = jend - jstart + 1
629
630      if (ido5d.eq.1.and.jtype.eq.2) then
631c
632c     The row side of the matrix: pure d
633c
634
635         i1 = 5 - jstart + 1
636
637c     when d   shell jstart = 5 , i1 = 1, jend = 10, jdim = 6
638c     when spd shell jstart = 1 , i1 = 5, jend = 10, jdim = 10
639c     single d shell is stored in H() as:
640c
641c     d1 d2 d3 d4 d5 d6 NOT s px py pz d1 d2 etc
642
643         do i=1,idim
644            dz2 = h(i1+2) - 0.5d0*(h(i1) + h(i1+1))
645            dx2y2 = r3ov2*(h(i1) - h(i1+1))
646            h(i1  ) = dz2
647            h(i1+1) = h(i1+4)
648            h(i1+2) = h(i1+5)
649            h(i1+4) = h(i1+3)
650            h(i1+3) = dx2y2
651            i1 = i1 + jdim
652         end do
653
654      endif
655
656      if (ido7f.eq.1.and.jtype.eq.3) then
657c
658c     The row side af the matrix: pure f
659c
660c     when f  shell jstart = 11 , jend = 20, jdim = 10
661c
662         i1 = 0
663         do i=1,idim
664
665            f0  = h(i1+3) - r2*(h(i1+6) + h(i1+9))
666            f1p = r4*(z1*h(i1+7) - h(i1+1) - z2*h(i1+4))
667            f1m = r4*(z1*h(i1+8) - h(i1+2) - z2*h(i1+5))
668            f2p = r3ov2*(h(i1+6) - h(i1+9))
669            f2m = h(i1+10)
670            f3p = r1*(h(i1+1) - z3*h(i1+4))
671            f3m = r1*(z3*h(i1+5) - h(i1+2))
672
673            h(i1+1) = f0
674            h(i1+2) = f1p
675            h(i1+3) = f1m
676            h(i1+4) = f2p
677            h(i1+5) = f2m
678            h(i1+6) = f3p
679            h(i1+7) = f3m
680
681            i1 = i1 + jdim
682
683         end do
684      endif
685
686      if (ido9g.eq.1.and.jtype.eq.4) then
687c
688c     The row side af the matrix: pure g
689c
690c     when g  shell jstart = 21 , jend = 35, jdim = 15
691c
692         i1 = 0
693         do i=1,idim
694
695          g0  =
696     &       0.375d0*h(i1+1)  + 0.375d0*  h(i1+2) +           h(i1+3)
697     &     -3.0d0*b2*h(i1+11) - 3.0d0*b2*h(i1+12) + 0.75d0*b2*h(i1+10)
698
699          g1p =
700     &  g1*(4.0d0*b1*h(i1+8)  - 3.0d0*b3*h(i1+14) - 3.0d0*b1* h(i1+5))
701
702          g1m =
703     &  g1*(4.0d0*b1*h(i1+9)  - 3.0d0*b3*h(i1+13) - 3.0d0*b1* h(i1+7))
704
705          g2p =
706     &  g2*(6.0d0*b2*h(i1+11) - 6.0d0*b2*h(i1+12) -           h(i1+1)
707     &    +          h(i1+2))
708
709          g2m =
710     &  g3*(6.0d0*b3*h(i1+15) - b1*      h(i1+4)  - b1*       h(i1+6))
711
712          g3p =
713     &  g4*(   b1*   h(i1+5)  - 3.0d0*b3*h(i1+14))
714
715          g3m =
716     &  g4*(3.0d0*b3*h(i1+13) -  b1*     h(i1+7))
717
718          g4p =
719     &  g5*(         h(i1+1)  - 6.0d0*b2*h(i1+10) +           h(i1+2))
720
721          g4m =
722     &  g6*b1*(      h(i1+4)  -          h(i1+6))
723
724            h(i1+1) = g0
725            h(i1+2) = g1p
726            h(i1+3) = g1m
727            h(i1+4) = g2p
728            h(i1+5) = g2m
729            h(i1+6) = g3p
730            h(i1+7) = g3m
731            h(i1+8) = g4p
732            h(i1+9) = g4m
733
734            i1 = i1 + jdim
735
736         end do
737      endif
738c
739c     the row size of the matrix has changed, so get rid of
740c     the superflous functions.
741c
742      if ((ido5d.eq.1.and.jtype.eq.2).or.
743     &    (ido7f.eq.1.and.jtype.eq.3).or.
744     &    (ido9g.eq.1.and.jtype.eq.4)) then
745
746         if (jtype.eq.2) then
747            jendp = 9
748         else if (jtype.eq.3) then
749c     when f   shell jstart = 11 , jend = 17, jpure = jdim = 7, jend = 17
750            jendp = 17
751         else if (jtype.eq.4) then
752c     when g   shell jstart = 21 , jendp = 29, jpure = jdim = 9, jend = 29
753            jendp = 29
754         endif
755         if (debug) print*,'row.old.new ',itype,jtype,jend,jendp
756         jrpure = jendp - jstart + 1
757
758         i1 = 0
759         i2 = 0
760
761         do i=1,idim
762            do j=1,jrpure
763               h(i2+j) = h(i1+j)
764            end do
765            i1 = i1 + jdim
766            i2 = i2 + jrpure
767         end do
768         jend = jendp
769         jdim = jrpure
770
771      endif
772c
773c     transformation at row side complete, start of column side
774c
775
776      if (ido5d.eq.1.and.itype.eq.2) then
777c
778c     The column side of the matrix: pure d
779c
780         i1 = (5-istart)*jdim + 1
781         do i=1,5
782            inc(i) = i*jdim
783         end do
784
785         do j=1,jdim
786            dz2 = h(i1+inc(2)) - 0.5d0*(h(i1) + h(i1+inc(1)))
787            dx2y2 = r3ov2*(h(i1) - h(i1+inc(1)))
788            h(i1       ) = dz2
789            h(i1+inc(1)) = h(i1+inc(4))
790            h(i1+inc(2)) = h(i1+inc(5))
791            h(i1+inc(4)) = h(i1+inc(3))
792            h(i1+inc(3)) = dx2y2
793            i1 = i1 + 1
794         end do
795
796      endif
797
798      if (ido7f.eq.1.and.itype.eq.3) then
799c
800c     The column side af the matrix: pure F
801c
802         i1 = 1
803         do i=1,9
804            inc(i) = i*jdim
805         end do
806
807         do j=1,jdim
808
809            f0  = h(i1+inc(2)) - r2*(h(i1+inc(5)) + h(i1+inc(8)))
810            f1p = r4*(z1*h(i1+inc(6)) - h(i1) - z2*h(i1+inc(3)))
811            f1m = r4*(z1*h(i1+inc(7)) - h(i1+inc(1))-z2*h(i1+inc(4)))
812            f2p = r3ov2*(h(i1+inc(5)) - h(i1+inc(8)))
813            f2m = h(i1+inc(9))
814            f3p = r1*(h(i1) - z3*h(i1+inc(3)))
815            f3m = r1*(z3*h(i1+inc(4)) - h(i1+inc(1)))
816
817            h(i1       ) = f0
818            h(i1+inc(1)) = f1p
819            h(i1+inc(2)) = f1m
820            h(i1+inc(3)) = f2p
821            h(i1+inc(4)) = f2m
822            h(i1+inc(5)) = f3p
823            h(i1+inc(6)) = f3m
824
825            i1 = i1 + 1
826         end do
827
828      endif
829
830      if (ido9g.eq.1.and.itype.eq.4) then
831c
832c     The column side af the matrix: pure G
833c
834         i1 = 1
835         do i=1,14
836            inc(i) = i*jdim
837         end do
838
839         do j=1,jdim
840
841          g0  =
842     &          0.375d0*h(i1)         + 0.375d0*  h(i1+inc(1))  +
843     &                  h(i1+inc(2))  - 3.0d0*b2* h(i1+inc(10)) -
844     &         3.0d0*b2*h(i1+inc(11)) + 0.75d0*b2*h(i1+inc(9))
845
846          g1p =
847     &     g1*(4.0d0*b1*h(i1+inc(7))  - 3.0d0*b3* h(i1+inc(13)) -
848     &         3.0d0*b1*h(i1+inc(4)))
849
850          g1m =
851     &     g1*(4.0d0*b1*h(i1+inc(8))  - 3.0d0*b3* h(i1+inc(12)) -
852     &         3.0d0*b1*h(i1+inc(6)))
853
854          g2p =
855     &     g2*(6.0d0*b2*h(i1+inc(10)) - 6.0d0*b2* h(i1+inc(11)) -
856     &                  h(i1)         +           h(i1+inc(1)))
857
858          g2m =
859     &     g3*(6.0d0*b3*h(i1+inc(14)) - b1*       h(i1+inc(3))  -
860     &               b1*h(i1+inc(5)))
861
862          g3p =
863     &     g4*(b1*      h(i1+inc(4))  - 3.0d0*b3* h(i1+inc(13)))
864
865          g3m =
866     &     g4*(3.0d0*b3*h(i1+inc(12)) - b1*       h(i1+inc(6)))
867
868          g4p =
869     &      g5*(        h(i1)         - 6.0d0*b2* h(i1+inc(9))  +
870     &                  h(i1+inc(1)))
871
872          g4m =
873     &      g6*b1*(     h(i1+inc(3))  -           h(i1+inc(5)) )
874
875            h(i1       ) = g0
876            h(i1+inc(1)) = g1p
877            h(i1+inc(2)) = g1m
878            h(i1+inc(3)) = g2p
879            h(i1+inc(4)) = g2m
880            h(i1+inc(5)) = g3p
881            h(i1+inc(6)) = g3m
882            h(i1+inc(7)) = g4p
883            h(i1+inc(8)) = g4m
884
885            i1 = i1 + 1
886         end do
887
888      endif
889      if ((ido5d.eq.1.and.itype.eq.2) .or.
890     &    (ido7f.eq.1.and.itype.eq.3) .or.
891     &    (ido9g.eq.1.and.itype.eq.4)) then
892
893         if (debug) print*,'colm.old ',itype,jtype,iend
894         if (itype.eq.2) then
895            iend = 9
896         else if (itype.eq.3) then
897            iend = 17
898         else if (itype.eq.4) then
899            iend = 29
900         endif
901         if (debug) print*,'colm.new ',itype,jtype,iend
902
903      endif
904
905      return
906      end
907
908      subroutine setcon
909      implicit double precision (a-h,o-z)
910      common /intcon/ d3,d5,d7,r1,r2,r3,r4,r3ov2,z1,z2,z3,
911     &                g1,g2,g3,g4,g5,g6
912      d3 = dsqrt(3.0d0)
913      d5 = dsqrt(5.0d0)
914      d7 = dsqrt(7.0d0)
915      z1 = 4.0d0/d5
916      z2 = 1.0d0/d5
917      z3 = 3.0d0/d5
918      r1 = 0.5d0*dsqrt(5.0d0/2.0d0)
919      r2 = 1.5d0/d5
920      r3ov2 = 0.5d0*dsqrt(3.0d0)
921      r3 = r3ov2
922      r4 = 0.5d0*dsqrt(3.0d0/2.0d0)
923
924      g1 = dsqrt(5.0d0/8.0d0)
925      g2 = dsqrt(5.0d0/16.0d0)
926      g3 = dsqrt(5.0d0/4.0d0)
927      g4 = dsqrt(35.0d0/8.0d0)
928      g5 = dsqrt(35.0d0/64.0d0)
929      g6 = dsqrt(35.0d0/4.0d0)
930
931      return
932      end
933
934      subroutine dipold(ntt,nor,dmao,focc,focb,vectrs,vectrb,p)
935      implicit double precision (a-h,o-z), integer (i-n)
936      parameter (numatm=2000)
937      parameter (numprm=1600)
938      parameter (numcex=numprm*3)
939
940c     dipole integrals
941
942      logical debug
943      common /moldat/ natoms, norbs, nelecs,nat(numatm)
944      common /charge/ dipo(3),ihasq,ihsdp,iqon,idipon
945      common /bounds/ ilow(5,5),iupp(5)
946      common /orbhlp/ mxorb,iuhf,ispd
947      common /coord / xyz(3,numatm)
948
949      integer shella,shelln,shellt,shellc,shladf,aos
950      common/b/exx(numcex),
951     &         c1(numcex),c2(numcex),c3(numprm),c4(numprm),c5(numprm),
952     &         shladf(numprm),gx(numprm),gy(numprm),gz(numprm),
953     &         jan(numprm),shella(numprm),shelln(numprm),shellt(numprm),
954     &         shellc(numprm),aos(numprm),nshell,maxtyp
955
956      common /indxe/  jx(35),jy(35),jz(35),ix(20),iy(20),iz(20)
957      common /slagau/ ihasd,isgau,ido5d,ido7f,ido9g,ihasg
958      double precision imprd
959      integer genaos
960
961      dimension focc(*),focb(*),vectrs(*),vectrb(*),nor(*)
962      dimension ca(35),cb(35),sx(36),sy(36),sz(36),dd(3,100),
963     &    ccx(330),ccy(330),ccz(330),p(*),dmao(ntt,3),
964     &    d12(100),dipon(3)
965
966      lind(i,j) = ((max0(i,j)*(max0(i,j)-1)/2)+min0(i,j))
967
968c iupp 1  4  10 20 35
969c      s  3p 6d 10f 15g
970c
971c ilow 0    1       2        3      4    <- shellt
972c---------------------------------------
973c 0    1(s) 1 (sp)  1 (spd)  1      1
974c 1    1    2 (p)   5       11     21
975c 2    1    1       5 (d)   11     21
976c 3    1    1       5       11 (f) 21
977c 4    1    1       5       11     21 (g)
978c---------------------------------------
979c ^
980c | shellc
981c
982
983      debug = .false.
984      cf = 2.5417463d0
985      idum = genaos(.false.)
986
987      pi = 4.d0*datan(1.d0)
988
989      do i=1,3
990         dipo(i) = 0.d0
991         dipon(i) = 0.d0
992      end do
993
994      do i=1,natoms
995         do j=1,3
996            dipon(j) = dipon(j) + dble(nat(i))*xyz(j,i)
997         end do
998      end do
999
1000      call denini
1001
1002      do ishell=1,nshell
1003
1004         xi = gx(ishell)
1005         yi = gy(ishell)
1006         zi = gz(ishell)
1007
1008         ifrst = shella(ishell)
1009         ilast = ifrst + shelln(ishell) - 1
1010         itype = shellt(ishell)
1011
1012         itype1 = itype + 1
1013         lidim = itype1 + 2
1014
1015         istart = ilow(itype1,shellc(ishell)+1)
1016         iend   = iupp(itype1)
1017         iendt  = iend
1018
1019         do jshell=1,ishell
1020
1021            xj = gx(jshell)
1022            yj = gy(jshell)
1023            zj = gz(jshell)
1024
1025            jfrst = shella(jshell)
1026            jlast = jfrst + shelln(jshell) - 1
1027            jtype = shellt(jshell)
1028
1029            jtype1 = jtype+1
1030            ljdim  = jtype+2
1031
1032            jstart = ilow(jtype1,shellc(jshell)+1)
1033            jend   = iupp(jtype1)
1034            jendt  = jend
1035            lpdim  = lidim + ljdim
1036
1037            xt = xi - xj
1038            yt = yi - yj
1039            zt = zi - zj
1040
1041            r2 = xt*xt + yt*yt + zt*zt
1042
1043            do l=1,100
1044               do k=1,3
1045                  dd(k,l) = 0.d0
1046               end do
1047            end do
1048
1049            irp = iendt - istart + 1
1050            jrp = jendt - jstart + 1
1051            lentqp = irp*jrp
1052
1053            if (ishell.eq.jshell) lentqp = (irp*(irp+1))/2
1054
1055            do igauss=ifrst,ilast
1056
1057               ei = exx(igauss)
1058               twoei = ei + ei
1059
1060               call fcij(itype,ifrst,igauss,shladf(ishell),ca)
1061
1062               do jgauss=jfrst,jlast
1063
1064                  ej = exx(jgauss)
1065                  twoej = ej + ej
1066
1067                  call fcij(jtype,jfrst,jgauss,shladf(jshell),cb)
1068
1069                  ep = ei + ej
1070
1071                  tx = (ei*xi + ej*xj) / ep
1072                  ty = (ei*yi + ej*yj) / ep
1073                  tz = (ei*zi + ej*zj) / ep
1074
1075                  xit = tx - xi
1076                  yit = ty - yi
1077                  zit = tz - zi
1078
1079                  xjt = tx - xj
1080                  yjt = ty - yj
1081                  zjt = tz - zj
1082
1083                  expont = (ej*ei*r2) / ep
1084                  pexp   = 0.d0
1085
1086                  if (expont.lt.600.d0) pexp = dexp(-expont)
1087
1088                  call dipint(lidim,ljdim,lpdim,
1089     &                  xit,xjt,ccx, yit,yjt,ccy, zit,zjt,ccz,
1090     &                  ei,ej,pi,pexp,sx,sy,sz)
1091
1092                  call dipao((ishell.eq.jshell),istart,iend,jstart,jend,
1093     &                 lidim,ljdim,ca,cb,
1094     &                 sx,sy,sz,xi,yi,zi,dd)
1095
1096               end do
1097            end do
1098
1099            do l=1,3
1100               call purdf(itype,jtype,istart,jstart,iend,jend,dd(l,1))
1101            end do
1102
1103            ii = 0
1104            ist = aos(ishell) - 1
1105            jst = aos(jshell) - 1
1106
1107            do i=istart,iend
1108               jnd = jend
1109               if (ishell.eq.jshell) jnd = i
1110
1111               do j=jstart,jnd
1112                  ii = ii + 1
1113                  d12(ii) = p((ist+i-1)*mxorb + (jst+j))
1114                  if (ist+i.ne.jst+j) d12(ii) = 2.0d0*d12(ii)
1115                  ll = lind(ist+i,jst+j)
1116                  do l=1,3
1117                     dmao(ll,l) =  dd(l,ii)
1118                  end do
1119               end do
1120
1121            end do
1122
1123            do l=1,3
1124               dipo(l) = dipo(l) - imprd(lentqp,d12,dd,l)
1125            end do
1126
1127         end do
1128      end do
1129
1130      if (debug) then
1131         write(6,'(a)') 'Dipole moment vector (Debye)'
1132         write(6,'(a,3f10.4)') 'elec dipole ',(cf*dipo(i),i=1,3)
1133         write(6,'(a,3f10.4)') 'nuc  dipole ',(cf*dipon(i),i=1,3)
1134
1135         dipt = 0.d0
1136         do l=1,3
1137            dipo(l) = dipo(l) + dipon(l)
1138            dipt = dipt + dipo(l)*dipo(l)
1139         end do
1140         dipt = dsqrt(dipt)
1141         write(6,'(a,3f10.4)') 'tot  dipole ',(cf*dipo(i),i=1,3)
1142         write(6,'(a,f10.4)') 'total dipole scalar ',cf*dipt
1143      endif
1144
1145      call setnor(nocc,norbs,nor,focc)
1146      call boys(nocc,nor,dmao,vectrs)
1147
1148      if (iuhf.ne.0) then
1149          call setnor(nocc,norbs,nor,focb)
1150          call boys(nocc,nor,dmao,vectrb)
1151      endif
1152
1153      return
1154      end
1155
1156      subroutine dipint(l1max,l2max,l12mxp,
1157     &                  ax,bx,ccx, ay,by,ccy, az,bz,ccz,
1158     &                  as,bs,pi,pexp,sx,sy,sz)
1159      implicit double precision (a-h,o-z), integer (i-n)
1160      dimension ccx(l12mxp,l2max,l1max)
1161      dimension ccy(l12mxp,l2max,l1max)
1162      dimension ccz(l12mxp,l2max,l1max)
1163      dimension sx(l2max,l1max),sy(l2max,l1max),sz(l2max,l1max)
1164      dimension s1(21)
1165
1166      do l1=1,l1max
1167         do l2=1,l2max
1168            do l12=1,l12mxp
1169               ccx(l12,l2,l1) = 0.d0
1170               ccy(l12,l2,l1) = 0.d0
1171               ccz(l12,l2,l1) = 0.d0
1172            end do
1173         end do
1174      end do
1175
1176      ccx(2,1,1) = 1.d0
1177      ccy(2,1,1) = 1.d0
1178      ccz(2,1,1) = 1.d0
1179
1180      do l1=1,l1max
1181         do l2=1,l2max
1182            iwmax = l1 + l2
1183            do iw=2,iwmax
1184
1185               if (l1.gt.1) then
1186
1187                   ccx(iw,l2,l1) =
1188     &                  ax*ccx(iw,l2,l1-1) + ccx(iw-1,l2,l1-1)
1189                   ccy(iw,l2,l1) =
1190     &                  ay*ccy(iw,l2,l1-1) + ccy(iw-1,l2,l1-1)
1191                   ccz(iw,l2,l1) =
1192     &                  az*ccz(iw,l2,l1-1) + ccz(iw-1,l2,l1-1)
1193
1194               elseif (l2.gt.1) then
1195
1196                   ccx(iw,l2,l1) =
1197     &                  bx*ccx(iw,l2-1,l1) + ccx(iw-1,l2-1,l1)
1198                   ccy(iw,l2,l1) =
1199     &                  by*ccy(iw,l2-1,l1) + ccy(iw-1,l2-1,l1)
1200                   ccz(iw,l2,l1) =
1201     &                  bz*ccz(iw,l2-1,l1) + ccz(iw-1,l2-1,l1)
1202
1203               endif
1204
1205            end do
1206         end do
1207      end do
1208
1209      temp = 1.d0/(2.d0*(as+bs))
1210      s1(1) = dsqrt(pi/(as+bs))
1211
1212      lim = l1max + l2max - 1
1213
1214      if (lim.ge.2) then
1215
1216         fact = 1.d0
1217
1218         do i=2,lim
1219            if (mod(i,2).eq.0) then
1220               s1(i) = 0.d0
1221            else
1222               s1(i) = s1(i-2)*temp*fact
1223               fact   = fact + 2.d0
1224            endif
1225         end do
1226
1227      endif
1228
1229      do l1=1,l1max
1230         do l2=1,l2max
1231
1232            x = 0.d0
1233            y = 0.d0
1234            z = 0.d0
1235            iwmax = l1 + l2 - 1
1236
1237            do iw=1,iwmax
1238               x = x + ccx(iw+1,l2,l1)*s1(iw)
1239               y = y + ccy(iw+1,l2,l1)*s1(iw)
1240               z = z + ccz(iw+1,l2,l1)*s1(iw)
1241            end do
1242
1243            sx(l2,l1) = x
1244            sy(l2,l1) = y
1245            sz(l2,l1) = z*pexp
1246
1247         end do
1248      end do
1249
1250      return
1251      end
1252
1253      subroutine dipao(ieqj,istart,iend,jstart,jend,
1254     &                 l1max,l2max,ca,cb,sx,sy,sz,xa,ya,za,dd)
1255      implicit double precision (a-h,o-z), integer (i-n)
1256      logical ieqj
1257      common /indxe/  lx(35),ly(35),lz(35),kx(20),ky(20),kz(20)
1258      dimension ca(*),cb(*),dd(3,*)
1259      dimension sx(l2max,l1max),sy(l2max,l1max),sz(l2max,l1max)
1260      dimension rxyz(3)
1261
1262      ii = 0
1263
1264      do i=istart,iend
1265
1266         ix = lx(i)
1267         iy = ly(i)
1268         iz = lz(i)
1269
1270         jnd = jend
1271
1272         if (ieqj) jnd = i
1273
1274         do j=jstart,jnd
1275
1276            jx = lx(j)
1277            jy = ly(j)
1278            jz = lz(j)
1279
1280            ii = ii + 1
1281            coef = ca(i)*cb(j)
1282
1283            rxyz(1) = (xa*sx(jx,ix) + sx(jx,ix+1))
1284     &                *sy(jy,iy)*sz(jz,iz)
1285            rxyz(2) = (ya*sy(jy,iy) + sy(jy,iy+1))
1286     &                *sx(jx,ix)*sz(jz,iz)
1287            rxyz(3) = (za*sz(jz,iz) + sz(jz,iz+1))
1288     &                *sx(jx,ix)*sy(jy,iy)
1289
1290            do l=1,3
1291               dd(l,ii) = dd(l,ii) + coef*rxyz(l)
1292            end do
1293
1294         end do
1295      end do
1296
1297      return
1298      end
1299
1300      double precision function imprd(n,a,b,l)
1301      implicit double precision (a-h,o-z), integer (i-n)
1302c     improduct of a and b.
1303      dimension a(*), b(3,*)
1304
1305      imprd = 0.d0
1306
1307      if (n.lt.1) return
1308
1309      do i=1,n
1310          imprd = imprd + a(i)*b(l,i)
1311      end do
1312
1313      return
1314      end
1315
1316      subroutine boyd(norb,nor,dmao,cc,nbasis,ntt,cl,
1317     &                iord,iir,rij,qpix,qpjx)
1318      implicit double precision (a-h,o-z), integer (i-n)
1319      parameter (numatm=2000)
1320      logical debug
1321
1322C     Boys Localization.  After QCPE program 354
1323C
1324C     NOrb   ... Number of orbitals to localize.
1325C     CC     ... Input and output orbitals (NBasis,).
1326C     CL     ... Scratch matrix (NBasis,NOrb).
1327
1328      common /orbhlp/ mxorb,iuhf,ispd
1329      common /moldat/ natoms, norbs, nelecs,nat(numatm)
1330      dimension cc(*), cl(nbasis,*), dmao(ntt,3), nor(nbasis,nbasis),
1331     &          iord(*), iir(*), rij(ntt,3), qpix(*), qpjx(*)
1332
1333      lind(i,j) = ((max0(i,j)*(max0(i,j)-1)/2)+min0(i,j))
1334      kind(j,i) = ((i-1)*mxorb + j)
1335
1336c essential that in above definition i and j have swapped places !!!
1337
1338
1339      debug = .false.
1340      cf = 2.5417463d0
1341
1342      do itry=1,3
1343
1344          do i=1,norb
1345              ii = i
1346
1347              do j=1,i
1348                  jj = j
1349                  sumx = 0.0d0
1350                  sumy = 0.0d0
1351                  sumz = 0.0d0
1352
1353                  do k=1,nbasis
1354                      do l=1,nbasis
1355                          cpr = cc(kind(k,ii))*cc(kind(l,jj))
1356                          sumx = sumx + cpr*dmao(lind(k,l),1)
1357                          sumy = sumy + cpr*dmao(lind(k,l),2)
1358                          sumz = sumz + cpr*dmao(lind(k,l),3)
1359                      end do
1360                  end do
1361
1362                  rij(lind(i,j),1) = sumx
1363                  rij(lind(i,j),2) = sumy
1364                  rij(lind(i,j),3) = sumz
1365
1366              end do
1367
1368          end do
1369
1370
1371      call caltra(norb,ntt,rij,tracia)
1372
1373      if (debug) write(6,'(a,f15.8)')
1374     &          ' Initial TraceA=',tracia*cf*cf
1375
1376      do i=1,nbasis
1377          do j=1,nbasis
1378             cl(i,j) = 0.0d0
1379          end do
1380      end do
1381
1382      do i=1,norb
1383          cl(i,i) = 1.0d0
1384      end do
1385
1386      iter = 0
1387      shift = datan(1.0d0)
1388
1389      change = 1.0d0
1390
1391       do while (iter.lt.150.and.change.gt.1.d-6)
1392
1393         if (change.ge.1.d-10) then
1394            change = 0.0d0
1395
1396            iter = iter + 1
1397
1398            do i=1,norb
1399                iir(i) = i
1400            end do
1401
1402            nnn = norb
1403
1404            do i=1,norb
1405                iii = irnumb(nnn)
1406                iord(i) = iir(iii)
1407                iir(iii) = iir(nnn)
1408                nnn = nnn - 1
1409            end do
1410
1411         endif
1412
1413         do iii=1,norb
1414
1415           i  = iord(iii)
1416           if (nor(i,i).eq.0) then
1417
1418             li = lind(i,i)
1419
1420             jm = 1
1421             rm = 0.0d0
1422             tm = 0.0d0
1423             sm = 0.0d0
1424             cm = 1.0d0
1425
1426             do j=1,norb
1427
1428               if (i.ne.j.and.nor(i,j).eq.0) then
1429
1430                 lj = lind(j,j)
1431                 ij = lind(i,j)
1432                 t  = 0.0d0
1433                 tx = 0.0d0
1434
1435                 do kk=1,3
1436                     t = t + 4.0d0*rij(ij,kk)**2 - rij(li,kk)**2
1437     &                  - rij(lj,kk)**2
1438     &                  + 2.0d0*rij(li,kk)*rij(lj,kk)
1439                     tx = tx + rij(ij,kk)*(rij(lj,kk)-rij(li,kk))
1440                 end do
1441
1442                 if (dabs(t).gt.1.d-10.or.dabs(tx).gt.1.d-10) then
1443
1444                    tx = 4.0d0*tx
1445                    t = datan2(tx,t)/4.0d0
1446                    sign = 1.0d0
1447                    if (t.gt.0.0d0) sign = -1.0d0
1448                    t = t + sign*shift
1449                    itim = 0
1450
1451                    s = dsin(t)
1452                    c = dcos(t)
1453                    rin = 0.0d0
1454                    itim = itim + 1
1455
1456                    do kk=1,3
1457                        qpi = c*c*rij(li,kk) + s*s*rij(lj,kk)
1458     &                        + 2.0d0*c*s*rij(ij,kk)
1459                        qpj = c*c*rij(lj,kk)+s*s*rij(li,kk)
1460     &                        - 2.0d0*c*s*rij(ij,kk)
1461                        rin = rin + qpi*qpi+qpj*qpj - rij(li,kk)**2
1462     &                            - rij(lj,kk)**2
1463                    end do
1464
1465                    ttest = dabs(t) - shift
1466
1467                    if (dabs(t).le.1.d-8.or.dabs(ttest).lt.1.d-8.or.
1468     &                  rin.ge.(-1.d-8)) then
1469                        if (rin.gt.rm) then
1470                            rm = rin
1471                            jm = j
1472                            sm = s
1473                            cm = c
1474                            tm = t
1475                        endif
1476                    else
1477
1478                       if (debug) then
1479                         write(6,'(a,i4,a,i4,a,f15.8,a,f15.8)')
1480     &               ' BoyLoc:  No rotation increases integrals for I=',
1481     &                    i,' J=',j,' Theta=',t,' Change=',rin
1482                       endif
1483
1484                    endif
1485
1486                 endif
1487
1488               endif
1489
1490             end do
1491
1492c end  loop over j
1493
1494             rin = rm
1495             j = jm
1496             s = sm
1497             c = cm
1498             t = tm
1499
1500             ij = lind(i,j)
1501             lj = lind(j,j)
1502
1503             if (nor(i,j).eq.0) then
1504
1505                change = change + t*t
1506
1507                do kk=1,3
1508                    qpi = c*c*rij(li,kk) + s*s*rij(lj,kk)
1509     &                                   + 2.0d0*c*s*rij(ij,kk)
1510                    qpj = c*c*rij(lj,kk) + s*s*rij(li,kk)
1511     &                                   - 2.0d0*c*s*rij(ij,kk)
1512                    qpij = (c*c-s*s)*rij(ij,kk) +
1513     &                      c*s*(rij(lj,kk)-rij(li,kk))
1514
1515                    do k=1,norb
1516                        if (i.ne.k.and.j.ne.k) then
1517                           ik = lind(i,k)
1518                           jk = lind(j,k)
1519                           qpix(k) = c*rij(ik,kk) + s*rij(jk,kk)
1520                           qpjx(k) = c*rij(jk,kk) - s*rij(ik,kk)
1521                        endif
1522                    end do
1523
1524                    do k=1,norb
1525                        if (i.ne.k.and.j.ne.k) then
1526                           ik = lind(i,k)
1527                           jk = lind(j,k)
1528                           rij(ik,kk) = qpix(k)
1529                           rij(jk,kk) = qpjx(k)
1530                        endif
1531                    end do
1532
1533                    rin = rin + qpi + qpj - rij(li,kk) - rij(lj,kk)
1534                    rij(li,kk) = qpi
1535                    rij(lj,kk) = qpj
1536                    rij(ij,kk) = qpij
1537
1538                end do
1539
1540                do k=1,norb
1541                    c1 =  c*cl(k,i) + s*cl(k,j)
1542                    c2 = -s*cl(k,i) + c*cl(k,j)
1543                    cl(k,i) = c1
1544                    cl(k,j) = c2
1545                end do
1546
1547             endif
1548
1549           endif
1550         end do
1551
1552c end loop over iii
1553
1554         change = dsqrt(2.0d0*change/dble(norb*(norb-1)))
1555
1556
1557         if (debug) then
1558
1559            write(6,'(a,i4,a,f15.6)')
1560     &              ' BoyLoc:  Iteration',iter,' Change=',change
1561
1562
1563         endif
1564
1565       end do
1566
1567      end do
1568
1569      if (iter.gt.150) then
1570            write(6,'(a)')
1571     &        ' Localization failed. after 3 tries of 150 iterations.'
1572      else
1573         if (change.le.1.d-6) then
1574           write(6,'(a,i4,a)')
1575     &          ' Localization complete after',iter,' iterations.'
1576         endif
1577      endif
1578
1579      do i=1,norb
1580
1581          do k=1,nbasis
1582             qpix(k) = 0.0d0
1583          end do
1584
1585          do j=1,norb
1586              jj = j
1587
1588              do k=1,nbasis
1589                  qpix(k) = qpix(k) + cc(kind(k,jj))*cl(j,i)
1590              end do
1591
1592          end do
1593
1594          do k=1,nbasis
1595             cl(k,i) = qpix(k)
1596          end do
1597
1598      end do
1599
1600c     put rotated orbitals back
1601
1602      do i=1,norb
1603          ii = i
1604
1605          do k=1,nbasis
1606             cc(kind(k,ii)) = cl(k,i)
1607          end do
1608
1609      end do
1610
1611      if (debug) then
1612
1613         print*,"rotated orbitals"
1614         print*," "
1615         call prev(cc,nbasis,nbasis,mxorb)
1616
1617         call caltra(norb,ntt,rij,tracfa)
1618         write(6,'(a,f15.8)') ' Final TraceA=',tracfa
1619
1620      endif
1621
1622      call caldrv(norb,ntt,rij,dmax,imax,jmax)
1623
1624      if (debug) then
1625          write(6,'(a,i3,a,i3,a,f15.8)')
1626     &    ' Largest derivative:  DDelta(',imax,',',jmax,
1627     &    ')=', dmax
1628      endif
1629
1630      if (dmax.gt.1.d-8) then
1631          write(6,'(a,i3,a,i3,a,f15.8)')
1632     &    ' Internal error in BoyLoc, DDelta(',imax,',',jmax,
1633     &    ')=',dmax
1634      endif
1635
1636      return
1637      end
1638
1639      integer function irnumb(max)
1640      implicit double precision (a-h,o-z), integer (i-n)
1641
1642      irnumb = int(random()*dble(max)+1.0d0)
1643
1644      return
1645      end
1646
1647      subroutine caltra(norb,ntt,rij,traca)
1648      implicit double precision (a-h,o-z), integer (i-n)
1649      dimension rij(ntt,3)
1650
1651      lind(i,j) = ((max0(i,j)*(max0(i,j)-1)/2)+min0(i,j))
1652
1653      traca = 0.0d0
1654
1655      ii = 0
1656      do i=1,norb
1657         ii = ii + i
1658         if (i.ne.1) then
1659         traca = traca + rij(ii,1)**2 + rij(ii,2)**2 + rij(ii,3)**2
1660         endif
1661      end do
1662
1663      return
1664      end
1665
1666      subroutine caldrv(norb,ntt,rij,dmax,imax,jmax)
1667      implicit double precision (a-h,o-z), integer (i-n)
1668      dimension rij(ntt,3)
1669
1670      lind(i,j) = ((max0(i,j)*(max0(i,j)-1)/2)+min0(i,j))
1671
1672      dmax = -1.0d0
1673
1674      do i=2,norb
1675
1676          li = lind(i,i)
1677
1678          do j=1,i-1
1679
1680              lj = lind(j,j)
1681              ij = lind(i,j)
1682              dt = 0.0d0
1683
1684              do kk=1,3
1685                  dt = dt + rij(ij,kk)*(rij(li,kk)-rij(lj,kk))
1686              end do
1687
1688              if (dabs(dt).gt.dmax) then
1689                  dmax = dabs(dt)
1690                  imax = i
1691                  jmax = j
1692              endif
1693
1694           end do
1695      end do
1696
1697      dmax = 2.0d0*dmax
1698
1699      return
1700      end
1701
1702      double precision function random()
1703      implicit double precision (a-h,o-z), integer (i-n)
1704      parameter (mplier=16807,modlus=2147483647,mobymp=127773,
1705     &           momdmp=2836)
1706      common  /seed/jseed,ifrst,nextn
1707c     mseed comes from alloc.c
1708      integer hvlue, lvlue, testv, nextn, mseed
1709
1710      if (ifrst.eq.0) then
1711        jseed = mseed()
1712        nextn = jseed
1713        ifrst = 1
1714      endif
1715
1716      hvlue = nextn / mobymp
1717      lvlue = mod(nextn, mobymp)
1718      testv = mplier*lvlue - momdmp*hvlue
1719
1720      if (testv.gt.0) then
1721        nextn = testv
1722      else
1723        nextn = testv + modlus
1724      endif
1725
1726      random = dble(nextn)/dble(modlus)
1727
1728      return
1729      end
1730
1731      integer function ncoree()
1732      implicit double precision (a-h,o-z), integer (i-n)
1733      parameter (numatm=2000)
1734      common /moldat/ natoms, norbs, nelecs,nat(numatm)
1735      dimension ncore(103)
1736      data ncore /2*0,2*1,6*1,2*5,6*5,2*9,10*9,9,5*14,
1737     &            2*18,10*18,18,5*23,2*27,14*27,10*34,
1738     &            34,5*39,2*43,14*43,50/
1739
1740      ncoree = 0
1741
1742      do i=1,natoms
1743          ncoree = ncoree + ncore(nat(i))
1744      end do
1745
1746      return
1747      end
1748
1749      subroutine setnor(nocc,nbasis,nor,focc)
1750      implicit double precision (a-h,o-z), integer (i-n)
1751      dimension focc(*),nor(nbasis,nbasis)
1752
1753      nocc = 0
1754      do i=1,nbasis
1755         if (focc(i).ne.0.d0) then
1756            nocc = nocc + 1
1757         endif
1758      end do
1759
1760      do i=1,nbasis
1761         do j=1,nbasis
1762            nor(i,j) = 1
1763         end do
1764      end do
1765
1766      mcore = ncoree()
1767      do i=mcore+1,nocc
1768         do j=mcore+1,nocc
1769            nor(i,j) = 0
1770         end do
1771      end do
1772
1773      return
1774      end
1775
1776      subroutine pregrd(npts1,npts2,npts3,xden,yden,zden)
1777      implicit double precision(a-h,o-z)
1778      common /plane/  px, py, pz, cx, cy, cz, r(3),v1(3),v2(3),iplat
1779      common /grdhlp/ mx3d,mx3d2
1780      dimension v3(3),xden(*),yden(*),zden(*)
1781
1782      v3(1) = cx
1783      v3(2) = cy
1784      v3(3) = cz
1785
1786      call vnrm(v3)
1787
1788      radius1 = 0.5d0 * r(1)
1789      radius2 = 0.5d0 * r(2)
1790      radius3 = 0.5d0 * r(3)
1791      rf1     = r(1) / dble(npts1-1)
1792      rf2     = r(2) / dble(npts2-1)
1793      rf3     = r(3) / dble(npts3-1)
1794
1795      do kc=1,npts3
1796          z2 = -radius3+dble(kc-1)*rf3
1797          do ic=1,npts1
1798             x2 = -radius1+dble(ic-1)*rf1
1799             do jc=1,npts2
1800                y2 = -radius2+dble(jc-1)*rf2
1801
1802                iadrs = (kc-1)*mx3d2 + ij
1803
1804                xden(iadrs)  = v1(1)*x2 + v2(1)*y2 + px - v3(1)*z2
1805                yden(iadrs)  = v1(2)*x2 + v2(2)*y2 + py - v3(2)*z2
1806                zden(iadrs)  = v1(3)*x2 + v2(3)*y2 + pz - v3(3)*z2
1807
1808             end do
1809          end do
1810      end do
1811
1812      return
1813      end
1814
1815