# [CP2K:2981] Re: Lowdin analysis

Teodoro Laino teodor... at gmail.com
Tue Dec 7 10:23:49 UTC 2010

```Hi Matthias,

I focused on this "purely technical issues" because of your original reply. "Just to recall it":

>>>>>>>> The MOLOPT basis sets include
>>>>>>>> relatively small exponents which might cause linear dependencies in
>>>>>>>> the overlap matrix. The Lowdin population analysis requires the
>>>>>>>> calculation of S**(1/2) which involves (in CP2K) a diagonalisation of
>>>>>>>> the overlap matrix. However, an ill-conditioned overlap matrix will
>>>>>>>> most likely result in a weird charge partition.

the fact that an ill-conditioned overlap matrix will most likely result in a weird charge partition IS a purely technical issue.
A more robust scheme? It's easy: just do an SVD of the C matrix. The condition number of S is the square of the condition number of C.
This means that a condition number of 10^-14 in S, is only 10^-7 in C. This is by far numerically more stable if you work in double precision.
Hope this helps you.

Best,
Teo

On Dec 7, 2010, at 11:16 AM, Matthias Krack wrote:

> Hi Teo,
>
> you are now referring to purely technical issues concerning the
> calculation of functions of S like S**(1/2), which is not (just to
> recall it) the original topic of this thread. However, I would
> certainly appreciate if you could come up with a more robust scheme to
> calculate S**(1/2) within CP2K, accurately. However, this won't change
> the pathological results of the Lowdin analysis using MOLOPT basis
> sets.
>
> cheers,
>
> Matthias
>
> On Dec 7, 10:51 am, Laino Teodoro <teodor... at gmail.com> wrote:
>> Hi Matthias,
>>
>> I do not agree : there several cases where diagonalization of S (done to compute functions of S) leads to huge numerical noise on the eigenvalues (due to almost linear dependencies of the C vectors) - with only 1-2 significant digits if you work in double precision - whether this is that case or not, I cannot say. I trust what you said: if you checked than for sure this is not that case.
>>
>> Nonetheless, for peculiar cases (numerically speaking) skipping the diagonalization of S is the way to go to compute functions of S.
>>
>> Best
>> Teo
>>
>> ---------------------------------------------
>> Teodoro Laino
>> Zurich Switzerland
>>
>> Contact info:
>> Tel.:    http://www.jajah.com/Teo
>> E-mail: teo... at laino.eu
>>              teodor... at gmail.com
>> ---------------------------------------------
>>
>> On 7 Dec 2010, at 10:36, Matthias Krack <matthia... at psi.ch> wrote:
>>
>>> Hi Teo,
>>
>>> On Dec 7, 10:10 am, Laino Teodoro <teodor... at gmail.com> wrote:
>>>> Hi Matthias,
>>
>>>> Ok - in this case, if you checked already before your first post, why mention linear dependencies and problems with
>>>> diagonalization?
>>
>>> Since such problems exist with MOLOPT basis sets certainly for larger
>>> and especially condensed systems. Like each basis set, MOLOPT basis
>>> sets have to be used with care.
>>
>>>> These are pure numerical issues that can be solved skipping the diagonalization of S to evaluate S**(1/2).
>>
>>> As I explained, finally this doesn't help, the problem of unphysical
>>> charges will persist, even if you are able to calculate S**(1/2) in
>>> any case.
>>
>>> cheers,
>>
>>> Matthias
>>
>>>> Best
>>>> Teo
>>
>>>> ---------------------------------------------
>>>> Teodoro Laino
>>>> Zurich Switzerland
>>
>>>> Contact info:
>>>> Tel.:    http://www.jajah.com/Teo
>>>> E-mail: teo... at laino.eu
>>>>              teodor... at gmail.com
>>>> ---------------------------------------------
>>
>>>> On 7 Dec 2010, at 09:54, Matthias Krack <matthia... at psi.ch> wrote:
>>
>>>>> Hi Teo,
>>
>>>>> I checked that already before my first reply. There is no numerical
>>>>> problem in this case, i.e. there is no warning about quenched
>>>>> eigenvalues.
>>
>>>>> Matthias
>>
>>>>> On Dec 7, 9:35 am, Teodoro Laino <teodor... at gmail.com> wrote:
>>>>>> Hi Matthias,
>>
>>>>>> maybe yes - maybe no -
>>>>>> I did not look into the condition number of that matrix -
>>>>>> A definitive answer (and not just assumptions like mines or yours) would  be confirmed or discarded by a more serious check.
>>
>>>>>> Best,
>>>>>> Teo
>>
>>>>>> On Dec 7, 2010, at 9:30 AM, Matthias Krack wrote:
>>
>>>>>>> Hi Teo,
>>
>>>>>>> I doubt that the applied computational scheme or the numerics are part
>>>>>>> of the problem, especially not in the case of a single water molecule,
>>>>>>> and thus any clever approach to evaluate S**(1/2) won't change
>>>>>>> anything. Maybe, my previous reply unintentionally implied that I
>>>>>>> would blame the MOLOPT basis set to be the only problem in this case.
>>>>>>> So let it put me more clearly: just the way how the Lowdin approach is
>>>>>>> partitioning the charge does not work properly, i.e. giving physically
>>>>>>> meaningful results, for all kinds of basis sets and the MOLOPT basis
>>>>>>> sets (but not its SR variants) belong to this group.
>>
>>>>>>> Best,
>>
>>>>>>> Matthias
>>
>>>>>>> On Dec 6, 8:03 pm, Teodoro Laino <teodor... at gmail.com> wrote:
>>>>>>>> Hi Matthias, Mahn,
>>
>>>>>>>> The real problem is not in the basis set itself but in how the population analysis is computed in CP2K.
>>>>>>>> Diagonalizing S to evaluate S**(1/2) is probably not the smartest thing to do.
>>
>>>>>>>> There are ways to evaluate S**(1/2), definitely numerically more stable, that would allow to have meaningful population analysis even with these "problematic" basis sets (although I still think that it's not the basis set itself to be problematic, rather the computational scheme used).
>>
>>>>>>>> Best,
>>>>>>>> Teo
>>
>>>>>>>> On Dec 5, 2010, at 3:06 PM, Matthias Krack wrote:
>>
>>>>>>>>> Hi Manh,
>>
>>>>>>>>> population analyses are known to be very sensible to the basis set
>>>>>>>>> choice. I guess also in this case the TZV2P MOLOPT basis set is the
>>>>>>>>> reason for the unexpected result. The MOLOPT basis sets include
>>>>>>>>> relatively small exponents which might cause linear dependencies in
>>>>>>>>> the overlap matrix. The Lowdin population analysis requires the
>>>>>>>>> calculation of S**(1/2) which involves (in CP2K) a diagonalisation of
>>>>>>>>> the overlap matrix. However, an ill-conditioned overlap matrix will
>>>>>>>>> most likely result in a weird charge partition. Therefore, I would
>>>>>>>>> suggest to employ different basis sets, e.g. the MOLOPT-SR or the DZVP/
>>>>>>>>> TZVP basis sets.
>>
>>>>>>>>> Best,
>>
>>>>>>>>> Matthias
>>
>>>>>>>>> On 3 Dez., 17:46, Manh <manht... at gmail.com> wrote:
>>>>>>>>>> Hello everyone,
>>
>>>>>>>>>> I calculate the Mulliken and Lowdin charges for a water molecule. The
>>>>>>>>>> results are:
>>
>>>>>>>>>> ----------------------------------------------------------------
>>
>>>>>>>>>>  MULLIKEN POPULATION ANALYSIS
>>
>>>>>>>>>>  # Atom  Element  Kind  Atomic population                Net charge
>>>>>>>>>>       1     H        1           0.780838                  0.219162
>>>>>>>>>>       2     H        1           0.780981                  0.219019
>>>>>>>>>>       3     O        2           6.438182                 -0.438182
>>>>>>>>>>  # Total charge                  8.000000                  0.000000
>>
>>>>>>>>>>  LOWDIN POPULATION ANALYSIS
>>
>>>>>>>>>>  # Atom  Element  Kind  Atomic population                Net charge
>>>>>>>>>>       1     H        1           1.210449                 -0.210449
>>>>>>>>>>       2     H        1           1.210473                 -0.210473
>>>>>>>>>>       3     O        2           5.579079                  0.420921
>>>>>>>>>>  # Total charge                  8.000000                  0.000000
>>
>>>>>>>>>> --------------------------------------------------------
>>
>>>>>>>>>> I think that the Lowdin analysis has some problem.
>>
>>>>>>>>>> Here is my input:
>>
>>>>>>>>>> ---------------------------------------------
>>
>>>>>>>>>> &FORCE_EVAL
>>>>>>>>>>   METHOD Quickstep
>>>>>>>>>>   &DFT
>>>>>>>>>>     BASIS_SET_FILE_NAME ../../../BASIS_MOLOPT
>>>>>>>>>>     POTENTIAL_FILE_NAME ../../../GTH_POTENTIALS
>>>>>>>>>>     RESTART_FILE_NAME H2O-RESTART.wfn.wfn
>>>>>>>>>>     &QS
>>>>>>>>>>       METHOD GPW
>>>>>>>>>>       EXTRAPOLATION ASPC
>>>>>>>>>>       EXTRAPOLATION_ORDER 3
>>>>>>>>>>     &END QS
>>>>>>>>>>     &MGRID
>>>>>>>>>>       CUTOFF 400
>>>>>>>>>>       NGRIDS 5
>>>>>>>>>>     &END
>>>>>>>>>>     &SCF
>>>>>>>>>>       MAX_SCF 20
>>>>>>>>>>       SCF_GUESS RESTART
>>>>>>>>>>       EPS_SCF 1.0E-6
>>>>>>>>>>       &OT
>>>>>>>>>>         PRECONDITIONER  FULL_SINGLE_INVERSE
>>>>>>>>>>         MINIMIZER  CG
>>>>>>>>>>       &END
>>>>>>>>>>       &OUTER_SCF
>>>>>>>>>>         MAX_SCF 10
>>>>>>>>>>         EPS_SCF 1.0E-6
>>>>>>>>>>       &END
>>>>>>>>>>       &PRINT
>>>>>>>>>>         &RESTART
>>>>>>>>>>           &EACH
>>>>>>>>>>             QS_SCF 0
>>>>>>>>>>             GEO_OPT 2
>>>>>>>>>>           &END
>>>>>>>>>>           FILENAME RESTART.wfn
>>>>>>>>>>         &END
>>>>>>>>>>         &RESTART_HISTORY OFF
>>>>>>>>>>         &END
>>>>>>>>>>       &END
>>>>>>>>>>     &END SCF
>>
>>>>>>>>>>     &XC
>>>>>>>>>>       &XC_FUNCTIONAL PBE
>>>>>>>>>>       &END XC_FUNCTIONAL
>>>>>>>>>>     &END XC
>>
>>>>>>>>>>     &PRINT
>>>>>>>>>>      &MULLIKEN
>>>>>>>>>>      &END MULLIKEN
>>
>>>>>>>>>>      &LOWDIN
>>>>>>>>>>      &END LOWDIN
>>>>>>>>>>     &END
>>
>>>>>>>>>>   &END DFT
>>>>>>>>>>   &SUBSYS
>>>>>>>>>>     &CELL
>>>>>>>>>>    ABC [angstrom]   30 30 30
>>>>>>>>>>     &END CELL
>>>>>>>>>>     &TOPOLOGY
>>>>>>>>>>       COORD_FILE_NAME water.xyz
>>>>>>>>>>       COORDINATE xyz
>>>>>>>>>>     &END
>>
>>>>>>>>>>     &KIND O
>>>>>>>>>>       BASIS_SET TZV2P-MOLOPT-GTH
>>>>>>>>>>       POTENTIAL GTH-PBE-q6
>>>>>>>>>>     &END KIND
>>>>>>>>>>     &KIND H
>>>>>>>>>>       BASIS_SET TZV2P-MOLOPT-GTH
>>>>>>>>>>       POTENTIAL GTH-PBE-q1
>>>>>>>>>>     &END KIND
>>
>>>>>>>>>>   &END SUBSYS
>>
>>>>>>>>>> &END FORCE_EVAL
>>
>>>>>>>>>> &GLOBAL
>>>>>>>>>>   PRINT_LEVEL LOW
>>>>>>>>>>   PROJECT H2O
>>>>>>>>>>   RUN_TYPE ENERGY
>>>>>>>>>>   WALLTIME 285000
>>>>>>>>>> &END GLOBAL
>>
>>>>>>>>>> &MOTION
>>>>>>>>>>     &GEO_OPT
>>>>>>>>>>      MAX_ITER 5000
>>>>>>>>>>      MAX_FORCE 0.00010
>>>>>>>>>>      OPTIMIZER BFGS
>>>>>>>>>>     &BFGS
>>>>>>>>>>     &END
>>>>>>>>>>   &END
>>>>>>>>>> &END
>>
>>>>>>>>>> ------------------------------------------------
>>>>>>>>>>  the water.xyz file is
>>>>>>>>>> ----------------------------------------
>>>>>>>>>>        3
>>>>>>>>>>  water
>>>>>>>>>>   H         8.4351137971        5.1620318678       12.1510352783
>>>>>>>>>>   H         8.5534822942        3.6347695084       12.1803075532
>>>>>>>>>>   O         7.9011326518        4.3528437024       12.1874564118
>>>>>>>>>> ------------------------------------------------
>>
>>>>>>>>>> Regards,
>>>>>>>>>> Manh
>>
>>>>>>>>> --
>>>>>>>>> You received this message because you are subscribed to the Google Groups "cp2k" group.
>>>>>>>>> To post to this group, send email to cp... at googlegroups.com.
>>>>>>>>> To unsubscribe from this group, send email to cp2k+uns... at googlegroups.com.
>>>>>>>>> For more options, visit this group athttp://groups.google.com/group/cp2k?hl=en.
>>
>>>>>>> --
>>>>>>> You received this message because you are subscribed to the Google Groups "cp2k" group.
>>>>>>> To post to this group, send email to cp... at googlegroups.com.
>>>>>>> To unsubscribe from this group, send email to cp2k+uns... at googlegroups.com.
>>>>>>> For more options, visit this group athttp://groups.google.com/group/cp2k?hl=en.
>>
>>>>> --
>>>>> You received this message because you are subscribed to the Google Groups "cp2k" group.
>>>>> To post to this group, send email to cp... at googlegroups.com.
>>>>> To unsubscribe from this group, send email to cp2k+uns... at googlegroups.com.
>>>>> For more options, visit this group athttp://groups.google.com/group/cp2k?hl=en.
>>
>>> --
>>> You received this message because you are subscribed to the Google Groups "cp2k" group.
>>> To post
>>
>> ...
>>