Integrated absolute spin density changes significantly during geometry optimization, and is not reproduced in a subsequent SCF run (no smearing)

[Ananth Govind Rajan]

Dear all,

I am using version 2.5.1 of CP2K for geometry optimizing my system of interest, which is a graphene nanoribbon. I have attached the input and output files for two separate cases: (a) a geometry optimization (geo_opt.inp and geo_opt.out) and (b) an SCF run with the optimized geometry (scf.inp and scf.out). Two XYZ files are also attached, the initial unoptimized file (zigzag_init.xyz), and the final optimized coordinates (optimized_zigzag_initial.xyz) 

As you will see, I am performing a spin polarized (UNRESTRICTED_KOHN_SHAM) calculation, with the MULTIPLICITY set to 1. The RELAX_MULTIPLICITY option is not turned on (it is set to zero by default). I am not using smearing. As the geometry optimization proceeds, the integrated absolute spin density (IASD) seems to increase significantly (from ~0 to ~18). Is this expected? I would imagine that the IASD is related to the specified multiplicity and therefore should not change significantly as optimization proceeds, because smearing is not turned on.

In fact, after geometry optimization finished, I ran a simple SCF run to compute the energy of the system using the optimized coordinates. Surprisingly, this energy is significantly different (~ 5 eV) than the energy obtained at the end of the geometry optimization step, despite having similar input files. Moreover, the IASD is also different at the end of the SCF run (~ 0), compared to at the end of the GEO_OPT run (~18).

Can anyone please advise what the issue might be? Why are the results so different between the SCF and GEO_OPT? Is it okay that the IASD changes during GEO_OPT?

Please let me know if you need any more information for diagnosing the issue.

Thank you very much

Best regards,

Ananth Govind Rajan

PhD Candidate

Department of Chemical Engineering, MIT

[Marcella Iannuzzi]

Dear Ananth,

the total integrated spin density can be different from zero, even if the multiplicity is stil one, this should not be a problem. 

After the geometry optimisation, did you try to restart also the wave function in addition to the coordinates?

I have noticed that in your scf.inp the SCF_GUESS is set to ATOMIC.

My guess is that the scf ends to be trapped in a metastable state ,due to the bad starting guess.

Kind regards,

Marcella 

[Ananth Govind Rajan]

Dear Marcella,

Thank you very much for your response.

I have now tried using the optimized coordinates along with SCF_GUESS set to RESTART and RANDOM. When I set it to RESTART, as you suggested, I am obtaining the same SCF energy as was at the end of the geometry optimization run, which is great. However, when I use SCF_GUESS as RANDOM, I am still getting a different SCF energy (the value is the same as what I got with SCF_GUESS ATOMIC).

This issue was causing problems when running NEB calculations, where if I was not using OPTIMIZE_END_POINTS, the end points were stuck in the metastable state you pointed out. I am now using OPTIMIZE_END_POINTS so that the end points relax to the final stable state. I will also try using SCF_GUESS RESTART for the end points in the NEB calculations.

On a related note, is smearing the only technically correct way to implement UKS if the ground state multiplicity is not know apriori? For example, it might be too computationally prohibitive to specify the multiplicity and check that the corresponding SCF energy is the lowest among various multiplicities for each image in an NEB calculation.

Thank you, and kind regards,

Ananth

[Marcella Iannuzzi]

Hi

this is not surprising, the RANDOM guess is even worse than the ATOMIC guess. 

And yes, if you have no clue on the multiplicity, the safer way is using the smearing.

Though this restricts to the diagonalisation scheme for the optimisation.  

As an alternative, you can run different NEBs at different fixed multiplicities and compare the energy profiles.

Kind regards,

Marcella

[Jannen]

Hello Marcella,

Can I ask a question follow up your explanation? How could the total integrated spin density be not zero when the multiplicity is 1? 

Thanks very much,

Jannen  

[Marcella Iannuzzi]

Dear Jannen,

the electrons with antiparallel spin may occupy different orbitals, having two singly occupied orbitals. 

Example two separate H atoms with antiparallel spin.

Regards

Marcella

[Jannen]

Dear Marcella,

Thanks so much for your response. Yes, the broken symmetry system. For the example, two separate H atoms, is the Integrated absolute spin density equal to  2, since there are two unpaired electrons if there is no spin contamination? Do you have any idea how is the Integrated absolute spin density calculated? 

I have a crystal structure that has 16 Co atoms, other atoms are C, H, O. The formal oxidation state of Co is 2+, so each Co is d7. When I do high spin calculation (all Co has 3 unpaired electrons and all are spin up), the Integrated absolute spin density is about 53.617 and S**2 is about 601.56 (ideally should be 600, since 48 unpaired electrons in total). 

I am confused by the Integrated absolute spin density. Why is it 53.617.  Does my system has strong spin contamination (S**2 is 1.56 larger than the ideal value)? 

I really appreciated your help!

Best regards,

Jannen