Dear Mr. Papior, it is still a bit strange: In both calculations with the mixing modifications, but one with the modified code, dH still jumps --- Without Code Modification ---- scf: 8 -1804868.778352 -1804868.732982 -1804868.733022 0.001146 -2.812871 0.312873 scf: 9 -1804868.749710 -1804868.735597 -1804868.735641 0.000361 -2.813709 0.234864 Geom step, scf iteration, dmax: 0 10 0.017262 redata: Recompute H after scf cycle = F siesta: DUscf = 413.344772 iscf Eharris(eV) E_KS(eV) FreeEng(eV) dDmax Ef(eV) dHmax(eV) scf: 1 -1804421.557777 -1804719.910852 -1804719.911018 0.050963 -2.853464 62.730089 scf: 2689464237.365204228548699.015671228548699.015350 51.191287-18.050496********** scf: 3 -1807335.880661 -1804834.918607 -1804834.918839 51.191233 -2.792301 28.033179 scf: 4 -1804850.731953 -1804845.865303 -1804845.865456 0.010372 -2.813543 24.742072 --- With Code Modification ---- scf: 8 -1804868.778352 -1804868.732982 -1804868.733022 0.001146 -2.812871 0.312873 scf: 9 -1804868.749710 -1804868.735597 -1804868.735641 0.000361 -2.813709 0.234864 Geom step, scf iteration, dmax: 0 10 0.017262 redata: Recompute H after scf cycle = F siesta: DUscf = 413.344772 iscf Eharris(eV) E_KS(eV) FreeEng(eV) dDmax Ef(eV) dHmax(eV) scf: 1 -1804421.557777 -1804719.910852 -1804719.911018 0.050963 -2.853464 62.730089 scf: 2 -1804742.932702 -1804731.529616 -1804731.529750 0.002985 -2.849860 60.125657 scf: 3 -1804889.097349 -1804864.557415 -1804864.557459 0.043795 -2.827967 11.121388 scf: 4 -1804867.632473 -1804866.689589 -1804866.689632 0.003100 -2.820711 8.032126 interestingly, the 2nd Jump is less severe and thereafter it converges quickly. So that modification is clearly useful. But the problem is still there. Concerning the mixing, Siesta prints: ps356562: cat GST.out | grep mixing mix.SCF: Pulay mixing = Pulay mix.SCF: Linear mixing weight = 0.020000 mix.SCF: Pulay mixing = Pulay mix.SCF: Linear mixing weight = 0.020000 mix.SCF: Spin-component mixing all mix.SCF: Pulay mixing = Pulay mix.SCF: Linear mixing weight = 0.000100 mix.SCF: Spin-component mixing all the last step is the restarted (soc from soc) one. Apparently it applies the linear mixing correctly. BUT: I also tried with permanent Linear mixing, and setting both the initial mixing weight and the later mixing weight to zero. However, this did not circumvent the initial dH jump (although, with the code modification, the 2nd dDmax jump to 51.* would probably not have happened) iscf Eharris(eV) E_KS(eV) FreeEng(eV) dDmax Ef(eV) dHmax(eV) scf: 1 -1804421.557777 -1804719.910852 -1804719.911018 0.050963 -2.853464 62.730089 scf: 2689464237.365204228548699.015671228548699.015350 51.191287-18.050496********** scf: 3228548699.015671228548699.015671228548699.015350 0.000000-18.050496********** scf: 4228548699.015671228548699.015671228548699.015350 0.000000-18.050496********** scf: 5228548699.015671228548699.015671228548699.015350 0.000000-18.050496********** scf: 6228548699.015671228548699.015671228548699.015350 0.000000-18.050496********** scf: 7228548699.015671228548699.015671228548699.015350 0.000000-18.050496********** scf: 8228548699.015671228548699.015671228548699.015350 0.000000-18.050496********** scf: 9228548699.015671228548699.015671228548699.015350 0.000000-18.050496********** Geom step, scf iteration, dmax: 0 10 5826.262685 Interestingly, the dDmax, Efermi and possibly dH does not change anymore after the initial 2 jumps, as it is supposed to be for zero mixing weight. But since the initial jumps are still there, i would guess that the problem either lies in the very initialization / reloading of the density matrix / Hamiltonian, or that the first mixing is not performed as intended. Is it possible to load H directly ?