scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Some more bug fixes to TDDFPT and projection algorithm. 

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@6490 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
obm 2010-03-15 10:44:08 +00:00
parent 676cd43906
commit 3ec28c1a27
6 changed files with 179 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
TDDFPT/src/lr_apply_liouvillian.f90
View File

@ -77,7 +77,7 @@ subroutine lr_apply_liouvillian( evc1, evc1_new, sevc1_new, interaction )
allocate( spsi1(npwx, nbnd) )
spsi1(:,:)=(0.0d0,0.0d0)
!
if( interaction ) then
if( interaction ) then !If true, the full L is calculated
!
call lr_calc_dens( evc1, .false. )
!
@ -86,14 +86,14 @@ subroutine lr_apply_liouvillian( evc1, evc1_new, sevc1_new, interaction )
dvrs(:,1)=0.0d0
call interpolate (dvrs(:,1),dvrss,-1)
else
dvrs(:,1)=rho_1(:,1)
!
!call lr_dv_of_drho(dvrs)
allocate( dvrs_temp(nrxx, nspin) )
dvrs_temp=CMPLX(dvrs,0.0d0) !OBM: This memory copy was hidden in lr_dv_of_drho, can it be avoided?
dvrs(:,1)=rho_1(:,1)
!
!call lr_dv_of_drho(dvrs)
allocate( dvrs_temp(nrxx, nspin) )
dvrs_temp=CMPLX(dvrs,0.0d0) !OBM: This memory copy was hidden in lr_dv_of_drho, can it be avoided?
call dv_of_drho(0,dvrs_temp,.false.)
dvrs=DBLE(dvrs_temp)
deallocate(dvrs_temp)
deallocate(dvrs_temp)
!
if ( okvan ) then
if ( tqr ) then
@ -104,7 +104,7 @@ subroutine lr_apply_liouvillian( evc1, evc1_new, sevc1_new, interaction )
endif
call add_paw_to_deeq(d_deeq)
!
call interpolate (dvrs(:,1),dvrss,-1)
call interpolate (dvrs(:,1),dvrss,-1)
endif
!
endif

29
TDDFPT/src/lr_charg_resp.f90
View File

@ -60,7 +60,7 @@ subroutine read_wT_beta_gamma_z()
inquire (file = filename, opened = exst)
if (.not.exst) call errore(' lr_main ','Lanczos coefficents can not be opened ',1)
!
WRITE(stdout,'("Reading Pre-calculated lanczos coefficents from ",A50)') filename
WRITE(stdout,'(/,/5x,"Reading Pre-calculated lanczos coefficents from ",A50)') filename
allocate(discard2(w_T_npol))
!
read(158,*,end=301,err=302) iter_restart
@ -73,7 +73,7 @@ subroutine read_wT_beta_gamma_z()
read(158,*,end=301,err=303) discard
!print *, discard
!
write(stdout,'("--------------Lanczos Matrix-------------------")')
!write(stdout,'("--------------Lanczos Matrix-------------------")')
do i=1,itermax
!
!print *, "Iter=",i
@ -85,12 +85,12 @@ subroutine read_wT_beta_gamma_z()
!print *, discard2(:)
!
enddo
print *, "closing file"
!print *, "closing file"
!
close(158)
!
deallocate(discard2)
print *, "starting broadcast"
!print *, "starting broadcast"
#ifdef __PARA
end if
call mp_barrier()
@ -98,6 +98,7 @@ subroutine read_wT_beta_gamma_z()
call mp_bcast (w_T_gamma_store(:), ionode_id)
#endif
!print *, "broadcast complete"
WRITE(stdout,'(5x,I8,1x,"steps succesfully read for polarization index",1x,I3)') itermax,LR_polarization
CALL stop_clock( 'post-processing' )
return
301 call errore ('read_beta_gamma_z', 'File is corrupted, no data', i )
@ -421,7 +422,7 @@ subroutine lr_dump_rho_tot_cube(rho,identifier)
ALLOCATE( rho_temp(dfftp%npp(1)+1) )
if (ionode) then
filename = trim(prefix) // "-" // identifier // "-pol" //trim(int_to_char(LR_polarization))// ".cube"
write(stdout,'(5X,"Writing Cube file for response charge density")')
write(stdout,'(/5X,"Writing Cube file for response charge density")')
!write(stdout, *) filename
!write(stdout,'(5X,"|rho|=",D15.8)') rho_sum
open (158, file = filename, form = 'formatted', status = 'replace', err=501)
@ -559,7 +560,7 @@ subroutine lr_dump_rho_tot_cube(rho,identifier)
!
filename = trim(prefix) // "-" // identifier // "-pol" //trim(int_to_char(LR_polarization))// ".cube"
write(stdout,'(5X,"Writing Cube file for response charge density")')
write(stdout,'(/5X,"Writing Cube file for response charge density")')
!write(stdout, *) filename
!write(stdout,'(5X,"|rho|=",D15.8)') rho_sum
open (158, file = filename, form = 'formatted', status = 'replace', err=501)
@ -580,7 +581,7 @@ subroutine lr_dump_rho_tot_cube(rho,identifier)
!C
!C ALL COORDINATES ARE GIVEN IN ATOMIC UNITS.
write(158,*) 'Cubfile created from TDDFPT calculation'
write(158,*) 'Cubefile created from TDDFPT calculation'
write(158,*) identifier
! origin is forced to (0.0,0.0,0.0)
write(158,'(I5,3F12.6)') nat, 0.0d0, 0.0d0, 0.0d0
@ -671,7 +672,7 @@ subroutine lr_dump_rho_tot_xyzd(rho,identifier)
ALLOCATE( kowner( nr3 ) )
if (ionode) then
filename = trim(prefix) // "-" // identifier // "-pol" //trim(int_to_char(LR_polarization))// ".xyzd"
write(stdout,'(5X,"Writing xyzd file for response charge density")')
write(stdout,'(/5X,"Writing xyzd file for response charge density")')
!write(stdout, *) filename
!write(stdout,'(5X,"|rho|=",D15.8)') rho_sum
open (158, file = filename, form = 'formatted', status = 'replace', err=501)
@ -767,7 +768,7 @@ subroutine lr_dump_rho_tot_xyzd(rho,identifier)
!
filename = trim(prefix) // "-" // identifier // "-pol" //trim(int_to_char(LR_polarization))// ".xyzd"
write(stdout,'(5X,"Writing xyzd file for response charge density")')
write(stdout,'(/5X,"Writing xyzd file for response charge density")')
!write(stdout, *) filename
!write(stdout,'(5X,"|rho|=",D15.8)') rho_sum
open (158, file = filename, form = 'formatted', status = 'replace', err=501)
@ -858,7 +859,7 @@ subroutine lr_dump_rho_tot_xcrys(rho, identifier)
!
!
filename = trim(prefix) // "-" // identifier // "-pol" //trim(int_to_char(LR_polarization))// ".xsf"
write(stdout,'(5X,"Writing xsf file for response charge density")')
write(stdout,'(/5X,"Writing xsf file for response charge density")')
!write(stdout, *) filename
open (158, file = filename, form = 'formatted', status = 'replace', err=501)
@ -983,7 +984,7 @@ subroutine lr_dump_rho_tot_xcrys(rho, identifier)
!
!
filename = trim(prefix) // "-" // identifier // "-pol" //trim(int_to_char(LR_polarization))// ".xsf"
write(stdout,'(5X,"Writing xsf file for response charge density")')
write(stdout,'(/5X,"Writing xsf file for response charge density")')
!write(stdout, *) filename
open (158, file = filename, form = 'formatted', status = 'replace', err=501)
@ -1110,7 +1111,7 @@ subroutine lr_dump_rho_tot_pxyd(rho,identifier)
!
filename = trim(prefix) // "-" // identifier // "-pol" //trim(int_to_char(LR_polarization))// ".pxyd"
write(stdout,'(5X,"Writing z plane averaged pxyd file for response charge density")')
write(stdout,'(/5X,"Writing z plane averaged pxyd file for response charge density")')
!write(stdout, *) filename
write(stdout,'(5X,"|rho|=",D15.8)') rho_sum
open (158, file = filename, form = 'formatted', status = 'replace', err=501)
@ -1186,7 +1187,7 @@ IMPLICIT none
ipol=1
endif
! The S term for ultrasoft calculation
IF ( okvan ) then
IF ( okvan) then
!
!
do ibnd_occ = 1, nbnd
@ -1272,7 +1273,7 @@ IMPLICIT none
!
!
!and finally
!and finally (note:parellization handled in dot product, each node has the copy of F)
!
F(ibnd_occ,ibnd_virt,ipol)=F(ibnd_occ,ibnd_virt,ipol)+2.0d0*SSUM*w_T(LR_iteration)
if (lr_verbosity>9) print *, "ibnd_occ=",ibnd_occ," ibnd_virt=",ibnd_virt," SSUM=",SSUM," w_T=",w_T(LR_iteration)," F=",F(ibnd_occ,ibnd_virt,ipol)

43
TDDFPT/src/lr_lanczos.f90
View File

@ -24,23 +24,20 @@ contains
use klist, only : nks,xk
use lr_variables, only : n_ipol, ltammd, itermax,&
evc1, evc1_new, sevc1_new, evc1_old, &
evc0, sevc0, d0psi, rho_1_tot,&
evc0, sevc0, d0psi, &
alpha_store, beta_store, gamma_store, zeta_store,&
restart_step, iunrestart, nwordrestart,&
charge_response, size_evc, LR_polarization, LR_iteration,&!,real_space
charge_response, size_evc, LR_polarization, LR_iteration,&!,real_space
test_case_no
use uspp, only : vkb, nkb, okvan
use wvfct, only : nbnd, npwx, npw
use control_flags, only : gamma_only,tqr
use becmod, only : bec_type, becp, calbec
use gvect, only : nrxx
!use real_beta, only : ccalbecr_gamma,s_psir,fft_orbital_gamma,bfft_orbital_gamma
USE realus, ONLY : real_space, fft_orbital_gamma, initialisation_level, &
bfft_orbital_gamma, calbec_rs_gamma, add_vuspsir_gamma, &
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
USE noncollin_module, ONLY : nspin_mag
!
implicit none
@ -60,7 +57,7 @@ contains
!
!integer :: n
!
logical :: exst
!
If (lr_verbosity > 5) THEN
WRITE(stdout,'("<lr_lanczos_one_step>")')
@ -309,25 +306,25 @@ contains
!
! Writing files for restart
!
if ( mod(LR_iteration,restart_step)==0 .or. LR_iteration==itermax ) then
!if ( mod(LR_iteration,restart_step)==0 .or. LR_iteration==itermax ) then
! !
! nwordrestart = 2 * nbnd * npwx * nks
! !
! call diropn ( iunrestart, 'restart_lanczos.'//trim(int_to_char(LR_polarization)), nwordrestart, exst)
!
nwordrestart = 2 * nbnd * npwx * nks
! call davcio(evc1(:,:,:,1),nwordrestart,iunrestart,1,1)
! call davcio(evc1(:,:,:,2),nwordrestart,iunrestart,2,1)
! call davcio(evc1_new(:,:,:,1),nwordrestart,iunrestart,3,1)
! call davcio(evc1_new(:,:,:,2),nwordrestart,iunrestart,4,1)
! !
! close( unit = iunrestart)
! if (charge_response == 2 ) then
! call diropn ( iunrestart, 'restart_lanczos-rho_tot.'//trim(int_to_char(LR_polarization)), 2*nrxx, exst)
! call davcio(rho_1_tot(:,:),2*nrxx*nspin_mag,iunrestart,1,1)
! close( unit = iunrestart)
! endif
!
call diropn ( iunrestart, 'restart_lanczos.'//trim(int_to_char(LR_polarization)), nwordrestart, exst)
!
call davcio(evc1(:,:,:,1),nwordrestart,iunrestart,1,1)
call davcio(evc1(:,:,:,2),nwordrestart,iunrestart,2,1)
call davcio(evc1_new(:,:,:,1),nwordrestart,iunrestart,3,1)
call davcio(evc1_new(:,:,:,2),nwordrestart,iunrestart,4,1)
!
close( unit = iunrestart)
if (charge_response == 2 ) then
call diropn ( iunrestart, 'restart_lanczos-rho_tot.'//trim(int_to_char(LR_polarization)), 2*nrxx, exst)
call davcio(rho_1_tot(:,:),2*nrxx*nspin_mag,iunrestart,1,1)
close( unit = iunrestart)
endif
!
end if
!end if
!
call stop_clock('one_step')
!

18
TDDFPT/src/lr_main.f90
View File

@ -227,9 +227,9 @@ PROGRAM lr_main
!call lr_dump_rho_tot_xcrys(rho_1_tot,"summed-rho")
endif
if (project) then
write(stdout,'(/,/5x,"Projection of virtual states")')
write(stdout,'(/,/5x,"Projection of virtual states for polarization direction",1x,i8)') LR_polarization
write(stdout,'(5x,"occ",1x,"con",8x,"Re(F)",14x,"Im(F)",12x,"|F|",5x"% presence")')
sum_F=SUM(abs(F(:,:,:)))
sum_F=SUM(abs(F(:,:,ip)))
sum_F=sum_F/100.0d0
do ibnd_occ=1,nbnd
do ibnd_virt=1,(nbnd_total-nbnd)
@ -252,6 +252,20 @@ PROGRAM lr_main
!If (lr_verbosity > 3) THEN
! CALL lr_calculate_spectrum()
!endif
if (project) then
write(stdout,'(/,/5x,"Summed projection of virtual states")')
write(stdout,'(5x,"occ",1x,"con",8x,"Re(F)",14x,"Im(F)",12x,"|F|",5x"% presence")')
sum_F=SUM(abs(F(:,:,1)))+SUM(abs(F(:,:,2)))+SUM(abs(F(:,:,3)))
sum_F=sum_F/100.0d0
do ibnd_occ=1,nbnd
do ibnd_virt=1,(nbnd_total-nbnd)
write(stdout,'(5x,i3,1x,i3,3x,E16.8,2X,E16.8,2X,E16.8,2X,"%",F5.2)') &
ibnd_occ,ibnd_virt,sum(F(ibnd_occ,ibnd_virt,:)),sum(abs(F(ibnd_occ,ibnd_virt,:))),&
sum(abs(F(ibnd_occ,ibnd_virt,:)))/sum_F
enddo
enddo
endif
!
! Deallocate pw variables
!

65
TDDFPT/src/lr_restart.f90
View File

@ -15,9 +15,8 @@ subroutine lr_restart(iter_restart,rflag)
use cell_base, only : tpiba2
use gvect, only : g
use io_files, only : tmp_dir, prefix
use lr_variables, only : itermax
use lr_variables, only : evc1, evc1_new, sevc1_new, rho_1_tot
use lr_variables, only : restart, nwordrestart, iunrestart
use lr_variables, only : itermax,evc1, evc1_new, sevc1_new, rho_1_tot ,&
restart, nwordrestart, iunrestart,project,nbnd_total,F
use wvfct, only : npw, igk, nbnd, g2kin, npwx
use lr_variables, only : beta_store, gamma_store, zeta_store, norm0!,real_space
use becmod, only : bec_type, becp, calbec
@ -45,7 +44,7 @@ subroutine lr_restart(iter_restart,rflag)
!
! local variables
!
integer :: i,ibnd
integer :: i,ibnd,ibnd_occ,ibnd_virt,temp
integer :: ik, ig, ip
logical :: exst
character(len=256) :: tempfile, filename
@ -95,35 +94,77 @@ subroutine lr_restart(iter_restart,rflag)
!
end if
!
!
!Ionode only reads
!
#ifdef __PARA
if (ionode) then
#endif
!
! Read and broadcast beta gamma zeta
!
open (158, file = tempfile, form = 'formatted', status = 'old')
!
read(158,*) iter_restart
read(158,*,end=301,err=303) iter_restart
!
if ( iter_restart .ge. itermax ) iter_restart = itermax
!
read(158,*) norm0(pol_index)
read(158,*,end=301,err=303) norm0(pol_index)
!
do i=1,iter_restart
!
read(158,*) beta_store(pol_index,i)
read(158,*) gamma_store(pol_index,i)
read(158,*) zeta_store (pol_index,:,i)
read(158,*,end=301,err=303) beta_store(pol_index,i)
read(158,*,end=301,err=303) gamma_store(pol_index,i)
read(158,*,end=301,err=303) zeta_store (pol_index,:,i)
!
end do
!
close(158)
!
#ifdef __PARA
end if
call mp_bcast (iter_restart, ionode_id)
call mp_bcast (norm0(pol_index), ionode_id)
call mp_bcast (beta_store(pol_index,:), ionode_id)
call mp_bcast (gamma_store(pol_index,:), ionode_id)
call mp_bcast (zeta_store(pol_index,:,:), ionode_id)
#endif
!
!
! Read projection
!
if (project) then
filename = trim(prefix) // ".projection." // trim(int_to_char(LR_polarization))
tempfile = trim(tmp_dir) // trim(filename)
!
!
open (158, file = tempfile, form = 'formatted', status = 'unknown')
!
read(158,*,end=301,err=303) temp
!
if (temp /= iter_restart) call errore ('lr_restart', 'Iteration mismatch reading projections', 1 )
!
read(158,*,end=301,err=303) temp !number of filled bands
!
if (temp /= nbnd) call errore ('lr_restart', 'NBND mismatch reading projections', 1 )
!
read(158,*,end=301,err=303) temp !total number of bands
!
if (temp /= nbnd_total) call errore ('lr_restart', 'Total number of bands mismatch reading projections', 1 )
!
do ibnd_occ=1,nbnd
do ibnd_virt=1,(nbnd_total-nbnd)
read(158,*,end=301,err=303) F(ibnd_occ,ibnd_virt,pol_index)
enddo
enddo
!
close(158)
endif
#ifdef __PARA
call mp_bcast (F, ionode_id)
#endif
#ifdef __PARA
end if
#endif
!
iter_restart = iter_restart + 1
@ -210,5 +251,7 @@ subroutine lr_restart(iter_restart,rflag)
!
!
return
301 call errore ('restart', 'A File is corrupted, file ended unexpectedly', 1 )
303 call errore ('restart', 'A File is corrupted, error in reading data', 1)
end subroutine lr_restart
!-----------------------------------------------------------------------

71
TDDFPT/src/lr_write_restart.f90
View File

@ -9,9 +9,14 @@ subroutine lr_write_restart()
!
use io_files, only : tmp_dir, prefix
use lr_variables, only : beta_store, gamma_store, zeta_store, norm0, &
LR_polarization, LR_iteration, n_ipol
LR_polarization, LR_iteration, n_ipol,F,project,&
evc1,evc1_new,iunrestart, nwordrestart, rho_1_tot, nbnd_total,&
charge_response,lr_verbosity
use wvfct, only : nbnd, npwx, npw
use gvect, only : nrxx
USE io_global, ONLY : ionode
USE lr_variables, ONLY : lr_verbosity
use klist, only : nks
USE noncollin_module, ONLY : nspin_mag
USE io_global, ONLY : stdout
!
implicit none
@ -22,18 +27,26 @@ subroutine lr_write_restart()
!
! local variables
!
integer :: i, pol_index
integer :: i, pol_index,ibnd_occ,ibnd_virt
character(len=256) :: tempfile, filename
logical :: exst
!
If (lr_verbosity > 5) THEN
WRITE(stdout,'("<lr_write_restart>")')
endif
!
!ionode only operations:
!
pol_index=1
if ( n_ipol /= 1 ) pol_index=LR_polarization
#ifdef __PARA
if (ionode) then
#endif
!
!Writing beta gamma and zeta
!
!
filename = trim(prefix) // ".beta_gamma_z." // trim(int_to_char(LR_polarization))
tempfile = trim(tmp_dir) // trim(filename)
@ -55,10 +68,58 @@ subroutine lr_write_restart()
!
close(158)
!
!Writing F
!
if (project) then
filename = trim(prefix) // ".projection." // trim(int_to_char(LR_polarization))
tempfile = trim(tmp_dir) // trim(filename)
!
!
open (158, file = tempfile, form = 'formatted', status = 'unknown')
!
write(158,*) LR_iteration
!
write(158,*) nbnd !number of filled bands
!
write(158,*) nbnd_total !total number of bands
!
do ibnd_occ=1,nbnd
do ibnd_virt=1,(nbnd_total-nbnd)
write(158,*) F(ibnd_occ,ibnd_virt,pol_index)
enddo
enddo
!
close(158)
endif
!
#ifdef __PARA
end if
#endif
!
!
! Writing wavefuncion files for restart
!
!
nwordrestart = 2 * nbnd * npwx * nks
!
call diropn ( iunrestart, 'restart_lanczos.'//trim(int_to_char(LR_polarization)), nwordrestart, exst)
!
call davcio(evc1(:,:,:,1),nwordrestart,iunrestart,1,1)
call davcio(evc1(:,:,:,2),nwordrestart,iunrestart,2,1)
call davcio(evc1_new(:,:,:,1),nwordrestart,iunrestart,3,1)
call davcio(evc1_new(:,:,:,2),nwordrestart,iunrestart,4,1)
!
close( unit = iunrestart)
!
! Writing charge response density for restart
!
if (charge_response == 2 ) then
call diropn ( iunrestart, 'restart_lanczos-rho_tot.'//trim(int_to_char(LR_polarization)), 2*nrxx*nspin_mag, exst)
call davcio(rho_1_tot(:,:),2*nrxx*nspin_mag,iunrestart,1,1)
close( unit = iunrestart)
endif
!
return
end subroutine lr_write_restart
!-----------------------------------------------------------------------