1c
2c $Id$
3c
4      subroutine get_density(g_dens,g_vecs,theory,scftype,nalpha,
5     *                       nbeta,nbfs,rtdb)
6      implicit none
7#include "errquit.fh"
8c
9c This routine gets alpha and beta densities and passes them back in g_dens.
10c
11#include "mafdecls.fh"
12#include "rtdb.fh"
13      integer g_dens(*)    ! [output] density handles
14      integer g_vecs(*)    ! [output] vector handles
15      character*3 theory   ! [input]  HF or DFT
16      character*4 scftype  ! [input]  RHF, UHF, or ROHF
17      integer nalpha       ! [input]  number of alpha occupieds
18      integer nbeta        ! [input]  number of beta occupieds
19      integer nbfs         ! [input]  number of basis functions
20      integer rtdb         ! [input]  RTDB handle
21c
22      integer l_evals, k_evals, l_occ, k_occ
23      double precision rhffact
24      character*255 movecs
25c
26      logical movecs_read
27      external movecs_read
28c
29c Temporarily disable ROHF
30c
31      if (scftype.eq.'ROHF')
32     *  call errquit('get_density: ROHF is not supported yet',555,
33     &       CAPMIS_ERR)
34c
35      if (theory.eq.'HF') then
36        if (.not. rtdb_cget(rtdb, 'scf:input vectors', 1, movecs))
37     *    call errquit('get_density: SCF MO vectors not defined',555,
38     &       RTDB_ERR)
39      else if (theory.eq.'DFT') then
40        if (.not. rtdb_cget(rtdb, 'dft:input vectors', 1, movecs))
41     *   call errquit('get_density: DFT MO vectors not defined',555,
42     &       RTDB_ERR)
43      else
44        call errquit('get_density:theory is unrecognized',555,
45     &       INPUT_ERR)
46      endif
47c
48c get additional space to store the eigenvalues and occupation numbers
49c
50      if (.not. ma_push_get(mt_dbl, nbfs,'MO evals', l_evals, k_evals))
51     *  call errquit('get_density: could not allocate l_evals',nbfs,
52     &       MA_ERR)
53      if (.not. ma_push_get(mt_dbl, nbfs,'occ. numbers', l_occ, k_occ))
54     *  call errquit('get_density: could not allocate l_occ',nbfs,
55     &       MA_ERR)
56c
57c read vectors
58c
59      if (.not. movecs_read (movecs, 1, dbl_mb(k_occ),
60     *  dbl_mb(k_evals), g_vecs(1)))
61     *  call errquit('get_density: could not read mo vectors', 555,
62     &       DISK_ERR)
63      if (scftype.eq.'UHF') then
64        if (.not. movecs_read (movecs, 2, dbl_mb(k_occ),
65     *    dbl_mb(k_evals), g_vecs(2)))
66     *    call errquit('get_density: could not read beta vectors', 555,
67     &       DISK_ERR)
68      endif
69c
70c free additional space
71c
72      if (.not.ma_pop_stack(l_occ))
73     *  call errquit('get_density:ma free occ',555, MA_ERR)
74      if (.not.ma_pop_stack(l_evals))
75     *  call errquit('get_density:ma free eval',555, MA_ERR)
76c
77c form the densities
78c
79      rhffact = 2.0
80      if (scftype.eq.'UHF') rhffact = 1.0
81      call ga_dgemm('n', 't', nbfs, nbfs, nalpha, rhffact,
82     *   g_vecs(1), g_vecs(1), 0.0d0, g_dens(1))
83      call ga_symmetrize(g_dens(1))
84      if (scftype.eq.'UHF') then
85        call ga_dgemm('n', 't', nbfs, nbfs, nbeta, rhffact,
86     *    g_vecs(2), g_vecs(2), 0.0d0, g_dens(2))
87        call ga_symmetrize(g_dens(2))
88      endif
89c
90      end
91