From 79f0fd08c04ed4a28bbd391f531f3d7f1ac64e79 Mon Sep 17 00:00:00 2001 From: degironc Date: Tue, 6 Nov 2007 10:26:07 +0000 Subject: [PATCH] more scf_mod changes git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@4398 c92efa57-630b-4861-b058-cf58834340f0 --- Modules/make.depend | 2 + PH/make.depend | 29 +++++++++ PP/make.depend | 1 - PP/pw2casino.f90 | 4 +- PW/electrons.f90 | 142 +++++++++++------------------------------ PW/make.depend | 6 ++ PW/potinit.f90 | 19 +++++- PW/read_file.f90 | 4 +- PW/restart_in_ions.f90 | 4 +- PW/sum_band.f90 | 58 ++++++++++++++++- PW/update_pot.f90 | 4 +- PW/v_of_rho.f90 | 18 +++--- 12 files changed, 164 insertions(+), 127 deletions(-) diff --git a/Modules/make.depend b/Modules/make.depend index cf5f32e82..f1426a663 100644 --- a/Modules/make.depend +++ b/Modules/make.depend @@ -89,6 +89,7 @@ ions_nose.o : ions_base.o ions_nose.o : kind.o ions_nose.o : parameters.o ions_nose.o : timestep.o +kplusg.o : kind.o metadyn_base.o : basic_algebra_routines.o metadyn_base.o : cell_base.o metadyn_base.o : constants.o @@ -255,6 +256,7 @@ uspp.o : kind.o uspp.o : parameters.o uspp.o : pseudo_types.o uspp.o : random_numbers.o +vms.o : kind.o vxc_t.o : functionals.o vxc_t.o : io_global.o vxc_t.o : kind.o diff --git a/PH/make.depend b/PH/make.depend index 79dcc9c19..25e659ed3 100644 --- a/PH/make.depend +++ b/PH/make.depend @@ -109,6 +109,7 @@ bcast_ph_input.o : ramanm.o bcast_ph_input1.o : ../Modules/mp.o bcast_ph_input1.o : ../PW/pwcom.o bcast_ph_input1.o : phcom.o +bicgstab.o : ../Modules/kind.o ccg_psi.o : ../Modules/kind.o cch_psi_all.o : ../Modules/kind.o cch_psi_all.o : ../PW/becmod.o @@ -350,6 +351,12 @@ dynmatrix.o : ../PW/noncol.o dynmatrix.o : ../PW/pwcom.o dynmatrix.o : phcom.o dynmatrix.o : ramanm.o +ec_rpa.o : ../Modules/io_files.o +ec_rpa.o : ../Modules/io_global.o +ec_rpa.o : ../Modules/kind.o +ec_rpa.o : ../Modules/random_numbers.o +ec_rpa.o : ../PW/pwcom.o +ec_rpa.o : phcom.o ef_shift.o : ../Modules/cell_base.o ef_shift.o : ../Modules/io_global.o ef_shift.o : ../Modules/kind.o @@ -391,6 +398,7 @@ find_mode_sym.o : ../PW/noncol.o find_mode_sym.o : ../PW/pwcom.o find_mode_sym.o : phcom.o generate_dynamical_matrix_c.o : ../Modules/kind.o +gl_weight.o : phcom.o gmressolve_all.o : ../Modules/kind.o h_psiq.o : ../Modules/kind.o h_psiq.o : ../Modules/wavefunctions.o @@ -480,6 +488,7 @@ phonon.o : phcom.o phonon.o : ramanm.o phq_init.o : ../Modules/atom.o phq_init.o : ../Modules/constants.o +phq_init.o : ../Modules/grid_paw_variables.o phq_init.o : ../Modules/io_files.o phq_init.o : ../Modules/io_global.o phq_init.o : ../Modules/ions_base.o @@ -488,6 +497,7 @@ phq_init.o : ../Modules/uspp.o phq_init.o : ../Modules/wavefunctions.o phq_init.o : ../PW/noncol.o phq_init.o : ../PW/pwcom.o +phq_init.o : ../PW/rad_paw_routines.o phq_init.o : phcom.o phq_readin.o : ../Modules/constants.o phq_readin.o : ../Modules/control_flags.o @@ -689,6 +699,17 @@ solve_e_nscf.o : ../Modules/wavefunctions.o solve_e_nscf.o : ../PW/becmod.o solve_e_nscf.o : ../PW/pwcom.o solve_e_nscf.o : phcom.o +solve_edv.o : ../Modules/check_stop.o +solve_edv.o : ../Modules/control_flags.o +solve_edv.o : ../Modules/io_files.o +solve_edv.o : ../Modules/io_global.o +solve_edv.o : ../Modules/ions_base.o +solve_edv.o : ../Modules/kind.o +solve_edv.o : ../Modules/uspp.o +solve_edv.o : ../Modules/wavefunctions.o +solve_edv.o : ../PW/becmod.o +solve_edv.o : ../PW/pwcom.o +solve_edv.o : phcom.o solve_linter.o : ../Modules/check_stop.o solve_linter.o : ../Modules/constants.o solve_linter.o : ../Modules/io_files.o @@ -782,6 +803,10 @@ transform_int_so.o : ../PW/pwcom.o transform_int_so.o : phcom.o trntnsc.o : ../Modules/kind.o trntnsr_3.o : ../Modules/kind.o +vc_dv.o : ../Modules/functionals.o +vc_dv.o : ../Modules/kind.o +vc_dv.o : ../PW/pwcom.o +vc_dv.o : phcom.o write_dyn_on_file.o : ../Modules/kind.o write_epsilon_and_zeu.o : ../Modules/io_global.o write_epsilon_and_zeu.o : ../Modules/kind.o @@ -829,6 +854,7 @@ allocate_part.o : ../include/f_defs.h allocate_phq.o : ../include/f_defs.h bcast_ph_input.o : ../include/f_defs.h bcast_ph_input1.o : ../include/f_defs.h +bicgstab.o : ../include/f_defs.h cch_psi_all.o : ../include/f_defs.h cgsolve_all.o : ../include/f_defs.h ch_psi_all.o : ../include/f_defs.h @@ -865,6 +891,7 @@ dynmat0.o : ../include/f_defs.h dynmat_us.o : ../include/f_defs.h dynmatcc.o : ../include/f_defs.h dynmatrix.o : ../include/f_defs.h +ec_rpa.o : ../include/f_defs.h ef_shift.o : ../include/f_defs.h el_opt.o : ../include/f_defs.h elphon.o : ../include/f_defs.h @@ -915,6 +942,7 @@ solve_e.o : ../include/f_defs.h solve_e2.o : ../include/f_defs.h solve_e_fpol.o : ../include/f_defs.h solve_e_nscf.o : ../include/f_defs.h +solve_edv.o : ../include/f_defs.h solve_linter.o : ../include/f_defs.h star_q.o : ../include/f_defs.h sym_and_write_zue.o : ../include/f_defs.h @@ -937,6 +965,7 @@ transform_int_nc.o : ../include/f_defs.h transform_int_so.o : ../include/f_defs.h trntnsc.o : ../include/f_defs.h trntnsr_3.o : ../include/f_defs.h +vc_dv.o : ../include/f_defs.h xk_wk_collect.o : ../include/f_defs.h zstar_eu.o : ../include/f_defs.h zstar_eu_us.o : ../include/f_defs.h diff --git a/PP/make.depend b/PP/make.depend index 4e92291b8..2e3baa2dd 100644 --- a/PP/make.depend +++ b/PP/make.depend @@ -108,7 +108,6 @@ dipole.o : ../Modules/mp.o dipole.o : ../Modules/mp_global.o dipole.o : ../PW/para.o dipole.o : ../PW/pwcom.o -dipole.o : ../PW/scf_mod.o do_initial_state.o : ../Modules/cell_base.o do_initial_state.o : ../Modules/io_global.o do_initial_state.o : ../Modules/ions_base.o diff --git a/PP/pw2casino.f90 b/PP/pw2casino.f90 index 7829c6d8d..9d3aad435 100644 --- a/PP/pw2casino.f90 +++ b/PP/pw2casino.f90 @@ -69,7 +69,7 @@ SUBROUTINE compute_casino nrxx, g, gg, ecutwfc, gcutm, nl, igtongl USE klist , ONLY: nks, nelec, xk USE lsda_mod, ONLY: lsda, nspin - USE scf, ONLY: rho, rho_core, rhog_core + USE scf, ONLY: rho, rho_core, rhog_core, tauk USE vlocal, ONLY: vloc, vnew, strf USE wvfct, ONLY: npw, npwx, nbnd, gamma_only, igk, g2kin, wg, et USE uspp, ONLY: nkb, vkb, dvan @@ -219,7 +219,7 @@ SUBROUTINE compute_casino ! ! compute hartree and xc contribution ! - CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, & + CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, tauk, & ehart, etxc, vtxc, etotefield, charge, vnew ) ! etot=(ek + (etxc-etxcc)+ehart+eloc+enl+ewld) diff --git a/PW/electrons.f90 b/PW/electrons.f90 index 6e01291d1..0be1ac55c 100644 --- a/PW/electrons.f90 +++ b/PW/electrons.f90 @@ -109,8 +109,8 @@ SUBROUTINE electrons() ! ... the charge density and of the HXC-potential ! type (scf_type) :: rhoin ! used to store rho_in of current/next iteration - COMPLEX(DP), ALLOCATABLE :: taukgnew(:,:) - REAL(DP), ALLOCATABLE :: tauknew(:,:) + COMPLEX(DP), ALLOCATABLE :: taukgin(:,:) + REAL(DP), ALLOCATABLE :: taukin(:,:) ! ! ... external functions ! @@ -208,6 +208,10 @@ SUBROUTINE electrons() becstep(:,:,:) = becsum(:,:,:) END IF call create_scf_type ( rhoin ) + IF ( dft_is_meta() ) then + ALLOCATE( taukin( nrxx, nspin ) ) + ALLOCATE( taukgin( ngm, nspin ) ) + END IF ! DO idum = 1, niter ! @@ -239,9 +243,11 @@ SUBROUTINE electrons() deband_hwf = delta_e() ! ! save input current density in rhoin - rhoin%of_r(:,:) = rho%of_r(:,:) - rhoin%of_g(:,:) = rho%of_g(:,:) - !rhoin = rho + rhoin = rho + if ( dft_is_meta() ) then + taukin(:,:) = tauk(:,:) + taukgin(:,:) = taukg(:,:) + end if ! scf_step: DO ! @@ -278,74 +284,10 @@ SUBROUTINE electrons() ! PAW : sum_band DOES NOT compute new one-center charges CALL sum_band() ! - ! ... bring output charge density (now stored in R-space in rho%of_r) - ! ... also to G-space (rho%of_g) for mixing - ! - DO is = 1, nspin - ! - psic(:) = rho%of_r(:,is) - ! - CALL cft3( psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1 ) - ! - rho%of_g(:,is) = psic(nl(:)) - ! - IF ( okvan .AND. tqr ) THEN - ! - ! ... in case the augmentation charges are computed in real space - ! ... we apply an FFT filter to the density in real space to - ! ... remove features that are not compatible with the FFT grid. - ! - psic(:) = ( 0.D0, 0.D0 ) - ! - psic(nl(:)) = rho%of_g(:,is) - ! - IF ( gamma_only ) psic(nlm(:)) = CONJG( rho%of_g(:,is) ) - ! - CALL cft3( psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1 ) - ! - rho%of_r(:,is) = psic(:) - ! - END IF - ! - END DO - ! - ! ... the same for tauk -> rhognew - ! - IF ( dft_is_meta()) THEN - ALLOCATE( taukgnew( ngm, nspin ) ) - DO is = 1, nspin - ! - psic(:) = tauk(:,is) - ! - CALL cft3( psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1 ) - ! - taukgnew(:,is) = psic(nl(:)) - ! - IF ( okvan .AND. tqr ) THEN - ! - ! ... in case the augmentation terms are computed in real space - ! ... we apply an FFT filter to the density in real space to - ! ... remove features that are not compatible with the FFT grid - ! - psic(:) = ( 0.D0, 0.D0 ) - ! - psic(nl(:)) = taukgnew(:,is) - ! - IF ( gamma_only ) psic(nlm(:)) = CONJG( taukgnew(:,is) ) - ! - CALL cft3( psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1 ) - ! - tauk(:,is) = psic(:) - ! - END IF - ! - END DO - END IF - ! ! ... the Harris-Weinert-Foulkes energy is computed here using only ! ... quantities obtained from the input density ! - hwf_energy = eband + deband_hwf + ( etxc - etxcc ) + ewld + ehart + demet + hwf_energy = eband + deband_hwf + (etxc - etxcc) + ewld + ehart + demet ! IF ( lda_plus_u ) THEN ! @@ -390,7 +332,7 @@ SUBROUTINE electrons() becnew(:,:,:) = becsum(:,:,:) END IF ! - CALL mix_rho( rho%of_g, rhoin%of_g, taukgnew, taukg, becnew, becstep, & + CALL mix_rho( rho%of_g, rhoin%of_g, taukg, taukgin, becnew, becstep, & nsnew, ns, mixing_beta, dr2, tr2_min, iter, nmix, conv_elec ) ! ! ... if convergence is achieved or if the self-consistency error @@ -398,19 +340,11 @@ SUBROUTINE electrons() ! ... (tr2_min), rhoin%of_g and rho%of_g are unchanged: rhoin%of_g ! ... contains the input density and rho%of_g contains the output ! ... density, both in G-space. - ! ... in the other cases rhoin%of_g now contains mixed charge density - ! ... (the new input density) in G-space. + ! ... In the other cases rhoin%of_g now contains mixed charge density + ! ... (the new input density) in G-space. In this case + ! NB: In this later case mix_rho leaves out-of-sync rho(in)%of_r and + ! rho(in)%of_g. This should be fixed as it is error-prone. ! - IF ( conv_elec ) THEN - ! - ! ... if convergence is achieved, rho%of_g and rho%of_r contain - ! ... the same charge density, one in G-space, the other in R-space - ! - IF ( dft_is_meta() ) taukg(:,:) = taukgnew(:,:) - ! - END IF - ! - IF ( dft_is_meta() ) DEALLOCATE( taukgnew ) IF ( okpaw ) DEALLOCATE (becnew) ! IF ( first .and. nat > 0) THEN @@ -438,11 +372,11 @@ SUBROUTINE electrons() not_converged_electrons : & IF ( .NOT. conv_elec ) THEN ! - ! ... bring the mixed charge density (rho%of_g) from G- to - ! ... R-space (rhoin%of_r) + ! ... synchronize R- and G- components of the mixed charge density ! DO is = 1, nspin ! + ! use psic as working array psic(:) = ( 0.D0, 0.D0 ) ! psic(nl(:)) = rhoin%of_g(:,is) @@ -455,30 +389,28 @@ SUBROUTINE electrons() ! END DO ! - ! the same for the kinetic energy density (tauknew) + ! the same for the kinetic energy density (taukin) ! IF ( dft_is_meta() ) THEN - ALLOCATE( tauknew( nrxx, nspin ) ) DO is = 1, nspin ! psic(:) = ( 0.D0, 0.D0 ) ! - psic(nl(:)) = taukg(:,is) + psic(nl(:)) = taukgin(:,is) ! - IF ( gamma_only ) psic(nlm(:)) = CONJG( taukg(:,is) ) + IF ( gamma_only ) psic(nlm(:)) = CONJG( taukgin(:,is) ) ! CALL cft3( psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1 ) ! - tauknew(:,is) = psic(:) + taukin(:,is) = psic(:) ! END DO - ! END IF ! ! ... no convergence yet: calculate new potential from mixed ! ... charge density (i.e. the new estimate) ! - CALL v_of_rho( rhoin%of_r, rhoin%of_g, rho_core, rhog_core, & + CALL v_of_rho( rhoin%of_r, rhoin%of_g, rho_core, rhog_core, taukin, & ehart, etxc, vtxc, etotefield, charge, vr ) ! ! ... estimate correction needed to have variational energy: @@ -490,16 +422,13 @@ SUBROUTINE electrons() ! descf = delta_escf() ! - ! ... now copy the mixed charge density in R-space and G-space in rho + ! ... now copy the mixed charge density in R- and G-space in rho ! - !rho%of_r(:,:) = rhoin%of_r(:,:) - !rho%of_g(:,:) = rhoin%of_g(:,:) rho = rhoin - ! - IF ( dft_is_meta() ) THEN - tauk(:,:) = tauknew(:,:) - DEALLOCATE( tauknew ) - END IF + if ( dft_is_meta() ) then + tauk(:,:) = taukin(:,:) + taukg(:,:) = taukgin(:,:) + end if ! ! ... compute PAW corrections to descf IF (okpaw) THEN @@ -537,7 +466,7 @@ SUBROUTINE electrons() ! vnew(:,:) = vr(:,:) ! - CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, & + CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, tauk, & ehart, etxc, vtxc, etotefield, charge, vr ) ! vnew(:,:) = vr(:,:) - vnew(:,:) @@ -673,7 +602,7 @@ SUBROUTINE electrons() #endif correction1c = (deband_PAW + descf_PAW + e_PAW) PRINT '(5x,A,f10.6,A)', 'PAW correction: ',correction1c, ', composed of: ' - PRINT '(5x,A,f10.6,A,f10.6,A,f12.6)',& + PRINT '(5x,A,f10.6,A,f10.6,A,f10.6)',& 'de_band: ',deband_PAW,', de_scf: ',descf_PAW,', 1-center E: ', e_PAW IF (really_do_paw) THEN etot = etot + correction1c @@ -694,7 +623,7 @@ SUBROUTINE electrons() IF ( first ) THEN ! fock0 = exxenergy2() - CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, & + CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, tauk, & ehart, etxc, vtxc, etotefield, charge, vr ) ! CALL set_vrs( vrs, vltot, vr, nrxx, nspin, doublegrid ) @@ -824,6 +753,11 @@ SUBROUTINE electrons() DEALLOCATE(deband_1ae_na, deband_1ps_na, descf_1ae_na, descf_1ps_na) #endif END IF + ! + IF ( dft_is_meta() ) then + DEALLOCATE( taukgin ) + DEALLOCATE( taukin ) + ENDIF call destroy_scf_type ( rhoin ) ! RETURN @@ -1049,7 +983,7 @@ SUBROUTINE electrons() IF ( dft_is_meta() ) THEN DO ipol = 1, nspin delta_escf = delta_escf - & - SUM( (tauknew(:,ipol)-tauk(:,ipol) )*kedtaur(:,ipol)) + SUM( (taukin(:,ipol)-tauk(:,ipol) )*kedtaur(:,ipol)) END DO END IF ! diff --git a/PW/make.depend b/PW/make.depend index c2ae8bb73..0b1f1802a 100644 --- a/PW/make.depend +++ b/PW/make.depend @@ -814,7 +814,9 @@ paw_xc.o : exx.o paw_xc.o : noncol.o paw_xc.o : pwcom.o potinit.o : ../Modules/cell_base.o +potinit.o : ../Modules/constants.o potinit.o : ../Modules/control_flags.o +potinit.o : ../Modules/functionals.o potinit.o : ../Modules/grid_paw_variables.o potinit.o : ../Modules/io_files.o potinit.o : ../Modules/io_global.o @@ -948,6 +950,7 @@ read_conf_from_file.o : pw_restart.o read_conf_from_file.o : pwcom.o read_file.o : ../Modules/cell_base.o read_file.o : ../Modules/constants.o +read_file.o : ../Modules/grid_paw_variables.o read_file.o : ../Modules/io_files.o read_file.o : ../Modules/ions_base.o read_file.o : ../Modules/kind.o @@ -956,9 +959,11 @@ read_file.o : ../Modules/uspp.o read_file.o : ../Modules/wavefunctions.o read_file.o : ../Modules/xml_io_base.o read_file.o : buffers.o +read_file.o : grid_paw_routines.o read_file.o : noncol.o read_file.o : pw_restart.o read_file.o : pwcom.o +read_file.o : rad_paw_routines.o read_file.o : scf_mod.o read_ncpp.o : ../Modules/functionals.o read_ncpp.o : ../Modules/kind.o @@ -1243,6 +1248,7 @@ sum_band.o : ../Modules/wavefunctions.o sum_band.o : buffers.o sum_band.o : noncol.o sum_band.o : pwcom.o +sum_band.o : realus.o sum_band.o : scf_mod.o sumkg.o : ../Modules/kind.o sumkt.o : ../Modules/kind.o diff --git a/PW/potinit.f90 b/PW/potinit.f90 index 09cd7583e..ea9adf070 100644 --- a/PW/potinit.f90 +++ b/PW/potinit.f90 @@ -24,6 +24,7 @@ SUBROUTINE potinit() ! ... In all cases the scf potential is recalculated and saved in vr ! USE kinds, ONLY : DP + USE constants, ONLY : pi USE io_global, ONLY : stdout USE cell_base, ONLY : alat, omega USE ions_base, ONLY : nat, ityp, ntyp => nsp @@ -34,8 +35,9 @@ SUBROUTINE potinit() ngm, gstart, nl, g, gg USE gsmooth, ONLY : doublegrid USE control_flags, ONLY : lscf - USE scf, ONLY : rho, rho_core, rhog_core, & + USE scf, ONLY : rho, rho_core, rhog_core, tauk, taukg, & vltot, vr, vrs + USE funct, ONLY : dft_is_meta USE wavefunctions_module, ONLY : psic USE ener, ONLY : ehart, etxc, vtxc USE ldaU, ONLY : niter_with_fixed_ns @@ -63,6 +65,7 @@ SUBROUTINE potinit() ! REAL(DP) :: charge ! the starting charge REAL(DP) :: etotefield ! + REAL(DP) :: fact INTEGER :: is, ios INTEGER :: ldim ! integer variable for I/O control LOGICAL :: exst @@ -217,9 +220,21 @@ SUBROUTINE potinit() ! END DO ! + if ( dft_is_meta()) then + ! ... define a starting (TF) guess for tauk and taukg + fact = (3.d0*pi*pi)**(2.0/3.0) + DO is = 1, nspin + tauk(:,is) = fact * abs(rho%of_r(:,is)*nspin)**(5.0/3.0)/nspin + psic(:) = tauk(:,is) + CALL cft3( psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1 ) + taukg(:,is) = psic(nl(:)) + END DO + ! + end if + ! ! ... compute the potential and store it in vr ! - CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, & + CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, tauk, & ehart, etxc, vtxc, etotefield, charge, vr ) ! ! ... define the total local potential (external+scf) diff --git a/PW/read_file.f90 b/PW/read_file.f90 index 44420685d..a032cb575 100644 --- a/PW/read_file.f90 +++ b/PW/read_file.f90 @@ -31,7 +31,7 @@ SUBROUTINE read_file() nl, gstart USE gsmooth, ONLY : ngms, nls, nrx1s, nr1s, nr2s, nr3s USE spin_orb, ONLY : so, lspinorb - USE scf, ONLY : rho, rho_core, rhog_core, vr + USE scf, ONLY : rho, rho_core, rhog_core, tauk, vr USE wavefunctions_module, ONLY : psic USE vlocal, ONLY : strf USE io_files, ONLY : tmp_dir, prefix, iunpun, nwordwfc, iunwfc @@ -212,7 +212,7 @@ SUBROUTINE read_file() ! ! ... recalculate the potential ! - CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, & + CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, tauk, & ehart, etxc, vtxc, etotefield, charge, vr ) ! ! ... reads the wavefunctions and writes them in 'distributed' form diff --git a/PW/restart_in_ions.f90 b/PW/restart_in_ions.f90 index a3a8581d3..ca8095574 100644 --- a/PW/restart_in_ions.f90 +++ b/PW/restart_in_ions.f90 @@ -18,7 +18,7 @@ subroutine restart_in_ions (iter, ik_, dr2) nrxx USE klist, ONLY: nks USE lsda_mod, ONLY: nspin - USE scf, ONLY : rho, rho_core, rhog_core + USE scf, ONLY : rho, rho_core, rhog_core, tauk USE control_flags, ONLY: restart, tr2, ethr USE vlocal, ONLY: vnew USE wvfct, ONLY: nbnd, et @@ -80,7 +80,7 @@ subroutine restart_in_ions (iter, ik_, dr2) ! recalculate etxc, vtxc, ehart, needed by stress calculation ! - CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, & + CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, tauk, & ehart, etxc, vtxc, etotefield, charge, psic ) ! ! restart procedure completed diff --git a/PW/sum_band.f90 b/PW/sum_band.f90 index a7803477b..7a6ef76b9 100644 --- a/PW/sum_band.f90 +++ b/PW/sum_band.f90 @@ -21,13 +21,15 @@ SUBROUTINE sum_band() USE control_flags, ONLY : diago_full_acc USE cell_base, ONLY : at, bg, omega, tpiba USE ions_base, ONLY : nat, ntyp => nsp, ityp - USE gvect, ONLY : nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, ngm, g + USE gvect, ONLY : nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, & + ngm, g, nl, nlm USE gsmooth, ONLY : nls, nlsm, nr1s, nr2s, nr3s, & nrx1s, nrx2s, nrx3s, nrxxs, doublegrid USE klist, ONLY : nks, nkstot, wk, xk, ngk USE ldaU, ONLY : lda_plus_U USE lsda_mod, ONLY : lsda, nspin, current_spin, isk - USE scf, ONLY : rho, tauk + USE realus, ONLY : tqr + USE scf, ONLY : rho, tauk, taukg USE symme, ONLY : nsym, s, ftau USE io_files, ONLY : iunwfc, nwordwfc, iunigk USE buffers, ONLY : get_buffer @@ -185,7 +187,57 @@ SUBROUTINE sum_band() ! #endif if (dft_is_meta() ) deallocate (kplusg) - ! + ! + ! ... synchronize rho%of_g to the calculated rho%of_r + ! + DO is = 1, nspin + ! + ! use psic as work array + psic(:) = rho%of_r(:,is) + CALL cft3( psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1 ) + rho%of_g(:,is) = psic(nl(:)) + ! + IF ( okvan .AND. tqr ) THEN + ! ... in case the augmentation charges are computed in real space + ! ... we apply an FFT filter to the density in real space to + ! ... remove features that are not compatible with the FFT grid. + ! + psic(:) = ( 0.D0, 0.D0 ) + psic(nl(:)) = rho%of_g(:,is) + IF ( gamma_only ) psic(nlm(:)) = CONJG( rho%of_g(:,is) ) + CALL cft3( psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1 ) + rho%of_r(:,is) = psic(:) + ! + END IF + ! + END DO + ! + ! ... the same for tauk and taukg + ! + IF ( dft_is_meta()) THEN + DO is = 1, nspin + ! + ! use psic as work array + psic(:) = tauk(:,is) + CALL cft3( psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1 ) + taukg(:,is) = psic(nl(:)) + ! + IF ( okvan .AND. tqr ) THEN + ! ... in case the augmentation charges are computed in real space + ! ... we apply an FFT filter to the density in real space to + ! ... remove features that are not compatible with the FFT grid. + ! + psic(:) = ( 0.D0, 0.D0 ) + psic(nl(:)) = taukg(:,is) + IF ( gamma_only ) psic(nlm(:)) = CONJG( taukg(:,is) ) + CALL cft3( psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1 ) + tauk(:,is) = psic(:) + ! + END IF + ! + END DO + END IF + ! CALL stop_clock( 'sum_band' ) ! RETURN diff --git a/PW/update_pot.f90 b/PW/update_pot.f90 index 5076722b5..a0e67fc06 100644 --- a/PW/update_pot.f90 +++ b/PW/update_pot.f90 @@ -192,7 +192,7 @@ SUBROUTINE extrapolate_charge( rho_extr ) USE gvect, ONLY : nrxx, ngm, g, gg, gstart, nr1, nr2, nr3, & nl, eigts1, eigts2, eigts3, nrx1, nrx2, nrx3 USE lsda_mod, ONLY : lsda, nspin - USE scf, ONLY : rho, rho_core, rhog_core, vr + USE scf, ONLY : rho, rho_core, rhog_core, tauk, vr USE wavefunctions_module, ONLY : psic USE control_flags, ONLY : alpha0, beta0 USE ener, ONLY : ehart, etxc, vtxc @@ -378,7 +378,7 @@ SUBROUTINE extrapolate_charge( rho_extr ) ! END DO ! - CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, & + CALL v_of_rho( rho%of_r, rho%of_g, rho_core, rhog_core, tauk, & ehart, etxc, vtxc, etotefield, charge, vr ) ! IF ( ABS( charge - nelec ) / charge > 1.D-7 ) THEN diff --git a/PW/v_of_rho.f90 b/PW/v_of_rho.f90 index adddd0d35..75e9bc829 100644 --- a/PW/v_of_rho.f90 +++ b/PW/v_of_rho.f90 @@ -6,7 +6,7 @@ ! or http://www.gnu.org/copyleft/gpl.txt . ! !---------------------------------------------------------------------------- -SUBROUTINE v_of_rho( rho, rhog, rho_core, rhog_core, & +SUBROUTINE v_of_rho( rho, rhog, rho_core, rhog_core, tauk, & ehart, etxc, vtxc, etotefield, charge, v ) !---------------------------------------------------------------------------- ! @@ -23,7 +23,7 @@ SUBROUTINE v_of_rho( rho, rhog, rho_core, rhog_core, & ! IMPLICIT NONE ! - REAL(DP), INTENT(IN) :: rho(nrxx,nspin), rho_core(nrxx) + REAL(DP), INTENT(IN) :: rho(nrxx,nspin), rho_core(nrxx), tauk(nrxx,nspin) ! input: the valence charge ! input: the core charge COMPLEX(DP), INTENT(IN) :: rhog(ngm,nspin), rhog_core(ngm) @@ -43,7 +43,7 @@ SUBROUTINE v_of_rho( rho, rhog, rho_core, rhog_core, & ! ... calculate exchange-correlation potential ! if (dft_is_meta()) then - call v_xc_meta( rho, rhog, rho_core, rhog_core, etxc, vtxc, v ) + call v_xc_meta( rho, rhog, rho_core, rhog_core, tauk, etxc, vtxc, v ) else CALL v_xc( rho, rhog, rho_core, rhog_core, etxc, vtxc, v ) endif @@ -70,7 +70,7 @@ SUBROUTINE v_of_rho( rho, rhog, rho_core, rhog_core, & ! END SUBROUTINE v_of_rho ! -SUBROUTINE v_xc_meta( rho, rhog, rho_core, rhog_core, etxc, vtxc, v ) +SUBROUTINE v_xc_meta( rho, rhog, rho_core, rhog_core, tauk, etxc, vtxc, v ) !---------------------------------------------------------------------------- ! ! ... Exchange-Correlation potential Vxc(r) from n(r) @@ -86,7 +86,7 @@ SUBROUTINE v_xc_meta( rho, rhog, rho_core, rhog_core, etxc, vtxc, v ) ! IMPLICIT NONE ! - REAL(DP), INTENT(IN) :: rho(nrxx,nspin), rho_core(nrxx) + REAL(DP), INTENT(IN) :: rho(nrxx,nspin), rho_core(nrxx), tauk(nrxx,nspin) ! the valence charge ! the core charge COMPLEX(DP), INTENT(IN) :: rhog(ngm,nspin), rhog_core(ngm) @@ -123,7 +123,7 @@ SUBROUTINE v_xc_meta( rho, rhog, rho_core, rhog_core, etxc, vtxc, v ) rhoneg = 0.D0 ! ! IF (get_igcx()==7.AND.get_igcc()==6) THEN - call v_xc_tpss( rho, rhog, rho_core, rhog_core, etxc, vtxc, v ) + call v_xc_tpss( rho, rhog, rho_core, rhog_core, tauk, etxc, vtxc, v ) ! ELSE ! CALL errore('v_xc_meta','wrong igcx and/or igcc',1) ! ENDIF @@ -131,7 +131,7 @@ SUBROUTINE v_xc_meta( rho, rhog, rho_core, rhog_core, etxc, vtxc, v ) RETURN END SUBROUTINE v_xc_meta ! -SUBROUTINE v_xc_tpss( rho, rhog, rho_core, rhog_core, etxc, vtxc, v ) +SUBROUTINE v_xc_tpss( rho, rhog, rho_core, rhog_core, tauk, etxc, vtxc, v ) ! =================== !-------------------------------------------------------------------- ! use gvecp, only: ng => ngm @@ -139,14 +139,14 @@ SUBROUTINE v_xc_tpss( rho, rhog, rho_core, rhog_core, etxc, vtxc, v ) USE gvect, ONLY : nrxx, nrx1,nrx2,nrx3,nr1,nr2,nr3, & g,nl,ngm USE gsmooth, ONLY : nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, nrxxs - USE scf, ONLY : kedtau,tauk, kedtaur + USE scf, ONLY : kedtau, kedtaur USE lsda_mod, ONLY : nspin USE cell_base, ONLY : omega, alat USE constants, ONLY : e2 IMPLICIT NONE ! ! input - REAL(DP),INTENT(IN) :: rho(nrxx,nspin), rho_core(nrxx) + REAL(DP),INTENT(IN) :: rho(nrxx,nspin), rho_core(nrxx), tauk(nrxx,nspin) COMPLEX(DP),INTENT(IN) :: rhog(ngm,nspin), rhog_core(ngm) REAL(DP),INTENT(OUT) :: etxc, vtxc, v(nrxx,nspin) ! integer nspin , nnr