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