(git:374b731)
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-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! **************************************************************************************************
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, 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
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) &
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
1791END 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:117
cell_types::pbc
Definition cell_types.F:92
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:145
pw_methods::pw_scale
Definition pw_methods.F:93
pw_methods::pw_transfer
Definition pw_methods.F:303
pw_methods::pw_zero
Definition pw_methods.F:82
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
qs_collocate_density::calculate_rho_resp_all
Definition qs_collocate_density.F:142
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::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_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_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
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:55
cp_control_types::dft_control_type
Definition cp_control_types.F:559
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:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
message_passing::mp_request_type
Definition message_passing.F:623
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:122
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:215
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56