|
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 qs_scf_post_gpw 00013 USE admm_types, ONLY: admm_type 00014 USE admm_utils, ONLY: admm_correct_for_eigenvalues,& 00015 admm_uncorrect_for_eigenvalues 00016 USE ai_onecenter, ONLY: sg_overlap 00017 USE atom_kind_orbitals, ONLY: calculate_atomic_density 00018 USE atomic_kind_types, ONLY: atomic_kind_type,& 00019 get_atomic_kind 00020 USE cell_types, ONLY: cell_type,& 00021 get_cell,& 00022 pbc 00023 USE cp_array_r_utils, ONLY: cp_1d_r_p_type 00024 USE cp_control_types, ONLY: dft_control_type,& 00025 tddfpt_control_type 00026 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 00027 copy_fm_to_dbcsr 00028 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix 00029 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type,& 00030 cp_dbcsr_type 00031 USE cp_ddapc_util, ONLY: get_ddapc 00032 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add 00033 USE cp_fm_diag, ONLY: choose_eigv_solver 00034 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00035 cp_fm_struct_release,& 00036 cp_fm_struct_type 00037 USE cp_fm_types, ONLY: cp_fm_create,& 00038 cp_fm_get_info,& 00039 cp_fm_init_random,& 00040 cp_fm_p_type,& 00041 cp_fm_release,& 00042 cp_fm_to_fm,& 00043 cp_fm_type 00044 USE cp_fm_vect, ONLY: cp_fm_vect_dealloc 00045 USE cp_output_handling, ONLY: cp_iter_string,& 00046 cp_p_file,& 00047 cp_print_key_finished_output,& 00048 cp_print_key_should_output,& 00049 cp_print_key_unit_nr 00050 USE cp_para_types, ONLY: cp_para_env_type 00051 USE cp_result_methods, ONLY: cp_results_erase,& 00052 get_results,& 00053 put_results 00054 USE cp_subsys_types, ONLY: cp_subsys_get,& 00055 cp_subsys_type 00056 USE et_coupling_types, ONLY: set_et_coupling_type 00057 USE f77_blas 00058 USE hfx_ri_gemopt, ONLY: geminal_optimize 00059 USE input_constants, ONLY: do_loc_both,& 00060 do_loc_homo,& 00061 do_loc_lumo,& 00062 ot_precond_full_all,& 00063 use_orb_basis_set 00064 USE input_section_types, ONLY: section_get_ival,& 00065 section_get_ivals,& 00066 section_get_lval,& 00067 section_get_rval,& 00068 section_vals_get,& 00069 section_vals_get_subs_vals,& 00070 section_vals_type,& 00071 section_vals_val_get 00072 USE kinds, ONLY: default_path_length,& 00073 default_string_length,& 00074 dp 00075 USE lapack, ONLY: lapack_sgesv 00076 USE mathconstants, ONLY: pi,& 00077 twopi 00078 USE molecular_states, ONLY: construct_molecular_states 00079 USE molecule_types_new, ONLY: molecule_type 00080 USE moments_utils, ONLY: get_reference_point 00081 USE orbital_pointers, ONLY: indso 00082 USE particle_list_types, ONLY: particle_list_type 00083 USE particle_types, ONLY: particle_type 00084 USE physcon, ONLY: angstrom,& 00085 debye,& 00086 evolt 00087 USE population_analyses, ONLY: lowdin_population_analysis,& 00088 mulliken_population_analysis 00089 USE preconditioner_types, ONLY: preconditioner_type 00090 USE pw_env_types, ONLY: pw_env_get,& 00091 pw_env_type 00092 USE pw_grids, ONLY: get_pw_grid_info 00093 USE pw_methods, ONLY: pw_axpy,& 00094 pw_copy,& 00095 pw_derive,& 00096 pw_integrate_function,& 00097 pw_scale,& 00098 pw_transfer,& 00099 pw_zero 00100 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00101 pw_pool_give_back_pw,& 00102 pw_pool_p_type,& 00103 pw_pool_type 00104 USE pw_types, ONLY: COMPLEXDATA1D,& 00105 REALDATA3D,& 00106 REALSPACE,& 00107 RECIPROCALSPACE,& 00108 pw_p_type,& 00109 pw_type 00110 USE qs_collocate_density, ONLY: calculate_rho_elec,& 00111 calculate_wavefunction 00112 USE qs_conductivity, ONLY: optical_conductivity 00113 USE qs_core_energies, ONLY: calculate_ecore 00114 USE qs_electric_field_gradient, ONLY: qs_efg_calc 00115 USE qs_energy_types, ONLY: qs_energy_type 00116 USE qs_environment_types, ONLY: get_qs_env,& 00117 qs_environment_type,& 00118 set_qs_env 00119 USE qs_epr_hyp, ONLY: qs_epr_hyp_calc 00120 USE qs_grid_atom, ONLY: grid_atom_type 00121 USE qs_ks_methods, ONLY: qs_ks_did_change,& 00122 qs_ks_update_qs_env 00123 USE qs_loc_methods, ONLY: qs_loc_driver 00124 USE qs_loc_molecules, ONLY: wfc_to_molecule 00125 USE qs_loc_types, ONLY: qs_loc_env_create,& 00126 qs_loc_env_destroy,& 00127 qs_loc_env_new_type 00128 USE qs_loc_utils, ONLY: loc_write_restart,& 00129 qs_loc_control_init,& 00130 qs_loc_env_init,& 00131 qs_loc_init,& 00132 retain_history 00133 USE qs_mo_methods, ONLY: calculate_orthonormality,& 00134 calculate_subspace_eigenvalues,& 00135 make_mo_eig 00136 USE qs_mo_types, ONLY: deallocate_mo_set,& 00137 duplicate_mo_set,& 00138 get_mo_set,& 00139 mo_set_p_type,& 00140 set_mo_occupation,& 00141 write_mo_set 00142 USE qs_moments, ONLY: qs_moment_berry_phase,& 00143 qs_moment_locop 00144 USE qs_ot_eigensolver, ONLY: ot_eigensolver 00145 USE qs_pdos, ONLY: calculate_projected_dos 00146 USE qs_resp, ONLY: resp_fit 00147 USE qs_rho_atom_types, ONLY: rho_atom_type 00148 USE qs_rho_types, ONLY: qs_rho_type 00149 USE qs_scf_types, ONLY: ot_method_nr,& 00150 qs_scf_env_type,& 00151 special_diag_method_nr 00152 USE realspace_grid_cube, ONLY: pw_to_cube 00153 USE s_square_methods, ONLY: compute_s_square 00154 USE scf_control_types, ONLY: scf_control_type 00155 USE stm_images, ONLY: th_stm_image 00156 USE termination, ONLY: stop_program 00157 USE timings, ONLY: timeset,& 00158 timestop 00159 USE wannier_states, ONLY: construct_wannier_states 00160 USE wannier_states_types, ONLY: wannier_centres_type 00161 USE xray_diffraction, ONLY: calculate_rhotot_elec_gspace,& 00162 xray_diffraction_spectrum 00163 #include "cp_common_uses.h" 00164 00165 IMPLICIT NONE 00166 PRIVATE 00167 00168 ! Global parameters 00169 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_post_gpw' 00170 PUBLIC :: scf_post_calculation_gpw,& 00171 qs_scf_post_occ_cubes,& 00172 qs_scf_post_unocc_cubes,& 00173 qs_scf_post_moments,& 00174 qs_scf_post_efg,& 00175 write_available_results,& 00176 write_mo_free_results 00177 00178 CONTAINS 00179 00180 ! ***************************************************************************** 00198 SUBROUTINE scf_post_calculation_gpw(dft_section, scf_env,qs_env, error) 00199 00200 TYPE(section_vals_type), POINTER :: dft_section 00201 TYPE(qs_scf_env_type), POINTER :: scf_env 00202 TYPE(qs_environment_type), POINTER :: qs_env 00203 TYPE(cp_error_type), INTENT(inout) :: error 00204 00205 CHARACTER(len=*), PARAMETER :: routineN = 'scf_post_calculation_gpw', 00206 routineP = moduleN//':'//routineN 00207 00208 INTEGER :: handle, homo, ispin, istat, min_lumos, n_rep, nhomo, nlumo, 00209 nlumo_stm, nlumo_tddft, nlumos, nmo, output_unit, unit_nr 00210 INTEGER, DIMENSION(:, :, :), POINTER :: marked_states 00211 LOGICAL :: check_write, compute_lumos, do_homo, do_mo_cubes, do_stm, 00212 do_wannier_cubes, failure, has_homo, has_lumo, loc_explicit, 00213 loc_print_explicit, my_localized_wfn, p_loc, p_loc_homo, p_loc_lumo 00214 REAL(dp) :: e_kin, e_kinAtt 00215 REAL(KIND=dp) :: gap, homo_lumo(2,2) 00216 REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues 00217 TYPE(admm_type), POINTER :: admm_env 00218 TYPE(atomic_kind_type), DIMENSION(:), 00219 POINTER :: atomic_kind_set 00220 TYPE(cp_1d_r_p_type), DIMENSION(:), 00221 POINTER :: occupied_evals, 00222 unoccupied_evals 00223 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00224 POINTER :: kinetic_m, ks_rmpv, 00225 ks_rmpv_aux_fit, matrix_s, 00226 mo_derivs 00227 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: homo_localized, 00228 lumo_localized, lumo_ptr, mo_loc_history, occupied_orbs, unoccupied_orbs 00229 TYPE(cp_fm_type), POINTER :: mo_coeff 00230 TYPE(cp_logger_type), POINTER :: logger 00231 TYPE(cp_para_env_type), POINTER :: para_env 00232 TYPE(cp_subsys_type), POINTER :: subsys 00233 TYPE(dft_control_type), POINTER :: dft_control 00234 TYPE(mo_set_p_type), DIMENSION(:), 00235 POINTER :: mos 00236 TYPE(molecule_type), POINTER :: molecule_set( : ) 00237 TYPE(particle_list_type), POINTER :: particles 00238 TYPE(particle_type), DIMENSION(:), 00239 POINTER :: particle_set 00240 TYPE(pw_env_type), POINTER :: pw_env 00241 TYPE(pw_p_type) :: wf_g, wf_r 00242 TYPE(pw_pool_p_type), DIMENSION(:), 00243 POINTER :: pw_pools 00244 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00245 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env_homo, 00246 qs_loc_env_lumo 00247 TYPE(qs_rho_type), POINTER :: rho 00248 TYPE(scf_control_type), POINTER :: scf_control 00249 TYPE(section_vals_type), POINTER :: input, loc_print_section, 00250 localize_section, print_key, 00251 stm_section 00252 TYPE(tddfpt_control_type), POINTER :: tddfpt_control 00253 00254 CALL timeset(routineN,handle) 00255 00256 ! Writes the data that is already available in qs_env 00257 para_env=>qs_env%para_env 00258 CALL write_available_results(qs_env,scf_env,error) 00259 00260 failure=.FALSE. 00261 my_localized_wfn = .FALSE. 00262 NULLIFY(admm_env, dft_control, pw_env, auxbas_pw_pool, pw_pools, mos, rho, & 00263 mo_coeff, ks_rmpv, ks_rmpv_aux_fit, matrix_s, qs_loc_env_homo,qs_loc_env_lumo, scf_control, & 00264 unoccupied_orbs, mo_eigenvalues, unoccupied_evals, molecule_set, mo_derivs,& 00265 tddfpt_control, subsys, particles, input, print_key, kinetic_m,marked_states) 00266 NULLIFY(homo_localized, lumo_localized,lumo_ptr) 00267 00268 has_homo=.FALSE. 00269 has_lumo=.FALSE. 00270 p_loc = .FALSE. 00271 p_loc_homo = .FALSE. 00272 p_loc_lumo= .FALSE. 00273 00274 00275 logger => cp_error_get_logger(error) 00276 output_unit= cp_logger_get_default_io_unit(logger) 00277 00278 CPPrecondition(ASSOCIATED(scf_env),cp_failure_level,routineP,error,failure) 00279 CPPrecondition(scf_env%ref_count>0,cp_failure_level,routineP,error,failure) 00280 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 00281 ! Here we start with data that needs a postprocessing... 00282 IF (.NOT. failure) THEN 00283 CALL get_qs_env(qs_env,dft_control=dft_control,molecule_set=molecule_set, & 00284 mos=mos,rho=rho,matrix_ks=ks_rmpv,scf_control=scf_control,matrix_s=matrix_s, & 00285 input=input, subsys=subsys, tddfpt_control=tddfpt_control, pw_env=pw_env,& 00286 particle_set=particle_set, atomic_kind_set=atomic_kind_set, error=error) 00287 CALL cp_subsys_get(subsys,particles=particles,error=error) 00288 00289 ! **** the kinetic energy 00290 IF (cp_print_key_should_output(logger%iter_info,input,& 00291 "DFT%PRINT%KINETIC_ENERGY",error=error)/=0) THEN 00292 CALL get_qs_env(qs_env,kinetic=kinetic_m,error=error) 00293 CPPrecondition(ASSOCIATED(kinetic_m),cp_failure_level,routineP,error,failure) 00294 CPPrecondition(ASSOCIATED(kinetic_m(1)%matrix),cp_failure_level,routineP,error,failure) 00295 e_kin=0.0_dp 00296 DO ispin=1,dft_control%nspins 00297 CALL calculate_ecore(h=kinetic_m(1)%matrix,& 00298 p=rho%rho_ao(ispin)%matrix,& 00299 ecore=e_kinAtt,& 00300 error=error) 00301 e_kin=e_kin+e_kinAtt 00302 END DO 00303 unit_nr = cp_print_key_unit_nr(logger,input,"DFT%PRINT%KINETIC_ENERGY",& 00304 extension=".Log",error=error) 00305 IF (unit_nr>0) THEN 00306 WRITE (unit_nr,'(T3,A,T55,F25.14)') "Electronic kinetic energy:",e_kin 00307 ENDIF 00308 CALL cp_print_key_finished_output(unit_nr,logger,input,& 00309 "DFT%PRINT%KINETIC_ENERGY", error=error) 00310 END IF 00311 ! Compute Atomic Charges 00312 CALL qs_scf_post_charges(input, logger, qs_env, rho, matrix_s, mos, error=error) 00313 00314 ! Moments of charge distribution 00315 CALL qs_scf_post_moments(input, logger, qs_env, output_unit, error=error) 00316 00317 ! Determine if we need to computer properties using the localized centers 00318 localize_section => section_vals_get_subs_vals(dft_section,"LOCALIZE",error=error) 00319 loc_print_section => section_vals_get_subs_vals(localize_section,"PRINT",error=error) 00320 CALL section_vals_get(localize_section, explicit=loc_explicit, error=error) 00321 CALL section_vals_get(loc_print_section, explicit=loc_print_explicit, error=error) 00322 00323 ! Print_keys controlled by localization 00324 IF(loc_print_explicit) THEN 00325 print_key => section_vals_get_subs_vals(loc_print_section,"MOLECULAR_DIPOLES",error=error) 00326 p_loc=BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file) 00327 print_key => section_vals_get_subs_vals(loc_print_section,"TOTAL_DIPOLE",error=error) 00328 p_loc=p_loc.OR.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file) 00329 print_key => section_vals_get_subs_vals(loc_print_section,"WANNIER_CENTERS",error=error) 00330 p_loc=p_loc.OR.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file) 00331 print_key => section_vals_get_subs_vals(loc_print_section,"WANNIER_SPREADS",error=error) 00332 p_loc=p_loc.OR.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file) 00333 print_key => section_vals_get_subs_vals(loc_print_section,"WANNIER_CUBES",error=error) 00334 p_loc=p_loc.OR.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file) 00335 print_key => section_vals_get_subs_vals(loc_print_section,"MOLECULAR_STATES",error=error) 00336 p_loc=p_loc.OR.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file) 00337 ELSE 00338 p_loc=.FALSE. 00339 END IF 00340 IF(loc_explicit) THEN 00341 p_loc_homo=(section_get_ival(localize_section,"STATES",error=error)==do_loc_homo.OR.& 00342 section_get_ival(localize_section,"STATES",error=error)==do_loc_both).AND.p_loc 00343 p_loc_lumo=(section_get_ival(localize_section,"STATES",error=error)==do_loc_lumo.OR.& 00344 section_get_ival(localize_section,"STATES",error=error)==do_loc_both).AND.p_loc 00345 CALL section_vals_val_get(localize_section,"LIST_UNOCCUPIED", n_rep_val=n_rep,error=error) 00346 ELSE 00347 p_loc_homo=.FALSE. 00348 p_loc_lumo=.FALSE. 00349 n_rep=0 00350 END IF 00351 00352 IF(n_rep==0.AND.p_loc_lumo)THEN 00353 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 00354 "No LIST_UNOCCUPIED was specified, therefore localization of unoccupied states will be skipped!"//& 00355 CPSourceFileRef,& 00356 only_ionode=.TRUE.) 00357 p_loc_lumo=.FALSE. 00358 END IF 00359 print_key => section_vals_get_subs_vals(loc_print_section,"WANNIER_STATES",error=error) 00360 p_loc=p_loc.OR.BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file) 00361 00362 00363 ! Control for STM 00364 stm_section => section_vals_get_subs_vals(input,"DFT%PRINT%STM", error=error) 00365 CALL section_vals_get(stm_section,explicit=do_stm,error=error) 00366 ! BTEST(cp_print_key_should_output(logger%iter_info,stm_section,error=error),cp_p_file) 00367 nlumo_stm = 0 00368 IF(do_stm) nlumo_stm = section_get_ival(stm_section,"NLUMO",error=error) 00369 00370 ! check for CUBES (MOs and WANNIERS) 00371 do_mo_cubes=BTEST(cp_print_key_should_output(logger%iter_info,dft_section,"PRINT%MO_CUBES",& 00372 error=error),cp_p_file) 00373 IF(loc_print_explicit) THEN 00374 do_wannier_cubes=BTEST(cp_print_key_should_output(logger%iter_info,loc_print_section,& 00375 "WANNIER_CUBES",error=error),cp_p_file) 00376 ELSE 00377 do_wannier_cubes=.FALSE. 00378 END IF 00379 nlumo=section_get_ival(dft_section,"PRINT%MO_CUBES%NLUMO",error=error) 00380 nhomo=section_get_ival(dft_section,"PRINT%MO_CUBES%NHOMO",error=error) 00381 nlumo_tddft = 0 00382 IF(dft_control%do_tddfpt_calculation) THEN 00383 nlumo_tddft=section_get_ival(dft_section,"TDDFPT%NLUMO",error=error) 00384 END IF 00385 00386 ! Setup the grids needed to compute a wavefunction given a vector.. 00387 IF ( ( ( do_mo_cubes .OR. do_wannier_cubes ).AND. (nlumo /= 0 .OR. nhomo /= 0 )) .OR. p_loc ) THEN 00388 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,& 00389 pw_pools=pw_pools,error=error) 00390 CALL pw_pool_create_pw(auxbas_pw_pool,wf_r%pw,& 00391 use_data = REALDATA3D,& 00392 in_space = REALSPACE, error=error) 00393 CALL pw_pool_create_pw(auxbas_pw_pool,wf_g%pw,& 00394 use_data = COMPLEXDATA1D,& 00395 in_space = RECIPROCALSPACE, error=error) 00396 END IF 00397 00398 ! Makes the MOs eigenstates, computes eigenvalues, write cubes 00399 IF ((do_mo_cubes .AND. nhomo /= 0) .OR. do_stm .OR. dft_control%do_tddfpt_calculation ) THEN 00400 IF (dft_control%restricted) THEN 00401 IF (output_unit>0) WRITE(output_unit,*) & 00402 " Unclear how we define MOs in the restricted case ... skipping" 00403 ELSE 00404 CALL get_qs_env(qs_env,mo_derivs=mo_derivs,error=error) 00405 IF(dft_control%do_admm) THEN 00406 CALL get_qs_env(qs_env,admm_env=admm_env, matrix_ks_aux_fit=ks_rmpv_aux_fit, error=error) 00407 CALL make_mo_eig(mos,dft_control%nspins,ks_rmpv,scf_control,mo_derivs,admm_env=admm_env,& 00408 ks_aux=ks_rmpv_aux_fit,error=error) 00409 ELSE 00410 CALL make_mo_eig(mos,dft_control%nspins,ks_rmpv,scf_control,mo_derivs,error=error) 00411 END IF 00412 END IF 00413 DO ispin=1,dft_control%nspins 00414 CALL get_mo_set(mo_set=mos(ispin)%mo_set, eigenvalues=mo_eigenvalues,homo=homo) 00415 homo_lumo(ispin,1)=mo_eigenvalues(homo) 00416 END DO 00417 has_homo=.TRUE. 00418 END IF 00419 IF (do_mo_cubes .AND. nhomo /= 0) THEN 00420 DO ispin=1,dft_control%nspins 00421 ! Prints the cube files of OCCUPIED ORBITALS 00422 CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff, & 00423 eigenvalues=mo_eigenvalues,homo=homo,nmo=nmo) 00424 CALL qs_scf_post_occ_cubes(input, dft_section, dft_control, logger, qs_env,& 00425 mo_coeff, wf_g, wf_r, particles, homo, ispin, error=error) 00426 END DO 00427 ENDIF 00428 00429 ! Initialize the localization environment, needed e.g. for wannier functions and molecular states 00430 ! Gets localization info for the occupied orbs 00431 ! - Possibly gets wannier functions 00432 ! - Possibly gets molecular states 00433 IF (p_loc_homo) THEN 00434 IF (qs_env%dft_control%restricted) THEN 00435 IF (output_unit>0) WRITE(output_unit,*) & 00436 " Unclear how we define MOs / localization in the restricted case ... skipping" 00437 ELSE 00438 ALLOCATE(occupied_orbs(dft_control%nspins),stat=istat) 00439 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00440 ALLOCATE(occupied_evals(dft_control%nspins),stat=istat) 00441 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00442 ALLOCATE(homo_localized(dft_control%nspins),stat=istat) 00443 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00444 DO ispin=1,dft_control%nspins 00445 CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff, & 00446 eigenvalues=mo_eigenvalues) 00447 occupied_orbs(ispin)%matrix=>mo_coeff 00448 occupied_evals(ispin)%array=>mo_eigenvalues 00449 CALL cp_fm_create(homo_localized(ispin)%matrix,occupied_orbs(ispin)%matrix%matrix_struct,error=error) 00450 CALL cp_fm_to_fm(occupied_orbs(ispin)%matrix,homo_localized(ispin)%matrix,error=error) 00451 END DO 00452 00453 CALL get_qs_env(qs_env,mo_loc_history=mo_loc_history,error=error) 00454 do_homo=.TRUE. 00455 00456 CALL qs_loc_env_create(qs_loc_env_homo,error=error) 00457 CALL qs_loc_control_init(qs_loc_env_homo,localize_section,do_homo=do_homo,error=error) 00458 CALL qs_loc_init(qs_env,qs_loc_env_homo,localize_section,homo_localized,do_homo,& 00459 do_mo_cubes,mo_loc_history=mo_loc_history,error=error) 00460 CALL get_localization_info(qs_env,qs_loc_env_homo,localize_section,homo_localized,& 00461 wf_r, wf_g,particles,occupied_orbs,occupied_evals,marked_states,error=error) 00462 00463 !retain the homo_localized for future use 00464 IF (qs_loc_env_homo%localized_wfn_control%use_history) THEN 00465 CALL retain_history(mo_loc_history,homo_localized,error) 00466 CALL set_qs_env(qs_env,mo_loc_history=mo_loc_history,error=error) 00467 ENDIF 00468 00469 !write restart for localization of occupied orbitals 00470 CALL loc_write_restart(qs_env,qs_loc_env_homo,loc_print_section,mos,& 00471 homo_localized, do_homo, error=error) 00472 CALL cp_fm_vect_dealloc(homo_localized,error) 00473 DEALLOCATE(occupied_orbs,stat=istat) 00474 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00475 DEALLOCATE(occupied_evals,stat=istat) 00476 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00477 ! Print Total Dipole if the localization has been performed 00478 IF (qs_loc_env_homo%do_localize) THEN 00479 CALL qs_scf_post_loc_dip(input, dft_control, qs_loc_env_homo, logger, qs_env, error) 00480 END IF 00481 END IF 00482 ENDIF 00483 00484 ! Gets the lumos, and eigenvalues for the lumos, and localize them if requested 00485 nlumo = MAX(nlumo, nlumo_stm,nlumo_tddft) 00486 compute_lumos=(do_mo_cubes .OR.dft_control%do_tddfpt_calculation .OR. do_stm ) .AND. nlumo .NE. 0 00487 compute_lumos=compute_lumos.OR.p_loc_lumo 00488 00489 IF (compute_lumos) THEN 00490 check_write=.TRUE. 00491 min_lumos=nlumo 00492 IF(nlumo==0)check_write=.FALSE. 00493 IF(p_loc_lumo)THEN 00494 do_homo=.FALSE. 00495 CALL qs_loc_env_create(qs_loc_env_lumo,error=error) 00496 CALL qs_loc_control_init(qs_loc_env_lumo,localize_section,do_homo=do_homo,error=error) 00497 min_lumos=MAX(MAXVAL(qs_loc_env_lumo%localized_wfn_control%loc_states(:,:)),nlumo) 00498 END IF 00499 00500 ALLOCATE(unoccupied_orbs(dft_control%nspins),STAT=istat) 00501 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00502 ALLOCATE(unoccupied_evals(dft_control%nspins),STAT=istat) 00503 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00504 CALL make_lumo(qs_env,scf_env,unoccupied_orbs,unoccupied_evals,min_lumos,nlumos,error) 00505 lumo_ptr=>unoccupied_orbs 00506 DO ispin=1,dft_control%nspins 00507 has_lumo=.TRUE. 00508 homo_lumo(ispin,2)=unoccupied_evals(ispin)%array(1) 00509 CALL get_mo_set(mo_set=mos(ispin)%mo_set, homo=homo) 00510 IF(check_write)THEN 00511 IF(p_loc_lumo.AND.nlumo.NE.-1)nlumos=MIN(nlumo,nlumos) 00512 ! Prints the cube files of UNOCCUPIED ORBITALS 00513 CALL qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env,& 00514 unoccupied_orbs, wf_g, wf_r, particles, nlumos, homo, ispin, error) 00515 END IF 00516 END DO 00517 00518 ! Save the info for tddfpt calculation 00519 IF (dft_control%do_tddfpt_calculation) THEN 00520 ALLOCATE(tddfpt_control%lumos_eigenvalues(nlumos,dft_control%nspins),stat=istat) 00521 DO ispin=1, dft_control%nspins 00522 tddfpt_control%lumos_eigenvalues(1:nlumos,ispin) = & 00523 unoccupied_evals(ispin)%array(1:nlumos) 00524 END DO 00525 tddfpt_control%lumos => unoccupied_orbs 00526 END IF 00527 00528 IF(p_loc_lumo)THEN 00529 ALLOCATE(lumo_localized(dft_control%nspins),STAT=istat) 00530 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00531 DO ispin=1,dft_control%nspins 00532 CALL cp_fm_create(lumo_localized(ispin)%matrix,unoccupied_orbs(ispin)%matrix%matrix_struct,error=error) 00533 CALL cp_fm_to_fm(unoccupied_orbs(ispin)%matrix,lumo_localized(ispin)%matrix,error=error) 00534 END DO 00535 CALL qs_loc_init(qs_env,qs_loc_env_lumo,localize_section,lumo_localized,do_homo,do_mo_cubes,& 00536 evals=unoccupied_evals, error=error) 00537 CALL qs_loc_env_init(qs_loc_env_lumo,qs_loc_env_lumo%localized_wfn_control,qs_env,& 00538 loc_coeff=unoccupied_orbs,error=error) 00539 CALL get_localization_info(qs_env,qs_loc_env_lumo,localize_section, & 00540 lumo_localized, wf_r, wf_g,particles,& 00541 unoccupied_orbs,unoccupied_evals,marked_states,error=error) 00542 CALL loc_write_restart(qs_env,qs_loc_env_lumo,loc_print_section,mos,homo_localized, do_homo,& 00543 evals=unoccupied_evals, error=error) 00544 lumo_ptr=>lumo_localized 00545 END IF 00546 ENDIF 00547 00548 IF (has_homo .AND. has_lumo) THEN 00549 IF (output_unit>0) WRITE(output_unit,*) " " 00550 DO ispin=1,dft_control%nspins 00551 IF(.NOT. scf_control%smear%do_smear) THEN 00552 gap = homo_lumo(ispin,2)-homo_lumo(ispin,1) 00553 IF (output_unit>0) WRITE(output_unit,'(T2,A,F12.6)')& 00554 "HOMO - LUMO gap [eV] :",gap*evolt 00555 END IF 00556 ENDDO 00557 ENDIF 00558 00559 ! Deallocate grids needed to compute wavefunctions 00560 IF ( ( ( do_mo_cubes .OR. do_wannier_cubes ).AND. (nlumo /= 0 .OR. nhomo /= 0 )) .OR. p_loc ) THEN 00561 CALL pw_pool_give_back_pw(auxbas_pw_pool,wf_r%pw, error=error) 00562 CALL pw_pool_give_back_pw(auxbas_pw_pool,wf_g%pw, error=error) 00563 END IF 00564 00565 ! Destroy the localization environment 00566 IF (p_loc_homo) CALL qs_loc_env_destroy(qs_loc_env_homo, error=error) 00567 IF (p_loc_lumo) CALL qs_loc_env_destroy(qs_loc_env_lumo, error=error) 00568 00569 00570 !stm images 00571 IF(do_stm)THEN 00572 CALL th_stm_image(qs_env, stm_section, particles, unoccupied_orbs, & 00573 unoccupied_evals, error) 00574 END IF 00575 00576 ! generate a mix of wfns, and write to a restart 00577 CALL wfn_mix(mos, particle_set, dft_section, atomic_kind_set, & 00578 lumo_ptr, scf_env, matrix_s, output_unit, marked_states,error) 00579 IF(p_loc_lumo)CALL cp_fm_vect_dealloc(lumo_localized,error) 00580 IF(ASSOCIATED(marked_states))THEN 00581 DEALLOCATE(marked_states,stat=istat) 00582 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00583 END IF 00584 00585 IF(compute_lumos) THEN 00586 ! Compute the optical conductivity (needs homos and lumos) 00587 CALL qs_scf_post_optc(input, dft_section, dft_control, logger, qs_env,& 00588 unoccupied_orbs, unoccupied_evals, output_unit, error) 00589 END IF 00590 00591 ! This is just a deallocation for printing MO_CUBES or TDDFPT 00592 IF(compute_lumos) THEN 00593 DO ispin=1,dft_control%nspins 00594 DEALLOCATE(unoccupied_evals(ispin)%array,stat=istat) 00595 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00596 IF (.NOT.dft_control%do_tddfpt_calculation) & 00597 CALL cp_fm_release(unoccupied_orbs(ispin)%matrix,error=error) 00598 ENDDO 00599 DEALLOCATE(unoccupied_evals,stat=istat) 00600 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00601 IF (.NOT.dft_control%do_tddfpt_calculation) THEN 00602 DEALLOCATE(unoccupied_orbs,stat=istat) 00603 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00604 END IF 00605 ENDIF 00606 00607 ! Print coherent X-ray diffraction spectrum 00608 CALL qs_scf_post_xray(input, dft_section, logger, qs_env, output_unit, error) 00609 00610 ! Calculation of Electric Field Gradients 00611 CALL qs_scf_post_efg(input, logger, qs_env, error) 00612 00613 ! Calculation of ET 00614 CALL qs_scf_post_et(input, qs_env, dft_control, error) 00615 00616 ! Calculation of EPR Hyperfine Coupling Tensors 00617 CALL qs_scf_post_epr(input, logger, qs_env, error) 00618 00619 ! Calculation of properties needed for BASIS_MOLOPT optimizations 00620 CALL qs_scf_post_molopt(input, logger, qs_env, error) 00621 00622 ! Optimize Geminals 00623 CALL qs_scf_post_gemopt(input, logger, qs_env, error) 00624 00625 ! Calculate ELF 00626 CALL qs_scf_post_elf(input, logger, qs_env, error) 00627 END IF 00628 CALL timestop(handle) 00629 00630 00631 END SUBROUTINE scf_post_calculation_gpw 00632 00633 ! ***************************************************************************** 00636 SUBROUTINE get_localization_info(qs_env,qs_loc_env,loc_section,mo_local,& 00637 wf_r, wf_g,particles,coeff,evals,marked_states,error) 00638 00639 TYPE(qs_environment_type), POINTER :: qs_env 00640 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 00641 TYPE(section_vals_type), POINTER :: loc_section 00642 TYPE(cp_fm_p_type), DIMENSION(:), 00643 POINTER :: mo_local 00644 TYPE(pw_p_type) :: wf_r, wf_g 00645 TYPE(particle_list_type), POINTER :: particles 00646 TYPE(cp_fm_p_type), DIMENSION(:), 00647 POINTER :: coeff 00648 TYPE(cp_1d_r_p_type), DIMENSION(:), 00649 POINTER :: evals 00650 INTEGER, DIMENSION(:, :, :), POINTER :: marked_states 00651 TYPE(cp_error_type), INTENT(inout) :: error 00652 00653 CHARACTER(len=*), PARAMETER :: routineN = 'get_localization_info', 00654 routineP = moduleN//':'//routineN 00655 00656 INTEGER :: handle, ispin, mystate, ns, 00657 output_unit, stat 00658 INTEGER, DIMENSION(:), POINTER :: lstates, marked_states_spin 00659 LOGICAL :: do_homo, failure 00660 REAL(KIND=dp), DIMENSION(:, :), POINTER :: scenter 00661 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00662 POINTER :: ks_rmpv, matrix_s 00663 TYPE(cp_logger_type), POINTER :: logger 00664 TYPE(dft_control_type), POINTER :: dft_control 00665 TYPE(mo_set_p_type), DIMENSION(:), 00666 POINTER :: mos 00667 TYPE(molecule_type), POINTER :: molecule_set( : ) 00668 TYPE(section_vals_type), POINTER :: loc_print_section 00669 TYPE(wannier_centres_type), 00670 DIMENSION(:), POINTER :: wc 00671 00672 CALL timeset(routineN,handle) 00673 failure=.FALSE. 00674 NULLIFY(mos, ks_rmpv, dft_control, loc_print_section, marked_states_spin,& 00675 matrix_s, molecule_set, scenter,wc) 00676 CALL get_qs_env(qs_env,mos=mos,matrix_ks=ks_rmpv,molecule_set=molecule_set,& 00677 dft_control=dft_control,matrix_s=matrix_s,error=error) 00678 logger => cp_error_get_logger(error) 00679 output_unit= cp_logger_get_default_io_unit(logger) 00680 loc_print_section => section_vals_get_subs_vals(loc_section,"PRINT",error=error) 00681 do_homo=qs_loc_env%localized_wfn_control%do_homo 00682 IF (BTEST(cp_print_key_should_output(logger%iter_info,loc_print_section,& 00683 "WANNIER_STATES",error=error),cp_p_file)) THEN 00684 CALL get_qs_env(qs_env=qs_env,WannierCentres=wc,error=error) 00685 IF (.NOT. ASSOCIATED(wc)) THEN 00686 ALLOCATE(wc(dft_control%nspins),stat=stat) 00687 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00688 ENDIF 00689 ENDIF 00690 DO ispin=1,dft_control%nspins 00691 ! 00692 IF(do_homo)THEN 00693 qs_loc_env%tag_mo="HOMO" 00694 ELSE 00695 qs_loc_env%tag_mo="LUMO" 00696 END IF 00697 00698 IF (qs_loc_env%do_localize) THEN 00699 ! Do the Real localization.. 00700 IF (output_unit>0.AND.do_homo) WRITE(output_unit,"(/,T2,A,I3)")& 00701 "LOCALIZATION| Computing localization properties "//& 00702 "for OCCUPIED ORBITALS. Spin:",ispin 00703 IF (output_unit>0.AND.(.NOT.do_homo)) WRITE(output_unit,"(/,T2,A,I3)")& 00704 "LOCALIZATION| Computing localization properties "//& 00705 "for UNOCCUPIED ORBITALS. Spin:",ispin 00706 00707 scenter => qs_loc_env%localized_wfn_control%centers_set(ispin)%array 00708 00709 CALL qs_loc_driver(qs_env,qs_loc_env,loc_section,loc_print_section,& 00710 myspin=ispin,ext_mo_coeff=mo_local(ispin)%matrix,error=error) 00711 00712 ! maps wfc to molecules, and compute the molecular dipoles if required 00713 IF (( BTEST(cp_print_key_should_output(logger%iter_info,loc_print_section,& 00714 "MOLECULAR_DIPOLES",error=error),cp_p_file) .OR. & 00715 BTEST(cp_print_key_should_output(logger%iter_info,loc_print_section,& 00716 "MOLECULAR_STATES",error=error),cp_p_file))) THEN 00717 CALL wfc_to_molecule(qs_env, qs_loc_env, loc_print_section, scenter,& 00718 molecule_set, dft_control%nspins, error) 00719 END IF 00720 00721 ! Compute the wannier states 00722 IF (BTEST(cp_print_key_should_output(logger%iter_info,loc_print_section,& 00723 "WANNIER_STATES",error=error),cp_p_file)) THEN 00724 ns=SIZE(qs_loc_env%localized_wfn_control%loc_states,1) 00725 IF (.NOT. ASSOCIATED(wc(ispin)%centres)) THEN 00726 ALLOCATE(wc(ispin)%WannierHamDiag(ns),stat=stat) 00727 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00728 ALLOCATE(wc(ispin)%centres(3,ns),stat=stat) 00729 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00730 ENDIF 00731 00732 wc(ispin)%centres(:,:)=scenter(1+(ispin-1)*3:ispin*3,:) 00733 lstates => qs_loc_env%localized_wfn_control%loc_states(:,ispin) 00734 CALL construct_wannier_states(molecule_set, mo_local(ispin)%matrix,& 00735 ks_rmpv(ispin)%matrix, qs_env, loc_print_section=loc_print_section,& 00736 WannierCentres=wc(ispin),ns=ns,states=lstates, error=error) 00737 ENDIF 00738 ! Compute the molecular states 00739 IF ( BTEST(cp_print_key_should_output(logger%iter_info,loc_print_section,& 00740 "MOLECULAR_STATES",error=error),cp_p_file)) THEN 00741 CALL construct_molecular_states(molecule_set, mo_local(ispin)%matrix, coeff(ispin)%matrix,& 00742 evals(ispin)%array,ks_rmpv(ispin)%matrix, matrix_s(1)%matrix, qs_env, wf_r, wf_g, & 00743 loc_print_section=loc_print_section, particles=particles, tag=TRIM(qs_loc_env%tag_mo),& 00744 marked_states=marked_states_spin,error=error) 00745 IF(ASSOCIATED(marked_states_spin))THEN 00746 IF(.NOT.ASSOCIATED(marked_states))THEN 00747 ALLOCATE(marked_states(SIZE(marked_states_spin),dft_control%nspins,2),stat=stat) 00748 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00749 END IF 00750 mystate=1 00751 IF(qs_loc_env%tag_mo=="LUMO")mystate=2 00752 marked_states(:,ispin,mystate)=marked_states_spin(:) 00753 DEALLOCATE(marked_states_spin,stat=stat) 00754 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 00755 END IF 00756 ENDIF 00757 END IF 00758 ENDDO 00759 IF (BTEST(cp_print_key_should_output(logger%iter_info,loc_print_section,& 00760 "WANNIER_STATES",error=error),cp_p_file)) THEN 00761 CALL set_qs_env(qs_env=qs_env,WannierCentres=wc,error=error) 00762 ENDIF 00763 00764 CALL timestop(handle) 00765 00766 END SUBROUTINE get_localization_info 00767 00768 ! ***************************************************************************** 00771 SUBROUTINE make_lumo(qs_env,scf_env,unoccupied_orbs,unoccupied_evals,nlumo,nlumos,error) 00772 00773 TYPE(qs_environment_type), POINTER :: qs_env 00774 TYPE(qs_scf_env_type), POINTER :: scf_env 00775 TYPE(cp_fm_p_type), DIMENSION(:), 00776 POINTER :: unoccupied_orbs 00777 TYPE(cp_1d_r_p_type), DIMENSION(:), 00778 POINTER :: unoccupied_evals 00779 INTEGER :: nlumo 00780 INTEGER, INTENT(OUT) :: nlumos 00781 TYPE(cp_error_type), INTENT(inout) :: error 00782 00783 CHARACTER(len=*), PARAMETER :: routineN = 'make_lumo', 00784 routineP = moduleN//':'//routineN 00785 00786 INTEGER :: handle, homo, ispin, istat, 00787 n, nao, nmo, output_unit 00788 LOGICAL :: failure 00789 TYPE(admm_type), POINTER :: admm_env 00790 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00791 POINTER :: ks_rmpv, ks_rmpv_aux_fit, 00792 matrix_s 00793 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 00794 TYPE(cp_fm_type), POINTER :: mo_coeff 00795 TYPE(cp_logger_type), POINTER :: logger 00796 TYPE(dft_control_type), POINTER :: dft_control 00797 TYPE(mo_set_p_type), DIMENSION(:), 00798 POINTER :: mos 00799 TYPE(preconditioner_type), POINTER :: local_preconditioner 00800 TYPE(scf_control_type), POINTER :: scf_control 00801 00802 CALL timeset(routineN,handle) 00803 00804 NULLIFY(mos,ks_rmpv,scf_control,dft_control,admm_env,ks_rmpv_aux_fit) 00805 CALL get_qs_env(qs_env,mos=mos,matrix_ks=ks_rmpv,scf_control=scf_control,& 00806 dft_control=dft_control,matrix_s=matrix_s,admm_env=admm_env,& 00807 matrix_ks_aux_fit=ks_rmpv_aux_fit, error=error) 00808 logger => cp_error_get_logger(error) 00809 output_unit= cp_logger_get_default_io_unit(logger) 00810 00811 DO ispin=1,dft_control%nspins 00812 NULLIFY(unoccupied_orbs(ispin)%matrix) 00813 NULLIFY(unoccupied_evals(ispin)%array) 00814 ! Always write eigenvalues 00815 IF (output_unit>0) WRITE(output_unit,*) " " 00816 IF (output_unit>0) WRITE(output_unit,*) " Lowest Eigenvalues of the unoccupied subspace spin ",ispin 00817 IF (output_unit>0) WRITE(output_unit,FMT='(1X,A)') "-----------------------------------------------------" 00818 CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff,homo=homo,nao=nao,nmo=nmo) 00819 CALL cp_fm_get_info(mo_coeff, nrow_global=n,error=error) 00820 nlumos=MAX(1,MIN(nlumo,nao-homo)) 00821 IF (nlumo==-1) nlumos=nao-homo 00822 ALLOCATE(unoccupied_evals(ispin)%array(nlumos),stat=istat) 00823 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 00824 CALL cp_fm_struct_create(fm_struct_tmp,para_env=qs_env%para_env,context=qs_env%blacs_env, & 00825 nrow_global=n,ncol_global=nlumos,error=error) 00826 CALL cp_fm_create(unoccupied_orbs(ispin)%matrix, fm_struct_tmp,name="lumos",error=error) 00827 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00828 CALL cp_fm_init_random(unoccupied_orbs(ispin)%matrix,nlumos,error=error) 00829 00830 ! the full_all preconditioner makes not much sense for lumos search 00831 NULLIFY(local_preconditioner) 00832 IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN 00833 local_preconditioner=>scf_env%ot_preconditioner(1)%preconditioner 00834 ! this one can for sure not be right (as it has to match a given C0) 00835 IF (local_preconditioner%in_use == ot_precond_full_all) THEN 00836 NULLIFY(local_preconditioner) 00837 ENDIF 00838 ENDIF 00839 00840 ! ** If we do ADMM, we add have to modify the kohn-sham matrix 00841 IF( dft_control%do_admm ) THEN 00842 CALL admm_correct_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix, & 00843 ks_rmpv_aux_fit(ispin)%matrix, error) 00844 END IF 00845 00846 CALL ot_eigensolver(matrix_h=ks_rmpv(ispin)%matrix,matrix_s=matrix_s(1)%matrix, & 00847 matrix_c_fm=unoccupied_orbs(ispin)%matrix, & 00848 matrix_orthogonal_space_fm=mo_coeff, & 00849 eps_gradient=scf_control%eps_lumos, & 00850 preconditioner=local_preconditioner, & 00851 iter_max=scf_control%max_iter_lumos,& 00852 size_ortho_space=nmo,error=error) 00853 00854 CALL calculate_subspace_eigenvalues(unoccupied_orbs(ispin)%matrix,ks_rmpv(ispin)%matrix,& 00855 unoccupied_evals(ispin)%array, scr=output_unit, & 00856 ionode=output_unit>0,error=error) 00857 00858 00859 ! ** If we do ADMM, we restore the original kohn-sham matrix 00860 IF( dft_control%do_admm ) THEN 00861 CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix, & 00862 ks_rmpv_aux_fit(ispin)%matrix, error) 00863 END IF 00864 00865 END DO 00866 00867 CALL timestop(handle) 00868 00869 END SUBROUTINE make_lumo 00870 00871 ! ***************************************************************************** 00874 SUBROUTINE wfn_mix(mos, particle_set, dft_section, atomic_kind_set, & 00875 unoccupied_orbs, scf_env, matrix_s, output_unit, marked_states,error) 00876 00877 TYPE(mo_set_p_type), DIMENSION(:), 00878 POINTER :: mos 00879 TYPE(particle_type), DIMENSION(:), 00880 POINTER :: particle_set 00881 TYPE(section_vals_type), POINTER :: dft_section 00882 TYPE(atomic_kind_type), DIMENSION(:), 00883 POINTER :: atomic_kind_set 00884 TYPE(cp_fm_p_type), DIMENSION(:), 00885 POINTER :: unoccupied_orbs 00886 TYPE(qs_scf_env_type), POINTER :: scf_env 00887 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00888 POINTER :: matrix_s 00889 INTEGER :: output_unit 00890 INTEGER, DIMENSION(:, :, :), POINTER :: marked_states 00891 TYPE(cp_error_type), INTENT(inout) :: error 00892 00893 CHARACTER(len=*), PARAMETER :: routineN = 'wfn_mix', 00894 routineP = moduleN//':'//routineN 00895 00896 INTEGER :: handle, i_rep, ispin, mark_ind, mark_number, n_rep, 00897 orig_mo_index, orig_spin_index, result_mo_index, result_spin_index 00898 LOGICAL :: explicit, failure, 00899 orig_is_virtual, overwrite_mos 00900 REAL(KIND=dp) :: orig_scale, orthonormality, 00901 result_scale 00902 TYPE(cp_fm_struct_type), POINTER :: fm_struct_vector 00903 TYPE(cp_fm_type), POINTER :: matrix_x, matrix_y 00904 TYPE(mo_set_p_type), DIMENSION(:), 00905 POINTER :: mos_new 00906 TYPE(section_vals_type), POINTER :: update_section, 00907 wfn_mix_section 00908 00909 failure=.FALSE. 00910 CALL timeset(routineN,handle) 00911 wfn_mix_section => section_vals_get_subs_vals(dft_section,"PRINT%WFN_MIX",error=error) 00912 CALL section_vals_get(wfn_mix_section,explicit=explicit, error=error) 00913 00914 ! only perform action if explicitly required 00915 IF (explicit) THEN 00916 00917 IF (output_unit>0) THEN 00918 WRITE(output_unit,'()') 00919 WRITE(output_unit,'(T2,A)') "Performing wfn mixing" 00920 WRITE(output_unit,'(T2,A)') "=====================" 00921 ENDIF 00922 00923 ALLOCATE(mos_new(SIZE(mos))) 00924 DO ispin=1,SIZE(mos) 00925 CALL duplicate_mo_set(mos_new(ispin)%mo_set,mos(ispin)%mo_set,error) 00926 ENDDO 00927 00928 ! a single vector matrix structure 00929 NULLIFY(fm_struct_vector) 00930 CALL cp_fm_struct_create(fm_struct_vector,template_fmstruct=mos(1)%mo_set%mo_coeff%matrix_struct, & 00931 ncol_global=1, error=error) 00932 CALL cp_fm_create(matrix_x,fm_struct_vector,name="x",error=error) 00933 CALL cp_fm_create(matrix_y,fm_struct_vector,name="y",error=error) 00934 CALL cp_fm_struct_release(fm_struct_vector,error) 00935 00936 update_section=>section_vals_get_subs_vals(wfn_mix_section,"UPDATE",error=error) 00937 CALL section_vals_get(update_section,n_repetition=n_rep,error=error) 00938 CALL section_vals_get(update_section,explicit=explicit, error=error) 00939 IF (.NOT. explicit) n_rep=0 00940 00941 DO i_rep=1,n_rep 00942 00943 CALL section_vals_val_get(update_section,"RESULT_MO_INDEX",i_rep_section=i_rep,i_val=result_mo_index,error=error) 00944 CALL section_vals_val_get(update_section,"RESULT_MARKED_STATE",i_rep_section=i_rep,i_val=mark_number,error=error) 00945 CALL section_vals_val_get(update_section,"RESULT_SPIN_INDEX",i_rep_section=i_rep,i_val=result_spin_index,error=error) 00946 CALL section_vals_val_get(update_section,"RESULT_SCALE",i_rep_section=i_rep,r_val=result_scale,error=error) 00947 00948 mark_ind=1 00949 IF(mark_number.GT.0)result_mo_index=marked_states(mark_number,result_spin_index,mark_ind) 00950 00951 CALL section_vals_val_get(update_section,"ORIG_MO_INDEX",i_rep_section=i_rep,i_val=orig_mo_index,error=error) 00952 CALL section_vals_val_get(update_section,"ORIG_MARKED_STATE",i_rep_section=i_rep,i_val=mark_number,error=error) 00953 CALL section_vals_val_get(update_section,"ORIG_SPIN_INDEX",i_rep_section=i_rep,i_val=orig_spin_index,error=error) 00954 CALL section_vals_val_get(update_section,"ORIG_SCALE",i_rep_section=i_rep,r_val=orig_scale,error=error) 00955 CALL section_vals_val_get(update_section,"ORIG_IS_VIRTUAL",i_rep_section=i_rep,l_val=orig_is_virtual,error=error) 00956 00957 IF(orig_is_virtual)mark_ind=2 00958 IF(mark_number.GT.0)orig_mo_index=marked_states(mark_number,orig_spin_index,mark_ind) 00959 00960 CALL section_vals_val_get(wfn_mix_section,"OVERWRITE_MOS",l_val=overwrite_mos,error=error) 00961 00962 ! first get a copy of the proper orig 00963 IF (.NOT. ORIG_IS_VIRTUAL) THEN 00964 CALL cp_fm_to_fm(mos(orig_spin_index)%mo_set%mo_coeff,matrix_x,1,mos(orig_spin_index)%mo_set%nmo-orig_mo_index+1,1) 00965 ELSE 00966 CALL cp_fm_to_fm(unoccupied_orbs(orig_spin_index)%matrix,matrix_x,1,orig_mo_index,1) 00967 ENDIF 00968 00969 ! now get a copy of the target 00970 CALL cp_fm_to_fm(mos_new(result_spin_index)%mo_set%mo_coeff,matrix_y, & 00971 1,mos_new(result_spin_index)%mo_set%nmo-result_mo_index+1,1) 00972 00973 ! properly mix 00974 CALL cp_fm_scale_and_add(result_scale,matrix_y,orig_scale,matrix_x,error) 00975 00976 ! and copy back in the result mos 00977 CALL cp_fm_to_fm(matrix_y,mos_new(result_spin_index)%mo_set%mo_coeff, & 00978 1,1,mos_new(result_spin_index)%mo_set%nmo-result_mo_index+1) 00979 00980 ENDDO 00981 00982 CALL cp_fm_release(matrix_x,error) 00983 CALL cp_fm_release(matrix_y,error) 00984 00985 IF(scf_env%method==special_diag_method_nr) THEN 00986 CALL calculate_orthonormality(orthonormality,mos,error=error) 00987 ELSE 00988 CALL calculate_orthonormality(orthonormality,mos,matrix_s(1)%matrix,error=error) 00989 END IF 00990 00991 IF (output_unit>0) THEN 00992 WRITE(output_unit,'()') 00993 WRITE(output_unit,'(T2,A,T61,E20.4)') & 00994 "Maximum deviation from MO S-orthonormality",orthonormality 00995 WRITE(output_unit,'(T2,A)') "Writing new MOs to file" 00996 ENDIF 00997 00998 ! *** Write WaveFunction restart file *** 00999 01000 DO ispin=1,SIZE(mos_new) 01001 IF(overwrite_mos)THEN 01002 CALL cp_fm_to_fm(mos_new(ispin)%mo_set%mo_coeff,mos(ispin)%mo_set%mo_coeff,error) 01003 IF(mos_new(1)%mo_set%use_mo_coeff_b)& 01004 CALL copy_fm_to_dbcsr(mos_new(ispin)%mo_set%mo_coeff,mos_new(ispin)%mo_set%mo_coeff_b,& 01005 error=error) 01006 END IF 01007 IF(mos(1)%mo_set%use_mo_coeff_b)& 01008 CALL copy_fm_to_dbcsr(mos_new(ispin)%mo_set%mo_coeff,mos(ispin)%mo_set%mo_coeff_b,& 01009 error=error) 01010 END DO 01011 CALL write_mo_set(mos_new,particle_set,dft_section=dft_section,scp=.FALSE., & 01012 atomic_kind_set=atomic_kind_set,error=error) 01013 01014 DO ispin=1,SIZE(mos_new) 01015 CALL deallocate_mo_set(mos_new(ispin)%mo_set,error) 01016 ENDDO 01017 DEALLOCATE(mos_new) 01018 01019 ENDIF 01020 01021 CALL timestop(handle) 01022 01023 END SUBROUTINE wfn_mix 01024 01025 ! ***************************************************************************** 01031 SUBROUTINE qs_scf_post_charges(input, logger, qs_env, rho, matrix_s, mos, error) 01032 TYPE(section_vals_type), POINTER :: input 01033 TYPE(cp_logger_type), POINTER :: logger 01034 TYPE(qs_environment_type), POINTER :: qs_env 01035 TYPE(qs_rho_type), POINTER :: rho 01036 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01037 POINTER :: matrix_s 01038 TYPE(mo_set_p_type), DIMENSION(:), 01039 POINTER :: mos 01040 TYPE(cp_error_type), INTENT(inout) :: error 01041 01042 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_charges', 01043 routineP = moduleN//':'//routineN 01044 01045 INTEGER :: handle, print_level, unit_nr 01046 LOGICAL :: print_it 01047 TYPE(particle_type), DIMENSION(:), 01048 POINTER :: particle_set 01049 TYPE(section_vals_type), POINTER :: density_fit_section, print_key 01050 01051 CALL timeset(routineN,handle) 01052 01053 NULLIFY(particle_set) 01054 CALL get_qs_env(qs_env=qs_env,& 01055 particle_set=particle_set,& 01056 error=error) 01057 01058 ! Compute the Mulliken charges 01059 print_key => section_vals_get_subs_vals(input,"DFT%PRINT%MULLIKEN", error=error) 01060 IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)) THEN 01061 unit_nr=cp_print_key_unit_nr(logger,input,"DFT%PRINT%MULLIKEN",extension=".mulliken",& 01062 middle_name="",log_filename=.FALSE.,error=error) 01063 print_level = 1 01064 CALL section_vals_val_get(print_key,"PRINT_GOP",l_val=print_it,error=error) 01065 IF (print_it) print_level = 2 01066 CALL section_vals_val_get(print_key,"PRINT_ALL",l_val=print_it,error=error) 01067 IF (print_it) print_level = 3 01068 CALL mulliken_population_analysis(qs_env,unit_nr,print_level,error) 01069 CALL cp_print_key_finished_output(unit_nr, logger,input,"DFT%PRINT%MULLIKEN",error=error) 01070 END IF 01071 01072 ! Compute the Lowdin charges 01073 print_key => section_vals_get_subs_vals(input,"DFT%PRINT%LOWDIN", error=error) 01074 IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)) THEN 01075 unit_nr = cp_print_key_unit_nr(logger,input,"DFT%PRINT%LOWDIN",extension=".lowdin",& 01076 log_filename=.FALSE.,error=error) 01077 print_level = 1 01078 CALL section_vals_val_get(print_key,"PRINT_GOP",l_val=print_it,error=error) 01079 IF (print_it) print_level = 2 01080 CALL section_vals_val_get(print_key,"PRINT_ALL",l_val=print_it,error=error) 01081 IF (print_it) print_level = 3 01082 CALL lowdin_population_analysis(qs_env,unit_nr,print_level,error) 01083 CALL cp_print_key_finished_output(unit_nr, logger,input,"DFT%PRINT%LOWDIN", error=error) 01084 END IF 01085 01086 ! Compute the RESP charges 01087 CALL resp_fit(qs_env,error) 01088 01089 ! Compute the Density Derived Atomic Point charges with the Bloechl scheme 01090 print_key => section_vals_get_subs_vals(input,"PROPERTIES%FIT_CHARGE", error=error) 01091 IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)) THEN 01092 unit_nr=cp_print_key_unit_nr(logger,input,"PROPERTIES%FIT_CHARGE",extension=".Fitcharge",& 01093 log_filename=.FALSE.,error=error) 01094 density_fit_section => section_vals_get_subs_vals(input,"DFT%DENSITY_FITTING", error=error) 01095 CALL get_ddapc(qs_env,.FALSE.,density_fit_section,iwc=unit_nr,error=error) 01096 CALL cp_print_key_finished_output(unit_nr, logger,input,"PROPERTIES%FIT_CHARGE", error=error) 01097 END IF 01098 01099 01100 CALL timestop(handle) 01101 01102 END SUBROUTINE qs_scf_post_charges 01103 01104 ! ***************************************************************************** 01110 SUBROUTINE qs_scf_post_loc_dip(input, dft_control, qs_loc_env, logger, qs_env, error) 01111 TYPE(section_vals_type), POINTER :: input 01112 TYPE(dft_control_type), POINTER :: dft_control 01113 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 01114 TYPE(cp_logger_type), POINTER :: logger 01115 TYPE(qs_environment_type), POINTER :: qs_env 01116 TYPE(cp_error_type), INTENT(inout) :: error 01117 01118 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_loc_dip', 01119 routineP = moduleN//':'//routineN 01120 01121 CHARACTER(LEN=default_string_length) :: description, 01122 descriptionThisDip, iter 01123 COMPLEX(KIND=dp) :: zeta 01124 COMPLEX(KIND=dp), DIMENSION(3) :: ggamma, zphase 01125 INTEGER :: handle, i, ispins, j, n_rep, 01126 reference, unit_nr 01127 LOGICAL :: do_berry, failure, 01128 first_time, ghost 01129 REAL(KIND=dp) :: charge_tot, theta, zeff, zwfc 01130 REAL(KIND=dp), DIMENSION(3) :: ci, dipole, dipole_old, gvec, 01131 rcc, ria 01132 REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point 01133 TYPE(cell_type), POINTER :: cell 01134 TYPE(particle_type), DIMENSION(:), 01135 POINTER :: particle_set 01136 TYPE(section_vals_type), POINTER :: print_key 01137 01138 CALL timeset(routineN,handle) 01139 01140 failure = .FALSE. 01141 print_key => section_vals_get_subs_vals(input,"DFT%LOCALIZE%PRINT%TOTAL_DIPOLE", error=error) 01142 IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,first_time=first_time,& 01143 error=error),cp_p_file)) THEN 01144 NULLIFY(cell, particle_set, ref_point) 01145 CALL get_qs_env(qs_env=qs_env,& 01146 cell=cell,& 01147 particle_set=particle_set,& 01148 error=error) 01149 01150 reference = section_get_ival(print_key,keyword_name="REFERENCE",error=error) 01151 CALL section_vals_val_get(print_key,"REF_POINT",r_vals=ref_point,error=error) 01152 CALL section_vals_val_get(print_key,"PERIODIC",l_val=do_berry,error=error) 01153 description='[DIPOLE]' 01154 descriptionThisDip='[TOTAL_DIPOLE]' 01155 CALL get_reference_point(rcc,qs_env=qs_env,reference=reference,ref_point=ref_point,error=error) 01156 01157 dipole=0.0_dp 01158 IF (do_berry) THEN 01159 rcc = pbc(rcc,cell) 01160 charge_tot = REAL(dft_control%charge,KIND=dp) 01161 ria = twopi * MATMUL(cell%h_inv, rcc) 01162 zphase = CMPLX(COS(ria), SIN(ria), KIND=dp)**charge_tot 01163 ggamma = CMPLX(1.0_dp,0.0_dp,KIND=dp) 01164 01165 ! Nuclear charges 01166 DO i=1,SIZE(particle_set) 01167 CALL get_atomic_kind(particle_set(i)%atomic_kind,ghost=ghost) 01168 IF (.NOT.ghost) THEN 01169 CALL get_atomic_kind(particle_set(i)%atomic_kind,core_charge=zeff) 01170 ria = pbc(particle_set(i)%r, cell) 01171 DO j = 1, 3 01172 gvec = twopi*cell%h_inv(j,:) 01173 theta = SUM(ria(:)*gvec(:)) 01174 zeta = CMPLX(COS(theta),SIN(theta),KIND=dp)**(zeff) 01175 ggamma(j) = ggamma(j) * zeta 01176 END DO 01177 END IF 01178 END DO 01179 01180 ! Charges of the wfc involved 01181 ! Warning, this assumes the same occupation for all states 01182 zwfc = 3.0_dp-REAL(dft_control%nspins,dp) 01183 01184 DO ispins=1,dft_control%nspins 01185 DO i=1,SIZE(qs_loc_env%localized_wfn_control%centers_set(ispins)%array,2) 01186 ria = pbc(qs_loc_env%localized_wfn_control%centers_set(ispins)%array(1:3,i),cell) 01187 DO j = 1, 3 01188 gvec = twopi*cell%h_inv(j,:) 01189 theta = SUM(ria(:)*gvec(:)) 01190 zeta = CMPLX(COS(theta),SIN(theta),KIND=dp)**(-zwfc) 01191 ggamma(j) = ggamma(j) * zeta 01192 END DO 01193 ENDDO 01194 ENDDO 01195 ggamma = ggamma * zphase 01196 ci = AIMAG(LOG(ggamma))/twopi 01197 dipole = MATMUL ( cell%hmat, ci ) 01198 ELSE 01199 ! Charges of the atoms involved 01200 DO i=1,SIZE(particle_set) 01201 CALL get_atomic_kind(particle_set(i)%atomic_kind,ghost=ghost) 01202 IF (.NOT.ghost) THEN 01203 CALL get_atomic_kind(particle_set(i)%atomic_kind,core_charge=zeff) 01204 ria = pbc(particle_set(i)%r, cell) 01205 dipole=dipole + zeff*(ria-rcc) 01206 END IF 01207 END DO 01208 01209 ! Charges of the wfc involved 01210 ! Warning, this assumes the same occupation for all states 01211 zwfc = 3.0_dp-REAL(dft_control%nspins,dp) 01212 01213 DO ispins=1,dft_control%nspins 01214 DO i=1,SIZE(qs_loc_env%localized_wfn_control%centers_set(ispins)%array,2) 01215 ria = pbc(qs_loc_env%localized_wfn_control%centers_set(ispins)%array(1:3,i),cell) 01216 dipole=dipole - zwfc * (ria-rcc) 01217 ENDDO 01218 ENDDO 01219 END IF 01220 01221 ! Print and possibly store results 01222 unit_nr=cp_print_key_unit_nr(logger,print_key,extension=".Dipole",& 01223 middle_name="TOTAL_DIPOLE",error=error) 01224 IF (unit_nr>0) THEN 01225 IF (first_time) THEN 01226 WRITE(unit=unit_nr,fmt="(A,T31,A,T88,A,T136,A)")& 01227 "# iter_level","dipole(x,y,z)[atomic units]",& 01228 "dipole(x,y,z)[debye]",& 01229 "delta_dipole(x,y,z)[atomic units]" 01230 END IF 01231 iter=cp_iter_string(logger%iter_info,error=error) 01232 CALL get_results(qs_env%results,descriptionThisDip,n_rep=n_rep,error=error) 01233 IF(n_rep==0)THEN 01234 dipole_old=0._dp 01235 ELSE 01236 CALL get_results(qs_env%results,descriptionThisDip,dipole_old,nval=n_rep,error=error) 01237 END IF 01238 IF (do_berry) THEN 01239 WRITE(unit=unit_nr,fmt="(a,9(es18.8))")& 01240 iter(1:15), dipole, dipole*debye, pbc(dipole-dipole_old,cell) 01241 ELSE 01242 WRITE(unit=unit_nr,fmt="(a,9(es18.8))")& 01243 iter(1:15), dipole, dipole*debye, (dipole-dipole_old) 01244 END IF 01245 END IF 01246 CALL cp_print_key_finished_output(unit_nr,logger,print_key,error=error) 01247 CALL cp_results_erase(qs_env%results,description,error=error) 01248 CALL put_results(qs_env%results,description,dipole,error) 01249 CALL cp_results_erase(qs_env%results,descriptionThisDip,error=error) 01250 CALL put_results(qs_env%results,descriptionThisDip,dipole,error) 01251 END IF 01252 01253 CALL timestop(handle) 01254 01255 END SUBROUTINE qs_scf_post_loc_dip 01256 01257 ! ***************************************************************************** 01263 SUBROUTINE qs_scf_post_occ_cubes(input, dft_section, dft_control, logger, qs_env,& 01264 mo_coeff, wf_g, wf_r, particles, homo, ispin, basis_set_id, error) 01265 TYPE(section_vals_type), POINTER :: input, dft_section 01266 TYPE(dft_control_type), POINTER :: dft_control 01267 TYPE(cp_logger_type), POINTER :: logger 01268 TYPE(qs_environment_type), POINTER :: qs_env 01269 TYPE(cp_fm_type), POINTER :: mo_coeff 01270 TYPE(pw_p_type) :: wf_g, wf_r 01271 TYPE(particle_list_type), POINTER :: particles 01272 INTEGER, INTENT(IN) :: homo, ispin 01273 INTEGER, INTENT(IN), OPTIONAL :: basis_set_id 01274 TYPE(cp_error_type), INTENT(inout) :: error 01275 01276 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_occ_cubes', 01277 routineP = moduleN//':'//routineN 01278 01279 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title 01280 INTEGER :: handle, ivector, 01281 my_basis_set_id, nhomo, 01282 unit_nr 01283 LOGICAL :: append_cube 01284 TYPE(atomic_kind_type), DIMENSION(:), 01285 POINTER :: atomic_kind_set 01286 TYPE(cell_type), POINTER :: cell 01287 TYPE(particle_type), DIMENSION(:), 01288 POINTER :: particle_set 01289 TYPE(pw_env_type), POINTER :: pw_env 01290 01291 CALL timeset(routineN,handle) 01292 01293 IF(PRESENT(basis_set_id)) THEN 01294 my_basis_set_id = basis_set_id 01295 ELSE 01296 my_basis_set_id = use_orb_basis_set 01297 END IF 01298 01299 IF (BTEST(cp_print_key_should_output(logger%iter_info,dft_section,"PRINT%MO_CUBES",& 01300 error=error),cp_p_file) .AND. section_get_lval(dft_section,"PRINT%MO_CUBES%WRITE_CUBE",error=error)) THEN 01301 nhomo=section_get_ival(dft_section,"PRINT%MO_CUBES%NHOMO",error=error) 01302 append_cube = section_get_lval(dft_section,"PRINT%MO_CUBES%APPEND",error=error) 01303 my_pos_cube="REWIND" 01304 IF(append_cube) THEN 01305 my_pos_cube="APPEND" 01306 END IF 01307 IF (nhomo==-1) nhomo=homo 01308 DO ivector=MAX(1,homo-nhomo+1),homo 01309 CALL get_qs_env(qs_env=qs_env,& 01310 atomic_kind_set=atomic_kind_set,& 01311 cell=cell,& 01312 particle_set=particle_set,& 01313 pw_env=pw_env,& 01314 error=error) 01315 CALL calculate_wavefunction(mo_coeff,ivector,wf_r,wf_g,atomic_kind_set,& 01316 cell,dft_control,particle_set,pw_env ,basis_set_id=my_basis_set_id, error=error) 01317 WRITE(filename,'(a4,I5.5,a1,I1.1)')"WFN_",ivector,"_",ispin 01318 unit_nr=cp_print_key_unit_nr(logger,input,"DFT%PRINT%MO_CUBES",extension=".cube",& 01319 middle_name=TRIM(filename),file_position=my_pos_cube,log_filename=.FALSE.,error=error) 01320 WRITE(title,*) "WAVEFUNCTION ",ivector," spin ",ispin," i.e. HOMO - ",ivector-homo 01321 CALL pw_to_cube(wf_r%pw,unit_nr,title,particles=particles,& 01322 stride=section_get_ivals(dft_section,"PRINT%MO_CUBES%STRIDE",error=error),& 01323 error=error) 01324 CALL cp_print_key_finished_output(unit_nr,logger,input,"DFT%PRINT%MO_CUBES",error=error) 01325 ENDDO 01326 END IF 01327 01328 CALL timestop(handle) 01329 01330 END SUBROUTINE qs_scf_post_occ_cubes 01331 01332 ! ***************************************************************************** 01338 SUBROUTINE qs_scf_post_unocc_cubes(input, dft_section, dft_control, logger, qs_env,& 01339 unoccupied_orbs, wf_g, wf_r, particles, nlumos, homo, ispin, error) 01340 TYPE(section_vals_type), POINTER :: input, dft_section 01341 TYPE(dft_control_type), POINTER :: dft_control 01342 TYPE(cp_logger_type), POINTER :: logger 01343 TYPE(qs_environment_type), POINTER :: qs_env 01344 TYPE(cp_fm_p_type), DIMENSION(:), 01345 POINTER :: unoccupied_orbs 01346 TYPE(pw_p_type) :: wf_g, wf_r 01347 TYPE(particle_list_type), POINTER :: particles 01348 INTEGER, INTENT(IN) :: nlumos, homo, ispin 01349 TYPE(cp_error_type), INTENT(inout) :: error 01350 01351 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_unocc_cubes', 01352 routineP = moduleN//':'//routineN 01353 01354 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title 01355 INTEGER :: handle, ivector, unit_nr 01356 LOGICAL :: append_cube 01357 TYPE(atomic_kind_type), DIMENSION(:), 01358 POINTER :: atomic_kind_set 01359 TYPE(cell_type), POINTER :: cell 01360 TYPE(particle_type), DIMENSION(:), 01361 POINTER :: particle_set 01362 TYPE(pw_env_type), POINTER :: pw_env 01363 01364 CALL timeset(routineN,handle) 01365 01366 IF (BTEST(cp_print_key_should_output(logger%iter_info,dft_section,"PRINT%MO_CUBES",error=error),cp_p_file)& 01367 .AND. section_get_lval(dft_section,"PRINT%MO_CUBES%WRITE_CUBE",error=error) ) THEN 01368 NULLIFY(atomic_kind_set, particle_set, pw_env, cell) 01369 append_cube = section_get_lval(dft_section,"PRINT%MO_CUBES%APPEND",error=error) 01370 my_pos_cube="REWIND" 01371 IF(append_cube) THEN 01372 my_pos_cube="APPEND" 01373 END IF 01374 DO ivector=1,nlumos 01375 CALL get_qs_env(qs_env=qs_env,& 01376 atomic_kind_set=atomic_kind_set,& 01377 cell=cell,& 01378 particle_set=particle_set,& 01379 pw_env=pw_env,& 01380 error=error) 01381 CALL calculate_wavefunction(unoccupied_orbs(ispin)%matrix, ivector, wf_r, wf_g, atomic_kind_set,& 01382 cell, dft_control, particle_set, pw_env, error=error) 01383 WRITE(filename,'(a4,I5.5,a1,I1.1)')"WFN_",homo+ivector,"_",ispin 01384 unit_nr=cp_print_key_unit_nr(logger,input,"DFT%PRINT%MO_CUBES",extension=".cube",& 01385 middle_name=TRIM(filename),file_position=my_pos_cube,log_filename=.FALSE.,error=error) 01386 WRITE(title,*) "WAVEFUNCTION ",homo+ivector," spin ",ispin," i.e. LUMO + ",ivector-1 01387 CALL pw_to_cube(wf_r%pw, unit_nr, title, particles=particles,& 01388 stride=section_get_ivals(dft_section,"PRINT%MO_CUBES%STRIDE",error=error),& 01389 error=error) 01390 CALL cp_print_key_finished_output(unit_nr,logger,input,"DFT%PRINT%MO_CUBES",error=error) 01391 ENDDO 01392 ENDIF 01393 01394 CALL timestop(handle) 01395 01396 END SUBROUTINE qs_scf_post_unocc_cubes 01397 01398 ! ***************************************************************************** 01404 SUBROUTINE qs_scf_post_optc(input, dft_section, dft_control, logger, qs_env,& 01405 unoccupied_orbs, unoccupied_evals, output_unit, error) 01406 TYPE(section_vals_type), POINTER :: input, dft_section 01407 TYPE(dft_control_type), POINTER :: dft_control 01408 TYPE(cp_logger_type), POINTER :: logger 01409 TYPE(qs_environment_type), POINTER :: qs_env 01410 TYPE(cp_fm_p_type), DIMENSION(:), 01411 POINTER :: unoccupied_orbs 01412 TYPE(cp_1d_r_p_type), DIMENSION(:), 01413 POINTER :: unoccupied_evals 01414 INTEGER, INTENT(IN) :: output_unit 01415 TYPE(cp_error_type), INTENT(inout) :: error 01416 01417 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_optc', 01418 routineP = moduleN//':'//routineN 01419 01420 CHARACTER(LEN=default_path_length) :: filename 01421 INTEGER :: handle, homo, ispin, nmo, 01422 unit_nr 01423 LOGICAL :: homoEQnmo 01424 REAL(KIND=dp) :: volume 01425 REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues 01426 TYPE(cell_type), POINTER :: cell 01427 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01428 POINTER :: matrix_s 01429 TYPE(cp_fm_type), POINTER :: mo_coeff 01430 TYPE(mo_set_p_type), DIMENSION(:), 01431 POINTER :: mos 01432 01433 CALL timeset(routineN,handle) 01434 01435 IF (BTEST(cp_print_key_should_output(logger%iter_info,input,& 01436 "DFT%PRINT%OPTICAL_CONDUCTIVITY",error=error),cp_p_file)) THEN 01437 NULLIFY(mos, cell, mo_coeff, mo_eigenvalues) 01438 CALL get_qs_env(qs_env, mos=mos, cell=cell, matrix_s=matrix_s,error=error) 01439 IF (output_unit>0) WRITE(output_unit,*) " " 01440 IF (output_unit>0) WRITE(output_unit,*) " Computing the optical conductivity " 01441 IF (output_unit>0) WRITE(output_unit,*) " Experimental version " 01442 IF (output_unit>0) WRITE(output_unit,*) " Check the code before believing results " 01443 homoEQnmo=.TRUE. 01444 DO ispin=1,dft_control%nspins 01445 CALL get_mo_set(mo_set=mos(ispin)%mo_set,homo=homo,nmo=nmo) 01446 IF (homo.NE.nmo) homoEQnmo=.FALSE. 01447 ENDDO 01448 IF (.NOT. homoEQnmo) THEN 01449 IF (output_unit>0) WRITE(output_unit,*) " homo.NE.nmo : skip optical conductivty " 01450 ELSE 01451 IF (.NOT.( BTEST(cp_print_key_should_output(logger%iter_info,dft_section,"PRINT%MO_CUBES",& 01452 error=error),cp_p_file) )) & 01453 CALL stop_program(routineN,moduleN,__LINE__,"Needs MO_CUBES to be activated") 01454 filename="CONDUCTIVITY" 01455 unit_nr=cp_print_key_unit_nr(logger,input,"DFT%PRINT%OPTICAL_CONDUCTIVITY",& 01456 extension=".data",middle_name=TRIM(filename),log_filename=.FALSE.,error=error) 01457 CALL get_cell(cell,deth=volume) 01458 01459 DO ispin=1,dft_control%nspins 01460 CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff, & 01461 eigenvalues=mo_eigenvalues) 01462 CALL optical_conductivity(matrix_s, mo_coeff, mo_eigenvalues, unoccupied_orbs(ispin)%matrix,& 01463 unoccupied_evals(ispin)%array, volume, unit_nr, error=error) 01464 ENDDO 01465 CALL cp_print_key_finished_output(unit_nr,logger,input,"DFT%PRINT%OPTICAL_CONDUCTIVITY",& 01466 error=error) 01467 IF (output_unit>0) WRITE(output_unit,*) " Results written to file ",filename 01468 ENDIF 01469 ENDIF 01470 01471 CALL timestop(handle) 01472 01473 END SUBROUTINE qs_scf_post_optc 01474 01475 ! ***************************************************************************** 01481 SUBROUTINE qs_scf_post_moments(input, logger, qs_env, output_unit, error) 01482 TYPE(section_vals_type), POINTER :: input 01483 TYPE(cp_logger_type), POINTER :: logger 01484 TYPE(qs_environment_type), POINTER :: qs_env 01485 INTEGER, INTENT(IN) :: output_unit 01486 TYPE(cp_error_type), INTENT(inout) :: error 01487 01488 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_moments', 01489 routineP = moduleN//':'//routineN 01490 01491 CHARACTER(LEN=default_path_length) :: filename 01492 INTEGER :: handle, maxmom, reference, 01493 unit_nr 01494 LOGICAL :: magnetic, periodic 01495 REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point 01496 TYPE(section_vals_type), POINTER :: print_key 01497 01498 CALL timeset(routineN,handle) 01499 01500 print_key => section_vals_get_subs_vals(section_vals=input,& 01501 subsection_name="DFT%PRINT%MOMENTS",error=error) 01502 01503 IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)) THEN 01504 01505 maxmom = section_get_ival(section_vals=input,& 01506 keyword_name="DFT%PRINT%MOMENTS%MAX_MOMENT",error=error) 01507 periodic = section_get_lval(section_vals=input,& 01508 keyword_name="DFT%PRINT%MOMENTS%PERIODIC",error=error) 01509 reference = section_get_ival(section_vals=input,& 01510 keyword_name="DFT%PRINT%MOMENTS%REFERENCE",error=error) 01511 magnetic = section_get_lval(section_vals=input,& 01512 keyword_name="DFT%PRINT%MOMENTS%MAGNETIC",error=error) 01513 NULLIFY ( ref_point ) 01514 CALL section_vals_val_get(input,"DFT%PRINT%MOMENTS%REF_POINT",r_vals=ref_point,error=error) 01515 unit_nr = cp_print_key_unit_nr(logger=logger,basis_section=input,& 01516 print_key_path="DFT%PRINT%MOMENTS",extension=".dat",& 01517 middle_name="moments",log_filename=.FALSE.,error=error) 01518 01519 IF (output_unit>0) THEN 01520 IF(unit_nr /= output_unit) THEN 01521 INQUIRE (UNIT=unit_nr,NAME=filename) 01522 WRITE (UNIT=output_unit,FMT="(/,T2,A,2(/,T3,A),/)")& 01523 "MOMENTS","The electric/magnetic moments are written to file:",& 01524 TRIM(filename) 01525 ELSE 01526 WRITE (UNIT=output_unit,FMT="(/,T2,A)") "ELECTRIC/MAGNETIC MOMENTS" 01527 END IF 01528 END IF 01529 01530 IF (periodic) THEN 01531 CALL qs_moment_berry_phase(qs_env,magnetic,maxmom,reference,ref_point,unit_nr,error) 01532 ELSE 01533 CALL qs_moment_locop(qs_env,magnetic,maxmom,reference,ref_point,unit_nr,error) 01534 END IF 01535 01536 CALL cp_print_key_finished_output(unit_nr=unit_nr,logger=logger,& 01537 basis_section=input,print_key_path="DFT%PRINT%MOMENTS",& 01538 error=error) 01539 END IF 01540 01541 CALL timestop(handle) 01542 01543 END SUBROUTINE qs_scf_post_moments 01544 01545 ! ***************************************************************************** 01551 SUBROUTINE qs_scf_post_xray(input,dft_section,logger,qs_env,output_unit,error) 01552 01553 TYPE(section_vals_type), POINTER :: input, dft_section 01554 TYPE(cp_logger_type), POINTER :: logger 01555 TYPE(qs_environment_type), POINTER :: qs_env 01556 INTEGER, INTENT(IN) :: output_unit 01557 TYPE(cp_error_type), INTENT(inout) :: error 01558 01559 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_post_xray', 01560 routineP = moduleN//':'//routineN 01561 01562 CHARACTER(LEN=default_path_length) :: filename 01563 INTEGER :: handle, unit_nr 01564 REAL(KIND=dp) :: q_max 01565 TYPE(section_vals_type), POINTER :: print_key 01566 01567 CALL timeset(routineN,handle) 01568 01569 print_key => section_vals_get_subs_vals(section_vals=input,& 01570 subsection_name="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM",& 01571 error=error) 01572 01573 IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),cp_p_file)) THEN 01574 q_max = section_get_rval(section_vals=dft_section,& 01575 keyword_name="PRINT%XRAY_DIFFRACTION_SPECTRUM%Q_MAX",& 01576 error=error) 01577 unit_nr = cp_print_key_unit_nr(logger=logger,& 01578 basis_section=input,& 01579 print_key_path="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM",& 01580 extension=".dat",& 01581 middle_name="xrd",& 01582 log_filename=.FALSE.,& 01583 error=error) 01584 IF (output_unit>0) THEN 01585 INQUIRE (UNIT=unit_nr,NAME=filename) 01586 WRITE (UNIT=output_unit,FMT="(/,/,T2,A)")& 01587 "X-RAY DIFFRACTION SPECTRUM" 01588 IF (unit_nr /= output_unit) THEN 01589 WRITE (UNIT=output_unit,FMT="(/,T3,A,/,/,T3,A,/)")& 01590 "The coherent X-ray diffraction spectrum is written to the file:",& 01591 TRIM(filename) 01592 END IF 01593 END IF 01594 CALL xray_diffraction_spectrum(qs_env=qs_env,& 01595 unit_number=unit_nr,& 01596 q_max=q_max,& 01597 error=error) 01598 CALL cp_print_key_finished_output(unit_nr=unit_nr,& 01599 logger=logger,& 01600 basis_section=input,& 01601 print_key_path="DFT%PRINT%XRAY_DIFFRACTION_SPECTRUM",& 01602 error=error) 01603 END IF 01604 01605 CALL timestop(handle) 01606 01607 END SUBROUTINE qs_scf_post_xray 01608 01609 ! ***************************************************************************** 01615 SUBROUTINE qs_scf_post_efg(input, logger, qs_env, error) 01616 TYPE(section_vals_type), POINTER :: input 01617 TYPE(cp_logger_type), POINTER :: logger 01618 TYPE(qs_environment_type), POINTER :: qs_env 01619 TYPE(cp_error_type), INTENT(inout) :: error 01620 01621 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_efg', 01622 routineP = moduleN//':'//routineN 01623 01624 INTEGER :: handle 01625 TYPE(section_vals_type), POINTER :: print_key 01626 01627 CALL timeset(routineN,handle) 01628 01629 print_key => section_vals_get_subs_vals(section_vals=input,& 01630 subsection_name="DFT%PRINT%ELECTRIC_FIELD_GRADIENT",& 01631 error=error) 01632 IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),& 01633 cp_p_file)) THEN 01634 CALL qs_efg_calc(qs_env=qs_env,error=error) 01635 END IF 01636 01637 CALL timestop(handle) 01638 01639 END SUBROUTINE qs_scf_post_efg 01640 01641 ! ***************************************************************************** 01647 SUBROUTINE qs_scf_post_et(input, qs_env, dft_control, error) 01648 TYPE(section_vals_type), POINTER :: input 01649 TYPE(qs_environment_type), POINTER :: qs_env 01650 TYPE(dft_control_type), POINTER :: dft_control 01651 TYPE(cp_error_type), INTENT(inout) :: error 01652 01653 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_et', 01654 routineP = moduleN//':'//routineN 01655 01656 INTEGER :: handle, ispin, stat 01657 LOGICAL :: do_et, failure 01658 TYPE(cp_fm_p_type), DIMENSION(:), 01659 POINTER :: my_mos 01660 TYPE(section_vals_type), POINTER :: et_section 01661 01662 CALL timeset(routineN,handle) 01663 failure=.FALSE. 01664 01665 do_et=.FALSE. 01666 et_section => section_vals_get_subs_vals(input,"PROPERTIES%ET_COUPLING",& 01667 error=error) 01668 CALL section_vals_get(et_section,explicit=do_et,error=error) 01669 IF(do_et)THEN 01670 IF(qs_env%et_coupling%first_run)THEN 01671 NULLIFY(my_mos) 01672 ALLOCATE(my_mos(dft_control%nspins),STAT=stat) 01673 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01674 ALLOCATE(qs_env%et_coupling%et_mo_coeff(dft_control%nspins),STAT=stat) 01675 CPPrecondition(stat==0,cp_failure_level,routineP,error,failure) 01676 DO ispin =1,dft_control%nspins 01677 NULLIFY(my_mos(ispin)%matrix) 01678 CALL cp_fm_create(matrix=my_mos(ispin)%matrix,& 01679 matrix_struct=qs_env%mos(ispin)%mo_set%mo_coeff%matrix_struct,& 01680 name="FIRST_RUN_COEFF"//TRIM(ADJUSTL(cp_to_string(ispin)))//"MATRIX",& 01681 error=error) 01682 CALL cp_fm_to_fm(qs_env%mos(ispin)%mo_set%mo_coeff,& 01683 my_mos(ispin)%matrix,error=error) 01684 END DO 01685 CALL set_et_coupling_type(qs_env%et_coupling,et_mo_coeff=my_mos,error=error) 01686 DEALLOCATE(my_mos) 01687 END IF 01688 END IF 01689 01690 CALL timestop(handle) 01691 01692 END SUBROUTINE qs_scf_post_et 01693 01694 01695 01696 ! ***************************************************************************** 01702 01703 SUBROUTINE qs_scf_post_elf(input, logger, qs_env, error) 01704 TYPE(section_vals_type), POINTER :: input 01705 TYPE(cp_logger_type), POINTER :: logger 01706 TYPE(qs_environment_type), POINTER :: qs_env 01707 TYPE(cp_error_type), INTENT(inout) :: error 01708 01709 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_elf', 01710 routineP = moduleN//':'//routineN 01711 01712 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube, title 01713 INTEGER :: handle, ispin, output_unit, 01714 stat, unit_nr 01715 LOGICAL :: append_cube, do_elf, failure, 01716 gapw 01717 REAL(dp) :: rho_cutoff 01718 TYPE(cp_subsys_type), POINTER :: subsys 01719 TYPE(dft_control_type), POINTER :: dft_control 01720 TYPE(particle_list_type), POINTER :: particles 01721 TYPE(pw_env_type), POINTER :: pw_env 01722 TYPE(pw_p_type), DIMENSION(:), POINTER :: elf_r 01723 TYPE(pw_pool_p_type), DIMENSION(:), 01724 POINTER :: pw_pools 01725 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 01726 TYPE(section_vals_type), POINTER :: elf_section 01727 01728 CALL timeset(routineN,handle) 01729 failure=.FALSE. 01730 output_unit= cp_logger_get_default_io_unit(logger) 01731 01732 do_elf=.FALSE. 01733 elf_section => section_vals_get_subs_vals(input,"DFT%PRINT%ELF_CUBE",error=error) 01734 CALL section_vals_get(elf_section,explicit=do_elf,error=error) 01735 IF(do_elf)THEN 01736 NULLIFY(dft_control, pw_env, auxbas_pw_pool, pw_pools, particles, subsys, elf_r) 01737 CALL get_qs_env(qs_env,dft_control=dft_control, pw_env=pw_env, subsys=subsys, error=error) 01738 CALL cp_subsys_get(subsys,particles=particles,error=error) 01739 01740 gapw=dft_control%qs_control%gapw 01741 IF(.NOT. gapw) THEN 01742 ! allocate 01743 ALLOCATE(elf_r(dft_control%nspins),stat=stat) 01744 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01745 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,& 01746 pw_pools=pw_pools,error=error) 01747 DO ispin = 1,dft_control%nspins 01748 CALL pw_pool_create_pw(auxbas_pw_pool,elf_r(ispin)%pw,& 01749 use_data = REALDATA3D,& 01750 in_space = REALSPACE, error=error) 01751 CALL pw_zero(elf_r(ispin)%pw, error=error) 01752 END DO 01753 01754 IF (output_unit>0) THEN 01755 WRITE (UNIT=output_unit,FMT="(/,T15,A,/)")& 01756 " ----- ELF is computed on the real space grid -----" 01757 END IF 01758 rho_cutoff=section_get_rval(elf_section,"density_cutoff",error=error) 01759 CALL qs_elf_calc(qs_env, elf_r, rho_cutoff, error=error) 01760 01761 ! write ELF into cube file 01762 append_cube = section_get_lval(elf_section,"APPEND",error=error) 01763 my_pos_cube="REWIND" 01764 IF(append_cube) THEN 01765 my_pos_cube="APPEND" 01766 END IF 01767 01768 DO ispin = 1,dft_control%nspins 01769 WRITE(filename,'(a5,I1.1)')"ELF_S",ispin 01770 WRITE(title,*) "ELF spin ", ispin 01771 unit_nr=cp_print_key_unit_nr(logger,input,"DFT%PRINT%ELF_CUBE",extension=".cube",& 01772 middle_name=TRIM(filename),file_position=my_pos_cube,log_filename=.FALSE.,error=error) 01773 IF (output_unit>0) THEN 01774 INQUIRE (UNIT=unit_nr,NAME=filename) 01775 WRITE (UNIT=output_unit,FMT="(/,T2,A,/,/,T2,A)")& 01776 "ELF is written in cube file format to the file:",& 01777 TRIM(filename) 01778 END IF 01779 01780 CALL pw_to_cube(elf_r(ispin)%pw,unit_nr,title,particles=particles,& 01781 stride=section_get_ivals(elf_section,"STRIDE",error=error),& 01782 error=error) 01783 CALL cp_print_key_finished_output(unit_nr,logger,input,"DFT%PRINT%ELF_CUBE",error=error) 01784 01785 CALL pw_pool_give_back_pw(auxbas_pw_pool,elf_r(ispin)%pw, error=error) 01786 END DO 01787 01788 ! deallocate 01789 DEALLOCATE(elf_r,stat=stat) 01790 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01791 01792 ELSE 01793 ! not implemented 01794 CALL cp_unimplemented_error(fromWhere=routineP, & 01795 message="ELF not implemented for GAPW calculations!!", & 01796 error=error, error_level=cp_warning_level) 01797 01798 END IF 01799 01800 END IF 01801 01802 CALL timestop(handle) 01803 01804 END SUBROUTINE qs_scf_post_elf 01805 01806 01807 SUBROUTINE qs_elf_calc(qs_env, elf_r, rho_cutoff, error) 01808 01809 TYPE(qs_environment_type), POINTER :: qs_env 01810 TYPE(pw_p_type), DIMENSION(:), POINTER :: elf_r 01811 REAL(kind=dp), INTENT(IN) :: rho_cutoff 01812 TYPE(cp_error_type), INTENT(inout) :: error 01813 01814 CHARACTER(len=*), PARAMETER :: routineN = 'qs_elf_calc', 01815 routineP = moduleN//':'//routineN 01816 REAL(KIND=dp), PARAMETER :: ELFCUT = 0.0001_dp, ELFEPS = 2.87E-05_dp, 01817 f18 = (1.0_dp/8.0_dp), f23 = (2.0_dp/3.0_dp), f53 = (5.0_dp/3.0_dp) 01818 01819 INTEGER :: handle, i, idir, ispin, j, k, 01820 nspin, stat 01821 INTEGER, DIMENSION(2, 3) :: bo 01822 INTEGER, DIMENSION(3, 3) :: nd 01823 LOGICAL :: deriv_pw, failure 01824 REAL(kind=dp) :: cfermi, dum, elf_kernel, 01825 norm_drho, rho_53, udvol 01826 TYPE(pw_env_type), POINTER :: pw_env 01827 TYPE(pw_p_type), DIMENSION(:), POINTER :: drho_g, drho_r, rho_r, tau_g, 01828 tau_r 01829 TYPE(pw_pool_p_type), DIMENSION(:), 01830 POINTER :: pw_pools 01831 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 01832 TYPE(pw_type), POINTER :: tmp_g 01833 TYPE(qs_rho_type), POINTER :: rho 01834 01835 ! for spin restricted systems 01836 01837 CALL timeset(routineN,handle) 01838 failure=.FALSE. 01839 01840 NULLIFY(rho,rho_r, drho_g, drho_r, tau_r, tau_g, pw_env, auxbas_pw_pool, pw_pools,tmp_g ) 01841 01842 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho, error=error) 01843 01844 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,& 01845 pw_pools=pw_pools,error=error) 01846 nspin = SIZE(rho%rho_r) 01847 bo = rho%rho_r(1)%pw%pw_grid%bounds_local 01848 cfermi = (3.0_dp/10.0_dp)*(pi*pi*3.0_dp)**f23 01849 01850 ALLOCATE(rho_r(nspin),STAT=stat) 01851 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01852 ALLOCATE(tau_r(nspin),STAT=stat) 01853 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01854 ALLOCATE(tau_g(nspin),STAT=stat) 01855 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01856 ALLOCATE(drho_r(3*nspin),STAT=stat) 01857 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01858 ALLOCATE(drho_g(3*nspin),STAT=stat) 01859 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01860 01861 DO ispin = 1,nspin 01862 rho_r(ispin)%pw => rho%rho_r(ispin)%pw 01863 IF (rho%tau_r_valid) THEN 01864 tau_r(ispin)%pw => rho%tau_r(ispin)%pw 01865 ELSE 01866 CALL pw_pool_create_pw(auxbas_pw_pool,tau_r(ispin)%pw,& 01867 use_data=REALDATA3D,in_space=REALSPACE,error=error) 01868 CALL pw_pool_create_pw(auxbas_pw_pool,tau_g(ispin)%pw,& 01869 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,error=error) 01870 CALL calculate_rho_elec(matrix_p=rho%rho_ao(ispin)%matrix,& 01871 rho=tau_r(ispin),& 01872 rho_gspace=tau_g(ispin),& 01873 total_rho=dum, & 01874 qs_env=qs_env, soft_valid=.FALSE., & 01875 compute_tau=.TRUE., error=error) 01876 END IF 01877 01878 IF(rho%drho_r_valid) THEN 01879 DO idir = 1,3 01880 drho_r(3*(ispin-1)+idir)%pw => rho%drho_r(3*(ispin-1)+idir)%pw 01881 END DO 01882 ELSE 01883 deriv_pw = .FALSE. 01884 ! deriv_pw = .TRUE. 01885 IF(deriv_pw) THEN 01886 nd = RESHAPE ((/1,0,0,0,1,0,0,0,1/),(/3,3/)) 01887 CALL pw_pool_create_pw(auxbas_pw_pool,tmp_g,& 01888 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,error=error) 01889 udvol = 1.0_dp/rho%rho_r(ispin)%pw%pw_grid%dvol 01890 DO idir = 1,3 01891 CALL pw_pool_create_pw(auxbas_pw_pool,drho_r(3*(ispin-1)+idir)%pw,& 01892 use_data=REALDATA3D,in_space=REALSPACE,error=error) 01893 CALL pw_transfer ( rho%rho_r(ispin)%pw, tmp_g , error=error) 01894 CALL pw_derive ( tmp_g, nd(:,idir) , error=error) 01895 CALL pw_transfer (tmp_g, drho_r(3*(ispin-1)+idir)%pw , error=error) 01896 ! CALL pw_scale(drho_r(3*(ispin-1)+idir)%pw,udvol,error=error) 01897 01898 END DO 01899 01900 ELSE 01901 DO idir = 1,3 01902 CALL pw_pool_create_pw(auxbas_pw_pool,drho_r(3*(ispin-1)+idir)%pw,& 01903 use_data=REALDATA3D,in_space=REALSPACE,error=error) 01904 CALL pw_pool_create_pw(auxbas_pw_pool,drho_g(3*(ispin-1)+idir)%pw,& 01905 use_data=COMPLEXDATA1D,in_space=RECIPROCALSPACE,error=error) 01906 CALL calculate_rho_elec(matrix_p=rho%rho_ao(ispin)%matrix,& 01907 rho=drho_r(3*(ispin-1)+idir),& 01908 rho_gspace=drho_g(3*(ispin-1)+idir),& 01909 total_rho=dum, & 01910 qs_env=qs_env, soft_valid=.FALSE., & 01911 compute_tau=.FALSE., compute_grad=.TRUE., idir =idir, error=error) 01912 01913 END DO 01914 END IF 01915 END IF 01916 01917 ! Calculate elf_r 01918 !$omp parallel do default(none) shared(bo,elf_r, ispin, drho_r,rho_r, tau_r, cfermi, rho_cutoff)& 01919 !$omp private(k,j,i, norm_drho, rho_53, elf_kernel) 01920 DO k = bo(1,3), bo(2,3) 01921 DO j = bo(1,2), bo(2,2) 01922 DO i = bo(1,1), bo(2,1) 01923 norm_drho = drho_r(3*(ispin-1)+1)%pw%cr3d(i,j,k)**2+& 01924 drho_r(3*(ispin-1)+2)%pw%cr3d(i,j,k)**2+& 01925 drho_r(3*(ispin-1)+3)%pw%cr3d(i,j,k)**2 01926 norm_drho = norm_drho/MAX( rho_r(ispin)%pw%cr3d(i,j,k),rho_cutoff) 01927 rho_53 = cfermi * MAX(rho_r(ispin)%pw%cr3d(i,j,k),rho_cutoff)**f53 01928 elf_kernel = (tau_r(ispin)%pw%cr3d(i,j,k) - f18*norm_drho) +2.87E-5_dp 01929 elf_kernel = (elf_kernel/rho_53)**2 01930 elf_r(ispin)%pw%cr3d(i,j,k) = 1.0_dp/(1.0_dp+elf_kernel) 01931 IF(elf_r(ispin)%pw%cr3d(i,j,k)<ELFCUT ) elf_r(ispin)%pw%cr3d(i,j,k) =0.0_dp 01932 END DO 01933 END DO 01934 END DO 01935 01936 IF (.NOT. rho%tau_r_valid) THEN 01937 CALL pw_pool_give_back_pw(auxbas_pw_pool,tau_r(ispin)%pw, error=error) 01938 CALL pw_pool_give_back_pw(auxbas_pw_pool,tau_g(ispin)%pw, error=error) 01939 END IF 01940 IF (.NOT. rho%drho_r_valid) THEN 01941 IF(deriv_pw) THEN 01942 CALL pw_pool_give_back_pw(auxbas_pw_pool,tmp_g, error=error) 01943 DO idir = 1,3 01944 CALL pw_pool_give_back_pw(auxbas_pw_pool,drho_r(3*(ispin-1)+idir)%pw, error=error) 01945 END DO 01946 ELSE 01947 DO idir = 1,3 01948 CALL pw_pool_give_back_pw(auxbas_pw_pool,drho_r(3*(ispin-1)+idir)%pw, error=error) 01949 CALL pw_pool_give_back_pw(auxbas_pw_pool,drho_g(3*(ispin-1)+idir)%pw, error=error) 01950 END DO 01951 END IF 01952 END IF 01953 END DO !ispin 01954 01955 DEALLOCATE(rho_r,STAT=stat) 01956 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01957 DEALLOCATE(tau_r,STAT=stat) 01958 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01959 DEALLOCATE(tau_g,STAT=stat) 01960 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01961 DEALLOCATE(drho_r,STAT=stat) 01962 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01963 DEALLOCATE(drho_g,STAT=stat) 01964 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01965 01966 CALL timestop(handle) 01967 01968 END SUBROUTINE qs_elf_calc 01969 01970 ! ***************************************************************************** 01980 SUBROUTINE qs_scf_post_molopt(input, logger, qs_env, error) 01981 TYPE(section_vals_type), POINTER :: input 01982 TYPE(cp_logger_type), POINTER :: logger 01983 TYPE(qs_environment_type), POINTER :: qs_env 01984 TYPE(cp_error_type), INTENT(inout) :: error 01985 01986 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_molopt', 01987 routineP = moduleN//':'//routineN 01988 01989 INTEGER :: handle, nao, unit_nr 01990 REAL(KIND=dp) :: S_cond_number 01991 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues 01992 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01993 POINTER :: matrix_s 01994 TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct 01995 TYPE(cp_fm_type), POINTER :: fm_s, fm_work, mo_coeff 01996 TYPE(mo_set_p_type), DIMENSION(:), 01997 POINTER :: mos 01998 TYPE(qs_energy_type), POINTER :: energy 01999 TYPE(section_vals_type), POINTER :: print_key 02000 02001 CALL timeset(routineN,handle) 02002 02003 print_key => section_vals_get_subs_vals(section_vals=input,& 02004 subsection_name="DFT%PRINT%BASIS_MOLOPT_QUANTITIES",& 02005 error=error) 02006 IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),& 02007 cp_p_file)) THEN 02008 02009 CALL get_qs_env(qs_env,energy=energy,matrix_s=matrix_s,mos=mos,error=error) 02010 02011 ! set up the two needed full matrices, using mo_coeff as a template 02012 CALL get_mo_set(mo_set=mos(1)%mo_set,mo_coeff=mo_coeff,nao=nao) 02013 CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct,& 02014 nrow_global=nao, ncol_global=nao,& 02015 template_fmstruct=mo_coeff%matrix_struct, error=error) 02016 CALL cp_fm_create(fm_s, matrix_struct=ao_ao_fmstruct,& 02017 name="fm_s", error=error) 02018 CALL cp_fm_create(fm_work, matrix_struct=ao_ao_fmstruct,& 02019 name="fm_work", error=error) 02020 CALL cp_fm_struct_release(ao_ao_fmstruct,error=error) 02021 ALLOCATE(eigenvalues(nao)) 02022 02023 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix,fm_s,error=error) 02024 CALL choose_eigv_solver(fm_s,fm_work,eigenvalues,error=error) 02025 02026 CALL cp_fm_release(fm_s,error=error) 02027 CALL cp_fm_release(fm_work,error=error) 02028 02029 S_cond_number=MAXVAL(ABS(eigenvalues))/MAX(MINVAL(ABS(eigenvalues)),EPSILON(0.0_dp)) 02030 02031 unit_nr=cp_print_key_unit_nr(logger,input,"DFT%PRINT%BASIS_MOLOPT_QUANTITIES",& 02032 extension=".molopt",error=error) 02033 02034 IF (unit_nr>0) THEN 02035 ! please keep this format fixed, needs to be grepable for molopt 02036 ! optimizations 02037 WRITE(unit_nr,'(T2,A28,2A25)') "","Tot. Ener.","S Cond. Numb." 02038 WRITE(unit_nr,'(T2,A28,2E25.17)') "BASIS_MOLOPT_QUANTITIES",energy%total,S_cond_number 02039 ENDIF 02040 02041 CALL cp_print_key_finished_output(unit_nr,logger,input,& 02042 "DFT%PRINT%BASIS_MOLOPT_QUANTITIES",error=error) 02043 02044 END IF 02045 02046 CALL timestop(handle) 02047 02048 END SUBROUTINE qs_scf_post_molopt 02049 02050 ! ***************************************************************************** 02056 SUBROUTINE qs_scf_post_epr(input, logger, qs_env, error) 02057 TYPE(section_vals_type), POINTER :: input 02058 TYPE(cp_logger_type), POINTER :: logger 02059 TYPE(qs_environment_type), POINTER :: qs_env 02060 TYPE(cp_error_type), INTENT(inout) :: error 02061 02062 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_epr', 02063 routineP = moduleN//':'//routineN 02064 02065 INTEGER :: handle 02066 TYPE(section_vals_type), POINTER :: print_key 02067 02068 CALL timeset(routineN,handle) 02069 02070 print_key => section_vals_get_subs_vals(section_vals=input,& 02071 subsection_name="DFT%PRINT%HYPERFINE_COUPLING_TENSOR",& 02072 error=error) 02073 IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),& 02074 cp_p_file)) THEN 02075 CALL qs_epr_hyp_calc(qs_env=qs_env,error=error) 02076 END IF 02077 02078 CALL timestop(handle) 02079 02080 END SUBROUTINE qs_scf_post_epr 02081 02082 ! ***************************************************************************** 02087 SUBROUTINE qs_scf_post_gemopt(input, logger, qs_env, error) 02088 TYPE(section_vals_type), POINTER :: input 02089 TYPE(cp_logger_type), POINTER :: logger 02090 TYPE(qs_environment_type), POINTER :: qs_env 02091 TYPE(cp_error_type), INTENT(inout) :: error 02092 02093 CHARACTER(len=*), PARAMETER :: routineN = 'qs_scf_post_gemopt', 02094 routineP = moduleN//':'//routineN 02095 02096 INTEGER :: handle, unit_nr 02097 LOGICAL :: explicit 02098 TYPE(section_vals_type), POINTER :: print_key 02099 02100 CALL timeset(routineN,handle) 02101 02102 02103 print_key => section_vals_get_subs_vals(section_vals=input,& 02104 subsection_name="DFT%PRINT%OPTIMIZE_GEMINALS",& 02105 error=error) 02106 unit_nr=cp_print_key_unit_nr(logger,input,"DFT%PRINT%OPTIMIZE_GEMINALS",& 02107 extension="",error=error) 02108 CALL section_vals_get(print_key,explicit=explicit, error=error) 02109 02110 IF ( explicit) THEN 02111 IF (BTEST(cp_print_key_should_output(logger%iter_info,print_key,error=error),& 02112 cp_p_file)) THEN 02113 IF (unit_nr>0) THEN 02114 WRITE(unit_nr,"(T2,79('-'),/,T2,'-',T30,A,T80,'-',/,T2,79('-'))") " Optimize Geminals " 02115 ENDIF 02116 ! 02117 CALL geminal_optimize(qs_env,print_key,unit_nr,error) 02118 ! 02119 IF (unit_nr>0) WRITE(unit_nr,"(T2,79('-'))") 02120 END IF 02121 END IF 02122 02123 CALL timestop(handle) 02124 02125 END SUBROUTINE qs_scf_post_gemopt 02126 02127 ! ***************************************************************************** 02135 02136 SUBROUTINE write_available_results(qs_env,scf_env,error) 02137 TYPE(qs_environment_type), POINTER :: qs_env 02138 TYPE(qs_scf_env_type), OPTIONAL, POINTER :: scf_env 02139 TYPE(cp_error_type), INTENT(inout) :: error 02140 02141 CHARACTER(len=*), PARAMETER :: routineN = 'write_available_results', 02142 routineP = moduleN//':'//routineN 02143 02144 INTEGER :: handle 02145 02146 CALL timeset(routineN,handle) 02147 CALL write_mo_dependent_results(qs_env,scf_env,error) 02148 CALL write_mo_free_results(qs_env,error) 02149 CALL timestop(handle) 02150 02151 END SUBROUTINE write_available_results 02152 02153 ! ***************************************************************************** 02161 SUBROUTINE write_mo_dependent_results(qs_env,scf_env,error) 02162 TYPE(qs_environment_type), POINTER :: qs_env 02163 TYPE(qs_scf_env_type), OPTIONAL, POINTER :: scf_env 02164 TYPE(cp_error_type), INTENT(inout) :: error 02165 02166 CHARACTER(len=*), PARAMETER :: routineN = 'write_mo_dependent_results', 02167 routineP = moduleN//':'//routineN 02168 02169 INTEGER :: handle, homo, ispin, nmo, 02170 output_unit 02171 LOGICAL :: all_equal, failure 02172 REAL(KIND=dp) :: maxocc, s_square, 02173 s_square_ideal, 02174 total_abs_spin_dens 02175 REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, 02176 occupation_numbers 02177 TYPE(admm_type), POINTER :: admm_env 02178 TYPE(atomic_kind_type), DIMENSION(:), 02179 POINTER :: atomic_kind_set 02180 TYPE(cell_type), POINTER :: cell 02181 TYPE(cp_dbcsr_p_type), DIMENSION(:), 02182 POINTER :: ks_rmpv, ks_rmpv_aux_fit, 02183 matrix_s 02184 TYPE(cp_dbcsr_type), POINTER :: mo_coeff_deriv 02185 TYPE(cp_fm_type), POINTER :: mo_coeff 02186 TYPE(cp_logger_type), POINTER :: logger 02187 TYPE(cp_para_env_type), POINTER :: para_env 02188 TYPE(cp_subsys_type), POINTER :: subsys 02189 TYPE(dft_control_type), POINTER :: dft_control 02190 TYPE(mo_set_p_type), DIMENSION(:), 02191 POINTER :: mos 02192 TYPE(molecule_type), POINTER :: molecule_set( : ) 02193 TYPE(particle_list_type), POINTER :: particles 02194 TYPE(particle_type), DIMENSION(:), 02195 POINTER :: particle_set 02196 TYPE(pw_env_type), POINTER :: pw_env 02197 TYPE(pw_p_type) :: wf_r 02198 TYPE(pw_pool_p_type), DIMENSION(:), 02199 POINTER :: pw_pools 02200 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 02201 TYPE(qs_energy_type), POINTER :: energy 02202 TYPE(qs_rho_type), POINTER :: rho 02203 TYPE(scf_control_type), POINTER :: scf_control 02204 TYPE(section_vals_type), POINTER :: dft_section, input 02205 02206 CALL timeset(routineN,handle) 02207 failure=.FALSE. 02208 NULLIFY(cell, dft_control, pw_env, auxbas_pw_pool, pw_pools, mo_coeff, & 02209 mo_coeff_deriv, mo_eigenvalues, mos, atomic_kind_set, & 02210 particle_set, rho, ks_rmpv, matrix_s, scf_control, dft_section, & 02211 molecule_set, input, particles, subsys) 02212 para_env=>qs_env%para_env 02213 logger => cp_error_get_logger(error) 02214 output_unit= cp_logger_get_default_io_unit(logger) 02215 02216 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 02217 IF (.NOT. failure) THEN 02218 CALL get_qs_env(qs_env,dft_control=dft_control,molecule_set=molecule_set, & 02219 mos=mos,atomic_kind_set=atomic_kind_set,particle_set=particle_set,& 02220 rho=rho,matrix_ks=ks_rmpv,matrix_ks_aux_fit=ks_rmpv_aux_fit,admm_env=admm_env,& 02221 scf_control=scf_control,matrix_s=matrix_s,& 02222 input=input,cell=cell,subsys=subsys,error=error) 02223 CALL cp_subsys_get(subsys,particles=particles,error=error) 02224 02225 ! *** if the dft_section tells you to do so, write last wavefunction to screen 02226 dft_section => section_vals_get_subs_vals(input,"DFT",error=error) 02227 IF(.NOT.qs_env%run_rtp)THEN 02228 IF (dft_control%nspins == 2) THEN 02229 CALL write_mo_set(mos(1)%mo_set,atomic_kind_set,particle_set,4,& 02230 dft_section,spin="ALPHA",last=.TRUE.,error=error) 02231 CALL write_mo_set(mos(2)%mo_set,atomic_kind_set,particle_set,4,& 02232 dft_section,spin="BETA",last=.TRUE.,error=error) 02233 ELSE 02234 CALL write_mo_set(mos(1)%mo_set,atomic_kind_set,particle_set,4,& 02235 dft_section,last=.TRUE.,error=error) 02236 END IF 02237 02238 ! *** at the end of scf print out the projected dos per kind 02239 IF (BTEST(cp_print_key_should_output(logger%iter_info,dft_section,"PRINT%PDOS",& 02240 error=error),cp_p_file) ) THEN 02241 DO ispin = 1,dft_control%nspins 02242 02243 ! ** If we do ADMM, we add have to modify the kohn-sham matrix 02244 IF( dft_control%do_admm ) THEN 02245 CALL admm_correct_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix, & 02246 ks_rmpv_aux_fit(ispin)%matrix, error) 02247 END IF 02248 IF(PRESENT(scf_env))THEN 02249 IF (scf_env%method == ot_method_nr) THEN 02250 CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff=mo_coeff, & 02251 eigenvalues=mo_eigenvalues) 02252 IF (ASSOCIATED(qs_env%mo_derivs)) THEN 02253 mo_coeff_deriv=>qs_env%mo_derivs(ispin)%matrix 02254 ELSE 02255 mo_coeff_deriv=>NULL() 02256 ENDIF 02257 02258 CALL calculate_subspace_eigenvalues(mo_coeff,ks_rmpv(ispin)%matrix,mo_eigenvalues, & 02259 do_rotation=.TRUE.,& 02260 co_rotate_dbcsr=mo_coeff_deriv,error=error) 02261 CALL set_mo_occupation(mo_set=mos(ispin)%mo_set,error=error) 02262 02263 END IF 02264 END IF 02265 IF(dft_control%nspins==2) THEN 02266 CALL calculate_projected_dos(mos(ispin)%mo_set,atomic_kind_set,& 02267 particle_set,qs_env, dft_section,ispin=ispin,error=error) 02268 ELSE 02269 CALL calculate_projected_dos(mos(ispin)%mo_set,atomic_kind_set,& 02270 particle_set,qs_env,dft_section,error=error) 02271 END IF 02272 02273 ! ** If we do ADMM, we add have to modify the kohn-sham matrix 02274 IF( dft_control%do_admm ) THEN 02275 CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_rmpv(ispin)%matrix, & 02276 ks_rmpv_aux_fit(ispin)%matrix, error) 02277 END IF 02278 02279 END DO 02280 ENDIF 02281 END IF 02282 02283 ! *** Integrated absolute spin density and spin contamination *** 02284 IF (dft_control%nspins.eq.2) THEN 02285 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error) 02286 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,& 02287 pw_pools=pw_pools,error=error) 02288 CALL pw_pool_create_pw(auxbas_pw_pool,wf_r%pw,& 02289 use_data = REALDATA3D,& 02290 in_space = REALSPACE, error=error) 02291 CALL pw_copy(rho%rho_r(1)%pw,wf_r%pw, error=error) 02292 CALL pw_axpy(rho%rho_r(2)%pw,wf_r%pw,alpha=-1._dp, error=error) 02293 total_abs_spin_dens=pw_integrate_function(wf_r%pw,oprt="ABS", error=error) 02294 IF (output_unit > 0) WRITE(UNIT=output_unit,FMT='(/,(T3,A,T61,F20.10))')& 02295 "Integrated absolute spin density : ",total_abs_spin_dens 02296 CALL pw_pool_give_back_pw(auxbas_pw_pool,wf_r%pw, error=error) 02297 ! 02298 ! XXX Fix Me XXX 02299 ! should be extended to the case where added MOs are present 02300 ! 02301 all_equal = .TRUE. 02302 DO ispin=1,dft_control%nspins 02303 CALL get_mo_set(mo_set=mos(ispin)%mo_set,& 02304 occupation_numbers=occupation_numbers,& 02305 homo=homo,& 02306 nmo=nmo,& 02307 maxocc=maxocc) 02308 IF (nmo > 0) THEN 02309 all_equal = all_equal.AND.& 02310 (ALL(occupation_numbers(1:homo) == maxocc).AND.& 02311 ALL(occupation_numbers(homo+1:nmo) == 0.0_dp)) 02312 END IF 02313 END DO 02314 IF (.NOT.all_equal) THEN 02315 IF (output_unit>0) WRITE(UNIT=output_unit,FMT="(T3,A)")& 02316 "WARNING: S**2 computation does not yet treat fractional occupied orbitals" 02317 ELSE 02318 CALL get_qs_env(qs_env=qs_env,& 02319 energy=energy,error=error) 02320 CALL compute_s_square(mos=mos, matrix_s=matrix_s, s_square=s_square,& 02321 s_square_ideal=s_square_ideal,error=error) 02322 IF (output_unit > 0) WRITE (UNIT=output_unit,FMT='(T3,A,T51,2F15.6)')& 02323 "Ideal and single determinant S**2 : ",s_square_ideal,s_square 02324 energy%s_square=s_square 02325 END IF 02326 ENDIF 02327 02328 END IF 02329 02330 CALL timestop(handle) 02331 02332 END SUBROUTINE write_mo_dependent_results 02333 02334 ! ***************************************************************************** 02341 SUBROUTINE write_mo_free_results(qs_env,error) 02342 TYPE(qs_environment_type), POINTER :: qs_env 02343 TYPE(cp_error_type), INTENT(inout) :: error 02344 02345 CHARACTER(len=*), PARAMETER :: routineN = 'write_mo_free_results', 02346 routineP = moduleN//':'//routineN 02347 02348 CHARACTER(len=1), DIMENSION(3) :: cdir = (/"x","y","z"/) 02349 CHARACTER(LEN=2) :: element_symbol 02350 CHARACTER(LEN=default_path_length) :: filename, my_pos_cube 02351 CHARACTER(LEN=default_string_length) :: name 02352 INTEGER :: handle, i, iat, id, ikind, 02353 iso, ispin, istat, iw, l, 02354 nd(3), ngto, niso, nkind, np, 02355 nr, output_unit, unit_nr 02356 LOGICAL :: append_cube, failure, 02357 print_total_density, 02358 write_ks, write_xc, 02359 xrd_interface 02360 REAL(KIND=dp) :: q_max, rho_hard, rho_soft, 02361 rho_total, rho_total_rspace, 02362 udvol, volume 02363 REAL(KIND=dp), ALLOCATABLE, 02364 DIMENSION(:, :) :: bfun 02365 REAL(KIND=dp), ALLOCATABLE, 02366 DIMENSION(:, :, :) :: aedens, ccdens, ppdens 02367 REAL(KIND=dp), DIMENSION(3) :: dr 02368 TYPE(atomic_kind_type), DIMENSION(:), 02369 POINTER :: atomic_kind_set 02370 TYPE(atomic_kind_type), POINTER :: atomic_kind 02371 TYPE(cell_type), POINTER :: cell 02372 TYPE(cp_dbcsr_p_type), DIMENSION(:), 02373 POINTER :: ks_rmpv, ks_rmpv_aux_fit, 02374 matrix_s 02375 TYPE(cp_logger_type), POINTER :: logger 02376 TYPE(cp_para_env_type), POINTER :: para_env 02377 TYPE(cp_subsys_type), POINTER :: subsys 02378 TYPE(dft_control_type), POINTER :: dft_control 02379 TYPE(grid_atom_type), POINTER :: grid_atom 02380 TYPE(particle_list_type), POINTER :: particles 02381 TYPE(particle_type), DIMENSION(:), 02382 POINTER :: particle_set 02383 TYPE(pw_env_type), POINTER :: pw_env 02384 TYPE(pw_p_type) :: aux_g, aux_r, 02385 rho_elec_gspace, 02386 rho_elec_rspace, wf_r 02387 TYPE(pw_p_type), POINTER :: rho0_s_gs, rho_core 02388 TYPE(pw_pool_p_type), DIMENSION(:), 02389 POINTER :: pw_pools 02390 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 02391 TYPE(qs_rho_type), POINTER :: rho 02392 TYPE(rho_atom_type), DIMENSION(:), 02393 POINTER :: rho_atom_set 02394 TYPE(rho_atom_type), POINTER :: rho_atom 02395 TYPE(section_vals_type), POINTER :: dft_section, input, xc_section 02396 02397 CALL timeset(routineN,handle) 02398 failure=.FALSE. 02399 NULLIFY(cell, dft_control, pw_env, auxbas_pw_pool, pw_pools, & 02400 atomic_kind_set, particle_set, rho, ks_rmpv, matrix_s, & 02401 dft_section, xc_section, input, particles, subsys) 02402 para_env=>qs_env%para_env 02403 logger => cp_error_get_logger(error) 02404 output_unit= cp_logger_get_default_io_unit(logger) 02405 02406 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 02407 IF (.NOT. failure) THEN 02408 CALL get_qs_env(qs_env,dft_control=dft_control, & 02409 atomic_kind_set=atomic_kind_set,particle_set=particle_set,& 02410 rho=rho,matrix_ks=ks_rmpv,matrix_ks_aux_fit=ks_rmpv_aux_fit,& 02411 matrix_s=matrix_s,input=input,cell=cell,subsys=subsys,error=error) 02412 dft_section => section_vals_get_subs_vals(input,"DFT",error=error) 02413 CALL cp_subsys_get(subsys,particles=particles,error=error) 02414 02415 ! Print the total density (electronic + core charge) 02416 02417 IF (BTEST(cp_print_key_should_output(logger%iter_info,input,& 02418 "DFT%PRINT%TOT_DENSITY_CUBE", error=error),cp_p_file)) THEN 02419 NULLIFY(rho_core,rho0_s_gs) 02420 append_cube = section_get_lval(input,"DFT%PRINT%TOT_DENSITY_CUBE%APPEND",error=error) 02421 my_pos_cube="REWIND" 02422 IF(append_cube) THEN 02423 my_pos_cube="APPEND" 02424 END IF 02425 02426 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,rho_core=rho_core,& 02427 rho0_s_gs=rho0_s_gs,error=error) 02428 CALL pw_env_get(pw_env, error=error) 02429 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error) 02430 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,& 02431 pw_pools=pw_pools,error=error) 02432 CALL pw_pool_create_pw(auxbas_pw_pool,wf_r%pw,& 02433 use_data = REALDATA3D,& 02434 in_space = REALSPACE, error=error) 02435 IF (dft_control%qs_control%gapw) THEN 02436 CALL pw_transfer(rho0_s_gs%pw,wf_r%pw,error=error) 02437 ELSE 02438 CALL pw_transfer(rho_core%pw,wf_r%pw,error=error) 02439 END IF 02440 DO ispin=1,dft_control%nspins 02441 CALL pw_axpy(rho%rho_r(ispin)%pw,wf_r%pw, error=error) 02442 END DO 02443 filename = "TOTAL_DENSITY" 02444 unit_nr=cp_print_key_unit_nr(logger,input,"DFT%PRINT%TOT_DENSITY_CUBE",& 02445 extension=".cube",middle_name=TRIM(filename),file_position=my_pos_cube,& 02446 log_filename=.FALSE.,error=error) 02447 CALL pw_to_cube(wf_r%pw,unit_nr,"TOTAL DENSITY",& 02448 particles=particles,& 02449 stride=section_get_ivals(dft_section,"PRINT%TOT_DENSITY_CUBE%STRIDE",error=error),& 02450 error=error) 02451 CALL cp_print_key_finished_output(unit_nr,logger,input,& 02452 "DFT%PRINT%TOT_DENSITY_CUBE",error=error) 02453 CALL pw_pool_give_back_pw(auxbas_pw_pool,wf_r%pw, error=error) 02454 END IF 02455 02456 ! Write cube file with electron density 02457 IF (BTEST(cp_print_key_should_output(logger%iter_info,input,& 02458 "DFT%PRINT%E_DENSITY_CUBE",error=error),cp_p_file)) THEN 02459 CALL section_vals_val_get(dft_section,& 02460 keyword_name="PRINT%E_DENSITY_CUBE%TOTAL_DENSITY",& 02461 l_val=print_total_density,& 02462 error=error) 02463 append_cube = section_get_lval(input,"DFT%PRINT%E_DENSITY_CUBE%APPEND",error=error) 02464 my_pos_cube="REWIND" 02465 IF(append_cube) THEN 02466 my_pos_cube="APPEND" 02467 END IF 02468 ! Write the info on core densities for the interface between cp2k and the XRD code 02469 ! together with the valence density they are used to compute the form factor (Fourier transform) 02470 xrd_interface = section_get_lval(input,"DFT%PRINT%E_DENSITY_CUBE%XRD_INTERFACE",error=error) 02471 IF(xrd_interface) THEN 02472 !cube file only contains soft density (GAPW) 02473 IF(dft_control%qs_control%gapw) print_total_density = .FALSE. 02474 02475 filename = "ELECTRON_DENSITY" 02476 unit_nr = cp_print_key_unit_nr(logger,input,"DFT%PRINT%E_DENSITY_CUBE",& 02477 extension=".xrd",middle_name=TRIM(filename),& 02478 file_position=my_pos_cube,log_filename=.FALSE.,error=error) 02479 ngto = section_get_ival(input,"DFT%PRINT%E_DENSITY_CUBE%NGAUSS",error=error) 02480 IF (output_unit>0) THEN 02481 INQUIRE (UNIT=unit_nr,NAME=filename) 02482 WRITE (UNIT=output_unit,FMT="(/,T2,A,/,/,T2,A)")& 02483 "The electron density (atomic part) is written to the file:",& 02484 TRIM(filename) 02485 END IF 02486 02487 IF (unit_nr>0) THEN 02488 xc_section => section_vals_get_subs_vals(input,"DFT%XC",error=error) 02489 WRITE(unit_nr,*) "Atomic (core) densities" 02490 WRITE(unit_nr,*) "Unit cell" 02491 WRITE(unit_nr,FMT="(3F20.12)") cell%hmat(1,1),cell%hmat(1,2),cell%hmat(1,3) 02492 WRITE(unit_nr,FMT="(3F20.12)") cell%hmat(2,1),cell%hmat(2,2),cell%hmat(2,3) 02493 WRITE(unit_nr,FMT="(3F20.12)") cell%hmat(3,1),cell%hmat(3,2),cell%hmat(3,3) 02494 ! calculate atomic density and core density 02495 nkind = SIZE(atomic_kind_set) 02496 WRITE(unit_nr,*) "Atomic types" 02497 WRITE(unit_nr,*) nkind 02498 ALLOCATE(ppdens(ngto,2,nkind),aedens(ngto,2,nkind),ccdens(ngto,2,nkind),stat=istat) 02499 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 02500 DO ikind=1,nkind 02501 atomic_kind => atomic_kind_set(ikind) 02502 CALL get_atomic_kind(atomic_kind=atomic_kind,name=name,element_symbol=element_symbol) 02503 WRITE(unit_nr,FMT="(I6,A10,A20)") ikind,TRIM(element_symbol),TRIM(name) 02504 WRITE(unit_nr,FMT="(I6)") ngto 02505 CALL calculate_atomic_density(ppdens(:,:,ikind),atomic_kind,ngto,iunit=output_unit,& 02506 confine=.TRUE.,error=error) 02507 CALL calculate_atomic_density(aedens(:,:,ikind),atomic_kind,ngto,iunit=output_unit,& 02508 allelectron=.TRUE.,confine=.TRUE.,error=error) 02509 WRITE(unit_nr,*) " Total density" 02510 WRITE(unit_nr,FMT="(2G24.12)") (aedens(i,1,ikind),aedens(i,2,ikind),i=1,ngto) 02511 ccdens(:,1,ikind) = aedens(:,1,ikind) 02512 ccdens(:,2,ikind) = 0._dp 02513 CALL project_function_a(ccdens(1:ngto,2,ikind),ccdens(1:ngto,1,ikind),& 02514 ppdens(1:ngto,2,ikind),ppdens(1:ngto,1,ikind),0,error) 02515 ccdens(:,2,ikind) = aedens(:,2,ikind) - ccdens(:,2,ikind) 02516 WRITE(unit_nr,*) " Core density" 02517 WRITE(unit_nr,FMT="(2G24.12)") (ccdens(i,1,ikind),ccdens(i,2,ikind),i=1,ngto) 02518 NULLIFY(atomic_kind) 02519 END DO 02520 END IF 02521 02522 IF (dft_control%qs_control%gapw) THEN 02523 CALL get_qs_env(qs_env=qs_env,rho_atom_set=rho_atom_set,error=error) 02524 02525 IF (unit_nr>0) THEN 02526 WRITE(unit_nr,*) "Coordinates and GAPW density" 02527 np = particles%n_els 02528 DO iat=1,np 02529 CALL get_atomic_kind(particles%els(iat)%atomic_kind,kind_number=ikind,& 02530 grid_atom=grid_atom) 02531 rho_atom => rho_atom_set(iat) 02532 nr = SIZE(rho_atom%rho_rad_h(1)%r_coef,1) 02533 niso = SIZE(rho_atom%rho_rad_h(1)%r_coef,2) 02534 02535 ALLOCATE(bfun(nr,niso),stat=istat) 02536 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 02537 bfun = 0._dp 02538 DO ispin = 1,dft_control%nspins 02539 bfun = bfun + rho_atom%rho_rad_h(ispin)%r_coef - rho_atom%rho_rad_s(ispin)%r_coef 02540 END DO 02541 ccdens(:,1,ikind) = ppdens(:,1,ikind) 02542 ccdens(:,2,ikind) = 0._dp 02543 WRITE(unit_nr,'(I10,I5,3f12.6)') iat,ikind,particles%els(iat)%r 02544 DO iso=1,niso 02545 l = indso(1,iso) 02546 CALL project_function_b(ccdens(:,2,ikind),ccdens(:,1,ikind),bfun(:,iso),grid_atom,l,error) 02547 WRITE(unit_nr,FMT="(3I6)") iso,l,ngto 02548 WRITE(unit_nr,FMT="(2G24.12)") (ccdens(i,1,ikind),ccdens(i,2,ikind),i=1,ngto) 02549 END DO 02550 DEALLOCATE(bfun,stat=istat) 02551 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 02552 END DO 02553 END IF 02554 ELSE 02555 IF (unit_nr>0) THEN 02556 WRITE(unit_nr,*) "Coordinates" 02557 np = particles%n_els 02558 DO iat=1,np 02559 CALL get_atomic_kind(particles%els(iat)%atomic_kind,kind_number=ikind) 02560 WRITE(unit_nr,'(I10,I5,3f12.6)') iat,ikind,particles%els(iat)%r 02561 END DO 02562 END IF 02563 END IF 02564 02565 IF (unit_nr>0) THEN 02566 DEALLOCATE(ppdens,aedens,ccdens,stat=istat) 02567 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 02568 END IF 02569 02570 CALL cp_print_key_finished_output(unit_nr,logger,input,& 02571 "DFT%PRINT%E_DENSITY_CUBE",& 02572 error=error) 02573 02574 END IF 02575 IF (dft_control%qs_control%gapw.AND.print_total_density) THEN 02576 ! Print total electronic density 02577 CALL get_qs_env(qs_env=qs_env,& 02578 pw_env=pw_env,& 02579 error=error) 02580 CALL pw_env_get(pw_env=pw_env,& 02581 auxbas_pw_pool=auxbas_pw_pool,& 02582 pw_pools=pw_pools,& 02583 error=error) 02584 CALL pw_pool_create_pw(pool=auxbas_pw_pool,& 02585 pw=rho_elec_rspace%pw,& 02586 use_data=REALDATA3D,& 02587 in_space=REALSPACE,& 02588 error=error) 02589 CALL pw_zero(rho_elec_rspace%pw,error=error) 02590 CALL pw_pool_create_pw(pool=auxbas_pw_pool,& 02591 pw=rho_elec_gspace%pw,& 02592 use_data=COMPLEXDATA1D,& 02593 in_space=RECIPROCALSPACE,& 02594 error=error) 02595 CALL pw_zero(rho_elec_gspace%pw,error=error) 02596 CALL get_pw_grid_info(pw_grid=rho_elec_gspace%pw%pw_grid,& 02597 dr=dr,& 02598 vol=volume,& 02599 error=error) 02600 q_max = SQRT(SUM((pi/dr(:))**2)) 02601 CALL calculate_rhotot_elec_gspace(qs_env=qs_env,& 02602 auxbas_pw_pool=auxbas_pw_pool,& 02603 rhotot_elec_gspace=rho_elec_gspace,& 02604 q_max=q_max,& 02605 rho_hard=rho_hard,& 02606 rho_soft=rho_soft,& 02607 error=error) 02608 rho_total = rho_hard + rho_soft 02609 CALL get_pw_grid_info(pw_grid=rho_elec_gspace%pw%pw_grid,& 02610 vol=volume,& 02611 error=error) 02612 CALL pw_transfer(rho_elec_gspace%pw,rho_elec_rspace%pw,debug=.FALSE.,error=error) 02613 rho_total_rspace = pw_integrate_function(rho_elec_rspace%pw,isign=-1,error=error)/volume 02614 filename = "TOTAL_ELECTRON_DENSITY" 02615 unit_nr = cp_print_key_unit_nr(logger,input,"DFT%PRINT%E_DENSITY_CUBE",& 02616 extension=".cube",middle_name=TRIM(filename),& 02617 file_position=my_pos_cube,log_filename=.FALSE.,error=error) 02618 IF (output_unit>0) THEN 02619 INQUIRE (UNIT=unit_nr,NAME=filename) 02620 WRITE (UNIT=output_unit,FMT="(/,T2,A,/,/,T2,A)")& 02621 "The total electron density is written in cube file format to the file:",& 02622 TRIM(filename) 02623 WRITE (UNIT=output_unit,FMT="(/,(T2,A,F20.10))")& 02624 "q(max) [1/Angstrom] :",q_max/angstrom,& 02625 "Soft electronic charge (G-space) :",rho_soft,& 02626 "Hard electronic charge (G-space) :",rho_hard,& 02627 "Total electronic charge (G-space):",rho_total,& 02628 "Total electronic charge (R-space):",rho_total_rspace 02629 END IF 02630 CALL pw_to_cube(rho_elec_rspace%pw,unit_nr,"TOTAL ELECTRON DENSITY",& 02631 particles=particles,& 02632 stride=section_get_ivals(dft_section,"PRINT%E_DENSITY_CUBE%STRIDE",error=error),& 02633 error=error) 02634 CALL cp_print_key_finished_output(unit_nr,logger,input,& 02635 "DFT%PRINT%E_DENSITY_CUBE",& 02636 error=error) 02637 ! Print total spin density for spin-polarized systems 02638 IF (dft_control%nspins > 1) THEN 02639 CALL pw_zero(rho_elec_gspace%pw,error=error) 02640 CALL pw_zero(rho_elec_rspace%pw,error=error) 02641 CALL calculate_rhotot_elec_gspace(qs_env=qs_env,& 02642 auxbas_pw_pool=auxbas_pw_pool,& 02643 rhotot_elec_gspace=rho_elec_gspace,& 02644 q_max=q_max,& 02645 rho_hard=rho_hard,& 02646 rho_soft=rho_soft,& 02647 fsign=-1.0_dp,& 02648 error=error) 02649 rho_total = rho_hard + rho_soft 02650 CALL pw_transfer(rho_elec_gspace%pw,rho_elec_rspace%pw,debug=.FALSE.,error=error) 02651 rho_total_rspace = pw_integrate_function(rho_elec_rspace%pw,isign=-1,error=error)/volume 02652 filename = "TOTAL_SPIN_DENSITY" 02653 unit_nr = cp_print_key_unit_nr(logger,input,"DFT%PRINT%E_DENSITY_CUBE",& 02654 extension=".cube",middle_name=TRIM(filename),& 02655 file_position=my_pos_cube,log_filename=.FALSE.,error=error) 02656 IF (output_unit>0) THEN 02657 INQUIRE (UNIT=unit_nr,NAME=filename) 02658 WRITE (UNIT=output_unit,FMT="(/,T2,A,/,/,T2,A)")& 02659 "The total spin density is written in cube file format to the file:",& 02660 TRIM(filename) 02661 WRITE (UNIT=output_unit,FMT="(/,(T2,A,F20.10))")& 02662 "q(max) [1/Angstrom] :",q_max/angstrom,& 02663 "Soft part of the spin density (G-space):",rho_soft,& 02664 "Hard part of the spin density (G-space):",rho_hard,& 02665 "Total spin density (G-space) :",rho_total,& 02666 "Total spin density (R-space) :",rho_total_rspace 02667 END IF 02668 CALL pw_to_cube(rho_elec_rspace%pw,unit_nr,"TOTAL SPIN DENSITY",& 02669 particles=particles,& 02670 stride=section_get_ivals(dft_section,"PRINT%E_DENSITY_CUBE%STRIDE",error=error),& 02671 error=error) 02672 END IF 02673 CALL pw_pool_give_back_pw(auxbas_pw_pool,rho_elec_gspace%pw,error=error) 02674 CALL pw_pool_give_back_pw(auxbas_pw_pool,rho_elec_rspace%pw,error=error) 02675 ELSE 02676 IF (dft_control%nspins > 1) THEN 02677 CALL get_qs_env(qs_env=qs_env,& 02678 pw_env=pw_env,& 02679 error=error) 02680 CALL pw_env_get(pw_env=pw_env,& 02681 auxbas_pw_pool=auxbas_pw_pool,& 02682 pw_pools=pw_pools,& 02683 error=error) 02684 CALL pw_pool_create_pw(pool=auxbas_pw_pool,& 02685 pw=rho_elec_rspace%pw,& 02686 use_data=REALDATA3D,& 02687 in_space=REALSPACE,& 02688 error=error) 02689 CALL pw_copy(rho%rho_r(1)%pw,rho_elec_rspace%pw, error=error) 02690 CALL pw_axpy(rho%rho_r(2)%pw,rho_elec_rspace%pw, error=error) 02691 filename = "ELECTRON_DENSITY" 02692 unit_nr = cp_print_key_unit_nr(logger,input,"DFT%PRINT%E_DENSITY_CUBE",& 02693 extension=".cube",middle_name=TRIM(filename),& 02694 file_position=my_pos_cube,log_filename=.FALSE.,error=error) 02695 IF (output_unit>0) THEN 02696 INQUIRE (UNIT=unit_nr,NAME=filename) 02697 WRITE (UNIT=output_unit,FMT="(/,T2,A,/,/,T2,A)")& 02698 "The sum of alpha and beta density is written in cube file format to the file:",& 02699 TRIM(filename) 02700 END IF 02701 CALL pw_to_cube(rho_elec_rspace%pw,unit_nr,"SUM OF ALPHA AND BETA DENSITY",& 02702 particles=particles,stride=section_get_ivals(dft_section,"PRINT%E_DENSITY_CUBE%STRIDE",& 02703 error=error),error=error) 02704 CALL cp_print_key_finished_output(unit_nr,logger,input,& 02705 "DFT%PRINT%E_DENSITY_CUBE",error=error) 02706 CALL pw_copy(rho%rho_r(1)%pw,rho_elec_rspace%pw, error=error) 02707 CALL pw_axpy(rho%rho_r(2)%pw,rho_elec_rspace%pw,alpha=-1.0_dp, error=error) 02708 filename = "SPIN_DENSITY" 02709 unit_nr = cp_print_key_unit_nr(logger,input,"DFT%PRINT%E_DENSITY_CUBE",& 02710 extension=".cube",middle_name=TRIM(filename),& 02711 file_position=my_pos_cube,log_filename=.FALSE.,error=error) 02712 IF (output_unit>0) THEN 02713 INQUIRE (UNIT=unit_nr,NAME=filename) 02714 WRITE (UNIT=output_unit,FMT="(/,T2,A,/,/,T2,A)")& 02715 "The spin density is written in cube file format to the file:",& 02716 TRIM(filename) 02717 END IF 02718 CALL pw_to_cube(rho_elec_rspace%pw,unit_nr,"SPIN DENSITY",& 02719 particles=particles,& 02720 stride=section_get_ivals(dft_section,"PRINT%E_DENSITY_CUBE%STRIDE",error=error),& 02721 error=error) 02722 CALL cp_print_key_finished_output(unit_nr,logger,input,& 02723 "DFT%PRINT%E_DENSITY_CUBE",& 02724 error=error) 02725 CALL pw_pool_give_back_pw(auxbas_pw_pool,rho_elec_rspace%pw,error=error) 02726 ELSE 02727 filename = "ELECTRON_DENSITY" 02728 unit_nr = cp_print_key_unit_nr(logger,input,"DFT%PRINT%E_DENSITY_CUBE",& 02729 extension=".cube",middle_name=TRIM(filename),& 02730 file_position=my_pos_cube,log_filename=.FALSE.,error=error) 02731 IF (output_unit>0) THEN 02732 INQUIRE (UNIT=unit_nr,NAME=filename) 02733 WRITE (UNIT=output_unit,FMT="(/,T2,A,/,/,T2,A)")& 02734 "The electron density is written in cube file format to the file:",& 02735 TRIM(filename) 02736 END IF 02737 CALL pw_to_cube(rho%rho_r(1)%pw,unit_nr,"ELECTRON DENSITY",& 02738 particles=particles,& 02739 stride=section_get_ivals(dft_section,"PRINT%E_DENSITY_CUBE%STRIDE",error=error),& 02740 error=error) 02741 CALL cp_print_key_finished_output(unit_nr,logger,input,& 02742 "DFT%PRINT%E_DENSITY_CUBE",& 02743 error=error) 02744 END IF ! nspins 02745 END IF ! total density for GAPW 02746 END IF ! print key 02747 02748 02749 ! Print the hartree potential 02750 IF (BTEST(cp_print_key_should_output(logger%iter_info,input,& 02751 "DFT%PRINT%V_HARTREE_CUBE",error=error),cp_p_file)) THEN 02752 02753 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error) 02754 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,error=error) 02755 CALL pw_pool_create_pw(auxbas_pw_pool,aux_r%pw,& 02756 use_data = REALDATA3D,& 02757 in_space = REALSPACE, error=error) 02758 02759 append_cube = section_get_lval(input,"DFT%PRINT%V_HARTREE_CUBE%APPEND",error=error) 02760 my_pos_cube="REWIND" 02761 IF(append_cube) THEN 02762 my_pos_cube="APPEND" 02763 END IF 02764 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error) 02765 CALL pw_env_get(pw_env, error=error) 02766 unit_nr=cp_print_key_unit_nr(logger,input,"DFT%PRINT%V_HARTREE_CUBE",& 02767 extension=".cube",middle_name="v_hartree",file_position=my_pos_cube,error=error) 02768 udvol = 1.0_dp/qs_env%ks_env%v_hartree_rspace%pw%pw_grid%dvol 02769 02770 CALL pw_copy(qs_env%ks_env%v_hartree_rspace%pw,aux_r%pw, error=error) 02771 CALL pw_scale(aux_r%pw,udvol,error=error) 02772 02773 CALL pw_to_cube(aux_r%pw,unit_nr,"HARTREE POTENTIAL",particles=particles,& 02774 stride=section_get_ivals(dft_section,& 02775 "PRINT%V_HARTREE_CUBE%STRIDE",error=error),& 02776 error=error) 02777 CALL cp_print_key_finished_output(unit_nr,logger,input,& 02778 "DFT%PRINT%V_HARTREE_CUBE",error=error) 02779 02780 CALL pw_pool_give_back_pw(auxbas_pw_pool,aux_r%pw, error=error) 02781 ENDIF 02782 02783 ! Print the Electrical Field Components 02784 IF (BTEST(cp_print_key_should_output(logger%iter_info,input,& 02785 "DFT%PRINT%EFIELD_CUBE",error=error),cp_p_file)) THEN 02786 02787 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error) 02788 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool,error=error) 02789 CALL pw_pool_create_pw(auxbas_pw_pool,aux_r%pw,& 02790 use_data = REALDATA3D,& 02791 in_space = REALSPACE, error=error) 02792 CALL pw_pool_create_pw(auxbas_pw_pool,aux_g%pw,& 02793 use_data = COMPLEXDATA1D,& 02794 in_space = RECIPROCALSPACE, error=error) 02795 02796 append_cube = section_get_lval(input,"DFT%PRINT%EFIELD_CUBE%APPEND",error=error) 02797 my_pos_cube="REWIND" 02798 IF(append_cube) THEN 02799 my_pos_cube="APPEND" 02800 END IF 02801 CALL get_qs_env(qs_env=qs_env,pw_env=pw_env,error=error) 02802 CALL pw_env_get(pw_env, error=error) 02803 udvol = 1.0_dp/qs_env%ks_env%v_hartree_rspace%pw%pw_grid%dvol 02804 DO id=1,3 02805 unit_nr=cp_print_key_unit_nr(logger,input,"DFT%PRINT%EFIELD_CUBE",& 02806 extension=".cube",middle_name="efield_"//cdir(id),file_position=my_pos_cube,error=error) 02807 02808 CALL pw_transfer(qs_env%ks_env%v_hartree_rspace%pw,aux_g%pw, error=error) 02809 nd=0 02810 nd(id)=1 02811 CALL pw_derive(aux_g%pw,nd,error=error) 02812 CALL pw_transfer(aux_g%pw,aux_r%pw, error=error) 02813 CALL pw_scale(aux_r%pw,udvol,error=error) 02814 02815 CALL pw_to_cube(aux_r%pw,& 02816 unit_nr,"ELECTRIC FIELD",particles=particles,& 02817 stride=section_get_ivals(dft_section,& 02818 "PRINT%EFIELD_CUBE%STRIDE",error=error),& 02819 error=error) 02820 CALL cp_print_key_finished_output(unit_nr,logger,input,& 02821 "DFT%PRINT%EFIELD_CUBE",error=error) 02822 END DO 02823 02824 CALL pw_pool_give_back_pw(auxbas_pw_pool,aux_r%pw, error=error) 02825 CALL pw_pool_give_back_pw(auxbas_pw_pool,aux_g%pw, error=error) 02826 END IF 02827 02828 ! Write the density matrices 02829 IF (BTEST(cp_print_key_should_output(logger%iter_info,input,& 02830 "DFT%PRINT%AO_MATRICES/DENSITY",error=error),cp_p_file)) THEN 02831 iw = cp_print_key_unit_nr(logger,input,"DFT%PRINT%AO_MATRICES/DENSITY",& 02832 extension=".Log",error=error) 02833 DO ispin=1,dft_control%nspins 02834 CALL cp_dbcsr_write_sparse_matrix(rho%rho_ao(ispin)%matrix,4,6,qs_env,& 02835 para_env,output_unit=iw,error=error) 02836 END DO 02837 CALL cp_print_key_finished_output(iw,logger,input,& 02838 "DFT%PRINT%AO_MATRICES/DENSITY",& 02839 error=error) 02840 END IF 02841 02842 ! Write the Kohn-Sham matrices 02843 write_ks=BTEST(cp_print_key_should_output(logger%iter_info,input,& 02844 "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX",error=error),cp_p_file) 02845 write_xc=BTEST(cp_print_key_should_output(logger%iter_info,input,& 02846 "DFT%PRINT%AO_MATRICES/MATRIX_VXC",error=error),cp_p_file) 02847 ! we need to update stuff before writing, potentially computing the matrix_vxc 02848 IF (write_ks .OR. write_xc) THEN 02849 IF (write_xc) qs_env%requires_matrix_vxc=.TRUE. 02850 CALL qs_ks_did_change(qs_env%ks_env,rho_changed=.TRUE.,error=error) 02851 CALL qs_ks_update_qs_env(qs_env%ks_env,qs_env=qs_env,error=error,& 02852 calculate_forces=.FALSE.,just_energy=.FALSE.) 02853 IF (write_xc) qs_env%requires_matrix_vxc=.FALSE. 02854 END IF 02855 02856 ! Write the Kohn-Sham matrices 02857 IF (write_ks) THEN 02858 iw = cp_print_key_unit_nr(logger,input,"DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX",& 02859 extension=".Log",error=error) 02860 DO ispin=1,dft_control%nspins 02861 CALL cp_dbcsr_write_sparse_matrix(ks_rmpv(ispin)%matrix,4,6,qs_env,& 02862 para_env,output_unit=iw,error=error) 02863 END DO 02864 CALL cp_print_key_finished_output(iw,logger,input,& 02865 "DFT%PRINT%AO_MATRICES/KOHN_SHAM_MATRIX",& 02866 error=error) 02867 END IF 02868 02869 ! Write the xc matrix 02870 IF (write_xc) THEN 02871 iw = cp_print_key_unit_nr(logger,input,"DFT%PRINT%AO_MATRICES/MATRIX_VXC",& 02872 extension=".Log",error=error) 02873 DO ispin=1,dft_control%nspins 02874 CALL cp_dbcsr_write_sparse_matrix(qs_env%matrix_vxc(ispin)%matrix,4,6,qs_env,& 02875 para_env,output_unit=iw,error=error) 02876 END DO 02877 CALL cp_print_key_finished_output(iw,logger,input,& 02878 "DFT%PRINT%AO_MATRICES/MATRIX_VXC",& 02879 error=error) 02880 END IF 02881 02882 END IF 02883 02884 CALL timestop(handle) 02885 02886 END SUBROUTINE write_mo_free_results 02887 02888 SUBROUTINE project_function_a(ca,a,cb,b,l,error) 02889 ! project function cb on ca 02890 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: ca 02891 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a, cb, b 02892 INTEGER, INTENT(IN) :: l 02893 TYPE(cp_error_type), INTENT(inout) :: error 02894 02895 CHARACTER(len=*), PARAMETER :: routineN = 'project_function_a', 02896 routineP = moduleN//':'//routineN 02897 02898 INTEGER :: info, istat, n 02899 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv 02900 LOGICAL :: failure = .FALSE. 02901 REAL(KIND=dp), ALLOCATABLE, 02902 DIMENSION(:, :) :: smat, tmat, v 02903 02904 n = SIZE(ca) 02905 ALLOCATE(smat(n,n),tmat(n,n),v(n,1),ipiv(n),stat=istat) 02906 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 02907 02908 CALL sg_overlap ( smat, l, a, a ) 02909 CALL sg_overlap ( tmat, l, a, b ) 02910 v(:,1) = MATMUL(tmat,cb) 02911 CALL lapack_sgesv ( n, 1, smat, n, ipiv, v, n, info ) 02912 CPPostcondition(info==0,cp_failure_level,routineP,error,failure) 02913 ca(:) = v(:,1) 02914 02915 DEALLOCATE(smat,tmat,v,ipiv,stat=istat) 02916 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02917 02918 END SUBROUTINE project_function_a 02919 02920 SUBROUTINE project_function_b(ca,a,bfun,grid_atom,l,error) 02921 ! project function f on ca 02922 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: ca 02923 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a, bfun 02924 TYPE(grid_atom_type), POINTER :: grid_atom 02925 INTEGER, INTENT(IN) :: l 02926 TYPE(cp_error_type), INTENT(inout) :: error 02927 02928 CHARACTER(len=*), PARAMETER :: routineN = 'project_function_b', 02929 routineP = moduleN//':'//routineN 02930 02931 INTEGER :: i, info, istat, n, nr 02932 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv 02933 LOGICAL :: failure = .FALSE. 02934 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: afun 02935 REAL(KIND=dp), ALLOCATABLE, 02936 DIMENSION(:, :) :: smat, v 02937 02938 n = SIZE(ca) 02939 nr = grid_atom%nr 02940 ALLOCATE(smat(n,n),v(n,1),ipiv(n),afun(nr),stat=istat) 02941 CPPrecondition(istat==0,cp_failure_level,routineP,error,failure) 02942 02943 CALL sg_overlap ( smat, l, a, a ) 02944 DO i=1,n 02945 afun(:) = grid_atom%rad(:)**l * EXP(-a(i)*grid_atom%rad2(:)) 02946 v(i,1) = SUM(afun(:)*bfun(:)*grid_atom%wr(:)) 02947 END DO 02948 CALL lapack_sgesv ( n, 1, smat, n, ipiv, v, n, info ) 02949 CPPostcondition(info==0,cp_failure_level,routineP,error,failure) 02950 ca(:) = v(:,1) 02951 02952 DEALLOCATE(smat,v,ipiv,afun,stat=istat) 02953 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02954 02955 END SUBROUTINE project_function_b 02956 02957 END MODULE qs_scf_post_gpw
1.7.3