|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00007 MODULE atom_electronic_structure 00008 USE atom_optimization, ONLY: atom_history_init,& 00009 atom_history_release,& 00010 atom_history_type,& 00011 atom_history_update,& 00012 atom_opt 00013 USE atom_output, ONLY: atom_print_energies,& 00014 atom_print_iteration,& 00015 atom_print_state 00016 USE atom_types, ONLY: & 00017 GTH_PSEUDO, NO_PSEUDO, atom_type, create_opgrid, create_opmat, & 00018 opgrid_type, opmat_type, release_opgrid, release_opmat 00019 USE atom_utils, ONLY: & 00020 atom_denmat, atom_density, atom_solve, atom_trace, ceri_contract, & 00021 coulomb_potential_analytic, coulomb_potential_numeric, eeri_contract, & 00022 err_matrix, exchange_numeric, exchange_semi_analytic, numpot_matrix, & 00023 slater_density, wigner_slater_functional 00024 USE atom_xc, ONLY: calculate_atom_vxc_lda,& 00025 calculate_atom_vxc_lsd 00026 USE f77_blas 00027 USE input_constants, ONLY: & 00028 do_analytic, do_dkh0_atom, do_dkh1_atom, do_dkh2_atom, do_dkh3_atom, & 00029 do_dkh4_atom, do_dkh5_atom, do_nonrel_atom, do_numeric, do_rhf_atom, & 00030 do_rks_atom, do_rohf_atom, do_semi_analytic, do_uhf_atom, do_uks_atom, & 00031 do_zoramp_atom 00032 USE input_section_types, ONLY: section_vals_get,& 00033 section_vals_get_subs_vals,& 00034 section_vals_type,& 00035 section_vals_val_get 00036 USE kinds, ONLY: dp 00037 USE timings, ONLY: timeset,& 00038 timestop 00039 #include "cp_common_uses.h" 00040 00041 IMPLICIT NONE 00042 PRIVATE 00043 PUBLIC :: calculate_atom 00044 00045 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_electronic_structure' 00046 00047 ! ***************************************************************************** 00048 00049 CONTAINS 00050 00051 ! ***************************************************************************** 00052 SUBROUTINE calculate_atom(atom,iw,noguess,error) 00053 TYPE(atom_type), POINTER :: atom 00054 INTEGER, INTENT(IN) :: iw 00055 LOGICAL, INTENT(IN), OPTIONAL :: noguess 00056 TYPE(cp_error_type), INTENT(INOUT) :: error 00057 00058 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_atom', 00059 routineP = moduleN//':'//routineN 00060 00061 INTEGER :: handle, method 00062 LOGICAL :: failure 00063 00064 failure = .FALSE. 00065 00066 CALL timeset(routineN,handle) 00067 00068 method = atom%method_type 00069 00070 SELECT CASE (method) 00071 CASE DEFAULT 00072 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00073 CASE (do_rks_atom) 00074 CALL calculate_atom_restricted(atom,iw,noguess,error) 00075 CASE (do_uks_atom) 00076 CALL calculate_atom_unrestricted(atom,iw,noguess,error) 00077 CASE (do_rhf_atom) 00078 CALL calculate_atom_restricted(atom,iw,noguess,error) 00079 CASE (do_uhf_atom) 00080 CALL calculate_atom_unrestricted(atom,iw,noguess,error) 00081 CASE (do_rohf_atom) 00082 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00083 END SELECT 00084 00085 CALL timestop(handle) 00086 00087 END SUBROUTINE calculate_atom 00088 00089 ! ***************************************************************************** 00090 SUBROUTINE calculate_atom_restricted(atom,iw,noguess,error) 00091 TYPE(atom_type), POINTER :: atom 00092 INTEGER, INTENT(IN) :: iw 00093 LOGICAL, INTENT(IN), OPTIONAL :: noguess 00094 TYPE(cp_error_type), INTENT(INOUT) :: error 00095 00096 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_atom_restricted', 00097 routineP = moduleN//':'//routineN 00098 00099 INTEGER :: handle, i, iter, l, max_iter, 00100 method, reltyp 00101 LOGICAL :: do_hfx, doguess, failure, 00102 need_x, need_xc 00103 REAL(KIND=dp) :: deps, eps_scf, hf_frac 00104 TYPE(atom_history_type) :: history 00105 TYPE(opgrid_type), POINTER :: cpot, density 00106 TYPE(opmat_type), POINTER :: fmat, hcore, jmat, kmat, xcmat 00107 TYPE(section_vals_type), POINTER :: hfx_sections, xc_section 00108 00109 failure = .FALSE. 00110 00111 CALL timeset(routineN,handle) 00112 00113 IF (PRESENT(noguess)) THEN 00114 doguess=.NOT.noguess 00115 ELSE 00116 doguess=.TRUE. 00117 END IF 00118 00119 hf_frac = 0._dp 00120 IF (ASSOCIATED(atom%xc_section)) THEN 00121 xc_section => atom%xc_section 00122 hfx_sections => section_vals_get_subs_vals(xc_section,"HF",error=error) 00123 CALL section_vals_get(hfx_sections,explicit=do_hfx,error=error) 00124 IF ( do_hfx ) THEN 00125 CALL section_vals_val_get(hfx_sections,"FRACTION", r_val=hf_frac, error=error) 00126 END IF 00127 ELSE 00128 NULLIFY(xc_section) 00129 do_hfx = .FALSE. 00130 END IF 00131 00132 method = atom%method_type 00133 max_iter = atom%optimization%max_iter 00134 eps_scf = atom%optimization%eps_scf 00135 00136 SELECT CASE (method) 00137 CASE DEFAULT 00138 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00139 CASE (do_rks_atom) 00140 need_x = do_hfx 00141 need_xc = .TRUE. 00142 CASE (do_uks_atom) 00143 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00144 CASE (do_rhf_atom) 00145 need_x = .TRUE. 00146 need_xc = .FALSE. 00147 hf_frac = 1._dp 00148 CASE (do_uhf_atom) 00149 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00150 CASE (do_rohf_atom) 00151 need_x = .TRUE. 00152 need_xc = .FALSE. 00153 hf_frac = 1._dp 00154 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00155 END SELECT 00156 00157 ! check for relativistic method 00158 reltyp = atom%relativistic 00159 00160 IF (iw>0) CALL atom_print_state(atom%state,iw,error) 00161 00162 NULLIFY(hcore) 00163 CALL create_opmat(hcore,atom%basis%nbas,error) 00164 ! Pseudopotentials 00165 SELECT CASE (atom%potential%ppot_type) 00166 CASE DEFAULT 00167 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00168 CASE (NO_PSEUDO) 00169 SELECT CASE (reltyp) 00170 CASE DEFAULT 00171 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00172 CASE (do_nonrel_atom) 00173 hcore%op = atom%integrals%kin - atom%zcore*atom%integrals%core 00174 CASE (do_zoramp_atom) 00175 hcore%op = atom%integrals%kin + atom%integrals%tzora - atom%zcore*atom%integrals%core 00176 CASE (do_dkh0_atom, do_dkh1_atom,do_dkh2_atom,do_dkh3_atom,do_dkh4_atom,do_dkh5_atom) 00177 hcore%op = atom%integrals%hdkh 00178 END SELECT 00179 CASE (GTH_PSEUDO) 00180 hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl 00181 END SELECT 00182 ! add confinement potential (not included in relativistic transformations) 00183 IF ( atom%potential%confinement ) THEN 00184 hcore%op = hcore%op + atom%potential%acon * atom%integrals%conf 00185 END IF 00186 00187 NULLIFY(fmat,jmat,kmat,xcmat) 00188 CALL create_opmat(fmat,atom%basis%nbas,error) 00189 CALL create_opmat(jmat,atom%basis%nbas,error) 00190 CALL create_opmat(kmat,atom%basis%nbas,error) 00191 CALL create_opmat(xcmat,atom%basis%nbas,error) 00192 00193 NULLIFY(density,cpot) 00194 CALL create_opgrid(density,atom%basis%grid,error) 00195 CALL create_opgrid(cpot,atom%basis%grid,error) 00196 00197 IF (doguess) THEN 00198 ! initial guess 00199 CALL slater_density(density%op,density%op,atom%z,atom%state,atom%basis%grid,error) 00200 CALL coulomb_potential_numeric(cpot%op,density%op,density%grid,error) 00201 CALL numpot_matrix(jmat%op,cpot%op,atom%basis,0,error) 00202 CALL wigner_slater_functional(density%op,cpot%op,error) 00203 CALL numpot_matrix(xcmat%op,cpot%op,atom%basis,0,error) 00204 fmat%op = hcore%op + jmat%op + xcmat%op 00205 CALL atom_solve(fmat%op,atom%integrals%utrans,atom%orbitals%wfn,atom%orbitals%ener,& 00206 atom%basis%nbas,atom%integrals%nne,atom%state%maxl_calc,error) 00207 END IF 00208 CALL atom_denmat(atom%orbitals%pmat,atom%orbitals%wfn,atom%basis%nbas,atom%state%occupation,& 00209 atom%state%maxl_occ,atom%state%maxn_occ,error) 00210 00211 ! wavefunction history 00212 NULLIFY(history%dmat,history%hmat) 00213 CALL atom_history_init (history,atom%optimization,fmat%op,error) 00214 00215 iter = 0 00216 DO !SCF Loop 00217 00218 ! Kinetic energy 00219 atom%energy%ekin = atom_trace(atom%integrals%kin,atom%orbitals%pmat,error) 00220 00221 ! Band energy 00222 atom%energy%eband = 0._dp 00223 DO l=0,3 00224 DO i=1,MIN(SIZE(atom%state%occupation,2),SIZE(atom%orbitals%ener,1)) 00225 atom%energy%eband = atom%energy%eband + atom%orbitals%ener(i,l)*atom%state%occupation(l,i) 00226 END DO 00227 END DO 00228 00229 ! Pseudopotential energy 00230 SELECT CASE (atom%potential%ppot_type) 00231 CASE DEFAULT 00232 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00233 CASE (NO_PSEUDO) 00234 atom%energy%eploc = 0._dp 00235 atom%energy%epnl = 0._dp 00236 CASE (GTH_PSEUDO) 00237 atom%energy%eploc = atom_trace(atom%integrals%core,atom%orbitals%pmat,error) 00238 atom%energy%epnl = atom_trace(atom%integrals%hnl,atom%orbitals%pmat,error) 00239 END SELECT 00240 atom%energy%epseudo = atom%energy%eploc + atom%energy%epnl 00241 00242 ! Core energy 00243 atom%energy%ecore = atom_trace(hcore%op,atom%orbitals%pmat,error) 00244 00245 ! Confinement energy 00246 IF ( atom%potential%confinement ) THEN 00247 atom%energy%econfinement = atom_trace(atom%integrals%conf,atom%orbitals%pmat,error) 00248 ELSE 00249 atom%energy%econfinement = 0._dp 00250 END IF 00251 00252 ! Hartree Term 00253 jmat%op = 0._dp 00254 SELECT CASE (atom%coulomb_integral_type) 00255 CASE DEFAULT 00256 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00257 CASE (do_analytic) 00258 CALL ceri_contract(jmat%op,atom%integrals%ceri,atom%orbitals%pmat,atom%integrals%n,error=error) 00259 CASE (do_semi_analytic) 00260 CALL coulomb_potential_analytic(cpot%op,atom%orbitals%pmat,atom%basis,atom%basis%grid,& 00261 atom%state%maxl_occ,error) 00262 CALL numpot_matrix(jmat%op,cpot%op,atom%basis,0,error) 00263 CASE (do_numeric) 00264 CALL atom_density(density%op,atom%orbitals%pmat,atom%basis,atom%state%maxl_occ,typ="RHO",error=error) 00265 CALL coulomb_potential_numeric(cpot%op,density%op,density%grid,error) 00266 CALL numpot_matrix(jmat%op,cpot%op,atom%basis,0,error) 00267 END SELECT 00268 atom%energy%ecoulomb = 0.5_dp * atom_trace(jmat%op,atom%orbitals%pmat,error) 00269 00270 ! Exchange Term 00271 IF (need_x) THEN 00272 kmat%op = 0._dp 00273 SELECT CASE (atom%exchange_integral_type) 00274 CASE DEFAULT 00275 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00276 CASE (do_analytic) 00277 CALL eeri_contract(kmat%op,atom%integrals%eeri,atom%orbitals%pmat,atom%integrals%n,error=error) 00278 CASE (do_semi_analytic) 00279 CALL exchange_semi_analytic(kmat%op,atom%state,atom%state%occupation,atom%orbitals%wfn,atom%basis,error=error) 00280 CASE (do_numeric) 00281 CALL exchange_numeric(kmat%op,atom%state,atom%state%occupation,atom%orbitals%wfn,atom%basis,error=error) 00282 END SELECT 00283 atom%energy%eexchange = hf_frac * 0.5_dp * atom_trace(kmat%op,atom%orbitals%pmat,error) 00284 kmat%op = hf_frac*kmat%op 00285 ELSE 00286 kmat%op = 0._dp 00287 atom%energy%eexchange = 0._dp 00288 END IF 00289 00290 ! XC 00291 IF (need_xc) THEN 00292 xcmat%op = 0._dp 00293 CALL calculate_atom_vxc_lda(xcmat,atom,xc_section,error) 00294 ELSE 00295 xcmat%op = 0._dp 00296 atom%energy%exc = 0._dp 00297 END IF 00298 00299 ! Zero this contribution 00300 atom%energy%elsd=0._dp 00301 00302 ! Total energy 00303 atom%energy%etot = atom%energy%ecore + atom%energy%ecoulomb + atom%energy%eexchange + atom%energy%exc 00304 00305 ! Potential energy 00306 atom%energy%epot = atom%energy%etot - atom%energy%ekin 00307 00308 ! Total HF/KS matrix 00309 fmat%op = hcore%op + jmat%op + kmat%op + xcmat%op 00310 00311 ! calculate error matrix 00312 CALL err_matrix(jmat%op,deps,fmat%op,atom%orbitals%pmat,atom%integrals%utrans,& 00313 atom%integrals%uptrans,atom%basis%nbas,atom%integrals%nne,error) 00314 00315 iter = iter + 1 00316 00317 IF ( iw > 0 ) THEN 00318 CALL atom_print_iteration(iter,deps,atom%energy%etot,iw,error) 00319 END IF 00320 00321 IF ( deps < eps_scf ) EXIT 00322 IF ( iter >= max_iter ) THEN 00323 IF ( iw > 0 ) THEN 00324 WRITE(iw,"(A)") " No convergence within maximum number of iterations " 00325 END IF 00326 EXIT 00327 END IF 00328 00329 ! update history container and extrapolate KS matrix 00330 CALL atom_history_update (history,fmat%op,jmat%op,error) 00331 CALL atom_opt (fmat%op,history,deps,error) 00332 00333 ! Solve HF/KS equations 00334 CALL atom_solve(fmat%op,atom%integrals%utrans,atom%orbitals%wfn,atom%orbitals%ener,& 00335 atom%basis%nbas,atom%integrals%nne,atom%state%maxl_calc,error) 00336 CALL atom_denmat(atom%orbitals%pmat,atom%orbitals%wfn,atom%basis%nbas,atom%state%occupation,& 00337 atom%state%maxl_occ,atom%state%maxn_occ,error) 00338 00339 END DO !SCF Loop 00340 00341 IF ( iw > 0 ) THEN 00342 CALL atom_print_energies(atom,iw,error) 00343 END IF 00344 00345 CALL atom_history_release(history,error) 00346 00347 CALL release_opmat(fmat,error) 00348 CALL release_opmat(jmat,error) 00349 CALL release_opmat(kmat,error) 00350 CALL release_opmat(xcmat,error) 00351 CALL release_opmat(hcore,error) 00352 00353 CALL release_opgrid(density,error) 00354 CALL release_opgrid(cpot,error) 00355 00356 CALL timestop(handle) 00357 00358 END SUBROUTINE calculate_atom_restricted 00359 00360 ! ***************************************************************************** 00361 SUBROUTINE calculate_atom_unrestricted(atom,iw,noguess,error) 00362 TYPE(atom_type), POINTER :: atom 00363 INTEGER, INTENT(IN) :: iw 00364 LOGICAL, INTENT(IN), OPTIONAL :: noguess 00365 TYPE(cp_error_type), INTENT(INOUT) :: error 00366 00367 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_atom_unrestricted', 00368 routineP = moduleN//':'//routineN 00369 00370 INTEGER :: handle, i, iter, k, l, 00371 max_iter, method, reltyp 00372 LOGICAL :: do_hfx, doguess, failure, 00373 lsdpot, need_x, need_xc 00374 REAL(KIND=dp) :: deps, depsa, depsb, eps_scf, 00375 hf_frac, ne, nm 00376 TYPE(atom_history_type) :: historya, historyb 00377 TYPE(opgrid_type), POINTER :: cpot, density, rhoa, rhob 00378 TYPE(opmat_type), POINTER :: fmata, fmatb, hcore, hlsd, 00379 jmat, kmata, kmatb, xcmata, 00380 xcmatb 00381 TYPE(section_vals_type), POINTER :: hfx_sections, xc_section 00382 00383 failure = .FALSE. 00384 00385 CALL timeset(routineN,handle) 00386 00387 IF (PRESENT(noguess)) THEN 00388 doguess=.NOT.noguess 00389 ELSE 00390 doguess=.TRUE. 00391 END IF 00392 00393 hf_frac = 0._dp 00394 IF (ASSOCIATED(atom%xc_section)) THEN 00395 xc_section => atom%xc_section 00396 hfx_sections => section_vals_get_subs_vals(xc_section,"HF",error=error) 00397 CALL section_vals_get(hfx_sections,explicit=do_hfx,error=error) 00398 IF ( do_hfx ) THEN 00399 CALL section_vals_val_get(hfx_sections,"FRACTION", r_val=hf_frac, error=error) 00400 END IF 00401 ELSE 00402 NULLIFY(xc_section) 00403 do_hfx = .FALSE. 00404 END IF 00405 00406 method = atom%method_type 00407 max_iter = atom%optimization%max_iter 00408 eps_scf = atom%optimization%eps_scf 00409 00410 SELECT CASE (method) 00411 CASE DEFAULT 00412 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00413 CASE (do_rks_atom) 00414 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00415 CASE (do_uks_atom) 00416 need_x = do_hfx 00417 need_xc = .TRUE. 00418 CASE (do_rhf_atom) 00419 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00420 CASE (do_uhf_atom) 00421 need_x = .TRUE. 00422 need_xc = .FALSE. 00423 hf_frac = 1._dp 00424 CASE (do_rohf_atom) 00425 need_x = .TRUE. 00426 need_xc = .FALSE. 00427 hf_frac = 1._dp 00428 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00429 END SELECT 00430 00431 ! set alpha and beta occupations 00432 atom%state%occa = 0._dp 00433 atom%state%occb = 0._dp 00434 DO l=0,3 00435 nm=REAL((2*l+1),KIND=dp) 00436 DO k=1,10 00437 ne = atom%state%occupation(l,k) 00438 IF(ne==0._dp) THEN !empty shell 00439 EXIT !assume there are no holes 00440 ELSEIF(ne==2._dp*nm) THEN !closed shell 00441 atom%state%occa(l,k)=nm 00442 atom%state%occb(l,k)=nm 00443 ELSEIF(atom%state%multiplicity==-2) THEN !High spin case 00444 atom%state%occa(l,k)=MIN(ne,nm) 00445 atom%state%occb(l,k)=MAX(0._dp,ne-nm) 00446 ELSE 00447 atom%state%occa(l,k)=0.5_dp*(ne+atom%state%multiplicity-1._dp) 00448 atom%state%occb(l,k)=ne-atom%state%occa(l,k) 00449 END IF 00450 END DO 00451 END DO 00452 ! check for relativistic method 00453 reltyp = atom%relativistic 00454 00455 IF (iw>0) CALL atom_print_state(atom%state,iw,error) 00456 00457 NULLIFY(hcore,hlsd) 00458 CALL create_opmat(hcore,atom%basis%nbas,error) 00459 CALL create_opmat(hlsd,atom%basis%nbas,error) 00460 hlsd%op = 0._dp 00461 ! Pseudopotentials 00462 lsdpot=.FALSE. 00463 SELECT CASE (atom%potential%ppot_type) 00464 CASE DEFAULT 00465 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00466 CASE (NO_PSEUDO) 00467 SELECT CASE (reltyp) 00468 CASE DEFAULT 00469 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00470 CASE (do_nonrel_atom) 00471 hcore%op = atom%integrals%kin - atom%zcore*atom%integrals%core 00472 CASE (do_zoramp_atom) 00473 hcore%op = atom%integrals%kin + atom%integrals%tzora - atom%zcore*atom%integrals%core 00474 CASE (do_dkh0_atom, do_dkh1_atom,do_dkh2_atom,do_dkh3_atom,do_dkh4_atom,do_dkh5_atom) 00475 hcore%op = atom%integrals%hdkh 00476 END SELECT 00477 CASE (GTH_PSEUDO) 00478 hcore%op = atom%integrals%kin + atom%integrals%core + atom%integrals%hnl 00479 IF(atom%potential%gth_pot%lsdpot) THEN 00480 lsdpot=.TRUE. 00481 hlsd%op = atom%integrals%clsd 00482 END IF 00483 END SELECT 00484 ! add confinement potential (not included in relativistic transformations) 00485 IF ( atom%potential%confinement ) THEN 00486 hcore%op = hcore%op + atom%potential%acon * atom%integrals%conf 00487 END IF 00488 00489 NULLIFY(fmata,fmatb,jmat,kmata,kmatb,xcmata,xcmatb) 00490 CALL create_opmat(fmata,atom%basis%nbas,error) 00491 CALL create_opmat(fmatb,atom%basis%nbas,error) 00492 CALL create_opmat(jmat,atom%basis%nbas,error) 00493 CALL create_opmat(kmata,atom%basis%nbas,error) 00494 CALL create_opmat(kmatb,atom%basis%nbas,error) 00495 CALL create_opmat(xcmata,atom%basis%nbas,error) 00496 CALL create_opmat(xcmatb,atom%basis%nbas,error) 00497 00498 NULLIFY(density,rhoa,rhob,cpot) 00499 CALL create_opgrid(density,atom%basis%grid,error) 00500 CALL create_opgrid(rhoa,atom%basis%grid,error) 00501 CALL create_opgrid(rhob,atom%basis%grid,error) 00502 CALL create_opgrid(cpot,atom%basis%grid,error) 00503 00504 IF (doguess) THEN 00505 ! initial guess 00506 CALL slater_density(rhoa%op,rhob%op,atom%z,atom%state,atom%basis%grid,error) 00507 density%op = rhoa%op + rhob%op 00508 CALL coulomb_potential_numeric(cpot%op,density%op,density%grid,error) 00509 CALL numpot_matrix(jmat%op,cpot%op,atom%basis,0,error) 00510 ! alpha spin 00511 density%op = 2._dp*rhoa%op 00512 CALL wigner_slater_functional(density%op,cpot%op,error) 00513 CALL numpot_matrix(xcmata%op,cpot%op,atom%basis,0,error) 00514 fmata%op = hcore%op + hlsd%op + jmat%op + xcmata%op 00515 CALL atom_solve(fmata%op,atom%integrals%utrans,atom%orbitals%wfna,atom%orbitals%enera,& 00516 atom%basis%nbas,atom%integrals%nne,atom%state%maxl_calc,error) 00517 ! beta spin 00518 density%op = 2._dp*rhob%op 00519 CALL wigner_slater_functional(density%op,cpot%op,error) 00520 CALL numpot_matrix(xcmatb%op,cpot%op,atom%basis,0,error) 00521 fmatb%op = hcore%op - hlsd%op + jmat%op + xcmatb%op 00522 CALL atom_solve(fmatb%op,atom%integrals%utrans,atom%orbitals%wfnb,atom%orbitals%enerb,& 00523 atom%basis%nbas,atom%integrals%nne,atom%state%maxl_calc,error) 00524 END IF 00525 CALL atom_denmat(atom%orbitals%pmata,atom%orbitals%wfna,atom%basis%nbas,atom%state%occa,& 00526 atom%state%maxl_occ,atom%state%maxn_occ,error) 00527 CALL atom_denmat(atom%orbitals%pmatb,atom%orbitals%wfnb,atom%basis%nbas,atom%state%occb,& 00528 atom%state%maxl_occ,atom%state%maxn_occ,error) 00529 atom%orbitals%pmat = atom%orbitals%pmata + atom%orbitals%pmatb 00530 00531 ! wavefunction history 00532 NULLIFY(historya%dmat,historya%hmat) 00533 CALL atom_history_init (historya,atom%optimization,fmata%op,error) 00534 NULLIFY(historyb%dmat,historyb%hmat) 00535 CALL atom_history_init (historyb,atom%optimization,fmatb%op,error) 00536 00537 iter = 0 00538 DO !SCF Loop 00539 00540 ! Kinetic energy 00541 atom%energy%ekin = atom_trace(atom%integrals%kin,atom%orbitals%pmat,error) 00542 00543 ! Band energy 00544 atom%energy%eband = 0._dp 00545 DO l=0,3 00546 DO i=1,MIN(SIZE(atom%state%occupation,2),SIZE(atom%orbitals%ener,1)) 00547 atom%energy%eband = atom%energy%eband + atom%orbitals%enera(i,l)*atom%state%occa(l,i) 00548 atom%energy%eband = atom%energy%eband + atom%orbitals%enerb(i,l)*atom%state%occb(l,i) 00549 END DO 00550 END DO 00551 00552 ! Pseudopotential energy 00553 SELECT CASE (atom%potential%ppot_type) 00554 CASE DEFAULT 00555 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00556 CASE (NO_PSEUDO) 00557 atom%energy%eploc = 0._dp 00558 atom%energy%epnl = 0._dp 00559 CASE (GTH_PSEUDO) 00560 atom%energy%eploc = atom_trace(atom%integrals%core,atom%orbitals%pmat,error) 00561 atom%energy%epnl = atom_trace(atom%integrals%hnl,atom%orbitals%pmat,error) 00562 END SELECT 00563 atom%energy%epseudo = atom%energy%eploc + atom%energy%epnl 00564 00565 ! Core energy 00566 atom%energy%ecore = atom_trace(hcore%op,atom%orbitals%pmat,error) 00567 00568 ! Confinement energy 00569 IF ( atom%potential%confinement ) THEN 00570 atom%energy%econfinement = atom_trace(atom%integrals%conf,atom%orbitals%pmat,error) 00571 ELSE 00572 atom%energy%econfinement = 0._dp 00573 END IF 00574 00575 ! Hartree Term 00576 jmat%op = 0._dp 00577 SELECT CASE (atom%coulomb_integral_type) 00578 CASE DEFAULT 00579 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00580 CASE (do_analytic) 00581 CALL ceri_contract(jmat%op,atom%integrals%ceri,atom%orbitals%pmat,atom%integrals%n,error=error) 00582 CASE (do_semi_analytic) 00583 CALL coulomb_potential_analytic(cpot%op,atom%orbitals%pmat,atom%basis,atom%basis%grid,& 00584 atom%state%maxl_occ,error) 00585 CALL numpot_matrix(jmat%op,cpot%op,atom%basis,0,error) 00586 CASE (do_numeric) 00587 CALL atom_density(density%op,atom%orbitals%pmat,atom%basis,atom%state%maxl_occ,typ="RHO",error=error) 00588 CALL coulomb_potential_numeric(cpot%op,density%op,density%grid,error) 00589 CALL numpot_matrix(jmat%op,cpot%op,atom%basis,0,error) 00590 END SELECT 00591 atom%energy%ecoulomb = 0.5_dp * atom_trace(jmat%op,atom%orbitals%pmat,error) 00592 00593 ! Exchange Term 00594 IF (need_x) THEN 00595 kmata%op = 0._dp 00596 kmatb%op = 0._dp 00597 SELECT CASE (atom%exchange_integral_type) 00598 CASE DEFAULT 00599 CPPostcondition(.FALSE., cp_failure_level, routineP, error, failure) 00600 CASE (do_analytic) 00601 CALL eeri_contract(kmata%op,atom%integrals%eeri,atom%orbitals%pmata,atom%integrals%n,error=error) 00602 CALL eeri_contract(kmatb%op,atom%integrals%eeri,atom%orbitals%pmatb,atom%integrals%n,error=error) 00603 CASE (do_semi_analytic) 00604 CALL exchange_semi_analytic(kmata%op,atom%state,atom%state%occa,atom%orbitals%wfna,atom%basis,error=error) 00605 CALL exchange_semi_analytic(kmatb%op,atom%state,atom%state%occb,atom%orbitals%wfnb,atom%basis,error=error) 00606 CASE (do_numeric) 00607 CALL exchange_numeric(kmata%op,atom%state,atom%state%occa,atom%orbitals%wfna,atom%basis,error=error) 00608 CALL exchange_numeric(kmatb%op,atom%state,atom%state%occb,atom%orbitals%wfnb,atom%basis,error=error) 00609 END SELECT 00610 atom%energy%eexchange = hf_frac * ( atom_trace(kmata%op,atom%orbitals%pmata,error) + & 00611 atom_trace(kmatb%op,atom%orbitals%pmatb,error) ) 00612 kmata%op = 2._dp*hf_frac*kmata%op 00613 kmatb%op = 2._dp*hf_frac*kmatb%op 00614 ELSE 00615 kmata%op = 0._dp 00616 kmatb%op = 0._dp 00617 atom%energy%eexchange = 0._dp 00618 END IF 00619 00620 ! XC 00621 IF (need_xc) THEN 00622 xcmata%op = 0._dp 00623 xcmatb%op = 0._dp 00624 CALL calculate_atom_vxc_lsd(xcmata,xcmatb,atom,xc_section,error) 00625 ELSE 00626 xcmata%op = 0._dp 00627 xcmatb%op = 0._dp 00628 atom%energy%exc = 0._dp 00629 END IF 00630 00631 IF(lsdpot) THEN 00632 atom%energy%elsd = atom_trace(hlsd%op,atom%orbitals%pmata,error) - & 00633 atom_trace(hlsd%op,atom%orbitals%pmatb,error) 00634 atom%energy%epseudo = atom%energy%epseudo + atom%energy%elsd 00635 atom%energy%ecore = atom%energy%ecore + atom%energy%elsd 00636 ELSE 00637 atom%energy%elsd=0._dp 00638 END IF 00639 00640 ! Total energy 00641 atom%energy%etot = atom%energy%ecore + atom%energy%ecoulomb + atom%energy%eexchange + atom%energy%exc 00642 00643 ! Potential energy 00644 atom%energy%epot = atom%energy%etot - atom%energy%ekin 00645 00646 ! Total HF/KS matrix 00647 fmata%op = hcore%op + hlsd%op + jmat%op + kmata%op + xcmata%op 00648 fmatb%op = hcore%op - hlsd%op + jmat%op + kmatb%op + xcmatb%op 00649 00650 ! calculate error matrix 00651 CALL err_matrix(xcmata%op,depsa,fmata%op,atom%orbitals%pmata,atom%integrals%utrans,& 00652 atom%integrals%uptrans,atom%basis%nbas,atom%integrals%nne,error) 00653 CALL err_matrix(xcmatb%op,depsb,fmatb%op,atom%orbitals%pmatb,atom%integrals%utrans,& 00654 atom%integrals%uptrans,atom%basis%nbas,atom%integrals%nne,error) 00655 deps=2._dp*MAX(depsa,depsb) 00656 00657 iter = iter + 1 00658 00659 IF ( iw > 0 ) THEN 00660 CALL atom_print_iteration(iter,deps,atom%energy%etot,iw,error) 00661 END IF 00662 00663 IF ( deps < eps_scf ) EXIT 00664 IF ( iter >= max_iter ) THEN 00665 IF ( iw > 0 ) THEN 00666 WRITE(iw,"(A)") " No convergence within maximum number of iterations " 00667 END IF 00668 EXIT 00669 END IF 00670 00671 ! update history container and extrapolate KS matrix 00672 CALL atom_history_update (historya,fmata%op,xcmata%op,error) 00673 CALL atom_history_update (historyb,fmatb%op,xcmatb%op,error) 00674 CALL atom_opt (fmata%op,historya,depsa,error) 00675 CALL atom_opt (fmatb%op,historyb,depsb,error) 00676 00677 ! Solve HF/KS equations 00678 CALL atom_solve(fmata%op,atom%integrals%utrans,atom%orbitals%wfna,atom%orbitals%enera,& 00679 atom%basis%nbas,atom%integrals%nne,atom%state%maxl_calc,error) 00680 CALL atom_denmat(atom%orbitals%pmata,atom%orbitals%wfna,atom%basis%nbas,atom%state%occa,& 00681 atom%state%maxl_occ,atom%state%maxn_occ,error) 00682 CALL atom_solve(fmatb%op,atom%integrals%utrans,atom%orbitals%wfnb,atom%orbitals%enerb,& 00683 atom%basis%nbas,atom%integrals%nne,atom%state%maxl_calc,error) 00684 CALL atom_denmat(atom%orbitals%pmatb,atom%orbitals%wfnb,atom%basis%nbas,atom%state%occb,& 00685 atom%state%maxl_occ,atom%state%maxn_occ,error) 00686 atom%orbitals%pmat = atom%orbitals%pmata + atom%orbitals%pmatb 00687 00688 END DO !SCF Loop 00689 00690 IF ( iw > 0 ) THEN 00691 CALL atom_print_energies(atom,iw,error) 00692 END IF 00693 00694 CALL atom_history_release(historya,error) 00695 CALL atom_history_release(historyb,error) 00696 00697 CALL release_opgrid(density,error) 00698 CALL release_opgrid(rhoa,error) 00699 CALL release_opgrid(rhob,error) 00700 CALL release_opgrid(cpot,error) 00701 00702 CALL release_opmat(fmata,error) 00703 CALL release_opmat(fmatb,error) 00704 CALL release_opmat(jmat,error) 00705 CALL release_opmat(kmata,error) 00706 CALL release_opmat(kmatb,error) 00707 CALL release_opmat(xcmata,error) 00708 CALL release_opmat(xcmatb,error) 00709 CALL release_opmat(hlsd,error) 00710 CALL release_opmat(hcore,error) 00711 00712 CALL timestop(handle) 00713 00714 END SUBROUTINE calculate_atom_unrestricted 00715 00716 ! ***************************************************************************** 00717 00718 END MODULE atom_electronic_structure
1.7.3