|
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 ! ***************************************************************************** 00018 MODULE qs_mo_types 00019 00020 USE atomic_kind_types, ONLY: atomic_kind_type,& 00021 get_atomic_kind,& 00022 get_atomic_kind_set 00023 USE basis_set_types, ONLY: get_gto_basis_set,& 00024 gto_basis_set_type 00025 USE cp_dbcsr_interface, ONLY: cp_dbcsr_copy,& 00026 cp_dbcsr_init_p,& 00027 cp_dbcsr_release_p 00028 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 00029 cp_dbcsr_copy_columns_hack 00030 USE cp_dbcsr_types, ONLY: cp_dbcsr_type 00031 USE cp_files, ONLY: close_file,& 00032 open_file 00033 USE cp_fm_pool_types, ONLY: cp_fm_pool_type,& 00034 fm_pool_create_fm 00035 USE cp_fm_types, ONLY: & 00036 cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_p_type, & 00037 cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, & 00038 cp_fm_type, cp_fm_write_unformatted 00039 USE cp_output_handling, ONLY: cp_p_file,& 00040 cp_print_key_finished_output,& 00041 cp_print_key_generate_filename,& 00042 cp_print_key_should_output,& 00043 cp_print_key_unit_nr 00044 USE cp_para_types, ONLY: cp_para_env_type 00045 USE fermi_utils, ONLY: FermiFixed,& 00046 FermiFixedDeriv 00047 USE input_constants, ONLY: smear_energy_window,& 00048 smear_fermi_dirac,& 00049 smear_list 00050 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00051 section_vals_type,& 00052 section_vals_val_get 00053 USE kahan_sum, ONLY: accurate_sum 00054 USE kinds, ONLY: default_path_length,& 00055 default_string_length,& 00056 dp 00057 USE message_passing, ONLY: mp_bcast 00058 USE orbital_pointers, ONLY: indco,& 00059 nco,& 00060 nso 00061 USE orbital_symbols, ONLY: cgf_symbol,& 00062 sgf_symbol 00063 USE orbital_transformation_matrices, ONLY: orbtramat 00064 USE particle_types, ONLY: particle_type 00065 USE physcon, ONLY: evolt 00066 USE qs_dftb_types, ONLY: qs_dftb_atom_type 00067 USE qs_dftb_utils, ONLY: get_dftb_atom_param 00068 USE scf_control_types, ONLY: smear_type 00069 USE scp_environment_types, ONLY: scp_environment_type 00070 USE scp_restarts, ONLY: read_aux_coeff_set,& 00071 write_scp_coeff_set 00072 USE string_utilities, ONLY: compress 00073 USE timings, ONLY: timeset,& 00074 timestop 00075 USE util, ONLY: sort 00076 USE xas_env_types, ONLY: get_xas_env,& 00077 xas_environment_type 00078 #include "cp_common_uses.h" 00079 00080 IMPLICIT NONE 00081 00082 PRIVATE 00083 00084 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_mo_types' 00085 00086 TYPE mo_set_type 00087 ! the actual MO coefficients as a matrix 00088 TYPE(cp_fm_type), POINTER :: mo_coeff 00089 TYPE(cp_dbcsr_type), POINTER :: mo_coeff_b 00090 ! we are using the dbcsr mo_coeff_b 00091 LOGICAL :: use_mo_coeff_b 00092 ! number of molecular orbitals (# cols in mo_coeff) 00093 INTEGER :: nmo 00094 ! number of atomic orbitals (# rows in mo_coeff) 00095 INTEGER :: nao 00096 ! occupation - eigenvalues of the nmo states (if eigenstates) 00097 REAL(KIND = dp), DIMENSION(:), POINTER :: eigenvalues,occupation_numbers 00098 ! maximum allowed occupation number of an MO (1 or 2) 00099 REAL(KIND = dp) :: maxocc 00100 ! number of electrons (taking occupation into account) 00101 INTEGER :: nelectron 00102 REAL(KIND=dp) :: n_el_f 00103 ! highest non-zero occupied orbital 00104 INTEGER :: homo 00105 ! lowest non maxocc occupied orbital (e.g. fractional or zero) 00106 INTEGER :: lfomo 00107 ! flag that indicates if the MOS have the same occupation number 00108 LOGICAL :: uniform_occupation 00109 ! the entropic energy contribution 00110 REAL(KIND=dp) :: kTS 00111 ! Fermi energy level 00112 REAL(KIND=dp) :: mu 00113 ! Threshold value for multiplicity change 00114 REAL(KIND=dp) :: flexible_electron_count 00115 END TYPE mo_set_type 00116 00117 TYPE mo_set_p_type 00118 TYPE(mo_set_type), POINTER :: mo_set 00119 END TYPE mo_set_p_type 00120 00121 PUBLIC :: mo_set_p_type,& 00122 mo_set_type 00123 00124 PUBLIC :: allocate_mo_set,& 00125 correct_mo_eigenvalues,& 00126 deallocate_mo_set,& 00127 get_mo_set,& 00128 init_mo_set,& 00129 wfn_restart_file_name,& 00130 read_mo_set,& 00131 read_mos_restart_low,& 00132 set_mo_occupation,& 00133 set_mo_set, & 00134 write_mo_set, & 00135 write_mo_set_low,& 00136 mo_set_restrict,& 00137 write_rt_mos_to_restart,& 00138 read_rt_mos_from_restart,& 00139 duplicate_mo_set 00140 00141 INTERFACE read_mo_set 00142 MODULE PROCEDURE read_mo_set_from_restart 00143 END INTERFACE 00144 00145 INTERFACE set_mo_occupation 00146 MODULE PROCEDURE set_mo_occupation_1,set_mo_occupation_2 00147 END INTERFACE 00148 00149 INTERFACE write_mo_set 00150 MODULE PROCEDURE write_mo_set_to_output_unit,write_mo_set_to_restart 00151 END INTERFACE 00152 00153 CONTAINS 00154 00155 ! ***************************************************************************** 00163 SUBROUTINE set_mo_occupation_3(mo_array,smear,error) 00164 00165 TYPE(mo_set_p_type), DIMENSION(:), 00166 POINTER :: mo_array 00167 TYPE(smear_type), POINTER :: smear 00168 TYPE(cp_error_type), INTENT(INOUT) :: error 00169 00170 CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mo_occupation_3', 00171 routineP = moduleN//':'//routineN 00172 00173 INTEGER :: all_nmo, handle, homo_a, 00174 homo_b, i, lfomo_a, lfomo_b, 00175 nmo_a, nmo_b, stat, xas_estate 00176 INTEGER, ALLOCATABLE, DIMENSION(:) :: all_index 00177 LOGICAL :: failure, is_large 00178 REAL(KIND=dp) :: all_nelec, kTS, mu, nelec_a, 00179 nelec_b, occ_estate 00180 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: all_eigval, all_occ 00181 REAL(KIND=dp), DIMENSION(:), POINTER :: eigval_a, eigval_b, occ_a, 00182 occ_b 00183 00184 failure = .FALSE. 00185 CPPrecondition(ASSOCIATED(mo_array),cp_failure_level,routineP,error,failure) 00186 CPPrecondition((SIZE(mo_array) == 2),cp_failure_level,routineP,error,failure) 00187 CALL timeset(routineN,handle) 00188 00189 NULLIFY(eigval_a,eigval_b,occ_a,occ_b) 00190 CALL get_mo_set(mo_set=mo_array(1)%mo_set,nmo=nmo_a,eigenvalues=eigval_a, & 00191 occupation_numbers=occ_a) 00192 CALL get_mo_set(mo_set=mo_array(2)%mo_set,nmo=nmo_b,eigenvalues=eigval_b, & 00193 occupation_numbers=occ_b) 00194 all_nmo = nmo_a+nmo_b 00195 ALLOCATE(all_eigval(all_nmo), STAT=stat) 00196 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00197 ALLOCATE(all_occ(all_nmo), STAT=stat) 00198 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00199 ALLOCATE(all_index(all_nmo), STAT=stat) 00200 00201 all_eigval(1:nmo_a) = eigval_a(1:nmo_a) 00202 all_eigval(nmo_a+1:all_nmo) = eigval_b(1:nmo_b) 00203 00204 CALL sort(all_eigval,all_nmo,all_index) 00205 00206 xas_estate = -1 00207 occ_estate = 0.0_dp 00208 00209 nelec_a=0.0_dp 00210 nelec_b=0.0_dp 00211 all_nelec=0.0_dp 00212 nelec_a=accurate_sum(occ_a(:)) 00213 nelec_b=accurate_sum(occ_b(:)) 00214 all_nelec = nelec_a+nelec_b 00215 00216 DO i = 1,all_nmo 00217 IF(all_index(i)<=nmo_a) THEN 00218 all_occ(i) = occ_a(all_index(i)) 00219 ELSE 00220 all_occ(i) = occ_b(all_index(i)-nmo_a) 00221 END IF 00222 END DO 00223 00224 CALL FermiFixed(all_occ,mu,kTS,all_eigval,all_nelec, & 00225 smear%electronic_temperature,1._dp,xas_estate,occ_estate) 00226 00227 is_large=ABS(MAXVAL(all_occ)-1.0_dp)> smear%eps_fermi_dirac 00228 ! this is not a real problem, but the temperature might be a bit large 00229 CALL cp_assert(.NOT.is_large,cp_warning_level,cp_assertion_failed,routineP,& 00230 "Fermi-Dirac smearing includes the first MO",& 00231 error,failure) 00232 00233 is_large=ABS(MINVAL(all_occ))> smear%eps_fermi_dirac 00234 CALL cp_assert(.NOT.is_large,cp_warning_level,cp_assertion_failed,routineP,& 00235 "Fermi-Dirac smearing includes the last MO => "//& 00236 "Add more MOs for proper smearing.",error,failure) 00237 00238 ! check that the total electron count is accurate 00239 is_large=(ABS(all_nelec - accurate_sum(all_occ(:))) > smear%eps_fermi_dirac*all_nelec) 00240 CALL cp_assert(.NOT.is_large,& 00241 cp_warning_level,cp_assertion_failed,routineP,& 00242 "Total number of electrons is not accurate",& 00243 error,failure) 00244 00245 DO i = 1,all_nmo 00246 IF(all_index(i)<=nmo_a) THEN 00247 occ_a(all_index(i)) = all_occ(i) 00248 eigval_a(all_index(i)) = all_eigval(i) 00249 ELSE 00250 occ_b(all_index(i)-nmo_a) = all_occ(i) 00251 eigval_b(all_index(i)-nmo_a) = all_eigval(i) 00252 END IF 00253 END DO 00254 00255 nelec_a=accurate_sum(occ_a(:)) 00256 nelec_b=accurate_sum(occ_b(:)) 00257 DO i =1,nmo_a 00258 IF (occ_a(i) < 1.0_dp) THEN 00259 lfomo_a=i 00260 EXIT 00261 END IF 00262 END DO 00263 DO i =1,nmo_b 00264 IF (occ_b(i) < 1.0_dp) THEN 00265 lfomo_b=i 00266 EXIT 00267 END IF 00268 END DO 00269 DO i = nmo_a,lfomo_a,-1 00270 IF (occ_a(i) > smear%eps_fermi_dirac) THEN 00271 homo_a=i 00272 EXIT 00273 END IF 00274 END DO 00275 DO i = nmo_b,lfomo_b,-1 00276 IF (occ_b(i) > smear%eps_fermi_dirac) THEN 00277 homo_b=i 00278 EXIT 00279 END IF 00280 END DO 00281 00282 CALL set_mo_set(mo_set=mo_array(1)%mo_set,kTS=kTS/2.0_dp,mu=mu,n_el_f=nelec_a,& 00283 lfomo=lfomo_a,homo=homo_a,uniform_occupation=.FALSE., error=error) 00284 CALL set_mo_set(mo_set=mo_array(2)%mo_set,kTS=kTS/2.0_dp,mu=mu,n_el_f=nelec_b,& 00285 lfomo=lfomo_b,homo=homo_b,uniform_occupation=.FALSE., error=error) 00286 00287 00288 CALL timestop(handle) 00289 00290 END SUBROUTINE set_mo_occupation_3 00291 ! ***************************************************************************** 00298 SUBROUTINE set_mo_occupation_2(mo_array,smear,eval_deriv,error) 00299 00300 TYPE(mo_set_p_type), DIMENSION(:), 00301 POINTER :: mo_array 00302 TYPE(smear_type), POINTER :: smear 00303 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00304 POINTER :: eval_deriv 00305 TYPE(cp_error_type), INTENT(INOUT) :: error 00306 00307 CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mo_occupation_2', 00308 routineP = moduleN//':'//routineN 00309 00310 INTEGER :: handle, i, lumo_a, lumo_b, 00311 multiplicity_new, 00312 multiplicity_old, nelec 00313 LOGICAL :: failure 00314 REAL(KIND=dp) :: threshold 00315 REAL(KIND=dp), DIMENSION(:), POINTER :: eigval_a, eigval_b 00316 TYPE(mo_set_type), POINTER :: mo_set_a, mo_set_b 00317 00318 failure = .FALSE. 00319 CPPrecondition(ASSOCIATED(mo_array),cp_failure_level,routineP,error,failure) 00320 mo_set_a => mo_array(1)%mo_set 00321 CPPrecondition(ASSOCIATED(mo_set_a),cp_failure_level,routineP,error,failure) 00322 ! Fall back for the case that we have only one MO set 00323 IF (SIZE(mo_array) == 1) THEN 00324 IF (PRESENT(eval_deriv)) THEN 00325 CALL set_mo_occupation_1(mo_set_a,smear=smear,eval_deriv=eval_deriv,error=error) 00326 ELSE 00327 CALL set_mo_occupation_1(mo_set_a,smear=smear,error=error) 00328 END IF 00329 RETURN 00330 END IF 00331 CPPrecondition((SIZE(mo_array) == 2),cp_failure_level,routineP,error,failure) 00332 mo_set_b => mo_array(2)%mo_set 00333 00334 IF(smear%do_smear .AND. smear%fixed_mag_mom<0.0_dp) THEN 00335 CPPrecondition(.NOT.(PRESENT(eval_deriv)),cp_failure_level,routineP,error,failure) 00336 CALL set_mo_occupation_3(mo_array,smear=smear,error=error) 00337 RETURN 00338 ELSEIF(smear%do_smear .AND. smear%fixed_mag_mom >=0.0) THEN 00339 nelec = mo_set_a%n_el_f + mo_set_b%n_el_f 00340 IF(ABS((mo_set_a%n_el_f-mo_set_b%n_el_f)-smear%fixed_mag_mom)>& 00341 smear%eps_fermi_dirac*nelec)THEN 00342 mo_set_a%n_el_f = nelec/2.0_dp + smear%fixed_mag_mom/2.0_dp 00343 mo_set_b%n_el_f = nelec/2.0_dp - smear%fixed_mag_mom/2.0_dp 00344 END IF 00345 CPPrecondition(.NOT.(PRESENT(eval_deriv)),cp_failure_level,routineP,error,failure) 00346 CALL set_mo_occupation_1(mo_set_a,smear=smear,error=error) 00347 CALL set_mo_occupation_1(mo_set_b,smear=smear,error=error) 00348 END IF 00349 00350 00351 IF (.NOT.((mo_set_a%flexible_electron_count > 0.0_dp).AND.& 00352 (mo_set_b%flexible_electron_count > 0.0_dp))) THEN 00353 IF (PRESENT(eval_deriv)) THEN 00354 CALL set_mo_occupation_1(mo_set_a,smear=smear,eval_deriv=eval_deriv,error=error) 00355 CALL set_mo_occupation_1(mo_set_b,smear=smear,eval_deriv=eval_deriv,error=error) 00356 ELSE 00357 CALL set_mo_occupation_1(mo_set_a,smear=smear,error=error) 00358 CALL set_mo_occupation_1(mo_set_b,smear=smear,error=error) 00359 END IF 00360 RETURN 00361 END IF 00362 00363 00364 CALL timeset(routineN,handle) 00365 00366 nelec = mo_set_a%nelectron + mo_set_b%nelectron 00367 00368 multiplicity_old = mo_set_a%nelectron - mo_set_b%nelectron + 1 00369 00370 CALL cp_assert((mo_set_a%nelectron < mo_set_a%nmo),& 00371 cp_warning_level,cp_assertion_failed,routineP,& 00372 "All alpha MOs are occupied. Add more alpha MOs to "//& 00373 "allow for a higher multiplicity",only_ionode=.TRUE.,& 00374 error=error) 00375 CALL cp_assert(((mo_set_b%nelectron < mo_set_b%nmo).OR.& 00376 (mo_set_b%nelectron == mo_set_a%nelectron)),& 00377 cp_warning_level,cp_assertion_failed,routineP,& 00378 "All beta MOs are occupied. Add more beta MOs to "//& 00379 "allow for a lower multiplicity",only_ionode=.TRUE.,& 00380 error=error) 00381 00382 eigval_a => mo_set_a%eigenvalues 00383 eigval_b => mo_set_b%eigenvalues 00384 00385 lumo_a = 1 00386 lumo_b = 1 00387 00388 ! Apply Aufbau principle 00389 DO i=1,nelec 00390 ! Threshold is needed to ensure a preference for alpha occupation in the case 00391 ! of degeneracy 00392 threshold = MAX(mo_set_a%flexible_electron_count,mo_set_b%flexible_electron_count) 00393 IF ((eigval_a(lumo_a) - threshold) < eigval_b(lumo_b)) THEN 00394 lumo_a = lumo_a + 1 00395 ELSE 00396 lumo_b = lumo_b + 1 00397 END IF 00398 IF (lumo_a > mo_set_a%nmo) THEN 00399 CALL cp_assert((i == nelec),& 00400 cp_warning_level,cp_assertion_failed,routineP,& 00401 "All alpha MOs are occupied. Add more alpha MOs to "//& 00402 "allow for a higher multiplicity",only_ionode=.TRUE.,& 00403 error=error) 00404 IF (i < nelec) THEN 00405 lumo_a = lumo_a - 1 00406 lumo_b = lumo_b + 1 00407 END IF 00408 END IF 00409 IF (lumo_b > mo_set_b%nmo) THEN 00410 CALL cp_assert((lumo_b >= lumo_a),& 00411 cp_warning_level,cp_assertion_failed,routineP,& 00412 "All beta MOs are occupied. Add more beta MOs to "//& 00413 "allow for a lower multiplicity",only_ionode=.TRUE.,& 00414 error=error) 00415 IF (i < nelec) THEN 00416 lumo_a = lumo_a + 1 00417 lumo_b = lumo_b - 1 00418 END IF 00419 END IF 00420 END DO 00421 00422 mo_set_a%homo = lumo_a - 1 00423 mo_set_b%homo = lumo_b - 1 00424 00425 IF (mo_set_b%homo > mo_set_a%homo) THEN 00426 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 00427 "More beta ("//& 00428 TRIM(ADJUSTL(cp_to_string(mo_set_b%homo)))//& 00429 ") than alpha ("//& 00430 TRIM(ADJUSTL(cp_to_string(mo_set_a%homo)))//& 00431 ") MOs are occupied. Resorting to low spin state",& 00432 only_ionode=.TRUE.,error=error) 00433 mo_set_a%homo = nelec/2 + MODULO(nelec,2) 00434 mo_set_b%homo = nelec/2 00435 END IF 00436 00437 mo_set_a%nelectron = mo_set_a%homo 00438 mo_set_b%nelectron = mo_set_b%homo 00439 multiplicity_new = mo_set_a%nelectron - mo_set_b%nelectron + 1 00440 00441 CALL cp_assert((multiplicity_new == multiplicity_old),& 00442 cp_note_level,cp_assertion_failed,routineP,& 00443 "Multiplicity changed from "//& 00444 TRIM(ADJUSTL(cp_to_string(multiplicity_old)))//" to "//& 00445 TRIM(ADJUSTL(cp_to_string(multiplicity_new))),& 00446 only_ionode=.TRUE.,error=error,failure=failure) 00447 00448 IF (PRESENT(eval_deriv)) THEN 00449 CALL set_mo_occupation_1(mo_set_a,smear=smear,eval_deriv=eval_deriv,error=error) 00450 CALL set_mo_occupation_1(mo_set_b,smear=smear,eval_deriv=eval_deriv,error=error) 00451 ELSE 00452 CALL set_mo_occupation_1(mo_set_a,smear=smear,error=error) 00453 CALL set_mo_occupation_1(mo_set_b,smear=smear,error=error) 00454 END IF 00455 00456 CALL timestop(handle) 00457 00458 END SUBROUTINE set_mo_occupation_2 00459 00460 ! ***************************************************************************** 00466 SUBROUTINE duplicate_mo_set(mo_set_new,mo_set_old,error) 00467 TYPE(mo_set_type), POINTER :: mo_set_new, mo_set_old 00468 TYPE(cp_error_type), INTENT(inout) :: error 00469 00470 CHARACTER(LEN=*), PARAMETER :: routineN = 'duplicate_mo_set', 00471 routineP = moduleN//':'//routineN 00472 00473 INTEGER :: istat, nmo 00474 LOGICAL :: failure 00475 00476 failure = .FALSE. 00477 00478 ALLOCATE (mo_set_new,STAT=istat) 00479 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00480 00481 mo_set_new%maxocc = mo_set_old%maxocc 00482 mo_set_new%nelectron = mo_set_old%nelectron 00483 mo_set_new%n_el_f = mo_set_old%n_el_f 00484 mo_set_new%nao = mo_set_old%nao 00485 mo_set_new%nmo = mo_set_old%nmo 00486 mo_set_new%homo = mo_set_old%homo 00487 mo_set_new%lfomo = mo_set_old%lfomo 00488 mo_set_new%uniform_occupation = mo_set_old%uniform_occupation 00489 mo_set_new%kTS = mo_set_old%kTS 00490 mo_set_new%mu = mo_set_old%mu 00491 mo_set_new%flexible_electron_count = mo_set_old%flexible_electron_count 00492 00493 nmo = mo_set_new%nmo 00494 00495 NULLIFY(mo_set_new%mo_coeff) 00496 CALL cp_fm_create(mo_set_new%mo_coeff,mo_set_old%mo_coeff%matrix_struct,error=error) 00497 CALL cp_fm_to_fm(mo_set_old%mo_coeff,mo_set_new%mo_coeff,error=error) 00498 00499 NULLIFY(mo_set_new%mo_coeff_b) 00500 IF(ASSOCIATED(mo_set_old%mo_coeff_b)) THEN 00501 CALL cp_dbcsr_init_p(mo_set_new%mo_coeff_b,error=error) 00502 CALL cp_dbcsr_copy(mo_set_new%mo_coeff_b,mo_set_old%mo_coeff_b,error=error) 00503 ENDIF 00504 mo_set_new%use_mo_coeff_b = mo_set_old%use_mo_coeff_b 00505 00506 ALLOCATE (mo_set_new%eigenvalues(nmo),STAT=istat) 00507 CPPostcondition((istat == 0),cp_failure_level,routineP,error,failure) 00508 mo_set_new%eigenvalues = mo_set_old%eigenvalues 00509 00510 ALLOCATE (mo_set_new%occupation_numbers(nmo),STAT=istat) 00511 CPPostcondition((istat == 0),cp_failure_level,routineP,error,failure) 00512 mo_set_new%occupation_numbers = mo_set_old%occupation_numbers 00513 00514 END SUBROUTINE duplicate_mo_set 00515 00516 ! ***************************************************************************** 00533 SUBROUTINE allocate_mo_set(mo_set,nao,nmo,nelectron,n_el_f,maxocc,& 00534 flexible_electron_count,error) 00535 00536 TYPE(mo_set_type), POINTER :: mo_set 00537 INTEGER, INTENT(IN) :: nao, nmo, nelectron 00538 REAL(KIND=dp), INTENT(IN) :: n_el_f, maxocc, 00539 flexible_electron_count 00540 TYPE(cp_error_type), INTENT(INOUT) :: error 00541 00542 CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_mo_set', 00543 routineP = moduleN//':'//routineN 00544 00545 INTEGER :: istat 00546 LOGICAL :: failure 00547 00548 failure = .FALSE. 00549 00550 IF (ASSOCIATED(mo_set)) CALL deallocate_mo_set(mo_set,error) 00551 00552 ALLOCATE (mo_set,STAT=istat) 00553 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00554 00555 mo_set%maxocc = maxocc 00556 mo_set%nelectron = nelectron 00557 mo_set%n_el_f = n_el_f 00558 mo_set%nao = nao 00559 mo_set%nmo = nmo 00560 mo_set%homo = 0 00561 mo_set%lfomo = 0 00562 mo_set%uniform_occupation = .TRUE. 00563 mo_set%kTS = 0.0_dp 00564 mo_set%mu = 0.0_dp 00565 mo_set%flexible_electron_count = flexible_electron_count 00566 00567 NULLIFY (mo_set%eigenvalues) 00568 NULLIFY (mo_set%occupation_numbers) 00569 NULLIFY (mo_set%mo_coeff) 00570 NULLIFY (mo_set%mo_coeff_b) 00571 mo_set%use_mo_coeff_b = .FALSE. 00572 00573 END SUBROUTINE allocate_mo_set 00574 00575 ! ***************************************************************************** 00587 SUBROUTINE init_mo_set(mo_set,fm_pool,name,error) 00588 00589 TYPE(mo_set_type), POINTER :: mo_set 00590 TYPE(cp_fm_pool_type), POINTER :: fm_pool 00591 CHARACTER(LEN=*), INTENT(in) :: name 00592 TYPE(cp_error_type), INTENT(inout) :: error 00593 00594 CHARACTER(LEN=*), PARAMETER :: routineN = 'init_mo_set', 00595 routineP = moduleN//':'//routineN 00596 00597 INTEGER :: istat, nao, nmo 00598 LOGICAL :: failure 00599 00600 failure = .FALSE. 00601 00602 CPPrecondition(ASSOCIATED(mo_set),cp_failure_level,routineP,error,failure) 00603 CPPrecondition(ASSOCIATED(fm_pool),cp_failure_level,routineP,error,failure) 00604 CPPrecondition(.NOT.ASSOCIATED(mo_set%eigenvalues),cp_failure_level,routineP,error,failure) 00605 CPPrecondition(.NOT.ASSOCIATED(mo_set%occupation_numbers),cp_failure_level,routineP,error,failure) 00606 CPPrecondition(.NOT.ASSOCIATED(mo_set%mo_coeff),cp_failure_level,routineP,error,failure) 00607 00608 CALL fm_pool_create_fm(fm_pool,mo_set%mo_coeff,name=name,error=error) 00609 CALL cp_fm_get_info(mo_set%mo_coeff,nrow_global=nao,ncol_global=nmo,error=error) 00610 CPPostcondition((nao >= mo_set%nao),cp_failure_level,routineP,error,failure) 00611 CPPostcondition((nmo >= mo_set%nmo),cp_failure_level,routineP,error,failure) 00612 00613 ALLOCATE (mo_set%eigenvalues(nmo),STAT=istat) 00614 CPPostcondition((istat == 0),cp_failure_level,routineP,error,failure) 00615 mo_set%eigenvalues(:) = 0.0_dp 00616 00617 ALLOCATE (mo_set%occupation_numbers(nmo),STAT=istat) 00618 CPPostcondition((istat == 0),cp_failure_level,routineP,error,failure) 00619 CALL set_mo_occupation_1(mo_set=mo_set,error=error) 00620 00621 END SUBROUTINE init_mo_set 00622 00623 ! ***************************************************************************** 00630 SUBROUTINE mo_set_restrict(mo_array,convert_dbcsr,error) 00631 TYPE(mo_set_p_type), DIMENSION(:), 00632 POINTER :: mo_array 00633 LOGICAL, INTENT(in), OPTIONAL :: convert_dbcsr 00634 TYPE(cp_error_type), INTENT(inout) :: error 00635 00636 CHARACTER(LEN=*), PARAMETER :: routineN = 'mo_set_restrict', 00637 routineP = moduleN//':'//routineN 00638 00639 INTEGER :: handle 00640 LOGICAL :: failure, my_convert_dbcsr 00641 00642 CALL timeset(routineN,handle) 00643 00644 failure = .FALSE. 00645 00646 my_convert_dbcsr = .FALSE. 00647 IF(PRESENT(convert_dbcsr)) my_convert_dbcsr = convert_dbcsr 00648 00649 CPPrecondition(ASSOCIATED(mo_array),cp_failure_level,routineP,error,failure) 00650 CPPrecondition(SIZE(mo_array).EQ.2,cp_failure_level,routineP,error,failure) 00651 CPPrecondition(mo_array(1)%mo_set%nmo>=mo_array(2)%mo_set%nmo,cp_failure_level,routineP,error,failure) 00652 00653 ! first nmo_beta orbitals are copied from alpha to beta 00654 IF (.NOT.failure) THEN 00655 IF(my_convert_dbcsr) THEN!fm->dbcsr 00656 CALL cp_dbcsr_copy_columns_hack(mo_array(2)%mo_set%mo_coeff_b,mo_array(1)%mo_set%mo_coeff_b,&!fm->dbcsr 00657 mo_array(2)%mo_set%nmo,1,1,&!fm->dbcsr 00658 para_env=mo_array(1)%mo_set%mo_coeff%matrix_struct%para_env,&!fm->dbcsr 00659 blacs_env=mo_array(1)%mo_set%mo_coeff%matrix_struct%context,error=error)!fm->dbcsr 00660 ELSE!fm->dbcsr 00661 CALL cp_fm_to_fm(mo_array(1)%mo_set%mo_coeff,mo_array(2)%mo_set%mo_coeff,mo_array(2)%mo_set%nmo) 00662 ENDIF 00663 END IF 00664 00665 CALL timestop(handle) 00666 00667 END SUBROUTINE mo_set_restrict 00668 00669 ! ***************************************************************************** 00675 SUBROUTINE correct_mo_eigenvalues(mo_set,level_shift) 00676 00677 TYPE(mo_set_type), POINTER :: mo_set 00678 REAL(KIND=dp), INTENT(IN) :: level_shift 00679 00680 CHARACTER(LEN=*), PARAMETER :: routineN = 'correct_mo_eigenvalues', 00681 routineP = moduleN//':'//routineN 00682 00683 INTEGER :: imo 00684 00685 IF (level_shift == 0.0_dp) RETURN 00686 00687 DO imo=mo_set%homo+1,mo_set%nmo 00688 mo_set%eigenvalues(imo) = mo_set%eigenvalues(imo) - level_shift 00689 END DO 00690 00691 END SUBROUTINE correct_mo_eigenvalues 00692 00693 ! ***************************************************************************** 00699 SUBROUTINE deallocate_mo_set(mo_set,error) 00700 00701 TYPE(mo_set_type), POINTER :: mo_set 00702 TYPE(cp_error_type), INTENT(inout) :: error 00703 00704 CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_mo_set', 00705 routineP = moduleN//':'//routineN 00706 00707 INTEGER :: istat 00708 LOGICAL :: failure 00709 00710 failure = .FALSE. 00711 00712 IF (ASSOCIATED(mo_set)) THEN 00713 IF (ASSOCIATED(mo_set%eigenvalues)) THEN 00714 DEALLOCATE (mo_set%eigenvalues,STAT=istat) 00715 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00716 END IF 00717 IF (ASSOCIATED(mo_set%occupation_numbers)) THEN 00718 DEALLOCATE (mo_set%occupation_numbers,STAT=istat) 00719 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00720 END IF 00721 CALL cp_fm_release(mo_set%mo_coeff,error=error) 00722 IF(ASSOCIATED(mo_set%mo_coeff_b)) CALL cp_dbcsr_release_p(mo_set%mo_coeff_b, error=error) 00723 DEALLOCATE (mo_set,STAT=istat) 00724 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00725 END IF 00726 00727 END SUBROUTINE deallocate_mo_set 00728 00729 ! ***************************************************************************** 00735 SUBROUTINE get_mo_set(mo_set,maxocc,homo,lfomo,nao,nelectron,n_el_f,nmo,& 00736 eigenvalues,occupation_numbers,mo_coeff,mo_coeff_b,& 00737 uniform_occupation,kTS,mu,flexible_electron_count) 00738 00739 TYPE(mo_set_type), POINTER :: mo_set 00740 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: maxocc 00741 INTEGER, INTENT(OUT), OPTIONAL :: homo, lfomo, nao, nelectron 00742 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: n_el_f 00743 INTEGER, INTENT(OUT), OPTIONAL :: nmo 00744 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00745 POINTER :: eigenvalues, 00746 occupation_numbers 00747 TYPE(cp_fm_type), OPTIONAL, POINTER :: mo_coeff 00748 TYPE(cp_dbcsr_type), OPTIONAL, POINTER :: mo_coeff_b 00749 LOGICAL, INTENT(OUT), OPTIONAL :: uniform_occupation 00750 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: kTS, mu, 00751 flexible_electron_count 00752 00753 IF (PRESENT(maxocc)) maxocc = mo_set%maxocc 00754 IF (PRESENT(homo)) homo = mo_set%homo 00755 IF (PRESENT(lfomo)) lfomo = mo_set%lfomo 00756 IF (PRESENT(nao)) nao = mo_set%nao 00757 IF (PRESENT(nelectron)) nelectron = mo_set%nelectron 00758 IF (PRESENT(n_el_f)) n_el_f = mo_set%n_el_f 00759 IF (PRESENT(nmo)) nmo = mo_set%nmo 00760 IF (PRESENT(eigenvalues)) eigenvalues => mo_set%eigenvalues 00761 IF (PRESENT(occupation_numbers)) THEN 00762 occupation_numbers => mo_set%occupation_numbers 00763 END IF 00764 IF (PRESENT(mo_coeff)) mo_coeff => mo_set%mo_coeff 00765 IF (PRESENT(mo_coeff_b)) mo_coeff_b => mo_set%mo_coeff_b 00766 IF (PRESENT(uniform_occupation)) uniform_occupation = mo_set%uniform_occupation 00767 IF (PRESENT(kTS)) kTS = mo_set%kTS 00768 IF (PRESENT(mu)) mu = mo_set%mu 00769 IF (PRESENT(flexible_electron_count)) flexible_electron_count = mo_set%flexible_electron_count 00770 00771 END SUBROUTINE get_mo_set 00772 00773 ! ***************************************************************************** 00783 SUBROUTINE set_mo_occupation_1(mo_set,smear,eval_deriv,xas_env, error) 00784 00785 TYPE(mo_set_type), POINTER :: mo_set 00786 TYPE(smear_type), OPTIONAL, POINTER :: smear 00787 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 00788 POINTER :: eval_deriv 00789 TYPE(xas_environment_type), OPTIONAL, 00790 POINTER :: xas_env 00791 TYPE(cp_error_type), INTENT(INOUT) :: error 00792 00793 CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mo_occupation_1', 00794 routineP = moduleN//':'//routineN 00795 00796 INTEGER :: handle, i, i_first, imo, nmo, 00797 nomo, stat, xas_estate 00798 LOGICAL :: equal_size, failure, is_large 00799 REAL(KIND=dp) :: e1, e2, edelta, edist, 00800 el_count, lengthscale, nelec, 00801 occ_estate, xas_nelectron 00802 REAL(KIND=dp), ALLOCATABLE, 00803 DIMENSION(:, :) :: dfde 00804 00805 CALL timeset(routineN,handle) 00806 00807 failure = .FALSE. 00808 00809 CPPrecondition(ASSOCIATED(mo_set),cp_failure_level,routineP,error,failure) 00810 CPPrecondition(ASSOCIATED(mo_set%eigenvalues),cp_failure_level,routineP,error,failure) 00811 CPPrecondition(ASSOCIATED(mo_set%occupation_numbers),cp_failure_level,routineP,error,failure) 00812 mo_set%occupation_numbers(:) = 0.0_dp 00813 00814 ! Quick return, if no electrons are available 00815 IF (mo_set%nelectron == 0) THEN 00816 CALL timestop(handle) 00817 RETURN 00818 END IF 00819 00820 IF (MODULO(mo_set%nelectron,INT(mo_set%maxocc)) == 0) THEN 00821 nomo = NINT(mo_set%nelectron/mo_set%maxocc) 00822 ! Initialize MO occupations 00823 mo_set%occupation_numbers(1:nomo) = mo_set%maxocc 00824 ELSE 00825 nomo = INT(mo_set%nelectron/mo_set%maxocc) + 1 00826 ! Initialize MO occupations 00827 mo_set%occupation_numbers(1:nomo-1) = mo_set%maxocc 00828 mo_set%occupation_numbers(nomo) = mo_set%nelectron -(nomo-1)* mo_set%maxocc 00829 END IF 00830 nmo = SIZE(mo_set%eigenvalues) 00831 xas_estate = -1 00832 occ_estate = 0.0_dp 00833 00834 CPPrecondition((nmo >= nomo),cp_failure_level,routineP,error,failure) 00835 CPPrecondition((SIZE(mo_set%occupation_numbers) == nmo),cp_failure_level,routineP,error,failure) 00836 00837 IF(PRESENT(xas_env) )THEN 00838 CALL get_xas_env(xas_env=xas_env,occ_estate=occ_estate,xas_estate=xas_estate,& 00839 xas_nelectron=xas_nelectron,error=error) 00840 mo_set%occupation_numbers(xas_estate) = occ_estate 00841 el_count = 0.0_dp 00842 DO i = 1,nomo 00843 el_count = el_count + mo_set%occupation_numbers(i) 00844 IF((el_count-xas_nelectron) >EPSILON(0.0_dp)) THEN 00845 nomo = i 00846 mo_set%occupation_numbers(nomo) = mo_set%occupation_numbers(nomo) - & 00847 (el_count-xas_nelectron) 00848 EXIT 00849 END IF 00850 END DO 00851 IF((xas_nelectron - el_count) > EPSILON(0.0_dp) ) THEN 00852 nomo = nomo + 1 00853 mo_set%occupation_numbers(nomo) = xas_nelectron - el_count 00854 END IF 00855 ENDIF 00856 ! zeros don't count as uniform 00857 !MK mo_set%uniform_occupation = ALL(mo_set%occupation_numbers==mo_set%maxocc) 00858 00859 mo_set%homo = nomo 00860 mo_set%lfomo = nomo + 1 00861 mo_set%mu = mo_set%eigenvalues(nomo) 00862 00863 ! Check consistency of the array lengths 00864 IF (PRESENT(eval_deriv)) THEN 00865 equal_size = (SIZE(mo_set%occupation_numbers,1) == SIZE(eval_deriv,1)) 00866 CPPrecondition(equal_size,cp_failure_level,routineP,error,failure) 00867 END IF 00868 00869 ! Quick return, if no smearing information is supplied (TO BE FIXED, smear should become non-optional...) 00870 IF (.NOT.PRESENT(smear)) THEN 00871 ! there is no dependence of the energy on the eigenvalues 00872 mo_set%uniform_occupation = .TRUE. 00873 IF (PRESENT(eval_deriv)) THEN 00874 eval_deriv = 0.0_dp 00875 END IF 00876 CALL timestop(handle) 00877 RETURN 00878 END IF 00879 00880 ! Check if proper eigenvalues are already available 00881 IF (smear%method /= smear_list) THEN 00882 IF ((ABS(mo_set%eigenvalues(1)) < 1.0E-12_dp).AND.& 00883 (ABS(mo_set%eigenvalues(nmo)) < 1.0E-12_dp)) THEN 00884 CALL timestop(handle) 00885 RETURN 00886 END IF 00887 END IF 00888 00889 ! Perform smearing 00890 IF (smear%do_smear) THEN 00891 IF(PRESENT(xas_env) ) THEN 00892 i_first = xas_estate + 1 00893 nelec = xas_nelectron 00894 ELSE 00895 i_first = 1 00896 ! nelec = REAL(mo_set%nelectron,dp) 00897 nelec = mo_set%n_el_f 00898 END IF 00899 SELECT CASE (smear%method) 00900 CASE (smear_fermi_dirac) 00901 IF (.NOT. PRESENT(eval_deriv)) THEN 00902 CALL FermiFixed(mo_set%occupation_numbers,mo_set%mu,mo_set%kTS,mo_set%eigenvalues,Nelec, & 00903 smear%electronic_temperature,mo_set%maxocc,xas_estate,occ_estate) 00904 ELSE 00905 ! could be a relatively large matrix, but one could get rid of it by never storing it 00906 ! we only need dE/df * df/de, one could equally parallelize over entries, this could become expensive 00907 ALLOCATE(dfde(nmo,nmo),stat=stat) 00908 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00909 ! lengthscale could become a parameter, but this is pretty good 00910 lengthscale=10*smear%electronic_temperature 00911 00912 CALL FermiFixedDeriv(dfde,mo_set%occupation_numbers,mo_set%mu,mo_set%kTS,mo_set%eigenvalues,Nelec, & 00913 smear%electronic_temperature,mo_set%maxocc,xas_estate,occ_estate,lengthscale) 00914 00915 ! deriv of E_{KS}-kT*S wrt to f_i 00916 eval_deriv=eval_deriv - mo_set%eigenvalues + mo_set%mu 00917 ! correspondingly the deriv of E_{KS}-kT*S wrt to e_i 00918 eval_deriv=MATMUL(TRANSPOSE(dfde),eval_deriv) 00919 00920 DEALLOCATE(dfde,stat=stat) 00921 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00922 END IF 00923 00924 ! Find the lowest fractional occupied MO (LFOMO) 00925 DO imo=i_first,nmo 00926 IF (mo_set%occupation_numbers(imo) < mo_set%maxocc) THEN 00927 mo_set%lfomo = imo 00928 EXIT 00929 END IF 00930 END DO 00931 is_large=ABS(MAXVAL(mo_set%occupation_numbers)-mo_set%maxocc)> smear%eps_fermi_dirac 00932 ! this is not a real problem, but the temperature might be a bit large 00933 CALL cp_assert(.NOT.is_large,cp_warning_level,cp_assertion_failed,routineP,& 00934 "Fermi-Dirac smearing includes the first MO",& 00935 error,failure) 00936 00937 ! Find the highest (fractional) occupied MO which will be now the HOMO 00938 DO imo=nmo,mo_set%lfomo,-1 00939 IF (mo_set%occupation_numbers(imo) > smear%eps_fermi_dirac) THEN 00940 mo_set%homo = imo 00941 EXIT 00942 END IF 00943 END DO 00944 is_large=ABS(MINVAL(mo_set%occupation_numbers))> smear%eps_fermi_dirac 00945 CALL cp_assert(.NOT.is_large,cp_warning_level,cp_assertion_failed,routineP,& 00946 "Fermi-Dirac smearing includes the last MO => "//& 00947 "Add more MOs for proper smearing.",error,failure) 00948 00949 ! check that the total electron count is accurate 00950 is_large=(ABS(nelec - accurate_sum(mo_set%occupation_numbers(:))) > smear%eps_fermi_dirac*nelec) 00951 CALL cp_assert(.NOT.is_large,& 00952 cp_warning_level,cp_assertion_failed,routineP,& 00953 "Total number of electrons is not accurate",& 00954 error,failure) 00955 00956 CASE (smear_energy_window) 00957 ! not implemented 00958 CPPrecondition(.NOT.PRESENT(eval_deriv),cp_failure_level,routineP,error,failure) 00959 00960 ! Define the energy window for the eigenvalues 00961 e1 = mo_set%eigenvalues(mo_set%homo) - 0.5_dp*smear%window_size 00962 CALL cp_assert((e1 > mo_set%eigenvalues(1)),cp_warning_level,cp_assertion_failed,routineP,& 00963 "Energy window for smearing includes the first MO",& 00964 error,failure) 00965 00966 e2 = mo_set%eigenvalues(mo_set%homo) + 0.5_dp*smear%window_size 00967 CALL cp_assert((e2 < mo_set%eigenvalues(nmo)),cp_warning_level,cp_assertion_failed,routineP,& 00968 "Energy window for smearing includes the last MO => "//& 00969 "Add more MOs for proper smearing.",error,failure) 00970 00971 ! Find the lowest fractional occupied MO (LFOMO) 00972 DO imo=i_first,nomo 00973 IF (mo_set%eigenvalues(imo) > e1) THEN 00974 mo_set%lfomo = imo 00975 EXIT 00976 END IF 00977 END DO 00978 00979 ! Find the highest fractional occupied (non-zero) MO which will be the HOMO 00980 DO imo=nmo,nomo,-1 00981 IF (mo_set%eigenvalues(imo) < e2) THEN 00982 mo_set%homo = imo 00983 EXIT 00984 END IF 00985 END DO 00986 00987 ! Get the number of electrons to be smeared 00988 edist = 0.0_dp 00989 nelec = 0.0_dp 00990 00991 DO imo=mo_set%lfomo,mo_set%homo 00992 nelec = nelec + mo_set%occupation_numbers(imo) 00993 edist = edist + ABS(e2 - mo_set%eigenvalues(imo)) 00994 END DO 00995 00996 ! Smear electrons inside the energy window 00997 DO imo=mo_set%lfomo,mo_set%homo 00998 edelta = ABS(e2 - mo_set%eigenvalues(imo)) 00999 mo_set%occupation_numbers(imo) = MIN(mo_set%maxocc,nelec*edelta/edist) 01000 nelec = nelec - mo_set%occupation_numbers(imo) 01001 edist = edist - edelta 01002 END DO 01003 01004 CASE (smear_list) 01005 equal_size = SIZE(mo_set%occupation_numbers,1)==SIZE(smear%list,1) 01006 CPPrecondition(equal_size,cp_failure_level,routineP,error,failure) 01007 mo_set%occupation_numbers = smear%list 01008 ! there is no dependence of the energy on the eigenvalues 01009 IF (PRESENT(eval_deriv)) THEN 01010 eval_deriv = 0.0_dp 01011 END IF 01012 ! most general case 01013 mo_set%lfomo=1 01014 mo_set%homo =nmo 01015 END SELECT 01016 01017 ! Check, if the smearing involves more than one MO 01018 IF (mo_set%lfomo == mo_set%homo) THEN 01019 mo_set%homo = nomo 01020 mo_set%lfomo = nomo + 1 01021 ELSE 01022 mo_set%uniform_occupation = .FALSE. 01023 END IF 01024 01025 END IF ! do smear 01026 01027 CALL timestop(handle) 01028 01029 END SUBROUTINE set_mo_occupation_1 01030 01031 ! ***************************************************************************** 01037 SUBROUTINE set_mo_set(mo_set,maxocc,homo,lfomo,nao,nelectron,n_el_f,nmo,& 01038 eigenvalues,occupation_numbers,uniform_occupation,& 01039 kTS,mu,flexible_electron_count,error) 01040 01041 TYPE(mo_set_type), POINTER :: mo_set 01042 REAL(KIND=dp), INTENT(IN), OPTIONAL :: maxocc 01043 INTEGER, INTENT(IN), OPTIONAL :: homo, lfomo, nao, nelectron 01044 REAL(KIND=dp), INTENT(IN), OPTIONAL :: n_el_f 01045 INTEGER, INTENT(IN), OPTIONAL :: nmo 01046 REAL(KIND=dp), DIMENSION(:), OPTIONAL, 01047 POINTER :: eigenvalues, 01048 occupation_numbers 01049 LOGICAL, INTENT(IN), OPTIONAL :: uniform_occupation 01050 REAL(KIND=dp), INTENT(IN), OPTIONAL :: kTS, mu, 01051 flexible_electron_count 01052 TYPE(cp_error_type), INTENT(INOUT) :: error 01053 01054 CHARACTER(LEN=*), PARAMETER :: routineN = 'set_mo_set', 01055 routineP = moduleN//':'//routineN 01056 01057 INTEGER :: istat 01058 LOGICAL :: failure 01059 01060 IF (PRESENT(maxocc)) mo_set%maxocc = maxocc 01061 IF (PRESENT(homo)) mo_set%homo = homo 01062 IF (PRESENT(lfomo)) mo_set%lfomo = lfomo 01063 IF (PRESENT(nao)) mo_set%nao = nao 01064 IF (PRESENT(nelectron)) mo_set%nelectron = nelectron 01065 IF (PRESENT(n_el_f)) mo_set%n_el_f = n_el_f 01066 IF (PRESENT(nmo)) mo_set%nmo = nmo 01067 IF (PRESENT(eigenvalues)) THEN 01068 IF (ASSOCIATED(mo_set%eigenvalues)) THEN 01069 DEALLOCATE(mo_set%eigenvalues,STAT=istat) 01070 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01071 END IF 01072 mo_set%eigenvalues => eigenvalues 01073 END IF 01074 IF (PRESENT(occupation_numbers)) THEN 01075 IF (ASSOCIATED(mo_set%occupation_numbers)) THEN 01076 DEALLOCATE(mo_set%occupation_numbers,STAT=istat) 01077 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01078 END IF 01079 mo_set%occupation_numbers => occupation_numbers 01080 END IF 01081 IF(PRESENT(uniform_occupation)) mo_set%uniform_occupation = uniform_occupation 01082 IF(PRESENT(kTS)) mo_set%kTS = kTS 01083 IF(PRESENT(mu)) mo_set%mu = mu 01084 IF(PRESENT(flexible_electron_count)) mo_set%flexible_electron_count = flexible_electron_count 01085 01086 END SUBROUTINE set_mo_set 01087 01088 ! ***************************************************************************** 01089 SUBROUTINE write_mo_set_to_restart(mo_array,particle_set,dft_section,scp,& 01090 scp_env,atomic_kind_set,error) 01091 01092 TYPE(mo_set_p_type), DIMENSION(:), 01093 POINTER :: mo_array 01094 TYPE(particle_type), DIMENSION(:), 01095 POINTER :: particle_set 01096 TYPE(section_vals_type), POINTER :: dft_section 01097 LOGICAL, INTENT(IN) :: scp 01098 TYPE(scp_environment_type), OPTIONAL, 01099 POINTER :: scp_env 01100 TYPE(atomic_kind_type), DIMENSION(:), 01101 OPTIONAL, POINTER :: atomic_kind_set 01102 TYPE(cp_error_type), INTENT(inout) :: error 01103 01104 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_mo_set_to_restart', 01105 routineP = moduleN//':'//routineN 01106 01107 CHARACTER(LEN=30), DIMENSION(2) :: 01108 keys = (/"SCF%PRINT%RESTART_HISTORY","SCF%PRINT%RESTART "/) 01109 INTEGER :: handle, ikey, ires, ispin 01110 LOGICAL :: failure 01111 TYPE(cp_logger_type), POINTER :: logger 01112 01113 CALL timeset(routineN,handle) 01114 failure = .FALSE. 01115 logger => cp_error_get_logger(error) 01116 01117 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01118 dft_section,keys(1),error=error),cp_p_file) .OR. & 01119 BTEST(cp_print_key_should_output(logger%iter_info,& 01120 dft_section,keys(2),error=error),cp_p_file) ) THEN 01121 01122 IF(mo_array(1)%mo_set%use_mo_coeff_b) THEN 01123 ! we are using the dbcsr mo_coeff 01124 ! we copy it to the fm for anycase 01125 DO ispin=1,SIZE(mo_array) 01126 IF(.not.ASSOCIATED(mo_array(ispin)%mo_set%mo_coeff_b)) THEN 01127 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 01128 ENDIF 01129 CALL copy_dbcsr_to_fm(mo_array(ispin)%mo_set%mo_coeff_b,& 01130 mo_array(ispin)%mo_set%mo_coeff,error=error)!fm->dbcsr 01131 ENDDO 01132 ENDIF 01133 01134 DO ikey=1,SIZE(keys) 01135 01136 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01137 dft_section,keys(ikey),error=error),cp_p_file)) THEN 01138 01139 ires = cp_print_key_unit_nr(logger,dft_section,keys(ikey),& 01140 extension=".wfn", file_status="REPLACE", file_action="WRITE",& 01141 do_backup=.TRUE., file_form="UNFORMATTED", error=error) 01142 01143 CALL write_mo_set_low(mo_array, particle_set=particle_set, ires=ires, error=error) 01144 01145 IF (scp) CALL write_scp_coeff_set(ires, scp_env, atomic_kind_set, particle_set, error) 01146 CALL cp_print_key_finished_output(ires,logger,dft_section,TRIM(keys(ikey)), error=error) 01147 END IF 01148 END DO 01149 END IF 01150 01151 CALL timestop(handle) 01152 01153 END SUBROUTINE write_mo_set_to_restart 01154 01155 ! ***************************************************************************** 01156 SUBROUTINE write_rt_mos_to_restart(mo_array,rt_mos,particle_set,dft_section,scp,& 01157 scp_env,atomic_kind_set,error) 01158 01159 TYPE(mo_set_p_type), DIMENSION(:), 01160 POINTER :: mo_array 01161 TYPE(cp_fm_p_type), DIMENSION(:), 01162 POINTER :: rt_mos 01163 TYPE(particle_type), DIMENSION(:), 01164 POINTER :: particle_set 01165 TYPE(section_vals_type), POINTER :: dft_section 01166 LOGICAL, INTENT(IN) :: scp 01167 TYPE(scp_environment_type), OPTIONAL, 01168 POINTER :: scp_env 01169 TYPE(atomic_kind_type), DIMENSION(:), 01170 OPTIONAL, POINTER :: atomic_kind_set 01171 TYPE(cp_error_type), INTENT(inout) :: error 01172 01173 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_rt_mos_to_restart', 01174 routineP = moduleN//':'//routineN 01175 01176 CHARACTER(LEN=43), DIMENSION(2) :: keys = (/ 01177 "REAL_TIME_PROPAGATION%PRINT%RESTART_HISTORY", 01178 "REAL_TIME_PROPAGATION%PRINT%RESTART "/) 01179 INTEGER :: handle, ikey, ires 01180 LOGICAL :: failure 01181 TYPE(cp_logger_type), POINTER :: logger 01182 01183 CALL timeset(routineN,handle) 01184 failure = .FALSE. 01185 logger => cp_error_get_logger(error) 01186 01187 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01188 dft_section,keys(1),error=error),cp_p_file) .OR. & 01189 BTEST(cp_print_key_should_output(logger%iter_info,& 01190 dft_section,keys(2),error=error),cp_p_file) ) THEN 01191 01192 DO ikey=1,SIZE(keys) 01193 01194 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01195 dft_section,keys(ikey),error=error),cp_p_file)) THEN 01196 01197 ires = cp_print_key_unit_nr(logger,dft_section,keys(ikey),& 01198 extension=".rtpwfn", file_status="REPLACE", file_action="WRITE",& 01199 do_backup=.TRUE., file_form="UNFORMATTED", error=error) 01200 01201 CALL write_mo_set_low(mo_array, rt_mos=rt_mos, particle_set=particle_set, ires=ires, error=error) 01202 IF (scp) CALL write_scp_coeff_set(ires, scp_env, atomic_kind_set, particle_set, error) 01203 CALL cp_print_key_finished_output(ires,logger,dft_section,TRIM(keys(ikey)), error=error) 01204 END IF 01205 END DO 01206 END IF 01207 01208 CALL timestop(handle) 01209 01210 END SUBROUTINE write_rt_mos_to_restart 01211 01212 ! ***************************************************************************** 01213 SUBROUTINE write_mo_set_low(mo_array, particle_set, ires, rt_mos, error) 01214 01215 TYPE(mo_set_p_type), DIMENSION(:), 01216 POINTER :: mo_array 01217 TYPE(particle_type), DIMENSION(:), 01218 POINTER :: particle_set 01219 INTEGER :: ires 01220 TYPE(cp_fm_p_type), DIMENSION(:), 01221 OPTIONAL, POINTER :: rt_mos 01222 TYPE(cp_error_type), INTENT(inout) :: error 01223 01224 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_mo_set_low', 01225 routineP = moduleN//':'//routineN 01226 01227 INTEGER :: handle, iatom, imat, iset, ishell, ispin, istat, lmax, lshell, 01228 max_block, nao, natom, nmo, nset, nset_max, nshell_max, nspin 01229 INTEGER, DIMENSION(:), POINTER :: nset_info, nshell 01230 INTEGER, DIMENSION(:, :), POINTER :: l, nshell_info 01231 INTEGER, DIMENSION(:, :, :), POINTER :: nso_info 01232 LOGICAL :: failure 01233 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 01234 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter 01235 01236 CALL timeset(routineN,handle) 01237 nspin = SIZE(mo_array) 01238 nao = mo_array(1)%mo_set%nao 01239 01240 IF (ires>0) THEN 01241 ! *** create some info about the basis set first *** 01242 natom = SIZE(particle_set,1) 01243 nset_max = 0 01244 nshell_max = 0 01245 01246 DO iatom=1,natom 01247 NULLIFY(orb_basis_set,dftb_parameter) 01248 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 01249 orb_basis_set=orb_basis_set,dftb_parameter=dftb_parameter) 01250 IF (ASSOCIATED(orb_basis_set)) THEN 01251 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01252 nset=nset,& 01253 nshell=nshell,& 01254 l=l) 01255 nset_max = MAX(nset_max,nset) 01256 DO iset=1,nset 01257 nshell_max = MAX(nshell_max,nshell(iset)) 01258 END DO 01259 ELSEIF (ASSOCIATED(dftb_parameter)) THEN 01260 CALL get_dftb_atom_param(dftb_parameter,lmax=lmax) 01261 nset_max = MAX(nset_max,1) 01262 nshell_max = MAX(nshell_max,lmax+1) 01263 ELSE 01264 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01265 "Unknown basis type. "//& 01266 CPSourceFileRef,& 01267 only_ionode=.TRUE.) 01268 END IF 01269 END DO 01270 01271 ALLOCATE (nso_info(nshell_max,nset_max,natom),STAT=istat) 01272 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01273 nso_info(:,:,:) = 0 01274 01275 ALLOCATE (nshell_info(nset_max,natom),STAT=istat) 01276 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01277 nshell_info(:,:) = 0 01278 01279 ALLOCATE (nset_info(natom),STAT=istat) 01280 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01281 nset_info(:) = 0 01282 01283 DO iatom=1,natom 01284 NULLIFY(orb_basis_set,dftb_parameter) 01285 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 01286 orb_basis_set=orb_basis_set,dftb_parameter=dftb_parameter) 01287 IF (ASSOCIATED(orb_basis_set)) THEN 01288 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01289 nset=nset,& 01290 nshell=nshell,& 01291 l=l) 01292 nset_info(iatom) = nset 01293 DO iset=1,nset 01294 nshell_info(iset,iatom)=nshell(iset) 01295 DO ishell=1,nshell(iset) 01296 lshell = l(ishell,iset) 01297 nso_info(ishell,iset,iatom) = nso(lshell) 01298 END DO 01299 END DO 01300 ELSEIF (ASSOCIATED(dftb_parameter)) THEN 01301 CALL get_dftb_atom_param(dftb_parameter,lmax=lmax) 01302 nset_info(iatom) = 1 01303 nshell_info(1,iatom)=lmax+1 01304 DO ishell=1,lmax+1 01305 lshell = ishell-1 01306 nso_info(ishell,1,iatom) = nso(lshell) 01307 END DO 01308 ELSE 01309 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01310 "Unknown basis set type. "//& 01311 CPSourceFileRef,& 01312 only_ionode=.TRUE.) 01313 END IF 01314 END DO 01315 01316 WRITE (ires) natom,nspin,nao,nset_max,nshell_max 01317 WRITE (ires) nset_info 01318 WRITE (ires) nshell_info 01319 WRITE (ires) nso_info 01320 01321 DEALLOCATE (nset_info,STAT=istat) 01322 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01323 01324 DEALLOCATE (nshell_info,STAT=istat) 01325 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01326 01327 DEALLOCATE (nso_info,STAT=istat) 01328 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01329 END IF 01330 01331 ! use the scalapack block size as a default for buffering columns 01332 CALL cp_fm_get_info(mo_array(1)%mo_set%mo_coeff,ncol_block=max_block,error=error) 01333 01334 01335 DO ispin=1,nspin 01336 nmo=mo_array(ispin)%mo_set%nmo 01337 IF ((ires>0).AND.(nmo > 0)) THEN 01338 WRITE (ires) nmo,& 01339 mo_array(ispin)%mo_set%homo,& 01340 mo_array(ispin)%mo_set%lfomo,& 01341 mo_array(ispin)%mo_set%nelectron 01342 WRITE (ires) mo_array(ispin)%mo_set%eigenvalues(1:nmo),& 01343 mo_array(ispin)%mo_set%occupation_numbers(1:nmo) 01344 END IF 01345 IF(PRESENT(rt_mos))THEN 01346 DO imat=2*ispin-1,2*ispin 01347 CALL cp_fm_write_unformatted(rt_mos(imat)%matrix,ires,error) 01348 END DO 01349 ELSE 01350 CALL cp_fm_write_unformatted(mo_array(ispin)%mo_set%mo_coeff,ires,error) 01351 END IF 01352 END DO 01353 01354 CALL timestop(handle) 01355 01356 END SUBROUTINE write_mo_set_low 01357 01358 ! ***************************************************************************** 01359 SUBROUTINE wfn_restart_file_name(filename,exist,section,logger,xas,rtp,error) 01360 CHARACTER(LEN=default_path_length), 01361 INTENT(OUT) :: filename 01362 LOGICAL, INTENT(OUT) :: exist 01363 TYPE(section_vals_type), POINTER :: section 01364 TYPE(cp_logger_type), POINTER :: logger 01365 LOGICAL, INTENT(IN), OPTIONAL :: xas, rtp 01366 TYPE(cp_error_type), INTENT(inout) :: error 01367 01368 INTEGER :: n_rep_val 01369 LOGICAL :: my_rtp, my_xas 01370 TYPE(section_vals_type), POINTER :: print_key 01371 01372 my_xas = .FALSE. 01373 my_rtp = .FALSE. 01374 IF(PRESENT(xas)) my_xas = xas 01375 IF(PRESENT(rtp)) my_rtp = rtp 01376 01377 exist = .FALSE. 01378 CALL section_vals_val_get(section,"WFN_RESTART_FILE_NAME",n_rep_val=n_rep_val,error=error) 01379 IF (n_rep_val>0) THEN 01380 CALL section_vals_val_get(section,"WFN_RESTART_FILE_NAME",c_val=filename,error=error) 01381 ELSE 01382 IF(my_xas) THEN 01383 ! try to read from the filename that is generated automatically from the printkey 01384 print_key => section_vals_get_subs_vals(section,"PRINT%RESTART",error=error) 01385 filename = cp_print_key_generate_filename(logger,print_key, & 01386 extension="",my_local=.FALSE., error=error) 01387 ELSE IF (my_rtp)THEN 01388 ! try to read from the filename that is generated automatically from the printkey 01389 print_key => section_vals_get_subs_vals(section,"REAL_TIME_PROPAGATION%PRINT%RESTART",error=error) 01390 filename = cp_print_key_generate_filename(logger,print_key, & 01391 extension=".rtpwfn",my_local=.FALSE., error=error) 01392 ELSE 01393 ! try to read from the filename that is generated automatically from the printkey 01394 print_key => section_vals_get_subs_vals(section,"SCF%PRINT%RESTART",error=error) 01395 filename = cp_print_key_generate_filename(logger,print_key, & 01396 extension=".wfn", my_local=.FALSE., error=error) 01397 END IF 01398 ENDIF 01399 IF(.NOT.my_xas) THEN 01400 INQUIRE(FILE=filename,exist=exist) 01401 END IF 01402 01403 END SUBROUTINE wfn_restart_file_name 01404 01405 ! ***************************************************************************** 01406 SUBROUTINE read_mo_set_from_restart(mo_array,atomic_kind_set,particle_set,& 01407 para_env,id_nr,multiplicity,dft_section,scp,scp_env,natom_mismatch,error) 01408 01409 TYPE(mo_set_p_type), DIMENSION(:), 01410 POINTER :: mo_array 01411 TYPE(atomic_kind_type), DIMENSION(:), 01412 POINTER :: atomic_kind_set 01413 TYPE(particle_type), DIMENSION(:), 01414 POINTER :: particle_set 01415 TYPE(cp_para_env_type), POINTER :: para_env 01416 INTEGER, INTENT(IN) :: id_nr, multiplicity 01417 TYPE(section_vals_type), POINTER :: dft_section 01418 LOGICAL, INTENT(IN) :: scp 01419 TYPE(scp_environment_type), POINTER :: scp_env 01420 LOGICAL, INTENT(OUT), OPTIONAL :: natom_mismatch 01421 TYPE(cp_error_type), INTENT(inout) :: error 01422 01423 CHARACTER(LEN=*), PARAMETER :: routineN = 'read_mo_set_from_restart', 01424 routineP = moduleN//':'//routineN 01425 01426 CHARACTER(LEN=default_path_length) :: file_name 01427 INTEGER :: group, handle, ispin, natom, 01428 nspin, restart_unit, source 01429 LOGICAL :: exist, failure 01430 TYPE(cp_logger_type), POINTER :: logger 01431 01432 CALL timeset(routineN,handle) 01433 logger => cp_error_get_logger(error) 01434 failure = .FALSE. 01435 01436 nspin = SIZE(mo_array) 01437 restart_unit = -1 01438 01439 group = para_env%group 01440 source = para_env%source 01441 01442 IF (para_env%ionode) THEN 01443 01444 natom = SIZE(particle_set,1) 01445 CALL wfn_restart_file_name(file_name,exist,dft_section,logger,error=error) 01446 IF (id_nr/=0) THEN 01447 ! Is it one of the backup files? 01448 file_name = TRIM(file_name)//".bak-"//ADJUSTL(cp_to_string(id_nr)) 01449 END IF 01450 01451 CALL open_file(file_name=file_name,& 01452 file_action="READ",& 01453 file_form="UNFORMATTED",& 01454 file_status="OLD",& 01455 unit_number=restart_unit) 01456 01457 END IF 01458 01459 CALL read_mos_restart_low (mo_array,para_env=para_env, particle_set=particle_set,& 01460 natom=natom,& 01461 rst_unit=restart_unit, multiplicity=multiplicity, natom_mismatch=natom_mismatch, & 01462 error=error) 01463 IF (PRESENT(natom_mismatch)) THEN 01464 ! read_mos_restart_low only the io_node returns natom_mismatch, must broadcast it 01465 CALL mp_bcast(natom_mismatch,source,group) 01466 IF (natom_mismatch) THEN 01467 IF (para_env%ionode) CALL close_file(unit_number=restart_unit) 01468 CALL timestop(handle) 01469 RETURN 01470 ENDIF 01471 ENDIF 01472 01473 IF (scp) CALL read_aux_coeff_set(restart_unit, scp_env%aux_coeff_set, para_env, error) 01474 01475 ! Close restart file 01476 IF (para_env%ionode) CALL close_file(unit_number=restart_unit) 01477 01478 DO ispin = 1,nspin 01479 CALL write_mo_set(mo_array(ispin)%mo_set,atomic_kind_set,particle_set,4,& 01480 dft_section,error=error) 01481 END DO 01482 01483 CALL timestop(handle) 01484 01485 END SUBROUTINE read_mo_set_from_restart 01486 ! ***************************************************************************** 01487 SUBROUTINE read_rt_mos_from_restart(mo_array,rt_mos,atomic_kind_set,particle_set,& 01488 para_env,id_nr,multiplicity,dft_section,error) 01489 01490 TYPE(mo_set_p_type), DIMENSION(:), 01491 POINTER :: mo_array 01492 TYPE(cp_fm_p_type), DIMENSION(:), 01493 POINTER :: rt_mos 01494 TYPE(atomic_kind_type), DIMENSION(:), 01495 POINTER :: atomic_kind_set 01496 TYPE(particle_type), DIMENSION(:), 01497 POINTER :: particle_set 01498 TYPE(cp_para_env_type), POINTER :: para_env 01499 INTEGER, INTENT(IN) :: id_nr, multiplicity 01500 TYPE(section_vals_type), POINTER :: dft_section 01501 TYPE(cp_error_type), INTENT(inout) :: error 01502 01503 CHARACTER(LEN=*), PARAMETER :: routineN = 'read_rt_mos_from_restart', 01504 routineP = moduleN//':'//routineN 01505 01506 CHARACTER(LEN=default_path_length) :: file_name 01507 INTEGER :: group, handle, ispin, natom, 01508 nspin, restart_unit, source 01509 LOGICAL :: exist, failure 01510 TYPE(cp_logger_type), POINTER :: logger 01511 01512 CALL timeset(routineN,handle) 01513 logger => cp_error_get_logger(error) 01514 failure = .FALSE. 01515 01516 nspin = SIZE(mo_array) 01517 restart_unit = -1 01518 01519 group = para_env%group 01520 source = para_env%source 01521 01522 IF (para_env%ionode) THEN 01523 01524 natom = SIZE(particle_set,1) 01525 CALL wfn_restart_file_name(file_name,exist,dft_section,logger,rtp=.TRUE.,error=error) 01526 IF (id_nr/=0) THEN 01527 ! Is it one of the backup files? 01528 file_name = TRIM(file_name)//".bak-"//ADJUSTL(cp_to_string(id_nr)) 01529 END IF 01530 01531 CALL open_file(file_name=file_name,& 01532 file_action="READ",& 01533 file_form="UNFORMATTED",& 01534 file_status="OLD",& 01535 unit_number=restart_unit) 01536 01537 END IF 01538 01539 CALL read_mos_restart_low (mo_array, rt_mos=rt_mos,para_env= para_env,& 01540 particle_set=particle_set, natom=natom,& 01541 rst_unit=restart_unit, multiplicity=multiplicity, error=error) 01542 01543 ! Close restart file 01544 IF (para_env%ionode) CALL close_file(unit_number=restart_unit) 01545 01546 DO ispin = 1,nspin 01547 CALL write_mo_set(mo_array(ispin)%mo_set,atomic_kind_set,particle_set,4,& 01548 dft_section,error=error) 01549 END DO 01550 01551 CALL timestop(handle) 01552 01553 END SUBROUTINE read_rt_mos_from_restart 01554 ! ***************************************************************************** 01560 SUBROUTINE read_mos_restart_low (mos, para_env, particle_set, natom, rst_unit, & 01561 multiplicity, rt_mos, natom_mismatch, error) 01562 01563 TYPE(mo_set_p_type), DIMENSION(:), 01564 POINTER :: mos 01565 TYPE(cp_para_env_type), POINTER :: para_env 01566 TYPE(particle_type), DIMENSION(:), 01567 POINTER :: particle_set 01568 INTEGER, INTENT(IN) :: natom, rst_unit 01569 INTEGER, INTENT(in), OPTIONAL :: multiplicity 01570 TYPE(cp_fm_p_type), DIMENSION(:), 01571 OPTIONAL, POINTER :: rt_mos 01572 LOGICAL, INTENT(OUT), OPTIONAL :: natom_mismatch 01573 TYPE(cp_error_type), INTENT(inout) :: error 01574 01575 CHARACTER(LEN=*), PARAMETER :: routineN = 'read_mos_restart_low', 01576 routineP = moduleN//':'//routineN 01577 01578 INTEGER :: group, homo, homo_read, i, iatom, imat, irow, iset, iset_read, 01579 ishell, ishell_read, iso, ispin, istat, lfomo_read, lmax, lshell, 01580 my_mult, nao, nao_read, natom_read, nelectron, nelectron_read, nmo, 01581 nmo_read, nnshell, nset, nset_max, nshell_max, nspin, nspin_read, 01582 offset_read, source 01583 INTEGER, DIMENSION(:), POINTER :: nset_info, nshell 01584 INTEGER, DIMENSION(:, :), POINTER :: l, nshell_info 01585 INTEGER, DIMENSION(:, :, :), POINTER :: nso_info, offset_info 01586 LOGICAL :: failure, minbas, natom_match, 01587 use_this 01588 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig_read, occ_read 01589 REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer, vecbuffer_read 01590 TYPE(cp_logger_type), POINTER :: logger 01591 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 01592 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter 01593 01594 logger => cp_error_get_logger(error) 01595 01596 nspin = SIZE(mos) 01597 nao = mos(1)%mo_set%nao 01598 my_mult = 0 01599 IF(PRESENT(multiplicity)) my_mult = multiplicity 01600 group = para_env%group 01601 source = para_env%source 01602 01603 IF (para_env%ionode) THEN 01604 READ (rst_unit) natom_read,nspin_read,nao_read,nset_max,nshell_max 01605 IF(PRESENT(rt_mos)) THEN 01606 IF (nspin_read /= nspin) THEN 01607 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01608 "To change nspin is not possible. "//& 01609 CPSourceFileRef,& 01610 only_ionode=.TRUE.) 01611 END IF 01612 ELSE 01613 ! we should allow for restarting with different spin settings 01614 IF (nspin_read /= nspin) THEN 01615 WRITE(cp_logger_get_default_unit_nr(logger),*) & 01616 "READ RESTART : WARNING : nspin is not equal " 01617 END IF 01618 ! this case needs fixing of homo/lfomo/nelec/occupations ... 01619 IF (nspin_read > nspin) THEN 01620 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01621 "Reducing nspin is not possible. "//& 01622 CPSourceFileRef,& 01623 only_ionode=.TRUE.) 01624 ENDIF 01625 END IF 01626 01627 natom_match = (natom_read == natom) 01628 01629 IF (natom_match) THEN ! actually do the read read 01630 01631 ! Let's make it possible to change the basis set 01632 ALLOCATE (nso_info(nshell_max,nset_max,natom_read),STAT=istat) 01633 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01634 ALLOCATE (nshell_info(nset_max,natom_read),STAT=istat) 01635 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01636 ALLOCATE (nset_info(natom_read),STAT=istat) 01637 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01638 ALLOCATE (offset_info(nshell_max,nset_max,natom_read),STAT=istat) 01639 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01640 01641 IF (nao_read /= nao) THEN 01642 WRITE(cp_logger_get_default_unit_nr(logger),*) & 01643 " READ RESTART : WARNING : DIFFERENT # AOs ",nao,nao_read 01644 IF(PRESENT(rt_mos))& 01645 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01646 "To change basis is not possible. "//& 01647 CPSourceFileRef,& 01648 only_ionode=.TRUE.) 01649 END IF 01650 01651 READ (rst_unit) nset_info 01652 READ (rst_unit) nshell_info 01653 READ (rst_unit) nso_info 01654 01655 i=1 01656 DO iatom=1,natom 01657 DO iset=1,nset_info(iatom) 01658 DO ishell=1,nshell_info(iset,iatom) 01659 offset_info(ishell,iset,iatom) = i 01660 i=i+nso_info(ishell,iset,iatom) 01661 END DO 01662 END DO 01663 END DO 01664 01665 ALLOCATE(vecbuffer_read(1,nao_read),STAT=istat) 01666 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01667 01668 END IF ! natom_match 01669 END IF ! ionode 01670 01671 ! make natom_match and natom_mismatch uniform across all nodes 01672 CALL mp_bcast(natom_match,source,group) 01673 IF (PRESENT(natom_mismatch)) natom_mismatch = .NOT. natom_match 01674 ! handle natom_match false 01675 IF (.NOT. natom_match) THEN 01676 IF (PRESENT(natom_mismatch)) THEN 01677 WRITE(cp_logger_get_default_unit_nr(logger),*) & 01678 " READ RESTART : WARNING : DIFFERENT natom, returning ",natom,natom_read 01679 RETURN 01680 ELSE 01681 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01682 "Incorrect number of atoms in restart file. "//& 01683 CPSourceFileRef,& 01684 only_ionode=.TRUE.) 01685 ENDIF 01686 ENDIF 01687 01688 CALL mp_bcast(nspin_read,source,group) 01689 01690 ALLOCATE (vecbuffer(1,nao),STAT=istat) 01691 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01692 01693 DO ispin=1,nspin 01694 01695 nmo=mos(ispin)%mo_set%nmo 01696 homo=mos(ispin)%mo_set%homo 01697 mos(ispin)%mo_set%eigenvalues(:) = 0.0_dp 01698 mos(ispin)%mo_set%occupation_numbers(:) = 0.0_dp 01699 CALL cp_fm_set_all(mos(ispin)%mo_set%mo_coeff,0.0_dp,error=error) 01700 01701 IF (para_env%ionode.AND.(nmo > 0)) THEN 01702 READ (rst_unit) nmo_read, homo_read, lfomo_read, nelectron_read 01703 ALLOCATE(eig_read(nmo_read), occ_read(nmo_read), STAT=istat) 01704 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01705 eig_read = 0.0_dp 01706 occ_read = 0.0_dp 01707 01708 nmo = MIN(nmo,nmo_read) 01709 CALL cp_assert((nmo_read >= nmo),cp_warning_level,cp_assertion_failed,routineP,& 01710 "The number of MOs on the restart unit is smaller than the number of "//& 01711 "the allocated MOs. The MO set will be padded with zeros!"//& 01712 CPSourceFileRef,& 01713 only_ionode=.TRUE.) 01714 CALL cp_assert((nmo_read<=nmo),cp_warning_level,cp_assertion_failed,routineP,& 01715 "The number of MOs on the restart unit is greater than the number of "//& 01716 "the allocated MOs. The read MO set will be truncated!"//& 01717 CPSourceFileRef,& 01718 only_ionode=.TRUE.) 01719 01720 READ (rst_unit) eig_read(1:nmo_read), occ_read(1:nmo_read) 01721 mos(ispin)%mo_set%eigenvalues(1:nmo) = eig_read(1:nmo) 01722 mos(ispin)%mo_set%occupation_numbers(1:nmo) = occ_read(1:nmo) 01723 DEALLOCATE(eig_read, occ_read, STAT=istat) 01724 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01725 01726 mos(ispin)%mo_set%homo = homo_read 01727 mos(ispin)%mo_set%lfomo = lfomo_read 01728 IF (homo_read > nmo) THEN 01729 IF(nelectron_read==mos(ispin)%mo_set%nelectron) THEN 01730 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 01731 "The number of occupied MOs on the restart unit is larger than "//& 01732 "the allocated MOs. The read MO set will be truncated and the occupation numbers recalculated!"//& 01733 CPSourceFileRef,& 01734 only_ionode=.TRUE.) 01735 CALL set_mo_occupation(mo_set=mos(ispin)%mo_set,error=error) 01736 ELSE 01737 ! can not make this a warning i.e. homo must be smaller than nmo 01738 ! otherwise e.g. set_mo_occupation will go out of bounds 01739 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01740 "Number of occupied MOs on restart unit larger than allocated MOs. "//& 01741 CPSourceFileRef,& 01742 only_ionode=.TRUE.) 01743 END IF 01744 END IF 01745 END IF 01746 01747 CALL mp_bcast(nmo,source,group) 01748 CALL mp_bcast(mos(ispin)%mo_set%homo,source,group) 01749 CALL mp_bcast(mos(ispin)%mo_set%lfomo,source,group) 01750 CALL mp_bcast(mos(ispin)%mo_set%nelectron,source,group) 01751 CALL mp_bcast(mos(ispin)%mo_set%eigenvalues,source,group) 01752 CALL mp_bcast(mos(ispin)%mo_set%occupation_numbers,source,group) 01753 IF(PRESENT(rt_mos))THEN 01754 DO imat=2*ispin-1,2*ispin 01755 DO i=1,nmo 01756 IF (para_env%ionode) THEN 01757 READ (rst_unit) vecbuffer 01758 ELSE 01759 vecbuffer(1,:) = 0.0_dp 01760 END IF 01761 CALL mp_bcast(vecbuffer,source,group) 01762 CALL cp_fm_set_submatrix(rt_mos(imat)%matrix,& 01763 vecbuffer,1,i,nao,1,transpose=.TRUE.,error=error) 01764 END DO 01765 END DO 01766 ELSE 01767 DO i=1,nmo 01768 IF (para_env%ionode) THEN 01769 READ (rst_unit) vecbuffer_read 01770 ! now, try to assign the read to the real vector 01771 ! in case the basis set changed this involves some guessing 01772 irow=1 01773 DO iatom=1,natom 01774 NULLIFY(orb_basis_set,dftb_parameter,l,nshell) 01775 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 01776 orb_basis_set=orb_basis_set,dftb_parameter=dftb_parameter) 01777 IF (ASSOCIATED(orb_basis_set)) THEN 01778 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 01779 nset=nset,& 01780 nshell=nshell,& 01781 l=l) 01782 minbas=.FALSE. 01783 ELSEIF (ASSOCIATED(dftb_parameter)) THEN 01784 CALL get_dftb_atom_param(dftb_parameter,lmax=lmax) 01785 nset = 1 01786 minbas=.TRUE. 01787 ELSE 01788 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01789 "Unknown basis set type. "//& 01790 CPSourceFileRef,& 01791 only_ionode=.TRUE.) 01792 END IF 01793 01794 use_this = .TRUE. 01795 iset_read = 1 01796 DO iset=1,nset 01797 ishell_read = 1 01798 IF(minbas) THEN 01799 nnshell = lmax+1 01800 ELSE 01801 nnshell = nshell(iset) 01802 END IF 01803 DO ishell=1,nnshell 01804 IF(minbas) THEN 01805 lshell = ishell-1 01806 ELSE 01807 lshell = l(ishell,iset) 01808 END IF 01809 IF (iset_read > nset_info(iatom)) use_this = .FALSE. 01810 IF (use_this) THEN ! avoids out of bound access of the lower line if false 01811 IF (nso(lshell) == nso_info(ishell_read,iset_read,iatom)) THEN 01812 offset_read=offset_info(ishell_read,iset_read,iatom) 01813 ishell_read=ishell_read+1 01814 IF (ishell_read > nshell_info(iset,iatom)) THEN 01815 ishell_read = 1 01816 iset_read = iset_read+1 01817 END IF 01818 ELSE 01819 use_this = .FALSE. 01820 END IF 01821 END IF 01822 DO iso=1,nso(lshell) 01823 IF (use_this) THEN 01824 IF (offset_read-1+iso.LT.1 .OR. offset_read-1+iso.GT.nao_read) THEN 01825 vecbuffer(1,irow)=0.0_dp 01826 ELSE 01827 vecbuffer(1,irow)=vecbuffer_read(1,offset_read-1+iso) 01828 END IF 01829 ELSE 01830 vecbuffer(1,irow) = 0.0_dp 01831 END IF 01832 irow = irow + 1 01833 END DO 01834 use_this = .TRUE. 01835 END DO 01836 END DO 01837 END DO 01838 01839 ELSE 01840 01841 vecbuffer(1,:) = 0.0_dp 01842 01843 END IF 01844 01845 CALL mp_bcast(vecbuffer,source,group) 01846 CALL cp_fm_set_submatrix(mos(ispin)%mo_set%mo_coeff,& 01847 vecbuffer,1,i,nao,1,transpose=.TRUE.,error=error) 01848 END DO 01849 END IF 01850 ! Skip extra MOs if there any 01851 IF (para_env%ionode) THEN 01852 !ignore nmo = 0 01853 IF(nmo>0) THEN 01854 DO i=nmo+1,nmo_read 01855 READ (rst_unit) vecbuffer_read 01856 END DO 01857 END IF 01858 END IF 01859 01860 IF(.NOT.PRESENT(rt_mos)) THEN 01861 IF(ispin==1 .AND. nspin_read<nspin) THEN 01862 01863 mos(ispin+1)%mo_set%homo = mos(ispin)%mo_set%homo 01864 mos(ispin+1)%mo_set%lfomo = mos(ispin)%mo_set%lfomo 01865 nelectron = mos(ispin)%mo_set%nelectron 01866 IF (my_mult.NE.1) THEN 01867 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01868 "Restarting an LSD calculation from an LDA wfn only works for multiplicity=1 (singlets)."//& 01869 CPSourceFileRef,& 01870 only_ionode=.TRUE.) 01871 END IF 01872 IF(mos(ispin+1)%mo_set%nelectron < 0)THEN 01873 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 01874 "LSD: too few electrons for this multiplisity. "//& 01875 CPSourceFileRef,& 01876 only_ionode=.TRUE.) 01877 END IF 01878 mos(ispin+1)%mo_set%eigenvalues = mos(ispin)%mo_set%eigenvalues 01879 mos(ispin)%mo_set%occupation_numbers = mos(ispin)%mo_set%occupation_numbers/2.0_dp 01880 mos(ispin+1)%mo_set%occupation_numbers = mos(ispin)%mo_set%occupation_numbers 01881 CALL cp_fm_to_fm(mos(ispin)%mo_set%mo_coeff,mos(ispin+1)%mo_set%mo_coeff,error=error) 01882 EXIT 01883 END IF 01884 END IF 01885 END DO ! ispin 01886 01887 DEALLOCATE(vecbuffer,STAT=istat) 01888 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01889 01890 IF (para_env%ionode) THEN 01891 DEALLOCATE(vecbuffer_read,STAT=istat) 01892 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01893 DEALLOCATE(offset_info,STAT=istat) 01894 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01895 DEALLOCATE(nso_info,STAT=istat) 01896 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01897 DEALLOCATE(nshell_info,STAT=istat) 01898 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01899 DEALLOCATE(nset_info,STAT=istat) 01900 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01901 END IF 01902 01903 END SUBROUTINE read_mos_restart_low 01904 01905 ! ***************************************************************************** 01917 SUBROUTINE write_mo_set_to_output_unit(mo_set,atomic_kind_set,particle_set,& 01918 before,dft_section,spin,last,error) 01919 01920 TYPE(mo_set_type), POINTER :: mo_set 01921 TYPE(atomic_kind_type), DIMENSION(:), 01922 POINTER :: atomic_kind_set 01923 TYPE(particle_type), DIMENSION(:), 01924 POINTER :: particle_set 01925 INTEGER, INTENT(IN) :: before 01926 TYPE(section_vals_type), POINTER :: dft_section 01927 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: spin 01928 LOGICAL, INTENT(IN), OPTIONAL :: last 01929 TYPE(cp_error_type), INTENT(inout) :: error 01930 01931 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_mo_set_to_output_unit', 01932 routineP = moduleN//':'//routineN 01933 01934 CHARACTER(LEN=12) :: symbol 01935 CHARACTER(LEN=12), DIMENSION(:), POINTER :: bcgf_symbol 01936 CHARACTER(LEN=16) :: fmtstr5, fmtstr6 01937 CHARACTER(LEN=2) :: element_symbol 01938 CHARACTER(LEN=2*default_string_length) :: name 01939 CHARACTER(LEN=22) :: fmtstr2 01940 CHARACTER(LEN=25) :: fmtstr1, fmtstr7 01941 CHARACTER(LEN=27) :: fmtstr4 01942 CHARACTER(LEN=38) :: fmtstr3 01943 CHARACTER(LEN=6), DIMENSION(:), POINTER :: bsgf_symbol 01944 INTEGER :: after, first_mo, from, iatom, icgf, ico, icol, imo, irow, 01945 iset, isgf, ishell, iso, istat, iw, jcol, last_mo, left, lmax, lshell, 01946 natom, ncgf, ncol, ncol_global, nrow_global, nset, nsgf, right, 01947 scf_step, to, width 01948 INTEGER, DIMENSION(:), POINTER :: mo_index_range, nshell 01949 INTEGER, DIMENSION(:, :), POINTER :: l 01950 LOGICAL :: failure, ionode, my_last, 01951 p_cart, p_eval, p_evec, 01952 p_occ, should_output 01953 REAL(KIND=dp) :: gap 01954 REAL(KIND=dp), DIMENSION(:, :), POINTER :: cmatrix, smatrix 01955 TYPE(cp_logger_type), POINTER :: logger 01956 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 01957 TYPE(qs_dftb_atom_type), POINTER :: dftb_parameter 01958 01959 NULLIFY (bcgf_symbol) 01960 NULLIFY (bsgf_symbol) 01961 NULLIFY (logger) 01962 NULLIFY (mo_index_range) 01963 NULLIFY (nshell) 01964 01965 logger => cp_error_get_logger(error) 01966 ionode = logger%para_env%mepos==logger%para_env%source 01967 failure = .FALSE. 01968 CALL section_vals_val_get(dft_section,"PRINT%MO%EIGENVALUES",l_val=p_eval,error=error) 01969 CALL section_vals_val_get(dft_section,"PRINT%MO%EIGENVECTORS",l_val=p_evec,error=error) 01970 CALL section_vals_val_get(dft_section,"PRINT%MO%OCCUPATION_NUMBERS",l_val=p_occ,error=error) 01971 CALL section_vals_val_get(dft_section,"PRINT%MO%CARTESIAN",l_val=p_cart,error=error) 01972 CALL section_vals_val_get(dft_section,"PRINT%MO%MO_INDEX_RANGE",i_vals=mo_index_range,error=error) 01973 CALL section_vals_val_get(dft_section,"PRINT%MO%NDIGITS",i_val=after,error=error) 01974 after = MIN(MAX(after,1),16) 01975 should_output = BTEST(cp_print_key_should_output(logger%iter_info,dft_section,& 01976 "PRINT%MO",error=error),cp_p_file) 01977 01978 IF ((.NOT.should_output).OR.(.NOT.(p_eval.OR.p_evec.OR.p_occ))) RETURN 01979 01980 IF (PRESENT(last)) THEN 01981 my_last = last 01982 ELSE 01983 my_last = .FALSE. 01984 END IF 01985 01986 scf_step = logger%iter_info%iteration(logger%iter_info%n_rlevel) - 1 01987 01988 IF (p_evec) THEN 01989 CALL cp_fm_get_info(mo_set%mo_coeff,& 01990 nrow_global=nrow_global,& 01991 ncol_global=ncol_global,error=error) 01992 ALLOCATE(smatrix(nrow_global,ncol_global),STAT=istat) 01993 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01994 CALL cp_fm_get_submatrix(mo_set%mo_coeff,smatrix,error=error) 01995 IF (.NOT.ionode) THEN 01996 DEALLOCATE(smatrix,STAT=istat) 01997 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01998 END IF 01999 END IF 02000 02001 iw = cp_print_key_unit_nr(logger,dft_section,"PRINT%MO",& 02002 ignore_should_output=should_output,& 02003 extension=".MOLog",error=error) 02004 02005 IF (iw > 0) THEN 02006 02007 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 02008 natom=natom,& 02009 ncgf=ncgf,& 02010 nsgf=nsgf) 02011 02012 ! Definition of the variable formats 02013 02014 fmtstr1 = "(/,T2,21X, ( X,I5, X))" 02015 fmtstr2 = "(T2,21X, (1X,F . ))" 02016 fmtstr3 = "(T2,I5,1X,I5,1X,A,1X,A6, (1X,F . ))" 02017 02018 width = before + after + 3 02019 ncol = INT(56/width) 02020 02021 right = MAX((after-2),1) 02022 left = width - right - 5 02023 02024 WRITE (UNIT=fmtstr1(11:12),FMT="(I2)") ncol 02025 WRITE (UNIT=fmtstr1(14:15),FMT="(I2)") left 02026 WRITE (UNIT=fmtstr1(21:22),FMT="(I2)") right 02027 02028 WRITE (UNIT=fmtstr2(9:10),FMT="(I2)") ncol 02029 WRITE (UNIT=fmtstr2(16:17),FMT="(I2)") width - 1 02030 WRITE (UNIT=fmtstr2(19:20),FMT="(I2)") after 02031 02032 WRITE (UNIT=fmtstr3(25:26),FMT="(I2)") ncol 02033 WRITE (UNIT=fmtstr3(32:33),FMT="(I2)") width - 1 02034 WRITE (UNIT=fmtstr3(35:36),FMT="(I2)") after 02035 02036 IF (p_evec) THEN 02037 02038 IF (p_cart) THEN 02039 02040 IF (my_last) THEN 02041 name = "MO EIGENVALUES, MO OCCUPATION NUMBERS, AND "//& 02042 "CARTESIAN MO EIGENVECTORS" 02043 ELSE 02044 WRITE (UNIT=name,FMT="(A,I6)")& 02045 "MO EIGENVALUES, MO OCCUPATION NUMBERS, AND "//& 02046 "CARTESIAN MO EIGENVECTORS AFTER SCF STEP",scf_step 02047 END IF 02048 02049 ALLOCATE(cmatrix(ncgf,ncgf),STAT=istat) 02050 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02051 02052 cmatrix = 0.0_dp 02053 02054 ! Transform spherical MOs to Cartesian MOs 02055 02056 icgf = 1 02057 isgf = 1 02058 DO iatom=1,natom 02059 NULLIFY (orb_basis_set,dftb_parameter) 02060 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 02061 orb_basis_set=orb_basis_set,& 02062 dftb_parameter=dftb_parameter) 02063 IF (ASSOCIATED(orb_basis_set)) THEN 02064 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 02065 nset=nset,& 02066 nshell=nshell,& 02067 l=l) 02068 DO iset=1,nset 02069 DO ishell=1,nshell(iset) 02070 lshell = l(ishell,iset) 02071 CALL dgemm("T","N",nco(lshell),nsgf,nso(lshell),1.0_dp,& 02072 orbtramat(lshell)%s2c,nso(lshell),& 02073 smatrix(isgf,1),nsgf,0.0_dp,& 02074 cmatrix(icgf,1),ncgf) 02075 icgf = icgf + nco(lshell) 02076 isgf = isgf + nso(lshell) 02077 END DO 02078 END DO 02079 ELSE IF (ASSOCIATED(dftb_parameter)) THEN 02080 CALL get_dftb_atom_param(dftb_parameter,lmax=lmax) 02081 DO ishell=1,lmax+1 02082 lshell = ishell-1 02083 CALL dgemm("T","N",nco(lshell),nsgf,nso(lshell),1.0_dp,& 02084 orbtramat(lshell)%s2c,nso(lshell),& 02085 smatrix(isgf,1),nsgf,0.0_dp,& 02086 cmatrix(icgf,1),ncgf) 02087 icgf = icgf + nco(lshell) 02088 isgf = isgf + nso(lshell) 02089 END DO 02090 ELSE 02091 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 02092 "Unknown basis set type. "//& 02093 CPSourceFileRef,& 02094 only_ionode=.TRUE.) 02095 END IF 02096 END DO ! iatom 02097 02098 ELSE 02099 02100 IF (my_last) THEN 02101 name = "MO EIGENVALUES, MO OCCUPATION NUMBERS, AND "//& 02102 "SPHERICAL MO EIGENVECTORS" 02103 ELSE 02104 WRITE (UNIT=name,FMT="(A,I6)")& 02105 "MO EIGENVALUES, MO OCCUPATION NUMBERS, AND "//& 02106 "SPHERICAL MO EIGENVECTORS AFTER SCF STEP",scf_step 02107 END IF 02108 02109 END IF ! p_cart 02110 02111 ELSE IF (p_occ.OR.p_eval) THEN 02112 02113 IF (my_last) THEN 02114 name = "MO EIGENVALUES AND MO OCCUPATION NUMBERS" 02115 ELSE 02116 WRITE (UNIT=name,FMT="(A,I6)")& 02117 "MO EIGENVALUES AND MO OCCUPATION NUMBERS AFTER "//& 02118 "SCF STEP",scf_step 02119 END IF 02120 02121 END IF ! p_evec 02122 02123 CALL compress(name) 02124 02125 ! Print headline 02126 02127 IF (PRESENT(spin)) THEN 02128 WRITE (UNIT=iw,FMT="(/,/,T2,A)") spin//" "//TRIM(name) 02129 ELSE 02130 WRITE (UNIT=iw,FMT="(/,/,T2,A)") TRIM(name) 02131 END IF 02132 02133 ! Check if only a subset of the MOs has to be printed 02134 IF ((mo_index_range(1) > 0).AND.& 02135 (mo_index_range(2) > 0).AND.& 02136 (mo_index_range(2) >= mo_index_range(1))) THEN 02137 first_mo = MAX(1,mo_index_range(1)) 02138 last_mo = MIN(mo_set%nmo,mo_index_range(2)) 02139 ELSE 02140 first_mo = 1 02141 last_mo = mo_set%nmo 02142 END IF 02143 02144 IF (p_evec) THEN 02145 02146 ! Print full MO information 02147 02148 DO icol=first_mo,last_mo,ncol 02149 02150 from = icol 02151 to = MIN((from+ncol-1),last_mo) 02152 02153 WRITE (UNIT=iw,FMT=fmtstr1) (jcol,jcol=from,to) 02154 WRITE (UNIT=iw,FMT=fmtstr2) (mo_set%eigenvalues(jcol),jcol=from,to) 02155 WRITE (UNIT=iw,FMT="(A)") "" 02156 02157 WRITE (UNIT=iw,FMT=fmtstr2) (mo_set%occupation_numbers(jcol),jcol=from,to) 02158 WRITE (UNIT=iw,FMT="(A)") "" 02159 02160 irow = 1 02161 02162 DO iatom=1,natom 02163 02164 IF (iatom /= 1) WRITE (UNIT=iw,FMT="(A)") "" 02165 02166 NULLIFY(orb_basis_set,dftb_parameter) 02167 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind,& 02168 element_symbol=element_symbol,& 02169 orb_basis_set=orb_basis_set,& 02170 dftb_parameter=dftb_parameter) 02171 02172 IF (p_cart) THEN 02173 02174 IF (ASSOCIATED(orb_basis_set)) THEN 02175 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 02176 nset=nset,& 02177 nshell=nshell,& 02178 l=l,& 02179 cgf_symbol=bcgf_symbol) 02180 02181 icgf = 1 02182 DO iset=1,nset 02183 DO ishell=1,nshell(iset) 02184 lshell = l(ishell,iset) 02185 DO ico=1,nco(lshell) 02186 WRITE (UNIT=iw,FMT=fmtstr3)& 02187 irow,iatom,ADJUSTR(element_symbol),bcgf_symbol(icgf),& 02188 (cmatrix(irow,jcol),jcol=from,to) 02189 icgf = icgf + 1 02190 irow = irow + 1 02191 END DO 02192 END DO 02193 END DO 02194 ELSE IF (ASSOCIATED(dftb_parameter)) THEN 02195 CALL get_dftb_atom_param(dftb_parameter,lmax=lmax) 02196 icgf = 1 02197 DO ishell=1,lmax+1 02198 lshell = ishell-1 02199 DO ico=1,nco(lshell) 02200 symbol = cgf_symbol(1,indco(1:3,icgf)) 02201 symbol(1:2) = " " 02202 WRITE (UNIT=iw,FMT=fmtstr3)& 02203 irow,iatom,ADJUSTR(element_symbol),symbol,& 02204 (cmatrix(irow,jcol),jcol=from,to) 02205 icgf = icgf + 1 02206 irow = irow + 1 02207 END DO 02208 END DO 02209 ELSE 02210 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 02211 "Unknown basis set type. "//& 02212 CPSourceFileRef,& 02213 only_ionode=.TRUE.) 02214 END IF 02215 02216 ELSE 02217 02218 IF (ASSOCIATED(orb_basis_set)) THEN 02219 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,& 02220 nset=nset,& 02221 nshell=nshell,& 02222 l=l,& 02223 sgf_symbol=bsgf_symbol) 02224 isgf = 1 02225 DO iset=1,nset 02226 DO ishell=1,nshell(iset) 02227 lshell = l(ishell,iset) 02228 DO iso=1,nso(lshell) 02229 WRITE (UNIT=iw,FMT=fmtstr3)& 02230 irow,iatom,ADJUSTR(element_symbol),bsgf_symbol(isgf),& 02231 (smatrix(irow,jcol),jcol=from,to) 02232 isgf = isgf + 1 02233 irow = irow + 1 02234 END DO 02235 END DO 02236 END DO 02237 ELSE IF (ASSOCIATED(dftb_parameter)) THEN 02238 CALL get_dftb_atom_param(dftb_parameter,lmax=lmax) 02239 isgf = 1 02240 DO ishell=1,lmax+1 02241 lshell = ishell-1 02242 DO iso=1,nso(lshell) 02243 symbol = sgf_symbol(1,lshell,-lshell+iso-1) 02244 symbol(1:2) = " " 02245 WRITE (UNIT=iw,FMT=fmtstr3)& 02246 irow,iatom,ADJUSTR(element_symbol),symbol,& 02247 (smatrix(irow,jcol),jcol=from,to) 02248 isgf = isgf + 1 02249 irow = irow + 1 02250 END DO 02251 END DO 02252 ELSE 02253 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,routineP,& 02254 "Unknown basis type. "//& 02255 CPSourceFileRef,& 02256 only_ionode=.TRUE.) 02257 END IF 02258 02259 END IF ! p_cart 02260 02261 END DO ! iatom 02262 02263 END DO ! icol 02264 02265 WRITE (UNIT=iw,FMT="(/)") 02266 02267 ! Release work storage 02268 02269 DEALLOCATE (smatrix,STAT=istat) 02270 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02271 IF (p_cart) THEN 02272 DEALLOCATE (cmatrix,STAT=istat) 02273 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02274 END IF 02275 02276 ELSE IF (p_occ.OR.p_eval) THEN 02277 02278 fmtstr4 = "(T2,I9,2X,F28. ,1X,F24. )" 02279 WRITE (UNIT=fmtstr4(15:16),FMT="(I2)") after 02280 WRITE (UNIT=fmtstr4(25:26),FMT="(I2)") after 02281 WRITE (UNIT=iw,FMT="(/,A)")& 02282 "# MO index MO eigenvalue [a.u.] MO occupation" 02283 DO imo=first_mo,last_mo 02284 WRITE (UNIT=iw,FMT=fmtstr4)& 02285 imo,mo_set%eigenvalues(imo),mo_set%occupation_numbers(imo) 02286 END DO 02287 fmtstr5 = "(A,T42,F24. ,/)" 02288 WRITE (UNIT=fmtstr5(12:13),FMT="(I2)") after 02289 WRITE (UNIT=iw,FMT=fmtstr5)& 02290 "# Sum",accurate_sum(mo_set%occupation_numbers(:)) 02291 02292 END IF ! p_evec 02293 02294 fmtstr6 = "(A,T17,F24. ,/)" 02295 WRITE (UNIT=fmtstr6(12:13),FMT="(I2)") after 02296 WRITE (UNIT=iw,FMT=fmtstr6) " Fermi energy:",mo_set%mu 02297 IF ((mo_set%uniform_occupation).AND.(last_mo > mo_set%homo)) THEN 02298 gap = mo_set%eigenvalues(mo_set%homo+1) -& 02299 mo_set%eigenvalues(mo_set%homo) 02300 fmtstr7 = "(A,T17,F24. ,A,F6.2,A,/)" 02301 WRITE (UNIT=fmtstr7(12:13),FMT="(I2)") after 02302 WRITE (UNIT=iw,FMT=fmtstr7)& 02303 " HOMO-LUMO gap:",gap," = ",gap*evolt," eV" 02304 END IF 02305 02306 END IF ! iw 02307 02308 CALL cp_print_key_finished_output(iw,logger,dft_section,"PRINT%MO",& 02309 ignore_should_output=should_output,& 02310 error=error) 02311 02312 END SUBROUTINE write_mo_set_to_output_unit 02313 02314 END MODULE qs_mo_types
1.7.3