Line data Source code
1 : module m_calc_cmt
2 :
3 : contains
4 112 : subroutine calc_cmt(atoms, cell, input, noco, nococonv, hybinp, hybdat, mpdata, kpts, &
5 112 : sym, zmat_ikp, jsp, ik, c_phase, cmt_out, submpi)
6 : use m_types
7 : use m_judft
8 : USE m_hyb_abcrot
9 : USE m_abcof
10 : use m_constants
11 : use m_trafo, only: waveftrafo_gen_cmt
12 : use m_io_hybrid
13 : use m_divide_most_evenly
14 : implicit none
15 : type(t_atoms), intent(in) :: atoms
16 : type(t_cell), intent(in) :: cell
17 : type(t_input), intent(in) :: input
18 : type(t_noco), intent(in) :: noco
19 : type(t_nococonv), intent(in) :: nococonv
20 : type(t_hybinp), intent(in) :: hybinp
21 : type(t_hybdat), intent(in) :: hybdat
22 : type(t_mpdata), intent(in) :: mpdata
23 : type(t_kpts), intent(in) :: kpts
24 : type(t_sym), intent(in) :: sym
25 :
26 : type(t_mat), intent(in), target :: zmat_ikp ! zmat of parent k-point
27 : ! not sure wether zmat works aswell
28 : integer, intent(in) :: jsp
29 : integer, intent(in) :: ik ! k-point
30 : complex, intent(in) :: c_phase(:)
31 :
32 : complex, intent(inout) :: cmt_out(:,:,:)
33 : type(t_hybmpi), intent(in), optional :: submpi
34 112 : complex, allocatable :: acof(:,:,:), bcof(:,:,:), ccof(:,:,:,:)
35 112 : complex, allocatable :: cmt(:,:,:)
36 448 : type(t_noco) :: nocoHyb
37 :
38 : integer :: ikp, nbands, ok(4) ! index of parent k-point
39 : integer :: iatom, iatom2, itype, indx, i, j, idum, iop, l, ll, lm, m, lm1, lm2
40 112 : integer :: map_lo(atoms%nlod)
41 112 : integer, allocatable :: start_idx(:), psize(:)
42 : integer :: my_psz, my_start, ierr
43 :
44 : REAL :: rdum, rfac
45 : complex :: cdum, cexp1, cexp2, cfac
46 112 : COMPLEX, ALLOCATABLE :: cmtTemp1(:,:), cmtTemp2(:,:)
47 112 : type(t_lapw) :: lapw_ik, lapw_ikp
48 112 : type(t_mat), target :: tmp
49 : type(t_mat), pointer :: mat_ptr
50 :
51 112 : call timestart("calc_cmt")
52 112 : ikp = kpts%bkp(ik)
53 112 : nbands = zmat_ikp%matsize2
54 :
55 112 : if(present(submpi)) then
56 24 : call divide_most_evenly(nbands, submpi%size, start_idx, psize)
57 24 : my_psz = psize(submpi%rank+1)
58 24 : my_start = start_idx(submpi%rank+1)
59 : else
60 88 : my_psz = nbands
61 88 : my_start = 1
62 : endif
63 :
64 112 : call timestart("alloc abccof")
65 299018 : allocate(acof(my_psz, 0:atoms%lmaxd*(atoms%lmaxd+2), atoms%nat), stat=ok(1), source=cmplx_0)
66 298906 : allocate(bcof(my_psz, 0:atoms%lmaxd*(atoms%lmaxd+2), atoms%nat), stat=ok(2), source=cmplx_0)
67 23368 : allocate(ccof(-atoms%llod:atoms%llod, my_psz, atoms%nlod, atoms%nat), stat=ok(3), source=cmplx_0)
68 611568 : allocate(cmt(nbands, hybdat%maxlmindx, atoms%nat), stat=ok(4), source=cmplx_0)
69 560 : if(any(ok /= 0)) then
70 0 : call judft_error("Error in allocating abcof arrays")
71 : endif
72 112 : call timestop("alloc abccof")
73 :
74 112 : CALL lapw_ik%init(input, noco, nococonv, kpts, atoms, sym, ik, cell)
75 112 : CALL lapw_ikp%init(input, noco, nococonv, kpts, atoms, sym, ikp, cell)
76 :
77 112 : lapw_ikp%nmat = lapw_ikp%nv(jsp) + atoms%nlotot
78 112 : if(my_psz /= nbands) then
79 0 : call tmp%init(zmat_ikp%l_real, zmat_ikp%matsize1, my_psz)
80 0 : if(tmp%l_real) then
81 0 : tmp%data_r = zmat_ikp%data_r(:,my_start:my_start+my_psz-1)
82 : else
83 0 : tmp%data_c = zmat_ikp%data_c(:,my_start:my_start+my_psz-1)
84 : endif
85 : mat_ptr => tmp
86 : else
87 : mat_ptr => zmat_ikp
88 : endif
89 :
90 112 : nocoHyb = noco
91 112 : IF (hybinp%l_hybrid .AND. noco%l_soc) nocoHyb%l_soc = .FALSE.
92 112 : CALL abcof(input, atoms, sym, cell, lapw_ikp, my_psz, hybdat%usdus, nocoHyb, nococonv, jsp, acof, bcof, ccof, mat_ptr)
93 :
94 112 : CALL hyb_abcrot(hybinp, atoms, my_psz, sym, acof, bcof, ccof)
95 :
96 112 : call timestart("copy to cmt")
97 : !$OMP parallel do default(none) private(iatom, itype, indx, l, ll, cdum, idum, map_lo, j, m, lm, i) &
98 112 : !$OMP shared(atoms, mpdata, cmt, acof, bcof, ccof, my_start, my_psz)
99 : DO iatom = 1,atoms%nat
100 : itype = atoms%itype(iatom)
101 :
102 : indx = 0
103 : DO l = 0, atoms%lmax(itype)
104 : ll = l*(l + 1)
105 : cdum = ImagUnit**l
106 :
107 : ! determine number of local orbitals with quantum number l
108 : ! map returns the number of the local orbital of quantum
109 : ! number l in the list of all local orbitals of the atom type
110 : idum = 0
111 : map_lo = 0
112 : IF (mpdata%num_radfun_per_l(l, itype) > 2) THEN
113 : DO j = 1, atoms%nlo(itype)
114 : IF (atoms%llo(j, itype) == l) THEN
115 : idum = idum + 1
116 : map_lo(idum) = j
117 : END IF
118 : END DO
119 : END IF
120 :
121 : DO m = -l, l
122 : lm = ll + m
123 : DO i = 1, mpdata%num_radfun_per_l(l, itype)
124 : indx = indx + 1
125 : IF (i == 1) THEN
126 : cmt(my_start:my_start+my_psz-1, indx, iatom) = cdum*acof(:, lm, iatom)
127 : ELSE IF (i == 2) THEN
128 : cmt(my_start:my_start+my_psz-1, indx, iatom) = cdum*bcof(:, lm, iatom)
129 : ELSE
130 : idum = i - 2
131 : cmt(my_start:my_start+my_psz-1, indx, iatom) = cdum*ccof(m, :, map_lo(idum), iatom)
132 : END IF
133 : END DO
134 : END DO
135 : END DO
136 : END DO
137 : !$OMP end parallel do
138 :
139 112 : call timestop("copy to cmt")
140 :
141 : #ifdef CPP_MPI
142 112 : if(my_psz /= nbands) then
143 0 : call timestart("allreduce cmt")
144 0 : do i = 1,atoms%nat
145 0 : call MPI_Allreduce(MPI_IN_PLACE, cmt(1,1,i), nbands*hybdat%maxlmindx, MPI_DOUBLE_COMPLEX, MPI_SUM, submpi%comm, ierr)
146 : enddo
147 0 : call timestop("allreduce cmt")
148 : endif
149 : #endif
150 :
151 : ! write cmt at irreducible k-points in direct-access file cmt
152 112 : if(ik <= kpts%nkpt) then
153 84 : call timestart("copy out cmt")
154 336 : call zcopy(size(cmt), cmt, 1, cmt_out, 1)
155 84 : call timestop("copy out cmt")
156 : else
157 28 : iop = kpts%bksym(ik)
158 :
159 : call waveftrafo_gen_cmt(cmt, c_phase, zmat_ikp%l_real, ikp, iop, atoms, &
160 28 : mpdata, hybinp, kpts, sym, nbands, cmt_out)
161 : endif
162 :
163 : ! calculate real values cmt for atoms that are inversion symmetric to each other
164 112 : rfac = 1.0 / SQRT(2.0)
165 112 : cfac = CMPLX(0.0,-1.0) / SQRT(2.0)
166 672 : ALLOCATE(cmtTemp1(SIZE(cmt,1),SIZE(cmt,2)), cmtTemp2(SIZE(cmt,1),SIZE(cmt,2)))
167 280 : DO iatom = 1,atoms%nat
168 168 : itype = atoms%itype(iatom)
169 168 : iatom2 = sym%invsatnr(iatom)
170 168 : IF(iatom2.EQ.0) iatom2 = iatom
171 :
172 280 : IF(iatom2.GT.iatom) THEN ! iatom is first of two atoms mapped to each other via inversion symmetry
173 0 : cmtTemp1 = CMPLX(0.0,0.0)
174 0 : cmtTemp2 = CMPLX(0.0,0.0)
175 0 : cexp1 = exp(CMPLX(0.0,-1.0)*tpi_const*dot_product(kpts%bkf(:, ik), atoms%taual(:, iatom)))
176 0 : cexp2 = exp(CMPLX(0.0,-1.0)*tpi_const*dot_product(kpts%bkf(:, ik), atoms%taual(:, iatom2)))
177 : ! Note: for the case k+q these phase calculations are slightly different then in the old version, where the resultin k+q point may be outside the BZ.
178 :
179 0 : lm1 = 0
180 0 : DO l = 0, atoms%lmax(itype)
181 0 : DO m = -l, l
182 0 : rdum = (-1)**(l + m)
183 0 : DO i = 1, mpdata%num_radfun_per_l(l, itype)
184 0 : lm1 = lm1 + 1
185 : ! lm index at l,-m
186 0 : lm2 = lm1 - 2 * m * mpdata%num_radfun_per_l(l, itype)
187 0 : cdum = rdum * cexp1 * cexp2
188 :
189 0 : cmtTemp1(:,lm1) = (cmt_out(:, lm1, iatom) + cdum*cmt_out(:, lm2, iatom2))*rfac
190 0 : cmtTemp2(:,lm1) = (cmt_out(:, lm1, iatom) - cdum*cmt_out(:, lm2, iatom2))*cfac
191 : END DO
192 : END DO
193 : END DO
194 0 : cmt_out(:,:,iatom) = cmtTemp1(:,:)
195 0 : cmt_out(:,:,iatom2) = cmtTemp2(:,:)
196 : END IF
197 : END DO
198 112 : DEALLOCATE (cmtTemp1, cmtTemp2)
199 :
200 112 : call timestop("calc_cmt")
201 112 : end subroutine calc_cmt
202 : end module m_calc_cmt
|
