Line data Source code
1 : !!****m* ABINIT/m_crystal
2 : !! NAME
3 : !! m_crystal
4 : !!
5 : !! FUNCTION
6 : !! Module containing the definition of the crystal_t data type and methods used to handle it.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG, YP, MJV, GA)
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_crystal
23 :
24 : use, intrinsic :: iso_c_binding
25 : use defs_basis
26 : use m_errors
27 : use m_abicore
28 : use m_atomdata
29 : use m_xmpi
30 : use m_nctk
31 : use netcdf
32 :
33 : use m_io_tools, only : file_exists
34 : use m_numeric_tools, only : set2unit
35 : use m_hide_lapack, only : matrginv
36 : use m_fstrings, only : int2char10, sjoin, yesno, itoa, strcat
37 : use m_matrix, only : mati3inv
38 : use m_symtk, only : sg_multable, symatm, print_symmetries
39 : use m_spgdata, only : spgdata
40 : use m_geometry, only : metric, xred2xcart, xcart2xred, remove_inversion, getspinrot, symredcart, normv
41 : use m_io_tools, only : open_file
42 : use m_pptools, only : write_xsf
43 :
44 : implicit none
45 :
46 : private
47 : !!***
48 :
49 : !----------------------------------------------------------------------
50 :
51 : !!****t* m_crystal/crystal_t
52 : !! NAME
53 : !! crystal_t
54 : !!
55 : !! FUNCTION
56 : !! Structure defining the unit cell (geometry, atomic positions and symmetry operations in real and reciprocal space)
57 : !!
58 : !! SOURCE
59 :
60 : type,public :: crystal_t
61 :
62 : !scalars
63 : !integer :: point_group ! Point group
64 : !integer :: bravais,crystsys ! Bravais lattice, Crystal system
65 : !integer :: nptsym ! No of point symmetries of the Bravais lattice
66 : !integer :: bravais(11) ! bravais(1)=iholohedry, bravais(2)=center
67 : ! bravais(3:11)=coordinates of rprim in the axes of the conventional
68 : ! bravais lattice (*2 if center/=0)
69 : !integer,pointer ptsymrel(:,:,:)
70 : !ptsymrel(3,3,nptsym)
71 : ! nptsym point-symmetry operations of the Bravais lattice in real space in terms of primitive translations.
72 :
73 : integer :: natom = -1
74 : ! Number of atoms
75 :
76 : integer :: nsym = -1
77 : ! Number of symmetry operations
78 :
79 : integer :: ntypat = -1
80 : ! Number of type of atoms
81 :
82 : integer :: nirredat = -1
83 : ! Number of irreducibel atoms
84 :
85 : integer :: npsp = -1
86 : ! No. of pseudopotentials
87 :
88 : integer :: space_group = -1
89 : ! Space group
90 :
91 : integer :: timrev = -1
92 : ! TODO BE CAREFUL here, as the convention used in abinit is different.
93 : ! 1 => do not use time-reversal symmetry.
94 : ! 2 => take advantage of time-reversal symmetry.
95 :
96 : real(dp) :: ucvol = zero
97 : ! Real space unit cell volume.
98 :
99 : logical :: use_antiferro
100 : ! .TRUE. if AFM symmetries are present and used.
101 :
102 : !arrays
103 : real(dp) :: angdeg(3) = zero
104 : ! Angles among rprim (degree).
105 :
106 : real(dp) :: gmet(3,3) = zero
107 : ! Reciprocal space metric ($\textrm{bohr}^{-2}$).
108 :
109 : real(dp) :: gprimd(3,3) = zero
110 : ! Dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
111 :
112 : real(dp) :: rmet(3,3) = zero
113 : ! Metric in real space.
114 :
115 : real(dp) :: rprimd(3,3) = zero
116 : ! Direct lattice vectors, Bohr units.
117 :
118 : integer,allocatable :: indsym(:,:,:)
119 : ! indsym(4,nsym,natom)
120 : ! indirect indexing array for atoms, see symatm.F90.
121 : ! $ R^{-1} (xred(:,iat)-\tau) = xred(:,iat_sym) + R_0 $
122 : ! * indsym(4, isym,iat) gives iat_sym in the original unit cell.
123 : ! * indsym(1:3,isym,iat) gives the lattice vector $R_0$.
124 :
125 : integer,allocatable :: symafm(:)
126 : ! symafm(nsym)
127 : ! (Anti)Ferromagnetic symmetries. +1/-1
128 :
129 : integer,allocatable :: symrec(:,:,:)
130 : ! symrec(3,3,nsym)
131 : ! Symmetry operation in reciprocal space (reduced coordinates)
132 :
133 : integer,allocatable :: symrel(:,:,:)
134 : ! symrel(3,3,nsym)
135 : ! Symmetry operations in direct space (reduced coordinates).
136 :
137 : real(dp),allocatable :: symrel_cart(:,:,:)
138 : ! symrel_cart(3,3,nsym)
139 : ! Symmetry operations in cartesian coordinates (same order as symrel)
140 :
141 : integer,allocatable :: atindx(:)
142 : integer,allocatable :: atindx1(:)
143 : ! atindx(natom), atindx1(natom)
144 : ! Index tables for atoms useful to treat atoms type after type.
145 :
146 : integer,allocatable :: typat(:)
147 : integer,allocatable :: nattyp(:)
148 : ! typat(natom), nattyp(ntypat)
149 : ! Type of each natom and number of atoms of each type.
150 :
151 : integer,allocatable :: irredatindx(:)
152 : ! Index of irreducible atoms
153 :
154 : real(dp),allocatable :: tnons(:,:)
155 : ! tnons(3,nsym)
156 : ! Fractional translations (reduced coordinates)
157 :
158 : real(dp),allocatable :: xcart(:,:)
159 : ! xcart(3,natom)
160 : ! Cartesian coordinates.
161 :
162 : real(dp),allocatable :: xred(:,:)
163 : ! xred(3,natom)
164 : ! Reduced coordinates.
165 :
166 : real(dp),allocatable :: spinrot(:,:)
167 : ! spinrot(4,nsym)
168 : ! spinor rotation matrices.
169 :
170 : real(dp),allocatable :: amu(:)
171 : ! amu(ntypat)
172 : ! mass of the atoms (atomic mass unit)
173 :
174 : real(dp),allocatable :: zion(:)
175 : ! zion(ntypat)
176 : ! Charge of the pseudo-ion
177 : ! (No of valence electrons needed to screen exactly the pseudopotential).
178 :
179 : !real(dp),allocatable :: znucltypat(:)
180 : ! znucltypat(ntypat)
181 : ! The atomic number of each type of atom (might be alchemy wrt psps)
182 :
183 : real(dp),allocatable :: znucl(:)
184 : ! znucl(npsp)
185 : ! Nuclear charge for each type of pseudopotential
186 :
187 : character(len=132),allocatable :: title(:)
188 : ! title(ntypat)
189 : ! The content of first line read from the psp file
190 :
191 : contains
192 :
193 : procedure :: init => crystal_init
194 : ! Main Creation method.
195 :
196 : procedure :: ncwrite => crystal_ncwrite
197 : ! Write the object in netcdf format
198 :
199 : procedure :: ncwrite_path => crystal_ncwrite_path
200 : ! Dump the object to netcdf file.
201 :
202 : procedure :: ncread => crystal_ncread
203 : ! Read the object from a netcdf file.
204 :
205 : procedure :: isymmorphic
206 : ! True if space group is symmorphic.
207 :
208 : procedure :: idx_spatial_inversion
209 : ! Return the index of the spatial inversion, 0 if not present.
210 :
211 : procedure :: isalchemical
212 : ! True if we are using alchemical pseudopotentials.
213 :
214 : procedure :: malloc => crystal_malloc
215 : ! Allocate memory.
216 :
217 : procedure :: free => crystal_free
218 : ! Free memory.
219 :
220 : procedure :: copy => crystal_copy
221 : ! Copy object.
222 :
223 : procedure :: bcast => crystal_bcast
224 : ! Master broadcasts data and others allocate their arrays.
225 :
226 : procedure :: new_without_symmetries => crystal_without_symmetries
227 : ! Return new object without symmetries (actually nsym = 1 and identity operation)
228 :
229 : procedure :: new_trinv_only => crystal_trinv_only
230 : ! Return new object without only identity, inversion & TR symmetries (if present)
231 :
232 : procedure :: get_point_group => crystal_point_group
233 : ! Return the symmetries of the point group of the crystal.
234 :
235 : procedure :: index_atoms => crystal_index_atoms
236 : ! Generate index table of atoms.
237 :
238 : procedure :: compute_sym => crystal_compute_sym
239 : ! Compute all symetries and construct tables.
240 :
241 : procedure :: compute_geometry => crystal_compute_geometry
242 : ! Compute the different metrics and the angle between primitive vectors.
243 :
244 : procedure :: symbol_type
245 : ! Return the atomic symbol from the itypat index.
246 :
247 : procedure :: symbol_iatom
248 : ! Return the atomic symbol from the iatom index.
249 :
250 : procedure :: adata_type
251 : ! Return atomic data from the itypat index.
252 :
253 : procedure :: compare => crystal_compare
254 : ! Compare two crystalline structures, write warning messages if they differ, return exit status
255 :
256 : procedure :: print => crystal_print
257 : ! Print dimensions and basic info stored in the object
258 :
259 : procedure :: print_abivars => crystal_print_abivars
260 : ! Print unit cell info in Abinit/abivars format
261 :
262 : procedure :: symmetrize_cart_vec3 => crystal_symmetrize_cart_vec3
263 : ! Symmetrize a 3d cartesian vector
264 :
265 : procedure :: symmetrize_cart_tens33 => crystal_symmetrize_cart_tens33
266 : ! Symmetrize a cartesian 3x3 tensor
267 :
268 : procedure :: get_redcart_qdirs => get_redcart_qdirs
269 : ! Return predefined list of 6 q-versors in reciprocal space reduced coordinates.
270 :
271 : procedure :: write_xsf_data => crystal_write_xsf_data
272 : ! Write array in real space in XSF format
273 :
274 : end type crystal_t
275 :
276 : public :: symbols_crystal ! Return an array with the atomic symbol: ["Sr","Ru","O1","O2","O3"]
277 : public :: prt_cif ! Print CIF file.
278 : public :: prtposcar ! output VASP style POSCAR and FORCES files.
279 : !!***
280 :
281 : CONTAINS !====================================================================================================
282 : !!***
283 :
284 : !!****f* m_crystal/crystal_init
285 : !! NAME
286 : !! crystal_init
287 : !!
288 : !! FUNCTION
289 : !! Initialize a crystal_t data type.
290 : !! Ideally the routine should work in two different modes:
291 : !! Either the symmetries are directly supplied or the space group
292 : !! is determined starting from the definition of the unit cell.
293 : !! Only the first method is implemented, the second one should be
294 : !! a wrapper for the symmetry finder library. To implement the
295 : !! second case I have to add additional entries in the object
296 : !! and I have also to pass an object describing the (optional) geometry builder.
297 : !!
298 : !! INPUTS
299 : !! natom=number of atom
300 : !! ntypat=number of type of atoms
301 : !! nsym=number of symmetry operations
302 : !! rprimd(3,3)=dimensional lattive vector (real space)
303 : !! typat(natom)=type of each atom
304 : !! xred(3,natom)=reduced coordinates of each atom
305 : !! symrel(3,3,nsym) [optional]=symmetry operations in real space
306 : !! space_group=Space group (0 if not available)
307 : !! tnons(3,nsym) [optional]=fractional Translations
308 : !! symafm(nsym) [optional]= ferromagnetic symmetries
309 : !! remove_inv [optional]= if .TRUE. the inversion is removed from the set of symmetries
310 : !! timrev ==2 => take advantage of time-reversal symmetry
311 : !! ==1 ==> do not use time-reversal symmetry
312 : !!
313 : !! OUTPUT
314 : !! Cryst<crystal_t>= the object completely initialized.
315 : !!
316 : !! TODO
317 : !! Add additional entries in the class:
318 : !! 1) Info on space and point group (generators?).
319 : !! 2) alchemy
320 : !! 3) masses and nuclear (pseudo&AE) charge
321 : !! 4) forces stresses, velocities.
322 : !! 5) constraints for the relaxation
323 : !! 6) Likely I will need also info on the electric field and berryopt
324 : !!
325 : !! SOURCE
326 :
327 31938 : subroutine crystal_init(cryst,amu,space_group,natom,npsp,ntypat,nsym,rprimd,typat,xred,&
328 31938 : zion,znucl,timrev,use_antiferro,remove_inv,title,&
329 95814 : symrel,tnons,symafm) ! Optional
330 :
331 : !Arguments ------------------------------------
332 : !scalars
333 : class(crystal_t),intent(inout) :: Cryst
334 : integer,intent(in) :: natom,ntypat,nsym,timrev,space_group,npsp
335 : logical,intent(in) :: remove_inv,use_antiferro
336 : !arrays
337 : integer,intent(in) :: typat(natom)
338 : integer,optional,intent(in) :: symrel(3,3,nsym),symafm(nsym)
339 : real(dp),intent(in) :: amu(ntypat),xred(3,natom),rprimd(3,3),zion(ntypat),znucl(npsp)
340 : real(dp),optional,intent(in) :: tnons(3,nsym)
341 : character(len=*),intent(in) :: title(ntypat)
342 :
343 : !Local variables-------------------------------
344 : !scalars
345 : integer :: pinv,nsym_noI
346 : !character(len=500) :: msg
347 : !arrays
348 31938 : integer,pointer :: symrel_noI(:,:,:)
349 31938 : real(dp),pointer :: tnons_noI(:,:)
350 : ! *************************************************************************
351 :
352 : !@crystal_t
353 31938 : Cryst%natom = natom
354 31938 : Cryst%ntypat = ntypat
355 31938 : Cryst%npsp = npsp
356 31938 : Cryst%space_group = space_group
357 31938 : Cryst%nsym = nsym
358 31938 : Cryst%timrev = timrev
359 31938 : Cryst%use_antiferro = use_antiferro
360 415194 : Cryst%rprimd = rprimd
361 :
362 31938 : call Cryst%free()
363 31938 : call Cryst%malloc()
364 :
365 107634 : Cryst%amu = amu
366 145847 : Cryst%typat = typat
367 391760 : Cryst%xred = xred
368 107634 : Cryst%zion = zion
369 107724 : Cryst%znucl = znucl
370 107634 : Cryst%title = title
371 :
372 31938 : call Cryst%compute_geometry()
373 :
374 31938 : call Cryst%index_atoms()
375 :
376 : ! TODO: Make this more elegant
377 31938 : if (PRESENT(symrel).and.PRESENT(tnons).and.PRESENT(symafm)) then
378 31938 : if (.not.remove_inv) then
379 : ! Just a copy
380 9726473 : Cryst%symrel=symrel
381 3036980 : Cryst%tnons=tnons
382 807149 : Cryst%symafm=symafm
383 : else
384 : ! Remove inversion, just to be compatible with old GW implementation
385 : ! TODO should be removed!
386 2 : call remove_inversion(nsym,symrel,tnons,nsym_noI,symrel_noI,tnons_noI,pinv)
387 2 : Cryst%nsym=nsym_noI
388 2 : ABI_SFREE(Cryst%symrel)
389 2 : ABI_SFREE(Cryst%symrec)
390 2 : ABI_SFREE(Cryst%tnons)
391 2 : ABI_SFREE(Cryst%symafm)
392 6 : ABI_MALLOC(Cryst%symrel,(3,3,nsym_noI))
393 4 : ABI_MALLOC(Cryst%symrec,(3,3,nsym_noI))
394 6 : ABI_MALLOC(Cryst%tnons,(3,nsym_noI))
395 6 : ABI_MALLOC(Cryst%symafm,(nsym_noI))
396 628 : Cryst%symrel=symrel_noI
397 196 : Cryst%tnons=tnons_noI
398 50 : if (ANY(symafm==-1)) then
399 0 : ABI_BUG('Solve the problem with inversion before adding ferromagnetic symmetries')
400 : end if
401 50 : Cryst%symafm=1
402 2 : ABI_FREE(symrel_noI)
403 4 : ABI_FREE(tnons_noI)
404 : end if
405 :
406 : else
407 : ! Find symmetries symrec,symrel,tnons,symafm
408 : ! TODO This should be a wrapper around the abinit library whose usage is not so straightforward
409 0 : ABI_BUG('NotImplememented: symrel, symrec and tnons should be specied')
410 : end if
411 :
412 : ! Compute all symetries and construct tables.
413 31938 : call Cryst%compute_sym()
414 :
415 127752 : end subroutine crystal_init
416 : !!***
417 :
418 : !----------------------------------------------------------------------
419 :
420 : !!****f* m_crystal/crystal_without_symmetries
421 : !! NAME
422 : !! crystal_without_symmetries
423 : !!
424 : !! FUNCTION
425 : !! Return new crystal_t object without symmetries (actually nsym = 1 and identity operation)
426 : !!
427 : !! INPUTS
428 : !!
429 : !! OUTPUT
430 : !!
431 : !! SOURCE
432 :
433 0 : type(crystal_t) function crystal_without_symmetries(self) result(new)
434 :
435 : !Arguments ------------------------------------
436 : class(crystal_t), intent(in) :: self
437 :
438 : !Local variables-------------------------------
439 : integer,parameter :: timrev1 = 1, new_symafm(1) = 1
440 : real(dp),parameter :: new_tnons(3,1) = zero
441 : ! *************************************************************************
442 :
443 : call new%init(self%amu, 1, self%natom, self%npsp, self%ntypat, 1, self%rprimd, self%typat, &
444 : self%xred, self%zion, self%znucl, timrev1, .False., .False., self%title, &
445 0 : symrel=identity_3d, tnons=new_tnons, symafm=new_symafm)
446 :
447 0 : end function crystal_without_symmetries
448 : !!***
449 :
450 : !----------------------------------------------------------------------
451 :
452 : !!****f* m_crystal/crystal_trinv_only
453 : !! NAME
454 : !! crystal_trinv_only
455 : !!
456 : !! FUNCTION
457 : ! ! Return new crystal_t without only identity, inversion and TR symmetries (if present)
458 : !!
459 : !! SOURCE
460 :
461 0 : type(crystal_t) function crystal_trinv_only(self) result(new)
462 :
463 : !Arguments ------------------------------------
464 : class(crystal_t), intent(in) :: self
465 :
466 : !Local variables-------------------------------
467 : integer,parameter :: timrev2 = 2
468 : real(dp),parameter :: new_tnons(3,2) = zero
469 : integer :: inv_idx, new_symafm(2), new_symrel(3,3,2)
470 : ! *************************************************************************
471 :
472 0 : inv_idx = self%idx_spatial_inversion()
473 0 : if (inv_idx == 0) then
474 : ! no spatial inversion
475 0 : new = self%new_without_symmetries(); new%timrev = timrev2
476 : else
477 : ! spatial inversion is present
478 0 : new_symrel(:,:,1) = identity_3d; new_symrel(:,:,2) = self%symrel(:,:,inv_idx)
479 0 : new_symafm(1) = 1; new_symafm(2) = self%symafm(inv_idx)
480 :
481 : call new%init(self%amu, 2, self%natom, self%npsp, self%ntypat, 2, self%rprimd, self%typat, &
482 : self%xred, self%zion, self%znucl, timrev2, .False., .False., self%title, &
483 0 : symrel=new_symrel, tnons=new_tnons, symafm=new_symafm)
484 : endif
485 :
486 0 : end function crystal_trinv_only
487 : !!***
488 :
489 : !----------------------------------------------------------------------
490 :
491 : !!****f* m_crystal/crystal_index_atoms
492 : !! NAME
493 : !! crystal_index_atoms
494 : !!
495 : !! FUNCTION
496 : !! Generate index table of atoms, in order for them to be used type after type.
497 : !!
498 : !! SOURCE
499 :
500 32339 : subroutine crystal_index_atoms(Cryst)
501 :
502 : !Arguments ------------------------------------
503 : class(crystal_t),intent(inout) :: Cryst
504 :
505 : !Local variables-------------------------------
506 : integer :: iat,indx,itypat
507 : ! *********************************************************************
508 :
509 32339 : indx=1
510 76644 : do itypat=1,Cryst%ntypat
511 44305 : Cryst%nattyp(itypat)=0
512 197888 : do iat=1,Cryst%natom
513 165549 : if (Cryst%typat(iat)==itypat) then
514 82757 : Cryst%atindx (iat )=indx
515 82757 : Cryst%atindx1(indx)=iat
516 82757 : indx=indx+1
517 82757 : Cryst%nattyp(itypat)=Cryst%nattyp(itypat)+1
518 : end if
519 : end do
520 : end do
521 :
522 32339 : end subroutine crystal_index_atoms
523 : !!***
524 :
525 : !----------------------------------------------------------------------
526 :
527 : !!****f* m_crystal/crystal_compute_sym
528 : !! NAME
529 : !! crystal_compute_sym
530 : !!
531 : !! FUNCTION
532 : !! Get symmetries in cartesian coordinates, construct rotation tables
533 : !! for atoms with and without spinor, and construct list of reducible atoms.
534 : !!
535 : !! SOURCE
536 :
537 32339 : subroutine crystal_compute_sym(Cryst)
538 :
539 : !Arguments ------------------------------------
540 : class(crystal_t),intent(inout) :: Cryst
541 :
542 : !Local variables-------------------------------
543 : !scalars
544 : integer :: iat,indx,isym
545 : real(dp) :: tolsym8
546 32339 : logical, allocatable :: irredat_tmp(:)
547 : !arrays
548 : integer :: symrec(3,3)
549 : ! *********************************************************************
550 :
551 : ! Get symmetries in reciprocal space
552 785523 : do isym=1,Cryst%nsym
553 753184 : call mati3inv(Cryst%symrel(:,:,isym),symrec)
554 9823731 : Cryst%symrec(:,:,isym)=symrec
555 : end do
556 :
557 : ! Get symmetries in cartesian coordinates
558 785523 : do isym =1,Cryst%nsym
559 753184 : call symredcart(Cryst%rprimd, Cryst%gprimd, Cryst%symrel_cart(:,:,isym), Cryst%symrel(:,:,isym))
560 : ! purify operations in cartesian coordinates.
561 9823731 : where (abs(Cryst%symrel_cart(:,:,isym)) < tol14)
562 : Cryst%symrel_cart(:,:,isym) = zero
563 : end where
564 : end do
565 :
566 : ! === Obtain a list of rotated atoms ===
567 : ! $ R^{-1} (xred(:,iat)-\tau) = xred(:,iat_sym) + R_0 $
568 : ! * indsym(4, isym,iat) gives iat_sym in the original unit cell.
569 : ! * indsym(1:3,isym,iat) gives the lattice vector $R_0$.
570 : !
571 32339 : tolsym8=tol8
572 32339 : call symatm(Cryst%indsym, Cryst%natom, Cryst%nsym, Cryst%symrec, Cryst%tnons, tolsym8, Cryst%typat, Cryst%xred)
573 :
574 : ! Rotations in spinor space
575 785523 : do isym=1,Cryst%nsym
576 785523 : call getspinrot(Cryst%rprimd, Cryst%spinrot(:,isym), Cryst%symrel(:,:,isym))
577 : end do
578 :
579 : ! Find list of irreducible atoms by using the indsym
580 97017 : ABI_MALLOC(irredat_tmp, (Cryst%natom))
581 115096 : irredat_tmp = .TRUE.
582 :
583 32339 : Cryst%nirredat = 0
584 115096 : do iat = 1,Cryst%natom
585 115096 : if(irredat_tmp(iat))then
586 65172 : Cryst%nirredat = Cryst%nirredat + 1
587 1037801 : do isym = 1,Cryst%nsym
588 1037801 : if (Cryst%indsym(4,isym,iat) /= iat)then
589 248167 : irredat_tmp(Cryst%indsym(4,isym,iat)) = .FALSE.
590 : endif
591 : enddo
592 : endif
593 : enddo
594 :
595 : ! Write indexes of irreducible atoms
596 97017 : ABI_MALLOC(Cryst%irredatindx, (Cryst%nirredat))
597 32339 : indx = 0
598 115096 : do iat = 1,Cryst%natom
599 115096 : if (irredat_tmp(iat)) then
600 65172 : indx = indx + 1
601 65172 : cryst%irredatindx(indx) = iat
602 : endif
603 : enddo
604 :
605 32339 : ABI_SFREE(irredat_tmp)
606 :
607 32339 : end subroutine crystal_compute_sym
608 : !!***
609 :
610 : !----------------------------------------------------------------------
611 :
612 : !!****f* m_crystal/crystal_compute_geometry
613 : !! NAME
614 : !! crystal_compute_geometry
615 : !!
616 : !! FUNCTION
617 : !! Compute the different metrics and the angle between primitive vectors.
618 : !! Also compute cartesian coordinates of atoms.
619 : !!
620 : !! SOURCE
621 :
622 32339 : subroutine crystal_compute_geometry(Cryst)
623 :
624 : !Arguments ------------------------------------
625 : class(crystal_t),intent(inout) :: Cryst
626 : ! *********************************************************************
627 :
628 32339 : call metric(Cryst%gmet,Cryst%gprimd,-1,Cryst%rmet,Cryst%rprimd,Cryst%ucvol)
629 :
630 32339 : Cryst%angdeg(1)=ACOS(Cryst%rmet(2,3)/SQRT(Cryst%rmet(2,2)*Cryst%rmet(3,3)))/two_pi*360.0d0
631 32339 : Cryst%angdeg(2)=ACOS(Cryst%rmet(1,3)/SQRT(Cryst%rmet(1,1)*Cryst%rmet(3,3)))/two_pi*360.0d0
632 32339 : Cryst%angdeg(3)=ACOS(Cryst%rmet(1,2)/SQRT(Cryst%rmet(1,1)*Cryst%rmet(2,2)))/two_pi*360.0d0
633 :
634 32339 : call xred2xcart(Cryst%natom,Cryst%rprimd,Cryst%xcart,Cryst%xred)
635 :
636 32339 : end subroutine crystal_compute_geometry
637 : !!***
638 :
639 : !----------------------------------------------------------------------
640 :
641 : !!****f* m_crystal/crystal_malloc
642 : !! NAME
643 : !! crystal_malloc
644 : !!
645 : !! FUNCTION
646 : !! Allocate the dynamic arrays in a crystal_t data type.
647 : !!
648 : !! SOURCE
649 :
650 32884 : subroutine crystal_malloc(Cryst)
651 :
652 : !Arguments ------------------------------------
653 : class(crystal_t),intent(inout) :: Cryst
654 :
655 : !Local variables-------------------------------
656 : integer :: ii
657 : ! *********************************************************************
658 :
659 : !integer
660 98652 : ABI_MALLOC(Cryst%typat,(Cryst%natom))
661 98652 : ABI_MALLOC(Cryst%xred,(3,Cryst%natom))
662 65768 : ABI_MALLOC(Cryst%xcart,(3,Cryst%natom))
663 98652 : ABI_MALLOC(Cryst%zion,(Cryst%ntypat))
664 98652 : ABI_MALLOC(Cryst%znucl,(Cryst%npsp))
665 65768 : ABI_MALLOC(Cryst%amu, (Cryst%ntypat))
666 :
667 98652 : ABI_MALLOC(Cryst%symrel,(3,3,Cryst%nsym))
668 65768 : ABI_MALLOC(Cryst%symrec,(3,3,Cryst%nsym))
669 98652 : ABI_MALLOC(Cryst%tnons,(3,Cryst%nsym))
670 98652 : ABI_MALLOC(Cryst%symafm,(Cryst%nsym))
671 98652 : ABI_MALLOC(Cryst%symrel_cart, (3, 3, Cryst%nsym))
672 131536 : ABI_MALLOC(Cryst%indsym,(4, Cryst%nsym, Cryst%natom))
673 :
674 65768 : ABI_MALLOC(Cryst%atindx,(Cryst%natom))
675 65768 : ABI_MALLOC(Cryst%atindx1,(Cryst%natom))
676 98652 : ABI_MALLOC(Cryst%nattyp,(Cryst%ntypat))
677 98652 : ABI_MALLOC(Cryst%spinrot, (4, Cryst%nsym))
678 :
679 98652 : ABI_MALLOC(Cryst%title,(Cryst%ntypat))
680 :
681 : ! nirredat must first be computed from indsym
682 : !ABI_MALLOC(Cryst%irredatindx,(Cryst%nirredat))
683 :
684 117047 : Cryst%typat = zero
685 369536 : Cryst%xred = zero
686 369536 : Cryst%xcart = zero
687 78115 : Cryst%zion = zero
688 78205 : Cryst%znucl = zero
689 78115 : Cryst%amu = zero
690 10063060 : Cryst%symrel = zero
691 10063060 : Cryst%symrec = zero
692 3119092 : Cryst%tnons = zero
693 804436 : Cryst%symafm = zero
694 10063060 : Cryst%symrel_cart = zero
695 8604972 : Cryst%indsym = zero
696 117047 : Cryst%atindx = zero
697 117047 : Cryst%atindx1 = zero
698 78115 : Cryst%nattyp = zero
699 3890644 : Cryst%spinrot = zero
700 :
701 78115 : do ii=1,Cryst%ntypat
702 78115 : Cryst%title(ii) = ''
703 : end do
704 :
705 32884 : end subroutine crystal_malloc
706 : !!***
707 :
708 : !----------------------------------------------------------------------
709 :
710 : !!****f* m_crystal/crystal_free
711 : !! NAME
712 : !! crystal_free
713 : !!
714 : !! FUNCTION
715 : !! Free dynamic memory in a crystal_t data type.
716 : !!
717 : !! SOURCE
718 :
719 66144 : subroutine crystal_free(Cryst)
720 :
721 : !Arguments ------------------------------------
722 : class(crystal_t),intent(inout) :: Cryst
723 : ! *********************************************************************
724 :
725 : !integer
726 66144 : ABI_SFREE(Cryst%indsym)
727 66144 : ABI_SFREE(Cryst%symafm)
728 66144 : ABI_SFREE(Cryst%symrec)
729 66144 : ABI_SFREE(Cryst%symrel)
730 66144 : ABI_SFREE(Cryst%symrel_cart)
731 66144 : ABI_SFREE(Cryst%atindx)
732 66144 : ABI_SFREE(Cryst%atindx1)
733 66144 : ABI_SFREE(Cryst%typat)
734 66144 : ABI_SFREE(Cryst%nattyp)
735 66144 : ABI_SFREE(Cryst%irredatindx)
736 :
737 : !real
738 66144 : ABI_SFREE(Cryst%tnons)
739 66144 : ABI_SFREE(Cryst%xcart)
740 66144 : ABI_SFREE(Cryst%xred)
741 66144 : ABI_SFREE(Cryst%zion)
742 66144 : ABI_SFREE(Cryst%znucl)
743 66144 : ABI_SFREE(Cryst%amu)
744 66144 : ABI_SFREE(Cryst%spinrot)
745 :
746 : !character
747 66144 : ABI_SFREE(Cryst%title)
748 :
749 66144 : end subroutine crystal_free
750 : !!***
751 :
752 : !----------------------------------------------------------------------
753 :
754 : !!****f* m_crystal/crystal_copy
755 : !! NAME
756 : !! crystal_copy
757 : !!
758 : !! FUNCTION
759 : !! Copy the object.
760 : !!
761 : !! OUTPUTS
762 : !! new = A new crystal instance
763 : !!
764 : !! SOURCE
765 :
766 529 : subroutine crystal_copy(Cryst, new)
767 :
768 : !Arguments ------------------------------------
769 : class(crystal_t),intent(in) :: Cryst
770 : class(crystal_t),intent(out) :: new
771 : ! *********************************************************************
772 :
773 : ! Copy dimensions, scalar variables, and static arrays
774 529 : new%natom = Cryst%natom
775 529 : new%nsym = Cryst%nsym
776 529 : new%ntypat = Cryst%ntypat
777 529 : new%nirredat = Cryst%nirredat
778 529 : new%npsp = Cryst%npsp
779 529 : new%space_group = Cryst%space_group
780 529 : new%timrev = Cryst%timrev
781 529 : new%ucvol = Cryst%ucvol
782 529 : new%use_antiferro = Cryst%use_antiferro
783 2116 : new%angdeg = Cryst%angdeg
784 6877 : new%gmet = Cryst%gmet
785 6877 : new%gprimd = Cryst%gprimd
786 6877 : new%rmet = Cryst%rmet
787 6877 : new%rprimd = Cryst%rprimd
788 :
789 : ! Allocate memory
790 529 : call new%malloc()
791 529 : if (allocated(Cryst%irredatindx)) then
792 1563 : ABI_MALLOC(new%irredatindx,(new%nirredat))
793 : end if
794 :
795 : ! Copy dynamic arrays
796 198482 : new%indsym = Cryst%indsym
797 18946 : new%symafm = Cryst%symafm
798 233602 : new%symrec = Cryst%symrec
799 233602 : new%symrel = Cryst%symrel
800 233602 : new%symrel_cart = Cryst%symrel_cart
801 2392 : new%atindx = Cryst%atindx
802 2392 : new%atindx1 = Cryst%atindx1
803 2392 : new%typat = Cryst%typat
804 1944 : new%nattyp = Cryst%nattyp
805 72610 : new%tnons = Cryst%tnons
806 6394 : new%xcart = Cryst%xcart
807 6394 : new%xred = Cryst%xred
808 90498 : new%spinrot = Cryst%spinrot
809 1944 : new%amu = Cryst%amu
810 1944 : new%zion = Cryst%zion
811 1944 : new%znucl = Cryst%znucl
812 1944 : new%title = Cryst%title
813 529 : if (allocated(Cryst%irredatindx)) then
814 2054 : new%irredatindx = Cryst%irredatindx
815 : end if
816 :
817 529 : end subroutine crystal_copy
818 : !!***
819 :
820 : !----------------------------------------------------------------------
821 :
822 : !!****f* m_crystal/crystal_bcast
823 : !! NAME
824 : !! crystal_bcast
825 : !!
826 : !! FUNCTION
827 : !! Master broadcasts data and others allocate their arrays.
828 : !!
829 : !! SOURCE
830 :
831 24 : subroutine crystal_bcast(Cryst, comm)
832 :
833 : !Arguments ------------------------------------
834 : class(crystal_t),intent(inout) :: Cryst
835 : integer, intent(in) :: comm
836 :
837 : !Local variables -------------------------
838 : integer, parameter :: master=0
839 : integer :: ierr
840 : ! *********************************************************************
841 :
842 24 : if (xmpi_comm_size(comm) == 1) return
843 :
844 : DBG_ENTER("COLL")
845 :
846 : ! Integers
847 24 : call xmpi_bcast(Cryst%natom, master, comm, ierr)
848 24 : call xmpi_bcast(Cryst%nsym, master, comm, ierr)
849 24 : call xmpi_bcast(Cryst%ntypat, master, comm, ierr)
850 24 : call xmpi_bcast(Cryst%nirredat, master, comm, ierr)
851 24 : call xmpi_bcast(Cryst%npsp, master, comm, ierr)
852 24 : call xmpi_bcast(Cryst%space_group, master, comm, ierr)
853 24 : call xmpi_bcast(Cryst%timrev, master, comm, ierr)
854 24 : call xmpi_bcast(Cryst%use_antiferro, master, comm, ierr)
855 :
856 24 : if (xmpi_comm_rank(comm) /= master) then
857 16 : call Cryst%free()
858 16 : call Cryst%malloc()
859 : end if
860 :
861 : ! Floats
862 24 : call xmpi_bcast(Cryst%ucvol, master, comm, ierr)
863 :
864 : ! Arrays
865 24 : call xmpi_bcast(Cryst%angdeg, master, comm, ierr)
866 24 : call xmpi_bcast(Cryst%gmet, master, comm, ierr)
867 24 : call xmpi_bcast(Cryst%gprimd, master, comm, ierr)
868 24 : call xmpi_bcast(Cryst%rmet, master, comm, ierr)
869 24 : call xmpi_bcast(Cryst%rprimd, master, comm, ierr)
870 24 : call xmpi_bcast(Cryst%indsym, master, comm, ierr)
871 24 : call xmpi_bcast(Cryst%symafm, master, comm, ierr)
872 24 : call xmpi_bcast(Cryst%symrec, master, comm, ierr)
873 24 : call xmpi_bcast(Cryst%symrel, master, comm, ierr)
874 24 : call xmpi_bcast(Cryst%symrel_cart, master, comm, ierr)
875 24 : call xmpi_bcast(Cryst%atindx, master, comm, ierr)
876 24 : call xmpi_bcast(Cryst%atindx1, master, comm, ierr)
877 24 : call xmpi_bcast(Cryst%typat, master, comm, ierr)
878 24 : call xmpi_bcast(Cryst%nattyp, master, comm, ierr)
879 24 : call xmpi_bcast(Cryst%tnons, master, comm, ierr)
880 24 : call xmpi_bcast(Cryst%xcart, master, comm, ierr)
881 24 : call xmpi_bcast(Cryst%xred, master, comm, ierr)
882 24 : call xmpi_bcast(Cryst%spinrot, master, comm, ierr)
883 24 : call xmpi_bcast(Cryst%amu, master, comm, ierr)
884 24 : call xmpi_bcast(Cryst%zion, master, comm, ierr)
885 24 : call xmpi_bcast(Cryst%znucl, master, comm, ierr)
886 24 : call xmpi_bcast(Cryst%title, master, comm, ierr)
887 :
888 : ! It is not always allocated on master node,
889 : ! and it can be computed afterward on each node.
890 : !call xmpi_bcast(Cryst%irredatindx, master, comm, ierr)
891 :
892 : DBG_EXIT("COLL")
893 :
894 : end subroutine crystal_bcast
895 : !!***
896 :
897 : !----------------------------------------------------------------------
898 :
899 : !!****f* m_crystal/crystal_compare
900 : !! NAME
901 : !! crystal_compare
902 : !!
903 : !! FUNCTION
904 : !! Compare two crystalline structures,
905 : !! write warning messages to stdout if they differ, return exit status
906 : !!
907 : !! INPUTS
908 : !! [header]=Optional header message.
909 : !!
910 : !! OUTPUT
911 : !!
912 : !! SOURCE
913 :
914 258 : integer function crystal_compare(self, other, header) result(ierr)
915 :
916 : !Arguments ------------------------------------
917 : class(crystal_t),intent(in) :: self, other
918 : character(len=*),optional,intent(in) :: header
919 :
920 : !Local variables-------------------------------
921 : !integer :: isym, iat, itypat
922 : ! *********************************************************************
923 :
924 258 : if (present(header)) call wrtout(std_out, header)
925 258 : ierr = 0
926 :
927 : ! Test basic dimensions and metadata.
928 258 : ABI_CHECK_IEQ_IERR(self%natom, other%natom, "Different natom" , ierr)
929 258 : ABI_CHECK_IEQ_IERR(self%ntypat, other%ntypat, "Different ntypat" , ierr)
930 258 : ABI_CHECK_IEQ_IERR(self%npsp, other%npsp, "Different npsp" , ierr)
931 258 : ABI_CHECK_IEQ_IERR(self%nsym, other%nsym, "Different nsym" , ierr)
932 258 : ABI_CHECK_IEQ_IERR(self%timrev, other%timrev, "Different timrev" , ierr)
933 :
934 258 : if (ierr /= 0) goto 10
935 : ! After this point, we know that basic dimensions agree with each other.
936 : ! Yes, I use GOTO and I'm proud of that!
937 :
938 : ! Check direct lattice
939 3354 : if (any(abs(self%rprimd - other%rprimd) > tol6)) then
940 0 : ABI_WARNING("Found critical diffs in rprimd lattice vectors.")
941 0 : ierr = ierr + 1
942 : end if
943 :
944 : ! Check Symmetries
945 145234 : if (any(self%symrel /= other%symrel)) then
946 0 : ABI_WARNING("Found critical diffs in symrel symmetries.")
947 0 : ierr = ierr + 1
948 : end if
949 44866 : if (any(abs(self%tnons - other%tnons) > tol3)) then
950 0 : ABI_WARNING("Found critical diffs in fractional translations tnons.")
951 0 : ierr = ierr + 1
952 : end if
953 258 : if (self%use_antiferro .neqv. other%use_antiferro) then
954 0 : ABI_WARNING("Different values of use_antiferro")
955 0 : ierr = ierr + 1
956 : end if
957 :
958 : ! Atoms
959 754 : if (any(self%typat /= other%typat)) then
960 0 : ABI_WARNING("Found critical diffs in typat.")
961 : !print *, self%typat; print *, other%typat
962 0 : ierr = ierr + 1
963 : end if
964 648 : if (any(abs(self%zion - other%zion) > tol3)) then
965 0 : ABI_WARNING("Found critical diffs in zion.")
966 0 : ierr = ierr + 1
967 : end if
968 648 : if (any(abs(self%znucl - other%znucl) > tol3)) then
969 0 : ABI_WARNING("Found critical diffs in znucl.")
970 0 : ierr = ierr + 1
971 : end if
972 648 : if (any(abs(self%amu - other%amu) > tol3)) then
973 0 : ABI_WARNING("Found critical diffs in amu.")
974 0 : ierr = ierr + 1
975 : end if
976 2242 : if (any(abs(self%xred - other%xred) > tol6)) then
977 0 : ABI_WARNING("Found critical diffs in xred.")
978 0 : ierr = ierr + 1
979 : end if
980 :
981 258 : if (ierr /= 0) goto 10
982 0 : return
983 :
984 : ! Print structure to aid debugging. Caller will handle exit status.
985 0 : 10 call wrtout(std_out, " Comparing crystal1 and crystal2 for possible differences before returning ierr /= 0!")
986 0 : call self%print(header="crystal1")
987 0 : call wrtout(std_out, "")
988 0 : call other%print(header="crystal2")
989 0 : call wrtout(std_out, "")
990 :
991 0 : end function crystal_compare
992 : !!***
993 :
994 : !----------------------------------------------------------------------
995 :
996 : !!****f* m_crystal/crystal_print
997 : !! NAME
998 : !! crystal_print
999 : !!
1000 : !! FUNCTION
1001 : !! Print the content of crystal_t data type
1002 : !!
1003 : !! INPUTS
1004 : !! Cryst<crystal_t>=The structure.
1005 : !! [unit]=Unit number for output. Defaults to std_out
1006 : !! [prtvol]=Verbosity level. If prtvol== -1, only lattice parameters are printed. Defaults to 0
1007 : !! [mode_paral]=Either "COLL" or "PERS"
1008 : !! [header]=String to be printed as header for additional info.
1009 : !!
1010 : !! OUTPUT
1011 : !! Only printing
1012 : !!
1013 : !! SOURCE
1014 :
1015 1232 : subroutine crystal_print(Cryst, header, unit, mode_paral, prtvol)
1016 :
1017 : !Arguments ------------------------------------
1018 : !scalars
1019 : class(crystal_t),intent(in) :: Cryst
1020 : integer,optional,intent(in) :: unit, prtvol
1021 : character(len=*),optional,intent(in) :: mode_paral
1022 : character(len=*),optional,intent(in) :: header
1023 :
1024 : !Local variables-------------------------------
1025 : integer :: my_unt,my_prtvol,nu,iatom, isym, ii, nsym, units(1)
1026 : character(len=4) :: my_mode
1027 : character(len=500) :: msg
1028 : ! *********************************************************************
1029 :
1030 694 : my_unt =std_out; if (PRESENT(unit )) my_unt =unit
1031 1232 : my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol
1032 1232 : my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral
1033 :
1034 2464 : units = [my_unt]
1035 :
1036 1232 : msg=' ==== Info on the Cryst% object ==== '
1037 1232 : if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
1038 1232 : call wrtout(my_unt, sjoin(ch10, msg), my_mode)
1039 :
1040 1232 : write(msg,'(a)')' Real(R)+Recip(G) space primitive vectors, cartesian coordinates (Bohr,Bohr^-1):'
1041 1232 : call wrtout(my_unt,msg,my_mode)
1042 4928 : do nu=1,3
1043 : write(msg,'(1x,a,i1,a,3f11.7,2x,a,i1,a,3f11.7)')&
1044 14784 : 'R(',nu,')=',Cryst%rprimd(:,nu)+tol10, &
1045 18480 : 'G(',nu,')=',Cryst%gprimd(:,nu)+tol10 ! tol10 is used to be consistent with metric.F90
1046 4928 : call wrtout(my_unt,msg,my_mode)
1047 : end do
1048 :
1049 1232 : write(msg,'(a,1p,e15.7,a)')' Unit cell volume ucvol=',Cryst%ucvol+tol10,' bohr^3'
1050 1232 : call wrtout(my_unt,msg,my_mode)
1051 :
1052 1232 : write(msg,'(a,3es16.8,a)')' Angles (23,13,12)=',Cryst%angdeg(1:3),' degrees'
1053 1232 : call wrtout(my_unt,msg,my_mode)
1054 :
1055 1232 : if (Cryst%timrev==1) then
1056 198 : msg = ' Time-reversal symmetry is not present '
1057 1034 : else if (Cryst%timrev==2) then
1058 1034 : msg = ' Time-reversal symmetry is present '
1059 : else
1060 0 : ABI_BUG(sjoin('Wrong value for timrev:', itoa(cryst%timrev)))
1061 : end if
1062 1232 : call wrtout(my_unt,msg,my_mode)
1063 1232 : if (my_prtvol == -1) return
1064 :
1065 1229 : if (my_prtvol > 0) then
1066 0 : call print_symmetries(units, Cryst%nsym, Cryst%symrel, Cryst%tnons, Cryst%symafm)
1067 0 : if (Cryst%use_antiferro) call wrtout(my_unt,' System has magnetic symmetries ',my_mode)
1068 :
1069 : ! Print indsym using the same format as in symatm
1070 0 : nsym = cryst%nsym
1071 0 : do iatom=1,cryst%natom
1072 0 : write(msg, '(a,i0,a)' )' symatm: atom number ',iatom,' is reached starting at atom'
1073 0 : call wrtout(std_out, msg)
1074 0 : do ii=1,(nsym-1)/24+1
1075 0 : if (cryst%natom<100) then
1076 0 : write(msg, '(1x,24i3)' ) (cryst%indsym(4,isym,iatom),isym=1+(ii-1)*24,min(nsym,ii*24))
1077 : else
1078 0 : write(msg, '(1x,24i6)' ) (cryst%indsym(4,isym,iatom),isym=1+(ii-1)*24,min(nsym,ii*24))
1079 : end if
1080 0 : call wrtout(std_out, msg)
1081 : end do
1082 : end do
1083 :
1084 : end if
1085 :
1086 1229 : call wrtout(my_unt, " Reduced atomic positions [iatom, xred, symbol]:", my_mode)
1087 4242 : do iatom=1,cryst%natom
1088 3013 : write(msg,"(i5,a,2x,3f11.7,2x,a)")iatom,")",cryst%xred(:,iatom), cryst%symbol_type(cryst%typat(iatom))
1089 4242 : call wrtout(my_unt,msg,my_mode)
1090 : end do
1091 :
1092 1232 : end subroutine crystal_print
1093 : !!***
1094 :
1095 : !!****f* m_crystal/crystal_print_abivars
1096 : !! NAME
1097 : !! crystal_print_abivars
1098 : !!
1099 : !! FUNCTION
1100 : !! Print unit cell info in Abinit/abivars format
1101 : !!
1102 : !! INPUTS
1103 : !! unit=Output unit
1104 : !!
1105 : !! OUTPUT
1106 : !! Only printing
1107 : !!
1108 : !! SOURCE
1109 :
1110 0 : subroutine crystal_print_abivars(cryst, unit)
1111 :
1112 : !Arguments ------------------------------------
1113 : !scalars
1114 : class(crystal_t),intent(in) :: cryst
1115 : integer,intent(in) :: unit
1116 :
1117 : !Local variables-------------------------------
1118 : integer :: iatom, ii
1119 : !character(len=500) :: fmt
1120 : ! *********************************************************************
1121 :
1122 0 : if (unit == dev_null) return
1123 :
1124 : ! Write variables using standard Abinit input format.
1125 0 : write(unit, "(/,/,a)")" # Abinit variables"
1126 0 : write(unit, "(a)")" acell 1.0 1.0 1.0"
1127 0 : write(unit, "(a)")" rprimd"
1128 0 : do ii=1,3
1129 0 : write(unit, "(3(f11.7,1x))")cryst%rprimd(:, ii)
1130 : end do
1131 0 : write(unit, "(a, i0)")" natom ", cryst%natom
1132 0 : write(unit, "(a, i0)")" ntypat ", cryst%ntypat
1133 0 : write(unit, strcat("(a, ", itoa(cryst%natom), "(i0,1x))")) " typat ", cryst%typat
1134 0 : write(unit, strcat("(a, ", itoa(cryst%npsp), "(f5.1,1x))")) " znucl ", cryst%znucl
1135 0 : write(unit, "(a)")" xred"
1136 0 : do iatom=1,cryst%natom
1137 0 : write(unit,"(1x, 3f11.7,2x,2a)")cryst%xred(:,iatom), " # ", cryst%symbol_type(cryst%typat(iatom))
1138 : end do
1139 :
1140 : ! Write variables using the abivars format supported by structure variable.
1141 : !write(unit, "(/,/,a)")" # Abivars format (external file with structure variable)"
1142 : !write(unit, "(a)")" acell 1.0 1.0 1.0"
1143 : !write(unit, "(a)")" rprimd"
1144 : !do ii=1,3
1145 : ! write(unit, "(1x, 3(f11.7,1x))")cryst%rprimd(:, ii)
1146 : !end do
1147 : !write(unit, "(a, i0)")" natom ", cryst%natom
1148 : !write(unit, "(a)")" xred_symbols"
1149 : !do iatom=1,cryst%natom
1150 : ! write(unit,"(1x, 3f11.7,2x,a)")cryst%xred(:,iatom), cryst%symbol_type(cryst%typat(iatom))
1151 : !end do
1152 :
1153 : end subroutine crystal_print_abivars
1154 : !!***
1155 :
1156 : !----------------------------------------------------------------------
1157 :
1158 : !!****f* m_crystal/symbols_crystal
1159 : !!
1160 : !! NAME
1161 : !! symbols_crystal
1162 : !!
1163 : !! FUNCTION
1164 : !! Return a array with the symbol of each atoms with indexation e.g.
1165 : !! ["Sr","Ru","O1","O2","O3"]
1166 : !!
1167 : !! INPUTS
1168 : !! natom = number of atoms
1169 : !! ntypat = number of typat
1170 : !! npsp = number of pseudopotentials
1171 : !! znucl = Nuclear charge for each type of pseudopotential
1172 : !!
1173 : !! OUTPUT
1174 : !! symbols = array with the symbol of each atoms
1175 : !!
1176 : !! SOURCE
1177 :
1178 104 : subroutine symbols_crystal(natom, ntypat, npsp, symbols, typat, znucl)
1179 :
1180 : !Arguments ------------------------------------
1181 : !scalars
1182 : integer,intent(in) :: natom,ntypat,npsp
1183 : !arrays
1184 : real(dp),intent(in):: znucl(npsp)
1185 : integer,intent(in) :: typat(natom)
1186 : character(len=5),intent(out) :: symbols(natom)
1187 : character(len=3) :: powerchar
1188 :
1189 : !Local variables-------------------------------
1190 : !scalar
1191 : integer :: ia,ii,itypat,jj
1192 : ! *************************************************************************
1193 :
1194 : ! Fill the symbols array
1195 624 : do ia=1,natom
1196 624 : symbols(ia) = adjustl(znucl2symbol(znucl(typat(ia))))
1197 : end do
1198 104 : itypat = 0
1199 416 : do itypat =1,ntypat
1200 : ii = 0
1201 1872 : do ia=1,natom
1202 1872 : if(typat(ia)==itypat) then
1203 520 : ii = ii + 1
1204 : end if
1205 : end do
1206 416 : if(ii>1)then
1207 104 : jj=1
1208 624 : do ia=1,natom
1209 624 : if(typat(ia)==itypat) then
1210 312 : write(powerchar,'(I0)') jj
1211 312 : symbols(ia) = trim(symbols(ia))//trim(powerchar)
1212 312 : jj=jj+1
1213 : end if
1214 : end do
1215 : end if
1216 : end do
1217 :
1218 104 : end subroutine symbols_crystal
1219 : !!***
1220 :
1221 : !----------------------------------------------------------------------
1222 :
1223 : !!****f* m_crystal/idx_spatial_inversion
1224 : !! NAME
1225 : !! idx_spatial_inversion
1226 : !!
1227 : !! FUNCTION
1228 : !! Return the index of the spatial inversion, 0 if not present
1229 : !!
1230 : !! SOURCE
1231 :
1232 18 : pure integer function idx_spatial_inversion(Cryst) result(inv_idx)
1233 :
1234 : !Arguments ------------------------------------
1235 : class(crystal_t),intent(in) :: Cryst
1236 :
1237 : !Local variables-------------------------------
1238 : integer :: isym
1239 : ! *************************************************************************
1240 :
1241 18 : inv_idx=0
1242 427 : do isym=1,cryst%nsym
1243 677 : if (all(cryst%symrel(:,:,isym) == inversion_3d)) then
1244 18 : inv_idx=isym; return
1245 : end if
1246 : end do
1247 :
1248 : end function idx_spatial_inversion
1249 : !!***
1250 :
1251 : !----------------------------------------------------------------------
1252 :
1253 : !!****f* m_crystal/isymmorphic
1254 : !! NAME
1255 : !! isymmorphic
1256 : !!
1257 : !! FUNCTION
1258 : !! Returns .TRUE. if space group is symmorphic, i.e. all fractional translations are zero.
1259 : !!
1260 : !! SOURCE
1261 :
1262 18636 : pure logical function isymmorphic(Cryst) result(ans)
1263 :
1264 : !Arguments ------------------------------------
1265 : class(crystal_t),intent(in) :: Cryst
1266 : ! *************************************************************************
1267 :
1268 1030544 : ans = ALL(ABS(Cryst%tnons) < tol6)
1269 :
1270 18636 : end function isymmorphic
1271 : !!***
1272 :
1273 : !----------------------------------------------------------------------
1274 :
1275 : !!****f* m_crystal/isalchemical
1276 : !! NAME
1277 : !! isalchemical
1278 : !!
1279 : !! FUNCTION
1280 : !! Returns .TRUE. if we are using alchemical pseudopotentials
1281 : !!
1282 : !! SOURCE
1283 :
1284 18636 : pure logical function isalchemical(Cryst) result(ans)
1285 :
1286 : !Arguments ------------------------------------
1287 : class(crystal_t),intent(in) :: Cryst
1288 : ! *************************************************************************
1289 :
1290 18636 : ans = (Cryst%npsp /= Cryst%ntypat)
1291 :
1292 18636 : end function isalchemical
1293 : !!***
1294 :
1295 : !----------------------------------------------------------------------
1296 :
1297 : !!****f* m_crystal/adata_type
1298 : !! NAME
1299 : !! adata_type
1300 : !!
1301 : !! FUNCTION
1302 : !! Return atomic data from the itypat index
1303 : !!
1304 : !! SOURCE
1305 :
1306 18574 : type(atomdata_t) function adata_type(crystal, itypat) result(atom)
1307 :
1308 : !Arguments ------------------------------------
1309 : class(crystal_t),intent(in) :: crystal
1310 : integer,intent(in) :: itypat
1311 : ! *************************************************************************
1312 :
1313 18574 : call atomdata_from_znucl(atom, crystal%znucl(itypat))
1314 :
1315 18574 : end function adata_type
1316 : !!***
1317 :
1318 : !----------------------------------------------------------------------
1319 :
1320 : !!****f* m_crystal/symbol_type
1321 : !! NAME
1322 : !! symbol_type
1323 : !!
1324 : !! FUNCTION
1325 : !! Return the atomic symbol from the itypat index
1326 : !!
1327 : !! SOURCE
1328 :
1329 18574 : function symbol_type(crystal, itypat) result(symbol)
1330 :
1331 : !Arguments ------------------------------------
1332 : !scalars
1333 : class(crystal_t),intent(in) :: crystal
1334 : integer,intent(in) :: itypat
1335 : character(len=2) :: symbol
1336 :
1337 : !Local variables-------------------------------
1338 : type(atomdata_t) :: atom
1339 : ! *************************************************************************
1340 :
1341 18574 : atom = crystal%adata_type(itypat)
1342 18574 : symbol = atom%symbol
1343 :
1344 18574 : end function symbol_type
1345 : !!***
1346 :
1347 : !----------------------------------------------------------------------
1348 :
1349 : !!****f* m_crystal/symbol_iatom
1350 : !! NAME
1351 : !! symbol_iatom
1352 : !!
1353 : !! FUNCTION
1354 : !! Return the atomic symbol from the iatom index
1355 : !!
1356 : !! SOURCE
1357 :
1358 15561 : function symbol_iatom(crystal, iatom) result(symbol)
1359 :
1360 : !Arguments ------------------------------------
1361 : !scalars
1362 : class(crystal_t),intent(in) :: crystal
1363 : integer,intent(in) :: iatom
1364 : character(len=2) :: symbol
1365 : ! *************************************************************************
1366 :
1367 15561 : symbol = crystal%symbol_type(crystal%typat(iatom))
1368 :
1369 15561 : end function symbol_iatom
1370 : !!***
1371 :
1372 : !----------------------------------------------------------------------
1373 :
1374 : !!****f* m_crystal/crystal_point_group
1375 : !! NAME
1376 : !! crystal_point_group
1377 : !!
1378 : !! FUNCTION
1379 : !! Return the symmetries of the point group of the crystal.
1380 : !!
1381 : !! INPUTS
1382 : !! [include_timrev]=If True, time-reversal symmetry is included in the point group unless
1383 : !! the system has spatial inversion. Default: False
1384 : !!
1385 : !! OUTPUT
1386 : !! ptg_nsym=Number of symmetries in the point group
1387 : !! ptg_symrel(3,3,ptg_nsym)=Rotations in real space
1388 : !! ptg_symrec(3,3,ptg_nsym)=Rotations in reciprocal space
1389 : !! has_inversion=True if spatial inversion is present in the point group.
1390 : !!
1391 : !! SOURCE
1392 :
1393 6 : subroutine crystal_point_group(cryst, ptg_nsym, ptg_symrel, ptg_symrec, has_inversion, include_timrev)
1394 :
1395 : !Arguments ------------------------------------
1396 : !scalars
1397 : class(crystal_t),intent(in) :: cryst
1398 : integer,intent(out) :: ptg_nsym
1399 : logical,intent(out) :: has_inversion
1400 : logical,optional,intent(in) :: include_timrev
1401 : !arrays
1402 : integer,allocatable,intent(out) :: ptg_symrel(:,:,:),ptg_symrec(:,:,:)
1403 :
1404 : !Local variables-------------------------------
1405 : !scalars
1406 : integer :: isym, search, tmp_nsym, ierr
1407 : logical :: found, my_include_timrev, debug
1408 : !arrays
1409 12 : integer :: work_symrel(3,3,cryst%nsym)
1410 6 : integer,allocatable :: symafm(:)
1411 : ! *************************************************************************
1412 :
1413 6 : my_include_timrev = .False.; if (present(include_timrev)) my_include_timrev = include_timrev
1414 :
1415 78 : tmp_nsym = 1; work_symrel(:,:,1) = cryst%symrel(:,:,1)
1416 288 : do isym=2,cryst%nsym
1417 282 : if (cryst%symafm(isym) == -1) cycle
1418 7050 : do search=1,tmp_nsym
1419 11124 : found = all(work_symrel(:,:,search) == cryst%symrel(:,:,isym))
1420 7050 : if (found) exit
1421 : end do
1422 288 : if (.not. found) then
1423 282 : tmp_nsym = tmp_nsym + 1
1424 3666 : work_symrel(:,:,tmp_nsym) = cryst%symrel(:,:,isym)
1425 : end if
1426 : end do
1427 :
1428 6 : has_inversion = .False.
1429 12 : do isym=1,tmp_nsym
1430 84 : if (all(work_symrel(:,:,isym) == inversion_3d) ) then
1431 6 : has_inversion = .True.; exit
1432 : end if
1433 : end do
1434 :
1435 : ! Now we know the symmetries of the point group.
1436 6 : ptg_nsym = tmp_nsym; if (.not. has_inversion .and. my_include_timrev) ptg_nsym = 2 * tmp_nsym
1437 18 : ABI_MALLOC(ptg_symrel, (3, 3, ptg_nsym))
1438 12 : ABI_MALLOC(ptg_symrec, (3, 3, ptg_nsym))
1439 :
1440 3750 : ptg_symrel(:,:,1:tmp_nsym) = work_symrel(:,:,1:tmp_nsym)
1441 294 : do isym=1,tmp_nsym
1442 294 : call mati3inv(ptg_symrel(:,:,isym), ptg_symrec(:,:,isym))
1443 : end do
1444 :
1445 6 : if (.not. has_inversion .and. my_include_timrev) then
1446 0 : ptg_symrel(:,:,tmp_nsym+1:) = -work_symrel(:,:,1:tmp_nsym)
1447 0 : do isym=tmp_nsym+1,ptg_nsym
1448 0 : call mati3inv(ptg_symrel(:,:,isym), ptg_symrec(:,:,isym))
1449 : end do
1450 : end if
1451 :
1452 6 : debug = .False.
1453 : if (debug) then
1454 : ABI_MALLOC(symafm, (ptg_nsym))
1455 : symafm = 1
1456 : call sg_multable(ptg_nsym, symafm, ptg_symrel, ierr)
1457 : ABI_CHECK(ierr == 0, "point group is not a group! See messages above")
1458 : ABI_FREE(symafm)
1459 : end if
1460 :
1461 6 : end subroutine crystal_point_group
1462 : !!***
1463 :
1464 : !----------------------------------------------------------------------
1465 :
1466 : !!****f* m_crystal/crystal_ncwrite
1467 : !! NAME
1468 : !! crystal_ncwrite
1469 : !!
1470 : !! FUNCTION
1471 : !! Output system geometry to a file, using the NETCDF file format and ETSF I/O.
1472 : !! Data are taken from the crystal_t object.
1473 : !!
1474 : !! INPUTS
1475 : !! cryst<crystal_t>=Object defining the unit cell and its symmetries.
1476 : !! ncid=NC file handle.
1477 : !!
1478 : !! OUTPUT
1479 : !! Only writing
1480 : !!
1481 : !! NOTES
1482 : !! Alchemy not treated, since crystal should be initialized at the beginning of the run.
1483 : !!
1484 : !! SOURCE
1485 :
1486 18636 : integer function crystal_ncwrite(cryst, ncid) result(ncerr)
1487 :
1488 : !Arguments ------------------------------------
1489 : !scalars
1490 : class(crystal_t),intent(in) :: cryst
1491 : integer,intent(in) :: ncid
1492 :
1493 : !Local variables-------------------------------
1494 : !scalars
1495 : integer :: itypat
1496 : character(len=500) :: msg
1497 : character(len=etsfio_charlen) :: symmorphic
1498 : type(atomdata_t) :: atom
1499 : !arrays
1500 37272 : character(len=2) :: symbols(cryst%ntypat)
1501 37272 : character(len=80) :: psp_desc(cryst%ntypat),symbols_long(cryst%ntypat)
1502 : ! *************************************************************************
1503 :
1504 : ! TODO alchemy not treated correctly by ETSF_IO specs.
1505 18636 : if (cryst%isalchemical()) then
1506 : write(msg,"(3a)")&
1507 70 : "Alchemical crystals are not fully supported by the netcdf format",ch10,&
1508 140 : "Important parameters (e.g. znucl, symbols) are not written with the correct value"
1509 70 : ABI_WARNING(msg)
1510 : end if
1511 :
1512 18636 : symmorphic = yesno(cryst%isymmorphic())
1513 :
1514 : ! Define dimensions.
1515 : ! npsp added in v9.
1516 : ncerr = nctk_def_dims(ncid, [ &
1517 : nctkdim_t("complex", 2), nctkdim_t("symbol_length", 2),&
1518 : nctkdim_t("character_string_length", 80), nctkdim_t("number_of_cartesian_directions", 3),&
1519 : nctkdim_t("number_of_reduced_dimensions", 3), nctkdim_t("number_of_vectors", 3),&
1520 : nctkdim_t("number_of_atoms", cryst%natom), nctkdim_t("number_of_atom_species", cryst%ntypat),&
1521 : nctkdim_t("number_of_atom_pseudopotentials", cryst%npsp),&
1522 204996 : nctkdim_t("number_of_symmetry_operations", cryst%nsym)], defmode=.True.)
1523 18636 : NCF_CHECK(ncerr)
1524 :
1525 : ! Define variables
1526 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
1527 74544 : "space_group", "time_reversal", "use_antiferromagnetic_symmetries"])
1528 18636 : NCF_CHECK(ncerr)
1529 :
1530 : ncerr = nctk_def_arrays(ncid, [ &
1531 : ! Atomic structure and symmetry operations
1532 : nctkarr_t("primitive_vectors", "dp", "number_of_cartesian_directions, number_of_vectors"), &
1533 : nctkarr_t("reduced_symmetry_matrices", "int", &
1534 : "number_of_reduced_dimensions, number_of_reduced_dimensions, number_of_symmetry_operations"), &
1535 : nctkarr_t("reduced_symmetry_translations", "dp", "number_of_reduced_dimensions, number_of_symmetry_operations"), &
1536 : nctkarr_t("atom_species", "int", "number_of_atoms"), &
1537 : nctkarr_t("reduced_atom_positions", "dp", "number_of_reduced_dimensions, number_of_atoms"), &
1538 : nctkarr_t("atomic_numbers", "dp", "number_of_atom_species"), &
1539 : nctkarr_t("atom_species_names", "char", "character_string_length, number_of_atom_species"), &
1540 : nctkarr_t("chemical_symbols", "char", "symbol_length, number_of_atom_species"), &
1541 : nctkarr_t('atomic_mass_units', "dp", "number_of_atom_species"), &
1542 : ! Atomic information.
1543 : nctkarr_t("valence_charges", "dp", "number_of_atom_species"), & ! NB: This variable is not written if alchemical
1544 : nctkarr_t("pseudopotential_types", "char", "character_string_length, number_of_atom_species") &
1545 223632 : ])
1546 18636 : NCF_CHECK(ncerr)
1547 :
1548 : ! Some variables require the "symmorphic" attribute.
1549 18636 : NCF_CHECK(nf90_put_att(ncid, vid("reduced_symmetry_matrices"), "symmorphic", symmorphic))
1550 18636 : NCF_CHECK(nf90_put_att(ncid, vid("reduced_symmetry_translations"), "symmorphic", symmorphic))
1551 :
1552 : ! At this point we have an ETSF-compliant file. Add additional data for internal use in abinit.
1553 : ncerr = nctk_def_arrays(ncid, [ &
1554 : nctkarr_t('symafm', "int", "number_of_symmetry_operations"), &
1555 : nctkarr_t('symrel_cart', "dp", "three, three, number_of_symmetry_operations"), &
1556 : nctkarr_t('indsym', "int", "four, number_of_symmetry_operations, number_of_atoms") &
1557 74544 : ])
1558 18636 : NCF_CHECK(ncerr)
1559 :
1560 : ! Set-up atomic symbols.
1561 42668 : do itypat=1,cryst%ntypat
1562 24032 : call atomdata_from_znucl(atom, cryst%znucl(itypat))
1563 24032 : symbols(itypat) = atom%symbol
1564 24032 : write(symbols_long(itypat),'(a2,a78)') symbols(itypat),REPEAT(CHAR(0),78)
1565 : write(psp_desc(itypat),'(2a)') &
1566 790344 : cryst%title(itypat)(1:MIN(80,LEN_TRIM(cryst%title(itypat)))),REPEAT(CHAR(0),MAX(0,80-LEN_TRIM(cryst%title(itypat))))
1567 : end do
1568 :
1569 : ! Write data.
1570 18636 : NCF_CHECK(nctk_set_datamode(ncid))
1571 18636 : NCF_CHECK(nf90_put_var(ncid, vid("space_group"), cryst%space_group))
1572 18636 : NCF_CHECK(nf90_put_var(ncid, vid("primitive_vectors"), cryst%rprimd))
1573 18636 : NCF_CHECK(nf90_put_var(ncid, vid("reduced_symmetry_matrices"), cryst%symrel))
1574 18636 : NCF_CHECK(nf90_put_var(ncid, vid("reduced_symmetry_translations"), cryst%tnons))
1575 18636 : NCF_CHECK(nf90_put_var(ncid, vid("atom_species"), cryst%typat))
1576 18636 : NCF_CHECK(nf90_put_var(ncid, vid("reduced_atom_positions"), cryst%xred))
1577 18636 : NCF_CHECK(nf90_put_var(ncid, vid("atomic_numbers"), cryst%znucl(1:cryst%ntypat)))
1578 18636 : NCF_CHECK(nf90_put_var(ncid, vid("atom_species_names"), symbols_long))
1579 18636 : NCF_CHECK(nf90_put_var(ncid, vid("chemical_symbols"), symbols))
1580 18636 : NCF_CHECK(nf90_put_var(ncid, vid('atomic_mass_units'), cryst%amu))
1581 18636 : NCF_CHECK(nf90_put_var(ncid, vid("pseudopotential_types"), psp_desc))
1582 18636 : if (cryst%npsp == cryst%ntypat) then
1583 18566 : NCF_CHECK(nf90_put_var(ncid, vid("valence_charges"), cryst%zion))
1584 : end if
1585 :
1586 18636 : NCF_CHECK(nf90_put_var(ncid, vid("symafm"), cryst%symafm))
1587 18636 : NCF_CHECK(nf90_put_var(ncid, vid("symrel_cart"), cryst%symrel_cart))
1588 18636 : NCF_CHECK(nf90_put_var(ncid, vid("indsym"), cryst%indsym))
1589 :
1590 : ! Variables pertaining to the symmetry of the wavefunctions.
1591 : ! Note that these variables will be used in crystal_compare
1592 18636 : NCF_CHECK(nf90_put_var(ncid, vid("time_reversal"), cryst%timrev))
1593 :
1594 18636 : if (cryst%use_antiferro) then
1595 119 : NCF_CHECK(nf90_put_var(ncid, vid("use_antiferromagnetic_symmetries"), 1))
1596 : else
1597 18517 : NCF_CHECK(nf90_put_var(ncid, vid("use_antiferromagnetic_symmetries"), 0))
1598 : end if
1599 :
1600 : contains
1601 354014 : integer function vid(vname)
1602 : character(len=*),intent(in) :: vname
1603 354014 : vid = nctk_idname(ncid, vname)
1604 : end function vid
1605 :
1606 : end function crystal_ncwrite
1607 : !!***
1608 :
1609 : !----------------------------------------------------------------------
1610 :
1611 : !!****f* m_crystal/crystal_ncwrite_path
1612 : !! NAME
1613 : !! crystal_ncwrite_path
1614 : !!
1615 : !! FUNCTION
1616 : !! Output system geometry to a file, using the NETCDF file format and ETSF I/O.
1617 : !!
1618 : !! INPUTS
1619 : !! crystal<crystal_t>=Object defining the unit cell and its symmetries.
1620 : !! path=filename
1621 : !!
1622 : !! OUTPUT
1623 : !! Only writing
1624 : !!
1625 : !! SOURCE
1626 :
1627 0 : integer function crystal_ncwrite_path(crystal, path) result(ncerr)
1628 :
1629 : !Arguments ------------------------------------
1630 : !scalars
1631 : character(len=*),intent(in) :: path
1632 : class(crystal_t),intent(in) :: crystal
1633 :
1634 : !Local variables-------------------------------
1635 : integer :: ncid
1636 : ! *************************************************************************
1637 :
1638 0 : ncerr = nf90_noerr
1639 0 : if (file_exists(path)) then
1640 0 : NCF_CHECK(nctk_open_modify(ncid, path, xmpi_comm_self))
1641 : else
1642 0 : ncerr = nctk_open_create(ncid, path, xmpi_comm_self)
1643 0 : NCF_CHECK_MSG(ncerr, sjoin("creating:", path))
1644 : end if
1645 :
1646 0 : NCF_CHECK(crystal_ncwrite(crystal, ncid))
1647 0 : NCF_CHECK(nf90_close(ncid))
1648 :
1649 0 : end function crystal_ncwrite_path
1650 : !!***
1651 :
1652 : !----------------------------------------------------------------------
1653 :
1654 : !!****f* m_crystal/crystal_ncread
1655 : !! NAME
1656 : !! crystal_ncread
1657 : !!
1658 : !! FUNCTION
1659 : !! Read the crystal object from a NETCDF file.
1660 : !!
1661 : !! INPUTS
1662 : !! cryst<crystal_t>=Object defining the unit cell and its symmetries.
1663 : !! ncid=NC file handle.
1664 : !!
1665 : !! OUTPUT
1666 : !! crystal
1667 : !!
1668 : !! SOURCE
1669 :
1670 401 : subroutine crystal_ncread(cryst, ncid)
1671 :
1672 : !Arguments ------------------------------------
1673 : !scalars
1674 : class(crystal_t),intent(inout) :: cryst
1675 : integer,intent(in) :: ncid
1676 :
1677 : !Local variables ------------------------------------
1678 : integer :: use_antiferro
1679 : ! *************************************************************************
1680 :
1681 : ! ---------------
1682 : ! Read dimensions
1683 : ! ---------------
1684 401 : NCF_CHECK(nctk_get_dim(ncid, "number_of_atoms", cryst%natom))
1685 401 : NCF_CHECK(nctk_get_dim(ncid, "number_of_atom_species", cryst%ntypat))
1686 401 : NCF_CHECK(nctk_get_dim(ncid, "number_of_atom_pseudopotentials", cryst%npsp))
1687 401 : NCF_CHECK(nctk_get_dim(ncid, "number_of_symmetry_operations", cryst%nsym))
1688 :
1689 : ! ---------------
1690 : ! Allocate memory
1691 : ! ---------------
1692 401 : call cryst%malloc()
1693 :
1694 : ! ------------
1695 : ! read scalars
1696 : ! ------------
1697 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "space_group"), cryst%space_group))
1698 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "time_reversal"), cryst%timrev))
1699 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "use_antiferromagnetic_symmetries"), use_antiferro))
1700 401 : cryst%use_antiferro = .False.
1701 401 : if (use_antiferro /= 0) cryst%use_antiferro = .True.
1702 :
1703 : ! -----------
1704 : ! read arrays
1705 : ! -----------
1706 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "primitive_vectors"), cryst%rprimd))
1707 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "reduced_symmetry_matrices"), cryst%symrel))
1708 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "reduced_symmetry_translations"), cryst%tnons))
1709 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "atom_species"), cryst%typat))
1710 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "reduced_atom_positions"), cryst%xred))
1711 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "atomic_numbers"), cryst%znucl(1:cryst%ntypat)))
1712 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "atomic_mass_units"), cryst%amu))
1713 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "symafm"), cryst%symafm))
1714 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "symrel_cart"), cryst%symrel_cart))
1715 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "indsym"), cryst%indsym))
1716 :
1717 401 : if (cryst%npsp == cryst%ntypat) then
1718 401 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "valence_charges"), cryst%zion))
1719 : end if
1720 :
1721 : ! Ignore those
1722 : !NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "pseudopotential_types"), psp_desc))
1723 : !NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "atom_species_names"), symbols_long))
1724 : !NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "chemical_symbols"), symbols))
1725 :
1726 : ! -----------------------
1727 : ! Complete initialization
1728 : ! -----------------------
1729 :
1730 401 : call cryst%compute_geometry()
1731 401 : call cryst%index_atoms()
1732 401 : call cryst%compute_sym()
1733 :
1734 401 : end subroutine crystal_ncread
1735 : !!***
1736 :
1737 : !----------------------------------------------------------------------
1738 :
1739 : !!****f* m_crystal/prt_cif
1740 : !! NAME
1741 : !! prt_cif
1742 : !!
1743 : !! FUNCTION
1744 : !! print out CIF format file
1745 : !!
1746 : !! INPUTS
1747 : !!
1748 : !! OUTPUT
1749 : !!
1750 : !! SOURCE
1751 :
1752 4 : subroutine prt_cif(brvltt, ciffname, natom, nsym, ntypat, rprimd, &
1753 2 : spgaxor, spgroup, spgorig, symrel, tnon, typat, xred, znucl)
1754 :
1755 : !Arguments ------------------------------------
1756 : !scalars
1757 : integer,intent(in) :: natom, ntypat, nsym
1758 : integer, intent(in) :: brvltt, spgaxor, spgroup, spgorig
1759 : !arrays
1760 : integer, intent(in) :: typat(natom)
1761 : integer, intent(in) :: symrel(3,3,nsym)
1762 : character(len=*), intent(in) :: ciffname
1763 : real(dp), intent(in) :: tnon(3,nsym)
1764 : real(dp), intent(in) :: rprimd(3,3)
1765 : real(dp), intent(in) :: xred(3,natom)
1766 : real(dp), intent(in) :: znucl(ntypat)
1767 :
1768 : !Local variables -------------------------------
1769 : !scalars
1770 : integer :: unitcif, iatom, isym, sporder, itypat, nat_this_type
1771 : real(dp) :: ucvol
1772 : type(atomdata_t) :: atom
1773 : !arrays
1774 : character(len=80) :: tmpstring
1775 : character(len=1) :: brvsb
1776 : character(len=15) :: intsb,ptintsb,ptschsb,schsb
1777 : character(len=35) :: intsbl
1778 : character(len=10) :: str_nat_type
1779 : character(len=100) :: chemformula
1780 : character(len=500) :: msg
1781 : real(dp) :: angle(3), gprimd(3,3), rmet(3,3), gmet(3,3)
1782 :
1783 : !*************************************************************************
1784 :
1785 : ! open file in append mode xlf and other compilers refuse append mode
1786 2 : if (open_file(ciffname,msg,newunit=unitcif) /=0) then
1787 0 : ABI_WARNING(msg)
1788 0 : return
1789 : end if
1790 :
1791 : ! print title for dataset
1792 2 : write (unitcif,'(a)') 'data_set'
1793 :
1794 : ! print cell parameters a,b,c, angles, volume
1795 2 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1796 2 : angle(1)=acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3)))/two_pi*360.0_dp
1797 2 : angle(2)=acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3)))/two_pi*360.0_dp
1798 2 : angle(3)=acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2)))/two_pi*360.0_dp
1799 :
1800 2 : write (unitcif,'(a,E20.10)') '_cell_length_a ', sqrt(rmet(1,1))*Bohr_Ang
1801 2 : write (unitcif,'(a,E20.10)') '_cell_length_b ', sqrt(rmet(2,2))*Bohr_Ang
1802 2 : write (unitcif,'(a,E20.10)') '_cell_length_c ', sqrt(rmet(3,3))*Bohr_Ang
1803 2 : write (unitcif,'(a,E20.10)') '_cell_angle_alpha ', angle(1)
1804 2 : write (unitcif,'(a,E20.10)') '_cell_angle_beta ', angle(2)
1805 2 : write (unitcif,'(a,E20.10)') '_cell_angle_gamma ', angle(3)
1806 2 : write (unitcif,'(a,E20.10)') '_cell_volume ', ucvol*(Bohr_Ang)**3
1807 :
1808 : ! print reduced positions
1809 2 : write (unitcif,'(a)') 'loop_'
1810 2 : write (unitcif,'(a,E20.10)') ' _atom_site_label '
1811 2 : write (unitcif,'(a,E20.10)') ' _atom_site_fract_x '
1812 2 : write (unitcif,'(a,E20.10)') ' _atom_site_fract_y '
1813 2 : write (unitcif,'(a,E20.10)') ' _atom_site_fract_z '
1814 5 : do iatom = 1, natom
1815 3 : call atomdata_from_znucl(atom,znucl(typat(iatom)))
1816 5 : write (unitcif,'(2a,3E20.10)') ' ', atom%symbol, xred(:,iatom)
1817 : end do
1818 :
1819 : !other specs in CIF dictionary which may be useful:
1820 : !GEOM_BOND GEOM_ANGLE GEOM_TORSION
1821 :
1822 : ! print chemical composition in simplest form
1823 2 : chemformula = "'"
1824 4 : do itypat = 1, ntypat
1825 2 : nat_this_type = 0
1826 5 : do iatom = 1, natom
1827 5 : if (typat(iatom) == itypat) nat_this_type = nat_this_type+1
1828 : end do
1829 2 : call atomdata_from_znucl(atom,znucl(itypat))
1830 2 : call int2char10(nat_this_type, str_nat_type)
1831 6 : chemformula = trim(chemformula) // atom%symbol // trim(str_nat_type) // " "
1832 : end do
1833 2 : chemformula = trim(chemformula) // "'"
1834 2 : write (unitcif,'(2a)') '_chemical_formula_analytical ', chemformula
1835 :
1836 : !FIXME: check that brvltt is correctly used here - is it equal to bravais(1) in the invars routines?
1837 2 : if (brvltt==1) then
1838 0 : write (unitcif,'(a)') '_symmetry_cell_setting triclinic'
1839 : else if(brvltt==2)then
1840 0 : write (unitcif,'(a)') '_symmetry_cell_setting monoclinic'
1841 : else if(brvltt==3)then
1842 0 : write (unitcif,'(a)') '_symmetry_cell_setting orthorhombic'
1843 : else if(brvltt==4)then
1844 0 : write (unitcif,'(a)') '_symmetry_cell_setting tetragonal'
1845 : else if(brvltt==5)then
1846 0 : write (unitcif,'(a)') '_symmetry_cell_setting rhombohedral'
1847 : else if(brvltt==6)then
1848 0 : write (unitcif,'(a)') '_symmetry_cell_setting hexagonal'
1849 : else if(brvltt==7)then
1850 0 : write (unitcif,'(a)') '_symmetry_cell_setting cubic'
1851 : end if
1852 :
1853 2 : call spgdata(brvsb,intsb,intsbl,ptintsb,ptschsb,schsb,spgaxor,spgroup,sporder,spgorig)
1854 :
1855 : ! print symmetry operations
1856 2 : write (unitcif,'(a,I6)') "_symmetry_Int_Tables_number ", spgroup
1857 2 : write (unitcif,'(5a)') "_symmetry_space_group_name_H-M '", brvsb, " ", trim(intsb), "'"
1858 2 : write (unitcif,'(a)') ''
1859 2 : write (unitcif,'(a)') 'loop_'
1860 2 : write (unitcif,'(a)') ' _symmetry_equiv_pos_as_xyz '
1861 74 : do isym = 1, nsym
1862 72 : call symrel2string(symrel(:,:,isym), tnon(:,isym), tmpstring)
1863 74 : write (unitcif,'(2a)') ' ', trim(tmpstring)
1864 : end do
1865 :
1866 2 : close(unitcif)
1867 :
1868 : end subroutine prt_cif
1869 : !!***
1870 :
1871 : !!****f* m_crystal/symrel2string
1872 : !! NAME
1873 : !! symrel2string
1874 : !!
1875 : !! FUNCTION
1876 : !!
1877 : !! INPUTS
1878 : !!
1879 : !! OUTPUT
1880 : !!
1881 : !! NOTES
1882 : !!
1883 : !! SOURCE
1884 :
1885 72 : subroutine symrel2string(symrel1, tnon, string)
1886 :
1887 : !Arguments ------------------------------------
1888 : !scalars
1889 : integer, intent(in) :: symrel1(3,3)
1890 : real(dp), intent(in) :: tnon(3)
1891 : character(len=80), intent(out) :: string
1892 :
1893 : !Local variables -------------------------
1894 : !scalars
1895 : integer :: i1,i2
1896 : character(len=1) :: xyz(3)
1897 :
1898 : ! *********************************************************************
1899 :
1900 72 : xyz(1) = 'x'
1901 72 : xyz(2) = 'y'
1902 72 : xyz(3) = 'z'
1903 :
1904 72 : string = ''
1905 288 : do i1=1,3
1906 216 : if (abs(tnon(i1)) > tol10) then
1907 : ! find fraction 1/n for tnon, otherwise do not know what to print
1908 36 : if (abs(one-two*tnon(i1)) < tol10) string = trim(string)//'1/2'
1909 36 : if (abs(one+two*tnon(i1)) < tol10) string = trim(string)//'-1/2'
1910 :
1911 36 : if (abs(one-three*tnon(i1)) < tol10) string = trim(string)//'1/3'
1912 36 : if (abs(one+three*tnon(i1)) < tol10) string = trim(string)//'-1/3'
1913 36 : if (abs(two-three*tnon(i1)) < tol10) string = trim(string)//'2/3'
1914 36 : if (abs(two+three*tnon(i1)) < tol10) string = trim(string)//'-2/3'
1915 :
1916 36 : if (abs(one-six*tnon(i1)) < tol10) string = trim(string)//'1/6'
1917 36 : if (abs(one+six*tnon(i1)) < tol10) string = trim(string)//'-1/6'
1918 36 : if (abs(five-six*tnon(i1)) < tol10) string = trim(string)//'5/6'
1919 36 : if (abs(five+six*tnon(i1)) < tol10) string = trim(string)//'-5/6'
1920 : end if
1921 864 : do i2=1,3
1922 : ! FIXME: check if this is correct ordering for symrel(i1,i2) looks ok
1923 648 : if (symrel1(i1,i2) == 1) string = trim(string)//'+'//xyz(i2)
1924 864 : if (symrel1(i1,i2) == -1) string = trim(string)//'-'//xyz(i2)
1925 : end do
1926 288 : if (i1 /= 3) string = trim(string)//','
1927 : end do
1928 :
1929 72 : end subroutine symrel2string
1930 : !!***
1931 :
1932 : !!****f* m_crystal/prtposcar
1933 : !! NAME
1934 : !! prtposcar
1935 : !!
1936 : !! FUNCTION
1937 : !! output VASP style POSCAR and FORCES files for use with frozen phonon codes, like
1938 : !! PHON from Dario Alfe' or frophon
1939 : !! IMPORTANT: the order of atoms is fixed such that typat is re-grouped.
1940 : !! First typat=1 then typat=2, etc...
1941 : !! Only master should call this routine in MPI-mode.
1942 : !!
1943 : !! INPUTS
1944 : !! fcart = forces on atoms in cartesian coordinates
1945 : !! natom = number of atoms
1946 : !! ntypat = number of types of atoms
1947 : !! rprimd = lattice vectors for the primitive cell
1948 : !! typat = type for each of the natom atoms
1949 : !! ucvol = unit cell volume
1950 : !! xred = reduced positions of the atoms
1951 : !! znucl = nuclear charge of each atomic type
1952 : !!
1953 : !! OUTPUTS
1954 : !! Only files written
1955 : !!
1956 : !! SOURCE
1957 :
1958 1 : subroutine prtposcar(fcart, fnameradix, natom, ntypat, rprimd, typat, ucvol, xred, znucl)
1959 :
1960 : !Arguments ------------------------------------
1961 : !scalars
1962 : integer, intent(in) :: natom, ntypat
1963 : real(dp), intent(in) :: ucvol
1964 : !arrays
1965 : integer, intent(in) :: typat(natom)
1966 : real(dp), intent(in) :: fcart(3,natom)
1967 : real(dp), intent(in) :: rprimd(3,3)
1968 : real(dp), intent(in) :: xred(3,natom)
1969 : real(dp), intent(in) :: znucl(ntypat)
1970 : character(len=fnlen), intent(in) :: fnameradix
1971 :
1972 : !Local variables-------------------------------
1973 : !scalars
1974 : integer :: iatom, itypat, iout
1975 : type(atomdata_t) :: atom
1976 : ! arrays
1977 2 : integer :: natoms_this_type(ntypat)
1978 : character(len=2) :: symbol
1979 : character(len=7) :: natoms_this_type_str
1980 : character(len=100) :: chem_formula, natoms_all_types, atomsstring_all_types
1981 : character(len=500) :: msg
1982 : !************************************************************************
1983 :
1984 : ! Output POSCAR file for positions, atom types etc
1985 1 : if (open_file(trim(fnameradix)//"_POSCAR", msg, newunit=iout) /= 0) then
1986 0 : ABI_ERROR(msg)
1987 : end if
1988 :
1989 3 : natoms_this_type = 0
1990 3 : do itypat=1,ntypat
1991 21 : do iatom=1,natom
1992 20 : if (typat(iatom) == itypat) natoms_this_type(itypat) = natoms_this_type(itypat) + 1
1993 : end do
1994 : end do
1995 :
1996 1 : chem_formula = ""
1997 1 : atomsstring_all_types = " "
1998 3 : do itypat=1, ntypat
1999 2 : call atomdata_from_znucl(atom, znucl(itypat))
2000 2 : symbol = atom%symbol
2001 2 : if (natoms_this_type(itypat) < 10) then
2002 2 : write(natoms_this_type_str, '(I1)') natoms_this_type(itypat)
2003 0 : else if (natoms_this_type(itypat) < 100) then
2004 0 : write(natoms_this_type_str, '(I2)') natoms_this_type(itypat)
2005 0 : else if (natoms_this_type(itypat) < 1000) then
2006 0 : write(natoms_this_type_str, '(I3)') natoms_this_type(itypat)
2007 : end if
2008 2 : chem_formula = trim(chem_formula) // symbol // trim(natoms_this_type_str)
2009 5 : atomsstring_all_types = trim(atomsstring_all_types) // " " // symbol
2010 : end do
2011 :
2012 1 : write (iout,'(3a,E24.14)') "ABINIT generated POSCAR file. Chemical formula is ",trim(chem_formula), &
2013 2 : " Volume (AA^3) = ", ucvol*Bohr_Ang*Bohr_Ang*Bohr_Ang
2014 :
2015 1 : write (iout,'(E10.1)') 1.0_dp
2016 4 : write (iout,'(3E24.14,1x)') Bohr_Ang*rprimd(:,1) ! (angstr? bohr?)
2017 4 : write (iout,'(3E24.14,1x)') Bohr_Ang*rprimd(:,2)
2018 4 : write (iout,'(3E24.14,1x)') Bohr_Ang*rprimd(:,3)
2019 1 : write (iout, '(a)') atomsstring_all_types
2020 :
2021 1 : natoms_all_types = " "
2022 3 : do itypat=1, ntypat
2023 2 : write (natoms_this_type_str, '(I7)') natoms_this_type(itypat)
2024 3 : natoms_all_types = trim(natoms_all_types) // " " // trim(natoms_this_type_str)
2025 : end do
2026 :
2027 1 : write(iout,'(a)') trim(natoms_all_types)
2028 1 : write(iout,'(a)') "Direct"
2029 :
2030 3 : do itypat=1, ntypat
2031 21 : do iatom=1,natom
2032 18 : if (typat(iatom) /= itypat) cycle
2033 20 : write (iout,'(3(E24.14,1x))') xred(:,iatom)
2034 : end do
2035 : end do
2036 1 : close(iout)
2037 :
2038 : ! output FORCES file for forces in same order as positions above
2039 1 : if (open_file(trim(fnameradix)//"_FORCES", msg, newunit=iout) /= 0 ) then
2040 0 : ABI_ERROR(msg)
2041 : end if
2042 :
2043 : !ndisplacements
2044 : !iatom_displaced displacement_red_coord(3)
2045 : !forces_cart_ev_Angstr(3)
2046 : !...
2047 : !<repeat for other displaced atoms>
2048 1 : write(iout,'(I7)') 1
2049 1 : write(iout,'(a)') '1 0 0 0 ! TO BE FILLED IN '
2050 3 : do itypat=1, ntypat
2051 21 : do iatom=1,natom
2052 18 : if (typat(iatom) /= itypat) cycle
2053 47 : write(iout,'(3(E24.14,1x))') Ha_eV/Bohr_Ang*fcart(:,iatom)
2054 : end do
2055 : end do
2056 :
2057 1 : close(iout)
2058 :
2059 1 : end subroutine prtposcar
2060 : !!***
2061 :
2062 : !----------------------------------------------------------------------
2063 :
2064 : !!****f* m_crystal/crystal_symmetrize_cart_vec3
2065 : !! NAME
2066 : !! crystal_symmetrize_cart_vec3
2067 : !!
2068 : !! FUNCTION
2069 : !! Use spatial and time-reversal symmetry (TR) to symmetrize a Cartesian vector of real elements.
2070 : !!
2071 : !! INPUTS
2072 : !! v(3)=Vector in Cartesian coordinates.
2073 : !! time_opt=Prefactor that defines how the vectors transforms under TR. Usually +1 or -1
2074 : !! Note that TR is used only if cryst%timrev == 2. time_opt = 0 disables TR for testing purposes.
2075 : !!
2076 : !! SOURCE
2077 :
2078 4495 : function crystal_symmetrize_cart_vec3(cryst, v, time_opt) result(vsum)
2079 :
2080 : !Arguments ------------------------------------
2081 : class(crystal_t),intent(in) :: cryst
2082 : real(dp),intent(in) :: v(3)
2083 : integer,intent(in) :: time_opt
2084 : real(dp) :: vsum(3)
2085 :
2086 : !Local variables-------------------------------
2087 : integer :: isym, itime, nsym_sum
2088 : real(dp) :: vsym(3), tsign
2089 : ! *************************************************************************
2090 :
2091 17980 : vsum = zero; nsym_sum = 0
2092 13485 : do itime=1,cryst%timrev
2093 8990 : tsign = 1
2094 8990 : if (itime == cryst%timrev) then
2095 4495 : if (time_opt == 0) cycle
2096 0 : tsign = time_opt
2097 : end if
2098 224750 : do isym=1, cryst%nsym
2099 215760 : nsym_sum = nsym_sum + 1
2100 5609760 : vsym = matmul(cryst%symrel_cart(:,:,isym), v) * tsign
2101 872030 : vsum = vsum + vsym
2102 : end do
2103 : end do
2104 17980 : vsum = vsum / nsym_sum
2105 :
2106 4495 : end function crystal_symmetrize_cart_vec3
2107 : !!***
2108 :
2109 : !----------------------------------------------------------------------
2110 :
2111 : !!****f* m_crystal/crystal_symmetrize_cart_tens33
2112 : !! NAME
2113 : !! crystal_symmetrize_cart_tens33
2114 : !!
2115 : !! FUNCTION
2116 : !! Symmetrize a cartesian 3x3 tensor
2117 : !! Use spatial and time-reversal symmetry (TR) to symmetrize a Cartesian 3x3 tensor of real elements.
2118 : !!
2119 : !! INPUTS
2120 : !! v(3)=Vector in Cartesian coordinates.
2121 : !! time_opt=Prefacator that defines how the vectors transforms under TR. Usually +1 or -1
2122 : !! Note that TR is used only if cryst%timrev == 2. time_opt = 0 disables TR for testing purposes.
2123 : !!
2124 : !! SOURCE
2125 :
2126 62537 : function crystal_symmetrize_cart_tens33(cryst, t, time_opt) result(tsum)
2127 :
2128 : !Arguments ------------------------------------
2129 : class(crystal_t),intent(in) :: cryst
2130 : real(dp),intent(in) :: t(3,3)
2131 : integer,intent(in) :: time_opt
2132 : real(dp) :: tsum(3,3)
2133 :
2134 : !Local variables-------------------------------
2135 : integer :: isym, itime, nsym_sum
2136 : real(dp) :: tsym(3,3), tsign
2137 : ! *************************************************************************
2138 :
2139 812981 : tsum = zero; nsym_sum = 0
2140 :
2141 187611 : do itime=1,cryst%timrev
2142 125074 : tsign = 1
2143 125074 : if (itime == cryst%timrev) then
2144 62537 : if (time_opt == 0) cycle
2145 0 : tsign = time_opt
2146 : end if
2147 1991722 : do isym=1, cryst%nsym
2148 1866648 : nsym_sum = nsym_sum + 1
2149 324796752 : tsym = tsign * matmul((cryst%symrel_cart(:,:,isym)), matmul(t, transpose(cryst%symrel_cart(:,:,isym))))
2150 24391498 : tsum = tsum + tsym
2151 : end do
2152 : end do
2153 :
2154 812981 : tsum = tsum / nsym_sum
2155 :
2156 62537 : end function crystal_symmetrize_cart_tens33
2157 : !!***
2158 :
2159 : !----------------------------------------------------------------------
2160 :
2161 : !!****f* m_crystal/get_recart_qdirs
2162 : !! NAME
2163 : !! get_recart_qdirs
2164 : !!
2165 : !! FUNCTION
2166 : !! Return predefined list of 6 q-versors in reciprocal space reduced coordinates.
2167 : !! First 3 entries are along the recip. space lattice vectors, then along the Cartesian axis x,y,z.
2168 : !! The optional qlen argument, can be used to rescale the vectors. Default: 1
2169 : !!
2170 : !! INPUTS
2171 : !!
2172 : !! SOURCE
2173 :
2174 29 : subroutine get_redcart_qdirs(cryst, nq, qdirs, qlen)
2175 :
2176 : !Arguments ------------------------------------
2177 : class(crystal_t),intent(in) :: cryst
2178 : integer,intent(out) :: nq
2179 : real(dp),allocatable,intent(out) :: qdirs(:,:)
2180 : real(dp),optional,intent(in) :: qlen
2181 :
2182 : !Local variables-------------------------------
2183 : integer :: iq
2184 : real(dp) :: qred2cart(3,3), qcart2red(3,3)
2185 : ! *************************************************************************
2186 :
2187 377 : qred2cart = two_pi * cryst%gprimd
2188 29 : qcart2red = qred2cart
2189 29 : call matrginv(qcart2red, 3, 3)
2190 :
2191 29 : nq = 6
2192 29 : ABI_MALLOC(qdirs, (3, nq))
2193 116 : qdirs(:,1) = [one, zero, zero] ! (100)
2194 116 : qdirs(:,2) = [zero, one, zero] ! (010)
2195 116 : qdirs(:,3) = [zero, zero, one] ! (001)
2196 493 : qdirs(:,4) = matmul(qcart2red, [one, zero, zero]) ! (x)
2197 493 : qdirs(:,5) = matmul(qcart2red, [zero, one, zero]) ! (y)
2198 493 : qdirs(:,6) = matmul(qcart2red, [zero, zero, one]) ! (z)
2199 :
2200 : ! normalization
2201 203 : do iq=1,nq
2202 725 : qdirs(:,iq) = qdirs(:,iq) / normv(qdirs(:,iq), cryst%gmet, "G")
2203 : end do
2204 :
2205 29 : if (present(qlen)) qdirs = qlen * qdirs
2206 :
2207 29 : end subroutine get_redcart_qdirs
2208 : !!***
2209 :
2210 : !----------------------------------------------------------------------
2211 :
2212 : !!****f* m_crystal/crystal_write_xsf_data
2213 : !! NAME
2214 : !! crystal_write_xsf_data
2215 : !!
2216 : !! FUNCTION
2217 : !! Write array in real space in XSF format
2218 : !!
2219 : !! INPUTS
2220 : !! filepath: Path to output file.
2221 : !! n1, n2, n3: Dimensions of input array.
2222 : !! datagrid(n1,n2,n3)=Array in real space.
2223 : !!
2224 : !! SOURCE
2225 :
2226 0 : subroutine crystal_write_xsf_data(cryst, filepath, n1, n2, n3, datagrid)
2227 :
2228 : !Arguments ------------------------------------
2229 : class(crystal_t),intent(in) :: cryst
2230 : character(len=*),intent(in) :: filepath
2231 : integer,intent(in) :: n1, n2, n3
2232 : real(dp),intent(in) :: datagrid(n1*n2*n3)
2233 :
2234 : !Local variables-------------------------------
2235 : integer,parameter :: realrecip = 0
2236 : real(dp),parameter :: origin(3) = zero
2237 : ! *************************************************************************
2238 :
2239 : call write_xsf(filepath, n1, n2, n3, datagrid, cryst%rprimd, origin, &
2240 0 : cryst%natom, cryst%ntypat, cryst%typat, cryst%xcart, cryst%znucl, realrecip)
2241 :
2242 0 : end subroutine crystal_write_xsf_data
2243 : !!***
2244 :
2245 1868235 : end module m_crystal
2246 : !!***
|
