Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Space Group Symmetry Module (version 1.0, January 16, 2020)
10 : !> \par History
11 : !> Pierre-André Cazade [pcazade] 01.2020 - University of Limerick
12 : !> \author Pierre-André Cazade (first version)
13 : ! **************************************************************************************************
14 : MODULE space_groups
15 : USE atomic_kind_types, ONLY: get_atomic_kind
16 : USE bibliography, ONLY: Togo2018,&
17 : cite_reference
18 : USE cell_methods, ONLY: cell_create,&
19 : init_cell
20 : USE cell_types, ONLY: cell_copy,&
21 : cell_type,&
22 : real_to_scaled,&
23 : scaled_to_real
24 : USE cp_subsys_types, ONLY: cp_subsys_get,&
25 : cp_subsys_type
26 : USE gopt_f_types, ONLY: gopt_f_type
27 : USE input_constants, ONLY: default_cell_direct_id,&
28 : default_cell_geo_opt_id,&
29 : default_cell_md_id,&
30 : default_cell_method_id,&
31 : default_minimization_method_id,&
32 : default_ts_method_id
33 : USE input_section_types, ONLY: section_vals_type,&
34 : section_vals_val_get
35 : USE kinds, ONLY: dp
36 : USE mathlib, ONLY: det_3x3,&
37 : inv_3x3,&
38 : jacobi
39 : USE particle_list_types, ONLY: particle_list_type
40 : USE physcon, ONLY: pascal
41 : USE space_groups_types, ONLY: cleanup_spgr_type,&
42 : spgr_type
43 : USE spglib_f08, ONLY: spg_get_international,&
44 : spg_get_multiplicity,&
45 : spg_get_pointgroup,&
46 : spg_get_schoenflies,&
47 : spg_get_symmetry
48 : USE string_utilities, ONLY: strlcpy_c2f
49 : #include "../base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 :
53 : PRIVATE
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'space_groups'
56 :
57 : PUBLIC :: spgr_create, identify_space_group, spgr_find_equivalent_atoms
58 : PUBLIC :: spgr_apply_rotations_coord, spgr_apply_rotations_force, print_spgr
59 : PUBLIC :: spgr_apply_rotations_stress, spgr_write_stress_tensor
60 :
61 : CONTAINS
62 :
63 : ! **************************************************************************************************
64 : !> \brief routine creates the space group structure
65 : !> \param scoor ...
66 : !> \param types ...
67 : !> \param cell ...
68 : !> \param gopt_env ...
69 : !> \param eps_symmetry ...
70 : !> \param pol ...
71 : !> \param ranges ...
72 : !> \param nparticle ...
73 : !> \param n_atom ...
74 : !> \param n_core ...
75 : !> \param n_shell ...
76 : !> \param iunit ...
77 : !> \param print_atoms ...
78 : !> \par History
79 : !> 01.2020 created [pcazade]
80 : !> \author Pierre-André Cazade (first version)
81 : ! **************************************************************************************************
82 14 : SUBROUTINE spgr_create(scoor, types, cell, gopt_env, eps_symmetry, pol, ranges, &
83 : nparticle, n_atom, n_core, n_shell, iunit, print_atoms)
84 :
85 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: scoor
86 : INTEGER, DIMENSION(:), INTENT(IN) :: types
87 : TYPE(cell_type), INTENT(IN), POINTER :: cell
88 : TYPE(gopt_f_type), INTENT(IN), POINTER :: gopt_env
89 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: eps_symmetry
90 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: pol
91 : INTEGER, DIMENSION(:, :), INTENT(IN), OPTIONAL :: ranges
92 : INTEGER, INTENT(IN), OPTIONAL :: nparticle, n_atom, n_core, n_shell
93 : INTEGER, INTENT(IN) :: iunit
94 : LOGICAL, INTENT(IN) :: print_atoms
95 :
96 : CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_create'
97 : #ifdef __SPGLIB
98 : CHARACTER(LEN=1000) :: buffer
99 : INTEGER :: ierr, nchars, nop, tra_mat(3, 3)
100 : #endif
101 : INTEGER :: handle, i, j, n_sr_rep
102 14 : INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp_types
103 : LOGICAL :: spglib
104 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_coor
105 : TYPE(spgr_type), POINTER :: spgr
106 :
107 14 : CALL timeset(routineN, handle)
108 :
109 14 : spgr => gopt_env%spgr
110 14 : CPASSERT(ASSOCIATED(spgr))
111 :
112 14 : CALL cleanup_spgr_type(spgr)
113 :
114 : !..total number of particles (atoms plus shells)
115 14 : IF (PRESENT(nparticle)) THEN
116 14 : CPASSERT(nparticle == SIZE(scoor, 2))
117 14 : spgr%nparticle = nparticle
118 : ELSE
119 0 : spgr%nparticle = SIZE(scoor, 2)
120 : END IF
121 :
122 14 : IF (PRESENT(n_atom)) THEN
123 14 : spgr%n_atom = n_atom
124 0 : ELSE IF (PRESENT(n_core)) THEN
125 0 : spgr%n_atom = spgr%nparticle - n_core
126 0 : ELSE IF (PRESENT(n_shell)) THEN
127 0 : spgr%n_atom = spgr%nparticle - n_shell
128 : ELSE
129 0 : spgr%n_atom = spgr%nparticle
130 : END IF
131 :
132 14 : IF (PRESENT(n_core)) THEN
133 14 : spgr%n_core = n_core
134 0 : ELSE IF (PRESENT(n_shell)) THEN
135 0 : spgr%n_core = n_shell
136 : END IF
137 :
138 14 : IF (PRESENT(n_shell)) THEN
139 14 : spgr%n_shell = n_shell
140 0 : ELSE IF (PRESENT(n_core)) THEN
141 0 : spgr%n_shell = n_core
142 : END IF
143 :
144 14 : IF (.NOT. (spgr%nparticle == (spgr%n_atom + spgr%n_shell))) THEN
145 0 : CPABORT("spgr_create: nparticle not equal to natom + nshell.")
146 : END IF
147 :
148 14 : spgr%nparticle_sym = spgr%nparticle
149 14 : spgr%n_atom_sym = spgr%n_atom
150 14 : spgr%n_core_sym = spgr%n_core
151 14 : spgr%n_shell_sym = spgr%n_shell
152 :
153 14 : spgr%iunit = iunit
154 14 : spgr%print_atoms = print_atoms
155 :
156 : ! accuracy for symmetry
157 14 : IF (PRESENT(eps_symmetry)) THEN
158 14 : spgr%eps_symmetry = eps_symmetry
159 : END IF
160 :
161 : ! vector to test reduced symmetry
162 14 : IF (PRESENT(pol)) THEN
163 14 : spgr%pol(1) = pol(1)
164 14 : spgr%pol(2) = pol(2)
165 14 : spgr%pol(3) = pol(3)
166 : END IF
167 :
168 42 : ALLOCATE (spgr%lat(spgr%nparticle))
169 170 : spgr%lat = .TRUE.
170 :
171 14 : IF (PRESENT(ranges)) THEN
172 0 : n_sr_rep = SIZE(ranges, 2)
173 0 : DO i = 1, n_sr_rep
174 0 : DO j = ranges(1, i), ranges(2, i)
175 0 : spgr%lat(j) = .FALSE.
176 0 : spgr%nparticle_sym = spgr%nparticle_sym - 1
177 0 : IF (j <= spgr%n_atom) THEN
178 0 : spgr%n_atom_sym = spgr%n_atom_sym - 1
179 0 : ELSE IF (j > spgr%n_atom .AND. j <= spgr%nparticle) THEN
180 0 : spgr%n_core_sym = spgr%n_core_sym - 1
181 0 : spgr%n_shell_sym = spgr%n_shell_sym - 1
182 : ELSE
183 0 : CPABORT("Symmetry exclusion range larger than actual number of particles.")
184 : END IF
185 : END DO
186 : END DO
187 : END IF
188 :
189 70 : ALLOCATE (tmp_coor(3, spgr%n_atom_sym), tmp_types(spgr%n_atom_sym))
190 :
191 14 : j = 0
192 122 : DO i = 1, spgr%n_atom
193 122 : IF (spgr%lat(i)) THEN
194 108 : j = j + 1
195 432 : tmp_coor(:, j) = scoor(:, i)
196 108 : tmp_types(j) = types(i)
197 : END IF
198 : END DO
199 :
200 : !..set cell values
201 14 : NULLIFY (spgr%cell_ref)
202 14 : CALL cell_create(spgr%cell_ref)
203 14 : CALL cell_copy(cell, spgr%cell_ref, tag="CELL_OPT_REF")
204 14 : SELECT CASE (gopt_env%type_id)
205 : CASE (default_minimization_method_id, default_ts_method_id)
206 0 : CALL init_cell(spgr%cell_ref, hmat=cell%hmat)
207 : CASE (default_cell_method_id)
208 28 : SELECT CASE (gopt_env%cell_method_id)
209 : CASE (default_cell_direct_id)
210 14 : CALL init_cell(spgr%cell_ref, hmat=gopt_env%h_ref)
211 : CASE (default_cell_geo_opt_id)
212 0 : CPABORT("SPACE_GROUP_SYMMETRY should not be invoked during the cell step.")
213 : CASE (default_cell_md_id)
214 0 : CPABORT("SPACE_GROUP_SYMMETRY is not compatible with md.")
215 : CASE DEFAULT
216 14 : CPABORT("SPACE_GROUP_SYMMETRY invoked with an unknown optimization method.")
217 : END SELECT
218 : CASE DEFAULT
219 14 : CPABORT("SPACE_GROUP_SYMMETRY is not compatible with md.")
220 : END SELECT
221 :
222 : ! atom types
223 42 : ALLOCATE (spgr%atype(spgr%nparticle))
224 170 : spgr%atype(1:spgr%nparticle) = types(1:spgr%nparticle)
225 :
226 14 : spgr%n_operations = 0
227 :
228 : #ifdef __SPGLIB
229 14 : spglib = .TRUE.
230 14 : CALL cite_reference(Togo2018)
231 : spgr%space_group_number = spg_get_international(spgr%international_symbol, TRANSPOSE(cell%hmat), tmp_coor, tmp_types, &
232 182 : spgr%n_atom_sym, eps_symmetry)
233 14 : buffer = ''
234 14 : nchars = strlcpy_c2f(buffer, spgr%international_symbol)
235 14 : spgr%international_symbol = buffer(1:nchars)
236 14 : IF (spgr%space_group_number == 0) THEN
237 0 : CPABORT("Symmetry Library SPGLIB failed, most likely due a problem with the coordinates.")
238 0 : spglib = .FALSE.
239 : ELSE
240 : nop = spg_get_multiplicity(TRANSPOSE(cell%hmat), tmp_coor, tmp_types, &
241 14 : spgr%n_atom_sym, eps_symmetry)
242 70 : ALLOCATE (spgr%rotations(3, 3, nop), spgr%translations(3, nop))
243 56 : ALLOCATE (spgr%eqatom(nop, spgr%nparticle))
244 42 : ALLOCATE (spgr%lop(nop))
245 14 : spgr%n_operations = nop
246 838 : spgr%lop = .TRUE.
247 : ierr = spg_get_symmetry(spgr%rotations, spgr%translations, nop, TRANSPOSE(cell%hmat), &
248 182 : tmp_coor, tmp_types, spgr%n_atom_sym, eps_symmetry)
249 : ! Schoenflies Symbol
250 : ierr = spg_get_schoenflies(spgr%schoenflies, TRANSPOSE(cell%hmat), tmp_coor, tmp_types, &
251 182 : spgr%n_atom_sym, eps_symmetry)
252 14 : buffer = ''
253 14 : nchars = strlcpy_c2f(buffer, spgr%schoenflies)
254 14 : spgr%schoenflies = buffer(1:nchars)
255 :
256 : ! Point Group
257 14 : tra_mat = 0
258 : ierr = spg_get_pointgroup(spgr%pointgroup_symbol, tra_mat, &
259 14 : spgr%rotations, spgr%n_operations)
260 14 : buffer = ''
261 14 : nchars = strlcpy_c2f(buffer, spgr%pointgroup_symbol)
262 14 : spgr%pointgroup_symbol = buffer(1:nchars)
263 : END IF
264 : #else
265 : CPABORT("Symmetry library SPGLIB not available")
266 : spglib = .FALSE.
267 : #endif
268 14 : spgr%symlib = spglib
269 :
270 14 : DEALLOCATE (tmp_coor, tmp_types)
271 :
272 14 : CALL timestop(handle)
273 :
274 14 : END SUBROUTINE spgr_create
275 :
276 : ! **************************************************************************************************
277 : !> \brief routine indentifies the space group and finds rotation matrices.
278 : !> \param subsys ...
279 : !> \param geo_section ...
280 : !> \param gopt_env ...
281 : !> \param iunit ...
282 : !> \par History
283 : !> 01.2020 created [pcazade]
284 : !> \author Pierre-André Cazade (first version)
285 : !> \note rotation matrices innclude translations and translation symmetry:
286 : !> it works with supercells as well.
287 : ! **************************************************************************************************
288 14 : SUBROUTINE identify_space_group(subsys, geo_section, gopt_env, iunit)
289 :
290 : TYPE(cp_subsys_type), INTENT(IN), POINTER :: subsys
291 : TYPE(section_vals_type), INTENT(IN), POINTER :: geo_section
292 : TYPE(gopt_f_type), INTENT(IN), POINTER :: gopt_env
293 : INTEGER, INTENT(IN) :: iunit
294 :
295 : CHARACTER(LEN=*), PARAMETER :: routineN = 'identify_space_group'
296 :
297 : INTEGER :: handle, i, k, n_atom, n_core, n_shell, &
298 : n_sr_rep, nparticle, shell_index
299 14 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atype
300 14 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ranges
301 14 : INTEGER, DIMENSION(:), POINTER :: tmp
302 : LOGICAL :: print_atoms
303 : REAL(KIND=dp) :: eps_symmetry
304 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: scoord
305 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: pol
306 : TYPE(cell_type), POINTER :: cell
307 : TYPE(particle_list_type), POINTER :: core_particles, particles, &
308 : shell_particles
309 : TYPE(spgr_type), POINTER :: spgr
310 :
311 14 : CALL timeset(routineN, handle)
312 :
313 14 : n_sr_rep = 0
314 14 : nparticle = 0
315 : n_atom = 0
316 14 : n_core = 0
317 14 : n_shell = 0
318 :
319 14 : NULLIFY (particles)
320 14 : NULLIFY (core_particles)
321 14 : NULLIFY (shell_particles)
322 :
323 : NULLIFY (cell)
324 14 : cell => subsys%cell
325 14 : CPASSERT(ASSOCIATED(cell))
326 :
327 14 : CALL cp_subsys_get(subsys, particles=particles, shell_particles=shell_particles, core_particles=core_particles)
328 :
329 14 : CPASSERT(ASSOCIATED(particles))
330 14 : n_atom = particles%n_els
331 : ! Check if we have other kinds of particles in this subsystem
332 14 : IF (ASSOCIATED(shell_particles)) THEN
333 4 : n_shell = shell_particles%n_els
334 4 : CPASSERT(ASSOCIATED(core_particles))
335 4 : n_core = subsys%core_particles%n_els
336 : ! The same number of shell and core particles is assumed
337 4 : CPASSERT(n_core == n_shell)
338 10 : ELSE IF (ASSOCIATED(core_particles)) THEN
339 : ! This case should not occur at the moment
340 0 : CPABORT("Core particles should not be defined without corresponding shell particles.")
341 : ELSE
342 : n_core = 0
343 : n_shell = 0
344 : END IF
345 :
346 14 : nparticle = n_atom + n_shell
347 70 : ALLOCATE (scoord(3, nparticle), atype(nparticle))
348 122 : DO i = 1, n_atom
349 108 : shell_index = particles%els(i)%shell_index
350 122 : IF (shell_index == 0) THEN
351 60 : CALL real_to_scaled(scoord(1:3, i), particles%els(i)%r(1:3), cell)
352 60 : CALL get_atomic_kind(atomic_kind=particles%els(i)%atomic_kind, kind_number=atype(i))
353 : ELSE
354 48 : CALL real_to_scaled(scoord(1:3, i), core_particles%els(shell_index)%r(1:3), cell)
355 48 : CALL get_atomic_kind(atomic_kind=core_particles%els(shell_index)%atomic_kind, kind_number=atype(i))
356 48 : k = n_atom + shell_index
357 48 : CALL real_to_scaled(scoord(1:3, k), shell_particles%els(shell_index)%r(1:3), cell)
358 48 : CALL get_atomic_kind(atomic_kind=shell_particles%els(shell_index)%atomic_kind, kind_number=atype(k))
359 : END IF
360 : END DO
361 :
362 14 : CALL section_vals_val_get(geo_section, "SPGR_PRINT_ATOMS", l_val=print_atoms)
363 14 : CALL section_vals_val_get(geo_section, "EPS_SYMMETRY", r_val=eps_symmetry)
364 14 : CALL section_vals_val_get(geo_section, "SYMM_REDUCTION", r_vals=pol)
365 14 : CALL section_vals_val_get(geo_section, "SYMM_EXCLUDE_RANGE", n_rep_val=n_sr_rep)
366 14 : IF (n_sr_rep > 0) THEN
367 0 : ALLOCATE (ranges(2, n_sr_rep))
368 0 : DO i = 1, n_sr_rep
369 0 : CALL section_vals_val_get(geo_section, "SYMM_EXCLUDE_RANGE", i_rep_val=i, i_vals=tmp)
370 0 : ranges(:, i) = tmp(:)
371 : END DO
372 : CALL spgr_create(scoord, atype, cell, gopt_env, eps_symmetry=eps_symmetry, pol=pol(1:3), &
373 : ranges=ranges, nparticle=nparticle, n_atom=n_atom, &
374 0 : n_core=n_core, n_shell=n_shell, iunit=iunit, print_atoms=print_atoms)
375 0 : DEALLOCATE (ranges)
376 : ELSE
377 : CALL spgr_create(scoord, atype, cell, gopt_env, eps_symmetry=eps_symmetry, pol=pol(1:3), &
378 : nparticle=nparticle, n_atom=n_atom, &
379 14 : n_core=n_core, n_shell=n_shell, iunit=iunit, print_atoms=print_atoms)
380 : END IF
381 :
382 : NULLIFY (spgr)
383 14 : spgr => gopt_env%spgr
384 :
385 14 : CALL spgr_find_equivalent_atoms(spgr, scoord)
386 14 : CALL spgr_reduce_symm(spgr)
387 14 : CALL spgr_rotations_subset(spgr)
388 :
389 14 : DEALLOCATE (scoord, atype)
390 :
391 14 : CALL timestop(handle)
392 :
393 56 : END SUBROUTINE identify_space_group
394 :
395 : ! **************************************************************************************************
396 : !> \brief routine indentifies the equivalent atoms for each rotation matrix.
397 : !> \param spgr ...
398 : !> \param scoord ...
399 : !> \par History
400 : !> 01.2020 created [pcazade]
401 : !> \author Pierre-André Cazade (first version)
402 : ! **************************************************************************************************
403 14 : SUBROUTINE spgr_find_equivalent_atoms(spgr, scoord)
404 :
405 : TYPE(spgr_type), INTENT(INOUT), POINTER :: spgr
406 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
407 : INTENT(IN) :: scoord
408 :
409 : CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_find_equivalent_atoms'
410 :
411 : INTEGER :: handle, i, ia, ib, ir, j, natom, nop, &
412 : nshell
413 : REAL(KIND=dp) :: diff
414 : REAL(KIND=dp), DIMENSION(3) :: rb, ri, ro, tr
415 : REAL(KIND=dp), DIMENSION(3, 3) :: rot
416 :
417 14 : CALL timeset(routineN, handle)
418 :
419 14 : nop = spgr%n_operations
420 14 : natom = spgr%n_atom
421 14 : nshell = spgr%n_shell
422 :
423 14 : IF (.NOT. (spgr%nparticle == (natom + nshell))) THEN
424 0 : CPABORT("spgr_find_equivalent_atoms: nparticle not equal to natom + nshell.")
425 : END IF
426 :
427 170 : DO ia = 1, spgr%nparticle
428 3906 : spgr%eqatom(:, ia) = ia
429 : END DO
430 :
431 14 : !$OMP PARALLEL DO PRIVATE (ia,ib,ir,ri,rb,ro,rot,tr,diff) SHARED (spgr,scoord,natom,nop) DEFAULT(NONE)
432 : DO ia = 1, natom
433 : IF (.NOT. spgr%lat(ia)) CYCLE
434 : ri(1:3) = scoord(1:3, ia)
435 : DO ir = 1, nop
436 : rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
437 : tr(1:3) = spgr%translations(1:3, ir)
438 : DO ib = 1, natom
439 : IF (.NOT. spgr%lat(ib)) CYCLE
440 : rb(1:3) = scoord(1:3, ib)
441 : ro(1) = REAL(rot(1, 1), dp)*rb(1) + REAL(rot(2, 1), dp)*rb(2) + REAL(rot(3, 1), dp)*rb(3) + tr(1)
442 : ro(2) = REAL(rot(1, 2), dp)*rb(1) + REAL(rot(2, 2), dp)*rb(2) + REAL(rot(3, 2), dp)*rb(3) + tr(2)
443 : ro(3) = REAL(rot(1, 3), dp)*rb(1) + REAL(rot(2, 3), dp)*rb(2) + REAL(rot(3, 3), dp)*rb(3) + tr(3)
444 : ro(1) = ro(1) - REAL(NINT(ro(1) - ri(1)), dp)
445 : ro(2) = ro(2) - REAL(NINT(ro(2) - ri(2)), dp)
446 : ro(3) = ro(3) - REAL(NINT(ro(3) - ri(3)), dp)
447 : diff = NORM2(ri(:) - ro(:))
448 : IF ((diff < spgr%eps_symmetry) .AND. (spgr%atype(ia) == spgr%atype(ib))) THEN
449 : spgr%eqatom(ir, ia) = ib
450 : EXIT
451 : END IF
452 : END DO
453 : END DO
454 : END DO
455 : !$OMP END PARALLEL DO
456 :
457 14 : !$OMP PARALLEL DO PRIVATE (i,j,ia,ib,ir,ri,rb,ro,rot,tr,diff) SHARED (spgr,scoord,natom,nshell,nop) DEFAULT(NONE)
458 : DO i = 1, nshell
459 : ia = natom + i
460 : IF (.NOT. spgr%lat(ia)) CYCLE
461 : ri(1:3) = scoord(1:3, ia)
462 : DO ir = 1, nop
463 : rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
464 : tr(1:3) = spgr%translations(1:3, ir)
465 : DO j = 1, nshell
466 : ib = natom + j
467 : IF (.NOT. spgr%lat(ib)) CYCLE
468 : rb(1:3) = scoord(1:3, ib)
469 : ro(1) = REAL(rot(1, 1), dp)*rb(1) + REAL(rot(2, 1), dp)*rb(2) + REAL(rot(3, 1), dp)*rb(3) + tr(1)
470 : ro(2) = REAL(rot(1, 2), dp)*rb(1) + REAL(rot(2, 2), dp)*rb(2) + REAL(rot(3, 2), dp)*rb(3) + tr(2)
471 : ro(3) = REAL(rot(1, 3), dp)*rb(1) + REAL(rot(2, 3), dp)*rb(2) + REAL(rot(3, 3), dp)*rb(3) + tr(3)
472 : ro(1) = ro(1) - REAL(NINT(ro(1) - ri(1)), dp)
473 : ro(2) = ro(2) - REAL(NINT(ro(2) - ri(2)), dp)
474 : ro(3) = ro(3) - REAL(NINT(ro(3) - ri(3)), dp)
475 : diff = NORM2(ri(:) - ro(:))
476 : IF ((diff < spgr%eps_symmetry) .AND. (spgr%atype(ia) == spgr%atype(ib))) THEN
477 : spgr%eqatom(ir, ia) = ib
478 : EXIT
479 : END IF
480 : END DO
481 : END DO
482 : END DO
483 : !$OMP END PARALLEL DO
484 :
485 14 : CALL timestop(handle)
486 :
487 14 : END SUBROUTINE spgr_find_equivalent_atoms
488 :
489 : ! **************************************************************************************************
490 : !> \brief routine looks for operations compatible with efield
491 : !> \param spgr ...
492 : !> \par History
493 : !> 01.2020 created [pcazade]
494 : !> \author Pierre-André Cazade (first version)
495 : ! **************************************************************************************************
496 14 : SUBROUTINE spgr_reduce_symm(spgr)
497 :
498 : TYPE(spgr_type), INTENT(INOUT), POINTER :: spgr
499 :
500 : CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_reduce_symm'
501 :
502 : INTEGER :: handle, ia, ib, ir, ja, jb, nop, nops, &
503 : nparticle
504 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x, xold
505 : REAL(KIND=dp), DIMENSION(3) :: ri, ro
506 : REAL(KIND=dp), DIMENSION(3, 3) :: rot
507 :
508 14 : CALL timeset(routineN, handle)
509 :
510 14 : nop = spgr%n_operations
511 14 : nparticle = spgr%nparticle
512 56 : ALLOCATE (x(3*nparticle), xold(3*nparticle))
513 482 : x = 0.0_dp
514 170 : DO ia = 1, nparticle
515 156 : ja = 3*(ia - 1)
516 156 : x(ja + 1) = x(ja + 1) + spgr%pol(1)
517 156 : x(ja + 2) = x(ja + 2) + spgr%pol(2)
518 170 : x(ja + 3) = x(ja + 3) + spgr%pol(3)
519 : END DO
520 482 : xold(:) = x(:)
521 :
522 : nops = 0
523 838 : DO ir = 1, nop
524 12032 : x = 0.d0
525 824 : spgr%lop(ir) = .TRUE.
526 10712 : rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
527 4560 : DO ia = 1, nparticle
528 3736 : IF (.NOT. spgr%lat(ia)) CYCLE
529 3736 : ja = 3*(ia - 1)
530 14944 : ri(1:3) = xold(ja + 1:ja + 3)
531 3736 : ro(1) = REAL(rot(1, 1), dp)*ri(1) + REAL(rot(2, 1), dp)*ri(2) + REAL(rot(3, 1), dp)*ri(3)
532 3736 : ro(2) = REAL(rot(1, 2), dp)*ri(1) + REAL(rot(2, 2), dp)*ri(2) + REAL(rot(3, 2), dp)*ri(3)
533 3736 : ro(3) = REAL(rot(1, 3), dp)*ri(1) + REAL(rot(2, 3), dp)*ri(2) + REAL(rot(3, 3), dp)*ri(3)
534 15768 : x(ja + 1:ja + 3) = ro(1:3)
535 : END DO
536 4560 : DO ia = 1, nparticle
537 3736 : IF (.NOT. spgr%lat(ia)) CYCLE
538 3736 : ib = spgr%eqatom(ir, ia)
539 3736 : ja = 3*(ia - 1)
540 3736 : jb = 3*(ib - 1)
541 14944 : ro = x(jb + 1:jb + 3) - xold(ja + 1:ja + 3)
542 : spgr%lop(ir) = (spgr%lop(ir) .AND. (ABS(ro(1)) < spgr%eps_symmetry) &
543 : .AND. (ABS(ro(2)) < spgr%eps_symmetry) &
544 4560 : .AND. (ABS(ro(3)) < spgr%eps_symmetry))
545 : END DO
546 838 : IF (spgr%lop(ir)) nops = nops + 1
547 : END DO
548 :
549 14 : spgr%n_reduced_operations = nops
550 :
551 14 : DEALLOCATE (x, xold)
552 14 : CALL timestop(handle)
553 :
554 14 : END SUBROUTINE spgr_reduce_symm
555 :
556 : ! **************************************************************************************************
557 : !> \brief routine looks for unique rotations
558 : !> \param spgr ...
559 : !> \par History
560 : !> 01.2020 created [pcazade]
561 : !> \author Pierre-André Cazade (first version)
562 : ! **************************************************************************************************
563 :
564 14 : SUBROUTINE spgr_rotations_subset(spgr)
565 :
566 : TYPE(spgr_type), INTENT(INOUT), POINTER :: spgr
567 :
568 : CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_rotations_subset'
569 :
570 : INTEGER :: handle, i, j
571 : INTEGER, DIMENSION(3, 3) :: d
572 14 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: mask
573 :
574 14 : CALL timeset(routineN, handle)
575 :
576 42 : ALLOCATE (mask(spgr%n_operations))
577 838 : mask = .TRUE.
578 :
579 838 : DO i = 1, spgr%n_operations
580 838 : IF (.NOT. spgr%lop(i)) mask(i) = .FALSE.
581 : END DO
582 :
583 824 : DO i = 1, spgr%n_operations - 1
584 810 : IF (.NOT. mask(i)) CYCLE
585 32552 : DO j = i + 1, spgr%n_operations
586 32300 : IF (.NOT. mask(j)) CYCLE
587 243932 : d(:, :) = spgr%rotations(:, :, j) - spgr%rotations(:, :, i)
588 244742 : IF (SUM(ABS(d)) == 0) mask(j) = .FALSE.
589 : END DO
590 : END DO
591 :
592 14 : spgr%n_operations_subset = 0
593 838 : DO i = 1, spgr%n_operations
594 838 : IF (mask(i)) spgr%n_operations_subset = spgr%n_operations_subset + 1
595 : END DO
596 :
597 42 : ALLOCATE (spgr%rotations_subset(3, 3, spgr%n_operations_subset))
598 :
599 14 : j = 0
600 838 : DO i = 1, spgr%n_operations
601 838 : IF (mask(i)) THEN
602 248 : j = j + 1
603 3224 : spgr%rotations_subset(:, :, j) = spgr%rotations(:, :, i)
604 : END IF
605 : END DO
606 :
607 14 : DEALLOCATE (mask)
608 14 : CALL timestop(handle)
609 :
610 14 : END SUBROUTINE spgr_rotations_subset
611 :
612 : ! **************************************************************************************************
613 : !> \brief routine applies the rotation matrices to the coordinates.
614 : !> \param spgr ...
615 : !> \param coord ...
616 : !> \par History
617 : !> 01.2020 created [pcazade]
618 : !> \author Pierre-André Cazade (first version)
619 : ! **************************************************************************************************
620 46 : SUBROUTINE spgr_apply_rotations_coord(spgr, coord)
621 :
622 : TYPE(spgr_type), INTENT(IN), POINTER :: spgr
623 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: coord
624 :
625 : CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_apply_rotations_coord'
626 :
627 : INTEGER :: handle, ia, ib, ir, ja, jb, nop, nops, &
628 : nparticle
629 46 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cold
630 : REAL(KIND=dp), DIMENSION(3) :: rf, ri, rn, ro, tr
631 : REAL(KIND=dp), DIMENSION(3, 3) :: rot
632 :
633 46 : CALL timeset(routineN, handle)
634 :
635 138 : ALLOCATE (cold(SIZE(coord)))
636 2074 : cold(:) = coord(:)
637 :
638 46 : nop = spgr%n_operations
639 46 : nparticle = spgr%nparticle
640 46 : nops = spgr%n_reduced_operations
641 :
642 46 : !$OMP PARALLEL DO PRIVATE (ia,ib,ja,jb,ir,ri,ro,rf,rn,rot,tr) SHARED (spgr,coord,nparticle,nop,nops) DEFAULT(NONE)
643 : DO ia = 1, nparticle
644 : IF (.NOT. spgr%lat(ia)) CYCLE
645 : ja = 3*(ia - 1)
646 : CALL real_to_scaled(rf(1:3), coord(ja + 1:ja + 3), spgr%cell_ref)
647 : rn(1:3) = 0.d0
648 : DO ir = 1, nop
649 : IF (.NOT. spgr%lop(ir)) CYCLE
650 : ib = spgr%eqatom(ir, ia)
651 : rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
652 : tr(1:3) = spgr%translations(1:3, ir)
653 : jb = 3*(ib - 1)
654 : CALL real_to_scaled(ri(1:3), coord(jb + 1:jb + 3), spgr%cell_ref)
655 : ro(1) = REAL(rot(1, 1), dp)*ri(1) + REAL(rot(2, 1), dp)*ri(2) + REAL(rot(3, 1), dp)*ri(3) + tr(1)
656 : ro(2) = REAL(rot(1, 2), dp)*ri(1) + REAL(rot(2, 2), dp)*ri(2) + REAL(rot(3, 2), dp)*ri(3) + tr(2)
657 : ro(3) = REAL(rot(1, 3), dp)*ri(1) + REAL(rot(2, 3), dp)*ri(2) + REAL(rot(3, 3), dp)*ri(3) + tr(3)
658 : ro(1) = ro(1) - REAL(NINT(ro(1) - rf(1)), dp)
659 : ro(2) = ro(2) - REAL(NINT(ro(2) - rf(2)), dp)
660 : ro(3) = ro(3) - REAL(NINT(ro(3) - rf(3)), dp)
661 : rn(1:3) = rn(1:3) + ro(1:3)
662 : END DO
663 : rn = rn/REAL(nops, dp)
664 : CALL scaled_to_real(coord(ja + 1:ja + 3), rn(1:3), spgr%cell_ref)
665 : END DO
666 : !$OMP END PARALLEL DO
667 :
668 46 : DEALLOCATE (cold)
669 46 : CALL timestop(handle)
670 :
671 46 : END SUBROUTINE spgr_apply_rotations_coord
672 :
673 : ! **************************************************************************************************
674 : !> \brief routine applies the rotation matrices to the forces.
675 : !> \param spgr ...
676 : !> \param force ...
677 : !> \par History
678 : !> 01.2020 created [pcazade]
679 : !> \author Pierre-André Cazade (first version)
680 : ! **************************************************************************************************
681 57 : SUBROUTINE spgr_apply_rotations_force(spgr, force)
682 :
683 : TYPE(spgr_type), INTENT(IN), POINTER :: spgr
684 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: force
685 :
686 : CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_apply_rotations_force'
687 :
688 : INTEGER :: handle, ia, ib, ir, ja, jb, nop, nops, &
689 : nparticle
690 57 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fold
691 : REAL(KIND=dp), DIMENSION(3) :: ri, rn, ro
692 : REAL(KIND=dp), DIMENSION(3, 3) :: rot
693 :
694 57 : CALL timeset(routineN, handle)
695 :
696 171 : ALLOCATE (fold(SIZE(force)))
697 2943 : fold(:) = force(:)
698 :
699 57 : nop = spgr%n_operations
700 57 : nparticle = spgr%nparticle
701 57 : nops = spgr%n_reduced_operations
702 :
703 57 : !$OMP PARALLEL DO PRIVATE (ia,ib,ja,jb,ir,ri,ro,rn,rot) SHARED (spgr,force,nparticle,nop,nops) DEFAULT(NONE)
704 : DO ia = 1, nparticle
705 : IF (.NOT. spgr%lat(ia)) CYCLE
706 : ja = 3*(ia - 1)
707 : rn(1:3) = 0.d0
708 : DO ir = 1, nop
709 : IF (.NOT. spgr%lop(ir)) CYCLE
710 : ib = spgr%eqatom(ir, ia)
711 : rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
712 : jb = 3*(ib - 1)
713 : CALL real_to_scaled(ri(1:3), force(jb + 1:jb + 3), spgr%cell_ref)
714 : ro(1) = REAL(rot(1, 1), dp)*ri(1) + REAL(rot(2, 1), dp)*ri(2) + REAL(rot(3, 1), dp)*ri(3)
715 : ro(2) = REAL(rot(1, 2), dp)*ri(1) + REAL(rot(2, 2), dp)*ri(2) + REAL(rot(3, 2), dp)*ri(3)
716 : ro(3) = REAL(rot(1, 3), dp)*ri(1) + REAL(rot(2, 3), dp)*ri(2) + REAL(rot(3, 3), dp)*ri(3)
717 : rn(1:3) = rn(1:3) + ro(1:3)
718 : END DO
719 : rn = rn/REAL(nops, dp)
720 : CALL scaled_to_real(force(ja + 1:ja + 3), rn(1:3), spgr%cell_ref)
721 : END DO
722 : !$OMP END PARALLEL DO
723 :
724 57 : DEALLOCATE (fold)
725 57 : CALL timestop(handle)
726 :
727 57 : END SUBROUTINE spgr_apply_rotations_force
728 :
729 : ! **************************************************************************************************
730 : !> \brief ...
731 : !> \param roti ...
732 : !> \param roto ...
733 : !> \param nop ...
734 : !> \param h1 ...
735 : !> \param h2 ...
736 : ! **************************************************************************************************
737 20 : SUBROUTINE spgr_change_basis(roti, roto, nop, h1, h2)
738 :
739 : INTEGER, DIMENSION(:, :, :) :: roti
740 : REAL(KIND=dp), DIMENSION(:, :, :) :: roto
741 : INTEGER :: nop
742 : REAL(KIND=dp), DIMENSION(3, 3) :: h1, h2
743 :
744 : CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_change_basis'
745 :
746 : INTEGER :: handle, ir
747 : REAL(KIND=dp), DIMENSION(3, 3) :: h1ih2, h2ih1, ih1, ih2, r, s
748 :
749 20 : CALL timeset(routineN, handle)
750 :
751 20 : ih1 = inv_3x3(h1)
752 20 : ih2 = inv_3x3(h2)
753 800 : h2ih1 = MATMUL(h2, ih1)
754 800 : h1ih2 = MATMUL(h1, ih2)
755 :
756 196 : DO ir = 1, nop
757 2288 : r(:, :) = roti(:, :, ir)
758 7040 : s = MATMUL(h2ih1, r)
759 7040 : r = MATMUL(s, h1ih2)
760 2308 : roto(:, :, ir) = r(:, :)
761 : END DO
762 :
763 20 : CALL timestop(handle)
764 :
765 20 : END SUBROUTINE spgr_change_basis
766 :
767 : ! **************************************************************************************************
768 : !> \brief routine applies the rotation matrices to the stress tensor.
769 : !> \param spgr ...
770 : !> \param cell ...
771 : !> \param stress ...
772 : !> \par History
773 : !> 01.2020 created [pcazade]
774 : !> \author Pierre-André Cazade (first version)
775 : ! **************************************************************************************************
776 20 : SUBROUTINE spgr_apply_rotations_stress(spgr, cell, stress)
777 :
778 : TYPE(spgr_type), INTENT(IN), POINTER :: spgr
779 : TYPE(cell_type), INTENT(IN), POINTER :: cell
780 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT) :: stress
781 :
782 : CHARACTER(LEN=*), PARAMETER :: routineN = 'spgr_apply_rotations_stress'
783 :
784 : INTEGER :: handle, i, ir, j, k, l, nop
785 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: roto
786 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat1, hmat2, r, stin
787 :
788 20 : CALL timeset(routineN, handle)
789 :
790 260 : hmat1 = TRANSPOSE(cell%hmat)
791 :
792 20 : hmat2 = 0d0
793 20 : hmat2(1, 1) = 1.d0
794 20 : hmat2(2, 2) = 1.d0
795 20 : hmat2(3, 3) = 1.d0
796 :
797 20 : nop = spgr%n_operations_subset
798 :
799 60 : ALLOCATE (roto(3, 3, nop))
800 :
801 20 : CALL spgr_change_basis(spgr%rotations_subset, roto, spgr%n_operations_subset, hmat1, hmat2)
802 :
803 20 : stin = stress
804 20 : stress = 0.d0
805 196 : DO ir = 1, nop
806 2288 : r(:, :) = roto(:, :, ir)
807 724 : DO i = 1, 3
808 2288 : DO j = 1, 3
809 6864 : DO k = 1, 3
810 20592 : DO l = 1, 3
811 19008 : stress(i, j) = stress(i, j) + (r(k, i)*r(l, j)*stin(k, l))
812 : END DO
813 : END DO
814 : END DO
815 : END DO
816 : END DO
817 260 : stress = stress/REAL(nop, dp)
818 :
819 20 : DEALLOCATE (roto)
820 :
821 20 : CALL timestop(handle)
822 :
823 20 : END SUBROUTINE spgr_apply_rotations_stress
824 :
825 : ! **************************************************************************************************
826 : !> \brief routine prints Space Group Information.
827 : !> \param spgr ...
828 : !> \par History
829 : !> 01.2020 created [pcazade]
830 : !> \author Pierre-André Cazade (first version)
831 : ! **************************************************************************************************
832 14 : SUBROUTINE print_spgr(spgr)
833 :
834 : TYPE(spgr_type), INTENT(IN), POINTER :: spgr
835 :
836 : INTEGER :: i, j
837 :
838 14 : IF (spgr%iunit > 0) THEN
839 7 : WRITE (spgr%iunit, '(/,T2,A,A)') "----------------------------------------", &
840 14 : "---------------------------------------"
841 7 : WRITE (spgr%iunit, "(T2,A,T25,A,T77,A)") "----", "SPACE GROUP SYMMETRY INFORMATION", "----"
842 7 : WRITE (spgr%iunit, '(T2,A,A)') "----------------------------------------", &
843 14 : "---------------------------------------"
844 7 : IF (spgr%symlib) THEN
845 7 : WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| SPACE GROUP NUMBER:", &
846 14 : spgr%space_group_number
847 7 : WRITE (spgr%iunit, '(T2,A,T70,A11)') "SPGR| INTERNATIONAL SYMBOL:", &
848 14 : TRIM(ADJUSTR(spgr%international_symbol))
849 7 : WRITE (spgr%iunit, '(T2,A,T75,A6)') "SPGR| POINT GROUP SYMBOL:", &
850 14 : TRIM(ADJUSTR(spgr%pointgroup_symbol))
851 7 : WRITE (spgr%iunit, '(T2,A,T74,A7)') "SPGR| SCHOENFLIES SYMBOL:", &
852 14 : TRIM(ADJUSTR(spgr%schoenflies))
853 7 : WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF SYMMETRY OPERATIONS:", &
854 14 : spgr%n_operations
855 7 : WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF UNIQUE ROTATIONS:", &
856 14 : spgr%n_operations_subset
857 7 : WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF REDUCED SYMMETRY OPERATIONS:", &
858 14 : spgr%n_reduced_operations
859 7 : WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF PARTICLES AND SYMMETRIZED PARTICLES:", &
860 14 : spgr%nparticle, spgr%nparticle_sym
861 7 : WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF ATOMS AND SYMMETRIZED ATOMS:", &
862 14 : spgr%n_atom, spgr%n_atom_sym
863 7 : WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF CORES AND SYMMETRIZED CORES:", &
864 14 : spgr%n_core, spgr%n_core_sym
865 7 : WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF SHELLS AND SYMMETRIZED SHELLS:", &
866 14 : spgr%n_shell, spgr%n_shell_sym
867 7 : IF (spgr%print_atoms) THEN
868 1 : WRITE (spgr%iunit, *) "SPGR| ACTIVE REDUCED SYMMETRY OPERATIONS:", spgr%lop
869 1 : WRITE (spgr%iunit, '(/,T2,A,A)') "----------------------------------------", &
870 2 : "---------------------------------------"
871 1 : WRITE (spgr%iunit, '(T2,A,T34,A,T77,A)') "----", "EQUIVALENT ATOMS", "----"
872 1 : WRITE (spgr%iunit, '(T2,A,A)') "----------------------------------------", &
873 2 : "---------------------------------------"
874 25 : DO i = 1, spgr%nparticle
875 121 : DO j = 1, spgr%n_operations
876 96 : WRITE (spgr%iunit, '(T2,A,T52,I8,I8,I8)') "SPGR| ATOM | SYMMETRY OPERATION | EQUIVALENT ATOM", &
877 216 : i, j, spgr%eqatom(j, i)
878 : END DO
879 : END DO
880 1 : WRITE (spgr%iunit, '(T2,A,A)') "----------------------------------------", &
881 2 : "---------------------------------------"
882 5 : DO i = 1, spgr%n_operations
883 : WRITE (spgr%iunit, '(T2,A,T46,i4,T51,3I10,/,T51,3I10,/,T51,3I10)') &
884 16 : "SPGR| SYMMETRY OPERATION #:", i, (spgr%rotations(j, :, i), j=1, 3)
885 5 : WRITE (spgr%iunit, '(T51,3F10.5)') spgr%translations(:, i)
886 : END DO
887 : END IF
888 : ELSE
889 0 : WRITE (spgr%iunit, "(T2,A)") "SPGLIB for Crystal Symmetry Information determination is not availale"
890 : END IF
891 : END IF
892 :
893 14 : END SUBROUTINE print_spgr
894 :
895 : ! **************************************************************************************************
896 : !> \brief Variable precision output of the symmetrized stress tensor
897 : !>
898 : !> \param stress tensor ...
899 : !> \param spgr ...
900 : !> \par History
901 : !> 07.2020 adapted to spgr [pcazade]
902 : !> \author MK (26.08.2010).
903 : ! **************************************************************************************************
904 20 : SUBROUTINE spgr_write_stress_tensor(stress, spgr)
905 :
906 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: stress
907 : TYPE(spgr_type), INTENT(IN), POINTER :: spgr
908 :
909 : REAL(KIND=dp), DIMENSION(3) :: eigval
910 : REAL(KIND=dp), DIMENSION(3, 3) :: eigvec, stress_tensor
911 :
912 260 : stress_tensor(:, :) = stress(:, :)*pascal*1.0E-9_dp
913 :
914 20 : IF (spgr%iunit > 0) THEN
915 : WRITE (UNIT=spgr%iunit, FMT='(/,T2,A)') &
916 11 : 'SPGR STRESS| Symmetrized stress tensor [GPa]'
917 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T19,3(19X,A1))') &
918 11 : 'SPGR STRESS|', 'x', 'y', 'z'
919 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
920 11 : 'SPGR STRESS| x', stress_tensor(1, 1:3)
921 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
922 11 : 'SPGR STRESS| y', stress_tensor(2, 1:3)
923 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
924 11 : 'SPGR STRESS| z', stress_tensor(3, 1:3)
925 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T66,ES20.11)') &
926 11 : 'SPGR STRESS| 1/3 Trace', (stress_tensor(1, 1) + &
927 : stress_tensor(2, 2) + &
928 22 : stress_tensor(3, 3))/3.0_dp
929 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T66,ES20.11)') &
930 11 : 'SPGR STRESS| Determinant', det_3x3(stress_tensor(1:3, 1), &
931 : stress_tensor(1:3, 2), &
932 22 : stress_tensor(1:3, 3))
933 11 : eigval(:) = 0.0_dp
934 11 : eigvec(:, :) = 0.0_dp
935 11 : CALL jacobi(stress_tensor, eigval, eigvec)
936 : WRITE (UNIT=spgr%iunit, FMT='(/,T2,A)') &
937 11 : 'SPGR STRESS| Eigenvectors and eigenvalues of the symmetrized stress tensor [GPa]'
938 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T19,3(1X,I19))') &
939 11 : 'SPGR STRESS|', 1, 2, 3
940 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
941 11 : 'SPGR STRESS| Eigenvalues', eigval(1:3)
942 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,F19.12))') &
943 11 : 'SPGR STRESS| x', eigvec(1, 1:3)
944 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,F19.12))') &
945 11 : 'SPGR STRESS| y', eigvec(2, 1:3)
946 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,F19.12))') &
947 11 : 'SPGR STRESS| z', eigvec(3, 1:3)
948 : END IF
949 :
950 20 : END SUBROUTINE spgr_write_stress_tensor
951 :
952 : END MODULE space_groups
|
