diff --git a/.github/workflows/dockerWF.yml b/.github/workflows/dockerWF.yml index e37f001c4f..0598c2c600 100644 --- a/.github/workflows/dockerWF.yml +++ b/.github/workflows/dockerWF.yml @@ -19,10 +19,10 @@ jobs: strategy: fail-fast: false matrix: - # fedora38 (gcc 13) + # fedora39 (gcc 13) # rockylinux8 variant: - - fedora38 + - fedora39 - rocky8 steps: - uses: actions/checkout@v3 diff --git a/.github/workflows/macWF.yml b/.github/workflows/macWF.yml index 3ba17ad760..0743819898 100644 --- a/.github/workflows/macWF.yml +++ b/.github/workflows/macWF.yml @@ -22,7 +22,7 @@ jobs: variant: [ "" , "+allmodules" ] # see https://github.community/t/how-to-conditionally-include-exclude-items-in-matrix-eg-based-on-branch/16853/6 for possible exclusions env: - PYVERS: "py37 py38 py39 py310" + PYVERS: "py39 py310" steps: - uses: actions/checkout@v3 - uses: actions/cache@v3 diff --git a/docker/Makefile b/docker/Makefile index e3675eff78..20246e9a40 100644 --- a/docker/Makefile +++ b/docker/Makefile @@ -1,11 +1,11 @@ -.PHONY: ubuntu plumed2.tgz clean fedora38 rocky8 +.PHONY: ubuntu plumed2.tgz clean fedora39 rocky8 ubuntu: plumed2.tgz docker build -t plumed . -fedora38: plumed2.tgz - docker build -t plumed -f fedora38 . +fedora39: plumed2.tgz + docker build -t plumed -f fedora39 . rocky8: plumed2.tgz docker build -t plumed -f rocky8 . diff --git a/docker/fedora38 b/docker/fedora39 similarity index 99% rename from docker/fedora38 rename to docker/fedora39 index a52c8ade91..95e77e39dd 100644 --- a/docker/fedora38 +++ b/docker/fedora39 @@ -1,4 +1,4 @@ -FROM fedora:38 +FROM fedora:39 # note: at variance with centos7, here we have to explicitly install gcc # mdtraj 1.9.7 does not compile with python 3.11 unless installed from source diff --git a/patches/qespresso-7.2.config b/patches/qespresso-7.2.config new file mode 100644 index 0000000000..69d5dab5ab --- /dev/null +++ b/patches/qespresso-7.2.config @@ -0,0 +1,28 @@ + + +function plumed_preliminary_test(){ +# check if the README.md contains the word ESPRESSO and if qe has been already configured + grep -q ESPRESSO README.md 1>/dev/null 2>/dev/null && test -f make.inc +} + +function plumed_before_patch(){ +cp make.inc make.inc.plumedbck +PWD=`pwd` +echo "include ${PWD}/Plumed.inc ">make.inc +awk '{if($1=="QELIBS" && $2=="="){sub("=","= $(PLUMED_LOAD)"); print}else{print }}' make.inc.plumedbck >> make.inc +} + +function plumed_after_revert(){ + mv make.inc.plumedbck make.inc +} + +function plumed_patch_info(){ + echo "" + echo "For more information on Quantum Espresso you should visit http://www.quantum-espresso.org" + echo "To apply this patch configure Quantum Espresso by running ./configure first." + echo "The newer CMake installation workflow is not supported yet." + echo "To enable PLUMED on md runs use pw.x -plumed < md.in > md.out." + echo "A fixed PLUMED input file name 'plumed.dat' is used." + echo "This patch was kindly provided by Ralf Meyer, email: meyer.ralf(at)yahoo.com" +} + 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..ed5ae7a26d --- /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"); + 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 diff --git a/regtest/basic/rt-make-threads/main.cpp b/regtest/basic/rt-make-threads/main.cpp index 4e10210ac3..815821f8fe 100644 --- a/regtest/basic/rt-make-threads/main.cpp +++ b/regtest/basic/rt-make-threads/main.cpp @@ -11,23 +11,50 @@ using namespace PLMD; void run(std::ostream & os, const std::string & name,std::function f,unsigned nthreads=4,unsigned nrepeats=10){ + + // vector containing possible error messages from threads + std::vector msgs(nthreads); + + // wrapper function that catch exceptions and store their messages + auto g=[&](int i){ + try { + f(i); + } catch(const std::exception& e) { + char buffer[1024]; + std::sprintf(buffer,"(thread %d)\n",i); + msgs[i]=std::string(buffer)+e.what(); + } + }; + os<<"Test "<0) msg+=msgs[j]; + // one could propagate the exception with plumed_error()<0) os<<"failed with error "< threads; - for(unsigned j=0;j0) msg+=msgs[j]; + if(msg.length()>0) os<<"failed with error "<&2 echo 'WARNING: using a legacy ActionRegister.h include path, please use <<#include "core/ActionRegister.h">>' sed -i.bak 's%^#include ".*/ActionRegister.h"%#include "core/ActionRegister.h"%g' "${tmpfile}" @@ -64,14 +77,14 @@ do exit 1 } - rm -f ${tmpfile} ${tmpfile}.bak ${tmpfile%.cpp} objs="$objs $obj" done -if test "$PLUMED_IS_INSTALLED" = yes ; then - eval "$link_installed" "$PLUMED_MKLIB_LDFLAGS" -o "$lib" "$objs" -else - eval "$link_uninstalled" "$PLUMED_MKLIB_LDFLAGS" -o "$lib" "$objs" +link_command="$link_uninstalled" + +if test "$PLUMED_IS_INSTALLED" = yes; then + link_command="$link_installed" fi +eval "$link_command" "$PLUMED_MKLIB_LDFLAGS" $objs -o "$lib" diff --git a/user-doc/tutorials/a-advanced-methods.txt b/user-doc/tutorials/a-advanced-methods.txt new file mode 100644 index 0000000000..ab5203b583 --- /dev/null +++ b/user-doc/tutorials/a-advanced-methods.txt @@ -0,0 +1,504 @@ +a-advanced-methods.txt/** +\page advanced-methods Advanced Methods in MD 2023: Metadynamics simulations with PLUMED + +\section advanced-methods-aims Aim + +The aim of this tutorial is to train users to perform and analyze metadynamics simulations with PLUMED. +This tutorial has been prepared by Max Bonomi (stealing a lot of material from other tutorials) for +the Advanced Methods in MD 2023 workshop, + held in Copenhagen (Denmark) on December 11-12, 2023. + +\section advanced-methods-objectives Objectives + +Once this tutorial is completed users will be able to: + +- Write the PLUMED input file to perform metadynamics simulations. +- Calculate the free energy as a function of the metadynamics collective variables. +- Unbias metadynamics simulations. +- Estimate the error in the reconstructed free energies using block analysis. +- Assess the convergence of metadynamics simulations. + +\section advanced-methods-resources Resources + +The \tarball{advanced-methods} for this tutorial contains the following files: +- `diala.pdb`: a PDB file for alanine dipeptide in vacuo. +- `topol.tpr`: a GROMACS run (binary) file to perform MD simulations of alanine dipeptide. +- `do_block_fes.py`: a python script to perform error analysis of metadynamics simulations. + +After dowloading the compressed archive to your local machine, you can unpack it using the following command: + +\verbatim +tar xvzf advanced-methods.tar.gz +\endverbatim + +Once unpacked, all the files can be found in the `advanced-methods` directory. To keep things clean, +it is recommended to run each exercise in a separate sub-directory that you can create inside `advanced-methods`. + +\note This tutorial has been tested with PLUMED version 2.8.2 and GROMACS version 2019.6. + +\section advanced-methods-intro Introduction + +PLUMED can be used to compute collective variables (CVs) on a pre-calculated +trajectory. However, PLUMED is most often use to add forces on the CVs during a MD simulation, for example, +in order to accelerate sampling. To this aim, we have implemented a variety of possible biases acting on CVs. +The complete documentation for all the biasing methods available in PLUMED can be found at the \ref Bias page. +In the following we will learn how to use PLUMED to perform and analyze a metadynamics simulation. +Here you can find a brief recap of the metadynamics theory. + +\hidden{Summary of theory} + +In metadynamics, an external history-dependent bias potential is constructed in the space of +a few selected degrees of freedom \f$ \vec{s}({q})\f$, generally called collective variables (CVs) \cite metad. +This potential is built as a sum of Gaussian kernels deposited along the trajectory in the CVs space: + +\f[ +V(\vec{s},t) = \sum_{ k \tau < t} W(k \tau) +\exp\left( +-\sum_{i=1}^{d} \frac{(s_i-s_i({q}(k \tau)))^2}{2\sigma_i^2} +\right). +\f] + +where \f$ \tau \f$ is the Gaussian deposition stride, +\f$ \sigma_i \f$ the width of the Gaussian for the \f$i\f$th CV, and \f$ W(k \tau) \f$ the +height of the Gaussian. The effect of the metadynamics bias potential is to push the system away +from local minima into visiting new regions of the phase space. Furthermore, in the long +time limit, the bias potential converges to minus the free energy as a function of the CVs: + +\f[ +V(\vec{s},t\rightarrow \infty) = -F(\vec{s}) + C. +\f] + +In standard metadynamics, Gaussian kernels of constant height are added for the entire course of a +simulation. As a result, the system is eventually pushed to explore high free-energy regions +and the estimate of the free energy calculated from the bias potential oscillates around +the real value. +In well-tempered metadynamics \cite Barducci:2008, the height of the Gaussian +is decreased with simulation time according to: + +\f[ + W (k \tau ) = W_0 \exp \left( -\frac{V(\vec{s}({q}(k \tau)),k \tau)}{k_B\Delta T} \right ), +\f] + +where \f$ W_0 \f$ is an initial Gaussian height, \f$ \Delta T \f$ an input parameter +with the dimension of a temperature, and \f$ k_B \f$ the Boltzmann constant. +With this rescaling of the Gaussian height, the bias potential smoothly converges in the long time limit, +but it does not fully compensate the underlying free energy: + +\f[ +V(\vec{s},t\rightarrow \infty) = -\frac{\Delta T}{T+\Delta T}F(\vec{s}) + C. +\f] + +where \f$ T \f$ is the temperature of the system. +In the long time limit, the CVs thus sample an ensemble +at a temperature \f$ T+\Delta T \f$ which is higher than the system temperature \f$ T \f$. +The parameter \f$ \Delta T \f$ can be chosen to regulate the extent of free-energy exploration: + \f$ \Delta T = 0\f$ corresponds to standard MD, \f$ \Delta T \rightarrow\infty\f$ to standard +metadynamics. In well-tempered metadynamics literature and in PLUMED, you will often encounter +the term "bias factor" which is the ratio between the temperature of the CVs (\f$ T+\Delta T \f$) +and the system temperature (\f$ T \f$): + +\f[ +\gamma = \frac{T+\Delta T}{T}. +\f] + +The bias factor should thus be carefully chosen in order for the relevant free-energy barriers to be crossed +efficiently in the time scale of the simulation. + +Additional information can be found in the several review papers on metadynamics +\cite gerv-laio09review \cite WCMS:WCMS31 \cite WCMS:WCMS1103 \cite bussi2015free. + +\endhidden + +We will play with a toy system, alanine dipeptide simulated in vacuo using the AMBER99SB-ILDN +force field (see Fig. \ref advanced-methods-ala-fig). +This rather simple molecule is useful to understand data analysis and free-energy methods. +This system is a nice example because it presents two metastable states separated by a high free-energy barrier. +It is conventional use to characterize the two states in terms of Ramachandran dihedral angles, which are denoted with \f$ \phi \f$ +(phi) and \f$ \psi \f$ (psi) in Fig. \ref advanced-methods-transition-fig. + +\anchor advanced-methods-ala-fig +\image html belfast-2-ala.png "The molecule of the day: alanine dipeptide." + +\anchor advanced-methods-transition-fig +\image html belfast-2-transition.png "Two metastable states of alanine dipeptide are characterized by their Ramachandran dihedral angles." + +\section advanced-methods-ex Exercises + +\subsection advanced-methods-ex-1 Exercise 1: My first metadynamics simulation + +In this exercise we will setup and perform a well-tempered metadynamics run using the backbone dihedral \f$ \phi \f$ +as collective variable. During the calculation, we will also monitor the behavior of the other backbone dihedral \f$ \psi \f$. + +Here you can find a sample `plumed.dat` file that you can use as a template. +Whenever you see an highlighted \highlight{FILL} string, this is a string that you must replace. + +\plumedfile +# Activate MOLINFO functionalities +MOLINFO STRUCTURE=diala.pdb + +# Compute the backbone dihedral angle phi, defined by atoms C-N-CA-C +# you might want to use MOLINFO shortcuts +phi: TORSION ATOMS=__FILL__ +# Compute the backbone dihedral angle psi, defined by atoms N-CA-C-N +# here also you might want to use MOLINFO shortcuts +psi: TORSION ATOMS=__FILL__ + +# Activate well-tempered metadynamics in phi +metad: __FILL__ ARG=__FILL__ ... +# Deposit a Gaussian every 500 time steps, with initial height equal to 1.2 kJ/mol + PACE=500 HEIGHT=1.2 +# The bias factor should be wisely chosen + BIASFACTOR=__FILL__ +# Gaussian width (sigma) should be chosen based on CV fluctuation in unbiased run + SIGMA=__FILL__ +# Gaussians will be written to file and also stored on grid + FILE=HILLS GRID_MIN=-pi GRID_MAX=pi +... + +# Print both collective variables on COLVAR file every 10 steps +PRINT ARG=__FILL__ FILE=COLVAR STRIDE=__FILL__ +\endplumedfile + +The syntax for the command \ref METAD is simple. +The directive is followed by a keyword `ARG` followed by the labels of the CVs +on which the metadynamics bias potential will act. +The keyword `PACE` determines the stride of Gaussian deposition in number of time steps, +while the keyword `HEIGHT` specifies the height of the Gaussian. For each CVs, one has +to specify the width of the Gaussian by using the keyword `SIGMA`. Gaussian will be written +to the file indicated by the keyword `FILE`. + +In this example, the bias potential will be stored on a grid, whose boundaries are specified by the keywords `GRID_MIN` and `GRID_MAX`. +Notice that you can provide either the number of bins for every CV (`GRID_BIN`) or +the desired grid spacing (`GRID_SPACING`). In case you provide both PLUMED will use +the most conservative choice (highest number of bins) for each dimension. +In case you do not provide any information about bin size (neither `GRID_BIN` nor `GRID_SPACING`) +and if Gaussian width is fixed, PLUMED will use 1/5 of the Gaussian width as grid spacing. +This default choice should be reasonable for most applications. + +Once your `plumed.dat` file is complete, you can run a 10-ns long metadynamics simulations with the following command: + +\verbatim +> gmx mdrun -s topol.tpr -nsteps 5000000 -plumed plumed.dat +\endverbatim + +During the metadynamics simulation, PLUMED will create two files, named `COLVAR` and `HILLS`. +The `COLVAR` file contains all the information specified by the \ref PRINT command, in this case +the value of the backbone dihedrals \f$ \phi \f$ and \f$ \psi \f$ every 10 steps of simulation. +We can use `gnuplot` to visualize the behavior of the metadynamics CV \f$ \phi \f$ during the simulation: + +\verbatim +gnuplot> p "COLVAR" u 1:2 +\endverbatim + +\anchor advanced-methods-phi-fig +\image html munster-metad-phi.png "Time evolution of the metadynamics CV during the first 2 ns of a metadynamics simulation of alanine dipeptide in vacuum." + +By inspecting Figure \ref advanced-methods-phi-fig, we can see that the system is initialized in one of the two metastable +states of alanine dipeptide. After a while (t=0.1 ns), the system is pushed +by the metadynamics bias potential to visit the other local minimum. As the simulation continues, +the bias potential fills the underlying free-energy landscape, and the system is able to diffuse in the +entire phase space. + +The `HILLS` file contains a list of the Gaussian kernels deposited along the simulation. +If we give a look at the header of this file, we can find relevant information about its content: + +\verbatim +#! FIELDS time phi sigma_phi height biasf +#! SET multivariate false +#! SET min_phi -pi +#! SET max_phi pi +\endverbatim + +The line starting with `FIELDS` tells us what is displayed in the various columns of the `HILLS` file: +the simulation time, the instantaneous value of \f$ \phi \f$, the Gaussian width and height, and the bias factor. +We can use the `HILLS` file to visualize the decrease of the Gaussian height during the simulation, +according to the well-tempered recipe: + +\anchor advanced-methods-phihills-fig +\image html munster-metad-phihills.png "Time evolution of the Gaussian height." + +If we look carefully at the scale of the y-axis, we will notice that in the beginning the value +of the Gaussian height is higher than the initial height specified in the input file, which should be 1.2 kJ/mol. +In fact, this column reports the height of the Gaussian scaled by the pre-factor that +in well-tempered metadynamics relates the bias potential to the free energy. + +\warning The fact that the Gaussian height is decreasing to zero should not be used as a measure of convergence +of your metadynamics simulation! + +\subsection advanced-methods-ex-2 Exercise 2: Estimating the free energy as a function of the metadynamics CVs + +One can estimate the free energy as a function of the metadynamics CVs directly from the metadynamics +bias potential. In order to do so, the utility \ref sum_hills can be used to sum the Gaussian kernels +deposited during the simulation and stored in the `HILLS` file. + +To calculate the free energy as a function of \f$ \phi \f$, it is sufficient to use the following command line: + +\verbatim +plumed sum_hills --hills HILLS +\endverbatim + +The command above generates a file called `fes.dat` in which the free-energy surface as function +of \f$ \phi \f$ is calculated on a regular grid. One can modify the default name for the free-energy file, +as well as the boundaries and bin size of the grid, by using the following \ref sum_hills options: + +\verbatim +--outfile - specify the outputfile for sumhills +--min - the lower bounds for the grid +--max - the upper bounds for the grid +--bin - the number of bins for the grid +--spacing - grid spacing, alternative to the number of bins +\endverbatim + +The result should look like this: + +\anchor advanced-methods-metad-phifes-fig +\image html munster-metad-phifes.png "Estimate of the free energy as a function of the dihedral phi from a 10ns-long well-tempered metadynamics simulation." + +To give a preliminary assessment of the convergence of a metadynamics simulation, one can calculate the estimate of the free energy as a function +of simulation time. At convergence, the reconstructed profiles should be similar. +The \ref sum_hills option `--stride` should be used to give an estimate of the free energy every `N` Gaussian kernels deposited, and +the option `--mintozero` can be used to align the profiles by setting the global minimum to zero. +If we use the following command line: + +\verbatim +plumed sum_hills --hills HILLS --stride 100 --mintozero +\endverbatim + +one free energy is calculated every 100 Gaussian kernels deposited, and the global minimum is set to zero in all profiles. +The resulting plot should look like the following: + +\anchor advanced-methods-metad-phifest-fig +\image html munster-metad-phifest.png "Estimates of the free energy as a function of the dihedral phi calculated every 100 Gaussian kernels deposited." + +These two qualitative observations: +1. the system is diffusing rapidly in the entire CV space (Figure \ref advanced-methods-phi-fig) +2. the estimated free energy does not significantly change as a function of time (Figure \ref advanced-methods-metad-phifest-fig) + +suggest that the simulation __might__ be converged. + +\warning The two conditions listed above are necessary, but not sufficient to declare convergence. +For a quantitative analysis of the convergence of metadynamics simulations, please have a look below at \ref advanced-methods-ex-4. + +\subsection advanced-methods-ex-3 Exercise 3: Reweighting (unbiasing) a metadynamics simulation + +In the previous exercise we biased \f$\phi\f$ and computed the free energy as a function of +the same variable directly from the metadynamics bias potential using the \ref sum_hills utility. +However, in many cases you might decide which variable should be analyzed _after_ +having performed a metadynamics simulation. For example, +you might want to calculate the free energy as a function of CVs other than those +biased during the metadynamics simulation, such as the dihedral \f$ \psi \f$. +At variance with standard MD simulations, you cannot simply calculate histograms of other variables directly from your metadynamics trajectory, +because the presence of the metadynamics bias potential has altered the statistical weight of each frame. +To remove the effect of this bias and thus be able to calculate properties of the system in the unbiased ensemble, +you must reweight (unbias) your simulation. + +There are multiple ways to calculate the correct +statistical weight of each frame in your metadynamics trajectory and thus to reweight your simulation. +For example: + +1. weights can be calculated by considering the time-dependence of the metadynamics bias + potential \cite Tiwary_jp504920s; +2. weights can be calculated using the metadynamics bias potential obtained at the end of the + simulation and assuming a constant bias during the entire course of the simulation \cite Branduardi:2012dl. + +In this exercise we will use the second method, which resembles the umbrella-sampling reweighting approach. +In order to compute the weights we will use the \ref driver tool. + +First of all, you need to prepare a `plumed_reweight.dat` file that is identical to the one you used +for running your metadynamics simulation except for a few modifications. First, you +need to add the keyword `RESTART=YES` to the \ref METAD command. +This will make this action behave as if PLUMED was restarting, i.e. PLUMED will +read from the `HILLS` file the Gaussians that have previously been accumulated. +Second, you need to set the Gaussian `HEIGHT` to zero and the `PACE` to a large number. +This will actually avoid adding new Gaussians (and even if they are added they will have +zero height). Finally, you need to modify the \ref PRINT statement so that you +write every frame and that, in addition to `phi` and `psi`, +you also write `metad.bias`. You might also want to change the name of the output file to `COLVAR_REWEIGHT`. + +\plumedfile +# Activate MOLINFO functionalities +MOLINFO STRUCTURE=diala.pdb + +__FILL__ # here goes the definitions of the phi and psi CVs + +# Activate well-tempered metadynamics in phi +metad: __FILL__ ARG=__FILL__ ... +# Deposit a Gaussian every 10000000 time steps (never!), with initial height equal to 0.0 kJ/mol + PACE=10000000 HEIGHT=0.0 # <- this is the new stuff! +# The bias factor should be wisely chosen + BIASFACTOR=__FILL__ +# Gaussian width (sigma) should be chosen based on CV fluctuation in unbiased run + SIGMA=__FILL__ +# Gaussians will be written to file and also stored on grid + FILE=HILLS GRID_MIN=-pi GRID_MAX=pi +# Say that METAD should be restarting (= reading an existing HILLS file) + RESTART=YES # <- this is the new stuff! +... + +# Print out the values of phi, psi and the metadynamics bias potential +# Make sure you print out the 3 variables in the specified order at every step +PRINT ARG=__FILL__ FILE=COLVAR_REWEIGHT STRIDE=__FILL__ # <- also change this one! +\endplumedfile + +Then run the \ref driver tool using this command: + +\verbatim +> plumed driver --mf_xtc traj_comp.xtc --plumed plumed_reweight.dat --kt 2.494339 +\endverbatim + +Notice that you have to specify the value of \f$k_BT\f$ in energy units. While running your simulation +this information was communicated by the MD code. + +As a result, PLUMED will produce a new `COLVAR_REWEIGHT` file with one additional column containing the +metadynamics bias potential \f$ V(s) \f$ calculated using all the Gaussians deposited along the entire trajectory. +The beginning of the file should look like this: + +\verbatim +#! FIELDS time phi psi metad.bias +#! SET min_phi -pi +#! SET max_phi pi +#! SET min_psi -pi +#! SET max_psi pi + 0.000000 -1.497988 0.273498 110.625670 + 1.000000 -1.449714 0.576594 110.873141 + 2.000000 -1.209587 0.831417 109.742353 + 3.000000 -1.475975 1.279726 110.752327 +\endverbatim + +You can easily obtain the weight \f$ w \f$ of each frame using the expression \f$w\propto\exp\left(\frac{V(s)}{k_BT}\right)\f$ +(umbrella-sampling-like reweighting). At this point, you can read the `COLVAR_REWEIGHT` file using, for example, a python script +and compute a weighted histogram. Alternatively, if you want PLUMED to do the weighted histograms for you, you can add the following +lines at the end of the `plumed_reweight.dat` file: + +\plumedfile +# Use the metadynamics bias as argument +as: REWEIGHT_BIAS ARG=__FILL__ + +# Calculate histograms of phi and psi dihedrals every 50 steps +# using the weights obtained from the metadynamics bias potentials (umbrella-sampling-like reweighting) +# Look at the manual to understand the parameters of the HISTOGRAM action! +hhphi: HISTOGRAM ARG=phi STRIDE=50 GRID_MIN=-pi GRID_MAX=pi GRID_BIN=50 BANDWIDTH=0.05 LOGWEIGHTS=as +hhpsi: HISTOGRAM ARG=psi STRIDE=50 GRID_MIN=-pi GRID_MAX=pi GRID_BIN=50 BANDWIDTH=0.05 LOGWEIGHTS=as + +# Convert histograms h(s) to free energies F(s) = -kBT * log(h(s)) +ffphi: CONVERT_TO_FES GRID=hhphi +ffpsi: CONVERT_TO_FES GRID=hhpsi + +# Print out the free energies F(s) to file once the entire trajectory is processed +DUMPGRID GRID=ffphi FILE=ffphi.dat +DUMPGRID GRID=ffpsi FILE=ffpsi.dat +\endplumedfile + +and plot the result using `gnuplot`: + +\verbatim +gnuplot> p "ffphi.dat" u 1:2 w lp +gnuplot> p "ffpsi.dat" u 1:2 w lp +\endverbatim + +You can now compare the free energies as a function of \f$ \phi \f$ calculated: +1. directly from the metadynamics bias potential using \ref sum_hills as done in \ref advanced-methods-ex-2; +2. using the reweighting procedure introduced in this exercise. + +The results should be identical (see Fig. \ref advanced-methods-fescomp-fig). + +\anchor advanced-methods-fescomp-fig +\image html master-ISDD-2-fescomp-fig.png "Comparison between the free energy as a function of the dihedral phi calculated from the metadynamics bias potential (bias) and by reweighting (rew)". + + +\subsection advanced-methods-ex-4 Exercise 4: Estimating the error in free energies using block-analysis + +In the previous exercise, we calculated the _final_ bias \f$ V(s) \f$ on the entire metadynamics trajectory and we used +this quantity to calculate the correct statistical weight of each frame that we need to reweight the biased simulation. +In this exercise, we will see how this information can be used to calculate the error in +the reconstructed free energies and assess whether our simulation is converged or not. +Let's first calculate the un-biasing weights \f$w\propto\exp\left(\frac{V(s)}{k_BT}\right)\f$ from the `COLVAR_REWEIGHT` +file obtained at the end of \ref advanced-methods-ex-3: + +\verbatim +# Find maximum value of bias to avoid numerical errors when calculating the un-biasing weights +bmax=`awk 'BEGIN{max=0.}{if($1!="#!" && $4>max)max=$4}END{print max}' COLVAR_REWEIGHT` + +# Print phi values and un-biasing (un-normalized) weights +awk '{if($1!="#!") print $2,exp(($4-bmax)/kbt)}' kbt=2.494339 bmax=$bmax COLVAR_REWEIGHT > phi.weight +\endverbatim + +If you inspect the `phi.weight` file, you will see that each line contains +the value of the dihedral \f$ \phi \f$ along with the corresponding (un-normalized) weight \f$ w \f$ for each frame of the +metadynamics trajectory: + +\verbatim +0.907347 0.0400579 +0.814296 0.0169656 +1.118951 0.0651276 +1.040781 0.0714174 +1.218571 0.0344903 +1.090823 0.0700568 +1.130800 0.0622998 +\endverbatim + +At this point we can apply the block-analysis technique (for more info about the theory, +have a look at \ref trieste-2) to calculate the average free energy across the blocks +and the error as a function of block size. For your convenience, you can use the `do_block_fes.py` python +script to read the `phi.weight` file and produce the desired output. +We use a bash loop to test block sizes ranging from 1 to 1000: + +\verbatim +# Arguments of do_block_fes.py +# - input file with CV value and weight for each frame of the trajectory: phi.weight +# - number of CVs: 1 +# - CV range (min, max): (-3.141593, 3.141593) +# - # points in output free energy: 51 +# - kBT (kJoule/mol): 2.494339 +# - Block size: 1<=i<=1000 (every 10) +# +for i in `seq 1 10 1000`; do python3 do_block_fes.py phi.weight 1 -3.141593 3.141593 51 2.494339 $i; done +\endverbatim + +For each value of block size `N`, you will obtain a separate `fes.N.dat` file, containing the value +of the \f$ \phi \f$ variable on a grid, the average free energy across the blocks with its associated error (in kJ/mol) +on each point of the grid: + +\verbatim + -3.141593 23.184653 0.080659 + -3.018393 17.264462 0.055181 + -2.895194 13.360259 0.047751 + -2.771994 10.772696 0.043548 + -2.648794 9.403544 0.042022 +\endverbatim + +Finally, we can calculate the average error along each free-energy profile as a function of the block size: + +\verbatim +for i in `seq 1 10 1000`; do a=`awk '{tot+=$3}END{print tot/NR}' fes.$i.dat`; echo $i $a; done > err.blocks +\endverbatim + +and visualize it using `gnuplot`: + +\verbatim +gnuplot> p "err.blocks" u 1:2 w lp +\endverbatim + +As expected, the error increases with the block size until it reaches a plateau in correspondence of a dimension +of the block that exceeds the correlation between data points (Fig. \ref advanced-methods-block-phi). + +\anchor advanced-methods-block-phi +\image html trieste-4-block-phi.png "Block analysis of a metadynamics simulation using phi as CV" + +__What can we learn from this analysis about the convergence of the metadynamics simulation?__ + +\section advanced-methods-conclusions Conclusions + +In summary, in this tutorial you should have learned how to use PLUMED to: +- Setup and run a metadynamics calculation. +- Compute free energies from the metadynamics bias potential using the \ref sum_hills utility. +- Reweight a metadynamics simulation. +- Calculate errors and assess convergence. + +*/ + +link: @subpage advanced-methods + +description: This tutorial explains how to use PLUMED to run metadynamics simulations + +additional-files: advanced-methods diff --git a/user-doc/tutorials/advanced-methods/diala.pdb b/user-doc/tutorials/advanced-methods/diala.pdb new file mode 100644 index 0000000000..2adf83e019 --- /dev/null +++ b/user-doc/tutorials/advanced-methods/diala.pdb @@ -0,0 +1,24 @@ +ATOM 1 1HH3 ACE 1 15.700 13.120 -0.270 1.00 0.00 +ATOM 2 CH3 ACE 1 14.670 12.910 -0.560 1.00 0.00 +ATOM 3 2HH3 ACE 1 14.080 12.680 0.330 1.00 0.00 +ATOM 4 3HH3 ACE 1 14.550 12.130 -1.310 1.00 0.00 +ATOM 5 C ACE 1 14.140 14.180 -1.210 1.00 0.00 +ATOM 6 O ACE 1 12.990 14.170 -1.650 1.00 0.00 +ATOM 7 N ALA 2 14.940 15.200 -1.540 1.00 0.00 +ATOM 8 H ALA 2 15.870 15.140 -1.150 1.00 0.00 +ATOM 9 CA ALA 2 14.600 16.550 -1.940 1.00 0.00 +ATOM 10 HA ALA 2 13.610 16.860 -1.590 1.00 0.00 +ATOM 11 CB ALA 2 15.650 17.430 -1.260 1.00 0.00 +ATOM 12 HB1 ALA 2 16.670 17.190 -1.550 1.00 0.00 +ATOM 13 HB2 ALA 2 15.420 18.460 -1.520 1.00 0.00 +ATOM 14 HB3 ALA 2 15.490 17.340 -0.190 1.00 0.00 +ATOM 15 C ALA 2 14.330 16.660 -3.430 1.00 0.00 +ATOM 16 O ALA 2 14.330 17.750 -4.010 1.00 0.00 +ATOM 17 N NME 3 14.120 15.510 -4.080 1.00 0.00 +ATOM 18 H NME 3 13.940 14.630 -3.620 1.00 0.00 +ATOM 19 CH3 NME 3 13.800 15.620 -5.490 1.00 0.00 +ATOM 20 1HH3 NME 3 13.140 14.830 -5.860 1.00 0.00 +ATOM 21 2HH3 NME 3 13.250 16.540 -5.710 1.00 0.00 +ATOM 22 3HH3 NME 3 14.690 15.590 -6.120 1.00 0.00 +TER +ENDMDL diff --git a/user-doc/tutorials/advanced-methods/do_block_fes.py b/user-doc/tutorials/advanced-methods/do_block_fes.py new file mode 100644 index 0000000000..c54d808beb --- /dev/null +++ b/user-doc/tutorials/advanced-methods/do_block_fes.py @@ -0,0 +1,185 @@ +import math +import sys +import numpy as np + +# Arguments of do_block_fes.py +# - FILE: input file, 1 column per CV + weights (optional) +# - NCV: number of CVs +# - *MIN: minimum value of CV +# - *MAX: max value of CV +# - *NBIN: # points in output free energy +# - KBT: temperature in energy units (kJoule/mol) +# - N: Block size +# +# * = repeat this block for each CV +# Example with 2 CVs: +# python3 do_block_fes.py phi_psi_w.dat 2 -3.141593 3.141593 50 -3.141593 3.141593 50 2.494339 100 +# +# +# Author: Max Bonomi (mbonomi@pasteur.fr) +# +# +# useful functions +# nD indexes from 1D index +def get_indexes_from_index(index, nbin): + indexes = [] + # get first index + indexes.append(index%nbin[0]) + # loop + kk = index + for i in range(1, len(nbin)-1): + kk = ( kk - indexes[i-1] ) / nbin[i-1] + indexes.append(kk%nbin[i]) + if(len(nbin)>=2): + indexes.append( ( kk - indexes[len(nbin)-2] ) / nbin[len(nbin) -2] ) + return tuple(indexes) + +# nD indexes from values +def get_indexes_from_cvs(cvs, gmin, dx, nbin): + idx = [] + for i in range(0, len(cvs)): + j = int( round( ( cvs[i] - gmin[i] ) / dx[i] ) ) + # check boundaries + if(j>=nbin[i]): + print("Point outside grid, check boundaries!") + exit() + idx.append(j) + return tuple(idx) + +# 1D index from values +def get_index_from_cvs(cvs, gmin, dx, nbin): + # get nD indices from value + idx = get_indexes_from_cvs(cvs, gmin, dx, nbin) + # transform in 1D index + i = idx[-1] + for j in range(len(nbin)-1,0,-1): + i = i*nbin[j-1]+idx[j-1] + return i + +# grid points from nD indexes +def get_points_from_indexes(idx, gmin, dx): + xs = [] + for i in range(0, len(idx)): + xs.append(gmin[i] + float(idx[i]) * dx[i]) + return xs + +# read CV/weight file and create numpy arrays +def read_file(filename,gmin,dx,nbin): + # read file and store lists + cvs=[]; ws=[] + # number of cv + ncv = len(gmin) + for lines in open(filename, "r").readlines(): + riga = lines.strip().split() + # check format + if(len(riga)!=ncv and len(riga)!=ncv+1): + print (filename,"is in the wrong format!") + exit() + # read CVs + cv = [] + for i in range(0, ncv): cv.append(float(riga[i])) + # get index in flattened array + idx = get_index_from_cvs(cv, gmin, dx, nbin) + # read weight, if present + if(len(riga)==ncv+1): + w = float(riga[ncv]) + else: w = 1.0 + # store into cv and weight lists + cvs.append(idx) + ws.append(w) + # return numpy arrays + return np.array(cvs),np.array(ws) + +# 1) READ INPUT parameters +# FILE with CV trajectory and (optionally) weights +FILENAME_ = sys.argv[1] +# number of CVs +NCV_ = int(sys.argv[2]) +# read minimum, maximum and number of bins for FES grid +gmin = []; gmax = []; nbin = [] +for i in range(0, NCV_): + i0 = 3*i + 3 + gmin.append(float(sys.argv[i0])) + gmax.append(float(sys.argv[i0+1])) + nbin.append(int(sys.argv[i0+2])) +# read KBT_ +KBT_ = float(sys.argv[3*NCV_+3]) +# block size +BSIZE_ = int(sys.argv[-1]) + +# 2) SETUP +# define bin sizes +dx = [] +for i in range(0, NCV_): + dx.append( (gmax[i]-gmin[i])/float(nbin[i]-1) ) +# total numbers of bins +nbins = 1 +for i in range(0, len(nbin)): nbins *= nbin[i] +# read file and store arrays +cv, w = read_file(FILENAME_, gmin, dx, nbin) +# total number of data points +ndata = cv.shape[0] +# number of blocks +nblock = int(ndata/BSIZE_) +# prepare numpy arrays for histogram and normalization +histo = np.zeros((nbins,nblock)) +norm = np.zeros(nblock) + +# 3) FILL IN HISTOGRAM ARRAY +for iblock in range(0, nblock): + # define range + i0 = iblock * BSIZE_ + i1 = i0 + BSIZE_ + # cycle on points in the block + for i in range(i0, i1): + # update histogram + histo[cv[i],iblock] += w[i] + # calculate normalization of the block + norm[iblock] = np.sum(w[i0:i1]) + # normalize block + histo[:,iblock] /= norm[iblock] + +# 4) CALCULATE WEIGHTED AVERAGE AND VARIANCE +# now we calculate weighted average across blocks +ave = np.sum(histo*norm, axis=1) / np.sum(norm) +avet = np.transpose(np.tile(ave, (nblock,1))) +# and variance +var = np.sum(np.power( norm * (histo-avet), 2), axis=1) / np.power(np.sum(norm), 2) + +# 5) PRINT FES + ERROR +log = open("fes."+str(BSIZE_)+".dat", "w") +# this is needed to add a blank line +xs_old = [] +for i in range(0, nbins): + # get the indexes in the multi-dimensional grid + idx = get_indexes_from_index(i, nbin) + # get values for grid point + xs = get_points_from_indexes(idx, gmin, dx) + # add a blank line for gnuplot + if(i == 0): + xs_old = xs[:] + else: + flag = 0 + for j in range(1,len(xs)): + if(xs[j] != xs_old[j]): + flag = 1 + xs_old = xs[:] + if (flag == 1): log.write("\n") + # print grid point + for x in xs: + log.write("%12.6lf " % x) + # calculate fes and error + try: + # fes + fes = -KBT_ * math.log(ave[i]) + # variance fes + varf = math.pow( KBT_ / ave[i], 2.0) * var[i] + # error fes + errf = math.sqrt(varf) + # printout + log.write(" %12.6lf %12.6lf\n" % (fes, errf)) + except ValueError: + log.write(" %12s %12s\n" % ("Inf","Inf")) + except OverflowError: + log.write(" %12s %12s\n" % ("Inf","Inf")) +log.close() diff --git a/user-doc/tutorials/advanced-methods/topol.tpr b/user-doc/tutorials/advanced-methods/topol.tpr new file mode 100644 index 0000000000..9bbf71ba8d Binary files /dev/null and b/user-doc/tutorials/advanced-methods/topol.tpr differ