(git:98357aa)
Loading...
Searching...
No Matches
qs_resp.F
Go to the documentation of this file.
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 provides a resp fit for gas phase systems
10!> \par History
11!> created
12!> Dorothea Golze [06.2012] (1) extension to periodic systems
13!> (2) re-structured the code
14!> \author Joost VandeVondele (02.2007)
15! **************************************************************************************************
16MODULE qs_resp
17 USE atomic_charges, ONLY: print_atomic_charges
18 USE atomic_kind_types, ONLY: atomic_kind_type,&
19 get_atomic_kind
20 USE bibliography, ONLY: campana2009,&
21 golze2015,&
22 rappe1992,&
23 cite_reference
24 USE cell_types, ONLY: cell_type,&
25 get_cell,&
26 pbc,&
27 use_perd_none,&
28 use_perd_xyz
29 USE cp_control_types, ONLY: dft_control_type
30 USE cp_log_handling, ONLY: cp_get_default_logger,&
31 cp_logger_type
32 USE cp_output_handling, ONLY: cp_p_file,&
33 cp_print_key_finished_output,&
34 cp_print_key_generate_filename,&
35 cp_print_key_should_output,&
36 cp_print_key_unit_nr
37 USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
38 USE cp_units, ONLY: cp_unit_from_cp2k,&
39 cp_unit_to_cp2k
40 USE input_constants, ONLY: do_resp_minus_x_dir,&
41 do_resp_minus_y_dir,&
42 do_resp_minus_z_dir,&
43 do_resp_x_dir,&
44 do_resp_y_dir,&
45 do_resp_z_dir,&
46 use_cambridge_vdw_radii,&
47 use_uff_vdw_radii
48 USE input_section_types, ONLY: section_get_ivals,&
49 section_get_lval,&
50 section_vals_get,&
51 section_vals_get_subs_vals,&
52 section_vals_type,&
53 section_vals_val_get
54 USE kahan_sum, ONLY: accurate_sum
55 USE kinds, ONLY: default_path_length,&
56 default_string_length,&
57 dp
58 USE machine, ONLY: m_flush
59 USE mathconstants, ONLY: pi
60 USE memory_utilities, ONLY: reallocate
61 USE message_passing, ONLY: mp_para_env_type,&
62 mp_request_type
63 USE particle_list_types, ONLY: particle_list_type
64 USE particle_types, ONLY: particle_type
65 USE periodic_table, ONLY: get_ptable_info
66 USE pw_env_types, ONLY: pw_env_get,&
67 pw_env_type
68 USE pw_methods, ONLY: pw_copy,&
69 pw_scale,&
70 pw_transfer,&
71 pw_zero
72 USE pw_poisson_methods, ONLY: pw_poisson_solve
73 USE pw_poisson_types, ONLY: pw_poisson_type
74 USE pw_pool_types, ONLY: pw_pool_type
75 USE pw_types, ONLY: pw_c1d_gs_type,&
76 pw_r3d_rs_type
77 USE qs_collocate_density, ONLY: calculate_rho_resp_all,&
78 calculate_rho_resp_single
79 USE qs_environment_types, ONLY: get_qs_env,&
80 qs_environment_type,&
81 set_qs_env
82 USE qs_kind_types, ONLY: qs_kind_type
83 USE qs_subsys_types, ONLY: qs_subsys_get,&
84 qs_subsys_type
85 USE uff_vdw_radii_table, ONLY: get_uff_vdw_radius
86#include "./base/base_uses.f90"
87
88 IMPLICIT NONE
89
90 PRIVATE
91
92! *** Global parameters ***
93
94 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_resp'
95
96 PUBLIC :: resp_fit
97
98 TYPE resp_type
99 LOGICAL :: equal_charges = .false., itc = .false., &
100 molecular_sys = .false., rheavies = .false., &
101 use_repeat_method = .false.
102 INTEGER :: nres = -1, ncons = -1, &
103 nrest_sec = -1, ncons_sec = -1, &
104 npoints = -1, stride(3) = -1, my_fit = -1, &
105 npoints_proc = -1, &
106 auto_vdw_radii_table = -1
107 INTEGER, DIMENSION(:), POINTER :: atom_surf_list => null()
108 INTEGER, DIMENSION(:, :), POINTER :: fitpoints => null()
109 REAL(KIND=dp) :: rheavies_strength = -1.0_dp, &
110 length = -1.0_dp, eta = -1.0_dp, &
111 sum_vhartree = -1.0_dp, offset = -1.0_dp
112 REAL(KIND=dp), DIMENSION(3) :: box_hi = -1.0_dp, box_low = -1.0_dp
113 REAL(KIND=dp), DIMENSION(:), POINTER :: rmin_kind => null(), &
114 rmax_kind => null()
115 REAL(KIND=dp), DIMENSION(:), POINTER :: range_surf => null()
116 REAL(KIND=dp), DIMENSION(:), POINTER :: rhs => null()
117 REAL(KIND=dp), DIMENSION(:), POINTER :: sum_vpot => null()
118 REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix => null()
119 END TYPE resp_type
120
121 TYPE resp_p_type
122 TYPE(resp_type), POINTER :: p_resp => null()
123 END TYPE resp_p_type
124
125CONTAINS
126
127! **************************************************************************************************
128!> \brief performs resp fit and generates RESP charges
129!> \param qs_env the qs environment
130! **************************************************************************************************
131 SUBROUTINE resp_fit(qs_env)
132 TYPE(qs_environment_type), POINTER :: qs_env
133
134 CHARACTER(len=*), PARAMETER :: routinen = 'resp_fit'
135
136 INTEGER :: handle, info, my_per, natom, nvar, &
137 output_unit
138 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
139 LOGICAL :: has_resp
140 REAL(kind=dp), DIMENSION(:), POINTER :: rhs_to_save
141 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
142 TYPE(cell_type), POINTER :: cell
143 TYPE(cp_logger_type), POINTER :: logger
144 TYPE(dft_control_type), POINTER :: dft_control
145 TYPE(particle_list_type), POINTER :: particles
146 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
147 TYPE(qs_subsys_type), POINTER :: subsys
148 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys
149 TYPE(resp_type), POINTER :: resp_env
150 TYPE(section_vals_type), POINTER :: cons_section, input, poisson_section, &
151 resp_section, rest_section
152
153 CALL timeset(routinen, handle)
154
155 NULLIFY (logger, atomic_kind_set, cell, subsys, particles, particle_set, input, &
156 resp_section, cons_section, rest_section, poisson_section, resp_env, rep_sys)
157
158 cpassert(ASSOCIATED(qs_env))
159
160 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, input=input, &
161 subsys=subsys, particle_set=particle_set, cell=cell)
162 resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
163 CALL section_vals_get(resp_section, explicit=has_resp)
164
165 IF (has_resp) THEN
166 logger => cp_get_default_logger()
167 poisson_section => section_vals_get_subs_vals(input, "DFT%POISSON")
168 CALL section_vals_val_get(poisson_section, "PERIODIC", i_val=my_per)
169 CALL create_resp_type(resp_env, rep_sys)
170 !initialize the RESP fitting, get all the keywords
171 CALL init_resp(resp_env, rep_sys, subsys, atomic_kind_set, &
172 cell, resp_section, cons_section, rest_section)
173
174 !print info
175 CALL print_resp_parameter_info(qs_env, resp_env, rep_sys, my_per)
176
177 CALL qs_subsys_get(subsys, particles=particles)
178 natom = particles%n_els
179 nvar = natom + resp_env%ncons
180
181 CALL resp_allocate(resp_env, natom, nvar)
182 ALLOCATE (ipiv(nvar))
183 ipiv = 0
184
185 ! calculate the matrix and the vector rhs
186 SELECT CASE (my_per)
187 CASE (use_perd_none)
188 CALL calc_resp_matrix_nonper(qs_env, resp_env, atomic_kind_set, particles, cell, &
189 resp_env%matrix, resp_env%rhs, natom)
190 CASE (use_perd_xyz)
191 CALL cite_reference(golze2015)
192 IF (resp_env%use_repeat_method) CALL cite_reference(campana2009)
193 CALL calc_resp_matrix_periodic(qs_env, resp_env, rep_sys, particles, cell, natom)
194 CASE DEFAULT
195 CALL cp_abort(__location__, &
196 "RESP charges only implemented for nonperiodic systems"// &
197 " or XYZ periodicity!")
198 END SELECT
199
200 output_unit = cp_print_key_unit_nr(logger, resp_section, "PRINT%PROGRAM_RUN_INFO", &
201 extension=".resp")
202 IF (output_unit > 0) THEN
203 WRITE (output_unit, '(T3,A,T69,I12)') "Number of fitting points "// &
204 "found: ", resp_env%npoints
205 WRITE (output_unit, '()')
206 END IF
207
208 !adding restraints and constraints
209 CALL add_restraints_and_constraints(qs_env, resp_env, rest_section, &
210 subsys, natom, cons_section, particle_set)
211
212 !solve system for the values of the charges and the lagrangian multipliers
213 CALL dgetrf(nvar, nvar, resp_env%matrix, nvar, ipiv, info)
214 cpassert(info == 0)
215
216 CALL dgetrs('N', nvar, 1, resp_env%matrix, nvar, ipiv, resp_env%rhs, nvar, info)
217 cpassert(info == 0)
218
219 IF (resp_env%use_repeat_method) resp_env%offset = resp_env%rhs(natom + 1)
220 CALL print_resp_charges(qs_env, resp_env, output_unit, natom)
221 CALL print_fitting_points(qs_env, resp_env)
222 CALL print_pot_from_resp_charges(qs_env, resp_env, particles, natom, output_unit)
223
224 ! In case of density functional embedding we need to save the charges to qs_env
225 NULLIFY (dft_control)
226 CALL get_qs_env(qs_env, dft_control=dft_control)
227 IF (dft_control%qs_control%ref_embed_subsys) THEN
228 ALLOCATE (rhs_to_save(SIZE(resp_env%rhs)))
229 rhs_to_save = resp_env%rhs
230 CALL set_qs_env(qs_env, rhs=rhs_to_save)
231 END IF
232
233 DEALLOCATE (ipiv)
234 CALL resp_dealloc(resp_env, rep_sys)
235 CALL cp_print_key_finished_output(output_unit, logger, resp_section, &
236 "PRINT%PROGRAM_RUN_INFO")
237
238 END IF
239
240 CALL timestop(handle)
241
242 END SUBROUTINE resp_fit
243
244! **************************************************************************************************
245!> \brief creates the resp_type structure
246!> \param resp_env the resp environment
247!> \param rep_sys structure for repeating input sections defining fit points
248! **************************************************************************************************
249 SUBROUTINE create_resp_type(resp_env, rep_sys)
250 TYPE(resp_type), POINTER :: resp_env
251 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys
252
253 IF (ASSOCIATED(resp_env)) CALL resp_dealloc(resp_env, rep_sys)
254 ALLOCATE (resp_env)
255
256 NULLIFY (resp_env%matrix, &
257 resp_env%fitpoints, &
258 resp_env%rmin_kind, &
259 resp_env%rmax_kind, &
260 resp_env%rhs, &
261 resp_env%sum_vpot)
262
263 resp_env%equal_charges = .false.
264 resp_env%itc = .false.
265 resp_env%molecular_sys = .false.
266 resp_env%rheavies = .false.
267 resp_env%use_repeat_method = .false.
268
269 resp_env%box_hi = 0.0_dp
270 resp_env%box_low = 0.0_dp
271
272 resp_env%ncons = 0
273 resp_env%ncons_sec = 0
274 resp_env%nres = 0
275 resp_env%nrest_sec = 0
276 resp_env%npoints = 0
277 resp_env%npoints_proc = 0
278 resp_env%auto_vdw_radii_table = use_cambridge_vdw_radii
279
280 END SUBROUTINE create_resp_type
281
282! **************************************************************************************************
283!> \brief allocates the resp
284!> \param resp_env the resp environment
285!> \param natom ...
286!> \param nvar ...
287! **************************************************************************************************
288 SUBROUTINE resp_allocate(resp_env, natom, nvar)
289 TYPE(resp_type), POINTER :: resp_env
290 INTEGER, INTENT(IN) :: natom, nvar
291
292 IF (.NOT. ASSOCIATED(resp_env%matrix)) THEN
293 ALLOCATE (resp_env%matrix(nvar, nvar))
294 END IF
295 IF (.NOT. ASSOCIATED(resp_env%rhs)) THEN
296 ALLOCATE (resp_env%rhs(nvar))
297 END IF
298 IF (.NOT. ASSOCIATED(resp_env%sum_vpot)) THEN
299 ALLOCATE (resp_env%sum_vpot(natom))
300 END IF
301 resp_env%matrix = 0.0_dp
302 resp_env%rhs = 0.0_dp
303 resp_env%sum_vpot = 0.0_dp
304
305 END SUBROUTINE resp_allocate
306
307! **************************************************************************************************
308!> \brief deallocates the resp_type structure
309!> \param resp_env the resp environment
310!> \param rep_sys structure for repeating input sections defining fit points
311! **************************************************************************************************
312 SUBROUTINE resp_dealloc(resp_env, rep_sys)
313 TYPE(resp_type), POINTER :: resp_env
314 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys
315
316 INTEGER :: i
317
318 IF (ASSOCIATED(resp_env)) THEN
319 IF (ASSOCIATED(resp_env%matrix)) THEN
320 DEALLOCATE (resp_env%matrix)
321 END IF
322 IF (ASSOCIATED(resp_env%rhs)) THEN
323 DEALLOCATE (resp_env%rhs)
324 END IF
325 IF (ASSOCIATED(resp_env%sum_vpot)) THEN
326 DEALLOCATE (resp_env%sum_vpot)
327 END IF
328 IF (ASSOCIATED(resp_env%fitpoints)) THEN
329 DEALLOCATE (resp_env%fitpoints)
330 END IF
331 IF (ASSOCIATED(resp_env%rmin_kind)) THEN
332 DEALLOCATE (resp_env%rmin_kind)
333 END IF
334 IF (ASSOCIATED(resp_env%rmax_kind)) THEN
335 DEALLOCATE (resp_env%rmax_kind)
336 END IF
337 DEALLOCATE (resp_env)
338 END IF
339 IF (ASSOCIATED(rep_sys)) THEN
340 DO i = 1, SIZE(rep_sys)
341 DEALLOCATE (rep_sys(i)%p_resp%atom_surf_list)
342 DEALLOCATE (rep_sys(i)%p_resp)
343 END DO
344 DEALLOCATE (rep_sys)
345 END IF
346
347 END SUBROUTINE resp_dealloc
348
349! **************************************************************************************************
350!> \brief initializes the resp fit. Getting the parameters
351!> \param resp_env the resp environment
352!> \param rep_sys structure for repeating input sections defining fit points
353!> \param subsys ...
354!> \param atomic_kind_set ...
355!> \param cell parameters related to the simulation cell
356!> \param resp_section resp section
357!> \param cons_section constraints section, part of resp section
358!> \param rest_section restraints section, part of resp section
359! **************************************************************************************************
360 SUBROUTINE init_resp(resp_env, rep_sys, subsys, atomic_kind_set, &
361 cell, resp_section, cons_section, rest_section)
362
363 TYPE(resp_type), POINTER :: resp_env
364 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys
365 TYPE(qs_subsys_type), POINTER :: subsys
366 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
367 TYPE(cell_type), POINTER :: cell
368 TYPE(section_vals_type), POINTER :: resp_section, cons_section, rest_section
369
370 CHARACTER(len=*), PARAMETER :: routinen = 'init_resp'
371
372 INTEGER :: handle, i, nrep
373 INTEGER, DIMENSION(:), POINTER :: atom_list_cons, my_stride
374 LOGICAL :: explicit
375 TYPE(section_vals_type), POINTER :: slab_section, sphere_section
376
377 CALL timeset(routinen, handle)
378
379 NULLIFY (atom_list_cons, my_stride, sphere_section, slab_section)
380
381 ! get the subsections
382 sphere_section => section_vals_get_subs_vals(resp_section, "SPHERE_SAMPLING")
383 slab_section => section_vals_get_subs_vals(resp_section, "SLAB_SAMPLING")
384 cons_section => section_vals_get_subs_vals(resp_section, "CONSTRAINT")
385 rest_section => section_vals_get_subs_vals(resp_section, "RESTRAINT")
386
387 ! get the general keywords
388 CALL section_vals_val_get(resp_section, "INTEGER_TOTAL_CHARGE", &
389 l_val=resp_env%itc)
390 IF (resp_env%itc) resp_env%ncons = resp_env%ncons + 1
391
392 CALL section_vals_val_get(resp_section, "RESTRAIN_HEAVIES_TO_ZERO", &
393 l_val=resp_env%rheavies)
394 IF (resp_env%rheavies) THEN
395 CALL section_vals_val_get(resp_section, "RESTRAIN_HEAVIES_STRENGTH", &
396 r_val=resp_env%rheavies_strength)
397 END IF
398 CALL section_vals_val_get(resp_section, "STRIDE", i_vals=my_stride)
399 IF (SIZE(my_stride) /= 1 .AND. SIZE(my_stride) /= 3) THEN
400 CALL cp_abort(__location__, "STRIDE keyword can accept only 1 (the same for X,Y,Z) "// &
401 "or 3 values. Correct your input file.")
402 END IF
403 IF (SIZE(my_stride) == 1) THEN
404 DO i = 1, 3
405 resp_env%stride(i) = my_stride(1)
406 END DO
407 ELSE
408 resp_env%stride = my_stride(1:3)
409 END IF
410 CALL section_vals_val_get(resp_section, "WIDTH", r_val=resp_env%eta)
411
412 ! get if the user wants to use REPEAT method
413 CALL section_vals_val_get(resp_section, "USE_REPEAT_METHOD", &
414 l_val=resp_env%use_repeat_method)
415 IF (resp_env%use_repeat_method) THEN
416 resp_env%ncons = resp_env%ncons + 1
417 ! restrain heavies should be off
418 resp_env%rheavies = .false.
419 END IF
420
421 ! get and set the parameters for molecular (non-surface) systems
422 ! this must come after the repeat settings being set
423 CALL get_parameter_molecular_sys(resp_env, sphere_section, cell, &
424 atomic_kind_set)
425
426 ! get the parameter for periodic/surface systems
427 CALL section_vals_get(slab_section, explicit=explicit, n_repetition=nrep)
428 IF (explicit) THEN
429 IF (resp_env%molecular_sys) THEN
430 CALL cp_abort(__location__, &
431 "You can only use either SPHERE_SAMPLING or SLAB_SAMPLING, but "// &
432 "not both.")
433 END IF
434 ALLOCATE (rep_sys(nrep))
435 DO i = 1, nrep
436 ALLOCATE (rep_sys(i)%p_resp)
437 NULLIFY (rep_sys(i)%p_resp%range_surf, rep_sys(i)%p_resp%atom_surf_list)
438 CALL section_vals_val_get(slab_section, "RANGE", r_vals=rep_sys(i)%p_resp%range_surf, &
439 i_rep_section=i)
440 CALL section_vals_val_get(slab_section, "LENGTH", r_val=rep_sys(i)%p_resp%length, &
441 i_rep_section=i)
442 CALL section_vals_val_get(slab_section, "SURF_DIRECTION", &
443 i_rep_section=i, i_val=rep_sys(i)%p_resp%my_fit)
444 IF (any(rep_sys(i)%p_resp%range_surf < 0.0_dp)) THEN
445 cpabort("Numbers in RANGE in SLAB_SAMPLING cannot be negative.")
446 END IF
447 IF (rep_sys(i)%p_resp%length <= epsilon(0.0_dp)) THEN
448 cpabort("Parameter LENGTH in SLAB_SAMPLING has to be larger than zero.")
449 END IF
450 !list of atoms specifying the surface
451 CALL build_atom_list(slab_section, subsys, rep_sys(i)%p_resp%atom_surf_list, rep=i)
452 END DO
453 END IF
454
455 ! get the parameters for the constraint and restraint sections
456 CALL section_vals_get(cons_section, explicit=explicit)
457 IF (explicit) THEN
458 CALL section_vals_get(cons_section, n_repetition=resp_env%ncons_sec)
459 DO i = 1, resp_env%ncons_sec
460 CALL section_vals_val_get(cons_section, "EQUAL_CHARGES", &
461 l_val=resp_env%equal_charges, explicit=explicit)
462 IF (.NOT. explicit) cycle
463 CALL build_atom_list(cons_section, subsys, atom_list_cons, i)
464 !instead of using EQUAL_CHARGES the constraint sections could be repeated
465 resp_env%ncons = resp_env%ncons + SIZE(atom_list_cons) - 2
466 DEALLOCATE (atom_list_cons)
467 END DO
468 END IF
469 CALL section_vals_get(rest_section, explicit=explicit)
470 IF (explicit) THEN
471 CALL section_vals_get(rest_section, n_repetition=resp_env%nrest_sec)
472 END IF
473 resp_env%ncons = resp_env%ncons + resp_env%ncons_sec
474 resp_env%nres = resp_env%nres + resp_env%nrest_sec
475
476 CALL timestop(handle)
477
478 END SUBROUTINE init_resp
479
480! **************************************************************************************************
481!> \brief getting the parameters for nonperiodic/non-surface systems
482!> \param resp_env the resp environment
483!> \param sphere_section input section setting parameters for sampling
484!> fitting in spheres around the atom
485!> \param cell parameters related to the simulation cell
486!> \param atomic_kind_set ...
487! **************************************************************************************************
488 SUBROUTINE get_parameter_molecular_sys(resp_env, sphere_section, cell, &
489 atomic_kind_set)
490
491 TYPE(resp_type), POINTER :: resp_env
492 TYPE(section_vals_type), POINTER :: sphere_section
493 TYPE(cell_type), POINTER :: cell
494 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
495
496 CHARACTER(LEN=2) :: symbol
497 CHARACTER(LEN=default_string_length) :: missing_rmax, missing_rmin
498 CHARACTER(LEN=default_string_length), &
499 DIMENSION(:), POINTER :: tmpstringlist
500 INTEGER :: ikind, j, kind_number, n_rmax_missing, &
501 n_rmin_missing, nkind, nrep_rmax, &
502 nrep_rmin, z
503 LOGICAL :: explicit, has_rmax, has_rmin
504 LOGICAL, ALLOCATABLE, DIMENSION(:) :: rmax_is_set, rmin_is_set
505 REAL(kind=dp) :: auto_rmax_scale, auto_rmin_scale, rmax, &
506 rmin
507 REAL(kind=dp), DIMENSION(3, 3) :: hmat
508 TYPE(atomic_kind_type), POINTER :: atomic_kind
509
510 nrep_rmin = 0
511 nrep_rmax = 0
512 nkind = SIZE(atomic_kind_set)
513
514 has_rmin = .false.
515 has_rmax = .false.
516
517 CALL section_vals_get(sphere_section, explicit=explicit)
518 IF (explicit) THEN
519 resp_env%molecular_sys = .true.
520 CALL section_vals_val_get(sphere_section, "AUTO_VDW_RADII_TABLE", &
521 i_val=resp_env%auto_vdw_radii_table)
522 CALL section_vals_val_get(sphere_section, "AUTO_RMIN_SCALE", r_val=auto_rmin_scale)
523 CALL section_vals_val_get(sphere_section, "AUTO_RMAX_SCALE", r_val=auto_rmax_scale)
524 CALL section_vals_val_get(sphere_section, "RMIN", explicit=has_rmin, r_val=rmin)
525 CALL section_vals_val_get(sphere_section, "RMAX", explicit=has_rmax, r_val=rmax)
526 CALL section_vals_val_get(sphere_section, "RMIN_KIND", n_rep_val=nrep_rmin)
527 CALL section_vals_val_get(sphere_section, "RMAX_KIND", n_rep_val=nrep_rmax)
528 ALLOCATE (resp_env%rmin_kind(nkind))
529 ALLOCATE (resp_env%rmax_kind(nkind))
530 resp_env%rmin_kind = 0.0_dp
531 resp_env%rmax_kind = 0.0_dp
532 ALLOCATE (rmin_is_set(nkind))
533 ALLOCATE (rmax_is_set(nkind))
534 rmin_is_set = .false.
535 rmax_is_set = .false.
536 ! define rmin_kind and rmax_kind to predefined vdW radii
537 DO ikind = 1, nkind
538 atomic_kind => atomic_kind_set(ikind)
539 CALL get_atomic_kind(atomic_kind, &
540 element_symbol=symbol, &
541 kind_number=kind_number, &
542 z=z)
543 SELECT CASE (resp_env%auto_vdw_radii_table)
544 CASE (use_cambridge_vdw_radii)
545 CALL get_ptable_info(symbol, vdw_radius=resp_env%rmin_kind(kind_number))
546 rmin_is_set(kind_number) = .true.
547 CASE (use_uff_vdw_radii)
548 CALL cite_reference(rappe1992)
549 CALL get_uff_vdw_radius(z, radius=resp_env%rmin_kind(kind_number), &
550 found=rmin_is_set(kind_number))
551 CASE DEFAULT
552 CALL get_ptable_info(symbol, vdw_radius=resp_env%rmin_kind(kind_number))
553 rmin_is_set(kind_number) = .true.
554 END SELECT
555 IF (rmin_is_set(kind_number)) THEN
556 resp_env%rmin_kind(kind_number) = cp_unit_to_cp2k(resp_env%rmin_kind(kind_number), &
557 "angstrom")
558 resp_env%rmin_kind(kind_number) = resp_env%rmin_kind(kind_number)*auto_rmin_scale
559 ! set RMAX_KIND accourding by scaling RMIN_KIND
560 resp_env%rmax_kind(kind_number) = &
561 max(resp_env%rmin_kind(kind_number), &
562 resp_env%rmin_kind(kind_number)*auto_rmax_scale)
563 rmax_is_set(kind_number) = .true.
564 END IF
565 END DO
566 ! if RMIN or RMAX are present, overwrite the rmin_kind(:) and
567 ! rmax_kind(:) to those values
568 IF (has_rmin) THEN
569 resp_env%rmin_kind = rmin
570 rmin_is_set = .true.
571 END IF
572 IF (has_rmax) THEN
573 resp_env%rmax_kind = rmax
574 rmax_is_set = .true.
575 END IF
576 ! if RMIN_KIND's or RMAX_KIND's are present, overwrite the
577 ! rmin_kinds(:) or rmax_kind(:) to those values
578 DO j = 1, nrep_rmin
579 CALL section_vals_val_get(sphere_section, "RMIN_KIND", i_rep_val=j, &
580 c_vals=tmpstringlist)
581 DO ikind = 1, nkind
582 atomic_kind => atomic_kind_set(ikind)
583 CALL get_atomic_kind(atomic_kind, element_symbol=symbol, kind_number=kind_number)
584 IF (trim(tmpstringlist(2)) == trim(symbol)) THEN
585 READ (tmpstringlist(1), *) resp_env%rmin_kind(kind_number)
586 resp_env%rmin_kind(kind_number) = &
587 cp_unit_to_cp2k(resp_env%rmin_kind(kind_number), &
588 "angstrom")
589 rmin_is_set(kind_number) = .true.
590 END IF
591 END DO
592 END DO
593 DO j = 1, nrep_rmax
594 CALL section_vals_val_get(sphere_section, "RMAX_KIND", i_rep_val=j, &
595 c_vals=tmpstringlist)
596 DO ikind = 1, nkind
597 atomic_kind => atomic_kind_set(ikind)
598 CALL get_atomic_kind(atomic_kind, element_symbol=symbol, kind_number=kind_number)
599 IF (trim(tmpstringlist(2)) == trim(symbol)) THEN
600 READ (tmpstringlist(1), *) resp_env%rmax_kind(kind_number)
601 resp_env%rmax_kind(kind_number) = cp_unit_to_cp2k(resp_env%rmax_kind(kind_number), &
602 "angstrom")
603 rmax_is_set(kind_number) = .true.
604 END IF
605 END DO
606 END DO
607 ! check if rmin and rmax are set for each kind
608 n_rmin_missing = 0
609 n_rmax_missing = 0
610 missing_rmin = ""
611 missing_rmax = ""
612 DO ikind = 1, nkind
613 atomic_kind => atomic_kind_set(ikind)
614 CALL get_atomic_kind(atomic_kind, &
615 element_symbol=symbol, &
616 kind_number=kind_number)
617 IF (.NOT. rmin_is_set(kind_number)) THEN
618 n_rmin_missing = n_rmin_missing + 1
619 missing_rmin = trim(missing_rmin)//" "//trim(symbol)//","
620 END IF
621 IF (.NOT. rmax_is_set(kind_number)) THEN
622 n_rmax_missing = n_rmax_missing + 1
623 missing_rmax = trim(missing_rmax)//" "//trim(symbol)//","
624 END IF
625 END DO
626 IF (n_rmin_missing > 0) THEN
627 CALL cp_warn(__location__, &
628 "RMIN for the following elements are missing: "// &
629 trim(missing_rmin)// &
630 " please set these values manually using "// &
631 "RMIN_KIND in SPHERE_SAMPLING section")
632 END IF
633 IF (n_rmax_missing > 0) THEN
634 CALL cp_warn(__location__, &
635 "RMAX for the following elements are missing: "// &
636 trim(missing_rmax)// &
637 " please set these values manually using "// &
638 "RMAX_KIND in SPHERE_SAMPLING section")
639 END IF
640 IF (n_rmin_missing > 0 .OR. &
641 n_rmax_missing > 0) THEN
642 cpabort("Insufficient data for RMIN or RMAX")
643 END IF
644
645 CALL get_cell(cell=cell, h=hmat)
646 resp_env%box_hi = [hmat(1, 1), hmat(2, 2), hmat(3, 3)]
647 resp_env%box_low = 0.0_dp
648 CALL section_vals_val_get(sphere_section, "X_HI", explicit=explicit)
649 IF (explicit) CALL section_vals_val_get(sphere_section, "X_HI", &
650 r_val=resp_env%box_hi(1))
651 CALL section_vals_val_get(sphere_section, "X_LOW", explicit=explicit)
652 IF (explicit) CALL section_vals_val_get(sphere_section, "X_LOW", &
653 r_val=resp_env%box_low(1))
654 CALL section_vals_val_get(sphere_section, "Y_HI", explicit=explicit)
655 IF (explicit) CALL section_vals_val_get(sphere_section, "Y_HI", &
656 r_val=resp_env%box_hi(2))
657 CALL section_vals_val_get(sphere_section, "Y_LOW", explicit=explicit)
658 IF (explicit) CALL section_vals_val_get(sphere_section, "Y_LOW", &
659 r_val=resp_env%box_low(2))
660 CALL section_vals_val_get(sphere_section, "Z_HI", explicit=explicit)
661 IF (explicit) CALL section_vals_val_get(sphere_section, "Z_HI", &
662 r_val=resp_env%box_hi(3))
663 CALL section_vals_val_get(sphere_section, "Z_LOW", explicit=explicit)
664 IF (explicit) CALL section_vals_val_get(sphere_section, "Z_LOW", &
665 r_val=resp_env%box_low(3))
666
667 DEALLOCATE (rmin_is_set)
668 DEALLOCATE (rmax_is_set)
669 END IF
670
671 END SUBROUTINE get_parameter_molecular_sys
672
673! **************************************************************************************************
674!> \brief building atom lists for different sections of RESP
675!> \param section input section
676!> \param subsys ...
677!> \param atom_list list of atoms for restraints, constraints and fit point
678!> sampling for slab-like systems
679!> \param rep input section can be repeated, this param defines for which
680!> repetition of the input section the atom_list is built
681! **************************************************************************************************
682 SUBROUTINE build_atom_list(section, subsys, atom_list, rep)
683
684 TYPE(section_vals_type), POINTER :: section
685 TYPE(qs_subsys_type), POINTER :: subsys
686 INTEGER, DIMENSION(:), POINTER :: atom_list
687 INTEGER, INTENT(IN), OPTIONAL :: rep
688
689 CHARACTER(len=*), PARAMETER :: routinen = 'build_atom_list'
690
691 INTEGER :: atom_a, atom_b, handle, i, irep, j, &
692 max_index, n_var, num_atom
693 INTEGER, DIMENSION(:), POINTER :: indexes
694 LOGICAL :: index_in_range
695
696 CALL timeset(routinen, handle)
697
698 NULLIFY (indexes)
699 irep = 1
700 IF (PRESENT(rep)) irep = rep
701
702 CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
703 n_rep_val=n_var)
704 num_atom = 0
705 DO i = 1, n_var
706 CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
707 i_rep_val=i, i_vals=indexes)
708 num_atom = num_atom + SIZE(indexes)
709 END DO
710 ALLOCATE (atom_list(num_atom))
711 atom_list = 0
712 num_atom = 1
713 DO i = 1, n_var
714 CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
715 i_rep_val=i, i_vals=indexes)
716 atom_list(num_atom:num_atom + SIZE(indexes) - 1) = indexes(:)
717 num_atom = num_atom + SIZE(indexes)
718 END DO
719 !check atom list
720 num_atom = num_atom - 1
721 CALL qs_subsys_get(subsys, nparticle=max_index)
722 cpassert(SIZE(atom_list) /= 0)
723 index_in_range = (maxval(atom_list) <= max_index) &
724 .AND. (minval(atom_list) > 0)
725 cpassert(index_in_range)
726 DO i = 1, num_atom
727 DO j = i + 1, num_atom
728 atom_a = atom_list(i)
729 atom_b = atom_list(j)
730 IF (atom_a == atom_b) THEN
731 cpabort("There are atoms doubled in atom list for RESP.")
732 END IF
733 END DO
734 END DO
735
736 CALL timestop(handle)
737
738 END SUBROUTINE build_atom_list
739
740! **************************************************************************************************
741!> \brief build matrix and vector for nonperiodic RESP fitting
742!> \param qs_env the qs environment
743!> \param resp_env the resp environment
744!> \param atomic_kind_set ...
745!> \param particles ...
746!> \param cell parameters related to the simulation cell
747!> \param matrix coefficient matrix of the linear system of equations
748!> \param rhs vector of the linear system of equations
749!> \param natom number of atoms
750! **************************************************************************************************
751 SUBROUTINE calc_resp_matrix_nonper(qs_env, resp_env, atomic_kind_set, particles, &
752 cell, matrix, rhs, natom)
753
754 TYPE(qs_environment_type), POINTER :: qs_env
755 TYPE(resp_type), POINTER :: resp_env
756 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
757 TYPE(particle_list_type), POINTER :: particles
758 TYPE(cell_type), POINTER :: cell
759 REAL(kind=dp), DIMENSION(:, :), POINTER :: matrix
760 REAL(kind=dp), DIMENSION(:), POINTER :: rhs
761 INTEGER, INTENT(IN) :: natom
762
763 CHARACTER(len=*), PARAMETER :: routinen = 'calc_resp_matrix_nonper'
764
765 INTEGER :: bo(2, 3), gbo(2, 3), handle, i, ikind, &
766 jx, jy, jz, k, kind_number, l, m, &
767 nkind, now, np(3), p
768 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: not_in_range
769 REAL(kind=dp) :: delta, dh(3, 3), dvol, r(3), rmax, rmin, &
770 vec(3), vec_pbc(3), vj
771 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dist
772 REAL(kind=dp), DIMENSION(3, 3) :: hmat, hmat_inv
773 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
774 TYPE(pw_r3d_rs_type), POINTER :: v_hartree_pw
775
776 CALL timeset(routinen, handle)
777
778 NULLIFY (particle_set, v_hartree_pw)
779 delta = 1.0e-13_dp
780
781 CALL get_cell(cell=cell, h=hmat, h_inv=hmat_inv)
782
783 IF (.NOT. cell%orthorhombic) THEN
784 CALL cp_abort(__location__, &
785 "Nonperiodic solution for RESP charges only"// &
786 " implemented for orthorhombic cells!")
787 END IF
788 IF (.NOT. resp_env%molecular_sys) THEN
789 CALL cp_abort(__location__, &
790 "Nonperiodic solution for RESP charges (i.e. nonperiodic"// &
791 " Poisson solver) can only be used with section SPHERE_SAMPLING")
792 END IF
793 IF (resp_env%use_repeat_method) THEN
794 CALL cp_abort(__location__, &
795 "REPEAT method only reasonable for periodic RESP fitting")
796 END IF
797 CALL get_qs_env(qs_env, particle_set=particle_set, v_hartree_rspace=v_hartree_pw)
798
799 bo = v_hartree_pw%pw_grid%bounds_local
800 gbo = v_hartree_pw%pw_grid%bounds
801 np = v_hartree_pw%pw_grid%npts
802 dh = v_hartree_pw%pw_grid%dh
803 dvol = v_hartree_pw%pw_grid%dvol
804 nkind = SIZE(atomic_kind_set)
805
806 ALLOCATE (dist(natom))
807 ALLOCATE (not_in_range(natom, 2))
808
809 ! store fitting points to calculate the RMS and RRMS later
810 IF (.NOT. ASSOCIATED(resp_env%fitpoints)) THEN
811 now = 1000
812 ALLOCATE (resp_env%fitpoints(3, now))
813 ELSE
814 now = SIZE(resp_env%fitpoints, 2)
815 END IF
816
817 DO jz = bo(1, 3), bo(2, 3)
818 DO jy = bo(1, 2), bo(2, 2)
819 DO jx = bo(1, 1), bo(2, 1)
820 IF (.NOT. (modulo(jz, resp_env%stride(3)) == 0)) cycle
821 IF (.NOT. (modulo(jy, resp_env%stride(2)) == 0)) cycle
822 IF (.NOT. (modulo(jx, resp_env%stride(1)) == 0)) cycle
823 !bounds bo reach from -np/2 to np/2. shift of np/2 so that r(1,1,1)=(0,0,0)
824 l = jx - gbo(1, 1)
825 k = jy - gbo(1, 2)
826 p = jz - gbo(1, 3)
827 r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
828 r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
829 r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
830 IF (r(3) < resp_env%box_low(3) .OR. r(3) > resp_env%box_hi(3)) cycle
831 IF (r(2) < resp_env%box_low(2) .OR. r(2) > resp_env%box_hi(2)) cycle
832 IF (r(1) < resp_env%box_low(1) .OR. r(1) > resp_env%box_hi(1)) cycle
833 ! compute distance from the grid point to all atoms
834 not_in_range = .false.
835 DO i = 1, natom
836 vec = r - particles%els(i)%r
837 vec_pbc(1) = vec(1) - hmat(1, 1)*anint(hmat_inv(1, 1)*vec(1))
838 vec_pbc(2) = vec(2) - hmat(2, 2)*anint(hmat_inv(2, 2)*vec(2))
839 vec_pbc(3) = vec(3) - hmat(3, 3)*anint(hmat_inv(3, 3)*vec(3))
840 dist(i) = sqrt(sum(vec_pbc**2))
841 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, &
842 kind_number=kind_number)
843 DO ikind = 1, nkind
844 IF (ikind == kind_number) THEN
845 rmin = resp_env%rmin_kind(ikind)
846 rmax = resp_env%rmax_kind(ikind)
847 EXIT
848 END IF
849 END DO
850 IF (dist(i) < rmin + delta) not_in_range(i, 1) = .true.
851 IF (dist(i) > rmax - delta) not_in_range(i, 2) = .true.
852 END DO
853 ! check if the point is sufficiently close and far. if OK, we can use
854 ! the point for fitting, add/subtract 1.0E-13 to get rid of rounding errors when shifting atoms
855 IF (any(not_in_range(:, 1)) .OR. all(not_in_range(:, 2))) cycle
856 resp_env%npoints_proc = resp_env%npoints_proc + 1
857 IF (resp_env%npoints_proc > now) THEN
858 now = 2*now
859 CALL reallocate(resp_env%fitpoints, 1, 3, 1, now)
860 END IF
861 resp_env%fitpoints(1, resp_env%npoints_proc) = jx
862 resp_env%fitpoints(2, resp_env%npoints_proc) = jy
863 resp_env%fitpoints(3, resp_env%npoints_proc) = jz
864 ! correct for the fact that v_hartree is scaled by dvol, and has the opposite sign
865 IF (qs_env%qmmm) THEN
866 ! If it's a QM/MM run let's remove the contribution of the MM potential out of the Hartree pot
867 vj = -v_hartree_pw%array(jx, jy, jz)/dvol + qs_env%ks_qmmm_env%v_qmmm_rspace%array(jx, jy, jz)
868 ELSE
869 vj = -v_hartree_pw%array(jx, jy, jz)/dvol
870 END IF
871 dist(:) = 1.0_dp/dist(:)
872
873 DO i = 1, natom
874 DO m = 1, natom
875 matrix(m, i) = matrix(m, i) + 2.0_dp*dist(i)*dist(m)
876 END DO
877 rhs(i) = rhs(i) + 2.0_dp*vj*dist(i)
878 END DO
879 END DO
880 END DO
881 END DO
882
883 resp_env%npoints = resp_env%npoints_proc
884 CALL v_hartree_pw%pw_grid%para%group%sum(resp_env%npoints)
885 CALL v_hartree_pw%pw_grid%para%group%sum(matrix)
886 CALL v_hartree_pw%pw_grid%para%group%sum(rhs)
887 !weighted sum
888 matrix = matrix/resp_env%npoints
889 rhs = rhs/resp_env%npoints
890
891 DEALLOCATE (dist)
892 DEALLOCATE (not_in_range)
893
894 CALL timestop(handle)
895
896 END SUBROUTINE calc_resp_matrix_nonper
897
898! **************************************************************************************************
899!> \brief build matrix and vector for periodic RESP fitting
900!> \param qs_env the qs environment
901!> \param resp_env the resp environment
902!> \param rep_sys structure for repeating input sections defining fit points
903!> \param particles ...
904!> \param cell parameters related to the simulation cell
905!> \param natom number of atoms
906! **************************************************************************************************
907 SUBROUTINE calc_resp_matrix_periodic(qs_env, resp_env, rep_sys, particles, cell, &
908 natom)
909
910 TYPE(qs_environment_type), POINTER :: qs_env
911 TYPE(resp_type), POINTER :: resp_env
912 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys
913 TYPE(particle_list_type), POINTER :: particles
914 TYPE(cell_type), POINTER :: cell
915 INTEGER, INTENT(IN) :: natom
916
917 CHARACTER(len=*), PARAMETER :: routinen = 'calc_resp_matrix_periodic'
918
919 INTEGER :: handle, i, ip, j, jx, jy, jz
920 INTEGER, DIMENSION(3) :: periodic
921 REAL(kind=dp) :: normalize_factor
922 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: vpot
923 TYPE(mp_para_env_type), POINTER :: para_env
924 TYPE(pw_c1d_gs_type) :: rho_ga, va_gspace
925 TYPE(pw_env_type), POINTER :: pw_env
926 TYPE(pw_poisson_type), POINTER :: poisson_env
927 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
928 TYPE(pw_r3d_rs_type) :: va_rspace
929
930 CALL timeset(routinen, handle)
931
932 NULLIFY (pw_env, para_env, auxbas_pw_pool, poisson_env)
933
934 CALL get_cell(cell=cell, periodic=periodic)
935
936 IF (.NOT. all(periodic /= 0)) THEN
937 CALL cp_abort(__location__, &
938 "Periodic solution for RESP (with periodic Poisson solver)"// &
939 " can only be obtained with a cell that has XYZ periodicity")
940 END IF
941
942 CALL get_qs_env(qs_env, pw_env=pw_env, para_env=para_env)
943
944 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
945 poisson_env=poisson_env)
946 CALL auxbas_pw_pool%create_pw(rho_ga)
947 CALL auxbas_pw_pool%create_pw(va_gspace)
948 CALL auxbas_pw_pool%create_pw(va_rspace)
949
950 !get fitting points and store them in resp_env%fitpoints
951 CALL get_fitting_points(qs_env, resp_env, rep_sys, particles=particles, &
952 cell=cell)
953 ALLOCATE (vpot(resp_env%npoints_proc, natom))
954 normalize_factor = sqrt((resp_env%eta/pi)**3)
955
956 DO i = 1, natom
957 !collocate gaussian for each atom
958 CALL pw_zero(rho_ga)
959 CALL calculate_rho_resp_single(rho_ga, qs_env, resp_env%eta, i)
960 !calculate potential va and store the part needed for fitting in vpot
961 CALL pw_zero(va_gspace)
962 CALL pw_poisson_solve(poisson_env, rho_ga, vhartree=va_gspace)
963 CALL pw_zero(va_rspace)
964 CALL pw_transfer(va_gspace, va_rspace)
965 CALL pw_scale(va_rspace, normalize_factor)
966 DO ip = 1, resp_env%npoints_proc
967 jx = resp_env%fitpoints(1, ip)
968 jy = resp_env%fitpoints(2, ip)
969 jz = resp_env%fitpoints(3, ip)
970 vpot(ip, i) = va_rspace%array(jx, jy, jz)
971 END DO
972 END DO
973
974 CALL va_gspace%release()
975 CALL va_rspace%release()
976 CALL rho_ga%release()
977
978 DO i = 1, natom
979 DO j = 1, natom
980 ! calculate matrix
981 resp_env%matrix(i, j) = resp_env%matrix(i, j) + 2.0_dp*sum(vpot(:, i)*vpot(:, j))
982 END DO
983 ! calculate vector resp_env%rhs
984 CALL calculate_rhs(qs_env, resp_env, resp_env%rhs(i), vpot(:, i))
985 END DO
986
987 CALL para_env%sum(resp_env%matrix)
988 CALL para_env%sum(resp_env%rhs)
989 !weighted sum
990 resp_env%matrix = resp_env%matrix/resp_env%npoints
991 resp_env%rhs = resp_env%rhs/resp_env%npoints
992
993 ! REPEAT stuff
994 IF (resp_env%use_repeat_method) THEN
995 ! sum over selected points of single Gaussian potential vpot
996 DO i = 1, natom
997 resp_env%sum_vpot(i) = 2.0_dp*accurate_sum(vpot(:, i))/resp_env%npoints
998 END DO
999 CALL para_env%sum(resp_env%sum_vpot)
1000 CALL para_env%sum(resp_env%sum_vhartree)
1001 resp_env%sum_vhartree = 2.0_dp*resp_env%sum_vhartree/resp_env%npoints
1002 END IF
1003
1004 DEALLOCATE (vpot)
1005 CALL timestop(handle)
1006
1007 END SUBROUTINE calc_resp_matrix_periodic
1008
1009! **************************************************************************************************
1010!> \brief get RESP fitting points for the periodic fitting
1011!> \param qs_env the qs environment
1012!> \param resp_env the resp environment
1013!> \param rep_sys structure for repeating input sections defining fit points
1014!> \param particles ...
1015!> \param cell parameters related to the simulation cell
1016! **************************************************************************************************
1017 SUBROUTINE get_fitting_points(qs_env, resp_env, rep_sys, particles, cell)
1018
1019 TYPE(qs_environment_type), POINTER :: qs_env
1020 TYPE(resp_type), POINTER :: resp_env
1021 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys
1022 TYPE(particle_list_type), POINTER :: particles
1023 TYPE(cell_type), POINTER :: cell
1024
1025 CHARACTER(len=*), PARAMETER :: routinen = 'get_fitting_points'
1026
1027 INTEGER :: bo(2, 3), gbo(2, 3), handle, i, iatom, &
1028 ikind, in_x, in_y, in_z, jx, jy, jz, &
1029 k, kind_number, l, m, natom, nkind, &
1030 now, p
1031 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: not_in_range
1032 REAL(kind=dp) :: delta, dh(3, 3), r(3), rmax, rmin, &
1033 vec_pbc(3)
1034 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dist
1035 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1036 TYPE(mp_para_env_type), POINTER :: para_env
1037 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1038 TYPE(pw_r3d_rs_type), POINTER :: v_hartree_pw
1039
1040 CALL timeset(routinen, handle)
1041
1042 NULLIFY (atomic_kind_set, v_hartree_pw, para_env, particle_set)
1043 delta = 1.0e-13_dp
1044
1045 CALL get_qs_env(qs_env, &
1046 particle_set=particle_set, &
1047 atomic_kind_set=atomic_kind_set, &
1048 para_env=para_env, &
1049 v_hartree_rspace=v_hartree_pw)
1050
1051 bo = v_hartree_pw%pw_grid%bounds_local
1052 gbo = v_hartree_pw%pw_grid%bounds
1053 dh = v_hartree_pw%pw_grid%dh
1054 natom = SIZE(particles%els)
1055 nkind = SIZE(atomic_kind_set)
1056
1057 IF (.NOT. ASSOCIATED(resp_env%fitpoints)) THEN
1058 now = 1000
1059 ALLOCATE (resp_env%fitpoints(3, now))
1060 ELSE
1061 now = SIZE(resp_env%fitpoints, 2)
1062 END IF
1063
1064 ALLOCATE (dist(natom))
1065 ALLOCATE (not_in_range(natom, 2))
1066
1067 !every proc gets another bo, grid is distributed
1068 DO jz = bo(1, 3), bo(2, 3)
1069 IF (.NOT. (modulo(jz, resp_env%stride(3)) == 0)) cycle
1070 DO jy = bo(1, 2), bo(2, 2)
1071 IF (.NOT. (modulo(jy, resp_env%stride(2)) == 0)) cycle
1072 DO jx = bo(1, 1), bo(2, 1)
1073 IF (.NOT. (modulo(jx, resp_env%stride(1)) == 0)) cycle
1074 !bounds gbo reach from -np/2 to np/2. shift of np/2 so that r(1,1,1)=(0,0,0)
1075 l = jx - gbo(1, 1)
1076 k = jy - gbo(1, 2)
1077 p = jz - gbo(1, 3)
1078 r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
1079 r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
1080 r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
1081 IF (resp_env%molecular_sys) THEN
1082 not_in_range = .false.
1083 DO m = 1, natom
1084 vec_pbc = pbc(r, particles%els(m)%r, cell)
1085 dist(m) = sqrt(sum(vec_pbc**2))
1086 CALL get_atomic_kind(atomic_kind=particle_set(m)%atomic_kind, &
1087 kind_number=kind_number)
1088 DO ikind = 1, nkind
1089 IF (ikind == kind_number) THEN
1090 rmin = resp_env%rmin_kind(ikind)
1091 rmax = resp_env%rmax_kind(ikind)
1092 EXIT
1093 END IF
1094 END DO
1095 IF (dist(m) < rmin + delta) not_in_range(m, 1) = .true.
1096 IF (dist(m) > rmax - delta) not_in_range(m, 2) = .true.
1097 END DO
1098 IF (any(not_in_range(:, 1)) .OR. all(not_in_range(:, 2))) cycle
1099 ELSE
1100 DO i = 1, SIZE(rep_sys)
1101 DO m = 1, SIZE(rep_sys(i)%p_resp%atom_surf_list)
1102 in_z = 0
1103 in_y = 0
1104 in_x = 0
1105 iatom = rep_sys(i)%p_resp%atom_surf_list(m)
1106 SELECT CASE (rep_sys(i)%p_resp%my_fit)
1107 CASE (do_resp_x_dir, do_resp_y_dir, do_resp_z_dir)
1108 vec_pbc = pbc(particles%els(iatom)%r, r, cell)
1109 CASE (do_resp_minus_x_dir, do_resp_minus_y_dir, do_resp_minus_z_dir)
1110 vec_pbc = pbc(r, particles%els(iatom)%r, cell)
1111 END SELECT
1112 SELECT CASE (rep_sys(i)%p_resp%my_fit)
1113 !subtract delta=1.0E-13 to get rid of rounding errors when shifting atoms
1114 CASE (do_resp_x_dir, do_resp_minus_x_dir)
1115 IF (abs(vec_pbc(3)) < rep_sys(i)%p_resp%length - delta) in_z = 1
1116 IF (abs(vec_pbc(2)) < rep_sys(i)%p_resp%length - delta) in_y = 1
1117 IF (vec_pbc(1) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
1118 vec_pbc(1) < rep_sys(i)%p_resp%range_surf(2) - delta) in_x = 1
1119 CASE (do_resp_y_dir, do_resp_minus_y_dir)
1120 IF (abs(vec_pbc(3)) < rep_sys(i)%p_resp%length - delta) in_z = 1
1121 IF (vec_pbc(2) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
1122 vec_pbc(2) < rep_sys(i)%p_resp%range_surf(2) - delta) in_y = 1
1123 IF (abs(vec_pbc(1)) < rep_sys(i)%p_resp%length - delta) in_x = 1
1124 CASE (do_resp_z_dir, do_resp_minus_z_dir)
1125 IF (vec_pbc(3) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
1126 vec_pbc(3) < rep_sys(i)%p_resp%range_surf(2) - delta) in_z = 1
1127 IF (abs(vec_pbc(2)) < rep_sys(i)%p_resp%length - delta) in_y = 1
1128 IF (abs(vec_pbc(1)) < rep_sys(i)%p_resp%length - delta) in_x = 1
1129 END SELECT
1130 IF (in_z*in_y*in_x == 1) EXIT
1131 END DO
1132 IF (in_z*in_y*in_x == 1) EXIT
1133 END DO
1134 IF (in_z*in_y*in_x == 0) cycle
1135 END IF
1136 resp_env%npoints_proc = resp_env%npoints_proc + 1
1137 IF (resp_env%npoints_proc > now) THEN
1138 now = 2*now
1139 CALL reallocate(resp_env%fitpoints, 1, 3, 1, now)
1140 END IF
1141 resp_env%fitpoints(1, resp_env%npoints_proc) = jx
1142 resp_env%fitpoints(2, resp_env%npoints_proc) = jy
1143 resp_env%fitpoints(3, resp_env%npoints_proc) = jz
1144 END DO
1145 END DO
1146 END DO
1147
1148 resp_env%npoints = resp_env%npoints_proc
1149 CALL para_env%sum(resp_env%npoints)
1150
1151 DEALLOCATE (dist)
1152 DEALLOCATE (not_in_range)
1153
1154 CALL timestop(handle)
1155
1156 END SUBROUTINE get_fitting_points
1157
1158! **************************************************************************************************
1159!> \brief calculate vector rhs
1160!> \param qs_env the qs environment
1161!> \param resp_env the resp environment
1162!> \param rhs vector
1163!> \param vpot single gaussian potential
1164! **************************************************************************************************
1165 SUBROUTINE calculate_rhs(qs_env, resp_env, rhs, vpot)
1166
1167 TYPE(qs_environment_type), POINTER :: qs_env
1168 TYPE(resp_type), POINTER :: resp_env
1169 REAL(kind=dp), INTENT(INOUT) :: rhs
1170 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: vpot
1171
1172 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rhs'
1173
1174 INTEGER :: handle, ip, jx, jy, jz
1175 REAL(kind=dp) :: dvol
1176 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: vhartree
1177 TYPE(pw_r3d_rs_type), POINTER :: v_hartree_pw
1178
1179 CALL timeset(routinen, handle)
1180
1181 NULLIFY (v_hartree_pw)
1182 CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_pw)
1183 dvol = v_hartree_pw%pw_grid%dvol
1184 ALLOCATE (vhartree(resp_env%npoints_proc))
1185 vhartree = 0.0_dp
1186
1187 !multiply v_hartree and va_rspace and calculate the vector rhs
1188 !taking into account that v_hartree has opposite site; remove v_qmmm
1189 DO ip = 1, resp_env%npoints_proc
1190 jx = resp_env%fitpoints(1, ip)
1191 jy = resp_env%fitpoints(2, ip)
1192 jz = resp_env%fitpoints(3, ip)
1193 vhartree(ip) = -v_hartree_pw%array(jx, jy, jz)/dvol
1194 IF (qs_env%qmmm) THEN
1195 !taking into account that v_qmmm has also opposite sign
1196 vhartree(ip) = vhartree(ip) + qs_env%ks_qmmm_env%v_qmmm_rspace%array(jx, jy, jz)
1197 END IF
1198 rhs = rhs + 2.0_dp*vhartree(ip)*vpot(ip)
1199 END DO
1200
1201 IF (resp_env%use_repeat_method) THEN
1202 resp_env%sum_vhartree = accurate_sum(vhartree)
1203 END IF
1204
1205 DEALLOCATE (vhartree)
1206
1207 CALL timestop(handle)
1208
1209 END SUBROUTINE calculate_rhs
1210
1211! **************************************************************************************************
1212!> \brief print the atom coordinates and the coordinates of the fitting points
1213!> to an xyz file
1214!> \param qs_env the qs environment
1215!> \param resp_env the resp environment
1216! **************************************************************************************************
1217 SUBROUTINE print_fitting_points(qs_env, resp_env)
1218
1219 TYPE(qs_environment_type), POINTER :: qs_env
1220 TYPE(resp_type), POINTER :: resp_env
1221
1222 CHARACTER(len=*), PARAMETER :: routinen = 'print_fitting_points'
1223
1224 CHARACTER(LEN=2) :: element_symbol
1225 CHARACTER(LEN=default_path_length) :: filename
1226 INTEGER :: gbo(2, 3), handle, i, iatom, ip, jx, jy, &
1227 jz, k, l, my_pos, nobjects, &
1228 output_unit, p
1229 INTEGER, DIMENSION(:), POINTER :: tmp_npoints, tmp_size
1230 INTEGER, DIMENSION(:, :), POINTER :: tmp_points
1231 REAL(kind=dp) :: conv, dh(3, 3), r(3)
1232 TYPE(cp_logger_type), POINTER :: logger
1233 TYPE(mp_para_env_type), POINTER :: para_env
1234 TYPE(mp_request_type), DIMENSION(6) :: req
1235 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1236 TYPE(pw_r3d_rs_type), POINTER :: v_hartree_pw
1237 TYPE(section_vals_type), POINTER :: input, print_key, resp_section
1238
1239 CALL timeset(routinen, handle)
1240
1241 NULLIFY (para_env, input, logger, resp_section, print_key, particle_set, tmp_size, &
1242 tmp_points, tmp_npoints, v_hartree_pw)
1243
1244 CALL get_qs_env(qs_env, input=input, para_env=para_env, &
1245 particle_set=particle_set, v_hartree_rspace=v_hartree_pw)
1246 conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
1247 gbo = v_hartree_pw%pw_grid%bounds
1248 dh = v_hartree_pw%pw_grid%dh
1249 nobjects = SIZE(particle_set) + resp_env%npoints
1250
1251 resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1252 print_key => section_vals_get_subs_vals(resp_section, "PRINT%COORD_FIT_POINTS")
1253 logger => cp_get_default_logger()
1254 output_unit = cp_print_key_unit_nr(logger, resp_section, &
1255 "PRINT%COORD_FIT_POINTS", &
1256 extension=".xyz", &
1257 file_status="REPLACE", &
1258 file_action="WRITE", &
1259 file_form="FORMATTED")
1260
1261 IF (btest(cp_print_key_should_output(logger%iter_info, &
1262 resp_section, "PRINT%COORD_FIT_POINTS"), &
1263 cp_p_file)) THEN
1264 IF (output_unit > 0) THEN
1265 filename = cp_print_key_generate_filename(logger, &
1266 print_key, extension=".xyz", &
1267 my_local=.false.)
1268 WRITE (unit=output_unit, fmt="(I12,/)") nobjects
1269 DO iatom = 1, SIZE(particle_set)
1270 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
1271 element_symbol=element_symbol)
1272 WRITE (unit=output_unit, fmt="(A,1X,3F10.5)") element_symbol, &
1273 particle_set(iatom)%r(1:3)*conv
1274 END DO
1275 !printing points of proc which is doing the output (should be proc 0)
1276 DO ip = 1, resp_env%npoints_proc
1277 jx = resp_env%fitpoints(1, ip)
1278 jy = resp_env%fitpoints(2, ip)
1279 jz = resp_env%fitpoints(3, ip)
1280 l = jx - gbo(1, 1)
1281 k = jy - gbo(1, 2)
1282 p = jz - gbo(1, 3)
1283 r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
1284 r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
1285 r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
1286 r(:) = r(:)*conv
1287 WRITE (unit=output_unit, fmt="(A,2X,3F10.5)") "X", r(1), r(2), r(3)
1288 END DO
1289 END IF
1290
1291 ALLOCATE (tmp_size(1))
1292 ALLOCATE (tmp_npoints(1))
1293
1294 !sending data of all other procs to proc which makes the output (proc 0)
1295 IF (output_unit > 0) THEN
1296 my_pos = para_env%mepos
1297 DO i = 1, para_env%num_pe
1298 IF (my_pos == i - 1) cycle
1299 CALL para_env%irecv(msgout=tmp_size, source=i - 1, &
1300 request=req(1))
1301 CALL req(1)%wait()
1302 ALLOCATE (tmp_points(3, tmp_size(1)))
1303 CALL para_env%irecv(msgout=tmp_points, source=i - 1, &
1304 request=req(3))
1305 CALL req(3)%wait()
1306 CALL para_env%irecv(msgout=tmp_npoints, source=i - 1, &
1307 request=req(5))
1308 CALL req(5)%wait()
1309 DO ip = 1, tmp_npoints(1)
1310 jx = tmp_points(1, ip)
1311 jy = tmp_points(2, ip)
1312 jz = tmp_points(3, ip)
1313 l = jx - gbo(1, 1)
1314 k = jy - gbo(1, 2)
1315 p = jz - gbo(1, 3)
1316 r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
1317 r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
1318 r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
1319 r(:) = r(:)*conv
1320 WRITE (unit=output_unit, fmt="(A,2X,3F10.5)") "X", r(1), r(2), r(3)
1321 END DO
1322 DEALLOCATE (tmp_points)
1323 END DO
1324 ELSE
1325 tmp_size(1) = SIZE(resp_env%fitpoints, 2)
1326 !para_env%source should be 0
1327 CALL para_env%isend(msgin=tmp_size, dest=para_env%source, &
1328 request=req(2))
1329 CALL req(2)%wait()
1330 CALL para_env%isend(msgin=resp_env%fitpoints, dest=para_env%source, &
1331 request=req(4))
1332 CALL req(4)%wait()
1333 tmp_npoints(1) = resp_env%npoints_proc
1334 CALL para_env%isend(msgin=tmp_npoints, dest=para_env%source, &
1335 request=req(6))
1336 CALL req(6)%wait()
1337 END IF
1338
1339 DEALLOCATE (tmp_size)
1340 DEALLOCATE (tmp_npoints)
1341 END IF
1342
1343 CALL cp_print_key_finished_output(output_unit, logger, resp_section, &
1344 "PRINT%COORD_FIT_POINTS")
1345
1346 CALL timestop(handle)
1347
1348 END SUBROUTINE print_fitting_points
1349
1350! **************************************************************************************************
1351!> \brief add restraints and constraints
1352!> \param qs_env the qs environment
1353!> \param resp_env the resp environment
1354!> \param rest_section input section for restraints
1355!> \param subsys ...
1356!> \param natom number of atoms
1357!> \param cons_section input section for constraints
1358!> \param particle_set ...
1359! **************************************************************************************************
1360 SUBROUTINE add_restraints_and_constraints(qs_env, resp_env, rest_section, &
1361 subsys, natom, cons_section, particle_set)
1362
1363 TYPE(qs_environment_type), POINTER :: qs_env
1364 TYPE(resp_type), POINTER :: resp_env
1365 TYPE(section_vals_type), POINTER :: rest_section
1366 TYPE(qs_subsys_type), POINTER :: subsys
1367 INTEGER, INTENT(IN) :: natom
1368 TYPE(section_vals_type), POINTER :: cons_section
1369 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1370
1371 CHARACTER(len=*), PARAMETER :: routinen = 'add_restraints_and_constraints'
1372
1373 INTEGER :: handle, i, k, m, ncons_v, z
1374 INTEGER, DIMENSION(:), POINTER :: atom_list_cons, atom_list_res
1375 LOGICAL :: explicit_coeff
1376 REAL(kind=dp) :: my_atom_coef(2), strength, TARGET
1377 REAL(kind=dp), DIMENSION(:), POINTER :: atom_coef
1378 TYPE(dft_control_type), POINTER :: dft_control
1379
1380 CALL timeset(routinen, handle)
1381
1382 NULLIFY (atom_coef, atom_list_res, atom_list_cons, dft_control)
1383
1384 CALL get_qs_env(qs_env, dft_control=dft_control)
1385
1386 !*** add the restraints
1387 DO i = 1, resp_env%nrest_sec
1388 CALL section_vals_val_get(rest_section, "TARGET", i_rep_section=i, r_val=TARGET)
1389 CALL section_vals_val_get(rest_section, "STRENGTH", i_rep_section=i, r_val=strength)
1390 CALL build_atom_list(rest_section, subsys, atom_list_res, i)
1391 CALL section_vals_val_get(rest_section, "ATOM_COEF", i_rep_section=i, explicit=explicit_coeff)
1392 IF (explicit_coeff) THEN
1393 CALL section_vals_val_get(rest_section, "ATOM_COEF", i_rep_section=i, r_vals=atom_coef)
1394 cpassert(SIZE(atom_list_res) == SIZE(atom_coef))
1395 END IF
1396 DO m = 1, SIZE(atom_list_res)
1397 IF (explicit_coeff) THEN
1398 DO k = 1, SIZE(atom_list_res)
1399 resp_env%matrix(atom_list_res(m), atom_list_res(k)) = &
1400 resp_env%matrix(atom_list_res(m), atom_list_res(k)) + &
1401 atom_coef(m)*atom_coef(k)*2.0_dp*strength
1402 END DO
1403 resp_env%rhs(atom_list_res(m)) = resp_env%rhs(atom_list_res(m)) + &
1404 2.0_dp*TARGET*strength*atom_coef(m)
1405 ELSE
1406 resp_env%matrix(atom_list_res(m), atom_list_res(m)) = &
1407 resp_env%matrix(atom_list_res(m), atom_list_res(m)) + &
1408 2.0_dp*strength
1409 resp_env%rhs(atom_list_res(m)) = resp_env%rhs(atom_list_res(m)) + &
1410 2.0_dp*TARGET*strength
1411 END IF
1412 END DO
1413 DEALLOCATE (atom_list_res)
1414 END DO
1415
1416 ! if heavies are restrained to zero, add these as well
1417 IF (resp_env%rheavies) THEN
1418 DO i = 1, natom
1419 CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, z=z)
1420 IF (z /= 1) THEN
1421 resp_env%matrix(i, i) = resp_env%matrix(i, i) + 2.0_dp*resp_env%rheavies_strength
1422 END IF
1423 END DO
1424 END IF
1425
1426 !*** add the constraints
1427 ncons_v = 0
1428 ncons_v = ncons_v + natom
1429
1430 ! REPEAT charges: treat the offset like a constraint
1431 IF (resp_env%use_repeat_method) THEN
1432 ncons_v = ncons_v + 1
1433 resp_env%matrix(1:natom, ncons_v) = resp_env%sum_vpot(1:natom)
1434 resp_env%matrix(ncons_v, 1:natom) = resp_env%sum_vpot(1:natom)
1435 resp_env%matrix(ncons_v, ncons_v) = 2.0_dp
1436 resp_env%rhs(ncons_v) = resp_env%sum_vhartree
1437 END IF
1438
1439 ! total charge constraint
1440 IF (resp_env%itc) THEN
1441 ncons_v = ncons_v + 1
1442 resp_env%matrix(1:natom, ncons_v) = 1.0_dp
1443 resp_env%matrix(ncons_v, 1:natom) = 1.0_dp
1444 resp_env%rhs(ncons_v) = dft_control%charge
1445 END IF
1446
1447 ! explicit constraints
1448 DO i = 1, resp_env%ncons_sec
1449 CALL build_atom_list(cons_section, subsys, atom_list_cons, i)
1450 IF (.NOT. resp_env%equal_charges) THEN
1451 ncons_v = ncons_v + 1
1452 CALL section_vals_val_get(cons_section, "ATOM_COEF", i_rep_section=i, r_vals=atom_coef)
1453 CALL section_vals_val_get(cons_section, "TARGET", i_rep_section=i, r_val=TARGET)
1454 cpassert(SIZE(atom_list_cons) == SIZE(atom_coef))
1455 DO m = 1, SIZE(atom_list_cons)
1456 resp_env%matrix(atom_list_cons(m), ncons_v) = atom_coef(m)
1457 resp_env%matrix(ncons_v, atom_list_cons(m)) = atom_coef(m)
1458 END DO
1459 resp_env%rhs(ncons_v) = TARGET
1460 ELSE
1461 my_atom_coef(1) = 1.0_dp
1462 my_atom_coef(2) = -1.0_dp
1463 DO k = 2, SIZE(atom_list_cons)
1464 ncons_v = ncons_v + 1
1465 resp_env%matrix(atom_list_cons(1), ncons_v) = my_atom_coef(1)
1466 resp_env%matrix(ncons_v, atom_list_cons(1)) = my_atom_coef(1)
1467 resp_env%matrix(atom_list_cons(k), ncons_v) = my_atom_coef(2)
1468 resp_env%matrix(ncons_v, atom_list_cons(k)) = my_atom_coef(2)
1469 resp_env%rhs(ncons_v) = 0.0_dp
1470 END DO
1471 END IF
1472 DEALLOCATE (atom_list_cons)
1473 END DO
1474 CALL timestop(handle)
1475
1476 END SUBROUTINE add_restraints_and_constraints
1477
1478! **************************************************************************************************
1479!> \brief print input information
1480!> \param qs_env the qs environment
1481!> \param resp_env the resp environment
1482!> \param rep_sys structure for repeating input sections defining fit points
1483!> \param my_per ...
1484! **************************************************************************************************
1485 SUBROUTINE print_resp_parameter_info(qs_env, resp_env, rep_sys, my_per)
1486
1487 TYPE(qs_environment_type), POINTER :: qs_env
1488 TYPE(resp_type), POINTER :: resp_env
1489 TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys
1490 INTEGER, INTENT(IN) :: my_per
1491
1492 CHARACTER(len=*), PARAMETER :: routinen = 'print_resp_parameter_info'
1493
1494 CHARACTER(len=2) :: symbol
1495 INTEGER :: handle, i, ikind, kind_number, nkinds, &
1496 output_unit
1497 REAL(kind=dp) :: conv, eta_conv
1498 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1499 TYPE(cp_logger_type), POINTER :: logger
1500 TYPE(section_vals_type), POINTER :: input, resp_section
1501
1502 CALL timeset(routinen, handle)
1503 NULLIFY (logger, input, resp_section)
1504
1505 CALL get_qs_env(qs_env, &
1506 input=input, &
1507 atomic_kind_set=atomic_kind_set)
1508 resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1509 logger => cp_get_default_logger()
1510 output_unit = cp_print_key_unit_nr(logger, resp_section, "PRINT%PROGRAM_RUN_INFO", &
1511 extension=".resp")
1512 nkinds = SIZE(atomic_kind_set)
1513
1514 conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
1515 IF (.NOT. my_per == use_perd_none) THEN
1516 eta_conv = cp_unit_from_cp2k(resp_env%eta, "angstrom", power=-2)
1517 END IF
1518
1519 IF (output_unit > 0) THEN
1520 WRITE (output_unit, '(/,1X,A,/)') "STARTING RESP FIT"
1521 IF (resp_env%use_repeat_method) THEN
1522 WRITE (output_unit, '(T3,A)') &
1523 "Fit the variance of the potential (REPEAT method)."
1524 END IF
1525 IF (.NOT. resp_env%equal_charges) THEN
1526 WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons_sec
1527 ELSE
1528 IF (resp_env%itc) THEN
1529 WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons - 1
1530 ELSE
1531 WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons
1532 END IF
1533 END IF
1534 WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit restraints: ", resp_env%nrest_sec
1535 WRITE (output_unit, '(T3,A,T80,A)') "Constrain total charge ", merge("T", "F", resp_env%itc)
1536 WRITE (output_unit, '(T3,A,T80,A)') "Restrain heavy atoms ", merge("T", "F", resp_env%rheavies)
1537 IF (resp_env%rheavies) THEN
1538 WRITE (output_unit, '(T3,A,T71,F10.6)') "Heavy atom restraint strength: ", &
1539 resp_env%rheavies_strength
1540 END IF
1541 WRITE (output_unit, '(T3,A,T66,3I5)') "Stride: ", resp_env%stride
1542 IF (resp_env%molecular_sys) THEN
1543 WRITE (output_unit, '(T3,A)') &
1544 "------------------------------------------------------------------------------"
1545 WRITE (output_unit, '(T3,A)') "Using sphere sampling"
1546 WRITE (output_unit, '(T3,A,T46,A,T66,A)') &
1547 "Element", "RMIN [angstrom]", "RMAX [angstrom]"
1548 DO ikind = 1, nkinds
1549 CALL get_atomic_kind(atomic_kind=atomic_kind_set(ikind), &
1550 kind_number=kind_number, &
1551 element_symbol=symbol)
1552 WRITE (output_unit, '(T3,A,T51,F10.5,T71,F10.5)') &
1553 symbol, &
1554 resp_env%rmin_kind(kind_number)*conv, &
1555 resp_env%rmax_kind(kind_number)*conv
1556 END DO
1557 IF (my_per == use_perd_none) THEN
1558 WRITE (output_unit, '(T3,A,T51,3F10.5)') "Box min [angstrom]: ", resp_env%box_low(1:3)*conv
1559 WRITE (output_unit, '(T3,A,T51,3F10.5)') "Box max [angstrom]: ", resp_env%box_hi(1:3)*conv
1560 END IF
1561 WRITE (output_unit, '(T3,A)') &
1562 "------------------------------------------------------------------------------"
1563 ELSE
1564 WRITE (output_unit, '(T3,A)') &
1565 "------------------------------------------------------------------------------"
1566 WRITE (output_unit, '(T3,A)') "Using slab sampling"
1567 WRITE (output_unit, '(2X,A,F10.5)') "Index of atoms defining the surface: "
1568 DO i = 1, SIZE(rep_sys)
1569 IF (i > 1 .AND. all(rep_sys(i)%p_resp%atom_surf_list == rep_sys(1)%p_resp%atom_surf_list)) EXIT
1570 WRITE (output_unit, '(7X,10I6)') rep_sys(i)%p_resp%atom_surf_list
1571 END DO
1572 DO i = 1, SIZE(rep_sys)
1573 IF (i > 1 .AND. all(rep_sys(i)%p_resp%range_surf == rep_sys(1)%p_resp%range_surf)) EXIT
1574 WRITE (output_unit, '(T3,A,T61,2F10.5)') &
1575 "Range for sampling above the surface [angstrom]:", &
1576 rep_sys(i)%p_resp%range_surf(1:2)*conv
1577 END DO
1578 DO i = 1, SIZE(rep_sys)
1579 IF (i > 1 .AND. rep_sys(i)%p_resp%length == rep_sys(1)%p_resp%length) EXIT
1580 WRITE (output_unit, '(T3,A,T71,F10.5)') "Length of sampling box above each"// &
1581 " surface atom [angstrom]: ", rep_sys(i)%p_resp%length*conv
1582 END DO
1583 WRITE (output_unit, '(T3,A)') &
1584 "------------------------------------------------------------------------------"
1585 END IF
1586 IF (.NOT. my_per == use_perd_none) THEN
1587 WRITE (output_unit, '(T3,A,T71,F10.5)') "Width of Gaussian charge"// &
1588 " distribution [angstrom^-2]: ", eta_conv
1589 END IF
1590 CALL m_flush(output_unit)
1591 END IF
1592 CALL cp_print_key_finished_output(output_unit, logger, resp_section, &
1593 "PRINT%PROGRAM_RUN_INFO")
1594
1595 CALL timestop(handle)
1596
1597 END SUBROUTINE print_resp_parameter_info
1598
1599! **************************************************************************************************
1600!> \brief print RESP charges to an extra file or to the normal output file
1601!> \param qs_env the qs environment
1602!> \param resp_env the resp environment
1603!> \param output_runinfo ...
1604!> \param natom number of atoms
1605! **************************************************************************************************
1606 SUBROUTINE print_resp_charges(qs_env, resp_env, output_runinfo, natom)
1607
1608 TYPE(qs_environment_type), POINTER :: qs_env
1609 TYPE(resp_type), POINTER :: resp_env
1610 INTEGER, INTENT(IN) :: output_runinfo, natom
1611
1612 CHARACTER(len=*), PARAMETER :: routinen = 'print_resp_charges'
1613
1614 CHARACTER(LEN=default_path_length) :: filename
1615 INTEGER :: handle, output_file
1616 TYPE(cp_logger_type), POINTER :: logger
1617 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1618 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1619 TYPE(section_vals_type), POINTER :: input, print_key, resp_section
1620
1621 CALL timeset(routinen, handle)
1622
1623 NULLIFY (particle_set, qs_kind_set, input, logger, resp_section, print_key)
1624
1625 CALL get_qs_env(qs_env, input=input, particle_set=particle_set, &
1626 qs_kind_set=qs_kind_set)
1627
1628 resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1629 print_key => section_vals_get_subs_vals(resp_section, &
1630 "PRINT%RESP_CHARGES_TO_FILE")
1631 logger => cp_get_default_logger()
1632
1633 IF (btest(cp_print_key_should_output(logger%iter_info, &
1634 resp_section, "PRINT%RESP_CHARGES_TO_FILE"), &
1635 cp_p_file)) THEN
1636 output_file = cp_print_key_unit_nr(logger, resp_section, &
1637 "PRINT%RESP_CHARGES_TO_FILE", &
1638 extension=".resp", &
1639 file_status="REPLACE", &
1640 file_action="WRITE", &
1641 file_form="FORMATTED")
1642 IF (output_file > 0) THEN
1643 filename = cp_print_key_generate_filename(logger, &
1644 print_key, extension=".resp", &
1645 my_local=.false.)
1646 CALL print_atomic_charges(particle_set, qs_kind_set, output_file, title="RESP charges:", &
1647 atomic_charges=resp_env%rhs(1:natom))
1648 IF (output_runinfo > 0) WRITE (output_runinfo, '(2X,A,/)') "PRINTED RESP CHARGES TO FILE"
1649 END IF
1650
1651 CALL cp_print_key_finished_output(output_file, logger, resp_section, &
1652 "PRINT%RESP_CHARGES_TO_FILE")
1653 ELSE
1654 CALL print_atomic_charges(particle_set, qs_kind_set, output_runinfo, title="RESP charges:", &
1655 atomic_charges=resp_env%rhs(1:natom))
1656 END IF
1657
1658 CALL timestop(handle)
1659
1660 END SUBROUTINE print_resp_charges
1661
1662! **************************************************************************************************
1663!> \brief print potential generated by RESP charges to file
1664!> \param qs_env the qs environment
1665!> \param resp_env the resp environment
1666!> \param particles ...
1667!> \param natom number of atoms
1668!> \param output_runinfo ...
1669! **************************************************************************************************
1670 SUBROUTINE print_pot_from_resp_charges(qs_env, resp_env, particles, natom, output_runinfo)
1671
1672 TYPE(qs_environment_type), POINTER :: qs_env
1673 TYPE(resp_type), POINTER :: resp_env
1674 TYPE(particle_list_type), POINTER :: particles
1675 INTEGER, INTENT(IN) :: natom, output_runinfo
1676
1677 CHARACTER(len=*), PARAMETER :: routinen = 'print_pot_from_resp_charges'
1678
1679 CHARACTER(LEN=default_path_length) :: my_pos_cube
1680 INTEGER :: handle, ip, jx, jy, jz, unit_nr
1681 LOGICAL :: append_cube, mpi_io
1682 REAL(kind=dp) :: dvol, normalize_factor, rms, rrms, &
1683 sum_diff, sum_hartree, udvol
1684 TYPE(cp_logger_type), POINTER :: logger
1685 TYPE(mp_para_env_type), POINTER :: para_env
1686 TYPE(pw_c1d_gs_type) :: rho_resp, v_resp_gspace
1687 TYPE(pw_env_type), POINTER :: pw_env
1688 TYPE(pw_poisson_type), POINTER :: poisson_env
1689 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1690 TYPE(pw_r3d_rs_type) :: aux_r, v_resp_rspace
1691 TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
1692 TYPE(section_vals_type), POINTER :: input, print_key, resp_section
1693
1694 CALL timeset(routinen, handle)
1695
1696 NULLIFY (auxbas_pw_pool, logger, pw_env, poisson_env, input, print_key, &
1697 para_env, resp_section, v_hartree_rspace)
1698 CALL get_qs_env(qs_env, &
1699 input=input, &
1700 para_env=para_env, &
1701 pw_env=pw_env, &
1702 v_hartree_rspace=v_hartree_rspace)
1703 normalize_factor = sqrt((resp_env%eta/pi)**3)
1704 resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1705 print_key => section_vals_get_subs_vals(resp_section, &
1706 "PRINT%V_RESP_CUBE")
1707 logger => cp_get_default_logger()
1708
1709 !*** calculate potential generated from RESP charges
1710 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1711 poisson_env=poisson_env)
1712
1713 CALL auxbas_pw_pool%create_pw(rho_resp)
1714 CALL auxbas_pw_pool%create_pw(v_resp_gspace)
1715 CALL auxbas_pw_pool%create_pw(v_resp_rspace)
1716
1717 CALL pw_zero(rho_resp)
1718 CALL calculate_rho_resp_all(rho_resp, resp_env%rhs, natom, &
1719 resp_env%eta, qs_env)
1720 CALL pw_zero(v_resp_gspace)
1721 CALL pw_poisson_solve(poisson_env, rho_resp, &
1722 vhartree=v_resp_gspace)
1723 CALL pw_zero(v_resp_rspace)
1724 CALL pw_transfer(v_resp_gspace, v_resp_rspace)
1725 dvol = v_resp_rspace%pw_grid%dvol
1726 CALL pw_scale(v_resp_rspace, dvol)
1727 CALL pw_scale(v_resp_rspace, -normalize_factor)
1728 ! REPEAT: correct for offset, take into account that potentials have reverse sign
1729 ! and are scaled by dvol
1730 IF (resp_env%use_repeat_method) THEN
1731 v_resp_rspace%array(:, :, :) = v_resp_rspace%array(:, :, :) - resp_env%offset*dvol
1732 END IF
1733 CALL v_resp_gspace%release()
1734 CALL rho_resp%release()
1735
1736 !***now print the v_resp_rspace%pw to a cube file if requested
1737 IF (btest(cp_print_key_should_output(logger%iter_info, resp_section, &
1738 "PRINT%V_RESP_CUBE"), cp_p_file)) THEN
1739 CALL auxbas_pw_pool%create_pw(aux_r)
1740 append_cube = section_get_lval(resp_section, "PRINT%V_RESP_CUBE%APPEND")
1741 my_pos_cube = "REWIND"
1742 IF (append_cube) THEN
1743 my_pos_cube = "APPEND"
1744 END IF
1745 mpi_io = .true.
1746 unit_nr = cp_print_key_unit_nr(logger, resp_section, &
1747 "PRINT%V_RESP_CUBE", &
1748 extension=".cube", &
1749 file_position=my_pos_cube, &
1750 mpi_io=mpi_io)
1751 udvol = 1.0_dp/dvol
1752 CALL pw_copy(v_resp_rspace, aux_r)
1753 CALL pw_scale(aux_r, udvol)
1754 CALL cp_pw_to_cube(aux_r, unit_nr, "RESP POTENTIAL", particles=particles, &
1755 stride=section_get_ivals(resp_section, &
1756 "PRINT%V_RESP_CUBE%STRIDE"), &
1757 mpi_io=mpi_io)
1758 CALL cp_print_key_finished_output(unit_nr, logger, resp_section, &
1759 "PRINT%V_RESP_CUBE", mpi_io=mpi_io)
1760 CALL auxbas_pw_pool%give_back_pw(aux_r)
1761 END IF
1762
1763 !*** RMS and RRMS
1764 sum_diff = 0.0_dp
1765 sum_hartree = 0.0_dp
1766 rms = 0.0_dp
1767 rrms = 0.0_dp
1768 DO ip = 1, resp_env%npoints_proc
1769 jx = resp_env%fitpoints(1, ip)
1770 jy = resp_env%fitpoints(2, ip)
1771 jz = resp_env%fitpoints(3, ip)
1772 sum_diff = sum_diff + (v_hartree_rspace%array(jx, jy, jz) - &
1773 v_resp_rspace%array(jx, jy, jz))**2
1774 sum_hartree = sum_hartree + v_hartree_rspace%array(jx, jy, jz)**2
1775 END DO
1776 CALL para_env%sum(sum_diff)
1777 CALL para_env%sum(sum_hartree)
1778 rms = sqrt(sum_diff/resp_env%npoints)
1779 rrms = sqrt(sum_diff/sum_hartree)
1780 IF (output_runinfo > 0) THEN
1781 WRITE (output_runinfo, '(2X,A,T69,ES12.5)') "Root-mean-square (RMS) "// &
1782 "error of RESP fit:", rms
1783 WRITE (output_runinfo, '(2X,A,T69,ES12.5,/)') "Relative root-mean-square "// &
1784 "(RRMS) error of RESP fit:", rrms
1785 END IF
1786
1787 CALL v_resp_rspace%release()
1788
1789 CALL timestop(handle)
1790
1791 END SUBROUTINE print_pot_from_resp_charges
1792
1793END MODULE qs_resp
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
cell_types::pbc
Definition cell_types.F:104
kahan_sum::accurate_sum
Definition kahan_sum.F:41
memory_utilities::reallocate
Definition memory_utilities.F:27
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
qs_collocate_density::calculate_rho_resp_all
Definition qs_collocate_density.F:143
atomic_charges
simple routine to print charges for all atomic charge methods (currently mulliken,...
Definition atomic_charges.F:14
atomic_charges::print_atomic_charges
subroutine, public print_atomic_charges(particle_set, qs_kind_set, scr, title, electronic_charges, atomic_charges)
generates a unified output format for atomic charges
Definition atomic_charges.F:48
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::golze2015
integer, save, public golze2015
Definition bibliography.F:36
bibliography::rappe1992
integer, save, public rappe1992
Definition bibliography.F:36
bibliography::campana2009
integer, save, public campana2009
Definition bibliography.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::use_perd_xyz
integer, parameter, public use_perd_xyz
Definition cell_types.F:42
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:233
cell_types::use_perd_none
integer, parameter, public use_perd_none
Definition cell_types.F:42
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_generate_filename
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
Definition cp_output_handling.F:754
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_realspace_grid_cube
A wrapper around pw_to_cube() which accepts particle_list_type.
Definition cp_realspace_grid_cube.F:12
cp_realspace_grid_cube::cp_pw_to_cube
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, zeff, stride, max_file_size_mb, zero_tails, silent, mpi_io)
...
Definition cp_realspace_grid_cube.F:48
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1251
cp_units::cp_unit_to_cp2k
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition cp_units.F:1222
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::do_resp_minus_x_dir
integer, parameter, public do_resp_minus_x_dir
Definition input_constants.F:1154
input_constants::do_resp_minus_y_dir
integer, parameter, public do_resp_minus_y_dir
Definition input_constants.F:1154
input_constants::do_resp_z_dir
integer, parameter, public do_resp_z_dir
Definition input_constants.F:1154
input_constants::do_resp_minus_z_dir
integer, parameter, public do_resp_minus_z_dir
Definition input_constants.F:1154
input_constants::use_uff_vdw_radii
integer, parameter, public use_uff_vdw_radii
Definition input_constants.F:1154
input_constants::do_resp_y_dir
integer, parameter, public do_resp_y_dir
Definition input_constants.F:1154
input_constants::do_resp_x_dir
integer, parameter, public do_resp_x_dir
Definition input_constants.F:1154
input_constants::use_cambridge_vdw_radii
integer, parameter, public use_cambridge_vdw_radii
Definition input_constants.F:1154
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_get_ivals
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
Definition input_section_types.F:1001
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition input_section_types.F:1017
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::get_ptable_info
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition periodic_table.F:85
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_poisson_methods
Definition pw_poisson_methods.F:13
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::calculate_rho_resp_single
subroutine, public calculate_rho_resp_single(rho_gb, qs_env, eta, iatom_in)
collocate a single Gaussian on the grid for periodic RESP fitting
Definition qs_collocate_density.F:1408
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, mimic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1120
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_resp
provides a resp fit for gas phase systems
Definition qs_resp.F:16
qs_resp::resp_fit
subroutine, public resp_fit(qs_env)
performs resp fit and generates RESP charges
Definition qs_resp.F:132
qs_subsys_types
types that represent a quickstep subsys
Definition qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition qs_subsys_types.F:137
uff_vdw_radii_table
provides a table for UFF vdW radii: Rappe et al. J. Am. Chem. Soc. 114, 10024 (1992)
Definition uff_vdw_radii_table.F:15
uff_vdw_radii_table::get_uff_vdw_radius
pure subroutine, public get_uff_vdw_radius(z, radius, found)
get UFF vdW radius for a given element
Definition uff_vdw_radii_table.F:163
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
message_passing::mp_request_type
Definition message_passing.F:610
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:123
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56