(git:b279b6b)
qs_resp.F
Go to the documentation of this file.
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 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 ! **************************************************************************************************
16 MODULE 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, itc, &
100  molecular_sys, rheavies, &
101  use_repeat_method
102  INTEGER :: nres, ncons, &
103  nrest_sec, ncons_sec, &
104  npoints, stride(3), my_fit, &
105  npoints_proc, &
106  auto_vdw_radii_table
107  INTEGER, DIMENSION(:), POINTER :: atom_surf_list
108  INTEGER, DIMENSION(:, :), POINTER :: fitpoints
109  REAL(KIND=dp) :: rheavies_strength, &
110  length, eta, &
111  sum_vhartree, offset
112  REAL(KIND=dp), DIMENSION(3) :: box_hi, box_low
113  REAL(KIND=dp), DIMENSION(:), POINTER :: rmin_kind, &
114  rmax_kind
115  REAL(KIND=dp), DIMENSION(:), POINTER :: range_surf
116  REAL(KIND=dp), DIMENSION(:), POINTER :: rhs
117  REAL(KIND=dp), DIMENSION(:), POINTER :: sum_vpot
118  REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix
119  END TYPE resp_type
120 
121  TYPE resp_p_type
122  TYPE(resp_type), POINTER :: p_resp
123  END TYPE resp_p_type
124 
125 CONTAINS
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) &
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  IF (SIZE(my_stride) == 1) THEN
403  DO i = 1, 3
404  resp_env%stride(i) = my_stride(1)
405  END DO
406  ELSE
407  resp_env%stride = my_stride(1:3)
408  END IF
409  CALL section_vals_val_get(resp_section, "WIDTH", r_val=resp_env%eta)
410 
411  ! get if the user wants to use REPEAT method
412  CALL section_vals_val_get(resp_section, "USE_REPEAT_METHOD", &
413  l_val=resp_env%use_repeat_method)
414  IF (resp_env%use_repeat_method) THEN
415  resp_env%ncons = resp_env%ncons + 1
416  ! restrain heavies should be off
417  resp_env%rheavies = .false.
418  END IF
419 
420  ! get and set the parameters for molecular (non-surface) systems
421  ! this must come after the repeat settings being set
422  CALL get_parameter_molecular_sys(resp_env, sphere_section, cell, &
423  atomic_kind_set)
424 
425  ! get the parameter for periodic/surface systems
426  CALL section_vals_get(slab_section, explicit=explicit, n_repetition=nrep)
427  IF (explicit) THEN
428  IF (resp_env%molecular_sys) THEN
429  CALL cp_abort(__location__, &
430  "You can only use either SPHERE_SAMPLING or SLAB_SAMPLING, but "// &
431  "not both.")
432  END IF
433  ALLOCATE (rep_sys(nrep))
434  DO i = 1, nrep
435  ALLOCATE (rep_sys(i)%p_resp)
436  NULLIFY (rep_sys(i)%p_resp%range_surf, rep_sys(i)%p_resp%atom_surf_list)
437  CALL section_vals_val_get(slab_section, "RANGE", r_vals=rep_sys(i)%p_resp%range_surf, &
438  i_rep_section=i)
439  CALL section_vals_val_get(slab_section, "LENGTH", r_val=rep_sys(i)%p_resp%length, &
440  i_rep_section=i)
441  CALL section_vals_val_get(slab_section, "SURF_DIRECTION", &
442  i_rep_section=i, i_val=rep_sys(i)%p_resp%my_fit)
443  IF (any(rep_sys(i)%p_resp%range_surf < 0.0_dp)) THEN
444  cpabort("Numbers in RANGE in SLAB_SAMPLING cannot be negative.")
445  END IF
446  IF (rep_sys(i)%p_resp%length <= epsilon(0.0_dp)) THEN
447  cpabort("Parameter LENGTH in SLAB_SAMPLING has to be larger than zero.")
448  END IF
449  !list of atoms specifying the surface
450  CALL build_atom_list(slab_section, subsys, rep_sys(i)%p_resp%atom_surf_list, rep=i)
451  END DO
452  END IF
453 
454  ! get the parameters for the constraint and restraint sections
455  CALL section_vals_get(cons_section, explicit=explicit)
456  IF (explicit) THEN
457  CALL section_vals_get(cons_section, n_repetition=resp_env%ncons_sec)
458  DO i = 1, resp_env%ncons_sec
459  CALL section_vals_val_get(cons_section, "EQUAL_CHARGES", &
460  l_val=resp_env%equal_charges, explicit=explicit)
461  IF (.NOT. explicit) cycle
462  CALL build_atom_list(cons_section, subsys, atom_list_cons, i)
463  !instead of using EQUAL_CHARGES the constraint sections could be repeated
464  resp_env%ncons = resp_env%ncons + SIZE(atom_list_cons) - 2
465  DEALLOCATE (atom_list_cons)
466  END DO
467  END IF
468  CALL section_vals_get(rest_section, explicit=explicit)
469  IF (explicit) THEN
470  CALL section_vals_get(rest_section, n_repetition=resp_env%nrest_sec)
471  END IF
472  resp_env%ncons = resp_env%ncons + resp_env%ncons_sec
473  resp_env%nres = resp_env%nres + resp_env%nrest_sec
474 
475  CALL timestop(handle)
476 
477  END SUBROUTINE init_resp
478 
479 ! **************************************************************************************************
480 !> \brief getting the parameters for nonperiodic/non-surface systems
481 !> \param resp_env the resp environment
482 !> \param sphere_section input section setting parameters for sampling
483 !> fitting in spheres around the atom
484 !> \param cell parameters related to the simulation cell
485 !> \param atomic_kind_set ...
486 ! **************************************************************************************************
487  SUBROUTINE get_parameter_molecular_sys(resp_env, sphere_section, cell, &
488  atomic_kind_set)
489 
490  TYPE(resp_type), POINTER :: resp_env
491  TYPE(section_vals_type), POINTER :: sphere_section
492  TYPE(cell_type), POINTER :: cell
493  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
494 
495  CHARACTER(LEN=2) :: symbol
496  CHARACTER(LEN=default_string_length) :: missing_rmax, missing_rmin
497  CHARACTER(LEN=default_string_length), &
498  DIMENSION(:), POINTER :: tmpstringlist
499  INTEGER :: ikind, j, kind_number, n_rmax_missing, &
500  n_rmin_missing, nkind, nrep_rmax, &
501  nrep_rmin, z
502  LOGICAL :: explicit, has_rmax, has_rmin
503  LOGICAL, ALLOCATABLE, DIMENSION(:) :: rmax_is_set, rmin_is_set
504  REAL(kind=dp) :: auto_rmax_scale, auto_rmin_scale, rmax, &
505  rmin
506  REAL(kind=dp), DIMENSION(3, 3) :: hmat
507  TYPE(atomic_kind_type), POINTER :: atomic_kind
508 
509  nrep_rmin = 0
510  nrep_rmax = 0
511  nkind = SIZE(atomic_kind_set)
512 
513  has_rmin = .false.
514  has_rmax = .false.
515 
516  CALL section_vals_get(sphere_section, explicit=explicit)
517  IF (explicit) THEN
518  resp_env%molecular_sys = .true.
519  CALL section_vals_val_get(sphere_section, "AUTO_VDW_RADII_TABLE", &
520  i_val=resp_env%auto_vdw_radii_table)
521  CALL section_vals_val_get(sphere_section, "AUTO_RMIN_SCALE", r_val=auto_rmin_scale)
522  CALL section_vals_val_get(sphere_section, "AUTO_RMAX_SCALE", r_val=auto_rmax_scale)
523  CALL section_vals_val_get(sphere_section, "RMIN", explicit=has_rmin, r_val=rmin)
524  CALL section_vals_val_get(sphere_section, "RMAX", explicit=has_rmax, r_val=rmax)
525  CALL section_vals_val_get(sphere_section, "RMIN_KIND", n_rep_val=nrep_rmin)
526  CALL section_vals_val_get(sphere_section, "RMAX_KIND", n_rep_val=nrep_rmax)
527  ALLOCATE (resp_env%rmin_kind(nkind))
528  ALLOCATE (resp_env%rmax_kind(nkind))
529  resp_env%rmin_kind = 0.0_dp
530  resp_env%rmax_kind = 0.0_dp
531  ALLOCATE (rmin_is_set(nkind))
532  ALLOCATE (rmax_is_set(nkind))
533  rmin_is_set = .false.
534  rmax_is_set = .false.
535  ! define rmin_kind and rmax_kind to predefined vdW radii
536  DO ikind = 1, nkind
537  atomic_kind => atomic_kind_set(ikind)
538  CALL get_atomic_kind(atomic_kind, &
539  element_symbol=symbol, &
540  kind_number=kind_number, &
541  z=z)
542  SELECT CASE (resp_env%auto_vdw_radii_table)
543  CASE (use_cambridge_vdw_radii)
544  CALL get_ptable_info(symbol, vdw_radius=resp_env%rmin_kind(kind_number))
545  rmin_is_set(kind_number) = .true.
546  CASE (use_uff_vdw_radii)
547  CALL cite_reference(rappe1992)
548  CALL get_uff_vdw_radius(z, radius=resp_env%rmin_kind(kind_number), &
549  found=rmin_is_set(kind_number))
550  CASE DEFAULT
551  CALL get_ptable_info(symbol, vdw_radius=resp_env%rmin_kind(kind_number))
552  rmin_is_set(kind_number) = .true.
553  END SELECT
554  IF (rmin_is_set(kind_number)) THEN
555  resp_env%rmin_kind(kind_number) = cp_unit_to_cp2k(resp_env%rmin_kind(kind_number), &
556  "angstrom")
557  resp_env%rmin_kind(kind_number) = resp_env%rmin_kind(kind_number)*auto_rmin_scale
558  ! set RMAX_KIND accourding by scaling RMIN_KIND
559  resp_env%rmax_kind(kind_number) = &
560  max(resp_env%rmin_kind(kind_number), &
561  resp_env%rmin_kind(kind_number)*auto_rmax_scale)
562  rmax_is_set(kind_number) = .true.
563  END IF
564  END DO
565  ! if RMIN or RMAX are present, overwrite the rmin_kind(:) and
566  ! rmax_kind(:) to those values
567  IF (has_rmin) THEN
568  resp_env%rmin_kind = rmin
569  rmin_is_set = .true.
570  END IF
571  IF (has_rmax) THEN
572  resp_env%rmax_kind = rmax
573  rmax_is_set = .true.
574  END IF
575  ! if RMIN_KIND's or RMAX_KIND's are present, overwrite the
576  ! rmin_kinds(:) or rmax_kind(:) to those values
577  DO j = 1, nrep_rmin
578  CALL section_vals_val_get(sphere_section, "RMIN_KIND", i_rep_val=j, &
579  c_vals=tmpstringlist)
580  DO ikind = 1, nkind
581  atomic_kind => atomic_kind_set(ikind)
582  CALL get_atomic_kind(atomic_kind, element_symbol=symbol, kind_number=kind_number)
583  IF (trim(tmpstringlist(2)) == trim(symbol)) THEN
584  READ (tmpstringlist(1), *) resp_env%rmin_kind(kind_number)
585  resp_env%rmin_kind(kind_number) = &
586  cp_unit_to_cp2k(resp_env%rmin_kind(kind_number), &
587  "angstrom")
588  rmin_is_set(kind_number) = .true.
589  END IF
590  END DO
591  END DO
592  DO j = 1, nrep_rmax
593  CALL section_vals_val_get(sphere_section, "RMAX_KIND", i_rep_val=j, &
594  c_vals=tmpstringlist)
595  DO ikind = 1, nkind
596  atomic_kind => atomic_kind_set(ikind)
597  CALL get_atomic_kind(atomic_kind, element_symbol=symbol, kind_number=kind_number)
598  IF (trim(tmpstringlist(2)) == trim(symbol)) THEN
599  READ (tmpstringlist(1), *) resp_env%rmax_kind(kind_number)
600  resp_env%rmax_kind(kind_number) = cp_unit_to_cp2k(resp_env%rmax_kind(kind_number), &
601  "angstrom")
602  rmax_is_set(kind_number) = .true.
603  END IF
604  END DO
605  END DO
606  ! check if rmin and rmax are set for each kind
607  n_rmin_missing = 0
608  n_rmax_missing = 0
609  missing_rmin = ""
610  missing_rmax = ""
611  DO ikind = 1, nkind
612  atomic_kind => atomic_kind_set(ikind)
613  CALL get_atomic_kind(atomic_kind, &
614  element_symbol=symbol, &
615  kind_number=kind_number)
616  IF (.NOT. rmin_is_set(kind_number)) THEN
617  n_rmin_missing = n_rmin_missing + 1
618  missing_rmin = trim(missing_rmin)//" "//trim(symbol)//","
619  END IF
620  IF (.NOT. rmax_is_set(kind_number)) THEN
621  n_rmax_missing = n_rmax_missing + 1
622  missing_rmax = trim(missing_rmax)//" "//trim(symbol)//","
623  END IF
624  END DO
625  IF (n_rmin_missing > 0) THEN
626  CALL cp_warn(__location__, &
627  "RMIN for the following elements are missing: "// &
628  trim(missing_rmin)// &
629  " please set these values manually using "// &
630  "RMIN_KIND in SPHERE_SAMPLING section")
631  END IF
632  IF (n_rmax_missing > 0) THEN
633  CALL cp_warn(__location__, &
634  "RMAX for the following elements are missing: "// &
635  trim(missing_rmax)// &
636  " please set these values manually using "// &
637  "RMAX_KIND in SPHERE_SAMPLING section")
638  END IF
639  IF (n_rmin_missing > 0 .OR. &
640  n_rmax_missing > 0) THEN
641  cpabort("Insufficient data for RMIN or RMAX")
642  END IF
643 
644  CALL get_cell(cell=cell, h=hmat)
645  resp_env%box_hi = (/hmat(1, 1), hmat(2, 2), hmat(3, 3)/)
646  resp_env%box_low = 0.0_dp
647  CALL section_vals_val_get(sphere_section, "X_HI", explicit=explicit)
648  IF (explicit) CALL section_vals_val_get(sphere_section, "X_HI", &
649  r_val=resp_env%box_hi(1))
650  CALL section_vals_val_get(sphere_section, "X_LOW", explicit=explicit)
651  IF (explicit) CALL section_vals_val_get(sphere_section, "X_LOW", &
652  r_val=resp_env%box_low(1))
653  CALL section_vals_val_get(sphere_section, "Y_HI", explicit=explicit)
654  IF (explicit) CALL section_vals_val_get(sphere_section, "Y_HI", &
655  r_val=resp_env%box_hi(2))
656  CALL section_vals_val_get(sphere_section, "Y_LOW", explicit=explicit)
657  IF (explicit) CALL section_vals_val_get(sphere_section, "Y_LOW", &
658  r_val=resp_env%box_low(2))
659  CALL section_vals_val_get(sphere_section, "Z_HI", explicit=explicit)
660  IF (explicit) CALL section_vals_val_get(sphere_section, "Z_HI", &
661  r_val=resp_env%box_hi(3))
662  CALL section_vals_val_get(sphere_section, "Z_LOW", explicit=explicit)
663  IF (explicit) CALL section_vals_val_get(sphere_section, "Z_LOW", &
664  r_val=resp_env%box_low(3))
665 
666  DEALLOCATE (rmin_is_set)
667  DEALLOCATE (rmax_is_set)
668  END IF
669 
670  END SUBROUTINE get_parameter_molecular_sys
671 
672 ! **************************************************************************************************
673 !> \brief building atom lists for different sections of RESP
674 !> \param section input section
675 !> \param subsys ...
676 !> \param atom_list list of atoms for restraints, constraints and fit point
677 !> sampling for slab-like systems
678 !> \param rep input section can be repeated, this param defines for which
679 !> repetition of the input section the atom_list is built
680 ! **************************************************************************************************
681  SUBROUTINE build_atom_list(section, subsys, atom_list, rep)
682 
683  TYPE(section_vals_type), POINTER :: section
684  TYPE(qs_subsys_type), POINTER :: subsys
685  INTEGER, DIMENSION(:), POINTER :: atom_list
686  INTEGER, INTENT(IN), OPTIONAL :: rep
687 
688  CHARACTER(len=*), PARAMETER :: routinen = 'build_atom_list'
689 
690  INTEGER :: atom_a, atom_b, handle, i, irep, j, &
691  max_index, n_var, num_atom
692  INTEGER, DIMENSION(:), POINTER :: indexes
693  LOGICAL :: index_in_range
694 
695  CALL timeset(routinen, handle)
696 
697  NULLIFY (indexes)
698  irep = 1
699  IF (PRESENT(rep)) irep = rep
700 
701  CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
702  n_rep_val=n_var)
703  num_atom = 0
704  DO i = 1, n_var
705  CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
706  i_rep_val=i, i_vals=indexes)
707  num_atom = num_atom + SIZE(indexes)
708  END DO
709  ALLOCATE (atom_list(num_atom))
710  atom_list = 0
711  num_atom = 1
712  DO i = 1, n_var
713  CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
714  i_rep_val=i, i_vals=indexes)
715  atom_list(num_atom:num_atom + SIZE(indexes) - 1) = indexes(:)
716  num_atom = num_atom + SIZE(indexes)
717  END DO
718  !check atom list
719  num_atom = num_atom - 1
720  CALL qs_subsys_get(subsys, nparticle=max_index)
721  cpassert(SIZE(atom_list) /= 0)
722  index_in_range = (maxval(atom_list) <= max_index) &
723  .AND. (minval(atom_list) > 0)
724  cpassert(index_in_range)
725  DO i = 1, num_atom
726  DO j = i + 1, num_atom
727  atom_a = atom_list(i)
728  atom_b = atom_list(j)
729  IF (atom_a == atom_b) &
730  cpabort("There are atoms doubled in atom list for RESP.")
731  END DO
732  END DO
733 
734  CALL timestop(handle)
735 
736  END SUBROUTINE build_atom_list
737 
738 ! **************************************************************************************************
739 !> \brief build matrix and vector for nonperiodic RESP fitting
740 !> \param qs_env the qs environment
741 !> \param resp_env the resp environment
742 !> \param atomic_kind_set ...
743 !> \param particles ...
744 !> \param cell parameters related to the simulation cell
745 !> \param matrix coefficient matrix of the linear system of equations
746 !> \param rhs vector of the linear system of equations
747 !> \param natom number of atoms
748 ! **************************************************************************************************
749  SUBROUTINE calc_resp_matrix_nonper(qs_env, resp_env, atomic_kind_set, particles, &
750  cell, matrix, rhs, natom)
751 
752  TYPE(qs_environment_type), POINTER :: qs_env
753  TYPE(resp_type), POINTER :: resp_env
754  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
755  TYPE(particle_list_type), POINTER :: particles
756  TYPE(cell_type), POINTER :: cell
757  REAL(kind=dp), DIMENSION(:, :), POINTER :: matrix
758  REAL(kind=dp), DIMENSION(:), POINTER :: rhs
759  INTEGER, INTENT(IN) :: natom
760 
761  CHARACTER(len=*), PARAMETER :: routinen = 'calc_resp_matrix_nonper'
762 
763  INTEGER :: bo(2, 3), gbo(2, 3), handle, i, ikind, &
764  jx, jy, jz, k, kind_number, l, m, &
765  nkind, now, np(3), p
766  LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: not_in_range
767  REAL(kind=dp) :: delta, dh(3, 3), dvol, r(3), rmax, rmin, &
768  vec(3), vec_pbc(3), vj
769  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dist
770  REAL(kind=dp), DIMENSION(3, 3) :: hmat, hmat_inv
771  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
772  TYPE(pw_r3d_rs_type), POINTER :: v_hartree_pw
773 
774  CALL timeset(routinen, handle)
775 
776  NULLIFY (particle_set, v_hartree_pw)
777  delta = 1.0e-13_dp
778 
779  CALL get_cell(cell=cell, h=hmat, h_inv=hmat_inv)
780 
781  IF (.NOT. cell%orthorhombic) THEN
782  CALL cp_abort(__location__, &
783  "Nonperiodic solution for RESP charges only"// &
784  " implemented for orthorhombic cells!")
785  END IF
786  IF (.NOT. resp_env%molecular_sys) THEN
787  CALL cp_abort(__location__, &
788  "Nonperiodic solution for RESP charges (i.e. nonperiodic"// &
789  " Poisson solver) can only be used with section SPHERE_SAMPLING")
790  END IF
791  IF (resp_env%use_repeat_method) THEN
792  CALL cp_abort(__location__, &
793  "REPEAT method only reasonable for periodic RESP fitting")
794  END IF
795  CALL get_qs_env(qs_env, particle_set=particle_set, v_hartree_rspace=v_hartree_pw)
796 
797  bo = v_hartree_pw%pw_grid%bounds_local
798  gbo = v_hartree_pw%pw_grid%bounds
799  np = v_hartree_pw%pw_grid%npts
800  dh = v_hartree_pw%pw_grid%dh
801  dvol = v_hartree_pw%pw_grid%dvol
802  nkind = SIZE(atomic_kind_set)
803 
804  ALLOCATE (dist(natom))
805  ALLOCATE (not_in_range(natom, 2))
806 
807  ! store fitting points to calculate the RMS and RRMS later
808  IF (.NOT. ASSOCIATED(resp_env%fitpoints)) THEN
809  now = 1000
810  ALLOCATE (resp_env%fitpoints(3, now))
811  ELSE
812  now = SIZE(resp_env%fitpoints, 2)
813  END IF
814 
815  DO jz = bo(1, 3), bo(2, 3)
816  DO jy = bo(1, 2), bo(2, 2)
817  DO jx = bo(1, 1), bo(2, 1)
818  IF (.NOT. (modulo(jz, resp_env%stride(3)) == 0)) cycle
819  IF (.NOT. (modulo(jy, resp_env%stride(2)) == 0)) cycle
820  IF (.NOT. (modulo(jx, resp_env%stride(1)) == 0)) cycle
821  !bounds bo reach from -np/2 to np/2. shift of np/2 so that r(1,1,1)=(0,0,0)
822  l = jx - gbo(1, 1)
823  k = jy - gbo(1, 2)
824  p = jz - gbo(1, 3)
825  r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
826  r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
827  r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
828  IF (r(3) < resp_env%box_low(3) .OR. r(3) > resp_env%box_hi(3)) cycle
829  IF (r(2) < resp_env%box_low(2) .OR. r(2) > resp_env%box_hi(2)) cycle
830  IF (r(1) < resp_env%box_low(1) .OR. r(1) > resp_env%box_hi(1)) cycle
831  ! compute distance from the grid point to all atoms
832  not_in_range = .false.
833  DO i = 1, natom
834  vec = r - particles%els(i)%r
835  vec_pbc(1) = vec(1) - hmat(1, 1)*anint(hmat_inv(1, 1)*vec(1))
836  vec_pbc(2) = vec(2) - hmat(2, 2)*anint(hmat_inv(2, 2)*vec(2))
837  vec_pbc(3) = vec(3) - hmat(3, 3)*anint(hmat_inv(3, 3)*vec(3))
838  dist(i) = sqrt(sum(vec_pbc**2))
839  CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, &
840  kind_number=kind_number)
841  DO ikind = 1, nkind
842  IF (ikind == kind_number) THEN
843  rmin = resp_env%rmin_kind(ikind)
844  rmax = resp_env%rmax_kind(ikind)
845  EXIT
846  END IF
847  END DO
848  IF (dist(i) < rmin + delta) not_in_range(i, 1) = .true.
849  IF (dist(i) > rmax - delta) not_in_range(i, 2) = .true.
850  END DO
851  ! check if the point is sufficiently close and far. if OK, we can use
852  ! the point for fitting, add/subtract 1.0E-13 to get rid of rounding errors when shifting atoms
853  IF (any(not_in_range(:, 1)) .OR. all(not_in_range(:, 2))) cycle
854  resp_env%npoints_proc = resp_env%npoints_proc + 1
855  IF (resp_env%npoints_proc > now) THEN
856  now = 2*now
857  CALL reallocate(resp_env%fitpoints, 1, 3, 1, now)
858  END IF
859  resp_env%fitpoints(1, resp_env%npoints_proc) = jx
860  resp_env%fitpoints(2, resp_env%npoints_proc) = jy
861  resp_env%fitpoints(3, resp_env%npoints_proc) = jz
862  ! correct for the fact that v_hartree is scaled by dvol, and has the opposite sign
863  IF (qs_env%qmmm) THEN
864  ! If it's a QM/MM run let's remove the contribution of the MM potential out of the Hartree pot
865  vj = -v_hartree_pw%array(jx, jy, jz)/dvol + qs_env%ks_qmmm_env%v_qmmm_rspace%array(jx, jy, jz)
866  ELSE
867  vj = -v_hartree_pw%array(jx, jy, jz)/dvol
868  END IF
869  dist(:) = 1.0_dp/dist(:)
870 
871  DO i = 1, natom
872  DO m = 1, natom
873  matrix(m, i) = matrix(m, i) + 2.0_dp*dist(i)*dist(m)
874  END DO
875  rhs(i) = rhs(i) + 2.0_dp*vj*dist(i)
876  END DO
877  END DO
878  END DO
879  END DO
880 
881  resp_env%npoints = resp_env%npoints_proc
882  CALL v_hartree_pw%pw_grid%para%group%sum(resp_env%npoints)
883  CALL v_hartree_pw%pw_grid%para%group%sum(matrix)
884  CALL v_hartree_pw%pw_grid%para%group%sum(rhs)
885  !weighted sum
886  matrix = matrix/resp_env%npoints
887  rhs = rhs/resp_env%npoints
888 
889  DEALLOCATE (dist)
890  DEALLOCATE (not_in_range)
891 
892  CALL timestop(handle)
893 
894  END SUBROUTINE calc_resp_matrix_nonper
895 
896 ! **************************************************************************************************
897 !> \brief build matrix and vector for periodic RESP fitting
898 !> \param qs_env the qs environment
899 !> \param resp_env the resp environment
900 !> \param rep_sys structure for repeating input sections defining fit points
901 !> \param particles ...
902 !> \param cell parameters related to the simulation cell
903 !> \param natom number of atoms
904 ! **************************************************************************************************
905  SUBROUTINE calc_resp_matrix_periodic(qs_env, resp_env, rep_sys, particles, cell, &
906  natom)
907 
908  TYPE(qs_environment_type), POINTER :: qs_env
909  TYPE(resp_type), POINTER :: resp_env
910  TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys
911  TYPE(particle_list_type), POINTER :: particles
912  TYPE(cell_type), POINTER :: cell
913  INTEGER, INTENT(IN) :: natom
914 
915  CHARACTER(len=*), PARAMETER :: routinen = 'calc_resp_matrix_periodic'
916 
917  INTEGER :: handle, i, ip, j, jx, jy, jz
918  INTEGER, DIMENSION(3) :: periodic
919  REAL(kind=dp) :: normalize_factor
920  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: vpot
921  TYPE(mp_para_env_type), POINTER :: para_env
922  TYPE(pw_c1d_gs_type) :: rho_ga, va_gspace
923  TYPE(pw_env_type), POINTER :: pw_env
924  TYPE(pw_poisson_type), POINTER :: poisson_env
925  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
926  TYPE(pw_r3d_rs_type) :: va_rspace
927 
928  CALL timeset(routinen, handle)
929 
930  NULLIFY (pw_env, para_env, auxbas_pw_pool, poisson_env)
931 
932  CALL get_cell(cell=cell, periodic=periodic)
933 
934  IF (.NOT. all(periodic /= 0)) THEN
935  CALL cp_abort(__location__, &
936  "Periodic solution for RESP (with periodic Poisson solver)"// &
937  " can only be obtained with a cell that has XYZ periodicity")
938  END IF
939 
940  CALL get_qs_env(qs_env, pw_env=pw_env, para_env=para_env)
941 
942  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
943  poisson_env=poisson_env)
944  CALL auxbas_pw_pool%create_pw(rho_ga)
945  CALL auxbas_pw_pool%create_pw(va_gspace)
946  CALL auxbas_pw_pool%create_pw(va_rspace)
947 
948  !get fitting points and store them in resp_env%fitpoints
949  CALL get_fitting_points(qs_env, resp_env, rep_sys, particles=particles, &
950  cell=cell)
951  ALLOCATE (vpot(resp_env%npoints_proc, natom))
952  normalize_factor = sqrt((resp_env%eta/pi)**3)
953 
954  DO i = 1, natom
955  !collocate gaussian for each atom
956  CALL pw_zero(rho_ga)
957  CALL calculate_rho_resp_single(rho_ga, qs_env, resp_env%eta, i)
958  !calculate potential va and store the part needed for fitting in vpot
959  CALL pw_zero(va_gspace)
960  CALL pw_poisson_solve(poisson_env, rho_ga, vhartree=va_gspace)
961  CALL pw_zero(va_rspace)
962  CALL pw_transfer(va_gspace, va_rspace)
963  CALL pw_scale(va_rspace, normalize_factor)
964  DO ip = 1, resp_env%npoints_proc
965  jx = resp_env%fitpoints(1, ip)
966  jy = resp_env%fitpoints(2, ip)
967  jz = resp_env%fitpoints(3, ip)
968  vpot(ip, i) = va_rspace%array(jx, jy, jz)
969  END DO
970  END DO
971 
972  CALL va_gspace%release()
973  CALL va_rspace%release()
974  CALL rho_ga%release()
975 
976  DO i = 1, natom
977  DO j = 1, natom
978  ! calculate matrix
979  resp_env%matrix(i, j) = resp_env%matrix(i, j) + 2.0_dp*sum(vpot(:, i)*vpot(:, j))
980  END DO
981  ! calculate vector resp_env%rhs
982  CALL calculate_rhs(qs_env, resp_env, resp_env%rhs(i), vpot(:, i))
983  END DO
984 
985  CALL para_env%sum(resp_env%matrix)
986  CALL para_env%sum(resp_env%rhs)
987  !weighted sum
988  resp_env%matrix = resp_env%matrix/resp_env%npoints
989  resp_env%rhs = resp_env%rhs/resp_env%npoints
990 
991  ! REPEAT stuff
992  IF (resp_env%use_repeat_method) THEN
993  ! sum over selected points of single Gaussian potential vpot
994  DO i = 1, natom
995  resp_env%sum_vpot(i) = 2.0_dp*accurate_sum(vpot(:, i))/resp_env%npoints
996  END DO
997  CALL para_env%sum(resp_env%sum_vpot)
998  CALL para_env%sum(resp_env%sum_vhartree)
999  resp_env%sum_vhartree = 2.0_dp*resp_env%sum_vhartree/resp_env%npoints
1000  END IF
1001 
1002  DEALLOCATE (vpot)
1003  CALL timestop(handle)
1004 
1005  END SUBROUTINE calc_resp_matrix_periodic
1006 
1007 ! **************************************************************************************************
1008 !> \brief get RESP fitting points for the periodic fitting
1009 !> \param qs_env the qs environment
1010 !> \param resp_env the resp environment
1011 !> \param rep_sys structure for repeating input sections defining fit points
1012 !> \param particles ...
1013 !> \param cell parameters related to the simulation cell
1014 ! **************************************************************************************************
1015  SUBROUTINE get_fitting_points(qs_env, resp_env, rep_sys, particles, cell)
1016 
1017  TYPE(qs_environment_type), POINTER :: qs_env
1018  TYPE(resp_type), POINTER :: resp_env
1019  TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys
1020  TYPE(particle_list_type), POINTER :: particles
1021  TYPE(cell_type), POINTER :: cell
1022 
1023  CHARACTER(len=*), PARAMETER :: routinen = 'get_fitting_points'
1024 
1025  INTEGER :: bo(2, 3), gbo(2, 3), handle, i, iatom, &
1026  ikind, in_x, in_y, in_z, jx, jy, jz, &
1027  k, kind_number, l, m, natom, nkind, &
1028  now, p
1029  LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: not_in_range
1030  REAL(kind=dp) :: delta, dh(3, 3), r(3), rmax, rmin, &
1031  vec_pbc(3)
1032  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: dist
1033  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1034  TYPE(mp_para_env_type), POINTER :: para_env
1035  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1036  TYPE(pw_r3d_rs_type), POINTER :: v_hartree_pw
1037 
1038  CALL timeset(routinen, handle)
1039 
1040  NULLIFY (atomic_kind_set, v_hartree_pw, para_env, particle_set)
1041  delta = 1.0e-13_dp
1042 
1043  CALL get_qs_env(qs_env, &
1044  particle_set=particle_set, &
1045  atomic_kind_set=atomic_kind_set, &
1046  para_env=para_env, &
1047  v_hartree_rspace=v_hartree_pw)
1048 
1049  bo = v_hartree_pw%pw_grid%bounds_local
1050  gbo = v_hartree_pw%pw_grid%bounds
1051  dh = v_hartree_pw%pw_grid%dh
1052  natom = SIZE(particles%els)
1053  nkind = SIZE(atomic_kind_set)
1054 
1055  IF (.NOT. ASSOCIATED(resp_env%fitpoints)) THEN
1056  now = 1000
1057  ALLOCATE (resp_env%fitpoints(3, now))
1058  ELSE
1059  now = SIZE(resp_env%fitpoints, 2)
1060  END IF
1061 
1062  ALLOCATE (dist(natom))
1063  ALLOCATE (not_in_range(natom, 2))
1064 
1065  !every proc gets another bo, grid is distributed
1066  DO jz = bo(1, 3), bo(2, 3)
1067  IF (.NOT. (modulo(jz, resp_env%stride(3)) == 0)) cycle
1068  DO jy = bo(1, 2), bo(2, 2)
1069  IF (.NOT. (modulo(jy, resp_env%stride(2)) == 0)) cycle
1070  DO jx = bo(1, 1), bo(2, 1)
1071  IF (.NOT. (modulo(jx, resp_env%stride(1)) == 0)) cycle
1072  !bounds gbo reach from -np/2 to np/2. shift of np/2 so that r(1,1,1)=(0,0,0)
1073  l = jx - gbo(1, 1)
1074  k = jy - gbo(1, 2)
1075  p = jz - gbo(1, 3)
1076  r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
1077  r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
1078  r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
1079  IF (resp_env%molecular_sys) THEN
1080  not_in_range = .false.
1081  DO m = 1, natom
1082  vec_pbc = pbc(r, particles%els(m)%r, cell)
1083  dist(m) = sqrt(sum(vec_pbc**2))
1084  CALL get_atomic_kind(atomic_kind=particle_set(m)%atomic_kind, &
1085  kind_number=kind_number)
1086  DO ikind = 1, nkind
1087  IF (ikind == kind_number) THEN
1088  rmin = resp_env%rmin_kind(ikind)
1089  rmax = resp_env%rmax_kind(ikind)
1090  EXIT
1091  END IF
1092  END DO
1093  IF (dist(m) < rmin + delta) not_in_range(m, 1) = .true.
1094  IF (dist(m) > rmax - delta) not_in_range(m, 2) = .true.
1095  END DO
1096  IF (any(not_in_range(:, 1)) .OR. all(not_in_range(:, 2))) cycle
1097  ELSE
1098  DO i = 1, SIZE(rep_sys)
1099  DO m = 1, SIZE(rep_sys(i)%p_resp%atom_surf_list)
1100  in_z = 0
1101  in_y = 0
1102  in_x = 0
1103  iatom = rep_sys(i)%p_resp%atom_surf_list(m)
1104  SELECT CASE (rep_sys(i)%p_resp%my_fit)
1105  CASE (do_resp_x_dir, do_resp_y_dir, do_resp_z_dir)
1106  vec_pbc = pbc(particles%els(iatom)%r, r, cell)
1107  CASE (do_resp_minus_x_dir, do_resp_minus_y_dir, do_resp_minus_z_dir)
1108  vec_pbc = pbc(r, particles%els(iatom)%r, cell)
1109  END SELECT
1110  SELECT CASE (rep_sys(i)%p_resp%my_fit)
1111  !subtract delta=1.0E-13 to get rid of rounding errors when shifting atoms
1112  CASE (do_resp_x_dir, do_resp_minus_x_dir)
1113  IF (abs(vec_pbc(3)) < rep_sys(i)%p_resp%length - delta) in_z = 1
1114  IF (abs(vec_pbc(2)) < rep_sys(i)%p_resp%length - delta) in_y = 1
1115  IF (vec_pbc(1) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
1116  vec_pbc(1) < rep_sys(i)%p_resp%range_surf(2) - delta) in_x = 1
1117  CASE (do_resp_y_dir, do_resp_minus_y_dir)
1118  IF (abs(vec_pbc(3)) < rep_sys(i)%p_resp%length - delta) in_z = 1
1119  IF (vec_pbc(2) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
1120  vec_pbc(2) < rep_sys(i)%p_resp%range_surf(2) - delta) in_y = 1
1121  IF (abs(vec_pbc(1)) < rep_sys(i)%p_resp%length - delta) in_x = 1
1122  CASE (do_resp_z_dir, do_resp_minus_z_dir)
1123  IF (vec_pbc(3) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
1124  vec_pbc(3) < rep_sys(i)%p_resp%range_surf(2) - delta) in_z = 1
1125  IF (abs(vec_pbc(2)) < rep_sys(i)%p_resp%length - delta) in_y = 1
1126  IF (abs(vec_pbc(1)) < rep_sys(i)%p_resp%length - delta) in_x = 1
1127  END SELECT
1128  IF (in_z*in_y*in_x == 1) EXIT
1129  END DO
1130  IF (in_z*in_y*in_x == 1) EXIT
1131  END DO
1132  IF (in_z*in_y*in_x == 0) cycle
1133  END IF
1134  resp_env%npoints_proc = resp_env%npoints_proc + 1
1135  IF (resp_env%npoints_proc > now) THEN
1136  now = 2*now
1137  CALL reallocate(resp_env%fitpoints, 1, 3, 1, now)
1138  END IF
1139  resp_env%fitpoints(1, resp_env%npoints_proc) = jx
1140  resp_env%fitpoints(2, resp_env%npoints_proc) = jy
1141  resp_env%fitpoints(3, resp_env%npoints_proc) = jz
1142  END DO
1143  END DO
1144  END DO
1145 
1146  resp_env%npoints = resp_env%npoints_proc
1147  CALL para_env%sum(resp_env%npoints)
1148 
1149  DEALLOCATE (dist)
1150  DEALLOCATE (not_in_range)
1151 
1152  CALL timestop(handle)
1153 
1154  END SUBROUTINE get_fitting_points
1155 
1156 ! **************************************************************************************************
1157 !> \brief calculate vector rhs
1158 !> \param qs_env the qs environment
1159 !> \param resp_env the resp environment
1160 !> \param rhs vector
1161 !> \param vpot single gaussian potential
1162 ! **************************************************************************************************
1163  SUBROUTINE calculate_rhs(qs_env, resp_env, rhs, vpot)
1164 
1165  TYPE(qs_environment_type), POINTER :: qs_env
1166  TYPE(resp_type), POINTER :: resp_env
1167  REAL(kind=dp), INTENT(INOUT) :: rhs
1168  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: vpot
1169 
1170  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_rhs'
1171 
1172  INTEGER :: handle, ip, jx, jy, jz
1173  REAL(kind=dp) :: dvol
1174  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: vhartree
1175  TYPE(pw_r3d_rs_type), POINTER :: v_hartree_pw
1176 
1177  CALL timeset(routinen, handle)
1178 
1179  NULLIFY (v_hartree_pw)
1180  CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_pw)
1181  dvol = v_hartree_pw%pw_grid%dvol
1182  ALLOCATE (vhartree(resp_env%npoints_proc))
1183  vhartree = 0.0_dp
1184 
1185  !multiply v_hartree and va_rspace and calculate the vector rhs
1186  !taking into account that v_hartree has opposite site; remove v_qmmm
1187  DO ip = 1, resp_env%npoints_proc
1188  jx = resp_env%fitpoints(1, ip)
1189  jy = resp_env%fitpoints(2, ip)
1190  jz = resp_env%fitpoints(3, ip)
1191  vhartree(ip) = -v_hartree_pw%array(jx, jy, jz)/dvol
1192  IF (qs_env%qmmm) THEN
1193  !taking into account that v_qmmm has also opposite sign
1194  vhartree(ip) = vhartree(ip) + qs_env%ks_qmmm_env%v_qmmm_rspace%array(jx, jy, jz)
1195  END IF
1196  rhs = rhs + 2.0_dp*vhartree(ip)*vpot(ip)
1197  END DO
1198 
1199  IF (resp_env%use_repeat_method) THEN
1200  resp_env%sum_vhartree = accurate_sum(vhartree)
1201  END IF
1202 
1203  DEALLOCATE (vhartree)
1204 
1205  CALL timestop(handle)
1206 
1207  END SUBROUTINE calculate_rhs
1208 
1209 ! **************************************************************************************************
1210 !> \brief print the atom coordinates and the coordinates of the fitting points
1211 !> to an xyz file
1212 !> \param qs_env the qs environment
1213 !> \param resp_env the resp environment
1214 ! **************************************************************************************************
1215  SUBROUTINE print_fitting_points(qs_env, resp_env)
1216 
1217  TYPE(qs_environment_type), POINTER :: qs_env
1218  TYPE(resp_type), POINTER :: resp_env
1219 
1220  CHARACTER(len=*), PARAMETER :: routinen = 'print_fitting_points'
1221 
1222  CHARACTER(LEN=2) :: element_symbol
1223  CHARACTER(LEN=default_path_length) :: filename
1224  INTEGER :: gbo(2, 3), handle, i, iatom, ip, jx, jy, &
1225  jz, k, l, my_pos, nobjects, &
1226  output_unit, p
1227  INTEGER, DIMENSION(:), POINTER :: tmp_npoints, tmp_size
1228  INTEGER, DIMENSION(:, :), POINTER :: tmp_points
1229  REAL(kind=dp) :: conv, dh(3, 3), r(3)
1230  TYPE(cp_logger_type), POINTER :: logger
1231  TYPE(mp_para_env_type), POINTER :: para_env
1232  TYPE(mp_request_type), DIMENSION(6) :: req
1233  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1234  TYPE(pw_r3d_rs_type), POINTER :: v_hartree_pw
1235  TYPE(section_vals_type), POINTER :: input, print_key, resp_section
1236 
1237  CALL timeset(routinen, handle)
1238 
1239  NULLIFY (para_env, input, logger, resp_section, print_key, particle_set, tmp_size, &
1240  tmp_points, tmp_npoints, v_hartree_pw)
1241 
1242  CALL get_qs_env(qs_env, input=input, para_env=para_env, &
1243  particle_set=particle_set, v_hartree_rspace=v_hartree_pw)
1244  conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
1245  gbo = v_hartree_pw%pw_grid%bounds
1246  dh = v_hartree_pw%pw_grid%dh
1247  nobjects = SIZE(particle_set) + resp_env%npoints
1248 
1249  resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1250  print_key => section_vals_get_subs_vals(resp_section, "PRINT%COORD_FIT_POINTS")
1251  logger => cp_get_default_logger()
1252  output_unit = cp_print_key_unit_nr(logger, resp_section, &
1253  "PRINT%COORD_FIT_POINTS", &
1254  extension=".xyz", &
1255  file_status="REPLACE", &
1256  file_action="WRITE", &
1257  file_form="FORMATTED")
1258 
1259  IF (btest(cp_print_key_should_output(logger%iter_info, &
1260  resp_section, "PRINT%COORD_FIT_POINTS"), &
1261  cp_p_file)) THEN
1262  IF (output_unit > 0) THEN
1263  filename = cp_print_key_generate_filename(logger, &
1264  print_key, extension=".xyz", &
1265  my_local=.false.)
1266  WRITE (unit=output_unit, fmt="(I12,/)") nobjects
1267  DO iatom = 1, SIZE(particle_set)
1268  CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
1269  element_symbol=element_symbol)
1270  WRITE (unit=output_unit, fmt="(A,1X,3F10.5)") element_symbol, &
1271  particle_set(iatom)%r(1:3)*conv
1272  END DO
1273  !printing points of proc which is doing the output (should be proc 0)
1274  DO ip = 1, resp_env%npoints_proc
1275  jx = resp_env%fitpoints(1, ip)
1276  jy = resp_env%fitpoints(2, ip)
1277  jz = resp_env%fitpoints(3, ip)
1278  l = jx - gbo(1, 1)
1279  k = jy - gbo(1, 2)
1280  p = jz - gbo(1, 3)
1281  r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
1282  r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
1283  r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
1284  r(:) = r(:)*conv
1285  WRITE (unit=output_unit, fmt="(A,2X,3F10.5)") "X", r(1), r(2), r(3)
1286  END DO
1287  END IF
1288 
1289  ALLOCATE (tmp_size(1))
1290  ALLOCATE (tmp_npoints(1))
1291 
1292  !sending data of all other procs to proc which makes the output (proc 0)
1293  IF (output_unit > 0) THEN
1294  my_pos = para_env%mepos
1295  DO i = 1, para_env%num_pe
1296  IF (my_pos == i - 1) cycle
1297  CALL para_env%irecv(msgout=tmp_size, source=i - 1, &
1298  request=req(1))
1299  CALL req(1)%wait()
1300  ALLOCATE (tmp_points(3, tmp_size(1)))
1301  CALL para_env%irecv(msgout=tmp_points, source=i - 1, &
1302  request=req(3))
1303  CALL req(3)%wait()
1304  CALL para_env%irecv(msgout=tmp_npoints, source=i - 1, &
1305  request=req(5))
1306  CALL req(5)%wait()
1307  DO ip = 1, tmp_npoints(1)
1308  jx = tmp_points(1, ip)
1309  jy = tmp_points(2, ip)
1310  jz = tmp_points(3, ip)
1311  l = jx - gbo(1, 1)
1312  k = jy - gbo(1, 2)
1313  p = jz - gbo(1, 3)
1314  r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
1315  r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
1316  r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
1317  r(:) = r(:)*conv
1318  WRITE (unit=output_unit, fmt="(A,2X,3F10.5)") "X", r(1), r(2), r(3)
1319  END DO
1320  DEALLOCATE (tmp_points)
1321  END DO
1322  ELSE
1323  tmp_size(1) = SIZE(resp_env%fitpoints, 2)
1324  !para_env%source should be 0
1325  CALL para_env%isend(msgin=tmp_size, dest=para_env%source, &
1326  request=req(2))
1327  CALL req(2)%wait()
1328  CALL para_env%isend(msgin=resp_env%fitpoints, dest=para_env%source, &
1329  request=req(4))
1330  CALL req(4)%wait()
1331  tmp_npoints(1) = resp_env%npoints_proc
1332  CALL para_env%isend(msgin=tmp_npoints, dest=para_env%source, &
1333  request=req(6))
1334  CALL req(6)%wait()
1335  END IF
1336 
1337  DEALLOCATE (tmp_size)
1338  DEALLOCATE (tmp_npoints)
1339  END IF
1340 
1341  CALL cp_print_key_finished_output(output_unit, logger, resp_section, &
1342  "PRINT%COORD_FIT_POINTS")
1343 
1344  CALL timestop(handle)
1345 
1346  END SUBROUTINE print_fitting_points
1347 
1348 ! **************************************************************************************************
1349 !> \brief add restraints and constraints
1350 !> \param qs_env the qs environment
1351 !> \param resp_env the resp environment
1352 !> \param rest_section input section for restraints
1353 !> \param subsys ...
1354 !> \param natom number of atoms
1355 !> \param cons_section input section for constraints
1356 !> \param particle_set ...
1357 ! **************************************************************************************************
1358  SUBROUTINE add_restraints_and_constraints(qs_env, resp_env, rest_section, &
1359  subsys, natom, cons_section, particle_set)
1360 
1361  TYPE(qs_environment_type), POINTER :: qs_env
1362  TYPE(resp_type), POINTER :: resp_env
1363  TYPE(section_vals_type), POINTER :: rest_section
1364  TYPE(qs_subsys_type), POINTER :: subsys
1365  INTEGER, INTENT(IN) :: natom
1366  TYPE(section_vals_type), POINTER :: cons_section
1367  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1368 
1369  CHARACTER(len=*), PARAMETER :: routinen = 'add_restraints_and_constraints'
1370 
1371  INTEGER :: handle, i, k, m, ncons_v, z
1372  INTEGER, DIMENSION(:), POINTER :: atom_list_cons, atom_list_res
1373  LOGICAL :: explicit_coeff
1374  REAL(kind=dp) :: my_atom_coef(2), strength, TARGET
1375  REAL(kind=dp), DIMENSION(:), POINTER :: atom_coef
1376  TYPE(dft_control_type), POINTER :: dft_control
1377 
1378  CALL timeset(routinen, handle)
1379 
1380  NULLIFY (atom_coef, atom_list_res, atom_list_cons, dft_control)
1381 
1382  CALL get_qs_env(qs_env, dft_control=dft_control)
1383 
1384  !*** add the restraints
1385  DO i = 1, resp_env%nrest_sec
1386  CALL section_vals_val_get(rest_section, "TARGET", i_rep_section=i, r_val=TARGET)
1387  CALL section_vals_val_get(rest_section, "STRENGTH", i_rep_section=i, r_val=strength)
1388  CALL build_atom_list(rest_section, subsys, atom_list_res, i)
1389  CALL section_vals_val_get(rest_section, "ATOM_COEF", i_rep_section=i, explicit=explicit_coeff)
1390  IF (explicit_coeff) THEN
1391  CALL section_vals_val_get(rest_section, "ATOM_COEF", i_rep_section=i, r_vals=atom_coef)
1392  cpassert(SIZE(atom_list_res) == SIZE(atom_coef))
1393  END IF
1394  DO m = 1, SIZE(atom_list_res)
1395  IF (explicit_coeff) THEN
1396  DO k = 1, SIZE(atom_list_res)
1397  resp_env%matrix(atom_list_res(m), atom_list_res(k)) = &
1398  resp_env%matrix(atom_list_res(m), atom_list_res(k)) + &
1399  atom_coef(m)*atom_coef(k)*2.0_dp*strength
1400  END DO
1401  resp_env%rhs(atom_list_res(m)) = resp_env%rhs(atom_list_res(m)) + &
1402  2.0_dp*TARGET*strength*atom_coef(m)
1403  ELSE
1404  resp_env%matrix(atom_list_res(m), atom_list_res(m)) = &
1405  resp_env%matrix(atom_list_res(m), atom_list_res(m)) + &
1406  2.0_dp*strength
1407  resp_env%rhs(atom_list_res(m)) = resp_env%rhs(atom_list_res(m)) + &
1408  2.0_dp*TARGET*strength
1409  END IF
1410  END DO
1411  DEALLOCATE (atom_list_res)
1412  END DO
1413 
1414  ! if heavies are restrained to zero, add these as well
1415  IF (resp_env%rheavies) THEN
1416  DO i = 1, natom
1417  CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, z=z)
1418  IF (z .NE. 1) THEN
1419  resp_env%matrix(i, i) = resp_env%matrix(i, i) + 2.0_dp*resp_env%rheavies_strength
1420  END IF
1421  END DO
1422  END IF
1423 
1424  !*** add the constraints
1425  ncons_v = 0
1426  ncons_v = ncons_v + natom
1427 
1428  ! REPEAT charges: treat the offset like a constraint
1429  IF (resp_env%use_repeat_method) THEN
1430  ncons_v = ncons_v + 1
1431  resp_env%matrix(1:natom, ncons_v) = resp_env%sum_vpot(1:natom)
1432  resp_env%matrix(ncons_v, 1:natom) = resp_env%sum_vpot(1:natom)
1433  resp_env%matrix(ncons_v, ncons_v) = 2.0_dp
1434  resp_env%rhs(ncons_v) = resp_env%sum_vhartree
1435  END IF
1436 
1437  ! total charge constraint
1438  IF (resp_env%itc) THEN
1439  ncons_v = ncons_v + 1
1440  resp_env%matrix(1:natom, ncons_v) = 1.0_dp
1441  resp_env%matrix(ncons_v, 1:natom) = 1.0_dp
1442  resp_env%rhs(ncons_v) = dft_control%charge
1443  END IF
1444 
1445  ! explicit constraints
1446  DO i = 1, resp_env%ncons_sec
1447  CALL build_atom_list(cons_section, subsys, atom_list_cons, i)
1448  IF (.NOT. resp_env%equal_charges) THEN
1449  ncons_v = ncons_v + 1
1450  CALL section_vals_val_get(cons_section, "ATOM_COEF", i_rep_section=i, r_vals=atom_coef)
1451  CALL section_vals_val_get(cons_section, "TARGET", i_rep_section=i, r_val=TARGET)
1452  cpassert(SIZE(atom_list_cons) == SIZE(atom_coef))
1453  DO m = 1, SIZE(atom_list_cons)
1454  resp_env%matrix(atom_list_cons(m), ncons_v) = atom_coef(m)
1455  resp_env%matrix(ncons_v, atom_list_cons(m)) = atom_coef(m)
1456  END DO
1457  resp_env%rhs(ncons_v) = TARGET
1458  ELSE
1459  my_atom_coef(1) = 1.0_dp
1460  my_atom_coef(2) = -1.0_dp
1461  DO k = 2, SIZE(atom_list_cons)
1462  ncons_v = ncons_v + 1
1463  resp_env%matrix(atom_list_cons(1), ncons_v) = my_atom_coef(1)
1464  resp_env%matrix(ncons_v, atom_list_cons(1)) = my_atom_coef(1)
1465  resp_env%matrix(atom_list_cons(k), ncons_v) = my_atom_coef(2)
1466  resp_env%matrix(ncons_v, atom_list_cons(k)) = my_atom_coef(2)
1467  resp_env%rhs(ncons_v) = 0.0_dp
1468  END DO
1469  END IF
1470  DEALLOCATE (atom_list_cons)
1471  END DO
1472  CALL timestop(handle)
1473 
1474  END SUBROUTINE add_restraints_and_constraints
1475 
1476 ! **************************************************************************************************
1477 !> \brief print input information
1478 !> \param qs_env the qs environment
1479 !> \param resp_env the resp environment
1480 !> \param rep_sys structure for repeating input sections defining fit points
1481 !> \param my_per ...
1482 ! **************************************************************************************************
1483  SUBROUTINE print_resp_parameter_info(qs_env, resp_env, rep_sys, my_per)
1484 
1485  TYPE(qs_environment_type), POINTER :: qs_env
1486  TYPE(resp_type), POINTER :: resp_env
1487  TYPE(resp_p_type), DIMENSION(:), POINTER :: rep_sys
1488  INTEGER, INTENT(IN) :: my_per
1489 
1490  CHARACTER(len=*), PARAMETER :: routinen = 'print_resp_parameter_info'
1491 
1492  CHARACTER(len=2) :: symbol
1493  INTEGER :: handle, i, ikind, kind_number, nkinds, &
1494  output_unit
1495  REAL(kind=dp) :: conv, eta_conv
1496  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1497  TYPE(cp_logger_type), POINTER :: logger
1498  TYPE(section_vals_type), POINTER :: input, resp_section
1499 
1500  CALL timeset(routinen, handle)
1501  NULLIFY (logger, input, resp_section)
1502 
1503  CALL get_qs_env(qs_env, &
1504  input=input, &
1505  atomic_kind_set=atomic_kind_set)
1506  resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1507  logger => cp_get_default_logger()
1508  output_unit = cp_print_key_unit_nr(logger, resp_section, "PRINT%PROGRAM_RUN_INFO", &
1509  extension=".resp")
1510  nkinds = SIZE(atomic_kind_set)
1511 
1512  conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
1513  IF (.NOT. my_per == use_perd_none) THEN
1514  eta_conv = cp_unit_from_cp2k(resp_env%eta, "angstrom", power=-2)
1515  END IF
1516 
1517  IF (output_unit > 0) THEN
1518  WRITE (output_unit, '(/,1X,A,/)') "STARTING RESP FIT"
1519  IF (resp_env%use_repeat_method) THEN
1520  WRITE (output_unit, '(T3,A)') &
1521  "Fit the variance of the potential (REPEAT method)."
1522  END IF
1523  IF (.NOT. resp_env%equal_charges) THEN
1524  WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons_sec
1525  ELSE
1526  IF (resp_env%itc) THEN
1527  WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons - 1
1528  ELSE
1529  WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons
1530  END IF
1531  END IF
1532  WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit restraints: ", resp_env%nrest_sec
1533  WRITE (output_unit, '(T3,A,T80,A)') "Constrain total charge ", merge("T", "F", resp_env%itc)
1534  WRITE (output_unit, '(T3,A,T80,A)') "Restrain heavy atoms ", merge("T", "F", resp_env%rheavies)
1535  IF (resp_env%rheavies) THEN
1536  WRITE (output_unit, '(T3,A,T71,F10.6)') "Heavy atom restraint strength: ", &
1537  resp_env%rheavies_strength
1538  END IF
1539  WRITE (output_unit, '(T3,A,T66,3I5)') "Stride: ", resp_env%stride
1540  IF (resp_env%molecular_sys) THEN
1541  WRITE (output_unit, '(T3,A)') &
1542  "------------------------------------------------------------------------------"
1543  WRITE (output_unit, '(T3,A)') "Using sphere sampling"
1544  WRITE (output_unit, '(T3,A,T46,A,T66,A)') &
1545  "Element", "RMIN [angstrom]", "RMAX [angstrom]"
1546  DO ikind = 1, nkinds
1547  CALL get_atomic_kind(atomic_kind=atomic_kind_set(ikind), &
1548  kind_number=kind_number, &
1549  element_symbol=symbol)
1550  WRITE (output_unit, '(T3,A,T51,F10.5,T71,F10.5)') &
1551  symbol, &
1552  resp_env%rmin_kind(kind_number)*conv, &
1553  resp_env%rmax_kind(kind_number)*conv
1554  END DO
1555  IF (my_per == use_perd_none) THEN
1556  WRITE (output_unit, '(T3,A,T51,3F10.5)') "Box min [angstrom]: ", resp_env%box_low(1:3)*conv
1557  WRITE (output_unit, '(T3,A,T51,3F10.5)') "Box max [angstrom]: ", resp_env%box_hi(1:3)*conv
1558  END IF
1559  WRITE (output_unit, '(T3,A)') &
1560  "------------------------------------------------------------------------------"
1561  ELSE
1562  WRITE (output_unit, '(T3,A)') &
1563  "------------------------------------------------------------------------------"
1564  WRITE (output_unit, '(T3,A)') "Using slab sampling"
1565  WRITE (output_unit, '(2X,A,F10.5)') "Index of atoms defining the surface: "
1566  DO i = 1, SIZE(rep_sys)
1567  IF (i > 1 .AND. all(rep_sys(i)%p_resp%atom_surf_list == rep_sys(1)%p_resp%atom_surf_list)) EXIT
1568  WRITE (output_unit, '(7X,10I6)') rep_sys(i)%p_resp%atom_surf_list
1569  END DO
1570  DO i = 1, SIZE(rep_sys)
1571  IF (i > 1 .AND. all(rep_sys(i)%p_resp%range_surf == rep_sys(1)%p_resp%range_surf)) EXIT
1572  WRITE (output_unit, '(T3,A,T61,2F10.5)') &
1573  "Range for sampling above the surface [angstrom]:", &
1574  rep_sys(i)%p_resp%range_surf(1:2)*conv
1575  END DO
1576  DO i = 1, SIZE(rep_sys)
1577  IF (i > 1 .AND. rep_sys(i)%p_resp%length == rep_sys(1)%p_resp%length) EXIT
1578  WRITE (output_unit, '(T3,A,T71,F10.5)') "Length of sampling box above each"// &
1579  " surface atom [angstrom]: ", rep_sys(i)%p_resp%length*conv
1580  END DO
1581  WRITE (output_unit, '(T3,A)') &
1582  "------------------------------------------------------------------------------"
1583  END IF
1584  IF (.NOT. my_per == use_perd_none) THEN
1585  WRITE (output_unit, '(T3,A,T71,F10.5)') "Width of Gaussian charge"// &
1586  " distribution [angstrom^-2]: ", eta_conv
1587  END IF
1588  CALL m_flush(output_unit)
1589  END IF
1590  CALL cp_print_key_finished_output(output_unit, logger, resp_section, &
1591  "PRINT%PROGRAM_RUN_INFO")
1592 
1593  CALL timestop(handle)
1594 
1595  END SUBROUTINE print_resp_parameter_info
1596 
1597 ! **************************************************************************************************
1598 !> \brief print RESP charges to an extra file or to the normal output file
1599 !> \param qs_env the qs environment
1600 !> \param resp_env the resp environment
1601 !> \param output_runinfo ...
1602 !> \param natom number of atoms
1603 ! **************************************************************************************************
1604  SUBROUTINE print_resp_charges(qs_env, resp_env, output_runinfo, natom)
1605 
1606  TYPE(qs_environment_type), POINTER :: qs_env
1607  TYPE(resp_type), POINTER :: resp_env
1608  INTEGER, INTENT(IN) :: output_runinfo, natom
1609 
1610  CHARACTER(len=*), PARAMETER :: routinen = 'print_resp_charges'
1611 
1612  CHARACTER(LEN=default_path_length) :: filename
1613  INTEGER :: handle, output_file
1614  TYPE(cp_logger_type), POINTER :: logger
1615  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1616  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1617  TYPE(section_vals_type), POINTER :: input, print_key, resp_section
1618 
1619  CALL timeset(routinen, handle)
1620 
1621  NULLIFY (particle_set, qs_kind_set, input, logger, resp_section, print_key)
1622 
1623  CALL get_qs_env(qs_env, input=input, particle_set=particle_set, &
1624  qs_kind_set=qs_kind_set)
1625 
1626  resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1627  print_key => section_vals_get_subs_vals(resp_section, &
1628  "PRINT%RESP_CHARGES_TO_FILE")
1629  logger => cp_get_default_logger()
1630 
1631  IF (btest(cp_print_key_should_output(logger%iter_info, &
1632  resp_section, "PRINT%RESP_CHARGES_TO_FILE"), &
1633  cp_p_file)) THEN
1634  output_file = cp_print_key_unit_nr(logger, resp_section, &
1635  "PRINT%RESP_CHARGES_TO_FILE", &
1636  extension=".resp", &
1637  file_status="REPLACE", &
1638  file_action="WRITE", &
1639  file_form="FORMATTED")
1640  IF (output_file > 0) THEN
1641  filename = cp_print_key_generate_filename(logger, &
1642  print_key, extension=".resp", &
1643  my_local=.false.)
1644  CALL print_atomic_charges(particle_set, qs_kind_set, output_file, title="RESP charges:", &
1645  atomic_charges=resp_env%rhs(1:natom))
1646  IF (output_runinfo > 0) WRITE (output_runinfo, '(2X,A,/)') "PRINTED RESP CHARGES TO FILE"
1647  END IF
1648 
1649  CALL cp_print_key_finished_output(output_file, logger, resp_section, &
1650  "PRINT%RESP_CHARGES_TO_FILE")
1651  ELSE
1652  CALL print_atomic_charges(particle_set, qs_kind_set, output_runinfo, title="RESP charges:", &
1653  atomic_charges=resp_env%rhs(1:natom))
1654  END IF
1655 
1656  CALL timestop(handle)
1657 
1658  END SUBROUTINE print_resp_charges
1659 
1660 ! **************************************************************************************************
1661 !> \brief print potential generated by RESP charges to file
1662 !> \param qs_env the qs environment
1663 !> \param resp_env the resp environment
1664 !> \param particles ...
1665 !> \param natom number of atoms
1666 !> \param output_runinfo ...
1667 ! **************************************************************************************************
1668  SUBROUTINE print_pot_from_resp_charges(qs_env, resp_env, particles, natom, output_runinfo)
1669 
1670  TYPE(qs_environment_type), POINTER :: qs_env
1671  TYPE(resp_type), POINTER :: resp_env
1672  TYPE(particle_list_type), POINTER :: particles
1673  INTEGER, INTENT(IN) :: natom, output_runinfo
1674 
1675  CHARACTER(len=*), PARAMETER :: routinen = 'print_pot_from_resp_charges'
1676 
1677  CHARACTER(LEN=default_path_length) :: my_pos_cube
1678  INTEGER :: handle, ip, jx, jy, jz, unit_nr
1679  LOGICAL :: append_cube, mpi_io
1680  REAL(kind=dp) :: dvol, normalize_factor, rms, rrms, &
1681  sum_diff, sum_hartree, udvol
1682  TYPE(cp_logger_type), POINTER :: logger
1683  TYPE(mp_para_env_type), POINTER :: para_env
1684  TYPE(pw_c1d_gs_type) :: rho_resp, v_resp_gspace
1685  TYPE(pw_env_type), POINTER :: pw_env
1686  TYPE(pw_poisson_type), POINTER :: poisson_env
1687  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1688  TYPE(pw_r3d_rs_type) :: aux_r, v_resp_rspace
1689  TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
1690  TYPE(section_vals_type), POINTER :: input, print_key, resp_section
1691 
1692  CALL timeset(routinen, handle)
1693 
1694  NULLIFY (auxbas_pw_pool, logger, pw_env, poisson_env, input, print_key, &
1695  para_env, resp_section, v_hartree_rspace)
1696  CALL get_qs_env(qs_env, &
1697  input=input, &
1698  para_env=para_env, &
1699  pw_env=pw_env, &
1700  v_hartree_rspace=v_hartree_rspace)
1701  normalize_factor = sqrt((resp_env%eta/pi)**3)
1702  resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1703  print_key => section_vals_get_subs_vals(resp_section, &
1704  "PRINT%V_RESP_CUBE")
1705  logger => cp_get_default_logger()
1706 
1707  !*** calculate potential generated from RESP charges
1708  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1709  poisson_env=poisson_env)
1710 
1711  CALL auxbas_pw_pool%create_pw(rho_resp)
1712  CALL auxbas_pw_pool%create_pw(v_resp_gspace)
1713  CALL auxbas_pw_pool%create_pw(v_resp_rspace)
1714 
1715  CALL pw_zero(rho_resp)
1716  CALL calculate_rho_resp_all(rho_resp, resp_env%rhs, natom, &
1717  resp_env%eta, qs_env)
1718  CALL pw_zero(v_resp_gspace)
1719  CALL pw_poisson_solve(poisson_env, rho_resp, &
1720  vhartree=v_resp_gspace)
1721  CALL pw_zero(v_resp_rspace)
1722  CALL pw_transfer(v_resp_gspace, v_resp_rspace)
1723  dvol = v_resp_rspace%pw_grid%dvol
1724  CALL pw_scale(v_resp_rspace, dvol)
1725  CALL pw_scale(v_resp_rspace, -normalize_factor)
1726  ! REPEAT: correct for offset, take into account that potentials have reverse sign
1727  ! and are scaled by dvol
1728  IF (resp_env%use_repeat_method) THEN
1729  v_resp_rspace%array(:, :, :) = v_resp_rspace%array(:, :, :) - resp_env%offset*dvol
1730  END IF
1731  CALL v_resp_gspace%release()
1732  CALL rho_resp%release()
1733 
1734  !***now print the v_resp_rspace%pw to a cube file if requested
1735  IF (btest(cp_print_key_should_output(logger%iter_info, resp_section, &
1736  "PRINT%V_RESP_CUBE"), cp_p_file)) THEN
1737  CALL auxbas_pw_pool%create_pw(aux_r)
1738  append_cube = section_get_lval(resp_section, "PRINT%V_RESP_CUBE%APPEND")
1739  my_pos_cube = "REWIND"
1740  IF (append_cube) THEN
1741  my_pos_cube = "APPEND"
1742  END IF
1743  mpi_io = .true.
1744  unit_nr = cp_print_key_unit_nr(logger, resp_section, &
1745  "PRINT%V_RESP_CUBE", &
1746  extension=".cube", &
1747  file_position=my_pos_cube, &
1748  mpi_io=mpi_io)
1749  udvol = 1.0_dp/dvol
1750  CALL pw_copy(v_resp_rspace, aux_r)
1751  CALL pw_scale(aux_r, udvol)
1752  CALL cp_pw_to_cube(aux_r, unit_nr, "RESP POTENTIAL", particles=particles, &
1753  stride=section_get_ivals(resp_section, &
1754  "PRINT%V_RESP_CUBE%STRIDE"), &
1755  mpi_io=mpi_io)
1756  CALL cp_print_key_finished_output(unit_nr, logger, resp_section, &
1757  "PRINT%V_RESP_CUBE", mpi_io=mpi_io)
1758  CALL auxbas_pw_pool%give_back_pw(aux_r)
1759  END IF
1760 
1761  !*** RMS and RRMS
1762  sum_diff = 0.0_dp
1763  sum_hartree = 0.0_dp
1764  rms = 0.0_dp
1765  rrms = 0.0_dp
1766  DO ip = 1, resp_env%npoints_proc
1767  jx = resp_env%fitpoints(1, ip)
1768  jy = resp_env%fitpoints(2, ip)
1769  jz = resp_env%fitpoints(3, ip)
1770  sum_diff = sum_diff + (v_hartree_rspace%array(jx, jy, jz) - &
1771  v_resp_rspace%array(jx, jy, jz))**2
1772  sum_hartree = sum_hartree + v_hartree_rspace%array(jx, jy, jz)**2
1773  END DO
1774  CALL para_env%sum(sum_diff)
1775  CALL para_env%sum(sum_hartree)
1776  rms = sqrt(sum_diff/resp_env%npoints)
1777  rrms = sqrt(sum_diff/sum_hartree)
1778  IF (output_runinfo > 0) THEN
1779  WRITE (output_runinfo, '(2X,A,T69,ES12.5)') "Root-mean-square (RMS) "// &
1780  "error of RESP fit:", rms
1781  WRITE (output_runinfo, '(2X,A,T69,ES12.5,/)') "Relative root-mean-square "// &
1782  "(RRMS) error of RESP fit:", rrms
1783  END IF
1784 
1785  CALL v_resp_rspace%release()
1786 
1787  CALL timestop(handle)
1788 
1789  END SUBROUTINE print_pot_from_resp_charges
1790 
1791 END MODULE qs_resp
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
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:117
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:28
bibliography::golze2015
integer, save, public golze2015
Definition: bibliography.F:43
bibliography::rappe1992
integer, save, public rappe1992
Definition: bibliography.F:43
bibliography::campana2009
integer, save, public campana2009
Definition: bibliography.F:43
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:195
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:803
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:704
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:1013
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, stride, zero_tails, silent, mpi_io)
...
Definition: cp_realspace_grid_cube.F:45
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:1179
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:1150
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_resp_minus_x_dir
integer, parameter, public do_resp_minus_x_dir
Definition: input_constants.F:1086
input_constants::do_resp_minus_y_dir
integer, parameter, public do_resp_minus_y_dir
Definition: input_constants.F:1086
input_constants::do_resp_z_dir
integer, parameter, public do_resp_z_dir
Definition: input_constants.F:1086
input_constants::do_resp_minus_z_dir
integer, parameter, public do_resp_minus_z_dir
Definition: input_constants.F:1086
input_constants::use_uff_vdw_radii
integer, parameter, public use_uff_vdw_radii
Definition: input_constants.F:1086
input_constants::do_resp_y_dir
integer, parameter, public do_resp_y_dir
Definition: input_constants.F:1086
input_constants::do_resp_x_dir
integer, parameter, public do_resp_x_dir
Definition: input_constants.F:1086
input_constants::use_cambridge_vdw_radii
integer, parameter, public use_cambridge_vdw_radii
Definition: input_constants.F:1086
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:989
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:724
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:697
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:1040
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition: input_section_types.F:1005
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:106
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:1398
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, 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, 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, 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, rhs)
Set the QUICKSTEP environment.
Definition: qs_environment_types.F:1035
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:134
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