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