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 : !> JGH [04042007] code refactoring
11 : ! **************************************************************************************************
12 : MODULE virial_methods
13 :
14 : USE atomic_kind_list_types, ONLY: atomic_kind_list_type
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind
17 : USE cell_types, ONLY: cell_type
18 : USE cp_subsys_types, ONLY: cp_subsys_get,&
19 : cp_subsys_type
20 : USE cp_units, ONLY: cp_unit_from_cp2k
21 : USE distribution_1d_types, ONLY: distribution_1d_type
22 : USE kinds, ONLY: default_string_length,&
23 : dp
24 : USE mathlib, ONLY: det_3x3,&
25 : jacobi
26 : USE message_passing, ONLY: mp_comm_type,&
27 : mp_para_env_type
28 : USE particle_list_types, ONLY: particle_list_type
29 : USE particle_types, ONLY: particle_type
30 : USE virial_types, ONLY: virial_type
31 : #include "./base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 :
35 : PRIVATE
36 : PUBLIC:: one_third_sum_diag, virial_evaluate, virial_pair_force, virial_update, &
37 : write_stress_tensor, write_stress_tensor_components
38 :
39 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'virial_methods'
40 :
41 : CONTAINS
42 : ! **************************************************************************************************
43 : !> \brief Updates the virial given the virial and subsys
44 : !> \param virial ...
45 : !> \param subsys ...
46 : !> \param para_env ...
47 : !> \par History
48 : !> none
49 : !> \author Teodoro Laino [tlaino] - 03.2008 - Zurich University
50 : ! **************************************************************************************************
51 6108 : SUBROUTINE virial_update(virial, subsys, para_env)
52 :
53 : TYPE(virial_type), INTENT(INOUT) :: virial
54 : TYPE(cp_subsys_type), POINTER :: subsys
55 : TYPE(mp_para_env_type), POINTER :: para_env
56 :
57 : TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
58 : TYPE(distribution_1d_type), POINTER :: local_particles
59 : TYPE(particle_list_type), POINTER :: particles
60 :
61 : CALL cp_subsys_get(subsys, local_particles=local_particles, atomic_kinds=atomic_kinds, &
62 6108 : particles=particles)
63 :
64 : CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles, &
65 6108 : virial, para_env)
66 :
67 6108 : END SUBROUTINE virial_update
68 :
69 : ! **************************************************************************************************
70 : !> \brief Computes the kinetic part of the pressure tensor and updates
71 : !> the full VIRIAL (PV)
72 : !> \param atomic_kind_set ...
73 : !> \param particle_set ...
74 : !> \param local_particles ...
75 : !> \param virial ...
76 : !> \param igroup ...
77 : !> \par History
78 : !> none
79 : !> \author CJM
80 : ! **************************************************************************************************
81 60665 : SUBROUTINE virial_evaluate(atomic_kind_set, particle_set, local_particles, &
82 : virial, igroup)
83 :
84 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
85 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
86 : TYPE(distribution_1d_type), POINTER :: local_particles
87 : TYPE(virial_type), INTENT(INOUT) :: virial
88 :
89 : CLASS(mp_comm_type), INTENT(IN) :: igroup
90 :
91 : CHARACTER(LEN=*), PARAMETER :: routineN = "virial_evaluate"
92 :
93 : INTEGER :: handle, i, iparticle, iparticle_kind, &
94 : iparticle_local, j, nparticle_kind, &
95 : nparticle_local
96 : REAL(KIND=dp) :: mass
97 : TYPE(atomic_kind_type), POINTER :: atomic_kind
98 :
99 60665 : IF (virial%pv_availability) THEN
100 19626 : CALL timeset(routineN, handle)
101 19626 : NULLIFY (atomic_kind)
102 19626 : nparticle_kind = SIZE(atomic_kind_set)
103 255138 : virial%pv_kinetic = 0.0_dp
104 78504 : DO i = 1, 3
105 196260 : DO j = 1, i
106 336816 : DO iparticle_kind = 1, nparticle_kind
107 219060 : atomic_kind => atomic_kind_set(iparticle_kind)
108 219060 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
109 219060 : nparticle_local = local_particles%n_el(iparticle_kind)
110 6481092 : DO iparticle_local = 1, nparticle_local
111 6144276 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
112 : virial%pv_kinetic(i, j) = virial%pv_kinetic(i, j) + &
113 6363336 : mass*particle_set(iparticle)%v(i)*particle_set(iparticle)%v(j)
114 : END DO
115 : END DO
116 176634 : virial%pv_kinetic(j, i) = virial%pv_kinetic(i, j)
117 : END DO
118 : END DO
119 :
120 19626 : CALL igroup%sum(virial%pv_kinetic)
121 :
122 : ! total virial
123 255138 : virial%pv_total = virial%pv_virial + virial%pv_kinetic + virial%pv_constraint
124 :
125 19626 : CALL timestop(handle)
126 : END IF
127 :
128 60665 : END SUBROUTINE virial_evaluate
129 :
130 : ! **************************************************************************************************
131 : !> \brief Computes the contribution to the stress tensor from two-body
132 : !> pair-wise forces
133 : !> \param pv_virial ...
134 : !> \param f0 ...
135 : !> \param force ...
136 : !> \param rab ...
137 : !> \par History
138 : !> none
139 : !> \author JGH
140 : ! **************************************************************************************************
141 76953943 : PURE SUBROUTINE virial_pair_force(pv_virial, f0, force, rab)
142 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: pv_virial
143 : REAL(KIND=dp), INTENT(IN) :: f0
144 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: force, rab
145 :
146 : INTEGER :: i, j
147 :
148 307815772 : DO i = 1, 3
149 1000401259 : DO j = 1, 3
150 923447316 : pv_virial(i, j) = pv_virial(i, j) + f0*force(i)*rab(j)
151 : END DO
152 : END DO
153 :
154 76953943 : END SUBROUTINE virial_pair_force
155 :
156 : ! **************************************************************************************************
157 : !> \brief ...
158 : !> \param virial ...
159 : !> \param iw ...
160 : !> \param cell ...
161 : !> \param unit_string ...
162 : !> \par History
163 : !> - Revised virial components (14.10.2020, MK)
164 : !> \author JGH
165 : ! **************************************************************************************************
166 113 : SUBROUTINE write_stress_tensor_components(virial, iw, cell, unit_string)
167 :
168 : TYPE(virial_type), INTENT(IN) :: virial
169 : INTEGER, INTENT(IN) :: iw
170 : TYPE(cell_type), POINTER :: cell
171 : CHARACTER(LEN=default_string_length), INTENT(IN) :: unit_string
172 :
173 : CHARACTER(LEN=*), PARAMETER :: fmt = "(T2,A,T41,2(1X,ES19.11))"
174 :
175 : REAL(KIND=dp) :: fconv
176 : REAL(KIND=dp), DIMENSION(3, 3) :: pv
177 :
178 113 : IF (iw > 0) THEN
179 113 : CPASSERT(ASSOCIATED(cell))
180 : ! Conversion factor
181 113 : fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, TRIM(ADJUSTL(unit_string)))
182 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
183 113 : "STRESS| Stress tensor components (GPW/GAPW) ["//TRIM(ADJUSTL(unit_string))//"]"
184 : WRITE (UNIT=iw, FMT="(T2,A,T52,A,T70,A)") &
185 113 : "STRESS|", "1/3 Trace", "Determinant"
186 1469 : pv = fconv*virial%pv_overlap
187 : WRITE (UNIT=iw, FMT=fmt) &
188 113 : "STRESS| Overlap", one_third_sum_diag(pv), det_3x3(pv)
189 1469 : pv = fconv*virial%pv_ekinetic
190 : WRITE (UNIT=iw, FMT=fmt) &
191 113 : "STRESS| Kinetic energy", one_third_sum_diag(pv), det_3x3(pv)
192 1469 : pv = fconv*virial%pv_ppl
193 : WRITE (UNIT=iw, FMT=fmt) &
194 113 : "STRESS| Local pseudopotential/core", one_third_sum_diag(pv), det_3x3(pv)
195 1469 : pv = fconv*virial%pv_ppnl
196 : WRITE (UNIT=iw, FMT=fmt) &
197 113 : "STRESS| Nonlocal pseudopotential", one_third_sum_diag(pv), det_3x3(pv)
198 1469 : pv = fconv*virial%pv_ecore_overlap
199 : WRITE (UNIT=iw, FMT=fmt) &
200 113 : "STRESS| Core charge overlap", one_third_sum_diag(pv), det_3x3(pv)
201 1469 : pv = fconv*virial%pv_ehartree
202 : WRITE (UNIT=iw, FMT=fmt) &
203 113 : "STRESS| Hartree", one_third_sum_diag(pv), det_3x3(pv)
204 1469 : pv = fconv*virial%pv_exc
205 : WRITE (UNIT=iw, FMT=fmt) &
206 113 : "STRESS| Exchange-correlation", one_third_sum_diag(pv), det_3x3(pv)
207 1469 : pv = fconv*virial%pv_exx
208 : WRITE (UNIT=iw, FMT=fmt) &
209 113 : "STRESS| Exact exchange (EXX)", one_third_sum_diag(pv), det_3x3(pv)
210 1469 : pv = fconv*virial%pv_vdw
211 : WRITE (UNIT=iw, FMT=fmt) &
212 113 : "STRESS| vdW correction", one_third_sum_diag(pv), det_3x3(pv)
213 1469 : pv = fconv*virial%pv_mp2
214 : WRITE (UNIT=iw, FMT=fmt) &
215 113 : "STRESS| Moller-Plesset (MP2)", one_third_sum_diag(pv), det_3x3(pv)
216 1469 : pv = fconv*virial%pv_nlcc
217 : WRITE (UNIT=iw, FMT=fmt) &
218 113 : "STRESS| Nonlinear core correction (NLCC)", one_third_sum_diag(pv), det_3x3(pv)
219 1469 : pv = fconv*virial%pv_gapw
220 : WRITE (UNIT=iw, FMT=fmt) &
221 113 : "STRESS| Local atomic parts (GAPW)", one_third_sum_diag(pv), det_3x3(pv)
222 1469 : pv = fconv*virial%pv_lrigpw
223 : WRITE (UNIT=iw, FMT=fmt) &
224 113 : "STRESS| Resolution-of-the-identity (LRI)", one_third_sum_diag(pv), det_3x3(pv)
225 : pv = fconv*(virial%pv_overlap + virial%pv_ekinetic + virial%pv_ppl + virial%pv_ppnl + &
226 : virial%pv_ecore_overlap + virial%pv_ehartree + virial%pv_exc + &
227 : virial%pv_exx + virial%pv_vdw + virial%pv_mp2 + virial%pv_nlcc + &
228 1469 : virial%pv_gapw + virial%pv_lrigpw)
229 : WRITE (UNIT=iw, FMT=fmt) &
230 113 : "STRESS| Sum of components", one_third_sum_diag(pv), det_3x3(pv)
231 1469 : pv = fconv*virial%pv_virial
232 : WRITE (UNIT=iw, FMT=fmt) &
233 113 : "STRESS| Total", one_third_sum_diag(pv), det_3x3(pv)
234 : END IF
235 :
236 113 : END SUBROUTINE write_stress_tensor_components
237 :
238 : ! **************************************************************************************************
239 : !> \brief ...
240 : !> \param a ...
241 : !> \return ...
242 : ! **************************************************************************************************
243 1770 : PURE FUNCTION one_third_sum_diag(a) RESULT(p)
244 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
245 : REAL(KIND=dp) :: p
246 :
247 1770 : p = (a(1, 1) + a(2, 2) + a(3, 3))/3.0_dp
248 1770 : END FUNCTION one_third_sum_diag
249 :
250 : ! **************************************************************************************************
251 : !> \brief Print stress tensor to output file
252 : !>
253 : !> \param pv_virial ...
254 : !> \param iw ...
255 : !> \param cell ...
256 : !> \param unit_string ...
257 : !> \param numerical ...
258 : !> \author MK (26.08.2010)
259 : ! **************************************************************************************************
260 6145 : SUBROUTINE write_stress_tensor(pv_virial, iw, cell, unit_string, numerical)
261 :
262 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: pv_virial
263 : INTEGER, INTENT(IN) :: iw
264 : TYPE(cell_type), POINTER :: cell
265 : CHARACTER(LEN=default_string_length), INTENT(IN) :: unit_string
266 : LOGICAL, INTENT(IN) :: numerical
267 :
268 : REAL(KIND=dp) :: fconv
269 : REAL(KIND=dp), DIMENSION(3) :: eigval
270 : REAL(KIND=dp), DIMENSION(3, 3) :: eigvec, stress_tensor
271 :
272 6145 : IF (iw > 0) THEN
273 6145 : CPASSERT(ASSOCIATED(cell))
274 : ! Conversion factor
275 6145 : fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, TRIM(ADJUSTL(unit_string)))
276 79885 : stress_tensor(:, :) = fconv*pv_virial(:, :)
277 6145 : IF (numerical) THEN
278 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
279 13 : "STRESS| Numerical stress tensor ["//TRIM(ADJUSTL(unit_string))//"]"
280 : ELSE
281 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
282 6132 : "STRESS| Analytical stress tensor ["//TRIM(ADJUSTL(unit_string))//"]"
283 : END IF
284 : WRITE (UNIT=iw, FMT="(T2,A,T14,3(19X,A1))") &
285 6145 : "STRESS|", "x", "y", "z"
286 : WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
287 6145 : "STRESS| x", stress_tensor(1, 1:3)
288 : WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
289 6145 : "STRESS| y", stress_tensor(2, 1:3)
290 : WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
291 6145 : "STRESS| z", stress_tensor(3, 1:3)
292 : WRITE (UNIT=iw, FMT="(T2,A,T61,ES20.11)") &
293 6145 : "STRESS| 1/3 Trace", (stress_tensor(1, 1) + &
294 : stress_tensor(2, 2) + &
295 12290 : stress_tensor(3, 3))/3.0_dp
296 : WRITE (UNIT=iw, FMT="(T2,A,T61,ES20.11)") &
297 6145 : "STRESS| Determinant", det_3x3(stress_tensor(1:3, 1), &
298 : stress_tensor(1:3, 2), &
299 12290 : stress_tensor(1:3, 3))
300 6145 : eigval(:) = 0.0_dp
301 6145 : eigvec(:, :) = 0.0_dp
302 6145 : CALL jacobi(stress_tensor, eigval, eigvec)
303 6145 : IF (numerical) THEN
304 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
305 : "STRESS| Eigenvectors and eigenvalues of the numerical stress tensor ["// &
306 13 : TRIM(ADJUSTL(unit_string))//"]"
307 : ELSE
308 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
309 : "STRESS| Eigenvectors and eigenvalues of the analytical stress tensor ["// &
310 6132 : TRIM(ADJUSTL(unit_string))//"]"
311 : END IF
312 : WRITE (UNIT=iw, FMT="(T2,A,T14,3(1X,I19))") &
313 6145 : "STRESS|", 1, 2, 3
314 : WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,ES19.11))") &
315 6145 : "STRESS| Eigenvalues", eigval(1:3)
316 : WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,F19.12))") &
317 6145 : "STRESS| x", eigvec(1, 1:3)
318 : WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,F19.12))") &
319 6145 : "STRESS| y", eigvec(2, 1:3)
320 : WRITE (UNIT=iw, FMT="(T2,A,T21,3(1X,F19.12))") &
321 6145 : "STRESS| z", eigvec(3, 1:3)
322 : END IF
323 :
324 6145 : END SUBROUTINE write_stress_tensor
325 :
326 : END MODULE virial_methods
|
