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 dm_ls_scf, ONLY: ls_scf
23 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
24 : section_vals_type
25 : USE kinds, ONLY: dp
26 : USE kpoint_types, ONLY: kpoint_type
27 : USE machine, ONLY: m_walltime
28 : USE message_passing, ONLY: mp_para_env_type
29 : USE qs_core_energies, ONLY: calculate_ptrace
30 : USE qs_energy_types, ONLY: qs_energy_type
31 : USE qs_environment_types, ONLY: get_qs_env,&
32 : qs_environment_type,&
33 : set_qs_env
34 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
35 : USE qs_ks_types, ONLY: qs_ks_did_change,&
36 : qs_ks_env_type
37 : USE qs_mo_types, ONLY: mo_set_type
38 : USE qs_nonscf_utils, ONLY: qs_nonscf_print_summary
39 : USE qs_rho_types, ONLY: qs_rho_get,&
40 : qs_rho_type
41 : USE qs_scf, ONLY: init_scf_loop
42 : USE qs_scf_initialization, ONLY: qs_scf_env_initialize
43 : USE qs_scf_loop_utils, ONLY: qs_scf_new_mos,&
44 : qs_scf_new_mos_kp
45 : USE qs_scf_output, ONLY: qs_scf_loop_print,&
46 : qs_scf_write_mos
47 : USE qs_scf_post_scf, ONLY: qs_scf_compute_properties
48 : USE qs_scf_types, ONLY: qs_scf_env_type
49 : USE qs_wf_history_methods, ONLY: wfi_update
50 : USE scf_control_types, ONLY: scf_control_type
51 : #include "./base/base_uses.f90"
52 :
53 : IMPLICIT NONE
54 :
55 : PRIVATE
56 :
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_nonscf'
58 :
59 : PUBLIC :: nonscf
60 :
61 : CONTAINS
62 :
63 : ! **************************************************************************************************
64 : !> \brief Find solution to HC=SCE
65 : !> \param qs_env the qs_environment where to perform the scf procedure
66 : !> \par History
67 : !> none
68 : !> \author JGH
69 : !> \note
70 : ! **************************************************************************************************
71 2010 : SUBROUTINE nonscf(qs_env)
72 : TYPE(qs_environment_type), POINTER :: qs_env
73 :
74 : TYPE(dft_control_type), POINTER :: dft_control
75 : TYPE(qs_scf_env_type), POINTER :: scf_env
76 : TYPE(scf_control_type), POINTER :: scf_control
77 :
78 2010 : CALL get_qs_env(qs_env, dft_control=dft_control)
79 :
80 2010 : IF (dft_control%qs_control%do_ls_scf) THEN
81 : ! Density matrix based solver
82 :
83 66 : CALL ls_scf(qs_env, nonscf=.TRUE.)
84 :
85 : ELSE
86 : ! Wavefunction based solver
87 :
88 1944 : CALL get_qs_env(qs_env, scf_env=scf_env, scf_control=scf_control)
89 1944 : IF (.NOT. ASSOCIATED(scf_env)) THEN
90 714 : CALL qs_scf_env_initialize(qs_env, scf_env)
91 714 : CALL set_qs_env(qs_env, scf_env=scf_env)
92 : ELSE
93 1230 : CALL qs_scf_env_initialize(qs_env, scf_env)
94 : END IF
95 :
96 1944 : CALL do_nonscf(qs_env, scf_env, scf_control)
97 :
98 : ! add the converged wavefunction to the wavefunction history
99 1944 : IF (ASSOCIATED(qs_env%wf_history)) THEN
100 1944 : 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 1944 : CALL qs_scf_compute_properties(qs_env)
105 :
106 : END IF
107 :
108 2010 : 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 1944 : 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, iounit, ispin
129 : LOGICAL :: diis_step, do_kpoints
130 : REAL(KIND=dp) :: pc_ener, qmmm_el, t1, t2, tdiag
131 : TYPE(cp_logger_type), POINTER :: logger
132 1944 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrixkp_ks, rho_ao_kp
133 : TYPE(dft_control_type), POINTER :: dft_control
134 : TYPE(kpoint_type), POINTER :: kpoints
135 1944 : 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 1944 : CALL timeset(routineN, handle)
143 :
144 1944 : t1 = m_walltime()
145 :
146 1944 : logger => cp_get_default_logger()
147 1944 : iounit = 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 1944 : para_env=para_env)
159 :
160 4116 : DO ispin = 1, dft_control%nspins
161 4116 : CPASSERT(.NOT. mos(ispin)%use_mo_coeff_b)
162 : END DO
163 :
164 1944 : dft_section => section_vals_get_subs_vals(input, "DFT")
165 1944 : scf_section => section_vals_get_subs_vals(dft_section, "SCF")
166 1944 : CALL init_scf_loop(scf_env=scf_env, qs_env=qs_env, scf_section=scf_section)
167 :
168 : ! Calculate KS matrix
169 1944 : CALL qs_ks_update_qs_env(qs_env, just_energy=.FALSE., calculate_forces=.FALSE.)
170 :
171 : ! print 'heavy weight' or relatively expensive quantities
172 1944 : CALL qs_scf_loop_print(qs_env, scf_env, para_env)
173 :
174 : ! Diagonalization
175 1944 : IF (do_kpoints) THEN
176 : ! kpoints
177 238 : CALL qs_scf_new_mos_kp(qs_env, scf_env, scf_control, diis_step)
178 : ELSE
179 : ! Gamma points only
180 1706 : 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 1944 : CALL qs_scf_write_mos(qs_env, scf_env, final_mos=.TRUE.)
185 :
186 : ! copy density matrix
187 1944 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
188 4116 : DO ispin = 1, dft_control%nspins
189 34746 : DO img = 1, SIZE(rho_ao_kp, 2)
190 32802 : CALL dbcsr_copy(rho_ao_kp(ispin, img)%matrix, scf_env%p_mix_new(ispin, img)%matrix)
191 : END DO
192 : END DO
193 :
194 1944 : CALL qs_ks_did_change(ks_env, rho_changed=.TRUE., potential_changed=.TRUE.)
195 :
196 : ! band energy : Tr(PH)
197 1944 : CALL get_qs_env(qs_env, matrix_ks_kp=matrixkp_ks)
198 1944 : CALL calculate_ptrace(matrixkp_ks, rho_ao_kp, energy%band, dft_control%nspins)
199 : ! core energy : Tr(Ph)
200 1944 : energy%total = energy%total - energy%core
201 1944 : CALL get_qs_env(qs_env, matrix_h_kp=matrix_h)
202 1944 : CALL calculate_ptrace(matrix_h, rho_ao_kp, energy%core, dft_control%nspins)
203 :
204 1944 : IF (qs_env%qmmm) THEN
205 : ! Compute QM/MM Energy
206 336 : CPASSERT(SIZE(matrixkp_ks, 2) == 1)
207 672 : DO ispin = 1, dft_control%nspins
208 : CALL dbcsr_dot(qs_env%ks_qmmm_env%matrix_h(1)%matrix, &
209 336 : matrixkp_ks(ispin, 1)%matrix, qmmm_el)
210 672 : energy%qmmm_el = energy%qmmm_el + qmmm_el
211 : END DO
212 336 : pc_ener = qs_env%ks_qmmm_env%pc_ener
213 336 : energy%qmmm_el = energy%qmmm_el + pc_ener
214 : ELSE
215 1608 : energy%qmmm_el = 0.0_dp
216 : END IF
217 :
218 1944 : t2 = m_walltime()
219 1944 : tdiag = t2 - t1
220 :
221 1944 : CALL qs_nonscf_print_summary(qs_env, tdiag, scf_env%nelectron, iounit)
222 :
223 1944 : CALL timestop(handle)
224 :
225 1944 : END SUBROUTINE do_nonscf
226 :
227 : END MODULE qs_nonscf
|
