(git:b279b6b)
qs_ks_qmmm_methods.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \author T.Laino
10 ! **************************************************************************************************
12  USE cp_control_types, ONLY: dft_control_type
14  cp_logger_type
17  USE cube_utils, ONLY: cube_info_type,&
19  USE d3_poly, ONLY: init_d3_poly_module
20  USE input_constants, ONLY: do_qmmm_gauss,&
22  USE input_section_types, ONLY: section_vals_type
23  USE kinds, ONLY: dp
24  USE pw_env_types, ONLY: pw_env_get,&
26  USE pw_methods, ONLY: pw_axpy,&
27  pw_integral_ab
28  USE pw_pool_types, ONLY: pw_pool_p_type,&
29  pw_pool_type
30  USE pw_types, ONLY: pw_r3d_rs_type
31  USE qmmm_types_low, ONLY: qmmm_env_qm_type
32  USE qs_environment_types, ONLY: get_qs_env,&
33  qs_environment_type,&
35  USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
36 #include "./base/base_uses.f90"
37 
38  IMPLICIT NONE
39 
40  PRIVATE
41 
42  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_qmmm_methods'
43 
46 
47 CONTAINS
48 
49 ! **************************************************************************************************
50 !> \brief Initialize the ks_qmmm_env
51 !> \param qs_env ...
52 !> \param qmmm_env ...
53 !> \par History
54 !> 05.2004 created [tlaino]
55 !> \author Teodoro Laino
56 ! **************************************************************************************************
57  SUBROUTINE ks_qmmm_env_rebuild(qs_env, qmmm_env)
58  TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
59  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
60 
61  TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env
62 
63  NULLIFY (ks_qmmm_env)
64  CALL get_qs_env(qs_env=qs_env, &
65  ks_qmmm_env=ks_qmmm_env)
66 
67  ! *** allocate the ks_qmmm env if not allocated yet!**
68  IF (.NOT. ASSOCIATED(ks_qmmm_env)) THEN
69  ALLOCATE (ks_qmmm_env)
70  CALL qs_ks_qmmm_create(ks_qmmm_env=ks_qmmm_env, qs_env=qs_env, &
71  qmmm_env=qmmm_env)
72  CALL set_qs_env(qs_env=qs_env, ks_qmmm_env=ks_qmmm_env)
73  END IF
74  END SUBROUTINE ks_qmmm_env_rebuild
75 
76 ! **************************************************************************************************
77 !> \brief allocates and initializes the given ks_qmmm_env.
78 !> \param ks_qmmm_env the ks_qmmm env to be initialized
79 !> \param qs_env the qs environment
80 !> \param qmmm_env ...
81 !> \par History
82 !> 05.2004 created [tlaino]
83 !> \author Teodoro Laino
84 ! **************************************************************************************************
85  SUBROUTINE qs_ks_qmmm_create(ks_qmmm_env, qs_env, qmmm_env)
86  TYPE(qs_ks_qmmm_env_type), INTENT(OUT) :: ks_qmmm_env
87  TYPE(qs_environment_type), POINTER :: qs_env
88  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
89 
90  CHARACTER(len=*), PARAMETER :: routinen = 'qs_ks_qmmm_create'
91 
92  INTEGER :: handle, igrid
93  TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
94  TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pools
95  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
96 
97  CALL timeset(routinen, handle)
98 
99  NULLIFY (ks_qmmm_env%pw_env, &
100  ks_qmmm_env%cube_info)
101  NULLIFY (auxbas_pw_pool)
102  CALL get_qs_env(qs_env=qs_env, &
103  pw_env=ks_qmmm_env%pw_env)
104  CALL pw_env_get(ks_qmmm_env%pw_env, auxbas_pw_pool=auxbas_pw_pool)
105  CALL pw_env_retain(ks_qmmm_env%pw_env)
106 
107  ks_qmmm_env%pc_ener = 0.0_dp
108  ks_qmmm_env%n_evals = 0
109 
110  CALL auxbas_pw_pool%create_pw(ks_qmmm_env%v_qmmm_rspace)
111 
112  IF ((qmmm_env%qmmm_coupl_type == do_qmmm_gauss) .OR. (qmmm_env%qmmm_coupl_type == do_qmmm_swave)) THEN
113  CALL init_d3_poly_module() ! a fairly arbitrary but sufficient spot to do this
114  CALL pw_env_get(ks_qmmm_env%pw_env, pw_pools=pools)
115  ALLOCATE (cube_info(SIZE(pools)))
116  DO igrid = 1, SIZE(pools)
117  CALL init_cube_info(cube_info(igrid), &
118  pools(igrid)%pool%pw_grid%dr(:), &
119  pools(igrid)%pool%pw_grid%dh(:, :), &
120  pools(igrid)%pool%pw_grid%dh_inv(:, :), &
121  pools(igrid)%pool%pw_grid%orthorhombic, &
122  qmmm_env%maxRadius(igrid))
123  END DO
124  ks_qmmm_env%cube_info => cube_info
125  END IF
126  NULLIFY (ks_qmmm_env%matrix_h)
127  !
128 
129  CALL timestop(handle)
130 
131  END SUBROUTINE qs_ks_qmmm_create
132 
133 ! **************************************************************************************************
134 !> \brief Computes the contribution to the total energy of the QM/MM
135 !> electrostatic coupling
136 !> \param qs_env ...
137 !> \param rho ...
138 !> \param v_qmmm ...
139 !> \param qmmm_energy ...
140 !> \par History
141 !> 05.2004 created [tlaino]
142 !> \author Teodoro Laino
143 ! **************************************************************************************************
144  SUBROUTINE qmmm_calculate_energy(qs_env, rho, v_qmmm, qmmm_energy)
145  TYPE(qs_environment_type), POINTER :: qs_env
146  TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: rho
147  TYPE(pw_r3d_rs_type), INTENT(IN) :: v_qmmm
148  REAL(kind=dp), INTENT(INOUT) :: qmmm_energy
149 
150  CHARACTER(len=*), PARAMETER :: routinen = 'qmmm_calculate_energy'
151 
152  INTEGER :: handle, ispin, output_unit
153  TYPE(cp_logger_type), POINTER :: logger
154  TYPE(dft_control_type), POINTER :: dft_control
155  TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs
156  TYPE(section_vals_type), POINTER :: input
157 
158  CALL timeset(routinen, handle)
159  NULLIFY (dft_control, input, logger)
160  logger => cp_get_default_logger()
161 
162  CALL get_qs_env(qs_env=qs_env, &
163  dft_control=dft_control, &
164  input=input)
165 
166  output_unit = cp_print_key_unit_nr(logger, input, "QMMM%PRINT%PROGRAM_RUN_INFO", &
167  extension=".qmmmLog")
168  IF (output_unit > 0) WRITE (unit=output_unit, fmt="(T3,A)") &
169  "Adding QM/MM electrostatic potential to the Kohn-Sham potential."
170 
171  qmmm_energy = 0.0_dp
172  DO ispin = 1, dft_control%nspins
173  qmmm_energy = qmmm_energy + pw_integral_ab(rho(ispin), v_qmmm)
174  END DO
175  IF (dft_control%qs_control%gapw) THEN
176  CALL get_qs_env(qs_env=qs_env, &
177  rho0_s_rs=rho0_s_rs)
178  cpassert(ASSOCIATED(rho0_s_rs))
179  qmmm_energy = qmmm_energy + pw_integral_ab(rho0_s_rs, v_qmmm)
180  END IF
181 
182  CALL cp_print_key_finished_output(output_unit, logger, input, &
183  "QMMM%PRINT%PROGRAM_RUN_INFO")
184 
185  CALL timestop(handle)
186  END SUBROUTINE qmmm_calculate_energy
187 
188 ! **************************************************************************************************
189 !> \brief Modify the hartree potential in order to include the QM/MM correction
190 !> \param v_hartree ...
191 !> \param v_qmmm ...
192 !> \param scale ...
193 !> \par History
194 !> 05.2004 created [tlaino]
195 !> \author Teodoro Laino
196 ! **************************************************************************************************
197  SUBROUTINE qmmm_modify_hartree_pot(v_hartree, v_qmmm, scale)
198  TYPE(pw_r3d_rs_type), INTENT(INOUT) :: v_hartree
199  TYPE(pw_r3d_rs_type), INTENT(IN) :: v_qmmm
200  REAL(kind=dp), INTENT(IN) :: scale
201 
202  CALL pw_axpy(v_qmmm, v_hartree, v_qmmm%pw_grid%dvol*scale)
203 
204  END SUBROUTINE qmmm_modify_hartree_pot
205 
206 END MODULE qs_ks_qmmm_methods
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
for a given dr()/dh(r) this will provide the bounds to be used if one wants to go over a sphere-subre...
Definition: cube_utils.F:18
subroutine, public init_cube_info(info, dr, dh, dh_inv, ortho, max_radius)
...
Definition: cube_utils.F:212
Routines to efficiently handle dense polynomial in 3 variables up to a given degree....
Definition: d3_poly.F:23
subroutine, public init_d3_poly_module()
initialization of the cache, is called by functions in this module that use cached values
Definition: d3_poly.F:74
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_qmmm_swave
integer, parameter, public do_qmmm_gauss
objects that represent the structure of input sections and the data contained in an input section
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
container for various plainwaves related things
Definition: pw_env_types.F:14
subroutine, public pw_env_retain(pw_env)
retains the pw_env (see doc/ReferenceCounting.html)
Definition: pw_env_types.F:161
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
subroutine, public qmmm_calculate_energy(qs_env, rho, v_qmmm, qmmm_energy)
Computes the contribution to the total energy of the QM/MM electrostatic coupling.
subroutine, public qmmm_modify_hartree_pot(v_hartree, v_qmmm, scale)
Modify the hartree potential in order to include the QM/MM correction.
subroutine, public ks_qmmm_env_rebuild(qs_env, qmmm_env)
Initialize the ks_qmmm_env.