From d10c76ea8aa3f3a98cc2d31fc41c81d821ac3bc4 Mon Sep 17 00:00:00 2001 From: Dmitry Ryndyk Date: Mon, 30 Mar 2026 10:38:21 +0200 Subject: [PATCH] NEGF: restart for NEGF SCF (#5023) --- src/common/cp_iter_types.F | 10 +- src/input_cp2k_negf.F | 34 +++- src/negf_control_types.F | 7 +- src/negf_env_types.F | 16 +- src/negf_io.F | 72 ++++--- src/negf_methods.F | 390 +++++++++++++++++++++++++++++++++---- 6 files changed, 451 insertions(+), 78 deletions(-) diff --git a/src/common/cp_iter_types.F b/src/common/cp_iter_types.F index 83df883e0c..f8b9426f3e 100644 --- a/src/common/cp_iter_types.F +++ b/src/common/cp_iter_types.F @@ -27,7 +27,7 @@ MODULE cp_iter_types LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE. ! When adding a new iteration level PLEASE update the following list with the proper name! - CHARACTER(LEN=default_path_length), PARAMETER, PUBLIC, DIMENSION(18) :: each_possible_labels = [ & + CHARACTER(LEN=default_path_length), PARAMETER, PUBLIC, DIMENSION(19) :: each_possible_labels = [ & "__ROOT__ ", & "JUST_ENERGY ", & "POWELL_OPT ", & @@ -45,9 +45,10 @@ MODULE cp_iter_types "REPLICA_EVAL ", & "BSSE ", & "SHELL_OPT ", & - "TDDFT_SCF "] + "TDDFT_SCF ", & + "NEGF_SCF "] - CHARACTER(LEN=default_path_length), PARAMETER, PUBLIC, DIMENSION(18) :: each_desc_labels = [ & + CHARACTER(LEN=default_path_length), PARAMETER, PUBLIC, DIMENSION(19) :: each_desc_labels = [ & "Iteration level for __ROOT__ (fictitious iteration level) ", & "Iteration level for an ENERGY/ENERGY_FORCE calculation. ", & "Iteration level for POWELL based optimization steps. ", & @@ -65,7 +66,8 @@ MODULE cp_iter_types "Iteration level for the evaluation of the Replica Environment ", & "Iteration level for the Basis Set Superposition Error (BSSE) calculation ", & "Iteration level for the Shell-Core distances optimization steps ", & - "Iteration level for the Time-Dependent Density Functional Theory SCF steps. "] + "Iteration level for the Time-Dependent Density Functional Theory SCF steps. ", & + "Iteration level for the NEGF SCF steps. "] ! ************************************************************************************************** !> \brief contains the information about the current state of the program diff --git a/src/input_cp2k_negf.F b/src/input_cp2k_negf.F index 73a861b58e..849a7830e1 100644 --- a/src/input_cp2k_negf.F +++ b/src/input_cp2k_negf.F @@ -209,12 +209,6 @@ CONTAINS n_var=1, type_of_var=integer_t, default_i_val=30) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="UPDATE_HS", & - description="Change the atomic Hamiltonian during NEGF self-consistent cycle."// & - " If .FALSE., only the density matrix is updated.", & - default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) - CALL section_add_keyword(section, keyword) - CALL keyword_release(keyword) ! PRINT subsection CALL section_create(subsection, __LOCATION__, "PRINT", "Printing of information during the NEGF.", & @@ -237,14 +231,40 @@ CONTAINS description="Controls the dumping of the restart files during NEGF SCF."// & " By default keeps a short history of three restarts.", & print_level=low_print_level, common_iter_levels=3, & - each_iter_names=s2a("NEGF_SCF"), each_iter_values=[20], & + each_iter_names=s2a("NEGF_SCF"), each_iter_values=[1], & add_last=add_last_numeric, filename="RESTART") + CALL keyword_create(keyword, __LOCATION__, name="BACKUP_COPIES", & + description="Specifies the maximum number of backup copies.", & + usage="BACKUP_COPIES {int}", & + default_i_val=3) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) CALL section_add_subsection(subsection, print_key) CALL section_release(print_key) CALL section_add_subsection(section, subsection) CALL section_release(subsection) + ! SCF subsection + CALL section_create(subsection, __LOCATION__, "SCF", "Additional parameters for an NEGF SCF run.", & + repeats=.FALSE.) + + CALL keyword_create(keyword, __LOCATION__, name="UPDATE_HS", & + description="Change the atomic Hamiltonian during NEGF self-consistent cycle."// & + " If .FALSE., only the density matrix is updated.", & + default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(subsection, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="RESTART_SCF", & + description="Restart SCF from files saved as defined in NEGF%PRINT%RESTART", & + default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(subsection, keyword) + CALL keyword_release(keyword) + + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + END SUBROUTINE create_negf_section ! ************************************************************************************************** diff --git a/src/negf_control_types.F b/src/negf_control_types.F index d287ee76e1..8d9c718239 100644 --- a/src/negf_control_types.F +++ b/src/negf_control_types.F @@ -88,7 +88,9 @@ MODULE negf_control_types !> to read/write H and S from/to file LOGICAL :: read_write_HS = .FALSE. !> to update the atomic Hamiltonian during NEGF self-consistent cycle - LOGICAL :: update_HS = .FALSE. + LOGICAL :: update_HS = .TRUE. + !> if NEGF SCF id restart from saved files + LOGICAL :: restart_scf = .TRUE. !> if dft of entire system is done LOGICAL :: is_dft_entire = .FALSE. !> if restart from files is really done @@ -374,7 +376,8 @@ CONTAINS CALL section_vals_val_get(negf_section, "V_SHIFT_OFFSET", r_val=negf_control%v_shift_offset) CALL section_vals_val_get(negf_section, "V_SHIFT_MAX_ITERS", i_val=negf_control%v_shift_maxiters) - CALL section_vals_val_get(negf_section, "UPDATE_HS", l_val=negf_control%update_HS) + CALL section_vals_val_get(negf_section, "SCF%UPDATE_HS", l_val=negf_control%update_HS) + CALL section_vals_val_get(negf_section, "SCF%RESTART_SCF", l_val=negf_control%restart_scf) ! check consistency IF (negf_control%eta < 0.0_dp) THEN diff --git a/src/negf_env_types.F b/src/negf_env_types.F index e370b181e5..b9fb915dee 100644 --- a/src/negf_env_types.F +++ b/src/negf_env_types.F @@ -334,7 +334,10 @@ CONTAINS CALL negf_env_scatt_read_write_hs(force_env, para_env, negf_env, sub_env, negf_control, negf_section, log_unit, & is_dft_entire=is_dft_entire) ELSE - IF (.NOT. is_dft_entire) CALL qs_energies(qs_env, consistent_energies=.FALSE., calc_forces=.FALSE.) + IF (.NOT. is_dft_entire) THEN + CALL qs_energies(qs_env, consistent_energies=.FALSE., calc_forces=.FALSE.) + is_dft_entire = .TRUE. + END IF ! extract device-related matrix blocks CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env) END IF @@ -549,8 +552,10 @@ CONTAINS IF (exist_all) THEN negf_control%contacts(icontact)%is_restart = .TRUE. - - IF (log_unit > 0) WRITE (log_unit, '(/,T2,A)') "All restart files exist." + IF (log_unit > 0) THEN + WRITE (log_unit, '(/,T2,A)') "User requested to read the Hamiltonian and overlap matrices from files." + WRITE (log_unit, '(T2,A)') "All restart files exist." + END IF ! ++ create matrices: s_00, s_01, h_00, h_01 IF (para_env%is_source()) THEN @@ -1112,7 +1117,10 @@ CONTAINS negf_control%is_restart = .TRUE. - IF (log_unit > 0) WRITE (log_unit, '(/,T2,A)') "All restart files exist." + IF (log_unit > 0) THEN + WRITE (log_unit, '(/,T2,A)') "User requested to read the Hamiltonian and overlap matrices from files." + WRITE (log_unit, '(T2,A)') "All restart files exist." + END IF ! ++ create matrices: s_s, s_sc, h_s, h_sc IF (para_env%is_source()) THEN diff --git a/src/negf_io.F b/src/negf_io.F index 639f4f98b3..88061d2ccf 100644 --- a/src/negf_io.F +++ b/src/negf_io.F @@ -15,6 +15,7 @@ MODULE negf_io open_file USE cp_log_handling, ONLY: cp_logger_generate_filename,& cp_logger_type + USE cp_output_handling, ONLY: cp_iter_string USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type,& section_vals_val_get @@ -51,20 +52,22 @@ CONTAINS !> \param s ... !> \param hc ... !> \param sc ... +!> \param h_scf ... !> \par History !> * 12.2025 created [Dmitry Ryndyk] ! ************************************************************************************************** - SUBROUTINE negf_restart_file_name(filename, exist, negf_section, logger, icontact, ispin, h00, h01, s00, s01, h, s, hc, sc) + SUBROUTINE negf_restart_file_name(filename, exist, negf_section, logger, icontact, ispin, h00, h01, & + s00, s01, h, s, hc, sc, h_scf) CHARACTER(LEN=default_path_length), INTENT(OUT) :: filename LOGICAL, INTENT(OUT) :: exist TYPE(section_vals_type), POINTER :: negf_section TYPE(cp_logger_type), POINTER :: logger INTEGER, INTENT(IN), OPTIONAL :: icontact, ispin - LOGICAL, INTENT(IN), OPTIONAL :: h00, h01, s00, s01, h, s, hc, sc + LOGICAL, INTENT(IN), OPTIONAL :: h00, h01, s00, s01, h, s, hc, sc, h_scf CHARACTER(len=default_string_length) :: middle_name, string1, string2 - LOGICAL :: my_h, my_h00, my_h01, my_hc, my_s, & - my_s00, my_s01, my_sc + LOGICAL :: my_h, my_h00, my_h01, my_h_scf, my_hc, & + my_s, my_s00, my_s01, my_sc TYPE(section_vals_type), POINTER :: contact_section, print_key my_h00 = .FALSE. @@ -83,6 +86,8 @@ CONTAINS IF (PRESENT(hc)) my_hc = hc my_sc = .FALSE. IF (PRESENT(sc)) my_sc = sc + my_h_scf = .FALSE. + IF (PRESENT(h_scf)) my_h_scf = h_scf exist = .FALSE. @@ -99,8 +104,8 @@ CONTAINS ELSE middle_name = "N"//TRIM(string1)//"-H00-S"//TRIM(string2) END IF - filename = negf_elecrodes_generate_filename(logger, print_key, middle_name=middle_name, & - extension=".hs", my_local=.FALSE.) + filename = negf_generate_filename(logger, print_key, middle_name=middle_name, & + extension=".hs", my_local=.FALSE.) END IF IF (my_h01) THEN @@ -109,20 +114,20 @@ CONTAINS ELSE middle_name = "N"//TRIM(string1)//"-H01-S"//TRIM(string2) END IF - filename = negf_elecrodes_generate_filename(logger, print_key, middle_name=middle_name, & - extension=".hs", my_local=.FALSE.) + filename = negf_generate_filename(logger, print_key, middle_name=middle_name, & + extension=".hs", my_local=.FALSE.) END IF IF (my_s00) THEN middle_name = "N"//TRIM(string1)//"-S00" - filename = negf_elecrodes_generate_filename(logger, print_key, middle_name=middle_name, & - extension=".hs", my_local=.FALSE.) + filename = negf_generate_filename(logger, print_key, middle_name=middle_name, & + extension=".hs", my_local=.FALSE.) END IF IF (my_s01) THEN middle_name = "N"//TRIM(string1)//"-S01" - filename = negf_elecrodes_generate_filename(logger, print_key, middle_name=middle_name, & - extension=".hs", my_local=.FALSE.) + filename = negf_generate_filename(logger, print_key, middle_name=middle_name, & + extension=".hs", my_local=.FALSE.) END IF ! try to read from the filename that is generated automatically from the printkey @@ -134,14 +139,14 @@ CONTAINS ELSE middle_name = "Hs-S"//TRIM(string2) END IF - filename = negf_elecrodes_generate_filename(logger, print_key, middle_name=middle_name, & - extension=".hs", my_local=.FALSE.) + filename = negf_generate_filename(logger, print_key, middle_name=middle_name, & + extension=".hs", my_local=.FALSE.) END IF IF (my_s) THEN middle_name = "Ss" - filename = negf_elecrodes_generate_filename(logger, print_key, middle_name=middle_name, & - extension=".hs", my_local=.FALSE.) + filename = negf_generate_filename(logger, print_key, middle_name=middle_name, & + extension=".hs", my_local=.FALSE.) END IF IF (my_hc) THEN @@ -150,14 +155,22 @@ CONTAINS ELSE middle_name = "Hsc-N"//TRIM(string1)//"-S"//TRIM(string2) END IF - filename = negf_elecrodes_generate_filename(logger, print_key, middle_name=middle_name, & - extension=".hs", my_local=.FALSE.) + filename = negf_generate_filename(logger, print_key, middle_name=middle_name, & + extension=".hs", my_local=.FALSE.) END IF IF (my_sc) THEN middle_name = "Ssc-N"//TRIM(string1) - filename = negf_elecrodes_generate_filename(logger, print_key, middle_name=middle_name, & - extension=".hs", my_local=.FALSE.) + filename = negf_generate_filename(logger, print_key, middle_name=middle_name, & + extension=".hs", my_local=.FALSE.) + END IF + + ! try to read from the filename that is generated automatically from the printkey + print_key => section_vals_get_subs_vals(negf_section, "PRINT%RESTART") + + IF (my_h_scf) THEN + filename = negf_generate_filename(logger, print_key, & + extension="", my_local=.FALSE., iter_string=.TRUE.) END IF INQUIRE (FILE=filename, exist=exist) @@ -175,17 +188,18 @@ CONTAINS !> \param extension extension to be applied to the filename (including the ".") !> \param my_local if the unit should be local to this task, or global to the !> program (defaults to false). +!> \param iter_string ... !> \return ... !> \par History !> * 12.2025 created [Dmitry Ryndyk] ! ************************************************************************************************** - FUNCTION negf_elecrodes_generate_filename(logger, print_key, middle_name, extension, & - my_local) RESULT(filename) + FUNCTION negf_generate_filename(logger, print_key, middle_name, extension, & + my_local, iter_string) RESULT(filename) TYPE(cp_logger_type), POINTER :: logger TYPE(section_vals_type), POINTER :: print_key CHARACTER(len=*), INTENT(IN), OPTIONAL :: middle_name CHARACTER(len=*), INTENT(IN) :: extension - LOGICAL, INTENT(IN) :: my_local + LOGICAL, INTENT(IN), OPTIONAL :: my_local, iter_string CHARACTER(len=default_path_length) :: filename CHARACTER(len=default_path_length) :: outName, outPath, postfix, root @@ -238,13 +252,21 @@ CONTAINS root = outPath(1:my_ind1)//my_middle_name(2:LEN_TRIM(my_middle_name)) END IF - postfix = extension + IF (PRESENT(iter_string) .AND. iter_string) THEN + ! use the cp_iter_string as a postfix + postfix = "-"//TRIM(cp_iter_string(logger%iter_info, print_key=print_key, for_file=.TRUE.)) + IF (TRIM(postfix) == "-") postfix = "" + ! and add the extension + postfix = TRIM(postfix)//extension + ELSE + postfix = extension + END IF ! and let the logger generate the filename CALL cp_logger_generate_filename(logger, res=filename, & root=root, postfix=postfix, local=my_local) - END FUNCTION negf_elecrodes_generate_filename + END FUNCTION negf_generate_filename ! ************************************************************************************************** !> \brief Prints full matrix to a file. diff --git a/src/negf_methods.F b/src/negf_methods.F index 482aeddbdf..8d8c55a554 100644 --- a/src/negf_methods.F +++ b/src/negf_methods.F @@ -37,16 +37,15 @@ MODULE negf_methods cp_fm_struct_type USE cp_fm_types, ONLY: & cp_fm_add_to_element, cp_fm_copy_general, cp_fm_create, cp_fm_get_info, & - cp_fm_get_submatrix, cp_fm_release, cp_fm_set_all, cp_fm_to_fm, cp_fm_type + cp_fm_get_submatrix, cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, & + cp_fm_type USE cp_log_handling, ONLY: cp_get_default_logger,& cp_logger_get_default_io_unit,& cp_logger_type - USE cp_output_handling, ONLY: cp_p_file,& - cp_print_key_finished_output,& - cp_print_key_should_output,& - cp_print_key_unit_nr,& - debug_print_level,& - high_print_level + USE cp_output_handling, ONLY: & + cp_add_iter_level, cp_iterate, cp_p_file, cp_print_key_finished_output, & + cp_print_key_should_output, cp_print_key_unit_nr, cp_rm_iter_level, debug_print_level, & + high_print_level USE cp_subsys_types, ONLY: cp_subsys_type USE force_env_types, ONLY: force_env_get,& force_env_p_type,& @@ -57,7 +56,8 @@ MODULE negf_methods USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type,& section_vals_val_get - USE kinds, ONLY: default_string_length,& + USE kinds, ONLY: default_path_length,& + default_string_length,& dp USE kpoint_types, ONLY: get_kpoint_info,& kpoint_type @@ -96,6 +96,8 @@ MODULE negf_methods simpsonrule_type,& sr_shape_arc,& sr_shape_linear + USE negf_io, ONLY: negf_read_matrix_from_file,& + negf_restart_file_name USE negf_matrix_utils, ONLY: invert_cell_to_index,& negf_copy_fm_submat_to_dbcsr,& negf_copy_sym_dbcsr_to_fm_submat @@ -291,7 +293,8 @@ CONTAINS ! self-consistent density ! ----------------------- - CALL converge_density(negf_env, negf_control, sub_env, qs_env, negf_control%v_shift, base_contact=1, log_unit=log_unit) + CALL converge_density(negf_env, negf_control, sub_env, negf_section, qs_env, negf_control%v_shift, & + base_contact=1, log_unit=log_unit) ! restart.hs ! ---------- @@ -1056,6 +1059,7 @@ CONTAINS !> \param negf_env NEGF environment !> \param negf_control NEGF control !> \param sub_env NEGF parallel (sub)group environment +!> \param negf_section ... !> \param qs_env QuickStep environment !> \param v_shift shift in Hartree potential !> \param base_contact index of the reference contact @@ -1064,10 +1068,11 @@ CONTAINS !> * 06.2017 created [Sergey Chulkov] !> * 11.2025 modified [Dmitry Ryndyk] ! ************************************************************************************************** - SUBROUTINE converge_density(negf_env, negf_control, sub_env, qs_env, v_shift, base_contact, log_unit) + SUBROUTINE converge_density(negf_env, negf_control, sub_env, negf_section, qs_env, v_shift, base_contact, log_unit) TYPE(negf_env_type), INTENT(inout) :: negf_env TYPE(negf_control_type), POINTER :: negf_control TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env + TYPE(section_vals_type), POINTER :: negf_section TYPE(qs_environment_type), POINTER :: qs_env REAL(kind=dp), INTENT(in) :: v_shift INTEGER, INTENT(in) :: base_contact, log_unit @@ -1076,18 +1081,22 @@ CONTAINS REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*EPSILON(0.0_dp) TYPE(cp_fm_type), PARAMETER :: fm_dummy = cp_fm_type() + CHARACTER(len=100) :: sfmt + CHARACTER(LEN=default_path_length) :: filebase, filename COMPLEX(kind=dp) :: lbound_cpath, ubound_cpath, ubound_lpath - INTEGER :: handle, icontact, image, ispin, & - iter_count, nao, ncontacts, nimages, & - nspins - LOGICAL :: do_kpoints + INTEGER :: handle, i, icontact, image, ispin, & + iter_count, j, nao, ncol, ncontacts, & + nimages, nrow, nspins, print_unit + LOGICAL :: do_kpoints, exist REAL(kind=dp) :: delta, iter_delta, mu_base, nelectrons, & nelectrons_diff, t1, t2, temperature, & trace, v_base + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: target_m TYPE(cp_blacs_env_type), POINTER :: blacs_env TYPE(cp_fm_struct_type), POINTER :: fm_struct TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: rho_ao_delta_fm, rho_ao_new_fm TYPE(cp_fm_type), POINTER :: matrix_s_fm + TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_initial_kp, matrix_ks_qs_kp, & rho_ao_initial_kp, rho_ao_new_kp, & rho_ao_qs_kp @@ -1100,6 +1109,8 @@ CONTAINS TYPE(qs_rho_type), POINTER :: rho_struct TYPE(qs_subsys_type), POINTER :: subsys + logger => cp_get_default_logger() + ncontacts = SIZE(negf_control%contacts) ! the current subroutine works for the general case as well, but the Poisson solver does not IF (ncontacts > 2) THEN @@ -1129,10 +1140,10 @@ CONTAINS t1 = m_walltime() ! need a globally distributed overlap matrix in order to compute integration errors + CALL cp_fm_get_info(negf_env%s_s, nrow_global=nao) IF (sub_env%ngroups > 1) THEN NULLIFY (matrix_s_fm, fm_struct) - CALL cp_fm_get_info(negf_env%s_s, nrow_global=nao) CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env) ALLOCATE (matrix_s_fm) @@ -1156,9 +1167,115 @@ CONTAINS v_base = negf_control%contacts(base_contact)%v_external mu_base = negf_control%contacts(base_contact)%fermi_level - v_base - ! keep the initial charge density matrix and Kohn-Sham matrix + ! keep the initial charge density matrix CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_qs_kp) + ALLOCATE (target_m(nao, nao)) + ALLOCATE (rho_ao_delta_fm(nspins), rho_ao_new_fm(nspins)) + DO ispin = 1, nspins + CALL cp_fm_create(rho_ao_delta_fm(ispin), fm_struct) + CALL cp_fm_create(rho_ao_new_fm(ispin), fm_struct) + END DO + + IF (negf_control%restart_scf) THEN + IF (para_env%is_source()) THEN + CALL negf_restart_file_name(filebase, exist, negf_section, logger, h_scf=.TRUE.) + END IF + CALL para_env%bcast(filebase) + IF (nspins == 1) THEN + filename = TRIM(filebase)//'.hs' + INQUIRE (FILE=filename, exist=exist) + IF (.NOT. exist) THEN + CALL cp_warn(__LOCATION__, & + "User requested to read the KS matrix from the file named: "// & + TRIM(filename)//". This file does not exist. The initial KS matrix will be used.") + ELSE + IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename, target_m) + CALL para_env%bcast(target_m) + CALL cp_fm_set_submatrix(negf_env%h_s(1), target_m) + IF (log_unit > 0) WRITE (log_unit, '(T2,A)') " H_s is read from "//TRIM(filename) + END IF + filename = TRIM(filebase)//'.rho' + INQUIRE (FILE=filename, exist=exist) + IF (.NOT. exist) THEN + CALL cp_warn(__LOCATION__, & + "User requested to read the density matrix from the file named: "// & + TRIM(filename)//". This file does not exist. The initial density matrix will be used.") + ELSE + IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename, target_m) + CALL para_env%bcast(target_m) + CALL cp_fm_set_submatrix(rho_ao_delta_fm(1), target_m) + IF (log_unit > 0) WRITE (log_unit, '(T2,A)') " rho_s is read from "//TRIM(filename) + CALL negf_copy_fm_submat_to_dbcsr(fm=rho_ao_delta_fm(1), & + matrix=rho_ao_qs_kp(1, 1)%matrix, & + atomlist_row=negf_control%atomlist_S_screening, & + atomlist_col=negf_control%atomlist_S_screening, & + subsys=subsys) + END IF + END IF + IF (nspins == 2) THEN + filename = TRIM(filebase)//'-S1.hs' + INQUIRE (FILE=filename, exist=exist) + IF (.NOT. exist) THEN + CALL cp_warn(__LOCATION__, & + "User requested to read the KS matrix from the file named: "// & + TRIM(filename)//". This file does not exist. The initial KS matrix will be used.") + ELSE + IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename, target_m) + CALL para_env%bcast(target_m) + CALL cp_fm_set_submatrix(negf_env%h_s(1), target_m) + IF (log_unit > 0) WRITE (log_unit, '(T2,A)') " H_s is read from "//TRIM(filename) + END IF + filename = TRIM(filebase)//'-S2.hs' + INQUIRE (FILE=filename, exist=exist) + IF (.NOT. exist) THEN + CALL cp_warn(__LOCATION__, & + "User requested to read the KS matrix from the file named: "// & + TRIM(filename)//". This file does not exist. The initial KS matrix will be used.") + ELSE + IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename, target_m) + CALL para_env%bcast(target_m) + CALL cp_fm_set_submatrix(negf_env%h_s(2), target_m) + IF (log_unit > 0) WRITE (log_unit, '(T2,A)') " H_s is read from "//TRIM(filename) + END IF + filename = TRIM(filebase)//'-S1.rho' + INQUIRE (FILE=filename, exist=exist) + IF (.NOT. exist) THEN + CALL cp_warn(__LOCATION__, & + "User requested to read the density matrix from the file named: "// & + TRIM(filename)//". This file does not exist. The initial density matrix will be used.") + ELSE + IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename, target_m) + CALL para_env%bcast(target_m) + CALL cp_fm_set_submatrix(rho_ao_delta_fm(1), target_m) + IF (log_unit > 0) WRITE (log_unit, '(T2,A)') " rho_s is read from "//TRIM(filename) + CALL negf_copy_fm_submat_to_dbcsr(fm=rho_ao_delta_fm(1), & + matrix=rho_ao_qs_kp(1, 1)%matrix, & + atomlist_row=negf_control%atomlist_S_screening, & + atomlist_col=negf_control%atomlist_S_screening, & + subsys=subsys) + END IF + filename = TRIM(filebase)//'-S2.rho' + INQUIRE (FILE=filename, exist=exist) + IF (.NOT. exist) THEN + CALL cp_warn(__LOCATION__, & + "User requested to read the density matrix from the file named: "// & + TRIM(filename)//". This file does not exist. The initial density matrix will be used.") + ELSE + IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename, target_m) + CALL para_env%bcast(target_m) + CALL cp_fm_set_submatrix(rho_ao_delta_fm(2), target_m) + IF (log_unit > 0) WRITE (log_unit, '(T2,A)') " rho_s is read from "//TRIM(filename) + CALL negf_copy_fm_submat_to_dbcsr(fm=rho_ao_delta_fm(2), & + matrix=rho_ao_qs_kp(2, 1)%matrix, & + atomlist_row=negf_control%atomlist_S_screening, & + atomlist_col=negf_control%atomlist_S_screening, & + subsys=subsys) + END IF + END IF + CALL qs_rho_update_rho(rho_struct, qs_env=qs_env) + END IF + NULLIFY (matrix_ks_initial_kp, rho_ao_initial_kp, rho_ao_new_kp) CALL dbcsr_allocate_matrix_set(matrix_ks_initial_kp, nspins, nimages) CALL dbcsr_allocate_matrix_set(rho_ao_initial_kp, nspins, nimages) @@ -1179,27 +1296,18 @@ CONTAINS ! extract the reference density matrix blocks nelectrons = 0.0_dp - ALLOCATE (rho_ao_delta_fm(nspins), rho_ao_new_fm(nspins)) DO ispin = 1, nspins - CALL cp_fm_create(rho_ao_delta_fm(ispin), fm_struct) - CALL cp_fm_create(rho_ao_new_fm(ispin), fm_struct) - END DO - IF (.NOT. negf_control%is_restart) THEN - DO ispin = 1, nspins - CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_qs_kp(ispin, 1)%matrix, & - fm=rho_ao_delta_fm(ispin), & - atomlist_row=negf_control%atomlist_S_screening, & - atomlist_col=negf_control%atomlist_S_screening, & - subsys=subsys, mpi_comm_global=para_env, & - do_upper_diag=.TRUE., do_lower=.TRUE.) + CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_qs_kp(ispin, 1)%matrix, & + fm=rho_ao_delta_fm(ispin), & + atomlist_row=negf_control%atomlist_S_screening, & + atomlist_col=negf_control%atomlist_S_screening, & + subsys=subsys, mpi_comm_global=para_env, & + do_upper_diag=.TRUE., do_lower=.TRUE.) - CALL cp_fm_trace(rho_ao_delta_fm(ispin), matrix_s_fm, trace) - nelectrons = nelectrons + trace - END DO - negf_env%nelectrons = nelectrons - ELSE - nelectrons = negf_env%nelectrons - END IF + CALL cp_fm_trace(rho_ao_delta_fm(ispin), matrix_s_fm, trace) + nelectrons = nelectrons + trace + END DO + negf_env%nelectrons = nelectrons ! mixing storage allocation IF (negf_env%mixing_method >= gspace_mixing_nr) THEN @@ -1231,11 +1339,13 @@ CONTAINS REAL(negf_control%delta_npoles, kind=dp)*twopi*temperature, kind=dp) ALLOCATE (g_surf_circular(nspins), g_surf_linear(nspins), g_surf_nonequiv(nspins)) + CALL cp_add_iter_level(logger%iter_info, "NEGF_SCF") !--- main SCF cycle -------------------------------------------------------------------! DO iter_count = 1, negf_control%max_scf ! compute an updated density matrix CALL integration_status_reset(stats) + CALL cp_iterate(logger%iter_info, last=.FALSE., iter_nr=iter_count) DO ispin = 1, nspins ! closed contour: residuals @@ -1392,8 +1502,8 @@ CONTAINS rho_struct, para_env, iter_count) END IF + ! update KS-matrix IF (negf_control%update_HS) THEN - ! update KS-matrix CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.) ! extract blocks from the updated Kohn-Sham matrix DO ispin = 1, nspins @@ -1406,8 +1516,216 @@ CONTAINS END DO END IF + ! Write the HS restart files + IF (nspins == 1) THEN + CALL cp_fm_get_submatrix(negf_env%h_s(1), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension=".hs", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + END IF + IF (nspins == 2) THEN + CALL cp_fm_get_submatrix(negf_env%h_s(1), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension="-S1.hs", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + CALL cp_fm_get_submatrix(negf_env%h_s(2), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension="-S2.hs", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + END IF + + ! Write the rho restart files + IF (nspins == 1) THEN + CALL cp_fm_get_submatrix(rho_ao_new_fm(1), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension=".rho", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + END IF + IF (nspins == 2) THEN + CALL cp_fm_get_submatrix(rho_ao_new_fm(1), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension="-S1.rho", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + CALL cp_fm_get_submatrix(rho_ao_new_fm(2), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension="-S2.rho", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + END IF + END DO + ! Write the final HS restart files + CALL cp_iterate(logger%iter_info, last=.TRUE., iter_nr=iter_count) + IF (nspins == 1) THEN + CALL cp_fm_get_submatrix(negf_env%h_s(1), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension=".hs", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + END IF + IF (nspins == 2) THEN + CALL cp_fm_get_submatrix(negf_env%h_s(1), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension="-S1.hs", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + CALL cp_fm_get_submatrix(negf_env%h_s(1), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension="-S2.hs", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + END IF + + ! Write the final rho restart files + IF (nspins == 1) THEN + CALL cp_fm_get_submatrix(rho_ao_new_fm(1), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension=".rho", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + END IF + IF (nspins == 2) THEN + CALL cp_fm_get_submatrix(rho_ao_new_fm(1), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension="-S1.rho", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + CALL cp_fm_get_submatrix(rho_ao_new_fm(2), target_m) + IF (para_env%is_source() .AND. BTEST(cp_print_key_should_output(logger%iter_info, & + negf_section, 'PRINT%RESTART'), cp_p_file)) THEN + print_unit = cp_print_key_unit_nr(logger, negf_section, 'PRINT%RESTART', & + extension="-S2.rho", file_status="REPLACE", file_action="WRITE", & + do_backup=.TRUE., file_form="FORMATTED") + nrow = SIZE(target_m, 1) + ncol = SIZE(target_m, 2) + WRITE (sfmt, "('(',i0,'(E15.5))')") ncol + WRITE (print_unit, *) nrow, ncol + DO i = 1, nrow + WRITE (print_unit, sfmt) (target_m(i, j), j=1, ncol) + END DO + CALL cp_print_key_finished_output(print_unit, logger, negf_section, 'PRINT%RESTART') + END IF + END IF + + DEALLOCATE (target_m) + CALL cp_rm_iter_level(logger%iter_info, level_name="NEGF_SCF") + !--------------------------------------------------------------------------------------! IF (log_unit > 0) THEN