Branch data Line data Source code
1 : : !!****m* ABINIT/m_geometry
2 : : !! NAME
3 : : !! m_geometry
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains basic tools to operate on vectors expressed in reduced coordinates.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2008-2026 ABINIT group (MG, MT, FJ, TRangel, DCA, XG, AHR, DJA, DRH)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : MODULE m_geometry
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_errors
27 : : use m_atomdata
28 : : use m_sort
29 : :
30 : : use m_io_tools, only : open_file
31 : : use m_numeric_tools, only : uniformrandom, isinteger, set2unit
32 : : use m_matrix, only : mati3inv, mati3det, matr3inv
33 : : use m_symtk, only : symdet
34 : : use m_hide_lapack, only : matr3eigval
35 : : use m_pptools, only : prmat
36 : : use m_numeric_tools, only : wrap2_pmhalf
37 : : use m_hide_lapack, only : matrginv
38 : :
39 : : implicit none
40 : :
41 : : private
42 : :
43 : : public :: normv ! Norm of vector(s) in reduced coordinates either in real or reciprocal space.
44 : : public :: vdotw ! Scalar product between two reduced vectors either in real or reciprocal space.
45 : : public :: acrossb ! Cross product of two 3-vectors.
46 : : public :: wigner_seitz ! Find the grid of points falling inside the Wigner-Seitz cell.
47 : : public :: phdispl_cart2red ! Calculate the displacement vectors for all branches in reduced coordinates.
48 : : public :: phdispl_cart2red_nmodes ! Calculate the displacement vectors for nmodes in reduced coordinates.
49 : : public :: getspinrot ! Compute the components of the spinor rotation matrix
50 : : public :: spinrot_cmat ! Construct 2x2 complex matrix representing rotation operator in spin-space.
51 : : public :: rotmat ! Finds the rotation matrix.
52 : : public :: fixsym ! Check that iatfix does not break symmetry.
53 : : public :: det3r ! Compute determinant of a 3x3 real matrix
54 : : public :: metric ! Compute metric matrices.
55 : : public :: mkradim ! Make rprim and acell from rprimd
56 : : public :: mkrdim ! Make rprimd from acell from rprim
57 : : public :: chkrprimd ! Test if {rprim,acell,rprimd} are consistent
58 : : public :: chkdilatmx ! check if dilatation of unit cell is consistent with initial G-sphere
59 : : public :: xcart2xred ! From cart coords to reduced
60 : : public :: xred2xcart ! From reduced coords to cart.
61 : : public :: gred2fcart ! Convert reduced gradients into cartesian forces
62 : : public :: fcart2gred ! Convert cartesian forces into reduced gradients
63 : : public :: bonds_lgth_angles ! Write GEO file
64 : : public :: randomcellpos ! Creates unit cell with random atomic positions.
65 : : public :: ioniondist ! Compute ion-ion distances
66 : : public :: dist2 ! Calculates the distance of v1 and v2 in a crystal by repeating the unit cell
67 : : public :: shellstruct ! Calculates shell structure (multiplicities, radii)
68 : : public :: remove_inversion ! Remove the inversion symmetry and improper rotations
69 : : public :: reduce2primitive ! Find real space primitive vectors from non-primitive ones and set of translations
70 : : public :: symredcart ! Convert a symmetry operation from reduced coordinates (integers) to cart coords (reals)
71 : : public :: strainsym ! Symmetrize the strain tensor.
72 : : public :: stresssym ! Symmetrize the stress tensor.
73 : : public :: stress_voigt_to_mat! Build 3x3 symmetric stress tensor from stress vector in Voigt notation.
74 : : public :: strconv ! Convert from symmetric storage mode in reduced coords to cart coords.
75 : : public :: littlegroup_pert ! Determines the set of symmetries that leaves a perturbation invariant.
76 : : public :: irreducible_set_pert ! Determines a set of perturbations that form a basis
77 : : public :: wedge_basis ! compute rprimd x gprimd vectors needed for generalized cross product
78 : : public :: wedge_product ! compute wedge product given wedge basis
79 : : public :: d3lwsym
80 : : public :: sylwtens ! Determines the set of irreductible elements of the spatial-dispersion tensors
81 : : public :: cart2spinaxis ! Compute the rotation matrix from cartesian to spinaxis coordinate
82 : : public :: vcart2ylm ! Convert Cartesian vector to spherical coordinates for Y_lm
83 : :
84 : : interface normv
85 : : module procedure normv_rdp_vector
86 : : module procedure normv_int_vector
87 : : !module procedure normv_int_vector_array ! WARNING for the time being, do not use these 2 procedures,
88 : : !module procedure normv_rdp_vector_array ! sunstudio12 is not able to resolve which sub should be called.
89 : : end interface normv
90 : :
91 : : interface vdotw
92 : : module procedure vdotw_rr_vector
93 : : module procedure vdotw_rc_vector
94 : : end interface vdotw
95 : :
96 : : CONTAINS !===========================================================
97 : : !!***
98 : :
99 : : !!****f* m_geometry/normv_rdp_vector
100 : : !! NAME
101 : : !! normv_rdp_vector
102 : : !!
103 : : !! FUNCTION
104 : : !! Compute the norm of a vector expressed in reduced coordinates using the metric met.
105 : : !! The result is multiplied by 2pi in case of a vector in reciprocal space
106 : : !! to take into account the correct normalisation of the reciprocal lattice vectors
107 : : !!
108 : : !! INPUTS
109 : : !! xv(3)=Vector in reduced coordinates
110 : : !! met(3,3)=Metric tensor
111 : : !! space=Character defining whether we are working in real (r|R) or reciprocal space (g|G)
112 : : !!
113 : : !! OUTPUT
114 : : !! normv_rdp_vector=norm of xv
115 : : !!
116 : : !! NOTES
117 : : !! The routine is able to deal both with a single vector as well as arrays of vectors.
118 : : !! Versions for integer and real vectors are provided.
119 : : !!
120 : : !! SOURCE
121 : :
122 : 1762457911 : function normv_rdp_vector(xv,met,space) result(res)
123 : :
124 : : !Arguments ------------------------------------
125 : : !scalars
126 : : real(dp) :: res
127 : : character(len=1),intent(in) :: space
128 : : !arrays
129 : : real(dp),intent(in) :: met(3,3),xv(3)
130 : : ! *************************************************************************
131 : :
132 : : res = (xv(1)*met(1,1)*xv(1) + xv(2)*met(2,2)*xv(2) + xv(3)*met(3,3)*xv(3) &
133 : 1762457911 : & +two*(xv(1)*met(1,2)*xv(2) + xv(1)*met(1,3)*xv(3) + xv(2)*met(2,3)*xv(3)) )
134 : :
135 : 2782949 : select case (space)
136 : : case ('r','R')
137 : 2782949 : res=SQRT(res)
138 : : case ('g','G')
139 : 1759674962 : res=two_pi*SQRT(res)
140 : : case default
141 [ + + - ]: 1762457911 : ABI_BUG('Wrong value for space')
142 : : end select
143 : :
144 : 1762457911 : end function normv_rdp_vector
145 : : !!***
146 : :
147 : : !----------------------------------------------------------------------
148 : :
149 : : !!****f* m_geometry/normv_int_vector
150 : : !! NAME
151 : : !! normv_int_vector
152 : : !!
153 : : !! FUNCTION
154 : : !! Returns the norm of an integer 3D vector expressed in reduced coordinates.
155 : : !! either in real or reciprocal space. In the later case the factor 2pi has
156 : : !! to be included, due to the conventions used in abinit to define the reciprocal lattice.
157 : : !!
158 : : !! INPUTS
159 : : !!
160 : : !! OUTPUT
161 : : !!
162 : : !! SOURCE
163 : :
164 : 1118484 : function normv_int_vector(xv, met, space) result(res)
165 : :
166 : : !Arguments ------------------------------------
167 : : !scalars
168 : : real(dp) :: res
169 : : character(len=1),intent(in) :: space
170 : : !arrays
171 : : real(dp),intent(in) :: met(3,3)
172 : : integer,intent(in) :: xv(3)
173 : : ! *************************************************************************
174 : :
175 : : res = ( xv(1)*met(1,1)*xv(1) + xv(2)*met(2,2)*xv(2) + xv(3)*met(3,3)*xv(3) &
176 : 1118484 : & +two*( xv(1)*met(1,2)*xv(2) + xv(1)*met(1,3)*xv(3) + xv(2)*met(2,3)*xv(3)) )
177 : :
178 : 0 : select case (space)
179 : : case ('r','R')
180 : 0 : res=SQRT(res)
181 : : case ('g','G')
182 : 1118484 : res=two_pi*SQRT(res)
183 : : case default
184 [ - + - ]: 1118484 : ABI_BUG('Wrong value for space')
185 : : end select
186 : :
187 : 1118484 : end function normv_int_vector
188 : : !!***
189 : :
190 : : !----------------------------------------------------------------------
191 : :
192 : : !!****f* m_geometry/normv_int_vector_array
193 : : !! NAME
194 : : !! normv_int_vector_array
195 : : !!
196 : : !! FUNCTION
197 : : !! Returns the norm of an array of integer 3D vectors expressed in reduced coordinates.
198 : : !! either in real or reciprocal space. In the later case the factor 2pi has
199 : : !! to be included, due to the conventions used in abinit to define the reciprocal lattice.
200 : : !!
201 : : !! INPUTS
202 : : !!
203 : : !! OUTPUT
204 : : !!
205 : : !! SOURCE
206 : :
207 : : function normv_int_vector_array(xv,met,space) result(res)
208 : :
209 : : !Arguments ------------------------------------
210 : : !scalars
211 : : character(len=1),intent(in) :: space
212 : : !arrays
213 : : real(dp),intent(in) :: met(3,3)
214 : : integer,intent(in) :: xv(:,:)
215 : : !this awful trick is needed to avoid problems with abilint
216 : : real(dp) :: res(SIZE(xv(1,:)))
217 : : ! *************************************************************************
218 : :
219 : : res(:) = ( xv(1,:)*met(1,1)*xv(1,:) + xv(2,:)*met(2,2)*xv(2,:) + xv(3,:)*met(3,3)*xv(3,:) &
220 : : +two*(xv(1,:)*met(1,2)*xv(2,:) + xv(1,:)*met(1,3)*xv(3,:) + xv(2,:)*met(2,3)*xv(3,:)) )
221 : :
222 : : select case (space)
223 : : case ('r','R')
224 : : res(:)=SQRT(res(:))
225 : : case ('g','G')
226 : : res(:)=two_pi*SQRT(res(:))
227 : : case default
228 : : ABI_BUG('Wrong value for space')
229 : : end select
230 : :
231 : : end function normv_int_vector_array
232 : : !!***
233 : :
234 : : !----------------------------------------------------------------------
235 : :
236 : : !!****f* m_geometry/normv_rdp_vector_array
237 : : !! NAME
238 : : !! normv_rdp_vector_array
239 : : !!
240 : : !! FUNCTION
241 : : !! Returns the norm of an array of real 3D vectors expressed in reduced coordinates.
242 : : !! either in real or reciprocal space. In the later case the factor 2pi has
243 : : !! to be included, due to the conventions used in abinit to define the reciprocal lattice.
244 : : !!
245 : : !! INPUTS
246 : : !!
247 : : !! OUTPUT
248 : : !!
249 : : !! SOURCE
250 : :
251 : : function normv_rdp_vector_array(xv,met,space) result(res)
252 : :
253 : : !Arguments ------------------------------------
254 : : !scalars
255 : : character(len=1),intent(in) :: space
256 : : !arrays
257 : : real(dp),intent(in) :: met(3,3)
258 : : real(dp),intent(in) :: xv(:,:)
259 : : !this awful trick is needed to avoid problems with abilint
260 : : real(dp) :: res(SIZE(xv(1,:)))
261 : : !real(dp) :: res(SIZE(xv,DIM=2))
262 : : ! *************************************************************************
263 : :
264 : : res(:) = ( xv(1,:)*met(1,1)*xv(1,:) + xv(2,:)*met(2,2)*xv(2,:) + xv(3,:)*met(3,3)*xv(3,:) &
265 : : & +two*(xv(1,:)*met(1,2)*xv(2,:) + xv(1,:)*met(1,3)*xv(3,:) + xv(2,:)*met(2,3)*xv(3,:)) )
266 : :
267 : : select case (space)
268 : : case ('r','R')
269 : : res(:)=SQRT(res(:))
270 : : case ('g','G')
271 : : res(:)=two_pi*SQRT(res)
272 : : case default
273 : : ABI_BUG('Wrong value for space')
274 : : end select
275 : :
276 : : end function normv_rdp_vector_array
277 : : !!***
278 : :
279 : : !----------------------------------------------------------------------
280 : :
281 : : !!****f* m_geometry/vdotw_rr_vector
282 : : !! NAME
283 : : !! vdotw_rr_vector
284 : : !!
285 : : !! FUNCTION
286 : : !! Compute the scalar product between two vectors expressed in reduced coordinates
287 : : !! The result is multiplied by (2pi)**2 in case of vectors in reciprocal space
288 : : !! to take into account the correct normalisation of the reciprocal lattice vectors
289 : : !!
290 : : !! INPUTS
291 : : !! xv(3),xw(3)=Vectors in reduced coordinates
292 : : !! met(3,3)=Metric tensor
293 : : !! space=Character defining whether we are working in real (r) or reciprocal space (g)
294 : : !!
295 : : !! OUTPUT
296 : : !! res=scalar product of xv and xw
297 : : !!
298 : : !! SOURCE
299 : :
300 : 0 : real(dp) function vdotw_rr_vector(xv,xw,met,space) result(res)
301 : :
302 : : !Arguments ------------------------------------
303 : : character(len=1),intent(in) :: space
304 : : !arrays
305 : : real(dp),intent(in) :: met(3,3),xv(3),xw(3)
306 : : ! *************************************************************************
307 : :
308 : : res = ( met(1,1)* xv(1)*xw(1) &
309 : : +met(2,2)* xv(2)*xw(2) &
310 : : +met(3,3)* xv(3)*xw(3) &
311 : : +met(1,2)*(xv(1)*xw(2) + xv(2)*xw(1)) &
312 : : +met(1,3)*(xv(1)*xw(3) + xv(3)*xw(1)) &
313 : 0 : +met(2,3)*(xv(2)*xw(3) + xv(3)*xw(2)) )
314 : :
315 : : select case (space)
316 : : case ('r', 'R')
317 : 0 : return
318 : : case ('g', 'G')
319 : 0 : res= res * (two_pi**2)
320 : : case default
321 [ # # # ]: 0 : ABI_BUG('Wrong value for space')
322 : : end select
323 : :
324 : : end function vdotw_rr_vector
325 : : !!***
326 : :
327 : : !----------------------------------------------------------------------
328 : :
329 : : !!****f* m_geometry/vdotw_rc_vector
330 : : !! NAME
331 : : !! vdotw_rc_vector
332 : : !!
333 : : !! FUNCTION
334 : : !! Compute the scalar product between two vectors expressed in reduced coordinates
335 : : !! First vector is real, the second one is complex.
336 : : !! The result is multiplied by (2pi)**2 in case of vectors in reciprocal space
337 : : !! to take into account the correct normalisation of the reciprocal lattice vectors
338 : : !!
339 : : !! INPUTS
340 : : !! xv(3),xw(3)=Vectors in reduced coordinates
341 : : !! met(3,3)=Metric tensor
342 : : !! space=Character defining whether we are working in real (r) or reciprocal space (g)
343 : : !!
344 : : !! OUTPUT
345 : : !! res=complex scalar product of xv and xw
346 : : !!
347 : : !! SOURCE
348 : :
349 : 78278 : complex(dp) function vdotw_rc_vector(xv, xw, met, space) result(res)
350 : :
351 : : !Arguments ------------------------------------
352 : : character(len=1),intent(in) :: space
353 : : !arrays
354 : : real(dp),intent(in) :: met(3,3),xv(3)
355 : : complex(dp),intent(in) :: xw(3)
356 : : ! *************************************************************************
357 : :
358 : : res = ( met(1,1)* xv(1)*xw(1) &
359 : : +met(2,2)* xv(2)*xw(2) &
360 : : +met(3,3)* xv(3)*xw(3) &
361 : : +met(1,2)*(xv(1)*xw(2) + xv(2)*xw(1)) &
362 : : +met(1,3)*(xv(1)*xw(3) + xv(3)*xw(1)) &
363 : 78278 : +met(2,3)*(xv(2)*xw(3) + xv(3)*xw(2)) )
364 : :
365 : : select case (space)
366 : : case ('r', 'R')
367 : 78278 : return
368 : : case ('g', 'G')
369 : 78278 : res= res * (two_pi**2)
370 : : case default
371 [ - + - ]: 78278 : ABI_BUG('Wrong value for space')
372 : : end select
373 : :
374 : : end function vdotw_rc_vector
375 : : !!***
376 : :
377 : : !----------------------------------------------------------------------
378 : :
379 : : !!****f* m_geometry/acrossb
380 : : !! NAME
381 : : !! acrossb
382 : : !!
383 : : !! FUNCTION
384 : : !! Calculates the cross product of two 3-vectors
385 : : !!
386 : : !! INPUTS
387 : : !! a(3): real(dp) vector
388 : : !! b(3): real(dp) vector
389 : : !!
390 : : !! OUTPUT
391 : : !! c(3): real(dp) vector = a X b
392 : : !!
393 : : !! SOURCE
394 : :
395 : 5484182 : subroutine acrossb(a,b,c)
396 : :
397 : : !Arguments ---------------------------------------------
398 : : !arrays
399 : : real(dp),intent(in) :: a(3),b(3)
400 : : real(dp),intent(out) :: c(3)
401 : : ! *********************************************************************
402 : :
403 : 5484182 : c(1) = a(2)*b(3) - a(3)*b(2)
404 : 5484182 : c(2) = -a(1)*b(3) + a(3)*b(1)
405 : 5484182 : c(3) = a(1)*b(2) - b(1)*a(2)
406 : :
407 : 5484182 : end subroutine acrossb
408 : : !!***
409 : :
410 : : !!****f* m_geometry/wedge_basis
411 : : !! NAME
412 : : !! wedge_basis
413 : : !!
414 : : !! FUNCTION
415 : : !! Calculates the basis vectors a ^ a* for a in rprimd and
416 : : !! a* in gprimd, needed for some generalized cross products
417 : : !!
418 : : !! INPUTS
419 : : !! rprimd(3,3) : real(dp) matrix
420 : : !! gprimd(3,3) : real(dp) matrix
421 : : !! normalize,optional : whether to normalize the output vectors
422 : : !!
423 : : !! OUTPUT
424 : : !! wedge(3,3,3) : 9 basis vectors of rprimd ^ gprimd
425 : : !!
426 : : !! SOURCE
427 : :
428 : 110 : subroutine wedge_basis(gprimd,rprimd,wedge,normalize)
429 : :
430 : : !Arguments ---------------------------------------------
431 : : ! scalars
432 : : logical,optional,intent(in) :: normalize
433 : : !arrays
434 : : real(dp),intent(in) :: gprimd(3,3),rprimd(3,3)
435 : : real(dp),intent(out) :: wedge(3,3,3)
436 : :
437 : : ! local
438 : : !scalars
439 : : integer :: igprimd, irprimd
440 : : real(dp) :: nfac
441 : : logical :: nvec
442 : : ! *********************************************************************
443 : :
444 [ - + ]: 110 : if(present(normalize)) then
445 : 0 : nvec = normalize
446 : : else
447 : : nvec = .FALSE.
448 : : end if
449 : :
450 [ + + ]: 440 : do irprimd = 1, 3
451 [ + + ]: 1430 : do igprimd = 1, 3
452 : 990 : wedge(1,irprimd,igprimd) = rprimd(2,irprimd)*gprimd(3,igprimd) - rprimd(3,irprimd)*gprimd(2,igprimd)
453 : 990 : wedge(2,irprimd,igprimd) = rprimd(3,irprimd)*gprimd(1,igprimd) - rprimd(1,irprimd)*gprimd(3,igprimd)
454 : 1320 : wedge(3,irprimd,igprimd) = rprimd(1,irprimd)*gprimd(2,igprimd) - rprimd(2,irprimd)*gprimd(1,igprimd)
455 : : end do
456 : : end do
457 : :
458 [ - + ]: 110 : if (nvec) then
459 [ # # ]: 0 : do irprimd = 1, 3
460 [ # # ]: 0 : do igprimd = 1, 3
461 [ # # # # : 0 : if(any(abs(wedge(1:3,irprimd,igprimd)).GT.tol8)) then
# # ]
462 [ # # # # : 0 : nfac = NORM2(wedge(1:3,irprimd,igprimd))
# # ]
463 [ # # ]: 0 : wedge(1:3,irprimd,igprimd) = wedge(1:3,irprimd,igprimd)/nfac
464 : : end if
465 : : end do
466 : : end do
467 : : end if
468 : :
469 : 110 : end subroutine wedge_basis
470 : : !!***
471 : :
472 : : !!****f* m_geometry/wedge_product
473 : : !! NAME
474 : : !! wedge_product
475 : : !!
476 : : !! FUNCTION
477 : : !! Calculates the wedge product u^w, given the wedge product basis a^b
478 : : !! typically u=(u1 a + u2 b + u3 c) and w = (w1 a* + w2 b* + w3 c*)
479 : : !!
480 : : !! INPUTS
481 : : !! u(3) :: real(dp) input vector
482 : : !! v(3) :: real(dp) input vector
483 : : !! wedgebasis(3,3,3) :: real(dp) input matrix
484 : : !!
485 : : !! OUTPUT
486 : : !! produv(3) :: real(dp) output vector
487 : : !!
488 : : !! SOURCE
489 : :
490 : 14464612 : subroutine wedge_product(produv,u,v,wedgebasis)
491 : :
492 : : !Arguments ---------------------------------------------
493 : : !arrays
494 : : real(dp),intent(in) :: u(3),v(3),wedgebasis(3,3,3)
495 : : real(dp),intent(out) :: produv(3)
496 : :
497 : : ! local
498 : : integer :: igprimd, irprimd
499 : : ! *********************************************************************
500 : :
501 : 14464612 : produv(:) = zero
502 [ + + ]: 57858448 : do irprimd = 1, 3
503 [ + + ]: 188039956 : do igprimd = 1, 3
504 [ + + ]: 564119868 : produv(1:3) = produv(1:3) + u(irprimd)*v(igprimd)*wedgebasis(1:3,irprimd,igprimd)
505 : : end do
506 : : end do
507 : :
508 : 14464612 : end subroutine wedge_product
509 : : !!***
510 : :
511 : : !!****f* m_geometry/wigner_seitz
512 : : !! NAME
513 : : !! wigner_seitz
514 : : !!
515 : : !! FUNCTION
516 : : !! Calculates a grid of points that falls inside of (and eventually on the surface of)
517 : : !! the Wigner-Seitz supercell centered on the origin of the B lattice with primitive
518 : : !! translations nmonkh(1)*a_1+nmonkh(2)*a_2+nmonkh(3)*a_3.
519 : : !! Subroutine taken from the Wannier90 code.
520 : : !! Modified by MG to fulfil abinit coding rules. API slightly changed wrt the wannier90 version.
521 : : !!
522 : : !! COPYRIGHT
523 : : !! Copyright (C) 2007 Jonathan Yates, Arash Mostofi,
524 : : !! Young-Su Lee, Nicola Marzari, Ivo Souza, David Vanderbilt.
525 : : !! This file is distributed under the terms of the
526 : : !! GNU General Public License, see ~abinit/COPYING
527 : : !! or http://www.gnu.org/copyleft/gpl.txt .
528 : : !!
529 : : !! INPUTS
530 : : !! center(3)=The Wigner-Seitz cell is centered on this point in reduced coordinates.
531 : : !! rmet(3,3)=Real space metric ($\textrm{bohr}^{2}$).
532 : : !! kptrlatt(3)=Values defining the supercell.
533 : : !! prtvol=If different from 0 print out the points falling inside the W-S cell and the corresponding weights.
534 : : !! lmax(3)=see Notes below.
535 : : !!
536 : : !! OUTPUT
537 : : !! npts=number of points falling inside the Wigner-Seitz cell
538 : : !! irvec(3,npts)=Reduced coordinated of the points inside the W-S cell (sorted by length)
539 : : !! ndegen(npts)=Weights associated to each point.
540 : : !! rmods(npts)=length of the irvec
541 : : !!
542 : : !! SIDE EFFECTS
543 : : !! irvec, ndegen and rmods are allocated with the correct
544 : : !! size inside the routine and returned to the caller.
545 : : !!
546 : : !! NOTES
547 : : !! The Wannier functions live in a supercell of the real space unit cell.
548 : : !! This supercell is mp_grid unit cells long in each direction
549 : : !! The algorithm loops over grid points r on a unit cell that is 8 times larger than this
550 : : !! primitive supercell.
551 : : !! One of these points is in the W-S cell if it is closer to center(:)
552 : : !! than any of the other points R where R are the translation vectors of the supercell.
553 : : !! In the end npts contains the total number of grid points that have been found in the Wigner-Seitz cell
554 : : !! The number of lattice vectors R along each direction of the supercell is defined by lmax.
555 : : !!
556 : : !! SOURCE
557 : :
558 : 15 : subroutine wigner_seitz(center, lmax, kptrlatt, rmet, npts, irvec, ndegen, rmods, prtvol)
559 : :
560 : : !Arguments ------------------------------------
561 : : !scalars
562 : : integer,optional,intent(in) :: prtvol
563 : : integer,intent(out) :: npts
564 : : !arrays
565 : : integer,intent(in) :: kptrlatt(3,3),lmax(3)
566 : : integer,allocatable,intent(out) :: irvec(:,:),ndegen(:)
567 : : real(dp),intent(in) :: center(3),rmet(3,3)
568 : : real(dp),allocatable,intent(out) :: rmods(:)
569 : :
570 : : !Local variables-------------------------------
571 : : !scalars
572 : : integer :: in1,in2,in3,l1,l2,l3,ii,icount,n1,n2,n3
573 : : integer :: l0,l1_max,l2_max,l3_max,nl,verbose,mm1,mm2,mm3,ir
574 : : real(dp),parameter :: TOL_DIST=tol7
575 : : real(dp) :: tot,dist_min
576 : : character(len=500) :: msg
577 : : !arrays
578 : 15 : integer,allocatable :: iperm(:), swap2(:,:), swap1(:)
579 : : real(dp),parameter :: gammak(3) = zero
580 : : real(dp) :: diff(3)
581 : 15 : real(dp),allocatable :: dist(:)
582 : : ! *************************************************************************
583 : :
584 [ + - ]: 15 : verbose = 0; if (present(prtvol)) verbose = prtvol
585 : :
586 : : if (kptrlatt(1,2) /= 0 .or. kptrlatt(2,1) /= 0 .or. &
587 : : kptrlatt(1,3) /= 0 .or. kptrlatt(3,1) /= 0 .or. &
588 [ + - + - : 15 : kptrlatt(2,3) /= 0 .or. kptrlatt(3,2) /= 0 ) then
+ - + - +
- - + ]
589 : 0 : ABI_ERROR('Off-diagonal elements of kptrlatt must be zero')
590 : : end if
591 : :
592 : 15 : n1 = kptrlatt(1,1); n2 = kptrlatt(2,2); n3 = kptrlatt(3,3)
593 : 15 : l1_max = lmax(1); l2_max = lmax(2); l3_max = lmax(3)
594 : :
595 : 15 : nl = (2*l1_max+1)*(2*l2_max+1)*(2*l3_max+1)
596 : 15 : l0 = 1+l1_max*(1+(2*l2_max+1)**2+(2*l3_max+1)) ! Index of the origin.
597 [ - + + - : 45 : ABI_MALLOC(dist, (nl))
+ - + - -
+ - + ]
598 : :
599 : : ! Allocate with maximum size
600 : 15 : mm1 = 2 * n1 + 1
601 : 15 : mm2 = 2 * n2 + 1
602 : 15 : mm3 = 2 * n3 + 1
603 [ - + + - : 45 : ABI_MALLOC(irvec, (3, mm1*mm2*mm3))
+ - + - -
+ - + -
+ ]
604 [ - + + - : 45 : ABI_MALLOC(ndegen, (mm1*mm2*mm3))
+ - + - -
+ - + -
+ ]
605 : :
606 : 15 : npts = 0
607 [ + + ]: 118 : do in1=-n1,n1
608 [ + + ]: 917 : do in2=-n2,n2
609 [ + + ]: 7869 : do in3=-n3,n3
610 : :
611 : : ! Loop over the nl points R. R=0 corresponds to l1=l2=l3=1, or icount=l0
612 : 6967 : icount = 0
613 [ + + ]: 41802 : do l1=-l1_max,l1_max
614 [ + + ]: 215977 : do l2=-l2_max,l2_max
615 [ + + ]: 1079885 : do l3=-l3_max,l3_max
616 : : ! Calculate |r - R -r0|^2.
617 : 870875 : diff(1) = in1 - l1 * n1 - center(1)
618 : 870875 : diff(2) = in2 - l2 * n2 - center(2)
619 : 870875 : diff(3) = in3 - l3 * n3 - center(3)
620 : 870875 : icount = icount+1
621 [ + + + + : 14108175 : dist(icount) = DOT_PRODUCT(diff, MATMUL(rmet, diff))
+ + ]
622 : : end do
623 : : end do
624 : : end do
625 : :
626 [ + - - + : 884809 : dist_min = MINVAL(dist)
- - + + ]
627 : :
628 [ + + ]: 7766 : if (ABS(dist(l0) - dist_min) < TOL_DIST) then
629 : 903 : npts = npts + 1
630 : 903 : ndegen (npts) = 0
631 [ + + ]: 113778 : do ii=1,nl
632 [ + + ]: 113778 : if (ABS(dist(ii) - dist_min) < TOL_DIST) ndegen(npts) = ndegen(npts) + 1
633 : : end do
634 [ + + ]: 3612 : irvec(:, npts) = [in1, in2, in3]
635 : : end if
636 : : end do !in3
637 : : end do !in2
638 : : end do !in1
639 : :
640 [ + + ]: 15 : if (verbose >= 1) then
641 : 3 : write(msg,'(a,i0)')' lattice points in Wigner-Seitz supercell: ',npts
642 : 3 : call wrtout(std_out, msg)
643 [ + + ]: 262 : do ii=1,npts
644 : 259 : write(msg,'(a,3(i3,1x),a,i0)')' vector: ', irvec(:,ii),', degeneracy: ', ndegen(ii)
645 : 262 : call wrtout(std_out, msg)
646 : : end do
647 : : end if
648 : :
649 : : ! Check the "sum rule"
650 : 15 : tot = zero
651 [ + + ]: 918 : do ii=1,npts
652 : 918 : tot = tot + one/ndegen(ii)
653 : : end do
654 [ - + ]: 15 : if (ABS(tot-(n1*n2*n3)) > tol8) then
655 : 0 : write(msg,'(a,es16.8,a,i0)')'Something wrong in the generation of the WS mesh: tot: ',tot,' /= n1*n2*n3: ',n1*n2*n3
656 : 0 : ABI_ERROR(msg)
657 : : end if
658 : :
659 [ - + ]: 15 : ABI_FREE(dist)
660 : :
661 : : ! Reallocate ndegen and irvec with correct size and sort by norm
662 [ - + + - : 45 : ABI_MALLOC(swap2, (3, npts))
+ - + - -
+ - + ]
663 [ + + + + ]: 3627 : swap2(:,:) = irvec(1:3, 1:npts)
664 [ - + ]: 15 : ABI_FREE(irvec)
665 : :
666 : 15 : call sort_gvecs(npts, gammak, rmet, swap2, out_gvec=irvec, iperm=iperm)
667 [ - + ]: 15 : ABI_FREE(swap2)
668 : :
669 [ - + + - : 45 : ABI_MALLOC(swap1, (npts))
+ - + - -
+ - + ]
670 [ + + ]: 918 : swap1(:) = ndegen(1:npts)
671 [ + - - + : 30 : ABI_REMALLOC(ndegen, (npts))
+ - - + -
+ - + ]
672 : :
673 [ + + ]: 918 : do ir=1,npts
674 : 918 : ndegen(ir) = swap1(iperm(ir))
675 : : end do
676 : :
677 [ - + + - : 45 : ABI_MALLOC(rmods, (npts))
+ - + - -
+ - + -
+ ]
678 [ + + ]: 918 : do ir=1,npts
679 [ + + + + : 17172 : rmods(ir) = sqrt(dot_product(irvec(:,ir), matmul(rmet, irvec(:,ir))))
+ + + - +
+ ]
680 : : end do
681 : :
682 [ - + ]: 15 : ABI_FREE(swap1)
683 [ - + ]: 15 : ABI_FREE(iperm)
684 : :
685 : 15 : end subroutine wigner_seitz
686 : : !!***
687 : :
688 : : !----------------------------------------------------------------------
689 : :
690 : : !!****f* m_geometry/phdispl_cart2red
691 : : !! NAME
692 : : !! phdispl_cart2red
693 : : !!
694 : : !! FUNCTION
695 : : !! Calculates the displacement vectors for all branches in reduced coordinates.
696 : : !! $ displ_red = displ_cart \cdot gprimd $ for each phonon branch.
697 : : !!
698 : : !! INPUTS
699 : : !! natom=Number of atoms.
700 : : !! gprimd(3,3)=Dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
701 : : !! displ_cart(2,3*natom,3*natom)=Phonon displacement in Cartesian coordinates.
702 : : !!
703 : : !! OUTPUT
704 : : !! displ_red(2,3*natom,3*natom)=Phonon displacement in reduded coordinates.
705 : : !!
706 : : !! SOURCE
707 : :
708 : 12085 : pure subroutine phdispl_cart2red(natom, gprimd, displ_cart, displ_red)
709 : :
710 : : !Arguments ------------------------------------
711 : : !scalars
712 : : integer,intent(in) :: natom
713 : : !arrays
714 : : real(dp),intent(in) :: gprimd(3,3)
715 : : real(dp),intent(in) :: displ_cart(2,3*natom,3*natom)
716 : : real(dp),intent(out) :: displ_red(2,3*natom,3*natom)
717 : : ! *************************************************************************
718 : :
719 : 12085 : call phdispl_cart2red_nmodes(natom, 3*natom, gprimd, displ_cart, displ_red)
720 : :
721 : 12085 : end subroutine phdispl_cart2red
722 : : !!***
723 : :
724 : : !!****f* m_geometry/phdispl_cart2red_nmodes
725 : : !! NAME
726 : : !! phdispl_cart2red_nmodes
727 : : !!
728 : : !! FUNCTION
729 : : !! Similar to phdispl_cart2red but operates on nmodes instead of 3*natom
730 : : !!
731 : : !! SOURCE
732 : :
733 : 13471 : pure subroutine phdispl_cart2red_nmodes(natom, nmodes, gprimd, displ_cart, displ_red)
734 : :
735 : : !Arguments ------------------------------------
736 : : !scalars
737 : : integer,intent(in) :: natom, nmodes
738 : : !arrays
739 : : real(dp),intent(in) :: gprimd(3,3)
740 : : real(dp),intent(in) :: displ_cart(2,3*natom, nmodes)
741 : : real(dp),intent(out) :: displ_red(2,3*natom, nmodes)
742 : :
743 : : !Local variables-------------------------
744 : : integer :: jbranch,iatom,idir,ibranch,kdir,k1
745 : : ! *************************************************************************
746 : :
747 [ + + + + : 866119 : displ_red = zero
+ + ]
748 : :
749 [ + + ]: 67675 : do jbranch=1,nmodes
750 : : !
751 [ + + ]: 156391 : do iatom=1,natom
752 [ + + ]: 409068 : do idir=1,3
753 : 266148 : ibranch=idir+3*(iatom-1)
754 [ + + ]: 1153308 : do kdir=1,3
755 : 798444 : k1 = kdir+3*(iatom-1)
756 : : ! WARNING: could be non-transpose of rprimd matrix : to be checked.
757 : : ! 23 june 2004: rprimd becomes gprimd. could be gprim and then multiply by acell...
758 : : ! Nope, checked and ok with gprimd 24 jun 2004
759 : 798444 : displ_red(1,ibranch,jbranch) = displ_red(1,ibranch,jbranch) + gprimd(kdir,idir) * displ_cart(1,k1,jbranch)
760 : 1064592 : displ_red(2,ibranch,jbranch) = displ_red(2,ibranch,jbranch) + gprimd(kdir,idir) * displ_cart(2,k1,jbranch)
761 : :
762 : : end do ! kdir
763 : : end do ! idir
764 : : end do ! iatom
765 : : end do ! jbranch
766 : :
767 : 13471 : end subroutine phdispl_cart2red_nmodes
768 : : !!***
769 : :
770 : : !----------------------------------------------------------------------
771 : :
772 : : !!****f* m_geometry/getspinrot
773 : : !! NAME
774 : : !! getspinrot
775 : : !!
776 : : !! FUNCTION
777 : : !! From the symmetry matrix symrel expressed in the coordinate system rprimd,
778 : : !! compute the components of the spinor rotation matrix.
779 : : !!
780 : : !! INPUTS
781 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
782 : : !! symrel(3,3)=symmetry operation in real space in terms of primitive translations rprimd
783 : : !!
784 : : !! OUTPUT
785 : : !! spinrot(4)=components of the spinor rotation matrix:
786 : : !!
787 : : !! spinrot(1)=$\cos \phi / 2$
788 : : !! spinrot(2)=$\sin \phi / 2 \times u_x$
789 : : !! spinrot(3)=$\sin \phi / 2 \times u_y$
790 : : !! spinrot(4)=$\sin \phi / 2 \times u_z$
791 : : !!
792 : : !! where $\phi$ is the angle of rotation, and $(u_x,u_y,u_z)$ is the normalized direction of the rotation axis
793 : : !!
794 : : !! NOTES
795 : : !! Only the proper part of the symmetry operation is taken into account:
796 : : !! pure rotations, while the inversion part is taken away, if present.
797 : : !! as inversion acts on spinors without affecting the spin index.
798 : : !!
799 : : !! The whole collection of symmetry matrices is call symrel(3,3,nsym)
800 : : !! symrel1 contains just one of those matrices symrel1(3,3)
801 : : !!
802 : : !! SOURCE
803 : :
804 : 759308 : subroutine getspinrot(rprimd, spinrot, symrel)
805 : :
806 : : !Arguments ------------------------------------
807 : : !arrays
808 : : integer,intent(in) :: symrel(3,3)
809 : : real(dp),intent(in) :: rprimd(3,3)
810 : : real(dp),intent(out) :: spinrot(4)
811 : :
812 : : !Local variables-------------------------------
813 : : !scalars
814 : : integer :: det
815 : : real(dp) :: cos_phi,norminv,phi,scprod,sin_phi
816 : : !character(len=500) :: msg
817 : : !arrays
818 : : integer :: identity(3,3),symrel1(3,3)
819 : : real(dp) :: axis(3),coord(3,3),coordinvt(3,3),matr1(3,3),matr2(3,3)
820 : : real(dp) :: rprimd_invt(3,3),vecta(3),vectb(3),vectc(3)
821 : : !**************************************************************************
822 : :
823 : 759308 : symrel1(:,:) = symrel(:,:)
824 : :
825 : : ! Compute determinant of the matrix
826 : 759308 : call mati3det(symrel1, det)
827 : :
828 : : ! Produce a rotation from an improper symmetry
829 [ + + + + : 5259116 : if (det==-1) symrel1(:,:) = -symrel1(:,:)
+ + ]
830 : :
831 : : ! Test the possibility of the unit matrix
832 : 759308 : identity(:,:)=0; identity(1,1)=1; identity(2,2)=1; identity(3,3)=1
833 : :
834 [ + + + + : 9871004 : if (sum((symrel1(:,:) - identity(:,:))**2)/=0) then
+ + ]
835 : :
836 : : ! Transform symmetry matrix in the system defined by rprimd
837 : 707516 : call matr3inv(rprimd, rprimd_invt)
838 [ + + + + ]: 9197708 : coord=TRANSPOSE(rprimd_invt)
839 : 707516 : call matr3inv(coord,coordinvt)
840 [ + + + + : 36790832 : matr1(:,:) = MATMUL(symrel1,coord)
+ + + + +
+ ]
841 [ + + + + : 28300640 : matr2(:,:) = MATMUL(TRANSPOSE(coordinvt),matr1)
+ + ]
842 : :
843 : : ! Find the eigenvector with unit eigenvalue of the rotation matrix in cartesian coordinate, matr2
844 : 707516 : matr1(:,:)=matr2(:,:)
845 : 707516 : matr1(1,1)=matr1(1,1)-one
846 : 707516 : matr1(2,2)=matr1(2,2)-one
847 : 707516 : matr1(3,3)=matr1(3,3)-one
848 : :
849 : : ! Compute the axis of rotation and the cos and sin of rotation angle
850 [ + + + + ]: 2830064 : if(DOT_PRODUCT(matr1(:,1),matr1(:,1)) < tol8 )then
851 : : ! The first direction is the axis
852 : 94180 : axis(1)=one ; axis(2)=zero ; axis(3)=zero
853 : 94180 : cos_phi=matr2(2,2)
854 : 94180 : sin_phi=matr2(3,2)
855 [ + + + + ]: 2453344 : else if(DOT_PRODUCT(matr1(:,2),matr1(:,2)) < tol8 )then
856 : : ! The second direction is the axis
857 : 90747 : axis(1)=zero ; axis(2)=one ; axis(3)=zero
858 : 90747 : cos_phi=matr2(3,3)
859 : 90747 : sin_phi=matr2(1,3)
860 : : else
861 : : ! In this case, try use the first and second vector to build the
862 : : ! rotation axis: compute their cross product
863 : 522589 : axis(1)=matr1(2,1)*matr1(3,2)-matr1(2,2)*matr1(3,1)
864 : 522589 : axis(2)=matr1(3,1)*matr1(1,2)-matr1(3,2)*matr1(1,1)
865 : 522589 : axis(3)=matr1(1,1)*matr1(2,2)-matr1(1,2)*matr1(2,1)
866 : : ! Then, try to normalize it
867 [ + + ]: 2090356 : scprod=DOT_PRODUCT(axis(:), axis(:))
868 [ + + ]: 522589 : if(scprod<tol8)then
869 : : ! The first and second vectors were linearly dependent
870 : : ! Thus, use the first and third vectors
871 : 71103 : axis(1)=matr1(2,1)*matr1(3,3)-matr1(2,3)*matr1(3,1)
872 : 71103 : axis(2)=matr1(3,1)*matr1(1,3)-matr1(3,3)*matr1(1,1)
873 : 71103 : axis(3)=matr1(1,1)*matr1(2,3)-matr1(1,3)*matr1(2,1)
874 : : ! Normalize the vector
875 [ + + ]: 284412 : scprod=DOT_PRODUCT(axis(:), axis(:))
876 [ - + ]: 71103 : if(scprod < tol8)then
877 : 0 : ABI_BUG('Cannot find the rotation axis.')
878 : : end if
879 : : end if
880 : 522589 : norminv=one/sqrt(scprod)
881 [ + + ]: 2090356 : axis(:)=axis(:)*norminv
882 : :
883 : : ! Project the axis vector out of the first unit vector,
884 : : ! and renormalize the projected vector
885 : : ! (the first vector cannot be the axis, as tested before)
886 : 522589 : vecta(1)=one-axis(1)**2
887 : 522589 : vecta(2)=-axis(1)*axis(2)
888 : 522589 : vecta(3)=-axis(1)*axis(3)
889 [ + + ]: 2090356 : scprod=DOT_PRODUCT(vecta(:),vecta(:))
890 : 522589 : norminv=one/sqrt(scprod)
891 [ + + ]: 2090356 : vecta(:)=vecta(:)*norminv
892 : : ! Rotate the vector A, to get vector B
893 [ + + ]: 2090356 : vectb(:)=matr2(:,1)*vecta(1)+matr2(:,2)*vecta(2)+matr2(:,3)*vecta(3)
894 : : ! Get dot product of vectors A and B, giving cos of the rotation angle
895 [ + + ]: 2090356 : cos_phi=DOT_PRODUCT(vecta(:),vectb(:))
896 : : ! Compute the cross product of the axis and vector A
897 : 522589 : vectc(1)=axis(2)*vecta(3)-axis(3)*vecta(2)
898 : 522589 : vectc(2)=axis(3)*vecta(1)-axis(1)*vecta(3)
899 : 522589 : vectc(3)=axis(1)*vecta(2)-axis(2)*vecta(1)
900 : : ! Get dot product of vectors B and C, giving sin of the rotation angle
901 [ + + ]: 2090356 : sin_phi=DOT_PRODUCT(vectb(:),vectc(:))
902 : : end if
903 : :
904 : : ! Get the rotation angle, then the parameters of the spinor rotation
905 : : ! Here, treat possible inaccurate values of the cosine of phi
906 [ + - ]: 707516 : if(cos_phi> one-tol8 )cos_phi= one-tol8
907 : : if(cos_phi<-(one-tol8))cos_phi=-(one-tol8)
908 : 707516 : phi=acos(cos_phi)
909 [ + + ]: 707516 : if(sin_phi<zero)phi=-phi
910 : : ! Rectify the angle, such that its absolute values corresponds to 180, 120, 90, 60, or 0 degrees
911 : 707516 : phi=(nint(six*phi/pi))/six*pi
912 : : ! Compute components of the spinor matrix
913 : 707516 : spinrot(1)=cos(half*phi)
914 : 707516 : spinrot(2)=axis(1)*sin(half*phi)
915 : 707516 : spinrot(3)=axis(2)*sin(half*phi)
916 : 707516 : spinrot(4)=axis(3)*sin(half*phi)
917 : :
918 : : else
919 : :
920 : : ! Here, the case of the unit matrix
921 : : axis(:)=zero
922 : 51792 : phi=zero
923 : 51792 : spinrot(1)=one
924 : 51792 : spinrot(2)=zero
925 : 51792 : spinrot(3)=zero
926 : 51792 : spinrot(4)=zero
927 : :
928 : : end if ! the case of the identity matrix
929 : :
930 : : !write(std_out,*)' getspinrot :'; write(std_out,*)' symre =',symrel(:,:)
931 : : !write(std_out,*)' symrel1 =',symrel1(:,:); write(std_out,*)' rprimd =',rprimd(:,:)
932 : : !write(std_out,*)' matr2 =',matr2(:,:); write(std_out,*)' matr1 =',matr1(:,:)
933 : : !write(std_out,*)' phi (degree)=',phi*180._dp/pi; write(std_out,'(a,3d16.6)' )' axis=',axis(:)
934 : : !write(std_out,*)' vecta=',vecta(:)
935 : : !stop
936 : :
937 : 759308 : end subroutine getspinrot
938 : : !!***
939 : :
940 : : !!****f* m_geometry/spinrot_cmat
941 : : !! NAME
942 : : !! spinrot_cmat
943 : : !!
944 : : !! FUNCTION
945 : : !! Construct 2x2 complex matrix representing the rotation operator in spin-space.
946 : : !!
947 : : !! INPUTS
948 : : !! spinrot(4)=components of the spinor rotation matrix computed by getspinrot
949 : : !!
950 : : !! OUTPUT
951 : : !! spinrot(2,2)=Rotation matrix (complex array)
952 : : !!
953 : : !! SOURCE
954 : :
955 [ # # ]: 0 : pure function spinrot_cmat(spinrot)
956 : :
957 : : !Arguments ------------------------------------
958 : : real(dp),intent(in) :: spinrot(4)
959 : : complex(dp) :: spinrot_cmat(2,2)
960 : : ! *************************************************************************
961 : :
962 : : ! Build rotation matrix from spinrot:
963 : : !
964 : : ! ( cos(phi/2) + i n_z sin(phi/2), (+n_y + i n_x) sin(phi/2) )
965 : : ! ( (-n_y + i n_x) sin(phi/2) , cos(phi/2) - i n_z sin(phi/2) )
966 : :
967 : : ! spinrot(1)=cos(half*phi)
968 : : ! spinrot(2)=axis(1)*sin(half*phi)
969 : : ! spinrot(3)=axis(2)*sin(half*phi)
970 : : ! spinrot(4)=axis(3)*sin(half*phi)
971 : :
972 : : ! Rotation in spinor space (same equations as in wfconv)
973 : : ! TODO: Be careful here as wfconv uses symrel^T to map k-points (listkk)
974 : : ! thus the inverse of the corresponding symrec.
975 : : ! This may explain why all the terms with sin(phi/2) change sign (phi --> -phi)
976 : :
977 : 0 : spinrot_cmat(1,1) = spinrot(1) + j_dpc*spinrot(4)
978 : 0 : spinrot_cmat(1,2) = spinrot(3) + j_dpc*spinrot(2)
979 : 0 : spinrot_cmat(2,1) =-spinrot(3) + j_dpc*spinrot(2)
980 : 0 : spinrot_cmat(2,2) = spinrot(1) - j_dpc*spinrot(4)
981 : :
982 : : ! My equation
983 : : !spinrot_cmat(1,1) = spinrot(1) - j_dpc*spinrot(4)
984 : : !spinrot_cmat(1,2) =-spinrot(3) - j_dpc*spinrot(2)
985 : : !spinrot_cmat(2,1) = spinrot(3) - j_dpc*spinrot(2)
986 : : !spinrot_cmat(2,2) = spinrot(1) + j_dpc*spinrot(4)
987 : :
988 : : end function spinrot_cmat
989 : : !!***
990 : :
991 : : !----------------------------------------------------------------------
992 : :
993 : : !!****f* m_geometry/rotmat
994 : : !! NAME
995 : : !! rotmat
996 : : !!
997 : : !! FUNCTION
998 : : !! Finds the rotation matrix.
999 : : !!
1000 : : !! INPUTS
1001 : : !! xaxis(3)= vectors defining the x axis
1002 : : !! zaxis(3)= vectors defining the z axis
1003 : : !!
1004 : : !! OUTPUT
1005 : : !! inversion_flag = flag that indicates that an inversion operation
1006 : : !! on the coordinate system should be done
1007 : : !! umat(3,3)= matrix that rotates the x=(1 0 0) and z=(0 0 1) to the new
1008 : : !! values defined in xaxis and zaxis
1009 : : !!
1010 : : !! NOTES
1011 : : !! Here I set that the axe x is originally at the 1 0 0 direction and z is originally 0 0 1.
1012 : : !! So calling rotmat(x',z') will find the rotation
1013 : : !! matrix for the case in which we rotate the x and z
1014 : : !! axes from their default values to x' and z'.
1015 : : !!
1016 : : !! SOURCE
1017 : :
1018 : 4 : subroutine rotmat(xaxis, zaxis, inversion_flag, umat)
1019 : :
1020 : : !Arguments ------------------------------------
1021 : : !scalars
1022 : : integer,intent(out) :: inversion_flag
1023 : : !arrays
1024 : : real(dp),intent(in) :: xaxis(3),zaxis(3)
1025 : : real(dp),intent(out) :: umat(3,3)
1026 : :
1027 : : !Local variables-------------------------------
1028 : : !scalars
1029 : : real(dp) :: cosine,xmod,zmod
1030 : : character(len=500) :: msg
1031 : : !arrays
1032 : : real(dp) :: yaxis(3)
1033 : : ! *************************************************************************
1034 : :
1035 [ + + + + : 16 : xmod = NORM2(xaxis(:))
- + ]
1036 [ + + + + : 16 : zmod = NORM2(zaxis(:))
- + ]
1037 : :
1038 [ - + ]: 4 : if(xmod < 1.d-8)then
1039 : : write(msg,'(a,a,a,i0)')&
1040 : 0 : 'The module of the xaxis should be greater than 1.d-8,',ch10,'however, |xaxis|=',xmod
1041 : 0 : ABI_BUG(msg)
1042 : : end if
1043 : :
1044 [ - + ]: 4 : if(zmod < 1.d-8)then
1045 : 0 : write(msg,'(a,a,a,i0)')'The module of the zaxis should be greater than 1.d-8,',ch10,'however, |zaxis|=',zmod
1046 : 0 : ABI_ERROR(msg)
1047 : : end if
1048 : :
1049 : : !verify that both axis are perpendicular
1050 [ + + ]: 16 : cosine = DOT_PRODUCT(xaxis,zaxis)/(xmod*zmod)
1051 : :
1052 [ - + ]: 4 : if(abs(cosine) > 1.d-8)then
1053 : 0 : write(msg,'(a,a,a,i6)')'xaxis and zaxis should be perpendicular,',ch10,'however, cosine=',cosine
1054 : 0 : ABI_BUG(msg)
1055 : : end if
1056 : :
1057 : : !new y axis as cross product
1058 : 4 : yaxis(1) = (zaxis(2)*xaxis(3) - xaxis(2)*zaxis(3))/(xmod*zmod)
1059 : 4 : yaxis(2) = (zaxis(3)*xaxis(1) - xaxis(3)*zaxis(1))/(xmod*zmod)
1060 : 4 : yaxis(3) = (zaxis(1)*xaxis(2) - xaxis(1)*zaxis(2))/(xmod*zmod)
1061 : :
1062 : : !hack to allow inversion operation on coordinate transformation
1063 : : !uses unlikely large but legal values of proj_x and/or proj_z
1064 : : !to flag inversion
1065 : 4 : inversion_flag=0
1066 [ + - - + ]: 4 : if(xmod>10._dp .or. zmod>10._dp) then
1067 : 0 : inversion_flag=1
1068 : : write(msg, '(4a)' )&
1069 : 0 : 'inversion operation will be appended to axis transformation',ch10,&
1070 : 0 : 'Action: If you did not intend this, make |z|<10 and |x|<10 ',ch10
1071 : 0 : call wrtout(std_out,msg)
1072 : : end if
1073 : :
1074 [ + + ]: 16 : umat(1,:) = xaxis(:)/xmod
1075 [ + + ]: 16 : umat(2,:) = yaxis(:)
1076 [ + + ]: 16 : umat(3,:) = zaxis(:)/zmod
1077 : :
1078 : 4 : end subroutine rotmat
1079 : : !!***
1080 : :
1081 : : !!****f* m_geometry/fixsym
1082 : : !! NAME
1083 : : !! fixsym
1084 : : !!
1085 : : !! FUNCTION
1086 : : !! Using input indsym which tells which atoms are related by symmetry,
1087 : : !! check that iatfix consistently fixes (freezes) all atoms which are
1088 : : !! related by symmetry, i.e. that iatfix does not break symmetry.
1089 : : !!
1090 : : !! INPUTS
1091 : : !! iatfix(3,natom)=integer array with 1 in every position for which
1092 : : !! the atom is to be kept fixed
1093 : : !! NOTE that this is not the input data structure for iatfix but it is
1094 : : !! the internal data structure used through most of the subroutines
1095 : : !! indsym(4,nsym,natom)=indirect indexing array for symmetrically related
1096 : : !! atoms; 4th element is label of symmetrically related atom
1097 : : !! natom=number of atoms
1098 : : !! nsym=number of symmetries (should be > 1 when this is called)
1099 : : !!
1100 : : !! OUTPUT
1101 : : !! (only checking)
1102 : : !!
1103 : : !! NOTE
1104 : : !! Stops execution with an error message if iatfix breaks symmetry.
1105 : : !!
1106 : : !! SOURCE
1107 : :
1108 : 4548 : subroutine fixsym(iatfix,indsym,natom,nsym)
1109 : :
1110 : : !Arguments ------------------------------------
1111 : : !scalars
1112 : : integer,intent(in) :: natom,nsym
1113 : : !arrays
1114 : : integer,intent(in) :: iatfix(3,natom),indsym(4,nsym,natom)
1115 : :
1116 : : !Local variables-------------------------------
1117 : : !scalars
1118 : : integer :: iatom,isym,jatom
1119 : : character(len=500) :: msg
1120 : : ! *************************************************************************
1121 : :
1122 [ + - ]: 4548 : if (nsym > 1) then
1123 [ + + ]: 14548 : do iatom=1,natom
1124 [ + + ]: 263950 : do isym=1,nsym
1125 : : ! jatom is the label of a symmetrically related atom
1126 : 249402 : jatom=indsym(4,isym,iatom)
1127 : : ! Thus the atoms jatom and iatom must be fixed along the same directions
1128 : : if (iatfix(1,jatom) /= iatfix(1,iatom) .or. &
1129 [ + - + - : 249402 : iatfix(2,jatom) /= iatfix(2,iatom) .or. &
- + ]
1130 : 10000 : iatfix(3,jatom) /= iatfix(3,iatom)) then
1131 : : write(msg, '(a,i0,a,i0,7a)' )&
1132 : 0 : 'Atom number: ',jatom,' is symmetrically equivalent to atom number: ',iatom,',',ch10,&
1133 : 0 : 'but according to iatfix, iatfixx, iatfixy and iatfixz, they',ch10,&
1134 : 0 : 'are not fixed along the same directions, which is forbidden.',ch10,&
1135 : 0 : 'Action: modify either the symmetry or iatfix(x,y,z) and resubmit.'
1136 : 0 : ABI_ERROR(msg)
1137 : : end if
1138 : : end do
1139 : : end do
1140 : : end if
1141 : :
1142 : 4548 : end subroutine fixsym
1143 : : !!***
1144 : :
1145 : : !!****f* m_geometry/det3r
1146 : : !! NAME
1147 : : !! det3r
1148 : : !!
1149 : : !! FUNCTION
1150 : : !! Compute determinant of a 3x3 real matrix
1151 : : !!
1152 : : !! SOURCE
1153 : :
1154 : 0 : pure real(dp) function det3r(rprimd)
1155 : :
1156 : : !Arguments ------------------------------------
1157 : : real(dp),intent(in) :: rprimd(3,3)
1158 : : ! *************************************************************************
1159 : :
1160 : : ! Compute unit cell volume
1161 : : det3r = rprimd(1,1)*(rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3))+&
1162 : : rprimd(2,1)*(rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3))+&
1163 : 0 : rprimd(3,1)*(rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3))
1164 : :
1165 : 0 : end function det3r
1166 : : !!***
1167 : :
1168 : : !!****f* m_geometry/metric
1169 : : !! NAME
1170 : : !! metric
1171 : : !!
1172 : : !! FUNCTION
1173 : : !! Compute first dimensional primitive translation vectors in reciprocal space
1174 : : !! gprimd from rprimd, and eventually writes out.
1175 : : !! Then, computes metrics for real and recip space rmet and gmet using length
1176 : : !! dimensional primitive translation vectors in columns of rprimd(3,3) and gprimd(3,3).
1177 : : !! gprimd is the inverse transpose of rprimd.
1178 : : !! i.e. $ rmet_{i,j}= \sum_k ( rprimd_{k,i}*rprimd_{k,j} ) $
1179 : : !! $ gmet_{i,j}= \sum_k ( gprimd_{k,i}*gprimd_{k,j} ) $
1180 : : !! Also computes unit cell volume ucvol in $\textrm{bohr}^3$
1181 : : !!
1182 : : !! INPUTS
1183 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1184 : : !! iout=unit number of output file. If iout<0, do not write output.
1185 : : !!
1186 : : !! OUTPUT
1187 : : !! gmet(3,3)=reciprocal space metric ($\textrm{bohr}^{-2}$).
1188 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
1189 : : !! rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).
1190 : : !! ucvol=unit cell volume ($\textrm{bohr}^{3}$).
1191 : : !!
1192 : : !! SOURCE
1193 : :
1194 : 1922801 : subroutine metric(gmet, gprimd, iout, rmet, rprimd, ucvol)
1195 : :
1196 : : !Arguments ------------------------------------
1197 : : !scalars
1198 : : integer,intent(in) :: iout
1199 : : real(dp),intent(out) :: ucvol
1200 : : !arrays
1201 : : real(dp),intent(in) :: rprimd(3,3)
1202 : : real(dp),intent(out) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1203 : :
1204 : : !Local variables-------------------------------
1205 : : !scalars
1206 : : integer :: nu
1207 : : character(len=500) :: msg
1208 : : !arrays
1209 : : real(dp) :: angle(3)
1210 : : ! *************************************************************************
1211 : :
1212 : : ! Compute unit cell volume
1213 : : ucvol=rprimd(1,1)*(rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3))+&
1214 : : rprimd(2,1)*(rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3))+&
1215 : 1922801 : rprimd(3,1)*(rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3))
1216 : : !ucvol = det3r(rprimd)
1217 : :
1218 : : ! Check that the input primitive translations are not linearly dependent (and none is zero); i.e. ucvol~=0
1219 : : ! Also ask that the mixed product is positive.
1220 [ - + ]: 1922801 : if (abs(ucvol)<tol12) then
1221 : : !write(std_out,*)"rprimd",rprimd,"ucvol",ucvol
1222 : : write(msg,'(6a,3(a,3es16.6,a))')&
1223 : 0 : 'Input rprim and acell gives vanishing unit cell volume.',ch10,&
1224 : 0 : 'This indicates linear dependency between primitive lattice vectors',ch10,&
1225 : 0 : 'Action: correct either rprim or acell in input file.', ch10, &
1226 : 0 : 'Rprimd =',rprimd(:,1),ch10,&
1227 : 0 : ' ',rprimd(:,2),ch10,&
1228 : 0 : ' ',rprimd(:,3),ch10
1229 : 0 : ABI_ERROR(msg)
1230 : : end if
1231 [ - + ]: 1922801 : if (ucvol<zero)then
1232 : : write(msg,'(2a,3(a,3es16.6,a),7a)')&
1233 : 0 : 'Current rprimd gives negative (R1 x R2) . R3 . ',ch10,&
1234 : 0 : 'Rprimd =',rprimd(:,1),ch10,&
1235 : 0 : ' ',rprimd(:,2),ch10,&
1236 : 0 : ' ',rprimd(:,3),ch10,&
1237 : 0 : 'Action: if the cell size and shape are fixed (optcell==0),',ch10,&
1238 : 0 : ' exchange two of the input rprim vectors;',ch10,&
1239 : 0 : ' if you are optimizing the cell size and shape (optcell/=0),',ch10,&
1240 : 0 : ' maybe the move was too large, and you might try to decrease strprecon.'
1241 : 0 : ABI_ERROR(msg)
1242 : : end if
1243 : :
1244 : : ! Generate gprimd
1245 : 1922801 : call matr3inv(rprimd, gprimd)
1246 : :
1247 : : ! Write out rprimd, gprimd and ucvol
1248 [ + + ]: 1922801 : if (iout>=0) then
1249 : 7909 : write(msg,'(2a)')' Real(R)+Recip(G) ','space primitive vectors, cartesian coordinates (Bohr,Bohr^-1):'
1250 : 7909 : call wrtout(iout,msg)
1251 [ + + ]: 31636 : do nu=1,3
1252 : : write(msg, '(1x,a,i1,a,3f11.7,2x,a,i1,a,3f11.7)' ) &
1253 [ + + ]: 94908 : 'R(',nu,')=',rprimd(:,nu)+tol10,&
1254 [ + + ]: 118635 : 'G(',nu,')=',gprimd(:,nu)+tol10
1255 : 31636 : call wrtout(iout,msg)
1256 : : end do
1257 : 7909 : write(msg,'(a,1p,e15.7,a)') ' Unit cell volume ucvol=',ucvol+tol10,' bohr^3'
1258 : 7909 : call wrtout(iout,msg)
1259 : 7909 : call wrtout(std_out,msg)
1260 : : end if
1261 : :
1262 : : ! Compute real space metric.
1263 [ + + + + : 76912040 : rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
+ + ]
1264 : :
1265 : : ! Compute reciprocal space metric.
1266 [ + + + + : 76912040 : gmet = MATMUL(TRANSPOSE(gprimd),gprimd)
+ + ]
1267 : :
1268 : : ! Write out the angles
1269 [ + + ]: 1922801 : if (iout>=0) then
1270 : 7909 : angle(1)=acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3)))/two_pi*360.0d0
1271 : 7909 : angle(2)=acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3)))/two_pi*360.0d0
1272 : 7909 : angle(3)=acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2)))/two_pi*360.0d0
1273 : 7909 : write(msg, '(a,3es16.8,a)' )' Angles (23,13,12)=',angle(1:3),' degrees'
1274 : 7909 : call wrtout(iout,msg)
1275 : 7909 : call wrtout(std_out,msg)
1276 : : end if
1277 : :
1278 : 1922801 : end subroutine metric
1279 : : !!***
1280 : :
1281 : : !!****f* m_geometry/mkradim
1282 : : !! NAME
1283 : : !! mkradim
1284 : : !!
1285 : : !! FUNCTION
1286 : : !! Not so trivial subroutine to make dimensionless real space
1287 : : !! primitive translations rprim(3,3) from dimensional rprimd(3).
1288 : : !! also make acell(3).
1289 : : !!
1290 : : !! INPUTS
1291 : : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
1292 : : !! where: rprimd(i,j)=rprim(i,j)*acell(j)
1293 : : !!
1294 : : !! OUTPUT
1295 : : !! acell(3)=unit cell length scales (bohr)
1296 : : !! rprim(3,3)=dimensionless real space primitive translations
1297 : : !!
1298 : : !! SOURCE
1299 : :
1300 : 513 : subroutine mkradim(acell,rprim,rprimd)
1301 : :
1302 : : !Arguments ------------------------------------
1303 : : !arrays
1304 : : real(dp),intent(out) :: acell(3),rprim(3,3)
1305 : : real(dp),intent(in) :: rprimd(3,3)
1306 : :
1307 : : !Local variables-------------------------------
1308 : : !scalars
1309 : : integer :: ii,jj
1310 : : real(dp) :: rprim_maxabs
1311 : : ! *************************************************************************
1312 : :
1313 : : !Use a representation based on normalised rprim vectors
1314 [ + + ]: 2052 : do ii=1,3
1315 [ + + + + : 6156 : acell(ii)=NORM2(rprimd(:,ii))
+ + ]
1316 [ + + ]: 6669 : rprim(:,ii)=rprimd(:,ii)/acell(ii)
1317 : : end do
1318 : :
1319 : : !Suppress meaningless values
1320 [ + + + + : 6669 : rprim_maxabs=maxval(abs(rprim))
+ + + - -
+ ]
1321 [ + + ]: 2052 : do ii=1,3
1322 [ + + ]: 6669 : do jj=1,3
1323 [ + + ]: 6156 : if(abs(rprim(ii,jj))<tol12*rprim_maxabs)rprim(ii,jj)=zero
1324 : : enddo
1325 : : enddo
1326 : :
1327 : 513 : end subroutine mkradim
1328 : : !!***
1329 : :
1330 : : !!****f* m_geometry/chkrprimd
1331 : : !!
1332 : : !! NAME
1333 : : !! chkrprimd
1334 : : !!
1335 : : !! FUNCTION
1336 : : !! Test if {rprim,acell,rprimd} are consistent
1337 : : !! It means that rprimd can be reconstructed from the rprim and acell
1338 : : !! Output a message if is not the case
1339 : : !!
1340 : : !! INPUTS
1341 : : !!
1342 : : !! OUTPUT
1343 : : !! (only writing)
1344 : : !!
1345 : : !! SOURCE
1346 : :
1347 : 0 : subroutine chkrprimd(acell,rprim,rprimd,iout)
1348 : :
1349 : : !Arguments ------------------------------------
1350 : : !scalars
1351 : : integer,intent(in) :: iout
1352 : : !arrays
1353 : : real(dp),intent(in) :: rprim(3,3)
1354 : : real(dp),intent(in) :: rprimd(3,3)
1355 : : real(dp),intent(in) :: acell(3)
1356 : :
1357 : : !Local variables-------------------------------
1358 : : !scalars
1359 : : integer :: ii,jj
1360 : : !arrays
1361 : : real(dp) :: rprimd_test(3,3)
1362 : : logical :: equal
1363 : :
1364 : : ! ***********************************************************
1365 : :
1366 : : !###########################################################
1367 : : !### 1. Compute rprimd from rprim and acell
1368 [ # # ]: 0 : do ii=1,3
1369 [ # # ]: 0 : rprimd_test(ii,1:3)=rprim(ii,1:3)*acell(1:3)
1370 : : end do
1371 : :
1372 : :
1373 : : !###########################################################
1374 : : !### 2. Compare rprimd and rprimd_test
1375 : :
1376 : : equal=.TRUE.
1377 [ # # ]: 0 : do ii=1,3
1378 [ # # ]: 0 : do jj=1,3
1379 [ # # ]: 0 : if (abs(rprimd_test(ii,jj)-rprimd(ii,jj))>1.E-12) then
1380 : 0 : equal=.FALSE.
1381 : : end if
1382 : : end do
1383 : : end do
1384 : :
1385 [ # # ]: 0 : if (equal)then
1386 : 0 : write(iout,*) 'chkrprimd: rprimd is consistent'
1387 : : else
1388 : 0 : write(iout,*) 'chkrprimd: rprimd is NOT consistent ERROR'
1389 : : end if
1390 : :
1391 : 0 : end subroutine chkrprimd
1392 : : !!***
1393 : :
1394 : : !!****f* m_geometry/chkdilatmx
1395 : : !! NAME
1396 : : !! chkdilatmx
1397 : : !!
1398 : : !! FUNCTION
1399 : : !! Check whether the new rprimd does not give a too large number
1400 : : !! of plane waves, compared to the one booked for rprimd, taking
1401 : : !! into account the maximal dilatation dilatmx. Actually check whether
1402 : : !! the new Fermi sphere is inside the old one, dilated.
1403 : : !!
1404 : : !! INPUTS
1405 : : !! chkdilatmx_ = if 1, will prevent to have any vector outside the Fermi sphere, possibly
1406 : : !! by rescaling (three times at most), and then stopping the execution
1407 : : !! if 0, simply send a warning, but continues execution
1408 : : !! dilatmx = maximal dilatation factor (usually the input variable)
1409 : : !! rprimd = new primitive vectors
1410 : : !! rprimd_orig = original primitive vectors (usually the input variable)
1411 : : !!
1412 : : !! OUTPUT
1413 : : !! dilatmx_errmsg=Empty string if calculation can continue.
1414 : : !! If the calculation cannot continue, dilatmx_errmsg will contain
1415 : : !! the message that should be reported in the output file.
1416 : : !!
1417 : : !! Client code should handle a possible problem with the following test:
1418 : : !!
1419 : : !! if (LEN_TRIM(dilatmx_errmsg) then
1420 : : !! dump dilatmx_errmsg to the main output file.
1421 : : !! handle_error
1422 : : !! end if
1423 : : !!
1424 : : !! SOURCE
1425 : :
1426 : 11540 : subroutine chkdilatmx(chkdilatmx_,dilatmx,rprimd,rprimd_orig,dilatmx_errmsg)
1427 : :
1428 : : !Arguments ------------------------------------
1429 : : !scalars
1430 : : integer,intent(in) :: chkdilatmx_
1431 : : real(dp),intent(in) :: dilatmx
1432 : : character(len=500),intent(out) :: dilatmx_errmsg
1433 : : !arrays
1434 : : real(dp),intent(inout) :: rprimd(3,3)
1435 : : real(dp),intent(in) :: rprimd_orig(3,3)
1436 : :
1437 : : !Local variables-------------------------------
1438 : : !scalars
1439 : : real(dp) :: alpha,dilatmx_new
1440 : : !arrays
1441 : : real(dp) :: eigval(3),gprimd_orig(3,3),met(3,3),old_to_new(3,3)
1442 : : character(len=500) :: msg
1443 : : ! *************************************************************************
1444 : :
1445 : : !Generates gprimd
1446 : 11540 : call matr3inv(rprimd_orig,gprimd_orig)
1447 : :
1448 : : !Find the matrix that transform an original xcart to xred, then to the new xcart
1449 [ + + + + : 461600 : old_to_new(:,:) = MATMUL(rprimd, TRANSPOSE(gprimd_orig))
+ + ]
1450 : :
1451 : : !The largest increase in length will be obtained thanks
1452 : : !to the diagonalization of the corresponding metric matrix :
1453 : : !it is the square root of its largest eigenvalue.
1454 [ + + + + : 461600 : met(:,:) = MATMUL(TRANSPOSE(old_to_new),old_to_new)
+ + ]
1455 : :
1456 : 11540 : call matr3eigval(eigval,met)
1457 : :
1458 [ - + + - : 57700 : dilatmx_new=sqrt(maxval(eigval(:)))
+ + ]
1459 : :
1460 : 11540 : dilatmx_errmsg = ""
1461 [ + + ]: 11540 : if(dilatmx_new>dilatmx+tol6)then
1462 : :
1463 : : ! MJV 2014 07 22: correct rprim to maximum jump allowed by dilatmx
1464 : : ! XG 20171011 : eigenvalues of "old_to_old" tensor are of course the unity !
1465 : :
1466 [ + + ]: 31 : if(chkdilatmx_/=0)then
1467 : 17 : alpha = (dilatmx - one) / (dilatmx_new - one)
1468 : : ! for safety, only 90 percent of max jump
1469 : 17 : alpha = 0.9_dp * alpha
1470 : :
1471 [ + + + + ]: 221 : rprimd = alpha * rprimd + (one - alpha) * rprimd_orig
1472 : :
1473 : : write(dilatmx_errmsg,'(3a,es16.6,4a,es16.6,2a,es16.6,a)')&
1474 : 17 : 'The new primitive vectors rprimd (an evolving quantity)',ch10,&
1475 : 17 : 'are too large with respect to the old rprimd and the accompanying dilatmx: ',dilatmx,ch10,&
1476 : 17 : 'This large change of unit cell parameters is not allowed by the present value of dilatmx.',ch10,&
1477 : 17 : 'An adequate value would have been dilatmx_new= ',dilatmx_new,ch10,&
1478 : 34 : 'Calculation continues with limited jump, by rescaling the projected move by the factor: ',alpha,'.'
1479 : : else
1480 : : write(msg, '(3a,es16.6,2a,es16.6,2a)' )&
1481 : 14 : 'The new primitive vectors rprimd (an evolving quantity)',ch10,&
1482 : 14 : 'are too large, given the initial rprimd and the accompanying dilatmx: ',dilatmx,ch10,&
1483 : 14 : 'An adequate value would have been dilatmx_new= ',dilatmx_new,ch10,&
1484 : 28 : 'As chkdilatmx=0, assume experienced user. Execution will continue.'
1485 : 14 : ABI_WARNING(msg)
1486 : : end if
1487 : :
1488 : : end if
1489 : :
1490 : 11540 : end subroutine chkdilatmx
1491 : : !!***
1492 : :
1493 : : !!****f* m_geometry/mkrdim
1494 : : !! NAME
1495 : : !! mkrdim
1496 : : !!
1497 : : !! FUNCTION
1498 : : !! Trivial subroutine to make dimensional real space
1499 : : !! primitive translations from length scales acell(3)
1500 : : !! and dimensionless translations rprim(3,3).
1501 : : !!
1502 : : !! INPUTS
1503 : : !! acell(3)=unit cell length scales (bohr)
1504 : : !! rprim(3,3)=dimensionless real space primitive translations
1505 : : !!
1506 : : !! OUTPUT
1507 : : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
1508 : : !! where: rprimd(i,j)=rprim(i,j)*acell(j)
1509 : : !!
1510 : : !! SOURCE
1511 : :
1512 : 91148 : subroutine mkrdim(acell,rprim,rprimd)
1513 : :
1514 : : !Arguments ------------------------------------
1515 : : !arrays
1516 : : real(dp),intent(in) :: acell(3),rprim(3,3)
1517 : : real(dp),intent(out) :: rprimd(3,3)
1518 : :
1519 : : !Local variables-------------------------------
1520 : : integer :: ii,jj
1521 : : ! *************************************************************************
1522 : :
1523 [ + + ]: 364592 : do ii=1,3
1524 [ + + ]: 1184924 : do jj=1,3
1525 : 1093776 : rprimd(ii,jj)=rprim(ii,jj)*acell(jj)
1526 : : end do
1527 : : end do
1528 : :
1529 : 91148 : end subroutine mkrdim
1530 : : !!***
1531 : :
1532 : : !!****f* m_geometry/xcart2xred
1533 : : !! NAME
1534 : : !! xcart2xred
1535 : : !!
1536 : : !! FUNCTION
1537 : : !! Convert from cartesian coordinates xcart(3,natom) in bohr to
1538 : : !! dimensionless reduced coordinates xred(3,natom) by using
1539 : : !! xred(mu,ia)=gprimd(1,mu)*xcart(1,ia)
1540 : : !! +gprimd(2,mu)*xcart(2,ia)
1541 : : !! +gprimd(3,mu)*xcart(3,ia)
1542 : : !! where gprimd is the inverse of rprimd
1543 : : !! Note that the reverse operation is done by xred2xcart
1544 : : !!
1545 : : !! INPUTS
1546 : : !! natom=number of atoms in unit cell
1547 : : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
1548 : : !! xcart(3,natom)=cartesian coordinates of atoms (bohr)
1549 : : !!
1550 : : !! OUTPUT
1551 : : !! xred(3,natom)=dimensionless reduced coordinates of atoms
1552 : : !!
1553 : : !! SOURCE
1554 : :
1555 : 24436 : subroutine xcart2xred(natom,rprimd,xcart,xred)
1556 : :
1557 : : !Arguments ------------------------------------
1558 : : !scalars
1559 : : integer,intent(in) :: natom
1560 : : !arrays
1561 : : real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
1562 : : real(dp),intent(out) :: xred(3,natom)
1563 : :
1564 : : !Local variables-------------------------------
1565 : : !scalars
1566 : : integer :: iatom,mu
1567 : : !arrays
1568 : : real(dp) :: gprimd(3,3)
1569 : : ! *************************************************************************
1570 : :
1571 : 24436 : call matr3inv(rprimd,gprimd)
1572 [ + + ]: 9288542 : do iatom=1,natom
1573 [ + + ]: 37080860 : do mu=1,3
1574 : 37056424 : xred(mu,iatom)= gprimd(1,mu)*xcart(1,iatom)+gprimd(2,mu)*xcart(2,iatom)+gprimd(3,mu)*xcart(3,iatom)
1575 : : end do
1576 : : end do
1577 : :
1578 : 24436 : end subroutine xcart2xred
1579 : : !!***
1580 : :
1581 : : !!****f* m_geometry/xred2xcart
1582 : : !! NAME
1583 : : !! xred2xcart
1584 : : !!
1585 : : !! FUNCTION
1586 : : !! Convert from dimensionless reduced coordinates xred(3,natom)
1587 : : !! to cartesian coordinates xcart(3,natom) in bohr by using
1588 : : !! xcart(mu,ia)=rprimd(mu,1)*xred(1,ia)
1589 : : !! +rprimd(mu,2)*xred(2,ia)
1590 : : !! +rprimd(mu,3)*xred(3,ia)
1591 : : !! Note that the reverse operation is done by xcart2xred.F90
1592 : : !!
1593 : : !! INPUTS
1594 : : !! natom=number of atoms in unit cell
1595 : : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
1596 : : !! xred(3,natom)=dimensionless reduced coordinates of atoms
1597 : : !!
1598 : : !! OUTPUT
1599 : : !! xcart(3,natom)=cartesian coordinates of atoms (bohr)
1600 : : !!
1601 : : !! SOURCE
1602 : :
1603 : 133188 : subroutine xred2xcart(natom, rprimd, xcart, xred)
1604 : :
1605 : : !Arguments ------------------------------------
1606 : : !scalars
1607 : : integer,intent(in) :: natom
1608 : : !arrays
1609 : : real(dp),intent(in) :: rprimd(3,3),xred(3,natom)
1610 : : real(dp),intent(out) :: xcart(3,natom)
1611 : :
1612 : : !Local variables-------------------------------
1613 : : integer :: iatom,mu
1614 : : ! *************************************************************************
1615 : :
1616 [ + + ]: 28107953 : do iatom=1,natom
1617 [ + + ]: 112032248 : do mu=1,3
1618 : 111899060 : xcart(mu,iatom)=rprimd(mu,1)*xred(1,iatom)+rprimd(mu,2)*xred(2,iatom)+rprimd(mu,3)*xred(3,iatom)
1619 : : end do
1620 : : end do
1621 : :
1622 : 133188 : end subroutine xred2xcart
1623 : : !!***
1624 : :
1625 : : !!****f* m_geometry/gred2fcart
1626 : : !! NAME
1627 : : !! gred2fcart
1628 : : !!
1629 : : !! FUNCTION
1630 : : !! Convert reduced forces into cartesian forces
1631 : : !!
1632 : : !! INPUTS
1633 : : !! gred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
1634 : : !! natom=Number of atoms in the unitary cell
1635 : : !! Favgz_null=TRUE if the average cartesian force has to be set to zero
1636 : : !! FALSE if it is set to zero only in x,y directions (not z)
1637 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1)
1638 : : !!
1639 : : !! OUTPUT
1640 : : !! fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)
1641 : : !!
1642 : : !! NOTES
1643 : : !! Unlike gred, fcart has been corrected by enforcing
1644 : : !! the translational symmetry, namely that the sum of force
1645 : : !! on all atoms is zero (except is a slab is used)
1646 : : !!
1647 : : !! SOURCE
1648 : :
1649 : 24016 : subroutine gred2fcart(favg,Favgz_null,fcart,gred,gprimd,natom)
1650 : :
1651 : : !Arguments ------------------------------------
1652 : : !scalars
1653 : : integer,intent(in) :: natom
1654 : : logical :: Favgz_null
1655 : : !arrays
1656 : : real(dp),intent(out) :: fcart(3,natom)
1657 : : real(dp),intent(in) :: gred(3,natom)
1658 : : real(dp),intent(in) :: gprimd(3,3)
1659 : : real(dp),intent(out) :: favg(3)
1660 : :
1661 : : !Local variables-------------------------------
1662 : : integer :: iatom,mu
1663 : : ! *************************************************************************
1664 : :
1665 : : !Note conversion to cartesian coordinates (bohr) AND
1666 : : !negation to make a force out of a gradient
1667 : 24016 : favg(:)=zero
1668 [ + + ]: 90147 : do iatom=1,natom
1669 [ + + ]: 288540 : do mu=1,3
1670 : : fcart(mu,iatom)= - (gprimd(mu,1)*gred(1,iatom)+&
1671 : : & gprimd(mu,2)*gred(2,iatom)+&
1672 : 198393 : & gprimd(mu,3)*gred(3,iatom))
1673 : 264524 : favg(mu)=favg(mu)+fcart(mu,iatom)
1674 : : end do
1675 : : end do
1676 : :
1677 : : !Subtract off average force from each force component
1678 : : !to avoid spurious drifting of atoms across cell.
1679 [ + + ]: 96064 : favg(:)=favg(:)/dble(natom)
1680 [ + + ]: 24016 : if(.not.Favgz_null) favg(3)=zero
1681 [ + + ]: 90147 : do iatom=1,natom
1682 [ + + ]: 288540 : fcart(:,iatom)=fcart(:,iatom)-favg(:)
1683 : : end do
1684 : :
1685 : 24016 : end subroutine gred2fcart
1686 : : !!***
1687 : :
1688 : : !!****f* m_geometry/fcart2gred
1689 : : !!
1690 : : !! NAME
1691 : : !! fcart2gred
1692 : : !!
1693 : : !! FUNCTION
1694 : : !! Convert cartesian forces into reduced forces
1695 : : !!
1696 : : !! INPUTS
1697 : : !! fcart(3,natom)=forces in cartesian coordinates (Ha/Bohr)
1698 : : !! natom=Number of atoms in the unitary cell
1699 : : !! rprimd(3,3)=dimensional primitive
1700 : : !!
1701 : : !! OUTPUT
1702 : : !! gred(3,natom)=symmetrized grtn = d(etotal)/d(xred)
1703 : : !!
1704 : : !! NOTES
1705 : : !! Unlike gred, fcart has been corrected by enforcing
1706 : : !! the translational symmetry, namely that the sum of force
1707 : : !! on all atoms is zero.
1708 : : !!
1709 : : !! SOURCE
1710 : :
1711 : 33159 : subroutine fcart2gred(fcart,gred,rprimd,natom)
1712 : :
1713 : : !Arguments ------------------------------------
1714 : : !scalars
1715 : : integer,intent(in) :: natom
1716 : : !arrays
1717 : : real(dp),intent(in) :: fcart(3,natom)
1718 : : real(dp),intent(out) :: gred(3,natom)
1719 : : real(dp),intent(in) :: rprimd(3,3)
1720 : :
1721 : : !Local variables-------------------------------
1722 : : integer :: iatom,mu
1723 : : ! *************************************************************************
1724 : :
1725 : : !MT, april 2012: the coding was not consistent with gred2fcart
1726 [ + + ]: 18244530 : do iatom=1,natom
1727 [ + + ]: 72878643 : do mu=1,3
1728 : : gred(mu,iatom)= - (rprimd(1,mu)*fcart(1,iatom)+&
1729 : : & rprimd(2,mu)*fcart(2,iatom)+&
1730 : 72845484 : & rprimd(3,mu)*fcart(3,iatom))
1731 : : end do
1732 : : end do
1733 : :
1734 : : !Previous version
1735 : : !do iatom=1,natom
1736 : : !do mu=1,3
1737 : : !gred(mu,iatom)= - (rprimd(mu,1)*fcart(1,iatom)+&
1738 : : !& rprimd(mu,2)*fcart(2,iatom)+&
1739 : : !& rprimd(mu,3)*fcart(3,iatom))
1740 : : !end do
1741 : : !end do
1742 : :
1743 : 33159 : end subroutine fcart2gred
1744 : : !!***
1745 : :
1746 : : !!****f* m_geometry/bonds_lgth_angles
1747 : : !! NAME
1748 : : !! bonds_lgth_angles
1749 : : !!
1750 : : !! FUNCTION
1751 : : !! From list of coordinates and primitive translations, output
1752 : : !! a list of bonds lengths and bond angles.
1753 : : !!
1754 : : !! INPUTS
1755 : : !! coordn = maximum coordination number to be taken into account
1756 : : !! fnameabo_app_geo=name of file for _GEO data
1757 : : !! natom = number of atoms in unit cell
1758 : : !! ntypat = number of types of atoms in unit cell.
1759 : : !! rprimd(3,3) = real space dimensional primitive translations (bohr)
1760 : : !! typat(natom) = type integer for each atom in cell
1761 : : !! znucl(ntypat)= real(dp), atomic number of atom type
1762 : : !! xred(3,natom)= reduced coordinates of atoms
1763 : : !!
1764 : : !! OUTPUT
1765 : : !! data written in file fnameabo_app_geo
1766 : : !!
1767 : : !! NOTES
1768 : : !! The tolerance tol8 aims at giving a machine-independent ordering.
1769 : : !! (this trick is used in bonds.f, listkk.f, prtrhomxmn.f and rsiaf9.f)
1770 : : !!
1771 : : !! SOURCE
1772 : :
1773 : 45 : subroutine bonds_lgth_angles(coordn,fnameabo_app_geo,natom,ntypat,rprimd,typat,xred,znucl)
1774 : :
1775 : : !Arguments ------------------------------------
1776 : : !scalars
1777 : : integer,intent(in) :: coordn,natom,ntypat
1778 : : character(len=*),intent(in) :: fnameabo_app_geo
1779 : : !arrays
1780 : : integer,intent(in) :: typat(natom)
1781 : : real(dp),intent(in) :: rprimd(3,3),znucl(ntypat)
1782 : : real(dp),intent(inout) :: xred(3,natom)
1783 : :
1784 : : !Local variables-------------------------------
1785 : : !scalars
1786 : : integer :: done,ia,ib,ic,ii,ineighb,jneighb,mneighb,mu,ndig,nu,t1,t2,t3,tmax,temp_unit
1787 : : real(dp) :: adotb,asq,bsq,co,length,sq,thdeg
1788 : : !real(dp)u1,u2,u3,v1,v2,v3
1789 : : character(len=500) :: msg
1790 : : type(atomdata_t) :: atom
1791 : : !arrays
1792 : 45 : integer,allocatable :: list_neighb(:,:,:)
1793 : : real(dp) :: bab(3),bac(3),dif(3),rmet(3,3)
1794 : 45 : real(dp),allocatable :: sqrlength(:),xcart(:,:)
1795 : 45 : character(len=8),allocatable :: iden(:)
1796 : : ! *************************************************************************
1797 : :
1798 : : !Initialize the file
1799 : 45 : write(msg, '(3a)' )' bonds_lgth_angles : about to open file ',trim(fnameabo_app_geo),ch10
1800 : 45 : call wrtout(std_out,msg); call wrtout(ab_out,msg)
1801 : :
1802 [ - + ]: 45 : if (open_file(fnameabo_app_geo,msg,newunit=temp_unit,status='unknown',form='formatted') /= 0) then
1803 : 0 : ABI_ERROR(msg)
1804 : : end if
1805 : 45 : rewind(temp_unit)
1806 : :
1807 : 45 : write(msg, '(a,a)' ) ch10,' ABINIT package : GEO file '
1808 : 45 : call wrtout(temp_unit,msg)
1809 : :
1810 : : !Compute maximum number of neighbors is the neighbor list,
1811 : : !from the indicative coordination number
1812 : : !Note : the following formula includes next nearest neighbors, but not others
1813 : 45 : mneighb=1+coordn+coordn*(coordn-1)
1814 : :
1815 : 45 : write(msg, '(a,a,i2,a,a,i4,a,a,a,i4,a)' ) ch10,&
1816 : 45 : & ' Maximal coordination number, as estimated by the user : ',coordn,ch10,&
1817 : 45 : & ' giving a maximum of ',coordn*coordn,&
1818 : 45 : & ' nearest neighbors and next nearest neighbors, ',ch10,&
1819 : 45 : & ' and ',(coordn*(coordn-1))/2,&
1820 : 90 : & ' distinct angles between nearest neighbors'
1821 : 45 : call wrtout(temp_unit,msg)
1822 : :
1823 : : !Compute metric tensor in real space rmet
1824 [ + + ]: 180 : do nu=1,3
1825 [ + + ]: 585 : do mu=1,3
1826 : : rmet(mu,nu)=rprimd(1,mu)*rprimd(1,nu)+&
1827 : : & rprimd(2,mu)*rprimd(2,nu)+&
1828 : 540 : & rprimd(3,mu)*rprimd(3,nu)
1829 : : end do
1830 : : end do
1831 : :
1832 : 45 : write(msg, '(a,a)' )ch10,' Primitive vectors of the periodic cell (bohr)'
1833 : 45 : call wrtout(temp_unit,msg)
1834 [ + + ]: 180 : do nu=1,3
1835 : 135 : write(msg, '(1x,a,i1,a,3f10.5)' ) ' R(',nu,')=',rprimd(:,nu)
1836 : 180 : call wrtout(temp_unit,msg)
1837 : : end do
1838 : :
1839 : 45 : write(msg, '(a,a)' ) ch10,&
1840 : 90 : & ' Atom list Reduced coordinates Cartesian coordinates (bohr)'
1841 : 45 : call wrtout(temp_unit,msg)
1842 : :
1843 : : !Set up a list of character identifiers for all atoms : iden(ia)
1844 [ - + + - : 135 : ABI_MALLOC(iden,(natom))
+ - + - -
+ - + ]
1845 [ + + ]: 271 : iden(:)=' '
1846 [ + + ]: 271 : do ia=1,natom
1847 : 226 : ndig=int(log10(dble(ia)+0.5d0))+1
1848 : 226 : call atomdata_from_znucl(atom,znucl(typat(ia)))
1849 [ + + ]: 226 : if(ndig==1) write(iden(ia), '(a,a,i1,a)' ) atom%symbol,'(',ia,') '
1850 [ + + ]: 226 : if(ndig==2) write(iden(ia), '(a,a,i2,a)' ) atom%symbol,'(',ia,') '
1851 [ - + ]: 226 : if(ndig==3) write(iden(ia), '(a,a,i3,a)' ) atom%symbol,'(',ia,') '
1852 [ - + ]: 226 : if(ndig==4) write(iden(ia), '(a,a,i4,a)' ) atom%symbol,'(',ia,')'
1853 [ - + ]: 497 : if(ndig>4)then
1854 : 0 : close(temp_unit)
1855 : : write(msg, '(a,i8,a,a)' )&
1856 : 0 : 'bonds_lgth_angles cannot handle more than 9999 atoms, while natom=',natom,ch10,&
1857 : 0 : 'Action: decrease natom, or contact ABINIT group.'
1858 : 0 : ABI_BUG(msg)
1859 : : end if
1860 : : end do
1861 : :
1862 : : !Compute cartesian coordinates, and print reduced and cartesian coordinates
1863 : : !then print coordinates in angstrom, with the format needed for xmol
1864 [ - + + - : 135 : ABI_MALLOC(xcart,(3,natom))
+ - + - -
+ - + ]
1865 : 45 : call xred2xcart(natom,rprimd,xcart,xred)
1866 : :
1867 [ + + ]: 271 : do ia=1,natom
1868 : : write(msg, '(a,a,3f10.5,a,3f10.5)' ) &
1869 [ + + ]: 1130 : ' ',iden(ia),(xred(ii,ia)+tol10,ii=1,3),&
1870 [ + - + - : 2034 : ' ',(xcart(ii,ia)+tol10,ii=1,3)
+ + ]
1871 : 271 : call wrtout(temp_unit,msg)
1872 : : end do
1873 : :
1874 : 45 : write(msg, '(a,a,a,a,i4,a)' )ch10,&
1875 : 90 : ' XMOL data : natom, followed by cartesian coordinates in Angstrom',ch10,ch10,natom,ch10
1876 : 45 : call wrtout(temp_unit,msg)
1877 : :
1878 [ + + ]: 271 : do ia=1,natom
1879 : 226 : call atomdata_from_znucl(atom,znucl(typat(ia)))
1880 [ + + ]: 904 : write(msg, '(a,a,3f10.5)' )' ',atom%symbol,xcart(1:3,ia)*Bohr_Ang
1881 : 497 : call wrtout(temp_unit,msg)
1882 : : end do
1883 : :
1884 [ - + ]: 45 : ABI_FREE(xcart)
1885 : :
1886 [ - + + - : 225 : ABI_MALLOC(list_neighb,(0:mneighb+1,4,2))
+ - + - +
- + - - +
- + ]
1887 [ - + + - : 135 : ABI_MALLOC(sqrlength,(0:mneighb+1))
+ - + - -
+ - + ]
1888 : :
1889 : : !Compute list of neighbors
1890 [ + + ]: 271 : do ia=1,natom
1891 : :
1892 : 226 : write(msg, '(a,a,a,a,a,a,a,a,a)' ) ch10,'===========',&
1893 : 226 : & '=====================================================================',&
1894 : 452 : & ch10,' ',iden(ia),ch10,ch10,' Bond lengths '
1895 : 226 : call wrtout(temp_unit,msg)
1896 : :
1897 : : ! Search other atoms for bonds, but must proceed
1898 : : ! in such a way to consider a search box sufficiently large,
1899 : : ! so increase the size of the search box until the
1900 : : ! final bond length list do not change
1901 [ - + ]: 648 : do tmax=0,5
1902 : :
1903 : : ! Set initial list of neighbors to zero,
1904 : : ! and initial square of bond lengths to a very large number.
1905 : : ! Note that the dimension is larger than neighb to ease
1906 : : ! the later sorting : neighbors 0 and neighb+1 are non-existent, while
1907 : : ! neighbor 1 will be the atom itself ...
1908 [ + + + + ]: 43848 : list_neighb(0:mneighb+1,1:4,1)=0
1909 [ + + ]: 10152 : sqrlength(1:mneighb+1)=huge(zero)
1910 : 648 : sqrlength(0)=-1.0d0
1911 : :
1912 : : ! Here search on all atoms inside the box defined by tmax
1913 [ + + ]: 10728 : do ib=1,natom
1914 [ + + ]: 40848 : do t3=-tmax,tmax
1915 [ + + ]: 157000 : do t2=-tmax,tmax
1916 [ + + ]: 656560 : do t1=-tmax,tmax
1917 : 509640 : dif(1)=xred(1,ia)-(xred(1,ib)+dble(t1))
1918 : 509640 : dif(2)=xred(2,ia)-(xred(2,ib)+dble(t2))
1919 : 509640 : dif(3)=xred(3,ia)-(xred(3,ib)+dble(t3))
1920 : 509640 : sq=rsdot(dif(1),dif(2),dif(3),dif(1),dif(2),dif(3),rmet)
1921 : :
1922 : : ! Insert the atom at the proper place in the neighbor list.
1923 [ + - ]: 958450 : do ineighb=mneighb,0,-1
1924 : : ! Note the tolerance
1925 [ + + ]: 841650 : if(sq+tol8>sqrlength(ineighb))then
1926 : 509640 : sqrlength(ineighb+1)=sq
1927 : 509640 : list_neighb(ineighb+1,1,1)=ib
1928 : 509640 : list_neighb(ineighb+1,2,1)=t1
1929 : 509640 : list_neighb(ineighb+1,3,1)=t2
1930 : 509640 : list_neighb(ineighb+1,4,1)=t3
1931 : : ! DEBUG
1932 : : ! if(ineighb/=mneighb)then
1933 : : ! write(std_out,*)' '
1934 : : ! do ii=1,mneighb
1935 : : ! write(std_out,*)ii,sqrlength(ii)
1936 : : ! end do
1937 : : ! end if
1938 : : ! ENDDEBUG
1939 : 509640 : exit
1940 : : else
1941 : 332010 : sqrlength(ineighb+1)=sqrlength(ineighb)
1942 [ + + ]: 1660050 : list_neighb(ineighb+1,1:4,1)=list_neighb(ineighb,1:4,1)
1943 : : end if
1944 : : end do
1945 : :
1946 : : end do
1947 : : end do
1948 : : end do
1949 : : ! end ib loop:
1950 : : end do
1951 : :
1952 : : ! Now, check that the box defined by tmax was large enough :
1953 : : ! require the present and old lists to be the same
1954 : 648 : done=0
1955 : :
1956 [ + + ]: 648 : if(tmax>0)then
1957 : : done=1
1958 [ + + ]: 6306 : do ineighb=1,mneighb
1959 : : ! DEBUG
1960 : : ! write(std_out,'(5i5,f12.5)' )ineighb,list_neighb(ineighb,1:4,1),&
1961 : : ! & sqrlength(ineighb)
1962 : : ! write(std_out,'(5i5)' )ineighb,list_neighb(ineighb,1:4,2)
1963 : : ! ENDDEBUG
1964 : : if( list_neighb(ineighb,1,1)/=list_neighb(ineighb,1,2) .or. &
1965 : : & list_neighb(ineighb,2,1)/=list_neighb(ineighb,2,2) .or. &
1966 [ + + + - : 5884 : & list_neighb(ineighb,3,1)/=list_neighb(ineighb,3,2) .or. &
+ + + + ]
1967 : 422 : & list_neighb(ineighb,4,1)/=list_neighb(ineighb,4,2) )then
1968 : 2531 : done=0
1969 : : end if
1970 : : end do
1971 : : end if
1972 : :
1973 : : ! If done==1, then one can exit the loop : the correct list of
1974 : : ! neighbors is contained in list_neighb(1:neighb,1:4,1),
1975 : : ! with the first neighbor being the atom itself
1976 [ + + ]: 422 : if(done==1)exit
1977 : :
1978 : : ! If the work is not done, while tmax==5, then there is a problem .
1979 [ - + ]: 422 : if(tmax==5)then
1980 : 0 : close(temp_unit)
1981 : : write(msg, '(2a)' )&
1982 : 0 : & 'Did not succeed to generate a reliable list of bonds ',&
1983 : 0 : & 'since tmax is exceeded.'
1984 : 0 : ABI_BUG(msg)
1985 : : end if
1986 : :
1987 : : ! Copy the new list into the old list.
1988 [ + + + + ]: 25872 : list_neighb(1:mneighb,1:4,2)=list_neighb(1:mneighb,1:4,1)
1989 : :
1990 : : ! Loop on tmax (note that there are exit instruction inside the loop)
1991 : : end do
1992 : :
1993 : :
1994 : :
1995 : : ! Output the bond list
1996 [ + + ]: 2972 : do ineighb=2,mneighb
1997 : 2746 : ib=list_neighb(ineighb,1,1)
1998 : 2746 : length=sqrt(sqrlength(ineighb))
1999 : : write(msg, '(a,a,a,a,3i2,t27,a,f10.5,a,f9.5,a)' )&
2000 [ + - + - ]: 2746 : & ' ',trim(iden(ia)),' - ',trim(iden(ib)),&
2001 : 2746 : & list_neighb(ineighb,2:4,1),'bond length is ',&
2002 : 5492 : & length,' bohr ( or ',Bohr_Ang*length,' Angst.)'
2003 : 2972 : call wrtout(temp_unit,msg)
2004 : : end do
2005 : :
2006 : : ! Output the angle list
2007 [ + + ]: 271 : if(coordn>1)then
2008 : :
2009 : 168 : write(msg, '(a,a)' ) ch10,' Bond angles '
2010 : 168 : call wrtout(temp_unit,msg)
2011 : :
2012 [ + + ]: 672 : do ineighb=2,coordn
2013 [ + + ]: 1680 : do jneighb=ineighb+1,coordn+1
2014 : :
2015 : 1008 : ib=list_neighb(ineighb,1,1)
2016 : 1008 : ic=list_neighb(jneighb,1,1)
2017 [ + + ]: 4032 : do mu=1,3
2018 : 3024 : bab(mu)=xred(mu,ib)+dble(list_neighb(ineighb,1+mu,1))-xred(mu,ia)
2019 : 4032 : bac(mu)=xred(mu,ic)+dble(list_neighb(jneighb,1+mu,1))-xred(mu,ia)
2020 : : end do
2021 : 1008 : asq=rsdot(bab(1),bab(2),bab(3),bab(1),bab(2),bab(3),rmet)
2022 : 1008 : bsq=rsdot(bac(1),bac(2),bac(3),bac(1),bac(2),bac(3),rmet)
2023 : 1008 : adotb=rsdot(bab(1),bab(2),bab(3),bac(1),bac(2),bac(3),rmet)
2024 : 1008 : co=adotb/sqrt(asq*bsq)
2025 [ - + ]: 1008 : if( abs(co)-1.0d0 >= 0.0d0 )then
2026 [ # # ]: 0 : if( abs(co)-1.0d0 <= 1.0d-12 )then
2027 : : ! Allows for a small numerical inaccuracy
2028 : 0 : thdeg=0.0d0
2029 [ # # ]: 0 : if(co < 0.0d0) thdeg=180.0d0
2030 : : else
2031 : 0 : ABI_BUG('the evaluation of the angle is wrong.')
2032 : : end if
2033 : : else
2034 : 1008 : thdeg=acos(co)*180.d0*piinv
2035 : : end if
2036 : :
2037 : : write(msg, '(a,a,3i2,a,a,a,a,3i2,t44,a,f13.5,a)' )&
2038 [ + - ]: 1008 : & ' ',trim(iden(ib)),list_neighb(ineighb,2:4,1),' - ',&
2039 [ + - + - ]: 1008 : & trim(iden(ia)),' - ',trim(iden(ic)),&
2040 : 2016 : & list_neighb(jneighb,2:4,1),'bond angle is ',thdeg,' degrees '
2041 : 1512 : call wrtout(temp_unit,msg)
2042 : : end do
2043 : : end do
2044 : :
2045 : : end if
2046 : : end do ! End big ia loop:
2047 : :
2048 [ - + ]: 45 : ABI_FREE(iden)
2049 [ - + ]: 45 : ABI_FREE(list_neighb)
2050 [ - + ]: 45 : ABI_FREE(sqrlength)
2051 : :
2052 : 45 : close(temp_unit)
2053 : :
2054 : : contains
2055 : :
2056 : 512664 : function rsdot(u1,u2,u3,v1,v2,v3,rmet)
2057 : :
2058 : : real(dp) :: rsdot
2059 : : real(dp),intent(in) :: u1,u2,u3,v1,v2,v3
2060 : : real(dp),intent(in) :: rmet(3,3)
2061 : : rsdot=rmet(1,1)*u1*v1+rmet(2,1)*u2*v1+&
2062 : : & rmet(3,1)*u3*v1+rmet(1,2)*u1*v2+rmet(2,2)*u2*v2+&
2063 : 512664 : & rmet(3,2)*u3*v2+rmet(1,3)*u1*v3+rmet(2,3)*u2*v3+rmet(3,3)*u3*v3
2064 : 512664 : end function rsdot
2065 : :
2066 : : end subroutine bonds_lgth_angles
2067 : : !!***
2068 : :
2069 : : !!****f* m_geometry/randomcellpos
2070 : : !! NAME
2071 : : !! randomcellpos
2072 : : !!
2073 : : !! FUNCTION
2074 : : !! This subroutine creates a unit cell with random atomic positions. It is
2075 : : !! assumed that the cell parameters are given and fixed. Several methods are
2076 : : !! used to generate the cell.
2077 : : !!
2078 : : !! INPUTS
2079 : : !! natom=number of atoms
2080 : : !! npsp=number of pseudopotentials (needed for the dimension of znucl)
2081 : : !! ntypat=number of type of atoms
2082 : : !! random_atpos=input variable
2083 : : !! 0 no generation of random atomic potision
2084 : : !! 1 completely random atomic potisions
2085 : : !! 2 random atomic positions, avoiding too close atoms
2086 : : !! (prevent coming closer than a fraction of the sum of covalent radii)
2087 : : !! 3 same than 2 but also generates the rprim and acell randomly
2088 : : !! within some given ranges (angles between 50 and 130)
2089 : : !! ratsph(1:ntypat)=radius of the atomic sphere
2090 : : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
2091 : : !! typat(1:natom)= input variable giving the type of each atom
2092 : : !! znucl(1:npsp)=nuclear number of atom as specified in psp file
2093 : : !!
2094 : : !! OUTPUT
2095 : : !! xred(3,natom)=reduced dimensionless atomic coordinates
2096 : : !!
2097 : : !! SIDE EFFECTS
2098 : : !!
2099 : : !! NOTES
2100 : : !!
2101 : : !! SOURCE
2102 : :
2103 : 7286 : subroutine randomcellpos(natom,npsp,ntypat,random_atpos,ratsph,rprim,rprimd,typat,xred,znucl,acell)
2104 : :
2105 : : !Arguments ------------------------------------
2106 : : !scalars
2107 : : integer,intent(in) :: natom,npsp,ntypat,random_atpos
2108 : : !arrays
2109 : : integer, intent(in) :: typat(natom)
2110 : : real(dp),intent(in) :: ratsph(ntypat)
2111 : : real(dp), intent(inout) :: rprim(3,3)
2112 : : real(dp), intent(inout) :: rprimd(3,3)
2113 : : real(dp), intent(inout) :: xred(3,natom)
2114 : : real(dp), intent(in) :: znucl(npsp)
2115 : : real(dp), intent(inout) :: acell(3)
2116 : :
2117 : : !Local variables-------------------------------
2118 : : integer :: iatom=0,ii,idum=-20
2119 : : real(dp) :: rij(3), rijd(3), radiuscovi, radiuscovj, dist, rati, ratj, angdeg(3)
2120 : : real(dp) :: cosang,aa,cc,a2
2121 : : character(len=500) :: msg
2122 : : type(atomdata_t) :: atom
2123 : : ! *************************************************************************
2124 : :
2125 : : !DEBUG
2126 : : !For the time being, print rprimd to keep it as an argument, in spite of abirule checking.
2127 : : !write (std_out,*) ' randomcellpos : enter'
2128 : : !write(std_out,*)' rprimd=',rprimd
2129 : : !write(std_out,*)' znucl=',znucl
2130 : : !write(std_out,*)' typat=',typat
2131 : : !write(std_out,*)' random_atpos=',random_atpos
2132 : : !ENDDEBUG
2133 : :
2134 [ + + - + ]: 7286 : if(random_atpos==2 .and. npsp/=ntypat)then
2135 : : write(msg, '(a,i5,2a,i5,a,i5,4a)' )&
2136 : 0 : & 'Input variable random_atpos= ',random_atpos,ch10,&
2137 : 0 : & 'However, the number of pseudopotentials ',npsp,', is not equal to the number of type of atoms ',ntypat,ch10,&
2138 : 0 : & 'The use of alchemical mixing cannot be combined with the constraint based on the mixing of covalent radii.',ch10,&
2139 : 0 : & 'Action: switch to another value of random_atpos.'
2140 : 0 : ABI_ERROR(msg)
2141 : : end if
2142 : :
2143 : : !random_atpos = 0 Default value, no random initialisation
2144 : : !random_atpos = 1 Fully random (Is it really useful ???)
2145 : : !random_atpos = 2 Random, but the sum of the two covalent radii is
2146 : : !less than the interatomic distance
2147 : : !random_atpos = 3 Random, but the sum of the two (other type of)
2148 : : !radii is less than the interatomic distance
2149 : : !random_atpos = 4 Random, but the sum of the two pseudopotential
2150 : : !radii is less than the interatomic distance
2151 : : !random_atpos = 5 Random, but the interatomic distance must be bigger
2152 : : !than the sum of
2153 : : !some input variable (well, instead of defining a new variable, why
2154 : : !not use ratsph ?)
2155 : : !Right now we are not using a factor for the tested distance.. something to be done, after a new variable has been defined
2156 : :
2157 [ + + ]: 7286 : if (random_atpos /= 0) then
2158 : : select case (random_atpos)
2159 : : case (1)
2160 [ # # ]: 0 : do ii=1,natom
2161 : 0 : xred(1,ii)=uniformrandom(idum)
2162 : 0 : xred(2,ii)=uniformrandom(idum)
2163 : 0 : xred(3,ii)=uniformrandom(idum)
2164 : : end do
2165 : : case (2)
2166 : 10 : iatom=0
2167 : 82 : do
2168 : 82 : iatom=iatom+1
2169 : 82 : xred(1,iatom)=uniformrandom(idum)
2170 : 82 : xred(2,iatom)=uniformrandom(idum)
2171 : 82 : xred(3,iatom)=uniformrandom(idum)
2172 : 82 : call atomdata_from_znucl(atom,znucl(typat(iatom)))
2173 : 82 : radiuscovi = atom%rcov
2174 [ + + ]: 366 : do ii=1,iatom-1
2175 [ + + ]: 1144 : rij=xred(:,iatom)-xred(:,ii)
2176 : : ! periodic boundary conditions
2177 [ + + ]: 1144 : rij = rij - 0.5
2178 [ + + ]: 1144 : rij = rij - anint (rij)
2179 : : ! coming back to cube between (0,1)
2180 [ + + ]: 1144 : rij = rij + 0.5
2181 : : ! convert reduced coordinates to cartesian coordinates
2182 : 286 : call xred2xcart(1,rprimd,rijd,rij)
2183 [ + + ]: 1144 : dist=dot_product(rijd,rijd)
2184 : 286 : call atomdata_from_znucl(atom,znucl(typat(ii)))
2185 : 286 : radiuscovj = atom%rcov
2186 [ + + ]: 366 : if (dist<(radiuscovj+radiuscovi)) then
2187 : 2 : iatom = iatom -1
2188 : 2 : EXIT
2189 : : end if
2190 : : end do
2191 [ + + ]: 82 : if (iatom>=natom) EXIT
2192 : : end do
2193 : : case(3)
2194 : 0 : iatom=0
2195 : 0 : do
2196 : 0 : iatom=iatom+1
2197 : 0 : xred(1,iatom)=uniformrandom(idum)
2198 : 0 : xred(2,iatom)=uniformrandom(idum)
2199 : 0 : xred(3,iatom)=uniformrandom(idum)
2200 : 0 : call atomdata_from_znucl(atom,znucl(typat(iatom)))
2201 : 0 : radiuscovi = atom%rcov
2202 [ # # ]: 0 : do ii=1,iatom-1
2203 [ # # ]: 0 : rij=xred(:,iatom)-xred(:,ii)
2204 : : ! periodic boundary conditions
2205 [ # # ]: 0 : rij = rij - 0.5
2206 [ # # ]: 0 : rij = rij - anint (rij)
2207 : : ! coming back to cube between (0,1)
2208 [ # # ]: 0 : rij = rij + 0.5
2209 : : ! convert reduced coordinates to cartesian coordinates
2210 : 0 : call xred2xcart(1,rprimd,rijd,rij)
2211 [ # # ]: 0 : dist=dot_product(rijd,rijd)
2212 : 0 : call atomdata_from_znucl(atom,znucl(typat(ii)))
2213 : 0 : radiuscovj = atom%rcov
2214 [ # # ]: 0 : if (dist<(radiuscovj+radiuscovi)) then
2215 : 0 : iatom = iatom -1
2216 : 0 : EXIT
2217 : : end if
2218 : : end do
2219 [ # # ]: 0 : if (iatom>=natom) EXIT
2220 : : end do
2221 [ # # ]: 0 : do ii=1,3
2222 : : ! generates cells with angles between 60 and 120 degrees
2223 : 0 : angdeg(ii)=60_dp+uniformrandom(idum)*60.0_dp
2224 : : end do
2225 [ # # ]: 0 : if (angdeg(1)+angdeg(2)+angdeg(3)>360._dp) then
2226 : 0 : angdeg(3)=360._dp-angdeg(1)-angdeg(2)
2227 : : end if
2228 : : ! check if angles are between the limits and create rprim
2229 : : if( abs(angdeg(1)-angdeg(2))<tol12 .and. &
2230 [ # # # # : 0 : & abs(angdeg(2)-angdeg(3))<tol12 .and. &
# # ]
2231 : : & abs(angdeg(1)-90._dp)+abs(angdeg(2)-90._dp)+abs(angdeg(3)-90._dp)>tol12 )then
2232 : : ! Treat the case of equal angles (except all right angles) :
2233 : : ! generates trigonal symmetry wrt third axis
2234 : 0 : cosang=cos(pi*angdeg(1)/180.0_dp)
2235 : 0 : a2=2.0_dp/3.0_dp*(1.0_dp-cosang)
2236 : 0 : aa=sqrt(a2)
2237 : 0 : cc=sqrt(1.0_dp-a2)
2238 : 0 : rprim(1,1)=aa ; rprim(2,1)=0.0_dp ; rprim(3,1)=cc
2239 : 0 : rprim(1,2)=-0.5_dp*aa ; rprim(2,2)= sqrt(3.0_dp)*0.5_dp*aa ; rprim(3,2)=cc
2240 : 0 : rprim(1,3)=-0.5_dp*aa ; rprim(2,3)=-sqrt(3.0_dp)*0.5_dp*aa ; rprim(3,3)=cc
2241 : : ! DEBUG
2242 : : ! write(std_out,*)' ingeo : angdeg=',angdeg(1:3)
2243 : : ! write(std_out,*)' ingeo : aa,cc=',aa,cc
2244 : : ! ENDDEBUG
2245 : : else
2246 : : ! Treat all the other cases
2247 : 0 : rprim(:,:)=0.0_dp
2248 : 0 : rprim(1,1)=1.0_dp
2249 : 0 : rprim(1,2)=cos(pi*angdeg(3)/180.0_dp)
2250 : 0 : rprim(2,2)=sin(pi*angdeg(3)/180.0_dp)
2251 : 0 : rprim(1,3)=cos(pi*angdeg(2)/180.0_dp)
2252 : 0 : rprim(2,3)=(cos(pi*angdeg(1)/180.0_dp)-rprim(1,2)*rprim(1,3))/rprim(2,2)
2253 : 0 : rprim(3,3)=sqrt(1.0_dp-rprim(1,3)**2-rprim(2,3)**2)
2254 : : end if
2255 : : ! generate acell
2256 : 0 : aa=zero
2257 [ # # ]: 0 : do ii=1,npsp
2258 : 0 : aa=znucl(ii)
2259 : : end do
2260 [ # # ]: 0 : do ii=1,3
2261 : 0 : acell(ii)=aa+uniformrandom(idum)*4.0
2262 : : end do
2263 : 0 : call mkrdim(acell,rprim,rprimd)
2264 : : case(4)
2265 : 0 : write(std_out,*) 'Not implemented yet'
2266 : : case(5)
2267 : 0 : iatom=0
2268 [ - + - - : 10 : do
- - ]
2269 : 0 : iatom=iatom+1
2270 : 0 : xred(1,iatom)=uniformrandom(idum)
2271 : 0 : xred(2,iatom)=uniformrandom(idum)
2272 : 0 : xred(3,iatom)=uniformrandom(idum)
2273 : 0 : rati=ratsph(typat(iatom))
2274 [ # # ]: 0 : do ii=1,iatom-1
2275 : 0 : ratj=ratsph(typat(ii))
2276 : : ! apply periodic boundary conditions
2277 [ # # ]: 0 : rij=(xred(:,iatom)-xred(:,ii))-0.5
2278 [ # # ]: 0 : rij = rij - ANINT ( rij )
2279 [ # # ]: 0 : rij = rij + 0.5
2280 : 0 : call xred2xcart(natom,rprimd,rijd,rij)
2281 [ # # ]: 0 : dist=dot_product(rijd,rijd)
2282 [ # # ]: 0 : if (dist<(rati+ratj)) EXIT
2283 : : end do
2284 [ # # ]: 0 : if (iatom==natom) EXIT
2285 [ # # ]: 0 : if (ii<(iatom-1)) iatom=iatom-1
2286 : : end do
2287 : : end select
2288 : : end if
2289 : :
2290 : 7286 : end subroutine randomcellpos
2291 : : !!***
2292 : :
2293 : : !!****f* m_geometry/shellstruct
2294 : : !! NAME
2295 : : !! shellstruct
2296 : : !!
2297 : : !! FUNCTION
2298 : : !! Calculates shell structure (multiplicities, radii) of an atomic configuration
2299 : : !!
2300 : : !! INPUTS
2301 : : !! natom=number of atoms in unit cell
2302 : : !! xred=reduced coordinates of atoms
2303 : : !! rprimd=unit cell vectors
2304 : : !! magv = magnetic ordering of atoms given as 1 and -1, if not given fm is assumed
2305 : : !! atp = atom on which the perturbation was done
2306 : : !!
2307 : : !! OUTPUT
2308 : : !! sdisv(nat)= distance of each shell to central atom (only the first nsh entries are relevant)
2309 : : !! nsh= number of shells
2310 : : !! mult(nat) = number of atoms on shell (only the first nsh entries are relevant)
2311 : : !!
2312 : : !! SOURCE
2313 : :
2314 : 55 : subroutine shellstruct(xred,rprimd,natom,magv,distv,smult,sdisv,nsh,atp,prtvol)
2315 : :
2316 : : !Arguments ------------------------------------
2317 : : !scalars
2318 : : integer,intent(in) :: natom
2319 : : integer,intent(in),optional :: atp
2320 : : integer,intent(in),optional :: prtvol
2321 : : integer,intent(out) :: nsh
2322 : : !arrays
2323 : : real(dp),intent(in) :: rprimd(3,3)
2324 : : real(dp),intent(in) :: xred(3,natom)
2325 : : integer,intent(out) :: smult(natom)
2326 : : integer,intent(in),optional :: magv(natom)
2327 : : real(dp),intent(out) :: sdisv(natom)
2328 : : real(dp),intent(out) :: distv(natom)
2329 : :
2330 : : !Local variables-------------------------------
2331 : : !scalars
2332 : : integer :: iatom,atpp,ish,prtvoll
2333 : : character(len=500) :: msg
2334 : : real(dp),parameter :: rndfact=10000_dp
2335 : : !arrays
2336 : 110 : integer :: iperm(natom),jperm(natom)
2337 : 110 : real(dp) :: distvh(natom,natom)
2338 : 110 : real(dp) :: magvv(natom)
2339 : : ! *************************************************************************
2340 : :
2341 [ + - ]: 55 : if (present(magv)) then
2342 [ + + ]: 4510 : magvv=magv
2343 : : else
2344 [ # # # # : 0 : magvv=(/ (1, iatom=1,natom) /)
# # ]
2345 : : end if
2346 : :
2347 [ + - ]: 55 : if (present(atp)) then
2348 : 55 : atpp=atp
2349 : : else
2350 : 0 : atpp=1
2351 : : end if
2352 : :
2353 [ + - ]: 55 : if (present(prtvol)) then
2354 : 55 : prtvoll=prtvol
2355 : : else
2356 : : prtvoll=1
2357 : : end if
2358 : :
2359 : : !DEBUB
2360 : 55 : write(std_out,*)'shellstruct start'
2361 : : !END DEBUG
2362 : :
2363 : : !Calculate ionic distances
2364 [ + - + + : 4510 : call ioniondist(natom,rprimd,xred,distvh,1,magv=int(magvv),atp=atpp)
- + ]
2365 [ + + ]: 4510 : distv=distvh(1,:)
2366 : :
2367 [ - + ]: 55 : if (prtvol>2) then
2368 : 0 : write(std_out,'(a)')' shellstruct ionic distances in cell (distv) : '
2369 : 0 : call prmat(distv(1:natom),1,natom,1,std_out)
2370 : : end if
2371 : :
2372 [ + + + + ]: 8965 : iperm=(/ (iatom, iatom=1,natom ) /)
2373 [ + + ]: 4510 : jperm=iperm
2374 [ + + ]: 4510 : distv=anint(distv*rndfact)/rndfact
2375 : : !Sort distances
2376 : 55 : call sort_dp(natom,distv,iperm,10d-5)
2377 : 55 : call sort_int(natom,iperm,jperm)
2378 : :
2379 [ + + ]: 4510 : smult=0
2380 [ + + + + : 5005 : sdisv=dot_product(rprimd(1,:),rprimd(1,:))+dot_product(rprimd(2,:),rprimd(2,:))+dot_product(rprimd(3,:),rprimd(3,:))
+ + + + ]
2381 : :
2382 : 55 : nsh=1
2383 : 55 : smult(1)=1
2384 : 55 : sdisv(1)=distv(1)
2385 : :
2386 [ + + ]: 4455 : do iatom=2,natom
2387 [ + - ]: 37025 : do ish=1,natom
2388 [ + + ]: 36970 : if (distv(iatom)>sdisv(ish)) then
2389 : : cycle
2390 [ + + ]: 4400 : else if (distv(iatom)==sdisv(ish)) then
2391 : 3920 : smult(ish)=smult(ish)+1
2392 : 3920 : exit
2393 [ + - ]: 480 : else if (distv(iatom)<sdisv(ish)) then
2394 [ + + ]: 65760 : smult(ish+1:natom)=smult(ish:natom-1)
2395 [ + + ]: 65760 : sdisv(ish+1:natom)=sdisv(ish:natom-1)
2396 : 480 : smult(ish)=1
2397 : 480 : sdisv(ish)=distv(iatom)
2398 : 480 : nsh=nsh+1
2399 : 480 : exit
2400 : : end if
2401 : : end do
2402 : : end do
2403 : :
2404 [ + - + + : 8965 : distv=(/ ( distv(jperm(iatom)),iatom=1,natom ) /)
+ + ]
2405 : :
2406 [ - + ]: 55 : if (prtvoll>2) then
2407 : 0 : write(msg,'(a,i4,a)')' shellstruct found ',nsh,' shells at distances (sdisv) '
2408 : 0 : call wrtout(std_out,msg)
2409 : 0 : call prmat(sdisv(1:nsh),1,nsh,1,std_out)
2410 : 0 : write(msg,fmt='(a,150i4)')' and multiplicities (smult) ', smult(1:nsh)
2411 : 0 : call wrtout(std_out,msg)
2412 : : end if
2413 : :
2414 : : !DEBUB
2415 : 55 : write(std_out,*)'shellstruct leave'
2416 : : !END DEBUG
2417 : :
2418 [ + - ]: 55 : end subroutine shellstruct
2419 : : !!***
2420 : :
2421 : : !!****f* m_geometry/ioniondist
2422 : : !! NAME
2423 : : !! ioniondist
2424 : : !!
2425 : : !! FUNCTION
2426 : : !! Compute ion-ion distances
2427 : : !!
2428 : : !! INPUTS
2429 : : !! natom= number of atoms in unit cell
2430 : : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
2431 : : !! xred(3,natom)=dimensionless reduced coordinates of atoms
2432 : : !! inm(natom,natom)=index (m,n) of the atom
2433 : : !! option= 1 output ion-ion distances / 2 output ordering of ion-ion
2434 : : !! distances / 3 output variables in varlist
2435 : : !! according to ion-ion distances * magnetic ordering
2436 : : !! magv magnetic ordering of atoms given also 1 and -1, if not
2437 : : !! given fm is assumed
2438 : : !! varlist=List of variables
2439 : : !! magv(natom)= magnetic ordering of atoms
2440 : : !! atp=atom on which the perturbation was done
2441 : : !!
2442 : : !! OUTPUT
2443 : : !!
2444 : : !! SOURCE
2445 : :
2446 : 165 : subroutine ioniondist(natom,rprimd,xred,inm,option,varlist,magv,atp,prtvol)
2447 : :
2448 : : !Arguments ------------------------------------
2449 : : !scalars
2450 : : integer,intent(in) :: natom,option
2451 : : integer,intent(in),optional :: atp !atom on which the perturbation was done
2452 : : !arrays
2453 : : real(dp),intent(in) :: rprimd(3,3)
2454 : : real(dp),intent(in) :: xred(3,natom)
2455 : : real(dp),intent(out) :: inm(natom,natom)
2456 : : integer,intent(in),optional :: magv(natom)
2457 : : real(dp),intent(in),optional :: varlist(natom)
2458 : : integer,intent(in),optional :: prtvol
2459 : :
2460 : : !Local variables-------------------------------
2461 : : !scalars
2462 : : integer :: iatom,jatom,katom,kdum,atpp,prtvoll
2463 : : !character(len=500) :: msg
2464 : : !arrays
2465 : 330 : integer :: interq(natom)
2466 : 330 : real(dp) :: hxcart(3,natom),distm(natom,natom)
2467 : 165 : real(dp) :: magvv(natom)
2468 : : ! *************************************************************************
2469 : :
2470 [ + - + + : 227370 : hxcart=matmul(rprimd,xred)
+ + + + +
+ + + ]
2471 [ + - + + : 26895 : interq=(/(iatom,iatom=1,natom)/)
+ + ]
2472 [ + + + + ]: 2251995 : inm=0
2473 : :
2474 [ + - ]: 165 : if (present(magv)) then
2475 [ + + ]: 13530 : magvv=magv
2476 : : else
2477 [ # # # # ]: 0 : magvv=(/ (1, iatom=1,natom) /)
2478 : : end if
2479 : :
2480 [ + - ]: 165 : if (present(atp)) then
2481 : 165 : atpp=atp
2482 : : else
2483 : : atpp=1
2484 : : end if
2485 : :
2486 [ + + ]: 165 : if (present(prtvol)) then
2487 : 110 : prtvoll=prtvol
2488 : : else
2489 : : prtvoll=1
2490 : : end if
2491 : :
2492 [ + + + - ]: 165 : if (option==3.and.(.not.present(varlist))) then
2493 : 0 : call wrtout(std_out,'ioniondist error: option=3 but no variable list provided for symmetrization')
2494 : : return
2495 : : end if
2496 : :
2497 : : !call wrtout(std_out,' ioniondist start ')
2498 : :
2499 [ + + + + ]: 2251995 : distm=0
2500 : 165 : katom=atpp-1
2501 [ + + ]: 13530 : do iatom=1,natom
2502 : 13365 : katom=katom+1
2503 [ - + ]: 13365 : if (katom > natom) katom=1
2504 : 13365 : distm(iatom,iatom)=0
2505 [ + + ]: 1139445 : do jatom=iatom,natom
2506 : 1125915 : distm(iatom,jatom)=dist2(xred(:,katom),xred(:,jatom),rprimd,1)*magvv(katom)*magvv(jatom)
2507 : 1139280 : distm(jatom,iatom)=distm(iatom,jatom)
2508 : : end do
2509 : : end do
2510 : :
2511 [ - + ]: 165 : if (prtvoll>=3) then
2512 : 0 : call wrtout(std_out,'ioniondist: ionic distances:')
2513 : 0 : call prmat(distm,natom,natom,natom,std_out)
2514 : : end if
2515 : :
2516 [ + + + + ]: 2251995 : distm=anint(distm*10000_dp)/10000_dp ! rounding needed else distm(iatom,jatom)/= distm(1,kdum) sometimes fails
2517 : :
2518 [ + + ]: 13530 : do iatom=1,natom
2519 [ + + ]: 13530 : if (option==1) then
2520 [ + + ]: 750610 : inm(iatom,:)=distm(iatom,:)
2521 : : else
2522 [ + + ]: 759520 : do jatom=iatom,natom
2523 : : kdum=1
2524 [ + - + + ]: 21266760 : do while ( (kdum <= natom) .and. (distm(iatom,jatom)/= distm(1,kdum)) )
2525 : 21266760 : kdum=kdum+1
2526 : : end do
2527 [ - + ]: 750610 : if (option==2) then
2528 : 0 : inm(iatom,jatom)=interq(kdum)
2529 [ + - ]: 750610 : else if (option==3) then
2530 : 750610 : inm(iatom,jatom)=varlist(kdum)
2531 : : end if
2532 : 759520 : inm(jatom,iatom)=inm(iatom,jatom)
2533 : : end do
2534 : : end if
2535 : : end do
2536 : :
2537 [ - + ]: 165 : if (prtvoll==2) then
2538 : 0 : call wrtout(std_out,'ioniondist: symmetrized matrix:')
2539 : 0 : call prmat(distm,1,natom,natom,std_out)
2540 [ - + ]: 165 : else if (prtvoll>=3) then
2541 : 0 : call wrtout(std_out,'ioniondist: symmetrized matrix:')
2542 : 0 : call prmat(distm,natom,natom,natom,std_out)
2543 : : end if
2544 : :
2545 [ + - ]: 165 : end subroutine ioniondist
2546 : : !!***
2547 : :
2548 : : !!****f* m_geometry/dist2
2549 : : !! NAME
2550 : : !! dist2
2551 : : !!
2552 : : !! FUNCTION
2553 : : !! Calculates the distance of v1 and v2 in a crystal by repeating the unit cell
2554 : : !!
2555 : : !! INPUTS
2556 : : !! v1,v2
2557 : : !! rprimd: dimensions of the unit cell. if not given 1,0,0/0,1,0/0,0,1 is assumed
2558 : : !! option: 0 v1, v2 given in cartesian coordinates (default)
2559 : : !! 1 v1,v2 given in reduced coordinates
2560 : : !! -1 v1 and v2 are supposed equal, and the routine returns the length of the smallest Bravais lattice vector
2561 : : !!
2562 : : !! OUTPUT
2563 : : !! dist2
2564 : : !!
2565 : : !! SOURCE
2566 : :
2567 : 1127887 : function dist2(v1,v2,rprimd,option)
2568 : :
2569 : : !Arguments ------------------------------------
2570 : : !scalars
2571 : : integer,intent(in),optional :: option
2572 : : real(dp) :: dist2
2573 : : !arrays
2574 : : real(dp),intent(in),optional :: rprimd(3,3)
2575 : : real(dp),intent(in) :: v1(3),v2(3)
2576 : :
2577 : : !Local variables-------------------------------
2578 : : !scalars
2579 : : integer :: i1,i2,i3,opt,s1,s2,s3
2580 : : real(dp):: min2,norm2,ucvol
2581 : : !arrays
2582 : : integer :: limits(3)
2583 : : real(dp) :: corner(3),dred(3),dtot(3),dv(3),dwrap(3),sh(3)
2584 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
2585 : : real(dp) :: vprimd(3,3)
2586 : : ! *************************************************************************
2587 : :
2588 [ - + ]: 1127887 : if (.not.PRESENT(rprimd)) then
2589 : 0 : vprimd=reshape((/1,0,0, 0,1,0, 0,0,1/),(/3,3/))
2590 : : else
2591 : 1127887 : vprimd=rprimd
2592 : : end if
2593 : :
2594 : 1127887 : call metric(gmet,gprimd,-1,rmet,vprimd,ucvol)
2595 : :
2596 [ + + ]: 4511548 : dv(:)=v2(:)-v1(:)
2597 : :
2598 : : !If in cartesian coordinates, need to be transformed to reduced coordinates.
2599 : 1127887 : opt=0
2600 [ + - ]: 1127887 : if(present(option))then
2601 : 1127887 : opt=option
2602 : : end if
2603 [ - + ]: 1127887 : if(opt==0)then
2604 [ # # ]: 0 : dred(:)=gprimd(1,:)*dv(1)+gprimd(2,:)*dv(2)+gprimd(3,:)*dv(3)
2605 [ + + ]: 1127887 : else if(opt==1)then
2606 : 1127833 : dred(:)=dv(:)
2607 [ + - ]: 54 : else if(opt==-1)then
2608 : 54 : dred(:)=zero
2609 : : end if
2610 : :
2611 : : !Wrap in the ]-1/2,1/2] interval
2612 : 1127887 : call wrap2_pmhalf(dred(1),dwrap(1),sh(1))
2613 : 1127887 : call wrap2_pmhalf(dred(2),dwrap(2),sh(2))
2614 : 1127887 : call wrap2_pmhalf(dred(3),dwrap(3),sh(3))
2615 : :
2616 : : !Compute the limits of the parallelipipedic box that contains the Wigner-Seitz cell
2617 : : !The reduced coordinates of the corners of the Wigner-Seitz cell are computed (multiplied by two)
2618 : : !Then, the maximal values of these reduced coordinates are stored.
2619 : 1127887 : limits(:)=0
2620 [ + + ]: 3383661 : do s1=-1,1,2
2621 [ + + ]: 7895209 : do s2=-1,1,2
2622 : 6767322 : do s3=-1,1,2
2623 [ + + ]: 36092384 : corner(:)=gmet(:,1)*s1*rmet(1,1)+gmet(:,2)*s2*rmet(2,2)+gmet(:,3)*s3*rmet(3,3)
2624 [ + - ]: 9023096 : limits(1)=max(limits(1),ceiling(abs(corner(1))+tol14))
2625 [ + - ]: 9023096 : limits(2)=max(limits(2),ceiling(abs(corner(2))+tol14))
2626 [ + - + + ]: 9023096 : limits(3)=max(limits(3),ceiling(abs(corner(3))+tol14))
2627 : : end do
2628 : : end do
2629 : : end do
2630 : :
2631 : : !Use all relevant primitive real space lattice vectors to find the minimal difference vector
2632 : 1127887 : min2=huge(zero)
2633 [ + + ]: 12852786 : do i1=-limits(1),limits(1)
2634 : 11724899 : dtot(1)=dwrap(1)+i1
2635 [ + + ]: 174928345 : do i2=-limits(2),limits(2)
2636 : 162075559 : dtot(2)=dwrap(2)+i2
2637 [ + + ]: 2742768725 : do i3=-limits(3),limits(3)
2638 [ + + + + ]: 2731043826 : if(opt/=-1.or.i1/=0.or.i2/=0.or.i3/=0)then
2639 : 2568968213 : dtot(3)=dwrap(3)+i3
2640 : : norm2=dtot(1)*rmet(1,1)*dtot(1)+dtot(2)*rmet(2,2)*dtot(2)+dtot(3)*rmet(3,3)*dtot(3)+&
2641 : 2568968213 : & 2*(dtot(1)*rmet(1,2)*dtot(2)+dtot(2)*rmet(2,3)*dtot(3)+dtot(3)*rmet(3,1)*dtot(1))
2642 : 2568968213 : min2=min(norm2,min2)
2643 : : endif
2644 : : end do
2645 : : end do
2646 : : end do
2647 : 1127887 : dist2=sqrt(min2)
2648 : :
2649 : 1127887 : end function dist2
2650 : : !!***
2651 : :
2652 : : !!****f* m_geometry/remove_inversion
2653 : : !! NAME
2654 : : !! remove_inversion
2655 : : !!
2656 : : !! FUNCTION
2657 : : !! Remove the inversion symmetry from a symmetry set as well
2658 : : !! all the improper rotations (if present)
2659 : : !!
2660 : : !! INPUTS
2661 : : !! nsym=initial number of symmetries
2662 : : !! symrel(3,3,nsym)=Initial set of symmetry operarations in real space
2663 : : !! tnons(3,nsym)=Initial fractional translations
2664 : : !!
2665 : : !! OUTPUT
2666 : : !! nsym_out=Number of symmetries in the set without improper rotation
2667 : : !! symrel_out(:,:) [pointer] = output symmetries without improper rotations
2668 : : !! tnons_out(:) [pointer] = fractional translations associated to symrel_out
2669 : : !! pinv=-1 if the inversion has been removed, 1 otherwise
2670 : : !!
2671 : : !! NOTES
2672 : : !! Note the use of pointers, memory is allocated inside the procedure and passed back
2673 : : !! to the caller. Thus memory deallocation is relegated to the caller. To be on the safe side
2674 : : !! the pointers should be nullified before entering.
2675 : : !!
2676 : : !! SOURCE
2677 : :
2678 : 15 : subroutine remove_inversion(nsym,symrel,tnons,nsym_out,symrel_out,tnons_out,pinv)
2679 : :
2680 : : !Arguments ------------------------------------
2681 : : !scalars
2682 : : integer,intent(in) :: nsym
2683 : : integer,intent(out) :: nsym_out,pinv
2684 : : !arrays
2685 : : integer,intent(in) :: symrel(3,3,nsym)
2686 : : integer,pointer :: symrel_out(:,:,:)
2687 : : real(dp),intent(in) :: tnons(3,nsym)
2688 : : real(dp),pointer :: tnons_out(:,:)
2689 : :
2690 : : !Local variables-------------------------------
2691 : : !scalars
2692 : : integer :: is,is2,is_discarded,is_inv,is_retained,nsym2
2693 : : logical :: found
2694 : : character(len=500) :: msg
2695 : : !arrays
2696 : 30 : integer :: determinant(nsym),inversion(3,3),symrel2(3,3,nsym)
2697 : 30 : real(dp) :: dtnons(3),tnons2(3,nsym)
2698 : : ! *********************************************************************
2699 : :
2700 : 15 : ABI_WARNING('Removing inversion related symmetrie from initial set')
2701 : :
2702 : : ! Find the occurrence of the inversion symmetry.
2703 [ + + + + ]: 210 : call set2unit(inversion) ; inversion=-inversion
2704 : :
2705 : 15 : is_inv=0; found=.FALSE.
2706 [ + + - + ]: 331 : do while (is_inv<nsym .and. .not.found)
2707 [ - + + + : 513 : is_inv=is_inv+1; found=ALL(symrel(:,:,is_inv)==inversion)
+ + ]
2708 : : end do
2709 [ - + ]: 15 : if (found) then
2710 : 0 : write(msg,'(a,i3)')' The inversion is symmetry operation no. ',is_inv
2711 : : else
2712 : 15 : write(msg,'(a)')' The inversion was not found in the symmetries list.'
2713 : : end if
2714 : 15 : call wrtout(std_out,msg)
2715 : :
2716 : : ! Find the symmetries that are related through the inversion symmetry
2717 : 15 : call symdet(determinant,nsym,symrel)
2718 : 15 : nsym2=0
2719 [ + + ]: 316 : do is=1,nsym-1
2720 [ + + ]: 3906 : do is2=is+1,nsym
2721 : :
2722 [ + + ]: 14360 : dtnons(:)=tnons(:,is2)-tnons(:,is)-tnons(:,is_inv)
2723 [ + - + + : 6125 : found=ALL(symrel(:,:,is)==-symrel(:,:,is2)).and.isinteger(dtnons,tol8)
+ + - + -
- ]
2724 : :
2725 : 301 : if (found) then
2726 : 0 : nsym2=nsym2+1
2727 : : ! Retain symmetries with positive determinant
2728 [ # # # # : 0 : if (ALL(tnons(:,is2)<tol8).and.ALL(tnons(:,is)<tol8)) then
# # # # #
# ]
2729 : 0 : is_retained=is2 ; is_discarded=is
2730 [ # # ]: 0 : if (determinant(is)==1) then
2731 : 0 : is_retained=is ; is_discarded=is2
2732 : : end if
2733 [ # # # # : 0 : else if (ALL(tnons(:,is2)<tol8)) then
# # ]
2734 : 0 : is_retained=is2 ; is_discarded=is
2735 : : else
2736 : 0 : is_retained=is ; is_discarded=is2
2737 : : end if
2738 : :
2739 [ # # # # ]: 0 : symrel2(:,:,nsym2)=symrel(:,:,is_retained)
2740 [ # # ]: 0 : tnons2 (:,nsym2)=tnons (:,is_retained)
2741 : : write(msg,'(a,i3,a,i3,3a,i3,a)')&
2742 : 0 : & ' Symmetry operations no. ',is,' and no. ',is2,&
2743 : 0 : & ' are related through the inversion.',ch10,&
2744 : 0 : & ' Symmetry operation no. ',is_discarded,' will be suppressed.'
2745 : 0 : call wrtout(std_out,msg)
2746 : : end if ! found
2747 : :
2748 : : end do !is2
2749 : : end do !is
2750 : :
2751 [ - + - - ]: 15 : if (nsym2/=(nsym/2).or.nsym==1) then
2752 : 15 : call wrtout(std_out, ' Program uses the original set of symmetries ')
2753 : 15 : nsym_out=nsym
2754 [ - + + - : 45 : ABI_MALLOC(symrel_out,(3,3,nsym))
+ - + - -
+ - + ]
2755 [ - + + - : 45 : ABI_MALLOC(tnons_out,(3,nsym))
+ - + - -
+ - + ]
2756 [ + + + + : 4123 : symrel_out(:,:,:)=symrel(:,:,1:nsym)
+ + ]
2757 [ + + + + ]: 1279 : tnons_out(:,:)=tnons(:,1:nsym)
2758 : 15 : pinv=1
2759 : : else
2760 : 0 : write(msg,'(a)')' Inversion related operations have been suppressed from symmetries list.'
2761 : 0 : call wrtout(std_out,msg)
2762 : 0 : nsym_out=nsym2
2763 [ # # # # : 0 : ABI_MALLOC(symrel_out,(3,3,nsym2))
# # # # #
# # # ]
2764 [ # # # # : 0 : ABI_MALLOC(tnons_out,(3,nsym2))
# # # # #
# # # ]
2765 [ # # # # : 0 : symrel_out(:,:,:)=symrel2(:,:,1:nsym2)
# # ]
2766 [ # # # # ]: 0 : tnons_out(:,:)=tnons(:,1:nsym2)
2767 : 0 : pinv=-1
2768 : : end if
2769 : :
2770 : 15 : end subroutine remove_inversion
2771 : : !!***
2772 : :
2773 : : !!****f* m_geometry/reduce2primitive
2774 : : !! NAME
2775 : : !! reduce2primitive
2776 : : !!
2777 : : !! FUNCTION
2778 : : !! Find real space primitive vectors from non-primitive ones and the set of non-integer translations
2779 : : !! that leave the system invariant
2780 : : !!
2781 : : !! INPUTS
2782 : : !! ntranslat=number of translations
2783 : : !! rprimd(3,3)=dimensional non-primitive vectors in real space (bohr)
2784 : : !! tolsym=tolerance for the symmetry operations
2785 : : !! translations(3,ntranslat)=translation vectors, in reduced coordinates
2786 : : !!
2787 : : !! OUTPUT
2788 : : !! rprimd_primitive(3,3)=dimensional primitive vectors in real space (bohr)
2789 : : !!
2790 : : !! SOURCE
2791 : :
2792 : 117 : subroutine reduce2primitive(ntranslat, rprimd, rprimd_primitive, tolsym, translations)
2793 : :
2794 : : !Arguments ------------------------------------
2795 : : !scalars
2796 : : integer,intent(in) :: ntranslat
2797 : : real(dp),intent(in) :: tolsym
2798 : : !arrays
2799 : : real(dp),intent(in) :: rprimd(3,3),translations(3,ntranslat)
2800 : : real(dp),intent(out) :: rprimd_primitive(3,3)
2801 : :
2802 : : !Local variables-------------------------------
2803 : : !scalars
2804 : : integer :: idir,itentative,itrans,replace
2805 : : character(len=500) :: msg
2806 : : !arrays
2807 : 234 : real(dp) :: trans_cart(3,ntranslat),trans_red(3,ntranslat)
2808 : : !**************************************************************************
2809 : :
2810 : : !These translations should form the primitive lattice when combined with the non-primitive vectors.
2811 : : !Each translation, in reduced coordinates, should be constituted of rational numbers.
2812 : : !They should pave the non-primitive cell homogeneously. The issue is to replace
2813 : : !at least one (or more) of the non-primitive vectors by one (or more) selected translation vectors among the list.
2814 : : !All translation vectors should be an integer linear combination of the vectors of the new basis.
2815 : :
2816 : 117 : rprimd_primitive(:,:)=rprimd(:,:)
2817 : :
2818 : : !First, the reduced coordinates of translation vectors are transferred to the [0,1[ interval
2819 [ + + + + ]: 1773 : trans_red(:,1:ntranslat)=translations(:,1:ntranslat)-nint(translations(:,1:ntranslat)-tolsym)
2820 : :
2821 : : !Then, one of the translation vectors with the smallest non-zero first coordinate will replace the first vector, if any.
2822 : : !Similarly for the three directions.
2823 [ + + ]: 468 : do idir=1,3
2824 : : replace=0
2825 [ + + ]: 1593 : do itrans=1,ntranslat
2826 [ + + ]: 1593 : if(trans_red(idir,itrans)>tolsym)then
2827 [ + + ]: 377 : if(replace==0)then
2828 : : replace=1 ; itentative=itrans
2829 : : else
2830 [ - + ]: 175 : if(trans_red(idir,itentative)>trans_red(idir,itrans)+tolsym)then
2831 : 202 : itentative=itrans
2832 : : endif
2833 : : endif
2834 : : endif
2835 : : enddo
2836 [ + + ]: 468 : if(replace==1)then
2837 : : ! Change the trans vectors to cartesian coordinates, using "old" rprimd
2838 : 202 : call xred2xcart(ntranslat,rprimd_primitive,trans_cart,trans_red)
2839 : : ! Replace rprimd vector with index idir by the selected trans_cart vector
2840 [ + + ]: 808 : rprimd_primitive(:,idir)=trans_cart(:,itentative)
2841 : : ! Change the translation vectors to new reduced coordinates using updated rprimd_primitive
2842 : 202 : call xcart2xred(ntranslat,rprimd_primitive,trans_cart,trans_red)
2843 : : !Transfer to the [0,1[ interval
2844 [ + + + + ]: 3234 : trans_red(:,1:ntranslat)=trans_red(:,1:ntranslat)-nint(trans_red(:,1:ntranslat)-tolsym)
2845 : : endif
2846 : : enddo ! idir
2847 : :
2848 : : !Now, check that all translation vectors have zero reduced coordinates.
2849 [ + + ]: 531 : do itrans=1,ntranslat
2850 [ + + ]: 1773 : do idir=1,3
2851 [ - + ]: 1656 : if (abs(trans_red(idir,itrans))>tolsym) then
2852 : : write(msg,'(5a)')&
2853 : 0 : 'Did not succeed to find primitive cell from non-primitive one.',ch10,&
2854 : 0 : 'Indeed, there remains a non-vanishing pure translation after reduction.',ch10,&
2855 : 0 : 'Action: this is a bug, contact ABINIT group. Then, use a primitive cell in your input file.'
2856 : 0 : ABI_ERROR(msg)
2857 : : end if
2858 : : enddo
2859 : : enddo
2860 : :
2861 : 117 : end subroutine reduce2primitive
2862 : : !!***
2863 : :
2864 : : !!****f* m_geometry/symredcart
2865 : : !! NAME
2866 : : !! symredcart
2867 : : !!
2868 : : !! FUNCTION
2869 : : !! Convert a symmetry operation from reduced coordinates (integers)
2870 : : !! to cartesian coordinates (reals). Can operate in real or reciprocal space
2871 : : !!
2872 : : !! INPUTS
2873 : : !! symred(3,3)=symmetry matrice in reduced coordinates (integers) (real or reciprocal space)
2874 : : !! aprim(3,3)=real or reciprocal space dimensional primitive translations (see below)
2875 : : !! bprim(3,3)=real or reciprocal space dimensional primitive translations (see below)
2876 : : !!
2877 : : !! OUTPUT
2878 : : !! symcart(3,3)=symmetry matrice in cartesian coordinates (reals)
2879 : : !!
2880 : : !! NOTES
2881 : : !! When aprim=rprimd and bprim=gprimd, the routine operates in real space (on a real space symmetry)
2882 : : !! When aprim=gprimd and bprim=rprimd, the routine operates in reciprocal space (on a real space symmetry)
2883 : : !!
2884 : : !! SOURCE
2885 : :
2886 : 801971 : subroutine symredcart(aprim,bprim,symcart,symred)
2887 : :
2888 : : !Arguments ------------------------------------
2889 : : !arrays
2890 : : integer,intent(in) :: symred(3,3)
2891 : : real(dp),intent(in) :: aprim(3,3),bprim(3,3)
2892 : : real(dp),intent(out) :: symcart(3,3)
2893 : :
2894 : : !Local variables-------------------------------
2895 : : !scalars
2896 : : integer :: ii,jj,kk
2897 : : real(dp) :: symtmp
2898 : : !arrays
2899 : : real(dp) :: work(3,3)
2900 : : ! *************************************************************************
2901 : :
2902 : 801971 : work=zero
2903 [ + + ]: 3207884 : do kk=1,3
2904 [ + + ]: 10425623 : do jj=1,3
2905 : 7217739 : symtmp=dble(symred(jj,kk))
2906 [ + + ]: 31276869 : do ii=1,3
2907 : 28870956 : work(ii,jj)=work(ii,jj)+bprim(ii,kk)*symtmp
2908 : : end do
2909 : : end do
2910 : : end do
2911 : :
2912 : : ! work = bprim * symred^T
2913 : :
2914 : 801971 : symcart=zero
2915 [ + + ]: 3207884 : do kk=1,3
2916 [ + + ]: 10425623 : do jj=1,3
2917 : 7217739 : symtmp=work(jj,kk)
2918 [ + + ]: 31276869 : do ii=1,3
2919 : : ! symcart = aprim * work^T = aprim * symred * bprim^T
2920 : 28870956 : symcart(ii,jj)=symcart(ii,jj)+aprim(ii,kk)*symtmp
2921 : : end do
2922 : : end do
2923 : : end do
2924 : :
2925 : 801971 : end subroutine symredcart
2926 : : !!***
2927 : :
2928 : : !!****f* m_geometry/strainsym
2929 : : !! NAME
2930 : : !! strainsym
2931 : : !!
2932 : : !! FUNCTION
2933 : : !! For given order of point group, symmetrizes the strain tensor,
2934 : : !! then produce primitive vectors based on the symmetrized strain.
2935 : : !!
2936 : : !! INPUTS
2937 : : !! nsym=order of group.
2938 : : !! rprimd(3,3)= primitive vectors, to be symmetrized
2939 : : !! rprimd0(3,3)= reference primitive vectors, already symmetrized
2940 : : !! symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
2941 : : !!
2942 : : !! OUTPUT
2943 : : !! rprimd_symm(3,3)= symmetrized primitive vectors
2944 : : !!
2945 : : !! SOURCE
2946 : :
2947 : 478 : subroutine strainsym(nsym,rprimd0,rprimd,rprimd_symm,symrel)
2948 : :
2949 : : use m_linalg_interfaces
2950 : :
2951 : : !Arguments ------------------------------------
2952 : : !scalars
2953 : : integer,intent(in) :: nsym
2954 : : !arrays
2955 : : integer,intent(in) :: symrel(3,3,nsym)
2956 : : real(dp),intent(in) :: rprimd(3,3),rprimd0(3,3)
2957 : : real(dp),intent(out) :: rprimd_symm(3,3)
2958 : :
2959 : : !Local variables-------------------------------
2960 : : !scalars
2961 : : integer :: isym
2962 : : !arrays
2963 : : integer :: symrel_it(3,3)
2964 : : real(dp) :: rprimd0_inv(3,3),strain(3,3),strain_symm(3,3),tmp_mat(3,3),symrel_db(3,3)
2965 : : !**************************************************************************
2966 : :
2967 : : !copy initial rprimd input and construct inverse
2968 : 478 : rprimd0_inv = rprimd0
2969 : 478 : call matrginv(rprimd0_inv,3,3)
2970 : :
2971 : : !define strain as rprimd = strain * rprimd0 (in cartesian frame)
2972 : : !so strain = rprimd * rprimd0^{-1}
2973 : : !transform to triclinic frame with rprimd0^{-1} * strain * rprimd0
2974 : : !giving strain as rprimd0^{-1} * rprimd
2975 : 478 : call dgemm('N','N',3,3,3,one,rprimd0_inv,3,rprimd,3,zero,strain,3)
2976 : :
2977 : : !loop over symmetry elements to obtain symmetrized strain matrix
2978 : 478 : strain_symm = zero
2979 [ + + ]: 8236 : do isym = 1, nsym
2980 : : ! this loop accumulates symrel^{-1}*strain*symrel into strain_symm
2981 : : ! mati3inv gives the inverse transpose of symrel
2982 : 7758 : call mati3inv(symrel(:,:,isym),symrel_it)
2983 [ + + + + ]: 100854 : symrel_db = dble(symrel(:,:,isym))
2984 : 7758 : call dgemm('N','N',3,3,3,one,strain,3,symrel_db,3,zero,tmp_mat,3)
2985 [ + + + + ]: 100854 : symrel_db = dble(symrel_it)
2986 : 8236 : call dgemm('T','N',3,3,3,one,symrel_db,3,tmp_mat,3,one,strain_symm,3)
2987 : : end do
2988 : :
2989 : : !normalize by number of symmetry operations
2990 [ + + + + ]: 6214 : strain_symm = strain_symm/dble(nsym)
2991 : :
2992 : : !this step is equivalent to r_new = r_old * strain * r_old^{-1} * r_old,
2993 : : !that is, convert strain back to cartesian frame and then multiply by r_old,
2994 : : !to get the r_new primitive vectors
2995 : :
2996 : 478 : call dgemm('N','N',3,3,3,one,rprimd0,3,strain_symm,3,zero,rprimd_symm,3)
2997 : :
2998 : 478 : end subroutine strainsym
2999 : : !!***
3000 : :
3001 : : !!****f* m_geometry/stresssym
3002 : : !! NAME
3003 : : !! stresssym
3004 : : !!
3005 : : !! FUNCTION
3006 : : !! For given order of point group, symmetrizes the stress tensor,
3007 : : !! in symmetrized storage mode and cartesian coordinates, using input
3008 : : !! 3x3 symmetry operators in reduced coordinates.
3009 : : !! symmetrized tensor replaces input tensor.
3010 : : !!
3011 : : !! INPUTS
3012 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1)
3013 : : !! nsym=order of group.
3014 : : !! sym(3,3,nsym)=symmetry operators (usually symrec=expressed in terms
3015 : : !! of action on reciprocal lattice primitive translations); integers.
3016 : : !!
3017 : : !! SIDE EFFECTS
3018 : : !! stress(6)=stress tensor, in cartesian coordinates, in symmetric storage mode
3019 : : !!
3020 : : !! SOURCE
3021 : :
3022 : 19996 : subroutine stresssym(gprimd,nsym,stress,sym)
3023 : :
3024 : : !Arguments ------------------------------------
3025 : : !scalars
3026 : : integer,intent(in) :: nsym
3027 : : !arrays
3028 : : integer,intent(in) :: sym(3,3,nsym)
3029 : : real(dp),intent(in) :: gprimd(3,3)
3030 : : real(dp),intent(inout) :: stress(6)
3031 : :
3032 : : !Local variables-------------------------------
3033 : : !scalars
3034 : : integer :: ii,isym,mu,nu
3035 : : real(dp) :: summ,tmp
3036 : : !arrays
3037 : : real(dp) :: rprimd(3,3),rprimdt(3,3),strfrac(6),tensor(3,3),tt(3,3)
3038 : : !*************************************************************************
3039 : :
3040 : : !Obtain matrix of real space dimensional primitive translations
3041 : : !(inverse transpose of gprimd), and its transpose
3042 : 19996 : call matr3inv(gprimd,rprimd)
3043 [ + + + + ]: 259948 : rprimdt=transpose(rprimd)
3044 : :
3045 : : !Compute stress tensor in reduced coordinates
3046 : : ! strfrac = rprimd^T * stress * rprimd
3047 : 19996 : call strconv(stress,rprimdt,strfrac)
3048 : :
3049 : : !Switch to full storage mode
3050 : 19996 : tensor(1,1)=strfrac(1)
3051 : 19996 : tensor(2,2)=strfrac(2)
3052 : 19996 : tensor(3,3)=strfrac(3)
3053 : 19996 : tensor(3,2)=strfrac(4)
3054 : 19996 : tensor(3,1)=strfrac(5)
3055 : 19996 : tensor(2,1)=strfrac(6)
3056 : 19996 : tensor(2,3)=tensor(3,2)
3057 : 19996 : tensor(1,3)=tensor(3,1)
3058 : 19996 : tensor(1,2)=tensor(2,1)
3059 : :
3060 : : ! these loops are useless - trivial action:
3061 : : ! tt = tensor / dble(nsym)
3062 : : ! tensor = zero
3063 [ + + ]: 79984 : do nu=1,3
3064 [ + + ]: 259948 : do mu=1,3
3065 : 179964 : tt(mu,nu)=tensor(mu,nu)/dble(nsym)
3066 : 239952 : tensor(mu,nu)=0.0_dp
3067 : : end do
3068 : : end do
3069 : :
3070 : : !loop over all symmetry operations:
3071 : : ! tensor = symrec * tt * symrec^T = symrec * rprimd^T * input * rprimd symrec^T
3072 : : ! TODO: this should be replaced by a little BLAS call or two
3073 [ + + ]: 352018 : do isym=1,nsym
3074 [ + + ]: 1348084 : do mu=1,3
3075 [ + + ]: 4316286 : do nu=1,3
3076 : : summ=0._dp
3077 [ + + ]: 11952792 : do ii=1,3
3078 : : tmp=tt(ii,1)*sym(nu,1,isym)+tt(ii,2)*sym(nu,2,isym)+&
3079 : 8964594 : & tt(ii,3)*sym(nu,3,isym)
3080 : 11952792 : summ=summ+sym(mu,ii,isym)*tmp
3081 : : end do
3082 : 3984264 : tensor(mu,nu)=tensor(mu,nu)+summ
3083 : : end do
3084 : : end do
3085 : : end do
3086 : :
3087 : : !Switch back to symmetric storage mode
3088 : 19996 : strfrac(1)=tensor(1,1)
3089 : 19996 : strfrac(2)=tensor(2,2)
3090 : 19996 : strfrac(3)=tensor(3,3)
3091 : 19996 : strfrac(4)=tensor(3,2)
3092 : 19996 : strfrac(5)=tensor(3,1)
3093 : 19996 : strfrac(6)=tensor(2,1)
3094 : :
3095 : : !Convert back stress tensor (symmetrized) in cartesian coordinates
3096 : : ! stress = gprimd * symrec * rprimd^T * input * rprimd symrec^T * gprimd^T
3097 : : ! symrec_cart = gprimd * symrec * rprimd^T
3098 : : ! sym_cart = symrec_cart^-1 ^T = rprimd * sym * gprimd^T
3099 : 19996 : call strconv(strfrac,gprimd,stress)
3100 : :
3101 : 19996 : end subroutine stresssym
3102 : : !!***
3103 : :
3104 : : !!****f* m_geometry/stress_voigt_to_mat
3105 : : !! NAME
3106 : : !! stress_voigt_to_mat
3107 : : !!
3108 : : !! FUNCTION
3109 : : !! Build 3x3 symmetric stress tensor from stress vector in Voigt notation.
3110 : : !!
3111 : : !! INPUTS
3112 : : !!
3113 : : !! OUTPUT
3114 : : !!
3115 : : !! SOURCE
3116 : :
3117 : 5505 : subroutine stress_voigt_to_mat(stress6, stress_mat)
3118 : :
3119 : : real(dp),intent(in) :: stress6(6)
3120 : : real(dp),intent(out) :: stress_mat(3,3)
3121 : :
3122 : 5505 : stress_mat(1,1) = stress6(1)
3123 : 5505 : stress_mat(2,2) = stress6(2)
3124 : 5505 : stress_mat(3,3) = stress6(3)
3125 : 5505 : stress_mat(2,3) = stress6(4)
3126 : 5505 : stress_mat(3,2) = stress6(4)
3127 : 5505 : stress_mat(1,3) = stress6(5)
3128 : 5505 : stress_mat(3,1) = stress6(5)
3129 : 5505 : stress_mat(1,2) = stress6(6)
3130 : 5505 : stress_mat(2,1) = stress6(6)
3131 : :
3132 : 5505 : end subroutine stress_voigt_to_mat
3133 : : !!***
3134 : :
3135 : : !!****f* m_geometry/strconv
3136 : : !! NAME
3137 : : !! strconv
3138 : : !!
3139 : : !! FUNCTION
3140 : : !! If original gprimd is input, convert from symmetric storage mode
3141 : : !! 3x3 tensor in reduced coordinates "frac" to symmetric storage mode
3142 : : !! symmetric tensor in cartesian coordinates "cart".
3143 : : !!
3144 : : !! INPUTS
3145 : : !! frac(6)=3x3 tensor in symmetric storage mode, reduced coordinates
3146 : : !! gprimd(3,3)=reciprocal space dimensional primitive translations (bohr^-1)
3147 : : !!
3148 : : !! OUTPUT
3149 : : !! cart(6)=symmetric storage mode for symmetric 3x3 tensor in cartesian coords.
3150 : : !!
3151 : : !! NOTES
3152 : : !! $cart(i,j)=G(i,a) G(j,b) frac(a,b)$
3153 : : !! "Symmetric" storage mode for 3x3 tensor is 6 element array with
3154 : : !! elements 11, 22, 33, 32, 31, and 21.
3155 : : !! "cart" may be same array as "frac".
3156 : : !! If rprimd transpose is input instead of gprimd, then convert tensor
3157 : : !! in cartesian coordinates to reduced coordinates
3158 : : !!
3159 : : !! SOURCE
3160 : :
3161 : 936385 : subroutine strconv(frac,gprimd,cart)
3162 : :
3163 : : !Arguments ------------------------------------
3164 : : !arrays
3165 : : real(dp),intent(in) :: frac(6),gprimd(3,3)
3166 : : real(dp),intent(inout) :: cart(6) ! alias of frac !vz_i
3167 : :
3168 : : !Local variables-------------------------------
3169 : : !scalars
3170 : : integer :: ii,jj
3171 : : !arrays
3172 : : real(dp) :: work1(3,3),work2(3,3)
3173 : : ! *************************************************************************
3174 : :
3175 : 936385 : work1(1,1)=frac(1)
3176 : 936385 : work1(2,2)=frac(2)
3177 : 936385 : work1(3,3)=frac(3)
3178 : 936385 : work1(3,2)=frac(4) ; work1(2,3)=frac(4)
3179 : 936385 : work1(3,1)=frac(5) ; work1(1,3)=frac(5)
3180 : 936385 : work1(2,1)=frac(6) ; work1(1,2)=frac(6)
3181 : :
3182 : : ! TODO: these are matmuls, replace or get BLAS
3183 : : ! work2 = work1 * gprimd^T
3184 [ + + ]: 3745540 : do ii=1,3
3185 [ + + ]: 11236620 : work2(:,ii)=zero
3186 [ + + ]: 12173005 : do jj=1,3
3187 [ + + ]: 36519015 : work2(:,ii)=work2(:,ii)+gprimd(ii,jj)*work1(:,jj)
3188 : : end do
3189 : : end do
3190 : :
3191 : : ! work1 = gprimd * work2 = gprimd * input * gprimd^T
3192 [ + + ]: 3745540 : do ii=1,3
3193 [ + + ]: 11236620 : work1(ii,:)=zero
3194 [ + + ]: 12173005 : do jj=1,3
3195 [ + + ]: 36519015 : work1(ii,:)=work1(ii,:)+gprimd(ii,jj)*work2(jj,:)
3196 : : end do
3197 : : end do
3198 : :
3199 : 936385 : cart(1)=work1(1,1)
3200 : 936385 : cart(2)=work1(2,2)
3201 : 936385 : cart(3)=work1(3,3)
3202 : 936385 : cart(4)=work1(2,3)
3203 : 936385 : cart(5)=work1(1,3)
3204 : 936385 : cart(6)=work1(1,2)
3205 : :
3206 : 936385 : end subroutine strconv
3207 : : !!***
3208 : :
3209 : : !!****f* m_geometry/littlegroup_pert
3210 : : !!
3211 : : !! NAME
3212 : : !! littlegroup_pert
3213 : : !!
3214 : : !! FUNCTION
3215 : : !! If syuse==0 and abs(rfmeth)==2, determines the set of symmetries that leaves a perturbation invariant.
3216 : : !! (Actually, all symmetries that leaves a q-wavevector invariant should be used to reduce the number
3217 : : !! of k-points for all perturbations. Unfortunately, one has to take into account the sign reversal of the
3218 : : !! perturbation under the symmetry operations, which makes GS routines not usable for the respfn code.
3219 : : !! The intermediate choice was to select only those that keep also the perturbation invariant.
3220 : : !! Note that the wavevector of the perturbation must also be invariant,
3221 : : !! a translation vector in real space is NOT allowed ).
3222 : : !!
3223 : : !! INPUTS
3224 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space (bohr**-1)
3225 : : !! idir=direction of the perturbation
3226 : : !! indsym(4,nsym,natom)=indirect indexing of atom labels--see subroutine symatm for definition (if nsym>1)
3227 : : !! iout=if non-zero, output on unit iout
3228 : : !! ipert=characteristics of the perturbation
3229 : : !! natom= number of atoms
3230 : : !! nsym=number of space group symmetries
3231 : : !! rfmeth =
3232 : : !! 1 or -1 if non-stationary block
3233 : : !! 2 or -2 if stationary block
3234 : : !! 3 or -3 if third order derivatives
3235 : : !! positive if symmetries are used to set elements to zero whenever possible, negative to prevent this to happen.
3236 : : !! symq(4,2,nsym)= Table computed by littlegroup_q.
3237 : : !! three first numbers define the G vector;
3238 : : !! fourth number is zero if the q-vector is not preserved, is 1 otherwise
3239 : : !! second index is one without time-reversal symmetry, two with time-reversal symmetry
3240 : : !! symafm(nsym)=(anti)ferromagnetic part of the symmetry operations
3241 : : !! symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
3242 : : !! symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
3243 : : !! syuse= flag to use the symmetries or not. If 0 usei it, if 1 do not use it.
3244 : : !! tnons(3,nsym)=nonsymmorphic translations of space group in terms
3245 : : !! of real space primitive translations (may be 0)
3246 : : !! [unit]=By default the routine writes to std_out and this is very annoying if we are inside a big loop.
3247 : : !! Use unit=dev_null or a negative integer to disable writing.
3248 : : !!
3249 : : !! OUTPUT
3250 : : !! nsym1 =number of space group symmetries that leaves the perturbation invariant
3251 : : !! symaf1(nsym1)=(anti)ferromagnetic part of the corresponding symmetry operations
3252 : : !! symrl1(3,3,nsym1)=corresponding 3x3 matrices of the group symmetries (real space)
3253 : : !! tnons1(3,nsym1)=corresponding nonsymmorphic translations of space group in terms
3254 : : !! of real space primitive translations (may be 0)!!
3255 : : !!
3256 : : !! SOURCE
3257 : :
3258 : 4129 : subroutine littlegroup_pert(gprimd,idir,indsym,iout,ipert,natom,nsym,nsym1, &
3259 : 4129 : & rfmeth,symafm,symaf1,symq,symrec,symrel,symrl1,syuse,tnons,tnons1, &
3260 : : & unit) ! Optional
3261 : :
3262 : : !Arguments -------------------------------
3263 : : !scalars
3264 : : integer,intent(in) :: idir,iout,ipert,natom,nsym,rfmeth,syuse
3265 : : integer,intent(in),optional :: unit
3266 : : integer,intent(out) :: nsym1
3267 : : !arrays
3268 : : integer,intent(in) :: indsym(4,nsym,natom),symafm(nsym),symq(4,2,nsym)
3269 : : integer,intent(in) :: symrec(3,3,nsym),symrel(3,3,nsym)
3270 : : integer,intent(out) :: symaf1(nsym),symrl1(3,3,nsym)
3271 : : real(dp),intent(in) :: gprimd(3,3),tnons(3,nsym)
3272 : : real(dp),intent(out) :: tnons1(3,nsym)
3273 : :
3274 : : !Local variables -------------------------
3275 : : !scalars
3276 : : integer :: idir1,ii,istr,isym,jj,nsym_test,tok,ount
3277 : : character(len=500) :: msg
3278 : : !arrays
3279 : : integer :: sym_test(3,3,2)
3280 : : real(dp) :: str_test(6)
3281 : : ! *********************************************************************
3282 : :
3283 [ + + ]: 4129 : ount = std_out; if (present(unit)) ount = unit
3284 : :
3285 : 4129 : nsym1=0
3286 [ + + + + : 4129 : if((ipert==natom+3 .or. ipert==natom+4) .and. syuse==0 .and. abs(rfmeth)==2) then
+ - + - ]
3287 : : ! Strain perturbation section
3288 : : ! Use ground state routine which symmetrizes cartesian stress as a quick
3289 : : ! and dirty test for the invariance of the strain (ipert,idir) under
3290 : : ! each candidate symmetry
3291 : : ! I am presently assuming that translations are acceptable because I dont
3292 : : ! see why not.
3293 : :
3294 : 156 : istr=3*(ipert-natom-3)+idir
3295 : 156 : nsym_test=2
3296 : : ! Store identity as first element for test
3297 [ + + + + ]: 2028 : sym_test(:,:,1)=0
3298 : 156 : sym_test(1,1,1)=1; sym_test(2,2,1)=1; sym_test(3,3,1)=1
3299 [ + + ]: 3686 : do isym=1,nsym
3300 [ + + + + ]: 45890 : sym_test(:,:,2)=symrec(:,:,isym)
3301 : 3530 : str_test(:)=0.0_dp
3302 : 3530 : str_test(istr)=1.0_dp
3303 : 3530 : call stresssym(gprimd,nsym_test,str_test,sym_test)
3304 [ + + ]: 3686 : if(abs(str_test(istr)-1.0_dp)<tol8)then
3305 : : ! The test has been successful !
3306 : 1004 : nsym1=nsym1+1
3307 : 1004 : symaf1(nsym1)=symafm(isym)
3308 [ + + ]: 4016 : do ii=1,3
3309 : 3012 : tnons1(ii,nsym1)=tnons(ii,isym)
3310 [ + + ]: 13052 : do jj=1,3
3311 : 12048 : symrl1(ii,jj,nsym1)=symrel(ii,jj,isym)
3312 : : end do
3313 : : end do
3314 : : end if
3315 : : end do
3316 : :
3317 [ + + + - : 3973 : else if(ipert>natom .or. syuse/=0 .or. abs(rfmeth)/=2)then
+ - ]
3318 : :
3319 : : ! Not yet coded for d/dk or electric field perturbations
3320 : 716 : nsym1=1
3321 [ + + ]: 2864 : do ii=1,3
3322 : 2148 : tnons1(ii,1)=0._dp
3323 : 2148 : symaf1(1)=1
3324 [ + + ]: 9308 : do jj=1,3
3325 : 6444 : symrl1(ii,jj,1)=0
3326 [ + + ]: 8592 : if(ii==jj)symrl1(ii,jj,1)=1
3327 : : end do
3328 : : end do
3329 : :
3330 : : else
3331 : :
3332 [ + + ]: 116642 : do isym=1,nsym
3333 : : ! Check that the symmetry operation preserves the wavevector
3334 : : ! (a translation is NOT allowed)
3335 : : if(symq(4,1,isym)==1 .and.&
3336 : : & symq(1,1,isym)==0 .and.&
3337 [ + + + + : 113385 : & symq(2,1,isym)==0 .and.&
+ + + + ]
3338 : 3257 : & symq(3,1,isym)==0 )then
3339 : : ! Check that the symmetry operation preserves the atom
3340 [ + + ]: 29871 : if(ipert==indsym(4,isym,ipert))then
3341 : : ! Check if the direction is preserved
3342 : : tok=1
3343 [ + + ]: 111740 : do idir1=1,3
3344 [ + + + + : 83805 : if((idir1==idir.and.symrec(idir,idir1,isym)/=1) .or.&
+ + + + ]
3345 : 27935 : & (idir1/=idir.and.symrec(idir,idir1,isym)/=0))then
3346 : 83805 : tok=0
3347 : : end if
3348 : : end do
3349 [ + + ]: 27935 : if(tok==1)then
3350 : : ! All the tests have been successful !
3351 : 5649 : nsym1=nsym1+1
3352 : 5649 : symaf1(nsym1)=symafm(isym)
3353 [ + + ]: 22596 : do ii=1,3
3354 : 16947 : tnons1(ii,nsym1)=tnons(ii,isym)
3355 [ + + ]: 73437 : do jj=1,3
3356 : 67788 : symrl1(ii,jj,nsym1)=symrel(ii,jj,isym)
3357 : : end do
3358 : : end do
3359 : : end if
3360 : :
3361 : : end if
3362 : : end if
3363 : : end do
3364 : : end if
3365 : :
3366 [ - + ]: 4129 : if (nsym1<1) then
3367 : 0 : write(msg,'(a,i0,a)')' The number of selected symmetries should be > 0, while it is nsym= ',nsym1,'.'
3368 : 0 : ABI_BUG(msg)
3369 : : end if
3370 : :
3371 [ + + ]: 4129 : if (nsym1 /= 1) then
3372 [ + + + - ]: 1586 : if (iout /= ount .and. iout > 0) then
3373 : 704 : write(msg,'(a,i5,a)')' Found ',nsym1,' symmetries that leave the perturbation invariant.'
3374 : 704 : call wrtout(iout,msg)
3375 : : end if
3376 : 1586 : write(msg,'(a,i5,a)')' littlegroup_pert: found ',nsym1,' symmetries that leave the perturbation invariant: '
3377 : 1586 : call wrtout(ount,msg)
3378 : : else
3379 [ + + + - ]: 2543 : if (iout /= ount .and. iout > 0) then
3380 : 1577 : write(msg,'(a,a)')' The set of symmetries contains',' only one element for this perturbation.'
3381 : 1577 : call wrtout(iout,msg)
3382 : : end if
3383 : 2543 : write(msg,'(a)')' littlegroup_pert: only one element in the set of symmetries for this perturbation:'
3384 : 2543 : call wrtout(ount,msg)
3385 : : end if
3386 : :
3387 [ + + ]: 4129 : if (ount > 0) then
3388 [ + + ]: 6812 : do isym=1,nsym1
3389 : 4401 : write(msg, '(9i4)' )((symrl1(ii,jj,isym),ii=1,3),jj=1,3)
3390 : 6812 : call wrtout(ount,msg)
3391 : : end do
3392 : : end if
3393 : :
3394 : 4129 : end subroutine littlegroup_pert
3395 : : !!***
3396 : :
3397 : : !!****f* ABINIT/irreducible_set_pert
3398 : : !! NAME
3399 : : !! irreducible_set_pert
3400 : : !!
3401 : : !! FUNCTION
3402 : : !! Determines a set of perturbations that form a basis
3403 : : !! in that, using symmetry, they can be used to generate
3404 : : !! all other perturbations that are asked to be calculated (target).
3405 : : !!
3406 : : !! INPUTS
3407 : : !! indsym(4,nsym,natom)=indirect indexing array described above: for each
3408 : : !! isym,iatom, fourth element is label of atom into which iatom is sent by
3409 : : !! INVERSE of symmetry operation isym; first three elements are the primitive
3410 : : !! translations which must be subtracted after the transformation to get back
3411 : : !! to the original unit cell.
3412 : : !! mpert =maximum number of iper
3413 : : !! natom= number of atoms
3414 : : !! nsym=number of space group symmetries
3415 : : !! rfdir(3)=direction for the perturbations
3416 : : !! rfpert(mpert)=information on the perturbations
3417 : : !! symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
3418 : : !! symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
3419 : : !! symq(4,2,nsym)= (integer) three first numbers define the G vector;
3420 : : !! fourth number is 0 if the q-vector is not preserved, is 1 otherwise
3421 : : !! second index is one without time-reversal symmetry, two with time-reversal symmetry
3422 : : !!
3423 : : !! OUTPUT
3424 : : !! pertsy(3,mpert)= the target perturbation is described by the two last indices (idir, and ipert),
3425 : : !! the value is 0, 1 or -1, see notes.
3426 : : !!
3427 : : !! NOTES
3428 : : !! Output will be in the pertsy array,
3429 : : !! 0 for non-target perturbations
3430 : : !! 1 for basis perturbations
3431 : : !! -1 for perturbations that can be found from basis perturbations
3432 : : !!
3433 : : !! SOURCE
3434 : :
3435 : 1523 : subroutine irreducible_set_pert(indsym,mpert,natom,nsym,pertsy,rfdir,rfpert,symq,symrec,symrel)
3436 : :
3437 : : !Arguments -------------------------------
3438 : : !scalars
3439 : : integer,intent(in) :: mpert,natom,nsym
3440 : : !arrays
3441 : : integer,intent(in) :: indsym(4,nsym,natom),rfdir(3),rfpert(mpert)
3442 : : integer,intent(in) :: symq(4,2,nsym),symrec(3,3,nsym),symrel(3,3,nsym)
3443 : : integer,intent(out) :: pertsy(3,mpert)
3444 : :
3445 : : !Local variables -------------------------
3446 : : !scalars
3447 : : integer :: found,idir1,idisy1,ii,ipert1,ipesy1,isign,isym,itirev,jj
3448 : : !arrays
3449 : : integer :: sym1(3,3)
3450 : : ! *********************************************************************
3451 : :
3452 : : !Zero pertsy
3453 [ + + + + ]: 53899 : pertsy(:,:)=0
3454 : :
3455 [ + + ]: 14617 : do ipert1=1,mpert
3456 [ + + ]: 53899 : do idir1=1,3
3457 [ + + + + ]: 52376 : if(rfpert(ipert1)==1.and.rfdir(idir1)==1)then
3458 : : ! write(std_out,*)' for candidate idir =',idir1,' ipert = ',ipert1
3459 : :
3460 : : ! Loop on all symmetries, including time-reversal
3461 [ + + ]: 255631 : do isym=1,nsym
3462 [ + + ]: 674042 : do itirev=1,2
3463 : 458776 : isign=3-2*itirev
3464 : :
3465 [ + + ]: 666454 : if(symq(4,itirev,isym)/=0)then
3466 : :
3467 : 182588 : found=1
3468 : :
3469 : : ! Here select the symmetric of ipert1
3470 [ + + ]: 182588 : if(ipert1<=natom)then
3471 : 137899 : ipesy1=indsym(4,isym,ipert1)
3472 [ + + ]: 551596 : do ii=1,3
3473 [ + + ]: 1792687 : do jj=1,3
3474 : 1654788 : sym1(ii,jj)=symrec(ii,jj,isym)
3475 : : end do
3476 : : end do
3477 [ + + ]: 44689 : else if(ipert1==(natom+2))then
3478 : 50136 : ipesy1=ipert1
3479 [ + + ]: 50136 : do ii=1,3
3480 [ + + ]: 162942 : do jj=1,3
3481 : 150408 : sym1(ii,jj)=symrel(ii,jj,isym)
3482 : : end do
3483 : : end do
3484 : : else
3485 : : found=0
3486 : : end if
3487 : :
3488 : : ! Now that a symmetric perturbation has been obtained,
3489 : : ! including the expression of the symmetry matrix, see
3490 : : ! if the symmetric perturbations are available
3491 : : if( found==1 ) then
3492 : :
3493 [ + + ]: 385804 : do idisy1=1,3
3494 [ + + ]: 385804 : if(sym1(idir1,idisy1)/=0)then
3495 [ + + ]: 186141 : if(pertsy(idisy1,ipesy1)==0)then
3496 : : found=0
3497 : : exit
3498 : : end if
3499 : : end if
3500 : : end do
3501 : : end if
3502 : :
3503 : : ! Now, if still found, then it is a symmetric
3504 : : ! of some linear combination of existing perturbations
3505 [ + + ]: 182588 : if(found==1)then
3506 : :
3507 : : ! DEBUG
3508 : : ! write(std_out,*)' all found ! isym, isign= ',isym,isign
3509 : : ! write(std_out,1010)((sym1(ii,jj),ii=1,3),jj=1,3)
3510 : : ! write(std_out,1010)((sym2(ii,jj),ii=1,3),jj=1,3)
3511 : : ! write(std_out,*)sumr,sumi
3512 : : ! 1010 format(9i4)
3513 : : ! ENDDEBUG
3514 : :
3515 : 40365 : pertsy(idir1,ipert1)=-1
3516 : 40365 : exit ! Exit loop on symmetry operations
3517 : :
3518 : : end if
3519 : :
3520 : : end if ! End loop on all symmetries + time-reversal
3521 : : end do
3522 : : end do
3523 : :
3524 : : ! Now that all symmetries have been examined,
3525 : : ! if still not symmetric of a linear combination
3526 : : ! of basis perturbations, then it is a basis perturbation
3527 [ + + ]: 7588 : if(pertsy(idir1,ipert1)/=-1) pertsy(idir1,ipert1)=1
3528 : : ! write(std_out,'(a,3i5)' ) ' irreducible_set_pert :',idir1,ipert1,pertsy(idir1,ipert1)
3529 : :
3530 : : end if ! End big loop on all elements
3531 : : end do
3532 : : end do
3533 : :
3534 : 1523 : end subroutine irreducible_set_pert
3535 : : !!***
3536 : :
3537 : : !!****f* m_dynmat/d3lwsym
3538 : : !! NAME
3539 : : !! d3lwsym
3540 : : !!
3541 : : !! FUNCTION
3542 : : !! Given a set of calculated elements of the 3DTE matrix,
3543 : : !! build (nearly) all the other matrix elements that can be build using symmetries.
3544 : : !!
3545 : : !! INPUTS
3546 : : !! has_strain = if .true. i2pert includes strain perturbation
3547 : : !! indsym(4,nsym,natom)=indirect indexing array : for each
3548 : : !! isym,iatom, fourth element is label of atom into which iatom is sent by
3549 : : !! INVERSE of symmetry operation isym; first three elements are the primitive
3550 : : !! translations which must be subtracted after the transformation to get back
3551 : : !! to the original unit cell.
3552 : : !! mpert =maximum number of ipert
3553 : : !! natom= number of atoms
3554 : : !! nsym=number of space group symmetries
3555 : : !! symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal reduced space)
3556 : : !! symrel(3,3,nsym)=3x3 matrices of the group symmetries (real reduced space)
3557 : : !! symrel_cart(3,3,nsym)=3x3 matrices of the group symmetries (real cartesian space)
3558 : : !!
3559 : : !! SIDE EFFECTS
3560 : : !! Input/Output
3561 : : !! blkflg(3,mpert,3,mpert,3,mpert)= matrix that indicates if an
3562 : : !! element of d3 is available (1 if available, 0 otherwise)
3563 : : !! d3(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE
3564 : : !!
3565 : : !! SOURCE
3566 : :
3567 : : !subroutine d3lwsym(blkflg,d3,has_strain,indsym,mpert,natom,nsym,symrec,symrel,symrel_cart)
3568 : 17 : subroutine d3lwsym(blkflg,d3,indsym,mpert,natom,nsym,symrec,symrel)
3569 : :
3570 : : !Arguments -------------------------------
3571 : : !scalars
3572 : : integer,intent(in) :: mpert,natom,nsym
3573 : : ! logical,intent(in) :: has_strain
3574 : : !arrays
3575 : : integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym)
3576 : : integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert)
3577 : : real(dp),intent(inout) :: d3(2,3,mpert,3,mpert,3,mpert)
3578 : : ! real(dp),intent(in) :: symrel_cart(3,3,nsym)
3579 : :
3580 : : !Local variables -------------------------
3581 : : !scalars
3582 : : integer :: found,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,idisy1,idisy2,idisy3
3583 : : integer :: ipesy1,ipesy2,ipesy3,isym,ithree
3584 : : !integer :: istr,i2dir_a,i2dir_b,disy2_a,idisy2_b
3585 : : real(dp) :: sumi,sumr
3586 : : logical :: is_strain
3587 : : !arrays
3588 : : ! integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/)
3589 : : integer :: sym1(3,3),sym2(3,3),sym3(3,3)
3590 : : ! integer :: strflg(3,mpert,3,3,3,mpert),strflg_car(3,mpert,3,3,3,mpert)
3591 : : ! real(dp) :: d3str(2,3,mpert,3,3,3,mpert)
3592 : : ! *********************************************************************
3593 : :
3594 : : !First, take into account the permutations symmetry of
3595 : : !(i1pert,i1dir) and (i2pert,i2dir)
3596 [ + + ]: 259 : do i1pert = 1, mpert
3597 [ + + ]: 4183 : do i2pert = 1, mpert
3598 [ + + ]: 75694 : do i3pert = 1, mpert
3599 : :
3600 [ + + ]: 290036 : do i1dir = 1, 3
3601 [ + + ]: 929864 : do i2dir = 1, 3
3602 [ + + ]: 2789592 : do i3dir = 1, 3
3603 : :
3604 [ + + + + ]: 1931256 : if ((blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1).and. &
3605 : 643752 : (blkflg(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert)/=1)) then
3606 : :
3607 : : d3(1,i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = &
3608 : 842 : d3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
3609 : : d3(2,i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = &
3610 : 842 : -d3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
3611 : :
3612 : 842 : blkflg(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = 1
3613 : :
3614 : : end if
3615 : :
3616 : : end do
3617 : : end do
3618 : : end do
3619 : :
3620 : : end do
3621 : : end do
3622 : : end do
3623 : :
3624 : : !For strain perturbation we need an array with the two strain indexes
3625 : : ! if (has_strain) then
3626 : : ! strflg=0
3627 : : ! d3str=zero
3628 : : ! do i3pert=1, mpert
3629 : : ! do i3dir=1,3
3630 : : ! do i2pert=natom+3,natom+4
3631 : : ! do i2dir=1,3
3632 : : ! if (i2pert==natom+3) istr=i2dir
3633 : : ! if (i2pert==natom+4) istr=3+i2dir
3634 : : ! i2dir_a=idx(2*istr-1); i2dir_b=idx(2*istr)
3635 : : ! do i1pert=1,mpert
3636 : : ! do i1dir=1,3
3637 : : ! if (blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then
3638 : : ! strflg(i1dir,i1pert,i2dir_a,i2dir_b,i3dir,i3pert)=1
3639 : : ! d3str(:,i1dir,i1pert,i2dir_a,i2dir_b,i3dir,i3pert)= &
3640 : : ! & d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
3641 : : ! if (i2pert==natom+4) then
3642 : : ! strflg(i1dir,i1pert,i2dir_b,i2dir_a,i3dir,i3pert)=1
3643 : : ! d3str(:,i1dir,i1pert,i2dir_b,i2dir_a,i3dir,i3pert)= &
3644 : : ! & d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)
3645 : : ! end if
3646 : : ! end if
3647 : : ! end do
3648 : : ! end do
3649 : : ! end do
3650 : : ! end do
3651 : : ! end do
3652 : : ! end do
3653 : : ! end if
3654 : :
3655 : : !Big Big Loop : symmetrize three times, because
3656 : : !of some cases in which one element is not yet available
3657 : : !at the first pass, and even at the second one !
3658 : :
3659 [ + + ]: 68 : do ithree=1,3
3660 : :
3661 : : ! Loop over perturbations
3662 [ + + ]: 794 : do i1pert = 1, mpert
3663 [ + + ]: 12549 : do i2pert = 1, mpert
3664 : : is_strain=.false.
3665 [ + + ]: 227082 : do i3pert = 1, mpert
3666 : :
3667 [ + + ]: 870108 : do i1dir = 1, 3
3668 [ + + ]: 2789592 : do i2dir = 1, 3
3669 [ + + ]: 8368776 : do i3dir = 1, 3
3670 : :
3671 : : ! Will get element (idir1,ipert1,idir2,ipert2)
3672 : : ! so this element should not yet be present ...
3673 [ + + ]: 7725024 : if(blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=1)then
3674 : :
3675 [ + + ]: 17360586 : d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 0_dp
3676 : :
3677 [ + + ]: 36156726 : do isym = 1, nsym
3678 : :
3679 : 30369864 : found = 1
3680 : :
3681 [ + + ]: 30369864 : if (i1pert <= natom) then
3682 : 7337334 : ipesy1 = indsym(4,isym,i1pert)
3683 [ + + + + ]: 95385342 : sym1(:,:) = symrec(:,:,isym)
3684 [ + + ]: 23032530 : else if (i1pert == natom + 2) then
3685 : 2634504 : ipesy1 = i1pert
3686 [ + + + + ]: 34248552 : sym1(:,:) = symrel(:,:,isym)
3687 : : else
3688 : : found = 0
3689 : : end if
3690 : :
3691 [ + + ]: 30369864 : if (i2pert <= natom) then
3692 : 7337334 : ipesy2 = indsym(4,isym,i2pert)
3693 [ + + + + ]: 95385342 : sym2(:,:) = symrec(:,:,isym)
3694 [ + + ]: 23032530 : else if (i2pert == natom + 2) then
3695 : 2634504 : ipesy2 = i2pert
3696 [ + + + + ]: 34248552 : sym2(:,:) = symrel(:,:,isym)
3697 [ + + + + ]: 20398026 : else if (i2pert == natom + 3.or. i2pert == natom + 4) then
3698 : : !TODO: Symmetries on strain perturbation do not work yet.
3699 : : found = 0
3700 : : is_strain=.true.
3701 : :
3702 : : ! ipesy2 = i2pert
3703 : : ! sym2(:,:) = NINT(symrel_cart(:,:,isym))
3704 : : ! if (i2pert==natom+3) istr=i2dir
3705 : : ! if (i2pert==natom+4) istr=3+i2dir
3706 : : ! i2dir_a=idx(2*istr-1); i2dir_b=idx(2*istr)
3707 : : else
3708 : 15071142 : found = 0
3709 : : end if
3710 : :
3711 [ + + ]: 30369864 : if (i3pert <= natom) then
3712 : 7457832 : ipesy3 = indsym(4,isym,i3pert)
3713 [ + + + + ]: 96951816 : sym3(:,:) = symrec(:,:,isym)
3714 [ + + + + ]: 22912032 : else if (i3pert == natom + 2.or.i3pert == natom + 8) then
3715 : 5142504 : ipesy3 = i3pert
3716 [ + + + + ]: 66852552 : sym3(:,:) = symrel(:,:,isym)
3717 : : else
3718 : : found = 0
3719 : : end if
3720 : :
3721 : 30369864 : sumr = 0_dp ; sumi = 0_dp;
3722 [ + + ]: 30369864 : if (.not.is_strain) then
3723 [ + + ]: 100171920 : do idisy1 = 1, 3
3724 [ + + ]: 325558740 : do idisy2 = 1, 3
3725 [ + + ]: 976676220 : do idisy3 = 1, 3
3726 : :
3727 : : if ((sym1(i1dir,idisy1) /=0).and.(sym2(i2dir,idisy2) /=0) &
3728 [ + + + + : 901547280 : & .and.(sym3(i3dir,idisy3) /=0)) then
+ + ]
3729 : :
3730 [ + + ]: 81492718 : if (blkflg(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 1) then
3731 : :
3732 : : sumr = sumr + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*&
3733 : 2909318 : & sym3(i3dir,idisy3)*d3(1,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3)
3734 : : sumi = sumi + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*&
3735 : 2909318 : & sym3(i3dir,idisy3)*d3(2,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3)
3736 : :
3737 : : else
3738 : :
3739 : : found = 0
3740 : :
3741 : : end if
3742 : :
3743 : : end if
3744 : :
3745 : : end do
3746 : : end do
3747 : : end do
3748 : : else
3749 : : ! do idisy1 = 1, 3
3750 : : ! !do idisy2_a = 1, 3
3751 : : ! ! do idisy2_b = 1, 3
3752 : : ! do idisy2 = 1, 3
3753 : : ! if (ipesy2==natom+3) istr=idisy2
3754 : : ! if (ipesy2==natom+4) istr=3+idisy2
3755 : : ! idisy2_a=idx(2*istr-1); idisy2_b=idx(2*istr)
3756 : : ! do idisy3 = 1, 3
3757 : : !
3758 : : ! if ((sym1(i1dir,idisy1) /=0).and.(sym2(i2dir_a,idisy2_a) /=0) &
3759 : : !& .and.(sym2(i2dir_b,idisy2_b) /=0).and.(sym3(i3dir,idisy3) /=0)) then
3760 : : !
3761 : : ! if (strflg(idisy1,ipesy1,idisy2_a,idisy2_b,idisy3,ipesy3) == 1) then
3762 : : !
3763 : : ! sumr = sumr + sym1(i1dir,idisy1)*sym2(i2dir_a,idisy2_a)* &
3764 : : !& sym2(i2dir_b,idisy2_b)*sym3(i3dir,idisy3)*&
3765 : : !& d3str(1,idisy1,ipesy1,idisy2_a,idisy2_b,idisy3,ipesy3)
3766 : : ! sumi = sumi + sym1(i1dir,idisy1)*sym2(i2dir_a,idisy2_b)*&
3767 : : !& sym2(i2dir_b,idisy2_b)*sym3(i3dir,idisy3)*&
3768 : : !& d3str(2,idisy1,ipesy1,idisy2_a,idisy2_b,idisy3,ipesy3)
3769 : : !
3770 : : ! else
3771 : : !
3772 : : ! found = 0
3773 : : !
3774 : : ! end if
3775 : : !
3776 : : ! end if
3777 : : !
3778 : : ! end do
3779 : : ! !end do
3780 : : ! end do
3781 : : ! end do
3782 : : end if
3783 : :
3784 [ + + ]: 36156726 : if (found == 1) then
3785 : 14544 : d3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr
3786 : 14544 : d3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi
3787 : 14544 : blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
3788 : : end if
3789 : :
3790 : : end do ! isym
3791 : :
3792 : : end if ! blkflg
3793 : :
3794 : : ! Close loop over perturbations
3795 : : end do
3796 : : end do
3797 : : end do
3798 : : end do
3799 : : end do
3800 : : end do
3801 : :
3802 : : end do ! close loop over ithree
3803 : :
3804 : 17 : end subroutine d3lwsym
3805 : : !!***
3806 : :
3807 : : !!****f* m_dynmat/sylwtens
3808 : : !!
3809 : : !! NAME
3810 : : !! sylwtens
3811 : : !!
3812 : : !! FUNCTION
3813 : : !! Determines the set of irreductible elements of the non-linear
3814 : : !! optical susceptibility and Raman tensors
3815 : : !!
3816 : : !! INPUTS
3817 : : !! has_strain = if .true. i2pert includes strain perturbation
3818 : : !! indsym(4,nsym,natom)=indirect indexing array described above: for each
3819 : : !! isym,iatom, fourth element is label of atom into which iatom is sent by
3820 : : !! INVERSE of symmetry operation isym; first three elements are the primitive
3821 : : !! translations which must be subtracted after the transformation to get back
3822 : : !! to the original unit cell.
3823 : : !! mpert =maximum number of ipert
3824 : : !! natom= number of atoms
3825 : : !! nsym=number of space group symmetries
3826 : : !! symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal reduced space)
3827 : : !! symrel(3,3,nsym)=3x3 matrices of the group symmetries (real reduced space)
3828 : : !! symrel_cart(3,3,nsym)=3x3 matrices of the group symmetries (real cartesian space)
3829 : : !!
3830 : : !! OUTPUT
3831 : : !! (see side effects)
3832 : : !!
3833 : : !! SIDE EFFECTS
3834 : : !! rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations
3835 : : !! that have to be computed
3836 : : !! At the input :
3837 : : !! 1 -> element has to be computed explicitly
3838 : : !! At the output :
3839 : : !! 1 -> element has to be computed explicitly
3840 : : !! -1 -> use symmetry operations to obtain the corresponding element
3841 : : !! -2 -> element is zero by symmetry
3842 : : !!
3843 : : !! SOURCE
3844 : :
3845 : 33 : subroutine sylwtens(indsym,mpert,natom,nsym,rfpert,symrec,symrel)
3846 : :
3847 : : !Arguments -------------------------------
3848 : : !scalars
3849 : : integer,intent(in) :: mpert,natom,nsym
3850 : : !arrays
3851 : : integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym)
3852 : : integer,intent(inout) :: rfpert(3,mpert,3,mpert,3,mpert)
3853 : : ! real(dp),intent(in) :: symrel_cart(3,3,nsym)
3854 : :
3855 : : !Local variables -------------------------
3856 : : !scalars
3857 : : integer :: flag,found,i1dir,i1dir_,i1pert,i1pert_,i2dir,i2dir_,i2pert,i2pert_
3858 : : ! integer :: i2dir_a,i2dir_b
3859 : : integer :: i3dir,i3dir_,i3pert,i3pert_,idisy1,idisy2,idisy3,ipesy1,ipesy2
3860 : : integer :: ipesy3,isym
3861 : : ! integer :: istr,idisy2_a,idisy2_b
3862 : : logical :: is_strain, is_timdisp
3863 : : ! real(dp) :: flag_dp
3864 : : !arrays
3865 : : ! integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/)
3866 : : integer :: sym1(3,3),sym2(3,3),sym3(3,3)
3867 : 33 : integer,allocatable :: pertsy(:,:,:,:,:,:)
3868 : : !***********************************************************************
3869 : :
3870 [ - + + - : 231 : ABI_MALLOC(pertsy,(3,mpert,3,mpert,3,mpert))
+ - - + +
- + - - +
+ - + - +
- - + +
- ]
3871 [ + + + + : 3406033 : pertsy(:,:,:,:,:,:) = 0
+ + + + +
+ + + ]
3872 : :
3873 : : !Loop over perturbations
3874 : :
3875 [ + + ]: 445 : do i1pert_ = 1, mpert
3876 [ + + ]: 6227 : do i2pert_ = 1, mpert
3877 : : is_strain=.false.
3878 [ + + ]: 98832 : do i3pert_ = 1, mpert
3879 : : is_timdisp=.false.
3880 : :
3881 [ + + ]: 376334 : do i1dir_ = 1, 3
3882 [ + + ]: 1204294 : do i2dir_ = 1, 3
3883 [ + + ]: 3612882 : do i3dir_ = 1, 3
3884 : :
3885 : 2501226 : i1pert = (mpert - i1pert_ + 1)
3886 [ + + ]: 2501226 : if (i1pert <= natom) i1pert = natom + 1 - i1pert
3887 : 2501226 : i2pert = (mpert - i2pert_ + 1)
3888 [ + + ]: 2501226 : if (i2pert <= natom) i2pert = natom + 1 - i2pert
3889 : 2501226 : i3pert = (mpert - i3pert_ + 1)
3890 [ + + ]: 2501226 : if (i3pert <= natom) i3pert = natom + 1 - i3pert
3891 : :
3892 [ + + ]: 2501226 : if (i1pert <= natom) then
3893 : : i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_
3894 [ + + ]: 1775250 : else if (i2pert <= natom) then
3895 : : i1dir = i2dir_ ; i2dir = i1dir_ ; i3dir = i3dir_
3896 [ + + ]: 1284606 : else if (i3pert <= natom) then
3897 : : i1dir = i3dir_ ; i2dir = i2dir_ ; i3dir = i1dir_
3898 : : else
3899 : 1669194 : i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_
3900 : : end if
3901 : :
3902 [ + + ]: 3334968 : if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) /= 0) then
3903 : :
3904 : : ! Loop over all symmetries
3905 : :
3906 : : flag = 0
3907 [ + + ]: 131904 : do isym = 1, nsym
3908 : :
3909 : 129012 : found = 1
3910 : :
3911 : : ! Select the symmetric element of i1pert,i2pert,i3pert
3912 : :
3913 [ + + ]: 129012 : if (i1pert <= natom) then
3914 : 62218 : ipesy1 = indsym(4,isym,i1pert)
3915 [ + + + + ]: 808834 : sym1(:,:) = symrec(:,:,isym)
3916 [ + + ]: 66794 : else if (i1pert == natom + 2) then
3917 : 51230 : ipesy1 = i1pert
3918 [ + + + + ]: 665990 : sym1(:,:) = symrel(:,:,isym)
3919 : : else
3920 : : found = 0
3921 : : end if
3922 : :
3923 [ + + ]: 129012 : if (i2pert <= natom) then
3924 : 74302 : ipesy2 = indsym(4,isym,i2pert)
3925 [ + + + + ]: 965926 : sym2(:,:) = symrec(:,:,isym)
3926 [ + + ]: 54710 : else if (i2pert == natom + 2) then
3927 : 266 : ipesy2 = i2pert
3928 [ + + + + ]: 3458 : sym2(:,:) = symrel(:,:,isym)
3929 [ + + + + ]: 54444 : else if (i2pert == natom + 3.or. i2pert == natom + 4) then
3930 : : ! !TODO: Symmetries on strain perturbation do not work yet.
3931 : : found = 0
3932 : : is_strain=.true.
3933 : : !
3934 : : ! ipesy2 = i2pert
3935 : : ! sym2(:,:) = NINT(symrel_cart(:,:,isym))
3936 : : ! if (i2pert==natom+3) istr=i2dir
3937 : : ! if (i2pert==natom+4) istr=3+i2dir
3938 : : ! i2dir_a=idx(2*istr-1); i2dir_b=idx(2*istr)
3939 : : else
3940 : 12 : found = 0
3941 : : end if
3942 : :
3943 [ + + ]: 129012 : if (i3pert == natom + 8) then
3944 : 128984 : ipesy3 = i3pert
3945 [ + + + + ]: 1676792 : sym3(:,:) = symrel(:,:,isym)
3946 [ - + ]: 28 : else if (i3pert == natom + 9) then
3947 : : is_timdisp=.true.
3948 : : found = 0
3949 : : else
3950 : : found = 0
3951 : : end if
3952 : :
3953 : :
3954 : : ! See if the symmetric element is available and check if some
3955 : : ! of the elements may be zero. In the latter case, they do not need
3956 : : ! to be computed.
3957 : :
3958 [ + + ]: 128984 : if (.not.is_timdisp.and..not.is_strain) then
3959 : : if ((flag /= -1).and.&
3960 [ + + + + ]: 74552 : & (ipesy1==i1pert).and.(ipesy2==i2pert).and.(ipesy3==i3pert)) then
3961 : 37368 : flag = sym1(i1dir,i1dir)*sym2(i2dir,i2dir)*sym3(i3dir,i3dir)
3962 : : end if
3963 : :
3964 [ + + ]: 298208 : do idisy1 = 1, 3
3965 [ + + ]: 969176 : do idisy2 = 1, 3
3966 [ + + ]: 2907528 : do idisy3 = 1, 3
3967 : :
3968 [ + + + + : 2012904 : if ((sym1(i1dir,idisy1) /= 0).and.(sym2(i2dir,idisy2) /= 0).and.&
+ + ]
3969 : : & (sym3(i3dir,idisy3) /= 0)) then
3970 [ + + ]: 262548 : if (pertsy(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 0) then
3971 : 2012904 : found = 0
3972 : : ! exit ! exit loop over symmetries
3973 : : end if
3974 : : end if
3975 : :
3976 : :
3977 [ + + + + : 2012904 : if ((flag == -1).and.&
+ + ]
3978 : 670968 : & ((idisy1/=i1dir).or.(idisy2/=i2dir).or.(idisy3/=i3dir))) then
3979 [ + + + + : 42972 : if ((sym1(i1dir,idisy1)/=0).and.(sym2(i2dir,idisy2)/=0).and.&
+ + ]
3980 : : & (sym3(i3dir,idisy3)/=0)) then
3981 : 2012904 : flag = 0
3982 : : end if
3983 : : end if
3984 : :
3985 : : end do
3986 : : end do
3987 : : end do
3988 : : ! else
3989 : : ! if ((flag_dp /= -1).and.&
3990 : : !& (ipesy1==i1pert).and.(ipesy2==i2pert).and.(ipesy3==i3pert)) then
3991 : : ! flag = sym1(i1dir,i1dir)*sym2(i2dir_a,i2dir_a)* &
3992 : : ! & sym2(i2dir_b,i2dir_b)*sym3(i3dir,i3dir)
3993 : : ! end if
3994 : : !
3995 : : ! do idisy1 = 1, 3
3996 : : ! do idisy2 = 1, 3
3997 : : ! if (ipesy2==natom+3) istr=idisy2
3998 : : ! if (ipesy2==natom+4) istr=3+idisy2
3999 : : ! idisy2_a=idx(2*istr-1); idisy2_b=idx(2*istr)
4000 : : ! do idisy3 = 1, 3
4001 : : !
4002 : : ! if ((sym1(i1dir,idisy1) /= 0).and.(sym2(i2dir_a,idisy2_a) /= 0).and.&
4003 : : !& (sym2(i2dir_b,idisy2_b) /= 0).and.(sym3(i3dir,idisy3) /= 0)) then
4004 : : ! if (pertsy(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 0) then
4005 : : ! found = 0
4006 : : !! exit ! exit loop over symmetries
4007 : : ! end if
4008 : : ! end if
4009 : : !
4010 : : !
4011 : : ! if ((flag == -1).and.&
4012 : : !& ((idisy1/=i1dir).or.(idisy2_a/=i2dir_a).or.(idisy2_b/=i2dir_b).or.(idisy3/=i3dir))) then
4013 : : ! if ((sym1(i1dir,idisy1)/=0).and.(sym2(i2dir_a,idisy2_a)/=0).and.&
4014 : : !& (sym2(i2dir_b,idisy2_b)/=0).and.(sym3(i3dir,idisy3)/=0)) then
4015 : : ! flag = 0
4016 : : ! end if
4017 : : ! end if
4018 : : !
4019 : : ! end do
4020 : : ! end do
4021 : : ! end do
4022 : : end if
4023 : :
4024 [ + + ]: 129012 : if (found == 1) then
4025 : 28177 : pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -1
4026 : : end if
4027 : :
4028 : : ! In case a symmetry operation only changes the sign of an
4029 : : ! element, this element has to be equal to zero
4030 : :
4031 [ + + ]: 131904 : if (flag == -1) then
4032 : 646 : pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -2
4033 : 646 : exit
4034 : : end if
4035 : :
4036 : : end do ! close loop on symmetries
4037 : :
4038 : : ! If the element i1pert,i2pert,i3pert is not symmetric
4039 : : ! to a basis element, it is a basis element
4040 : :
4041 [ + + ]: 3538 : if (pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) > -1) then
4042 : 1551 : pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1
4043 : : end if
4044 : :
4045 : : end if ! rfpert /= 0
4046 : :
4047 : : end do ! close loop over perturbations
4048 : : end do
4049 : : end do
4050 : : end do
4051 : : end do
4052 : : end do
4053 : :
4054 : : !Now, take into account the permutation of (i1pert,i1dir)
4055 : : !and (i2pert,i2dir)
4056 : :
4057 : :
4058 [ + + ]: 445 : do i1pert = 1, mpert
4059 [ + + ]: 6227 : do i2pert = 1, mpert
4060 [ + + ]: 98832 : do i3pert = 1, mpert
4061 : :
4062 [ + + ]: 376334 : do i1dir = 1, 3
4063 [ + + ]: 1204294 : do i2dir = 1, 3
4064 [ + + ]: 3612882 : do i3dir = 1, 3
4065 : :
4066 [ + + ]: 3334968 : if ((i1pert /= i2pert).or.(i1dir /= i2dir)) then
4067 : :
4068 [ + + + + ]: 2449188 : if ((pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) == 1).and.&
4069 : : (pertsy(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) == 1)) then
4070 : 10 : pertsy(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = -1
4071 : : end if
4072 : :
4073 : : end if
4074 : :
4075 : : end do
4076 : : end do
4077 : : end do
4078 : :
4079 : : end do
4080 : : end do
4081 : : end do
4082 : :
4083 [ + + + + : 3406033 : rfpert(:,:,:,:,:,:) = pertsy(:,:,:,:,:,:)
+ + + + +
+ + + ]
4084 : :
4085 [ - + ]: 33 : ABI_FREE(pertsy)
4086 : :
4087 : 33 : end subroutine sylwtens
4088 : : !!***
4089 : : !!****f* m_geometry/vcart2ylm
4090 : : !! NAME
4091 : : !! vcart2ylm
4092 : : !!
4093 : : !! FUNCTION
4094 : : !! Convert a 3D Cartesian vector into spherical coordinates (r, theta, phi)
4095 : : !! suitable for spherical harmonics calculations (Y_lm).
4096 : : !! Angles are returned in degrees.
4097 : : !!
4098 : : !! INPUTS
4099 : : !! vector(3) = Cartesian vector (x, y, z)
4100 : : !!
4101 : : !! OUTPUTS
4102 : : !! length = radial distance r = sqrt(x^2+y^2+z^2)
4103 : : !! theta = polar angle (from +z axis), in degrees
4104 : : !! phi = azimuthal angle (from +x axis in xy-plane), in degrees
4105 : : !!
4106 : : !! NOTES
4107 : : !! - If the vector magnitude is very small (<1e-9), theta and phi are set to 0.
4108 : : !! - Uses the physics/Y_lm convention: theta = polar, phi = azimuth.
4109 : : !! - phi is computed using atan2 to account for the correct quadrant.
4110 : : !! - This routine assumes input vector is real(8).
4111 : : !! - Can be easily extended to arrays of vectors.
4112 : :
4113 : 8900 : subroutine vcart2ylm(vector, length, theta, phi)
4114 : :
4115 : : !Arguments ---------------------------------------------
4116 : : !arrays
4117 : : real(8),intent(in) :: vector(3)
4118 : : !scalars
4119 : : real(8),intent(out) :: length, theta, phi
4120 : :
4121 : : ! Local
4122 : : real(8):: pi
4123 : :
4124 : 8900 : pi=4.0d0*datan(1.0d0)
4125 : : ! Compute spherical coordinates
4126 : 8900 : length = sqrt(vector(1)**2+vector(2)**2+vector(3)**2)
4127 : :
4128 [ + + ]: 8900 : if (length > tol6) then
4129 : 8318 : theta = acos(vector(3) / length) * 180.d0 / pi
4130 [ + + + + ]: 8318 : if (abs(vector(1)) > tol6 .or. abs(vector(2)) > tol6) then
4131 : 7495 : phi = atan2(vector(2), vector(1)) * 180.d0 / pi
4132 [ + + ]: 7495 : if (phi<-179) phi=180
4133 : : else
4134 : 823 : phi = 0.d0
4135 : : end if
4136 : : else
4137 : 582 : theta = 0.d0
4138 : 582 : phi = 0.d0
4139 : : end if
4140 : :
4141 : 8900 : end subroutine vcart2ylm
4142 : : !!***
4143 : :
4144 : : !!****f* m_geometry/cart2spinaxis
4145 : : !! NAME
4146 : : !! cart2spinaxis
4147 : : !!
4148 : : !! FUNCTION
4149 : : !! Compute the rotation matrix R = Rz(alpha)*Ry(beta) and rotate a vector in
4150 : : !! cartesian coordinate to spinaxis coordinates
4151 : : !!
4152 : : !! INPUTS
4153 : : !! alpha=Euler angle for rotation around z-axis
4154 : : !! beta=Euler angle for rotation around y-axis
4155 : : !! vin(3)=vector in the cartesian coordinate
4156 : : !!
4157 : : !! OUTPUT
4158 : : !! R(3,3)=rotation matrix from cartesian to spinaxis coordinates
4159 : : !! vout(3)=vector in spinaxis coordinate
4160 : : !!
4161 : : !! SOURCE
4162 : :
4163 : 2 : subroutine cart2spinaxis(alpha, beta, R, vin, vout)
4164 : :
4165 : : !Arguments -------------------------------
4166 : : !scalars
4167 : : real(dp),intent(in) :: alpha, beta
4168 : : !arrays
4169 : : real(dp),intent(out) :: R(3,3)
4170 : : real(dp),optional,intent(in) :: vin(3)
4171 : : real(dp),optional,intent(out) :: vout(3)
4172 : :
4173 : : !Local variables -------------------------
4174 : : !scalars
4175 : : real(dp) :: sb, cb, sa, ca
4176 : : !***********************************************************************
4177 : :
4178 : 2 : sb = sin(beta); cb = cos(beta)
4179 : 2 : sa = sin(alpha); ca = cos(alpha)
4180 : :
4181 : 2 : R(1,1) = cb*ca; R(2,1) = -sa; R(3,1) = sb*ca
4182 : 2 : R(1,2) = cb*sa; R(2,2) = ca; R(3,2) = sb*sa
4183 : 2 : R(1,3) = -sb; R(2,3) = zero; R(3,3) = cb
4184 : :
4185 [ - + ]: 2 : if (present(vin) .and. present(vout)) then
4186 [ # # # # ]: 0 : vout(:) = matmul(R, vin)
4187 : : end if
4188 : :
4189 : 2 : end subroutine cart2spinaxis
4190 : : !!***
4191 : :
4192 : 903 : end module m_geometry
4193 : : !!***
|
