Geometry optimization with constrained DFT

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.

-John

Hi John,

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.

thanks,
Jason

Hi John,

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.

thanks,
Jason

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.

Hi John,

Do you know what would cause the Becke populations not to add up to the constrained value?

thanks,
Jason

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).

Hi John,

From my testing it seems that the issues are related to the use of becke_shift universal_density and bragg_slater. It seems I can do relaxations without them but using either universal_density or bragg_slater the energy goes up really fast and I get strange geometry relaxation. On the other hand it seems becke_shift universal_density(or bragg_slater possibly) is needed to get reasonable binding energies at a fixed geometry. So now I don’t really know what to do and my project is a bit stuck.

I have a few possible suggestions.

First, I do not think that the quality of the CDFT gradients should be sensitive to the particular way in which the atom-centered grid is weighted. The constrained value will be different, but the differences in the XC quadrature grid between BECKE_SHIFT = NONE, BRAGG_SLATER, or UNIVERSAL_DENSITY occur at the level of the DFT quadrature and don’t have anything to do with the CDFT code.

Second, ensure that you turn integral symmetry off for CDFT jobs (SYMMETRY = FALSE), as the constraints allow the electron density to violate point-group symmetries. This is not the default according to ticket #1916. Could be that in one case, the violations of point-group symmetry are small or negligible, but in the other case (with the re-weighted quadrature) they are not.

Lastly, it is possible as well that CDFT constraints generate very unphysical electron distributions that could lead to erratic behavior. For a softer constraint, it might be advisable to try applying a charge and a spin constraint within the CDFT-CI framework, as CDFT-CI corrects for the constraints based on intuition via the construction of promolecule densities. CDFT simply takes the integer constraint and runs with it, which can be satisfactory in cases where the charges are spatially distant from one another, but in cases where the electron clouds of donor and acceptor overlap significantly the softer CDFT-CI constraint is likely the better option. See sec. 2.5 of Kaduk, Kowalczyk, and Van Voorhis Chem. Rev. 2012 (DOI: 10.1021/cr200148b).

So I have tried increasing the grid to XC_GRID 3 and setting SYMMETRY FALSE but I still can not find a case were doing AIMD with the CDFT and BECKE_SHIFT UNIVERSAL_DENSITY does not have a nearly continuously increasing internal energy and erratic temperature. I am even trying very small time steps of 2 or 1 atomic units. I have not tried the CDFT-CI because I am still not completely clear on the theory or how I would implement it.

Do you happen to know of an example of using CDFT with geometry optimization or AIMD that does not exhibit erratic behavior?

thanks,
Jason

I posted your query on our internal developer page and got the following example (from Yuezhi Mao). Ethylene dimer with a charge on one monomer. Energy goes down, step-by-step in the optimization.

$molecule
1 2
C 0.00000 0.00000 1.00000
C 0.00000 1.33200 1.00000
H -0.92878 -0.57430 1.00000
H 0.92878 -0.57430 1.00000
H 0.92878 1.90630 1.00000
H -0.92878 1.90630 1.00000
C 0.00000 -0.00000 5.00000
C -0.00000 1.33200 5.00000
H 0.92878 -0.57430 5.00000
H -0.92878 -0.57430 5.00000
H -0.92878 1.90630 5.00000
H 0.92878 1.90630 5.00000
$end

$rem
JOBTYPE OPT
METHOD wB97X-D
BASIS 6-31G(D)
UNRESTRICTED TRUE
MAX_SCF_CYCLES 100
SCF_ALGORITHM DIIS
SCF_CONVERGENCE 8
THRESH 14
SYMMETRY FALSE
SYM_IGNORE TRUE
SCF_FINAL_PRINT 1
CDFT TRUE
CDFT_THRESH 7
MEM_TOTAL 8000
MEM_STATIC 2000
SCF_GUESS SAD
$end

$cdft
-1.0
1.0 1 6
1.0
1.0 1 6 s
$end

I also agree that doing the geometry optimization for a single diabatic state under the CDFT-CI framework might help since that uses a promolecule scheme to determine the “softened” constraint values rather than using the integer constraint values specified by the user. Based on my previous testing, this should work correctly for pure GGA (e.g. BLYP) or global hybrid GGA (e.g. B3LYP) functionals while it doesn’t correctly support meta-GGA and RSH functionals yet.

Here is an input for the ethylene dimer example above, which is quite similar (the functional has been changed to B3LYP since wB97X-D is buggy for CDFT-CI):

$molecule
1 2
C 0.00000 0.00000 1.00000
C 0.00000 1.33200 1.00000
H -0.92878 -0.57430 1.00000
H 0.92878 -0.57430 1.00000
H 0.92878 1.90630 1.00000
H -0.92878 1.90630 1.00000
C 0.00000 -0.00000 5.00000
C -0.00000 1.33200 5.00000
H 0.92878 -0.57430 5.00000
H -0.92878 -0.57430 5.00000
H -0.92878 1.90630 5.00000
H 0.92878 1.90630 5.00000
$end

$rem
JOBTYPE OPT
METHOD B3LYP
DFT_D D3_ZERO
BASIS 6-31G(D)
UNRESTRICTED TRUE
MAX_SCF_CYCLES 100
SCF_ALGORITHM DIIS
SCF_CONVERGENCE 8
THRESH 14
SYMMETRY FALSE
SYM_IGNORE TRUE
SCF_FINAL_PRINT 1
CDFTCI TRUE
CDFT_THRESH 7
MEM_TOTAL 8000
MEM_STATIC 2000
SCF_GUESS SAD
$end

$cdft
-1.0
1.0 1 6
1.0
1.0 1 6 s
$end