using the code, I think I was able to differentiate the forces but had a few questions.

[mike marchywka]

I got a result but can't tell if it is reasonable if anyone cares to comment it would be appreciated. These are attempts to find force constants without perturbing the electron density - only the diagonal elements would be presumed non-zero. I started trying to code the perturbation approach in section 5.5 of DFT++ guide but ran into coding problems and still not real confident of this more primitive result yet.

As I indicated in prior thread, I was trying to find the vibrational modes of methane as a starting point to test my code. I first did this with brute force array of displaced atoms and found modes from the forces reported by jdftx. I have started hacking into the code with the help of the DFT++ Derivation guide.
I got the diagonal elements from the force constant matrix approximated by differeniating he expressions internal to JDFTX for the reported forces. This appears to be largely multiplying the force expressions by wavevector and then the gInfo.dV.
I could not actually determine this latter scale factor but it gives the right answer and is at least plausible. In the Vnl code, I ended up throwing daggers around until I seemed to get the right result as my matrix differentiation is not great but ended up with this,

diagMatrix cdiag = diag(MVdagC * Fq * dagger(dV2[k]^Cq)+(Mdagger(dV[k]^Cq)Fq*dagger(dV[k]^Cq)));

      (curves)[atom][k] += gInfo.dV*2.0*Cq.qnum->weight * cdiag[p+atom*nProj];

The ewald component is apparently easy to do analyitically but a bit tedious.

These are xx,yy,zz components of force constants for the 4 hydorgrns followed by the
carbon in various pieces,
MJM ewald CURVE 5.34505e-05 1.2585e-06 -0.00588062
MJM ewald CURVE -0.00556138 1.74345e-06 0.0019106
MJM ewald CURVE 0.00275628 0.00480066 0.00198636
MJM ewald CURVE 0.00275256 -0.0048036 0.00198423
MJM ewald CURVE -9.23539e-07 -5.45571e-08 -5.73132e-07

MJM local CURVE 0.117898 0.117981 0.0323116
MJM local CURVE 0.0414687 0.118016 0.1089
MJM local CURVE 0.0987764 0.0605658 0.107689
MJM local CURVE 0.0988283 0.0604951 0.107711
MJM local CURVE 0.350265 0.353225 0.349813

MJM MIKE curveNL h 0.00377684 0.00379081 0.00299545 1
MJM MIKE curveNL h 0.00267793 0.00404087 0.00384243 1
MJM MIKE curveNL h 0.00369407 0.00331586 0.00412333 1
MJM MIKE curveNL h 0.00346529 0.00309643 0.00368302 1
MJM MIKE curveNL c 0.256989 0.2553 0.257415 1

MJM NET TOTAL CURVES h 0.121729 0.121773 0.0294264 1
MJM NET TOTAL CURVES h 0.0385853 0.122059 0.114653 1
MJM NET TOTAL CURVES h 0.105227 0.0686824 0.113799 1
MJM NET TOTAL CURVES h 0.105046 0.0587879 0.113378 1
MJM NET TOTAL CURVES c 0.607252 0.608524 0.607227 1

These compare somewhat with the diagonal elements of he brute force technique remembering that the coordinates here are different with each ""x" being the direction to nearest neighbor. The carbon seems a bit high with hydrogens a bit low but hopefully when I add the perturbations to get delta-n that will work out.

cat methane_phonon_gga_nat2.tcl | ./tcl_tool | grep FCONST | awk 'BEGIN{i=3;} {print $i; ++i}'
-0.342744
-0.0464313
-0.0463531

-0.342823
-0.0464951
-0.0467174

-0.342959
-0.0464457
-0.0467845

-0.342923
-0.0464237
-0.0467629

-0.542804
-0.518814
-0.528373

Thanks again.

[Ravishankar Sundararaman]

Hey Mike,

It's really nice to see that you've started implementing density-functional perturbation theory into the code! We have always wanted to have that functionality, but haven't found the time so far to implement it.

Would you be interested in formalizing this work and committing it as a feature of the code? I can probably help you plan out the command structure etc. for actually implementing something like this in a non-hacky way.

Thanks!
Shankar

[mike marchywka]

Well, thanks I'd like that but I'd like to get some idea I know what I'm doing first :) Still not sure I got any of those derivatives right except that when you add them up they seem about right based on the brute force approach. I was going to try to work through the perturbation theory next but keep hitting coding issues- I may be approaching this incorrectly. I tried to make a copy of "Everything" and then thought I could modify it but I have to go find all the things that don't copy and it is turning into a huge mess.

But yes all this stuff I've kluged in in places where I though it would be least obtrusive so I don't have to back out my only working version. It looks like from the DFT++ docs circa 2000 you have put a lot of work into various versions that are designed for modern architectures already, I assumed from your pubs you had the perturbation and related code internally and just had not released it yet.

Thanks, happy to contribute.

[mike marchywka]

FWIW, I haven't forgotten about this but it has taken a lot more time than I expected. I had to back up a bit and work with Stark effect in CO instead of going right to vibrations and still looking for more simple systems for testing. I generated a constant electric field for input to JDFTX as an external potential and then found the energy shifts component by component. I could then calculate the same components using the perturbation code and play with it until they matched. I finally was able to get coarse agreement, about 10 percent ( see below ) and think there is some hope I can eventually get useful code. I had to get more papers and plow through basics and keep retreating a bit but as they say "I'm learning a lot." LOL

This is about where I am now, the energy components are named after the things
in Energies more or less, column 2 is the name, col 5 is the shift from JDFTX calculations with and without external field and col 6 is the ratio of the perturbation result to column5. ( these are taken along a plateau
after many iterations but it still does not "converge" as I have
the energy limit set way down ),

[1] "ei.eh 3 3 -0.00118182 1.00468768509587"
[1] "ei.eloc 8 3 0.00136951 1.0826792064315"
[1] "ei.enl 9 4 6.82395e-05 1.07171213153672"
[1] "ei.exce 18 3 0.000132755 1.03205905615608"
[1] "ei.kecombo 24 4 -0.000387906 0.974916603507035"

kecombo is the kinetic energy, a bit below target with enl and eloc being a bit high. I think I can make this better but getting this far took some time.

I've been making graphs of these values as the minimizer iterates and
just trying to understand it right now. I will probably try to get the
energy levels next and see how those look and then try to fix some of
the other kluges. This seems to have the potential to be very fast
if I can move all the zero-order ( constant ) calculations out of the
iteration loop and fix other basic problems.

I do have code to write DataRptr arrays to an xplor file which apparently works in Pymol( I could visiualize electron density ) but the coordinate transforms are not complete.

Thanks again, definitely intresting.

[Deniz Gunceler]

That looks really promising Mike.

How did you generate the constant electric field? With a sawtooth potential? How big are the electric fields you are using?

[mike marchywka]

Thanks, I'll have to go back and check that as I finally set the strength to get level shifts I thought were small. But yes it should just be an external potential file generated from code that I also include with my jdftx build
and then apply internally in the perturbation code.

I probably should have started with LDA too but I'm hoping to figure this out with the TPSS :)