From 7b022ce276d678235c61247d9430b53d7fbcde34 Mon Sep 17 00:00:00 2001 From: sbraccia Date: Mon, 10 Jan 2005 06:56:14 +0000 Subject: [PATCH] Order in wavefunctions extrapolation made independent from the order in potential extrapolation. The input keyword "potential_extrapolation" has been substituted by two separate keywords "pot_extrapolation" and "wfc_extrapolation". Default values are still 'atomic' for the former and 'none' for the latter. Documentation updated. C.S. git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1551 c92efa57-630b-4861-b058-cf58834340f0 --- Doc/INPUT_FPMD | 6 +- Doc/INPUT_PW | 28 ++-- Gamma/phcg.f90 | 7 +- Modules/control_flags.f90 | 3 +- Modules/input_parameters.f90 | 9 +- Modules/read_namelists.f90 | 44 ++++--- PW/close_files.f90 | 1 - PW/compute_scf.f90 | 18 ++- PW/hinit1.f90 | 67 +++++----- PW/input.f90 | 39 ++++-- PW/move_ions.f90 | 6 +- PW/openfil.f90 | 1 - PW/update_pot.f90 | 234 ++++++++++++++++++--------------- examples/example03/run_example | 3 +- examples/example17/run_example | 41 +++--- 15 files changed, 288 insertions(+), 219 deletions(-) diff --git a/Doc/INPUT_FPMD b/Doc/INPUT_FPMD index cb28b9835..f5f28a36e 100644 --- a/Doc/INPUT_FPMD +++ b/Doc/INPUT_FPMD @@ -692,13 +692,9 @@ INTEGER :: upscale ! This variable is NOT used in FPMD - CHARACTER(LEN=80) :: potential_extrapolation - ! This variable is used only by PWSCF - ! NOT used in FPMD - NAMELIST / ions / ion_dynamics, ion_radius, ion_damping, ion_positions, & ion_velocities, ion_temperature, tempw, fnosep, tranp, amprp, greasp, tolp, & - ion_nstepe, ion_maxstep, upscale, potential_extrapolation + ion_nstepe, ion_maxstep, upscale, pot_extrapolation, wfc_extrapolation =----------------------------------------------------------------------------=! diff --git a/Doc/INPUT_PW b/Doc/INPUT_PW index 574501880..d9ddd6668 100644 --- a/Doc/INPUT_PW +++ b/Doc/INPUT_PW @@ -480,19 +480,23 @@ nraise INTEGER (default = 100 ) This keyword is NOT used in the case of variable cell calculations. -potential_extrapolation - CHARACTER - used to extrapolate the potential and the wave-functions - from preceding ionic step(s) - - 'none': no extrapolation - 'atomic': extrapolate the potential as if it was a sum of - atomic-like orbitals (default for calculation='relax') - 'wfc': extrapolate the potential as above - extrapolate wave-functions with first-order formula - (default for calcualtion='md' and calcualtion='neb') - 'wfc2': as above, with second order formula +pot_extrapolation + CHARACTER ( default = "atomic" ) + used to extrapolate the potential from preceding ionic steps + 'none' : no extrapolation + 'atomic' : extrapolate the potential as if it was a sum of + atomic-like orbitals + 'first_rder' : extrapolate the potential with first-order + formula + 'second_order': as above, with second order formula +wfc_extrapolation + CHARACTER ( default = "none" ) + used to extrapolate the wavefunctions from preceding ionic steps + 'none' : no extrapolation + 'first_rder' : extrapolate the wave-functions with first-order + formula + 'second_order': as above, with second order formula ! ! ... keywords used only in BFGS calculations diff --git a/Gamma/phcg.f90 b/Gamma/phcg.f90 index 1c703caa3..488c317af 100644 --- a/Gamma/phcg.f90 +++ b/Gamma/phcg.f90 @@ -388,8 +388,8 @@ SUBROUTINE newscf USE funct USE io_files, ONLY : iunigk, iunwfc, input_drho, output_drho USE control_flags, ONLY : restart, reduce_io, lscf, istep, iprint, & - order, david, max_cg_iter, isolve, tr2, & - ethr, mixing_beta, nmix, niter + pot_order, wfc_order, david, max_cg_iter, & + isolve, tr2, ethr, mixing_beta, nmix, niter ! IMPLICIT NONE INTEGER :: iter @@ -408,7 +408,8 @@ SUBROUTINE newscf qcutz=0.0 istep=1 iprint=10000 - order=1 + pot_order=1 + wfc_order=0 input_drho=' ' output_drho=' ' ! diff --git a/Modules/control_flags.f90 b/Modules/control_flags.f90 index b4d0d19e5..2630a53ec 100644 --- a/Modules/control_flags.f90 +++ b/Modules/control_flags.f90 @@ -216,7 +216,8 @@ diis_buff, &! dimension of the buffer in diis diis_ndim, &! dimension of reduced basis in DIIS history, &! number of old steps available for potential updating - order ! type of potential updating ( see update_pot ) + pot_order, &! type of potential updating ( see update_pot ) + wfc_order ! type of wavefunctions updating ( see update_pot ) ! LOGICAL, PUBLIC :: & lfixatom, &! if .TRUE. some atom is kept fixed diff --git a/Modules/input_parameters.f90 b/Modules/input_parameters.f90 index 12eb8b4ec..5b49afdd2 100644 --- a/Modules/input_parameters.f90 +++ b/Modules/input_parameters.f90 @@ -907,8 +907,9 @@ MODULE input_parameters REAL(dbl) :: upscale = 0.0d0 ! This variable is NOT used in FPMD - CHARACTER(LEN=80) :: potential_extrapolation = 'default' - ! This variable is used only by PWSCF + CHARACTER(LEN=80) :: pot_extrapolation = 'default', & + wfc_extrapolation = 'default' + ! These variables are used only by PWSCF ! NOT used in FPMD ! @@ -1047,8 +1048,8 @@ MODULE input_parameters NAMELIST / ions / ion_dynamics, ion_radius, ion_damping, ion_positions, & ion_velocities, ion_temperature, tempw, fnosep, tranp, amprp, greasp, & - tolp, ion_nstepe, ion_maxstep, upscale, potential_extrapolation, & - delta_t, nraise, & + tolp, ion_nstepe, ion_maxstep, upscale, pot_extrapolation, & + wfc_extrapolation, delta_t, nraise, & num_of_images, CI_scheme, opt_scheme, first_last_opt, use_multistep, & reset_vel, write_save, damp, temp_req, ds, k_max, k_min, path_thr, & init_num_of_images, & diff --git a/Modules/read_namelists.f90 b/Modules/read_namelists.f90 index be6333a8d..a637a09e3 100644 --- a/Modules/read_namelists.f90 +++ b/Modules/read_namelists.f90 @@ -333,21 +333,22 @@ MODULE read_namelists_module ! ! ... ( 'nose' | 'not_controlled' | 'rescaling' ) ! - ion_temperature = 'not_controlled' - tempw = 300.D0 - fnosep = 1.D0 - tranp = .FALSE. - amprp = 0.D0 - greasp = 1.D0 - tolp = 100.D0 - ion_nstepe = 1 - ion_maxstep = 100 - delta_t = 1.D0 - nraise = 100 - upscale = 10 - potential_extrapolation = 'default' + ion_temperature = 'not_controlled' + tempw = 300.D0 + fnosep = 1.D0 + tranp = .FALSE. + amprp = 0.D0 + greasp = 1.D0 + tolp = 100.D0 + ion_nstepe = 1 + ion_maxstep = 100 + delta_t = 1.D0 + nraise = 100 + upscale = 10 + pot_extrapolation = 'atomic' + wfc_extrapolation = 'none' ! - ! ... defaults for "path" optimizations + ! ... defaults for "path" optimisations ! ! ... general input variables ! @@ -746,7 +747,8 @@ MODULE read_namelists_module CALL mp_bcast( delta_t, ionode_id ) CALL mp_bcast( nraise, ionode_id ) CALL mp_bcast( upscale, ionode_id ) - CALL mp_bcast( potential_extrapolation, ionode_id ) + CALL mp_bcast( pot_extrapolation, ionode_id ) + CALL mp_bcast( wfc_extrapolation, ionode_id ) ! ! ... "path" variables broadcast ! @@ -1456,16 +1458,24 @@ MODULE read_namelists_module ! ! ... "path" optimizations ! - IF( prog == 'PW' ) startingpot = 'atomic' - IF( prog == 'PW' ) startingwfc = 'atomic' + IF ( prog == 'PW' ) THEN + ! + startingpot = 'atomic' + startingwfc = 'atomic' + ! + END IF IF( prog == 'FP' ) THEN + ! electron_dynamics = 'sd' ion_dynamics = 'none' cell_dynamics = 'none' + ! END IF IF( prog == 'CP' ) THEN + ! electron_dynamics = 'damp' ion_dynamics = 'damp' + ! END IF CASE DEFAULT CALL errore( sub_name,' calculation '// & diff --git a/PW/close_files.f90 b/PW/close_files.f90 index 82e958250..736294a73 100644 --- a/PW/close_files.f90 +++ b/PW/close_files.f90 @@ -11,7 +11,6 @@ SUBROUTINE close_files() ! ! ... Close all files and synchronize processes for a new scf calculation. ! - USE control_flags, ONLY : order USE ldaU, ONLY : lda_plus_u USE io_files, ONLY : prefix, iunwfc, iunigk, iunat USE mp_global, ONLY : intra_image_comm diff --git a/PW/compute_scf.f90 b/PW/compute_scf.f90 index 7ac3118a1..2db78b7bd 100644 --- a/PW/compute_scf.f90 +++ b/PW/compute_scf.f90 @@ -21,7 +21,7 @@ SUBROUTINE compute_scf( N_in, N_fin, stat ) diago_thr_init USE constants, ONLY : e2 USE control_flags, ONLY : lneb, lsmd, conv_elec, istep, & - history, alpha0, beta0, ethr, order + history, alpha0, beta0, ethr, pot_order USE check_stop, ONLY : check_stop_now USE vlocal, ONLY : strf USE cell_base, ONLY : bg, alat @@ -213,7 +213,7 @@ SUBROUTINE compute_scf( N_in, N_fin, stat ) CALL mp_bcast( history, ionode_id, intra_image_comm ) CALL mp_bcast( tauold, ionode_id, intra_image_comm ) ! - IF ( conv_elec .AND. order > 0 .AND. history > 0 ) THEN + IF ( conv_elec .AND. history > 0 ) THEN ! ! ... potential and wavefunctions are extrapolated only if ! ... we are starting a new self-consistency (scf on the @@ -231,11 +231,15 @@ SUBROUTINE compute_scf( N_in, N_fin, stat ) CALL mp_bcast( alpha0, ionode_id, intra_image_comm ) CALL mp_bcast( beta0, ionode_id, intra_image_comm ) ! - ! ... structure factors of the old positions are computed - ! ... (needed for the old atomic charge) - ! - CALL struc_fact( nat, tauold(:,:,1), nsp, ityp, ngm, g, bg, & - nr1, nr2, nr3, strf, eigts1, eigts2, eigts3 ) + IF ( pot_order > 0 ) THEN + ! + ! ... structure factors of the old positions are computed + ! ... (needed for the old atomic charge) + ! + CALL struc_fact( nat, tauold(:,:,1), nsp, ityp, ngm, g, bg, & + nr1, nr2, nr3, strf, eigts1, eigts2, eigts3 ) + ! + END IF ! CALL update_pot() ! diff --git a/PW/hinit1.f90 b/PW/hinit1.f90 index 0576d4fe4..9dcd3ff3c 100644 --- a/PW/hinit1.f90 +++ b/PW/hinit1.f90 @@ -1,15 +1,16 @@ ! -! Copyright (C) 2001 PWSCF group +! Copyright (C) 2001-2004 PWSCF group ! This file is distributed under the terms of the ! GNU General Public License. See the file `License' ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! ! -!----------------------------------------------------------------------- -subroutine hinit1 - !----------------------------------------------------------------------- - ! Atomic configuration dependent hamiltonian initialization +!---------------------------------------------------------------------------- +SUBROUTINE hinit1() + !---------------------------------------------------------------------------- + ! + ! ... Atomic configuration dependent hamiltonian initialization ! USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau USE cell_base, ONLY : at, bg, omega, tpiba2 @@ -21,44 +22,50 @@ subroutine hinit1 USE lsda_mod, ONLY : nspin USE scf, ONLY : vrs, vltot, vr USE vlocal, ONLY : strf - USE control_flags, ONLY : order + USE control_flags, ONLY : pot_order ! - implicit none + IMPLICIT NONE ! - ! update the potential ! - call update_pot + ! ... update the potential ! - ! initialize structure factor array if it has not already been calculat - ! update_pot ( this is done if order > 0 ) + CALL update_pot() ! - if (order.eq.0) then - if (lmovecell) call scale_h - call struc_fact (nat, tau, ntyp, ityp, ngm, g, bg, nr1, nr2, & - nr3, strf, eigts1, eigts2, eigts3) + ! ... initialize structure factor array if it has not already been + ! ... calculated update_pot ( this is done if order > 0 ) + ! + IF ( pot_order == 0 ) THEN ! - ! calculate the core charge (if any) for the nonlinear core correction + IF ( lmovecell ) CALL scale_h() ! - call set_rhoc - endif + CALL struc_fact( nat, tau, ntyp, ityp, ngm, g, bg, & + nr1, nr2, nr3, strf, eigts1, eigts2, eigts3 ) + ! + ! ... calculate the core charge (if any) for the nonlinear + ! ... core correction + ! + CALL set_rhoc() + ! + END IF ! - ! calculate the total local potential + ! ... calculate the total local potential ! - call setlocal + CALL setlocal() ! - ! define the total local potential (external+scf) + ! ... define the total local potential (external+scf) ! - call set_vrs (vrs, vltot, vr, nrxx, nspin, doublegrid) + CALL set_vrs( vrs, vltot, vr, nrxx, nspin, doublegrid ) ! - ! update the D matrix + ! ... update the D matrix ! - call newd + CALL newd() ! - ! and recalculate the products of the S with the atomic wfcs used in LDA+U - ! calculations + ! ... and recalculate the products of the S with the atomic wfcs used + ! ... in LDA+U calculations ! - if (lda_plus_u) call orthoatwfc - - return -end subroutine hinit1 + IF ( lda_plus_u ) CALL orthoatwfc() + ! + RETURN + ! +END SUBROUTINE hinit1 diff --git a/PW/input.f90 b/PW/input.f90 index 1fdda9b96..3edf6a369 100644 --- a/PW/input.f90 +++ b/PW/input.f90 @@ -112,7 +112,7 @@ SUBROUTINE iosys() ! USE control_flags, ONLY : diis_ndim, isolve, & max_cg_iter, diis_buff, david, imix, nmix, & - iverbosity, tr2, niter, order, & + iverbosity, tr2, niter, pot_order, wfc_order, & upscale_ => upscale, & mixing_beta_ => mixing_beta, & nstep_ => nstep, & @@ -203,7 +203,7 @@ SUBROUTINE iosys() ! USE input_parameters, ONLY : ion_dynamics, ion_positions, ion_temperature, & tolp, tempw, delta_t, nraise, upscale, & - potential_extrapolation, & + pot_extrapolation, wfc_extrapolation, & num_of_images, path_thr, CI_scheme, opt_scheme, & reset_vel, use_multistep, first_last_opt, damp, & init_num_of_images, temp_req, k_max, k_min, ds, & @@ -604,19 +604,38 @@ SUBROUTINE iosys() tr2 = conv_thr niter = electron_maxstep ! - SELECT CASE ( TRIM( potential_extrapolation ) ) + SELECT CASE ( TRIM( pot_extrapolation ) ) CASE ( 'none' ) - order = 0 + pot_order = 0 CASE ( 'atomic' ) - order = 1 - CASE ( 'wfc' ) - order = 2 - CASE ( 'wfc2' ) - order = 3 + pot_order = 1 + CASE ( 'first_order' ) + pot_order = 2 + CASE ( 'second_order' ) + pot_order = 3 CASE DEFAULT - order = 1 + pot_order = 1 END SELECT ! + SELECT CASE ( TRIM( wfc_extrapolation ) ) + CASE ( 'none' ) + wfc_order = 0 + CASE ( 'first_order' ) + wfc_order = 2 + CASE ( 'second_order' ) + wfc_order = 3 + CASE DEFAULT + wfc_order = 0 + END SELECT + ! + IF ( wfc_order > 0 .AND. noncolin ) THEN + ! + CALL errore( ' iosys ', & + & ' wfc extrapolation not implemented in the' // & + & ' noncollinear case', -1 ) + ! + END IF + ! IF ( occupations == 'fixed' .AND. nspin == 2 .AND. lscf ) THEN CALL errore( ' iosys ', & & ' fixed occupations and lsda not implemented ', 1 ) diff --git a/PW/move_ions.f90 b/PW/move_ions.f90 index a1abc1eb0..2d0241464 100644 --- a/PW/move_ions.f90 +++ b/PW/move_ions.f90 @@ -355,7 +355,7 @@ SUBROUTINE find_alpha_and_beta( nat, tau, tauold, alpha0, beta0 ) USE constants, ONLY : eps16 USE kinds, ONLY : DP USE io_global, ONLY : stdout - USE control_flags, ONLY : order, history + USE control_flags, ONLY : history ! IMPLICIT NONE ! @@ -364,11 +364,11 @@ SUBROUTINE find_alpha_and_beta( nat, tau, tauold, alpha0, beta0 ) REAL(KIND=DP) :: a11, a12, a21, a22, b1, b2, c, det ! ! - IF ( MIN( history, order ) < 2 ) THEN + IF ( history < 2 ) THEN ! RETURN ! - ELSE IF ( MIN( history, order ) == 2 ) THEN + ELSE IF ( history == 2 ) THEN ! alpha0 = 1.D0 beta0 = 0.D0 diff --git a/PW/openfil.f90 b/PW/openfil.f90 index 6e75a214e..1f1e1d361 100644 --- a/PW/openfil.f90 +++ b/PW/openfil.f90 @@ -18,7 +18,6 @@ SUBROUTINE openfil() USE io_global, ONLY : stdout USE basis, ONLY : natomwfc, startingwfc USE wvfct, ONLY : nbnd, npwx - USE control_flags, ONLY : order, lneb USE ldaU, ONLY : lda_plus_U USE io_files, ONLY : prefix, iunpun, iunat, iunwfc, iunigk, & nwordwfc, nwordatwfc diff --git a/PW/update_pot.f90 b/PW/update_pot.f90 index 5b081275d..1765068a5 100644 --- a/PW/update_pot.f90 +++ b/PW/update_pot.f90 @@ -14,37 +14,54 @@ SUBROUTINE update_pot() !---------------------------------------------------------------------------- ! - ! ... update potential, use the integer variable order to decide the way + ! ... update the potential extrapolating the charge density and extrapolates + ! ... the wave-functions ! - ! ... order = 0 copy the old potential (nothing is done) + ! ... charge density extrapolation : ! - ! ... order = 1 subtract old atomic charge density and sum the new + ! ... pot_order = 0 copy the old potential (nothing is done) + ! + ! ... pot_order = 1 subtract old atomic charge density and sum the new ! ... if dynamics is done the routine extrapolates also ! ... the difference between the the scf charge and the ! ... atomic one, ! - ! ... order = 2 extrapolate the wavefunctions: + ! ... pot_order = 2 first order extrapolation : + ! + ! ... rho(t+dt) = 2*rho(t) - rho(t-dt) + ! + ! ... pot_order = 3 second order extrapolation : + ! + ! ... rho(t+dt) = rho(t) + + ! ... + alpha0*( rho(t) - rho(t-dt) ) + ! ... + beta0* ( rho(t-dt) - rho(t-2*dt) ) + ! + ! + ! ... wave-functions extrapolation : + ! + ! ... wfc_order = 0 nothing is done + ! + ! ... wfc_order = 2 first order extrapolation : ! ! ... |psi(t+dt)> = 2*|psi(t)> - |psi(t-dt)> ! - ! ... order = 3 extrapolate the wavefunctions with the second-order - ! ... formula: + ! ... wfc_order = 3 second order extrapolation : ! - ! ... |psi(t+dt)> = |psi(t) + + ! ... |psi(t+dt)> = |psi(t)> + ! ... + alpha0*( |psi(t)> - |psi(t-dt)> ) ! ... + beta0* ( |psi(t-dt)> - |psi(t-2*dt)> ) ! - ! ... where alpha0 and beta0 are calculated in - ! ... "find_alpha_and_beta()" so that |tau'-tau(t+dt)| is - ! ... minimum; - ! ... tau' and tau(t+dt) are respectively the atomic positions - ! ... at time t+dt and the extrapolated one: ! - ! ... tau(t+dt) = tau(t) + alpha0*( tau(t) - tau(t-dt) ) - ! ... + beta0*( tau(t-dt) -tau(t-2*dt) ) + ! ... alpha0 and beta0 are calculated in "find_alpha_and_beta()" so that + ! ... |tau'-tau(t+dt)| is minimum; + ! ... tau' and tau(t+dt) are respectively the atomic positions at time + ! ... t+dt and the extrapolated one: + ! + ! ... tau(t+dt) = tau(t) + alpha0*( tau(t) - tau(t-dt) ) + ! ... + beta0*( tau(t-dt) -tau(t-2*dt) ) ! ! - USE control_flags, ONLY : order, history + USE control_flags, ONLY : pot_order, wfc_order, history USE io_files, ONLY : prefix, tmp_dir, nd_nmbr USE io_global, ONLY : ionode, ionode_id USE mp, ONLY : mp_bcast @@ -52,13 +69,13 @@ SUBROUTINE update_pot() ! IMPLICIT NONE ! - INTEGER :: rho_order, wfc_order + INTEGER :: rho_extr, wfc_extr LOGICAL :: exists ! ! CALL start_clock( 'update_pot' ) ! - IF ( order == 0 ) THEN + IF ( pot_order == 0 .AND. wfc_order == 0 ) THEN ! CALL stop_clock( 'update_pot' ) ! @@ -69,47 +86,47 @@ SUBROUTINE update_pot() ! ... determines the maximum effective order of the extrapolation on the ! ... basis of the files that are really available ! - rho_order = MIN( 1, history ) + rho_extr = MIN( 1, history, pot_order ) ! INQUIRE( FILE = TRIM( tmp_dir ) // & & TRIM( prefix ) // '.oldrho', EXIST = exists ) ! IF ( exists ) THEN ! - rho_order = MIN( 2, history ) + rho_extr = MIN( 2, history, pot_order ) ! INQUIRE( FILE = TRIM( tmp_dir ) // & & TRIM( prefix ) // '.old2rho', EXIST = exists ) ! - IF ( exists ) rho_order = MIN( 3, history ) + IF ( exists ) rho_extr = MIN( 3, history, pot_order ) ! END IF ! - CALL extrapolate_charge( rho_order ) + IF ( pot_order > 0 ) CALL extrapolate_charge( rho_extr ) ! IF ( ionode ) THEN ! - wfc_order = MIN( 1, history, order ) + wfc_extr = MIN( 1, history, wfc_order ) ! INQUIRE( FILE = TRIM( tmp_dir ) // & & TRIM( prefix ) // '.oldwfc' // nd_nmbr, EXIST = exists ) ! IF ( exists ) THEN ! - wfc_order = MIN( 2, history, order ) + wfc_extr = MIN( 2, history, wfc_order ) ! INQUIRE( FILE = TRIM( tmp_dir ) // & & TRIM( prefix ) // '.old2wfc' // nd_nmbr , EXIST = exists ) ! - IF ( exists ) wfc_order = MIN( 3, history, order ) + IF ( exists ) wfc_extr = MIN( 3, history, wfc_order ) ! END IF ! END IF ! - CALL mp_bcast( wfc_order, ionode_id, intra_image_comm ) + CALL mp_bcast( wfc_extr, ionode_id, intra_image_comm ) ! - IF ( order >= 2 ) CALL extrapolate_wfcs( wfc_order ) + IF ( wfc_order > 0 ) CALL extrapolate_wfcs( wfc_extr ) ! CALL stop_clock( 'update_pot' ) ! @@ -119,31 +136,32 @@ END SUBROUTINE update_pot ! ! !---------------------------------------------------------------------------- -SUBROUTINE extrapolate_charge( rho_order ) +SUBROUTINE extrapolate_charge( rho_extr ) !---------------------------------------------------------------------------- ! - USE io_global, ONLY : stdout - USE kinds, ONLY : DP - USE cell_base, ONLY : omega, bg, alat - USE ions_base, ONLY : nat, tau, nsp, ityp - 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, vr - USE control_flags, ONLY : alpha0, beta0, imix - USE ener, ONLY : ehart, etxc, vtxc - USE extfield, ONLY : etotefield - USE cellmd, ONLY : lmovecell, omega_old - USE vlocal, ONLY : strf + USE constants, ONLY : eps32 + USE io_global, ONLY : stdout + USE kinds, ONLY : DP + USE cell_base, ONLY : omega, bg, alat + USE ions_base, ONLY : nat, tau, nsp, ityp + 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, vr + USE control_flags, ONLY : alpha0, beta0, imix + USE ener, ONLY : ehart, etxc, vtxc + USE extfield, ONLY : etotefield + USE cellmd, ONLY : lmovecell, omega_old + USE vlocal, ONLY : strf USE noncollin_module, ONLY : noncolin USE noncollin_module, ONLY : factlist, pointlist, pointnum, mcons,& i_cons, lambda, vtcon, report - USE io_files, ONLY : prefix - USE klist, ONLY : nelec + USE io_files, ONLY : prefix + USE klist, ONLY : nelec ! IMPLICIT NONE ! - INTEGER, INTENT(IN) :: rho_order + INTEGER, INTENT(IN) :: rho_extr ! REAL(KIND=DP), ALLOCATABLE :: work(:), work1(:) ! work is the difference between charge density and atomic charge @@ -154,34 +172,32 @@ SUBROUTINE extrapolate_charge( rho_order ) INTEGER :: ir, is ! ! - IF ( rho_order == 0 ) RETURN + IF ( rho_extr == 0 ) RETURN ! ALLOCATE( work(nrxx) ) ! work(:) = 0.D0 ! - ! ... if rho_order = 1 update the potential subtracting to the charge density - ! ... the "old" atomic charge and summing the new one - ! - WRITE( stdout,'(/5X,"NEW-OLD atomic charge density approx. for the potential")' ) - ! ! ... in the lsda case the magnetization will follow rigidly the density ! ... keeping fixed the value of zeta = mag / rho_tot. - ! ... zeta is set here and put in rho(* ??? while rho(*,1) will contain the + ! ... zeta is set here and put in rho(:,2) while rho(:,1) will contain the ! ... total valence charge ! IF ( lsda ) CALL rho2zeta( rho, rho_core, nrxx, nspin, 1 ) ! - IF (noncolin) THEN + IF ( noncolin ) THEN + ! DO is = 2, nspin - DO ir = 1, nrxx - IF ( rho(ir,1) .GT. 1.D-30 ) THEN - rho(ir,is) = rho(ir,is) / rho(ir,1) - ELSE - rho(ir,is) = 0.D0 - END IF + ! + WHERE( rho(:,1) > eps32 ) ! - END DO + rho(:,is) = rho(:,is) / rho(:,1) + ! + ELSE WHERE + ! + rho(:,is) = 0.D0 + ! + END WHERE ! END DO ! @@ -198,11 +214,21 @@ SUBROUTINE extrapolate_charge( rho_order ) ! ... extrapolate the difference between the atomic charge and ! ... the self-consistent one ! - IF ( rho_order == 1 ) THEN + IF ( rho_extr == 1 ) THEN + ! + ! ... if rho_extr = 1 update the potential subtracting to the charge + ! ... density the "old" atomic charge and summing the + ! ... new one + ! + WRITE( stdout, & + '(/5X,"NEW-OLD atomic charge density approx. for the potential")' ) ! CALL io_pot( + 1, TRIM( prefix )//'.oldrho', rho, 1 ) ! - ELSE IF ( rho_order == 2 ) THEN + ELSE IF ( rho_extr == 2 ) THEN + ! + WRITE( UNIT = stdout, & + FMT = '(/5X,"first order charge density extrapolation")' ) ! ! ... oldrho -> work ! @@ -214,11 +240,14 @@ SUBROUTINE extrapolate_charge( rho_order ) CALL io_pot( + 1, TRIM( prefix )//'.oldrho', rho, 1 ) CALL io_pot( + 1, TRIM( prefix )//'.old2rho', work, 1 ) ! - ! ... alpha0 has been calculated in move_ions + ! ... extrapolation ! - rho(:,1) = rho(:,1) + alpha0 * ( rho(:,1) - work(:) ) + rho(:,1) = 2.D0 * rho(:,1) - work(:) ! - ELSE IF ( rho_order == 3 ) THEN + ELSE IF ( rho_extr == 3 ) THEN + ! + WRITE( UNIT = stdout, & + FMT = '(/5X,"second order charge density extrapolation")' ) ! ALLOCATE( work1(nrxx) ) ! @@ -266,25 +295,33 @@ SUBROUTINE extrapolate_charge( rho_order ) ! IF ( lsda ) CALL rho2zeta( rho, rho_core, nrxx, nspin, -1 ) ! - IF (noncolin) THEN + IF ( noncolin ) THEN + ! DO is = 2, nspin - DO ir = 1, nrxx - IF ( rho(ir,1) .GT. 1.D-30 ) THEN - rho(ir,is) = rho(ir,is) * rho(ir,1) - ELSE - rho(ir,is) = 0.D0 - END IF - END DO + ! + WHERE( rho(:,1) > eps32 ) + ! + rho(:,is) = rho(:,is) * rho(:,1) + ! + ELSE WHERE + ! + rho(:,is) = 0.D0 + ! + END WHERE + ! END DO - - CALL v_of_rho_nc (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, & - nrxx, nl, ngm, gstart, nspin, g, gg, alat, omega, & - ehart, etxc, vtxc, charge, vr, lambda, vtcon, i_cons, mcons, & - pointlist, pointnum, factlist, nat, nsp, ityp) + ! + CALL v_of_rho_nc( rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, & + nl, ngm, gstart, nspin, g, gg, alat, omega, ehart, & + etxc, vtxc, charge, vr, lambda, vtcon, i_cons, mcons, & + pointlist, pointnum, factlist, nat, nsp, ityp ) + ! ELSE + ! CALL v_of_rho( rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, & - nrxx, nl, ngm, gstart, nspin, g, gg, alat, omega, & - ehart, etxc, vtxc, etotefield, charge, vr ) + nrxx, nl, ngm, gstart, nspin, g, gg, alat, omega, & + ehart, etxc, vtxc, etotefield, charge, vr ) + ! END IF ! IF ( ABS( charge - nelec ) / charge > 1.D-7 ) THEN @@ -299,12 +336,6 @@ SUBROUTINE extrapolate_charge( rho_order ) ! END IF ! - ! ... write potential (and rho) on file - ! - IF ( imix >= 0 ) CALL io_pot( + 1, TRIM( prefix )//'.rho', rho, nspin ) - ! - CALL io_pot( + 1, TRIM( prefix )//'.pot', vr, nspin ) - ! DEALLOCATE( work ) ! RETURN @@ -313,7 +344,7 @@ END SUBROUTINE extrapolate_charge ! ! !----------------------------------------------------------------------- -SUBROUTINE extrapolate_wfcs( wfc_order ) +SUBROUTINE extrapolate_wfcs( wfc_extr ) !----------------------------------------------------------------------- ! ! ... This routine extrapolate the wfc's after a "parallel alignment" @@ -323,7 +354,7 @@ SUBROUTINE extrapolate_wfcs( wfc_order ) USE io_global, ONLY : stdout USE kinds, ONLY : DP USE klist, ONLY : nks - USE control_flags, ONLY : isolve, alpha0, beta0, order + USE control_flags, ONLY : isolve, alpha0, beta0, wfc_order USE wvfct, ONLY : nbnd, npw, npwx, igk USE io_files, ONLY : nwordwfc, iunigk, iunwfc, iunoldwfc, & iunoldwfc2, prefix @@ -332,7 +363,7 @@ SUBROUTINE extrapolate_wfcs( wfc_order ) ! IMPLICIT NONE ! - INTEGER, INTENT(IN) :: wfc_order + INTEGER, INTENT(IN) :: wfc_extr ! INTEGER :: j, i, ik, zero_ew, lwork, info ! do-loop variables @@ -355,14 +386,7 @@ SUBROUTINE extrapolate_wfcs( wfc_order ) LOGICAL :: exst ! ! - IF (noncolin) call errore('extrapolate_wfcs', & - 'not implemented in the noncollinear case',-1) - - IF ( wfc_order == 0 .or. noncolin) THEN - ! - RETURN - ! - ELSE IF ( wfc_order == 1 ) THEN + IF ( wfc_extr == 1 ) THEN ! CALL diropn( iunoldwfc, TRIM( prefix ) // '.oldwfc', nwordwfc, exst ) ! @@ -377,17 +401,17 @@ SUBROUTINE extrapolate_wfcs( wfc_order ) ! CLOSE( UNIT = iunoldwfc, STATUS = 'KEEP' ) ! - ELSE IF ( wfc_order == 2 ) THEN + ELSE IF ( wfc_extr == 2 ) THEN ! CALL diropn( iunoldwfc, TRIM( prefix ) // '.oldwfc', nwordwfc, exst ) ! - IF ( order > 2 ) & + IF ( wfc_order > 2 ) & CALL diropn( iunoldwfc2, TRIM( prefix ) // '.old2wfc', nwordwfc, exst ) ! ALLOCATE( evcold(npwx,nbnd) ) ! WRITE( UNIT = stdout, & - FMT = '(5X,"Extrapolating wave-functions (first order) ...")' ) + FMT = '(5X,"first order wave-functions extrapolation")' ) ! lwork = 5 * nbnd ! @@ -459,13 +483,11 @@ SUBROUTINE extrapolate_wfcs( wfc_order ) ! ... extrapolate the wfc's (note that evcold contains wavefcts ! ... at (t) and evc contains wavefcts at (t-dt) ) ! - ! evc = 2.D0 * evcold - evc <= this was the previous recipe - ! - evc = evcold + alpha0 * ( evcold - evc ) + evc = 2.D0 * evcold - evc ! ! ... move the files: "old" -> "old1" and "now" -> "old" ! - IF ( order > 2 ) THEN + IF ( wfc_order > 2 ) THEN ! CALL davcio( evcold, nwordwfc, iunoldwfc, ik, - 1 ) CALL davcio( evcold, nwordwfc, iunoldwfc2, ik, + 1 ) @@ -491,12 +513,12 @@ SUBROUTINE extrapolate_wfcs( wfc_order ) DEALLOCATE( evcold ) ! CLOSE( UNIT = iunoldwfc, STATUS = 'KEEP' ) - IF ( order > 2 ) & + IF ( wfc_order > 2 ) & CLOSE( UNIT = iunoldwfc2, STATUS = 'KEEP' ) ! ELSE ! - ! ... case : wfc_order = 3 + ! ... case : wfc_extr = 3 ! CALL diropn( iunoldwfc, TRIM( prefix ) // '.oldwfc', nwordwfc, exst ) CALL diropn( iunoldwfc2, TRIM( prefix ) // '.old2wfc', nwordwfc, exst ) @@ -504,7 +526,7 @@ SUBROUTINE extrapolate_wfcs( wfc_order ) ALLOCATE( evcold(npwx,nbnd) ) ! WRITE( UNIT = stdout, & - FMT = '(5X,"Extrapolating wave-functions (second order) ...")' ) + FMT = '(5X,"second order wave-functions extrapolation")' ) ! lwork = 5 * nbnd ! @@ -574,7 +596,7 @@ SUBROUTINE extrapolate_wfcs( wfc_order ) CALL davcio( evc, nwordwfc, iunoldwfc, ik, - 1 ) ! ! ... extrapolate the wfc's, - ! ... if wfc_order == 3 use the second order extrapolation formula + ! ... if wfc_extr == 3 use the second order extrapolation formula ! ... alpha0 and beta0 are calculated in "move_ions" ! evc = ( 1 + alpha0 ) * evcold + ( beta0 - alpha0 ) * evc diff --git a/examples/example03/run_example b/examples/example03/run_example index dfa0bfa28..9d3726804 100755 --- a/examples/example03/run_example +++ b/examples/example03/run_example @@ -192,7 +192,8 @@ cat > al001.mm.in << EOF / &ions ion_dynamics='damp', - potential_extrapolation='wfc2' + pot_extrapolation='second_order' + wfc_extrapolation='second-order' / ATOMIC_SPECIES Al 26.98 Al.vbc.UPF diff --git a/examples/example17/run_example b/examples/example17/run_example index f6ca8db8b..f66f030c0 100755 --- a/examples/example17/run_example +++ b/examples/example17/run_example @@ -101,11 +101,13 @@ cat > H2+H.in << EOF mixing_beta = 0.3D0, / &IONS - num_of_images = 7, - k_max = 0.3D0, - k_min = 0.2D0, - potential_extrapolation = "wfc2", - path_thr = 0.1D0, + ds = 2.D0, + num_of_images = 7, + k_max = 0.3D0, + k_min = 0.2D0, + pot_extrapolation = "second_order", + wfc_extrapolation = "second_order", + path_thr = 0.1D0, / ATOMIC_SPECIES H 1.00794 HUSPBE.RRKJ3 @@ -160,11 +162,13 @@ cat > symmetric_H2+H.in << EOF mixing_beta = 0.3D0, / &IONS - num_of_images = 8, - k_max = 0.3D0, - k_min = 0.2D0, - potential_extrapolation = "wfc2", - path_thr = 0.2D0, + ds = 2.D0, + num_of_images = 8, + k_max = 0.3D0, + k_min = 0.2D0, + pot_extrapolation = "second_order", + wfc_extrapolation = "second_order", + path_thr = 0.2D0, / ATOMIC_SPECIES H 1.00794 HUSPBE.RRKJ3 @@ -226,14 +230,15 @@ cat > asymmetric_H2+H.in << EOF mixing_beta = 0.3D0, / &IONS - num_of_images = 8, - reset_vel = .TRUE., - ds = 1.D0, - k_max = 0.3D0, - k_min = 0.2D0, - CI_scheme = "manual", - potential_extrapolation = "wfc2", - path_thr = 0.05D0, + num_of_images = 8, + reset_vel = .TRUE., + ds = 2.D0, + k_max = 0.3D0, + k_min = 0.2D0, + CI_scheme = "manual", + pot_extrapolation = "second_order", + wfc_extrapolation = "second_order", + path_thr = 0.05D0, / CLIMBING_IMAGES 5