|
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 ! ***************************************************************************** 00012 MODULE xas_methods 00013 00014 USE ai_overlap_new, ONLY: overlap 00015 USE array_types, ONLY: array_i1d_obj,& 00016 array_new,& 00017 array_nullify,& 00018 array_release 00019 USE atomic_kind_types, ONLY: atomic_kind_type,& 00020 get_atomic_kind 00021 USE basis_set_types, ONLY: & 00022 allocate_sto_basis_set, create_gto_from_sto_basis, & 00023 deallocate_sto_basis_set, get_gto_basis_set, gto_basis_set_type, & 00024 init_orb_basis_set, set_sto_basis_set, srules, sto_basis_set_type 00025 USE cell_types, ONLY: cell_type,& 00026 pbc 00027 USE cp_array_r_utils, ONLY: cp_2d_r_p_type 00028 USE cp_control_types, ONLY: dft_control_type 00029 USE cp_dbcsr_interface, ONLY: cp_dbcsr_copy,& 00030 cp_dbcsr_create,& 00031 cp_dbcsr_init,& 00032 cp_dbcsr_set 00033 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 00034 cp_dbcsr_alloc_block_from_nbl,& 00035 cp_dbcsr_allocate_matrix_set,& 00036 cp_dbcsr_sm_fm_multiply 00037 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00038 USE cp_external_control, ONLY: external_control 00039 USE cp_fm_basic_linalg, ONLY: cp_fm_gemm 00040 USE cp_fm_pool_types, ONLY: fm_pool_create_fm 00041 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00042 cp_fm_struct_release,& 00043 cp_fm_struct_type 00044 USE cp_fm_types, ONLY: cp_fm_create,& 00045 cp_fm_get_element,& 00046 cp_fm_get_submatrix,& 00047 cp_fm_p_type,& 00048 cp_fm_set_all,& 00049 cp_fm_set_submatrix,& 00050 cp_fm_to_fm,& 00051 cp_fm_type 00052 USE cp_output_handling, ONLY: cp_p_file,& 00053 cp_print_key_finished_output,& 00054 cp_print_key_should_output,& 00055 cp_print_key_unit_nr 00056 USE cp_para_types, ONLY: cp_para_env_type 00057 USE dbcsr_types, ONLY: dbcsr_distribution_obj,& 00058 dbcsr_type_antisymmetric 00059 USE dbcsr_util, ONLY: convert_offsets_to_sizes 00060 USE input_constants, ONLY: & 00061 do_loc_none, state_loc_list, state_loc_range, xas_1s_type, & 00062 xas_2p_type, xas_2s_type, xas_dip_len2, xas_dip_vel, xas_dscf, & 00063 xas_scf_general, xas_tp_fh, xas_tp_hh, xas_tp_xfh, xas_tp_xhh, & 00064 xes_tp_val 00065 USE input_section_types, ONLY: section_get_lval,& 00066 section_vals_get_subs_vals,& 00067 section_vals_type,& 00068 section_vals_val_get 00069 USE kinds, ONLY: default_string_length,& 00070 dp 00071 USE mathconstants, ONLY: twopi 00072 USE memory_utilities, ONLY: reallocate 00073 USE orbital_pointers, ONLY: ncoset 00074 USE particle_types, ONLY: get_particle_set,& 00075 particle_type 00076 USE periodic_table, ONLY: ptable 00077 USE physcon, ONLY: evolt 00078 USE qs_density_mixing_types, ONLY: direct_mixing_nr,& 00079 mixing_storage_create 00080 USE qs_diis, ONLY: qs_diis_b_clear,& 00081 qs_diis_b_create 00082 USE qs_environment_types, ONLY: get_qs_env,& 00083 qs_environment_type,& 00084 set_qs_env 00085 USE qs_loc_methods, ONLY: qs_loc_driver,& 00086 qs_print_cubes 00087 USE qs_loc_types, ONLY: localized_wfn_control_type,& 00088 qs_loc_env_create,& 00089 qs_loc_env_new_type,& 00090 qs_loc_env_release 00091 USE qs_loc_utils, ONLY: qs_loc_control_init,& 00092 qs_loc_env_init,& 00093 set_loc_centers,& 00094 set_loc_wfn_lists 00095 USE qs_matrix_pools, ONLY: mpools_get,& 00096 qs_matrix_pools_type 00097 USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues 00098 USE qs_mo_types, ONLY: get_mo_set,& 00099 mo_set_p_type,& 00100 set_mo_set,& 00101 write_mo_set 00102 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 00103 USE qs_operators_ao, ONLY: p_xyz_ao,& 00104 rRc_xyz_ao 00105 USE qs_scf, ONLY: init_scf_run,& 00106 scf_env_cleanup 00107 USE qs_scf_types, ONLY: general_diag_method_nr,& 00108 ot_method_nr,& 00109 qs_scf_env_type 00110 USE scf_control_types, ONLY: init_smear,& 00111 read_smear_section,& 00112 scf_control_type 00113 USE termination, ONLY: stop_program 00114 USE timings, ONLY: timeset,& 00115 timestop 00116 USE xas_control, ONLY: xas_control_type 00117 USE xas_env_types, ONLY: get_xas_env,& 00118 set_xas_env,& 00119 xas_env_create,& 00120 xas_env_release,& 00121 xas_environment_type 00122 USE xas_restart, ONLY: xas_read_restart 00123 USE xas_tp_scf, ONLY: xas_do_tp_scf,& 00124 xes_scf_once 00125 #include "cp_common_uses.h" 00126 00127 IMPLICIT NONE 00128 PRIVATE 00129 00130 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_methods' 00131 00132 ! *** Public subroutines *** 00133 00134 PUBLIC :: xas 00135 00136 CONTAINS 00137 00138 ! ***************************************************************************** 00159 SUBROUTINE xas(qs_env, dft_control, error) 00160 00161 TYPE(qs_environment_type), POINTER :: qs_env 00162 TYPE(dft_control_type), POINTER :: dft_control 00163 TYPE(cp_error_type), INTENT(inout) :: error 00164 00165 CHARACTER(LEN=*), PARAMETER :: routineN = 'xas', 00166 routineP = moduleN//':'//routineN 00167 00168 INTEGER :: handle, homo, i, iat, iatom, ispin, istat, method0, 00169 my_homo(2), my_nelectron(2), nao, nelectron, nexc_atoms, nexc_search, 00170 nmo, nspins, output_unit, state_to_be_excited 00171 INTEGER, DIMENSION(:), POINTER :: state_of_atom 00172 LOGICAL :: ch_method_flags, converged, failure, my_uocc(2), should_stop, 00173 skip_scf, transition_potential 00174 REAL(dp) :: maxocc, occ_estate, 00175 xas_nelectron 00176 REAL(dp), DIMENSION(:), POINTER :: eigenvalues, 00177 occupation_numbers 00178 REAL(dp), DIMENSION(:, :), POINTER :: vecbuffer 00179 TYPE(atomic_kind_type), DIMENSION(:), 00180 POINTER :: atomic_kind_set 00181 TYPE(cell_type), POINTER :: cell 00182 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00183 POINTER :: matrix_ks, op_sm, ostrength_sm 00184 TYPE(cp_fm_p_type), DIMENSION(:), 00185 POINTER :: groundstate_coeff 00186 TYPE(cp_fm_type), POINTER :: all_vectors, excvec_coeff, 00187 mo_coeff 00188 TYPE(cp_logger_type), POINTER :: logger 00189 TYPE(cp_para_env_type), POINTER :: para_env 00190 TYPE(mo_set_p_type), DIMENSION(:), 00191 POINTER :: mos 00192 TYPE(particle_type), DIMENSION(:), 00193 POINTER :: particle_set 00194 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 00195 TYPE(qs_scf_env_type), POINTER :: scf_env 00196 TYPE(scf_control_type), POINTER :: scf_control 00197 TYPE(section_vals_type), POINTER :: dft_section, loc_section, 00198 print_loc_section, 00199 scf_section, xas_section 00200 TYPE(xas_control_type), POINTER :: xas_control 00201 TYPE(xas_environment_type), POINTER :: xas_env 00202 00203 CALL timeset(routineN,handle) 00204 00205 failure = .FALSE. 00206 transition_potential = .FALSE. 00207 skip_scf = .FALSE. 00208 converged = .TRUE. 00209 should_stop = .FALSE. 00210 ch_method_flags = .FALSE. 00211 00212 NULLIFY(logger) 00213 logger => cp_error_get_logger(error) 00214 output_unit= cp_logger_get_default_io_unit(logger) 00215 00216 NULLIFY(xas_env, groundstate_coeff, ostrength_sm, op_sm) 00217 NULLIFY(excvec_coeff, state_of_atom, qs_loc_env, cell) 00218 NULLIFY(occupation_numbers,matrix_ks) 00219 NULLIFY(all_vectors,state_of_atom,xas_control) 00220 NULLIFY(vecbuffer,op_sm) 00221 NULLIFY(dft_section, xas_section, loc_section, print_loc_section) 00222 dft_section => section_vals_get_subs_vals(qs_env%input,"DFT",error=error) 00223 scf_section => section_vals_get_subs_vals(dft_section,"SCF",error=error) 00224 xas_section => section_vals_get_subs_vals(dft_section,"XAS",error=error) 00225 loc_section => section_vals_get_subs_vals(xas_section,"LOCALIZE",error=error) 00226 print_loc_section => section_vals_get_subs_vals(loc_section,"PRINT",error=error) 00227 00228 xas_control => dft_control%xas_control 00229 output_unit = cp_print_key_unit_nr(logger,xas_section,"PRINT%PROGRAM_RUN_INFO",& 00230 extension=".Log",error=error) 00231 IF (output_unit>0) THEN 00232 WRITE (UNIT=output_unit,FMT="(/,T3,A,/,T25,A,/,T3,A,/)")& 00233 REPEAT("=",77),& 00234 "START CORE LEVEL SPECTROSCOPY CALCULATION",& 00235 REPEAT("=",77) 00236 END IF 00237 00238 ! Create the xas environment 00239 CALL get_qs_env(qs_env,xas_env=xas_env,error=error) 00240 IF (.NOT.ASSOCIATED(xas_env)) THEN 00241 IF (output_unit>0) THEN 00242 WRITE (UNIT=output_unit,FMT="(/,T10,A,/)")& 00243 "Create and initialize the xas environment" 00244 END IF 00245 CALL xas_env_create(xas_env, error=error) 00246 CALL xas_env_init(xas_env, xas_control, qs_env, xas_section, logger, error=error) 00247 CALL set_qs_env(qs_env,xas_env=xas_env, error=error) 00248 CALL xas_env_release(xas_env,error=error) 00249 CALL get_qs_env(qs_env,xas_env=xas_env,error=error) 00250 END IF 00251 00252 ! Initialize the type of calculation 00253 NULLIFY(atomic_kind_set, scf_control, mos, para_env, particle_set) 00254 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, & 00255 cell = cell, scf_control=scf_control,& 00256 matrix_ks=matrix_ks,mos=mos, para_env=para_env, & 00257 particle_set=particle_set ,error=error) 00258 00259 ! The eigenstate of the KS Hamiltonian are nedeed 00260 NULLIFY(mo_coeff,eigenvalues) 00261 IF(scf_control%use_ot) THEN 00262 IF (output_unit>0) THEN 00263 WRITE (UNIT=output_unit,FMT="(/,T10,A,/)")& 00264 "Get eigenstates and eigenvalues from ground state MOs" 00265 END IF 00266 DO ispin = 1,dft_control%nspins 00267 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff,nao=nao,nmo=nmo,& 00268 eigenvalues=eigenvalues,homo=homo) 00269 CALL calculate_subspace_eigenvalues(mo_coeff,& 00270 matrix_ks(ispin)%matrix,eigenvalues, & 00271 do_rotation=.TRUE.,error=error) 00272 END DO 00273 END IF 00274 00275 ! Set initial occupation numbers, and store the original ones 00276 my_homo = 0 00277 my_nelectron = 0 00278 DO ispin = 1,dft_control%nspins 00279 CALL get_mo_set(mos(ispin)%mo_set,nelectron=my_nelectron(ispin),maxocc=maxocc,& 00280 homo=my_homo(ispin),uniform_occupation=my_uocc(ispin)) 00281 END DO 00282 00283 nspins = dft_control%nspins 00284 transition_potential = (xas_control%xas_method==xas_tp_hh).OR.& 00285 (xas_control%xas_method==xas_tp_fh).OR.& 00286 (xas_control%xas_method==xas_tp_xhh).OR.& 00287 (xas_control%xas_method==xas_tp_xfh).OR.& 00288 (xas_control%xas_method==xas_dscf) 00289 IF(nspins==1 .AND. transition_potential) THEN 00290 CALL stop_program(routineN,moduleN,__LINE__,& 00291 "XAS with TP method requires LSD calculations") 00292 END IF 00293 00294 ! Set of states among which there is the state to be excited 00295 CALL get_mo_set(mos(1)%mo_set,nao=nao,homo=homo) 00296 IF(xas_control%nexc_search < 0) xas_control%nexc_search = homo 00297 nexc_search = xas_control%nexc_search 00298 00299 CALL get_xas_env(xas_env=xas_env,& 00300 state_of_atom=state_of_atom,all_vectors=all_vectors,& 00301 groundstate_coeff=groundstate_coeff,excvec_coeff=excvec_coeff,& 00302 nexc_atoms=nexc_atoms,occ_estate=occ_estate,xas_nelectron=xas_nelectron,error=error) 00303 00304 CALL set_xas_env(xas_env=xas_env,nexc_search=nexc_search,error=error) 00305 00306 ! SCF for only XES using occupied core and empty homo (only one SCF) 00307 ! Probably better not to do the localization in this case, but only single out the 00308 ! core orbital for the specific atom for which the spectrum is computed 00309 IF(xas_control%xas_method==xes_tp_val .AND. & 00310 xas_control%xes_core_occupation==1.0_dp) THEN 00311 IF (output_unit>0) WRITE(UNIT=output_unit,FMT='(/,/,T10,A)') & 00312 "START Core Level Spectroscopy Calculation for the Emission Spectrum" 00313 IF (output_unit>0) WRITE(UNIT=output_unit,FMT='(T10,A,/,A)')& 00314 "The core state is fully occupied and the homo is empty",& 00315 " (final state of the core hole decay). Only one SCF is needed (not one per atom)" 00316 skip_scf = .TRUE. 00317 00318 CALL set_xas_env(xas_env=xas_env,xas_estate=1,error=error) 00319 CALL xes_scf_once(qs_env,xas_env,scf_section,converged,should_stop,error=error) 00320 00321 IF(converged .AND. .NOT. should_stop) THEN 00322 00323 IF (output_unit>0) WRITE(UNIT=output_unit,FMT='(/,T10,A,I6)') & 00324 "SCF with empty homo converged " 00325 ELSE 00326 IF (output_unit>0) WRITE(UNIT=output_unit,FMT='(/,T10,A,I6)') & 00327 "SCF with empty homo NOT converged" 00328 GOTO 1000 00329 END IF 00330 00331 END IF 00332 00333 !Define the qs_loc_env : to find centers, spread and possibly localize them 00334 CALL get_xas_env(xas_env=xas_env, qs_loc_env=qs_loc_env,error=error) 00335 IF (qs_loc_env%do_localize) THEN 00336 IF(output_unit>0) THEN 00337 WRITE (UNIT=output_unit,FMT="(/,T10,A,/)")& 00338 "Localize a sub-set of MOs with spin alpha,"//& 00339 " to better identify the core states that have to be excited" 00340 IF(& 00341 qs_loc_env%localized_wfn_control%set_of_states == state_loc_range ) THEN 00342 WRITE (UNIT=output_unit,FMT="( A , I7, A, I7)") " The sub-set contains states from ",& 00343 qs_loc_env%localized_wfn_control%lu_bound_states(1,1), " to ",& 00344 qs_loc_env%localized_wfn_control%lu_bound_states(2,1) 00345 ELSEIF (qs_loc_env%localized_wfn_control%set_of_states==state_loc_list) THEN 00346 WRITE (UNIT=output_unit,FMT="( A )") " The sub-set contains states given in the input list" 00347 END IF 00348 00349 END IF 00350 CALL qs_loc_driver(qs_env,qs_loc_env,loc_section,print_loc_section,myspin=1,error=error) 00351 END IF 00352 00353 ! Assign the character of the selected core states 00354 ! through the overlap with atomic-like states 00355 CALL cls_assign_core_states(xas_control,xas_env,qs_loc_env%localized_wfn_control,& 00356 qs_env,error=error) 00357 00358 IF(skip_scf) THEN 00359 CALL get_mo_set(mos(1)%mo_set,mo_coeff=mo_coeff) 00360 CALL cp_fm_to_fm(mo_coeff,all_vectors,ncol=nexc_search,& 00361 source_start=1,target_start=1) 00362 END IF 00363 00364 ALLOCATE(vecbuffer(1,nao),STAT=istat) 00365 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00366 ALLOCATE (op_sm(3),STAT=istat) 00367 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00368 00369 ! copy the coefficients of the mos in a temporary fm with the right structure 00370 IF(transition_potential) THEN 00371 CPPrecondition(ASSOCIATED(groundstate_coeff),cp_failure_level,routineP,error,failure) 00372 DO ispin = 1,nspins 00373 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff, nao=nao, nmo=nmo) 00374 CALL cp_fm_to_fm(mo_coeff,groundstate_coeff(ispin)%matrix,nmo,1,1) 00375 END DO 00376 00377 ! Calculate the operator 00378 CALL get_xas_env(xas_env=xas_env,ostrength_sm=ostrength_sm,error=error) 00379 DO i = 1,3 00380 NULLIFY(op_sm(i)%matrix) 00381 op_sm(i)%matrix => ostrength_sm(i)%matrix 00382 END DO 00383 IF(xas_control%dipole_form==xas_dip_vel) THEN 00384 CALL p_xyz_ao(op_sm,qs_env,error=error) 00385 END IF 00386 END IF 00387 ! for td-dft-xas the oscillator strength should be obtained from the linear response orbitals 00388 00389 ! DO SCF if required 00390 NULLIFY(scf_env) 00391 CALL get_qs_env(qs_env,scf_env=scf_env,error=error) 00392 ! OT cannot be used with different occupation numbers in the core 00393 ! Therefore in case of XAS the general_diag_method_nr is imposed 00394 ch_method_flags = .FALSE. 00395 IF( occ_estate .LT. 1.0_dp .AND. scf_env%method == ot_method_nr) THEN 00396 CPPostcondition(xas_control%scf_method==xas_scf_general,cp_warning_level,routineP,error,failure) 00397 method0 = scf_env%method 00398 scf_env%method = general_diag_method_nr 00399 ch_method_flags = .TRUE. 00400 END IF 00401 00402 !some further setup depending on the mixing method 00403 scf_env%mixing_method=xas_env%mixing_method 00404 scf_env%p_mix_alpha = xas_env%mixing_store%alpha 00405 scf_env%skip_diis = .FALSE. 00406 IF(xas_env%mixing_method>direct_mixing_nr) THEN 00407 scf_env%skip_diis = .TRUE. 00408 IF(xas_env%mixing_store%beta==0.0_dp) THEN 00409 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00410 routineP,"CLS: Mixing employing the Kerker damping factor needs BETA /= 0.0",error,failure) 00411 END IF 00412 END IF 00413 IF (xas_env%mixing_method==direct_mixing_nr) THEN 00414 IF(xas_env%eps_diis<qs_env%scf_control%eps_scf) THEN 00415 scf_env%skip_diis = .TRUE. 00416 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 00417 "CLS: the DIIS scheme is disabled, since EPS_DIIS < EPS_SCF",& 00418 error,failure) 00419 END IF 00420 END IF 00421 !dbg 00422 ! write(*,*) xas_env%mixing_method, xas_env%mixing_store%alpha 00423 ! STOP 00424 !dbg 00425 00426 DO iat = 1,nexc_atoms 00427 iatom = xas_env%exc_atoms(iat) 00428 ! determine which state has to be excited in the global list 00429 state_to_be_excited = state_of_atom(iat) 00430 00431 ! Take the state_to_be_excited vector from the full set and copy into excvec_coeff 00432 CALL get_mo_set(mos(1)%mo_set, & 00433 occupation_numbers=occupation_numbers,homo=homo,nmo=nmo,nelectron=nelectron) 00434 IF(xas_control%xas_method==xas_dscf .OR. xas_control%xas_method==xas_tp_xhh & 00435 .OR. xas_control%xas_method==xas_tp_xfh) THEN 00436 CALL cp_assert(nmo>(nelectron+1),cp_failure_level,cp_assertion_failed,& 00437 routineP,"CLS: the required method needs added_mos to the ground state",error,failure) 00438 END IF 00439 occupation_numbers(1:homo) = 1.0_dp 00440 ! If the restart file for this atom exists, the mos and the 00441 ! occupation numbers are overwritten 00442 ! It is necessary that the restart is for the same xas method 00443 ! otherwise the number of electrons and the occupation numbers 00444 ! may not be consistent 00445 IF(xas_control%xas_restart) THEN 00446 CALL xas_read_restart(xas_env,xas_section, qs_env, xas_control%xas_method, iatom,& 00447 state_to_be_excited,error=error) 00448 END IF 00449 00450 CALL set_xas_env(xas_env=xas_env,xas_estate=state_to_be_excited,error=error) 00451 CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff) 00452 CPPrecondition(ASSOCIATED(excvec_coeff),cp_failure_level,routineP,error,failure) 00453 CALL cp_fm_get_submatrix(mo_coeff,vecbuffer,1,state_to_be_excited,& 00454 nao,1,transpose=.TRUE.,error=error) 00455 CALL cp_fm_set_submatrix(excvec_coeff,vecbuffer,1,1,& 00456 nao,1,transpose=.TRUE.,error=error) 00457 00458 IF (transition_potential) THEN 00459 IF(.NOT. skip_scf) THEN 00460 IF (output_unit>0) THEN 00461 IF(xas_control%xas_method==xas_dscf) THEN 00462 WRITE(UNIT=output_unit,FMT='(/,/,T10,A,I6)') & 00463 "START DeltaSCF for the first excited state from the core state of ATOM ", iatom 00464 ELSE 00465 WRITE(UNIT=output_unit,FMT='(/,/,T10,A,I6)') & 00466 "Start Core Level Spectroscopy Calculation with TP approach for ATOM ", iatom 00467 WRITE(UNIT=output_unit,FMT='(T10,A,f10.4)') "Occupation of the core orbital",& 00468 occ_estate 00469 IF(xas_nelectron==nelectron) THEN 00470 WRITE(UNIT=output_unit,FMT='(T10,A,f10.4)') & 00471 "Occupation of lowest virtual orbitals to keep the same number of electrons"//& 00472 " as in the groundstate" 00473 END IF 00474 END IF 00475 END IF 00476 00477 CALL init_scf_run(scf_env=scf_env,qs_env=qs_env,& 00478 scf_section=scf_section, do_xas_tp=.TRUE., error=error) 00479 ! if using mo_coeff_b then copy to fm 00480 DO ispin=1,SIZE(mos)!fm->dbcsr 00481 IF(mos(1)%mo_set%use_mo_coeff_b)THEN!fm->dbcsr 00482 !write(*,*) routinen//' copy_dbcsr_to_fm',__LINE__ 00483 CALL copy_dbcsr_to_fm(mos(ispin)%mo_set%mo_coeff_b,mos(ispin)%mo_set%mo_coeff,error=error)!fm->dbcsr 00484 ENDIF!fm->dbcsr 00485 ENDDO!fm->dbcsr 00486 00487 IF(.NOT.scf_env%skip_diis) THEN 00488 IF (.NOT.ASSOCIATED(scf_env%scf_diis_buffer)) THEN 00489 CALL qs_diis_b_create(scf_env%scf_diis_buffer,nbuffer=scf_control%max_diis,error=error) 00490 END IF 00491 CALL qs_diis_b_clear(scf_env%scf_diis_buffer,error=error) 00492 END IF 00493 00494 CALL xas_do_tp_scf(dft_control,xas_env,iatom,scf_env,qs_env,& 00495 xas_section,scf_section,converged,should_stop,error=error) 00496 00497 CALL external_control(should_stop,"CLS",target_time=qs_env%target_time,& 00498 start_time=qs_env%start_time,error=error) 00499 IF(should_stop)THEN 00500 CALL scf_env_cleanup(scf_env,qs_env=qs_env,error=error) 00501 GO TO 1000 00502 END IF 00503 IF(iat == nexc_atoms) THEN 00504 CALL scf_env_cleanup(scf_env,qs_env=qs_env,error=error) 00505 END IF 00506 00507 END IF 00508 ! SCF DONE 00509 00510 ! *** Write last wavefunction to screen *** 00511 DO ispin=1,dft_control%nspins 00512 CALL write_mo_set(mos(ispin)%mo_set,atomic_kind_set,particle_set,4,& 00513 dft_section,error=error) 00514 ENDDO 00515 00516 ELSE 00517 ! Core level spectroscopy by TDDFT is not yet implemented 00518 ! the states defined by the rotation are the ground state orbitals 00519 ! the initial state from which I excite should be localized 00520 ! I take the excitations from lumo to nmo 00521 END IF 00522 00523 IF(converged) THEN 00524 CALL cls_calculate_spectrum(xas_control,xas_env,qs_env,xas_section,& 00525 iatom,error=error) 00526 ELSE 00527 IF (output_unit>0) WRITE(UNIT=output_unit,FMT='(/,/,T10,A,I6)') & 00528 "SCF with core hole NOT converged for ATOM ", iatom 00529 END IF 00530 00531 IF(.NOT.skip_scf) THEN 00532 ! Reset the initial core orbitals. 00533 ! The valence orbitals are taken from the last SCF, 00534 ! it should be a better initial guess 00535 DO ispin = 1,nspins 00536 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff,& 00537 nao=nao, nmo=nmo) 00538 CALL cp_fm_to_fm(groundstate_coeff(ispin)%matrix,mo_coeff,nmo,1,1) 00539 END DO 00540 END IF 00541 00542 END DO 00543 00544 1000 CONTINUE ! END of Calculation 00545 00546 ! Set back the original method for the optimize 00547 IF(ch_method_flags) THEN 00548 scf_env%method = method0 00549 END IF 00550 00551 ! Release what has to be released 00552 IF(ASSOCIATED(vecbuffer)) THEN 00553 DEALLOCATE(vecbuffer,STAT=istat) 00554 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00555 DEALLOCATE(op_sm,STAT=istat) 00556 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00557 END IF 00558 00559 00560 DO ispin = 1,dft_control%nspins 00561 CALL set_mo_set(mos(ispin)%mo_set, homo=my_homo(ispin),& 00562 uniform_occupation=my_uocc(ispin), error=error) 00563 END DO 00564 00565 xas_env % xas_estate = -1 00566 00567 IF (output_unit>0) THEN 00568 WRITE (UNIT=output_unit,FMT="(/,T3,A,/,T25,A,/,T3,A,/)")& 00569 REPEAT("=",77),& 00570 "END CORE LEVEL SPECTROSCOPY CALCULATION",& 00571 REPEAT("=",77) 00572 END IF 00573 00574 CALL cp_print_key_finished_output(output_unit,logger,xas_section,& 00575 "PRINT%PROGRAM_RUN_INFO",error=error) 00576 CALL timestop(handle) 00577 00578 END SUBROUTINE xas 00579 00580 ! ***************************************************************************** 00591 SUBROUTINE xas_env_init(xas_env, xas_control, qs_env, xas_section, logger, error) 00592 00593 TYPE(xas_environment_type), POINTER :: xas_env 00594 TYPE(xas_control_type) :: xas_control 00595 TYPE(qs_environment_type), POINTER :: qs_env 00596 TYPE(section_vals_type), POINTER :: xas_section 00597 TYPE(cp_logger_type), POINTER :: logger 00598 TYPE(cp_error_type), INTENT(inout) :: error 00599 00600 CHARACTER(len=*), PARAMETER :: routineN = 'xas_env_init', 00601 routineP = moduleN//':'//routineN 00602 00603 CHARACTER(LEN=default_string_length) :: name_sto 00604 INTEGER :: homo, i, iat, iatom, ik, ikind, ispin, istat, j, l, lfomo, 00605 n_mo(2), n_rep, nao, natom, ncubes, nelectron, nexc_atoms, nexc_search, 00606 nj, nk, nkind, nmo, nmoloc(2), norb(0:3), nsgf_gto, nsgf_sto, nspins, 00607 nvirtual, nvirtual2, nwork 00608 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, kind_type_tmp, 00609 kind_z_tmp, last_sgf 00610 INTEGER, DIMENSION(4, 7) :: ne 00611 INTEGER, DIMENSION(:), POINTER :: bounds, list, lq, nq, rbs 00612 LOGICAL :: failure, ihavethis 00613 REAL(dp) :: eps_diis, nele, occ_estate, 00614 zatom 00615 REAL(dp), DIMENSION(:), POINTER :: sto_zet 00616 TYPE(array_i1d_obj) :: row_blk_sizes 00617 TYPE(atomic_kind_type), DIMENSION(:), 00618 POINTER :: atomic_kind_set 00619 TYPE(atomic_kind_type), POINTER :: atomic_kind 00620 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00621 POINTER :: matrix_s 00622 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 00623 TYPE(cp_fm_type), POINTER :: mo_coeff 00624 TYPE(cp_para_env_type), POINTER :: para_env 00625 TYPE(dbcsr_distribution_obj), POINTER :: dbcsr_dist 00626 TYPE(dft_control_type), POINTER :: dft_control 00627 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 00628 TYPE(mo_set_p_type), DIMENSION(:), 00629 POINTER :: mos 00630 TYPE(neighbor_list_set_p_type), 00631 DIMENSION(:), POINTER :: sab_orb 00632 TYPE(particle_type), DIMENSION(:), 00633 POINTER :: particle_set 00634 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 00635 TYPE(qs_matrix_pools_type), POINTER :: mpools 00636 TYPE(qs_scf_env_type), POINTER :: scf_env 00637 TYPE(scf_control_type), POINTER :: scf_control 00638 TYPE(section_vals_type), POINTER :: loc_section, mixing_section, 00639 smear_section 00640 TYPE(sto_basis_set_type), POINTER :: sto_basis_set 00641 00642 failure=.FALSE. 00643 00644 n_mo(1:2) = 0 00645 CPPrecondition(ASSOCIATED(xas_env),cp_failure_level,routineP,error,failure) 00646 00647 IF(.NOT. failure) THEN 00648 00649 NULLIFY( atomic_kind_set, dft_control, scf_control, matrix_s, mos, mpools ) 00650 NULLIFY( para_env, particle_set, scf_env) 00651 NULLIFY(qs_loc_env) 00652 NULLIFY(sab_orb) 00653 CALL get_qs_env(qs_env=qs_env, & 00654 atomic_kind_set=atomic_kind_set, & 00655 dft_control = dft_control, & 00656 scf_control = scf_control, & 00657 mpools=mpools,& 00658 matrix_s=matrix_s, mos=mos, & 00659 para_env=para_env, particle_set=particle_set,& 00660 scf_env=scf_env, & 00661 sab_orb=sab_orb,& 00662 dbcsr_dist=dbcsr_dist,& 00663 error=error) 00664 00665 nexc_search = xas_control%nexc_search 00666 IF(nexc_search < 0) THEN 00667 ! ground state occupation 00668 CALL get_mo_set(mos(1)%mo_set,nmo=nmo,lfomo=lfomo) 00669 nexc_search = lfomo -1 00670 END IF 00671 nexc_atoms = xas_control%nexc_atoms 00672 ALLOCATE(xas_env%exc_atoms(nexc_atoms),STAT=istat) 00673 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00674 xas_env%exc_atoms = xas_control%exc_atoms 00675 eps_diis = xas_control%eps_diis 00676 IF(eps_diis .LT. 0.0_dp) eps_diis = scf_control%eps_diis 00677 CALL set_xas_env(xas_env=xas_env,eps_diis=eps_diis,nexc_search=nexc_search,& 00678 nexc_atoms=nexc_atoms,error=error) 00679 00680 CALL mpools_get(mpools, ao_mo_fm_pools= xas_env%ao_mo_fm_pools,error=error) 00681 00682 NULLIFY(mo_coeff) 00683 CALL get_mo_set(mos(1)%mo_set,nao=nao,homo=homo,nmo=nmo,mo_coeff=mo_coeff, nelectron=nelectron) 00684 00685 nvirtual2 = 0 00686 IF(xas_control%added_mos .GT. 0) THEN 00687 nvirtual2 = MIN(xas_control%added_mos,nao-nmo) 00688 xas_env%unoccupied_eps = xas_control%eps_added 00689 xas_env%unoccupied_max_iter = xas_control%max_iter_added 00690 END IF 00691 nvirtual = nmo + nvirtual2 00692 00693 n_mo(1:2) = nmo 00694 00695 ALLOCATE(xas_env%centers_wfn(3,nexc_search),STAT=istat) 00696 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00697 ALLOCATE(xas_env%atom_of_state(nexc_search),STAT=istat) 00698 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00699 ALLOCATE(xas_env%type_of_state(nexc_search),STAT=istat) 00700 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00701 ALLOCATE(xas_env%state_of_atom(nexc_atoms),STAT=istat) 00702 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00703 ALLOCATE(xas_env%mykind_of_atom(nexc_atoms),STAT=istat) 00704 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00705 nkind = SIZE(atomic_kind_set,1) 00706 ALLOCATE(xas_env%mykind_of_kind(nkind),STAT=istat) 00707 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00708 xas_env%mykind_of_kind = 0 00709 00710 ! create a new matrix structure nao x 1 00711 NULLIFY(tmp_fm_struct) 00712 CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=nao,& 00713 ncol_global=1,para_env=para_env,context=mo_coeff%matrix_struct%context,error=error) 00714 CALL cp_fm_create (xas_env%excvec_coeff, tmp_fm_struct ,error=error) 00715 CALL cp_fm_struct_release ( tmp_fm_struct ,error=error) 00716 00717 NULLIFY(tmp_fm_struct) 00718 CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=1,& 00719 ncol_global=nexc_search,para_env=para_env,& 00720 context=mo_coeff%matrix_struct%context,error=error) 00721 CALL cp_fm_create (xas_env%excvec_overlap, tmp_fm_struct ,error=error) 00722 CALL cp_fm_struct_release ( tmp_fm_struct ,error=error) 00723 00724 nspins = SIZE(mos,1) 00725 00726 ! initialize operators for the calculation of the oscillator strengts 00727 IF (xas_control%xas_method==xas_tp_hh) THEN 00728 occ_estate = 0.5_dp 00729 nele = REAL(nelectron,dp) - 0.5_dp 00730 ELSEIF(xas_control%xas_method==xas_tp_xhh) THEN 00731 occ_estate = 0.5_dp 00732 nele = REAL(nelectron,dp) 00733 ELSEIF(xas_control%xas_method==xas_tp_fh) THEN 00734 occ_estate = 0.0_dp 00735 nele = REAL(nelectron,dp) - 1.0_dp 00736 ELSEIF(xas_control%xas_method==xas_tp_xfh) THEN 00737 occ_estate = 0.0_dp 00738 nele = REAL(nelectron,dp) 00739 ELSEIF(xas_control%xas_method==xes_tp_val) THEN 00740 occ_estate = xas_control%xes_core_occupation 00741 nele = REAL(nelectron,dp)-xas_control%xes_core_occupation 00742 ELSEIF(xas_control%xas_method==xas_dscf) THEN 00743 occ_estate = 0.0_dp 00744 nele = REAL(nelectron,dp) 00745 ENDIF 00746 CALL set_xas_env(xas_env=xas_env,occ_estate=occ_estate,xas_nelectron=nele,& 00747 nvirtual2=nvirtual2,nvirtual=nvirtual,error=error) 00748 00749 ! Initialize the list of orbitals for cube files printing 00750 IF (BTEST(cp_print_key_should_output(logger%iter_info,xas_section,& 00751 "PRINT%CLS_FUNCTION_CUBES",error=error),cp_p_file)) THEN 00752 NULLIFY(bounds,list) 00753 CALL section_vals_val_get(xas_section,& 00754 "PRINT%CLS_FUNCTION_CUBES%CUBES_LU_BOUNDS",& 00755 i_vals=bounds,error=error) 00756 ncubes = bounds(2) - bounds(1) + 1 00757 IF(ncubes > 0 ) THEN 00758 ALLOCATE( xas_control%list_cubes(ncubes),STAT=istat) 00759 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00760 00761 DO ik = 1,ncubes 00762 xas_control%list_cubes(ik) = bounds(1) + (ik-1) 00763 END DO 00764 END IF 00765 00766 IF(.NOT. ASSOCIATED(xas_control%list_cubes)) THEN 00767 CALL section_vals_val_get(xas_section, & 00768 "PRINT%CLS_FUNCTION_CUBES%CUBES_LIST",& 00769 n_rep_val=n_rep,error=error) 00770 ncubes = 0 00771 DO ik = 1,n_rep 00772 NULLIFY(list) 00773 CALL section_vals_val_get(xas_section,& 00774 "PRINT%CLS_FUNCTION_CUBES%CUBES_LIST",& 00775 i_rep_val=ik,i_vals=list,error=error) 00776 IF(ASSOCIATED(list)) THEN 00777 CALL reallocate(xas_control%list_cubes,1,ncubes+ SIZE(list)) 00778 DO i = 1, SIZE(list) 00779 xas_control%list_cubes(i+ncubes) = list(i) 00780 END DO 00781 ncubes = ncubes + SIZE(list) 00782 END IF 00783 END DO ! ik 00784 END IF 00785 00786 IF(.NOT. ASSOCIATED(xas_control%list_cubes)) THEN 00787 ncubes = MAX(10,xas_control%added_mos/10) 00788 ncubes = MIN(ncubes,xas_control%added_mos) 00789 ALLOCATE( xas_control%list_cubes(ncubes),STAT=istat) 00790 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00791 DO ik = 1,ncubes 00792 xas_control%list_cubes(ik) = homo + ik 00793 END DO 00794 END IF 00795 ELSE 00796 NULLIFY(xas_control%list_cubes) 00797 ENDIF 00798 00799 NULLIFY(tmp_fm_struct) 00800 nwork = nvirtual 00801 CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=nao,ncol_global=nwork,& 00802 para_env=para_env,context=mo_coeff%matrix_struct%context,error=error) 00803 CALL cp_fm_create (xas_env%fm_work,tmp_fm_struct,error=error) 00804 CALL cp_fm_struct_release ( tmp_fm_struct ,error=error) 00805 00806 ALLOCATE (xas_env%groundstate_coeff(nspins), STAT=istat) 00807 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00808 DO ispin = 1,nspins 00809 NULLIFY(xas_env%groundstate_coeff(ispin)%matrix) 00810 CALL get_mo_set(mos(ispin)%mo_set,nao=nao,nmo=nmo) 00811 CALL fm_pool_create_fm(xas_env%ao_mo_fm_pools(ispin)%pool,& 00812 xas_env%groundstate_coeff(ispin)%matrix,& 00813 name="xas_env%mo0"//TRIM(ADJUSTL(cp_to_string(ispin))),error=error) 00814 END DO ! ispin 00815 00816 NULLIFY(tmp_fm_struct) 00817 CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=1,& 00818 ncol_global=nvirtual,para_env=para_env,& 00819 context=mo_coeff%matrix_struct%context,error=error) 00820 ALLOCATE (xas_env%dip_fm_set(2,3),STAT=istat) 00821 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00822 DO i = 1,3 00823 DO j = 1,2 00824 NULLIFY(xas_env%dip_fm_set(j,i)%matrix) 00825 CALL cp_fm_create (xas_env%dip_fm_set(j,i)%matrix, tmp_fm_struct ,error=error) 00826 END DO 00827 END DO 00828 CALL cp_fm_struct_release ( tmp_fm_struct ,error=error) 00829 00830 !Array to store all the eigenstates: occupied and the required not occupied 00831 IF(nvirtual2 .GT. 0) THEN 00832 ALLOCATE(xas_env%unoccupied_evals(nvirtual2), STAT=istat) 00833 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00834 NULLIFY(tmp_fm_struct) 00835 CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=nao,& 00836 ncol_global=nvirtual2,& 00837 para_env=para_env,context=mo_coeff%matrix_struct%context,error=error) 00838 NULLIFY(xas_env%unoccupied_orbs) 00839 CALL cp_fm_create (xas_env%unoccupied_orbs,tmp_fm_struct,error=error) 00840 CALL cp_fm_struct_release ( tmp_fm_struct ,error=error) 00841 END IF 00842 00843 NULLIFY(tmp_fm_struct) 00844 CALL cp_fm_struct_create(tmp_fm_struct,nrow_global=nao, & 00845 ncol_global=nvirtual,& 00846 para_env=para_env,context=mo_coeff%matrix_struct%context,error=error) 00847 NULLIFY(xas_env%all_vectors) 00848 CALL cp_fm_create (xas_env%all_vectors,tmp_fm_struct,error=error) 00849 CALL cp_fm_struct_release ( tmp_fm_struct ,error=error) 00850 00851 ! Array to store all the energies needed for the spectrum 00852 ALLOCATE(xas_env%all_evals(nvirtual), STAT=istat) 00853 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00854 00855 IF(xas_control%dipole_form==xas_dip_len2) THEN 00856 CALL cp_dbcsr_allocate_matrix_set(xas_env%ostrength_sm,3,error=error) 00857 DO i = 1,3 00858 ALLOCATE(xas_env%ostrength_sm(i)%matrix) 00859 CALL cp_dbcsr_init(xas_env%ostrength_sm(i)%matrix,error=error) 00860 CALL cp_dbcsr_copy(xas_env%ostrength_sm(i)%matrix,matrix_s(1)%matrix,& 00861 "xas_env%ostrength_sm"//& 00862 "-"//TRIM(ADJUSTL(cp_to_string(i))),error=error) 00863 CALL cp_dbcsr_set(xas_env%ostrength_sm(i)%matrix,0.0_dp,error=error) 00864 END DO 00865 ELSEIF(xas_control%dipole_form==xas_dip_vel) THEN 00866 ! 00867 ! prepare for allocation 00868 natom = SIZE(particle_set,1) 00869 ALLOCATE (first_sgf(natom),STAT=istat) 00870 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00871 ALLOCATE (last_sgf(natom),STAT=istat) 00872 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00873 CALL get_particle_set(particle_set=particle_set,& 00874 first_sgf=first_sgf,& 00875 last_sgf=last_sgf,error=error) 00876 ALLOCATE (rbs(natom), STAT=istat) 00877 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00878 CALL convert_offsets_to_sizes (first_sgf, rbs, last_sgf) 00879 CALL array_nullify (row_blk_sizes) 00880 CALL array_new (row_blk_sizes, rbs, gift=.TRUE.) 00881 DEALLOCATE (first_sgf,STAT=istat) 00882 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00883 DEALLOCATE (last_sgf,STAT=istat) 00884 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00885 ! 00886 ! 00887 CALL cp_dbcsr_allocate_matrix_set(xas_env%ostrength_sm,3,error=error) 00888 ALLOCATE(xas_env%ostrength_sm(1)%matrix) 00889 CALL cp_dbcsr_init(xas_env%ostrength_sm(1)%matrix,error=error) 00890 CALL cp_dbcsr_create(matrix=xas_env%ostrength_sm(1)%matrix, & 00891 name="xas_env%ostrength_sm", & 00892 dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric,& 00893 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, & 00894 nblks=0, nze=0, mutable_work=.TRUE., & 00895 error=error) 00896 CALL cp_dbcsr_alloc_block_from_nbl(xas_env%ostrength_sm(1)%matrix,sab_orb,error=error) 00897 CALL cp_dbcsr_set(xas_env%ostrength_sm(1)%matrix,0.0_dp,error=error) 00898 DO i = 2,3 00899 ALLOCATE(xas_env%ostrength_sm(i)%matrix) 00900 CALL cp_dbcsr_init(xas_env%ostrength_sm(i)%matrix,error=error) 00901 CALL cp_dbcsr_copy(xas_env%ostrength_sm(i)%matrix,xas_env%ostrength_sm(1)%matrix,& 00902 "xas_env%ostrength_sm-"//TRIM(ADJUSTL(cp_to_string(i))),error=error) 00903 CALL cp_dbcsr_set(xas_env%ostrength_sm(i)%matrix,0.0_dp,error=error) 00904 END DO 00905 00906 CALL array_release (row_blk_sizes) 00907 END IF 00908 00909 !Define the qs_loc_env : to find centers, spread and possibly localize them 00910 IF(.NOT.(ASSOCIATED(xas_env%qs_loc_env))) THEN 00911 CALL qs_loc_env_create(qs_loc_env,error=error) 00912 CALL set_xas_env(xas_env=xas_env, qs_loc_env=qs_loc_env,error=error) 00913 CALL qs_loc_env_release(qs_loc_env,error=error) 00914 CALL get_xas_env(xas_env=xas_env, qs_loc_env=qs_loc_env,error=error) 00915 loc_section => section_vals_get_subs_vals(xas_section,"LOCALIZE",error=error) 00916 CALL qs_loc_control_init(qs_loc_env,loc_section,do_homo=.TRUE.,& 00917 do_xas=.TRUE.,nloc_xas=nexc_search,error=error) 00918 IF(.NOT. qs_loc_env%do_localize) THEN 00919 qs_loc_env%localized_wfn_control%localization_method = do_loc_none 00920 ELSE 00921 nmoloc = qs_loc_env%localized_wfn_control%nloc_states 00922 CALL set_loc_wfn_lists(qs_loc_env%localized_wfn_control,nmoloc,n_mo,nspins,error=error) 00923 CALL set_loc_centers(qs_loc_env%localized_wfn_control,nmoloc,nspins,error=error) 00924 CALL qs_loc_env_init(qs_loc_env,qs_loc_env%localized_wfn_control,& 00925 qs_env,myspin=1,do_localize=qs_loc_env%do_localize,error=error) 00926 END IF 00927 END IF 00928 00929 ! Initialize Mixing scheme 00930 mixing_section => section_vals_get_subs_vals(xas_section,"MIXING",error=error) 00931 CALL section_vals_val_get(mixing_section,"METHOD",i_val=xas_env%mixing_method,error=error) 00932 CALL mixing_storage_create(xas_env%mixing_store, mixing_section, xas_env%mixing_method, & 00933 dft_control%qs_control%cutoff, error=error) 00934 00935 CALL init_smear(xas_env%smear,error) 00936 smear_section => section_vals_get_subs_vals(xas_section,"SMEAR",error=error) 00937 CALL read_smear_section(xas_env%smear,smear_section,error=error) 00938 00939 !Type of state 00940 norb = 0 00941 ALLOCATE(nq(1),lq(1),sto_zet(1),STAT=istat) 00942 IF( xas_control%state_type == xas_1s_type) THEN 00943 norb(0) = 1 00944 nq(1) = 1 00945 lq(1) = 0 00946 ELSEIF( xas_control%state_type == xas_2s_type ) THEN 00947 norb(0) = 2 00948 nq(1) = 2 00949 lq(1) = 0 00950 ELSEIF( xas_control%state_type == xas_2p_type ) THEN 00951 norb(0) = 2 00952 norb(1) = 1 00953 nq(1) = 2 00954 lq(1) = 1 00955 ELSE 00956 CALL stop_program(routineN,moduleN,__LINE__,& 00957 "XAS type of state not implemented") 00958 END IF 00959 00960 ALLOCATE(kind_type_tmp(nkind),STAT=istat) 00961 ALLOCATE(kind_z_tmp(nkind),STAT=istat) 00962 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00963 kind_type_tmp = 0 00964 kind_z_tmp = 0 00965 nk=0 00966 DO iat = 1,nexc_atoms 00967 iatom = xas_env%exc_atoms(iat) 00968 NULLIFY(atomic_kind) 00969 atomic_kind => particle_set(iatom)%atomic_kind 00970 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00971 kind_number=ikind,zeff=zatom) 00972 ihavethis = .FALSE. 00973 DO ik = 1,nk 00974 IF(ikind==kind_type_tmp(ik)) THEN 00975 ihavethis = .TRUE. 00976 xas_env%mykind_of_atom(iat) = ik 00977 EXIT 00978 END IF 00979 END DO 00980 IF(.NOT. ihavethis) THEN 00981 nk = nk +1 00982 kind_type_tmp(nk) = ikind 00983 kind_z_tmp(nk) = INT(zatom) 00984 xas_env%mykind_of_atom(iat) = nk 00985 xas_env%mykind_of_kind(ikind) = nk 00986 END IF 00987 END DO ! iat 00988 00989 ALLOCATE(xas_env%my_gto_basis(nk),STAT=istat) 00990 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00991 ALLOCATE(xas_env%stogto_overlap(nk),STAT=istat) 00992 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00993 DO ik = 1,nk 00994 NULLIFY(xas_env%my_gto_basis(ik)%gto_basis_set,sto_basis_set) 00995 ne = 0 00996 DO l=1,4 !lq(1)+1 00997 nj = 2*(l-1)+1 00998 DO i=l, 7! nq(1) 00999 ne(l,i) = ptable(kind_z_tmp(ik))%e_conv(l-1) - 2*nj*(i-l) 01000 ne(l,i) = MAX(ne(l,i),0) 01001 ne(l,i) = MIN(ne(l,i),2*nj) 01002 END DO 01003 END DO 01004 01005 sto_zet(1) = srules(kind_z_tmp(ik),ne,nq(1),lq(1)) 01006 CALL allocate_sto_basis_set(sto_basis_set,error) 01007 name_sto='xas_tmp_sto' 01008 CALL set_sto_basis_set(sto_basis_set,nshell=1,nq=nq,& 01009 lq=lq,zet=sto_zet,name=name_sto) 01010 CALL create_gto_from_sto_basis(sto_basis_set,& 01011 xas_env%my_gto_basis(ik)%gto_basis_set,xas_control%ngauss,error=error) 01012 CALL deallocate_sto_basis_set(sto_basis_set,error) 01013 xas_env%my_gto_basis(ik)%gto_basis_set%norm_type = 2 01014 CALL init_orb_basis_set(xas_env%my_gto_basis(ik)%gto_basis_set,error=error) 01015 01016 atomic_kind => atomic_kind_set(kind_type_tmp(ik)) 01017 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01018 orb_basis_set=orb_basis_set) 01019 01020 CALL get_gto_basis_set(gto_basis_set=orb_basis_set,nsgf=nsgf_gto) 01021 CALL get_gto_basis_set(gto_basis_set=xas_env%my_gto_basis(ik)%gto_basis_set,nsgf=nsgf_sto) 01022 ALLOCATE(xas_env%stogto_overlap(ik)%array(nsgf_sto,nsgf_gto),STAT=istat) 01023 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01024 01025 CALL calc_stogto_overlap(xas_env%my_gto_basis(ik)%gto_basis_set,orb_basis_set,& 01026 xas_env%stogto_overlap(ik)%array,error=error) 01027 END DO 01028 01029 DEALLOCATE(nq,lq,sto_zet,STAT=istat) 01030 DEALLOCATE(kind_type_tmp,kind_z_tmp,STAT=istat) 01031 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01032 END IF ! failure 01033 01034 END SUBROUTINE xas_env_init 01035 01036 ! ***************************************************************************** 01050 SUBROUTINE cls_calculate_spectrum(xas_control,xas_env,qs_env,xas_section,& 01051 iatom,error) 01052 01053 TYPE(xas_control_type) :: xas_control 01054 TYPE(xas_environment_type), POINTER :: xas_env 01055 TYPE(qs_environment_type), POINTER :: qs_env 01056 TYPE(section_vals_type), POINTER :: xas_section 01057 INTEGER, INTENT(IN) :: iatom 01058 TYPE(cp_error_type), INTENT(inout) :: error 01059 01060 CHARACTER(LEN=*), PARAMETER :: routineN = 'cls_calculate_spectrum', 01061 routineP = moduleN//':'//routineN 01062 01063 INTEGER :: homo, i, istat, lfomo, nabs, 01064 nmo, nvirtual, output_unit, 01065 xas_estate 01066 LOGICAL :: append_cube, failure, length 01067 REAL(dp) :: rc(3) 01068 REAL(dp), DIMENSION(:), POINTER :: all_evals 01069 REAL(dp), DIMENSION(:, :), POINTER :: sp_ab, sp_em 01070 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01071 POINTER :: op_sm, ostrength_sm 01072 TYPE(cp_fm_p_type), DIMENSION(:, :), 01073 POINTER :: dip_fm_set 01074 TYPE(cp_fm_type), POINTER :: all_vectors, excvec_coeff, 01075 fm_work 01076 TYPE(cp_logger_type), POINTER :: logger 01077 TYPE(mo_set_p_type), DIMENSION(:), 01078 POINTER :: mos 01079 TYPE(particle_type), DIMENSION(:), 01080 POINTER :: particle_set 01081 01082 failure =.FALSE. 01083 NULLIFY(logger) 01084 logger => cp_error_get_logger(error) 01085 output_unit= cp_logger_get_default_io_unit(logger) 01086 01087 NULLIFY(fm_work,ostrength_sm,op_sm, dip_fm_set) 01088 NULLIFY(all_evals,all_vectors,excvec_coeff) 01089 NULLIFY(mos,particle_set,sp_em,sp_ab) 01090 ALLOCATE (op_sm(3),STAT=istat) 01091 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01092 01093 CALL get_qs_env(qs_env=qs_env, & 01094 mos=mos,particle_set=particle_set,error=error) 01095 01096 CALL get_mo_set(mos(1)%mo_set, homo=homo,lfomo=lfomo,nmo=nmo) 01097 CALL get_xas_env(xas_env=xas_env,all_vectors=all_vectors,xas_estate=xas_estate,& 01098 all_evals=all_evals,dip_fm_set=dip_fm_set,excvec_coeff=excvec_coeff,& 01099 fm_work=fm_work, ostrength_sm=ostrength_sm,nvirtual=nvirtual,error=error) 01100 01101 nabs = nvirtual -lfomo + 1 01102 ALLOCATE (sp_em(6,homo),STAT=istat) 01103 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01104 ALLOCATE (sp_ab(6,nabs),STAT=istat) 01105 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01106 CPPrecondition(ASSOCIATED(excvec_coeff),cp_failure_level,routineP,error,failure) 01107 01108 IF(.NOT. failure) THEN 01109 01110 IF(.NOT. xas_control%xas_method == xas_dscf) THEN 01111 ! Calculate the spectrum 01112 IF(xas_control%dipole_form==xas_dip_len2) THEN 01113 rc(1:3) = particle_set(iatom)%r(1:3) 01114 DO i = 1,3 01115 NULLIFY(op_sm(i)%matrix) 01116 op_sm(i)%matrix => ostrength_sm(i)%matrix 01117 END DO 01118 CALL rRc_xyz_ao(op_sm,qs_env,rc,order=1,minimum_image=.TRUE.,error=error) 01119 CALL spectrum_dip_vel(dip_fm_set,op_sm,mos,excvec_coeff,& 01120 all_vectors,all_evals,fm_work,& 01121 sp_em,sp_ab,xas_estate,nvirtual,error=error) 01122 DO i = 1,SIZE(ostrength_sm,1) 01123 CALL cp_dbcsr_set(ostrength_sm(i)%matrix,0.0_dp,error=error) 01124 END DO 01125 ELSE 01126 DO i = 1,3 01127 NULLIFY(op_sm(i)%matrix) 01128 op_sm(i)%matrix => ostrength_sm(i)%matrix 01129 END DO 01130 CALL spectrum_dip_vel(dip_fm_set,op_sm,mos,excvec_coeff,& 01131 all_vectors,all_evals,fm_work,& 01132 sp_em,sp_ab,xas_estate,nvirtual,error=error) 01133 END IF 01134 END IF 01135 01136 CALL get_mo_set(mos(1)%mo_set, lfomo=lfomo) 01137 ! writw thw spectrum, if the file exists it is appended 01138 IF( .NOT. xas_control%xas_method == xas_dscf) THEN 01139 length = (.NOT. xas_control%dipole_form==xas_dip_vel) 01140 CALL xas_write(sp_em,sp_ab, xas_estate, & 01141 xas_section, iatom, lfomo, length=length, error=error) 01142 END IF 01143 01144 DEALLOCATE(sp_em,STAT=istat) 01145 DEALLOCATE(sp_ab,STAT=istat) 01146 01147 IF(BTEST(cp_print_key_should_output(logger%iter_info,xas_section,& 01148 "PRINT%CLS_FUNCTION_CUBES",error=error),cp_p_file)) THEN 01149 append_cube= section_get_lval(xas_section,"PRINT%CLS_FUNCTION_CUBES%APPEND",error=error) 01150 CALL xas_print_cubes(xas_control,qs_env,xas_section,mos,all_vectors,& 01151 iatom,append_cube,error=error) 01152 END IF 01153 01154 END IF ! failure 01155 01156 DEALLOCATE (op_sm,STAT=istat) 01157 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01158 01159 END SUBROUTINE cls_calculate_spectrum 01160 01161 ! ***************************************************************************** 01175 SUBROUTINE xas_write(sp_em, sp_ab, estate, xas_section, iatom, & 01176 lfomo, length, error) 01177 01178 REAL(dp), DIMENSION(:, :), POINTER :: sp_em, sp_ab 01179 INTEGER, INTENT(IN) :: estate 01180 TYPE(section_vals_type), POINTER :: xas_section 01181 INTEGER, INTENT(IN) :: iatom, lfomo 01182 LOGICAL, INTENT(IN) :: length 01183 TYPE(cp_error_type), INTENT(inout) :: error 01184 01185 CHARACTER(len=*), PARAMETER :: routineN = 'xas_write', 01186 routineP = moduleN//':'//routineN 01187 01188 CHARACTER(LEN=default_string_length) :: mittle_ab, mittle_em, my_act, 01189 my_pos 01190 INTEGER :: i, istate, out_sp_ab, 01191 out_sp_em 01192 LOGICAL :: failure 01193 REAL(dp) :: ene2 01194 TYPE(cp_logger_type), POINTER :: logger 01195 01196 NULLIFY(logger) 01197 logger => cp_error_get_logger(error) 01198 01199 failure = .FALSE. 01200 01201 my_pos = "APPEND" 01202 my_act = "WRITE" 01203 01204 mittle_em = "xes_at"//TRIM(ADJUSTL(cp_to_string(iatom))) 01205 out_sp_em = cp_print_key_unit_nr(logger,xas_section,"PRINT%XES_SPECTRUM",& 01206 extension=".spectrum", file_position=my_pos, file_action=my_act,& 01207 file_form="FORMATTED", middle_name=TRIM(mittle_em), & 01208 error=error) 01209 01210 IF(out_sp_em>0) THEN 01211 WRITE(out_sp_em,'(A,I6,A,I6,A,I6)') " Emission spectrum for atom ", iatom,& 01212 ", index of excited core MO is", estate, ", # of lines ", SIZE(sp_em,2) 01213 ene2=1.0_dp 01214 DO istate = estate,SIZE(sp_em,2) 01215 IF(length) ene2 = sp_em(1,istate)*sp_em(1,istate) 01216 WRITE(out_sp_em,'(I6,5F16.8,F10.5)') istate, sp_em(1,istate)*evolt, & 01217 sp_em(2,istate)*ene2, sp_em(3,istate)*ene2,& 01218 sp_em(4,istate)*ene2, sp_em(5,istate)*ene2,sp_em(6,istate) 01219 END DO 01220 END IF 01221 CALL cp_print_key_finished_output(out_sp_em,logger,xas_section,& 01222 "PRINT%XES_SPECTRUM", error=error) 01223 01224 mittle_ab = "xas_at"//TRIM(ADJUSTL(cp_to_string(iatom))) 01225 out_sp_ab = cp_print_key_unit_nr(logger,xas_section,"PRINT%XAS_SPECTRUM",& 01226 extension=".spectrum", file_position=my_pos, file_action=my_act,& 01227 file_form="FORMATTED", middle_name=TRIM(mittle_ab), & 01228 error=error) 01229 01230 IF(out_sp_ab>0) THEN 01231 WRITE(out_sp_ab,'(A,I6,A,I6,A,I6)') " Absorption spectrum for atom ", iatom,& 01232 ", index of excited core MO is", estate, ", # of lines ", SIZE(sp_ab,2) 01233 ene2=1.0_dp 01234 DO i = 1,SIZE(sp_ab,2) 01235 istate = lfomo -1 + i 01236 IF(length) ene2 = sp_ab(1,i)*sp_ab(1,i) 01237 WRITE(out_sp_ab,'(I6,5F16.8,F10.5)') istate, sp_ab(1,i)*evolt, & 01238 sp_ab(2,i)*ene2, sp_ab(3,i)*ene2,& 01239 sp_ab(4,i)*ene2, sp_ab(5,i)*ene2,sp_ab(6,i) 01240 END DO 01241 END IF 01242 01243 CALL cp_print_key_finished_output(out_sp_ab,logger,xas_section,& 01244 "PRINT%XAS_SPECTRUM", error=error) 01245 01246 END SUBROUTINE xas_write 01247 01248 ! ***************************************************************************** 01260 SUBROUTINE xas_print_cubes(xas_control,qs_env,xas_section,& 01261 mos,all_vectors,iatom,append_cube, error) 01262 01263 TYPE(xas_control_type) :: xas_control 01264 TYPE(qs_environment_type), POINTER :: qs_env 01265 TYPE(section_vals_type), POINTER :: xas_section 01266 TYPE(mo_set_p_type), DIMENSION(:), 01267 POINTER :: mos 01268 TYPE(cp_fm_type), POINTER :: all_vectors 01269 INTEGER, INTENT(IN) :: iatom 01270 LOGICAL, INTENT(IN) :: append_cube 01271 TYPE(cp_error_type), INTENT(inout) :: error 01272 01273 CHARACTER(len=*), PARAMETER :: routineN = 'xas_print_cubes', 01274 routineP = moduleN//':'//routineN 01275 01276 CHARACTER(LEN=default_string_length) :: my_mittle, my_pos 01277 INTEGER :: homo, istat, istate0, nspins, 01278 nstates 01279 LOGICAL :: failure 01280 REAL(dp), DIMENSION(:, :), POINTER :: centers 01281 TYPE(section_vals_type), POINTER :: print_key 01282 01283 failure = .FALSE. 01284 nspins = SIZE(mos) 01285 01286 print_key => section_vals_get_subs_vals(xas_section,"PRINT%CLS_FUNCTION_CUBES",error=error) 01287 my_mittle = 'at'//TRIM(ADJUSTL(cp_to_string(iatom))) 01288 nstates = SIZE(xas_control%list_cubes,1) 01289 01290 IF(xas_control%do_centers) THEN 01291 ! one might like to calculate the centers of the xas orbital (without localizing them) 01292 ELSE 01293 ALLOCATE(centers(6,nstates),STAT=istat) 01294 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01295 centers = 0.0_dp 01296 END IF 01297 01298 CALL get_mo_set(mos(1)%mo_set, homo=homo) 01299 istate0 = 0 01300 01301 my_pos = "REWIND" 01302 IF (append_cube) THEN 01303 my_pos = "APPEND" 01304 END IF 01305 01306 CALL qs_print_cubes(qs_env,all_vectors,nstates,xas_control%list_cubes,& 01307 centers,print_key,my_mittle,state0=istate0,file_position=my_pos,error=error) 01308 01309 DEALLOCATE(centers,STAT=istat) 01310 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01311 01312 END SUBROUTINE xas_print_cubes 01313 01314 ! ***************************************************************************** 01337 SUBROUTINE spectrum_dip_len(fm_set,op_sm,mos,excvec,& 01338 all_vectors,all_evals,fm_work,cell, sp_em, sp_ab,estate, nstate,error) 01339 01340 TYPE(cp_fm_p_type), DIMENSION(:, :), 01341 POINTER :: fm_set 01342 TYPE(cp_dbcsr_p_type), DIMENSION(:, :), 01343 POINTER :: op_sm 01344 TYPE(mo_set_p_type), DIMENSION(:), 01345 POINTER :: mos 01346 TYPE(cp_fm_type), POINTER :: excvec, all_vectors 01347 REAL(dp), DIMENSION(:), POINTER :: all_evals 01348 TYPE(cp_fm_type), POINTER :: fm_work 01349 TYPE(cell_type), POINTER :: cell 01350 REAL(dp), DIMENSION(:, :), POINTER :: sp_em, sp_ab 01351 INTEGER, INTENT(IN) :: estate, nstate 01352 TYPE(cp_error_type), INTENT(inout) :: error 01353 01354 CHARACTER(LEN=*), PARAMETER :: routineN = 'spectrum_dip_len', 01355 routineP = moduleN//':'//routineN 01356 01357 COMPLEX(KIND=dp) :: z 01358 INTEGER :: homo, i, i_abs, istate, j, 01359 lfomo, nao, nmo 01360 LOGICAL :: failure 01361 REAL(dp) :: ene_f, ene_i, imagpart, 01362 ra(3), realpart 01363 REAL(dp), DIMENSION(:), POINTER :: eigenvalues, 01364 occupation_numbers 01365 01366 failure = .FALSE. 01367 01368 CPPrecondition(ASSOCIATED(fm_set),cp_failure_level,routineP,error,failure) 01369 CPPrecondition(ASSOCIATED(mos),cp_failure_level,routineP,error,failure) 01370 01371 NULLIFY(eigenvalues, occupation_numbers) 01372 IF(.NOT. failure) THEN 01373 CALL get_mo_set(mos(1)%mo_set, eigenvalues=eigenvalues, occupation_numbers=occupation_numbers,& 01374 nao=nao, nmo=nmo, homo =homo, lfomo=lfomo) 01375 DO i=1,SIZE(fm_set,2) 01376 DO j = 1,SIZE(fm_set,1) 01377 CPPrecondition(ASSOCIATED(fm_set(j,i)%matrix),cp_failure_level,routineP,error,failure) 01378 CALL cp_fm_set_all(fm_set(j,i)%matrix, 0.0_dp, error=error) 01379 CALL cp_fm_set_all(fm_work, 0.0_dp, error=error) 01380 CALL cp_dbcsr_sm_fm_multiply(op_sm(j,i)%matrix,all_vectors,fm_work,ncol=nstate,error=error) 01381 CALL cp_fm_gemm("T","N",1,nstate,nao,1.0_dp,excvec,& 01382 fm_work,0.0_dp, fm_set(j,i)%matrix,b_first_col=1,error=error) 01383 END DO 01384 END DO 01385 01386 sp_em = 0.0_dp 01387 sp_ab = 0.0_dp 01388 ene_i = eigenvalues(estate) 01389 DO istate = 1,nstate 01390 ene_f = all_evals(istate) 01391 DO i = 1,3 01392 CALL cp_fm_get_element(fm_set(1,i)%matrix,1,istate,realpart) 01393 CALL cp_fm_get_element(fm_set(2,i)%matrix,1,istate,imagpart) 01394 z = CMPLX(realpart,imagpart,dp) 01395 ra(i) = ( cell % hmat ( i, i ) / twopi ) * AIMAG ( LOG ( z ) ) 01396 END DO 01397 IF(istate<=homo) THEN 01398 sp_em(1,istate) = ene_f - ene_i 01399 sp_em(2,istate) = ra(1)*ra(1) 01400 sp_em(3,istate) = ra(2)*ra(2) 01401 sp_em(4,istate) = ra(3)*ra(3) 01402 sp_em(5,istate) = ra(1)*ra(1)+ra(2)*ra(2)+ra(3)*ra(3) 01403 sp_em(6,istate) = occupation_numbers(istate) 01404 END IF 01405 IF(istate>=lfomo ) THEN 01406 i_abs = istate - lfomo + 1 01407 sp_ab(1,i_abs) = ene_f - ene_i 01408 sp_ab(2,i_abs) = ra(1)*ra(1) 01409 sp_ab(3,i_abs) = ra(2)*ra(2) 01410 sp_ab(4,i_abs) = ra(3)*ra(3) 01411 sp_ab(5,i_abs) = ra(1)*ra(1)+ra(2)*ra(2)+ra(3)*ra(3) 01412 IF(istate<=nmo) sp_ab(6,i_abs) = occupation_numbers(istate) 01413 END IF 01414 END DO 01415 END IF 01416 01417 END SUBROUTINE spectrum_dip_len 01418 01419 ! ***************************************************************************** 01440 SUBROUTINE spectrum_dip_vel(fm_set,op_sm,mos,excvec, & 01441 all_vectors,all_evals,fm_work,sp_em, sp_ab,estate, nstate,error) 01442 01443 TYPE(cp_fm_p_type), DIMENSION(:, :), 01444 POINTER :: fm_set 01445 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01446 POINTER :: op_sm 01447 TYPE(mo_set_p_type), DIMENSION(:), 01448 POINTER :: mos 01449 TYPE(cp_fm_type), POINTER :: excvec, all_vectors 01450 REAL(dp), DIMENSION(:), POINTER :: all_evals 01451 TYPE(cp_fm_type), POINTER :: fm_work 01452 REAL(dp), DIMENSION(:, :), POINTER :: sp_em, sp_ab 01453 INTEGER, INTENT(IN) :: estate, nstate 01454 TYPE(cp_error_type), INTENT(inout) :: error 01455 01456 CHARACTER(LEN=*), PARAMETER :: routineN = 'spectrum_dip_vel', 01457 routineP = moduleN//':'//routineN 01458 01459 INTEGER :: homo, i, i_abs, istate, 01460 lfomo, nao, nmo 01461 LOGICAL :: failure 01462 REAL(dp) :: dip(3), ene_f, ene_i 01463 REAL(dp), DIMENSION(:), POINTER :: eigenvalues, 01464 occupation_numbers 01465 01466 failure = .FALSE. 01467 01468 CPPrecondition(ASSOCIATED(fm_set),cp_failure_level,routineP,error,failure) 01469 CPPrecondition(ASSOCIATED(mos),cp_failure_level,routineP,error,failure) 01470 NULLIFY(eigenvalues, occupation_numbers) 01471 01472 IF(.NOT. failure) THEN 01473 CALL get_mo_set(mos(1)%mo_set, eigenvalues=eigenvalues,occupation_numbers=occupation_numbers,& 01474 nao=nao, nmo=nmo, homo =homo, lfomo=lfomo) 01475 01476 DO i=1,SIZE(fm_set,2) 01477 CPPrecondition(ASSOCIATED(fm_set(1,i)%matrix),cp_failure_level,routineP,error,failure) 01478 CALL cp_fm_set_all(fm_set(1,i)%matrix, 0.0_dp, error=error) 01479 CALL cp_fm_set_all(fm_work, 0.0_dp, error=error) 01480 CALL cp_dbcsr_sm_fm_multiply(op_sm(i)%matrix,all_vectors,fm_work,ncol=nstate,error=error) 01481 CALL cp_fm_gemm("T","N",1,nstate,nao,1.0_dp,excvec,& 01482 fm_work,0.0_dp, fm_set(1,i)%matrix,b_first_col=1,error=error) 01483 END DO 01484 01485 sp_em = 0.0_dp 01486 sp_ab = 0.0_dp 01487 ene_i = eigenvalues(estate) 01488 01489 DO istate = 1,nstate 01490 ene_f = all_evals(istate) 01491 01492 DO i = 1,3 01493 CALL cp_fm_get_element(fm_set(1,i)%matrix,1,istate,dip(i)) 01494 END DO 01495 IF(istate<=homo) THEN 01496 sp_em(1,istate) = ene_f - ene_i 01497 sp_em(2,istate) = dip(1) 01498 sp_em(3,istate) = dip(2) 01499 sp_em(4,istate) = dip(3) 01500 sp_em(5,istate) = dip(1)*dip(1)+dip(2)*dip(2)+dip(3)*dip(3) 01501 sp_em(6,istate) = occupation_numbers(istate) 01502 END IF 01503 IF(istate>=lfomo ) THEN 01504 i_abs = istate - lfomo + 1 01505 sp_ab(1,i_abs) = ene_f - ene_i 01506 sp_ab(2,i_abs) = dip(1) 01507 sp_ab(3,i_abs) = dip(2) 01508 sp_ab(4,i_abs) = dip(3) 01509 sp_ab(5,i_abs) = dip(1)*dip(1)+dip(2)*dip(2)+dip(3)*dip(3) 01510 IF(istate<=nmo) sp_ab(6,i_abs) = occupation_numbers(istate) 01511 END IF 01512 01513 END DO 01514 END IF 01515 01516 END SUBROUTINE spectrum_dip_vel 01517 01518 ! ***************************************************************************** 01519 SUBROUTINE calc_stogto_overlap(base_a, base_b, matrix,error) 01520 01521 TYPE(gto_basis_set_type), POINTER :: base_a, base_b 01522 REAL(dp), DIMENSION(:, :), POINTER :: matrix 01523 TYPE(cp_error_type), INTENT(inout) :: error 01524 01525 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_stogto_overlap', 01526 routineP = moduleN//':'//routineN 01527 01528 INTEGER :: iset, istat, jset, ldai, ldsab, maxcoa, maxcob, maxl, maxla, 01529 maxlb, ncoa, ncob, nseta, nsetb, nsgfa, nsgfb, sgfa, sgfb 01530 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, 01531 lb_min, npgfa, npgfb, 01532 nsgfa_set, nsgfb_set 01533 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 01534 LOGICAL :: failure 01535 REAL(dp) :: rab(3) 01536 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work 01537 REAL(dp), ALLOCATABLE, 01538 DIMENSION(:, :, :) :: ai_work 01539 REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, 01540 zeta, zetb 01541 01542 failure = .FALSE. 01543 01544 NULLIFY(la_max,la_min,lb_max,lb_min) 01545 NULLIFY(npgfa,npgfb,nsgfa_set,nsgfb_set) 01546 NULLIFY(first_sgfa,first_sgfb) 01547 NULLIFY(rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb) 01548 01549 CALL get_gto_basis_set(gto_basis_set=base_a,nsgf=nsgfa,nsgf_set=nsgfa_set,lmax=la_max,& 01550 lmin=la_min,npgf=npgfa,pgf_radius=rpgfa, & 01551 sphi=sphi_a,zet=zeta,first_sgf=first_sgfa, & 01552 maxco=maxcoa,nset=nseta, maxl=maxla) 01553 01554 CALL get_gto_basis_set(gto_basis_set=base_b,nsgf=nsgfb,nsgf_set=nsgfb_set, lmax=lb_max,& 01555 lmin=lb_min,npgf=npgfb,pgf_radius=rpgfb, & 01556 sphi=sphi_b,zet=zetb,first_sgf=first_sgfb, & 01557 maxco=maxcob,nset=nsetb,maxl=maxlb) 01558 rab = 0.0_dp 01559 ! Initialize and allocate 01560 matrix = 0.0_dp 01561 01562 ldsab = MAX(maxcoa,maxcob,nsgfa,nsgfb) 01563 maxl = MAX(maxla,maxlb) 01564 ldai = ncoset(maxl) 01565 01566 ALLOCATE(sab(ldsab,ldsab),STAT=istat) 01567 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01568 ALLOCATE(work(ldsab,ldsab),STAT=istat) 01569 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01570 ALLOCATE(ai_work(ldai,ldai,1),STAT=istat) 01571 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01572 01573 DO iset=1,nseta 01574 01575 ncoa = npgfa(iset)*ncoset(la_max(iset)) 01576 sgfa = first_sgfa(1,iset) 01577 01578 DO jset=1,nsetb 01579 01580 ncob = npgfb(jset)*ncoset(lb_max(jset)) 01581 sgfb = first_sgfb(1,jset) 01582 01583 ai_work = 0.0_dp 01584 01585 CALL overlap(la_max(iset),la_min(iset),npgfa(iset),& 01586 rpgfa(:,iset),zeta(:,iset),& 01587 lb_max(jset),lb_min(jset),npgfb(jset),& 01588 rpgfb(:,jset),zetb(:,jset),rab,0.0_dp,sab,0,.FALSE.,ai_work,ldai) 01589 CALL dgemm("N","N",ncoa,nsgfb_set(jset),ncob,& 01590 1.0_dp,sab(1,1),SIZE(sab,1),& 01591 sphi_b(1,sgfb),SIZE(sphi_b,1),& 01592 0.0_dp,work(1,1),SIZE(work,1)) 01593 CALL dgemm("T","N",nsgfa_set(iset),nsgfb_set(jset),ncoa,& 01594 1.0_dp,sphi_a(1,sgfa),SIZE(sphi_a,1),& 01595 work(1,1),SIZE(work,1),& 01596 1.0_dp,matrix(sgfa,sgfb), SIZE(matrix,1)) 01597 01598 END DO ! jset 01599 END DO ! iset 01600 01601 DEALLOCATE(sab,work,STAT=istat) 01602 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01603 01604 END SUBROUTINE calc_stogto_overlap 01605 01606 ! ***************************************************************************** 01621 SUBROUTINE cls_assign_core_states(xas_control,xas_env,localized_wfn_control,qs_env,error) 01622 01623 TYPE(xas_control_type) :: xas_control 01624 TYPE(xas_environment_type), POINTER :: xas_env 01625 TYPE(localized_wfn_control_type), 01626 POINTER :: localized_wfn_control 01627 TYPE(qs_environment_type), POINTER :: qs_env 01628 TYPE(cp_error_type), INTENT(inout) :: error 01629 01630 CHARACTER(LEN=*), PARAMETER :: routineN = 'cls_assign_core_states', 01631 routineP = moduleN//':'//routineN 01632 01633 INTEGER :: homo, i, iat, iatom, ikind, isgf, istat, istate, j, my_kind, 01634 nao, natom, nexc_atoms, nexc_search, output_unit 01635 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf 01636 INTEGER, DIMENSION(3) :: perd0 01637 INTEGER, DIMENSION(:), POINTER :: atom_of_state, 01638 mykind_of_kind, 01639 state_of_atom, 01640 state_of_mytype, type_of_state 01641 LOGICAL :: failure 01642 REAL(dp) :: component, dist, distmin, 01643 maxocc, ra(3), rac(3), rc(3) 01644 REAL(dp), DIMENSION(:), POINTER :: max_overlap, sto_state_overlap 01645 REAL(dp), DIMENSION(:, :), POINTER :: centers_wfn 01646 REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer 01647 TYPE(atomic_kind_type), POINTER :: atomic_kind 01648 TYPE(cell_type), POINTER :: cell 01649 TYPE(cp_2d_r_p_type), DIMENSION(:), 01650 POINTER :: stogto_overlap 01651 TYPE(cp_fm_type), POINTER :: mo_coeff 01652 TYPE(cp_logger_type), POINTER :: logger 01653 TYPE(mo_set_p_type), DIMENSION(:), 01654 POINTER :: mos 01655 TYPE(particle_type), DIMENSION(:), 01656 POINTER :: particle_set 01657 01658 failure = .FALSE. 01659 NULLIFY(cell,mos,particle_set) 01660 NULLIFY(atom_of_state,centers_wfn,mykind_of_kind,state_of_atom) 01661 NULLIFY(stogto_overlap,type_of_state,max_overlap) 01662 NULLIFY(state_of_mytype,type_of_state,sto_state_overlap) 01663 01664 NULLIFY(logger) 01665 logger => cp_error_get_logger(error) 01666 output_unit= cp_logger_get_default_io_unit(logger) 01667 01668 CALL get_qs_env(qs_env=qs_env, cell=cell, mos=mos, particle_set=particle_set,error=error) 01669 01670 ! The Berry operator can be used only for periodic systems 01671 ! If an isolated system is used the periodicity is overimposed 01672 perd0(1:3) = cell%perd(1:3) 01673 cell%perd(1:3) = 1 01674 01675 CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff, maxocc=maxocc, nao=nao, homo=homo) 01676 01677 CALL get_xas_env(xas_env=xas_env,& 01678 centers_wfn=centers_wfn,atom_of_state=atom_of_state,& 01679 mykind_of_kind=mykind_of_kind,& 01680 type_of_state=type_of_state, state_of_atom=state_of_atom,& 01681 stogto_overlap=stogto_overlap,nexc_atoms=nexc_atoms, nexc_search=nexc_search,error=error) 01682 01683 ! scratch array for the state 01684 ALLOCATE(vecbuffer(1,nao),STAT=istat) 01685 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01686 natom = SIZE(particle_set) 01687 01688 ALLOCATE (first_sgf(natom),STAT=istat) 01689 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01690 CALL get_particle_set(particle_set=particle_set, first_sgf=first_sgf,error=error) 01691 ALLOCATE (sto_state_overlap(nexc_search),STAT=istat) 01692 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01693 ALLOCATE (max_overlap(natom),STAT=istat) 01694 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01695 max_overlap = 0.0_dp 01696 ALLOCATE (state_of_mytype(natom),STAT=istat) 01697 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01698 state_of_mytype = 0 01699 atom_of_state = 0 01700 01701 DO istate = 1,nexc_search 01702 centers_wfn(1,istate) = localized_wfn_control%centers_set(1)%array(1,istate) 01703 centers_wfn(2,istate) = localized_wfn_control%centers_set(1)%array(2,istate) 01704 centers_wfn(3,istate) = localized_wfn_control%centers_set(1)%array(3,istate) 01705 01706 ! Assign the state to the closest atom 01707 distmin = 100.0_dp 01708 01709 DO iat = 1,nexc_atoms 01710 iatom = xas_control%exc_atoms(iat) 01711 ra(1:3) = particle_set(iatom)%r(1:3) 01712 rc(1:3) = centers_wfn(1:3,istate) 01713 rac = pbc(ra,rc,cell) 01714 dist = rac(1)*rac(1) + rac(2)*rac(2) + rac(3)*rac(3) 01715 01716 IF(dist < distmin) THEN 01717 atom_of_state(istate) = iatom 01718 distmin = dist 01719 END IF 01720 END DO 01721 IF(atom_of_state(istate) /= 0) THEN 01722 !Character of the state 01723 CALL cp_fm_get_submatrix(mo_coeff,vecbuffer,1,istate,& 01724 nao,1,transpose=.TRUE.,error=error) 01725 01726 iatom = atom_of_state(istate) 01727 01728 NULLIFY(atomic_kind) 01729 atomic_kind => particle_set(iatom)%atomic_kind 01730 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01731 kind_number=ikind) 01732 01733 my_kind = mykind_of_kind(ikind) 01734 01735 sto_state_overlap(istate) = 0.0_dp 01736 DO i = 1,SIZE(stogto_overlap(my_kind)%array,1) 01737 component = 0.0_dp 01738 DO j = 1,SIZE(stogto_overlap(my_kind)%array,2) 01739 isgf = first_sgf(iatom) + j - 1 01740 component = component + stogto_overlap(my_kind)%array(i,j)*vecbuffer(1,isgf) 01741 END DO 01742 sto_state_overlap(istate) = sto_state_overlap(istate) + & 01743 component * component 01744 END DO 01745 IF(sto_state_overlap(istate)>max_overlap(iatom)) THEN 01746 state_of_mytype(iatom) = istate 01747 max_overlap(iatom) = sto_state_overlap(istate) 01748 END IF 01749 END IF 01750 END DO ! istate 01751 01752 ! In the set of states, assign the index of the state to be excited for iatom 01753 IF (output_unit>0) THEN 01754 WRITE (UNIT=output_unit,FMT="(/,T10,A,/)")& 01755 "List the atoms to be excited and the relative of MOs index " 01756 END IF 01757 DO iat = 1,nexc_atoms 01758 iatom = xas_env%exc_atoms(iat) 01759 state_of_atom(iat) = state_of_mytype(iatom) 01760 IF (output_unit>0) THEN 01761 WRITE(UNIT=output_unit,FMT="(T10,A,I6,T32,A,I6)") 'Atom: ',iatom ,& 01762 "MO index", state_of_atom(iat) 01763 END IF 01764 IF(state_of_atom(iat)==0 .OR. state_of_atom(iat) .GT. homo) THEN 01765 CALL stop_program(routineN,moduleN,__LINE__,& 01766 "A wrong state has been selected for excitation, check the Wannier centers") 01767 END IF 01768 END DO 01769 01770 ! Set back the correct periodicity 01771 cell%perd(1:3) = perd0(1:3) 01772 01773 DEALLOCATE(vecbuffer,STAT=istat) 01774 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01775 DEALLOCATE (first_sgf,STAT=istat) 01776 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01777 DEALLOCATE (sto_state_overlap,STAT=istat) 01778 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01779 DEALLOCATE (max_overlap,STAT=istat) 01780 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01781 DEALLOCATE (state_of_mytype,STAT=istat) 01782 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 01783 01784 END SUBROUTINE cls_assign_core_states 01785 01786 END MODULE xas_methods
1.7.3