Branch data Line data Source code
1 : : !!****m* ABINIT/m_symfind
2 : : !! NAME
3 : : !! m_symfind
4 : : !!
5 : : !! FUNCTION
6 : : !! Symmetry finder high-level API.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2000-2026 ABINIT group (XG, RC)
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_symfind
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_abicore
27 : : use m_symlist
28 : :
29 : :
30 : : use m_matrix, only : mati3inv, matr3inv, mati3det
31 : : use m_symtk, only : chkprimit, symrelrot, symdet, symcharac, holocell, symatm, &
32 : : smallprim, print_symmetries, sg_multable, symmetrize_tnons, symmetrize_xred
33 : : use m_geometry, only : acrossb, xred2xcart
34 : : use m_spgdata, only : getptgroupma, symptgroup, spgdata
35 : :
36 : : implicit none
37 : :
38 : : private
39 : : !!***
40 : :
41 : : public :: symfind ! From the symmetries of the Bravais lattice,
42 : : ! select those that leave invariant the system, and generate tnons. Not always robust.
43 : : public :: symfind_expert ! Wrap symfind to provide robust determination of the symmetries,
44 : : ! for which resymmetrization of atomic positions and tnons is needed.
45 : : public :: symanal ! Find the space group from the list of symmetries and lattice parameters
46 : : public :: symbrav ! Determine the Bravais information from the list of symmetry operations, and the lattice vectors.
47 : : public :: symlatt ! Find the Bravais lattice and its symmetry operations (ptsymrel).
48 : : ! From the unit cell vectors (rprimd) and the corresponding metric tensor.
49 : :
50 : : contains
51 : : !!***
52 : :
53 : : !!****f* m_symfind/symfind
54 : : !! NAME
55 : : !! symfind
56 : : !!
57 : : !! FUNCTION
58 : : !! Symmetry finder.
59 : : !! From the symmetries of the Bravais lattice (ptsymrel),
60 : : !! select those that leave invariant the system, and generate
61 : : !! the corresponding tnons vectors and symafm information.
62 : : !! Unlike symfind_expert, does NOT resymmetrize atomic positions and tnons for more robust determination of the symmetries.
63 : : !! The algorithm is explained in T.G. Worlton and J.L. Warren, Comp. Phys. Comm. 3, 88 (1972) [[cite:Worlton1972]]
64 : : !!
65 : : !! INPUTS
66 : : !! chrgat(natom) (optional)=target charge for each atom. Not always used, it depends on the value of constraint_kind
67 : : !! invardir_red (optional)=reduced coordinates of an invariant direction (only acting with symrel - not tnons)
68 : : !! invaraxial_red (optional)=reduced coordinates of an axial vector, similar to invardir_red, but transforms with an additional
69 : : !! deteminant factor under symrel operations
70 : : !! invar_z (optional)= if 1, the z direction must stay invariant for all symrel applied ;
71 : : !! if 2, z must stay invariant and also there cannot be any associated tnons along z.
72 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space
73 : : !! msym=default maximal number of symmetries
74 : : !! natom=number of atoms in cell.
75 : : !! nptsym=number of point symmetries of the Bravais lattice
76 : : !! nspden= number of spin-density components. When 4, the three components of spinat are taken into account, instead of only z-component.
77 : : !! nucdipmom(3,natom) (optional) array of nuclear dipole moments
78 : : !! ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
79 : : !! of the Bravais lattice in real space in terms
80 : : !! of primitive translations.
81 : : !! spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
82 : : !! tolsym=tolerance for the symmetries
83 : : !! typat(natom)=integer identifying type of atom.
84 : : !! use_inversion=1 if inversion and improper rotations can be included in set of symmetries
85 : : !! xred(3,natom)=reduced coordinates of atoms in terms of real space
86 : : !! primitive translations
87 : : !!
88 : : !! OUTPUT
89 : : !! ierr (optional)=if non-zero, the symmetry operations do not form a group
90 : : !! nsym=actual number of symmetries
91 : : !! symafm(1:msym)=(anti)ferromagnetic part of nsym symmetry operations
92 : : !! symrel(3,3,1:msym)= nsym symmetry operations in real space in terms
93 : : !! of primitive translations
94 : : !! tnons(3,1:msym)=nonsymmorphic translations for each symmetry (would
95 : : !! be 0 0 0 each for a symmorphic space group)
96 : : !!
97 : : !! SOURCE
98 : :
99 : 18163 : subroutine symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
100 : 18163 : prtvol, ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
101 : 35898 : chrgat,ierr,nucdipmom,invardir_red,invaraxial_red,invar_z) ! Optional
102 : :
103 : : !Arguments ------------------------------------
104 : : !scalars
105 : : integer,intent(in) :: msym,natom,nptsym,nspden,use_inversion
106 : : integer,intent(in) :: prtvol
107 : : integer,optional,intent(in) :: invar_z
108 : : integer,optional,intent(out) :: ierr
109 : : integer,intent(out) :: nsym
110 : : real(dp),intent(in) :: tolsym
111 : : !arrays
112 : : integer,intent(in) :: ptsymrel(3,3,msym),typat(natom)
113 : : integer,intent(inout) :: symafm(msym),symrel(3,3,msym) !vz_i
114 : : real(dp),intent(in) :: gprimd(3,3),spinat(3,natom),xred(3,natom)
115 : : real(dp),optional,intent(in) :: invardir_red(3),invaraxial_red(3),chrgat(natom)
116 : : real(dp),optional, intent(in) :: nucdipmom(3,natom)
117 : : real(dp),intent(inout) :: tnons(3,msym) !vz_i
118 : :
119 : : !Local variables-------------------------------
120 : : !scalars
121 : : integer :: found3,foundcl,iatom,iatom0,iatom1,iatom2,iatom3,iclass,iclass0,ierr_,ii
122 : : integer :: isym,jj,kk,natom0,nclass,ntrial,printed,trialafm,trialok
123 : : integer :: mm(3,3), detR
124 : : real(dp) :: det,diff1,diff2,diff3,diffr1,diffr2,diffr3,ndnorm,nucdipmomcl2,nucdipmomcl20
125 : : real(dp) :: spinat2,spinatcl2,spinatcl20,tolsym2
126 : : ! TRUE if antiferro symmetries are used with non-collinear magnetism.
127 : : integer :: afm_noncoll=1 !For nspden==4. If 1, all symops are permitted ; if 0 symafm must be 1.
128 : : !For nspden=4. If noncoll_orthorhombic1, require the symmetry operations to be a subset of the orthorhombic symmetries, except if all spinat=0..
129 : : integer :: noncoll_orthorhombic=0
130 : : logical :: test_sameabsspin,test_samechrg
131 : : logical :: test_samenucdipmom
132 : : character(len=500) :: msg
133 : : !arrays
134 : 18163 : integer,allocatable :: class(:,:),natomcl(:),typecl(:)
135 : : real(dp) :: diff(3),invardir_red_rot(3),invaraxial_red_rot(3),hand2(3),hand3(3),ndtest(3),rprimd(3,3),spinat0(3),xred0(3)
136 : : !real(dp) :: symnucdipmom2(3)
137 : : real(dp) :: symnucdipmom2cart(3,3),symnucdipmom2red(3,3)
138 : : real(dp) :: symspinat1(3),symspinat2(3),symxred2(3),trialnons(3)
139 : 18163 : real(dp),allocatable :: chrgat_(:)
140 : 18163 : real(dp),allocatable :: chrgatcl(:)
141 : 18163 : real(dp),allocatable :: local_nucdipmom(:,:,:),nucdipmomcl(:,:),nucdipmomred(:,:,:)
142 : 18163 : real(dp),allocatable :: spinatcl(:,:),spinatred(:,:)
143 : : !**************************************************************************
144 : :
145 : : !write(std_out,'(a)')' m_symfind%symfind : enter '; call flush(std_out)
146 : :
147 [ - + + - : 54489 : ABI_MALLOC(local_nucdipmom,(3,3,natom))
+ - + - -
+ - + ]
148 [ + + + + : 601031 : local_nucdipmom(:,:,:) = zero
+ + ]
149 [ + + ]: 18163 : if(present(nucdipmom)) then
150 [ + + + + ]: 194669 : local_nucdipmom(1:3,1,:) = nucdipmom(1:3,:)
151 : : end if
152 : : !DEBUG
153 : : ! if (prtvol>1) msg="remove me later"
154 : : ! write(std_out,*)' symfind : enter'
155 : : ! call flush(6)
156 : : ! write(std_out,*)' ptsymrel matrices are :'
157 : : ! do isym=1,nptsym
158 : : ! write(std_out,'(i4,4x,9i4)' )isym,ptsymrel(:,:,isym)
159 : : ! end do
160 : : ! write(std_out,*)' symfind : natom=',natom
161 : : ! do iatom=1,natom
162 : : ! write(std_out,*)' atom number',iatom
163 : : ! write(std_out,*)' typat =',typat(iatom)
164 : : ! write(std_out,*)' spinat =',spinat(:,iatom)
165 : : ! write(std_out,*)' xred =',xred(:,iatom)
166 : : ! if (ANY(ABS(local_nucdipmom)>tol8)) then
167 : : ! write(std_out,*)' nucdipmom = ',local_nucdipmom(:,1,iatom)
168 : : ! end if
169 : : ! end do
170 : : ! write(std_out,*)' '
171 : : ! call flush(6)
172 : : !ENDDEBUG
173 : :
174 : : ABI_UNUSED(prtvol)
175 : :
176 : : !Find the number of classes of atoms (type, chrg and spinat must be identical,
177 : : !spinat might differ by a sign, if aligned with the z direction, or,
178 : : ! type and nucdipmom must be identical)
179 : : !natomcl(iclass) will contain the number of atoms in the class
180 : : !typecl(iclass) will contain the type of the atoms in the class
181 : : !chrgcl(iclass) will contain the charge of the atoms in the class
182 : : !spinatcl(1:3,iclass) will contain the spinat of the atoms in the class
183 : : !class(1:natomclass(iclass),iclass) will contain the index of the
184 : : !atoms belonging to the class
185 [ + - - + : 72652 : ABI_MALLOC(class,(natom+3,natom))
+ - - + +
- + - + -
- + - + ]
186 [ - + + - : 54489 : ABI_MALLOC(natomcl,(natom))
+ - + - -
+ - + ]
187 [ - + + - : 36326 : ABI_MALLOC(typecl,(natom))
- + - + ]
188 [ - + + - : 54489 : ABI_MALLOC(chrgat_,(natom))
+ - + - -
+ - + ]
189 [ - + + - : 36326 : ABI_MALLOC(chrgatcl,(natom))
- + - + ]
190 [ - + + - : 54489 : ABI_MALLOC(spinatcl,(3,natom))
+ - + - -
+ - + ]
191 [ - + + - : 36326 : ABI_MALLOC(nucdipmomcl,(3,natom))
- + - + ]
192 : :
193 : 18163 : tolsym2=tolsym**2
194 : :
195 [ + + ]: 62999 : chrgat_(:)=zero
196 [ + + ]: 18163 : if(present(chrgat))then
197 [ + + ]: 62129 : chrgat_(:)=chrgat(:)
198 : : endif
199 : :
200 : : ! for each nuclear dipole we need a local right handed coord system, so we can
201 : : ! test later for whether a symmetry operation preserves the circulation induced
202 : : ! by the dipole
203 [ + + ]: 62999 : do iatom=1, natom
204 [ + + ]: 179344 : ndnorm=sqrt(DOT_PRODUCT(local_nucdipmom(1:3,1,iatom),local_nucdipmom(1:3,1,iatom)))
205 : :
206 : : ! if nuclear dipole has effectively zero size, move on to the next atom
207 [ + + ]: 44836 : if (ndnorm < tol8) cycle
208 : :
209 : : ! for testing purposes, we care only about direction so renormalize to unity
210 [ + + ]: 512 : local_nucdipmom(1:3,1,iatom) = local_nucdipmom(1:3,1,iatom)/ndnorm
211 : :
212 : : ! make a random vector, each component is (0,1]
213 : 128 : call random_number(ndtest)
214 : :
215 : : ! vector 2 is constructed to be orthogonal to original nuclear dipole moment vector
216 : 128 : call acrossb(local_nucdipmom(1:3,1,iatom),ndtest(1:3),local_nucdipmom(1:3,2,iatom))
217 : :
218 : : ! vector 3 is orthogonal to 1 and 2, and 1,2,3 form a right-handed set
219 : 62999 : call acrossb(local_nucdipmom(1:3,1,iatom),local_nucdipmom(1:3,2,iatom),local_nucdipmom(1:3,3,iatom))
220 : : end do
221 : :
222 : : ! need rprimd later to transform back to cart coords
223 : 18163 : call matr3inv(gprimd,rprimd)
224 : :
225 : : !write(std_out,'(a)')' m_symfind%symfind : before initialise with the first atom '; call flush(std_out)
226 : :
227 : : !Initialise with the first atom
228 : 18163 : nclass=1
229 : 18163 : natomcl(1)=1
230 : 18163 : typecl(1)=typat(1)
231 : 18163 : chrgatcl(1)=chrgat_(1)
232 [ + + ]: 72652 : spinatcl(:,1)=spinat(:,1)
233 [ + + ]: 72652 : nucdipmomcl(:,1)=local_nucdipmom(:,1,1)
234 : 18163 : class(1,1)=1
235 [ + + ]: 18163 : if(natom>1)then
236 [ + + ]: 40435 : do iatom=2,natom
237 : : ! DEBUG
238 : : ! write(std_out,*)' '
239 : : ! write(std_out,*)' symfind : examine iatom=',iatom
240 : : ! ENDDEBUG
241 : 39033 : foundcl=0
242 [ + + ]: 39033 : do iclass=1,nclass
243 : : ! Compare the typat, chrg and spinat of atom iatom with existing ones.
244 : : ! At this stage, admit either identical spinat, or spin-flip spinat.
245 [ + + ]: 39033 : if( typat(iatom)==typecl(iclass)) then
246 : 20654 : test_samechrg= (abs(chrgat_(iatom)-chrgatcl(iclass))<tolsym)
247 [ + + ]: 20654 : if(nspden/=4)then
248 : 20010 : test_sameabsspin=(abs(abs(spinat(3,iatom))-abs(spinatcl(3,iclass)))<tolsym)
249 : : else if(nspden==4)then
250 : 644 : spinat2 =spinat(1,iatom)**2+spinat(2,iatom)**2+spinat(3,iatom)**2
251 : 644 : spinatcl2=spinatcl(1,iclass)**2+spinatcl(2,iclass)**2+spinatcl(3,iclass)**2
252 : 644 : test_sameabsspin=abs(spinat2-spinatcl2)<tolsym
253 : : endif
254 : : test_samenucdipmom= &
255 : : & abs(local_nucdipmom(1,1,iatom)-nucdipmomcl(1,iclass))<tolsym .and. &
256 : : & abs(local_nucdipmom(2,1,iatom)-nucdipmomcl(2,iclass))<tolsym .and. &
257 [ + - + - : 20654 : & abs(local_nucdipmom(3,1,iatom)-nucdipmomcl(3,iclass))<tolsym
- + ]
258 : : ! note in the following test, m_chkinp/chkinp has already prevented nucdipmom to be
259 : : ! nonzero when spinat is nonzero
260 [ + + + - ]: 20654 : if( test_samechrg .and. test_sameabsspin .and. test_samenucdipmom ) then
261 : : ! write(std_out,*)' symfind : find it belongs to class iclass=',iclass
262 : : ! write(std_out,*)' symfind : spinat(:,iatom)=',spinat(:,iatom)
263 : : ! write(std_out,*)' symfind : spinatcl(:,iclass)=',spinatcl(:,iclass)
264 : : ! write(std_out,*)' symfind : test_sameabsspin=',test_sameabsspin
265 : : ! write(std_out,*)' '
266 : 20572 : natomcl(iclass)=natomcl(iclass)+1
267 : 20572 : class(natomcl(iclass),iclass)=iatom
268 : : foundcl=1
269 : : exit
270 : : end if
271 : : end if
272 : : end do
273 : : ! If no class with these characteristics exist, create one
274 : 13762 : if(foundcl==0)then
275 : 6101 : nclass=nclass+1
276 : 6101 : natomcl(nclass)=1
277 : 6101 : typecl(nclass)=typat(iatom)
278 : 6101 : chrgatcl(nclass)=chrgat_(iatom)
279 [ + + ]: 24404 : spinatcl(:,nclass)=spinat(:,iatom)
280 [ + + ]: 24404 : nucdipmomcl(:,nclass)=local_nucdipmom(:,1,iatom)
281 : 6101 : class(1,nclass)=iatom
282 : : end if
283 : : end do
284 : : end if
285 : :
286 : : !write(std_out,*)' '
287 : : !write(std_out,*)' symfind : found ',nclass,' nclass of atoms'
288 : : !do iclass=1,nclass
289 : : !write(std_out,*)' class number',iclass
290 : : !write(std_out,*)' natomcl =',natomcl(iclass)
291 : : !write(std_out,*)' typecl =',typecl(iclass)
292 : : !write(std_out,*)' spinatcl=',spinatcl(:,iclass)
293 : : !write(std_out,*)' class =',(class(iatom,iclass),iatom=1,natomcl(iclass))
294 : : !end do
295 : : !write(std_out,*)' '
296 : :
297 : : !write(std_out,'(a)')' m_symfind%symfind : before select the class '; call flush(std_out)
298 : :
299 : : !Select the class with the least number of atoms, and non-zero spinat if any
300 : : !It is important to select a magnetic class of atom, if any, otherwise
301 : : !the determination of the initial (inclusive) set of symmetries takes only
302 : : !non-magnetic symmetries, and not both magnetic and non-magnetic ones, see later.
303 : : !On the contrary, the chrgat_ data does not play any role, it is invariant upon atomic-centered symmetries
304 : 18163 : iclass0=1
305 : 18163 : natom0=natomcl(1)
306 : 18163 : spinatcl20=spinatcl(1,1)**2+spinatcl(2,1)**2+spinatcl(3,1)**2
307 : 18163 : nucdipmomcl20=nucdipmomcl(1,1)**2+nucdipmomcl(2,1)**2+nucdipmomcl(3,1)**2
308 [ + + ]: 18163 : if(nclass>1)then
309 [ + + ]: 11834 : do iclass=2,nclass
310 : 6101 : spinatcl2=spinatcl(1,iclass)**2+spinatcl(2,iclass)**2+spinatcl(3,iclass)**2
311 : 6101 : nucdipmomcl2=nucdipmomcl(1,iclass)**2+nucdipmomcl(2,iclass)**2+nucdipmomcl(3,iclass)**2
312 : : if( (natomcl(iclass)<natom0 &
313 : : & .and. .not. (spinatcl20>tolsym .and. spinatcl2<tolsym) &
314 : : & .and. .not. (nucdipmomcl20>tolsym .and. nucdipmomcl2<tolsym) ) &
315 : : & .or. (spinatcl20<tolsym .and. spinatcl2>tolsym) &
316 [ + + + + : 11834 : & .or. (nucdipmomcl20<tolsym .and. nucdipmomcl2>tolsym)) then
- + - + -
- + + + +
+ + + + ]
317 : 6101 : iclass0=iclass
318 : 6101 : natom0=natomcl(iclass)
319 : 6101 : spinatcl20=spinatcl2
320 : 6101 : nucdipmomcl20=nucdipmomcl2
321 : : end if
322 : : end do
323 : : end if
324 : :
325 : 18163 : printed=0
326 : :
327 : : !If non-collinear spinat have to be used, transfer them in reduced coordinates
328 [ + + ]: 18163 : if (nspden==4) then
329 [ - + + - : 1092 : ABI_MALLOC(spinatred,(3,natom))
- + - + ]
330 [ + + ]: 1890 : do iatom=1,natom
331 [ + + ]: 5922 : do ii=1,3
332 : 5376 : spinatred(1:3,iatom)=MATMUL(TRANSPOSE(gprimd),spinat(1:3,iatom))
333 : : end do
334 : : end do
335 : : end if
336 : :
337 : : !write(std_out,*)' '
338 : : !write(std_out,*)' symfind : has selected iclass0=',iclass0
339 : : !write(std_out,*)' # iatom xred spinat (spinatred if nspden=4) '
340 : : !do iatom0=1,natomcl(iclass0)
341 : : !iatom=class(iatom0,iclass0)
342 : : !if(nspden/=4)then
343 : : ! write(std_out,'(2i4,6f10.4)' )iatom0,iatom,xred(:,iatom),spinat(:,iatom)
344 : : !else if(nspden==4)then
345 : : ! write(std_out,'(2i4,9f10.4)' )iatom0,iatom,xred(:,iatom),spinat(:,iatom),spinatred(:,iatom)
346 : : !endif
347 : : !end do
348 : : !write(std_out,*)' '
349 : :
350 : : !represent nuclear dipole moments in reduced coords
351 [ - + + - : 36326 : ABI_MALLOC(nucdipmomred,(3,3,natom))
- + - + ]
352 [ + + ]: 62999 : do iatom=1,natom
353 [ + + ]: 197507 : do ii=1,3
354 : 179344 : nucdipmomred(1:3,ii,iatom)=MATMUL(TRANSPOSE(gprimd),local_nucdipmom(1:3,ii,iatom))
355 : : end do
356 : : end do
357 : :
358 : : !write(std_out,'(a)')' m_symfind%symfind : before big loop '; call flush(std_out)
359 : :
360 : : !Big loop over each symmetry operation of the Bravais lattice
361 : 18163 : nsym=0
362 [ + + ]: 761883 : do isym=1,nptsym
363 : :
364 : : !write(std_out,'(a,i4)')' m_symfind%symfind : enter loop isym=',isym; call flush(std_out)
365 : :
366 [ + + ]: 743720 : if(present(invardir_red))then
367 : : ! ji: Check whether symmetry operation leaves invardir_red invariant
368 : : invardir_red_rot(:) = ptsymrel(:,1,isym)*invardir_red(1) + &
369 : : & ptsymrel(:,2,isym)*invardir_red(2) + &
370 [ + + ]: 2939488 : & ptsymrel(:,3,isym)*invardir_red(3)
371 [ + + ]: 2939488 : diff(:)=invardir_red(:)-invardir_red_rot(:)
372 [ + + ]: 734872 : if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 ) cycle
373 : : endif
374 : :
375 : : ! check whether symmetry operation leaves invaraxial_red invariant (axial vector)
376 [ + + ]: 740760 : if(present(invaraxial_red))then
377 : : !! CALL MatrixHyb_getDet(symref_hyb(isym), detR)
378 [ + + + + ]: 9514856 : mm(:,:) = ptsymrel(:,:,isym)
379 : 731912 : call mati3det(mm,detR)
380 : : invaraxial_red_rot(:) = detR * (&
381 : : ptsymrel(:,1,isym)*invaraxial_red(1) + &
382 : : & ptsymrel(:,2,isym)*invaraxial_red(2) + &
383 [ + + ]: 2927648 : & ptsymrel(:,3,isym)*invaraxial_red(3) )
384 [ + + ]: 2927648 : diff(:)=invaraxial_red(:)-invaraxial_red_rot(:)
385 : :
386 [ + + ]: 731912 : if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 ) cycle
387 : : endif
388 : :
389 : : !write(std_out,'(a,i4)')' m_symfind%symfind : 1'; call flush(std_out)
390 : :
391 : : det=ptsymrel(1,1,isym)*ptsymrel(2,2,isym)*ptsymrel(3,3,isym)+&
392 : : & ptsymrel(2,1,isym)*ptsymrel(3,2,isym)*ptsymrel(1,3,isym)+&
393 : : & ptsymrel(1,2,isym)*ptsymrel(2,3,isym)*ptsymrel(3,1,isym) - &
394 : : & (ptsymrel(3,1,isym)*ptsymrel(2,2,isym)*ptsymrel(1,3,isym)+&
395 : : & ptsymrel(2,1,isym)*ptsymrel(1,2,isym)*ptsymrel(3,3,isym)+&
396 : 739976 : & ptsymrel(3,2,isym)*ptsymrel(2,3,isym)*ptsymrel(1,1,isym))
397 [ + + + + ]: 739976 : if(use_inversion==0 .and. det==-1) cycle
398 : :
399 : : !write(std_out,'(a,i4)')' m_symfind%symfind : 2'; call flush(std_out)
400 : :
401 : : ! jellium slab and spatially varying chemical potential cases:
402 : : ! (actually, an inversion symmetry/mirror plane perpendicular to z symmetry operation might still be allowed... TO BE DONE !)
403 [ + + ]: 729612 : if(present(invar_z))then
404 [ + + ]: 720772 : if (invar_z/=0) then
405 : : ! check whether symmetry operation produce a rotation only in the xy plane
406 : : if (ptsymrel(1,3,isym)/=0 .or. ptsymrel(2,3,isym)/=0 .or. &
407 [ + + + + : 1856 : ptsymrel(3,1,isym)/=0 .or. ptsymrel(3,2,isym)/=0 ) cycle
- + - + ]
408 : : ! check whether symmetry operation does not change the z
409 [ + + ]: 832 : if( ptsymrel(3,3,isym)/=1 ) cycle
410 : : end if
411 : : end if
412 : :
413 : : !write(std_out,'(a,i4)')' m_symfind%symfind : 3'; call flush(std_out)
414 : :
415 : : ! If noncoll_orthorhombic=1, require orthorhombic operations of symmetries, except if spinat=0.
416 [ + + - + ]: 728172 : if (nspden==4 .and. noncoll_orthorhombic==1)then
417 [ # # # # : 0 : if(sum(abs(spinat(:,:)))>tol14)then
# # ]
418 : : if( ptsymrel(1,3,isym)/=0 .or. ptsymrel(2,3,isym)/=0 .or. &
419 : : & ptsymrel(1,2,isym)/=0 .or. ptsymrel(3,2,isym)/=0 .or. &
420 [ # # # # : 0 : & ptsymrel(2,1,isym)/=0 .or. ptsymrel(3,2,isym)/=0 ) cycle
# # # # #
# ]
421 : : endif
422 : : endif
423 : :
424 : : !write(std_out,'(a,i4)')' m_symfind%symfind : 4'; call flush(std_out)
425 : :
426 : : ! Select a tentative set of associated translations
427 : : ! First compute the symmetric of the first atom in the smallest class,
428 : : ! using the point symmetry, and also the symmetric of spinat(red).
429 : 728172 : iatom0=class(1,iclass0)
430 : : xred0(:)=ptsymrel(:,1,isym)*xred(1,iatom0)+ &
431 : : ptsymrel(:,2,isym)*xred(2,iatom0)+ &
432 [ + + ]: 2912688 : ptsymrel(:,3,isym)*xred(3,iatom0)
433 [ + + ]: 728172 : if (nspden/=4) then
434 [ + + ]: 2857744 : spinat0(:)=spinat(:,iatom0)
435 : : else
436 : : spinat0(:)=det*(&
437 : : ptsymrel(:,1,isym)*spinatred(1,iatom0)+ &
438 : : & ptsymrel(:,2,isym)*spinatred(2,iatom0)+ &
439 [ + + ]: 54944 : & ptsymrel(:,3,isym)*spinatred(3,iatom0))
440 : : ! spinat should be treated as an axial vector
441 : : ! i.e. the improper part of a symm. op. has no effect on spinat
442 : : endif
443 : :
444 : : !write(std_out,'(a,i4)')' m_symfind%symfind : 5'; call flush(std_out)
445 : :
446 : : ! From the set of possible images, deduce tentative translations,
447 : : ! and magnetic factor then test whether it send each atom on a symmetric one
448 : : ntrial=0
449 [ + + ]: 2027991 : do ii=1,natom0
450 : : !write(std_out,'(a,2i4)')' symfind : loop isym,ii=',isym,ii
451 : 1281656 : iatom1=class(ii,iclass0)
452 : :
453 : : ! The tentative translation is found
454 [ + + ]: 5126624 : trialnons(:)=xred(:,iatom1)-xred0(:)
455 : : ! Compare the spinat vectors
456 [ + + ]: 1281656 : if (nspden/=4) then
457 [ + + ]: 5045488 : symspinat1(:)=spinat(:,iatom1)
458 : : else
459 [ + + ]: 81136 : symspinat1(:)=spinatred(:,iatom1)
460 : : end if
461 : :
462 : : !write(std_out,'(a,6f10.4)')' symspinat1,spinat0=',symspinat1(:),spinat0(:)
463 : :
464 : 1281656 : trialafm=1
465 [ + + + + ]: 5126624 : if(sum(abs(symspinat1(:)-spinat0(:)))>tolsym)then
466 : 16966 : trialafm=-1
467 [ + + + - ]: 16966 : if(nspden==4 .and. afm_noncoll==0)cycle
468 [ + + + + ]: 67864 : if(sum(abs(symspinat1(:)+spinat0(:)))>tolsym)cycle
469 : : endif
470 : :
471 [ + + - + ]: 5087232 : if(sum(abs(local_nucdipmom(:,1,iatom1)-local_nucdipmom(:,1,iatom0)))>tolsym)then
472 : : write(msg,'(3a,3i5)')&
473 : 0 : 'Problem with matching the nuclear dipole moment within a class.',ch10,&
474 : 0 : 'isym,iatom0,iatom1=',isym,iatom0,iatom1
475 : 0 : ABI_ERROR(msg)
476 : : end if
477 : :
478 : : ! jellium slab case: check whether symmetry operation has no translational
479 : : ! component along z
480 [ + + ]: 1271808 : if(present(invar_z))then
481 [ + + + - ]: 1258596 : if( invar_z==2 .and. abs(trialnons(3)) > tolsym ) cycle
482 : : endif
483 : : trialok=1
484 : :
485 : : ! write(std_out, '(a,i3,a,i3,a,i3,a,3f12.4,i3)') ' Try isym=',isym,' sending iatom0 ',iatom0,' to iatom1 ',iatom1,' with trialnons(:),trialafm =',trialnons(:),trialafm
486 : :
487 : : ! Loop over all classes, then all atoms in the class,
488 : : ! to find whether they have a symmetric
489 [ + + ]: 2086166 : do iclass=1,nclass
490 [ + + ]: 4186758 : do jj=1,natomcl(iclass)
491 : :
492 : 3372400 : iatom2=class(jj,iclass)
493 : : ! Generate the tentative symmetric position of iatom2
494 : : symxred2(:)=ptsymrel(:,1,isym)*xred(1,iatom2)+ &
495 : : & ptsymrel(:,2,isym)*xred(2,iatom2)+ &
496 [ + + ]: 13489600 : & ptsymrel(:,3,isym)*xred(3,iatom2)+ trialnons(:)
497 : : ! Generate the tentative symmetric spinat of iatom2
498 [ + + ]: 3372400 : if (nspden/=4) then
499 [ + + ]: 13391744 : symspinat2(:)=trialafm*spinat(:,iatom2)
500 : : else
501 : : symspinat2(:)=trialafm*det*(ptsymrel(:,1,isym)*spinatred(1,iatom2)+ &
502 : : & ptsymrel(:,2,isym)*spinatred(2,iatom2)+ &
503 [ + + ]: 97856 : & ptsymrel(:,3,isym)*spinatred(3,iatom2))
504 : : end if
505 : : ! Generate the tentative symmetric nucdipmom of iatom2
506 [ + + ]: 13489600 : do kk = 1, 3
507 : : symnucdipmom2red(:,kk)=ptsymrel(:,1,isym)*nucdipmomred(1,kk,iatom2)+ &
508 : : & ptsymrel(:,2,isym)*nucdipmomred(2,kk,iatom2)+ &
509 [ + + ]: 40468800 : & ptsymrel(:,3,isym)*nucdipmomred(3,kk,iatom2)
510 : : ! transform back to cart coords for final comparison to nucdipmom
511 [ + + + + : 165247600 : symnucdipmom2cart(:,kk)=MATMUL(rprimd,symnucdipmom2red(:,kk))
+ + ]
512 : : end do
513 : :
514 : : ! write(std_out,'(a,i4,a,3f8.4,a,3f8.4,a,3f8.4)')&
515 : : ! ' Test iatom2=',iatom2,' at xred=',xred(:,iatom2),'. Is sent to',symxred2(:),' with symspinat2=',symspinat2(:)
516 : : ! write(std_out,'(a,3f8.4)')' and nucdipmom2=',symnucdipmom2cart(:,1)
517 : :
518 : : ! Check whether there exists an atom of the same class at the
519 : : ! same location, with the correct spinat and nuclear dipole moment circulation
520 [ + + ]: 28461684 : do kk=1,natomcl(iclass)
521 : :
522 : 27780217 : found3=1
523 : 27780217 : iatom3=class(kk,iclass)
524 : : ! Check the location
525 : 27780217 : diffr1=xred(1,iatom3)-symxred2(1)
526 : 27780217 : diff1=diffr1-nint(diffr1)
527 [ + + ]: 27780217 : if(diff1**2>tolsym2)then
528 : : found3=0 ; cycle
529 : : else
530 : 6111514 : diffr2=xred(2,iatom3)-symxred2(2)
531 : 6111514 : diff2=diffr2-nint(diffr2)
532 [ + + ]: 6111514 : if(diff2**2>tolsym2)then
533 : : found3=0 ; cycle
534 : : else
535 : 3310919 : diffr3=xred(3,iatom3)-symxred2(3)
536 : 3310919 : diff3=diffr3-nint(diffr3)
537 [ + + ]: 3310919 : if( (diff1**2+diff2**2+diff3**2) > tolsym**2 )then
538 : : found3=0 ; cycle
539 : : endif
540 : : endif
541 : : endif
542 : : ! Check the spinat
543 [ + + ]: 2693485 : if (nspden/=4) then
544 [ + + ]: 10693884 : diff(:)=spinat(:,iatom3)-symspinat2(:)
545 : : else
546 [ + + ]: 80056 : diff(:)=spinatred(:,iatom3)-symspinat2(:)
547 : : end if
548 [ + + ]: 2693485 : if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 )then
549 : : found3=0
550 : : cycle
551 : : endif
552 : : ! Check the nucdipmom
553 : : ! hand3 gives original circulation sense of nuclear dipole
554 : 2693293 : call acrossb(local_nucdipmom(1:3,2,iatom3),local_nucdipmom(1:3,3,iatom3),hand3)
555 : :
556 : : ! hand2 gives circulation sense of tentative, symmetry equivalent nuclear dipole
557 : 2693293 : call acrossb(symnucdipmom2cart(1:3,2),symnucdipmom2cart(1:3,3),hand2)
558 : :
559 [ + + ]: 10773172 : diff(:)=hand3(:)-hand2(:)
560 [ + + + + : 10766508 : if( any(abs(diff)>tolsym) )found3=0
+ + ]
561 : :
562 : : !diff(:)=symnucdipmom2cart(:,1) - local_nucdipmom(:,1,iatom3)
563 : : !if ( any(abs(diff)>tolsym) ) found3=0
564 : :
565 : 681467 : if(found3==1)exit
566 : : end do ! End loop over iatom3
567 : :
568 [ + + ]: 4186758 : if(found3==0)then
569 : : trialok=0
570 : : exit
571 : : end if
572 : : end do ! End loop over iatom2
573 : :
574 [ + + ]: 2086166 : if(trialok==0)exit
575 : : end do ! End loop over all classes
576 : :
577 : : !write(std_out,*)' For trial isym=',isym,', trialok = ',trialok
578 : : !write(std_out,*)' '
579 : :
580 [ + + ]: 2015528 : if(trialok==1)then
581 : 590341 : nsym=nsym+1
582 [ - + ]: 590341 : if(nsym>msym)then
583 : : write(msg,'(a,i0,2a,i0,4a)')&
584 : 0 : 'The number of symmetries (including non-symmorphic translations) is: ', nsym, ch10,&
585 : 0 : 'is larger than maxnsym: ',msym,ch10,&
586 : 0 : 'Action: increase maxnsym in the input, or take a cell that is primitive, ',ch10,&
587 : 0 : 'or at least smaller than the present one.'
588 : 0 : ABI_ERROR(msg)
589 : : end if
590 : : ntrial=ntrial+1
591 [ + + + + ]: 7674433 : symrel(:,:,nsym)=ptsymrel(:,:,isym)
592 : 590341 : symafm(nsym)=trialafm
593 [ + + ]: 2361364 : tnons(:,nsym)=trialnons(:)-nint(trialnons(:)-tolsym)
594 : : end if
595 : :
596 : : end do ! End the loop on tentative translations
597 : : end do ! End big loop over each symmetry operation of the Bravais lattice
598 : :
599 : : !write(std_out,'(a)')' m_symfind%symfind : after big loop, will call ABI_FREE '; call flush(std_out)
600 : :
601 : 18163 : ABI_FREE(class)
602 : 18163 : ABI_FREE(natomcl)
603 : 18163 : ABI_FREE(chrgat_)
604 : 18163 : ABI_FREE(chrgatcl)
605 : 18163 : ABI_FREE(spinatcl)
606 : 18163 : ABI_FREE(typecl)
607 [ - + ]: 18163 : ABI_FREE(local_nucdipmom)
608 : 18163 : ABI_FREE(nucdipmomcl)
609 [ - + ]: 18163 : ABI_FREE(nucdipmomred)
610 [ + + ]: 18163 : if (nspden==4) then
611 [ - + ]: 546 : ABI_FREE(spinatred)
612 : : end if
613 : :
614 : : !write(std_out,'(a,i6)')' m_symfind%symfind : call sg_multable, nsym= ',nsym; call flush(std_out)
615 : :
616 : : ! The algorithm in sg_multable is still cubic in nsym, so avoid calling it uselessly when nsym is too large
617 [ + + + - ]: 18163 : if(present(ierr) .or. nsym<=384)then
618 : 18163 : call sg_multable(nsym, symafm, symrel, ierr_, tnons=tnons, tnons_tol=tolsym)
619 : : else
620 : 0 : ierr_=0
621 : : endif
622 : :
623 : : !write(std_out,'(a)')' m_symfind%symfind : call print_symmetries, ierr_= ',ierr_; call flush(std_out)
624 : :
625 [ + + + + ]: 18164 : if (ierr_/=0) call print_symmetries([std_out], nsym, symrel, tnons, symafm)
626 : :
627 [ + + ]: 18163 : if(.not.present(ierr))then
628 [ - + ]: 9160 : ABI_CHECK(ierr_==0,"Error in group closure")
629 : : else
630 : 9003 : ierr=ierr_
631 : : endif
632 : :
633 : : !write(msg,'(a,I0,es16.6,a)')' symfind : exit, nsym, tolsym=',nsym,tolsym,ch10
634 : : !write(msg,'(2a)') trim(msg),' symrel matrices, symafm and tnons are :'
635 : : !call wrtout(std_out,msg)
636 : : !do isym=1,nsym
637 : : ! write(msg,'(i4,4x,3i4,2x,3i4,2x,3i4,4x,i4,4x,3f8.4)' ) isym,symrel(:,:,isym),symafm(isym),tnons(:,isym)
638 : : ! call wrtout(std_out,msg)
639 : : !end do
640 : : !stop
641 : :
642 : : !write(std_out,'(a)')' m_symfind%symfind : exit '; call flush(std_out)
643 : :
644 [ + + + + : 54275 : end subroutine symfind
- + ]
645 : : !!***
646 : :
647 : : !!****f* m_symfind/symfind_expert
648 : : !! NAME
649 : : !! symfind_expert
650 : : !!
651 : : !! FUNCTION
652 : : !! Symmetry finder, with an added layer of robustness compared to symfind,
653 : : !! and for which resymmetrization of atomic positions and tnons is needed..
654 : : !! From the symmetries of the Bravais lattice (ptsymrel),
655 : : !! select those that leave invariant the system, and generate
656 : : !! the corresponding tnons vectors and symafm information.
657 : : !! Unlike symfind_expert, does NOT resymmetrize atomic positions and tnons for more robust determination of the symmetries.
658 : : !! The algorithm is explained in T.G. Worlton and J.L. Warren, Comp. Phys. Comm. 3, 88 (1972) [[cite:Worton1972]]
659 : : !!
660 : : !! INPUTS
661 : : !! chrgat(natom) (optional)=target charge for each atom. Not always used, it depends on the value of constraint_kind
662 : : !! invardir_red (optional)=reduced coordinates of an invariant direction (only acting with symrel - not tnons)
663 : : !! invaraxial_red (optional)=reduced coordinates of an axial vector, similar to invardir_red, but transforms with an additional
664 : : !! deteminant factor under symrel operations
665 : : !! invar_z (optional)= if 1, the z direction must stay invariant for all symrel applied ;
666 : : !! if 2, z must stay invariant and also there cannot be any associated tnons along z.
667 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space
668 : : !! msym=default maximal number of symmetries
669 : : !! natom=number of atoms in cell.
670 : : !! nptsym=number of point symmetries of the Bravais lattice
671 : : !! nspden= number of spin-density components. When 4, the three components of spinat are taken into account, instead of only z-component.
672 : : !! nucdipmom(3,natom) (optional) array of nuclear dipole moments
673 : : ! pawspnorb=flag: 1 if spin-orbit coupling is activated
674 : : !! ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
675 : : !! of the Bravais lattice in real space in terms of primitive translations.
676 : : !! spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
677 : : !! tolsym=tolerance for the symmetries
678 : : !! typat(natom)=integer identifying type of atom.
679 : : !! usepaw= 0 for non paw calculation; =1 for paw calculation
680 : : !!
681 : : !! OUTPUT
682 : : !! nsym=actual number of symmetries
683 : : !! symafm(1:msym)=(anti)ferromagnetic part of nsym symmetry operations
684 : : !! symrel(3,3,1:msym)= nsym symmetry operations in real space in terms of primitive translations
685 : : !! tnons(3,1:msym)=nonsymmorphic translations for each symmetry (would be 0 0 0 each for a symmorphic space group)
686 : : !!
687 : : !! SIDE EFFECTS
688 : : !! xred(3,natom)=reduced coordinates of atoms in terms of real space
689 : : !! primitive translations. Might be changed during the resymmetrization.
690 : : !!
691 : : !! SOURCE
692 : :
693 : 9002 : subroutine symfind_expert(gprimd,msym,natom,nptsym,nspden,nsym,&
694 : 9002 : pawspnorb,prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,usepaw,xred,&
695 : : chrgat,nucdipmom,invardir_red,invaraxial_red,invar_z) ! Optional - although for the time being all are required ...
696 : :
697 : : !Arguments ------------------------------------
698 : : !scalars
699 : : integer,intent(in) :: msym,natom,nptsym,nspden,pawspnorb,usepaw
700 : : integer,intent(in) :: prtvol
701 : : integer,optional,intent(in) :: invar_z
702 : : integer,intent(out) :: nsym
703 : : real(dp),intent(in) :: tolsym
704 : : !arrays
705 : : integer,intent(in) :: ptsymrel(3,3,msym),typat(natom)
706 : : integer,intent(inout) :: symafm(msym),symrel(3,3,msym) !vz_i
707 : : real(dp),intent(in) :: gprimd(3,3),spinat(3,natom)
708 : : real(dp),intent(inout) :: xred(3,natom)
709 : : real(dp),optional,intent(in) :: invardir_red(3),invaraxial_red(3),chrgat(natom)
710 : : real(dp),optional, intent(in) :: nucdipmom(3,natom)
711 : : real(dp),intent(inout) :: tnons(3,msym) !vz_i
712 : :
713 : : !Local variables-------------------------------
714 : : !scalars
715 : : integer, save :: print_comment_tolsym=1
716 : : integer :: fixed_mismatch,mismatch_fft_tnons
717 : : integer :: ierr,isym,use_inversion
718 : : character(len=1000) :: msg
719 : : !arrays
720 : 9002 : integer,allocatable :: indsym(:,:,:),symrec(:,:,:)
721 : 9002 : real(dp),allocatable :: tnons_new(:,:)
722 : : !**************************************************************************
723 : :
724 : : ! write(std_out,*)' m_symfind%symfind_expert : enter '
725 : :
726 : 9002 : use_inversion=1
727 [ + + + + : 9002 : if (usepaw == 1 .and. (nspden==4.or.pawspnorb>0)) then
+ + ]
728 : 239 : ABI_COMMENT("Removing inversion and improper rotations from initial space group because of PAW + SOC")
729 : : ! MMignolet: PAW can be used with inversion, however it results in seg faults in the dmft code. To enable when this is fixed...
730 : 239 : use_inversion=0
731 : : end if
732 : :
733 : : ! write(std_out,*)' m_symfind%symfind_expert : before call symfind (1) '
734 : :
735 : : call symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
736 : : prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
737 : 9002 : chrgat=chrgat,nucdipmom=nucdipmom,ierr=ierr,invardir_red=invardir_red,invaraxial_red=invaraxial_red,invar_z=invar_z)
738 : :
739 : : ! write(std_out,*)' m_symfind%symfind_expert : after call symfind (1) '
740 : :
741 : : !If the group closure is not obtained, which should be exceptional, try with a larger tolsym (three times larger)
742 [ + + ]: 9002 : if(ierr/=0)then
743 : 1 : ABI_WARNING('Will try to obtain group closure by using a tripled tolsym.')
744 : : call symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
745 : : prtvol,ptsymrel,spinat,symafm,symrel,tnons,three*tolsym,typat,use_inversion,xred,&
746 : 1 : chrgat=chrgat,nucdipmom=nucdipmom,ierr=ierr,invardir_red=invardir_red,invaraxial_red=invaraxial_red,invar_z=invar_z)
747 [ - + ]: 1 : ABI_CHECK(ierr==0,"Error in group closure")
748 : 1 : ABI_WARNING('Succeeded to obtain group closure by using a tripled tolsym.')
749 : : endif
750 : :
751 : : ! If the tolerance on symmetries is bigger than 1.e-8, symmetrize tnons for gliding or screw operations,
752 : : ! symmetrize the atomic positions and recompute the symmetry operations
753 [ + + ]: 9002 : if(tolsym>1.00001e-8)then
754 : :
755 : 8946 : call symmetrize_tnons(nsym,symrel,tnons,tolsym)
756 [ - + + - : 35784 : ABI_MALLOC(indsym,(4,natom,nsym))
- + + - +
- + - - +
- + ]
757 [ - + + - : 26838 : ABI_MALLOC(symrec,(3,3,nsym))
+ - + - -
+ + - ]
758 [ + + ]: 299428 : do isym=1,nsym
759 : 299428 : call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
760 : : end do
761 : 8946 : call symatm(indsym,natom,nsym,symrec,tnons,tolsym,typat,xred)
762 : 8946 : call symmetrize_xred(natom,nsym,symrel,tnons,xred,indsym=indsym)
763 [ - + ]: 8946 : ABI_FREE(indsym)
764 [ - + ]: 8946 : ABI_FREE(symrec)
765 : :
766 [ + + ]: 8946 : if(print_comment_tolsym==1)then
767 : : write(msg,'(a,es12.3,18a)')&
768 : 1370 : 'The tolerance on symmetries =',tolsym,' is bigger than 1.0e-8.',ch10,&
769 : 1370 : 'In order to avoid spurious effects, the atomic coordinates have been',ch10,&
770 : 1370 : 'symmetrized before storing them in the dataset internal variable.',ch10,&
771 : 1370 : 'So, do not be surprised by the fact that your input variables (xcart, xred, ...)',ch10,&
772 : 1370 : 'do not correspond exactly to the ones echoed by ABINIT, the latter being used to do the calculations.',ch10,&
773 : 1370 : 'This is not a problem per se.',ch10,&
774 : 1370 : 'Still, in order to avoid this symmetrization (e.g. for specific debugging/development),',&
775 : 1370 : ' decrease tolsym to 1.0e-8 or lower,',ch10,&
776 : 1370 : 'or (much preferred) use input primitive vectors that are accurate to better than 1.0e-8.',ch10,&
777 : 2740 : 'This message will only be printed once, even if there are other datasets where tolsym is bigger than 1.0e-8.'
778 : 1370 : ABI_COMMENT(msg)
779 : 1370 : print_comment_tolsym=0
780 : : endif
781 : :
782 : : !write(std_out,*)' m_symfind%symfind_expert : before call symfind (3) '
783 : :
784 : : call symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
785 : : prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
786 : 8946 : chrgat=chrgat,nucdipmom=nucdipmom,invardir_red=invardir_red,invaraxial_red=invaraxial_red,invar_z=invar_z)
787 : :
788 : : ! write(std_out,*)' m_symfind%symfind_expert : after call symfind (3) '
789 : :
790 : : !Needs one more resymmetrization, for the tnons
791 [ - + + - : 26838 : ABI_MALLOC(tnons_new,(3,nsym))
+ - + - -
+ - + ]
792 : :
793 : : call symmetrize_xred(natom,nsym,symrel,tnons,xred,&
794 : 8946 : fixed_mismatch=fixed_mismatch,mismatch_fft_tnons=mismatch_fft_tnons,tnons_new=tnons_new,tolsym=tolsym)
795 [ + + + + ]: 1170874 : tnons(:,1:nsym)=tnons_new(:,:)
796 [ - + ]: 26838 : ABI_FREE(tnons_new)
797 : : end if ! tolsym >1.00001e-8
798 : :
799 : : ! write(std_out,*)' m_symfind%symfind_expert : exit '
800 : :
801 : 9002 : end subroutine symfind_expert
802 : : !!***
803 : :
804 : : !!****f* m_symfind/symanal
805 : : !! NAME
806 : : !! symanal
807 : : !!
808 : : !! FUNCTION
809 : : !! Find the space group, Bravais lattice, including Shubnikov characteristics
810 : : !! from the list of symmetries (including magnetic characteristics), and lattice parameters
811 : : !! Warning: the recognition of the space group might not yet work for the
812 : : !! Shubnikov group of type IV
813 : : !!
814 : : !! INPUTS
815 : : !! chkprim= if 1 then stop if the cell is not primitive
816 : : !! msym=default maximal number of symmetries
817 : : !! nsym=actual number of symmetries
818 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
819 : : !! symafm(1:msym)=(anti)ferromagnetic part of symmetry operations
820 : : !! symrel(3,3,1:msym)=symmetry operations in real space in terms
821 : : !! of primitive translations
822 : : !! tnons(3,1:msym)=nonsymmorphic translations for symmetry operations
823 : : !! tolsym=tolerance for the symmetry operations
824 : : !! [verbose]= if true, will list the symmetry operation labels
825 : : !!
826 : : !! OUTPUT
827 : : !! bravais(11)=characteristics of Bravais lattice (see symlatt.F90)
828 : : !! genafm(3)=magnetic translation generator (in case of Shubnikov group type IV)
829 : : !! ptgroupma = magnetic point group number
830 : : !! spgroup=symmetry space group
831 : : !!
832 : : !! SOURCE
833 : :
834 : 12028 : subroutine symanal(bravais,chkprim,genafm,msym,nsym,ptgroupma,rprimd,spgroup,symafm,symrel,tnons,tolsym, &
835 : : verbose) ! optional
836 : :
837 : : !Arguments ------------------------------------
838 : : !scalars
839 : : integer,intent(in) :: chkprim,msym,nsym
840 : : integer,intent(out) :: ptgroupma,spgroup
841 : : real(dp),intent(in) :: tolsym
842 : : logical,optional,intent(in) :: verbose
843 : : !arrays
844 : : integer,intent(out) :: bravais(11)
845 : : integer,intent(in) :: symafm(msym),symrel(3,3,msym)
846 : : real(dp),intent(in) :: rprimd(3,3), tnons(3,msym)
847 : : real(dp),intent(out) :: genafm(3)
848 : :
849 : : !Local variables-------------------------------
850 : : !scalars
851 : : integer, parameter :: maxsym=192
852 : : ! In this routine, maxsym is used either to determine the ptsymrel from the rprimd (routine symlatt),
853 : : ! so, it might be up to 192 = 4*48 for FCC, and also to define the maximum number of symmetry operation labels,
854 : : ! but only in case the cell is primitive, which gives the same upper bound. Thus in this routine, msym might
855 : : ! be equal to nsym.
856 : : integer :: iholohedry_nomagn,isym,isym_nomagn,multi, nptsym,nsym_nomagn,shubnikov
857 : : logical :: verbose_
858 : : character(len=5) :: ptgroup,ptgroupha
859 : : character(len=500) :: msg
860 : : !arrays
861 : : integer :: identity(3,3)
862 : 12028 : integer,allocatable :: ptsymrel(:,:,:),symrel_nomagn(:,:,:)
863 : 12028 : real(dp),allocatable :: tnons_nomagn(:,:)
864 : : character(len=128) :: labels(maxsym)
865 : : ! *************************************************************************
866 : :
867 : : !write(std_out,*)' symanal : enter'
868 : : !write(std_out,*)' symanal : chkprim =',chkprim
869 : : !write(std_out,*)' symanal : nsym=',nsym
870 : : !do isym=1,nsym
871 : : ! write(std_out,*)' symanal : symrel=',symrel(1:3,1:3,isym)
872 : : !enddo
873 : :
874 : 12028 : verbose_=.false.
875 [ + + ]: 12028 : if (present(verbose)) verbose_=verbose
876 : :
877 : : !This routine finds the Bravais characteristics, without actually looking at the symmetry operations.
878 [ - + ]: 12028 : ABI_MALLOC(ptsymrel,(3,3,maxsym))
879 : 12028 : call symlatt(bravais,dev_null,maxsym,nptsym,ptsymrel,rprimd,tolsym)
880 [ - + ]: 12028 : ABI_FREE(ptsymrel)
881 : :
882 : : ! Check whether the cell is primitive or not.
883 : 12028 : call chkprimit(chkprim,multi,nsym,symafm,symrel)
884 : :
885 : 12028 : spgroup=0 ; ptgroupma=0 ; genafm(:)=zero
886 : :
887 [ + + ]: 12028 : if (multi>1) then ! Modify bravais if the cell is not primitive ; no determination of the space group
888 : 289 : bravais(1)=-bravais(1)
889 : : else
890 : :
891 : : ! The cell is primitive, so that the space group can be
892 : : ! determined. Need to distinguish Fedorov and Shubnikov groups.
893 : : ! Do not distinguish Shubnikov types I and II.
894 : : ! Also identify genafm, in case of Shubnikov type IV
895 : 11739 : identity(:,:)=reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/))
896 : 11739 : shubnikov=1
897 [ + + ]: 299029 : do isym=1,nsym
898 [ + + ]: 299029 : if(symafm(isym)==-1)then
899 : 5756 : shubnikov=3
900 [ + + + + : 74828 : if(sum(abs(symrel(:,:,isym)-identity(:,:)))==0)then
+ + ]
901 : 611 : shubnikov=4
902 [ + + ]: 2444 : genafm(:)=tnons(:,isym)
903 : : !write(std_out,*)' isym=',isym
904 : : !write(std_out,*)' symrel(:,:,isym)',symrel(:,:,isym)
905 : : !write(std_out,*)' tnons(:,isym)',tnons(:,isym)
906 : : !write(std_out,*)' symafm(isym)',symafm(isym)
907 : : exit
908 : : end if
909 : : end if
910 : : end do
911 : :
912 [ + + ]: 11739 : if(shubnikov/=1)then
913 [ + + ]: 1515 : if(shubnikov==3)write(msg, '(a)' )' Shubnikov space group type III'
914 [ + + ]: 1515 : if(shubnikov==4)write(msg, '(a)' )' Shubnikov space group type IV'
915 : 1515 : call wrtout(std_out,msg)
916 : : end if
917 : :
918 [ + + ]: 11739 : if(shubnikov==1 .or. shubnikov==3)then
919 : : ! Find the correct Bravais characteristics and point group
920 : : ! Should also be used for Shubnikov groups of type IV ...
921 : 11128 : call symbrav(bravais,msym,nsym,ptgroup,rprimd,symrel,tolsym)
922 : :
923 : : ! Find the space group
924 : 11128 : call symspgr(bravais,labels,nsym,spgroup,symrel,tnons,tolsym)
925 : :
926 [ + + ]: 11128 : if (verbose_) then
927 [ + + ]: 2275 : do isym=1,nsym
928 [ + - ]: 2203 : write(msg,'(a,i3,2a)')' symanal : the symmetry operation no. ',isym,' is ',trim(labels(isym))
929 : 2275 : call wrtout(std_out,msg)
930 : : enddo
931 : : endif
932 : :
933 : : end if
934 : :
935 [ + + ]: 11128 : if (shubnikov/=1) then
936 : :
937 : : ! Determine nonmagnetic symmetry operations
938 : 1515 : nsym_nomagn=nsym/2
939 [ - + + - : 4545 : ABI_MALLOC(symrel_nomagn,(3,3,nsym_nomagn))
+ - + - -
+ - + ]
940 [ - + + - : 4545 : ABI_MALLOC(tnons_nomagn,(3,nsym_nomagn))
+ - + - -
+ - + ]
941 : 1515 : isym_nomagn=0
942 [ + + ]: 23791 : do isym=1,nsym
943 [ + + ]: 23791 : if(symafm(isym)==1)then
944 : 11138 : isym_nomagn=isym_nomagn+1
945 [ + + + + ]: 144794 : symrel_nomagn(:,:,isym_nomagn)=symrel(:,:,isym)
946 [ + + ]: 44552 : tnons_nomagn(:,isym_nomagn)=tnons(:,isym)
947 : : end if
948 : : end do
949 : :
950 [ + + ]: 1515 : if (shubnikov==3) then
951 : :
952 : : ! write(std_out,*)' symanal : will enter symbrav with halved symmetry set'
953 : : ! write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
954 : : ! do isym=1,nsym_nomagn
955 : : ! write(std_out,'(i3,2x,9i3,3es12.2,i3)')isym,symrel_nomagn(:,:,isym),tnons_nomagn(:,isym)
956 : : ! end do
957 : :
958 : : ! Find the point group of the halved symmetry set
959 : 904 : call symptgroup(iholohedry_nomagn,nsym_nomagn,ptgroupha,symrel_nomagn)
960 : :
961 : : ! Deduce the magnetic point group (ptgroupma) from ptgroup and ptgroupha
962 : 904 : call getptgroupma(ptgroup,ptgroupha,ptgroupma)
963 : :
964 : 611 : else if(shubnikov==4)then
965 : :
966 : : ! Find the Fedorov space group of the halved symmetry set
967 : 611 : call symspgr(bravais,labels,nsym_nomagn,spgroup,symrel_nomagn,tnons_nomagn,tolsym)
968 : :
969 : : ! The magnetic translation generator genafm has already been determined
970 : : ! write(std_out,*)' genafm =',genafm, ' spgroup=',spgroup
971 [ - + ]: 611 : if (verbose_) then
972 : 0 : write(msg, '(a)' )' Select only the non-magnetic symmetry operations '
973 : 0 : call wrtout(std_out,msg)
974 : :
975 [ # # ]: 0 : do isym=1,nsym
976 [ # # ]: 0 : if(symafm(isym)==1)then
977 : 0 : isym_nomagn=isym_nomagn+1
978 [ # # ]: 0 : write(msg,'(a,i3,2a)')' symspgr : the symmetry operation no. ',isym,' is ',trim(labels(isym_nomagn))
979 : 0 : call wrtout(std_out,msg)
980 : : endif
981 : : enddo
982 : : endif
983 : :
984 : : end if
985 : :
986 [ - + ]: 1515 : ABI_FREE(symrel_nomagn)
987 [ - + ]: 1515 : ABI_FREE(tnons_nomagn)
988 : : end if ! Shubnikov groups
989 : :
990 : : end if
991 : :
992 : 12028 : end subroutine symanal
993 : : !!***
994 : :
995 : : !!****f* m_symfind/symbrav
996 : : !! NAME
997 : : !! symbrav
998 : : !!
999 : : !! FUNCTION
1000 : : !! From the list of symmetry operations, and the lattice vectors,
1001 : : !! determine the Bravais information (including the holohedry, the centering,
1002 : : !! the coordinate of the primitive vectors in the conventional vectors),
1003 : : !! as well as the point group.
1004 : : !!
1005 : : !! INPUTS
1006 : : !! msym=dimension of symrel
1007 : : !! nsym=actual number of symmetries
1008 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1009 : : !! symrel(3,3,msym)=symmetry operations in real space in terms
1010 : : !! of primitive translations
1011 : : !! tolsym=tolerance for the symmetries
1012 : : !!
1013 : : !! OUTPUT
1014 : : !! bravais(11): bravais(1)=iholohedry
1015 : : !! bravais(2)=center
1016 : : !! bravais(3:11)=coordinates of rprimd in the axes
1017 : : !! of the conventional bravais lattice (*2 if center/=0)
1018 : : !! ptgroup=symmetry point group
1019 : : !! [axis(3)]=Invariant axis in the conventional vector coordinates
1020 : : !! Set to (/0,0,0/) if the lattice belongs to the same holohedry as the lattice+atoms (+electric field + ...).
1021 : : !!
1022 : : !! SOURCE
1023 : :
1024 : 11128 : subroutine symbrav(bravais,msym,nsym,ptgroup,rprimd,symrel,tolsym,axis)
1025 : :
1026 : : !Arguments ------------------------------------
1027 : : !scalars
1028 : : integer,intent(in) :: msym,nsym
1029 : : real(dp),intent(in) :: tolsym
1030 : : character(len=5),intent(out) :: ptgroup
1031 : : !arrays
1032 : : integer,intent(in) :: symrel(3,3,msym)
1033 : : integer,optional,intent(out) :: axis(3)
1034 : : integer,intent(out) :: bravais(11)
1035 : : real(dp),intent(in) :: rprimd(3,3)
1036 : :
1037 : : !Local variables-------------------------------
1038 : : !scalars
1039 : : integer :: iaxis,ii,bravais1now,ideform,iholohedry,invariant,isym
1040 : : integer :: jaxis,next_stage,nptsym,problem,maxsym
1041 : : integer, parameter :: naxes_ortho=22, naxes_hexa=7
1042 : : real(dp) :: norm,scprod
1043 : : character(len=500) :: msg
1044 : : !arrays
1045 : : integer :: identity(3,3),axis_trial(3),hexa_axes(3,naxes_hexa),ortho_axes(3,naxes_ortho)
1046 : 11128 : integer,allocatable :: ptsymrel(:,:,:),symrelconv(:,:,:)
1047 : : real(dp) :: axes(3,3),axis_cart(3),axis_red(3)
1048 : : real(dp) :: rprimdconv(3,3),rprimdtry(3,3),rprimdnow(3,3)
1049 : : real(dp) :: rprimdconv_invt(3,3)
1050 : : !**************************************************************************
1051 : :
1052 : : !DEBUG
1053 : : !write(std_out,*)' symbrav : enter '
1054 : : !call flush(std_out)
1055 : : !ENDDEBUG
1056 : :
1057 : : identity(:,:)=0
1058 : 11128 : identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
1059 : :
1060 : 11128 : ortho_axes(:,:)=0
1061 : 11128 : ortho_axes(1,1)=1
1062 : 11128 : ortho_axes(2,2)=1
1063 : 11128 : ortho_axes(3,3)=1
1064 [ + + ]: 44512 : ortho_axes(:,4)=(/0,1,1/)
1065 [ + + ]: 44512 : ortho_axes(:,5)=(/1,0,1/)
1066 [ + + ]: 44512 : ortho_axes(:,6)=(/1,1,0/)
1067 [ + + ]: 44512 : ortho_axes(:,7)=(/0,1,-1/)
1068 [ + + ]: 44512 : ortho_axes(:,8)=(/-1,0,1/)
1069 [ + + ]: 44512 : ortho_axes(:,9)=(/1,-1,0/)
1070 [ + + ]: 44512 : ortho_axes(:,10)=(/0,1,2/)
1071 [ + + ]: 44512 : ortho_axes(:,11)=(/2,0,1/)
1072 [ + + ]: 44512 : ortho_axes(:,12)=(/1,2,0/)
1073 [ + + ]: 44512 : ortho_axes(:,13)=(/1,1,1/)
1074 [ + + ]: 44512 : ortho_axes(:,14)=(/-1,1,1/)
1075 [ + + ]: 44512 : ortho_axes(:,15)=(/1,-1,1/)
1076 [ + + ]: 44512 : ortho_axes(:,16)=(/1,1,-1/)
1077 [ + + ]: 44512 : ortho_axes(:,17)=(/2,1,1/)
1078 [ + + ]: 44512 : ortho_axes(:,18)=(/1,2,1/)
1079 [ + + ]: 44512 : ortho_axes(:,19)=(/1,1,2/)
1080 [ + + ]: 44512 : ortho_axes(:,20)=(/2,1,-1/)
1081 [ + + ]: 44512 : ortho_axes(:,21)=(/-1,2,1/)
1082 [ + + ]: 44512 : ortho_axes(:,22)=(/1,-1,2/)
1083 : :
1084 : 11128 : hexa_axes(:,:)=0
1085 : 11128 : hexa_axes(1,1)=1
1086 : 11128 : hexa_axes(2,2)=1
1087 : 11128 : hexa_axes(3,3)=1
1088 [ + + ]: 44512 : hexa_axes(:,4)=(/1,-1,0/)
1089 [ + + ]: 44512 : hexa_axes(:,5)=(/2,1,0/)
1090 [ + + ]: 44512 : hexa_axes(:,6)=(/1,1,0/)
1091 [ + + ]: 44512 : hexa_axes(:,7)=(/1,2,0/)
1092 : :
1093 : : !Determine the point group from the list of symmetry operations.
1094 : : !Also determine the holohedry, up to one undeterminacy : hR versus hP
1095 : 11128 : call symptgroup(iholohedry,nsym,ptgroup,symrel)
1096 : :
1097 : : !DEBUG
1098 : : !write(std_out,*)' symbrav, after symptgroup: nsym=',nsym
1099 : : !call flush(std_out)
1100 : : !write(std_out,*)' symbrav: symrel='
1101 : : !do isym=1,nsym
1102 : : ! write(std_out,'(9i4)')symrel(:,:,isym)
1103 : : !enddo
1104 : : !write(std_out,*)' symbrav: iholohedry=',iholohedry
1105 : : !call flush(std_out)
1106 : : !ENDDEBUG
1107 : :
1108 : : !Loop over trial deformations
1109 : : !This is needed in case the Bravais lattice determination from the lattice vectors
1110 : : !has a higher holohedry than the real one, in which the symmetry
1111 : : !operations for the atoms (or electric field, etc) are taken into account
1112 : 11128 : iaxis=0
1113 : 11128 : invariant=0
1114 : 11128 : next_stage=0
1115 : 11128 : rprimdnow(:,:)=rprimd(:,:)
1116 : 11128 : rprimdtry(:,:)=rprimd(:,:)
1117 [ - + + - : 33384 : ABI_MALLOC(symrelconv,(3,3,nsym))
+ - + - -
+ + - ]
1118 : :
1119 : : !At most will have to try naxes_ortho*5 deformations (naxes_ortho axes, five stages)
1120 : : !First, test whether the current recognition of Bravais lattice is problematic (iholohedry differs from bravais(1)).
1121 : : !Then, if there is a problem, test different deformations of rprimd, one after the other.
1122 : : !For each try, there is a new rprimdtry from the different set of deformation, that is generated later in the loop
1123 : : !Also bravais1now and rprimdnow might changed (and progressively lowered).
1124 : : !The latter change induces at most 5 stages for the computation (cubic->tetragonal->orthorhombic->monoclinic->triclinic).
1125 : : !After an upgrade of bravais1now and rprimdnow, one has to restart the full set of deformations.
1126 : : !Not sure that this procedure resolves all cases, but seems to work on >40000 inaccurate POSCAR files.
1127 [ + - ]: 15626 : do ideform=1,naxes_ortho*5
1128 : :
1129 : 15626 : maxsym=max(192,msym)
1130 [ - + + - : 46878 : ABI_MALLOC(ptsymrel,(3,3,maxsym))
+ - - + -
+ ]
1131 : 15626 : call symlatt(bravais,std_out,maxsym,nptsym,ptsymrel,rprimdtry,tolsym)
1132 [ - + ]: 15626 : ABI_FREE(ptsymrel)
1133 : :
1134 : : !DEBUG
1135 : : !write(std_out,*)' symbrav: inside loop with ideform,iaxis=',ideform,iaxis
1136 : : !write(std_out,'(a,9f12.6)')' rprimdtry=',rprimdtry(:,:)
1137 : : !write(std_out,'(a,2i4)')' bravais(1:2)=',bravais(1:2)
1138 : : !call flush(std_out)
1139 : : !ENDDEBUG
1140 : :
1141 : :
1142 : : ! Examine the agreement with bravais(1)
1143 : : ! Warning : might change Bravais lattice hR to hP, if hexagonal axes
1144 : 15626 : problem=0
1145 : 23617 : select case (bravais(1))
1146 : : case (7)
1147 [ + + ]: 7991 : if(iholohedry<6)problem=1
1148 [ - + ]: 7991 : if(iholohedry==6)problem=2
1149 : : case (6)
1150 [ + + ]: 669 : if(iholohedry<4)problem=1
1151 [ - + ]: 669 : if(iholohedry==7 .or. iholohedry==4)problem=2
1152 : : ! Here, change hR into hP
1153 [ + + ]: 669 : if(iholohedry==5)iholohedry=6
1154 : : case (5)
1155 [ + + ]: 806 : if(iholohedry<4)problem=1
1156 [ - + ]: 806 : if(iholohedry==7 .or. iholohedry==6 .or. iholohedry==4)problem=2
1157 : : case (4)
1158 [ + + ]: 2042 : if(iholohedry<4)problem=1
1159 [ - + ]: 2042 : if(iholohedry>4)problem=2
1160 : : case (3)
1161 [ + + ]: 898 : if(iholohedry<3)problem=1
1162 [ - + ]: 898 : if(iholohedry>3)problem=2
1163 : : case (2)
1164 [ + + ]: 1775 : if(iholohedry<2)problem=1
1165 [ - + ]: 1775 : if(iholohedry>2)problem=2
1166 : : case (1)
1167 [ + + + + : 15626 : if(iholohedry>1)problem=2
+ + + - -
+ ]
1168 : : end select
1169 : :
1170 : : ! This is the usual situation, in which the lattice belong to the same holohedry
1171 : : ! as the lattice+atoms (+electric field + ...)
1172 [ + + ]: 15626 : if(problem==0)exit
1173 : :
1174 [ - + ]: 4498 : if(problem==2)then
1175 [ # # ]: 0 : if(iaxis==0)then
1176 : : write(msg, '(3a,i3,3a,i3,7a)' )&
1177 : 0 : 'The Bravais lattice determined only from the primitive',ch10,&
1178 : 0 : 'vectors (rprim or angdeg), bravais(1)=',bravais(1),', is not compatible',ch10,&
1179 : 0 : 'with the real one, iholohedry=',iholohedry,', obtained by taking into',ch10,&
1180 : 0 : 'account the symmetry operations. This might be due to an insufficient',ch10,&
1181 : 0 : 'number of digits in the specification of rprim (at least 10),',ch10,&
1182 : 0 : 'or to an erroneous rprim or angdeg. If this is not the case, then ...'
1183 : 0 : ABI_BUG(msg)
1184 : : end if
1185 [ # # ]: 0 : if(iaxis==1)then
1186 : : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1187 : 0 : 'Could not succeed to determine the bravais lattice',ch10,&
1188 : 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1189 : 0 : 'bravais(1)=',bravais(1),ch10,&
1190 : 0 : 'iholohedry=',iholohedry
1191 : 0 : ABI_BUG(msg)
1192 : : end if
1193 : : ! Try to increase tolsym to find the Bravais lattice.
1194 : : maxsym=max(192,msym)
1195 [ # # # # : 0 : ABI_MALLOC(ptsymrel,(3,3,maxsym))
# # # # ]
1196 : : !DEBUG
1197 : : ! write(6,*)' symbrav : will call symlatt, 3*tolsym=',3*tolsym
1198 : : !ENDDEBUG
1199 : 0 : call symlatt(bravais,std_out,maxsym,nptsym,ptsymrel,rprimdtry,3*tolsym)
1200 [ # # ]: 0 : ABI_FREE(ptsymrel)
1201 [ # # ]: 0 : if(bravais(1)==iholohedry)then
1202 : : ! Succeeded
1203 : : exit
1204 : : else
1205 : : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1206 : 0 : 'Could not succeed to determine the bravais lattice, even after considering a larger tolsym',ch10,&
1207 : 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1208 : 0 : 'bravais(1)=',bravais(1),ch10,&
1209 : 0 : 'iholohedry=',iholohedry
1210 : 0 : ABI_BUG(msg)
1211 : : end if
1212 : : end if
1213 : :
1214 [ + - ]: 4498 : if(problem==1)then ! One is left with the problem=1 case, basically iholohedry is lower than bravais(1)
1215 [ + + ]: 4498 : if(iaxis==0)then
1216 : : write(msg, '(a,a,a,i3,a,a,a,i3,a,a,a)' )&
1217 : 2934 : 'The Bravais lattice determined only from the primitive',ch10,&
1218 : 2934 : 'vectors, bravais(1)=',bravais(1),', is more symmetric',ch10,&
1219 : 2934 : 'than the real one, iholohedry=',iholohedry,', obtained by taking into',ch10,&
1220 : 5868 : 'account the atomic positions. Start deforming the primitive vector set.'
1221 : 2934 : ABI_COMMENT(msg)
1222 : : next_stage=1
1223 : : else if(iaxis/=0)then
1224 [ + + ]: 1564 : if(bravais(1)<bravais1now)then
1225 : : write(msg, '(3a,i3,3a,i3,2a)' )&
1226 : 1349 : 'The Bravais lattice determined from modified primitive',ch10,&
1227 : 1349 : 'vectors, bravais(1)=',bravais(1),', has a lower symmetry than before,',ch10,&
1228 : 1349 : 'but is still more symmetric than the real one, iholohedry=',iholohedry,ch10,&
1229 : 2698 : 'obtained by taking into account the atomic positions.'
1230 : 1349 : ABI_COMMENT(msg)
1231 : : next_stage=1
1232 [ - + ]: 215 : else if(iaxis==1)then
1233 : : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1234 : 0 : 'Could not succeed to determine the bravais lattice',ch10,&
1235 : 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1236 : 0 : 'bravais(1)=',bravais(1),ch10,&
1237 : 0 : 'iholohedry=',iholohedry
1238 : 0 : ABI_BUG(msg)
1239 : : end if
1240 : : end if
1241 : : end if ! problem==1
1242 : :
1243 : : ! One is here when problem=1 (iholohedry < bravais(1)) and either
1244 : : ! - iaxis=0 (no deformation has been tried yet),
1245 : : ! - some deformation iaxis has been tried giving bravais(1), but iholohedry < bravais(1) < bravais1now
1246 : : ! - some deformation iaxis has been tried giving bravais(1), but iholohedry < bravais(1) = bravais1now and iaxis/=1 .
1247 : : ! Also, note that next_stage is still 0 when bravais(1)=bravais1now .
1248 : : ! The loop has been ended (so, the search failed) when bravais(1)=bravais1now and iaxis==1.
1249 : :
1250 : : if(next_stage==1)then
1251 : 4283 : bravais1now=bravais(1)
1252 : 4283 : rprimdnow(:,:)=rprimdtry(:,:)
1253 : : ! Generate the symmetry operations in the conventional vector coordinates
1254 [ + + ]: 17132 : rprimdconv(:,1)=bravais(3:5)
1255 [ + + ]: 17132 : rprimdconv(:,2)=bravais(6:8)
1256 [ + + ]: 17132 : rprimdconv(:,3)=bravais(9:11)
1257 : 4283 : axes(:,:)=zero
1258 : 4283 : axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one
1259 [ + + + + : 211269 : symrelconv(:,:,1:nsym)=symrel(:,:,1:nsym)
+ + ]
1260 : 4283 : call symrelrot(nsym,rprimdconv,axes,symrelconv,tolsym)
1261 [ + + ]: 4283 : if(bravais(1)/=6)then
1262 : 4266 : iaxis=naxes_ortho+1
1263 : : else
1264 : 17 : iaxis=naxes_hexa+1
1265 : : end if
1266 : : next_stage=0
1267 : : !DEBUG
1268 : : ! write(std_out,*)' symbrav: next stage, bravais(1), bravais(2) and symrelconv'
1269 : : ! write(std_out,'(a,2i4)')' bravais(1:2)=',bravais(1:2)
1270 : : ! write(std_out,'(a,9f12.6)')' rprimdconv=',rprimdconv(:,:)
1271 : : ! do isym=1,nsym
1272 : : ! write(std_out,'(9i4)')symrelconv(:,:,isym)
1273 : : ! enddo
1274 : : ! call flush(std_out)
1275 : : !ENDDEBUG
1276 : : end if
1277 : :
1278 : : ! Go to the next iaxis that will be left invariant
1279 : 4498 : iaxis=iaxis-1
1280 [ + - ]: 25008 : do jaxis=iaxis,1,-1
1281 [ + + ]: 25008 : if(bravais(1)/=6)then
1282 [ + + ]: 99956 : axis_trial(:)=ortho_axes(:,jaxis)
1283 : : else
1284 [ + + ]: 76 : axis_trial(:)=hexa_axes(:,jaxis)
1285 : : end if
1286 : : ! DEBUG
1287 : : ! write(std_out,*)' symbrav : ixaxis, trial jaxis=',iaxis,jaxis
1288 : : ! write(std_out,*)' axis_trial=',axis_trial
1289 : : ! ENDDEBUG
1290 : 25008 : invariant=1
1291 : : ! Examine whether all symmetry operations leave the axis invariant (might be reversed, though)
1292 [ + + ]: 225362 : do isym=1,nsym
1293 [ - + + + : 1803186 : if(sum(abs(matmul(symrelconv(:,:,isym),axis_trial)+(-axis_trial(:))))/=0 .and. &
+ + + + ]
1294 [ + + + + : 7037398 : sum(abs(matmul(symrelconv(:,:,isym),axis_trial)+axis_trial(:)))/=0 )invariant=0
+ + + + +
+ + + + -
+ - ]
1295 : : end do
1296 [ + + ]: 25008 : if(invariant==1)then
1297 : 4498 : iaxis=jaxis
1298 : : !DEBUG
1299 : : ! write(msg, '(2a,i3)' )ch10,' symbrav : found invariant axis, jaxis=',jaxis
1300 : : ! call wrtout(std_out,msg)
1301 : : !ENDDEBUG
1302 : 4498 : exit
1303 : : end if
1304 : : end do
1305 : :
1306 [ - + ]: 4498 : if(invariant==0)then
1307 : : ! Not a single axis was invariant with respect to all operations ?!
1308 : : ! do isym=1,nsym; write(std_out, '(a,10i4)' )' isym,symrelconv=',isym,symrelconv(:,:,isym); enddo
1309 : : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1310 : 0 : 'Could not succeed to determine the bravais lattice (not a single invariant)',ch10,&
1311 : 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1312 : 0 : 'bravais(1)=',bravais(1),ch10,&
1313 : 0 : 'iholohedry=',iholohedry
1314 : 0 : ABI_BUG(msg)
1315 : : end if
1316 : :
1317 : 4498 : call matr3inv(rprimdconv,rprimdconv_invt)
1318 : : axis_red(:)=axis_trial(1)*rprimdconv_invt(1,:)+ &
1319 : : & axis_trial(2)*rprimdconv_invt(2,:)+ &
1320 [ + + ]: 17992 : & axis_trial(3)*rprimdconv_invt(3,:)
1321 : : axis_cart(:)=axis_red(1)*rprimdnow(:,1)+ &
1322 : : & axis_red(2)*rprimdnow(:,2)+ &
1323 [ + + ]: 17992 : & axis_red(3)*rprimdnow(:,3)
1324 [ + + ]: 17992 : norm=sum(axis_cart(:)**2)
1325 : : !DEBUG
1326 : : ! write(6,*)' axis_trial =',axis_trial
1327 : : ! write(6,*)' axis_red =',axis_red
1328 : : ! write(6,*)' axis_cart =',axis_cart
1329 : : ! write(6,*)' rprimdnow=',rprimdnow
1330 : : !ENDDEBUG
1331 : : ! Expand by a uniform, quite arbitrary, dilatation, along the invariant axis
1332 : : ! Note : make these dilatation different, according to ideform
1333 : : ! XG 20151221 : Still, the interplay between the size of the deformation and the tolsym is not easy to address.
1334 : : ! Indeed the deformation must be sufficiently large to be perceived by symlatt as a real breaking of the
1335 : : ! symmetry of the lattice. In order to deal with all the small values od tolsym, it has been set at a minimum of tol3,
1336 : : ! but it must also be larger than tolsym. Moreover, for some axis choice, the deformation is not aligned with the axis, decreasing
1337 : : ! the effective deformation length. An additional factor of three is thus included, actually increased to six just to be sure...
1338 [ + + ]: 44746 : do ii=1,3
1339 : 13494 : scprod=axis_cart(1)*rprimdnow(1,ii)+axis_cart(2)*rprimdnow(2,ii)+axis_cart(3)*rprimdnow(3,ii)
1340 [ + + ]: 58474 : rprimdtry(:,ii)=rprimdnow(:,ii)+ideform*(max(tol3,six*tolsym)-tol6)*scprod/norm*axis_cart(:)
1341 : : end do
1342 : :
1343 : : !DEBUG
1344 : : ! write(6,*)' rprimdtry=',rprimdtry
1345 : : !ENDDEBUG
1346 : :
1347 : : end do ! ideform
1348 : :
1349 [ - + ]: 11128 : if(bravais(1)/=iholohedry)then
1350 : : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1351 : 0 : 'Despite efforts, Could not succeed to determine the bravais lattice :',ch10,&
1352 : 0 : 'bravais(1)=',bravais(1),ch10,&
1353 : 0 : 'iholohedry=',iholohedry
1354 : 0 : ABI_BUG(msg)
1355 : : end if
1356 : :
1357 [ - + ]: 11128 : ABI_FREE(symrelconv)
1358 : :
1359 [ - + ]: 11128 : if (PRESENT(axis)) then ! Return symmetry axis.
1360 : 0 : axis=(/0,0,0/)
1361 [ # # ]: 0 : if (iaxis/=0) then
1362 [ # # ]: 0 : if(bravais(1)/=6)then
1363 [ # # ]: 0 : axis=ortho_axes(:,iaxis)
1364 : : else
1365 [ # # ]: 0 : axis=hexa_axes(:,iaxis)
1366 : : end if
1367 : : end if
1368 : : end if
1369 : :
1370 : : !DEBUG
1371 : : !write(std_out,'(a)')' symbrav : exit '
1372 : : !ENDDEBUG
1373 : :
1374 : 11128 : end subroutine symbrav
1375 : : !!***
1376 : :
1377 : : !!****f* m_symfind/symspgr
1378 : : !! NAME
1379 : : !! symspgr
1380 : : !!
1381 : : !! FUNCTION
1382 : : !! Find the type of each symmetry operation (calling symcharac):
1383 : : !! proper symmetries 1,2,2_1,3,3_1,3_2,4,4_1,4_2,4_3,6,6_1,...6_5
1384 : : !! improper symmetries -1,m,a,b,c,d,n,g,-3,-4,-6 ,
1385 : : !! Then, build an array with the number of such operations.
1386 : : !! Then, call symlist to identify the space group.
1387 : : !!
1388 : : !! INPUTS
1389 : : !! bravais(11): bravais(1)=iholohedry
1390 : : !! bravais(2)=center
1391 : : !! bravais(3:11)=coordinates of rprimd in the axes
1392 : : !! of the conventional bravais lattice (*2 if center/=0)
1393 : : !! nsym=actual number of symmetries
1394 : : !! symrel(3,3,nsym)= nsym symmetry operations in real space in terms
1395 : : !! of primitive translations
1396 : : !! tnons(3,nsym)=nonsymmorphic translations for each symmetry (would
1397 : : !! be 0 0 0 each for a symmorphic space group)
1398 : : !!
1399 : : !! OUTPUT
1400 : : !! labels(maxsym=192)= labels of the symmetry operations
1401 : : !! spgroup=symmetry space group number
1402 : : !!
1403 : : !! NOTES
1404 : : !! It is assumed that the symmetry operations will be entered in the
1405 : : !! symrel tnons arrays, for the PRIMITIVE cell. The matrix of transformation
1406 : : !! from the primitive cell to the conventional cell is described
1407 : : !! in the array "bravais" (see symlatt.F90).
1408 : : !! The present routine first make the transformation from the
1409 : : !! primitive coordinates to the conventional ones, then eventually
1410 : : !! generate additional symmetries, taking into account the
1411 : : !! centering translations.
1412 : : !! Then, the order and determinant of each symmetry operation
1413 : : !! is determined.
1414 : : !!
1415 : : !! For proper symmetries (rotations), the
1416 : : !! associated translation is also determined.
1417 : : !! However, left or right handed screw rotations are
1418 : : !! not (presently) distinguished, and will be attributed equally
1419 : : !! to left or right.
1420 : : !!
1421 : : !! For the detailed description of the labelling of the axes, see symaxes.f and symplanes.f
1422 : : !!
1423 : : !! SOURCE
1424 : :
1425 : 11739 : subroutine symspgr(bravais,labels,nsym,spgroup,symrel,tnons,tolsym)
1426 : :
1427 : : use m_numeric_tools, only : OPERATOR(.x.)
1428 : :
1429 : : !Arguments ------------------------------------
1430 : : !scalars
1431 : : integer,intent(in) :: nsym
1432 : : integer,intent(out) :: spgroup
1433 : : real(dp),intent(in) :: tolsym
1434 : : !arrays
1435 : : integer,intent(in) :: bravais(11),symrel(3,3,nsym)
1436 : : real(dp),intent(in) :: tnons(3,nsym)
1437 : : character(len=128),intent(out) :: labels(192) ! 192 = maxsym
1438 : :
1439 : : !Local variables-------------------------------
1440 : : !scalars
1441 : : ! logical,parameter :: verbose=.FALSE.
1442 : : integer :: additional_info,brvltt,center,direction=0,found,iholohedry,ii
1443 : : integer :: ishift,isym,jj,nshift,nsymconv,spgaxor,spgorig,sporder
1444 : : character(len=1) :: brvsb
1445 : : character(len=15) :: intsb,ptintsb,ptschsb,schsb
1446 : : character(len=35) :: intsbl
1447 : : character(len=500) :: msg
1448 : : !arrays
1449 : : integer :: ivec1(3), ivec2(3)
1450 : : integer :: n_axes(31),n_axest(31),prime(5),test_direction(3),symrel_uni(3,3)
1451 : : integer :: uniaxis(3),uniaxis_try(3)
1452 : 11739 : integer,allocatable :: determinant(:),symrelconv(:,:,:),t_axes(:)
1453 : : real(dp) :: axes(3,3),rprimdconv(3,3),vect(3,3)
1454 : 11739 : real(dp),allocatable :: shift(:,:),tnonsconv(:,:)
1455 : : !**************************************************************************
1456 : :
1457 : : DBG_ENTER("COLL")
1458 : :
1459 : : !Initialize brvltt, from bravais(2) and bravais(1)
1460 : 11739 : center=bravais(2)
1461 : 11739 : iholohedry=bravais(1)
1462 : 11739 : brvltt=1
1463 [ + + ]: 11739 : if(center==-1)brvltt=2 ! Inner centering
1464 [ + + ]: 11739 : if(center==-3)brvltt=3 ! Face centering
1465 [ + + ]: 11739 : if(center==1)brvltt=5 ! A-Face centering
1466 [ + + ]: 11739 : if(center==2)brvltt=6 ! B-Face centering
1467 [ + + ]: 11739 : if(center==3)brvltt=4 ! C-Face centering
1468 [ + + ]: 11739 : if(iholohedry==5)brvltt=7 ! Rhombohedral
1469 : :
1470 : : !Produce the symmetry operations, in the axis of the conventional cell
1471 : 11739 : nsymconv=nsym
1472 [ + + ]: 11739 : if(center/=0)nsymconv=2*nsymconv
1473 [ + + ]: 11739 : if(center==-3)nsymconv=4*nsym
1474 [ - + + - : 35217 : ABI_MALLOC(symrelconv,(3,3,nsymconv))
+ - + - -
+ - + ]
1475 [ - + + - : 35217 : ABI_MALLOC(tnonsconv,(3,nsymconv))
+ - + - -
+ - + ]
1476 : :
1477 : : !Produce symrel and tnons in conventional axes,
1478 : : !name them symrelconv and tnonsconv
1479 [ + + ]: 46956 : rprimdconv(:,1)=bravais(3:5)
1480 [ + + ]: 46956 : rprimdconv(:,2)=bravais(6:8)
1481 [ + + ]: 46956 : rprimdconv(:,3)=bravais(9:11)
1482 : :
1483 [ + + + + : 79803 : if(center/=0)rprimdconv(:,:)=rprimdconv(:,:)*half
+ + ]
1484 : :
1485 : 11739 : axes(:,:)=zero
1486 : 11739 : axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one
1487 [ + + + + : 3764852 : symrelconv(:,:,1:nsym)=symrel(:,:,1:nsym)
+ + ]
1488 : : !Note that the number of symmetry operations is still nsym
1489 : 11739 : call symrelrot(nsym,rprimdconv,axes,symrelconv,tolsym)
1490 : :
1491 : 11739 : call xred2xcart(nsym,rprimdconv,tnonsconv,tnons)
1492 : : !Gives the associated translation, with components in the interval ]-0.5,0.5] .
1493 [ + + + + ]: 1166543 : tnonsconv(:,1:nsym)=tnonsconv(:,1:nsym)-nint(tnonsconv(:,1:nsym)-tol6)
1494 : :
1495 : : !If the Bravais lattice is centered, duplicate or quadruplicate
1496 : : !the number of symmetry operations, using the Bravais lattice shifts
1497 : 11739 : nshift=1
1498 [ + + ]: 11739 : if(center/=0)nshift=2
1499 [ + + ]: 11739 : if(center==-3)nshift=4
1500 [ - + + - : 23478 : ABI_MALLOC(shift,(3,nshift))
- + - + ]
1501 [ + + ]: 46956 : shift(:,1)=zero
1502 [ + + ]: 11739 : if(center/=0 .and. center/=-3)then
1503 [ + + ]: 5188 : shift(:,2)=half
1504 [ + + ]: 1297 : if(center==1)shift(1,2)=zero
1505 [ + + ]: 1297 : if(center==2)shift(2,2)=zero
1506 [ + + ]: 1297 : if(center==3)shift(3,2)=zero
1507 [ + + ]: 10442 : else if(center==-3)then
1508 [ + + ]: 17500 : shift(:,2)=half ; shift(1,2)=zero
1509 [ + + ]: 17500 : shift(:,3)=half ; shift(2,3)=zero
1510 [ + + ]: 17500 : shift(:,4)=half ; shift(3,4)=zero
1511 : : end if ! center/=0 or -3
1512 [ + + ]: 11739 : if(nshift/=1)then
1513 [ + + ]: 20094 : do ishift=2,nshift
1514 [ + - + + : 14844718 : symrelconv(:,:,(ishift-1)*nsym+1:ishift*nsym)=symrelconv(:,:,1:nsym)
+ + + + +
+ + + +
+ ]
1515 [ + + ]: 590490 : do isym=1,nsym
1516 [ + + ]: 2296006 : tnonsconv(:,(ishift-1)*nsym+isym)=tnonsconv(:,isym)+shift(:,ishift)
1517 : : end do
1518 : : end do ! ishift
1519 : : end if ! nshift/=1
1520 : :
1521 : : !At this stage, all the symmetry operations are available,
1522 : : !expressed in the conventional axis, and also include
1523 : : !the Bravais lattive translations, and associated operations...
1524 : :
1525 : 11739 : n_axes(:)=0
1526 : :
1527 [ - + + - : 35217 : ABI_MALLOC(determinant,(nsymconv))
+ - + - -
+ - + ]
1528 : :
1529 : : !Get the determinant
1530 : 11739 : call symdet(determinant,nsymconv,symrelconv)
1531 : :
1532 : : !Get the order of each the symmetry operation, as well as the maximal order
1533 : : !Also, examine whether each symmetry operation is the inversion, or a root of the inversion (like -3)
1534 : : !Decide which kind of point symmetry operation it is
1535 : : !Finally assign tnonsconv order and decide the space symmetry operation
1536 : :
1537 [ - + + - : 23478 : ABI_MALLOC(t_axes,(nsymconv))
- + - + ]
1538 : :
1539 [ + + ]: 870836 : do isym=1,nsymconv
1540 : :
1541 : : ! Note : nsymconv might be bigger than 192, but only for non-primitive cells, in which case labels will not be echoed anywhere.
1542 : : ! 192 is the fixed dimension of labels, so this avoids possible memory problems.
1543 : : call symcharac(center, determinant(isym), iholohedry, isym, labels(mod(isym-1,192)+1), &
1544 : 859097 : symrelconv(:,:,isym), tnonsconv(:,isym), t_axes(isym))
1545 [ - + ]: 859097 : if (t_axes(isym) == -1) then
1546 : 0 : write(msg, '(a,a,i3,a,3(a,3i4,a),a,3es22.12,a,a,3es22.12)' )ch10,&
1547 : 0 : ' symspgr: problem with isym=',isym,ch10,&
1548 : 0 : ' symrelconv(:,1,isym)=',symrelconv(:,1,isym),ch10,&
1549 : 0 : ' symrelconv(:,2,isym)=',symrelconv(:,2,isym),ch10,&
1550 : 0 : ' symrelconv(:,3,isym)=',symrelconv(:,3,isym),ch10,&
1551 : 0 : ' tnonsconv(:,isym)=',tnonsconv(:,isym)
1552 : 0 : call wrtout(std_out,msg)
1553 : 0 : write(msg, '(a,i0,2a)' )'The space symmetry operation number',isym,ch10,'is not a (translated) root of unity'
1554 : 0 : ABI_BUG(msg)
1555 [ - + ]: 859097 : else if (t_axes(isym) == -2) then
1556 : 0 : write(msg, '(a,i0,a)' )'The symmetry operation number ',isym,' is not a root of unity'
1557 : 0 : ABI_BUG(msg)
1558 : : end if
1559 : :
1560 : 870836 : n_axes(t_axes(isym))=n_axes(t_axes(isym))+1
1561 : : end do ! isym=1,nsymconv
1562 : :
1563 [ + + - + ]: 375648 : if (sum(n_axes)-nsymconv/=0) then
1564 : : write(msg, '(7a)' )&
1565 : 0 : 'Not all the symmetries have been recognized. ',ch10,&
1566 : 0 : 'This might be due either to an error in the input file',ch10,&
1567 : 0 : 'or to a BUG in ABINIT',ch10,&
1568 : 0 : 'Please contact the ABINIT group.'
1569 : 0 : ABI_WARNING(msg)
1570 : : end if
1571 : :
1572 : : !write(std_out,*)' symspgr : brvltt,nsymconv=',brvltt,nsymconv
1573 : : !write(std_out,*)' n_axes(1:10)=',n_axes(1:10)
1574 : : !write(std_out,*)' n_axes(11:20)=',n_axes(11:20)
1575 : : !write(std_out,*)' n_axes(21:31)=',n_axes(21:31)
1576 : :
1577 : : !Treat cases in which the space group cannot be identified on the
1578 : : !basis of n_axes one need additional information
1579 [ + + ]: 11739 : if(brvltt==1)then
1580 : : ! If the bravais lattice is primitive
1581 [ + + ]: 5255 : if(nsymconv==4)then
1582 : 309 : n_axest=(/0,0,0,0,0,0,0,1,1,0, 0,0,0,0,0,2,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0/)
1583 [ + + + + ]: 9888 : if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 27 (Pcc2) or 32 (Pba2)
1584 : 16 : write(std_out,*)' symspgr: 27 or 32'
1585 : 16 : additional_info=2
1586 : : ! Select binary axis
1587 [ + + ]: 80 : do isym=1,nsymconv
1588 [ + + ]: 80 : if(t_axes(isym)==8)then
1589 : : ! Find direction of binary axis
1590 [ + - ]: 16 : if(symrelconv(1,1,isym)==1)direction=1
1591 [ + - ]: 16 : if(symrelconv(2,2,isym)==1)direction=2
1592 [ + - ]: 16 : if(symrelconv(3,3,isym)==1)direction=3
1593 : : end if
1594 : : end do
1595 : : ! Examine the projection of the translation vector of the a, b or c mirror planes onto the binary axis
1596 [ + + ]: 80 : do isym=1,nsymconv
1597 [ + + ]: 80 : if(t_axes(isym)==16)then
1598 [ + + ]: 32 : if(abs(tnonsconv(direction,isym))>tol8)additional_info=1
1599 : : end if
1600 : : end do
1601 : : end if
1602 [ + + ]: 4946 : else if(nsymconv==8)then
1603 : 881 : n_axest=(/0,0,0,0,1,0,0,1,1,0, 0,0,0,0,1,2,0,0,0,2, 0,0,0,0,0,0,0,0,0,0,0/)
1604 [ + + + + ]: 28192 : if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 55 (Pbam) or 57 (Pbcm)
1605 : 26 : write(std_out,*)' symspgr: 55 or 57'
1606 : 26 : additional_info=1
1607 : : ! Select mirror plane m
1608 [ + + ]: 234 : do isym=1,nsymconv
1609 [ + + ]: 234 : if(t_axes(isym)==15)then
1610 : : ! Find direction of mirror plane
1611 [ - + ]: 26 : if(symrelconv(1,1,isym)==-1)direction=1
1612 [ - + ]: 26 : if(symrelconv(2,2,isym)==-1)direction=2
1613 [ + - ]: 26 : if(symrelconv(3,3,isym)==-1)direction=3
1614 : : end if
1615 : : end do
1616 : : ! Examine the projection of the translation vector of the a, b, or c mirror planes
1617 : : ! onto the binary axis
1618 [ + + ]: 234 : do isym=1,nsymconv
1619 [ + + ]: 234 : if(t_axes(isym)==16)then
1620 [ + + ]: 52 : if(abs(tnonsconv(direction,isym))>tol8)additional_info=2
1621 : : end if
1622 : : end do
1623 : : end if
1624 : 881 : n_axest=(/0,0,0,0,1,0,0,1,1,0, 0,0,0,0,0,2,0,1,0,2, 0,0,0,0,0,0,0,0,0,0,0/)
1625 [ + + + + ]: 28192 : if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 56 (Pccn) or 60 (Pbcn)
1626 : 26 : write(std_out,*)' symspgr: 56 or 60'
1627 : 26 : additional_info=1
1628 : : ! Select mirror plane n
1629 [ + + ]: 234 : do isym=1,nsymconv
1630 [ + + ]: 234 : if(t_axes(isym)==18)then
1631 : : ! Find direction of mirror plane
1632 [ - + ]: 26 : if(symrelconv(1,1,isym)==-1)direction=1
1633 [ - + ]: 26 : if(symrelconv(2,2,isym)==-1)direction=2
1634 [ + - ]: 26 : if(symrelconv(3,3,isym)==-1)direction=3
1635 : : end if
1636 : : end do
1637 : : ! Examine the projection of the translation vector of the a, b, or c mirror planes
1638 : : ! onto the binary axis
1639 [ + + ]: 234 : do isym=1,nsymconv
1640 [ + + ]: 234 : if(t_axes(isym)==16)then
1641 [ + + ]: 52 : if(abs(tnonsconv(direction,isym))<tol8)additional_info=2
1642 : : end if
1643 : : end do
1644 : : end if
1645 : : end if
1646 [ + + ]: 6484 : else if(brvltt==2)then
1647 : : ! In the few next lines, use additional_info as a flag
1648 : 698 : additional_info=0
1649 : : ! If the bravais lattice is inner-centered
1650 [ + + ]: 698 : if(nsymconv==8)then
1651 : : ! Test spgroup 23 (I222) or 24 (I2_{1}2_{1}2_{1})
1652 : 105 : n_axest=(/0,0,0,0,0,0,1,1,3,0, 0,0,0,0,0,0,0,0,0,3, 0,0,0,0,0,0,0,0,0,0,0/)
1653 [ + + + + ]: 3360 : if(sum((n_axes-n_axest)**2)==0) additional_info=1
1654 [ + + ]: 593 : else if(nsymconv==24)then
1655 : : ! Test spgroup 197 (I23) or 199 (I2_{1}3)
1656 : 2 : n_axest=(/0,0,0,0,0,0,1,1,3,16, 0,0,0,0,0,0,0,0,0,3, 0,0,0,0,0,0,0,0,0,0,0/)
1657 [ + + + - ]: 64 : if(sum((n_axes-n_axest)**2)==0) additional_info=1
1658 : : end if
1659 [ + + ]: 698 : if(additional_info==1)then
1660 : 8 : write(std_out,*)' symspgr: (23 or 24) or (197 or 199)'
1661 : : ! Select the three binary axes (they might be 2 or 2_1 !)
1662 : 8 : test_direction(:)=0
1663 [ + + ]: 104 : do isym=1,nsymconv
1664 [ + + ]: 104 : if(t_axes(isym)==20)then
1665 : : ! Find direction of axis
1666 [ + + ]: 96 : do direction=1,3
1667 [ + + ]: 96 : if(symrelconv(direction,direction,isym)==1)then
1668 : 24 : test_direction(direction)=1
1669 [ - + ]: 24 : if(abs(tnonsconv(direction,isym))<tol8)then
1670 [ # # ]: 0 : vect(:,direction)=tnonsconv(:,isym)
1671 : : else
1672 [ + + ]: 96 : vect(:,direction)=tnonsconv(:,isym)+half
1673 : : end if
1674 [ + + ]: 96 : vect(:,direction)=vect(:,direction)-nint(vect(:,direction)-tol8)
1675 : 24 : vect(direction,direction)=zero
1676 : : end if
1677 : : end do ! direction=1,3
1678 : : end if ! if binary axis
1679 : : end do ! isym
1680 [ + - - + : 8 : if(test_direction(1)/=1 .or. test_direction(2)/=1 .and. test_direction(3)/=1)then
- - ]
1681 : : write(msg, '(5a,3i4)' )&
1682 : 0 : 'For space groups 23, 24, 197 or 197, the three binary axes',ch10,&
1683 : 0 : 'are not equally partitioned along the x, y and z directions',ch10,&
1684 : 0 : 'test_direction(1:3)=',test_direction(:)
1685 : 0 : ABI_BUG(msg)
1686 : : end if
1687 : 8 : additional_info=1
1688 [ + + + - : 8 : if(abs(vect(1,2)-vect(1,3))>tol8 .or. abs(vect(2,1)-vect(2,3))>tol8 .or. &
- + ]
1689 : 4 : abs(vect(3,1)-vect(3,2))>tol8) additional_info=2
1690 : : end if ! additional information are needed
1691 : : end if ! brvltt==1
1692 : :
1693 [ + + ]: 11739 : if (brvltt==0 .or. brvltt==1) then ! Primitive
1694 : 5255 : call symlist_prim(additional_info,nsymconv,n_axes,spgroup)
1695 [ + + ]: 6484 : else if(brvltt==2)then
1696 : 698 : call symlist_bcc(additional_info,nsymconv,n_axes,spgroup)
1697 [ + + ]: 5786 : else if(brvltt==3)then
1698 : 4375 : call symlist_fcc(nsymconv,n_axes,spgroup)
1699 : : else
1700 : 1411 : call symlist_others(brvltt,nsymconv,n_axes,spgroup)
1701 : : end if
1702 : :
1703 [ + + ]: 11739 : if(spgroup==0) then
1704 : : write(msg, '(5a)' )&
1705 : 17 : 'Could not find the space group.',ch10,&
1706 : 17 : 'This often happens when the user selects a restricted set of symmetries ',ch10,&
1707 : 34 : 'in the input file, instead of letting the code automatically find symmetries.'
1708 : 17 : ABI_WARNING(msg)
1709 : : end if
1710 : :
1711 : 11739 : spgorig=1 ; spgaxor=1
1712 : 11739 : call spgdata(brvsb,intsb,intsbl,ptintsb,ptschsb,schsb,spgaxor,spgroup,sporder,spgorig)
1713 : :
1714 [ + + ]: 11739 : if(spgroup/=0)then
1715 : 11722 : write(msg, '(a,i4,2x,a,a,a,a,a)' ) ' symspgr: spgroup=',spgroup,trim(brvsb),trim(intsb),' (=',trim(schsb),')'
1716 : 11722 : call wrtout(std_out,msg)
1717 : : end if
1718 : :
1719 [ + + ]: 11739 : if(bravais(1)==7)then
1720 : 5622 : write(msg, '(a)' ) ' symspgr: optical characteristics = isotropic '
1721 : 5622 : call wrtout(std_out,msg)
1722 [ + + ]: 6117 : else if(bravais(1)==4 .or. bravais(1)==5 .or. bravais(1)==6)then
1723 : 3060 : write(msg, '(a)' ) ' symspgr: optical characteristics = uniaxial '
1724 : 3060 : call wrtout(std_out,msg)
1725 : : ! Identify the first symmetry operation that is order 3, 4 or 6
1726 : 3060 : found=0
1727 [ + + ]: 16596 : do isym=1,nsym
1728 [ + + + + ]: 232019 : if( minval( abs( t_axes(isym)-(/10,12,14,22,23,24,25,26,27,28,29,30,31/) ))==0) then
1729 : : ! Proper rotations
1730 : : found=1 ; exit
1731 : :
1732 [ + + + + ]: 58124 : else if( minval( abs( t_axes(isym)-(/1,2,3/) ))==0) then
1733 : : ! Improper symmetry operations
1734 : 995 : found=-1 ; exit
1735 : : end if
1736 : : end do
1737 [ + + ]: 3060 : if(found==-1 .or. found==1)then
1738 [ + + + + ]: 39455 : symrel_uni=symrel(:,:,isym)
1739 [ + + + + : 14975 : if(found==-1)symrel_uni=-symrel_uni
+ + ]
1740 : : ! Now, symrel_uni is a rotation of order 3, 4, 6, for which the axis must be identified
1741 : : ! It is actually the only eigenvector with eigenvalue 1. It can be found by cross products
1742 : : ! Subtract the unit matrix.
1743 [ + + ]: 12140 : do ii=1,3
1744 : 12140 : symrel_uni(ii,ii)=symrel_uni(ii,ii)-1
1745 : : end do
1746 : 3650 : found=0
1747 [ - + ]: 3650 : do ii=1,3
1748 [ + + ]: 3650 : jj=ii+1 ; if(jj==4)jj=1
1749 : : ! Cross product
1750 [ + + + + ]: 25550 : ivec1 = symrel_uni(ii,:); ivec2 = symrel_uni(jj,:)
1751 : 3650 : uniaxis = ivec1 .x. ivec2
1752 [ + + + + ]: 14600 : if(sum(uniaxis**2)/=0)then
1753 : : found=1 ; exit
1754 : : end if
1755 : : end do
1756 [ + - ]: 3035 : if(found==1)then
1757 : : ! Try to reduce the length, by an integer factor (try only primes 2, 3, 5, 7, 11)
1758 : 3035 : prime=(/2,3,5,7,11/)
1759 : 3035 : ii=1
1760 [ + + ]: 19780 : do while (ii<6)
1761 [ + + ]: 66980 : uniaxis_try=uniaxis/prime(ii)
1762 [ + + + + ]: 70015 : if(sum(abs(uniaxis_try*prime(ii)-uniaxis))==0)then
1763 : 1570 : uniaxis=uniaxis_try
1764 : : else
1765 : 15175 : ii=ii+1
1766 : : end if
1767 : : end do
1768 : 3035 : write(msg, '(a,3i4)' ) ' Optical axis (in reduced coordinates, real space ) :',uniaxis
1769 : : end if
1770 : : end if
1771 : : if(found==0)then
1772 : 25 : write(msg, '(a)' ) ' However, the axis has not been found. Sorry for this.'
1773 : : end if
1774 : 3060 : call wrtout(std_out,msg)
1775 : : end if
1776 : :
1777 [ - + ]: 11739 : ABI_FREE(determinant)
1778 [ - + ]: 11739 : ABI_FREE(shift)
1779 [ - + ]: 11739 : ABI_FREE(symrelconv)
1780 [ - + ]: 11739 : ABI_FREE(tnonsconv)
1781 [ - + ]: 11739 : ABI_FREE(t_axes)
1782 : :
1783 : : DBG_EXIT("COLL")
1784 : :
1785 : 11739 : end subroutine symspgr
1786 : : !!***
1787 : :
1788 : : !!****f* m_symfind/symlatt
1789 : : !! NAME
1790 : : !! symlatt
1791 : : !!
1792 : : !! FUNCTION
1793 : : !! From the unit cell vectors (rprimd) and the corresponding metric tensor,
1794 : : !! find the Bravais lattice and its symmetry operations (ptsymrel).
1795 : : !! 1) Find the shortest possible primitive vectors for the lattice
1796 : : !! 2) Determines the holohedral group of the lattice, and the
1797 : : !! axes to be used for the conventional cell
1798 : : !! (this is a delicate part, in which the centering of the
1799 : : !! reduced cell must be taken into account)
1800 : : !! The idea is to determine the basis vectors of the conventional
1801 : : !! cell from the reduced cell basis vectors.
1802 : : !! 3) Generate the symmetry operations of the holohedral group
1803 : : !!
1804 : : !! INPUTS
1805 : : !! iout=unit number of output file
1806 : : !! msym=default maximal number of symmetries. WARNING : cannot be simply set to nsym, because
1807 : : !! the number of symmetries found here will likely be bigger than sym !
1808 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1809 : : !! tolsym=tolerance for the symmetries
1810 : : !!
1811 : : !! OUTPUT
1812 : : !! bravais(11): bravais(1)=iholohedry
1813 : : !! bravais(2)=center
1814 : : !! bravais(3:11)=coordinates of rprim in the axes
1815 : : !! of the conventional bravais lattice (*2 if center/=0)
1816 : : !! nptsym=number of point symmetries of the Bravais lattice
1817 : : !! ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
1818 : : !! of the Bravais lattice in real space in terms of primitive translations.
1819 : : !!
1820 : : !! NOTES
1821 : : !! WARNING: bravais(1) might be given a negative value in another
1822 : : !! routine, if the cell is non-primitive.
1823 : : !! The holohedral groups are numbered as follows
1824 : : !! (see international tables for crystallography (1983), p. 13)
1825 : : !! iholohedry=1 triclinic 1bar
1826 : : !! iholohedry=2 monoclinic 2/m
1827 : : !! iholohedry=3 orthorhombic mmm
1828 : : !! iholohedry=4 tetragonal 4/mmm
1829 : : !! iholohedry=5 trigonal 3bar m
1830 : : !! iholohedry=6 hexagonal 6/mmm
1831 : : !! iholohedry=7 cubic m3bar m
1832 : : !! Centering
1833 : : !! center=0 no centering
1834 : : !! center=-1 body-centered
1835 : : !! center=-3 face-centered
1836 : : !! center=1 A-face centered
1837 : : !! center=2 B-face centered
1838 : : !! center=3 C-face centered
1839 : : !!
1840 : : !! SOURCE
1841 : :
1842 : 40278 : subroutine symlatt(bravais,iout,msym,nptsym,ptsymrel,rprimd,tolsym)
1843 : :
1844 : : !Arguments ------------------------------------
1845 : : !scalars
1846 : : integer,intent(in) :: iout,msym
1847 : : integer,intent(out) :: nptsym
1848 : : real(dp),intent(in) :: tolsym
1849 : : !arrays
1850 : : integer,intent(out) :: bravais(11),ptsymrel(3,3,msym)
1851 : : real(dp),intent(in) :: rprimd(3,3)
1852 : :
1853 : : !Local variables-------------------------------
1854 : : !scalars
1855 : : integer,parameter :: mgen=4
1856 : : integer :: center,fact,found,foundc,ia,iaxis1,iaxis2
1857 : : integer :: isign1,isign2,ib,icase,igen,iholohedry,ii,index,isym
1858 : : integer :: itrial,jj,jsym,ngen=0,orthogonal,sign12,sign13,sign23,sumsign
1859 : : real(dp) :: determinant,norm2a,norm2b,norm2c,norm2trial,reduceda,reducedb,sca
1860 : : real(dp) :: scalarprod,scb,trace,trace_best,val
1861 : : character(len=500) :: msg
1862 : : !arrays
1863 : : integer,parameter :: list_holo(7)=(/7,6,4,3,5,2,1/)
1864 : : integer :: ang90(3),equal(3),gen(3,3,mgen),gen2xy(3,3),gen2y(3,3),gen2z(3,3)
1865 : : integer :: gen3(3,3),gen6(3,3),icoord(3,3),identity(3,3),nvecta(3),nvectb(3)
1866 : : integer :: order(mgen)
1867 : : real(dp) :: axes(3,3),axesinvt(3,3),axes_best(3,3),axes_try(3,3)
1868 : : real(dp) :: cell_base(3,3),coord(3,3),metmin(3,3)
1869 : : real(dp) :: minim(3,3),scprods(3,3),vecta(3),vectb(3),vectc(3),vin1(3),vin2(3),vext(3)
1870 : : !**************************************************************************
1871 : :
1872 : : !DEBUG
1873 : : !write(std_out,'(a,es14.6)') ' m_symfind%symlatt : enter, tolsym= ',tolsym
1874 : : !call flush(std_out)
1875 : : !ENDDEBUG
1876 : :
1877 : 40278 : identity(:,:)=0 ; identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
1878 : 40278 : nvecta(1)=2 ; nvectb(1)=3
1879 : 40278 : nvecta(2)=1 ; nvectb(2)=3
1880 : 40278 : nvecta(3)=1 ; nvectb(3)=2
1881 : :
1882 : : !--------------------------------------------------------------------------
1883 : : !Reduce the input vectors to a set of minimal vectors
1884 : 40278 : call smallprim(metmin,minim,rprimd)
1885 : :
1886 : : !DEBUG
1887 : : !write(std_out,*)' symlatt : minim(:,1)=',minim(:,1)
1888 : : !write(std_out,*)' symlatt : minim(:,2)=',minim(:,2)
1889 : : !write(std_out,*)' symlatt : minim(:,3)=',minim(:,3)
1890 : : !call flush(std_out)
1891 : : !ENDDEBUG
1892 : :
1893 : : !--------------------------------------------------------------------------
1894 : : !Examine the angles and vector lengths
1895 : 40278 : ang90(:)=0
1896 [ + + ]: 40278 : if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1
1897 [ + + ]: 40278 : if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1
1898 [ + + ]: 40278 : if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1
1899 : 40278 : equal(:)=0
1900 [ + + ]: 40278 : if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1
1901 [ + + ]: 40278 : if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1
1902 [ + + ]: 40278 : if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1
1903 : :
1904 : : !DEBUG
1905 : : !write(std_out,*)' ang90=',ang90(:)
1906 : : !write(std_out,*)' equal=',equal(:)
1907 : : !call flush(std_out)
1908 : : !ENDDEBUG
1909 : :
1910 : : !-----------------------------------------------------------------------
1911 : : !Identification of the centering
1912 : :
1913 : 40278 : foundc=0
1914 : : !Default values
1915 : 40278 : fact=1 ; center=0
1916 : 40278 : cell_base(:,:)=minim(:,:)
1917 : :
1918 : : !Examine each holohedral group
1919 : : !This search is ordered : should not be happy with tetragonal,
1920 : : !while there is FCC ...
1921 [ + + ]: 126813 : do index=1,6
1922 : :
1923 : : ! If the holohedry is already found, exit
1924 [ + + ]: 123013 : if(foundc==1)exit
1925 : :
1926 : : ! Initialize the target holohedry
1927 : 86535 : iholohedry=list_holo(index)
1928 : :
1929 : : ! DEBUG
1930 : : ! write(std_out,*)
1931 : : ! write(std_out,*)' symlatt : trial holohedry',iholohedry
1932 : : ! ENDDEBUG
1933 : :
1934 : 86535 : orthogonal=0
1935 [ + + + + ]: 86535 : if(iholohedry==7 .or. iholohedry==4 .or. iholohedry==3)orthogonal=1
1936 : :
1937 : : ! Now, will examine different working hypothesis.
1938 : : ! The set of these hypothesis is thought to cover all possible cases ...
1939 : :
1940 : : ! Working hypothesis : the basis is orthogonal
1941 [ + + + + ]: 86535 : if(ang90(1)+ang90(2)+ang90(3)==3 .and. orthogonal==1)then
1942 : 21486 : fact=1 ; center=0
1943 : 21486 : cell_base(:,:)=minim(:,:)
1944 : : ! Checks that the basis vectors are OK for the target holohedry
1945 : 21486 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1946 : : end if
1947 : :
1948 : : ! Select one trial direction
1949 [ + + ]: 303739 : do itrial=1,3
1950 : :
1951 : : ! If the holohedry is already found, exit
1952 [ + + ]: 211869 : if(foundc==1)exit
1953 : :
1954 : 176926 : ia=nvecta(itrial) ; ib=nvectb(itrial)
1955 : :
1956 : : ! This is in case of hexagonal holohedry
1957 [ + - + + : 176926 : if(foundc==0 .and. iholohedry==6 .and. ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then
+ + + + +
+ ]
1958 : 7406 : reduceda=metmin(ib,ia)/metmin(ia,ia)
1959 : 7406 : fact=1 ; center=0
1960 [ + + ]: 7406 : if(abs(reduceda+0.5d0)<tolsym)then
1961 [ + + ]: 8520 : cell_base(:,1)=minim(:,ia)
1962 [ + + ]: 8520 : cell_base(:,2)=minim(:,ib)
1963 [ + + ]: 8520 : cell_base(:,3)=minim(:,itrial)
1964 : : !DEBUG
1965 : : ! write(std_out,*)' cell_base(:,1)=',cell_base(:,1)
1966 : : ! write(std_out,*)' cell_base(:,2)=',cell_base(:,2)
1967 : : ! write(std_out,*)' cell_base(:,3)=',cell_base(:,3)
1968 : : !ENDDEBUG
1969 : : ! Checks that the basis vectors are OK for the target holohedry
1970 : 2130 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1971 [ + + ]: 5276 : else if(abs(reduceda-0.5d0)<tolsym)then
1972 [ + + ]: 276 : cell_base(:,1)=minim(:,ia)
1973 [ + + ]: 276 : cell_base(:,2)=minim(:,ib)-minim(:,ia)
1974 [ + + ]: 276 : cell_base(:,3)=minim(:,itrial)
1975 : : ! Checks that the basis vectors are OK for the target holohedry
1976 : 69 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1977 : : end if
1978 : : end if
1979 : :
1980 : : ! Working hypothesis : the conventional cell is orthogonal,
1981 : : ! and the two other vectors are axes of the conventional cell
1982 [ + + + + : 176926 : if(foundc==0 .and. orthogonal==1 .and. ang90(itrial)==1)then
+ + ]
1983 : :
1984 : : ! Compute the reduced coordinate of trial vector in the basis
1985 : : ! of the two other vectors
1986 : 37612 : reduceda=metmin(itrial,ia)/metmin(ia,ia)
1987 : 37612 : reducedb=metmin(itrial,ib)/metmin(ib,ib)
1988 [ + + ]: 150448 : cell_base(:,ia)=minim(:,ia)
1989 [ + + ]: 150448 : cell_base(:,ib)=minim(:,ib)
1990 [ + + + + : 37612 : if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. &
+ + ]
1991 [ + + ]: 37612 : & ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym) )then
1992 [ + + ]: 4602 : if(abs(abs(reduceda)-0.5d0)<tolsym)center=ib
1993 [ + + ]: 4602 : if(abs(abs(reducedb)-0.5d0)<tolsym)center=ia
1994 : 4602 : fact=2
1995 : : cell_base(:,itrial)= &
1996 [ + + ]: 18408 : & (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0
1997 : 4602 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1998 [ + + + + ]: 33010 : else if( abs(abs(reduceda)-0.5d0)<tolsym .and.&
1999 : : & abs(abs(reducedb)-0.5d0)<tolsym ) then
2000 : 553 : fact=2 ; center=-1
2001 : : cell_base(:,itrial)= &
2002 [ + + ]: 2212 : & (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0
2003 : 553 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2004 : : end if
2005 : : end if
2006 : :
2007 : : ! Working hypothesis : the conventional cell is orthogonal, and
2008 : : ! the trial vector is one of the future axes, and the face perpendicular to it is centered
2009 : : if(foundc==0 .and. iholohedry==3 .and. &
2010 [ + + + + : 176926 : & ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then
+ + + + +
+ ]
2011 : 460 : fact=2 ; center=itrial
2012 [ + + ]: 1840 : cell_base(:,ia)=minim(:,ia)+minim(:,ib)
2013 [ + + ]: 1840 : cell_base(:,ib)=minim(:,ia)-minim(:,ib)
2014 [ + + ]: 1840 : cell_base(:,itrial)=minim(:,itrial)
2015 : : ! Checks that the basis vectors are OK for the target holohedry
2016 : 460 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2017 : : end if
2018 : :
2019 : : ! DEBUG
2020 : : ! write(std_out,*)' after test_b, foundc=',foundc
2021 : : ! ENDDEBUG
2022 : :
2023 : : ! Working hypothesis : the conventional cell is orthogonal, and
2024 : : ! the trial vector is one of the future axes
2025 [ + + + + ]: 176926 : if(foundc==0 .and. orthogonal==1)then
2026 : : ! Compute the projection of the two other vectors on the trial vector
2027 : 107035 : reduceda=metmin(itrial,ia)/metmin(itrial,itrial)
2028 : 107035 : reducedb=metmin(itrial,ib)/metmin(itrial,itrial)
2029 : : ! If both projections are half-integer, one might have found an axis
2030 [ + + + + ]: 107035 : if( abs(abs(reduceda)-0.5d0)<tolsym .and.&
2031 : : & abs(abs(reducedb)-0.5d0)<tolsym ) then
2032 [ + + ]: 96608 : vecta(:)=minim(:,ia)-reduceda*minim(:,itrial)
2033 [ + + ]: 96608 : vectb(:)=minim(:,ib)-reducedb*minim(:,itrial)
2034 : 24152 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2035 : 24152 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2036 : 24152 : scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
2037 : : ! Note the order of selection : body-centered is prefered
2038 : : ! over face centered, which is correct for the tetragonal case
2039 [ + + ]: 24152 : if(abs(norm2a-norm2b)<tolsym*half*(norm2a+norm2b))then
2040 : : ! The lattice is body centered
2041 : 17279 : fact=2 ; center=-1
2042 [ + + ]: 69116 : cell_base(:,ia)=vecta(:)+vectb(:)
2043 [ + + ]: 69116 : cell_base(:,ib)=vecta(:)-vectb(:)
2044 [ + + ]: 69116 : cell_base(:,itrial)=minim(:,itrial)
2045 : 17279 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2046 [ + - ]: 6873 : else if(abs(scalarprod)<tolsym*half*(norm2a+norm2b))then
2047 : : ! The lattice is face centered
2048 : 0 : fact=2 ; center=-3
2049 [ # # ]: 0 : cell_base(:,ia)=2.0d0*vecta(:)
2050 [ # # ]: 0 : cell_base(:,ib)=2.0d0*vectb(:)
2051 [ # # ]: 0 : cell_base(:,itrial)=minim(:,itrial)
2052 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2053 : : end if
2054 : : end if
2055 : : end if
2056 : :
2057 : : ! DEBUG
2058 : : ! write(std_out,*)' after test_c, foundc=',foundc
2059 : : ! ENDDEBUG
2060 : :
2061 : : ! Working hypothesis : the conventional cell is orthogonal,
2062 : : ! and body centered with no basis vector being an axis,
2063 : : ! in which case the basis vectors must be equal (even for orthorhombic)
2064 : : if(foundc==0 .and. orthogonal==1 .and. &
2065 [ + + + + : 176926 : & equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then
+ + + + +
- ]
2066 : : ! Compute the combination of the two other vectors
2067 [ + + ]: 93436 : vecta(:)=minim(:,ia)+minim(:,ib)
2068 [ + + ]: 93436 : vectb(:)=minim(:,ia)-minim(:,ib)
2069 : 23359 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2070 : 23359 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2071 : : ! Project the trial vector on the first of the two vectors
2072 : : reduceda=( minim(1,itrial)*vecta(1)+ &
2073 : : & minim(2,itrial)*vecta(2)+ &
2074 : 23359 : & minim(3,itrial)*vecta(3) )/norm2a
2075 : : reducedb=( minim(1,itrial)*vectb(1)+ &
2076 : : & minim(2,itrial)*vectb(2)+ &
2077 : 23359 : & minim(3,itrial)*vectb(3) )/norm2b
2078 [ + + ]: 23359 : if( abs(abs(reduceda)-0.5d0)<tolsym )then
2079 : : ! The first vector is an axis
2080 : 3653 : fact=2 ; center=-1
2081 [ + + ]: 14612 : cell_base(:,ia)=vecta(:)
2082 [ + + ]: 14612 : vecta(:)=minim(:,itrial)-reduceda*vecta(:)
2083 [ + + ]: 14612 : vectb(:)=0.5d0*vectb(:)
2084 [ + + ]: 14612 : cell_base(:,ib)=vecta(:)+vectb(:)
2085 [ + + ]: 14612 : cell_base(:,itrial)=vecta(:)-vectb(:)
2086 : 3653 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2087 [ + + ]: 19706 : else if( abs(abs(reducedb)-0.5d0)<tolsym )then
2088 : : ! The second vector is an axis
2089 : 84 : fact=2 ; center=-1
2090 [ + + ]: 336 : cell_base(:,ib)=vectb(:)
2091 [ + + ]: 336 : vectb(:)=minim(:,itrial)-reducedb*vectb(:)
2092 [ + + ]: 336 : vecta(:)=0.5d0*vecta(:)
2093 [ + + ]: 336 : cell_base(:,ia)=vectb(:)+vecta(:)
2094 [ + + ]: 336 : cell_base(:,itrial)=vectb(:)-vecta(:)
2095 : 84 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2096 : : end if
2097 : : end if
2098 : :
2099 : : ! Working hypothesis : the conventional cell is orthogonal,
2100 : : ! and face centered, in the case where two minimal vectors are equal
2101 [ + + + + : 176926 : if(foundc==0 .and. orthogonal==1 .and. equal(itrial)==1 ) then
+ + ]
2102 : : ! Compute the combination of these two vectors
2103 [ + + ]: 157048 : vecta(:)=minim(:,ia)+minim(:,ib)
2104 [ + + ]: 157048 : vectb(:)=minim(:,ia)-minim(:,ib)
2105 : 39262 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2106 : 39262 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2107 : : ! Project the trial vector on the two vectors
2108 : : reduceda=( minim(1,itrial)*vecta(1)+ &
2109 : : & minim(2,itrial)*vecta(2)+ &
2110 : 39262 : & minim(3,itrial)*vecta(3) )/norm2a
2111 : : reducedb=( minim(1,itrial)*vectb(1)+ &
2112 : : & minim(2,itrial)*vectb(2)+ &
2113 : 39262 : & minim(3,itrial)*vectb(3) )/norm2b
2114 [ + + + + : 39262 : if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. &
+ + ]
2115 [ + + ]: 39262 : & ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym) )then
2116 : 978 : fact=2 ; center=-3
2117 : : cell_base(:,itrial)= &
2118 [ + + ]: 3912 : & (minim(:,itrial)-reduceda*vecta(:)-reducedb*vectb(:) )*2.0d0
2119 [ + + ]: 3912 : cell_base(:,ia)=vecta(:)
2120 [ + + ]: 3912 : cell_base(:,ib)=vectb(:)
2121 : 978 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2122 : : end if
2123 : : end if
2124 : :
2125 : : ! Working hypothesis : the conventional cell is orthogonal,
2126 : : ! face centered, but no two vectors are on the same "square"
2127 [ + + + + ]: 176926 : if(foundc==0 .and. orthogonal==1)then
2128 : : ! Compute the combination of these two vectors
2129 [ + + ]: 419860 : vecta(:)=minim(:,ia)+minim(:,ib)
2130 [ + + ]: 419860 : vectb(:)=minim(:,ia)-minim(:,ib)
2131 : 104965 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2132 : 104965 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2133 : : ! The trial vector length must be equal to one of these lengths
2134 [ + + ]: 104965 : if(abs(metmin(itrial,itrial)-norm2a)<tolsym*norm2a)then
2135 : 1211 : fact=2 ; center=-3
2136 [ + + ]: 4844 : cell_base(:,ia)=vecta(:)+minim(:,itrial)
2137 [ + + ]: 4844 : cell_base(:,ib)=vecta(:)-minim(:,itrial)
2138 : : ! Project vectb perpendicular to cell_base(:,ia) and cell_base(:,ib)
2139 : 1211 : norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2
2140 : 1211 : norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2
2141 : : reduceda=( cell_base(1,ia)*vectb(1)+ &
2142 : : & cell_base(2,ia)*vectb(2)+ &
2143 : 1211 : & cell_base(3,ia)*vectb(3) )/norm2a
2144 : : reducedb=( cell_base(1,ib)*vectb(1)+ &
2145 : : & cell_base(2,ib)*vectb(2)+ &
2146 : 1211 : & cell_base(3,ib)*vectb(3) )/norm2b
2147 [ + + + - ]: 1211 : if( abs(abs(reduceda)-0.5d0)<tolsym .and. &
2148 : : & abs(abs(reducedb)-0.5d0)<tolsym )then
2149 [ # # ]: 0 : cell_base(:,itrial)=vectb(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib)
2150 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2151 : : end if
2152 [ + + ]: 103754 : else if(abs(metmin(itrial,itrial)-norm2b)<tolsym*norm2b)then
2153 : 17405 : fact=2 ; center=-3
2154 [ + + ]: 69620 : cell_base(:,ia)=vectb(:)+minim(:,itrial)
2155 [ + + ]: 69620 : cell_base(:,ib)=vectb(:)-minim(:,itrial)
2156 : : ! Project vecta perpendicular to cell_base(:,ia) and cell_base(:,ib)
2157 : 17405 : norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2
2158 : 17405 : norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2
2159 : : reduceda=( cell_base(1,ia)*vecta(1)+ &
2160 : : & cell_base(2,ia)*vecta(2)+ &
2161 : 17405 : & cell_base(3,ia)*vecta(3) )/norm2a
2162 : : reducedb=( cell_base(1,ib)*vecta(1)+ &
2163 : : & cell_base(2,ib)*vecta(2)+ &
2164 : 17405 : & cell_base(3,ib)*vecta(3) )/norm2b
2165 [ + + + + ]: 17405 : if( abs(abs(reduceda)-0.5d0)<tolsym .and. &
2166 : : & abs(abs(reducedb)-0.5d0)<tolsym )then
2167 [ + + ]: 66572 : cell_base(:,itrial)=vecta(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib)
2168 : 16643 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2169 : : end if
2170 : : end if
2171 : : end if
2172 : :
2173 : : ! Working hypothesis : the cell is rhombohedral, and
2174 : : ! the three minimal vectors have same length and same absolute scalar product
2175 : : if(foundc==0 .and. iholohedry==5 .and. &
2176 [ + + + + : 176926 : & equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then
+ + + + +
- ]
2177 [ + + ]: 347 : if( abs(abs(metmin(1,2))-abs(metmin(1,3)))<tolsym*metmin(1,1) .and. &
2178 [ + - ]: 347 : & abs(abs(metmin(1,2))-abs(metmin(2,3)))<tolsym*metmin(2,2) )then
2179 : 323 : fact=1 ; center=0
2180 : 323 : cell_base(:,:)=minim(:,:)
2181 : : ! One might have to change the sign of one of the vectors
2182 : 323 : sign12=1 ; sign13=1 ; sign23=1
2183 [ + + ]: 323 : if(metmin(1,2)<0.0d0)sign12=-1
2184 [ + + ]: 323 : if(metmin(1,3)<0.0d0)sign13=-1
2185 [ + + ]: 323 : if(metmin(2,3)<0.0d0)sign23=-1
2186 : 323 : sumsign=sign12+sign13+sign23
2187 [ + + ]: 323 : if(sumsign==-1)then
2188 [ - + - - ]: 50 : if(sign12==1)cell_base(:,3)=-cell_base(:,3)
2189 [ - + - - ]: 50 : if(sign13==1)cell_base(:,2)=-cell_base(:,2)
2190 [ + - + + ]: 200 : if(sign23==1)cell_base(:,1)=-cell_base(:,1)
2191 [ + + ]: 273 : else if(sumsign==1)then
2192 [ - + - - ]: 4 : if(sign12==-1)cell_base(:,3)=-cell_base(:,3)
2193 [ - + - - ]: 4 : if(sign13==-1)cell_base(:,2)=-cell_base(:,2)
2194 [ + - + + ]: 16 : if(sign23==-1)cell_base(:,1)=-cell_base(:,1)
2195 : : end if
2196 : 323 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2197 : : end if
2198 : : end if
2199 : :
2200 : : ! DEBUG
2201 : : ! write(std_out,*)' after test_3a, foundc=',foundc
2202 : : ! write(std_out,*)' after test_3a, itrial=',itrial
2203 : : ! write(std_out,*)' after test_3a, equal(:)=',equal(:)
2204 : : ! ENDDEBUG
2205 : :
2206 : : ! Working hypothesis : the cell is rhombohedral, one vector
2207 : : ! is parallel to the trigonal axis
2208 [ + + + + : 176926 : if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 )then
+ + ]
2209 [ + + + + ]: 20888 : vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib)
2210 : 2984 : norm2trial=minim(1,itrial)**2+minim(2,itrial)**2+minim(3,itrial)**2
2211 : : reduceda=( minim(1,itrial)*vecta(1)+ &
2212 : : & minim(2,itrial)*vecta(2)+ &
2213 : 2984 : & minim(3,itrial)*vecta(3) )/norm2trial
2214 : : reducedb=( minim(1,itrial)*vectb(1)+ &
2215 : : & minim(2,itrial)*vectb(2)+ &
2216 : 2984 : & minim(3,itrial)*vectb(3) )/norm2trial
2217 : : ! DEBUG
2218 : : ! write(std_out,*)' reduceda,reducedb=',reduceda,reducedb
2219 : : ! ENDDEBUG
2220 [ + + - + ]: 2984 : if(abs(abs(reduceda)-1.0d0/3.0d0)<tolsym .and. &
2221 : : & abs(abs(reducedb)-1.0d0/3.0d0)<tolsym ) then
2222 : : ! Possibly change of sign to make positive the scalar product with
2223 : : ! the vector parallel to the trigonal axis
2224 [ + + + + ]: 82 : if(reduceda<zero)vecta(:)=-vecta(:)
2225 [ + + + + ]: 130 : if(reducedb<zero)vectb(:)=-vectb(:)
2226 : : ! Projection on the orthogonal plane
2227 [ + + ]: 148 : vecta(:)=vecta(:)-abs(reduceda)*cell_base(:,itrial)
2228 [ + + ]: 148 : vectb(:)=vectb(:)-abs(reducedb)*cell_base(:,itrial)
2229 : : ! These two vectors should have an angle of 120 degrees
2230 : 37 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2231 : 37 : scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
2232 : : ! DEBUG
2233 : : ! write(std_out,*)' norm2a,scalarprod=',norm2a,scalarprod
2234 : : ! ENDDEBUG
2235 [ + + ]: 37 : if(abs(two*scalarprod+norm2a)<tolsym*norm2a)then
2236 : 12 : fact=1 ; center=0
2237 [ - + - - ]: 12 : if(scalarprod>0.0d0)vectb(:)=-vectb(:)
2238 : : ! Now vecta and vectb have an angle of 120 degrees
2239 [ + + ]: 48 : cell_base(:,1)=cell_base(:,itrial)/3.0d0+vecta(:)
2240 [ + + ]: 48 : cell_base(:,2)=cell_base(:,itrial)/3.0d0+vectb(:)
2241 [ + + ]: 48 : cell_base(:,3)=cell_base(:,itrial)/3.0d0-vecta(:)-vectb(:)
2242 : : ! DEBUG
2243 : : ! write(std_out,*)' cell_base(:,1)=',cell_base(:,1)
2244 : : ! write(std_out,*)' cell_base(:,2)=',cell_base(:,2)
2245 : : ! write(std_out,*)' cell_base(:,3)=',cell_base(:,3)
2246 : : ! ENDDEBUG
2247 : 12 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2248 : : end if
2249 : : end if
2250 : : end if
2251 : :
2252 : : ! Working hypothesis : the cell is rhombohedral, one vector
2253 : : ! is in the plane perpendicular to the trigonal axis
2254 [ + + + + : 176926 : if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then
+ + ]
2255 [ + + ]: 11888 : vecta(:)=minim(:,ia)+minim(:,ib)
2256 [ + + ]: 11888 : vectb(:)=minim(:,ia)-minim(:,ib)
2257 : 2972 : norm2trial=cell_base(1,itrial)**2+cell_base(2,itrial)**2+cell_base(3,itrial)**2
2258 : 2972 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2259 : 2972 : norm2b=vecta(1)**2+vecta(2)**2+vecta(3)**2
2260 : : reduceda=( cell_base(1,itrial)*vecta(1)+ &
2261 : : & cell_base(2,itrial)*vecta(2)+ &
2262 : 2972 : & cell_base(3,itrial)*vecta(3) )/norm2trial
2263 : : reducedb=( cell_base(1,itrial)*vectb(1)+ &
2264 : : & cell_base(2,itrial)*vectb(2)+ &
2265 : 2972 : & cell_base(3,itrial)*vectb(3) )/norm2trial
2266 [ + + - + ]: 2972 : if(abs(norm2trial-norm2a)<tolsym*norm2a .and. &
2267 : : & abs(abs(2*reduceda)-norm2trial)<tolsym*norm2trial )then
2268 : 0 : fact=1 ; center=0
2269 [ # # ]: 0 : cell_base(:,1)=minim(:,ia)
2270 [ # # ]: 0 : cell_base(:,2)=-minim(:,ib)
2271 [ # # ]: 0 : cell_base(:,3)=-minim(:,ib)+2*reduceda*minim(:,itrial)
2272 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2273 [ + + - + ]: 2972 : else if (abs(norm2trial-norm2b)<tolsym*norm2b .and. &
2274 : : & abs(abs(2*reducedb)-norm2trial)<tolsym*norm2trial )then
2275 : 0 : fact=1 ; center=0
2276 [ # # ]: 0 : cell_base(:,1)=minim(:,ia)
2277 [ # # ]: 0 : cell_base(:,2)=minim(:,ib)
2278 [ # # ]: 0 : cell_base(:,3)=minim(:,ib)+2*reducedb*minim(:,itrial)
2279 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2280 : : end if
2281 : : end if
2282 : :
2283 : : ! Working hypothesis : the cell is rhombohedral, two vectors
2284 : : ! are in the plane perpendicular to the trigonal axis
2285 [ + + + + : 176926 : if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then
+ + ]
2286 [ + + + + ]: 20804 : vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib)
2287 : 2972 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2288 : 2972 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2289 : 2972 : scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
2290 [ + + ]: 2972 : if(abs(abs(2*scalarprod)-norm2a)<tolsym*norm2a)then
2291 : : ! This is in order to have 120 angle between vecta and vectb
2292 [ + + + + ]: 1957 : if(scalarprod>0.0d0)vectb(:)=-vectb(:)
2293 : : reduceda=( cell_base(1,itrial)*vecta(1)+ &
2294 : : & cell_base(2,itrial)*vecta(2)+ &
2295 : 1246 : & cell_base(3,itrial)*vecta(3) )/norm2a
2296 : : reducedb=( cell_base(1,itrial)*vectb(1)+ &
2297 : : & cell_base(2,itrial)*vectb(2)+ &
2298 : 1246 : & cell_base(3,itrial)*vectb(3) )/norm2b
2299 : 1246 : fact=1 ; center=0
2300 [ + + ]: 4984 : cell_base(:,1)=minim(:,itrial)
2301 [ + + + + ]: 1246 : if(abs(reduceda-0.5d0)<tolsym .and. abs(reducedb)<tolsym )then
2302 [ + + ]: 900 : cell_base(:,2)=minim(:,itrial)-vecta(:)
2303 [ + + ]: 900 : cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:)
2304 : 225 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2305 [ + + + - ]: 1021 : else if(abs(reduceda-0.5d0)<tolsym.and. abs(reducedb+0.5d0)<tolsym )then
2306 [ + + ]: 24 : cell_base(:,2)=minim(:,itrial)-vecta(:)
2307 [ + + ]: 24 : cell_base(:,3)=minim(:,itrial)+vectb(:)
2308 : 6 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2309 [ - + - - ]: 1015 : else if(abs(reduceda)<tolsym .and. abs(reducedb+0.5d0)<tolsym )then
2310 [ # # ]: 0 : cell_base(:,2)=minim(:,itrial)+vectb(:)
2311 [ # # ]: 0 : cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:)
2312 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2313 [ + - + + ]: 1015 : else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb)<tolsym)then
2314 [ + + ]: 24 : cell_base(:,2)=minim(:,itrial)+vecta(:)
2315 [ + + ]: 24 : cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:)
2316 : 6 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2317 [ + - + - ]: 1009 : else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb-0.5d0)<tolsym)then
2318 [ + + ]: 4036 : cell_base(:,2)=minim(:,itrial)+vecta(:)
2319 [ + + ]: 4036 : cell_base(:,3)=minim(:,itrial)-vectb(:)
2320 : 1009 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2321 [ # # # # ]: 0 : else if(abs(reduceda)<tolsym .and. abs(reducedb-0.5d0)<tolsym )then
2322 [ # # ]: 0 : cell_base(:,2)=minim(:,itrial)-vectb(:)
2323 [ # # ]: 0 : cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:)
2324 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2325 : : end if
2326 : : end if
2327 : : end if
2328 : :
2329 : : ! Working hypothesis : monoclinic holohedry, primitive. Then, two angles are 90 degrees
2330 : : if(foundc==0 .and. iholohedry==2 .and. &
2331 [ + + + + : 176926 : & ang90(ia)==1 .and. ang90(ib)==1 ) then
+ + + + ]
2332 : 297 : fact=1 ; center=0
2333 [ + + ]: 1188 : cell_base(:,1)=minim(:,ia)
2334 [ + + ]: 1188 : cell_base(:,2)=minim(:,itrial)
2335 [ + + ]: 1188 : cell_base(:,3)=minim(:,ib)
2336 : : ! Checks that the basis vectors are OK for the target holohedry
2337 : 297 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2338 : : end if
2339 : :
2340 : : ! Monoclinic holohedry, one-face-centered cell
2341 : : ! Working hypothesis, two vectors have equal length.
2342 [ + + ]: 1061556 : do icase=1,5
2343 [ + + + + : 1061556 : if(foundc==0 .and. iholohedry==2 .and. equal(itrial)==1 ) then
+ + ]
2344 [ + + ]: 15492 : vecta(:)=cell_base(:,ia)+cell_base(:,ib)
2345 [ + + ]: 15492 : vectb(:)=cell_base(:,ia)-cell_base(:,ib)
2346 : 3873 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2347 : 3873 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2348 : : ! The minim(:,trial) vector belongs to the
2349 : : ! plane parallel to the cell_base(:,ia),cell_base(:,ib) plane
2350 : : ! In that plane, must try minim(:,itrial),
2351 : : ! as well as the 4 different combinations of
2352 : : ! minim(:,itrial) with the vectors in the plane
2353 [ + + + + ]: 7170 : if(icase==1)vectc(:)=minim(:,itrial)
2354 [ + + + + ]: 5964 : if(icase==2)vectc(:)=minim(:,itrial)+cell_base(:,ia)
2355 [ + + + + ]: 5958 : if(icase==3)vectc(:)=minim(:,itrial)+cell_base(:,ib)
2356 [ + + + + ]: 5946 : if(icase==4)vectc(:)=minim(:,itrial)-cell_base(:,ia)
2357 [ + + + + ]: 5946 : if(icase==5)vectc(:)=minim(:,itrial)-cell_base(:,ib)
2358 : 3873 : norm2c=vectc(1)**2+vectc(2)**2+vectc(3)**2
2359 : : sca=vectc(1)*vecta(1)+&
2360 : : & vectc(2)*vecta(2)+&
2361 : 3873 : & vectc(3)*vecta(3)
2362 : : scb=vectc(1)*vectb(1)+&
2363 : : & vectc(2)*vectb(2)+&
2364 : 3873 : & vectc(3)*vectb(3)
2365 : : ! DEBUG
2366 : : ! write(std_out,*)' symlatt : test iholohedry=2, sca,scb=',sca,scb
2367 : : ! ENDDEBUG
2368 [ + + + + ]: 3873 : if(abs(sca)<tolsym*sqrt(norm2c*norm2a) .or. abs(scb)<tolsym*sqrt(norm2c*norm2b))then
2369 : 408 : fact=2 ; center=3
2370 : : ! The itrial direction is centered
2371 [ + + ]: 1632 : cell_base(:,3)=vectc(:)
2372 [ + + ]: 408 : if(abs(sca)<tolsym*sqrt(norm2c*norm2a))then
2373 [ + + ]: 560 : cell_base(:,2)=vecta(:)
2374 [ + + ]: 560 : cell_base(:,1)=vectb(:)
2375 : 140 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2376 [ - + ]: 268 : else if(abs(scb)<tolsym*sqrt(norm2c*norm2b))then
2377 [ + + ]: 1072 : cell_base(:,2)=vectb(:)
2378 [ + + ]: 1072 : cell_base(:,1)=vecta(:)
2379 : 268 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2380 : : end if
2381 : : end if
2382 : : end if
2383 : : end do ! icase=1,5
2384 : :
2385 : : ! Monoclinic holohedry, one-face-centered cell, but non equivalent.
2386 : : ! This case, one pair of vectors is orthogonal
2387 [ + + + + : 176926 : if(foundc==0 .and. iholohedry==2 .and. ang90(itrial)==1) then
+ + ]
2388 [ + + ]: 4704 : vecta(:)=minim(:,ia)
2389 [ + + ]: 4704 : vectb(:)=minim(:,ib)
2390 : 1176 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2391 : 1176 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2392 : : ! Project the trial vector on the two vectors
2393 : : reduceda=( minim(1,itrial)*vecta(1)+ &
2394 : : & minim(2,itrial)*vecta(2)+ &
2395 : 1176 : & minim(3,itrial)*vecta(3) )/norm2a
2396 : : reducedb=( minim(1,itrial)*vectb(1)+ &
2397 : : & minim(2,itrial)*vectb(2)+ &
2398 : 1176 : & minim(3,itrial)*vectb(3) )/norm2b
2399 [ + + + + ]: 1176 : if(abs(abs(reduceda)-0.5d0)<tolsym .or. abs(abs(reducedb)-0.5d0)<tolsym) then
2400 : 759 : fact=2 ; center=3
2401 [ + + ]: 759 : if(abs(abs(reduceda)-0.5d0)<tolsym)then
2402 [ + + ]: 20 : cell_base(:,2)=vecta(:)
2403 [ + + ]: 20 : cell_base(:,3)=vectb(:)
2404 [ + + ]: 20 : cell_base(:,1)=2*(minim(:,itrial)-reduceda*vecta(:))
2405 : 5 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2406 [ - + ]: 754 : else if(abs(abs(reducedb)-0.5d0)<tolsym)then
2407 [ + + ]: 3016 : cell_base(:,2)=vectb(:)
2408 [ + + ]: 3016 : cell_base(:,3)=vecta(:)
2409 [ + + ]: 3016 : cell_base(:,1)=2*(minim(:,itrial)-reducedb*vectb(:))
2410 : 754 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2411 : : end if
2412 : : end if
2413 : : end if
2414 : :
2415 : : ! Monoclinic holohedry, one-face-centered cell, but non equivalent.
2416 : : ! This case, no pair of vectors is orthogonal, no pair of vector of equal lentgh
2417 [ + + + + ]: 263461 : if(foundc==0 .and. iholohedry==2)then
2418 : : ! Try to find a vector that belongs to the mediator plane, or the binary vector.
2419 : : ! There must be one such vector, if centered monoclinic and no pair of vectors of equal length,
2420 : : ! either among the three vectors, or among one of their differences or sums.
2421 : : ! And there must be, among the two other vectors, one vector whose projection
2422 : : ! on this vector is half the length of this vector.
2423 [ + + ]: 29696 : vecta(:)=minim(:,ia)
2424 [ + + ]: 29696 : vectb(:)=minim(:,ib)
2425 : : ! Try the different possibilities for the vector on which the projection will be half ...
2426 [ + + ]: 32759 : do ii=1,5
2427 [ + + + + ]: 50388 : if(ii==1)vectc(:)=minim(:,itrial)
2428 [ + + + + ]: 48129 : if(ii==2)vectc(:)=minim(:,itrial)+vecta(:)
2429 [ + + + + ]: 42174 : if(ii==3)vectc(:)=minim(:,itrial)-vecta(:)
2430 [ + + + + ]: 42174 : if(ii==4)vectc(:)=minim(:,itrial)+vectb(:)
2431 [ + + + + ]: 42063 : if(ii==5)vectc(:)=minim(:,itrial)-vectb(:)
2432 : 28116 : norm2trial=vectc(1)**2+vectc(2)**2+vectc(3)**2
2433 : : ! Project the two vectors on the trial vector
2434 : 28116 : reduceda=( vectc(1)*vecta(1)+vectc(2)*vecta(2)+vectc(3)*vecta(3) )/norm2trial
2435 : 28116 : reducedb=( vectc(1)*vectb(1)+vectc(2)*vectb(2)+vectc(3)*vectb(3) )/norm2trial
2436 : 28116 : found=0
2437 [ + + ]: 28116 : if(abs(abs(reduceda)-0.5d0)<tolsym)then
2438 : 2055 : vin1(:)=vectc(:)
2439 [ + + ]: 8220 : vin2(:)=2.0d0*(vecta(:)-reduceda*vectc(:))
2440 : 2055 : vext(:)=vectb(:)
2441 : : found=1
2442 [ + + ]: 26061 : else if(abs(abs(reducedb)-0.5d0)<tolsym)then
2443 : 726 : vin1(:)=vectc(:)
2444 [ + + ]: 2904 : vin2(:)=2.0d0*(vectb(:)-reduceda*vectc(:))
2445 : 726 : vext(:)=vecta(:)
2446 : : found=1
2447 : : end if
2448 : 4643 : if(found==1)exit
2449 : : end do
2450 : : ! Now, vin1 and vin2 are perpendicular to each other, and in the plane that contains the binary vector.
2451 : : ! One of them must be the binary vector if any.
2452 : : ! On the other hand, vext is out-of-plane. Might belong to the mediator plane or not.
2453 : : ! If C monoclinc, then the projection of this vext on the binary vector will be either 0 or +1/2 or -1/2.
2454 : : ! The binary axis must be stored in cell_base(:,2) for conventional C-cell
2455 [ + + ]: 7424 : if(found==1)then
2456 : 2781 : found=0
2457 : :
2458 : : ! Test vin1 being the binary axis
2459 : 2781 : norm2trial=vin1(1)**2+vin1(2)**2+vin1(3)**2
2460 : 2781 : reduceda=(vext(1)*vin1(1)+vext(2)*vin1(2)+vext(3)*vin1(3))/norm2trial
2461 [ + + ]: 2781 : if(abs(reduceda)<tolsym)then ! vin1 is the binary axis and vext is in the mediator plane
2462 : 5060 : found=1
2463 [ + + ]: 5060 : cell_base(:,1)=vin2(:)
2464 [ + + ]: 5060 : cell_base(:,2)=vin1(:)
2465 [ + + ]: 5060 : cell_base(:,3)=vext(:)
2466 [ + + ]: 1516 : else if(abs(abs(reduceda)-0.5d0)<tolsym)then ! vin1 is the binary axis and vext has +1/2 or -1/2 as projection
2467 : 144 : found=1
2468 [ + + ]: 144 : cell_base(:,1)=vin2(:)
2469 [ + + ]: 144 : cell_base(:,2)=vin1(:)
2470 [ + + ]: 144 : cell_base(:,3)=vext(:)-reduceda*vin1(:)+vin2(:)*half
2471 : : else
2472 : : ! Test vin2 being the binary axis
2473 : 1480 : norm2trial=vin2(1)**2+vin2(2)**2+vin2(3)**2
2474 : 1480 : reduceda=(vext(1)*vin2(1)+vext(2)*vin2(2)+vext(3)*vin2(3))/norm2trial
2475 [ + + ]: 1480 : if(abs(reduceda)<tolsym)then ! vin2 is the binary axis and vext is in the mediator plane
2476 : 328 : found=1
2477 [ + + ]: 328 : cell_base(:,1)=vin1(:)
2478 [ + + ]: 328 : cell_base(:,2)=vin2(:)
2479 [ + + ]: 328 : cell_base(:,3)=vext(:)
2480 [ - + ]: 1398 : else if(abs(abs(reduceda)-0.5d0)<tolsym)then ! vin2 is the binary axis and vext has +1/2 or -1/2 as projection
2481 : 0 : found=1
2482 [ # # ]: 0 : cell_base(:,1)=vin1(:)
2483 [ # # ]: 0 : cell_base(:,2)=vin2(:)
2484 [ # # ]: 0 : cell_base(:,3)=vext(:)-reduceda*vin2(:)+vin1(:)*half
2485 : : end if
2486 : : end if
2487 : :
2488 : : if(found==1)then
2489 : 1383 : fact=2 ; center=3
2490 : 1383 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2491 : : end if
2492 : : end if
2493 : : end if
2494 : :
2495 : : end do ! Do-loop on three different directions
2496 : : end do ! Do-loop on different target holohedries
2497 : :
2498 [ + + ]: 40278 : if(foundc==0)then
2499 : 1649 : iholohedry=1 ; fact=1 ; center=0
2500 : 1649 : cell_base(:,:)=minim(:,:)
2501 : : end if
2502 : :
2503 : : !DEBUG
2504 : : !write(std_out,*)' symlatt : done with centering tests, foundc=',foundc
2505 : : !write(std_out,*)' center=',center
2506 : : !write(std_out,*)' iholohedry=',iholohedry
2507 : : !call flush(std_out)
2508 : : !ENDDEBUG
2509 : :
2510 : : !--------------------------------------------------------------------------
2511 : : !Final check on the Bravais lattice, using the basis vectors
2512 : :
2513 : : !Recompute the metric tensor
2514 [ + + ]: 40278 : if(foundc==1)then
2515 [ + + ]: 154516 : do ii=1,3
2516 : : metmin(:,ii)=cell_base(1,:)*cell_base(1,ii)+&
2517 : : & cell_base(2,:)*cell_base(2,ii)+&
2518 [ + + ]: 502177 : & cell_base(3,:)*cell_base(3,ii)
2519 : : end do
2520 : : end if
2521 : :
2522 : : !Examine the angles and vector lengths
2523 : 40278 : ang90(:)=0
2524 [ + + ]: 40278 : if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1
2525 [ + + ]: 40278 : if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1
2526 [ + + ]: 40278 : if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1
2527 : 40278 : equal(:)=0
2528 [ + + ]: 40278 : if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1
2529 [ + + ]: 40278 : if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1
2530 [ + + ]: 40278 : if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1
2531 : :
2532 : : !DEBUG
2533 : : !write(std_out,*)' symlatt : recompute the metric tensor '
2534 : : !write(std_out,*)' ang90=',ang90
2535 : : !write(std_out,*)' equal=',equal
2536 : : !call flush(std_out)
2537 : : !ENDDEBUG
2538 : :
2539 : : !The axes will be aligned with the previously determined
2540 : : !basis vectors, EXCEPT for the tetragonal cell, see later
2541 : 40278 : axes(:,:)=cell_base(:,:)
2542 : :
2543 : 40278 : found=0
2544 : : !Check orthogonal conventional cells
2545 [ + + ]: 40278 : if(ang90(1)+ang90(2)+ang90(3)==3)then
2546 : :
2547 : : ! Cubic system
2548 [ + + ]: 32701 : if(equal(1)+equal(2)+equal(3)==3)then
2549 : : ! However, one-face centered is not admitted
2550 [ + + + - ]: 24669 : if(center==0 .or. center==-1 .or. center==-3)then
2551 : 16739 : iholohedry=7 ; found=1
2552 [ + + ]: 7930 : if(center==0)then
2553 : 6739 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cP (primitive cubic)'
2554 [ + + ]: 17930 : else if(center==-1)then
2555 : 1191 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cI (body-centered cubic)'
2556 : : else if(center==-3)then
2557 : 16739 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cF (face-centered cubic)'
2558 : : end if
2559 : : end if
2560 : : end if
2561 : :
2562 : : ! Tetragonal system
2563 [ + + + + : 8032 : if(found==0 .and. &
+ + ]
2564 : : & (equal(1)==1 .or. equal(2)==1 .or. equal(3)==1) )then
2565 : : ! However, one-face centered or face-centered is not admitted
2566 [ + + ]: 5431 : if(center==0 .or. center==-1)then
2567 : 5353 : iholohedry=4 ; found=1
2568 [ + + ]: 5353 : if(equal(1)==1)then
2569 [ + + + + : 24560 : axes(:,3)=cell_base(:,1) ; axes(:,1)=cell_base(:,2) ; axes(:,2)=cell_base(:,3)
+ + ]
2570 [ + + ]: 2897 : else if(equal(2)==1)then
2571 [ + + + + : 1260 : axes(:,3)=cell_base(:,2) ; axes(:,2)=cell_base(:,1) ; axes(:,1)=cell_base(:,3)
+ + ]
2572 [ + - ]: 2771 : else if(equal(3)==1)then
2573 : 2771 : axes(:,:)=cell_base(:,:)
2574 : : end if
2575 [ + + ]: 5353 : if(center==0)then
2576 : 4747 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is tP (primitive tetragonal)'
2577 : : else if(center==-1)then
2578 : 606 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is tI (body-centered tetragonal)'
2579 : : end if
2580 : : end if
2581 : : end if
2582 : :
2583 : : ! Orthorhombic system
2584 : : if(found==0)then
2585 : 2679 : iholohedry=3 ; found=1
2586 : 2679 : axes(:,:)=cell_base(:,:)
2587 [ + + ]: 2679 : if(center==0)then
2588 : 1751 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oP (primitive orthorhombic)'
2589 [ + + ]: 928 : else if(center==-1)then
2590 : 302 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oI (body-centered orthorhombic)'
2591 [ + + ]: 626 : else if(center==1 .or. center==2 .or. center==3)then
2592 : 526 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oC (one-face-centered orthorhombic)'
2593 [ + - ]: 100 : else if(center==-3)then
2594 : 100 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oF (face-centered orthorhombic)'
2595 : : end if
2596 : : end if
2597 : :
2598 : : else
2599 : :
2600 : : ! Hexagonal system
2601 [ + + + + : 7577 : if(found==0 .and. ang90(1)==1 .and. ang90(2)==1 .and. equal(3)==1 .and. (2*metmin(2,1)+metmin(1,1))<tolsym*metmin(1,1))then
+ - + - ]
2602 : 2196 : iholohedry=6 ; found=1
2603 : 2196 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is hP (primitive hexagonal)'
2604 : : end if
2605 : :
2606 : : ! Rhombohedral system
2607 : : if(found==0 .and. equal(1)+equal(2)+equal(3)==3 .and. &
2608 [ + + + - : 5381 : & abs(metmin(2,1)-metmin(3,2))<tolsym*metmin(2,2) .and. &
+ - ]
2609 : : & abs(metmin(2,1)-metmin(3,1))<tolsym*metmin(1,1) )then
2610 : 1581 : iholohedry=5 ; found=1
2611 : 1581 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is hR (rhombohedral)'
2612 : : end if
2613 : :
2614 : : ! Monoclinic system
2615 [ + + ]: 3800 : if(found==0 .and. ang90(1)+ang90(2)+ang90(3)==2 )then
2616 : 2151 : iholohedry=2 ; found=1
2617 [ + + ]: 2151 : if(center==0)then
2618 : 297 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is mP (primitive monoclinic)'
2619 [ + - ]: 1854 : else if(center==3)then
2620 : 1854 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is mC (one-face-centered monoclinic)'
2621 : : end if
2622 : : end if
2623 : :
2624 : : ! Triclinic system
2625 : : if(found==0)then
2626 : 1649 : iholohedry=1 ; found=1
2627 : 1649 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is aP (primitive triclinic)'
2628 : : end if
2629 : :
2630 : : end if
2631 : :
2632 : 40278 : call wrtout(iout,msg)
2633 : :
2634 : : !DEBUG
2635 : : !write(std_out,*)' symlatt : after checking conventional orthogonal cell '
2636 : : !call flush(std_out)
2637 : : !ENDDEBUG
2638 : :
2639 : : !--------------------------------------------------------------------------
2640 : : !Make sure that axes form a right-handed coordinate system
2641 : : !(Note : this should be done in the body of the routine,
2642 : : !by making changes that leave the sign of the mixed product of the three
2643 : : !vectors invariant)
2644 : : determinant=axes(1,1)*axes(2,2)*axes(3,3) &
2645 : : & +axes(1,2)*axes(2,3)*axes(3,1) &
2646 : : & +axes(1,3)*axes(3,2)*axes(2,1) &
2647 : : & -axes(1,1)*axes(3,2)*axes(2,3) &
2648 : : & -axes(1,3)*axes(2,2)*axes(3,1) &
2649 : 40278 : & -axes(1,2)*axes(2,1)*axes(3,3)
2650 [ + + ]: 40278 : if(determinant<zero)then
2651 [ + + + + ]: 71162 : axes(:,:)=-axes(:,:)
2652 : : end if
2653 : :
2654 : : !DEBUG
2655 : : !write(std_out,'(a,i4)')' symlatt : before itrial do loop, iholohedry= ',iholohedry
2656 : : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',&
2657 : : !& rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
2658 : : !call flush(std_out)
2659 : : !ENDDEBUG
2660 : :
2661 : : !--------------------------------------------------------------------------
2662 : : !Prefer symmetry axes on the same side as the primitive axes,
2663 : : !when the changes are allowed
2664 [ + - ]: 63878 : do itrial=1,100
2665 : :
2666 : : ! DEBUG
2667 : : ! write(std_out,'(a)')' '
2668 : : ! write(std_out,'(a,i5)')' symlatt : itrial do loop, itrial= ',itrial
2669 : : ! write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',&
2670 : : ! & axes(:,1),ch10,axes(:,2),ch10,axes(:,3)
2671 : : ! call flush(std_out)
2672 : : ! ENDDEBUG
2673 : :
2674 [ + + ]: 255512 : do ia=1,3
2675 : : scprods(ia,:)=axes(1,ia)*rprimd(1,:)+&
2676 : : & axes(2,ia)*rprimd(2,:)+&
2677 [ + + ]: 766536 : & axes(3,ia)*rprimd(3,:)
2678 [ + + ]: 766536 : norm2trial=sum(axes(:,ia)**2)
2679 [ + + ]: 830414 : scprods(ia,:)=scprods(ia,:)/sqrt(norm2trial)
2680 : : end do
2681 [ + + ]: 255512 : do ia=1,3
2682 [ + + ]: 766536 : norm2trial=sum(rprimd(:,ia)**2)
2683 [ + + ]: 830414 : scprods(:,ia)=scprods(:,ia)/sqrt(norm2trial)
2684 : : end do
2685 : :
2686 : : !DEBUG
2687 : : ! write(std_out,'(a,3f12.6)')' diagonal scalar products ',scprods(1,1),scprods(2,2),scprods(3,3)
2688 : : ! call flush(std_out)
2689 : : !ENDDEBUG
2690 : :
2691 : : ! One should now try all the generators of the
2692 : : ! proper rotations of each Bravais lattice, coupled with change of
2693 : : ! signs of each vector. This is not done systematically in what follows ...
2694 : : ! Here, the third axis is left unchanged
2695 [ + + ]: 63878 : if(iholohedry/=5)then
2696 [ + + + + ]: 59476 : if(scprods(1,1)<-tolsym .and. scprods(2,2)<-tolsym)then
2697 [ + + + + ]: 16450 : axes(:,1)=-axes(:,1) ; axes(:,2)=-axes(:,2)
2698 : : cycle
2699 : : end if
2700 : : end if
2701 : : ! The first (or second) axis is left unchanged
2702 [ + + ]: 61528 : if(iholohedry/=5 .and. iholohedry/=6)then
2703 [ + + + + ]: 54930 : if(scprods(2,2)<-tolsym .and. scprods(3,3)<-tolsym)then
2704 [ + + + + ]: 9093 : axes(:,2)=-axes(:,2) ; axes(:,3)=-axes(:,3)
2705 : : cycle
2706 : : end if
2707 [ + + + + ]: 53631 : if(scprods(1,1)<-tolsym .and. scprods(3,3)<-tolsym)then
2708 [ + + + + ]: 3318 : axes(:,1)=-axes(:,1) ; axes(:,3)=-axes(:,3)
2709 : : cycle
2710 : : end if
2711 : : end if
2712 : : ! Permutation of the three axis
2713 [ + + ]: 59755 : if(iholohedry==5 .or. iholohedry==7)then
2714 : 45727 : trace=scprods(1,1)+scprods(2,2)+scprods(3,3)
2715 [ + + ]: 45727 : if(trace+tolsym< scprods(1,2)+scprods(2,3)+scprods(3,1))then
2716 [ + + + + ]: 3591 : vecta(:)=axes(:,1) ; axes(:,1)=axes(:,3)
2717 [ + + + + ]: 3591 : axes(:,3)=axes(:,2); axes(:,2)=vecta(:)
2718 : : cycle
2719 : : end if
2720 [ + + ]: 45214 : if(trace+tolsym < scprods(1,3)+scprods(2,1)+scprods(3,2))then
2721 [ + + + + ]: 123648 : vecta(:)=axes(:,1) ; axes(:,1)=axes(:,2)
2722 [ + + + + ]: 123648 : axes(:,2)=axes(:,3); axes(:,3)=vecta(:)
2723 : : cycle
2724 : : end if
2725 : : ! This case is observed when the three new vectors
2726 : : ! are pointing opposite to the three original vectors
2727 : : ! One takes their opposite, then switch two of them, then process
2728 : : ! them again in the loop
2729 [ + + + + : 358150 : if(sum(scprods(:,:))<-tolsym)then
+ + ]
2730 [ + + ]: 5200 : axes(:,1)=-axes(:,1)
2731 [ + + ]: 5200 : vecta(:)=-axes(:,2)
2732 [ + + ]: 5200 : axes(:,2)=-axes(:,3)
2733 [ + + ]: 5200 : axes(:,3)=vecta(:)
2734 : : cycle
2735 : : end if
2736 : : end if
2737 : :
2738 : : ! Actually, for iholohedry==7, can test specifically all possibilities
2739 : : ! and take the best one.
2740 : : ! Not activated, because changing the order of symmetries in many tests !
2741 : : if(iholohedry==7 .and. .false.)then
2742 : : ! if(iholohedry==7)then
2743 : :
2744 : : !DEBUG
2745 : : !write(std_out,'(a,a)')ch10,' enter search of all possibilities for iholohedry==7 '
2746 : : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',&
2747 : : !& axes(:,1),ch10,axes(:,2),ch10,axes(:,3)
2748 : : !call flush(std_out)
2749 : : !ENDDEBUG
2750 : :
2751 : : do iaxis1=1,3
2752 : : axes_try(:,1)=axes(:,iaxis1)
2753 : : do iaxis2=1,2
2754 : : if(iaxis1==1)axes_try(:,2)=axes(:,1+iaxis2)
2755 : : if(iaxis1==2)axes_try(:,2)=axes(:,2*iaxis2-1)
2756 : : if(iaxis1==3)axes_try(:,2)=axes(:,iaxis2)
2757 : : if(iaxis2==1)axes_try(:,3)=axes(:,3)
2758 : : if(iaxis2==2)axes_try(:,3)=axes(:,1)
2759 : : if(iaxis1==1.and.iaxis2==2)axes_try(:,3)=axes(:,2)
2760 : : if(iaxis1==3.and.iaxis2==1)axes_try(:,3)=axes(:,2)
2761 : : do isign1=1,-1,-2
2762 : : axes_try(:,1)=-axes_try(:,1)
2763 : : do isign2=1,-1,-2
2764 : : axes_try(:,2)=-axes_try(:,2)
2765 : : determinant=axes_try(1,1)*axes_try(2,2)*axes_try(3,3) &
2766 : : & +axes_try(1,2)*axes_try(2,3)*axes_try(3,1) &
2767 : : & +axes_try(1,3)*axes_try(3,2)*axes_try(2,1) &
2768 : : & -axes_try(1,1)*axes_try(3,2)*axes_try(2,3) &
2769 : : & -axes_try(1,3)*axes_try(2,2)*axes_try(3,1) &
2770 : : & -axes_try(1,2)*axes_try(2,1)*axes_try(3,3)
2771 : : if(determinant<zero)axes_try(:,3)=-axes_try(:,3)
2772 : : do ia=1,3
2773 : : scprods(ia,:)=axes_try(1,ia)*rprimd(1,:)+&
2774 : : & axes_try(2,ia)*rprimd(2,:)+&
2775 : : & axes_try(3,ia)*rprimd(3,:)
2776 : : norm2trial=sum(axes_try(:,ia)**2)
2777 : : scprods(ia,:)=scprods(ia,:)/sqrt(norm2trial)
2778 : : end do
2779 : : do ia=1,3
2780 : : norm2trial=sum(rprimd(:,ia)**2)
2781 : : scprods(:,ia)=scprods(:,ia)/sqrt(norm2trial)
2782 : : end do
2783 : : trace=scprods(1,1)+scprods(2,2)+scprods(3,3)
2784 : : if(iaxis1==1.and.iaxis2==1.and.isign1==1.and.isign2==1)then
2785 : : trace_best=trace
2786 : : axes_best=axes_try
2787 : : else if (trace>trace_best+tolsym)then
2788 : : trace_best=trace
2789 : : axes_best=axes_try
2790 : : endif
2791 : : enddo ! isign2
2792 : : enddo ! isign1
2793 : : enddo ! iaxes2
2794 : : enddo ! iaxes1
2795 : : axes=axes_best
2796 : : endif ! iholohedry=7
2797 : 0 : exit
2798 : : end do
2799 : :
2800 : : !--------------------------------------------------------------------------
2801 : :
2802 : : !DEBUG
2803 : : !write(std_out,'(a,a)')ch10,' after order/sign optimization do-loop '
2804 : : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',&
2805 : : !& rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
2806 : : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',&
2807 : : !& axes(:,1),ch10,axes(:,2),ch10,axes(:,3)
2808 : : !call flush(std_out)
2809 : : !ENDDEBUG
2810 : :
2811 : : !Compute the coordinates of rprimd in the system defined by axes(:,:)
2812 : 40278 : call matr3inv(axes,axesinvt)
2813 [ + + ]: 161112 : do ii=1,3
2814 : : coord(:,ii)=rprimd(1,ii)*axesinvt(1,:)+ &
2815 : : & rprimd(2,ii)*axesinvt(2,:)+ &
2816 [ + + ]: 523614 : & rprimd(3,ii)*axesinvt(3,:)
2817 : : end do
2818 : :
2819 : : !Check that the coordinates are integers, or half-integer in
2820 : : !the case there is a centering, and generate integer coordinates
2821 [ + + ]: 161112 : do ii=1,3
2822 [ + + ]: 523614 : do jj=1,3
2823 : 362502 : val=coord(ii,jj)*fact
2824 [ - + ]: 362502 : if(abs(val-nint(val))>fact*two*tolsym)then
2825 : : write(msg,'(4a,a,3es18.10,a,a,3es18.10,a,a,3es18.10,a,a,i4)')&
2826 : 0 : 'One of the coordinates of rprimd in axes is non-integer,',ch10,&
2827 : 0 : 'or non-half-integer (if centering), within 2*tolsym.',ch10,&
2828 : 0 : 'coord=',coord(:,1),ch10,&
2829 : 0 : ' ',coord(:,2),ch10,&
2830 : 0 : ' ',coord(:,3),ch10,&
2831 : 0 : 'fact=',fact
2832 : 0 : ABI_BUG(msg)
2833 : : end if
2834 : 483336 : icoord(ii,jj)=nint(val)
2835 : : end do
2836 : : end do
2837 : :
2838 : : !Store the bravais lattice characteristics
2839 : 40278 : bravais(1)=iholohedry
2840 : 40278 : bravais(2)=center
2841 [ + + ]: 161112 : bravais(3:5)=icoord(1:3,1)
2842 [ + + ]: 161112 : bravais(6:8)=icoord(1:3,2)
2843 [ + + ]: 161112 : bravais(9:11)=icoord(1:3,3)
2844 : :
2845 : : !--------------------------------------------------------------------------
2846 : : !Initialize the set of symmetries
2847 : : !Bravais lattices are always invariant under identity and inversion
2848 : :
2849 : : !Identity and inversion
2850 [ + + + + : 1006950 : ptsymrel(:,:,1)=identity(:,:) ; ptsymrel(:,:,2)=-identity(:,:)
+ + + + ]
2851 : 40278 : nptsym=2
2852 : :
2853 : : !Keep this for IFCv70 compiler
2854 : : if(nptsym/=2)then
2855 : : write(msg,'(a,a,a,a)')ch10,&
2856 : : ' symlatt : BUG -',ch10,&
2857 : : ' Crazy error, compiler bug '
2858 : : call wrtout(std_out,msg)
2859 : : end if
2860 : :
2861 : : !--------------------------------------------------------------------------
2862 : : !Initialize some generators
2863 : : !gen6 is defined in a coordinated system with gamma=120 degrees
2864 : 40278 : gen6(:,:)=0 ; gen6(3,3)=1 ; gen6(1,1)=1 ; gen6(1,2)=-1 ; gen6(2,1)=1
2865 : 40278 : gen3(:,:)=0 ; gen3(1,2)=1 ; gen3(2,3)=1 ; gen3(3,1)=1
2866 : 40278 : gen2xy(:,:)=0 ; gen2xy(2,1)=1 ; gen2xy(1,2)=1; gen2xy(3,3)=1
2867 : 40278 : gen2y(:,:)=0 ; gen2y(1,1)=-1; gen2y(2,2)=1 ; gen2y(3,3)=-1
2868 : 40278 : gen2z(:,:)=0 ; gen2z(1,1)=-1; gen2z(2,2)=-1; gen2z(3,3)=1
2869 : :
2870 : : !--------------------------------------------------------------------------
2871 : :
2872 : : !Define the generators for each holohedry (inversion is already included)
2873 [ + + ]: 40278 : if(iholohedry==6)then
2874 : 2196 : ngen=2
2875 [ + + + + ]: 28548 : gen(:,:,1)=gen2xy(:,:) ; order(1)=2
2876 [ + + + + ]: 28548 : gen(:,:,2)=gen6(:,:) ; order(2)=6
2877 [ + + ]: 38082 : else if(iholohedry==5)then
2878 : 1581 : ngen=2
2879 [ + + + + ]: 20553 : gen(:,:,1)=gen2xy(:,:) ; order(1)=2
2880 [ + + + + ]: 20553 : gen(:,:,2)=gen3(:,:) ; order(2)=3
2881 : : else
2882 [ + + + + ]: 474513 : gen(:,:,1)=gen2y(:,:) ; order(1)=2
2883 [ + + + + ]: 474513 : gen(:,:,2)=gen2z(:,:) ; order(2)=2
2884 [ + + + + ]: 474513 : gen(:,:,3)=gen2xy(:,:) ; order(3)=2
2885 [ + + + + ]: 474513 : gen(:,:,4)=gen3(:,:) ; order(4)=3
2886 [ + + ]: 36501 : if(iholohedry<=4)ngen=iholohedry-1
2887 [ + + ]: 36501 : if(iholohedry==7)ngen=4
2888 : : end if
2889 : :
2890 : : !Build the point symmetry operations from generators, in the reduced system
2891 : : !of coordinates defined by axes(:,:)
2892 [ + + ]: 40278 : if(ngen/=0)then
2893 [ + + ]: 168427 : do igen=1,ngen
2894 [ + + ]: 1424012 : do isym=1+nptsym,order(igen)*nptsym
2895 : 1294214 : jsym=isym-nptsym
2896 [ + + ]: 5306654 : do ii=1,3
2897 : : ptsymrel(:,ii,isym)=gen(:,1,igen)*ptsymrel(1,ii,jsym)+ &
2898 : : & gen(:,2,igen)*ptsymrel(2,ii,jsym)+ &
2899 [ + + ]: 16824782 : & gen(:,3,igen)*ptsymrel(3,ii,jsym)
2900 : : end do
2901 : : end do
2902 : 168427 : nptsym=order(igen)*nptsym
2903 : :
2904 : : end do
2905 : : end if
2906 : :
2907 : : !--------------------------------------------------------------------------
2908 : :
2909 : : !Transform symmetry matrices in the system defined by rprimd
2910 : 40278 : call symrelrot(nptsym,axes,rprimd,ptsymrel,tolsym)
2911 : :
2912 : : !DEBUG
2913 : : !write(std_out,'(a)') ' symlatt : exit '
2914 : : !call flush(std_out)
2915 : : !stop
2916 : : !ENDDEBUG
2917 : :
2918 : 40278 : end subroutine symlatt
2919 : : !!***
2920 : :
2921 [ - + ]: 200354 : end module m_symfind
2922 : : !!***
|
