Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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_dot,&
18 : dbcsr_p_type
19 : USE cp_log_handling, ONLY: cp_get_default_logger,&
20 : cp_logger_get_default_io_unit,&
21 : cp_logger_type
22 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
23 : section_vals_type
24 : USE kinds, ONLY: dp
25 : USE kpoint_types, ONLY: kpoint_type
26 : USE machine, ONLY: m_walltime
27 : USE message_passing, ONLY: mp_para_env_type
28 : USE qs_core_energies, ONLY: calculate_ptrace
29 : USE qs_energy_types, ONLY: qs_energy_type
30 : USE qs_environment_types, ONLY: get_qs_env,&
31 : qs_environment_type,&
32 : set_qs_env
33 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
34 : USE qs_ks_types, ONLY: qs_ks_did_change,&
35 : qs_ks_env_type
36 : USE qs_mo_types, ONLY: mo_set_type
37 : USE qs_rho_types, ONLY: qs_rho_get,&
38 : qs_rho_type
39 : USE qs_scf, ONLY: init_scf_loop
40 : USE qs_scf_initialization, ONLY: qs_scf_env_initialize
41 : USE qs_scf_loop_utils, ONLY: qs_scf_new_mos,&
42 : qs_scf_new_mos_kp
43 : USE qs_scf_output, ONLY: qs_scf_loop_print,&
44 : qs_scf_print_summary,&
45 : qs_scf_write_mos
46 : USE qs_scf_post_scf, ONLY: qs_scf_compute_properties
47 : USE qs_scf_types, ONLY: qs_scf_env_type
48 : USE qs_wf_history_methods, ONLY: wfi_update
49 : USE scf_control_types, ONLY: scf_control_type
50 : #include "./base/base_uses.f90"
51 :
52 : IMPLICIT NONE
53 :
54 : PRIVATE
55 :
56 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_nonscf'
57 :
58 : PUBLIC :: nonscf
59 :
60 : CONTAINS
61 :
62 : ! **************************************************************************************************
63 : !> \brief Find solution to HC=SCE
64 : !> \param qs_env the qs_environment where to perform the scf procedure
65 : !> \par History
66 : !> none
67 : !> \author JGH
68 : !> \note
69 : ! **************************************************************************************************
70 576 : SUBROUTINE nonscf(qs_env)
71 : TYPE(qs_environment_type), POINTER :: qs_env
72 :
73 : TYPE(dft_control_type), POINTER :: dft_control
74 : TYPE(qs_scf_env_type), POINTER :: scf_env
75 : TYPE(scf_control_type), POINTER :: scf_control
76 :
77 : CALL get_qs_env(qs_env, scf_env=scf_env, scf_control=scf_control, &
78 576 : dft_control=dft_control)
79 :
80 576 : IF (dft_control%qs_control%do_ls_scf) THEN
81 : ! Density matrix based solver
82 :
83 0 : CPABORT("NOT AVAILABLE")
84 : !CALL ls_scf(qs_env=qs_env)
85 :
86 : ELSE
87 : ! Wavefunction based solver
88 :
89 576 : IF (.NOT. ASSOCIATED(scf_env)) THEN
90 58 : CALL qs_scf_env_initialize(qs_env, scf_env)
91 58 : CALL set_qs_env(qs_env, scf_env=scf_env)
92 : ELSE
93 518 : CALL qs_scf_env_initialize(qs_env, scf_env)
94 : END IF
95 :
96 576 : CALL do_nonscf(qs_env, scf_env, scf_control)
97 :
98 : ! add the converged wavefunction to the wavefunction history
99 576 : IF (ASSOCIATED(qs_env%wf_history)) THEN
100 576 : CALL wfi_update(qs_env%wf_history, qs_env=qs_env, dt=1.0_dp)
101 : END IF
102 :
103 : ! compute properties that depend on the wavefunction
104 576 : CALL qs_scf_compute_properties(qs_env)
105 :
106 : END IF
107 :
108 576 : END SUBROUTINE nonscf
109 :
110 : ! **************************************************************************************************
111 : !> \brief Solve KS equation for fixed potential
112 : !> \param qs_env ...
113 : !> \param scf_env the scf_env where to perform the scf procedure
114 : !> \param scf_control ...
115 : !> \par History
116 : !> none
117 : !> \author JGH
118 : !> \note
119 : ! **************************************************************************************************
120 576 : SUBROUTINE do_nonscf(qs_env, scf_env, scf_control)
121 :
122 : TYPE(qs_environment_type), POINTER :: qs_env
123 : TYPE(qs_scf_env_type), POINTER :: scf_env
124 : TYPE(scf_control_type), POINTER :: scf_control
125 :
126 : CHARACTER(LEN=*), PARAMETER :: routineN = 'do_nonscf'
127 :
128 : INTEGER :: handle, img, ispin, output_unit
129 : LOGICAL :: diis_step, do_kpoints
130 : REAL(KIND=dp) :: pc_ener, qmmm_el, t1, t2
131 : TYPE(cp_logger_type), POINTER :: logger
132 576 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_ks, rho_ao_kp
133 : TYPE(dft_control_type), POINTER :: dft_control
134 : TYPE(kpoint_type), POINTER :: kpoints
135 576 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
136 : TYPE(mp_para_env_type), POINTER :: para_env
137 : TYPE(qs_energy_type), POINTER :: energy
138 : TYPE(qs_ks_env_type), POINTER :: ks_env
139 : TYPE(qs_rho_type), POINTER :: rho
140 : TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
141 :
142 576 : CALL timeset(routineN, handle)
143 :
144 576 : t1 = m_walltime()
145 :
146 576 : logger => cp_get_default_logger()
147 576 : output_unit = cp_logger_get_default_io_unit(logger)
148 :
149 : CALL get_qs_env(qs_env=qs_env, &
150 : energy=energy, &
151 : ks_env=ks_env, &
152 : rho=rho, &
153 : mos=mos, &
154 : input=input, &
155 : dft_control=dft_control, &
156 : do_kpoints=do_kpoints, &
157 : kpoints=kpoints, &
158 576 : para_env=para_env)
159 :
160 1152 : DO ispin = 1, dft_control%nspins
161 1152 : CPASSERT(.NOT. mos(ispin)%use_mo_coeff_b)
162 : END DO
163 :
164 576 : dft_section => section_vals_get_subs_vals(input, "DFT")
165 576 : scf_section => section_vals_get_subs_vals(dft_section, "SCF")
166 576 : CALL init_scf_loop(scf_env=scf_env, qs_env=qs_env, scf_section=scf_section)
167 :
168 : ! Calculate KS matrix
169 576 : CALL qs_ks_update_qs_env(qs_env, just_energy=.FALSE., calculate_forces=.FALSE.)
170 :
171 : ! print 'heavy weight' or relatively expensive quantities
172 576 : CALL qs_scf_loop_print(qs_env, scf_env, para_env)
173 :
174 : ! Diagonalization
175 576 : IF (do_kpoints) THEN
176 : ! kpoints
177 100 : CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step)
178 : ELSE
179 : ! Gamma points only
180 476 : CALL qs_scf_new_mos(qs_env, scf_env, scf_control, scf_section, diis_step, .FALSE.)
181 : END IF
182 :
183 : ! Print requested MO information (can be computationally expensive with OT)
184 576 : CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.TRUE.)
185 :
186 : ! copy density matrix
187 576 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
188 1152 : DO ispin = 1, dft_control%nspins
189 3722 : DO img = 1, SIZE(rho_ao_kp, 2)
190 3146 : CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, scf_env%p_mix_new(ispin, img)%matrix)
191 : END DO
192 : END DO
193 576 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE., potential_changed=.TRUE.)
194 :
195 : ! core energy : Tr(PH)
196 576 : CALL get_qs_env(qs_env, matrix_ks_kp=matrixkp_ks)
197 576 : IF (qs_env%qmmm) THEN
198 : ! Compute QM/MM Energy
199 336 : CPASSERT(SIZE(matrixkp_ks, 2) == 1)
200 672 : DO ispin = 1, dft_control%nspins
201 : CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
202 336 : matrixkp_ks(ispin, 1)%matrix, qmmm_el)
203 672 : energy%qmmm_el = energy%qmmm_el + qmmm_el
204 : END DO
205 336 : pc_ener = qs_env%ks_qmmm_env%pc_ener
206 336 : energy%qmmm_el = energy%qmmm_el + pc_ener
207 : ELSE
208 240 : energy%qmmm_el = 0.0_dp
209 : END IF
210 576 : energy%total = energy%total - energy%core
211 576 : CALL calculate_ptrace(matrixkp_ks, rho_ao_kp, energy%core, dft_control%nspins)
212 576 : energy%total = energy%total + energy%core + energy%qmmm_el
213 :
214 576 : t2 = m_walltime()
215 :
216 576 : IF (output_unit > 0) THEN
217 : WRITE (UNIT=output_unit, &
218 : FMT="(T2,A,T40,A,F10.2,T61,F20.10)") &
219 288 : "Diagonalization", "Time:", t2 - t1, energy%total
220 : END IF
221 :
222 576 : CALL qs_scf_print_summary(output_unit, qs_env)
223 :
224 576 : CALL timestop(handle)
225 :
226 576 : END SUBROUTINE do_nonscf
227 :
228 : END MODULE qs_nonscf
|
