Modifying orbital occupancy

Dear Q-Chem community,

I’m new to Q-Chem and working on a system with two identical monomers separated by 20 Å. I’ve ensured that the ground-state calculations show degenerate energies for orbitals localized on each monomer. Now, I want to modify the orbital occupancy by transferring one electron from the HOMO on one monomer to the LUMO on the other (hence the charge of each monomer isn’t neutral anymore).
I tried to use the $swap_occupied_virtual feature by applying it only to the alpha electrons as follows:
$molecule
0 1
O -0.2341000 -2.0616800 1.5520400
O -6.2766000 1.5508000 -1.4368900
O -5.2904700 -2.1657900 1.0047900
N -5.7678900 -0.3078200 -0.2092700
C -1.6572700 0.6890800 -0.3489500
C -3.9890500 1.1869900 -0.9413700
C -3.0403700 0.3596900 -0.3055500
C -0.7178100 -0.1576600 0.2941600
C -2.5403000 -1.6250300 0.9862700
H -2.9304900 -2.5140200 1.4988000
C -3.4571300 -0.8005500 0.3687000
C -1.2769800 1.8523400 -1.0559500
H -0.2094400 2.1020300 -1.1152900
C -3.5831400 2.3205900 -1.6156700
H -4.3550800 2.9341700 -2.0979400
C -1.1731600 -1.3061000 0.9526600
C -5.4196300 0.8565400 -0.9004400
C -2.2196300 2.6471500 -1.6734000
H -1.9004300 3.5433800 -2.2198600
C -4.8779500 -1.1633100 0.4287400
C -0.6657400 -3.2197000 2.2216000
H -1.3656500 -2.9759400 3.0466800
H 0.2398300 -3.6891700 2.6370300
H -1.1634500 -3.9307300 1.5306900
C -7.1621800 -0.6701400 -0.1435400
H -7.3171500 -1.6643500 -0.5975400
H -7.4733700 -0.6959200 0.8798800
H -7.7381600 0.0540800 -0.6808000
H 0.2885700 0.0702000 0.2767200
H 20.2094500 -2.1020300 -1.1152900
C 23.9890600 -1.1869900 -0.9413700
C 25.4196400 -0.8565400 -0.9004300
H 19.7601100 3.6891400 2.6370200
O 26.2766100 -1.5508000 -1.4368800
C 21.2769900 -1.8523400 -1.0559500
C 22.2196400 -2.6471500 -1.6734000
H 21.9004400 -3.5433800 -2.2198500
C 23.5831500 -2.3205900 -1.6156700
H 24.3550900 -2.9341700 -2.0979400
O 20.2340900 2.0616800 1.5520400
O 25.2904700 2.1658000 1.0048000
N 25.7679000 0.3078200 -0.2092700
C 20.7178200 0.1576600 0.2941600
C 23.4571300 0.8005500 0.3687000
C 23.0403700 -0.3596900 -0.3055500
C 21.1731600 1.3061000 0.9526600
C 21.6572800 -0.6890800 -0.3489500
C 24.8779500 1.1633100 0.4287400
C 22.5403000 1.6250300 0.9862700
H 22.9304900 2.5140200 1.4988000
C 20.6656900 3.2196900 2.2216100
H 21.3656000 2.9759500 3.0467000
H 21.1633900 3.9307400 1.5307000
C 27.1621900 0.6701500 -0.1435300
H 27.3171600 1.6643600 -0.5975400
H 27.4733700 0.6959400 0.8799000
H 27.7381800 -0.0540700 -0.6807800
H 19.7114400 -0.0702100 0.2767200
$end

$rem
BASIS = 6-31G*
GUI = 2
METHOD = B3LYP
$end

@@@

$molecule
read
$end

$rem
BASIS = 6-31G*
GUI = 2
JOB_TYPE = SP
METHOD = B3LYP
SCF_CONVERGENCE = 8
SCF_GUESS = READ
$end

$swap_occupied_virtual
alpha 126 127
$end

Where alpha 126 is the HOMO, localized on the left monomer and alpha 127 is the LUMO, localized on the right monomer.

Inspecting the output, I found out that Q-Chem essentially just swapped the orbitals 127 with 128 (which was localized on the left monomer, as 126, but with the same shape and energy of the 127), while the 126 orbital has the same energy and shape as before. Moreover, the charge of the monomers remained neutral.

My goal is to freeze the involved orbitals and observe an actual electron transfer resulting in non-neutral monomers. Any advice on how to achieve this would be greatly appreciated. Additionally, if you know of any resources that explain the swapping feature in more detail beyond the usual manual, I would love to hear about them.

Thanks in advance for your help!

Best,
Arianna Cantarella, University of Parma.

P.S. I also tried constrained DFT to enforce charge separation but couldn’t get the calculations to converge, also enforcing the convergence with different methods.

Have you verified that the MOs are truly localized? On a symmetric dimer system you could easily get MOs that are equal linear combinations of left and right orbitals.

Having verified that, the way that I would do this is with a $occupied section,

$occupied
1:125 127
1:126
$end

where the first line is for alpha (transfer electron from 126 → 127) and the second line is for beta (normal aufbau occupation). You can use this with SCF_GUESS = READ following the regular SCF that gets you the MOs, and set MOM_START = 1 to turn on maximum overlap method. (I would also set MOM_METHOD = IMOM.)

2 Likes

Thank you very much for your helpful suggestion.

I’d like to explore this subject further, but I haven’t been able to find a detailed explanation of the algorithm in the Q-Chem manual. Do you know where to find any references or documentation where the procedure for the occupation constraint is described in more detail?

It’s not an occupation “constraint”, it’s a non-aufbau initial guess with a modified SCF procedure (MOM) that attempts to converge a non-aufbau solution. For $occupied, see here:
https://manual.q-chem.com/latest/sect_occupied.html
For MOM, see here:
https://manual.q-chem.com/latest/sect_mom.html

1 Like