Branch data Line data Source code
1 : :
2 : : #if defined HAVE_CONFIG_H
3 : : #include "config.h"
4 : : #endif
5 : :
6 : : #include "abi_common.h"
7 : :
8 : : module m_tdep_sym
9 : :
10 : : use defs_basis
11 : : use m_abicore
12 : : use m_errors
13 : : use m_xmpi
14 : : use m_matrix, only : mati3inv
15 : : use m_symtk, only : symatm
16 : : use m_symfind, only : symfind, symanal, symlatt
17 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
18 : : use m_tdep_latt, only : Lattice_type, tdep_make_inbox
19 : :
20 : : implicit none
21 : :
22 : : type,public :: Symmetries_type
23 : :
24 : : integer :: msym
25 : : integer :: nptsym
26 : : integer :: nsym
27 : : integer :: spgroup
28 : : integer, allocatable :: ptsymrel(:,:,:)
29 : : integer, allocatable :: symrec(:,:,:)
30 : : integer, allocatable :: symrel(:,:,:)
31 : : integer, allocatable :: symafm(:)
32 : : integer, allocatable :: indsym(:,:,:)
33 : : double precision, allocatable :: S_ref(:,:,:,:)
34 : : double precision, allocatable :: S_inv(:,:,:,:)
35 : : double precision, allocatable :: tnons(:,:)
36 : : double precision, allocatable :: xred_zero(:,:)
37 : :
38 : : end type Symmetries_type
39 : :
40 : : public :: tdep_make_sym
41 : : public :: tdep_SearchS_1at
42 : : public :: tdep_SearchS_2at
43 : : public :: tdep_SearchS_3at
44 : : public :: tdep_SearchS_4at
45 : : public :: tdep_destroy_sym
46 : : public :: tdep_calc_indsym2
47 : :
48 : : contains
49 : :
50 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
51 [ - + - + : 44 : subroutine tdep_make_sym(Invar,Lattice,MPIdata,Sym)
- + - + -
+ - + - +
- + ]
52 : :
53 : : type(Symmetries_type),intent(out) :: Sym
54 : : type(atdep_dataset_type),intent(inout) :: Invar
55 : : type(Lattice_type),intent(in) :: Lattice
56 : : type(MPI_enreg_type), intent(in) :: MPIdata
57 : :
58 : : integer :: nspden,use_inversion,chkprim
59 : : integer :: ptgroupma,isym,ii,jj,iatom_unitcell
60 : : integer :: bravais(11)
61 : 44 : integer, allocatable :: symrel(:,:,:),symafm(:)
62 : : double precision :: genafm(3)
63 : : double precision :: temp1(3,1), temp2(3,1)
64 : 44 : double precision, allocatable :: spinat(:,:)
65 : 44 : double precision, allocatable :: tnons(:,:)
66 [ - + ]: 44 : double precision, allocatable :: tmp1(:,:)
67 : :
68 : :
69 : : ! Compute all the symetries coming from the bravais lattice
70 : : ! The routine used is symlatt (from Abinit code)
71 : : ! GA: FIXME msym needs to be an input variable
72 : 44 : Sym%msym=1000 !msym needs to be very large due to non-primitive cell calculations
73 [ - + + + : 572044 : ABI_MALLOC(Sym%ptsymrel,(3,3,Sym%msym)) ; Sym%ptsymrel(:,:,:)=0
+ + + + ]
74 : 44 : call symlatt(Invar%bravais,std_out,Sym%msym,Sym%nptsym,Sym%ptsymrel,Lattice%rprimdt,tol8)
75 : 44 : write(Invar%stdout,'(a,1x,11(i4,1x))')' bravais=',Invar%bravais(:)
76 : 44 : Sym%nsym=Sym%nptsym
77 : :
78 : : ! Initialize all the symetries using the first atom at (0.0;0.0;0.0)
79 [ - + + - : 556 : ABI_MALLOC(Sym%xred_zero,(3,Invar%natom_unitcell)) ; Sym%xred_zero(:,:)=0.d0
+ - + - -
+ - + - +
+ + + + ]
80 [ + + ]: 150 : do iatom_unitcell=1,Invar%natom_unitcell
81 [ + + ]: 424 : temp1(:,1)=Invar%xred_unitcell(:,iatom_unitcell)-Invar%xred_unitcell(:,1)
82 [ + + ]: 424 : do ii=1,3
83 [ + + ]: 318 : if (temp1(ii,1).lt.0.d0) then
84 : 66 : jj=int(temp1(ii,1)+1.0)
85 [ + - ]: 252 : else if (temp1(ii,1).ge.0.d0) then
86 : 252 : jj=int(temp1(ii,1)+0.0)
87 : : end if
88 : 424 : temp2(ii,1)=temp1(ii,1)-real(jj)
89 : : end do
90 [ + + ]: 468 : Sym%xred_zero(:,iatom_unitcell)=temp2(:,1)
91 : : end do
92 : :
93 : : ! Calculation of the (non-symmorphic) translations
94 : 44 : nspden=1
95 [ - + + - : 512 : ABI_MALLOC(spinat,(3,Invar%natom_unitcell)); spinat(:,:)=0.d0
- + - + +
+ + + ]
96 : 44 : use_inversion=1
97 [ - + + - : 572132 : ABI_MALLOC(symrel ,(3,3,Sym%msym)) ; symrel (:,:,:)=0
+ - + - -
+ - + + +
+ + + + ]
98 [ - + + - : 176132 : ABI_MALLOC(tnons ,(3,Sym%msym)) ; tnons (:,:)=0.d0
+ - + - -
+ - + + +
+ + ]
99 [ - + + - : 44132 : ABI_MALLOC(symafm ,(Sym%msym)) ; symafm (:)=1
+ - + - -
+ - + +
+ ]
100 : : call symfind(Lattice%gprimd,Sym%msym,Invar%natom_unitcell,Sym%nptsym,nspden,Sym%nsym,&
101 : 44 : & 0,Sym%ptsymrel,spinat,symafm,symrel,tnons,tol8,Invar%typat_unitcell,use_inversion,Sym%xred_zero)
102 [ - + ]: 44 : ABI_FREE(spinat)
103 [ - + + - : 21218 : ABI_MALLOC(Sym%symrel,(3,3,Sym%nsym)) ; Sym%symrel(:,:,:)=0
+ - + - -
+ - + - +
+ + + + +
+ ]
104 [ - + + - : 6620 : ABI_MALLOC(Sym%tnons ,(3,Sym%nsym)) ; Sym%tnons (:,:)=0.d0
+ - + - -
+ - + - +
+ + + + ]
105 [ - + + - : 1754 : ABI_MALLOC(Sym%symafm ,(Sym%nsym)) ; Sym%symafm (:)=1
+ - + - -
+ - + - +
+ + ]
106 [ + + + + : 21130 : Sym%symrel(:,:,:)=symrel(:,:,1:Sym%nsym)
+ + ]
107 [ + + + + ]: 6532 : Sym%tnons (:,:)=tnons (:,1:Sym%nsym)
108 [ + + ]: 1666 : Sym%symafm (:)=symafm (1:Sym%nsym)
109 [ - + ]: 44 : ABI_FREE(symrel)
110 [ - + ]: 44 : ABI_FREE(tnons)
111 [ - + ]: 44 : ABI_FREE(symafm)
112 [ + + ]: 44 : if (Sym%nptsym.ne.Sym%nsym) then
113 : 3 : write(Invar%stdlog,'(a,i4,a,i4)') 'WARNING: nsym=',Sym%nsym,' is not equal to nptsym=',Sym%nptsym
114 : 3 : write(Invar%stdlog,'(a)') ' Symrel='
115 [ + + ]: 293 : do isym=1,Sym%nsym
116 : 290 : write(Invar%stdlog,*) ' For sym=',isym
117 : 290 : write(Invar%stdlog,*) Sym%symrel(1,:,isym)
118 : 290 : write(Invar%stdlog,*) Sym%symrel(2,:,isym)
119 : 293 : write(Invar%stdlog,*) Sym%symrel(3,:,isym)
120 : : end do
121 : 3 : write(Invar%stdlog,'(a)') ' Ptsymrel='
122 [ + + ]: 111 : do isym=1,Sym%nptsym
123 : 108 : write(Invar%stdlog,*) ' For sym=',isym
124 : 108 : write(Invar%stdlog,*) Sym%ptsymrel(1,:,isym)
125 : 108 : write(Invar%stdlog,*) Sym%ptsymrel(2,:,isym)
126 : 111 : write(Invar%stdlog,*) Sym%ptsymrel(3,:,isym)
127 : : end do
128 : : end if
129 : :
130 : : ! Calculation of symrec
131 [ - + + - : 21218 : ABI_MALLOC(Sym%symrec,(3,3,Sym%nsym)); Sym%symrec(:,:,:)=0
+ - + - -
+ - + - +
+ + + + +
+ ]
132 [ + + ]: 1666 : do isym=1,Sym%nsym
133 : 1666 : call mati3inv(Sym%symrel(:,:,isym),Sym%symrec(:,:,isym))
134 : : end do
135 : : ! Calculation of spgroup
136 : 44 : chkprim=0
137 : 44 : genafm(:)=0.d0
138 : : ptgroupma=0
139 : 44 : Sym%spgroup=0
140 : 44 : call symanal(bravais,chkprim,genafm,Sym%msym,Sym%nsym,ptgroupma,Lattice%rprimdt,Sym%spgroup,Sym%symafm,Sym%symrel,Sym%tnons,tol8)
141 : :
142 : : ! Transform S_ref in cartesian coordinates
143 [ - + + - : 21262 : ABI_MALLOC(Sym%S_ref,(3,3,Sym%nsym,2)) ; Sym%S_ref(:,:,:,1)=real(Sym%symrel(:,:,1:Sym%nsym))
+ - + - +
- - + - +
- + + + +
+ + + ]
144 [ - + + - : 21174 : ABI_MALLOC(Sym%S_inv,(3,3,Sym%nsym,2)) ; Sym%S_inv(:,:,:,1)=zero
- + - + -
+ + + + +
+ + ]
145 [ - + + + : 572 : ABI_MALLOC(tmp1,(3,3)); tmp1(:,:)=0.d0
+ + ]
146 [ + - ]: 44 : if (MPIdata%iam_master) open(unit=75,file=trim(Invar%output_prefix)//'_sym.dat')
147 [ + + ]: 1666 : do isym=1,Sym%nsym
148 [ + - ]: 1622 : if (MPIdata%iam_master) then
149 : 1622 : write(75,*) ' '
150 : 1622 : write(75,*) 'For isym=',isym
151 : 1622 : write(75,*) 'In reduced coordinates:'
152 : 1622 : write(75,'(3(f5.2,1x))') Sym%S_ref(1,1,isym,1),Sym%S_ref(1,2,isym,1),Sym%S_ref(1,3,isym,1)
153 : 1622 : write(75,'(3(f5.2,1x))') Sym%S_ref(2,1,isym,1),Sym%S_ref(2,2,isym,1),Sym%S_ref(2,3,isym,1)
154 : 1622 : write(75,'(3(f5.2,1x))') Sym%S_ref(3,1,isym,1),Sym%S_ref(3,2,isym,1),Sym%S_ref(3,3,isym,1)
155 : : end if
156 : 1622 : call DGEMM('N','N',3,3,3,1.d0,Lattice%rprimdt,3,Sym%S_ref(:,:,isym,1),3,0.d0,tmp1,3)
157 : 1622 : call DGEMM('N','N',3,3,3,1.d0,tmp1,3,Lattice%rprimdtm1,3,0.d0,Sym%S_ref(:,:,isym,1),3)
158 [ + + ]: 6488 : do ii=1,3
159 [ + + ]: 21086 : do jj=1,3
160 [ + + ]: 19464 : if (abs(Sym%S_ref(ii,jj,isym,1)).lt.tol8) Sym%S_ref(ii,jj,isym,1)=zero
161 : : end do
162 : : end do
163 : : ! Inversion of S_ref (which is equivalent to the transposition)
164 [ + + ]: 6488 : do ii=1,3
165 [ + + ]: 21086 : do jj=1,3
166 : 19464 : Sym%S_inv(ii,jj,isym,1)=Sym%S_ref(jj,ii,isym,1)
167 : : end do
168 : : end do
169 [ + - ]: 1666 : if (MPIdata%iam_master) then
170 : 1622 : write(75,*) 'In cartesian coordinates:'
171 : 1622 : write(75,'(3(f5.2,1x))') Sym%S_ref(1,1,isym,1),Sym%S_ref(1,2,isym,1),Sym%S_ref(1,3,isym,1)
172 : 1622 : write(75,'(3(f5.2,1x))') Sym%S_ref(2,1,isym,1),Sym%S_ref(2,2,isym,1),Sym%S_ref(2,3,isym,1)
173 : 1622 : write(75,'(3(f5.2,1x))') Sym%S_ref(3,1,isym,1),Sym%S_ref(3,2,isym,1),Sym%S_ref(3,3,isym,1)
174 : : end if
175 : : end do
176 [ + - ]: 44 : if (MPIdata%iam_master) close(75)
177 : 44 : write(Invar%stdout,'(a)') ' See the sym.dat file'
178 : :
179 : : ! We verify that S.S^T = Id
180 [ + + + + ]: 572 : tmp1(:,:)=zero
181 [ + + ]: 1666 : do isym=1,Sym%nsym
182 : 1622 : call DGEMM('T','N',3,3,3,1.d0,Sym%S_ref(:,:,isym,1),3,Sym%S_ref(:,:,isym,1),3,0.d0,tmp1,3)
183 [ + + ]: 6532 : do ii=1,3
184 [ + + ]: 21086 : do jj=1,3
185 [ + + + - : 19464 : if ((ii/=jj.and.abs(tmp1(ii,jj)).gt.tol8).or.(ii==jj.and.abs(tmp1(ii,jj)-1.d0).gt.tol8)) then
+ + - + ]
186 : 0 : write(Invar%stdout,'(a)') ' STOP : the matrix is not orthogonal : '
187 : 0 : write(Invar%stdout,'(a,1x,i3,1x,i3,1x,f15.10)') ' ii, jj, sym(ii,jj)=',ii,jj,tmp1(ii,jj)
188 : 0 : ABI_ERROR('The matrix is not orthogonal')
189 : : end if
190 : : end do
191 : : end do
192 : : end do
193 : :
194 [ - + ]: 44 : ABI_FREE(tmp1)
195 : :
196 : 44 : end subroutine tdep_make_sym
197 : :
198 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
199 : :
200 : 44 : subroutine tdep_SearchS_1at(Invar,MPIdata,Sym,xred_ideal)
201 : :
202 : : type(atdep_dataset_type), intent(in) :: Invar
203 : : type(Symmetries_type), intent(inout) :: Sym
204 : : type(MPI_enreg_type), intent(in) :: MPIdata
205 : : double precision, intent(in) :: xred_ideal(3,Invar%natom)
206 : :
207 : : integer :: isym,jatom,mu
208 : : integer :: iatom_unitcell,jatom_unitcell,iatom
209 : 88 : integer :: vecti(3),vectj(3),vectsym(4,Sym%nsym)
210 : 44 : integer, allocatable :: indsym2(:,:,:,:)
211 : : double precision :: temp3(3,1)
212 : 88 : double precision :: tmpi(3,Invar%natom),tmpj(3,Invar%natom),tmp_store(3,Invar%natom_unitcell)
213 : :
214 : 44 : write(Invar%stdout,*) ' '
215 : 44 : write(Invar%stdout,*) '#############################################################################'
216 : 44 : write(Invar%stdout,*) '###################### Find the symetry operations ##########################'
217 : 44 : write(Invar%stdout,*) '#################### (connecting the atoms together) ########################'
218 : 44 : write(Invar%stdout,*) '#############################################################################'
219 : :
220 : : ! write(Invar%stdout,'(a)') 'Begining SearchS'
221 : : ! TODO: Initialize the local variables at 0
222 : :
223 : : ! Compute the indsym fundamental quantity
224 : : ! === Obtain a list of rotated atoms ===
225 : : ! $ R^{-1} (xred(:,iat)-\tau) = xred(:,iat_sym) + R_0 $
226 : : ! * indsym(4, isym,iat) gives iat_sym in the original unit cell.
227 : : ! * indsym(1:3,isym,iat) gives the lattice vector $R_0$.
228 : 44 : write(Invar%stdout,'(a)') ' Search the matrix transformation going from (k) to (i)...'
229 [ - + + - : 1207538 : ABI_MALLOC(Sym%indsym,(4,Sym%nsym,Invar%natom)); Sym%indsym(:,:,:)=zero
- + + - +
- + - - +
- + - + +
+ + + +
+ ]
230 : : call symatm(Sym%indsym(:,:,1:Invar%natom_unitcell),Invar%natom_unitcell,Sym%nsym,&
231 : 44 : & Sym%symrec,Sym%tnons,tol8,Invar%typat_unitcell,Sym%xred_zero)
232 : :
233 : : ! Store the positions of the atoms in the motif
234 [ + + ]: 150 : do iatom=1,Invar%natom_unitcell
235 : 150 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,xred_ideal(:,iatom),1,0.d0,tmp_store(:,iatom),1)
236 : : end do
237 : : ! Write the Indsym of the atoms included in the (reference) unitcell (i.e.: the motif)
238 [ + + + - ]: 44 : if (Invar%debug.and.MPIdata%iam_master) then
239 : 2 : open(unit=40,file=trim(Invar%output_prefix)//'_Indsym-unitcell.dat')
240 [ + + ]: 12 : do iatom=1,Invar%natom_unitcell
241 : 10 : write(40,*) '=========================================='
242 : 10 : write(40,'(a,i4,a,3(f10.5,1x))') 'For iatom=',iatom,' with xred (supercell)=',xred_ideal(:,iatom)
243 [ + + ]: 172 : do isym=1,Sym%nsym
244 : 160 : write(40,'(a,i2,a,i4,a,3(i4,1x),a,i2,a,3(f10.5,1x))') ' indsym(isym=',isym,',',iatom,')=',&
245 : 330 : & Sym%indsym(1:3,isym,iatom),'|iat=',Sym%indsym(4,isym,iatom),'| with tnons=',Sym%tnons(:,isym)
246 : : end do
247 : : end do
248 : 2 : close(40)
249 : : end if
250 : :
251 : : ! Search the matrix transformation going from (k,l) to (i,j)
252 : 44 : write(Invar%stdout,'(a)') ' Search the matrix transformation going from (k,l) to (i,j)...'
253 [ - + + - : 404752562 : ABI_MALLOC(indsym2,(8,Sym%nsym,Invar%natom,Invar%natom)); indsym2(:,:,:,:)=0
- + + - -
+ + - + -
+ - - + -
+ + + + +
+ + + + ]
254 [ + + + + ]: 27572 : tmpi(:,:)=0.d0
255 [ + + + + ]: 27572 : tmpj(:,:)=0.d0
256 [ + + ]: 6926 : do iatom=1,Invar%natom
257 : : ! For a single iatom
258 : 6882 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,xred_ideal(:,iatom),1,0.d0,tmpi(:,iatom),1)
259 : 6882 : iatom_unitcell=mod(iatom-1,Invar%natom_unitcell)+1
260 [ + + ]: 27528 : vecti(:)=nint(tmpi(:,iatom)-tmp_store(:,iatom_unitcell))
261 [ + + ]: 247022 : do isym=1,Sym%nsym
262 [ + + + + ]: 69595616 : vectsym(:,:)=0
263 [ + + ]: 960384 : do mu=1,3 ! Apply inverse transformation to original coordinates. Note transpose of symrec.
264 : 960384 : vectsym(mu,isym) = Sym%symrec(1,mu,isym)*vecti(1)+Sym%symrec(2,mu,isym)*vecti(2)+Sym%symrec(3,mu,isym)*vecti(3)
265 : : end do
266 [ + + ]: 1200480 : Sym%indsym(1:4,isym,iatom)=Sym%indsym(1:4,isym,iatom_unitcell)+vectsym(1:4,isym)
267 [ + + ]: 45055106 : do jatom=1,Invar%natom
268 [ + + ]: 224280736 : indsym2(1:4,isym,iatom,jatom)=Sym%indsym(1:4,isym,iatom_unitcell)+vectsym(1:4,isym)
269 : : end do
270 : : end do
271 : : end do
272 [ + + + - ]: 44 : if (Invar%debug.and.MPIdata%iam_master) then
273 : 2 : open(unit=40,file=trim(Invar%output_prefix)//'_Indsym-supercell.dat')
274 [ + + ]: 226 : do iatom=1,Invar%natom
275 : 224 : write(40,*) '=========================================='
276 : 224 : write(40,'(a,i4,a,3(f10.5,1x))') 'For iatom=',iatom,' with xred (supercell)=',xred_ideal(:,iatom)
277 [ + + ]: 7136 : do isym=1,Sym%nsym
278 : 6912 : write(40,'(a,i2,a,i4,a,3(i4,1x),a,i2,a,3(f10.5,1x))') ' indsym(isym=',isym,',',iatom,')=',&
279 : 14048 : & Sym%indsym(1:3,isym,iatom),'|iat=',Sym%indsym(4,isym,iatom),'| with tnons=',Sym%tnons(:,isym)
280 : : end do
281 : 226 : write(40,'(a,i4)') ' '
282 : : end do
283 : 2 : close(40)
284 : : end if
285 : :
286 : : ! For a couple of (iatom,jatom). The (iatom,jatom) vector depends on the position of iatom (due to PBC)
287 [ + + ]: 6926 : do iatom=1,Invar%natom
288 [ + + ]: 1479234 : do jatom=1,Invar%natom
289 [ + + ]: 5889232 : temp3(:,1)=xred_ideal(:,jatom)-xred_ideal(:,iatom)
290 : 1472308 : call tdep_make_inbox(temp3(:,1),1,1d-3,temp3(:,1))
291 [ + + ]: 5889232 : temp3(:,1)=xred_ideal(:,iatom)+temp3(:,1)
292 : 1472308 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,temp3(:,1),1,0.d0,tmpj(:,jatom),1)
293 : 1472308 : jatom_unitcell=mod(jatom-1,Invar%natom_unitcell)+1
294 [ + + ]: 5889232 : vectj(:)=nint(tmpj(:,jatom)-tmp_store(:,jatom_unitcell))
295 [ + + ]: 46287318 : do isym=1,Sym%nsym
296 [ + + + + ]: 11296166848 : vectsym(:,:)=0
297 [ + + ]: 179232512 : do mu=1,3 ! Apply inverse transformation to original coordinates. Note transpose of symrec.
298 : 179232512 : vectsym(mu,isym) = Sym%symrec(1,mu,isym)*vectj(1)+Sym%symrec(2,mu,isym)*vectj(2)+Sym%symrec(3,mu,isym)*vectj(3)
299 : : end do
300 [ + + ]: 225512948 : indsym2(5:8,isym,iatom,jatom)=Sym%indsym(1:4,isym,jatom_unitcell)+vectsym(1:4,isym)
301 : : end do
302 : : end do
303 : : end do
304 [ + + + - ]: 44 : if (Invar%debug.and.MPIdata%iam_master) then
305 : 2 : open(unit=40,file=trim(Invar%output_prefix)//'_Indsym-2atoms.dat')
306 [ + + ]: 226 : do iatom=1,Invar%natom
307 : 224 : write(40,*) '=========================================='
308 : 224 : write(40,'(a,i4,a,3(f10.5,1x))') 'For iatom=',iatom,' with xred (supercell)=',xred_ideal(:,iatom)
309 [ + + ]: 25824 : do jatom=1,Invar%natom
310 : 25600 : write(40,*) ' =========================================='
311 : 25600 : write(40,'(a,i4,a,3(f10.5,1x))') ' For jatom=',jatom,' with xred (supercell)=',xred_ideal(:,jatom)
312 [ + + ]: 885984 : do isym=1,Sym%nsym
313 : 860160 : write(40,'(a,i2,a,i4,a,i4,a,3(i4,1x),a,i2,a,3(i4,1x),a,i2,a)') ' indsym2(isym=',isym,',',iatom,',',jatom,')=',&
314 : 860160 : & indsym2(1:3,isym,iatom,jatom),&
315 : 1745920 : & '|iat=',indsym2(4,isym,iatom,jatom),'|',indsym2(5:7,isym,iatom,jatom),'|iat=',indsym2(8,isym,iatom,jatom),'|'
316 : : end do
317 : : end do
318 : 226 : write(40,'(a,i4)') ' '
319 : : end do
320 : 2 : close(40)
321 : : end if
322 [ - + ]: 44 : ABI_FREE(indsym2)
323 : 44 : write(Invar%stdout,'(a)') ' See the Indsym*.dat files (if debug)'
324 : :
325 : 44 : end subroutine tdep_SearchS_1at
326 : :
327 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
328 : 306801 : subroutine tdep_SearchS_2at(Invar,iatom,jatom,eatom,fatom,Isym2at,Sym,xred_ideal)
329 : :
330 : : integer, intent(in) :: iatom,jatom,eatom,fatom
331 : : type(atdep_dataset_type),intent(in) :: Invar
332 : : type(Symmetries_type),intent(in) :: Sym
333 : : integer, intent(inout) :: Isym2at(Invar%natom,Invar%natom,2)
334 : : double precision, intent(in) :: xred_ideal(3,Invar%natom)
335 : :
336 : : integer :: isym,ee,ff,ii
337 : : integer :: iatom_unitcell,jatom_unitcell
338 : : integer :: vecti(3),vectj(3),latte(3),lattf(3),indsym2(8)
339 : 613602 : double precision :: temp(3),tmp_store(3,Invar%natom_unitcell)
340 : 613602 : double precision :: tmp(3,Invar%natom)
341 : : logical :: ok
342 : :
343 : : ! Search the couple of atoms (e,f) obtained through the transformation (R+t)
344 : : ! starting from (i,j)
345 : : ! In that case, note that two cases are possible:
346 : : ! - The bond is just transformed, so indsym(4,e)=i and indsym(4,f)=j
347 : : ! - The bond is transformed and reversed, so indsym(4,e)=j and indsym(4,f)=i
348 : :
349 : : ! Store the positions of the atoms in the motif
350 [ + + ]: 931887 : do ii=1,Invar%natom_unitcell
351 : 931887 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,xred_ideal(:,ii),1,0.d0,tmp_store(:,ii),1)
352 : : end do
353 : :
354 : : ! Search the atom equivalent to iatom in the (reference) unitcell
355 [ + + ]: 77606881 : do ii=1,Invar%natom
356 : 77606881 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,xred_ideal(:,ii),1,0.d0,tmp(:,ii),1)
357 : : end do
358 : : ! Note that in the (present) particular case : iatom_unitcell=iatom and vecti(:)=zero
359 : 306801 : iatom_unitcell=mod(iatom-1,Invar%natom_unitcell)+1
360 [ + + ]: 1227204 : vecti(:)=nint(tmp(:,iatom)-tmp(:,iatom_unitcell))
361 : :
362 : : ! Search the atom equivalent to jatom in the (reference) unitcell.
363 : : ! jatom can be outside the box (according to the value of iatom)
364 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
365 [ + + ]: 1227204 : temp(:)=xred_ideal(:,jatom)-xred_ideal(:,iatom)
366 : 306801 : call tdep_make_inbox(temp,1,1d-3,temp)
367 [ + + ]: 1227204 : temp(:)=xred_ideal(:,iatom)+temp(:)
368 : 306801 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,temp(:),1,0.d0,tmp(:,jatom),1)
369 : 306801 : jatom_unitcell=mod(jatom-1,Invar%natom_unitcell)+1
370 [ + + ]: 1227204 : vectj(:)=nint(tmp(:,jatom)-tmp_store(:,jatom_unitcell))
371 : :
372 : : ! To understand the meaning of "latt", see SearchS_1at
373 : 306801 : ok=.false.
374 [ + + ]: 5232538 : do isym=1,Sym%nsym
375 : 5178841 : indsym2(:)=0
376 : 5178841 : call tdep_calc_indsym2(Invar,eatom,fatom,indsym2,isym,Sym,xred_ideal)
377 : 5178841 : ee=indsym2(4)
378 : 5178841 : ff=indsym2(8)
379 [ + + ]: 20715364 : latte(:)=indsym2(1:3)
380 [ + + ]: 20715364 : lattf(:)=indsym2(5:7)
381 [ + + ]: 5178841 : if (ee==iatom_unitcell.and.ff==jatom_unitcell) then
382 : : !FB if (Invar%debug) then
383 : : !FB write(Invar%stdout,*) 'For ee=iatom and ff=jatom, isym=',isym
384 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vecti(1),' - ',vectj(1),' - ',latte(1),' + ',lattf(1)
385 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vecti(2),' - ',vectj(2),' - ',latte(2),' + ',lattf(2)
386 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vecti(3),' - ',vectj(3),' - ',latte(3),' + ',lattf(3)
387 : : !FB end if
388 [ + + + + ]: 10628184 : if (sum(abs((vecti(:)-latte(:))-(vectj(:)-lattf(:)))).lt.tol8) then
389 : 253104 : Isym2at(eatom,fatom,1)=isym
390 : 253104 : Isym2at(eatom,fatom,2)=1
391 : : ok=.true.
392 : : end if
393 : : end if
394 : 53697 : if (ok) exit
395 : : end do
396 : :
397 : 306801 : end subroutine tdep_SearchS_2at
398 : :
399 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
400 : 36986 : subroutine tdep_SearchS_3at(Invar,iatom,jatom,katom,eatom,fatom,gatom,Isym3at,Sym,xred_ideal)
401 : :
402 : : integer, intent(in) :: iatom,jatom,katom,eatom,fatom,gatom
403 : : type(atdep_dataset_type),intent(in) :: Invar
404 : : type(Symmetries_type),intent(in) :: Sym
405 : : integer, intent(inout) :: Isym3at(2)
406 : : double precision, intent(in) :: xred_ideal(3,Invar%natom)
407 : :
408 : : integer :: isym,ee,ff,gg,ii
409 : : integer :: iatom_unitcell,jatom_unitcell,katom_unitcell
410 : : integer :: vecti(3),vectj(3),vectk(3),indsym2(8)
411 : : integer :: lattef(3),lattfg(3),lattge(3),lattfe(3),lattgf(3),latteg(3)
412 : 73972 : double precision :: temp(3),tmp_store(3,Invar%natom_unitcell)
413 : 73972 : double precision :: tmp(3,Invar%natom)
414 : : logical :: ok
415 : :
416 : : ! Search the couple of atoms (e,f) obtained through the transformation (R+t)
417 : : ! starting from (i,j)
418 : : ! In that case, note that two cases are possible:
419 : : ! - The bond is just transformed, so indsym(4,e)=i and indsym(4,f)=j
420 : : ! - The bond is transformed and reversed, so indsym(4,e)=j and indsym(4,f)=i
421 : :
422 : : ! Store the positions of the atoms in the motif
423 [ + + ]: 124446 : do ii=1,Invar%natom_unitcell
424 : 124446 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,xred_ideal(:,ii),1,0.d0,tmp_store(:,ii),1)
425 : : end do
426 : :
427 : : ! Search the atom equivalent to iatom in the (reference) unitcell
428 [ + + ]: 6147570 : do ii=1,Invar%natom
429 : 6147570 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,xred_ideal(:,ii),1,0.d0,tmp(:,ii),1)
430 : : end do
431 : : ! Note that in the (present) particular case : iatom_unitcell=iatom and vecti(:)=zero
432 : 36986 : iatom_unitcell=mod(iatom-1,Invar%natom_unitcell)+1
433 [ + + ]: 147944 : vecti(:)=nint(tmp(:,iatom)-tmp(:,iatom_unitcell))
434 : :
435 : : ! Search the atom equivalent to jatom in the (reference) unitcell.
436 : : ! jatom can be outside the box (according to the value of iatom)
437 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
438 [ + + ]: 147944 : temp(:)=xred_ideal(:,jatom)-xred_ideal(:,iatom)
439 : 36986 : call tdep_make_inbox(temp,1,1d-3,temp)
440 [ + + ]: 147944 : temp(:)=xred_ideal(:,iatom)+temp(:)
441 : 36986 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,temp(:),1,0.d0,tmp(:,jatom),1)
442 : 36986 : jatom_unitcell=mod(jatom-1,Invar%natom_unitcell)+1
443 [ + + ]: 147944 : vectj(:)=nint(tmp(:,jatom)-tmp_store(:,jatom_unitcell))
444 : :
445 : : ! Search the atom equivalent to katom in the (reference) unitcell.
446 : : ! katom can be outside the box (according to the value of iatom)
447 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
448 [ + + ]: 147944 : temp(:)=xred_ideal(:,katom)-xred_ideal(:,iatom)
449 : 36986 : call tdep_make_inbox(temp,1,1d-3,temp)
450 [ + + ]: 147944 : temp(:)=xred_ideal(:,iatom)+temp(:)
451 : 36986 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,temp(:),1,0.d0,tmp(:,katom),1)
452 : 36986 : katom_unitcell=mod(katom-1,Invar%natom_unitcell)+1
453 [ + + ]: 147944 : vectk(:)=nint(tmp(:,katom)-tmp_store(:,katom_unitcell))
454 : :
455 : : ! To understand the meaning of "latt", see SearchS_1at
456 : 36986 : ok=.false.
457 [ + + ]: 907139 : do isym=1,Sym%nsym
458 : : ! TODO : A CHECKER!!!!!!!!!!!!!!!!!
459 : 892539 : indsym2(:)=zero
460 : 892539 : call tdep_calc_indsym2(Invar,eatom,fatom,indsym2,isym,Sym,xred_ideal)
461 : 892539 : ee=indsym2(4)
462 [ + + ]: 3570156 : lattef(:)=indsym2(1:3)
463 [ + + ]: 3570156 : lattfe(:)=indsym2(5:7)
464 : :
465 : 892539 : indsym2(:)=zero
466 : 892539 : call tdep_calc_indsym2(Invar,fatom,gatom,indsym2,isym,Sym,xred_ideal)
467 : 892539 : ff=indsym2(4)
468 [ + + ]: 3570156 : lattfg(:)=indsym2(1:3)
469 [ + + ]: 3570156 : lattgf(:)=indsym2(5:7)
470 : :
471 : 892539 : indsym2(:)=zero
472 : 892539 : call tdep_calc_indsym2(Invar,gatom,eatom,indsym2,isym,Sym,xred_ideal)
473 : 892539 : gg=indsym2(4)
474 [ + + ]: 3570156 : lattge(:)=indsym2(1:3)
475 [ + + ]: 3570156 : latteg(:)=indsym2(5:7)
476 [ + + + + ]: 892539 : if (ee==iatom_unitcell.and.ff==jatom_unitcell.and.gg==katom_unitcell) then
477 : : !FB if (Invar%debug) then
478 : : !FB write(Invar%stdout,*) 'For ee=iatom, ff=jatom and gg=katom, isym=',isym
479 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vecti(1),' - ',vectj(1),' - ',lattef(1),' + ',lattfe(1)
480 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vecti(2),' - ',vectj(2),' - ',lattef(2),' + ',lattfe(2)
481 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vecti(3),' - ',vectj(3),' - ',lattef(3),' + ',lattfe(3)
482 : : !FB write(Invar%stdout,*) ' '
483 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectj(1),' - ',vectk(1),' - ',lattfg(1),' + ',lattgf(1)
484 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectj(2),' - ',vectk(2),' - ',lattfg(2),' + ',lattgf(2)
485 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectj(3),' - ',vectk(3),' - ',lattfg(3),' + ',lattgf(3)
486 : : !FB write(Invar%stdout,*) ' '
487 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectk(1),' - ',vecti(1),' - ',lattge(1),' + ',latteg(1)
488 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectk(2),' - ',vecti(2),' - ',lattge(2),' + ',latteg(2)
489 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectk(3),' - ',vecti(3),' - ',lattge(3),' + ',latteg(3)
490 : : !FB write(Invar%stdout,*) ' '
491 : : !FB end if
492 : : if ((sum(abs((vecti(:)-lattef(:))-(vectj(:)-lattfe(:)))).lt.tol8).and.&
493 [ + + + + : 1812760 : & (sum(abs((vectj(:)-lattfg(:))-(vectk(:)-lattgf(:)))).lt.tol8).and.&
+ + + + +
- ]
494 : : & (sum(abs((vectk(:)-lattge(:))-(vecti(:)-latteg(:)))).lt.tol8)) then
495 : 22386 : Isym3at(1)=isym
496 : 22386 : Isym3at(2)=1
497 : : ok=.true.
498 : : end if
499 : : end if
500 : 14600 : if (ok) exit
501 : : end do
502 : :
503 : 36986 : end subroutine tdep_SearchS_3at
504 : :
505 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
506 : 66824 : subroutine tdep_SearchS_4at(Invar,iatom,jatom,katom,latom,eatom,fatom,gatom,hatom,Isym4at,Sym,xred_ideal)
507 : :
508 : : integer, intent(in) :: iatom,jatom,katom,latom,eatom,fatom,gatom,hatom
509 : : type(atdep_dataset_type),intent(in) :: Invar
510 : : type(Symmetries_type),intent(in) :: Sym
511 : : integer, intent(inout) :: Isym4at(2)
512 : : double precision, intent(in) :: xred_ideal(3,Invar%natom)
513 : :
514 : : integer :: isym,ee,ff,gg,hh,ii
515 : : integer :: iatom_unitcell,jatom_unitcell,katom_unitcell,latom_unitcell
516 : : integer :: vecti(3),vectj(3),vectk(3),vectl(3),indsym2(8)
517 : : integer :: lattef(3),lattfe(3)
518 : : integer :: latteg(3),lattge(3)
519 : : integer :: latteh(3),latthe(3)
520 : : integer :: lattfg(3),lattgf(3)
521 : : integer :: lattfh(3),latthf(3)
522 : : integer :: lattgh(3),latthg(3)
523 : 133648 : double precision :: temp(3),tmp_store(3,Invar%natom_unitcell)
524 : 133648 : double precision :: tmp(3,Invar%natom)
525 : : logical :: ok
526 : :
527 : : ! Search the couple of atoms (e,f) obtained through the transformation (R+t)
528 : : ! starting from (i,j)
529 : : ! In that case, note that two cases are possible:
530 : : ! - The bond is just transformed, so indsym(4,e)=i and indsym(4,f)=j
531 : : ! - The bond is transformed and reversed, so indsym(4,e)=j and indsym(4,f)=i
532 : :
533 : : ! Store the positions of the atoms in the motif
534 [ + + ]: 226052 : do ii=1,Invar%natom_unitcell
535 : 226052 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,xred_ideal(:,ii),1,0.d0,tmp_store(:,ii),1)
536 : : end do
537 : :
538 : : ! Search the atom equivalent to iatom in the (reference) unitcell
539 [ + + ]: 10948880 : do ii=1,Invar%natom
540 : 10948880 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,xred_ideal(:,ii),1,0.d0,tmp(:,ii),1)
541 : : end do
542 : : ! Note that in the (present) particular case : iatom_unitcell=iatom and vecti(:)=zero
543 : 66824 : iatom_unitcell=mod(iatom-1,Invar%natom_unitcell)+1
544 [ + + ]: 267296 : vecti(:)=nint(tmp(:,iatom)-tmp(:,iatom_unitcell))
545 : :
546 : : ! Search the atom equivalent to jatom in the (reference) unitcell.
547 : : ! jatom can be outside the box (according to the value of iatom)
548 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
549 [ + + ]: 267296 : temp(:)=xred_ideal(:,jatom)-xred_ideal(:,iatom)
550 : 66824 : call tdep_make_inbox(temp,1,1d-3,temp)
551 [ + + ]: 267296 : temp(:)=xred_ideal(:,iatom)+temp(:)
552 : 66824 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,temp(:),1,0.d0,tmp(:,jatom),1)
553 : 66824 : jatom_unitcell=mod(jatom-1,Invar%natom_unitcell)+1
554 [ + + ]: 267296 : vectj(:)=nint(tmp(:,jatom)-tmp_store(:,jatom_unitcell))
555 : :
556 : : ! Search the atom equivalent to katom in the (reference) unitcell.
557 : : ! katom can be outside the box (according to the value of iatom)
558 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
559 [ + + ]: 267296 : temp(:)=xred_ideal(:,katom)-xred_ideal(:,iatom)
560 : 66824 : call tdep_make_inbox(temp,1,1d-3,temp)
561 [ + + ]: 267296 : temp(:)=xred_ideal(:,iatom)+temp(:)
562 : 66824 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,temp(:),1,0.d0,tmp(:,katom),1)
563 : 66824 : katom_unitcell=mod(katom-1,Invar%natom_unitcell)+1
564 [ + + ]: 267296 : vectk(:)=nint(tmp(:,katom)-tmp_store(:,katom_unitcell))
565 : :
566 : : ! Search the atom equivalent to latom in the (reference) unitcell.
567 : : ! latom can be outside the box (according to the value of iatom)
568 : : ! so we use the inbox procedure to put the distance within [-0.5,0.5[
569 [ + + ]: 267296 : temp(:)=xred_ideal(:,latom)-xred_ideal(:,iatom)
570 : 66824 : call tdep_make_inbox(temp,1,1d-3,temp)
571 [ + + ]: 267296 : temp(:)=xred_ideal(:,iatom)+temp(:)
572 : 66824 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,temp(:),1,0.d0,tmp(:,latom),1)
573 : 66824 : latom_unitcell=mod(latom-1,Invar%natom_unitcell)+1
574 [ + + ]: 267296 : vectl(:)=nint(tmp(:,latom)-tmp_store(:,latom_unitcell))
575 : :
576 : : ! To understand the meaning of "latt", see SearchS_1at
577 : 66824 : ok=.false.
578 [ + + ]: 2069687 : do isym=1,Sym%nsym
579 : : ! TODO : A CHECKER!!!!!!!!!!!!!!!!!
580 : 2033675 : indsym2(:)=zero
581 : 2033675 : call tdep_calc_indsym2(Invar,eatom,fatom,indsym2,isym,Sym,xred_ideal)
582 : 2033675 : ee=indsym2(4)
583 : 2033675 : ff=indsym2(8)
584 [ + + ]: 8134700 : lattef(:)=indsym2(1:3)
585 [ + + ]: 8134700 : lattfe(:)=indsym2(5:7)
586 : :
587 : 2033675 : indsym2(:)=zero
588 : 2033675 : call tdep_calc_indsym2(Invar,eatom,gatom,indsym2,isym,Sym,xred_ideal)
589 : 2033675 : gg=indsym2(8)
590 [ + + ]: 8134700 : latteg(:)=indsym2(1:3)
591 [ + + ]: 8134700 : lattge(:)=indsym2(5:7)
592 : :
593 : 2033675 : indsym2(:)=zero
594 : 2033675 : call tdep_calc_indsym2(Invar,eatom,hatom,indsym2,isym,Sym,xred_ideal)
595 : 2033675 : hh=indsym2(8)
596 [ + + ]: 8134700 : latteh(:)=indsym2(1:3)
597 [ + + ]: 8134700 : latthe(:)=indsym2(5:7)
598 : :
599 : 2033675 : indsym2(:)=zero
600 : 2033675 : call tdep_calc_indsym2(Invar,fatom,gatom,indsym2,isym,Sym,xred_ideal)
601 [ + + ]: 8134700 : lattfg(:)=indsym2(1:3)
602 [ + + ]: 8134700 : lattgf(:)=indsym2(5:7)
603 : :
604 : 2033675 : indsym2(:)=zero
605 : 2033675 : call tdep_calc_indsym2(Invar,fatom,hatom,indsym2,isym,Sym,xred_ideal)
606 [ + + ]: 8134700 : lattfh(:)=indsym2(1:3)
607 [ + + ]: 8134700 : latthf(:)=indsym2(5:7)
608 : :
609 : 2033675 : indsym2(:)=zero
610 : 2033675 : call tdep_calc_indsym2(Invar,gatom,hatom,indsym2,isym,Sym,xred_ideal)
611 [ + + ]: 8134700 : lattgh(:)=indsym2(1:3)
612 [ + + ]: 8134700 : latthg(:)=indsym2(5:7)
613 : :
614 [ + + + + ]: 2033675 : if (ee==iatom_unitcell.and.ff==jatom_unitcell.and.gg==katom_unitcell.and.hh==latom_unitcell) then
615 : : !FB if (Invar%debug) then
616 : : !FB write(Invar%stdout,*) 'For indef=iatom, indfg=jatom and indge=katom, isym=',isym
617 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vecti(1),' - ',vectj(1),' - ',lattef(1),' + ',lattfe(1)
618 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vecti(2),' - ',vectj(2),' - ',lattef(2),' + ',lattfe(2)
619 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vecti(3),' - ',vectj(3),' - ',lattef(3),' + ',lattfe(3)
620 : : !FB write(Invar%stdout,*) ' '
621 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectj(1),' - ',vectk(1),' - ',lattfg(1),' + ',lattgf(1)
622 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectj(2),' - ',vectk(2),' - ',lattfg(2),' + ',lattgf(2)
623 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectj(3),' - ',vectk(3),' - ',lattfg(3),' + ',lattgf(3)
624 : : !FB write(Invar%stdout,*) ' '
625 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectk(1),' - ',vecti(1),' - ',lattge(1),' + ',latteg(1)
626 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectk(2),' - ',vecti(2),' - ',lattge(2),' + ',latteg(2)
627 : : !FB write(Invar%stdout,'(4(a,i4))') ' + ',vectk(3),' - ',vecti(3),' - ',lattge(3),' + ',latteg(3)
628 : : !FB write(Invar%stdout,*) ' '
629 : : !FB end if
630 : : if ((sum(abs((vecti(:)-lattef(:))-(vectj(:)-lattfe(:)))).lt.tol8).and.&
631 : : & (sum(abs((vecti(:)-latteg(:))-(vectk(:)-lattge(:)))).lt.tol8).and.&
632 : : & (sum(abs((vecti(:)-latteh(:))-(vectl(:)-latthe(:)))).lt.tol8).and.&
633 : : & (sum(abs((vectj(:)-lattfg(:))-(vectk(:)-lattgf(:)))).lt.tol8).and.&
634 [ + + + + : 5712882 : & (sum(abs((vectj(:)-lattfh(:))-(vectl(:)-latthf(:)))).lt.tol8).and.&
+ + + + +
+ + + + +
+ + + - ]
635 : : & (sum(abs((vectk(:)-lattgh(:))-(vectl(:)-latthg(:)))).lt.tol8)) then
636 : 30812 : Isym4at(1)=isym
637 : 30812 : Isym4at(2)=1
638 : : ok=.true.
639 : : end if
640 : : end if
641 : 36012 : if (ok) exit
642 : : end do
643 : :
644 : 66824 : end subroutine tdep_SearchS_4at
645 : :
646 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
647 : 20058508 : subroutine tdep_calc_indsym2(Invar,iatom,jatom,indsym2,isym,Sym,xred_ideal)
648 : :
649 : : integer, intent(in) :: iatom,jatom
650 : : type(atdep_dataset_type),intent(in) :: Invar
651 : : type(Symmetries_type),intent(in) :: Sym
652 : : integer, intent(out) :: indsym2(8)
653 : : double precision, intent(in) :: xred_ideal(3,Invar%natom)
654 : :
655 : : integer :: isym,mu,ii
656 : : integer :: iatom_unitcell,jatom_unitcell
657 : 40117016 : integer :: vecti(3),vectj(3),vectsym(4,Sym%nsym)
658 : 40117016 : double precision :: tmpi(3,Invar%natom),tmpj(3,Invar%natom),temp3(3,1),tmp_store(3,Invar%natom_unitcell)
659 : :
660 : : ! Store the positions of the atoms in the motif
661 [ + + ]: 66414480 : do ii=1,Invar%natom_unitcell
662 : 66414480 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,xred_ideal(:,ii),1,0.d0,tmp_store(:,ii),1)
663 : : end do
664 : :
665 : : ! Search the matrix transformation going from (k,l) to (i,j)
666 [ + + + + ]: 14104154820 : tmpi(:,:)=0.d0
667 [ + + + + ]: 14104154820 : tmpj(:,:)=0.d0
668 : : ! For a single iatom
669 : 20058508 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,xred_ideal(:,iatom),1,0.d0,tmpi(:,iatom),1)
670 : 20058508 : iatom_unitcell=mod(iatom-1,Invar%natom_unitcell)+1
671 [ + + ]: 80234032 : vecti(:)=nint(tmpi(:,iatom)-tmp_store(:,iatom_unitcell))
672 [ + + + + ]: 4720379728 : vectsym(:,:)=0
673 [ + + ]: 80234032 : do mu=1,3 ! Apply inverse transformation to original coordinates. Note transpose of symrec.
674 : 80234032 : vectsym(mu,isym) = Sym%symrec(1,mu,isym)*vecti(1)+Sym%symrec(2,mu,isym)*vecti(2)+Sym%symrec(3,mu,isym)*vecti(3)
675 : : end do
676 [ + + ]: 100292540 : indsym2(1:4)=Sym%indsym(1:4,isym,iatom_unitcell)+vectsym(1:4,isym)
677 : :
678 : : ! For a couple of (iatom,jatom). The (iatom,jatom) vector depends on the position of iatom (due to PBC)
679 : : !FB if (Invar%debug) write(Invar%stdout,*) '=========================================='
680 : : !FB if (Invar%debug) write(Invar%stdout,'(a,i4,a,3(f10.5,x))') 'For iatom=',iatom,' with xred=',xred_ideal(:,iatom)
681 : : !FB if (Invar%debug) write(Invar%stdout,*) ' =========================================='
682 : : !FB if (Invar%debug) write(Invar%stdout,'(a,i4,a,3(f10.5,x))') ' For jatom=',jatom,' with xred=',xred_ideal(:,jatom)
683 [ + + ]: 80234032 : temp3(:,1)=xred_ideal(:,jatom)-xred_ideal(:,iatom)
684 : 20058508 : call tdep_make_inbox(temp3(:,1),1,1d-3,temp3(:,1))
685 [ + + ]: 80234032 : temp3(:,1)=xred_ideal(:,iatom)+temp3(:,1)
686 : 20058508 : call DGEMV('T',3,3,1.d0,Invar%multiplicity(:,:),3,temp3(:,1),1,0.d0,tmpj(:,jatom),1)
687 : 20058508 : jatom_unitcell=mod(jatom-1,Invar%natom_unitcell)+1
688 [ + + ]: 80234032 : vectj(:)=nint(tmpj(:,jatom)-tmp_store(:,jatom_unitcell))
689 [ + + + + ]: 4720379728 : vectsym(:,:)=0
690 [ + + ]: 80234032 : do mu=1,3 ! Apply inverse transformation to original coordinates. Note transpose of symrec.
691 : 80234032 : vectsym(mu,isym) = Sym%symrec(1,mu,isym)*vectj(1)+Sym%symrec(2,mu,isym)*vectj(2)+Sym%symrec(3,mu,isym)*vectj(3)
692 : : end do
693 [ + + ]: 100292540 : indsym2(5:8)=Sym%indsym(1:4,isym,jatom_unitcell)+vectsym(1:4,isym)
694 : : !FB if (Invar%debug) then
695 : : !FB write(Invar%stdout,'(a,i2,a,i4,a,i4,a,3(i4,x),a,i2,a,3(i4,x),a,i2,a)') ' indsym2(isym=',isym,',',iatom,',',jatom,')=',indsym2(1:3),&
696 : : !FB& '|iat=',indsym2(4),'|',indsym2(5:7),'|iat=',indsym2(8),'|'
697 : : !FB end if
698 : :
699 : 20058508 : end subroutine tdep_calc_indsym2
700 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
701 : :
702 : 44 : subroutine tdep_destroy_sym(Sym)
703 : :
704 : : type(Symmetries_type),intent(inout) :: Sym
705 : :
706 [ - + ]: 44 : ABI_FREE(Sym%ptsymrel)
707 [ - + ]: 44 : ABI_FREE(Sym%S_ref)
708 [ - + ]: 44 : ABI_FREE(Sym%S_inv)
709 [ - + ]: 44 : ABI_FREE(Sym%xred_zero)
710 [ - + ]: 44 : ABI_FREE(Sym%tnons)
711 [ - + ]: 44 : ABI_FREE(Sym%symafm)
712 [ - + ]: 44 : ABI_FREE(Sym%symrec)
713 [ - + ]: 44 : ABI_FREE(Sym%symrel)
714 [ - + ]: 44 : ABI_FREE(Sym%indsym)
715 : :
716 : 44 : end subroutine tdep_destroy_sym
717 : :
718 : : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
719 [ # # # # : 0 : end module m_tdep_sym
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
|
