(git:1f285aa)
ec_efield_local.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 Calculates the energy contribution and the mo_derivative of
10 !> a static electric field (nonperiodic)
11 !> \par History
12 !> Adjusted from qs_efield_local
13 !> \author JGH (10.2019)
14 ! **************************************************************************************************
16  USE ai_moments, ONLY: dipole_force
17  USE atomic_kind_types, ONLY: atomic_kind_type,&
20  USE basis_set_types, ONLY: gto_basis_set_p_type,&
21  gto_basis_set_type
22  USE cell_types, ONLY: cell_type,&
23  pbc
24  USE cp_control_types, ONLY: dft_control_type
25  USE dbcsr_api, ONLY: dbcsr_add,&
26  dbcsr_copy,&
27  dbcsr_get_block_p,&
28  dbcsr_p_type,&
29  dbcsr_set
30  USE ec_env_types, ONLY: energy_correction_type
31  USE kinds, ONLY: dp
32  USE message_passing, ONLY: mp_para_env_type
33  USE orbital_pointers, ONLY: ncoset
34  USE particle_types, ONLY: particle_type
35  USE qs_energy_types, ONLY: qs_energy_type
36  USE qs_environment_types, ONLY: get_qs_env,&
37  qs_environment_type
38  USE qs_force_types, ONLY: qs_force_type
39  USE qs_kind_types, ONLY: get_qs_kind,&
40  qs_kind_type
45  neighbor_list_iterator_p_type,&
47  neighbor_list_set_p_type
48  USE qs_period_efield_types, ONLY: efield_berry_type,&
51 #include "./base/base_uses.f90"
52 
53  IMPLICIT NONE
54 
55  PRIVATE
56 
57  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_efield_local'
58 
59  ! *** Public subroutines ***
60 
62 
63 ! **************************************************************************************************
64 
65 CONTAINS
66 
67 ! **************************************************************************************************
68 
69 ! **************************************************************************************************
70 !> \brief ...
71 !> \param qs_env ...
72 !> \param ec_env ...
73 !> \param calculate_forces ...
74 ! **************************************************************************************************
75  SUBROUTINE ec_efield_local_operator(qs_env, ec_env, calculate_forces)
76 
77  TYPE(qs_environment_type), POINTER :: qs_env
78  TYPE(energy_correction_type), POINTER :: ec_env
79  LOGICAL, INTENT(IN) :: calculate_forces
80 
81  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_efield_local_operator'
82 
83  INTEGER :: handle
84  REAL(dp), DIMENSION(3) :: rpoint
85  TYPE(dft_control_type), POINTER :: dft_control
86 
87  CALL timeset(routinen, handle)
88 
89  NULLIFY (dft_control)
90  CALL get_qs_env(qs_env, dft_control=dft_control)
91 
92  IF (dft_control%apply_efield) THEN
93  rpoint = 0.0_dp
94  CALL ec_efield_integrals(qs_env, ec_env, rpoint)
95  CALL ec_efield_mo_derivatives(qs_env, ec_env, rpoint, calculate_forces)
96  END IF
97 
98  CALL timestop(handle)
99 
100  END SUBROUTINE ec_efield_local_operator
101 
102 ! **************************************************************************************************
103 !> \brief ...
104 !> \param qs_env ...
105 !> \param ec_env ...
106 !> \param rpoint ...
107 ! **************************************************************************************************
108  SUBROUTINE ec_efield_integrals(qs_env, ec_env, rpoint)
109 
110  TYPE(qs_environment_type), POINTER :: qs_env
111  TYPE(energy_correction_type), POINTER :: ec_env
112  REAL(dp), DIMENSION(3), INTENT(IN) :: rpoint
113 
114  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_efield_integrals'
115 
116  INTEGER :: handle, i
117  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s
118  TYPE(efield_berry_type), POINTER :: efield, efieldref
119 
120  CALL timeset(routinen, handle)
121 
122  CALL get_qs_env(qs_env=qs_env, efield=efieldref)
123  efield => ec_env%efield
124  CALL init_efield_matrices(efield)
125  matrix_s => ec_env%matrix_s(:, 1)
126  ALLOCATE (dipmat(3))
127  DO i = 1, 3
128  ALLOCATE (dipmat(i)%matrix)
129  CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'DIP MAT')
130  CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
131  END DO
132  CALL build_local_moment_matrix(qs_env, dipmat, 1, rpoint, basis_type="HARRIS")
133  CALL set_efield_matrices(efield=efield, dipmat=dipmat)
134  ec_env%efield => efield
135 
136  CALL timestop(handle)
137 
138  END SUBROUTINE ec_efield_integrals
139 
140 ! **************************************************************************************************
141 !> \brief ...
142 !> \param qs_env ...
143 !> \param ec_env ...
144 !> \param rpoint ...
145 !> \param calculate_forces ...
146 ! **************************************************************************************************
147  SUBROUTINE ec_efield_mo_derivatives(qs_env, ec_env, rpoint, calculate_forces)
148  TYPE(qs_environment_type), POINTER :: qs_env
149  TYPE(energy_correction_type), POINTER :: ec_env
150  REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rpoint
151  LOGICAL :: calculate_forces
152 
153  CHARACTER(LEN=*), PARAMETER :: routinen = 'ec_efield_mo_derivatives'
154 
155  INTEGER :: atom_a, atom_b, handle, i, ia, iatom, icol, idir, ikind, irow, iset, ispin, &
156  jatom, jkind, jset, ldab, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
157  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
158  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
159  npgfb, nsgfa, nsgfb
160  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
161  LOGICAL :: found, trans
162  REAL(dp) :: charge, dab, fdir
163  REAL(dp), DIMENSION(3) :: ci, fieldpol, ra, rab, rac, rbc, ria
164  REAL(dp), DIMENSION(3, 3) :: forcea, forceb
165  REAL(dp), DIMENSION(:, :), POINTER :: p_block_a, p_block_b, pblock, pmat, work
166  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
167  REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
168  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
169  TYPE(cell_type), POINTER :: cell
170  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_ks
171  TYPE(dft_control_type), POINTER :: dft_control
172  TYPE(efield_berry_type), POINTER :: efield
173  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
174  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
175  TYPE(mp_para_env_type), POINTER :: para_env
176  TYPE(neighbor_list_iterator_p_type), &
177  DIMENSION(:), POINTER :: nl_iterator
178  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
179  POINTER :: sab_orb
180  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
181  TYPE(qs_energy_type), POINTER :: energy
182  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
183  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
184  TYPE(qs_kind_type), POINTER :: qs_kind
185 
186  CALL timeset(routinen, handle)
187 
188  CALL get_qs_env(qs_env, dft_control=dft_control, cell=cell, particle_set=particle_set)
189  CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
190  energy=energy, para_env=para_env, sab_orb=sab_orb)
191 
192  efield => ec_env%efield
193 
194  fieldpol = dft_control%efield_fields(1)%efield%polarisation* &
195  dft_control%efield_fields(1)%efield%strength
196 
197  ! nuclear contribution
198  natom = SIZE(particle_set)
199  IF (calculate_forces) THEN
200  CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force)
201  CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
202  END IF
203  ci = 0.0_dp
204  DO ia = 1, natom
205  CALL get_atomic_kind(particle_set(ia)%atomic_kind, kind_number=ikind)
206  CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
207  ria = particle_set(ia)%r - rpoint
208  ria = pbc(ria, cell)
209  ci(:) = ci(:) + charge*ria(:)
210  IF (calculate_forces) THEN
211  IF (para_env%mepos == 0) THEN
212  iatom = atom_of_kind(ia)
213  DO idir = 1, 3
214  force(ikind)%efield(idir, iatom) = force(ikind)%efield(idir, iatom) - fieldpol(idir)*charge
215  END DO
216  END IF
217  END IF
218  END DO
219 
220  IF (ec_env%should_update) THEN
221  ec_env%efield_nuclear = -sum(ci(:)*fieldpol(:))
222  ! Update KS matrix
223  matrix_ks => ec_env%matrix_h(:, 1)
224  dipmat => efield%dipmat
225  DO ispin = 1, SIZE(matrix_ks)
226  DO idir = 1, 3
227  CALL dbcsr_add(matrix_ks(ispin)%matrix, dipmat(idir)%matrix, &
228  alpha_scalar=1.0_dp, beta_scalar=fieldpol(idir))
229  END DO
230  END DO
231  END IF
232 
233  ! forces from the efield contribution
234  IF (calculate_forces) THEN
235  nkind = SIZE(qs_kind_set)
236  natom = SIZE(particle_set)
237 
238  ALLOCATE (basis_set_list(nkind))
239  DO ikind = 1, nkind
240  qs_kind => qs_kind_set(ikind)
241  CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type="HARRIS")
242  IF (ASSOCIATED(basis_set_a)) THEN
243  basis_set_list(ikind)%gto_basis_set => basis_set_a
244  ELSE
245  NULLIFY (basis_set_list(ikind)%gto_basis_set)
246  END IF
247  END DO
248  !
249  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
250  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
251  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
252  iatom=iatom, jatom=jatom, r=rab)
253  basis_set_a => basis_set_list(ikind)%gto_basis_set
254  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
255  basis_set_b => basis_set_list(jkind)%gto_basis_set
256  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
257  ! basis ikind
258  first_sgfa => basis_set_a%first_sgf
259  la_max => basis_set_a%lmax
260  la_min => basis_set_a%lmin
261  npgfa => basis_set_a%npgf
262  nseta = basis_set_a%nset
263  nsgfa => basis_set_a%nsgf_set
264  rpgfa => basis_set_a%pgf_radius
265  set_radius_a => basis_set_a%set_radius
266  sphi_a => basis_set_a%sphi
267  zeta => basis_set_a%zet
268  ! basis jkind
269  first_sgfb => basis_set_b%first_sgf
270  lb_max => basis_set_b%lmax
271  lb_min => basis_set_b%lmin
272  npgfb => basis_set_b%npgf
273  nsetb = basis_set_b%nset
274  nsgfb => basis_set_b%nsgf_set
275  rpgfb => basis_set_b%pgf_radius
276  set_radius_b => basis_set_b%set_radius
277  sphi_b => basis_set_b%sphi
278  zetb => basis_set_b%zet
279 
280  atom_a = atom_of_kind(iatom)
281  atom_b = atom_of_kind(jatom)
282 
283  ra(:) = particle_set(iatom)%r(:) - rpoint(:)
284  rac(:) = pbc(ra(:), cell)
285  rbc(:) = rac(:) + rab(:)
286  dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
287 
288  IF (iatom <= jatom) THEN
289  irow = iatom
290  icol = jatom
291  trans = .false.
292  ELSE
293  irow = jatom
294  icol = iatom
295  trans = .true.
296  END IF
297 
298  fdir = 2.0_dp
299  IF (iatom == jatom .AND. dab < 1.e-10_dp) fdir = 1.0_dp
300 
301  ! density matrix
302  NULLIFY (p_block_a)
303  CALL dbcsr_get_block_p(ec_env%matrix_p(1, 1)%matrix, irow, icol, p_block_a, found)
304  IF (.NOT. found) cycle
305  IF (SIZE(ec_env%matrix_p, 1) > 1) THEN
306  NULLIFY (p_block_b)
307  CALL dbcsr_get_block_p(ec_env%matrix_p(2, 1)%matrix, irow, icol, p_block_b, found)
308  cpassert(found)
309  END IF
310  forcea = 0.0_dp
311  forceb = 0.0_dp
312 
313  DO iset = 1, nseta
314  ncoa = npgfa(iset)*ncoset(la_max(iset))
315  sgfa = first_sgfa(1, iset)
316  DO jset = 1, nsetb
317  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
318  ncob = npgfb(jset)*ncoset(lb_max(jset))
319  sgfb = first_sgfb(1, jset)
320  ! Calculate the primitive integrals (da|O|b) and (a|O|db)
321  ldab = max(ncoa, ncob)
322  ALLOCATE (work(ldab, ldab), pmat(ncoa, ncob))
323  ! Decontract P matrix block
324  pmat = 0.0_dp
325  DO i = 1, SIZE(ec_env%matrix_p, 1)
326  IF (i == 1) THEN
327  pblock => p_block_a
328  ELSE
329  pblock => p_block_b
330  END IF
331  IF (.NOT. ASSOCIATED(pblock)) cycle
332  IF (trans) THEN
333  CALL dgemm("N", "T", ncoa, nsgfb(jset), nsgfa(iset), &
334  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
335  pblock(sgfb, sgfa), SIZE(pblock, 1), &
336  0.0_dp, work(1, 1), ldab)
337  ELSE
338  CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
339  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
340  pblock(sgfa, sgfb), SIZE(pblock, 1), &
341  0.0_dp, work(1, 1), ldab)
342  END IF
343  CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
344  1.0_dp, work(1, 1), ldab, &
345  sphi_b(1, sgfb), SIZE(sphi_b, 1), &
346  1.0_dp, pmat(1, 1), ncoa)
347  END DO
348 
349  CALL dipole_force(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
350  lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
351  1, rac, rbc, pmat, forcea, forceb)
352 
353  DEALLOCATE (work, pmat)
354  END DO
355  END DO
356 
357  DO idir = 1, 3
358  force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) &
359  + fdir*fieldpol(idir)*forcea(idir, 1:3)
360  force(jkind)%efield(1:3, atom_b) = force(jkind)%efield(1:3, atom_b) &
361  + fdir*fieldpol(idir)*forceb(idir, 1:3)
362  END DO
363 
364  END DO
365  CALL neighbor_list_iterator_release(nl_iterator)
366  DEALLOCATE (basis_set_list)
367  END IF
368 
369  CALL timestop(handle)
370 
371  END SUBROUTINE ec_efield_mo_derivatives
372 
373 END MODULE ec_efield_local
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition: ai_moments.F:17
subroutine, public dipole_force(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, pab, forcea, forceb)
This returns the derivative of the dipole integrals [a|x|b], with respect to the position of the prim...
Definition: ai_moments.F:1667
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.
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...
Calculates the energy contribution and the mo_derivative of a static electric field (nonperiodic)
subroutine, public ec_efield_local_operator(qs_env, ec_env, calculate_forces)
...
subroutine, public ec_efield_integrals(qs_env, ec_env, rpoint)
...
Types needed for a for a Energy Correction.
Definition: ec_env_types.F:14
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
Define the data structure for the particle information.
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.
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition: qs_moments.F:14
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Definition: qs_moments.F:558
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
type for berry phase efield matrices. At the moment only used for cosmat and sinmat
subroutine, public set_efield_matrices(efield, sinmat, cosmat, dipmat)
...
subroutine, public init_efield_matrices(efield)
...