diff --git a/src/COMMON/m_inpmodules.f90 b/src/COMMON/m_inpmodules.f90 index 1d17a41..c5b5771 100644 --- a/src/COMMON/m_inpmodules.f90 +++ b/src/COMMON/m_inpmodules.f90 @@ -546,7 +546,7 @@ module potential_inp ! rovr(novrx,0:nphx) - r for overlap shell ! Added by Fer ! Used to correct the excitation energy for chemical shifts - integer ChSh_Type + integer ChSh_Type integer configtype !KJ 12-2010 : which method for choosing atomic configuration? double precision corval_emin !KJ 12-2012 defines energy window for search for core-valence separation energy. @@ -565,6 +565,7 @@ module potential_inp integer negrid INTEGER iscfth ! TS 07/2020 integer nmu ! cja 10/15/20 + real*8 custom_chsh contains subroutine potential_allocate @@ -611,6 +612,8 @@ subroutine potential_write write(3,20) ChSh_Type write(3,10) 'ConfigType:' write(3,20) configtype + WRITE(3,10) 'ChSh (in eV):' + WRITE(3,*) custom_chsh ! Electronic temperature for thermal scf routine write(3,10) 'Temperature (in eV):' write(3,*) scf_temperature, scf_thermal_vxc @@ -629,8 +632,8 @@ end subroutine potential_write subroutine potential_read integer ititle,ip,iph,iovr - call potential_allocate - call potential_init + call potential_allocate + call potential_init open (file=filename, unit=3, status='old') read(3,*) ; read(3,*) mpot, nph, ntitle, ihole, ipr1, iafolp, ixc, ispec, iscfxc read(3,*) ; read(3,*) nmix, nohole, jumprm, inters, nscmt, icoul, lfms1, iunf @@ -646,12 +649,13 @@ subroutine potential_read read(3,*) ; read(3,*) (novr(iph), iph=0,nph) read(3,*) do iph = 0, nph - do iovr = 1, novr(iph) - read(3,*) iphovr(iovr, iph), nnovr(iovr,iph), rovr(iovr,iph) - enddo - enddo + do iovr = 1, novr(iph) + read(3,*) iphovr(iovr, iph), nnovr(iovr,iph), rovr(iovr,iph) + enddo + enddo read(3,*) ; read(3,*) ChSh_Type read(3,*,end=55) ; read(3,*,end=55) configtype + read(3,*) ; read(3,*) custom_chsh read(3,*) ; read(3,*) scf_temperature, scf_thermal_vxc read(3,*) ; read(3,*) iscfth, xntol, nmu read(3,*) ; read(3,*) negrid, emaxscf diff --git a/src/DEP/xsph.mk b/src/DEP/xsph.mk index 172933e..9885e7a 100644 --- a/src/DEP/xsph.mk +++ b/src/DEP/xsph.mk @@ -57,7 +57,9 @@ xsphSRC = \ ./TDLDA/ellfun.f90 ./TDLDA/yzktd.f90 ./MATH/conv.f90 \ ./TDLDA/kkchi.f90 ./TDLDA/dmscf.f90 ./XSPH/phase.f90 \ ./XSPH/phase_h.f90 ./XSPH/wphase.f90 ./XSPH/wrxsph.f90 \ -./XSPH/axafs.f90 ./MATH/determ.f90 +./XSPH/axafs.f90 ./MATH/determ.f90 \ +./EXCH/rgw.f90 ./EXCH/imgw.f90 + xsph_MODULESRC = \ ./PAR/m_par.f90 ./COMMON/m_constants.f90 ./KSPACE/m_struct.f90 \ ./KSPACE/m_controls.f90 ./KSPACE/m_boundaries.f90 ./KSPACE/m_kklist.f90 \ @@ -73,3 +75,4 @@ xsph_MODULESRC = \ ./IOMODS/m_iomod.f90 ./COMMON/m_stkets.f90 ./FMS/m_fms.f90 \ ./COMMON/m_t3j.f90 ./SELF/m_SelfEnergy.f90 ./TDLDA/getchi0.f90 \ ./COMMON/m_kinds.f90 ./POT/m_mtdp.f90 ./POT/m_atomicpotio.f90 + diff --git a/src/EXCH/imgw.f90 b/src/EXCH/imgw.f90 new file mode 100644 index 0000000..83c7d1a --- /dev/null +++ b/src/EXCH/imgw.f90 @@ -0,0 +1,643 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Information about last revision of $RCSfile: rhl.f90,v $: +! $Revision: 1.4 $ +! $Author: hebhop $ +! $Date: 2010/02/23 23:52:06 $ +! +! v1 - Mixed weights, lost track of origin +! v2 - Using local training weights - scipy 1.7.1 +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +SUBROUTINE imgw_v1(rs, xk, temp, eim) + ! temperature in unit of hartree + ! xk in unit of a.u. + ! eim in unit of hatree + USE constants + IMPLICIT NONE + REAL(8), INTENT(IN) :: rs, xk, temp + REAL(8), INTENT(OUT) :: eim + + INTEGER :: mrs + REAL(8), DIMENSION(:) :: rsleft(24), tright(6), rsright(6), scaling(4) + REAL(8), DIMENSION(27) :: left1, left2, left3, left4, p + REAL(8), DIMENSION(36) :: right1, right2, right3, right4, qq + REAL(8) :: vleft, vright, weight, k_boundary, t, k, rkf, ef, eim_kf + ! REAL(8), PARAMETER :: fa = 1.9191582926775128d0, pi = 3.1415926535897932384626433832795d0 + + DATA scaling / 0.012163698621989414d0, 0.06690755186151459d0, 0.35402213905982843d0, 0.6540650734032093d0 / + DATA left1 / -0.09730116611517951d0, 0.019811787581018432d0, -0.9107593983619665d0, & + 0.04429457958312912d0, -0.002829886983994427d0, 0.4679177248651782d0, & + -0.12920151858760917d0, -0.3185960944710708d0, -0.6853260400838367d0, & + 0.05864888772460627d0, 0.19062026200270787d0, 0.3320705055412537d0, & + 0.22569358685714208d0, 0.24878999466001417d0, 1.6620823965743552d0, & + -0.10285278556823735d0, -0.1581723331764946d0, -0.838851955728067d0, & + 0.021481438539048946d0, 0.18375504201186213d0, 0.2646365593837374d0, & + 0.028118057019198332d0, 0.3878911095790498d0, 0.12429020135523595d0, & + -0.04923156950114449d0, -0.5504701432426441d0, -0.4160785576258503d0 / + + DATA right1 / -1.1281702352664078d0, 15.913917121155512d0, -15.781763286731191d0, & + -8.196227222731984d0, 49.57397322108513d0, -54.859186787268406d0, & + 7.148275167140558d0, -37.08415785218802d0, 41.63729110321263d0, & + -12.218014372748724d0, 59.87283626607546d0, -67.71261977276231d0, & + -17.246302998638555d0, 82.30740207451787d0, -94.09633980972686d0, & + 2.1606534565867697d0, 5.411632668505843d0, -12.179092411480768d0, & + 17.994004196294703d0, 9.123347236345582d0, 8.606772473958118d0, & + -58.47312852780925d0, 1.0625113379641786d0, -48.71757053540676d0, & + -355.91252983782226d0, 190.63518944485625d0, -310.49378130241234d0, & + 214.80355604369302d0, -131.7421202745962d0, 118.66215483704421d0, & + -429.60939281012014d0, 286.4172905542458d0, -315.45204472938093d0, & + -1.4907763646708778d0, 8.076166718277117d0, -1.614001151428006d0 / + + DATA left2 / -0.1601209222500568d0, -1.1689316143699513d0, 1.0558126886303d0, & + 0.11309902121772979d0, 0.5849459844670819d0, -0.5726049588912363d0, & + -0.36416056624418897d0, -0.86890695165242d0, 0.8088218974342717d0, & + 0.23222343055157713d0, 0.38696774972725145d0, -0.41242408747826753d0, & + 0.5499117076571199d0, 1.8790287734597042d0, -1.730582839999667d0, & + -0.36043076429733906d0, -0.8785085669745413d0, 0.9058703916880849d0, & + -0.014018914957450708d0, 0.6041754397450181d0, -0.3781817796528243d0, & + 0.015415318386368074d0, 0.7431802071380266d0, -0.47248889385247694d0, & + 0.002648081567479078d0, -1.3461391390947097d0, 0.8485070219155748d0 / + DATA right2 / 1.2124314405967718d0, 4.362653904918381d0, 0.0644870840259468d0, & + -6.052826011729766d0, -14.582800518248705d0, -3.0475157298312276d0, & + 4.891627448962329d0, 15.535535274868563d0, 8.71073073332997d0, & + 1.6916902132193923d0, 6.129100021914905d0, 12.986683851023564d0, & + -4.795324539850384d0, 9.780199586281656d0, -68.77445618921504d0, & + 2.1223240396942145d0, -0.9853858615672741d0, 0.07285272264221357d0, & + 3.040897422460228d0, -1.624057453264639d0, 1.2388597307831737d0, & + 3.0808972467314835d0, -3.2825285078006523d0, 1.8322531925879286d0, & + 7.538853613459758d0, -10.689962213288286d0, 4.685103645264265d0, & + -3.149847495870372d0, 5.846712789516093d0, -2.6993422306528916d0, & + -16.785211916080595d0, 26.283869943955267d0, -11.124802704730993d0, & + -0.25041311445387615d0, 2.9450940007931896d0, -0.28594019092691975d0 / + DATA left3 / -0.4888688812237603d0, 0.28539850046718035d0, -0.08273287486482286d0, & + 0.2656325902645167d0, -0.1849691784345592d0, 0.052504778245743705d0, & + -0.7222924620244626d0, 0.3576236963517782d0, -0.07483691492465128d0, & + 0.3952865822390625d0, -0.22927563599528483d0, 0.0494082188549157d0, & + 1.2045675641067852d0, -0.6302953954558612d0, 0.1535481192430782d0, & + -0.6593883225209146d0, 0.4092545459441451d0, -0.10005795779991933d0, & + 0.2187487514974954d0, -0.003790727935850534d0, 0.00153321067844436d0, & + 0.3222868901764764d0, -0.03137964574105422d0, 0.0009654443979920085d0, & + -0.5338036685775893d0, 0.03076617000441788d0, -0.002236088599183545d0 / + DATA right3 / -2.4652385454861045d0, 2.2709369946116476d0, -0.016335383516986d0, & + 6.407346282173194d0, -3.5241610699720565d0, -0.938052150275685d0, & + -5.330201111232691d0, 0.2097150467245504d0, 2.3913144246314464d0, & + -2.759643208969734d0, -4.6749207208828905d0, 4.846048477877578d0, & + 1.274936238844436d0, -8.99534035230194d0, 6.344175089142765d0, & + 1.5942536840899906d0, -0.7080734263680196d0, 0.06020037998104392d0, & + 1.9858949542281634d0, -0.2126674702243764d0, 0.4372020277869658d0, & + 3.0259804599009854d0, -0.1865921940484292d0, 1.2872381342191748d0, & + 6.2597474728625455d0, -0.6447257049655256d0, 5.661178863033677d0, & + 0.9511711716444875d0, -5.071556705823834d0, 11.077219657751685d0, & + -5.235149214896328d0, 2.958138301791755d0, -11.085413113349038d0, & + 0.6646470406851569d0, 0.40897370602960853d0, 0.18063699128754584d0 / + + DATA left4 / -0.2784040622721631d0, -0.015923208658958306d0, 0.010936406080942379d0, & + 0.15177687739002343d0, -0.013671299652304182d0, -0.001955551467937442d0, & + -0.6602235069910223d0, 0.22243064479560798d0, -0.026186885954275728d0, & + 0.3507468887735972d0, -0.14026926391698746d0, 0.018100444392833998d0, & + 0.9916213260703696d0, -0.2500523604133429d0, 0.024459854621137585d0, & + -0.5327649552961942d0, 0.17914527978686218d0, -0.02142230136888505d0, & + 0.2837717099371378d0, -0.04225811704879709d0, 0.006162296703193555d0, & + 0.3957333770673754d0, -0.08131130660031956d0, 0.009023081733350145d0, & + -0.6749994220174723d0, 0.12072704401870088d0, -0.015089642050132496d0 / + DATA right4 / -6.000006016596565d0, -0.5674645360422231d0, -0.08194210870224415d0, & + -20.09957970932959d0, -7.6736716302976475d0, -0.9082009774103044d0, & + -9.721151171770373d0, -10.391407207511312d0, -3.7162057539502245d0, & + -3.375549352685843d0, 1.6362545698898001d0, -3.2904412671133585d0, & + 7.399267039799666d0, -8.563505107372434d0, 4.015324424866241d0, & + -1.0294921526228644d0, 1.1312226646636925d0, -0.05235612645649135d0, & + -0.1766141100011096d0, -0.11009354980440053d0, -0.12319707443052594d0, & + 0.008069772186224945d0, 0.1689695741639606d0, 0.15614385532021963d0, & + 0.272368412409354d0, -0.5726722195398571d0, -0.5296986853892816d0, & + -4.436330713599893d0, 3.7741343244695202d0, 6.344404987856783d0, & + -2.2339091003314797d0, 3.2818227500365795d0, 3.7788440868255257d0, & + 0.2676830746051518d0, -0.10237475806456857d0, -0.020508952084653007d0 / + + !------------------------------------------------------------ + ! Preprocessing + !------------------------------------------------------------ + rkf = fa/rs + ef = (rkf**2.d0)/2.d0 + k = xk / rkf + t = temp / ef + + !------------------------------------------------------------ + ! Check for interpolation or extrapolation + !------------------------------------------------------------ + IF (t.GT.2.d0) THEN + PRINT*, "ERROR: Local electronic temperature is higher interpolation grid (T >= 2 TF)." + PRINT*, "T = ", t, temp + PRINT*, "rs = ", rs + CALL par_stop + ENDIF + + IF (rs.GT.20.d0) THEN + PRINT*, "ERROR: rs > 20" + CALL par_stop + ENDIF + + IF (rs.LT.1.d-2) THEN + PRINT*, "ERROR: rs < 0.01" + CALL par_stop + ENDIF + + ! Set to fermi level if below fermi level + IF (k.LT.1.00001d0) THEN + k = 1.00001d0 + ENDIF + + ! Determine which region of rs to use + IF (rs.LT.0.2) THEN + p = left1 + qq = right1 + mrs = 1 + ELSEIF (rs.LT.1.0) THEN + p = left2 + qq = right2 + mrs = 2 + ELSEIF (rs.LT.5.0) THEN + p = left3 + qq = right3 + mrs = 3 + ELSE + p = left4 + qq = right4 + mrs = 4 + ENDIF + + !-------------------- + ! Imaginary part + !-------------------- + CALL im_cusp_v1(rs, t, k_boundary) + CALL im_kf_v1(rs, t, eim_kf) + + ! Left hand side + rsleft(1) = p(1) * rs + p(2)*rs**(3.d0/2.d0) + p(3)*rs**2.d0 + rsleft(2) = p(4) * rs + p(5)*rs**(3.d0/2.d0) + p(6)*rs**2.d0 + rsleft(3) = p(7) * rs + p(8)*rs**(3.d0/2.d0) + p(9)*rs**2.d0 + rsleft(4) = p(10) * rs + p(11)*rs**(3.d0/2.d0) + p(12)*rs**2.d0 + rsleft(5) = p(13) * rs + p(14)*rs**(3.d0/2.d0) + p(15)*rs**2.d0 + rsleft(6) = p(16) * rs + p(17)*rs**(3.d0/2.d0) + p(18)*rs**2.d0 + rsleft(7) = p(19) * rs + p(20)*rs**(3.d0/2.d0) + p(21)*rs**2.d0 + rsleft(8) = p(22) * rs + p(23)*rs**(3.d0/2.d0) + p(24)*rs**2.d0 + rsleft(9) = p(25) * rs + p(26)*rs**(3.d0/2.d0) + p(27)*rs**2.d0 + + vleft = rsleft(1)*k*t & + + rsleft(2)*k*(t**1.5d0) & + + rsleft(3)*(k**2.d0)*t & + + rsleft(4)*(k**2.d0)*(t**1.5d0) & + + rsleft(5)*(k**1.5d0)*t & + + rsleft(6)*(k**1.5d0)*(t**1.5d0) & + + rsleft(7)*k & + + rsleft(8)*(k**2.d0) & + + rsleft(9)*(k**1.5d0) & + + eim_kf + IF (vleft.LT.0.d0) THEN + vleft = 0.d0 + ENDIF + ! print*, rsleft(1), rsleft(2), rsleft(3), rsleft(4), rsleft(5), rsleft(6) + ! print*, rsleft(7), rsleft(8), rsleft(9) + + + ![-0.46372935] [0.20346769] [0.39460677] [-0.68524943] [0.32048492] [0.51356274] + ![1.15491929] [-0.52727977] [-0.90495454] + + + ! Right hand side + rsright(1) = qq(1) * rs + qq(2)*rs**(3.d0/2.d0) + qq(3)*rs**2.d0 + rsright(2) = qq(4) * rs + qq(5)*rs**(3.d0/2.d0) + qq(6)*rs**2.d0 + rsright(3) = qq(7) * rs + qq(8)*rs**(3.d0/2.d0) + qq(9)*rs**2.d0 + rsright(4) = qq(10) * rs + qq(11)*rs**(3.d0/2.d0) + qq(12)*rs**2.d0 + rsright(5) = qq(13) * rs + qq(14)*rs**(3.d0/2.d0) + qq(15)*rs**2.d0 + rsright(6) = qq(16) * rs + qq(17)*rs**(3.d0/2.d0) + qq(18)*rs**2.d0 + + tright(1) = qq(19) + qq(20)*t + qq(21)*t**2.d0 + tright(2) = qq(22) + qq(23)*t + qq(24)*t**2.d0 + tright(3) = qq(25) + qq(26)*t + qq(27)*t**2.d0 + tright(5) = qq(28) + qq(29)*t + qq(30)*t**2.d0 + tright(4) = qq(31) + qq(32)*t + qq(33)*t**2.d0 + tright(6) = qq(34) + qq(35)*t + qq(36)*t**2.d0 + ! print*, tright(1), tright(2), tright(3), tright(4), tright(5), tright(6) + ! print*, rsright(1), rsright(2), rsright(3), rsright(4), rsright(5), rsright(6) + + ! [1.99086217] [3.03543934] [6.25979544] [-5.23141634] [0.9495449] [0.66449063] tright + ! [1.03885912] [-0.11115822] [-1.25343129] [-0.43254481] [1.26940609] [1.36394528] rsright + ! [3.30936812] [7.80174861] [27.61501152] [-43.66052506] [35.11693639] [2.20514242] + ! [1.04151751] [-0.11174835] [-1.25312891] [-0.43187571] [1.26780354] [1.36359809] + + vright = rsright(1)*tright(1)/k & + + rsright(2)*tright(2)/(k**2.d0) & + + rsright(3)*tright(3)/(k**3.d0) & + + rsright(4)*tright(4)/(k**4.d0) & + + rsright(5)*tright(5)/(k**5.d0) & + + rsright(6)*tright(6)/(k**0.5d0) + IF (vright.LT.0.d0) THEN + vright = 0.d0 + ENDIF + ! Interpolate between regions + weight = 1.d0/(1.d0 + EXP(10.d0*(k - 1.1d0*k_boundary))) + eim = weight*vleft + (1.d0-weight)*vright*scaling(mrs) + + ! print*, "Reduced temperature = ", t, " rs = ", rs + ! Scale back + eim = -eim * ef + RETURN +END SUBROUTINE imgw_v1 + +SUBROUTINE im_kf_v1(rs, t, kf) + ! Imaginary part of quasiparticle self-energy at fermi level + IMPLICIT NONE + REAL(8), INTENT(IN) :: rs, t + REAL(8), INTENT(OUT) :: kf + REAL(8), DIMENSION(9) :: factors1, factors2, factors3, factors4, factors + REAL(8) :: A, B, C + INTEGER :: mrs + DATA factors1 / -0.09094323574018673d0, 1.089160117026636d0, -0.7292298662705564d0, & + 0.006462381104206332d0, -0.019802206542834953d0, 0.01113010261742388d0, & + 0.04903214800885628d0, -0.6506372471873646d0, 0.5274238755315628d0 / + DATA factors2 / 0.4316369508779085d0, -0.17245974727080232d0, 0.009809795774884911d0, & + 3.5768469939869205d-5, -0.004287037992540818d0, 0.0020415490627589093d0, & + -0.47317800570803087d0, 0.610089125206505d0, -0.211092319733901d0 / + DATA factors3 / -0.009515617697671723d0, 0.09982092613155993d0, 0.07964035978307292d0, & + -0.030979381407273337d0, 0.41236562921436426d0, -0.33077114773922256d0, & + 0.00322951237168846d0, -0.03565817188719173d0, 0.004253783508400203d0 / + DATA factors4 / -0.2759688536430732d0, 0.7430961585847626d0, -0.297181644446135d0, & + 0.26653826932442615d0, -0.30590572296821034d0, 0.08998256117324012d0, & + 0.04277863363447566d0, -0.13113824718924794d0, 0.0601847393212351d0 / + + ! 00 -> 1 + ! 01 -> 2 + ! 10 -> 3 + ! 11 -> 4 + mrs = 1 + IF (rs.LE.0.2d0) THEN + mrs = mrs + 0 + ELSE + mrs = mrs + 2 + ENDIF + + IF (t.LE.0.5d0) THEN + mrs = mrs + 0 + ELSE + mrs = mrs + 1 + ENDIF + + ! Shift from 0-index to 1-index + SELECT CASE (mrs) + CASE (1) + factors = factors1 + CASE (2) + factors = factors2 + CASE (3) + factors = factors3 + CASE (4) + factors = factors4 + END SELECT + + A = factors(1)*t + factors(2)*t**1.5d0 + factors(3)*t**2.d0 + B = factors(4)*t + factors(5)*t**1.5d0 + factors(6)*t**2.d0 + C = factors(7)*t + factors(8)*t**1.5d0 + factors(9)*t**2.d0 + kf = A*rs + B*rs**0.5d0 + C*rs**1.5d0 + RETURN +END SUBROUTINE im_kf_v1 + +SUBROUTINE im_cusp_v1(rs, t, k) + ! Position of minimum for imaginary quasiparticle self-energy + ! + ! Inputs: + ! rs - radius + ! t - temperature in unit of fermi + ! Output: + ! k - temperature in unit of fermi + ! + IMPLICIT NONE + REAL(8), INTENT(IN) :: rs, t + REAL(8), INTENT(OUT) :: k + REAL(8), PARAMETER :: alpha = 1.1495192084389452d0 + REAL(8), PARAMETER :: beta = 0.033376683130363946d0 + REAL(8), DIMENSION(6) :: factors + REAL(8) :: A, B + DATA factors / 0.1450147977529337d0, 0.233244548502498d0, -0.017661811059191113d0, & + 1.3723845833398063d0, -0.05395611403472079d0, -0.02022335357291743d0 / + + A = factors(1) + factors(2)*t + factors(3)*t**2.d0 + B = factors(4) + factors(5)*t + factors(6)*t**2.d0 + k = (1.d0 + TANH((rs - 0.5d0)/alpha))*(beta*rs + A) + B + + RETURN +END SUBROUTINE im_cusp_v1 + +SUBROUTINE imgw_v2(rs, xk, temp, eim) + ! temperature in unit of hartree + ! xk in unit of a.u. + ! eim in unit of hatree + ! USE constants + IMPLICIT NONE + REAL(8), INTENT(IN) :: rs, xk, temp + REAL(8), INTENT(OUT) :: eim + + INTEGER :: mrs + REAL(8), DIMENSION(:) :: rsleft(24), tright(6), rsright(6), scaling(4) + REAL(8), DIMENSION(27) :: left1, left2, left3, left4, p + REAL(8), DIMENSION(36) :: right1, right2, right3, right4, qq + REAL(8) :: vleft, vright, weight, k_boundary, t, k, rkf, ef, eim_kf + REAL(8), PARAMETER :: fa = 1.9191582926775128d0, pi = 3.1415926535897932384626433832795d0 + + DATA scaling / 0.011358051321446327, 0.0637710259641874, 0.3458217538503771, 0.6459817122849331 / + DATA left1 / -0.3996084877917776d0, 1.5619375484723628d0, -2.8394175870705833d0, & + 0.23209925292414174d0, -0.9654187419185194d0, 1.6712478624220377d0, & + -0.4874561959141273d0, 1.5037256449081757d0, -2.9596237738238536d0, & + 0.2809690251217274d0, -0.9428606121254859d0, 1.7436729543186533d0, & + 0.8875815885158318d0, -3.122410046101803d0, 5.873364460085901d0, & + -0.5136592905099191d0, 1.9416086553623595d0, -3.4584841020954467d0, & + 0.10522063267277822d0, -0.23502107642632086d0, 0.7931114188357979d0, & + 0.12741075179508662d0, -0.11312575631205433d0, 0.7600031082923175d0, & + -0.232860065381581d0, 0.3717152558382517d0, -1.582843856594346d0 / + + DATA right1 / -1.051699097983817d0, 15.7605480356678d0, -15.233090644952417d0, & + 3.578781012760475d0, -21.73202416174778d0, 23.81498829488864d0, & + -7.88089609855281d0, 40.88236714341268d0, -45.525707747104185d0, & + -30.002393040691693d0, 146.7394852005877d0, -164.67171231292798d0, & + -15.605129680157468d0, 74.23528764158678d0, -84.30470772187726d0, & + 2.062510163377484d0, 6.089922176277204d0, -13.268159816006849d0, & + 17.321743110637776d0, 10.192118590735618d0, 8.7641751541308d0, & + 128.82200744826002d0, 4.464680955918343d0, 112.21833288947664d0, & + 313.0137985207077d0, -156.72890989858854d0, 277.7869396298057d0, & + -168.9585087316882d0, 105.49447512610465d0, -122.59474478201288d0, & + 227.52417460894688d0, -126.2597384129594d0, 118.75117989962092d0, & + -1.5051768093475526d0, 8.562414898102029d0, -1.693502453897044d0 / + + DATA left2 / -0.09165761253230811d0, -1.0100877689661594d0, 0.7905437610016559d0, & + 0.08529550258124016d0, 0.46021875472047963d0, -0.3980814806818628d0, & + -0.13887859413506645d0, -1.2985935331860234d0, 0.9885016763968162d0, & + 0.11079913844671392d0, 0.6123501259132299d0, -0.5028960917095654d0, & + 0.22234507917902402d0, 2.291507964580495d0, -1.756313542111056d0, & + -0.1918368990456516d0, -1.0612123498169421d0, 0.8865352374221388d0, & + -0.09934055662710367d0, 0.826682845998962d0, -0.5077049296467399d0, & + -0.08334894751470627d0, 1.002357219687475d0, -0.6237860934520787d0, & + 0.18796116873132362d0, -1.831089556300899d0, 1.1312699585795796d0 / + DATA right2 / 1.0111657040642084d0, 4.064326038687206d0, 0.16846258022713467d0, & + -3.5506429566055306d0, -7.990788334751925d0, -2.182518930294372d0, & + 6.246803206576854d0, 12.98741573261871d0, 10.367308772292617d0, & + 5.636778874087016d0, 10.690116333976434d0, 30.739773671316126d0, & + -0.7190194276269423d0, 0.9545234353855702d0, -14.129097658240307d0, & + 2.102500838214577d0, -0.9431144829469508d0, 0.05131256788772643d0, & + 3.314841881874237d0, -1.7817457593618435d0, 1.3493106467838916d0, & + 5.228419733993317d0, -5.588526786097672d0, 3.09535168989133d0, & + 7.031182163777919d0, -10.012834657452048d0, 4.2963979842778866d0, & + -6.7570100154263395d0, 10.639863170117842d0, -4.431157054856769d0, & + -12.496869829024636d0, 23.95309512752951d0, -11.191038511186278d0, & + -0.25511793974662095d0, 3.105525998799246d0, -0.2700075247878443d0 / + DATA left3 / -0.26704300240729906d0, -0.021244137631247734d0, 0.005698796489891547d0, & + 0.11283604277220974d0, 0.019045454633044828d0, -0.0054277807029545185d0, & + -0.5469201480659264d0, 0.1390065934386878d0, -0.01601625505901977d0, & + 0.2708598173868438d0, -0.07878943179251854d0, 0.009344615815916932d0, & + 0.8111313207199795d0, -0.12157599964456617d0, 0.013189802361680087d0, & + -0.3834160670049347d0, 0.06336570626257887d0, -0.005829363923190015d0, & + 0.2125698651055652d0, 0.007932536482860863d0, -0.00124853550536081d0, & + 0.3193597622536932d0, -0.02167195555951595d0, -0.0015187500633450784d0, & + -0.5247600432423136d0, 0.009563140232016603d0, 0.0029171635166231604d0 / + DATA right3 / -0.3150344671029144d0, 7.363395001078842d0, 0.18630501915445655d0, & + -0.9576430005524192d0, -2.8271106778713824d0, -0.8599630632534063d0, & + -2.2992847874530837d0, -3.534208989980021d0, -2.319658720592565d0, & + 2.115642134468694d0, 6.272665143976334d0, 5.974450461799069d0, & + 0.5888717390998074d0, -1.7755382467572316d0, -1.0209269767285236d0, & + 1.4707685845122445d0, -0.24926081938255257d0, 0.0306479062551625d0, & + 0.347339469478298d0, -0.18901400041978014d0, 0.044021750409341934d0, & + 2.2022033947039525d0, -1.6328495350124697d0, 0.12767043506549883d0, & + -3.465557297038497d0, 2.8964610624546245d0, 0.531207342268555d0, & + -2.7075622949220635d0, 2.1708189359224606d0, 0.9901363088912304d0, & + -7.6563481562115285d0, 7.7215240277806485d0, 5.1441415569288385d0, & + 0.08361335594468824d0, 0.5330486534456583d0, 0.08731307883881419d0 / + + DATA left4 / -0.7255764908622497d0, 0.3093792864396717d0, -0.04972994251238999d0, & + 0.43196107390671407d0, -0.21694226100958927d0, 0.03587075561538992d0, & + -0.8331809868322648d0, 0.33645009481627325d0, -0.046566693017091876d0, & + 0.46642805756248357d0, -0.2174406467324864d0, 0.03192758597637695d0, & + 1.5214282268545376d0, -0.6182276335453512d0, 0.09193920949179708d0, & + -0.8769564937958089d0, 0.4187362571018319d0, -0.06527303630214303d0, & + 0.2591505164264935d0, -0.018219696625697477d0, 0.0014599894833166724d0, & + 0.37702649496173785d0, -0.06106371579984695d0, 0.0048997611081584d0, & + -0.6284002621066553d0, 0.073893129187332d0, -0.00579591281241008d0 / + DATA right4 / 3.8526675952729166d0, -0.5569276374977795d0, 0.04528019859397363d0, & + 17.166502327516778d0, -0.45002686665649694d0, -0.18820682694326124d0, & + -24.462015831364145d0, -6.159224437087954d0, 1.8399672551799358d0, & + -6.261590646081136d0, -11.807307405220511d0, 2.748734636468226d0, & + -3.7126314337964885d0, 14.060405013130548d0, -2.6174447194118313d0, & + -1.1921369427480966d0, 1.3239339182810987d0, -0.05785069917390398d0, & + 0.2647118302555023d0, 0.20610301216238663d0, 0.2333740612693166d0, & + 0.3617041524382783d0, -0.285361076924459d0, -0.23845121580912843d0, & + 1.5877270765978815d0, -0.5343997363603614d0, -0.341836713115954d0, & + -6.038348097520621d0, 1.0453428353642356d0, 0.6260596337936863d0, & + -7.133728728079829d0, 0.11396959493126725d0, 0.04687586216904349d0, & + 0.289235549542202d0, -0.06477313647733891d0, 0.010361014288148208d0 / + + !------------------------------------------------------------ + ! Preprocessing + !------------------------------------------------------------ + rkf = fa/rs + ef = (rkf**2.d0)/2.d0 + k = xk / rkf + t = temp / ef + + !------------------------------------------------------------ + ! Check for interpolation or extrapolation + !------------------------------------------------------------ + IF (t.GT.2.d0) THEN + PRINT*, "ERROR: Local electronic temperature is higher interpolation grid (T >= 2 TF)." + PRINT*, "T = ", t, temp + PRINT*, "rs = ", rs + ! CALL par_stop + ENDIF + + IF (rs.GT.20.d0) THEN + PRINT*, "ERROR: rs > 20 is outside of fit range" + CALL par_stop + ENDIF + + IF (rs.LT.1.d-2) THEN + PRINT*, "ERROR: rs < 0.01 is outside of fit range" + CALL par_stop + ENDIF + + ! Set to fermi level if below fermi level + IF (k.LT.1.00001d0) THEN + k = 1.00001d0 + ENDIF + + ! Determine which region of rs to use + IF (rs.LT.0.2) THEN + p = left1 + qq = right1 + mrs = 1 + ELSEIF (rs.LT.1.0) THEN + p = left2 + qq = right2 + mrs = 2 + ELSEIF (rs.LT.5.0) THEN + p = left3 + qq = right3 + mrs = 3 + ELSE + p = left4 + qq = right4 + mrs = 4 + ENDIF + + !-------------------- + ! Imaginary part + !-------------------- + CALL im_cusp_v2(rs, t, k_boundary) + CALL im_kf_v2(rs, t, eim_kf) + + ! Left hand side + rsleft(1) = p(1) * rs + p(2)*rs**(3.d0/2.d0) + p(3)*rs**2.d0 + rsleft(2) = p(4) * rs + p(5)*rs**(3.d0/2.d0) + p(6)*rs**2.d0 + rsleft(3) = p(7) * rs + p(8)*rs**(3.d0/2.d0) + p(9)*rs**2.d0 + rsleft(4) = p(10) * rs + p(11)*rs**(3.d0/2.d0) + p(12)*rs**2.d0 + rsleft(5) = p(13) * rs + p(14)*rs**(3.d0/2.d0) + p(15)*rs**2.d0 + rsleft(6) = p(16) * rs + p(17)*rs**(3.d0/2.d0) + p(18)*rs**2.d0 + rsleft(7) = p(19) * rs + p(20)*rs**(3.d0/2.d0) + p(21)*rs**2.d0 + rsleft(8) = p(22) * rs + p(23)*rs**(3.d0/2.d0) + p(24)*rs**2.d0 + rsleft(9) = p(25) * rs + p(26)*rs**(3.d0/2.d0) + p(27)*rs**2.d0 + + vleft = rsleft(1)*k*t & + + rsleft(2)*k*(t**1.5d0) & + + rsleft(3)*(k**2.d0)*t & + + rsleft(4)*(k**2.d0)*(t**1.5d0) & + + rsleft(5)*(k**1.5d0)*t & + + rsleft(6)*(k**1.5d0)*(t**1.5d0) & + + rsleft(7)*k & + + rsleft(8)*(k**2.d0) & + + rsleft(9)*(k**1.5d0) & + + eim_kf + + IF (vleft.LT.0.d0) THEN + vleft = 0.d0 + ENDIF + + ! Right hand side + rsright(1) = qq(1) * rs + qq(2)*rs**(3.d0/2.d0) + qq(3)*rs**2.d0 + rsright(2) = qq(4) * rs + qq(5)*rs**(3.d0/2.d0) + qq(6)*rs**2.d0 + rsright(3) = qq(7) * rs + qq(8)*rs**(3.d0/2.d0) + qq(9)*rs**2.d0 + rsright(4) = qq(10) * rs + qq(11)*rs**(3.d0/2.d0) + qq(12)*rs**2.d0 + rsright(5) = qq(13) * rs + qq(14)*rs**(3.d0/2.d0) + qq(15)*rs**2.d0 + rsright(6) = qq(16) * rs + qq(17)*rs**(3.d0/2.d0) + qq(18)*rs**2.d0 + + tright(1) = qq(19) + qq(20)*t + qq(21)*t**2.d0 + tright(2) = qq(22) + qq(23)*t + qq(24)*t**2.d0 + tright(3) = qq(25) + qq(26)*t + qq(27)*t**2.d0 + tright(4) = qq(28) + qq(29)*t + qq(30)*t**2.d0 + tright(5) = qq(31) + qq(32)*t + qq(33)*t**2.d0 + tright(6) = qq(34) + qq(35)*t + qq(36)*t**2.d0 + + vright = rsright(1)*tright(1)/k & + + rsright(2)*tright(2)/(k**2.d0) & + + rsright(3)*tright(3)/(k**3.d0) & + + rsright(4)*tright(4)/(k**4.d0) & + + rsright(5)*tright(5)/(k**5.d0) & + + rsright(6)*tright(6)/SQRT(k) + + IF (vright.LT.0.d0) THEN + vright = 0.d0 + ENDIF + + ! Interpolate between regions + weight = 1.d0/(1.d0 + EXP(15.d0*(k - 1.1d0*k_boundary))) + eim = weight*vleft + (1.d0-weight)*vright*scaling(mrs) + + ! Scale back + eim = -eim * ef + RETURN +END SUBROUTINE imgw_v2 + +SUBROUTINE im_kf_v2(rs, t, kf) + ! Imaginary part of quasiparticle self-energy at fermi level + IMPLICIT NONE + REAL(8), INTENT(IN) :: rs, t + REAL(8), INTENT(OUT) :: kf + REAL(8), DIMENSION(12) :: factors1, factors2, factors3, factors4, factors + REAL(8) :: A, B, C, D + INTEGER :: mrs + DATA factors1 / 0.006126619079058038d0, -0.022846382754034916d0, 0.015045114528861647d0, & + -0.14013765402506528d0, 1.732411629425611d0, -1.527915604967462d0, & + 0.2041869214818972d0, -2.6592424440301077d0, 3.0210152215271107d0, & + -0.159621537197849d0, 2.269255518220703d0, -2.8208884235019287d0 / + DATA factors2 / -0.00383652380506507d0, 0.002372632021212689d0, -0.0006936859623145972d0, & + 0.477670622519175d0, -0.2516296967247694d0, 0.04232619703799145d0, & + -0.6070288526414329d0, 0.8402894704206122d0, -0.3056394124643722d0, & + 0.25197834782489975d0, -0.43335938835885957d0, 0.17798795242083698d0 / + DATA factors3 / -0.02526320970813707d0, 0.35287667298530234d0, -0.340763616688629d0, & + 0.0035392012435348998d0, -0.05837406966292795d0, 0.417344190284107d0, & + -0.00024626463664617225d0, 0.00753983071122524d0, -0.1067418540332355d0, & + -1.4336855836569554e-06, -0.00010468083863930219d0, 0.0022898442682024144d0 / + DATA factors4 / 0.23639874118322862d0, -0.3495821472113405d0, 0.12933931275338106d0, & + -0.22490690598335195d0, 0.8170914266911294d0, -0.36385883500314564d0, & + 0.022055208809747342d0, -0.16116895879846668d0, 0.08724548039489805d0, & + 0.0007796226273447165d0, 0.0011297583383124885d0, -0.001018029772211792d0/ + + ! 00 -> 1 + ! 01 -> 2 + ! 10 -> 3 + ! 11 -> 4 + mrs = 1 + IF (rs.LE.0.2d0) THEN + mrs = mrs + 0 + ELSE + mrs = mrs + 2 + ENDIF + + IF (t.LE.0.5d0) THEN + mrs = mrs + 0 + ELSE + mrs = mrs + 1 + ENDIF + + ! Shift from 0-index to 1-index + SELECT CASE (mrs) + CASE (1) + factors = factors1 + CASE (2) + factors = factors2 + CASE (3) + factors = factors3 + CASE (4) + factors = factors4 + END SELECT + + A = factors(1)*t + factors(2)*t**1.5d0 + factors(3)*t**2.d0 + B = factors(4)*t + factors(5)*t**1.5d0 + factors(6)*t**2.d0 + C = factors(7)*t + factors(8)*t**1.5d0 + factors(9)*t**2.d0 + D = factors(10)*t + factors(11)*t**1.5d0 + factors(12)*t**2.d0 + kf = A*rs**0.5d0 + B*rs + C*rs**1.5d0 + D*rs**2.5d0 + IF (kf.LT.0.d0) THEN + kf = 0.d0 + ENDIF + RETURN +END SUBROUTINE im_kf_v2 + +SUBROUTINE im_cusp_v2(rs, t, k) + ! Position of minimum for imaginary quasiparticle self-energy + ! + ! Inputs: + ! rs - radius + ! t - temperature in unit of fermi + ! Output: + ! k - temperature in unit of fermi + ! + IMPLICIT NONE + REAL(8), INTENT(IN) :: rs, t + REAL(8), INTENT(OUT) :: k + REAL(8), PARAMETER :: alpha = 1.24469774d0 + REAL(8), PARAMETER :: beta = 0.03303176d0 + REAL(8), DIMENSION(6) :: factors + REAL(8) :: A, B + DATA factors / 0.15104746d0, 0.23333948d0, -0.01594679d0, 1.36174536d0, -0.05237112d0, -0.02221209d0 / + + A = factors(1) + factors(2)*t + factors(3)*t**2.d0 + B = factors(4) + factors(5)*t + factors(6)*t**2.d0 + k = (1.d0 + TANH((rs - 0.5d0)/alpha))*(beta*rs + A) + B + + RETURN +END SUBROUTINE im_cusp_v2 + \ No newline at end of file diff --git a/src/EXCH/rgw.f90 b/src/EXCH/rgw.f90 new file mode 100644 index 0000000..08c1af8 --- /dev/null +++ b/src/EXCH/rgw.f90 @@ -0,0 +1,242 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! Information about last revision of $RCSfile: rhl.f90,v $: +! $Revision: 1.4 $ +! $Author: hebhop $ +! $Date: 2010/02/23 23:52:06 $ +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +SUBROUTINE rgw(rs, xk, temp, erl) + USE constants + IMPLICIT NONE + REAL(8), INTENT(IN) :: rs, xk, temp + REAL(8), INTENT(OUT) :: erl + + INTEGER, PARAMETER :: nlt = 5, nlrs = 4, nrt = 2, nrrs = 3 + REAL(8), DIMENSION(:) :: tleft(5), rsleft(4), tright(3), rsright(3) + REAL(8), DIMENSION(27) :: left1, left2, left3, left4 + REAL(8), DIMENSION(15) :: right1, right2, right3, right4 + REAL(8) :: vleft, vright, weight, k_boundary, eee, left0, t, k, wp, rkf, ef + + REAL(8), DIMENSION(27) :: p + REAL(8), DIMENSION(15) :: qq + ! DATA p / 2.5344908519150277648180d0, -0.3404592250356968663461d0, 0.2090924044894818256690d0, & + ! -6.2531413812700469279093d0, -1.2882761371134692218732d0, -0.1545959102722999634061d0, & + ! 6.3092788968560729045976d0, -0.3712659829797280708341d0, 0.4561735015743797472254d0, & + ! 1.0517044701577911158807d0, 0.9980885389251320516379d0, -0.1262961858879391718880d0, & + ! -1.5864423552937068073732d0, -0.0703517047897053715566d0, 0.0012095316268887058939d0, & + ! 0.6265240382967289800575d0, -0.0345491408824314047421d0, 0.0475813589973176492842d0, & + ! 1.1214446081724578352379d0, 0.0000158118170277364426d0, 0.0378391859065193569833d0, & + ! 0.5760739892638068182507d0, -0.0591688129460346545763d0, 0.0696403645308140417658d0, & + ! -0.0209456339452350516483d0, 0.0320039948622821771029d0, -0.0007924053800554207690d0/ + + ! DATA qq / 78637966.0787817537784576416016d0, 311217296.7055425047874450683594d0, -12127688.5634898133575916290283d0, & + ! -0.1313637765537457313680d0, -0.2846930273653715537385d0, 0.0139082160012960597689d0, & + ! -0.7582546438957535883674d0, 3.7042943840493092189092d0, -0.2725344901748169568201d0, & + ! 3.0528050686602790086965d0, 0.7089354835299107593372d0, -0.6289703483429308628061d0, & + ! 0.4129987887383732569901d0, 0.2055510771521086854641d0, -0.1132810050944537821893d0 / + + + DATA left1/ 1.2154134208095004d0, -4.508816240359177d0, 11.97698761986273d0, & + 1.029009842611702d0, -7.243415266109083d0, 23.33202585152756d0, & + 1.2369761501564556d0, -6.223807061499287d0, 18.14943199920946d0, & + 0.40358858199494435d0, -4.354830869637849d0, 15.501627868057078d0, & + -6.525885979048827d0, 3.252201256859724d0, -0.6080473609125667d0, & + -8.622666469061805d0, 4.821019804126895d0, -0.9930420373493627d0, & + 12.073120786108186d0, -6.4928164712145415d0, 1.2909333493367048d0, & + 3.917630502317266d0, -2.2250639617884422d0, 0.4664176928931159d0, & + -0.0002913680954701001d0, 0.001102119504261212d0, -0.0003347159507350774d0 / + + DATA right1/ 62798400.694744214d0, 79885552.49571729d0, -61337388.52634566d0, & + 0.5346554263875789d0, -0.5651101777324324d0, 1.8738775250643693d0, & + 0.3610699572007695d0, -1.5637471619292d0, 4.107115105254552d0, & + -0.5644568492314478d0, -0.44104734176745714d0, 0.07004901166054882d0, & + -0.07218918308832134d0, 0.8775972530665443d0, 0.03470894995029971d0 / + + DATA left2/0.42478487081400573d0, -2.765565312702711d0, 1.413712045371367d0, & + 1.8008959412556513d0, -5.511781319389452d0, 2.875404181280661d0, & + -1.1165589022045965d0, 4.155085805524251d0, -2.150770107359026d0, & + 1.2917078801078052d0, -3.614898998919843d0, 1.9022101205756106d0, & + 5.350869687205111d0, -0.5026845175588389d0, -0.15396574514270045d0, & + 9.373849864823613d0, -1.6438575214764517d0, 0.004280818867507963d0, & + 11.694218235303818d0, -1.6291144624856047d0, -0.16101555512766996d0, & + -4.868323364472416d0, 0.9470904674414855d0, -0.04702950176365168d0, & + -0.012640460224829784d0, 0.019161157710872564d0, -0.004262966329670912d0 / + + DATA right2/ 8506976.524928845d0, 49450609.38150841d0, -19148117.040330067d0, & + 2.125260762220303d0, 0.39526176892125203d0, -0.2637442881384847d0, & + 2.0976711275164974d0, -0.5177622167149698d0, 0.9646134330597216d0, & + -0.3334561771389293d0, -0.20556723599763893d0, 0.07739791566506975d0, & + 0.22005329730563267d0, 0.3222403436910038d0, -0.09636278516858485d0 / + + DATA left3/ 2.5344908519150278d0, -0.34045922503569687d0, 0.20909240448948183d0, & + -6.253141381270047d0, -1.2882761371134692d0, -0.15459591027229996d0, & + 6.309278896856073d0, -0.37126598297972807d0, 0.45617350157437975d0, & + 1.0517044701577911d0, 0.998088538925132d0, -0.12629618588793917d0, & + -1.5864423552937068d0, -0.07035170478970537d0, 0.001209531626888706d0, & + 0.626524038296729d0, -0.034549140882431405d0, 0.04758135899731765d0, & + 1.1214446081724578d0, 1.5811817027736443d-5, 0.03783918590651936d0, & + 0.5760739892638068d0, -0.059168812946034655d0, 0.06964036453081404d0, & + -0.02094563394523505d0, 0.03200399486228218d0, -0.0007924053800554208d0 / + + DATA right3/ 78637966.07878175d0, 311217296.7055425d0, -12127688.563489813d0, & + -0.13136377655374573d0, -0.28469302736537155d0, 0.01390821600129606d0, & + -0.7582546438957536d0, 3.704294384049309d0, -0.27253449017481696d0, & + 3.052805068660279d0, 0.7089354835299108d0, -0.6289703483429309d0, & + 0.41299878873837326d0, 0.20555107715210869d0, -0.11328100509445378d0 / + + DATA left4/ 0.2993158634964539d0, 1.625125795930575d0, -0.2287001580618493d0, & + 1.9351156358534447d0, 1.486385504665352d0, -0.27231997267857966d0, & + 1.2396680043134514d0, 1.9448397131933621d0, -0.31003687759808185d0, & + 2.328977751230831d0, 0.8835261791434674d0, -0.19774034154529302d0, & + -1.3242329783616744d0, 0.0950061930026322d0, 0.01668566331718798d0, & + -1.3819322069998992d0, 0.29087403267083034d0, -0.007844189680559794d0, & + 1.828253099747466d0, -0.2607012137868943d0, -0.007553184872345978d0, & + 0.4075947748403373d0, -0.10769700661487398d0, 0.0055513391816981806d0, & + -0.5402740829426358d0, 0.058999510608476025d0, 0.004566792900479818d0 / + + DATA right4/ 25343208.713695824d0, -28305801.47662847d0, -74437595.52635065d0, & + 0.7558224839914135d0, -1.0500056608761195d0, 0.1563053901402796d0, & + -1.7177311660206924d0, 1.7866351771464843d0, -0.2807943980619071d0, & + 2.6934109989894948d0, -1.2072374495844942d0, 0.5334795894434334d0, & + 3.1140946102853815d0, -1.6792786833930664d0, 0.9473416778998698d0 / + + !------------------------------------------------------------ + ! Preprocessing + !------------------------------------------------------------ + rkf = fa/rs + ef = rkf**2.d0/2.d0 + wp = SQRT(3.d0/rs**3.d0) + eee = -pi*wp/(4.d0*rkf*ef) + k = xk / rkf + t = temp / ef + + !------------------------------------------------------------ + ! Check for interpolation or extrapolation + !------------------------------------------------------------ + IF (t.GT.2.d0) THEN + PRINT*, "ERROR: Local electronic temperature is higher interpolation grid (T > 2 TF)." + PRINT*, "T = ", t, temp + PRINT*, "rs = ", rs + CALL par_stop + ENDIF + + IF (rs.GT.20.d0) THEN + PRINT*, "ERROR: rs > 20 is outside of fit range" + CALL par_stop + ENDIF + + IF (rs.LT.1.d-2) THEN + PRINT*, "ERROR: rs < 0.01 is outside of fit range" + CALL par_stop + ENDIF + + IF (k.GT.20.d0) THEN + PRINT*, "ERROR: k is higher interpolation grid (k > 20kF)." + PRINT*, "T = ", t, temp + PRINT*, "rs = ", rs + PRINT*, "k = ", k + CALL par_stop + ENDIF + + ! Set to fermi level if below fermi level + IF (k.LT.1.00001d0) THEN + k = 1.00001d0 + ENDIF + + ! Determine which region of rs to use + IF (rs.LT.0.2) THEN + p = left1 + qq = right1 + ELSEIF (rs.LT.1.0) THEN + p = left2 + qq = right2 + ELSEIF (rs.LT.5.0) THEN + p = left3 + qq = right3 + ELSE + p = left4 + qq = right4 + ENDIF + + !-------------------- + ! Real part + !-------------------- + CALL real_cusp(rs, t, k_boundary) + + ! Reference + + + ! Left hand side + rsleft(1) = p(1) * rs + p(2)*rs**(3.d0/2.d0) + p(3)*rs**2.d0 + rsleft(2) = p(4) * rs + p(5)*rs**(3.d0/2.d0) + p(6)*rs**2.d0 + rsleft(3) = p(7) * rs + p(8)*rs**(3.d0/2.d0) + p(9)*rs**2.d0 + rsleft(4) = p(10) * rs + p(11)*rs**(3.d0/2.d0) + p(12)*rs**2.d0 + tleft(1) = p(13) + p(14)*t + p(15)*t**2.d0 + tleft(2) = p(16) + p(17)*t + p(18)*t**2.d0 + tleft(3) = p(19) + p(20)*t + p(21)*t**2.d0 + tleft(4) = p(22) + p(23)*t + p(24)*t**2.d0 + tleft(5) = p(25) + p(26)*t + p(27)*t**2.d0 + vleft = tleft(1)*rsleft(1)*k & + + tleft(2)*rsleft(2)*k**(2.d0) & + + tleft(3)*rsleft(3)*k**(1.5d0) & + + tleft(4)*rsleft(4)*k**(2.5d0) & + + tleft(5) + + left0 = tleft(1)*rsleft(1)*k_boundary & + + tleft(2)*rsleft(2)*k_boundary**(2.d0) & + + tleft(3)*rsleft(3)*k_boundary**(1.5d0) & + + tleft(4)*rsleft(4)*k_boundary**(2.5d0) & + + tleft(5) + + ! Right hand side + rsright(1) = qq(1) * rs + qq(2)*rs**(3.d0/2.d0) + qq(3)*rs**2.d0 + rsright(2) = qq(4) * rs + qq(5)*rs**(3.d0/2.d0) + qq(6)*rs**2.d0 + rsright(3) = qq(7) * rs + qq(8)*rs**(3.d0/2.d0) + qq(9)*rs**2.d0 + tright(2) = qq(10) + qq(11)*t + qq(12)*t**2.d0 + tright(3) = qq(13) + qq(14)*t + qq(15)*t**2.d0 + tright(1) = k_boundary*(left0 - tright(2)*rsright(2)/k_boundary**2 - tright(3)*rsright(3)/k_boundary**(3))/(rsright(1)*eee) + vright = eee*tright(1)*rsright(1)/k + tright(2)*rsright(2)/k**2 + tright(3)*rsright(3)/k**(3) + + + ! Interpolation + weight = 1.d0/(1.d0 + EXP(15.d0*(k - k_boundary))) + erl = weight*vleft + (1.d0-weight)*vright + + ! Scale back + erl = erl * ef + RETURN +END + +SUBROUTINE real_cusp(rs, t, k) + ! Position of minimum for real sigma + ! + ! Inputs: + ! rs - radius + ! t - temperature in unit of fermi + ! Output: + ! k - temperature in unit of fermi + ! + IMPLICIT NONE + REAL(8), INTENT(IN) :: rs, t + REAL(8), INTENT(OUT) :: k + REAL(8), PARAMETER :: alpha = 6.4226400052591925d0 + REAL(8), PARAMETER :: beta = 0.06579409601454807d0 + REAL(8), DIMENSION(9) :: factors + REAL(8) :: A, B, C + DATA factors / 11.789533109217384d0, 25.56619569021551d0, -8.718509700959657d0, & + 0.07937961365220357d0, 0.1642978703589978d0, -0.0019130964169998144d0, & + 1.4999860733980672d0, -3.021030697663565d-6, -0.0005684185963082104d0 / + + A = factors(1) + factors(2)*t + factors(3)*t**2.d0 + B = factors(4) + factors(5)*t + factors(6)*t**2.d0 + C = factors(7) + factors(8)*t + factors(9)*t**2.d0 + + IF (rs.LE.1.0d0) THEN + k = 1.5d0 + ELSE + k = (1.d0 + TANH(alpha*rs-A)) * (beta*LOG(rs)**2.d0 + B) + C + ENDIF + + RETURN +END SUBROUTINE + \ No newline at end of file diff --git a/src/EXCH/xcpot.f90 b/src/EXCH/xcpot.f90 index 96a288d..b0d5e70 100644 --- a/src/EXCH/xcpot.f90 +++ b/src/EXCH/xcpot.f90 @@ -5,479 +5,683 @@ ! $Date: 2012/03/29 22:52:37 $ ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine xcpot (iph, ie, index, lreal, ifirst, jri, em, xmu, & - & vtot, vvalgs, densty, dmag, denval, & - & eref, v, vval, iPl, WpCorr, Gamma, AmpFac, EGap, & - & vxcrmu, vxcimu, gsrel, vvxcrm, vvxcim) - - USE IOMOD - use dimsmod, only: nrptx, MxPole - use constants - implicit double precision (a-h, o-z) -! calculate self-energy correction -! first coded j. mustre de leon -! last modified a.ankudinov 1996 for non-local self-energies -! Ankudinov, Rehr, J. Physique IV, vol. 7, C2-121, (1997). - -! INPUT -! iph, ie used only for debug and labels. -! index 0 Hedin-Lunqvist + const real & imag part -! 1 Dirac-Hara + const real & imag part -! 2 ground state + const real & imag part -! 3 Dirac-Hara + HL imag part + const real & imag part -! 4 See rdinp for comment -! lreal not equal zero for real self energy -! ifirst first entry flag, set to zero before first call for -! each unique potential, see vxcrmu and vxcimu below -! jri index of first interstitial point in current -! Loucks r grid -! em current energy grid point -! xmu fermi level -! vi0 const imag part to subtract from potential -! gamach core hole lifetime -! vtot(nr) total potential (coulomb and gs exchange corr) -! vvalgs(nr) total coulomb + gs xc potential from valence electrons -! densty(nr) electron density -! dmag(nr) density magnetization -! denval(nr) valence electron density -! iPl Control for many pole self energy (Josh) -! -! OUTPUT -! eref complex energy reference for current energy -! v(nr) complex potential including energy dep xc -! vval(nr) as above,but xc from valence electrons only -! em current energy -! -! WORKSPACE -! vxcrmu and vxcimu are calculated only on first entry for a -! particular unique potential, re-used on subsequent entries. -! vxcrmu(nr) real part of xc at fermi level -! vxcimu(nr) imag part of xc at fermi level -! gsrel(nr) ratio of gs xc potentials with and without magnetization -! vvxcrm(nr) real part of xc at fermi level from valence electrons -! vvxcim(nr) imag part of xc at fermi level from valence electrons - - - - dimension vtot(nrptx), vvalgs(nrptx), densty(nrptx) - dimension dmag(nrptx), denval(nrptx) - complex*16 em, eref, v(nrptx), vval(nrptx) - dimension vxcrmu(nrptx), vxcimu(nrptx) - dimension vvxcrm(nrptx), vvxcim(nrptx), gsrel(nrptx) - -! Josh added variables: -! ZRnrm - renormalization constant -! csig - control for using many pole self energy -! NRPts - Number of points to use inside atom. -! Other points are linearly interpolated. -! WpCorr - Array of frequencies for many pole self energy. -! rsTmp - Temp var for rs -! WpTmp - Temp var for Wp -! AmpFac - g_i (pole strengths) -! Gamma - pole broadening -! RsInt - Rs in the intersitial -! DRs - RsCore - RsInt -! delrHL - Re[delta sigma] for many pole self energy -! deliHL - Im[delta sigma] for many pole self energy -! Rs1(NRPts) - Array of Rs points for interpolation - - character*512 slog - logical csig - integer NRPts, NPoles - parameter (tol=0.0004) - parameter (NRPts=10) - complex*16 delta, deltav, ZRnrm, cmu, SigF(NRPts), deltaHL(NRPts) - double precision WpCorr(MxPole), rsTmp, WpTmp, & - & AmpFac(MxPole), Gamma(MxPole) - double precision RsInt, DRs, delrHL(NRPts), deliHL(NRPts), & - & Rs1(NRPts), EGap, RsMin, RsMax - character(LEN=10) ColumnLabels(20) - save SigF -! Josh END - -! First calculate vxc to correct the local momentum dispersion -! relation, delta = vxc(e,k) - vxc(mu,k), and -! p^2 = k^2 -mu + kf^2 - delta. -! In jr theory, v(e,r) = vcoul(r) + vxc(e,r) = -! = vcoul(r) + vxcgs(r) + delta(e,r). - - -! at jri potential is smooth continuation of potential to r(jri) -! at this point potential jumps to interstitial value at jri+1 - csig=.false. - jri1 = jri + 1 - nmax=1 - nul=0 - ibp = index / 10 - ixc = mod(index,10) - ixcTmp=ixc - ZRNrm = 1.d0 - if((iPl.gt.0).and.(ixc.eq.0)) then - csig=.true. - end if - if (ixc .eq. 2 .or. dble(em).le.xmu) then - do 10 i = 1, jri1 - v(i) = vtot(i) - vval(i) = vvalgs(i) - 10 continue -! Ground state exchange, no self energy calculation - goto 888 - endif -! Josh - Added CSigma to calculate HL Sigma with broadening and such. -! Calculate Rs at the core and interstitial densities. - if (densty(jri1).le.0) then - RsInt =10 - else - RsInt = (3 / (4*pi*densty(jri1))) ** third - endif - if (densty(1).le.0) then - rscore =101.d0 - else - rscore = (3 / (4*pi*densty(1))) ** third - endif - if (MAXVAL(densty(1:jri1)).le.0) then - RsMin =RsInt*1.d-3 - else - RsMin = (3 / (4*pi*MAXVAL(densty(1:jri1)))) ** third - endif - if (MINVAL(densty(1:jri1)).le.0) then - RsMax=RsInt*2.d0 - else - RsMax = (3 / (4*pi*MINVAL(densty(1:jri1)))) ** third - endif - !RsMax = RsInt - - DRs = (RsMax-RsMin)/(NRPts-2) - IF(iPl.EQ.2) THEN - Rs1(1) = RsMin - DO i= 2, NRPts - IF(RsMax.GT.RsInt) THEN - IF(Rs1(i-1).LT.RsInt) THEN - Rs1(i) = RsMin + DBLE(i-1)*DRs - IF(Rs1(i).GT.RsInt) Rs1(i) = RsInt - ELSE - IF(i.GT.2) THEN - Rs1(i)=RsMin+DBLE(i-2)*DRs - END IF - END IF - ELSE - Rs1(i) = RsMin + DBLE(i-1)*DRs - END IF - END DO +SUBROUTINE xcpot(iph, ie, index, lreal, ifirst, jri, em_in, xmu, & + & vtot, vvalgs, densty, dmag, denval, & + & eref, v, vval, iPl, WpCorr, Gamma, AmpFac, EGap, & + & vxcrmu, vxcimu, gsrel, vvxcrm, vvxcim, fname) + + USE IOMOD + USE dimsmod, only: nrptx, MxPole + USE constants + USE potential_inp, ONLY: scf_temperature + IMPLICIT NONE + + ! implicit double precision (a-h, o-z) + ! calculate self-energy correction + ! first coded j. mustre de leon + ! last modified a.ankudinov 1996 for non-local self-energies + ! Ankudinov, Rehr, J. Physique IV, vol. 7, C2-121, (1997). + + ! INPUT + ! iph, ie used only for debug and labels. + ! index 0 Hedin-Lunqvist + const real & imag part + ! 1 Dirac-Hara + const real & imag part + ! 2 ground state + const real & imag part + ! 3 Dirac-Hara + HL imag part + const real & imag part + ! 4 See rdinp for comment + ! lreal not equal zero for real self energy + ! ifirst first entry flag, set to zero before first call for + ! each unique potential, see vxcrmu and vxcimu below + ! jri index of first interstitial point in current + ! Loucks r grid + ! em current energy grid point + ! xmu fermi level + ! vi0 const imag part to subtract from potential + ! gamach core hole lifetime + ! vtot(nr) total potential (coulomb and gs exchange corr) + ! vvalgs(nr) total coulomb + gs xc potential from valence electrons + ! densty(nr) electron density + ! dmag(nr) density magnetization + ! denval(nr) valence electron density + ! iPl Control for many pole self energy (Josh) + INTEGER, INTENT(INOUT) :: iph, ie, index, lreal, jri, iPl + INTEGER, INTENT(INOUT) :: ifirst + COMPLEX*16, INTENT(INOUT) :: em_in + REAL*8, INTENT(INOUT) :: xmu + REAL*8, DIMENSION(nrptx), INTENT(INOUT) :: vtot, vvalgs, dmag, denval, densty + ! OUTPUT + ! eref complex energy reference for current energy + ! v(nr) complex potential including energy dep xc + ! vval(nr) as above,but xc from valence electrons only + ! em current energy + COMPLEX*16, INTENT(INOUT) :: eref + COMPLEX*16, DIMENSION(nrptx), INTENT(INOUT) :: v, vval + + ! WORKSPACE + ! vxcrmu and vxcimu are calculated only on first entry for a + ! particular unique potential, re-used on subsequent entries. + ! vxcrmu(nr) real part of xc at fermi level + ! vxcimu(nr) imag part of xc at fermi level + ! gsrel(nr) ratio of gs xc potentials with and without magnetization + ! vvxcrm(nr) real part of xc at fermi level from valence electrons + ! vvxcim(nr) imag part of xc at fermi level from valence electrons + REAL*8, DIMENSION(nrptx), INTENT(INOUT) :: gsrel + REAL*8, DIMENSION(nrptx), INTENT(INOUT) :: vxcrmu, vxcimu + REAL*8, DIMENSION(nrptx), INTENT(INOUT) :: vvxcrm, vvxcim + CHARACTER(LEN=50) :: fname + + REAL*8 :: xkpp, xkm2, xkm, xk2, xk, xfval, xf, xfm + REAL*8 :: vxcvi, vxcvr, vxci, vxcr, rsm, rsval, rs, rscore + REAL*8 :: delr, delvi, delvr, deli, del1r, tol + COMPLEX*16 :: em + INTEGER :: nul, nmax, jri1, niter, ixc, ixcTmp, i, ibp + ! Josh added variables: + ! ZRnrm - renormalization constant + ! csig - control for using many pole self energy + ! NRPts - Number of points to use inside atom. + ! Other points are linearly interpolated. + ! WpCorr - Array of frequencies for many pole self energy. + ! rsTmp - Temp var for rs + ! WpTmp - Temp var for Wp + ! AmpFac - g_i (pole strengths) + ! Gamma - pole broadening + ! RsInt - Rs in the intersitial + ! DRs - RsCore - RsInt + ! delrHL - Re[delta sigma] for many pole self energy + ! deliHL - Im[delta sigma] for many pole self energy + ! Rs1(NRPts) - Array of Rs points for interpolation + + CHARACTER*512 slog + LOGICAL csig + INTEGER NRPts, NPoles + PARAMETER (tol=0.0004) + PARAMETER (NRPts=10) + COMPLEX*16 :: delta, deltav, ZRnrm, cmu, SigF(NRPts), deltaHL(NRPts) + DOUBLE PRECISION WpCorr(MxPole), rsTmp, WpTmp, AmpFac(MxPole), Gamma(MxPole) + DOUBLE PRECISION RsInt, DRs, delrHL(NRPts), deliHL(NRPts), Rs1(NRPts), EGap, RsMin, RsMax + CHARACTER(LEN=10) ColumnLabels(20) + SAVE SigF + ! Josh END + + ! First calculate vxc to correct the local momentum dispersion + ! relation, delta = vxc(e,k) - vxc(mu,k), and + ! p^2 = k^2 -mu + kf^2 - delta. + ! In jr theory, v(e,r) = vcoul(r) + vxc(e,r) = + ! = vcoul(r) + vxcgs(r) + delta(e,r). + + + ! at jri potential is smooth continuation of potential to r(jri) + ! at this point potential jumps to interstitial value at jri+1 + csig = .FALSE. + jri1 = jri + 1 + nmax = 1 + nul = 0 + ibp = index / 10 + ixc = MOD(index, 10) + ixcTmp = ixc + ZRNrm = 1.d0 + IF ((iPl.GT.0).AND.(ixc.EQ.0)) THEN + csig = .TRUE. + ENDIF + + ! IF (ixc .EQ. 2 .OR. DBLE(em).LE.xmu) THEN + ! DO i = 1, jri1 + ! v(i) = vtot(i) + ! vval(i) = vvalgs(i) + ! ENDDO + ! ! Ground state exchange, no self energy calculation + ! GOTO 888 + ! ENDIF + + ! TTS - Apply self-energy calculation below Fermi level at finite temperature + em = em_in + IF (index.EQ.2) THEN + DO i = 1, jri1 + v(i) = vtot(i) + vval(i) = vvalgs(i) + ENDDO + ! Ground state exchange, no self energy calculation + GOTO 888 + ELSE!IF ((index.EQ.6).OR.(index.EQ.7)) THEN + ! For finite temperature self-energies + IF (DBLE(em_in).LE.xmu) THEN + ! Use self energy at xmu for E < xmu + IF (xmu.LT.0.0d0) THEN + em = DCMPLX(xmu*0.9999d0, DIMAG(em_in)) + ELSE + em = DCMPLX(xmu*1.0001d0, DIMAG(em_in)) + ENDIF + ENDIF + ENDIF + IF (index.EQ.-2) index = 6 ! JJK - calculate self-energy if temperature dependent ground state + + ! Josh - Added CSigma to calculate HL Sigma with broadening and such. + ! Calculate Rs at the core and interstitial densities. + IF (densty(jri1).LE.0) THEN + RsInt =10 + ELSE + RsInt = (3 / (4*pi*densty(jri1))) ** third + ENDIF + IF (densty(1).LE.0) THEN + rscore =101.d0 + ELSE + rscore = (3 / (4*pi*densty(1))) ** third + ENDIF + IF (MAXVAL(densty(1:jri1)).LE.0) THEN + RsMin =RsInt*1.d-3 + ELSE + RsMin = (3 / (4*pi*MAXVAL(densty(1:jri1)))) ** third + ENDIF + IF (MINVAL(densty(1:jri1)).LE.0) THEN + RsMax=RsInt*2.d0 + ELSE + RsMax = (3 / (4*pi*MINVAL(densty(1:jri1)))) ** third + ENDIF + + !RsMax = RsInt + + DRs = (RsMax-RsMin)/(NRPts-2) + IF (iPl.EQ.2) THEN + Rs1(1) = RsMin + DO i= 2, NRPts + IF (RsMax.GT.RsInt) THEN + IF (Rs1(i-1).LT.RsInt) THEN + Rs1(i) = RsMin + DBLE(i-1)*DRs + IF (Rs1(i).GT.RsInt) Rs1(i) = RsInt + ELSE + IF (i.GT.2) THEN + Rs1(i) = RsMin + DBLE(i-2)*DRs + END IF + END IF + ELSE + Rs1(i) = RsMin + DBLE(i-1)*DRs + END IF + END DO + ELSE + Rs1(:) = RsInt + END IF + !WRITE(66,*) Rs1 + ! Now calculate delta sigma as a function of Rs and energy + IF (csig) THEN + ! Count the number of poles + DO i = 1, MxPole + IF (WpCorr(i).LT.-1.d0) THEN + NPoles = i - 1 + EXIT + ENDIF + ENDDO + ! Self energy at the fermi level is calculated once only since it is independent of + ! energy. + IF (ifirst.EQ.0) THEN + cmu = xmu*1.00001d0 + DO i = 1, NRPts + SigF(i) = 0.d0 + !Rs1(i)=RsMin+DBLE(i-1)*DRs + + ! If iPl = 1, use r independent sigma with parameters calculated + ! from the interstitial density, i.e. bulk self-energy + ! If iPl = 2, use Sigma(r) = Sigma[Wp(r)*Wp/Wp(RInt)] + IF ((iPl.EQ.2).OR.(i.EQ.NRPts)) THEN + ! The next line will not work with ipl = 2. Get rid of ipl = 2 in future. + WpCorr(:) = WpCorr(:)*SQRT(3.d0/Rs1(i)**3) + CALL CSigZ(cmu,xmu,Rs1(i),SigF(i),ZRnrm,WpCorr,Gamma,AmpFac,EGap,NPoles,.TRUE.,.FALSE.) + WpCorr(:) = WpCorr(:)/SQRT(3.d0/Rs1(i)**3) + ENDIF + ENDDO + IF (iPl.EQ.1) THEN + SigF(1:NRPts-1) = SigF(NRPts) + ENDIF + ENDIF + + DO i = 1, NRPts + deltaHL(i) = 0.d0 + !Rs1(i)=RsMin+DBLE(i-1)*DRs + ! If iPl = 1, use r independent sigma with parameters calculated + ! from the interstitial density, i.e. bulk self-energy + ! If iPl = 2, use Sigma(r) = Sigma[Wp(r)*Wp/Wp(RInt)] + IF ((iPl.EQ.2).OR.(i.EQ.NRPts)) THEN + WpCorr(:) = WpCorr(:)*SQRT(3.d0/Rs1(i)**3) + CALL CSigZ(em,xmu,Rs1(i),deltaHL(i),ZRnrm,WpCorr,Gamma,AmpFac,EGap,NPoles,.TRUE.,.FALSE.) + WpCorr(:) = WpCorr(:)/SQRT(3.d0/Rs1(i)**3) + deltaHL(i) = ZRnrm*(deltaHL(i) - SigF(i)) + ENDIF + ENDDO + IF (iPl.EQ.1) THEN + deltaHL(1:NRPts-1) = deltaHL(NRPts) + ENDIF + delrHL(:) = DBLE(deltaHL(:)) + deliHL(:) = DIMAG(deltaHL(:)) + ENDIF + + ! END Josh + ! Add the self energy correction + DO i = jri1, 1, -1 + niter = 0 + IF (densty(i).LE.0) THEN + rs = 10 + ELSE + rs = (3 / (4*pi*densty(i))) ** third + ENDIF + IF(ixc.EQ.-2) rs = RsInt ! DEBUG - JK use only constant imaginary part at the interstitial for + ! temperature dependent ground state calculations. + ! write(22,*) 1.d0*exp(dble(i)*0.01), densty(i) + ! Josh - If csigma is turned on, interpolate onto rs. + ! Then skip to 15 (skip other calculations and self + ! consistency) + ! Test with constant SE. + IF (.FALSE.) THEN + delrHL(:) = 0.d0/hart + deliHL(:) = -5.d0/hart + ENDIF + + IF (csig) THEN + IF (iPl.EQ.2) THEN + IF ((rs.LT.RsMin).OR.(rs.GT.RsMax)) THEN + delrHL = 0.d0 + ELSE + CALL terp(Rs1, delrHL, NRPts, 1, rs, delr) + CALL terp(Rs1, deliHL, NRPts, 1, rs, deli) + END IF ELSE - Rs1(:) = RsInt + delr = delrHL(1) + deli = deliHL(1) END IF - !WRITE(66,*) Rs1 -! Now calculate delta sigma as a function of Rs and energy - if (csig) then - ! Count the number of poles - do i = 1, MxPole - if(WpCorr(i).lt.-1.d0) then - NPoles = i - 1 - EXIT - end if - end do - ! Self energy at the fermi level is calculated once only since it is independent of - ! energy. - if(ifirst .eq. 0) then - cmu = xmu*1.00001d0 - do i= 1, NRPts - SigF(i) = 0.d0 - !Rs1(i)=RsMin+DBLE(i-1)*DRs - - ! If iPl = 1, use r independent sigma with parameters calculated - ! from the interstitial density, i.e. bulk self-energy - ! If iPl = 2, use Sigma(r) = Sigma[Wp(r)*Wp/Wp(RInt)] - if((iPl.eq.2).or.(i.eq.NRPts)) then - ! The next line will not work with ipl = 2. Get rid of ipl = 2 in future. - WpCorr(:) = WpCorr(:)*SQRT(3.d0/Rs1(i)**3) - call CSigZ(cmu,xmu,Rs1(i),SigF(i),ZRnrm,WpCorr,Gamma,AmpFac,EGap,NPoles,.TRUE.,.FALSE.) - WpCorr(:) = WpCorr(:)/SQRT(3.d0/Rs1(i)**3) - end if - end do - if (iPl.eq.1) then - SigF(1:NRPts-1) = SigF(NRPts) - end if - end if - - do i= 1, NRPts - deltaHL(i) = 0.d0 - !Rs1(i)=RsMin+DBLE(i-1)*DRs - ! If iPl = 1, use r independent sigma with parameters calculated - ! from the interstitial density, i.e. bulk self-energy - ! If iPl = 2, use Sigma(r) = Sigma[Wp(r)*Wp/Wp(RInt)] - if((iPl.eq.2).or.(i.eq.NRPts)) then - WpCorr(:) = WpCorr(:)*SQRT(3.d0/Rs1(i)**3) - call CSigZ(em,xmu,Rs1(i),deltaHL(i),ZRnrm,WpCorr,Gamma,AmpFac,EGap,NPoles,.TRUE.,.FALSE.) - WpCorr(:) = WpCorr(:)/SQRT(3.d0/Rs1(i)**3) - deltaHL(i) = ZRnrm*(deltaHL(i) - SigF(i)) - end if - end do - if (iPl.eq.1) then - deltaHL(1:NRPts-1) = deltaHL(NRPts) - end if - delrHL(:) = DBLE(deltaHL(:)) - deliHL(:) = DIMAG(deltaHL(:)) - - end if - -! END Josh - -! Add the self energy correction - do 20 i = jri1,1,-1 - niter = 0 - if (densty(i).le.0) then - rs =10 - else - rs = (3 / (4*pi*densty(i))) ** third - endif -! write(22,*) 1.d0*exp(dble(i)*0.01), densty(i) -! Josh - If csigma is turned on, interpolate onto rs. -! Then skip to 15 (skip other calculations and self -! consistency) - ! Test with constant SE. - if(.FALSE.) THEN - delrHL(:) = 0.d0/hart - deliHL(:) = -5.d0/hart - end if - - if(csig) then - IF(iPl.EQ.2) THEN - IF((rs.LT.RsMin).OR.(rs.GT.RsMax)) THEN - delrHL = 0.d0 - ELSE - call terp (Rs1, delrHL, NRPts, 1, rs, delr) - call terp (Rs1, deliHL, NRPts, 1, rs, deli) - END IF - ELSE - delr = delrHL(1) - deli = deliHL(1) - END IF - goto 15 - end if -! END Josh - -! xf = 1.9191.../rs - xf = fa / rs - rsm = rs / (1+dmag(i))**third - xfm = fa / rsm - - if (ixc.eq.5) then - if ( denval(i) .gt. 0.00001) then - rsval = (3 / (4*pi*denval(i))) ** third - if (rsval.gt.10.0) rsval=10.0 - else - rsval = 10.0 - endif - xfval = fa / rsval - elseif (ixc.ge.6) then - if (densty(i) .le. denval(i) ) then - rscore = 101.0 - else - rscore = (3 / (4*pi*(densty(i)-denval(i)))) ** third - endif - endif - - if (ifirst .eq. 0) then -! vxc_mu indep of energy, calc only once -! Calculate vxc at fermi level e = mu, j.m. 1/12/89 - xk = xf * 1.00001 - gsrel(i) = 1.0d0 - if (ixc .lt. 5) then - call sigma(ixc, ibp,rs,rscore,xk,vxcrmu(i),vxcimu(i)) - if (index .eq. 0) then -! do not need 4 following lines for gs difference in potential -! xmag = 1.0d0+ dmag(i) -! call vbh(rs,xmag,v1) -! call vbh(rs, 1.0d0,v0) -! if (v0 .ne. 0) gsrel(i) = v1/v0 - endif - else - call sigma(nul,ibp, rs, rscore,xk,vxcrmu(i),vxcimu(i)) - endif - if (ixc.eq.5 ) then - xkpp = xfval * 1.00001 - call sigma & - & (ixc, ibp, rsval, rscore, xkpp, vvxcrm(i),vvxcim(i)) - if (ixc.eq.5 .and. i.eq.jri1) then - vvxcrm(jri1) = vxcrmu(jri1) - vvxcim(jri1) = vxcimu(jri1) - endif - elseif (ixc .ge. 6) then - call sigma & - & (ixc, ibp, rs, rscore, xk, vvxcrm(i), vvxcim(i)) - if (ixc.eq.6 .and. i.eq.jri1) then - vvxcrm(jri1) = vxcrmu(jri1) - vvxcim(jri1) = vxcimu(jri1) - endif - else - vvxcrm(i) = 0.0d0 - vvxcim(i) = 0.0d0 - endif - endif - -! xk2 is the local momentum squared, p^2 = k^2 - 2*mu + kf^2, -! k^2 represents energy measured from vacuum. -! See formula 2.15 in Lee and Beni's paper with the last 2 -! terms neglected. (complete reference?) - xk2 = 2 * (dble(em) - xmu) + xf**2 - xk = sqrt(xk2) - xkm2 = 2 * (dble(em) - xmu) + xfm**2 -! quick fix - if (xkm2.lt.0) xkm2=xk2 - xkm = sqrt(xkm2) - -! find \delta_1 - if (ixc .lt. 5) then - call sigma (ixc, ibp, rs, rscore, xk, vxcr, vxci) - else - call sigma (nul, ibp, rs, rscore, xk, vxcr, vxci) - endif - del1r = gsrel(i) * (vxcr - vxcrmu(i)) - -! Correct local momentum according to the formula -! p^2 = k^2 - 2*mu + kf^2 - 2*delta. Note that imag part -! of delta is ignored, since xk2 is a real quantity. - -! find xk(em) by iterative solution of dyson equation - 50 continue - xk2 = 2*(dble(em) - xmu - del1r) + xf**2 - if (xk2 .lt. 0) then - write(slog,'(1pe13.5, 3i8, a)') xk2, i, ie, iph, ' xk2, i, ie, iph' - call wlog(slog) - call wlog(' em, xf**2, xmu, delta') - write(slog,'(1p, 5e13.5)') dble(em), xf**2, xmu, del1r - call wlog(slog) - call par_stop('XCPOT-2') - endif - xk = sqrt (xk2) - -! calculate \delta_2 and \delta_v,2 with the corrected -! local momentum - call sigma (ixc, ibp, rs, rscore, xk, vxcr, vxci) -! delta corrected calculated with new local momentum - delr = gsrel(i) * (vxcr - vxcrmu(i)) - deli = vxci-vxcimu(i) - - if (ixc.ge.5 .and. i.eq.jri1 .and. xk.gt.xf) then - if (ixc.eq.5 .or. ixc.eq.6) then - delvr = delr - delvi = deli - endif - endif - - if (niter.lt.nmax) then - del1r=delr - niter=niter+1 - go to 50 - endif - - if (ixc .ge. 5 .and. i.lt.jri1 .and. xk.gt.xf) then - if (ixc.eq.5) then - xkpp=sqrt(xk**2-xf**2+xfval**2) - call sigma (ixc, ibp, rsval,rscore,xkpp,vxcvr,vxcvi) - else - call sigma (ixc, ibp, rs, rscore, xk, vxcvr, vxcvi) - endif - delvr = vxcvr-vvxcrm(i) - delvi = vxcvi-vvxcim(i) - endif - -! Josh - Skip SC loop if CSigma is called. CSigma calculates self consistently. - 15 continue - - delta = dcmplx(delr,deli) - -! Josh - write out delta sigma at interstitial level to sigma.dat. - if(i.eq.jri1) then -! write(45,'(X,20e14.6)') (DBLE(em) - xmu)*hart, delr*hart, & -! & deli*hart, DBLE(ZRnrm), DIMAG(ZRnrm), & -! & SQRT(DBLE(ZRnrm)**2+DIMAG(ZRnrm)**2), & -! & ATAN2(DIMAG(ZRnrm),DBLE(ZRnrm)), & -! & SQRT(DBLE(em-xmu)/2.d0)/ABS(deli)*bohr - ColumnLabels(:) = ' ' - ColumnLabels(1) = 'E-Mu' - ColumnLabels(2) = 'Re[Sigma]' - ColumnLabels(3) = 'Im[Sigma]' - ColumnLabels(4) = 'Re[Z]' - ColumnLabels(5) = 'Im[Z]' - ColumnLabels(6) = '|Z|' - ColumnLabels(7) = 'Phase[Z]' - ColumnLabels(8) = 'IMFP' - CALL WriteData('mpse.dat', & ! Specify file name. - & Double1 = (DBLE(em) - xmu)*hart, & ! Specify 1st col. - & DComplex2 = delta*hart, & ! Specify 2nd col. - & DComplex3 = ZRnrm, & ! Specify 3nd col. - & Double4 = SQRT(DBLE(ZRnrm)**2+DIMAG(ZRnrm)**2), & ! Specify 4rd col. - & Double5 = ATAN2(DIMAG(ZRnrm),DBLE(ZRnrm)), & ! 5th col. - & Double6 = SQRT(DBLE(em-xmu)/2.d0)/ABS(deli)*bohr, & ! 6th col - & ColumnLabels = ColumnLabels) - end if -! Josh END - - if (ixc .eq. 5) delta = dcmplx(delr,delvi) - v(i) = vtot(i) + delta - if (ixc .ge. 5) then - deltav = dcmplx(delvr,delvi) - vval(i) = vvalgs(i) + deltav - endif - 20 continue - 25 continue - - ifirst = 1 - -! Reference the potential with respect to mt potential, ie, -! first interstitial point. v(jri1) = 0 - -! Note that the reference does not contain the core hole lifetime -! since the total atomic potential should have it. However in the -! perturbation deltav = v - vmt it cancels out. -! ( deltav = vat - igamma - (vatmt-igamma) ). - - 888 eref = v(jri1) - do 910 i = 1, jri1 - 910 v(i) = v(i) - eref - if (ixc.ge.5) then - do 920 i = 1, jri1 - 920 vval(i) = vval(i) - eref - else - do 930 i = 1, jri1 - 930 vval(i) = v(i) - endif - -! Real self energy, zero imag part - if (lreal.gt.0) then - do 950 i = 1, jri1 - v(i) = dble(v(i)) - if (ixc.gt.4) vval(i) = dble(vval(i)) - 950 continue - eref = dble(eref) - endif - - - return - end - - subroutine sigma (ixc, ibp, rs, rscore, xk, vr, vi) - implicit double precision (a-h, o-z) - - if ((ixc.eq.0 .or. ixc.ge.5) .and. ibp .eq. 0) then - call rhl (rs, xk, vr, vi) - elseif ((ixc.eq.0.or. ixc.ge.5) .and. ibp .eq. 1) then - call rhlbp (rs, xk, vr, vi) - elseif (ixc .eq. 1) then - vi = 0 - call edp(rs,xk,vr) - elseif (ixc .eq. 3) then - call edp(rs,xk,vr) - call imhl (rs,xk,vi,icusp) - endif - - if (ixc .ge. 6) then - call edp(rscore,xk,vrp) - vr = vr - vrp - endif - - return - end - + GOTO 15 + ENDIF + ! END Josh + + ! xf = 1.9191.../rs + xf = fa / rs + rsm = rs / (1+dmag(i))**third + xfm = fa / rsm + + ! IF (ixc.EQ.5) THEN + ! IF ( denval(i).GT.0.00001) THEN + ! rsval = (3 / (4*pi*denval(i))) ** third + ! IF (rsval.GT.10.0) rsval=10.0 + ! ELSE + ! rsval = 10.0 + ! ENDIF + ! xfval = fa / rsval + ! ELSEIF (ixc.GE.6) THEN + ! IF (densty(i).LE.denval(i)) THEN + ! rscore = 101.0 + ! ELSE + ! rscore = (3 / (4*pi*(densty(i)-denval(i)))) ** third + ! ENDIF + ! ENDIF + + IF ((index.EQ.5).OR.(index.EQ.15)) THEN + IF ( denval(i).GT.0.00001) THEN + rsval = (3 / (4*pi*denval(i))) ** third + IF (rsval.GT.10.0) rsval=10.0 + ELSE + rsval = 10.0 + ENDIF + xfval = fa / rsval + ELSEIF (index.EQ.9) THEN + IF (densty(i).LE.denval(i)) THEN + rscore = 101.0 + ELSE + rscore = (3 / (4*pi*(densty(i)-denval(i)))) ** third + ENDIF + ENDIF + + ! IF (ifirst.EQ.0) THEN + ! ! vxc_mu indep of energy, calc only once + ! ! Calculate vxc at fermi level e = mu, j.m. 1/12/89 + ! xk = xf * 1.00001 + ! gsrel(i) = 1.0d0 + ! IF (ixc.LT.5) THEN + ! CALL sigma(ixc, ibp,rs,rscore,xk,vxcrmu(i),vxcimu(i),.TRUE.) + ! IF (index.EQ.0) THEN + ! ! do not need 4 following lines for gs difference in potential + ! ! xmag = 1.0d0+ dmag(i) + ! ! call vbh(rs,xmag,v1) + ! ! call vbh(rs, 1.0d0,v0) + ! ! if (v0 .ne. 0) gsrel(i) = v1/v0 + ! ENDIF + ! ELSE + ! CALL sigma(nul,ibp, rs, rscore,xk,vxcrmu(i),vxcimu(i),.TRUE.) + ! ENDIF + ! IF (ixc.EQ.5) THEN + ! xkpp = xfval * 1.00001 + ! CALL sigma(ixc, ibp, rsval, rscore, xkpp, vvxcrm(i),vvxcim(i),.FALSE.) + ! IF (ixc.EQ.5 .AND. i.EQ.jri1) THEN + ! vvxcrm(jri1) = vxcrmu(jri1) + ! vvxcim(jri1) = vxcimu(jri1) + ! ENDIF + ! ELSEIF (ixc.GE.6) THEN + ! CALL sigma(ixc, ibp, rs, rscore, xk, vvxcrm(i), vvxcim(i),.FALSE.) + ! IF (ixc.EQ.6.AND.i.EQ.jri1) THEN + ! vvxcrm(jri1) = vxcrmu(jri1) + ! vvxcim(jri1) = vxcimu(jri1) + ! ENDIF + ! ELSE + ! vvxcrm(i) = 0.0d0 + ! vvxcim(i) = 0.0d0 + ! ENDIF + ! ENDIF + + IF (ifirst.EQ.0) THEN + ! vxc_mu indep of energy, calc only once + ! Calculate vxc at fermi level e = mu, j.m. 1/12/89 + xk = xf * 1.00001 + gsrel(i) = 1.0d0 + IF ((index.EQ.0).OR.(index.EQ.1).OR.(index.EQ.2).OR.(index.EQ.3).OR.(index.EQ.6).OR.(index.EQ.7) & + .OR.(index.EQ.11).OR.(index.EQ.13)) THEN + CALL sigma(index, ibp,rs, rscore, xk, vxcrmu(i), vxcimu(i), xmu, .TRUE.) + ! Only get the real part + ELSE + ! Use index = 0 -> HL + ! index = 5, 15, 9 + CALL sigma(nul, ibp, rs, rscore, xk, vxcrmu(i), vxcimu(i), xmu, .FALSE.) + ENDIF + + IF ((index.EQ.5).OR.(index.EQ.15)) THEN + xkpp = xfval * 1.00001 + CALL sigma(index, ibp, rsval, rscore, xkpp, vvxcrm(i),vvxcim(i), xmu, .FALSE.) + IF (((index.EQ.5).OR.(index.EQ.15)).AND. i.EQ.jri1) THEN + vvxcrm(jri1) = vxcrmu(jri1) + vvxcim(jri1) = vxcimu(jri1) + ENDIF + ELSEIF (index.EQ.9) THEN + ! index = 9 + CALL sigma(index, ibp, rs, rscore, xk, vvxcrm(i), vvxcim(i), xmu, .FALSE.) + IF (index.EQ.9.AND.i.EQ.jri1) THEN + vvxcrm(jri1) = vxcrmu(jri1) + vvxcim(jri1) = vxcimu(jri1) + ENDIF + ELSE + ! index = 0, 1, 11, 2, 3, 13, 6, 7 + vvxcrm(i) = 0.0d0 + vvxcim(i) = 0.0d0 + ENDIF + ENDIF + + ! xk2 is the local momentum squared, p^2 = k^2 - 2*mu + kf^2, + ! k^2 represents energy measured from vacuum. + ! See formula 2.15 in Lee and Beni's paper with the last 2 + ! terms neglected. Phys Rev B 15 (6) 2862 (1977) + xk2 = 2 * (DBLE(em) - xmu) + xf**2 + xk = SQRT(xk2) + xkm2 = 2 * (DBLE(em) - xmu) + xfm**2 + ! quick fix + IF (xkm2.LT.0) xkm2 = xk2 + xkm = SQRT(xkm2) + + ! find \delta_1 + ! IF (ixc .lt. 5) THEN + ! CALL sigma(ixc, ibp, rs, rscore, xk, vxcr, vxci,.FALSE.) + ! ELSE + ! CALL sigma(nul, ibp, rs, rscore, xk, vxcr, vxci,.FALSE.) + ! ENDIF + IF ((index.EQ.0).OR.(index.EQ.1).OR.(index.EQ.2).OR.(index.EQ.3).OR.(index.EQ.6).OR.(index.EQ.7) & + .OR.(index.EQ.11).OR.(index.EQ.13)) THEN + CALL sigma(index, ibp, rs, rscore, xk, vxcr, vxci, xmu, .FALSE.) + ELSE + ! index = 5, 15, 9, (4 ?) + CALL sigma(nul, ibp, rs, rscore, xk, vxcr, vxci, xmu, .FALSE.) + ENDIF + del1r = gsrel(i) * (vxcr - vxcrmu(i)) + + ! Correct local momentum according to the formula + ! p^2 = k^2 - 2*mu + kf^2 - 2*delta. Note that imag part + ! of delta is ignored, since xk2 is a real quantity. + + ! find xk(em) by iterative solution of dyson equation + 50 CONTINUE + + xk2 = 2*(DBLE(em) - xmu - del1r) + xf**2 + IF (xk2.LT.0) THEN + WRITE(slog,'(1pe13.5, 3i8, a)') xk2, i, ie, iph, ' xk2, i, ie, iph' + CALL wlog(slog) + CALL wlog(' em, xf**2, xmu, delta') + WRITE(slog,'(1p, 5e13.5)') DBLE(em), xf**2, xmu, del1r + CALL wlog(slog) + CALL par_stop('XCPOT-2') + ENDIF + xk = SQRT(xk2) + + ! calculate \delta_2 and \delta_v,2 with the corrected local momentum + CALL sigma(index, ibp, rs, rscore, xk, vxcr, vxci, xmu, .FALSE.) + + ! delta corrected calculated with new local momentum + delr = gsrel(i) * (vxcr - vxcrmu(i)) + deli = vxci - vxcimu(i) + + IF (((index.EQ.5).OR.(index.EQ.15).OR.(index.EQ.9)).AND.(i.EQ.jri1).AND.(xk.GT.xf)) THEN + delvr = delr + delvi = deli + ENDIF + + IF (niter.LT.nmax) THEN + del1r=delr + niter=niter+1 + GOTO 50 + ENDIF + + ! IF (ixc .GE. 5 .AND. i.LT.jri1 .AND. xk.GT.xf) THEN + ! IF (ixc.EQ.5) THEN + ! xkpp = SQRT(xk**2-xf**2+xfval**2) + ! CALL sigma(ixc, ibp, rsval,rscore,xkpp,vxcvr,vxcvi,.FALSE.) + ! ELSE + ! CALL sigma(ixc, ibp, rs, rscore, xk, vxcvr, vxcvi,.FALSE.) + ! ENDIF + ! delvr = vxcvr-vvxcrm(i) + ! delvi = vxcvi-vvxcim(i) + ! ENDIF + IF ( ((index.EQ.5).OR.(index.EQ.15).OR.(index.EQ.9)).AND.(i.LT.jri1).AND.(xk.GT.xf) ) THEN + IF ((index.EQ.5).OR.(index.EQ.15)) THEN + xkpp = SQRT(xk**2-xf**2+xfval**2) + CALL sigma(index, ibp, rsval, rscore, xkpp, vxcvr, vxcvi, xmu, .FALSE.) + ELSE + ! index = 9 + CALL sigma(index, ibp, rs, rscore, xk, vxcvr, vxcvi, xmu, .FALSE.) + ENDIF + delvr = vxcvr-vvxcrm(i) + delvi = vxcvi-vvxcim(i) + ENDIF + + ! Josh - Skip SC loop if CSigma is called. CSigma calculates self consistently. + 15 CONTINUE + + delta = DCMPLX(delr, deli) + + ! Josh - write out delta sigma at interstitial level to sigma.dat. + ! TTS - Replaced em with em_in + IF (i.EQ.jri1) THEN + ! write(45,'(X,20e14.6)') (DBLE(em) - xmu)*hart, delr*hart, & + ! & deli*hart, DBLE(ZRnrm), DIMAG(ZRnrm), & + ! & SQRT(DBLE(ZRnrm)**2+DIMAG(ZRnrm)**2), & + ! & ATAN2(DIMAG(ZRnrm),DBLE(ZRnrm)), & + ! & SQRT(DBLE(em-xmu)/2.d0)/ABS(deli)*bohr + ColumnLabels(:) = ' ' + ColumnLabels(1) = 'E-Mu' + ColumnLabels(2) = 'Re[Sigma]' + ColumnLabels(3) = 'Im[Sigma]' + ColumnLabels(4) = 'Re[Z]' + ColumnLabels(5) = 'Im[Z]' + ColumnLabels(6) = '|Z|' + ColumnLabels(7) = 'Phase[Z]' + ColumnLabels(8) = 'IMFP' + CALL WriteData(TRIM(fname)//'.dat', & ! Specify file name. + & Double1 = (DBLE(em_in) - xmu)*hart, & ! Specify 1st col. + & DComplex2 = delta*hart, & ! Specify 2nd col. + & DComplex3 = ZRnrm, & ! Specify 3nd col. + & Double4 = SQRT(DBLE(ZRnrm)**2+DIMAG(ZRnrm)**2), & ! Specify 4rd col. + & Double5 = ATAN2(DIMAG(ZRnrm),DBLE(ZRnrm)), & ! 5th col. + & Double6 = SQRT(DBLE(em-xmu)/2.d0)/ABS(deli)*bohr, & ! 6th col + & ColumnLabels = ColumnLabels) + ENDIF + ! Josh END + + ! IF (ixc.EQ.5) delta = dcmplx(delr,delvi) + IF ((index.EQ.5).OR.(index.EQ.15)) delta = DCMPLX(delr, delvi) + + v(i) = vtot(i) + delta + + IF ((index.EQ.5).OR.(index.EQ.15).OR.(index.EQ.9)) THEN + deltav = DCMPLX(delvr, delvi) + vval(i) = vvalgs(i) + deltav + ENDIF + ENDDO + + ifirst = 1 + + ! Reference the potential with respect to mt potential, ie, + ! first interstitial point. v(jri1) = 0 + + ! Note that the reference does not contain the core hole lifetime + ! since the total atomic potential should have it. However in the + ! perturbation deltav = v - vmt it cancels out. + ! ( deltav = vat - igamma - (vatmt-igamma) ). + +888 IF (ifirst.eq.0) THEN + +ENDIF + +eref = v(jri1) + +DO i = 1, jri1 + v(i) = v(i) - eref +ENDDO + +IF ((index.EQ.5).OR.(index.EQ.15).OR.(index.EQ.9)) THEN + DO i = 1, jri1 + vval(i) = vval(i) - eref + ENDDO +ELSE + ! index = 0, 1, 11, 2, 3, 13, 6, 7 + DO i = 1, jri1 + vval(i) = v(i) + ENDDO +ENDIF + +! Real self energy, zero imag part +IF (lreal.GT.0) THEN + DO i = 1, jri1 + v(i) = DBLE(v(i)) + ! IF (ixc.GT.4) vval(i) = DBLE(vval(i)) + IF ((index.EQ.5).OR.(index.EQ.9).OR.(index.EQ.15)) vval(i) = DBLE(vval(i)) + ENDDO + eref = DBLE(eref) +ENDIF +! We need to restore index to original value +IF (ixc.EQ.-2) index = -2 +RETURN +END SUBROUTINE xcpot + +SUBROUTINE sigma(ixc, ibp, rs, rscore, xk, vr, vi, xmu, realonly) + ! USE m_COHSEX, ONLY: gw_exact + USE potential_inp, ONLY: scf_temperature + USE constants, ONLY: hart + IMPLICIT NONE + ! Computes the self-energy given radius rs and momentum xk + ! ixc - index + ! ibp - Deprecated + ! rs - Radius + ! rscore - Radius at core + ! vr - Real part in a.u. + ! vi - Imaginary part in a.u. + ! realonly - Return only the real part of self-energy + INTEGER, INTENT(IN) :: ixc, ibp + REAL(8), INTENT(IN) :: rs, rscore, xk, xmu + REAL(8), INTENT(OUT) :: vr, vi + LOGICAL, INTENT(IN) :: realonly + + INTEGER :: icusp + REAL(8) :: temperature, vrp + + SELECT CASE(ixc) + CASE(0) + ! Hedin-Lundqvist + constant imaginary + CALL rhl(rs, xk, vr, vi) + CASE(1) + ! Dirac-Hara + constant imaginary + vi = 0.d0 + CALL edp(rs, xk, vr) + CASE(2) + ! Ground state + const real & imag part + ! This seems to be redundant since we never call sigma + vi = 0.d0 + vr = 0.d0 + ! temperature = scf_temperature / hart + ! IF (realonly) THEN + ! vi = 0.d0 + ! ELSE + ! CALL imgw_v2(rs, xk, temperature, vi) + ! ENDIF + CASE(3) + ! Dirac-Hara + HL imag part + constant imaginary + CALL edp(rs, xk, vr) + CALL imhl(rs, xk, vi, icusp) + CASE(5) + ! Partially nonlocal: Dirac-Fock for core + HL for valence electrons + constant imaginary + ! TODO: Why is this the same as HL (case 0) ? Update documentation ? + CALL rhl(rs, xk, vr, vi) + CASE(6) + ! Finite temperature RPA GW + temperature = scf_temperature / hart + CALL rgw(rs, xk, temperature, vr) + IF (realonly) THEN + vi = 0.d0 + ELSE + CALL imgw_v2(rs, xk, temperature, vi) + ENDIF + CASE(7) + ! Temperature independent RPA GW (For testing purposes) + temperature = 0.d0 + CALL rgw(rs, xk, temperature, vr) + IF (realonly) THEN + vi = 0.d0 + ELSE + CALL imgw_v2(rs, xk, temperature, vi) + ENDIF + CASE(9) + ! TODO: What is this ? + ! This is the original ixc > 6 + ! IF (ixc.GE.6) THEN + CALL edp(rscore, xk, vrp) + vr = vr - vrp + ! ENDIF + CASE(10) + ! Same as ixc=0 with broadened plasmon HL selfenergy + ! using interpolation for the real and imaginary part. + CALL rhlbp(rs, xk, vr, vi) + CASE(13) + ! Same as ixc=3 with broadened plasmon HL selfenergy + ! using interpolation for the real and imaginary part. + CALL edp(rs, xk, vr) + CALL imhl(rs, xk, vi, icusp) + CASE(15) + ! Same as ixc=5 with broadened plasmon HL selfenergy + ! using interpolation for the real and imaginary part. + CALL rhlbp (rs, xk, vr, vi) + CASE DEFAULT + vr = 0.d0 + vi = 0.d0 + END SELECT + ! original script + ! if ((ixc.eq.0 .or. ixc.ge.5) .and. ibp .eq. 0) then + ! call rhl (rs, xk, vr, vi) + ! elseif ((ixc.eq.0.or. ixc.ge.5) .and. ibp .eq. 1) then + ! call rhlbp (rs, xk, vr, vi) + ! elseif (ixc .eq. 1) then + ! vi = 0 + ! call edp(rs,xk,vr) + ! elseif (ixc .eq. 3) then + ! call edp(rs,xk,vr) + ! call imhl (rs,xk,vi,icusp) + ! endif + + ! if (ixc .ge. 6) then + ! call edp(rscore,xk,vrp) + ! vr = vr - vrp + ! endif + + ! call gw_exact(rs, xk, vr, vi) + RETURN +END SUBROUTINE sigma + + diff --git a/src/FF2X/ff2xmu.f90 b/src/FF2X/ff2xmu.f90 index 44d7ffc..09a4d62 100644 --- a/src/FF2X/ff2xmu.f90 +++ b/src/FF2X/ff2xmu.f90 @@ -192,247 +192,262 @@ subroutine ff2xmu (ispec, ipr4, idwopt, critcw, s02, sig2g, & ! make combined title ntitle=ntitle-nheadold - do 120 ihead = 1, nhead - 120 title(ntitle+ihead) = head(ihead) - ntitle = ntitle + nhead - nheadold=nhead - - ! write feffnnnn.dat - if (ipr4.ge.3) then - call feffdt(ntotal,ip,nptot,ntitle,title,ne,npot, & - & ihole, iorder, ilinit, rnrmav, xmu, edge, potlbl, & - & iz,phc,ck,xk,index, & - & nleg,deg,nepts,reff,crit,ipot,rat,achi,phchi) - end if - - ! If there is a vicorr, will need a mean free path factor xlam0. - ! Use it as chi(ie) * exp (2 * reff * xlam0) - ! ckp is ck' = ck prime. - if (abs(vicorr) .ge. eps4) then - do 180 ie = 1, ne - ckp = sqrt (ck(ie)**2 + coni*2*vicorr) - xlam0 = aimag(ck(ie)) - dimag(ckp) - do 170 ipath = 1, nptot - achi(ie,ipath) = achi(ie,ipath) * exp (2 * reff(ipath) * xlam0) - 170 continue - 180 continue - endif - - ! k'**2 = k**2 + vr. If there is no real correction - ! (vrcorr = 0), these two grids will be the same. - ! k' is value for output, k is value used for - ! interpolations with original grid. - ! edge = -17.918/hart ! 1000K - ! edge = -17.937/hart ! 300K - ! vrcorr shifts the edge and the k grid - if (abs(vrcorr) .gt. eps4) edge = edge - vrcorr - - - ! ik0 is index at fermi level - do i = 1, ne - temp = xk(i)*abs(xk(i)) + 2*vrcorr - if (temp.ge. 0) then - xkp(i) = sqrt(temp) - else - xkp(i) = - sqrt(-temp) - endif - enddo - - - dwcorr = .false. - if (tk .gt. 1.0e-3) dwcorr = .true. - - ! Open chi.dat and xmu.dat (output) and start headers - if (iabs.eq.nabs) then - open (unit=3, file=f1, status='unknown', iostat=ios) !KJ changed chi.dat to f1 - call chopen (ios, f1, 'ff2chi') !KJ id. 1-06 - open (unit=8, file=f2, status='unknown', iostat=ios) !KJ changed xmu.dat to f2 - call chopen (ios, f2, 'ff2chi') !KJ id. - - ! write miscellaneous staff into headers !KJ corrected typo - call wrhead (8, ntitle, title, dwcorr, s02, & - & tk, thetad, sig2g, alphat, vrcorr, vicorr, critcw) - - ! also write information on the screen - if (alphat .gt. zero) then - write(slog,322) alphat - 322 format (' 1st and 3rd cumulants, alphat = ', 1pe20.4) - if(write_to_screen) call wlog(slog) - endif - if (abs(vrcorr).ge.eps4 .or. abs(vicorr).ge.eps4) then - write(slog,343) vrcorr*hart, vicorr*hart - 343 format (' Energy zero shift, vr, vi ', 1p, 2e14.5) - if(write_to_screen) call wlog(slog) - endif - - write(slog,370) critcw - if(write_to_screen) call wlog(slog) - 370 format (' Use all paths with cw amplitude ratio', f7.2, '%') - if (dwcorr) then - write(slog,380) s02, tk, thetad, sig2g - if(write_to_screen) call wlog(slog) - else - write(slog,381) s02, sig2g - if(write_to_screen) call wlog(slog) - endif - 380 format(' S02', f7.3, ' Temp', f8.2, ' Debye temp', f8.2, & - & ' Global sig2', f9.5) - 381 format(' S02', f7.3, ' Global sig2', f9.5) - endif - - - ! make chi and sum it - cchi(:) = 0 - do ik = 1, ne - cchi(ik)= s02 * gtr(ik) - enddo - - - ! add Debye-Waller factors - call dwadd (ntotal, nptot, idwopt, ip, index, crit, critcw, sig2g,& - & sig2u, dwcorr, rnrmav, nleg, deg, reff, iz, ipot, rat,tk,thetad,& - & alphat, thetae, mbconv, s02, ne1, ck, achi, phchi, ne, xk, xkp, & - & xkp, cchi, iabs, nabs, ispec, ipr4, ntitle, & - & title, vrcorr, vicorr, nused) - - ! read or initialize chia - result of configuration average - if (iabs.eq.1) then - chia(:) = 0 - else - open (unit=1, file='chia.bin', status='old', access='sequential', form='unformatted', iostat=ios) - do ie = 1,ne ; read(1) chia(ie) ; enddo - close (unit=1, status='delete') - endif - - if(iabs.eq.1) then - ! compare grids in xsect.bin and feff.bin - do i = 1, nxsec - del = xk(i)**2 - xkxs(i)**2 - ! JK - added extra test below because very large k fails the first test even if points are relatively very close - if (abs(del) .gt. eps4 .and. (abs(del)/xk(i)**2) .gt.eps4) then - call wlog(' Emesh in feff.bin and xsect.bin different.') - call par_stop('FF2XMU-1') - endif - enddo - endif - - ! add contribution from an absorber iabs - ! present scheme assumes that xsec is the same for all iabs. - do ik = 1, ne - chia(ik) = chia(ik) + cchi(ik)/ nabs - enddo - if (iabs.lt.nabs) then - ! save chia in chia.bin for averaging - open (unit=1, file='chia.bin', status='unknown', access='sequential', form='unformatted', iostat=ios) - do ie=1,ne ; write(1) chia(ie) ; enddo - close(unit=1) - endif - - if (iabs.eq.nabs) then - ! The loop over absorbers is finished. Write out the results. - write(8,600) coment, nused, ntotal - 600 format ( a2, 1x, i4, '/', i4, ' paths used') - 610 format ( a2, 1x, 71('-')) - - - rchtot(:) = dimag (chia(:)) ! = Im gtr - ! prepare the output grid omega - efermi = edge + omega(1) - dble(emxs(1)) - - ! do convolution with excitation spectrum - if (mbconv .gt. 0) then - wp = wp / 2. - call exconv (omega, ne1, efermi, s02p, erelax, wp, xsnorm) ! from xsect.bin - call exconv (omega, ne1, efermi, s02p, erelax, wp, rchtot) ! Im gtr - endif - - ! normalize to xsec at 50 eV above edge - ! and prepare the output energy grid omega - ! Tun: modify normalization at 100+kT eV above edge - ! edg50 = efermi + 50/hart - edg50 = efermi + 50/hart - if (ispec.eq.2) edg50 = efermi - call terp (omega, xsnorm, ne1, 1, edg50, xsedge) - ! IF (electronic_temperature.GT.0.d0) CALL terp (omega, xsnorm, ne1/2, 1, edg50, xsedge) - if (absolu.eq.1) xsedge=dble(1) !KJ 3-06 don't normalize - write(8,660) coment, xsedge - 660 format (a2, ' xsedge+ 50, used to normalize mu ', 1pe20.4) - write(8,610) coment - write(8,670) coment - 670 format (a2,' omega e k mu mu0 chi @#') - - if (.not.cross) then !KJ I added this block 1-06 - kxsec(:)=xsec(:) ! from xsect.bin - else - kxsec(:)=dcmplx(0,0) - endif !KJ end - - ! So, at this point, all we have is kxsec and xsnorm (from xsect.bin) and rchtot = Im gtr - - ! do correction using brouder method - vi0 = 0 - - write(slog,'(a,f10.5,a)') "electronic_temperature = ", electronic_temperature, " (eV)" - call wlog(slog) - IF (electronic_temperature.GT.0) THEN - call thermal_xscorr(ispec,emxs, ne1, ne, ik0, kxsec,xsnorm,chia,vrcorr, vi0, cchi, electronic_temperature) !KJ changed xsec to kxsec 1-06 - ELSE - call xscorr(ispec,emxs, ne1, ne, ik0, kxsec,xsnorm,chia,vrcorr, vi0, cchi) - ENDIF - - do ie=1,ne1 - rchtot(ie)=dimag(kxsec(ie)+xsnorm(ie)*chia(ie)+cchi(ie)) !KJ id. - enddo - - chia(:) = 0 - - if (cross) then !KJ bugfix 05/2010 to avoid NaN for crossterms (kxsec=chia=0) - cchi(:)=(0.d0,0.d0) - else - IF (electronic_temperature.GT.0) THEN ! Atomic part - call thermal_xscorr(ispec,emxs, ne1, ne, ik0, kxsec,xsnorm,chia,vrcorr, vi0, cchi, electronic_temperature) !KJ changed xsec to kxsec 1-06 - ELSE - call xscorr(ispec, emxs, ne1, ne, ik0, kxsec,xsnorm,chia,vrcorr, vi0, cchi) - ENDIF - endif - - do ie = 1, ne1 - cchi(ie) = dimag(kxsec(ie)+cchi(ie)) * coni+rchtot(ie) !KJ id. - enddo - - if (vicorr.gt.eps4 .and. ntotal.eq.0) then - ! add correction due to vicorr - call conv(omega,cchi,ne1,vicorr) - ! call conv(omega,xsec,ne1,vicorr) - endif - - do ie = 1, ne1 - em0 = dble(emxs(ie)) - xsec0 = dimag(cchi(ie)) - rchtot(ie) = dble (cchi(ie)) - chi0 = (rchtot(ie) - xsec0) - if (siggk) then - dwcoeff = exp(-(sig_gk*xkp(ie)/bohr)**2) !KJ energy dependent global DW factor for SIGGK card - chi0 = chi0 * dwcoeff - rchtot(ie) = xsec0 + chi0 - endif - write(8,700) omega(ie)*hart, em0*hart, xkp(ie)/bohr, rchtot(ie)/xsedge, xsec0/xsedge, chi0/xsedge + do ihead = 1, nhead + title(ntitle+ihead) = head(ihead) + enddo + + ntitle = ntitle + nhead + nheadold=nhead + + ! write feffnnnn.dat + if (ipr4.ge.3) then + call feffdt(ntotal,ip,nptot,ntitle,title,ne,npot, & + & ihole, iorder, ilinit, rnrmav, xmu, edge, potlbl, & + & iz,phc,ck,xk,index, & + & nleg,deg,nepts,reff,crit,ipot,rat,achi,phchi) + end if + + ! If there is a vicorr, will need a mean free path factor xlam0. + ! Use it as chi(ie) * exp (2 * reff * xlam0) + ! ckp is ck' = ck prime. + if (abs(vicorr) .ge. eps4) then + do ie = 1, ne + ckp = sqrt (ck(ie)**2 + coni*2*vicorr) + xlam0 = aimag(ck(ie)) - dimag(ckp) + do ipath = 1, nptot + achi(ie,ipath) = achi(ie,ipath) * exp (2 * reff(ipath) * xlam0) + enddo + enddo + endif + + ! k'**2 = k**2 + vr. If there is no real correction + ! (vrcorr = 0), these two grids will be the same. + ! k' is value for output, k is value used for + ! interpolations with original grid. + ! edge = -17.918/hart ! 1000K + ! edge = -17.937/hart ! 300K + ! vrcorr shifts the edge and the k grid + if (abs(vrcorr) .gt. eps4) edge = edge - vrcorr + + + ! ik0 is index at fermi level + do i = 1, ne + temp = xk(i)*abs(xk(i)) + 2*vrcorr + if (temp.ge. 0) then + xkp(i) = sqrt(temp) + else + xkp(i) = - sqrt(-temp) + endif + enddo + + + dwcorr = .false. + if (tk .gt. 1.0e-3) dwcorr = .true. + + ! Open chi.dat and xmu.dat (output) and start headers + if (iabs.eq.nabs) then + open (unit=3, file=f1, status='unknown', iostat=ios) !KJ changed chi.dat to f1 + call chopen (ios, f1, 'ff2chi') !KJ id. 1-06 + open (unit=8, file=f2, status='unknown', iostat=ios) !KJ changed xmu.dat to f2 + call chopen (ios, f2, 'ff2chi') !KJ id. + + ! write miscellaneous staff into headers !KJ corrected typo + call wrhead (8, ntitle, title, dwcorr, s02, & + & tk, thetad, sig2g, alphat, vrcorr, vicorr, critcw) + + ! also write information on the screen + if (alphat .gt. zero) then + write(slog,322) alphat + 322 format (' 1st and 3rd cumulants, alphat = ', 1pe20.4) + if(write_to_screen) call wlog(slog) + endif + if (abs(vrcorr).ge.eps4 .or. abs(vicorr).ge.eps4) then + write(slog,343) vrcorr*hart, vicorr*hart + 343 format (' Energy zero shift, vr, vi ', 1p, 2e14.5) + if(write_to_screen) call wlog(slog) + endif + + write(slog,370) critcw + if (write_to_screen) call wlog(slog) + 370 format (' Use all paths with cw amplitude ratio', f7.2, '%') + if (dwcorr) then + write(slog,380) s02, tk, thetad, sig2g + if(write_to_screen) call wlog(slog) + else + write(slog,381) s02, sig2g + if(write_to_screen) call wlog(slog) + endif + 380 format(' S02', f7.3, ' Temp', f8.2, ' Debye temp', f8.2, & + & ' Global sig2', f9.5) + 381 format(' S02', f7.3, ' Global sig2', f9.5) + endif + + + ! make chi and sum it + cchi(:) = 0 + do ik = 1, ne + cchi(ik)= s02 * gtr(ik) + enddo + + + ! add Debye-Waller factors + call dwadd (ntotal, nptot, idwopt, ip, index, crit, critcw, sig2g,& + & sig2u, dwcorr, rnrmav, nleg, deg, reff, iz, ipot, rat,tk,thetad,& + & alphat, thetae, mbconv, s02, ne1, ck, achi, phchi, ne, xk, xkp, & + & xkp, cchi, iabs, nabs, ispec, ipr4, ntitle, & + & title, vrcorr, vicorr, nused) + + ! read or initialize chia - result of configuration average + if (iabs.eq.1) then + chia(:) = 0 + else + open (unit=1, file='chia.bin', status='old', access='sequential', form='unformatted', iostat=ios) + do ie = 1,ne + read(1) chia(ie) + enddo + close (unit=1, status='delete') + endif + + if(iabs.eq.1) then + ! compare grids in xsect.bin and feff.bin + do i = 1, nxsec + del = xk(i)**2 - xkxs(i)**2 + ! JK - added extra test below because very large k fails the first test even if points are relatively very close + if (abs(del) .gt. eps4 .and. (abs(del)/xk(i)**2) .gt.eps4) then + call wlog(' Emesh in feff.bin and xsect.bin different.') + call par_stop('FF2XMU-1') + endif + enddo + endif + + ! add contribution from an absorber iabs + ! present scheme assumes that xsec is the same for all iabs. + do ik = 1, ne + chia(ik) = chia(ik) + cchi(ik)/ nabs + enddo + if (iabs.lt.nabs) then + ! save chia in chia.bin for averaging + open (unit=1, file='chia.bin', status='unknown', access='sequential', form='unformatted', iostat=ios) + do ie=1,ne + write(1) chia(ie) + enddo + close(unit=1) + endif + + if (iabs.eq.nabs) then + ! The loop over absorbers is finished. Write out the results. + write(8,600) coment, nused, ntotal + 600 format ( a2, 1x, i4, '/', i4, ' paths used') + 610 format ( a2, 1x, 71('-')) + + rchtot(:) = dimag (chia(:)) ! = Im gtr + ! prepare the output grid omega + efermi = edge + omega(1) - dble(emxs(1)) + + write(8, 680) coment, efermi*hart + write(8, 690) coment, electronic_temperature + 680 FORMAT(a2,' Chemical potential = ', 1pe20.8, ' eV') + 690 FORMAT(a2,' Electronic temperature = ', 1pe20.8, ' eV') + + ! do convolution with excitation spectrum + if (mbconv .gt. 0) then + wp = wp / 2. + call exconv (omega, ne1, efermi, s02p, erelax, wp, xsnorm) ! from xsect.bin + call exconv (omega, ne1, efermi, s02p, erelax, wp, rchtot) ! Im gtr + endif + + ! normalize to xsec at 50 eV above edge + ! and prepare the output energy grid omega + ! Tun: modify normalization at 100+kT eV above edge + ! edg50 = efermi + 50/hart + edg50 = efermi + 50/hart + if (ispec.eq.2) edg50 = efermi + call terp (omega, xsnorm, ne1, 1, edg50, xsedge) + ! IF (electronic_temperature.GT.0.d0) CALL terp (omega, xsnorm, ne1/2, 1, edg50, xsedge) + if (absolu.eq.1) xsedge=dble(1) !KJ 3-06 don't normalize + write(8,660) coment, xsedge + 660 format (a2, ' xsedge+ 50, used to normalize mu ', 1pe20.4) + write(8,610) coment + write(8,670) coment + 670 format (a2,' omega e k mu mu0 chi @#') + + if (.not.cross) then !KJ I added this block 1-06 + kxsec(:)=xsec(:) ! from xsect.bin + else + kxsec(:)=dcmplx(0,0) + endif !KJ end + + ! So, at this point, all we have is kxsec and xsnorm (from xsect.bin) and rchtot = Im gtr + + ! do correction using brouder method + vi0 = 0 + + write(slog,'(a,f10.5,a)') "electronic_temperature = ", electronic_temperature, " (eV)" + call wlog(slog) + IF (electronic_temperature.GT.0) THEN + call thermal_xscorr(ispec,emxs, ne1, ne, ik0, kxsec,xsnorm,chia,vrcorr, vi0, cchi, electronic_temperature) !KJ changed xsec to kxsec 1-06 + ELSE + call xscorr(ispec,emxs, ne1, ne, ik0, kxsec,xsnorm,chia,vrcorr, vi0, cchi) + ENDIF + + OPEN(UNIT=1149, file="chia.dat", status='REPLACE') + do ie = 1, ne + WRITE(1149, '(20E20.10E3)') DBLE(chia(ie)), DIMAG(chia(ie)) + enddo + close(1149) + + OPEN(UNIT=1149, file="rchtot1.dat", status='REPLACE') + do ie=1,ne1 + rchtot(ie)=dimag(kxsec(ie)+xsnorm(ie)*chia(ie)+cchi(ie)) !KJ id. + WRITE(1149, '(20E20.10E3)') DBLE(omega(ie))*hart, dimag(kxsec(ie)), dimag(xsnorm(ie)*chia(ie)), dimag(cchi(ie)), rchtot(ie) + enddo + close(1149) + chia(:) = 0 + + if (cross) then !KJ bugfix 05/2010 to avoid NaN for crossterms (kxsec=chia=0) + cchi(:)=(0.d0,0.d0) + else + IF (electronic_temperature.GT.0) THEN ! Atomic part + call thermal_xscorr(ispec,emxs, ne1, ne, ik0, kxsec,xsnorm,chia,vrcorr, vi0, cchi, electronic_temperature) !KJ changed xsec to kxsec 1-06 + ELSE + call xscorr(ispec, emxs, ne1, ne, ik0, kxsec,xsnorm,chia,vrcorr, vi0, cchi) + ENDIF + endif + + do ie = 1, ne1 + cchi(ie) = dimag(kxsec(ie)+cchi(ie)) * coni+rchtot(ie) !KJ id. + enddo + + if (vicorr.gt.eps4 .and. ntotal.eq.0) then + ! add correction due to vicorr + call conv(omega,cchi,ne1,vicorr) + ! call conv(omega,xsec,ne1,vicorr) + endif + + do ie = 1, ne1 + em0 = dble(emxs(ie)) + xsec0 = dimag(cchi(ie)) + rchtot(ie) = dble (cchi(ie)) + chi0 = (rchtot(ie) - xsec0) + if (siggk) then + dwcoeff = exp(-(sig_gk*xkp(ie)/bohr)**2) !KJ energy dependent global DW factor for SIGGK card + chi0 = chi0 * dwcoeff + rchtot(ie) = xsec0 + chi0 + endif + write(8,700) omega(ie)*hart, em0*hart, xkp(ie)/bohr, rchtot(ie)/xsedge, xsec0/xsedge, chi0/xsedge ! if you want f'' at the output in el. units use next line ! rchtot(ie)*omega(ie)*prefac, xsec0*omega(ie)*prefac, chi0 ! with prefac = alpinv / 4 / pi /bohr**2 - 700 format (1x, 3f20.10, 1p, 3e20.10) + 700 format (1x, 3f20.10, 1p, 3e20.10) - enddo + enddo - close (unit=8) - close (unit=3, status='delete') - endif + close (unit=8) + close (unit=3, status='delete') + endif ! for if (iabs=abs); or the last absorber - - enddo !KJ of my iip=ipmin,ipmax,ipstep loop 1-06 - - - return - end + enddo !KJ of my iip=ipmin,ipmax,ipstep loop 1-06 + return +end diff --git a/src/FF2X/m_thermal_xscorr.f90 b/src/FF2X/m_thermal_xscorr.f90 index d23155b..a5ba47e 100644 --- a/src/FF2X/m_thermal_xscorr.f90 +++ b/src/FF2X/m_thermal_xscorr.f90 @@ -33,7 +33,7 @@ SUBROUTINE thermal_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & COMPLEX*16 :: ec(nex), fc(nex), xmu(nex), pole COMPLEX*16 :: xmu1(nex), xmu2(nex), xmu3(nex), xmu4(nex) REAL(8), PARAMETER :: eps4 = 1.0d-4 - LOGICAL, PARAMETER :: print_out = .TRUE. + LOGICAL, PARAMETER :: print_out = .FALSE. EXTERNAL :: lorenz, astep COMPLEX*16 :: lorenz REAL(8) :: astep @@ -58,6 +58,7 @@ SUBROUTINE thermal_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & ne4 = ne1 efermi = DBLE(emxs(ne-extra)) + xloss = DIMAG(emxs(ne4+1)) cheight = DIMAG(emxs(1)) @@ -67,25 +68,25 @@ SUBROUTINE thermal_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & xmu(ie) = xsec(ie) + xsnorm(ie)*chia(ie) ENDDO - - OPEN(UNIT=155, file='curve.dat', status='REPLACE') - OPEN(UNIT=150, file="prexmu.dat", status='REPLACE') - OPEN(UNIT=149, file="raw.dat", status='REPLACE') - OPEN(UNIT=153, file="leg1.dat", status='REPLACE') - OPEN(UNIT=151, file='residue.dat', status='REPLACE') - OPEN(UNIT=152, file='contour.dat', status='REPLACE') - OPEN(UNIT=1543, file='ratio.dat', status='REPLACE') - WRITE(149,*) "Temperature (Hatree)", tk - WRITE(149,*) "Electronic Temperature (eV)", electronic_temperature - WRITE(149,*) "xloss = ", xloss, " Hatree" - WRITE(149,*) "Chemical potential = ", efermi*hart, " eV with shift", vrcorr*hart - WRITE(149,*) "Number of poles = ", npole - WRITE(149, *) 'Omega(Hart) \t Re CCHI \t Im CCHI \t 1-Fermi \t Re xmu0 \t Im xmu0' - - DO ie =1, ne - WRITE(155, '(I4,20E20.10E3)') ie, DBLE(emxs(ie)), DIMAG(emxs(ie)) - ENDDO - CLOSE(155) + IF (print_out) THEN + ! OPEN(UNIT=155, file='curve.dat', status='REPLACE') + OPEN(UNIT=150, file="prexmu.dat", status='REPLACE') + OPEN(UNIT=149, file="raw.dat", status='REPLACE') + OPEN(UNIT=153, file="leg1.dat", status='REPLACE') + ! OPEN(UNIT=151, file='residue.dat', status='REPLACE') + OPEN(UNIT=152, file='contour.dat', status='REPLACE') + ! OPEN(UNIT=1543, file='ratio.dat', status='REPLACE') + WRITE(149,*) "Temperature (Hatree)", tk + WRITE(149,*) "Electronic Temperature (eV)", electronic_temperature + WRITE(149,*) "xloss = ", xloss, " Hatree" + WRITE(149,*) "Chemical potential = ", efermi*hart, " eV with shift", vrcorr*hart + WRITE(149,*) "Number of poles = ", npole + WRITE(149, *) 'Omega(Hart) \t Re CCHI \t Im CCHI \t 1-Fermi \t Re xmu0 \t Im xmu0' + ENDIF + ! DO ie =1, ne + ! WRITE(155, '(I4,20E20.10E3)') ie, DBLE(emxs(ie)), DIMAG(emxs(ie)) + ! ENDDO + ! CLOSE(155) !---------- Construct the integration contour @@ -137,7 +138,7 @@ SUBROUTINE thermal_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & IF (ispec.ne.2) cchi(ie) = xmu0 - cchi(ie) ! Absorption dummy = 1-dummy z1 = omega(ie) - WRITE(149, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)), DBLE(dummy), DBLE(xmu0), DIMAG(xmu0), DBLE(fermiDirac(z1,tk,efermi)) + ! WRITE(149, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)), DBLE(dummy), DBLE(xmu0), DIMAG(xmu0), DBLE(fermiDirac(z1,tk,efermi)) !!!!!!!!!!!!!!!!!!!!!! !! RESIDUE of Fermi !! @@ -154,7 +155,7 @@ SUBROUTINE thermal_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & residue = -residue*(2.d0*coni*tk) IF (ispec.ne.2) residue = - residue ! Absorption - WRITE(151, '(20E20.10E3)') DBLE(omega(ie)), DBLE(residue), DIMAG(residue), DBLE(efermi) + !WRITE(151, '(20E20.10E3)') DBLE(omega(ie)), DBLE(residue), DIMAG(residue), DBLE(efermi) !!!!!!!!!!!!!!! !! First leg !! @@ -162,7 +163,7 @@ SUBROUTINE thermal_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & ! Integrate from 0 to cheight with x in R-space. corr1 = 0.d0 - WRITE(153, '(20E20.10E3)') DBLE(omega(ie)), DBLE(corr1), DIMAG(corr1) + IF (print_out) WRITE(153, '(20E20.10E3)') DBLE(omega(ie)), DBLE(corr1), DIMAG(corr1) !!!!!!!!!!!!!!!!! !! Second leg !! @@ -171,36 +172,41 @@ SUBROUTINE thermal_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & CALL cchi_eff(emxs(1:ne4), xmu2, ne4, efermi, tk, xmu0, xloss, omega(ie), ispec, corr2) ! CALL fermilorenH(emxs(1:ne4), ne4, efermi, tk, xloss, omega(ie), ispec, dummy2) ! CALL fermilorenTH(omega(ie),tk,efermi,xloss,cheight,dummy2) - ! WRITE(152, '(20E20.10E3)') DBLE(omega(ie)), DBLE(corr2), DIMAG(corr2), DBLE(xmu2(ie)-xmu0), DIMAG(xmu2(ie)-xmu0) + IF (print_out) WRITE(152, '(20E20.10E3)') DBLE(omega(ie)), DBLE(corr2), DIMAG(corr2), DBLE(xmu2(ie)-xmu0), DIMAG(xmu2(ie)-xmu0), xmu0 ! WRITE(1543, '(20E20.10E3)') DBLE(omega(ie)), DBLE(corr2), DIMAG(corr2), & ! & DBLE(cchi(ie)), DIMAG(cchi(ie)), & ! & DBLE(dummy), DBLE(xmu0), DIMAG(xmu0), DBLE(dummy2), DIMAG(dummy2) - WRITE(1543, '(20E20.10E3)') DBLE(omega(ie)), DBLE(1-fermiDirac(omega(ie)+coni*0.d0,tk,efermi)), & - & DBLE(dummy),& - & DBLE(cchi(ie)), DIMAG(cchi(ie)), & - & DBLE(corr2), DIMAG(corr2), & - & DBLE(cauchy(omega(ie)+coni*0.0d0,efermi,xloss)), & - & DIMAG(xmu2(ie)-xmu0) + ! WRITE(1543, '(20E20.10E3)') DBLE(omega(ie)), DBLE(1-fermiDirac(omega(ie)+coni*0.d0,tk,efermi)), & + ! & DBLE(dummy),& + ! & DBLE(cchi(ie)), DIMAG(cchi(ie)), & + ! & DBLE(corr2), DIMAG(corr2), & + ! & DBLE(cauchy(omega(ie)+coni*0.0d0,efermi,xloss)), & + ! & DIMAG(xmu2(ie)-xmu0) !!!!!!!!!!!!!!!!!!!!! !! Combine results !! !!!!!!!!!!!!!!!!!!!!! corr = corr1 + corr2 cchi(ie) = cchi(ie) + corr + residue - WRITE(150, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)) + ! WRITE(149, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)) ! Return the result of convolution minus bare value cchi(ie) = cchi(ie) - xmu2(ie) + ! WRITE(150, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)) + IF (print_out) WRITE(149, '(20E20.10E3)') DBLE(omega(ie)*hart), DIMAG(cchi(ie)), DIMAG(xsec(ie)), DIMAG(xsnorm(ie)*chia(ie)), DIMAG(xsec(ie)+xsnorm(ie)*chia(ie)+cchi(ie)) + + ENDDO !------- End of convolution - CLOSE(149) - CLOSE(150) - CLOSE(151) - CLOSE(152) - CLOSE(153) - CLOSE(1543) - + IF (print_out) THEN + CLOSE(149) + CLOSE(150) + ! CLOSE(151) + CLOSE(152) + CLOSE(153) + ! CLOSE(1543) + ENDIF RETURN - END SUBROUTINE + END SUBROUTINE thermal_xscorr SUBROUTINE fermilorenT(omega, tk, efermi, xloss, result) ! Perform convolution between fermi function @@ -539,8 +545,8 @@ SUBROUTINE cchi_eff(emxs, xmu2, ne4, efermi, tk, xmu0, xloss, omega, ispec, resu ! Find where does xmu +/- window_terp lies on the contour ind_m3 = binarysearch(DBLE(efermi_m3), DBLE(emxs), ne4)-1 ind_p3 = binarysearch(DBLE(efermi_p3), DBLE(emxs), ne4) - - IF (ind_m3.EQ.0) BUG=.TRUE. + ! IF ((ind_m3.EQ.0).AND.(ABS( DBLE(efermi_m3) - DBLE(emxs(1))).GT.1e-16)) BUG=.TRUE. + IF ((ind_m3.EQ.0).AND.(.NOT.(DBLE(efermi_m3).EQ.DBLE(emxs(1))))) BUG=.TRUE. IF (BUG) THEN ! Need to be careful regarding the window_size_terp ! efermi_m3 must be greater than emxs(1) @@ -738,7 +744,7 @@ SUBROUTINE sommerfeld_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & COMPLEX*16 :: ec(nex), fc(nex), xmu(nex), pole COMPLEX*16 :: xmu1(nex), xmu2(nex), xmu3(nex), xmu4(nex) REAL(8), PARAMETER :: eps4 = 1.0d-4 - LOGICAL, PARAMETER :: print_out = .TRUE. + LOGICAL, PARAMETER :: print_out = .FALSE. EXTERNAL :: lorenz, astep COMPLEX*16 :: lorenz, astep @@ -771,24 +777,26 @@ SUBROUTINE sommerfeld_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & xmu(ie) = xsec(ie) + xsnorm(ie)*chia(ie) ENDDO - OPEN(UNIT=155, file='curve.dat', status='REPLACE') - OPEN(UNIT=150, file="prexmu.dat", status='REPLACE') - OPEN(UNIT=149, file="raw.dat", status='REPLACE') - OPEN(UNIT=153, file="leg1.dat", status='REPLACE') - OPEN(UNIT=151, file='residue.dat', status='REPLACE') - OPEN(UNIT=152, file='contour.dat', status='REPLACE') - ! OPEN(UNIT=154, file='lorentzian.dat', status='REPLACE') - WRITE(149,*) "Temperature (Hatree)", tk - WRITE(149,*) "Electronic Temperature (eV)", electronic_temperature - WRITE(149,*) "xloss = ", xloss, " Hatree" - WRITE(149,*) "Chemical potential = ", efermi*hart, " eV" - WRITE(149,*) "Number of poles = ", npole - WRITE(149, *) 'Omega(Hart) \t Re CCHI \t Im CCHI \t 1-Fermi \t Re xmu0 \t Im xmu0' - - DO ie =1, ne - WRITE(155, '(20E20.10E3)') ie, DBLE(emxs(ie)), DIMAG(emxs(ie)) - ENDDO - CLOSE(155) + IF (print_out) THEN + OPEN(UNIT=155, file='curve.dat', status='REPLACE') + OPEN(UNIT=150, file="prexmu.dat", status='REPLACE') + OPEN(UNIT=149, file="raw.dat", status='REPLACE') + OPEN(UNIT=153, file="leg1.dat", status='REPLACE') + OPEN(UNIT=151, file='residue.dat', status='REPLACE') + OPEN(UNIT=152, file='contour.dat', status='REPLACE') + ! OPEN(UNIT=154, file='lorentzian.dat', status='REPLACE') + WRITE(149,*) "Temperature (Hatree)", tk + WRITE(149,*) "Electronic Temperature (eV)", electronic_temperature + WRITE(149,*) "xloss = ", xloss, " Hatree" + WRITE(149,*) "Chemical potential = ", efermi*hart, " eV" + WRITE(149,*) "Number of poles = ", npole + WRITE(149, *) 'Omega(Hart) \t Re CCHI \t Im CCHI \t 1-Fermi \t Re xmu0 \t Im xmu0' + + DO ie =1, ne + WRITE(155, '(20E20.10E3)') ie, DBLE(emxs(ie)), DIMAG(emxs(ie)) + ENDDO + CLOSE(155) + ENDIF !---------- Construct the integration contour nc = 0 ! DO ie = 1, ne2 @@ -840,7 +848,7 @@ SUBROUTINE sommerfeld_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & cchi(ie) = xmu0*cchi(ie) IF (ispec.ne.2) cchi(ie) = xmu0 - cchi(ie) ! Absorption dummy = 1-dummy - WRITE(149, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)), DBLE(dummy), DBLE(xmu0), DIMAG(xmu0) + IF (print_out) WRITE(149, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)), DBLE(dummy), DBLE(xmu0), DIMAG(xmu0) !!!!!!!!!!!!!!!!!!!!!! !! RESIDUE of Fermi !! @@ -856,7 +864,7 @@ SUBROUTINE sommerfeld_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & residue = -residue*(2.d0*coni*tk) IF (ispec.ne.2) residue = - residue ! Absorption - WRITE(151, '(20E20.10E3)') DBLE(omega(ie)), DBLE(residue), DIMAG(residue), DBLE(efermi) + IF (print_out) WRITE(151, '(20E20.10E3)') DBLE(omega(ie)), DBLE(residue), DIMAG(residue), DBLE(efermi) !!!!!!!!!!!!!!!! !! Correction !! @@ -880,7 +888,7 @@ SUBROUTINE sommerfeld_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & IF (ispec.eq.2) corr1 = -corr1 - WRITE(153, '(20E20.10E3)') DBLE(omega(ie)), DBLE(corr1), DIMAG(corr1) + IF (print_out) WRITE(153, '(20E20.10E3)') DBLE(omega(ie)), DBLE(corr1), DIMAG(corr1) !!!!!!!!!!!!!!!!! !! Second leg !! @@ -902,25 +910,27 @@ SUBROUTINE sommerfeld_xscorr(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & ENDDO IF (ispec.EQ.2) corr2 = -corr2 - WRITE(152, '(20E20.10E3)') DBLE(omega(ie)), DBLE(corr2), DIMAG(corr2) + IF (print_out) WRITE(152, '(20E20.10E3)') DBLE(omega(ie)), DBLE(corr2), DIMAG(corr2) !!!!!!!!!!!!!!!!!!!!! !! Combine results !! !!!!!!!!!!!!!!!!!!!!! corr = corr1 + corr2 cchi(ie) = cchi(ie) + corr + residue - WRITE(150, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)) + IF (print_out) WRITE(150, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)) ! Return the result of convolution minus bare value cchi(ie) = cchi(ie) - xmu2(ie) ENDDO !------- End of convolution - CLOSE(149) - CLOSE(150) - CLOSE(151) - CLOSE(152) - CLOSE(153) - ! CLOSE(154) + IF (print_out) THEN + CLOSE(149) + CLOSE(150) + CLOSE(151) + CLOSE(152) + CLOSE(153) + ! CLOSE(154) + ENDIF RETURN END SUBROUTINE @@ -944,6 +954,7 @@ SUBROUTINE sommerfeld_xscorr0(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & COMPLEX*16 :: corr, corr1, corr2, correction, leg3, residue COMPLEX*16 :: ec(nex), fc(nex), xmu(nex), pole REAL(8), PARAMETER :: eps4 = 1.0d-4 + LOGICAL, PARAMETER :: print_out = .FALSE. EXTERNAL :: lorenz, astep REAL(8) :: astep COMPLEX*16 :: lorenz @@ -967,21 +978,22 @@ SUBROUTINE sommerfeld_xscorr0(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & ! Only happens at low temperature IF (MOD(xloss, tk*pi).LT.eps4) PRINT*, "WARNING: xloss close to pole" - OPEN(UNIT=150, file="prexmu.dat", status='REPLACE') - OPEN(UNIT=149, file="raw.dat", status='REPLACE') - OPEN(UNIT=153, file="leg1.dat", status='REPLACE') - OPEN(UNIT=151, file='residue.dat', status='REPLACE') - OPEN(UNIT=152, file='contour.dat', status='REPLACE') - OPEN(UNIT=155, file='curve.dat', status='REPLACE') - OPEN(UNIT=156, file='correction.dat', status='REPLACE') - ! OPEN(UNIT=154, file='lorentzian.dat', status='REPLACE') - WRITE(149,*) "Temperature (Hatree)", tk - WRITE(149,*) "Electronic Temperature (eV)", electronic_temperature - WRITE(149,*) "xloss = ", xloss, " Hatree" - WRITE(149,*) "Chemical potential = ", efermi*hart, " eV" - WRITE(149,*) "Number of poles = ", npole - WRITE(149, *) 'Omega(Hart) \t Re CCHI \t Im CCHI \t 1-Fermi \t Re xmu0 \t Im xmu0' - + IF (print_out) THEN + OPEN(UNIT=150, file="prexmu.dat", status='REPLACE') + OPEN(UNIT=149, file="raw.dat", status='REPLACE') + OPEN(UNIT=153, file="leg1.dat", status='REPLACE') + OPEN(UNIT=151, file='residue.dat', status='REPLACE') + OPEN(UNIT=152, file='contour.dat', status='REPLACE') + OPEN(UNIT=155, file='curve.dat', status='REPLACE') + OPEN(UNIT=156, file='correction.dat', status='REPLACE') + ! OPEN(UNIT=154, file='lorentzian.dat', status='REPLACE') + WRITE(149,*) "Temperature (Hatree)", tk + WRITE(149,*) "Electronic Temperature (eV)", electronic_temperature + WRITE(149,*) "xloss = ", xloss, " Hatree" + WRITE(149,*) "Chemical potential = ", efermi*hart, " eV" + WRITE(149,*) "Number of poles = ", npole + WRITE(149, *) 'Omega(Hart) \t Re CCHI \t Im CCHI \t 1-Fermi \t Re xmu0 \t Im xmu0' + ENDIF DO ie = 1, ne xmu(ie) = xsec(ie) + xsnorm(ie)*chia(ie) ENDDO @@ -1016,7 +1028,7 @@ SUBROUTINE sommerfeld_xscorr0(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & nc = nc + 1 ec(nc) = emxs(ne1+ie) fc(nc) = xmu(ne1+ie) - WRITE(155, '(20E20.10E3)') DBLE(ec(nc)), DIMAG(ec(nc)), DBLE(fc(nc)), DIMAG(fc(nc)) + IF (print_out) WRITE(155, '(20E20.10E3)') DBLE(ec(nc)), DIMAG(ec(nc)), DBLE(fc(nc)), DIMAG(fc(nc)) ENDIF ENDDO @@ -1026,11 +1038,11 @@ SUBROUTINE sommerfeld_xscorr0(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & if (abs(vrcorr).gt.eps4) then ec(nc) = efermi + coni*xloss fc(nc) = bb * xmu(ik0) - WRITE(155, '(20E20.10E3)') DBLE(ec(nc)), DIMAG(ec(nc)), DBLE(fc(nc)), DIMAG(fc(nc)) + IF (print_out) WRITE(155, '(20E20.10E3)') DBLE(ec(nc)), DIMAG(ec(nc)), DBLE(fc(nc)), DIMAG(fc(nc)) else ec(nc) = emxs(ik0) fc(nc) = xmu(ik0) - WRITE(155, '(20E20.10E3)') DBLE(ec(nc)), DIMAG(ec(nc)), DBLE(fc(nc)), DIMAG(fc(nc)) + IF (print_out) WRITE(155, '(20E20.10E3)') DBLE(ec(nc)), DIMAG(ec(nc)), DBLE(fc(nc)), DIMAG(fc(nc)) endif IF (ispec.NE.2) THEN @@ -1039,7 +1051,7 @@ SUBROUTINE sommerfeld_xscorr0(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & nc = nc+1 ec(nc) = emxs(ie) fc(nc) = xmu(ie) - WRITE(155, '(20E20.10E3)') DBLE(ec(nc)), DIMAG(ec(nc)), DBLE(fc(nc)), DIMAG(fc(nc)) + IF (print_out) WRITE(155, '(20E20.10E3)') DBLE(ec(nc)), DIMAG(ec(nc)), DBLE(fc(nc)), DIMAG(fc(nc)) ENDIF ENDDO ELSE @@ -1048,7 +1060,7 @@ SUBROUTINE sommerfeld_xscorr0(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & nc = nc+1 ec(nc) = emxs(ie) fc(nc) = xmu(ie) - WRITE(155, '(20E20.10E3)') DBLE(ec(nc)), DIMAG(ec(nc)), DBLE(fc(nc)), DIMAG(fc(nc)) + IF (print_out) WRITE(155, '(20E20.10E3)') DBLE(ec(nc)), DIMAG(ec(nc)), DBLE(fc(nc)), DIMAG(fc(nc)) ENDIF ENDDO ENDIF @@ -1079,7 +1091,7 @@ SUBROUTINE sommerfeld_xscorr0(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & cchi(ie) = xmu0*cchi(ie) IF (ispec.ne.2) cchi(ie) = xmu0 - cchi(ie) ! Absorption dummy = 1-dummy - WRITE(149, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)), DBLE(dummy), DBLE(xmu0), DIMAG(xmu0) + IF (print_out) WRITE(149, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)), DBLE(dummy), DBLE(xmu0), DIMAG(xmu0) !!!!!!!!!!!!!!!!!!!!!!!!!!! !! sommerfeld correction !! @@ -1114,7 +1126,7 @@ SUBROUTINE sommerfeld_xscorr0(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & w1 = DIMAG(ec(1)) corr = corr + lorenz0(xloss,w1,dele)*ff(1)*coni*w1 leg3 = lorenz0(xloss,w1,dele) - WRITE(151, '(20E20.10E3)') DBLE(omega(ie)), DBLE(leg3), DIMAG(leg3) + IF (print_out) WRITE(151, '(20E20.10E3)') DBLE(omega(ie)), DBLE(leg3), DIMAG(leg3) ! Horizontal Contour DO ic = 1, nc-1 @@ -1134,8 +1146,8 @@ SUBROUTINE sommerfeld_xscorr0(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & IF (ispec.eq.2) corr = -corr cchi(ie) = cchi(ie) + corr + correction !+ correction + residue - WRITE(152, '(20E20.10E3)') DBLE(omega(ie)), DBLE(corr), DIMAG(corr) - WRITE(150, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)) + IF (print_out) WRITE(152, '(20E20.10E3)') DBLE(omega(ie)), DBLE(corr), DIMAG(corr) + IF (print_out) WRITE(150, '(20E20.10E3)') DBLE(omega(ie)), DBLE(cchi(ie)), DIMAG(cchi(ie)) @@ -1144,13 +1156,15 @@ SUBROUTINE sommerfeld_xscorr0(ispec, emxs, ne1, ne, ik0, xsec, xsnorm, chia, & cchi(ie) = cchi(ie) - xmu(ie) ENDDO !------- End of convolution - CLOSE(149) - CLOSE(150) - CLOSE(151) - CLOSE(152) - CLOSE(153) - CLOSE(155) - CLOSE(156) + IF (print_out) THEN + CLOSE(149) + CLOSE(150) + CLOSE(151) + CLOSE(152) + CLOSE(153) + CLOSE(155) + CLOSE(156) + ENDIF ! end of cycle over frequency points if (abs(vrcorr).gt.eps4) then diff --git a/src/LDOS/ldossub.f90 b/src/LDOS/ldossub.f90 index 50a83a0..c62bbe0 100644 --- a/src/LDOS/ldossub.f90 +++ b/src/LDOS/ldossub.f90 @@ -117,7 +117,8 @@ subroutine ldos ( nph, edens, edenvl, dmag, vtot, vvalgs, & call fixvar (rmt(iph),edens(1,iph),vtot(1,iph),dmag(1,iph), & & vint, rhoint, dx, rgrd, jumprm, & & vjump, ri, vtotph, dum, dmagx) - if (mod(ixc,10) .ge.5) then + !if (mod(ixc,10) .ge.5) then + if ((mod(ixc,10).ge.5).and.(ixc.ne.6).and.(ixc.ne.7)) then ! TTS if (jumprm .gt. 0) jumprm = 2 call fixvar (rmt(iph), edenvl(1,iph), vvalgs(1,iph), & & dmag(1,iph), vint, rhoint, dx, rgrd , jumprm, & @@ -132,7 +133,8 @@ subroutine ldos ( nph, edens, edenvl, dmag, vtot, vvalgs, & do i = 1, jri1 vtotph(i) = vtotph(i) - eref(1) enddo - if (ixc.ge.5) then + ! if (ixc.ge.5) then + if ((ixc.ge.5).and.(ixc.ne.6).and.(ixc.ne.7)) then ! TTS do i = 1, jri1 vvalph(i) = vvalph(i) - eref(1) enddo diff --git a/src/LDOS/ldossub_h.f90 b/src/LDOS/ldossub_h.f90 index 198c46e..d99b71f 100755 --- a/src/LDOS/ldossub_h.f90 +++ b/src/LDOS/ldossub_h.f90 @@ -123,7 +123,8 @@ subroutine ldos_h ( nph, edens_sp, edenvl_sp, dmag, vtot_sp, vvalgs_sp, & edens(:,iph)=edens_sp(:,iph,is) edenvl(:,iph)=edenvl_sp(:,iph,is) call fixvar (rmt(iph),edens(1,iph),vtot(1,iph),dmag(1,iph), vint, rhoint, dx, rgrd, jumprm, vjump, ri, vtotph, dum, dmagx) - if (mod(ixc,10) .ge.5) then + ! if (mod(ixc,10) .ge.5) then + if ((mod(ixc,10).ge.5).and.(ixc.ne.6).and.(ixc.ne.7)) then ! TTS if (jumprm .gt. 0) jumprm = 2 call fixvar (rmt(iph), edenvl(1,iph), vvalgs(1,iph), dmag(1,iph), vint, rhoint, dx, rgrd , jumprm, vjump, ri, vvalph, dum, dmagx) if (jumprm .gt. 0) jumprm = 1 @@ -141,7 +142,8 @@ subroutine ldos_h ( nph, edens_sp, edenvl_sp, dmag, vtot_sp, vvalgs_sp, & do i = 1, jri1 vtotph(i) = vtotph(i) - eref(1) enddo - if (ixc.ge.5) then + ! if (ixc.ge.5) then + if ((ixc.ge.5).and.(ixc.ne.6).and.(ixc.ne.7)) then ! TTS do i = 1, jri1 vvalph(i) = vvalph(i) - eref(1) enddo diff --git a/src/LDOS/ldossub_h_unrolled.f90 b/src/LDOS/ldossub_h_unrolled.f90 index 04fca10..5ed3bd7 100644 --- a/src/LDOS/ldossub_h_unrolled.f90 +++ b/src/LDOS/ldossub_h_unrolled.f90 @@ -112,7 +112,8 @@ subroutine ldos_h_unrolled ( nph, edens_sp, edenvl_sp, dmag, vtot_sp, vvalgs_sp, call fixvar (rmt(iph),edens(1,iph),vtot(1,iph),dmag(1,iph), & vint, rhoint, dx, rgrd, jumprm, vjump, ri, & vtotph, dum, dmagx) - if (mod(ixc,10) .ge.5) then + ! if (mod(ixc,10) .ge.5) then + if ((mod(ixc,10).ge.5).and.(ixc.ne.6).and.(ixc.ne.7)) then ! TTS if (jumprm .gt. 0) jumprm = 2 call fixvar (rmt(iph), edenvl(1,iph), vvalgs(1,iph), & dmag(1,iph), vint, rhoint, dx, rgrd , & @@ -128,7 +129,8 @@ subroutine ldos_h_unrolled ( nph, edens_sp, edenvl_sp, dmag, vtot_sp, vvalgs_sp, eref_sp(1:ne,is) = eref(1) eref(1:ne) = eref(1) vtotph(1:jri1) = vtotph(1:jri1) - eref(1) - if (ixc.ge.5) then + ! if (ixc.ge.5) then + if ((ixc.ge.5).and.(ixc.ne.6).and.(ixc.ne.7)) then ! TTS vvalph(1:jri1) = vvalph(1:jri1) - eref(1) else vvalph(1:jri1) = vtotph(1:jri1) diff --git a/src/LDOS/rhol.f90 b/src/LDOS/rhol.f90 index 139449a..285d2e6 100644 --- a/src/LDOS/rhol.f90 +++ b/src/LDOS/rhol.f90 @@ -25,6 +25,8 @@ subroutine rhol ( dx, x0, ri, ne, em, & ! 3 Dirac-Hara + HL imag part + const real & imag part ! 5 Dirac-Fock exchange with core electrons + ! ixc=0 for valence electron density + ! 6 Finite temperature RPA + ! 7 Zero-temperature RPA ! rmt r muffin tin ! rnrm r norman ! vtot(nr) total potential, including gsxc, final state @@ -93,7 +95,8 @@ subroutine rhol ( dx, x0, ri, ne, em, & IE_LOOP: do ie = 1, ne ! p2 is (complex momentum)**2 referenced to energy dep xc p2 = em(ie) - eref - if (mod(ixc,10) .lt. 5) then + ! if (mod(ixc,10) .lt. 5) then + if ((mod(ixc,10).lt.5).OR.(ixc.EQ.6).OR.(ixc.EQ.7)) then ! TTS ncycle = 0 else ncycle = 3 diff --git a/src/LDOS/rhol_h.f90 b/src/LDOS/rhol_h.f90 index d67b8d2..a04f7ed 100755 --- a/src/LDOS/rhol_h.f90 +++ b/src/LDOS/rhol_h.f90 @@ -107,7 +107,8 @@ subroutine rhol_h (iph, dx, x0, ri, ne, em, & do ie = 1, ne ! p2 is (complex momentum)**2 referenced to energy dep xc p2 = em(ie) - eref - if (mod(ixc,10) .lt. 5) then + ! if (mod(ixc,10) .lt. 5) then + if ((mod(ixc,10).lt.5).OR.(ixc.EQ.6).OR.(ixc.EQ.7)) then ! TTS ncycle = 0 else ncycle = 3 diff --git a/src/LDOS/rhol_h_step1.f90 b/src/LDOS/rhol_h_step1.f90 index 63bf158..efdff74 100644 --- a/src/LDOS/rhol_h_step1.f90 +++ b/src/LDOS/rhol_h_step1.f90 @@ -87,7 +87,8 @@ subroutine rhol_h_step1 (iph, dx, x0, ri, ne, em, ixc, rmt, rnrm, do ie = 1, ne ! p2 is (complex momentum)**2 referenced to energy dep xc p2 = em(ie) - eref - if (mod(ixc,10) .lt. 5) then + ! if (mod(ixc,10) .lt. 5) then + if ((mod(ixc,10).lt.5).OR.(ixc.EQ.6).OR.(ixc.EQ.7)) then ! TTS ncycle = 0 else ncycle = 3 diff --git a/src/LDOS/rhol_h_step2.f90 b/src/LDOS/rhol_h_step2.f90 index 88ad394..c480e0a 100644 --- a/src/LDOS/rhol_h_step2.f90 +++ b/src/LDOS/rhol_h_step2.f90 @@ -89,7 +89,8 @@ subroutine rhol_h_step2 (iph, dx, x0, ri, ne, em, ixc, rmt, rnrm, & do ie = 1, ne ! p2 is (complex momentum)**2 referenced to energy dep xc p2 = em(ie) - eref - if (mod(ixc,10) .lt. 5) then + ! if (mod(ixc,10) .lt. 5) then + if ((mod(ixc,10).lt.5).OR.(ixc.EQ.6).OR.(ixc.EQ.7)) then ! TTS ncycle = 0 else ncycle = 3 diff --git a/src/PAR/parallel.src b/src/PAR/parallel.src index 824b06c..48330c8 100644 --- a/src/PAR/parallel.src +++ b/src/PAR/parallel.src @@ -288,6 +288,18 @@ return end + subroutine par_bcast_double_scalar(buf,count,source) +! ************************************************** +! Call mpi_bcast for double arrays +! ************************************************** + integer count,source + double precision buf + include 'mpif.h' + call MPI_BCAST(buf,count,MPI_DOUBLE_PRECISION,source, & + & MPI_COMM_WORLD,ierrorflag) + return + end + subroutine MPE_DECOMP1D( n, num_procs, myid, s, e ) ! ****************************************************** ! A routine for producing a decomposition of a 1-d diff --git a/src/PAR/sequential.src b/src/PAR/sequential.src index dc28f92..63f31a2 100644 --- a/src/PAR/sequential.src +++ b/src/PAR/sequential.src @@ -221,6 +221,15 @@ return end + subroutine par_bcast_double_scalar(buf,count,source) +! ************************************************** +! Call mpi_bcast for double arrays +! ************************************************** + integer count, source + double precision buf + return + end + subroutine MPE_DECOMP1D( n, num_procs, myid, s, e ) ! ****************************************************** ! A routine for producing a decomposition of a 1-d diff --git a/src/POT/corval.f90 b/src/POT/corval.f90 index d12538b..b3fa708 100644 --- a/src/POT/corval.f90 +++ b/src/POT/corval.f90 @@ -138,7 +138,8 @@ subroutine corval ( ecv, xnvmu, eorb, norb, xnval, kappa, rgrd, & do 500 iph = 0, nph call fixvar (rmt(iph),edens(1,iph),vtot(1,iph),dmag0, vint, rhoint, dx, rgrd, jumprm, vjump, ri, vtotph, dum, dmagx) - if (mod(ixc,10) .ge.5) then + ! if (mod(ixc,10) .ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.15)) THEN ! Replace mod(ixc,10).ge.5 if (jumprm .gt. 0) jumprm = 2 call fixvar (rmt(iph), edenvl(1,iph), vvalgs(1,iph), dmag0, & vint, rhoint, dx, rgrd, jumprm, vjump, ri, vvalph, dum, dmagx) @@ -151,7 +152,8 @@ subroutine corval ( ecv, xnvmu, eorb, norb, xnval, kappa, rgrd, & eref = vtotph(jri1) do 40 i = 1, jri1 40 vtotph(i) = vtotph(i) - eref - if (ixc.ge.5) then + ! if (ixc.ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.10).OR.(ixc.EQ.13).OR.(ixc.EQ.15)) THEN ! Replace ixc.ge.5 do 50 i = 1, jri1 50 vvalph(i) = vvalph(i) - eref else diff --git a/src/POT/istprm.f90 b/src/POT/istprm.f90 index 5ff5050..2cd0d8c 100644 --- a/src/POT/istprm.f90 +++ b/src/POT/istprm.f90 @@ -315,6 +315,10 @@ subroutine istprm ( nph, nat, iphat, rat, iatph, xnatph, & call vbh( rs, xmag, vvbh ) case(12) call pz_vxc( rs, vvbh ) + case default + write(slog, 911) 'Invalid choice of xc potential ', iscfxc + 911 format (a, i3) + call wlog(slog) end select ! else @@ -332,14 +336,16 @@ subroutine istprm ( nph, nat, iphat, rat, iatph, xnatph, & ! end if vtot(i,iph) = vclap(i,iph) + vvbh - if (mod(ixc,10).eq.5) then + ! if (mod(ixc,10).eq.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.15)) THEN rsval = 10.0 if (edenvl(i,iph) .gt. 0.00001) rsval = (edenvl(i,iph)/3)**(-third) if (rsval.gt.10.0) rsval = 10.0 xmagvl = 1.0 + idmag * dmag(i,iph) * edens(i,iph) / edenvl(i,iph) call vbh(rsval,xmagvl,vvbhvl) vvalgs(i,iph) = vclap(i,iph) + vvbhvl - elseif (mod(ixc,10) .ge. 6) then + ! elseif (mod(ixc,10) .ge. 6) then + ELSEIF (ixc.EQ.9) THEN if (edens(i,iph).le.edenvl(i,iph)) then rscore =101.0 else @@ -397,7 +403,8 @@ subroutine istprm ( nph, nat, iphat, rat, iatph, xnatph, & call ovp2mt(nph, edens, 0, qtotel, ri, xnatph, lnear, inrm, imt, rnrm, rmt, cmovp, ipiv, rhoint,inters) rhoint = 4*pi * rhoint / volint - if (ixc.ge.5) then + ! if (ixc.ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.10).OR.(ixc.EQ.13).OR.(ixc.EQ.15)) THEN ! Replace ixc.ge.5 ! find valence potential inside mt sphere (vintvl -dummy) call ovp2mt(nph, vvalgs, 1, qtotel, ri, xnatph, lnear, inrm, imt, rnrm, rmt, cmovp, ipiv, vintvl,inters) endif diff --git a/src/POT/m_sommerfeld_scf.f90 b/src/POT/m_sommerfeld_scf.f90 index f63cbe8..250b10b 100644 --- a/src/POT/m_sommerfeld_scf.f90 +++ b/src/POT/m_sommerfeld_scf.f90 @@ -51,7 +51,7 @@ MODULE m_sommerfeld_scf use constants use atoms_inp, only: nat, nph, iatph, iphat, rat use potential_inp, only: ixc, xnatph, xion, iunf, iz, ihole, lmaxsc, & - nohole, nscmt, icoul, ca1, rfms1, lfms1, jumprm, scf_temperature, xntol, nmu + nohole, nscmt, icoul, ca1, rfms1, lfms1, jumprm, scf_temperature, xntol, nmu, emaxscf, negrid implicit none private @@ -73,7 +73,8 @@ MODULE m_sommerfeld_scf DOUBLE PRECISION, ALLOCATABLE :: history_xntot(:), history_xmu(:) COMPLEX*16, ALLOCATABLE :: rhoe1(:,:,:), rhore1(:,:,:) ! points for d/dx DOUBLE PRECISION:: lower_bound, upper_bound - COMPLEX*16 :: dstep ! for d/dx + ! COMPLEX*16 :: dstep ! for d/dx + REAL*8 :: dstep CONTAINS subroutine sommerfeld_scf_init(ecv0, mu, iscmt) @@ -94,7 +95,8 @@ subroutine sommerfeld_scf_init(ecv0, mu, iscmt) end do ! ne = 80 - ne = 220 + ! ne = 220 + ne = negrid call sommerfeld_scf_energies_init() end subroutine @@ -120,6 +122,7 @@ subroutine sommerfeld_scf_energies_init() eimmax = 0.15d0 ! 4.05 eV same as grids.f90 ! eimmax = 0.30d0 emax = mu0 + 0.5 + emax = emaxscf e1 = 2*pi*kT np = 1 @@ -277,7 +280,7 @@ subroutine sommerfeld_scf_main(iscmt, ecv, vclap, edens, edenvl, vtot, vvalgs, & dstep = 0.001d0/hart ! step size of 0.02 V ! dstep = 1e-6 endif - call par_bcast_double(dstep, 1, 0) + call par_bcast_double_scalar(dstep, 1, 0) call par_barrier ! Compute points needed for derivatives @@ -343,7 +346,7 @@ subroutine sommerfeld_scf_main(iscmt, ecv, vclap, edens, edenvl, vtot, vvalgs, & call par_barrier call par_bcast_int(converged, 1, 0) call par_bcast_int(imu, 1, 0) - call par_bcast_double(xmunew, 1, 0) + call par_bcast_double_scalar(xmunew, 1, 0) if (converged.eq.1) then ok = .true. @@ -360,7 +363,7 @@ subroutine sommerfeld_scf_main(iscmt, ecv, vclap, edens, edenvl, vtot, vvalgs, & ! the rest of this routine is quick and is performed by all nodes - call par_bcast_double(xntot, 1, 0) + call par_bcast_double_scalar(xntot, 1, 0) call par_bcast_double(xnmues, (lx+1)*(nphx+1), 0) call par_bcast_double(rhoval, 251*(nphx+1), 0) diff --git a/src/POT/m_thermal_scf.f90 b/src/POT/m_thermal_scf.f90 index b89f26a..b7fdc00 100644 --- a/src/POT/m_thermal_scf.f90 +++ b/src/POT/m_thermal_scf.f90 @@ -68,7 +68,6 @@ module m_thermal_scf DOUBLE PRECISION, ALLOCATABLE :: history_xntot(:), history_xmu(:) real*8, parameter :: eimmax = 0.15d0 REAL*8 :: lower_bound, upper_bound - CONTAINS subroutine thscf_init(ecv0, mu, iscmt) @@ -326,7 +325,7 @@ subroutine thscf_main(iscmt, ecv, vclap, edens, edenvl, vtot, vvalgs, & call par_bcast_int(converged, 1, 0) call par_bcast_int(imu, 1, 0) - call par_bcast_double(xmunew, 1, 0) + call par_bcast_double_scalar(xmunew, 1, 0) if (converged.eq.1) then ok = .true. @@ -342,7 +341,7 @@ subroutine thscf_main(iscmt, ecv, vclap, edens, edenvl, vtot, vvalgs, & end do ! the rest of this routine is quick and is performed by all nodes - call par_bcast_double(xntot, 1, 0) + call par_bcast_double_scalar(xntot, 1, 0) call par_bcast_double(xnmues, (lx+1)*(nphx+1), 0) call par_bcast_double(rhoval, 251*(nphx+1), 0) @@ -360,28 +359,31 @@ subroutine thscf_main(iscmt, ecv, vclap, edens, edenvl, vtot, vvalgs, & if (diff.gt.13.1 .or. (il.eq.2 .and. diff.gt. 9.1) .or. & & (il.eq.1 .and. diff.gt.5.1) .or. (il.eq.0 .and. diff.gt.1.95)) then ! the difference in number of atoms is too large + ! print*, "Difference = ", diff, " computed = ", xnmues(il,ip), " ref = ", xnvmu(il,ip) call wlog (' Found bad counts.') write (slog,311) xnvmu(il,ip) 311 format(' Occupation number in getorb is ', f9.3) call wlog(slog) call wlog (' Will repeat this iteration. ') if (iscmt.gt.1) ok = .false. + ! call par_stop endif end do end do ! III. if ok, then update output values using Broyden algorithm. otherwise, scf loop will restart if (ok) then + xmu = xmunew ! find rhoval via Broyden algorithm call broydn( iscmt, ca1, nph, xnvmu, nr05 , xnatph, rnrm, qnrm, edenvl, rhoval, dq) ! xnatph,qrnm,edenvl,rhoval => rhoval,dq ! calculate new vclap - overlap coulomb potential call coulom (icoul, nph, nr05 , rhoval, edenvl, edens, nat, rat, iatph, iphat, rnrm, dq, iz, vclap) ! everything => vclap - + ! update array edens do iph=0,nph - do ir=1,nr05 (iph) + do ir=1, nr05 (iph) edens(ir,iph)=edens(ir,iph)-edenvl(ir,iph)+rhoval(ir,iph) end do do ir=nr05 (iph)+1,251 @@ -408,7 +410,9 @@ SUBROUTINE bracketing_method(current_step, xnferm, xmunew, xmuprev) integer :: iOrder(current_step), i_left, i_right, imu, ip double precision :: sorted_xmu(current_step), sorted_xntot(current_step) character(512) :: slog - + write(slog,*) "Calling Bracketing " + call wlog(slog) + ! call par_stop sorted_xmu(:) = -100.d0 sorted_xntot(:) = -1.d0 CALL qsortd(iOrder, current_step, history_xmu(1:current_step)) @@ -640,7 +644,7 @@ SUBROUTINE integrate_rhos_lt(nr05,xmu,xntot,xnmues,rhoval,print_flag) xnmues = 0 rhoval = 0 dx = window_size_terp*temp - dxmu = (2*dx)/nxmu + dxmu = (2.d0*dx)/nxmu xmu_m3 = xmu - dx + coni*DIMAG(energies(11)) xmu_p3 = xmu + dx + coni*DIMAG(energies(11)) @@ -719,30 +723,6 @@ SUBROUTINE integrate_rhos_lt(nr05,xmu,xntot,xnmues,rhoval,print_flag) ENDDO ENDIF - ! OPEN(unit=99, file='RHOE_TH.dat') - ! OPEN(unit=1233, file="rhoe_th.dat") - ! OPEN(unit=899, file='FERMI_TH.dat') - ! DO ie = 1, ne_terp - ! f = fermiF(energies_terp(ie), temp, xmu) - ! IF (ie > ne_terp) f = -2*pi*coni - ! WRITE(899,'(20E20.10E3)') DBLE(energies_terp(ie)*hart), DBLE(f) - ! WRITE(99,'(20E20.10E3)') DBLE(energies_terp(ie)*hart), DIMAG(energies_terp(ie)), DBLE(rhoe_terp(0,0,ie)*f), DIMAG(rhoe_terp(0,0,ie)*f), & - ! & DBLE(rhoe_terp(1,0,ie)*f), DIMAG(rhoe_terp(1,0,ie)*f), DBLE(rhoe_terp(2,0,ie)*f), DIMAG(rhoe_terp(2,0,ie)*f) - ! - ! ENDDO - ! DO ie=1,ne - ! f = fermiF(energies(ie), temp, xmu) - ! f = 1.d0 - ! WRITE(1233,'(20E20.10E3)') DBLE(energies(ie)*hart), DIMAG(energies(ie)), DBLE(rhoe(0,0,ie)*f), DIMAG(rhoe(0,0,ie)*f), & - ! & DBLE(rhoe(1,0,ie)*f), DIMAG(rhoe(1,0,ie)*f), DBLE(rhoe(2,0,ie)*f), DIMAG(rhoe(2,0,ie)*f) - ! ENDDO - ! WRITE(99,*) - ! WRITE(899,*) - ! CLOSE(99) - ! CLOSE(899) - ! CLOSE(1233) - - ! I. Integrate density over energy, fermi distribution ! The first calculated point has finite imaginary part. So, include contribution from ! leg between real axis and this point by approximating the integrand as constant @@ -756,7 +736,13 @@ SUBROUTINE integrate_rhos_lt(nr05,xmu,xntot,xnmues,rhoval,print_flag) f = fermiF(ee, temp, xmu) fp = fermiF(ep, temp, xmu) - + IF (ie.GT.10) THEN + ! No imaginary part for contour between two poles + f = fermiF(DBLE(ee)+coni*0.d0, temp, xmu) + fp = fermiF(DBLE(ep)+coni*0.d0, temp, xmu) + de = DBLE(de) + ENDIF + DO iph = 0, nph DO il = 0, lx if (iunf.EQ.0 .AND. il.GT.2) CYCLE @@ -964,13 +950,9 @@ COMPLEX*16 FUNCTION fermiF(x,t,mu) IMPLICIT none COMPLEX*16, INTENT(IN) :: x REAL(8), INTENT(IN) :: t, mu - REAL(8) :: b - COMPLEX*16 :: a IF (t.GT.0.d0) THEN - a = (x-mu)/t - b = DIMAG((x-mu)/t) - if (DBLE(x-mu)/t.le.5e2) then - fermiF = 1/(1+EXP((x-mu)/t)) + if (DBLE(x-mu)/t.LE.MAXEXPONENT(DBLE(x))) THEN + fermiF = 1.d0/(1.d0+EXP((x-mu)/t)) else fermiF = 0.d0 endif diff --git a/src/POT/potsub.f90 b/src/POT/potsub.f90 index b35cbf0..b70df0f 100644 --- a/src/POT/potsub.f90 +++ b/src/POT/potsub.f90 @@ -22,7 +22,7 @@ subroutine pot !KJ put everything in modules 7-09 ! for unique potentials specifed by atoms and overlap cards. implicit none !double precision (a-h, o-z) -!Changed the dimensions to 40 to account for superheavy elements. Pavlo Baranov 07/2016 + !Changed the dimensions to 40 to account for superheavy elements. Pavlo Baranov 07/2016 ! Notes: ! nat number of atoms in problem ! nph number of unique potentials @@ -158,10 +158,10 @@ subroutine pot !KJ put everything in modules 7-09 enddo CALL ReadAtomicPots(nph, iz(0:nph), ihole, rho, dmag(:,0:nph+1), rhoval, vcoul, dgc0, & - & dpc0, dgc(:,:,0:nph+1), dpc(:,:,0:nph+1), adgc(:,:,0:nph+1), adpc(:,:,0:nph+1), & - & erelax, emu, xnvmu, xnval(:,0:nph+1), norb, eorb, drho, dvcoul, iphat, & - & ratTmp, iatph(0:nph), novr(0:nph), iphovr, nnovr, rovr, nat, edens, & - & edenvl, vclap, rnrm(0:nph), kappa(:,0:nph+1), iorb(:,0:nph+1), s02) + & dpc0, dgc(:,:,0:nph+1), dpc(:,:,0:nph+1), adgc(:,:,0:nph+1), adpc(:,:,0:nph+1), & + & erelax, emu, xnvmu, xnval(:,0:nph+1), norb, eorb, drho, dvcoul, iphat, & + & ratTmp, iatph(0:nph), novr(0:nph), iphovr, nnovr, rovr, nat, edens, & + & edenvl, vclap, rnrm(0:nph), kappa(:,0:nph+1), iorb(:,0:nph+1), s02) edensTrack(:,0)=edens(:,0) !atomic overlap for ipot=0 goes into iteration 0 of density tracker @@ -344,7 +344,8 @@ subroutine pot !KJ put everything in modules 7-09 ! need to update vxcval in case the core-valence separation was ! made in subroutine corval. Need vxcval only for nonlocal exchange. - if (mod(ixc,10).ge.5) then + ! if (mod(ixc,10).ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.15)) THEN ! Replace mod(ixc,10).ge.5 call istprm (nph, nat, iphat, rat, iatph, xnatph, & & novr, iphovr, nnovr, rovr, folp, folpx, iafolp, & & edens, edenvl, idmag, & @@ -353,7 +354,7 @@ subroutine pot !KJ put everything in modules 7-09 & rnrmav, qtotel, inters, totvol) xmunew = xmu endif - if(master) then + if (master) then write(slog,131) ecv*hart call wlog(slog) write(slog,130) xmu*hart @@ -364,6 +365,7 @@ subroutine pot !KJ put everything in modules 7-09 131 format('Core-valence separation energy: ecv= ',f9.3,' eV') 132 format(' ca1= ',f9.3) endif + dosTrack=0.d0 ! do first nmix iterations with mixing scheme. Need for f-elements. @@ -565,104 +567,102 @@ subroutine pot !KJ put everything in modules 7-09 if (nohole.gt.0) then ! testing new final state potential - do 220 j = 1,251 - 220 edens(j,0) = edens(j,0) - drho(j) - - ! notice that vclap is actually for the next iteration - ! in SCMT loop, thus vclap may be wrong if self-consistency has not been reached - do 230 j = 1,251 - 230 vclap(j,0) = vclap(j,0) - dvcoul(j) - - call istprm (nph, nat, iphat, rat, iatph, xnatph, & - & novr, iphovr, nnovr, rovr, folp, folpx, iafolp, & - & edens, edenvl, idmag, & - & dmag, vclap, vtot, vvalgs, imt, inrm, rmt, rnrm, & - & ixc, rhoint,vint, rs, xf, xmu, xmunew, & - & rnrmav, qtotel, inters, totvol) - endif - - ! correct the excitation energy - ! emu = emu -vclap(1,0) + vcoul(1,0) done also above - ! emu = emu+xmu should be done in principle but leads - ! to worse estimate of edge position. fix later. ala - - if (ipr1 .ge. 2 .and. master) call wpot (nph, edens, imt, inrm, rho, vclap, vcoul, vtot, ntitle, title) - - - ! Debug: Fer - ! Changing eorb by hand for now - ! eorb(1,0) = -14.0050405749 - ! eorb(1,0) = -14.1374268761 - - ! Debug: Fer - ! write(6,fmt='(a,f16.10)') ' eorb in pot: ', eorb(1,0) - - ! Debug: Fer - ! Write out some of the stuff that goes into pot.bin - ! print *, 'xmu: ', xmu - ! print *, 'vint: ', vint - ! print *, 'eorb: ', eorb(1,0) - ! print *, 'emu: ', emu - ! print *, 'erelax: ', erelax - ! print *, 'iorb: ', iorb - - if ( (ChSh_Type == 1) .or. (ChSh_Type == 3) ) then - - ! Set the target center - iCenter = 0 - sh_iz = iz(iCenter) - sh_ihole = ihole - sh_rmt = rmt(iCenter) - sh_jri = imt(iCenter)+1 - sh_dx = dx - sh_ri = ri - sh_p2f = xmu-vint - sh_edge = xmu - sh_vxc(1:251) = vtot(1:251,iCenter) - vint - sh_dgcn(1:251,:) = dgc(:,:,iCenter) - sh_dpcn(1:251,:) = dpc(:,:,iCenter) - ! sh_adgc = adgc(:,:,iCenter) - ! sh_adpc = adpc(:,:,iCenter) - sh_adgc = 0.0d0 - sh_adpc = 0.0d0 - sh_eorb = eorb(:,iCenter) - sh_kappa = kappa(:,iCenter) - - call correorb(sh_iz, sh_ihole, sh_rmt, sh_jri, & - sh_dx,sh_ri, sh_p2f, sh_edge, sh_vxc, & - sh_dgcn, sh_dpcn, sh_adgc, sh_adpc, & - sh_eorb, & - sh_neg, sh_eng, sh_rhoj, sh_kappa, sh_norbp, 0) !KJ iph=0 - e_chsh = -sh_eng(1,1) - write(6,fmt='(a,f12.8)') 'e_chsh: ', e_chsh - emu = e_chsh + erelax - end if - - ! write stuff into pot.bin - if (master) call wrpot (nph, ntitle, title, rnrmav, xmu, vint, rhoint, & - & emu, s02, erelax, wp, ecv,rs,xf, qtotel, & - & imt, rmt, inrm, rnrm, folp, folpx, xnatph, & - & dgc0, dpc0, dgc, dpc, adgc, adpc, & - & edens, vclap, vtot, edenvl, vvalgs, dmag, xnval, & - & eorb(1,0), kappa(1,0), iorb, qnrm, xnmues, nhtmp, & - & ihole, inters, totvol, iafolp, xion, iunf, iz, jumprm,ipr1) !KJ added ipr1 3-2012 - - ! write misc.dat - if (ipr1 .ge. 1 .and. master) then - open (unit=1, file='misc.dat', status='unknown', iostat=ios) - call chopen (ios, 'misc.dat', 'potph') - call wthead(1, ntitle, title) - close (unit=1) - endif - - 400 call par_barrier - - OPEN(UNIT=1993, FILE='chemical.dat',STATUS='REPLACE', POSITION="APPEND") - WRITE(1993,*) scf_temperature, scf_temperature/0.000086173423, xmu*hart - CLOSE(1993) - - ! Deallocate local variables - deallocate(xnmues,xnvmu,gtr,xrhoce,xrhole,yrhoce,yrhole) - - return - end + do j = 1,251 + edens(j,0) = edens(j,0) - drho(j) + enddo + + ! notice that vclap is actually for the next iteration + ! in SCMT loop, thus vclap may be wrong if self-consistency has not been reached + do j = 1,251 + vclap(j,0) = vclap(j,0) - dvcoul(j) + enddo + + call istprm (nph, nat, iphat, rat, iatph, xnatph, & + & novr, iphovr, nnovr, rovr, folp, folpx, iafolp, & + & edens, edenvl, idmag, & + & dmag, vclap, vtot, vvalgs, imt, inrm, rmt, rnrm, & + & ixc, rhoint,vint, rs, xf, xmu, xmunew, & + & rnrmav, qtotel, inters, totvol) + endif + + ! correct the excitation energy + ! emu = emu -vclap(1,0) + vcoul(1,0) done also above + ! emu = emu+xmu should be done in principle but leads + ! to worse estimate of edge position. fix later. ala + + if (ipr1 .ge. 2 .and. master) call wpot (nph, edens, imt, inrm, rho, vclap, vcoul, vtot, ntitle, title) + + + ! Debug: Fer + ! Changing eorb by hand for now + ! eorb(1,0) = -14.0050405749 + ! eorb(1,0) = -14.1374268761 + + ! Debug: Fer + ! write(6,fmt='(a,f16.10)') ' eorb in pot: ', eorb(1,0) + + ! Debug: Fer + ! Write out some of the stuff that goes into pot.bin + ! print *, 'xmu: ', xmu + ! print *, 'vint: ', vint + ! print *, 'eorb: ', eorb(1,0) + ! print *, 'emu: ', emu + ! print *, 'erelax: ', erelax + ! print *, 'iorb: ', iorb + + if ( (ChSh_Type == 1) .or. (ChSh_Type == 3) ) then + + ! Set the target center + iCenter = 0 + sh_iz = iz(iCenter) + sh_ihole = ihole + sh_rmt = rmt(iCenter) + sh_jri = imt(iCenter)+1 + sh_dx = dx + sh_ri = ri + sh_p2f = xmu-vint + sh_edge = xmu + sh_vxc(1:251) = vtot(1:251,iCenter) - vint + sh_dgcn(1:251,:) = dgc(:,:,iCenter) + sh_dpcn(1:251,:) = dpc(:,:,iCenter) + ! sh_adgc = adgc(:,:,iCenter) + ! sh_adpc = adpc(:,:,iCenter) + sh_adgc = 0.0d0 + sh_adpc = 0.0d0 + sh_eorb = eorb(:,iCenter) + sh_kappa = kappa(:,iCenter) + + call correorb(sh_iz, sh_ihole, sh_rmt, sh_jri, & + sh_dx,sh_ri, sh_p2f, sh_edge, sh_vxc, & + sh_dgcn, sh_dpcn, sh_adgc, sh_adpc, & + sh_eorb, & + sh_neg, sh_eng, sh_rhoj, sh_kappa, sh_norbp, 0) !KJ iph=0 + e_chsh = -sh_eng(1,1) + write(6,fmt='(a,f12.8)') 'e_chsh: ', e_chsh + emu = e_chsh + erelax + end if + + ! write stuff into pot.bin + if (master) call wrpot (nph, ntitle, title, rnrmav, xmu, vint, rhoint, & + & emu, s02, erelax, wp, ecv,rs,xf, qtotel, & + & imt, rmt, inrm, rnrm, folp, folpx, xnatph, & + & dgc0, dpc0, dgc, dpc, adgc, adpc, & + & edens, vclap, vtot, edenvl, vvalgs, dmag, xnval, & + & eorb(1,0), kappa(1,0), iorb, qnrm, xnmues, nhtmp, & + & ihole, inters, totvol, iafolp, xion, iunf, iz, jumprm,ipr1) !KJ added ipr1 3-2012 + + ! write misc.dat + if (ipr1 .ge. 1 .and. master) then + open (unit=1, file='misc.dat', status='unknown', iostat=ios) + call chopen (ios, 'misc.dat', 'potph') + call wthead(1, ntitle, title) + close (unit=1) + endif + + 400 call par_barrier + + ! Deallocate local variables + deallocate(xnmues,xnvmu,gtr,xrhoce,xrhole,yrhoce,yrhole) + + return +end diff --git a/src/POT/rhofmslie.f90 b/src/POT/rhofmslie.f90 index 543158c..39f851a 100644 --- a/src/POT/rhofmslie.f90 +++ b/src/POT/rhofmslie.f90 @@ -76,7 +76,8 @@ subroutine rhofmslie(edens, edenvl, vtot, vvalgs, rmt, rnrm, rhoint, vint, & !c use spin-unpolarized case to get SCF. set dmagx to zero !c may want to replace dmag0 with dmag(1,iph) for spin-dependent extension of SCF procedure. call fixvar (rmt(iph),edens(1,iph),vtot(1,iph),dmag0, vint, rhoint, dx, rgrd, jumprm, vjump, ri, vtotph, dum, dmagx) !makes vtotph,dum,dmagx out of vtot,edens,dmag0 - if (mod(ixc,10) .ge.5) then + ! if (mod(ixc,10) .ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.15)) THEN ! Replace mod(ixc,10).ge.5 if (jumprm .gt. 0) jumprm = 2 call fixvar (rmt(iph), edenvl(1,iph), vvalgs(1,iph), dmag0, vint, rhoint, dx, rgrd , jumprm, vjump, ri, vvalph, dum, dmagx) if (jumprm .gt. 0) jumprm = 1 @@ -89,7 +90,8 @@ subroutine rhofmslie(edens, edenvl, vtot, vvalgs, rmt, rnrm, rhoint, vint, & do i = 1, jri1 vtotph(i) = vtotph(i) - eref end do - if (ixc.ge.5) then + ! if (ixc.ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.10).OR.(ixc.EQ.13).OR.(ixc.EQ.15)) THEN ! Replace ixc.ge.5 do i = 1, jri1 vvalph(i) = vvalph(i) - eref end do diff --git a/src/POT/rholie.f90 b/src/POT/rholie.f90 index b028f51..d54ebab 100644 --- a/src/POT/rholie.f90 +++ b/src/POT/rholie.f90 @@ -100,7 +100,8 @@ subroutine rholie ( ri05, nr05, dx, x0, ri, em, & ! p2 is 0.5*(complex momentum)**2 referenced to energy dep xc ! need hartree units for dfovrg p2 = em - eref - if (mod(ixc,10) .lt. 5) then +! if (mod(ixc,10) .lt. 5) then + IF ((ixc.EQ.0).OR.(ixc.EQ.1).OR.(ixc.EQ.2).OR.(ixc.EQ.3).OR.(ixc.EQ.6).OR.(ixc.EQ.7)) THEN ! Replace mod(ixc,10) .lt. 5 ncycle = 0 else ncycle = 3 diff --git a/src/POT/scmt.f90 b/src/POT/scmt.f90 index d38286e..6f9ab8f 100644 --- a/src/POT/scmt.f90 +++ b/src/POT/scmt.f90 @@ -120,7 +120,8 @@ subroutine scmt (iscmt, ecv, nph, nat, vclap, edens, edenvl, vtot, vvalgs, rmt, !c use spin-unpolarized case to get SCF. set dmagx to zero !c may want to replace dmag0 with dmag(1,iph) for spin-dependent extension of SCF procedure. call fixvar (rmt(iph),edens(1,iph),vtot(1,iph),dmag0, vint, rhoint, dx, rgrd, jumprm, vjump, ri, vtotph, dum, dmagx) !makes vtotph,dum,dmagx out of vtot,edens,dmag0 - if (mod(ixc,10) .ge.5) then + ! if (mod(ixc,10) .ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.15)) THEN ! Replace mod(ixc,10).ge.5 if (jumprm .gt. 0) jumprm = 2 call fixvar (rmt(iph), edenvl(1,iph), vvalgs(1,iph), dmag0, & vint, rhoint, dx, rgrd , jumprm, vjump, ri, vvalph, dum, dmagx) @@ -133,7 +134,8 @@ subroutine scmt (iscmt, ecv, nph, nat, vclap, edens, edenvl, vtot, vvalgs, rmt, eref = vtotph(jri1) do 40 i = 1, jri1 40 vtotph(i) = vtotph(i) - eref - if (ixc.ge.5) then + ! if (ixc.ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.10).OR.(ixc.EQ.13).OR.(ixc.EQ.15)) THEN ! Replace ixc.ge.5 do 50 i = 1, jri1 50 vvalph(i) = vvalph(i) - eref else diff --git a/src/POT/scmtmp.f90 b/src/POT/scmtmp.f90 index b7eee14..2a66882 100644 --- a/src/POT/scmtmp.f90 +++ b/src/POT/scmtmp.f90 @@ -164,7 +164,8 @@ subroutine scmtmp (npr, iscmt, ecv, nph, nat, vclap, & !c may want to replace dmag0 with dmag(1,iph) for spin-dependent !c extension of SCF procedure. call fixvar (rmt(iph),edens(1,iph),vtot(1,iph),dmag0,vint, rhoint, dx, rgrd, jumprm, vjump, ri, vtotph, dum, dmagx) - if (mod(ixc,10) .ge.5) then + ! if (mod(ixc,10) .ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.15)) THEN ! Replace mod(ixc,10).ge.5 if (jumprm .gt. 0) jumprm = 2 call fixvar (rmt(iph),edenvl(1,iph),vvalgs(1,iph),dmag0, vint, rhoint, dx, rgrd , jumprm, vjump, ri, vvalph, dum, dmagx) if (jumprm .gt. 0) jumprm = 1 @@ -176,7 +177,8 @@ subroutine scmtmp (npr, iscmt, ecv, nph, nat, vclap, & eref = vtotph(jri1) do 40 i = 1, jri1 40 vtotph(i) = vtotph(i) - eref - if (ixc.ge.5) then + ! if (ixc.ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.10).OR.(ixc.EQ.13).OR.(ixc.EQ.15)) THEN ! Replace ixc.ge.5 do 50 i = 1, jri1 50 vvalph(i) = vvalph(i) - eref else diff --git a/src/RDINP/rdinp.f90 b/src/RDINP/rdinp.f90 index 8afdbb9..744c7e6 100755 --- a/src/RDINP/rdinp.f90 +++ b/src/RDINP/rdinp.f90 @@ -324,16 +324,23 @@ program rdinp ! ixc=3 Dirac-Hara + HL imag part + const real & imag part ! ixc=5 partially nonlocal: Dirac-Fock for core + HL for ! valence electrons, + const real & imag part -! ixc=10 same as ixc=0 with broadened plasmon HL selfenergy -! ixc=13 same as ixc=3 with broadened plasmon HL selfenergy -! ixc=15 same as ixc=5 with broadened plasmon HL selfenergy +! ixc=6 Finite temperature RPA GW +! ixc=9 Undocumented options, it does the following (not sure whats the purpose) +! For unnocupied states: +! Use HL for the valence electrons: +! Use HL - DH (correlation only) for core electrons. +! For occupied states: +! Use VBH(rs) - DH(rs_core, k_fermi(rs_core) +! ixc=10 same as ixc=0 with broadened plasmon HL self-energy +! ixc=13 same as ixc=3 with broadened plasmon HL self-energy +! ixc=15 same as ixc=5 with broadened plasmon HL self-energy ! vr0 is const imag part of potential ! vi0 is const imag part of potential ! Default is HL. (ixc=0, vr0=0, vi0=0, ixc0 = 2) vr0=0.0 vi0=0.0 read(words(2),20,err=900) ixc - if (ixc.lt.0) ixc=0 !default + if (ixc.lt.0) ixc=0 !default ! if (nwords.ge.3) (read(words(3),30,err=900) vr0 read(words(3),30,err=900) vr0 ! if (nwords.ge.4) read(words(4),30,err=900) vi0 @@ -1069,6 +1076,7 @@ program rdinp elseif (itok .eq. 76) then ! Added by Fer read(words(2),*,err=900) ChSh_Type + if (nwords.gt.2) read(words(3),*,err=900) custom_chsh ! TTS 05/22 custom core level shift elseif (itok .eq. 77) then ! DIMS card ! First is nclusx, second is lx @@ -1174,7 +1182,7 @@ program rdinp IF(iph.gt.nphxhardlimit) then call wlog("iph > nphxhardlimit in feff.inp") call wlog(TRIM(ADJUSTL(line))) - call par_stop + call par_stop('') END IF read(words(3),30) NumDens(iph) elseif (itok .eq. 85) then diff --git a/src/XSPH/phase.f90 b/src/XSPH/phase.f90 index 9ac56c9..beef9af 100644 --- a/src/XSPH/phase.f90 +++ b/src/XSPH/phase.f90 @@ -59,7 +59,7 @@ subroutine phase (iph, dx, x0, ri, ne, ne1, ne3, em, & complex*16 eref(nex) complex*16 ph(nex,-ltot:ltot) integer ispin - + character(len=50) :: fname ! work space for xcpot dimension vxcrmu(nrptx), vxcimu(nrptx), gsrel(nrptx) dimension vvxcrm(nrptx), vvxcim(nrptx) @@ -162,8 +162,9 @@ subroutine phase (iph, dx, x0, ri, ne, ne1, ne3, em, & ColumnLabels(:) = ' ' ColumnLabels(1) = 'Rs' ColumnLabels(1) = 'wp (eV)' -!KJ write(*,*) 'writing mpse.dat in PHASE' - CALL WriteData('mpse.dat',Double1 = (3 / (4*pi*edens(jri+1))) ** third, & +!KJ write(*,*) 'writing mpse.dat in PHASE' + fname='mpse' + CALL WriteData(TRIM(fname)//'.dat',Double1 = (3 / (4*pi*edens(jri+1))) ** third, & & Double2 = SQRT(3.d0/((3 / (4*pi*edens(jri+1))) ** third)**3)*hart, & & ColumnLabels = ColumnLabels, WriteDataInHeader = .TRUE., & & Headers = (/ 'This file contains information about the' //& @@ -179,7 +180,7 @@ subroutine phase (iph, dx, x0, ri, ne, ne1, ne3, em, & & em(ie), xmu, & & vtot, vvalgs, edens, dmag, edenvl, & & eref(ie), v, vval, iPl, WpCorr, Gamma, AmpFac, EGap, & - & vxcrmu, vxcimu, gsrel, vvxcrm, vvxcim) + & vxcrmu, vxcimu, gsrel, vvxcrm, vvxcim, fname) if (dble(em(ie)).lt.-10.d0 .or. dble(em(ie)) .gt.3.d2) goto 220 ! p2 is (complex momentum)**2 referenced to energy dep xc @@ -197,7 +198,8 @@ subroutine phase (iph, dx, x0, ri, ne, ne1, ne3, em, & call besjn (xkmt, jl, nl) call besjn (xkmtEC, jlEC, nlEC) - if (mod(ixc,10) .lt. 5) then + ! if (mod(ixc,10) .lt. 5) then + IF ((ixc.EQ.0).OR.(ixc.EQ.1).OR.(ixc.EQ.2).OR.(ixc.EQ.3).OR.(ixc.EQ.6).OR.(ixc.EQ.7)) THEN ! Replace mod(ixc,10) .lt. 5 ncycle = 0 else ncycle = 3 @@ -310,7 +312,7 @@ subroutine phase (iph, dx, x0, ri, ne, ne1, ne3, em, & 220 continue ! Josh - Close mpse.dat, rl.dat IF((iph.EQ.0).AND.PrintRl) CALL CloseFl('rl.dat') - CALL CloseFl('mpse.dat') + CALL CloseFl(TRIM(fname)//'.dat') ! Josh END do ie = max(1,ne12+1), ne !KJ Feb 14 - added "max(1" to respect array bounds diff --git a/src/XSPH/phase_h.f90 b/src/XSPH/phase_h.f90 index 1716157..1143e53 100644 --- a/src/XSPH/phase_h.f90 +++ b/src/XSPH/phase_h.f90 @@ -198,7 +198,8 @@ subroutine phase_h (iph, dx, x0, ri, ne, ne1, ne3, em, & call besjn (xkmt, jl, nl) call besjn (xkmtEC, jlEC, nlEC) - if (mod(ixc,10) .lt. 5) then + ! if (mod(ixc,10) .lt. 5) then + IF ((ixc.EQ.0).OR.(ixc.EQ.1).OR.(ixc.EQ.2).OR.(ixc.EQ.3).OR.(ixc.EQ.6).OR.(ixc.EQ.7)) THEN ! Replace mod(ixc,10) .lt. 5 ncycle = 0 else ncycle = 3 diff --git a/src/XSPH/rholat.f90 b/src/XSPH/rholat.f90 index 1768799..14fcdb9 100644 --- a/src/XSPH/rholat.f90 +++ b/src/XSPH/rholat.f90 @@ -144,7 +144,8 @@ subroutine rholat ( icount, dx, x0, ri, em, & ! p2 is 0.5*(complex momentum)**2 referenced to energy dep xc ! need hartree units for dfovrg p2 = em - eref - if (mod(ixc,10) .lt. 5) then + ! if (mod(ixc,10) .lt. 5) then + IF ((ixc.EQ.0).OR.(ixc.EQ.1).OR.(ixc.EQ.2).OR.(ixc.EQ.3).OR.(ixc.EQ.6).OR.(ixc.EQ.7)) THEN ! Replace mod(ixc,10) .lt. 5 ncycle = 0 else ncycle = 3 diff --git a/src/XSPH/rholsz.f90 b/src/XSPH/rholsz.f90 index 4798909..5b2e9d7 100644 --- a/src/XSPH/rholsz.f90 +++ b/src/XSPH/rholsz.f90 @@ -100,7 +100,8 @@ subroutine rholsz ( dx, x0, ri, em, & ! p2 is 0.5*(complex momentum)**2 referenced to energy dep xc ! need hartree units for dfovrg p2 = em - eref - if (mod(ixc,10) .lt. 5) then + ! if (mod(ixc,10) .lt. 5) then + IF ((ixc.EQ.0).OR.(ixc.EQ.1).OR.(ixc.EQ.2).OR.(ixc.EQ.3).OR.(ixc.EQ.6).OR.(ixc.EQ.7)) THEN ! Replace mod(ixc,10) .lt. 5 ncycle = 0 else ncycle = 3 diff --git a/src/XSPH/szlz.f90 b/src/XSPH/szlz.f90 index d26a969..26ac521 100644 --- a/src/XSPH/szlz.f90 +++ b/src/XSPH/szlz.f90 @@ -111,7 +111,8 @@ subroutine szlz (ispin, ecv, nph, nat, rgrd, nohole, rfms2, lfms2,& call fixvar (rmt(iph),edens(1,iph),vtot(1,iph),dmag(1,iph), & & vint, rhoint, dx, rgrd, jumprm, & & vjump, ri, vtotph, dum, dmagx) - if (mod(ixc,10) .ge.5) then + ! if (mod(ixc,10) .ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.15)) THEN ! Replace mod(ixc,10).ge.5 if (jumprm .gt. 0) jumprm = 2 call fixvar (rmt(iph), edenvl(1,iph), vvalgs(1,iph), & & dmag(1,iph), vint, rhoint, dx, rgrd , jumprm, & @@ -125,7 +126,8 @@ subroutine szlz (ispin, ecv, nph, nat, rgrd, nohole, rfms2, lfms2,& eref = vtotph(jri1) do 40 i = 1, jri1 40 vtotph(i) = vtotph(i) - eref - if (ixc.ge.5) then + ! if (ixc.ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.10).OR.(ixc.EQ.13).OR.(ixc.EQ.15)) THEN ! Replace ixc.ge.5 do 50 i = 1, jri1 50 vvalph(i) = vvalph(i) - eref else diff --git a/src/XSPH/xsect.f90 b/src/XSPH/xsect.f90 index 7d5adeb..623cf10 100644 --- a/src/XSPH/xsect.f90 +++ b/src/XSPH/xsect.f90 @@ -7,754 +7,755 @@ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Josh - argument iPl has been added to arguments of xsect ! Josh - added NPoles and Eps0 3/8/2010 - subroutine xsect (ipr2, dx, x0, ri, ne, ne1, ik0, em, edge, & - ihole, emu, corr, dgc0, dpc0, jnew, & - ixc, lreal, rmt, rnrm, xmu, & - vi0, iPl, NPoles, Eps0, EGap, gamach, & - vtot, vvalgs, edens, dmag, edenvl, & - dgcn, dpcn, adgc, adpc, xsec, xsnorm, rkk, & - iz, xion, iunf, xnval, & - izstd, ifxc, eorb, kappa, iorb, l2lp, & - ipol, ispin, le2, angks, ptz, iph) !KJ - - USE IOMOD - use constants - use DimsMod, only: nex, nrptx, lx, MxPole, nspx=>nspu - USE SelfEnergyMod -! right now the same self-energy is used for calculation of the central atom part (xsec) and dipole m.e. for -! scattering (rkk). You may want to run xsect separately for xsec and for rkk, if you want to use different self-energy -! for central and scattering parts. ala. fix later - - !implicit double precision (a-h, o-z) - implicit none -!Changed the dimensions to 40 to account for superheavy elements. Pavlo Baranov 07/2016 - integer, intent(in) :: ihole, ipr2, ne, ne1, ik0, l2lp, le2, ispin, ipol, iPl, & - ixc, lreal, iunf, izstd, iz, ifxc - integer, intent(inout) :: NPoles, iph !KJ I hope that is as intended?? - real*8, intent(in) :: dx, x0, rmt, rnrm, xmu, vi0, gamach, corr, xion, emu -! INPUT -! dx, x0, ri(nr) -! Loucks r-grid, ri=exp((i-1)*dx-x0) -! ne, em(ne) number of energy points, real energy grid -! edge chemical potential (energy for k=0) -! ihole hole code -! emu position of chemical potential in absorption specrum -! dgc0(nr) dirac upper component, ground state hole orbital -! dpc0(nr) dirac lower component, ground state hole orbital -! ixc 0 Hedin-Lunqist + const real & imag part -! 1 Dirac-Hara + const real & imag part -! 2 ground state + const real & imag part -! 3 Dirac-Hara + HL imag part + const real & imag part -! 5 Dirac-Fock exchange with core electrons + -! ixc=0 for valence electron density -! lreal logical, true for real phase shifts only -! rmt r muffin tin -! xmu fermi level -! vi0 const imag part to add to complex potential -! gamach core hole lifetime -! vtot(nr) total potential, including gsxc, final state -! edens(nr) density, hole orbital, final state -! dmag(251) density magnetization -! edenvl valence charge density -! dgcn(dpcn) large (small) dirac components for central atom -! adgc(adpc) their development coefficients -! -! OUTPUT -! xsec(ne) atomic absorption cross section to multiply \chi -! (atomic background for XMCD) -! xsnorm(ne) atomic absorption cross section (norm for XMCD) -! rkk(ne, 8) normalized reduced matrix elements for construction -! of termination matrix in genfmt. - - complex*16, intent(in) :: ptz(-1:1, -1:1) - complex*16, intent(in) :: em(nex) - real*8, intent(in) :: ri(nrptx), vtot(nrptx), edens(nrptx),dmag(nrptx) - real*8, intent(in) :: dgc0(nrptx), dpc0(nrptx), vvalgs(nrptx), edenvl(nrptx) - real*8, intent(in) :: dgcn(nrptx,41), dpcn(nrptx,41), eorb(41) - real*8, intent(in) :: adgc(10,41), adpc(10,41), xnval(41) - integer, intent(in) :: iorb(-5:4),kappa(41) - complex*16, intent(out) :: rkk(nex, 8), xsec(nex) - real*8, intent(out) :: xsnorm(nex) - - integer kiind(8), lind(8) !KJ renamed kind to kiind to avoid conflict with reserved name and its use in modules - real*8 xp(nrptx), xq(nrptx) - -! work space for xcpot - real*8 vxcrmu(nrptx), vxcimu(nrptx), gsrel(nrptx) - real*8 vvxcrm(nrptx), vvxcim(nrptx) - -! work space for fovrg - complex*16 p(nrptx), q(nrptx), pn(nrptx), qn(nrptx), fscf(nrptx) - complex*16 pp(nrptx), qp(nrptx), pnp(nrptx), qnp(nrptx) -! storage for calculation of cross term (SPIN 1 only) - complex*16 xrcold(nrptx) , xncold(nrptx), yvec(nrptx,1) - - complex*16 p2, ck, xkmt, xkmtp - complex*16 pu, qu, dum1, factor - complex*16 xfnorm, xirf, xirf1 - complex*16 temp, aa, bb, cc, rkk1, rkk0, phold - complex*16 phx(8), ph0 - complex*16 eref, xm1, xm2, xm3, xm4 - - complex*16 jl,jlp1,nl,nlp1 - complex*16 v(nrptx), vval(nrptx) - complex*16 xrc(nrptx), xnc(nrptx) - character*512 slog - logical ltrace - complex*16 xrhoce(nex), xrhopr(nex), chia(nex), cchi(nex) - real*8 omega1(nex), bf(0:2, nrptx) - - real*8 pat(nrptx),qat(nrptx) - complex*16 intr(nrptx),var(nrptx) -! to pass energy levels and projected DOS - integer neg(41) - real*8 eng(nex, 41), rhoj(nex,41) -! Josh - Added iPl switch for PLASMON card -! - and WpCorr = Wi/Wp, Gamma, AmpFac -! - to describe Im[eps^-1] - integer ipole, NData - real*8, allocatable :: Energy(:), Loss(:) - real*8 WpCorr(MxPole), Gamma(MxPole), AmpFac(MxPole), Eps0, EGap - character(LEN=10) ColumnLabels(20) -! Josh END - integer imt, jri, jri1, inrm, jnrm, i, index, ifirst, ie, ncycle, ilast, ll, & - kinit, linit, isp, mult, kx, ks, id, ifl, j, ilp, lfin, ic3p, ic3, iold, & - ikap, ind, kdif, ios, i0, k1, kdif1, jproj, isel, irr, matsize, maxsize, & - norbp, jnew - real*8 xinorm, del, angks, p2f, edge, omega, xk0, x, sinx, cosx, ww, wse, & - dum, sign, ratio, ratiop, xsedge, edg50, prefac, vrcorr, vicorr, xnorm, sfun - - complex*16, allocatable :: bmat(:,:,:,:,:,:) - logical,parameter :: correction_josh = .false. - logical CorrectOrbitalEnergies ! Set to true until this is done, then set to false. -! Debug: FDV -! dimension sh_eng(nex, 41) +subroutine xsect (ipr2, dx, x0, ri, ne, ne1, ik0, em, edge, & + ihole, emu, corr, dgc0, dpc0, jnew, & + ixc, lreal, rmt, rnrm, xmu, & + vi0, iPl, NPoles, Eps0, EGap, gamach, & + vtot, vvalgs, edens, dmag, edenvl, & + dgcn, dpcn, adgc, adpc, xsec, xsnorm, rkk, & + iz, xion, iunf, xnval, & + izstd, ifxc, eorb, kappa, iorb, l2lp, & + ipol, ispin, le2, angks, ptz, iph) !KJ + + USE IOMOD + use constants + use DimsMod, only: nex, nrptx, lx, MxPole, nspx=>nspu + USE SelfEnergyMod + ! right now the same self-energy is used for calculation of the central atom part (xsec) and dipole m.e. for + ! scattering (rkk). You may want to run xsect separately for xsec and for rkk, if you want to use different self-energy + ! for central and scattering parts. ala. fix later + + !implicit double precision (a-h, o-z) + implicit none + !Changed the dimensions to 40 to account for superheavy elements. Pavlo Baranov 07/2016 + integer, intent(in) :: ihole, ipr2, ne, ne1, ik0, l2lp, le2, ispin, ipol, iPl, & + ixc, lreal, iunf, izstd, iz, ifxc + integer, intent(inout) :: NPoles, iph !KJ I hope that is as intended?? + real*8, intent(in) :: dx, x0, rmt, rnrm, xmu, vi0, gamach, corr, xion, emu + character(len=50) :: fname + ! INPUT + ! dx, x0, ri(nr) + ! Loucks r-grid, ri=exp((i-1)*dx-x0) + ! ne, em(ne) number of energy points, real energy grid + ! edge chemical potential (energy for k=0) + ! ihole hole code + ! emu position of chemical potential in absorption specrum + ! dgc0(nr) dirac upper component, ground state hole orbital + ! dpc0(nr) dirac lower component, ground state hole orbital + ! ixc 0 Hedin-Lunqist + const real & imag part + ! 1 Dirac-Hara + const real & imag part + ! 2 ground state + const real & imag part + ! 3 Dirac-Hara + HL imag part + const real & imag part + ! 5 Dirac-Fock exchange with core electrons + + ! ixc=0 for valence electron density + ! lreal logical, true for real phase shifts only + ! rmt r muffin tin + ! xmu fermi level + ! vi0 const imag part to add to complex potential + ! gamach core hole lifetime + ! vtot(nr) total potential, including gsxc, final state + ! edens(nr) density, hole orbital, final state + ! dmag(251) density magnetization + ! edenvl valence charge density + ! dgcn(dpcn) large (small) dirac components for central atom + ! adgc(adpc) their development coefficients + + ! OUTPUT + ! xsec(ne) atomic absorption cross section to multiply \chi + ! (atomic background for XMCD) + ! xsnorm(ne) atomic absorption cross section (norm for XMCD) + ! rkk(ne, 8) normalized reduced matrix elements for construction + ! of termination matrix in genfmt. + + complex*16, intent(in) :: ptz(-1:1, -1:1) + complex*16, intent(in) :: em(nex) + real*8, intent(in) :: ri(nrptx), vtot(nrptx), edens(nrptx),dmag(nrptx) + real*8, intent(in) :: dgc0(nrptx), dpc0(nrptx), vvalgs(nrptx), edenvl(nrptx) + real*8, intent(in) :: dgcn(nrptx,41), dpcn(nrptx,41), eorb(41) + real*8, intent(in) :: adgc(10,41), adpc(10,41), xnval(41) + integer, intent(in) :: iorb(-5:4),kappa(41) + complex*16, intent(out) :: rkk(nex, 8), xsec(nex) + real*8, intent(out) :: xsnorm(nex) + + integer kiind(8), lind(8) !KJ renamed kind to kiind to avoid conflict with reserved name and its use in modules + real*8 xp(nrptx), xq(nrptx) + + ! work space for xcpot + real*8 vxcrmu(nrptx), vxcimu(nrptx), gsrel(nrptx) + real*8 vvxcrm(nrptx), vvxcim(nrptx) + + ! work space for fovrg + complex*16 p(nrptx), q(nrptx), pn(nrptx), qn(nrptx), fscf(nrptx) + complex*16 pp(nrptx), qp(nrptx), pnp(nrptx), qnp(nrptx) + ! storage for calculation of cross term (SPIN 1 only) + complex*16 xrcold(nrptx) , xncold(nrptx), yvec(nrptx,1) + + complex*16 p2, ck, xkmt, xkmtp + complex*16 pu, qu, dum1, factor + complex*16 xfnorm, xirf, xirf1 + complex*16 temp, aa, bb, cc, rkk1, rkk0, phold + complex*16 phx(8), ph0 + complex*16 eref, xm1, xm2, xm3, xm4, erefs(ne) + + complex*16 jl,jlp1,nl,nlp1 + complex*16 v(nrptx), vval(nrptx) + complex*16 xrc(nrptx), xnc(nrptx) + character*512 slog + logical ltrace + complex*16 xrhoce(nex), xrhopr(nex), chia(nex), cchi(nex) + real*8 omega1(nex), bf(0:2, nrptx) + + real*8 pat(nrptx),qat(nrptx) + complex*16 intr(nrptx),var(nrptx) + ! to pass energy levels and projected DOS + integer neg(41) + real*8 eng(nex, 41), rhoj(nex,41) + ! Josh - Added iPl switch for PLASMON card + ! - and WpCorr = Wi/Wp, Gamma, AmpFac + ! - to describe Im[eps^-1] + integer ipole, NData + real*8, allocatable :: Energy(:), Loss(:) + real*8 WpCorr(MxPole), Gamma(MxPole), AmpFac(MxPole), Eps0, EGap + character(LEN=10) ColumnLabels(20) + ! Josh END + integer imt, jri, jri1, inrm, jnrm, i, index, ifirst, ie, ncycle, ilast, ll, & + kinit, linit, isp, mult, kx, ks, id, ifl, j, ilp, lfin, ic3p, ic3, iold, & + ikap, ind, kdif, ios, i0, k1, kdif1, jproj, isel, irr, matsize, maxsize, & + norbp, jnew + real*8 xinorm, del, angks, p2f, edge, omega, xk0, x, sinx, cosx, ww, wse, & + dum, sign, ratio, ratiop, xsedge, edg50, prefac, vrcorr, vicorr, xnorm, sfun + + complex*16, allocatable :: bmat(:,:,:,:,:,:) + logical,parameter :: correction_josh = .false. + logical CorrectOrbitalEnergies ! Set to true until this is done, then set to false. + ! Debug: FDV + ! dimension sh_eng(nex, 41) - CorrectOrbitalEnergies = .TRUE. - allocate(bmat(-lx:lx,0:1,8, -lx:lx,0:1,8)) - - call setkap(ihole, kinit, linit) - -! set imt and jri (use general Loucks grid) -! rmt is between imt and jri (see function ii(r) in file xx.f) - imt = (log(rmt) + x0) / dx + 1 - jri = imt+1 - jri1 = jri+1 - if (jri1 .gt. nrptx) call par_stop('jri .gt. nrptx in phase') - -! Debug: FDV -! call correorb(iz, ihole, rmt, jri, dx,ri, p2f,edge, v, dgcn, dpcn, adgc, adpc, eorb, neg, sh_eng, rhoj, kappa, norbp) -! e_chsh = -sh_eng(1,1) -! write(6,fmt='(a,f12.8)') 'e_chsh: ', e_chsh - -! nesvi: define jnrm - inrm = (log(rnrm) + x0) / dx + 1 - jnrm = inrm + 1 -!write(*,*) 'x0,dx,rnrm,inrm,jnrm',x0,dx,rnrm,inrm,jnrm -! We'll need later to normalize dipole matrix elements -! . NB, dgc and dpc are r*wave_fn, so use '0' in somm to get integral psi**2 r**2 dr. -! Square the dgc0 and dpc0 arrays before integrating. -! == xinorm. -! dgc and dpc should be normalized =1, check this here - do i = 1, nrptx - xp(i) = dpc0(i)**2 - xq(i) = dgc0(i)**2 - enddo -! nb, xinorm is used for exponent on input to somm - xinorm = 2*linit + 2 - call somm (ri, xp, xq, dx, xinorm, 0, jnrm) - del = abs (abs(xinorm) - 1) - if (del .gt. 1.e-2) then - write(slog,'(a,i8,1p2e13.5)') ' ihole, xinorm ', ihole , xinorm - call wlog(slog) -! if using real phase shifts, don't expect great results - if (lreal.lt.2) then - call wlog(' There may be convergence problems.') - call wlog(' Xinorm should be 1. If you set the RGRID, minor ' // & - 'interpolation errors that will not affect final results may occur.') - endif + CorrectOrbitalEnergies = .TRUE. + allocate(bmat(-lx:lx,0:1,8, -lx:lx,0:1,8)) + + call setkap(ihole, kinit, linit) + + ! set imt and jri (use general Loucks grid) + ! rmt is between imt and jri (see function ii(r) in file xx.f) + imt = (log(rmt) + x0) / dx + 1 + jri = imt+1 + jri1 = jri+1 + if (jri1 .gt. nrptx) call par_stop('jri .gt. nrptx in phase') + + ! Debug: FDV + ! call correorb(iz, ihole, rmt, jri, dx,ri, p2f,edge, v, dgcn, dpcn, adgc, adpc, eorb, neg, sh_eng, rhoj, kappa, norbp) + ! e_chsh = -sh_eng(1,1) + ! write(6,fmt='(a,f12.8)') 'e_chsh: ', e_chsh + + ! nesvi: define jnrm + inrm = (log(rnrm) + x0) / dx + 1 + jnrm = inrm + 1 + !write(*,*) 'x0,dx,rnrm,inrm,jnrm',x0,dx,rnrm,inrm,jnrm + ! We'll need later to normalize dipole matrix elements + ! . NB, dgc and dpc are r*wave_fn, so use '0' in somm to get integral psi**2 r**2 dr. + ! Square the dgc0 and dpc0 arrays before integrating. + ! == xinorm. + ! dgc and dpc should be normalized =1, check this here + do i = 1, nrptx + xp(i) = dpc0(i)**2 + xq(i) = dgc0(i)**2 + enddo + ! nb, xinorm is used for exponent on input to somm + xinorm = 2*linit + 2 + call somm (ri, xp, xq, dx, xinorm, 0, jnrm) + del = abs (abs(xinorm) - 1) + if (del .gt. 1.e-2) then + write(slog,'(a,i8,1p2e13.5)') ' ihole, xinorm ', ihole , xinorm + call wlog(slog) + ! if using real phase shifts, don't expect great results + if (lreal.lt.2) then + call wlog(' There may be convergence problems.') + call wlog(' Xinorm should be 1. If you set the RGRID, minor ' // & + 'interpolation errors that will not affect final results may occur.') + endif + endif + + ! use ixc for testing + index = ixc + ! Always use ground state self energy for xsection, quick fix + ! JJR, Jan 93 + ! change for testing broadened plasmon pole 6/93 + ! index = 2 + ! ALA found that it is better to use index=ixc and real part of + ! self-energy for atomic xsection. 12/96 + ltrace = .true. + call bcoef(kinit, ipol, ptz, le2, ltrace, ispin, angks, kiind, lind, bmat) + ! set spin index to use bmat + isp = 0 + if (ispin.eq.1) isp = nspx - 1 + + ! zero rkk and phx + phx(:)=dcmplx(0) + rkk(:,:)=dcmplx(0) + + ifirst = 0 + ! Josh - if PLASMON card is set, and using HL exc, + ! - read loss function from loss.dat and find poles etc. + ! Josh - Changing to read directly from loss.dat : 3/9/2010 + IF ( (iPl.gt.0).and.(ixc.eq.0) ) THEN + WpCorr(:) = -1.d30 + CALL OpenFl('loss.dat', FileStatus = 'OLD') + NData = NumberOfLines('loss.dat') + + ALLOCATE(Energy(NData), Loss(NData)) + + CALL ReadArrayData('loss.dat', Double1 = Energy, Double2 = Loss) + + !NPoles = 100 ! 100 poles should be enough for anything. Can add a card for this later. + !Eps0 = -2.d0 ! This will not do anything. Will add a card to set Eps0 later. + CALL MkExc(Energy, Loss, Eps0, WpCorr, AmpFac, NPoles) + + Gamma(:) = Gamma(:)/hart + WpCorr(:) = (WpCorr(:)/hart) / SQRT(3.d0/((3 / (4*pi*edens(jri+1))) ** third)**3) + CALL CloseFl('loss.dat') + DEALLOCATE(Energy,Loss) + END IF + ! Write wp as calculated from density to sigma.dat + ColumnLabels(:) = ' ' + ColumnLabels(1) = 'Rs' + ColumnLabels(2) = 'wp (eV)' + fname = 'mpse_atomic' + CALL WriteData(TRIM(fname)//'.dat', & + Double1 = (3 / (4*pi*edens(jri+1))) ** third, & + Double2 = SQRT(3.d0/((3 / (4*pi*edens(jri+1))) ** third)**3)*hart, & + ColumnLabels = ColumnLabels, WriteDataInHeader = .TRUE., & + Headers = (/ 'This file contains information about the self-energy.' /)) + + ! Josh END + ! ********************* START OF BIG LOOP OVER ENERGY MESH ************** + do ie = 1, ne + iph = 0 + ! Josh - xcpot now has new arguments: + ! - iPl, WpCorr, Gamma, AmpFac, EGap + call xcpot(iph, ie, index, lreal, ifirst, jri, em(ie), xmu, vtot, vvalgs, edens, dmag, edenvl, & + eref, v, vval, iPl, WpCorr, Gamma, AmpFac, EGap, vxcrmu, vxcimu, gsrel, vvxcrm, vvxcim, fname) + ! set the method to calculate atomic cross section + ! p2 is (complex momentum)**2 referenced to energy dep xc + p2 = em(ie) - eref + p2f = edge - dble(eref) + ck = sqrt (2*p2 + (p2*alphfs)**2) + xkmt = rmt * ck + IF ((index.LT.5).OR.(index.EQ.10).OR.(index.EQ.11).OR.(index.EQ.12).OR.(index.EQ.13).OR.(index.EQ.14).OR.(index.EQ.6).OR.(index.EQ.7)) THEN + ncycle = 0 + ELSE + ! fix later. maybe ncycle can be less + ncycle = 3 + ENDIF + ! if (mod(index,10) .lt. 5) then + ! ncycle = 0 + ! else + ! ! fix later. maybe ncycle can be less + ! ncycle = 3 + ! endif + omega = (dble(em(ie)) - edge) + emu + omega = max (omega, 0.001d0 / hart) + ! nesvi: add omega1(ie)- need it later + omega1(ie) = omega + + ! remember the bessel functions for multipole matrix elements + xk0 = omega * alphfs + ilast = jnrm+6 + if (ilast.lt.jnew) ilast = jnew + if (ilast.gt.nrptx) ilast = nrptx + do i = 1, ilast + temp = xk0 * ri(i) + if (abs(temp).lt.1.d0) then + ! use series expansion + do ll = 0,2 + call bjnser(temp,ll, xirf, dum1,1) + bf(ll,i) = dble(xirf) + enddo + else + ! use formula + x = dble(temp) + sinx = sin(x) + cosx = cos(x) + bf(0,i) = sinx/x + bf(1,i) = sinx/x**2 - cosx/x + bf(2,i) = sinx*(3/x**3-1/x) - 3*cosx/x**2 + endif + enddo + + ! notice for spin-dep case xsnorm and xsec are spin-dep + ! and kept separately (see call to xsect in subroutine xsph) + xsnorm(ie) = 0 + xsec(ie) = 0 + if (dble(em(ie)).lt.-10.d0) cycle !goto 400 + if (dimag(p2).le.0.d0 .and. dble(p2).le.0.d0) cycle !goto 400 + + ! matrix elements for multipole (E1,E2,M1) transitions + ! The terms up to (k/c)^2 for absorption are kept. + ! L3 edge: kdif=1 (3d5/2) kdif=2 (3d3/2), kdif=3(4s) + ! L2 edge: kdif=1 (no transition), 2 (4s), 3 (3d3/2) + do 350 mult = 0, 2 + if (mult.eq.0) then + ! E1 transitions + kx = 1 + ks = 2 + else + ! M1 transitions + kx = 1 + ks = 6 + ! E2 transitions + if (mult.eq.2) kx = 2 + endif + ! skip unnecessary calculations + if (mult.gt.0 .and. (mult.ne.le2)) goto 350 + + ! set ilast larger than jri for better interpolation for pu + ! also need 5 points after jri for irregular solution + ilast = jnrm + 6 + if (ilast.lt.jnew) ilast = jnew + + ! calculate screened dipole field + ww = dble(emu+p2-edge) + + if (mult.eq.0 .and. izstd.gt.0) then + ! In f' calculations this gets skipped by line 262 + ! Change to logical first. + !if (ie.eq.1) then + if(CorrectOrbitalEnergies) then + CorrectOrbitalEnergies = .FALSE. + call correorb(iz, ihole, rmt, jri, dx,ri, p2f,edge, v, dgcn, & + dpcn, adgc, adpc, eorb, neg, eng, rhoj, kappa, norbp, iph) !KJ + if(correction_josh) eng(1,:) = eng(1,:) + p2f - edge ! Josh added this line so that + ! correorb calculates orbital energies relative to absolute zero. + end if + maxsize = 1 + matsize = 0 + sfun = 1.d0 + call phiscf (ifxc, rmt, ilast, jri, p2, p2f, emu, dx, ri, v, & + edens, dgcn, dpcn, adgc, adpc, iz, ihole, neg, eng, rhoj, & + kappa, norbp, fscf, yvec, maxsize, matsize, sfun, iph) !KJ + wse = dble(p2-eng(1,ihole)) + else + fscf(:) = 1.d0 + wse = ww + endif + if(.false.) then + write(*,*) 'ie',ie + write(*,*) 'p2',p2 + write(*,*) 'p2f',p2f + write(*,*) 'em(ie)',em(ie) + write(*,*) 'eref',eref + write(*,*) 'alphfs',alphfs + write(*,*) 'ck',ck + write(*,*) 'omega',omega + write(*,*) 'eng(1,ihole)',eng(1,ihole) + write(*,*) 'edge',edge + write(*,*) 'emu',emu + write(*,*) 'ww',ww + write(*,*) 'wse',wse endif -! use ixc for testing - index = ixc -! Always use ground state self energy for xsection, quick fix -! JJR, Jan 93 -! change for testing broadened plasmon pole 6/93 -! index = 2 -! ALA found that it is better to use index=ixc and real part of -! self-energy for atomic xsection. 12/96 - ltrace = .true. - call bcoef(kinit, ipol, ptz, le2, ltrace, ispin, angks, kiind, lind, bmat) -! set spin index to use bmat - isp = 0 - if (ispin.eq.1) isp = nspx - 1 - -! zero rkk and phx - phx(:)=dcmplx(0) - rkk(:,:)=dcmplx(0) - - ifirst = 0 -! Josh - if PLASMON card is set, and using HL exc, -! - read loss function from loss.dat and find poles etc. - ! Josh - Changing to read directly from loss.dat : 3/9/2010 - IF ( (iPl.gt.0).and.(ixc.eq.0) ) THEN - WpCorr(:) = -1.d30 - CALL OpenFl('loss.dat', FileStatus = 'OLD') - NData = NumberOfLines('loss.dat') - - ALLOCATE(Energy(NData), Loss(NData)) - - CALL ReadArrayData('loss.dat', Double1 = Energy, Double2 = Loss) - - !NPoles = 100 ! 100 poles should be enough for anything. Can add a card for this later. - !Eps0 = -2.d0 ! This will not do anything. Will add a card to set Eps0 later. - CALL MkExc(Energy, Loss, Eps0, WpCorr, AmpFac, NPoles) - - Gamma(:) = Gamma(:)/hart - WpCorr(:) = (WpCorr(:)/hart) / SQRT(3.d0/((3 / (4*pi*edens(jri+1))) ** third)**3) - CALL CloseFl('loss.dat') - DEALLOCATE(Energy,Loss) - END IF - ! Write wp as calculated from density to sigma.dat - ColumnLabels(:) = ' ' - ColumnLabels(1) = 'Rs' - ColumnLabels(2) = 'wp (eV)' - CALL WriteData('mpse.dat', & - Double1 = (3 / (4*pi*edens(jri+1))) ** third, & - Double2 = SQRT(3.d0/((3 / (4*pi*edens(jri+1))) ** third)**3)*hart, & - ColumnLabels = ColumnLabels, WriteDataInHeader = .TRUE., & - Headers = (/ 'This file contains information about the self-energy.' /)) - -! Josh END - - -! ********************* START OF BIG LOOP OVER ENERGY MESH ************** - do ie = 1, ne - iph = 0 -! Josh - xcpot now has new arguments: -! - iPl, WpCorr, Gamma, AmpFac, EGap - call xcpot (iph, ie, index, lreal, ifirst, jri, em(ie), xmu, vtot, vvalgs, edens, dmag, edenvl, & - eref, v, vval, iPl, WpCorr, Gamma, AmpFac, EGap, vxcrmu, vxcimu, gsrel, vvxcrm, vvxcim) - -! set the method to calculate atomic cross section -! p2 is (complex momentum)**2 referenced to energy dep xc - p2 = em(ie) - eref - p2f = edge - dble(eref) - ck = sqrt (2*p2 + (p2*alphfs)**2) - xkmt = rmt * ck - if (mod(index,10) .lt. 5) then - ncycle = 0 - else -! fix later. maybe ncycle can be less - ncycle = 3 - endif - omega = (dble(em(ie)) - edge) + emu - omega = max (omega, 0.001d0 / hart) -! nesvi: add omega1(ie)- need it later - omega1(ie) = omega - -! remember the bessel functions for multipole matrix elements - xk0 = omega * alphfs - ilast = jnrm+6 - if (ilast.lt.jnew) ilast = jnew - if (ilast.gt.nrptx) ilast = nrptx - do i = 1, ilast - temp = xk0 * ri(i) - if (abs(temp).lt.1.d0) then -! use series expansion - do ll = 0,2 - call bjnser(temp,ll, xirf, dum1,1) - bf(ll,i) = dble(xirf) + + ! ww = 1 + ! ww = wse / ww + ww = sqrt(wse/ww) + + do 300 kdif = -kx, kx + if (omega.le.0.0) goto 300 + ind = kdif + ks + ikap = kiind (ind) + if (ikap .eq. 0) goto 300 + ! use l2lp =0 to include both transitions l-->l+/-1 + ! if (l2lp.ne.0) only dipole transitions are calculated. + ! l-->l+1 transitions + if (l2lp.eq. 1 .and. ((kinit.lt.0 .and. ind.ge.3) .or. (kinit.gt.0 .and. ind.ne.3)) ) goto 300 + ! l-->l-1 transitions + if (l2lp.eq.-1 .and. ((kinit.lt.0 .and. ind.ne.3) .or. (kinit.gt.0 .and. ind.ge.3)) ) goto 300 + + iold = 0 + ic3=0 + ! start cycle do ic3=0,1 + ! return for ic3=1 calculations only for |ispin|=1 + ! where the central atom cross terms are needed + 100 continue + irr = -1 + ! ic3p=1 to test K2Cr2O7 L3 XES + ic3p = ic3 + call dfovrg ( ncycle, ikap, rmt, ilast, jri, p2, dx, ri, v, vval, dgcn, dpcn, adgc, adpc, & + xnval, pu, qu, p, q, iz, ihole, xion, iunf, irr, ic3p, iph) !KJ iph + lfin = lind (ind) + ilp = lfin - 1 + if (ikap .lt. 0) ilp = lfin + 1 + call exjlnl (xkmt, lfin, jl, nl) + call exjlnl (xkmt, ilp, jlp1, nlp1) + call phamp(rmt,pu,qu, ck, jl,nl,jlp1,nlp1, ikap, ph0,temp) + + sign = -1.0 + if (ikap.gt.0) sign = 1.0 + factor = ck*alphfs + factor = sign * factor/(1+sqrt(1+factor**2)) + dum1 = 1/ sqrt(1+factor**2) + xfnorm = 1 / temp *dum1 + ! normalization factor + ! xfnorm = dum1*rmt*(jl*cos(delta) - nl*sin(delta))/ Rl(rmt) + ! dum1 is relativistic correction to normalization + ! normalize regular solution + do i = 1,ilast + p(i)=p(i)*xfnorm + q(i)=q(i)*xfnorm + enddo + + ! calculate xirf including fscf - TDLDA result + do id = 1, 2 + if((izstd.le.0).AND.(id.eq.2)) EXIT ! Josh added to fix glitch in xmcd. + if (id.eq.1) then + do j = 1,ilast + pp(j) = p(j)*dble(fscf(j)) + qp(j) = q(j)*dble(fscf(j)) enddo + else ! id.eq.2 + do j = 1,ilast + pp(j) = p(j)*dimag(fscf(j)) + qp(j) = q(j)*dimag(fscf(j)) + enddo + endif + ifl = 1 + if (izstd.gt.0) ifl = -1 + xirf1 = 0 + call radint(ifl, mult, bf, kinit, dgc0,dpc0, ikap, pp,qp, pn,qn,ri,dx, ilast,iold, xrc,xnc, xrcold,xncold, xirf1) + ! if (ifl.lt.0) xirf1 = xirf1 * xk0 * ww + if (ifl.lt.0) xirf1 = xirf1 * xk0 + if (id.eq.1) then + xirf = xirf1 else -! use formula - x = dble(temp) - sinx = sin(x) - cosx = cos(x) - bf(0,i) = sinx/x - bf(1,i) = sinx/x**2 - cosx/x - bf(2,i) = sinx*(3/x**3-1/x) - 3*cosx/x**2 + if (abs(xirf) .eq. 0.d0) then + xirf = xirf1 + elseif (abs(xirf1) .eq. 0.d0) then + xirf = xirf + elseif (abs(xirf1) .lt. abs(xirf)) then + dum = abs(xirf1) / abs(xirf) + xirf = xirf * sqrt(1.d0 + dum**2) + else + dum = abs(xirf) / abs(xirf1) + xirf = xirf1 * sqrt(1.d0 + dum**2) + endif endif enddo -! notice for spin-dep case xsnorm and xsec are spin-dep -! and kept separately (see call to xsect in subroutine xsph) - xsnorm(ie) = 0 - xsec(ie) = 0 - if (dble(em(ie)).lt.-10.d0) cycle !goto 400 - if (dimag(p2).le.0.d0 .and. dble(p2).le.0.d0) cycle !goto 400 - -! matrix elements for multipole (E1,E2,M1) transitions -! The terms up to (k/c)^2 for absorption are kept. -! L3 edge: kdif=1 (3d5/2) kdif=2 (3d3/2), kdif=3(4s) -! L2 edge: kdif=1 (no transition), 2 (4s), 3 (3d3/2) - do 350 mult = 0, 2 - if (mult.eq.0) then -! E1 transitions - kx = 1 - ks = 2 - else -! M1 transitions - kx = 1 - ks = 6 -! E2 transitions - if (mult.eq.2) kx = 2 - endif -! skip unnecessary calculations - if (mult.gt.0 .and. (mult.ne.le2)) goto 350 - -! set ilast larger than jri for better interpolation for pu -! also need 5 points after jri for irregular solution - ilast = jnrm + 6 - if (ilast.lt.jnew) ilast = jnew - -!c calculate screened dipole field - ww = dble(emu+p2-edge) - - if (mult.eq.0 .and. izstd.gt.0) then - ! In f' calculations this gets skipped by line 262 - ! Change to logical first. - !if (ie.eq.1) then - if(CorrectOrbitalEnergies) then - CorrectOrbitalEnergies = .FALSE. - call correorb(iz, ihole, rmt, jri, dx,ri, p2f,edge, v, dgcn, & - dpcn, adgc, adpc, eorb, neg, eng, rhoj, kappa, norbp, iph) !KJ - if(correction_josh) eng(1,:) = eng(1,:) + p2f - edge ! Josh added this line so that - ! correorb calculates orbital energies relative to absolute zero. - end if - maxsize = 1 - matsize = 0 - sfun = 1.d0 - call phiscf (ifxc, rmt, ilast, jri, p2, p2f, emu, dx, ri, v, & - edens, dgcn, dpcn, adgc, adpc, iz, ihole, neg, eng, rhoj, & - kappa, norbp, fscf, yvec, maxsize, matsize, sfun, iph) !KJ - wse = dble(p2-eng(1,ihole)) - else - fscf(:) = 1.d0 - wse = ww + ! note that for real potential xirf is real or reduced matrix + ! element for dipole transition is pure imaginary. + if (ic3.eq.0) then + ! can remember only E2 or M1 matrix elements + if (mult.eq.0 .or. le2.eq.mult) then + rkk(ie,ind)=xirf + phx(ind) = ph0 endif -if(.false.) then -write(*,*) 'ie',ie -write(*,*) 'p2',p2 -write(*,*) 'p2f',p2f -write(*,*) 'em(ie)',em(ie) -write(*,*) 'eref',eref -write(*,*) 'alphfs',alphfs -write(*,*) 'ck',ck -write(*,*) 'omega',omega -write(*,*) 'eng(1,ihole)',eng(1,ihole) -write(*,*) 'edge',edge -write(*,*) 'emu',emu -write(*,*) 'ww',ww -write(*,*) 'wse',wse -endif - - -! ww = 1 -! ww = wse / ww - ww = sqrt(wse/ww) - - do 300 kdif = -kx, kx - if (omega.le.0.0) goto 300 - ind = kdif + ks - ikap = kiind (ind) - if (ikap .eq. 0) goto 300 -! use l2lp =0 to include both transitions l-->l+/-1 -! if (l2lp.ne.0) only dipole transitions are calculated. -! l-->l+1 transitions - if (l2lp.eq. 1 .and. ((kinit.lt.0 .and. ind.ge.3) .or. (kinit.gt.0 .and. ind.ne.3)) ) goto 300 -! l-->l-1 transitions - if (l2lp.eq.-1 .and. ((kinit.lt.0 .and. ind.ne.3) .or. (kinit.gt.0 .and. ind.ge.3)) ) goto 300 - - iold = 0 - ic3=0 -! start cycle do ic3=0,1 -! return for ic3=1 calculations only for |ispin|=1 -! where the central atom cross terms are needed - 100 continue - irr = -1 -! ic3p=1 to test K2Cr2O7 L3 XES - ic3p = ic3 - call dfovrg ( ncycle, ikap, rmt, ilast, jri, p2, dx, ri, v, vval, dgcn, dpcn, adgc, adpc, & - xnval, pu, qu, p, q, iz, ihole, xion, iunf, irr, ic3p, iph) !KJ iph - lfin = lind (ind) - ilp = lfin - 1 - if (ikap .lt. 0) ilp = lfin + 1 - call exjlnl (xkmt, lfin, jl, nl) - call exjlnl (xkmt, ilp, jlp1, nlp1) - call phamp(rmt,pu,qu, ck, jl,nl,jlp1,nlp1, ikap, ph0,temp) - - sign = -1.0 - if (ikap.gt.0) sign = 1.0 - factor = ck*alphfs - factor = sign * factor/(1+sqrt(1+factor**2)) - dum1 = 1/ sqrt(1+factor**2) - xfnorm = 1 / temp *dum1 -! normalization factor -! xfnorm = dum1*rmt*(jl*cos(delta) - nl*sin(delta))/ Rl(rmt) -! dum1 is relativistic correction to normalization -! normalize regular solution - do i = 1,ilast - p(i)=p(i)*xfnorm - q(i)=q(i)*xfnorm - enddo + ! for f' want to include both E2 and M1 for xsnorm and xsec + ! but now only one of them is included (fix later) + xsnorm(ie)=xsnorm(ie) + ( dble(xirf)**2 + dimag(xirf)**2 )/ (2*kx+1) + aa = - coni*xirf**2 + xsec(ie) = xsec(ie) - aa * bmat(0,isp,ind, 0,isp,ind) + elseif (iold.eq.1) then + rkk1=xirf + phold = ph0 + elseif (iold.eq.2) then + rkk0=xirf + endif -! calculate xirf including fscf - TDLDA result - do id = 1, 2 - if((izstd.le.0).AND.(id.eq.2)) EXIT ! Josh added to fix glitch in xmcd. - if (id.eq.1) then - do j = 1,ilast - pp(j) = p(j)*dble(fscf(j)) - qp(j) = q(j)*dble(fscf(j)) - enddo - else ! id.eq.2 - do j = 1,ilast - pp(j) = p(j)*dimag(fscf(j)) - qp(j) = q(j)*dimag(fscf(j)) - enddo - endif - ifl = 1 - if (izstd.gt.0) ifl = -1 - xirf1 = 0 - call radint(ifl, mult, bf, kinit, dgc0,dpc0, ikap, pp,qp, pn,qn,ri,dx, ilast,iold, xrc,xnc, xrcold,xncold, xirf1) -! if (ifl.lt.0) xirf1 = xirf1 * xk0 * ww - if (ifl.lt.0) xirf1 = xirf1 * xk0 - if (id.eq.1) then - xirf = xirf1 - else - if (abs(xirf) .eq. 0.d0) then - xirf = xirf1 - elseif (abs(xirf1) .eq. 0.d0) then - xirf = xirf - elseif (abs(xirf1) .lt. abs(xirf)) then - dum = abs(xirf1) / abs(xirf) - xirf = xirf * sqrt(1.d0 + dum**2) - else - dum = abs(xirf) / abs(xirf1) - xirf = xirf1 * sqrt(1.d0 + dum**2) - endif - endif - enddo + ! get irregular solution and atomic cross-section xsec + ! find irregular solution + if(dimag(em(ie)).gt.0.d0) then + irr = 1 + ! set pu, qu - initial condition for irregular solution + pu = (nl*cos(ph0)+jl*sin(ph0)) *rmt * dum1 + qu=(nlp1*cos(ph0)+jlp1*sin(ph0))*factor *rmt * dum1 + + ! test on bessel functions + ! if (ikap.gt.0) print*,'test1',xkmt**2*(jl*nlp1-nl*jlp1) + + call dfovrg (ncycle, ikap, rmt, ilast, jri, p2, dx, ri, v,vval, dgcn, dpcn, adgc, adpc, & + xnval, pu, qu, pn, qn, iz, ihole, xion, iunf, irr, ic3p, iph) !KJ iph + ! set N- irregular solution , which is outside + ! N=(nlp1*cos(ph0)+jlp1*sin(ph0))*factor *rmt * dum1 + ! N = i*R - H*exp(i*ph0) + temp = exp(coni*ph0) + do i = 1, ilast + pn(i) = coni * p(i) - temp * pn(i) + qn(i) = coni * q(i) - temp * qn(i) + enddo + else + pn(1:ilast) = 0 + qn(1:ilast) = 0 + endif -! note that for real potential xirf is real or reduced matrix -! element for dipole transition is pure imaginary. - if (ic3.eq.0) then -! can remember only E2 or M1 matrix elements - if (mult.eq.0 .or. le2.eq.mult) then - rkk(ie,ind)=xirf - phx(ind) = ph0 - endif -! for f' want to include both E2 and M1 for xsnorm and xsec -! but now only one of them is included (fix later) - xsnorm(ie)=xsnorm(ie) + ( dble(xirf)**2 + dimag(xirf)**2 )/ (2*kx+1) - aa = - coni*xirf**2 - xsec(ie) = xsec(ie) - aa * bmat(0,isp,ind, 0,isp,ind) - elseif (iold.eq.1) then - rkk1=xirf - phold = ph0 - elseif (iold.eq.2) then - rkk0=xirf - endif + ! combine regular and irregular solution into the central atom absorption coefficient xsec (mu = dimag(xsec)) + ! thus for real energy dimag(xsec)=xsnorm + ! also include TDLDA effects + do id = 1, 2 + if((izstd.le.0).AND.(id.eq.2)) EXIT ! Josh - added to avoid glitch in XMCD. + if (id.eq.1) then + do j = 1,ilast + pp(j) = p(j)*dble(fscf(j)) + qp(j) = q(j)*dble(fscf(j)) + pnp(j) = pn(j)*dble(fscf(j)) + qnp(j) = qn(j)*dble(fscf(j)) + enddo + else + do j = 1,ilast + pp(j) = p(j)*dimag(fscf(j)) + qp(j) = q(j)*dimag(fscf(j)) + pnp(j) = pn(j)*dimag(fscf(j)) + qnp(j) = qn(j)*dimag(fscf(j)) + enddo + endif -! get irregular solution and atomic cross-section xsec -! find irregular solution - if(dimag(em(ie)).gt.0.d0) then - irr = 1 -! set pu, qu - initial condition for irregular solution - pu = (nl*cos(ph0)+jl*sin(ph0)) *rmt * dum1 - qu=(nlp1*cos(ph0)+jlp1*sin(ph0))*factor *rmt * dum1 - -! test on bessel functions -! if (ikap.gt.0) print*,'test1',xkmt**2*(jl*nlp1-nl*jlp1) - - call dfovrg (ncycle, ikap, rmt, ilast, jri, p2, dx, ri, v,vval, dgcn, dpcn, adgc, adpc, & - xnval, pu, qu, pn, qn, iz, ihole, xion, iunf, irr, ic3p, iph) !KJ iph -! set N- irregular solution , which is outside -! N=(nlp1*cos(ph0)+jlp1*sin(ph0))*factor *rmt * dum1 -! N = i*R - H*exp(i*ph0) - temp = exp(coni*ph0) - do i = 1, ilast - pn(i) = coni * p(i) - temp * pn(i) - qn(i) = coni * q(i) - temp * qn(i) - enddo + ! TDLDA theory is written for the r-form of matrix elements + ! so one might want to use ifl=-1,-2 for these calculations + ! on the other hand want ifl=1,2 for DANES calculations + ! since it is more reliable at high energies and gives better results for Cu test. + ifl = 2 + if (izstd.gt.0) ifl = -2 + + call radint(ifl,mult, bf, kinit, dgc0, dpc0, ikap, pp, qp, & + pnp, qnp, ri,dx, ilast,iold, xrc, xnc, xrcold, xncold, xirf1) + if (ifl.lt.0) xirf1 = xirf1 * xk0**2 * ww**2 + if (id.eq.1) then + xirf = xirf1 + isel=0 + else + if (abs(xirf) .eq. 0.d0) then + xirf = xirf1 + isel=1 + elseif (abs(xirf1) .eq. 0.d0) then + xirf = xirf + isel=2 + elseif (abs(xirf1) .lt. abs(xirf)) then + dum = abs(xirf1) / abs(xirf) + xirf = xirf * sqrt(1.d0 + dum**2) + isel=3 else - pn(1:ilast) = 0 - qn(1:ilast) = 0 + dum = abs(xirf) / abs(xirf1) + xirf = xirf1 * sqrt(1.d0 + dum**2) + isel=4 endif + endif + !write(*,'(3i4,2e12.4,i4,6e12.4)') ie,id,ic3,xirf1,isel,xirf,ww,xk0 + enddo + if (ic3.eq.0) xsec(ie) = xsec(ie) - xirf * bmat(0,isp,ind, 0,isp,ind) + + + ! ------start of density of states part------------- + ! added by nesvi 07/12/00 + ! Calculate rhoc00 and rho_projected on the same grid as xsect. Need this to calculate the smooth atomic ratio rho_0/mu_0 or rho_proj/mu_0. + ! The atomic functions are normalized to 1 inside Norman radius. This procedure can be called 'Renormalized atomic sphere method". + ! It gives very reasonable numbers for hole counts. In order to get Mulliken counts one can extend integration limits till very + ! large R, but it's currently not recommended because of the problems with the wave function's tails above Rnm. + jproj = iorb(ikap) + if (jproj.eq.0 .and. ikap.lt.0) jproj = iorb(-ikap-1) + kdif1 = -1 + if(kinit.gt.0) kdif1 = 1 + + if (kdif .eq. kdif1 .and. ic3.eq.0 .and. jproj.gt.0) then + ! calculate rhoc00 (rho_0) -! combine regular and irregular solution into the central atom absorption coefficient xsec (mu = dimag(xsec)) -! thus for real energy dimag(xsec)=xsnorm -! also include TDLDA effects - do id = 1, 2 - if((izstd.le.0).AND.(id.eq.2)) EXIT ! Josh - added to avoid glitch in XMCD. - if (id.eq.1) then - do j = 1,ilast - pp(j) = p(j)*dble(fscf(j)) - qp(j) = q(j)*dble(fscf(j)) - pnp(j) = pn(j)*dble(fscf(j)) - qnp(j) = qn(j)*dble(fscf(j)) - enddo - else - do j = 1,ilast - pp(j) = p(j)*dimag(fscf(j)) - qp(j) = q(j)*dimag(fscf(j)) - pnp(j) = pn(j)*dimag(fscf(j)) - qnp(j) = qn(j)*dimag(fscf(j)) - enddo - endif - -! TDLDA theory is written for the r-form of matrix elements -! so one might want to use ifl=-1,-2 for these calculations -! on the other hand want ifl=1,2 for DANES calculations -! since it is more reliable at high energies and gives better results for Cu test. - ifl = 2 - if (izstd.gt.0) ifl = -2 - - call radint(ifl,mult, bf, kinit, dgc0, dpc0, ikap, pp, qp, & - pnp, qnp, ri,dx, ilast,iold, xrc, xnc, xrcold, xncold, xirf1) - if (ifl.lt.0) xirf1 = xirf1 * xk0**2 * ww**2 - if (id.eq.1) then - xirf = xirf1 - isel=0 - else - if (abs(xirf) .eq. 0.d0) then - xirf = xirf1 - isel=1 - elseif (abs(xirf1) .eq. 0.d0) then - xirf = xirf - isel=2 - elseif (abs(xirf1) .lt. abs(xirf)) then - dum = abs(xirf1) / abs(xirf) - xirf = xirf * sqrt(1.d0 + dum**2) - isel=3 - else - dum = abs(xirf) / abs(xirf1) - xirf = xirf1 * sqrt(1.d0 + dum**2) - isel=4 - endif - endif - !write(*,'(3i4,2e12.4,i4,6e12.4)') ie,id,ic3,xirf1,isel,xirf,ww,xk0 - enddo - if (ic3.eq.0) xsec(ie) = xsec(ie) - xirf * bmat(0,isp,ind, 0,isp,ind) + temp = (2*lfin+1.0d0) / (1+factor**2) /pi *4*ck /hart + do i = 1, ilast + xrc(i) = pn(i)*p(i) - coni*p(i)*p(i) + qn(i)*q(i) - coni*q(i)*q(i) + enddo + xirf = 1 + ! integration is till Norman radius, not Rmt as in xsect + i0 = jnrm + 1 + call csomm2 (ri, xrc, dx, xirf, rnrm, i0) + xrhoce(ie) = - xirf * temp + + ! calculate rho_projected: + ! pat, qat - atomic functions that we make projection on. + do i=1,nrptx + pat(i) = dgcn(i,jproj) + qat(i) = dpcn(i,jproj) + enddo -! ------start of density of states part------------- -! added by nesvi 07/12/00 -! Calculate rhoc00 and rho_projected on the same grid as xsect. Need this to calculate the smooth atomic ratio rho_0/mu_0 or rho_proj/mu_0. -! The atomic functions are normalized to 1 inside Norman radius. This procedure can be called 'Renormalized atomic sphere method". -! It gives very reasonable numbers for hole counts. In order to get Mulliken counts one can extend integration limits till very -! large R, but it's currently not recommended because of the problems with the wave function's tails above Rnm. - + ! normalize pat and qat in the Norman radius sphere: =1, (renormalized atomic sphere method) + do i = 1, ilast + xp(i) = pat(i)**2 + qat(i)**2 + xq(i) = 0 + enddo + !nb, xinorm is used for exponent on input to somm + xinorm = 2*lfin + 2 + call somm2 (ri, xp, dx, xinorm, rnrm, 0, i0) + ! call somm (ri, xp, xq, dx, xinorm, 0, jnrm) + + xinorm = sqrt(xinorm) + do i=1,nrptx + pat(i) = pat(i) / xinorm + qat(i) = qat(i) / xinorm + enddo + + !calculate overlap integral between f and atomic function (integral Rl(r)*Psi_at(r)dr from 0 till r'). + !intr(i) is that overlap integral. Later it will be multiplied by pr(i)*Psi_at(r') and integrated till Norman radius. + do i=1,ilast + var(i)=pat(i)*p(i)+qat(i)*q(i) + ! factor of 2 -integration r< r> -->2 r r' + enddo - jproj = iorb(ikap) - if (jproj.eq.0 .and. ikap.lt.0) jproj = iorb(-ikap-1) - kdif1 = -1 - if(kinit.gt.0) kdif1 = 1 - - if (kdif .eq. kdif1 .and. ic3.eq.0 .and. jproj.gt.0) then -! calculate rhoc00 (rho_0) - - temp = (2*lfin+1.0d0) / (1+factor**2) /pi *4*ck /hart - do i = 1, ilast - xrc(i) = pn(i)*p(i) - coni*p(i)*p(i) + qn(i)*q(i) - coni*q(i)*q(i) - enddo - xirf = 1 -! integration is till Norman radius, not Rmt as in xsect - i0 = jnrm + 1 - call csomm2 (ri, xrc, dx, xirf, rnrm, i0) - xrhoce(ie) = - xirf * temp - -! calculate rho_projected: - -! pat, qat - atomic functions that we make projection on. - do i=1,nrptx - pat(i) = dgcn(i,jproj) - qat(i) = dpcn(i,jproj) - enddo - -! normalize pat and qat in the Norman radius sphere: =1, (renormalized atomic sphere method) - do i = 1, ilast - xp(i) = pat(i)**2 + qat(i)**2 - xq(i) = 0 - enddo - !nb, xinorm is used for exponent on input to somm - xinorm = 2*lfin + 2 - call somm2 (ri, xp, dx, xinorm, rnrm, 0, i0) -! call somm (ri, xp, xq, dx, xinorm, 0, jnrm) + ! integration by trapezoid method + intr(1)=var(1)*ri(1) + do i=2,ilast + intr(i)=intr(i-1)+ (var(i)+var(i-1))*(ri(i)-ri(i-1)) + enddo + + ! now calculate rho_projected - xrhopr + temp = (2*lfin+1.0d0) / (1+factor**2) /pi *4*ck /hart + ! temp = abs(ikap) / (1+factor**2) /pi *4*ck /hart + do i = 1, ilast + xrc(i) = pn(i)*pat(i)*intr(i)+ qn(i)*qat(i)*intr(i) + xrc(i) = xrc(i) - coni*(p(i)*pat(i)*intr(i) + q(i)*qat(i)*intr(i)) + enddo + + xirf = 1 + call csomm2 (ri, xrc, dx, xirf, rnrm, i0) + xrhopr(ie) = - xirf * temp - xinorm = sqrt(xinorm) - do i=1,nrptx - pat(i) = pat(i) / xinorm - qat(i) = qat(i) / xinorm - enddo - - !calculate overlap integral between f and atomic function (integral Rl(r)*Psi_at(r)dr from 0 till r'). - !intr(i) is that overlap integral. Later it will be multiplied by pr(i)*Psi_at(r') and integrated till Norman radius. - do i=1,ilast - var(i)=pat(i)*p(i)+qat(i)*q(i) -! factor of 2 -integration r< r> -->2 r r' - enddo - -! integration by trapezoid method - intr(1)=var(1)*ri(1) - do i=2,ilast - intr(i)=intr(i-1)+ (var(i)+var(i-1))*(ri(i)-ri(i-1)) - enddo - -! now calculate rho_projected - xrhopr - temp = (2*lfin+1.0d0) / (1+factor**2) /pi *4*ck /hart -! temp = abs(ikap) / (1+factor**2) /pi *4*ck /hart - do i = 1, ilast - xrc(i) = pn(i)*pat(i)*intr(i)+ qn(i)*qat(i)*intr(i) - xrc(i) = xrc(i) - coni*(p(i)*pat(i)*intr(i) + q(i)*qat(i)*intr(i)) - enddo - - xirf = 1 - call csomm2 (ri, xrc, dx, xirf, rnrm, i0) - xrhopr(ie) = - xirf * temp + endif + ! ----------end of density of states part--- + + + if (iold.gt.0) then + ! calculate cross term contribution to XMCD + ! in both cases coupling between neighbors + ! need to remove SO interaction (ic3=1) in order to avoid unphysical peak in Gd XMCD. a.l. ankudinov + k1 = ind - 1 + if (k1.ge.1 .and.k1.le.8) then + if (lind(k1).eq.lind(ind) .and. lind(k1).gt.0) then + aa = exp( coni*(ph0 - phold)) + bb = 1/aa + cc = - ( bmat(0,isp,k1, 0,isp,ind) + bmat(0,isp,ind, 0,isp,k1) ) / 2.d0 + xsec(ie) = xsec(ie) - coni * rkk1 * rkk0 * (bb+aa) * cc + ! combine regular and irregular solution into the central atom absorption coefficient (mu=dimag(xsec)) + ! thus for real energy dimag(xsec)=xsnorm + call radint (3, mult, bf, kinit, dgc0, dpc0, ikap, p, q, pn, & + qn, ri, dx, ilast, iold, xrc, xnc, xrcold, xncold, xirf) + xsec(ie) = xsec(ie) + xirf * cc * bb + call radint (4, mult, bf, kinit, dgc0, dpc0, ikap, p, q, pn, & + qn, ri, dx, ilast, iold, xrc, xnc, xrcold, xncold, xirf) + xsec(ie) = xsec(ie) + xirf * cc * aa endif -! ----------end of density of states part--- - - - if (iold.gt.0) then -! calculate cross term contribution to XMCD -! in both cases coupling between neighbors -! need to remove SO interaction (ic3=1) in order to avoid unphysical peak in Gd XMCD. a.l. ankudinov - k1 = ind - 1 - if (k1.ge.1 .and.k1.le.8) then - if (lind(k1).eq.lind(ind) .and. lind(k1).gt.0) then - aa = exp( coni*(ph0 - phold)) - bb = 1/aa - cc = - ( bmat(0,isp,k1, 0,isp,ind) + bmat(0,isp,ind, 0,isp,k1) ) / 2.d0 - xsec(ie) = xsec(ie) - coni * rkk1 * rkk0 * (bb+aa) * cc -! combine regular and irregular solution into the central atom absorption coefficient (mu=dimag(xsec)) -! thus for real energy dimag(xsec)=xsnorm - call radint (3, mult, bf, kinit, dgc0, dpc0, ikap, p, q, pn, & - qn, ri, dx, ilast, iold, xrc, xnc, xrcold, xncold, xirf) - xsec(ie) = xsec(ie) + xirf * cc * bb - - call radint (4, mult, bf, kinit, dgc0, dpc0, ikap, p, q, pn, & - qn, ri, dx, ilast, iold, xrc, xnc, xrcold, xncold, xirf) - xsec(ie) = xsec(ie) + xirf * cc * aa - endif - endif - endif -! end of |ispin=1| cross term calculations - -! prepare for ic3=1 cross term calculations if needed - if (ic3.eq.0 .and. abs(ispin).eq.1) then - iold = 0 - ! Josh - Write equivalent in if elseif (easier to see what is happening). - ! if orbital angular momentum (lind) is not zero - if (lind(ind).gt.0) then - if(ind.eq.1) then - ! If this is the first transition of this angular momentum, save xrc and xnc from radint. - k1 = ind + 1 - if (lind(k1).eq.lind(ind)) iold = 1 - else - ! Otherwise, use the old result since xrc and xnc should be the same. - k1 = ind - 1 - if(lind(k1).eq.lind(ind)) iold = 2 - end if + endif + endif + ! end of |ispin=1| cross term calculations + + ! prepare for ic3=1 cross term calculations if needed + if (ic3.eq.0 .and. abs(ispin).eq.1) then + iold = 0 + ! Josh - Write equivalent in if elseif (easier to see what is happening). + ! if orbital angular momentum (lind) is not zero + if (lind(ind).gt.0) then + if(ind.eq.1) then + ! If this is the first transition of this angular momentum, save xrc and xnc from radint. + k1 = ind + 1 + if (lind(k1).eq.lind(ind)) iold = 1 + else + ! Otherwise, use the old result since xrc and xnc should be the same. + k1 = ind - 1 + if(lind(k1).eq.lind(ind)) iold = 2 end if - !if (ind.lt.8 .and. lind(ind).gt.0) then - ! k1 = ind + 1 - ! if (lind(k1).eq.lind(ind)) iold = 1 - !endif - !if (ind.gt.1 .and. lind(ind).gt.0) then - ! k1 = ind - 1 - ! if (lind(k1).eq.lind(ind)) iold = 2 - !endif -! need to remove SO interaction to calculate cross term -! big effect for Gd XMCD calculations - if (iold.gt.0) then -! repeat calculation for current kdif with SO turned off - ic3 = 1 - goto 100 - endif - endif + end if + !if (ind.lt.8 .and. lind(ind).gt.0) then + ! k1 = ind + 1 + ! if (lind(k1).eq.lind(ind)) iold = 1 + !endif + !if (ind.gt.1 .and. lind(ind).gt.0) then + ! k1 = ind - 1 + ! if (lind(k1).eq.lind(ind)) iold = 2 + !endif + ! need to remove SO interaction to calculate cross term + ! big effect for Gd XMCD calculations + if (iold.gt.0) then + ! repeat calculation for current kdif with SO turned off + ic3 = 1 + goto 100 + endif + endif - 300 continue - 350 continue - if (omega.gt.0.0) then -! prefac = (8 * pi / 3) * alphfs * omega -- nonrelativistic -! relativistic is (for alpha form) - prefac = 4 * pi * alpinv / omega * bohr**2 - xsnorm(ie) = xsnorm(ie) * prefac * 2*abs(ck) - xnorm= sqrt( xsnorm(ie) ) - xsec(ie) = xsec(ie) * prefac* 2*ck - -! put complex sqrt(prefactor) into reduced matrix elements rkk - ck = sqrt ( prefac * (2*ck)) -! guarantee that we have the right root - if (dimag(ck) .lt. 0) ck = -ck -! add central atom phase shift here. - do kdif = 1 , 8 - rkk(ie,kdif)= rkk(ie,kdif) * ck/xnorm * exp(coni*phx(kdif)) - enddo - endif + 300 continue + 350 continue + if (omega.gt.0.0) then + ! prefac = (8 * pi / 3) * alphfs * omega -- nonrelativistic + ! relativistic is (for alpha form) + prefac = 4 * pi * alpinv / omega * bohr**2 + xsnorm(ie) = xsnorm(ie) * prefac * 2*abs(ck) + xnorm= sqrt( xsnorm(ie) ) + xsec(ie) = xsec(ie) * prefac* 2*ck + + ! put complex sqrt(prefactor) into reduced matrix elements rkk + ck = sqrt ( prefac * (2*ck)) + ! guarantee that we have the right root + if (dimag(ck) .lt. 0) ck = -ck + ! add central atom phase shift here. + do kdif = 1 , 8 + rkk(ie,kdif)= rkk(ie,kdif) * ck/xnorm * exp(coni*phx(kdif)) enddo -! ********************* END OF BIG LOOP OVER ENERGY MESH ************** - - - CALL CloseFl('mpse.dat') - -! *** optionally write ratio.dat and ratiop.dat ***** - if (ipr2.ge.3) then -! calculate mu_0/rho_0 for XMCD normalization. - chia(1:ne) = 0 - vrcorr = 0 - vicorr = 0 - call xscorr(1, em, ne1, ne, ik0, xrhoce,xsnorm,chia,vrcorr, vicorr, cchi) - do ie = 1, ne1 - xrhoce(ie) = coni* dimag(xrhoce(ie)+cchi(ie)) - enddo - call xscorr(1, em, ne1, ne, ik0, xrhopr,xsnorm,chia,vrcorr, vicorr, cchi) - do ie = 1, ne1 - xrhopr(ie) = coni* dimag(xrhopr(ie)+cchi(ie)) - enddo - call xscorr(1, em, ne1, ne, ik0, xsec,xsnorm,chia,vrcorr, vicorr, cchi) - do ie = 1, ne1 - cchi(ie) = coni* dimag(xsec(ie)+cchi(ie)) - enddo - open(unit=3,file='ratio.dat',status='unknown', iostat=ios) - open(unit=4,file='ratiop.dat',status='unknown', iostat=ios) -! normalize to xsec at 50 ev above edge - edg50 = emu +50.0 / hart - call terp (omega1, xsnorm, ne1, 1, edg50, xsedge) - write(3,440) xsedge, emu * hart - 440 format ('# Normalization factor:', e12.4,' Angstrom**2. Fermi level at ', f7.1, ' eV.') - write(3,450) - 450 format ('# Energy rho_0 mu_0 rho_0/mu_0 ') - write(4,440) xsedge, emu * hart - write(4,455) - 455 format ('# Energy rho_proj mu_0 rho_proj/mu_0',' mu_deloc ') - - do ie=1,ne1 - if (dimag(cchi(ie)).eq.0.d0 .and. ie.lt.ik0) then - cchi(ie)=cchi(ik0) - xrhoce(ie)=xrhoce(ik0) - xrhopr(ie)=xrhopr(ik0) - endif - ratio = dimag(xrhoce(ie)) / dimag(cchi(ie)) * xsedge - ratiop = dimag(xrhopr(ie)) / dimag(cchi(ie)) * xsedge - - write(3,460) dble(em(ie))*hart, dimag(xrhoce(ie)), dimag(cchi(ie))/xsedge, ratio*corr -! corr is the ratio N_av/N_j, responsible for difference in -! counts due to variation of wave function due to spin-orbit - 460 format(f12.6, 2x, e12.6,2x,e12.6,2x,e12.6,1x,e12.6) - write(4,465) dble(em(ie))*hart, dimag(xrhopr(ie)), dimag(cchi(ie))/xsedge, ratiop, dimag(xrhoce(ie)-xrhopr(ie))/ratio -! also write contribution to mu_0 from delocalized states defined as -! (rho-rho_proj)/ratio - 465 format(f12.6, 2x, e12.6,2x,e12.6,2x,e12.6,1x,e12.6,2x,e12.6) + endif + enddo + ! ********************* END OF BIG LOOP OVER ENERGY MESH ************** + + CALL CloseFl(TRIM(fname)//'.dat') + + ! *** optionally write ratio.dat and ratiop.dat ***** + if (ipr2.ge.3) then + ! calculate mu_0/rho_0 for XMCD normalization. + chia(1:ne) = 0 + vrcorr = 0 + vicorr = 0 + call xscorr(1, em, ne1, ne, ik0, xrhoce,xsnorm,chia,vrcorr, vicorr, cchi) + do ie = 1, ne1 + xrhoce(ie) = coni* dimag(xrhoce(ie)+cchi(ie)) + enddo + call xscorr(1, em, ne1, ne, ik0, xrhopr,xsnorm,chia,vrcorr, vicorr, cchi) + do ie = 1, ne1 + xrhopr(ie) = coni* dimag(xrhopr(ie)+cchi(ie)) + enddo + call xscorr(1, em, ne1, ne, ik0, xsec,xsnorm,chia,vrcorr, vicorr, cchi) + do ie = 1, ne1 + cchi(ie) = coni* dimag(xsec(ie)+cchi(ie)) + enddo + open(unit=3,file='ratio.dat',status='unknown', iostat=ios) + open(unit=4,file='ratiop.dat',status='unknown', iostat=ios) + ! normalize to xsec at 50 ev above edge + edg50 = emu +50.0 / hart + call terp (omega1, xsnorm, ne1, 1, edg50, xsedge) + write(3,440) xsedge, emu * hart + 440 format ('# Normalization factor:', e12.4,' Angstrom**2. Fermi level at ', f7.1, ' eV.') + write(3,450) + 450 format ('# Energy rho_0 mu_0 rho_0/mu_0 ') + write(4,440) xsedge, emu * hart + write(4,455) + 455 format ('# Energy rho_proj mu_0 rho_proj/mu_0',' mu_deloc ') + + do ie=1,ne1 + if (dimag(cchi(ie)).eq.0.d0 .and. ie.lt.ik0) then + cchi(ie)=cchi(ik0) + xrhoce(ie)=xrhoce(ik0) + xrhopr(ie)=xrhopr(ik0) + endif + ratio = dimag(xrhoce(ie)) / dimag(cchi(ie)) * xsedge + ratiop = dimag(xrhopr(ie)) / dimag(cchi(ie)) * xsedge + + write(3,460) dble(em(ie))*hart, dimag(xrhoce(ie)), dimag(cchi(ie))/xsedge, ratio*corr + ! corr is the ratio N_av/N_j, responsible for difference in + ! counts due to variation of wave function due to spin-orbit + 460 format(f12.6, 2x, e12.6,2x,e12.6,2x,e12.6,1x,e12.6) + write(4,465) dble(em(ie))*hart, dimag(xrhopr(ie)), dimag(cchi(ie))/xsedge, ratiop, dimag(xrhoce(ie)-xrhopr(ie))/ratio + ! also write contribution to mu_0 from delocalized states defined as + ! (rho-rho_proj)/ratio + 465 format(f12.6, 2x, e12.6,2x,e12.6,2x,e12.6,1x,e12.6,2x,e12.6) - enddo - close(unit=3) - close(unit=4) - endif - - deallocate(bmat) - - return - end + enddo + close(unit=3) + close(unit=4) + endif + + deallocate(bmat) + return +end diff --git a/src/XSPH/xsphsub.f90 b/src/XSPH/xsphsub.f90 index 58834a9..e658574 100644 --- a/src/XSPH/xsphsub.f90 +++ b/src/XSPH/xsphsub.f90 @@ -113,7 +113,7 @@ subroutine xsph ! NRIXS variables !KJ integer indmax - integer iri + integer iri, ixctmp ! real*8, allocatable :: xnmues(:,:) @@ -207,9 +207,15 @@ subroutine xsph e_chsh = -sh_eng(1,1) write(6,fmt='(a,f20.10)') 'e_chsh: ', e_chsh emu = e_chsh + erelax - end if + IF (ChSh_Type == 4) THEN + ! Use core level shift from user + e_chsh = custom_chsh/hart + write(6,fmt='(a,f20.10)') 'e_chsh: ', e_chsh + emu = e_chsh + erelax + ENDIF + ! Write potentials for feffFP IF(.FALSE.) THEN OPEN(unit = 33, file = 'vtot.dat', status = 'replace') @@ -543,7 +549,8 @@ subroutine xsph call fixvar (rmt(0), edens(1,0), vtot(1,0), dmag(1,0), vint, rhoint, dx, dxnew, jumprm, vjump, ri, vtotph, rhoph, dmagx) END IF call fixdsx (iph, dx, dxnew, dgc, dpc, dgcn, dpcn) - if (mod(ixc,10) .ge. 5) then + ! if (mod(ixc,10) .ge. 5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.15)) THEN ! Replace mod(ixc,10).ge.5 if (jumprm .gt. 0) jumprm = 2 call fixvar (rmt(0), edenvl(1,0), vvalgs(1,0), dmag(1,0), vint, rhoint, dx, dxnew, jumprm, vjump, ri, vvalph, rhphvl, dmagx) if (jumprm .gt. 0) jumprm = 1 @@ -573,9 +580,16 @@ subroutine xsph ! Josh - added argument iPl to control many pole self energy NPolesTmp = NPoles ! write(*,*) 'Breakdown at xsect' + ! JJK - Quick fix for finite temperature calculations. If using ground state, + ! still need to add finite temperature part at the fermi level since there is an imaginary part. + IF((electronic_temperature.GT.0.d0).AND.(ixc0.EQ.2)) THEN + ixctmp = -2 + ELSE + ixctmp = ixc0 + END IF call xsect (ipr2, dxnew, x0, ri, ne, ne1, ik0, em, edge, & ihole, emu, corr, dgcx, dpcx, jnew, & - ixc0, lreal, rmt(0), rnrm(0), xmu, vi0, iPl, NPolesTmp, Eps0, EGap, & !iPl new JK, NPoles and Eps0 new JJK 3/9/2010 + ixctmp, lreal, rmt(0), rnrm(0), xmu, vi0, iPl, NPolesTmp, Eps0, EGap, & !iPl new JK, NPoles and Eps0 new JJK 3/9/2010 gamach, vtotph, vvalph, rhoph, dmagx, rhphvl, & dgcn, dpcn, adgc(1,1,iph), adpc(1,1,iph), xsec(1,isp), & xsnorm(1,isp), rkk(1:nex,1:nq,1:kfinmax,isp), iz(0), xion(0), iunf, & !KJ 12/10 @@ -658,7 +672,8 @@ subroutine xsph ! write(*,*) ! write(*,*) - if (mod(ixc,10) .ge.5) then + ! if (mod(ixc,10) .ge.5) then + IF ((ixc.EQ.5).OR.(ixc.EQ.9).OR.(ixc.EQ.15)) THEN ! Replace mod(ixc,10).ge.5 if (jumprm .gt. 0) jumprm = 2 call fixvar (rmt(iph), edenvl(1,iph), vvalgs(1,iph), dmag(1,iph), & vint, rhoint, dx, dxnew, jumprm, vjump, ri, vvalph, rhphvl, dmagx) @@ -693,7 +708,9 @@ subroutine xsph 300 continue ! end of big loop over spins -! write main output to xsect.dat + ! JK - Shift imaginary part of self-energy by MIN(eref) since we want the fermi function broadened as well in XAS. + em(ne1+1:2*ne1) = em(ne1+1:2*ne1) - coni*DIMAG(eref(ne1+1,1)) + ! write main output to xsect.dat 340 format (e17.9, 4e13.5) if (abs(ispin).ne.1 .or. nspx.eq.1) then do ie = 1, ne