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 ...
10 : !> \author ...
11 : ! *****************************************************************************
12 : MODULE qs_basis_gradient
13 :
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind,&
16 : get_atomic_kind_set
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
19 : dbcsr_p_type,&
20 : dbcsr_set,&
21 : dbcsr_type
22 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
23 : USE cp_fm_types, ONLY: cp_fm_type
24 : USE cp_log_handling, ONLY: cp_get_default_logger,&
25 : cp_logger_type
26 : USE input_section_types, ONLY: section_vals_type
27 : USE kinds, ONLY: dp
28 : USE message_passing, ONLY: mp_para_env_type
29 : USE particle_types, ONLY: particle_type
30 : USE qs_core_hamiltonian, ONLY: build_core_hamiltonian_matrix
31 : USE qs_density_matrices, ONLY: calculate_density_matrix
32 : USE qs_energy_matrix_w, ONLY: qs_energies_compute_w
33 : USE qs_environment_types, ONLY: get_qs_env,&
34 : qs_environment_type
35 : USE qs_force_types, ONLY: allocate_qs_force,&
36 : deallocate_qs_force,&
37 : qs_force_type,&
38 : replicate_qs_force,&
39 : zero_qs_force
40 : USE qs_kind_types, ONLY: get_qs_kind,&
41 : qs_kind_type
42 : USE qs_ks_methods, ONLY: qs_ks_allocate_basics,&
43 : qs_ks_update_qs_env
44 : USE qs_ks_types, ONLY: get_ks_env,&
45 : qs_ks_env_type,&
46 : set_ks_env
47 : USE qs_mixing_utils, ONLY: mixing_allocate
48 : USE qs_mo_methods, ONLY: make_basis_sm
49 : USE qs_mo_types, ONLY: get_mo_set,&
50 : mo_set_type
51 : USE qs_neighbor_lists, ONLY: build_qs_neighbor_lists
52 : USE qs_rho_methods, ONLY: qs_rho_update_rho
53 : USE qs_rho_types, ONLY: qs_rho_get,&
54 : qs_rho_set,&
55 : qs_rho_type
56 : USE qs_scf_types, ONLY: qs_scf_env_type
57 : USE qs_subsys_types, ONLY: qs_subsys_set,&
58 : qs_subsys_type
59 : USE qs_update_s_mstruct, ONLY: qs_env_update_s_mstruct
60 : #include "./base/base_uses.f90"
61 :
62 : IMPLICIT NONE
63 :
64 : PRIVATE
65 :
66 : ! *** Global parameters ***
67 :
68 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_basis_gradient'
69 :
70 : ! *** Public subroutines ***
71 :
72 : PUBLIC :: qs_basis_center_gradient, qs_update_basis_center_pos, &
73 : return_basis_center_gradient_norm
74 :
75 : CONTAINS
76 :
77 : ! *****************************************************************************
78 : ! for development of floating basis functions
79 : ! *****************************************************************************
80 : !> \brief ...
81 : !> \param qs_env ...
82 : ! **************************************************************************************************
83 0 : SUBROUTINE qs_basis_center_gradient(qs_env)
84 :
85 : TYPE(qs_environment_type), POINTER :: qs_env
86 :
87 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_basis_center_gradient'
88 :
89 : INTEGER :: handle, i, iatom, ikind, img, ispin, &
90 : natom, nimg, nkind, nspin
91 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, natom_of_kind
92 : LOGICAL :: floating, has_unit_metric
93 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient
94 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
95 : TYPE(cp_logger_type), POINTER :: logger
96 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, matrix_w_kp
97 : TYPE(dbcsr_type), POINTER :: matrix
98 : TYPE(dft_control_type), POINTER :: dft_control
99 : TYPE(mp_para_env_type), POINTER :: para_env
100 0 : TYPE(qs_force_type), DIMENSION(:), POINTER :: basis_force, force, qs_force
101 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
102 : TYPE(qs_ks_env_type), POINTER :: ks_env
103 : TYPE(qs_scf_env_type), POINTER :: scf_env
104 : TYPE(qs_subsys_type), POINTER :: subsys
105 :
106 0 : CALL timeset(routineN, handle)
107 0 : NULLIFY (logger)
108 0 : logger => cp_get_default_logger()
109 :
110 : ! get the gradient array
111 0 : CALL get_qs_env(qs_env, scf_env=scf_env, natom=natom)
112 0 : IF (ASSOCIATED(scf_env%floating_basis%gradient)) THEN
113 0 : gradient => scf_env%floating_basis%gradient
114 0 : CPASSERT(SIZE(gradient) == 3*natom)
115 : ELSE
116 0 : ALLOCATE (gradient(3, natom))
117 0 : scf_env%floating_basis%gradient => gradient
118 : END IF
119 0 : gradient = 0.0_dp
120 :
121 : ! init the force environment
122 0 : CALL get_qs_env(qs_env, force=force, subsys=subsys, atomic_kind_set=atomic_kind_set)
123 0 : IF (ASSOCIATED(force)) THEN
124 0 : qs_force => force
125 : ELSE
126 0 : NULLIFY (qs_force)
127 : END IF
128 : ! Allocate the force data structure
129 0 : nkind = SIZE(atomic_kind_set)
130 0 : ALLOCATE (natom_of_kind(nkind))
131 0 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
132 0 : CALL allocate_qs_force(basis_force, natom_of_kind)
133 0 : DEALLOCATE (natom_of_kind)
134 0 : CALL qs_subsys_set(subsys, force=basis_force)
135 0 : CALL zero_qs_force(basis_force)
136 :
137 : ! get atom mapping
138 0 : ALLOCATE (atom_of_kind(natom), kind_of(natom))
139 0 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
140 :
141 : ! allocate energy weighted density matrices, if needed
142 0 : CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
143 0 : IF (.NOT. has_unit_metric) THEN
144 0 : CALL get_qs_env(qs_env, matrix_w_kp=matrix_w_kp)
145 0 : IF (.NOT. ASSOCIATED(matrix_w_kp)) THEN
146 0 : NULLIFY (matrix_w_kp)
147 0 : CALL get_qs_env(qs_env, ks_env=ks_env, matrix_s_kp=matrix_s, dft_control=dft_control)
148 0 : nspin = dft_control%nspins
149 0 : nimg = dft_control%nimages
150 0 : matrix => matrix_s(1, 1)%matrix
151 0 : CALL dbcsr_allocate_matrix_set(matrix_w_kp, nspin, nimg)
152 0 : DO ispin = 1, nspin
153 0 : DO img = 1, nimg
154 0 : ALLOCATE (matrix_w_kp(ispin, img)%matrix)
155 0 : CALL dbcsr_copy(matrix_w_kp(ispin, img)%matrix, matrix, name="W MATRIX")
156 0 : CALL dbcsr_set(matrix_w_kp(ispin, img)%matrix, 0.0_dp)
157 : END DO
158 : END DO
159 0 : CALL set_ks_env(ks_env, matrix_w_kp=matrix_w_kp)
160 : END IF
161 : END IF
162 : ! time to compute the w matrix
163 0 : CALL qs_energies_compute_w(qs_env, .TRUE.)
164 :
165 : ! core hamiltonian forces
166 0 : CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.TRUE.)
167 : ! Compute grid-based forces
168 0 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.TRUE.)
169 :
170 : ! replicate forces
171 0 : CALL replicate_qs_force(basis_force, para_env)
172 0 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
173 0 : DO iatom = 1, natom
174 0 : ikind = kind_of(iatom)
175 0 : i = atom_of_kind(iatom)
176 0 : CALL get_qs_kind(qs_kind_set(ikind), floating=floating)
177 0 : IF (floating) gradient(1:3, iatom) = -basis_force(ikind)%total(1:3, i)
178 : END DO
179 : ! clean up force environment and reinitialize qs_force
180 0 : IF (ASSOCIATED(basis_force)) CALL deallocate_qs_force(basis_force)
181 0 : CALL qs_subsys_set(subsys, force=qs_force)
182 0 : CALL get_qs_env(qs_env, ks_env=ks_env)
183 0 : CALL set_ks_env(ks_env, forces_up_to_date=.FALSE.)
184 :
185 0 : DEALLOCATE (atom_of_kind, kind_of)
186 :
187 0 : CALL timestop(handle)
188 :
189 0 : END SUBROUTINE qs_basis_center_gradient
190 :
191 : ! *****************************************************************************
192 : !> \brief ... returns the norm of the gradient vector, taking only floating
193 : !> components into account
194 : !> \param qs_env ...
195 : !> \return ...
196 : ! **************************************************************************************************
197 0 : FUNCTION return_basis_center_gradient_norm(qs_env) RESULT(norm)
198 :
199 : TYPE(qs_environment_type), POINTER :: qs_env
200 : REAL(KIND=dp) :: norm
201 :
202 : INTEGER :: iatom, ikind, natom, nfloat
203 : LOGICAL :: floating
204 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient
205 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
206 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
207 : TYPE(qs_scf_env_type), POINTER :: scf_env
208 :
209 0 : norm = 0.0_dp
210 0 : CALL get_qs_env(qs_env, scf_env=scf_env, particle_set=particle_set, qs_kind_set=qs_kind_set)
211 0 : gradient => scf_env%floating_basis%gradient
212 0 : natom = SIZE(particle_set)
213 0 : nfloat = 0
214 0 : DO iatom = 1, natom
215 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
216 0 : CALL get_qs_kind(qs_kind_set(ikind), floating=floating)
217 0 : IF (floating) THEN
218 0 : nfloat = nfloat + 1
219 0 : norm = norm + SUM(ABS(gradient(1:3, iatom)))
220 : END IF
221 : END DO
222 0 : IF (nfloat > 0) THEN
223 0 : norm = norm/(3.0_dp*REAL(nfloat, KIND=dp))
224 : END IF
225 :
226 0 : END FUNCTION return_basis_center_gradient_norm
227 :
228 : ! *****************************************************************************
229 : !> \brief move atoms with kind float according to gradient
230 : !> \param qs_env ...
231 : ! **************************************************************************************************
232 0 : SUBROUTINE qs_update_basis_center_pos(qs_env)
233 :
234 : TYPE(qs_environment_type), POINTER :: qs_env
235 :
236 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_update_basis_center_pos'
237 :
238 : INTEGER :: handle, iatom, ikind, natom
239 : LOGICAL :: floating
240 : REAL(KIND=dp) :: alpha
241 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gradient
242 : TYPE(cp_logger_type), POINTER :: logger
243 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
244 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
245 : TYPE(qs_scf_env_type), POINTER :: scf_env
246 :
247 0 : CALL timeset(routineN, handle)
248 0 : NULLIFY (logger)
249 0 : logger => cp_get_default_logger()
250 :
251 : ! update positions
252 0 : CALL get_qs_env(qs_env, scf_env=scf_env, particle_set=particle_set, qs_kind_set=qs_kind_set)
253 0 : gradient => scf_env%floating_basis%gradient
254 0 : natom = SIZE(particle_set)
255 0 : alpha = 0.50_dp
256 0 : DO iatom = 1, natom
257 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
258 0 : CALL get_qs_kind(qs_kind_set(ikind), floating=floating)
259 0 : IF (floating) THEN
260 0 : particle_set(iatom)%r(1:3) = particle_set(iatom)%r(1:3) + alpha*gradient(1:3, iatom)
261 : END IF
262 : END DO
263 :
264 0 : CALL qs_basis_reinit_energy(qs_env)
265 :
266 0 : CALL timestop(handle)
267 :
268 0 : END SUBROUTINE qs_update_basis_center_pos
269 :
270 : ! *****************************************************************************
271 : !> \brief rebuilds the structures after a floating basis update
272 : !> \param qs_env ...
273 : !> \par History
274 : !> 05.2016 created [JGH]
275 : !> \author JGH
276 : ! **************************************************************************************************
277 0 : SUBROUTINE qs_basis_reinit_energy(qs_env)
278 : TYPE(qs_environment_type), POINTER :: qs_env
279 :
280 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_basis_reinit_energy'
281 :
282 : INTEGER :: handle, ispin, nmo
283 : LOGICAL :: ks_is_complex
284 : TYPE(cp_fm_type), POINTER :: mo_coeff
285 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_kp
286 : TYPE(dft_control_type), POINTER :: dft_control
287 0 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
288 : TYPE(mp_para_env_type), POINTER :: para_env
289 : TYPE(qs_ks_env_type), POINTER :: ks_env
290 : TYPE(qs_rho_type), POINTER :: rho
291 : TYPE(qs_scf_env_type), POINTER :: scf_env
292 : TYPE(section_vals_type), POINTER :: input
293 :
294 0 : CALL timeset(routineN, handle)
295 :
296 0 : NULLIFY (input, para_env, ks_env)
297 : ! rebuild neighbor lists
298 0 : CALL get_qs_env(qs_env, input=input, para_env=para_env, ks_env=ks_env)
299 : CALL build_qs_neighbor_lists(qs_env, para_env, molecular=.FALSE., &
300 0 : force_env_section=input)
301 : ! update core hamiltonian
302 0 : CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.FALSE.)
303 : ! update structures
304 0 : CALL qs_env_update_s_mstruct(qs_env)
305 : ! KS matrices
306 0 : CALL get_ks_env(ks_env, complex_ks=ks_is_complex)
307 0 : CALL qs_ks_allocate_basics(qs_env, is_complex=ks_is_complex)
308 : ! reinit SCF task matrices
309 0 : NULLIFY (scf_env)
310 0 : CALL get_qs_env(qs_env, scf_env=scf_env, dft_control=dft_control)
311 0 : IF (scf_env%mixing_method > 0) THEN
312 : CALL mixing_allocate(qs_env, scf_env%mixing_method, scf_env%p_mix_new, &
313 : scf_env%p_delta, dft_control%nspins, &
314 0 : scf_env%mixing_store)
315 : ELSE
316 0 : NULLIFY (scf_env%p_mix_new)
317 : END IF
318 0 : CALL get_qs_env(qs_env, mos=mos, rho=rho, matrix_s_kp=matrix_s_kp)
319 0 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
320 0 : DO ispin = 1, SIZE(mos)
321 0 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
322 : ! reorthogonalize MOs
323 0 : CALL make_basis_sm(mo_coeff, nmo, matrix_s_kp(1, 1)%matrix)
324 : ! update density matrix
325 0 : CALL calculate_density_matrix(mos(ispin), rho_ao_kp(ispin, 1)%matrix)
326 : END DO
327 : CALL qs_rho_set(rho, rho_r_valid=.FALSE., drho_r_valid=.FALSE., rho_g_valid=.FALSE., &
328 0 : drho_g_valid=.FALSE., tau_r_valid=.FALSE., tau_g_valid=.FALSE.)
329 0 : CALL qs_rho_update_rho(rho, qs_env)
330 :
331 0 : CALL timestop(handle)
332 :
333 0 : END SUBROUTINE qs_basis_reinit_energy
334 :
335 : END MODULE qs_basis_gradient
|