Spin contamination for calculations including iron

Clicked A Few Times

8:00:00 AM PDT - Tue, Jun 21st 2016

Hi all

I am running DFT calculation using nwchem6.6 on CentOs including Fe atoms. I am experienced with the Turbomole program, but I am now in a position where I am using NWChem. I have been trying to run calculations with iron containing heme radical complexes, but I am running into convergence problems and spin-contamination. I am now doing some tests with iron-oxide radicals.

I have tried many combinations of RODFT, ROHF, SCF and DFT. I have tried to optimize with different charges, spin states and reusing vectors. I still have large spin contamination which I never experienced to this degree with Turbomole.

      Expectation value of S2:
      <S2> =      1.3408 (Exact =     0.7500)

I have found an input script for an iron molecule (fe2dft.nw) in the NWChem FAQ. I tried to run this and found that this script also produces spin contamination:

Expectation value of S2:
      <S2> =      8.7503 (Exact =     8.7500)
     Expectation value of S2:
      <S2> =      8.7503 (Exact =     3.7500)
     Expectation value of S2:
      <S2> =      5.0004 (Exact =     0.0000)

I am wondering if it is an inherent problem with the NWChem code. I have not been able to find any publications that use NWchem to study heme complexes and I would like to know if anybody is running into similar problems.

Clicked A Few Times

8:09:31 AM PDT - Tue, Jun 21st 2016

This is one of my input files:

start
title "name"
geometry
 load ironoxide.xyz
end
basis
 Fe  library ahlrichs_vtz
 H   library 6-31G*
 C   library 6-31G*
 N   library 6-31G*
 O   library 6-31G*
 F   library 6-31G*
 S   library 6-31G*
 Cl  library 6-31G*
end
driver
 xyz geo
 maxiter 1000
end
geometry adjust
 zcoord
 end
end
charge -3

dft
 xc b3lyp
 iterations 1000
 mult 2
 vectors atomic output fe.mo
end

task dft optimize

unset dft

charge -1

dft
 xc b3lyp
 mult 2
 iterations 1000
 vectors fe.mo
end

task dft optimize

ironoxide.xyz

     2
 geometry
 Fe                    0.00000000     0.00000000     0.38403467
 O                     0.00000000     0.00000000    -1.24308515

I get the least spin contamination after optimizing with charge -3 and changing to charge -1. However, the spin contamination increases during optimization and slight distortions in the input geometry (~0.01 Å) changes increases the spin contamination even more.

     Expectation value of S2:
      <S2> =      1.7552 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.7562 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.7568 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.7572 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.7573 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.7573 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.2067 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      0.7767 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      0.8908 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.1086 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.2432 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.2975 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.3333 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.3400 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.3407 (Exact =     0.7500)
     Expectation value of S2:
      <S2> =      1.3408 (Exact =     0.7500)

Are there any suggestions on how I can decrease the spin contamination?

Clicked A Few Times

7:14:51 AM PST - Fri, Feb 17th 2017

any suggestions or updates with respect to this? Are there at least a member who could potentially run the fe2dft.nw and check if is will always produce spin contamination? I have pasted the script below:

start fe2

geometry noautoz 
  Fe  0 0 0 
  Fe 0 0 3.35
symmetry c2v
end

geometry Fe
  Fe  0 0 0
symmetry c2v
end

basis
  Fe library sto-3g
end

set atomscf:tags_z Fe
set atomscf:z      3.
dft
convergence energy 1d-5 density 1d-4 gradient 1d-3
mulliken
end

set geometry Fe
charge +3
title "Fe spin up"
dft
mult 6
vectors input atomic output feup.mos
end
task dft energy ignore

title "Fe spin down"
dft
mult -6
vectors input atomic output fedown.mos
end
task dft energy ignore

charge +6
unset geometry
title "Fe2 anti-f"
dft
iterations 999
odft
mult 1
vectors input fragment feup.mos fedown.mos output fe2.mos
convergence ncyds 999 energy 1d-6 density 1d-5 gradient 1d-4
end
task dft

Forum Regular

10:47:22 AM PST - Fri, Feb 17th 2017
There is a small bug in the code when it comes to calculating the exact value of S2 for cases where the number of beta electrons is greater than the number of alpha electrons. The code always calculates the exact S2 as 0.5(N_alpha - N_beta)(0.5*(N_alpha - N_beta) + 1). This formula is only going to be correct if the number of alpha electrons (N_alpha) is greater than the number of beta electrons (N_beta). This is why in the fe2dft.nw example the exact value goes from 8.75 to 3.75 when clearly they should be the same for multiplicities of 6 and -6. The last exact S2 value for this input is 0.0 as would be expected from the simple formula given above, but you have created an open-shell singlet and the calculated value of S2 is 5.0 as would be expected given the arrangement of the electrons. So in this case you are getting essentially no spin contamination. For your iron oxide example, when I optimized the geometry starting in the -1 state (i.e. I did not first optimize in the -3 state), my calculated S2 was 0.7745 for the optimized geometry. This seems like a reasonable amount of spin contamination to me. Doing as you suggest (first -3 then -1) I end up with a different structure with a significant amount of spin contamination. DFT does not often play well with transition metals so I would carefully check your results to see if you are getting reasonable orderings of the molecular orbitals.

Forum >> NWChem's corner >> Running NWChem