LCOV - code coverage report
Current view: top level - src - qs_resp.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 763 824 92.6 %
Date: 2024-03-28 07:31:50 Functions: 16 18 88.9 %

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

Generated by: LCOV version 1.15