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