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 -- ---------------------------------------------------------------------- Kuilin Lu, Ph.D student College of Materials Science and Engineering Hunan University Changsha 410082 China E-Mail: lukuilin at gmail.com ---------------------------------------------------------------------- -------------- next part -------------- A non-text attachment was scrubbed... Name: figure.png Type: image/png Size: 5702 bytes Desc: not available Url : http://www.democritos.it/pipermail/pw_forum/attachments/20120416/5c15d0fa/attachment.png