Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for Quickstep NON-SCF run.
10 : !> \par History
11 : !> - initial setup [JGH, 2024]
12 : !> \author JGH (13.05.2024)
13 : ! **************************************************************************************************
14 : MODULE qs_nonscf
15 : USE cp_control_types, ONLY: dft_control_type
16 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
17 : dbcsr_p_type
18 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
19 : USE cp_log_handling, ONLY: cp_get_default_logger,&
20 : cp_logger_get_default_io_unit,&
21 : cp_logger_type
22 : USE dm_ls_scf, ONLY: ls_scf
23 : USE input_constants, ONLY: hfit_relative_entropy
24 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
25 : section_vals_type
26 : USE kinds, ONLY: dp
27 : USE kpoint_types, ONLY: kpoint_type
28 : USE machine, ONLY: m_walltime
29 : USE message_passing, ONLY: mp_para_env_type
30 : USE qs_core_energies, ONLY: calculate_ptrace
31 : USE qs_energy_types, ONLY: qs_energy_type
32 : USE qs_environment_types, ONLY: get_qs_env,&
33 : qs_environment_type,&
34 : set_qs_env
35 : USE qs_harris_methods, ONLY: harris_direct_density_matrix_energy,&
36 : harris_relative_entropy_reconstruction
37 : USE qs_harris_types, ONLY: harris_type
38 : USE qs_harris_utils, ONLY: harris_density_update
39 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
40 : USE qs_ks_types, ONLY: qs_ks_did_change,&
41 : qs_ks_env_type
42 : USE qs_mo_types, ONLY: mo_set_type
43 : USE qs_nonscf_utils, ONLY: qs_nonscf_print_summary
44 : USE qs_rho_methods, ONLY: qs_rho_update_rho
45 : USE qs_rho_types, ONLY: qs_rho_get,&
46 : qs_rho_type
47 : USE qs_scf, ONLY: init_scf_loop
48 : USE qs_scf_initialization, ONLY: qs_scf_env_initialize
49 : USE qs_scf_loop_utils, ONLY: qs_scf_new_mos,&
50 : qs_scf_new_mos_kp
51 : USE qs_scf_output, ONLY: qs_scf_loop_print,&
52 : qs_scf_write_mos
53 : USE qs_scf_post_scf, ONLY: qs_scf_compute_properties
54 : USE qs_scf_types, ONLY: qs_scf_env_type
55 : USE qs_wf_history_methods, ONLY: wfi_update
56 : USE scf_control_types, ONLY: scf_control_type
57 : #include "./base/base_uses.f90"
58 :
59 : IMPLICIT NONE
60 :
61 : PRIVATE
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_nonscf'
64 :
65 : PUBLIC :: nonscf
66 :
67 : CONTAINS
68 :
69 : ! **************************************************************************************************
70 : !> \brief Find solution to HC=SCE
71 : !> \param qs_env the qs_environment where to perform the scf procedure
72 : !> \par History
73 : !> none
74 : !> \author JGH
75 : !> \note
76 : ! **************************************************************************************************
77 2916 : SUBROUTINE nonscf(qs_env)
78 : TYPE(qs_environment_type), POINTER :: qs_env
79 :
80 : TYPE(dft_control_type), POINTER :: dft_control
81 : TYPE(harris_type), POINTER :: harris_env
82 : TYPE(qs_rho_type), POINTER :: rho
83 : TYPE(qs_scf_env_type), POINTER :: scf_env
84 : TYPE(scf_control_type), POINTER :: scf_control
85 :
86 2916 : NULLIFY (harris_env, rho)
87 2916 : CALL get_qs_env(qs_env, dft_control=dft_control, harris_env=harris_env, rho=rho)
88 :
89 2916 : IF (ASSOCIATED(harris_env)) THEN
90 2916 : IF (harris_env%fit_method == hfit_relative_entropy .AND. &
91 : .NOT. harris_env%density_fit_ready) THEN
92 8 : CALL harris_density_update(qs_env, harris_env)
93 8 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
94 8 : CALL harris_relative_entropy_reconstruction(qs_env)
95 : END IF
96 2916 : IF (harris_env%direct_density_matrix_energy) THEN
97 8 : CALL harris_density_update(qs_env, harris_env)
98 8 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
99 8 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE., potential_changed=.TRUE.)
100 8 : CALL harris_direct_density_matrix_energy(qs_env)
101 8 : RETURN
102 : END IF
103 : END IF
104 :
105 2908 : IF (dft_control%qs_control%do_ls_scf) THEN
106 : ! Density matrix based solver
107 :
108 72 : CALL ls_scf(qs_env, nonscf=.TRUE.)
109 :
110 : ELSE
111 : ! Wavefunction based solver
112 :
113 2836 : CALL get_qs_env(qs_env, scf_env=scf_env, scf_control=scf_control)
114 2836 : IF (.NOT. ASSOCIATED(scf_env)) THEN
115 772 : CALL qs_scf_env_initialize(qs_env, scf_env)
116 772 : CALL set_qs_env(qs_env, scf_env=scf_env)
117 : ELSE
118 2064 : CALL qs_scf_env_initialize(qs_env, scf_env)
119 : END IF
120 :
121 2836 : CALL do_nonscf(qs_env, scf_env, scf_control)
122 :
123 : ! add the converged wavefunction to the wavefunction history
124 2836 : IF (ASSOCIATED(qs_env%wf_history)) THEN
125 2836 : CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
126 : END IF
127 :
128 : ! compute properties that depend on the wavefunction
129 2836 : CALL qs_scf_compute_properties(qs_env)
130 :
131 : END IF
132 :
133 : END SUBROUTINE nonscf
134 :
135 : ! **************************************************************************************************
136 : !> \brief Solve KS equation for fixed potential
137 : !> \param qs_env ...
138 : !> \param scf_env the scf_env where to perform the scf procedure
139 : !> \param scf_control ...
140 : !> \par History
141 : !> none
142 : !> \author JGH
143 : !> \note
144 : ! **************************************************************************************************
145 2836 : SUBROUTINE do_nonscf(qs_env, scf_env, scf_control)
146 :
147 : TYPE(qs_environment_type), POINTER :: qs_env
148 : TYPE(qs_scf_env_type), POINTER :: scf_env
149 : TYPE(scf_control_type), POINTER :: scf_control
150 :
151 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_nonscf'
152 :
153 : INTEGER :: handle, img, iounit, ispin
154 : LOGICAL :: diis_step, do_kpoints
155 : REAL(KIND=dp) :: pc_ener, qmmm_el, t1, t2, tdiag
156 : TYPE(cp_logger_type), POINTER :: logger
157 2836 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrixkp_ks, rho_ao_kp
158 : TYPE(dft_control_type), POINTER :: dft_control
159 : TYPE(kpoint_type), POINTER :: kpoints
160 2836 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
161 : TYPE(mp_para_env_type), POINTER :: para_env
162 : TYPE(qs_energy_type), POINTER :: energy
163 : TYPE(qs_ks_env_type), POINTER :: ks_env
164 : TYPE(qs_rho_type), POINTER :: rho
165 : TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
166 :
167 2836 : CALL timeset(routineN, handle)
168 :
169 2836 : t1 = m_walltime()
170 :
171 2836 : logger => cp_get_default_logger()
172 2836 : iounit = cp_logger_get_default_io_unit(logger)
173 :
174 : CALL get_qs_env(qs_env=qs_env, &
175 : energy=energy, &
176 : ks_env=ks_env, &
177 : rho=rho, &
178 : mos=mos, &
179 : input=input, &
180 : dft_control=dft_control, &
181 : do_kpoints=do_kpoints, &
182 : kpoints=kpoints, &
183 2836 : para_env=para_env)
184 :
185 6044 : DO ispin = 1, dft_control%nspins
186 6044 : CPASSERT(.NOT. mos(ispin)%use_mo_coeff_b)
187 : END DO
188 :
189 2836 : dft_section => section_vals_get_subs_vals(input, "DFT")
190 2836 : scf_section => section_vals_get_subs_vals(dft_section, "SCF")
191 2836 : CALL init_scf_loop(scf_env=scf_env, qs_env=qs_env, scf_section=scf_section)
192 :
193 : ! Calculate KS matrix
194 2836 : CALL qs_ks_update_qs_env(qs_env, just_energy=.FALSE., calculate_forces=.FALSE.)
195 :
196 : ! print 'heavy weight' or relatively expensive quantities
197 2836 : CALL qs_scf_loop_print(qs_env, scf_env, para_env)
198 :
199 : ! Diagonalization
200 2836 : IF (do_kpoints) THEN
201 : ! kpoints
202 566 : CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step)
203 : ELSE
204 : ! Gamma points only
205 2270 : CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, .FALSE.)
206 : END IF
207 :
208 : ! Print requested MO information (can be computationally expensive with OT)
209 2836 : CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.TRUE.)
210 :
211 : ! copy density matrix
212 2836 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
213 6044 : DO ispin = 1, dft_control%nspins
214 44606 : DO img = 1, SIZE(rho_ao_kp, 2)
215 41770 : CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, scf_env%p_mix_new(ispin, img)%matrix)
216 : END DO
217 : END DO
218 :
219 2836 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE., potential_changed=.TRUE.)
220 :
221 : ! band energy : Tr(PH)
222 2836 : CALL get_qs_env(qs_env, matrix_ks_kp=matrixkp_ks)
223 2836 : CALL calculate_ptrace(matrixkp_ks, rho_ao_kp, energy%band, dft_control%nspins, .TRUE.)
224 : ! core energy : Tr(Ph)
225 2836 : energy%total = energy%total - energy%core
226 2836 : CALL get_qs_env(qs_env, matrix_h_kp=matrix_h)
227 2836 : CALL calculate_ptrace(matrix_h, rho_ao_kp, energy%core, dft_control%nspins)
228 :
229 2836 : IF (qs_env%qmmm) THEN
230 : ! Compute QM/MM Energy
231 336 : CPASSERT(SIZE(matrixkp_ks, 2) == 1)
232 672 : DO ispin = 1, dft_control%nspins
233 : CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
234 336 : matrixkp_ks(ispin, 1)%matrix, qmmm_el)
235 672 : energy%qmmm_el = energy%qmmm_el + qmmm_el
236 : END DO
237 336 : pc_ener = qs_env%ks_qmmm_env%pc_ener
238 336 : energy%qmmm_el = energy%qmmm_el + pc_ener
239 : ELSE
240 2500 : energy%qmmm_el = 0.0_dp
241 : END IF
242 :
243 2836 : t2 = m_walltime()
244 2836 : tdiag = t2 - t1
245 :
246 2836 : CALL qs_nonscf_print_summary(qs_env, tdiag, scf_env%nelectron, iounit)
247 :
248 2836 : CALL timestop(handle)
249 :
250 2836 : END SUBROUTINE do_nonscf
251 :
252 : END MODULE qs_nonscf
|
