Line data Source code
1 : MODULE m_mapatom
2 : use m_juDFT
3 : !*******************************************************************
4 : ! determines the group operation which maps the representive
5 : ! atom into its equivalent atoms c.l.fu
6 : !*******************************************************************
7 : CONTAINS
8 80 : SUBROUTINE mapatom(sym,atoms,cell,input,noco,gfinp)
9 : !
10 : ! if (l_f) setup multab,invtab,invarop,invarind for force_a12 & 21
11 : !***********************************************************************
12 : ! the contribution to the hamiltonian and to the overlap matrix in a
13 : ! system with inversion symmetry from one muffin tin is the complex
14 : ! conjugate of the contribution from the "invers" muffin tin. this fact
15 : ! can be exploited to save cpu-time in hssphn. Thus, it is nessessary to
16 : ! know whether an atom can be mapped onto an equivalent atom via 3d
17 : ! inversion. where both atoms have to belong to the same unit cell, i.e.
18 : ! the are not related to each other by a lattice translation. therefore,
19 : ! an array invsatom is set up.
20 : ! invsatom(natom) =
21 : ! 0 if the atom cannot be mapped onto an eqivalent atom via inversion
22 : ! 1 if the atom can be mapped onto an eqivalent atom via inversion, and
23 : ! has a smaller atom index than the related atom
24 : ! 2 if the atom can be mapped onto an eqivalent atom via inversion, and
25 : ! has a bigger atom index than the related atom
26 : ! p.kurz aug. 1996
27 : !***********************************************************************
28 :
29 : USE m_types_sym
30 : use m_types_atoms
31 : use m_types_cell
32 : use m_types_input
33 : use m_types_noco
34 : use m_types_gfinp
35 : USE m_constants
36 : USE m_socsym
37 :
38 : IMPLICIT NONE
39 :
40 : TYPE(t_sym),INTENT(INOUT):: sym
41 : TYPE(t_atoms),INTENT(IN) :: atoms
42 : TYPE(t_cell),INTENT(IN) :: cell
43 : TYPE(t_input),INTENT(IN) :: input
44 : TYPE(t_noco),INTENT(IN) :: noco
45 : TYPE(t_gfinp),INTENT(IN) :: gfinp
46 :
47 : ! .. Local Scalars ..
48 : REAL s3,norm
49 : INTEGER i,icount,j,j1,j2,j3,jop,n,na,nat,natp,iop,nat1,nat2,nb,na_r
50 : INTEGER k,ij,n1,n2,ix,iy,iz,na2
51 : REAL, PARAMETER :: del = 1.0e-4
52 :
53 : ! .. Local Arrays ..
54 : INTEGER mt(3,3),mp(3,3)
55 : REAL aamat(3,3),sum_tau_lat(3),sum_taual(3)
56 : REAL gam(3),gaminv(3),gamr(3),sr(3),ttau(3),gammap(3)
57 80 : LOGICAL error(sym%nop)
58 :
59 : ! CALL dotset(&
60 : ! & cell,&
61 : ! & aamat, )
62 3200 : aamat=matmul(transpose(cell%amat),cell%amat)
63 :
64 80 : IF (ALLOCATED(sym%ngopr)) deallocate(sym%ngopr)
65 80 : if (allocated(sym%invsatnr)) deallocate(sym%invsatnr)
66 80 : if (allocated(sym%invarop)) deallocate(sym%invarop)
67 80 : if (allocated(sym%invarind)) deallocate(sym%invarind)
68 80 : if (allocated(sym%invsat)) deallocate(sym%invsat)
69 :
70 240 : ALLOCATE(sym%invsatnr(atoms%nat))
71 320 : ALLOCATE(sym%invarop(atoms%nat,sym%nop))
72 160 : ALLOCATE(sym%invarind(atoms%nat))
73 160 : ALLOCATE(sym%ngopr(atoms%nat))
74 160 : ALLOCATE(sym%invsat(atoms%nat))
75 :
76 :
77 80 : IF (noco%l_soc) THEN ! check once more here...
78 : CALL soc_sym(&
79 : & sym%nop,sym%mrot,noco%theta_inp,noco%phi_inp,cell%amat,&
80 16 : & error)
81 : ELSE
82 1213 : error(:) = .false.
83 : ENDIF
84 :
85 80 : WRITE (oUnit,FMT=8000)
86 : 8000 FORMAT (/,/,5x,'group operations on equivalent atoms:')
87 217 : DO n = 1,atoms%ntype
88 137 : nat1 = atoms%firstAtom(n)
89 137 : nat2 = nat1 + atoms%neq(n) - 1
90 137 : sym%ngopr(nat1) = 1
91 : !+gu
92 137 : na_r = nat1
93 398 : DO na = nat1,nat2
94 181 : IF (sym%ntypsy(na).NE.sym%ntypsy(na_r)) na_r = na
95 : !-gu
96 724 : DO i = 1,3
97 724 : gam(i) = atoms%taual(i,na)
98 : END DO
99 181 : sym%invarind(na) = 0
100 181 : icount = 0
101 2507 : DO jop = 1,sym%nop
102 9304 : DO i = 1,3
103 6978 : gamr(i) = 0.
104 27912 : DO j = 1,3
105 27912 : gamr(i) = gamr(i) + sym%mrot(i,j,jop)*gam(j)
106 : END DO
107 9304 : gamr(i) = gamr(i) + sym%tau(i,jop)
108 : END DO
109 9304 : DO i = 1,3
110 6978 : gaminv(i) = gamr(i) - atoms%taual(i,na)
111 9304 : gamr(i) = gamr(i) - atoms%taual(i,nat1) ! cf local_sym
112 : END DO
113 2326 : IF (icount.EQ.0) THEN
114 2010 : DO j3 = -2,2
115 1675 : sr(3) = gamr(3) + real(j3)
116 10385 : DO j2 = -2,2
117 8375 : sr(2) = gamr(2) + real(j2)
118 51925 : DO j1 = -2,2
119 41875 : sr(1) = gamr(1) + real(j1)
120 670000 : s3 = sqrt(dot_product(matmul(sr,aamat),sr))
121 50250 : IF ((s3.LT.del).AND.(.not.error(jop))) THEN
122 181 : icount = icount + 1
123 181 : sym%ngopr(na) = jop
124 : END IF
125 : END DO
126 : END DO
127 : END DO
128 : END IF
129 : !
130 : ! search for operations which leave taual invariant
131 : !
132 2507 : IF (input%l_f.OR.(atoms%n_denmat+gfinp%n.GT.0)) THEN
133 4524 : DO j3 = -2,2
134 3770 : sr(3) = gaminv(3) + real(j3)
135 24946 : DO j2 = -2,2
136 18850 : sr(2) = gaminv(2) + real(j2)
137 116870 : DO j1 = -2,2
138 94250 : sr(1) = gaminv(1) + real(j1)
139 1508000 : s3 = sqrt(dot_product(matmul(sr,aamat),sr))
140 113100 : IF (s3.LT.del) THEN
141 678 : sym%invarind(na) = sym%invarind(na) + 1
142 678 : sym%invarop(na,sym%invarind(na)) = jop
143 : END IF
144 : END DO
145 : END DO
146 : END DO
147 : ENDIF
148 : !
149 : ! end of operations
150 : ENDDO
151 181 : IF (icount.LE.0) THEN
152 0 : write(oUnit,*) "Mapping failed for atom:",nat1
153 0 : write(oUnit,*) "No of symmetries tested:",sym%nop
154 0 : CALL juDFT_error("mapatom",calledby="mapatom")
155 : ENDIF
156 318 : WRITE (oUnit,FMT=8010) nat1,na,sym%ngopr(na)
157 : 8010 FORMAT (5x,'atom',i5,' can be mapped into atom',i5,&
158 : & ' through group operation',i4)
159 : !
160 : ! end of equivalent atoms
161 : ENDDO
162 : ! end of different types of atoms
163 : ENDDO
164 :
165 : !------------------------- FORCE PART -------------------------------
166 : !+gu this is the remainder of spgset necessary for force calculations
167 : !
168 80 : IF (input%l_f.OR.(atoms%n_denmat+gfinp%n.GT.0)) THEN
169 :
170 : WRITE (oUnit,FMT=&
171 26 : & '(//,"list of operations which leave taual invariant",/)')
172 88 : DO na = 1,nat2
173 62 : WRITE (oUnit,FMT='("atom nr.",i3,3x,(t14,"ops are:",24i3))') na,&
174 882 : & (sym%invarop(na,nb),nb=1,sym%invarind(na))
175 : END DO
176 :
177 : ENDIF
178 :
179 : !IF(gfinp%n > 0) THEN ! TODO: Bind this to a switch gfinp%n > 0 .OR. juPhon%l_dos
180 : ! Create mapped_atom array
181 : ! Store, where atoms are mapped to when symmetry operations are applied
182 2827 : IF ( .NOT. ALLOCATED(sym%mapped_atom)) ALLOCATE(sym%mapped_atom(sym%nop,atoms%nat),source=0)
183 :
184 217 : DO n = 1 , atoms%ntype
185 398 : DO nat = atoms%firstAtom(n), atoms%firstAtom(n) + atoms%neq(n) - 1
186 :
187 2644 : DO iop = 1, sym%nop
188 58150 : gamr = matmul(sym%mrot(:,:,iop), atoms%taual(:,nat))
189 9304 : gamr = gamr + sym%tau(:,iop)
190 :
191 2326 : icount = 0
192 6057 : DO natp = atoms%firstAtom(n), atoms%firstAtom(n) + atoms%neq(n) - 1
193 14200 : gammap = gamr - atoms%taual(:,natp)
194 5876 : IF (icount.EQ.0) THEN
195 17652 : DO j3 = -2,2
196 14710 : sr(3) = gammap(3) + real(j3)
197 91202 : DO j2 = -2,2
198 73550 : sr(2) = gammap(2) + real(j2)
199 456010 : DO j1 = -2,2
200 367750 : sr(1) = gammap(1) + real(j1)
201 5884000 : s3 = sqrt(dot_product(matmul(sr,aamat),sr))
202 441300 : IF ((s3.LT.del).AND.(.not.error(iop))) THEN
203 2290 : icount = icount + 1
204 2290 : sym%mapped_atom(iop,nat) = natp
205 : END IF
206 : END DO
207 : END DO
208 : END DO
209 : END IF
210 : ENDDO
211 : ENDDO
212 : ENDDO
213 : ENDDO
214 :
215 80 : WRITE (oUnit,FMT='(//,"list of atoms, which are mapped to by symmetry operations",/)')
216 261 : DO na = 1,atoms%nat
217 2587 : WRITE (oUnit,FMT='("atom nr.",i3,3x,(t14,"mapped atoms are:",24i3))') na,(sym%mapped_atom(nb,na),nb=1,sym%nop)
218 : END DO
219 : !ENDIF
220 :
221 :
222 : !------------------------- FORCE PART ENDS --------------------------
223 : !
224 : ! check closure ; note that: {R|t} tau = R^{-1} tau - R^{-1} t
225 : !
226 : !---> loop over all operations
227 : !
228 80 : WRITE (oUnit,FMT=8040)
229 : 8040 FORMAT (/,/,' multiplication table',/,/)
230 45426 : sym%multab = 0
231 1286 : DO j=1,sym%nop
232 :
233 : !---> multiply {R_j|t_j}{R_i|t_i}
234 45426 : DO i=1,sym%nop
235 2957380 : mp = matmul( sym%mrot(:,:,j) , sym%mrot(:,:,i) )
236 1235920 : ttau = sym%tau(:,j) + matmul( sym%mrot(:,:,j) , sym%tau(:,i) )
237 176560 : ttau = ttau - anint( ttau - 1.e-7 )
238 :
239 : !---> determine which operation this is
240 2001244 : DO k=1,sym%nop
241 9054396 : IF( all( mp(:,:) == sym%mrot(:,:,k) ) .AND.&
242 44140 : & ALL( abs( ttau(:)-sym%tau(:,k) ) < 1.e-7 ) ) THEN
243 44140 : IF (sym%multab(j,i) .EQ. 0 ) THEN
244 44140 : sym%multab(j,i) = k
245 44140 : IF (k .EQ. 1) sym%invtab(j)=i
246 : ELSE
247 0 : WRITE(oUnit,'(" Symmetry error: multiple ops")')
248 0 : CALL juDFT_error("Multiple ops",calledby="mapatom")
249 : ENDIF
250 : ENDIF
251 : ENDDO
252 :
253 45346 : IF (sym%multab(j,i).EQ.0) THEN
254 0 : WRITE (oUnit,'(" Group not closed")')
255 0 : WRITE (oUnit,'(" j , i =",2i4)') j,i
256 : CALL juDFT_error("mapatom: group not closed",calledby&
257 0 : & ="mapatom")
258 : ENDIF
259 : ENDDO
260 : ENDDO
261 :
262 1286 : DO n1 = 1,sym%nop
263 45426 : WRITE (oUnit,FMT=8060) (sym%multab(n1,n2),n2=1,sym%nop)
264 : END DO
265 : 8060 FORMAT (1x,48i3)
266 80 : WRITE (oUnit,FMT='(//," inverse operations",//)')
267 1286 : DO n1 = 1,sym%nop
268 1286 : WRITE (oUnit,FMT=8060) n1,sym%invtab(n1)
269 : END DO
270 :
271 261 : DO na = 1,atoms%nat
272 181 : sym%invsat(na) = 0
273 261 : sym%invsatnr(na) = 0
274 : END DO
275 :
276 80 : IF (.not.(noco%l_soc.and.atoms%n_u+atoms%n_opc>0) .and. atoms%n_hia==0) THEN
277 76 : IF (sym%invs) THEN
278 42 : WRITE (oUnit,FMT=*)
279 119 : DO n = 1,atoms%ntype
280 77 : nat1 = atoms%firstAtom(n)
281 77 : nat2 = nat1 + atoms%neq(n) - 1
282 156 : DO na = nat1,nat2 - 1
283 114 : IF (sym%invsat(na).EQ.0.AND..NOT.noco%l_noco) THEN
284 80 : naloop:DO na2 = na + 1,nat2
285 196 : DO i = 1,3
286 196 : sum_taual(i) = atoms%taual(i,na) + atoms%taual(i,na2)
287 : END DO
288 229 : DO ix = -2,2
289 180 : sum_tau_lat(1) = sum_taual(1) + real(ix)
290 1002 : DO iy = -2,2
291 835 : sum_tau_lat(2) = sum_taual(2) + real(iy)
292 5066 : DO iz = -2,2
293 4113 : sum_tau_lat(3) = sum_taual(3) + real(iz)
294 65808 : norm = sqrt(dot_product(matmul(sum_tau_lat,aamat),sum_tau_lat))
295 4917 : IF (norm.LT.del) THEN
296 31 : sym%invsat(na) = 1
297 31 : sym%invsat(na2) = 2
298 31 : sym%invsatnr(na) = na2
299 31 : sym%invsatnr(na2) = na
300 31 : WRITE (oUnit,FMT=9000) n,na,na2
301 31 : cycle naloop
302 : END IF
303 : END DO
304 : END DO
305 : END DO
306 : END DO naloop
307 : END IF
308 : END DO
309 : END DO
310 42 : WRITE (oUnit,FMT=*) sym%invsat
311 : 9000 FORMAT ('atom type',i3,': atom',i3,' can be mapped into atom',i3,&
312 : & ' via 3d inversion')
313 : END IF
314 : END IF
315 :
316 80 : END SUBROUTINE mapatom
317 46466 : END MODULE m_mapatom
|
