|
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 ! ***************************************************************************** 00040 MODULE qs_scf 00041 USE array_types, ONLY: array_i1d_obj,& 00042 array_release 00043 USE atomic_kind_types, ONLY: atomic_kind_type,& 00044 get_atomic_kind,& 00045 get_atomic_kind_set,& 00046 set_atomic_kind 00047 USE cp_control_types, ONLY: dft_control_type 00048 USE cp_dbcsr_interface, ONLY: & 00049 cp_create_bl_distribution, cp_dbcsr_col_block_sizes, cp_dbcsr_copy, & 00050 cp_dbcsr_create, cp_dbcsr_distribution, cp_dbcsr_distribution_release, & 00051 cp_dbcsr_init, cp_dbcsr_init_p, cp_dbcsr_name, & 00052 cp_dbcsr_row_block_sizes, cp_dbcsr_set 00053 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 00054 copy_fm_to_dbcsr,& 00055 cp_dbcsr_alloc_block_from_nbl,& 00056 cp_dbcsr_allocate_matrix_set,& 00057 cp_dbcsr_deallocate_matrix_set,& 00058 cp_dbcsr_sm_fm_multiply 00059 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix,& 00060 write_fm_with_basis_info 00061 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type,& 00062 cp_dbcsr_type 00063 USE cp_external_control, ONLY: external_control 00064 USE cp_fm_basic_linalg, ONLY: cp_fm_triangular_invert 00065 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose 00066 USE cp_fm_diag, ONLY: cp_fm_power 00067 USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,& 00068 fm_pool_get_el_struct 00069 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00070 cp_fm_struct_get,& 00071 cp_fm_struct_release,& 00072 cp_fm_struct_type 00073 USE cp_fm_types, ONLY: cp_fm_create,& 00074 cp_fm_get_info,& 00075 cp_fm_p_type,& 00076 cp_fm_release,& 00077 cp_fm_to_fm,& 00078 cp_fm_type 00079 USE cp_output_handling, ONLY: cp_add_iter_level,& 00080 cp_iterate,& 00081 cp_p_file,& 00082 cp_print_key_finished_output,& 00083 cp_print_key_should_output,& 00084 cp_print_key_unit_nr,& 00085 cp_rm_iter_level 00086 USE cp_para_types, ONLY: cp_para_env_type 00087 USE dbcsr_methods, ONLY: dbcsr_distribution_mp,& 00088 dbcsr_distribution_new,& 00089 dbcsr_distribution_row_dist,& 00090 dbcsr_mp_npcols 00091 USE dbcsr_types, ONLY: dbcsr_distribution_obj,& 00092 dbcsr_type_no_symmetry,& 00093 dbcsr_type_real_default,& 00094 dbcsr_type_symmetric 00095 USE f77_blas 00096 USE harris_env_types, ONLY: harris_env_type 00097 USE harris_functional, ONLY: harris_eigenvalue_calculation,& 00098 harris_eigenvalue_trace_KS_Pmix,& 00099 harris_energy_correction,& 00100 harris_postprocessing 00101 USE input_constants, ONLY: & 00102 broy_mix, broy_mix_new, cholesky_inverse, cholesky_off, & 00103 densities_guess, diag_block_davidson, diag_block_krylov, diag_ot, & 00104 diag_standard, direct_p_mix, history_guess, kerker_mix, multisec_mix, & 00105 no_mix, ot_precond_full_all, ot_precond_full_single, & 00106 ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, & 00107 ot_precond_sparse_diag, outer_scf_none, outer_scf_scp, plus_u_lowdin, & 00108 pulay_mix, wfi_frozen_method_nr, wfi_use_guess_method_nr 00109 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00110 section_vals_type,& 00111 section_vals_val_get 00112 USE kahan_sum, ONLY: accurate_sum 00113 USE kinds, ONLY: default_string_length,& 00114 dp 00115 USE machine, ONLY: m_flush,& 00116 m_walltime 00117 USE mp2, ONLY: mp2_main 00118 USE particle_types, ONLY: particle_type 00119 USE physcon, ONLY: evolt,& 00120 kcalmol 00121 USE preconditioner, ONLY: prepare_preconditioner,& 00122 restart_preconditioner 00123 USE preconditioner_types, ONLY: destroy_preconditioner 00124 USE pw_env_types, ONLY: pw_env_get 00125 USE pw_pool_types, ONLY: pw_pool_give_back_pw,& 00126 pw_pool_type 00127 USE qmmm_image_charge, ONLY: calculate_image_matrix,& 00128 print_image_coefficients 00129 USE qs_block_davidson_types, ONLY: block_davidson_allocate,& 00130 block_davidson_deallocate,& 00131 block_davidson_env_create 00132 USE qs_charges_types, ONLY: qs_charges_type 00133 USE qs_density_mixing_types, ONLY: & 00134 broyden_mixing_new_nr, broyden_mixing_nr, direct_mixing_nr, & 00135 gspace_mixing_nr, mixing_storage_create, mixing_storage_release, & 00136 mixing_storage_type, multisecant_mixing_nr, no_mixing_nr, & 00137 pulay_mixing_nr 00138 USE qs_diis, ONLY: qs_diis_b_clear,& 00139 qs_diis_b_create 00140 USE qs_energy_types, ONLY: qs_energy_type 00141 USE qs_environment_types, ONLY: get_qs_env,& 00142 qs_env_reorthogonalize_vectors,& 00143 qs_environment_type,& 00144 set_qs_env 00145 USE qs_gspace_mixing, ONLY: gspace_mixing,& 00146 mixing_allocate,& 00147 mixing_init,& 00148 self_consistency_check 00149 USE qs_initial_guess, ONLY: calculate_first_density_matrix 00150 USE qs_ks_methods, ONLY: qs_ks_create,& 00151 qs_ks_did_change,& 00152 qs_ks_update_qs_env 00153 USE qs_ks_scp_methods, ONLY: qs_ks_scp_did_change,& 00154 qs_ks_scp_update 00155 USE qs_ks_types, ONLY: qs_ks_env_type,& 00156 qs_ks_release 00157 USE qs_matrix_pools, ONLY: mpools_get 00158 USE qs_mo_methods, ONLY: calculate_density_matrix,& 00159 calculate_magnitude,& 00160 calculate_orthonormality,& 00161 make_basis_simple,& 00162 make_basis_sm 00163 USE qs_mo_types, ONLY: get_mo_set,& 00164 init_mo_set,& 00165 mo_set_p_type,& 00166 set_mo_occupation,& 00167 set_mo_set,& 00168 write_mo_set 00169 USE qs_ot, ONLY: qs_ot_new_preconditioner 00170 USE qs_ot_scf, ONLY: ot_scf_destroy,& 00171 ot_scf_init,& 00172 ot_scf_mini,& 00173 ot_scf_read_input 00174 USE qs_outer_scf, ONLY: outer_loop_extrapolate,& 00175 outer_loop_gradient,& 00176 outer_loop_optimize,& 00177 outer_loop_update_qs_env,& 00178 outer_loop_variables_count 00179 USE qs_rho_atom_types, ONLY: rho_atom_type 00180 USE qs_rho_methods, ONLY: duplicate_rho_type,& 00181 qs_rho_update_rho 00182 USE qs_rho_types, ONLY: qs_rho_type 00183 USE qs_scf_diagonalization, ONLY: diag_subspace_allocate,& 00184 do_block_davidson_diag,& 00185 do_block_krylov_diag,& 00186 do_general_diag,& 00187 do_ot_diag,& 00188 do_roks_diag,& 00189 do_scf_diag_subspace,& 00190 do_special_diag 00191 USE qs_scf_lanczos, ONLY: krylov_space_allocate 00192 USE qs_scf_methods, ONLY: scf_env_density_mixing 00193 USE qs_scf_post_dftb, ONLY: scf_post_calculation_dftb 00194 USE qs_scf_post_gpw, ONLY: scf_post_calculation_gpw 00195 USE qs_scf_post_se, ONLY: scf_post_calculation_se 00196 USE qs_scf_types, ONLY: & 00197 block_davidson_diag_method_nr, block_krylov_diag_method_nr, & 00198 diag_subspace_env_create, general_diag_method_nr, krylov_space_create, & 00199 ot_diag_method_nr, ot_method_nr, qs_scf_env_type, scf_env_create, & 00200 scf_env_release, special_diag_method_nr 00201 USE qs_wf_history_methods, ONLY: wfi_extrapolate,& 00202 wfi_get_method_label,& 00203 wfi_update 00204 USE scf_control_types, ONLY: scf_control_type,& 00205 smear_type 00206 USE scp_coeff_types, ONLY: aux_coeff_set_type 00207 USE scp_density_methods, ONLY: update_rhoscp 00208 USE scp_energy_types, ONLY: scp_energy_type 00209 USE scp_environment_types, ONLY: get_scp_env,& 00210 scp_environment_type,& 00211 set_scp_env 00212 USE termination, ONLY: stop_program 00213 USE timings, ONLY: timeset,& 00214 timestop 00215 USE xas_env_types, ONLY: xas_environment_type 00216 USE xas_restart, ONLY: xas_initialize_rho 00217 #include "cp_common_uses.h" 00218 00219 IMPLICIT NONE 00220 00221 PRIVATE 00222 00223 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf' 00224 00225 PUBLIC :: scf, scf_env_cleanup, scf_env_do_scf,& 00226 init_scf_run, init_scf_loop, & 00227 qs_scf_print_summary 00228 00229 CONTAINS 00230 00231 ! ***************************************************************************** 00241 SUBROUTINE scf(qs_env,error) 00242 TYPE(qs_environment_type), POINTER :: qs_env 00243 TYPE(cp_error_type), INTENT(INOUT) :: error 00244 00245 CHARACTER(len=*), PARAMETER :: routineN = 'scf', 00246 routineP = moduleN//':'//routineN 00247 00248 INTEGER :: max_scf_tmp 00249 LOGICAL :: converged, failure, 00250 outer_scf_loop, should_stop 00251 LOGICAL, SAVE :: first_step_flag = .TRUE. 00252 TYPE(cp_logger_type), POINTER :: logger 00253 TYPE(dft_control_type), POINTER :: dft_control 00254 TYPE(qs_rho_type), POINTER :: rho 00255 TYPE(qs_scf_env_type), POINTER :: scf_env 00256 TYPE(scf_control_type), POINTER :: scf_control 00257 TYPE(section_vals_type), POINTER :: dft_section, input, 00258 mixing_section, scf_section 00259 00260 NULLIFY(scf_env) 00261 failure=.FALSE. 00262 logger => cp_error_get_logger(error) 00263 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 00264 IF (.NOT.failure) THEN 00265 CALL get_qs_env(qs_env,scf_env=scf_env,error=error,input=input, & 00266 dft_control=dft_control,scf_control=scf_control) 00267 00268 dft_section => section_vals_get_subs_vals(input,"DFT",error=error) 00269 scf_section => section_vals_get_subs_vals(dft_section,"SCF",error=error) 00270 00271 IF (.NOT.ASSOCIATED(scf_env)) THEN ! i.e. for MD this is associated on the second step (it so seems) 00272 CALL scf_env_create(scf_env, error=error) 00273 CALL set_qs_env(qs_env,scf_env=scf_env,error=error) 00274 CALL scf_env_release(scf_env,error=error) 00275 00276 ! set some of the methods that might be used in this SCF. 00277 ! this might not yet be the ideal place to set this kind of stuff (who knows?) 00278 CALL get_qs_env(qs_env,scf_env=scf_env,scf_control=scf_control, & 00279 dft_control=dft_control, error=error) 00280 00281 IF (scf_control%use_diag) THEN 00282 IF ( scf_control%diagonalization%method == diag_standard ) THEN 00283 scf_env%method=general_diag_method_nr 00284 IF (qs_env%has_unit_metric) THEN 00285 scf_env%method=special_diag_method_nr 00286 END IF 00287 ELSEIF ( scf_control%diagonalization%method == diag_ot ) THEN 00288 scf_env%method=ot_diag_method_nr 00289 IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%semi_empirical) THEN 00290 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00291 routineP,"DFTB and SE not possible with OT diagonalization",error,failure) 00292 END IF 00293 ELSEIF ( scf_control%diagonalization%method == diag_block_krylov ) THEN 00294 scf_env%method=block_krylov_diag_method_nr 00295 IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%semi_empirical) THEN 00296 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00297 routineP,"DFTB and SE not possible with block diagonalization",error,failure) 00298 END IF 00299 CALL krylov_space_create(scf_env%krylov_space, scf_section,error=error) 00300 ELSEIF ( scf_control%diagonalization%method == diag_block_davidson ) THEN 00301 scf_env%method=block_davidson_diag_method_nr 00302 IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%semi_empirical) THEN 00303 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00304 routineP,"DFTB and SE not possible with block diagonalization",error,failure) 00305 END IF 00306 CALL block_davidson_env_create(scf_env%block_davidson_env,dft_control%nspins,& 00307 scf_section,error=error) 00308 ELSE 00309 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00310 routineP,"Unknown diagonalization method",error,failure) 00311 END IF 00312 IF(scf_control%do_diag_sub) THEN 00313 scf_env%do_diag_sub = .TRUE. 00314 CALL diag_subspace_env_create(scf_env%subspace_env,scf_section,& 00315 dft_control%qs_control%cutoff,error=error) 00316 END IF 00317 ELSEIF (scf_control%use_ot) THEN 00318 scf_env%method=ot_method_nr 00319 ELSE 00320 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00321 routineP,"OT or DIAGONALIZATION have to be set",error,failure) 00322 END IF 00323 00324 SELECT CASE(scf_control%mixing_method) 00325 CASE(no_mix) 00326 scf_env%mixing_method=no_mixing_nr 00327 scf_env%p_mix_alpha = 1.0_dp 00328 CASE(direct_p_mix,kerker_mix,pulay_mix,broy_mix,broy_mix_new,multisec_mix) 00329 scf_env%mixing_method = scf_control%mixing_method 00330 mixing_section => section_vals_get_subs_vals(scf_section,"MIXING",error=error) 00331 CALL mixing_storage_create(scf_env%mixing_store, mixing_section, scf_env%mixing_method, & 00332 dft_control%qs_control%cutoff, error=error) 00333 CASE DEFAULT 00334 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00335 routineP,"Unknown mixing method",error,failure) 00336 END SELECT 00337 00338 ! Disable DIIS for OT and g-space density mixing methods 00339 IF( scf_env%method==ot_method_nr ) THEN 00340 ! No mixing is used with OT 00341 scf_env%mixing_method=no_mixing_nr 00342 scf_env%p_mix_alpha = 1.0_dp 00343 scf_env%skip_diis = .TRUE. 00344 END IF 00345 00346 IF (scf_control%use_diag .AND. scf_env%mixing_method==no_mixing_nr) THEN 00347 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00348 routineP,"Diagonalization procedures without mixing are not recommendable",error,failure) 00349 END IF 00350 00351 IF(scf_env%mixing_method>direct_mixing_nr) THEN 00352 scf_env%skip_diis = .TRUE. 00353 scf_env%p_mix_alpha = scf_env%mixing_store%alpha 00354 IF(scf_env%mixing_store%beta==0.0_dp) THEN 00355 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00356 routineP,"Mixing employing the Kerker damping factor needs BETA /= 0.0",error,failure) 00357 END IF 00358 END IF 00359 00360 IF (scf_env%mixing_method==direct_mixing_nr) THEN 00361 scf_env%p_mix_alpha = scf_env%mixing_store%alpha 00362 IF(qs_env%scf_control%eps_diis<qs_env%scf_control%eps_scf) THEN 00363 scf_env%skip_diis = .TRUE. 00364 CALL cp_assert(.FALSE.,cp_warning_level,cp_assertion_failed,routineP,& 00365 "the DIIS scheme is disabled, since EPS_DIIS < EPS_SCF",& 00366 error,failure) 00367 END IF 00368 END IF 00369 00370 CALL section_vals_val_get(scf_section,"CHOLESKY",i_val=scf_env%cholesky_method,error=error) 00371 ELSE 00372 ! Reallocate mixing store, if the g space grid (cell) has changed 00373 SELECT CASE (scf_env%mixing_method) 00374 CASE (kerker_mix,pulay_mix,broy_mix,broy_mix_new,multisec_mix) 00375 IF (ASSOCIATED(scf_env%mixing_store)) THEN 00376 CALL get_qs_env(qs_env,rho=rho,error=error) 00377 ! The current mixing_store data structure does not allow for an unique 00378 ! grid comparison, but the probability that the 1d lengths of the old and 00379 ! the new grid are accidentily equal is rather low 00380 IF (SIZE(rho%rho_g(1)%pw%pw_grid%gsq) /= SIZE(scf_env%mixing_store%rhoin(1)%cc)) THEN 00381 CALL mixing_storage_release(scf_env%mixing_store,error=error) 00382 mixing_section => section_vals_get_subs_vals(scf_section,"MIXING",error=error) 00383 CALL mixing_storage_create(scf_env%mixing_store,mixing_section,scf_env%mixing_method,& 00384 dft_control%qs_control%cutoff,error=error) 00385 END IF 00386 END IF 00387 END SELECT 00388 END IF 00389 00390 IF (qs_env%scf_control%max_scf > 0) THEN 00391 00392 CALL init_scf_run(scf_env=scf_env, qs_env=qs_env, scf_section=scf_section, error=error) 00393 00394 IF ((qs_env%scf_control%density_guess .EQ. history_guess) .AND. (first_step_flag)) THEN 00395 max_scf_tmp = qs_env%scf_control%max_scf 00396 qs_env%scf_control%max_scf = 1 00397 outer_scf_loop = qs_env%scf_control%outer_scf%have_scf 00398 qs_env%scf_control%outer_scf%have_scf = .FALSE. 00399 END IF 00400 CALL scf_env_do_scf(scf_env=scf_env, qs_env=qs_env, & 00401 converged=converged, should_stop=should_stop, error=error) 00402 00403 ! *** add the converged wavefunction to the wavefunction history 00404 IF ((ASSOCIATED(qs_env%wf_history)) .AND. & 00405 ((qs_env%scf_control%density_guess .NE. history_guess) .OR. & 00406 (.NOT. first_step_flag))) THEN 00407 CALL wfi_update(qs_env%wf_history,qs_env=qs_env,dt=1.0_dp, error=error) 00408 ELSE IF ((qs_env%scf_control%density_guess .EQ. history_guess) .AND. & 00409 (first_step_flag)) THEN 00410 qs_env%scf_control%max_scf = max_scf_tmp 00411 qs_env%scf_control%outer_scf%have_scf = outer_scf_loop 00412 first_step_flag = .FALSE. 00413 END IF 00414 00415 ! *** compute properties that depend on the converged wavefunction 00416 IF (.NOT.(should_stop)) THEN 00417 IF (dft_control%qs_control%semi_empirical) THEN 00418 CALL scf_post_calculation_se (dft_section, scf_env, qs_env, error) 00419 ELSEIF (dft_control%qs_control%dftb) THEN 00420 CALL scf_post_calculation_dftb (dft_section, scf_env, qs_env, error) 00421 ELSEIF (dft_control%qs_control%scptb) THEN 00422 ! Do Nothing 00423 ELSEIF (dft_control%qs_control%ofgpw) THEN 00424 ! Do Nothing 00425 ELSE 00426 CALL scf_post_calculation_gpw(dft_section, scf_env, qs_env, error) 00427 END IF 00428 00429 ! Compute MP2 energy 00430 IF (ASSOCIATED(qs_env%mp2_env)) THEN 00431 CALL mp2_main(qs_env=qs_env,error=error) 00432 ENDIF 00433 00434 END IF 00435 00436 END IF 00437 00438 00439 ! *** cleanup 00440 CALL scf_env_cleanup(scf_env,qs_env=qs_env,error=error) 00441 00442 END IF 00443 00444 END SUBROUTINE scf 00445 00446 ! ***************************************************************************** 00461 SUBROUTINE scf_env_do_scf(scf_env,qs_env,converged,should_stop,error) 00462 00463 TYPE(qs_scf_env_type), POINTER :: scf_env 00464 TYPE(qs_environment_type), POINTER :: qs_env 00465 LOGICAL, INTENT(OUT) :: converged, should_stop 00466 TYPE(cp_error_type), INTENT(inout) :: error 00467 00468 CHARACTER(LEN=*), PARAMETER :: routineN = 'scf_env_do_scf', 00469 routineP = moduleN//':'//routineN 00470 00471 INTEGER :: handle, handle2, ispin, 00472 iter_count, output_unit, 00473 total_steps 00474 LOGICAL :: diis_step, do_std_diag, energy_only, exit_inner_loop, failure, 00475 gapw, gapw_xc, harris_energy, harris_flag, has_unit_metric, 00476 just_energy, outer_loop_converged, scp_dft, scp_nddo, skip_diag_sub, 00477 use_jacobi 00478 REAL(KIND=dp) :: outer_loop_eps, t1, t2 00479 TYPE(atomic_kind_type), DIMENSION(:), 00480 POINTER :: atomic_kind_set 00481 TYPE(aux_coeff_set_type), POINTER :: aux_coeff_set 00482 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00483 POINTER :: matrix_ks, matrix_s 00484 TYPE(cp_dbcsr_type), POINTER :: ks_scp, pscp 00485 TYPE(cp_logger_type), POINTER :: logger 00486 TYPE(cp_para_env_type), POINTER :: para_env 00487 TYPE(dft_control_type), POINTER :: dft_control 00488 TYPE(harris_env_type), POINTER :: harris_env 00489 TYPE(mo_set_p_type), DIMENSION(:), 00490 POINTER :: mos 00491 TYPE(particle_type), DIMENSION(:), 00492 POINTER :: particle_set 00493 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00494 TYPE(qs_charges_type), POINTER :: qs_charges 00495 TYPE(qs_energy_type), POINTER :: energy 00496 TYPE(qs_ks_env_type), POINTER :: ks_env 00497 TYPE(qs_rho_type), POINTER :: rho, rho_xc 00498 TYPE(scf_control_type), POINTER :: scf_control 00499 TYPE(scp_environment_type), POINTER :: scp_env 00500 TYPE(section_vals_type), POINTER :: dft_section, harris_section, 00501 input, scf_section 00502 00503 CALL timeset(routineN,handle) 00504 00505 failure=.FALSE. 00506 converged=.TRUE. 00507 00508 NULLIFY(dft_control,matrix_s,matrix_ks,rho,rho_xc,energy, & 00509 scf_control,logger,qs_charges,ks_env,mos,atomic_kind_set, & 00510 particle_set,harris_env,dft_section,input,& 00511 scf_section, para_env) 00512 NULLIFY ( aux_coeff_set, scp_env, pscp, ks_scp) 00513 00514 CPPrecondition(ASSOCIATED(scf_env),cp_failure_level,routineP,error,failure) 00515 CPPrecondition(scf_env%ref_count>0,cp_failure_level,routineP,error,failure) 00516 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 00517 CPPrecondition(qs_env%ref_count>0,cp_failure_level,routineP,error,failure) 00518 para_env=>qs_env%para_env 00519 00520 logger => cp_error_get_logger(error) 00521 t1 = m_walltime() 00522 00523 CALL get_qs_env(qs_env=qs_env,& 00524 matrix_s=matrix_s,energy=energy,& 00525 particle_set=particle_set,& 00526 qs_charges=qs_charges,& 00527 ks_env=ks_env, & 00528 harris_env=harris_env,& 00529 atomic_kind_set=atomic_kind_set,& 00530 matrix_ks=matrix_ks,rho=rho,rho_xc=rho_xc,mos=mos, & 00531 input=input, dft_control=dft_control, scf_control=scf_control, & 00532 error=error) 00533 00534 ! Defining SCP logicals 00535 scp_dft = dft_control%scp 00536 scp_nddo = dft_control%qs_control%se_control%scp 00537 00538 dft_section => section_vals_get_subs_vals(input,"DFT",error=error) 00539 scf_section => section_vals_get_subs_vals(dft_section,"SCF",error=error) 00540 harris_section => section_vals_get_subs_vals(dft_section,"QS%HARRIS", error=error) 00541 00542 output_unit=cp_print_key_unit_nr(logger,scf_section,"PRINT%PROGRAM_RUN_INFO",& 00543 extension=".scfLog",error=error) 00544 00545 IF (output_unit>0) WRITE (UNIT=output_unit,FMT="(/,/,T2,A)") & 00546 "SCF WAVEFUNCTION OPTIMIZATION" 00547 00548 ! short cut flags setting the different methods 00549 gapw = dft_control%qs_control%gapw 00550 gapw_xc = dft_control%qs_control%gapw_xc 00551 ! **** SCP 00552 IF ( ( scp_dft ) .OR. ( scp_nddo ) ) THEN 00553 CALL get_qs_env ( qs_env = qs_env, scp_env = scp_env, error=error ) 00554 END IF 00555 ! Assign SCP-DFT structures 00556 IF ( scp_dft ) THEN 00557 CALL get_scp_env ( scp_env = scp_env, aux_coeff_set = aux_coeff_set, error = error ) 00558 ENDIF 00559 ! Assign SCP-NDDO structures 00560 IF ( scp_nddo ) THEN 00561 CALL get_scp_env ( scp_env = scp_env, pscp=pscp, ks_scp=ks_scp, error = error ) 00562 ENDIF 00563 ! **** SCP 00564 harris_flag = qs_env%use_harris 00565 CALL get_qs_env(qs_env,has_unit_metric=has_unit_metric,error=error) 00566 00567 IF ((output_unit > 0).AND.(.NOT.scf_control%use_ot)) THEN 00568 WRITE (UNIT=output_unit,& 00569 FMT="(/,T3,A,T12,A,T31,A,T39,A,T59,A,T75,A,/,T3,A)")& 00570 "Step","Update method","Time","Convergence","Total energy","Change",& 00571 REPEAT("-",78) 00572 END IF 00573 CALL cp_add_iter_level(logger%iter_info,"QS_SCF",error=error) 00574 ! *** outer loop of the scf, can treat other variables, 00575 ! *** such as lagrangian multipliers 00576 scf_env%outer_scf%iter_count=0 00577 iter_count = 0 00578 total_steps = 0 00579 energy%tot_old = 0.0_dp 00580 scf_outer_loop: DO 00581 00582 CALL init_scf_loop(scf_env=scf_env, qs_env=qs_env, & 00583 scf_section=scf_section, error=error) 00584 00585 ! ****** Setting Some SCP flags******* 00586 ! check if we are doing SCP: 00587 ! **** qs_ot_env%settings%scp_dft,qs_ot_env%settings%ks ARE BOTH .TRUE. here ***** 00588 IF ( scp_dft ) THEN 00589 ! check if we are doing outer SCF optimization of SCP coeffs: 00590 IF ( scf_control%outer_scf%type==outer_scf_none ) THEN 00591 scf_env%qs_ot_env(1)%settings%ks = .TRUE. 00592 !*** qs_ot_env%settings%scp_dft*** may be set .FALSE. on input 00593 ! and we will preserve this if there is no outer_scf 00594 ELSEIF ( scf_control%outer_scf%type==outer_scf_scp ) THEN 00595 ! IF outer scf count is even, then only optimize SCP, else optimize KS 00596 IF ( scf_env%outer_scf%iter_count == 0 ) THEN 00597 scf_env%qs_ot_env(1)%settings%ks = .TRUE. 00598 scf_env%qs_ot_env(1)%settings%scp_dft = .FALSE. 00599 ELSEIF ( MOD ( scf_env%outer_scf%iter_count, 2 ) /= 0 ) THEN 00600 scf_env%qs_ot_env(1)%settings%ks = .FALSE. 00601 scf_env%qs_ot_env(1)%settings%scp_dft = .TRUE. 00602 ELSE 00603 scf_env%qs_ot_env(1)%settings%ks = .TRUE. 00604 scf_env%qs_ot_env(1)%settings%scp_dft = .FALSE. 00605 END IF 00606 ENDIF 00607 END IF 00608 IF ( scp_nddo ) THEN 00609 ! check if we are doing outer SCF optimization of SCP coeffs: 00610 IF ( scf_control%outer_scf%type==outer_scf_none ) THEN 00611 scf_env%qs_ot_env(1)%settings%ks = .TRUE. 00612 !*** qs_ot_env%settings%scp_nddo*** may be set .FALSE. on input 00613 ! and we will preserve this if there is no outer_scf 00614 ELSEIF ( scf_control%outer_scf%type==outer_scf_scp ) THEN 00615 ! IF outer scf count is even, then only optimize SCP, else optimize KS 00616 IF ( scf_env%outer_scf%iter_count == 0 ) THEN 00617 scf_env%qs_ot_env(1)%settings%ks = .TRUE. 00618 scf_env%qs_ot_env(1)%settings%scp_nddo = .FALSE. 00619 ELSEIF ( MOD ( scf_env%outer_scf%iter_count, 2 ) /= 0 ) THEN 00620 scf_env%qs_ot_env(1)%settings%ks = .FALSE. 00621 scf_env%qs_ot_env(1)%settings%scp_nddo = .TRUE. 00622 ELSE 00623 scf_env%qs_ot_env(1)%settings%ks = .TRUE. 00624 scf_env%qs_ot_env(1)%settings%scp_nddo = .FALSE. 00625 END IF 00626 ENDIF 00627 END IF 00628 00629 00630 ! Some flags needed to be set at the beginning of the loop 00631 00632 diis_step = .FALSE. 00633 use_jacobi = .FALSE. 00634 energy_only = .FALSE. 00635 just_energy = .FALSE. 00636 00637 ! SCF loop, optimisation of the wfn coefficients 00638 ! qs_env%rho%rho_r and qs_env%rho%rho_g should be up to date here 00639 00640 scf_env%iter_count = 0 00641 exit_inner_loop = .FALSE. 00642 scf_loop: DO 00643 00644 CALL timeset(routineN//"_inner_loop",handle2) 00645 00646 scf_env%iter_count = scf_env%iter_count + 1 00647 iter_count = iter_count + 1 00648 CALL cp_iterate(logger%iter_info,last=.FALSE.,iter_nr=iter_count,error=error) 00649 00650 IF (output_unit > 0) CALL m_flush(output_unit) 00651 00652 total_steps = total_steps + 1 00653 just_energy = energy_only 00654 ! ***SCP 00655 IF ( scp_dft ) CALL qs_ks_scp_update ( qs_env, just_energy=just_energy, error=error ) 00656 ! ***SCP 00657 CALL qs_ks_update_qs_env(ks_env,qs_env=qs_env,calculate_forces=.FALSE.,& 00658 just_energy=just_energy,error=error) 00659 00660 ! print 'heavy weight' or relatively expensive quantities 00661 CALL qs_scf_loop_print(qs_env,scf_env,para_env,error) 00662 00663 scf_env%iter_param = 0.0_dp 00664 00665 ! this takes known energy and derivatives and produces 00666 ! new wfns and or density matrix 00667 SELECT CASE (scf_env%method) 00668 CASE DEFAULT 00669 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00670 routineP,"unknown scf method: "//& 00671 cp_to_string(scf_env%method),error,failure) 00672 CASE(general_diag_method_nr) 00673 IF (dft_control%roks) THEN 00674 CALL do_roks_diag(scf_env,mos,matrix_ks,matrix_s,& 00675 scf_control,scf_section,diis_step,& 00676 has_unit_metric,error) 00677 ELSE 00678 CALL do_general_diag(scf_env,mos,matrix_ks,& 00679 matrix_s,scf_control,scf_section, & 00680 diis_step,use_jacobi,error=error) 00681 IF(scf_env%do_diag_sub) THEN 00682 skip_diag_sub = (scf_env%subspace_env%eps_diag_sub > 0.0_dp) .AND. & 00683 (scf_env%iter_count==1 .OR. scf_env%iter_delta>scf_env%subspace_env%eps_diag_sub ) 00684 IF( .NOT. skip_diag_sub) THEN 00685 CALL do_scf_diag_subspace(qs_env,scf_env,scf_env%subspace_env,mos,rho,& 00686 ks_env,scf_section,scf_control,error=error) 00687 END IF 00688 END IF 00689 END IF 00690 CASE(special_diag_method_nr) 00691 IF (dft_control%roks) THEN 00692 CALL do_roks_diag(scf_env,mos,matrix_ks,matrix_s,& 00693 scf_control,scf_section,diis_step,& 00694 has_unit_metric,error) 00695 ELSE 00696 CALL do_special_diag(scf_env,mos,matrix_ks,& 00697 scf_control,scf_section, & 00698 diis_step,error) 00699 END IF 00700 CASE(ot_diag_method_nr) 00701 CALL cp_assert(.NOT.(dft_control%roks),cp_failure_level,cp_assertion_failed,& 00702 routineP,"ROKS with OT diagonalization not possible",error,failure) 00703 CALL do_ot_diag(scf_env,mos,matrix_ks,matrix_s,& 00704 scf_control,scf_section,diis_step,error) 00705 CASE(block_krylov_diag_method_nr) 00706 CALL cp_assert(.NOT.(dft_control%roks),cp_failure_level,cp_assertion_failed,& 00707 routineP,"ROKS with block PF diagonalization not possible",error,failure) 00708 do_std_diag = (scf_env%krylov_space%eps_std_diag > 0.0_dp) 00709 IF(do_std_diag .AND.& 00710 (scf_env%iter_count==1 .OR. scf_env%iter_delta>scf_env%krylov_space%eps_std_diag)) THEN 00711 IF(scf_env%krylov_space%always_check_conv) THEN 00712 CALL do_block_krylov_diag(scf_env,mos,matrix_ks,& 00713 scf_control, scf_section, check_moconv_only=.TRUE., error=error) 00714 END IF 00715 CALL do_general_diag(scf_env,mos,matrix_ks,& 00716 matrix_s,scf_control,scf_section, diis_step,use_jacobi,error=error) 00717 ELSE 00718 CALL do_block_krylov_diag(scf_env,mos,matrix_ks, & 00719 scf_control, scf_section, error=error) 00720 END IF 00721 IF(scf_env%do_diag_sub) THEN 00722 skip_diag_sub = (scf_env%subspace_env%eps_diag_sub > 0.0_dp) .AND. & 00723 (scf_env%iter_count==1 .OR. scf_env%iter_delta>scf_env%subspace_env%eps_diag_sub ) 00724 IF( .NOT. skip_diag_sub) THEN 00725 CALL do_scf_diag_subspace(qs_env,scf_env,scf_env%subspace_env,mos,rho,& 00726 ks_env,scf_section, scf_control,error=error) 00727 END IF 00728 END IF 00729 CASE(block_davidson_diag_method_nr) 00730 00731 CALL do_block_davidson_diag(qs_env,scf_env,mos,matrix_ks,matrix_s,scf_control,& 00732 scf_section,.FALSE.,error=error) 00733 00734 CASE(ot_method_nr) 00735 ! ***SCP 00736 IF ( scp_dft ) THEN 00737 CALL qs_scf_loop_do_ot(qs_env,scf_env,scf_control%smear,mos,rho,& 00738 qs_env%mo_derivs,energy%total, & 00739 matrix_s, aux_coeff_set=aux_coeff_set, energy_only=energy_only, & 00740 has_unit_metric=has_unit_metric,error=error) 00741 ELSEIF ( scp_nddo ) THEN 00742 CALL get_scp_env ( scp_env = scp_env, ks_scp=ks_scp, error = error ) 00743 CALL qs_scf_loop_do_ot(qs_env,scf_env,scf_control%smear,mos,rho,& 00744 qs_env%mo_derivs,energy%total, & 00745 matrix_s, pscp=pscp, fscp=ks_scp, energy_only=energy_only, & 00746 has_unit_metric=has_unit_metric,error=error) 00747 CALL set_scp_env ( scp_env = scp_env, pscp = pscp, error = error ) 00748 ELSE 00749 CALL qs_scf_loop_do_ot(qs_env,scf_env,scf_control%smear,mos,rho,& 00750 qs_env%mo_derivs,energy%total, & 00751 matrix_s, energy_only=energy_only,has_unit_metric=has_unit_metric,error=error) 00752 ENDIF 00753 ! ***SCP 00754 END SELECT 00755 00756 ! another heavy weight print object, print controlled by dft_section 00757 IF (SIZE(mos) > 1) THEN 00758 CALL write_mo_set(mos(1)%mo_set,atomic_kind_set,particle_set,4,& 00759 dft_section,spin="ALPHA",last=.FALSE.,error=error) 00760 CALL write_mo_set(mos(2)%mo_set,atomic_kind_set,particle_set,4,& 00761 dft_section,spin="BETA",last=.FALSE.,error=error) 00762 ELSE 00763 CALL write_mo_set(mos(1)%mo_set,atomic_kind_set,particle_set,4,& 00764 dft_section,last=.FALSE.,error=error) 00765 END IF 00766 00767 ! ** calculation of the harris energy correction *** 00768 IF (harris_flag) THEN 00769 CALL harris_energy_correction(qs_env, harris_env, para_env=para_env, fast=.TRUE., error=error) 00770 IF (scf_env%method .NE. ot_method_nr) THEN 00771 CALL harris_eigenvalue_trace_KS_Pmix(scf_env, qs_env, harris_env,error=error) 00772 ELSE 00773 CALL harris_eigenvalue_calculation(qs_env=qs_env, harris_env=harris_env,error=error) 00774 END IF 00775 00776 CALL section_vals_val_get(harris_section, "HARRIS_ENERGY",l_val=harris_energy,error=error) 00777 IF ((harris_energy)) THEN 00778 energy%total = harris_env%harris_energy%Eharris 00779 END IF 00780 END IF 00781 00782 SELECT CASE(scf_env%mixing_method) 00783 CASE(direct_mixing_nr) 00784 CALL scf_env_density_mixing(scf_env%p_mix_new,& 00785 scf_env%mixing_store, rho%rho_ao, para_env, scf_env%iter_delta, scf_env%iter_count, & 00786 diis=diis_step, error=error) 00787 CASE(gspace_mixing_nr,pulay_mixing_nr,broyden_mixing_nr,& 00788 broyden_mixing_new_nr,multisecant_mixing_nr) 00789 ! Compute the difference p_out-p_in 00790 CALL self_consistency_check(qs_env%rho%rho_ao,scf_env%p_delta,para_env,scf_env%p_mix_new,& 00791 delta= scf_env%iter_delta, error=error) 00792 CASE(no_mixing_nr) 00793 CASE DEFAULT 00794 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 00795 routineP,"unknown scf mixing method: "//& 00796 cp_to_string(scf_env%mixing_method),error,failure) 00797 END SELECT 00798 00799 t2 = m_walltime() 00800 IF ((output_unit > 0).AND.scf_env%print_iter_line) THEN 00801 IF (just_energy) THEN 00802 WRITE (UNIT=output_unit,& 00803 FMT="(T2,I5,1X,A,T20,E8.2,1X,F6.1,16X,F20.10)")& 00804 scf_env%iter_count,TRIM(scf_env%iter_method),scf_env%iter_param,& 00805 t2 - t1,energy%total 00806 ELSE 00807 IF ((ABS(scf_env%iter_delta) < 1.0E-8_dp).OR.& 00808 (ABS(scf_env%iter_delta) >= 1.0E5_dp)) THEN 00809 WRITE (UNIT=output_unit,& 00810 FMT="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,ES14.4,1X,F20.10,1X,ES9.2)")& 00811 scf_env%iter_count,TRIM(scf_env%iter_method),scf_env%iter_param,& 00812 t2 - t1,scf_env%iter_delta,energy%total,energy%total-energy%tot_old 00813 ELSE 00814 WRITE (UNIT=output_unit,& 00815 FMT="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,F14.8,1X,F20.10,1X,ES9.2)")& 00816 scf_env%iter_count,TRIM(scf_env%iter_method),scf_env%iter_param,& 00817 t2 - t1,scf_env%iter_delta,energy%total,energy%total-energy%tot_old 00818 END IF 00819 END IF 00820 END IF 00821 IF (.NOT.just_energy) energy%tot_old = energy%total 00822 IF (harris_flag) THEN 00823 CALL harris_postprocessing(harris_env, error=error) 00824 END IF 00825 00826 ! ** can we exit the SCF loop ? 00827 CALL external_control(should_stop,"SCF",target_time=qs_env%target_time, & 00828 start_time=qs_env%start_time,error=error) 00829 IF (scf_env%iter_delta < scf_control%eps_scf) THEN 00830 IF (output_unit>0) THEN 00831 WRITE(UNIT=output_unit,FMT="(/,T3,A,I5,A/)")& 00832 "*** SCF run converged in ",scf_env%iter_count," steps ***" 00833 END IF 00834 exit_inner_loop = .TRUE. 00835 ELSE IF (should_stop.OR. scf_env%iter_count >= scf_control%max_scf) THEN 00836 IF (output_unit>0) THEN 00837 WRITE(UNIT=output_unit,FMT="(/,T3,A,/)")& 00838 "*** SCF run NOT converged ***" 00839 END IF 00840 converged=.FALSE. 00841 exit_inner_loop = .TRUE. 00842 END IF 00843 ! ** In case we decide to exit we perform few more check to see if this one 00844 ! ** is really the last SCF step 00845 IF (exit_inner_loop) THEN 00846 ! ********* Need to reset qs_ot_env%settings%scp_dft = .TRUE. to ensure ********* 00847 ! consistent deallocation in cleanup_scf_loop 00848 IF ( scp_dft .AND. scf_control%outer_scf%type==outer_scf_scp ) & 00849 scf_env%qs_ot_env(1)%settings%scp_dft=.TRUE. 00850 ! ********* Need to reset qs_ot_env%settings%scp_dft = .TRUE. to ensure ********* 00851 ! consistent deallocation in cleanup_scf_loop 00852 IF ( scp_nddo .AND. scf_control%outer_scf%type==outer_scf_scp ) & 00853 scf_env%qs_ot_env(1)%settings%scp_nddo=.TRUE. 00854 00855 CALL cleanup_scf_loop(scf_env,error) 00856 00857 ! now, print out energies and charges corresponding to the obtained wfn 00858 ! (this actually is not 100% consistent at this point)! 00859 IF ( scp_dft ) CALL qs_scf_scp_print_summary ( output_unit, scp_env, qs_env%qmmm, error ) 00860 CALL qs_scf_print_summary(output_unit,rho,qs_charges,energy,scf_env%nelectron, & 00861 dft_control,qs_env%qmmm,qs_env,gapw,gapw_xc,error) 00862 00863 IF (harris_flag) THEN 00864 energy%total = harris_env%harris_energy%Eharris 00865 END IF 00866 00867 ! *** mixing methods need to undo mixing of the density matrix 00868 ! (restore original density) *** 00869 IF (scf_env%mixing_method>0) THEN 00870 SELECT CASE(scf_env%mixing_method) 00871 CASE(direct_mixing_nr) 00872 CALL scf_env_density_mixing(scf_env%p_mix_new,scf_env%mixing_store,& 00873 rho%rho_ao,para_env,scf_env%iter_delta,& 00874 scf_env%iter_count,diis=diis_step,& 00875 invert=.TRUE.,error=error) 00876 DO ispin=1,dft_control%nspins 00877 CALL cp_dbcsr_copy(rho%rho_ao(ispin)%matrix,scf_env%p_mix_new(ispin)%matrix,& 00878 name=TRIM(cp_dbcsr_name(rho%rho_ao(ispin)%matrix)),error=error) 00879 END DO 00880 CASE(gspace_mixing_nr,pulay_mixing_nr,broyden_mixing_nr,& 00881 broyden_mixing_new_nr,multisecant_mixing_nr) 00882 DO ispin=1,dft_control%nspins 00883 CALL cp_dbcsr_copy(rho%rho_ao(ispin)%matrix,scf_env%p_mix_new(ispin)%matrix,& 00884 name=TRIM(cp_dbcsr_name(rho%rho_ao(ispin)%matrix)),error=error) 00885 END DO 00886 END SELECT 00887 ENDIF 00888 00889 ! *** update rspace rho since the mo changed 00890 ! *** this might not always be needed (i.e. no post calculation / no forces ) 00891 ! *** but guarantees that rho and wfn are consistent at this point 00892 CALL qs_rho_update_rho(rho, qs_env=qs_env, error=error) 00893 CALL qs_ks_did_change(ks_env,rho_changed=.TRUE.,error=error) 00894 ! ***SCP 00895 IF ( scp_dft ) THEN 00896 CALL update_rhoscp ( qs_env = qs_env, error = error ) 00897 CALL qs_ks_scp_did_change(qs_env, rho_changed=.TRUE.,error=error) 00898 ENDIF 00899 ! ***SCP 00900 00901 outer_loop_converged=.TRUE. 00902 IF (scf_control%outer_scf%have_scf) THEN 00903 ! We have an outer SCF loop... 00904 scf_env%outer_scf%iter_count=scf_env%outer_scf%iter_count+1 00905 outer_loop_converged=.FALSE. 00906 00907 CALL outer_loop_gradient(qs_env,scf_env,error) 00908 outer_loop_eps=SQRT(SUM(scf_env%outer_scf%gradient(:,scf_env%outer_scf%iter_count)**2))/ & 00909 SIZE(scf_env%outer_scf%gradient,1) 00910 IF (outer_loop_eps<scf_control%outer_scf%eps_scf) outer_loop_converged=.TRUE. 00911 END IF 00912 ! ** Let's tag the last SCF cycle so we can print informations only of the last step 00913 IF (outer_loop_converged.OR.& 00914 scf_env%outer_scf%iter_count>scf_control%outer_scf%max_scf .OR.& 00915 should_stop) THEN 00916 CALL cp_iterate(logger%iter_info,last=.TRUE.,iter_nr=iter_count,error=error) 00917 END IF 00918 END IF 00919 00920 ! *** Write WaveFunction restart file *** 00921 CALL write_mo_set(mos,particle_set,dft_section=dft_section,scp=scp_dft,& 00922 scp_env=scp_env,& 00923 atomic_kind_set=atomic_kind_set,error=error) 00924 00925 ! *** Exit if we have finished with the SCF inner loop *** 00926 IF (exit_inner_loop) THEN 00927 CALL timestop(handle2) 00928 EXIT scf_loop 00929 END IF 00930 00931 IF (.NOT.BTEST(cp_print_key_should_output(logger%iter_info,& 00932 scf_section,"PRINT%ITERATION_INFO/TIME_CUMUL",error=error),cp_p_file)) & 00933 t1 = m_walltime() 00934 00935 ! ** mixing methods have the new density matrix in p_mix_new ** 00936 IF (scf_env%mixing_method > 0) THEN 00937 DO ispin=1,dft_control%nspins 00938 CALL cp_dbcsr_copy(rho%rho_ao(ispin)%matrix,scf_env%p_mix_new(ispin)%matrix,& 00939 name=TRIM(cp_dbcsr_name(rho%rho_ao(ispin)%matrix)),error=error) 00940 END DO 00941 END IF 00942 00943 ! ** update qs_env%rho 00944 CALL qs_rho_update_rho(rho, qs_env=qs_env, error=error) 00945 ! ** Density mixing through density matrix or on the reciprocal space grid (exclusive) 00946 IF(scf_env%mixing_method>=gspace_mixing_nr) THEN 00947 CALL gspace_mixing(qs_env, scf_env, scf_env%mixing_store, rho, qs_env%para_env, error=error) 00948 END IF 00949 CALL qs_ks_did_change(ks_env,rho_changed=.TRUE.,error=error) 00950 00951 ! *** SCP 00952 IF ( scp_dft ) THEN 00953 CALL update_rhoscp ( qs_env = qs_env, error = error ) 00954 CALL qs_ks_scp_did_change(qs_env, rho_changed=.TRUE.,error=error) 00955 ENDIF 00956 ! *** SCP 00957 CALL timestop(handle2) 00958 00959 END DO scf_loop 00960 00961 IF (.NOT. scf_control%outer_scf%have_scf) EXIT scf_outer_loop 00962 00963 ! ** In case we use the OUTER SCF loop let's print some info.. 00964 IF (output_unit>0) WRITE(output_unit,'(/,T3,A,I4,A,E10.2,A,F22.10)') & 00965 "outer SCF iter = ",scf_env%outer_scf%iter_count, & 00966 " RMS gradient = ",outer_loop_eps," energy =",energy%total 00967 00968 IF (outer_loop_converged) THEN 00969 IF (output_unit>0) WRITE(output_unit,'(T3,A,I4,A,I4,A,/)') & 00970 "outer SCF loop converged in",scf_env%outer_scf%iter_count,& 00971 " iterations or ",total_steps," steps" 00972 EXIT scf_outer_loop 00973 END IF 00974 00975 IF (scf_env%outer_scf%iter_count>scf_control%outer_scf%max_scf & 00976 .OR. should_stop ) THEN 00977 IF (output_unit>0) WRITE(output_unit,'(T3,A,I4,A,I4,A,/)') & 00978 "outer SCF loop FAILED to converge after ", & 00979 scf_env%outer_scf%iter_count," iterations or ",total_steps," steps" 00980 EXIT scf_outer_loop 00981 END IF 00982 00983 CALL outer_loop_optimize(scf_env,scf_control,error) 00984 CALL outer_loop_update_qs_env(qs_env,scf_env,error) 00985 CALL qs_ks_did_change(ks_env,potential_changed=.TRUE.,error=error) 00986 ! *** SCP 00987 IF ( scp_dft ) CALL qs_ks_scp_did_change(qs_env, potential_changed=.TRUE.,error=error) 00988 ! *** SCP 00989 00990 END DO scf_outer_loop 00991 00992 00993 !if needed copy mo_coeff dbcsr->fm for later use in post_scf!fm->dbcsr 00994 DO ispin=1,SIZE(mos)!fm -> dbcsr 00995 IF(mos(ispin)%mo_set%use_mo_coeff_b) THEN!fm->dbcsr 00996 IF (.NOT.ASSOCIATED(mos(ispin)%mo_set%mo_coeff_b)) &!fm->dbcsr 00997 CALL stop_program(routineN,moduleN,__LINE__,"mo_coeff_b is not allocated")!fm->dbcsr 00998 CALL copy_dbcsr_to_fm(mos(ispin)%mo_set%mo_coeff_b, &!fm->dbcsr 00999 mos(ispin)%mo_set%mo_coeff,error=error)!fm -> dbcsr 01000 ENDIF!fm->dbcsr 01001 ENDDO!fm -> dbcsr 01002 01003 01004 IF(qs_env%dft_control%qs_control%becke_restraint)THEN 01005 CALL pw_env_get(qs_env%pw_env,auxbas_pw_pool=auxbas_pw_pool,error=error) 01006 CALL pw_pool_give_back_pw(auxbas_pw_pool,& 01007 qs_env%dft_control%qs_control%becke_control%becke_pot%pw,error=error) 01008 qs_env%dft_control%qs_control%becke_control%need_pot=.TRUE. 01009 END IF 01010 CALL cp_rm_iter_level(logger%iter_info,level_name="QS_SCF",error=error) 01011 CALL timestop(handle) 01012 01013 END SUBROUTINE scf_env_do_scf 01014 01015 ! ***************************************************************************** 01021 SUBROUTINE qs_scf_loop_do_ot(qs_env,scf_env,smear,mos,rho,mo_derivs,total_energy,& 01022 matrix_s,aux_coeff_set,pscp,fscp,energy_only, & 01023 has_unit_metric,error) 01024 01025 TYPE(qs_environment_type), POINTER :: qs_env 01026 TYPE(qs_scf_env_type), POINTER :: scf_env 01027 TYPE(smear_type), POINTER :: smear 01028 TYPE(mo_set_p_type), DIMENSION(:), 01029 POINTER :: mos 01030 TYPE(qs_rho_type), POINTER :: rho 01031 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01032 POINTER :: mo_derivs 01033 REAL(KIND=dp), INTENT(IN) :: total_energy 01034 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01035 POINTER :: matrix_s 01036 TYPE(aux_coeff_set_type), OPTIONAL, 01037 POINTER :: aux_coeff_set 01038 TYPE(cp_dbcsr_type), OPTIONAL, POINTER :: pscp, fscp 01039 LOGICAL, INTENT(INOUT) :: energy_only 01040 LOGICAL, INTENT(IN) :: has_unit_metric 01041 TYPE(cp_error_type), INTENT(inout) :: error 01042 01043 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_loop_do_ot', 01044 routineP = moduleN//':'//routineN 01045 01046 INTEGER :: handle, ispin 01047 LOGICAL :: failure 01048 TYPE(cp_dbcsr_type), POINTER :: orthogonality_metric 01049 01050 CALL timeset(routineN,handle) 01051 01052 failure = .FALSE. 01053 01054 IF (has_unit_metric) THEN 01055 NULLIFY(orthogonality_metric) 01056 ELSE 01057 orthogonality_metric=>matrix_s(1)%matrix 01058 END IF 01059 01060 ! in case of LSD the first spin qs_ot_env will drive the minimization 01061 ! in the case of a restricted calculation, it will make sure the spin orbitals are equal 01062 01063 CALL ot_scf_mini(mos,mo_derivs,smear,orthogonality_metric, & 01064 total_energy, aux_coeff_set, & 01065 pscp,fscp,energy_only,scf_env%iter_delta, & 01066 scf_env%qs_ot_env,qs_env%input,error=error) 01067 01068 01069 DO ispin=1,SIZE(mos) 01070 CALL set_mo_occupation(mo_set=mos(ispin)%mo_set,& 01071 smear=smear,& 01072 error=error) 01073 ENDDO 01074 01075 DO ispin=1,SIZE(mos) 01076 CALL calculate_density_matrix(mos(ispin)%mo_set,& 01077 rho%rho_ao(ispin)%matrix,& 01078 use_dbcsr=.TRUE.,& 01079 error=error) 01080 END DO 01081 01082 scf_env%iter_method=scf_env%qs_ot_env(1)%OT_METHOD_FULL 01083 scf_env%iter_param=scf_env%qs_ot_env(1)%ds_min 01084 qs_env%broyden_adaptive_sigma=scf_env%qs_ot_env(1)%broyden_adaptive_sigma 01085 01086 CALL timestop(handle) 01087 01088 END SUBROUTINE qs_scf_loop_do_ot 01089 01090 ! ***************************************************************************** 01096 SUBROUTINE qs_scf_scp_print_summary(output_unit,scp_env,qmmm,error) 01097 INTEGER, INTENT(IN) :: output_unit 01098 TYPE(scp_environment_type), POINTER :: scp_env 01099 LOGICAL, INTENT(IN) :: qmmm 01100 TYPE(cp_error_type), INTENT(inout) :: error 01101 01102 TYPE(scp_energy_type), POINTER :: energy 01103 01104 CALL get_scp_env ( scp_env=scp_env,& 01105 energy=energy, error=error) 01106 IF (output_unit>0) THEN 01107 IF ( qmmm ) WRITE (UNIT=output_unit,FMT="((T3,A,T55,F25.14))")& 01108 "SCP Hartree (SCP density, QMMM potential) ",energy % e_scp_qmmm 01109 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T55,F25.14))")& 01110 "SCP Hartree (SCP,SCP) ",energy % e_scp_hartree,& 01111 "SCP Hartree Self (SCP,SCP) ",energy%e_scp_self,& 01112 "SCP Hartree (KS,SCP) ",energy % e_scp_ks,& 01113 "SCP Hartree Self (KS,SCP) ",energy%e_scp_ks_self,& 01114 "SCP Hartree Self Core (KS core,SCP) ",energy%e_scp_core,& 01115 "SCP Polarization Kernel ",energy % e_scp_kernel, & 01116 "SCP TOTAL ",energy % e_scp_total 01117 01118 CALL m_flush(output_unit) 01119 END IF 01120 END SUBROUTINE qs_scf_scp_print_summary 01121 01122 ! ***************************************************************************** 01128 SUBROUTINE qs_scf_print_summary(output_unit,rho,qs_charges,energy,nelectron_total,& 01129 dft_control,qmmm,qs_env,gapw,gapw_xc,error) 01130 01131 INTEGER, INTENT(IN) :: output_unit 01132 TYPE(qs_rho_type), POINTER :: rho 01133 TYPE(qs_charges_type), POINTER :: qs_charges 01134 TYPE(qs_energy_type), POINTER :: energy 01135 INTEGER, INTENT(IN) :: nelectron_total 01136 TYPE(dft_control_type), POINTER :: dft_control 01137 LOGICAL, INTENT(IN) :: qmmm 01138 TYPE(qs_environment_type), POINTER :: qs_env 01139 LOGICAL, INTENT(IN) :: gapw, gapw_xc 01140 TYPE(cp_error_type), INTENT(inout) :: error 01141 01142 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_print_summary', 01143 routineP = moduleN//':'//routineN 01144 01145 INTEGER :: handle, ispin 01146 REAL(kind=dp) :: tot1_h, tot1_s 01147 01148 CALL timeset(routineN,handle) 01149 01150 IF (output_unit>0) THEN 01151 IF(.NOT.(dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR.& 01152 dft_control%qs_control%scptb)) THEN 01153 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T41,2F20.10))")& 01154 "Electronic density on regular grids: ",& 01155 accurate_sum(rho%tot_rho_r),& 01156 accurate_sum(rho%tot_rho_r) + nelectron_total ,& 01157 "Core density on regular grids:",& 01158 qs_charges%total_rho_core_rspace,& 01159 qs_charges%total_rho_core_rspace - REAL(nelectron_total+dft_control%charge,dp) 01160 IF (gapw) THEN 01161 tot1_h = qs_charges%total_rho1_hard(1) 01162 tot1_s = qs_charges%total_rho1_soft(1) 01163 DO ispin=2,dft_control%nspins 01164 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin) 01165 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin) 01166 END DO 01167 WRITE (UNIT=output_unit,FMT="((T3,A,T41,2F20.10))")& 01168 "Hard and soft densities (Lebedev):",& 01169 tot1_h, tot1_s 01170 WRITE (UNIT=output_unit,FMT="(T3,A,T41,F20.10)")& 01171 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ",& 01172 accurate_sum(rho%tot_rho_r)+ tot1_h - tot1_s ,& 01173 "Total charge density (r-space): ",& 01174 accurate_sum(rho%tot_rho_r)+ tot1_h - tot1_s & 01175 + qs_charges%total_rho_core_rspace,& 01176 "Total Rho_soft + Rho0_soft (g-space):",& 01177 qs_charges%total_rho_gspace 01178 ELSE 01179 WRITE (UNIT=output_unit,FMT="(T3,A,T41,F20.10)")& 01180 "Total charge density on r-space grids: ",& 01181 accurate_sum(rho%tot_rho_r)+& 01182 qs_charges%total_rho_core_rspace,& 01183 "Total charge density g-space grids: ",& 01184 qs_charges%total_rho_gspace 01185 END IF 01186 END IF 01187 IF (dft_control%qs_control%semi_empirical) THEN 01188 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01189 "Core-core repulsion energy [eV]: ",energy%core_overlap*evolt,& 01190 "Core Hamiltonian energy [eV]: ",energy%core*evolt,& 01191 "Two-electron integral energy [eV]: ",energy%hartree*evolt,& 01192 "Electronic energy [eV]: ",& 01193 (energy%core+0.5_dp*energy%hartree)*evolt 01194 ELSEIF (dft_control%qs_control%dftb) THEN 01195 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01196 "Core Hamiltonian energy: ",energy%core,& 01197 "Repulsive potential energy: ",energy%repulsive,& 01198 "Electronic energy: ",energy%hartree,& 01199 "Dispersion energy: ",energy%dispersion 01200 ELSEIF (dft_control%qs_control%scptb) THEN 01201 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01202 "Repulsive pair potential energy: ",energy%repulsive,& 01203 "Zeroth order Hamiltonian energy: ",energy%core,& 01204 "Kinetic energy: ",energy%kinetic,& 01205 "Charge fluctuation (SCP) energy: ",energy%hartree,& 01206 "London dispersion energy: ",energy%dispersion 01207 ELSE 01208 IF( dft_control%do_admm ) THEN 01209 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01210 "Overlap energy of the core charge distribution:",energy%core_overlap,& 01211 "Self energy of the core charge distribution: ",energy%core_self,& 01212 "Core Hamiltonian energy: ",energy%core,& 01213 "Hartree energy: ",energy%hartree,& 01214 "Exchange-correlation energy: ",energy%exc+energy%exc_aux_fit 01215 ELSE 01216 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01217 "Overlap energy of the core charge distribution:",energy%core_overlap,& 01218 "Self energy of the core charge distribution: ",energy%core_self,& 01219 "Core Hamiltonian energy: ",energy%core,& 01220 "Hartree energy: ",energy%hartree,& 01221 "Exchange-correlation energy: ",energy%exc 01222 END IF 01223 IF (energy%e_hartree /= 0.0_dp) & 01224 WRITE (UNIT=output_unit,FMT="(T3,A,/,T3,A,T56,F25.14)")& 01225 "Coulomb Electron-Electron Interaction Energy ",& 01226 "- Already included in the total Hartree term ",energy%e_hartree 01227 IF (energy%ex /= 0.0_dp)& 01228 WRITE (UNIT=output_unit,FMT="(T3,A,T56,F25.14)")& 01229 "Hartree-Fock Exchange energy: ",energy%ex 01230 IF (energy%dispersion /= 0.0_dp)& 01231 WRITE (UNIT=output_unit,FMT="(T3,A,T56,F25.14)")& 01232 "Dispersion energy: ",energy%dispersion 01233 IF(gapw) THEN 01234 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01235 "GAPW| Exc from hard and soft atomic rho1: ",energy%exc1,& 01236 "GAPW| local Eh = 1 center integrals: ",energy%hartree_1c 01237 END IF 01238 IF(gapw_xc) THEN 01239 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01240 "GAPW_XC| Exc from hard and soft atomic rho1: ",energy%exc1 01241 END IF 01242 END IF 01243 IF (dft_control%smear) THEN 01244 WRITE (UNIT=output_unit,FMT="((T3,A,T56,F25.14))")& 01245 "Electronic entropic energy:",energy%kTS 01246 WRITE (UNIT=output_unit,FMT="((T3,A,T56,F25.14))")& 01247 "Fermi energy:",energy%efermi 01248 END IF 01249 IF (dft_control%dft_plus_u) THEN 01250 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01251 "DFT+U energy:",energy%dft_plus_u 01252 END IF 01253 IF (qmmm) THEN 01254 WRITE (UNIT=output_unit,FMT="(T3,A,T61,F20.10)")& 01255 "QM/MM Electrostatic energy: ",energy%qmmm_el 01256 IF(qs_env%qmmm_env_qm%image_charge) THEN 01257 WRITE (UNIT=output_unit,FMT="(T3,A,T61,F20.10)")& 01258 "QM/MM image charge energy: ",energy%image_charge 01259 ENDIF 01260 END IF 01261 IF (dft_control%qs_control%mulliken_restraint) THEN 01262 WRITE (UNIT=output_unit,FMT="(T3,A,T61,F20.10)")& 01263 "Mulliken restraint energy: ",energy%mulliken 01264 END IF 01265 IF (dft_control%qs_control%semi_empirical) THEN 01266 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01267 "Total energy [eV]: ",energy%total*evolt 01268 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01269 "Atomic reference energy [eV]: ",energy%core_self*evolt,& 01270 "Heat of formation [kcal/mol]: ",& 01271 (energy%total+energy%core_self)*kcalmol 01272 ELSE 01273 WRITE (UNIT=output_unit,FMT="(/,(T3,A,T56,F25.14))")& 01274 "Total energy: ",energy%total 01275 END IF 01276 IF (qmmm) THEN 01277 IF(qs_env%qmmm_env_qm%image_charge) THEN 01278 CALL print_image_coefficients(qs_env%image_coeff,qs_env,error) 01279 ENDIF 01280 ENDIF 01281 CALL m_flush(output_unit) 01282 END IF 01283 01284 CALL timestop(handle) 01285 01286 END SUBROUTINE qs_scf_print_summary 01287 01288 ! ***************************************************************************** 01293 SUBROUTINE qs_scf_loop_print(qs_env,scf_env,para_env,error) 01294 TYPE(qs_environment_type), POINTER :: qs_env 01295 TYPE(qs_scf_env_type), POINTER :: scf_env 01296 TYPE(cp_para_env_type), POINTER :: para_env 01297 TYPE(cp_error_type), INTENT(inout) :: error 01298 01299 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_loop_print', 01300 routineP = moduleN//':'//routineN 01301 01302 INTEGER :: handle, ispin, iw 01303 REAL(KIND=dp) :: mo_mag_max, mo_mag_min, 01304 orthonormality 01305 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01306 POINTER :: matrix_ks, matrix_p, matrix_s 01307 TYPE(cp_logger_type), POINTER :: logger 01308 TYPE(dft_control_type), POINTER :: dft_control 01309 TYPE(mo_set_p_type), DIMENSION(:), 01310 POINTER :: mos 01311 TYPE(qs_rho_type), POINTER :: rho 01312 TYPE(section_vals_type), POINTER :: dft_section, input, 01313 scf_section 01314 01315 logger => cp_error_get_logger(error) 01316 CALL timeset(routineN,handle) 01317 01318 CALL get_qs_env(qs_env=qs_env,input=input, dft_control=dft_control, & 01319 error=error) 01320 01321 dft_section => section_vals_get_subs_vals(input,"DFT",error=error) 01322 scf_section => section_vals_get_subs_vals(dft_section,"SCF",error=error) 01323 01324 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, & 01325 matrix_ks=matrix_ks,rho=rho,mos=mos, & 01326 error=error) 01327 01328 matrix_p => rho%rho_ao 01329 01330 DO ispin=1,dft_control%nspins 01331 01332 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01333 dft_section,"PRINT%AO_MATRICES/DENSITY",error=error),cp_p_file)) THEN 01334 iw = cp_print_key_unit_nr(logger,dft_section,"PRINT%AO_MATRICES/DENSITY",& 01335 extension=".Log",error=error) 01336 CALL cp_dbcsr_write_sparse_matrix(matrix_p(ispin)%matrix,4,6,qs_env,para_env,& 01337 output_unit=iw,error=error) 01338 CALL cp_print_key_finished_output(iw,logger,dft_section,& 01339 "PRINT%AO_MATRICES/DENSITY", error=error) 01340 END IF 01341 01342 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01343 dft_section,"PRINT%AO_MATRICES/KOHN_SHAM_MATRIX",error=error),cp_p_file)) THEN 01344 iw = cp_print_key_unit_nr(logger,dft_section,"PRINT%AO_MATRICES/KOHN_SHAM_MATRIX",& 01345 extension=".Log",error=error) 01346 IF (dft_control%qs_control%semi_empirical) THEN 01347 CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin)%matrix,4,6,qs_env,para_env,& 01348 scale=evolt,output_unit=iw,error=error) 01349 ELSE 01350 CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin)%matrix,4,6,qs_env,para_env,& 01351 output_unit=iw,error=error) 01352 END IF 01353 CALL cp_print_key_finished_output(iw,logger,dft_section,& 01354 "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", error=error) 01355 END IF 01356 01357 ENDDO 01358 01359 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01360 scf_section,"PRINT%MO_ORTHONORMALITY",error=error),cp_p_file)) THEN 01361 IF(scf_env%method==special_diag_method_nr) THEN 01362 CALL calculate_orthonormality(orthonormality,mos,error=error) 01363 ELSE 01364 CALL calculate_orthonormality(orthonormality,mos,matrix_s(1)%matrix,error=error) 01365 END IF 01366 iw=cp_print_key_unit_nr(logger,scf_section,"PRINT%MO_ORTHONORMALITY",& 01367 extension=".scfLog",error=error) 01368 IF (iw>0) THEN 01369 WRITE(iw,'(T8,A,T61,E20.4)') & 01370 " Maximum deviation from MO S-orthonormality",orthonormality 01371 ENDIF 01372 CALL cp_print_key_finished_output(iw,logger,scf_section,& 01373 "PRINT%MO_ORTHONORMALITY", error=error) 01374 ENDIF 01375 01376 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01377 scf_section,"PRINT%MO_MAGNITUDE",error=error),cp_p_file)) THEN 01378 CALL calculate_magnitude(mos,mo_mag_min,mo_mag_max,error=error) 01379 iw=cp_print_key_unit_nr(logger,scf_section,"PRINT%MO_MAGNITUDE",& 01380 extension=".scfLog",error=error) 01381 IF (iw>0) THEN 01382 WRITE(iw,'(T8,A,T41,2E20.4)') & 01383 " Minimum/Maximum MO magnitude ",mo_mag_min,mo_mag_max 01384 ENDIF 01385 CALL cp_print_key_finished_output(iw,logger,scf_section,& 01386 "PRINT%MO_MAGNITUDE", error=error) 01387 ENDIF 01388 01389 CALL timestop(handle) 01390 01391 END SUBROUTINE qs_scf_loop_print 01392 01393 ! ***************************************************************************** 01400 SUBROUTINE init_scf_loop(scf_env,qs_env,scf_section,error) 01401 01402 TYPE(qs_scf_env_type), POINTER :: scf_env 01403 TYPE(qs_environment_type), POINTER :: qs_env 01404 TYPE(section_vals_type), POINTER :: scf_section 01405 TYPE(cp_error_type), INTENT(inout) :: error 01406 01407 CHARACTER(LEN=*), PARAMETER :: routineN = 'init_scf_loop', 01408 routineP = moduleN//':'//routineN 01409 01410 INTEGER :: handle, ispin, nmo, 01411 number_of_OT_envs, stat 01412 LOGICAL :: do_rotation, failure, 01413 has_unit_metric, scp_dft, 01414 scp_nddo 01415 TYPE(aux_coeff_set_type), POINTER :: aux_coeff_set 01416 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01417 POINTER :: matrix_ks, matrix_s 01418 TYPE(cp_dbcsr_type), POINTER :: ks_scp, orthogonality_metric, 01419 pscp 01420 TYPE(cp_fm_type), POINTER :: mo_coeff 01421 TYPE(dft_control_type), POINTER :: dft_control 01422 TYPE(mo_set_p_type), DIMENSION(:), 01423 POINTER :: mos 01424 TYPE(scf_control_type), POINTER :: scf_control 01425 TYPE(scp_environment_type), POINTER :: scp_env 01426 01427 CALL timeset(routineN,handle) 01428 01429 NULLIFY(scf_control,matrix_s,matrix_ks,dft_control,mos,mo_coeff) 01430 ! **SCP 01431 NULLIFY ( scp_env, aux_coeff_set, pscp, ks_scp ) 01432 ! **SCP 01433 01434 failure=.FALSE. 01435 01436 CPPrecondition(ASSOCIATED(scf_env),cp_failure_level,routineP,error,failure) 01437 CPPrecondition(scf_env%ref_count>0,cp_failure_level,routineP,error,failure) 01438 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 01439 CPPrecondition(qs_env%ref_count>0,cp_failure_level,routineP,error,failure) 01440 01441 CALL get_qs_env(qs_env=qs_env,& 01442 scf_control=scf_control,& 01443 dft_control=dft_control,& 01444 mos=mos,matrix_ks=matrix_ks,& 01445 matrix_s=matrix_s, error=error) 01446 01447 scp_dft = dft_control%scp 01448 scp_nddo = dft_control%qs_control%se_control%scp 01449 01450 ! if using mo_coeff_b then copy to fm 01451 DO ispin=1,SIZE(mos)!fm->dbcsr 01452 IF(mos(1)%mo_set%use_mo_coeff_b)THEN!fm->dbcsr 01453 CALL copy_dbcsr_to_fm(mos(ispin)%mo_set%mo_coeff_b,mos(ispin)%mo_set%mo_coeff,error=error)!fm->dbcsr 01454 ENDIF!fm->dbcsr 01455 ENDDO!fm->dbcsr 01456 01457 ! this just guarantees that all mo_occupations match the eigenvalues, if smear 01458 DO ispin=1,dft_control%nspins 01459 CALL set_mo_occupation(mo_set=mos(ispin)%mo_set,& 01460 smear=scf_control%smear, error=error) 01461 ENDDO 01462 01463 SELECT CASE (scf_env%method) 01464 CASE DEFAULT 01465 01466 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 01467 routineP,"unknown scf method method:"//& 01468 cp_to_string(scf_env%method),error,failure) 01469 01470 CASE (general_diag_method_nr,special_diag_method_nr,block_krylov_diag_method_nr) 01471 IF(.NOT.scf_env%skip_diis) THEN 01472 IF (.NOT.ASSOCIATED(scf_env%scf_diis_buffer)) THEN 01473 CALL qs_diis_b_create(scf_env%scf_diis_buffer,nbuffer=scf_control%max_diis,error=error) 01474 END IF 01475 CALL qs_diis_b_clear(scf_env%scf_diis_buffer,error=error) 01476 END IF 01477 01478 CASE (ot_diag_method_nr) 01479 01480 IF(.NOT.scf_env%skip_diis) THEN 01481 IF (.NOT.ASSOCIATED(scf_env%scf_diis_buffer)) THEN 01482 CALL qs_diis_b_create(scf_env%scf_diis_buffer,nbuffer=scf_control%max_diis,error=error) 01483 END IF 01484 CALL qs_diis_b_clear(scf_env%scf_diis_buffer,error=error) 01485 END IF 01486 01487 ! disable DFTB and SE for now 01488 IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%semi_empirical) THEN 01489 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 01490 routineP,"DFTB and SE not available with OT/DIAG",error,failure) 01491 END IF 01492 01493 ! if an old preconditioner is still around (i.e. outer SCF is active), 01494 ! remove it if this could be worthwhile 01495 CALL restart_preconditioner(qs_env,scf_env%ot_preconditioner,& 01496 scf_control%diagonalization%ot_settings%preconditioner_type,& 01497 dft_control%nspins,error=error) 01498 01499 CALL prepare_preconditioner(qs_env,mos,matrix_ks,matrix_s,scf_env%ot_preconditioner,& 01500 scf_control%diagonalization%ot_settings%preconditioner_type,& 01501 scf_control%diagonalization%ot_settings%precond_solver_type,& 01502 scf_control%diagonalization%ot_settings%energy_gap,dft_control%nspins,error=error) 01503 CASE (block_davidson_diag_method_nr) 01504 ! Preconditioner initialized within the loop, when required 01505 CASE (ot_method_nr) 01506 01507 CALL get_qs_env(qs_env,has_unit_metric=has_unit_metric,error=error) 01508 01509 ! reortho the wavefunctions if we are having an outer scf and 01510 ! this is not the first iteration 01511 ! this is useful to avoid the build-up of numerical noise 01512 ! however, we can not play this trick if restricted (don't mix non-equivalent orbs) 01513 IF(scf_control%do_outer_scf_reortho) THEN 01514 IF (scf_control%outer_scf%have_scf .AND. .NOT. dft_control%restricted) THEN 01515 IF (scf_env%outer_scf%iter_count>0) THEN 01516 DO ispin=1,dft_control%nspins 01517 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, nmo=nmo) 01518 IF (has_unit_metric) THEN 01519 CALL make_basis_simple(mo_coeff,nmo,error=error) 01520 ELSE 01521 CALL make_basis_sm(mo_coeff,nmo,matrix_s(1)%matrix,error=error) 01522 ENDIF 01523 ENDDO 01524 ENDIF 01525 ENDIF 01526 ELSE 01527 ! dont need any dirty trick for the numerically stable irac algorithm. 01528 ENDIF 01529 01530 IF (.NOT.ASSOCIATED(scf_env%qs_ot_env)) THEN 01531 01532 ! restricted calculations require just one set of OT orbitals 01533 number_of_OT_envs=dft_control%nspins 01534 IF (dft_control%restricted) number_of_OT_envs=1 01535 01536 ALLOCATE(scf_env%qs_ot_env(number_of_OT_envs),stat=stat) 01537 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01538 01539 ! XXX Joost XXX should disentangle reading input from this part 01540 CALL ot_scf_read_input(scf_env%qs_ot_env,scf_section,error) 01541 01542 ! keep a note that we are restricted 01543 IF (dft_control%restricted) THEN 01544 scf_env%qs_ot_env(1)%restricted=.TRUE. 01545 ! requires rotation 01546 CALL cp_assert(scf_env%qs_ot_env(1)%settings%do_rotation,cp_fatal_level,& 01547 cp_assertion_failed,routineP,& 01548 "Restricted calculation with OT requires orbital rotation. Please "//& 01549 "activate the OT%ROTATION keyword! "//& 01550 CPSourceFileRef) 01551 ELSE 01552 scf_env%qs_ot_env(:)%restricted=.FALSE. 01553 ENDIF 01554 01555 ! might need the KS matrix to init properly 01556 ! **SCP 01557 ! Updates density and potential due to SCP 01558 IF ( scp_dft ) CALL qs_ks_scp_update ( qs_env, just_energy=.FALSE., error=error ) 01559 ! **SCP 01560 CALL qs_ks_update_qs_env(qs_env%ks_env,& 01561 qs_env=qs_env,& 01562 calculate_forces=.FALSE.,& 01563 just_energy=.FALSE.,& 01564 error=error) 01565 01566 ! if an old preconditioner is still around (i.e. outer SCF is active), 01567 ! remove it if this could be worthwhile 01568 CALL restart_preconditioner(qs_env,scf_env%ot_preconditioner,& 01569 scf_env%qs_ot_env(1)%settings%preconditioner_type,& 01570 dft_control%nspins,error=error) 01571 01572 ! 01573 ! preconditioning still needs to be done correctly with has_unit_metric 01574 ! notice that a big part of the preconditioning (S^-1) is fine anyhow 01575 ! 01576 IF (has_unit_metric) THEN 01577 NULLIFY(orthogonality_metric) 01578 ELSE 01579 orthogonality_metric=>matrix_s(1)%matrix 01580 ENDIF 01581 IF(dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb) THEN 01582 scf_env%qs_ot_env(1)%settings%mixed_precision=.FALSE. 01583 END IF 01584 !new 01585 CALL prepare_preconditioner(qs_env,mos,matrix_ks,matrix_s,scf_env%ot_preconditioner,& 01586 scf_env%qs_ot_env(1)%settings%preconditioner_type,& 01587 scf_env%qs_ot_env(1)%settings%precond_solver_type,& 01588 scf_env%qs_ot_env(1)%settings%energy_gap,dft_control%nspins,& 01589 has_unit_metric=has_unit_metric,& 01590 mixed_precision=scf_env%qs_ot_env(1)%settings%mixed_precision,error=error) 01591 !new 01592 IF(dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb) THEN 01593 ! we have to treat the semi-empirical methods separately 01594 IF ( scp_nddo ) THEN 01595 CALL get_qs_env ( qs_env = qs_env, scp_env = scp_env, error = error ) 01596 CALL get_scp_env ( scp_env = scp_env, pscp=pscp, ks_scp=ks_scp, error = error ) 01597 CALL ot_scf_init(mos,orthogonality_metric, pscp=pscp, fscp=ks_scp, & 01598 broyden_adaptive_sigma=qs_env%broyden_adaptive_sigma,& 01599 qs_ot_env=scf_env%qs_ot_env, matrix_ks=matrix_ks(1)%matrix, & 01600 error=error) 01601 ELSE 01602 CALL ot_scf_init(mo_array = mos, matrix_s = orthogonality_metric, & 01603 broyden_adaptive_sigma=qs_env%broyden_adaptive_sigma,& 01604 qs_ot_env = scf_env%qs_ot_env,matrix_ks=matrix_ks(1)%matrix,& 01605 error=error) 01606 END IF 01607 ELSE 01608 ! **** SCP 01609 IF ( scp_dft ) THEN 01610 CALL get_qs_env ( qs_env = qs_env, scp_env = scp_env, error = error ) 01611 CALL get_scp_env ( scp_env = scp_env, aux_coeff_set = aux_coeff_set, error = error ) 01612 CALL ot_scf_init(mos,orthogonality_metric, aux_coeff_set=aux_coeff_set, & 01613 broyden_adaptive_sigma=qs_env%broyden_adaptive_sigma,& 01614 qs_ot_env=scf_env%qs_ot_env, matrix_ks=matrix_ks(1)%matrix, & 01615 error=error) 01616 ELSE 01617 CALL ot_scf_init(mo_array = mos, matrix_s = orthogonality_metric, & 01618 broyden_adaptive_sigma=qs_env%broyden_adaptive_sigma,& 01619 qs_ot_env = scf_env%qs_ot_env, matrix_ks=matrix_ks(1)%matrix, & 01620 error=error) 01621 ENDIF 01622 ! **** SCP 01623 END IF 01624 01625 SELECT CASE(scf_env%qs_ot_env(1)%settings%preconditioner_type) 01626 CASE(ot_precond_none) 01627 CASE(ot_precond_full_all,ot_precond_full_single_inverse) 01628 ! this will rotate the MOs to be eigen states, which is not compatible with rotation 01629 ! e.g. mo_derivs here do not yet include potentially different occupations numbers 01630 do_rotation=scf_env%qs_ot_env(1)%settings%do_rotation 01631 CPPrecondition(.NOT.do_rotation,cp_failure_level,routineP,error,failure) 01632 DO ispin=1,SIZE(scf_env%qs_ot_env) 01633 CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin),& 01634 scf_env%ot_preconditioner(ispin)%preconditioner,error=error) 01635 ENDDO 01636 CASE(ot_precond_s_inverse,ot_precond_sparse_diag,ot_precond_full_single) 01637 DO ispin=1,SIZE(scf_env%qs_ot_env) 01638 CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin),& 01639 scf_env%ot_preconditioner(1)%preconditioner,error=error) 01640 ENDDO 01641 CASE DEFAULT 01642 DO ispin=1,SIZE(scf_env%qs_ot_env) 01643 CALL qs_ot_new_preconditioner(scf_env%qs_ot_env(ispin),& 01644 scf_env%ot_preconditioner(1)%preconditioner,error=error) 01645 ENDDO 01646 END SELECT 01647 ENDIF 01648 01649 ! if we have non-uniform occupations we should be using rotation 01650 do_rotation=scf_env%qs_ot_env(1)%settings%do_rotation 01651 DO ispin=1,SIZE(mos) 01652 IF (.NOT. mos(ispin)%mo_set%uniform_occupation) THEN 01653 CPPrecondition(do_rotation,cp_failure_level,routineP,error,failure) 01654 ENDIF 01655 ENDDO 01656 END SELECT 01657 01658 ! another safety check 01659 IF (dft_control%low_spin_roks) THEN 01660 CPPrecondition(scf_env%method==ot_method_nr,cp_failure_level,routineP,error,failure) 01661 do_rotation=scf_env%qs_ot_env(1)%settings%do_rotation 01662 CPPrecondition(do_rotation,cp_failure_level,routineP,error,failure) 01663 ENDIF 01664 01665 CALL timestop(handle) 01666 01667 END SUBROUTINE init_scf_loop 01668 01669 ! ***************************************************************************** 01676 SUBROUTINE init_scf_run(scf_env, qs_env, scf_section, do_xas_tp, error) 01677 01678 TYPE(qs_scf_env_type), POINTER :: scf_env 01679 TYPE(qs_environment_type), POINTER :: qs_env 01680 TYPE(section_vals_type), POINTER :: scf_section 01681 LOGICAL, INTENT(IN), OPTIONAL :: do_xas_tp 01682 TYPE(cp_error_type), INTENT(inout) :: error 01683 01684 CHARACTER(LEN=*), PARAMETER :: routineN = 'init_scf_run', 01685 routineP = moduleN//':'//routineN 01686 01687 INTEGER :: handle, homo, ikind, ispin, 01688 iw, nao, ndep, 01689 nelectron_spin, nmo, 01690 output_unit 01691 LOGICAL :: dft_plus_u_atom, failure, gth_potential_present, 01692 init_u_ramping_each_scf, my_transition_potential, s_minus_half_available 01693 REAL(KIND=dp) :: u_ramping 01694 TYPE(atomic_kind_type), DIMENSION(:), 01695 POINTER :: atomic_kind_set 01696 TYPE(atomic_kind_type), POINTER :: atomic_kind 01697 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01698 POINTER :: matrix_h, matrix_ks, matrix_s 01699 TYPE(cp_logger_type), POINTER :: logger 01700 TYPE(cp_para_env_type), POINTER :: para_env 01701 TYPE(dft_control_type), POINTER :: dft_control 01702 TYPE(mo_set_p_type), DIMENSION(:), 01703 POINTER :: mos 01704 TYPE(qs_rho_type), POINTER :: rho 01705 TYPE(scf_control_type), POINTER :: scf_control 01706 01707 CALL timeset(routineN,handle) 01708 01709 NULLIFY(scf_control, atomic_kind_set, matrix_h, matrix_s, matrix_ks, & 01710 dft_control, mos, rho) 01711 failure=.FALSE. 01712 my_transition_potential = .FALSE. 01713 IF(PRESENT(do_xas_tp)) my_transition_potential = do_xas_tp 01714 logger => cp_error_get_logger(error) 01715 output_unit = cp_print_key_unit_nr(logger,scf_section,"PRINT%PROGRAM_RUN_INFO",& 01716 extension=".scfLog",error=error) 01717 01718 CPPrecondition(ASSOCIATED(scf_env),cp_failure_level,routineP,error,failure) 01719 CPPrecondition(scf_env%ref_count>0,cp_failure_level,routineP,error,failure) 01720 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 01721 CPPrecondition(qs_env%ref_count>0,cp_failure_level,routineP,error,failure) 01722 NULLIFY(para_env) 01723 para_env=>qs_env%para_env 01724 01725 s_minus_half_available = .FALSE. 01726 01727 CALL scf_env_check_i_alloc(scf_env=scf_env, qs_env=qs_env,& 01728 do_xas_tp=my_transition_potential,error=error) 01729 01730 CALL get_qs_env(qs_env=qs_env,& 01731 scf_control=scf_control,& 01732 dft_control=dft_control,& 01733 atomic_kind_set=atomic_kind_set,& 01734 mos=mos,matrix_ks=matrix_ks, rho=rho,& 01735 matrix_h=matrix_h,matrix_s=matrix_s,& 01736 error=error) 01737 01738 ! some special checks to eliminate the first problems with restricted 01739 ! should move earlier I think 01740 IF (dft_control%restricted) THEN 01741 IF (scf_env%method .NE. ot_method_nr) THEN 01742 IF (output_unit>0) WRITE(output_unit,*) "OT only for restricted" 01743 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 01744 ENDIF 01745 CPPrecondition(dft_control%nspins.EQ.2,cp_failure_level,routineP,error,failure) 01746 ENDIF 01747 01748 IF((scf_env%method == ot_method_nr) .AND. & 01749 (scf_control%added_mos(1)+scf_control%added_mos(2)) > 0) THEN 01750 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 01751 routineP,"OT with ADDED_MOS/=0 not implemented", error,failure) 01752 END IF 01753 01754 01755 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 01756 gth_potential_present=gth_potential_present) 01757 01758 ! updates the total number of electrons 01759 CALL get_qs_env(qs_env,nelectron_total=scf_env%nelectron,error=error) 01760 01761 ! print occupation numbers 01762 IF (output_unit > 0) THEN 01763 DO ispin=1,dft_control%nspins 01764 CALL get_mo_set(mo_set=mos(ispin)%mo_set,& 01765 homo=homo, & 01766 nelectron=nelectron_spin, & 01767 nao=nao,& 01768 nmo=nmo) 01769 IF (dft_control%nspins > 1) THEN 01770 WRITE (UNIT=output_unit,FMT="(/,T2,A,I2)") "Spin",ispin 01771 END IF 01772 WRITE (UNIT=output_unit,FMT="(/,(T2,A,T71,I10))")& 01773 "Number of electrons:",nelectron_spin,& 01774 "Number of occupied orbitals:",homo,& 01775 "Number of molecular orbitals:",nmo 01776 END DO 01777 WRITE (UNIT=output_unit,FMT="(/,T2,A,T71,I10)")& 01778 "Number of orbital functions:",nao 01779 END IF 01780 CALL cp_print_key_finished_output(output_unit,logger,scf_section,& 01781 "PRINT%PROGRAM_RUN_INFO", error=error) 01782 01783 ! calc ortho matrix 01784 ndep = 0 01785 IF (scf_env%method /= ot_method_nr .AND. & 01786 scf_env%method /= ot_diag_method_nr .AND. & 01787 scf_env%method /= special_diag_method_nr .AND.& 01788 scf_env%method /= block_davidson_diag_method_nr) THEN 01789 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix,scf_env%ortho,error=error) 01790 IF (scf_env%cholesky_method>cholesky_off) THEN 01791 CALL cp_fm_cholesky_decompose(scf_env%ortho,error=error) 01792 IF(scf_env%cholesky_method==cholesky_inverse) THEN 01793 CALL cp_fm_to_fm(scf_env%ortho,scf_env%ortho_m1,error=error) 01794 CALL cp_fm_triangular_invert(scf_env%ortho_m1,error=error) 01795 END IF 01796 ELSE 01797 CALL cp_fm_power(scf_env%ortho,scf_env%scf_work2,-0.5_dp,& 01798 scf_control%eps_eigval,ndep,error=error) 01799 s_minus_half_available = .TRUE. 01800 END IF 01801 01802 IF (BTEST(cp_print_key_should_output(logger%iter_info,& 01803 qs_env%input,"DFT%PRINT%AO_MATRICES/ORTHO",error=error),cp_p_file)) THEN 01804 iw = cp_print_key_unit_nr(logger,qs_env%input,"DFT%PRINT%AO_MATRICES/ORTHO",& 01805 extension=".Log",error=error) 01806 CALL write_fm_with_basis_info(scf_env%ortho,4,6,qs_env,para_env,output_unit=iw,error=error) 01807 CALL cp_print_key_finished_output(iw,logger,qs_env%input,& 01808 "DFT%PRINT%AO_MATRICES/ORTHO", error=error) 01809 END IF 01810 END IF 01811 01812 CALL get_mo_set(mo_set=mos(1)%mo_set,nao=nao) 01813 01814 ! DFT+U methods based on Lowdin charges need S^(1/2) 01815 IF (dft_control%dft_plus_u) THEN 01816 IF (dft_control%plus_u_method_id == plus_u_lowdin) THEN 01817 IF (s_minus_half_available) THEN 01818 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix,scf_env%ortho,scf_env%s_half,& 01819 nao,error=error) 01820 ELSE 01821 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix,scf_env%s_half,error=error) 01822 CALL cp_fm_power(scf_env%s_half,scf_env%scf_work2,0.5_dp,& 01823 scf_control%eps_eigval,ndep,error=error) 01824 END IF 01825 END IF 01826 DO ikind=1,SIZE(atomic_kind_set) 01827 atomic_kind => atomic_kind_set(ikind) 01828 CALL get_atomic_kind(atomic_kind=atomic_kind,& 01829 dft_plus_u_atom=dft_plus_u_atom,& 01830 u_ramping=u_ramping,& 01831 init_u_ramping_each_scf=init_u_ramping_each_scf) 01832 IF (dft_plus_u_atom.AND.(u_ramping /= 0.0_dp)) THEN 01833 IF (init_u_ramping_each_scf) THEN 01834 CALL set_atomic_kind(atomic_kind=atomic_kind,& 01835 u_minus_j=0.0_dp) 01836 END IF 01837 END IF 01838 END DO 01839 END IF 01840 01841 output_unit=cp_print_key_unit_nr(logger,scf_section,"PRINT%PROGRAM_RUN_INFO",& 01842 extension=".scfLog",error=error) 01843 IF (output_unit > 0) THEN 01844 WRITE (UNIT=output_unit,FMT="(T2,A,T71,I10)")& 01845 "Number of independent orbital functions:",nao - ndep 01846 END IF 01847 CALL cp_print_key_finished_output(output_unit,logger,scf_section,& 01848 "PRINT%PROGRAM_RUN_INFO", error=error) 01849 01850 ! extrapolate outer loop variables 01851 IF (scf_control%outer_scf%have_scf) THEN 01852 CALL outer_loop_extrapolate(qs_env,error) 01853 ENDIF 01854 01855 ! initializes rho and the mos 01856 IF( my_transition_potential) THEN 01857 ! initialize the density with the localized mos 01858 CALL xas_initialize_rho(qs_env,error=error) 01859 ELSE 01860 CALL scf_env_initial_rho_setup(scf_env,qs_env=qs_env,& 01861 scf_section=scf_section, error=error) 01862 END IF 01863 01864 ! Frozen density approximation 01865 IF (ASSOCIATED(qs_env%wf_history)) THEN 01866 IF (qs_env%wf_history%interpolation_method_nr==wfi_frozen_method_nr) THEN 01867 IF (.NOT. ASSOCIATED(qs_env%wf_history%past_states(1)%snapshot)) THEN 01868 CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp, & 01869 error=error) 01870 CALL duplicate_rho_type(rho_input=rho, & 01871 rho_output=qs_env%wf_history%past_states(1)%snapshot%rho_frozen, & 01872 qs_env=qs_env, error=error) 01873 END IF 01874 END IF 01875 END IF 01876 01877 !image charge method, calculate image_matrix 01878 IF(qs_env%qmmm) THEN 01879 IF(qs_env%qmmm_env_qm%image_charge) THEN 01880 IF(.NOT.qs_env%qmmm_env_qm%image_charge_pot%coeff_iterative) & 01881 CALL calculate_image_matrix(image_matrix=qs_env%image_matrix,& 01882 ipiv=qs_env%ipiv,qs_env=qs_env,qmmm_env=qs_env%qmmm_env_qm,& 01883 error=error) 01884 ENDIF 01885 ENDIF 01886 01887 CALL timestop(handle) 01888 01889 END SUBROUTINE init_scf_run 01890 01891 ! ***************************************************************************** 01902 SUBROUTINE scf_env_check_i_alloc(scf_env,qs_env,do_xas_tp,error) 01903 TYPE(qs_scf_env_type), POINTER :: scf_env 01904 TYPE(qs_environment_type), POINTER :: qs_env 01905 LOGICAL, INTENT(IN) :: do_xas_tp 01906 TYPE(cp_error_type), INTENT(inout) :: error 01907 01908 CHARACTER(LEN=*), PARAMETER :: routineN = 'scf_env_check_i_alloc', 01909 routineP = moduleN//':'//routineN 01910 01911 CHARACTER(LEN=default_string_length) :: headline 01912 INTEGER :: handle, homo, ispin, 01913 mixing_method, nao, nhistory, 01914 nmo, nrow_block, nvariables, 01915 stat 01916 LOGICAL :: failure, gth_potential_present 01917 TYPE(array_i1d_obj) :: col_blk_size, col_dist 01918 TYPE(atomic_kind_type), DIMENSION(:), 01919 POINTER :: atomic_kind_set 01920 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01921 POINTER :: matrix_ks, matrix_ks_aux_fit, 01922 matrix_s, mo_derivs 01923 TYPE(cp_dbcsr_type), POINTER :: mo_coeff_b 01924 TYPE(cp_fm_p_type), DIMENSION(:), 01925 POINTER :: mo_derivs_aux_fit 01926 TYPE(cp_fm_pool_p_type), DIMENSION(:), 01927 POINTER :: ao_mo_fm_pools, 01928 ao_mo_fm_pools_aux_fit 01929 TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct, ao_mo_fmstruct 01930 TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_aux_fit 01931 TYPE(cp_para_env_type), POINTER :: para_env 01932 TYPE(dbcsr_distribution_obj) :: dist 01933 TYPE(dft_control_type), POINTER :: dft_control 01934 TYPE(mixing_storage_type), POINTER :: mixing_store 01935 TYPE(mo_set_p_type), DIMENSION(:), 01936 POINTER :: mos, mos_aux_fit 01937 TYPE(qs_ks_env_type), POINTER :: ks_env 01938 TYPE(scf_control_type), POINTER :: scf_control 01939 TYPE(xas_environment_type), POINTER :: xas_env 01940 01941 CALL timeset(routineN,handle) 01942 01943 NULLIFY(matrix_ks, ao_mo_fm_pools, matrix_s, ao_mo_fmstruct, ao_ao_fmstruct,& 01944 dft_control, mos, ks_env,& 01945 ao_mo_fm_pools_aux_fit, matrix_ks_aux_fit) 01946 NULLIFY(atomic_kind_set, mo_derivs, mo_coeff, scf_control, & 01947 mo_coeff_aux_fit, mo_derivs_aux_fit, para_env, xas_env) 01948 01949 failure=.FALSE. 01950 01951 CPPrecondition(ASSOCIATED(scf_env),cp_failure_level,routineP,error,failure) 01952 CPPrecondition(scf_env%ref_count>0,cp_failure_level,routineP,error,failure) 01953 01954 CALL get_qs_env(qs_env=qs_env,& 01955 dft_control=dft_control,& 01956 mos=mos,& 01957 mos_aux_fit=mos_aux_fit,& 01958 matrix_ks=matrix_ks,& 01959 ks_env=ks_env,& 01960 atomic_kind_set=atomic_kind_set,& 01961 matrix_s=matrix_s,& 01962 para_env=para_env, xas_env=xas_env,error=error) 01963 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& 01964 gth_potential_present=gth_potential_present) 01965 CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools,& 01966 error=error) 01967 IF( dft_control%do_admm) THEN 01968 CALL mpools_get(qs_env%mpools_aux_fit, ao_mo_fm_pools=ao_mo_fm_pools_aux_fit,& 01969 error=error) 01970 END IF 01971 01972 01973 ! *** finish initialization of the MOs *** 01974 CPPrecondition(ASSOCIATED(mos),cp_failure_level,routineP,error,failure) 01975 IF (.NOT.failure) THEN 01976 DO ispin=1,SIZE(mos) 01977 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff=mo_coeff,mo_coeff_b=mo_coeff_b,homo=homo) 01978 IF (.NOT.ASSOCIATED(mo_coeff)) THEN 01979 CALL init_mo_set(mos(ispin)%mo_set,& 01980 ao_mo_fm_pools(ispin)%pool,& 01981 name="qs_env"//TRIM(ADJUSTL(cp_to_string(qs_env%id_nr)))//& 01982 "%mo"//TRIM(ADJUSTL(cp_to_string(ispin))),& 01983 error=error) 01984 END IF 01985 IF(.not.ASSOCIATED(mo_coeff_b)) THEN 01986 CALL cp_fm_get_info(mos(ispin)%mo_set%mo_coeff,ncol_global=nmo,error=error) 01987 ! 01988 ! check if we need the mo_coeff_b 01989 mos(ispin)%mo_set%use_mo_coeff_b = qs_env%scf_control%use_ot .OR.& 01990 scf_env%method .EQ. ot_diag_method_nr .OR. & 01991 scf_env%method .EQ. block_davidson_diag_method_nr 01992 CALL cp_create_bl_distribution (col_dist, col_blk_size, nmo, &! fm->dbcsr 01993 dbcsr_mp_npcols(dbcsr_distribution_mp(cp_dbcsr_distribution(matrix_s(1)%matrix))))! fm->dbcsr 01994 CALL dbcsr_distribution_new (dist, dbcsr_distribution_mp (cp_dbcsr_distribution(matrix_s(1)%matrix)),&! fm->dbcsr 01995 dbcsr_distribution_row_dist(cp_dbcsr_distribution(matrix_s(1)%matrix)), col_dist)! fm->dbcsr 01996 CALL cp_dbcsr_init_p(mos(ispin)%mo_set%mo_coeff_b,error=error) 01997 CALL cp_dbcsr_create(mos(ispin)%mo_set%mo_coeff_b, "mo_coeff_b", dist, dbcsr_type_no_symmetry,&! fm->dbcsr 01998 cp_dbcsr_row_block_sizes(matrix_s(1)%matrix), col_blk_size, 0, 0, dbcsr_type_real_default,&!fm->dbcsr 01999 error=error)! fm->dbcsr 02000 CALL cp_dbcsr_distribution_release (dist)! fm->dbcsr 02001 CALL array_release (col_blk_size)! fm->dbcsr 02002 CALL array_release (col_dist)! fm->dbcsr 02003 ENDIF 02004 ! very first tests for xas 02005 IF(do_xas_tp .AND. ispin==1) THEN 02006 CALL set_mo_occupation(mo_set=mos(ispin)%mo_set,smear=xas_env%smear,& 02007 xas_env=xas_env,error=error) 02008 CALL set_mo_set(mos(ispin)%mo_set,& 02009 uniform_occupation=.FALSE.,error=error) 02010 END IF 02011 END DO 02012 END IF 02013 02014 ! *** finish initialization of the MOs for ADMM *** 02015 IF(dft_control%do_admm) THEN 02016 CPPrecondition(ASSOCIATED(mos_aux_fit),cp_failure_level,routineP,error,failure) 02017 IF (.NOT.failure) THEN 02018 DO ispin=1,SIZE(mos_aux_fit) 02019 CALL get_mo_set(mos_aux_fit(ispin)%mo_set,mo_coeff=mo_coeff_aux_fit,homo=homo) 02020 IF (.NOT.ASSOCIATED(mo_coeff_aux_fit)) THEN 02021 CALL init_mo_set(mos_aux_fit(ispin)%mo_set,& 02022 ao_mo_fm_pools_aux_fit(ispin)%pool,& 02023 name="qs_env"//TRIM(ADJUSTL(cp_to_string(qs_env%id_nr)))//& 02024 "%mo_aux_fit"//TRIM(ADJUSTL(cp_to_string(ispin))),& 02025 error=error) 02026 END IF 02027 IF(do_xas_tp .AND. ispin==1) THEN 02028 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 02029 END IF 02030 END DO 02031 END IF 02032 END IF 02033 02034 ! *** get the mo_derivs OK if needed *** 02035 IF (qs_env%requires_mo_derivs) THEN 02036 CALL get_qs_env(qs_env,mo_derivs=mo_derivs,error=error) 02037 IF (.NOT.ASSOCIATED(mo_derivs)) THEN 02038 IF (dft_control%restricted) THEN ! right now, there might be more mos than needed derivs 02039 ALLOCATE(mo_derivs(1)) 02040 CALL get_mo_set(mos(1)%mo_set,mo_coeff_b=mo_coeff_b) 02041 NULLIFY(mo_derivs(1)%matrix) 02042 CALL cp_dbcsr_init_p(mo_derivs(1)%matrix,error=error) 02043 CALL cp_dbcsr_create(mo_derivs(1)%matrix, "mo_derivs",& 02044 cp_dbcsr_distribution(mo_coeff_b), dbcsr_type_no_symmetry,& 02045 cp_dbcsr_row_block_sizes(mo_coeff_b), cp_dbcsr_col_block_sizes(mo_coeff_b),& 02046 0, 0, dbcsr_type_real_default, error=error) 02047 ELSE 02048 ALLOCATE(mo_derivs(dft_control%nspins)) 02049 DO ispin=1,dft_control%nspins 02050 CALL get_mo_set(mos(ispin)%mo_set,mo_coeff_b=mo_coeff_b) 02051 NULLIFY(mo_derivs(ispin)%matrix) 02052 CALL cp_dbcsr_init_p(mo_derivs(ispin)%matrix,error=error) 02053 CALL cp_dbcsr_create(mo_derivs(ispin)%matrix, "mo_derivs",& 02054 cp_dbcsr_distribution(mo_coeff_b), dbcsr_type_no_symmetry,& 02055 cp_dbcsr_row_block_sizes(mo_coeff_b), cp_dbcsr_col_block_sizes(mo_coeff_b),& 02056 0, 0, dbcsr_type_real_default, error=error) 02057 ENDDO 02058 ENDIF 02059 CALL set_qs_env(qs_env,mo_derivs=mo_derivs,error=error) 02060 ENDIF 02061 02062 ELSE 02063 ! nothing should be done 02064 ENDIF 02065 02066 ! *** get the mo_derivs OK if needed *** 02067 IF( dft_control%do_admm) THEN 02068 IF (qs_env%requires_mo_derivs) THEN 02069 CALL get_qs_env(qs_env,mo_derivs_aux_fit=mo_derivs_aux_fit,error=error) 02070 IF (.NOT.ASSOCIATED(mo_derivs_aux_fit)) THEN 02071 IF (dft_control%restricted) THEN ! right now, there might be more mos than needed derivs 02072 ALLOCATE(mo_derivs_aux_fit(1)) 02073 CALL get_mo_set(mos_aux_fit(1)%mo_set,mo_coeff=mo_coeff_aux_fit) 02074 CALL cp_fm_create(mo_derivs_aux_fit(1)%matrix,mo_coeff_aux_fit%matrix_struct,error=error) 02075 ELSE 02076 ALLOCATE(mo_derivs_aux_fit(dft_control%nspins)) 02077 DO ispin=1,dft_control%nspins 02078 CALL get_mo_set(mos_aux_fit(ispin)%mo_set,mo_coeff=mo_coeff_aux_fit) 02079 CALL cp_fm_create(mo_derivs_aux_fit(ispin)%matrix,mo_coeff_aux_fit%matrix_struct,error=error) 02080 ENDDO 02081 ENDIF 02082 CALL set_qs_env(qs_env,mo_derivs_aux_fit=mo_derivs_aux_fit,error=error) 02083 ENDIF 02084 ELSE 02085 ! nothing should be done 02086 ENDIF 02087 END IF 02088 02089 ! *** Allocate the distributed SCF matrices *** 02090 IF ((.NOT.ASSOCIATED(scf_env%scf_work1)).OR.& 02091 (.NOT.ASSOCIATED(scf_env%scf_work2)).OR.& 02092 (.NOT.ASSOCIATED(scf_env%ortho)).OR.& 02093 (.NOT.ASSOCIATED(scf_env%s_half))) THEN 02094 02095 ao_mo_fmstruct => fm_pool_get_el_struct(ao_mo_fm_pools(1)%pool,& 02096 error=error) 02097 CALL cp_fm_struct_get(ao_mo_fmstruct, nrow_block=nrow_block,& 02098 error=error) 02099 CALL get_mo_set(mos(1)%mo_set,nao=nao) 02100 CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct,& 02101 nrow_block=nrow_block,ncol_block=nrow_block,& 02102 nrow_global=nao, ncol_global=nao,& 02103 template_fmstruct=ao_mo_fmstruct, error=error) 02104 02105 IF((scf_env%method /= ot_method_nr) .AND. & 02106 (scf_env%method /= block_davidson_diag_method_nr) ) THEN 02107 IF (.NOT.ASSOCIATED(scf_env%scf_work1)) THEN 02108 ALLOCATE(scf_env%scf_work1(dft_control%nspins), stat=stat) 02109 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02110 DO ispin=1,SIZE(scf_env%scf_work1) 02111 NULLIFY(scf_env%scf_work1(ispin)%matrix) 02112 CALL cp_fm_create(scf_env%scf_work1(ispin)%matrix,& 02113 matrix_struct=ao_ao_fmstruct,& 02114 name="SCF"//TRIM(ADJUSTL(cp_to_string(scf_env%id_nr)))//& 02115 "WORK_MATRIX-1-"//TRIM(ADJUSTL(cp_to_string(ispin))),& 02116 error=error) 02117 ENDDO 02118 END IF 02119 IF ((.NOT.ASSOCIATED(scf_env%ortho)).AND.& 02120 (scf_env%method /= ot_diag_method_nr).AND.& 02121 (scf_env%method /= special_diag_method_nr)) THEN 02122 CALL cp_fm_create(scf_env%ortho,& 02123 matrix_struct=ao_ao_fmstruct,& 02124 name="SCF"//TRIM(ADJUSTL(cp_to_string(scf_env%id_nr)))//& 02125 "ORTHO_MATRIX",& 02126 error=error) 02127 IF(scf_env%cholesky_method==cholesky_inverse) THEN 02128 CALL cp_fm_create(scf_env%ortho_m1, matrix_struct=ao_ao_fmstruct,& 02129 name="SCF"//TRIM(ADJUSTL(cp_to_string(scf_env%id_nr)))//& 02130 "ORTHO_MATRIX-1",& 02131 error=error) 02132 END IF 02133 END IF 02134 IF ((.NOT.ASSOCIATED(scf_env%scf_work2))) THEN 02135 CALL cp_fm_create(scf_env%scf_work2,& 02136 matrix_struct=ao_ao_fmstruct,& 02137 name="SCF"//TRIM(ADJUSTL(cp_to_string(scf_env%id_nr)))//& 02138 "WORK_MATRIX-2",& 02139 error=error) 02140 END IF 02141 END IF 02142 02143 IF (dft_control%dft_plus_u) THEN 02144 IF (dft_control%plus_u_method_id == plus_u_lowdin) THEN 02145 IF (.NOT.ASSOCIATED(scf_env%scf_work2)) THEN 02146 CALL cp_fm_create(scf_env%scf_work2,& 02147 matrix_struct=ao_ao_fmstruct,& 02148 name="SCF"//TRIM(ADJUSTL(cp_to_string(scf_env%id_nr)))//& 02149 "WORK_MATRIX-2",& 02150 error=error) 02151 END IF 02152 IF (.NOT.ASSOCIATED(scf_env%s_half)) THEN 02153 CALL cp_fm_create(scf_env%s_half,& 02154 matrix_struct=ao_ao_fmstruct,& 02155 name="S**(1/2) MATRIX",& 02156 error=error) 02157 END IF 02158 END IF 02159 END IF 02160 CALL cp_fm_struct_release(ao_ao_fmstruct,error=error) 02161 END IF 02162 02163 ! *** Allocate matrix_ks and put it in the QS environment *** 02164 02165 IF (.NOT.ASSOCIATED(matrix_ks)) THEN 02166 CALL cp_dbcsr_allocate_matrix_set(matrix_ks, dft_control%nspins, error) 02167 DO ispin=1,dft_control%nspins 02168 ALLOCATE(matrix_ks(ispin)%matrix) 02169 CALL cp_dbcsr_init(matrix_ks(ispin)%matrix,error=error) 02170 IF (dft_control%nspins > 1) THEN 02171 IF (ispin == 1) THEN 02172 headline="KOHN-SHAM MATRIX FOR ALPHA SPIN" 02173 ELSE 02174 headline="KOHN-SHAM MATRIX FOR BETA SPIN" 02175 END IF 02176 ELSE 02177 headline="KOHN-SHAM MATRIX" 02178 END IF 02179 CALL cp_dbcsr_create(matrix=matrix_ks(ispin)%matrix,& 02180 name=TRIM(headline)//" (SCF ENV "//TRIM(ADJUSTL(cp_to_string(scf_env%id_nr)))//")",& 02181 dist=cp_dbcsr_distribution(matrix_s(1)%matrix), matrix_type=dbcsr_type_symmetric,& 02182 row_blk_size=cp_dbcsr_row_block_sizes(matrix_s(1)%matrix),& 02183 col_blk_size=cp_dbcsr_col_block_sizes(matrix_s(1)%matrix),& 02184 nblks=0, nze=0, error=error) 02185 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks(ispin)%matrix,qs_env%sab_orb,error=error) 02186 CALL cp_dbcsr_set(matrix_ks(ispin)%matrix,0.0_dp,error=error) 02187 ENDDO 02188 02189 CALL set_qs_env(qs_env=qs_env,& 02190 matrix_ks=matrix_ks,error=error) 02191 02192 END IF 02193 02194 ! *** Allocate matrix_ks_aux_fit if requested and put it in the QS environment *** 02195 02196 IF( dft_control%do_admm) THEN 02197 IF (.NOT.ASSOCIATED(matrix_ks_aux_fit)) THEN 02198 CALL cp_dbcsr_allocate_matrix_set(matrix_ks_aux_fit, dft_control%nspins, error) 02199 DO ispin=1,dft_control%nspins 02200 ALLOCATE(matrix_ks_aux_fit(ispin)%matrix) 02201 CALL cp_dbcsr_init(matrix_ks_aux_fit(ispin)%matrix,error=error) 02202 02203 CALL cp_dbcsr_create(matrix=matrix_ks_aux_fit(ispin)%matrix, & 02204 name="SCF"//TRIM(ADJUSTL(cp_to_string(scf_env%id_nr)))//& 02205 "KOHN-SHAM_MATRIX for ADMM", & 02206 dist=cp_dbcsr_distribution(qs_env%matrix_s_aux_fit(1)%matrix), matrix_type=dbcsr_type_symmetric, & 02207 row_blk_size=cp_dbcsr_row_block_sizes(qs_env%matrix_s_aux_fit(1)%matrix),& 02208 col_blk_size=cp_dbcsr_col_block_sizes(qs_env%matrix_s_aux_fit(1)%matrix),& 02209 nblks=0, nze=0, error=error) 02210 CALL cp_dbcsr_alloc_block_from_nbl(matrix_ks_aux_fit(ispin)%matrix,qs_env%sab_aux_fit,& 02211 error=error) 02212 CALL cp_dbcsr_set(matrix_ks_aux_fit(ispin)%matrix,0.0_dp,error=error) 02213 ENDDO 02214 02215 CALL set_qs_env(qs_env=qs_env,& 02216 matrix_ks_aux_fit=matrix_ks_aux_fit,error=error) 02217 END IF 02218 END IF 02219 02220 ! *** allocate the mixing storage 02221 IF(do_xas_tp) THEN 02222 mixing_method = xas_env%mixing_method 02223 mixing_store => xas_env%mixing_store 02224 ELSE 02225 mixing_method = scf_env%mixing_method 02226 mixing_store => scf_env%mixing_store 02227 END IF 02228 IF (mixing_method>0) THEN 02229 CALL mixing_allocate(qs_env,mixing_method,mixing_store,dft_control%nspins,error) 02230 ELSE 02231 NULLIFY(scf_env%p_mix_new) 02232 END IF 02233 02234 ! *** allocate the ks env ** 02235 IF (.NOT.ASSOCIATED(ks_env)) THEN 02236 CALL qs_ks_create(ks_env,qs_env=qs_env,error=error) 02237 CALL set_qs_env(qs_env, ks_env=ks_env,error=error) 02238 CALL qs_ks_release(ks_env,error=error) 02239 END IF 02240 02241 ! If there is an outer scf loop allocate the space for the variables 02242 CALL get_qs_env(qs_env=qs_env,scf_control=scf_control, dft_control=dft_control,error=error) 02243 IF (scf_control%outer_scf%have_scf) THEN 02244 nhistory = scf_control%outer_scf%max_scf + 1 02245 nvariables = outer_loop_variables_count(scf_control,error) 02246 ALLOCATE(scf_env%outer_scf%variables(nvariables,nhistory),STAT=stat) 02247 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02248 ALLOCATE(scf_env%outer_scf%count(nhistory),STAT=stat) 02249 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02250 scf_env%outer_scf%count=0 02251 ALLOCATE(scf_env%outer_scf%gradient(nvariables,nhistory),STAT=stat) 02252 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02253 ALLOCATE(scf_env%outer_scf%energy(nhistory),STAT=stat) 02254 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02255 ENDIF 02256 02257 IF(scf_env%method==block_krylov_diag_method_nr ) THEN 02258 CALL krylov_space_allocate(scf_env%krylov_space, scf_control, mos, error=error) 02259 END IF 02260 IF(scf_env%do_diag_sub ) THEN 02261 CALL diag_subspace_allocate(scf_env%subspace_env, qs_env, mos, error=error) 02262 END IF 02263 IF(scf_env%method==block_davidson_diag_method_nr ) THEN 02264 DO ispin=1,dft_control%nspins 02265 CALL block_davidson_allocate(scf_env%block_davidson_env(ispin), mos(ispin)%mo_set, error=error) 02266 END DO 02267 END IF 02268 02269 CALL timestop(handle) 02270 02271 END SUBROUTINE scf_env_check_i_alloc 02272 02273 ! ***************************************************************************** 02283 SUBROUTINE scf_env_initial_rho_setup(scf_env, qs_env, scf_section, error) 02284 TYPE(qs_scf_env_type), POINTER :: scf_env 02285 TYPE(qs_environment_type), POINTER :: qs_env 02286 TYPE(section_vals_type), POINTER :: scf_section 02287 TYPE(cp_error_type), INTENT(inout) :: error 02288 02289 CHARACTER(len=*), PARAMETER :: routineN = 'scf_env_initial_rho_setup', 02290 routineP = moduleN//':'//routineN 02291 02292 INTEGER :: extrapolation_method_nr, 02293 handle, ispin, nmo, 02294 output_unit 02295 LOGICAL :: failure, orthogonal_wf 02296 TYPE(cp_fm_type), POINTER :: mo_coeff 02297 TYPE(cp_logger_type), POINTER :: logger 02298 TYPE(rho_atom_type), DIMENSION(:), 02299 POINTER :: rho_atom 02300 02301 CALL timeset(routineN,handle) 02302 failure=.FALSE. 02303 NULLIFY(mo_coeff ) 02304 logger => cp_error_get_logger(error) 02305 CPPrecondition(ASSOCIATED(scf_env),cp_failure_level,routineP,error,failure) 02306 CPPrecondition(scf_env%ref_count>0,cp_failure_level,routineP,error,failure) 02307 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 02308 02309 IF (.NOT.failure) THEN 02310 extrapolation_method_nr=wfi_use_guess_method_nr 02311 IF (ASSOCIATED(qs_env%wf_history)) THEN 02312 CALL wfi_extrapolate(qs_env%wf_history, & 02313 qs_env=qs_env, dt=1.0_dp, & 02314 extrapolation_method_nr=extrapolation_method_nr,& 02315 orthogonal_wf=orthogonal_wf, error=error) 02316 ! wfi_use_guess_method_nr the wavefunctions are not yet initialized 02317 IF ((.NOT.orthogonal_wf).AND.& 02318 (scf_env%method == ot_method_nr).AND.& 02319 (.NOT.(extrapolation_method_nr == wfi_use_guess_method_nr))) THEN 02320 DO ispin=1,SIZE(qs_env%mos) 02321 CALL get_mo_set(qs_env%mos(ispin)%mo_set, & 02322 mo_coeff=mo_coeff, nmo=nmo) 02323 CALL qs_env_reorthogonalize_vectors(qs_env, & 02324 v_matrix=mo_coeff, n_col=nmo,& 02325 error=error) 02326 CALL set_mo_occupation(mo_set=qs_env%mos(ispin)%mo_set, & 02327 smear=qs_env%scf_control%smear, error=error) 02328 END DO 02329 END IF 02330 END IF 02331 output_unit=cp_print_key_unit_nr(logger,scf_section,"PRINT%PROGRAM_RUN_INFO",& 02332 extension=".scfLog",error=error) 02333 IF (output_unit>0) THEN 02334 WRITE (UNIT=output_unit,FMT="(/,T2,A)")& 02335 "Extrapolation method: "//& 02336 TRIM(wfi_get_method_label(extrapolation_method_nr,error=error)) 02337 END IF 02338 02339 CALL cp_print_key_finished_output(output_unit,logger,scf_section,& 02340 "PRINT%PROGRAM_RUN_INFO", error=error) 02341 IF (extrapolation_method_nr==wfi_use_guess_method_nr) THEN 02342 CALL calculate_first_density_matrix(scf_env=scf_env,qs_env=qs_env,error=error) 02343 IF (.NOT.(qs_env%scf_control%density_guess==densities_guess)) THEN 02344 CALL qs_rho_update_rho(qs_env%rho,qs_env=qs_env, error=error) 02345 CALL qs_ks_did_change(qs_env%ks_env,rho_changed=.TRUE.,& 02346 error=error) 02347 END IF 02348 END IF 02349 02350 ! Some preparation for the mixing 02351 IF(qs_env%scf_env%mixing_method>1) THEN 02352 IF(qs_env % dft_control % qs_control%gapw) THEN 02353 CALL get_qs_env(qs_env=qs_env,& 02354 rho_atom_set=rho_atom,error=error) 02355 CALL mixing_init(qs_env%scf_env%mixing_method,qs_env%rho,qs_env%scf_env%mixing_store,& 02356 qs_env%para_env,rho_atom=rho_atom,error=error) 02357 ELSE 02358 CALL mixing_init(qs_env%scf_env%mixing_method,qs_env%rho,qs_env%scf_env%mixing_store,& 02359 qs_env%para_env,error=error) 02360 END IF 02361 END IF 02362 02363 ! *** SCP 02364 IF ( qs_env % dft_control % scp ) THEN 02365 CALL update_rhoscp ( qs_env = qs_env, error = error ) 02366 CALL qs_ks_scp_did_change ( qs_env=qs_env, rho_changed = .TRUE., error = error ) 02367 END IF 02368 ! *** SCP 02369 02370 DO ispin=1,SIZE(qs_env%mos)!fm->dbcsr 02371 IF(qs_env%mos(ispin)%mo_set%use_mo_coeff_b) THEN 02372 CALL copy_fm_to_dbcsr(qs_env%mos(ispin)%mo_set%mo_coeff,& 02373 qs_env%mos(ispin)%mo_set%mo_coeff_b,error=error)!fm->dbcsr 02374 ENDIF 02375 ENDDO!fm->dbcsr 02376 02377 END IF 02378 02379 CALL timestop(handle) 02380 02381 END SUBROUTINE scf_env_initial_rho_setup 02382 02383 ! ***************************************************************************** 02392 SUBROUTINE scf_env_cleanup(scf_env,qs_env,error) 02393 TYPE(qs_scf_env_type), POINTER :: scf_env 02394 TYPE(qs_environment_type), POINTER :: qs_env 02395 TYPE(cp_error_type), INTENT(inout) :: error 02396 02397 CHARACTER(len=*), PARAMETER :: routineN = 'scf_env_cleanup', 02398 routineP = moduleN//':'//routineN 02399 02400 INTEGER :: handle, ispin, stat 02401 LOGICAL :: failure 02402 02403 CALL timeset(routineN,handle) 02404 02405 failure=.FALSE. 02406 02407 CPPrecondition(ASSOCIATED(scf_env),cp_failure_level,routineP,error,failure) 02408 CPPrecondition(scf_env%ref_count>0,cp_failure_level,routineP,error,failure) 02409 CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure) 02410 IF (.NOT.failure) THEN 02411 02412 ! *** Release SCF work storage *** 02413 02414 IF (ASSOCIATED(scf_env%scf_work1)) THEN 02415 DO ispin=1,SIZE(scf_env%scf_work1) 02416 CALL cp_fm_release(scf_env%scf_work1(ispin)%matrix,error=error) 02417 ENDDO 02418 DEALLOCATE(scf_env%scf_work1) 02419 ENDIF 02420 IF (ASSOCIATED(scf_env%scf_work2)) CALL cp_fm_release(scf_env%scf_work2,error) 02421 IF (ASSOCIATED(scf_env%ortho)) CALL cp_fm_release(scf_env%ortho,error=error) 02422 IF (ASSOCIATED(scf_env%ortho_m1)) CALL cp_fm_release(scf_env%ortho_m1,error=error) 02423 02424 IF (ASSOCIATED(scf_env%p_mix_new)) THEN 02425 CALL cp_dbcsr_deallocate_matrix_set(scf_env%p_mix_new,error=error) 02426 END IF 02427 02428 IF (ASSOCIATED(scf_env%p_delta)) THEN 02429 CALL cp_dbcsr_deallocate_matrix_set(scf_env%p_delta,error=error) 02430 END IF 02431 02432 ! *** method dependent cleanup 02433 SELECT CASE(scf_env%method) 02434 CASE(ot_method_nr) 02435 DO ispin=1,SIZE(scf_env%ot_preconditioner) 02436 CALL destroy_preconditioner(scf_env%ot_preconditioner(ispin)%preconditioner,error=error) 02437 DEALLOCATE(scf_env%ot_preconditioner(ispin)%preconditioner) 02438 ENDDO 02439 DEALLOCATE(scf_env%ot_preconditioner,stat=stat) 02440 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02441 CASE(ot_diag_method_nr) 02442 ! 02443 CASE(general_diag_method_nr) 02444 ! 02445 CASE(special_diag_method_nr) 02446 ! 02447 CASE(block_krylov_diag_method_nr) 02448 CASE(block_davidson_diag_method_nr) 02449 CALL block_davidson_deallocate(scf_env%block_davidson_env,error=error) 02450 CASE DEFAULT 02451 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 02452 routineP,"unknown scf method method:"//& 02453 cp_to_string(scf_env%method),error,failure) 02454 END SELECT 02455 02456 IF (ASSOCIATED(scf_env%outer_scf%variables)) THEN 02457 DEALLOCATE(scf_env%outer_scf%variables,stat=stat) 02458 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02459 ENDIF 02460 IF (ASSOCIATED(scf_env%outer_scf%count)) THEN 02461 DEALLOCATE(scf_env%outer_scf%count,stat=stat) 02462 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02463 ENDIF 02464 IF (ASSOCIATED(scf_env%outer_scf%gradient)) THEN 02465 DEALLOCATE(scf_env%outer_scf%gradient,stat=stat) 02466 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02467 ENDIF 02468 IF (ASSOCIATED(scf_env%outer_scf%energy)) THEN 02469 DEALLOCATE(scf_env%outer_scf%energy,stat=stat) 02470 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02471 ENDIF 02472 02473 IF (scf_env%do_diag_sub) THEN 02474 IF(ASSOCIATED(scf_env%subspace_env%p_matrix_store)) THEN 02475 CALL cp_dbcsr_deallocate_matrix_set(scf_env%subspace_env%p_matrix_store,error=error) 02476 END IF 02477 IF(ASSOCIATED(scf_env%subspace_env%p_matrix_mix)) THEN 02478 CALL cp_dbcsr_deallocate_matrix_set(scf_env%subspace_env%p_matrix_mix,error=error) 02479 END IF 02480 END IF 02481 02482 END IF 02483 CALL timestop(handle) 02484 02485 END SUBROUTINE scf_env_cleanup 02486 02487 ! ***************************************************************************** 02492 SUBROUTINE cleanup_scf_loop(scf_env,error) 02493 TYPE(qs_scf_env_type), POINTER :: scf_env 02494 TYPE(cp_error_type), INTENT(inout) :: error 02495 02496 CHARACTER(len=*), PARAMETER :: routineN = 'cleanup_scf_loop', 02497 routineP = moduleN//':'//routineN 02498 02499 INTEGER :: handle, ispin, stat 02500 LOGICAL :: failure 02501 02502 CALL timeset(routineN,handle) 02503 02504 failure=.FALSE. 02505 02506 CPPrecondition(ASSOCIATED(scf_env),cp_failure_level,routineP,error,failure) 02507 CPPrecondition(scf_env%ref_count>0,cp_failure_level,routineP,error,failure) 02508 IF (.NOT.failure) THEN 02509 ! *** method dependent cleanup 02510 SELECT CASE(scf_env%method) 02511 CASE(ot_method_nr) 02512 DO ispin=1,SIZE(scf_env%qs_ot_env) 02513 CALL ot_scf_destroy(scf_env%qs_ot_env(ispin),error=error) 02514 ENDDO 02515 DEALLOCATE(scf_env%qs_ot_env,stat=stat) 02516 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02517 CASE(ot_diag_method_nr) 02518 ! 02519 CASE(general_diag_method_nr) 02520 ! 02521 CASE(special_diag_method_nr) 02522 ! 02523 CASE(block_krylov_diag_method_nr,block_davidson_diag_method_nr) 02524 CASE DEFAULT 02525 CALL cp_assert(.FALSE.,cp_failure_level,cp_assertion_failed,& 02526 routineP,"unknown scf method method:"//& 02527 cp_to_string(scf_env%method),error,failure) 02528 END SELECT 02529 02530 END IF 02531 02532 CALL timestop(handle) 02533 02534 END SUBROUTINE cleanup_scf_loop 02535 02536 END MODULE qs_scf
1.7.3