Merge pull request #17 from cclib/dev/calculation_methods

Dev/calculation methods

Author: berquist

docs/make.jl

@@ -16,6 +16,7 @@ makedocs(;
     pages=[
         "Home" => "index.md",
         "How to parse files" => "io.md",
+        "Additional Analyses" => "calculation.md",
         "AtomsBase Integration" => "atombase.md"
     ],
 )

docs/src/calculation.md

@@ -0,0 +1,85 @@
+# Additional Analyses
+
+Cclib also allows to further analyse calculation ouputs.
+
+# C squared population analysis (CSPA)
+**CSPA** can be used to determine and interpret the electron density of a molecule. The contribution of the a-th atomic orbital to the i-th molecular orbital can be written in terms of the molecular orbital coefficients:
+$$\Phi_{ai} = \frac{c^2_{ai}}{\sum_k c^2_{ki}}$$
+```Julia
+# Example calculation files can be found in the test folder of the main branch.
+julia> using Cclib
+julia> aoresults, fragresults, fragcharges = cspa("./Trp_polar.fchk")
+```
+
+* ``aoresults``: a three dimensional array with spin, molecular orbital, and atomic orbitals as the axes, so that ``aoresults[1, 46, 1]`` gives the contribution of the 1st atomic orbital to the 46th alpha/restricted molecular orbital,
+* ``fragresults``: a three dimensional array with spin, molecular orbital, and atoms as the axes, so that ``fragresults[1, 24, 5]`` gives the contribution of the 5th fragment orbitals to the 24th beta molecular orbital)
+* ``fragcharges``: a vector with the number of (partial) electrons in each fragment, so that ``fragcharges[3]`` gives the number of electrons in the 3rd fragment.
+
+# Mulliken population analysis (MPA)
+MPA can be used to determine and interpret the electron density of a molecule. The contribution of the a-th atomic orbital to the i-th molecular orbital in this method is written in terms of the molecular orbital coefficients, c, and the overlap matrix, S:
+$$\Phi_{ai} = \sum_b c_{ai} c_{bi} S_{ab}$$
+```Julia
+julia> using Cclib
+julia> aoresults, fragresults, fragcharges = mpa("./Trp_polar.fchk")
+```
+
+* ``aoresults``: a three dimensional array with spin, molecular orbital, and atomic orbitals as the axes, so that ``aoresults[1, 46, 1]`` gives the contribution of the 1st atomic orbital to the 46th alpha/restricted molecular orbital,
+* ``fragresults``: a three dimensional array with spin, molecular orbital, and atoms as the axes, so that ``fragresults[1, 24, 5]`` gives the contribution of the 5th fragment orbitals to the 24th beta molecular orbital)
+* ``fragcharges``: a vector with the number of (partial) electrons in each fragment, so that ``fragcharges[3]`` gives the number of electrons in the 3rd fragment.
+
+# Löwdin Population Analysis
+```Julia
+julia> using Cclib
+julia> aoresults, fragresults, fragcharges = lpa("./Trp_polar.fchk")
+```
+
+* ``aoresults``: a three dimensional array with spin, molecular orbital, and atomic orbitals as the axes, so that ``aoresults[1, 46, 1]`` gives the contribution of the 1st atomic orbital to the 46th alpha/restricted molecular orbital,
+* ``fragresults``: a three dimensional array with spin, molecular orbital, and atoms as the axes, so that ``fragresults[1, 24, 5]`` gives the contribution of the 5th fragment orbitals to the 24th beta molecular orbital)
+* ``fragcharges``: a vector with the number of (partial) electrons in each fragment, so that ``fragcharges[3]`` gives the number of electrons in the 3rd fragment.
+
+# Bickelhaupt Population Analysis
+The Bickelhaupt class available from cclib.method performs Bickelhaupt population analysis that has been proposed in *Organometallics* 1996, 15, 13, 2923–2931. [doi:10.1021/om950966x](https://pubs.acs.org/doi/abs/10.1021/om950966x)
+
+The contribution of the a-th atomic orbital to the i-th molecular orbital in this method is written in terms of the molecular orbital coefficients, c, and the overlap matrix, S:
+
+$$\Phi_{ai,\alpha} = \sum_b w_{ab,\alpha} c_{ai,\alpha} c_{bi,\alpha} S_{ab}$$
+
+where the weights $w_{ab}$ that are applied on the Mulliken atomic orbital contributions are defined as following when the coefficients of the molecular orbitals are substituted into equation 11 in the original article.
+
+$$w_{ab,\alpha} = 2 \frac{\sum_k c_{ak,\alpha}^2}{\sum_i c_{ai,\alpha}^2 + \sum_j c_{bj,\alpha}^2}$$
+
+In case of unrestricted calculations, $\alpha$ charges and $\beta$ charges are each determined to obtain total charge. In restricted calculations, $\alpha$ subscript can be ignored since the coefficients are equivalent for both spin orbitals.
+
+The weights are introduced to replace the somewhat arbitrary partitioning of off-diagonal charges in the Mulliken population analysis, which divides the off-diagonal charges identically to both atoms. Bickelhaupt population analysis instead divides the off-diagonal elements based on the relative magnitude of diagonal elements.
+```Julia
+julia> using Cclib
+julia> aoresults, fragresults, fragcharges = bpa("./Trp_polar.fchk")
+```
+* ``aoresults``: a three dimensional array with spin, molecular orbital, and atomic orbitals as the axes, so that ``aoresults[1, 46, 1]`` gives the contribution of the 1st atomic orbital to the 46th alpha/restricted molecular orbital,
+* ``fragresults``: a three dimensional array with spin, molecular orbital, and atoms as the axes, so that ``fragresults[1, 24, 5]`` gives the contribution of the 5th fragment orbitals to the 24th beta molecular orbital)
+* ``fragcharges``: a vector with the number of (partial) electrons in each fragment, so that ``fragcharges[3]`` gives the number of electrons in the 3rd fragment.
+
+# Density Matrix calculation
+Calculates the electron density matrix
+```Julia
+julia> using Cclib
+julia> result = density("./Trp_polar.fchk")
+```
+Returns an array with three axes. The first axis is for the spin contributions, the second and the third axes for the density matrix, which follows standard definition.
+# Mayer’s Bond Orders (MBO)
+Calculates Mayer's bond orders
+```Julia
+julia> using Cclib
+julia> result = mbo("./Trp_polar.fchk")
+```
+Returns an array with three axes. The first axis is for contributions of each spin to the MBO, while the second and the third correspond to the indices of the atoms.
+# Charge Decomposition Analysis
+The Charge Decomposition Analysis (CDA) as developed by Gernot Frenking et al. is used to study the donor-acceptor interactions of a molecule in terms of two user-specified fragments.
+```Julia
+julia> using Cclib
+julia> donations, bdonations, repulsions = cda("BH3CO-sp.log", "BH3.log", "CO.log")
+```
+Returns donations, bdonations (back donations), and repulsions attributes.
+These attributes are simply lists of 1-dimensional arrays corresponding to the restricted or alpha/beta molecular orbitals of the entire molecule.
+
+<!-- # Bader’s QTAIM -->

src/Cclib.jl

@@ -7,6 +7,7 @@ using UnitfulAtomic
 using Logging
 
 include("config.jl")
+include("calculation_methods.jl")
 include("functions.jl")
 include("ab_integration.jl")
 

src/calculation_methods.jl

@@ -0,0 +1,191 @@
+#
+# Calculation methods for qchem outputs
+#
+export cspa
+export mpa
+export density
+export lpa
+export bpa
+export density
+export mbo
+export cda
+export bader
+
+# dimensions of some outputs are NxMxK, however python return dimension N as a list, not an array
+function expand(x)
+    N = length(x)
+    x = pyconvert(Array{Float64}, x[0])
+    return reshape(x, (N, size(x)...))
+end
+
+"""
+    cspa(file::String)
+
+C-Squared Population Analysis (CSPA)
+# Arguments
+- `file::String`: Cclib-supported output file
+# Returns
+Tuple with 3 elements:
+- `aoresults`: a three dimensional array with spin, molecular orbital, and atomic orbitals as the axes, so that ``aoresults[1][46][1]`` gives the contribution of the 1st atomic orbital to the 46th alpha/restricted molecular orbital,
+- `fragresults`: a three dimensional array with spin, molecular orbital, and atoms as the axes, so that ``fragresults[1][24][5]`` gives the contribution of the 5th fragment orbitals to the 24th beta molecular orbital)
+- `fragcharges`: a vector with the number of (partial) electrons in each fragment, so that ``fragcharges[3]`` gives the number of electrons in the 3rd fragment.
+"""
+function cspa(file::String)
+    #TODO: implement the optional args from docs
+    data = cclib[].io.ccread(file)
+    mol = cclib[].method.CSPA(data)
+    mol.calculate()
+    aoresults = mol.__dict__["aoresults"] |> expand
+    fragresults = mol.__dict__["fragresults"] |> expand
+    fragcharges = pyconvert(Array{Float64}, mol.__dict__["fragcharges"])
+    return aoresults, fragresults, fragcharges
+end
+
+"""
+    mpa(file::String)
+
+Mulliken Population Analysis
+# Arguments
+- `file::String`: Cclib-supported output file
+# Returns
+Tuple with 3 elements:
+- `aoresults`: a three dimensional array with spin, molecular orbital, and atomic orbitals as the axes, so that ``aoresults[1][46][1]`` gives the contribution of the 1st atomic orbital to the 46th alpha/restricted molecular orbital,
+- `fragresults`: a three dimensional array with spin, molecular orbital, and atoms as the axes, so that ``fragresults[1][24][5]`` gives the contribution of the 5th fragment orbitals to the 24th beta molecular orbital)
+- `fragcharges`: a vector with the number of (partial) electrons in each fragment, so that ``fragcharges[3]`` gives the number of electrons in the 3rd fragment.
+"""
+function mpa(file::String)
+    data = cclib[].io.ccread(file)
+    mol = cclib[].method.MPA(data)
+    mol.calculate()
+    aoresults = mol.__dict__["aoresults"] |> expand
+    fragresults = mol.__dict__["fragresults"] |> expand
+    fragcharges = pyconvert(Array{Float64}, mol.__dict__["fragcharges"])
+    return aoresults, fragresults, fragcharges
+end
+
+"""
+    lpa(file::String)
+
+Lowdin Population Analysis
+# Arguments
+- `file::String`: Cclib-supported output file
+# Returns
+Tuple with 3 elements:
+- `aoresults`: a three dimensional array with spin, molecular orbital, and atomic orbitals as the axes, so that ``aoresults[1][46][1]`` gives the contribution of the 1st atomic orbital to the 46th alpha/restricted molecular orbital,
+- `fragresults`: a three dimensional array with spin, molecular orbital, and atoms as the axes, so that ``fragresults[1][24][5]`` gives the contribution of the 5th fragment orbitals to the 24th beta molecular orbital)
+- `fragcharges`: a vector with the number of (partial) electrons in each fragment, so that ``fragcharges[3]`` gives the number of electrons in the 3rd fragment.
+"""
+function lpa(file::String)
+    data = cclib[].io.ccread(file)
+    mol = cclib[].method.LPA(data)
+    mol.calculate()
+    aoresults = mol.__dict__["aoresults"] |> expand
+    fragresults = mol.__dict__["fragresults"] |> expand
+    fragcharges = pyconvert(Array{Float64}, mol.__dict__["fragcharges"])
+    return aoresults, fragresults, fragcharges
+end
+
+"""
+    bpa(file::String)
+
+Bickelhaupt Population Analysis
+# Arguments
+- `file::String`: Cclib-supported output file
+# Returns
+Tuple with 3 elements:
+- `aoresults`: a three dimensional array with spin, molecular orbital, and atomic orbitals as the axes, so that ``aoresults[1][46][1]`` gives the contribution of the 1st atomic orbital to the 46th alpha/restricted molecular orbital,
+- `fragresults`: a three dimensional array with spin, molecular orbital, and atoms as the axes, so that ``fragresults[1][24][5]`` gives the contribution of the 5th fragment orbitals to the 24th beta molecular orbital)
+- `fragcharges`: a vector with the number of (partial) electrons in each fragment, so that ``fragcharges[3]`` gives the number of electrons in the 3rd fragment.
+"""
+function bpa(file::String)
+    data = cclib[].io.ccread(file)
+    mol = cclib[].method.Bickelhaupt(data)
+    mol.calculate()
+    aoresults = mol.__dict__["aoresults"] |> expand
+    fragresults = mol.__dict__["fragresults"] |> expand
+    fragcharges = pyconvert(Array{Float64}, mol.__dict__["fragcharges"])
+    return aoresults, fragresults, fragcharges
+end
+
+"""
+    density(file::String)
+
+Density matrix calculation
+# Arguments
+- `file::String`: Cclib-supported output file
+# Returns
+- An array with three axes. The first axis is for the spin contributions,
+the second and the third axes for the density matrix, which follows standard definition.
+"""
+function density(file::String)
+    data = cclib[].io.ccread(file)
+    mol = cclib[].method.Density(data)
+    mol.calculate()
+    return pyconvert(Array{Float64}, mol.__dict__["density"])
+end
+
+"""
+    mbo(file::String)
+
+Calculate Mayer's bond orders
+# Arguments
+- `file::String`: Cclib-supported output file
+# Returns
+- An array with three axes. The first axis is for contributions of each spin to the MBO,
+while the second and third correspond to the indices of the atoms.
+"""
+function mbo(file::String)
+    data = cclib[].io.ccread(file)
+    mol = cclib[].method.MBO(data)
+    mol.calculate()
+    return pyconvert(Array{Float64}, mol.__dict__["fragresults"])
+end
+
+"""
+    cda(file::String)
+
+Charge decomposition analysis
+# Arguments
+- `mol::String`: Cclib-supported output file
+- `frag1::String`: Cclib-supported output file
+- `frag2::String`: Cclib-supported output file
+# Returns
+- Tuple (donations, backdonations, repulsions)
+donations, bdonations (back donations), and repulsions attributes.
+These attributes are simply lists of 1-dimensional arrays corresponding to the restricted or alpha/beta molecular orbitals of the entire molecule.
+"""
+function cda(mol::String, frag1::String, frag2::String)
+    mol = cclib[].io.ccread(mol)
+    frag1 = cclib[].io.ccread(frag1)
+    frag2 = cclib[].io.ccread(frag2)
+    mol = cclib[].method.CDA(mol)
+    mol.calculate([frag1, frag2])
+    donations = mol.__dict__["donations"] |> expand
+    bdonations = mol.__dict__["bdonations"] |> expand
+    repulsions = mol.__dict__["repulsions"] |> expand
+    return donations, bdonations, repulsions
+end
+
+# TODO: this requires PyQuante which doens't seem to work
+# when trying to install it into a conda env
+# """
+#     bader(file::string)
+
+# bader's qtaim charges calculation
+# # arguments
+# - `file::string`: cclib-supported output file
+# # returns
+# tuple (donations, backdonations, repulsions)
+# """
+# function bader(file::String, vol::Vector{Vector{Float64}})
+#     data = cclib[].io.ccread(file)
+#     vol = pytuple(pytuple(i) for i in vol)
+#     vol = cclib[].method.volume(vol...)
+#     mol = cclib[].method.bader(data, vol)
+#     mol.calculate()
+#     return mol
+#     # aoresults = mol.__dict__["aoresults"] |> expand
+#     # fragresults = mol.__dict__["fragresults"] |> expand
+#     # fragcharges = pyconvert(array{float64}, mol.__dict__["fragcharges"])
+#     # return aoresults, fragresults, fragcharges
+# end