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 Utility subroutines for CDFT calculations
10 : !> \par History
11 : !> separated from et_coupling [03.2017]
12 : !> \author Nico Holmberg [03.2017]
13 : ! **************************************************************************************************
14 : MODULE qs_cdft_utils
15 : USE ao_util, ONLY: exp_radius_very_extended
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind
18 : USE bibliography, ONLY: Becke1988b,&
19 : Holmberg2017,&
20 : Holmberg2018,&
21 : cite_reference
22 : USE cell_types, ONLY: cell_type,&
23 : pbc
24 : USE cp_control_types, ONLY: dft_control_type,&
25 : qs_control_type
26 : USE cp_log_handling, ONLY: cp_get_default_logger,&
27 : cp_logger_type,&
28 : cp_to_string
29 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
30 : cp_print_key_unit_nr
31 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
32 : USE cp_units, ONLY: cp_unit_from_cp2k
33 : USE grid_api, ONLY: GRID_FUNC_AB,&
34 : collocate_pgf_product
35 : USE hirshfeld_methods, ONLY: create_shape_function
36 : USE hirshfeld_types, ONLY: create_hirshfeld_type,&
37 : hirshfeld_type,&
38 : set_hirshfeld_info
39 : USE input_constants, ONLY: &
40 : becke_cutoff_element, becke_cutoff_global, cdft_charge_constraint, &
41 : outer_scf_becke_constraint, outer_scf_cdft_constraint, outer_scf_hirshfeld_constraint, &
42 : outer_scf_none, radius_user, shape_function_gaussian
43 : USE input_section_types, ONLY: section_get_ivals,&
44 : section_vals_get,&
45 : section_vals_get_subs_vals,&
46 : section_vals_type,&
47 : section_vals_val_get
48 : USE kinds, ONLY: default_path_length,&
49 : dp
50 : USE memory_utilities, ONLY: reallocate
51 : USE message_passing, ONLY: mp_para_env_type
52 : USE outer_scf_control_types, ONLY: outer_scf_read_parameters
53 : USE particle_list_types, ONLY: particle_list_type
54 : USE particle_types, ONLY: particle_type
55 : USE pw_env_types, ONLY: pw_env_get,&
56 : pw_env_type
57 : USE pw_methods, ONLY: pw_zero
58 : USE pw_pool_types, ONLY: pw_pool_type
59 : USE qs_cdft_types, ONLY: becke_constraint_type,&
60 : cdft_control_type,&
61 : cdft_group_type,&
62 : hirshfeld_constraint_type
63 : USE qs_environment_types, ONLY: get_qs_env,&
64 : qs_environment_type
65 : USE qs_kind_types, ONLY: get_qs_kind,&
66 : qs_kind_type
67 : USE qs_scf_output, ONLY: qs_scf_cdft_constraint_info
68 : USE qs_subsys_types, ONLY: qs_subsys_get,&
69 : qs_subsys_type
70 : USE realspace_grid_types, ONLY: realspace_grid_type,&
71 : rs_grid_zero,&
72 : transfer_rs2pw
73 : #include "./base/base_uses.f90"
74 :
75 : IMPLICIT NONE
76 :
77 : PRIVATE
78 :
79 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cdft_utils'
80 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
81 :
82 : ! *** Public subroutines ***
83 : PUBLIC :: becke_constraint_init, read_becke_section, read_cdft_control_section
84 : PUBLIC :: hfun_scale, hirshfeld_constraint_init, cdft_constraint_print, &
85 : cdft_print_hirshfeld_density, cdft_print_weight_function
86 :
87 : CONTAINS
88 :
89 : ! **************************************************************************************************
90 : !> \brief Initializes the Becke constraint environment
91 : !> \param qs_env the qs_env where to build the constraint
92 : !> \par History
93 : !> Created 01.2007 [fschiff]
94 : !> Extended functionality 12/15-12/16 [Nico Holmberg]
95 : ! **************************************************************************************************
96 202 : SUBROUTINE becke_constraint_init(qs_env)
97 : TYPE(qs_environment_type), POINTER :: qs_env
98 :
99 : CHARACTER(len=*), PARAMETER :: routineN = 'becke_constraint_init'
100 :
101 : CHARACTER(len=2) :: element_symbol
102 : INTEGER :: atom_a, bounds(2), handle, i, iatom, iex, igroup, ikind, ip, ithread, iw, j, &
103 : jatom, katom, natom, nkind, npme, nthread, numexp, unit_nr
104 : INTEGER, DIMENSION(2, 3) :: bo
105 202 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores, stride
106 : LOGICAL :: build, in_memory, mpi_io
107 202 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_constraint
108 : REAL(KIND=dp) :: alpha, chi, coef, eps_cavity, ircov, &
109 : jrcov, radius, uij
110 : REAL(KIND=dp), DIMENSION(3) :: cell_v, dist_vec, r, r1, ra
111 202 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii_list
112 202 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
113 202 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
114 : TYPE(becke_constraint_type), POINTER :: becke_control
115 : TYPE(cdft_control_type), POINTER :: cdft_control
116 202 : TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
117 : TYPE(cell_type), POINTER :: cell
118 : TYPE(cp_logger_type), POINTER :: logger
119 : TYPE(dft_control_type), POINTER :: dft_control
120 : TYPE(hirshfeld_type), POINTER :: cavity_env
121 : TYPE(mp_para_env_type), POINTER :: para_env
122 : TYPE(particle_list_type), POINTER :: particles
123 202 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
124 : TYPE(pw_env_type), POINTER :: pw_env
125 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
126 202 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
127 : TYPE(qs_subsys_type), POINTER :: subsys
128 : TYPE(realspace_grid_type), POINTER :: rs_cavity
129 : TYPE(section_vals_type), POINTER :: cdft_constraint_section
130 :
131 202 : NULLIFY (cores, stride, atom_list, cell, para_env, dft_control, &
132 202 : particle_set, logger, cdft_constraint_section, qs_kind_set, &
133 202 : particles, subsys, pab, pw_env, rs_cavity, cavity_env, &
134 202 : auxbas_pw_pool, atomic_kind_set, group, radii_list, cdft_control)
135 404 : logger => cp_get_default_logger()
136 202 : CALL timeset(routineN, handle)
137 : CALL get_qs_env(qs_env, &
138 : cell=cell, &
139 : particle_set=particle_set, &
140 : natom=natom, &
141 : dft_control=dft_control, &
142 202 : para_env=para_env)
143 202 : cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
144 202 : iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog")
145 202 : cdft_control => dft_control%qs_control%cdft_control
146 202 : becke_control => cdft_control%becke_control
147 202 : group => cdft_control%group
148 202 : in_memory = .FALSE.
149 202 : IF (cdft_control%save_pot) THEN
150 72 : in_memory = becke_control%in_memory
151 : END IF
152 202 : IF (becke_control%cavity_confine) THEN
153 522 : ALLOCATE (is_constraint(natom))
154 174 : is_constraint = .FALSE.
155 516 : DO i = 1, cdft_control%natoms
156 : ! Notice that here is_constraint=.TRUE. also for dummy atoms to properly compute their Becke charges
157 : ! A subsequent check (atom_in_group) ensures that the gradients of these dummy atoms are correct
158 516 : is_constraint(cdft_control%atoms(i)) = .TRUE.
159 : END DO
160 : END IF
161 202 : eps_cavity = becke_control%eps_cavity
162 : ! Setup atomic radii for adjusting cell boundaries
163 202 : IF (becke_control%adjust) THEN
164 118 : IF (.NOT. ASSOCIATED(becke_control%radii)) THEN
165 94 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
166 94 : IF (.NOT. SIZE(atomic_kind_set) == SIZE(becke_control%radii_tmp)) THEN
167 : CALL cp_abort(__LOCATION__, &
168 : "Length of keyword BECKE_CONSTRAINT\ATOMIC_RADII does not "// &
169 0 : "match number of atomic kinds in the input coordinate file.")
170 : END IF
171 282 : ALLOCATE (becke_control%radii(SIZE(atomic_kind_set)))
172 282 : becke_control%radii(:) = becke_control%radii_tmp(:)
173 94 : DEALLOCATE (becke_control%radii_tmp)
174 : END IF
175 : END IF
176 : ! Setup cutoff scheme
177 202 : IF (.NOT. ASSOCIATED(becke_control%cutoffs)) THEN
178 162 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
179 486 : ALLOCATE (becke_control%cutoffs(natom))
180 288 : SELECT CASE (becke_control%cutoff_type)
181 : CASE (becke_cutoff_global)
182 390 : becke_control%cutoffs(:) = becke_control%rglobal
183 : CASE (becke_cutoff_element)
184 36 : IF (.NOT. SIZE(atomic_kind_set) == SIZE(becke_control%cutoffs_tmp)) THEN
185 : CALL cp_abort(__LOCATION__, &
186 : "Length of keyword BECKE_CONSTRAINT\ELEMENT_CUTOFFS does not "// &
187 0 : "match number of atomic kinds in the input coordinate file.")
188 : END IF
189 108 : DO ikind = 1, SIZE(atomic_kind_set)
190 72 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=katom, atom_list=atom_list)
191 202 : DO iatom = 1, katom
192 94 : atom_a = atom_list(iatom)
193 166 : becke_control%cutoffs(atom_a) = becke_control%cutoffs_tmp(ikind)
194 : END DO
195 : END DO
196 198 : DEALLOCATE (becke_control%cutoffs_tmp)
197 : END SELECT
198 : END IF
199 : ! Zero weight functions
200 432 : DO igroup = 1, SIZE(group)
201 432 : CALL pw_zero(group(igroup)%weight)
202 : END DO
203 202 : IF (cdft_control%atomic_charges) THEN
204 310 : DO iatom = 1, cdft_control%natoms
205 310 : CALL pw_zero(cdft_control%charge(iatom))
206 : END DO
207 : END IF
208 : ! Allocate storage for cell adjustment coefficients and needed distance vectors
209 202 : build = .FALSE.
210 202 : IF (becke_control%adjust .AND. .NOT. ASSOCIATED(becke_control%aij)) THEN
211 376 : ALLOCATE (becke_control%aij(natom, natom))
212 94 : build = .TRUE.
213 : END IF
214 202 : IF (becke_control%vector_buffer%store_vectors) THEN
215 606 : ALLOCATE (becke_control%vector_buffer%distances(natom))
216 606 : ALLOCATE (becke_control%vector_buffer%distance_vecs(3, natom))
217 388 : IF (in_memory) ALLOCATE (becke_control%vector_buffer%pair_dist_vecs(3, natom, natom))
218 404 : ALLOCATE (becke_control%vector_buffer%position_vecs(3, natom))
219 : END IF
220 808 : ALLOCATE (becke_control%vector_buffer%R12(natom, natom))
221 : ! Calculate pairwise distances between each atom pair
222 808 : DO i = 1, 3
223 808 : cell_v(i) = cell%hmat(i, i)
224 : END DO
225 438 : DO iatom = 1, natom - 1
226 708 : DO jatom = iatom + 1, natom
227 1080 : r = particle_set(iatom)%r
228 1080 : r1 = particle_set(jatom)%r
229 1080 : DO i = 1, 3
230 810 : r(i) = MODULO(r(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
231 1080 : r1(i) = MODULO(r1(i), cell%hmat(i, i)) - cell%hmat(i, i)/2._dp
232 : END DO
233 1080 : dist_vec = (r - r1) - ANINT((r - r1)/cell_v)*cell_v
234 : ! Store pbc corrected position and pairwise distance vectors for later reuse
235 270 : IF (becke_control%vector_buffer%store_vectors) THEN
236 1080 : becke_control%vector_buffer%position_vecs(:, iatom) = r(:)
237 876 : IF (iatom == 1 .AND. jatom == natom) becke_control%vector_buffer%position_vecs(:, jatom) = r1(:)
238 270 : IF (in_memory) THEN
239 248 : becke_control%vector_buffer%pair_dist_vecs(:, iatom, jatom) = dist_vec(:)
240 248 : becke_control%vector_buffer%pair_dist_vecs(:, jatom, iatom) = -dist_vec(:)
241 : END IF
242 : END IF
243 1080 : becke_control%vector_buffer%R12(iatom, jatom) = NORM2(dist_vec)
244 270 : becke_control%vector_buffer%R12(jatom, iatom) = becke_control%vector_buffer%R12(iatom, jatom)
245 : ! Set up heteronuclear cell partitioning using user defined radii
246 506 : IF (build) THEN
247 150 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, kind_number=ikind)
248 150 : ircov = becke_control%radii(ikind)
249 150 : CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, kind_number=ikind)
250 150 : jrcov = becke_control%radii(ikind)
251 150 : IF (ircov /= jrcov) THEN
252 122 : chi = ircov/jrcov
253 122 : uij = (chi - 1.0_dp)/(chi + 1.0_dp)
254 122 : becke_control%aij(iatom, jatom) = uij/(uij**2 - 1.0_dp)
255 122 : IF (becke_control%aij(iatom, jatom) > 0.5_dp) THEN
256 0 : becke_control%aij(iatom, jatom) = 0.5_dp
257 122 : ELSE IF (becke_control%aij(iatom, jatom) < -0.5_dp) THEN
258 0 : becke_control%aij(iatom, jatom) = -0.5_dp
259 : END IF
260 : ELSE
261 28 : becke_control%aij(iatom, jatom) = 0.0_dp
262 : END IF
263 : ! Note change of sign
264 150 : becke_control%aij(jatom, iatom) = -becke_control%aij(iatom, jatom)
265 : END IF
266 : END DO
267 : END DO
268 : ! Dump some additional information about the calculation
269 202 : IF (cdft_control%first_iteration) THEN
270 162 : IF (iw > 0) THEN
271 : WRITE (iw, '(/,T3,A)') &
272 81 : '----------------------- Becke atomic parameters ------------------------'
273 81 : IF (becke_control%adjust) THEN
274 : WRITE (iw, '(T3,A)') &
275 47 : 'Atom Element Cutoff (angstrom) CDFT Radius (angstrom)'
276 155 : DO iatom = 1, natom
277 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol, &
278 108 : kind_number=ikind)
279 108 : ircov = cp_unit_from_cp2k(becke_control%radii(ikind), "angstrom")
280 : WRITE (iw, "(i6,T15,A2,T37,F8.3,T67,F8.3)") &
281 108 : iatom, ADJUSTR(element_symbol), cp_unit_from_cp2k(becke_control%cutoffs(iatom), "angstrom"), &
282 371 : ircov
283 : END DO
284 : ELSE
285 : WRITE (iw, '(T3,A)') &
286 34 : 'Atom Element Cutoff (angstrom)'
287 105 : DO iatom = 1, natom
288 71 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol)
289 : WRITE (iw, "(i7,T15,A2,T37,F8.3)") &
290 105 : iatom, ADJUSTR(element_symbol), cp_unit_from_cp2k(becke_control%cutoffs(iatom), "angstrom")
291 : END DO
292 : END IF
293 : WRITE (iw, '(T3,A)') &
294 81 : '------------------------------------------------------------------------'
295 : WRITE (iw, '(/,T3,A,T60)') &
296 81 : '----------------------- Becke group definitions ------------------------'
297 172 : DO igroup = 1, SIZE(group)
298 91 : IF (igroup > 1) WRITE (iw, '(T3,A)') ' '
299 : WRITE (iw, '(T5,A,I5,A,I5)') &
300 91 : 'Atomic group', igroup, ' of ', SIZE(group)
301 91 : WRITE (iw, '(T5,A)') 'Atom Element Coefficient'
302 349 : DO ip = 1, SIZE(group(igroup)%atoms)
303 177 : iatom = group(igroup)%atoms(ip)
304 177 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol)
305 268 : WRITE (iw, '(i8,T16,A2,T23,F8.3)') iatom, ADJUSTR(element_symbol), group(igroup)%coeff(ip)
306 : END DO
307 : END DO
308 : WRITE (iw, '(T3,A)') &
309 81 : '------------------------------------------------------------------------'
310 : END IF
311 162 : cdft_control%first_iteration = .FALSE.
312 : END IF
313 : ! Setup cavity confinement using spherical Gaussians
314 202 : IF (becke_control%cavity_confine) THEN
315 174 : cavity_env => becke_control%cavity_env
316 174 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, pw_env=pw_env, qs_kind_set=qs_kind_set)
317 174 : CPASSERT(ASSOCIATED(qs_kind_set))
318 174 : nkind = SIZE(qs_kind_set)
319 : ! Setup the Gaussian shape function
320 174 : IF (.NOT. ASSOCIATED(cavity_env%kind_shape_fn)) THEN
321 138 : IF (ASSOCIATED(becke_control%radii)) THEN
322 276 : ALLOCATE (radii_list(SIZE(becke_control%radii)))
323 276 : DO ikind = 1, SIZE(becke_control%radii)
324 276 : IF (cavity_env%use_bohr) THEN
325 4 : radii_list(ikind) = becke_control%radii(ikind)
326 : ELSE
327 180 : radii_list(ikind) = cp_unit_from_cp2k(becke_control%radii(ikind), "angstrom")
328 : END IF
329 : END DO
330 : END IF
331 : CALL create_shape_function(cavity_env, qs_kind_set, atomic_kind_set, &
332 : radius=becke_control%rcavity, &
333 138 : radii_list=radii_list)
334 138 : IF (ASSOCIATED(radii_list)) THEN
335 92 : DEALLOCATE (radii_list)
336 : END IF
337 : END IF
338 : ! Form cavity by summing isolated Gaussian densities over constraint atoms
339 174 : NULLIFY (rs_cavity)
340 174 : CALL pw_env_get(pw_env, auxbas_rs_grid=rs_cavity, auxbas_pw_pool=auxbas_pw_pool)
341 174 : CALL rs_grid_zero(rs_cavity)
342 174 : ALLOCATE (pab(1, 1))
343 174 : nthread = 1
344 174 : ithread = 0
345 482 : DO ikind = 1, SIZE(atomic_kind_set)
346 308 : numexp = cavity_env%kind_shape_fn(ikind)%numexp
347 308 : IF (numexp <= 0) CYCLE
348 308 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=katom, atom_list=atom_list)
349 924 : ALLOCATE (cores(katom))
350 616 : DO iex = 1, numexp
351 308 : alpha = cavity_env%kind_shape_fn(ikind)%zet(iex)
352 308 : coef = cavity_env%kind_shape_fn(ikind)%coef(iex)
353 308 : npme = 0
354 688 : cores = 0
355 688 : DO iatom = 1, katom
356 688 : IF (rs_cavity%desc%parallel .AND. .NOT. rs_cavity%desc%distributed) THEN
357 : ! replicated realspace grid, split the atoms up between procs
358 380 : IF (MODULO(iatom, rs_cavity%desc%group_size) == rs_cavity%desc%my_pos) THEN
359 190 : npme = npme + 1
360 190 : cores(npme) = iatom
361 : END IF
362 : ELSE
363 0 : npme = npme + 1
364 0 : cores(npme) = iatom
365 : END IF
366 : END DO
367 806 : DO j = 1, npme
368 190 : iatom = cores(j)
369 190 : atom_a = atom_list(iatom)
370 190 : pab(1, 1) = coef
371 190 : IF (becke_control%vector_buffer%store_vectors) THEN
372 760 : ra(:) = becke_control%vector_buffer%position_vecs(:, atom_a) + cell_v(:)/2._dp
373 : ELSE
374 0 : ra(:) = pbc(particle_set(atom_a)%r, cell)
375 : END IF
376 498 : IF (is_constraint(atom_a)) THEN
377 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
378 : ra=ra, rb=ra, rp=ra, zetp=alpha, &
379 : eps=dft_control%qs_control%eps_rho_rspace, &
380 : pab=pab, o1=0, o2=0, & ! without map_consistent
381 171 : prefactor=1.0_dp, cutoff=0.0_dp)
382 :
383 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
384 : [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, &
385 : pab, 0, 0, rs_cavity, &
386 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
387 171 : use_subpatch=.TRUE., subpatch_pattern=0)
388 : END IF
389 : END DO
390 : END DO
391 790 : DEALLOCATE (cores)
392 : END DO
393 174 : DEALLOCATE (pab)
394 174 : CALL auxbas_pw_pool%create_pw(becke_control%cavity)
395 174 : CALL transfer_rs2pw(rs_cavity, becke_control%cavity)
396 : ! Grid points where the Gaussian density falls below eps_cavity are ignored
397 : ! We can calculate the smallest/largest values along z-direction outside
398 : ! which the cavity is zero at every point (x, y)
399 : ! If gradients are needed storage needs to be allocated only for grid points within
400 : ! these bounds
401 174 : IF (in_memory .OR. cdft_control%save_pot) THEN
402 64 : CALL hfun_zero(becke_control%cavity%array, eps_cavity, just_bounds=.TRUE., bounds=bounds)
403 : ! Save bounds (first nonzero grid point indices)
404 640 : bo = group(1)%weight%pw_grid%bounds_local
405 64 : IF (bounds(2) < bo(2, 3)) THEN
406 8 : bounds(2) = bounds(2) - 1
407 : ELSE
408 56 : bounds(2) = bo(2, 3)
409 : END IF
410 64 : IF (bounds(1) > bo(1, 3)) THEN
411 : ! In the special case bounds(1) == bounds(2) == bo(2, 3), after this check
412 : ! bounds(1) > bounds(2) and the subsequent gradient allocation (:, :, :, bounds(1):bounds(2))
413 : ! will correctly allocate a 0-sized array
414 8 : bounds(1) = bounds(1) + 1
415 : ELSE
416 56 : bounds(1) = bo(1, 3)
417 : END IF
418 302 : becke_control%confine_bounds = bounds
419 : END IF
420 : ! Optional printing of cavity (meant for testing, so options currently hardcoded...)
421 174 : IF (becke_control%print_cavity) THEN
422 2 : CALL hfun_zero(becke_control%cavity%array, eps_cavity, just_bounds=.FALSE.)
423 2 : ALLOCATE (stride(3))
424 8 : stride = [2, 2, 2]
425 2 : mpi_io = .TRUE.
426 : ! Note PROGRAM_RUN_INFO section neeeds to be active!
427 : unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", &
428 : middle_name="BECKE_CAVITY", &
429 : extension=".cube", file_position="REWIND", &
430 2 : log_filename=.FALSE., mpi_io=mpi_io)
431 2 : IF (para_env%is_source() .AND. unit_nr < 1) THEN
432 : CALL cp_abort(__LOCATION__, &
433 0 : "Please turn on PROGRAM_RUN_INFO to print cavity")
434 : END IF
435 2 : CALL get_qs_env(qs_env, subsys=subsys)
436 2 : CALL qs_subsys_get(subsys, particles=particles)
437 2 : CALL cp_pw_to_cube(becke_control%cavity, unit_nr, "CAVITY", particles=particles, stride=stride, mpi_io=mpi_io)
438 2 : CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io)
439 2 : DEALLOCATE (stride)
440 : END IF
441 : END IF
442 202 : IF (ALLOCATED(is_constraint)) THEN
443 174 : DEALLOCATE (is_constraint)
444 : END IF
445 202 : CALL timestop(handle)
446 :
447 404 : END SUBROUTINE becke_constraint_init
448 :
449 : ! **************************************************************************************************
450 : !> \brief reads the input parameters specific to Becke-based CDFT constraints
451 : !> \param cdft_control the cdft_control which holds the Becke control type
452 : !> \param becke_section the input section containing Becke constraint information
453 : !> \par History
454 : !> Created 01.2007 [fschiff]
455 : !> Merged Becke into CDFT 09.2018 [Nico Holmberg]
456 : !> \author Nico Holmberg [09.2018]
457 : ! **************************************************************************************************
458 254 : SUBROUTINE read_becke_section(cdft_control, becke_section)
459 :
460 : TYPE(cdft_control_type), INTENT(INOUT) :: cdft_control
461 : TYPE(section_vals_type), POINTER :: becke_section
462 :
463 : INTEGER :: j
464 : LOGICAL :: exists
465 254 : REAL(KIND=dp), DIMENSION(:), POINTER :: rtmplist
466 : TYPE(becke_constraint_type), POINTER :: becke_control
467 :
468 254 : NULLIFY (rtmplist)
469 254 : becke_control => cdft_control%becke_control
470 0 : CPASSERT(ASSOCIATED(becke_control))
471 :
472 : ! Atomic size corrections
473 254 : CALL section_vals_val_get(becke_section, "ADJUST_SIZE", l_val=becke_control%adjust)
474 254 : IF (becke_control%adjust) THEN
475 156 : CALL section_vals_val_get(becke_section, "ATOMIC_RADII", explicit=exists)
476 156 : IF (.NOT. exists) CPABORT("Keyword ATOMIC_RADII is missing.")
477 156 : CALL section_vals_val_get(becke_section, "ATOMIC_RADII", r_vals=rtmplist)
478 156 : CPASSERT(SIZE(rtmplist) > 0)
479 468 : ALLOCATE (becke_control%radii_tmp(SIZE(rtmplist)))
480 624 : DO j = 1, SIZE(rtmplist)
481 468 : becke_control%radii_tmp(j) = rtmplist(j)
482 : END DO
483 : END IF
484 :
485 : ! Cutoff scheme
486 254 : CALL section_vals_val_get(becke_section, "CUTOFF_TYPE", i_val=becke_control%cutoff_type)
487 412 : SELECT CASE (becke_control%cutoff_type)
488 : CASE (becke_cutoff_global)
489 158 : CALL section_vals_val_get(becke_section, "GLOBAL_CUTOFF", r_val=becke_control%rglobal)
490 : CASE (becke_cutoff_element)
491 96 : CALL section_vals_val_get(becke_section, "ELEMENT_CUTOFF", r_vals=rtmplist)
492 96 : CPASSERT(SIZE(rtmplist) > 0)
493 288 : ALLOCATE (becke_control%cutoffs_tmp(SIZE(rtmplist)))
494 542 : DO j = 1, SIZE(rtmplist)
495 288 : becke_control%cutoffs_tmp(j) = rtmplist(j)
496 : END DO
497 : END SELECT
498 :
499 : ! Gaussian cavity confinement
500 254 : CALL section_vals_val_get(becke_section, "CAVITY_CONFINE", l_val=becke_control%cavity_confine)
501 254 : CALL section_vals_val_get(becke_section, "SHOULD_SKIP", l_val=becke_control%should_skip)
502 254 : CALL section_vals_val_get(becke_section, "IN_MEMORY", l_val=becke_control%in_memory)
503 254 : IF (cdft_control%becke_control%cavity_confine) THEN
504 222 : CALL section_vals_val_get(becke_section, "CAVITY_SHAPE", i_val=becke_control%cavity_shape)
505 222 : IF (becke_control%cavity_shape == radius_user .AND. .NOT. becke_control%adjust) THEN
506 : CALL cp_abort(__LOCATION__, &
507 0 : "Activate keyword ADJUST_SIZE to use cavity shape USER.")
508 : END IF
509 222 : CALL section_vals_val_get(becke_section, "CAVITY_RADIUS", r_val=becke_control%rcavity)
510 222 : CALL section_vals_val_get(becke_section, "EPS_CAVITY", r_val=becke_control%eps_cavity)
511 222 : CALL section_vals_val_get(becke_section, "CAVITY_PRINT", l_val=becke_control%print_cavity)
512 222 : CALL section_vals_val_get(becke_section, "CAVITY_USE_BOHR", l_val=becke_control%use_bohr)
513 222 : IF (.NOT. cdft_control%becke_control%use_bohr) THEN
514 220 : becke_control%rcavity = cp_unit_from_cp2k(becke_control%rcavity, "angstrom")
515 : END IF
516 222 : CALL create_hirshfeld_type(becke_control%cavity_env)
517 : CALL set_hirshfeld_info(becke_control%cavity_env, &
518 : shape_function_type=shape_function_gaussian, iterative=.FALSE., &
519 : radius_type=becke_control%cavity_shape, &
520 222 : use_bohr=becke_control%use_bohr)
521 : END IF
522 :
523 254 : CALL cite_reference(Becke1988b)
524 :
525 254 : END SUBROUTINE read_becke_section
526 :
527 : ! **************************************************************************************************
528 : !> \brief reads the input parameters needed to define CDFT constraints
529 : !> \param cdft_control the object which holds the CDFT control type
530 : !> \param cdft_control_section the input section containing CDFT constraint information
531 : !> \author Nico Holmberg [09.2018]
532 : ! **************************************************************************************************
533 282 : SUBROUTINE read_constraint_definitions(cdft_control, cdft_control_section)
534 :
535 : TYPE(cdft_control_type), INTENT(INOUT) :: cdft_control
536 : TYPE(section_vals_type), INTENT(INOUT), POINTER :: cdft_control_section
537 :
538 : INTEGER :: i, j, jj, k, n_rep, natoms, nvar, &
539 : tot_natoms
540 282 : INTEGER, DIMENSION(:), POINTER :: atomlist, dummylist, tmplist
541 : LOGICAL :: exists, is_duplicate
542 282 : REAL(KIND=dp), DIMENSION(:), POINTER :: rtmplist
543 : TYPE(section_vals_type), POINTER :: group_section
544 :
545 282 : NULLIFY (tmplist, rtmplist, atomlist, dummylist, group_section)
546 :
547 564 : group_section => section_vals_get_subs_vals(cdft_control_section, "ATOM_GROUP")
548 282 : CALL section_vals_get(group_section, n_repetition=nvar, explicit=exists)
549 282 : IF (.NOT. exists) CPABORT("Section ATOM_GROUP is missing.")
550 1150 : ALLOCATE (cdft_control%group(nvar))
551 282 : tot_natoms = 0
552 : ! Parse all ATOM_GROUP sections
553 586 : DO k = 1, nvar
554 : ! First determine how much storage is needed
555 304 : natoms = 0
556 304 : CALL section_vals_val_get(group_section, "ATOMS", i_rep_section=k, n_rep_val=n_rep)
557 608 : DO j = 1, n_rep
558 304 : CALL section_vals_val_get(group_section, "ATOMS", i_rep_section=k, i_rep_val=j, i_vals=tmplist)
559 304 : IF (SIZE(tmplist) < 1) THEN
560 0 : CPABORT("Each ATOM_GROUP must contain at least 1 atom.")
561 : END IF
562 608 : natoms = natoms + SIZE(tmplist)
563 : END DO
564 912 : ALLOCATE (cdft_control%group(k)%atoms(natoms))
565 912 : ALLOCATE (cdft_control%group(k)%coeff(natoms))
566 304 : NULLIFY (cdft_control%group(k)%weight)
567 304 : NULLIFY (cdft_control%group(k)%integrated)
568 304 : tot_natoms = tot_natoms + natoms
569 : ! Now parse
570 304 : jj = 0
571 608 : DO j = 1, n_rep
572 304 : CALL section_vals_val_get(group_section, "ATOMS", i_rep_section=k, i_rep_val=j, i_vals=tmplist)
573 1174 : DO i = 1, SIZE(tmplist)
574 566 : jj = jj + 1
575 870 : cdft_control%group(k)%atoms(jj) = tmplist(i)
576 : END DO
577 : END DO
578 304 : CALL section_vals_val_get(group_section, "COEFF", i_rep_section=k, n_rep_val=n_rep)
579 304 : jj = 0
580 608 : DO j = 1, n_rep
581 304 : CALL section_vals_val_get(group_section, "COEFF", i_rep_section=k, i_rep_val=j, r_vals=rtmplist)
582 1174 : DO i = 1, SIZE(rtmplist)
583 566 : jj = jj + 1
584 566 : IF (jj > natoms) CPABORT("Length of keywords ATOMS and COEFF must match.")
585 566 : IF (ABS(rtmplist(i)) /= 1.0_dp) CPABORT("Keyword COEFF accepts only values +/-1.0")
586 870 : cdft_control%group(k)%coeff(jj) = rtmplist(i)
587 : END DO
588 : END DO
589 304 : IF (jj < natoms) CPABORT("Length of keywords ATOMS and COEFF must match.")
590 : CALL section_vals_val_get(group_section, "CONSTRAINT_TYPE", i_rep_section=k, &
591 304 : i_val=cdft_control%group(k)%constraint_type)
592 : CALL section_vals_val_get(group_section, "FRAGMENT_CONSTRAINT", i_rep_section=k, &
593 304 : l_val=cdft_control%group(k)%is_fragment_constraint)
594 1194 : IF (cdft_control%group(k)%is_fragment_constraint) cdft_control%fragment_density = .TRUE.
595 : END DO
596 : ! Create a list containing all constraint atoms
597 846 : ALLOCATE (atomlist(tot_natoms))
598 848 : atomlist = -1
599 282 : jj = 0
600 586 : DO k = 1, nvar
601 1152 : DO j = 1, SIZE(cdft_control%group(k)%atoms)
602 566 : is_duplicate = .FALSE.
603 1376 : DO i = 1, jj + 1
604 1376 : IF (cdft_control%group(k)%atoms(j) == atomlist(i)) THEN
605 : is_duplicate = .TRUE.
606 : EXIT
607 : END IF
608 : END DO
609 870 : IF (.NOT. is_duplicate) THEN
610 528 : jj = jj + 1
611 528 : atomlist(jj) = cdft_control%group(k)%atoms(j)
612 : END IF
613 : END DO
614 : END DO
615 282 : CALL reallocate(atomlist, 1, jj)
616 : CALL section_vals_val_get(cdft_control_section, "ATOMIC_CHARGES", &
617 282 : l_val=cdft_control%atomic_charges)
618 : ! Parse any dummy atoms (no constraint, just charges)
619 282 : IF (cdft_control%atomic_charges) THEN
620 116 : group_section => section_vals_get_subs_vals(cdft_control_section, "DUMMY_ATOMS")
621 116 : CALL section_vals_get(group_section, explicit=exists)
622 116 : IF (exists) THEN
623 : ! First determine how many atoms there are
624 2 : natoms = 0
625 2 : CALL section_vals_val_get(group_section, "ATOMS", n_rep_val=n_rep)
626 4 : DO j = 1, n_rep
627 2 : CALL section_vals_val_get(group_section, "ATOMS", i_rep_val=j, i_vals=tmplist)
628 2 : IF (SIZE(tmplist) < 1) THEN
629 0 : CPABORT("DUMMY_ATOMS must contain at least 1 atom.")
630 : END IF
631 4 : natoms = natoms + SIZE(tmplist)
632 : END DO
633 6 : ALLOCATE (dummylist(natoms))
634 : ! Now parse
635 2 : jj = 0
636 4 : DO j = 1, n_rep
637 2 : CALL section_vals_val_get(group_section, "ATOMS", i_rep_val=j, i_vals=tmplist)
638 6 : DO i = 1, SIZE(tmplist)
639 2 : jj = jj + 1
640 4 : dummylist(jj) = tmplist(i)
641 : END DO
642 : END DO
643 : ! Check for duplicates
644 4 : DO j = 1, natoms
645 4 : DO i = j + 1, natoms
646 2 : IF (dummylist(i) == dummylist(j)) THEN
647 0 : CPABORT("Duplicate atoms defined in section DUMMY_ATOMS.")
648 : END IF
649 : END DO
650 : END DO
651 : ! Check that a dummy atom is not included in any ATOM_GROUP
652 6 : DO j = 1, SIZE(atomlist)
653 6 : DO i = 1, SIZE(dummylist)
654 4 : IF (dummylist(i) == atomlist(j)) THEN
655 : CALL cp_abort(__LOCATION__, &
656 0 : "Duplicate atoms defined in sections ATOM_GROUP and DUMMY_ATOMS.")
657 : END IF
658 : END DO
659 : END DO
660 : END IF
661 : END IF
662 : ! Join dummy atoms and constraint atoms into one list
663 282 : IF (ASSOCIATED(dummylist)) THEN
664 2 : cdft_control%natoms = SIZE(atomlist) + SIZE(dummylist)
665 : ELSE
666 280 : cdft_control%natoms = SIZE(atomlist)
667 : END IF
668 846 : ALLOCATE (cdft_control%atoms(cdft_control%natoms))
669 564 : ALLOCATE (cdft_control%is_constraint(cdft_control%natoms))
670 750 : IF (cdft_control%atomic_charges) ALLOCATE (cdft_control%charge(cdft_control%natoms))
671 810 : cdft_control%atoms(1:SIZE(atomlist)) = atomlist
672 282 : IF (ASSOCIATED(dummylist)) THEN
673 4 : cdft_control%atoms(1 + SIZE(atomlist):) = dummylist
674 2 : DEALLOCATE (dummylist)
675 : END IF
676 812 : cdft_control%is_constraint = .FALSE.
677 810 : cdft_control%is_constraint(1:SIZE(atomlist)) = .TRUE.
678 282 : DEALLOCATE (atomlist)
679 : ! Get constraint potential definitions from input
680 846 : ALLOCATE (cdft_control%strength(nvar))
681 564 : ALLOCATE (cdft_control%value(nvar))
682 564 : ALLOCATE (cdft_control%target(nvar))
683 282 : CALL section_vals_val_get(cdft_control_section, "STRENGTH", r_vals=rtmplist)
684 282 : IF (SIZE(rtmplist) /= nvar) THEN
685 : CALL cp_abort(__LOCATION__, &
686 : "The length of keyword STRENGTH is incorrect. "// &
687 : "Expected "//TRIM(ADJUSTL(cp_to_string(nvar)))// &
688 : " value(s), got "// &
689 0 : TRIM(ADJUSTL(cp_to_string(SIZE(rtmplist))))//" value(s).")
690 : END IF
691 586 : DO j = 1, nvar
692 586 : cdft_control%strength(j) = rtmplist(j)
693 : END DO
694 282 : CALL section_vals_val_get(cdft_control_section, "TARGET", r_vals=rtmplist)
695 282 : IF (SIZE(rtmplist) /= nvar) THEN
696 : CALL cp_abort(__LOCATION__, &
697 : "The length of keyword TARGET is incorrect. "// &
698 : "Expected "//TRIM(ADJUSTL(cp_to_string(nvar)))// &
699 : " value(s), got "// &
700 0 : TRIM(ADJUSTL(cp_to_string(SIZE(rtmplist))))//" value(s).")
701 : END IF
702 586 : DO j = 1, nvar
703 586 : cdft_control%target(j) = rtmplist(j)
704 : END DO
705 : ! Read fragment constraint definitions
706 282 : IF (cdft_control%fragment_density) THEN
707 : CALL section_vals_val_get(cdft_control_section, "FRAGMENT_A_FILE_NAME", &
708 10 : c_val=cdft_control%fragment_a_fname)
709 : CALL section_vals_val_get(cdft_control_section, "FRAGMENT_B_FILE_NAME", &
710 10 : c_val=cdft_control%fragment_b_fname)
711 : CALL section_vals_val_get(cdft_control_section, "FRAGMENT_A_SPIN_FILE", &
712 10 : c_val=cdft_control%fragment_a_spin_fname)
713 : CALL section_vals_val_get(cdft_control_section, "FRAGMENT_B_SPIN_FILE", &
714 10 : c_val=cdft_control%fragment_b_spin_fname)
715 : CALL section_vals_val_get(cdft_control_section, "FLIP_FRAGMENT_A", &
716 10 : l_val=cdft_control%flip_fragment(1))
717 : CALL section_vals_val_get(cdft_control_section, "FLIP_FRAGMENT_B", &
718 10 : l_val=cdft_control%flip_fragment(2))
719 : END IF
720 :
721 282 : END SUBROUTINE read_constraint_definitions
722 :
723 : ! **************************************************************************************************
724 : !> \brief reads the input parameters needed for CDFT with OT
725 : !> \param qs_control the qs_control which holds the CDFT control type
726 : !> \param cdft_control_section the input section for CDFT
727 : !> \author Nico Holmberg [12.2015]
728 : ! **************************************************************************************************
729 564 : SUBROUTINE read_cdft_control_section(qs_control, cdft_control_section)
730 : TYPE(qs_control_type), INTENT(INOUT) :: qs_control
731 : TYPE(section_vals_type), POINTER :: cdft_control_section
732 :
733 : INTEGER :: k, nvar
734 : LOGICAL :: exists
735 : TYPE(cdft_control_type), POINTER :: cdft_control
736 : TYPE(section_vals_type), POINTER :: becke_constraint_section, group_section, &
737 : hirshfeld_constraint_section, &
738 : outer_scf_section, print_section
739 :
740 282 : NULLIFY (outer_scf_section, hirshfeld_constraint_section, becke_constraint_section, &
741 282 : print_section, group_section)
742 282 : cdft_control => qs_control%cdft_control
743 282 : CPASSERT(ASSOCIATED(cdft_control))
744 282 : group_section => section_vals_get_subs_vals(cdft_control_section, "ATOM_GROUP")
745 282 : CALL section_vals_get(group_section, n_repetition=nvar, explicit=exists)
746 :
747 : CALL section_vals_val_get(cdft_control_section, "TYPE_OF_CONSTRAINT", &
748 282 : i_val=qs_control%cdft_control%type)
749 :
750 282 : IF (cdft_control%type /= outer_scf_none) THEN
751 : CALL section_vals_val_get(cdft_control_section, "REUSE_PRECOND", &
752 282 : l_val=cdft_control%reuse_precond)
753 : CALL section_vals_val_get(cdft_control_section, "PRECOND_FREQ", &
754 282 : i_val=cdft_control%precond_freq)
755 : CALL section_vals_val_get(cdft_control_section, "MAX_REUSE", &
756 282 : i_val=cdft_control%max_reuse)
757 : CALL section_vals_val_get(cdft_control_section, "PURGE_HISTORY", &
758 282 : l_val=cdft_control%purge_history)
759 : CALL section_vals_val_get(cdft_control_section, "PURGE_FREQ", &
760 282 : i_val=cdft_control%purge_freq)
761 : CALL section_vals_val_get(cdft_control_section, "PURGE_OFFSET", &
762 282 : i_val=cdft_control%purge_offset)
763 : CALL section_vals_val_get(cdft_control_section, "COUNTER", &
764 282 : i_val=cdft_control%ienergy)
765 282 : print_section => section_vals_get_subs_vals(cdft_control_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION")
766 282 : CALL section_vals_get(print_section, explicit=cdft_control%print_weight)
767 :
768 282 : outer_scf_section => section_vals_get_subs_vals(cdft_control_section, "OUTER_SCF")
769 282 : CALL outer_scf_read_parameters(cdft_control%constraint_control, outer_scf_section)
770 282 : IF (cdft_control%constraint_control%have_scf) THEN
771 282 : IF (cdft_control%constraint_control%type /= outer_scf_cdft_constraint) THEN
772 0 : CPABORT("Unsupported CDFT constraint.")
773 : END IF
774 : ! Constraint definitions
775 282 : CALL read_constraint_definitions(cdft_control, cdft_control_section)
776 : ! Constraint-specific initializations
777 536 : SELECT CASE (cdft_control%type)
778 : CASE (outer_scf_becke_constraint)
779 254 : becke_constraint_section => section_vals_get_subs_vals(cdft_control_section, "BECKE_CONSTRAINT")
780 254 : CALL section_vals_get(becke_constraint_section, explicit=exists)
781 254 : IF (.NOT. exists) CPABORT("BECKE_CONSTRAINT section is missing.")
782 530 : DO k = 1, nvar
783 530 : NULLIFY (cdft_control%group(k)%gradients)
784 : END DO
785 254 : CALL read_becke_section(cdft_control, becke_constraint_section)
786 : CASE (outer_scf_hirshfeld_constraint)
787 28 : hirshfeld_constraint_section => section_vals_get_subs_vals(cdft_control_section, "HIRSHFELD_CONSTRAINT")
788 28 : CALL section_vals_get(hirshfeld_constraint_section, explicit=exists)
789 28 : IF (.NOT. exists) CPABORT("HIRSHFELD_CONSTRAINT section is missing.")
790 56 : DO k = 1, nvar
791 28 : NULLIFY (cdft_control%group(k)%gradients_x)
792 28 : NULLIFY (cdft_control%group(k)%gradients_y)
793 56 : NULLIFY (cdft_control%group(k)%gradients_z)
794 : END DO
795 28 : CALL read_hirshfeld_constraint_section(cdft_control, hirshfeld_constraint_section)
796 : CASE DEFAULT
797 564 : CPABORT("Unknown constraint type.")
798 : END SELECT
799 :
800 282 : CALL cite_reference(Holmberg2017)
801 282 : CALL cite_reference(Holmberg2018)
802 : ELSE
803 0 : qs_control%cdft = .FALSE.
804 : END IF
805 : ELSE
806 0 : qs_control%cdft = .FALSE.
807 : END IF
808 :
809 282 : END SUBROUTINE read_cdft_control_section
810 :
811 : ! **************************************************************************************************
812 : !> \brief reads the input parameters needed for Hirshfeld constraint
813 : !> \param cdft_control the cdft_control which holds the Hirshfeld constraint
814 : !> \param hirshfeld_section the input section for a Hirshfeld constraint
815 : ! **************************************************************************************************
816 28 : SUBROUTINE read_hirshfeld_constraint_section(cdft_control, hirshfeld_section)
817 : TYPE(cdft_control_type), INTENT(INOUT) :: cdft_control
818 : TYPE(section_vals_type), POINTER :: hirshfeld_section
819 :
820 : LOGICAL :: exists
821 28 : REAL(KIND=dp), DIMENSION(:), POINTER :: rtmplist
822 : TYPE(hirshfeld_constraint_type), POINTER :: hirshfeld_control
823 :
824 28 : NULLIFY (rtmplist)
825 28 : hirshfeld_control => cdft_control%hirshfeld_control
826 0 : CPASSERT(ASSOCIATED(hirshfeld_control))
827 :
828 28 : CALL section_vals_val_get(hirshfeld_section, "SHAPE_FUNCTION", i_val=hirshfeld_control%shape_function)
829 28 : CALL section_vals_val_get(hirshfeld_section, "GAUSSIAN_SHAPE", i_val=hirshfeld_control%gaussian_shape)
830 28 : CALL section_vals_val_get(hirshfeld_section, "GAUSSIAN_RADIUS", r_val=hirshfeld_control%radius)
831 28 : CALL section_vals_val_get(hirshfeld_section, "USE_BOHR", l_val=hirshfeld_control%use_bohr)
832 28 : CALL section_vals_val_get(hirshfeld_section, "USE_ATOMIC_CUTOFF", l_val=hirshfeld_control%use_atomic_cutoff)
833 28 : CALL section_vals_val_get(hirshfeld_section, "PRINT_DENSITY", l_val=hirshfeld_control%print_density)
834 28 : CALL section_vals_val_get(hirshfeld_section, "EPS_CUTOFF", r_val=hirshfeld_control%eps_cutoff)
835 28 : CALL section_vals_val_get(hirshfeld_section, "ATOMIC_CUTOFF", r_val=hirshfeld_control%atomic_cutoff)
836 :
837 28 : IF (.NOT. hirshfeld_control%use_bohr) THEN
838 28 : hirshfeld_control%radius = cp_unit_from_cp2k(hirshfeld_control%radius, "angstrom")
839 : END IF
840 :
841 28 : IF (hirshfeld_control%shape_function == shape_function_gaussian .AND. &
842 : hirshfeld_control%gaussian_shape == radius_user) THEN
843 0 : CALL section_vals_val_get(hirshfeld_section, "ATOMIC_RADII", explicit=exists)
844 0 : IF (.NOT. exists) CPABORT("Keyword ATOMIC_RADII is missing.")
845 0 : CALL section_vals_val_get(hirshfeld_section, "ATOMIC_RADII", r_vals=rtmplist)
846 0 : CPASSERT(SIZE(rtmplist) > 0)
847 0 : ALLOCATE (hirshfeld_control%radii(SIZE(rtmplist)))
848 0 : hirshfeld_control%radii(:) = rtmplist
849 : END IF
850 :
851 28 : CALL create_hirshfeld_type(hirshfeld_control%hirshfeld_env)
852 : CALL set_hirshfeld_info(hirshfeld_control%hirshfeld_env, &
853 : shape_function_type=hirshfeld_control%shape_function, &
854 : iterative=.FALSE., &
855 : radius_type=hirshfeld_control%gaussian_shape, &
856 28 : use_bohr=hirshfeld_control%use_bohr)
857 :
858 28 : END SUBROUTINE read_hirshfeld_constraint_section
859 :
860 : ! **************************************************************************************************
861 : !> \brief Calculate fout = fun1/fun2 or fout = fun1*fun2
862 : !> \param fout the output 3D potential
863 : !> \param fun1 the first input 3D potential
864 : !> \param fun2 the second input 3D potential
865 : !> \param divide logical that decides whether to divide or multiply the input potentials
866 : !> \param small customisable parameter to determine lower bound of division
867 : ! **************************************************************************************************
868 40 : SUBROUTINE hfun_scale(fout, fun1, fun2, divide, small)
869 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: fout
870 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: fun1, fun2
871 : LOGICAL, INTENT(IN) :: divide
872 : REAL(KIND=dp), INTENT(IN) :: small
873 :
874 : INTEGER :: i1, i2, i3, n1, n2, n3
875 :
876 40 : n1 = SIZE(fout, 1)
877 40 : n2 = SIZE(fout, 2)
878 40 : n3 = SIZE(fout, 3)
879 40 : CPASSERT(n1 == SIZE(fun1, 1))
880 40 : CPASSERT(n2 == SIZE(fun1, 2))
881 40 : CPASSERT(n3 == SIZE(fun1, 3))
882 40 : CPASSERT(n1 == SIZE(fun2, 1))
883 40 : CPASSERT(n2 == SIZE(fun2, 2))
884 40 : CPASSERT(n3 == SIZE(fun2, 3))
885 :
886 40 : IF (divide) THEN
887 1640 : DO i3 = 1, n3
888 65640 : DO i2 = 1, n2
889 1409600 : DO i1 = 1, n1
890 1408000 : IF (fun2(i1, i2, i3) > small) THEN
891 1163532 : fout(i1, i2, i3) = fun1(i1, i2, i3)/fun2(i1, i2, i3)
892 : ELSE
893 180468 : fout(i1, i2, i3) = 0.0_dp
894 : END IF
895 : END DO
896 : END DO
897 : END DO
898 : ELSE
899 0 : DO i3 = 1, n3
900 0 : DO i2 = 1, n2
901 0 : DO i1 = 1, n1
902 0 : fout(i1, i2, i3) = fun1(i1, i2, i3)*fun2(i1, i2, i3)
903 : END DO
904 : END DO
905 : END DO
906 : END IF
907 :
908 40 : END SUBROUTINE hfun_scale
909 :
910 : ! **************************************************************************************************
911 : !> \brief Determine confinement bounds along confinement dir (hardcoded to be z)
912 : !> and optionally zero entries below a given threshold
913 : !> \param fun input 3D potential (real space)
914 : !> \param th threshold for screening values
915 : !> \param just_bounds if the bounds should be computed without zeroing values
916 : !> \param bounds the confinement bounds: fun is nonzero only between these values along 3rd dimension
917 : ! **************************************************************************************************
918 66 : SUBROUTINE hfun_zero(fun, th, just_bounds, bounds)
919 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: fun
920 : REAL(KIND=dp), INTENT(IN) :: th
921 : LOGICAL :: just_bounds
922 : INTEGER, OPTIONAL :: bounds(2)
923 :
924 : INTEGER :: i1, i2, i3, lb, n1, n2, n3, nzeroed, &
925 : nzeroed_inner, ub
926 : LOGICAL :: lb_final, ub_final
927 :
928 66 : n1 = SIZE(fun, 1)
929 66 : n2 = SIZE(fun, 2)
930 66 : n3 = SIZE(fun, 3)
931 66 : IF (just_bounds) THEN
932 64 : CPASSERT(PRESENT(bounds))
933 : lb = 1
934 : lb_final = .FALSE.
935 : ub_final = .FALSE.
936 : END IF
937 :
938 2898 : DO i3 = 1, n3
939 2832 : IF (just_bounds) nzeroed = 0
940 20466 : DO i2 = 1, n2
941 20306 : IF (just_bounds) nzeroed_inner = 0
942 518013 : DO i1 = 1, n1
943 520685 : IF (fun(i1, i2, i3) < th) THEN
944 446635 : IF (just_bounds) THEN
945 433707 : nzeroed_inner = nzeroed_inner + 1
946 : ELSE
947 12928 : fun(i1, i2, i3) = 0.0_dp
948 : END IF
949 : ELSE
950 53744 : IF (just_bounds) EXIT
951 : END IF
952 : END DO
953 23138 : IF (just_bounds) THEN
954 17106 : IF (nzeroed_inner < n1) EXIT
955 14434 : nzeroed = nzeroed + nzeroed_inner
956 : END IF
957 : END DO
958 2898 : IF (just_bounds) THEN
959 2752 : IF (nzeroed == (n2*n1)) THEN
960 80 : IF (.NOT. lb_final) THEN
961 : lb = i3
962 56 : ELSE IF (.NOT. ub_final) THEN
963 8 : ub = i3
964 8 : ub_final = .TRUE.
965 : END IF
966 : ELSE
967 : IF (.NOT. lb_final) lb_final = .TRUE.
968 : IF (ub_final) ub_final = .FALSE. ! Safeguard against "holes"
969 : END IF
970 : END IF
971 : END DO
972 66 : IF (just_bounds) THEN
973 64 : IF (.NOT. ub_final) ub = n3
974 64 : bounds(1) = lb
975 64 : bounds(2) = ub
976 192 : bounds = bounds - (n3/2) - 1
977 : END IF
978 :
979 66 : END SUBROUTINE hfun_zero
980 :
981 : ! **************************************************************************************************
982 : !> \brief Initializes Gaussian Hirshfeld constraints
983 : !> \param qs_env the qs_env where to build the constraint
984 : !> \author Nico Holmberg (09.2018)
985 : ! **************************************************************************************************
986 22 : SUBROUTINE hirshfeld_constraint_init(qs_env)
987 : TYPE(qs_environment_type), POINTER :: qs_env
988 :
989 : CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_constraint_init'
990 :
991 : CHARACTER(len=2) :: element_symbol
992 : INTEGER :: handle, iat, iatom, igroup, ikind, ip, &
993 : iw, natom, nkind
994 22 : INTEGER, DIMENSION(:), POINTER :: atom_list
995 : REAL(KIND=dp) :: zeff
996 22 : REAL(KIND=dp), DIMENSION(:), POINTER :: radii_list
997 22 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
998 : TYPE(atomic_kind_type), POINTER :: atomic_kind
999 : TYPE(cdft_control_type), POINTER :: cdft_control
1000 22 : TYPE(cdft_group_type), DIMENSION(:), POINTER :: group
1001 : TYPE(cp_logger_type), POINTER :: logger
1002 : TYPE(dft_control_type), POINTER :: dft_control
1003 : TYPE(hirshfeld_constraint_type), POINTER :: hirshfeld_control
1004 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
1005 22 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1006 22 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1007 : TYPE(section_vals_type), POINTER :: print_section
1008 :
1009 22 : NULLIFY (cdft_control, hirshfeld_control, hirshfeld_env, qs_kind_set, atomic_kind_set, &
1010 22 : radii_list, dft_control, group, atomic_kind, atom_list)
1011 22 : CALL timeset(routineN, handle)
1012 :
1013 22 : logger => cp_get_default_logger()
1014 22 : print_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
1015 22 : iw = cp_print_key_unit_nr(logger, print_section, "PROGRAM_RUN_INFO", extension=".cdftLog")
1016 :
1017 22 : CALL get_qs_env(qs_env, dft_control=dft_control)
1018 22 : cdft_control => dft_control%qs_control%cdft_control
1019 22 : hirshfeld_control => cdft_control%hirshfeld_control
1020 22 : hirshfeld_env => hirshfeld_control%hirshfeld_env
1021 :
1022 : ! Setup the Hirshfeld shape function
1023 22 : IF (.NOT. ASSOCIATED(hirshfeld_env%kind_shape_fn)) THEN
1024 : hirshfeld_env => hirshfeld_control%hirshfeld_env
1025 22 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
1026 22 : CPASSERT(ASSOCIATED(qs_kind_set))
1027 22 : nkind = SIZE(qs_kind_set)
1028 : ! Parse atomic radii for setting up Gaussian shape function
1029 22 : IF (ASSOCIATED(hirshfeld_control%radii)) THEN
1030 0 : IF (.NOT. SIZE(atomic_kind_set) == SIZE(hirshfeld_control%radii)) THEN
1031 : CALL cp_abort(__LOCATION__, &
1032 : "Length of keyword HIRSHFELD_CONSTRAINT\ATOMIC_RADII does not "// &
1033 0 : "match number of atomic kinds in the input coordinate file.")
1034 : END IF
1035 :
1036 0 : ALLOCATE (radii_list(SIZE(hirshfeld_control%radii)))
1037 0 : DO ikind = 1, SIZE(hirshfeld_control%radii)
1038 0 : IF (hirshfeld_control%use_bohr) THEN
1039 0 : radii_list(ikind) = hirshfeld_control%radii(ikind)
1040 : ELSE
1041 0 : radii_list(ikind) = cp_unit_from_cp2k(hirshfeld_control%radii(ikind), "angstrom")
1042 : END IF
1043 : END DO
1044 : END IF
1045 : ! radius/radii_list parameters are optional for shape_function_density
1046 : CALL create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, &
1047 : radius=hirshfeld_control%radius, &
1048 22 : radii_list=radii_list)
1049 22 : IF (ASSOCIATED(radii_list)) DEALLOCATE (radii_list)
1050 : END IF
1051 :
1052 : ! Atomic reference charges (Mulliken not supported)
1053 22 : IF (.NOT. ASSOCIATED(hirshfeld_env%charges)) THEN
1054 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
1055 22 : nkind=nkind, natom=natom)
1056 66 : ALLOCATE (hirshfeld_env%charges(natom))
1057 66 : DO ikind = 1, nkind
1058 44 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
1059 44 : atomic_kind => atomic_kind_set(ikind)
1060 44 : CALL get_atomic_kind(atomic_kind, atom_list=atom_list)
1061 154 : DO iat = 1, SIZE(atom_list)
1062 44 : iatom = atom_list(iat)
1063 88 : hirshfeld_env%charges(iatom) = zeff
1064 : END DO
1065 : END DO
1066 : END IF
1067 :
1068 : ! Print some additional information about the calculation on the first iteration
1069 22 : IF (cdft_control%first_iteration) THEN
1070 22 : IF (iw > 0) THEN
1071 12 : group => cdft_control%group
1072 12 : CALL get_qs_env(qs_env, particle_set=particle_set)
1073 12 : IF (ASSOCIATED(hirshfeld_control%radii)) THEN
1074 : WRITE (iw, '(T3,A)') &
1075 0 : 'Atom Element Gaussian radius (angstrom)'
1076 0 : DO iatom = 1, natom
1077 0 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol)
1078 : WRITE (iw, "(i7,T15,A2,T37,F8.3)") &
1079 0 : iatom, ADJUSTR(element_symbol), cp_unit_from_cp2k(hirshfeld_control%radii(iatom), "angstrom")
1080 : END DO
1081 : WRITE (iw, '(T3,A)') &
1082 0 : '------------------------------------------------------------------------'
1083 : END IF
1084 : WRITE (iw, '(/,T3,A,T60)') &
1085 12 : '----------------------- CDFT group definitions -------------------------'
1086 24 : DO igroup = 1, SIZE(group)
1087 12 : IF (igroup > 1) WRITE (iw, '(T3,A)') ' '
1088 : WRITE (iw, '(T5,A,I5,A,I5)') &
1089 12 : 'Atomic group', igroup, ' of ', SIZE(group)
1090 12 : WRITE (iw, '(T5,A)') 'Atom Element Coefficient'
1091 47 : DO ip = 1, SIZE(group(igroup)%atoms)
1092 23 : iatom = group(igroup)%atoms(ip)
1093 23 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, element_symbol=element_symbol)
1094 35 : WRITE (iw, '(i8,T16,A2,T23,F8.3)') iatom, ADJUSTR(element_symbol), group(igroup)%coeff(ip)
1095 : END DO
1096 : END DO
1097 : WRITE (iw, '(T3,A)') &
1098 12 : '------------------------------------------------------------------------'
1099 : END IF
1100 22 : cdft_control%first_iteration = .FALSE.
1101 : END IF
1102 :
1103 : ! Radii no longer needed
1104 22 : IF (ASSOCIATED(hirshfeld_control%radii)) DEALLOCATE (hirshfeld_control%radii)
1105 22 : CALL timestop(handle)
1106 :
1107 22 : END SUBROUTINE hirshfeld_constraint_init
1108 :
1109 : ! **************************************************************************************************
1110 : !> \brief Prints information about CDFT constraints
1111 : !> \param qs_env the qs_env where to build the constraint
1112 : !> \param electronic_charge the CDFT charges
1113 : !> \par History
1114 : !> Created 9.2018 [Nico Holmberg]
1115 : ! **************************************************************************************************
1116 3410 : SUBROUTINE cdft_constraint_print(qs_env, electronic_charge)
1117 : TYPE(qs_environment_type), POINTER :: qs_env
1118 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: electronic_charge
1119 :
1120 : CHARACTER(len=2) :: element_symbol
1121 : INTEGER :: iatom, ikind, iw, jatom
1122 : REAL(kind=dp) :: tc(2), zeff
1123 : TYPE(cdft_control_type), POINTER :: cdft_control
1124 : TYPE(cp_logger_type), POINTER :: logger
1125 : TYPE(dft_control_type), POINTER :: dft_control
1126 3410 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1127 3410 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1128 : TYPE(section_vals_type), POINTER :: cdft_constraint_section
1129 :
1130 3410 : NULLIFY (cdft_constraint_section, logger, particle_set, dft_control, qs_kind_set)
1131 6820 : logger => cp_get_default_logger()
1132 :
1133 : CALL get_qs_env(qs_env, &
1134 : particle_set=particle_set, &
1135 : dft_control=dft_control, &
1136 3410 : qs_kind_set=qs_kind_set)
1137 3410 : CPASSERT(ASSOCIATED(qs_kind_set))
1138 :
1139 3410 : cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
1140 3410 : iw = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", extension=".cdftLog")
1141 3410 : cdft_control => dft_control%qs_control%cdft_control
1142 :
1143 : ! Print constraint information
1144 3410 : CALL qs_scf_cdft_constraint_info(iw, cdft_control)
1145 :
1146 : ! Print weight function(s) to cube file(s) whenever weight is (re)built
1147 3410 : IF (cdft_control%print_weight .AND. cdft_control%need_pot) THEN
1148 2 : CALL cdft_print_weight_function(qs_env)
1149 : END IF
1150 :
1151 : ! Print atomic CDFT charges
1152 3410 : IF (iw > 0 .AND. cdft_control%atomic_charges) THEN
1153 799 : IF (.NOT. cdft_control%fragment_density) THEN
1154 794 : IF (dft_control%nspins == 1) THEN
1155 : WRITE (iw, '(/,T3,A)') &
1156 0 : '-------------------------------- CDFT atomic charges --------------------------------'
1157 : WRITE (iw, '(T3,A,A)') &
1158 0 : '#Atom Element Is_constraint', ' Core charge Population (total)'// &
1159 0 : ' Net charge'
1160 0 : tc = 0.0_dp
1161 0 : DO iatom = 1, cdft_control%natoms
1162 0 : jatom = cdft_control%atoms(iatom)
1163 : CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, &
1164 : element_symbol=element_symbol, &
1165 0 : kind_number=ikind)
1166 0 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
1167 : WRITE (iw, "(i7,T15,A2,T23,L10,T39,F8.3,T61,F8.3,T81,F8.3)") &
1168 0 : jatom, ADJUSTR(element_symbol), cdft_control%is_constraint(iatom), zeff, electronic_charge(iatom, 1), &
1169 0 : (zeff - electronic_charge(iatom, 1))
1170 0 : tc(1) = tc(1) + (zeff - electronic_charge(iatom, 1))
1171 : END DO
1172 0 : WRITE (iw, '(/,T3,A,T81,F8.3,/)') "Total Charge: ", tc(1)
1173 : ELSE
1174 : WRITE (iw, '(/,T3,A)') &
1175 794 : '------------------------------------------ CDFT atomic charges -------------------------------------------'
1176 : WRITE (iw, '(T3,A,A)') &
1177 794 : '#Atom Element Is_constraint', ' Core charge Population (alpha, beta)'// &
1178 1588 : ' Net charge Spin population'
1179 794 : tc = 0.0_dp
1180 2343 : DO iatom = 1, cdft_control%natoms
1181 1549 : jatom = cdft_control%atoms(iatom)
1182 : CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, &
1183 : element_symbol=element_symbol, &
1184 1549 : kind_number=ikind)
1185 1549 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
1186 : WRITE (iw, "(i7,T15,A2,T23,L10,T39,F8.3,T53,F8.3,T67,F8.3,T81,F8.3,T102,F8.3)") &
1187 1549 : jatom, ADJUSTR(element_symbol), &
1188 1549 : cdft_control%is_constraint(iatom), &
1189 1549 : zeff, electronic_charge(iatom, 1), electronic_charge(iatom, 2), &
1190 1549 : (zeff - electronic_charge(iatom, 1) - electronic_charge(iatom, 2)), &
1191 3098 : electronic_charge(iatom, 1) - electronic_charge(iatom, 2)
1192 1549 : tc(1) = tc(1) + (zeff - electronic_charge(iatom, 1) - electronic_charge(iatom, 2))
1193 3892 : tc(2) = tc(2) + (electronic_charge(iatom, 1) - electronic_charge(iatom, 2))
1194 : END DO
1195 794 : WRITE (iw, '(/,T3,A,T81,F8.3,T102,F8.3/)') "Total Charge and Spin Moment: ", tc(1), tc(2)
1196 : END IF
1197 : ELSE
1198 8 : IF (ALL(cdft_control%group(:)%constraint_type == cdft_charge_constraint)) THEN
1199 : WRITE (iw, '(/,T3,A)') &
1200 3 : '-------------------------------- CDFT atomic charges --------------------------------'
1201 3 : IF (dft_control%nspins == 1) THEN
1202 : WRITE (iw, '(T3,A,A)') &
1203 0 : '#Atom Element Is_constraint', ' Fragment charge Population (total)'// &
1204 0 : ' Net charge'
1205 : ELSE
1206 : WRITE (iw, '(T3,A,A)') &
1207 3 : '#Atom Element Is_constraint', ' Fragment charge Population (alpha, beta)'// &
1208 6 : ' Net charge'
1209 : END IF
1210 3 : tc = 0.0_dp
1211 7 : DO iatom = 1, cdft_control%natoms
1212 4 : jatom = cdft_control%atoms(iatom)
1213 : CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, &
1214 : element_symbol=element_symbol, &
1215 4 : kind_number=ikind)
1216 7 : IF (dft_control%nspins == 1) THEN
1217 : WRITE (iw, "(i7,T15,A2,T23,L10,T43,F8.3,T65,F8.3,T81,F8.3)") &
1218 0 : jatom, ADJUSTR(element_symbol), &
1219 0 : cdft_control%is_constraint(iatom), &
1220 0 : cdft_control%charges_fragment(iatom, 1), &
1221 0 : electronic_charge(iatom, 1), &
1222 : (electronic_charge(iatom, 1) - &
1223 0 : cdft_control%charges_fragment(iatom, 1))
1224 : tc(1) = tc(1) + (electronic_charge(iatom, 1) - &
1225 0 : cdft_control%charges_fragment(iatom, 1))
1226 : ELSE
1227 : WRITE (iw, "(i7,T15,A2,T23,L10,T43,F8.3,T57,F8.3,T69,F8.3,T81,F8.3)") &
1228 4 : jatom, ADJUSTR(element_symbol), &
1229 4 : cdft_control%is_constraint(iatom), &
1230 4 : cdft_control%charges_fragment(iatom, 1), &
1231 4 : electronic_charge(iatom, 1), electronic_charge(iatom, 2), &
1232 : (electronic_charge(iatom, 1) + electronic_charge(iatom, 2) - &
1233 8 : cdft_control%charges_fragment(iatom, 1))
1234 : tc(1) = tc(1) + (electronic_charge(iatom, 1) + electronic_charge(iatom, 2) - &
1235 4 : cdft_control%charges_fragment(iatom, 1))
1236 : END IF
1237 : END DO
1238 3 : WRITE (iw, '(/,T3,A,T81,F8.3,/)') "Total Charge: ", tc(1)
1239 : ELSE
1240 : WRITE (iw, '(/,T3,A)') &
1241 2 : '------------------------------------------ CDFT atomic charges -------------------------------------------'
1242 : WRITE (iw, '(T3,A,A)') &
1243 2 : '#Atom Element Is_constraint', ' Fragment charge/spin moment'// &
1244 4 : ' Population (alpha, beta) Net charge/spin moment'
1245 2 : tc = 0.0_dp
1246 5 : DO iatom = 1, cdft_control%natoms
1247 3 : jatom = cdft_control%atoms(iatom)
1248 : CALL get_atomic_kind(atomic_kind=particle_set(jatom)%atomic_kind, &
1249 : element_symbol=element_symbol, &
1250 3 : kind_number=ikind)
1251 : WRITE (iw, "(i7,T15,A2,T22,L10,T40,F8.3,T52,F8.3,T66,F8.3,T78,F8.3,T90,F8.3,T102,F8.3)") &
1252 3 : jatom, ADJUSTR(element_symbol), &
1253 3 : cdft_control%is_constraint(iatom), &
1254 3 : cdft_control%charges_fragment(iatom, 1), &
1255 3 : cdft_control%charges_fragment(iatom, 2), &
1256 3 : electronic_charge(iatom, 1), electronic_charge(iatom, 2), &
1257 : (electronic_charge(iatom, 1) + electronic_charge(iatom, 2) - &
1258 3 : cdft_control%charges_fragment(iatom, 1)), &
1259 : (electronic_charge(iatom, 1) - electronic_charge(iatom, 2) - &
1260 6 : cdft_control%charges_fragment(iatom, 2))
1261 : tc(1) = tc(1) + (electronic_charge(iatom, 1) + electronic_charge(iatom, 2) - &
1262 3 : cdft_control%charges_fragment(iatom, 1))
1263 : tc(2) = tc(2) + (electronic_charge(iatom, 1) - electronic_charge(iatom, 2) - &
1264 8 : cdft_control%charges_fragment(iatom, 2))
1265 : END DO
1266 2 : WRITE (iw, '(/,T3,A,T90,F8.3,T102,F8.3/)') "Total Charge and Spin Moment: ", tc(1), tc(2)
1267 : END IF
1268 : END IF
1269 : END IF
1270 :
1271 3410 : END SUBROUTINE cdft_constraint_print
1272 :
1273 : ! **************************************************************************************************
1274 : !> \brief Prints CDFT weight functions to cube files
1275 : !> \param qs_env ...
1276 : ! **************************************************************************************************
1277 2 : SUBROUTINE cdft_print_weight_function(qs_env)
1278 : TYPE(qs_environment_type), POINTER :: qs_env
1279 :
1280 : CHARACTER(LEN=default_path_length) :: middle_name
1281 : INTEGER :: igroup, unit_nr
1282 : LOGICAL :: mpi_io
1283 : TYPE(cdft_control_type), POINTER :: cdft_control
1284 : TYPE(cp_logger_type), POINTER :: logger
1285 : TYPE(dft_control_type), POINTER :: dft_control
1286 : TYPE(mp_para_env_type), POINTER :: para_env
1287 : TYPE(particle_list_type), POINTER :: particles
1288 : TYPE(qs_subsys_type), POINTER :: subsys
1289 : TYPE(section_vals_type), POINTER :: cdft_constraint_section
1290 :
1291 2 : NULLIFY (cdft_constraint_section, logger, particles, dft_control, &
1292 2 : para_env, subsys, cdft_control)
1293 2 : logger => cp_get_default_logger()
1294 :
1295 2 : CALL get_qs_env(qs_env, subsys=subsys, para_env=para_env, dft_control=dft_control)
1296 2 : CALL qs_subsys_get(subsys, particles=particles)
1297 2 : cdft_control => dft_control%qs_control%cdft_control
1298 2 : cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
1299 :
1300 4 : DO igroup = 1, SIZE(cdft_control%group)
1301 2 : mpi_io = .TRUE.
1302 2 : middle_name = "cdft_weight_"//TRIM(ADJUSTL(cp_to_string(igroup)))
1303 : unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", &
1304 : middle_name=middle_name, &
1305 : extension=".cube", file_position="REWIND", &
1306 2 : log_filename=.FALSE., mpi_io=mpi_io)
1307 : ! Note PROGRAM_RUN_INFO section neeeds to be active!
1308 2 : IF (para_env%is_source() .AND. unit_nr < 1) THEN
1309 : CALL cp_abort(__LOCATION__, &
1310 0 : "Please turn on PROGRAM_RUN_INFO to print CDFT weight function.")
1311 : END IF
1312 :
1313 : CALL cp_pw_to_cube(cdft_control%group(igroup)%weight, &
1314 : unit_nr, &
1315 : "CDFT Weight Function", &
1316 : particles=particles, &
1317 : stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE"), &
1318 2 : mpi_io=mpi_io)
1319 4 : CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io)
1320 : END DO
1321 :
1322 2 : END SUBROUTINE cdft_print_weight_function
1323 :
1324 : ! **************************************************************************************************
1325 : !> \brief Prints Hirshfeld weight function and promolecule density
1326 : !> \param qs_env ...
1327 : ! **************************************************************************************************
1328 0 : SUBROUTINE cdft_print_hirshfeld_density(qs_env)
1329 : TYPE(qs_environment_type), POINTER :: qs_env
1330 :
1331 : CHARACTER(LEN=default_path_length) :: middle_name
1332 : INTEGER :: iatom, igroup, unit_nr
1333 : LOGICAL :: mpi_io
1334 : TYPE(cdft_control_type), POINTER :: cdft_control
1335 : TYPE(cp_logger_type), POINTER :: logger
1336 : TYPE(dft_control_type), POINTER :: dft_control
1337 : TYPE(mp_para_env_type), POINTER :: para_env
1338 : TYPE(particle_list_type), POINTER :: particles
1339 : TYPE(pw_env_type), POINTER :: pw_env
1340 : TYPE(qs_subsys_type), POINTER :: subsys
1341 : TYPE(section_vals_type), POINTER :: cdft_constraint_section
1342 :
1343 0 : NULLIFY (cdft_constraint_section, logger, particles, dft_control, &
1344 0 : para_env, subsys, cdft_control, pw_env)
1345 0 : logger => cp_get_default_logger()
1346 :
1347 0 : CALL get_qs_env(qs_env, subsys=subsys, para_env=para_env, dft_control=dft_control, pw_env=pw_env)
1348 0 : CALL qs_subsys_get(subsys, particles=particles)
1349 0 : cdft_control => dft_control%qs_control%cdft_control
1350 0 : cdft_constraint_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%CDFT")
1351 :
1352 0 : mpi_io = .TRUE.
1353 :
1354 0 : DO igroup = 1, SIZE(cdft_control%group)
1355 :
1356 0 : middle_name = "hw_rho_total"//TRIM(ADJUSTL(cp_to_string(igroup)))
1357 : unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io, &
1358 0 : file_position="REWIND", middle_name=middle_name, extension=".cube")
1359 :
1360 : CALL cp_pw_to_cube(cdft_control%hw_rho_total, unit_nr, "CDFT Weight Function", mpi_io=mpi_io, &
1361 0 : particles=particles, stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE"))
1362 :
1363 0 : CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io)
1364 :
1365 : END DO
1366 :
1367 0 : DO igroup = 1, SIZE(cdft_control%group)
1368 :
1369 0 : middle_name = "hw_rho_total_constraint_"//TRIM(ADJUSTL(cp_to_string(igroup)))
1370 : unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io, &
1371 0 : file_position="REWIND", middle_name=middle_name, extension=".cube")
1372 :
1373 : CALL cp_pw_to_cube(cdft_control%group(igroup)%hw_rho_total_constraint, unit_nr, &
1374 : "CDFT Weight Function", mpi_io=mpi_io, particles=particles, &
1375 0 : stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE"))
1376 :
1377 0 : CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io)
1378 :
1379 : END DO
1380 :
1381 0 : DO igroup = 1, SIZE(cdft_control%group)
1382 0 : DO iatom = 1, (cdft_control%natoms)
1383 :
1384 0 : middle_name = "hw_rho_atomic_"//TRIM(ADJUSTL(cp_to_string(iatom)))
1385 : unit_nr = cp_print_key_unit_nr(logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io, &
1386 0 : file_position="REWIND", middle_name=middle_name, extension=".cube")
1387 :
1388 : CALL cp_pw_to_cube(cdft_control%group(igroup)%hw_rho_atomic(iatom), unit_nr, &
1389 : "CDFT Weight Function", mpi_io=mpi_io, particles=particles, &
1390 0 : stride=section_get_ivals(cdft_constraint_section, "PROGRAM_RUN_INFO%WEIGHT_FUNCTION%STRIDE"))
1391 :
1392 0 : CALL cp_print_key_finished_output(unit_nr, logger, cdft_constraint_section, "PROGRAM_RUN_INFO", mpi_io=mpi_io)
1393 :
1394 : END DO
1395 : END DO
1396 :
1397 0 : END SUBROUTINE cdft_print_hirshfeld_density
1398 :
1399 : END MODULE qs_cdft_utils
|
