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