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