1C 2C zfock_cs_build.F 3C 4C Builds closed shell complex Fock matrix from a complex density 5C matrix. 6C 7C A comment on notation. The "z" in g_zfock means this GA is 8C complex data type. The corresponding real and imaginary parts are 9C called g_fockre and g_fockim, which are both real data types, so 10C no "z". 11C 12C All quantities are in the atomic orbital (AO) basis. 13C 14 15 subroutine zfock_cs_build (params, g_zdens, energies, g_zfock) 16 implicit none 17 18#include "errquit.fh" 19#include "mafdecls.fh" 20#include "stdio.fh" 21#include "global.fh" 22#include "msgids.fh" 23#include "util.fh" 24#include "cdft.fh" 25#include "matutils.fh" 26#include "zora.fh" 27#include "case.fh" 28#include "rt_tddft.fh" 29 30 31C == Inputs == 32 type(rt_params_t), intent(in) :: params 33 integer, intent(in) :: g_zdens 34 35 36C == Outputs == 37 type(rt_energies_t), intent(out) :: energies 38 integer, intent(in) :: g_zfock 39 40 41C == Parameters == 42 character(*), parameter :: pname = "zfock_cs_build: " 43 44 45C == Variables == 46 logical need_exactexch 47 logical need_dftxc 48 integer g_densre !real part of dens mat 49 integer g_densim !imag part of dens mat 50 integer g_fockre !real part of fock mat 51 integer g_fockim !imag part of fock mat 52 integer g_v1, g_v2 !potentials--results of each piece of F 53 double precision ener1, ener2 !energies-results of each piece of F 54 double precision dftExc(2) 55 logical status 56 double precision elapsed, elapsed_tot 57 integer g_zfock_mo, g_zcapmo ! for MO damping 58 59C == External routines == 60 logical xc_gotxc 61 external xc_gotxc 62 external sandbox_complex 63 64 65 call rt_tddft_cs_confirm (params,'zfock_cs_build.F') 66 67 68C == Initializations == 69 if (params%prof) call prof_start (elapsed_tot) 70 71 72 if (params%nExc .ne. 1) 73 $ call errquit (pname//"might not work for nExc /= 1",0,0) 74 75 if (oskel) 76 $ call errquit (pname//"not working with symmetry", 0, 0) 77 78 if (nbf_ao .ne. params%nbf_ao) call errquit (pname// 79 $ "nbf_ao in params /= nbf_ao in cdft header", 0, 0) 80 81 82C 83C If asked for (usually a "static" calculation), dont build but 84C instead use stored Fock matrix. 85C 86 if (params%skip_fock) then 87 call ga_copy (params%g_zfock_ao_init(1), g_zfock) 88 call rt_tddft_print_notice ("Static CS Fock matrix") 89 return 90 endif 91 92 93 call ga_zero (g_zfock) 94 need_exactexch = abs(xfac(1)).gt.1d-8 95 need_dftxc = xc_gotxc() 96 97 98 if (.not. ga_create(mt_dbl, nbf_ao, nbf_ao, 99 $ "Re[Dens]", 0, 0, g_densre)) 100 $ call errquit ("couldnt create g_densre", 0, GA_ERR) 101 if (.not. ga_duplicate(g_densre, g_densim, "Im[Dens]")) 102 $ call errquit ("couldnt duplicate g_densre", 0, GA_ERR) 103 if (.not. ga_duplicate(g_densre, g_fockre, "Re[F]")) 104 $ call errquit ("couldnt duplicate g_densre", 0, GA_ERR) 105 if (.not. ga_duplicate(g_densre, g_fockim, "Im[F]")) 106 $ call errquit ("couldnt duplicate g_densre", 0, GA_ERR) 107 if (.not. ga_duplicate(g_densre, g_v1, "V1")) 108 $ call errquit ("couldnt duplicate g_densre", 0, GA_ERR) 109 if (.not. ga_duplicate(g_densre, g_v2, "V2")) 110 $ call errquit ("couldnt duplicate g_densre", 0, GA_ERR) 111 112 113C == Extract real and imag parts of density matrix == 114 if (params%prof) call prof_start (elapsed) 115 call ga_zero (g_densre) 116 call ga_zero (g_densim) 117 call convert_z2d (g_zdens, g_densre, g_densim) 118 if (params%prof) call prof_end (elapsed, "Fock CS z2d") 119 120 121C 122C == Compute complex Fock matrix == 123C 124C For each piece we compute the energy and potential, then 125C accumulate the result in the real or imag part of F. Note that 126C the only piece that involves the imag part of the density matrix 127C is the exact exchange. We also only call the DFT XC routine if we 128C need it, i.e., anything but pure Hartree-Fock. 129C 130 energies%core = 0d0 131 energies%coul = 0d0 132 energies%xc(1) = 0d0 133 energies%xc(2) = 0d0 134 135 136 call ga_zero (g_fockre) 137 call ga_zero (g_fockim) 138 139 140C == Standard core (kinetic+potential) == 141 call zfock_cs_core (params, g_densre, g_v1) 142 call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre) 143C energies%core = ga_ddot (g_densre, g_v1) 144 145 146C == Scalar ZORA core == 147 if (params%lzorasf) then 148C call zfock_cs_core_scalarzora (params, g_densre, g_v1) 149 call ga_add (1d0, params%g_zora_sf(1), 1d0, g_fockre, g_fockre) 150C energies%core = energies%core + ga_ddot (g_densre, g_v1) 151 endif 152 153 energies%core = ga_ddot (g_densre, g_fockre) 154 155 156C == Standard coul + standard exch == 157 if ((need_exactexch).and.(.not.cam_exch).and.(.not.cdfit)) then 158 call zfock_cs_coul_exchre(params, g_densre, ener1, ener2, g_v1) 159 call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre) 160 energies%coul = ener1 161 energies%xc(1) = ener2 162 energies%xc(2) = 0d0 163 164 call zfock_cs_exchim (params, g_densim, ener1, g_v1) 165 call ga_add (1d0, g_v1, 1d0, g_fockim, g_fockim) 166 energies%xc(1) = energies%xc(1) + ener1 167 168 169C == Standard coul + CAM exchange == 170 elseif ((need_exactexch).and.(cam_exch).and.(.not.cdfit)) then 171 call zfock_cs_coul (params, g_densre, ener1, g_v1) 172 call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre) 173 energies%coul = ener1 174 175 call zfock_cs_exch (params, g_densre, g_densim, .true., 176 $ ener1, g_v1, g_v2) 177 call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre) 178 call ga_add (1d0, g_v2, 1d0, g_fockim, g_fockim) 179 energies%xc(1) = ener1 180 energies%xc(2) = 0d0 181 182 183C == CD fit coul + standard exch == 184 elseif ((need_exactexch).and.(.not. cam_exch).and.(cdfit)) then 185 call zfock_cs_coulcdfit (params, g_densre, ener1, g_v1) 186 call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre) 187 energies%coul = ener1 188 189 call zfock_cs_exch (params, g_densre, g_densim, .false., 190 $ ener1, g_v1, g_v2) 191 call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre) 192 call ga_add (1d0, g_v2, 1d0, g_fockim, g_fockim) 193 energies%xc(1) = ener1 194 195 196C == CD fit coul + CAM exch == 197 elseif ((need_exactexch).and.(cam_exch).and.(cdfit)) then 198 call zfock_cs_coulcdfit (params, g_densre, ener1, g_v1) 199 call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre) 200 energies%coul = ener1 201 202 call zfock_cs_exch (params, g_densre, g_densim, .true., 203 $ ener1, g_v1, g_v2) 204 call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre) 205 call ga_add (1d0, g_v2, 1d0, g_fockim, g_fockim) 206 energies%xc(1) = ener1 207 energies%xc(2) = 0d0 208 209 210C == Standard coul + no exchange == 211 elseif ((.not.need_exactexch).and.(.not.cdfit)) then 212 energies%xc(1) = 0d0 213 energies%xc(2) = 0d0 214 215 call zfock_cs_coul (params, g_densre, ener1, g_v1) 216 call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre) 217 energies%coul = ener1 218 219 220C == CD fit coul + no exchange == 221 elseif ((.not.need_exactexch).and.(cdfit)) then 222 energies%xc(1) = 0d0 223 energies%xc(2) = 0d0 224 225 call zfock_cs_coulcdfit (params, g_densre, ener1, g_v1) 226 call ga_add (1d0, g_v1, 1d0, g_fockre, g_fockre) 227 energies%coul = ener1 228 229 else 230 call errquit (pname//"Unsupported combination of flags", 0, 0) 231 endif 232 233 234C == DFT XC == 235 if (need_dftxc) then 236 dftExc(1) = 0d0 237 dftExc(2) = 0d0 238 call zfock_cs_dftxc (params, g_densre, dftExc, g_v1) 239 call ga_dadd (1d0, g_v1, 1d0, g_fockre, g_fockre) 240 energies%xc(1) = energies%xc(1) + dftExc(1) 241 energies%xc(2) = energies%xc(2) + dftExc(2) 242 endif 243 244 245C 246C Fudge symmetry 247C 248 if (oskel) call errquit (pname//"not working with oskel", 0, 0) 249 250 251 252C == Pack real and imag parts of Fock matrix back into g_zfock == 253C No need to repack density matrix, as its value should 254C not have changed. Perhaps double check this?? 255 256 if (params%prof) call prof_start (elapsed) 257 call convert_d2z (1d0, g_fockre, 1d0, g_fockim, g_zfock) 258 if (params%prof) call prof_end (elapsed, "Fock CS d2z") 259 260 261 262C 263C Apply molecular orbital complex absorbing potential (MO CAP) 264C 265 if (params%mocap_active) then 266 if (.not. ga_create(mt_dcpl, params%ns_mo, params%ns_mo, 267 $ "Fock MO tmp" , 0, 0, g_zfock_mo)) 268 $ call errquit(pname//"GA alloc failed",0, GA_ERR) 269 270 if (.not. ga_create(mt_dcpl, params%ns_mo, params%ns_mo, 271 $ "Fock MO tmp" , 0, 0, g_zcapmo)) 272 $ call errquit(pname//"GA alloc failed",0, GA_ERR) 273 274 call canorg_trans (params,"F","AO->MO", g_zfock, g_zfock_mo) 275 call rt_tddft_mocap (params, g_zfock_mo, g_zcapmo) 276 call ga_zadd (z1, g_zfock_mo, z1, g_zcapmo, g_zfock_mo) 277 call canorg_trans (params,"F","MO->AO", g_zfock_mo, g_zfock) 278 call rt_tddft_print_notice ("Applied closed shell MO CAP") 279 280 if (.not. ga_destroy(g_zfock_mo)) 281 $ call errquit(pname//"GA destroy failed",0, GA_ERR) 282 283 if (.not. ga_destroy(g_zcapmo)) 284 $ call errquit(pname//"GA destroy failed",0, GA_ERR) 285 endif 286 287 288C == Spatial complex absorbing boundary conditions (CAP) == 289 if (params%cap_active) then 290 call ga_zadd (z1, g_zfock, z1, params%g_zcap, g_zfock) 291 call rt_tddft_print_notice ("Applied spatial CAP to "// 292 $ "restricted Fock matrix") 293 endif 294 295 296C 297C Extra checks, disabled for speed. 298C 299 if (params%checklvl .ge. 2) then 300 if (.not. mat_is_hermitian (g_zfock, params%tol_zero)) 301 $ call errquit ("F not hermitian in fock builder", 0, 0) 302 if (.not. mat_is_symmetric (g_fockim, "A", params%tol_zero)) 303 $ call errquit ("Im[F] not antisym in fock builder", 0, 0) 304 if (energies%xc(2) > 1d-8) 305 $ call errquit (pname//"Exc(2) /= 0?", 0, 0) 306 endif 307 308 status = .true. 309 status=status.and.ga_destroy(g_densre) 310 status=status.and.ga_destroy(g_densim) 311 status=status.and.ga_destroy(g_fockre) 312 status=status.and.ga_destroy(g_fockim) 313 status=status.and.ga_destroy(g_v1) 314 status=status.and.ga_destroy(g_v2) 315 316 if (.not.status) 317 $ call errquit (pname//"couldnt free arrays", 0, GA_ERR) 318 319 if (params%prof) call prof_end (elapsed_tot, 320 $ "Fock CS total build") 321 322 end subroutine 323 324c $Id$ 325