[cclib-devel] Lowdin and generalized wavefunction projection population analysis

[Karol L. <kar...@kn...>]

Having added the LPA module for calculating Lowdin Population Analyses, I 
would like to comment on it a bit. First of all, it's a nice alternative to 
Mulliken charges. Second, it allows the calculation to be done using any 
given exponent over the overlap matrix:

q_{A} = sum_{a in A} (S^{x} P S^{1-x})_{aa}

where a are the atomic orbitals that belong to atom or group A, S is the 
overlap matrix, and P is the density matrix. When x=1.0 or x=0.0, the 
equation simplifies to Mulliken charges. The Lowdin analysis assumes x=0.5 
(default in the new code).

Take a look:

Python 2.5 (r25:51908, Apr 30 2007, 15:03:13)
[GCC 3.4.6 (Debian 3.4.6-5)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import cclib
>>> cclib.__version__
'0.7'
>>> cclib.method.mpa.__revision__
'$Revision: 627 $'
>>> mol = cclib.parser.ccopen("dvb_sp.out")
>>> mol.parse()
[Gaussian dvb_sp.out INFO] Creating attribute charge: 0
[Gaussian dvb_sp.out INFO] Creating attribute mult: 1
[Gaussian dvb_sp.out INFO] Creating attribute atomnos[]
[Gaussian dvb_sp.out INFO] Creating attribute natom: 20
[Gaussian dvb_sp.out INFO] Creating attribute atomcoords[]
[Gaussian dvb_sp.out INFO] Creating attribute nbasis: 60
[Gaussian dvb_sp.out INFO] Creating attribute aooverlaps[]
[Gaussian dvb_sp.out INFO] Creating attribute nmo: 60
[Gaussian dvb_sp.out INFO] Creating attribute scftargets[]
[Gaussian dvb_sp.out INFO] Creating attribute scfvalues[]
[Gaussian dvb_sp.out INFO] Creating attribute scfenergies[]
[Gaussian dvb_sp.out INFO] Creating attribute mosyms[]
[Gaussian dvb_sp.out INFO] Creating attribute homos[]
[Gaussian dvb_sp.out INFO] Creating attribute moenergies[]
[Gaussian dvb_sp.out INFO] Creating attribute aonames[]
[Gaussian dvb_sp.out INFO] Creating attribute mocoeffs[]
[Gaussian dvb_sp.out INFO] Creating attribute coreelectrons[]
[Gaussian dvb_sp.out INFO] Creating attribute coreelectrons[]
>>> mpa = cclib.method.MPA(mol)
>>> lpa = cclib.method.LPA(mol)
>>> mpa.calculate()
[MPA Gaussian log file dvb_sp.out INFO] Creating attribute aoresults: 
[array[2]]
[MPA Gaussian log file dvb_sp.out INFO] Saving partitioned results in 
fragresults: [array[2]]
[MPA Gaussian log file dvb_sp.out INFO] Creating fragcharges: array[1]
True
>>> lpa.calculate()
[LPA Gaussian log file dvb_sp.out INFO] Creating attribute aoresults: 
[array[2]]
[LPA Gaussian log file dvb_sp.out INFO] Saving partitioned results in 
fragresults: [array[2]]
[LPA Gaussian log file dvb_sp.out INFO] Creating fragcharges: array[1]
True
>>> mpa.fragcharges
array([ 6.0045089 ,  6.07712116,  6.07668927,  6.0045089 ,  6.07712116,
        0.92078088,  0.92230013,  0.92078088,  6.07636862,  6.15486853,
        0.92367176,  0.92311728,  0.92066214,  6.07636862,  0.92311728,
        6.15486853,  0.92367176,  0.92066214,  6.07668927,  0.92230013])
>>> lpa.fragcharges
array([ 5.99633692,  6.04422133,  6.04398657,  5.99633692,  6.04422133,
        0.9563139 ,  0.95803503,  0.9563139 ,  6.03894201,  6.09166045,
        0.95680625,  0.95882711,  0.9549591 ,  6.03894201,  0.95882711,
        6.09166045,  0.95680625,  0.9549591 ,  6.04398657,  0.95803503])
>>> lpa.calculate(x=0.0)
[LPA Gaussian log file dvb_sp.out INFO] Creating attribute aoresults: 
[array[2]]
[LPA Gaussian log file dvb_sp.out INFO] Saving partitioned results in 
fragresults: [array[2]]
[LPA Gaussian log file dvb_sp.out INFO] Creating fragcharges: array[1]
True
>>> lpa.fragcharges
array([ 6.0045089 ,  6.07712116,  6.07668927,  6.0045089 ,  6.07712116,
        0.92078088,  0.92230013,  0.92078088,  6.07636862,  6.15486853,
        0.92367176,  0.92311728,  0.92066214,  6.07636862,  0.92311728,
        6.15486853,  0.92367176,  0.92066214,  6.07668927,  0.92230013])

Notice how when x=0.0 is passed the charges are the same as the ones obtained 
using the current MPA class. It might be nice to integrate the two with a 
more general class, but not fully since when 0<x<1 two roots of the overlap 
matrix need to be calculated which takes some time.

- Karol

-- 
written by Karol Langner
Sun May 20 00:31:36 CEST 2007