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 : !> \par History
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE fist_efield_methods
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_result_methods, ONLY: cp_results_erase,&
18 : put_results
19 : USE cp_result_types, ONLY: cp_result_type
20 : USE fist_efield_types, ONLY: fist_efield_type
21 : USE fist_environment_types, ONLY: fist_env_get,&
22 : fist_environment_type
23 : USE input_section_types, ONLY: section_get_ival,&
24 : section_vals_type,&
25 : section_vals_val_get
26 : USE kinds, ONLY: default_string_length,&
27 : dp
28 : USE mathconstants, ONLY: twopi,&
29 : z_one,&
30 : z_zero
31 : USE moments_utils, ONLY: get_reference_point
32 : USE particle_types, ONLY: particle_type
33 : USE physcon, ONLY: debye
34 : #include "./base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_efield_methods'
39 :
40 : PRIVATE
41 :
42 : PUBLIC :: fist_dipole, fist_efield_energy_force
43 :
44 : ! **************************************************************************************************
45 :
46 : CONTAINS
47 :
48 : ! **************************************************************************************************
49 : !> \brief ...
50 : !> \param qenergy ...
51 : !> \param qforce ...
52 : !> \param qpv ...
53 : !> \param atomic_kind_set ...
54 : !> \param particle_set ...
55 : !> \param cell ...
56 : !> \param efield ...
57 : !> \param use_virial ...
58 : !> \param iunit ...
59 : !> \param charges ...
60 : ! **************************************************************************************************
61 116 : SUBROUTINE fist_efield_energy_force(qenergy, qforce, qpv, atomic_kind_set, particle_set, cell, &
62 : efield, use_virial, iunit, charges)
63 : REAL(KIND=dp), INTENT(OUT) :: qenergy
64 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: qforce
65 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: qpv
66 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
67 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
68 : TYPE(cell_type), POINTER :: cell
69 : TYPE(fist_efield_type), POINTER :: efield
70 : LOGICAL, INTENT(IN), OPTIONAL :: use_virial
71 : INTEGER, INTENT(IN), OPTIONAL :: iunit
72 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
73 :
74 : COMPLEX(KIND=dp) :: zeta
75 : COMPLEX(KIND=dp), DIMENSION(3) :: ggamma
76 : INTEGER :: i, ii, iparticle_kind, iw, j
77 116 : INTEGER, DIMENSION(:), POINTER :: atom_list
78 : LOGICAL :: use_charges, virial
79 : REAL(KIND=dp) :: q, theta
80 : REAL(KIND=dp), DIMENSION(3) :: ci, dfilter, di, dipole, fieldpol, fq, &
81 : gvec, ria, tmp
82 : TYPE(atomic_kind_type), POINTER :: atomic_kind
83 :
84 116 : qenergy = 0.0_dp
85 1508 : qforce = 0.0_dp
86 116 : qpv = 0.0_dp
87 :
88 116 : use_charges = .FALSE.
89 116 : IF (PRESENT(charges)) THEN
90 116 : IF (ASSOCIATED(charges)) use_charges = .TRUE.
91 : END IF
92 :
93 : IF (PRESENT(iunit)) THEN
94 116 : iw = iunit
95 : ELSE
96 116 : iw = -1
97 : END IF
98 :
99 116 : IF (PRESENT(use_virial)) THEN
100 116 : virial = use_virial
101 : ELSE
102 : virial = .FALSE.
103 : END IF
104 :
105 464 : fieldpol = efield%polarisation
106 812 : fieldpol = fieldpol/SQRT(DOT_PRODUCT(fieldpol, fieldpol))
107 464 : fieldpol = -fieldpol*efield%strength
108 :
109 464 : dfilter = efield%dfilter
110 :
111 116 : dipole = 0.0_dp
112 464 : ggamma = z_one
113 348 : DO iparticle_kind = 1, SIZE(atomic_kind_set)
114 232 : atomic_kind => atomic_kind_set(iparticle_kind)
115 232 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q, atom_list=atom_list)
116 : ! TODO parallelization over atoms (local_particles)
117 696 : DO i = 1, SIZE(atom_list)
118 348 : ii = atom_list(i)
119 1392 : ria = particle_set(ii)%r(:)
120 1392 : ria = pbc(ria, cell)
121 348 : IF (use_charges) q = charges(ii)
122 1392 : DO j = 1, 3
123 4176 : gvec = twopi*cell%h_inv(j, :)
124 4176 : theta = SUM(ria(:)*gvec(:))
125 1044 : zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(-q)
126 1392 : ggamma(j) = ggamma(j)*zeta
127 : END DO
128 1624 : qforce(1:3, ii) = q
129 : END DO
130 : END DO
131 :
132 464 : IF (ALL(REAL(ggamma, KIND=dp) /= 0.0_dp)) THEN
133 464 : tmp = AIMAG(ggamma)/REAL(ggamma, KIND=dp)
134 464 : ci = ATAN(tmp)
135 1856 : dipole = MATMUL(cell%hmat, ci)/twopi
136 : END IF
137 :
138 116 : IF (efield%displacement) THEN
139 : ! E = (omega/8Pi)(D - 4Pi*P)^2
140 152 : di = dipole/cell%deth
141 152 : DO i = 1, 3
142 114 : theta = fieldpol(i) - 2._dp*twopi*di(i)
143 114 : qenergy = qenergy + dfilter(i)*theta**2
144 152 : fq(i) = -dfilter(i)*theta
145 : END DO
146 38 : qenergy = 0.25_dp*cell%deth/twopi*qenergy
147 152 : DO i = 1, SIZE(qforce, 2)
148 494 : qforce(1:3, i) = fq(1:3)*qforce(1:3, i)
149 : END DO
150 : ELSE
151 : ! E = -omega*E*P
152 312 : qenergy = SUM(fieldpol*dipole)
153 312 : DO i = 1, SIZE(qforce, 2)
154 1014 : qforce(1:3, i) = fieldpol(1:3)*qforce(1:3, i)
155 : END DO
156 : END IF
157 :
158 116 : IF (virial) THEN
159 6 : DO iparticle_kind = 1, SIZE(atomic_kind_set)
160 4 : atomic_kind => atomic_kind_set(iparticle_kind)
161 4 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list)
162 12 : DO i = 1, SIZE(atom_list)
163 6 : ii = atom_list(i)
164 24 : ria = particle_set(ii)%r(:)
165 24 : ria = pbc(ria, cell)
166 28 : DO j = 1, 3
167 78 : qpv(j, 1:3) = qpv(j, 1:3) + qforce(j, ii)*ria(1:3)
168 : END DO
169 : END DO
170 : END DO
171 : ! Stress tensor for constant D needs further investigation
172 2 : IF (efield%displacement) THEN
173 0 : CPABORT("Stress Tensor for constant D simulation is not working")
174 : END IF
175 : END IF
176 :
177 116 : END SUBROUTINE fist_efield_energy_force
178 : ! **************************************************************************************************
179 : !> \brief Evaluates the Dipole of a classical charge distribution(point-like)
180 : !> possibly using the berry phase formalism
181 : !> \param fist_env ...
182 : !> \param print_section ...
183 : !> \param atomic_kind_set ...
184 : !> \param particle_set ...
185 : !> \param cell ...
186 : !> \param unit_nr ...
187 : !> \param charges ...
188 : !> \par History
189 : !> [01.2006] created
190 : !> [12.2007] tlaino - University of Zurich - debug and extended
191 : !> \author Teodoro Laino
192 : ! **************************************************************************************************
193 42408 : SUBROUTINE fist_dipole(fist_env, print_section, atomic_kind_set, particle_set, &
194 : cell, unit_nr, charges)
195 : TYPE(fist_environment_type), POINTER :: fist_env
196 : TYPE(section_vals_type), POINTER :: print_section
197 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
198 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
199 : TYPE(cell_type), POINTER :: cell
200 : INTEGER, INTENT(IN) :: unit_nr
201 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
202 :
203 : CHARACTER(LEN=default_string_length) :: description, dipole_type
204 : COMPLEX(KIND=dp) :: dzeta, dzphase(3), zeta, zphase(3)
205 : COMPLEX(KIND=dp), DIMENSION(3) :: dggamma, ggamma
206 : INTEGER :: i, iparticle_kind, j, reference
207 21204 : INTEGER, DIMENSION(:), POINTER :: atom_list
208 : LOGICAL :: do_berry, use_charges
209 : REAL(KIND=dp) :: charge_tot, ci(3), dci(3), dipole(3), dipole_deriv(3), drcc(3), dria(3), &
210 : dtheta, gvec(3), q, rcc(3), ria(3), theta, tmp(3), via(3)
211 21204 : REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
212 : TYPE(atomic_kind_type), POINTER :: atomic_kind
213 : TYPE(cp_result_type), POINTER :: results
214 :
215 21204 : NULLIFY (atomic_kind)
216 : ! Reference point
217 42408 : reference = section_get_ival(print_section, keyword_name="DIPOLE%REFERENCE")
218 21204 : NULLIFY (ref_point)
219 21204 : description = '[DIPOLE]'
220 21204 : CALL section_vals_val_get(print_section, "DIPOLE%REF_POINT", r_vals=ref_point)
221 21204 : CALL section_vals_val_get(print_section, "DIPOLE%PERIODIC", l_val=do_berry)
222 21204 : use_charges = .FALSE.
223 21204 : IF (PRESENT(charges)) THEN
224 21204 : IF (ASSOCIATED(charges)) use_charges = .TRUE.
225 : END IF
226 :
227 21204 : CALL get_reference_point(rcc, drcc, fist_env=fist_env, reference=reference, ref_point=ref_point)
228 :
229 : ! Dipole deriv will be the derivative of the Dipole(dM/dt=\sum e_j v_j)
230 21204 : dipole_deriv = 0.0_dp
231 21204 : dipole = 0.0_dp
232 21204 : IF (do_berry) THEN
233 21204 : dipole_type = "periodic (Berry phase)"
234 84816 : rcc = pbc(rcc, cell)
235 21204 : charge_tot = 0._dp
236 21204 : IF (use_charges) THEN
237 2644 : charge_tot = SUM(charges)
238 : ELSE
239 2235036 : DO i = 1, SIZE(particle_set)
240 2214460 : atomic_kind => particle_set(i)%atomic_kind
241 2214460 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q)
242 2235036 : charge_tot = charge_tot + q
243 : END DO
244 : END IF
245 339264 : ria = twopi*MATMUL(cell%h_inv, rcc)
246 84816 : zphase = CMPLX(COS(ria), SIN(ria), dp)**charge_tot
247 :
248 339264 : dria = twopi*MATMUL(cell%h_inv, drcc)
249 84816 : dzphase = charge_tot*CMPLX(-SIN(ria), COS(ria), dp)**(charge_tot - 1.0_dp)*dria
250 :
251 84816 : ggamma = z_one
252 21204 : dggamma = z_zero
253 80920 : DO iparticle_kind = 1, SIZE(atomic_kind_set)
254 59716 : atomic_kind => atomic_kind_set(iparticle_kind)
255 59716 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q, atom_list=atom_list)
256 :
257 2297396 : DO i = 1, SIZE(atom_list)
258 8865904 : ria = particle_set(atom_list(i))%r(:)
259 8865904 : ria = pbc(ria, cell)
260 8865904 : via = particle_set(atom_list(i))%v(:)
261 2216476 : IF (use_charges) q = charges(atom_list(i))
262 8925620 : DO j = 1, 3
263 26597712 : gvec = twopi*cell%h_inv(j, :)
264 26597712 : theta = SUM(ria(:)*gvec(:))
265 26597712 : dtheta = SUM(via(:)*gvec(:))
266 6649428 : zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(-q)
267 6649428 : dzeta = -q*CMPLX(-SIN(theta), COS(theta), KIND=dp)**(-q - 1.0_dp)*dtheta
268 6649428 : dggamma(j) = dggamma(j)*zeta + ggamma(j)*dzeta
269 8865904 : ggamma(j) = ggamma(j)*zeta
270 : END DO
271 : END DO
272 : END DO
273 84816 : dggamma = dggamma*zphase + ggamma*dzphase
274 84816 : ggamma = ggamma*zphase
275 84816 : IF (ALL(REAL(ggamma, KIND=dp) /= 0.0_dp)) THEN
276 84816 : tmp = AIMAG(ggamma)/REAL(ggamma, KIND=dp)
277 84816 : ci = ATAN(tmp)
278 : dci = (1.0_dp/(1.0_dp + tmp**2))* &
279 84816 : (AIMAG(dggamma)*REAL(ggamma, KIND=dp) - AIMAG(ggamma)*REAL(dggamma, KIND=dp))/(REAL(ggamma, KIND=dp))**2
280 :
281 339264 : dipole = MATMUL(cell%hmat, ci)/twopi
282 339264 : dipole_deriv = MATMUL(cell%hmat, dci)/twopi
283 : END IF
284 21204 : CALL fist_env_get(fist_env=fist_env, results=results)
285 21204 : CALL cp_results_erase(results, description)
286 21204 : CALL put_results(results, description, dipole)
287 : ELSE
288 0 : dipole_type = "non-periodic"
289 0 : DO i = 1, SIZE(particle_set)
290 0 : atomic_kind => particle_set(i)%atomic_kind
291 0 : ria = particle_set(i)%r(:) ! no pbc(particle_set(i)%r(:),cell) so that the total dipole
292 : ! is the sum of the molecular dipoles
293 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q)
294 0 : IF (use_charges) q = charges(i)
295 0 : dipole = dipole - q*(ria - rcc)
296 0 : dipole_deriv(:) = dipole_deriv(:) - q*(particle_set(i)%v(:) - drcc)
297 : END DO
298 0 : CALL fist_env_get(fist_env=fist_env, results=results)
299 0 : CALL cp_results_erase(results, description)
300 0 : CALL put_results(results, description, dipole)
301 : END IF
302 21204 : IF (unit_nr > 0) THEN
303 : WRITE (unit_nr, '(/,T2,A,T31,A50)') &
304 10855 : 'MM_DIPOLE| Dipole type', ADJUSTR(TRIM(dipole_type))
305 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
306 10855 : 'MM_DIPOLE| Moment [a.u.]', dipole(1:3)
307 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
308 43420 : 'MM_DIPOLE| Moment [Debye]', dipole(1:3)*debye
309 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
310 10855 : 'MM_DIPOLE| Derivative [a.u.]', dipole_deriv(1:3)
311 : END IF
312 :
313 21204 : END SUBROUTINE fist_dipole
314 :
315 : END MODULE fist_efield_methods
|
