(git:b279b6b)
qmmm_se_forces.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 !> \brief Calculation of the derivative of the QMMM Hamiltonian integral
10 !> matrix <a|\sum_i q_i|b> for semi-empirical methods
11 !> \author Teodoro Laino - 04.2007 [tlaino]
12 ! **************************************************************************************************
14  USE atomic_kind_types, ONLY: atomic_kind_type,&
16  USE cell_types, ONLY: cell_type,&
17  pbc
18  USE cp_control_types, ONLY: dft_control_type
19  USE dbcsr_api, ONLY: dbcsr_get_block_p,&
20  dbcsr_p_type
21  USE input_constants, ONLY: &
24  USE kinds, ONLY: dp
25  USE message_passing, ONLY: mp_para_env_type
27  USE particle_types, ONLY: particle_type
28  USE qmmm_types_low, ONLY: qmmm_env_qm_type,&
29  qmmm_pot_p_type,&
30  qmmm_pot_type
32  USE qs_environment_types, ONLY: get_qs_env,&
33  qs_environment_type
34  USE qs_kind_types, ONLY: get_qs_kind,&
35  qs_kind_type
36  USE qs_ks_qmmm_types, ONLY: qs_ks_qmmm_env_type
37  USE qs_rho_types, ONLY: qs_rho_get,&
38  qs_rho_type
41  drotnuc
43  se_int_control_type,&
44  se_taper_type,&
47  semi_empirical_type,&
51 #include "./base/base_uses.f90"
52 
53  IMPLICIT NONE
54 
55  PRIVATE
56 
57  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_se_forces'
58  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
59  PUBLIC :: deriv_se_qmmm_matrix
60 
61 CONTAINS
62 
63 ! **************************************************************************************************
64 !> \brief Constructs the derivative w.r.t. 1-el semi-empirical hamiltonian
65 !> QMMM terms
66 !> \param qs_env ...
67 !> \param qmmm_env ...
68 !> \param particles_mm ...
69 !> \param mm_cell ...
70 !> \param para_env ...
71 !> \param calc_force ...
72 !> \param Forces ...
73 !> \param Forces_added_charges ...
74 !> \author Teodoro Laino 04.2007 [created]
75 ! **************************************************************************************************
76  SUBROUTINE deriv_se_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
77  calc_force, Forces, Forces_added_charges)
78 
79  TYPE(qs_environment_type), POINTER :: qs_env
80  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env
81  TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
82  TYPE(cell_type), POINTER :: mm_cell
83  TYPE(mp_para_env_type), POINTER :: para_env
84  LOGICAL, INTENT(in), OPTIONAL :: calc_force
85  REAL(kind=dp), DIMENSION(:, :), POINTER :: forces, forces_added_charges
86 
87  CHARACTER(len=*), PARAMETER :: routinen = 'deriv_se_qmmm_matrix'
88 
89  INTEGER :: handle, i, iatom, ikind, iqm, ispin, &
90  itype, natom, natorb_a, nkind, &
91  number_qm_atoms
92  INTEGER, DIMENSION(:), POINTER :: list
93  LOGICAL :: anag, defined, found
94  REAL(kind=dp) :: delta, enuclear
95  REAL(kind=dp), DIMENSION(:, :), POINTER :: forces_qm, p_block_a
96  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
97  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
98  TYPE(dft_control_type), POINTER :: dft_control
99  TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
100  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
101  TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
102  TYPE(qs_rho_type), POINTER :: rho
103  TYPE(se_int_control_type) :: se_int_control
104  TYPE(se_taper_type), POINTER :: se_taper
105  TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm
106 
107  CALL timeset(routinen, handle)
108  IF (calc_force) THEN
109  NULLIFY (rho, atomic_kind_set, qs_kind_set, se_taper)
110  NULLIFY (se_kind_a, se_kind_mm, particles_qm)
111  CALL get_qs_env(qs_env=qs_env, &
112  rho=rho, &
113  se_taper=se_taper, &
114  atomic_kind_set=atomic_kind_set, &
115  qs_kind_set=qs_kind_set, &
116  ks_qmmm_env=ks_qmmm_env_loc, &
117  dft_control=dft_control, &
118  particle_set=particles_qm, &
119  natom=number_qm_atoms)
120  SELECT CASE (dft_control%qs_control%method_id)
123  ! Go on with the calculation..
124  CASE DEFAULT
125  ! Otherwise stop..
126  cpabort("Method not available")
127  END SELECT
128  anag = dft_control%qs_control%se_control%analytical_gradients
129  delta = dft_control%qs_control%se_control%delta
130  ! Setup SE integral control type
132  se_int_control, shortrange=.false., do_ewald_r3=.false., &
133  do_ewald_gks=.false., integral_screening=dft_control%qs_control%se_control%integral_screening, &
134  max_multipole=do_multipole_none, pc_coulomb_int=.false.)
135 
136  ! Create a fake semi-empirical type to handle the classical atom
137  ALLOCATE (forces_qm(3, number_qm_atoms))
138  CALL semi_empirical_create(se_kind_mm)
139  CALL se_param_set_default(se_kind_mm, 0, do_method_pchg)
140  itype = get_se_type(se_kind_mm%typ)
141  nkind = SIZE(atomic_kind_set)
142  enuclear = 0.0_dp
143  forces_qm = 0.0_dp
144  CALL qs_rho_get(rho, rho_ao=matrix_p)
145 
146  DO ispin = 1, dft_control%nspins
147  iqm = 0
148  kinds: DO ikind = 1, nkind
149  CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=list)
150  CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
151  CALL get_se_param(se_kind_a, &
152  defined=defined, &
153  natorb=natorb_a)
154  IF (.NOT. defined .OR. natorb_a < 1) cycle
155  atoms: DO i = 1, SIZE(list)
156  iqm = iqm + 1
157  iatom = list(i)
158  ! Give back block
159  NULLIFY (p_block_a)
160  CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
161  row=iatom, col=iatom, block=p_block_a, found=found)
162 
163  IF (ASSOCIATED(p_block_a)) THEN
164  ! Expand derivative of geometrical factors
165  CALL deriv_se_qmmm_matrix_low(p_block_a, &
166  se_kind_a, &
167  se_kind_mm, &
168  qmmm_env%Potentials, &
169  particles_mm, &
170  qmmm_env%mm_atom_chrg, &
171  qmmm_env%mm_atom_index, &
172  mm_cell, &
173  iatom, &
174  itype, &
175  forces, &
176  forces_qm(:, iqm), &
177  se_taper, &
178  se_int_control, &
179  anag, &
180  delta, &
181  qmmm_env%spherical_cutoff, &
182  particles_qm)
183  ! Possibly added charges
184  IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
185  CALL deriv_se_qmmm_matrix_low(p_block_a, &
186  se_kind_a, &
187  se_kind_mm, &
188  qmmm_env%added_charges%potentials, &
189  qmmm_env%added_charges%added_particles, &
190  qmmm_env%added_charges%mm_atom_chrg, &
191  qmmm_env%added_charges%mm_atom_index, &
192  mm_cell, &
193  iatom, &
194  itype, &
195  forces_added_charges, &
196  forces_qm(:, iqm), &
197  se_taper, &
198  se_int_control, &
199  anag, &
200  delta, &
201  qmmm_env%spherical_cutoff, &
202  particles_qm)
203  END IF
204  END IF
205  END DO atoms
206  END DO kinds
207  END DO
208  cpassert(iqm == number_qm_atoms)
209  ! Transfer QM gradients to the QM particles..
210  CALL para_env%sum(forces_qm)
211  iqm = 0
212  DO ikind = 1, nkind
213  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
214  CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
215  CALL get_se_param(se_kind_a, &
216  defined=defined, &
217  natorb=natorb_a)
218  IF (.NOT. defined .OR. natorb_a < 1) cycle
219  DO i = 1, SIZE(list)
220  iqm = iqm + 1
221  iatom = qmmm_env%qm_atom_index(list(i))
222  particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + forces_qm(:, iqm)
223  END DO
224  END DO
225  ! MM forces will be handled directly from the QMMM module in the same way
226  ! as for GPW/GAPW methods
227  DEALLOCATE (forces_qm)
228  CALL semi_empirical_release(se_kind_mm)
229 
230  END IF
231  CALL timestop(handle)
232  END SUBROUTINE deriv_se_qmmm_matrix
233 
234 ! **************************************************************************************************
235 !> \brief Low Level : Computes derivatives of the 1-el semi-empirical QMMM
236 !> hamiltonian block w.r.t. MM and QM coordinates
237 !> \param p_block_a ...
238 !> \param se_kind_a ...
239 !> \param se_kind_mm ...
240 !> \param potentials ...
241 !> \param particles_mm ...
242 !> \param mm_charges ...
243 !> \param mm_atom_index ...
244 !> \param mm_cell ...
245 !> \param IndQM ...
246 !> \param itype ...
247 !> \param forces ...
248 !> \param forces_qm ...
249 !> \param se_taper ...
250 !> \param se_int_control ...
251 !> \param anag ...
252 !> \param delta ...
253 !> \param qmmm_spherical_cutoff ...
254 !> \param particles_qm ...
255 !> \author Teodoro Laino 04.2007 [created]
256 ! **************************************************************************************************
257  SUBROUTINE deriv_se_qmmm_matrix_low(p_block_a, se_kind_a, se_kind_mm, &
258  potentials, particles_mm, mm_charges, mm_atom_index, &
259  mm_cell, IndQM, itype, forces, forces_qm, se_taper, &
260  se_int_control, anag, delta, qmmm_spherical_cutoff, particles_qm)
261 
262  REAL(kind=dp), DIMENSION(:, :), POINTER :: p_block_a
263  TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm
264  TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER :: potentials
265  TYPE(particle_type), DIMENSION(:), POINTER :: particles_mm
266  REAL(kind=dp), DIMENSION(:), POINTER :: mm_charges
267  INTEGER, DIMENSION(:), POINTER :: mm_atom_index
268  TYPE(cell_type), POINTER :: mm_cell
269  INTEGER, INTENT(IN) :: indqm, itype
270  REAL(kind=dp), DIMENSION(:, :), POINTER :: forces
271  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: forces_qm
272  TYPE(se_taper_type), POINTER :: se_taper
273  TYPE(se_int_control_type), INTENT(IN) :: se_int_control
274  LOGICAL, INTENT(IN) :: anag
275  REAL(kind=dp), INTENT(IN) :: delta, qmmm_spherical_cutoff(2)
276  TYPE(particle_type), DIMENSION(:), POINTER :: particles_qm
277 
278  CHARACTER(len=*), PARAMETER :: routinen = 'deriv_se_qmmm_matrix_low'
279 
280  INTEGER :: handle, i1, i1l, i2, imm, imp, indmm, &
281  ipot, j1, j1l
282  REAL(kind=dp) :: rt1, rt2, rt3, sph_chrg_factor
283  REAL(kind=dp), DIMENSION(3) :: denuc, force_ab, r_pbc, rij
284  REAL(kind=dp), DIMENSION(3, 45) :: de1b
285  TYPE(qmmm_pot_type), POINTER :: pot
286 
287  CALL timeset(routinen, handle)
288  ! Loop Over MM atoms - parallelization over MM atoms...
289  ! Loop over Pot stores atoms with the same charge
290  mainlooppot: DO ipot = 1, SIZE(potentials)
291  pot => potentials(ipot)%Pot
292  ! Loop over atoms belonging to this type
293  loopmm: DO imp = 1, SIZE(pot%mm_atom_index)
294  imm = pot%mm_atom_index(imp)
295  indmm = mm_atom_index(imm)
296  r_pbc = pbc(particles_mm(indmm)%r - particles_qm(indqm)%r, mm_cell)
297  rt1 = r_pbc(1)
298  rt2 = r_pbc(2)
299  rt3 = r_pbc(3)
300  rij = (/rt1, rt2, rt3/)
301  se_kind_mm%zeff = mm_charges(imm)
302  ! Computes the screening factor for the spherical cutoff
303  IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
304  CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
305  se_kind_mm%zeff = se_kind_mm%zeff*sph_chrg_factor
306  END IF
307  IF (abs(se_kind_mm%zeff) <= epsilon(0.0_dp)) cycle
308  ! Integrals derivatives involving QM - MM atoms
309  CALL drotnuc(se_kind_a, se_kind_mm, rij, itype=itype, de1b=de1b, &
310  se_int_control=se_int_control, anag=anag, delta=delta, &
311  se_taper=se_taper)
312  CALL dcorecore(se_kind_a, se_kind_mm, rij, itype=itype, denuc=denuc, &
313  se_int_control=se_int_control, anag=anag, delta=delta, &
314  se_taper=se_taper)
315  ! Nucler - Nuclear term
316  force_ab(1:3) = -denuc(1:3)
317  ! Force contribution from the QMMM Hamiltonian
318  i2 = 0
319  DO i1l = 1, se_kind_a%natorb
320  i1 = se_orbital_pointer(i1l)
321  DO j1l = 1, i1l - 1
322  j1 = se_orbital_pointer(j1l)
323  i2 = i2 + 1
324  force_ab = force_ab - 2.0_dp*de1b(:, i2)*p_block_a(i1, j1)
325  END DO
326  j1 = se_orbital_pointer(j1l)
327  i2 = i2 + 1
328  force_ab = force_ab - de1b(:, i2)*p_block_a(i1, j1)
329  END DO
330  ! The array of QM forces are really the forces
331  forces_qm(:) = forces_qm(:) - force_ab
332  ! The one of MM atoms are instead gradients
333  forces(:, imm) = forces(:, imm) - force_ab
334  END DO loopmm
335  END DO mainlooppot
336  CALL timestop(handle)
337  END SUBROUTINE deriv_se_qmmm_matrix_low
338 
339 END MODULE qmmm_se_forces
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Handles all functions related to the CELL.
Definition: cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_method_pchg
integer, parameter, public do_method_pdg
integer, parameter, public do_method_pnnl
integer, parameter, public do_method_rm1
integer, parameter, public do_method_pm3
integer, parameter, public do_method_mndo
integer, parameter, public do_method_mndod
integer, parameter, public do_method_am1
integer, parameter, public do_method_pm6fm
integer, parameter, public do_method_pm6
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
Interface to the message passing library MPI.
Multipole structure: for multipole (fixed and induced) in FF based MD.
integer, parameter, public do_multipole_none
Define the data structure for the particle information.
Calculation of the derivative of the QMMM Hamiltonian integral matrix <a|\sum_i q_i|b> for semi-empir...
subroutine, public deriv_se_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, calc_force, Forces, Forces_added_charges)
Constructs the derivative w.r.t. 1-el semi-empirical hamiltonian QMMM terms.
subroutine, public spherical_cutoff_factor(spherical_cutoff, rij, factor)
Computes a spherical cutoff factor for the QMMM interactions.
Definition: qmmm_util.F:615
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.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
Arrays of parameters used in the semi-empirical calculations \References Everywhere in this module TC...
integer, dimension(9), public se_orbital_pointer
Set of wrappers for semi-empirical analytical/numerical Integrals routines.
subroutine, public dcorecore(sepi, sepj, rij, denuc, itype, delta, anag, se_int_control, se_taper)
wrapper for numerical/analytical routines
subroutine, public drotnuc(sepi, sepj, rij, de1b, de2a, itype, delta, anag, se_int_control, se_taper)
wrapper for numerical/analytical routines
Definition of the semi empirical parameter types.
subroutine, public semi_empirical_create(sep)
Allocate semi-empirical type.
subroutine, public setup_se_int_control_type(se_int_control, shortrange, do_ewald_r3, do_ewald_gks, integral_screening, max_multipole, pc_coulomb_int)
Setup the Semiempirical integral control type.
subroutine, public get_se_param(sep, name, typ, defined, z, zeff, natorb, eheat, beta, sto_exponents, uss, upp, udd, uff, alp, eisol, gss, gsp, gpp, gp2, acoul, nr, de, ass, asp, app, hsp, gsd, gpd, gdd, ppddg, dpddg, ngauss)
Get info from the semi-empirical type.
subroutine, public semi_empirical_release(sep)
Deallocate the semi-empirical type.
Working with the semi empirical parameter types.
integer function, public get_se_type(se_method)
Gives back the unique semi_empirical METHOD type.
subroutine, public se_param_set_default(sep, z, method)
Initialize parameter for a semi_empirival type.