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 : MODULE m_genMTBasis
7 :
8 : CONTAINS
9 :
10 7248 : SUBROUTINE genMTBasis(atoms,enpara,vTot,fmpi,iType,jspin,usdus,f,g,flo,hub1data,l_writeArg)
11 : USE m_types
12 : USE m_constants
13 : USE m_radfun
14 : USE m_radflo
15 : USE m_find_enpara
16 : USE m_radsra
17 : USE m_intgr, ONLY : intgr0
18 : !$ use omp_lib
19 :
20 : IMPLICIT NONE
21 :
22 : TYPE(t_atoms), INTENT(IN) :: atoms
23 : TYPE(t_enpara), INTENT(IN) :: enpara
24 : TYPE(t_potden), INTENT(IN) :: vTot
25 : TYPE(t_mpi), INTENT(IN) :: fmpi
26 : TYPE(t_usdus), INTENT(INOUT) :: usdus
27 :
28 : INTEGER, INTENT(IN) :: iType
29 : INTEGER, INTENT(IN) :: jspin
30 :
31 : TYPE(t_hub1data), OPTIONAL, INTENT(IN) :: hub1data
32 :
33 : REAL, INTENT(INOUT) :: f(atoms%jmtd,2,0:atoms%lmaxd)
34 : REAL, INTENT(INOUT) :: g(atoms%jmtd,2,0:atoms%lmaxd)
35 : REAL, INTENT(INOUT) :: flo(atoms%jmtd,2,atoms%nlod)
36 : LOGICAL, OPTIONAL, INTENT(IN) :: l_writeArg
37 :
38 :
39 : INTEGER :: l, nodeu, noded, iLO, nLOsForL, iLOForL, nMTBasisFcts, iFunA, iFunB, iFunC
40 : REAL :: wronk, overlap, energy, e_lo, e_up, us, dus, coreNorm
41 :
42 :
43 : LOGICAL :: l_write,l_hia,l_performSpinavg
44 7248 : REAL :: vrTmp(atoms%jmtd)
45 : INTEGER :: i, info, addFunA, addFunB, iState, nQuantumNumber
46 :
47 7248 : REAL, ALLOCATABLE :: coreFun(:,:)
48 7248 : REAL, ALLOCATABLE :: funValues(:,:,:)
49 7248 : REAL, ALLOCATABLE :: bestRadFun(:,:)
50 7248 : REAL, ALLOCATABLE :: mtOverlapMat(:,:), invMTOverlapMat(:,:)
51 7248 : REAL, ALLOCATABLE :: rWorkArray(:), projVector(:), productVals(:)
52 7248 : INTEGER, ALLOCATABLE :: ipiv(:)
53 :
54 : INTEGER :: largestCoreMainQuantumNumbers(0:3)
55 :
56 7248 : l_write = .TRUE.
57 7248 : IF(PRESENT(l_writeArg)) l_write = l_writeArg
58 :
59 7248 : l_performSpinavg = .false.
60 7248 : IF(PRESENT(hub1data)) l_performSpinavg = hub1data%l_performSpinavg
61 :
62 7248 : l_write=l_write .and. fmpi%irank==0
63 2660 : !$ l_write = l_write .and. omp_get_num_threads()==1
64 :
65 :
66 1523 : IF (l_write) WRITE (oUnit,FMT=8000) iType
67 :
68 75176 : DO l = 0,atoms%lmax(iType)
69 : !Check if the orbital is to be treated with Hubbard 1
70 67928 : l_hia=.FALSE.
71 67928 : DO i = atoms%n_u+1, atoms%n_u+atoms%n_hia
72 67928 : IF(atoms%lda_u(i)%atomType.EQ.itype.AND.atoms%lda_u(i)%l.EQ.l) THEN
73 0 : l_hia=.TRUE.
74 : ENDIF
75 : ENDDO
76 : !In the case of a spin-polarized calculation with Hubbard 1 we want to treat
77 : !the correlated orbitals with a non-spin-polarized basis
78 67928 : IF((l_hia.AND.SIZE(vTot%mt,4).GT.1 .AND. l_performSpinavg).or.atoms%l_nonpolbas(itype)) THEN
79 1528128 : vrTmp = (vTot%mt(:,0,iType,1) + vTot%mt(:,0,iType,2))/2.0
80 : ELSE
81 49564608 : vrTmp = vTot%mt(:,0,iType,jspin)
82 : ENDIF
83 :
84 : CALL radfun(l,iType,jspin,enpara%el0(l,iType,jspin),vrTmp,atoms,&
85 67928 : f(1,1,l),g(1,1,l),usdus,nodeu,noded,wronk)
86 75176 : IF (l_write) THEN
87 14059 : WRITE (oUnit,FMT=8010) l,enpara%el0(l,iType,jspin),usdus%us(l,iType,jspin),usdus%dus(l,iType,jspin),&
88 28118 : nodeu,usdus%uds(l,iType,jspin),usdus%duds(l,iType,jspin),noded,usdus%ddn(l,iType,jspin),wronk
89 : END IF
90 : END DO
91 :
92 : ! Generate the extra wavefunctions for the local orbitals, if there are any.
93 7248 : IF (atoms%nlo(iType).GE.1) THEN
94 : CALL radflo(atoms,iType,jspin,enpara%ello0(1,1,jspin),vrtmp,f,g,fmpi,&
95 4888 : usdus,usdus%uuilon(1,1,jspin),usdus%duilon(1,1,jspin),usdus%ulouilopn(1,1,1,jspin),flo)
96 : END IF
97 :
98 36240 : largestCoreMainQuantumNumbers(:) = -1
99 45428 : DO iState = 1, atoms%econf(iType)%num_core_states
100 38180 : nQuantumNumber = atoms%econf(iType)%nprnc(iState)
101 38180 : l = atoms%econf(iType)%get_state_l(iState)
102 45428 : IF (nQuantumNumber.GT.largestCoreMainQuantumNumbers(l)) largestCoreMainQuantumNumbers(l) = nQuantumNumber
103 : END DO
104 :
105 : ! As quality measure: Calculate and write out overlap in MT spheres for each l for which LOs exist
106 7248 : IF (l_write) THEN
107 :
108 1523 : WRITE(oUnit,*) ''
109 1523 : WRITE(oUnit,'(a,i5,a)') 'Representation of radial basis functions, core electron states in terms of the other radial basis functions for atom type ',iType,':'
110 1523 : WRITE(oUnit,*) ''
111 :
112 15582 : DO l = 0, atoms%lmax(iType)
113 14059 : nLOsForL = 0
114 27164 : DO iLO = 1, atoms%nlo(iType)
115 27164 : If (atoms%llo(iLO, iType).EQ.l) nLOsForL = nLOsForL + 1
116 : END DO
117 14059 : nMTBasisFcts = nLOsForL + 2
118 : ! Fill radial function array for each function
119 56236 : ALLOCATE (funValues(atoms%jri(iType),2,nMTBasisFcts))
120 56236 : ALLOCATE (mtOverlapMat(nMTBasisFcts,nMTBasisFcts))
121 56236 : ALLOCATE (invMTOverlapMat(nMTBasisFcts-1,nMTBasisFcts-1))
122 42177 : ALLOCATE (rWorkArray(4*nMTBasisFcts))
123 42177 : ALLOCATE (ipiv(nMTBasisFcts))
124 42177 : ALLOCATE (bestRadFun(atoms%jri(iType),2))
125 28118 : ALLOCATE (coreFun(atoms%jri(iType),2))
126 42177 : ALLOCATE (projVector(nMTBasisFcts-1))
127 42177 : ALLOCATE (productVals(atoms%jri(iType)))
128 :
129 19490083 : funValues(:,:,1) = f(1:atoms%jri(iType),:,l)
130 19490083 : funValues(:,:,2) = g(1:atoms%jri(iType),:,l)
131 : iLOForL = 0
132 27164 : DO iLO = 1, atoms%nlo(iType)
133 27164 : If (atoms%llo(iLO, iType).EQ.l) THEN
134 1387 : iLOForL = iLOForL + 1
135 2133631 : funValues(:,:,iLOForL+2) = flo(1:atoms%jri(iType),:,iLO)
136 : END IF
137 : END DO
138 : ! calculate MT overlap
139 43564 : DO iFunA = 1, nMTBasisFcts
140 106743 : DO iFunB = 1, nMTBasisFcts
141 44288722 : productVals = 0.0
142 44288722 : DO i = 1, atoms%jri(iType)
143 44288722 : productVals(i) = funValues(i,1,iFunA)*funValues(i,1,iFunB)+funValues(i,2,iFunA)*funValues(i,2,iFunB)
144 : END DO
145 63179 : CALL intgr0(productVals,atoms%rmsh(1,iType),atoms%dx(iType),atoms%jri(iType),overlap)
146 92684 : mtOverlapMat(iFunA, iFunB) = overlap
147 : END DO
148 : END DO
149 :
150 14059 : IF(nLOsForL.GT.0) THEN
151 5536 : DO iFunC = 1, nMTBasisFcts
152 16628 : DO iFunA = 1, nMTBasisFcts
153 12475 : addFunA = 0
154 12475 : IF(iFunA.EQ.iFunC) CYCLE
155 8322 : IF(iFunA.GT.iFunC) addFunA = -1
156 8322 : projVector(iFunA+addFunA) = mtOverlapMat(iFunA,iFunC)
157 37489 : DO iFunB = 1, nMTBasisFcts
158 25014 : addFunB = 0
159 25014 : IF(iFunB.EQ.iFunC) CYCLE
160 16692 : IF(iFunB.GT.iFunC) addFunB = -1
161 37489 : invMTOverlapMat(iFunA+addFunA,iFunB+addFunB) = mtOverlapMat(iFunA,iFunB)
162 : END DO
163 : END DO
164 16628 : ipiv = 0
165 4153 : CALL DGETRF(nMTBasisFcts-1,nMTBasisFcts-1,invMTOverlapMat,nMTBasisFcts-1,ipiv(1:nMTBasisFcts-1),info)
166 4153 : CALL DGETRI(nMTBasisFcts-1,invMTOverlapMat,nMTBasisFcts-1,ipiv(1:nMTBasisFcts-1),rWorkArray,SIZE(rWorkArray),info)
167 45811 : projVector(:) = MATMUL(invMTOverlapMat(:,:),projVector(:))
168 6388757 : bestRadFun = 0.0
169 16628 : DO iFunA = 1, nMTBasisFcts
170 12475 : addFunA = 0
171 12475 : IF(iFunA.EQ.iFunC) CYCLE
172 8322 : IF(iFunA.GT.iFunC) addFunA = -1
173 : ! Calculate best representation of iFunC with all other basis functions.
174 12805939 : bestRadFun(:,:) = bestRadFun(:,:) + projVector(iFunA+addFunA)*funValues(:,:,iFunA)
175 : END DO
176 : ! Subtract iFunC from its best representation
177 6388757 : bestRadFun(:,:) = bestRadFun(:,:) - funValues(:,:,iFunC)
178 : ! Calculate norm of resulting function as representation error
179 3192302 : productVals(:) = bestRadFun(:,1)*bestRadFun(:,1) + bestRadFun(:,2)*bestRadFun(:,2)
180 : overlap = 0.0
181 4153 : CALL intgr0(productVals,atoms%rmsh(1,iType),atoms%dx(iType),atoms%jri(iType),overlap)
182 4153 : overlap = SQRT(overlap / mtOverlapMat(iFunC, iFunC))
183 5536 : WRITE(oUnit,'(a,i2,a,i2,a,f12.8)') 'Representation error for l = ', l, ', radial basis function ', iFunC, ': ', overlap
184 : END DO
185 : END IF
186 :
187 14059 : DEALLOCATE (projVector)
188 42177 : ALLOCATE(projVector(nMTBasisFcts))
189 :
190 14059 : IF (l.LE.3) THEN
191 6092 : IF (largestCoreMainQuantumNumbers(l).GT.0) THEN
192 4212526 : coreFun = 0.0
193 2994 : energy = find_enpara(.TRUE.,l,iType,jspin,largestCoreMainQuantumNumbers(l),atoms,vrTmp,e_lo,e_up,.TRUE.)
194 2994 : CALL radsra(energy,l,vrTmp(:),atoms%rmsh(1,iType),atoms%dx(iType),atoms%jri(iType),c_light(1.0), us,dus,nodeu,coreFun(:,1),coreFun(:,2))
195 10279 : DO iFunA = 1, nMTBasisFcts
196 5208890 : productVals = 0.0
197 5208890 : DO i = 1, atoms%jri(iType)
198 5208890 : productVals(i) = funValues(i,1,iFunA)*coreFun(i,1)+funValues(i,2,iFunA)*coreFun(i,2)
199 : END DO
200 7285 : CALL intgr0(productVals,atoms%rmsh(1,iType),atoms%dx(iType),atoms%jri(iType),overlap)
201 10279 : projVector(iFunA) = overlap
202 : END DO
203 10279 : ipiv = 0
204 2994 : CALL DGETRF(nMTBasisFcts,nMTBasisFcts,mtOverlapMat,nMTBasisFcts,ipiv,info)
205 2994 : CALL DGETRI(nMTBasisFcts,mtOverlapMat,nMTBasisFcts,ipiv,rWorkArray,SIZE(rWorkArray),info)
206 43318 : projVector(:) = MATMUL(mtOverlapMat(:,:),projVector(:))
207 4212526 : bestRadFun = 0.0
208 10279 : DO iFunA = 1, nMTBasisFcts
209 10428059 : bestRadFun(:,:) = bestRadFun(:,:) + projVector(iFunA)*funValues(:,:,iFunA)
210 : END DO
211 4212526 : bestRadFun(:,:) = bestRadFun(:,:) - coreFun(:,:)
212 2104766 : productVals(:) = bestRadFun(:,1)*bestRadFun(:,1) + bestRadFun(:,2)*bestRadFun(:,2)
213 : overlap = 0.0
214 2994 : CALL intgr0(productVals,atoms%rmsh(1,iType),atoms%dx(iType),atoms%jri(iType),overlap)
215 2104766 : productVals(:) = coreFun(:,1)*coreFun(:,1) + coreFun(:,2)*coreFun(:,2)
216 2994 : CALL intgr0(productVals,atoms%rmsh(1,iType),atoms%dx(iType),atoms%jri(iType),coreNorm)
217 2994 : overlap = SQRT(overlap / coreNorm)
218 5988 : WRITE(oUnit,'(a,i2,a,i2,a,f12.8)') 'Valence basis representation error for core state with n = ', largestCoreMainQuantumNumbers(l), ', l = ', l, ': ', overlap
219 : END IF
220 : END IF
221 :
222 14059 : DEALLOCATE (productVals)
223 14059 : DEALLOCATE (projVector)
224 14059 : DEALLOCATE (coreFun)
225 14059 : DEALLOCATE (bestRadFun)
226 14059 : DEALLOCATE (ipiv)
227 14059 : DEALLOCATE (rWorkArray)
228 14059 : DEALLOCATE (invMTOverlapMat)
229 14059 : DEALLOCATE (mtOverlapMat)
230 15582 : DEALLOCATE (funValues)
231 : END DO
232 : END IF
233 :
234 : 8000 FORMAT (1x,/,/,' wavefunction parameters for atom type',i5,':',&
235 : /,t32,'radial function',t79,'energy derivative',/,t3,&
236 : 'l',t8,'energy',t26,'value',t39,'derivative',t53,&
237 : 'nodes',t68,'value',t81,'derivative',t95,'nodes',t107,&
238 : 'norm',t119,'wronskian')
239 : 8010 FORMAT (i3,f10.5,2 (5x,1p,2e16.7,i5),1p,2e16.7)
240 :
241 7248 : END SUBROUTINE genMTBasis
242 :
243 7147 : END MODULE m_genMTBasis
|
