|
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 ! ***************************************************************************** 00016 MODULE dft_plus_u 00017 00018 USE atomic_kind_types, ONLY: atomic_kind_type,& 00019 get_atomic_kind,& 00020 get_atomic_kind_set,& 00021 set_atomic_kind 00022 USE basis_set_types, ONLY: get_gto_basis_set,& 00023 gto_basis_set_type 00024 USE bibliography, ONLY: Dudarev1997,& 00025 Dudarev1998,& 00026 cite_reference 00027 USE cp_control_types, ONLY: dft_control_type 00028 USE cp_dbcsr_interface, ONLY: cp_dbcsr_get_block_diag,& 00029 cp_dbcsr_get_block_p,& 00030 cp_dbcsr_init,& 00031 cp_dbcsr_iterator_blocks_left,& 00032 cp_dbcsr_iterator_next_block,& 00033 cp_dbcsr_iterator_start,& 00034 cp_dbcsr_iterator_stop,& 00035 cp_dbcsr_set 00036 USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,& 00037 cp_dbcsr_deallocate_matrix,& 00038 cp_dbcsr_plus_fm_fm_t,& 00039 cp_dbcsr_sm_fm_multiply 00040 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix,& 00041 write_fm_with_basis_info 00042 USE cp_dbcsr_types, ONLY: cp_dbcsr_iterator,& 00043 cp_dbcsr_p_type,& 00044 cp_dbcsr_type 00045 USE cp_fm_basic_linalg, ONLY: cp_fm_gemm,& 00046 cp_fm_transpose 00047 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00048 cp_fm_struct_release,& 00049 cp_fm_struct_type 00050 USE cp_fm_types, ONLY: cp_fm_create,& 00051 cp_fm_release,& 00052 cp_fm_type 00053 USE cp_output_handling, ONLY: cp_p_file,& 00054 cp_print_key_finished_output,& 00055 cp_print_key_should_output,& 00056 cp_print_key_unit_nr 00057 USE f77_blas 00058 USE input_constants, ONLY: low_print_level,& 00059 plus_u_lowdin,& 00060 plus_u_mulliken,& 00061 plus_u_mulliken_charges 00062 USE input_section_types, ONLY: section_vals_type 00063 USE kinds, ONLY: default_string_length,& 00064 dp 00065 USE mathlib, ONLY: jacobi 00066 USE message_passing, ONLY: mp_sum 00067 USE orbital_symbols, ONLY: sgf_symbol 00068 USE particle_types, ONLY: get_particle_set,& 00069 particle_type 00070 USE physcon, ONLY: evolt 00071 USE qs_energy_types, ONLY: qs_energy_type 00072 USE qs_environment_types, ONLY: get_qs_env,& 00073 qs_environment_type 00074 USE qs_rho_types, ONLY: qs_rho_type 00075 USE qs_scf_types, ONLY: qs_scf_env_type 00076 USE timings, ONLY: timeset,& 00077 timestop 00078 #include "cp_common_uses.h" 00079 00080 IMPLICIT NONE 00081 00082 PRIVATE 00083 00084 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dft_plus_u' 00085 00086 PUBLIC :: plus_u 00087 00088 CONTAINS 00089 ! ***************************************************************************** 00100 SUBROUTINE plus_u(qs_env,matrix_h,matrix_w,error) 00101 00102 TYPE(qs_environment_type), POINTER :: qs_env 00103 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00104 OPTIONAL, POINTER :: matrix_h, matrix_w 00105 TYPE(cp_error_type), INTENT(INOUT) :: error 00106 00107 CHARACTER(LEN=*), PARAMETER :: routineN = 'plus_u', 00108 routineP = moduleN//':'//routineN 00109 00110 INTEGER :: handle, output_unit, 00111 print_level 00112 LOGICAL :: failure, orthonormal_basis, 00113 should_output 00114 TYPE(cp_logger_type), POINTER :: logger 00115 TYPE(section_vals_type), POINTER :: input 00116 00117 CALL timeset(routineN,handle) 00118 00119 failure = .FALSE. 00120 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 00121 00122 NULLIFY (input) 00123 00124 logger => cp_error_get_logger(error) 00125 00126 CALL get_qs_env(qs_env=qs_env,& 00127 input=input,& 00128 error=error) 00129 00130 CALL cite_reference(Dudarev1997) 00131 CALL cite_reference(Dudarev1998) 00132 00133 ! Later we could save here some time, if the method in use has this property 00134 ! which then has to be figured out here 00135 00136 orthonormal_basis = .FALSE. 00137 00138 ! Setup print control 00139 00140 print_level = logger%iter_info%print_level 00141 should_output = (BTEST(cp_print_key_should_output(logger%iter_info,input,& 00142 "DFT%PRINT%PLUS_U",error=error),cp_p_file).AND.& 00143 (.NOT.PRESENT(matrix_w))) 00144 output_unit = cp_print_key_unit_nr(logger,input,"DFT%PRINT%PLUS_U",& 00145 extension=".plus_u",& 00146 ignore_should_output=should_output,& 00147 log_filename=.FALSE.,& 00148 error=error) 00149 00150 ! Select DFT+U method 00151 00152 SELECT CASE (qs_env%dft_control%plus_u_method_id) 00153 CASE (plus_u_lowdin) 00154 IF (orthonormal_basis) THEN 00155 ! For an orthonormal basis the Lowdin method and the Mulliken method 00156 ! are equivalent 00157 CALL mulliken(qs_env,orthonormal_basis,matrix_h,& 00158 should_output,output_unit,print_level,error) 00159 ELSE 00160 CALL lowdin(qs_env,matrix_h,matrix_w,& 00161 should_output,output_unit,print_level,error) 00162 END IF 00163 CASE (plus_u_mulliken) 00164 CALL mulliken(qs_env,orthonormal_basis,matrix_h,& 00165 should_output,output_unit,print_level,error) 00166 CASE (plus_u_mulliken_charges) 00167 CALL mulliken_charges(qs_env,orthonormal_basis,matrix_h,matrix_w,& 00168 should_output,output_unit,print_level,error) 00169 CASE DEFAULT 00170 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 00171 "Invalid DFT+U method requested",only_ionode=.TRUE.) 00172 END SELECT 00173 00174 CALL cp_print_key_finished_output(output_unit,logger,input,"DFT%PRINT%PLUS_U",& 00175 ignore_should_output=should_output,& 00176 error=error) 00177 00178 CALL timestop(handle) 00179 00180 END SUBROUTINE plus_u 00181 00182 ! ***************************************************************************** 00209 SUBROUTINE lowdin(qs_env,matrix_h,matrix_w,should_output,output_unit,& 00210 print_level,error) 00211 00212 TYPE(qs_environment_type), POINTER :: qs_env 00213 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00214 OPTIONAL, POINTER :: matrix_h, matrix_w 00215 LOGICAL, INTENT(IN) :: should_output 00216 INTEGER, INTENT(IN) :: output_unit, print_level 00217 TYPE(cp_error_type), INTENT(INOUT) :: error 00218 00219 CHARACTER(LEN=*), PARAMETER :: routineN = 'lowdin', 00220 routineP = moduleN//':'//routineN 00221 00222 CHARACTER(LEN=10) :: spin_info 00223 CHARACTER(LEN=6), ALLOCATABLE, 00224 DIMENSION(:) :: symbol 00225 CHARACTER(LEN=default_string_length) :: atomic_kind_name 00226 INTEGER :: atom_a, handle, i, i0, iatom, ikind, iorb, isb, iset, isgf, 00227 ishell, ispin, j, jsb, jset, jsgf, jshell, lu, m, max_scf, n, natom, 00228 natom_of_kind, nkind, norb, nsb, nsbsize, nset, nsgf, nsgf_kind, nspin, 00229 stat 00230 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom 00231 INTEGER, DIMENSION(1) :: iloc 00232 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals 00233 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf 00234 LOGICAL :: debug, dft_plus_u_atom, failure, fm_work1_local_alloc, 00235 fm_work2_local_alloc, found, just_energy, smear 00236 LOGICAL, ALLOCATABLE, DIMENSION(:) :: orb_occ 00237 REAL(KIND=dp) :: eps_scf, eps_u_ramping, 00238 fspin, occ, trq, trq2, 00239 u_minus_j, u_minus_j_target, 00240 u_ramping 00241 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: q_eigval 00242 REAL(KIND=dp), ALLOCATABLE, 00243 DIMENSION(:, :) :: q_eigvec, q_matrix, q_work 00244 REAL(KIND=dp), DIMENSION(:, :), POINTER :: q_block, v_block 00245 TYPE(atomic_kind_type), DIMENSION(:), 00246 POINTER :: atomic_kind_set 00247 TYPE(atomic_kind_type), POINTER :: atomic_kind 00248 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00249 POINTER :: matrix_p, matrix_s 00250 TYPE(cp_dbcsr_type), POINTER :: sm_h, sm_p, sm_q, sm_s, sm_v, 00251 sm_w 00252 TYPE(cp_fm_struct_type), POINTER :: fmstruct 00253 TYPE(cp_fm_type), POINTER :: fm_s_half, fm_work1, fm_work2 00254 TYPE(dft_control_type), POINTER :: dft_control 00255 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 00256 TYPE(particle_type), DIMENSION(:), 00257 POINTER :: particle_set 00258 TYPE(qs_energy_type), POINTER :: energy 00259 TYPE(qs_rho_type), POINTER :: rho 00260 TYPE(qs_scf_env_type), POINTER :: scf_env 00261 00262 CALL timeset(routineN,handle) 00263 00264 debug = .FALSE. ! Set to .TRUE. to print debug information 00265 00266 NULLIFY (atom_list) 00267 NULLIFY (atomic_kind) 00268 NULLIFY (atomic_kind_set) 00269 NULLIFY (dft_control) 00270 NULLIFY (energy) 00271 NULLIFY (first_sgf) 00272 NULLIFY (fm_s_half) 00273 NULLIFY (fm_work1) 00274 NULLIFY (fm_work2) 00275 NULLIFY (fmstruct) 00276 NULLIFY (matrix_p) 00277 NULLIFY (matrix_s) 00278 NULLIFY (l) 00279 NULLIFY (last_sgf) 00280 NULLIFY (nshell) 00281 NULLIFY (orb_basis_set) 00282 NULLIFY (orbitals) 00283 NULLIFY (particle_set) 00284 NULLIFY (q_block) 00285 NULLIFY (rho) 00286 NULLIFY (scf_env) 00287 NULLIFY (sm_h) 00288 NULLIFY (sm_p) 00289 NULLIFY (sm_q) 00290 NULLIFY (sm_s) 00291 NULLIFY (sm_v) 00292 NULLIFY (sm_w) 00293 NULLIFY (v_block) 00294 00295 failure = .FALSE. 00296 00297 smear = .FALSE. 00298 max_scf = -1 00299 eps_scf = 1.0E30_dp 00300 00301 CALL get_qs_env(qs_env=qs_env,& 00302 atomic_kind_set=atomic_kind_set,& 00303 dft_control=dft_control,& 00304 energy=energy,& 00305 matrix_s=matrix_s,& 00306 particle_set=particle_set,& 00307 rho=rho,& 00308 scf_env=scf_env,& 00309 error=error) 00310 00311 CPPrecondition(ASSOCIATED(atomic_kind_set),cp_failure_level,routineP,error,failure) 00312 CPPrecondition(ASSOCIATED(dft_control),cp_failure_level,routineP,error,failure) 00313 CPPrecondition(ASSOCIATED(energy),cp_failure_level,routineP,error,failure) 00314 CPPrecondition(ASSOCIATED(matrix_s),cp_failure_level,routineP,error,failure) 00315 CPPrecondition(ASSOCIATED(particle_set),cp_failure_level,routineP,error,failure) 00316 CPPrecondition(ASSOCIATED(rho),cp_failure_level,routineP,error,failure) 00317 00318 sm_s => matrix_s(1)%matrix ! Overlap matrix in sparse format 00319 matrix_p => rho%rho_ao ! Density matrices in sparse format 00320 00321 energy%dft_plus_u = 0.0_dp 00322 00323 nspin = dft_control%nspins 00324 00325 IF (nspin == 2) THEN 00326 fspin = 1.0_dp 00327 ELSE 00328 fspin = 0.5_dp 00329 END IF 00330 00331 ! Get the total number of atoms, contracted spherical Gaussian basis 00332 ! functions, and atomic kinds 00333 00334 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 00335 natom=natom,& 00336 nsgf=nsgf) 00337 00338 nkind = SIZE(atomic_kind_set) 00339 00340 ALLOCATE (first_sgf_atom(natom),STAT=stat) 00341 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00342 first_sgf_atom(:) = 0 00343 00344 CALL get_particle_set(particle_set=particle_set,& 00345 first_sgf=first_sgf_atom,& 00346 error=error) 00347 00348 IF (PRESENT(matrix_h).OR.PRESENT(matrix_w)) THEN 00349 just_energy = .FALSE. 00350 ELSE 00351 just_energy = .TRUE. 00352 END IF 00353 00354 ! Retrieve S^(1/2) from the SCF environment 00355 00356 fm_s_half => scf_env%s_half 00357 CPPrecondition(ASSOCIATED(fm_s_half),cp_failure_level,routineP,error,failure) 00358 00359 ! Try to retrieve (full) work matrices from the SCF environment and reuse 00360 ! them, if available 00361 00362 fm_work1_local_alloc = .FALSE. 00363 fm_work2_local_alloc = .FALSE. 00364 00365 IF (ASSOCIATED(scf_env%scf_work1)) THEN 00366 IF (ASSOCIATED(scf_env%scf_work1(1)%matrix)) THEN 00367 fm_work1 => scf_env%scf_work1(1)%matrix 00368 END IF 00369 END IF 00370 00371 fm_work2 => scf_env%scf_work2 00372 00373 IF ((.NOT.ASSOCIATED(fm_work1)).OR.& 00374 (.NOT.ASSOCIATED(fm_work2))) THEN 00375 CALL cp_fm_struct_create(fmstruct=fmstruct,& 00376 para_env=qs_env%para_env,& 00377 context=qs_env%blacs_env,& 00378 nrow_global=nsgf,& 00379 ncol_global=nsgf,& 00380 error=error) 00381 IF (.NOT.ASSOCIATED(fm_work1)) THEN 00382 CALL cp_fm_create(matrix=fm_work1,& 00383 matrix_struct=fmstruct,& 00384 name="FULL WORK MATRIX 1",& 00385 error=error) 00386 fm_work1_local_alloc = .TRUE. 00387 END IF 00388 IF (.NOT.ASSOCIATED(fm_work2)) THEN 00389 CALL cp_fm_create(matrix=fm_work2,& 00390 matrix_struct=fmstruct,& 00391 name="FULL WORK MATRIX 2",& 00392 error=error) 00393 fm_work2_local_alloc = .TRUE. 00394 END IF 00395 CALL cp_fm_struct_release(fmstruct=fmstruct,error=error) 00396 END IF 00397 00398 ! Create local block diagonal matrices 00399 00400 ALLOCATE (sm_q,STAT=stat) 00401 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00402 CALL cp_dbcsr_init(sm_q,error=error) 00403 CALL cp_dbcsr_get_block_diag(sm_s,sm_q,error=error) 00404 00405 ALLOCATE (sm_v,STAT=stat) 00406 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00407 CALL cp_dbcsr_init(sm_v,error=error) 00408 CALL cp_dbcsr_get_block_diag(sm_s,sm_v,error=error) 00409 00410 ! Loop over all spins 00411 00412 DO ispin=1,nspin 00413 00414 IF (PRESENT(matrix_h)) THEN 00415 ! Hamiltonian matrix for spin ispin in sparse format 00416 sm_h => matrix_h(ispin)%matrix 00417 ELSE 00418 NULLIFY (sm_h) 00419 END IF 00420 00421 IF (PRESENT(matrix_w)) THEN 00422 ! Energy weighted density matrix for spin ispin in sparse format 00423 sm_w => matrix_w(ispin)%matrix 00424 ELSE 00425 NULLIFY (sm_w) 00426 END IF 00427 00428 sm_p => matrix_p(ispin)%matrix ! Density matrix for spin ispin in sparse format 00429 00430 CALL cp_dbcsr_set(sm_q,0.0_dp,error=error) 00431 CALL cp_dbcsr_set(sm_v,0.0_dp,error=error) 00432 00433 ! Calculate S^(1/2)*P*S^(1/2) as a full matrix (Lowdin) 00434 00435 CALL cp_dbcsr_sm_fm_multiply(sm_p,fm_s_half,fm_work1,nsgf,error=error) 00436 CALL cp_fm_gemm(transa="N",& 00437 transb="N",& 00438 m=nsgf,& 00439 n=nsgf,& 00440 k=nsgf,& 00441 alpha=1.0_dp,& 00442 matrix_a=fm_s_half,& 00443 matrix_b=fm_work1,& 00444 beta=0.0_dp,& 00445 matrix_c=fm_work2,& 00446 error=error) 00447 IF (debug) THEN 00448 CALL cp_dbcsr_write_sparse_matrix(sm_p,4,6,qs_env,qs_env%para_env,& 00449 output_unit=output_unit,& 00450 error=error) 00451 CALL write_fm_with_basis_info(fm_s_half,4,6,qs_env,qs_env%para_env,& 00452 output_unit=output_unit,& 00453 error=error) 00454 CALL write_fm_with_basis_info(fm_work2,4,6,qs_env,qs_env%para_env,& 00455 output_unit=output_unit,& 00456 error=error) 00457 END IF ! debug 00458 00459 ! Copy occupation matrix to sparse matrix format, finally we are only 00460 ! interested in the diagonal (atomic) blocks, i.e. the previous full 00461 ! matrix product is not the most efficient choice, anyway. 00462 00463 CALL copy_fm_to_dbcsr(fm_work2,sm_q,keep_sparsity=.TRUE.,error=error) 00464 00465 ! E[DFT+U] = E[DFT] + E[U] 00466 ! = E[DFT] + (U - J)*(Tr(q) - Tr(q*q))/2 00467 00468 ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U] 00469 ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j) 00470 ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j)) 00471 00472 ! Loop over all atomic kinds 00473 00474 DO ikind=1,nkind 00475 00476 atomic_kind => atomic_kind_set(ikind) 00477 00478 ! Load the required atomic kind data 00479 00480 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00481 atom_list=atom_list,& 00482 dft_plus_u_atom=dft_plus_u_atom,& 00483 l_of_dft_plus_u=lu,& 00484 name=atomic_kind_name,& 00485 natom=natom_of_kind,& 00486 nsgf=nsgf_kind,& 00487 orb_basis_set=orb_basis_set,& 00488 u_minus_j=u_minus_j,& 00489 u_minus_j_target=u_minus_j_target,& 00490 u_ramping=u_ramping,& 00491 eps_u_ramping=eps_u_ramping,& 00492 orbitals=orbitals,& 00493 eps_scf=eps_scf,& 00494 max_scf=max_scf,& 00495 smear=smear) 00496 00497 ! Check, if the atoms of this atomic kind need a DFT+U correction 00498 00499 IF (.NOT.ASSOCIATED(orb_basis_set)) CYCLE 00500 IF (.NOT.dft_plus_u_atom) CYCLE 00501 IF (lu < 0) CYCLE 00502 00503 ! Apply U ramping if requested 00504 00505 IF ((ispin == 1).AND.(u_ramping > 0.0_dp)) THEN 00506 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN 00507 u_minus_j = MIN(u_minus_j + u_ramping,u_minus_j_target) 00508 CALL set_atomic_kind(atomic_kind=atomic_kind,u_minus_j=u_minus_j) 00509 END IF 00510 IF (should_output.AND.(output_unit > 0)) THEN 00511 WRITE (UNIT=output_unit,FMT="(T3,A,3X,A,F0.3,A)")& 00512 "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)),& 00513 "U(eff) = ",u_minus_j*evolt," eV" 00514 END IF 00515 END IF 00516 00517 IF (u_minus_j == 0.0_dp) CYCLE 00518 00519 ! Load the required Gaussian basis set data 00520 00521 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 00522 first_sgf=first_sgf,& 00523 l=l,& 00524 last_sgf=last_sgf,& 00525 nset=nset,& 00526 nshell=nshell) 00527 00528 ! Count the relevant shell blocks of this atomic kind 00529 00530 nsb = 0 00531 DO iset=1,nset 00532 DO ishell=1,nshell(iset) 00533 IF (l(ishell,iset) == lu) nsb = nsb + 1 00534 END DO 00535 END DO 00536 00537 nsbsize = (2*lu + 1) 00538 n = nsb*nsbsize 00539 00540 ALLOCATE (q_matrix(n,n),STAT=stat) 00541 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00542 q_matrix(:,:) = 0.0_dp 00543 00544 ! Print headline if requested 00545 00546 IF (should_output.AND.(print_level > low_print_level)) THEN 00547 IF (output_unit > 0) THEN 00548 ALLOCATE (symbol(nsbsize),STAT=stat) 00549 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00550 DO m=-lu,lu 00551 symbol(lu+m+1) = sgf_symbol(0,lu,m) 00552 END DO 00553 IF (nspin > 1) THEN 00554 WRITE (UNIT=spin_info,FMT="(A8,I2)") " of spin",ispin 00555 ELSE 00556 spin_info = "" 00557 END IF 00558 WRITE (UNIT=output_unit,FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))")& 00559 "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ",ikind,& 00560 ": "//TRIM(atomic_kind_name),& 00561 "Atom Shell ",(ADJUSTR(symbol(i)),i=1,nsbsize)," Trace" 00562 DEALLOCATE (symbol,STAT=stat) 00563 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00564 END IF 00565 END IF 00566 00567 ! Loop over all atoms of the current atomic kind 00568 00569 DO iatom=1,natom_of_kind 00570 00571 atom_a = atom_list(iatom) 00572 00573 q_matrix(:,:) = 0.0_dp 00574 00575 ! Get diagonal block 00576 00577 CALL cp_dbcsr_get_block_p(matrix=sm_q,& 00578 row=atom_a,& 00579 col=atom_a,& 00580 block=q_block,& 00581 found=found) 00582 00583 IF (ASSOCIATED(q_block)) THEN 00584 00585 ! Calculate energy contribution to E(U) 00586 00587 i = 0 00588 DO iset=1,nset 00589 DO ishell=1,nshell(iset) 00590 IF (l(ishell,iset) /= lu) CYCLE 00591 DO isgf=first_sgf(ishell,iset),last_sgf(ishell,iset) 00592 i = i + 1 00593 j = 0 00594 DO jset=1,nset 00595 DO jshell=1,nshell(jset) 00596 IF (l(jshell,jset) /= lu) CYCLE 00597 DO jsgf=first_sgf(jshell,jset),last_sgf(jshell,jset) 00598 j = j + 1 00599 q_matrix(i,j) = q_block(isgf,jsgf) 00600 END DO ! next contracted spherical Gaussian function "jsgf" 00601 END DO ! next shell "jshell" 00602 END DO ! next shell set "jset" 00603 END DO ! next contracted spherical Gaussian function "isgf" 00604 END DO ! next shell "ishell" 00605 END DO ! next shell set "iset" 00606 00607 ! Perform the requested manipulations of the (initial) orbital occupations 00608 00609 IF (ASSOCIATED(orbitals)) THEN 00610 IF ((qs_env%scf_env%iter_delta >= eps_scf).OR.& 00611 ((qs_env%scf_env%outer_scf%iter_count == 0).AND.& 00612 (qs_env%scf_env%iter_count <= max_scf))) THEN 00613 ALLOCATE (orb_occ(nsbsize),STAT=stat) 00614 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00615 ALLOCATE (q_eigval(n),STAT=stat) 00616 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00617 q_eigval(:) = 0.0_dp 00618 ALLOCATE (q_eigvec(n,n),STAT=stat) 00619 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00620 q_eigvec(:,:) = 0.0_dp 00621 norb = SIZE(orbitals) 00622 CALL jacobi(q_matrix,q_eigval,q_eigvec) 00623 q_matrix(:,:) = 0.0_dp 00624 DO isb=1,nsb 00625 trq = 0.0_dp 00626 DO i=(isb-1)*nsbsize+1,isb*nsbsize 00627 trq = trq + q_eigval(i) 00628 END DO 00629 IF (smear) THEN 00630 occ = trq/REAL(norb,KIND=dp) 00631 ELSE 00632 occ = 1.0_dp/fspin 00633 END IF 00634 orb_occ(:) = .FALSE. 00635 iloc = MAXLOC(q_eigvec(:,isb*nsbsize)) 00636 jsb = INT((iloc(1)-1)/nsbsize) + 1 00637 i = 0 00638 i0 = (jsb-1)*nsbsize + 1 00639 iorb = -1000 00640 DO j=i0,jsb*nsbsize 00641 i = i + 1 00642 IF (i > norb) THEN 00643 DO m=-lu,lu 00644 IF (.NOT.orb_occ(lu+m+1)) THEN 00645 iorb = i0 + lu + m 00646 orb_occ(lu+m+1) = .TRUE. 00647 END IF 00648 END DO 00649 ELSE 00650 iorb = i0 + lu + orbitals(i) 00651 orb_occ(lu+orbitals(i)+1) = .TRUE. 00652 END IF 00653 CPPostcondition((iorb /= -1000),cp_failure_level,routineP,error,failure) 00654 iloc = MAXLOC(q_eigvec(iorb,:)) 00655 q_eigval(iloc(1)) = MIN(occ,trq) 00656 q_matrix(:,iloc(1)) = q_eigval(iloc(1))*q_eigvec(:,iloc(1)) ! backtransform left 00657 trq = trq - q_eigval(iloc(1)) 00658 END DO 00659 END DO 00660 q_matrix = MATMUL(q_matrix,TRANSPOSE(q_eigvec)) ! backtransform right 00661 DEALLOCATE (orb_occ,STAT=stat) 00662 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00663 DEALLOCATE (q_eigval,STAT=stat) 00664 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00665 DEALLOCATE (q_eigvec,STAT=stat) 00666 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00667 END IF 00668 END IF ! orbitals associated 00669 00670 trq = 0.0_dp 00671 trq2 = 0.0_dp 00672 00673 DO i=1,n 00674 trq = trq + q_matrix(i,i) 00675 DO j=1,n 00676 trq2 = trq2 + q_matrix(i,j)*q_matrix(j,i) 00677 END DO 00678 END DO 00679 00680 trq = fspin*trq 00681 trq2 = fspin*fspin*trq2 00682 00683 ! Calculate energy contribution to E(U) 00684 00685 energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin 00686 00687 ! Calculate potential V(U) = dE(U)/dq 00688 00689 IF (.NOT.just_energy) THEN 00690 00691 CALL cp_dbcsr_get_block_p(matrix=sm_v,& 00692 row=atom_a,& 00693 col=atom_a,& 00694 block=v_block,& 00695 found=found) 00696 CPPostcondition(ASSOCIATED(v_block),cp_failure_level,routineP,error,failure) 00697 00698 i = 0 00699 DO iset=1,nset 00700 DO ishell=1,nshell(iset) 00701 IF (l(ishell,iset) /= lu) CYCLE 00702 DO isgf=first_sgf(ishell,iset),last_sgf(ishell,iset) 00703 i = i + 1 00704 j = 0 00705 DO jset=1,nset 00706 DO jshell=1,nshell(jset) 00707 IF (l(jshell,jset) /= lu) CYCLE 00708 DO jsgf=first_sgf(jshell,jset),last_sgf(jshell,jset) 00709 j = j + 1 00710 IF (isgf == jsgf) THEN 00711 v_block(isgf,isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i,i)) 00712 ELSE 00713 v_block(isgf,jsgf) = -u_minus_j*fspin*q_matrix(j,i) 00714 END IF 00715 END DO ! next contracted spherical Gaussian function "jsgf" 00716 END DO ! next shell "jshell" 00717 END DO ! next shell set "jset" 00718 END DO ! next contracted spherical Gaussian function "isgf" 00719 END DO ! next shell "ishell" 00720 END DO ! next shell set "iset" 00721 00722 END IF ! not just energy 00723 00724 END IF ! q_block associated 00725 00726 ! Consider print requests 00727 00728 IF (should_output.AND.(print_level > low_print_level)) THEN 00729 CALL mp_sum(q_matrix,qs_env%para_env%group) 00730 IF (output_unit > 0) THEN 00731 ALLOCATE (q_work(nsb,nsbsize),STAT=stat) 00732 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00733 q_work(:,:) = 0.0_dp 00734 DO isb=1,nsb 00735 j = 0 00736 DO i=(isb-1)*nsbsize+1,isb*nsbsize 00737 j = j + 1 00738 q_work(isb,j) = q_matrix(i,i) 00739 END DO 00740 END DO 00741 DO isb=1,nsb 00742 WRITE (UNIT=output_unit,FMT="(T3,I6,2X,I6,2X,10F8.3)")& 00743 atom_a,isb,q_work(isb,:),SUM(q_work(isb,:)) 00744 END DO 00745 WRITE (UNIT=output_unit,FMT="(T12,A,2X,10F8.3)")& 00746 "Total",(SUM(q_work(:,i)),i=1,nsbsize),SUM(q_work) 00747 WRITE (UNIT=output_unit,FMT="(A)") "" 00748 DEALLOCATE (q_work,STAT=stat) 00749 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00750 IF (debug) THEN 00751 ! Print the DFT+U occupation matrix 00752 WRITE (UNIT=output_unit,FMT="(T9,70I10)") (i,i=1,n) 00753 DO i=1,n 00754 WRITE (UNIT=output_unit,FMT="(T3,I6,70F10.6)") i,q_matrix(i,:) 00755 END DO 00756 ! Print the eigenvalues and eigenvectors of the occupation matrix 00757 ALLOCATE (q_eigval(n),STAT=stat) 00758 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00759 q_eigval(:) = 0.0_dp 00760 ALLOCATE (q_eigvec(n,n),STAT=stat) 00761 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00762 q_eigvec(:,:) = 0.0_dp 00763 CALL jacobi(q_matrix,q_eigval,q_eigvec) 00764 WRITE (UNIT=output_unit,FMT="(/,T9,70I10)") (i,i=1,n) 00765 WRITE (UNIT=output_unit,FMT="(T9,71F10.6)") (q_eigval(i),i=1,n),& 00766 SUM(q_eigval(1:n)) 00767 DO i=1,n 00768 WRITE (UNIT=output_unit,FMT="(T3,I6,70F10.6)") i,q_eigvec(i,:) 00769 END DO 00770 DEALLOCATE (q_eigval,STAT=stat) 00771 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00772 DEALLOCATE (q_eigvec,STAT=stat) 00773 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00774 END IF ! debug 00775 END IF 00776 IF (debug) THEN 00777 ! Print the full atomic occupation matrix block 00778 ALLOCATE (q_work(nsgf_kind,nsgf_kind),STAT=stat) 00779 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00780 q_work(:,:) = 0.0_dp 00781 IF (ASSOCIATED(q_block)) q_work(:,:) = q_block(:,:) 00782 CALL mp_sum(q_work,qs_env%para_env%group) 00783 IF (output_unit > 0) THEN 00784 norb = SIZE(q_work,1) 00785 WRITE (UNIT=output_unit,FMT="(/,T9,200I10)") (i,i=1,norb) 00786 DO i=1,norb 00787 WRITE (UNIT=output_unit,FMT="(T3,I6,200F10.6)") i,q_work(i,:) 00788 END DO 00789 ALLOCATE (q_eigval(norb),STAT=stat) 00790 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00791 q_eigval(:) = 0.0_dp 00792 ALLOCATE (q_eigvec(norb,norb),STAT=stat) 00793 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00794 q_eigvec(:,:) = 0.0_dp 00795 CALL jacobi(q_work,q_eigval,q_eigvec) 00796 WRITE (UNIT=output_unit,FMT="(/,T9,200I10)") (i,i=1,norb) 00797 WRITE (UNIT=output_unit,FMT="(T9,201F10.6)") (q_eigval(i),i=1,norb),& 00798 SUM(q_eigval(1:norb)) 00799 DO i=1,norb 00800 WRITE (UNIT=output_unit,FMT="(T3,I6,200F10.6)") i,q_eigvec(i,:) 00801 END DO 00802 DEALLOCATE (q_eigval,STAT=stat) 00803 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00804 DEALLOCATE (q_eigvec,STAT=stat) 00805 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00806 END IF 00807 DEALLOCATE (q_work,STAT=stat) 00808 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00809 END IF ! debug 00810 END IF ! should output 00811 00812 END DO ! next atom "iatom" of atomic kind "ikind" 00813 00814 IF (ALLOCATED(q_matrix)) THEN 00815 DEALLOCATE (q_matrix,STAT=stat) 00816 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 00817 END IF 00818 END DO ! next atomic kind "ikind" 00819 00820 ! Add V(i,j)[U] to V(i,j)[DFT] 00821 00822 IF (ASSOCIATED(sm_h)) THEN 00823 00824 CALL cp_dbcsr_sm_fm_multiply(sm_v,fm_s_half,fm_work1,nsgf,error=error) 00825 CALL cp_fm_transpose(fm_work1,fm_work2,error=error) 00826 CALL cp_dbcsr_plus_fm_fm_t(sm_h,fm_s_half,fm_work2,nsgf,error=error) 00827 00828 END IF ! An update of the Hamiltonian matrix is requested 00829 00830 ! Calculate the contribution (non-Pulay part) to the derivatives 00831 ! w.r.t. the nuclear positions 00832 00833 IF (ASSOCIATED(sm_w)) THEN 00834 00835 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 00836 "Forces are not yet fully implemented for DFT+U method LOWDIN",& 00837 only_ionode=.TRUE.) 00838 00839 END IF ! W matrix update requested 00840 00841 END DO ! next spin "ispin" 00842 00843 ! Collect the energy contributions from all processes 00844 00845 CALL mp_sum(energy%dft_plus_u,qs_env%para_env%group) 00846 00847 CALL cp_assert(energy%dft_plus_u >= 0.0_dp,cp_warning_level,& 00848 cp_assertion_failed,routineP,& 00849 "DFT+U energy contibution is negative possibly due "//& 00850 "to unphysical Lowdin charges. Check your input, "//& 00851 "if this warning persists or try a different method!",& 00852 only_ionode=.TRUE.) 00853 00854 ! Release (local) full matrices 00855 00856 NULLIFY (fm_s_half) 00857 IF (fm_work1_local_alloc) THEN 00858 CALL cp_fm_release(matrix=fm_work1,error=error) 00859 END IF 00860 IF (fm_work2_local_alloc) THEN 00861 CALL cp_fm_release(matrix=fm_work2,error=error) 00862 END IF 00863 00864 ! Release (local) sparse matrices 00865 00866 CALL cp_dbcsr_deallocate_matrix(sm_q,error=error) 00867 CALL cp_dbcsr_deallocate_matrix(sm_v,error=error) 00868 00869 CALL timestop(handle) 00870 00871 END SUBROUTINE lowdin 00872 00873 ! ***************************************************************************** 00900 SUBROUTINE mulliken(qs_env,orthonormal_basis,matrix_h,should_output,& 00901 output_unit,print_level,error) 00902 00903 TYPE(qs_environment_type), POINTER :: qs_env 00904 LOGICAL, INTENT(IN) :: orthonormal_basis 00905 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00906 OPTIONAL, POINTER :: matrix_h 00907 LOGICAL, INTENT(IN) :: should_output 00908 INTEGER, INTENT(IN) :: output_unit, print_level 00909 TYPE(cp_error_type), INTENT(INOUT) :: error 00910 00911 CHARACTER(LEN=*), PARAMETER :: routineN = 'mulliken', 00912 routineP = moduleN//':'//routineN 00913 00914 CHARACTER(LEN=10) :: spin_info 00915 CHARACTER(LEN=6), ALLOCATABLE, 00916 DIMENSION(:) :: symbol 00917 CHARACTER(LEN=default_string_length) :: atomic_kind_name 00918 INTEGER :: atom_a, handle, i, i0, iatom, ikind, iorb, isb, iset, isgf, 00919 ishell, ispin, j, jsb, jset, jsgf, jshell, lu, m, max_scf, n, natom, 00920 natom_of_kind, nkind, norb, nsb, nsbsize, nset, nsgf_kind, nspin, stat 00921 INTEGER, DIMENSION(1) :: iloc 00922 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell, orbitals 00923 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf 00924 LOGICAL :: debug, dft_plus_u_atom, 00925 failure, found, just_energy, 00926 smear 00927 LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_plus_u_kind, orb_occ 00928 REAL(KIND=dp) :: eps_scf, eps_u_ramping, 00929 fspin, occ, trq, trq2, 00930 u_minus_j, u_minus_j_target, 00931 u_ramping 00932 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: q_eigval 00933 REAL(KIND=dp), ALLOCATABLE, 00934 DIMENSION(:, :) :: q_eigvec, q_matrix, q_work 00935 REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, q_block, 00936 s_block, v_block 00937 TYPE(atomic_kind_type), DIMENSION(:), 00938 POINTER :: atomic_kind_set 00939 TYPE(atomic_kind_type), POINTER :: atomic_kind, kind_a 00940 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00941 POINTER :: matrix_p, matrix_s 00942 TYPE(cp_dbcsr_type), POINTER :: sm_h, sm_p, sm_q, sm_s, sm_v 00943 TYPE(dft_control_type), POINTER :: dft_control 00944 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 00945 TYPE(particle_type), DIMENSION(:), 00946 POINTER :: particle_set 00947 TYPE(qs_energy_type), POINTER :: energy 00948 TYPE(qs_rho_type), POINTER :: rho 00949 00950 CALL timeset(routineN,handle) 00951 00952 debug = .FALSE. ! Set to .TRUE. to print debug information 00953 00954 NULLIFY (atom_list) 00955 NULLIFY (atomic_kind) 00956 NULLIFY (atomic_kind_set) 00957 NULLIFY (dft_control) 00958 NULLIFY (energy) 00959 NULLIFY (first_sgf) 00960 NULLIFY (h_block) 00961 NULLIFY (matrix_p) 00962 NULLIFY (matrix_s) 00963 NULLIFY (l) 00964 NULLIFY (last_sgf) 00965 NULLIFY (nshell) 00966 NULLIFY (orb_basis_set) 00967 NULLIFY (p_block) 00968 NULLIFY (particle_set) 00969 NULLIFY (q_block) 00970 NULLIFY (rho) 00971 NULLIFY (s_block) 00972 NULLIFY (orbitals) 00973 NULLIFY (sm_h) 00974 NULLIFY (sm_p) 00975 NULLIFY (sm_q) 00976 NULLIFY (sm_s) 00977 NULLIFY (sm_v) 00978 NULLIFY (v_block) 00979 00980 failure = .FALSE. 00981 00982 smear = .FALSE. 00983 max_scf = -1 00984 eps_scf = 1.0E30_dp 00985 00986 CALL get_qs_env(qs_env=qs_env,& 00987 atomic_kind_set=atomic_kind_set,& 00988 dft_control=dft_control,& 00989 energy=energy,& 00990 particle_set=particle_set,& 00991 rho=rho,& 00992 error=error) 00993 00994 CPPrecondition(ASSOCIATED(atomic_kind_set),cp_failure_level,routineP,error,failure) 00995 CPPrecondition(ASSOCIATED(dft_control),cp_failure_level,routineP,error,failure) 00996 CPPrecondition(ASSOCIATED(energy),cp_failure_level,routineP,error,failure) 00997 CPPrecondition(ASSOCIATED(particle_set),cp_failure_level,routineP,error,failure) 00998 CPPrecondition(ASSOCIATED(rho),cp_failure_level,routineP,error,failure) 00999 01000 IF (orthonormal_basis) THEN 01001 NULLIFY (sm_s) 01002 ELSE 01003 ! Get overlap matrix in sparse format 01004 CALL get_qs_env(qs_env=qs_env,& 01005 matrix_s=matrix_s,& 01006 error=error) 01007 CPPrecondition(ASSOCIATED(matrix_s),cp_failure_level,routineP,error,failure) 01008 sm_s => matrix_s(1)%matrix 01009 END IF 01010 01011 ! Get density matrices in sparse format 01012 01013 matrix_p => rho%rho_ao 01014 01015 energy%dft_plus_u = 0.0_dp 01016 01017 nspin = dft_control%nspins 01018 01019 IF (nspin == 2) THEN 01020 fspin = 1.0_dp 01021 ELSE 01022 fspin = 0.5_dp 01023 END IF 01024 01025 ! Get the total number of atoms, contracted spherical Gaussian basis 01026 ! functions, and atomic kinds 01027 01028 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 01029 natom=natom) 01030 01031 nkind = SIZE(atomic_kind_set) 01032 01033 ALLOCATE (is_plus_u_kind(nkind),STAT=stat) 01034 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01035 is_plus_u_kind(:) = .FALSE. 01036 01037 IF (PRESENT(matrix_h)) THEN 01038 just_energy = .FALSE. 01039 ELSE 01040 just_energy = .TRUE. 01041 END IF 01042 01043 ! Loop over all spins 01044 01045 DO ispin=1,nspin 01046 01047 IF (PRESENT(matrix_h)) THEN 01048 ! Hamiltonian matrix for spin ispin in sparse format 01049 sm_h => matrix_h(ispin)%matrix 01050 ELSE 01051 NULLIFY (sm_h) 01052 END IF 01053 01054 ! Get density matrix for spin ispin in sparse format 01055 01056 sm_p => matrix_p(ispin)%matrix 01057 01058 IF (.NOT.ASSOCIATED(sm_q)) THEN 01059 ALLOCATE (sm_q,STAT=stat) 01060 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01061 CALL cp_dbcsr_init(sm_q,error=error) 01062 CALL cp_dbcsr_get_block_diag(sm_p,sm_q,error=error) 01063 END IF 01064 CALL cp_dbcsr_set(sm_q,0.0_dp,error=error) 01065 01066 IF (.NOT.ASSOCIATED(sm_v)) THEN 01067 ALLOCATE (sm_v,STAT=stat) 01068 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01069 CALL cp_dbcsr_init(sm_v,error=error) 01070 CALL cp_dbcsr_get_block_diag(sm_p,sm_v,error=error) 01071 END IF 01072 CALL cp_dbcsr_set(sm_v,0.0_dp,error=error) 01073 01074 DO iatom=1,natom 01075 01076 CALL cp_dbcsr_get_block_p(matrix=sm_p,& 01077 row=iatom,& 01078 col=iatom,& 01079 block=p_block,& 01080 found=found) 01081 01082 IF (.NOT.ASSOCIATED(p_block)) CYCLE 01083 01084 CALL cp_dbcsr_get_block_p(matrix=sm_q,& 01085 row=iatom,& 01086 col=iatom,& 01087 block=q_block,& 01088 found=found) 01089 CPPostcondition(ASSOCIATED(q_block),cp_failure_level,routineP,error,failure) 01090 01091 IF (orthonormal_basis) THEN 01092 ! S is the unit matrix 01093 DO isgf=1,SIZE(q_block,1) 01094 q_block(isgf,isgf) = p_block(isgf,isgf) 01095 END DO 01096 ELSE 01097 CALL cp_dbcsr_get_block_p(matrix=sm_s,& 01098 row=iatom,& 01099 col=iatom,& 01100 block=s_block,& 01101 found=found) 01102 CPPostcondition(ASSOCIATED(s_block),cp_failure_level,routineP,error,failure) 01103 ! Exploit that P and S are symmetric 01104 DO jsgf=1,SIZE(p_block,2) 01105 DO isgf=1,SIZE(p_block,1) 01106 q_block(isgf,jsgf) = p_block(isgf,jsgf)*s_block(isgf,jsgf) 01107 END DO 01108 END DO 01109 END IF ! orthonormal basis set 01110 01111 END DO ! next atom "iatom" 01112 01113 ! E[DFT+U] = E[DFT] + E[U] 01114 ! = E[DFT] + (U - J)*(Tr(q) - Tr(q*q))/2 01115 01116 ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U] 01117 ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j) 01118 ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j)) 01119 01120 ! Loop over all atomic kinds 01121 01122 DO ikind=1,nkind 01123 01124 atomic_kind => atomic_kind_set(ikind) 01125 01126 ! Load the required atomic kind data 01127 01128 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01129 atom_list=atom_list,& 01130 dft_plus_u_atom=dft_plus_u_atom,& 01131 l_of_dft_plus_u=lu,& 01132 name=atomic_kind_name,& 01133 natom=natom_of_kind,& 01134 nsgf=nsgf_kind,& 01135 orb_basis_set=orb_basis_set,& 01136 u_minus_j=u_minus_j,& 01137 u_minus_j_target=u_minus_j_target,& 01138 u_ramping=u_ramping,& 01139 eps_u_ramping=eps_u_ramping,& 01140 orbitals=orbitals,& 01141 eps_scf=eps_scf,& 01142 max_scf=max_scf,& 01143 smear=smear) 01144 01145 ! Check, if the atoms of this atomic kind need a DFT+U correction 01146 01147 IF (.NOT.ASSOCIATED(orb_basis_set)) CYCLE 01148 IF (.NOT.dft_plus_u_atom) CYCLE 01149 IF (lu < 0) CYCLE 01150 01151 ! Apply U ramping if requested 01152 01153 IF ((ispin == 1).AND.(u_ramping > 0.0_dp)) THEN 01154 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN 01155 u_minus_j = MIN(u_minus_j + u_ramping,u_minus_j_target) 01156 CALL set_atomic_kind(atomic_kind=atomic_kind,u_minus_j=u_minus_j) 01157 END IF 01158 IF (should_output.AND.(output_unit > 0)) THEN 01159 WRITE (UNIT=output_unit,FMT="(T3,A,3X,A,F0.3,A)")& 01160 "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)),& 01161 "U(eff) = ",u_minus_j*evolt," eV" 01162 END IF 01163 END IF 01164 01165 IF (u_minus_j == 0.0_dp) CYCLE 01166 01167 is_plus_u_kind(ikind) = .TRUE. 01168 01169 ! Load the required Gaussian basis set data 01170 01171 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01172 first_sgf=first_sgf,& 01173 l=l,& 01174 last_sgf=last_sgf,& 01175 nset=nset,& 01176 nshell=nshell) 01177 01178 ! Count the relevant shell blocks of this atomic kind 01179 01180 nsb = 0 01181 DO iset=1,nset 01182 DO ishell=1,nshell(iset) 01183 IF (l(ishell,iset) == lu) nsb = nsb + 1 01184 END DO 01185 END DO 01186 01187 nsbsize = (2*lu + 1) 01188 n = nsb*nsbsize 01189 01190 ALLOCATE (q_matrix(n,n),STAT=stat) 01191 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01192 q_matrix(:,:) = 0.0_dp 01193 01194 ! Print headline if requested 01195 01196 IF (should_output.AND.(print_level > low_print_level)) THEN 01197 IF (output_unit > 0) THEN 01198 ALLOCATE (symbol(nsbsize),STAT=stat) 01199 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01200 DO m=-lu,lu 01201 symbol(lu+m+1) = sgf_symbol(0,lu,m) 01202 END DO 01203 IF (nspin > 1) THEN 01204 WRITE (UNIT=spin_info,FMT="(A8,I2)") " of spin",ispin 01205 ELSE 01206 spin_info = "" 01207 END IF 01208 WRITE (UNIT=output_unit,FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))")& 01209 "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ",ikind,& 01210 ": "//TRIM(atomic_kind_name),& 01211 "Atom Shell ",(ADJUSTR(symbol(i)),i=1,nsbsize)," Trace" 01212 DEALLOCATE (symbol,STAT=stat) 01213 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01214 END IF 01215 END IF 01216 01217 ! Loop over all atoms of the current atomic kind 01218 01219 DO iatom=1,natom_of_kind 01220 01221 atom_a = atom_list(iatom) 01222 01223 q_matrix(:,:) = 0.0_dp 01224 01225 ! Get diagonal block 01226 01227 CALL cp_dbcsr_get_block_p(matrix=sm_q,& 01228 row=atom_a,& 01229 col=atom_a,& 01230 block=q_block,& 01231 found=found) 01232 01233 ! Calculate energy contribution to E(U) 01234 01235 IF (ASSOCIATED(q_block)) THEN 01236 01237 i = 0 01238 DO iset=1,nset 01239 DO ishell=1,nshell(iset) 01240 IF (l(ishell,iset) /= lu) CYCLE 01241 DO isgf=first_sgf(ishell,iset),last_sgf(ishell,iset) 01242 i = i + 1 01243 j = 0 01244 DO jset=1,nset 01245 DO jshell=1,nshell(jset) 01246 IF (l(jshell,jset) /= lu) CYCLE 01247 DO jsgf=first_sgf(jshell,jset),last_sgf(jshell,jset) 01248 j = j + 1 01249 q_matrix(i,j) = q_block(isgf,jsgf) 01250 END DO ! next contracted spherical Gaussian function "jsgf" 01251 END DO ! next shell "jshell" 01252 END DO ! next shell set "jset" 01253 END DO ! next contracted spherical Gaussian function "isgf" 01254 END DO ! next shell "ishell" 01255 END DO ! next shell set "iset" 01256 01257 ! Perform the requested manipulations of the (initial) orbital occupations 01258 01259 IF (ASSOCIATED(orbitals)) THEN 01260 IF ((qs_env%scf_env%iter_delta >= eps_scf).OR.& 01261 ((qs_env%scf_env%outer_scf%iter_count == 0).AND.& 01262 (qs_env%scf_env%iter_count <= max_scf))) THEN 01263 ALLOCATE (orb_occ(nsbsize),STAT=stat) 01264 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01265 ALLOCATE (q_eigval(n),STAT=stat) 01266 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01267 q_eigval(:) = 0.0_dp 01268 ALLOCATE (q_eigvec(n,n),STAT=stat) 01269 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01270 q_eigvec(:,:) = 0.0_dp 01271 norb = SIZE(orbitals) 01272 CALL jacobi(q_matrix,q_eigval,q_eigvec) 01273 q_matrix(:,:) = 0.0_dp 01274 DO isb=1,nsb 01275 trq = 0.0_dp 01276 DO i=(isb-1)*nsbsize+1,isb*nsbsize 01277 trq = trq + q_eigval(i) 01278 END DO 01279 IF (smear) THEN 01280 occ = trq/REAL(norb,KIND=dp) 01281 ELSE 01282 occ = 1.0_dp/fspin 01283 END IF 01284 orb_occ(:) = .FALSE. 01285 iloc = MAXLOC(q_eigvec(:,isb*nsbsize)) 01286 jsb = INT((iloc(1)-1)/nsbsize) + 1 01287 i = 0 01288 i0 = (jsb-1)*nsbsize + 1 01289 iorb = -1000 01290 DO j=i0,jsb*nsbsize 01291 i = i + 1 01292 IF (i > norb) THEN 01293 DO m=-lu,lu 01294 IF (.NOT.orb_occ(lu+m+1)) THEN 01295 iorb = i0 + lu + m 01296 orb_occ(lu+m+1) = .TRUE. 01297 END IF 01298 END DO 01299 ELSE 01300 iorb = i0 + lu + orbitals(i) 01301 orb_occ(lu+orbitals(i)+1) = .TRUE. 01302 END IF 01303 CPPostcondition((iorb /= -1000),cp_failure_level,routineP,error,failure) 01304 iloc = MAXLOC(q_eigvec(iorb,:)) 01305 q_eigval(iloc(1)) = MIN(occ,trq) 01306 q_matrix(:,iloc(1)) = q_eigval(iloc(1))*q_eigvec(:,iloc(1)) ! backtransform left 01307 trq = trq - q_eigval(iloc(1)) 01308 END DO 01309 END DO 01310 q_matrix = MATMUL(q_matrix,TRANSPOSE(q_eigvec)) ! backtransform right 01311 DEALLOCATE (orb_occ,STAT=stat) 01312 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01313 DEALLOCATE (q_eigval,STAT=stat) 01314 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01315 DEALLOCATE (q_eigvec,STAT=stat) 01316 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01317 END IF 01318 END IF ! orbitals associated 01319 01320 trq = 0.0_dp 01321 trq2 = 0.0_dp 01322 01323 DO i=1,n 01324 trq = trq + q_matrix(i,i) 01325 DO j=1,n 01326 trq2 = trq2 + q_matrix(i,j)*q_matrix(j,i) 01327 END DO 01328 END DO 01329 01330 trq = fspin*trq 01331 trq2 = fspin*fspin*trq2 01332 01333 ! Calculate energy contribution to E(U) 01334 01335 energy%dft_plus_u = energy%dft_plus_u + 0.5_dp*u_minus_j*(trq - trq2)/fspin 01336 01337 ! Calculate potential V(U) = dE(U)/dq 01338 01339 IF (.NOT.just_energy) THEN 01340 01341 CALL cp_dbcsr_get_block_p(matrix=sm_v,& 01342 row=atom_a,& 01343 col=atom_a,& 01344 block=v_block,& 01345 found=found) 01346 CPPostcondition(ASSOCIATED(v_block),cp_failure_level,routineP,error,failure) 01347 01348 i = 0 01349 DO iset=1,nset 01350 DO ishell=1,nshell(iset) 01351 IF (l(ishell,iset) /= lu) CYCLE 01352 DO isgf=first_sgf(ishell,iset),last_sgf(ishell,iset) 01353 i = i + 1 01354 j = 0 01355 DO jset=1,nset 01356 DO jshell=1,nshell(jset) 01357 IF (l(jshell,jset) /= lu) CYCLE 01358 DO jsgf=first_sgf(jshell,jset),last_sgf(jshell,jset) 01359 j = j + 1 01360 IF (isgf == jsgf) THEN 01361 v_block(isgf,isgf) = u_minus_j*(0.5_dp - fspin*q_matrix(i,i)) 01362 ELSE 01363 v_block(isgf,jsgf) = -u_minus_j*fspin*q_matrix(j,i) 01364 END IF 01365 END DO ! next contracted spherical Gaussian function "jsgf" 01366 END DO ! next shell "jshell" 01367 END DO ! next shell set "jset" 01368 END DO ! next contracted spherical Gaussian function "isgf" 01369 END DO ! next shell "ishell" 01370 END DO ! next shell set "iset" 01371 01372 END IF ! not just energy 01373 01374 END IF ! q_block associated 01375 01376 ! Consider print requests 01377 01378 IF (should_output.AND.(print_level > low_print_level)) THEN 01379 CALL mp_sum(q_matrix,qs_env%para_env%group) 01380 IF (output_unit > 0) THEN 01381 ALLOCATE (q_work(nsb,nsbsize),STAT=stat) 01382 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01383 q_work(:,:) = 0.0_dp 01384 DO isb=1,nsb 01385 j = 0 01386 DO i=(isb-1)*nsbsize+1,isb*nsbsize 01387 j = j + 1 01388 q_work(isb,j) = q_matrix(i,i) 01389 END DO 01390 END DO 01391 DO isb=1,nsb 01392 WRITE (UNIT=output_unit,FMT="(T3,I6,2X,I6,2X,10F8.3)")& 01393 atom_a,isb,q_work(isb,:),SUM(q_work(isb,:)) 01394 END DO 01395 WRITE (UNIT=output_unit,FMT="(T12,A,2X,10F8.3)")& 01396 "Total",(SUM(q_work(:,i)),i=1,nsbsize),SUM(q_work) 01397 WRITE (UNIT=output_unit,FMT="(A)") "" 01398 DEALLOCATE (q_work,STAT=stat) 01399 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01400 IF (debug) THEN 01401 ! Print the DFT+U occupation matrix 01402 WRITE (UNIT=output_unit,FMT="(T9,70I10)") (i,i=1,n) 01403 DO i=1,n 01404 WRITE (UNIT=output_unit,FMT="(T3,I6,70F10.6)") i,q_matrix(i,:) 01405 END DO 01406 ! Print the eigenvalues and eigenvectors of the occupation matrix 01407 ALLOCATE (q_eigval(n),STAT=stat) 01408 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01409 q_eigval(:) = 0.0_dp 01410 ALLOCATE (q_eigvec(n,n),STAT=stat) 01411 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01412 q_eigvec(:,:) = 0.0_dp 01413 CALL jacobi(q_matrix,q_eigval,q_eigvec) 01414 WRITE (UNIT=output_unit,FMT="(/,T9,70I10)") (i,i=1,n) 01415 WRITE (UNIT=output_unit,FMT="(T9,71F10.6)") (q_eigval(i),i=1,n),& 01416 SUM(q_eigval(1:n)) 01417 DO i=1,n 01418 WRITE (UNIT=output_unit,FMT="(T3,I6,70F10.6)") i,q_eigvec(i,:) 01419 END DO 01420 DEALLOCATE (q_eigval,STAT=stat) 01421 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01422 DEALLOCATE (q_eigvec,STAT=stat) 01423 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01424 END IF ! debug 01425 END IF 01426 IF (debug) THEN 01427 ! Print the full atomic occupation matrix block 01428 ALLOCATE (q_work(nsgf_kind,nsgf_kind),STAT=stat) 01429 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01430 q_work(:,:) = 0.0_dp 01431 IF (ASSOCIATED(q_block)) q_work(:,:) = q_block(:,:) 01432 CALL mp_sum(q_work,qs_env%para_env%group) 01433 IF (output_unit > 0) THEN 01434 norb = SIZE(q_work,1) 01435 WRITE (UNIT=output_unit,FMT="(/,T9,200I10)") (i,i=1,norb) 01436 DO i=1,norb 01437 WRITE (UNIT=output_unit,FMT="(T3,I6,200F10.6)") i,q_work(i,:) 01438 END DO 01439 ALLOCATE (q_eigval(norb),STAT=stat) 01440 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01441 q_eigval(:) = 0.0_dp 01442 ALLOCATE (q_eigvec(norb,norb),STAT=stat) 01443 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01444 q_eigvec(:,:) = 0.0_dp 01445 CALL jacobi(q_work,q_eigval,q_eigvec) 01446 WRITE (UNIT=output_unit,FMT="(/,T9,200I10)") (i,i=1,norb) 01447 WRITE (UNIT=output_unit,FMT="(T9,201F10.6)") (q_eigval(i),i=1,norb),& 01448 SUM(q_eigval(1:norb)) 01449 DO i=1,norb 01450 WRITE (UNIT=output_unit,FMT="(T3,I6,200F10.6)") i,q_eigvec(i,:) 01451 END DO 01452 DEALLOCATE (q_eigval,STAT=stat) 01453 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01454 DEALLOCATE (q_eigvec,STAT=stat) 01455 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01456 END IF 01457 DEALLOCATE (q_work,STAT=stat) 01458 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01459 END IF ! debug 01460 END IF ! should output 01461 01462 END DO ! next atom "iatom" of atomic kind "ikind" 01463 01464 IF (ALLOCATED(q_matrix)) THEN 01465 DEALLOCATE (q_matrix,STAT=stat) 01466 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01467 END IF 01468 01469 END DO ! next atomic kind "ikind" 01470 01471 ! Add V(i,j)[U] to V(i,j)[DFT] 01472 01473 IF (ASSOCIATED(sm_h)) THEN 01474 01475 DO ikind=1,nkind 01476 01477 IF (.NOT.is_plus_u_kind(ikind)) CYCLE 01478 01479 kind_a => atomic_kind_set(ikind) 01480 01481 CALL get_atomic_kind(atomic_kind=kind_a,& 01482 atom_list=atom_list,& 01483 natom=natom_of_kind) 01484 01485 DO iatom=1,natom_of_kind 01486 01487 atom_a = atom_list(iatom) 01488 01489 CALL cp_dbcsr_get_block_p(matrix=sm_h,& 01490 row=atom_a,& 01491 col=atom_a,& 01492 block=h_block,& 01493 found=found) 01494 01495 IF (.NOT.ASSOCIATED(h_block)) CYCLE 01496 01497 CALL cp_dbcsr_get_block_p(matrix=sm_v,& 01498 row=atom_a,& 01499 col=atom_a,& 01500 block=v_block,& 01501 found=found) 01502 CPPostcondition(ASSOCIATED(v_block),cp_failure_level,routineP,error,failure) 01503 01504 IF (orthonormal_basis) THEN 01505 DO isgf=1,SIZE(h_block,1) 01506 h_block(isgf,isgf) = h_block(isgf,isgf) + v_block(isgf,isgf) 01507 END DO 01508 ELSE 01509 CALL cp_dbcsr_get_block_p(matrix=sm_s,& 01510 row=atom_a,& 01511 col=atom_a,& 01512 block=s_block,& 01513 found=found) 01514 CPPostcondition(ASSOCIATED(s_block),cp_failure_level,routineP,error,failure) 01515 DO jsgf=1,SIZE(h_block,2) 01516 DO isgf=1,SIZE(h_block,1) 01517 h_block(isgf,jsgf) = h_block(isgf,jsgf) + v_block(isgf,jsgf)*s_block(isgf,jsgf) 01518 END DO 01519 END DO 01520 END IF ! orthonormal basis set 01521 01522 END DO ! next atom "iatom" of atomic kind "ikind" 01523 01524 END DO ! Next atomic kind "ikind" 01525 01526 END IF ! An update of the Hamiltonian matrix is requested 01527 01528 END DO ! next spin "ispin" 01529 01530 ! Collect the energy contributions from all processes 01531 01532 CALL mp_sum(energy%dft_plus_u,qs_env%para_env%group) 01533 01534 CALL cp_assert(energy%dft_plus_u >= 0.0_dp,cp_warning_level,& 01535 cp_assertion_failed,routineP,& 01536 "DFT+U energy contibution is negative possibly due "//& 01537 "to unphysical Mulliken charges. Check your input, "//& 01538 "if this warning persists or try a different method!",& 01539 only_ionode=.TRUE.) 01540 01541 CALL cp_dbcsr_deallocate_matrix(sm_q,error=error) 01542 CALL cp_dbcsr_deallocate_matrix(sm_v,error=error) 01543 01544 CALL timestop(handle) 01545 01546 END SUBROUTINE mulliken 01547 01548 ! ***************************************************************************** 01581 SUBROUTINE mulliken_charges(qs_env,orthonormal_basis,matrix_h,matrix_w,& 01582 should_output,output_unit,print_level,error) 01583 01584 TYPE(qs_environment_type), POINTER :: qs_env 01585 LOGICAL, INTENT(IN) :: orthonormal_basis 01586 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01587 OPTIONAL, POINTER :: matrix_h, matrix_w 01588 LOGICAL, INTENT(IN) :: should_output 01589 INTEGER, INTENT(IN) :: output_unit, print_level 01590 TYPE(cp_error_type), INTENT(INOUT) :: error 01591 01592 CHARACTER(LEN=*), PARAMETER :: routineN = 'mulliken_charges', 01593 routineP = moduleN//':'//routineN 01594 01595 CHARACTER(LEN=10) :: spin_info 01596 CHARACTER(LEN=6), ALLOCATABLE, 01597 DIMENSION(:) :: symbol 01598 CHARACTER(LEN=default_string_length) :: atomic_kind_name 01599 INTEGER :: atom_a, blk, handle, i, iatom, ikind, isb, iset, isgf, ishell, 01600 ispin, jatom, jsgf, lu, m, natom, natom_of_kind, nkind, nsb, nset, 01601 nsgf, nspin, sgf, stat 01602 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf_atom 01603 INTEGER, DIMENSION(:), POINTER :: atom_list, nshell 01604 INTEGER, DIMENSION(:, :), POINTER :: first_sgf, l, last_sgf 01605 LOGICAL :: dft_plus_u_atom, failure, 01606 found, just_energy 01607 REAL(KIND=dp) :: eps_u_ramping, fspin, q, 01608 u_minus_j, u_minus_j_target, 01609 u_ramping, v 01610 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dEdq, trps 01611 REAL(KIND=dp), ALLOCATABLE, 01612 DIMENSION(:, :) :: q_ii 01613 REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, p_block, s_block, 01614 w_block 01615 TYPE(atomic_kind_type), DIMENSION(:), 01616 POINTER :: atomic_kind_set 01617 TYPE(atomic_kind_type), POINTER :: atomic_kind 01618 TYPE(cp_dbcsr_iterator) :: iter 01619 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01620 POINTER :: matrix_p, matrix_s 01621 TYPE(cp_dbcsr_type), POINTER :: sm_h, sm_p, sm_s, sm_w 01622 TYPE(dft_control_type), POINTER :: dft_control 01623 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 01624 TYPE(particle_type), DIMENSION(:), 01625 POINTER :: particle_set 01626 TYPE(qs_energy_type), POINTER :: energy 01627 TYPE(qs_rho_type), POINTER :: rho 01628 01629 CALL timeset(routineN,handle) 01630 01631 NULLIFY (atom_list) 01632 NULLIFY (atomic_kind) 01633 NULLIFY (atomic_kind_set) 01634 NULLIFY (dft_control) 01635 NULLIFY (energy) 01636 NULLIFY (first_sgf) 01637 NULLIFY (h_block) 01638 NULLIFY (matrix_p) 01639 NULLIFY (matrix_s) 01640 NULLIFY (l) 01641 NULLIFY (last_sgf) 01642 NULLIFY (nshell) 01643 NULLIFY (orb_basis_set) 01644 NULLIFY (p_block) 01645 NULLIFY (particle_set) 01646 NULLIFY (rho) 01647 NULLIFY (s_block) 01648 NULLIFY (sm_h) 01649 NULLIFY (sm_p) 01650 NULLIFY (sm_s) 01651 NULLIFY (w_block) 01652 01653 failure = .FALSE. 01654 01655 CALL get_qs_env(qs_env=qs_env,& 01656 atomic_kind_set=atomic_kind_set,& 01657 dft_control=dft_control,& 01658 energy=energy,& 01659 particle_set=particle_set,& 01660 rho=rho,& 01661 error=error) 01662 01663 CPPrecondition(ASSOCIATED(atomic_kind_set),cp_failure_level,routineP,error,failure) 01664 CPPrecondition(ASSOCIATED(dft_control),cp_failure_level,routineP,error,failure) 01665 CPPrecondition(ASSOCIATED(energy),cp_failure_level,routineP,error,failure) 01666 CPPrecondition(ASSOCIATED(particle_set),cp_failure_level,routineP,error,failure) 01667 CPPrecondition(ASSOCIATED(rho),cp_failure_level,routineP,error,failure) 01668 01669 IF (orthonormal_basis) THEN 01670 NULLIFY (sm_s) 01671 ELSE 01672 ! Get overlap matrix in sparse format 01673 CALL get_qs_env(qs_env=qs_env,& 01674 matrix_s=matrix_s,& 01675 error=error) 01676 CPPrecondition(ASSOCIATED(matrix_s),cp_failure_level,routineP,error,failure) 01677 sm_s => matrix_s(1)%matrix 01678 END IF 01679 01680 ! Get density matrices in sparse format 01681 01682 matrix_p => rho%rho_ao 01683 01684 energy%dft_plus_u = 0.0_dp 01685 01686 nspin = dft_control%nspins 01687 01688 IF (nspin == 2) THEN 01689 fspin = 1.0_dp 01690 ELSE 01691 fspin = 0.5_dp 01692 END IF 01693 01694 ! Get the total number of atoms, contracted spherical Gaussian basis 01695 ! functions, and atomic kinds 01696 01697 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 01698 natom=natom,& 01699 nsgf=nsgf) 01700 01701 nkind = SIZE(atomic_kind_set) 01702 01703 ALLOCATE (first_sgf_atom(natom),STAT=stat) 01704 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01705 first_sgf_atom(:) = 0 01706 01707 CALL get_particle_set(particle_set=particle_set,& 01708 first_sgf=first_sgf_atom,& 01709 error=error) 01710 01711 ALLOCATE (trps(nsgf),STAT=stat) 01712 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01713 trps(:) = 0.0_dp 01714 01715 IF (PRESENT(matrix_h).OR.PRESENT(matrix_w)) THEN 01716 ALLOCATE (dEdq(nsgf),STAT=stat) 01717 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01718 just_energy = .FALSE. 01719 ELSE 01720 just_energy = .TRUE. 01721 END IF 01722 01723 ! Loop over all spins 01724 01725 DO ispin=1,nspin 01726 01727 IF (PRESENT(matrix_h)) THEN 01728 ! Hamiltonian matrix for spin ispin in sparse format 01729 sm_h => matrix_h(ispin)%matrix 01730 dEdq(:) = 0.0_dp 01731 ELSE 01732 NULLIFY (sm_h) 01733 END IF 01734 01735 IF (PRESENT(matrix_w)) THEN 01736 ! Energy weighted density matrix for spin ispin in sparse format 01737 sm_w => matrix_w(ispin)%matrix 01738 dEdq(:) = 0.0_dp 01739 ELSE 01740 NULLIFY (sm_w) 01741 END IF 01742 01743 sm_p => matrix_p(ispin)%matrix ! Density matrix for spin ispin in sparse format 01744 01745 ! Calculate Trace(P*S) assuming symmetric matrices 01746 01747 trps(:) = 0.0_dp 01748 01749 CALL cp_dbcsr_iterator_start(iter,sm_p) 01750 01751 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 01752 01753 CALL cp_dbcsr_iterator_next_block(iter,iatom,jatom,p_block,blk) 01754 01755 IF (orthonormal_basis) THEN 01756 01757 IF (iatom /= jatom) CYCLE 01758 01759 IF (ASSOCIATED(p_block)) THEN 01760 sgf = first_sgf_atom(iatom) 01761 DO isgf=1,SIZE(p_block,1) 01762 trps(sgf) = trps(sgf) + p_block(isgf,isgf) 01763 sgf = sgf + 1 01764 END DO 01765 END IF 01766 01767 ELSE 01768 01769 CALL cp_dbcsr_get_block_p(matrix=sm_s,& 01770 row=iatom,& 01771 col=jatom,& 01772 block=s_block,& 01773 found=found) 01774 CPPostcondition(ASSOCIATED(s_block),cp_failure_level,routineP,error,failure) 01775 01776 sgf = first_sgf_atom(jatom) 01777 DO jsgf=1,SIZE(p_block,2) 01778 DO isgf=1,SIZE(p_block,1) 01779 trps(sgf) = trps(sgf) + p_block(isgf,jsgf)*s_block(isgf,jsgf) 01780 END DO 01781 sgf = sgf + 1 01782 END DO 01783 01784 IF (iatom /= jatom) THEN 01785 sgf = first_sgf_atom(iatom) 01786 DO isgf=1,SIZE(p_block,1) 01787 DO jsgf=1,SIZE(p_block,2) 01788 trps(sgf) = trps(sgf) + p_block(isgf,jsgf)*s_block(isgf,jsgf) 01789 END DO 01790 sgf = sgf + 1 01791 END DO 01792 END IF 01793 01794 END IF ! orthonormal basis set 01795 01796 END DO ! next atom "iatom" 01797 01798 CALL cp_dbcsr_iterator_stop(iter) 01799 01800 CALL mp_sum(trps,qs_env%para_env%group) 01801 01802 ! q <- Trace(PS) 01803 01804 ! E[DFT+U] = E[DFT] + E[U] 01805 ! = E[DFT] + (U - J)*(q - q**2))/2 01806 01807 ! V(i,j)[DFT+U] = V(i,j)[DFT] + V(i,j)[U] 01808 ! = dE[DFT]/dP(i,j) + dE[U]/dP(i,j) 01809 ! = dE[DFT]/dP(i,j) + (dE(U)/dq)*(dq/dP(i,j)) 01810 01811 ! Loop over all atomic kinds 01812 01813 DO ikind=1,nkind 01814 01815 atomic_kind => atomic_kind_set(ikind) 01816 01817 ! Load the required atomic kind data 01818 01819 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01820 atom_list=atom_list,& 01821 dft_plus_u_atom=dft_plus_u_atom,& 01822 l_of_dft_plus_u=lu,& 01823 name=atomic_kind_name,& 01824 natom=natom_of_kind,& 01825 orb_basis_set=orb_basis_set,& 01826 u_minus_j=u_minus_j,& 01827 u_minus_j_target=u_minus_j_target,& 01828 u_ramping=u_ramping,& 01829 eps_u_ramping=eps_u_ramping) 01830 01831 ! Check, if this atom needs a DFT+U correction 01832 01833 IF (.NOT.ASSOCIATED(orb_basis_set)) CYCLE 01834 IF (.NOT.dft_plus_u_atom) CYCLE 01835 IF (lu < 0) CYCLE 01836 01837 ! Apply U ramping if requested 01838 01839 IF ((ispin == 1).AND.(u_ramping > 0.0_dp)) THEN 01840 IF (qs_env%scf_env%iter_delta <= eps_u_ramping) THEN 01841 u_minus_j = MIN(u_minus_j + u_ramping,u_minus_j_target) 01842 CALL set_atomic_kind(atomic_kind=atomic_kind,u_minus_j=u_minus_j) 01843 END IF 01844 IF (should_output.AND.(output_unit > 0)) THEN 01845 WRITE (UNIT=output_unit,FMT="(T3,A,3X,A,F0.3,A)")& 01846 "Kind name: "//TRIM(ADJUSTL(atomic_kind_name)),& 01847 "U(eff) = ",u_minus_j*evolt," eV" 01848 END IF 01849 END IF 01850 01851 IF (u_minus_j == 0.0_dp) CYCLE 01852 01853 ! Load the required Gaussian basis set data 01854 01855 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01856 first_sgf=first_sgf,& 01857 l=l,& 01858 last_sgf=last_sgf,& 01859 nset=nset,& 01860 nshell=nshell) 01861 01862 ! Count the relevant shell blocks of this atomic kind 01863 01864 nsb = 0 01865 DO iset=1,nset 01866 DO ishell=1,nshell(iset) 01867 IF (l(ishell,iset) == lu) nsb = nsb + 1 01868 END DO 01869 END DO 01870 01871 ALLOCATE (q_ii(nsb,2*lu+1),STAT=stat) 01872 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01873 01874 ! Print headline if requested 01875 01876 IF (should_output.AND.(print_level > low_print_level)) THEN 01877 IF (output_unit > 0) THEN 01878 ALLOCATE (symbol(2*lu+1),STAT=stat) 01879 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01880 DO m=-lu,lu 01881 symbol(lu+m+1) = sgf_symbol(0,lu,m) 01882 END DO 01883 IF (nspin > 1) THEN 01884 WRITE (UNIT=spin_info,FMT="(A8,I2)") " of spin",ispin 01885 ELSE 01886 spin_info = "" 01887 END IF 01888 WRITE (UNIT=output_unit,FMT="(/,T3,A,I0,A,/,/,T5,A,10(2X,A6))")& 01889 "DFT+U occupations"//TRIM(spin_info)//" for the atoms of atomic kind ",ikind,& 01890 ": "//TRIM(atomic_kind_name),& 01891 "Atom Shell ",(ADJUSTR(symbol(i)),i=1,2*lu+1)," Trace" 01892 DEALLOCATE (symbol,STAT=stat) 01893 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01894 END IF 01895 END IF 01896 01897 ! Loop over all atoms of the current atomic kind 01898 01899 DO iatom=1,natom_of_kind 01900 01901 atom_a = atom_list(iatom) 01902 01903 q_ii(:,:) = 0.0_dp 01904 01905 ! Get diagonal block 01906 01907 CALL cp_dbcsr_get_block_p(matrix=sm_p,& 01908 row=atom_a,& 01909 col=atom_a,& 01910 block=p_block,& 01911 found=found) 01912 01913 ! Calculate E(U) and dE(U)/dq 01914 01915 IF (ASSOCIATED(p_block)) THEN 01916 01917 sgf = first_sgf_atom(atom_a) 01918 01919 isb = 0 01920 DO iset=1,nset 01921 DO ishell=1,nshell(iset) 01922 IF (l(ishell,iset) == lu) THEN 01923 isb = isb + 1 01924 i = 0 01925 DO isgf=first_sgf(ishell,iset),last_sgf(ishell,iset) 01926 q = fspin*trps(sgf) 01927 i = i + 1 01928 q_ii(isb,i) = q 01929 energy%dft_plus_u = energy%dft_plus_u +& 01930 0.5_dp*u_minus_j*(q - q**2)/fspin 01931 IF (.NOT.just_energy) THEN 01932 dEdq(sgf) = dEdq(sgf) + u_minus_j*(0.5_dp - q) 01933 END IF 01934 sgf = sgf + 1 01935 END DO ! next contracted spherical Gaussian function "isgf" 01936 ELSE 01937 sgf = sgf + last_sgf(ishell,iset) - first_sgf(ishell,iset) + 1 01938 END IF ! angular momentum requested for DFT+U correction 01939 END DO ! next shell "ishell" 01940 END DO ! next shell set "iset" 01941 01942 END IF ! this process is the owner of the sparse matrix block? 01943 01944 ! Consider print requests 01945 01946 IF (should_output.AND.(print_level > low_print_level)) THEN 01947 CALL mp_sum(q_ii,qs_env%para_env%group) 01948 IF (output_unit > 0) THEN 01949 DO isb=1,nsb 01950 WRITE (UNIT=output_unit,FMT="(T3,I6,2X,I6,2X,10F8.3)")& 01951 atom_a,isb,q_ii(isb,:),SUM(q_ii(isb,:)) 01952 END DO 01953 WRITE (UNIT=output_unit,FMT="(T12,A,2X,10F8.3)")& 01954 "Total",(SUM(q_ii(:,i)),i=1,2*lu+1),SUM(q_ii) 01955 WRITE (UNIT=output_unit,FMT="(A)") "" 01956 END IF 01957 END IF ! should output 01958 01959 END DO ! next atom "iatom" of atomic kind "ikind" 01960 01961 IF (ALLOCATED(q_ii)) THEN 01962 DEALLOCATE (q_ii,STAT=stat) 01963 CPPostcondition((stat == 0),cp_failure_level,routineP,error,failure) 01964 END IF 01965 01966 END DO ! next atomic kind "ikind" 01967 01968 IF (.NOT.just_energy) THEN 01969 CALL mp_sum(dEdq,qs_env%para_env%group) 01970 END IF 01971 01972 ! Add V(i,j)[U] to V(i,j)[DFT] 01973 01974 IF (ASSOCIATED(sm_h)) THEN 01975 01976 CALL cp_dbcsr_iterator_start(iter,sm_h) 01977 01978 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 01979 01980 CALL cp_dbcsr_iterator_next_block(iter,iatom,jatom,h_block,blk) 01981 01982 IF (orthonormal_basis) THEN 01983 01984 IF (iatom /= jatom) CYCLE 01985 01986 IF (ASSOCIATED(h_block)) THEN 01987 sgf = first_sgf_atom(iatom) 01988 DO isgf=1,SIZE(h_block,1) 01989 h_block(isgf,isgf) = h_block(isgf,isgf) + dEdq(sgf) 01990 sgf = sgf + 1 01991 END DO 01992 END IF 01993 01994 ELSE 01995 01996 ! Request katom just to check for consistent sparse matrix pattern 01997 01998 CALL cp_dbcsr_get_block_p(matrix=sm_s,& 01999 row=iatom,& 02000 col=jatom,& 02001 block=s_block,& 02002 found=found) 02003 CPPostcondition(ASSOCIATED(s_block),cp_failure_level,routineP,error,failure) 02004 02005 ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation 02006 02007 sgf = first_sgf_atom(iatom) 02008 02009 DO isgf=1,SIZE(h_block,1) 02010 IF (dEdq(sgf) /= 0.0_dp) THEN 02011 v = 0.5_dp*dEdq(sgf) 02012 DO jsgf=1,SIZE(h_block,2) 02013 h_block(isgf,jsgf) = h_block(isgf,jsgf) + v*s_block(isgf,jsgf) 02014 END DO 02015 END IF 02016 sgf = sgf + 1 02017 END DO 02018 02019 sgf = first_sgf_atom(jatom) 02020 02021 DO jsgf=1,SIZE(h_block,2) 02022 IF (dEdq(sgf) /= 0.0_dp) THEN 02023 v = 0.5_dp*dEdq(sgf) 02024 DO isgf=1,SIZE(h_block,1) 02025 h_block(isgf,jsgf) = h_block(isgf,jsgf) + v*s_block(isgf,jsgf) 02026 END DO 02027 END IF 02028 sgf = sgf + 1 02029 END DO 02030 02031 END IF ! orthonormal basis set 02032 02033 END DO ! Next atom "iatom" 02034 02035 CALL cp_dbcsr_iterator_stop(iter) 02036 02037 END IF ! An update of the Hamiltonian matrix is requested 02038 02039 ! Calculate the contribution (non-Pulay part) to the derivatives 02040 ! w.r.t. the nuclear positions, which requires an update of the 02041 ! energy weighted density W. 02042 02043 IF (ASSOCIATED(sm_w).AND.(.NOT.orthonormal_basis)) THEN 02044 02045 CALL cp_dbcsr_iterator_start(iter,sm_p) 02046 02047 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 02048 02049 CALL cp_dbcsr_iterator_next_block(iter,iatom,jatom,p_block,blk) 02050 02051 ! Skip the diagonal blocks of the W matrix 02052 02053 IF (iatom == jatom) CYCLE 02054 02055 ! Request katom just to check for consistent sparse matrix patterns 02056 02057 CALL cp_dbcsr_get_block_p(matrix=sm_w,& 02058 row=iatom,& 02059 col=jatom,& 02060 block=w_block,& 02061 found=found) 02062 CPPostcondition(ASSOCIATED(w_block),cp_failure_level,routineP,error,failure) 02063 02064 ! Consider the symmetric form 1/2*(P*S + S*P) for the calculation 02065 02066 sgf = first_sgf_atom(iatom) 02067 02068 DO isgf=1,SIZE(w_block,1) 02069 IF (dEdq(sgf) /= 0.0_dp) THEN 02070 v = -0.5_dp*dEdq(sgf) 02071 DO jsgf=1,SIZE(w_block,2) 02072 w_block(isgf,jsgf) = w_block(isgf,jsgf) + v*p_block(isgf,jsgf) 02073 END DO 02074 END IF 02075 sgf = sgf + 1 02076 END DO 02077 02078 sgf = first_sgf_atom(jatom) 02079 02080 DO jsgf=1,SIZE(w_block,2) 02081 IF (dEdq(sgf) /= 0.0_dp) THEN 02082 v = -0.5_dp*dEdq(sgf) 02083 DO isgf=1,SIZE(w_block,1) 02084 w_block(isgf,jsgf) = w_block(isgf,jsgf) + v*p_block(isgf,jsgf) 02085 END DO 02086 END IF 02087 sgf = sgf + 1 02088 END DO 02089 02090 END DO ! next block node "jatom" 02091 02092 CALL cp_dbcsr_iterator_stop(iter) 02093 02094 END IF ! W matrix update requested 02095 02096 END DO ! next spin "ispin" 02097 02098 ! Collect the energy contributions from all processes 02099 02100 CALL mp_sum(energy%dft_plus_u,qs_env%para_env%group) 02101 02102 CALL cp_assert(energy%dft_plus_u >= 0.0_dp,cp_warning_level,& 02103 cp_assertion_failed,routineP,& 02104 "DFT+U energy contibution is negative possibly due "//& 02105 "to unphysical Mulliken charges. Check your input, "//& 02106 "if this warning persists or try a different method!",& 02107 only_ionode=.TRUE.) 02108 02109 ! Release local work storage 02110 02111 IF (ALLOCATED(first_sgf_atom)) THEN 02112 DEALLOCATE (first_sgf_atom,STAT=stat) 02113 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02114 END IF 02115 02116 IF (ALLOCATED(trps)) THEN 02117 DEALLOCATE (trps,STAT=stat) 02118 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02119 END IF 02120 02121 IF (ALLOCATED(dEdq)) THEN 02122 DEALLOCATE (dEdq,STAT=stat) 02123 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02124 END IF 02125 02126 CALL timestop(handle) 02127 02128 END SUBROUTINE mulliken_charges 02129 02130 END MODULE dft_plus_u
1.7.3