Home Page for Computational Physics 2006
(Cornell Phys 480, Phys 680 and Astro 690)
(Cornell Phys 480, Phys 680 and Astro 690)
Ewald: -0.0885440851224253 H, Electronic: -0.356725655980680 H.
Ewald: 0.3135897639081198 H, Electronic: -1.4495422745929 H
Ewald: -31.9038 H, Electronic: 1.0561 H
I've computed the exchange correllation energy as per equation (17) in
http://www.muchomas.lassp.cornell.edu/P480/Psets/ps8/node34.html
I think there is a mistake, exc(I n) should be simply exc(n), as it is in equation (18). That corrected my problem.
Yes, that was a typo! It should always be exc(n)! (n is already
a set of sample values, not expansion coefficients)
-Prof. Arias
I found one minor typo in node10.html, when you say to debug "after having just run your ``sch.m'' " but you mean dft.m.
Also, in node11.html, it might be better to use different colors for the different lines, instead of having 4 of them be blue.
...
Finally, for the hydrogen molecule, I calculated the total energy for
each of the separations you gave, but then you say to see the files
here for an analysis:
/home/Faculty/muchomas/Minicourse/H2
But there is no "Minicourse" directory ... are these files somewhere else?
Thanks for letting me know about these things!
For the first item, I'll update the problem set. For the last item, I'd deleted the directory by mistake. It is now restored.
For the second item, everyone should replace the current semilogy(...) statement in "dft.m" with
semilogy([1:length(Elm)],Elm-43.3371147782040,'-;lm;', ... [1:length(Epclm)],Epclm-43.3371147782040,'-;pclm;', ... [1:length(Ecg1)],Ecg1-43.3371147782040,'-;pccg-FR;', ... [1:length(Ecg2)],Ecg2-43.3371147782040,'-;pccg-PR;', ... [1:length(Ecg3)],Ecg3-43.3371147782040,'-;pccg-HS;');
-Prof. Arias
This movie is from a web based quantum mechanics course which has a variety of other interesting animations.
Hi Prof. Tomas Arias:
OK. My first question is about the behavior of the wave function. I believe that what we got for psi[0] is some wave function (I'm taking the Schrodinger's equation). If that is correct, why it changed so much over small energy difference? My psi[0] goes basically the same trend as you described on the blackboard this morning, but for the -0.4999999999999x, it goes up to e+60 order (N=2000, the N=200 case is the same with what you wrote on the blackboard), although it did flip sign between "3 and 4". Is it normal? I could not figure a reasonable explanation for that. Is it because that wave function goes to infinity at point r=0? Or is there any real physical reasons?
The second one is about end points. We deal with the r[N]=dr[N]=0, then everything is ok. But when I added r[0]=dr[0]=0, the values of psi[0] changed. Still, I don't think it will affect the result if you stare at the change of variable part. Or, I got something wrong in my routines?
Two very good questions!
My apologies for not answering earlier. My email inbox sometimes gets
quite backed up. So that you know, usually, I reserve Monday mornings for
catching up.
For the first question, the large values of psi are normal. Also, you
are understanding things correctly. Bascially, the large values of
psi come
from the fact that psi->0 at infinity. psi is exponentially decaying
out at
infinity, so when we integrate backwards, it is exponentially growing.
The
numbers are much larger when we have more points basically, because
our grid
is reaching out further towards infinity. All of these large numbers,
of
course, will go away in the end when we normalize the wave functions,
the
process we talked about in lecture.
For your second question, I am pleased that you are thinking carefully
about
the end points. To handle the infinities there, one should always
take the
appropriate limits and make sure that the choices you make for r[] and
dr[]
at the end points make the limits work appropriately. For u->1
(r->infinity), the wave functions are usually exponentially decaying and
that decay kills off any algebraic growth in terms like r and dr. So,
the
limit of the final products taken as u->1 will be zero, and so setting
r[N]=dr[N]=0 to avoid Nan's is fine.
For u->0 (r->0), things are a little more subtle. Analytically, at
small r,
we have that psi is proportional to r, but also we hake that
T=(E-V(r)) goes
like 1/r, and finally dr goes like 1 (because r=u/(1-u)=u+u^2+u^3+...).
Thus, the product T*psi*dr which appears in the code goes to some
finite limit as r->0, and so setting dr[0]=0 (which will treate the
limit
effectively as zero) isn't quite correct. A nice way to fix this
uncertainty, is to try the map r=u^2/(1-u). Because then, near u->0, we
have r=u^2, dr=2*u=2*sqrt(r). Thus, T*psi*dr will go like
(1/r)*(r)*(sqrt(r))->0, and then it is fair to set dr[0]=0 so that the
value
which the computer evaluates at r=0 is the correct value of zero.
Cheers,
-Prof. Arias
Solutions: [Postscript, PDF].
Solutions: [Postscript, PDF].
Solutions: [Postscript, PDF].
Hint: If your return value from getg() is acting wierd, don't forget to #include "physics.h" at the top of your main programs (and include all of the function prototypes from physics.c into physics.h)!!! Several students have had this problem and this fixed it, so it may help you too!
Solutions: [Postscript, PDF].