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 : !> \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 118 : 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 118 : 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
82 : TYPE(atomic_kind_type), POINTER :: atomic_kind
83 :
84 118 : qenergy = 0.0_dp
85 1526 : qforce = 0.0_dp
86 118 : qpv = 0.0_dp
87 :
88 118 : use_charges = .FALSE.
89 118 : IF (PRESENT(charges)) THEN
90 118 : IF (ASSOCIATED(charges)) use_charges = .TRUE.
91 : END IF
92 :
93 : IF (PRESENT(iunit)) THEN
94 118 : iw = iunit
95 : ELSE
96 118 : iw = -1
97 : END IF
98 :
99 118 : IF (PRESENT(use_virial)) THEN
100 118 : virial = use_virial
101 : ELSE
102 : virial = .FALSE.
103 : END IF
104 :
105 472 : fieldpol = efield%polarisation
106 826 : fieldpol = fieldpol/NORM2(fieldpol)
107 472 : fieldpol = -fieldpol*efield%strength
108 :
109 472 : dfilter = efield%dfilter
110 :
111 : dipole = 0.0_dp
112 472 : ggamma = z_one
113 354 : DO iparticle_kind = 1, SIZE(atomic_kind_set)
114 236 : atomic_kind => atomic_kind_set(iparticle_kind)
115 236 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q, atom_list=atom_list)
116 : ! TODO parallelization over atoms (local_particles)
117 706 : DO i = 1, SIZE(atom_list)
118 352 : ii = atom_list(i)
119 1408 : ria = particle_set(ii)%r(:)
120 1408 : ria = pbc(ria, cell)
121 352 : IF (use_charges) q = charges(ii)
122 1408 : DO j = 1, 3
123 4224 : gvec = twopi*cell%h_inv(j, :)
124 4224 : theta = SUM(ria(:)*gvec(:))
125 1056 : zeta = CMPLX(COS(q*theta), SIN(q*theta), KIND=dp)
126 1408 : ggamma(j) = ggamma(j)*zeta
127 : END DO
128 1644 : qforce(1:3, ii) = q
129 : END DO
130 : END DO
131 :
132 472 : ci = ATAN2(AIMAG(ggamma), REAL(ggamma, KIND=dp))
133 1888 : dipole = MATMUL(cell%hmat, ci)/twopi
134 :
135 118 : IF (efield%displacement) THEN
136 : ! E = (omega/8Pi)(D - 4Pi*P)^2
137 152 : di = dipole/cell%deth
138 152 : DO i = 1, 3
139 114 : theta = fieldpol(i) + 2._dp*twopi*di(i)
140 114 : qenergy = qenergy + dfilter(i)*theta**2
141 152 : fq(i) = -dfilter(i)*theta
142 : END DO
143 38 : qenergy = 0.25_dp*cell%deth/twopi*qenergy
144 152 : DO i = 1, SIZE(qforce, 2)
145 494 : qforce(1:3, i) = fq(1:3)*qforce(1:3, i)
146 : END DO
147 : ELSE
148 : ! E = -omega*E*P
149 320 : qenergy = SUM(fieldpol*dipole)
150 318 : DO i = 1, SIZE(qforce, 2)
151 1032 : qforce(1:3, i) = -fieldpol(1:3)*qforce(1:3, i)
152 : END DO
153 : END IF
154 :
155 118 : IF (virial) THEN
156 6 : DO iparticle_kind = 1, SIZE(atomic_kind_set)
157 4 : atomic_kind => atomic_kind_set(iparticle_kind)
158 4 : CALL get_atomic_kind(atomic_kind=atomic_kind, atom_list=atom_list)
159 12 : DO i = 1, SIZE(atom_list)
160 6 : ii = atom_list(i)
161 24 : ria = particle_set(ii)%r(:)
162 24 : ria = pbc(ria, cell)
163 28 : DO j = 1, 3
164 78 : qpv(j, 1:3) = qpv(j, 1:3) + qforce(j, ii)*ria(1:3)
165 : END DO
166 : END DO
167 : END DO
168 : ! Stress tensor for constant D needs further investigation
169 2 : IF (efield%displacement) THEN
170 0 : CPABORT("Stress Tensor for constant D simulation is not working")
171 : END IF
172 : END IF
173 :
174 118 : END SUBROUTINE fist_efield_energy_force
175 : ! **************************************************************************************************
176 : !> \brief Evaluates the Dipole of a classical charge distribution(point-like)
177 : !> possibly using the berry phase formalism
178 : !> \param fist_env ...
179 : !> \param print_section ...
180 : !> \param atomic_kind_set ...
181 : !> \param particle_set ...
182 : !> \param cell ...
183 : !> \param unit_nr ...
184 : !> \param charges ...
185 : !> \par History
186 : !> [01.2006] created
187 : !> [12.2007] tlaino - University of Zurich - debug and extended
188 : !> \author Teodoro Laino
189 : ! **************************************************************************************************
190 40972 : SUBROUTINE fist_dipole(fist_env, print_section, atomic_kind_set, particle_set, &
191 : cell, unit_nr, charges)
192 : TYPE(fist_environment_type), POINTER :: fist_env
193 : TYPE(section_vals_type), POINTER :: print_section
194 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
195 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
196 : TYPE(cell_type), POINTER :: cell
197 : INTEGER, INTENT(IN) :: unit_nr
198 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
199 :
200 : CHARACTER(LEN=default_string_length) :: description, dipole_type
201 : COMPLEX(KIND=dp) :: dzeta, dzphase(3), zeta, zphase(3)
202 : COMPLEX(KIND=dp), DIMENSION(3) :: dggamma, ggamma
203 : INTEGER :: i, iparticle_kind, j, reference
204 20486 : INTEGER, DIMENSION(:), POINTER :: atom_list
205 : LOGICAL :: do_berry, use_charges
206 : REAL(KIND=dp) :: charge_tot, ci(3), dci(3), dipole(3), &
207 : dipole_deriv(3), drcc(3), dria(3), &
208 : dtheta, gvec(3), q, rcc(3), ria(3), &
209 : theta, via(3)
210 20486 : REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
211 : TYPE(atomic_kind_type), POINTER :: atomic_kind
212 : TYPE(cp_result_type), POINTER :: results
213 :
214 20486 : NULLIFY (atomic_kind)
215 : ! Reference point
216 40972 : reference = section_get_ival(print_section, keyword_name="DIPOLE%REFERENCE")
217 20486 : NULLIFY (ref_point)
218 20486 : description = '[DIPOLE]'
219 20486 : CALL section_vals_val_get(print_section, "DIPOLE%REF_POINT", r_vals=ref_point)
220 20486 : CALL section_vals_val_get(print_section, "DIPOLE%PERIODIC", l_val=do_berry)
221 20486 : use_charges = .FALSE.
222 20486 : IF (PRESENT(charges)) THEN
223 20486 : IF (ASSOCIATED(charges)) use_charges = .TRUE.
224 : END IF
225 :
226 20486 : CALL get_reference_point(rcc, drcc, fist_env=fist_env, reference=reference, ref_point=ref_point)
227 :
228 : ! Dipole deriv will be the derivative of the Dipole(dM/dt=\sum e_j v_j)
229 20486 : dipole_deriv = 0.0_dp
230 20486 : dipole = 0.0_dp
231 20486 : IF (do_berry) THEN
232 20486 : dipole_type = "periodic (Berry phase)"
233 81944 : rcc = pbc(rcc, cell)
234 20486 : charge_tot = 0._dp
235 20486 : IF (use_charges) THEN
236 2644 : charge_tot = SUM(charges)
237 : ELSE
238 2202356 : DO i = 1, SIZE(particle_set)
239 2182498 : atomic_kind => particle_set(i)%atomic_kind
240 2182498 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q)
241 2202356 : charge_tot = charge_tot + q
242 : END DO
243 : END IF
244 327776 : ria = twopi*MATMUL(cell%h_inv, rcc)
245 81944 : zphase = CMPLX(COS(charge_tot*ria), -SIN(charge_tot*ria), KIND=dp)
246 :
247 327776 : dria = twopi*MATMUL(cell%h_inv, drcc)
248 81944 : dzphase = -charge_tot*CMPLX(SIN(charge_tot*ria), COS(charge_tot*ria), KIND=dp)*dria
249 :
250 81944 : ggamma = z_one
251 20486 : dggamma = z_zero
252 86664 : DO iparticle_kind = 1, SIZE(atomic_kind_set)
253 66178 : atomic_kind => atomic_kind_set(iparticle_kind)
254 66178 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q, atom_list=atom_list)
255 :
256 2271178 : DO i = 1, SIZE(atom_list)
257 8738056 : ria = particle_set(atom_list(i))%r(:)
258 8738056 : ria = pbc(ria, cell)
259 8738056 : via = particle_set(atom_list(i))%v(:)
260 2184514 : IF (use_charges) q = charges(atom_list(i))
261 8804234 : DO j = 1, 3
262 26214168 : gvec = twopi*cell%h_inv(j, :)
263 26214168 : theta = SUM(ria(:)*gvec(:))
264 26214168 : dtheta = SUM(via(:)*gvec(:))
265 6553542 : zeta = CMPLX(COS(q*theta), SIN(q*theta), KIND=dp)
266 6553542 : dzeta = q*CMPLX(-SIN(q*theta), COS(q*theta), KIND=dp)*dtheta
267 6553542 : dggamma(j) = dggamma(j)*zeta + ggamma(j)*dzeta
268 8738056 : ggamma(j) = ggamma(j)*zeta
269 : END DO
270 : END DO
271 : END DO
272 81944 : dggamma = dggamma*zphase + ggamma*dzphase
273 81944 : ggamma = ggamma*zphase
274 81944 : ci = ATAN2(AIMAG(ggamma), REAL(ggamma, KIND=dp))
275 : dci = (REAL(ggamma, KIND=dp)*AIMAG(dggamma) - &
276 81944 : AIMAG(ggamma)*REAL(dggamma, KIND=dp))/ABS(ggamma)**2
277 :
278 327776 : dipole = MATMUL(cell%hmat, ci)/twopi
279 327776 : dipole_deriv = MATMUL(cell%hmat, dci)/twopi
280 20486 : CALL fist_env_get(fist_env=fist_env, results=results)
281 20486 : CALL cp_results_erase(results, description)
282 20486 : CALL put_results(results, description, dipole)
283 : ELSE
284 0 : dipole_type = "non-periodic"
285 0 : DO i = 1, SIZE(particle_set)
286 0 : atomic_kind => particle_set(i)%atomic_kind
287 0 : ria = particle_set(i)%r(:) ! no pbc(particle_set(i)%r(:),cell) so that the total dipole
288 : ! is the sum of the molecular dipoles
289 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q)
290 0 : IF (use_charges) q = charges(i)
291 0 : dipole = dipole + q*(ria - rcc)
292 0 : dipole_deriv(:) = dipole_deriv(:) + q*(particle_set(i)%v(:) - drcc)
293 : END DO
294 0 : CALL fist_env_get(fist_env=fist_env, results=results)
295 0 : CALL cp_results_erase(results, description)
296 0 : CALL put_results(results, description, dipole)
297 : END IF
298 20486 : IF (unit_nr > 0) THEN
299 : WRITE (unit_nr, '(/,T2,A,T31,A50)') &
300 10496 : 'MM_DIPOLE| Dipole type', ADJUSTR(TRIM(dipole_type))
301 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
302 10496 : 'MM_DIPOLE| Moment [a.u.]', dipole(1:3)
303 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
304 41984 : 'MM_DIPOLE| Moment [Debye]', dipole(1:3)*debye
305 : WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
306 10496 : 'MM_DIPOLE| Derivative [a.u.]', dipole_deriv(1:3)
307 : END IF
308 :
309 20486 : END SUBROUTINE fist_dipole
310 :
311 : END MODULE fist_efield_methods
|
