|
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 ! ***************************************************************************** 00013 MODULE xas_tp_scf 00014 USE atomic_kind_types, ONLY: atomic_kind_type,& 00015 get_atomic_kind 00016 USE cell_types, ONLY: cell_type,& 00017 pbc 00018 USE cp_array_r_utils, ONLY: cp_2d_r_p_type 00019 USE cp_control_types, ONLY: dft_control_type 00020 USE cp_dbcsr_interface, ONLY: cp_dbcsr_copy 00021 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00022 USE cp_external_control, ONLY: external_control 00023 USE cp_fm_types, ONLY: cp_fm_get_submatrix,& 00024 cp_fm_init_random,& 00025 cp_fm_set_submatrix,& 00026 cp_fm_to_fm,& 00027 cp_fm_type 00028 USE cp_output_handling, ONLY: cp_add_iter_level,& 00029 cp_iterate,& 00030 cp_p_file,& 00031 cp_print_key_finished_output,& 00032 cp_print_key_should_output,& 00033 cp_print_key_unit_nr,& 00034 cp_rm_iter_level 00035 USE cp_para_types, ONLY: cp_para_env_type 00036 USE f77_blas 00037 USE input_constants, ONLY: ot_precond_full_kinetic,& 00038 ot_precond_solver_default,& 00039 xas_dscf,& 00040 xas_scf_general 00041 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00042 section_vals_type 00043 USE kinds, ONLY: dp 00044 USE machine, ONLY: m_flush,& 00045 m_walltime 00046 USE message_passing, ONLY: mp_sync 00047 USE particle_types, ONLY: get_particle_set,& 00048 particle_type 00049 USE preconditioner, ONLY: make_preconditioner 00050 USE preconditioner_types, ONLY: destroy_preconditioner,& 00051 init_preconditioner,& 00052 preconditioner_type 00053 USE qs_charges_types, ONLY: qs_charges_type 00054 USE qs_density_mixing_types, ONLY: broyden_mixing_new_nr,& 00055 broyden_mixing_nr,& 00056 direct_mixing_nr,& 00057 gspace_mixing_nr,& 00058 multisecant_mixing_nr,& 00059 no_mixing_nr,& 00060 pulay_mixing_nr 00061 USE qs_energy_types, ONLY: qs_energy_type 00062 USE qs_environment_types, ONLY: get_qs_env,& 00063 qs_environment_type 00064 USE qs_gspace_mixing, ONLY: gspace_mixing,& 00065 self_consistency_check 00066 USE qs_ks_methods, ONLY: qs_ks_did_change,& 00067 qs_ks_update_qs_env 00068 USE qs_ks_types, ONLY: qs_ks_env_type 00069 USE qs_loc_methods, ONLY: qs_loc_driver 00070 USE qs_loc_types, ONLY: get_qs_loc_env,& 00071 localized_wfn_control_type,& 00072 qs_loc_env_new_type 00073 USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues 00074 USE qs_mo_types, ONLY: get_mo_set,& 00075 mo_set_p_type,& 00076 set_mo_occupation,& 00077 set_mo_set 00078 USE qs_ot_eigensolver, ONLY: ot_eigensolver 00079 USE qs_rho_methods, ONLY: qs_rho_update_rho 00080 USE qs_rho_types, ONLY: qs_rho_type 00081 USE qs_scf, ONLY: init_scf_run,& 00082 qs_scf_print_summary,& 00083 scf_env_cleanup,& 00084 scf_env_do_scf 00085 USE qs_scf_diagonalization, ONLY: do_general_diag 00086 USE qs_scf_methods, ONLY: scf_env_density_mixing 00087 USE qs_scf_types, ONLY: qs_scf_env_type 00088 USE scf_control_types, ONLY: scf_control_type 00089 USE timings, ONLY: timeset,& 00090 timestop 00091 USE xas_control, ONLY: xas_control_type 00092 USE xas_env_types, ONLY: get_xas_env,& 00093 set_xas_env,& 00094 xas_environment_type 00095 USE xas_restart, ONLY: xas_write_restart 00096 #include "cp_common_uses.h" 00097 00098 IMPLICIT NONE 00099 PRIVATE 00100 00101 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_tp_scf' 00102 00103 ! *** Public subroutines *** 00104 00105 PUBLIC :: xas_do_tp_scf, xes_scf_once 00106 00107 CONTAINS 00108 00109 ! ***************************************************************************** 00122 SUBROUTINE xas_do_tp_scf (dft_control,xas_env,iatom,scf_env,qs_env,& 00123 xas_section,scf_section,converged,should_stop,error) 00124 00125 TYPE(dft_control_type), POINTER :: dft_control 00126 TYPE(xas_environment_type), POINTER :: xas_env 00127 INTEGER, INTENT(IN) :: iatom 00128 TYPE(qs_scf_env_type), POINTER :: scf_env 00129 TYPE(qs_environment_type), POINTER :: qs_env 00130 TYPE(section_vals_type), POINTER :: xas_section, scf_section 00131 LOGICAL, INTENT(OUT) :: converged, should_stop 00132 TYPE(cp_error_type), INTENT(inout) :: error 00133 00134 CHARACTER(LEN=*), PARAMETER :: routineN = 'xas_do_tp_scf', 00135 routineP = moduleN//':'//routineN 00136 00137 INTEGER :: handle, handle2, ispin, 00138 iter_count, output_unit 00139 LOGICAL :: diis_step, energy_only, 00140 exit_loop, failure, gapw, 00141 use_jacobi 00142 REAL(KIND=dp) :: t1, t2 00143 TYPE(atomic_kind_type), DIMENSION(:), 00144 POINTER :: atomic_kind_set 00145 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00146 POINTER :: matrix_ks, matrix_s 00147 TYPE(cp_logger_type), POINTER :: logger 00148 TYPE(cp_para_env_type), POINTER :: para_env 00149 TYPE(mo_set_p_type), DIMENSION(:), 00150 POINTER :: mos 00151 TYPE(qs_charges_type), POINTER :: qs_charges 00152 TYPE(qs_energy_type), POINTER :: energy 00153 TYPE(qs_ks_env_type), POINTER :: ks_env 00154 TYPE(qs_rho_type), POINTER :: rho 00155 TYPE(scf_control_type), POINTER :: scf_control 00156 TYPE(xas_control_type), POINTER :: xas_control 00157 00158 CALL timeset(routineN,handle) 00159 NULLIFY(xas_control,matrix_s,matrix_ks,para_env) 00160 NULLIFY(rho,energy,scf_control,logger, ks_env,mos,atomic_kind_set) 00161 NULLIFY(qs_charges) 00162 00163 logger => cp_error_get_logger(error) 00164 t1 = m_walltime() 00165 failure=.FALSE. 00166 converged = .TRUE. 00167 00168 CPPrecondition(ASSOCIATED(xas_env),cp_failure_level,routineP,error,failure) 00169 CPPrecondition(xas_env%ref_count>0,cp_failure_level,routineP,error,failure) 00170 CPPrecondition(ASSOCIATED(scf_env),cp_failure_level,routineP,error,failure) 00171 CPPrecondition(scf_env%ref_count>0,cp_failure_level,routineP,error,failure) 00172 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 00173 CPPrecondition(qs_env%ref_count>0,cp_failure_level,routineP,error,failure) 00174 00175 CALL get_qs_env(qs_env=qs_env,& 00176 atomic_kind_set=atomic_kind_set,& 00177 scf_control=scf_control,& 00178 matrix_s=matrix_s,energy=energy,& 00179 qs_charges=qs_charges,& 00180 ks_env=ks_env,para_env=para_env,& 00181 error=error) 00182 00183 energy_only = .FALSE. 00184 output_unit=cp_print_key_unit_nr(logger,xas_section,"PRINT%PROGRAM_RUN_INFO",& 00185 extension=".xasLog",error=error) 00186 IF (output_unit>0) THEN 00187 WRITE (UNIT=output_unit,FMT="(/,/,T2,A)") "XAS_TP_SCF WAVEFUNCTION OPTIMIZATION" 00188 END IF 00189 00190 ! GAPW method must be used 00191 gapw = dft_control%qs_control%gapw 00192 CPPrecondition(gapw,cp_failure_level,routineP,error,failure) 00193 xas_control => dft_control%xas_control 00194 00195 CALL cp_add_iter_level(logger%iter_info,"XAS_SCF",error=error) 00196 00197 CALL get_qs_env(qs_env,matrix_ks=matrix_ks,rho=rho,mos=mos,error=error) 00198 00199 iter_count = 0 00200 diis_step = .FALSE. 00201 use_jacobi = .FALSE. 00202 00203 IF (output_unit>0) THEN 00204 WRITE (UNIT=output_unit,& 00205 FMT="(/,T3,A,T12,A,T31,A,T40,A,T60,A,T75,A/,T3,A)")& 00206 "Step","Update method","Time","Convergence","Total energy","Change",& 00207 REPEAT("-",77) 00208 END IF 00209 00210 ! *** SCF loop *** 00211 00212 energy%tot_old = 0.0_dp 00213 scf_loop: DO 00214 CALL timeset(routineN//"_inner_loop",handle2) 00215 00216 exit_loop = .FALSE. 00217 IF (output_unit > 0) CALL m_flush(output_unit) 00218 00219 iter_count = iter_count + 1 00220 CALL cp_iterate(logger%iter_info,last=.FALSE.,iter_nr=iter_count,error=error) 00221 00222 ! ** here qs_env%rho%rho_r and qs_env%rho%rho_g should be up to date 00223 00224 CALL qs_ks_update_qs_env(ks_env,qs_env=qs_env,error=error,& 00225 calculate_forces=.FALSE.,just_energy=energy_only) 00226 00227 scf_env%mixing_store%alpha = xas_env%mixing_store%alpha 00228 scf_env%mixing_store%iter_method = xas_env%mixing_store%iter_method 00229 00230 SELECT CASE (xas_control%scf_method) 00231 CASE DEFAULT 00232 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00233 routineP,"unknown scf method method for core level spectroscopy"//& 00234 cp_to_string(xas_control%scf_method),error,failure) 00235 00236 CASE(xas_scf_general) ! diagonalisation (default) 00237 scf_env%iter_count=iter_count 00238 CALL do_general_diag(scf_env,mos,matrix_ks,& 00239 matrix_s,scf_control,scf_section, & 00240 diis_step,use_jacobi,xas_env,error) 00241 00242 END SELECT 00243 00244 SELECT CASE(xas_env%mixing_method) 00245 CASE(direct_mixing_nr) 00246 CALL scf_env_density_mixing(scf_env%p_mix_new,& 00247 xas_env%mixing_store, rho%rho_ao, para_env, scf_env%iter_delta, scf_env%iter_count, & 00248 diis=diis_step, error=error) 00249 CASE(gspace_mixing_nr,pulay_mixing_nr,broyden_mixing_nr,& 00250 broyden_mixing_new_nr,multisecant_mixing_nr) 00251 ! Compute the difference p_out-p_in 00252 CALL self_consistency_check(qs_env%rho%rho_ao,scf_env%p_delta,para_env,scf_env%p_mix_new,& 00253 delta= scf_env%iter_delta, error=error) 00254 CASE(no_mixing_nr) 00255 CASE DEFAULT 00256 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00257 routineP,"unknown scf mixing method: "//& 00258 cp_to_string(xas_env%mixing_method),error,failure) 00259 END SELECT 00260 00261 t2 = m_walltime() 00262 00263 IF (output_unit>0.and.scf_env%print_iter_line) THEN 00264 WRITE (UNIT=output_unit,& 00265 FMT="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,F14.8,1X,F20.10,1X,ES9.2)")& 00266 iter_count,TRIM(scf_env%iter_method),& 00267 scf_env%iter_param,t2 - t1,scf_env%iter_delta,energy%total,& 00268 energy%total-energy%tot_old 00269 00270 END IF 00271 energy%tot_old = energy%total 00272 00273 ! ** convergence check 00274 CALL external_control(should_stop,"XASSCF",target_time=qs_env%target_time,& 00275 start_time=qs_env%start_time,error=error) 00276 IF (scf_env%iter_delta < xas_control%eps_scf) THEN 00277 IF (output_unit>0) THEN 00278 WRITE(UNIT=output_unit,FMT="(/,T3,A,I5,A/)")& 00279 "*** SCF run converged in ",iter_count," steps ***" 00280 END IF 00281 exit_loop = .TRUE. 00282 ELSE IF (should_stop.OR. iter_count == xas_control%max_scf) THEN 00283 IF (output_unit>0) THEN 00284 WRITE(UNIT=output_unit,FMT="(/,T3,A,/)")& 00285 "*** SCF run NOT converged ***" 00286 END IF 00287 converged = .FALSE. 00288 exit_loop = .TRUE. 00289 END IF 00290 ! *** Exit if we have finished with the SCF inner loop *** 00291 IF (exit_loop) THEN 00292 ! now, print out energies and charges corresponding to the obtained wfn 00293 ! (this actually is not 100% consistent at this point)! 00294 CALL qs_scf_print_summary(output_unit,rho,qs_charges,energy,scf_env%nelectron, & 00295 dft_control,qs_env%qmmm,qs_env,.TRUE.,.FALSE.,error) 00296 CALL cp_iterate(logger%iter_info,last=.TRUE.,iter_nr=iter_count,error=error) 00297 END IF 00298 00299 ! ** Write restart file ** 00300 CALL xas_write_restart(xas_env, xas_section, qs_env, xas_control%xas_method,& 00301 iatom,error=error) 00302 00303 IF (exit_loop) THEN 00304 CALL timestop(handle2) 00305 EXIT scf_loop 00306 END IF 00307 00308 IF (.NOT.BTEST(cp_print_key_should_output(logger%iter_info,& 00309 xas_section,"PRINT%ITERATION_INFO/TIME_CUMUL",error=error),cp_p_file)) t1 = m_walltime() 00310 00311 ! *** mixing methods have the new density matrix in p_mix_new 00312 IF (xas_env%mixing_method > 0) THEN 00313 DO ispin=1,dft_control%nspins 00314 CALL cp_dbcsr_copy(rho%rho_ao(ispin)%matrix,scf_env%p_mix_new(ispin)%matrix,& 00315 error=error) 00316 END DO 00317 ENDIF 00318 00319 ! ** update qs_env%rho 00320 CALL qs_rho_update_rho(rho, qs_env=qs_env, error=error) 00321 IF(xas_env%mixing_method>=gspace_mixing_nr) THEN 00322 CALL gspace_mixing(qs_env, scf_env, xas_env%mixing_store, rho, qs_env%para_env, error=error) 00323 END IF 00324 00325 CALL qs_ks_did_change(ks_env,rho_changed=.TRUE.,error=error) 00326 CALL timestop(handle2) 00327 00328 END DO scf_loop 00329 00330 IF (output_unit>0) THEN 00331 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T55,F25.14))")& 00332 "Ionization potential of the excited atom: ",xas_env%IP_energy 00333 CALL m_flush(output_unit) 00334 END IF 00335 00336 CALL mp_sync(qs_env%para_env%group) 00337 CALL qs_ks_did_change(ks_env,rho_changed=.TRUE.,error=error) 00338 00339 CALL cls_prepare_states(xas_control,xas_env,qs_env,iatom,xas_section,output_unit,error=error) 00340 00341 CALL mp_sync(qs_env%para_env%group) 00342 00343 CALL cp_print_key_finished_output(output_unit,logger,xas_section,& 00344 "PRINT%PROGRAM_RUN_INFO", error=error) 00345 CALL cp_rm_iter_level(logger%iter_info,"XAS_SCF",error=error) 00346 00347 CALL timestop(handle) 00348 00349 END SUBROUTINE xas_do_tp_scf 00350 ! ***************************************************************************** 00362 SUBROUTINE cls_prepare_states(xas_control,xas_env,qs_env,iatom,xas_section,output_unit,error) 00363 00364 TYPE(xas_control_type) :: xas_control 00365 TYPE(xas_environment_type), POINTER :: xas_env 00366 TYPE(qs_environment_type), POINTER :: qs_env 00367 INTEGER, INTENT(IN) :: iatom 00368 TYPE(section_vals_type), POINTER :: xas_section 00369 INTEGER, INTENT(IN) :: output_unit 00370 TYPE(cp_error_type), INTENT(inout) :: error 00371 00372 CHARACTER(LEN=*), PARAMETER :: routineN = 'cls_prepare_states', 00373 routineP = moduleN//':'//routineN 00374 00375 INTEGER :: handle, i, ikind, isgf, ispin, istat, istate, j, my_kind, 00376 my_state, nao, natom, nexc_search, nmo, nvirtual2, uno_iter, xas_estate 00377 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf 00378 INTEGER, DIMENSION(:), POINTER :: mykind_of_kind 00379 LOGICAL :: failure 00380 REAL(dp), DIMENSION(:, :), POINTER :: centers_wfn 00381 REAL(KIND=dp) :: component, dist, max_overlap, 00382 ra(3), rac(3), rc(3), 00383 sto_state_overlap, uno_eps 00384 REAL(KIND=dp), DIMENSION(:), POINTER :: all_evals, eigenvalues, 00385 uno_evals 00386 REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer, vecbuffer2 00387 TYPE(atomic_kind_type), POINTER :: atomic_kind 00388 TYPE(cell_type), POINTER :: cell 00389 TYPE(cp_2d_r_p_type), DIMENSION(:), 00390 POINTER :: stogto_overlap 00391 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00392 POINTER :: matrix_ks, matrix_s 00393 TYPE(cp_fm_type), POINTER :: all_vectors, excvec_coeff, 00394 excvec_overlap, mo_coeff, 00395 uno_orbs 00396 TYPE(dft_control_type), POINTER :: dft_control 00397 TYPE(localized_wfn_control_type), 00398 POINTER :: localized_wfn_control 00399 TYPE(mo_set_p_type), DIMENSION(:), 00400 POINTER :: mos 00401 TYPE(particle_type), DIMENSION(:), 00402 POINTER :: particle_set 00403 TYPE(preconditioner_type), POINTER :: local_preconditioner 00404 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 00405 TYPE(section_vals_type), POINTER :: loc_section, print_loc_section 00406 00407 CALL timeset(routineN,handle) 00408 00409 NULLIFY(atomic_kind,dft_control,matrix_s,matrix_ks) 00410 NULLIFY(cell,particle_set,local_preconditioner,vecbuffer,vecbuffer2) 00411 NULLIFY(dft_control,loc_section,mos,mo_coeff,eigenvalues) 00412 NULLIFY(centers_wfn,mykind_of_kind,qs_loc_env,localized_wfn_control,stogto_overlap) 00413 NULLIFY(all_evals,all_vectors,excvec_coeff,excvec_overlap,uno_evals) 00414 00415 failure=.FALSE. 00416 00417 CPPrecondition(ASSOCIATED(xas_env),cp_failure_level,routineP,error,failure) 00418 00419 CALL get_qs_env(qs_env=qs_env,& 00420 cell=cell,dft_control=dft_control,& 00421 matrix_ks=matrix_ks, matrix_s=matrix_s,mos=mos,& 00422 particle_set=particle_set,error=error) 00423 00424 ! Some elements from the xas_env 00425 CALL get_xas_env(xas_env=xas_env,& 00426 all_vectors=all_vectors,all_evals=all_evals,& 00427 excvec_coeff=excvec_coeff,& 00428 nvirtual2=nvirtual2,xas_estate=xas_estate,& 00429 excvec_overlap=excvec_overlap,nexc_search=nexc_search,error=error) 00430 CPPrecondition(ASSOCIATED(excvec_overlap),cp_failure_level,routineP,error,failure) 00431 00432 CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff,nao=nao,& 00433 eigenvalues=eigenvalues) 00434 00435 ALLOCATE(vecbuffer(1,nao),STAT=istat) 00436 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00437 vecbuffer = 0.0_dp 00438 ALLOCATE(vecbuffer2(1,nao),STAT=istat) 00439 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00440 vecbuffer2 = 0.0_dp 00441 natom=SIZE(particle_set,1) 00442 ALLOCATE (first_sgf(natom),STAT=istat) 00443 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00444 CALL get_particle_set(particle_set=particle_set, first_sgf=first_sgf,& 00445 error=error) 00446 ALLOCATE(centers_wfn(3,nexc_search),STAT=istat) 00447 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00448 centers_wfn=0.0_dp 00449 00450 ! Possible only for emission only 00451 IF(.NOT.xas_control%scf_method==xas_scf_general) THEN 00452 IF (output_unit>0) THEN 00453 WRITE (UNIT=output_unit,FMT="(/,/,T2,A)") " Eigenstates are derived "//& 00454 "from the MOs optimized by OT. Follows localization of the core states"//& 00455 " to identify the excited orital. " 00456 END IF 00457 00458 CALL get_xas_env(xas_env=xas_env, & 00459 mykind_of_kind=mykind_of_kind, qs_loc_env=qs_loc_env,& 00460 stogto_overlap=stogto_overlap,error=error) 00461 CALL get_qs_loc_env(qs_loc_env=qs_loc_env,& 00462 localized_wfn_control=localized_wfn_control,error=error) 00463 loc_section => section_vals_get_subs_vals(xas_section,"LOCALIZE",error=error) 00464 print_loc_section => section_vals_get_subs_vals(loc_section,"PRINT",error=error) 00465 CALL qs_loc_driver(qs_env,qs_loc_env,loc_section,print_loc_section,myspin=1,error=error) 00466 ra(1:3) = particle_set(iatom)%r(1:3) 00467 00468 NULLIFY(atomic_kind) 00469 atomic_kind => particle_set(iatom)%atomic_kind 00470 CALL get_atomic_kind(atomic_kind=atomic_kind,& 00471 kind_number=ikind) 00472 my_kind = mykind_of_kind(ikind) 00473 00474 00475 CALL cp_fm_get_submatrix(mo_coeff,vecbuffer2,1,my_state,& 00476 nao,1,transpose=.TRUE.,error=error) 00477 00478 ! Rotate the wfn to get the eigenstate of the KS hamiltonian 00479 ! Only ispin=1 should be needed 00480 DO ispin=1,dft_control%nspins 00481 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff,nmo=nmo,& 00482 eigenvalues=eigenvalues) 00483 CALL calculate_subspace_eigenvalues(mo_coeff,& 00484 matrix_ks(ispin)%matrix,eigenvalues, & 00485 do_rotation=.TRUE.,error=error) 00486 END DO ! ispin 00487 00488 !Search for the core state to be excited 00489 max_overlap = 0.0_dp 00490 DO istate = 1,nexc_search 00491 centers_wfn(1,istate) = localized_wfn_control%centers_set(1)%array(1,istate) 00492 centers_wfn(2,istate) = localized_wfn_control%centers_set(1)%array(2,istate) 00493 centers_wfn(3,istate) = localized_wfn_control%centers_set(1)%array(3,istate) 00494 00495 rc(1:3) = centers_wfn(1:3,istate) 00496 rac = pbc(ra,rc,cell) 00497 dist = rac(1)*rac(1) + rac(2)*rac(2) + rac(3)*rac(3) 00498 00499 IF(dist < 1.0_dp) THEN 00500 CALL cp_fm_get_submatrix(mo_coeff,vecbuffer,1,istate,& 00501 nao,1,transpose=.TRUE.,error=error) 00502 sto_state_overlap=0.0_dp 00503 DO i = 1,SIZE(stogto_overlap(my_kind)%array,1) 00504 component = 0.0_dp 00505 DO j = 1,SIZE(stogto_overlap(my_kind)%array,2) 00506 isgf = first_sgf(iatom) + j - 1 00507 component = component + stogto_overlap(my_kind)%array(i,j)*vecbuffer(1,isgf) 00508 END DO ! j size 00509 sto_state_overlap = sto_state_overlap + & 00510 component * component 00511 END DO ! i size 00512 IF(sto_state_overlap .GT. max_overlap) THEN 00513 max_overlap = sto_state_overlap 00514 my_state = istate 00515 END IF 00516 END IF 00517 xas_estate = my_state 00518 END DO ! istate 00519 00520 CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff) 00521 CALL cp_fm_get_submatrix(mo_coeff,vecbuffer,1,xas_estate,& 00522 nao,1,transpose=.TRUE.,error=error) 00523 CALL cp_fm_set_submatrix(excvec_coeff,vecbuffer2,1,1,& 00524 nao,1,transpose=.TRUE.,error=error) 00525 ! 00526 END IF 00527 00528 CALL mp_sync(qs_env%para_env%group) 00529 !Calculate the virtual states from the KS matrix matrix_ks(1) 00530 IF(nvirtual2 .GT. 0) THEN 00531 NULLIFY(mo_coeff) 00532 CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff,nmo=nmo) 00533 IF (output_unit>0) THEN 00534 WRITE (UNIT=output_unit,FMT="(/,/,T2,A,I5,A,I6,A)") " Calculation of ", nvirtual2,& 00535 " additional virtual states of the subspace complementary to the "//& 00536 " lowest ", nmo, " states" 00537 END IF 00538 00539 NULLIFY(uno_orbs,uno_evals,local_preconditioner) 00540 ALLOCATE(local_preconditioner,STAT=istat) 00541 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00542 CALL init_preconditioner(local_preconditioner,para_env=qs_env%para_env,& 00543 blacs_env=qs_env%blacs_env,error=error) 00544 00545 CALL make_preconditioner(local_preconditioner, & 00546 precon_type=ot_precond_full_kinetic, & 00547 solver_type=ot_precond_solver_default, & 00548 matrix_h=matrix_ks(1)%matrix, & 00549 matrix_s=matrix_s(1)%matrix,& 00550 matrix_t=qs_env%kinetic(1)%matrix, & 00551 convert_precond_to_dbcsr=.TRUE.,& 00552 mo_set=mos(1)%mo_set,energy_gap=0.2_dp,error=error) 00553 00554 CALL get_xas_env(xas_env=xas_env,unoccupied_orbs=uno_orbs,& 00555 unoccupied_evals=uno_evals,unoccupied_eps=uno_eps,unoccupied_max_iter=uno_iter,error=error) 00556 CALL cp_fm_init_random(uno_orbs,nvirtual2,error=error) 00557 00558 CALL ot_eigensolver(matrix_h=matrix_ks(1)%matrix,matrix_s=matrix_s(1)%matrix, & 00559 matrix_c_fm=uno_orbs,matrix_orthogonal_space_fm=mo_coeff,& 00560 preconditioner=local_preconditioner,eps_gradient=uno_eps,& 00561 iter_max=uno_iter,size_ortho_space=nmo,error=error) 00562 00563 CALL calculate_subspace_eigenvalues(uno_orbs,matrix_ks(1)%matrix,& 00564 uno_evals,do_rotation=.TRUE.,error=error) 00565 CALL destroy_preconditioner(local_preconditioner,error=error) 00566 00567 DEALLOCATE(local_preconditioner,STAT=istat) 00568 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00569 END IF 00570 00571 CALL mp_sync(qs_env%para_env%group) 00572 ! Prapare arrays for the calculation of the spectrum 00573 IF(.NOT. xas_control%xas_method == xas_dscf) THEN 00574 ! Copy the final vectors in the array 00575 NULLIFY(all_vectors,all_evals) 00576 CALL get_xas_env(xas_env=xas_env,all_vectors=all_vectors,& 00577 all_evals=all_evals,error=error) 00578 CALL get_mo_set(mos(1)%mo_set, eigenvalues=eigenvalues,mo_coeff=mo_coeff,& 00579 nmo=nmo) 00580 00581 CALL cp_fm_to_fm(mo_coeff,all_vectors,ncol=nmo,& 00582 source_start=1,target_start=1) 00583 DO istate = 1,nmo 00584 all_evals(istate) = eigenvalues(istate) 00585 ENDDO 00586 IF(nvirtual2 .GT. 0) THEN 00587 CALL cp_fm_to_fm(uno_orbs,all_vectors,ncol=nvirtual2,& 00588 source_start=1,target_start=1+nmo) 00589 DO istate = 1,nvirtual2 00590 all_evals(istate+nmo) = uno_evals(istate) 00591 END DO 00592 END IF 00593 END IF 00594 00595 DEALLOCATE(vecbuffer,STAT=istat) 00596 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00597 DEALLOCATE(vecbuffer2,STAT=istat) 00598 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00599 DEALLOCATE(centers_wfn,first_sgf,STAT=istat) 00600 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00601 00602 CALL timestop(handle) 00603 00604 END SUBROUTINE cls_prepare_states 00605 00606 ! ***************************************************************************** 00617 SUBROUTINE xes_scf_once(qs_env,xas_env,scf_section,converged,should_stop,error) 00618 00619 TYPE(qs_environment_type), POINTER :: qs_env 00620 TYPE(xas_environment_type), POINTER :: xas_env 00621 TYPE(section_vals_type), POINTER :: scf_section 00622 LOGICAL, INTENT(OUT) :: converged, should_stop 00623 TYPE(cp_error_type), INTENT(inout) :: error 00624 00625 CHARACTER(LEN=*), PARAMETER :: routineN = 'xes_scf_once', 00626 routineP = moduleN//':'//routineN 00627 00628 INTEGER :: handle, ispin, istate, 00629 nelectron, nmo, nvirtual, 00630 nvirtual2, output_unit 00631 LOGICAL :: failure 00632 REAL(KIND=dp), DIMENSION(:), POINTER :: all_evals, eigenvalues 00633 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00634 POINTER :: matrix_ks 00635 TYPE(cp_fm_type), POINTER :: all_vectors, mo_coeff 00636 TYPE(cp_logger_type), POINTER :: logger 00637 TYPE(cp_para_env_type), POINTER :: para_env 00638 TYPE(dft_control_type), POINTER :: dft_control 00639 TYPE(mo_set_p_type), DIMENSION(:), 00640 POINTER :: mos 00641 TYPE(qs_scf_env_type), POINTER :: scf_env 00642 TYPE(scf_control_type), POINTER :: scf_control 00643 00644 failure = .FALSE. 00645 NULLIFY(dft_control,scf_control,matrix_ks,mos,para_env) 00646 NULLIFY(logger) 00647 logger => cp_error_get_logger(error) 00648 output_unit= cp_logger_get_default_io_unit(logger) 00649 00650 CALL timeset(routineN,handle) 00651 00652 CPPrecondition(ASSOCIATED(xas_env),cp_failure_level,routineP,error,failure) 00653 IF(.NOT. failure) THEN 00654 CALL get_qs_env(qs_env=qs_env,dft_control=dft_control, scf_control=scf_control,& 00655 matrix_ks=matrix_ks,mos=mos,para_env=para_env,error=error) 00656 00657 CALL get_mo_set(mos(1)%mo_set,nelectron=nelectron) 00658 nelectron = nelectron - 1 00659 CALL set_mo_set(mos(1)%mo_set,nelectron=nelectron,error=error) 00660 CALL set_mo_occupation(mo_set=mos(1)%mo_set,smear=xas_env%smear,& 00661 error=error) 00662 00663 00664 NULLIFY(scf_env) 00665 CALL get_qs_env(qs_env,scf_env=scf_env,error=error) 00666 CALL init_scf_run(scf_env=scf_env,qs_env=qs_env,& 00667 scf_section=scf_section, do_xas_tp=.TRUE., error=error) 00668 00669 CALL scf_env_do_scf(scf_env=scf_env, qs_env=qs_env, & 00670 converged=converged, should_stop=should_stop, error=error) 00671 CALL scf_env_cleanup(scf_env,qs_env=qs_env,error=error) 00672 00673 ! The eigenstate of the KS Hamiltonian are nedeed 00674 NULLIFY(mo_coeff,eigenvalues) 00675 IF(scf_control%use_ot) THEN 00676 IF (output_unit>0) THEN 00677 WRITE (UNIT=output_unit,FMT="(/,T10,A,/)")& 00678 "Get eigenstates and eigenvalues from ground state MOs" 00679 END IF 00680 DO ispin = 1,SIZE(mos) 00681 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=mo_coeff,& 00682 eigenvalues=eigenvalues) 00683 CALL calculate_subspace_eigenvalues(mo_coeff,& 00684 matrix_ks(ispin)%matrix,eigenvalues, & 00685 do_rotation=.TRUE.,error=error) 00686 END DO 00687 END IF 00688 NULLIFY(all_vectors,all_evals) 00689 CALL get_xas_env(xas_env=xas_env,all_vectors=all_vectors,& 00690 all_evals=all_evals,nvirtual2=nvirtual2,error=error) 00691 CALL get_mo_set(mos(1)%mo_set, eigenvalues=eigenvalues,mo_coeff=mo_coeff, nmo=nmo) 00692 00693 CALL cp_fm_to_fm(mo_coeff,all_vectors,ncol=nmo,& 00694 source_start=1,target_start=1) 00695 DO istate = 1,nmo 00696 all_evals(istate) = eigenvalues(istate) 00697 ENDDO 00698 00699 IF(nvirtual2/=0) THEN 00700 IF (output_unit>0) THEN 00701 WRITE (UNIT=output_unit,FMT="(/,T10,A,/)")& 00702 "WARNING: for this XES calculation additional unoccupied MOs are not needed" 00703 END IF 00704 nvirtual2=0 00705 nvirtual = nmo 00706 CALL set_xas_env(xas_env=xas_env,nvirtual=nvirtual,nvirtual2=nvirtual2,error=error) 00707 END IF 00708 00709 nelectron = nelectron + 1 00710 CALL set_mo_set(mos(1)%mo_set,nelectron=nelectron,error=error) 00711 END IF ! failure 00712 CALL timestop(handle) 00713 00714 END SUBROUTINE xes_scf_once 00715 00716 END MODULE xas_tp_scf
1.7.3