DFTB charge equilibration discrepancy

[Nir]

Hi All,

I have been benchmarking the DFTB implementation in CP2K against a
version of the DFTB code from Prof. Frauenheim's group (www.dftb.org)
and have run into some problems. FYI, the DFTB code I used force
matches to DFTB+ to within 1%. Specifically, for this problem we are
interested in systems containing C, N, O, and H. For the benchmarking
tests I'm talking about, both codes were run using charge
equilibration (scc) and Ewald sums, without dispersion, and using
identical parameter files. Similar very low SCF convergence criteria
were used for both codes.

For most systems I tested (e. g., CH4, H2O, CN, diamond), the two
codes compared very well and the forces and computed energies matched
to within ~1%. However, the forces and energies were considerably off
for systems containing C-O, N-H, and N-O bonds (e. g., CO, NH3 and
NO2, respectively). I hacked into both codes and saw that the
repulsive forces are virtually identical for both CP2K and the DFTB
code, but that for some reason they equilibrate to different charges
for the above systems. Any help or insight into solving this problem
would be greatly appreciated. I uploaded a tar file called
nir_dftb.tar of sample input and output files for a single NH3
molecule containing the following files (with hopefully self-
explanatory names):

cp2k.nh3_single.inp
cp2k.nh3_singe.out
cp2k.charges.dat
DFTB.nh3_single.out
DFTB.FRC_nh3_single.DAT
DFTB.CHR_nh3_single.dat

In case it matters, I ran cp2k.popt on a single node on a Linux-x86-64
machine compiled with the intel fortran compiler.

Thanks for your help!

Nir

[Juerg Hutter]

Hi Nir

there might be several reasons for this problem.
To get more information could you do the following tests:

- run both programs without EWALD (and PERIODIC NONE in CP2K)
in order to verify that we agree on the basics

- Converge all EWALD parameters in both codes in order to
exclude problems with different implementations

When we tested the implementation we reached agreements
with the reference code (DMON) far better than 1%

best regards

Juerg

----------------------------------------------------------
Juerg Hutter Phone : ++41 44 635 4491
Physical Chemistry Institute FAX : ++41 44 635 6838
University of Zurich E-mail: hut...@pci.uzh.ch
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland
----------------------------------------------------------

[Nir]

Hi Juerg,

> there might be several reasons for this problem.
> To get more information could you do the following tests:
>
> - run both programs without EWALD (and PERIODIC NONE in CP2K)
>    in order to verify that we agree on the basics

Thank you for your speedy reply. Unfortunately turning off the Ewald
sums in both codes (and using PERIODIC NONE in CP2K) did not produce
agreement. I also tried doing the same by using the "cluster" option
with DFTB+ without success. Is there anything else you can suggest to
try? Again, DFTB/DFTB+ and CP2K agree very well for systems such as
CO, H2O and CH4, but disagree for CO and NO. I'm using the scc input
files provided by CP2K for all three codes.

Thanks again,

Nir

[Juerg Hutter]

Hi

so it's not Ewald, good.
Without more information it is still wild guessing:

Could it be an odd number of electrons? What happens if
you are using "LSD" for NO?

Could it be that your DFTB+ run is using a l-dependent
hardness by default? CP2K has only a single hardness
value per atom implemented.
Or is DFTB+ taking another hardness value from the
input file, CP2K is using the s-value.

regards

Juerg

----------------------------------------------------------
Juerg Hutter Phone : ++41 44 635 4491
Physical Chemistry Institute FAX : ++41 44 635 6838
University of Zurich E-mail: hut...@pci.uzh.ch
Winterthurerstrasse 190
CH-8057 Zurich, Switzerland
----------------------------------------------------------

[Nir]

Hi Juerg,

> Could it be an odd number of electrons? What happens if
> you are using "LSD" for NO?

I don't think an odd number of electrons is causing the problem. For
example, if I add a -1 charge to NO, the codes still disagree.
However, contrary to what I wrote earlier, DFTB+ and CP2K do agree
for CO and NH3. They still disagree for CO2, NO2, and for systems
such as nitromethane and TATB. As a wrote earlier, I get excellent
comparisons for H2O, diamond, methane, and LN2. It's just for CO2 and
systems with NO2 that I observe a discrepancy. The energies for these
systems tend to be off by ~0.6 Ha, and the forces are off by anywhere
from 10-100%. I also tried turning off the scc charge equilibration
for both codes. The agreement improved slightly, but was still off
overall. Another thing is that CP2K tends to produce lower energies
than DFTB+. I haven't quantified this at all, but that seemed to be
what was happening when I eyeballed the results.

> Could it be that your DFTB+ run is using a l-dependent
> hardness by default? CP2K has only a single hardness
> value per atom implemented.
> Or is DFTB+ taking another hardness value from the
> input file, CP2K is using the s-value.

DFTB+ has an input flag where the code will only use the Hubbard U for
the s-shell. I assume this is the same thing as you are talking
about, above. I used that flag for the comparisons I'm talking about.

Any other thoughts on other things I could look into?

Thanks for the help,

Nir

[Juerg Hutter]

Hi

>
> I don't think an odd number of electrons is causing the problem. For
> example, if I add a -1 charge to NO, the codes still disagree.
> However, contrary to what I wrote earlier, DFTB+ and CP2K do agree
> for CO and NH3. They still disagree for CO2, NO2, and for systems
> such as nitromethane and TATB. As a wrote earlier, I get excellent
> comparisons for H2O, diamond, methane, and LN2. It's just for CO2 and
> systems with NO2 that I observe a discrepancy. The energies for these
> systems tend to be off by ~0.6 Ha, and the forces are off by anywhere
> from 10-100%. I also tried turning off the scc charge equilibration
> for both codes. The agreement improved slightly, but was still off
> overall. Another thing is that CP2K tends to produce lower energies

So the problem is not in SCC! In this case you can just print out the
Hamiltonian and Overlap matrices with both codes and look for differences.
From there it should be easy to identify the problem.

Juerg