1c $Id$
2c
3c (atw) Modified to directly create the distributed AO-density
4c       from the atomic density blocks
5c (7/3/94)
6c
7c Routines changed: (see guess_dens1.F)
8c                    guess_dens
9c                    guess_mem
10c                    denat
11c                    datoms
12c                    creded
13c                    pdfded
14c
15c
16
17
18
19
20
21
22
23#ifndef ADRIANS_CRAP
24      subroutine guess_dens(rtdb, basis, g_dens, g_old_dens, iprint)
25      implicit none
26c
27#include "mafdecls.h"
28#include "global.h"
29#include "tcgmsg.h"
30c
31c     arguments
32c
33      integer rtdb, basis, g_dens, g_old_dens, iprint
34c
35c     local variables
36c
37      integer nshell, iproc, nproc, mem1, max1e
38      integer natoms, nbf, nprim, maxprim, max_l, max_sh_bf, max_at_bf
39c
40      integer l_dens, l_scr
41      integer k_dens, k_scr
42      logical status
43c
44      integer iread, iwrite
45      common /iofile/ iread,iwrite
46c
47c     Global array g_dens     = Density Matrix
48c     Global array g_old_dens = old Density Matrix
49c
50c
51      call ga_zero(g_dens)
52      call ga_zero(g_old_dens)
53c
54c     Get info about the basis set
55c
56      call gto_info(basis, natoms, nshell, nbf, nprim, maxprim,
57     $     max_l, max_sh_bf, max_at_bf)
58c
59      iwrite = 6
60c
61      iproc = nodeid()
62      nproc = nnodes()
63c     iproc = ga_nodeid()
64c     nproc = ga_nnodes()
65c
66c     allocate necessary local temporary arrays on the stack
67c
68c     l_dens ... buffer for density
69c     l_scr ... workspace for guess routines
70c
71c     k_* are the offsets corrsponding to the l_* handles
72c
73c
74c     guess_mem call ...
75c
76      call guess_mem(mem1, max1e, nbf)
77c
78      status = MA_push_get(MT_DBL, max1e, 'guess_dens: dens',
79     +                     l_dens, k_dens)
80      status = MA_push_get(MT_DBL, mem1,  'guess_dens: scr',
81     +                     l_scr, k_scr)
82c
83c
84      call denat(dbl_mb(k_dens), nbf, natoms, iprint, iproc,
85     +           dbl_mb(k_scr),  mem1)
86      call square(dbl_mb(k_scr), dbl_mb(k_dens), nbf, nbf)
87c
88      call ga_put(g_dens,1,nbf,1,nbf,dbl_mb(k_scr),nbf)
89      call ga_put(g_old_dens,1,nbf,1,nbf,dbl_mb(k_scr),nbf)
90c
91c     chop stack at first item allocated
92c
93      status = MA_pop_stack(l_scr)
94      status = MA_pop_stack(l_dens)
95c
96      end
97      subroutine guess_mem(mem1, max1e, nbf)
98c
99      implicit none
100c
101      integer mem1, nbf, max1e
102c
103c..   calculate memory requirements for atomic guess routines
104c..
105      integer nb, no, ntr, nsq
106      integer i10, ipcap, iqcap, ifc, ifo, is, iu, it
107      integer ih, idc, idos, idt, idold, iss, ic, icopn, ismin
108      integer iqmin, itransf, icc
109c..
110c..    core partitioning
111c..
112      max1e = nbf*(nbf+1)/2
113c
114c     allow for maximum of 100 bfns on any given atom
115c
116      nb = 100
117      no = 50
118      ntr = nb*(nb+1)/2
119      nsq = nb * nb
120c
121c
122      i10 = 1
123      ipcap = i10 + max1e
124      iqcap = ipcap + ntr
125      ifc = iqcap + ntr
126      ifo = ifc + ntr
127      is = ifo + ntr
128      iu = is + ntr
129      it = iu + ntr
130      ih = it + ntr
131      idc = ih + ntr
132      idos = idc + ntr
133      idt = idos + ntr
134      idold = idt + ntr
135      iss = idold + ntr
136      ic = iss + ntr
137c
138      icopn = ic + nsq
139      ismin = icopn + nsq
140      iqmin = ismin + nb * no
141      itransf = iqmin + nb * no
142      icc = itransf + nsq
143      mem1 = icc + nsq - 1
144c
145c     NOTE: later inserts required for pseudopotentials
146c
147      end
148      subroutine denat(dens,nbf,nat,iprint,iproc,q,memq)
149c
150      implicit none
151c
152      integer nat, nbf
153      integer iprint, memq, iproc
154c
155c..   get starting vectors from superposition of atomic densities
156c..
157c..   routines involved are (all present in this order) :
158c..   denat  : general control routine
159c..   datoms : translate orbital info, call atomscf , gather d-matrix
160c..   atomd,tracd,trafsd,cmergd,oeigd,teigd,densid,denmad, ..
161c..   .. hamild,outpud,shalfd,atcond,starcd,creded,pdfded, ..
162c..   .. jacod,orderd,tramad
163c..   oeigd has been changed to pick up the 1-elec ints
164c..   atcond has been changed to allow for effective nuclear charge
165c..   creded/pdfded have been adapted to yield directly a d-matrix
166crz   xpsnld: implemented for use of non local pseudo potentials
167c..
168c...    start all types after one cycle closed shell in routine denscf
169c
170      real *8 q(memq),dens(*)
171c
172      integer iread, iwrite
173      common/iofile/ iread,iwrite
174c
175      logical oprinv
176      integer l2
177      integer nb, no, ntr, nsq
178      integer i10, ipcap, iqcap, ifc, ifo, is, iu, it
179      integer ih, idc, idos, idt, idold, iss, ic, icopn, ismin
180      integer iqmin, itransf, icc, last
181c
182      if (iproc.eq.0) write (iwrite,6010)
183c..
184c..    core partitioning
185c..
186      oprinv = iprint.eq.2
187      l2 = nbf*(nbf+1)/2
188c
189c..   zero total density matrix
190c..
191C      call vclr(dens,1,l2)
192      call dfill(l2,0.d0,dens,1)
193
194c
195c     allow for maximum of 100 bfns on any given atom
196c
197      nb = 100
198      no = 50
199      ntr = nb*(nb+1)/2
200      nsq = nb * nb
201c
202c
203      i10 = 1
204      ipcap = i10 + l2
205      iqcap = ipcap + ntr
206      ifc = iqcap + ntr
207      ifo = ifc + ntr
208      is = ifo + ntr
209      iu = is + ntr
210      it = iu + ntr
211      ih = it + ntr
212      idc = ih + ntr
213      idos = idc + ntr
214      idt = idos + ntr
215      idold = idt + ntr
216      iss = idold + ntr
217      ic = iss + ntr
218c
219      icopn = ic + nsq
220      ismin = icopn + nsq
221      iqmin = ismin + nb * no
222      itransf = iqmin + nb * no
223      icc = itransf + nsq
224      last = icc + nsq
225c
226c     NOTE: later inserts required for pseudopotentials
227c
228c..     get core
229c
230c..    to supply to the atomic scf routines to cover
231c..    pseudo-potential calculations
232c..     d in  dens    /  workspace in i10 (exact fit)
233c      pseudo corrections in i15 .. NOW REMOVED
234c..    **** not symmetry adapted ****
235c..
236c..
237c..   now loop over the atoms / do atomic scf and gather d-matrix
238c..
239      call datoms(q(i10),dens,oprinv,iproc,
240     +    q(ipcap), q(iqcap), q(ifc), q(ifo), q(is), q(iu),
241     +    q(it), q(ih), q(idc), q(idos), q(idt), q(idold), q(iss) ,
242     +    q(ic), q(icopn), q(ismin), q(iqmin), q(itransf), q(icc) ,
243     +    nb)
244c
245c..   print if requested
246c
247      if (oprinv.and.iproc.eq.0) then
248         write (iwrite,6020)
249         call writel(dens,nbf,.false.)
250      end if
251c
252      return
253 6010 format (/,' initial guess orbitals generated by ',
254     +        'superposition of atomic densities',/)
255 6020 format (//30x,28('-'),/,30x,'   initial guess density   ',/,30x,
256     +        28('-')//)
257      end
258      subroutine datoms(hatom,d,oprin,iproc,
259     + pcap, qcap, fc, fo, s, u, t, h, dc, dos, dt, dold, ss,
260     + cvec, copn, smin, qmin, transf, cc, nbb)
261c
262      implicit none
263      logical osatod
264c
265c...   subroutine to coordinate atom-scf calls and d-matrix gathering
266c...   for atomic startup
267c...   h,t  : full h-matrix and t-matrix to supply to atom
268c...   d    : full density matrix as return parameter
269c...
270c...   **note** data is transferred to atom directly via common/cguess/
271c
272c
273      integer nbb, iproc
274      real *8 hatom(*), d(*)
275      real *8 pcap(*), qcap(*), fc(*), fo(*), s(*), u(*), t(*)
276      real *8 h(*), dc(*), dos(*), dt(*), dold(*), ss(*)
277      real *8 cvec(*), copn(*), smin(nbb,*), qmin(nbb,*),transf(*),cc(*)
278      logical oprin
279c
280c...   commons where information comes from :
281c
282#include 'int.h'
283c
284      integer iread, iwrite
285      common /iofile/ iread,iwrite
286      integer kmin, kmax, nuct
287      common /iguess/ kmin(mxshell), kmax(mxshell), nuct(mxnat)
288c
289c...   common where info goes to :
290c
291      integer nb, no
292      parameter (nb=100, no=50)
293      integer nsym, nbas, ncsh, nosh, nccup, nsht, nitscf
294      integer nqn, n1, nconv, nbc, nbct, nstrt, ifcont
295      real *8 zn, zeta, eps, cin, vir, energ, ajmn, damp, cont
296      common /cguess/nsym,nbas(5),ncsh(5),nosh(5),nccup(6),nsht,nitscf,
297     + zn,n1(6),nqn(nb),zeta(nb),eps(no),cin,vir,energ,ajmn(24),damp,
298     + nconv,nbc(5),cont(nb),nbct,nstrt(nb),ifcont
299c
300      integer ic(6,nb),iiloc(nb,6),iisch(nb,6)
301      logical odone(mxnat)
302c
303      integer i, ii, ispdf, iat, j, k, iorb
304      integer kk, mini, maxi
305      integer kkzc, kh, isymax, is, if
306      real *8 pi32, toteng, ee, fac, znps
307c
308c..    we need xy for d and xyz for f in gathering h-ints
309c..    so set proper offsets in ioffhp (see do 140)
310c
311      integer ioffhp(4)
312      data ioffhp/0,0,3,9/
313c
314      integer maxtyp
315      parameter (maxtyp = 6)
316      integer minf(maxtyp),maxf(maxtyp)
317      data minf  / 1, 2,  5, 11, 21, 1 /
318      data maxf  / 1, 4, 10, 20, 35, 4 /
319c
320c..
321      data pi32/5.56832799683170d0/
322c
323      do i = 1,nshell
324        ii = ktype(i)
325        kmin(i) = minf(ii)
326        kmax(i) = maxf(ii)
327      enddo
328c
329      toteng = 0.0d0
330      do i = 1 , nat
331         nuct(i) = 1
332         odone(i) = .false.
333      enddo
334c
335c...  loop over atoms like adapt does
336c
337      do iat = 1 , nat
338c
339c ... eliminate ghost/dummies/point charges
340c
341         if (nuct(iat).ne.1) then
342            odone(iat) = .true.
343            go to 110
344c..
345c..   check if we have already treated this one or if it is of same
346c..   type as the one we just did
347         else if (iat.gt.1) then
348            if (odone(iat)) then
349               odone(iat) = .true.
350               go to 110
351            else if (osatod(iat,ic,iiloc,iisch,nbb)) then
352               go to 100
353            end if
354         end if
355c
356c...  gather  shell / symmetry info
357c
358         do i = 1 , 4
359            nbc(i) = 0
360         enddo
361c
362c   nbc  # shell's / symmetry
363c   iisch  contains index of shell
364c   iiloc  contains position of starting ao of shell in "real" world
365c   translate to 1 (s)
366c
367         do ii = 1 , nshell
368            i = katom(ii)
369            if (i.eq.iat) then
370               mini = kmin(ii)
371               maxi = kmax(ii)
372               kk = ktype(ii)
373               if (kk.eq.6) kk = 2
374               do iorb = mini , maxi
375                  if (iorb.eq.1) then
376                     nbc(1) = nbc(1) + 1
377                     iisch(nbc(1),1) = ii
378                     iiloc(nbc(1),1) = kloc(ii)
379                  else if (iorb.eq.2 .or. iorb.eq.5 .or. iorb.eq.11)
380     +                     then
381c..  translate to 2 (p) 3(d) or  4(f)
382                     ispdf = kk
383                     nbc(ispdf) = nbc(ispdf) + 1
384                     iisch(nbc(ispdf),ispdf) = ii
385                     iiloc(nbc(ispdf),ispdf) = kloc(ii) + iorb - mini
386                  end if
387               enddo
388            end if
389         enddo
390c..
391c..     we gathered symmetry/shell info ; now get the real thing
392c..
393         kkzc = 0
394         kh = 0
395         isymax = 0
396         do ispdf = 1 , 4
397c..      nbas = total # primitives for this symmetry
398            nbas(ispdf) = 0
399            if (nbc(ispdf).gt.0) isymax = ispdf
400            do j = 1 , nbc(ispdf)
401               ii = iisch(j,ispdf)
402               is = kstart(ii)
403               if = is + kng(ii) - 1
404c..      ic = # number of primitives /contracted /symmetry
405               ic(ispdf,j) = kng(ii)
406               nbas(ispdf) = nbas(ispdf) + kng(ii)
407c..      gather the primitives / watch the subtle use of 2-dim cspd
408               do k = is , if
409                  kkzc = kkzc + 1
410                  zeta(kkzc) = ex(k)
411                  cont(kkzc) = cspd(k,ispdf)
412c...     get contraction coeff's as we are used to
413                  ee = 2*zeta(kkzc)
414                  fac = pi32/(ee*sqrt(ee))
415                  if (ispdf.eq.2) then
416                     fac = 0.5d0*fac/ee
417                  else if (ispdf.eq.3) then
418                     fac = 0.75d0*fac/(ee*ee)
419                  else if (ispdf.eq.4) then
420                     fac = 1.875d0*fac/(ee**3)
421                  end if
422                  cont(kkzc) = cont(kkzc)*sqrt(fac)
423               enddo
424c...     in the pseudopotential case, we would be involved in
425c...     gathering  the h integrals for the contracted ao's
426c...     so that they are added in at the right time (in atomd)
427c...     use proper offset to use pure d or f functions (ioffhp)
428c...     only the comments remain ...
429               do k = 1 , j
430                  kh = kh + 1
431                  hatom(kh) = 0.0d0
432               enddo
433c...
434            enddo
435         enddo
436c..
437c..     all prepared call  atomd
438c..     zeta,cont,nbas,nbc,nbas,ic,zn are passed via cguess
439c..     energ and the density matrix dt are received via cguess
440c..     note zn is the real nuclear charge / znps is the effective charg
441c..
442         zn = zan(iat)
443         znps = zan(iat)
444c
445         call atomd(oprin,iwrite,znps,ic,isymax,hatom,
446     +    pcap, qcap, fc, fo, s, u, t, h, dc, dos, dt, dold, ss,
447     +    cvec, copn, smin, qmin, transf, cc, nbb)
448c..
449 100     toteng = toteng + energ
450c..
451c..      now add density-matrix to the molecular d-matrix
452c..
453         call creded(d,dt,iiloc,nbb)
454c..
455         odone(iat) = .true.
456 110  continue
457      enddo
458c..
459      if (iproc.eq.0) write (iwrite,6010) toteng
460c..
461      return
462 6010 format (/1x,'***** total atomic energy ',f17.8,' *** ')
463      end
464#endif /* ADRIANS_CRAP */
465
466
467
468#ifndef ADRIANS_CRAP
469      subroutine creded(d,dt,iiloc,nbb)
470c
471      implicit none
472c
473c..   routine to merge atomic density matrix into molecular one
474c..   is called straight after atom, so all info is still in
475c..   common /cguess/ : i.e. dt,nbc
476c
477      integer nbb
478      real *8 d(*),dt(*)
479      integer iiloc(nbb,6)
480c
481      integer nb, no, nsym, nbas, ncsh, nosh, nccup, nsht, nitscf, n1
482      integer nqn, nconv, nbc, nbct, nstrt, ifcont
483      real *8 zn, zeta, eps, cin, vir, energ, ajmn, damp, cont
484      parameter (nb=100, no=50)
485      common /cguess/nsym,nbas(5),ncsh(5),nosh(5),nccup(6),nsht,nitscf,
486     + zn,n1(6),nqn(nb),zeta(nb),eps(no),cin,vir,energ,ajmn(24),damp,
487     + nconv,nbc(5),cont(nb),nbct,nstrt(nb),ifcont
488c..
489      real *8 dmult(145)
490c..
491crz   order of d's and f's in GAMESS is completely different
492crz   from MOLECULE !!
493c
494c     dmult consists of :
495c     nothing for s;
496c      p functions (1 line)
497c      d functions (next 4 lines)
498c      f functions (last lines)
499c      -0.670820393 = -0.3*sqrt(5) (used for f)
500c
501      data dmult/ 1.0d0, 3*0.0d0, 1.0d0, 3*0.0d0, 1.0d0,
502     x  1.0d0, -0.5d0, -0.5d0, 3*0.0d0,
503     x -0.5d0,  1.0d0, -0.5d0, 3*0.0d0,
504     x -0.5d0, -0.5d0,  1.0d0, 3*0.0d0,
505     x  3*0.0d0, 1.0d0, 6*0.0d0, 1.0d0, 6*0.0d0, 1.0d0,
506     x  1.0d0,   4*0.0d0, -0.670820393d0, 0.0d0,-0.670820393d0, 2*0.0d0,
507     x  0.0d0,1.0d0,0.0d0,-0.670820393d0, 4*0.0d0,-0.670820393d0, 0.0d0,
508     x2*0.0d0,1.0d0,0.0d0,-0.670820393d0, 0.0d0,-0.670820393d0, 3*0.0d0,
509     x   0.0d0, -0.670820393d0,  0.0d0,  1.2d0, 4*0.0d0, -0.3d0,  0.0d0,
510     x 2*0.0d0, -0.670820393d0,  0.0d0,  1.2d0,  0.0d0, -0.3d0, 3*0.0d0,
511     x  -0.670820393d0,         4*0.0d0, 1.2d0,  0.0d0, -0.3d0, 2*0.0d0,
512     x 2*0.0d0, -0.670820393d0,  0.0d0, -0.3d0,  0.0d0,  1.2d0, 3*0.0d0,
513     x  -0.670820393d0,     4*0.0d0,    -0.3d0,  0.0d0,  1.2d0, 2*0.0d0,
514     x   0.0d0, -0.670820393d0, 0.0d0, -0.3d0, 4*0.0d0,  1.2d0,   0.0d0,
515     x 9*0.0d0,  1.0d0/
516c.......................................................................
517c         among the variables used are:
518c
519c              nsym       - highest l-quantum no. used in atomic calc.
520c              dt         - atomic density matrix
521c              d          - area for final molecular density matrix.
522c.......................................................................
523c
524      integer k, l, m
525      integer lm, nbci, noff, kdim, kmone, kpoint
526      real *8 factor
527c..
528c..    triangle statement function
529c..
530      integer itrian, i, j
531      itrian(i,j) = max0(i,j)*(max0(i,j)-1)/2 + min0(i,j)
532
533c..
534      lm = 0
535      do k = 1 , nsym
536         nbci = nbc(k)
537         if (k.eq.1) then
538c.......................................................................
539c
540c       s-orbitals, simple distribution of matrix elements.
541c
542c.......................................................................
543            do l = 1 , nbci
544               do m = 1 , l
545                  noff = itrian(iiloc(l,1),iiloc(m,1))
546                  lm = lm + 1
547                  d(noff) = dt(lm)
548               enddo
549            enddo
550         else if (k.le.4) then
551            kdim = k*(k+1)/2
552            kmone = k - 1
553            kpoint = kmone*(kmone+1)*(((6*kmone+24)*kmone+26)*kmone+4)
554     +               /120
555            factor = 1.d0/(k+k-1)
556            call pdfded(k,kdim,d,dt,factor,dmult(kpoint),lm,nbci,
557     +                  iiloc(1,k))
558         end if
559      enddo
560c..
561      return
562      end
563#endif /* ADRIANS_CRAP */
564
565
566#ifndef ADRIANS_CRAP
567
568      subroutine pdfded(k,kdim,d,dhelp,factor,dmult,lm,nbci,iiloc)
569      implicit none
570c.......................................................................
571c
572c     routine for distributing atomic densities to the molecular
573c     density matrix. note that the number of orbitals differ in
574c     the molecular and atomic case for d and f orbitals. transformat-
575c     ion matrices are provided in dmult(*,*). to work out these tables
576c     real atomic orbitals are needed. for f functions these are:
577c              sqrt(1/60)*(2zzz-3zyy-3zxx)
578c              sqrt(1/40)*(4zzy-yyy-xxy)
579c              sqrt(1/40)*(4zzx-xxx-xyy)
580c                    xyz
581c              sqrt(1/4)*(xxz-yyz)
582c              sqrt(1/24)*(3xxy-yyy)
583c              sqrt(1/24)*(3xyy-xxx)
584c     normalization of primitives is given by (xyz:xyz)=1, (xxy:xxy)=3
585c     (xxx:xxx)=15.
586c
587c.......................................................................
588c..
589      integer k, kdim, lm, nbci
590      real *8 d(*),dhelp(*),dmult(kdim,kdim)
591      real *8 factor
592      integer iiloc(nbci)
593c
594      integer l, na, m, nb
595      integer lmsave, nbrang, noff
596      real *8 delem
597c
598      integer i, j, itria2
599      itria2(i,j,na) = (max0(i,j)+na-1)*(max0(i,j)+na-2)/2 + min0(i,j)
600c..
601c..
602      do l = 1 , nbci
603         lmsave = lm
604         do na = 1 , kdim
605            lm = lmsave
606            do m = 1 , l
607               noff = itria2(iiloc(m),iiloc(l),na) - 1
608               lm = lm + 1
609               delem = dhelp(lm)*factor
610               nbrang = kdim
611               if (m.eq.l) nbrang = na
612               do nb = 1 , nbrang
613                  noff = noff + 1
614                  d(noff) = delem*dmult(na,nb)
615               enddo
616            enddo
617         enddo
618      enddo
619c..
620      return
621      end
622      subroutine guess_dens(geom, basis, g_dens)
623      implicit none
624#include "mafdecls.fh"
625      integer geom, basis       ! [input] handles
626      integer g_dens            ! [input] GA returns superposition of AO dens
627c
628      integer memscr            ! Size of scratch required in doubles
629      integer l_scr,k_scr
630      logical status
631c
632c     Return in g_dens a full square matrix with atomic HF densities
633c     along the diagonal and zeroes elsewhere
634c
635      call ga_zero(g_dens)
636c
637      call guess_mem(memscr)
638
639#ifdef DEBUG
640      if (ga_nodeid().eq.0) then
641         print*,'Maximum primitives:',nprim
642         print*,'Scratch space:',memscr
643         call flush_output()
644      endif
645#endif
646
647      status = ma_push_get(MT_DBL,memscr,'guess scratch',l_scr,k_scr)
648      call denat(geom,basis,1,dbl_mb(k_scr),memscr,g_dens)
649      status = ma_pop_stack(l_scr)
650
651      return
652      end
653