Does anyone have any tips for doing geometry optimizations with constrained DFT? It seems that trying to constrain some molecules to oxidized states on a carbon surface model and relaxing geometry is difficult. The total energy seems to not track smoothly to a lower energy and frequently goes up. Is this just a matter of needing higher thresholds all the way around than for a non constrained optimization or is there a simpler fix? I know this post is a little vague but any tips on geometry optimizations with constrained DFT are greatly appreciated.
This may be a bug in the constrained DFT code. We have reported it as a ticket (#2347). For a test job, analytic gradients do not agree with finite-difference results (IDERIV=0).
Just an update to say that we are looking into this. Errors between analytic and finite-diff gradient are not huge (~ 1e-4 in a.u.), but larger than I’m accustomed to seeing, and we see the behavior that you describe, i.e., energy does not go down consistently in some cases. Bottom line is that I think you can continue using the gradients for now, but you are not imagining things or doing something wrong, so far as I can tell.
I’m no longer 100% convinced this is a bug (though still investigating). As a consequence of looking into this, however, Evgeny @ Q-Chem discovered 3 undocumented job control variables, which we will add to the manual for the 5.3.2 bug-fix release but I can tell you here. They are:
CDFT_THRESH ==> threshold for convergence of CDFT weights
CDFT_MAXITER ==> max number of iterations for converging the CDFT weight
CDFT_PRINT ==> enable printing of CDFT iterations
The most important one in the present context is the first one, which defaults to 5 and that was giving erratic behavior. A better choice may be to set it to 8 for geometry optimizations, which is the default value for SCF_CONVERGENCE in the case of a geometry optimization. For the example that we were looking at, this makes things look a lot better. Please try setting CDFT_THRESH = 8 in your $rem section and (if you don’t mind) reporting back whether or not this solves your problem.
I tried using the higher threshold but still the total energy climbs quite quickly for my test case. I will try another test case and see what that does.
You might find it useful that it seems that everything works okay if doing a MD run but it is when doing an OPT run the energy keeps rising.
Is your molecule top secret or do you mind posting an input file that I could submit as a bug ticket? That way the Q-Chem scientific staff (of which I am not) could take a look. And probably I would look too.
Do you know what would cause the Becke populations not to add up to the constrained value?
Poor quality grid. Note there is an XC_GRID for c-DFT even if you use Hartree-Fock theory. Try something like XC_GRID = 000099000590, i.e., (99,590). That is “high quality” so I would hope is good enough. The “kitchen sink” grid that we use as a stand-in for grid limit is (250,974).