From 99e9255211e976ee80c83db002074b979e33ef7e Mon Sep 17 00:00:00 2001 From: Daniele Rapetti <5535617+Iximiel@users.noreply.github.com> Date: Tue, 28 May 2024 09:55:08 +0200 Subject: [PATCH] updated qe patches from higher versions --- .../PW/src/plugin_ext_forces.f90 | 13 +- patches/qespresso-7.0.diff/Modules/Makefile | 175 ++++++ .../Modules/Makefile.preplumed | 171 ++++++ patches/qespresso-7.0.diff/PW/src/forces.f90 | 486 +++++++++++++++ .../PW/src/forces.f90.preplumed | 484 +++++++++++++++ .../PW/src/plugin_ext_forces.f90 | 74 +++ .../PW/src/plugin_ext_forces.f90.preplumed | 23 + .../PW/src/plugin_initialization.f90 | 64 ++ .../src/plugin_initialization.f90.preplumed | 21 + .../qespresso-7.0.diff/PW/src/run_pwscf.f90 | 532 ++++++++++++++++ .../PW/src/run_pwscf.f90.preplumed | 518 ++++++++++++++++ patches/qespresso-7.2.diff/Modules/Makefile | 249 ++++++++ .../Modules/Makefile.preplumed | 244 ++++++++ patches/qespresso-7.2.diff/PW/src/forces.f90 | 532 ++++++++++++++++ .../PW/src/forces.f90.preplumed | 535 ++++++++++++++++ .../PW/src/plugin_ext_forces.f90 | 74 +++ .../PW/src/plugin_ext_forces.f90.preplumed | 23 + .../PW/src/plugin_initialization.f90 | 64 ++ .../src/plugin_initialization.f90.preplumed | 21 + .../qespresso-7.2.diff/PW/src/run_pwscf.f90 | 569 ++++++++++++++++++ .../PW/src/run_pwscf.f90.preplumed | 560 +++++++++++++++++ 21 files changed, 5428 insertions(+), 4 deletions(-) create mode 100644 patches/qespresso-7.0.diff/Modules/Makefile create mode 100644 patches/qespresso-7.0.diff/Modules/Makefile.preplumed create mode 100644 patches/qespresso-7.0.diff/PW/src/forces.f90 create mode 100644 patches/qespresso-7.0.diff/PW/src/forces.f90.preplumed create mode 100644 patches/qespresso-7.0.diff/PW/src/plugin_ext_forces.f90 create mode 100644 patches/qespresso-7.0.diff/PW/src/plugin_ext_forces.f90.preplumed create mode 100644 patches/qespresso-7.0.diff/PW/src/plugin_initialization.f90 create mode 100644 patches/qespresso-7.0.diff/PW/src/plugin_initialization.f90.preplumed create mode 100644 patches/qespresso-7.0.diff/PW/src/run_pwscf.f90 create mode 100644 patches/qespresso-7.0.diff/PW/src/run_pwscf.f90.preplumed create mode 100644 patches/qespresso-7.2.diff/Modules/Makefile create mode 100644 patches/qespresso-7.2.diff/Modules/Makefile.preplumed create mode 100644 patches/qespresso-7.2.diff/PW/src/forces.f90 create mode 100644 patches/qespresso-7.2.diff/PW/src/forces.f90.preplumed create mode 100644 patches/qespresso-7.2.diff/PW/src/plugin_ext_forces.f90 create mode 100644 patches/qespresso-7.2.diff/PW/src/plugin_ext_forces.f90.preplumed create mode 100644 patches/qespresso-7.2.diff/PW/src/plugin_initialization.f90 create mode 100644 patches/qespresso-7.2.diff/PW/src/plugin_initialization.f90.preplumed create mode 100644 patches/qespresso-7.2.diff/PW/src/run_pwscf.f90 create mode 100644 patches/qespresso-7.2.diff/PW/src/run_pwscf.f90.preplumed diff --git a/patches/qespresso-6.2.diff/PW/src/plugin_ext_forces.f90 b/patches/qespresso-6.2.diff/PW/src/plugin_ext_forces.f90 index ac9925473c..4df0bbaab7 100644 --- a/patches/qespresso-6.2.diff/PW/src/plugin_ext_forces.f90 +++ b/patches/qespresso-6.2.diff/PW/src/plugin_ext_forces.f90 @@ -18,19 +18,21 @@ SUBROUTINE plugin_ext_forces() USE plugin_flags ! USE cell_base, ONLY : alat, at - USE ions_base, ONLY : tau, nat,amass + USE ions_base, ONLY : tau, nat, amass, ityp USE force_mod, ONLY : force,sigma USE control_flags, ONLY : istep USE ener, ONLY : etot ! IMPLICIT NONE ! - INTEGER:: i,j + INTEGER:: i,j,ia REAL(DP) :: at_plumed(3,3) REAL(DP) :: virial(3,3) REAL(DP) :: volume REAL(DP), ALLOCATABLE :: tau_plumed(:,:) + REAL(DP) :: masses_plumed(nat) ! + masses_plumed = 0.0_DP IF(use_plumed) then IF(ionode)THEN at_plumed=alat*at; ! the cell, rescaled properly @@ -43,9 +45,12 @@ SUBROUTINE plugin_ext_forces() -at_plumed(1,2)*at_plumed(3,3)*at_plumed(2,1) & -at_plumed(1,3)*at_plumed(3,1)*at_plumed(2,2) virial=-sigma*volume - + ! the masses in QE are stored per type, see q-e//Modules/ions_base.f90 + do ia=1,nat + masses_plumed(ia)=amass(ityp(ia)) + end do CALL plumed_f_gcmd("setStep"//char(0),istep) - CALL plumed_f_gcmd("setMasses"//char(0),amass) + CALL plumed_f_gcmd("setMasses"//char(0),masses_plumed) CALL plumed_f_gcmd("setForces"//char(0),force) CALL plumed_f_gcmd("setPositions"//char(0),tau_plumed) CALL plumed_f_gcmd("setBox"//char(0),at_plumed) diff --git a/patches/qespresso-7.0.diff/Modules/Makefile b/patches/qespresso-7.0.diff/Modules/Makefile new file mode 100644 index 0000000000..36b065e4f8 --- /dev/null +++ b/patches/qespresso-7.0.diff/Modules/Makefile @@ -0,0 +1,175 @@ +#/a Makefile for Modules + +include ../make.inc + +# location of needed modules +MODFLAGS=$(BASEMOD_FLAGS) \ + $(MOD_FLAG)../ELPA/src + +# list of modules + +MODULES = \ +additional_kpoints.o \ +autopilot.o \ +basic_algebra_routines.o \ +becmod.o \ +bfgs_module.o \ +bspline.o \ +bz_form.o \ +cell_base.o \ +check_stop.o \ +command_line_options.o \ +compute_dipole.o \ +constants.o \ +constraints_module.o \ +control_flags.o \ +coulomb_vcut.o \ +dist.o \ +electrons_base.o \ +environment.o \ +fd_gradient.o \ +fft_base.o \ +fft_rho.o \ +fsockets.o \ +funct.o \ +generate_function.o \ +gradutils.o \ +gvecw.o \ +input_parameters.o \ +invmat.o \ +io_files.o \ +io_global.o \ +ions_base.o \ +kind.o \ +lmdif.o \ +mdiis.o \ +mm_dispersion.o \ +mp_bands.o \ +mp_exx.o \ +mp_global.o \ +mp_images.o \ +mp_pools.o \ +mp_wave.o \ +mp_world.o \ +noncol.o \ +open_close_input_file.o \ +parameters.o \ +parser.o \ +plugin_flags.o \ +plugin_arguments.o \ +plugin_variables.o \ +pw_dot.o \ +qmmm.o \ +random_numbers.o \ +read_cards.o \ +read_input.o \ +read_namelists.o \ +read_pseudo.o \ +recvec.o \ +recvec_subs.o \ +run_info.o \ +space_group.o \ +set_para_diag.o \ +set_signal.o \ +set_vdw_corr.o \ +setqf.o \ +timestep.o\ +tsvdw.o\ +mbdlib.o\ +version.o \ +wannier_gw.o\ +wannier_new.o \ +wavefunctions.o \ +ws_base.o \ +xc_vdW_DF.o \ +xc_rVV10.o \ +io_base.o \ +qes_types_module.o \ +qes_libs_module.o \ +qes_write_module.o \ +qes_read_module.o \ +qes_reset_module.o \ +qes_init_module.o \ +qes_bcast_module.o \ +qexsd.o \ +qexsd_copy.o \ +qexsd_init.o \ +qexsd_input.o \ +hdf5_qe.o\ +qeh5_module.o\ +fox_init_module.o \ +xsf.o \ +wyckoff.o \ +wypos.o \ +zvscal.o \ +wave_gauge.o \ +plumed.o + +# list of subroutines and functions (not modules) previously found in flib/clib + +OBJS = \ +atom_weight.o \ +capital.o \ +cryst_to_car.o \ +expint.o \ +generate_k_along_lines.o \ +has_xml.o \ +inpfile.o \ +int_to_char.o \ +latgen.o \ +linpack.o \ +matches.o \ +plot_io.o \ +radial_gradients.o \ +rgen.o \ +recips.o \ +remove_tot_torque.o \ +set_hubbard_l.o \ +set_hubbard_n.o \ +sort.o \ +trimcheck.o \ +test_input_file.o \ +date_and_tim.o \ +volume.o \ +wgauss.o \ +w0gauss.o \ +w1gauss.o \ +deviatoric.o \ +customize_signals.o \ +qmmm_aux.o \ +sockets.o \ +stack.o + +# GPU versions of modules +MODULES += \ + wavefunctions_gpu.o \ + becmod_gpu.o \ + becmod_subs_gpu.o \ + cuda_subroutines.o \ + random_numbers_gpu.o + +TLDEPS= libfox libutil libla libfft librxc + +all : libqemod.a + +## The following is needed only for lapack compiled from sources + +dlamch.o : dlamch.f + $(F77) $(FFLAGS_NOOPT) -c $< + +libqemod.a: $(MODULES) $(OBJS) + $(AR) $(ARFLAGS) $@ $? + $(RANLIB) $@ + +tldeps : + if test -n "$(TLDEPS)" ; then \ + ( cd ../.. ; $(MAKE) $(TLDEPS) || exit 1 ) ; fi + + +clean : + - /bin/rm -f *.o *.a *.d *.i *~ *_tmp.f90 *.mod *.L plumed.f90 + +plumed.f90: + cp $(PLUMED_FORTRAN) plumed.f90 + +include make.depend diff --git a/patches/qespresso-7.0.diff/Modules/Makefile.preplumed b/patches/qespresso-7.0.diff/Modules/Makefile.preplumed new file mode 100644 index 0000000000..3cbdc0ccee --- /dev/null +++ b/patches/qespresso-7.0.diff/Modules/Makefile.preplumed @@ -0,0 +1,171 @@ +#/a Makefile for Modules + +include ../make.inc + +# location of needed modules +MODFLAGS=$(BASEMOD_FLAGS) \ + $(MOD_FLAG)../ELPA/src + +# list of modules + +MODULES = \ +additional_kpoints.o \ +autopilot.o \ +basic_algebra_routines.o \ +becmod.o \ +bfgs_module.o \ +bspline.o \ +bz_form.o \ +cell_base.o \ +check_stop.o \ +command_line_options.o \ +compute_dipole.o \ +constants.o \ +constraints_module.o \ +control_flags.o \ +coulomb_vcut.o \ +dist.o \ +electrons_base.o \ +environment.o \ +fd_gradient.o \ +fft_base.o \ +fft_rho.o \ +fsockets.o \ +funct.o \ +generate_function.o \ +gradutils.o \ +gvecw.o \ +input_parameters.o \ +invmat.o \ +io_files.o \ +io_global.o \ +ions_base.o \ +kind.o \ +lmdif.o \ +mdiis.o \ +mm_dispersion.o \ +mp_bands.o \ +mp_exx.o \ +mp_global.o \ +mp_images.o \ +mp_pools.o \ +mp_wave.o \ +mp_world.o \ +noncol.o \ +open_close_input_file.o \ +parameters.o \ +parser.o \ +plugin_flags.o \ +plugin_arguments.o \ +plugin_variables.o \ +pw_dot.o \ +qmmm.o \ +random_numbers.o \ +read_cards.o \ +read_input.o \ +read_namelists.o \ +read_pseudo.o \ +recvec.o \ +recvec_subs.o \ +run_info.o \ +space_group.o \ +set_para_diag.o \ +set_signal.o \ +set_vdw_corr.o \ +setqf.o \ +timestep.o\ +tsvdw.o\ +mbdlib.o\ +version.o \ +wannier_gw.o\ +wannier_new.o \ +wavefunctions.o \ +ws_base.o \ +xc_vdW_DF.o \ +xc_rVV10.o \ +io_base.o \ +qes_types_module.o \ +qes_libs_module.o \ +qes_write_module.o \ +qes_read_module.o \ +qes_reset_module.o \ +qes_init_module.o \ +qes_bcast_module.o \ +qexsd.o \ +qexsd_copy.o \ +qexsd_init.o \ +qexsd_input.o \ +hdf5_qe.o\ +qeh5_module.o\ +fox_init_module.o \ +xsf.o \ +wyckoff.o \ +wypos.o \ +zvscal.o \ +wave_gauge.o + +# list of subroutines and functions (not modules) previously found in flib/clib + +OBJS = \ +atom_weight.o \ +capital.o \ +cryst_to_car.o \ +expint.o \ +generate_k_along_lines.o \ +has_xml.o \ +inpfile.o \ +int_to_char.o \ +latgen.o \ +linpack.o \ +matches.o \ +plot_io.o \ +radial_gradients.o \ +rgen.o \ +recips.o \ +remove_tot_torque.o \ +set_hubbard_l.o \ +set_hubbard_n.o \ +sort.o \ +trimcheck.o \ +test_input_file.o \ +date_and_tim.o \ +volume.o \ +wgauss.o \ +w0gauss.o \ +w1gauss.o \ +deviatoric.o \ +customize_signals.o \ +qmmm_aux.o \ +sockets.o \ +stack.o + +# GPU versions of modules +MODULES += \ + wavefunctions_gpu.o \ + becmod_gpu.o \ + becmod_subs_gpu.o \ + cuda_subroutines.o \ + random_numbers_gpu.o + +TLDEPS= libfox libutil libla libfft librxc + +all : libqemod.a + +## The following is needed only for lapack compiled from sources + +dlamch.o : dlamch.f + $(F77) $(FFLAGS_NOOPT) -c $< + +libqemod.a: $(MODULES) $(OBJS) + $(AR) $(ARFLAGS) $@ $? + $(RANLIB) $@ + +tldeps : + if test -n "$(TLDEPS)" ; then \ + ( cd ../.. ; $(MAKE) $(TLDEPS) || exit 1 ) ; fi + + +clean : + - /bin/rm -f *.o *.a *.d *.i *~ *_tmp.f90 *.mod *.L + +include make.depend diff --git a/patches/qespresso-7.0.diff/PW/src/forces.f90 b/patches/qespresso-7.0.diff/PW/src/forces.f90 new file mode 100644 index 0000000000..4a9d2e7d6f --- /dev/null +++ b/patches/qespresso-7.0.diff/PW/src/forces.f90 @@ -0,0 +1,486 @@ +! +! Copyright (C) 2001-2011 Quantum ESPRESSO 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 forces() + !---------------------------------------------------------------------------- + !! This routine is a driver routine which computes the forces + !! acting on the atoms. The complete expression of the forces + !! contains many parts which are computed by different routines: + ! + !! - force_lc: local potential contribution + !! - force_us: non-local potential contribution + !! - (esm_)force_ew: (ESM) electrostatic ewald term + !! - force_cc: nonlinear core correction contribution + !! - force_corr: correction term for incomplete self-consistency + !! - force_hub: contribution due to the Hubbard term; + !! - force_london: Grimme DFT+D dispersion forces + !! - force_d3: Grimme-D3 (DFT-D3) dispersion forces + !! - force_xdm: XDM dispersion forces + !! - more terms from external electric fields, Martyna-Tuckerman, etc. + ! + USE kinds, ONLY : DP + USE io_global, ONLY : stdout + USE cell_base, ONLY : at, bg, alat, omega + USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau, zv, amass, extfor, atm + USE fft_base, ONLY : dfftp + USE gvect, ONLY : ngm, gstart, ngl, igtongl, igtongl_d, g, gg, & + g_d, gcutm + USE lsda_mod, ONLY : nspin + USE symme, ONLY : symvector + USE vlocal, ONLY : strf, vloc + USE force_mod, ONLY : force, sumfor + USE scf, ONLY : rho + USE ions_base, ONLY : if_pos + USE ldaU, ONLY : lda_plus_u, U_projection + USE extfield, ONLY : tefield, forcefield, gate, forcegate, relaxz + USE control_flags, ONLY : gamma_only, remove_rigid_rot, textfor, & + iverbosity, llondon, ldftd3, lxdm, ts_vdw, & + mbd_vdw, lforce => tprnfor + USE plugin_flags + USE bp, ONLY : lelfield, gdir, l3dstring, efield_cart, & + efield_cry,efield + USE uspp, ONLY : okvan + USE martyna_tuckerman, ONLY : do_comp_mt, wg_corr_force + USE london_module, ONLY : force_london + USE dftd3_api, ONLY : get_atomic_number, dftd3_calc + USE dftd3_qe, ONLY : dftd3_pbc_gdisp, dftd3 + + USE xdm_module, ONLY : force_xdm + USE tsvdw_module, ONLY : FtsvdW + USE libmbd_interface, ONLY : FmbdvdW + USE esm, ONLY : do_comp_esm, esm_bc, esm_force_ew + USE qmmm, ONLY : qmmm_mode + ! + USE control_flags, ONLY : use_gpu + USE device_fbuff_m, ONLY : dev_buf + USE device_memcpy_m, ONLY : dev_memcpy + ! + IMPLICIT NONE + ! + REAL(DP), ALLOCATABLE :: forcenl(:,:), & + forcelc(:,:), & + forcecc(:,:), & + forceion(:,:), & + force_disp(:,:), & + force_d3(:,:), & + force_disp_xdm(:,:), & + force_mt(:,:), & + forcescc(:,:), & + forces_bp_efield(:,:),& + forceh(:,:) + ! nonlocal, local, core-correction, ewald, scf correction terms, and hubbard + ! + ! aux is used to store a possible additional density + ! now defined in real space + ! + COMPLEX(DP), ALLOCATABLE :: auxg(:), auxr(:) + ! + REAL(DP) :: sumscf, sum_mm + REAL(DP), PARAMETER :: eps = 1.e-12_dp + INTEGER :: ipol, na + ! counter on polarization + ! counter on atoms + ! + REAL(DP) :: latvecs(3,3) + INTEGER :: atnum(1:nat) + REAL(DP) :: stress_dftd3(3,3) + ! + ! TODO: get rid of this !!!! Use standard method for duplicated global data + REAL(DP), POINTER :: vloc_d (:, :) + INTEGER :: ierr +#if defined(__CUDA) + attributes(DEVICE) :: vloc_d +#endif + ! + force(:,:) = 0.D0 + ! + ! Early return if all forces to be set to zero + ! + IF ( ALL( if_pos == 0 ) ) RETURN + ! + CALL start_clock( 'forces' ) + ! Cleanup scratch space used in previous SCF iterations. + ! This will reduce memory footprint. + CALL dev_buf%reinit(ierr) + IF (ierr .ne. 0) CALL infomsg('forces', 'Cannot reset GPU buffers! Some buffers still locked.') + ! + ! + ALLOCATE( forcenl(3,nat), forcelc(3,nat), forcecc(3,nat), & + forceh(3,nat), forceion(3,nat), forcescc(3,nat) ) + ! + forcescc(:,:) = 0.D0 + forceh(:,:) = 0.D0 + ! + ! ... The nonlocal contribution is computed here + ! + call start_clock('frc_us') + IF (.not. use_gpu) CALL force_us( forcenl ) + IF ( use_gpu) CALL force_us_gpu( forcenl ) + call stop_clock('frc_us') + ! + ! ... The local contribution + ! + CALL start_clock('frc_lc') + IF (.not. use_gpu) & ! On the CPU + CALL force_lc( nat, tau, ityp, alat, omega, ngm, ngl, igtongl, & + g, rho%of_r(:,1), dfftp%nl, gstart, gamma_only, vloc, & + forcelc ) + IF ( use_gpu) THEN ! On the GPU + ! move these data to the GPU + CALL dev_buf%lock_buffer(vloc_d, (/ ngl, ntyp /) , ierr) + IF (ierr /= 0) CALL errore( 'forces', 'cannot allocate buffers', -1 ) + CALL dev_memcpy(vloc_d, vloc) + CALL force_lc_gpu( nat, tau, ityp, alat, omega, ngm, ngl, igtongl_d, & + g_d, rho%of_r(:,1), dfftp%nl_d, gstart, gamma_only, vloc_d, & + forcelc ) + CALL dev_buf%release_buffer(vloc_d, ierr) + END IF + call stop_clock('frc_lc') + ! + ! ... The NLCC contribution + ! + call start_clock('frc_cc') + IF (.not. use_gpu) CALL force_cc( forcecc ) + IF ( use_gpu) CALL force_cc_gpu( forcecc ) + ! + call stop_clock('frc_cc') + + ! ... The Hubbard contribution + ! (included by force_us if using beta as local projectors) + ! + IF (.not. use_gpu) THEN + IF ( lda_plus_u .AND. U_projection.NE.'pseudo' ) CALL force_hub( forceh ) + ELSE + IF ( lda_plus_u .AND. U_projection.NE.'pseudo' ) CALL force_hub_gpu( forceh ) + ENDIF + ! + ! ... The ionic contribution is computed here + ! + IF( do_comp_esm ) THEN + CALL esm_force_ew( forceion ) + ELSE + CALL force_ew( alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, & + gg, ngm, gstart, gamma_only, gcutm, strf, forceion ) + ENDIF + ! + ! ... the semi-empirical dispersion correction + ! + IF ( llondon ) THEN + ! + ALLOCATE( force_disp(3,nat) ) + force_disp(:,:) = 0.0_DP + force_disp = force_london( alat , nat , ityp , at , bg , tau ) + ! + ENDIF + ! + ! ... The Grimme-D3 dispersion correction + ! + IF ( ldftd3 ) THEN + ! + CALL start_clock('force_dftd3') + ALLOCATE( force_d3(3, nat) ) + force_d3(:,:) = 0.0_DP + latvecs(:,:) = at(:,:)*alat + tau(:,:) = tau(:,:)*alat + atnum(:) = get_atomic_number(atm(ityp(:))) + CALL dftd3_pbc_gdisp( dftd3, tau, atnum, latvecs, & + force_d3, stress_dftd3 ) + force_d3 = -2.d0*force_d3 + tau(:,:) = tau(:,:)/alat + CALL stop_clock('force_dftd3') + ENDIF + ! + ! + IF (lxdm) THEN + ALLOCATE( force_disp_xdm(3,nat) ) + force_disp_xdm = 0._dp + force_disp_xdm = force_xdm(nat) + ENDIF + ! + ! ... The SCF contribution + ! + call start_clock('frc_scc') + ! Cleanup scratch space again, next subroutines uses a lot of memory. + ! In an ideal world this should be done only if really needed (TODO). + CALL dev_buf%reinit(ierr) + IF (ierr .ne. 0) CALL errore('forces', 'Cannot reset GPU buffers! Buffers still locked: ', abs(ierr)) + ! + IF ( .not. use_gpu ) CALL force_corr( forcescc ) + IF ( use_gpu ) CALL force_corr_gpu( forcescc ) + call stop_clock('frc_scc') + ! + IF (do_comp_mt) THEN + ! + ALLOCATE( force_mt(3,nat) ) + CALL wg_corr_force( .TRUE., omega, nat, ntyp, ityp, ngm, g, tau, zv, strf, & + rho%of_g(:,1), force_mt ) + ENDIF + ! + ! ... call void routine for user define/ plugin patches on internal forces + ! + CALL plugin_int_forces() + ! + ! ... Berry's phase electric field terms + ! + IF (lelfield) THEN + ALLOCATE( forces_bp_efield(3,nat) ) + forces_bp_efield(:,:) = 0.d0 + IF (.NOT.l3dstring) THEN + IF (okvan) CALL forces_us_efield( forces_bp_efield, gdir, efield ) + CALL forces_ion_efield( forces_bp_efield, gdir, efield ) + ELSE + IF (okvan) THEN + DO ipol = 1, 3 + CALL forces_us_efield( forces_bp_efield, ipol, efield_cry(ipol) ) + ENDDO + ENDIF + DO ipol = 1, 3 + CALL forces_ion_efield( forces_bp_efield, ipol, efield_cart(ipol) ) + ENDDO + ENDIF + ENDIF + ! + ! ... here we sum all the contributions and compute the total force acting + ! ... on the crystal + ! + DO ipol = 1, 3 + ! + sumfor = 0.D0 + ! + DO na = 1, nat + ! + force(ipol,na) = force(ipol,na) + & + forcenl(ipol,na) + & + forceion(ipol,na) + & + forcelc(ipol,na) + & + forcecc(ipol,na) + & + forceh(ipol,na) + & + forcescc(ipol,na) + ! + IF ( llondon ) force(ipol,na) = force(ipol,na) + force_disp(ipol,na) + IF ( ldftd3 ) force(ipol,na) = force(ipol,na) + force_d3(ipol,na) + IF ( lxdm ) force(ipol,na) = force(ipol,na) + force_disp_xdm(ipol,na) + ! factor 2 converts from Ha to Ry a.u. + ! the IF condition is to avoid double counting + IF ( mbd_vdw ) THEN + force(ipol, na) = force(ipol, na) + 2.0_dp*FmbdvdW(ipol, na) + ELSE IF ( ts_vdw ) THEN + force(ipol, na) = force(ipol, na) + 2.0_dp*FtsvdW(ipol, na) + ENDIF + IF ( tefield ) force(ipol,na) = force(ipol,na) + forcefield(ipol,na) + IF ( gate ) force(ipol,na) = force(ipol,na) + forcegate(ipol,na) ! TB + IF (lelfield) force(ipol,na) = force(ipol,na) + forces_bp_efield(ipol,na) + IF (do_comp_mt) force(ipol,na) = force(ipol,na) + force_mt(ipol,na) + ! + sumfor = sumfor + force(ipol,na) + ! + ENDDO + ! + !TB + IF ((gate.AND.relaxz).AND.(ipol==3)) WRITE( stdout, '("Total force in z direction = 0 disabled")') + ! + IF ( (do_comp_esm .AND. ( esm_bc /= 'pbc' )).OR.(gate.AND.relaxz) ) THEN + ! + ! ... impose total force along xy = 0 + ! + DO na = 1, nat + IF ( ipol /= 3) force(ipol,na) = force(ipol,na) & + - sumfor / DBLE( nat ) + ENDDO + ! + ELSEIF ( qmmm_mode < 0 ) THEN + ! + ! ... impose total force = 0 except in a QM-MM calculation + ! + DO na = 1, nat + force(ipol,na) = force(ipol,na) - sumfor / DBLE( nat ) + ENDDO + ! + ENDIF + ! + ENDDO + ! + ! ... resymmetrize (should not be needed, but ...) + ! + CALL symvector( nat, force ) + ! + IF ( remove_rigid_rot ) & + CALL remove_tot_torque( nat, tau, amass(ityp(:)), force ) + ! + IF( textfor ) force(:,:) = force(:,:) + extfor(:,:) + ! + ! ... call void routine for user define/ plugin patches on external forces + ! + ! plugin should be called after stress has been computed + ! look in run_pwscf.f90 + !CALL plugin_ext_forces() + ! + ! ... write on output the forces + ! + WRITE( stdout, '(/,5x,"Forces acting on atoms (cartesian axes, Ry/au):", / )') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), force(:,na) + ENDDO + ! + ! ... forces on fixed coordinates are set to zero ( C.S. 15/10/2003 ) + ! + force(:,:) = force(:,:) * DBLE( if_pos ) + forcescc(:,:) = forcescc(:,:) * DBLE( if_pos ) + ! + IF ( iverbosity > 0 ) THEN + IF ( do_comp_mt ) THEN + WRITE( stdout, '(5x,"The Martyna-Tuckerman correction term to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( force_mt(ipol,na), ipol = 1, 3 ) + ENDDO + END IF + ! + WRITE( stdout, '(5x,"The non-local contrib. to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcenl(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The ionic contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forceion(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The local contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcelc(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The core correction contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcecc(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The Hubbard contrib. to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forceh(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The SCF correction term to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcescc(ipol,na), ipol = 1, 3 ) + ENDDO + ! + IF ( llondon) THEN + WRITE( stdout, '(/,5x,"Dispersion contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_disp(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + IF ( ldftd3 ) THEN + WRITE( stdout, '(/,5x,"DFT-D3 dispersion contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_d3(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + IF (lxdm) THEN + WRITE( stdout, '(/,5x,"XDM contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_disp_xdm(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + ! again, as above, if condition is to avoid redundant printing + IF ( mbd_vdw ) THEN + WRITE( stdout, '(/,5x, "MBD contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (2.0d0*FmbdvdW(ipol, na), ipol = 1, 3) + ENDDO + ELSE IF ( ts_vdw ) THEN + WRITE( stdout, '(/,5x, "TS-VDW contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (2.0d0*FtsvdW(ipol, na), ipol = 1, 3) + ENDDO + ENDIF + + ! + ! TB gate forces + IF ( gate ) THEN + WRITE( stdout, '(/,5x,"Gate contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (forcegate(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + END IF + ! + sumfor = 0.D0 + sumscf = 0.D0 + ! + DO na = 1, nat + ! + sumfor = sumfor + force(1,na)**2 + force(2,na)**2 + force(3,na)**2 + sumscf = sumscf + forcescc(1,na)**2 + forcescc(2,na)**2+ forcescc(3,na)**2 + ! + ENDDO + ! + sumfor = SQRT( sumfor ) + sumscf = SQRT( sumscf ) + ! + WRITE( stdout, '(/5x,"Total force = ",F12.6,5X, & + & "Total SCF correction = ",F12.6)') sumfor, sumscf + ! + IF ( llondon .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_disp(1,na)**2 + force_disp(2,na)**2 + force_disp(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "Total Dispersion Force = ",F12.6)') sum_mm + ! + END IF + ! + IF ( ldftd3 .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_d3(1,na)**2 + force_d3(2,na)**2 + force_d3(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "DFT-D3 dispersion Force = ",F12.6)') sum_mm + ! + END IF + ! + IF ( lxdm .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_disp_xdm(1,na)**2 + force_disp_xdm(2,na)**2 + force_disp_xdm(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "Total XDM Force = ",F12.6)') sum_mm + ! + END IF + ! + DEALLOCATE( forcenl, forcelc, forcecc, forceh, forceion, forcescc ) + IF ( llondon ) DEALLOCATE( force_disp ) + IF ( ldftd3 ) DEALLOCATE( force_d3 ) + IF ( lxdm ) DEALLOCATE( force_disp_xdm ) + IF ( lelfield ) DEALLOCATE( forces_bp_efield ) + IF(ALLOCATED(force_mt)) DEALLOCATE( force_mt ) + ! + ! FIXME: what is the following line good for? + ! + lforce = .TRUE. + ! + CALL stop_clock( 'forces' ) + ! + IF ( ( sumfor < 10.D0*sumscf ) .AND. ( sumfor > nat*eps ) ) & + WRITE( stdout,'(5x,"SCF correction compared to forces is large: ", & + & "reduce conv_thr to get better values")') + + RETURN + ! +9035 FORMAT(5X,'atom ',I4,' type ',I2,' force = ',3F14.8) + ! +END SUBROUTINE forces diff --git a/patches/qespresso-7.0.diff/PW/src/forces.f90.preplumed b/patches/qespresso-7.0.diff/PW/src/forces.f90.preplumed new file mode 100644 index 0000000000..11c1928181 --- /dev/null +++ b/patches/qespresso-7.0.diff/PW/src/forces.f90.preplumed @@ -0,0 +1,484 @@ +! +! Copyright (C) 2001-2011 Quantum ESPRESSO 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 forces() + !---------------------------------------------------------------------------- + !! This routine is a driver routine which computes the forces + !! acting on the atoms. The complete expression of the forces + !! contains many parts which are computed by different routines: + ! + !! - force_lc: local potential contribution + !! - force_us: non-local potential contribution + !! - (esm_)force_ew: (ESM) electrostatic ewald term + !! - force_cc: nonlinear core correction contribution + !! - force_corr: correction term for incomplete self-consistency + !! - force_hub: contribution due to the Hubbard term; + !! - force_london: Grimme DFT+D dispersion forces + !! - force_d3: Grimme-D3 (DFT-D3) dispersion forces + !! - force_xdm: XDM dispersion forces + !! - more terms from external electric fields, Martyna-Tuckerman, etc. + ! + USE kinds, ONLY : DP + USE io_global, ONLY : stdout + USE cell_base, ONLY : at, bg, alat, omega + USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau, zv, amass, extfor, atm + USE fft_base, ONLY : dfftp + USE gvect, ONLY : ngm, gstart, ngl, igtongl, igtongl_d, g, gg, & + g_d, gcutm + USE lsda_mod, ONLY : nspin + USE symme, ONLY : symvector + USE vlocal, ONLY : strf, vloc + USE force_mod, ONLY : force, sumfor + USE scf, ONLY : rho + USE ions_base, ONLY : if_pos + USE ldaU, ONLY : lda_plus_u, U_projection + USE extfield, ONLY : tefield, forcefield, gate, forcegate, relaxz + USE control_flags, ONLY : gamma_only, remove_rigid_rot, textfor, & + iverbosity, llondon, ldftd3, lxdm, ts_vdw, & + mbd_vdw, lforce => tprnfor + USE plugin_flags + USE bp, ONLY : lelfield, gdir, l3dstring, efield_cart, & + efield_cry,efield + USE uspp, ONLY : okvan + USE martyna_tuckerman, ONLY : do_comp_mt, wg_corr_force + USE london_module, ONLY : force_london + USE dftd3_api, ONLY : get_atomic_number, dftd3_calc + USE dftd3_qe, ONLY : dftd3_pbc_gdisp, dftd3 + + USE xdm_module, ONLY : force_xdm + USE tsvdw_module, ONLY : FtsvdW + USE libmbd_interface, ONLY : FmbdvdW + USE esm, ONLY : do_comp_esm, esm_bc, esm_force_ew + USE qmmm, ONLY : qmmm_mode + ! + USE control_flags, ONLY : use_gpu + USE device_fbuff_m, ONLY : dev_buf + USE device_memcpy_m, ONLY : dev_memcpy + ! + IMPLICIT NONE + ! + REAL(DP), ALLOCATABLE :: forcenl(:,:), & + forcelc(:,:), & + forcecc(:,:), & + forceion(:,:), & + force_disp(:,:), & + force_d3(:,:), & + force_disp_xdm(:,:), & + force_mt(:,:), & + forcescc(:,:), & + forces_bp_efield(:,:),& + forceh(:,:) + ! nonlocal, local, core-correction, ewald, scf correction terms, and hubbard + ! + ! aux is used to store a possible additional density + ! now defined in real space + ! + COMPLEX(DP), ALLOCATABLE :: auxg(:), auxr(:) + ! + REAL(DP) :: sumscf, sum_mm + REAL(DP), PARAMETER :: eps = 1.e-12_dp + INTEGER :: ipol, na + ! counter on polarization + ! counter on atoms + ! + REAL(DP) :: latvecs(3,3) + INTEGER :: atnum(1:nat) + REAL(DP) :: stress_dftd3(3,3) + ! + ! TODO: get rid of this !!!! Use standard method for duplicated global data + REAL(DP), POINTER :: vloc_d (:, :) + INTEGER :: ierr +#if defined(__CUDA) + attributes(DEVICE) :: vloc_d +#endif + ! + force(:,:) = 0.D0 + ! + ! Early return if all forces to be set to zero + ! + IF ( ALL( if_pos == 0 ) ) RETURN + ! + CALL start_clock( 'forces' ) + ! Cleanup scratch space used in previous SCF iterations. + ! This will reduce memory footprint. + CALL dev_buf%reinit(ierr) + IF (ierr .ne. 0) CALL infomsg('forces', 'Cannot reset GPU buffers! Some buffers still locked.') + ! + ! + ALLOCATE( forcenl(3,nat), forcelc(3,nat), forcecc(3,nat), & + forceh(3,nat), forceion(3,nat), forcescc(3,nat) ) + ! + forcescc(:,:) = 0.D0 + forceh(:,:) = 0.D0 + ! + ! ... The nonlocal contribution is computed here + ! + call start_clock('frc_us') + IF (.not. use_gpu) CALL force_us( forcenl ) + IF ( use_gpu) CALL force_us_gpu( forcenl ) + call stop_clock('frc_us') + ! + ! ... The local contribution + ! + CALL start_clock('frc_lc') + IF (.not. use_gpu) & ! On the CPU + CALL force_lc( nat, tau, ityp, alat, omega, ngm, ngl, igtongl, & + g, rho%of_r(:,1), dfftp%nl, gstart, gamma_only, vloc, & + forcelc ) + IF ( use_gpu) THEN ! On the GPU + ! move these data to the GPU + CALL dev_buf%lock_buffer(vloc_d, (/ ngl, ntyp /) , ierr) + IF (ierr /= 0) CALL errore( 'forces', 'cannot allocate buffers', -1 ) + CALL dev_memcpy(vloc_d, vloc) + CALL force_lc_gpu( nat, tau, ityp, alat, omega, ngm, ngl, igtongl_d, & + g_d, rho%of_r(:,1), dfftp%nl_d, gstart, gamma_only, vloc_d, & + forcelc ) + CALL dev_buf%release_buffer(vloc_d, ierr) + END IF + call stop_clock('frc_lc') + ! + ! ... The NLCC contribution + ! + call start_clock('frc_cc') + IF (.not. use_gpu) CALL force_cc( forcecc ) + IF ( use_gpu) CALL force_cc_gpu( forcecc ) + ! + call stop_clock('frc_cc') + + ! ... The Hubbard contribution + ! (included by force_us if using beta as local projectors) + ! + IF (.not. use_gpu) THEN + IF ( lda_plus_u .AND. U_projection.NE.'pseudo' ) CALL force_hub( forceh ) + ELSE + IF ( lda_plus_u .AND. U_projection.NE.'pseudo' ) CALL force_hub_gpu( forceh ) + ENDIF + ! + ! ... The ionic contribution is computed here + ! + IF( do_comp_esm ) THEN + CALL esm_force_ew( forceion ) + ELSE + CALL force_ew( alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, & + gg, ngm, gstart, gamma_only, gcutm, strf, forceion ) + ENDIF + ! + ! ... the semi-empirical dispersion correction + ! + IF ( llondon ) THEN + ! + ALLOCATE( force_disp(3,nat) ) + force_disp(:,:) = 0.0_DP + force_disp = force_london( alat , nat , ityp , at , bg , tau ) + ! + ENDIF + ! + ! ... The Grimme-D3 dispersion correction + ! + IF ( ldftd3 ) THEN + ! + CALL start_clock('force_dftd3') + ALLOCATE( force_d3(3, nat) ) + force_d3(:,:) = 0.0_DP + latvecs(:,:) = at(:,:)*alat + tau(:,:) = tau(:,:)*alat + atnum(:) = get_atomic_number(atm(ityp(:))) + CALL dftd3_pbc_gdisp( dftd3, tau, atnum, latvecs, & + force_d3, stress_dftd3 ) + force_d3 = -2.d0*force_d3 + tau(:,:) = tau(:,:)/alat + CALL stop_clock('force_dftd3') + ENDIF + ! + ! + IF (lxdm) THEN + ALLOCATE( force_disp_xdm(3,nat) ) + force_disp_xdm = 0._dp + force_disp_xdm = force_xdm(nat) + ENDIF + ! + ! ... The SCF contribution + ! + call start_clock('frc_scc') + ! Cleanup scratch space again, next subroutines uses a lot of memory. + ! In an ideal world this should be done only if really needed (TODO). + CALL dev_buf%reinit(ierr) + IF (ierr .ne. 0) CALL errore('forces', 'Cannot reset GPU buffers! Buffers still locked: ', abs(ierr)) + ! + IF ( .not. use_gpu ) CALL force_corr( forcescc ) + IF ( use_gpu ) CALL force_corr_gpu( forcescc ) + call stop_clock('frc_scc') + ! + IF (do_comp_mt) THEN + ! + ALLOCATE( force_mt(3,nat) ) + CALL wg_corr_force( .TRUE., omega, nat, ntyp, ityp, ngm, g, tau, zv, strf, & + rho%of_g(:,1), force_mt ) + ENDIF + ! + ! ... call void routine for user define/ plugin patches on internal forces + ! + CALL plugin_int_forces() + ! + ! ... Berry's phase electric field terms + ! + IF (lelfield) THEN + ALLOCATE( forces_bp_efield(3,nat) ) + forces_bp_efield(:,:) = 0.d0 + IF (.NOT.l3dstring) THEN + IF (okvan) CALL forces_us_efield( forces_bp_efield, gdir, efield ) + CALL forces_ion_efield( forces_bp_efield, gdir, efield ) + ELSE + IF (okvan) THEN + DO ipol = 1, 3 + CALL forces_us_efield( forces_bp_efield, ipol, efield_cry(ipol) ) + ENDDO + ENDIF + DO ipol = 1, 3 + CALL forces_ion_efield( forces_bp_efield, ipol, efield_cart(ipol) ) + ENDDO + ENDIF + ENDIF + ! + ! ... here we sum all the contributions and compute the total force acting + ! ... on the crystal + ! + DO ipol = 1, 3 + ! + sumfor = 0.D0 + ! + DO na = 1, nat + ! + force(ipol,na) = force(ipol,na) + & + forcenl(ipol,na) + & + forceion(ipol,na) + & + forcelc(ipol,na) + & + forcecc(ipol,na) + & + forceh(ipol,na) + & + forcescc(ipol,na) + ! + IF ( llondon ) force(ipol,na) = force(ipol,na) + force_disp(ipol,na) + IF ( ldftd3 ) force(ipol,na) = force(ipol,na) + force_d3(ipol,na) + IF ( lxdm ) force(ipol,na) = force(ipol,na) + force_disp_xdm(ipol,na) + ! factor 2 converts from Ha to Ry a.u. + ! the IF condition is to avoid double counting + IF ( mbd_vdw ) THEN + force(ipol, na) = force(ipol, na) + 2.0_dp*FmbdvdW(ipol, na) + ELSE IF ( ts_vdw ) THEN + force(ipol, na) = force(ipol, na) + 2.0_dp*FtsvdW(ipol, na) + ENDIF + IF ( tefield ) force(ipol,na) = force(ipol,na) + forcefield(ipol,na) + IF ( gate ) force(ipol,na) = force(ipol,na) + forcegate(ipol,na) ! TB + IF (lelfield) force(ipol,na) = force(ipol,na) + forces_bp_efield(ipol,na) + IF (do_comp_mt) force(ipol,na) = force(ipol,na) + force_mt(ipol,na) + ! + sumfor = sumfor + force(ipol,na) + ! + ENDDO + ! + !TB + IF ((gate.AND.relaxz).AND.(ipol==3)) WRITE( stdout, '("Total force in z direction = 0 disabled")') + ! + IF ( (do_comp_esm .AND. ( esm_bc /= 'pbc' )).OR.(gate.AND.relaxz) ) THEN + ! + ! ... impose total force along xy = 0 + ! + DO na = 1, nat + IF ( ipol /= 3) force(ipol,na) = force(ipol,na) & + - sumfor / DBLE( nat ) + ENDDO + ! + ELSEIF ( qmmm_mode < 0 ) THEN + ! + ! ... impose total force = 0 except in a QM-MM calculation + ! + DO na = 1, nat + force(ipol,na) = force(ipol,na) - sumfor / DBLE( nat ) + ENDDO + ! + ENDIF + ! + ENDDO + ! + ! ... resymmetrize (should not be needed, but ...) + ! + CALL symvector( nat, force ) + ! + IF ( remove_rigid_rot ) & + CALL remove_tot_torque( nat, tau, amass(ityp(:)), force ) + ! + IF( textfor ) force(:,:) = force(:,:) + extfor(:,:) + ! + ! ... call void routine for user define/ plugin patches on external forces + ! + CALL plugin_ext_forces() + ! + ! ... write on output the forces + ! + WRITE( stdout, '(/,5x,"Forces acting on atoms (cartesian axes, Ry/au):", / )') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), force(:,na) + ENDDO + ! + ! ... forces on fixed coordinates are set to zero ( C.S. 15/10/2003 ) + ! + force(:,:) = force(:,:) * DBLE( if_pos ) + forcescc(:,:) = forcescc(:,:) * DBLE( if_pos ) + ! + IF ( iverbosity > 0 ) THEN + IF ( do_comp_mt ) THEN + WRITE( stdout, '(5x,"The Martyna-Tuckerman correction term to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( force_mt(ipol,na), ipol = 1, 3 ) + ENDDO + END IF + ! + WRITE( stdout, '(5x,"The non-local contrib. to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcenl(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The ionic contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forceion(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The local contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcelc(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The core correction contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcecc(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The Hubbard contrib. to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forceh(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The SCF correction term to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcescc(ipol,na), ipol = 1, 3 ) + ENDDO + ! + IF ( llondon) THEN + WRITE( stdout, '(/,5x,"Dispersion contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_disp(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + IF ( ldftd3 ) THEN + WRITE( stdout, '(/,5x,"DFT-D3 dispersion contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_d3(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + IF (lxdm) THEN + WRITE( stdout, '(/,5x,"XDM contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_disp_xdm(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + ! again, as above, if condition is to avoid redundant printing + IF ( mbd_vdw ) THEN + WRITE( stdout, '(/,5x, "MBD contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (2.0d0*FmbdvdW(ipol, na), ipol = 1, 3) + ENDDO + ELSE IF ( ts_vdw ) THEN + WRITE( stdout, '(/,5x, "TS-VDW contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (2.0d0*FtsvdW(ipol, na), ipol = 1, 3) + ENDDO + ENDIF + + ! + ! TB gate forces + IF ( gate ) THEN + WRITE( stdout, '(/,5x,"Gate contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (forcegate(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + END IF + ! + sumfor = 0.D0 + sumscf = 0.D0 + ! + DO na = 1, nat + ! + sumfor = sumfor + force(1,na)**2 + force(2,na)**2 + force(3,na)**2 + sumscf = sumscf + forcescc(1,na)**2 + forcescc(2,na)**2+ forcescc(3,na)**2 + ! + ENDDO + ! + sumfor = SQRT( sumfor ) + sumscf = SQRT( sumscf ) + ! + WRITE( stdout, '(/5x,"Total force = ",F12.6,5X, & + & "Total SCF correction = ",F12.6)') sumfor, sumscf + ! + IF ( llondon .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_disp(1,na)**2 + force_disp(2,na)**2 + force_disp(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "Total Dispersion Force = ",F12.6)') sum_mm + ! + END IF + ! + IF ( ldftd3 .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_d3(1,na)**2 + force_d3(2,na)**2 + force_d3(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "DFT-D3 dispersion Force = ",F12.6)') sum_mm + ! + END IF + ! + IF ( lxdm .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_disp_xdm(1,na)**2 + force_disp_xdm(2,na)**2 + force_disp_xdm(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "Total XDM Force = ",F12.6)') sum_mm + ! + END IF + ! + DEALLOCATE( forcenl, forcelc, forcecc, forceh, forceion, forcescc ) + IF ( llondon ) DEALLOCATE( force_disp ) + IF ( ldftd3 ) DEALLOCATE( force_d3 ) + IF ( lxdm ) DEALLOCATE( force_disp_xdm ) + IF ( lelfield ) DEALLOCATE( forces_bp_efield ) + IF(ALLOCATED(force_mt)) DEALLOCATE( force_mt ) + ! + ! FIXME: what is the following line good for? + ! + lforce = .TRUE. + ! + CALL stop_clock( 'forces' ) + ! + IF ( ( sumfor < 10.D0*sumscf ) .AND. ( sumfor > nat*eps ) ) & + WRITE( stdout,'(5x,"SCF correction compared to forces is large: ", & + & "reduce conv_thr to get better values")') + + RETURN + ! +9035 FORMAT(5X,'atom ',I4,' type ',I2,' force = ',3F14.8) + ! +END SUBROUTINE forces diff --git a/patches/qespresso-7.0.diff/PW/src/plugin_ext_forces.f90 b/patches/qespresso-7.0.diff/PW/src/plugin_ext_forces.f90 new file mode 100644 index 0000000000..5d4762fb42 --- /dev/null +++ b/patches/qespresso-7.0.diff/PW/src/plugin_ext_forces.f90 @@ -0,0 +1,74 @@ +! +! Copyright (C) 2001-2009 Quantum ESPRESSO 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 plugin_ext_forces() + !---------------------------------------------------------------------------- + ! + ! + USE mp, ONLY : mp_bcast + USE mp_images, ONLY : intra_image_comm + USE io_global, ONLY : stdout, ionode, ionode_id + USE kinds, ONLY : DP + ! + USE plugin_flags, ONLY : use_plumed + ! + USE cell_base, ONLY : alat, at + USE ions_base, ONLY : tau, nat, amass, ityp + USE force_mod, ONLY : force,sigma + USE control_flags, ONLY : istep + USE ener, ONLY : etot + + USE plumed_module, ONLY: plumed_f_gcmd + ! + IMPLICIT NONE + ! + INTEGER:: i,j,ia + REAL(DP) :: at_plumed(3,3) + REAL(DP) :: virial(3,3) + REAL(DP) :: volume + REAL(DP), ALLOCATABLE :: tau_plumed(:,:) + REAL(DP) :: masses_plumed(nat) + ! + masses_plumed = 0.0_DP + IF(use_plumed) then + IF(ionode)THEN + at_plumed=alat*at; ! the cell, rescaled properly + allocate(tau_plumed(3,nat)) + tau_plumed=alat*tau + volume=+at_plumed(1,1)*at_plumed(2,2)*at_plumed(3,3) & + +at_plumed(1,2)*at_plumed(2,3)*at_plumed(3,1) & + +at_plumed(1,3)*at_plumed(2,1)*at_plumed(3,2) & + -at_plumed(1,1)*at_plumed(3,2)*at_plumed(2,3) & + -at_plumed(1,2)*at_plumed(3,3)*at_plumed(2,1) & + -at_plumed(1,3)*at_plumed(3,1)*at_plumed(2,2) + virial=-sigma*volume + + ! the masses in QE are stored per type, see q-e//Modules/ions_base.f90 + do ia=1,nat + masses_plumed(ia)=amass(ityp(ia)) + end do + + CALL plumed_f_gcmd("setStep"//char(0),istep) + CALL plumed_f_gcmd("setMasses"//char(0),masses_plumed) + CALL plumed_f_gcmd("setForces"//char(0),force) + CALL plumed_f_gcmd("setPositions"//char(0),tau_plumed) + CALL plumed_f_gcmd("setBox"//char(0),at_plumed) + CALL plumed_f_gcmd("setVirial"//char(0),virial) + CALL plumed_f_gcmd("setEnergy"//char(0),etot) + CALL plumed_f_gcmd("calc"//char(0),0) + + sigma=-virial/volume + + deallocate(tau_plumed) + ENDIF + CALL mp_bcast(force, ionode_id, intra_image_comm) + CALL mp_bcast(sigma, ionode_id, intra_image_comm) + ENDIF + ! + ! +END SUBROUTINE plugin_ext_forces diff --git a/patches/qespresso-7.0.diff/PW/src/plugin_ext_forces.f90.preplumed b/patches/qespresso-7.0.diff/PW/src/plugin_ext_forces.f90.preplumed new file mode 100644 index 0000000000..b184498b58 --- /dev/null +++ b/patches/qespresso-7.0.diff/PW/src/plugin_ext_forces.f90.preplumed @@ -0,0 +1,23 @@ +! +! Copyright (C) 2001-2009 Quantum ESPRESSO 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 plugin_ext_forces() + !---------------------------------------------------------------------------- + ! + ! + USE mp, ONLY : mp_bcast + USE mp_images, ONLY : intra_image_comm + USE io_global, ONLY : stdout, ionode, ionode_id + USE kinds, ONLY : DP + ! + USE plugin_flags + ! + IMPLICIT NONE + ! + ! +END SUBROUTINE plugin_ext_forces diff --git a/patches/qespresso-7.0.diff/PW/src/plugin_initialization.f90 b/patches/qespresso-7.0.diff/PW/src/plugin_initialization.f90 new file mode 100644 index 0000000000..0b41434353 --- /dev/null +++ b/patches/qespresso-7.0.diff/PW/src/plugin_initialization.f90 @@ -0,0 +1,64 @@ +! +! Copyright (C) 2010 Quantum ESPRESSO 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 plugin_initialization() + !---------------------------------------------------------------------------- + ! + USE io_global, ONLY : stdout, ionode + USE kinds, ONLY : DP + USE io_files, ONLY : tmp_dir + ! + USE plugin_flags, ONLY : use_plumed + ! + USE ions_base, ONLY : amass, ityp, nat + ! + USE dynamics_module, ONLY : dt + USE constants, ONLY : au_ps + + USE plumed_module, ONLY : plumed_f_installed, plumed_f_gcreate, plumed_f_gcmd + ! + ! + IMPLICIT NONE + ! + INTEGER :: na + INTEGER :: plumedavailable + REAL*8 :: energyUnits,lengthUnits,timeUnits + ! + IF(use_plumed) then + + CALL plumed_f_installed(plumedavailable) + + IF(plumedavailable<=0)THEN + write(stdout,*)"YOU ARE LOOKING FOR PLUMED BUT LOOKS LIKE IT IS NOT AVAILABLE: DO YOU HAVE IT IN YOUR LD_LIBRARY_PATH?" + STOP + ELSE + IF (ionode) THEN + + write(stdout,*)" CREATING PLUMED FROM THE PROGRAM" + call plumed_f_gcreate() + CALL plumed_f_gcmd("setRealPrecision"//char(0),8) + energyUnits=1312.75 ! Ry to kjoule mol + lengthUnits=0.0529177249 ! bohr to nm + timeUnits=2*au_ps ! internal time to ps + call plumed_f_gcmd("setMDEnergyUnits"//char(0),energyUnits) + call plumed_f_gcmd("setMDLengthUnits"//char(0),lengthUnits) + call plumed_f_gcmd("setMDTimeUnits"//char(0),timeUnits) + call plumed_f_gcmd("setPlumedDat"//char(0),"plumed.dat"//char(0)) + call plumed_f_gcmd("setLogFile"//char(0),"PLUMED.OUT"//char(0)) + call plumed_f_gcmd("setNatoms"//char(0),nat) + call plumed_f_gcmd("setMDEngine"//char(0),"qespresso"//char(0)); + call plumed_f_gcmd("setTimestep"//char(0),dt); + call plumed_f_gcmd("init"//char(0),0); + + + ENDIF + ENDIF + ENDIF + ! + ! +END SUBROUTINE plugin_initialization diff --git a/patches/qespresso-7.0.diff/PW/src/plugin_initialization.f90.preplumed b/patches/qespresso-7.0.diff/PW/src/plugin_initialization.f90.preplumed new file mode 100644 index 0000000000..f2c9c0097f --- /dev/null +++ b/patches/qespresso-7.0.diff/PW/src/plugin_initialization.f90.preplumed @@ -0,0 +1,21 @@ +! +! Copyright (C) 2010 Quantum ESPRESSO 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 plugin_initialization() + !---------------------------------------------------------------------------- + ! + USE io_global, ONLY : stdout, ionode + USE kinds, ONLY : DP + USE io_files, ONLY : tmp_dir + ! + USE plugin_flags + ! + IMPLICIT NONE + ! + ! +END SUBROUTINE plugin_initialization diff --git a/patches/qespresso-7.0.diff/PW/src/run_pwscf.f90 b/patches/qespresso-7.0.diff/PW/src/run_pwscf.f90 new file mode 100644 index 0000000000..e2883c098e --- /dev/null +++ b/patches/qespresso-7.0.diff/PW/src/run_pwscf.f90 @@ -0,0 +1,532 @@ +! +! Copyright (C) 2013-2020 Quantum ESPRESSO 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 run_pwscf( exit_status ) + !---------------------------------------------------------------------------- + !! Author: Paolo Giannozzi + !! License: GNU + !! Summary: Run an instance of the Plane Wave Self-Consistent Field code + ! + !! Run an instance of the Plane Wave Self-Consistent Field code + !! MPI initialization and input data reading is performed in the + !! calling code - returns in exit_status the exit code for pw.x, + !! returned in the shell. Values are: + !! * 0: completed successfully + !! * 1: an error has occurred (value returned by the errore() routine) + !! * 2-127: convergence error + !! * 2: scf convergence error + !! * 3: ion convergence error + !! * 128-255: code exited due to specific trigger + !! * 255: exit due to user request, or signal trapped, + !! or time > max_seconds + !! (note: in the future, check_stop_now could also return a value + !! to specify the reason of exiting, and the value could be used + !! to return a different value for different reasons) + ! + !! @Note + !! 10/01/17 Samuel Ponce: Add Ford documentation + !! @endnote + !! + ! + USE kinds, ONLY : DP + USE mp, ONLY : mp_bcast, mp_sum + USE io_global, ONLY : stdout, ionode, ionode_id + USE parameters, ONLY : ntypx, npk + USE upf_params, ONLY : lmaxx + USE cell_base, ONLY : fix_volume, fix_area + USE control_flags, ONLY : conv_elec, gamma_only, ethr, lscf, treinit_gvecs + USE control_flags, ONLY : conv_ions, istep, nstep, restart, lmd, lbfgs,& + lensemb, lforce=>tprnfor, tstress + USE control_flags, ONLY : io_level + USE cellmd, ONLY : lmovecell + USE command_line_options, ONLY : command_line + USE force_mod, ONLY : sigma, force + USE check_stop, ONLY : check_stop_init, check_stop_now + USE mp_images, ONLY : intra_image_comm + USE extrapolation, ONLY : update_file, update_pot + USE scf, ONLY : rho + USE lsda_mod, ONLY : nspin + USE fft_base, ONLY : dfftp + ! + USE plugin_flags, ONLY : use_plumed + ! + USE qmmm, ONLY : qmmm_initialization, qmmm_shutdown, & + qmmm_update_positions, qmmm_update_forces + USE qexsd_module, ONLY : qexsd_set_status + USE xc_lib, ONLY : xclib_dft_is, stop_exx + USE beef, ONLY : beef_energies + USE ldaU, ONLY : lda_plus_u + USE add_dmft_occ, ONLY : dmft + ! + USE device_fbuff_m, ONLY : dev_buf + USE plumed_module, ONLY : plumed_f_gfinalize + ! + IMPLICIT NONE + ! + INTEGER, INTENT(OUT) :: exit_status + !! Gives the exit status at the end + ! + LOGICAL, EXTERNAL :: matches + ! checks if first string is contained in the second + ! + ! ... local variables + ! + INTEGER :: idone + ! counter of electronic + ionic steps done in this run + INTEGER :: ions_status + ! ions_status = 3 not yet converged + ! ions_status = 2 converged, restart with nonzero magnetization + ! ions_status = 1 converged, final step with current cell needed + ! ions_status = 0 converged, exiting + ! + LOGICAL :: optimizer_failed = .FALSE. + ! + INTEGER :: ierr + ! collect error codes + ! + ions_status = 3 + exit_status = 0 + IF ( ionode ) WRITE( UNIT = stdout, FMT = 9010 ) ntypx, npk, lmaxx + ! + IF (ionode) CALL plugin_arguments() + CALL plugin_arguments_bcast( ionode_id, intra_image_comm ) + ! + ! ... needs to come before iosys() so some input flags can be + ! overridden without needing to write PWscf specific code. + ! + CALL qmmm_initialization() + ! + ! ... convert to internal variables + ! + CALL iosys() + ! + ! ... If executable names is "dist.x", compute atomic distances, angles, + ! ... nearest neighbors, write them to file "dist.out", exit + ! + IF ( matches('dist.x',command_line) ) THEN + IF (ionode) CALL run_dist( exit_status ) + RETURN + ENDIF + ! + IF ( gamma_only ) WRITE( UNIT = stdout, & + & FMT = '(/,5X,"gamma-point specific algorithms are used")' ) + ! + ! call to void routine for user defined / plugin patches initializations + ! + CALL plugin_initialization() + ! + CALL check_stop_init() + ! + CALL setup() + ! + CALL qmmm_update_positions() + ! + ! ... dry run: code will stop here if called with exit file present + ! ... useful for a quick and automated way to check input data + ! + IF ( nstep == 0 .OR. check_stop_now() ) THEN + CALL pre_init() + CALL data_structure( gamma_only ) + CALL summary() + CALL memory_report() + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config-init' ) + RETURN + ENDIF + ! + CALL init_run() + ! + IF ( check_stop_now() ) THEN + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config' ) + RETURN + ENDIF + ! + main_loop: DO idone = 1, nstep + ! + ! ... electronic self-consistency or band structure calculation + ! + IF ( .NOT. lscf) THEN + CALL non_scf() + ELSE + CALL electrons() + END IF + ! + ! ... code stopped by user or not converged + ! + IF ( check_stop_now() .OR. .NOT. conv_elec ) THEN + IF ( check_stop_now() ) exit_status = 255 + IF ( .NOT. conv_elec) THEN + IF (dmft) exit_status = 131 + ELSE + exit_status = 2 + ENDIF + CALL qexsd_set_status(exit_status) + CALL punch( 'config' ) + RETURN + ENDIF + ! + ! ... file in CASINO format written here if required + ! + IF ( lmd ) THEN + CALL pw2casino( istep ) + ELSE + CALL pw2casino( 0 ) + END IF + ! + ! ... ionic section starts here + ! + CALL start_clock( 'ions' ); !write(*,*)' start ions' ; FLUSH(6) + conv_ions = .TRUE. + ! + ! ... force calculation + ! + IF ( lforce ) CALL forces() + ! + ! ... stress calculation + ! + IF ( tstress ) CALL stress( sigma ) + ! + ! ... this is the right place for plugin forces: + ! + CALL plugin_ext_forces() + ! + IF ( lmd .OR. lbfgs ) THEN + ! + ! ... add information on this ionic step to xml file + ! + CALL add_qexsd_step( idone ) + ! + IF (fix_volume) CALL impose_deviatoric_stress( sigma ) + IF (fix_area) CALL impose_deviatoric_stress_2d( sigma ) + ! + ! ... save data needed for potential and wavefunction extrapolation + ! + CALL update_file() + ! + ! ... ionic step (for molecular dynamics or optimization) + ! + CALL move_ions ( idone, ions_status, optimizer_failed ) + conv_ions = ( ions_status == 0 ) .OR. & + ( ions_status == 1 .AND. treinit_gvecs ) + ! + IF ( xclib_dft_is('hybrid') ) CALL stop_exx() + ! + ! ... save restart information for the new configuration + ! + IF ( idone <= nstep .AND. .NOT. conv_ions ) THEN + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config-only' ) + END IF + ! + END IF + ! + CALL stop_clock( 'ions' ); !write(*,*)' stop ions' ; FLUSH(6) + ! + ! ... send out forces to MM code in QM/MM run + ! + CALL qmmm_update_forces( force, rho%of_r, nspin, dfftp ) + ! + ! ... exit condition (ionic convergence) is checked here + ! + IF ( conv_ions .OR. optimizer_failed ) EXIT main_loop + ! + ! ... receive new positions from MM code in QM/MM run + ! + CALL qmmm_update_positions() + ! + ! ... terms of the hamiltonian depending upon nuclear positions + ! ... are reinitialized here + ! + IF ( lmd .OR. lbfgs ) THEN + ! + IF ( ions_status == 1 ) THEN + ! + ! ... final scf calculation with G-vectors for final cell + ! + lbfgs=.FALSE.; lmd=.FALSE. + WRITE( UNIT = stdout, FMT=9020 ) + ! + CALL reset_gvectors( ) + ! + ! ... read atomic occupations for DFT+U(+V) + ! + IF ( lda_plus_u ) CALL read_ns() + ! + ELSE IF ( ions_status == 2 ) THEN + ! + ! ... check whether nonzero magnetization is real + ! + CALL reset_magn() + ! + ELSE + ! + IF ( treinit_gvecs ) THEN + ! + ! ... prepare for next step with freshly computed G vectors + ! + IF ( lmovecell) CALL scale_h() + CALL reset_gvectors ( ) + ! + ELSE + ! + ! ... update the wavefunctions, charge density, potential + ! ... update_pot initializes structure factor array as well + ! + CALL update_pot() + ! + ! ... re-initialize atomic position-dependent quantities + ! + CALL hinit1() + ! + END IF + ! + END IF + ! + ENDIF + ! ... Reset convergence threshold of iterative diagonalization for + ! ... the first scf iteration of each ionic step (after the first) + ! + ethr = 1.0D-6 + ! + CALL dev_buf%reinit( ierr ) + IF ( ierr .ne. 0 ) CALL infomsg( 'run_pwscf', 'Cannot reset GPU buffers! Some buffers still locked.' ) + ! + ENDDO main_loop + ! + ! Set correct exit_status + ! + IF ( .NOT. conv_ions .OR. optimizer_failed ) THEN + exit_status = 3 + ELSE + ! All good + exit_status = 0 + END IF + ! + ! ... save final data file + ! + CALL qexsd_set_status( exit_status ) + IF ( lensemb ) CALL beef_energies( ) + IF ( io_level > -2 ) CALL punch( 'all' ) + ! + ! finalize plumed + ! + IF(use_plumed) then + CALL plumed_f_gfinalize() + ENDIF + ! + CALL qmmm_shutdown() + ! + RETURN + ! +9010 FORMAT( /,5X,'Current dimensions of program PWSCF are:', & + & /,5X,'Max number of different atomic species (ntypx) = ',I2,& + & /,5X,'Max number of k-points (npk) = ',I6, & + & /,5X,'Max angular momentum in pseudopotentials (lmaxx) = ',i2) +9020 FORMAT( /,5X,'Final scf calculation at the relaxed structure.', & + & /,5X,'The G-vectors are recalculated for the final unit cell', & + & /,5X,'Results may differ from those at the preceding step.' ) + ! +END SUBROUTINE run_pwscf +! +! +!------------------------------------------------------------- +SUBROUTINE reset_gvectors( ) +!------------------------------------------------------------- + ! + !! Prepare a new scf calculation with newly recomputed grids, + !! restarting from scratch, not from available data of previous + !! steps (dimensions and file lengths will be different in general) + !! Useful as a check of variable-cell optimization: + !! once convergence is achieved, compare the final energy with the + !! energy computed with G-vectors and plane waves for the final cell + ! + USE io_global, ONLY : stdout + USE basis, ONLY : starting_wfc, starting_pot + USE fft_base, ONLY : dfftp + USE fft_base, ONLY : dffts + USE xc_lib, ONLY : xclib_dft_is + ! + IMPLICIT NONE + ! + ! ... get magnetic moments from previous run before charge is deleted + ! + CALL reset_starting_magnetization() + ! + ! ... clean everything (FIXME: clean only what has to be cleaned) + ! + CALL clean_pw( .FALSE. ) + CALL close_files(.TRUE.) + ! + IF (TRIM(starting_wfc) == 'file') starting_wfc = 'atomic+random' + starting_pot='atomic' + ! + ! ... re-set FFT grids and re-compute needed stuff (FIXME: which?) + ! + dfftp%nr1=0; dfftp%nr2=0; dfftp%nr3=0 + dffts%nr1=0; dffts%nr2=0; dffts%nr3=0 + ! + CALL init_run() + ! + ! ... re-set and re-initialize EXX-related stuff + ! + IF ( xclib_dft_is('hybrid') ) CALL reset_exx( ) + ! +END SUBROUTINE reset_gvectors +! +! +!------------------------------------------------------------- +SUBROUTINE reset_exx( ) +!------------------------------------------------------------- + USE fft_types, ONLY : fft_type_deallocate + USE exx_base, ONLY : exx_grid_init, exx_mp_init, exx_div_check, & + coulomb_fac, coulomb_done + USE exx, ONLY : dfftt, exx_fft_create, deallocate_exx + USE exx_band, ONLY : igk_exx + ! + IMPLICIT NONE + ! + ! ... re-set EXX-related stuff... + ! + IF (ALLOCATED(coulomb_fac) ) DEALLOCATE( coulomb_fac, coulomb_done ) + CALL deallocate_exx( ) + IF (ALLOCATED(igk_exx)) DEALLOCATE(igk_exx) + dfftt%nr1=0; dfftt%nr2=0; dfftt%nr3=0 + CALL fft_type_deallocate( dfftt ) ! FIXME: is this needed? + ! + ! ... re-compute needed EXX-related stuff + ! + CALL exx_grid_init( REINIT = .TRUE. ) + CALL exx_mp_init() + CALL exx_fft_create() + CALL exx_div_check() + ! +END SUBROUTINE reset_exx +! +! +!---------------------------------------------------------------- +SUBROUTINE reset_magn() + !---------------------------------------------------------------- + !! LSDA optimization: a final configuration with zero + !! absolute magnetization has been found and we check + !! if it is really the minimum energy structure by + !! performing a new scf iteration without any "electronic" history. + ! + USE io_global, ONLY : stdout + USE dfunct, ONLY : newd + ! + IMPLICIT NONE + ! + WRITE( UNIT = stdout, FMT = 9010 ) + WRITE( UNIT = stdout, FMT = 9020 ) + ! + ! ... re-initialize the potential (no need to re-initialize wavefunctions) + ! + CALL potinit() + CALL newd() + ! +9010 FORMAT( /5X,'lsda relaxation : a final configuration with zero', & + & /5X,' absolute magnetization has been found' ) +9020 FORMAT( /5X,'the program is checking if it is really ', & + & 'the minimum energy structure', & + & /5X,'by performing a new scf iteration ', & + & 'without any "electronic" history' ) + ! +END SUBROUTINE reset_magn +! +! +!------------------------------------------------------------------- +SUBROUTINE reset_starting_magnetization() + !------------------------------------------------------------------- + !! On input, the scf charge density is needed. + !! On output, new values for starting_magnetization, angle1, angle2 + !! estimated from atomic magnetic moments - to be used in last step. + ! + USE kinds, ONLY : DP + USE constants, ONLY : pi + USE ions_base, ONLY : nsp, ityp, nat + USE lsda_mod, ONLY : nspin, starting_magnetization + USE scf, ONLY : rho + USE noncollin_module, ONLY : noncolin, angle1, angle2, domag + ! + IMPLICIT NONE + ! + ! ... local variables + ! + INTEGER :: i, nt, iat + ! loop counter on species + ! number of atoms per species + ! loop counter on atoms + REAL(DP) :: norm_tot, norm_xy + ! modulus of atomic magnetization + ! xy-projection of atomic magnetization + REAL(DP) :: theta, phi + ! angle between magnetization and z-axis + ! angle between xy-magnetization and x-axis + REAL(DP), ALLOCATABLE :: r_loc(:) + ! auxiliary array for density + REAL(DP), ALLOCATABLE :: m_loc(:,:) + ! auxiliary array for magnetization + ! + IF ( (noncolin .AND. domag) .OR. nspin==2) THEN + ALLOCATE( r_loc(nat), m_loc(nspin-1,nat) ) + CALL get_locals( r_loc,m_loc, rho%of_r ) + ELSE + RETURN + ENDIF + ! + DO i = 1, nsp + ! + starting_magnetization(i) = 0.0_DP + angle1(i) = 0.0_DP + angle2(i) = 0.0_DP + nt = 0 + ! + DO iat = 1, nat + IF (ityp(iat) == i) THEN + nt = nt + 1 + IF (noncolin) THEN + norm_tot = SQRT(m_loc(1,iat)**2+m_loc(2,iat)**2+m_loc(3,iat)**2) + norm_xy = SQRT(m_loc(1,iat)**2+m_loc(2,iat)**2) + IF (norm_tot > 1.d-10) THEN + theta = ACOS(m_loc(3,iat)/norm_tot) + IF (norm_xy > 1.d-10) THEN + phi = ACOS(m_loc(1,iat)/norm_xy) + IF (m_loc(2,iat) < 0.d0) phi = - phi + ELSE + phi = 2.d0*pi + ENDIF + ELSE + theta = 2.d0*pi + phi = 2.d0*pi + ENDIF + angle1(i) = angle1(i) + theta + angle2(i) = angle2(i) + phi + starting_magnetization(i) = starting_magnetization(i) + & + norm_tot/r_loc(iat) + ELSE + starting_magnetization(i) = starting_magnetization(i) + & + m_loc(1,iat)/r_loc(iat) + ENDIF + ENDIF + ENDDO + ! + IF ( nt > 0 ) THEN + starting_magnetization(i) = starting_magnetization(i) / DBLE(nt) + angle1(i) = angle1(i) / DBLE(nt) + angle2(i) = angle2(i) / DBLE(nt) + ENDIF + ! + ENDDO + ! + DEALLOCATE( r_loc, m_loc ) + ! +END SUBROUTINE reset_starting_magnetization diff --git a/patches/qespresso-7.0.diff/PW/src/run_pwscf.f90.preplumed b/patches/qespresso-7.0.diff/PW/src/run_pwscf.f90.preplumed new file mode 100644 index 0000000000..b31516afac --- /dev/null +++ b/patches/qespresso-7.0.diff/PW/src/run_pwscf.f90.preplumed @@ -0,0 +1,518 @@ +! +! Copyright (C) 2013-2020 Quantum ESPRESSO 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 run_pwscf( exit_status ) + !---------------------------------------------------------------------------- + !! Author: Paolo Giannozzi + !! License: GNU + !! Summary: Run an instance of the Plane Wave Self-Consistent Field code + ! + !! Run an instance of the Plane Wave Self-Consistent Field code + !! MPI initialization and input data reading is performed in the + !! calling code - returns in exit_status the exit code for pw.x, + !! returned in the shell. Values are: + !! * 0: completed successfully + !! * 1: an error has occurred (value returned by the errore() routine) + !! * 2-127: convergence error + !! * 2: scf convergence error + !! * 3: ion convergence error + !! * 128-255: code exited due to specific trigger + !! * 255: exit due to user request, or signal trapped, + !! or time > max_seconds + !! (note: in the future, check_stop_now could also return a value + !! to specify the reason of exiting, and the value could be used + !! to return a different value for different reasons) + ! + !! @Note + !! 10/01/17 Samuel Ponce: Add Ford documentation + !! @endnote + !! + ! + USE kinds, ONLY : DP + USE mp, ONLY : mp_bcast, mp_sum + USE io_global, ONLY : stdout, ionode, ionode_id + USE parameters, ONLY : ntypx, npk + USE upf_params, ONLY : lmaxx + USE cell_base, ONLY : fix_volume, fix_area + USE control_flags, ONLY : conv_elec, gamma_only, ethr, lscf, treinit_gvecs + USE control_flags, ONLY : conv_ions, istep, nstep, restart, lmd, lbfgs,& + lensemb, lforce=>tprnfor, tstress + USE control_flags, ONLY : io_level + USE cellmd, ONLY : lmovecell + USE command_line_options, ONLY : command_line + USE force_mod, ONLY : sigma, force + USE check_stop, ONLY : check_stop_init, check_stop_now + USE mp_images, ONLY : intra_image_comm + USE extrapolation, ONLY : update_file, update_pot + USE scf, ONLY : rho + USE lsda_mod, ONLY : nspin + USE fft_base, ONLY : dfftp + USE qmmm, ONLY : qmmm_initialization, qmmm_shutdown, & + qmmm_update_positions, qmmm_update_forces + USE qexsd_module, ONLY : qexsd_set_status + USE xc_lib, ONLY : xclib_dft_is, stop_exx + USE beef, ONLY : beef_energies + USE ldaU, ONLY : lda_plus_u + USE add_dmft_occ, ONLY : dmft + ! + USE device_fbuff_m, ONLY : dev_buf + ! + IMPLICIT NONE + ! + INTEGER, INTENT(OUT) :: exit_status + !! Gives the exit status at the end + ! + LOGICAL, EXTERNAL :: matches + ! checks if first string is contained in the second + ! + ! ... local variables + ! + INTEGER :: idone + ! counter of electronic + ionic steps done in this run + INTEGER :: ions_status + ! ions_status = 3 not yet converged + ! ions_status = 2 converged, restart with nonzero magnetization + ! ions_status = 1 converged, final step with current cell needed + ! ions_status = 0 converged, exiting + ! + LOGICAL :: optimizer_failed = .FALSE. + ! + INTEGER :: ierr + ! collect error codes + ! + ions_status = 3 + exit_status = 0 + IF ( ionode ) WRITE( UNIT = stdout, FMT = 9010 ) ntypx, npk, lmaxx + ! + IF (ionode) CALL plugin_arguments() + CALL plugin_arguments_bcast( ionode_id, intra_image_comm ) + ! + ! ... needs to come before iosys() so some input flags can be + ! overridden without needing to write PWscf specific code. + ! + CALL qmmm_initialization() + ! + ! ... convert to internal variables + ! + CALL iosys() + ! + ! ... If executable names is "dist.x", compute atomic distances, angles, + ! ... nearest neighbors, write them to file "dist.out", exit + ! + IF ( matches('dist.x',command_line) ) THEN + IF (ionode) CALL run_dist( exit_status ) + RETURN + ENDIF + ! + IF ( gamma_only ) WRITE( UNIT = stdout, & + & FMT = '(/,5X,"gamma-point specific algorithms are used")' ) + ! + ! call to void routine for user defined / plugin patches initializations + ! + CALL plugin_initialization() + ! + CALL check_stop_init() + ! + CALL setup() + ! + CALL qmmm_update_positions() + ! + ! ... dry run: code will stop here if called with exit file present + ! ... useful for a quick and automated way to check input data + ! + IF ( nstep == 0 .OR. check_stop_now() ) THEN + CALL pre_init() + CALL data_structure( gamma_only ) + CALL summary() + CALL memory_report() + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config-init' ) + RETURN + ENDIF + ! + CALL init_run() + ! + IF ( check_stop_now() ) THEN + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config' ) + RETURN + ENDIF + ! + main_loop: DO idone = 1, nstep + ! + ! ... electronic self-consistency or band structure calculation + ! + IF ( .NOT. lscf) THEN + CALL non_scf() + ELSE + CALL electrons() + END IF + ! + ! ... code stopped by user or not converged + ! + IF ( check_stop_now() .OR. .NOT. conv_elec ) THEN + IF ( check_stop_now() ) exit_status = 255 + IF ( .NOT. conv_elec) THEN + IF (dmft) exit_status = 131 + ELSE + exit_status = 2 + ENDIF + CALL qexsd_set_status(exit_status) + CALL punch( 'config' ) + RETURN + ENDIF + ! + ! ... file in CASINO format written here if required + ! + IF ( lmd ) THEN + CALL pw2casino( istep ) + ELSE + CALL pw2casino( 0 ) + END IF + ! + ! ... ionic section starts here + ! + CALL start_clock( 'ions' ); !write(*,*)' start ions' ; FLUSH(6) + conv_ions = .TRUE. + ! + ! ... force calculation + ! + IF ( lforce ) CALL forces() + ! + ! ... stress calculation + ! + IF ( tstress ) CALL stress( sigma ) + ! + IF ( lmd .OR. lbfgs ) THEN + ! + ! ... add information on this ionic step to xml file + ! + CALL add_qexsd_step( idone ) + ! + IF (fix_volume) CALL impose_deviatoric_stress( sigma ) + IF (fix_area) CALL impose_deviatoric_stress_2d( sigma ) + ! + ! ... save data needed for potential and wavefunction extrapolation + ! + CALL update_file() + ! + ! ... ionic step (for molecular dynamics or optimization) + ! + CALL move_ions ( idone, ions_status, optimizer_failed ) + conv_ions = ( ions_status == 0 ) .OR. & + ( ions_status == 1 .AND. treinit_gvecs ) + ! + IF ( xclib_dft_is('hybrid') ) CALL stop_exx() + ! + ! ... save restart information for the new configuration + ! + IF ( idone <= nstep .AND. .NOT. conv_ions ) THEN + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config-only' ) + END IF + ! + END IF + ! + CALL stop_clock( 'ions' ); !write(*,*)' stop ions' ; FLUSH(6) + ! + ! ... send out forces to MM code in QM/MM run + ! + CALL qmmm_update_forces( force, rho%of_r, nspin, dfftp ) + ! + ! ... exit condition (ionic convergence) is checked here + ! + IF ( conv_ions .OR. optimizer_failed ) EXIT main_loop + ! + ! ... receive new positions from MM code in QM/MM run + ! + CALL qmmm_update_positions() + ! + ! ... terms of the hamiltonian depending upon nuclear positions + ! ... are reinitialized here + ! + IF ( lmd .OR. lbfgs ) THEN + ! + IF ( ions_status == 1 ) THEN + ! + ! ... final scf calculation with G-vectors for final cell + ! + lbfgs=.FALSE.; lmd=.FALSE. + WRITE( UNIT = stdout, FMT=9020 ) + ! + CALL reset_gvectors( ) + ! + ! ... read atomic occupations for DFT+U(+V) + ! + IF ( lda_plus_u ) CALL read_ns() + ! + ELSE IF ( ions_status == 2 ) THEN + ! + ! ... check whether nonzero magnetization is real + ! + CALL reset_magn() + ! + ELSE + ! + IF ( treinit_gvecs ) THEN + ! + ! ... prepare for next step with freshly computed G vectors + ! + IF ( lmovecell) CALL scale_h() + CALL reset_gvectors ( ) + ! + ELSE + ! + ! ... update the wavefunctions, charge density, potential + ! ... update_pot initializes structure factor array as well + ! + CALL update_pot() + ! + ! ... re-initialize atomic position-dependent quantities + ! + CALL hinit1() + ! + END IF + ! + END IF + ! + ENDIF + ! ... Reset convergence threshold of iterative diagonalization for + ! ... the first scf iteration of each ionic step (after the first) + ! + ethr = 1.0D-6 + ! + CALL dev_buf%reinit( ierr ) + IF ( ierr .ne. 0 ) CALL infomsg( 'run_pwscf', 'Cannot reset GPU buffers! Some buffers still locked.' ) + ! + ENDDO main_loop + ! + ! Set correct exit_status + ! + IF ( .NOT. conv_ions .OR. optimizer_failed ) THEN + exit_status = 3 + ELSE + ! All good + exit_status = 0 + END IF + ! + ! ... save final data file + ! + CALL qexsd_set_status( exit_status ) + IF ( lensemb ) CALL beef_energies( ) + IF ( io_level > -2 ) CALL punch( 'all' ) + ! + CALL qmmm_shutdown() + ! + RETURN + ! +9010 FORMAT( /,5X,'Current dimensions of program PWSCF are:', & + & /,5X,'Max number of different atomic species (ntypx) = ',I2,& + & /,5X,'Max number of k-points (npk) = ',I6, & + & /,5X,'Max angular momentum in pseudopotentials (lmaxx) = ',i2) +9020 FORMAT( /,5X,'Final scf calculation at the relaxed structure.', & + & /,5X,'The G-vectors are recalculated for the final unit cell', & + & /,5X,'Results may differ from those at the preceding step.' ) + ! +END SUBROUTINE run_pwscf +! +! +!------------------------------------------------------------- +SUBROUTINE reset_gvectors( ) +!------------------------------------------------------------- + ! + !! Prepare a new scf calculation with newly recomputed grids, + !! restarting from scratch, not from available data of previous + !! steps (dimensions and file lengths will be different in general) + !! Useful as a check of variable-cell optimization: + !! once convergence is achieved, compare the final energy with the + !! energy computed with G-vectors and plane waves for the final cell + ! + USE io_global, ONLY : stdout + USE basis, ONLY : starting_wfc, starting_pot + USE fft_base, ONLY : dfftp + USE fft_base, ONLY : dffts + USE xc_lib, ONLY : xclib_dft_is + ! + IMPLICIT NONE + ! + ! ... get magnetic moments from previous run before charge is deleted + ! + CALL reset_starting_magnetization() + ! + ! ... clean everything (FIXME: clean only what has to be cleaned) + ! + CALL clean_pw( .FALSE. ) + CALL close_files(.TRUE.) + ! + IF (TRIM(starting_wfc) == 'file') starting_wfc = 'atomic+random' + starting_pot='atomic' + ! + ! ... re-set FFT grids and re-compute needed stuff (FIXME: which?) + ! + dfftp%nr1=0; dfftp%nr2=0; dfftp%nr3=0 + dffts%nr1=0; dffts%nr2=0; dffts%nr3=0 + ! + CALL init_run() + ! + ! ... re-set and re-initialize EXX-related stuff + ! + IF ( xclib_dft_is('hybrid') ) CALL reset_exx( ) + ! +END SUBROUTINE reset_gvectors +! +! +!------------------------------------------------------------- +SUBROUTINE reset_exx( ) +!------------------------------------------------------------- + USE fft_types, ONLY : fft_type_deallocate + USE exx_base, ONLY : exx_grid_init, exx_mp_init, exx_div_check, & + coulomb_fac, coulomb_done + USE exx, ONLY : dfftt, exx_fft_create, deallocate_exx + USE exx_band, ONLY : igk_exx + ! + IMPLICIT NONE + ! + ! ... re-set EXX-related stuff... + ! + IF (ALLOCATED(coulomb_fac) ) DEALLOCATE( coulomb_fac, coulomb_done ) + CALL deallocate_exx( ) + IF (ALLOCATED(igk_exx)) DEALLOCATE(igk_exx) + dfftt%nr1=0; dfftt%nr2=0; dfftt%nr3=0 + CALL fft_type_deallocate( dfftt ) ! FIXME: is this needed? + ! + ! ... re-compute needed EXX-related stuff + ! + CALL exx_grid_init( REINIT = .TRUE. ) + CALL exx_mp_init() + CALL exx_fft_create() + CALL exx_div_check() + ! +END SUBROUTINE reset_exx +! +! +!---------------------------------------------------------------- +SUBROUTINE reset_magn() + !---------------------------------------------------------------- + !! LSDA optimization: a final configuration with zero + !! absolute magnetization has been found and we check + !! if it is really the minimum energy structure by + !! performing a new scf iteration without any "electronic" history. + ! + USE io_global, ONLY : stdout + USE dfunct, ONLY : newd + ! + IMPLICIT NONE + ! + WRITE( UNIT = stdout, FMT = 9010 ) + WRITE( UNIT = stdout, FMT = 9020 ) + ! + ! ... re-initialize the potential (no need to re-initialize wavefunctions) + ! + CALL potinit() + CALL newd() + ! +9010 FORMAT( /5X,'lsda relaxation : a final configuration with zero', & + & /5X,' absolute magnetization has been found' ) +9020 FORMAT( /5X,'the program is checking if it is really ', & + & 'the minimum energy structure', & + & /5X,'by performing a new scf iteration ', & + & 'without any "electronic" history' ) + ! +END SUBROUTINE reset_magn +! +! +!------------------------------------------------------------------- +SUBROUTINE reset_starting_magnetization() + !------------------------------------------------------------------- + !! On input, the scf charge density is needed. + !! On output, new values for starting_magnetization, angle1, angle2 + !! estimated from atomic magnetic moments - to be used in last step. + ! + USE kinds, ONLY : DP + USE constants, ONLY : pi + USE ions_base, ONLY : nsp, ityp, nat + USE lsda_mod, ONLY : nspin, starting_magnetization + USE scf, ONLY : rho + USE noncollin_module, ONLY : noncolin, angle1, angle2, domag + ! + IMPLICIT NONE + ! + ! ... local variables + ! + INTEGER :: i, nt, iat + ! loop counter on species + ! number of atoms per species + ! loop counter on atoms + REAL(DP) :: norm_tot, norm_xy + ! modulus of atomic magnetization + ! xy-projection of atomic magnetization + REAL(DP) :: theta, phi + ! angle between magnetization and z-axis + ! angle between xy-magnetization and x-axis + REAL(DP), ALLOCATABLE :: r_loc(:) + ! auxiliary array for density + REAL(DP), ALLOCATABLE :: m_loc(:,:) + ! auxiliary array for magnetization + ! + IF ( (noncolin .AND. domag) .OR. nspin==2) THEN + ALLOCATE( r_loc(nat), m_loc(nspin-1,nat) ) + CALL get_locals( r_loc,m_loc, rho%of_r ) + ELSE + RETURN + ENDIF + ! + DO i = 1, nsp + ! + starting_magnetization(i) = 0.0_DP + angle1(i) = 0.0_DP + angle2(i) = 0.0_DP + nt = 0 + ! + DO iat = 1, nat + IF (ityp(iat) == i) THEN + nt = nt + 1 + IF (noncolin) THEN + norm_tot = SQRT(m_loc(1,iat)**2+m_loc(2,iat)**2+m_loc(3,iat)**2) + norm_xy = SQRT(m_loc(1,iat)**2+m_loc(2,iat)**2) + IF (norm_tot > 1.d-10) THEN + theta = ACOS(m_loc(3,iat)/norm_tot) + IF (norm_xy > 1.d-10) THEN + phi = ACOS(m_loc(1,iat)/norm_xy) + IF (m_loc(2,iat) < 0.d0) phi = - phi + ELSE + phi = 2.d0*pi + ENDIF + ELSE + theta = 2.d0*pi + phi = 2.d0*pi + ENDIF + angle1(i) = angle1(i) + theta + angle2(i) = angle2(i) + phi + starting_magnetization(i) = starting_magnetization(i) + & + norm_tot/r_loc(iat) + ELSE + starting_magnetization(i) = starting_magnetization(i) + & + m_loc(1,iat)/r_loc(iat) + ENDIF + ENDIF + ENDDO + ! + IF ( nt > 0 ) THEN + starting_magnetization(i) = starting_magnetization(i) / DBLE(nt) + angle1(i) = angle1(i) / DBLE(nt) + angle2(i) = angle2(i) / DBLE(nt) + ENDIF + ! + ENDDO + ! + DEALLOCATE( r_loc, m_loc ) + ! +END SUBROUTINE reset_starting_magnetization diff --git a/patches/qespresso-7.2.diff/Modules/Makefile b/patches/qespresso-7.2.diff/Modules/Makefile new file mode 100644 index 0000000000..cf1e729f18 --- /dev/null +++ b/patches/qespresso-7.2.diff/Modules/Makefile @@ -0,0 +1,249 @@ +#/a Makefile for Modules + +include ../make.inc + +# location of needed modules +MODFLAGS=$(BASEMOD_FLAGS) \ + $(MOD_FLAG)../ELPA/src + +# list of modules + +MODULES = \ +additional_kpoints.o \ +autopilot.o \ +basic_algebra_routines.o \ +becmod.o \ +bfgs_module.o \ +bspline.o \ +bz_form.o \ +cell_base.o \ +check_stop.o \ +command_line_options.o \ +compute_dipole.o \ +constants.o \ +constraints_module.o \ +control_flags.o \ +coulomb_vcut.o \ +dist.o \ +electrons_base.o \ +environ_base_module.o \ +environment.o \ +extffield.o \ +fd_gradient.o \ +fft_base.o \ +fft_rho.o \ +fft_wave.o \ +fsockets.o \ +funct.o \ +generate_function.o \ +gradutils.o \ +gvecw.o \ +input_parameters.o \ +invmat.o \ +io_files.o \ +io_global.o \ +ions_base.o \ +kind.o \ +lmdif.o \ +mdiis.o \ +mm_dispersion.o \ +mp_bands.o \ +mp_exx.o \ +mp_global.o \ +mp_images.o \ +mp_pools.o \ +mp_wave.o \ +mp_world.o \ +noncol.o \ +open_close_input_file.o \ +parameters.o \ +parser.o \ +plugin_flags.o \ +plugin_arguments.o \ +plugin_variables.o \ +pw_dot.o \ +qmmm.o \ +random_numbers.o \ +read_cards.o \ +read_input.o \ +read_namelists.o \ +read_pseudo.o \ +recvec.o \ +recvec_subs.o \ +run_info.o \ +space_group.o \ +set_para_diag.o \ +set_signal.o \ +set_vdw_corr.o \ +setqf.o \ +timestep.o\ +tsvdw.o\ +mbdlib.o\ +version.o \ +wannier_gw.o\ +wannier_new.o \ +wavefunctions.o \ +ws_base.o \ +xc_vdW_DF.o \ +xc_rVV10.o \ +io_base.o \ +qes_types_module.o \ +qes_libs_module.o \ +qes_write_module.o \ +qes_read_module.o \ +qes_reset_module.o \ +qes_init_module.o \ +qes_bcast_module.o \ +qexsd.o \ +qexsd_copy.o \ +qexsd_init.o \ +qexsd_input.o \ +hdf5_qe.o\ +qeh5_module.o\ +fox_init_module.o \ +xsf.o \ +wyckoff.o \ +wypos.o \ +zvscal.o \ +wave_gauge.o \ +plumed.o + +# list of RISM's modules + +RISMLIB = \ +allocate_fft_3drism.o \ +chempot.o \ +chempot_lauerism.o \ +closure.o \ +corrdipole_laue.o \ +correctat0_vv.o \ +corrgxy0_laue.o \ +cryst_to_car_2d.o \ +data_structure_3drism.o \ +do_1drism.o \ +do_3drism.o \ +do_lauerism.o \ +eqn_1drism.o \ +eqn_3drism.o \ +eqn_lauedipole.o \ +eqn_lauegxy0.o \ +eqn_lauelong.o \ +eqn_lauerism.o \ +eqn_laueshort.o \ +eqn_lauevoid.o \ +err_rism.o \ +guess_3drism.o \ +init_1drism.o \ +init_3drism.o \ +input_1drism.o \ +input_3drism.o \ +io_rism_xml.o \ +lauefft.o \ +lauefft_subs.o \ +lj_forcefield.o \ +lj_solute.o \ +molecorr_vv.o \ +molebridge_vv.o \ +molecule_const.o \ +molecule_types.o \ +mp_rism.o \ +mp_swap_ax_rism.o \ +normalize_lauerism.o \ +plot_rism.o \ +potential_3drism.o \ +potential_esm.o \ +potential_vv.o \ +print_chempot_3drism.o \ +print_chempot_lauerism.o \ +print_chempot_vv.o \ +print_corr_vv.o \ +print_solvavg.o \ +radfft.o \ +read_mol.o \ +read_solv.o \ +recvec_3drism.o \ +rism.o \ +rism1d_facade.o \ +rism3d_facade.o \ +rms_residual.o \ +scale_fft_3drism.o \ +scale_fft_lauerism.o \ +solute.o \ +solvation_3drism.o \ +solvation_esm.o \ +solvation_force.o \ +solvation_lauerism.o \ +solvation_pbc.o \ +solvation_stress.o \ +solvavg.o \ +solvmol.o \ +summary_1drism.o \ +summary_3drism.o \ +suscept_g0.o \ +suscept_laue.o \ +suscept_laueint.o \ +suscept_vv.o \ +write_rism_type.o \ +xml_io_rism.o + +# list of subroutines and functions (not modules) previously found in flib/clib + +OBJS = \ +atom_weight.o \ +capital.o \ +cryst_to_car.o \ +expint.o \ +generate_k_along_lines.o \ +has_xml.o \ +inpfile.o \ +int_to_char.o \ +latgen.o \ +linpack.o \ +matches.o \ +plot_io.o \ +radial_gradients.o \ +rgen.o \ +recips.o \ +remove_tot_torque.o \ +sort.o \ +trimcheck.o \ +test_input_file.o \ +date_and_tim.o \ +volume.o \ +wgauss.o \ +w0gauss.o \ +w1gauss.o \ +deviatoric.o \ +customize_signals.o \ +qmmm_aux.o \ +sockets.o \ +stack.o + +# GPU versions of modules +MODULES += \ + wavefunctions_gpu.o \ + becmod_gpu.o \ + becmod_subs_gpu.o \ + cuda_subroutines.o \ + random_numbers_gpu.o + +TLDEPS= libfox libla libfft libutil libmbd librxc libupf + +all : libqemod.a + +libqemod.a: $(MODULES) $(OBJS) $(RISMLIB) + $(AR) $(ARFLAGS) $@ $? + $(RANLIB) $@ + +tldeps : + if test -n "$(TLDEPS)" ; then \ + ( cd ../.. ; $(MAKE) $(TLDEPS) || exit 1 ) ; fi + + +clean : + - /bin/rm -f *.o *.a *.d *.i *~ *_tmp.f90 *.mod *.L plumed.f90 + +plumed.f90: + cp $(PLUMED_FORTRAN) plumed.f90 + +include make.depend diff --git a/patches/qespresso-7.2.diff/Modules/Makefile.preplumed b/patches/qespresso-7.2.diff/Modules/Makefile.preplumed new file mode 100644 index 0000000000..acc568ec27 --- /dev/null +++ b/patches/qespresso-7.2.diff/Modules/Makefile.preplumed @@ -0,0 +1,244 @@ +#/a Makefile for Modules + +include ../make.inc + +# location of needed modules +MODFLAGS=$(BASEMOD_FLAGS) + +# list of modules + +MODULES = \ +additional_kpoints.o \ +autopilot.o \ +basic_algebra_routines.o \ +becmod.o \ +bfgs_module.o \ +bspline.o \ +bz_form.o \ +cell_base.o \ +check_stop.o \ +command_line_options.o \ +compute_dipole.o \ +constants.o \ +constraints_module.o \ +control_flags.o \ +coulomb_vcut.o \ +dist.o \ +electrons_base.o \ +environ_base_module.o \ +environment.o \ +extffield.o \ +fd_gradient.o \ +fft_base.o \ +fft_rho.o \ +fft_wave.o \ +fsockets.o \ +funct.o \ +generate_function.o \ +gradutils.o \ +gvecw.o \ +input_parameters.o \ +invmat.o \ +io_files.o \ +io_global.o \ +ions_base.o \ +kind.o \ +lmdif.o \ +mdiis.o \ +mm_dispersion.o \ +mp_bands.o \ +mp_exx.o \ +mp_global.o \ +mp_images.o \ +mp_pools.o \ +mp_wave.o \ +mp_world.o \ +noncol.o \ +open_close_input_file.o \ +parameters.o \ +parser.o \ +plugin_flags.o \ +plugin_arguments.o \ +plugin_variables.o \ +pw_dot.o \ +qmmm.o \ +random_numbers.o \ +read_cards.o \ +read_input.o \ +read_namelists.o \ +read_pseudo.o \ +recvec.o \ +recvec_subs.o \ +run_info.o \ +space_group.o \ +set_para_diag.o \ +set_signal.o \ +set_vdw_corr.o \ +setqf.o \ +timestep.o\ +tsvdw.o\ +mbdlib.o\ +version.o \ +wannier_gw.o\ +wannier_new.o \ +wavefunctions.o \ +ws_base.o \ +xc_vdW_DF.o \ +xc_rVV10.o \ +io_base.o \ +qes_types_module.o \ +qes_libs_module.o \ +qes_write_module.o \ +qes_read_module.o \ +qes_reset_module.o \ +qes_init_module.o \ +qes_bcast_module.o \ +qexsd.o \ +qexsd_copy.o \ +qexsd_init.o \ +qexsd_input.o \ +hdf5_qe.o\ +qeh5_module.o\ +fox_init_module.o \ +xsf.o \ +wyckoff.o \ +wypos.o \ +zvscal.o \ +wave_gauge.o + +# list of RISM's modules + +RISMLIB = \ +allocate_fft_3drism.o \ +chempot.o \ +chempot_lauerism.o \ +closure.o \ +corrdipole_laue.o \ +correctat0_vv.o \ +corrgxy0_laue.o \ +cryst_to_car_2d.o \ +data_structure_3drism.o \ +do_1drism.o \ +do_3drism.o \ +do_lauerism.o \ +eqn_1drism.o \ +eqn_3drism.o \ +eqn_lauedipole.o \ +eqn_lauegxy0.o \ +eqn_lauelong.o \ +eqn_lauerism.o \ +eqn_laueshort.o \ +eqn_lauevoid.o \ +err_rism.o \ +guess_3drism.o \ +init_1drism.o \ +init_3drism.o \ +input_1drism.o \ +input_3drism.o \ +io_rism_xml.o \ +lauefft.o \ +lauefft_subs.o \ +lj_forcefield.o \ +lj_solute.o \ +molecorr_vv.o \ +molebridge_vv.o \ +molecule_const.o \ +molecule_types.o \ +mp_rism.o \ +mp_swap_ax_rism.o \ +normalize_lauerism.o \ +plot_rism.o \ +potential_3drism.o \ +potential_esm.o \ +potential_vv.o \ +print_chempot_3drism.o \ +print_chempot_lauerism.o \ +print_chempot_vv.o \ +print_corr_vv.o \ +print_solvavg.o \ +radfft.o \ +read_mol.o \ +read_solv.o \ +recvec_3drism.o \ +rism.o \ +rism1d_facade.o \ +rism3d_facade.o \ +rms_residual.o \ +scale_fft_3drism.o \ +scale_fft_lauerism.o \ +solute.o \ +solvation_3drism.o \ +solvation_esm.o \ +solvation_force.o \ +solvation_lauerism.o \ +solvation_pbc.o \ +solvation_stress.o \ +solvavg.o \ +solvmol.o \ +summary_1drism.o \ +summary_3drism.o \ +suscept_g0.o \ +suscept_laue.o \ +suscept_laueint.o \ +suscept_vv.o \ +write_rism_type.o \ +xml_io_rism.o + +# list of subroutines and functions (not modules) previously found in flib/clib + +OBJS = \ +atom_weight.o \ +capital.o \ +cryst_to_car.o \ +expint.o \ +generate_k_along_lines.o \ +has_xml.o \ +inpfile.o \ +int_to_char.o \ +latgen.o \ +linpack.o \ +matches.o \ +plot_io.o \ +radial_gradients.o \ +rgen.o \ +recips.o \ +remove_tot_torque.o \ +sort.o \ +trimcheck.o \ +test_input_file.o \ +date_and_tim.o \ +volume.o \ +wgauss.o \ +w0gauss.o \ +w1gauss.o \ +deviatoric.o \ +customize_signals.o \ +qmmm_aux.o \ +sockets.o \ +stack.o + +# GPU versions of modules +MODULES += \ + wavefunctions_gpu.o \ + becmod_gpu.o \ + becmod_subs_gpu.o \ + cuda_subroutines.o \ + random_numbers_gpu.o + +TLDEPS= libfox libla libfft libutil libmbd librxc libupf + +all : libqemod.a + +libqemod.a: $(MODULES) $(OBJS) $(RISMLIB) + $(AR) $(ARFLAGS) $@ $? + $(RANLIB) $@ + +tldeps : + if test -n "$(TLDEPS)" ; then \ + ( cd ../.. ; $(MAKE) $(TLDEPS) || exit 1 ) ; fi + + +clean : + - /bin/rm -f *.o *.a *.d *.i *~ *_tmp.f90 *.mod *.L + +include make.depend diff --git a/patches/qespresso-7.2.diff/PW/src/forces.f90 b/patches/qespresso-7.2.diff/PW/src/forces.f90 new file mode 100644 index 0000000000..588222cb3d --- /dev/null +++ b/patches/qespresso-7.2.diff/PW/src/forces.f90 @@ -0,0 +1,532 @@ +! +! Copyright (C) 2001-2022 Quantum ESPRESSO 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 forces() + !---------------------------------------------------------------------------- + !! This routine is a driver routine which computes the forces + !! acting on the atoms. The complete expression of the forces + !! contains many parts which are computed by different routines: + ! + !! - force_lc: local potential contribution + !! - force_us: non-local potential contribution + !! - (esm_)force_ew: (ESM) electrostatic ewald term + !! - force_cc: nonlinear core correction contribution + !! - force_corr: correction term for incomplete self-consistency + !! - force_hub: contribution due to the Hubbard term; + !! - force_london: Grimme DFT+D dispersion forces + !! - force_d3: Grimme-D3 (DFT-D3) dispersion forces + !! - force_xdm: XDM dispersion forces + !! - more terms from external electric fields, Martyna-Tuckerman, etc. + !! - force_sol: contribution due to 3D-RISM + ! + USE kinds, ONLY : DP + USE io_global, ONLY : stdout + USE cell_base, ONLY : at, bg, alat, omega + USE ions_base, ONLY : nat, ntyp => nsp,nsp, ityp, tau, zv, amass, extfor, atm + USE gvect, ONLY : ngm, gstart, ngl, igtongl, g, gg, gcutm + USE lsda_mod, ONLY : nspin + USE symme, ONLY : symvector + USE vlocal, ONLY : strf, vloc + USE force_mod, ONLY : force, sumfor + USE scf, ONLY : rho + USE ions_base, ONLY : if_pos + USE ldaU, ONLY : lda_plus_u, Hubbard_projectors + USE extfield, ONLY : tefield, forcefield, gate, forcegate, relaxz + USE control_flags, ONLY : gamma_only, remove_rigid_rot, textfor, & + iverbosity, llondon, ldftd3, lxdm, ts_vdw, & + mbd_vdw, lforce => tprnfor, istep + USE bp, ONLY : lelfield, gdir, l3dstring, efield_cart, & + efield_cry,efield + USE uspp, ONLY : okvan + USE martyna_tuckerman, ONLY : do_comp_mt, wg_corr_force + USE london_module, ONLY : force_london + USE dftd3_api, ONLY : get_atomic_number + USE dftd3_qe, ONLY : dftd3_pbc_gdisp, dftd3 + + USE xdm_module, ONLY : force_xdm + USE tsvdw_module, ONLY : FtsvdW + USE libmbd_interface, ONLY : FmbdvdW + USE esm, ONLY : do_comp_esm, esm_bc, esm_force_ew + USE qmmm, ONLY : qmmm_mode + USE rism_module, ONLY : lrism, force_rism + USE extffield, ONLY : apply_extffield_PW + USE input_parameters, ONLY : nextffield + ! +#if defined(__CUDA) + USE device_fbuff_m, ONLY : dev_buf +#endif + ! +#if defined (__ENVIRON) + USE plugin_flags, ONLY : use_environ + USE environ_base_module, ONLY : calc_environ_force + USE environ_pw_module, ONLY : is_ms_gcs, run_ms_gcs +#endif +#if defined (__OSCDFT) + USE plugin_flags, ONLY : use_oscdft + USE oscdft_base, ONLY : oscdft_ctx + USE oscdft_forces_subs, ONLY : oscdft_apply_forces, oscdft_print_forces +#endif + ! + IMPLICIT NONE + ! + REAL(DP), ALLOCATABLE :: forcenl(:,:), & + forcelc(:,:), & + forcecc(:,:), & + forceion(:,:), & + force_disp(:,:), & + force_d3(:,:), & + force_disp_xdm(:,:), & + force_mt(:,:), & + forcescc(:,:), & + forces_bp_efield(:,:),& + forceh(:,:), & + force_sol(:,:) + ! nonlocal, local, core-correction, ewald, scf correction terms, and hubbard + ! + ! aux is used to store a possible additional density + ! now defined in real space + ! + COMPLEX(DP), ALLOCATABLE :: auxg(:), auxr(:) + ! + REAL(DP) :: sumscf, sum_mm + REAL(DP), PARAMETER :: eps = 1.e-12_dp + INTEGER :: ipol, na + ! counter on polarization + ! counter on atoms + ! + REAL(DP) :: latvecs(3,3) + INTEGER :: atnum(1:nat) + REAL(DP) :: stress_dftd3(3,3) + ! + INTEGER :: ierr + ! + force(:,:) = 0.D0 + ! + ! Early return if all forces to be set to zero + ! + IF ( ALL( if_pos == 0 ) ) RETURN + ! + CALL start_clock( 'forces' ) +#if defined(__CUDA) + ! Cleanup scratch space used in previous SCF iterations. + ! This will reduce memory footprint. + CALL dev_buf%reinit(ierr) + IF (ierr .ne. 0) CALL infomsg('forces', 'Cannot reset GPU buffers! Some buffers still locked.') +#endif + ! + ! + ALLOCATE( forcenl(3,nat), forcelc(3,nat), forcecc(3,nat), & + forceh(3,nat), forceion(3,nat), forcescc(3,nat) ) + ! + forcescc(:,:) = 0.D0 + forceh(:,:) = 0.D0 + ! + ! ... The nonlocal contribution is computed here + ! + call start_clock('frc_us') + CALL force_us( forcenl ) + call stop_clock('frc_us') + ! + ! ... The local contribution + ! + call start_clock('frc_lc') + CALL force_lc( nat, tau, ityp, ntyp, alat, omega, ngm, ngl, igtongl, & + g, rho%of_r(:,1), gstart, gamma_only, vloc, forcelc ) + call stop_clock('frc_lc') + ! + ! ... The NLCC contribution + ! + call start_clock('frc_cc') + CALL force_cc( forcecc ) + call stop_clock('frc_cc') + + ! ... The Hubbard contribution + ! (included by force_us if using beta as local projectors) + ! + IF ( lda_plus_u .AND. Hubbard_projectors.NE.'pseudo' ) CALL force_hub( forceh ) + ! + ! ... The ionic contribution is computed here + ! + IF( do_comp_esm ) THEN + CALL esm_force_ew( forceion ) + ELSE + CALL force_ew( alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, & + gg, ngm, gstart, gamma_only, gcutm, strf, forceion ) + ENDIF + ! + ! ... the semi-empirical dispersion correction + ! + IF ( llondon ) THEN + ! + ALLOCATE( force_disp(3,nat) ) + force_disp(:,:) = 0.0_DP + force_disp = force_london( alat , nat , ityp , at , bg , tau ) + ! + ENDIF + ! + ! ... The Grimme-D3 dispersion correction + ! + IF ( ldftd3 ) THEN + ! + CALL start_clock('force_dftd3') + ALLOCATE( force_d3(3, nat) ) + force_d3(:,:) = 0.0_DP + latvecs(:,:) = at(:,:)*alat + tau(:,:) = tau(:,:)*alat + atnum(:) = get_atomic_number(atm(ityp(:))) + CALL dftd3_pbc_gdisp( dftd3, tau, atnum, latvecs, & + force_d3, stress_dftd3 ) + force_d3 = -2.d0*force_d3 + tau(:,:) = tau(:,:)/alat + CALL stop_clock('force_dftd3') + ENDIF + ! + ! + IF (lxdm) THEN + ALLOCATE( force_disp_xdm(3,nat) ) + force_disp_xdm = 0._dp + force_disp_xdm = force_xdm(nat) + ENDIF + ! + ! ... The SCF contribution + ! + call start_clock('frc_scc') +#if defined(__CUDA) + ! Cleanup scratch space again, next subroutines uses a lot of memory. + ! In an ideal world this should be done only if really needed (TODO). + CALL dev_buf%reinit(ierr) + IF (ierr .ne. 0) CALL errore('forces', 'Cannot reset GPU buffers! Buffers still locked: ', abs(ierr)) +#endif + ! + CALL force_corr( forcescc ) + call stop_clock('frc_scc') + ! + IF (do_comp_mt) THEN + ! + ALLOCATE( force_mt(3,nat) ) + CALL wg_corr_force( .TRUE., omega, nat, ntyp, ityp, ngm, g, tau, zv, strf, & + rho%of_g(:,1), force_mt ) + ENDIF + ! + ! ... The solvation contribution (3D-RISM) + ! + IF (lrism) THEN + ALLOCATE ( force_sol ( 3 , nat ) ) + CALL force_rism( force_sol ) + END IF + ! + ! ... call void routine for user define/ plugin patches on internal forces + ! +#if defined(__LEGACY_PLUGINS) + CALL plugin_int_forces() +#endif +#if defined (__ENVIRON) + IF (use_environ) CALL calc_environ_force(force) +#endif +#if defined (__OSCDFT) + IF (use_oscdft) CALL oscdft_apply_forces(oscdft_ctx) +#endif + ! + ! ... Berry's phase electric field terms + ! + IF (lelfield) THEN + ALLOCATE( forces_bp_efield(3,nat) ) + forces_bp_efield(:,:) = 0.d0 + IF (.NOT.l3dstring) THEN + IF (okvan) CALL forces_us_efield( forces_bp_efield, gdir, efield ) + CALL forces_ion_efield( forces_bp_efield, gdir, efield ) + ELSE + IF (okvan) THEN + DO ipol = 1, 3 + CALL forces_us_efield( forces_bp_efield, ipol, efield_cry(ipol) ) + ENDDO + ENDIF + DO ipol = 1, 3 + CALL forces_ion_efield( forces_bp_efield, ipol, efield_cart(ipol) ) + ENDDO + ENDIF + ENDIF + ! + ! ... here we sum all the contributions and compute the total force acting + ! ... on the crystal + ! + DO ipol = 1, 3 + ! + sumfor = 0.D0 + ! + DO na = 1, nat + ! + force(ipol,na) = force(ipol,na) + & + forcenl(ipol,na) + & + forceion(ipol,na) + & + forcelc(ipol,na) + & + forcecc(ipol,na) + & + forceh(ipol,na) + & + forcescc(ipol,na) + ! + IF ( llondon ) force(ipol,na) = force(ipol,na) + force_disp(ipol,na) + IF ( ldftd3 ) force(ipol,na) = force(ipol,na) + force_d3(ipol,na) + IF ( lxdm ) force(ipol,na) = force(ipol,na) + force_disp_xdm(ipol,na) + ! factor 2 converts from Ha to Ry a.u. + ! the IF condition is to avoid double counting + IF ( mbd_vdw ) THEN + force(ipol, na) = force(ipol, na) + 2.0_dp*FmbdvdW(ipol, na) + ELSE IF ( ts_vdw ) THEN + force(ipol, na) = force(ipol, na) + 2.0_dp*FtsvdW(ipol, na) + ENDIF + IF ( tefield ) force(ipol,na) = force(ipol,na) + forcefield(ipol,na) + IF ( gate ) force(ipol,na) = force(ipol,na) + forcegate(ipol,na) ! TB + IF (lelfield) force(ipol,na) = force(ipol,na) + forces_bp_efield(ipol,na) + IF (do_comp_mt) force(ipol,na) = force(ipol,na) + force_mt(ipol,na) + IF ( lrism ) force(ipol,na) = force(ipol,na) + force_sol(ipol,na) + ! + sumfor = sumfor + force(ipol,na) + ! + ENDDO + ! + !TB + IF ((gate.AND.relaxz).AND.(ipol==3)) WRITE( stdout, '("Total force in z direction = 0 disabled")') + ! + IF ( (do_comp_esm .AND. ( esm_bc /= 'pbc' )).OR.(gate.AND.relaxz) ) THEN + ! + ! ... impose total force along xy = 0 + ! + DO na = 1, nat + IF ( ipol /= 3) force(ipol,na) = force(ipol,na) & + - sumfor / DBLE( nat ) + ENDDO + ! + ELSEIF ( qmmm_mode < 0 ) THEN + ! + ! ... impose total force = 0 except in a QM-MM calculation + ! + DO na = 1, nat + force(ipol,na) = force(ipol,na) - sumfor / DBLE( nat ) + ENDDO + ! + ENDIF + ! + ENDDO + ! + ! ... call run_extffield to apply external force fields on ions + ! + IF ( nextffield > 0 ) THEN + tau(:,:) = tau(:,:)*alat + CALL apply_extffield_PW(istep,nextffield,tau,force) + tau(:,:) = tau(:,:)/alat + END IF + ! + ! ... resymmetrize (should not be needed, but ...) + ! + CALL symvector( nat, force ) + ! + IF ( remove_rigid_rot ) & + CALL remove_tot_torque( nat, tau, amass(ityp(:)), force ) + ! + IF( textfor ) force(:,:) = force(:,:) + extfor(:,:) + ! + ! ... call void routine for user define/ plugin patches on external forces + ! +#if defined(__LEGACY_PLUGINS) + CALL plugin_ext_forces() +#endif +#if defined (__ENVIRON) + IF (use_environ) THEN + IF (is_ms_gcs()) CALL run_ms_gcs() + END IF +#endif + ! + ! ... write on output the forces + ! + WRITE( stdout, '(/,5x,"Forces acting on atoms (cartesian axes, Ry/au):", / )') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), force(:,na) + ENDDO + ! + ! ... forces on fixed coordinates are set to zero ( C.S. 15/10/2003 ) + ! + force(:,:) = force(:,:) * DBLE( if_pos ) + forcescc(:,:) = forcescc(:,:) * DBLE( if_pos ) + ! +!civn +! IF ( iverbosity > 0 ) THEN + IF ( .true. ) THEN +! + IF ( do_comp_mt ) THEN + WRITE( stdout, '(5x,"The Martyna-Tuckerman correction term to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( force_mt(ipol,na), ipol = 1, 3 ) + ENDDO + END IF + ! + WRITE( stdout, '(5x,"The non-local contrib. to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcenl(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The ionic contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forceion(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The local contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcelc(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The core correction contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcecc(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The Hubbard contrib. to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forceh(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The SCF correction term to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcescc(ipol,na), ipol = 1, 3 ) + ENDDO + ! + IF ( llondon) THEN + WRITE( stdout, '(/,5x,"Dispersion contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_disp(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + IF ( ldftd3 ) THEN + WRITE( stdout, '(/,5x,"DFT-D3 dispersion contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_d3(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + IF (lxdm) THEN + WRITE( stdout, '(/,5x,"XDM contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_disp_xdm(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + ! again, as above, if condition is to avoid redundant printing + IF ( mbd_vdw ) THEN + WRITE( stdout, '(/,5x, "MBD contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (2.0d0*FmbdvdW(ipol, na), ipol = 1, 3) + ENDDO + ELSE IF ( ts_vdw ) THEN + WRITE( stdout, '(/,5x, "TS-VDW contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (2.0d0*FtsvdW(ipol, na), ipol = 1, 3) + ENDDO + ENDIF + + ! + ! TB gate forces + IF ( gate ) THEN + WRITE( stdout, '(/,5x,"Gate contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (forcegate(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + IF ( lrism ) THEN + WRITE( stdout, '(/,5x,"3D-RISM Solvation contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_sol(ipol,na), ipol = 1, 3) + END DO + END IF + ! + END IF +#if defined (__OSCDFT) + IF (use_oscdft) CALL oscdft_print_forces(oscdft_ctx) +#endif + ! + sumfor = 0.D0 + sumscf = 0.D0 + ! + DO na = 1, nat + ! + sumfor = sumfor + force(1,na)**2 + force(2,na)**2 + force(3,na)**2 + sumscf = sumscf + forcescc(1,na)**2 + forcescc(2,na)**2+ forcescc(3,na)**2 + ! + ENDDO + ! + sumfor = SQRT( sumfor ) + sumscf = SQRT( sumscf ) + ! + WRITE( stdout, '(/5x,"Total force = ",F12.6,5X, & + & "Total SCF correction = ",F12.6)') sumfor, sumscf + ! + IF ( llondon .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_disp(1,na)**2 + force_disp(2,na)**2 + force_disp(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "Total Dispersion Force = ",F12.6)') sum_mm + ! + END IF + ! + IF ( ldftd3 .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_d3(1,na)**2 + force_d3(2,na)**2 + force_d3(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "DFT-D3 dispersion Force = ",F12.6)') sum_mm + ! + END IF + ! + IF ( lxdm .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_disp_xdm(1,na)**2 + force_disp_xdm(2,na)**2 + force_disp_xdm(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "Total XDM Force = ",F12.6)') sum_mm + ! + END IF + ! + IF ( lrism .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_sol(1,na)**2 + force_sol(2,na)**2 + force_sol(3,na)**2 + END DO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "Total 3D-RISM Solvation Force = ",F12.6)') sum_mm + ! + END IF + ! + DEALLOCATE( forcenl, forcelc, forcecc, forceh, forceion, forcescc ) + IF ( llondon ) DEALLOCATE( force_disp ) + IF ( ldftd3 ) DEALLOCATE( force_d3 ) + IF ( lxdm ) DEALLOCATE( force_disp_xdm ) + IF ( lelfield ) DEALLOCATE( forces_bp_efield ) + IF ( lrism ) DEALLOCATE( force_sol ) + IF(ALLOCATED(force_mt)) DEALLOCATE( force_mt ) + ! + ! FIXME: what is the following line good for? + ! + lforce = .TRUE. + ! + CALL stop_clock( 'forces' ) + ! + IF ( ( sumfor < 10.D0*sumscf ) .AND. ( sumfor > nat*eps ) ) & + WRITE( stdout,'(5x,"SCF correction compared to forces is large: ", & + & "reduce conv_thr to get better values")') + RETURN + ! +9035 FORMAT(5X,'atom ',I4,' type ',I2,' force = ',3F14.8) + ! +END SUBROUTINE forces diff --git a/patches/qespresso-7.2.diff/PW/src/forces.f90.preplumed b/patches/qespresso-7.2.diff/PW/src/forces.f90.preplumed new file mode 100644 index 0000000000..96f6bfc6b6 --- /dev/null +++ b/patches/qespresso-7.2.diff/PW/src/forces.f90.preplumed @@ -0,0 +1,535 @@ +! +! Copyright (C) 2001-2022 Quantum ESPRESSO 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 forces() + !---------------------------------------------------------------------------- + !! This routine is a driver routine which computes the forces + !! acting on the atoms. The complete expression of the forces + !! contains many parts which are computed by different routines: + ! + !! - force_lc: local potential contribution + !! - force_us: non-local potential contribution + !! - (esm_)force_ew: (ESM) electrostatic ewald term + !! - force_cc: nonlinear core correction contribution + !! - force_corr: correction term for incomplete self-consistency + !! - force_hub: contribution due to the Hubbard term; + !! - force_london: Grimme DFT+D dispersion forces + !! - force_d3: Grimme-D3 (DFT-D3) dispersion forces + !! - force_xdm: XDM dispersion forces + !! - more terms from external electric fields, Martyna-Tuckerman, etc. + !! - force_sol: contribution due to 3D-RISM + ! + USE kinds, ONLY : DP + USE io_global, ONLY : stdout + USE cell_base, ONLY : at, bg, alat, omega + USE ions_base, ONLY : nat, ntyp => nsp,nsp, ityp, tau, zv, amass, extfor, atm + USE gvect, ONLY : ngm, gstart, ngl, igtongl, g, gg, gcutm + USE lsda_mod, ONLY : nspin + USE symme, ONLY : symvector + USE vlocal, ONLY : strf, vloc + USE force_mod, ONLY : force, sumfor + USE scf, ONLY : rho + USE ions_base, ONLY : if_pos + USE ldaU, ONLY : lda_plus_u, Hubbard_projectors + USE extfield, ONLY : tefield, forcefield, gate, forcegate, relaxz + USE control_flags, ONLY : gamma_only, remove_rigid_rot, textfor, & + iverbosity, llondon, ldftd3, lxdm, ts_vdw, & + mbd_vdw, lforce => tprnfor, istep + USE bp, ONLY : lelfield, gdir, l3dstring, efield_cart, & + efield_cry,efield + USE uspp, ONLY : okvan + USE martyna_tuckerman, ONLY : do_comp_mt, wg_corr_force + USE london_module, ONLY : force_london + USE dftd3_api, ONLY : get_atomic_number + USE dftd3_qe, ONLY : dftd3_pbc_gdisp, dftd3 + + USE xdm_module, ONLY : force_xdm + USE tsvdw_module, ONLY : FtsvdW + USE libmbd_interface, ONLY : FmbdvdW + USE esm, ONLY : do_comp_esm, esm_bc, esm_force_ew + USE qmmm, ONLY : qmmm_mode + USE rism_module, ONLY : lrism, force_rism + USE extffield, ONLY : apply_extffield_PW + USE input_parameters, ONLY : nextffield + ! +#if defined(__CUDA) + USE device_fbuff_m, ONLY : dev_buf +#endif + ! +#if defined (__ENVIRON) + USE plugin_flags, ONLY : use_environ + USE environ_base_module, ONLY : calc_environ_force + USE environ_pw_module, ONLY : is_ms_gcs, run_ms_gcs +#endif +#if defined (__OSCDFT) + USE plugin_flags, ONLY : use_oscdft + USE oscdft_base, ONLY : oscdft_ctx + USE oscdft_forces_subs, ONLY : oscdft_apply_forces, oscdft_print_forces +#endif +#if defined(__LEGACY_PLUGINS) + USE plugin_flags, ONLY : plugin_ext_forces, plugin_int_forces +#endif + ! + IMPLICIT NONE + ! + REAL(DP), ALLOCATABLE :: forcenl(:,:), & + forcelc(:,:), & + forcecc(:,:), & + forceion(:,:), & + force_disp(:,:), & + force_d3(:,:), & + force_disp_xdm(:,:), & + force_mt(:,:), & + forcescc(:,:), & + forces_bp_efield(:,:),& + forceh(:,:), & + force_sol(:,:) + ! nonlocal, local, core-correction, ewald, scf correction terms, and hubbard + ! + ! aux is used to store a possible additional density + ! now defined in real space + ! + COMPLEX(DP), ALLOCATABLE :: auxg(:), auxr(:) + ! + REAL(DP) :: sumscf, sum_mm + REAL(DP), PARAMETER :: eps = 1.e-12_dp + INTEGER :: ipol, na + ! counter on polarization + ! counter on atoms + ! + REAL(DP) :: latvecs(3,3) + INTEGER :: atnum(1:nat) + REAL(DP) :: stress_dftd3(3,3) + ! + INTEGER :: ierr + ! + force(:,:) = 0.D0 + ! + ! Early return if all forces to be set to zero + ! + IF ( ALL( if_pos == 0 ) ) RETURN + ! + CALL start_clock( 'forces' ) +#if defined(__CUDA) + ! Cleanup scratch space used in previous SCF iterations. + ! This will reduce memory footprint. + CALL dev_buf%reinit(ierr) + IF (ierr .ne. 0) CALL infomsg('forces', 'Cannot reset GPU buffers! Some buffers still locked.') +#endif + ! + ! + ALLOCATE( forcenl(3,nat), forcelc(3,nat), forcecc(3,nat), & + forceh(3,nat), forceion(3,nat), forcescc(3,nat) ) + ! + forcescc(:,:) = 0.D0 + forceh(:,:) = 0.D0 + ! + ! ... The nonlocal contribution is computed here + ! + call start_clock('frc_us') + CALL force_us( forcenl ) + call stop_clock('frc_us') + ! + ! ... The local contribution + ! + call start_clock('frc_lc') + CALL force_lc( nat, tau, ityp, ntyp, alat, omega, ngm, ngl, igtongl, & + g, rho%of_r(:,1), gstart, gamma_only, vloc, forcelc ) + call stop_clock('frc_lc') + ! + ! ... The NLCC contribution + ! + call start_clock('frc_cc') + CALL force_cc( forcecc ) + call stop_clock('frc_cc') + + ! ... The Hubbard contribution + ! (included by force_us if using beta as local projectors) + ! + IF ( lda_plus_u .AND. Hubbard_projectors.NE.'pseudo' ) CALL force_hub( forceh ) + ! + ! ... The ionic contribution is computed here + ! + IF( do_comp_esm ) THEN + CALL esm_force_ew( forceion ) + ELSE + CALL force_ew( alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, & + gg, ngm, gstart, gamma_only, gcutm, strf, forceion ) + ENDIF + ! + ! ... the semi-empirical dispersion correction + ! + IF ( llondon ) THEN + ! + ALLOCATE( force_disp(3,nat) ) + force_disp(:,:) = 0.0_DP + force_disp = force_london( alat , nat , ityp , at , bg , tau ) + ! + ENDIF + ! + ! ... The Grimme-D3 dispersion correction + ! + IF ( ldftd3 ) THEN + ! + CALL start_clock('force_dftd3') + ALLOCATE( force_d3(3, nat) ) + force_d3(:,:) = 0.0_DP + latvecs(:,:) = at(:,:)*alat + tau(:,:) = tau(:,:)*alat + atnum(:) = get_atomic_number(atm(ityp(:))) + CALL dftd3_pbc_gdisp( dftd3, tau, atnum, latvecs, & + force_d3, stress_dftd3 ) + force_d3 = -2.d0*force_d3 + tau(:,:) = tau(:,:)/alat + CALL stop_clock('force_dftd3') + ENDIF + ! + ! + IF (lxdm) THEN + ALLOCATE( force_disp_xdm(3,nat) ) + force_disp_xdm = 0._dp + force_disp_xdm = force_xdm(nat) + ENDIF + ! + ! ... The SCF contribution + ! + call start_clock('frc_scc') +#if defined(__CUDA) + ! Cleanup scratch space again, next subroutines uses a lot of memory. + ! In an ideal world this should be done only if really needed (TODO). + CALL dev_buf%reinit(ierr) + IF (ierr .ne. 0) CALL errore('forces', 'Cannot reset GPU buffers! Buffers still locked: ', abs(ierr)) +#endif + ! + CALL force_corr( forcescc ) + call stop_clock('frc_scc') + ! + IF (do_comp_mt) THEN + ! + ALLOCATE( force_mt(3,nat) ) + CALL wg_corr_force( .TRUE., omega, nat, ntyp, ityp, ngm, g, tau, zv, strf, & + rho%of_g(:,1), force_mt ) + ENDIF + ! + ! ... The solvation contribution (3D-RISM) + ! + IF (lrism) THEN + ALLOCATE ( force_sol ( 3 , nat ) ) + CALL force_rism( force_sol ) + END IF + ! + ! ... call void routine for user define/ plugin patches on internal forces + ! +#if defined(__LEGACY_PLUGINS) + CALL plugin_int_forces() +#endif +#if defined (__ENVIRON) + IF (use_environ) CALL calc_environ_force(force) +#endif +#if defined (__OSCDFT) + IF (use_oscdft) CALL oscdft_apply_forces(oscdft_ctx) +#endif + ! + ! ... Berry's phase electric field terms + ! + IF (lelfield) THEN + ALLOCATE( forces_bp_efield(3,nat) ) + forces_bp_efield(:,:) = 0.d0 + IF (.NOT.l3dstring) THEN + IF (okvan) CALL forces_us_efield( forces_bp_efield, gdir, efield ) + CALL forces_ion_efield( forces_bp_efield, gdir, efield ) + ELSE + IF (okvan) THEN + DO ipol = 1, 3 + CALL forces_us_efield( forces_bp_efield, ipol, efield_cry(ipol) ) + ENDDO + ENDIF + DO ipol = 1, 3 + CALL forces_ion_efield( forces_bp_efield, ipol, efield_cart(ipol) ) + ENDDO + ENDIF + ENDIF + ! + ! ... here we sum all the contributions and compute the total force acting + ! ... on the crystal + ! + DO ipol = 1, 3 + ! + sumfor = 0.D0 + ! + DO na = 1, nat + ! + force(ipol,na) = force(ipol,na) + & + forcenl(ipol,na) + & + forceion(ipol,na) + & + forcelc(ipol,na) + & + forcecc(ipol,na) + & + forceh(ipol,na) + & + forcescc(ipol,na) + ! + IF ( llondon ) force(ipol,na) = force(ipol,na) + force_disp(ipol,na) + IF ( ldftd3 ) force(ipol,na) = force(ipol,na) + force_d3(ipol,na) + IF ( lxdm ) force(ipol,na) = force(ipol,na) + force_disp_xdm(ipol,na) + ! factor 2 converts from Ha to Ry a.u. + ! the IF condition is to avoid double counting + IF ( mbd_vdw ) THEN + force(ipol, na) = force(ipol, na) + 2.0_dp*FmbdvdW(ipol, na) + ELSE IF ( ts_vdw ) THEN + force(ipol, na) = force(ipol, na) + 2.0_dp*FtsvdW(ipol, na) + ENDIF + IF ( tefield ) force(ipol,na) = force(ipol,na) + forcefield(ipol,na) + IF ( gate ) force(ipol,na) = force(ipol,na) + forcegate(ipol,na) ! TB + IF (lelfield) force(ipol,na) = force(ipol,na) + forces_bp_efield(ipol,na) + IF (do_comp_mt) force(ipol,na) = force(ipol,na) + force_mt(ipol,na) + IF ( lrism ) force(ipol,na) = force(ipol,na) + force_sol(ipol,na) + ! + sumfor = sumfor + force(ipol,na) + ! + ENDDO + ! + !TB + IF ((gate.AND.relaxz).AND.(ipol==3)) WRITE( stdout, '("Total force in z direction = 0 disabled")') + ! + IF ( (do_comp_esm .AND. ( esm_bc /= 'pbc' )).OR.(gate.AND.relaxz) ) THEN + ! + ! ... impose total force along xy = 0 + ! + DO na = 1, nat + IF ( ipol /= 3) force(ipol,na) = force(ipol,na) & + - sumfor / DBLE( nat ) + ENDDO + ! + ELSEIF ( qmmm_mode < 0 ) THEN + ! + ! ... impose total force = 0 except in a QM-MM calculation + ! + DO na = 1, nat + force(ipol,na) = force(ipol,na) - sumfor / DBLE( nat ) + ENDDO + ! + ENDIF + ! + ENDDO + ! + ! ... call run_extffield to apply external force fields on ions + ! + IF ( nextffield > 0 ) THEN + tau(:,:) = tau(:,:)*alat + CALL apply_extffield_PW(istep,nextffield,tau,force) + tau(:,:) = tau(:,:)/alat + END IF + ! + ! ... resymmetrize (should not be needed, but ...) + ! + CALL symvector( nat, force ) + ! + IF ( remove_rigid_rot ) & + CALL remove_tot_torque( nat, tau, amass(ityp(:)), force ) + ! + IF( textfor ) force(:,:) = force(:,:) + extfor(:,:) + ! + ! ... call void routine for user define/ plugin patches on external forces + ! +#if defined(__LEGACY_PLUGINS) + CALL plugin_ext_forces() +#endif +#if defined (__ENVIRON) + IF (use_environ) THEN + IF (is_ms_gcs()) CALL run_ms_gcs() + END IF +#endif + ! + ! ... write on output the forces + ! + WRITE( stdout, '(/,5x,"Forces acting on atoms (cartesian axes, Ry/au):", / )') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), force(:,na) + ENDDO + ! + ! ... forces on fixed coordinates are set to zero ( C.S. 15/10/2003 ) + ! + force(:,:) = force(:,:) * DBLE( if_pos ) + forcescc(:,:) = forcescc(:,:) * DBLE( if_pos ) + ! +!civn +! IF ( iverbosity > 0 ) THEN + IF ( .true. ) THEN +! + IF ( do_comp_mt ) THEN + WRITE( stdout, '(5x,"The Martyna-Tuckerman correction term to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( force_mt(ipol,na), ipol = 1, 3 ) + ENDDO + END IF + ! + WRITE( stdout, '(5x,"The non-local contrib. to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcenl(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The ionic contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forceion(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The local contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcelc(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The core correction contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcecc(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The Hubbard contrib. to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forceh(ipol,na), ipol = 1, 3 ) + ENDDO + WRITE( stdout, '(5x,"The SCF correction term to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), ( forcescc(ipol,na), ipol = 1, 3 ) + ENDDO + ! + IF ( llondon) THEN + WRITE( stdout, '(/,5x,"Dispersion contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_disp(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + IF ( ldftd3 ) THEN + WRITE( stdout, '(/,5x,"DFT-D3 dispersion contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_d3(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + IF (lxdm) THEN + WRITE( stdout, '(/,5x,"XDM contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_disp_xdm(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + ! again, as above, if condition is to avoid redundant printing + IF ( mbd_vdw ) THEN + WRITE( stdout, '(/,5x, "MBD contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (2.0d0*FmbdvdW(ipol, na), ipol = 1, 3) + ENDDO + ELSE IF ( ts_vdw ) THEN + WRITE( stdout, '(/,5x, "TS-VDW contribution to forces")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (2.0d0*FtsvdW(ipol, na), ipol = 1, 3) + ENDDO + ENDIF + + ! + ! TB gate forces + IF ( gate ) THEN + WRITE( stdout, '(/,5x,"Gate contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (forcegate(ipol,na), ipol = 1, 3) + ENDDO + END IF + ! + IF ( lrism ) THEN + WRITE( stdout, '(/,5x,"3D-RISM Solvation contribution to forces:")') + DO na = 1, nat + WRITE( stdout, 9035) na, ityp(na), (force_sol(ipol,na), ipol = 1, 3) + END DO + END IF + ! + END IF +#if defined (__OSCDFT) + IF (use_oscdft) CALL oscdft_print_forces(oscdft_ctx) +#endif + ! + sumfor = 0.D0 + sumscf = 0.D0 + ! + DO na = 1, nat + ! + sumfor = sumfor + force(1,na)**2 + force(2,na)**2 + force(3,na)**2 + sumscf = sumscf + forcescc(1,na)**2 + forcescc(2,na)**2+ forcescc(3,na)**2 + ! + ENDDO + ! + sumfor = SQRT( sumfor ) + sumscf = SQRT( sumscf ) + ! + WRITE( stdout, '(/5x,"Total force = ",F12.6,5X, & + & "Total SCF correction = ",F12.6)') sumfor, sumscf + ! + IF ( llondon .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_disp(1,na)**2 + force_disp(2,na)**2 + force_disp(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "Total Dispersion Force = ",F12.6)') sum_mm + ! + END IF + ! + IF ( ldftd3 .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_d3(1,na)**2 + force_d3(2,na)**2 + force_d3(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "DFT-D3 dispersion Force = ",F12.6)') sum_mm + ! + END IF + ! + IF ( lxdm .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_disp_xdm(1,na)**2 + force_disp_xdm(2,na)**2 + force_disp_xdm(3,na)**2 + ENDDO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "Total XDM Force = ",F12.6)') sum_mm + ! + END IF + ! + IF ( lrism .AND. iverbosity > 0 ) THEN + ! + sum_mm = 0.D0 + DO na = 1, nat + sum_mm = sum_mm + & + force_sol(1,na)**2 + force_sol(2,na)**2 + force_sol(3,na)**2 + END DO + sum_mm = SQRT( sum_mm ) + WRITE ( stdout, '(/,5x, "Total 3D-RISM Solvation Force = ",F12.6)') sum_mm + ! + END IF + ! + DEALLOCATE( forcenl, forcelc, forcecc, forceh, forceion, forcescc ) + IF ( llondon ) DEALLOCATE( force_disp ) + IF ( ldftd3 ) DEALLOCATE( force_d3 ) + IF ( lxdm ) DEALLOCATE( force_disp_xdm ) + IF ( lelfield ) DEALLOCATE( forces_bp_efield ) + IF ( lrism ) DEALLOCATE( force_sol ) + IF(ALLOCATED(force_mt)) DEALLOCATE( force_mt ) + ! + ! FIXME: what is the following line good for? + ! + lforce = .TRUE. + ! + CALL stop_clock( 'forces' ) + ! + IF ( ( sumfor < 10.D0*sumscf ) .AND. ( sumfor > nat*eps ) ) & + WRITE( stdout,'(5x,"SCF correction compared to forces is large: ", & + & "reduce conv_thr to get better values")') + RETURN + ! +9035 FORMAT(5X,'atom ',I4,' type ',I2,' force = ',3F14.8) + ! +END SUBROUTINE forces diff --git a/patches/qespresso-7.2.diff/PW/src/plugin_ext_forces.f90 b/patches/qespresso-7.2.diff/PW/src/plugin_ext_forces.f90 new file mode 100644 index 0000000000..5d4762fb42 --- /dev/null +++ b/patches/qespresso-7.2.diff/PW/src/plugin_ext_forces.f90 @@ -0,0 +1,74 @@ +! +! Copyright (C) 2001-2009 Quantum ESPRESSO 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 plugin_ext_forces() + !---------------------------------------------------------------------------- + ! + ! + USE mp, ONLY : mp_bcast + USE mp_images, ONLY : intra_image_comm + USE io_global, ONLY : stdout, ionode, ionode_id + USE kinds, ONLY : DP + ! + USE plugin_flags, ONLY : use_plumed + ! + USE cell_base, ONLY : alat, at + USE ions_base, ONLY : tau, nat, amass, ityp + USE force_mod, ONLY : force,sigma + USE control_flags, ONLY : istep + USE ener, ONLY : etot + + USE plumed_module, ONLY: plumed_f_gcmd + ! + IMPLICIT NONE + ! + INTEGER:: i,j,ia + REAL(DP) :: at_plumed(3,3) + REAL(DP) :: virial(3,3) + REAL(DP) :: volume + REAL(DP), ALLOCATABLE :: tau_plumed(:,:) + REAL(DP) :: masses_plumed(nat) + ! + masses_plumed = 0.0_DP + IF(use_plumed) then + IF(ionode)THEN + at_plumed=alat*at; ! the cell, rescaled properly + allocate(tau_plumed(3,nat)) + tau_plumed=alat*tau + volume=+at_plumed(1,1)*at_plumed(2,2)*at_plumed(3,3) & + +at_plumed(1,2)*at_plumed(2,3)*at_plumed(3,1) & + +at_plumed(1,3)*at_plumed(2,1)*at_plumed(3,2) & + -at_plumed(1,1)*at_plumed(3,2)*at_plumed(2,3) & + -at_plumed(1,2)*at_plumed(3,3)*at_plumed(2,1) & + -at_plumed(1,3)*at_plumed(3,1)*at_plumed(2,2) + virial=-sigma*volume + + ! the masses in QE are stored per type, see q-e//Modules/ions_base.f90 + do ia=1,nat + masses_plumed(ia)=amass(ityp(ia)) + end do + + CALL plumed_f_gcmd("setStep"//char(0),istep) + CALL plumed_f_gcmd("setMasses"//char(0),masses_plumed) + CALL plumed_f_gcmd("setForces"//char(0),force) + CALL plumed_f_gcmd("setPositions"//char(0),tau_plumed) + CALL plumed_f_gcmd("setBox"//char(0),at_plumed) + CALL plumed_f_gcmd("setVirial"//char(0),virial) + CALL plumed_f_gcmd("setEnergy"//char(0),etot) + CALL plumed_f_gcmd("calc"//char(0),0) + + sigma=-virial/volume + + deallocate(tau_plumed) + ENDIF + CALL mp_bcast(force, ionode_id, intra_image_comm) + CALL mp_bcast(sigma, ionode_id, intra_image_comm) + ENDIF + ! + ! +END SUBROUTINE plugin_ext_forces diff --git a/patches/qespresso-7.2.diff/PW/src/plugin_ext_forces.f90.preplumed b/patches/qespresso-7.2.diff/PW/src/plugin_ext_forces.f90.preplumed new file mode 100644 index 0000000000..b184498b58 --- /dev/null +++ b/patches/qespresso-7.2.diff/PW/src/plugin_ext_forces.f90.preplumed @@ -0,0 +1,23 @@ +! +! Copyright (C) 2001-2009 Quantum ESPRESSO 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 plugin_ext_forces() + !---------------------------------------------------------------------------- + ! + ! + USE mp, ONLY : mp_bcast + USE mp_images, ONLY : intra_image_comm + USE io_global, ONLY : stdout, ionode, ionode_id + USE kinds, ONLY : DP + ! + USE plugin_flags + ! + IMPLICIT NONE + ! + ! +END SUBROUTINE plugin_ext_forces diff --git a/patches/qespresso-7.2.diff/PW/src/plugin_initialization.f90 b/patches/qespresso-7.2.diff/PW/src/plugin_initialization.f90 new file mode 100644 index 0000000000..0b41434353 --- /dev/null +++ b/patches/qespresso-7.2.diff/PW/src/plugin_initialization.f90 @@ -0,0 +1,64 @@ +! +! Copyright (C) 2010 Quantum ESPRESSO 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 plugin_initialization() + !---------------------------------------------------------------------------- + ! + USE io_global, ONLY : stdout, ionode + USE kinds, ONLY : DP + USE io_files, ONLY : tmp_dir + ! + USE plugin_flags, ONLY : use_plumed + ! + USE ions_base, ONLY : amass, ityp, nat + ! + USE dynamics_module, ONLY : dt + USE constants, ONLY : au_ps + + USE plumed_module, ONLY : plumed_f_installed, plumed_f_gcreate, plumed_f_gcmd + ! + ! + IMPLICIT NONE + ! + INTEGER :: na + INTEGER :: plumedavailable + REAL*8 :: energyUnits,lengthUnits,timeUnits + ! + IF(use_plumed) then + + CALL plumed_f_installed(plumedavailable) + + IF(plumedavailable<=0)THEN + write(stdout,*)"YOU ARE LOOKING FOR PLUMED BUT LOOKS LIKE IT IS NOT AVAILABLE: DO YOU HAVE IT IN YOUR LD_LIBRARY_PATH?" + STOP + ELSE + IF (ionode) THEN + + write(stdout,*)" CREATING PLUMED FROM THE PROGRAM" + call plumed_f_gcreate() + CALL plumed_f_gcmd("setRealPrecision"//char(0),8) + energyUnits=1312.75 ! Ry to kjoule mol + lengthUnits=0.0529177249 ! bohr to nm + timeUnits=2*au_ps ! internal time to ps + call plumed_f_gcmd("setMDEnergyUnits"//char(0),energyUnits) + call plumed_f_gcmd("setMDLengthUnits"//char(0),lengthUnits) + call plumed_f_gcmd("setMDTimeUnits"//char(0),timeUnits) + call plumed_f_gcmd("setPlumedDat"//char(0),"plumed.dat"//char(0)) + call plumed_f_gcmd("setLogFile"//char(0),"PLUMED.OUT"//char(0)) + call plumed_f_gcmd("setNatoms"//char(0),nat) + call plumed_f_gcmd("setMDEngine"//char(0),"qespresso"//char(0)); + call plumed_f_gcmd("setTimestep"//char(0),dt); + call plumed_f_gcmd("init"//char(0),0); + + + ENDIF + ENDIF + ENDIF + ! + ! +END SUBROUTINE plugin_initialization diff --git a/patches/qespresso-7.2.diff/PW/src/plugin_initialization.f90.preplumed b/patches/qespresso-7.2.diff/PW/src/plugin_initialization.f90.preplumed new file mode 100644 index 0000000000..f2c9c0097f --- /dev/null +++ b/patches/qespresso-7.2.diff/PW/src/plugin_initialization.f90.preplumed @@ -0,0 +1,21 @@ +! +! Copyright (C) 2010 Quantum ESPRESSO 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 plugin_initialization() + !---------------------------------------------------------------------------- + ! + USE io_global, ONLY : stdout, ionode + USE kinds, ONLY : DP + USE io_files, ONLY : tmp_dir + ! + USE plugin_flags + ! + IMPLICIT NONE + ! + ! +END SUBROUTINE plugin_initialization diff --git a/patches/qespresso-7.2.diff/PW/src/run_pwscf.f90 b/patches/qespresso-7.2.diff/PW/src/run_pwscf.f90 new file mode 100644 index 0000000000..b94431e9b1 --- /dev/null +++ b/patches/qespresso-7.2.diff/PW/src/run_pwscf.f90 @@ -0,0 +1,569 @@ +! +! Copyright (C) 2013-2020 Quantum ESPRESSO 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 run_pwscf( exit_status ) + !---------------------------------------------------------------------------- + !! Author: Paolo Giannozzi + !! License: GNU + !! Summary: Run an instance of the Plane Wave Self-Consistent Field code + ! + !! Run an instance of the Plane Wave Self-Consistent Field code + !! MPI initialization and input data reading is performed in the + !! calling code - returns in exit_status the exit code for pw.x, + !! returned in the shell. Values are: + !! * 0: completed successfully + !! * 1: an error has occurred (value returned by the errore() routine) + !! * 2-127: convergence error + !! * 2: scf convergence error + !! * 3: ion convergence error + !! * 128-255: code exited due to specific trigger + !! * 255: exit due to user request, or signal trapped, + !! or time > max_seconds + !! (note: in the future, check_stop_now could also return a value + !! to specify the reason of exiting, and the value could be used + !! to return a different value for different reasons) + ! + !! @Note + !! 10/01/17 Samuel Ponce: Add Ford documentation + !! @endnote + !! + ! + USE kinds, ONLY : DP + USE mp, ONLY : mp_bcast, mp_sum + USE io_global, ONLY : stdout, ionode, ionode_id + USE parameters, ONLY : ntypx, npk + USE upf_params, ONLY : lmaxx + USE cell_base, ONLY : fix_volume, fix_area + USE control_flags, ONLY : conv_elec, gamma_only, ethr, lscf, treinit_gvecs + USE control_flags, ONLY : conv_ions, istep, nstep, restart, lmd, lbfgs,& + lensemb, lforce=>tprnfor, tstress + USE control_flags, ONLY : io_level + USE cellmd, ONLY : lmovecell + USE command_line_options, ONLY : command_line + USE force_mod, ONLY : sigma, force + USE check_stop, ONLY : check_stop_init, check_stop_now + USE mp_images, ONLY : intra_image_comm + USE extrapolation, ONLY : update_file, update_pot + USE scf, ONLY : rho + USE lsda_mod, ONLY : nspin + USE fft_base, ONLY : dfftp + ! + USE plugin_flags, ONLY : use_plumed + ! + USE qmmm, ONLY : qmmm_initialization, qmmm_shutdown, & + qmmm_update_positions, qmmm_update_forces + USE qexsd_module, ONLY : qexsd_set_status + USE xc_lib, ONLY : xclib_dft_is, stop_exx, exx_is_active + USE beef, ONLY : beef_energies + USE ldaU, ONLY : lda_plus_u + USE add_dmft_occ, ONLY : dmft + USE extffield, ONLY : init_extffield, close_extffield + USE input_parameters, ONLY : nextffield + ! + USE device_fbuff_m, ONLY : dev_buf + ! +#if defined (__ENVIRON) + USE plugin_flags, ONLY : use_environ + USE environ_pw_module, ONLY : is_ms_gcs, init_ms_gcs +#endif +#if defined (__OSCDFT) + USE plugin_flags, ONLY : use_oscdft + USE oscdft_base, ONLY : oscdft_ctx + USE oscdft_functions, ONLY : oscdft_run_pwscf +#endif + ! + IMPLICIT NONE + ! + INTEGER, INTENT(OUT) :: exit_status + !! Gives the exit status at the end + ! + LOGICAL, EXTERNAL :: matches + ! checks if first string is contained in the second + ! + ! ... local variables + ! + INTEGER :: idone + ! counter of electronic + ionic steps done in this run + INTEGER :: ions_status + ! ions_status = 3 not yet converged + ! ions_status = 2 converged, restart with nonzero magnetization + ! ions_status = 1 converged, final step with current cell needed + ! ions_status = 0 converged, exiting + ! + LOGICAL :: optimizer_failed = .FALSE. + ! + INTEGER :: ierr + ! collect error codes + ! + ions_status = 3 + exit_status = 0 + IF ( ionode ) WRITE( UNIT = stdout, FMT = 9010 ) ntypx, npk, lmaxx + ! + IF (ionode) CALL plugin_arguments() + CALL plugin_arguments_bcast( ionode_id, intra_image_comm ) + ! + ! ... needs to come before iosys() so some input flags can be + ! overridden without needing to write PWscf specific code. + ! + CALL qmmm_initialization() + ! + ! ... convert to internal variables + ! + CALL iosys() + ! + ! ... If executable names is "dist.x", compute atomic distances, angles, + ! ... nearest neighbors, write them to file "dist.out", exit + ! + IF ( matches('dist.x',command_line) ) THEN + IF (ionode) CALL run_dist( exit_status ) + RETURN + ENDIF + ! + IF ( gamma_only ) WRITE( UNIT = stdout, & + & FMT = '(/,5X,"gamma-point specific algorithms are used")' ) + ! + ! call to void routine for user defined / plugin patches initializations + ! +#if defined(__LEGACY_PLUGINS) + CALL plugin_initialization() +#endif +#if defined (__ENVIRON) + IF (use_environ) THEN + IF (is_ms_gcs()) CALL init_ms_gcs() + END IF +#endif + ! + CALL check_stop_init() + ! + CALL setup() + ! + CALL qmmm_update_positions() + ! + ! ... dry run: code will stop here if called with exit file present + ! ... useful for a quick and automated way to check input data + ! + IF ( nstep == 0 .OR. check_stop_now() ) THEN + CALL pre_init() + CALL data_structure( gamma_only ) + CALL summary() + CALL memory_report() + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config-init' ) + RETURN + ENDIF + ! + CALL init_run() + ! + ! read external force fields parameters + ! + IF ( nextffield > 0 .AND. ionode) THEN + ! + CALL init_extffield( 'PW', nextffield ) + ! + END IF + ! + IF ( check_stop_now() ) THEN + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config' ) + RETURN + ENDIF + ! + main_loop: DO idone = 1, nstep + ! + ! ... electronic self-consistency or band structure calculation + ! +#if defined (__OSCDFT) + IF (use_oscdft) THEN + CALL oscdft_run_pwscf(oscdft_ctx) + ELSE +#endif + IF ( .NOT. lscf) THEN + CALL non_scf() + ELSE + CALL electrons() + END IF +#if defined (__OSCDFT) + END IF +#endif + ! + ! ... code stopped by user or not converged + ! + IF ( check_stop_now() .OR. .NOT. conv_elec ) THEN + IF ( check_stop_now() ) THEN + exit_status = 255 + ELSE + IF (dmft) THEN + exit_status = 131 + ELSE + exit_status = 2 + ENDIF + ENDIF + CALL qexsd_set_status(exit_status) + IF(exx_is_active()) then + CALL punch( 'all' ) + ELSE + CALL punch( 'config' ) + ENDIF + RETURN + ENDIF + ! + ! ... file in CASINO format written here if required + ! + IF ( lmd ) THEN + CALL pw2casino( istep ) + ELSE + CALL pw2casino( 0 ) + END IF + ! + ! ... ionic section starts here + ! + CALL start_clock( 'ions' ); !write(*,*)' start ions' ; FLUSH(6) + conv_ions = .TRUE. + ! + ! ... force calculation + ! + IF ( lforce ) CALL forces() + ! + ! ... stress calculation + ! + IF ( tstress ) CALL stress( sigma ) + ! + IF ( lmd .OR. lbfgs ) THEN + ! + ! ... add information on this ionic step to xml file + ! + CALL add_qexsd_step( idone ) + ! + IF (fix_volume) CALL impose_deviatoric_stress( sigma ) + IF (fix_area) CALL impose_deviatoric_stress_2d( sigma ) + ! + ! ... save data needed for potential and wavefunction extrapolation + ! + CALL update_file() + ! + ! ... ionic step (for molecular dynamics or optimization) + ! + CALL move_ions ( idone, ions_status, optimizer_failed ) + conv_ions = ( ions_status == 0 ) .OR. & + ( ions_status == 1 .AND. treinit_gvecs ) + ! + IF ( xclib_dft_is('hybrid') ) CALL stop_exx() + ! + ! ... save restart information for the new configuration + ! + IF ( idone <= nstep .AND. .NOT. conv_ions ) THEN + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config-only' ) + END IF + ! + END IF + ! + CALL stop_clock( 'ions' ); !write(*,*)' stop ions' ; FLUSH(6) + ! + ! ... send out forces to MM code in QM/MM run + ! + CALL qmmm_update_forces( force, rho%of_r, nspin, dfftp ) + ! + ! ... exit condition (ionic convergence) is checked here + ! + IF ( conv_ions .OR. optimizer_failed ) EXIT main_loop + ! + ! ... receive new positions from MM code in QM/MM run + ! + CALL qmmm_update_positions() + ! + ! ... terms of the hamiltonian depending upon nuclear positions + ! ... are reinitialized here + ! + IF ( lmd .OR. lbfgs ) THEN + ! + IF ( ions_status == 1 ) THEN + ! + ! ... final scf calculation with G-vectors for final cell + ! + lbfgs=.FALSE.; lmd=.FALSE. + WRITE( UNIT = stdout, FMT=9020 ) + ! + CALL reset_gvectors( ) + ! + ! ... read atomic occupations for DFT+U(+V) + ! + IF ( lda_plus_u ) CALL read_ns() + ! + ELSE IF ( ions_status == 2 ) THEN + ! + ! ... check whether nonzero magnetization is real + ! + CALL reset_magn() + ! + ELSE + ! + IF ( treinit_gvecs ) THEN + ! + ! ... prepare for next step with freshly computed G vectors + ! + IF ( lmovecell) CALL scale_h() + CALL reset_gvectors ( ) + ! + ELSE + ! + ! ... update the wavefunctions, charge density, potential + ! ... update_pot initializes structure factor array as well + ! + CALL update_pot() + ! + ! ... re-initialize atomic position-dependent quantities + ! + CALL hinit1() + ! + END IF + ! + END IF + ! + ENDIF + ! ... Reset convergence threshold of iterative diagonalization for + ! ... the first scf iteration of each ionic step (after the first) + ! + ethr = 1.0D-6 + ! + CALL dev_buf%reinit( ierr ) + IF ( ierr .ne. 0 ) CALL infomsg( 'run_pwscf', 'Cannot reset GPU buffers! Some buffers still locked.' ) + ! + ENDDO main_loop + ! + ! Set correct exit_status + ! + IF ( .NOT. conv_ions .OR. optimizer_failed ) THEN + exit_status = 3 + ELSE + ! All good + exit_status = 0 + END IF + ! + ! ... save final data file + ! + CALL qexsd_set_status( exit_status ) + IF ( lensemb ) CALL beef_energies( ) + IF ( io_level > -2 ) CALL punch( 'all' ) + ! + ! finalize plumed + ! + IF(use_plumed) then + CALL plumed_f_gfinalize() + ENDIF + ! + CALL qmmm_shutdown() + ! + RETURN + ! +9010 FORMAT( /,5X,'Current dimensions of program PWSCF are:', & + & /,5X,'Max number of different atomic species (ntypx) = ',I2,& + & /,5X,'Max number of k-points (npk) = ',I6, & + & /,5X,'Max angular momentum in pseudopotentials (lmaxx) = ',i2) +9020 FORMAT( /,5X,'Final scf calculation at the relaxed structure.', & + & /,5X,'The G-vectors are recalculated for the final unit cell', & + & /,5X,'Results may differ from those at the preceding step.' ) + ! +END SUBROUTINE run_pwscf +! +! +!------------------------------------------------------------- +SUBROUTINE reset_gvectors( ) +!------------------------------------------------------------- + ! + !! Prepare a new scf calculation with newly recomputed grids, + !! restarting from scratch, not from available data of previous + !! steps (dimensions and file lengths will be different in general) + !! Useful as a check of variable-cell optimization: + !! once convergence is achieved, compare the final energy with the + !! energy computed with G-vectors and plane waves for the final cell + ! + USE io_global, ONLY : stdout + USE basis, ONLY : starting_wfc, starting_pot + USE fft_base, ONLY : dfftp + USE fft_base, ONLY : dffts + USE xc_lib, ONLY : xclib_dft_is + ! + IMPLICIT NONE + ! + ! ... get magnetic moments from previous run before charge is deleted + ! + CALL reset_starting_magnetization() + ! + ! ... clean everything (FIXME: clean only what has to be cleaned) + ! + CALL clean_pw( .FALSE. ) + CALL close_files(.TRUE.) + ! + IF (TRIM(starting_wfc) == 'file') starting_wfc = 'atomic+random' + starting_pot='atomic' + ! + ! ... re-set FFT grids and re-compute needed stuff (FIXME: which?) + ! + dfftp%nr1=0; dfftp%nr2=0; dfftp%nr3=0 + dffts%nr1=0; dffts%nr2=0; dffts%nr3=0 + ! + CALL init_run() + ! + ! ... re-set and re-initialize EXX-related stuff + ! + IF ( xclib_dft_is('hybrid') ) CALL reset_exx( ) + ! +END SUBROUTINE reset_gvectors +! +! +!------------------------------------------------------------- +SUBROUTINE reset_exx( ) +!------------------------------------------------------------- + USE fft_types, ONLY : fft_type_deallocate + USE exx_base, ONLY : exx_grid_init, exx_mp_init, exx_div_check, & + coulomb_fac, coulomb_done + USE exx, ONLY : dfftt, exx_fft_create, deallocate_exx + USE exx_band, ONLY : igk_exx + ! + IMPLICIT NONE + ! + ! ... re-set EXX-related stuff... + ! + IF (ALLOCATED(coulomb_fac) ) DEALLOCATE( coulomb_fac, coulomb_done ) + CALL deallocate_exx( ) + IF (ALLOCATED(igk_exx)) DEALLOCATE(igk_exx) + dfftt%nr1=0; dfftt%nr2=0; dfftt%nr3=0 + CALL fft_type_deallocate( dfftt ) ! FIXME: is this needed? + ! + ! ... re-compute needed EXX-related stuff + ! + CALL exx_grid_init( REINIT = .TRUE. ) + CALL exx_mp_init() + CALL exx_fft_create() + CALL exx_div_check() + ! +END SUBROUTINE reset_exx +! +! +!---------------------------------------------------------------- +SUBROUTINE reset_magn() + !---------------------------------------------------------------- + !! LSDA optimization: a final configuration with zero + !! absolute magnetization has been found and we check + !! if it is really the minimum energy structure by + !! performing a new scf iteration without any "electronic" history. + ! + USE io_global, ONLY : stdout + USE dfunct, ONLY : newd + ! + IMPLICIT NONE + ! + WRITE( UNIT = stdout, FMT = 9010 ) + WRITE( UNIT = stdout, FMT = 9020 ) + ! + ! ... re-initialize the potential (no need to re-initialize wavefunctions) + ! + CALL potinit() + CALL newd() + ! +9010 FORMAT( /5X,'lsda relaxation : a final configuration with zero', & + & /5X,' absolute magnetization has been found' ) +9020 FORMAT( /5X,'the program is checking if it is really ', & + & 'the minimum energy structure', & + & /5X,'by performing a new scf iteration ', & + & 'without any "electronic" history' ) + ! +END SUBROUTINE reset_magn +! +! +!------------------------------------------------------------------- +SUBROUTINE reset_starting_magnetization() + !------------------------------------------------------------------- + !! On input, the scf charge density is needed. + !! On output, new values for starting_magnetization, angle1, angle2 + !! estimated from atomic magnetic moments - to be used in last step. + ! + USE kinds, ONLY : DP + USE constants, ONLY : pi + USE ions_base, ONLY : nsp, ityp, nat + USE lsda_mod, ONLY : nspin, starting_magnetization + USE scf, ONLY : rho + USE noncollin_module, ONLY : noncolin, angle1, angle2, domag + ! + IMPLICIT NONE + ! + ! ... local variables + ! + INTEGER :: i, nt, iat + ! loop counter on species + ! number of atoms per species + ! loop counter on atoms + REAL(DP) :: norm_tot, norm_xy + ! modulus of atomic magnetization + ! xy-projection of atomic magnetization + REAL(DP) :: theta, phi + ! angle between magnetization and z-axis + ! angle between xy-magnetization and x-axis + REAL(DP), ALLOCATABLE :: r_loc(:) + ! auxiliary array for density + REAL(DP), ALLOCATABLE :: m_loc(:,:) + ! auxiliary array for magnetization + ! + IF ( (noncolin .AND. domag) .OR. nspin==2) THEN + ALLOCATE( r_loc(nat), m_loc(nspin-1,nat) ) + CALL get_locals( r_loc,m_loc, rho%of_r ) + ELSE + RETURN + ENDIF + ! + DO i = 1, nsp + ! + starting_magnetization(i) = 0.0_DP + angle1(i) = 0.0_DP + angle2(i) = 0.0_DP + nt = 0 + ! + DO iat = 1, nat + IF (ityp(iat) == i) THEN + nt = nt + 1 + IF (noncolin) THEN + norm_tot = SQRT(m_loc(1,iat)**2+m_loc(2,iat)**2+m_loc(3,iat)**2) + norm_xy = SQRT(m_loc(1,iat)**2+m_loc(2,iat)**2) + IF (norm_tot > 1.d-10) THEN + theta = ACOS(m_loc(3,iat)/norm_tot) + IF (norm_xy > 1.d-10) THEN + phi = ACOS(m_loc(1,iat)/norm_xy) + IF (m_loc(2,iat) < 0.d0) phi = - phi + ELSE + phi = 2.d0*pi + ENDIF + ELSE + theta = 2.d0*pi + phi = 2.d0*pi + ENDIF + angle1(i) = angle1(i) + theta + angle2(i) = angle2(i) + phi + starting_magnetization(i) = starting_magnetization(i) + & + norm_tot/r_loc(iat) + ELSE + starting_magnetization(i) = starting_magnetization(i) + & + m_loc(1,iat)/r_loc(iat) + ENDIF + ENDIF + ENDDO + ! + IF ( nt > 0 ) THEN + starting_magnetization(i) = starting_magnetization(i) / DBLE(nt) + angle1(i) = angle1(i) / DBLE(nt) + angle2(i) = angle2(i) / DBLE(nt) + ENDIF + ! + ENDDO + ! + DEALLOCATE( r_loc, m_loc ) + ! +END SUBROUTINE reset_starting_magnetization diff --git a/patches/qespresso-7.2.diff/PW/src/run_pwscf.f90.preplumed b/patches/qespresso-7.2.diff/PW/src/run_pwscf.f90.preplumed new file mode 100644 index 0000000000..7bc03cbbc5 --- /dev/null +++ b/patches/qespresso-7.2.diff/PW/src/run_pwscf.f90.preplumed @@ -0,0 +1,560 @@ +! +! Copyright (C) 2013-2020 Quantum ESPRESSO 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 run_pwscf( exit_status ) + !---------------------------------------------------------------------------- + !! Author: Paolo Giannozzi + !! License: GNU + !! Summary: Run an instance of the Plane Wave Self-Consistent Field code + ! + !! Run an instance of the Plane Wave Self-Consistent Field code + !! MPI initialization and input data reading is performed in the + !! calling code - returns in exit_status the exit code for pw.x, + !! returned in the shell. Values are: + !! * 0: completed successfully + !! * 1: an error has occurred (value returned by the errore() routine) + !! * 2-127: convergence error + !! * 2: scf convergence error + !! * 3: ion convergence error + !! * 128-255: code exited due to specific trigger + !! * 255: exit due to user request, or signal trapped, + !! or time > max_seconds + !! (note: in the future, check_stop_now could also return a value + !! to specify the reason of exiting, and the value could be used + !! to return a different value for different reasons) + ! + !! @Note + !! 10/01/17 Samuel Ponce: Add Ford documentation + !! @endnote + !! + ! + USE kinds, ONLY : DP + USE mp, ONLY : mp_bcast, mp_sum + USE io_global, ONLY : stdout, ionode, ionode_id + USE parameters, ONLY : ntypx, npk + USE upf_params, ONLY : lmaxx + USE cell_base, ONLY : fix_volume, fix_area + USE control_flags, ONLY : conv_elec, gamma_only, ethr, lscf, treinit_gvecs + USE control_flags, ONLY : conv_ions, istep, nstep, restart, lmd, lbfgs,& + lensemb, lforce=>tprnfor, tstress + USE control_flags, ONLY : io_level + USE cellmd, ONLY : lmovecell + USE command_line_options, ONLY : command_line + USE force_mod, ONLY : sigma, force + USE check_stop, ONLY : check_stop_init, check_stop_now + USE mp_images, ONLY : intra_image_comm + USE extrapolation, ONLY : update_file, update_pot + USE scf, ONLY : rho + USE lsda_mod, ONLY : nspin + USE fft_base, ONLY : dfftp + USE qmmm, ONLY : qmmm_initialization, qmmm_shutdown, & + qmmm_update_positions, qmmm_update_forces + USE qexsd_module, ONLY : qexsd_set_status + USE xc_lib, ONLY : xclib_dft_is, stop_exx, exx_is_active + USE beef, ONLY : beef_energies + USE ldaU, ONLY : lda_plus_u + USE add_dmft_occ, ONLY : dmft + USE extffield, ONLY : init_extffield, close_extffield + USE input_parameters, ONLY : nextffield + ! + USE device_fbuff_m, ONLY : dev_buf + ! +#if defined (__ENVIRON) + USE plugin_flags, ONLY : use_environ + USE environ_pw_module, ONLY : is_ms_gcs, init_ms_gcs +#endif +#if defined (__OSCDFT) + USE plugin_flags, ONLY : use_oscdft + USE oscdft_base, ONLY : oscdft_ctx + USE oscdft_functions, ONLY : oscdft_run_pwscf +#endif + ! + IMPLICIT NONE + ! + INTEGER, INTENT(OUT) :: exit_status + !! Gives the exit status at the end + ! + LOGICAL, EXTERNAL :: matches + ! checks if first string is contained in the second + ! + ! ... local variables + ! + INTEGER :: idone + ! counter of electronic + ionic steps done in this run + INTEGER :: ions_status + ! ions_status = 3 not yet converged + ! ions_status = 2 converged, restart with nonzero magnetization + ! ions_status = 1 converged, final step with current cell needed + ! ions_status = 0 converged, exiting + ! + LOGICAL :: optimizer_failed = .FALSE. + ! + INTEGER :: ierr + ! collect error codes + ! + ions_status = 3 + exit_status = 0 + IF ( ionode ) WRITE( UNIT = stdout, FMT = 9010 ) ntypx, npk, lmaxx + ! + IF (ionode) CALL plugin_arguments() + CALL plugin_arguments_bcast( ionode_id, intra_image_comm ) + ! + ! ... needs to come before iosys() so some input flags can be + ! overridden without needing to write PWscf specific code. + ! + CALL qmmm_initialization() + ! + ! ... convert to internal variables + ! + CALL iosys() + ! + ! ... If executable names is "dist.x", compute atomic distances, angles, + ! ... nearest neighbors, write them to file "dist.out", exit + ! + IF ( matches('dist.x',command_line) ) THEN + IF (ionode) CALL run_dist( exit_status ) + RETURN + ENDIF + ! + IF ( gamma_only ) WRITE( UNIT = stdout, & + & FMT = '(/,5X,"gamma-point specific algorithms are used")' ) + ! + ! call to void routine for user defined / plugin patches initializations + ! +#if defined(__LEGACY_PLUGINS) + CALL plugin_initialization() +#endif +#if defined (__ENVIRON) + IF (use_environ) THEN + IF (is_ms_gcs()) CALL init_ms_gcs() + END IF +#endif + ! + CALL check_stop_init() + ! + CALL setup() + ! + CALL qmmm_update_positions() + ! + ! ... dry run: code will stop here if called with exit file present + ! ... useful for a quick and automated way to check input data + ! + IF ( nstep == 0 .OR. check_stop_now() ) THEN + CALL pre_init() + CALL data_structure( gamma_only ) + CALL summary() + CALL memory_report() + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config-init' ) + RETURN + ENDIF + ! + CALL init_run() + ! + ! read external force fields parameters + ! + IF ( nextffield > 0 .AND. ionode) THEN + ! + CALL init_extffield( 'PW', nextffield ) + ! + END IF + ! + IF ( check_stop_now() ) THEN + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config' ) + RETURN + ENDIF + ! + main_loop: DO idone = 1, nstep + ! + ! ... electronic self-consistency or band structure calculation + ! +#if defined (__OSCDFT) + IF (use_oscdft) THEN + CALL oscdft_run_pwscf(oscdft_ctx) + ELSE +#endif + IF ( .NOT. lscf) THEN + CALL non_scf() + ELSE + CALL electrons() + END IF +#if defined (__OSCDFT) + END IF +#endif + ! + ! ... code stopped by user or not converged + ! + IF ( check_stop_now() .OR. .NOT. conv_elec ) THEN + IF ( check_stop_now() ) THEN + exit_status = 255 + ELSE + IF (dmft) THEN + exit_status = 131 + ELSE + exit_status = 2 + ENDIF + ENDIF + CALL qexsd_set_status(exit_status) + IF(exx_is_active()) then + CALL punch( 'all' ) + ELSE + CALL punch( 'config' ) + ENDIF + RETURN + ENDIF + ! + ! ... file in CASINO format written here if required + ! + IF ( lmd ) THEN + CALL pw2casino( istep ) + ELSE + CALL pw2casino( 0 ) + END IF + ! + ! ... ionic section starts here + ! + CALL start_clock( 'ions' ); !write(*,*)' start ions' ; FLUSH(6) + conv_ions = .TRUE. + ! + ! ... force calculation + ! + IF ( lforce ) CALL forces() + ! + ! ... stress calculation + ! + IF ( tstress ) CALL stress( sigma ) + ! + IF ( lmd .OR. lbfgs ) THEN + ! + ! ... add information on this ionic step to xml file + ! + CALL add_qexsd_step( idone ) + ! + IF (fix_volume) CALL impose_deviatoric_stress( sigma ) + IF (fix_area) CALL impose_deviatoric_stress_2d( sigma ) + ! + ! ... save data needed for potential and wavefunction extrapolation + ! + CALL update_file() + ! + ! ... ionic step (for molecular dynamics or optimization) + ! + CALL move_ions ( idone, ions_status, optimizer_failed ) + conv_ions = ( ions_status == 0 ) .OR. & + ( ions_status == 1 .AND. treinit_gvecs ) + ! + IF ( xclib_dft_is('hybrid') ) CALL stop_exx() + ! + ! ... save restart information for the new configuration + ! + IF ( idone <= nstep .AND. .NOT. conv_ions ) THEN + exit_status = 255 + CALL qexsd_set_status( exit_status ) + CALL punch( 'config-only' ) + END IF + ! + END IF + ! + CALL stop_clock( 'ions' ); !write(*,*)' stop ions' ; FLUSH(6) + ! + ! ... send out forces to MM code in QM/MM run + ! + CALL qmmm_update_forces( force, rho%of_r, nspin, dfftp ) + ! + ! ... exit condition (ionic convergence) is checked here + ! + IF ( conv_ions .OR. optimizer_failed ) EXIT main_loop + ! + ! ... receive new positions from MM code in QM/MM run + ! + CALL qmmm_update_positions() + ! + ! ... terms of the hamiltonian depending upon nuclear positions + ! ... are reinitialized here + ! + IF ( lmd .OR. lbfgs ) THEN + ! + IF ( ions_status == 1 ) THEN + ! + ! ... final scf calculation with G-vectors for final cell + ! + lbfgs=.FALSE.; lmd=.FALSE. + WRITE( UNIT = stdout, FMT=9020 ) + ! + CALL reset_gvectors( ) + ! + ! ... read atomic occupations for DFT+U(+V) + ! + IF ( lda_plus_u ) CALL read_ns() + ! + ELSE IF ( ions_status == 2 ) THEN + ! + ! ... check whether nonzero magnetization is real + ! + CALL reset_magn() + ! + ELSE + ! + IF ( treinit_gvecs ) THEN + ! + ! ... prepare for next step with freshly computed G vectors + ! + IF ( lmovecell) CALL scale_h() + CALL reset_gvectors ( ) + ! + ELSE + ! + ! ... update the wavefunctions, charge density, potential + ! ... update_pot initializes structure factor array as well + ! + CALL update_pot() + ! + ! ... re-initialize atomic position-dependent quantities + ! + CALL hinit1() + ! + END IF + ! + END IF + ! + ENDIF + ! ... Reset convergence threshold of iterative diagonalization for + ! ... the first scf iteration of each ionic step (after the first) + ! + ethr = 1.0D-6 + ! + CALL dev_buf%reinit( ierr ) + IF ( ierr .ne. 0 ) CALL infomsg( 'run_pwscf', 'Cannot reset GPU buffers! Some buffers still locked.' ) + ! + ENDDO main_loop + ! + ! Set correct exit_status + ! + IF ( .NOT. conv_ions .OR. optimizer_failed ) THEN + exit_status = 3 + ELSE + ! All good + exit_status = 0 + END IF + ! + ! ... save final data file + ! + CALL qexsd_set_status( exit_status ) + IF ( lensemb ) CALL beef_energies( ) + IF ( io_level > -2 ) CALL punch( 'all' ) + ! + CALL qmmm_shutdown() + ! + RETURN + ! +9010 FORMAT( /,5X,'Current dimensions of program PWSCF are:', & + & /,5X,'Max number of different atomic species (ntypx) = ',I2,& + & /,5X,'Max number of k-points (npk) = ',I6, & + & /,5X,'Max angular momentum in pseudopotentials (lmaxx) = ',i2) +9020 FORMAT( /,5X,'Final scf calculation at the relaxed structure.', & + & /,5X,'The G-vectors are recalculated for the final unit cell', & + & /,5X,'Results may differ from those at the preceding step.' ) + ! +END SUBROUTINE run_pwscf +! +! +!------------------------------------------------------------- +SUBROUTINE reset_gvectors( ) +!------------------------------------------------------------- + ! + !! Prepare a new scf calculation with newly recomputed grids, + !! restarting from scratch, not from available data of previous + !! steps (dimensions and file lengths will be different in general) + !! Useful as a check of variable-cell optimization: + !! once convergence is achieved, compare the final energy with the + !! energy computed with G-vectors and plane waves for the final cell + ! + USE io_global, ONLY : stdout + USE basis, ONLY : starting_wfc, starting_pot + USE fft_base, ONLY : dfftp + USE fft_base, ONLY : dffts + USE xc_lib, ONLY : xclib_dft_is + ! + IMPLICIT NONE + ! + ! ... get magnetic moments from previous run before charge is deleted + ! + CALL reset_starting_magnetization() + ! + ! ... clean everything (FIXME: clean only what has to be cleaned) + ! + CALL clean_pw( .FALSE. ) + CALL close_files(.TRUE.) + ! + IF (TRIM(starting_wfc) == 'file') starting_wfc = 'atomic+random' + starting_pot='atomic' + ! + ! ... re-set FFT grids and re-compute needed stuff (FIXME: which?) + ! + dfftp%nr1=0; dfftp%nr2=0; dfftp%nr3=0 + dffts%nr1=0; dffts%nr2=0; dffts%nr3=0 + ! + CALL init_run() + ! + ! ... re-set and re-initialize EXX-related stuff + ! + IF ( xclib_dft_is('hybrid') ) CALL reset_exx( ) + ! +END SUBROUTINE reset_gvectors +! +! +!------------------------------------------------------------- +SUBROUTINE reset_exx( ) +!------------------------------------------------------------- + USE fft_types, ONLY : fft_type_deallocate + USE exx_base, ONLY : exx_grid_init, exx_mp_init, exx_div_check, & + coulomb_fac, coulomb_done + USE exx, ONLY : dfftt, exx_fft_create, deallocate_exx + USE exx_band, ONLY : igk_exx + ! + IMPLICIT NONE + ! + ! ... re-set EXX-related stuff... + ! + IF (ALLOCATED(coulomb_fac) ) DEALLOCATE( coulomb_fac, coulomb_done ) + CALL deallocate_exx( ) + IF (ALLOCATED(igk_exx)) DEALLOCATE(igk_exx) + dfftt%nr1=0; dfftt%nr2=0; dfftt%nr3=0 + CALL fft_type_deallocate( dfftt ) ! FIXME: is this needed? + ! + ! ... re-compute needed EXX-related stuff + ! + CALL exx_grid_init( REINIT = .TRUE. ) + CALL exx_mp_init() + CALL exx_fft_create() + CALL exx_div_check() + ! +END SUBROUTINE reset_exx +! +! +!---------------------------------------------------------------- +SUBROUTINE reset_magn() + !---------------------------------------------------------------- + !! LSDA optimization: a final configuration with zero + !! absolute magnetization has been found and we check + !! if it is really the minimum energy structure by + !! performing a new scf iteration without any "electronic" history. + ! + USE io_global, ONLY : stdout + USE dfunct, ONLY : newd + ! + IMPLICIT NONE + ! + WRITE( UNIT = stdout, FMT = 9010 ) + WRITE( UNIT = stdout, FMT = 9020 ) + ! + ! ... re-initialize the potential (no need to re-initialize wavefunctions) + ! + CALL potinit() + CALL newd() + ! +9010 FORMAT( /5X,'lsda relaxation : a final configuration with zero', & + & /5X,' absolute magnetization has been found' ) +9020 FORMAT( /5X,'the program is checking if it is really ', & + & 'the minimum energy structure', & + & /5X,'by performing a new scf iteration ', & + & 'without any "electronic" history' ) + ! +END SUBROUTINE reset_magn +! +! +!------------------------------------------------------------------- +SUBROUTINE reset_starting_magnetization() + !------------------------------------------------------------------- + !! On input, the scf charge density is needed. + !! On output, new values for starting_magnetization, angle1, angle2 + !! estimated from atomic magnetic moments - to be used in last step. + ! + USE kinds, ONLY : DP + USE constants, ONLY : pi + USE ions_base, ONLY : nsp, ityp, nat + USE lsda_mod, ONLY : nspin, starting_magnetization + USE scf, ONLY : rho + USE noncollin_module, ONLY : noncolin, angle1, angle2, domag + ! + IMPLICIT NONE + ! + ! ... local variables + ! + INTEGER :: i, nt, iat + ! loop counter on species + ! number of atoms per species + ! loop counter on atoms + REAL(DP) :: norm_tot, norm_xy + ! modulus of atomic magnetization + ! xy-projection of atomic magnetization + REAL(DP) :: theta, phi + ! angle between magnetization and z-axis + ! angle between xy-magnetization and x-axis + REAL(DP), ALLOCATABLE :: r_loc(:) + ! auxiliary array for density + REAL(DP), ALLOCATABLE :: m_loc(:,:) + ! auxiliary array for magnetization + ! + IF ( (noncolin .AND. domag) .OR. nspin==2) THEN + ALLOCATE( r_loc(nat), m_loc(nspin-1,nat) ) + CALL get_locals( r_loc,m_loc, rho%of_r ) + ELSE + RETURN + ENDIF + ! + DO i = 1, nsp + ! + starting_magnetization(i) = 0.0_DP + angle1(i) = 0.0_DP + angle2(i) = 0.0_DP + nt = 0 + ! + DO iat = 1, nat + IF (ityp(iat) == i) THEN + nt = nt + 1 + IF (noncolin) THEN + norm_tot = SQRT(m_loc(1,iat)**2+m_loc(2,iat)**2+m_loc(3,iat)**2) + norm_xy = SQRT(m_loc(1,iat)**2+m_loc(2,iat)**2) + IF (norm_tot > 1.d-10) THEN + theta = ACOS(m_loc(3,iat)/norm_tot) + IF (norm_xy > 1.d-10) THEN + phi = ACOS(m_loc(1,iat)/norm_xy) + IF (m_loc(2,iat) < 0.d0) phi = - phi + ELSE + phi = 2.d0*pi + ENDIF + ELSE + theta = 2.d0*pi + phi = 2.d0*pi + ENDIF + angle1(i) = angle1(i) + theta + angle2(i) = angle2(i) + phi + starting_magnetization(i) = starting_magnetization(i) + & + norm_tot/r_loc(iat) + ELSE + starting_magnetization(i) = starting_magnetization(i) + & + m_loc(1,iat)/r_loc(iat) + ENDIF + ENDIF + ENDDO + ! + IF ( nt > 0 ) THEN + starting_magnetization(i) = starting_magnetization(i) / DBLE(nt) + angle1(i) = angle1(i) / DBLE(nt) + angle2(i) = angle2(i) / DBLE(nt) + ENDIF + ! + ENDDO + ! + DEALLOCATE( r_loc, m_loc ) + ! +END SUBROUTINE reset_starting_magnetization