(git:ccc2433)
virial_methods.F
Go to the documentation of this file.
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  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  particles=particles)
61 
62  CALL virial_evaluate(atomic_kinds%els, particles%els, local_particles, &
63  virial, para_env)
64 
65  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  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  IF (virial%pv_availability) THEN
97  CALL timeset(routinen, handle)
98  NULLIFY (atomic_kind)
99  nparticle_kind = SIZE(atomic_kind_set)
100  virial%pv_kinetic = 0.0_dp
101  DO i = 1, 3
102  DO j = 1, i
103  DO iparticle_kind = 1, nparticle_kind
104  atomic_kind => atomic_kind_set(iparticle_kind)
105  CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
106  nparticle_local = local_particles%n_el(iparticle_kind)
107  DO iparticle_local = 1, nparticle_local
108  iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
109  virial%pv_kinetic(i, j) = virial%pv_kinetic(i, j) + &
110  mass*particle_set(iparticle)%v(i)*particle_set(iparticle)%v(j)
111  END DO
112  END DO
113  virial%pv_kinetic(j, i) = virial%pv_kinetic(i, j)
114  END DO
115  END DO
116 
117  CALL igroup%sum(virial%pv_kinetic)
118 
119  ! total virial
120  virial%pv_total = virial%pv_virial + virial%pv_kinetic + virial%pv_constraint
121 
122  CALL timestop(handle)
123  END IF
124 
125  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  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  DO i = 1, 3
146  DO j = 1, 3
147  pv_virial(i, j) = pv_virial(i, j) + f0*force(i)*rab(j)
148  END DO
149  END DO
150 
151  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  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  IF (iw > 0) THEN
173  cpassert(ASSOCIATED(cell))
174  ! Conversion factor to GPa
175  fconv = 1.0e-9_dp*pascal/cell%deth
176  WRITE (unit=iw, fmt='(/,T2,A)') &
177  'STRESS| Stress tensor components (GPW/GAPW) [GPa]'
178  WRITE (unit=iw, fmt='(T2,A,T52,A,T70,A)') &
179  'STRESS|', '1/3 Trace', 'Determinant'
180  pv = fconv*virial%pv_overlap
181  WRITE (unit=iw, fmt=fmt) &
182  'STRESS| Overlap', one_third_sum_diag(pv), det_3x3(pv)
183  pv = fconv*virial%pv_ekinetic
184  WRITE (unit=iw, fmt=fmt) &
185  'STRESS| Kinetic energy', one_third_sum_diag(pv), det_3x3(pv)
186  pv = fconv*virial%pv_ppl
187  WRITE (unit=iw, fmt=fmt) &
188  'STRESS| Local pseudopotential/core', one_third_sum_diag(pv), det_3x3(pv)
189  pv = fconv*virial%pv_ppnl
190  WRITE (unit=iw, fmt=fmt) &
191  'STRESS| Nonlocal pseudopotential', one_third_sum_diag(pv), det_3x3(pv)
192  pv = fconv*virial%pv_ecore_overlap
193  WRITE (unit=iw, fmt=fmt) &
194  'STRESS| Core charge overlap', one_third_sum_diag(pv), det_3x3(pv)
195  pv = fconv*virial%pv_ehartree
196  WRITE (unit=iw, fmt=fmt) &
197  'STRESS| Hartree', one_third_sum_diag(pv), det_3x3(pv)
198  pv = fconv*virial%pv_exc
199  WRITE (unit=iw, fmt=fmt) &
200  'STRESS| Exchange-correlation', one_third_sum_diag(pv), det_3x3(pv)
201  pv = fconv*virial%pv_exx
202  WRITE (unit=iw, fmt=fmt) &
203  'STRESS| Exact exchange (EXX)', one_third_sum_diag(pv), det_3x3(pv)
204  pv = fconv*virial%pv_vdw
205  WRITE (unit=iw, fmt=fmt) &
206  'STRESS| vdW correction', one_third_sum_diag(pv), det_3x3(pv)
207  pv = fconv*virial%pv_mp2
208  WRITE (unit=iw, fmt=fmt) &
209  'STRESS| Moller-Plesset (MP2)', one_third_sum_diag(pv), det_3x3(pv)
210  pv = fconv*virial%pv_nlcc
211  WRITE (unit=iw, fmt=fmt) &
212  'STRESS| Nonlinear core correction (NLCC)', one_third_sum_diag(pv), det_3x3(pv)
213  pv = fconv*virial%pv_gapw
214  WRITE (unit=iw, fmt=fmt) &
215  'STRESS| Local atomic parts (GAPW)', one_third_sum_diag(pv), det_3x3(pv)
216  pv = fconv*virial%pv_lrigpw
217  WRITE (unit=iw, fmt=fmt) &
218  '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  virial%pv_gapw + virial%pv_lrigpw)
223  WRITE (unit=iw, fmt=fmt) &
224  'STRESS| Sum of components', one_third_sum_diag(pv), det_3x3(pv)
225  pv = fconv*virial%pv_virial
226  WRITE (unit=iw, fmt=fmt) &
227  'STRESS| Total', one_third_sum_diag(pv), det_3x3(pv)
228  END IF
229 
230  END SUBROUTINE write_stress_tensor_components
231 
232 ! **************************************************************************************************
233 !> \brief ...
234 !> \param a ...
235 !> \return ...
236 ! **************************************************************************************************
237  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  p = (a(1, 1) + a(2, 2) + a(3, 3))/3.0_dp
242  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  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  IF (iw > 0) THEN
264  cpassert(ASSOCIATED(cell))
265  stress_tensor(:, :) = pv_virial(:, :)/cell%deth*pascal*1.0e-9_dp
266  IF (numerical) THEN
267  WRITE (unit=iw, fmt='(/,T2,A)') &
268  'STRESS| Numerical stress tensor [GPa]'
269  ELSE
270  WRITE (unit=iw, fmt='(/,T2,A)') &
271  'STRESS| Analytical stress tensor [GPa]'
272  END IF
273  WRITE (unit=iw, fmt='(T2,A,T14,3(19X,A1))') &
274  'STRESS|', 'x', 'y', 'z'
275  WRITE (unit=iw, fmt='(T2,A,T21,3(1X,ES19.11))') &
276  'STRESS| x', stress_tensor(1, 1:3)
277  WRITE (unit=iw, fmt='(T2,A,T21,3(1X,ES19.11))') &
278  'STRESS| y', stress_tensor(2, 1:3)
279  WRITE (unit=iw, fmt='(T2,A,T21,3(1X,ES19.11))') &
280  'STRESS| z', stress_tensor(3, 1:3)
281  WRITE (unit=iw, fmt='(T2,A,T61,ES20.11)') &
282  'STRESS| 1/3 Trace', (stress_tensor(1, 1) + &
283  stress_tensor(2, 2) + &
284  stress_tensor(3, 3))/3.0_dp
285  WRITE (unit=iw, fmt='(T2,A,T61,ES20.11)') &
286  'STRESS| Determinant', det_3x3(stress_tensor(1:3, 1), &
287  stress_tensor(1:3, 2), &
288  stress_tensor(1:3, 3))
289  eigval(:) = 0.0_dp
290  eigvec(:, :) = 0.0_dp
291  CALL jacobi(stress_tensor, eigval, eigvec)
292  IF (numerical) THEN
293  WRITE (unit=iw, fmt='(/,T2,A)') &
294  'STRESS| Eigenvectors and eigenvalues of the numerical stress tensor [GPa]'
295  ELSE
296  WRITE (unit=iw, fmt='(/,T2,A)') &
297  '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  'STRESS|', 1, 2, 3
301  WRITE (unit=iw, fmt='(T2,A,T21,3(1X,ES19.11))') &
302  'STRESS| Eigenvalues', eigval(1:3)
303  WRITE (unit=iw, fmt='(T2,A,T21,3(1X,F19.12))') &
304  'STRESS| x', eigvec(1, 1:3)
305  WRITE (unit=iw, fmt='(T2,A,T21,3(1X,F19.12))') &
306  'STRESS| y', eigvec(2, 1:3)
307  WRITE (unit=iw, fmt='(T2,A,T21,3(1X,F19.12))') &
308  'STRESS| z', eigvec(3, 1:3)
309  END IF
310 
311  END SUBROUTINE write_stress_tensor
312 
313 END MODULE virial_methods
314 
det_3x3
real(kind=dp) function det_3x3(a)
...
Definition: dumpdcd.F:1091
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
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
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:1479
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
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::pascal
real(kind=dp), parameter, public pascal
Definition: physcon.F:174
virial_methods
Definition: virial_methods.F:12
virial_methods::write_stress_tensor_components
subroutine, public write_stress_tensor_components(virial, iw, cell)
...
Definition: virial_methods.F:163
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:139
virial_methods::write_stress_tensor
subroutine, public write_stress_tensor(pv_virial, iw, cell, numerical)
Print stress tensor to output file.
Definition: virial_methods.F:254
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:81
virial_methods::one_third_sum_diag
pure real(kind=dp) function, public one_third_sum_diag(a)
...
Definition: virial_methods.F:238
virial_methods::virial_update
subroutine, public virial_update(virial, subsys, para_env)
Updates the virial given the virial and subsys.
Definition: virial_methods.F:51
virial_types
Definition: virial_types.F:13