core electron binding energy (cebe)

[Alessandro Motta]

dear all,

  I am trying to perform some calculation tests for the evaluation of the core electron binding energy within the delta-SCF approach. I used a simple C9H18 linear alkene molecule, focusing on the C 1s signal. I performed the neutral ground state job and the core ionized job, using the swap option in vector and the max_ovl option to avoid function collapse.

- First of all I obtained a 280 eV for the C 1s signal (using two dft functionals: b3lyp and pbe0) that seem lower then experimentally expected (290 eV??)

- Second: some core molecular orbitals are delocazized in two or more C centers. If I interested on a specific C atom, can I localized the molecular orbitals in any way? can I assign a specific basis set to that atom? can I check where core molecular orbitals are located?

- third: I tried to perform similar calculation in the scf module (instead of dft) using the lock keyword instead of the max_ovl, but, the wave function collapses. On the other hand, if I use the HFexch in the dft section, the job finishes correctly. 

- 4th: I used different basis set (6-311g**, 6-311+g** and a pcX-3, obtaining similar results. Anyone can give me some suggestions about the best basis set for this kind of calculations?

I attach two input files (ground state and ionized state) for more details

Thanks!

Alessandro

[Edoardo Aprà]

Second: some core molecular orbitals are localized in two or more C centers. If I interested on a specific C atom, can I localized the molecular orbitals in any way? can I assign a specific basis set to that atom? can I check where core molecular orbitals are located?

You could try to replace the C atom with N, run this N-substitute molecule with a +1 charge, then using this molecular orbitals (where the symmetry between C atoms should be broken with C being replaced  by N) as a starting guess for the core ionized case

Third: I tried to perform similar calculation in the scf module (instead of dft) using the lock keyword instead of the max_ovl, but, the wave function collapses. On the other hand, if I use the HFexch in the dft section, the job finishes correctly.

the DFT module and the SCF module use different solver for the SCF solution

[Alessandro Motta]

Dear Edoardo,

thank you for your suggestions! I have tryied to perform a job in two step for the ionized state calculation. You will find the input I prepared attached. Unfortunally, the second step does not easily converge. Please find attached the output for details. Anything wrong in the input? bad basis set o dft functional?

thanks in advance

Alessandro

[Edoardo Aprà]

The attached input file shows how to get the first 1s hole using the z+1 starting guess combined with the occup input option. https://nwchemgit.github.io/Density-Functional-Theory-for-Molecules.html#occup-controlling-the-occupations-of-molecular-orbitals

[Alessandro Motta]

I tried to use occup instead of swap, but I found similar convergence instability. Probably my input have something wrong. I attach both the new input with the keyword OCCUP and the relative output for details. Sorry for my poor ability with the setup of the jobs, but I am a newbie with nwchem

Alessandro

[Edoardo Aprà]

Dumb suggestion: could you try first the input I sent you, check it works on your installation and then modify it by steps.

[Edoardo Aprà]

The main problem with your last input was the use of projection in the scf guess. It is not needed.

 vectors input  C9H18N+1_pbe0.movecs output C9H18+1_pbe0.movecs

Since the number of basis function does not change between the z+1 case and "regular" +1 case, you do not need any projection. You do want to use the very same molecular orbitals.

See the attached modified input file (with minimal changes with your last one).

Please not also the single definition of the basis set.

basis "ao basis" spherical noprint
  H  library "6-311+g**"
  C  library "6-311+g**"
  N  library "6-311+g**"
end

[Edoardo Aprà]

Please ignore the cocc.db file that was sent by mistake

[Alessandro Motta]

I tried the input you attached in the previous message, but a basis mismatch error occurred in reading the vector at the beginning of the second step. Maybe for this reason (I don't rememeber) I decided to use the project option. Attached you will find the output. As you can see, the job at the secon step loads an atomic guess instead of the vector produced at the first step.

[Edoardo Aprà]

I have missed the basis set mismatch (I was running my test with a smaller basis).

To overcame this problem, you have to add the following line before the task line

set lindep:n_dep 0

For example

...

  noprint "final vectors analysis"
end

set lindep:n_dep 0

task dft

You need to re-run the neutral case with the same condition in order to get consistent energy differences.

[Alessandro Motta]

Thank Edoardo for all your suggestions and patience! now the job works well. In order to better understand the keywords used in this work I have two other questions: 1) what is the difference between occup and swap? I used both options obtaining results strongly different (about 10 eV in the C 1s binding energy). 2) can you suggest me any readings to deep inside the meaning of the lindep option?

thanks again and all the best

Alessandro

[Edoardo Aprà]

You can find the documentation for the swap and occup keywords on the NWChem website

https://nwchemgit.github.io/Density-Functional-Theory-for-Molecules.html#occup-controlling-the-occupations-of-molecular-orbitals

https://nwchemgit.github.io/Hartree-Fock-Theory-for-Molecules.html#vectors-swap-keyword

One difference between the two approaches is that while the swap keyword does not apply any constraint on the orbital occupation, the occup keyword enforces the orbital occupation set by input at every SCF cycle. My suggestion is to look at the final eigenvectors (aka molecular orbitals) printed at the end of each run and compare what you got for the two choices.

Suggest reading for lindep? I am afraid we have never really documented this option in the NWChem user manual (and this might be chance to do it ... users' contribution are welcome!)

You could read what I posted just a few days ago in another thread https://groups.google.com/g/nwchem-forum/c/FG892T6MLRU/m/HzShIvHcBwAJ

Two or more basis functions can be consider linearly dependent when they span the same region of space. This can result in SCF converge problems. Analysis of the eigenvectors of the S-1/2 matrix (where S is the overlap matrix) is used to detect linear dependencies: if there are eigenvalues  close to zero, the basis set goes through the process of canonical orthogonalization (as describe in Section 3.4.5 of Szabo & Ostlund "Modern Quantum Chemistry" book). This has net effect of a reduction of number of basis function used, compared to the original number set by input. 

By setting set lindep:n_dep 0 this orthogonalization process is skipped.

In your specific case, the basis set you chose had some diffuse function (i.e. with small exponents). This resulted in a a different number of linear dependencies between the z+1 case and +1 case, therefore causing the molecular orbitals guess to fail.

[Alessandro Motta]

very exaustive! thanks again

Alessandro

[Alessandro Motta]

Dear Edoardo, dear all,

  I am proceeding with my test calculations for the evaluation of the core electron ionization energy (cebe). Now, the input works well but some times I obtain wrong results. I noted that job with wrong binding energy shows also wrong <S2> value compared to the expected (exact) one. I don't know if the wrong <S2> is related to the wrong binding energy. I attached two jobs (input and output). In the first one (G_e_24ch_hfexch) the binding energy (291.5 eV) is a good value with respect to the experimental one and the <S2> (0.87) is near to the exact one (one uncoupled core electron, 0.75). In the second job (G_COH3_29ch_hfexch), the binding energy (288.2 eV) is too much lower with respect the experimental one and the <S2> (5.67) is also far away compared to the exact value (0.75). I tried to use other dft functionals (pbe0 and b3lyp) obtaining a <S2> values stable at 1.8. Are binding energy wrong value related to the <S2> wrong value? does exist a way to improve the <S2> value compared to the exact one, in order to improve the accuracy of the binding energy value ?

thanks a lot!

Alessandro

[Alessandro Motta]

dear all,

  in order to better understand the setting to use for cebe calculations, I realized that I used the keywords charge 1 and mult 2, in the core-hole calculation thinking that once extracted a core electron (using the keyword occup) the system was positively charged and with a multiplicity of 2. I saw into the water example of the occup section in the nwchem webpage that for the core-hole calculation the charge remains zero and the multiplicity remains  one. I can't well understand: if I set an hole in the core level with the occup keyword, how can the system remain neutral (charge 0) and he multiplicity remain 1? Where is the extracted electron gone? Sorry for my naive questions.

Alessandro

[Edoardo Aprà]

Once occup is applied, charge and multiplicity become whatever the new occupation numbers correspond to.

[Alessandro Motta]

OK clear, thanks! one more question: in order to avoid spin contamination I would to use rodft instead of odft for the core hole calculation, but I am not sure if the max_ovl (useful to avoid function collapse after that occup set the right position of the core hole) works when rodft is requested.

[Edoardo Aprà]

max_ovl is no longer needed when occup is used.

[Alessandro Motta]

I triyed to uses occup e rodft for the core hole calculation (in order to avoid spin contamination), but it seems that when rodft is requested, occup is not considered since the calculation converges to the ground state energy.

[Edoardo Aprà]

Yes. the RODFT code does not have the same functionality as the ODFT code. OCCUP is one of the missing features.

How bad is the spin-contamination you are getting with ODFT?

[Alessandro Motta]

it depends. Using b3lyp e pcX-1 as basis set, I found a <S2> value of about 1.0 (instead of 0.75) with small molecules (benzene, hexane and derivatives), but I found a <S2> in the range 1.5-2.1 with little graphene sheets (about 50 carbon atoms) with some oxygen functionalities (epoxide, hydroxyl, carbonyl..). I attach one output file as example of the entire calculation process (ground state/z+1 noscf/core hole).