Siesta

Code review comment for lp:~nickpapior/siesta/4.1-md-supercell

4.1-md-supercell
Merge into rel-4.1

Revision history for this message

Alberto Garcia (albertog) wrote on 2018-09-03:

> On 3 Sep 2018, at 10:21, Nick Papior <email address hidden> wrote:
>
> 1) Is it a calculation with folding that you want to initialize from? I.e. the nsc_read = (1,1,1) in this case?

Yes. Assume a user has run a calculation without k-points (and without an auxcell), and then wants to use the DM to initialize a
k-point calculation. This is what happens in the si64 example:

—————————————
Initializing Density Matrix...

Attempting to read DM from file... Succeeded...
DM from file:
<dSpData2D:IO-DM: si64.DM
  <sparsity:IO-DM: si64.DM
    nrows_g=256 nrows=136 sparsity=.4275 nnzs=28016, refcount: 1>
  <dData2D:(new from dSpData2D) n=28016 m=1, refcount: 1>
refcount: 1>
Note: For starting DM, Qtot, Tr[D*S] = 256.00000000 229.72632220
—————————————

This is a calculation with “normal” folding, i.e. only edge image interactions, but no direct diagonal images, so that S(i,i) = 1.0.

The missing norm is coming (or not coming…) from the zeros in the supercell DM.

Later on:

———————————
        iscf Eharris(eV) E_KS(eV) FreeEng(eV) dDmax Ef(eV) dHmax(eV)
   scf: 1 10293.156932 -1324.422347 -1324.455994 3.337586 -3.198437142.619752
timer: Routine,Calls,Time,% = IterSCF 1 4.132 80.76
   scf: 2 35111.572145 3756.559710 3756.553850 4.416191 -2.224851288.324694
   scf: 3 -1549.974918 444.670819 444.641579 3.270160 -2.128185164.985102
   scf: 4 -5371.294656 -2706.434665 -2706.485243 1.485187 -2.565191125.811367
   scf: 5 18837.075132 -728.224618 -728.241614 4.270983 -3.174556153.635637
   scf: 6 -8033.559239 -5034.599795 -5034.645848 2.447968 -2.942511 88.795877
   scf: 7 -7165.186747 -6839.938924 -6839.938924 2.852955 -2.675640 4.831900
   scf: 8 -6836.790121 -6839.363374 -6839.363374 0.078137 -2.668903 4.867
…..
————————

Note that, for k=0, even in extreme folding cases, the DM has a very simple structure, and is NOT folded. Only H and S are actually folded. ALL the image elements of the DM are the same as those in the unit-cell square matrix. This is the code in diagg:

                do j = 1,numd(io)
                  ind = listdptr(io) + j
                  jo = listd(ind)
                  jo = MODP(jo,nuotot) ! To allow auxiliary supercells
                  Dnew(ind,ispin) = Dnew(ind,ispin) + qei*paux(jo)
                  Enew(ind,ispin) = Enew(ind,ispin) + eei*paux(jo)
                enddo

So all the entries (arbitrary “ind”) have the same value as those in the “fold-to” corner of the rectangular DM.

To get back the missing norm approximately in the nsc_read= (1,1,1) case, one then could do:

foreach non-zero extended S entry “ind_new” in the new, full-nsc supercell:
find the io,jo entry in the read DM, and copy the entry to the new DM.

It might be easier to reverse the logic and just “un-fold” the read DM, with a loop very similar to the above.

Of course, it is simpler to run the k=0 calculation with “force-auxcell”, but this change looks simple, and will add functionality.

The handling of intermediate cases, when nsc_read is a “subset” of nsc, but not (1,1,1) is more complicated (assuming k-points have been used).
In this case the DM entries contain different phases for each k-point’s contribution.

> 2) I agree, this would be the best, to handle the overlap matrix, however, in cases where one wan't to create an initial DM from previously (separate) calculated DM it becomes too elaborate. In the end I think this is mainly a good idea from a perspective of doing manual interpolation and/or manual initialization.
>

The issues discussed above can be present even in a normal MD calculation.

> As for the nsc_read == 1, I can successfully run the script I send (privately) with nsc_read == 1, without problems?
> Secondly, where do you want the extra clarification to be put, in the dev-docs, or the source?
> Could you check out the new commit?

I did not get the script… I am also experimenting with a new idiom in FORD documentation for the internals of routines. More later.

> --
> https://code.launchpad.net/~nickpapior/siesta/4.1-md-supercell/+merge/354028
> You are requested to review the proposed merge of lp:~nickpapior/siesta/4.1-md-supercell into lp:siesta/4.1.

> On 3 Sep 2018, at 10:21, Nick Papior <nickpapior@gmail.com> wrote:
> 
> 1) Is it a calculation with folding that you want to initialize from? I.e. the nsc_read = (1,1,1) in this case?

Yes. Assume a user has run a calculation without k-points (and without an auxcell), and then wants to use the DM to initialize a
k-point calculation. This is what happens in the si64 example:

—————————————
Initializing Density Matrix...

Attempting to read DM from file... Succeeded...
DM from file:
<dSpData2D:IO-DM: si64.DM
  <sparsity:IO-DM: si64.DM
    nrows_g=256 nrows=136 sparsity=.4275 nnzs=28016, refcount: 1>
  <dData2D:(new from dSpData2D) n=28016 m=1, refcount: 1>
refcount: 1>
Note: For starting DM, Qtot, Tr[D*S] =        256.00000000        229.72632220
—————————————

This is a calculation with “normal” folding, i.e. only edge image interactions, but no direct diagonal images, so that S(i,i) = 1.0.

The missing norm is coming (or not coming…) from the zeros in the supercell DM.

Later on:

———————————
        iscf     Eharris(eV)        E_KS(eV)     FreeEng(eV)     dDmax    Ef(eV) dHmax(eV)
   scf:    1    10293.156932    -1324.422347    -1324.455994  3.337586 -3.198437142.619752
timer: Routine,Calls,Time,% = IterSCF        1       4.132  80.76
   scf:    2    35111.572145     3756.559710     3756.553850  4.416191 -2.224851288.324694
   scf:    3    -1549.974918      444.670819      444.641579  3.270160 -2.128185164.985102
   scf:    4    -5371.294656    -2706.434665    -2706.485243  1.485187 -2.565191125.811367
   scf:    5    18837.075132     -728.224618     -728.241614  4.270983 -3.174556153.635637
   scf:    6    -8033.559239    -5034.599795    -5034.645848  2.447968 -2.942511 88.795877
   scf:    7    -7165.186747    -6839.938924    -6839.938924  2.852955 -2.675640  4.831900
   scf:    8    -6836.790121    -6839.363374    -6839.363374  0.078137 -2.668903  4.867
…..
————————

Note that,  for k=0, even in extreme folding cases, the DM has a very simple structure, and is NOT folded. Only H and S are actually folded. ALL the image elements of the DM are the same as those in the unit-cell square matrix. This is the code in diagg:

do j = 1,numd(io)
                  ind = listdptr(io) + j
                  jo  = listd(ind)
                  jo = MODP(jo,nuotot) ! To allow auxiliary supercells                         
                  Dnew(ind,ispin) = Dnew(ind,ispin) + qei*paux(jo)
                  Enew(ind,ispin) = Enew(ind,ispin) + eei*paux(jo)
                enddo

So all the entries (arbitrary “ind”) have the same value as those in the “fold-to” corner of the rectangular DM.

To get back the missing norm approximately in the nsc_read= (1,1,1) case, one then could do:

foreach non-zero extended S entry “ind_new” in the new, full-nsc supercell:
       find the io,jo entry in the read DM, and copy the entry to the new DM.

It might be easier to reverse the logic and just “un-fold” the read DM, with a loop very similar to the above.

Of course, it is simpler to run the k=0 calculation with “force-auxcell”, but this change looks simple, and will add functionality.
      
The handling of intermediate cases, when nsc_read is a “subset” of nsc, but not (1,1,1) is more complicated (assuming k-points have been used).
In this case the DM entries contain different phases for each k-point’s contribution.

The issues discussed above can be present even in a normal MD calculation.

I did not get the script… I am also experimenting with a new idiom in FORD documentation for the internals of routines. More later.

> -- 
> https://code.launchpad.net/~nickpapior/siesta/4.1-md-supercell/+merge/354028
> You are requested to review the proposed merge of lp:~nickpapior/siesta/4.1-md-supercell into lp:siesta/4.1.

« Back to merge proposal