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