Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_chkmt
8 : USE m_juDFT
9 : !---------------------------------------------------------------------
10 : ! Check muffin tin radii and determine a reasonable choice for MTRs.
11 : ! Derive also other parameters for the input file, to provide some
12 : ! help in the out-file.
13 : ! GM'16
14 : !
15 : ! Note: This procedure is prepared to also return the next nearest
16 : ! neighbors. That data is stored in nearestNeighborAtomsSorted,
17 : ! nearestNeighborShiftsSorted, numNextNearestNeighbors.
18 : !---------------------------------------------------------------------
19 : CONTAINS
20 :
21 82 : SUBROUTINE chkmt(atoms,input,vacuum,cell ,l_test,&
22 0 : l_gga,noel, kmax,dtild,dvac1,lmax1,jri1,rmt1,dx1,&!optional, if l_gga and ... are present suggestions are calculated
23 2 : overlap)!this is optional, if present and l_test the routine returns the overlaps and does not stop
24 :
25 : USE m_types_fleurinput
26 : USE m_constants
27 : USE m_sort
28 : USE m_inv3
29 : USE m_juDFT
30 :
31 : IMPLICIT NONE
32 :
33 : ! .. Scalar Arguments ..
34 : TYPE(t_atoms),INTENT(IN) :: atoms
35 : TYPE(t_input),INTENT(IN) :: input
36 : TYPE(t_vacuum),INTENT(IN):: vacuum
37 : TYPE(t_cell),INTENT(IN) :: cell
38 :
39 : CHARACTER*3, INTENT (IN),OPTIONAL :: noel(atoms%ntype)
40 : LOGICAL, INTENT (IN),OPTIONAL :: l_gga
41 : LOGICAL, INTENT (IN) ::l_test
42 : REAL, INTENT (OUT),OPTIONAL :: kmax,dtild,dvac1
43 :
44 : ! .. Array Arguments ..
45 : INTEGER, INTENT (OUT),OPTIONAL :: lmax1(atoms%ntype),jri1(atoms%ntype)
46 : REAL, INTENT (OUT),OPTIONAL :: rmt1(atoms%ntype),dx1(atoms%ntype)
47 : REAL,OPTIONAL,INTENT(OUT):: overlap(0:atoms%ntype,atoms%ntype)
48 :
49 : ! .. Local Scalars ..
50 : INTEGER na,n
51 : INTEGER i,j,k,l,jri11,lmax11
52 : INTEGER maxCubeAtoms, iAtom, numAtoms, iNeighborAtom, identicalAtoms
53 : INTEGER typeA, typeB
54 : REAL dx11,rkm,sum_r,facA,facB
55 : REAL rmtMax, rmtMin, rmtMaxDefault, rmtDelta
56 : REAL rmtFac, cubeLength, amatAuxDet
57 : REAL maxSqrDist, dist, currentDist
58 : LOGICAL error, outOfBounds
59 :
60 : ! .. Local Arrays ..
61 : REAL t_rmt(0:103), minRmts(0:103)
62 : REAL amatAux(3,3), invAmatAux(3,3)
63 82 : REAL taualAux(3,atoms%nat), posAux(3,atoms%nat)
64 : REAL minPos(3), maxPos(3), pos(3), point(3), realCellPos(3)
65 : REAL offsetPos(3)
66 82 : REAL nearestAtomDists(atoms%ntype)
67 82 : INTEGER nearestAtoms(atoms%ntype)
68 82 : INTEGER sortedDistList(atoms%ntype)
69 : INTEGER minCubeIndex(3), maxCubeIndex(3), cubeIndex(3)
70 : INTEGER minCellIndices(3), maxCellIndices(3)
71 :
72 82 : INTEGER, ALLOCATABLE :: numAtomsInCubes(:,:,:)
73 82 : INTEGER, ALLOCATABLE :: atomRefsInCubes(:,:,:,:)
74 82 : INTEGER, ALLOCATABLE :: atomsInCubes(:,:,:,:)
75 82 : INTEGER, ALLOCATABLE :: atomShiftsInCubes(:,:,:,:,:)
76 82 : INTEGER, ALLOCATABLE :: refCubes(:,:)
77 82 : INTEGER, ALLOCATABLE :: nearestNeighborRefsSorted(:,:)
78 82 : INTEGER, ALLOCATABLE :: nearestNeighborAtomsSorted(:,:)
79 82 : INTEGER, ALLOCATABLE :: nearestNeighborShiftsSorted(:,:,:)
80 82 : INTEGER, ALLOCATABLE :: numNearestNeighbors(:)
81 82 : INTEGER, ALLOCATABLE :: numNextNearestNeighbors(:)
82 82 : INTEGER, ALLOCATABLE :: neighborAtomRefs(:)
83 82 : INTEGER, ALLOCATABLE :: neighborAtoms(:)
84 82 : INTEGER, ALLOCATABLE :: neighborShifts(:,:)
85 82 : INTEGER, ALLOCATABLE :: distIndexList(:)
86 82 : REAL, ALLOCATABLE :: posInCubes(:,:,:,:,:)
87 82 : REAL, ALLOCATABLE :: refPos(:,:)
88 82 : REAL, ALLOCATABLE :: nearestNeighborDists(:,:)
89 82 : REAL, ALLOCATABLE :: sqrDistances(:)
90 :
91 : ! Plan for this routine:
92 : ! 0. Do initializations and set some constants
93 : ! 1. Locally replace unit cell by an auxiliary unit cell with:
94 : ! a) all atoms within the unit cell
95 : ! b) basis vectors obtained by lattice reduction of the original cell.
96 : ! [not in 1st (this) version of routine. Can be implemented with the LLL algorithm when needed.]
97 : ! 2. Get minimal and maximal coordinates within auxiliary unit cell
98 : ! 3. Construct mesh of cubes covering the auxiliary unit cell and a boundary of width 2*rmtMax + rmtDelta
99 : ! 4. Fill mesh of cubes with atoms
100 : ! a) Store atoms in cubes and representative cube for each atom type
101 : ! 5. For each atom in auxiliary unit cell select cube and collect shortest distances to other atoms in neighborhood
102 : ! a) Sort distances and set MT radii for the atoms
103 : ! 6. Correct bad choices and set missing MT radii, vacuum distances, and other parameters
104 : ! 7. Test old MT radii
105 :
106 :
107 : ! 0. Do initializations and set some constants
108 98 : if (present(overlap)) overlap=0.0
109 :
110 82 : rmtMaxDefault = 2.8
111 82 : rmtMax = rmtMaxDefault
112 82 : rmtMin = 1.0
113 82 : IF (l_test) THEN
114 223 : rmtMax = MAX(rmtMax,MAXVAL(atoms%rmt(:)))
115 223 : rmtMin = MIN(rmtMin,MINVAL(atoms%rmt(:)))
116 : END IF
117 82 : rmtDelta = 0.3
118 82 : IF (input%film) THEN
119 : rmtFac = 0.95
120 : ELSE
121 70 : rmtFac = 0.975
122 : ENDIF
123 8610 : t_rmt(0:103) = 2.3 ! default value
124 656 : t_rmt(1) = 1.0 ; t_rmt(5:9) = 1.3 ; t_rmt(16:17) = 1.8
125 82 : cubeLength = 2*rmtMax+rmtDelta
126 82 : maxCubeAtoms = (FLOOR(cubeLength / (0.7*2.0*rmtMin)) + 1)**3
127 82 : error = .FALSE.
128 :
129 : ! 1. For the 1st version the auxiliary unit cell is just a copy of the original unit cell with
130 : ! all atoms within the cell.
131 :
132 328 : DO i = 1, 3
133 1066 : DO j = 1, 3
134 984 : amatAux(i,j) = cell%amat(i,j)
135 : END DO
136 : END DO
137 :
138 268 : DO i = 1, atoms%nat
139 186 : taualAux(1,i) = atoms%taual(1,i) - FLOOR(atoms%taual(1,i))
140 186 : taualAux(2,i) = atoms%taual(2,i) - FLOOR(atoms%taual(2,i))
141 186 : taualAux(3,i) = atoms%taual(3,i) - FLOOR(atoms%taual(3,i))
142 3058 : posAux(:,i) = MATMUL(amatAux,taualAux(:,i))
143 : END DO
144 :
145 : ! 2. Get minimal and maximal coordinates for auxiliary unit cell
146 :
147 82 : minPos = 0.0
148 82 : maxPos = 0.0
149 :
150 246 : DO i = 0, 1
151 574 : DO j = 0, 1
152 1148 : DO k = 0, 1
153 2952 : DO l = 1, 3
154 1968 : pos(l) = i*amatAux(l,1) + j*amatAux(l,2) + k*amatAux(l,3)
155 1968 : IF (pos(l).GT.maxPos(l)) maxPos(l) = pos(l)
156 2624 : IF (pos(l).LT.minPos(l)) minPos(l) = pos(l)
157 : END DO
158 : END DO
159 : END DO
160 : END DO
161 :
162 : ! 3. Construct cube mesh:
163 : ! In each dimension cube i covers the interval from i*cubeLength to (i+1)*cubeLength
164 : ! Each cube may cover up to maxCubeAtoms atoms. This should be set to a save size.
165 :
166 328 : DO i = 1, 3
167 246 : minPos(i) = minPos(i) - cubeLength
168 246 : maxPos(i) = maxPos(i) + cubeLength
169 246 : minCubeIndex(i) = FLOOR(minPos(i)/cubeLength)
170 328 : maxCubeIndex(i) = CEILING(maxPos(i)/cubeLength)
171 : END DO
172 :
173 : ALLOCATE (numAtomsInCubes(minCubeIndex(1):maxCubeIndex(1),&
174 : minCubeIndex(2):maxCubeIndex(2),&
175 410 : minCubeIndex(3):maxCubeIndex(3)))
176 : ALLOCATE (atomRefsInCubes(maxCubeAtoms,minCubeIndex(1):maxCubeIndex(1),&
177 : minCubeIndex(2):maxCubeIndex(2),&
178 492 : minCubeIndex(3):maxCubeIndex(3)))
179 : ALLOCATE (atomsInCubes(maxCubeAtoms,minCubeIndex(1):maxCubeIndex(1),&
180 : minCubeIndex(2):maxCubeIndex(2),&
181 410 : minCubeIndex(3):maxCubeIndex(3)))
182 : ALLOCATE (atomShiftsInCubes(3,maxCubeAtoms,minCubeIndex(1):maxCubeIndex(1),&
183 : minCubeIndex(2):maxCubeIndex(2),&
184 492 : minCubeIndex(3):maxCubeIndex(3)))
185 : ALLOCATE (posInCubes(3,maxCubeAtoms,minCubeIndex(1):maxCubeIndex(1),&
186 : minCubeIndex(2):maxCubeIndex(2),&
187 492 : minCubeIndex(3):maxCubeIndex(3)))
188 410 : ALLOCATE (refCubes(3,atoms%ntype),refPos(3,atoms%ntype))
189 574 : ALLOCATE (nearestNeighborRefsSorted(maxCubeAtoms,atoms%ntype),numNearestNeighbors(atoms%ntype),numNextNearestNeighbors(atoms%ntype))
190 246 : ALLOCATE (nearestNeighborAtomsSorted(maxCubeAtoms,atoms%ntype))
191 328 : ALLOCATE (nearestNeighborShiftsSorted(3,maxCubeAtoms,atoms%ntype))
192 328 : ALLOCATE (nearestNeighborDists(maxCubeAtoms,atoms%ntype))
193 :
194 12833 : numAtomsInCubes = 0
195 :
196 : ! 4. Fill mesh of cubes with atoms
197 : ! First obtain minimal and maximal indices for relevant unit cells
198 :
199 82 : minCellIndices = 0
200 82 : maxCellIndices = 0
201 :
202 82 : CALL inv3(amatAux,invAmatAux,amatAuxDet)
203 :
204 246 : DO i = 0, 1
205 574 : DO j = 0, 1
206 1148 : DO k = 0, 1
207 656 : point(:) = minPos(:)
208 656 : IF(i.EQ.1) point(1) = maxPos(1)
209 656 : IF(j.EQ.1) point(2) = maxPos(2)
210 656 : IF(k.EQ.1) point(3) = maxPos(3)
211 8528 : realCellPos(:) = matmul(invAmatAux,point(:))
212 2952 : DO l = 1, 3
213 1968 : IF(minCellIndices(l).GT.realCellPos(l)) THEN
214 297 : minCellIndices(l) = FLOOR(realCellPos(l))
215 : END IF
216 2624 : IF(maxCellIndices(l).LT.realCellPos(l)) THEN
217 974 : maxCellIndices(l) = FLOOR(realCellPos(l)) ! Is 'FLOOR' enough?
218 : END IF
219 : END DO
220 : END DO
221 : END DO
222 : END DO
223 :
224 : ! Store atoms in cubes and representative cube for each atom type
225 :
226 534 : DO i = minCellIndices(1), maxCellIndices(1)
227 3358 : DO j = minCellIndices(2), maxCellIndices(2)
228 22504 : DO k = minCellIndices(3), maxCellIndices(3)
229 76912 : DO l = 1, 3
230 76912 : offsetPos(l) = i*amatAux(l,1) + j*amatAux(l,2) + k*amatAux(l,3)
231 : END DO
232 : iAtom = 0
233 47067 : DO n = 1, atoms%ntype
234 74137 : DO na = 1, atoms%neq(n)
235 29894 : iAtom = iAtom + 1
236 119576 : pos(:) = posAux(:,iAtom) + offsetPos(:)
237 : outOfBounds = .FALSE.
238 119576 : DO l = 1, 3
239 89682 : cubeIndex(l) = FLOOR(pos(l)/cubeLength)
240 89682 : IF(cubeIndex(l).LT.minCubeIndex(l)) outOfBounds = .TRUE.
241 119576 : IF(cubeIndex(l).GT.maxCubeIndex(l)) outOfBounds = .TRUE.
242 : END DO
243 54909 : IF(.NOT.outOfBounds) THEN
244 : numAtomsInCubes(cubeIndex(1),cubeIndex(2),cubeIndex(3)) = &
245 11701 : numAtomsInCubes(cubeIndex(1),cubeIndex(2),cubeIndex(3)) + 1
246 11701 : numAtoms = numAtomsInCubes(cubeIndex(1),cubeIndex(2),cubeIndex(3))
247 11701 : IF(numAtoms.GT.maxCubeAtoms) THEN
248 0 : STOP 'ERROR: maxCubeAtoms is not large enough in chkmt.'
249 : END IF
250 11701 : atomRefsInCubes(numAtoms,cubeIndex(1),cubeIndex(2),cubeIndex(3)) = n
251 11701 : atomsInCubes(numAtoms,cubeIndex(1),cubeIndex(2),cubeIndex(3)) = iAtom
252 11701 : atomShiftsInCubes(1,numAtoms,cubeIndex(1),cubeIndex(2),cubeIndex(3)) = i - FLOOR(atoms%taual(1,na))
253 11701 : atomShiftsInCubes(2,numAtoms,cubeIndex(1),cubeIndex(2),cubeIndex(3)) = j - FLOOR(atoms%taual(2,na))
254 11701 : atomShiftsInCubes(3,numAtoms,cubeIndex(1),cubeIndex(2),cubeIndex(3)) = k - FLOOR(atoms%taual(3,na))
255 46804 : posInCubes(:,numAtoms,cubeIndex(1),cubeIndex(2),cubeIndex(3)) = pos(:)
256 11701 : IF((i.EQ.0).AND.(j.EQ.0).AND.(k.EQ.0).AND.(na.EQ.1)) THEN
257 564 : refCubes(:,n) = cubeIndex(:)
258 564 : refPos(:,n) = pos(:)
259 : END IF
260 : END IF
261 : END DO
262 : END DO
263 : END DO
264 : END DO
265 : END DO
266 :
267 : ! 5. For each atom type in auxiliary unit cell select cube and collect shortest distances
268 : ! to other atoms in neighborhood
269 :
270 82 : maxSqrDist = cubeLength**2
271 246 : ALLOCATE(sqrDistances(8*maxCubeAtoms)) ! Formally 27, but 8 should be enough due to maxSqrDist
272 246 : ALLOCATE(neighborAtomRefs(8*maxCubeAtoms))
273 164 : ALLOCATE(neighborAtoms(8*maxCubeAtoms))
274 246 : ALLOCATE(neighborShifts(3,8*maxCubeAtoms))
275 164 : ALLOCATE(distIndexList(8*maxCubeAtoms))
276 :
277 223 : DO n = 1, atoms%ntype
278 564 : cubeIndex(:) = refCubes(:,n)
279 148117 : neighborAtomRefs = 0
280 141 : iNeighborAtom = 0
281 141 : identicalAtoms = 0
282 564 : DO i = cubeIndex(1) - 1, cubeIndex(1) + 1
283 1833 : DO j = cubeIndex(2) - 1, cubeIndex(2) + 1
284 5499 : DO k = cubeIndex(3) - 1, cubeIndex(3) + 1
285 13529 : DO iAtom = 1, numAtomsInCubes(i,j,k)
286 : currentDist = (refPos(1,n) - posInCubes(1,iAtom,i,j,k))**2 + &
287 : (refPos(2,n) - posInCubes(2,iAtom,i,j,k))**2 + &
288 8453 : (refPos(3,n) - posInCubes(3,iAtom,i,j,k))**2
289 12260 : IF (currentDist.LT.0.000001) THEN
290 141 : identicalAtoms = identicalAtoms + 1
291 8312 : ELSE IF (currentDist.LT.maxSqrDist) THEN
292 1316 : iNeighborAtom = iNeighborAtom + 1
293 1316 : neighborAtomRefs(iNeighborAtom) = atomRefsInCubes(iAtom,i,j,k)
294 1316 : neighborAtoms(iNeighborAtom) = atomsInCubes(iAtom,i,j,k)
295 5264 : neighborShifts(:,iNeighborAtom) = atomShiftsInCubes(:,iAtom,i,j,k)
296 1316 : sqrDistances(iNeighborAtom) = currentDist
297 : END IF
298 : END DO
299 : END DO
300 : END DO
301 : END DO
302 141 : IF (identicalAtoms.GT.1) THEN
303 0 : WRITE(*,*) 'Position: ', refPos(:,n)
304 0 : CALL juDFT_error("Too many atoms at same position.",calledby ="chkmt")
305 : END IF
306 141 : numNearestNeighbors(n) = MIN(maxCubeAtoms,iNeighborAtom)
307 141 : CALL sort(distIndexList(:iNeighborAtom),sqrDistances(:iNeighborAtom))
308 1539 : DO i = 1, numNearestNeighbors(n)
309 1316 : nearestNeighborRefsSorted(i,n) = neighborAtomRefs(distIndexList(i))
310 1316 : nearestNeighborAtomsSorted(i,n) = neighborAtoms(distIndexList(i))
311 5264 : nearestNeighborShiftsSorted(:,i,n) = neighborShifts(:,distIndexList(i))
312 1457 : nearestNeighborDists(i,n) = SQRT(sqrDistances(distIndexList(i)))
313 : END DO
314 : END DO
315 :
316 82 : WRITE(ounit,*) ''
317 223 : DO i = 1, atoms%ntype
318 223 : IF(numNearestNeighbors(i).GE.1) THEN
319 119 : nearestAtoms(i) = nearestNeighborRefsSorted(1,i)
320 119 : nearestAtomDists(i) = nearestNeighborDists(1,i)
321 119 : WRITE(ounit,'(a,i6,a)') 'Nearest atoms, shifts for reference of atom type ', i, ":"
322 119 : WRITE(ounit,'(a)') ' atom shifts'
323 119 : WRITE(ounit,'(4i6)') nearestNeighborAtomsSorted(1,i), nearestNeighborShiftsSorted(:,1,i)
324 747 : DO j = 2, numNearestNeighbors(i)
325 709 : IF ((nearestNeighborDists(j,i) - nearestAtomDists(i)).GT.0.0001) THEN
326 81 : numNextNearestNeighbors(i) = j - 1
327 81 : EXIT
328 : END IF
329 628 : IF (j.EQ.numNearestNeighbors(i)) numNextNearestNeighbors(i) = j
330 666 : WRITE(ounit,'(4i6)') nearestNeighborAtomsSorted(j,i), nearestNeighborShiftsSorted(:,j,i)
331 : END DO
332 : ELSE
333 22 : nearestAtoms(i) = -1
334 22 : nearestAtomDists(i) = 5000.0 * cubeLength
335 22 : numNextNearestNeighbors(i) = 0
336 : END IF
337 : END DO
338 :
339 :
340 82 : IF (PRESENT(l_gga)) THEN
341 : ! Sort distances and set MT radii for the atoms
342 :
343 0 : CALL sort(sortedDistList,nearestAtomDists)
344 0 : rmt1 = -1.0
345 0 : minRmts = -1.0
346 0 : DO i = 1, atoms%ntype
347 0 : typeA = sortedDistList(i)
348 0 : typeB = nearestAtoms(typeA)
349 0 : IF(typeB.LT.0) CYCLE
350 0 : dist = nearestAtomDists(typeA)
351 0 : IF (dist.LT.0.5) THEN
352 0 : WRITE (*,*) "Distance between atoms too small!"
353 0 : WRITE (*,*) "atom type A: ", typeA
354 0 : WRITE (*,*) "atom type B: ", typeB
355 0 : WRITE (*,*) "distance: ", dist
356 0 : CALL juDFT_error("Distance between atoms too small!",calledby ="chkmt")
357 : END IF
358 0 : sum_r = 1.0 / ( t_rmt(atoms%nz(typeA)) + t_rmt(atoms%nz(typeB)) )
359 0 : facA = t_rmt(atoms%nz(typeA)) * sum_r
360 0 : facB = t_rmt(atoms%nz(typeB)) * sum_r
361 : ! Note: The result of this section may be slightly different from the old version
362 : ! iff the nearest atom is another atom of the same type
363 0 : IF (minRmts(atoms%nz(typeA)).LT.0.0) THEN
364 0 : IF (minRmts(atoms%nz(typeB)).LT.0.0) THEN
365 0 : minRmts(atoms%nz(typeA)) = rmtFac * dist * facA
366 0 : minRmts(atoms%nz(typeB)) = rmtFac * dist * facB
367 : ELSE
368 0 : minRmts(atoms%nz(typeA)) = rmtFac * (dist - minRmts(atoms%nz(typeB)))
369 : END IF
370 0 : ELSE IF (minRmts(atoms%nz(typeB)).LT.0.0) THEN
371 0 : minRmts(atoms%nz(typeB)) = rmtFac * (dist - minRmts(atoms%nz(typeA)))
372 : END IF
373 0 : minRmts(atoms%nz(typeA)) = min(minRmts(atoms%nz(typeA)),rmtMaxDefault) ! limit already here
374 0 : minRmts(atoms%nz(typeB)) = min(minRmts(atoms%nz(typeB)),rmtMaxDefault) ! to a reasonable value
375 : END DO
376 :
377 : ! 6. Correct bad choices and set missing MT radii, vacuum distances, and other parameters
378 :
379 0 : DO i = 1, atoms%ntype
380 0 : IF((minRmts(atoms%nz(i)).LT.0.0).OR.(minRmts(atoms%nz(i)).GE.rmtMaxDefault)) THEN
381 0 : minRmts(atoms%nz(i)) = rmtMaxDefault
382 : END IF
383 0 : rmt1(i) = minRmts(atoms%nz(i))
384 : END DO
385 :
386 : ! NOTE: The result of this section may be slightly different from the old version
387 : ! iff the old version would enlarge a MT sphere at this point.
388 : ! Also the old version does not propagate the changes of the MT radii to all
389 : ! atoms with the same atomic number
390 0 : DO i = 1, atoms%ntype
391 0 : DO j = 1, numNearestNeighbors(i)
392 0 : k = nearestNeighborRefsSorted(j,i)
393 0 : IF (rmt1(i)+rmt1(k).GE.nearestNeighborDists(j,i)) THEN
394 0 : minRmts(atoms%nz(i)) = MIN(rmtFac*nearestNeighborDists(j,i)/2.0,MIN(rmt1(i),minRmts(atoms%nz(i))))
395 0 : minRmts(atoms%nz(k)) = MIN(rmtFac*nearestNeighborDists(j,i)/2.0,MIN(rmt1(k),minRmts(atoms%nz(k))))
396 : END IF
397 : END DO
398 : END DO
399 :
400 0 : DO i = 1, atoms%ntype
401 0 : rmt1(i) = minRmts(atoms%nz(i))
402 : END DO
403 :
404 0 : WRITE (oUnit,*) '----------------------------------------------------'
405 0 : WRITE (oUnit,*) 'Suggested values for input: '
406 0 : WRITE (oUnit,*)
407 :
408 0 : dvac1 = 0.0
409 0 : IF (input%film) THEN
410 0 : iAtom = 0
411 0 : DO i = 1, atoms%ntype
412 0 : DO na = 1, atoms%neq(i)
413 0 : iAtom = iAtom + 1
414 :
415 0 : dvac1 = MAX(dvac1, ABS(atoms%pos(3,iAtom))+rmt1(i))
416 :
417 : END DO
418 : END DO
419 0 : dvac1 = 2.0 * (dvac1+0.3)
420 0 : dtild = dvac1 + 1.5 * MAXVAL(rmt1(:))
421 0 : WRITE (oUnit,'("vacuum distance dvac =",f10.5)') dvac1
422 0 : WRITE (oUnit,'("extra vac.dist. dtild=",f10.5)') dtild
423 : END IF
424 :
425 0 : rkm = 0.0
426 0 : WRITE (oUnit,*) 'Atom Z lmax jri rmt dx'
427 0 : DO n = 1, atoms%ntype
428 0 : IF (rmt1(n).LT.1.8) THEN
429 0 : lmax11 = 6
430 0 : ELSE IF (rmt1(n).LT.2.4) THEN
431 0 : lmax11 = 8
432 : ELSE
433 0 : lmax11 = 10
434 : END IF
435 0 : IF (l_gga) THEN
436 0 : jri11 = NINT(330*rmt1(n))
437 : ELSE
438 0 : jri11 = NINT(220*rmt1(n))
439 : END IF
440 0 : jri11 = NINT(jri11*0.5) * 2 + 1
441 0 : IF (atoms%nz(n) > 0) THEN
442 0 : dx11 = LOG(3200*atoms%nz(n)*rmt1(n))/(jri11-1)
443 : ELSE
444 0 : dx11 = LOG(3200*rmt1(n))/(jri11-1)
445 : ENDIF
446 0 : rkm = MAX(rkm, lmax11/rmt1(n))
447 0 : WRITE (oUnit,'(a3,i3,2i5,2f10.6)') noel(n),atoms%nz(n),lmax11,jri11,rmt1(n),dx11
448 0 : dx1(n) = dx11
449 0 : lmax1(n) = lmax11
450 0 : jri1(n) = jri11
451 : END DO
452 0 : WRITE (oUnit,'("k_max =",f8.5)') rkm
453 0 : WRITE (oUnit,'("G_max =",f8.5)') 3*rkm
454 0 : kmax = rkm
455 : ENDIF
456 :
457 : ! 7. Test old MT radii
458 :
459 82 : IF (l_test) THEN
460 82 : iAtom = 0
461 223 : DO i = 1, atoms%ntype
462 1457 : DO j = 1, numNearestNeighbors(i)
463 1316 : k = nearestNeighborRefsSorted(j,i)
464 1457 : IF (atoms%rmt(i)+atoms%rmt(k).GE.nearestNeighborDists(j,i)) THEN
465 0 : error = .TRUE.
466 0 : IF (PRESENT(overlap)) overlap(i,k)=atoms%rmt(i)+atoms%rmt(k)-nearestNeighborDists(j,i)
467 0 : WRITE(oUnit,240) i,k,nearestNeighborDists(j,i),atoms%rmt(i),atoms%rmt(k)
468 0 : IF (.NOT.PRESENT(overlap)) CALL juDFT_error("Overlapping MT-radii of two neighbours detected.",calledby ="chkmt")
469 : END IF
470 : END DO
471 223 : IF (input%film) THEN
472 58 : DO na = 1, atoms%neq(i)
473 35 : iAtom = iAtom + 1
474 :
475 35 : IF (((atoms%pos(3,iAtom)+atoms%rmt(i)).GT. vacuum%dvac/2.).OR.&
476 23 : ((atoms%pos(3,iAtom)-atoms%rmt(i)).LT.-vacuum%dvac/2.)) THEN
477 0 : error=.TRUE.
478 0 : WRITE(oUnit,241) i ,na
479 0 : IF (PRESENT(overlap)) overlap(0,i)=MAX(overlap(0,i),abs(atoms%pos(3,iAtom))+atoms%rmt(i)-vacuum%dvac/2.)
480 0 : WRITE(oUnit,*) atoms%pos(3,iAtom),atoms%rmt(i),vacuum%dvac/2.
481 0 : IF (.NOT.PRESENT(overlap)) CALL juDFT_error("MT is partly inside the vacuum.",calledby ="chkmt")
482 : ENDIF
483 :
484 : END DO
485 : END IF
486 : END DO
487 82 : IF (error.AND..NOT.PRESENT(overlap)) CALL juDFT_error("Error checking M.T. radii",calledby ="chkmt")
488 : END IF
489 :
490 82 : DEALLOCATE(nearestNeighborRefsSorted,nearestNeighborAtomsSorted,nearestNeighborShiftsSorted,numNearestNeighbors,numNextNearestNeighbors,nearestNeighborDists)
491 82 : DEALLOCATE(distIndexList,neighborAtomRefs,neighborAtoms,neighborShifts,sqrDistances)
492 82 : DEALLOCATE(numAtomsInCubes,atomRefsInCubes,atomsInCubes,atomShiftsInCubes,posInCubes,refCubes,refPos)
493 :
494 : 240 FORMAT('Error in muffin tin radii pair (',i5,',',i5,'):',3f10.5)
495 : 241 FORMAT(' error: atom ',i3,' # ',i3,'reaches out into vaccuum')
496 :
497 166 : END SUBROUTINE chkmt
498 : END MODULE m_chkmt
|
