search

shankar1729 / jdftx

JDFTx: software for joint density functional theory
http://jdftx.org
82 stars 54 forks source link

Regarding formation energy differences between JDFTx and QE #40

Closed ttgmichael closed 6 years ago

ttgmichael commented 6 years ago

Hey @shankar1729 ,

I've been comparing jDFTx energies with the ones I get from QE pwscf for the formation energy of CO and adsorption formation energies of CO on Pt(111) and Cu(111) and I get some deviation in numbers.

I tried switching between the GBRV pseudopotentials from the installation and the one at my cluster and they get the same results.

The XC functional I'm using is PBE. And I kept the elec/ rho energy cutoff to be 500eV / 5000eV. (in hartrees, 18.3746545014 183.746545014)

What I get with jDFTx is similar to the values in the example in the documentation;

BFGS: 0 15:12:29 -591.216686 8.6057 BFGS: 1 15:15:01 -591.705210 2.9871 BFGS: 2 15:17:30 -591.735280 1.8346 BFGS: 3 15:19:59 -591.748251 0.1889 BFGS: 4 15:22:29 -591.748399 0.0123

What I get with a very similar setting in QE is:

BFGSLineSearch: 0[ 0] 15:19:04 -593.902306 8.6049 BFGSLineSearch: 1[ 2] 15:22:24 -594.419610 1.7779 BFGSLineSearch: 2[ 4] 15:25:12 -594.434013 0.0597 BFGSLineSearch: 3[ 5] 15:26:24 -594.434028 0.0044

The E components in both codes seem to only agree on E_xc (I'm not sure how to connect the others from jDFTx to pwscf.) It appears Eewald disagrees between the two codes. I'm just multiplying jDFTx hartree energies by 2 to get the rydberg values from pwscf.

jDFTx: 

 Eewald =       11.1193940300900049
 EH =       28.7476069859467707
 Eloc =      -69.8132696518777323
 Enl =        3.1415972827461465
 Exc =       -5.2632552603605181
 Exc_core =        0.0987190434473193
 KE =       10.2228625447723616
 Etot =      -21.7463450252356481
 TS =        0.0000000000000001
 F =      -21.7463450252356481

Pwscf: The total energy is the sum of the following terms:

 one-electron contribution =     -98.23791950 Ry
 hartree contribution      =      50.17414508 Ry
 xc contribution           =     -10.52677357 Ry
 ewald contribution        =      14.90045088 Ry
 smearing contrib. (-TS)   =      -0.00000000 Ry

Could you point me to a potential area I may have missed in setting up the CO calculation?

my jdftx input looks like this:

lattice \ 22.676711 0.000000 0.000000 \ 0.000000 22.676711 0.000000 \ 0.000000 0.000000 25.038869

elec-cutoff 18.3746545014 183.746545014 spintype no-spin elec-smearing Cold 0.00459366362534 elec-ex-corr gga-x-pbe gga-c-pbe davidson-band-ratio 1.1 kpoint 0 0 0 1 kpoint-folding 1 1 1

coords-type cartesian ion C 11.338356 11.338356 11.338356 1 ion O 11.338356 11.338356 13.700513 1

ion-species ..//jdftx-1.4.2/psp/gbrv/c.upf ion-species ..//jdftx-1.4.2/psp/gbrv/o.upf

coulomb-interaction isolated

pw.inp looks like: &CONTROL calculation='relax', prefix='calc', pseudo_dir='/home/users/ttgmichael/psp/gbrv1.5pbe', outdir='.', tprnfor=.true., / &SYSTEM ibrav=0, celldm(1)=1.8897261245650618d0, nat=2, ntyp=2, ecutwfc=36.7493237909d0, ecutrho=367.493237909d0, nbnd=15, occupations='smearing', smearing='fd', degauss=0.00734986475817d0, input_dft='PBE', / &ELECTRONS diagonalization='david', conv_thr=7.34986475817e-08, mixing_beta=0.7d0, electron_maxstep=100, conv_thr=7.34986475817e-08, / &IONS ion_dynamics='ase3', / CELL_PARAMETERS 20.000000000000000d0 0.000000000000000d0 0.000000000000000d0 0.000000000000000d0 21.250000000000000d0 0.000000000000000d0 0.000000000000000d0 0.000000000000000d0 20.000000000000000d0 ATOMIC_SPECIES C1 12.011d0 C.UPF O1 15.9994d0 O.UPF ATOMIC_POSITIONS {crystal} C1 0.500000000000000d0 0.470588235294118d0 0.500000000000000d0 O1 0.500000000000000d0 0.529411764705882d0 0.500000000000000d0 K_POINTS automatic 1 1 1 0 0 0

shankar1729 commented 6 years ago

This is a convention difference between JDFTx and QE on how the internal (partial) core energies are handled. In JDFTx, the partial core energy (Exc_core) is subtracted from the total energy. Since this is a constant per atom, this does not affect energy differences. (We may switch this to make it more consistent with other codes in the next major version change (i.e. in 2.x).) You can see that the forces are very similar between the two codes.

The remaining components are probably disagreeing because of Coulomb truncation in JDFTx vs not in QE. For systems without net charges and dipoles, the total energy will be similar between these cases, but the breakup into Hartee, Ewald etc. will differ.

Best, Shankar

ttgmichael commented 6 years ago

Hey Shankar,

Thanks for the fast reply! I see there are different conventions between jDFTx and QE at the moment. Yes the forces are very similar (and so are the C-O bond distances across codes).

It seems like that covered the actual discrepancies in the output. The rest of this comment is just exploring the other cases metntioned: if coulomb truncation (the coulomb-interaction setting) in JDFTx play a role or if dipoles in the system play a role.

I set the cells of CO to be periodic, so that the input looks like:

lattice \ 22.676711 0.000000 0.000000 \ 0.000000 22.676711 0.000000 \ 0.000000 0.000000 25.038869

elec-cutoff 18.3746545014 183.746545014 spintype no-spin elec-smearing Cold 0.00459366362534 elec-ex-corr gga-x-pbe gga-c-pbe davidson-band-ratio 1.1 kpoint 0 0 0 1 kpoint-folding 1 1 1

coords-type cartesian ion C 11.338356 11.338356 11.338356 1 ion O 11.338356 11.338356 13.700513 1

ion-species /home/groups/suncat/mtt013/group-resources2/jdftx-1.4.2/psp/gbrv/c.upf ion-species /home/groups/suncat/mtt013/group-resources2/jdftx-1.4.2/psp/gbrv/o.upf

coulomb-interaction periodic

dump-name temp.$VAR dump End Forces Ecomponents

But the energies seem to be near identical against the coulomb-interaction isolated setting.

BFGS: 0 11:42:42 -591.216635 8.6052 BFGS: 1 11:46:08 -591.705154 2.9875 BFGS: 2 11:49:38 -591.735228 1.8353 BFGS: 3 11:53:04 -591.748208 0.1890

I also tried a system that doesn't have any dipoles (CH4) for both QE and jdftx with the same settings (PBE, GBRV v1.5, periodic cells, and 500eV cutoffs). The energies from the ase output file is:

jdftx: BFGS: 0 11:41:56 -220.131625 0.0042

pwscf: BFGSLineSearch: 0[ 0] 11:02:17 -221.041677 0.3222 BFGSLineSearch: 1[ 2] 11:03:45 -221.047685 0.0002

the E contributions between both codes for the CH4 case is:

jdftx:

 Eewald =        5.8709843421693897
 EH =       11.5705730020518800
 Eloc =      -28.1980744672991399
 Enl =        0.6701690148797548
 Exc =       -3.2322073647096881
 Exc_core =        0.0336749726267480
 KE =        5.1951952914789290
 Etot =       -8.0896852088021269
 TS =        0.0000000000000000
  F =       -8.0896852088021269

pwscf: 

 total energy              =     -16.24670589 Ry
 Harris-Foulkes estimate   =     -16.24670590 Ry
 estimated scf accuracy    <       0.00000002 Ry
 The total energy is the sum of the following terms:

 one-electron contribution =     -50.01871785 Ry
 hartree contribution      =      25.80242008 Ry
 xc contribution           =      -6.46395695 Ry
 ewald contribution        =      14.43354883 Ry
 smearing contrib. (-TS)   =      -0.00000000 Ry

But looking back to the energy differences between them it really does just look like a Exc_core difference for both the CO and CH4 calculations.

For CH4:

-16.24670589 Ry/2 - -8.0896852088021269 = 0.0336677 hartrees (Exc_core = 0.0336749726267480 hartrees)

hartree diff of 0.00000727262675 and for CO:

-43.69009601 Ry/2 - -21.7463450252356481 = 0.09870297 hartrees (Exc_core = 0.0987190434473193 hartrees)

hartree diff of 0.00001607344732

both hartree diffs are way too small when converted back into eV. So I think consistency is met with the subtraction of Ecore_xc from the total system in jdftx!

Thanks again!