1c
2c        construct the transition density
3         subroutine get_transden(spin,iroot,ipol,filename,basis,
4     &        g_movecs, g_tdens)
5c
6         implicit none
7c
8#include "errquit.fh"
9#include "global.fh"
10#include "tcgmsg.fh"
11#include "msgtypesf.h"
12#include "mafdecls.fh"
13#include "msgids.fh"
14#include "inp.fh"
15#include "util.fh"
16#include "stdio.fh"
17#include "bas.fh"
18#include "geom.fh"
19c
20         character*8 spin
21         integer iroot
22         integer ipol
23         character*256 filename
24         integer basis            ! AO basis set handle
25         integer g_movecs(ipol)   ! MO vectors
26         integer g_tdens(ipol)    ! Transition density matrix
27c
28         integer l_trials, k_trials
29c
30         integer i,m,n
31         logical tda
32         integer nroots
33         integer nocc(2)
34         integer nmo(2)
35         integer nfc(2)
36         integer nfv(2)
37         integer nov(2)
38         integer icntr,itmom
39         double precision r,cntr(3),tmom(20)
40         integer nbf_ao
41c
42         logical nodezero
43c
44c        CI vectors are written out as X+Y and X-Y
45         integer g_x(2),g_y(2)
46         integer g_temp(2) ! scratch space
47c
48         integer inntsize,ddblsize,logisize,ok
49c
50         character*32 pname
51         pname = 'get_transden: '
52c
53c        preliminaries
54         inntsize=MA_sizeof(MT_INT,1,MT_BYTE)
55         ddblsize=MA_sizeof(MT_DBL,1,MT_BYTE)
56         logisize=MA_sizeof(MT_LOG,1,MT_BYTE)
57         ok = 0
58         call ga_sync()
59         nodezero=(ga_nodeid().eq.0)
60         if (.not.bas_numbf(basis,nbf_ao))
61     &      call errquit(pname//'bas_numbf failed',0,0)
62c
63c        initialization
64         do i=1,ipol
65	   call ga_zero(g_tdens(i))
66         end do ! ipol
67         do icntr=1,3
68           cntr(icntr)=0.0d0
69         enddo
70         do itmom=1,20
71           tmom(itmom)=0.0d0
72         enddo
73c
74c        Read header information from civecs file
75         if (nodezero) then
76           open(unit=69,file=filename,form='unformatted',
77     &          status='unknown',err=1000)
78            rewind(69)
79            read(69,err=1001) tda
80            read(69,err=1001) ipol
81            read(69,err=1001) nroots
82            read(69,err=1001) nocc(1),nocc(2)
83            read(69,err=1001) nmo(1),nmo(2)
84            read(69,err=1001) nfc(1),nfc(2)
85            read(69,err=1001) nfv(1),nfv(2)
86            read(69,err=1001) nov(1),nov(2)
87            read(69,err=1001)
88c
89            if (ipol.eq.1) nocc(2)=0
90            if (ipol.eq.1) nmo(2)=0
91            if (ipol.eq.1) nfc(2)=0
92            if (ipol.eq.1) nfv(2)=0
93c
94            do i=1,ipol
95             nov(i)=(nmo(i)-nfv(i)-nocc(i))*(nocc(i)-nfc(i))
96            end do ! ipol
97            if (ipol.eq.1) nov(2)=0
98           close(unit=69,status='keep',err=1002) ! file
99           ok = 1
100         end if ! nodezero
101c
102c        broadcast status and variables to other nodes
103         call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0)
104         call ga_brdcst(Msg_Vec_NMO+MSGINT, tda, logisize, 0)
105         call ga_brdcst(Msg_Vec_NMO+MSGINT,ipol, inntsize, 0)
106         call ga_brdcst(Msg_Vec_NMO+MSGINT,nroots, inntsize, 0)
107         call ga_brdcst(Msg_Vec_NMO+MSGINT, nocc, inntsize*2, 0)
108         call ga_brdcst(Msg_Vec_NMO+MSGINT, nmo, inntsize*2, 0)
109         call ga_brdcst(Msg_Vec_NMO+MSGINT, nfc, inntsize*2, 0)
110         call ga_brdcst(Msg_Vec_NMO+MSGINT, nfv, inntsize*2, 0)
111         call ga_brdcst(Msg_Vec_NMO+MSGINT, nov, inntsize*2, 0)
112         call ga_sync()
113c
114c        Set up X, Y vectors and transition density matrices
115         do i=1,ipol
116            if (.not.ga_create(mt_dbl,nov(i),nroots,'X vector',
117     1        -1,-1,g_x(i))) call errquit
118     2       (pname//'failed to create g_x',0, GA_ERR)
119             call ga_zero(g_x(i))
120             if (.not.ga_create(mt_dbl,nbf_ao,nbf_ao,'temp',
121     1        -1,-1,g_temp(i))) call errquit
122     2       (pname//'failed to create g_temp',0, GA_ERR)
123             call ga_zero(g_temp(i))
124c
125             if (.not.tda) then
126              if (.not.ga_create(mt_dbl,nov(i),nroots,'Y vector',
127     1          -1,-1,g_y(i))) call errquit
128     2          (pname//'failed to create g_y',0, GA_ERR)
129              call ga_zero(g_y(i))
130             end if ! .not. tda
131         end do  ! ipol
132c
133c        Read remainder of the civecs file
134         if (nodezero) then
135            open(unit=69,file=filename,form='unformatted',
136     &                status='unknown',err=1000)
137            rewind(69)
138            read(69,err=1001) tda
139            read(69,err=1001) ipol
140            read(69,err=1001) nroots
141            read(69,err=1001) nocc(1),nocc(2)
142            read(69,err=1001) nmo(1),nmo(2)
143            read(69,err=1001) nfc(1),nfc(2)
144            read(69,err=1001) nfv(1),nfv(2)
145            read(69,err=1001) nov(1),nov(2)
146            read(69,err=1001)
147         end if ! nodezero
148c
149         do n = 1,nroots
150c
151            if (nodezero) then
152               read(69) r  ! energy of root
153               read(69) r  ! s2
154            end if
155c
156            do i=1,ipol
157c
158c            Allocate memory
159             if (.not.ma_push_get(mt_dbl,nov(i),"slice",l_trials,
160     &          k_trials))
161     &          call errquit(trim(pname)//"failed to alloc slice",0,0)
162c
163             if (.not.tda) then
164               if (nodezero) then
165                call sread(69,dbl_mb(k_trials),nov(i))
166                call ga_put(g_x(i),1,nov(i),n,n,dbl_mb(k_trials),nov(i))
167                call sread(69,dbl_mb(k_trials),nov(i))
168                call ga_put(g_y(i),1,nov(i),n,n,dbl_mb(k_trials),nov(i))
169               end if ! nodezero
170             else
171               if (nodezero) then
172                call sread(69,dbl_mb(k_trials),nov(i))
173                call ga_put(g_x(i),1,nov(i),n,n,dbl_mb(k_trials),nov(i))
174               end if ! nodezero
175             end if  !tda
176c
177c            Deallocate memory
178             if (.not.ma_pop_stack(l_trials))
179     &          call errquit(trim(pname)//"failed to pop stack",0,0)
180c
181            end do ! ipol
182         end do ! nroots
183c
184         if (nodezero) close(unit=69,status='keep',err=1002) ! file
185c
186         call ga_brdcst(Msg_Vec_Stat+MSGINT, ok, inntsize, 0)
187         call ga_sync()
188c
189c        Separate out X and Y vectors if the calculations are not TDA
190c        Nothing needs to be done for TDA as Y is zero
191         if (.not.tda) then
192           do i=1,ipol
193              call ga_add(0.5d0,g_x(i), 0.5d0,g_y(i),g_x(i))
194              call ga_add(1.0d0,g_x(i),-1.0d0,g_y(i),g_y(i))
195           enddo
196         end if
197c
198c        calculate X component of the transition density matrix
199         call tddft_transfm(iroot,g_x,g_movecs,nbf_ao,nocc,nmo,
200     &             nfc,nfv,ipol,g_temp,1)  ! x-transition density
201         do i = 1,ipol
202          call multipole_density(basis,cntr,3,g_temp(i),tmom,20)  ! transition moments
203          call ga_copy(g_temp(i),g_tdens(i))
204         end do
205c
206c        calculate Y component of the transition density matrix
207         if (.not.tda) then
208           do i = 1,ipol
209                call ga_zero(g_temp(i))
210           end do
211           call tddft_transfm(iroot,g_y,g_movecs,nbf_ao,nocc,nmo,
212     &             nfc,nfv,ipol,g_temp,1)  ! y-transition density
213c
214c          accumulate the Y component of the transition density matrix
215           do i = 1,ipol
216              call multipole_density(basis,cntr,3,g_temp(i),tmom,20)  ! transition moments
217              call ga_add(1.d0,g_tdens(i),1.d0,g_temp(i),g_tdens(i))
218           end do
219         end if  ! tda
220c
221         if (ipol.eq.1) then
222          do i=1,20
223            tmom(i)=tmom(i)*dsqrt(2.0d0)
224          enddo
225         end if
226c
227         if (ga_nodeid().eq.0) then
228                write(luout,*) " *** tmom(2)***: ", tmom(2)
229                write(luout,*) " *** tmom(3)***: ", tmom(3)
230                write(luout,*) " *** tmom(4)***: ", tmom(4)
231         end if
232c
233c        symmetrize the transition density matrix
234         do i = 1,ipol
235             call ga_symmetrize(g_tdens(i))
236         enddo
237c
238c        calculate total, spin density or individual components
239         If (ipol.eq.2) Then
240            If (Spin.eq.'TOTAL') Then
241              Call GA_dAdd(1.d0,g_tdens(1),1.d0,g_tdens(2),g_tdens(1))
242            Else If (Spin.eq.'SPINDENS') Then
243              Call GA_dAdd(1.d0,g_tdens(1),-1.d0,g_tdens(2),g_tdens(1))
244            Else If (Spin.eq.'ALPHA') Then
245            Else If (Spin.eq.'BETA') Then
246               Call GA_Copy(g_tdens(2),g_tdens(1))
247            End If
248         End If ! ipol check
249c
250c        cleanup
251         do i=1,ipol
252            if (.not.ga_destroy(g_x(i))) call errquit
253     2       (pname//'failed to destroy g_x',0, GA_ERR)
254            if (.not.ga_destroy(g_temp(i))) call errquit
255     2       (pname//'failed to destroy g_temp',0, GA_ERR)
256            if (.not.tda) then
257              if (.not.ga_destroy(g_y(i))) call errquit
258     2          (pname//'failed to destroy g_y',0, GA_ERR)
259            end if ! tda
260         enddo ! ipol
261c
262      return
263c
264 1000 call errquit(pname//'failed to open file',0,-1)
265 1001 call errquit(pname//'failed to read file',0,-1)
266 1002 call errquit(pname//'failed to close file',0,-1)
267c
268      end
269c $Id$
270