Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 = .FALSE., itc = .FALSE., &
100 : molecular_sys = .FALSE., rheavies = .FALSE., &
101 : use_repeat_method = .FALSE.
102 : INTEGER :: nres = -1, ncons = -1, &
103 : nrest_sec = -1, ncons_sec = -1, &
104 : npoints = -1, stride(3) = -1, my_fit = -1, &
105 : npoints_proc = -1, &
106 : auto_vdw_radii_table = -1
107 : INTEGER, DIMENSION(:), POINTER :: atom_surf_list => NULL()
108 : INTEGER, DIMENSION(:, :), POINTER :: fitpoints => NULL()
109 : REAL(KIND=dp) :: rheavies_strength = -1.0_dp, &
110 : length = -1.0_dp, eta = -1.0_dp, &
111 : sum_vhartree = -1.0_dp, offset = -1.0_dp
112 : REAL(KIND=dp), DIMENSION(3) :: box_hi = -1.0_dp, box_low = -1.0_dp
113 : REAL(KIND=dp), DIMENSION(:), POINTER :: rmin_kind => NULL(), &
114 : rmax_kind => NULL()
115 : REAL(KIND=dp), DIMENSION(:), POINTER :: range_surf => NULL()
116 : REAL(KIND=dp), DIMENSION(:), POINTER :: rhs => NULL()
117 : REAL(KIND=dp), DIMENSION(:), POINTER :: sum_vpot => NULL()
118 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix => NULL()
119 : END TYPE resp_type
120 :
121 : TYPE resp_p_type
122 : TYPE(resp_type), POINTER :: p_resp => NULL()
123 : END TYPE resp_p_type
124 :
125 : CONTAINS
126 :
127 : ! **************************************************************************************************
128 : !> \brief performs resp fit and generates RESP charges
129 : !> \param qs_env the qs environment
130 : ! **************************************************************************************************
131 9901 : 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 9901 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
139 : LOGICAL :: has_resp
140 9901 : REAL(KIND=dp), DIMENSION(:), POINTER :: rhs_to_save
141 9901 : 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 9901 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
147 : TYPE(qs_subsys_type), POINTER :: subsys
148 9901 : 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 9901 : CALL timeset(routineN, handle)
154 :
155 9901 : NULLIFY (logger, atomic_kind_set, cell, subsys, particles, particle_set, input, &
156 9901 : resp_section, cons_section, rest_section, poisson_section, resp_env, rep_sys)
157 :
158 9901 : CPASSERT(ASSOCIATED(qs_env))
159 :
160 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, input=input, &
161 9901 : subsys=subsys, particle_set=particle_set, cell=cell)
162 9901 : resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
163 9901 : CALL section_vals_get(resp_section, explicit=has_resp)
164 :
165 9901 : IF (has_resp) THEN
166 14 : logger => cp_get_default_logger()
167 14 : poisson_section => section_vals_get_subs_vals(input, "DFT%POISSON")
168 14 : CALL section_vals_val_get(poisson_section, "PERIODIC", i_val=my_per)
169 14 : 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 14 : cell, resp_section, cons_section, rest_section)
173 :
174 : !print info
175 14 : CALL print_resp_parameter_info(qs_env, resp_env, rep_sys, my_per)
176 :
177 14 : CALL qs_subsys_get(subsys, particles=particles)
178 14 : natom = particles%n_els
179 14 : nvar = natom + resp_env%ncons
180 :
181 14 : CALL resp_allocate(resp_env, natom, nvar)
182 42 : ALLOCATE (ipiv(nvar))
183 138 : ipiv = 0
184 :
185 : ! calculate the matrix and the vector rhs
186 4 : 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 4 : resp_env%matrix, resp_env%rhs, natom)
190 : CASE (use_perd_xyz)
191 10 : CALL cite_reference(Golze2015)
192 10 : IF (resp_env%use_repeat_method) CALL cite_reference(Campana2009)
193 10 : 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 14 : " or XYZ periodicity!")
198 : END SELECT
199 :
200 : output_unit = cp_print_key_unit_nr(logger, resp_section, "PRINT%PROGRAM_RUN_INFO", &
201 14 : extension=".resp")
202 14 : IF (output_unit > 0) THEN
203 : WRITE (output_unit, '(T3,A,T69,I12)') "Number of fitting points "// &
204 7 : "found: ", resp_env%npoints
205 7 : 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 14 : subsys, natom, cons_section, particle_set)
211 :
212 : !solve system for the values of the charges and the lagrangian multipliers
213 14 : CALL DGETRF(nvar, nvar, resp_env%matrix, nvar, ipiv, info)
214 14 : CPASSERT(info == 0)
215 :
216 14 : CALL DGETRS('N', nvar, 1, resp_env%matrix, nvar, ipiv, resp_env%rhs, nvar, info)
217 14 : CPASSERT(info == 0)
218 :
219 14 : IF (resp_env%use_repeat_method) resp_env%offset = resp_env%rhs(natom + 1)
220 14 : CALL print_resp_charges(qs_env, resp_env, output_unit, natom)
221 14 : CALL print_fitting_points(qs_env, resp_env)
222 14 : 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 14 : NULLIFY (dft_control)
226 14 : CALL get_qs_env(qs_env, dft_control=dft_control)
227 14 : IF (dft_control%qs_control%ref_embed_subsys) THEN
228 6 : ALLOCATE (rhs_to_save(SIZE(resp_env%rhs)))
229 28 : rhs_to_save = resp_env%rhs
230 2 : CALL set_qs_env(qs_env, rhs=rhs_to_save)
231 : END IF
232 :
233 14 : DEALLOCATE (ipiv)
234 14 : CALL resp_dealloc(resp_env, rep_sys)
235 : CALL cp_print_key_finished_output(output_unit, logger, resp_section, &
236 28 : "PRINT%PROGRAM_RUN_INFO")
237 :
238 : END IF
239 :
240 9901 : CALL timestop(handle)
241 :
242 9901 : 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 14 : 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 14 : IF (ASSOCIATED(resp_env)) CALL resp_dealloc(resp_env, rep_sys)
254 182 : 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 56 : resp_env%box_hi = 0.0_dp
270 56 : resp_env%box_low = 0.0_dp
271 :
272 14 : resp_env%ncons = 0
273 14 : resp_env%ncons_sec = 0
274 14 : resp_env%nres = 0
275 14 : resp_env%nrest_sec = 0
276 14 : resp_env%npoints = 0
277 14 : resp_env%npoints_proc = 0
278 14 : resp_env%auto_vdw_radii_table = use_cambridge_vdw_radii
279 :
280 14 : 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 14 : SUBROUTINE resp_allocate(resp_env, natom, nvar)
289 : TYPE(resp_type), POINTER :: resp_env
290 : INTEGER, INTENT(IN) :: natom, nvar
291 :
292 14 : IF (.NOT. ASSOCIATED(resp_env%matrix)) THEN
293 56 : ALLOCATE (resp_env%matrix(nvar, nvar))
294 : END IF
295 14 : IF (.NOT. ASSOCIATED(resp_env%rhs)) THEN
296 42 : ALLOCATE (resp_env%rhs(nvar))
297 : END IF
298 14 : IF (.NOT. ASSOCIATED(resp_env%sum_vpot)) THEN
299 42 : ALLOCATE (resp_env%sum_vpot(natom))
300 : END IF
301 1258 : resp_env%matrix = 0.0_dp
302 138 : resp_env%rhs = 0.0_dp
303 104 : resp_env%sum_vpot = 0.0_dp
304 :
305 14 : 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 14 : 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 14 : IF (ASSOCIATED(resp_env)) THEN
319 14 : IF (ASSOCIATED(resp_env%matrix)) THEN
320 14 : DEALLOCATE (resp_env%matrix)
321 : END IF
322 14 : IF (ASSOCIATED(resp_env%rhs)) THEN
323 14 : DEALLOCATE (resp_env%rhs)
324 : END IF
325 14 : IF (ASSOCIATED(resp_env%sum_vpot)) THEN
326 14 : DEALLOCATE (resp_env%sum_vpot)
327 : END IF
328 14 : IF (ASSOCIATED(resp_env%fitpoints)) THEN
329 14 : DEALLOCATE (resp_env%fitpoints)
330 : END IF
331 14 : IF (ASSOCIATED(resp_env%rmin_kind)) THEN
332 10 : DEALLOCATE (resp_env%rmin_kind)
333 : END IF
334 14 : IF (ASSOCIATED(resp_env%rmax_kind)) THEN
335 10 : DEALLOCATE (resp_env%rmax_kind)
336 : END IF
337 14 : DEALLOCATE (resp_env)
338 : END IF
339 14 : IF (ASSOCIATED(rep_sys)) THEN
340 8 : DO i = 1, SIZE(rep_sys)
341 4 : DEALLOCATE (rep_sys(i)%p_resp%atom_surf_list)
342 8 : DEALLOCATE (rep_sys(i)%p_resp)
343 : END DO
344 4 : DEALLOCATE (rep_sys)
345 : END IF
346 :
347 14 : 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 14 : 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 14 : INTEGER, DIMENSION(:), POINTER :: atom_list_cons, my_stride
374 : LOGICAL :: explicit
375 : TYPE(section_vals_type), POINTER :: slab_section, sphere_section
376 :
377 14 : CALL timeset(routineN, handle)
378 :
379 14 : NULLIFY (atom_list_cons, my_stride, sphere_section, slab_section)
380 :
381 : ! get the subsections
382 14 : sphere_section => section_vals_get_subs_vals(resp_section, "SPHERE_SAMPLING")
383 14 : slab_section => section_vals_get_subs_vals(resp_section, "SLAB_SAMPLING")
384 14 : cons_section => section_vals_get_subs_vals(resp_section, "CONSTRAINT")
385 14 : 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 14 : l_val=resp_env%itc)
390 14 : 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 14 : l_val=resp_env%rheavies)
394 14 : IF (resp_env%rheavies) THEN
395 : CALL section_vals_val_get(resp_section, "RESTRAIN_HEAVIES_STRENGTH", &
396 14 : r_val=resp_env%rheavies_strength)
397 : END IF
398 14 : CALL section_vals_val_get(resp_section, "STRIDE", i_vals=my_stride)
399 14 : 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 14 : "or 3 values. Correct your input file.")
402 14 : IF (SIZE(my_stride) == 1) THEN
403 48 : DO i = 1, 3
404 48 : resp_env%stride(i) = my_stride(1)
405 : END DO
406 : ELSE
407 16 : resp_env%stride = my_stride(1:3)
408 : END IF
409 14 : 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 14 : l_val=resp_env%use_repeat_method)
414 14 : IF (resp_env%use_repeat_method) THEN
415 4 : resp_env%ncons = resp_env%ncons + 1
416 : ! restrain heavies should be off
417 4 : 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 14 : atomic_kind_set)
424 :
425 : ! get the parameter for periodic/surface systems
426 14 : CALL section_vals_get(slab_section, explicit=explicit, n_repetition=nrep)
427 14 : IF (explicit) THEN
428 4 : IF (resp_env%molecular_sys) THEN
429 : CALL cp_abort(__LOCATION__, &
430 : "You can only use either SPHERE_SAMPLING or SLAB_SAMPLING, but "// &
431 0 : "not both.")
432 : END IF
433 16 : ALLOCATE (rep_sys(nrep))
434 8 : DO i = 1, nrep
435 52 : ALLOCATE (rep_sys(i)%p_resp)
436 4 : 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 4 : i_rep_section=i)
439 : CALL section_vals_val_get(slab_section, "LENGTH", r_val=rep_sys(i)%p_resp%length, &
440 4 : i_rep_section=i)
441 : CALL section_vals_val_get(slab_section, "SURF_DIRECTION", &
442 4 : i_rep_section=i, i_val=rep_sys(i)%p_resp%my_fit)
443 12 : IF (ANY(rep_sys(i)%p_resp%range_surf < 0.0_dp)) THEN
444 0 : CPABORT("Numbers in RANGE in SLAB_SAMPLING cannot be negative.")
445 : END IF
446 4 : IF (rep_sys(i)%p_resp%length <= EPSILON(0.0_dp)) THEN
447 0 : CPABORT("Parameter LENGTH in SLAB_SAMPLING has to be larger than zero.")
448 : END IF
449 : !list of atoms specifying the surface
450 8 : 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 14 : CALL section_vals_get(cons_section, explicit=explicit)
456 14 : IF (explicit) THEN
457 8 : CALL section_vals_get(cons_section, n_repetition=resp_env%ncons_sec)
458 22 : DO i = 1, resp_env%ncons_sec
459 : CALL section_vals_val_get(cons_section, "EQUAL_CHARGES", &
460 14 : l_val=resp_env%equal_charges, explicit=explicit)
461 14 : IF (.NOT. explicit) CYCLE
462 2 : CALL build_atom_list(cons_section, subsys, atom_list_cons, i)
463 : !instead of using EQUAL_CHARGES the constraint sections could be repeated
464 2 : resp_env%ncons = resp_env%ncons + SIZE(atom_list_cons) - 2
465 24 : DEALLOCATE (atom_list_cons)
466 : END DO
467 : END IF
468 14 : CALL section_vals_get(rest_section, explicit=explicit)
469 14 : IF (explicit) THEN
470 6 : CALL section_vals_get(rest_section, n_repetition=resp_env%nrest_sec)
471 : END IF
472 14 : resp_env%ncons = resp_env%ncons + resp_env%ncons_sec
473 14 : resp_env%nres = resp_env%nres + resp_env%nrest_sec
474 :
475 14 : CALL timestop(handle)
476 :
477 56 : 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 14 : 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 14 : 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 14 : 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 14 : nrep_rmin = 0
510 14 : nrep_rmax = 0
511 14 : nkind = SIZE(atomic_kind_set)
512 :
513 14 : has_rmin = .FALSE.
514 14 : has_rmax = .FALSE.
515 :
516 14 : CALL section_vals_get(sphere_section, explicit=explicit)
517 14 : IF (explicit) THEN
518 10 : resp_env%molecular_sys = .TRUE.
519 : CALL section_vals_val_get(sphere_section, "AUTO_VDW_RADII_TABLE", &
520 10 : i_val=resp_env%auto_vdw_radii_table)
521 10 : CALL section_vals_val_get(sphere_section, "AUTO_RMIN_SCALE", r_val=auto_rmin_scale)
522 10 : CALL section_vals_val_get(sphere_section, "AUTO_RMAX_SCALE", r_val=auto_rmax_scale)
523 10 : CALL section_vals_val_get(sphere_section, "RMIN", explicit=has_rmin, r_val=rmin)
524 10 : CALL section_vals_val_get(sphere_section, "RMAX", explicit=has_rmax, r_val=rmax)
525 10 : CALL section_vals_val_get(sphere_section, "RMIN_KIND", n_rep_val=nrep_rmin)
526 10 : CALL section_vals_val_get(sphere_section, "RMAX_KIND", n_rep_val=nrep_rmax)
527 30 : ALLOCATE (resp_env%rmin_kind(nkind))
528 20 : ALLOCATE (resp_env%rmax_kind(nkind))
529 38 : resp_env%rmin_kind = 0.0_dp
530 38 : resp_env%rmax_kind = 0.0_dp
531 30 : ALLOCATE (rmin_is_set(nkind))
532 20 : ALLOCATE (rmax_is_set(nkind))
533 38 : rmin_is_set = .FALSE.
534 38 : rmax_is_set = .FALSE.
535 : ! define rmin_kind and rmax_kind to predefined vdW radii
536 38 : DO ikind = 1, nkind
537 28 : atomic_kind => atomic_kind_set(ikind)
538 : CALL get_atomic_kind(atomic_kind, &
539 : element_symbol=symbol, &
540 : kind_number=kind_number, &
541 28 : z=z)
542 50 : SELECT CASE (resp_env%auto_vdw_radii_table)
543 : CASE (use_cambridge_vdw_radii)
544 22 : CALL get_ptable_info(symbol, vdw_radius=resp_env%rmin_kind(kind_number))
545 22 : rmin_is_set(kind_number) = .TRUE.
546 : CASE (use_uff_vdw_radii)
547 6 : CALL cite_reference(Rappe1992)
548 : CALL get_uff_vdw_radius(z, radius=resp_env%rmin_kind(kind_number), &
549 6 : found=rmin_is_set(kind_number))
550 : CASE DEFAULT
551 0 : CALL get_ptable_info(symbol, vdw_radius=resp_env%rmin_kind(kind_number))
552 28 : rmin_is_set(kind_number) = .TRUE.
553 : END SELECT
554 66 : 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 28 : "angstrom")
557 28 : 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 28 : resp_env%rmin_kind(kind_number)*auto_rmax_scale)
562 28 : 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 10 : IF (has_rmin) THEN
568 24 : resp_env%rmin_kind = rmin
569 24 : rmin_is_set = .TRUE.
570 : END IF
571 10 : IF (has_rmax) THEN
572 24 : resp_env%rmax_kind = rmax
573 24 : 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 10 : DO j = 1, nrep_rmin
578 : CALL section_vals_val_get(sphere_section, "RMIN_KIND", i_rep_val=j, &
579 0 : c_vals=tmpstringlist)
580 10 : DO ikind = 1, nkind
581 0 : atomic_kind => atomic_kind_set(ikind)
582 0 : CALL get_atomic_kind(atomic_kind, element_symbol=symbol, kind_number=kind_number)
583 0 : IF (TRIM(tmpstringlist(2)) == TRIM(symbol)) THEN
584 0 : 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 0 : "angstrom")
588 0 : rmin_is_set(kind_number) = .TRUE.
589 : END IF
590 : END DO
591 : END DO
592 10 : DO j = 1, nrep_rmax
593 : CALL section_vals_val_get(sphere_section, "RMAX_KIND", i_rep_val=j, &
594 0 : c_vals=tmpstringlist)
595 10 : DO ikind = 1, nkind
596 0 : atomic_kind => atomic_kind_set(ikind)
597 0 : CALL get_atomic_kind(atomic_kind, element_symbol=symbol, kind_number=kind_number)
598 0 : IF (TRIM(tmpstringlist(2)) == TRIM(symbol)) THEN
599 0 : 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 0 : "angstrom")
602 0 : 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 10 : n_rmin_missing = 0
608 10 : n_rmax_missing = 0
609 10 : missing_rmin = ""
610 10 : missing_rmax = ""
611 38 : DO ikind = 1, nkind
612 28 : atomic_kind => atomic_kind_set(ikind)
613 : CALL get_atomic_kind(atomic_kind, &
614 : element_symbol=symbol, &
615 28 : kind_number=kind_number)
616 28 : IF (.NOT. rmin_is_set(kind_number)) THEN
617 0 : n_rmin_missing = n_rmin_missing + 1
618 0 : missing_rmin = TRIM(missing_rmin)//" "//TRIM(symbol)//","
619 : END IF
620 66 : IF (.NOT. rmax_is_set(kind_number)) THEN
621 0 : n_rmax_missing = n_rmax_missing + 1
622 0 : missing_rmax = TRIM(missing_rmax)//" "//TRIM(symbol)//","
623 : END IF
624 : END DO
625 10 : 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 0 : "RMIN_KIND in SPHERE_SAMPLING section")
631 : END IF
632 10 : 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 0 : "RMAX_KIND in SPHERE_SAMPLING section")
638 : END IF
639 10 : IF (n_rmin_missing > 0 .OR. &
640 : n_rmax_missing > 0) THEN
641 0 : CPABORT("Insufficient data for RMIN or RMAX")
642 : END IF
643 :
644 10 : CALL get_cell(cell=cell, h=hmat)
645 40 : resp_env%box_hi = (/hmat(1, 1), hmat(2, 2), hmat(3, 3)/)
646 40 : resp_env%box_low = 0.0_dp
647 10 : CALL section_vals_val_get(sphere_section, "X_HI", explicit=explicit)
648 10 : IF (explicit) CALL section_vals_val_get(sphere_section, "X_HI", &
649 0 : r_val=resp_env%box_hi(1))
650 10 : CALL section_vals_val_get(sphere_section, "X_LOW", explicit=explicit)
651 10 : IF (explicit) CALL section_vals_val_get(sphere_section, "X_LOW", &
652 0 : r_val=resp_env%box_low(1))
653 10 : CALL section_vals_val_get(sphere_section, "Y_HI", explicit=explicit)
654 10 : IF (explicit) CALL section_vals_val_get(sphere_section, "Y_HI", &
655 0 : r_val=resp_env%box_hi(2))
656 10 : CALL section_vals_val_get(sphere_section, "Y_LOW", explicit=explicit)
657 10 : IF (explicit) CALL section_vals_val_get(sphere_section, "Y_LOW", &
658 0 : r_val=resp_env%box_low(2))
659 10 : CALL section_vals_val_get(sphere_section, "Z_HI", explicit=explicit)
660 10 : IF (explicit) CALL section_vals_val_get(sphere_section, "Z_HI", &
661 0 : r_val=resp_env%box_hi(3))
662 10 : CALL section_vals_val_get(sphere_section, "Z_LOW", explicit=explicit)
663 10 : IF (explicit) CALL section_vals_val_get(sphere_section, "Z_LOW", &
664 0 : r_val=resp_env%box_low(3))
665 :
666 10 : DEALLOCATE (rmin_is_set)
667 80 : DEALLOCATE (rmax_is_set)
668 : END IF
669 :
670 14 : 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 26 : 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 26 : INTEGER, DIMENSION(:), POINTER :: indexes
693 : LOGICAL :: index_in_range
694 :
695 26 : CALL timeset(routineN, handle)
696 :
697 26 : NULLIFY (indexes)
698 26 : irep = 1
699 26 : IF (PRESENT(rep)) irep = rep
700 :
701 : CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
702 26 : n_rep_val=n_var)
703 26 : num_atom = 0
704 52 : DO i = 1, n_var
705 : CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
706 26 : i_rep_val=i, i_vals=indexes)
707 52 : num_atom = num_atom + SIZE(indexes)
708 : END DO
709 78 : ALLOCATE (atom_list(num_atom))
710 100 : atom_list = 0
711 26 : num_atom = 1
712 52 : DO i = 1, n_var
713 : CALL section_vals_val_get(section, "ATOM_LIST", i_rep_section=irep, &
714 26 : i_rep_val=i, i_vals=indexes)
715 200 : atom_list(num_atom:num_atom + SIZE(indexes) - 1) = indexes(:)
716 52 : num_atom = num_atom + SIZE(indexes)
717 : END DO
718 : !check atom list
719 26 : num_atom = num_atom - 1
720 26 : CALL qs_subsys_get(subsys, nparticle=max_index)
721 26 : CPASSERT(SIZE(atom_list) /= 0)
722 : index_in_range = (MAXVAL(atom_list) <= max_index) &
723 200 : .AND. (MINVAL(atom_list) > 0)
724 0 : CPASSERT(index_in_range)
725 100 : DO i = 1, num_atom
726 236 : DO j = i + 1, num_atom
727 136 : atom_a = atom_list(i)
728 136 : atom_b = atom_list(j)
729 136 : IF (atom_a == atom_b) &
730 74 : CPABORT("There are atoms doubled in atom list for RESP.")
731 : END DO
732 : END DO
733 :
734 26 : CALL timestop(handle)
735 :
736 78 : 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 4 : 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 4 : 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 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dist
770 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat, hmat_inv
771 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
772 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_pw
773 :
774 4 : CALL timeset(routineN, handle)
775 :
776 4 : NULLIFY (particle_set, v_hartree_pw)
777 4 : delta = 1.0E-13_dp
778 :
779 4 : CALL get_cell(cell=cell, h=hmat, h_inv=hmat_inv)
780 :
781 4 : IF (.NOT. cell%orthorhombic) THEN
782 : CALL cp_abort(__LOCATION__, &
783 : "Nonperiodic solution for RESP charges only"// &
784 0 : " implemented for orthorhombic cells!")
785 : END IF
786 4 : IF (.NOT. resp_env%molecular_sys) THEN
787 : CALL cp_abort(__LOCATION__, &
788 : "Nonperiodic solution for RESP charges (i.e. nonperiodic"// &
789 0 : " Poisson solver) can only be used with section SPHERE_SAMPLING")
790 : END IF
791 4 : IF (resp_env%use_repeat_method) THEN
792 : CALL cp_abort(__LOCATION__, &
793 0 : "REPEAT method only reasonable for periodic RESP fitting")
794 : END IF
795 4 : CALL get_qs_env(qs_env, particle_set=particle_set, v_hartree_rspace=v_hartree_pw)
796 :
797 40 : bo = v_hartree_pw%pw_grid%bounds_local
798 40 : gbo = v_hartree_pw%pw_grid%bounds
799 16 : np = v_hartree_pw%pw_grid%npts
800 52 : dh = v_hartree_pw%pw_grid%dh
801 4 : dvol = v_hartree_pw%pw_grid%dvol
802 4 : nkind = SIZE(atomic_kind_set)
803 :
804 12 : ALLOCATE (dist(natom))
805 12 : ALLOCATE (not_in_range(natom, 2))
806 :
807 : ! store fitting points to calculate the RMS and RRMS later
808 4 : IF (.NOT. ASSOCIATED(resp_env%fitpoints)) THEN
809 4 : now = 1000
810 4 : ALLOCATE (resp_env%fitpoints(3, now))
811 : ELSE
812 0 : now = SIZE(resp_env%fitpoints, 2)
813 : END IF
814 :
815 184 : DO jz = bo(1, 3), bo(2, 3)
816 8284 : DO jy = bo(1, 2), bo(2, 2)
817 190530 : DO jx = bo(1, 1), bo(2, 1)
818 182250 : IF (.NOT. (MODULO(jz, resp_env%stride(3)) == 0)) CYCLE
819 60750 : IF (.NOT. (MODULO(jy, resp_env%stride(2)) == 0)) CYCLE
820 20250 : 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 6750 : l = jx - gbo(1, 1)
823 6750 : k = jy - gbo(1, 2)
824 6750 : p = jz - gbo(1, 3)
825 6750 : r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
826 6750 : r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
827 6750 : r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
828 6750 : IF (r(3) < resp_env%box_low(3) .OR. r(3) > resp_env%box_hi(3)) CYCLE
829 6750 : IF (r(2) < resp_env%box_low(2) .OR. r(2) > resp_env%box_hi(2)) CYCLE
830 6750 : 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 101250 : not_in_range = .FALSE.
833 47250 : DO i = 1, natom
834 162000 : vec = r - particles%els(i)%r
835 40500 : vec_pbc(1) = vec(1) - hmat(1, 1)*ANINT(hmat_inv(1, 1)*vec(1))
836 40500 : vec_pbc(2) = vec(2) - hmat(2, 2)*ANINT(hmat_inv(2, 2)*vec(2))
837 40500 : vec_pbc(3) = vec(3) - hmat(3, 3)*ANINT(hmat_inv(3, 3)*vec(3))
838 162000 : dist(i) = SQRT(SUM(vec_pbc**2))
839 : CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, &
840 40500 : kind_number=kind_number)
841 101250 : DO ikind = 1, nkind
842 101250 : IF (ikind == kind_number) THEN
843 40500 : rmin = resp_env%rmin_kind(ikind)
844 40500 : rmax = resp_env%rmax_kind(ikind)
845 40500 : EXIT
846 : END IF
847 : END DO
848 40500 : IF (dist(i) < rmin + delta) not_in_range(i, 1) = .TRUE.
849 87750 : 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 85830 : IF (ANY(not_in_range(:, 1)) .OR. ALL(not_in_range(:, 2))) CYCLE
854 72 : resp_env%npoints_proc = resp_env%npoints_proc + 1
855 72 : IF (resp_env%npoints_proc > now) THEN
856 0 : now = 2*now
857 0 : CALL reallocate(resp_env%fitpoints, 1, 3, 1, now)
858 : END IF
859 72 : resp_env%fitpoints(1, resp_env%npoints_proc) = jx
860 72 : resp_env%fitpoints(2, resp_env%npoints_proc) = jy
861 72 : 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 72 : 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 0 : vj = -v_hartree_pw%array(jx, jy, jz)/dvol + qs_env%ks_qmmm_env%v_qmmm_rspace%array(jx, jy, jz)
866 : ELSE
867 72 : vj = -v_hartree_pw%array(jx, jy, jz)/dvol
868 : END IF
869 504 : dist(:) = 1.0_dp/dist(:)
870 :
871 8604 : DO i = 1, natom
872 3024 : DO m = 1, natom
873 3024 : matrix(m, i) = matrix(m, i) + 2.0_dp*dist(i)*dist(m)
874 : END DO
875 182682 : rhs(i) = rhs(i) + 2.0_dp*vj*dist(i)
876 : END DO
877 : END DO
878 : END DO
879 : END DO
880 :
881 4 : resp_env%npoints = resp_env%npoints_proc
882 4 : CALL v_hartree_pw%pw_grid%para%group%sum(resp_env%npoints)
883 724 : CALL v_hartree_pw%pw_grid%para%group%sum(matrix)
884 76 : CALL v_hartree_pw%pw_grid%para%group%sum(rhs)
885 : !weighted sum
886 364 : matrix = matrix/resp_env%npoints
887 40 : rhs = rhs/resp_env%npoints
888 :
889 4 : DEALLOCATE (dist)
890 4 : DEALLOCATE (not_in_range)
891 :
892 4 : CALL timestop(handle)
893 :
894 4 : 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 10 : 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 10 : 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 10 : CALL timeset(routineN, handle)
929 :
930 10 : NULLIFY (pw_env, para_env, auxbas_pw_pool, poisson_env)
931 :
932 10 : CALL get_cell(cell=cell, periodic=periodic)
933 :
934 40 : IF (.NOT. ALL(periodic /= 0)) THEN
935 : CALL cp_abort(__LOCATION__, &
936 : "Periodic solution for RESP (with periodic Poisson solver)"// &
937 0 : " can only be obtained with a cell that has XYZ periodicity")
938 : END IF
939 :
940 10 : 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 10 : poisson_env=poisson_env)
944 10 : CALL auxbas_pw_pool%create_pw(rho_ga)
945 10 : CALL auxbas_pw_pool%create_pw(va_gspace)
946 10 : 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 10 : cell=cell)
951 40 : ALLOCATE (vpot(resp_env%npoints_proc, natom))
952 10 : normalize_factor = SQRT((resp_env%eta/pi)**3)
953 :
954 76 : DO i = 1, natom
955 : !collocate gaussian for each atom
956 66 : CALL pw_zero(rho_ga)
957 66 : 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 66 : CALL pw_zero(va_gspace)
960 66 : CALL pw_poisson_solve(poisson_env, rho_ga, vhartree=va_gspace)
961 66 : CALL pw_zero(va_rspace)
962 66 : CALL pw_transfer(va_gspace, va_rspace)
963 66 : CALL pw_scale(va_rspace, normalize_factor)
964 10659 : DO ip = 1, resp_env%npoints_proc
965 10583 : jx = resp_env%fitpoints(1, ip)
966 10583 : jy = resp_env%fitpoints(2, ip)
967 10583 : jz = resp_env%fitpoints(3, ip)
968 10649 : vpot(ip, i) = va_rspace%array(jx, jy, jz)
969 : END DO
970 : END DO
971 :
972 10 : CALL va_gspace%release()
973 10 : CALL va_rspace%release()
974 10 : CALL rho_ga%release()
975 :
976 76 : DO i = 1, natom
977 516 : DO j = 1, natom
978 : ! calculate matrix
979 55897 : 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 76 : CALL calculate_rhs(qs_env, resp_env, resp_env%rhs(i), vpot(:, i))
983 : END DO
984 :
985 1778 : CALL para_env%sum(resp_env%matrix)
986 186 : CALL para_env%sum(resp_env%rhs)
987 : !weighted sum
988 894 : resp_env%matrix = resp_env%matrix/resp_env%npoints
989 98 : resp_env%rhs = resp_env%rhs/resp_env%npoints
990 :
991 : ! REPEAT stuff
992 10 : IF (resp_env%use_repeat_method) THEN
993 : ! sum over selected points of single Gaussian potential vpot
994 32 : DO i = 1, natom
995 32 : resp_env%sum_vpot(i) = 2.0_dp*accurate_sum(vpot(:, i))/resp_env%npoints
996 : END DO
997 60 : CALL para_env%sum(resp_env%sum_vpot)
998 4 : CALL para_env%sum(resp_env%sum_vhartree)
999 4 : resp_env%sum_vhartree = 2.0_dp*resp_env%sum_vhartree/resp_env%npoints
1000 : END IF
1001 :
1002 10 : DEALLOCATE (vpot)
1003 10 : CALL timestop(handle)
1004 :
1005 10 : 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 10 : 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 10 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: not_in_range
1030 : REAL(KIND=dp) :: delta, dh(3, 3), r(3), rmax, rmin, &
1031 : vec_pbc(3)
1032 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dist
1033 10 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1034 : TYPE(mp_para_env_type), POINTER :: para_env
1035 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1036 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_pw
1037 :
1038 10 : CALL timeset(routineN, handle)
1039 :
1040 10 : NULLIFY (atomic_kind_set, v_hartree_pw, para_env, particle_set)
1041 10 : 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 10 : v_hartree_rspace=v_hartree_pw)
1048 :
1049 100 : bo = v_hartree_pw%pw_grid%bounds_local
1050 100 : gbo = v_hartree_pw%pw_grid%bounds
1051 130 : dh = v_hartree_pw%pw_grid%dh
1052 10 : natom = SIZE(particles%els)
1053 10 : nkind = SIZE(atomic_kind_set)
1054 :
1055 10 : IF (.NOT. ASSOCIATED(resp_env%fitpoints)) THEN
1056 10 : now = 1000
1057 10 : ALLOCATE (resp_env%fitpoints(3, now))
1058 : ELSE
1059 0 : now = SIZE(resp_env%fitpoints, 2)
1060 : END IF
1061 :
1062 30 : ALLOCATE (dist(natom))
1063 30 : ALLOCATE (not_in_range(natom, 2))
1064 :
1065 : !every proc gets another bo, grid is distributed
1066 350 : DO jz = bo(1, 3), bo(2, 3)
1067 340 : IF (.NOT. (MODULO(jz, resp_env%stride(3)) == 0)) CYCLE
1068 4338 : DO jy = bo(1, 2), bo(2, 2)
1069 4204 : IF (.NOT. (MODULO(jy, resp_env%stride(2)) == 0)) CYCLE
1070 31554 : DO jx = bo(1, 1), bo(2, 1)
1071 29642 : 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 11246 : l = jx - gbo(1, 1)
1074 11246 : k = jy - gbo(1, 2)
1075 11246 : p = jz - gbo(1, 3)
1076 11246 : r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
1077 11246 : r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
1078 11246 : r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
1079 11246 : IF (resp_env%molecular_sys) THEN
1080 154498 : not_in_range = .FALSE.
1081 71826 : DO m = 1, natom
1082 60980 : vec_pbc = pbc(r, particles%els(m)%r, cell)
1083 243920 : dist(m) = SQRT(SUM(vec_pbc**2))
1084 : CALL get_atomic_kind(atomic_kind=particle_set(m)%atomic_kind, &
1085 60980 : kind_number=kind_number)
1086 138114 : DO ikind = 1, nkind
1087 138114 : IF (ikind == kind_number) THEN
1088 60980 : rmin = resp_env%rmin_kind(ikind)
1089 60980 : rmax = resp_env%rmax_kind(ikind)
1090 60980 : EXIT
1091 : END IF
1092 : END DO
1093 60980 : IF (dist(m) < rmin + delta) not_in_range(m, 1) = .TRUE.
1094 132806 : IF (dist(m) > rmax - delta) not_in_range(m, 2) = .TRUE.
1095 : END DO
1096 121136 : IF (ANY(not_in_range(:, 1)) .OR. ALL(not_in_range(:, 2))) CYCLE
1097 : ELSE
1098 752 : DO i = 1, SIZE(rep_sys)
1099 3348 : DO m = 1, SIZE(rep_sys(i)%p_resp%atom_surf_list)
1100 2996 : in_z = 0
1101 2996 : in_y = 0
1102 2996 : in_x = 0
1103 2996 : iatom = rep_sys(i)%p_resp%atom_surf_list(m)
1104 5992 : SELECT CASE (rep_sys(i)%p_resp%my_fit)
1105 : CASE (do_resp_x_dir, do_resp_y_dir, do_resp_z_dir)
1106 2996 : 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 2996 : vec_pbc = pbc(r, particles%els(iatom)%r, cell)
1109 : END SELECT
1110 2996 : 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 0 : IF (ABS(vec_pbc(3)) < rep_sys(i)%p_resp%length - delta) in_z = 1
1114 0 : IF (ABS(vec_pbc(2)) < rep_sys(i)%p_resp%length - delta) in_y = 1
1115 0 : IF (vec_pbc(1) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
1116 0 : 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 0 : IF (ABS(vec_pbc(3)) < rep_sys(i)%p_resp%length - delta) in_z = 1
1119 0 : IF (vec_pbc(2) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
1120 0 : vec_pbc(2) < rep_sys(i)%p_resp%range_surf(2) - delta) in_y = 1
1121 0 : 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 2996 : IF (vec_pbc(3) > rep_sys(i)%p_resp%range_surf(1) + delta .AND. &
1124 196 : vec_pbc(3) < rep_sys(i)%p_resp%range_surf(2) - delta) in_z = 1
1125 2996 : IF (ABS(vec_pbc(2)) < rep_sys(i)%p_resp%length - delta) in_y = 1
1126 5992 : IF (ABS(vec_pbc(1)) < rep_sys(i)%p_resp%length - delta) in_x = 1
1127 : END SELECT
1128 3348 : IF (in_z*in_y*in_x == 1) EXIT
1129 : END DO
1130 752 : IF (in_z*in_y*in_x == 1) EXIT
1131 : END DO
1132 400 : IF (in_z*in_y*in_x == 0) CYCLE
1133 : END IF
1134 2044 : resp_env%npoints_proc = resp_env%npoints_proc + 1
1135 2044 : IF (resp_env%npoints_proc > now) THEN
1136 1 : now = 2*now
1137 1 : CALL reallocate(resp_env%fitpoints, 1, 3, 1, now)
1138 : END IF
1139 2044 : resp_env%fitpoints(1, resp_env%npoints_proc) = jx
1140 2044 : resp_env%fitpoints(2, resp_env%npoints_proc) = jy
1141 33846 : resp_env%fitpoints(3, resp_env%npoints_proc) = jz
1142 : END DO
1143 : END DO
1144 : END DO
1145 :
1146 10 : resp_env%npoints = resp_env%npoints_proc
1147 10 : CALL para_env%sum(resp_env%npoints)
1148 :
1149 10 : DEALLOCATE (dist)
1150 10 : DEALLOCATE (not_in_range)
1151 :
1152 10 : CALL timestop(handle)
1153 :
1154 10 : 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 66 : 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 66 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vhartree
1175 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_pw
1176 :
1177 66 : CALL timeset(routineN, handle)
1178 :
1179 66 : NULLIFY (v_hartree_pw)
1180 66 : CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_pw)
1181 66 : dvol = v_hartree_pw%pw_grid%dvol
1182 198 : ALLOCATE (vhartree(resp_env%npoints_proc))
1183 10649 : 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 10649 : DO ip = 1, resp_env%npoints_proc
1188 10583 : jx = resp_env%fitpoints(1, ip)
1189 10583 : jy = resp_env%fitpoints(2, ip)
1190 10583 : jz = resp_env%fitpoints(3, ip)
1191 10583 : vhartree(ip) = -v_hartree_pw%array(jx, jy, jz)/dvol
1192 10583 : IF (qs_env%qmmm) THEN
1193 : !taking into account that v_qmmm has also opposite sign
1194 0 : vhartree(ip) = vhartree(ip) + qs_env%ks_qmmm_env%v_qmmm_rspace%array(jx, jy, jz)
1195 : END IF
1196 10649 : rhs = rhs + 2.0_dp*vhartree(ip)*vpot(ip)
1197 : END DO
1198 :
1199 66 : IF (resp_env%use_repeat_method) THEN
1200 28 : resp_env%sum_vhartree = accurate_sum(vhartree)
1201 : END IF
1202 :
1203 66 : DEALLOCATE (vhartree)
1204 :
1205 66 : CALL timestop(handle)
1206 :
1207 132 : 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 28 : 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 14 : INTEGER, DIMENSION(:), POINTER :: tmp_npoints, tmp_size
1228 14 : 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 98 : TYPE(mp_request_type), DIMENSION(6) :: req
1233 14 : 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 14 : CALL timeset(routineN, handle)
1238 :
1239 14 : NULLIFY (para_env, input, logger, resp_section, print_key, particle_set, tmp_size, &
1240 14 : tmp_points, tmp_npoints, v_hartree_pw)
1241 :
1242 : CALL get_qs_env(qs_env, input=input, para_env=para_env, &
1243 14 : particle_set=particle_set, v_hartree_rspace=v_hartree_pw)
1244 14 : conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
1245 140 : gbo = v_hartree_pw%pw_grid%bounds
1246 182 : dh = v_hartree_pw%pw_grid%dh
1247 14 : nobjects = SIZE(particle_set) + resp_env%npoints
1248 :
1249 14 : resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1250 14 : print_key => section_vals_get_subs_vals(resp_section, "PRINT%COORD_FIT_POINTS")
1251 14 : 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 14 : file_form="FORMATTED")
1258 :
1259 14 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1260 : resp_section, "PRINT%COORD_FIT_POINTS"), &
1261 : cp_p_file)) THEN
1262 2 : IF (output_unit > 0) THEN
1263 : filename = cp_print_key_generate_filename(logger, &
1264 : print_key, extension=".xyz", &
1265 1 : my_local=.FALSE.)
1266 1 : WRITE (unit=output_unit, FMT="(I12,/)") nobjects
1267 7 : DO iatom = 1, SIZE(particle_set)
1268 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
1269 6 : element_symbol=element_symbol)
1270 6 : WRITE (UNIT=output_unit, FMT="(A,1X,3F10.5)") element_symbol, &
1271 31 : 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 101 : DO ip = 1, resp_env%npoints_proc
1275 100 : jx = resp_env%fitpoints(1, ip)
1276 100 : jy = resp_env%fitpoints(2, ip)
1277 100 : jz = resp_env%fitpoints(3, ip)
1278 100 : l = jx - gbo(1, 1)
1279 100 : k = jy - gbo(1, 2)
1280 100 : p = jz - gbo(1, 3)
1281 100 : r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
1282 100 : r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
1283 100 : r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
1284 400 : r(:) = r(:)*conv
1285 101 : WRITE (UNIT=output_unit, FMT="(A,2X,3F10.5)") "X", r(1), r(2), r(3)
1286 : END DO
1287 : END IF
1288 :
1289 2 : ALLOCATE (tmp_size(1))
1290 2 : ALLOCATE (tmp_npoints(1))
1291 :
1292 : !sending data of all other procs to proc which makes the output (proc 0)
1293 2 : IF (output_unit > 0) THEN
1294 1 : my_pos = para_env%mepos
1295 3 : DO i = 1, para_env%num_pe
1296 2 : IF (my_pos == i - 1) CYCLE
1297 : CALL para_env%irecv(msgout=tmp_size, source=i - 1, &
1298 1 : request=req(1))
1299 1 : CALL req(1)%wait()
1300 3 : ALLOCATE (tmp_points(3, tmp_size(1)))
1301 : CALL para_env%irecv(msgout=tmp_points, source=i - 1, &
1302 1 : request=req(3))
1303 1 : CALL req(3)%wait()
1304 : CALL para_env%irecv(msgout=tmp_npoints, source=i - 1, &
1305 1 : request=req(5))
1306 1 : CALL req(5)%wait()
1307 84 : DO ip = 1, tmp_npoints(1)
1308 83 : jx = tmp_points(1, ip)
1309 83 : jy = tmp_points(2, ip)
1310 83 : jz = tmp_points(3, ip)
1311 83 : l = jx - gbo(1, 1)
1312 83 : k = jy - gbo(1, 2)
1313 83 : p = jz - gbo(1, 3)
1314 83 : r(3) = p*dh(3, 3) + k*dh(3, 2) + l*dh(3, 1)
1315 83 : r(2) = p*dh(2, 3) + k*dh(2, 2) + l*dh(2, 1)
1316 83 : r(1) = p*dh(1, 3) + k*dh(1, 2) + l*dh(1, 1)
1317 332 : r(:) = r(:)*conv
1318 84 : WRITE (UNIT=output_unit, FMT="(A,2X,3F10.5)") "X", r(1), r(2), r(3)
1319 : END DO
1320 3 : DEALLOCATE (tmp_points)
1321 : END DO
1322 : ELSE
1323 1 : 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 1 : request=req(2))
1327 1 : CALL req(2)%wait()
1328 : CALL para_env%isend(msgin=resp_env%fitpoints, dest=para_env%source, &
1329 1 : request=req(4))
1330 1 : CALL req(4)%wait()
1331 1 : tmp_npoints(1) = resp_env%npoints_proc
1332 : CALL para_env%isend(msgin=tmp_npoints, dest=para_env%source, &
1333 1 : request=req(6))
1334 1 : CALL req(6)%wait()
1335 : END IF
1336 :
1337 2 : DEALLOCATE (tmp_size)
1338 2 : DEALLOCATE (tmp_npoints)
1339 : END IF
1340 :
1341 : CALL cp_print_key_finished_output(output_unit, logger, resp_section, &
1342 14 : "PRINT%COORD_FIT_POINTS")
1343 :
1344 14 : CALL timestop(handle)
1345 :
1346 28 : 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 14 : 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 14 : INTEGER, DIMENSION(:), POINTER :: atom_list_cons, atom_list_res
1373 : LOGICAL :: explicit_coeff
1374 : REAL(KIND=dp) :: my_atom_coef(2), strength, TARGET
1375 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: atom_coef
1376 : TYPE(dft_control_type), POINTER :: dft_control
1377 :
1378 14 : CALL timeset(routineN, handle)
1379 :
1380 14 : NULLIFY (atom_coef, atom_list_res, atom_list_cons, dft_control)
1381 :
1382 14 : CALL get_qs_env(qs_env, dft_control=dft_control)
1383 :
1384 : !*** add the restraints
1385 20 : DO i = 1, resp_env%nrest_sec
1386 6 : CALL section_vals_val_get(rest_section, "TARGET", i_rep_section=i, r_val=TARGET)
1387 6 : CALL section_vals_val_get(rest_section, "STRENGTH", i_rep_section=i, r_val=strength)
1388 6 : CALL build_atom_list(rest_section, subsys, atom_list_res, i)
1389 6 : CALL section_vals_val_get(rest_section, "ATOM_COEF", i_rep_section=i, explicit=explicit_coeff)
1390 6 : IF (explicit_coeff) THEN
1391 6 : CALL section_vals_val_get(rest_section, "ATOM_COEF", i_rep_section=i, r_vals=atom_coef)
1392 6 : CPASSERT(SIZE(atom_list_res) == SIZE(atom_coef))
1393 : END IF
1394 12 : DO m = 1, SIZE(atom_list_res)
1395 12 : IF (explicit_coeff) THEN
1396 12 : 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 12 : 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 6 : 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 0 : 2.0_dp*strength
1407 : resp_env%rhs(atom_list_res(m)) = resp_env%rhs(atom_list_res(m)) + &
1408 0 : 2.0_dp*TARGET*strength
1409 : END IF
1410 : END DO
1411 32 : DEALLOCATE (atom_list_res)
1412 : END DO
1413 :
1414 : ! if heavies are restrained to zero, add these as well
1415 14 : IF (resp_env%rheavies) THEN
1416 72 : DO i = 1, natom
1417 62 : CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, z=z)
1418 72 : IF (z .NE. 1) THEN
1419 30 : 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 14 : ncons_v = 0
1426 14 : ncons_v = ncons_v + natom
1427 :
1428 : ! REPEAT charges: treat the offset like a constraint
1429 14 : IF (resp_env%use_repeat_method) THEN
1430 4 : ncons_v = ncons_v + 1
1431 32 : resp_env%matrix(1:natom, ncons_v) = resp_env%sum_vpot(1:natom)
1432 32 : resp_env%matrix(ncons_v, 1:natom) = resp_env%sum_vpot(1:natom)
1433 4 : resp_env%matrix(ncons_v, ncons_v) = 2.0_dp
1434 4 : resp_env%rhs(ncons_v) = resp_env%sum_vhartree
1435 : END IF
1436 :
1437 : ! total charge constraint
1438 14 : IF (resp_env%itc) THEN
1439 14 : ncons_v = ncons_v + 1
1440 104 : resp_env%matrix(1:natom, ncons_v) = 1.0_dp
1441 104 : resp_env%matrix(ncons_v, 1:natom) = 1.0_dp
1442 14 : resp_env%rhs(ncons_v) = dft_control%charge
1443 : END IF
1444 :
1445 : ! explicit constraints
1446 28 : DO i = 1, resp_env%ncons_sec
1447 14 : CALL build_atom_list(cons_section, subsys, atom_list_cons, i)
1448 14 : IF (.NOT. resp_env%equal_charges) THEN
1449 12 : ncons_v = ncons_v + 1
1450 12 : CALL section_vals_val_get(cons_section, "ATOM_COEF", i_rep_section=i, r_vals=atom_coef)
1451 12 : CALL section_vals_val_get(cons_section, "TARGET", i_rep_section=i, r_val=TARGET)
1452 12 : CPASSERT(SIZE(atom_list_cons) == SIZE(atom_coef))
1453 36 : DO m = 1, SIZE(atom_list_cons)
1454 24 : resp_env%matrix(atom_list_cons(m), ncons_v) = atom_coef(m)
1455 36 : resp_env%matrix(ncons_v, atom_list_cons(m)) = atom_coef(m)
1456 : END DO
1457 12 : resp_env%rhs(ncons_v) = TARGET
1458 : ELSE
1459 2 : my_atom_coef(1) = 1.0_dp
1460 2 : my_atom_coef(2) = -1.0_dp
1461 6 : DO k = 2, SIZE(atom_list_cons)
1462 4 : ncons_v = ncons_v + 1
1463 4 : resp_env%matrix(atom_list_cons(1), ncons_v) = my_atom_coef(1)
1464 4 : resp_env%matrix(ncons_v, atom_list_cons(1)) = my_atom_coef(1)
1465 4 : resp_env%matrix(atom_list_cons(k), ncons_v) = my_atom_coef(2)
1466 4 : resp_env%matrix(ncons_v, atom_list_cons(k)) = my_atom_coef(2)
1467 6 : resp_env%rhs(ncons_v) = 0.0_dp
1468 : END DO
1469 : END IF
1470 28 : DEALLOCATE (atom_list_cons)
1471 : END DO
1472 14 : CALL timestop(handle)
1473 :
1474 14 : 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 14 : 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 14 : 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 14 : CALL timeset(routineN, handle)
1501 14 : NULLIFY (logger, input, resp_section)
1502 :
1503 : CALL get_qs_env(qs_env, &
1504 : input=input, &
1505 14 : atomic_kind_set=atomic_kind_set)
1506 14 : resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1507 14 : logger => cp_get_default_logger()
1508 : output_unit = cp_print_key_unit_nr(logger, resp_section, "PRINT%PROGRAM_RUN_INFO", &
1509 14 : extension=".resp")
1510 14 : nkinds = SIZE(atomic_kind_set)
1511 :
1512 14 : conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
1513 14 : IF (.NOT. my_per == use_perd_none) THEN
1514 10 : eta_conv = cp_unit_from_cp2k(resp_env%eta, "angstrom", power=-2)
1515 : END IF
1516 :
1517 14 : IF (output_unit > 0) THEN
1518 7 : WRITE (output_unit, '(/,1X,A,/)') "STARTING RESP FIT"
1519 7 : IF (resp_env%use_repeat_method) THEN
1520 : WRITE (output_unit, '(T3,A)') &
1521 2 : "Fit the variance of the potential (REPEAT method)."
1522 : END IF
1523 7 : IF (.NOT. resp_env%equal_charges) THEN
1524 6 : WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons_sec
1525 : ELSE
1526 1 : IF (resp_env%itc) THEN
1527 1 : WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons - 1
1528 : ELSE
1529 0 : WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit constraints: ", resp_env%ncons
1530 : END IF
1531 : END IF
1532 7 : WRITE (output_unit, '(T3,A,T75,I6)') "Number of explicit restraints: ", resp_env%nrest_sec
1533 7 : WRITE (output_unit, '(T3,A,T80,A)') "Constrain total charge ", MERGE("T", "F", resp_env%itc)
1534 9 : WRITE (output_unit, '(T3,A,T80,A)') "Restrain heavy atoms ", MERGE("T", "F", resp_env%rheavies)
1535 7 : IF (resp_env%rheavies) THEN
1536 5 : WRITE (output_unit, '(T3,A,T71,F10.6)') "Heavy atom restraint strength: ", &
1537 10 : resp_env%rheavies_strength
1538 : END IF
1539 7 : WRITE (output_unit, '(T3,A,T66,3I5)') "Stride: ", resp_env%stride
1540 7 : IF (resp_env%molecular_sys) THEN
1541 : WRITE (output_unit, '(T3,A)') &
1542 5 : "------------------------------------------------------------------------------"
1543 5 : WRITE (output_unit, '(T3,A)') "Using sphere sampling"
1544 : WRITE (output_unit, '(T3,A,T46,A,T66,A)') &
1545 5 : "Element", "RMIN [angstrom]", "RMAX [angstrom]"
1546 19 : DO ikind = 1, nkinds
1547 : CALL get_atomic_kind(atomic_kind=atomic_kind_set(ikind), &
1548 : kind_number=kind_number, &
1549 14 : element_symbol=symbol)
1550 : WRITE (output_unit, '(T3,A,T51,F10.5,T71,F10.5)') &
1551 14 : symbol, &
1552 14 : resp_env%rmin_kind(kind_number)*conv, &
1553 33 : resp_env%rmax_kind(kind_number)*conv
1554 : END DO
1555 5 : IF (my_per == use_perd_none) THEN
1556 8 : WRITE (output_unit, '(T3,A,T51,3F10.5)') "Box min [angstrom]: ", resp_env%box_low(1:3)*conv
1557 8 : 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 5 : "------------------------------------------------------------------------------"
1561 : ELSE
1562 : WRITE (output_unit, '(T3,A)') &
1563 2 : "------------------------------------------------------------------------------"
1564 2 : WRITE (output_unit, '(T3,A)') "Using slab sampling"
1565 2 : WRITE (output_unit, '(2X,A,F10.5)') "Index of atoms defining the surface: "
1566 4 : DO i = 1, SIZE(rep_sys)
1567 18 : IF (i > 1 .AND. ALL(rep_sys(i)%p_resp%atom_surf_list == rep_sys(1)%p_resp%atom_surf_list)) EXIT
1568 20 : WRITE (output_unit, '(7X,10I6)') rep_sys(i)%p_resp%atom_surf_list
1569 : END DO
1570 4 : DO i = 1, SIZE(rep_sys)
1571 6 : 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 2 : "Range for sampling above the surface [angstrom]:", &
1574 10 : rep_sys(i)%p_resp%range_surf(1:2)*conv
1575 : END DO
1576 4 : DO i = 1, SIZE(rep_sys)
1577 2 : 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 4 : " surface atom [angstrom]: ", rep_sys(i)%p_resp%length*conv
1580 : END DO
1581 : WRITE (output_unit, '(T3,A)') &
1582 2 : "------------------------------------------------------------------------------"
1583 : END IF
1584 7 : IF (.NOT. my_per == use_perd_none) THEN
1585 : WRITE (output_unit, '(T3,A,T71,F10.5)') "Width of Gaussian charge"// &
1586 5 : " distribution [angstrom^-2]: ", eta_conv
1587 : END IF
1588 7 : CALL m_flush(output_unit)
1589 : END IF
1590 : CALL cp_print_key_finished_output(output_unit, logger, resp_section, &
1591 14 : "PRINT%PROGRAM_RUN_INFO")
1592 :
1593 14 : CALL timestop(handle)
1594 :
1595 14 : 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 14 : 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 14 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1616 14 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1617 : TYPE(section_vals_type), POINTER :: input, print_key, resp_section
1618 :
1619 14 : CALL timeset(routineN, handle)
1620 :
1621 14 : 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 14 : qs_kind_set=qs_kind_set)
1625 :
1626 14 : resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1627 : print_key => section_vals_get_subs_vals(resp_section, &
1628 14 : "PRINT%RESP_CHARGES_TO_FILE")
1629 14 : logger => cp_get_default_logger()
1630 :
1631 14 : 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 0 : file_form="FORMATTED")
1640 0 : IF (output_file > 0) THEN
1641 : filename = cp_print_key_generate_filename(logger, &
1642 : print_key, extension=".resp", &
1643 0 : my_local=.FALSE.)
1644 : CALL print_atomic_charges(particle_set, qs_kind_set, output_file, title="RESP charges:", &
1645 0 : atomic_charges=resp_env%rhs(1:natom))
1646 0 : 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 0 : "PRINT%RESP_CHARGES_TO_FILE")
1651 : ELSE
1652 : CALL print_atomic_charges(particle_set, qs_kind_set, output_runinfo, title="RESP charges:", &
1653 14 : atomic_charges=resp_env%rhs(1:natom))
1654 : END IF
1655 :
1656 14 : CALL timestop(handle)
1657 :
1658 14 : 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 14 : 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 14 : CALL timeset(routineN, handle)
1693 :
1694 14 : NULLIFY (auxbas_pw_pool, logger, pw_env, poisson_env, input, print_key, &
1695 14 : 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 14 : v_hartree_rspace=v_hartree_rspace)
1701 14 : normalize_factor = SQRT((resp_env%eta/pi)**3)
1702 14 : resp_section => section_vals_get_subs_vals(input, "PROPERTIES%RESP")
1703 : print_key => section_vals_get_subs_vals(resp_section, &
1704 14 : "PRINT%V_RESP_CUBE")
1705 14 : 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 14 : poisson_env=poisson_env)
1710 :
1711 14 : CALL auxbas_pw_pool%create_pw(rho_resp)
1712 14 : CALL auxbas_pw_pool%create_pw(v_resp_gspace)
1713 14 : CALL auxbas_pw_pool%create_pw(v_resp_rspace)
1714 :
1715 14 : CALL pw_zero(rho_resp)
1716 : CALL calculate_rho_resp_all(rho_resp, resp_env%rhs, natom, &
1717 14 : resp_env%eta, qs_env)
1718 14 : CALL pw_zero(v_resp_gspace)
1719 : CALL pw_poisson_solve(poisson_env, rho_resp, &
1720 14 : vhartree=v_resp_gspace)
1721 14 : CALL pw_zero(v_resp_rspace)
1722 14 : CALL pw_transfer(v_resp_gspace, v_resp_rspace)
1723 14 : dvol = v_resp_rspace%pw_grid%dvol
1724 14 : CALL pw_scale(v_resp_rspace, dvol)
1725 14 : 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 14 : IF (resp_env%use_repeat_method) THEN
1729 101437 : v_resp_rspace%array(:, :, :) = v_resp_rspace%array(:, :, :) - resp_env%offset*dvol
1730 : END IF
1731 14 : CALL v_resp_gspace%release()
1732 14 : CALL rho_resp%release()
1733 :
1734 : !***now print the v_resp_rspace%pw to a cube file if requested
1735 14 : IF (BTEST(cp_print_key_should_output(logger%iter_info, resp_section, &
1736 : "PRINT%V_RESP_CUBE"), cp_p_file)) THEN
1737 2 : CALL auxbas_pw_pool%create_pw(aux_r)
1738 2 : append_cube = section_get_lval(resp_section, "PRINT%V_RESP_CUBE%APPEND")
1739 2 : my_pos_cube = "REWIND"
1740 2 : IF (append_cube) THEN
1741 0 : my_pos_cube = "APPEND"
1742 : END IF
1743 2 : 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 2 : mpi_io=mpi_io)
1749 2 : udvol = 1.0_dp/dvol
1750 2 : CALL pw_copy(v_resp_rspace, aux_r)
1751 2 : 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 2 : mpi_io=mpi_io)
1756 : CALL cp_print_key_finished_output(unit_nr, logger, resp_section, &
1757 2 : "PRINT%V_RESP_CUBE", mpi_io=mpi_io)
1758 2 : CALL auxbas_pw_pool%give_back_pw(aux_r)
1759 : END IF
1760 :
1761 : !*** RMS and RRMS
1762 14 : sum_diff = 0.0_dp
1763 14 : sum_hartree = 0.0_dp
1764 14 : rms = 0.0_dp
1765 14 : rrms = 0.0_dp
1766 2130 : DO ip = 1, resp_env%npoints_proc
1767 2116 : jx = resp_env%fitpoints(1, ip)
1768 2116 : jy = resp_env%fitpoints(2, ip)
1769 2116 : jz = resp_env%fitpoints(3, ip)
1770 : sum_diff = sum_diff + (v_hartree_rspace%array(jx, jy, jz) - &
1771 2116 : v_resp_rspace%array(jx, jy, jz))**2
1772 2130 : sum_hartree = sum_hartree + v_hartree_rspace%array(jx, jy, jz)**2
1773 : END DO
1774 14 : CALL para_env%sum(sum_diff)
1775 14 : CALL para_env%sum(sum_hartree)
1776 14 : rms = SQRT(sum_diff/resp_env%npoints)
1777 14 : rrms = SQRT(sum_diff/sum_hartree)
1778 14 : IF (output_runinfo > 0) THEN
1779 : WRITE (output_runinfo, '(2X,A,T69,ES12.5)') "Root-mean-square (RMS) "// &
1780 7 : "error of RESP fit:", rms
1781 : WRITE (output_runinfo, '(2X,A,T69,ES12.5,/)') "Relative root-mean-square "// &
1782 7 : "(RRMS) error of RESP fit:", rrms
1783 : END IF
1784 :
1785 14 : CALL v_resp_rspace%release()
1786 :
1787 14 : CALL timestop(handle)
1788 :
1789 14 : END SUBROUTINE print_pot_from_resp_charges
1790 :
1791 0 : END MODULE qs_resp
|
