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