Spin Coupling in Fe4S4 Cluster¶

Step 1: Create and pre-optimize the Fe₄ S₄ cubane model¶

Embedded into proteins and coordinated by cysteine ligands, iron-sulfur cubanes are often used by nature in electron transfer and catalysis. While their native structures can be extracted from protein data bank (PDB) files, we will use a model of the Fe₄ S₄ cubane in this example.

In this step are the instructions to build the system. If you just want the final system go to the end of this step, there you can copy-paste the coordinates.

Start ADFinput

Select the Structure tool → Polyhedra → Cube structure

Click anywhere in the empty structure drawing area

The cube built is constructed of carbon atoms. We will now change four of these atoms to iron (Fe), and next the remaining four atoms into sulfur (S).

Via Shift-click, select four carbon atoms in the corners of the cubane

Atoms → Change Atom Type → Fe: change the four atoms into Fe

Selection → Invert Selection

Atoms → Change Atom Type → S: change the other atoms into S

Now you should see the Fe₄ S₄ cubane in the structure drawing area of ADFinput. The proper coordination to Fe atoms is important in modeling their electronic structure. In proteins, iron-sulfur cubanes are coordinated by cysteine ligands to the Fe sites. Here, we will model these four cysteines by ^- SH ligands. The procedure to add these ligands is described below.

Select all Fe atoms (for example by inverting the selection)

Use the Atoms → Details (Color, Radius, Mass, ...) menu command

Change the number of connectors for all Fe atoms from 10 into 4

Switch to the Main panel

Add hydrogens to the (still selected) Fe atoms using the Atoms → Add Hydrogen menu command

Select the (newly added) hydrogens

Use the Atoms → Replace By Structure → Ligands → OH menu command:

This will replace hydrogen atoms into OH ligands.

Select one of the O atoms

Use the Select → Select Atoms Of Same Type command to select them all

Use the Atoms → Change Atom Type → S to change them into S

Instead of building the structure using the instructions above, you can also copy-paste the following coordinates:

Fe      0.92466000       0.92466000      -0.92466000
S       0.92466000       0.92466000       0.92466000
Fe      0.92466000      -0.92466000       0.92466000
S       0.92466000      -0.92466000      -0.92466000
S      -0.92466000       0.92466000      -0.92466000
Fe     -0.92466000       0.92466000       0.92466000
S      -0.92466000      -0.92466000       0.92466000
Fe     -0.92466000      -0.92466000      -0.92466000
H       2.61053000       2.61053000      -2.61053000
H       2.61053000      -2.61053000       2.61053000
H      -2.61053000       2.61053000       2.61053000
H      -2.61053000      -2.61053000      -2.61053000
S      -2.03318000      -2.03318000      -2.03318000
S      -2.03318000       2.03318000       2.03318000
S       2.03318000      -2.03318000       2.03318000
S       2.03318000       2.03318000      -2.03318000

Step 2: Optimize the structure with ADF¶

The next step is to optimize this structure. It is a difficult system, and the pre-optimizers will fail. So we will use ADF to optimize the geometry.

Set the charge to -2

Select the Geometry Optimization task

Click the Symmetrize button (the star) to check the symmetry (should be Td)

File → Run (enter FeS as name)

When ready, click OK to use the optimized geometry in ADFinput

SCM → Movie

In the Movie window: Graph → Energy

Click the Play button

Close the movie window: File → Close

Step 3: Obtain the solution for the high-spin (HS) state of the cubane¶

In Fe₄ S₄ systems, the iron sites are commonly high-spin ferrous (Fe³⁺ , S = 5/2) or ferric (Fe²⁺ , S = 2). For the present example, we will use the iron-sulfur cubane oxidation state where the two sites are ferric and the remaining two are ferrous. This oxidation level of Fe₄ S₄ is well-defined and occurs, for example, in rubredoxin and high-potential iron-sulfur proteins (HIPIPs). For our model system, Fe₄ S₄ (SH)₄ , this implies the total charge of -2.

The relative directions, or coupling, of the individual site spin vectors is a very important issue in obtaining the desired density functional solution in Fe₄ S₄ , as well as many other systems which display a multi-center radical character.

Within the common open-shell approach, the spin vectors are either parallel or anti-parallel. The case when all the spins are parallel is called high-spin (HS). Obtaining self-consisted field (SCF) solution for the HS case is normally simpler because the program does not need to resolve the ambiguity of distribution the sites spin vectors. While the ferromagnetic HS solution commonly does not correspond to the lowest energy electronic state, this solution can be used to obtain the electron density corresponding to the lower energy spin-coupled state. In this step, we will obtain the HS solution for the iron-sulfur cubane, which will be used later in the tutorial. The HS solution for the [Fe₄ S₄ (SH)₄ ]^2- model corresponds to S = 2 × 5/2 + 2 × 2 = 9.

Select the ADFinput window

Select the Single Point task

Keep the total charge at -2

Set the spin polarization to 18 (corresponds to S = 9)

Check the unrestricted box

File → Save As...

Enter FeS_HS as filename and save

File → Run

Click on the progress lines so that ADFtail window will pop up showing the progress of the job:

The logfile shows that the convergence has been obtained in about 15 SCF cycles.

Step 3: Couple the spins in Fe₄ S₄ using the SpinFlip option¶

While SCF solution for the high-spin (HS) state of a multi-center spin system can often be easily found, this solution does not necessarily correspond to the lowest energy state.

To generate the solution with the desired collinear spin arrangement, one option is to use the ‘SpinFlip’ concept that has been earlier introduced by L. Noodleman and coworkers. In this two step procedure:

• first the spin-unrestricted HS solution is generated with all the site spins ferromagnetically coupled in one direction (all spins up, \(\uparrow\) ).
• Next the α (\(\uparrow\) ) and β (\(\downarrow\) ) electron densities centered at the sites which are expected to be antiferromagnetically coupled (spins down, \(\downarrow\) ) to the total spin vector are exchanged for the earlier obtained HS solution, using the SpinFlip option, and the calculation is restarted.

Because this approach often results in lowering the electronic symmetry of the system while retaining its structural symmetry, a solution obtained in such way is often called the broken symmetry (BS) state. In such cases you will need to make sure that your BS calculation is done with lower symmetry.

The concept of SpinFlip and BS state can be illustrated considering our iron-sulfur cubane case with two ferrous (Fe³⁺ , S = 5/2) and two ferric (Fe²⁺ , S = 2) sites. One of the well-characterized BS states for this level of Fe₄ S₄ oxidation corresponds to S = (5/2 + 2) - (5/2 + 2) = 0, \(2 \uparrow :2 \downarrow\).

Select the ADFinput window with your FeS_HS calculation

Make sure the ‘Main’ panel is visible

Change the spin polarization from 18 to 0

This spin polarization setting corresponds to S = 0 zero spin of the BS electronic state.

Panel bar Model → Spin and Occupation

ADFlevels will show the levels from your previous calculation. As we are not changing the occupations this info is not needed for this tutorial.

Close the window showing the energy level diagram (ADFlevels ...)

Select the two (out of four) arbitrary Fe sites in the drawing area

Click + next to the Spin Flip on Restart For: line.

To achieve the desired BS solution, the SpinFlip procedure should be applied to 2 out of 4 Fe sites of Fe₄ S₄ . In the above example, we selected sites Fe(2) and Fe(7). This will instruct SpinFlip algorithm to interchange α (\(\uparrow\) ) and β (\(\downarrow\) ) electron densities associated with these two Fe sites when the job will be restarted, changing the HS \(4 \uparrow :0 \downarrow\) spin state to the BS \(2 \uparrow :2 \downarrow\) spin state.

The Spin Flip option only works when restarting, so set up the restart calculation from the HS results:

Panel bar Details → Files (Restart)

Click folder icon in front of the ‘Restart file:’ field,

Select the Fe_HS.t21 file

The above will instruct ADF to read the converged HS solution we obtained in the previous step of the tutorial. This solution has been saved in the Fe_HS.t21.

Panel bar Details → Symmetry

Set the symmetry symbol to NOSYM

File → Save As, save the job as FeS_BS_SPINFLIP

File → Run

As you see, our FeS_BS_SPINFLIP job converges in 21 cycles.

Step 4: Coupling the spins using the ModifyStartPotential option, use ARH SCF convergence method¶

There is an alternative to SpinFlip available in ADF, which is aimed to achieve a specific spin-coupled solution in a single calculation only. This is done using the MODIFYSTARTPOTENTIAL key in ADF: it allows you to create a spin-polarized potential at the very start of the calculation.

Please read the page of ADF user’s guide on MODIFYSTARTPOTENTIAL key prior to proceeding with this step of the tutorial. As follows from the MODIFYSTARTPOTENTIAL description, this key allows to control the ratio of spin-α and spin-β electrons associated with fragments via ‘alpha’ and ‘beta’ numbers. For the purpose of the present tutorial, we will consider the four Fe sites as fragments. Apparently, ‘alpha’ and ‘beta’ numbers will correspond to the number of spin-α and spin-β electrons, correspondingly, associated with a Fe site.

So follow these instructions to obtain the BS solution via the MODIFYSTARTPOPTENTIAL option:

Open your FeS_HS calculation in ADFinput

Change the spin polarization from 18 to 0

Panel bar Details → Symmetry

Set the symmetry to NOSYM

Panel bar Model → Spin and Occupations

Close the ADFlevels window

Select the four Fe atoms

Click the ‘+’ in front of ‘Modify Start Potential’

Set the alpha and beta occupations, for the four Fe atoms:

14 alpha and 10 beta

14 alpha and 9 beta

10 alpha and 14 beta

9 alpha and 14 beta

For the spin-up Fe³⁺ , S = 5/2, the alpha and beta numbers would be 14 and 9, correspondingly, and for the spin-up Fe²⁺ , S = 2, these numbers are 14 and 10. For the spin-down Fe sites, alpha and beta numbers should be apparently transposed. Also, our desired BS state for this level of Fe₄ S₄ oxidation corresponds to S = (5/2 + 2) - (5/2 + 2) = 0, \(2 \uparrow :2 \downarrow\).

Note that SCF procedure might be problematic.

You can play with several options to help the convergence. In the SCF panel (in Details) you can experiment with several methods like LISTi, A-DIIS, E-DIIS, ARH and the option to enable the old SCF algorithm. In the SCF convergence details panel further options are available, like mixing, level shifting, orbital freezing and DIIS details.

In this particular case, the default method works but it needs a lot of iterations. So no need to make changes.

File → Save As: save the job as FeS_BS_MODIFYSTARTPOTENTIAL

File → Run

Hopefully you will be able to converge the job, to the same energy and state as the SpinFlip job.

Thus both the SpinFlip and the ModifyStartPotential option allow you to obtain a desired Fe spin coupling pattern in the Fe₄ S₄ case. While SpinFlip is a two-step approach and ModifyStartPotential works as in a single step, the SpinFlip approach shows a better performance during SCF (much better and faster SCF convergence).

Step 5: View the spin density of the broken symmetry (BS) solutions¶

In the two previous steps of this tutorial, we have generated the broken symmetry (BS) solution for the Fe₄ S₄ cubane using alternatively the SpinFlip and ModifyStartPotential options. Here, we will analyze this BS solution viewing the Fe spin densities using ADFview, and confirm that the spin alignment of the iron sites is \(2 \uparrow :2 \downarrow\). This type of analysis can also be powerfully used during presentations and for scientific illustrations.

Select your FeS_BS_MODIFYSTARTPOTENTIAL calculation in ADFjobs

Use the SCM → View menu command to activate ADFview

You should see the [Fe₄ S₄ (SH)₄ ]^2- system in the ADFview window.

The spin-density is available as a short-cut in the Properties menu:

Properties → Spin Density

You should see now the two Fe ions surrounded by blue blobs (spin up, \(\uparrow\) ) and the other two by red blobs (spin down, \(\downarrow\)):

Change the iso value to 0.03

Note that you now have two visualization lines: one is a calculated field that actually calculates the difference between the alpha and beta spin density. The other is a double isosurface that visualized the calculated field. You can use these calculated fields for many purposed. For example, by changing the ‘-‘ into a ‘+’ you are visualizing the total density instead of the spin difference density.

In the same way you can also check the spin densities for the FeS_HS and FeS_BS_SPINFLIP results.