comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

NEGF: restart for NEGF SCF (#5023)

Browse source
This commit is contained in:
Dmitry Ryndyk 2026-03-30 10:38:21 +02:00 • committed by GitHub
parent ed08826bad
commit d10c76ea8a
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
6 changed files with 451 additions and 78 deletions
Download patch file Download diff file Expand all files Collapse all files

10
src/common/cp_iter_types.F
View file

@ -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

34
src/input_cp2k_negf.F
View file

@ -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
! **************************************************************************************************

7
src/negf_control_types.F
View file

@ -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

16
src/negf_env_types.F
View file

@ -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

72
src/negf_io.F
View file

@ -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.

390
src/negf_methods.F
View file

@ -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