[Q-e-developers] [Pw_forum] two program logic questions in TDDFPT

[Paolo Giannozzi]

Thank you for your remarks. I don't know about the first point
you raise. I notice however a dangerous construct:
   CALL zcopy(size_evc,evc1_new(:,:,:,1),1,evc1_new(:,:,:,2),1)
Old-style fortran routines should be called with old-style array
syntax, i.e.
   CALL zcopy(size_evc,evc1_new(1,1,1,1),1,evc1_new(1,1,1,2),1)
About the second point you raise, you are right: either the code 
is stopped when beta is very small, or a warning is issued and the
code goes on, but then beta and gamma must be properly calculated.
(CC: to the q-e-developers at qe-forge.org mailing list)

Paolo

On Mon, 2012-04-16 at 15:19 +0800, kuilin lu wrote:
> Dear QE users and developers:
>    I have two program logic questions in TDDFPT.  The questions are
> about lr_lanczos.f90.
> 
> -----------------------
> place 1:
> /-----------------------------------------------------------------------------------------------------------\
> place 1
> 273     IF(.not.ltammd) THEN
> 274     !
> 275     IF ( mod(LR_iteration,2)==0 ) THEN
> 276        CALL
> lr_apply_liouvillian(evc1(1,1,1,1),evc1_new(1,1,1,1),sevc1_new(1,1,1,1),.false.)
> 277        CALL
> lr_apply_liouvillian(evc1(1,1,1,2),evc1_new(1,1,1,2),sevc1_new(1,1,1,2),.true.)
> 278     ELSE
> 279        CALL
> lr_apply_liouvillian(evc1(1,1,1,1),evc1_new(1,1,1,1),sevc1_new(1,1,1,1),.true.)
> 280        CALL
> lr_apply_liouvillian(evc1(1,1,1,2),evc1_new(1,1,1,2),sevc1_new(1,1,1,2),.false.)
> 281     ENDIF
> 282     !
> 283     ELSE
> 284        CALL
> lr_apply_liouvillian(evc1(1,1,1,1),evc1_new(1,1,1,1),sevc1_new(1,1,1,1),.true.)
> 285        CALL zcopy(size_evc,evc1(:,:,:,1),1,evc1(:,:,:,2),1)
> !evc1(,1) = evc1(,2)
> 286        CALL
> zcopy(size_evc,evc1_new(:,:,:,1),1,evc1_new(:,:,:,2),1) !evc1_new(,1)
> = evc1_new(,2)
> 287        CALL
> zcopy(size_evc,evc1_new(:,:,:,1),1,evc1_new(:,:,:,2),1) !evc1_new(,1)
> = evc1_new(,2)
> 288
> 289     ENDIF
> \-----------------------------------------------------------------------------------------------------------/
> at line 287, it is
>   CALL zcopy(size_evc,evc1_new(:,:,:,1),1,evc1_new(:,:,:,2),1)
> !evc1_new(,1) = evc1_new(,2)
> should here be the following?
>   CALL zcopy(size_evc,sevc1_new(:,:,:,1),1,sevc1_new(:,:,:,2),1)
> !sevc1_new(,1) = sevc1_new(,2)
> 
>   since it is a waste to copy evc1_new(:,:,:,1) two times and missing
> setting value of sevc1_new(:,:,:,2)
> \
> 
> -----------------------
> place 2:
> /-----------------------------------------------------------------------------------------------------------\
> place 2
> 217     IF ( abs(beta)<1.0d-12 ) THEN
> 218        !
> 219        WRITE(stdout,'(5x,"lr_lanczos: Left and right Lanczos
> vectors are orthogonal, this is a violation of oblique projection")')
> 220        !
> 221     ELSEIF ( beta<0.0d0 ) THEN
> 222        !
> 223        beta=sqrt(-beta)
> 224        gamma=-beta
> 225        !
> 226     ELSEIF ( beta>0.0d0 ) THEN
> 227        !
> 228        beta=sqrt(beta)
> 229        gamma=beta
> 230        !
> 231     ENDIF
> \-----------------------------------------------------------------------------------------------------------
> /
> at line 217
> 
>   if beta is by accidently small, it just print a warning
> “lr_lanczos: Left and right Lanczos vectors are orthogonal, this is a
> violation of oblique projection”
> then exit the IF structure. So variable gamma is actually not
> computed, it will use the value of last step.  It might cause NaN
> error later since beta and gamma are not same in this step and
> different from lanczos algorithm
> (http://www.cs.utk.edu/~dongarra/etemplates/node245.html) which is
> mentioned in the code.
> 
>  Should here be the following?
> 
> !------
>     IF ( abs(beta)<1.0d-12 ) THEN
>        !
>        WRITE(stdout,'(5x,"lr_lanczos: Left and right Lanczos vectors
> are orthogonal, this is a violation of oblique projection")')
>        !
>     ENDIF
>     IF ( beta<0.0d0 ) THEN
>        !
>        beta=sqrt(-beta)
>        gamma=-beta
>        !
>     ELSEIF ( beta>0.0d0 ) THEN
>        !
>        beta=sqrt(beta)
>        gamma=beta
>        !
>     ENDIF
> !------
>     I met this problem when I wrote my program based on Siesta (which
> is a DFT software),  finding that  variable beta can be small, and get
> NaN later.  My  program could get reasonable result after using my
> changed logic. The test result is in accessory of this mail.  So I
> doubt there is some problem here. If the original logic is right,  is
> there some reason to do so? Thanks in advance.
>  \
> 
> Best Wishes
> Kuilin Lu
> _______________________________________________
> Pw_forum mailing list
> Pw_forum at pwscf.org
> http://www.democritos.it/mailman/listinfo/pw_forum
-- 
Paolo Giannozzi, IOM-Democritos and University of Udine, Italy