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 16 : 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 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp_types
103 : LOGICAL :: spglib
104 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_coor
105 : TYPE(spgr_type), POINTER :: spgr
106 :
107 16 : CALL timeset(routineN, handle)
108 :
109 16 : spgr => gopt_env%spgr
110 16 : CPASSERT(ASSOCIATED(spgr))
111 :
112 16 : CALL cleanup_spgr_type(spgr)
113 :
114 : !..total number of particles (atoms plus shells)
115 16 : IF (PRESENT(nparticle)) THEN
116 16 : CPASSERT(nparticle == SIZE(scoor, 2))
117 16 : spgr%nparticle = nparticle
118 : ELSE
119 0 : spgr%nparticle = SIZE(scoor, 2)
120 : END IF
121 :
122 16 : IF (PRESENT(n_atom)) THEN
123 16 : 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 16 : IF (PRESENT(n_core)) THEN
133 16 : 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 16 : IF (PRESENT(n_shell)) THEN
139 16 : 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 16 : 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 16 : spgr%nparticle_sym = spgr%nparticle
149 16 : spgr%n_atom_sym = spgr%n_atom
150 16 : spgr%n_core_sym = spgr%n_core
151 16 : spgr%n_shell_sym = spgr%n_shell
152 :
153 16 : spgr%iunit = iunit
154 16 : spgr%print_atoms = print_atoms
155 :
156 : ! accuracy for symmetry
157 16 : IF (PRESENT(eps_symmetry)) THEN
158 16 : spgr%eps_symmetry = eps_symmetry
159 : END IF
160 :
161 : ! vector to test reduced symmetry
162 16 : IF (PRESENT(pol)) THEN
163 16 : spgr%pol(1) = pol(1)
164 16 : spgr%pol(2) = pol(2)
165 16 : spgr%pol(3) = pol(3)
166 : END IF
167 :
168 48 : ALLOCATE (spgr%lat(spgr%nparticle))
169 190 : spgr%lat = .TRUE.
170 :
171 16 : 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 80 : ALLOCATE (tmp_coor(3, spgr%n_atom_sym), tmp_types(spgr%n_atom_sym))
190 :
191 16 : j = 0
192 142 : DO i = 1, spgr%n_atom
193 142 : IF (spgr%lat(i)) THEN
194 126 : j = j + 1
195 504 : tmp_coor(:, j) = scoor(:, i)
196 126 : tmp_types(j) = types(i)
197 : END IF
198 : END DO
199 :
200 : !..set cell values
201 16 : NULLIFY (spgr%cell_ref)
202 16 : CALL cell_create(spgr%cell_ref)
203 16 : CALL cell_copy(cell, spgr%cell_ref, tag="CELL_OPT_REF")
204 16 : 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 32 : SELECT CASE (gopt_env%cell_method_id)
209 : CASE (default_cell_direct_id)
210 16 : 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 16 : CPABORT("SPACE_GROUP_SYMMETRY invoked with an unknown optimization method.")
217 : END SELECT
218 : CASE DEFAULT
219 16 : CPABORT("SPACE_GROUP_SYMMETRY is not compatible with md.")
220 : END SELECT
221 :
222 : ! atom types
223 48 : ALLOCATE (spgr%atype(spgr%nparticle))
224 190 : spgr%atype(1:spgr%nparticle) = types(1:spgr%nparticle)
225 :
226 16 : spgr%n_operations = 0
227 :
228 : #ifdef __SPGLIB
229 16 : spglib = .TRUE.
230 16 : CALL cite_reference(Togo2018)
231 : spgr%space_group_number = spg_get_international(spgr%international_symbol, TRANSPOSE(cell%hmat), tmp_coor, tmp_types, &
232 208 : spgr%n_atom_sym, eps_symmetry)
233 16 : buffer = ''
234 16 : nchars = strlcpy_c2f(buffer, spgr%international_symbol)
235 16 : spgr%international_symbol = buffer(1:nchars)
236 16 : 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 16 : spgr%n_atom_sym, eps_symmetry)
242 80 : ALLOCATE (spgr%rotations(3, 3, nop), spgr%translations(3, nop))
243 64 : ALLOCATE (spgr%eqatom(nop, spgr%nparticle))
244 48 : ALLOCATE (spgr%lop(nop))
245 16 : spgr%n_operations = nop
246 852 : spgr%lop = .TRUE.
247 : ierr = spg_get_symmetry(spgr%rotations, spgr%translations, nop, TRANSPOSE(cell%hmat), &
248 208 : 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 208 : spgr%n_atom_sym, eps_symmetry)
252 16 : buffer = ''
253 16 : nchars = strlcpy_c2f(buffer, spgr%schoenflies)
254 16 : spgr%schoenflies = buffer(1:nchars)
255 :
256 : ! Point Group
257 16 : tra_mat = 0
258 : ierr = spg_get_pointgroup(spgr%pointgroup_symbol, tra_mat, &
259 16 : spgr%rotations, spgr%n_operations)
260 16 : buffer = ''
261 16 : nchars = strlcpy_c2f(buffer, spgr%pointgroup_symbol)
262 16 : spgr%pointgroup_symbol = buffer(1:nchars)
263 : END IF
264 : #else
265 : CPABORT("Symmetry library SPGLIB not available")
266 : spglib = .FALSE.
267 : #endif
268 16 : spgr%symlib = spglib
269 :
270 16 : DEALLOCATE (tmp_coor, tmp_types)
271 :
272 16 : CALL timestop(handle)
273 :
274 16 : 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 16 : 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 16 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atype
300 16 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ranges
301 16 : INTEGER, DIMENSION(:), POINTER :: tmp
302 : LOGICAL :: print_atoms
303 : REAL(KIND=dp) :: eps_symmetry
304 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: scoord
305 16 : 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 16 : CALL timeset(routineN, handle)
312 :
313 16 : n_sr_rep = 0
314 16 : nparticle = 0
315 : n_atom = 0
316 16 : n_core = 0
317 16 : n_shell = 0
318 :
319 16 : NULLIFY (particles)
320 16 : NULLIFY (core_particles)
321 16 : NULLIFY (shell_particles)
322 :
323 : NULLIFY (cell)
324 16 : cell => subsys%cell
325 16 : CPASSERT(ASSOCIATED(cell))
326 :
327 16 : CALL cp_subsys_get(subsys, particles=particles, shell_particles=shell_particles, core_particles=core_particles)
328 :
329 16 : CPASSERT(ASSOCIATED(particles))
330 16 : n_atom = particles%n_els
331 : ! Check if we have other kinds of particles in this subsystem
332 16 : 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 12 : 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 16 : nparticle = n_atom + n_shell
347 80 : ALLOCATE (scoord(3, nparticle), atype(nparticle))
348 142 : DO i = 1, n_atom
349 126 : shell_index = particles%els(i)%shell_index
350 142 : IF (shell_index == 0) THEN
351 78 : CALL real_to_scaled(scoord(1:3, i), particles%els(i)%r(1:3), cell)
352 78 : 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 16 : CALL section_vals_val_get(geo_section, "SPGR_PRINT_ATOMS", l_val=print_atoms)
363 16 : CALL section_vals_val_get(geo_section, "EPS_SYMMETRY", r_val=eps_symmetry)
364 16 : CALL section_vals_val_get(geo_section, "SYMM_REDUCTION", r_vals=pol)
365 16 : CALL section_vals_val_get(geo_section, "SYMM_EXCLUDE_RANGE", n_rep_val=n_sr_rep)
366 16 : 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 16 : n_core=n_core, n_shell=n_shell, iunit=iunit, print_atoms=print_atoms)
380 : END IF
381 :
382 : NULLIFY (spgr)
383 16 : spgr => gopt_env%spgr
384 :
385 16 : CALL spgr_find_equivalent_atoms(spgr, scoord)
386 16 : CALL spgr_reduce_symm(spgr)
387 16 : CALL spgr_rotations_subset(spgr)
388 :
389 16 : DEALLOCATE (scoord, atype)
390 :
391 16 : CALL timestop(handle)
392 :
393 64 : 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 16 : 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 16 : CALL timeset(routineN, handle)
418 :
419 16 : nop = spgr%n_operations
420 16 : natom = spgr%n_atom
421 16 : nshell = spgr%n_shell
422 :
423 16 : 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 190 : DO ia = 1, spgr%nparticle
428 4034 : spgr%eqatom(:, ia) = ia
429 : END DO
430 :
431 16 : !$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 16 : !$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 16 : CALL timestop(handle)
486 :
487 16 : 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 16 : 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 16 : 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 16 : CALL timeset(routineN, handle)
509 :
510 16 : nop = spgr%n_operations
511 16 : nparticle = spgr%nparticle
512 64 : ALLOCATE (x(3*nparticle), xold(3*nparticle))
513 538 : x = 0.0_dp
514 190 : DO ia = 1, nparticle
515 174 : ja = 3*(ia - 1)
516 174 : x(ja + 1) = x(ja + 1) + spgr%pol(1)
517 174 : x(ja + 2) = x(ja + 2) + spgr%pol(2)
518 190 : x(ja + 3) = x(ja + 3) + spgr%pol(3)
519 : END DO
520 538 : xold(:) = x(:)
521 :
522 : nops = 0
523 852 : DO ir = 1, nop
524 12368 : x = 0.d0
525 836 : spgr%lop(ir) = .TRUE.
526 10868 : rot(1:3, 1:3) = spgr%rotations(1:3, 1:3, ir)
527 4680 : DO ia = 1, nparticle
528 3844 : IF (.NOT. spgr%lat(ia)) CYCLE
529 3844 : ja = 3*(ia - 1)
530 15376 : ri(1:3) = xold(ja + 1:ja + 3)
531 3844 : ro(1) = REAL(rot(1, 1), dp)*ri(1) + REAL(rot(2, 1), dp)*ri(2) + REAL(rot(3, 1), dp)*ri(3)
532 3844 : ro(2) = REAL(rot(1, 2), dp)*ri(1) + REAL(rot(2, 2), dp)*ri(2) + REAL(rot(3, 2), dp)*ri(3)
533 3844 : ro(3) = REAL(rot(1, 3), dp)*ri(1) + REAL(rot(2, 3), dp)*ri(2) + REAL(rot(3, 3), dp)*ri(3)
534 16212 : x(ja + 1:ja + 3) = ro(1:3)
535 : END DO
536 4680 : DO ia = 1, nparticle
537 3844 : IF (.NOT. spgr%lat(ia)) CYCLE
538 3844 : ib = spgr%eqatom(ir, ia)
539 3844 : ja = 3*(ia - 1)
540 3844 : jb = 3*(ib - 1)
541 15376 : 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 4680 : .AND. (ABS(ro(3)) < spgr%eps_symmetry))
545 : END DO
546 852 : IF (spgr%lop(ir)) nops = nops + 1
547 : END DO
548 :
549 16 : spgr%n_reduced_operations = nops
550 :
551 16 : DEALLOCATE (x, xold)
552 16 : CALL timestop(handle)
553 :
554 16 : 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 16 : 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 16 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: mask
573 :
574 16 : CALL timeset(routineN, handle)
575 :
576 48 : ALLOCATE (mask(spgr%n_operations))
577 852 : mask = .TRUE.
578 :
579 852 : DO i = 1, spgr%n_operations
580 852 : IF (.NOT. spgr%lop(i)) mask(i) = .FALSE.
581 : END DO
582 :
583 836 : DO i = 1, spgr%n_operations - 1
584 820 : IF (.NOT. mask(i)) CYCLE
585 32594 : DO j = i + 1, spgr%n_operations
586 32330 : IF (.NOT. mask(j)) CYCLE
587 244322 : d(:, :) = spgr%rotations(:, :, j) - spgr%rotations(:, :, i)
588 245142 : IF (SUM(ABS(d)) == 0) mask(j) = .FALSE.
589 : END DO
590 : END DO
591 :
592 16 : spgr%n_operations_subset = 0
593 852 : DO i = 1, spgr%n_operations
594 852 : IF (mask(i)) spgr%n_operations_subset = spgr%n_operations_subset + 1
595 : END DO
596 :
597 48 : ALLOCATE (spgr%rotations_subset(3, 3, spgr%n_operations_subset))
598 :
599 16 : j = 0
600 852 : DO i = 1, spgr%n_operations
601 852 : IF (mask(i)) THEN
602 260 : j = j + 1
603 3380 : spgr%rotations_subset(:, :, j) = spgr%rotations(:, :, i)
604 : END IF
605 : END DO
606 :
607 16 : DEALLOCATE (mask)
608 16 : CALL timestop(handle)
609 :
610 16 : 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 78 : 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 78 : 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 78 : CALL timeset(routineN, handle)
634 :
635 234 : ALLOCATE (cold(SIZE(coord)))
636 3162 : cold(:) = coord(:)
637 :
638 78 : nop = spgr%n_operations
639 78 : nparticle = spgr%nparticle
640 78 : nops = spgr%n_reduced_operations
641 :
642 78 : !$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 78 : DEALLOCATE (cold)
669 78 : CALL timestop(handle)
670 :
671 78 : 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 125 : 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 125 : 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 125 : CALL timeset(routineN, handle)
695 :
696 375 : ALLOCATE (fold(SIZE(force)))
697 5255 : fold(:) = force(:)
698 :
699 125 : nop = spgr%n_operations
700 125 : nparticle = spgr%nparticle
701 125 : nops = spgr%n_reduced_operations
702 :
703 125 : !$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 125 : DEALLOCATE (fold)
725 125 : CALL timestop(handle)
726 :
727 125 : 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 32 : 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 32 : CALL timeset(routineN, handle)
750 :
751 32 : ih1 = inv_3x3(h1)
752 32 : ih2 = inv_3x3(h2)
753 1280 : h2ih1 = MATMUL(h2, ih1)
754 1280 : h1ih2 = MATMUL(h1, ih2)
755 :
756 280 : DO ir = 1, nop
757 3224 : r(:, :) = roti(:, :, ir)
758 9920 : s = MATMUL(h2ih1, r)
759 9920 : r = MATMUL(s, h1ih2)
760 3256 : roto(:, :, ir) = r(:, :)
761 : END DO
762 :
763 32 : CALL timestop(handle)
764 :
765 32 : 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 32 : 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 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: roto
786 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat1, hmat2, r, stin
787 :
788 32 : CALL timeset(routineN, handle)
789 :
790 416 : hmat1 = TRANSPOSE(cell%hmat)
791 :
792 32 : hmat2 = 0d0
793 32 : hmat2(1, 1) = 1.d0
794 32 : hmat2(2, 2) = 1.d0
795 32 : hmat2(3, 3) = 1.d0
796 :
797 32 : nop = spgr%n_operations_subset
798 :
799 96 : ALLOCATE (roto(3, 3, nop))
800 :
801 32 : CALL spgr_change_basis(spgr%rotations_subset, roto, spgr%n_operations_subset, hmat1, hmat2)
802 :
803 32 : stin = stress
804 32 : stress = 0.d0
805 280 : DO ir = 1, nop
806 3224 : r(:, :) = roto(:, :, ir)
807 1024 : DO i = 1, 3
808 3224 : DO j = 1, 3
809 9672 : DO k = 1, 3
810 29016 : DO l = 1, 3
811 26784 : 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 416 : stress = stress/REAL(nop, dp)
818 :
819 32 : DEALLOCATE (roto)
820 :
821 32 : CALL timestop(handle)
822 :
823 32 : 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 16 : SUBROUTINE print_spgr(spgr)
833 :
834 : TYPE(spgr_type), INTENT(IN), POINTER :: spgr
835 :
836 : INTEGER :: i, j
837 :
838 16 : IF (spgr%iunit > 0) THEN
839 8 : WRITE (spgr%iunit, '(/,T2,A,A)') "----------------------------------------", &
840 16 : "---------------------------------------"
841 8 : WRITE (spgr%iunit, "(T2,A,T25,A,T77,A)") "----", "SPACE GROUP SYMMETRY INFORMATION", "----"
842 8 : WRITE (spgr%iunit, '(T2,A,A)') "----------------------------------------", &
843 16 : "---------------------------------------"
844 8 : IF (spgr%symlib) THEN
845 8 : WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| SPACE GROUP NUMBER:", &
846 16 : spgr%space_group_number
847 8 : WRITE (spgr%iunit, '(T2,A,T70,A11)') "SPGR| INTERNATIONAL SYMBOL:", &
848 16 : TRIM(ADJUSTR(spgr%international_symbol))
849 8 : WRITE (spgr%iunit, '(T2,A,T75,A6)') "SPGR| POINT GROUP SYMBOL:", &
850 16 : TRIM(ADJUSTR(spgr%pointgroup_symbol))
851 8 : WRITE (spgr%iunit, '(T2,A,T74,A7)') "SPGR| SCHOENFLIES SYMBOL:", &
852 16 : TRIM(ADJUSTR(spgr%schoenflies))
853 8 : WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF SYMMETRY OPERATIONS:", &
854 16 : spgr%n_operations
855 8 : WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF UNIQUE ROTATIONS:", &
856 16 : spgr%n_operations_subset
857 8 : WRITE (spgr%iunit, '(T2,A,T73,I8)') "SPGR| NUMBER OF REDUCED SYMMETRY OPERATIONS:", &
858 16 : spgr%n_reduced_operations
859 8 : WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF PARTICLES AND SYMMETRIZED PARTICLES:", &
860 16 : spgr%nparticle, spgr%nparticle_sym
861 8 : WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF ATOMS AND SYMMETRIZED ATOMS:", &
862 16 : spgr%n_atom, spgr%n_atom_sym
863 8 : WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF CORES AND SYMMETRIZED CORES:", &
864 16 : spgr%n_core, spgr%n_core_sym
865 8 : WRITE (spgr%iunit, '(T2,A,T65,I8,I8)') "SPGR| NUMBER OF SHELLS AND SYMMETRIZED SHELLS:", &
866 16 : spgr%n_shell, spgr%n_shell_sym
867 8 : 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 16 : 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 32 : 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 416 : stress_tensor(:, :) = stress(:, :)*pascal*1.0E-9_dp
913 :
914 32 : IF (spgr%iunit > 0) THEN
915 : WRITE (UNIT=spgr%iunit, FMT='(/,T2,A)') &
916 17 : 'SPGR STRESS| Symmetrized stress tensor [GPa]'
917 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T19,3(19X,A1))') &
918 17 : 'SPGR STRESS|', 'x', 'y', 'z'
919 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
920 17 : 'SPGR STRESS| x', stress_tensor(1, 1:3)
921 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
922 17 : 'SPGR STRESS| y', stress_tensor(2, 1:3)
923 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
924 17 : 'SPGR STRESS| z', stress_tensor(3, 1:3)
925 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T66,ES20.11)') &
926 17 : 'SPGR STRESS| 1/3 Trace', (stress_tensor(1, 1) + &
927 : stress_tensor(2, 2) + &
928 34 : stress_tensor(3, 3))/3.0_dp
929 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T66,ES20.11)') &
930 17 : 'SPGR STRESS| Determinant', det_3x3(stress_tensor(1:3, 1), &
931 : stress_tensor(1:3, 2), &
932 34 : stress_tensor(1:3, 3))
933 17 : eigval(:) = 0.0_dp
934 17 : eigvec(:, :) = 0.0_dp
935 17 : CALL jacobi(stress_tensor, eigval, eigvec)
936 : WRITE (UNIT=spgr%iunit, FMT='(/,T2,A)') &
937 17 : 'SPGR STRESS| Eigenvectors and eigenvalues of the symmetrized stress tensor [GPa]'
938 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T19,3(1X,I19))') &
939 17 : 'SPGR STRESS|', 1, 2, 3
940 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,ES19.11))') &
941 17 : 'SPGR STRESS| Eigenvalues', eigval(1:3)
942 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,F19.12))') &
943 17 : 'SPGR STRESS| x', eigvec(1, 1:3)
944 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,F19.12))') &
945 17 : 'SPGR STRESS| y', eigvec(2, 1:3)
946 : WRITE (UNIT=spgr%iunit, FMT='(T2,A,T26,3(1X,F19.12))') &
947 17 : 'SPGR STRESS| z', eigvec(3, 1:3)
948 : END IF
949 :
950 32 : END SUBROUTINE spgr_write_stress_tensor
951 :
952 : END MODULE space_groups
|
