Line data Source code
1 : MODULE m_coredr
2 : CONTAINS
3 2 : SUBROUTINE coredr(input,atoms,seig, rho,sphhar, vrs, qints,rhc)
4 : ! *******************************************************
5 : ! ***** set up the core densities for compounds *****
6 : ! ***** for relativistic core *****
7 : ! *******************************************************
8 :
9 : USE m_etabinit
10 : USE m_spratm
11 : USE m_ccdnup
12 : USE m_cdn_io
13 : USE m_types
14 : IMPLICIT NONE
15 :
16 : TYPE(t_input),INTENT(IN) :: input
17 : TYPE(t_sphhar),INTENT(IN) :: sphhar
18 : TYPE(t_atoms),INTENT(IN) :: atoms
19 : !
20 : ! .. Scalar Arguments ..
21 : REAL seig
22 : ! ..
23 : ! .. Array Arguments ..
24 : REAL , INTENT (IN) :: vrs(atoms%jmtd,atoms%ntype,input%jspins)
25 : REAL, INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
26 : REAL, INTENT (OUT) :: rhc(atoms%msh,atoms%ntype,input%jspins),qints(atoms%ntype,input%jspins)
27 : ! ..
28 : ! .. Local Scalars ..
29 : REAL dxx,rnot,sume,t2,t2b,z,t1,rr,d,v1,v2
30 : INTEGER i,j,jatom,jspin,k,n,ncmsh
31 : LOGICAL exetab
32 : ! ..
33 : ! .. Local Arrays ..
34 2 : REAL br(atoms%jmtd,atoms%ntype),brd(atoms%msh),etab(100,atoms%ntype),&
35 2 : rhcs(atoms%jmtd,atoms%ntype,input%jspins),rhochr(atoms%msh),rhospn(atoms%msh),&
36 2 : tecs(atoms%ntype,input%jspins),vr(atoms%jmtd,atoms%ntype),vrd(atoms%msh)
37 2 : INTEGER nkmust(atoms%ntype),ntab(100,atoms%ntype),ltab(100,atoms%ntype)
38 :
39 : ! ..
40 608 : ntab(:,:) = -1 ; ltab(:,:) = -1 ; etab(:,:) = 0.0
41 : !
42 : ! setup potential and field
43 : !
44 2 : IF (input%jspins.EQ.1) THEN
45 0 : DO n = 1,atoms%ntype
46 0 : DO j = 1,atoms%jmtd
47 0 : vr(j,n) = vrs(j,n,1)
48 0 : br(j,n) = 0.0
49 : END DO
50 : END DO
51 : ELSE
52 4 : DO n = 1,atoms%ntype
53 1134 : DO j = 1,atoms%jmtd
54 1130 : vr(j,n) = (vrs(j,n,1)+vrs(j,n,input%jspins))/2.
55 1132 : br(j,n) = (vrs(j,n,input%jspins)-vrs(j,n,1))/2.
56 : END DO
57 : END DO
58 : END IF
59 : !
60 : ! setup eigenvalues
61 2 : exetab = .FALSE.
62 2 : INQUIRE (file='core.dat',exist=exetab)
63 2 : IF (exetab) THEN
64 1 : OPEN (58,file='core.dat',form='formatted',status='old')
65 1 : REWIND 58
66 2 : DO n = 1,atoms%ntype
67 1 : READ (58,FMT=*) nkmust(n)
68 20 : DO k = 1,nkmust(n)
69 18 : READ (58,FMT='(f12.6,2i3)') etab(k,n),ntab(k,n),&
70 37 : & ltab(k,n)
71 :
72 : END DO
73 : END DO
74 : ELSE
75 1 : OPEN (58,file='core.dat',form='formatted',status='new')
76 1 : CALL etabinit(atoms,input, vr, etab,ntab,ltab,nkmust)
77 : END IF
78 : !
79 2 : ncmsh = atoms%msh
80 2 : seig = 0.
81 : ! ---> set up densities
82 4 : DO jatom = 1,atoms%ntype
83 : !
84 1132 : DO j = 1,atoms%jri(jatom)
85 1130 : vrd(j) = vr(j,jatom)
86 1132 : brd(j) = br(j,jatom)
87 : END DO
88 :
89 2 : IF (input%l_core_confpot) THEN
90 : !---> linear extension of the potential with slope t1 / a.u.
91 2 : rr = atoms%rmt(jatom)
92 2 : d = EXP(atoms%dx(jatom))
93 2 : t1=0.125
94 : ! t2 = vrd(jri(jatom))/rr - rr*t1
95 : ! t2b = brd(jri(jatom))/rr - rr*t1
96 2 : t2 = vrs(atoms%jri(jatom),jatom,1) /rr - rr*t1
97 2 : t2b = vrs(atoms%jri(jatom),jatom,input%jspins)/rr - rr*t1
98 : ELSE
99 0 : t2 = vrd(atoms%jri(jatom))/ (atoms%jri(jatom)-ncmsh)
100 0 : t2b = brd(atoms%jri(jatom))/ (atoms%jri(jatom)-ncmsh)
101 : ENDIF
102 2 : IF (atoms%jri(jatom).LT.ncmsh) THEN
103 218 : DO i = atoms%jri(jatom) + 1,ncmsh
104 218 : IF (input%l_core_confpot) THEN
105 216 : rr = d*rr
106 216 : v1 = rr*( t2 + rr*t1 )
107 216 : v2 = rr*( t2b + rr*t1 )
108 216 : vrd(i) = 0.5*(v2 + v1)
109 216 : brd(i) = 0.5*(v2 - v1)
110 : ELSE
111 0 : vrd(i) = vrd(atoms%jri(jatom)) + t2* (i-atoms%jri(jatom))
112 0 : brd(i) = brd(atoms%jri(jatom)) + t2b* (i-atoms%jri(jatom))
113 : ENDIF
114 : END DO
115 : END IF
116 :
117 : ! rr = rmsh(1,jatom)
118 : ! do i =1, ncmsh
119 : ! rr = d*rr
120 : ! write(*,'(3f20.10)') rr,vrd(i),brd(i)
121 : ! enddo
122 :
123 : !
124 2 : rnot = atoms%rmsh(1,jatom)
125 2 : z = atoms%zatom(jatom)
126 2 : dxx = atoms%dx(jatom)
127 :
128 : CALL spratm(atoms%msh,vrd,brd,z,rnot,dxx,ncmsh,&
129 2 : etab(1,jatom),ntab(1,jatom),ltab(1,jatom), sume,rhochr,rhospn)
130 :
131 2 : seig = seig + atoms%neq(jatom)*sume
132 : !
133 : ! rho_up=2(ir) = (rhochr(ir) + rhospn(ir))*0.5
134 : ! rho_dw=1(ir) = (rhochr(ir) - rhospn(ir))*0.5
135 : !
136 2 : IF (input%jspins.EQ.2) THEN
137 1132 : DO j = 1,atoms%jri(jatom)
138 1130 : rhcs(j,jatom,input%jspins) = (rhochr(j)+rhospn(j))*0.5
139 1132 : rhcs(j,jatom,1) = (rhochr(j)-rhospn(j))*0.5
140 : END DO
141 : ELSE
142 0 : DO j = 1,atoms%jri(jatom)
143 0 : rhcs(j,jatom,1) = rhochr(j)
144 : END DO
145 : END IF
146 2 : IF (input%jspins.EQ.2) THEN
147 1348 : DO j = 1,atoms%msh
148 1346 : rhc(j,jatom,input%jspins) = (rhochr(j)+rhospn(j))*0.5
149 1348 : rhc(j,jatom,1) = (rhochr(j)-rhospn(j))*0.5
150 : ENDDO
151 : ELSE
152 0 : DO j = 1,atoms%msh
153 0 : rhc(j,jatom,1) = rhochr(j)
154 : END DO
155 : END IF
156 : !
157 : ! store atomic eigenvalues to file.58
158 2 : IF (jatom.EQ.1) REWIND 58
159 2 : WRITE (58,FMT=*) nkmust(jatom)
160 38 : DO k = 1,nkmust(jatom)
161 38 : WRITE (58,FMT='(f12.6,2i3)') etab(k,jatom),ntab(k,jatom), ltab(k,jatom)
162 : END DO
163 : !---->update spherical charge density rho with the core density.
164 6 : CALL ccdnup(atoms,sphhar,input,jatom, rho, sume,vrs,rhochr,rhospn, tecs,qints)
165 :
166 : END DO ! loop over atoms (jatom)
167 : !
168 : !----> store core charge densities
169 2 : CALL writeCoreDensity(input,atoms,rhcs,tecs,qints)
170 2 : END SUBROUTINE coredr
171 : END MODULE m_coredr
|
