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 : 18161 : subroutine symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
100 : 18161 : prtvol, ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
101 : 35894 : 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 : 18161 : 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 : 18161 : real(dp),allocatable :: chrgat_(:)
140 : 18161 : real(dp),allocatable :: chrgatcl(:)
141 : 18161 : real(dp),allocatable :: local_nucdipmom(:,:,:),nucdipmomcl(:,:),nucdipmomred(:,:,:)
142 : 18161 : real(dp),allocatable :: spinatcl(:,:),spinatred(:,:)
143 : : !**************************************************************************
144 : :
145 : : !write(std_out,'(a)')' m_symfind%symfind : enter '; call flush(std_out)
146 : :
147 [ - + + - : 54483 : ABI_MALLOC(local_nucdipmom,(3,3,natom))
+ - + - -
+ - + ]
148 [ + + + + : 601029 : local_nucdipmom(:,:,:) = zero
+ + ]
149 [ + + ]: 18161 : if(present(nucdipmom)) then
150 [ + + + + ]: 194667 : 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 [ + - - + : 72644 : ABI_MALLOC(class,(natom+3,natom))
+ - - + +
- + - + -
- + - + ]
186 [ - + + - : 54483 : ABI_MALLOC(natomcl,(natom))
+ - + - -
+ - + ]
187 [ - + + - : 36322 : ABI_MALLOC(typecl,(natom))
- + - + ]
188 [ - + + - : 54483 : ABI_MALLOC(chrgat_,(natom))
+ - + - -
+ - + ]
189 [ - + + - : 36322 : ABI_MALLOC(chrgatcl,(natom))
- + - + ]
190 [ - + + - : 54483 : ABI_MALLOC(spinatcl,(3,natom))
+ - + - -
+ - + ]
191 [ - + + - : 36322 : ABI_MALLOC(nucdipmomcl,(3,natom))
- + - + ]
192 : :
193 : 18161 : tolsym2=tolsym**2
194 : :
195 [ + + ]: 62997 : chrgat_(:)=zero
196 [ + + ]: 18161 : if(present(chrgat))then
197 [ + + ]: 62127 : 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 [ + + ]: 62997 : 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 [ + + ]: 480 : 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 : 120 : call random_number(ndtest)
214 : :
215 : : ! vector 2 is constructed to be orthogonal to original nuclear dipole moment vector
216 : 120 : 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 : 62997 : 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 : 18161 : 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 : 18161 : nclass=1
229 : 18161 : natomcl(1)=1
230 : 18161 : typecl(1)=typat(1)
231 : 18161 : chrgatcl(1)=chrgat_(1)
232 [ + + ]: 72644 : spinatcl(:,1)=spinat(:,1)
233 [ + + ]: 72644 : nucdipmomcl(:,1)=local_nucdipmom(:,1,1)
234 : 18161 : class(1,1)=1
235 [ + + ]: 18161 : 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 : 39037 : foundcl=0
242 [ + + ]: 39037 : 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 [ + + ]: 39037 : if( typat(iatom)==typecl(iclass)) then
246 : 20658 : test_samechrg= (abs(chrgat_(iatom)-chrgatcl(iclass))<tolsym)
247 [ + + ]: 20658 : 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 : 648 : spinat2 =spinat(1,iatom)**2+spinat(2,iatom)**2+spinat(3,iatom)**2
251 : 648 : spinatcl2=spinatcl(1,iclass)**2+spinatcl(2,iclass)**2+spinatcl(3,iclass)**2
252 : 648 : 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 [ + - + - : 20658 : & 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 [ + + + - ]: 20658 : 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 : 20576 : natomcl(iclass)=natomcl(iclass)+1
267 : 20576 : 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 : 13760 : if(foundcl==0)then
275 : 6099 : nclass=nclass+1
276 : 6099 : natomcl(nclass)=1
277 : 6099 : typecl(nclass)=typat(iatom)
278 : 6099 : chrgatcl(nclass)=chrgat_(iatom)
279 [ + + ]: 24396 : spinatcl(:,nclass)=spinat(:,iatom)
280 [ + + ]: 24396 : nucdipmomcl(:,nclass)=local_nucdipmom(:,1,iatom)
281 : 6099 : 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 : 18161 : iclass0=1
305 : 18161 : natom0=natomcl(1)
306 : 18161 : spinatcl20=spinatcl(1,1)**2+spinatcl(2,1)**2+spinatcl(3,1)**2
307 : 18161 : nucdipmomcl20=nucdipmomcl(1,1)**2+nucdipmomcl(2,1)**2+nucdipmomcl(3,1)**2
308 [ + + ]: 18161 : if(nclass>1)then
309 [ + + ]: 11830 : do iclass=2,nclass
310 : 6099 : spinatcl2=spinatcl(1,iclass)**2+spinatcl(2,iclass)**2+spinatcl(3,iclass)**2
311 : 6099 : 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 [ + + + + : 11830 : & .or. (nucdipmomcl20<tolsym .and. nucdipmomcl2>tolsym)) then
- + - + -
- + + + +
+ + + + ]
317 : 6099 : iclass0=iclass
318 : 6099 : natom0=natomcl(iclass)
319 : 6099 : spinatcl20=spinatcl2
320 : 6099 : nucdipmomcl20=nucdipmomcl2
321 : : end if
322 : : end do
323 : : end if
324 : :
325 : 18161 : printed=0
326 : :
327 : : !If non-collinear spinat have to be used, transfer them in reduced coordinates
328 [ + + ]: 18161 : if (nspden==4) then
329 [ - + + - : 1096 : ABI_MALLOC(spinatred,(3,natom))
- + - + ]
330 [ + + ]: 1900 : do iatom=1,natom
331 [ + + ]: 5956 : do ii=1,3
332 : 5408 : 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 [ - + + - : 36322 : ABI_MALLOC(nucdipmomred,(3,3,natom))
- + - + ]
352 [ + + ]: 62997 : do iatom=1,natom
353 [ + + ]: 197505 : 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 : 18161 : nsym=0
362 [ + + ]: 761785 : do isym=1,nptsym
363 : :
364 : : !write(std_out,'(a,i4)')' m_symfind%symfind : enter loop isym=',isym; call flush(std_out)
365 : :
366 [ + + ]: 743624 : 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 [ + + ]: 2939104 : & ptsymrel(:,3,isym)*invardir_red(3)
371 [ + + ]: 2939104 : diff(:)=invardir_red(:)-invardir_red_rot(:)
372 [ + + ]: 734776 : 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 [ + + ]: 740664 : if(present(invaraxial_red))then
377 : : !! CALL MatrixHyb_getDet(symref_hyb(isym), detR)
378 [ + + + + ]: 9513608 : mm(:,:) = ptsymrel(:,:,isym)
379 : 731816 : call mati3det(mm,detR)
380 : : invaraxial_red_rot(:) = detR * (&
381 : : ptsymrel(:,1,isym)*invaraxial_red(1) + &
382 : : & ptsymrel(:,2,isym)*invaraxial_red(2) + &
383 [ + + ]: 2927264 : & ptsymrel(:,3,isym)*invaraxial_red(3) )
384 [ + + ]: 2927264 : diff(:)=invaraxial_red(:)-invaraxial_red_rot(:)
385 : :
386 [ + + ]: 731816 : 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 : 739880 : & ptsymrel(3,2,isym)*ptsymrel(2,3,isym)*ptsymrel(1,1,isym))
397 [ + + + + ]: 739880 : 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 [ + + ]: 739872 : if(present(invar_z))then
404 [ + + ]: 731032 : 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 [ + + - + ]: 738432 : 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 : 738432 : iatom0=class(1,iclass0)
430 : : xred0(:)=ptsymrel(:,1,isym)*xred(1,iatom0)+ &
431 : : ptsymrel(:,2,isym)*xred(2,iatom0)+ &
432 [ + + ]: 2953728 : ptsymrel(:,3,isym)*xred(3,iatom0)
433 [ + + ]: 738432 : if (nspden/=4) then
434 [ + + ]: 2881216 : spinat0(:)=spinat(:,iatom0)
435 : : else
436 : : spinat0(:)=det*(&
437 : : ptsymrel(:,1,isym)*spinatred(1,iatom0)+ &
438 : : & ptsymrel(:,2,isym)*spinatred(2,iatom0)+ &
439 [ + + ]: 72512 : & 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 [ + + ]: 2050141 : do ii=1,natom0
450 : : !write(std_out,'(a,2i4)')' symfind : loop isym,ii=',isym,ii
451 : 1293548 : iatom1=class(ii,iclass0)
452 : :
453 : : ! The tentative translation is found
454 [ + + ]: 5174192 : trialnons(:)=xred(:,iatom1)-xred0(:)
455 : : ! Compare the spinat vectors
456 [ + + ]: 1293548 : if (nspden/=4) then
457 [ + + ]: 5073232 : symspinat1(:)=spinat(:,iatom1)
458 : : else
459 [ + + ]: 100960 : symspinat1(:)=spinatred(:,iatom1)
460 : : end if
461 : :
462 : : !write(std_out,'(a,6f10.4)')' symspinat1,spinat0=',symspinat1(:),spinat0(:)
463 : :
464 : 1293548 : trialafm=1
465 [ + + + + ]: 5174192 : if(sum(abs(symspinat1(:)-spinat0(:)))>tolsym)then
466 : 20536 : trialafm=-1
467 [ + + + - ]: 20536 : if(nspden==4 .and. afm_noncoll==0)cycle
468 [ + + + + ]: 82144 : if(sum(abs(symspinat1(:)+spinat0(:)))>tolsym)cycle
469 : : endif
470 : :
471 [ + + - + ]: 5122992 : 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 [ + + ]: 1280748 : if(present(invar_z))then
481 [ + + + - ]: 1267536 : 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 [ + + ]: 2103026 : do iclass=1,nclass
490 [ + + ]: 4212030 : do jj=1,natomcl(iclass)
491 : :
492 : 3389752 : 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 [ + + ]: 13559008 : & ptsymrel(:,3,isym)*xred(3,iatom2)+ trialnons(:)
497 : : ! Generate the tentative symmetric spinat of iatom2
498 [ + + ]: 3389752 : if (nspden/=4) then
499 [ + + ]: 13447776 : symspinat2(:)=trialafm*spinat(:,iatom2)
500 : : else
501 : : symspinat2(:)=trialafm*det*(ptsymrel(:,1,isym)*spinatred(1,iatom2)+ &
502 : : & ptsymrel(:,2,isym)*spinatred(2,iatom2)+ &
503 [ + + ]: 111232 : & ptsymrel(:,3,isym)*spinatred(3,iatom2))
504 : : end if
505 : : ! Generate the tentative symmetric nucdipmom of iatom2
506 [ + + ]: 13559008 : do kk = 1, 3
507 : : symnucdipmom2red(:,kk)=ptsymrel(:,1,isym)*nucdipmomred(1,kk,iatom2)+ &
508 : : & ptsymrel(:,2,isym)*nucdipmomred(2,kk,iatom2)+ &
509 [ + + ]: 40677024 : & ptsymrel(:,3,isym)*nucdipmomred(3,kk,iatom2)
510 : : ! transform back to cart coords for final comparison to nucdipmom
511 [ + + + + : 166097848 : 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 [ + + ]: 28490242 : do kk=1,natomcl(iclass)
521 : :
522 : 27802703 : found3=1
523 : 27802703 : iatom3=class(kk,iclass)
524 : : ! Check the location
525 : 27802703 : diffr1=xred(1,iatom3)-symxred2(1)
526 : 27802703 : diff1=diffr1-nint(diffr1)
527 [ + + ]: 27802703 : if(diff1**2>tolsym2)then
528 : : found3=0 ; cycle
529 : : else
530 : 6125266 : diffr2=xred(2,iatom3)-symxred2(2)
531 : 6125266 : diff2=diffr2-nint(diffr2)
532 [ + + ]: 6125266 : if(diff2**2>tolsym2)then
533 : : found3=0 ; cycle
534 : : else
535 : 3323323 : diffr3=xred(3,iatom3)-symxred2(3)
536 : 3323323 : diff3=diffr3-nint(diffr3)
537 [ + + ]: 3323323 : 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 [ + + ]: 2705287 : if (nspden/=4) then
544 [ + + ]: 10729836 : diff(:)=spinat(:,iatom3)-symspinat2(:)
545 : : else
546 [ + + ]: 91312 : diff(:)=spinatred(:,iatom3)-symspinat2(:)
547 : : end if
548 [ + + ]: 2705287 : 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 : 2705053 : 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 : 2705053 : call acrossb(symnucdipmom2cart(1:3,2),symnucdipmom2cart(1:3,3),hand2)
558 : :
559 [ + + ]: 10820212 : diff(:)=hand3(:)-hand2(:)
560 [ + + + + : 10812364 : 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 : 687539 : if(found3==1)exit
566 : : end do ! End loop over iatom3
567 : :
568 [ + + ]: 4212030 : if(found3==0)then
569 : : trialok=0
570 : : exit
571 : : end if
572 : : end do ! End loop over iatom2
573 : :
574 [ + + ]: 2103026 : 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 [ + + ]: 2024372 : if(trialok==1)then
581 : 593209 : nsym=nsym+1
582 [ - + ]: 593209 : 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 [ + + + + ]: 7711717 : symrel(:,:,nsym)=ptsymrel(:,:,isym)
592 : 593209 : symafm(nsym)=trialafm
593 [ + + ]: 2372836 : 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 : 18161 : ABI_FREE(class)
602 : 18161 : ABI_FREE(natomcl)
603 : 18161 : ABI_FREE(chrgat_)
604 : 18161 : ABI_FREE(chrgatcl)
605 : 18161 : ABI_FREE(spinatcl)
606 : 18161 : ABI_FREE(typecl)
607 [ - + ]: 18161 : ABI_FREE(local_nucdipmom)
608 : 18161 : ABI_FREE(nucdipmomcl)
609 [ - + ]: 18161 : ABI_FREE(nucdipmomred)
610 [ + + ]: 18161 : if (nspden==4) then
611 [ - + ]: 548 : 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 [ + + + - ]: 18161 : if(present(ierr) .or. nsym<=384)then
618 : 18161 : 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 [ + + + + ]: 18162 : if (ierr_/=0) call print_symmetries([std_out], nsym, symrel, tnons, symafm)
626 : :
627 [ + + ]: 18161 : if(.not.present(ierr))then
628 [ - + ]: 9159 : ABI_CHECK(ierr_==0,"Error in group closure")
629 : : else
630 : 9002 : 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 [ + + + + : 54269 : 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 : 9001 : subroutine symfind_expert(gprimd,msym,natom,nptsym,nspden,nsym,&
694 : 9001 : 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 : 9001 : integer,allocatable :: indsym(:,:,:),symrec(:,:,:)
721 : 9001 : real(dp),allocatable :: tnons_new(:,:)
722 : : !**************************************************************************
723 : :
724 : : ! write(std_out,*)' m_symfind%symfind_expert : enter '
725 : :
726 : 9001 : use_inversion=1
727 : : if (usepaw == 1 .and. (nspden==4.or.pawspnorb>0)) then
728 : : !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 : : ! problem should be fixed now
731 : : use_inversion=1
732 : : end if
733 : :
734 : : ! write(std_out,*)' m_symfind%symfind_expert : before call symfind (1) '
735 : :
736 : : call symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
737 : : prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
738 : 9001 : chrgat=chrgat,nucdipmom=nucdipmom,ierr=ierr,invardir_red=invardir_red,invaraxial_red=invaraxial_red,invar_z=invar_z)
739 : :
740 : : ! write(std_out,*)' m_symfind%symfind_expert : after call symfind (1) '
741 : :
742 : : !If the group closure is not obtained, which should be exceptional, try with a larger tolsym (three times larger)
743 [ + + ]: 9001 : if(ierr/=0)then
744 : 1 : ABI_WARNING('Will try to obtain group closure by using a tripled tolsym.')
745 : : call symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
746 : : prtvol,ptsymrel,spinat,symafm,symrel,tnons,three*tolsym,typat,use_inversion,xred,&
747 : 1 : chrgat=chrgat,nucdipmom=nucdipmom,ierr=ierr,invardir_red=invardir_red,invaraxial_red=invaraxial_red,invar_z=invar_z)
748 [ - + ]: 1 : ABI_CHECK(ierr==0,"Error in group closure")
749 : 1 : ABI_WARNING('Succeeded to obtain group closure by using a tripled tolsym.')
750 : : endif
751 : :
752 : : ! If the tolerance on symmetries is bigger than 1.e-8, symmetrize tnons for gliding or screw operations,
753 : : ! symmetrize the atomic positions and recompute the symmetry operations
754 [ + + ]: 9001 : if(tolsym>1.00001e-8)then
755 : :
756 : 8945 : call symmetrize_tnons(nsym,symrel,tnons,tolsym)
757 [ - + + - : 35780 : ABI_MALLOC(indsym,(4,natom,nsym))
- + + - +
- + - - +
- + ]
758 [ - + + - : 26835 : ABI_MALLOC(symrec,(3,3,nsym))
+ - + - -
+ + - ]
759 [ + + ]: 300861 : do isym=1,nsym
760 : 300861 : call mati3inv(symrel(:,:,isym),symrec(:,:,isym))
761 : : end do
762 : 8945 : call symatm(indsym,natom,nsym,symrec,tnons,tolsym,typat,xred)
763 : 8945 : call symmetrize_xred(natom,nsym,symrel,tnons,xred,indsym=indsym)
764 [ - + ]: 8945 : ABI_FREE(indsym)
765 [ - + ]: 8945 : ABI_FREE(symrec)
766 : :
767 [ + + ]: 8945 : if(print_comment_tolsym==1)then
768 : : write(msg,'(a,es12.3,18a)')&
769 : 1371 : 'The tolerance on symmetries =',tolsym,' is bigger than 1.0e-8.',ch10,&
770 : 1371 : 'In order to avoid spurious effects, the atomic coordinates have been',ch10,&
771 : 1371 : 'symmetrized before storing them in the dataset internal variable.',ch10,&
772 : 1371 : 'So, do not be surprised by the fact that your input variables (xcart, xred, ...)',ch10,&
773 : 1371 : 'do not correspond exactly to the ones echoed by ABINIT, the latter being used to do the calculations.',ch10,&
774 : 1371 : 'This is not a problem per se.',ch10,&
775 : 1371 : 'Still, in order to avoid this symmetrization (e.g. for specific debugging/development),',&
776 : 1371 : ' decrease tolsym to 1.0e-8 or lower,',ch10,&
777 : 1371 : 'or (much preferred) use input primitive vectors that are accurate to better than 1.0e-8.',ch10,&
778 : 2742 : 'This message will only be printed once, even if there are other datasets where tolsym is bigger than 1.0e-8.'
779 : 1371 : ABI_COMMENT(msg)
780 : 1371 : print_comment_tolsym=0
781 : : endif
782 : :
783 : : !write(std_out,*)' m_symfind%symfind_expert : before call symfind (3) '
784 : :
785 : : call symfind(gprimd,msym,natom,nptsym,nspden,nsym,&
786 : : prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,&
787 : 8945 : chrgat=chrgat,nucdipmom=nucdipmom,invardir_red=invardir_red,invaraxial_red=invaraxial_red,invar_z=invar_z)
788 : :
789 : : ! write(std_out,*)' m_symfind%symfind_expert : after call symfind (3) '
790 : :
791 : : !Needs one more resymmetrization, for the tnons
792 [ - + + - : 26835 : ABI_MALLOC(tnons_new,(3,nsym))
+ - + - -
+ - + ]
793 : :
794 : : call symmetrize_xred(natom,nsym,symrel,tnons,xred,&
795 : 8945 : fixed_mismatch=fixed_mismatch,mismatch_fft_tnons=mismatch_fft_tnons,tnons_new=tnons_new,tolsym=tolsym)
796 [ + + + + ]: 1176609 : tnons(:,1:nsym)=tnons_new(:,:)
797 [ - + ]: 26835 : ABI_FREE(tnons_new)
798 : : end if ! tolsym >1.00001e-8
799 : :
800 : : ! write(std_out,*)' m_symfind%symfind_expert : exit '
801 : :
802 : 9001 : end subroutine symfind_expert
803 : : !!***
804 : :
805 : : !!****f* m_symfind/symanal
806 : : !! NAME
807 : : !! symanal
808 : : !!
809 : : !! FUNCTION
810 : : !! Find the space group, Bravais lattice, including Shubnikov characteristics
811 : : !! from the list of symmetries (including magnetic characteristics), and lattice parameters
812 : : !! Warning: the recognition of the space group might not yet work for the
813 : : !! Shubnikov group of type IV
814 : : !!
815 : : !! INPUTS
816 : : !! chkprim= if 1 then stop if the cell is not primitive
817 : : !! msym=default maximal number of symmetries
818 : : !! nsym=actual number of symmetries
819 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
820 : : !! symafm(1:msym)=(anti)ferromagnetic part of symmetry operations
821 : : !! symrel(3,3,1:msym)=symmetry operations in real space in terms
822 : : !! of primitive translations
823 : : !! tnons(3,1:msym)=nonsymmorphic translations for symmetry operations
824 : : !! tolsym=tolerance for the symmetry operations
825 : : !! [verbose]= if true, will list the symmetry operation labels
826 : : !!
827 : : !! OUTPUT
828 : : !! bravais(11)=characteristics of Bravais lattice (see symlatt.F90)
829 : : !! genafm(3)=magnetic translation generator (in case of Shubnikov group type IV)
830 : : !! ptgroupma = magnetic point group number
831 : : !! spgroup=symmetry space group
832 : : !!
833 : : !! SOURCE
834 : :
835 : 12033 : subroutine symanal(bravais,chkprim,genafm,msym,nsym,ptgroupma,rprimd,spgroup,symafm,symrel,tnons,tolsym, &
836 : : verbose) ! optional
837 : :
838 : : !Arguments ------------------------------------
839 : : !scalars
840 : : integer,intent(in) :: chkprim,msym,nsym
841 : : integer,intent(out) :: ptgroupma,spgroup
842 : : real(dp),intent(in) :: tolsym
843 : : logical,optional,intent(in) :: verbose
844 : : !arrays
845 : : integer,intent(out) :: bravais(11)
846 : : integer,intent(in) :: symafm(msym),symrel(3,3,msym)
847 : : real(dp),intent(in) :: rprimd(3,3), tnons(3,msym)
848 : : real(dp),intent(out) :: genafm(3)
849 : :
850 : : !Local variables-------------------------------
851 : : !scalars
852 : : integer, parameter :: maxsym=192
853 : : ! In this routine, maxsym is used either to determine the ptsymrel from the rprimd (routine symlatt),
854 : : ! so, it might be up to 192 = 4*48 for FCC, and also to define the maximum number of symmetry operation labels,
855 : : ! but only in case the cell is primitive, which gives the same upper bound. Thus in this routine, msym might
856 : : ! be equal to nsym.
857 : : integer :: iholohedry_nomagn,isym,isym_nomagn,multi, nptsym,nsym_nomagn,shubnikov
858 : : logical :: verbose_
859 : : character(len=5) :: ptgroup,ptgroupha
860 : : character(len=500) :: msg
861 : : !arrays
862 : : integer :: identity(3,3)
863 : 12033 : integer,allocatable :: ptsymrel(:,:,:),symrel_nomagn(:,:,:)
864 : 12033 : real(dp),allocatable :: tnons_nomagn(:,:)
865 : : character(len=128) :: labels(maxsym)
866 : : ! *************************************************************************
867 : :
868 : : !write(std_out,*)' symanal : enter'
869 : : !write(std_out,*)' symanal : chkprim =',chkprim
870 : : !write(std_out,*)' symanal : nsym=',nsym
871 : : !do isym=1,nsym
872 : : ! write(std_out,*)' symanal : symrel=',symrel(1:3,1:3,isym)
873 : : !enddo
874 : :
875 : 12033 : verbose_=.false.
876 [ + + ]: 12033 : if (present(verbose)) verbose_=verbose
877 : :
878 : : !This routine finds the Bravais characteristics, without actually looking at the symmetry operations.
879 [ - + ]: 12033 : ABI_MALLOC(ptsymrel,(3,3,maxsym))
880 : 12033 : call symlatt(bravais,dev_null,maxsym,nptsym,ptsymrel,rprimd,tolsym)
881 [ - + ]: 12033 : ABI_FREE(ptsymrel)
882 : :
883 : : ! Check whether the cell is primitive or not.
884 : 12033 : call chkprimit(chkprim,multi,nsym,symafm,symrel)
885 : :
886 : 12033 : spgroup=0 ; ptgroupma=0 ; genafm(:)=zero
887 : :
888 [ + + ]: 12033 : if (multi>1) then ! Modify bravais if the cell is not primitive ; no determination of the space group
889 : 289 : bravais(1)=-bravais(1)
890 : : else
891 : :
892 : : ! The cell is primitive, so that the space group can be
893 : : ! determined. Need to distinguish Fedorov and Shubnikov groups.
894 : : ! Do not distinguish Shubnikov types I and II.
895 : : ! Also identify genafm, in case of Shubnikov type IV
896 : 11744 : identity(:,:)=reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/))
897 : 11744 : shubnikov=1
898 [ + + ]: 300418 : do isym=1,nsym
899 [ + + ]: 300418 : if(symafm(isym)==-1)then
900 : 5968 : shubnikov=3
901 [ + + + + : 77584 : if(sum(abs(symrel(:,:,isym)-identity(:,:)))==0)then
+ + ]
902 : 611 : shubnikov=4
903 [ + + ]: 2444 : genafm(:)=tnons(:,isym)
904 : : !write(std_out,*)' isym=',isym
905 : : !write(std_out,*)' symrel(:,:,isym)',symrel(:,:,isym)
906 : : !write(std_out,*)' tnons(:,isym)',tnons(:,isym)
907 : : !write(std_out,*)' symafm(isym)',symafm(isym)
908 : : exit
909 : : end if
910 : : end if
911 : : end do
912 : :
913 [ + + ]: 11744 : if(shubnikov/=1)then
914 [ + + ]: 1515 : if(shubnikov==3)write(msg, '(a)' )' Shubnikov space group type III'
915 [ + + ]: 1515 : if(shubnikov==4)write(msg, '(a)' )' Shubnikov space group type IV'
916 : 1515 : call wrtout(std_out,msg)
917 : : end if
918 : :
919 [ + + ]: 11744 : if(shubnikov==1 .or. shubnikov==3)then
920 : : ! Find the correct Bravais characteristics and point group
921 : : ! Should also be used for Shubnikov groups of type IV ...
922 : 11133 : call symbrav(bravais,msym,nsym,ptgroup,rprimd,symrel,tolsym)
923 : :
924 : : ! Find the space group
925 : 11133 : call symspgr(bravais,labels,nsym,spgroup,symrel,tnons,tolsym)
926 : :
927 [ + + ]: 11133 : if (verbose_) then
928 [ + + ]: 2275 : do isym=1,nsym
929 [ + - ]: 2203 : write(msg,'(a,i3,2a)')' symanal : the symmetry operation no. ',isym,' is ',trim(labels(isym))
930 : 2275 : call wrtout(std_out,msg)
931 : : enddo
932 : : endif
933 : :
934 : : end if
935 : :
936 [ + + ]: 11133 : if (shubnikov/=1) then
937 : :
938 : : ! Determine nonmagnetic symmetry operations
939 : 1515 : nsym_nomagn=nsym/2
940 [ - + + - : 4545 : ABI_MALLOC(symrel_nomagn,(3,3,nsym_nomagn))
+ - + - -
+ - + ]
941 [ - + + - : 4545 : ABI_MALLOC(tnons_nomagn,(3,nsym_nomagn))
+ - + - -
+ - + ]
942 : 1515 : isym_nomagn=0
943 [ + + ]: 24243 : do isym=1,nsym
944 [ + + ]: 24243 : if(symafm(isym)==1)then
945 : 11364 : isym_nomagn=isym_nomagn+1
946 [ + + + + ]: 147732 : symrel_nomagn(:,:,isym_nomagn)=symrel(:,:,isym)
947 [ + + ]: 45456 : tnons_nomagn(:,isym_nomagn)=tnons(:,isym)
948 : : end if
949 : : end do
950 : :
951 [ + + ]: 1515 : if (shubnikov==3) then
952 : :
953 : : ! write(std_out,*)' symanal : will enter symbrav with halved symmetry set'
954 : : ! write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)'
955 : : ! do isym=1,nsym_nomagn
956 : : ! write(std_out,'(i3,2x,9i3,3es12.2,i3)')isym,symrel_nomagn(:,:,isym),tnons_nomagn(:,isym)
957 : : ! end do
958 : :
959 : : ! Find the point group of the halved symmetry set
960 : 904 : call symptgroup(iholohedry_nomagn,nsym_nomagn,ptgroupha,symrel_nomagn)
961 : :
962 : : ! Deduce the magnetic point group (ptgroupma) from ptgroup and ptgroupha
963 : 904 : call getptgroupma(ptgroup,ptgroupha,ptgroupma)
964 : :
965 : 611 : else if(shubnikov==4)then
966 : :
967 : : ! Find the Fedorov space group of the halved symmetry set
968 : 611 : call symspgr(bravais,labels,nsym_nomagn,spgroup,symrel_nomagn,tnons_nomagn,tolsym)
969 : :
970 : : ! The magnetic translation generator genafm has already been determined
971 : : ! write(std_out,*)' genafm =',genafm, ' spgroup=',spgroup
972 [ - + ]: 611 : if (verbose_) then
973 : 0 : write(msg, '(a)' )' Select only the non-magnetic symmetry operations '
974 : 0 : call wrtout(std_out,msg)
975 : :
976 [ # # ]: 0 : do isym=1,nsym
977 [ # # ]: 0 : if(symafm(isym)==1)then
978 : 0 : isym_nomagn=isym_nomagn+1
979 [ # # ]: 0 : write(msg,'(a,i3,2a)')' symspgr : the symmetry operation no. ',isym,' is ',trim(labels(isym_nomagn))
980 : 0 : call wrtout(std_out,msg)
981 : : endif
982 : : enddo
983 : : endif
984 : :
985 : : end if
986 : :
987 [ - + ]: 1515 : ABI_FREE(symrel_nomagn)
988 [ - + ]: 1515 : ABI_FREE(tnons_nomagn)
989 : : end if ! Shubnikov groups
990 : :
991 : : end if
992 : :
993 : 12033 : end subroutine symanal
994 : : !!***
995 : :
996 : : !!****f* m_symfind/symbrav
997 : : !! NAME
998 : : !! symbrav
999 : : !!
1000 : : !! FUNCTION
1001 : : !! From the list of symmetry operations, and the lattice vectors,
1002 : : !! determine the Bravais information (including the holohedry, the centering,
1003 : : !! the coordinate of the primitive vectors in the conventional vectors),
1004 : : !! as well as the point group.
1005 : : !!
1006 : : !! INPUTS
1007 : : !! msym=dimension of symrel
1008 : : !! nsym=actual number of symmetries
1009 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1010 : : !! symrel(3,3,msym)=symmetry operations in real space in terms
1011 : : !! of primitive translations
1012 : : !! tolsym=tolerance for the symmetries
1013 : : !!
1014 : : !! OUTPUT
1015 : : !! bravais(11): bravais(1)=iholohedry
1016 : : !! bravais(2)=center
1017 : : !! bravais(3:11)=coordinates of rprimd in the axes
1018 : : !! of the conventional bravais lattice (*2 if center/=0)
1019 : : !! ptgroup=symmetry point group
1020 : : !! [axis(3)]=Invariant axis in the conventional vector coordinates
1021 : : !! Set to (/0,0,0/) if the lattice belongs to the same holohedry as the lattice+atoms (+electric field + ...).
1022 : : !!
1023 : : !! SOURCE
1024 : :
1025 : 11133 : subroutine symbrav(bravais,msym,nsym,ptgroup,rprimd,symrel,tolsym,axis)
1026 : :
1027 : : !Arguments ------------------------------------
1028 : : !scalars
1029 : : integer,intent(in) :: msym,nsym
1030 : : real(dp),intent(in) :: tolsym
1031 : : character(len=5),intent(out) :: ptgroup
1032 : : !arrays
1033 : : integer,intent(in) :: symrel(3,3,msym)
1034 : : integer,optional,intent(out) :: axis(3)
1035 : : integer,intent(out) :: bravais(11)
1036 : : real(dp),intent(in) :: rprimd(3,3)
1037 : :
1038 : : !Local variables-------------------------------
1039 : : !scalars
1040 : : integer :: iaxis,ii,bravais1now,ideform,iholohedry,invariant,isym
1041 : : integer :: jaxis,next_stage,nptsym,problem,maxsym
1042 : : integer, parameter :: naxes_ortho=22, naxes_hexa=7
1043 : : real(dp) :: norm,scprod
1044 : : character(len=500) :: msg
1045 : : !arrays
1046 : : integer :: identity(3,3),axis_trial(3),hexa_axes(3,naxes_hexa),ortho_axes(3,naxes_ortho)
1047 : 11133 : integer,allocatable :: ptsymrel(:,:,:),symrelconv(:,:,:)
1048 : : real(dp) :: axes(3,3),axis_cart(3),axis_red(3)
1049 : : real(dp) :: rprimdconv(3,3),rprimdtry(3,3),rprimdnow(3,3)
1050 : : real(dp) :: rprimdconv_invt(3,3)
1051 : : !**************************************************************************
1052 : :
1053 : : !DEBUG
1054 : : !write(std_out,*)' symbrav : enter '
1055 : : !call flush(std_out)
1056 : : !ENDDEBUG
1057 : :
1058 : : identity(:,:)=0
1059 : 11133 : identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
1060 : :
1061 : 11133 : ortho_axes(:,:)=0
1062 : 11133 : ortho_axes(1,1)=1
1063 : 11133 : ortho_axes(2,2)=1
1064 : 11133 : ortho_axes(3,3)=1
1065 [ + + ]: 44532 : ortho_axes(:,4)=(/0,1,1/)
1066 [ + + ]: 44532 : ortho_axes(:,5)=(/1,0,1/)
1067 [ + + ]: 44532 : ortho_axes(:,6)=(/1,1,0/)
1068 [ + + ]: 44532 : ortho_axes(:,7)=(/0,1,-1/)
1069 [ + + ]: 44532 : ortho_axes(:,8)=(/-1,0,1/)
1070 [ + + ]: 44532 : ortho_axes(:,9)=(/1,-1,0/)
1071 [ + + ]: 44532 : ortho_axes(:,10)=(/0,1,2/)
1072 [ + + ]: 44532 : ortho_axes(:,11)=(/2,0,1/)
1073 [ + + ]: 44532 : ortho_axes(:,12)=(/1,2,0/)
1074 [ + + ]: 44532 : ortho_axes(:,13)=(/1,1,1/)
1075 [ + + ]: 44532 : ortho_axes(:,14)=(/-1,1,1/)
1076 [ + + ]: 44532 : ortho_axes(:,15)=(/1,-1,1/)
1077 [ + + ]: 44532 : ortho_axes(:,16)=(/1,1,-1/)
1078 [ + + ]: 44532 : ortho_axes(:,17)=(/2,1,1/)
1079 [ + + ]: 44532 : ortho_axes(:,18)=(/1,2,1/)
1080 [ + + ]: 44532 : ortho_axes(:,19)=(/1,1,2/)
1081 [ + + ]: 44532 : ortho_axes(:,20)=(/2,1,-1/)
1082 [ + + ]: 44532 : ortho_axes(:,21)=(/-1,2,1/)
1083 [ + + ]: 44532 : ortho_axes(:,22)=(/1,-1,2/)
1084 : :
1085 : 11133 : hexa_axes(:,:)=0
1086 : 11133 : hexa_axes(1,1)=1
1087 : 11133 : hexa_axes(2,2)=1
1088 : 11133 : hexa_axes(3,3)=1
1089 [ + + ]: 44532 : hexa_axes(:,4)=(/1,-1,0/)
1090 [ + + ]: 44532 : hexa_axes(:,5)=(/2,1,0/)
1091 [ + + ]: 44532 : hexa_axes(:,6)=(/1,1,0/)
1092 [ + + ]: 44532 : hexa_axes(:,7)=(/1,2,0/)
1093 : :
1094 : : !Determine the point group from the list of symmetry operations.
1095 : : !Also determine the holohedry, up to one undeterminacy : hR versus hP
1096 : 11133 : call symptgroup(iholohedry,nsym,ptgroup,symrel)
1097 : :
1098 : : !DEBUG
1099 : : !write(std_out,*)' symbrav, after symptgroup: nsym=',nsym
1100 : : !call flush(std_out)
1101 : : !write(std_out,*)' symbrav: symrel='
1102 : : !do isym=1,nsym
1103 : : ! write(std_out,'(9i4)')symrel(:,:,isym)
1104 : : !enddo
1105 : : !write(std_out,*)' symbrav: iholohedry=',iholohedry
1106 : : !call flush(std_out)
1107 : : !ENDDEBUG
1108 : :
1109 : : !Loop over trial deformations
1110 : : !This is needed in case the Bravais lattice determination from the lattice vectors
1111 : : !has a higher holohedry than the real one, in which the symmetry
1112 : : !operations for the atoms (or electric field, etc) are taken into account
1113 : 11133 : iaxis=0
1114 : 11133 : invariant=0
1115 : 11133 : next_stage=0
1116 : 11133 : rprimdnow(:,:)=rprimd(:,:)
1117 : 11133 : rprimdtry(:,:)=rprimd(:,:)
1118 [ - + + - : 33399 : ABI_MALLOC(symrelconv,(3,3,nsym))
+ - + - -
+ + - ]
1119 : :
1120 : : !At most will have to try naxes_ortho*5 deformations (naxes_ortho axes, five stages)
1121 : : !First, test whether the current recognition of Bravais lattice is problematic (iholohedry differs from bravais(1)).
1122 : : !Then, if there is a problem, test different deformations of rprimd, one after the other.
1123 : : !For each try, there is a new rprimdtry from the different set of deformation, that is generated later in the loop
1124 : : !Also bravais1now and rprimdnow might changed (and progressively lowered).
1125 : : !The latter change induces at most 5 stages for the computation (cubic->tetragonal->orthorhombic->monoclinic->triclinic).
1126 : : !After an upgrade of bravais1now and rprimdnow, one has to restart the full set of deformations.
1127 : : !Not sure that this procedure resolves all cases, but seems to work on >40000 inaccurate POSCAR files.
1128 [ + - ]: 15629 : do ideform=1,naxes_ortho*5
1129 : :
1130 : 15629 : maxsym=max(192,msym)
1131 [ - + + - : 46887 : ABI_MALLOC(ptsymrel,(3,3,maxsym))
+ - - + -
+ ]
1132 : 15629 : call symlatt(bravais,std_out,maxsym,nptsym,ptsymrel,rprimdtry,tolsym)
1133 [ - + ]: 15629 : ABI_FREE(ptsymrel)
1134 : :
1135 : : !DEBUG
1136 : : !write(std_out,*)' symbrav: inside loop with ideform,iaxis=',ideform,iaxis
1137 : : !write(std_out,'(a,9f12.6)')' rprimdtry=',rprimdtry(:,:)
1138 : : !write(std_out,'(a,2i4)')' bravais(1:2)=',bravais(1:2)
1139 : : !call flush(std_out)
1140 : : !ENDDEBUG
1141 : :
1142 : :
1143 : : ! Examine the agreement with bravais(1)
1144 : : ! Warning : might change Bravais lattice hR to hP, if hexagonal axes
1145 : 15629 : problem=0
1146 : 23625 : select case (bravais(1))
1147 : : case (7)
1148 [ + + ]: 7996 : if(iholohedry<6)problem=1
1149 [ - + ]: 7996 : if(iholohedry==6)problem=2
1150 : : case (6)
1151 [ + + ]: 669 : if(iholohedry<4)problem=1
1152 [ - + ]: 669 : if(iholohedry==7 .or. iholohedry==4)problem=2
1153 : : ! Here, change hR into hP
1154 [ + + ]: 669 : if(iholohedry==5)iholohedry=6
1155 : : case (5)
1156 [ + + ]: 806 : if(iholohedry<4)problem=1
1157 [ - + ]: 806 : if(iholohedry==7 .or. iholohedry==6 .or. iholohedry==4)problem=2
1158 : : case (4)
1159 [ + + ]: 2050 : if(iholohedry<4)problem=1
1160 [ - + ]: 2050 : if(iholohedry>4)problem=2
1161 : : case (3)
1162 [ + + ]: 898 : if(iholohedry<3)problem=1
1163 [ - + ]: 898 : if(iholohedry>3)problem=2
1164 : : case (2)
1165 [ + + ]: 1771 : if(iholohedry<2)problem=1
1166 [ - + ]: 1771 : if(iholohedry>2)problem=2
1167 : : case (1)
1168 [ + + + + : 15629 : if(iholohedry>1)problem=2
+ + + - -
+ ]
1169 : : end select
1170 : :
1171 : : ! This is the usual situation, in which the lattice belong to the same holohedry
1172 : : ! as the lattice+atoms (+electric field + ...)
1173 [ + + ]: 15629 : if(problem==0)exit
1174 : :
1175 [ - + ]: 4496 : if(problem==2)then
1176 [ # # ]: 0 : if(iaxis==0)then
1177 : : write(msg, '(3a,i3,3a,i3,7a)' )&
1178 : 0 : 'The Bravais lattice determined only from the primitive',ch10,&
1179 : 0 : 'vectors (rprim or angdeg), bravais(1)=',bravais(1),', is not compatible',ch10,&
1180 : 0 : 'with the real one, iholohedry=',iholohedry,', obtained by taking into',ch10,&
1181 : 0 : 'account the symmetry operations. This might be due to an insufficient',ch10,&
1182 : 0 : 'number of digits in the specification of rprim (at least 10),',ch10,&
1183 : 0 : 'or to an erroneous rprim or angdeg. If this is not the case, then ...'
1184 : 0 : ABI_BUG(msg)
1185 : : end if
1186 [ # # ]: 0 : if(iaxis==1)then
1187 : : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1188 : 0 : 'Could not succeed to determine the bravais lattice',ch10,&
1189 : 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1190 : 0 : 'bravais(1)=',bravais(1),ch10,&
1191 : 0 : 'iholohedry=',iholohedry
1192 : 0 : ABI_BUG(msg)
1193 : : end if
1194 : : ! Try to increase tolsym to find the Bravais lattice.
1195 : : maxsym=max(192,msym)
1196 [ # # # # : 0 : ABI_MALLOC(ptsymrel,(3,3,maxsym))
# # # # ]
1197 : : !DEBUG
1198 : : ! write(6,*)' symbrav : will call symlatt, 3*tolsym=',3*tolsym
1199 : : !ENDDEBUG
1200 : 0 : call symlatt(bravais,std_out,maxsym,nptsym,ptsymrel,rprimdtry,3*tolsym)
1201 [ # # ]: 0 : ABI_FREE(ptsymrel)
1202 [ # # ]: 0 : if(bravais(1)==iholohedry)then
1203 : : ! Succeeded
1204 : : exit
1205 : : else
1206 : : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1207 : 0 : 'Could not succeed to determine the bravais lattice, even after considering a larger tolsym',ch10,&
1208 : 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1209 : 0 : 'bravais(1)=',bravais(1),ch10,&
1210 : 0 : 'iholohedry=',iholohedry
1211 : 0 : ABI_BUG(msg)
1212 : : end if
1213 : : end if
1214 : :
1215 [ + - ]: 4496 : if(problem==1)then ! One is left with the problem=1 case, basically iholohedry is lower than bravais(1)
1216 [ + + ]: 4496 : if(iaxis==0)then
1217 : : write(msg, '(a,a,a,i3,a,a,a,i3,a,a,a)' )&
1218 : 2934 : 'The Bravais lattice determined only from the primitive',ch10,&
1219 : 2934 : 'vectors, bravais(1)=',bravais(1),', is more symmetric',ch10,&
1220 : 2934 : 'than the real one, iholohedry=',iholohedry,', obtained by taking into',ch10,&
1221 : 5868 : 'account the atomic positions. Start deforming the primitive vector set.'
1222 : 2934 : ABI_COMMENT(msg)
1223 : : next_stage=1
1224 : : else if(iaxis/=0)then
1225 [ + + ]: 1562 : if(bravais(1)<bravais1now)then
1226 : : write(msg, '(3a,i3,3a,i3,2a)' )&
1227 : 1347 : 'The Bravais lattice determined from modified primitive',ch10,&
1228 : 1347 : 'vectors, bravais(1)=',bravais(1),', has a lower symmetry than before,',ch10,&
1229 : 1347 : 'but is still more symmetric than the real one, iholohedry=',iholohedry,ch10,&
1230 : 2694 : 'obtained by taking into account the atomic positions.'
1231 : 1347 : ABI_COMMENT(msg)
1232 : : next_stage=1
1233 [ - + ]: 215 : else if(iaxis==1)then
1234 : : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1235 : 0 : 'Could not succeed to determine the bravais lattice',ch10,&
1236 : 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1237 : 0 : 'bravais(1)=',bravais(1),ch10,&
1238 : 0 : 'iholohedry=',iholohedry
1239 : 0 : ABI_BUG(msg)
1240 : : end if
1241 : : end if
1242 : : end if ! problem==1
1243 : :
1244 : : ! One is here when problem=1 (iholohedry < bravais(1)) and either
1245 : : ! - iaxis=0 (no deformation has been tried yet),
1246 : : ! - some deformation iaxis has been tried giving bravais(1), but iholohedry < bravais(1) < bravais1now
1247 : : ! - some deformation iaxis has been tried giving bravais(1), but iholohedry < bravais(1) = bravais1now and iaxis/=1 .
1248 : : ! Also, note that next_stage is still 0 when bravais(1)=bravais1now .
1249 : : ! The loop has been ended (so, the search failed) when bravais(1)=bravais1now and iaxis==1.
1250 : :
1251 : : if(next_stage==1)then
1252 : 4281 : bravais1now=bravais(1)
1253 : 4281 : rprimdnow(:,:)=rprimdtry(:,:)
1254 : : ! Generate the symmetry operations in the conventional vector coordinates
1255 [ + + ]: 17124 : rprimdconv(:,1)=bravais(3:5)
1256 [ + + ]: 17124 : rprimdconv(:,2)=bravais(6:8)
1257 [ + + ]: 17124 : rprimdconv(:,3)=bravais(9:11)
1258 : 4281 : axes(:,:)=zero
1259 : 4281 : axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one
1260 [ + + + + : 217143 : symrelconv(:,:,1:nsym)=symrel(:,:,1:nsym)
+ + ]
1261 : 4281 : call symrelrot(nsym,rprimdconv,axes,symrelconv,tolsym)
1262 [ + + ]: 4281 : if(bravais(1)/=6)then
1263 : 4264 : iaxis=naxes_ortho+1
1264 : : else
1265 : 17 : iaxis=naxes_hexa+1
1266 : : end if
1267 : : next_stage=0
1268 : : !DEBUG
1269 : : ! write(std_out,*)' symbrav: next stage, bravais(1), bravais(2) and symrelconv'
1270 : : ! write(std_out,'(a,2i4)')' bravais(1:2)=',bravais(1:2)
1271 : : ! write(std_out,'(a,9f12.6)')' rprimdconv=',rprimdconv(:,:)
1272 : : ! do isym=1,nsym
1273 : : ! write(std_out,'(9i4)')symrelconv(:,:,isym)
1274 : : ! enddo
1275 : : ! call flush(std_out)
1276 : : !ENDDEBUG
1277 : : end if
1278 : :
1279 : : ! Go to the next iaxis that will be left invariant
1280 : 4496 : iaxis=iaxis-1
1281 [ + - ]: 25112 : do jaxis=iaxis,1,-1
1282 [ + + ]: 25112 : if(bravais(1)/=6)then
1283 [ + + ]: 100372 : axis_trial(:)=ortho_axes(:,jaxis)
1284 : : else
1285 [ + + ]: 76 : axis_trial(:)=hexa_axes(:,jaxis)
1286 : : end if
1287 : : ! DEBUG
1288 : : ! write(std_out,*)' symbrav : ixaxis, trial jaxis=',iaxis,jaxis
1289 : : ! write(std_out,*)' axis_trial=',axis_trial
1290 : : ! ENDDEBUG
1291 : 25112 : invariant=1
1292 : : ! Examine whether all symmetry operations leave the axis invariant (might be reversed, though)
1293 [ + + ]: 231370 : do isym=1,nsym
1294 [ - + + + : 1856322 : if(sum(abs(matmul(symrelconv(:,:,isym),axis_trial)+(-axis_trial(:))))/=0 .and. &
+ + + + ]
1295 [ + + + + : 7244142 : sum(abs(matmul(symrelconv(:,:,isym),axis_trial)+axis_trial(:)))/=0 )invariant=0
+ + + + +
+ + + + -
+ - ]
1296 : : end do
1297 [ + + ]: 25112 : if(invariant==1)then
1298 : 4496 : iaxis=jaxis
1299 : : !DEBUG
1300 : : ! write(msg, '(2a,i3)' )ch10,' symbrav : found invariant axis, jaxis=',jaxis
1301 : : ! call wrtout(std_out,msg)
1302 : : !ENDDEBUG
1303 : 4496 : exit
1304 : : end if
1305 : : end do
1306 : :
1307 [ - + ]: 4496 : if(invariant==0)then
1308 : : ! Not a single axis was invariant with respect to all operations ?!
1309 : : ! do isym=1,nsym; write(std_out, '(a,10i4)' )' isym,symrelconv=',isym,symrelconv(:,:,isym); enddo
1310 : : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1311 : 0 : 'Could not succeed to determine the bravais lattice (not a single invariant)',ch10,&
1312 : 0 : 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,&
1313 : 0 : 'bravais(1)=',bravais(1),ch10,&
1314 : 0 : 'iholohedry=',iholohedry
1315 : 0 : ABI_BUG(msg)
1316 : : end if
1317 : :
1318 : 4496 : call matr3inv(rprimdconv,rprimdconv_invt)
1319 : : axis_red(:)=axis_trial(1)*rprimdconv_invt(1,:)+ &
1320 : : & axis_trial(2)*rprimdconv_invt(2,:)+ &
1321 [ + + ]: 17984 : & axis_trial(3)*rprimdconv_invt(3,:)
1322 : : axis_cart(:)=axis_red(1)*rprimdnow(:,1)+ &
1323 : : & axis_red(2)*rprimdnow(:,2)+ &
1324 [ + + ]: 17984 : & axis_red(3)*rprimdnow(:,3)
1325 [ + + ]: 17984 : norm=sum(axis_cart(:)**2)
1326 : : !DEBUG
1327 : : ! write(6,*)' axis_trial =',axis_trial
1328 : : ! write(6,*)' axis_red =',axis_red
1329 : : ! write(6,*)' axis_cart =',axis_cart
1330 : : ! write(6,*)' rprimdnow=',rprimdnow
1331 : : !ENDDEBUG
1332 : : ! Expand by a uniform, quite arbitrary, dilatation, along the invariant axis
1333 : : ! Note : make these dilatation different, according to ideform
1334 : : ! XG 20151221 : Still, the interplay between the size of the deformation and the tolsym is not easy to address.
1335 : : ! Indeed the deformation must be sufficiently large to be perceived by symlatt as a real breaking of the
1336 : : ! symmetry of the lattice. In order to deal with all the small values od tolsym, it has been set at a minimum of tol3,
1337 : : ! but it must also be larger than tolsym. Moreover, for some axis choice, the deformation is not aligned with the axis, decreasing
1338 : : ! the effective deformation length. An additional factor of three is thus included, actually increased to six just to be sure...
1339 [ + + ]: 44746 : do ii=1,3
1340 : 13488 : scprod=axis_cart(1)*rprimdnow(1,ii)+axis_cart(2)*rprimdnow(2,ii)+axis_cart(3)*rprimdnow(3,ii)
1341 [ + + ]: 58448 : rprimdtry(:,ii)=rprimdnow(:,ii)+ideform*(max(tol3,six*tolsym)-tol6)*scprod/norm*axis_cart(:)
1342 : : end do
1343 : :
1344 : : !DEBUG
1345 : : ! write(6,*)' rprimdtry=',rprimdtry
1346 : : !ENDDEBUG
1347 : :
1348 : : end do ! ideform
1349 : :
1350 [ - + ]: 11133 : if(bravais(1)/=iholohedry)then
1351 : : write(msg, '(3a,3i3,2a,i3,2a,i3)' )&
1352 : 0 : 'Despite efforts, Could not succeed to determine the bravais lattice :',ch10,&
1353 : 0 : 'bravais(1)=',bravais(1),ch10,&
1354 : 0 : 'iholohedry=',iholohedry
1355 : 0 : ABI_BUG(msg)
1356 : : end if
1357 : :
1358 [ - + ]: 11133 : ABI_FREE(symrelconv)
1359 : :
1360 [ - + ]: 11133 : if (PRESENT(axis)) then ! Return symmetry axis.
1361 : 0 : axis=(/0,0,0/)
1362 [ # # ]: 0 : if (iaxis/=0) then
1363 [ # # ]: 0 : if(bravais(1)/=6)then
1364 [ # # ]: 0 : axis=ortho_axes(:,iaxis)
1365 : : else
1366 [ # # ]: 0 : axis=hexa_axes(:,iaxis)
1367 : : end if
1368 : : end if
1369 : : end if
1370 : :
1371 : : !DEBUG
1372 : : !write(std_out,'(a)')' symbrav : exit '
1373 : : !ENDDEBUG
1374 : :
1375 : 11133 : end subroutine symbrav
1376 : : !!***
1377 : :
1378 : : !!****f* m_symfind/symspgr
1379 : : !! NAME
1380 : : !! symspgr
1381 : : !!
1382 : : !! FUNCTION
1383 : : !! Find the type of each symmetry operation (calling symcharac):
1384 : : !! proper symmetries 1,2,2_1,3,3_1,3_2,4,4_1,4_2,4_3,6,6_1,...6_5
1385 : : !! improper symmetries -1,m,a,b,c,d,n,g,-3,-4,-6 ,
1386 : : !! Then, build an array with the number of such operations.
1387 : : !! Then, call symlist to identify the space group.
1388 : : !!
1389 : : !! INPUTS
1390 : : !! bravais(11): bravais(1)=iholohedry
1391 : : !! bravais(2)=center
1392 : : !! bravais(3:11)=coordinates of rprimd in the axes
1393 : : !! of the conventional bravais lattice (*2 if center/=0)
1394 : : !! nsym=actual number of symmetries
1395 : : !! symrel(3,3,nsym)= nsym symmetry operations in real space in terms
1396 : : !! of primitive translations
1397 : : !! tnons(3,nsym)=nonsymmorphic translations for each symmetry (would
1398 : : !! be 0 0 0 each for a symmorphic space group)
1399 : : !!
1400 : : !! OUTPUT
1401 : : !! labels(maxsym=192)= labels of the symmetry operations
1402 : : !! spgroup=symmetry space group number
1403 : : !!
1404 : : !! NOTES
1405 : : !! It is assumed that the symmetry operations will be entered in the
1406 : : !! symrel tnons arrays, for the PRIMITIVE cell. The matrix of transformation
1407 : : !! from the primitive cell to the conventional cell is described
1408 : : !! in the array "bravais" (see symlatt.F90).
1409 : : !! The present routine first make the transformation from the
1410 : : !! primitive coordinates to the conventional ones, then eventually
1411 : : !! generate additional symmetries, taking into account the
1412 : : !! centering translations.
1413 : : !! Then, the order and determinant of each symmetry operation
1414 : : !! is determined.
1415 : : !!
1416 : : !! For proper symmetries (rotations), the
1417 : : !! associated translation is also determined.
1418 : : !! However, left or right handed screw rotations are
1419 : : !! not (presently) distinguished, and will be attributed equally
1420 : : !! to left or right.
1421 : : !!
1422 : : !! For the detailed description of the labelling of the axes, see symaxes.f and symplanes.f
1423 : : !!
1424 : : !! SOURCE
1425 : :
1426 : 11744 : subroutine symspgr(bravais,labels,nsym,spgroup,symrel,tnons,tolsym)
1427 : :
1428 : : use m_numeric_tools, only : OPERATOR(.x.)
1429 : :
1430 : : !Arguments ------------------------------------
1431 : : !scalars
1432 : : integer,intent(in) :: nsym
1433 : : integer,intent(out) :: spgroup
1434 : : real(dp),intent(in) :: tolsym
1435 : : !arrays
1436 : : integer,intent(in) :: bravais(11),symrel(3,3,nsym)
1437 : : real(dp),intent(in) :: tnons(3,nsym)
1438 : : character(len=128),intent(out) :: labels(192) ! 192 = maxsym
1439 : :
1440 : : !Local variables-------------------------------
1441 : : !scalars
1442 : : ! logical,parameter :: verbose=.FALSE.
1443 : : integer :: additional_info,brvltt,center,direction=0,found,iholohedry,ii
1444 : : integer :: ishift,isym,jj,nshift,nsymconv,spgaxor,spgorig,sporder
1445 : : character(len=1) :: brvsb
1446 : : character(len=15) :: intsb,ptintsb,ptschsb,schsb
1447 : : character(len=35) :: intsbl
1448 : : character(len=500) :: msg
1449 : : !arrays
1450 : : integer :: ivec1(3), ivec2(3)
1451 : : integer :: n_axes(31),n_axest(31),prime(5),test_direction(3),symrel_uni(3,3)
1452 : : integer :: uniaxis(3),uniaxis_try(3)
1453 : 11744 : integer,allocatable :: determinant(:),symrelconv(:,:,:),t_axes(:)
1454 : : real(dp) :: axes(3,3),rprimdconv(3,3),vect(3,3)
1455 : 11744 : real(dp),allocatable :: shift(:,:),tnonsconv(:,:)
1456 : : !**************************************************************************
1457 : :
1458 : : DBG_ENTER("COLL")
1459 : :
1460 : : !Initialize brvltt, from bravais(2) and bravais(1)
1461 : 11744 : center=bravais(2)
1462 : 11744 : iholohedry=bravais(1)
1463 : 11744 : brvltt=1
1464 [ + + ]: 11744 : if(center==-1)brvltt=2 ! Inner centering
1465 [ + + ]: 11744 : if(center==-3)brvltt=3 ! Face centering
1466 [ + + ]: 11744 : if(center==1)brvltt=5 ! A-Face centering
1467 [ + + ]: 11744 : if(center==2)brvltt=6 ! B-Face centering
1468 [ + + ]: 11744 : if(center==3)brvltt=4 ! C-Face centering
1469 [ + + ]: 11744 : if(iholohedry==5)brvltt=7 ! Rhombohedral
1470 : :
1471 : : !Produce the symmetry operations, in the axis of the conventional cell
1472 : 11744 : nsymconv=nsym
1473 [ + + ]: 11744 : if(center/=0)nsymconv=2*nsymconv
1474 [ + + ]: 11744 : if(center==-3)nsymconv=4*nsym
1475 [ - + + - : 35232 : ABI_MALLOC(symrelconv,(3,3,nsymconv))
+ - + - -
+ - + ]
1476 [ - + + - : 35232 : ABI_MALLOC(tnonsconv,(3,nsymconv))
+ - + - -
+ - + ]
1477 : :
1478 : : !Produce symrel and tnons in conventional axes,
1479 : : !name them symrelconv and tnonsconv
1480 [ + + ]: 46976 : rprimdconv(:,1)=bravais(3:5)
1481 [ + + ]: 46976 : rprimdconv(:,2)=bravais(6:8)
1482 [ + + ]: 46976 : rprimdconv(:,3)=bravais(9:11)
1483 : :
1484 [ + + + + : 79856 : if(center/=0)rprimdconv(:,:)=rprimdconv(:,:)*half
+ + ]
1485 : :
1486 : 11744 : axes(:,:)=zero
1487 : 11744 : axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one
1488 [ + + + + : 3783031 : symrelconv(:,:,1:nsym)=symrel(:,:,1:nsym)
+ + ]
1489 : : !Note that the number of symmetry operations is still nsym
1490 : 11744 : call symrelrot(nsym,rprimdconv,axes,symrelconv,tolsym)
1491 : :
1492 : 11744 : call xred2xcart(nsym,rprimdconv,tnonsconv,tnons)
1493 : : !Gives the associated translation, with components in the interval ]-0.5,0.5] .
1494 [ + + + + ]: 1172140 : tnonsconv(:,1:nsym)=tnonsconv(:,1:nsym)-nint(tnonsconv(:,1:nsym)-tol6)
1495 : :
1496 : : !If the Bravais lattice is centered, duplicate or quadruplicate
1497 : : !the number of symmetry operations, using the Bravais lattice shifts
1498 : 11744 : nshift=1
1499 [ + + ]: 11744 : if(center/=0)nshift=2
1500 [ + + ]: 11744 : if(center==-3)nshift=4
1501 [ - + + - : 23488 : ABI_MALLOC(shift,(3,nshift))
- + - + ]
1502 [ + + ]: 46976 : shift(:,1)=zero
1503 [ + + ]: 11744 : if(center/=0 .and. center/=-3)then
1504 [ + + ]: 5204 : shift(:,2)=half
1505 [ + + ]: 1301 : if(center==1)shift(1,2)=zero
1506 [ + + ]: 1301 : if(center==2)shift(2,2)=zero
1507 [ + + ]: 1301 : if(center==3)shift(3,2)=zero
1508 [ + + ]: 10443 : else if(center==-3)then
1509 [ + + ]: 17500 : shift(:,2)=half ; shift(1,2)=zero
1510 [ + + ]: 17500 : shift(:,3)=half ; shift(2,3)=zero
1511 [ + + ]: 17500 : shift(:,4)=half ; shift(3,4)=zero
1512 : : end if ! center/=0 or -3
1513 [ + + ]: 11744 : if(nshift/=1)then
1514 [ + + ]: 20102 : do ishift=2,nshift
1515 [ + - + + : 14897658 : symrelconv(:,:,(ishift-1)*nsym+1:ishift*nsym)=symrelconv(:,:,1:nsym)
+ + + + +
+ + + +
+ ]
1516 [ + + ]: 592534 : do isym=1,nsym
1517 [ + + ]: 2304154 : tnonsconv(:,(ishift-1)*nsym+isym)=tnonsconv(:,isym)+shift(:,ishift)
1518 : : end do
1519 : : end do ! ishift
1520 : : end if ! nshift/=1
1521 : :
1522 : : !At this stage, all the symmetry operations are available,
1523 : : !expressed in the conventional axis, and also include
1524 : : !the Bravais lattive translations, and associated operations...
1525 : :
1526 : 11744 : n_axes(:)=0
1527 : :
1528 [ - + + - : 35232 : ABI_MALLOC(determinant,(nsymconv))
+ - + - -
+ - + ]
1529 : :
1530 : : !Get the determinant
1531 : 11744 : call symdet(determinant,nsymconv,symrelconv)
1532 : :
1533 : : !Get the order of each the symmetry operation, as well as the maximal order
1534 : : !Also, examine whether each symmetry operation is the inversion, or a root of the inversion (like -3)
1535 : : !Decide which kind of point symmetry operation it is
1536 : : !Finally assign tnonsconv order and decide the space symmetry operation
1537 : :
1538 [ - + + - : 23488 : ABI_MALLOC(t_axes,(nsymconv))
- + - + ]
1539 : :
1540 [ + + ]: 874275 : do isym=1,nsymconv
1541 : :
1542 : : ! Note : nsymconv might be bigger than 192, but only for non-primitive cells, in which case labels will not be echoed anywhere.
1543 : : ! 192 is the fixed dimension of labels, so this avoids possible memory problems.
1544 : : call symcharac(center, determinant(isym), iholohedry, isym, labels(mod(isym-1,192)+1), &
1545 : 862531 : symrelconv(:,:,isym), tnonsconv(:,isym), t_axes(isym))
1546 [ - + ]: 862531 : if (t_axes(isym) == -1) then
1547 : 0 : write(msg, '(a,a,i3,a,3(a,3i4,a),a,3es22.12,a,a,3es22.12)' )ch10,&
1548 : 0 : ' symspgr: problem with isym=',isym,ch10,&
1549 : 0 : ' symrelconv(:,1,isym)=',symrelconv(:,1,isym),ch10,&
1550 : 0 : ' symrelconv(:,2,isym)=',symrelconv(:,2,isym),ch10,&
1551 : 0 : ' symrelconv(:,3,isym)=',symrelconv(:,3,isym),ch10,&
1552 : 0 : ' tnonsconv(:,isym)=',tnonsconv(:,isym)
1553 : 0 : call wrtout(std_out,msg)
1554 : 0 : write(msg, '(a,i0,2a)' )'The space symmetry operation number',isym,ch10,'is not a (translated) root of unity'
1555 : 0 : ABI_BUG(msg)
1556 [ - + ]: 862531 : else if (t_axes(isym) == -2) then
1557 : 0 : write(msg, '(a,i0,a)' )'The symmetry operation number ',isym,' is not a root of unity'
1558 : 0 : ABI_BUG(msg)
1559 : : end if
1560 : :
1561 : 874275 : n_axes(t_axes(isym))=n_axes(t_axes(isym))+1
1562 : : end do ! isym=1,nsymconv
1563 : :
1564 [ + + - + ]: 375808 : if (sum(n_axes)-nsymconv/=0) then
1565 : : write(msg, '(7a)' )&
1566 : 0 : 'Not all the symmetries have been recognized. ',ch10,&
1567 : 0 : 'This might be due either to an error in the input file',ch10,&
1568 : 0 : 'or to a BUG in ABINIT',ch10,&
1569 : 0 : 'Please contact the ABINIT group.'
1570 : 0 : ABI_WARNING(msg)
1571 : : end if
1572 : :
1573 : : !write(std_out,*)' symspgr : brvltt,nsymconv=',brvltt,nsymconv
1574 : : !write(std_out,*)' n_axes(1:10)=',n_axes(1:10)
1575 : : !write(std_out,*)' n_axes(11:20)=',n_axes(11:20)
1576 : : !write(std_out,*)' n_axes(21:31)=',n_axes(21:31)
1577 : :
1578 : : !Treat cases in which the space group cannot be identified on the
1579 : : !basis of n_axes one need additional information
1580 [ + + ]: 11744 : if(brvltt==1)then
1581 : : ! If the bravais lattice is primitive
1582 [ + + ]: 5256 : if(nsymconv==4)then
1583 : 307 : 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/)
1584 [ + + + + ]: 9824 : if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 27 (Pcc2) or 32 (Pba2)
1585 : 16 : write(std_out,*)' symspgr: 27 or 32'
1586 : 16 : additional_info=2
1587 : : ! Select binary axis
1588 [ + + ]: 80 : do isym=1,nsymconv
1589 [ + + ]: 80 : if(t_axes(isym)==8)then
1590 : : ! Find direction of binary axis
1591 [ + - ]: 16 : if(symrelconv(1,1,isym)==1)direction=1
1592 [ + - ]: 16 : if(symrelconv(2,2,isym)==1)direction=2
1593 [ + - ]: 16 : if(symrelconv(3,3,isym)==1)direction=3
1594 : : end if
1595 : : end do
1596 : : ! Examine the projection of the translation vector of the a, b or c mirror planes onto the binary axis
1597 [ + + ]: 80 : do isym=1,nsymconv
1598 [ + + ]: 80 : if(t_axes(isym)==16)then
1599 [ + + ]: 32 : if(abs(tnonsconv(direction,isym))>tol8)additional_info=1
1600 : : end if
1601 : : end do
1602 : : end if
1603 [ + + ]: 4949 : else if(nsymconv==8)then
1604 : 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/)
1605 [ + + + + ]: 28192 : if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 55 (Pbam) or 57 (Pbcm)
1606 : 26 : write(std_out,*)' symspgr: 55 or 57'
1607 : 26 : additional_info=1
1608 : : ! Select mirror plane m
1609 [ + + ]: 234 : do isym=1,nsymconv
1610 [ + + ]: 234 : if(t_axes(isym)==15)then
1611 : : ! Find direction of mirror plane
1612 [ - + ]: 26 : if(symrelconv(1,1,isym)==-1)direction=1
1613 [ - + ]: 26 : if(symrelconv(2,2,isym)==-1)direction=2
1614 [ + - ]: 26 : if(symrelconv(3,3,isym)==-1)direction=3
1615 : : end if
1616 : : end do
1617 : : ! Examine the projection of the translation vector of the a, b, or c mirror planes
1618 : : ! onto the binary axis
1619 [ + + ]: 234 : do isym=1,nsymconv
1620 [ + + ]: 234 : if(t_axes(isym)==16)then
1621 [ + + ]: 52 : if(abs(tnonsconv(direction,isym))>tol8)additional_info=2
1622 : : end if
1623 : : end do
1624 : : end if
1625 : 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/)
1626 [ + + + + ]: 28192 : if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 56 (Pccn) or 60 (Pbcn)
1627 : 26 : write(std_out,*)' symspgr: 56 or 60'
1628 : 26 : additional_info=1
1629 : : ! Select mirror plane n
1630 [ + + ]: 234 : do isym=1,nsymconv
1631 [ + + ]: 234 : if(t_axes(isym)==18)then
1632 : : ! Find direction of mirror plane
1633 [ - + ]: 26 : if(symrelconv(1,1,isym)==-1)direction=1
1634 [ - + ]: 26 : if(symrelconv(2,2,isym)==-1)direction=2
1635 [ + - ]: 26 : if(symrelconv(3,3,isym)==-1)direction=3
1636 : : end if
1637 : : end do
1638 : : ! Examine the projection of the translation vector of the a, b, or c mirror planes
1639 : : ! onto the binary axis
1640 [ + + ]: 234 : do isym=1,nsymconv
1641 [ + + ]: 234 : if(t_axes(isym)==16)then
1642 [ + + ]: 52 : if(abs(tnonsconv(direction,isym))<tol8)additional_info=2
1643 : : end if
1644 : : end do
1645 : : end if
1646 : : end if
1647 [ + + ]: 6488 : else if(brvltt==2)then
1648 : : ! In the few next lines, use additional_info as a flag
1649 : 706 : additional_info=0
1650 : : ! If the bravais lattice is inner-centered
1651 [ + + ]: 706 : if(nsymconv==8)then
1652 : : ! Test spgroup 23 (I222) or 24 (I2_{1}2_{1}2_{1})
1653 : 113 : 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/)
1654 [ + + + + ]: 3616 : if(sum((n_axes-n_axest)**2)==0) additional_info=1
1655 [ + + ]: 593 : else if(nsymconv==24)then
1656 : : ! Test spgroup 197 (I23) or 199 (I2_{1}3)
1657 : 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/)
1658 [ + + + - ]: 64 : if(sum((n_axes-n_axest)**2)==0) additional_info=1
1659 : : end if
1660 [ + + ]: 706 : if(additional_info==1)then
1661 : 8 : write(std_out,*)' symspgr: (23 or 24) or (197 or 199)'
1662 : : ! Select the three binary axes (they might be 2 or 2_1 !)
1663 : 8 : test_direction(:)=0
1664 [ + + ]: 104 : do isym=1,nsymconv
1665 [ + + ]: 104 : if(t_axes(isym)==20)then
1666 : : ! Find direction of axis
1667 [ + + ]: 96 : do direction=1,3
1668 [ + + ]: 96 : if(symrelconv(direction,direction,isym)==1)then
1669 : 24 : test_direction(direction)=1
1670 [ - + ]: 24 : if(abs(tnonsconv(direction,isym))<tol8)then
1671 [ # # ]: 0 : vect(:,direction)=tnonsconv(:,isym)
1672 : : else
1673 [ + + ]: 96 : vect(:,direction)=tnonsconv(:,isym)+half
1674 : : end if
1675 [ + + ]: 96 : vect(:,direction)=vect(:,direction)-nint(vect(:,direction)-tol8)
1676 : 24 : vect(direction,direction)=zero
1677 : : end if
1678 : : end do ! direction=1,3
1679 : : end if ! if binary axis
1680 : : end do ! isym
1681 [ + - - + : 8 : if(test_direction(1)/=1 .or. test_direction(2)/=1 .and. test_direction(3)/=1)then
- - ]
1682 : : write(msg, '(5a,3i4)' )&
1683 : 0 : 'For space groups 23, 24, 197 or 197, the three binary axes',ch10,&
1684 : 0 : 'are not equally partitioned along the x, y and z directions',ch10,&
1685 : 0 : 'test_direction(1:3)=',test_direction(:)
1686 : 0 : ABI_BUG(msg)
1687 : : end if
1688 : 8 : additional_info=1
1689 [ + + + - : 8 : if(abs(vect(1,2)-vect(1,3))>tol8 .or. abs(vect(2,1)-vect(2,3))>tol8 .or. &
- + ]
1690 : 4 : abs(vect(3,1)-vect(3,2))>tol8) additional_info=2
1691 : : end if ! additional information are needed
1692 : : end if ! brvltt==1
1693 : :
1694 [ + + ]: 11744 : if (brvltt==0 .or. brvltt==1) then ! Primitive
1695 : 5256 : call symlist_prim(additional_info,nsymconv,n_axes,spgroup)
1696 [ + + ]: 6488 : else if(brvltt==2)then
1697 : 706 : call symlist_bcc(additional_info,nsymconv,n_axes,spgroup)
1698 [ + + ]: 5782 : else if(brvltt==3)then
1699 : 4375 : call symlist_fcc(nsymconv,n_axes,spgroup)
1700 : : else
1701 : 1407 : call symlist_others(brvltt,nsymconv,n_axes,spgroup)
1702 : : end if
1703 : :
1704 [ + + ]: 11744 : if(spgroup==0) then
1705 : : write(msg, '(5a)' )&
1706 : 17 : 'Could not find the space group.',ch10,&
1707 : 17 : 'This often happens when the user selects a restricted set of symmetries ',ch10,&
1708 : 34 : 'in the input file, instead of letting the code automatically find symmetries.'
1709 : 17 : ABI_WARNING(msg)
1710 : : end if
1711 : :
1712 : 11744 : spgorig=1 ; spgaxor=1
1713 : 11744 : call spgdata(brvsb,intsb,intsbl,ptintsb,ptschsb,schsb,spgaxor,spgroup,sporder,spgorig)
1714 : :
1715 [ + + ]: 11744 : if(spgroup/=0)then
1716 : 11727 : write(msg, '(a,i4,2x,a,a,a,a,a)' ) ' symspgr: spgroup=',spgroup,trim(brvsb),trim(intsb),' (=',trim(schsb),')'
1717 : 11727 : call wrtout(std_out,msg)
1718 : : end if
1719 : :
1720 [ + + ]: 11744 : if(bravais(1)==7)then
1721 : 5627 : write(msg, '(a)' ) ' symspgr: optical characteristics = isotropic '
1722 : 5627 : call wrtout(std_out,msg)
1723 [ + + ]: 6117 : else if(bravais(1)==4 .or. bravais(1)==5 .or. bravais(1)==6)then
1724 : 3068 : write(msg, '(a)' ) ' symspgr: optical characteristics = uniaxial '
1725 : 3068 : call wrtout(std_out,msg)
1726 : : ! Identify the first symmetry operation that is order 3, 4 or 6
1727 : 3068 : found=0
1728 [ + + ]: 16842 : do isym=1,nsym
1729 [ + + + + ]: 235463 : if( minval( abs( t_axes(isym)-(/10,12,14,22,23,24,25,26,27,28,29,30,31/) ))==0) then
1730 : : ! Proper rotations
1731 : : found=1 ; exit
1732 : :
1733 [ + + + + ]: 59188 : else if( minval( abs( t_axes(isym)-(/1,2,3/) ))==0) then
1734 : : ! Improper symmetry operations
1735 : 1023 : found=-1 ; exit
1736 : : end if
1737 : : end do
1738 [ + + ]: 3068 : if(found==-1 .or. found==1)then
1739 [ + + + + ]: 39559 : symrel_uni=symrel(:,:,isym)
1740 [ + + + + : 15319 : if(found==-1)symrel_uni=-symrel_uni
+ + ]
1741 : : ! Now, symrel_uni is a rotation of order 3, 4, 6, for which the axis must be identified
1742 : : ! It is actually the only eigenvector with eigenvalue 1. It can be found by cross products
1743 : : ! Subtract the unit matrix.
1744 [ + + ]: 12172 : do ii=1,3
1745 : 12172 : symrel_uni(ii,ii)=symrel_uni(ii,ii)-1
1746 : : end do
1747 : 3658 : found=0
1748 [ - + ]: 3658 : do ii=1,3
1749 [ + + ]: 3658 : jj=ii+1 ; if(jj==4)jj=1
1750 : : ! Cross product
1751 [ + + + + ]: 25606 : ivec1 = symrel_uni(ii,:); ivec2 = symrel_uni(jj,:)
1752 : 3658 : uniaxis = ivec1 .x. ivec2
1753 [ + + + + ]: 14632 : if(sum(uniaxis**2)/=0)then
1754 : : found=1 ; exit
1755 : : end if
1756 : : end do
1757 [ + - ]: 3043 : if(found==1)then
1758 : : ! Try to reduce the length, by an integer factor (try only primes 2, 3, 5, 7, 11)
1759 : 3043 : prime=(/2,3,5,7,11/)
1760 : 3043 : ii=1
1761 [ + + ]: 19828 : do while (ii<6)
1762 [ + + ]: 67140 : uniaxis_try=uniaxis/prime(ii)
1763 [ + + + + ]: 70183 : if(sum(abs(uniaxis_try*prime(ii)-uniaxis))==0)then
1764 : 1570 : uniaxis=uniaxis_try
1765 : : else
1766 : 15215 : ii=ii+1
1767 : : end if
1768 : : end do
1769 : 3043 : write(msg, '(a,3i4)' ) ' Optical axis (in reduced coordinates, real space ) :',uniaxis
1770 : : end if
1771 : : end if
1772 : : if(found==0)then
1773 : 25 : write(msg, '(a)' ) ' However, the axis has not been found. Sorry for this.'
1774 : : end if
1775 : 3068 : call wrtout(std_out,msg)
1776 : : end if
1777 : :
1778 [ - + ]: 11744 : ABI_FREE(determinant)
1779 [ - + ]: 11744 : ABI_FREE(shift)
1780 [ - + ]: 11744 : ABI_FREE(symrelconv)
1781 [ - + ]: 11744 : ABI_FREE(tnonsconv)
1782 [ - + ]: 11744 : ABI_FREE(t_axes)
1783 : :
1784 : : DBG_EXIT("COLL")
1785 : :
1786 : 11744 : end subroutine symspgr
1787 : : !!***
1788 : :
1789 : : !!****f* m_symfind/symlatt
1790 : : !! NAME
1791 : : !! symlatt
1792 : : !!
1793 : : !! FUNCTION
1794 : : !! From the unit cell vectors (rprimd) and the corresponding metric tensor,
1795 : : !! find the Bravais lattice and its symmetry operations (ptsymrel).
1796 : : !! 1) Find the shortest possible primitive vectors for the lattice
1797 : : !! 2) Determines the holohedral group of the lattice, and the
1798 : : !! axes to be used for the conventional cell
1799 : : !! (this is a delicate part, in which the centering of the
1800 : : !! reduced cell must be taken into account)
1801 : : !! The idea is to determine the basis vectors of the conventional
1802 : : !! cell from the reduced cell basis vectors.
1803 : : !! 3) Generate the symmetry operations of the holohedral group
1804 : : !!
1805 : : !! INPUTS
1806 : : !! iout=unit number of output file
1807 : : !! msym=default maximal number of symmetries. WARNING : cannot be simply set to nsym, because
1808 : : !! the number of symmetries found here will likely be bigger than sym !
1809 : : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
1810 : : !! tolsym=tolerance for the symmetries
1811 : : !!
1812 : : !! OUTPUT
1813 : : !! bravais(11): bravais(1)=iholohedry
1814 : : !! bravais(2)=center
1815 : : !! bravais(3:11)=coordinates of rprim in the axes
1816 : : !! of the conventional bravais lattice (*2 if center/=0)
1817 : : !! nptsym=number of point symmetries of the Bravais lattice
1818 : : !! ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
1819 : : !! of the Bravais lattice in real space in terms of primitive translations.
1820 : : !!
1821 : : !! NOTES
1822 : : !! WARNING: bravais(1) might be given a negative value in another
1823 : : !! routine, if the cell is non-primitive.
1824 : : !! The holohedral groups are numbered as follows
1825 : : !! (see international tables for crystallography (1983), p. 13)
1826 : : !! iholohedry=1 triclinic 1bar
1827 : : !! iholohedry=2 monoclinic 2/m
1828 : : !! iholohedry=3 orthorhombic mmm
1829 : : !! iholohedry=4 tetragonal 4/mmm
1830 : : !! iholohedry=5 trigonal 3bar m
1831 : : !! iholohedry=6 hexagonal 6/mmm
1832 : : !! iholohedry=7 cubic m3bar m
1833 : : !! Centering
1834 : : !! center=0 no centering
1835 : : !! center=-1 body-centered
1836 : : !! center=-3 face-centered
1837 : : !! center=1 A-face centered
1838 : : !! center=2 B-face centered
1839 : : !! center=3 C-face centered
1840 : : !!
1841 : : !! SOURCE
1842 : :
1843 : 40291 : subroutine symlatt(bravais,iout,msym,nptsym,ptsymrel,rprimd,tolsym)
1844 : :
1845 : : !Arguments ------------------------------------
1846 : : !scalars
1847 : : integer,intent(in) :: iout,msym
1848 : : integer,intent(out) :: nptsym
1849 : : real(dp),intent(in) :: tolsym
1850 : : !arrays
1851 : : integer,intent(out) :: bravais(11),ptsymrel(3,3,msym)
1852 : : real(dp),intent(in) :: rprimd(3,3)
1853 : :
1854 : : !Local variables-------------------------------
1855 : : !scalars
1856 : : integer,parameter :: mgen=4
1857 : : integer :: center,fact,found,foundc,ia,iaxis1,iaxis2
1858 : : integer :: isign1,isign2,ib,icase,igen,iholohedry,ii,index,isym
1859 : : integer :: itrial,jj,jsym,ngen=0,orthogonal,sign12,sign13,sign23,sumsign
1860 : : real(dp) :: determinant,norm2a,norm2b,norm2c,norm2trial,reduceda,reducedb,sca
1861 : : real(dp) :: scalarprod,scb,trace,trace_best,val
1862 : : character(len=500) :: msg
1863 : : !arrays
1864 : : integer,parameter :: list_holo(7)=(/7,6,4,3,5,2,1/)
1865 : : integer :: ang90(3),equal(3),gen(3,3,mgen),gen2xy(3,3),gen2y(3,3),gen2z(3,3)
1866 : : integer :: gen3(3,3),gen6(3,3),icoord(3,3),identity(3,3),nvecta(3),nvectb(3)
1867 : : integer :: order(mgen)
1868 : : real(dp) :: axes(3,3),axesinvt(3,3),axes_best(3,3),axes_try(3,3)
1869 : : real(dp) :: cell_base(3,3),coord(3,3),metmin(3,3)
1870 : : real(dp) :: minim(3,3),scprods(3,3),vecta(3),vectb(3),vectc(3),vin1(3),vin2(3),vext(3)
1871 : : !**************************************************************************
1872 : :
1873 : : !DEBUG
1874 : : !write(std_out,'(a,es14.6)') ' m_symfind%symlatt : enter, tolsym= ',tolsym
1875 : : !call flush(std_out)
1876 : : !ENDDEBUG
1877 : :
1878 : 40291 : identity(:,:)=0 ; identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1
1879 : 40291 : nvecta(1)=2 ; nvectb(1)=3
1880 : 40291 : nvecta(2)=1 ; nvectb(2)=3
1881 : 40291 : nvecta(3)=1 ; nvectb(3)=2
1882 : :
1883 : : !--------------------------------------------------------------------------
1884 : : !Reduce the input vectors to a set of minimal vectors
1885 : 40291 : call smallprim(metmin,minim,rprimd)
1886 : :
1887 : : !DEBUG
1888 : : !write(std_out,*)' symlatt : minim(:,1)=',minim(:,1)
1889 : : !write(std_out,*)' symlatt : minim(:,2)=',minim(:,2)
1890 : : !write(std_out,*)' symlatt : minim(:,3)=',minim(:,3)
1891 : : !call flush(std_out)
1892 : : !ENDDEBUG
1893 : :
1894 : : !--------------------------------------------------------------------------
1895 : : !Examine the angles and vector lengths
1896 : 40291 : ang90(:)=0
1897 [ + + ]: 40291 : if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1
1898 [ + + ]: 40291 : if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1
1899 [ + + ]: 40291 : if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1
1900 : 40291 : equal(:)=0
1901 [ + + ]: 40291 : if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1
1902 [ + + ]: 40291 : if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1
1903 [ + + ]: 40291 : if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1
1904 : :
1905 : : !DEBUG
1906 : : !write(std_out,*)' ang90=',ang90(:)
1907 : : !write(std_out,*)' equal=',equal(:)
1908 : : !call flush(std_out)
1909 : : !ENDDEBUG
1910 : :
1911 : : !-----------------------------------------------------------------------
1912 : : !Identification of the centering
1913 : :
1914 : 40291 : foundc=0
1915 : : !Default values
1916 : 40291 : fact=1 ; center=0
1917 : 40291 : cell_base(:,:)=minim(:,:)
1918 : :
1919 : : !Examine each holohedral group
1920 : : !This search is ordered : should not be happy with tetragonal,
1921 : : !while there is FCC ...
1922 [ + + ]: 126805 : do index=1,6
1923 : :
1924 : : ! If the holohedry is already found, exit
1925 [ + + ]: 123015 : if(foundc==1)exit
1926 : :
1927 : : ! Initialize the target holohedry
1928 : 86514 : iholohedry=list_holo(index)
1929 : :
1930 : : ! DEBUG
1931 : : ! write(std_out,*)
1932 : : ! write(std_out,*)' symlatt : trial holohedry',iholohedry
1933 : : ! ENDDEBUG
1934 : :
1935 : 86514 : orthogonal=0
1936 [ + + + + ]: 86514 : if(iholohedry==7 .or. iholohedry==4 .or. iholohedry==3)orthogonal=1
1937 : :
1938 : : ! Now, will examine different working hypothesis.
1939 : : ! The set of these hypothesis is thought to cover all possible cases ...
1940 : :
1941 : : ! Working hypothesis : the basis is orthogonal
1942 [ + + + + ]: 86514 : if(ang90(1)+ang90(2)+ang90(3)==3 .and. orthogonal==1)then
1943 : 21501 : fact=1 ; center=0
1944 : 21501 : cell_base(:,:)=minim(:,:)
1945 : : ! Checks that the basis vectors are OK for the target holohedry
1946 : 21501 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1947 : : end if
1948 : :
1949 : : ! Select one trial direction
1950 [ + + ]: 303629 : do itrial=1,3
1951 : :
1952 : : ! If the holohedry is already found, exit
1953 [ + + ]: 211786 : if(foundc==1)exit
1954 : :
1955 : 176824 : ia=nvecta(itrial) ; ib=nvectb(itrial)
1956 : :
1957 : : ! This is in case of hexagonal holohedry
1958 [ + - + + : 176824 : if(foundc==0 .and. iholohedry==6 .and. ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then
+ + + + +
+ ]
1959 : 7406 : reduceda=metmin(ib,ia)/metmin(ia,ia)
1960 : 7406 : fact=1 ; center=0
1961 [ + + ]: 7406 : if(abs(reduceda+0.5d0)<tolsym)then
1962 [ + + ]: 8520 : cell_base(:,1)=minim(:,ia)
1963 [ + + ]: 8520 : cell_base(:,2)=minim(:,ib)
1964 [ + + ]: 8520 : cell_base(:,3)=minim(:,itrial)
1965 : : !DEBUG
1966 : : ! write(std_out,*)' cell_base(:,1)=',cell_base(:,1)
1967 : : ! write(std_out,*)' cell_base(:,2)=',cell_base(:,2)
1968 : : ! write(std_out,*)' cell_base(:,3)=',cell_base(:,3)
1969 : : !ENDDEBUG
1970 : : ! Checks that the basis vectors are OK for the target holohedry
1971 : 2130 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1972 [ + + ]: 5276 : else if(abs(reduceda-0.5d0)<tolsym)then
1973 [ + + ]: 276 : cell_base(:,1)=minim(:,ia)
1974 [ + + ]: 276 : cell_base(:,2)=minim(:,ib)-minim(:,ia)
1975 [ + + ]: 276 : cell_base(:,3)=minim(:,itrial)
1976 : : ! Checks that the basis vectors are OK for the target holohedry
1977 : 69 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1978 : : end if
1979 : : end if
1980 : :
1981 : : ! Working hypothesis : the conventional cell is orthogonal,
1982 : : ! and the two other vectors are axes of the conventional cell
1983 [ + + + + : 176824 : if(foundc==0 .and. orthogonal==1 .and. ang90(itrial)==1)then
+ + ]
1984 : :
1985 : : ! Compute the reduced coordinate of trial vector in the basis
1986 : : ! of the two other vectors
1987 : 37640 : reduceda=metmin(itrial,ia)/metmin(ia,ia)
1988 : 37640 : reducedb=metmin(itrial,ib)/metmin(ib,ib)
1989 [ + + ]: 150560 : cell_base(:,ia)=minim(:,ia)
1990 [ + + ]: 150560 : cell_base(:,ib)=minim(:,ib)
1991 [ + + + + : 37640 : if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. &
+ + ]
1992 [ + + ]: 37640 : & ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym) )then
1993 [ + + ]: 4602 : if(abs(abs(reduceda)-0.5d0)<tolsym)center=ib
1994 [ + + ]: 4602 : if(abs(abs(reducedb)-0.5d0)<tolsym)center=ia
1995 : 4602 : fact=2
1996 : : cell_base(:,itrial)= &
1997 [ + + ]: 18408 : & (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0
1998 : 4602 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
1999 [ + + + + ]: 33038 : else if( abs(abs(reduceda)-0.5d0)<tolsym .and.&
2000 : : & abs(abs(reducedb)-0.5d0)<tolsym ) then
2001 : 569 : fact=2 ; center=-1
2002 : : cell_base(:,itrial)= &
2003 [ + + ]: 2276 : & (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0
2004 : 569 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2005 : : end if
2006 : : end if
2007 : :
2008 : : ! Working hypothesis : the conventional cell is orthogonal, and
2009 : : ! the trial vector is one of the future axes, and the face perpendicular to it is centered
2010 : : if(foundc==0 .and. iholohedry==3 .and. &
2011 [ + + + + : 176824 : & ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then
+ + + + +
+ ]
2012 : 460 : fact=2 ; center=itrial
2013 [ + + ]: 1840 : cell_base(:,ia)=minim(:,ia)+minim(:,ib)
2014 [ + + ]: 1840 : cell_base(:,ib)=minim(:,ia)-minim(:,ib)
2015 [ + + ]: 1840 : cell_base(:,itrial)=minim(:,itrial)
2016 : : ! Checks that the basis vectors are OK for the target holohedry
2017 : 460 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2018 : : end if
2019 : :
2020 : : ! DEBUG
2021 : : ! write(std_out,*)' after test_b, foundc=',foundc
2022 : : ! ENDDEBUG
2023 : :
2024 : : ! Working hypothesis : the conventional cell is orthogonal, and
2025 : : ! the trial vector is one of the future axes
2026 [ + + + + ]: 176824 : if(foundc==0 .and. orthogonal==1)then
2027 : : ! Compute the projection of the two other vectors on the trial vector
2028 : 106977 : reduceda=metmin(itrial,ia)/metmin(itrial,itrial)
2029 : 106977 : reducedb=metmin(itrial,ib)/metmin(itrial,itrial)
2030 : : ! If both projections are half-integer, one might have found an axis
2031 [ + + + + ]: 106977 : if( abs(abs(reduceda)-0.5d0)<tolsym .and.&
2032 : : & abs(abs(reducedb)-0.5d0)<tolsym ) then
2033 [ + + ]: 96656 : vecta(:)=minim(:,ia)-reduceda*minim(:,itrial)
2034 [ + + ]: 96656 : vectb(:)=minim(:,ib)-reducedb*minim(:,itrial)
2035 : 24164 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2036 : 24164 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2037 : 24164 : scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
2038 : : ! Note the order of selection : body-centered is prefered
2039 : : ! over face centered, which is correct for the tetragonal case
2040 [ + + ]: 24164 : if(abs(norm2a-norm2b)<tolsym*half*(norm2a+norm2b))then
2041 : : ! The lattice is body centered
2042 : 17279 : fact=2 ; center=-1
2043 [ + + ]: 69116 : cell_base(:,ia)=vecta(:)+vectb(:)
2044 [ + + ]: 69116 : cell_base(:,ib)=vecta(:)-vectb(:)
2045 [ + + ]: 69116 : cell_base(:,itrial)=minim(:,itrial)
2046 : 17279 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2047 [ + - ]: 6885 : else if(abs(scalarprod)<tolsym*half*(norm2a+norm2b))then
2048 : : ! The lattice is face centered
2049 : 0 : fact=2 ; center=-3
2050 [ # # ]: 0 : cell_base(:,ia)=2.0d0*vecta(:)
2051 [ # # ]: 0 : cell_base(:,ib)=2.0d0*vectb(:)
2052 [ # # ]: 0 : cell_base(:,itrial)=minim(:,itrial)
2053 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2054 : : end if
2055 : : end if
2056 : : end if
2057 : :
2058 : : ! DEBUG
2059 : : ! write(std_out,*)' after test_c, foundc=',foundc
2060 : : ! ENDDEBUG
2061 : :
2062 : : ! Working hypothesis : the conventional cell is orthogonal,
2063 : : ! and body centered with no basis vector being an axis,
2064 : : ! in which case the basis vectors must be equal (even for orthorhombic)
2065 : : if(foundc==0 .and. orthogonal==1 .and. &
2066 [ + + + + : 176824 : & equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then
+ + + + +
- ]
2067 : : ! Compute the combination of the two other vectors
2068 [ + + ]: 93148 : vecta(:)=minim(:,ia)+minim(:,ib)
2069 [ + + ]: 93148 : vectb(:)=minim(:,ia)-minim(:,ib)
2070 : 23287 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2071 : 23287 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2072 : : ! Project the trial vector on the first of the two vectors
2073 : : reduceda=( minim(1,itrial)*vecta(1)+ &
2074 : : & minim(2,itrial)*vecta(2)+ &
2075 : 23287 : & minim(3,itrial)*vecta(3) )/norm2a
2076 : : reducedb=( minim(1,itrial)*vectb(1)+ &
2077 : : & minim(2,itrial)*vectb(2)+ &
2078 : 23287 : & minim(3,itrial)*vectb(3) )/norm2b
2079 [ + + ]: 23287 : if( abs(abs(reduceda)-0.5d0)<tolsym )then
2080 : : ! The first vector is an axis
2081 : 3653 : fact=2 ; center=-1
2082 [ + + ]: 14612 : cell_base(:,ia)=vecta(:)
2083 [ + + ]: 14612 : vecta(:)=minim(:,itrial)-reduceda*vecta(:)
2084 [ + + ]: 14612 : vectb(:)=0.5d0*vectb(:)
2085 [ + + ]: 14612 : cell_base(:,ib)=vecta(:)+vectb(:)
2086 [ + + ]: 14612 : cell_base(:,itrial)=vecta(:)-vectb(:)
2087 : 3653 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2088 [ + + ]: 19634 : else if( abs(abs(reducedb)-0.5d0)<tolsym )then
2089 : : ! The second vector is an axis
2090 : 84 : fact=2 ; center=-1
2091 [ + + ]: 336 : cell_base(:,ib)=vectb(:)
2092 [ + + ]: 336 : vectb(:)=minim(:,itrial)-reducedb*vectb(:)
2093 [ + + ]: 336 : vecta(:)=0.5d0*vecta(:)
2094 [ + + ]: 336 : cell_base(:,ia)=vectb(:)+vecta(:)
2095 [ + + ]: 336 : cell_base(:,itrial)=vectb(:)-vecta(:)
2096 : 84 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2097 : : end if
2098 : : end if
2099 : :
2100 : : ! Working hypothesis : the conventional cell is orthogonal,
2101 : : ! and face centered, in the case where two minimal vectors are equal
2102 [ + + + + : 176824 : if(foundc==0 .and. orthogonal==1 .and. equal(itrial)==1 ) then
+ + ]
2103 : : ! Compute the combination of these two vectors
2104 [ + + ]: 156768 : vecta(:)=minim(:,ia)+minim(:,ib)
2105 [ + + ]: 156768 : vectb(:)=minim(:,ia)-minim(:,ib)
2106 : 39192 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2107 : 39192 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2108 : : ! Project the trial vector on the two vectors
2109 : : reduceda=( minim(1,itrial)*vecta(1)+ &
2110 : : & minim(2,itrial)*vecta(2)+ &
2111 : 39192 : & minim(3,itrial)*vecta(3) )/norm2a
2112 : : reducedb=( minim(1,itrial)*vectb(1)+ &
2113 : : & minim(2,itrial)*vectb(2)+ &
2114 : 39192 : & minim(3,itrial)*vectb(3) )/norm2b
2115 [ + + + + : 39192 : if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. &
+ + ]
2116 [ + + ]: 39192 : & ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym) )then
2117 : 986 : fact=2 ; center=-3
2118 : : cell_base(:,itrial)= &
2119 [ + + ]: 3944 : & (minim(:,itrial)-reduceda*vecta(:)-reducedb*vectb(:) )*2.0d0
2120 [ + + ]: 3944 : cell_base(:,ia)=vecta(:)
2121 [ + + ]: 3944 : cell_base(:,ib)=vectb(:)
2122 : 986 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2123 : : end if
2124 : : end if
2125 : :
2126 : : ! Working hypothesis : the conventional cell is orthogonal,
2127 : : ! face centered, but no two vectors are on the same "square"
2128 [ + + + + ]: 176824 : if(foundc==0 .and. orthogonal==1)then
2129 : : ! Compute the combination of these two vectors
2130 [ + + ]: 419628 : vecta(:)=minim(:,ia)+minim(:,ib)
2131 [ + + ]: 419628 : vectb(:)=minim(:,ia)-minim(:,ib)
2132 : 104907 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2133 : 104907 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2134 : : ! The trial vector length must be equal to one of these lengths
2135 [ + + ]: 104907 : if(abs(metmin(itrial,itrial)-norm2a)<tolsym*norm2a)then
2136 : 1211 : fact=2 ; center=-3
2137 [ + + ]: 4844 : cell_base(:,ia)=vecta(:)+minim(:,itrial)
2138 [ + + ]: 4844 : cell_base(:,ib)=vecta(:)-minim(:,itrial)
2139 : : ! Project vectb perpendicular to cell_base(:,ia) and cell_base(:,ib)
2140 : 1211 : norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2
2141 : 1211 : norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2
2142 : : reduceda=( cell_base(1,ia)*vectb(1)+ &
2143 : : & cell_base(2,ia)*vectb(2)+ &
2144 : 1211 : & cell_base(3,ia)*vectb(3) )/norm2a
2145 : : reducedb=( cell_base(1,ib)*vectb(1)+ &
2146 : : & cell_base(2,ib)*vectb(2)+ &
2147 : 1211 : & cell_base(3,ib)*vectb(3) )/norm2b
2148 [ + + + - ]: 1211 : if( abs(abs(reduceda)-0.5d0)<tolsym .and. &
2149 : : & abs(abs(reducedb)-0.5d0)<tolsym )then
2150 [ # # ]: 0 : cell_base(:,itrial)=vectb(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib)
2151 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2152 : : end if
2153 [ + + ]: 103696 : else if(abs(metmin(itrial,itrial)-norm2b)<tolsym*norm2b)then
2154 : 17405 : fact=2 ; center=-3
2155 [ + + ]: 69620 : cell_base(:,ia)=vectb(:)+minim(:,itrial)
2156 [ + + ]: 69620 : cell_base(:,ib)=vectb(:)-minim(:,itrial)
2157 : : ! Project vecta perpendicular to cell_base(:,ia) and cell_base(:,ib)
2158 : 17405 : norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2
2159 : 17405 : norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2
2160 : : reduceda=( cell_base(1,ia)*vecta(1)+ &
2161 : : & cell_base(2,ia)*vecta(2)+ &
2162 : 17405 : & cell_base(3,ia)*vecta(3) )/norm2a
2163 : : reducedb=( cell_base(1,ib)*vecta(1)+ &
2164 : : & cell_base(2,ib)*vecta(2)+ &
2165 : 17405 : & cell_base(3,ib)*vecta(3) )/norm2b
2166 [ + + + + ]: 17405 : if( abs(abs(reduceda)-0.5d0)<tolsym .and. &
2167 : : & abs(abs(reducedb)-0.5d0)<tolsym )then
2168 [ + + ]: 66572 : cell_base(:,itrial)=vecta(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib)
2169 : 16643 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2170 : : end if
2171 : : end if
2172 : : end if
2173 : :
2174 : : ! Working hypothesis : the cell is rhombohedral, and
2175 : : ! the three minimal vectors have same length and same absolute scalar product
2176 : : if(foundc==0 .and. iholohedry==5 .and. &
2177 [ + + + + : 176824 : & equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then
+ + + + +
- ]
2178 [ + - ]: 323 : if( abs(abs(metmin(1,2))-abs(metmin(1,3)))<tolsym*metmin(1,1) .and. &
2179 [ + - ]: 323 : & abs(abs(metmin(1,2))-abs(metmin(2,3)))<tolsym*metmin(2,2) )then
2180 : 323 : fact=1 ; center=0
2181 : 323 : cell_base(:,:)=minim(:,:)
2182 : : ! One might have to change the sign of one of the vectors
2183 : 323 : sign12=1 ; sign13=1 ; sign23=1
2184 [ + + ]: 323 : if(metmin(1,2)<0.0d0)sign12=-1
2185 [ + + ]: 323 : if(metmin(1,3)<0.0d0)sign13=-1
2186 [ + + ]: 323 : if(metmin(2,3)<0.0d0)sign23=-1
2187 : 323 : sumsign=sign12+sign13+sign23
2188 [ + + ]: 323 : if(sumsign==-1)then
2189 [ - + - - ]: 50 : if(sign12==1)cell_base(:,3)=-cell_base(:,3)
2190 [ - + - - ]: 50 : if(sign13==1)cell_base(:,2)=-cell_base(:,2)
2191 [ + - + + ]: 200 : if(sign23==1)cell_base(:,1)=-cell_base(:,1)
2192 [ + + ]: 273 : else if(sumsign==1)then
2193 [ - + - - ]: 4 : if(sign12==-1)cell_base(:,3)=-cell_base(:,3)
2194 [ - + - - ]: 4 : if(sign13==-1)cell_base(:,2)=-cell_base(:,2)
2195 [ + - + + ]: 16 : if(sign23==-1)cell_base(:,1)=-cell_base(:,1)
2196 : : end if
2197 : 323 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2198 : : end if
2199 : : end if
2200 : :
2201 : : ! DEBUG
2202 : : ! write(std_out,*)' after test_3a, foundc=',foundc
2203 : : ! write(std_out,*)' after test_3a, itrial=',itrial
2204 : : ! write(std_out,*)' after test_3a, equal(:)=',equal(:)
2205 : : ! ENDDEBUG
2206 : :
2207 : : ! Working hypothesis : the cell is rhombohedral, one vector
2208 : : ! is parallel to the trigonal axis
2209 [ + + + + : 176824 : if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 )then
+ + ]
2210 [ + + + + ]: 20706 : vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib)
2211 : 2958 : norm2trial=minim(1,itrial)**2+minim(2,itrial)**2+minim(3,itrial)**2
2212 : : reduceda=( minim(1,itrial)*vecta(1)+ &
2213 : : & minim(2,itrial)*vecta(2)+ &
2214 : 2958 : & minim(3,itrial)*vecta(3) )/norm2trial
2215 : : reducedb=( minim(1,itrial)*vectb(1)+ &
2216 : : & minim(2,itrial)*vectb(2)+ &
2217 : 2958 : & minim(3,itrial)*vectb(3) )/norm2trial
2218 : : ! DEBUG
2219 : : ! write(std_out,*)' reduceda,reducedb=',reduceda,reducedb
2220 : : ! ENDDEBUG
2221 [ + + - + ]: 2958 : if(abs(abs(reduceda)-1.0d0/3.0d0)<tolsym .and. &
2222 : : & abs(abs(reducedb)-1.0d0/3.0d0)<tolsym ) then
2223 : : ! Possibly change of sign to make positive the scalar product with
2224 : : ! the vector parallel to the trigonal axis
2225 [ + + + + ]: 82 : if(reduceda<zero)vecta(:)=-vecta(:)
2226 [ + + + + ]: 130 : if(reducedb<zero)vectb(:)=-vectb(:)
2227 : : ! Projection on the orthogonal plane
2228 [ + + ]: 148 : vecta(:)=vecta(:)-abs(reduceda)*cell_base(:,itrial)
2229 [ + + ]: 148 : vectb(:)=vectb(:)-abs(reducedb)*cell_base(:,itrial)
2230 : : ! These two vectors should have an angle of 120 degrees
2231 : 37 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2232 : 37 : scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
2233 : : ! DEBUG
2234 : : ! write(std_out,*)' norm2a,scalarprod=',norm2a,scalarprod
2235 : : ! ENDDEBUG
2236 [ + + ]: 37 : if(abs(two*scalarprod+norm2a)<tolsym*norm2a)then
2237 : 12 : fact=1 ; center=0
2238 [ - + - - ]: 12 : if(scalarprod>0.0d0)vectb(:)=-vectb(:)
2239 : : ! Now vecta and vectb have an angle of 120 degrees
2240 [ + + ]: 48 : cell_base(:,1)=cell_base(:,itrial)/3.0d0+vecta(:)
2241 [ + + ]: 48 : cell_base(:,2)=cell_base(:,itrial)/3.0d0+vectb(:)
2242 [ + + ]: 48 : cell_base(:,3)=cell_base(:,itrial)/3.0d0-vecta(:)-vectb(:)
2243 : : ! DEBUG
2244 : : ! write(std_out,*)' cell_base(:,1)=',cell_base(:,1)
2245 : : ! write(std_out,*)' cell_base(:,2)=',cell_base(:,2)
2246 : : ! write(std_out,*)' cell_base(:,3)=',cell_base(:,3)
2247 : : ! ENDDEBUG
2248 : 12 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2249 : : end if
2250 : : end if
2251 : : end if
2252 : :
2253 : : ! Working hypothesis : the cell is rhombohedral, one vector
2254 : : ! is in the plane perpendicular to the trigonal axis
2255 [ + + + + : 176824 : if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then
+ + ]
2256 [ + + ]: 11784 : vecta(:)=minim(:,ia)+minim(:,ib)
2257 [ + + ]: 11784 : vectb(:)=minim(:,ia)-minim(:,ib)
2258 : 2946 : norm2trial=cell_base(1,itrial)**2+cell_base(2,itrial)**2+cell_base(3,itrial)**2
2259 : 2946 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2260 : 2946 : norm2b=vecta(1)**2+vecta(2)**2+vecta(3)**2
2261 : : reduceda=( cell_base(1,itrial)*vecta(1)+ &
2262 : : & cell_base(2,itrial)*vecta(2)+ &
2263 : 2946 : & cell_base(3,itrial)*vecta(3) )/norm2trial
2264 : : reducedb=( cell_base(1,itrial)*vectb(1)+ &
2265 : : & cell_base(2,itrial)*vectb(2)+ &
2266 : 2946 : & cell_base(3,itrial)*vectb(3) )/norm2trial
2267 [ + + - + ]: 2946 : if(abs(norm2trial-norm2a)<tolsym*norm2a .and. &
2268 : : & abs(abs(2*reduceda)-norm2trial)<tolsym*norm2trial )then
2269 : 0 : fact=1 ; center=0
2270 [ # # ]: 0 : cell_base(:,1)=minim(:,ia)
2271 [ # # ]: 0 : cell_base(:,2)=-minim(:,ib)
2272 [ # # ]: 0 : cell_base(:,3)=-minim(:,ib)+2*reduceda*minim(:,itrial)
2273 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2274 [ + + - + ]: 2946 : else if (abs(norm2trial-norm2b)<tolsym*norm2b .and. &
2275 : : & abs(abs(2*reducedb)-norm2trial)<tolsym*norm2trial )then
2276 : 0 : fact=1 ; center=0
2277 [ # # ]: 0 : cell_base(:,1)=minim(:,ia)
2278 [ # # ]: 0 : cell_base(:,2)=minim(:,ib)
2279 [ # # ]: 0 : cell_base(:,3)=minim(:,ib)+2*reducedb*minim(:,itrial)
2280 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2281 : : end if
2282 : : end if
2283 : :
2284 : : ! Working hypothesis : the cell is rhombohedral, two vectors
2285 : : ! are in the plane perpendicular to the trigonal axis
2286 [ + + + + : 176824 : if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then
+ + ]
2287 [ + + + + ]: 20622 : vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib)
2288 : 2946 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2289 : 2946 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2290 : 2946 : scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3)
2291 [ + + ]: 2946 : if(abs(abs(2*scalarprod)-norm2a)<tolsym*norm2a)then
2292 : : ! This is in order to have 120 angle between vecta and vectb
2293 [ + + + + ]: 1957 : if(scalarprod>0.0d0)vectb(:)=-vectb(:)
2294 : : reduceda=( cell_base(1,itrial)*vecta(1)+ &
2295 : : & cell_base(2,itrial)*vecta(2)+ &
2296 : 1246 : & cell_base(3,itrial)*vecta(3) )/norm2a
2297 : : reducedb=( cell_base(1,itrial)*vectb(1)+ &
2298 : : & cell_base(2,itrial)*vectb(2)+ &
2299 : 1246 : & cell_base(3,itrial)*vectb(3) )/norm2b
2300 : 1246 : fact=1 ; center=0
2301 [ + + ]: 4984 : cell_base(:,1)=minim(:,itrial)
2302 [ + + + + ]: 1246 : if(abs(reduceda-0.5d0)<tolsym .and. abs(reducedb)<tolsym )then
2303 [ + + ]: 900 : cell_base(:,2)=minim(:,itrial)-vecta(:)
2304 [ + + ]: 900 : cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:)
2305 : 225 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2306 [ + + + - ]: 1021 : else if(abs(reduceda-0.5d0)<tolsym.and. abs(reducedb+0.5d0)<tolsym )then
2307 [ + + ]: 24 : cell_base(:,2)=minim(:,itrial)-vecta(:)
2308 [ + + ]: 24 : cell_base(:,3)=minim(:,itrial)+vectb(:)
2309 : 6 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2310 [ - + - - ]: 1015 : else if(abs(reduceda)<tolsym .and. abs(reducedb+0.5d0)<tolsym )then
2311 [ # # ]: 0 : cell_base(:,2)=minim(:,itrial)+vectb(:)
2312 [ # # ]: 0 : cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:)
2313 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2314 [ + - + + ]: 1015 : else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb)<tolsym)then
2315 [ + + ]: 24 : cell_base(:,2)=minim(:,itrial)+vecta(:)
2316 [ + + ]: 24 : cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:)
2317 : 6 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2318 [ + - + - ]: 1009 : else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb-0.5d0)<tolsym)then
2319 [ + + ]: 4036 : cell_base(:,2)=minim(:,itrial)+vecta(:)
2320 [ + + ]: 4036 : cell_base(:,3)=minim(:,itrial)-vectb(:)
2321 : 1009 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2322 [ # # # # ]: 0 : else if(abs(reduceda)<tolsym .and. abs(reducedb-0.5d0)<tolsym )then
2323 [ # # ]: 0 : cell_base(:,2)=minim(:,itrial)-vectb(:)
2324 [ # # ]: 0 : cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:)
2325 : 0 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2326 : : end if
2327 : : end if
2328 : : end if
2329 : :
2330 : : ! Working hypothesis : monoclinic holohedry, primitive. Then, two angles are 90 degrees
2331 : : if(foundc==0 .and. iholohedry==2 .and. &
2332 [ + + + + : 176824 : & ang90(ia)==1 .and. ang90(ib)==1 ) then
+ + + + ]
2333 : 299 : fact=1 ; center=0
2334 [ + + ]: 1196 : cell_base(:,1)=minim(:,ia)
2335 [ + + ]: 1196 : cell_base(:,2)=minim(:,itrial)
2336 [ + + ]: 1196 : cell_base(:,3)=minim(:,ib)
2337 : : ! Checks that the basis vectors are OK for the target holohedry
2338 : 299 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2339 : : end if
2340 : :
2341 : : ! Monoclinic holohedry, one-face-centered cell
2342 : : ! Working hypothesis, two vectors have equal length.
2343 [ + + ]: 1060944 : do icase=1,5
2344 [ + + + + : 1060944 : if(foundc==0 .and. iholohedry==2 .and. equal(itrial)==1 ) then
+ + ]
2345 [ + + ]: 15460 : vecta(:)=cell_base(:,ia)+cell_base(:,ib)
2346 [ + + ]: 15460 : vectb(:)=cell_base(:,ia)-cell_base(:,ib)
2347 : 3865 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2348 : 3865 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2349 : : ! The minim(:,trial) vector belongs to the
2350 : : ! plane parallel to the cell_base(:,ia),cell_base(:,ib) plane
2351 : : ! In that plane, must try minim(:,itrial),
2352 : : ! as well as the 4 different combinations of
2353 : : ! minim(:,itrial) with the vectors in the plane
2354 [ + + + + ]: 7138 : if(icase==1)vectc(:)=minim(:,itrial)
2355 [ + + + + ]: 5956 : if(icase==2)vectc(:)=minim(:,itrial)+cell_base(:,ia)
2356 [ + + + + ]: 5950 : if(icase==3)vectc(:)=minim(:,itrial)+cell_base(:,ib)
2357 [ + + + + ]: 5938 : if(icase==4)vectc(:)=minim(:,itrial)-cell_base(:,ia)
2358 [ + + + + ]: 5938 : if(icase==5)vectc(:)=minim(:,itrial)-cell_base(:,ib)
2359 : 3865 : norm2c=vectc(1)**2+vectc(2)**2+vectc(3)**2
2360 : : sca=vectc(1)*vecta(1)+&
2361 : : & vectc(2)*vecta(2)+&
2362 : 3865 : & vectc(3)*vecta(3)
2363 : : scb=vectc(1)*vectb(1)+&
2364 : : & vectc(2)*vectb(2)+&
2365 : 3865 : & vectc(3)*vectb(3)
2366 : : ! DEBUG
2367 : : ! write(std_out,*)' symlatt : test iholohedry=2, sca,scb=',sca,scb
2368 : : ! ENDDEBUG
2369 [ + + + + ]: 3865 : if(abs(sca)<tolsym*sqrt(norm2c*norm2a) .or. abs(scb)<tolsym*sqrt(norm2c*norm2b))then
2370 : 400 : fact=2 ; center=3
2371 : : ! The itrial direction is centered
2372 [ + + ]: 1600 : cell_base(:,3)=vectc(:)
2373 [ + + ]: 400 : if(abs(sca)<tolsym*sqrt(norm2c*norm2a))then
2374 [ + + ]: 528 : cell_base(:,2)=vecta(:)
2375 [ + + ]: 528 : cell_base(:,1)=vectb(:)
2376 : 132 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2377 [ - + ]: 268 : else if(abs(scb)<tolsym*sqrt(norm2c*norm2b))then
2378 [ + + ]: 1072 : cell_base(:,2)=vectb(:)
2379 [ + + ]: 1072 : cell_base(:,1)=vecta(:)
2380 : 268 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2381 : : end if
2382 : : end if
2383 : : end if
2384 : : end do ! icase=1,5
2385 : :
2386 : : ! Monoclinic holohedry, one-face-centered cell, but non equivalent.
2387 : : ! This case, one pair of vectors is orthogonal
2388 [ + + + + : 176824 : if(foundc==0 .and. iholohedry==2 .and. ang90(itrial)==1) then
+ + ]
2389 [ + + ]: 4712 : vecta(:)=minim(:,ia)
2390 [ + + ]: 4712 : vectb(:)=minim(:,ib)
2391 : 1178 : norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2
2392 : 1178 : norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2
2393 : : ! Project the trial vector on the two vectors
2394 : : reduceda=( minim(1,itrial)*vecta(1)+ &
2395 : : & minim(2,itrial)*vecta(2)+ &
2396 : 1178 : & minim(3,itrial)*vecta(3) )/norm2a
2397 : : reducedb=( minim(1,itrial)*vectb(1)+ &
2398 : : & minim(2,itrial)*vectb(2)+ &
2399 : 1178 : & minim(3,itrial)*vectb(3) )/norm2b
2400 [ + + + + ]: 1178 : if(abs(abs(reduceda)-0.5d0)<tolsym .or. abs(abs(reducedb)-0.5d0)<tolsym) then
2401 : 759 : fact=2 ; center=3
2402 [ + + ]: 759 : if(abs(abs(reduceda)-0.5d0)<tolsym)then
2403 [ + + ]: 20 : cell_base(:,2)=vecta(:)
2404 [ + + ]: 20 : cell_base(:,3)=vectb(:)
2405 [ + + ]: 20 : cell_base(:,1)=2*(minim(:,itrial)-reduceda*vecta(:))
2406 : 5 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2407 [ - + ]: 754 : else if(abs(abs(reducedb)-0.5d0)<tolsym)then
2408 [ + + ]: 3016 : cell_base(:,2)=vectb(:)
2409 [ + + ]: 3016 : cell_base(:,3)=vecta(:)
2410 [ + + ]: 3016 : cell_base(:,1)=2*(minim(:,itrial)-reducedb*vectb(:))
2411 : 754 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2412 : : end if
2413 : : end if
2414 : : end if
2415 : :
2416 : : ! Monoclinic holohedry, one-face-centered cell, but non equivalent.
2417 : : ! This case, no pair of vectors is orthogonal, no pair of vector of equal lentgh
2418 [ + + + + ]: 263338 : if(foundc==0 .and. iholohedry==2)then
2419 : : ! Try to find a vector that belongs to the mediator plane, or the binary vector.
2420 : : ! There must be one such vector, if centered monoclinic and no pair of vectors of equal length,
2421 : : ! either among the three vectors, or among one of their differences or sums.
2422 : : ! And there must be, among the two other vectors, one vector whose projection
2423 : : ! on this vector is half the length of this vector.
2424 [ + + ]: 29656 : vecta(:)=minim(:,ia)
2425 [ + + ]: 29656 : vectb(:)=minim(:,ib)
2426 : : ! Try the different possibilities for the vector on which the projection will be half ...
2427 [ + + ]: 32667 : do ii=1,5
2428 [ + + + + ]: 50282 : if(ii==1)vectc(:)=minim(:,itrial)
2429 [ + + + + ]: 47999 : if(ii==2)vectc(:)=minim(:,itrial)+vecta(:)
2430 [ + + + + ]: 42050 : if(ii==3)vectc(:)=minim(:,itrial)-vecta(:)
2431 [ + + + + ]: 42050 : if(ii==4)vectc(:)=minim(:,itrial)+vectb(:)
2432 [ + + + + ]: 41939 : if(ii==5)vectc(:)=minim(:,itrial)-vectb(:)
2433 : 28040 : norm2trial=vectc(1)**2+vectc(2)**2+vectc(3)**2
2434 : : ! Project the two vectors on the trial vector
2435 : 28040 : reduceda=( vectc(1)*vecta(1)+vectc(2)*vecta(2)+vectc(3)*vecta(3) )/norm2trial
2436 : 28040 : reducedb=( vectc(1)*vectb(1)+vectc(2)*vectb(2)+vectc(3)*vectb(3) )/norm2trial
2437 : 28040 : found=0
2438 [ + + ]: 28040 : if(abs(abs(reduceda)-0.5d0)<tolsym)then
2439 : 2061 : vin1(:)=vectc(:)
2440 [ + + ]: 8244 : vin2(:)=2.0d0*(vecta(:)-reduceda*vectc(:))
2441 : 2061 : vext(:)=vectb(:)
2442 : : found=1
2443 [ + + ]: 25979 : else if(abs(abs(reducedb)-0.5d0)<tolsym)then
2444 : 726 : vin1(:)=vectc(:)
2445 [ + + ]: 2904 : vin2(:)=2.0d0*(vectb(:)-reduceda*vectc(:))
2446 : 726 : vext(:)=vecta(:)
2447 : : found=1
2448 : : end if
2449 : 4627 : if(found==1)exit
2450 : : end do
2451 : : ! Now, vin1 and vin2 are perpendicular to each other, and in the plane that contains the binary vector.
2452 : : ! One of them must be the binary vector if any.
2453 : : ! On the other hand, vext is out-of-plane. Might belong to the mediator plane or not.
2454 : : ! If C monoclinc, then the projection of this vext on the binary vector will be either 0 or +1/2 or -1/2.
2455 : : ! The binary axis must be stored in cell_base(:,2) for conventional C-cell
2456 [ + + ]: 7414 : if(found==1)then
2457 : 2787 : found=0
2458 : :
2459 : : ! Test vin1 being the binary axis
2460 : 2787 : norm2trial=vin1(1)**2+vin1(2)**2+vin1(3)**2
2461 : 2787 : reduceda=(vext(1)*vin1(1)+vext(2)*vin1(2)+vext(3)*vin1(3))/norm2trial
2462 [ + + ]: 2787 : if(abs(reduceda)<tolsym)then ! vin1 is the binary axis and vext is in the mediator plane
2463 : 5052 : found=1
2464 [ + + ]: 5052 : cell_base(:,1)=vin2(:)
2465 [ + + ]: 5052 : cell_base(:,2)=vin1(:)
2466 [ + + ]: 5052 : cell_base(:,3)=vext(:)
2467 [ + + ]: 1524 : else if(abs(abs(reduceda)-0.5d0)<tolsym)then ! vin1 is the binary axis and vext has +1/2 or -1/2 as projection
2468 : 160 : found=1
2469 [ + + ]: 160 : cell_base(:,1)=vin2(:)
2470 [ + + ]: 160 : cell_base(:,2)=vin1(:)
2471 [ + + ]: 160 : cell_base(:,3)=vext(:)-reduceda*vin1(:)+vin2(:)*half
2472 : : else
2473 : : ! Test vin2 being the binary axis
2474 : 1484 : norm2trial=vin2(1)**2+vin2(2)**2+vin2(3)**2
2475 : 1484 : reduceda=(vext(1)*vin2(1)+vext(2)*vin2(2)+vext(3)*vin2(3))/norm2trial
2476 [ + + ]: 1484 : if(abs(reduceda)<tolsym)then ! vin2 is the binary axis and vext is in the mediator plane
2477 : 328 : found=1
2478 [ + + ]: 328 : cell_base(:,1)=vin1(:)
2479 [ + + ]: 328 : cell_base(:,2)=vin2(:)
2480 [ + + ]: 328 : cell_base(:,3)=vext(:)
2481 [ - + ]: 1402 : else if(abs(abs(reduceda)-0.5d0)<tolsym)then ! vin2 is the binary axis and vext has +1/2 or -1/2 as projection
2482 : 0 : found=1
2483 [ # # ]: 0 : cell_base(:,1)=vin1(:)
2484 [ # # ]: 0 : cell_base(:,2)=vin2(:)
2485 [ # # ]: 0 : cell_base(:,3)=vext(:)-reduceda*vin2(:)+vin1(:)*half
2486 : : end if
2487 : : end if
2488 : :
2489 : : if(found==1)then
2490 : 1385 : fact=2 ; center=3
2491 : 1385 : call holocell(cell_base,0,foundc,iholohedry,tolsym)
2492 : : end if
2493 : : end if
2494 : : end if
2495 : :
2496 : : end do ! Do-loop on three different directions
2497 : : end do ! Do-loop on different target holohedries
2498 : :
2499 [ + + ]: 40291 : if(foundc==0)then
2500 : 1643 : iholohedry=1 ; fact=1 ; center=0
2501 : 1643 : cell_base(:,:)=minim(:,:)
2502 : : end if
2503 : :
2504 : : !DEBUG
2505 : : !write(std_out,*)' symlatt : done with centering tests, foundc=',foundc
2506 : : !write(std_out,*)' center=',center
2507 : : !write(std_out,*)' iholohedry=',iholohedry
2508 : : !call flush(std_out)
2509 : : !ENDDEBUG
2510 : :
2511 : : !--------------------------------------------------------------------------
2512 : : !Final check on the Bravais lattice, using the basis vectors
2513 : :
2514 : : !Recompute the metric tensor
2515 [ + + ]: 40291 : if(foundc==1)then
2516 [ + + ]: 154592 : do ii=1,3
2517 : : metmin(:,ii)=cell_base(1,:)*cell_base(1,ii)+&
2518 : : & cell_base(2,:)*cell_base(2,ii)+&
2519 [ + + ]: 502424 : & cell_base(3,:)*cell_base(3,ii)
2520 : : end do
2521 : : end if
2522 : :
2523 : : !Examine the angles and vector lengths
2524 : 40291 : ang90(:)=0
2525 [ + + ]: 40291 : if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1
2526 [ + + ]: 40291 : if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1
2527 [ + + ]: 40291 : if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1
2528 : 40291 : equal(:)=0
2529 [ + + ]: 40291 : if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1
2530 [ + + ]: 40291 : if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1
2531 [ + + ]: 40291 : if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1
2532 : :
2533 : : !DEBUG
2534 : : !write(std_out,*)' symlatt : recompute the metric tensor '
2535 : : !write(std_out,*)' ang90=',ang90
2536 : : !write(std_out,*)' equal=',equal
2537 : : !call flush(std_out)
2538 : : !ENDDEBUG
2539 : :
2540 : : !The axes will be aligned with the previously determined
2541 : : !basis vectors, EXCEPT for the tetragonal cell, see later
2542 : 40291 : axes(:,:)=cell_base(:,:)
2543 : :
2544 : 40291 : found=0
2545 : : !Check orthogonal conventional cells
2546 [ + + ]: 40291 : if(ang90(1)+ang90(2)+ang90(3)==3)then
2547 : :
2548 : : ! Cubic system
2549 [ + + ]: 32724 : if(equal(1)+equal(2)+equal(3)==3)then
2550 : : ! However, one-face centered is not admitted
2551 [ + + + - ]: 24684 : if(center==0 .or. center==-1 .or. center==-3)then
2552 : 16739 : iholohedry=7 ; found=1
2553 [ + + ]: 7945 : if(center==0)then
2554 : 6754 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cP (primitive cubic)'
2555 [ + + ]: 17930 : else if(center==-1)then
2556 : 1191 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cI (body-centered cubic)'
2557 : : else if(center==-3)then
2558 : 16739 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cF (face-centered cubic)'
2559 : : end if
2560 : : end if
2561 : : end if
2562 : :
2563 : : ! Tetragonal system
2564 [ + + + + : 8040 : if(found==0 .and. &
+ + ]
2565 : : & (equal(1)==1 .or. equal(2)==1 .or. equal(3)==1) )then
2566 : : ! However, one-face centered or face-centered is not admitted
2567 [ + + ]: 5439 : if(center==0 .or. center==-1)then
2568 : 5361 : iholohedry=4 ; found=1
2569 [ + + ]: 5361 : if(equal(1)==1)then
2570 [ + + + + : 24560 : axes(:,3)=cell_base(:,1) ; axes(:,1)=cell_base(:,2) ; axes(:,2)=cell_base(:,3)
+ + ]
2571 [ + + ]: 2905 : else if(equal(2)==1)then
2572 [ + + + + : 1340 : axes(:,3)=cell_base(:,2) ; axes(:,2)=cell_base(:,1) ; axes(:,1)=cell_base(:,3)
+ + ]
2573 [ + - ]: 2771 : else if(equal(3)==1)then
2574 : 2771 : axes(:,:)=cell_base(:,:)
2575 : : end if
2576 [ + + ]: 5361 : if(center==0)then
2577 : 4747 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is tP (primitive tetragonal)'
2578 : : else if(center==-1)then
2579 : 614 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is tI (body-centered tetragonal)'
2580 : : end if
2581 : : end if
2582 : : end if
2583 : :
2584 : : ! Orthorhombic system
2585 : : if(found==0)then
2586 : 2679 : iholohedry=3 ; found=1
2587 : 2679 : axes(:,:)=cell_base(:,:)
2588 [ + + ]: 2679 : if(center==0)then
2589 : 1751 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oP (primitive orthorhombic)'
2590 [ + + ]: 928 : else if(center==-1)then
2591 : 302 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oI (body-centered orthorhombic)'
2592 [ + + ]: 626 : else if(center==1 .or. center==2 .or. center==3)then
2593 : 526 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oC (one-face-centered orthorhombic)'
2594 [ + - ]: 100 : else if(center==-3)then
2595 : 100 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oF (face-centered orthorhombic)'
2596 : : end if
2597 : : end if
2598 : :
2599 : : else
2600 : :
2601 : : ! Hexagonal system
2602 [ + + + + : 7567 : 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
+ - + - ]
2603 : 2196 : iholohedry=6 ; found=1
2604 : 2196 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is hP (primitive hexagonal)'
2605 : : end if
2606 : :
2607 : : ! Rhombohedral system
2608 : : if(found==0 .and. equal(1)+equal(2)+equal(3)==3 .and. &
2609 [ + + + - : 5371 : & abs(metmin(2,1)-metmin(3,2))<tolsym*metmin(2,2) .and. &
+ - ]
2610 : : & abs(metmin(2,1)-metmin(3,1))<tolsym*metmin(1,1) )then
2611 : 1581 : iholohedry=5 ; found=1
2612 : 1581 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is hR (rhombohedral)'
2613 : : end if
2614 : :
2615 : : ! Monoclinic system
2616 [ + + ]: 3790 : if(found==0 .and. ang90(1)+ang90(2)+ang90(3)==2 )then
2617 : 2147 : iholohedry=2 ; found=1
2618 [ + + ]: 2147 : if(center==0)then
2619 : 299 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is mP (primitive monoclinic)'
2620 [ + - ]: 1848 : else if(center==3)then
2621 : 1848 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is mC (one-face-centered monoclinic)'
2622 : : end if
2623 : : end if
2624 : :
2625 : : ! Triclinic system
2626 : : if(found==0)then
2627 : 1643 : iholohedry=1 ; found=1
2628 : 1643 : write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is aP (primitive triclinic)'
2629 : : end if
2630 : :
2631 : : end if
2632 : :
2633 : 40291 : call wrtout(iout,msg)
2634 : :
2635 : : !DEBUG
2636 : : !write(std_out,*)' symlatt : after checking conventional orthogonal cell '
2637 : : !call flush(std_out)
2638 : : !ENDDEBUG
2639 : :
2640 : : !--------------------------------------------------------------------------
2641 : : !Make sure that axes form a right-handed coordinate system
2642 : : !(Note : this should be done in the body of the routine,
2643 : : !by making changes that leave the sign of the mixed product of the three
2644 : : !vectors invariant)
2645 : : determinant=axes(1,1)*axes(2,2)*axes(3,3) &
2646 : : & +axes(1,2)*axes(2,3)*axes(3,1) &
2647 : : & +axes(1,3)*axes(3,2)*axes(2,1) &
2648 : : & -axes(1,1)*axes(3,2)*axes(2,3) &
2649 : : & -axes(1,3)*axes(2,2)*axes(3,1) &
2650 : 40291 : & -axes(1,2)*axes(2,1)*axes(3,3)
2651 [ + + ]: 40291 : if(determinant<zero)then
2652 [ + + + + ]: 71136 : axes(:,:)=-axes(:,:)
2653 : : end if
2654 : :
2655 : : !DEBUG
2656 : : !write(std_out,'(a,i4)')' symlatt : before itrial do loop, iholohedry= ',iholohedry
2657 : : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',&
2658 : : !& rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
2659 : : !call flush(std_out)
2660 : : !ENDDEBUG
2661 : :
2662 : : !--------------------------------------------------------------------------
2663 : : !Prefer symmetry axes on the same side as the primitive axes,
2664 : : !when the changes are allowed
2665 [ + - ]: 63889 : do itrial=1,100
2666 : :
2667 : : ! DEBUG
2668 : : ! write(std_out,'(a)')' '
2669 : : ! write(std_out,'(a,i5)')' symlatt : itrial do loop, itrial= ',itrial
2670 : : ! write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',&
2671 : : ! & axes(:,1),ch10,axes(:,2),ch10,axes(:,3)
2672 : : ! call flush(std_out)
2673 : : ! ENDDEBUG
2674 : :
2675 [ + + ]: 255556 : do ia=1,3
2676 : : scprods(ia,:)=axes(1,ia)*rprimd(1,:)+&
2677 : : & axes(2,ia)*rprimd(2,:)+&
2678 [ + + ]: 766668 : & axes(3,ia)*rprimd(3,:)
2679 [ + + ]: 766668 : norm2trial=sum(axes(:,ia)**2)
2680 [ + + ]: 830557 : scprods(ia,:)=scprods(ia,:)/sqrt(norm2trial)
2681 : : end do
2682 [ + + ]: 255556 : do ia=1,3
2683 [ + + ]: 766668 : norm2trial=sum(rprimd(:,ia)**2)
2684 [ + + ]: 830557 : scprods(:,ia)=scprods(:,ia)/sqrt(norm2trial)
2685 : : end do
2686 : :
2687 : : !DEBUG
2688 : : ! write(std_out,'(a,3f12.6)')' diagonal scalar products ',scprods(1,1),scprods(2,2),scprods(3,3)
2689 : : ! call flush(std_out)
2690 : : !ENDDEBUG
2691 : :
2692 : : ! One should now try all the generators of the
2693 : : ! proper rotations of each Bravais lattice, coupled with change of
2694 : : ! signs of each vector. This is not done systematically in what follows ...
2695 : : ! Here, the third axis is left unchanged
2696 [ + + ]: 63889 : if(iholohedry/=5)then
2697 [ + + + + ]: 59487 : if(scprods(1,1)<-tolsym .and. scprods(2,2)<-tolsym)then
2698 [ + + + + ]: 16450 : axes(:,1)=-axes(:,1) ; axes(:,2)=-axes(:,2)
2699 : : cycle
2700 : : end if
2701 : : end if
2702 : : ! The first (or second) axis is left unchanged
2703 [ + + ]: 61539 : if(iholohedry/=5 .and. iholohedry/=6)then
2704 [ + + + + ]: 54941 : if(scprods(2,2)<-tolsym .and. scprods(3,3)<-tolsym)then
2705 [ + + + + ]: 9079 : axes(:,2)=-axes(:,2) ; axes(:,3)=-axes(:,3)
2706 : : cycle
2707 : : end if
2708 [ + + + + ]: 53644 : if(scprods(1,1)<-tolsym .and. scprods(3,3)<-tolsym)then
2709 [ + + + + ]: 3318 : axes(:,1)=-axes(:,1) ; axes(:,3)=-axes(:,3)
2710 : : cycle
2711 : : end if
2712 : : end if
2713 : : ! Permutation of the three axis
2714 [ + + ]: 59768 : if(iholohedry==5 .or. iholohedry==7)then
2715 : 45742 : trace=scprods(1,1)+scprods(2,2)+scprods(3,3)
2716 [ + + ]: 45742 : if(trace+tolsym< scprods(1,2)+scprods(2,3)+scprods(3,1))then
2717 [ + + + + ]: 3591 : vecta(:)=axes(:,1) ; axes(:,1)=axes(:,3)
2718 [ + + + + ]: 3591 : axes(:,3)=axes(:,2); axes(:,2)=vecta(:)
2719 : : cycle
2720 : : end if
2721 [ + + ]: 45229 : if(trace+tolsym < scprods(1,3)+scprods(2,1)+scprods(3,2))then
2722 [ + + + + ]: 123648 : vecta(:)=axes(:,1) ; axes(:,1)=axes(:,2)
2723 [ + + + + ]: 123648 : axes(:,2)=axes(:,3); axes(:,3)=vecta(:)
2724 : : cycle
2725 : : end if
2726 : : ! This case is observed when the three new vectors
2727 : : ! are pointing opposite to the three original vectors
2728 : : ! One takes their opposite, then switch two of them, then process
2729 : : ! them again in the loop
2730 [ + + + + : 358345 : if(sum(scprods(:,:))<-tolsym)then
+ + ]
2731 [ + + ]: 5200 : axes(:,1)=-axes(:,1)
2732 [ + + ]: 5200 : vecta(:)=-axes(:,2)
2733 [ + + ]: 5200 : axes(:,2)=-axes(:,3)
2734 [ + + ]: 5200 : axes(:,3)=vecta(:)
2735 : : cycle
2736 : : end if
2737 : : end if
2738 : :
2739 : : ! Actually, for iholohedry==7, can test specifically all possibilities
2740 : : ! and take the best one.
2741 : : ! Not activated, because changing the order of symmetries in many tests !
2742 : : if(iholohedry==7 .and. .false.)then
2743 : : ! if(iholohedry==7)then
2744 : :
2745 : : !DEBUG
2746 : : !write(std_out,'(a,a)')ch10,' enter search of all possibilities for iholohedry==7 '
2747 : : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',&
2748 : : !& axes(:,1),ch10,axes(:,2),ch10,axes(:,3)
2749 : : !call flush(std_out)
2750 : : !ENDDEBUG
2751 : :
2752 : : do iaxis1=1,3
2753 : : axes_try(:,1)=axes(:,iaxis1)
2754 : : do iaxis2=1,2
2755 : : if(iaxis1==1)axes_try(:,2)=axes(:,1+iaxis2)
2756 : : if(iaxis1==2)axes_try(:,2)=axes(:,2*iaxis2-1)
2757 : : if(iaxis1==3)axes_try(:,2)=axes(:,iaxis2)
2758 : : if(iaxis2==1)axes_try(:,3)=axes(:,3)
2759 : : if(iaxis2==2)axes_try(:,3)=axes(:,1)
2760 : : if(iaxis1==1.and.iaxis2==2)axes_try(:,3)=axes(:,2)
2761 : : if(iaxis1==3.and.iaxis2==1)axes_try(:,3)=axes(:,2)
2762 : : do isign1=1,-1,-2
2763 : : axes_try(:,1)=-axes_try(:,1)
2764 : : do isign2=1,-1,-2
2765 : : axes_try(:,2)=-axes_try(:,2)
2766 : : determinant=axes_try(1,1)*axes_try(2,2)*axes_try(3,3) &
2767 : : & +axes_try(1,2)*axes_try(2,3)*axes_try(3,1) &
2768 : : & +axes_try(1,3)*axes_try(3,2)*axes_try(2,1) &
2769 : : & -axes_try(1,1)*axes_try(3,2)*axes_try(2,3) &
2770 : : & -axes_try(1,3)*axes_try(2,2)*axes_try(3,1) &
2771 : : & -axes_try(1,2)*axes_try(2,1)*axes_try(3,3)
2772 : : if(determinant<zero)axes_try(:,3)=-axes_try(:,3)
2773 : : do ia=1,3
2774 : : scprods(ia,:)=axes_try(1,ia)*rprimd(1,:)+&
2775 : : & axes_try(2,ia)*rprimd(2,:)+&
2776 : : & axes_try(3,ia)*rprimd(3,:)
2777 : : norm2trial=sum(axes_try(:,ia)**2)
2778 : : scprods(ia,:)=scprods(ia,:)/sqrt(norm2trial)
2779 : : end do
2780 : : do ia=1,3
2781 : : norm2trial=sum(rprimd(:,ia)**2)
2782 : : scprods(:,ia)=scprods(:,ia)/sqrt(norm2trial)
2783 : : end do
2784 : : trace=scprods(1,1)+scprods(2,2)+scprods(3,3)
2785 : : if(iaxis1==1.and.iaxis2==1.and.isign1==1.and.isign2==1)then
2786 : : trace_best=trace
2787 : : axes_best=axes_try
2788 : : else if (trace>trace_best+tolsym)then
2789 : : trace_best=trace
2790 : : axes_best=axes_try
2791 : : endif
2792 : : enddo ! isign2
2793 : : enddo ! isign1
2794 : : enddo ! iaxes2
2795 : : enddo ! iaxes1
2796 : : axes=axes_best
2797 : : endif ! iholohedry=7
2798 : 0 : exit
2799 : : end do
2800 : :
2801 : : !--------------------------------------------------------------------------
2802 : :
2803 : : !DEBUG
2804 : : !write(std_out,'(a,a)')ch10,' after order/sign optimization do-loop '
2805 : : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',&
2806 : : !& rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3)
2807 : : !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',&
2808 : : !& axes(:,1),ch10,axes(:,2),ch10,axes(:,3)
2809 : : !call flush(std_out)
2810 : : !ENDDEBUG
2811 : :
2812 : : !Compute the coordinates of rprimd in the system defined by axes(:,:)
2813 : 40291 : call matr3inv(axes,axesinvt)
2814 [ + + ]: 161164 : do ii=1,3
2815 : : coord(:,ii)=rprimd(1,ii)*axesinvt(1,:)+ &
2816 : : & rprimd(2,ii)*axesinvt(2,:)+ &
2817 [ + + ]: 523783 : & rprimd(3,ii)*axesinvt(3,:)
2818 : : end do
2819 : :
2820 : : !Check that the coordinates are integers, or half-integer in
2821 : : !the case there is a centering, and generate integer coordinates
2822 [ + + ]: 161164 : do ii=1,3
2823 [ + + ]: 523783 : do jj=1,3
2824 : 362619 : val=coord(ii,jj)*fact
2825 [ - + ]: 362619 : if(abs(val-nint(val))>fact*two*tolsym)then
2826 : : write(msg,'(4a,a,3es18.10,a,a,3es18.10,a,a,3es18.10,a,a,i4)')&
2827 : 0 : 'One of the coordinates of rprimd in axes is non-integer,',ch10,&
2828 : 0 : 'or non-half-integer (if centering), within 2*tolsym.',ch10,&
2829 : 0 : 'coord=',coord(:,1),ch10,&
2830 : 0 : ' ',coord(:,2),ch10,&
2831 : 0 : ' ',coord(:,3),ch10,&
2832 : 0 : 'fact=',fact
2833 : 0 : ABI_BUG(msg)
2834 : : end if
2835 : 483492 : icoord(ii,jj)=nint(val)
2836 : : end do
2837 : : end do
2838 : :
2839 : : !Store the bravais lattice characteristics
2840 : 40291 : bravais(1)=iholohedry
2841 : 40291 : bravais(2)=center
2842 [ + + ]: 161164 : bravais(3:5)=icoord(1:3,1)
2843 [ + + ]: 161164 : bravais(6:8)=icoord(1:3,2)
2844 [ + + ]: 161164 : bravais(9:11)=icoord(1:3,3)
2845 : :
2846 : : !--------------------------------------------------------------------------
2847 : : !Initialize the set of symmetries
2848 : : !Bravais lattices are always invariant under identity and inversion
2849 : :
2850 : : !Identity and inversion
2851 [ + + + + : 1007275 : ptsymrel(:,:,1)=identity(:,:) ; ptsymrel(:,:,2)=-identity(:,:)
+ + + + ]
2852 : 40291 : nptsym=2
2853 : :
2854 : : !Keep this for IFCv70 compiler
2855 : : if(nptsym/=2)then
2856 : : write(msg,'(a,a,a,a)')ch10,&
2857 : : ' symlatt : BUG -',ch10,&
2858 : : ' Crazy error, compiler bug '
2859 : : call wrtout(std_out,msg)
2860 : : end if
2861 : :
2862 : : !--------------------------------------------------------------------------
2863 : : !Initialize some generators
2864 : : !gen6 is defined in a coordinated system with gamma=120 degrees
2865 : 40291 : gen6(:,:)=0 ; gen6(3,3)=1 ; gen6(1,1)=1 ; gen6(1,2)=-1 ; gen6(2,1)=1
2866 : 40291 : gen3(:,:)=0 ; gen3(1,2)=1 ; gen3(2,3)=1 ; gen3(3,1)=1
2867 : 40291 : gen2xy(:,:)=0 ; gen2xy(2,1)=1 ; gen2xy(1,2)=1; gen2xy(3,3)=1
2868 : 40291 : gen2y(:,:)=0 ; gen2y(1,1)=-1; gen2y(2,2)=1 ; gen2y(3,3)=-1
2869 : 40291 : gen2z(:,:)=0 ; gen2z(1,1)=-1; gen2z(2,2)=-1; gen2z(3,3)=1
2870 : :
2871 : : !--------------------------------------------------------------------------
2872 : :
2873 : : !Define the generators for each holohedry (inversion is already included)
2874 [ + + ]: 40291 : if(iholohedry==6)then
2875 : 2196 : ngen=2
2876 [ + + + + ]: 28548 : gen(:,:,1)=gen2xy(:,:) ; order(1)=2
2877 [ + + + + ]: 28548 : gen(:,:,2)=gen6(:,:) ; order(2)=6
2878 [ + + ]: 38095 : else if(iholohedry==5)then
2879 : 1581 : ngen=2
2880 [ + + + + ]: 20553 : gen(:,:,1)=gen2xy(:,:) ; order(1)=2
2881 [ + + + + ]: 20553 : gen(:,:,2)=gen3(:,:) ; order(2)=3
2882 : : else
2883 [ + + + + ]: 474682 : gen(:,:,1)=gen2y(:,:) ; order(1)=2
2884 [ + + + + ]: 474682 : gen(:,:,2)=gen2z(:,:) ; order(2)=2
2885 [ + + + + ]: 474682 : gen(:,:,3)=gen2xy(:,:) ; order(3)=2
2886 [ + + + + ]: 474682 : gen(:,:,4)=gen3(:,:) ; order(4)=3
2887 [ + + ]: 36514 : if(iholohedry<=4)ngen=iholohedry-1
2888 [ + + ]: 36514 : if(iholohedry==7)ngen=4
2889 : : end if
2890 : :
2891 : : !Build the point symmetry operations from generators, in the reduced system
2892 : : !of coordinates defined by axes(:,:)
2893 [ + + ]: 40291 : if(ngen/=0)then
2894 [ + + ]: 168526 : do igen=1,ngen
2895 [ + + ]: 1424886 : do isym=1+nptsym,order(igen)*nptsym
2896 : 1295008 : jsym=isym-nptsym
2897 [ + + ]: 5309910 : do ii=1,3
2898 : : ptsymrel(:,ii,isym)=gen(:,1,igen)*ptsymrel(1,ii,jsym)+ &
2899 : : & gen(:,2,igen)*ptsymrel(2,ii,jsym)+ &
2900 [ + + ]: 16835104 : & gen(:,3,igen)*ptsymrel(3,ii,jsym)
2901 : : end do
2902 : : end do
2903 : 168526 : nptsym=order(igen)*nptsym
2904 : :
2905 : : end do
2906 : : end if
2907 : :
2908 : : !--------------------------------------------------------------------------
2909 : :
2910 : : !Transform symmetry matrices in the system defined by rprimd
2911 : 40291 : call symrelrot(nptsym,axes,rprimd,ptsymrel,tolsym)
2912 : :
2913 : : !DEBUG
2914 : : !write(std_out,'(a)') ' symlatt : exit '
2915 : : !call flush(std_out)
2916 : : !stop
2917 : : !ENDDEBUG
2918 : :
2919 : 40291 : end subroutine symlatt
2920 : : !!***
2921 : :
2922 [ - + ]: 206258 : end module m_symfind
2923 : : !!***
|
