(git:b77b4be)
Loading...
Searching...
No Matches
virial_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
41CONTAINS
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 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 particles=particles)
63
64 CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles, &
65 virial, para_env)
66
67 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 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 IF (virial%pv_availability) THEN
100 CALL timeset(routinen, handle)
101 NULLIFY (atomic_kind)
102 nparticle_kind = SIZE(atomic_kind_set)
103 virial%pv_kinetic = 0.0_dp
104 DO i = 1, 3
105 DO j = 1, i
106 DO iparticle_kind = 1, nparticle_kind
107 atomic_kind => atomic_kind_set(iparticle_kind)
108 CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
109 nparticle_local = local_particles%n_el(iparticle_kind)
110 DO iparticle_local = 1, nparticle_local
111 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
112 virial%pv_kinetic(i, j) = virial%pv_kinetic(i, j) + &
113 mass*particle_set(iparticle)%v(i)*particle_set(iparticle)%v(j)
114 END DO
115 END DO
116 virial%pv_kinetic(j, i) = virial%pv_kinetic(i, j)
117 END DO
118 END DO
119
120 CALL igroup%sum(virial%pv_kinetic)
121
122 ! total virial
123 virial%pv_total = virial%pv_virial + virial%pv_kinetic + virial%pv_constraint
124
125 CALL timestop(handle)
126 END IF
127
128 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 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 DO i = 1, 3
149 DO j = 1, 3
150 pv_virial(i, j) = pv_virial(i, j) + f0*force(i)*rab(j)
151 END DO
152 END DO
153
154 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 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 IF (iw > 0) THEN
179 cpassert(ASSOCIATED(cell))
180 ! Conversion factor
181 fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, trim(adjustl(unit_string)))
182 WRITE (unit=iw, fmt="(/,T2,A)") &
183 "STRESS| Stress tensor components (GPW/GAPW) ["//trim(adjustl(unit_string))//"]"
184 WRITE (unit=iw, fmt="(T2,A,T52,A,T70,A)") &
185 "STRESS|", "1/3 Trace", "Determinant"
186 pv = fconv*virial%pv_overlap
187 WRITE (unit=iw, fmt=fmt) &
188 "STRESS| Overlap", one_third_sum_diag(pv), det_3x3(pv)
189 pv = fconv*virial%pv_ekinetic
190 WRITE (unit=iw, fmt=fmt) &
191 "STRESS| Kinetic energy", one_third_sum_diag(pv), det_3x3(pv)
192 pv = fconv*virial%pv_ppl
193 WRITE (unit=iw, fmt=fmt) &
194 "STRESS| Local pseudopotential/core", one_third_sum_diag(pv), det_3x3(pv)
195 pv = fconv*virial%pv_ppnl
196 WRITE (unit=iw, fmt=fmt) &
197 "STRESS| Nonlocal pseudopotential", one_third_sum_diag(pv), det_3x3(pv)
198 pv = fconv*virial%pv_ecore_overlap
199 WRITE (unit=iw, fmt=fmt) &
200 "STRESS| Core charge overlap", one_third_sum_diag(pv), det_3x3(pv)
201 pv = fconv*virial%pv_ehartree
202 WRITE (unit=iw, fmt=fmt) &
203 "STRESS| Hartree", one_third_sum_diag(pv), det_3x3(pv)
204 pv = fconv*virial%pv_exc
205 WRITE (unit=iw, fmt=fmt) &
206 "STRESS| Exchange-correlation", one_third_sum_diag(pv), det_3x3(pv)
207 pv = fconv*virial%pv_exx
208 WRITE (unit=iw, fmt=fmt) &
209 "STRESS| Exact exchange (EXX)", one_third_sum_diag(pv), det_3x3(pv)
210 pv = fconv*virial%pv_vdw
211 WRITE (unit=iw, fmt=fmt) &
212 "STRESS| vdW correction", one_third_sum_diag(pv), det_3x3(pv)
213 pv = fconv*virial%pv_mp2
214 WRITE (unit=iw, fmt=fmt) &
215 "STRESS| Moller-Plesset (MP2)", one_third_sum_diag(pv), det_3x3(pv)
216 pv = fconv*virial%pv_nlcc
217 WRITE (unit=iw, fmt=fmt) &
218 "STRESS| Nonlinear core correction (NLCC)", one_third_sum_diag(pv), det_3x3(pv)
219 pv = fconv*virial%pv_gapw
220 WRITE (unit=iw, fmt=fmt) &
221 "STRESS| Local atomic parts (GAPW)", one_third_sum_diag(pv), det_3x3(pv)
222 pv = fconv*virial%pv_lrigpw
223 WRITE (unit=iw, fmt=fmt) &
224 "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 virial%pv_gapw + virial%pv_lrigpw)
229 WRITE (unit=iw, fmt=fmt) &
230 "STRESS| Sum of components", one_third_sum_diag(pv), det_3x3(pv)
231 pv = fconv*virial%pv_virial
232 WRITE (unit=iw, fmt=fmt) &
233 "STRESS| Total", one_third_sum_diag(pv), det_3x3(pv)
234 END IF
235
236 END SUBROUTINE write_stress_tensor_components
237
238! **************************************************************************************************
239!> \brief ...
240!> \param a ...
241!> \return ...
242! **************************************************************************************************
243 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 p = (a(1, 1) + a(2, 2) + a(3, 3))/3.0_dp
248 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 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 IF (iw > 0) THEN
273 cpassert(ASSOCIATED(cell))
274 ! Conversion factor
275 fconv = cp_unit_from_cp2k(1.0_dp/cell%deth, trim(adjustl(unit_string)))
276 stress_tensor(:, :) = fconv*pv_virial(:, :)
277 IF (numerical) THEN
278 WRITE (unit=iw, fmt="(/,T2,A)") &
279 "STRESS| Numerical stress tensor ["//trim(adjustl(unit_string))//"]"
280 ELSE
281 WRITE (unit=iw, fmt="(/,T2,A)") &
282 "STRESS| Analytical stress tensor ["//trim(adjustl(unit_string))//"]"
283 END IF
284 WRITE (unit=iw, fmt="(T2,A,T14,3(19X,A1))") &
285 "STRESS|", "x", "y", "z"
286 WRITE (unit=iw, fmt="(T2,A,T21,3(1X,ES19.11))") &
287 "STRESS| x", stress_tensor(1, 1:3)
288 WRITE (unit=iw, fmt="(T2,A,T21,3(1X,ES19.11))") &
289 "STRESS| y", stress_tensor(2, 1:3)
290 WRITE (unit=iw, fmt="(T2,A,T21,3(1X,ES19.11))") &
291 "STRESS| z", stress_tensor(3, 1:3)
292 WRITE (unit=iw, fmt="(T2,A,T61,ES20.11)") &
293 "STRESS| 1/3 Trace", (stress_tensor(1, 1) + &
294 stress_tensor(2, 2) + &
295 stress_tensor(3, 3))/3.0_dp
296 WRITE (unit=iw, fmt="(T2,A,T61,ES20.11)") &
297 "STRESS| Determinant", det_3x3(stress_tensor(1:3, 1), &
298 stress_tensor(1:3, 2), &
299 stress_tensor(1:3, 3))
300 eigval(:) = 0.0_dp
301 eigvec(:, :) = 0.0_dp
302 CALL jacobi(stress_tensor, eigval, eigvec)
303 IF (numerical) THEN
304 WRITE (unit=iw, fmt="(/,T2,A)") &
305 "STRESS| Eigenvectors and eigenvalues of the numerical stress tensor ["// &
306 trim(adjustl(unit_string))//"]"
307 ELSE
308 WRITE (unit=iw, fmt="(/,T2,A)") &
309 "STRESS| Eigenvectors and eigenvalues of the analytical stress tensor ["// &
310 trim(adjustl(unit_string))//"]"
311 END IF
312 WRITE (unit=iw, fmt="(T2,A,T14,3(1X,I19))") &
313 "STRESS|", 1, 2, 3
314 WRITE (unit=iw, fmt="(T2,A,T21,3(1X,ES19.11))") &
315 "STRESS| Eigenvalues", eigval(1:3)
316 WRITE (unit=iw, fmt="(T2,A,T21,3(1X,F19.12))") &
317 "STRESS| x", eigvec(1, 1:3)
318 WRITE (unit=iw, fmt="(T2,A,T21,3(1X,F19.12))") &
319 "STRESS| y", eigvec(2, 1:3)
320 WRITE (unit=iw, fmt="(T2,A,T21,3(1X,F19.12))") &
321 "STRESS| z", eigvec(3, 1:3)
322 END IF
323
324 END SUBROUTINE write_stress_tensor
325
326END MODULE virial_methods
mathlib::det_3x3
Definition mathlib.F:66
atomic_kind_list_types
represent a simple array based list of the given type
Definition atomic_kind_list_types.F:15
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:345
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1179
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::jacobi
subroutine, public jacobi(a, d, v)
Jacobi matrix diagonalization. The eigenvalues are returned in vector d and the eigenvectors are retu...
Definition mathlib.F:1468
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F:142
virial_methods::write_stress_tensor
subroutine, public write_stress_tensor(pv_virial, iw, cell, unit_string, numerical)
Print stress tensor to output file.
Definition virial_methods.F:261
virial_methods::virial_evaluate
subroutine, public virial_evaluate(atomic_kind_set, particle_set, local_particles, virial, igroup)
Computes the kinetic part of the pressure tensor and updates the full VIRIAL (PV)
Definition virial_methods.F:83
virial_methods::write_stress_tensor_components
subroutine, public write_stress_tensor_components(virial, iw, cell, unit_string)
...
Definition virial_methods.F:167
virial_methods::one_third_sum_diag
pure real(kind=dp) function, public one_third_sum_diag(a)
...
Definition virial_methods.F:244
virial_methods::virial_update
subroutine, public virial_update(virial, subsys, para_env)
Updates the virial given the virial and subsys.
Definition virial_methods.F:52
virial_types
Definition virial_types.F:13
atomic_kind_list_types::atomic_kind_list_type
represent a list of objects
Definition atomic_kind_list_types.F:45
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
virial_types::virial_type
Definition virial_types.F:26