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