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 1762519429 : 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 1762519429 : & +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 1759736480 : res=two_pi*SQRT(res)
140 : case default
141 1762519429 : ABI_BUG('Wrong value for space')
142 : end select
143 :
144 1762519429 : 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 1115268 : 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 1115268 : & +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 1115268 : res=two_pi*SQRT(res)
183 : case default
184 1115268 : ABI_BUG('Wrong value for space')
185 : end select
186 :
187 1115268 : 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 78114 : 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 78114 : +met(2,3)*(xv(2)*xw(3) + xv(3)*xw(2)) )
364 :
365 : select case (space)
366 : case ('r', 'R')
367 78114 : return
368 : case ('g', 'G')
369 78114 : res= res * (two_pi**2)
370 : case default
371 78114 : 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 5439078 : 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 5439078 : c(1) = a(2)*b(3) - a(3)*b(2)
404 5439078 : c(2) = -a(1)*b(3) + a(3)*b(1)
405 5439078 : c(3) = a(1)*b(2) - b(1)*a(2)
406 :
407 5439078 : 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 90 : 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 90 : if(present(normalize)) then
445 0 : nvec = normalize
446 : else
447 : nvec = .FALSE.
448 : end if
449 :
450 360 : do irprimd = 1, 3
451 1170 : do igprimd = 1, 3
452 810 : wedge(1,irprimd,igprimd) = rprimd(2,irprimd)*gprimd(3,igprimd) - rprimd(3,irprimd)*gprimd(2,igprimd)
453 810 : wedge(2,irprimd,igprimd) = rprimd(3,irprimd)*gprimd(1,igprimd) - rprimd(1,irprimd)*gprimd(3,igprimd)
454 1080 : wedge(3,irprimd,igprimd) = rprimd(1,irprimd)*gprimd(2,igprimd) - rprimd(2,irprimd)*gprimd(1,igprimd)
455 : end do
456 : end do
457 :
458 90 : 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 90 : 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 5420008 : 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 5420008 : produv(:) = zero
502 21680032 : do irprimd = 1, 3
503 70460104 : do igprimd = 1, 3
504 211380312 : produv(1:3) = produv(1:3) + u(irprimd)*v(igprimd)*wedgebasis(1:3,irprimd,igprimd)
505 : end do
506 : end do
507 :
508 5420008 : 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 753756 : 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 753756 : symrel1(:,:) = symrel(:,:)
824 :
825 : ! Compute determinant of the matrix
826 753756 : call mati3det(symrel1, det)
827 :
828 : ! Produce a rotation from an improper symmetry
829 5203668 : if (det==-1) symrel1(:,:) = -symrel1(:,:)
830 :
831 : ! Test the possibility of the unit matrix
832 753756 : identity(:,:)=0; identity(1,1)=1; identity(2,2)=1; identity(3,3)=1
833 :
834 9798828 : if (sum((symrel1(:,:) - identity(:,:))**2)/=0) then
835 :
836 : ! Transform symmetry matrix in the system defined by rprimd
837 702485 : call matr3inv(rprimd, rprimd_invt)
838 9132305 : coord=TRANSPOSE(rprimd_invt)
839 702485 : call matr3inv(coord,coordinvt)
840 36529220 : matr1(:,:) = MATMUL(symrel1,coord)
841 28099400 : matr2(:,:) = MATMUL(TRANSPOSE(coordinvt),matr1)
842 :
843 : ! Find the eigenvector with unit eigenvalue of the rotation matrix in cartesian coordinate, matr2
844 702485 : matr1(:,:)=matr2(:,:)
845 702485 : matr1(1,1)=matr1(1,1)-one
846 702485 : matr1(2,2)=matr1(2,2)-one
847 702485 : matr1(3,3)=matr1(3,3)-one
848 :
849 : ! Compute the axis of rotation and the cos and sin of rotation angle
850 2809940 : if(DOT_PRODUCT(matr1(:,1),matr1(:,1)) < tol8 )then
851 : ! The first direction is the axis
852 93604 : axis(1)=one ; axis(2)=zero ; axis(3)=zero
853 93604 : cos_phi=matr2(2,2)
854 93604 : sin_phi=matr2(3,2)
855 2435524 : else if(DOT_PRODUCT(matr1(:,2),matr1(:,2)) < tol8 )then
856 : ! The second direction is the axis
857 90153 : axis(1)=zero ; axis(2)=one ; axis(3)=zero
858 90153 : cos_phi=matr2(3,3)
859 90153 : 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 518728 : axis(1)=matr1(2,1)*matr1(3,2)-matr1(2,2)*matr1(3,1)
864 518728 : axis(2)=matr1(3,1)*matr1(1,2)-matr1(3,2)*matr1(1,1)
865 518728 : axis(3)=matr1(1,1)*matr1(2,2)-matr1(1,2)*matr1(2,1)
866 : ! Then, try to normalize it
867 2074912 : scprod=DOT_PRODUCT(axis(:), axis(:))
868 518728 : if(scprod<tol8)then
869 : ! The first and second vectors were linearly dependent
870 : ! Thus, use the first and third vectors
871 70465 : axis(1)=matr1(2,1)*matr1(3,3)-matr1(2,3)*matr1(3,1)
872 70465 : axis(2)=matr1(3,1)*matr1(1,3)-matr1(3,3)*matr1(1,1)
873 70465 : axis(3)=matr1(1,1)*matr1(2,3)-matr1(1,3)*matr1(2,1)
874 : ! Normalize the vector
875 281860 : scprod=DOT_PRODUCT(axis(:), axis(:))
876 70465 : if(scprod < tol8)then
877 0 : ABI_BUG('Cannot find the rotation axis.')
878 : end if
879 : end if
880 518728 : norminv=one/sqrt(scprod)
881 2074912 : 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 518728 : vecta(1)=one-axis(1)**2
887 518728 : vecta(2)=-axis(1)*axis(2)
888 518728 : vecta(3)=-axis(1)*axis(3)
889 2074912 : scprod=DOT_PRODUCT(vecta(:),vecta(:))
890 518728 : norminv=one/sqrt(scprod)
891 2074912 : vecta(:)=vecta(:)*norminv
892 : ! Rotate the vector A, to get vector B
893 2074912 : 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 2074912 : cos_phi=DOT_PRODUCT(vecta(:),vectb(:))
896 : ! Compute the cross product of the axis and vector A
897 518728 : vectc(1)=axis(2)*vecta(3)-axis(3)*vecta(2)
898 518728 : vectc(2)=axis(3)*vecta(1)-axis(1)*vecta(3)
899 518728 : 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 2074912 : 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 702485 : if(cos_phi> one-tol8 )cos_phi= one-tol8
907 : if(cos_phi<-(one-tol8))cos_phi=-(one-tol8)
908 702485 : phi=acos(cos_phi)
909 702485 : 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 702485 : phi=(nint(six*phi/pi))/six*pi
912 : ! Compute components of the spinor matrix
913 702485 : spinrot(1)=cos(half*phi)
914 702485 : spinrot(2)=axis(1)*sin(half*phi)
915 702485 : spinrot(3)=axis(2)*sin(half*phi)
916 702485 : spinrot(4)=axis(3)*sin(half*phi)
917 :
918 : else
919 :
920 : ! Here, the case of the unit matrix
921 : axis(:)=zero
922 51271 : phi=zero
923 51271 : spinrot(1)=one
924 51271 : spinrot(2)=zero
925 51271 : spinrot(3)=zero
926 51271 : 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 753756 : 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 4536 : 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 4536 : if (nsym > 1) then
1123 14498 : do iatom=1,natom
1124 261882 : do isym=1,nsym
1125 : ! jatom is the label of a symmetrically related atom
1126 247384 : 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 247384 : iatfix(2,jatom) /= iatfix(2,iatom) .or. &
1130 9962 : 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 4536 : 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 1915833 : 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 1915833 : 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 1915833 : 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 1915833 : 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 1915833 : call matr3inv(rprimd, gprimd)
1246 :
1247 : ! Write out rprimd, gprimd and ucvol
1248 1915833 : if (iout>=0) then
1249 7853 : write(msg,'(2a)')' Real(R)+Recip(G) ','space primitive vectors, cartesian coordinates (Bohr,Bohr^-1):'
1250 7853 : call wrtout(iout,msg)
1251 31412 : do nu=1,3
1252 : write(msg, '(1x,a,i1,a,3f11.7,2x,a,i1,a,3f11.7)' ) &
1253 94236 : 'R(',nu,')=',rprimd(:,nu)+tol10,&
1254 117795 : 'G(',nu,')=',gprimd(:,nu)+tol10
1255 31412 : call wrtout(iout,msg)
1256 : end do
1257 7853 : write(msg,'(a,1p,e15.7,a)') ' Unit cell volume ucvol=',ucvol+tol10,' bohr^3'
1258 7853 : call wrtout(iout,msg)
1259 7853 : call wrtout(std_out,msg)
1260 : end if
1261 :
1262 : ! Compute real space metric.
1263 76633320 : rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
1264 :
1265 : ! Compute reciprocal space metric.
1266 76633320 : gmet = MATMUL(TRANSPOSE(gprimd),gprimd)
1267 :
1268 : ! Write out the angles
1269 1915833 : if (iout>=0) then
1270 7853 : angle(1)=acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3)))/two_pi*360.0d0
1271 7853 : angle(2)=acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3)))/two_pi*360.0d0
1272 7853 : angle(3)=acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2)))/two_pi*360.0d0
1273 7853 : write(msg, '(a,3es16.8,a)' )' Angles (23,13,12)=',angle(1:3),' degrees'
1274 7853 : call wrtout(iout,msg)
1275 7853 : call wrtout(std_out,msg)
1276 : end if
1277 :
1278 1915833 : 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 90539 : 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 362156 : do ii=1,3
1524 1177007 : do jj=1,3
1525 1086468 : rprimd(ii,jj)=rprim(ii,jj)*acell(jj)
1526 : end do
1527 : end do
1528 :
1529 90539 : 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 24383 : 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 24383 : call matr3inv(rprimd,gprimd)
1572 9288385 : do iatom=1,natom
1573 37080391 : do mu=1,3
1574 37056008 : 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 24383 : 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 132608 : 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 28103778 : do iatom=1,natom
1617 112017288 : do mu=1,3
1618 111884680 : 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 132608 : 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 23761 : 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 23761 : favg(:)=zero
1668 89386 : do iatom=1,natom
1669 286261 : do mu=1,3
1670 : fcart(mu,iatom)= - (gprimd(mu,1)*gred(1,iatom)+&
1671 : & gprimd(mu,2)*gred(2,iatom)+&
1672 196875 : & gprimd(mu,3)*gred(3,iatom))
1673 262500 : 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 95044 : favg(:)=favg(:)/dble(natom)
1680 23761 : if(.not.Favgz_null) favg(3)=zero
1681 89386 : do iatom=1,natom
1682 286261 : fcart(:,iatom)=fcart(:,iatom)-favg(:)
1683 : end do
1684 :
1685 23761 : 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 7233 : 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 7233 : 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 7233 : 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 7233 : 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 17 : 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 34 : integer :: determinant(nsym),inversion(3,3),symrel2(3,3,nsym)
2697 34 : real(dp) :: dtnons(3),tnons2(3,nsym)
2698 : ! *********************************************************************
2699 :
2700 17 : ABI_WARNING('Removing inversion related symmetrie from initial set')
2701 :
2702 : ! Find the occurrence of the inversion symmetry.
2703 238 : call set2unit(inversion) ; inversion=-inversion
2704 :
2705 17 : is_inv=0; found=.FALSE.
2706 381 : do while (is_inv<nsym .and. .not.found)
2707 593 : is_inv=is_inv+1; found=ALL(symrel(:,:,is_inv)==inversion)
2708 : end do
2709 17 : if (found) then
2710 0 : write(msg,'(a,i3)')' The inversion is symmetry operation no. ',is_inv
2711 : else
2712 17 : write(msg,'(a)')' The inversion was not found in the symmetries list.'
2713 : end if
2714 17 : call wrtout(std_out,msg)
2715 :
2716 : ! Find the symmetries that are related through the inversion symmetry
2717 17 : call symdet(determinant,nsym,symrel)
2718 17 : nsym2=0
2719 364 : do is=1,nsym-1
2720 4506 : do is2=is+1,nsym
2721 :
2722 16568 : dtnons(:)=tnons(:,is2)-tnons(:,is)-tnons(:,is_inv)
2723 7077 : found=ALL(symrel(:,:,is)==-symrel(:,:,is2)).and.isinteger(dtnons,tol8)
2724 :
2725 347 : 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 17 : if (nsym2/=(nsym/2).or.nsym==1) then
2752 17 : call wrtout(std_out, ' Program uses the original set of symmetries ')
2753 17 : nsym_out=nsym
2754 51 : ABI_MALLOC(symrel_out,(3,3,nsym))
2755 51 : ABI_MALLOC(tnons_out,(3,nsym))
2756 4749 : symrel_out(:,:,:)=symrel(:,:,1:nsym)
2757 1473 : tnons_out(:,:)=tnons(:,1:nsym)
2758 17 : 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 17 : 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 789281 : 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 789281 : work=zero
2903 3157124 : do kk=1,3
2904 10260653 : do jj=1,3
2905 7103529 : symtmp=dble(symred(jj,kk))
2906 30781959 : do ii=1,3
2907 28414116 : 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 789281 : symcart=zero
2915 3157124 : do kk=1,3
2916 10260653 : do jj=1,3
2917 7103529 : symtmp=work(jj,kk)
2918 30781959 : do ii=1,3
2919 : ! symcart = aprim * work^T = aprim * symred * bprim^T
2920 28414116 : symcart(ii,jj)=symcart(ii,jj)+aprim(ii,kk)*symtmp
2921 : end do
2922 : end do
2923 : end do
2924 :
2925 789281 : 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 19939 : 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 19939 : call matr3inv(gprimd,rprimd)
3043 259207 : rprimdt=transpose(rprimd)
3044 :
3045 : !Compute stress tensor in reduced coordinates
3046 : ! strfrac = rprimd^T * stress * rprimd
3047 19939 : call strconv(stress,rprimdt,strfrac)
3048 :
3049 : !Switch to full storage mode
3050 19939 : tensor(1,1)=strfrac(1)
3051 19939 : tensor(2,2)=strfrac(2)
3052 19939 : tensor(3,3)=strfrac(3)
3053 19939 : tensor(3,2)=strfrac(4)
3054 19939 : tensor(3,1)=strfrac(5)
3055 19939 : tensor(2,1)=strfrac(6)
3056 19939 : tensor(2,3)=tensor(3,2)
3057 19939 : tensor(1,3)=tensor(3,1)
3058 19939 : tensor(1,2)=tensor(2,1)
3059 :
3060 : ! these loops are useless - trivial action:
3061 : ! tt = tensor / dble(nsym)
3062 : ! tensor = zero
3063 79756 : do nu=1,3
3064 259207 : do mu=1,3
3065 179451 : tt(mu,nu)=tensor(mu,nu)/dble(nsym)
3066 239268 : 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 348786 : do isym=1,nsym
3074 1335327 : do mu=1,3
3075 4275011 : do nu=1,3
3076 : summ=0._dp
3077 11838492 : do ii=1,3
3078 : tmp=tt(ii,1)*sym(nu,1,isym)+tt(ii,2)*sym(nu,2,isym)+&
3079 8878869 : & tt(ii,3)*sym(nu,3,isym)
3080 11838492 : summ=summ+sym(mu,ii,isym)*tmp
3081 : end do
3082 3946164 : 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 19939 : strfrac(1)=tensor(1,1)
3089 19939 : strfrac(2)=tensor(2,2)
3090 19939 : strfrac(3)=tensor(3,3)
3091 19939 : strfrac(4)=tensor(3,2)
3092 19939 : strfrac(5)=tensor(3,1)
3093 19939 : 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 19939 : call strconv(strfrac,gprimd,stress)
3100 :
3101 19939 : 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 5474 : 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 5474 : stress_mat(1,1) = stress6(1)
3123 5474 : stress_mat(2,2) = stress6(2)
3124 5474 : stress_mat(3,3) = stress6(3)
3125 5474 : stress_mat(2,3) = stress6(4)
3126 5474 : stress_mat(3,2) = stress6(4)
3127 5474 : stress_mat(1,3) = stress6(5)
3128 5474 : stress_mat(3,1) = stress6(5)
3129 5474 : stress_mat(1,2) = stress6(6)
3130 5474 : stress_mat(2,1) = stress6(6)
3131 :
3132 5474 : 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 931827 : 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 931827 : work1(1,1)=frac(1)
3176 931827 : work1(2,2)=frac(2)
3177 931827 : work1(3,3)=frac(3)
3178 931827 : work1(3,2)=frac(4) ; work1(2,3)=frac(4)
3179 931827 : work1(3,1)=frac(5) ; work1(1,3)=frac(5)
3180 931827 : 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 3727308 : do ii=1,3
3185 11181924 : work2(:,ii)=zero
3186 12113751 : do jj=1,3
3187 36341253 : 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 3727308 : do ii=1,3
3193 11181924 : work1(ii,:)=zero
3194 12113751 : do jj=1,3
3195 36341253 : work1(ii,:)=work1(ii,:)+gprimd(ii,jj)*work2(jj,:)
3196 : end do
3197 : end do
3198 :
3199 931827 : cart(1)=work1(1,1)
3200 931827 : cart(2)=work1(2,2)
3201 931827 : cart(3)=work1(3,3)
3202 931827 : cart(4)=work1(2,3)
3203 931827 : cart(5)=work1(1,3)
3204 931827 : cart(6)=work1(1,2)
3205 :
3206 931827 : 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 4115 : subroutine littlegroup_pert(gprimd,idir,indsym,iout,ipert,natom,nsym,nsym1, &
3259 4115 : & 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 4115 : ount = std_out; if (present(unit)) ount = unit
3284 :
3285 4115 : nsym1=0
3286 4115 : 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 3959 : 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 704 : nsym1=1
3321 2816 : do ii=1,3
3322 2112 : tnons1(ii,1)=0._dp
3323 2112 : symaf1(1)=1
3324 9152 : do jj=1,3
3325 6336 : symrl1(ii,jj,1)=0
3326 8448 : if(ii==jj)symrl1(ii,jj,1)=1
3327 : end do
3328 : end do
3329 :
3330 : else
3331 :
3332 116592 : 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 113337 : & symq(2,1,isym)==0 .and.&
3338 3255 : & symq(3,1,isym)==0 )then
3339 : ! Check that the symmetry operation preserves the atom
3340 29823 : if(ipert==indsym(4,isym,ipert))then
3341 : ! Check if the direction is preserved
3342 : tok=1
3343 111644 : do idir1=1,3
3344 83733 : if((idir1==idir.and.symrec(idir,idir1,isym)/=1) .or.&
3345 27911 : & (idir1/=idir.and.symrec(idir,idir1,isym)/=0))then
3346 83733 : tok=0
3347 : end if
3348 : end do
3349 27911 : if(tok==1)then
3350 : ! All the tests have been successful !
3351 5645 : nsym1=nsym1+1
3352 5645 : symaf1(nsym1)=symafm(isym)
3353 22580 : do ii=1,3
3354 16935 : tnons1(ii,nsym1)=tnons(ii,isym)
3355 73385 : do jj=1,3
3356 67740 : 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 4115 : 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 4115 : if (nsym1 /= 1) then
3372 1584 : if (iout /= ount .and. iout > 0) then
3373 702 : write(msg,'(a,i5,a)')' Found ',nsym1,' symmetries that leave the perturbation invariant.'
3374 702 : call wrtout(iout,msg)
3375 : end if
3376 1584 : write(msg,'(a,i5,a)')' littlegroup_pert: found ',nsym1,' symmetries that leave the perturbation invariant: '
3377 1584 : call wrtout(ount,msg)
3378 : else
3379 2531 : if (iout /= ount .and. iout > 0) then
3380 1565 : write(msg,'(a,a)')' The set of symmetries contains',' only one element for this perturbation.'
3381 1565 : call wrtout(iout,msg)
3382 : end if
3383 2531 : write(msg,'(a)')' littlegroup_pert: only one element in the set of symmetries for this perturbation:'
3384 2531 : call wrtout(ount,msg)
3385 : end if
3386 :
3387 4115 : if (ount > 0) then
3388 6782 : do isym=1,nsym1
3389 4385 : write(msg, '(9i4)' )((symrl1(ii,jj,isym),ii=1,3),jj=1,3)
3390 6782 : call wrtout(ount,msg)
3391 : end do
3392 : end if
3393 :
3394 4115 : 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 1518 : 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 53722 : pertsy(:,:)=0
3454 :
3455 14569 : do ipert1=1,mpert
3456 53722 : do idir1=1,3
3457 52204 : 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 255019 : do isym=1,nsym
3462 672538 : do itirev=1,2
3463 457742 : isign=3-2*itirev
3464 :
3465 664974 : if(symq(4,itirev,isym)/=0)then
3466 :
3467 181554 : found=1
3468 :
3469 : ! Here select the symmetric of ipert1
3470 181554 : if(ipert1<=natom)then
3471 137465 : ipesy1=indsym(4,isym,ipert1)
3472 549860 : do ii=1,3
3473 1787045 : do jj=1,3
3474 1649580 : sym1(ii,jj)=symrec(ii,jj,isym)
3475 : end do
3476 : end do
3477 44089 : 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 384588 : do idisy1=1,3
3494 384588 : if(sym1(idir1,idisy1)/=0)then
3495 185707 : 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 181554 : 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 40223 : pertsy(idir1,ipert1)=-1
3516 40223 : 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 7564 : 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 1518 : 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 8928 : 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 8928 : pi=4.0d0*datan(1.0d0)
4125 : ! Compute spherical coordinates
4126 8928 : length = sqrt(vector(1)**2+vector(2)**2+vector(3)**2)
4127 :
4128 8928 : if (length > tol6) then
4129 8374 : theta = acos(vector(3) / length) * 180.d0 / pi
4130 8374 : if (abs(vector(1)) > tol6 .or. abs(vector(2)) > tol6) then
4131 7481 : phi = atan2(vector(2), vector(1)) * 180.d0 / pi
4132 7481 : if (phi<-179) phi=180
4133 : else
4134 893 : phi = 0.d0
4135 : end if
4136 : else
4137 554 : theta = 0.d0
4138 554 : phi = 0.d0
4139 : end if
4140 :
4141 8928 : 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 : !!***
|
