Line data Source code
1 : MODULE m_vham
2 :
3 : !----------------------------------------------------------------------------------- !
4 : ! !
5 : ! MODULE: m_v_ham !
6 : ! !
7 : ! Author : Wejdan Beida !
8 : ! !
9 : ! Description: This module calculates the intersite occupation matrix !
10 : ! !
11 : ! s s --- --- IJ JI s I J s !
12 : !<Phsi |V_POT_MAT| Phis > = - > > V n <Psi |phi > <phi |Psi > !
13 : ! k k --- --- k L Lp k !
14 : ! I,J,s L,Lp !
15 : ! !
16 : ! where L={l,m},Lp={lp,mp} !
17 : ! !
18 : !------------------------------------------------------------------------------------!
19 :
20 :
21 : CONTAINS
22 :
23 0 : SUBROUTINE v_ham(input,usdus,atoms,kpts,cell,lapw,sym,noco,fmpi,nococonv,fjgj,den,jspin,kptindx,hmat)
24 :
25 : USE m_types
26 : USE m_constants
27 : USE m_juDFT
28 : USE m_hsmt_ab
29 : USE m_hsmt_fjgj
30 : USE m_ylm
31 : USE m_radsrd
32 :
33 : IMPLICIT NONE
34 :
35 : TYPE(t_input), INTENT(IN) :: input
36 : TYPE(t_usdus), INTENT(IN) :: usdus
37 : TYPE(t_atoms), INTENT(IN) :: atoms
38 : TYPE(t_kpts), INTENT(IN) :: kpts
39 : TYPE(t_cell), INTENT(IN) :: cell
40 : TYPE(t_lapw), INTENT(IN) :: lapw
41 : TYPE(t_sym), INTENT(IN) :: sym
42 : TYPE(t_noco), INTENT(IN) :: noco
43 : TYPE(t_mpi), INTENT(IN) :: fmpi
44 : TYPE(t_nococonv), INTENT(IN) :: nococonv
45 : TYPE(t_fjgj), INTENT(INOUT) :: fjgj
46 : TYPE(t_potden), INTENT(IN) :: den
47 : INTEGER, INTENT(IN) :: jspin,kptindx
48 : CLASS(t_mat), INTENT(INOUT) :: hmat
49 :
50 :
51 :
52 :
53 : INTEGER i_v,i_pair,natom1,latom1,ll1atom1,atom2,natom2,latom2,ll1atom2,matom1,matom2,lm1atom1,lm1atom2,iG1,iG2,abSizeG1,abSizeG2,indx_pair,counter, atom1index, atom2index
54 : COMPLEX c_0,c_pair, a1, b1, a2, b2, power_fac, exponent
55 : REAL norm1_W, norm2_W, V_inp
56 0 : COMPLEX, ALLOCATABLE :: abG1(:,:),abG2(:,:)
57 0 : COMPLEX, ALLOCATABLE :: c_0_pair(:,:,:)
58 0 : ALLOCATE(abG1(2*atoms%lmaxd*(atoms%lmaxd+2)+2,MAXVAL(lapw%nv)))
59 0 : ALLOCATE(abG2(2*atoms%lmaxd*(atoms%lmaxd+2)+2,MAXVAL(lapw%nv)))
60 :
61 0 : counter=0
62 0 : DO i_v = 1,atoms%n_v
63 0 : Do atom2=1,atoms%lda_v(i_v)%numOtherAtoms
64 0 : counter=counter+1
65 : ENDDO
66 : ENDDO
67 0 : ALLOCATE(c_0_pair(counter,MAXVAL(lapw%nv),MAXVAL(lapw%nv)))
68 :
69 0 : DO iG1=1,lapw%nv(jspin)
70 0 : DO iG2=1,lapw%nv(jspin)
71 0 : DO indx_pair=1, counter-1
72 0 : c_0_pair(indx_pair,iG1,iG2)=cmplx_0
73 : ENDDO
74 : ENDDO
75 : ENDDO
76 :
77 : i_pair=1
78 0 : DO i_v = 1,atoms%n_v
79 0 : V_inp=atoms%lda_v(i_v)%V / hartree_to_ev_const
80 0 : natom1=atoms%lda_v(i_v)%atomIndex
81 0 : latom1=atoms%lda_v(i_v)%thisAtomL
82 0 : ll1atom1=latom1*(latom1+1)
83 0 : norm1_W = usdus%ddn(latom1,atoms%itype(natom1),jspin)**0.5
84 0 : CALL fjgj%calculate(input,atoms,cell,lapw,noco,usdus,atoms%itype(natom1),jspin)
85 0 : CALL hsmt_ab(sym,atoms,noco,nococonv,jspin,1,atoms%itype(natom1),natom1,cell,lapw,fjgj,abG1,abSizeG1,.FALSE.)
86 0 : Do atom2=1,atoms%lda_v(i_v)%numOtherAtoms
87 0 : natom2=atoms%lda_v(i_v)%otherAtomIndices(atom2)
88 0 : latom2=atoms%lda_v(i_v)%otherAtomL
89 0 : ll1atom2=latom2*(latom2+1)
90 0 : norm2_W = usdus%ddn(latom2,atoms%itype(natom2),jspin)**0.5
91 0 : power_fac=(cmplx(0, -1)**latom1) *(cmplx(0, 1)**latom2)
92 0 : CALL fjgj%calculate(input,atoms,cell,lapw,noco,usdus,atoms%itype(natom2),jspin)
93 0 : CALL hsmt_ab(sym,atoms,noco,nococonv,jspin,1,atoms%itype(natom2),natom2,cell,lapw,fjgj,abG2,abSizeG2,.FALSE.)
94 0 : DO iG1=1,lapw%nv(jspin)
95 0 : Do iG2=1,lapw%nv(jspin)
96 0 : c_0=cmplx_0
97 0 : exponent=EXP(cmplx(0.0,tpi_const)*dot_product(atoms%lda_v(i_v)%atomShifts(:,atom2),(kpts%bk(:,kptindx)+lapw%gvec(:, iG2,jspin))))
98 0 : Do matom1=-latom1,latom1
99 0 : lm1atom1=ll1atom1+matom1
100 0 : Do matom2=-latom2,latom2
101 0 : lm1atom2=ll1atom2+matom2
102 0 : a1 = abG1(lm1atom1+1,iG1)
103 0 : b1 = abG1(lm1atom1+1+abSizeG1/2,iG1)
104 0 : a2 = abG2(lm1atom2+1,iG2)
105 0 : b2 = abG2(lm1atom2+1+abSizeG2/2,iG2)
106 : c_0 = c_0 - (V_inp)*(conjg(den%nIJ_llp_mmp(matom1,matom2,i_pair,jspin)))* exponent *power_fac * &
107 0 : (conjg(a1)*a2 + conjg(b1)*a2*norm1_W + conjg(a1)*b2*norm2_W + conjg(b1)*b2*norm2_W*norm1_W)
108 : ENDDO
109 : ENDDO
110 0 : c_0_pair(i_pair,iG1,iG2)=c_0
111 : ENDDO
112 : ENDDO
113 0 : i_pair=i_pair+1
114 : ENDDO
115 : ENDDO
116 :
117 0 : DO iG1=1,lapw%nv(jspin)
118 0 : Do iG2=1,lapw%nv(jspin)
119 : c_pair=cmplx_0
120 0 : DO indx_pair=1, i_pair-1
121 0 : c_pair=c_pair+c_0_pair(indx_pair,iG1,iG2)
122 : ENDDO
123 0 : IF(hmat%l_real) THEN
124 0 : hmat%data_r(iG1,iG2)=hmat%data_r(iG1,iG2)+REAL(c_pair)
125 : ELSE
126 0 : hmat%data_c(iG1,iG2)=hmat%data_c(iG1,iG2)+c_pair
127 : END IF
128 : ENDDO
129 : ENDDO
130 :
131 :
132 : !INTEGER i_v,i_pair,natom1,latom1,ll1atom1,atom2,natom2,latom2,ll1atom2,iG2,abSizeG1,abSizeG2,counter, atom1index, atom2index, i,j,k
133 : !COMPLEX power_fac, exponent
134 : !REAL norm1_W, norm2_W, V_inp
135 : !COMPLEX, ALLOCATABLE :: abG1(:,:),abG2(:,:),X1(:,:), X2(:,:), PotMat(:,:)
136 : !ALLOCATE(abG1(2*atoms%lmaxd*(atoms%lmaxd+2)+2,MAXVAL(lapw%nv)))
137 : !ALLOCATE(abG2(2*atoms%lmaxd*(atoms%lmaxd+2)+2,MAXVAL(lapw%nv)))
138 : !atom1index=0
139 : !atom2index=0
140 : !counter=0
141 : !DO i_v = 1,atoms%n_v
142 : ! Do atom2=1,atoms%lda_v(i_v)%numOtherAtoms
143 : ! counter=counter+1
144 : ! atom1index= atom1index + 2 * ( 2 * atoms%lda_v(i_v)%thisAtomL +1 )
145 : ! atom2index= atom2index + 2 * ( 2 * atoms%lda_v(i_v)%otherAtomL +1 )
146 : ! ENDDO
147 : !ENDDO
148 : !ALLOCATE(X1(atom1index,lapw%nv(jspin)))
149 : !ALLOCATE(X2(atom2index,lapw%nv(jspin)))
150 : !ALLOCATE(PotMat(atom1index,atom2index))
151 :
152 : !DO i=1,atom1index
153 : ! DO j=1,atom2index
154 : ! PotMat(i,j)= cmplx_0
155 : ! DO k =1, MAXVAL(lapw%nv)
156 : ! X1(i,k)=cmplx_0
157 : ! X2(j,k)=cmplx_0
158 : ! ENDDO
159 : ! ENDDO
160 : !ENDDO
161 :
162 : !i_pair=1
163 : !atom1index=0
164 : !atom2index=0
165 : !counter=0
166 : !DO i_v = 1,atoms%n_v
167 : ! V_inp=atoms%lda_v(i_v)%V / hartree_to_ev_const
168 : ! natom1=atoms%lda_v(i_v)%atomIndex
169 : ! latom1=atoms%lda_v(i_v)%thisAtomL
170 : ! ll1atom1=latom1*(latom1+1)
171 : ! norm1_W = usdus%ddn(latom1,atoms%itype(natom1),jspin)**0.5
172 : ! CALL fjgj%calculate(input,atoms,cell,lapw,noco,usdus,atoms%itype(natom1),jspin)
173 : ! CALL hsmt_ab(sym,atoms,noco,nococonv,jspin,1,atoms%itype(natom1),natom1,cell,lapw,fjgj,abG1,abSizeG1,.FALSE.)
174 : ! DO atom2=1,atoms%lda_v(i_v)%numOtherAtoms
175 : ! natom2=atoms%lda_v(i_v)%otherAtomIndices(atom2)
176 : ! latom2=atoms%lda_v(i_v)%otherAtomL
177 : ! ll1atom2=latom2*(latom2+1)
178 : ! norm2_W = usdus%ddn(latom2,atoms%itype(natom2),jspin)**0.5
179 : ! power_fac=(cmplx(0, -1)**latom1) *(cmplx(0, 1)**latom2)
180 : ! X1( atom1index +1 : atom1index + 2 * latom1 + 1,:) = conjg(abG1(ll1atom1-latom1+1:ll1atom1+latom1+1,:))*power_fac*V_inp
181 : ! X1( atom1index + 2*latom1 +1 +1: atom1index + 2*(2*latom1 +1) ,:)= conjg(abG1(ll1atom1-latom1+abSizeG1/2+1:ll1atom1+latom1+abSizeG1/2+1,:))*norm1_W*power_fac*V_inp
182 : ! CALL fjgj%calculate(input,atoms,cell,lapw,noco,usdus,atoms%itype(natom2),jspin)
183 : ! CALL hsmt_ab(sym,atoms,noco,nococonv,jspin,1,atoms%itype(natom2),natom2,cell,lapw,fjgj,abG2,abSizeG2,.FALSE.)
184 : ! DO iG2= 1,lapw%nv(jspin)
185 : ! exponent=EXP(cmplx(0.0,tpi_const)*dot_product(atoms%lda_v(i_v)%atomShifts(:,atom2),(kpts%bk(:,kptindx)+lapw%gvec(:, iG2,jspin))))
186 : ! X2( atom2index +1 : atom2index + 2 * latom2 + 1,:) = abG2(ll1atom2-latom2+1:ll1atom2+latom2+1,:)*exponent
187 : ! X2( atom2index + 2*latom2 +1 +1: atom2index + 2*(2*latom2 +1) ,:)= abG2(ll1atom2-latom2+abSizeG2/2+1:ll1atom2+latom2+abSizeG2/2+1,:)*exponent*norm2_W
188 : ! ENDDO
189 : ! PotMat(1+atom1index:atom1index+(2*latom1 +1),1+atom2index:atom2index+(2*latom2 +1)) = conjg(den%nIJ_llp_mmp(-latom1:latom1,-latom2:latom2,i_pair,jspin))
190 : ! PotMat(1+atom1index+(2*latom1 +1):atom1index+2*(2*latom1 +1),1+atom2index+(2*latom2 +1):atom2index+2*(2*latom2 +1)) = conjg(den%nIJ_llp_mmp(-latom1:latom1,-latom2:latom2,i_pair,jspin))
191 : ! atom2index = atom2index + 2*(2*latom2 +1)
192 : ! atom1index = atom1index + 2*(2*latom1 +1)
193 : ! i_pair=i_pair +1
194 : ! ENDDO
195 : !ENDDO
196 : !IF(hmat%l_real) THEN
197 : ! hmat%data_r=hmat%data_r - REAL(MATMUL(TRANSPOSE(X1), MATMUL(TRANSPOSE(PotMat),X2)))
198 : !ELSE
199 : ! hmat%data_c=hmat%data_c - MATMUL(TRANSPOSE(X1), MATMUL(TRANSPOSE(PotMat),X2))
200 : !END IF
201 :
202 : ! For debugging (END)
203 : !i_pair=1
204 : !WRITE(1000+kptindx,*) '=========='
205 : !DO i_v = 1,atoms%n_v
206 : ! natom1=atoms%lda_v(i_v)%atomIndex
207 : ! latom1=atoms%lda_v(i_v)%thisAtomL
208 : ! Do atom2=1,atoms%lda_v(i_v)%numOtherAtoms
209 : ! natom2=atoms%lda_v(i_v)%otherAtomIndices(atom2)
210 : ! latom2=atoms%lda_v(i_v)%otherAtomL
211 : ! DO matom1=-latom1,latom1
212 : ! DO matom2=-latom2,latom2
213 : ! WRITE(1000+kptindx,'(4i7,2f15.8)') matom2,matom1,i_Pair,jspin,den%nIJ_llp_mmp(matom1,matom2,i_pair,jspin)
214 : ! END DO
215 : ! END DO
216 : ! i_pair=i_pair+1
217 : ! END DO
218 : !END DO
219 : ! For debugging (END)
220 :
221 : ! T5=0.0
222 : ! T6=0.0
223 : ! call cpu_time(T5)
224 : ! call cpu_time(T6)
225 : ! WRITE(20,*) 'intenal loop call', T6-T5
226 0 : END SUBROUTINE v_ham
227 : END MODULE m_vham
|
