still have a problem with TPSS but could not post plots,
Joint Density Functional Theory
Brought to you by:
shankar1729
Menu
▾
▴
I was still trying to look at Exc(r) for H comparing PBE and TPSS with JDFTX AE
and an analytical calculation. I wanted to post these plots but could not figure out
how to do it on the forum. There was no obvious way to attach them.
Have you tried TPSS with just a simple AE atom? It is possible I
messed up the code...
Thanks.
Testing attachments ... There seems to be an attach button right next to post. Or is that only for the admins?
-Shankar
it seems to be an option for reply, still does not show up for
starting a topic. It seems to let me pick and image but not sure if it is uploading.
ok, I guess for reply it works but still can't find it for starting a topic.
Maybe I need to update the browser...
Anyway, it looks like the most relevant one posted. The Exc(r) I dumped for
TPSS from JDFT doing what should be an AE H atom seems to oscillate
but NB all prior caveats in first post. This graph should have 4 curves, each showing Exc(r) from proton energy in Ha with distance scaled. The green curve was JDFTX using PBE, red/blue were test code for analytic hydrogen with PBE and TPSS done using LIBXC called from "R". The oscillating curve is the TPSS result. 120 Ha cutoff.
What you are observing seems to be Nyquist frequency oscillations. If you apply a smoothing kernel that averages out the oscillations, it looks to me like the results would line up.
The magic of the plane-wave basis is that the final energy can be correct even though the potential is not resolvable on the grid - note that the energy is related to the integral of this against the electron density which is a smooth band-width limited function.
I am curious though - what is the quantity that you're plotting i.e. what is its definition / expression and what is the JDFTx code snippet that generates it?
I wrote some code to dump grids as either ssv or xplor format but the coord conversion was not complete. And yes I am concerned I added bugs after a year or so of hacking :) Anyway, this is the relevant code,
e.exCorr(get_nXC(1), &Vxc2, needKE, &tau, &Vtau2);
mjm_dump_ssv(StrTy("vxczed.ssv"), StrTy("my"), e.gInfo, Vxc2[0]);
All the other stuff came from a TCL tool to create hydrogen rho,sigma, and tau
from closed form result and then using "R" interface to libxc to get the
X and C energies added together to make these plots.
( not debugged but results seemed to line up with jdftx output so far,
the only reason I'm posting now is to get some idea of reaosnablenes of oscillations in TPSS ).
I thought the "R" tool would be nice for comparing XC values from many
different algorithms.
Thanks.
OK, the issue is most likely the resolvability of the AE tau on the Fourier grid. The gradient of the cusp'd exponential at the origin is particularly problematic. On a radial grid, the r=0 doesn't contribute elsewhere, but with the Fourier basis the non-differentiability at r=0 pollutes the derivative everywhere else, and the oscillations die out on a length scale inversely proportional to Gmax, the maximum reciprocal lattice vector in the basis. This is an issue for both grad n entering GGA and tau entering mGGA, but note that the Gmax for wavefunctions is typically 1/2 that for densities in plane wave electronic codes.
What I would recommend is using a pseudopotential with a very small core radius so that the cusp is eliminated and replaced by a little analytic piece at some distance from the origin, but then only looking at the results outside that core radius. I have attached an fhi file (made using opium) with a core radius of 0.25 bohrs, which is just soft enough to be resolvable at the 0.1 mEh energy scale at 120 Eh cutoff. Note that outside this core radius the potential is exactly 1/r. I can attach an even smaller core if you're interested, but you'll need to ramp the cutoff even higher.
ok, thanks. Is this gibbs phenomenon or something else? I should know the fourier problems by now but I will need to research. It occurred to me I did not run PBE with AE, I'm trying that now. Also note that the libxc "TPSS" had much reduced oscillations probably due to interface bug in tau convention.
Yep, basically Gibbs phenomenon.
Also, fwiw, here is the earlier GGA PP that had been working fine, IIRC it came from QE or Abinit site. I thought this was a bug at first but seems to vary with parameters we are discussing.
!
So anyway then I take the ssv xc file and graph against these things
in "R", mjmlibxc is my interface to libxc that as with all the other pieces is still being tested. ( none of this had been checked yet).
cat testH.R
dyn.load("foo/src/mjmlibxc.so")
library("mjmlibxc",lib.loc="./mjmlibxc.Rcheck")
df<-read.table("testH.txt",sep=" ",header=T)
str(df)
.Call("find_xc",df)
a<-.Call("libxcdo",df)
str(df)
str(a)
df
dfj<-read.table("vchpbe.txt",sep=" ",header=F)
dfjtpss<-read.table("vch.txt",sep=" ",header=F)
dfjtpssae<-read.table("vchtpssae.txt",sep=" ",header=F)
dfjtpsslibxcae<-read.table("vchtpsslibxcae.txt",sep=" ",header=F)
df2= df
df2$functional<- df2$functional+29
str(df2)
a<-.Call("libxcdo",df2)
df3= df
df3$functional<- df2$functional0+101
a<-.Call("libxcdo",df3)
df4= df
df4$functional<- df2$functional0+130
a<-.Call("libxcdo",df4)
...
source("testH.R")
plot(dfj$V1,col="green",xlim=c(0,100))
points(df$radius10+2,df3$drho+df4$drho,col="blue")
points(df$radius10+2,df$drho+df2$drho,col="red")
points(dfjtpss$V1,col="black")
?legend
?legend
legend("bottomright",c("PBE-analytic","TPSS-analytic","PBE-JDFTX-PBEAtom","TPSS-JDFTX-PBEAtom"),col=c("blue","red","green","black"))
?legend
legend("bottomright",c("PBE-analytic","TPSS-analytic","PBE-JDFTX-PBEAtom","TPSS-JDFTX-PBEAtom"),col=c("blue","red","green","black"),pch=1)
png("xcompare.png")
plot(dfj$V1,col="green",xlim=c(0,100))
points(df$radius*10+2,df3$drho+df4$drho,col="blue")
ok, that is better but than AE but did not create the apparent chaos from the PBE PP ( I guess knowing what you are doing helps LOL). It is interesting however that the PBE AE sets right on top of the PBE-PBE and they both match more or less my analytical H with PBE or TPSS ( similar precision of these two matching each other ). The lengend is a bit cryptic but the + is your latest PP with TPSS in JDFTX, the X is the AE TPSS and the triangle and green circles are PBE with AE and old PBE PP respectively. ( closed form H not shown).
So right now the oscillations appear to go with the JDFTX TPSS and can be minimzed with proper PP selection or made non-obvious with careless selection. The oscillations went to very low with the libxc TPSS which is probably like your TPSS except for a factor of 2 error in tau. While the original objective was to debug my new code, this should all be helpful in understanding the TPSS tau part too. Thanks.
!
If I make this change to ExCorr_internal_mGGA.h the problem goes away.
the real odd things was vtau,
head runs/h1/zed_dump_headvtau.ssv
0 0 0 +0.000000000000000000e+00
1 0 0 +0.000000000000000000e+00
2 0 0 +6.473175360293148550e-02
3 0 0 +0.000000000000000000e+00
4 0 0 +4.534857574820924664e-02
5 0 0 +0.000000000000000000e+00
6 0 0 +3.760358776751759197e-02
7 0 0 +0.000000000000000000e+00
8 0 0 +3.059191417873884794e-02
9 0 0 +0.000000000000000000e+00
One patch, which I have not entirely thought through but does eliminate the symptom is this. The cutoff value seems a bit high but also I am not sure you want to skip this if tau is out of range or not. I have a lot of plots and equations in a latex file, it is actually interesting but of no use now LOL.
ls -al jdftx/electronic/ExCorr_internal_mGGA.h recent/electronic/ExCorr_internal_mGGA.h
-rwxrwxrwx 1 marchywka marchywka 21521 2013-11-17 13:02 jdftx/electronic/ExCorr_internal_mGGA.h
-rw-rw-r-- 1 marchywka marchywka 20928 2013-03-05 09:23 recent/electronic/ExCorr_internal_mGGA.h
marchywka@marchywka-Latitude-E6510:~/d/jdft$
< // mjm this seems to be 1e-16 which is large...
< // if(tauTot<nCutoff) return;
< double mjm_zed=nTot*tauTot;
< double z_sigma = mjm_zed?(0.125/(mjm_zed)):0;
< //double z_sigma = 0.125/(nTot * tauTot);
I have attached a plot showing that this version agrees with PBE and analystical H to much better extent although HH bond length has not changed much. I am working on a low pass filter which accomplished similar results. This plot also shows PBE vs TPSS and you can still see some higher frequency noise of unknown significance. I suspect it does not matter in results but may be interesting for convergence of learning tool. Thanks.
I got a little lost in that last message - could you point out exactly what the code change you made was that fixed the oscillation? :)
Yes, I'd be happy to but I just wanted to make another attempt to verify it is not something really stupid I did while playing with the code :) I had worked through your TPSS and written my own and they seemed to agree on a few data points for things like Vrho but I can not just easily drop in mine for yours in this setting. The symptom is still a bit distant from
this change but I'm trying to piece it all together- squaring or dividing band limited signals is a bit interesting in fourier space. If you just decde
to set some points to zero, multiply by square wave, I guess that could produce some ringing etc.
Also note in the prior plot that if create Vxc(TPSS)/Vxc(PBE) the jdftx and analytical results seem to be offset, there are probably a bunch of things like offsets I have not fixed properly.
This code in question seems to use an arbitrary cutoff and avoids the exc evaluation based only on the value of tau.
Tau of course has a denomintator that can be zero but should have a well defined value AFAICT for any physical situation.
This is the code I have in ExCorr_internal_mGGA.h circa line 190
( but likely some chars will be taken as formatting commands.... )
// mjm this seems to be 1e-16 which is large...
// if(tauTot<ncutoff) return;="" double="" mjm_zed="nTot*tauTot;" z_sigma="mjm_zed?(0.125/(mjm_zed)):0;" *="" tautot);="" z="z_sigma" sigmatot;="" bool="" zoffrange="false;" mjm="" if="" (="" !mjm_zed)="" is="" now="" nan,="" may="" be="" a="" problem?="" if(z="">1.) { z=1.; zOffRange = true; }</ncutoff)>
I see, so basically you replaced:
with:
(and appropriately changed the next few lines to not overwrite that behavior). That sounds like a decent idea.
I think by the time I was done I only checked for zero denominator, I don't think you ever want to return at this point but instead as a last resort set z to some limit. There are a mix of observations and verbose codes as I'm still not sure on best way.
But yeah I can't tell from my prior post what that was upposed to be, 1/2 the code got taken as formatting commands....LOL.
See if this helps,
AFAICT if everything is ok there is no need for a cutoff but how this translated into the symnptoms is another matter. I was debating about trying to pass "z" instead of making it up from pieces. Quite interesting actually in PW basis similar to signal processing now :).
Anyway, after some more analysis as you probably guessed it is unlikely that was a significant factor, the oscillations would come and go and also occur with libxc intermittently.
However, it does appear that if you calculate Ztpss directly from psi and derivatives, without finding "n" first, Ztpss is almost assured to stay within range and the nyquist stuff reduces significantly at least in 1 electron case I have not closely examined limit as denomiator goes to zero but in any case where it is not exactly zero the sum should stay in range as it reduces to something like c+ Re(x)/|x|
I have some simulations with closed form hydrogen and various numerical approaches and this seems to work - z would stay in bounds. If I did the math right( and I have not checked it) , cases where electron density of one state overlaps with ked of another state should reduce the risk of out-of-bounds ( something like c+ Re(x)/(|x|+|a|). also ran a few cases with finite difference and lanczos derivatives.
Not sure any of this is of practical importance but I wanted to understand
how the various energy terms act and the kind of oscillation originally
described is quite confusing :)
Thanks