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_slomat
8 : !***********************************************************************
9 : ! updates the overlap matrix with the contributions from the local
10 : ! orbitals.
11 : ! p.kurz sept. 1996
12 : !***********************************************************************
13 : CONTAINS
14 11286 : SUBROUTINE slomat(input,atoms,sym,fmpi,lapw,cell,nococonv,ntyp,na,&
15 11286 : isp,ud, alo1,blo1,clo1,fjgj,&
16 : igSpinPr,igSpin,chi,smat,l_fullj,lapwq,fjgjq)
17 : !***********************************************************************
18 : ! locol stores the number of columns already processed; on parallel
19 : ! computers this decides, whether the LO-contribution is
20 : ! done on this node gb00
21 : !
22 : ! function legpol() at end of module
23 : !***********************************************************************
24 :
25 : USE m_constants,ONLY: fpi_const
26 : !USE m_types_mpimat
27 : USE m_types
28 : USE m_hsmt_fjgj
29 :
30 : IMPLICIT NONE
31 :
32 : TYPE(t_input),INTENT(IN) :: input
33 : TYPE(t_atoms),INTENT(IN) :: atoms
34 : TYPE(t_sym),INTENT(IN) :: sym
35 : TYPE(t_lapw),INTENT(IN),Target :: lapw
36 : TYPE(t_mpi),INTENT(IN) :: fmpi
37 : TYPE(t_cell),INTENT(IN) :: cell
38 : TYPE(t_nococonv),INTENT(IN) :: nococonv
39 : TYPE(t_fjgj),INTENT(IN),target :: fjgj
40 :
41 : ! Scalar Arguments
42 : INTEGER, INTENT (IN) :: na,ntyp
43 : INTEGER, INTENT (IN) :: igSpin,igSpinPr
44 : COMPLEX, INTENT (IN) :: chi
45 : INTEGER, INTENT(IN) :: isp
46 :
47 : ! Array Arguments
48 : REAL, INTENT (IN) :: alo1(atoms%nlod),blo1(atoms%nlod),clo1(atoms%nlod)
49 : TYPE(t_usdus),INTENT(IN) :: ud
50 : CLASS(t_mat),INTENT(INOUT) :: smat
51 : LOGICAL, INTENT(IN) :: l_fullj
52 :
53 : TYPE(t_lapw), OPTIONAL, TARGET, INTENT(IN) :: lapwq
54 : TYPE(t_fjgj), OPTIONAL, Target, INTENT(IN) :: fjgjq
55 :
56 : ! Local Scalars
57 : REAL :: con,dotp,fact1,fact2,fact3,fl2p1
58 : INTEGER :: invsfct,k ,l,lo,lop,lp,nkvec,nkvecp,kp,i
59 : INTEGER :: locol,lorow
60 : LOGICAL :: l_samelapw
61 :
62 11286 : COMPLEX, ALLOCATABLE :: cphPr(:), cph(:)
63 :
64 : TYPE(t_lapw),pointer :: lapwPr
65 : TYPE(t_fjgj),pointer :: fjgjPr
66 :
67 11286 : l_samelapw = .FALSE.
68 11286 : IF (.NOT.PRESENT(lapwq)) l_samelapw = .TRUE.
69 : IF (.NOT.l_samelapw) THEN
70 : lapwPr => lapwq
71 : fjgjPr => fjgjq
72 : ELSE
73 : lapwPr => lapw
74 : fjgjPr => fjgj
75 : END IF
76 :
77 56430 : ALLOCATE(cph(MAXVAL(lapw%nv)))
78 56430 : ALLOCATE(cphPr(MAXVAL(lapwPr%nv)))
79 :
80 : ! TODO: Introduce the logic for different lapw and the full rectangular
81 : ! instead of triangular construction...
82 :
83 11286 : CALL lapwPr%phase_factors(igSpinPr,atoms%taual(:,na),nococonv%qss,cphPr)
84 11286 : IF (l_samelapw.and.igspin==igSpinPr) THEN
85 1444464 : cph = cphPr
86 : ELSE
87 0 : CALL lapw%phase_factors(igSpin,atoms%taual(:,na),nococonv%qss,cph)
88 : END IF
89 :
90 11286 : IF ((sym%invsat(na) == 0) .OR. (sym%invsat(na) == 1)) THEN
91 : ! If this atom is the first of two atoms related by inversion, the
92 : ! contributions to the overlap matrix of both atoms are added simultaneously.
93 : ! Where it is made use of the fact, that the sum of these contributions
94 : ! is twice the real part of the contribution of each atom. Note, that in
95 : ! this case there are twice as many (2*(2*l+1)) k-vectors.
96 : ! (compare abccoflo and comments there).
97 : IF (sym%invsat(na) == 0) invsfct = 1
98 11286 : IF (sym%invsat(na) == 1) invsfct = 2
99 :
100 11286 : con = fpi_const/SQRT(cell%omtil)* ((atoms%rmt(ntyp))**2)/2.0
101 :
102 : !$acc kernels present(smat,smat%data_c,smat%data_r)&
103 : !$acc & copyin(fjgj,fjgj%fj,fjgj%gj,fjgjPr,fjgjPr%fj,fjgjPr%gj,l,lapw,lapw%kvec(:,:,na),lapwPr,lapwPr%kvec(:,:,na),ud,clo1(:),dotp,cph(:),cphPr(:),atoms,lapw%index_lo(:,na),lapw%gk(:,:,:)) &
104 : !$acc & copyin(lapwPr%index_lo(:,na),lapwPr%gk(:,:,:),ud%dulon(:,ntyp,isp),ud%ddn(:,ntyp,isp),ud%uulon(:,ntyp,isp),ud%uloulopn(:,:,ntyp,isp),blo1(:)) &
105 : !$acc & copyin(atoms,atoms%nlo(ntyp),lapw%nv(:),lapwPr%nv(:),atoms%llo(:,ntyp),alo1(:), fmpi, fmpi%n_size, fmpi%n_rank)&
106 : !$acc & default(none)
107 :
108 : !$acc loop gang private(fact1,fl2p1,l) independent
109 28724 : DO lo = 1,atoms%nlo(ntyp) !loop over all LOs for this atom
110 17438 : l = atoms%llo(lo,ntyp)
111 17438 : fl2p1 = (2*l+1)/fpi_const
112 : fact1 = (con**2) * fl2p1 * ( &
113 : & alo1(lo)*(alo1(lo) &
114 : & + 2*clo1(lo)*ud%uulon(lo,ntyp,isp)) + &
115 : & blo1(lo)*(blo1(lo)*ud%ddn(l,ntyp,isp) &
116 : & + 2*clo1(lo)*ud%dulon(lo,ntyp,isp)) + &
117 17438 : & clo1(lo)*clo1(lo) )
118 :
119 70144 : DO nkvec = 1,invsfct* (2*l+1) !Each LO can have several functions
120 : ! lop-kG part (l_fullj). |====|
121 41420 : IF (l_fullj) THEN
122 0 : lorow = lapwPr%nv(igSpinPr)+lapwPr%index_lo(lo,na)+nkvec
123 0 : kp = lapwPr%kvec(nkvec,lo,na)
124 : !$acc loop vector private(fact2,dotp) independent
125 0 : DO k = fmpi%n_rank + 1, lapw%nv(igSpin), fmpi%n_size
126 : fact2 = con * fl2p1 * ( &
127 : & fjgj%fj(k,l,isp,igSpin)*(alo1(lo) &
128 : & + clo1(lo)*ud%uulon(lo,ntyp,isp)) + &
129 : & fjgj%gj(k,l,isp,igSpin)*(blo1(lo)*ud%ddn(l,ntyp,isp) &
130 0 : & + clo1(lo)*ud%dulon(lo,ntyp,isp)) )
131 0 : dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
132 :
133 0 : IF (smat%l_real) THEN
134 : smat%data_r(lorow,(k-1) / fmpi%n_size + 1) = smat%data_r(lorow,(k-1) / fmpi%n_size + 1) &
135 : & + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
136 0 : & * CONJG(cphPr(kp)) * cph(k)
137 : ELSE
138 : smat%data_c(lorow,(k-1) / fmpi%n_size + 1) = smat%data_c(lorow,(k-1) / fmpi%n_size + 1) &
139 : & + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
140 0 : & * CONJG(cphPr(kp)) * cph(k)
141 : END IF
142 : END DO
143 : !$acc end loop
144 : END IF
145 41420 : locol = lapw%nv(igSpin)+lapw%index_lo(lo,na)+nkvec !this is the column of the matrix
146 58858 : IF (MOD(locol-1,fmpi%n_size) == fmpi%n_rank) THEN
147 26146 : locol=(locol-1)/fmpi%n_size+1 !this is the column in local storage!
148 26146 : k = lapw%kvec(nkvec,lo,na)
149 : ! Calculate the overlap matrix elements with the regular Flapw
150 : ! basis functions
151 : !$acc loop vector private(fact2,dotp,kp) independent
152 3109227 : DO kp = 1,lapwPr%nv(igSpinPr)
153 : fact2 = con * fl2p1 * ( &
154 : & fjgjPr%fj(kp,l,isp,igSpinPr)*(alo1(lo) &
155 : & + clo1(lo)*ud%uulon(lo,ntyp,isp)) + &
156 : & fjgjPr%gj(kp,l,isp,igSpinPr)*(blo1(lo)*ud%ddn(l,ntyp,isp) &
157 3083081 : & + clo1(lo)*ud%dulon(lo,ntyp,isp)) )
158 12332324 : dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
159 :
160 3109227 : IF (smat%l_real) THEN
161 : smat%data_r(kp,locol) = smat%data_r(kp,locol) &
162 : & + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
163 1122558 : & * CONJG(cphPr(kp)) * cph(k)
164 : ELSE
165 : smat%data_c(kp,locol) = smat%data_c(kp,locol) &
166 : & + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
167 1960523 : & * CONJG(cphPr(kp)) * cph(k)
168 : END IF
169 : END DO
170 : !$acc end loop
171 : ! Calculate the overlap matrix elements with other local orbitals
172 : ! of the same atom, if they have the same l
173 : !$acc loop vector private(lp,fact3,nkvecp,kp,lorow,dotp) independent
174 39764 : DO lop = 1, MERGE(lo-1,atoms%nlo(ntyp),igSpinPr==igSpin.AND..NOT.l_fullj)
175 13618 : IF (lop==lo) CYCLE !Do later
176 13618 : lp = atoms%llo(lop,ntyp)
177 39764 : IF (l == lp) THEN
178 : fact3 = con**2 * fl2p1 * ( &
179 : & alo1(lop)*(alo1(lo) &
180 : & + clo1(lo)*ud%uulon(lo,ntyp,isp)) + &
181 : & blo1(lop)*(blo1(lo)*ud%ddn(l,ntyp,isp) &
182 : & + clo1(lo)*ud%dulon(lo,ntyp,isp)) + &
183 : & clo1(lop)*(alo1(lo)*ud%uulon(lop,ntyp,isp) &
184 : & + blo1(lo)*ud%dulon(lop,ntyp,isp) &
185 60 : & + clo1(lo)*ud%uloulopn(lop,lo,ntyp,isp)) )
186 240 : DO nkvecp = 1,invsfct* (2*lp+1)
187 180 : kp = lapwPr%kvec(nkvecp,lop,na)
188 180 : lorow=lapwPr%nv(igSpinPr)+lapwPr%index_lo(lop,na)+nkvecp
189 720 : dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
190 :
191 240 : IF (smat%l_real) THEN
192 : smat%data_r(lorow,locol) = smat%data_r(lorow,locol) &
193 : & + chi * invsfct * fact3 * legpol(atoms%llo(lo,ntyp),dotp) &
194 180 : & * CONJG(cphPr(kp)) * cph(k)
195 : ELSE
196 : smat%data_c(lorow,locol) = smat%data_c(lorow,locol) &
197 : & + chi * invsfct * fact3 * legpol(atoms%llo(lo,ntyp),dotp) &
198 0 : & * CONJG(cphPr(kp)) * cph(k)
199 : END IF
200 : END DO
201 : END IF
202 : END DO
203 : ! Calculate the overlap matrix elements of one local
204 : ! orbital with itself
205 26146 : lop=lo
206 : !$acc loop vector private(kp,lorow,dotp) independent
207 75926 : DO nkvecp = 1,MERGE(nkvec,invsfct* (2*l+1),igSpinPr==igSpin.AND..NOT.l_fullj)
208 49780 : kp = lapwPr%kvec(nkvecp,lo,na)
209 49780 : lorow = lapwPr%nv(igSpinPr)+lapwPr%index_lo(lo,na)+nkvecp
210 199120 : dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
211 :
212 75926 : IF (smat%l_real) THEN
213 : smat%data_r(lorow,locol) = smat%data_r(lorow,locol) &
214 : & + chi * invsfct * fact1 * legpol(l,dotp) &
215 18047 : & * CONJG(cphPr(kp)) * cph(k)
216 : ELSE
217 : smat%data_c(lorow,locol) = smat%data_c(lorow,locol) &
218 : & + chi * invsfct * fact1 * legpol(l,dotp) &
219 31733 : & * CONJG(cphPr(kp)) * cph(k)
220 : END IF
221 : END DO
222 : !$acc end loop
223 : END IF ! mod(locol-1,n_size) = nrank
224 : END DO
225 : !$acc end loop
226 : END DO
227 : !$acc end loop
228 : !$acc end kernels
229 : END IF
230 :
231 11286 : END SUBROUTINE slomat
232 :
233 3133041 : PURE REAL FUNCTION legpol(l,arg)
234 : !$acc routine seq
235 :
236 : IMPLICIT NONE
237 :
238 : REAL,INTENT(IN) :: arg
239 : INTEGER,INTENT(IN):: l
240 :
241 : INTEGER :: lp
242 :
243 3133041 : REAL :: plegend(0:l)
244 :
245 3133041 : plegend(0) = 1.0
246 3133041 : IF (l.GE.1) THEN
247 2654253 : plegend(1) = arg
248 2870268 : DO lp = 1,l - 1
249 2870268 : plegend(lp+1) = (lp+lp+1)*arg*plegend(lp)/ (lp+1) -lp*plegend(lp-1)/ (lp+1)
250 : END DO
251 : END IF
252 3133041 : legpol = plegend(l)
253 3133041 : END FUNCTION legpol
254 :
255 : END MODULE m_slomat
|