Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : 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 214 : 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 214 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments
73 214 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h
74 : TYPE(dft_control_type), POINTER :: dft_control
75 :
76 214 : NULLIFY (dft_control)
77 214 : CALL timeset(routineN, handle)
78 :
79 : CALL get_qs_env(qs_env, &
80 : dft_control=dft_control, &
81 : matrix_h_kp=matrix_h, &
82 214 : matrix_s=matrix_s)
83 :
84 214 : NULLIFY (moments)
85 214 : CALL dbcsr_allocate_matrix_set(moments, 3)
86 856 : DO i = 1, 3
87 642 : ALLOCATE (moments(i)%matrix)
88 642 : CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
89 856 : CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
90 : END DO
91 :
92 214 : CALL build_local_moment_matrix(qs_env, moments, 1)
93 :
94 214 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
95 :
96 856 : DO i = 1, 3
97 1498 : DO image = 1, dft_control%nimages
98 1284 : CALL dbcsr_add(matrix_h(1, image)%matrix, moments(i)%matrix, 1.0_dp, field(i))
99 : END DO
100 : END DO
101 :
102 214 : CALL dbcsr_deallocate_matrix_set(moments)
103 :
104 214 : CALL timestop(handle)
105 :
106 214 : 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 986 : 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 986 : c = 137.03599962875_dp
129 986 : field = 0._dp
130 986 : nu = 0.0_dp
131 986 : nfield = SIZE(dft_control%efield_fields)
132 1972 : DO i = 1, nfield
133 986 : efield => dft_control%efield_fields(i)%efield
134 986 : IF (.NOT. efield%envelop_id == custom_env .AND. efield%wavelength > EPSILON(0.0_dp)) nu = c/(efield%wavelength) !in case of a custom efield we do not need nu
135 986 : strength = SQRT(efield%strength/(3.50944_dp*10.0_dp**16))
136 3944 : IF (DOT_PRODUCT(efield%polarisation, efield%polarisation) == 0) THEN
137 0 : pol(:) = 1.0_dp/3.0_dp
138 : ELSE
139 6902 : pol(:) = efield%polarisation(:)/(SQRT(DOT_PRODUCT(efield%polarisation, efield%polarisation)))
140 : END IF
141 1972 : IF (efield%envelop_id == constant_env) THEN
142 212 : IF (sim_step >= efield%envelop_i_vars(1) .AND. &
143 : (sim_step <= efield%envelop_i_vars(2) .OR. efield%envelop_i_vars(2) < 0)) THEN
144 : field = field + strength*COS(sim_time*nu*2.0_dp*pi + &
145 530 : efield%phase_offset*pi)*pol(:)
146 : END IF
147 774 : ELSE IF (efield%envelop_id == ramp_env) THEN
148 118 : IF (sim_step >= efield%envelop_i_vars(1) .AND. sim_step <= efield%envelop_i_vars(2)) &
149 52 : strength = strength*(sim_step - efield%envelop_i_vars(1))/(efield%envelop_i_vars(2) - efield%envelop_i_vars(1))
150 118 : IF (sim_step >= efield%envelop_i_vars(3) .AND. sim_step <= efield%envelop_i_vars(4)) &
151 60 : strength = strength*(efield%envelop_i_vars(4) - sim_step)/(efield%envelop_i_vars(4) - efield%envelop_i_vars(3))
152 118 : IF (sim_step > efield%envelop_i_vars(4) .AND. efield%envelop_i_vars(4) > 0) strength = 0.0_dp
153 118 : IF (sim_step <= efield%envelop_i_vars(1)) strength = 0.0_dp
154 : field = field + strength*COS(sim_time*nu*2.0_dp*pi + &
155 472 : efield%phase_offset*pi)*pol(:)
156 656 : ELSE IF (efield%envelop_id == gaussian_env) THEN
157 542 : env = EXP(-0.5_dp*((sim_time - efield%envelop_r_vars(1))/efield%envelop_r_vars(2))**2.0_dp)
158 : field = field + strength*env*COS(sim_time*nu*2.0_dp*pi + &
159 2168 : efield%phase_offset*pi)*pol(:)
160 114 : ELSE IF (efield%envelop_id == custom_env) THEN
161 114 : dt = efield%envelop_r_vars(1)
162 114 : IF (sim_time < (SIZE(efield%envelop_r_vars) - 2)*dt) THEN
163 : !make a linear interpolation between the two next points
164 62 : lower = FLOOR(sim_time/dt)
165 62 : upper = lower + 1
166 62 : strength = (efield%envelop_r_vars(lower + 2)*(upper*dt - sim_time) + efield%envelop_r_vars(upper + 2)*(sim_time - lower*dt))/dt
167 : ELSE
168 : strength = 0.0_dp
169 : END IF
170 456 : field = field + strength*pol(:)
171 : END IF
172 : END DO
173 :
174 986 : END SUBROUTINE make_field
175 :
176 : ! **************************************************************************************************
177 : !> \brief Computes the force and the energy due to a efield on the cores
178 : !> Note: In the velocity gauge, the energy term is not added because
179 : !> it would lead to an unbalanced energy (center of negative charge not
180 : !> involved in the electric energy in this gauge).
181 : !> \param qs_env ...
182 : !> \param calculate_forces ...
183 : !> \author Florian Schiffmann (02.09)
184 : ! **************************************************************************************************
185 :
186 16969 : SUBROUTINE calculate_ecore_efield(qs_env, calculate_forces)
187 : TYPE(qs_environment_type), POINTER :: qs_env
188 : LOGICAL, OPTIONAL :: calculate_forces
189 :
190 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ecore_efield'
191 :
192 : INTEGER :: atom_a, handle, iatom, ikind, natom, &
193 : nkind
194 16969 : INTEGER, DIMENSION(:), POINTER :: list
195 : LOGICAL :: my_force
196 : REAL(KIND=dp) :: efield_ener, zeff
197 : REAL(KIND=dp), DIMENSION(3) :: field, r
198 16969 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
199 : TYPE(cell_type), POINTER :: cell
200 : TYPE(dft_control_type), POINTER :: dft_control
201 16969 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
202 : TYPE(qs_energy_type), POINTER :: energy
203 16969 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
204 16969 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
205 :
206 16969 : NULLIFY (dft_control)
207 16969 : CALL timeset(routineN, handle)
208 16969 : CALL get_qs_env(qs_env, dft_control=dft_control)
209 16969 : IF (dft_control%apply_efield_field .OR. dft_control%apply_vector_potential) THEN
210 320 : my_force = .FALSE.
211 320 : IF (PRESENT(calculate_forces)) my_force = calculate_forces
212 :
213 : CALL get_qs_env(qs_env=qs_env, &
214 : atomic_kind_set=atomic_kind_set, &
215 : qs_kind_set=qs_kind_set, &
216 : energy=energy, &
217 : particle_set=particle_set, &
218 320 : cell=cell)
219 320 : efield_ener = 0.0_dp
220 320 : nkind = SIZE(atomic_kind_set)
221 320 : CALL make_field(dft_control, field, qs_env%sim_step, qs_env%sim_time)
222 :
223 718 : DO ikind = 1, SIZE(atomic_kind_set)
224 398 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list, natom=natom)
225 398 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
226 :
227 398 : natom = SIZE(list)
228 1436 : DO iatom = 1, natom
229 718 : IF (dft_control%apply_efield_field) THEN
230 506 : atom_a = list(iatom)
231 506 : r(:) = pbc(particle_set(atom_a)%r(:), cell)
232 2024 : efield_ener = efield_ener - zeff*DOT_PRODUCT(r, field)
233 : END IF
234 1116 : IF (my_force) THEN
235 408 : CALL get_qs_env(qs_env=qs_env, force=force)
236 1632 : force(ikind)%efield(:, iatom) = force(ikind)%efield(:, iatom) - field*zeff
237 : END IF
238 : END DO
239 :
240 : END DO
241 320 : IF (dft_control%apply_efield_field) energy%efield_core = efield_ener
242 : END IF
243 16969 : CALL timestop(handle)
244 16969 : END SUBROUTINE calculate_ecore_efield
245 : END MODULE efield_utils
|
