hsmt_sph.F90 29 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7
MODULE m_hsmt_sph
8 9
  USE m_juDFT
  IMPLICIT NONE
10
CONTAINS
Daniel Wortmann's avatar
Daniel Wortmann committed
11
  SUBROUTINE hsmt_sph(n,atoms,mpi,isp,input,noco,cell,iintsp,jintsp,chi,lapw,el,e_shift,usdus,fj,gj,smat,hmat)
12 13 14
    USE m_constants, ONLY : fpi_const,tpi_const
    USE m_types
    IMPLICIT NONE
15
    TYPE(t_input),INTENT(IN)      :: input
Daniel Wortmann's avatar
Daniel Wortmann committed
16
    TYPE(t_mpi),INTENT(IN)        :: mpi
17 18 19
    TYPE(t_noco),INTENT(IN)       :: noco
    TYPE(t_cell),INTENT(IN)       :: cell
    TYPE(t_atoms),INTENT(IN)      :: atoms
Daniel Wortmann's avatar
Daniel Wortmann committed
20 21 22
    TYPE(t_lapw),INTENT(IN)       :: lapw
    TYPE(t_usdus),INTENT(IN)      :: usdus
    CLASS(t_mat),INTENT(INOUT)     :: smat,hmat
23 24
    !     ..
    !     .. Scalar Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
25 26
    INTEGER, INTENT (IN) :: n,isp,iintsp,jintsp
    COMPLEX, INTENT(IN)  :: chi
27 28
    !     ..
    !     .. Array Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
29 30 31
    REAL,    INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntype,input%jspins)
    REAL,    INTENT (IN) :: e_shift(input%jspins)
    REAL,    INTENT (IN) :: fj(:,0:,:),gj(:,0:,:)
32 33
    !     ..
    !     .. Local Scalars ..
Daniel Wortmann's avatar
Daniel Wortmann committed
34 35
    REAL tnn(3), elall,fct,fjkiln,gjkiln,ddnln,ski(3)
    REAL apw_lo1,apw_lo2,apw1,w1
36

Daniel Wortmann's avatar
Daniel Wortmann committed
37 38
    COMPLEX capw1
    INTEGER kii,ki,kj,l,nn
39 40 41

    !     ..
    !     .. Local Arrays ..
Daniel Wortmann's avatar
Daniel Wortmann committed
42 43
    REAL fleg1(0:atoms%lmaxd),fleg2(0:atoms%lmaxd),fl2p1(0:atoms%lmaxd)     
    REAL fl2p1bt(0:atoms%lmaxd)
44 45
    REAL qssbti(3),qssbtj(3)
    REAL, ALLOCATABLE :: plegend(:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
46
    COMPLEX, ALLOCATABLE :: cph(:)
47 48
    LOGICAL apw(0:atoms%lmaxd)

Daniel Wortmann's avatar
Daniel Wortmann committed
49
    CALL timestart("spherical setup")
50

Daniel Wortmann's avatar
Daniel Wortmann committed
51
<<<<<<< HEAD
52 53 54 55 56
    ! for Spin-orbit...
    REAL, ALLOCATABLE :: dplegend(:,:)
    REAL, ALLOCATABLE :: cross_k(:,:)
    INTEGER :: j1,j2
    COMPLEX :: isigma_x(2,2),isigma_y(2,2),isigma_z(2,2)
57
    COMPLEX :: chi11so(2,2),chi21so(2,2),chi22so(2,2),angso(DIMENSION%nvd,2,2)
58

59 60
    REAL :: tmp,tmp1,tmp2,tmp3

61 62 63

    IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) ALLOCATE ( aahlp(hlpmsize),bbhlp(hlpmsize) )
    !     ..
64
    con1 = fpi_const/SQRT(cell%omtil)
Daniel Wortmann's avatar
Daniel Wortmann committed
65 66
=======
>>>>>>> hsmt_simple
67
    DO l = 0,atoms%lmaxd
68 69 70
       fleg1(l) = REAL(l+l+1)/REAL(l+1)
       fleg2(l) = REAL(l)/REAL(l+1)
       fl2p1(l) = REAL(l+l+1)/fpi_const
71 72
       fl2p1bt(l) = fl2p1(l)*0.5
    END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
73
<<<<<<< HEAD
74 75 76
    !---> calculate the overlap of the radial basis functions with different
    !---> spin directions.
    IF (noco%l_constr) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
77 78
       ALLOCATE(uun21(0:atoms%lmaxd,atoms%ntype),udn21(0:atoms%lmaxd,atoms%ntype),&
            dun21(0:atoms%lmaxd,atoms%ntype),ddn21(0:atoms%lmaxd,atoms%ntype) )
79 80 81 82 83 84 85 86 87 88
       CALL rad_ovlp(atoms,usdus,input,vr,el(0:,:,:), uun21,udn21,dun21,ddn21)
    ENDIF
    !---> loop over each atom type

    DO iintsp = 1,nintsp
       !$OMP parallel do DEFAULT(SHARED)&
       !$OMP PRIVATE(n,l,apw,lo,k,gs,fb,gb,ws,ff,gg)
       DO n = 1,atoms%ntype

          DO l = 0,atoms%lmax(n)
89
             apw(l) = .FALSE.
90
             DO lo = 1,atoms%nlo(n)
91
                IF (atoms%l_dulo(lo,n)) apw(l) = .TRUE.
92
             ENDDO
93
             IF ((input%l_useapw).AND.(atoms%lapw_l(n).GE.l)) apw(l) = .FALSE.
94

95 96
          ENDDO
          DO lo = 1,atoms%nlo(n)
97
             IF (atoms%l_dulo(lo,n)) apw(atoms%llo(lo,n)) = .TRUE.
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
          ENDDO

          DO k = 1,lapw%nv(iintsp)
             gs = lapw%rk(k,iintsp)*atoms%rmt(n)
             CALL sphbes(atoms%lmax(n),gs, fb)
             CALL dsphbs(atoms%lmax(n),gs,fb, gb)
             DO l = 0,atoms%lmax(n)
                !---> set up wronskians for the matching conditions for each ntype
                ws = con1/(usdus%uds(l,n,isp)*usdus%dus(l,n,isp)&
                     - usdus%us(l,n,isp)*usdus%duds(l,n,isp))
                ff = fb(l)
                gg = lapw%rk(k,iintsp)*gb(l)
                IF ( apw(l) ) THEN
                   fj(k,l,n,iintsp) = 1.0*con1 * ff / usdus%us(l,n,isp)
                   gj(k,l,n,iintsp) = 0.0d0
                ELSE
114
                   IF (noco%l_constr.OR.l_socfirst) THEN
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
                      !--->                 in a constrained calculation fj and gj are needed
                      !--->                 both local spin directions (isp) at the same time
                      fj(k,l,n,isp) = ws * ( usdus%uds(l,n,isp)*gg - usdus%duds(l,n,isp)*ff )
                      gj(k,l,n,isp) = ws * ( usdus%dus(l,n,isp)*ff - usdus%us(l,n,isp)*gg )
                   ELSE
                      !--->                 in a spin-spiral calculation fj and gj are needed
                      !--->                 both interstitial spin directions at the same time
                      !--->                 In all other cases iintsp runs from 1 to 1.
                      fj(k,l,n,iintsp) = ws * ( usdus%uds(l,n,isp)*gg - usdus%duds(l,n,isp)*ff )
                      gj(k,l,n,iintsp) = ws * ( usdus%dus(l,n,isp)*ff - usdus%us(l,n,isp)*gg )
                   ENDIF
                ENDIF
             ENDDO
          ENDDO ! k = 1, lapw%nv
       ENDDO    ! n = 1,atoms%ntype
       !$OMP end parallel do

    ENDDO       ! iintsp = 1,nintsp
    !
    IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
       !---> pk non-collinear
       !--->    initialize hamiltonian help array
137 138
       aahlp=CMPLX(0.0,0.0)
       bbhlp=CMPLX(0.0,0.0)
139 140 141 142 143 144 145 146
    ENDIF

    !$OMP PARALLEL  DEFAULT(shared)&
    !$OMP PRIVATE(kii,ki,nc,iii,kjmax,ski,kj,plegend,l,n1,n)&
    !$OMP PRIVATE(n0,rph,cph,nn,tnn,fjkiln,gjkiln)&
    !$OMP PRIVATE(w1,apw_lo1,apw_lo2,ddnln,elall,fct,ij,apw1)&
    !$OMP PRIVATE(cross_k,dplegend,chi,chi11,chi21,chi22,nsp,chj)&
    !$OMP PRIVATE(isigma_x,isigma_y,isigma_z,j1,j2,chi11so,chi21so,chi22so)&
Gregor Michalicek's avatar
Gregor Michalicek committed
147
    !$OMP PRIVATE(tmp1,tmp2,tmp3,tmp)&
148 149 150
    !$OMP PRIVATE(aawa,bbwa,capw1,ii) IF (.not.l_socfirst)
    !$     IF (l_socfirst) WRITE(*,*) "WARNING: first variation SOC does not work with OPENMP in hsmt_sph"
    !$     IF (l_socfirst) WRITE(*,*) "         switching off openmp parallelization"
151 152 153
    ALLOCATE(rph(DIMENSION%nvd,ab_dim))
    ALLOCATE(cph(DIMENSION%nvd,ab_dim))
    ALLOCATE(plegend(DIMENSION%nvd,0:atoms%lmaxd))
154
    IF (l_socfirst)THEN
155
       ALLOCATE ( dplegend(DIMENSION%nvd,0:atoms%lmaxd),cross_k(DIMENSION%nvd,3))
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178
       dplegend(:,0)=0.e0
       dplegend(:,1)=1.e0
    ENDIF

    plegend=0.0
    plegend(:,0)=1.0
    DO iintsp = 1,nintsp
       IF (iintsp.EQ.1) THEN
          qssbti = - noco%qss/2
       ELSE
          qssbti = + noco%qss/2
       ENDIF
       DO jintsp = 1,iintsp
          IF (jintsp.EQ.1) THEN
             qssbtj = - noco%qss/2
          ELSE
             qssbtj = + noco%qss/2
          ENDIF

          nc = 0                                    ! maybe IF (iintsp.EQ
          IF ( noco%l_noco .AND. (n_size.GT.1) ) THEN
             !--->       for EV-parallelization & noco ( see comments at top )
             lapw%nv_tot = lapw%nv(1) + lapw%nv(2)
179
             IF (noco%l_ss)  CALL juDFT_error("ev-parallelization & spin-spiral !",calledby ="hsmt_sph")
180 181 182 183
          ELSE
             lapw%nv_tot = lapw%nv(iintsp)
          ENDIF

Daniel Wortmann's avatar
Daniel Wortmann committed
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222
=======
    !$OMP PARALLEL DEFAULT(SHARED)&
    !$OMP PRIVATE(kii,ki,ski,kj,plegend,l)&
    !$OMP PRIVATE(cph,nn,tnn,fjkiln,gjkiln)&
    !$OMP PRIVATE(w1,apw_lo1,apw_lo2,ddnln,elall,fct,apw1)&
    !$OMP PRIVATE(capw1) 
    ALLOCATE(cph(MAXVAL(lapw%nv)))
    ALLOCATE(plegend(MAXVAL(lapw%nv),0:atoms%lmaxd))
    plegend=0.0
    plegend(:,0)=1.0
    qssbti=MERGE(- noco%qss/2,+ noco%qss/2,iintsp.EQ.1)
    qssbtj=MERGE(- noco%qss/2,+ noco%qss/2,jintsp.EQ.1)
    !$OMP  DO SCHEDULE(DYNAMIC,1)
    DO  ki =  mpi%n_rank+1, lapw%nv(iintsp), mpi%n_size
       kii=(ki-1)/mpi%n_size+1
       ski = lapw%gvec(:,ki,iintsp) + qssbti
       !--->       legendre polynomials
       DO kj = 1,ki
          plegend(kj,1) = DOT_PRODUCT(lapw%gk(:,kj,jintsp),lapw%gk(:,ki,iintsp))
       END DO
       DO l = 1,atoms%lmax(n) - 1
          plegend(:ki,l+1) = fleg1(l)*plegend(:ki,1)*plegend(:ki,l) - fleg2(l)*plegend(:ki,l-1)
       END DO
       !--->             set up phase factors
       cph = 0.0
       DO nn = SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n))
          tnn = tpi_const*atoms%taual(:,nn)
          DO kj = 1,ki
             cph(kj) = cph(kj) +&
                  CMPLX(COS(DOT_PRODUCT(ski-lapw%gvec(:,kj,jintsp)+qssbtj,tnn)),&
                  SIN(DOT_PRODUCT(lapw%gvec(:,kj,jintsp)+qssbtj-ski,tnn)))
          END DO
       END DO

       !--->          update overlap and l-diagonal hamiltonian matrix
       DO  l = 0,atoms%lmax(n)
          fjkiln = fj(ki,l,iintsp)
          gjkiln = gj(ki,l,iintsp)
>>>>>>> hsmt_simple
223
          !
Daniel Wortmann's avatar
Daniel Wortmann committed
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
          IF (input%l_useapw) THEN
             w1 = 0.5 * ( usdus%uds(l,n,isp)*usdus%dus(l,n,isp) + &
                  usdus%us(l,n,isp)*usdus%duds(l,n,isp) )
             apw_lo1 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( gjkiln * w1 +&
                  fjkiln * usdus%us(l,n,isp) * usdus%dus(l,n,isp) )
             apw_lo2 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( fjkiln * w1 +&
                  gjkiln * usdus%uds(l,n,isp) * usdus%duds(l,n,isp) )
             !
          ENDIF
          ddnln =  usdus%ddn(l,n,isp)
          elall = el(l,n,isp)
          IF (l<=atoms%lnonsph(n)) elall=elall-e_shift(isp)
          IF (smat%l_real) THEN
             DO kj = 1,ki
                fct  = plegend(kj,l)*fl2p1(l)*&
                     ( fjkiln*fj(kj,l,jintsp) + gjkiln*gj(kj,l,jintsp)*ddnln )
                smat%data_r(kj,kii)=smat%data_r(kj,kii)+REAL(cph(kj))*fct
                hmat%data_r(kj,kii)=hmat%data_r(kj,kii) + REAL(cph(kj)) * &
                     ( fct * elall + plegend(kj,l) * fl2p1bt(l) *&
                     ( fjkiln*gj(kj,l,jintsp) + gjkiln*fj(kj,l,jintsp) ) )
                !+APW
                IF (input%l_useapw) THEN
                   apw1 = REAL(cph(kj)) * plegend(kj,l)  * &
                        ( apw_lo1 * fj(kj,l,iintsp) + apw_lo2 * gj(kj,l,iintsp) )
                   hmat%data_r(kj,kii)=hmat%data_r(kj,kii) + apw1
249
                ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
250 251 252 253 254 255 256 257 258 259 260 261 262 263 264
                !-APW
             ENDDO
          ELSE
             DO kj = 1,ki
                fct  = chi*plegend(kj,l)*fl2p1(l)*&
                     ( fjkiln*fj(kj,l,jintsp) + gjkiln*gj(kj,l,jintsp)*ddnln )

                smat%data_c(kj,kii)=smat%data_c(kj,kii) + cph(kj)*fct
                hmat%data_c(kj,kii)=hmat%data_c(kj,kii) + cph(kj) * ( fct*elall &
                     + chi*plegend(kj,l)*fl2p1bt(l) * ( fjkiln*gj(kj,l,jintsp) + gjkiln*fj(kj,l,jintsp) ) )

                IF (input%l_useapw) THEN
                   capw1 = cph(kj)*plegend(kj,l)&
                        * ( apw_lo1 * fj(kj,l,iintsp) + apw_lo2 * gj(kj,l,iintsp) )
                   hmat%data_c(kj,kii)=hmat%data_c(kj,kii) + capw1
265 266
                ENDIF
             END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
267
<<<<<<< HEAD
268 269 270 271 272 273 274 275
             !--->       loop over equivalent atoms
             n1 = 0
             DO n = 1,atoms%ntype

                IF (noco%l_noco) THEN
                   !--->          pk non-collinear
                   !--->             set up the spinors of this atom within global
                   !--->             spin-coordinateframe
276
                   CALL hsmt_spinor(isp,n, noco, input,chi, chi11, chi21, chi22,l_socfirst,&
277 278 279 280 281 282 283 284 285 286
                        isigma,isigma_x,isigma_y,isigma_z,chi11so,chi21so,chi22so,angso,chj,cross_k)
                ENDIF
                !---> pk non-collinear
                n0 = n1 + 1
                n1 = n1 + atoms%neq(n)
                rph(:,1) = 0.0
                cph(:,1) = 0.0
                DO nn = n0,n1
                   tnn = tpi_const*atoms%taual(:,nn)
                   !--->             set up phase factors
287
                   !!$OMP SIMD PRIVATE(tmp1,tmp2,tmp3,tmp) 
288
                   DO kj = 1,kjmax
289 290 291 292 293 294 295
                      !rph(kj,1) = rph(kj,1) +&
                      !     COS(DOT_PRODUCT(ski-(/lapw%k1(kj,jintsp),lapw%k2(kj,jintsp),lapw%k3(kj,jintsp)/)+qssbtj,tnn))
                      tmp1 = (ski(1)-lapw%k1(kj,jintsp)+qssbtj(1)) * tnn(1)
                      tmp2 = (ski(2)-lapw%k2(kj,jintsp)+qssbtj(2)) * tnn(2)
                      tmp3 = (ski(3)-lapw%k3(kj,jintsp)+qssbtj(3)) * tnn(3)
                      tmp = cos(tmp1 + tmp2 + tmp3)
                      rph(kj,1) = rph(kj,1) + tmp
296
                   END DO
297 298 299
                   IF (.NOT.sym%invs) THEN
                      !--->                if the system does not posses inversion symmetry
                      !--->                the complex part of the exponential is needed.
300
                      !!$OMP SIMD PRIVATE(tmp1,tmp2,tmp3,tmp) 
301 302 303 304 305 306 307 308 309 310
                      DO kj = 1,kjmax
                      !   cph(kj,1) = cph(kj,1) +&
                      !        SIN(DOT_PRODUCT((/lapw%k1(kj,jintsp),lapw%k2(kj,jintsp),lapw%k3(kj,jintsp)/)+qssbtj-ski,tnn))
                      tmp1 = (lapw%k1(kj,jintsp)+qssbtj(1)-ski(1)) * tnn(1)
                      tmp2 = (lapw%k2(kj,jintsp)+qssbtj(2)-ski(2)) * tnn(2)
                      tmp3 = (lapw%k3(kj,jintsp)+qssbtj(3)-ski(3)) * tnn(3)
                      tmp = sin(tmp1 + tmp2 + tmp3)
                      cph(kj,1) = cph(kj,1) + tmp
                      END DO
                   ENDIF
311 312 313 314
                END DO

                !--->          update overlap and l-diagonal hamiltonian matrix
                DO  l = 0,atoms%lmax(n)
315
                   IF (noco%l_constr.OR.l_socfirst) THEN
316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334
                      fjkiln = fj(ki,l,n,isp)
                      gjkiln = gj(ki,l,n,isp)
                   ELSE
                      fjkiln = fj(ki,l,n,iintsp)
                      gjkiln = gj(ki,l,n,iintsp)
                   ENDIF
                   !
                   w1 = 0.5 * ( usdus%uds(l,n,isp)*usdus%dus(l,n,isp) + &
                        usdus%us(l,n,isp)*usdus%duds(l,n,isp) )
                   apw_lo1 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( gjkiln * w1 +&
                        fjkiln * usdus%us(l,n,isp) * usdus%dus(l,n,isp) )
                   apw_lo2 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( fjkiln * w1 +&
                        gjkiln * usdus%uds(l,n,isp) * usdus%duds(l,n,isp) )
                   !
                   ddnln =  usdus%ddn(l,n,isp)
                   elall = el(l,n,isp)

                   IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
                      !--->             pk non-collinear
335
                      IF (noco%l_constr.OR.l_socfirst) THEN
336 337 338 339
                         DO kj = 1,ki
                            fct  = plegend(kj,l)*fl2p1(l)*&
                                 ( fjkiln*fj(kj,l,n,isp) + gjkiln*gj(kj,l,n,isp)*ddnln )

340 341
                            bbwa(kj) = CMPLX(rph(kj,1),cph(kj,1))*fct
                            aawa(kj) = CMPLX(rph(kj,1),cph(kj,1)) * ( &
342 343
                                 fct*elall + plegend(kj,l)*fl2p1bt(l)*&
                                 ( fjkiln*gj(kj,l,n,isp) + gjkiln*fj(kj,l,n,isp) ) )
344 345 346 347 348
                            IF (input%l_useapw) THEN
                               capw1 = CMPLX(rph(kj,1),cph(kj,1))*plegend(kj,l)&
                                    * ( apw_lo1 * fj(kj,l,n,isp) + apw_lo2 * gj(kj,l,n,isp) )
                               aawa(kj) = aawa(kj) + capw1
                            ENDIF
349 350 351 352 353 354
                         ENDDO
                      ELSE
                         DO kj = 1,ki
                            fct  = plegend(kj,l)*fl2p1(l)*&
                                 ( fjkiln*fj(kj,l,n,jintsp) + gjkiln*gj(kj,l,n,jintsp)*ddnln )

355 356
                            bbwa(kj) = CMPLX(rph(kj,1),cph(kj,1))*fct
                            aawa(kj) = CMPLX(rph(kj,1),cph(kj,1)) * (&
357 358
                                 fct*elall + plegend(kj,l)*fl2p1bt(l)*&
                                 ( fjkiln*gj(kj,l,n,jintsp) + gjkiln*fj(kj,l,n,jintsp) ) )
359 360 361 362 363
                            IF (input%l_useapw) THEN
                               capw1 = CMPLX(rph(kj,1),cph(kj,1))*plegend(kj,l)&
                                    * ( apw_lo1 * fj(kj,l,n,jintsp) + apw_lo2 * gj(kj,l,n,jintsp) )
                               aawa(kj) = aawa(kj) + capw1
                            ENDIF
364 365 366 367 368
                         ENDDO
                      ENDIF
                      !+||
                      IF ( kii+1.LE.lapw%nv(1) ) THEN
                         !--->                spin-up spin-up part
369 370
                         CALL CPP_BLAS_caxpy(ki,chi11,bbwa,1,hamOvlp%b_c(iii+1),1)
                         CALL CPP_BLAS_caxpy(ki,chi11,aawa,1,hamOvlp%a_c(iii+1),1)
371 372 373 374 375 376
                         !--->                spin-down spin-up part, upper triangle.
                         !--->                the help array is used to allow vectorization on
                         !--->                Cray PVP systems. it is mapped onto the hamiltonian
                         !--->                matrix after the setup is complete.
                         DO kj = 1,ki - 1
                            ii = iii + kj
377 378
                            aahlp(ii)=aahlp(ii)+CONJG(aawa(kj))*chi21
                            bbhlp(ii)=bbhlp(ii)+CONJG(bbwa(kj))*chi21
379 380 381 382 383 384 385 386 387
                         END DO
                      ENDIF
                      IF ( (kii+1.GT.lapw%nv(1)).OR.(n_size.EQ.1) ) THEN
                         IF (n_size.EQ.1) THEN
                            ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2
                         ELSE
                            ii = iii
                         ENDIF
                         !--->                spin-down spin-up part, lower triangle
388 389
                         CALL CPP_BLAS_caxpy(ki,chi21,bbwa,1,hamOvlp%b_c(ii+1),1)
                         CALL CPP_BLAS_caxpy(ki,chi21,aawa,1,hamOvlp%a_c(ii+1),1)
390 391
                         !--->                spin-down spin-down part
                         ii = ii + lapw%nv(1)+atoms%nlotot
392 393
                         CALL CPP_BLAS_caxpy(ki,chi22,bbwa,1,hamOvlp%b_c(ii+1),1)
                         CALL CPP_BLAS_caxpy(ki,chi22,aawa,1,hamOvlp%a_c(ii+1),1)
394 395 396 397 398 399 400 401 402
                      ENDIF
                      !-||
                      !--->                when fj and gj are available for both local spins
                      !--->                (isp), add the contribution from the constraint
                      !--->                B-field.
                      IF (noco%l_constr .AND. (isp .EQ. 2)) THEN
                         DO nsp = 1,input%jspins
                            IF (nsp.EQ.1) THEN
                               DO kj = 1,lapw%nv(1)
403 404
                                  aawa(kj) = (-0.5)*CMPLX(noco%b_con(1,n),noco%b_con(2,n))*&
                                       CMPLX(rph(kj,1),cph(kj,1))*&
405 406 407 408 409 410 411 412
                                       plegend(kj,l)*fl2p1(l)*(&
                                       + fj(ki,l,n,2)*fj(kj,l,n,1)*uun21(l,n)&
                                       + fj(ki,l,n,2)*gj(kj,l,n,1)*udn21(l,n)&
                                       + gj(ki,l,n,2)*fj(kj,l,n,1)*dun21(l,n)&
                                       + gj(ki,l,n,2)*gj(kj,l,n,1)*ddn21(l,n))
                               ENDDO
                            ELSE
                               DO kj = 1,lapw%nv(1)
413 414
                                  aawa(kj) = (-0.5)*CMPLX(noco%b_con(1,n),-noco%b_con(2,n))*&
                                       CMPLX(rph(kj,1),cph(kj,1))*&
415 416 417 418 419 420 421 422 423 424
                                       plegend(kj,l)*fl2p1(l)*(&
                                       + fj(ki,l,n,1)*fj(kj,l,n,2)*uun21(l,n)&
                                       + fj(ki,l,n,1)*gj(kj,l,n,2)*dun21(l,n)&
                                       + gj(ki,l,n,1)*fj(kj,l,n,2)*udn21(l,n)&
                                       + gj(ki,l,n,1)*gj(kj,l,n,2)*ddn21(l,n))
                               ENDDO
                            ENDIF
                            !--->                  spin-up spin-up part
                            ii = (ki-1)*(ki)/2
                            DO kj = 1,ki
425
                               hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chj(nsp,1,1,n)
426 427 428 429 430
                            ENDDO
                            !--->                  spin-down spin-down part
                            ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2 + &
                                 lapw%nv(1)+atoms%nlotot
                            DO kj = 1,ki
431
                               hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chj(nsp,2,2,n)
432 433 434 435
                            ENDDO
                            !--->                  spin-down spin-up part
                            ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2
                            DO kj = 1,lapw%nv(1)
436
                               hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chj(nsp,2,1,n)
437 438 439 440 441 442 443 444 445 446 447
                            ENDDO
                         ENDDO
                      ENDIF

                      !                 Add spin-orbit coupling
                      IF (isp.EQ.2) THEN
                         IF ((l.GT.0).AND.l_socfirst) THEN

                            DO j1 = 1,input%jspins
                               DO j2 = 1,input%jspins
                                  DO kj = 1,lapw%nv(1)
448
                                     aawa(kj)=CMPLX(rph(kj,1),cph(kj,1))*(&
449 450 451 452 453 454 455 456 457 458 459
                                          + fj(ki,l,n,j1)*fj(kj,l,n,j2)*rsoc%rsopp(n,l,j1,j2)&
                                          + fj(ki,l,n,j1)*gj(kj,l,n,j2)*rsoc%rsopdp(n,l,j1,j2)&
                                          + gj(ki,l,n,j1)*fj(kj,l,n,j2)*rsoc%rsoppd(n,l,j1,j2)&
                                          + gj(ki,l,n,j1)*gj(kj,l,n,j2)*rsoc%rsopdpd(n,l,j1,j2))&
                                          *dplegend(kj,l)*fl2p1(l)*angso(kj,j1,j2)
                                  ENDDO
                                  IF (n_size.EQ.1) THEN

                                     !--->                  spin-up spin-up part
                                     ii = (ki-1)*(ki)/2
                                     DO kj = 1,ki
460
                                        hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chi11so(j1,j2)
461 462 463 464 465
                                     ENDDO
                                     !--->                  spin-down spin-down part
                                     ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2 +&
                                          lapw%nv(1)+atoms%nlotot
                                     DO kj = 1,ki
466
                                        hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chi22so(j1,j2)
467 468 469 470
                                     ENDDO
                                     !--->                  spin-down spin-up part
                                     ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2
                                     DO kj = 1,lapw%nv(1)
471
                                        hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chi21so(j1,j2)
472 473 474 475 476 477
                                     ENDDO

                                  ELSE  ! eigenvalue parallelization

                                     IF ( kii+1.LE.lapw%nv(1) ) THEN
                                        !--->                    spin-up spin-up part
478
                                        CALL CPP_BLAS_caxpy(ki,chi11so(j1,j2),aawa,1, hamOvlp%a_c(iii+1),1)
479 480 481 482

                                        !--->                    spin-down spin-up part, upper triangle.
                                        DO kj = 1,ki - 1
                                           ii = iii + kj
483
                                           aahlp(ii) = aahlp(ii) + CONJG(aawa(kj))*chi21so(j2,j1)
484 485 486 487 488
                                        END DO
                                     ENDIF
                                     IF  (kii+1.GT.lapw%nv(1)) THEN
                                        ii = iii
                                        !--->                    spin-down spin-up part, lower triangle
489
                                        CALL CPP_BLAS_caxpy(ki,chi21so(j1,j2),aawa,1, hamOvlp%a_c(ii+1),1)
490 491
                                        !--->                    spin-down spin-down part
                                        ii = ii + lapw%nv(1)+atoms%nlotot
492
                                        CALL CPP_BLAS_caxpy(ki,chi22so(j1,j2),aawa,1, hamOvlp%a_c(ii+1),1)
493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511
                                     ENDIF

                                  ENDIF ! eigenvalue par.

                               ENDDO ! j2
                            ENDDO    ! j1
                         ENDIF       ! ( l > 0 ) & socfirst
                      ENDIF          ! (isp = 2)
                      !               End spin-orbit
                   ELSEIF ( noco%l_noco .AND. noco%l_ss ) THEN
                      IF ( iintsp.EQ.2 .AND. jintsp.EQ.1 ) THEN
                         kjmax = lapw%nv(1)
                      ELSE
                         kjmax = ki
                      ENDIF
                      DO kj = 1,kjmax
                         fct  = plegend(kj,l)*fl2p1(l)* ( fjkiln*fj(kj,l,n,jintsp) +&
                              gjkiln*gj(kj,l,n,jintsp)*ddnln )

512 513
                         bbwa(kj) = CMPLX(rph(kj,1),cph(kj,1))*fct
                         aawa(kj) = CMPLX(rph(kj,1),cph(kj,1)) * ( &
514 515 516 517 518 519 520
                              fct*elall + plegend(kj,l)*fl2p1bt(l)*&
                              ( fjkiln*gj(kj,l,n,jintsp) + gjkiln*fj(kj,l,n,jintsp) ) )
                      ENDDO
                      IF ( iintsp.EQ.1 .AND. jintsp.EQ.1 ) THEN
                         !--->                   spin-up spin-up part
                         ii = (ki-1)*(ki)/2
                         DO kj = 1,ki
521 522
                            hamOvlp%b_c(ii+kj) = hamOvlp%b_c(ii+kj) + bbwa(kj)*chi11
                            hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chi11
523 524 525 526 527 528
                         ENDDO
                      ELSEIF ( iintsp.EQ.2 .AND. jintsp.EQ.2 ) THEN
                         !--->                   spin-down spin-down part
                         ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2 +&
                              lapw%nv(1)+atoms%nlotot
                         DO kj = 1,ki
529 530
                            hamOvlp%b_c(ii+kj) = hamOvlp%b_c(ii+kj) + bbwa(kj)*chi22
                            hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chi22
531 532 533 534 535
                         ENDDO
                      ELSE
                         !--->                   spin-down spin-up part
                         ii = (lapw%nv(1)+atoms%nlotot+ki-1)*(lapw%nv(1)+atoms%nlotot+ki)/2
                         DO kj = 1,lapw%nv(1)
536 537
                            hamOvlp%b_c(ii+kj) = hamOvlp%b_c(ii+kj) + bbwa(kj)*chi21
                            hamOvlp%a_c(ii+kj) = hamOvlp%a_c(ii+kj) + aawa(kj)*chi21
538 539 540 541
                         ENDDO
                      ENDIF
                      !--->             pk non-collinear
                   ELSE
542 543 544 545
                      IF (l_real) THEN
                         DO kj = 1,ki
                            fct  = plegend(kj,l)*fl2p1(l)*&
                                 ( fjkiln*fj(kj,l,n,jintsp) + gjkiln*gj(kj,l,n,jintsp)*ddnln )
546

547
                            ij = iii + kj
548 549
                            hamOvlp%b_r(ij) = hamOvlp%b_r(ij) + rph(kj,1) * fct
                            hamOvlp%a_r(ij) = hamOvlp%a_r(ij) + rph(kj,1) * ( fct * elall + plegend(kj,l) * fl2p1bt(l) *&
550 551 552 553 554
                                 ( fjkiln*gj(kj,l,n,jintsp) + gjkiln*fj(kj,l,n,jintsp) ) )
                            !+APW
                            IF (input%l_useapw) THEN
                               apw1 = rph(kj,1) * plegend(kj,l)  * &
                                    ( apw_lo1 * fj(kj,l,n,iintsp) + apw_lo2 * gj(kj,l,n,iintsp) )
555
                               hamOvlp%a_r(ij) = hamOvlp%a_r(ij) + apw1
556 557 558 559 560 561 562 563 564
                            ENDIF
                            !-APW
                         ENDDO
                      ELSE
                         DO kj = 1,ki
                            fct  = plegend(kj,l)*fl2p1(l)*&
                                 ( fjkiln*fj(kj,l,n,jintsp) + gjkiln*gj(kj,l,n,jintsp)*ddnln )

                            ij = iii + kj
565 566
                            hamOvlp%b_c(ij) = hamOvlp%b_c(ij) + CMPLX(rph(kj,1),cph(kj,1))*fct
                            hamOvlp%a_c(ij) = hamOvlp%a_c(ij) + CMPLX(rph(kj,1),cph(kj,1)) * (fct*elall + plegend(kj,l)*fl2p1bt(l) *&
567 568 569 570 571
                                 ( fjkiln*gj(kj,l,n,jintsp) + gjkiln*fj(kj,l,n,jintsp) ) )
                            IF (input%l_useapw) THEN

                               capw1 = CMPLX(rph(kj,1),cph(kj,1))*plegend(kj,l)&
                                    * ( apw_lo1 * fj(kj,l,n,iintsp) + apw_lo2 * gj(kj,l,n,iintsp) )
572
                               hamOvlp%a_c(ij) = hamOvlp%a_c(ij) + capw1
573 574 575
                            ENDIF
                         END DO
                      ENDIF
576 577 578
                   ENDIF

                   !--->          end loop over l
579
                ENDDO
580 581

                !--->       end loop over atom types (ntype)
582
             ENDDO
583 584 585

             !--->    end loop over ki

586
          ENDDO
587 588 589 590
          !$OMP END DO
          !---> end loops over interstitial spins
       ENDDO ! jintsp = 1,iintsp
    ENDDO   ! iintsp = 1,nintsp
Daniel Wortmann's avatar
Daniel Wortmann committed
591 592 593 594 595 596 597 598 599
=======
          ENDIF
          !--->          end loop over l
       ENDDO
       !--->    end loop over ki
    ENDDO
    !$OMP END DO
    !--->       end loop over atom types (ntype)
>>>>>>> hsmt_simple
600
    DEALLOCATE(plegend)
Daniel Wortmann's avatar
Daniel Wortmann committed
601
    DEALLOCATE(cph)
602
    !$OMP END PARALLEL
Daniel Wortmann's avatar
Daniel Wortmann committed
603
    CALL timestop("spherical setup")
604 605 606 607

    RETURN
  END SUBROUTINE hsmt_sph
END MODULE m_hsmt_sph