(git:58e3e09)
qmmm_tb_coulomb.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 QMMM Coulomb contributions in TB
10 !> \author JGH
11 ! **************************************************************************************************
13  USE atomic_kind_types, ONLY: atomic_kind_type,&
15  USE atprop_types, ONLY: atprop_type
16  USE cell_types, ONLY: cell_type,&
17  get_cell
18  USE cp_control_types, ONLY: dft_control_type
19  USE dbcsr_api, ONLY: dbcsr_add,&
20  dbcsr_get_block_p,&
21  dbcsr_iterator_blocks_left,&
22  dbcsr_iterator_next_block,&
23  dbcsr_iterator_start,&
24  dbcsr_iterator_stop,&
25  dbcsr_iterator_type,&
26  dbcsr_p_type
27  USE distribution_1d_types, ONLY: distribution_1d_type
29  ewald_environment_type
31  USE ewald_pw_types, ONLY: ewald_pw_type
32  USE kinds, ONLY: dp
33  USE mathconstants, ONLY: oorootpi,&
34  pi
35  USE message_passing, ONLY: mp_para_env_type
36  USE particle_types, ONLY: particle_type
37  USE pw_poisson_types, ONLY: do_ewald_ewald,&
39  do_ewald_pme,&
41  USE qmmm_types_low, ONLY: qmmm_env_qm_type
42  USE qs_energy_types, ONLY: qs_energy_type
43  USE qs_environment_types, ONLY: get_qs_env,&
44  qs_environment_type
45  USE qs_force_types, ONLY: qs_force_type
46  USE qs_rho_types, ONLY: qs_rho_get,&
47  qs_rho_type
48  USE virial_types, ONLY: virial_type
49 #include "./base/base_uses.f90"
50 
51  IMPLICIT NONE
52 
53  PRIVATE
54 
55  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_tb_coulomb'
56 
57  PUBLIC :: build_tb_coulomb_qmqm
58 
59 CONTAINS
60 
61 ! **************************************************************************************************
62 !> \brief ...
63 !> \param qs_env ...
64 !> \param ks_matrix ...
65 !> \param rho ...
66 !> \param mcharge ...
67 !> \param energy ...
68 !> \param calculate_forces ...
69 !> \param just_energy ...
70 ! **************************************************************************************************
71  SUBROUTINE build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, &
72  calculate_forces, just_energy)
73 
74  TYPE(qs_environment_type), INTENT(IN) :: qs_env
75  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
76  TYPE(qs_rho_type), POINTER :: rho
77  REAL(dp), DIMENSION(:) :: mcharge
78  TYPE(qs_energy_type), POINTER :: energy
79  LOGICAL, INTENT(in) :: calculate_forces, just_energy
80 
81  CHARACTER(len=*), PARAMETER :: routinen = 'build_tb_coulomb_qmqm'
82 
83  INTEGER :: atom_i, atom_j, blk, ewald_type, handle, &
84  i, ia, iatom, ikind, jatom, jkind, &
85  natom, nmat
86  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
87  INTEGER, DIMENSION(3) :: periodic
88  LOGICAL :: found, use_virial
89  REAL(kind=dp) :: alpha, deth, dfr, dr, fi, fr, gmij
90  REAL(kind=dp), DIMENSION(3) :: fij, rij
91  REAL(kind=dp), DIMENSION(:, :), POINTER :: dsblock, gmcharge, ksblock, ksblock_2, &
92  pblock, sblock
93  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
94  TYPE(atprop_type), POINTER :: atprop
95  TYPE(cell_type), POINTER :: cell, mm_cell
96  TYPE(dbcsr_iterator_type) :: iter
97  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p, matrix_s
98  TYPE(dft_control_type), POINTER :: dft_control
99  TYPE(distribution_1d_type), POINTER :: local_particles
100  TYPE(ewald_environment_type), POINTER :: ewald_env
101  TYPE(ewald_pw_type), POINTER :: ewald_pw
102  TYPE(mp_para_env_type), POINTER :: para_env
103  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
104  TYPE(qmmm_env_qm_type), POINTER :: qmmm_env_qm
105  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
106  TYPE(virial_type), POINTER :: virial
107 
108  CALL timeset(routinen, handle)
109 
110  NULLIFY (matrix_p, matrix_s, virial, atprop, dft_control)
111 
112  use_virial = .false.
113 
114  IF (calculate_forces) THEN
115  nmat = 4
116  ELSE
117  nmat = 1
118  END IF
119 
120  natom = SIZE(mcharge)
121  ALLOCATE (gmcharge(natom, nmat))
122  gmcharge = 0._dp
123 
124  CALL get_qs_env(qs_env, &
125  particle_set=particle_set, &
126  cell=cell, &
127  virial=virial, &
128  atprop=atprop, &
129  dft_control=dft_control)
130 
131  IF (calculate_forces) THEN
132  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
133  END IF
134 
135  ! Qm-QM long range correction for QMMM calculations
136  ! no atomic energy evaluation
137  cpassert(.NOT. atprop%energy)
138  ! no stress tensor possible for QMMM
139  cpassert(.NOT. use_virial)
140  qmmm_env_qm => qs_env%qmmm_env_qm
141  ewald_env => qmmm_env_qm%ewald_env
142  ewald_pw => qmmm_env_qm%ewald_pw
143  CALL get_qs_env(qs_env=qs_env, super_cell=mm_cell)
144  CALL get_cell(cell=mm_cell, periodic=periodic, deth=deth)
145  CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
146  gmcharge = 0.0_dp
147  ! direct sum for overlap and local correction
148  CALL get_qs_env(qs_env=qs_env, &
149  atomic_kind_set=atomic_kind_set, &
150  local_particles=local_particles, &
151  force=force, para_env=para_env)
152  DO ikind = 1, SIZE(local_particles%n_el)
153  DO ia = 1, local_particles%n_el(ikind)
154  iatom = local_particles%list(ikind)%array(ia)
155  DO jatom = 1, iatom - 1
156  rij = particle_set(iatom)%r - particle_set(jatom)%r
157  ! no pbc(rij,mm_cell) at this point, we assume that QM particles are
158  ! inside QM box and QM box << MM box
159  dr = sqrt(sum(rij(:)**2))
160  ! local (unit cell) correction 1/R
161  gmcharge(iatom, 1) = gmcharge(iatom, 1) - mcharge(jatom)/dr
162  gmcharge(jatom, 1) = gmcharge(jatom, 1) - mcharge(iatom)/dr
163  DO i = 2, nmat
164  gmcharge(iatom, i) = gmcharge(iatom, i) - rij(i - 1)*mcharge(jatom)/dr**3
165  gmcharge(jatom, i) = gmcharge(jatom, i) + rij(i - 1)*mcharge(iatom)/dr**3
166  END DO
167  ! overlap correction
168  fr = erfc(alpha*dr)/dr
169  gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)*fr
170  gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)*fr
171  IF (nmat > 1) THEN
172  dfr = -2._dp*alpha*exp(-alpha*alpha*dr*dr)*oorootpi/dr - fr/dr
173  dfr = -dfr/dr
174  DO i = 2, nmat
175  gmcharge(iatom, i) = gmcharge(iatom, i) - rij(i - 1)*mcharge(jatom)*dfr
176  gmcharge(jatom, i) = gmcharge(jatom, i) + rij(i - 1)*mcharge(iatom)*dfr
177  END DO
178  END IF
179  END DO
180  END DO
181  END DO
182 
183  SELECT CASE (ewald_type)
184  CASE DEFAULT
185  cpabort("Invalid Ewald type")
186  CASE (do_ewald_none)
187  cpabort("Not allowed with DFTB")
188  CASE (do_ewald_ewald)
189  cpabort("Standard Ewald not implemented in DFTB")
190  CASE (do_ewald_pme)
191  cpabort("PME not implemented in DFTB")
192  CASE (do_ewald_spme)
193  CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, mm_cell, &
194  gmcharge, mcharge, calculate_forces, virial, use_virial, atprop)
195  END SELECT
196  !
197  CALL para_env%sum(gmcharge(:, 1))
198  !
199  ! add self charge interaction and background charge contribution
200  gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:)
201  IF (any(periodic(:) == 1)) THEN
202  gmcharge(:, 1) = gmcharge(:, 1) - pi/alpha**2/deth
203  END IF
204  !
205  energy%qmmm_el = energy%qmmm_el + 0.5_dp*sum(mcharge(:)*gmcharge(:, 1))
206  !
207  IF (calculate_forces) THEN
208  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
209  kind_of=kind_of, &
210  atom_of_kind=atom_of_kind)
211  END IF
212  !
213  IF (.NOT. just_energy) THEN
214  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
215  CALL qs_rho_get(rho, rho_ao=matrix_p)
216 
217  IF (calculate_forces .AND. SIZE(matrix_p) == 2) THEN
218  CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
219  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
220  END IF
221 
222  CALL dbcsr_iterator_start(iter, ks_matrix(1, 1)%matrix)
223  DO WHILE (dbcsr_iterator_blocks_left(iter))
224  CALL dbcsr_iterator_next_block(iter, iatom, jatom, ksblock, blk)
225  NULLIFY (sblock, ksblock_2)
226  IF (SIZE(ks_matrix, 1) > 1) THEN
227  CALL dbcsr_get_block_p(matrix=ks_matrix(2, 1)%matrix, &
228  row=iatom, col=jatom, block=ksblock_2, found=found)
229  END IF
230  CALL dbcsr_get_block_p(matrix=matrix_s(1)%matrix, &
231  row=iatom, col=jatom, block=sblock, found=found)
232  gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
233  ksblock = ksblock - gmij*sblock
234  IF (SIZE(ks_matrix, 1) > 1) ksblock_2 = ksblock_2 - gmij*sblock
235  IF (calculate_forces) THEN
236  ikind = kind_of(iatom)
237  atom_i = atom_of_kind(iatom)
238  jkind = kind_of(jatom)
239  atom_j = atom_of_kind(jatom)
240  NULLIFY (pblock)
241  CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
242  row=iatom, col=jatom, block=pblock, found=found)
243  DO i = 1, 3
244  NULLIFY (dsblock)
245  CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
246  row=iatom, col=jatom, block=dsblock, found=found)
247  fi = -2.0_dp*gmij*sum(pblock*dsblock)
248  force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
249  force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
250  fij(i) = fi
251  END DO
252  END IF
253  END DO
254  CALL dbcsr_iterator_stop(iter)
255  IF (calculate_forces .AND. SIZE(matrix_p) == 2) THEN
256  CALL dbcsr_add(matrix_p(1)%matrix, matrix_p(2)%matrix, &
257  alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
258  END IF
259  END IF
260 
261  DEALLOCATE (gmcharge)
262 
263  CALL timestop(handle)
264 
265  END SUBROUTINE build_tb_coulomb_qmqm
266 
267 END MODULE qmmm_tb_coulomb
268 
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Holds information on atomic properties.
Definition: atprop_types.F:14
Handles all functions related to the CELL.
Definition: cell_types.F:15
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition: cell_types.F:195
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Calculation of Ewald contributions in DFTB.
subroutine, public tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, gmcharge, mcharge, calculate_forces, virial, use_virial, atprop)
...
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public oorootpi
real(kind=dp), parameter, public pi
Interface to the message passing library MPI.
Define the data structure for the particle information.
functions related to the poisson solver on regular grids
integer, parameter, public do_ewald_pme
integer, parameter, public do_ewald_ewald
integer, parameter, public do_ewald_none
integer, parameter, public do_ewald_spme
Calculation of QMMM Coulomb contributions in TB.
subroutine, public build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
...
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.
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