Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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_api, ONLY: dbcsr_add,&
20 : dbcsr_copy,&
21 : dbcsr_p_type,&
22 : dbcsr_set
23 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
24 : dbcsr_deallocate_matrix_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 : twopi
32 : USE particle_types, ONLY: particle_type
33 : USE physcon, ONLY: c_light_au
34 : USE qs_energy_types, ONLY: qs_energy_type
35 : USE qs_environment_types, ONLY: get_qs_env,&
36 : qs_environment_type
37 : USE qs_force_types, ONLY: qs_force_type
38 : USE qs_kind_types, ONLY: get_qs_kind,&
39 : qs_kind_type
40 : USE qs_moments, ONLY: build_local_moment_matrix
41 : #include "./base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 :
45 : PRIVATE
46 :
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'efield_utils'
48 :
49 : ! *** Public subroutines ***
50 :
51 : PUBLIC :: efield_potential_lengh_gauge, &
52 : calculate_ecore_efield, &
53 : make_field
54 :
55 : CONTAINS
56 :
57 : ! **************************************************************************************************
58 : !> \brief Replace the original implementation of the electric-electronic
59 : !> interaction in the length gauge. This calculation is no longer done in
60 : !> the grid but using matrices to match the velocity gauge implementation.
61 : !> Note: The energy is stored in energy%core and computed later on.
62 : !> \param qs_env ...
63 : !> \author Guillaume Le Breton (02.23)
64 : ! **************************************************************************************************
65 :
66 214 : SUBROUTINE efield_potential_lengh_gauge(qs_env)
67 :
68 : TYPE(qs_environment_type), POINTER :: qs_env
69 :
70 : CHARACTER(len=*), PARAMETER :: routineN = 'efield_potential_lengh_gauge'
71 :
72 : INTEGER :: handle, i, image
73 : REAL(kind=dp) :: field(3)
74 214 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments
75 214 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h
76 : TYPE(dft_control_type), POINTER :: dft_control
77 :
78 214 : NULLIFY (dft_control)
79 214 : CALL timeset(routineN, handle)
80 :
81 : CALL get_qs_env(qs_env, &
82 : dft_control=dft_control, &
83 : matrix_h_kp=matrix_h, &
84 214 : matrix_s=matrix_s)
85 :
86 214 : NULLIFY (moments)
87 214 : CALL dbcsr_allocate_matrix_set(moments, 3)
88 856 : DO i = 1, 3
89 642 : ALLOCATE (moments(i)%matrix)
90 642 : CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
91 856 : CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
92 : END DO
93 :
94 214 : CALL build_local_moment_matrix(qs_env, moments, 1)
95 :
96 214 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
97 :
98 856 : DO i = 1, 3
99 1498 : DO image = 1, dft_control%nimages
100 1284 : CALL dbcsr_add(matrix_h(1, image)%matrix, moments(i)%matrix, 1.0_dp, field(i))
101 : END DO
102 : END DO
103 :
104 214 : CALL dbcsr_deallocate_matrix_set(moments)
105 :
106 214 : CALL timestop(handle)
107 :
108 214 : END SUBROUTINE efield_potential_lengh_gauge
109 :
110 : ! **************************************************************************************************
111 : !> \brief computes the amplitude of the efield within a given envelop
112 : !> \param dft_control ...
113 : !> \param field ...
114 : !> \param sim_step ...
115 : !> \param sim_time ...
116 : !> \author Florian Schiffmann (02.09)
117 : ! **************************************************************************************************
118 :
119 1204 : SUBROUTINE make_field(dft_control, field, sim_step, sim_time)
120 : TYPE(dft_control_type), INTENT(IN) :: dft_control
121 : REAL(dp), INTENT(OUT) :: field(3)
122 : INTEGER, INTENT(IN) :: sim_step
123 : REAL(KIND=dp), INTENT(IN) :: sim_time
124 :
125 : INTEGER :: i, lower, nfield, upper
126 : REAL(dp) :: E_0, env, nu, pol(3), strength
127 : REAL(KIND=dp) :: dt
128 : TYPE(efield_type), POINTER :: efield
129 :
130 1204 : field = 0._dp
131 1204 : nfield = SIZE(dft_control%efield_fields)
132 2408 : DO i = 1, nfield
133 1204 : efield => dft_control%efield_fields(i)%efield
134 1204 : nu = 0.0_dp
135 1204 : IF (.NOT. efield%envelop_id == custom_env .AND. efield%wavelength > EPSILON(0.0_dp)) THEN
136 1090 : nu = c_light_au/(efield%wavelength) !in case of a custom efield we do not need nu
137 : END IF
138 1204 : E_0 = efield%amplitude
139 :
140 4816 : IF (DOT_PRODUCT(efield%polarisation, efield%polarisation) == 0) THEN
141 0 : pol(:) = 1.0_dp/3.0_dp
142 : ELSE
143 8428 : pol(:) = efield%polarisation(:)/(SQRT(DOT_PRODUCT(efield%polarisation, efield%polarisation)))
144 : END IF
145 :
146 1204 : SELECT CASE (efield%envelop_id)
147 : CASE (constant_env)
148 212 : IF (sim_step >= efield%envelop_i_vars(1) .AND. &
149 118 : (sim_step <= efield%envelop_i_vars(2) .OR. efield%envelop_i_vars(2) < 0)) THEN
150 : field = field + E_0*COS(sim_time*nu*twopi + &
151 530 : efield%phase_offset*pi)*pol(:)
152 : END IF
153 : CASE (ramp_env)
154 118 : IF (sim_step >= efield%envelop_i_vars(1) .AND. sim_step <= efield%envelop_i_vars(2)) &
155 52 : strength = E_0*(sim_step - efield%envelop_i_vars(1))/(efield%envelop_i_vars(2) - efield%envelop_i_vars(1))
156 118 : IF (sim_step >= efield%envelop_i_vars(3) .AND. sim_step <= efield%envelop_i_vars(4)) &
157 60 : strength = E_0*(efield%envelop_i_vars(4) - sim_step)/(efield%envelop_i_vars(4) - efield%envelop_i_vars(3))
158 118 : IF (sim_step > efield%envelop_i_vars(4) .AND. efield%envelop_i_vars(4) > 0) strength = 0.0_dp
159 118 : IF (sim_step <= efield%envelop_i_vars(1)) strength = 0.0_dp
160 : field = field + strength*COS(sim_time*nu*twopi + &
161 472 : efield%phase_offset*pi)*pol(:)
162 : CASE (gaussian_env)
163 760 : env = EXP(-0.5_dp*((sim_time - efield%envelop_r_vars(1))/efield%envelop_r_vars(2))**2.0_dp)
164 : field = field + E_0*env*COS(sim_time*nu*twopi + &
165 3040 : efield%phase_offset*pi)*pol(:)
166 : CASE (custom_env)
167 114 : dt = efield%envelop_r_vars(1)
168 114 : IF (sim_time < (SIZE(efield%envelop_r_vars) - 2)*dt) THEN
169 : !make a linear interpolation between the two next points
170 62 : lower = FLOOR(sim_time/dt)
171 62 : upper = lower + 1
172 : strength = (efield%envelop_r_vars(lower + 2)*(upper*dt - sim_time) &
173 62 : + efield%envelop_r_vars(upper + 2)*(sim_time - lower*dt))/dt
174 : ELSE
175 : strength = 0.0_dp
176 : END IF
177 1660 : field = field + strength*pol(:)
178 : CASE default
179 : END SELECT
180 : END DO
181 :
182 1204 : END SUBROUTINE make_field
183 :
184 : ! **************************************************************************************************
185 : !> \brief Computes the force and the energy due to a efield on the cores
186 : !> Note: In the velocity gauge, the energy term is not added because
187 : !> it would lead to an unbalanced energy (center of negative charge not
188 : !> involved in the electric energy in this gauge).
189 : !> \param qs_env ...
190 : !> \param calculate_forces ...
191 : !> \author Florian Schiffmann (02.09)
192 : ! **************************************************************************************************
193 :
194 18797 : SUBROUTINE calculate_ecore_efield(qs_env, calculate_forces)
195 : TYPE(qs_environment_type), POINTER :: qs_env
196 : LOGICAL, OPTIONAL :: calculate_forces
197 :
198 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ecore_efield'
199 :
200 : INTEGER :: atom_a, handle, iatom, ikind, natom, &
201 : nkind
202 18797 : INTEGER, DIMENSION(:), POINTER :: list
203 : LOGICAL :: my_force
204 : REAL(KIND=dp) :: efield_ener, zeff
205 : REAL(KIND=dp), DIMENSION(3) :: field, r
206 18797 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
207 : TYPE(cell_type), POINTER :: cell
208 : TYPE(dft_control_type), POINTER :: dft_control
209 18797 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
210 : TYPE(qs_energy_type), POINTER :: energy
211 18797 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
212 18797 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
213 :
214 18797 : NULLIFY (dft_control)
215 18797 : CALL timeset(routineN, handle)
216 18797 : CALL get_qs_env(qs_env, dft_control=dft_control)
217 18797 : IF (dft_control%apply_efield_field .OR. dft_control%apply_vector_potential) THEN
218 322 : my_force = .FALSE.
219 322 : IF (PRESENT(calculate_forces)) my_force = calculate_forces
220 :
221 : CALL get_qs_env(qs_env=qs_env, &
222 : atomic_kind_set=atomic_kind_set, &
223 : qs_kind_set=qs_kind_set, &
224 : energy=energy, &
225 : particle_set=particle_set, &
226 322 : cell=cell)
227 322 : efield_ener = 0.0_dp
228 322 : nkind = SIZE(atomic_kind_set)
229 322 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
230 :
231 722 : DO ikind = 1, SIZE(atomic_kind_set)
232 400 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list, natom=natom)
233 400 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
234 :
235 400 : natom = SIZE(list)
236 1444 : DO iatom = 1, natom
237 722 : IF (dft_control%apply_efield_field) THEN
238 510 : atom_a = list(iatom)
239 510 : r(:) = pbc(particle_set(atom_a)%r(:), cell)
240 2040 : efield_ener = efield_ener - zeff*DOT_PRODUCT(r, field)
241 : END IF
242 1122 : IF (my_force) THEN
243 408 : CALL get_qs_env(qs_env=qs_env, force=force)
244 1632 : force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) - field*zeff
245 : END IF
246 : END DO
247 :
248 : END DO
249 322 : IF (dft_control%apply_efield_field) energy%efield_core = efield_ener
250 : END IF
251 18797 : CALL timestop(handle)
252 18797 : END SUBROUTINE calculate_ecore_efield
253 : END MODULE efield_utils
|
