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