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