(git:b195825)
efield_utils.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 all routins needed for a nonperiodic electric field
10 ! **************************************************************************************************
11 
12 MODULE efield_utils
13  USE atomic_kind_types, ONLY: atomic_kind_type,&
14  get_atomic_kind
15  USE cell_types, ONLY: cell_type,&
16  pbc
17  USE cp_control_types, ONLY: dft_control_type,&
18  efield_type
19  USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
20  dbcsr_deallocate_matrix_set
21  USE dbcsr_api, ONLY: dbcsr_add,&
22  dbcsr_copy,&
23  dbcsr_p_type,&
24  dbcsr_set
25  USE input_constants, ONLY: constant_env,&
26  custom_env,&
27  gaussian_env,&
28  ramp_env
29  USE kinds, ONLY: dp
30  USE mathconstants, ONLY: pi
31  USE particle_types, ONLY: particle_type
32  USE qs_energy_types, ONLY: qs_energy_type
33  USE qs_environment_types, ONLY: get_qs_env,&
34  qs_environment_type
35  USE qs_force_types, ONLY: qs_force_type
36  USE qs_kind_types, ONLY: get_qs_kind,&
37  qs_kind_type
38  USE qs_moments, ONLY: build_local_moment_matrix
39 #include "./base/base_uses.f90"
40 
41  IMPLICIT NONE
42 
43  PRIVATE
44 
45  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'efield_utils'
46 
47 ! *** Public subroutines ***
48 
49  PUBLIC :: efield_potential_lengh_gauge, &
50  calculate_ecore_efield, &
51  make_field
52 
53 CONTAINS
54 
55 ! **************************************************************************************************
56 !> \brief Replace the original implementation of the electric-electronic
57 !> interaction in the length gauge. This calculation is no longer done in
58 !> the grid but using matrices to match the velocity gauge implementation.
59 !> Note: The energy is stored in energy%core and computed later on.
60 !> \param qs_env ...
61 !> \author Guillaume Le Breton (02.23)
62 ! **************************************************************************************************
63 
64  SUBROUTINE efield_potential_lengh_gauge(qs_env)
65 
66  TYPE(qs_environment_type), POINTER :: qs_env
67 
68  CHARACTER(len=*), PARAMETER :: routinen = 'efield_potential_lengh_gauge'
69 
70  INTEGER :: handle, i, image
71  REAL(kind=dp) :: field(3)
72  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments
73  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h
74  TYPE(dft_control_type), POINTER :: dft_control
75 
76  NULLIFY (dft_control)
77  CALL timeset(routinen, handle)
78 
79  CALL get_qs_env(qs_env, &
80  dft_control=dft_control, &
81  matrix_h_kp=matrix_h, &
82  matrix_s=matrix_s)
83 
84  NULLIFY (moments)
85  CALL dbcsr_allocate_matrix_set(moments, 3)
86  DO i = 1, 3
87  ALLOCATE (moments(i)%matrix)
88  CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
89  CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
90  END DO
91 
92  CALL build_local_moment_matrix(qs_env, moments, 1)
93 
94  CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
95 
96  DO i = 1, 3
97  DO image = 1, dft_control%nimages
98  CALL dbcsr_add(matrix_h(1, image)%matrix, moments(i)%matrix, 1.0_dp, field(i))
99  END DO
100  END DO
101 
102  CALL dbcsr_deallocate_matrix_set(moments)
103 
104  CALL timestop(handle)
105 
106  END SUBROUTINE efield_potential_lengh_gauge
107 
108 ! **************************************************************************************************
109 !> \brief computes the amplitude of the efield within a given envelop
110 !> \param dft_control ...
111 !> \param field ...
112 !> \param sim_step ...
113 !> \param sim_time ...
114 !> \author Florian Schiffmann (02.09)
115 ! **************************************************************************************************
116 
117  SUBROUTINE make_field(dft_control, field, sim_step, sim_time)
118  TYPE(dft_control_type), INTENT(IN) :: dft_control
119  REAL(dp), INTENT(OUT) :: field(3)
120  INTEGER, INTENT(IN) :: sim_step
121  REAL(kind=dp), INTENT(IN) :: sim_time
122 
123  INTEGER :: i, lower, nfield, upper
124  REAL(dp) :: c, env, nu, pol(3), strength
125  REAL(kind=dp) :: dt
126  TYPE(efield_type), POINTER :: efield
127 
128  c = 137.03599962875_dp
129  field = 0._dp
130  nfield = SIZE(dft_control%efield_fields)
131  DO i = 1, nfield
132  efield => dft_control%efield_fields(i)%efield
133  IF (.NOT. efield%envelop_id == custom_env) nu = c/(efield%wavelength) !in case of a custom efield we do not need nu
134  strength = sqrt(efield%strength/(3.50944_dp*10.0_dp**16))
135  IF (dot_product(efield%polarisation, efield%polarisation) == 0) THEN
136  pol(:) = 1.0_dp/3.0_dp
137  ELSE
138  pol(:) = efield%polarisation(:)/(sqrt(dot_product(efield%polarisation, efield%polarisation)))
139  END IF
140  IF (efield%envelop_id == constant_env) THEN
141  IF (sim_step .GE. efield%envelop_i_vars(1) .AND. &
142  (sim_step .LE. efield%envelop_i_vars(2) .OR. efield%envelop_i_vars(2) .LT. 0)) THEN
143  field = field + strength*cos(sim_time*nu*2.0_dp*pi + &
144  efield%phase_offset*pi)*pol(:)
145  END IF
146  ELSE IF (efield%envelop_id == ramp_env) THEN
147  IF (sim_step .GE. efield%envelop_i_vars(1) .AND. sim_step .LE. efield%envelop_i_vars(2)) &
148  strength = strength*(sim_step - efield%envelop_i_vars(1))/(efield%envelop_i_vars(2) - efield%envelop_i_vars(1))
149  IF (sim_step .GE. efield%envelop_i_vars(3) .AND. sim_step .LE. efield%envelop_i_vars(4)) &
150  strength = strength*(efield%envelop_i_vars(4) - sim_step)/(efield%envelop_i_vars(4) - efield%envelop_i_vars(3))
151  IF (sim_step .GT. efield%envelop_i_vars(4) .AND. efield%envelop_i_vars(4) .GT. 0) strength = 0.0_dp
152  IF (sim_step .LE. efield%envelop_i_vars(1)) strength = 0.0_dp
153  field = field + strength*cos(sim_time*nu*2.0_dp*pi + &
154  efield%phase_offset*pi)*pol(:)
155  ELSE IF (efield%envelop_id == gaussian_env) THEN
156  env = exp(-0.5_dp*((sim_time - efield%envelop_r_vars(1))/efield%envelop_r_vars(2))**2.0_dp)
157  field = field + strength*env*cos(sim_time*nu*2.0_dp*pi + &
158  efield%phase_offset*pi)*pol(:)
159  ELSE IF (efield%envelop_id == custom_env) THEN
160  dt = efield%envelop_r_vars(1)
161  IF (sim_time .LT. (SIZE(efield%envelop_r_vars) - 2)*dt) THEN
162  !make a linear interpolation between the two next points
163  lower = floor(sim_time/dt)
164  upper = lower + 1
165  strength = (efield%envelop_r_vars(lower + 2)*(upper*dt - sim_time) + efield%envelop_r_vars(upper + 2)*(sim_time - lower*dt))/dt
166  ELSE
167  strength = 0.0_dp
168  END IF
169  field = field + strength*pol(:)
170  END IF
171  END DO
172 
173  END SUBROUTINE make_field
174 
175 ! **************************************************************************************************
176 !> \brief Computes the force and the energy due to a efield on the cores
177 !> Note: In the velocity gauge, the energy term is not added because
178 !> it would lead to an unbalanced energy (center of negative charge not
179 !> involved in the electric energy in this gauge).
180 !> \param qs_env ...
181 !> \param calculate_forces ...
182 !> \author Florian Schiffmann (02.09)
183 ! **************************************************************************************************
184 
185  SUBROUTINE calculate_ecore_efield(qs_env, calculate_forces)
186  TYPE(qs_environment_type), POINTER :: qs_env
187  LOGICAL, OPTIONAL :: calculate_forces
188 
189  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_ecore_efield'
190 
191  INTEGER :: atom_a, handle, iatom, ikind, natom, &
192  nkind
193  INTEGER, DIMENSION(:), POINTER :: list
194  LOGICAL :: my_force
195  REAL(kind=dp) :: efield_ener, zeff
196  REAL(kind=dp), DIMENSION(3) :: field, r
197  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
198  TYPE(cell_type), POINTER :: cell
199  TYPE(dft_control_type), POINTER :: dft_control
200  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
201  TYPE(qs_energy_type), POINTER :: energy
202  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
203  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
204 
205  NULLIFY (dft_control)
206  CALL timeset(routinen, handle)
207  CALL get_qs_env(qs_env, dft_control=dft_control)
208  IF (dft_control%apply_efield_field .OR. dft_control%apply_vector_potential) THEN
209  my_force = .false.
210  IF (PRESENT(calculate_forces)) my_force = calculate_forces
211 
212  CALL get_qs_env(qs_env=qs_env, &
213  atomic_kind_set=atomic_kind_set, &
214  qs_kind_set=qs_kind_set, &
215  energy=energy, &
216  particle_set=particle_set, &
217  cell=cell)
218  efield_ener = 0.0_dp
219  nkind = SIZE(atomic_kind_set)
220  CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
221 
222  DO ikind = 1, SIZE(atomic_kind_set)
223  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list, natom=natom)
224  CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
225 
226  natom = SIZE(list)
227  DO iatom = 1, natom
228  IF (dft_control%apply_efield_field) THEN
229  atom_a = list(iatom)
230  r(:) = pbc(particle_set(atom_a)%r(:), cell)
231  efield_ener = efield_ener - zeff*dot_product(r, field)
232  END IF
233  IF (my_force) THEN
234  CALL get_qs_env(qs_env=qs_env, force=force)
235  force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) - field*zeff
236  END IF
237 ! END IF
238  END DO
239 
240  END DO
241  IF (dft_control%apply_efield_field) energy%efield_core = efield_ener
242 ! energy%efield_core = efield_ener
243  END IF
244  CALL timestop(handle)
245  END SUBROUTINE calculate_ecore_efield
246 END MODULE efield_utils
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
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.
Definition: atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
efield_utils
all routins needed for a nonperiodic electric field
Definition: efield_utils.F:12
efield_utils::make_field
subroutine, public make_field(dft_control, field, sim_step, sim_time)
computes the amplitude of the efield within a given envelop
Definition: efield_utils.F:118
efield_utils::calculate_ecore_efield
subroutine, public calculate_ecore_efield(qs_env, calculate_forces)
Computes the force and the energy due to a efield on the cores Note: In the velocity gauge,...
Definition: efield_utils.F:186
efield_utils::efield_potential_lengh_gauge
subroutine, public efield_potential_lengh_gauge(qs_env)
Replace the original implementation of the electric-electronic interaction in the length gauge....
Definition: efield_utils.F:65
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::ramp_env
integer, parameter, public ramp_env
Definition: input_constants.F:931
input_constants::constant_env
integer, parameter, public constant_env
Definition: input_constants.F:931
input_constants::gaussian_env
integer, parameter, public gaussian_env
Definition: input_constants.F:931
input_constants::custom_env
integer, parameter, public custom_env
Definition: input_constants.F:931
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
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.
Definition: qs_environment_types.F:502
qs_force_types
Definition: qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
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.
Definition: qs_kind_types.F:451
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition: qs_moments.F:14
qs_moments::build_local_moment_matrix
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Definition: qs_moments.F:558