1c 2c Determine the HOMO-LUMO gap 3c 4 subroutine calc_homolumogap(k_eval,nelec,rlshift,homo,lumo,gap) 5c 6 implicit none 7c 8#include "errquit.fh" 9#include "stdio.fh" 10#include "mafdecls.fh" 11c 12 integer k_eval(2) ! Pointers to the eigenvalue arrays. Only (1) is used 13 integer nelec 14 double precision rlshift 15 double precision homo 16 double precision lumo 17 double precision gap 18c 19c Determine the HOMO, LUMO and gap 20c 21 homo = Dbl_MB(k_eval(1)+nelec-1) ! Extract the contents of the eigenvalue array 22 lumo = Dbl_MB(k_eval(1)+nelec) 23 gap = min(gap, (lumo-homo-rlshift)) 24c 25 return 26 end 27c 28c Transform Fock To Orthonormal 29c 30 subroutine trans_fock_to_ortho(g_tmp2,nbf_mo,g_sm12,g_fockso) 31c 32 implicit none 33c 34#include "errquit.fh" 35#include "global.fh" 36#include "consts.fh" 37c 38 integer g_tmp2 ! temp scratch array 39 integer nbf_mo ! number of molecular orbitals 40 integer g_sm12 ! S^(-1/2) 41 integer g_fockso(2) ! 1: real, 2: imag 42c 43c Real part of the Fock matrix 44c 45 call ga_zero(g_tmp2) 46 call ga_dgemm('T', 'N', nbf_mo, nbf_mo, nbf_mo, one, 47 & g_sm12, g_fockso(1), zero, g_tmp2) 48 call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one, 49 & g_tmp2, g_sm12, zero, g_fockso(1)) 50c 51c Imag part of the Fock matrix 52c 53 call ga_zero(g_tmp2) 54 call ga_dgemm('T', 'N', nbf_mo, nbf_mo, nbf_mo, one, 55 & g_sm12, g_fockso(2), zero, g_tmp2) 56 call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one, 57 & g_tmp2, g_sm12, zero, g_fockso(2)) 58c 59 return 60 end 61c 62c Diagonalize complex Fock Matrix 63c 64 subroutine diag_fock(nbf_mo,ia,g_fockso,ibuff,g_moso,iwork,irwork, 65 & k_eval,trace,llwork,info) 66c 67 implicit none 68c 69#include "errquit.fh" 70#include "global.fh" 71#include "mafdecls.fh" 72#include "stdio.fh" 73#include "util.fh" 74c 75 integer nbf_mo 76 integer ia 77 integer g_fockso(2) 78 integer ibuff 79 integer g_moso(2) 80 integer iwork 81 integer irwork 82 integer k_eval(2) 83 double precision trace 84 integer llwork 85 integer info 86c 87 integer i,j,i1 88 double precision ddot 89 external ddot 90c 91c Prepare arrays for diagonalization 92c 93 trace = 0.d0 94 do i = 1, nbf_mo 95 do j = 1, nbf_mo 96 DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1))=dcmplx(0.0, 0.0) 97 enddo 98 enddo 99 do i = 1, nbf_mo 100 call ga_get(g_fockso(1), 1,i, i,i, dbl_mb(ibuff),1) 101 do j=1,i 102 DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1))= 103 = dcmplx(dbl_mb(ibuff+j-1),0d0) 104 enddo 105 call ga_get(g_fockso(2), 1,i, i,i, dbl_mb(ibuff),1) 106 do j=1,i 107 DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1))= 108 $ DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1)) 109 $ +dcmplx(0d0,dbl_mb(ibuff+j-1)) 110 enddo 111 enddo 112 call ga_zero(g_moso(1)) 113 call ga_zero(g_moso(2)) 114c 115c Call the diagonalizer (complex diagonalizer) 116c 117 call zheev( 'V', 'U', nbf_mo, DCpl_mb(ia), nbf_mo, 118 $ Dbl_mb(k_eval(1)), 119 $ DCpl_mb(iwork), llwork, Dbl_mb(irwork), info ) 120 do i = 1, nbf_mo 121 do j = 1, nbf_mo 122 dbl_mb(ibuff+j-1)=0.0d0 123 dbl_mb(ibuff+j-1)=dble(DCpl_mb(ia+nbf_mo*(i-1)+(j-1))) 124 enddo 125 i1=i 126 call ga_put(g_moso(1),1,nbf_mo,i1,i1,dbl_mb(ibuff),1) 127 trace = ddot(nbf_mo,dbl_mb(ibuff),1,dbl_mb(ibuff),1) 128 do j = 1, nbf_mo 129 dbl_mb(ibuff+j-1)=0.0d0 130 dbl_mb(ibuff+j-1)= 131 $ dimag(dcmplx(DCpl_mb(ia+nbf_mo*(i-1)+(j-1)))) 132 enddo 133 i1=i 134 call ga_put(g_moso(2),1,nbf_mo,i1,i1,dbl_mb(ibuff),1) 135 trace = ddot(nbf_mo,dbl_mb(ibuff),1,dbl_mb(ibuff),1) 136 enddo 137c 138 return 139 end 140c 141c Back-transform eigenvectors with S^-1/2. 142c 143 subroutine trans_vec_to_ao(nbf_mo,g_sm12,g_moso,g_fockso) 144c 145 implicit none 146c 147#include "errquit.fh" 148#include "global.fh" 149#include "stdio.fh" 150#include "consts.fh" 151c 152 integer nbf_mo 153 integer g_sm12 ! contains S^-1/2 154 integer g_moso(2) ! MO vecs 1: real, 2: imag 155 integer g_fockso(2) ! being used as scratch 156c 157c Back-transform eigenvectors with S^-1/2. 158c 159 call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one, 160 & g_sm12, g_moso(1), zero, g_fockso(1)) 161 call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one, 162 & g_sm12, g_moso(2), zero, g_fockso(2)) 163c 164c Transfer into MO vec arrays 165c 166 call ga_zero(g_moso(1)) 167 call ga_zero(g_moso(2)) 168 call ga_copy(g_fockso(1), g_moso(1)) 169 call ga_copy(g_fockso(2), g_moso(2)) 170c 171 return 172 end 173c 174c Calculate S powers: S,S^(-1/2),S^(1/2),S^(-1) 175c 176 subroutine calc_s_powers(g_scr,g_tmp,nbf_ao,toll_s,svals,g_svecs, 177 & g_sp1,g_sm12,g_sp12,g_sm1) 178c 179 implicit none 180c 181#include "errquit.fh" 182#include "global.fh" 183#include "stdio.fh" 184#include "consts.fh" 185c 186 integer iw 187 integer g_scr 188 integer g_tmp 189 integer nbf_ao 190 double precision toll_s 191 double precision svals(*) 192 integer g_svecs 193 integer g_sp1 194 integer g_sm12 195 integer g_sp12 196 integer g_sm1 197c 198c Calculate S^1: g_sp1 199c 200 iw = 1 201 call ga_zero(g_scr) 202 call diis_bld12_so(toll_s, svals, g_svecs, g_scr, 203 & g_tmp, nbf_ao, iw) 204 call ga_zero(g_sp1) 205 call ga_fock_sf(g_scr, g_sp1, nbf_ao) ! Map onto larger array 206c 207c Calculate S^-1/2: g_sm12 208c 209 iw = 2 210 call ga_zero(g_scr) 211 call diis_bld12_so(toll_s, svals, g_svecs, g_scr, 212 & g_tmp, nbf_ao, iw) 213 call ga_zero(g_sm12) 214 call ga_fock_sf(g_scr, g_sm12, nbf_ao) ! Map onto larger array 215c 216c Calculate S^+1/2: g_sp12 217c 218 iw = 3 219 call ga_zero(g_scr) 220 call diis_bld12_so(toll_s, svals, g_svecs, g_scr, 221 & g_tmp, nbf_ao, iw) 222 call ga_zero(g_sp12) 223 call ga_fock_sf(g_scr, g_sp12, nbf_ao) ! Map onto larger array 224c 225c Calculate S^-1: g_sm1 226c 227 iw = 4 228 call ga_zero(g_scr) 229 call diis_bld12_so(toll_s, svals, g_svecs, g_scr, 230 & g_tmp, nbf_ao, iw) 231 call ga_zero(g_sm1) 232 call ga_fock_sf(g_scr, g_sm1, nbf_ao) ! Map onto larger array 233c 234 return 235 end 236c 237c Print energies 238c 239 subroutine print_energies(etnew,enuc,ecore,ecoul,exc, 240 & nexc,rho_n,dft_time) 241c 242 implicit none 243c 244#include "global.fh" 245#include "errquit.fh" 246#include "stdio.fh" 247c 248 double precision etnew 249 double precision enuc 250 double precision ecore 251 double precision ecoul 252 double precision exc(2) 253 double precision nexc 254 double precision rho_n 255 double precision dft_time 256c 257c Print 258c 259 if (nexc.le.1) then 260 write(luout,222)etnew+enuc,ecore,ecoul,exc(1),enuc 261 else 262 write(luout,223)etnew+enuc,ecore,ecoul,exc(1),exc(2),enuc 263 end if 264 write(luout,2222) rho_n 265 write(luout,2223) dft_time 266c 267 222 format(// 268 & ' Total SO-DFT energy =', f22.12/ 269 & ' One electron energy =', f22.12/ 270 & ' Coulomb energy =', f22.12/ 271 & ' Exchange-Corr. energy =', f22.12/ 272 & ' Nuclear repulsion energy =', f22.12/) 273c 274 223 format(// 275 & ' Total SO-DFT energy =', f22.12/ 276 & ' One electron energy =', f22.12/ 277 & ' Coulomb energy =', f22.12/ 278 & ' Exchange energy =', f22.12/ 279 & ' Correlation energy =', f22.12/ 280 & ' Nuclear repulsion energy =', f22.12/) 281c 282 2222 format(' Numeric. integr. density =', f22.12/) 283 2223 format(' Total iterative time =', f9.1,'s'//) 284c 285 return 286 end 287c 288c Level shifting is implemented here (similarity 289c transformation before standard eigensolver). Note, 290c levelshifting is appropriate once a transformation 291c is available which makes the resulting Fock matrix 292c diagonally dominant, e.g., in an approximate MO basis. 293c Also note, there are many matrix multiplies with S^+-1/2 294c which are redundant if one is sure that the former basis 295c is orthonormal. 296c 297 subroutine levelshift_fock(nbf_mo,ntotocc,g_tmp2,g_sp12,g_moso, 298 & g_scr2,g_fockso,rlshift) 299c 300 implicit none 301c 302#include "global.fh" 303#include "errquit.fh" 304#include "stdio.fh" 305#include "consts.fh" 306c 307 integer nbf_mo 308 integer ntotocc 309 integer g_tmp2 310 integer g_sp12 311 integer g_moso(2) 312 integer g_scr2 313 integer g_fockso(2) 314 double precision rlshift 315c 316 integer j 317 integer me 318 integer nproc 319c 320c Preliminaries 321c 322 me = ga_nodeid() 323 nproc = ga_nnodes() 324c 325c Build a matrix which is diagonal in the "MO" rep, 326c back-transform, and shift the current Fock matrix 327c 328c Use S^+1/2 (g_sp12) * old movecs (as a transform). 329c 330c Real part 331 call ga_zero(g_tmp2) 332 call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one, 333 & g_sp12, g_moso(1), zero, g_tmp2) 334 call ga_copy(g_tmp2, g_moso(1)) 335 336c Imag part 337 call ga_zero(g_tmp2) 338 call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one, 339 & g_sp12, g_moso(2), zero, g_tmp2) 340 call ga_copy(g_tmp2, g_moso(2)) 341c 342c Build diagonal matrix for the shift 343 call ga_zero(g_tmp2) 344 do j = nTotOcc+1+me, nbf_mo, nproc 345 call ga_put(g_tmp2, j, j, j, j, rlshift, 1) 346 enddo 347c 348c Transform this into "AO" basis and add to current 349c Fock matrix 350c 351c Real part 352 call ga_zero(g_scr2) ! used as a work area 353 call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one, 354 & g_moso(1), g_tmp2, zero, g_scr2) 355 call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, one, 356 & g_scr2, g_moso(1), one, g_fockso(1)) 357 call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one, 358 & g_moso(2), g_tmp2, zero, g_scr2) 359 call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, one, 360 & g_scr2, g_moso(2), one, g_fockso(1)) 361c 362c Imag part 363 call ga_zero(g_scr2) ! used as a work area 364 call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one, 365 & g_moso(1), g_tmp2, zero, g_scr2) 366 call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, mone, 367 & g_scr2, g_moso(2), one, g_fockso(2)) 368 call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one, 369 & g_moso(2), g_tmp2, zero, g_scr2) 370 call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, one, 371 & g_scr2, g_moso(1), one, g_fockso(2)) 372c 373 return 374 end 375c $Id$ 376