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