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_qpwtonmt
8 : !***************************************************************
9 : ! This subroutine calculates the lattice harmonic expansions
10 : ! for the plane wave density inside the spheres.
11 : !
12 : ! Stefan Bl"ugel , IFF, Nov. 1997
13 : !***************************************************************
14 : CONTAINS
15 636 : SUBROUTINE qpw_to_nmt(sphhar,atoms,stars,sym,cell ,fmpi,jspin,l_cutoff,qpwc,rho)
16 : !
17 : USE m_constants
18 : USE m_phasy1
19 : USE m_sphbes
20 : USE m_types
21 :
22 :
23 : IMPLICIT NONE
24 : ! ..
25 : ! .. Scalar Arguments ..
26 : TYPE(t_sphhar),INTENT(IN) :: sphhar
27 : TYPE(t_atoms),INTENT(IN) :: atoms
28 : TYPE(t_stars),INTENT(IN) :: stars
29 : TYPE(t_cell),INTENT(IN) :: cell
30 : TYPE(t_sym),INTENT(IN) :: sym
31 :
32 : TYPE(t_mpi),INTENT(IN) :: fmpi
33 :
34 : INTEGER, INTENT (IN) :: jspin,l_cutoff
35 : ! ..
36 : ! .. Array Arguments ..
37 : COMPLEX, INTENT (IN) :: qpwc(stars%ng3)
38 : REAL, INTENT (INOUT) :: rho(:,0:,:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
39 : !-odim
40 : !+odim
41 : ! ..
42 : ! .. Local Scalars ..
43 : REAL, PARAMETER :: zero = 0.0
44 : COMPLEX, PARAMETER :: czero = (0.0,0.0)
45 : INTEGER in,j,jl,j0,jm,k,l,lh,m,n,nd,nrm,n1,n2,na,lm
46 : REAL d0,gr,r0,rr
47 : COMPLEX cp,sm,cprr2
48 : ! ..
49 : ! .. Local Arrays ..
50 636 : COMPLEX pylm( (atoms%lmaxd+1)**2,atoms%ntype ) !,bsl_c(atoms%jmtd,0:lmaxd)
51 636 : REAL bsl(0:atoms%lmaxd),bsl_r(atoms%jmtd,0:atoms%lmaxd),bsl_i(atoms%jmtd,0:atoms%lmaxd)
52 636 : INTEGER mr(atoms%ntype),lmx(atoms%ntype),lmxx(atoms%ntype),ntypsy_o(atoms%ntype)
53 : ! ..
54 636 : !$ REAL,ALLOCATABLE :: rho_tmp(:,:,:)
55 :
56 :
57 : !----> cut-off l-expansion of non-spherical charge contribution
58 : ! from coretails of neighboring atom for l> l_cutoff
59 : !
60 1858 : DO n = 1,atoms%ntype
61 1222 : na = atoms%firstAtom(n)
62 1222 : lmx(n) = MIN( atoms%lmax(n) , l_cutoff )
63 1858 : ntypsy_o(n) = sym%ntypsy(na)
64 : END DO
65 : !
66 : !----> identify atoms with the same radial mesh
67 : !
68 : j0 = 0
69 : r0 = zero
70 : d0 = zero
71 : nrm= 0
72 1858 : DO n = 1 , atoms%ntype
73 1222 : IF (.NOT.(atoms%jri(n).EQ.j0 .AND. atoms%rmsh(1,n).EQ.r0 &
74 : & .AND. atoms%dx(n).EQ.d0)) THEN
75 936 : j0 = atoms%jri(n)
76 936 : r0 = atoms%rmsh(1,n)
77 936 : d0 = atoms%dx(n)
78 936 : nrm= nrm + 1
79 936 : lmxx(nrm) = lmx(n)
80 : END IF
81 1222 : mr(nrm)=n
82 1858 : lmxx(nrm) = MAX( lmx(n) , lmx(nrm) )
83 : END DO
84 : !
85 : !=====> Loop over the g-vectors
86 : !
87 : ! ----> g=0 component
88 636 : CALL timestart("qpw_to_nmt")
89 636 : IF (fmpi%irank == 0) THEN
90 318 : cp = qpwc(1)*stars%nstr(1)
91 929 : DO n = 1 , atoms%ntype
92 423406 : DO j = 1,atoms%jri(n)
93 422477 : rr = atoms%rmsh(j,n)*atoms%rmsh(j,n)
94 423088 : rho(j,0,n,jspin) = rho(j,0,n,jspin) + rr*sfp_const*cp
95 : ENDDO
96 : ENDDO
97 : ELSE
98 15296783 : rho(:,:,:,jspin) = 0.0
99 : ENDIF
100 :
101 : ! ----> g.ne.0 components
102 : !
103 : ! g=|=0 terms : \sum_{g =|= 0} 4 \pi i^l \rho_int(g) r_i^{2} \times
104 : ! j_{l} (gr_i) \exp{ig\xi_i} Y^*_{lm} (g)
105 : !
106 : ! ----> calculate structure constant for each atom
107 : ! (4pi*i**l/nop(3)*sum(R){exp(iRG(taual-taur)*conjg(ylm(RG))
108 : !
109 : !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,cp,pylm,n1,in,n2,j) &
110 : !$OMP PRIVATE(cprr2,gr,bsl,jl,bsl_r,bsl_i,n,nd,lh,l,sm,jm,m,lm) &
111 636 : !$OMP PRIVATE(rho_tmp)
112 :
113 :
114 : !$ ALLOCATE(rho_tmp(size(rho,1),0:size(rho,2)-1,size(rho,3)))
115 : !$ rho_tmp=0
116 : !$OMP DO
117 : DO k = fmpi%irank+2,stars%ng3,fmpi%isize
118 : cp = qpwc(k)*stars%nstr(k)
119 :
120 : CALL phasy1(atoms,stars,sym,cell,k,pylm)
121 :
122 : !
123 : n1 = 1
124 : DO in = 1 , nrm
125 : n2 = mr(in)
126 : bsl_r = 0.0
127 : bsl_i = 0.0
128 : DO j = 1,atoms%jri(n1)
129 : cprr2 = cp*atoms%rmsh(j,n1)*atoms%rmsh(j,n1)
130 : gr = stars%sk3(k)*atoms%rmsh(j,n1)
131 : CALL sphbes(lmxx(in),gr,bsl)
132 : DO jl=0,lmxx(in)
133 : bsl_r(j,jl) = bsl(jl) * REAL(cprr2)
134 : bsl_i(j,jl) = bsl(jl) * AIMAG(cprr2)
135 : END DO
136 : END DO
137 : DO n = n1,n2
138 : nd = ntypsy_o(n)
139 : DO lh = 0,sphhar%nlh(nd)
140 : l = sphhar%llh(lh,nd)
141 : IF ( l.LE.lmx(n) ) THEN
142 : sm = czero
143 : DO jm = 1,sphhar%nmem(lh,nd)
144 : m = sphhar%mlh(jm,lh,nd)
145 : lm = l*(l+1) + m + 1
146 : sm = sm + CONJG(sphhar%clnu(jm,lh,nd)) *pylm(lm,n)
147 : END DO
148 : !$ if (.false.) THEN
149 : rho(:,lh,n,jspin) = rho(:,lh,n,jspin) + bsl_r(:,l) * REAL(sm)
150 : rho(:,lh,n,jspin) = rho(:,lh,n,jspin) - bsl_i(:,l) * AIMAG(sm)
151 : !$ endif
152 : !$ rho_tmp(:,lh,n)=rho_tmp(:,lh,n)+bsl_r(:,l)*real(sm)
153 : !$ rho_tmp(:,lh,n)=rho_tmp(:,lh,n)-bsl_i(:,l)*aimag(sm)
154 : END IF
155 : END DO
156 : END DO
157 : n1 = n2 + 1
158 : ENDDO
159 : ENDDO
160 : !$OMP END DO
161 : !$OMP CRITICAL
162 : !$ rho(:,:,:,jspin)=rho(:,:,:,jspin)+rho_tmp
163 : !$OMP END CRITICAL
164 : !$ DEALLOCATE(rho_tmp)
165 : !$OMP END PARALLEL
166 : !
167 636 : CALL timestop("qpw_to_nmt")
168 : !
169 636 : END SUBROUTINE qpw_to_nmt
170 : END MODULE m_qpwtonmt
|
