Line data Source code
1 : MODULE m_mcdinit
2 : CONTAINS
3 0 : SUBROUTINE mcd_init(atoms,banddos,input,vr,g,f,mcd,itype,jspin)
4 :
5 : !-----------------------------------------------------------------------
6 : !
7 : ! For a given atom-type 'itype' look, whether a core state is in the
8 : ! energy interval [emcd_lo,emcd_up] and, if found, calculate the
9 : ! MCD-matrix elements 'm_mcd'.
10 : !
11 : !-----------------------------------------------------------------------
12 :
13 : USE m_nabla
14 : USE m_dr2fdr
15 : USE m_constants, ONLY : c_light
16 : !USE m_setcor
17 : USE m_differ
18 : USE m_types
19 : use m_types_mcd
20 : IMPLICIT NONE
21 :
22 :
23 : TYPE(t_input),INTENT(IN) :: input
24 : TYPE(t_atoms),INTENT(IN) :: atoms
25 : TYPE(t_banddos),INTENT(IN) :: banddos
26 : TYPE(t_mcd),INTENT(INOUT) :: mcd
27 :
28 : INTEGER, PARAMETER :: l_max = 3
29 :
30 : ! Arguments ...
31 :
32 : INTEGER, INTENT (IN) :: itype
33 : INTEGER, INTENT (IN) :: jspin
34 : REAL, INTENT (IN) :: vr(atoms%jmtd,atoms%ntype,input%jspins)
35 : REAL, INTENT (IN) :: f(atoms%jmtd,2,0:atoms%lmaxd,input%jspins)
36 : REAL, INTENT (IN) :: g(atoms%jmtd,2,0:atoms%lmaxd,input%jspins)
37 :
38 : ! Locals ...
39 :
40 : INTEGER kap,mue,iri,l,ispin,i,icore,korb,nst,n_core,ierr,n_dos
41 : REAL c,t2,e,fj,fl,fn ,d,ms,rn
42 0 : INTEGER kappa(maxval(atoms%econf%num_states)),nprnc(maxval(atoms%econf%num_states)),l_core(maxval(atoms%econf%num_states))
43 0 : REAL vrd(atoms%msh),occ(maxval(atoms%econf%num_states),2),a(atoms%msh),b(atoms%msh),j_core(maxval(atoms%econf%num_states)),e_mcd1(maxval(atoms%econf%num_states))
44 0 : REAL gv1(atoms%jmtd)
45 0 : REAL, ALLOCATABLE :: gc(:,:,:),fc(:,:,:)
46 0 : REAL, ALLOCATABLE :: gv(:,:,:,:),fv(:,:,:,:),dgv(:,:,:,:)
47 :
48 : !-----------------------------------------------------------------------
49 :
50 0 : if (.not.any(banddos%dos_atom(atoms%firstAtom(itype):atoms%firstAtom(itype)+atoms%neq(itype)-1))) return
51 :
52 0 : c = c_light(1.0)
53 0 : ALLOCATE ( gc(atoms%jri(itype),atoms%econf(itype)%num_core_states,input%jspins) )
54 0 : ALLOCATE ( fc(atoms%jri(itype),atoms%econf(itype)%num_core_states,input%jspins) )
55 :
56 : ! core setup
57 0 : n_dos=banddos%map_atomtype(itype)
58 0 : mcd%ncore(n_dos) = 0
59 0 : CALL atoms%econf(itype)%get_core(nst,nprnc,kappa,occ)
60 :
61 0 : DO ispin = jspin, jspin
62 :
63 : ! extend core potential
64 :
65 0 : DO iri = 1, atoms%jri(itype)
66 0 : vrd(iri) = vr(iri,itype,ispin)
67 : ENDDO
68 0 : t2 = vrd(atoms%jri(itype)) / (atoms%jri(itype) - atoms%msh)
69 0 : DO iri = atoms%jri(itype) + 1, atoms%msh
70 0 : vrd(iri) = vrd(atoms%jri(itype)) + t2* ( iri-atoms%jri(itype) )
71 : ENDDO
72 :
73 : ! calculate core
74 :
75 0 : n_core = 0
76 0 : DO korb = 1, atoms%econf(itype)%num_core_states
77 0 : IF (occ(korb,1).GT.0) THEN
78 0 : fn = nprnc(korb)
79 0 : fj = iabs(kappa(korb)) - .5e0
80 0 : fl = fj + (.5e0)*isign(1,kappa(korb))
81 0 : e = -2* (atoms%zatom(itype)/ (fn+fl))**2
82 0 : d = EXP(atoms%dx(itype))
83 0 : rn = atoms%rmsh(1,itype)*( d**(atoms%msh-1) )
84 : CALL differ(fn,fl,fj,c,atoms%zatom(itype),atoms%dx(itype),atoms%rmsh(1,itype),&
85 0 : rn,d,atoms%msh,vrd, e, a,b,ierr)
86 0 : IF (ierr/=0) CALL juDFT_error("error in core-levels", calledby="mcd_init")
87 0 : IF ( (e.LE.mcd%emcd_up).AND.(e.GE.mcd%emcd_lo) ) THEN
88 0 : WRITE(*,*) 'good ev = ',e
89 0 : n_core = n_core + 1
90 0 : j_core(n_core) = fj
91 0 : l_core(n_core) = NINT( fl )
92 0 : e_mcd1(n_core) = e
93 0 : DO iri = 1, atoms%jri(itype)
94 0 : gc(iri,n_core,ispin) = a(iri)
95 0 : fc(iri,n_core,ispin) = b(iri)
96 : ENDDO
97 : ENDIF
98 : ENDIF
99 : ENDDO
100 :
101 : ENDDO
102 :
103 : !-----------------------------------------------------------------------
104 :
105 0 : IF (n_core.GT.0) THEN
106 :
107 0 : ALLOCATE ( gv(atoms%jri(itype),0:l_max,input%jspins,2) )
108 0 : ALLOCATE (dgv(atoms%jri(itype),0:l_max,input%jspins,2) )
109 0 : ALLOCATE ( fv(atoms%jri(itype),0:l_max,input%jspins,2) )
110 0 : DO i = 1, 2
111 0 : DO iri = 3*(n_dos-1)+1 , 3*(n_dos-1)+3
112 0 : DO l = 1, (l_max+1)**2
113 0 : DO icore = 1, maxval(atoms%econf%num_states)
114 0 : mcd%m_mcd(icore,l,iri,i) = CMPLX(0.0,0.0)
115 : ENDDO
116 : ENDDO
117 : ENDDO
118 : ENDDO
119 : !
120 : ! bring LAPW wavefunctions in a proper form:
121 : !
122 0 : DO ispin = jspin, jspin
123 0 : ms = ispin - 1.5
124 0 : DO l = 0, l_max
125 0 : DO iri = 1, atoms%jri(itype)
126 0 : gv(iri,l,ispin,1) = f(iri,1,l,ispin) ! large component of u
127 0 : fv(iri,l,ispin,1) = f(iri,2,l,ispin) ! small .
128 0 : gv(iri,l,ispin,2) = g(iri,1,l,ispin) ! large component of u
129 0 : fv(iri,l,ispin,2) = g(iri,2,l,ispin) ! small
130 : ENDDO
131 0 : gv1(:) = atoms%rmsh(:,itype) * gv(:,l,ispin,1)
132 : CALL dr2fdr(& ! deriative of u (large)&
133 0 : gv1,atoms%rmsh(1,itype),atoms%jri(itype), dgv(1,l,ispin,1) )
134 0 : gv1(:) = atoms%rmsh(:,itype) * gv(:,l,ispin,2) ! .
135 : CALL dr2fdr(& ! deriative of u (large)&
136 0 : gv1,atoms%rmsh(1,itype),atoms%jri(itype), dgv(1,l,ispin,2) )
137 : ENDDO
138 : !
139 : !
140 : !
141 0 : DO icore = 1, n_core
142 :
143 0 : DO i = 1, 2
144 : ! write(*,*) j_core(icore),l_core(icore),l_max,ms
145 : CALL nabla(n_dos,icore,atoms%jri(itype),atoms%dx(itype),maxval(atoms%econf%num_states),atoms%ntype,&
146 : j_core(icore),l_core(icore),l_max,ms,atoms%rmsh(:,itype),gc(:,icore,ispin),&
147 0 : gv(:,0:,ispin,i),dgv(:,0:,ispin,i), mcd%m_mcd(:,:,:,i) )
148 : ENDDO
149 :
150 0 : DO i = 1, 2*icore*l_core(icore)
151 0 : mcd%ncore(n_dos) = mcd%ncore(n_dos) + 1
152 0 : IF (mcd%ncore(n_dos)>maxval(atoms%econf%num_states)) CALL juDFT_error("maxval(atoms%econf%num_states) too small" ,calledby ="mcd_init")
153 0 : mcd%e_mcd(n_dos,ispin,mcd%ncore(n_dos)) = e_mcd1(icore)
154 0 : IF (e_mcd1(icore).GT.mcd%maxE_mcd) mcd%maxE_mcd = e_mcd1(icore)
155 : ENDDO
156 : ENDDO
157 : ENDDO
158 :
159 0 : DEALLOCATE (gv,fv,dgv)
160 : ENDIF
161 0 : DEALLOCATE (gc,fc)
162 :
163 :
164 : ! DO i = 1, 2
165 : ! DO iri = 3*(itype-1)+1 , 3*(itype-1)+3
166 : ! write (*,*) iri
167 : ! DO icore = 1, mcd%ncore(itype)
168 : ! write (*,'(10f10.5)') (mcd%m_mcd(icore,l,iri,i),l=1,9)
169 : ! ENDDO
170 : ! ENDDO
171 : ! ENDDO
172 : END SUBROUTINE mcd_init
173 : END MODULE m_mcdinit
|
