From 0d4d58b532b7c7d30acca5f50772c4946b9d1596 Mon Sep 17 00:00:00 2001 From: obm Date: Sat, 3 Apr 2010 01:27:40 +0000 Subject: [PATCH] Bugfixes for AIX machines and intel compilers git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@6593 c92efa57-630b-4861-b058-cf58834340f0 --- TDDFPT/Doc/TODO | 9 ++ TDDFPT/src/lr_apply_liouvillian.f90 | 3 +- TDDFPT/src/lr_dvpsi_e.f90 | 91 +++++++++++++-- TDDFPT/src/lr_lanczos.f90 | 85 +++++++++++--- TDDFPT/src/lr_main.f90 | 33 ++++-- TDDFPT/src/lr_ortho.f90 | 169 +++++++++++----------------- TDDFPT/src/lr_read_wf.f90 | 25 +++- TDDFPT/src/lr_variables.f90 | 130 +++++++++++++++++++++ 8 files changed, 399 insertions(+), 146 deletions(-) create mode 100644 TDDFPT/Doc/TODO diff --git a/TDDFPT/Doc/TODO b/TDDFPT/Doc/TODO new file mode 100644 index 000000000..935ab9d54 --- /dev/null +++ b/TDDFPT/Doc/TODO @@ -0,0 +1,9 @@ +1) Possible optimizations +1.1 revc0 each other band is not used in gamma point case +2) Feature enhancements +2.1 Exact exchange +2.2 PAW +2.3 Spin +3) Known bugs +3.1 Precision is not achieved in xl compiler + diff --git a/TDDFPT/src/lr_apply_liouvillian.f90 b/TDDFPT/src/lr_apply_liouvillian.f90 index fc04d7867..6c588389d 100755 --- a/TDDFPT/src/lr_apply_liouvillian.f90 +++ b/TDDFPT/src/lr_apply_liouvillian.f90 @@ -45,7 +45,8 @@ subroutine lr_apply_liouvillian( evc1, evc1_new, sevc1_new, interaction ) ! implicit none ! - complex(kind=dp) :: evc1(npwx,nbnd,nks), evc1_new(npwx,nbnd,nks), sevc1_new(npwx,nbnd,nks) + complex(kind=dp),intent(in) :: evc1(npwx,nbnd,nks) + complex(kind=dp),intent(out) :: evc1_new(npwx,nbnd,nks), sevc1_new(npwx,nbnd,nks) logical, intent(in) :: interaction ! ! Local variables diff --git a/TDDFPT/src/lr_dvpsi_e.f90 b/TDDFPT/src/lr_dvpsi_e.f90 index 1620a8b24..788a33566 100755 --- a/TDDFPT/src/lr_dvpsi_e.f90 +++ b/TDDFPT/src/lr_dvpsi_e.f90 @@ -52,12 +52,15 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi) USE lr_variables, ONLY : lr_verbosity, evc0 USE io_global, ONLY : stdout - +!DEBUG + USE lr_variables, ONLY: check_all_bands_gamma, check_density_gamma,check_vector_gamma,check_vector_f + ! + ! implicit none ! integer, intent(IN) :: ipol, ik ! - complex(kind=dp) :: dvpsi(npwx,nbnd) + complex(kind=dp),intent(out) :: dvpsi(npwx,nbnd) real(kind=dp) :: atnorm ! @@ -112,9 +115,7 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi) !OBM!!!! eprec is also calculated on the fly for each k point allocate(eprec(nbnd)) !OBM!!!! - evc(:,:)=evc0(:,:,ik) - !print *, "evc", evc(1:3,1) - !print *,"npw_k(",ik,")=",npw_k(ik) + evc(:,:)=evc0(:,:,ik) if (nkb > 0) then allocate (dvkb (npwx, nkb), dvkb1(npwx, nkb)) dvkb (:,:) = (0.d0, 0.d0) @@ -125,7 +126,12 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi) gk (2, ig) = (xk (2, ik) + g (2, igk (ig) ) ) * tpiba gk (3, ig) = (xk (3, ik) + g (3, igk (ig) ) ) * tpiba g2kin (ig) = gk (1, ig) **2 + gk (2, ig) **2 + gk (3, ig) **2 - enddo + enddo + if (lr_verbosity > 10 ) then + write(stdout,'("lr_dvpsi_e g2kin:",F15.8)') SUM(g2kin(:)) + endif + + ! ! this is the kinetic contribution to [H,x]: -2i (k+G)_ipol * psi ! @@ -153,7 +159,12 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi) ! ! ! enddo ! print *, "lr_dvpsi_e d0psi kinetic contribution", obm_debug - + if (lr_verbosity > 10 ) then + write(stdout,'("lr_dvpsi_e d0psi kinetic contribution:")') + do ibnd=1,nbnd + call check_vector_gamma(d0psi(:,ibnd)) + enddo + endif !!obm_debug @@ -168,7 +179,38 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi) ! and this is the contribution from nonlocal pseudopotentials ! call gen_us_dj (ik, dvkb) - call gen_us_dy (ik, at (1, ipol), dvkb1) + call gen_us_dy (ik, at (1, ipol), dvkb1) + if (lr_verbosity > 10 ) then + write(stdout,'("lr_dvpsi_e dvkb:")') + jkb=0 + do nt = 1, ntyp + do na = 1, nat + if (nt == ityp (na)) then + do ikb = 1, nh (nt) + jkb = jkb + 1 + call check_vector_f(dvkb(:,jkb)) + enddo + endif + enddo + enddo + endif + if (lr_verbosity > 10 ) then + write(stdout,'("lr_dvpsi_e dvkb1:")') + jkb=0 + do nt = 1, ntyp + do na = 1, nat + if (nt == ityp (na)) then + do ikb = 1, nh (nt) + jkb = jkb + 1 + call check_vector_f(dvkb1(:,jkb)) + enddo + endif + enddo + enddo + endif + + + do ig = 1, npw_k(ik) if (g2kin (ig) < 1.0d-10) then gk (1, ig) = 0.d0 @@ -201,9 +243,23 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi) deallocate (gk) !OBM!!!be careful, from bwalker, why?!!!!! work(:,:)=(0.0d0,1.0d0)*work(:,:) - !OBM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !print *, "nonlocal contribution"!, d0psi(1:3,1) - !CALL lr_normalise( d0psi(:,:), anorm) + !OBM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if (lr_verbosity > 10 ) then + write(stdout,'("lr_dvpsi_e non-local contribution:")') + jkb=0 + do nt = 1, ntyp + do na = 1, nat + if (nt == ityp (na)) then + do ikb = 1, nh (nt) + jkb = jkb + 1 + call check_vector_gamma(work(:,jkb)) + enddo + endif + enddo + enddo + endif + + if(gamma_only) then call lr_dvpsi_e_gamma() else if (noncolin) then @@ -219,7 +275,18 @@ subroutine lr_dvpsi_e(ik,ipol,dvpsi) ! ! ! enddo ! print *, "lr_dvpsi_e norm of dvpsi being returned=", obm_debug - ! + ! + if (lr_verbosity > 10 ) then + write(stdout,'("lr_dvpsi_e d0psi after lr_dvpsi_e case specific calc:")') + do ibnd=1,nbnd + call check_vector_gamma(d0psi(:,ibnd)) + enddo + write(stdout,'("lr_dvpsi_e dvpsii after lr_dvpsi_e case specific calc:")') + do ibnd=1,nbnd + call check_vector_gamma(dvpsi(:,ibnd)) + enddo + endif + !!obm_debug IF (nkb > 0) THEN diff --git a/TDDFPT/src/lr_lanczos.f90 b/TDDFPT/src/lr_lanczos.f90 index 4c3334fd1..ca0d0bd41 100755 --- a/TDDFPT/src/lr_lanczos.f90 +++ b/TDDFPT/src/lr_lanczos.f90 @@ -38,7 +38,8 @@ contains v_loc_psir, s_psir_gamma, igk_k,npw_k, real_space_debug USE lr_variables, ONLY : lr_verbosity, charge_response use charg_resp, only : w_T_beta_store,w_T,lr_calc_F - + !Debugging + USE lr_variables, ONLY: check_all_bands_gamma, check_density_gamma,check_vector_gamma ! implicit none ! @@ -51,7 +52,7 @@ contains ! ! Local variables ! - real(kind=dp) :: alpha, beta, gamma + real(kind=dp) :: alpha, beta, gamma, temp ! complex(kind=dp) :: zeta ! @@ -61,14 +62,14 @@ contains ! If (lr_verbosity > 5) THEN WRITE(stdout,'("")') + endif + If (lr_verbosity > 10) THEN print *, "Real space = ", real_space print *, "Real space debug ", real_space_debug print *, "TQR = ", tqr endif ! call start_clock('one_step') - - !print *, "norm of d0psi", lr_dot(d0psi(1,1,1,1),d0psi(1,1,1,1)) pol_index=1 if ( n_ipol /= 1 ) pol_index=LR_polarization ! @@ -117,7 +118,28 @@ contains evc1_new(:,:,:,2)=evc1_new(:,:,:,1) sevc1_new(:,:,:,2)=sevc1_new(:,:,:,1) endif - LR_iteration=1 + LR_iteration=1 + !Debugging + if (lr_verbosity >10) then + !write(stdout,'("evc1_new(1), first step")') + !do ibnd=1,nbnd + ! call check_vector_gamma(evc1_new(:,ibnd,1,1)) + !enddo + !write(stdout,'("sevc1_new(1), first step")') + !do ibnd=1,nbnd + ! call check_vector_gamma(sevc1_new(:,ibnd,1,1)) + !enddo + !write(stdout,'("evc1_new(2), first step")') + !do ibnd=1,nbnd + ! call check_vector_gamma(evc1_new(:,ibnd,1,2)) + !enddo + !write(stdout,'("sevc1_new(2), first step")') + !do ibnd=1,nbnd + ! call check_vector_gamma(sevc1_new(:,ibnd,1,2)) + !enddo + temp=dble(lr_dot(evc1_new(1,1,1,1),sevc1_new(1,1,1,2))) + write(stdout,'(" first step",E15.8)') temp + endif endif ! ! The lanczos algorithm starts here @@ -125,12 +147,6 @@ contains ! ! Left and right vectors are orthogonalised wrto ground state wf - !print *, "norm of evc1,1 before lr_ortho", lr_dot(evc1_new(1,1,1,1),evc1_new(1,1,1,1)) - !print *, "norm of evc1,2 before lr_ortho", lr_dot(evc1_new(1,1,1,2),evc1_new(1,1,1,2)) - !print *, "norm of evc0 before lr_ortho", lr_dot(evc0(1,1,1),evc0(1,1,1)) - !print *, "norm of sevc0 before lr_ortho", lr_dot(sevc0(1,1,1),sevc0(1,1,1)) - !print *,"nks=", nks - !OBM: Notice that here "orthogonalization" is not strictly the true word, as the norm of the vectors change !This is due to how the uspp scheme is implemented, the beta are evc1_new(left).sevc1_new(right), that is, !a mixing of two vectors, thus the resultant vector from belov should be devoid from S, which affects the norm @@ -139,10 +155,12 @@ contains call lr_ortho(evc1_new(:,:,ik,1), evc0(:,:,ik), ik, ik, sevc0(:,:,ik),.true.) call lr_ortho(evc1_new(:,:,ik,2), evc0(:,:,ik), ik, ik, sevc0(:,:,ik),.true.) enddo - - !print *, "norm of evc1,1 after lr_ortho", lr_dot(evc1_new(1,1,1,1),evc1_new(1,1,1,1)) - !print *, "norm of evc1,2 after lr_ortho", lr_dot(evc1_new(1,1,1,1),evc1_new(1,1,1,2)) - + if (lr_verbosity >10) then + temp=dble(lr_dot(evc1_new(1,1,1,1),sevc0(1,1,1))) + write(stdout,'("",E15.8)') temp + temp=dble(lr_dot(evc1_new(1,1,1,2),sevc0(1,1,1))) + write(stdout,'("",E15.8)') temp + endif ! ! By construction =0 should be 0, forcing this both conserves resources and increases stability @@ -215,7 +233,6 @@ contains end if !print *, "norm of sevc1,1 after spsi", lr_dot(sevc1_new(1,1,1,1),sevc1_new(1,1,1,1)) !print *, "norm of sevc1,2 after spsi", lr_dot(sevc1_new(1,1,1,1),sevc1_new(1,1,1,2)) - !Resume the LR ! ! Orthogonality requirement as proposed by Y. Saad beta=sqrt(|qdash.pdash|) gamma=sign(qdash.pdash)*beta @@ -276,6 +293,18 @@ contains evc1_new(:,:,:,:)=(0.0d0,0.0d0) sevc1_new(:,:,:,:)=(0.0d0,0.0d0) ! + if (lr_verbosity >10) then + write(stdout,'("evc1(1), rotate")') + do ibnd=1,nbnd + call check_vector_gamma(evc1(:,ibnd,1,1)) + enddo + write(stdout,'("evc1(2), rotate")') + do ibnd=1,nbnd + call check_vector_gamma(evc1(:,ibnd,1,2)) + enddo + endif + + ! ! if(.not.ltammd) then @@ -294,7 +323,17 @@ contains call zcopy(size_evc,evc1_new(:,:,:,1),1,evc1_new(:,:,:,2),1) !evc1_new(,1) = evc1_new(,2) call zcopy(size_evc,evc1_new(:,:,:,1),1,evc1_new(:,:,:,2),1) !evc1_new(,1) = evc1_new(,2) - end if + end if + if (lr_verbosity >10) then + write(stdout,'("evc1(1), apply L")') + do ibnd=1,nbnd + call check_vector_gamma(evc1_new(:,ibnd,1,1)) + enddo + write(stdout,'("evc1(2), apply L")') + do ibnd=1,nbnd + call check_vector_gamma(evc1_new(:,ibnd,1,2)) + enddo + endif ! ! qdash(i+1)=f(q(i))-gamma*q(i-1) ! pdash(i+1)=f(p(i))-beta*p(i-1) @@ -302,7 +341,17 @@ contains ! !OBM BLAS call zaxpy(size_evc,-cmplx(gamma,0.0d0,kind=dp),evc1_old(:,:,:,1),1,evc1_new(:,:,:,1),1) - call zaxpy(size_evc,-cmplx(beta,0.0d0,kind=dp),evc1_old(:,:,:,2),1,evc1_new(:,:,:,2),1) + call zaxpy(size_evc,-cmplx(beta,0.0d0,kind=dp),evc1_old(:,:,:,2),1,evc1_new(:,:,:,2),1) + if (lr_verbosity >10) then + write(stdout,'("evc1(1), final")') + do ibnd=1,nbnd + call check_vector_gamma(evc1_new(:,ibnd,1,1)) + enddo + write(stdout,'("evc1(2), final")') + do ibnd=1,nbnd + call check_vector_gamma(evc1_new(:,ibnd,1,2)) + enddo + endif ! ! Writing files for restart ! diff --git a/TDDFPT/src/lr_main.f90 b/TDDFPT/src/lr_main.f90 index b0bff9645..a427e82b3 100755 --- a/TDDFPT/src/lr_main.f90 +++ b/TDDFPT/src/lr_main.f90 @@ -34,14 +34,15 @@ PROGRAM lr_main USE control_ph, ONLY : nbnd_occ use wvfct, only : nbnd - + !Debugging + USE lr_variables, ONLY: check_all_bands_gamma, check_density_gamma,check_vector_gamma ! IMPLICIT NONE ! ! Local variables ! !INTEGER :: iter, ip, op - INTEGER :: ip,pol_index,ibnd_occ,ibnd_virt + INTEGER :: ip,pol_index,ibnd_occ,ibnd_virt,ibnd INTEGER :: iter_restart,iteration LOGICAL :: rflag CHARACTER (len=9) :: code = 'TDDFPT' @@ -89,8 +90,6 @@ PROGRAM lr_main endif !OBM_DEBUG ! - !if (restart) call test_restart() - ! ! !Initialisation of degauss/openshell related stuff ! @@ -138,7 +137,7 @@ PROGRAM lr_main CALL lr_read_d0psi() else CALL lr_solve_e() - endif + endif ! ! Set up initial stuff for derivatives ! @@ -200,10 +199,28 @@ PROGRAM lr_main ! write(stdout,'(/5x,"Starting Lanczos loop",1x,i8)') LR_polarization ! - END IF - ! - CALL sd0psi() + END IF + if (lr_verbosity >10) then + write(stdout,'("d0psi")') + do ibnd=1,nbnd + call check_vector_gamma(d0psi(:,ibnd,1,pol_index)) + enddo + endif + + ! + CALL sd0psi() !after this d0psi is Sd0psi ! + if (lr_verbosity >10) then + write(stdout,'("initial evc1")') + do ibnd=1,nbnd + call check_vector_gamma(evc1(:,ibnd,1,1)) + enddo + write(stdout,'("initial sd0psi")') + do ibnd=1,nbnd + call check_vector_gamma(d0psi(:,ibnd,1,pol_index)) + enddo + endif + ! lancz_loop1 : DO iteration = iter_restart, itermax ! diff --git a/TDDFPT/src/lr_ortho.f90 b/TDDFPT/src/lr_ortho.f90 index b78238d1f..c41415e46 100755 --- a/TDDFPT/src/lr_ortho.f90 +++ b/TDDFPT/src/lr_ortho.f90 @@ -46,7 +46,7 @@ COMPLEX(DP), INTENT(INOUT) :: dvpsi(npwx*npol,nbnd) COMPLEX(DP), INTENT(IN) :: sevc(npwx*npol,nbnd) ! work space allocated by ! the calling routine (was called dpsi) !real(kind=dp), intent(IN) :: lr_alpha_pv !This is calculated manually in tddfpt -logical, optional :: inverse !if .true. |dvspi> = |dvpsi> - |evq> instead of |dvspi> = |dvpsi> - |sevc> +logical, intent(in):: inverse !if .true. |dvspi> = |dvpsi> - |evq> instead of |dvspi> = |dvpsi> - |sevc> logical:: inverse_mode @@ -57,11 +57,11 @@ CALL start_clock ('lr_ortho') If (lr_verbosity > 5) WRITE(stdout,'("")') ! - if (.not. present(inverse)) then - inverse_mode=.false. - else + !if (.not. present(inverse)) then + ! inverse_mode=.false. + !else inverse_mode=inverse - endif + !endif if (gamma_only) then ! call lr_ortho_gamma() @@ -205,116 +205,73 @@ contains ! if (lgauss) then call errore ('lr_ortho', "degauss with gamma point algorithms",1) - ! - ! metallic case - ! - ps = (0.d0, 0.d0) - if (inverse_mode) then - CALL ZGEMM( 'C', 'N', nbnd, nbnd_occ (ikk), npw_k(ikk), (1.d0,0.d0), & - sevc, npwx, dvpsi, npwx, (0.d0,0.d0), ps, nbnd ) - else - CALL ZGEMM( 'C', 'N', nbnd, nbnd_occ (ikk), npw_k(ikk), (1.d0,0.d0), & - evq, npwx, dvpsi, npwx, (0.d0,0.d0), ps, nbnd ) - endif - ! - DO ibnd = 1, nbnd_occ (ikk) - wg1 = wgauss ((ef-et(ibnd,ikk)) / degauss, ngauss) - w0g = w0gauss((ef-et(ibnd,ikk)) / degauss, ngauss) / degauss - DO jbnd = 1, nbnd - wgp = wgauss ( (ef - et (jbnd, ikq) ) / degauss, ngauss) - deltae = et (jbnd, ikq) - et (ibnd, ikk) - theta = wgauss (deltae / degauss, 0) - wwg = wg1 * (1.d0 - theta) + wgp * theta - IF (jbnd <= nbnd_occ (ikq) ) THEN - IF (abs (deltae) > 1.0d-5) THEN - wwg = wwg + alpha_pv * theta * (wgp - wg1) / deltae - ELSE - ! - ! if the two energies are too close takes the limit - ! of the 0/0 ratio - ! - wwg = wwg - alpha_pv * theta * w0g - ENDIF - ENDIF - ! - ps(jbnd,ibnd) = wwg * ps(jbnd,ibnd) - ! - ENDDO - call DSCAL (2*npw_k(ikk), wg1, dvpsi(1,ibnd), 1) - END DO - nbnd_eff=nbnd ELSE ! ! insulators ! ps = ! in old version it was ps = - ps = 0.d0 - if (inverse_mode) then - CALL DGEMM( 'C', 'N', nbnd, nbnd ,2*npw_k(1), & - 2.d0, sevc, 2*npwx, dvpsi, 2*npwx, & - 0.d0, ps, nbnd ) - else - CALL DGEMM( 'C', 'N', nbnd, nbnd ,2*npw_k(1), & - 2.d0, evq, 2*npwx, dvpsi, 2*npwx, & - 0.d0, ps, nbnd ) - endif - nbnd_eff=nbnd - if (gstart == 2) then - if (inverse_mode) then - CALL DGER( nbnd, nbnd, -1.D0, sevc, 2*npwx, dvpsi, 2*npwx, ps, nbnd ) - else - CALL DGER( nbnd, nbnd, -1.D0, evq, 2*npwx, dvpsi, 2*npwx, ps, nbnd ) - endif - endif - END IF + ps = 0.d0 + if (inverse_mode) then + CALL DGEMM( 'C', 'N', nbnd, nbnd ,2*npw_k(1), & + 2.d0, sevc, 2*npwx, dvpsi, 2*npwx, & + 0.d0, ps, nbnd ) + !ps = 2* + else + CALL DGEMM( 'C', 'N', nbnd, nbnd ,2*npw_k(1), & + 2.d0, evq, 2*npwx, dvpsi, 2*npwx, & + 0.d0, ps, nbnd ) + !ps = 2* + endif + nbnd_eff=nbnd + if (gstart == 2) then + if (inverse_mode) then + CALL DGER( nbnd, nbnd, -1.D0, sevc, 2*npwx, dvpsi, 2*npwx, ps, nbnd ) + !PS = PS - sevc*dvpsi + else + CALL DGER( nbnd, nbnd, -1.D0, evq, 2*npwx, dvpsi, 2*npwx, ps, nbnd ) + !PS = PS - evc*dvpsi + endif + endif + END IF #ifdef __PARA - call mp_sum(ps(:,1:nbnd_eff),intra_pool_comm) + call mp_sum(ps(:,:),intra_pool_comm) #endif - ! in the original dpsi was used as a storage for sevc, since in - ! tddfpt we have it stored in memory as sevc0 this part is obsolote - !! - !! dpsi is used as work space to store S|evc> - !! - !IF (noncolin) THEN - ! IF (okvan) CALL calbec ( npw_k(ikk), vkb, evq, becp_nc, nbnd_eff ) - !ELSE - ! IF (okvan) CALL calbec ( npwq, vkb, evq, becp, nbnd_eff) - !ENDIF - !CALL s_psi (npwx, npwq, nbnd_eff, evq, dpsi) + ! in the original dpsi was used as a storage for sevc, since in + ! tddfpt we have it stored in memory as sevc0 this part is obsolote + !! + !! dpsi is used as work space to store S|evc> + !! + !IF (noncolin) THEN + ! IF (okvan) CALL calbec ( npw_k(ikk), vkb, evq, becp_nc, nbnd_eff ) + !ELSE + ! IF (okvan) CALL calbec ( npwq, vkb, evq, becp, nbnd_eff) + !ENDIF + !CALL s_psi (npwx, npwq, nbnd_eff, evq, dpsi) - ps_c = cmplx(ps, 0.d0, dp) - !! - !! |dvspi> = -(|dvpsi> - S|evq>) - !! - !OBM!!! changed to |dvspi> = |dvpsi> - |sevc> - if (lgauss) then - ! - ! metallic case - ! - if (inverse_mode) then - CALL ZGEMM( 'N', 'N', npw_k(1), nbnd_occ(1), nbnd, & - (1.d0,0.d0), evq, npwx, ps_c, nbnd, (-1.0d0,0.d0), & - dvpsi, npwx ) - else - CALL ZGEMM( 'N', 'N', npw_k(1), nbnd_occ(1), nbnd, & - (1.d0,0.d0), sevc, npwx, ps_c, nbnd, (-1.0d0,0.d0), & - dvpsi, npwx ) - endif - ELSE - ! - ! Insulators: note that nbnd_occ(ikk)=nbnd_occ(ikq) in an insulator - ! - if (inverse_mode) then - CALL ZGEMM( 'N', 'N', npw_k(1), nbnd, nbnd, & - (-1.d0,0.d0), evq, npwx, ps_c, nbnd, (1.0d0,0.d0), & - dvpsi, npwx ) - else - CALL ZGEMM( 'N', 'N', npw_k(1), nbnd, nbnd, & - (-1.d0,0.d0), sevc, npwx, ps_c, nbnd, (1.0d0,0.d0), & - dvpsi, npwx ) - endif - ENDIF + ps_c = cmplx(ps, 0.d0, dp) + !! + !! |dvspi> = -(|dvpsi> - S|evq>) + !! + !OBM!!! changed to |dvspi> = |dvpsi> - |sevc> + if (lgauss) then + !errore ? + ELSE + ! + ! Insulators: note that nbnd_occ(ikk)=nbnd_occ(ikq) in an insulator + ! + if (inverse_mode) then + CALL ZGEMM( 'N', 'N', npw_k(1), nbnd, nbnd, & + (-1.d0,0.d0), evq, npwx, ps_c, nbnd, (1.0d0,0.d0), & + dvpsi, npwx ) + !dvpsi=dvpsi-|evq> + else + CALL ZGEMM( 'N', 'N', npw_k(1), nbnd, nbnd, & + (-1.d0,0.d0), sevc, npwx, ps_c, nbnd, (1.0d0,0.d0), & + dvpsi, npwx ) + !dvpsi=dvpsi-|sevc> + endif + ENDIF DEALLOCATE(ps) DEALLOCATE(ps_c) end subroutine lr_ortho_gamma diff --git a/TDDFPT/src/lr_read_wf.f90 b/TDDFPT/src/lr_read_wf.f90 index 882049be5..57520b897 100755 --- a/TDDFPT/src/lr_read_wf.f90 +++ b/TDDFPT/src/lr_read_wf.f90 @@ -73,6 +73,8 @@ subroutine normal_read() ! !The usual way of reading wavefunctions ! +USE lr_variables, ONLY: check_all_bands_gamma, check_density_gamma,& + check_vector_gamma IMPLICIT NONE nwordwfc = 2 * nbnd * npwx @@ -287,6 +289,14 @@ subroutine normal_read() ! !print * , "evc0 ",evc0(1:3,1,1) ! + if (lr_verbosity >10) then + call check_all_bands_gamma(evc0(:,:,1),sevc0(:,:,1),nbnd,nbnd) + write(stdout,'("evc0")') + do ibnd=1,nbnd + call check_vector_gamma(evc0(:,ibnd,1)) + enddo + call check_density_gamma(revc0(:,:,1),nbnd) + endif ! !OBM!!! debug--- !CALL lr_normalise( evc0(:,:,1), obm_debug) @@ -306,6 +316,8 @@ subroutine virt_read() ! USE control_ph, ONLY : nbnd_occ use gvect, only : nrxx +USE lr_variables, ONLY: check_all_bands_gamma, check_density_gamma,& + check_vector_gamma IMPLICIT NONE complex(kind=dp), allocatable :: evc_all(:,:,:) complex(kind=dp), allocatable :: sevc_all(:,:,:) @@ -423,7 +435,7 @@ use gvect, only : nrxx ! Inverse fourier transform of evc_all ! ! - revc0=(0.0d0,0.0d0) + revc_all=(0.0d0,0.0d0) ! if ( gamma_only ) then ! @@ -481,8 +493,11 @@ use gvect, only : nrxx ! nbnd=nbnd_occ(1) ! + evc0=(0.0d0,0.0d0) evc0(:,:,:)=evc_all(:,1:nbnd,:) + sevc0=(0.0d0,0.0d0) sevc0(:,:,:)=sevc_all(:,1:nbnd,:) + revc0=(0.0d0,0.0d0) revc0(:,:,:)=revc_all(:,1:nbnd,:) if (nkb>0) then if (gamma_only) then @@ -506,6 +521,14 @@ use gvect, only : nrxx endif endif endif + if (lr_verbosity >10) then + call check_all_bands_gamma(evc_all(:,:,1),sevc_all(:,:,1),nbnd_total,nbnd) + call check_density_gamma(revc_all(:,:,1),nbnd) + write(stdout,'("evc0")') + do ibnd=1,nbnd + call check_vector_gamma(evc0(:,ibnd,1)) + enddo + endif if (nkb>0) then if (gamma_only) then deallocate(becp1_all) diff --git a/TDDFPT/src/lr_variables.f90 b/TDDFPT/src/lr_variables.f90 index a03605fb8..d70c3f0a9 100755 --- a/TDDFPT/src/lr_variables.f90 +++ b/TDDFPT/src/lr_variables.f90 @@ -170,5 +170,135 @@ module lr_variables ! !real(kind=dp) :: broadening !Broadening integer :: plot_type ! Dumps rho as: 1=xyzd 2=xsf 3=cube + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !! Debugging subroutines + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + contains +!---------------------------------------------------------------------------- +SUBROUTINE check_vector_gamma (x) +! Checks the inner product for a given vector, and its imaginary and real +! component +! input, evc +! output : screen output + USE mp_global, ONLY : inter_pool_comm, intra_pool_comm + use mp, only : mp_sum + use realus, only : npw_k + use gvect, only : gstart + use io_global, only : stdout + +IMPLICIT NONE + !input/output + complex(kind=dp),intent(in) :: x(:) + !local + real(kind=dp) :: temp_gamma + real(kind=dp), external :: DDOT + + temp_gamma = 2.D0*DDOT(2*npw_k(1),x(:),1,x(:),1) + if (gstart==2) temp_gamma = temp_gamma - dble(x(1))*dble(x(1)) +#ifdef __PARA + call mp_sum(temp_gamma, intra_pool_comm) +#endif + write(stdout,'(" = ",E15.8)') temp_gamma + END SUBROUTINE check_vector_gamma + +!---------------------------------------------------------------------------- +SUBROUTINE check_vector_f (x) +! Checks the inner product for a given vector, and its imaginary and real +! component +! input, evc +! output : screen output + USE mp_global, ONLY : inter_pool_comm, intra_pool_comm + use mp, only : mp_sum + use realus, only : npw_k + use gvect, only : gstart + use io_global, only : stdout + +IMPLICIT NONE + !input/output + complex(kind=dp),intent(in) :: x(:) + !local + complex(kind=dp) :: temp_f + complex(kind=dp), external :: ZDOTC + + temp_f = ZDOTC(npw_k(1),x(:),1,x(:),1) +#ifdef __PARA + call mp_sum(temp_f, intra_pool_comm) +#endif + write(stdout,'(" = ",2E15.8,1X)') temp_f + END SUBROUTINE check_vector_f +!---------------------------------------------------------------------------- +SUBROUTINE check_all_bands_gamma (x,sx,nbnd1,nbnd2) +! Checks all bands of given KS states for orthoganilty +! input, evc and sevc +! output : screen output + USE mp_global, ONLY : inter_pool_comm, intra_pool_comm + use mp, only : mp_sum + use realus, only : npw_k + use io_global, only : stdout + use gvect, only : gstart +IMPLICIT NONE + !input/output + integer,intent(in) :: nbnd1,nbnd2 !Total number of bands for x and sx + complex(kind=dp),intent(in) :: x(:,:), sx(:,:) + !local + integer :: ibnd, jbnd + real(kind=dp) :: temp_gamma + real(kind=dp), external :: DDOT + + do ibnd=1,nbnd1 + do jbnd=ibnd,nbnd2 + ! + temp_gamma = 2.D0*DDOT(2*npw_k(1),x(:,ibnd),1,sx(:,jbnd),1) + if (gstart==2) temp_gamma = temp_gamma - dble(x(1,ibnd))*dble(sx(1,jbnd)) +#ifdef __PARA + call mp_sum(temp_gamma, intra_pool_comm) +#endif + write(stdout,'(" =",E15.8)') ibnd,jbnd,temp_gamma + enddo + enddo + END SUBROUTINE check_all_bands_gamma +!---------------------------------------------------------------------------- +SUBROUTINE check_density_gamma (rx,nbnd) +! Checks the contirbution of a given function transformed into real space +! input, revc +! output : screen output + USE mp_global, ONLY : inter_pool_comm, intra_pool_comm + use mp, only : mp_sum + use realus, only : npw_k + use wvfct, only : wg + use gvect, only : nrxx + use io_global, only : stdout + use cell_base, only : omega + +IMPLICIT NONE + !input/output + integer,intent(in) :: nbnd !Total number of bands for x and sx + complex(kind=dp),intent(in) :: rx(:,:) + !local + integer :: ibnd + real(kind=dp) :: temp_gamma,w1,w2 + + do ibnd=1,nbnd,2 + w1=wg(ibnd,1)/omega + ! + if(ibnd