From 2712fa6496f039f945573cc44050db1bb5d82321 Mon Sep 17 00:00:00 2001 From: Gregor Michalicek Date: Wed, 27 Apr 2016 18:28:54 +0200 Subject: [PATCH] Rewrote chkmt.f90. This commit introduces the new chkmt.f90. Besides some exceptions this new routine should yield the same results as the old one. With a runtime scaling of O(N_atoms) it overcomes the bad O(N_atoms^4) scaling of the old pendant. --- cmake/Files_and_Targets.txt | 3 +- global/chkmt.f90 | 591 +++++++++++++++++++++++------------- main/fleur_init.F90 | 1 + 3 files changed, 386 insertions(+), 209 deletions(-) diff --git a/cmake/Files_and_Targets.txt b/cmake/Files_and_Targets.txt index 64fb7ce2..862939bd 100644 --- a/cmake/Files_and_Targets.txt +++ b/cmake/Files_and_Targets.txt @@ -36,7 +36,8 @@ global/constants.f inpgen/element.f inpgen/atom_input.f inpgen/crystal.f inpgen/ inpgen/atom_sym.f inpgen/generator.f inpgen/read_record.f inpgen/soc_or_ssdw.f inpgen/symproperties.f inpgen/bravais_symm.f inpgen/set_atom_core.f inpgen/spg_gen.f global/triang.f inpgen/closure.f inpgen/lapw_input.f inpgen/struct_input.f inpgen/write_struct.f - global/hybridmix.f io/calculator.f global/ss_sym.f global/soc_sym.f math/inv3.f io/rw_symfile.f +global/hybridmix.f io/calculator.f global/ss_sym.f global/soc_sym.f math/inv3.f io/rw_symfile.f +global/sort.f ) set(inpgen_F90 io/xsf_io.f90 diff --git a/global/chkmt.f90 b/global/chkmt.f90 index 2e82c4dc..52527431 100644 --- a/global/chkmt.f90 +++ b/global/chkmt.f90 @@ -1,12 +1,12 @@ MODULE m_chkmt - use m_juDFT + USE m_juDFT private public chkmt - INTEGER,PARAMETER:: MAX_CHECK=100 !--------------------------------------------------------------------- ! Check muffin tin radii and determine a reasonable choice for MTRs. ! Derive also other parameters for the input file, to provide some -! help in the out-file. gb`02 +! help in the out-file. +! GM'16 !--------------------------------------------------------------------- CONTAINS SUBROUTINE chkmt(& @@ -15,6 +15,8 @@ & kmax,dtild,dvac1,lmax1,jri1,rmt1,dx1) USE m_types + USE m_sort + USE m_inv3 IMPLICIT NONE ! .. ! .. Scalar Arguments .. @@ -32,229 +34,402 @@ REAL, INTENT (OUT) :: rmt1(atoms%ntype),dx1(atoms%ntype) ! .. ! .. Local Scalars .. - INTEGER na,n,nna,nn,n1,nn1,k1,k2,k3 - INTEGER i,j,jri11,lmax11 - REAL sum,sss,xmin,dx11,rkm,fac,sum_r,fac_1,fac_2 - LOGICAL error - REAL distvec(3),vec(3,26) + INTEGER na,n + INTEGER i,j,k,l,jri11,lmax11 + INTEGER maxCubeAtoms, iAtom, numAtoms, iNeighborAtom, identicalAtoms + INTEGER typeA, typeB + REAL dx11,rkm,sum_r,facA,facB + REAL rmtMax, rmtMin, rmtMaxDefault, rmtDelta + REAL rmtFac, cubeLength, amatAuxDet + REAL maxSqrDist, dist, currentDist + LOGICAL error, outOfBounds ! .. ! .. Local Arrays .. - INTEGER minni(2) - REAL dist(atoms%ntype,atoms%ntype),dist1(atoms%ntype,atoms%ntype),t_rmt(0:103) + REAL t_rmt(0:103), minRmts(0:103) + REAL amatAux(3,3), invAmatAux(3,3) + REAL taualAux(3,atoms%nat), posAux(3,atoms%nat) + REAL minPos(3), maxPos(3), pos(3), point(3), realCellPos(3) + REAL offsetPos(3) + REAL nearestAtomDists(atoms%ntype) + INTEGER nearestAtoms(atoms%ntype) + INTEGER sortedDistList(atoms%ntype) + INTEGER minCubeIndex(3), maxCubeIndex(3), cubeIndex(3) + INTEGER minCellIndices(3), maxCellIndices(3) - IF (l_test.and.atoms%ntype>MAX_CHECK) THEN - write(6,*) "No test for MT radii performed" - write(6,*) "Change MAX_CHECK in chkmt if needed" - return - ENDIF -! -! typical muffin-tin radii -! - i=0 - DO k1=-1,1 - DO k2=-1,1 - DO k3=-1,1 - if (k1==0.and.k2==0.and.k3==0) cycle - i=i+1 - vec(:,i)=matmul(cell%amat,(/k1,k2,k3/)) - ENDDO - ENDDO - ENDDO - t_rmt(0:103) = 2.3 ! default value - t_rmt(1) = 1.0 ; t_rmt(5:9) = 1.3 ; t_rmt(16:17) = 1.8 + INTEGER, ALLOCATABLE :: numAtomsInCubes(:,:,:) + INTEGER, ALLOCATABLE :: atomRefsInCubes(:,:,:,:) + INTEGER, ALLOCATABLE :: refCubes(:,:) + INTEGER, ALLOCATABLE :: nearestNeighbors(:,:) + INTEGER, ALLOCATABLE :: numNearestNeighbors(:) + INTEGER, ALLOCATABLE :: neighborAtoms(:) + INTEGER, ALLOCATABLE :: distIndexList(:) + REAL, ALLOCATABLE :: posInCubes(:,:,:,:,:) + REAL, ALLOCATABLE :: refPos(:,:) + REAL, ALLOCATABLE :: nearestNeighborDists(:,:) + REAL, ALLOCATABLE :: sqrDistances(:) + +! Plan for this routine: +! 0. Do initializations and set some constants +! 1. Locally replace unit cell by an auxiliary unit cell with: +! a) all atoms within the unit cell +! b) basis vectors obtained by lattice reduction of the original cell. +! [not in 1st (this) version of routine. Can be implemented with the LLL algorithm when needed.] +! 2. Get minimal and maximal coordinates within auxiliary unit cell +! 3. Construct mesh of cubes covering the auxiliary unit cell and a boundary of width 2*rmtMax + rmtDelta +! 4. Fill mesh of cubes with atoms +! a) Store atoms in cubes and representative cube for each atom type +! 5. For each atom in auxiliary unit cell select cube and collect shortest distances to other atoms in neighborhood +! a) Sort distances and set MT radii for the atoms +! 6. Correct bad choices and set missing MT radii, vacuum distances, and other parameters +! 7. Test old MT radii - error=.false. - dist(:,:) = 9.99e19 - !$OMP PARALLEL DO DEFAULT(NONE) & - !$OMP SHARED(atoms,dist,error,l_test,input,oneD,vacuum,vec)& - !$OMP PRIVATE(n,n1,na,nn,nn1,nna,sss,distvec,k1) - DO n=1,atoms%ntype - DO n1=1,atoms%neq(n) - na=n1+sum(atoms%neq(:n-1)) -! -! check distance to other atoms: -! - DO nn=n,atoms%ntype - DO nn1=1,atoms%neq(nn) - nna=nn1+sum(atoms%neq(:nn-1)) - sss=1.E19 - distvec=atoms%pos(:,na)-atoms%pos(:,nna) - if (nna.ne.na) sss=distvec(1)**2+distvec(2)**2+distvec(3)**2 - DO k1=1,26 - sss=min(sss,(distvec(1)+vec(1,k1))**2+(distvec(2)+vec(2,k1))**2+(distvec(3)+vec(3,k1))**2) - ENDDO - dist(n,nn) = min( dist(n,nn),sqrt(sss) ) - dist(nn,n) = dist(n,nn) - IF ( sss.LE.(atoms%rmt(nn)+atoms%rmt(n))**2 ) THEN - error=.true. - IF (l_test)& - & WRITE(6,240) nn,nn1,(atoms%pos(i,nna),i=1,3),atoms%rmt(nn),& - & n ,n1 ,(atoms%pos(i,na),i=1,3),atoms%rmt(n ) - ENDIF - - ENDDO ! nn1 - ENDDO ! nn -! -! distance to vacuum -! - IF (input%film) THEN - IF (oneD%odd%d1) THEN - IF ((sqrt(atoms%pos(1,na)**2+atoms%pos(2,na)**2)+& - & atoms%rmt(n)).GT.vacuum%dvac/2.) THEN - error=.true. - WRITE(6,241) n ,n1 - WRITE(6,*) sqrt(atoms%pos(1,na)**2+atoms%pos(2,na)**2),& - & atoms%rmt(n),vacuum%dvac/2. - END IF - ELSE - IF ( ( (atoms%pos(3,na)+atoms%rmt(n) ).GT. vacuum%dvac/2.).OR.& - & ( (atoms%pos(3,na)-atoms%rmt(n) ).LT.-vacuum%dvac/2.) ) THEN - error=.true. - WRITE(6,241) n ,n1 - WRITE(6,*) atoms%pos(3,na),atoms%rmt(n),vacuum%dvac/2. - ENDIF - ENDIF - END IF - ENDDO - ENDDO - !$OMP END PARALLEL DO - -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -! DO n = 1, ntype -! WRITE (*,'(12f12.6)') dist(:,n) -! ENDDO - dist1 = dist - rmt1(:) = 999. - WRITE (6,*) '----------------------------------------------------' - WRITE (6,*) 'Suggested values for input: ' - WRITE (6,*) +! 0. Do initializations and set some constants + + rmtMaxDefault = 2.8 + rmtMax = rmtMaxDefault + rmtMin = 1.0 + IF (l_test) THEN + rmtMax = MAX(rmtMax,MAXVAL(atoms%rmt(:))) + rmtMin = MIN(rmtMin,MINVAL(atoms%rmt(:))) + END IF + rmtDelta = 0.3 IF (input%film) THEN - fac = 0.95 + rmtFac = 0.95 ELSE - fac = 0.975 + rmtFac = 0.975 ENDIF + t_rmt(0:103) = 2.3 ! default value + t_rmt(1) = 1.0 ; t_rmt(5:9) = 1.3 ; t_rmt(16:17) = 1.8 + cubeLength = 2*rmtMax+rmtDelta + maxCubeAtoms = (FLOOR(cubeLength / (0.8*2.0*rmtMin)) + 1)**3 + error = .FALSE. - minni = minloc(dist) ! minni(1) and minni(2) are the indices of the closest atoms - xmin = minval(dist) ! xmin is their distance - - DO WHILE ( xmin < 999.0 ) - - sum_r = 1.0 / ( t_rmt(atoms%nz(minni(1))) + t_rmt(atoms%nz(minni(2))) ) - fac_1 = t_rmt(atoms%nz(minni(1))) * sum_r - fac_2 = t_rmt(atoms%nz(minni(2))) * sum_r - - IF (rmt1(minni(1)) > 990.) THEN ! if not set, determine MTR - IF (rmt1(minni(2)) > 990.) THEN ! both not set, choose in between - rmt1(minni(1)) = fac * xmin * fac_1 ! / 2 - rmt1(minni(2)) = fac * xmin * fac_2 ! / 2 - ELSE - rmt1(minni(1)) = fac * ( xmin - rmt1(minni(2)) ) - IF (2*rmt1(minni(1)).GT.dist(minni(1),minni(1))) THEN - rmt1(minni(1)) = fac * dist(minni(1),minni(1)) / 2 - ENDIF - ENDIF - ELSEIF (rmt1(minni(2)) > 990.) THEN - rmt1(minni(2)) = fac * ( xmin - rmt1(minni(1)) ) - IF (2*rmt1(minni(2)).GT.dist(minni(2),minni(2))) THEN - rmt1(minni(2)) = fac * dist(minni(2),minni(2)) / 2 - ENDIF - ENDIF - - dist(minni(1),minni(1)) = 999.0 - dist(minni(2),minni(1)) = 999.0 - dist(minni(1),minni(2)) = 999.0 - dist(minni(2),minni(2)) = 999.0 - - DO j = 1, 2 - DO n = 1, atoms%ntype - IF (atoms%nz(n) == atoms%nz(minni(j))) THEN - IF (rmt1(n) > 990.) THEN - rmt1(n) = rmt1(minni(j)) - ENDIF - ENDIF - ENDDO - ENDDO - - minni = minloc(dist) - xmin = minval(dist) - - ENDDO +! 1. For the 1st version the auxiliary unit cell is just a copy of the original unit cell with +! all atoms within the cell. - dvac1 = 0.0 - rkm = 0.0 - na = 0 - WRITE (6,230) - DO n= 1,atoms%ntype -! -!--> determine M.T. radii -! - DO j= 1,atoms%ntype - dist(j,n) = rmt1(n)+rmt1(j) - IF ( dist1(j,n)-dist(j,n) < 0.0 ) THEN - WRITE(*,*) j,n,dist1(j,n)-dist(j,n) - rmt1(n) = fac * dist1(j,n) / 2. - rmt1(j) = fac * dist1(j,n) / 2. - ENDIF - ENDDO - IF (input%film) THEN - DO nn = 1, atoms%neq(n) - na = na + 1 - IF (oneD%odd%d1) THEN - dvac1 = max( dvac1, sqrt(atoms%pos(1,na)**2+atoms%pos(2,na)**2)& - & +rmt1(n) ) + DO i = 1, 3 + DO j = 1, 3 + amatAux(i,j) = cell%amat(i,j) + END DO + END DO + + DO i = 1, atoms%nat + taualAux(1,i) = atoms%taual(1,i) - FLOOR(atoms%taual(1,i)) + taualAux(2,i) = atoms%taual(2,i) - FLOOR(atoms%taual(2,i)) + taualAux(3,i) = atoms%taual(3,i) - FLOOR(atoms%taual(3,i)) + posAux(:,i) = matmul(amatAux,taualAux(:,i)) + END DO + +! 2. Get minimal and maximal coordinates for auxiliary unit cell + + minPos = 0.0 + maxPos = 0.0 + + DO i = 0, 1 + DO j = 0, 1 + DO k = 0, 1 + DO l = 1, 3 + pos(l) = i*amatAux(l,1) + j*amatAux(l,2) + k*amatAux(l,3) + IF (pos(l).GT.maxPos(l)) maxPos(l) = pos(l) + IF (pos(l).LT.minPos(l)) minPos(l) = pos(l) + END DO + END DO + END DO + END DO + +! 3. Construct cube mesh: +! In each dimension cube i covers the interval from i*cubeLength to (i+1)*cubeLength +! Each cube may cover up to maxCubeAtoms atoms. This should be set to a save size. + + DO i = 1, 3 + minPos(i) = minPos(i) - cubeLength + maxPos(i) = maxPos(i) + cubeLength + minCubeIndex(i) = FLOOR(minPos(i)/cubeLength) + maxCubeIndex(i) = CEILING(maxPos(i)/cubeLength) + END DO + + ALLOCATE (numAtomsInCubes(minCubeIndex(1):maxCubeIndex(1),& + minCubeIndex(2):maxCubeIndex(2),& + minCubeIndex(3):maxCubeIndex(3))) + ALLOCATE (atomRefsInCubes(maxCubeAtoms,minCubeIndex(1):maxCubeIndex(1),& + minCubeIndex(2):maxCubeIndex(2),& + minCubeIndex(3):maxCubeIndex(3))) + ALLOCATE (posInCubes(3,maxCubeAtoms,minCubeIndex(1):maxCubeIndex(1),& + minCubeIndex(2):maxCubeIndex(2),& + minCubeIndex(3):maxCubeIndex(3))) + ALLOCATE (refCubes(3,atoms%ntype),refPos(3,atoms%ntype)) + ALLOCATE (nearestNeighbors(maxCubeAtoms,atoms%ntype),numNearestNeighbors(atoms%ntype)) + ALLOCATE (nearestNeighborDists(maxCubeAtoms,atoms%ntype)) + + numAtomsInCubes = 0 + +! 4. Fill mesh of cubes with atoms +! First obtain minimal and maximal indices for relevant unit cells + + minCellIndices = 0 + maxCellIndices = 0 + + CALL inv3(amatAux,invAmatAux,amatAuxDet) + + DO i = 0, 1 + DO j = 0, 1 + DO k = 0, 1 + point(:) = minPos(:) + IF(i.EQ.1) point(1) = maxPos(1) + IF(j.EQ.1) point(2) = maxPos(2) + IF(k.EQ.1) point(3) = maxPos(3) + realCellPos(:) = matmul(invAmatAux,point(:)) + DO l = 1, 3 + IF(minCellIndices(l).GT.realCellPos(l)) THEN + minCellIndices(l) = FLOOR(realCellPos(l)) + END IF + IF(maxCellIndices(l).LT.realCellPos(l)) THEN + maxCellIndices(l) = FLOOR(realCellPos(l)) ! Is 'FLOOR' enough? + END IF + END DO + END DO + END DO + END DO + +! Store atoms in cubes and representative cube for each atom type + + DO i = minCellIndices(1), maxCellIndices(1) + DO j = minCellIndices(2), maxCellIndices(2) + DO k = minCellIndices(3), maxCellIndices(3) + DO l = 1, 3 + offsetPos(l) = i*amatAux(l,1) + j*amatAux(l,2) + k*amatAux(l,3) + END DO + iAtom = 0 + DO n = 1, atoms%ntype + DO na = 1, atoms%neq(n) + iAtom = iAtom + 1 + pos(:) = posAux(:,iAtom) + offsetPos(:) + outOfBounds = .FALSE. + DO l = 1, 3 + cubeIndex(l) = FLOOR(pos(l)/cubeLength) + IF(cubeIndex(l).LT.minCubeIndex(l)) outOfBounds = .TRUE. + IF(cubeIndex(l).GT.maxCubeIndex(l)) outOfBounds = .TRUE. + END DO + IF(.NOT.outOfBounds) THEN + numAtomsInCubes(cubeIndex(1),cubeIndex(2),cubeIndex(3)) = & + numAtomsInCubes(cubeIndex(1),cubeIndex(2),cubeIndex(3)) + 1 + numAtoms = numAtomsInCubes(cubeIndex(1),cubeIndex(2),cubeIndex(3)) + IF(numAtoms.GT.maxCubeAtoms) THEN + STOP 'ERROR: maxCubeAtoms is not large enough in chkmt.' + END IF + atomRefsInCubes(numAtoms,cubeIndex(1),cubeIndex(2),cubeIndex(3)) = n + posInCubes(:,numAtoms,cubeIndex(1),cubeIndex(2),cubeIndex(3)) = pos(:) + IF((i.EQ.0).AND.(j.EQ.0).AND.(k.EQ.0).AND.(na.EQ.1)) THEN + refCubes(:,n) = cubeIndex(:) + refPos(:,n) = pos(:) + END IF + END IF + END DO + END DO + END DO + END DO + END DO + +! 5. For each atom type in auxiliary unit cell select cube and collect shortest distances +! to other atoms in neighborhood + + maxSqrDist = cubeLength**2 + ALLOCATE(sqrDistances(8*maxCubeAtoms)) ! Formally 27, but 8 should be enough due to maxSqrDist + ALLOCATE(neighborAtoms(8*maxCubeAtoms)) + ALLOCATE(distIndexList(8*maxCubeAtoms)) + + DO n = 1, atoms%ntype + cubeIndex(:) = refCubes(:,n) + neighborAtoms = 0 + iNeighborAtom = 0 + identicalAtoms = 0 + DO i = cubeIndex(1) - 1, cubeIndex(1) + 1 + DO j = cubeIndex(2) - 1, cubeIndex(2) + 1 + DO k = cubeIndex(3) - 1, cubeIndex(3) + 1 + DO iAtom = 1, numAtomsInCubes(i,j,k) + currentDist = (refPos(1,n) - posInCubes(1,iAtom,i,j,k))**2 + & + (refPos(2,n) - posInCubes(2,iAtom,i,j,k))**2 + & + (refPos(3,n) - posInCubes(3,iAtom,i,j,k))**2 + IF (currentDist.LT.0.000001) THEN + identicalAtoms = identicalAtoms + 1 + ELSE IF (currentDist.LT.maxSqrDist) THEN + iNeighborAtom = iNeighborAtom + 1 + neighborAtoms(iNeighborAtom) = atomRefsInCubes(iAtom,i,j,k) + sqrDistances(iNeighborAtom) = currentDist + END IF + END DO + END DO + END DO + END DO + IF (identicalAtoms.GT.1) THEN + WRITE(*,*) 'Position: ', refPos(:,n) + STOP 'Error: Too many atoms at same position.' + END IF + numNearestNeighbors(n) = MIN(maxCubeAtoms,iNeighborAtom) + CALL sort(iNeighborAtom,sqrDistances,distIndexList) + DO i = 1, numNearestNeighbors(n) + nearestNeighbors(i,n) = neighborAtoms(distIndexList(i)) + nearestNeighborDists(i,n) = SQRT(sqrDistances(distIndexList(i))) + END DO + END DO + + DO i = 1, atoms%ntype + IF(numNearestNeighbors(i).GE.1) THEN + nearestAtoms(i) = nearestNeighbors(1,i) + nearestAtomDists(i) = nearestNeighborDists(1,i) + ELSE + nearestAtoms(i) = -1 + nearestAtomDists(i) = 5000.0 * cubeLength + END IF + END DO + +! Sort distances and set MT radii for the atoms + + CALL sort(atoms%ntype,nearestAtomDists,sortedDistList) + rmt1 = -1.0 + minRmts = -1.0 + DO i = 1, atoms%ntype + typeA = sortedDistList(i) + typeB = nearestAtoms(typeA) + IF(typeB.LT.0) CYCLE + dist = nearestAtomDists(typeA) + sum_r = 1.0 / ( t_rmt(atoms%nz(typeA)) + t_rmt(atoms%nz(typeB)) ) + facA = t_rmt(atoms%nz(typeA)) * sum_r + facB = t_rmt(atoms%nz(typeB)) * sum_r + ! Note: The result of this section may be slightly different from the old version + ! iff the nearest atom is another atom of the same type + IF (minRmts(atoms%nz(typeA)).LT.0.0) THEN + IF (minRmts(atoms%nz(typeB)).LT.0.0) THEN + minRmts(atoms%nz(typeA)) = rmtFac * dist * facA + minRmts(atoms%nz(typeB)) = rmtFac * dist * facB ELSE - dvac1 = max( dvac1, abs(atoms%pos(3,na))+rmt1(n) ) + minRmts(atoms%nz(typeA)) = rmtFac * (dist - minRmts(atoms%nz(typeB))) END IF - ENDDO - ENDIF -! -!--> calculate jri1, dx1 and lmax1 -! - IF (rmt1(n).LT.1.8) THEN - lmax11 = 6 - ELSEIF (rmt1(n).LT.2.4) THEN - lmax11 = 8 - ELSEIF (rmt1(n).LT.2.8) THEN - lmax11 = 10 - ELSE - WRITE (6,'("Atom Nr.",i3,"( ",a3,") has a M.T. radius of",& - & f8.4)') n,noel(n),rmt1(n) - WRITE (6,'("that was truncated to 2.8")') - rmt1(n) = 2.8 - lmax11 = 10 - ENDIF - IF (l_gga) THEN - jri11 = nint( 330*rmt1(n) ) - ELSE - jri11 = nint( 220*rmt1(n) ) - ENDIF - jri11 = nint( jri11*0.5 ) * 2 + 1 - dx11 = log(3200*atoms%nz(n)*rmt1(n))/(jri11-1) - rkm = max( rkm , lmax11/rmt1(n) ) - - WRITE (6,9070) noel(n),atoms%nz(n),lmax11,jri11,rmt1(n),dx11 - 9070 FORMAT (a3,i3,2i5,2f10.6) - dx1(n) = dx11 ; lmax1(n) = lmax11 ; jri1(n) = jri11 - - ENDDO ! loop over atom types - 230 FORMAT ('Atom Z lmax jri rmt dx') + ELSE IF (minRmts(atoms%nz(typeB)).LT.0.0) THEN + minRmts(atoms%nz(typeB)) = rmtFac * (dist - minRmts(atoms%nz(typeA))) + END IF + END DO + +! 6. Correct bad choices and set missing MT radii, vacuum distances, and other parameters + + DO i = 1, atoms%ntype + IF((minRmts(atoms%nz(i)).LT.0.0).OR.(minRmts(atoms%nz(i)).GE.rmtMaxDefault)) THEN + minRmts(atoms%nz(i)) = rmtMaxDefault + END IF + rmt1(i) = minRmts(atoms%nz(i)) + END DO + + ! NOTE: The result of this section may be slightly different from the old version + ! iff the old version would enlarge a MT sphere at this point. + ! Also the old version does not propagate the changes of the MT radii to all + ! atoms with the same atomic number + DO i = 1, atoms%ntype + DO j = 1, numNearestNeighbors(i) + k = nearestNeighbors(j,i) + IF (rmt1(i)+rmt1(k).GE.nearestNeighborDists(j,i)) THEN + minRmts(atoms%nz(i)) = MIN(rmtFac*nearestNeighborDists(j,i)/2.0,MIN(rmt1(i),minRmts(atoms%nz(i)))) + minRmts(atoms%nz(k)) = MIN(rmtFac*nearestNeighborDists(j,i)/2.0,MIN(rmt1(k),minRmts(atoms%nz(k)))) + END IF + END DO + END DO + + DO i = 1, atoms%ntype + rmt1(i) = minRmts(atoms%nz(i)) + END DO + + WRITE (6,*) '----------------------------------------------------' + WRITE (6,*) 'Suggested values for input: ' + WRITE (6,*) + dvac1 = 0.0 IF (input%film) THEN - dvac1 = 2* (dvac1+0.3) - dtild = dvac1 + 1.5 * maxval( rmt1(:) ) - WRITE (6,'("vacuum distance dvac =",f10.5)') dvac1 - WRITE (6,'("extra vac.dist. dtild=",f10.5)') dtild - ENDIF + iAtom = 0 + DO i = 1, atoms%ntype + DO na = 1, atoms%neq(i) + iAtom = iAtom + 1 + IF (oneD%odd%d1) THEN + dvac1 = MAX(dvac1, SQRT(atoms%pos(1,iAtom)**2+atoms%pos(2,iAtom)**2)+rmt1(i)) + ELSE + dvac1 = MAX(dvac1, ABS(atoms%pos(3,iAtom))+rmt1(i)) + END IF + END DO + END DO + dvac1 = 2.0 * (dvac1+0.3) + dtild = dvac1 + 1.5 * MAXVAL(rmt1(:)) + WRITE (6,'("vacuum distance dvac =",f10.5)') dvac1 + WRITE (6,'("extra vac.dist. dtild=",f10.5)') dtild + END IF + + rkm = 0.0 + WRITE (6,*) 'Atom Z lmax jri rmt dx' + DO n = 1, atoms%ntype + IF (rmt1(n).LT.1.8) THEN + lmax11 = 6 + ELSE IF (rmt1(n).LT.2.4) THEN + lmax11 = 8 + ELSE + lmax11 = 10 + END IF + IF (l_gga) THEN + jri11 = NINT(330*rmt1(n)) + ELSE + jri11 = NINT(220*rmt1(n)) + END IF + jri11 = NINT(jri11*0.5) * 2 + 1 + dx11 = LOG(3200*atoms%nz(n)*rmt1(n))/(jri11-1) + rkm = MAX(rkm, lmax11/rmt1(n)) + WRITE (6,'(a3,i3,2i5,2f10.6)') noel(n),atoms%nz(n),lmax11,jri11,rmt1(n),dx11 + dx1(n) = dx11 + lmax1(n) = lmax11 + jri1(n) = jri11 + END DO WRITE (6,'("k_max =",f8.5)') rkm WRITE (6,'("G_max =",f8.5)') 3*rkm kmax = rkm - WRITE (6,*) '----------------------------------------------------' - - IF ( error.AND.l_test ) CALL juDFT_error& - & ("Error checking M.T. radii",calledby ="chkmt") - 240 FORMAT(' error in muffin tin radii , pair ',2(/,2i5,4f10.5)) +! 7. Test old MT radii + + IF (l_test) THEN + iAtom = 0 + DO i = 1, atoms%ntype + DO j = 1, numNearestNeighbors(i) + k = nearestNeighbors(j,i) + IF (atoms%rmt(i)+atoms%rmt(k).GE.nearestNeighborDists(j,i)) THEN + error = .TRUE. + WRITE(6,240) i,k,nearestNeighborDists(j,i),atoms%rmt(i),atoms%rmt(k) + END IF + END DO + IF (input%film) THEN + DO na = 1, atoms%neq(i) + iAtom = iAtom + 1 + IF (oneD%odd%d1) THEN + IF ((sqrt(atoms%pos(1,iAtom)**2+atoms%pos(2,iAtom)**2)+& + atoms%rmt(i)).GT.vacuum%dvac/2.) THEN + error=.TRUE. + WRITE(6,241) i ,na + WRITE(6,*) sqrt(atoms%pos(1,iAtom)**2+atoms%pos(2,iAtom)**2),& + atoms%rmt(i),vacuum%dvac/2. + END IF + ELSE + IF (((atoms%pos(3,iAtom)+atoms%rmt(i)).GT. vacuum%dvac/2.).OR.& + ((atoms%pos(3,iAtom)-atoms%rmt(i)).LT.-vacuum%dvac/2.)) THEN + error=.TRUE. + WRITE(6,241) i ,na + WRITE(6,*) atoms%pos(3,iAtom),atoms%rmt(i),vacuum%dvac/2. + ENDIF + ENDIF + END DO + END IF + END DO + IF (error) CALL juDFT_error("Error checking M.T. radii",calledby ="chkmt") + END IF + + DEALLOCATE(nearestNeighbors,numNearestNeighbors,nearestNeighborDists) + DEALLOCATE(distIndexList,neighborAtoms,sqrDistances) + DEALLOCATE(numAtomsInCubes,atomRefsInCubes,posInCubes,refCubes,refPos) + + 240 FORMAT('Error in muffin tin radii pair (',i5,',',i5,'):',3f10.5) 241 FORMAT(' error: atom ',i3,' # ',i3,'reaches out into vaccuum') END SUBROUTINE chkmt diff --git a/main/fleur_init.F90 b/main/fleur_init.F90 index 2ffed66b..7e787b4e 100644 --- a/main/fleur_init.F90 +++ b/main/fleur_init.F90 @@ -96,6 +96,7 @@ oneD%odd%k3 = 0 ; oneD%odd%M =0 ; oneD%odd%nn2d = 1 ENDIF !-odim + atoms%nat = atoms%natd ! This is preliminary. The value of nat changes later. ALLOCATE ( atoms%nz(atoms%ntypd),atoms%relax(3,atoms%ntypd),atoms%nlhtyp(atoms%ntype)) ALLOCATE ( sphhar%clnu(sphhar%memd,0:sphhar%nlhd,sphhar%ntypsd),stars%ustep(stars%n3d) ) ALLOCATE ( stars%ig(-stars%k1d:stars%k1d,-stars%k2d:stars%k2d,-stars%k3d:stars%k3d),stars%ig2(stars%n3d),stars%igz(stars%n3d) ) -- GitLab