dephasing in yambo_nl

[sangallidavide]

The correct expression should be
sum_n |u_n(t)><u_n(t) | - |u_n(0)><u_n(t) | [1]
An alternative expression is
sum_n |u_n(t)><u_n(t) | - |u_n(0)><u_n(0) | [2]

I had a look to the way it is coded into yambo_nl and it looks something in between the two.
Indeed the collision integral is define with 3 band indexes I( i1, i2, ib), and the ib index is used to apply the dephasing only to band "ib".

This means it uses |u_n(t)><u_n(t) | u_b(t) > = |u_n(t)> delta_{n,b}
As a consequence (first part of the dephasing)
sum_{jn}<u0_i |u_n(t)><u_n(t) | u0_j> <u0_j| u_b(t) > = sum_{j}<u0_i |u_b(t)><u_b(t) | u0_j> <u0_j| u_b(t) > ( = sum_{j}<u0_i |u_b(t)> )
So far all good.

However the same trick is used for the term |u0_n><u0_n |, i.e. implying that |u0_n><u0_n | u_b(t) > = |u0_n> delta_{n,b}
This is not true for |u0_n><u0_n | [2], but it would for |u0_n><u_n(t) | [1]

Then in yambo_nl there is implemented
-sum_{jn}<u0_i |u0_n><u0_n | u0_j> <u0_j| u_b(t) > = -sum_{j}<u0_i |u0_b><u0_b | u0_j> <u0_j| u_b(t) > = -\delta_{i,b} <u0_b| u_b(t) > [2nl]
in place of
-sum_{jn}<u0_i |u0_n><u0_n | u0_j> <u0_j| u_b(t) > = <u0_i| u_b(t) > [2theo]

the expression with |u0_n><u_n(t) | would be
-sum_{jn}<u0_i |u0_n><u_n(t) | u0_j> <u0_j| u_b(t) > = -sum_{j}<u0_i |u0_b><u_b(t) | u0_j> <u0_j| u_b(t) > = -\delta_{i,b} [1code]

I have the impression that the one implemented [2nl] is quite similar to [1code]

(duplicate of https://github.com/yambo-code/yambo-devel/issues/674)

[sangallidavide]

Comment from Ignacio:

This may be very basic, or it may add something, I don't know.

For me, the issue is in the summation over bands that is in the definition of the dephasing term.

1. If you define the dephasing term as in the paper:
   sum_n ( |u_n(t)><u_n(t)| - |u_n(0)><u_n(0)| ),
   then you cannot get to what is coded, for the reasons Davide explained.
   1.a. The second part of the dephasing in KS space is the identity in the occupied block and zero everywhere else (Gamma0_{i,j} =sum_{n} delta_{i,n} delta_{j,n} with n running over occupied bands). [2theo]
   1.b. Meanwhile, what is coded has an additional delta_{n,b}, i.e., just 1 in i=j=b when the EOM is for band |u_b(t)>. [2nl]
   1.c Then, we start thinking that the last bra should actually be time-dependent, to force that delta appearing, but that also brings about an extra <u_b(t) | u_j(0)> for all j at row i=b of matrix Gamma0. This turns into delta_{i,b} upon multiplying with the corresponding vector, as Davide says in [1code].

Nothing new so far.

2. Alternatively, you can define the dephasing as in the paper but without the summation over occupied bands. Instead, you'd define it for the particular band that you are integrating in the EOM, i.e., in the EOM for band |u_b(t)>, you'd have:
   |u_b(t)><u_b(t)| - |u_b(0)><u_b(0)| [2 no sum] = [2nl].
   2.a. I believe this, in KS space, is exactly what is coded. The additional delta_{n,b} is implied when deciding not to use the summation in the definition. Also, it has the property of acting as a restoring force when |u_b(t)> departs from |u_b(0)>.

I ended up in this line of thinking because I actually tried [1code] and it didn't give good results, importantly, not conserving the norm of the vectors.

Now, the reason why my take may be very basic is that I ignore the derivation of this dephasing term, and maybe the summation over bands arises from such derivation. Or maybe it is just not correct that the Hamiltonian is the same for every band at a given kpt, but the dephasing is different.

[sangallidavide]

Comment from myself

Addendum on the damping in yambo_nl.

At present the damping gives asymmetric excited state density between the valence and the conduction channel.
(See below.)

To fix this I think it is enough to damp in the same way the unoccupied and the empty states.
This can be done in the "old damping" by removing these lines

   !                                                                    
   ! Remove the equilibrium contribution                                
   !                                                                    
!$omp parallel do default(shared), private(ib)                          
   do ib=NL_bands(1),E%nbf(i_sp_pol)                                    
     I_relax(ib,ib,ib,ik_mem,i_sp_pol)=I_relax(ib,ib,ib,ik_mem,i_sp_pol)+cI*2._SP/RAD_LifeTime
   enddo                                                                
!$omp end parallel do                

and similarly in the new damping by avoiding the "cycle" command

     do ib=1,E%nbf(i_sp_pol)                                            
       if(ib==i1) cycle                                                 
       I_relax(i1,i1,ib,ik_mem,i_sp_pol) = -cI*Gamma_ij                 
     enddo                  

After these changes I_relax must then be applied to psi-psi_0

---

Asymmetric excited density example:

• Case with damping
  violet= excited density in valence
  green= excited density in conduction
  

#
#    Time[fs]           d(N_hole-N_elec)   dN_hole            dN_elec
#
      0.00000000         0.00000000         0.00000000         0.00000000
     0.10000000         0.14176308E-11     0.32022329E-10     0.30604698E-10
     0.20000000         0.33824048E-11     0.41433135E-10     0.38050730E-10
     0.30000000         0.54020576E-11     0.40367099E-10     0.34965041E-10
     0.40000000         0.67837464E-11     0.29289293E-10     0.22505547E-10
     0.50000000         0.74766757E-11     0.24822366E-10     0.17345691E-10
  ...

      6.30000000        0.157054311E-6     0.516333760E-6     0.359279449E-6
      6.40000000        0.168251346E-6     0.485183143E-6     0.316931798E-6
      6.50000000        0.174301182E-6     0.406980470E-6     0.232679289E-6
      6.60000000        0.177467175E-6     0.422291212E-6     0.244824037E-6
      6.70000000        0.184106704E-6     0.520605685E-6     0.336498981E-6
      6.80000000        0.194450167E-6     0.556460077E-6     0.362009910E-6
      6.90000000        0.202171760E-6     0.481935406E-6     0.279763647E-6
      7.00000000        0.204394094E-6     0.421227031E-6     0.216832938E-6
      7.10000000        0.206058105E-6     0.470370302E-6     0.264312197E-6

• Case without damping
  

#
#    Time[fs]           d(N_hole-N_elec)   dN_hole            dN_elec
#
      0.00000000         0.00000000         0.00000000         0.00000000
     0.10000000         0.10080449E-15     0.10547119E-14     0.95390739E-15
     0.20000000         0.56955534E-15     0.24147351E-13     0.23577795E-13
     0.30000000        -0.74781671E-15     0.11074475E-12     0.11149256E-12
     0.40000000        -0.18246052E-14     0.26151303E-12     0.26333764E-12
...
      14.300000        -0.11504262E-12     0.25065744E-08     0.25066895E-08
      14.400000        -0.11543977E-12     0.18098391E-08     0.18099546E-08
      14.500000        -0.11544604E-12     0.15258592E-08     0.15259746E-08
      14.600000        -0.11600221E-12     0.18552923E-08     0.18554083E-08
      14.700000        -0.11579392E-12     0.25238999E-08     0.25240157E-08
      14.800000        -0.11723663E-12     0.30024183E-08     0.30025355E-08
      14.900000        -0.11848214E-12     0.29152621E-08     0.29153806E-08
      15.000000        -0.11968546E-12     0.23239206E-08     0.23240403E-08
      15.100000        -0.11990438E-12     0.16763167E-08     0.16764366E-08