LCOV - code coverage report
Current view: top level - global - chkmt.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 215 241 89.2 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

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

Generated by: LCOV version 1.13