Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

first try at implementing bond-order calculations #672

Merged
merged 5 commits into from
Mar 1, 2024
Merged

Conversation

zerothi
Copy link
Owner

@zerothi zerothi commented Jan 21, 2024

This fixes #507 by adding the Wiberg and Mayer implementations.

I am not fully sure the Mayer is correct, since the expressions are PS_abPS_ba, but due to symmetry, I would expect this to be PS_ab ** 2.

@tfrederiksen would you please have a look at whether this makes sense? I am a bit puzzled about the Mayer implementation, it just looks the same as Wiberg, except the overlap matrix, and perhaps Mayer should be the default.

  • Closes Bond orders from DensityMatrix #507
  • Added tests for new/changed functions?
  • Ran isort . and black . at top-level
  • Documentation for functionality in docs/
  • Changes documented in CHANGELOG.md

Copy link

codecov bot commented Jan 21, 2024

Codecov Report

Attention: Patch coverage is 84.73282% with 20 lines in your changes are missing coverage. Please review.

Project coverage is 75.71%. Comparing base (2836e38) to head (b65bbe8).
Report is 3 commits behind head on main.

Files Patch % Lines
src/sisl/physics/densitymatrix.py 81.73% 19 Missing ⚠️
src/sisl/_core/sparse.py 92.30% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main     #672      +/-   ##
==========================================
+ Coverage   75.68%   75.71%   +0.02%     
==========================================
  Files         398      398              
  Lines       50418    50547     +129     
==========================================
+ Hits        38161    38270     +109     
- Misses      12257    12277      +20     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

src/sisl/io/siesta/_help.py Fixed Show fixed Hide fixed
src/sisl/sparse.py Fixed Show fixed Hide fixed
src/sisl/physics/densitymatrix.py Fixed Show fixed Hide fixed
@tfrederiksen
Copy link
Contributor

Thanks @zerothi, apologies for this late response. I just took a look now and tested the methods for a benzene molecule. This is what I got:

import sisl

dm = sisl.get_sile('benzene/RUN.fdf').read_density_matrix()

def calcBO(method):
    print(f"method={method}:")
    bo = dm.bond_order(method=method)
    for i in range(bo.na):
        print(bo.atoms[i].symbol, end=" ")
        #for j in range(bo.na):
        for j in range(i + 1):
            print(f"{bo[i, j]:5.2f}", end=" ")
        print()

calcBO("mayer")
calcBO("wiberg")
calcBO("bond+anti")

that produces

method=mayer:
C  0.42 
C  0.01  0.39 
C  0.01  0.00  0.40 
C  0.00  0.01  0.00  0.40 
C  0.00  0.00  0.01  0.00  0.39 
C  0.00  0.00  0.00  0.01  0.01  0.42 
H  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
H  0.00  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
H  0.00  0.00  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
H  0.00  0.00  0.00  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
H  0.00  0.00  0.00  0.00  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
H  0.00  0.00  0.00  0.00  0.00  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
method=wiberg:
C 25.27 
C  6.24 25.27 
C  6.25  0.71 25.27 
C  0.70  6.24  0.26 25.27 
C  0.71  0.26  6.24  0.71 25.27 
C  0.26  0.71  0.70  6.25  6.24 25.27 
H  5.26  0.65  0.65  0.01  0.01  0.00  8.78 
H  0.65  5.27  0.01  0.65  0.00  0.01  0.18  8.78 
H  0.65  0.01  5.27  0.00  0.65  0.01  0.18  0.00  8.78 
H  0.01  0.65  0.00  5.27  0.01  0.65  0.00  0.18  0.00  8.78 
H  0.01  0.00  0.65  0.01  5.27  0.65  0.00  0.00  0.18  0.00  8.78 
H  0.00  0.01  0.01  0.65  0.65  5.26  0.00  0.00  0.00  0.18  0.18  8.78 
method=bond+anti:
C  1.13 
C  0.27  1.13 
C  0.27  0.01  1.13 
C  0.01  0.27 -0.00  1.13 
C  0.01 -0.00  0.27  0.01  1.13 
C -0.00  0.01  0.01  0.27  0.27  1.13 
H  0.20 -0.00 -0.00  0.00  0.00  0.00  0.42 
H -0.00  0.20  0.00 -0.00  0.00  0.00 -0.00  0.42 
H -0.00  0.00  0.20  0.00 -0.00  0.00 -0.00 -0.00  0.42 
H  0.00 -0.00  0.00  0.20  0.00 -0.00 -0.00 -0.00  0.00  0.42 
H  0.00  0.00 -0.00  0.00  0.20 -0.00 -0.00  0.00 -0.00 -0.00  0.42 
H  0.00  0.00  0.00 -0.00 -0.00  0.20  0.00 -0.00 -0.00 -0.00 -0.00  0.42

I don't quite understand this as for benzene I would expect to find (quoting from Wikipedia/Bond_order):

In benzene, the delocalized molecular orbitals contain 6 pi electrons over six carbons, essentially yielding half a pi bond together with the sigma bond for each pair of carbon atoms, giving a calculated bond order of 1.5 (one and a half bond).

None of the three methods provide approx. the value 1.5 for the C-C bonds. Do you understand what is going on?

@tfrederiksen
Copy link
Contributor

tfrederiksen commented Feb 2, 2024

PS If I read the DM from a spin-polarized calculation I get completely different results for wiberg (although benzene is closed shell). mayer and bond+anti appear to be identical to the unpolarized calculation reported above.

method=mayer:
C  0.42 
C  0.01  0.39 
C  0.01  0.00  0.40 
C  0.00  0.01  0.00  0.40 
C  0.00  0.00  0.01  0.00  0.39 
C  0.00  0.00  0.00  0.01  0.01  0.42 
H  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
H  0.00  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
H  0.00  0.00  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
H  0.00  0.00  0.00  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
H  0.00  0.00  0.00  0.00  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
H  0.00  0.00  0.00  0.00  0.00  0.02  0.00  0.00  0.00  0.00  0.00  0.26 
method=wiberg:
C  0.14 
C  0.05  0.15 
C  0.05  0.00  0.15 
C  0.00  0.05  0.01  0.15 
C  0.00  0.01  0.05  0.00  0.15 
C  0.01  0.00  0.00  0.05  0.05  0.14 
H  0.06  0.00  0.00  0.00  0.00  0.00  0.07 
H  0.00  0.06  0.00  0.00  0.00  0.00  0.00  0.07 
H  0.00  0.00  0.06  0.00  0.00  0.00  0.00  0.00  0.07 
H  0.00  0.00  0.00  0.06  0.00  0.00  0.00  0.00  0.00  0.07 
H  0.00  0.00  0.00  0.00  0.06  0.00  0.00  0.00  0.00  0.00  0.07 
H  0.00  0.00  0.00  0.00  0.00  0.06  0.00  0.00  0.00  0.00  0.00  0.07 
method=bond+anti:
C  1.13 
C  0.27  1.13 
C  0.27  0.01  1.13 
C  0.01  0.27 -0.00  1.13 
C  0.01 -0.00  0.27  0.01  1.13 
C -0.00  0.01  0.01  0.27  0.27  1.13 
H  0.20 -0.00 -0.00  0.00  0.00  0.00  0.42 
H -0.00  0.20  0.00 -0.00  0.00  0.00 -0.00  0.42 
H -0.00  0.00  0.20  0.00 -0.00  0.00 -0.00 -0.00  0.42 
H  0.00 -0.00  0.00  0.20  0.00 -0.00 -0.00 -0.00  0.00  0.42 
H  0.00  0.00 -0.00  0.00  0.20 -0.00 -0.00  0.00 -0.00 -0.00  0.42 
H  0.00  0.00  0.00 -0.00 -0.00  0.20  0.00 -0.00 -0.00 -0.00 -0.00  0.42 

@tfrederiksen
Copy link
Contributor

tfrederiksen commented Feb 2, 2024

As a reference, running ORCA for benzene produces these Mayer charges and bond orders in the output file:

  NA   - Mulliken gross atomic population
  ZA   - Total nuclear charge
  QA   - Mulliken gross atomic charge
  VA   - Mayer's total valence
  BVA  - Mayer's bonded valence
  FA   - Mayer's free valence

  ATOM       NA         ZA         QA         VA         BVA        FA
  0 C      6.1010     6.0000    -0.1010     3.9514     3.9514    -0.0000
  1 C      6.1012     6.0000    -0.1012     3.9515     3.9515    -0.0000
  2 C      6.1014     6.0000    -0.1014     3.9515     3.9515    -0.0000
  3 C      6.1014     6.0000    -0.1014     3.9515     3.9515    -0.0000
  4 C      6.1012     6.0000    -0.1012     3.9515     3.9515     0.0000
  5 C      6.1010     6.0000    -0.1010     3.9514     3.9514    -0.0000
  6 H      0.8988     1.0000     0.1012     0.9652     0.9652    -0.0000
  7 H      0.8988     1.0000     0.1012     0.9652     0.9652     0.0000
  8 H      0.8988     1.0000     0.1012     0.9652     0.9652    -0.0000
  9 H      0.8988     1.0000     0.1012     0.9652     0.9652     0.0000
 10 H      0.8988     1.0000     0.1012     0.9652     0.9652     0.0000
 11 H      0.8988     1.0000     0.1012     0.9652     0.9652    -0.0000

  Mayer bond orders larger than 0.100000
B(  0-C ,  1-C ) :   1.4171 B(  0-C ,  2-C ) :   1.4171 B(  0-C ,  6-H ) :   0.9584
B(  1-C ,  3-C ) :   1.4172 B(  1-C ,  7-H ) :   0.9584 B(  2-C ,  4-C ) :   1.4172
B(  2-C ,  8-H ) :   0.9583 B(  3-C ,  5-C ) :   1.4171 B(  3-C ,  9-H ) :   0.9584
B(  4-C ,  5-C ) :   1.4171 B(  4-C , 10-H ) :   0.9584 B(  5-C , 11-H ) :   0.9584

I think these numbers make: ~1.5 for C-C resonant bonds and ~1 for C-H single bonds.

@zerothi
Copy link
Owner Author

zerothi commented Feb 13, 2024

Thanks @tfrederiksen I have found the problem (it had to do with how the formula should be understood, MM vs simple multiplications).

I just have to generalize it for supercells so it works for chains etc.

src/sisl/physics/densitymatrix.py Fixed Show fixed Hide fixed

from .sparse import SparseOrbitalBZSpin
from .spin import Spin

__all__ = ["DensityMatrix"]


def _get_density(DM, orthogonal, what="sum"):

Check notice

Code scanning / CodeQL

Explicit returns mixed with implicit (fall through) returns Note

Mixing implicit and explicit returns may indicate an error as implicit returns always return None.
@zerothi
Copy link
Owner Author

zerothi commented Feb 27, 2024

@tfrederiksen could you have a look again. I believe it is better now, I get this with Benzene and a Cchain:

Benzene:
C  1.39  0.09  0.10  0.09  1.39  0.99  0.07  0.01  0.01  0.01  0.07 
C  1.39  0.09  0.10  0.09  0.07  0.99  0.07  0.01  0.01  0.01 
C  1.39  0.09  0.10  0.01  0.07  0.99  0.07  0.01  0.01 
C  1.39  0.09  0.01  0.01  0.07  0.99  0.07  0.01 
C  1.39  0.01  0.01  0.01  0.07  0.99  0.07 
C  0.07  0.01  0.01  0.01  0.07  0.99 
H  0.01  0.00  0.00  0.00  0.01 
H  0.01  0.00  0.00  0.00 
H  0.01  0.00  0.00 
H  0.01  0.00 
H  0.01 
H 


Cchain:
method=mayer:
C  0.01  0.07  0.13  1.84  1.84  0.13  0.07  0.01 

@@ -57,24 +57,44 @@
__all__ = ["SparseCSR", "ispmatrix", "ispmatrixd"]


def _rows_and_cols(csr, rows=None):
"""Retrieve the sparse patterns rows and columns from the sparse matrix
def _to_coo(csr, data: bool = True, rows=None):

Check notice

Code scanning / CodeQL

Returning tuples with varying lengths Note

_to_coo returns
tuple of size 2
and
tuple of size 3
.
src/sisl/io/siesta/_help.py Fixed Show fixed Hide fixed
# the matrix is 0. Hence we just neglect that contribution.
j = latt.sc_index(sc)
r = r + Al[i] @ Bl[j]
except:

Check notice

Code scanning / CodeQL

Except block handles 'BaseException' Note

Except block directly handles BaseException.
This fixes #507 by adding the Wiberg and Mayer implementations.

I am not fully sure the Mayer is correct, since the
expressions are PS_abPS_ba, but due to symmetry, I would
expect this to be PS_ab ** 2.

fixed bond+antibond method, factor 0.5

Signed-off-by: Nick Papior <[email protected]>
Still needs a bit refinement, but it should be
there now.

Signed-off-by: Nick Papior <[email protected]>
Also corrected Wiberg (which gives bad results)

Signed-off-by: Nick Papior <[email protected]>
While the mulliken bond-order can easily be calculated
from the mulliken output, this will make it simpler.

Signed-off-by: Nick Papior <[email protected]>
While it is a bit weird to have bond_order
return orbital-order, I think it is not necessary
to create a new method to retrieve the bond-order.

Better to keep it simple.

Signed-off-by: Nick Papior <[email protected]>
@zerothi
Copy link
Owner Author

zerothi commented Mar 1, 2024

I'll merge this, I think it looks correct, @tfrederiksen if/when you get the time, a check from your side would be grateful. For now I think it is OK to put in and let people play with it :)

Thanks!

@zerothi zerothi merged commit 6844de2 into main Mar 1, 2024
8 checks passed
@zerothi zerothi deleted the 507-bond-order branch March 1, 2024 21:09
@tfrederiksen
Copy link
Contributor

I finally got to test this and it looks good to me too.

For benzene, I now get this:

benzene/RUN.fdf method=mayer
C  3.38 
C  1.38  3.38 
C  1.38  0.08  3.38 
C  0.08  1.38  0.10  3.38 
C  0.08  0.10  1.38  0.08  3.38 
C  0.10  0.08  0.08  1.38  1.38  3.38 
H  1.01  0.05  0.05  0.01  0.01  0.01  1.30 
H  0.05  1.01  0.01  0.05  0.01  0.01  0.01  1.30 
H  0.05  0.01  1.01  0.01  0.05  0.01  0.01  0.00  1.30 
H  0.01  0.05  0.01  1.01  0.01  0.05  0.00  0.01  0.00  1.30 
H  0.01  0.01  0.05  0.01  1.01  0.05  0.00  0.00  0.01  0.00  1.30 
H  0.01  0.01  0.01  0.05  0.05  1.01  0.00  0.00  0.00  0.01  0.01  1.30 
benzene/RUN.fdf method=wiberg
C  0.58 
C  0.19  0.58 
C  0.19  0.02  0.58 
C  0.02  0.19  0.04  0.58 
C  0.02  0.04  0.19  0.02  0.58 
C  0.04  0.02  0.02  0.19  0.19  0.58 
H  0.22  0.01  0.01  0.00  0.00  0.00  0.27 
H  0.01  0.23  0.00  0.01  0.00  0.00  0.00  0.27 
H  0.01  0.00  0.22  0.00  0.01  0.00  0.00  0.00  0.27 
H  0.00  0.01  0.00  0.22  0.00  0.01  0.00  0.00  0.00  0.27 
H  0.00  0.00  0.01  0.00  0.23  0.01  0.00  0.00  0.00  0.00  0.27 
H  0.00  0.00  0.00  0.01  0.01  0.22  0.00  0.00  0.00  0.00  0.00  0.27 
benzene/RUN.fdf method=mulliken
C  4.52 
C  1.08  4.52 
C  1.08  0.04  4.52 
C  0.04  1.08 -0.00  4.52 
C  0.04 -0.00  1.08  0.04  4.52 
C -0.00  0.04  0.04  1.08  1.08  4.52 
H  0.80 -0.01 -0.01  0.00  0.00  0.00  1.68 
H -0.01  0.80  0.00 -0.01  0.00  0.00 -0.00  1.68 
H -0.01  0.00  0.80  0.00 -0.01  0.00 -0.00 -0.00  1.68 
H  0.00 -0.01  0.00  0.80  0.00 -0.01 -0.00 -0.00  0.00  1.68 
H  0.00  0.00 -0.01  0.00  0.80 -0.01 -0.00  0.00 -0.00 -0.00  1.68 
H  0.00  0.00  0.00 -0.01 -0.01  0.80  0.00 -0.00 -0.00 -0.00 -0.00  1.68 

With the :spin option I just get zeros, which I suppose is the correct result for this closed-shell molecule.

For [2]-triangulene, which has a S=1/2 ground state, I find among the carbons:

2T/RUN.fdf method=mayer
C  3.38 
C  1.38  3.41 
C  1.38  0.07  3.41 
C  0.07  1.27  0.06  3.58 
C  0.07  0.06  1.27  0.06  3.58 
C  0.02  0.07  0.00  1.27  0.01  3.41 
C  0.06  0.06  0.06  1.26  1.26  0.06  3.63 
C  0.02  0.00  0.07  0.01  1.27  0.00  0.06  3.41 
C  0.00  0.02  0.01  0.07  0.00  1.38  0.06  0.01  3.38 
C  0.00  0.01  0.01  0.06  0.06  0.06  1.26  0.06  0.07  3.58 
C  0.00  0.01  0.02  0.00  0.07  0.01  0.06  1.38  0.00  0.07  3.38 
C  0.01  0.00  0.00  0.06  0.01  0.07  0.06  0.00  1.38  1.27  0.02  3.41 
C  0.01  0.00  0.00  0.01  0.06  0.00  0.06  0.07  0.02  1.27  1.38  0.07  3.41 
2T/RUN.fdf method=wiberg
C  0.58 
C  0.19  0.62 
C  0.19  0.02  0.62 
C  0.02  0.16  0.03  0.45 
C  0.02  0.03  0.16  0.01  0.45 
C  0.01  0.02  0.00  0.16  0.01  0.62 
C  0.03  0.01  0.01  0.15  0.15  0.01  0.48 
C  0.01  0.00  0.02  0.01  0.16  0.00  0.01  0.62 
C  0.00  0.01  0.01  0.02  0.00  0.19  0.03  0.01  0.58 
C  0.00  0.01  0.01  0.01  0.01  0.03  0.15  0.03  0.02  0.45 
C  0.00  0.01  0.01  0.00  0.02  0.01  0.03  0.19  0.00  0.02  0.58 
C  0.01  0.00  0.00  0.03  0.01  0.02  0.01  0.00  0.19  0.16  0.01  0.62 
C  0.01  0.00  0.00  0.01  0.03  0.00  0.01  0.02  0.01  0.16  0.19  0.02  0.62 
2T/RUN.fdf method=mulliken
C  4.52 
C  1.09  4.60 
C  1.09  0.03  4.60 
C  0.03  1.03  0.00  4.74 
C  0.03  0.00  1.03  0.02  4.74 
C  0.00  0.02 -0.00  1.03  0.00  4.60 
C  0.00  0.02  0.02  1.01  1.01  0.02  4.81 
C  0.00 -0.00  0.02  0.00  1.03  0.00  0.02  4.60 
C  0.00  0.00  0.00  0.03 -0.00  1.09  0.00  0.00  4.52 
C -0.00  0.00  0.00  0.02  0.02  0.00  1.01  0.00  0.03  4.74 
C  0.00  0.00  0.00 -0.00  0.03  0.00  0.00  1.09  0.00  0.03  4.52 
C  0.00 -0.00  0.00  0.00  0.00  0.03  0.02 -0.00  1.09  1.03  0.00  4.60 
C  0.00  0.00 -0.00  0.00  0.00 -0.00  0.02  0.03  0.00  1.03  1.09  0.02  4.60 

and some nonzero "spin-z" components:

2T/RUN.fdf method=mayer:spin
C  0.00 
C  0.00  0.03 
C  0.00  0.02  0.03 
C  0.00  0.00  0.00  0.00 
C  0.00  0.00  0.00  0.00  0.00 
C  0.00  0.02  0.02  0.00  0.00  0.03 
C  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
C  0.00  0.02  0.02  0.00  0.00  0.02  0.00  0.03 
C  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
C  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
C  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
C  0.00  0.02  0.02  0.00  0.00  0.02  0.00  0.02  0.00  0.00  0.00  0.03 
C  0.00  0.02  0.02  0.00  0.00  0.02  0.00  0.02  0.00  0.00  0.00  0.02  0.03 
2T/RUN.fdf method=wiberg:spin
C  0.00 
C  0.00  0.01 
C  0.00  0.01  0.01 
C  0.00  0.00  0.00  0.00 
C  0.00  0.00  0.00  0.00  0.00 
C  0.00  0.01  0.01  0.00  0.00  0.01 
C  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
C  0.00  0.01  0.01  0.00  0.00  0.01  0.00  0.01 
C  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
C  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
C  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00 
C  0.00  0.01  0.01  0.00  0.00  0.01  0.00  0.01  0.00  0.00  0.00  0.01 
C  0.00  0.01  0.01  0.00  0.00  0.01  0.00  0.01  0.00  0.00  0.00  0.01  0.01 
2T/RUN.fdf method=mulliken:spin
C -0.08 
C  0.01  0.38 
C  0.01 -0.01  0.38 
C  0.00  0.01 -0.00 -0.07 
C  0.00 -0.00  0.01  0.00 -0.07 
C -0.00 -0.01  0.00  0.01  0.00  0.38 
C -0.00 -0.00 -0.00 -0.00 -0.00 -0.00  0.02 
C -0.00  0.00 -0.01  0.00  0.01 -0.00 -0.00  0.38 
C  0.00 -0.00  0.00  0.00  0.00  0.01 -0.00  0.00 -0.08 
C  0.00  0.00  0.00  0.00  0.00 -0.00 -0.00 -0.00  0.00 -0.07 
C  0.00  0.00 -0.00  0.00  0.00  0.00 -0.00  0.01  0.00  0.00 -0.08 
C  0.00  0.00 -0.00 -0.00  0.00 -0.01 -0.00  0.00  0.01  0.01 -0.00  0.38 
C  0.00 -0.00  0.00  0.00 -0.00  0.00 -0.00 -0.01 -0.00  0.01  0.01 -0.01  0.38 

Not sure what to take away from the latter, though.

@zerothi
Copy link
Owner Author

zerothi commented Mar 8, 2024

Thanks for following up. Yeah, I am not sure about the spin details. But I guess they can be used to correlate differences between two structures/states.

I would assume it means something is the spin configuration is a localized or delocalized quantity.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Bond orders from DensityMatrix
2 participants