diff --git a/ChangeLog b/ChangeLog new file mode 100644 index 0000000000..e900dd1038 --- /dev/null +++ b/ChangeLog @@ -0,0 +1,219 @@ +2010-09-26 16:03:43 r7 Andrea Marini + + Bug in K_Haydock.F fixed + +2010-09-21 09:43:22 r6 Andrea Marini + + Version 3.2.2 + + * configure include/version.inc int_modules/.objects int_modules/mod_xc2y.F p2y/p2y_db1.F a2y/a2y_db1.F config/Makefile.in config/configure.ac xc_functionals/el_current.F xc_functionals/XC_potential_driver.F xc_functionals/xc_lda_driver.F xc_functionals/XC_lda_driver.F xc_functionals/XC_libxc_driver.F xc_functionals/spin_density.F xc_functionals/.objects qp/QP_driver.F qp/XCo_local.F qp/QP_real_axis.F qp/QP_newton.F qp/bracket.F qp/QP_life_transitions.F qp/QP_secant.F qp/secant.F qp/QP_W2Sc.F qp/.objects common/FREQUENCIES_reset.F common/FREQUENCIES_Green_Function.F common/.objects bz_ops/k_ibz2bz.F pol_function/X_os.F pol_function/O_driver.F pol_function/X_em1.F pol_function/Dipole_driver.F pol_function/Dipole_transverse.F interface/init_load.F interface/init.F interface/barriers.F setup/G_shells_finder.F tddft/tddft_alda_g_space.F io/variables_BS.F io/ioDB1.F io/io_HF_and_locXC.F io/ioX.F io/ioQP.F io/io_header.F io/ioBSS_Haydock.F io/ioOSTNTS.F io/.objects io/io_DIPOLES.F modules/mod_BS.F modules/mod_R_lattice.F modules/mod_fields.F modules/mod_functions.F modules/mod_libxc_int.F modules/mod_drivers.F modules/mod_xc_functionals.F modules/mod_X_output.F modules/std_presets.F modules/mod_X.F modules/mod_IO.F modules/mod_com.F modules/.objects modules/mod_global_XC.F bse/K_Haydock.F bse/K_BSmat_by_V.F bse/K_diagonalization.F bse/K_solvers.F bse/K_output_file.F bse/K.F ChangeLog ypp/ypp_init.F ypp/plot_cube.F ypp/excitons_sort_and_report.F ypp/k_grids.F ypp/.objects ypp/excitons_driver.F ypp/mod_YPP.F ypp/plot_check_and_launch.F ypp/electrons_driver.F ypp/ypp_init_load.F driver/yambo_driver.F driver/driver.c driver/yambo.h driver/ypp.h sbin/make_makefile.sh.in + + Additions: + + - *** BSE/TDDFT beyond the TDA *** + It is now possible to calculate TDDFT and BSE response functions beyond the + Tamm-Dancoff approximation. The antiresonant and coupling blocks of the + Hamiltonians can be computed. + The final Dyson equation is, then. solved using diagonalization or the + quasi-Hermitian Haydock algorithm proposed in Nano Letters, 9, 2820 (2009) + + - *** LIBXC support *** + This means that Yambo can ALWAYS read GS (abinit) databases generated using libxc. + If, in addition, yambo is linked against LIBXC (see --with-libxc-include/--with-libxc-lib options + of configure) this is used to evalute DFT xc terms. + + - *** REAL AXIS GW *** + It is now possible to calculate the GW self-energy on the real-axis. This means no plasmon-pole + approx and the possibility of calculating electronic lifetimes. + In addition the secant (recursive) solver of the Dyson equation has been added. + + - *** CUBE format *** + Added support plot 3D isosurfaces in Gaussian cube format. + + Patch sent by: The Yambo Team + +2010-07-09 10:05:07 r5 Andrea Marini + + Version 3.2.2 + + * configure include/version.inc config/Makefile.in config/configure.ac interface/init_load.F sbin/svn2cl.pl + + Bugs: + - In src/interface/init_load.F a too long description ( > 100 Chars) was causing seg-fault + + Changes: + - Tests option removed from Makefile. To be soon replaced by the use of the test-suite branch. + + Patch sent by: Andrea + +2010-07-02 14:56:07 r4 Andrea Marini + + Version 3.2.2 + + * configure include/version.inc e2y/e2y_wf.F e2y/e2y_db1.F e2y/e2y_i.F e2y/e2y_kb_pp.F e2y/etsf_data.F int_modules/mod_wf2y.F int_modules/mod_com2y.F p2y/qexml_v3.1.F p2y/qexml_v4.0.F p2y/mod_pw_data.F p2y/qexml_v3.2.F p2y/p2y_db1.F p2y/p2y_i.F p2y/pw_errore.F p2y/mod_pw_export.F p2y/p2y_wf.F p2y/mod_p2y.F a2y/a2y_db1.F a2y/.objects a2y/a2y_i.F a2y/hdr_io.F a2y/a2y_wf.F a2y/a2y_KSS_file_name.F config/libxc.m4 config/report.in config/Makefile.in config/acx_report.m4 config/fftw.m4 config/acx.m4 config/acx_cpp.m4 config/acx_mpi.m4 config/netcdf_f90.m4 config/acx_fortran_flags.m4 config/configure.ac config/ax_f90_module_flag.m4 config/setup.in config/acx_get_fc_version.m4 doc/QuickGuidedTour.pdf doc/sample/bulk_silicon/SAVE doc/sample/bulk_silicon/SAVE/database_DB1_NETCDF_format.gz doc/sample/bulk_silicon/SAVE/database_WF_NETCDF_format.gz doc/README lib/lapack/dlacpy.f lib/lapack/slacpy.f lib/lapack/.objects xc_functionals/el_density.F xc_functionals/xcxalp.F xc_functionals/.objects xc_functionals/xc_lda_driver.F xc_functionals/xc_rpa_kp.F xc_functionals/mod_xc_costants.F xc_functionals/xcspol.F xc_functionals/el_magnetization.F qp/QP_driver.F qp/DFT_Vxc.F qp/QP_XX_Vxc.F qp/.objects qp/XCo_local.F qp/QP_of.F qp/XCo_driver.F qp/XCo_Hartree_Fock.F qp/QP_descriptions.F qp/QP_report_and_write.F qp/QP_ppa_chosex.F qp/QP_ppa_cohsex.F qp/QP_newton.F qp_ctl/QP_apply_interpolate.F qp_ctl/QP_apply_stretch.F qp_ctl/mod_QP_CTL.F qp_ctl/QP_apply_global_stretch.F qp_ctl/QP_apply_DB_interpolation.F qp_ctl/.objects qp_ctl/QP_load_DB.F qp_ctl/QP_check_if_corrected.F qp_ctl/QP_fit_DB_values.F qp_ctl/QP_apply_fit.F qp_ctl/QP_apply_db.F qp_ctl/QP_apply_done.F external_c/c_printing.c common/pol_fit.F common/OCCUPATIONS_Extend.F common/extend_occupations.F common/freqs_setup.F common/kk.F common/QP_state_table_setup.F common/FREQUENCIES_damping.F common/OCCUPATIONS_Fermi.F common/FREQUENCIES_coarse_grid.F common/QP_state_extract.F common/Kramers_Kronig.F common/eval_G_minus_G.F common/damping.F common/coarse_grid.F common/.objects common/crystal_lattice.F common/fermi_level.F common/FREQUENCIES_setup.F bz_ops/k_the_nearest.F bz_ops/k_ibz2bz.F bz_ops/k_sym2sym.F bz_ops/k_lattice.F bz_ops/k_reduce.F bz_ops/.objects bz_ops/k_expand.F bz_ops/bz_samp_indexes.F parser/Gclose.F parser/parser_lib.F parser/ps_convert.F parser/init_convert.F parser/mod_parser_m.F parser/mod_itm.F wf_and_fft/wf_load.F wf_and_fft/WF_load.F wf_and_fft/space_inv_by_wf.F wf_and_fft/fft_setup.F wf_and_fft/.objects wf_and_fft/scatter_Bamp.F wf_and_fft/fft_3d.F wf_and_fft/WF_spatial_invertion.F wf_and_fft/scatter_Gamp.F pol_function/X_bare_RIM_GreenF.F pol_function/X_pre_setup.F pol_function/X_eh_setup.F pol_function/Dipole_kb_pp_comp.F pol_function/O_select_q_and_G.F pol_function/X_simple_GreenF.F pol_function/O_eels.F pol_function/X_delta_part.F pol_function/X_os.F pol_function/Dipole_kb_Ylm.F pol_function/O_driver.F pol_function/X_em1.F pol_function/Dipole_kb_sum.F pol_function/X_Gf.F pol_function/X_bare_RIM_Gf.F pol_function/.objects pol_function/Dipole_driver.F pol_function/Dipole_transverse.F pol_function/X_s.F pol_function/X_drude.F setup setup/QP_external_corrections.F setup/check_periodic_directions.F setup/setup.F setup/.objects setup/G_shells_finder.F setup/QP_state_table_setup.F coulomb/cutoff_sphere.F coulomb/cutoff_box.F coulomb/cutoff_driver.F coulomb/bessel_F2.F coulomb/bessel_F3.F coulomb/bessel_F4.F coulomb/bessel_J0.F coulomb/bessel_F5.F coulomb/bessel_F6.F coulomb/cutoff_cylinder.F coulomb/cutoff_test.F coulomb/rim.F coulomb/col_driver.F coulomb/rim_integrate.F interface/check_periodic_directions.F interface/parse_K_list.F interface/setup.F interface/QP_ctl_load.F interface/init_load.F interface/.objects interface/QP_init.F interface/init.F interface/QP_ctl_switch.F interface/G_shells_finder.F interface/stop_now.F interface/barriers.F interface/init_q_pts.F communicate/report_energies.F communicate/REPORT_Energies.F communicate/REPORT_Occupations.F communicate/.objects communicate/section.F communicate/fermi_msg.F communicate/acknowledge_yambo.F communicate/job_setup.F tddft/tddft_alda_r_space.F tddft/.objects tddft/tddft_do_X_W_typs.F tddft/tddft_alda_g_space.F io/variables_BS.F io/ioBS.F io/io_elemental.F io/io_HF_and_locXC.F io/ioDB1.F io/ioX.F io/IO_mute.F io/ioGROT.F io/ioWF.F io/ver_is_gt_or_eq.F io/IO_and_Messaging_switch.F io/ioQP.F io/ioXXVXC.F io/ioQINDX.F io/io_header.F io/ioRIM.F io/ioCOL_CUT.F io/ioBSS_Haydock.F io/Fragments_Synchronize.F io/ioOSTNTS.F io/ioBSS_diago.F io/ioKB_PP.F io/.objects io/ioE_RIM.F io/io_bulk.F io/ioDB1_reload.F io/IO_Setup.F io/Fragments_Restart.F modules/mod_electrons.F modules/mod_logo.F modules/mod_frequency.F modules/mod_zeros.F modules/mod_units.F modules/mod_BS.F modules/mod_TDDFT.F modules/mod_R_lattice.F modules/mod_memory.F modules/mod_functions.F modules/mod_wave_func.F modules/mod_drivers.F modules/mod_xc_functionals.F modules/mod_ACFDT.F modules/mod_fragments.F modules/mod_timing.F modules/mod_X_output.F modules/mod_vec_operate.F modules/mod_pseudo.F modules/std_presets.F modules/mod_par_indexes.F modules/mod_stderr.F modules/mod_par_proc.F modules/mod_pars.F modules/mod_X.F modules/mod_IO.F modules/mod_com.F modules/mod_FFT.F modules/mod_wrapper.F modules/mod_QP.F modules/mod_matrix_operate.F modules/.objects modules/mod_D_lattice.F modules/mod_global_XC.F modules/mod_collision.F bse/K_Haydock.F bse/K_exchange.F bse/K_scattering.F bse/K_dump_to_o_file.F bse/K_filling.F bse/Kernel_by_V.F bse/K_BSmat_by_V.F bse/.objects bse/K_driver.F bse/K_diagonalization.F bse/K_eh_setup.F bse/K_solvers.F bse/K_output_file.F bse/K.F ypp/ypp_init.F ypp/plot_xcrysden.F ypp/excitons_sort_and_report.F ypp/ypp_i.F ypp/k_grids.F ypp/.objects ypp/plot_gnuplot.F ypp/excitons_driver.F ypp/mod_YPP.F ypp/plot_check_and_launch.F ypp/electrons_driver.F ypp/ypp_init_load.F driver/p2y.h driver/a2y.h driver/editor.h.in driver/codever.h.in driver/e2y.h driver/yambo_driver.F driver/driver.c driver/yambo.h driver/ypp.h sbin/make_message.pl.in sbin/make_makefile.sh.in sbin/objects_debug.sh.in + + Additions: + - COHSEX support added to the GPL version + - Wrapper module for BLAS subroutines => Many (not all) caxpy,cgemm, cdotc calls replaced + - Added configure support for SUN compiler (-M instead of -I for module path; -I for INCLUDE file path) tested on x86_64. + - "-g" option to p2y in order to force p2y to use ALL RL-vectors. This means that the default p2y + behaviour has been changed. From now on p2y will write only the wavefuntion RL-vectors unless the + "-g" option is used. + - Added configure support for g95 under x86_64; for Open64 under x86_64. + - MacOSX support in configure + - Time-Dependent HF support added + - Large transferred momenta (outside IBZ) now possible + - Electronic DOS can now be calculated using ypp. + + Bugs: + - In k_lattice.F: fixed wrong definition of K-grid size in specific (and very few) cases + - Fixed compilation and configure problems with pgf95 + - Wrong DB loading in BSE when using PPA screened interaction + - CPP flags fixed + - Wrong precompiler options used to link SCALAPACK not fixed. + + Changes: + - Routines involved in the evaluation of the electronic occupations renamed. + - Fixed symmetries setup in interfaces + - Diagonalization DB I/O added in K_diagonalization.F (it has been removed by mistake) + - Global renaming: chosex -> cohsex + - Definition and use of database fragments deeply restyled. + - Gfortran support in configure + - **GLOBAL** restyling of messaging and output file naming in O_driver and + K_solvers + - **VERBOSITY** options changed. Use yambo -H to learn more. + + Patch sent by: Andrea + +2010-04-15 12:04:10 r3 svn + + Yambo GPL Release 3 version 3.2.1 (oldrev448) + + Additions: + - XC functional now recognized by p2y. + - Added self-detection of lattice geometry (devel) + - Support for Intel Fortrant v. 11 + - Q databases I/O to apply previously calculated + Q corrections on-fly + - During I/O databases are searched in all user provide + paths (-I and -O options). + - check in p2y of non symmorphic symmetries + - Acknowledge message at the end of each run + - A2Y support for Abinit 5.7.3 + - Added warning in case optical oscillators are loaded + to perform an optical calculation without including + non-local pseudo commutator + - In configure added : --enable-msgs-comps to exapand the compilation messages + --enable-options-check to allow the check of the options + passed to driver.c + By using --disable-options-check it is possible to solve the + annoying problem caused by some mpirun scripts (-p option) + - Added DOUBLE PRECISION slatec sources + - O_eels.F and check_periodic_directions.F. Used to allow the user + to specify non periodic directions. In case those are given + the polarizability is calculated by using KK transformations + even in the TDA. + - Added Time-Rev check using WFs to all branches. + - Added check for dlaran when using external lapack lib. + In case dlaran is missing an additional lib (lib/local) + is compiled and linked. + - Added check for C/Fortran precompiler in configure. + - Changed -with--netcdf with --with-netcdf-include AND + --with-netcdf-lib. Both must be passed to configure to + activate NETCDF support. + - Added V_testing verbosity level. + - X_os.F restiled to be processed by yamboo.pl. + - Verbosity to ypp. + + Changes: + - fermi_level.F and k_lattice.F definitively replaced + with the devel versions. + - Ifort optimization option "-ip" added + - Overall syncronization of trunk with and +SIN branches + - Removed different make clean*. Now only available "make clean" + and "make clean_all". In the first case the configuration files + are not removed and the code can be compiled without running + configure again. + - slatec routines used for polynomial fit replaced with + single/double versions. + - Input file name set to "(none)" when the file is not found/written + - Added set_real_printed_length rotuine in STDERR module in order + to tune on-fly the format of the output/report file for reals + - el_density.F moved in mod_electrons.F + - exc_properties renamed in excitons_driver.F + - Fixed use of "implicit none" in many modules + - plot_driver renamed in electrons_driver.F + - To avoid long preprocessor commands in driver.c, introduced two + new global variables (_YAMBO_MAIN _Y_MAIN) + - Deep restyling of ypp command-line options. + - "-xT" Intel compiler option replaced with "-xW" + - removed V_debug intger from mod_itm + - removed old interfaces/int_modules/mod_interfaces.F file + + Bugs: + - Double precision build fixed. + - Fixed wrong input file variables when using "-g n -p p" + command line options. + - Fixed minor bug in X_O_transverse + - Added in the devel/04_03_2009 a new test version of fermi_level.F + in response to the problem posted in the forum: + http://www.yambo-code.org/forum/viewtopic.php?f=12&t=18 + - Bug in K_driver.F. Fixed thanks to John awson + - Fix BU in ioKB_. Wrong VARIABE_NAME. + - In ypp (Claudio) + - C test failed when using IBM "XF" based compilers + - Wrong linking of BAS/AACK + - In ct_etime.c wrong declaration of tt variable in double + precision builds. + - qp_ctl: wrong transfer of gap correction obtained after a fit + - mod_IO: use of optional "subfolder" even if not passed + - ypp_i : wrong use of ioWF (not declared as part if wave_func module) + - Wrong dimension of c1vars in mod_itm was causing weird + behaviour when compiled using pgf90 + - Wrong use of integers in mod function in K_filling.F/K.F + - Wrong transferred momentum vector written to output files + - Unrecognized atomic elements (to fix) changed with '--' + in setup.F + - Bug fix in self-evaluation of the X energy range + +2010-04-15 12:00:03 r2 svn + + Yambo GPL Release 2 version 3.2.0 (oldrev315) + + Additions: + - Optimization of ifort 10 on x86_64 platforms + - qexml_v4.0.F added + - Removed -ifile option to ypp (appeared twice) + - doc/sample/bulk_silicon/experiment.dat + - ver_is_gt_or_eq function moved out of mod_IO so to add + optional revision argument. + - doc/ dir in GPL branch. + - the QuickGuidedTour.pdf in the doc/ dir. + - Added the sample core databases of + silicon used on the CPC paper and in the QuickGuidedTour.pdf + - Support for cygwin in Makefile (still to test) + - No other runlevel than setup is now permitted when + either G-shells or K-pts DBs are not present. + - Support for PWscf 4.0 in p2y + Changes: + - In GPL release, fermi level search permitted using a + finite temperature. + - Remove spurious option in ypp.h/ypp.h_gpl + - Cleaned the acx.m4,acx_mpi.m4 and configure.ac. Removed call + to preprocessor inquire, as it is already detected at the beginning. + - Changed tags in ioX.F to permit calculations above 999 Q-points. + - CPP replaced with "gcc -E" . This to avoid compatibility problems + with cpp. + Bugs: + - Wrong use of CPPFLAGS in configure fixed. This bug was preventing + some compilation with MPI support. + - In QP_driver QP_solver defined "by hand" to initialize the Dyson solver + - In GPL procedure to initialize the Dyson solver + - Preload of WF DB in ypp. Fragmentation parameters were not read at the + start-up => NETCDF I/O error in ypp for bands-fragmented DBs. + +2010-04-15 11:55:40 r1 svn + + Yambo GPL Release 1, version 3.1.2 (oldrev300) + + diff --git a/config/Makefile.in b/config/Makefile.in index 8f5f27f0f5..f873ecb03c 100644 --- a/config/Makefile.in +++ b/config/Makefile.in @@ -23,10 +23,11 @@ cpu = @build_cpu@ os = @build_os@ mpi = @mpi_cpp@ netcdf = @dnetcdf@ +libxc = @dlibxc@ scalapack = @dscalapack@ precision = @dp_cpp@ fft = @FFT_CPP@ -xcpp = @dnetcdf@ @mpi_cpp@ @FFT_CPP@ @dscalapack@ @dp_cpp@ +xcpp = @dnetcdf@ @dlibxc@ @mpi_cpp@ @FFT_CPP@ @dscalapack@ @dp_cpp@ debug = @enable_debug@ do_blas = @compile_blas@ do_lapack = @compile_lapack@ @@ -98,6 +99,17 @@ ypp: libs @LIBS2DO="$(YPPLIBS)" ; $(mksrc) @X2DO="ypp" ;XPATH="ypp";XLIBS="$(YPPLIBS)";$(mkx) # +# TESTs +# +tests: yambo + @if ! test -f ../branches/test-suite/run-tests.pl ; then echo "TESTs folder not found"; fi + @if test -f ../branches/test-suite/run-tests.pl ; then \ + cd ../branches/test-suite/; ./run-tests.pl -t b -test "Al111 01_init 02_eels 03_eps_drude 04_HF" -bindir=${PWD}/bin/ ; \ + cd ../branches/test-suite/; ./run-tests.pl -t b -test "C6H3Cl3 01_init 02_ALDA_diago 03_ALDA_haydo 04_QP 05_BSE 06_QP_RIM" -bindir=${PWD}/bin/; \ + cd ../branches/test-suite/; ./run-tests.pl -t b -test "hBN 01_init 02_HF 03_QP 04_BSE" -bindir=${PWD}/bin/; \ + cd ../branches/test-suite/; ./run-tests.pl -t b -test "LiF 01_init 02_RPA_no_LF 04_alda_g_space" -bindir=${PWD}/bin/; \ + fi +# clean: clean_fast clean_fast: @$(objects_clean) @@ -132,7 +144,7 @@ define mklib if test ! -f "$(libdir)/lib$$ldir.a" ; \ echo " " ; \ echo ">>>[Making $$ldir]<<<" ; \ - then ./sbin/make_makefile.sh $$DIR2GO/$$ldir lib$$ldir.a .objects l $(precision) $$ADF ; \ + then ./sbin/make_makefile.sh $$DIR2GO/$$ldir lib$$ldir.a .objects l $(precision) $(xcpp) $$ADF ; \ cd $$DIR2GO/$$ldir ; $(make) VPATH=$$DIR2GO/$$ldir || exit "$$?" ; cd ../../; fi \ done endef @@ -161,7 +173,7 @@ define objects_clean endef define lib_mod_clean find . \( -name '*.a' -o -name '*.mod' \) -type f -print | \ - grep -v netcdf | grep -v iotk | grep -v typesize | xargs rm -f + grep -v netcdf | grep -v xc_ | grep -v libxc | grep -v iotk | grep -v typesize | xargs rm -f echo "[CLEAN] Libraries ... done" echo "[CLEAN] Modules ... done" endef diff --git a/config/configure.ac b/config/configure.ac index 72a52b4db2..dd2a323c9e 100644 --- a/config/configure.ac +++ b/config/configure.ac @@ -19,11 +19,11 @@ # Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, # MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. # -AC_INIT(Yambo, 3.2.2 r.633 , yambo@yambo-code.org) +AC_INIT(Yambo, 3.2.3 r.676 , yambo@yambo-code.org) SVERSION="3" SPATCHLEVEL="2" -SSUBLEVEL="2" -SREVISION="633" +SSUBLEVEL="3" +SREVISION="676" AC_SUBST(SVERSION) AC_SUBST(SPATCHLEVEL) AC_SUBST(SSUBLEVEL) diff --git a/configure b/configure index 8cc7790224..07b3f5d7d1 100755 --- a/configure +++ b/configure @@ -1,6 +1,6 @@ #! /bin/sh # Guess values for system-dependent variables and create Makefiles. -# Generated by GNU Autoconf 2.63 for Yambo 3.2.2 r.633 . +# Generated by GNU Autoconf 2.63 for Yambo 3.2.3 r.676 . # # Report bugs to . # @@ -596,8 +596,8 @@ SHELL=${CONFIG_SHELL-/bin/sh} # Identity of this package. PACKAGE_NAME='Yambo' PACKAGE_TARNAME='yambo' -PACKAGE_VERSION='3.2.2 r.633 ' -PACKAGE_STRING='Yambo 3.2.2 r.633 ' +PACKAGE_VERSION='3.2.3 r.676 ' +PACKAGE_STRING='Yambo 3.2.3 r.676 ' PACKAGE_BUGREPORT='yambo@yambo-code.org' ac_unique_file="driver/driver.c" @@ -1370,7 +1370,7 @@ if test "$ac_init_help" = "long"; then # Omit some internal or obsolete options to make the list less imposing. # This message is too long to be a string in the A/UX 3.1 sh. cat <<_ACEOF -\`configure' configures Yambo 3.2.2 r.633 to adapt to many kinds of systems. +\`configure' configures Yambo 3.2.3 r.676 to adapt to many kinds of systems. Usage: $0 [OPTION]... [VAR=VALUE]... @@ -1435,7 +1435,7 @@ fi if test -n "$ac_init_help"; then case $ac_init_help in - short | recursive ) echo "Configuration of Yambo 3.2.2 r.633 :";; + short | recursive ) echo "Configuration of Yambo 3.2.3 r.676 :";; esac cat <<\_ACEOF @@ -1556,7 +1556,7 @@ fi test -n "$ac_init_help" && exit $ac_status if $ac_init_version; then cat <<\_ACEOF -Yambo configure 3.2.2 r.633 +Yambo configure 3.2.3 r.676 generated by GNU Autoconf 2.63 Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001, @@ -1570,7 +1570,7 @@ cat >config.log <<_ACEOF This file contains any messages produced by compilers while running configure, to aid debugging if configure makes a mistake. -It was created by Yambo $as_me 3.2.2 r.633 , which was +It was created by Yambo $as_me 3.2.3 r.676 , which was generated by GNU Autoconf 2.63. Invocation command line was $ $0 $@ @@ -1940,8 +1940,8 @@ ac_compiler_gnu=$ac_cv_c_compiler_gnu SVERSION="3" SPATCHLEVEL="2" -SSUBLEVEL="2" -SREVISION="633" +SSUBLEVEL="3" +SREVISION="676" @@ -13643,7 +13643,7 @@ exec 6>&1 # report actual input values of CONFIG_FILES etc. instead of their # values after options handling. ac_log=" -This file was extended by Yambo $as_me 3.2.2 r.633 , which was +This file was extended by Yambo $as_me 3.2.3 r.676 , which was generated by GNU Autoconf 2.63. Invocation command line was CONFIG_FILES = $CONFIG_FILES @@ -13693,7 +13693,7 @@ Report bugs to ." _ACEOF cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1 ac_cs_version="\\ -Yambo config.status 3.2.2 r.633 +Yambo config.status 3.2.3 r.676 configured by $0, generated by GNU Autoconf 2.63, with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`\\" diff --git a/driver/driver.c b/driver/driver.c index 83fb3f1bfa..720f2ca11e 100644 --- a/driver/driver.c +++ b/driver/driver.c @@ -67,6 +67,9 @@ typedef struct #if defined _a2y #include "a2y.h" #endif +#if defined _c2y + #include "c2y.h" +#endif #if defined _f2y #include "f2y.h" #endif @@ -268,6 +271,14 @@ main(int argc, char *argv[]) F90_FUNC(ypp_i,YPP_I)( rnstr2,&lni,inf,&iif,id,&iid,od,&iod,com_dir,&icd,js,&ijs,&np,&pid); #endif +#if defined _c2y + /* + Running the Fortran c2y driver + =========================================================================== + */ + F90_FUNC(c2y_i,C2Y_I)( + rnstr2,&lni,inf,&iif,id,&iid,od,&iod,com_dir,&icd,js,&ijs,&np,&pid); +#endif #if defined _a2y /* Running the Fortran a2y driver diff --git a/driver/yambo.h b/driver/yambo.h index db7aec9b7b..ebcc21f82a 100644 --- a/driver/yambo.h +++ b/driver/yambo.h @@ -36,23 +36,27 @@ {"help", "h","Short Help",0,0,0,0}, {"lhelp", "H","Long Help",0,0,0,0}, {"jobstr", "J","Job string identifier",0,0,1,0}, - {"infver", "V","Input file verbosity [opt=RL,kpt,sc,qp,io,gen,resp]",0,0,1,0}, + {"infver", "V","Input file verbosity",0,0,1,0}, + {"DESC", " ","opt=RL,kpt,sc,qp,io,gen,resp,rt,all",0,0,0,0}, {"ifile", "F","Input file",0,0,1,0}, - {"idir", "I","Core I/O directory",0,0,1,0}, - {"odir", "O","Additional I/O directory",0,0,1,0}, - {"cdir", "C","Communications I/O directory",0,0,1,0}, + {"idir", "I","Core I/O directory",0,0,1,0}, + {"odir", "O","Additional I/O directory",0,0,1,0}, + {"cdir", "C","Communications I/O directory",0,0,1,0}, {"nompi", "N","Skip MPI initialization",0,0,0,0}, {"dbpr", "D","DataBases properties",0,0,0,0}, {"dbfrag", "S","DataBases fragmentation",0,0,0,0}, {"setup", "i","Initialization",0,0,0,0}, {"optics", "o","Optics [opt=(c)hi/(b)se/(t)dhf]",0,0,1,0}, - {"tddft", "t","The TDDFTs [opt=(a)LDA/(l)RC]",0,0,1,0}, - {"rim_cut","c","Coulomb interaction",0,0,0,0}, - {"HF_and_locXC", "x","Hartree-Fock Self-energy and local XC",0,0,0,0}, - {"em1s", "b","Static Inverse Dielectric Matrix",0,0,0,0}, - {"gwapprx","p","GW approximations [opt=(p)PA/c(HOSEX)]",0,0,1,0}, - {"gw0", "g","Dyson Equation solver [opt=n(ewton)]",0,0,1,0}, - {"bss", "y","BSE solver [opt=h/d]",0,0,1,0}, + {"tddft", "t","The TDDFTs [opt=(a)LDA/(l)RC]",0,0,1,0}, + {"rim_cut","c","Coulomb interaction",0,0,0,0}, + {"HF_and_locXC", "x","Hartree-Fock Self-energy and local XC",0,0,0,0}, + {"em1d", "d","Dynamical Inverse Dielectric Matrix",0,0,0,0}, + {"em1s", "b","Static Inverse Dielectric Matrix",0,0,0,0}, + {"gwapprx","p","GW approximations [opt=(p)PA/c(HOSEX)]",0,0,1,0}, + {"gw0", "g","Dyson Equation solver",0,0,1,0}, + {"DESC", " ","opt=n(ewton)/s(ecant)/g(reen)",0,0,0,0}, + {"life", "l","GoWo Quasiparticle lifetimes",0,0,0,0}, + {"bss", "y","BSE solver [opt=h/d/i/t]",0,0,1,0}, {NULL,NULL,NULL,0,0,0,0} }; char *tool="yambo"; diff --git a/driver/yambo_driver.F b/driver/yambo_driver.F index 6e8e9ec9e4..649149bdb0 100644 --- a/driver/yambo_driver.F +++ b/driver/yambo_driver.F @@ -162,6 +162,12 @@ integer function yambo_driver(instr,lnstr,inf,iinf,ind,& ! if (l_pp_em1) i_err=X_em1(Xen,Xk,q,X(4),Xw(4),.false.) ! + ! Dynamical + !----------- + ! + l_dyn_em1= (.not.l_ppa.and..not.l_cohsex) .and. (l_em1d.or.(l_gw0.and.l_el_corr)) + ! + if (l_dyn_em1) i_err=X_em1(Xen,Xk,q,X(3),Xw(3),.false.) ! ! EXTENDED COLLISIONS (HF & COHSEX) !=================================== @@ -188,6 +194,10 @@ integer function yambo_driver(instr,lnstr,inf,iinf,ind,& ! else ! + ! Real Axis + !============ + ! + call QP_driver(X(3),Xen,Xk,en,k,q,Xw(3)) ! endif ! diff --git a/driver/ypp.h b/driver/ypp.h index 2b1d6028ca..2f6d36dfcb 100644 --- a/driver/ypp.h +++ b/driver/ypp.h @@ -20,7 +20,7 @@ MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. */ /* - Driver declaration + Driver */ #if defined _FORTRAN_US int ypp_i @@ -36,14 +36,14 @@ {"help", "h","Short Help",0,0,0,0}, {"lhelp", "H","Long Help",0,0,0,0}, {"jobstr","J","Job string identifier",0,0,1,0}, - {"infver", "V","Input file verbosity [opt=gen]",0,0,1,0}, + {"infver", "V","Input file verbosity [opt=gen,qp]",0,0,1,0}, {"ifile", "F","Input file",0,0,1,0}, {"idir", "I","Core I/O directory",0,0,1,0}, {"odir", "O","Additional I/O directory",0,0,1,0}, - {"cdir", "C","Communications I/O directory",0,0,1,0}, + {"cdir", "C","Communications I/O directory",0,0,1,0}, {"nompi", "N","Skip MPI initialization",0,0,0,0}, {"dbfrag","S","DataBases fragmentation",0,0,0,0}, - {"bzgrids","k","BZ Grid generator [(k)pt,(q)pt,(l)ongitudinal]",0,0,1,0}, + {"bzgrids","k","BZ Grid generator [(k)pt,(q)pt]",0,0,1,0}, {"excitons", "e","Excitons [(s)ort,(a)mplitude,(w)ave]",0,0,1,0}, {"electrons","s","Electrons [(w)ave,(d)ensity,do(s)]",0,0,1,0}, {"freehole","f","Free hole position [excitons plot]",0,0,0,0}, diff --git a/include/version.inc b/include/version.inc index b129123e3a..ae6c9bba46 100644 --- a/include/version.inc +++ b/include/version.inc @@ -1,4 +1,4 @@ code_version(1)=3 code_version(2)=2 -code_version(3)=2 -code_revision=633 +code_version(3)=3 +code_revision=676 diff --git a/interfaces/a2y/a2y_db1.F b/interfaces/a2y/a2y_db1.F index 7b9ae7db55..7e947106a2 100644 --- a/interfaces/a2y/a2y_db1.F +++ b/interfaces/a2y/a2y_db1.F @@ -33,9 +33,8 @@ subroutine a2y_db1(en,k,KSS_file_name) use wave_func, ONLY: wf_nc_k, wf_igk,wf_ncx,wf_ng use vec_operate, ONLY: cross_product,v_is_zero use defs_datatypes, ONLY: hdr_type, wffile_type - use xc_functionals, ONLY: GS_xc_KIND,GS_xc_FUNCTIONAL,XC_USED_LIBXC - use xc_functionals, ONLY: XC_EXCHANGE_CORRELATION,XC_EXCHANGE,XC_LDA_C_PW,XC_LDA_C_PZ,& -& XC_LDA_C_WIGNER,XC_LDA_C_HL,XC_LDA_C_XALPHA,NONE,XC_LDA_X + use xc_functionals, ONLY: GS_xc_KIND,GS_xc_FUNCTIONAL + use mod_xc2y, ONLY: XC_Abinit2Yambo,XC_Libxc2Yambo use mod_com2y, ONLY: print_interface_dimensions,symmetries_check_and_load,& & alat_mult_factor,artificial_spin_pol implicit none @@ -74,37 +73,11 @@ subroutine a2y_db1(en,k,KSS_file_name) ! XC KIND/FUNCTIONAL !=================== ! - select case (ahdr%ixc) - case (0) - GS_xc_FUNCTIONAL=NONE - GS_xc_KIND=NONE - case (1,7) - GS_xc_FUNCTIONAL=XC_LDA_C_PW - GS_xc_KIND=XC_EXCHANGE_CORRELATION - case (2) - GS_xc_FUNCTIONAL=XC_LDA_C_PZ - GS_xc_KIND=XC_EXCHANGE_CORRELATION - case (3) - GS_xc_FUNCTIONAL=XC_LDA_C_PZ - GS_xc_KIND=XC_EXCHANGE_CORRELATION - case (4) - GS_xc_FUNCTIONAL=XC_LDA_C_WIGNER - GS_xc_KIND=XC_EXCHANGE_CORRELATION - case (5) - GS_xc_FUNCTIONAL=XC_LDA_C_HL - GS_xc_KIND=XC_EXCHANGE_CORRELATION - case (6) - GS_xc_FUNCTIONAL=XC_LDA_C_XALPHA - GS_xc_KIND=XC_EXCHANGE_CORRELATION - case (8) - GS_xc_FUNCTIONAL=XC_EXCHANGE - GS_xc_KIND=XC_LDA_X - case (9:) - call warning ('GGA not supported: switching to LDA PW') - GS_xc_FUNCTIONAL = XC_LDA_C_PW - GS_xc_KIND = XC_EXCHANGE_CORRELATION - end select - ! + if (ahdr%ixc.ge. 0) then + call XC_Abinit2Yambo(ahdr%ixc,GS_xc_FUNCTIONAL,GS_xc_KIND) + else !negative ahdr%ixc means abinit used libxc + call XC_Libxc2Yambo(ahdr%ixc,GS_xc_FUNCTIONAL,GS_xc_KIND) + end if ! n_spinor=1 n_sp_pol=1 diff --git a/interfaces/int_modules/.objects b/interfaces/int_modules/.objects index 3648bcdf31..971b01c781 100644 --- a/interfaces/int_modules/.objects +++ b/interfaces/int_modules/.objects @@ -1 +1 @@ -objs = mod_com2y.o mod_wf2y.o +objs = mod_com2y.o mod_wf2y.o mod_xc2y.o diff --git a/interfaces/int_modules/mod_xc2y.F b/interfaces/int_modules/mod_xc2y.F new file mode 100644 index 0000000000..57b23a5d3e --- /dev/null +++ b/interfaces/int_modules/mod_xc2y.F @@ -0,0 +1,193 @@ +! +! Copyright (C) 2000-2010 Myrta Gruning and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +module mod_xc2y + ! + use xc_functionals + implicit none + ! + integer, parameter :: N_Abinit_XC = 28, N_LibXC_LDA= 20 +! integer, parameter :: N_PW_X = 9, N_PW_C = 12, N_PW_GX = 12,& +! & N_PW_GC = 9 + ! + integer, parameter, dimension(N_Abinit_XC) :: XC_A2Y =& +& (/XC_LDA_XC_TETER93, & !ixc= 1 + XC_LDA_C_PZ, & !ixc= 2 + XC_LDA_C_PZ, & !ixc= 3 (???) + XC_LDA_C_WIGNER, & !ixc= 4 + XC_LDA_C_HL, & !ixc= 5 + XC_LDA_C_XALPHA, & !ixc= 6 + XC_LDA_C_PW, & !ixc= 7 + XC_LDA_X, & !ixc= 8 + XC_LDA_C_PW, & !ixc= 9 (???) + XC_NOT_AVAILABLE, & !ixc= 10 + XC_NOT_AVAILABLE, & !ixc= 11 + XC_NOT_AVAILABLE, & !ixc= 12 + XC_NOT_AVAILABLE, & !ixc= 13 + XC_NOT_AVAILABLE, &!ixc= 14 (???) + XC_NOT_AVAILABLE, &!ixc= 15 (???) + XC_NOT_AVAILABLE, &!ixc= 16 + XC_NOT_AVAILABLE, &!ixc= 17 + XC_NOT_AVAILABLE, &!ixc= 18 + XC_NOT_AVAILABLE, &!ixc= 19 + XC_NOT_AVAILABLE, &!ixc= 20 + XC_NOT_AVAILABLE, &!ixc= 21 + XC_NOT_AVAILABLE, &!ixc= 22 + XC_NOT_AVAILABLE, &!ixc= 23 (???) + XC_NOT_AVAILABLE, &!ixc= 24 + XC_NOT_AVAILABLE, &!ixc= 25 + XC_NOT_AVAILABLE, &!ixc= 26 + XC_NOT_AVAILABLE, &!ixc= 27 + XC_NOT_AVAILABLE/) !ixc= 28 + integer, parameter, dimension(N_LibXC_LDA) :: LDA_Libxc2Y =& + & (/XC_LDA_X, & !1 + XC_LDA_C_WIGNER, & !2 + XC_LDA_C_RPA, & !3 + XC_LDA_C_HL, & !4 + XC_LDA_C_GL, & !5 + XC_LDA_C_XALPHA, & !6 + XC_LDA_C_VWN, & !7 + XC_LDA_C_VWN_RPA,& !8 + XC_LDA_C_PZ, & !9 + XC_LDA_C_PZ_MOD, & !10 + XC_LDA_C_OB_PZ, & !11 + XC_LDA_C_PW, & !12 + XC_LDA_C_PW_MOD, & !13 + XC_LDA_C_OB_PW, & !14 + XC_NOT_AVAILABLE,& !15 + XC_NOT_AVAILABLE,& !16 + XC_LDA_C_vBH, & !17 + XC_NOT_AVAILABLE,& !18 + XC_NOT_AVAILABLE,& !19 + XC_LDA_XC_TETER93/) !20 +contains + subroutine XC_Abinit2Yambo(abint_func,yambo_func,yambo_kind) + integer abint_func, yambo_func,yambo_kind + ! + yambo_func = XC_NOT_AVAILABLE + yambo_kind = NONE + ! + if (abint_func == 0) then + yambo_func = NONE + else + if (abint_func.le.N_Abinit_XC) yambo_func = XC_A2Y(abint_func) + yambo_kind = XC_EXCHANGE + if (yambo_func.ne.XC_LDA_X) yambo_kind = yambo_kind + XC_CORRELATION + end if + ! + end subroutine XC_Abinit2Yambo + ! + subroutine XC_Libxc2Yambo(libxc_func,yambo_func,yambo_kind) + integer libxc_func,func(2),yambo_func,yambo_kind,ii,ixc + integer, parameter :: factor = 1000 + ! + ! This serves when abinit was used with libxc + ! in that case, ixc < 0 + ! + ixc = libxc_func + if (libxc_func<0) ixc = -libxc_func + ! + func(1) = ixc/factor + func(2) = ixc-func(1)*factor + yambo_func = XC_NOT_AVAILABLE + yambo_kind = NONE + ! + do ii=1,2 + select case (func(ii)) + case(0) + if (yambo_func==NONE.or.yambo_func==XC_NOT_AVAILABLE)& + & yambo_func = NONE + case(1) + yambo_func = LDA_Libxc2Y(func(ii)) + yambo_kind = yambo_kind + XC_EXCHANGE + case(2:N_LibXC_LDA) + yambo_func = LDA_Libxc2Y(func(ii)) + yambo_kind = yambo_kind + XC_CORRELATION + if (yambo_func==XC_LDA_XC_TETER93) yambo_kind = XC_EXCHANGE_CORRELATION + case DEFAULT + yambo_func = XC_NOT_AVAILABLE + yambo_kind = NONE + end select + if (yambo_func==XC_NOT_AVAILABLE) EXIT + end do + ! + end subroutine XC_Libxc2Yambo + ! +end module mod_xc2y +! +! FOR PWSCF: FOR LATER USE +! +! character(len=6),parameter,dimension(N_PW_X) :: PW_X=& +! & (/' NOX', ' SLA', ' SL1', ' RXC', ' OEP', ' HF', 'PB0X', 'B3LP', ' KZK'/) +! character(len=6),parameter,dimension(N_PW_C) :: PW_C=& +! & (/ ' NOC', ' PZ', ' VWN', ' LYP', ' PW', ' WIG', ' HL', ' OBZ', & +! ' OBW', ' GL' , 'B3LP', ' KZK' /) +! character(len=6),parameter,dimension(N_PW_GX) :: PW_GX=& +! & (/ 'NOGX', ' B88', ' GGX', ' PBX', ' RPB', 'HCTH', 'OPTX',& +! 'META', 'PB0X', 'B3LP',' PSX', ' WCX'/) +! character(len=6),parameter,dimension(N_PW_GC) :: PW_GC=& +! & (/ 'NOGC', ' P86', ' GGC', 'BLYP', ' PBC', 'HCTH', 'META',& +! 'B3LP', ' PSC'/) +! integer,parameter,dimension(N_PW_X) :: PWX2LIBXC=& +! & (/ NOXC, XC_LDA_X, XC_LDA_X*XC_FACTOR+XC_LDA_C_XALPHA, XC_NOT_AVAILABLE, & +! & XC_NOT_AVAILABLE, XC_NOT_AVAILABLE, XC_NOT_AVAILABLE, XC_NOT_AVAILABLE,& +! XC_NOT_AVAILABLE/) +! integer,parameter,dimension(N_PW_C) :: PWC2LIBXC=& +! & (/ NOXC, XC_LDA_C_PZ, XC_LDA_C_VWN, XC_GGA_C_LYP, XC_LDA_C_PW,& +! & XC_LDA_C_WIGNER, XC_LDA_C_HL, XC_LDA_C_OB_PZ, XC_LDA_C_OB_PW,XC_LDA_C_GL,& +! & XC_NOT_AVAILABLE, XC_NOT_AVAILABLE/) +! integer,parameter,dimension(N_PW_GX) :: PWGX2LIBXC=& +! & (/ NOXC, XC_GGA_X_B88, XC_GGA_X_PW91, XC_GGA_X_PBE, XC_GGA_X_PBE_R,& +! & XC_GGA_XC_HCTH_93, XC_GGA_X_OPTX, XC_MGGA_X_TPSS, XC_HYB_GGA_XC_PBEH,& +! & XC_HYB_GGA_XC_B3LYP, XC_GGA_X_PBE_SOL, XC_GGA_X_WC/) +! integer,parameter,dimension(N_PW_GC) :: PWGC2LIBXC=& +! & (/ NOXC, XC_GGA_C_P86, XC_GGA_C_PW91, XC_GGA_C_LYP, XC_GGA_C_PBE, XC_GGA_XC_HCTH_93,& +! & XC_MGGA_C_TPSS, XC_HYB_GGA_XC_B3LYP, XC_GGA_C_PBE_SOL/) +! end module + + +! character(len=6),parameter,dimension(N_PW_XC) :: PW_XC=& +!& (/ 'BP', 'BLYP', 'PBE', 'REVPBE', 'PBESOL', 'HCTH', 'OLYP', & +!& 'TPSS', 'WC', 'PBE0', 'B3LYP'/) +! integer,parameter,dimension(N_PW_XC) :: PWXC2LIBXC=& +!& (/ XC_GGA_X_B88*XC_FACTOR+XC_GGA_C_P86,& +!& XC_GGA_X_PW91*XC_FACTOR+XC_GGA_C_PW91,& +!& XC_GGA_X_PBE*XC_FACTOR+XC_GGA_C_PBE,& +!& XC_GGA_X_PBE_R*XC_FACTOR+XC_GGA_C_PBE,& +!& XC_GGA_X_PBE_SOL*XC_FACTOR+XC_GGA_C_PBE_SOL,& +!& XC_GGA_XC_HCTH_120,& +!& XC_GGA_X_OPTX*XC_FACTOR+XC_GGA_C_LYP,& +!& XC_MGGA_X_TPSS*XC_FACTOR+XC_MGGA_C_TPSS,& +!& XC_GGA_X_WC*XC_FACTOR+XC_GGA_C_PBE,& +!& XC_HYB_GGA_XC_B3LYP/) +! +! For LATER USE: ABINIT GGAs +! + ! XC_GGA_X_PBE*XC_FACTOR+XC_GGA_C_PBE, & !ixc= 11 + ! XC_GGA_X_PBE*XC_FACTOR, & !ixc= 12 + ! XC_GGA_XC_LB*XC_FACTOR, & !ixc= 13 + ! XC_GGA_X_PBE_R*XC_FACTOR+XC_GGA_C_PBE,&!ixc= 14 (???) + ! XC_GGA_X_RPBE*XC_FACTOR+XC_GGA_C_PBE, &!ixc= 15 (???) + ! XC_GGA_XC_HCTH_93*XC_FACTOR, &!ixc= 16 + ! XC_GGA_XC_HCTH_120*XC_FACTOR, &!ixc= 17 + ! XC_GGA_X_WC*XC_FACTOR, &!ixc= 23 (???) + ! XC_GGA_XC_HCTH_147*XC_FACTOR, &!ixc= 26 + ! XC_GGA_XC_HCTH_407*XC_FACTOR, &!ixc= 27 diff --git a/interfaces/p2y/p2y_db1.F b/interfaces/p2y/p2y_db1.F index bdac8ae503..95f8f57ffc 100644 --- a/interfaces/p2y/p2y_db1.F +++ b/interfaces/p2y/p2y_db1.F @@ -30,7 +30,6 @@ subroutine p2y_db1(en,k) implicit none type(levels), intent(out) :: en ! Energies type(bz_samp), intent(out) :: k ! K/Q points - !---------------------------------------------------------------------* ! Read dimensions * !---------------------------------------------------------------------* diff --git a/sbin/make_makefile.sh.in b/sbin/make_makefile.sh.in index 2d5f43aea8..5ae9639a11 100755 --- a/sbin/make_makefile.sh.in +++ b/sbin/make_makefile.sh.in @@ -24,6 +24,7 @@ os="@build_os@" cpp="@CPP@" cppflags="@C_AS_CPP_FLAGS@" ECHO_N="@ECHO_N@" +PREFIX="@MKMF_PREFIX@" # if [ $# = 0 ] ; then echo $0 "dir target objectfile mode(l/x) Dflag1 Dflag2 Dflag3 ..." @@ -77,7 +78,7 @@ case $target in llibs="$llibs \$(lscalapack) \$(llapack) \$(lblacs) \$(lblas)" llibs="$llibs \$(lnetcdf) \$(llibxc) -lslatec \$(llocal) \$(lfftw) -lm $llmpi" ;; - a2y|elk2y) + a2y|elk2y|c2y) llibs="-lint_modules $llibs \$(lscalapack) \$(llapack) \$(lblacs) \$(lblas) \$(lnetcdf) \$(llibxc) \$(llocal) \$(lfftw) -lm $llmpi" ;; p2y*) diff --git a/src/bse/K.F b/src/bse/K.F index f0d211b1d6..673a0c1172 100644 --- a/src/bse/K.F +++ b/src/bse/K.F @@ -259,7 +259,7 @@ subroutine K(iq,Ken,Xk,q,X,Xw,W_bss) ! if (BS_K_is_ALDA) then allocate(F_xc(fft_size)) - call xc_lda_driver(Ken,Xk,WF_KIND,WF_xc_functional,2) + call XC_potential_driver(Ken,Xk,WF_KIND,WF_xc_functional,2) endif ! if (BS_res_K_corr) then @@ -472,6 +472,15 @@ subroutine K(iq,Ken,Xk,q,X,Xw,W_bss) if (BS_res_K_exchange) H_res_x=cdotc(BS_n_g_exch,O2x(1,icv2),1,O1x(1,icv1),1) #endif ! + ! :::Exchange (coupling)::: + ! +#if defined _DOUBLE + if (BS_cpl_K_exchange) H_cpl_x=zdotu(BS_n_g_exch,& +& O2x( g_rot(inv_index,1:BS_n_g_exch),icv2),1,O1x(1,icv1),1) +#else + if (BS_cpl_K_exchange) H_cpl_x=cdotu(BS_n_g_exch,& +& O2x( g_rot(inv_index,1:BS_n_g_exch),icv2),1,O1x(1,icv1),1) +#endif ! ! :::ALDA (resonant)::: ! @@ -479,6 +488,11 @@ subroutine K(iq,Ken,Xk,q,X,Xw,W_bss) & tddft_alda_R_space(iq,(/ic1,ic2/),(/ik1,ik2,ik1,ik2/),& & (/iv1,iv2/),(/is1,is2,is1,is2/),1) ! + ! :::ALDA (coupling)::: + ! + if (BS_K_is_ALDA) H_cpl_x=H_cpl_x+& +& tddft_alda_R_space(iq,(/ic1,ic2/),(/ik1,ik2,ik1,ik2/),& +& (/iv1,iv2/),(/is1,is2,is1,is2/),2) ! ! :::Correlation (resonant)::: ! @@ -520,14 +534,52 @@ subroutine K(iq,Ken,Xk,q,X,Xw,W_bss) ! endif ! + ! :::Correlation (coupling)::: ! - ! Impose the kernel to be hermitian + if (BS_cpl_K_corr.and.i_sp1==i_sp2) then + ! + iOcv=O_table(ic1-BS_bands(1)+1,is1, iv2-BS_bands(1)+1,is2,i_sp1) + iOvc=O_table(iv1-BS_bands(1)+1,is1p,ic2-BS_bands(1)+1,is2p,i_sp2) + inv_s1=sop_inv(is1) + is1xis2=sop_tab(sop_inv(is1),is2) + ! + forall(i1=1:BS_n_g_W) Ocv(i1)=& + & BS_O(g_rot(inv_s1,G_m_G(g_rot(iqs,i1),ig0)),iOcv)*O_phase(is1xis2,iv2,i_sp1) + forall(i1=1:BS_n_g_W) Ovc(i1)=& + & BS_O(g_rot(inv_s1,G_m_G(g_rot(iqs,i1),ig0)),iOvc)*O_phase(is1xis2,ic2,i_sp2) + ! + if (is1>nsym/(i_time_rev+1)) Ocv=conjg(Ocv) + if (is1>nsym/(i_time_rev+1)) Ovc=conjg(Ovc) + if (BS_W_is_diagonal) then + ! + forall(i1=1:BS_n_g_W) Cc(i1)=Ocv(i1)*BS_W(i1,1,iq_W) + ! + else + ! + do i1=1,BS_n_g_W +#if defined _DOUBLE + Cc(i1)=zdotu(BS_n_g_W,Ocv,1,BS_W(1,i1,iq_W),1) +#else + Cc(i1)=cdotu(BS_n_g_W,Ocv,1,BS_W(1,i1,iq_W),1) +#endif + enddo + endif + ! +#if defined _DOUBLE + H_cpl_c=zdotc(BS_n_g_W,Ovc,1,Cc,1)*4.*pi +#else + H_cpl_c=cdotc(BS_n_g_W,Ovc,1,Cc,1)*4.*pi +#endif + ! + endif + ! + ! Impose the resonant part of the kernel to be hermitian ! if (ik1==ik2.and.icv1==icv2) H_res_c=real(H_res_c) ! BS_mat(icv1,icv2)=H_res_x*real(spin_occ)*Co-H_res_c ! - ! + if (BS_K_coupling) BS_cpl_mat(icv1,icv2)=H_cpl_x*real(spin_occ)*Co-H_cpl_c enddo enddo call pp_redux_wait(BS_mat) @@ -579,8 +631,10 @@ subroutine K(iq,Ken,Xk,q,X,Xw,W_bss) call mem_est("O_RES_WS") deallocate(BS_O) call mem_est("BS_O") - ! - ! + if (BS_cpl_K_corr) then + deallocate(Ovc,Ocv) + call mem_est("O_CPL_WS") + endif endif ! if (l_bs_fxc) then diff --git a/src/bse/K_BSmat_by_V.F b/src/bse/K_BSmat_by_V.F index dd42d01eaf..180facee82 100644 --- a/src/bse/K_BSmat_by_V.F +++ b/src/bse/K_BSmat_by_V.F @@ -160,6 +160,32 @@ subroutine K_BSmat_by_V(iq,iter,Vi,Vo) & Vo(frag_pointer(1)+1),1) #endif ! + if (BS_K_coupling) then + ! + ! Calculate respectively: the diagonal anti-resonant, coupling, anti-coupling part + ! +#if defined _DOUBLE + call zgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik1),(-1._SP,0._SP),& +& conjg(BS_mat),BS_blk_dim(ik1),Vi(frag_pointer(1)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+BS_K_dim+1),1) + call zgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik1),(1._SP,0._SP),& +& BS_cpl_mat,BS_blk_dim(ik1),Vi(frag_pointer(1)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+1),1) + call zgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik1),(-1._SP,0._SP),& +& conjg(BS_cpl_mat),BS_blk_dim(ik1),Vi(frag_pointer(1)+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+BS_K_dim+1),1) +#else + call cgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik1),(-1._SP,0._SP),& +& conjg(BS_mat),BS_blk_dim(ik1),Vi(frag_pointer(1)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+BS_K_dim+1),1) + call cgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik1),(1._SP,0._SP),& +& BS_cpl_mat,BS_blk_dim(ik1),Vi(frag_pointer(1)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+1),1) + call cgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik1),(-1._SP,0._SP),& +& conjg(BS_cpl_mat),BS_blk_dim(ik1),Vi(frag_pointer(1)+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+BS_K_dim+1),1) +#endif + end if ! end if ! @@ -175,6 +201,11 @@ subroutine K_BSmat_by_V(iq,iter,Vi,Vo) & i2=BS_k_and_row_restart(2)+1:BS_k_and_row_restart(2)+BS_blk_dim(ik2)) & & BS_mat(i2,i1)=conjg(BS_mat(i1,i2)) ! + if (BS_K_coupling) then + forall(i1=BS_k_and_row_restart(1)+1:BS_k_and_row_restart(1)+BS_blk_dim(ik1),& +& i2=BS_k_and_row_restart(2)+1:BS_k_and_row_restart(2)+BS_blk_dim(ik2)) & +& BS_mat(i2,BS_K_dim+i1)=BS_mat(i1,BS_K_dim+i2) + endif ! else ! @@ -194,6 +225,50 @@ subroutine K_BSmat_by_V(iq,iter,Vi,Vo) & Vo(frag_pointer(2)+1),1) #endif ! + if (BS_K_coupling) then + ! + ! Calculate respectively: the off-diagonal anti-resonant, coupling, anti-coupling part + ! +#if defined _DOUBLE + call zgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik2),(-1._SP,0._SP),& +& conjg(BS_mat),BS_blk_dim(ik1),Vi(frag_pointer(2)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+BS_K_dim+1),1) + call zgemv('c',BS_blk_dim(ik1),BS_blk_dim(ik2),(-1._SP,0._SP),& +& conjg(BS_mat),BS_blk_dim(ik1),Vi(frag_pointer(1)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(2)+BS_K_dim+1),1) + call zgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik2),(1._SP,0._SP),& +& BS_cpl_mat,BS_blk_dim(ik1),Vi(frag_pointer(2)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+1),1) + call zgemv('t',BS_blk_dim(ik1),BS_blk_dim(ik2),(1._SP,0._SP),& +& BS_cpl_mat,BS_blk_dim(ik1),Vi(frag_pointer(1)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(2)+1),1) + call zgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik2),(-1._SP,0._SP),& +& conjg(BS_cpl_mat),BS_blk_dim(ik1),Vi(frag_pointer(2)+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+BS_K_dim+1),1) + call zgemv('t',BS_blk_dim(ik1),BS_blk_dim(ik2),(-1._SP,0._SP),& +& conjg(BS_cpl_mat),BS_blk_dim(ik1),Vi(frag_pointer(1)+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(2)+BS_K_dim+1),1) +#else + call cgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik2),(-1._SP,0._SP),& +& conjg(BS_mat),BS_blk_dim(ik1),Vi(frag_pointer(2)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+BS_K_dim+1),1) + call cgemv('c',BS_blk_dim(ik1),BS_blk_dim(ik2),(-1._SP,0._SP),& +& conjg(BS_mat),BS_blk_dim(ik1),Vi(frag_pointer(1)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(2)+BS_K_dim+1),1) + call cgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik2),(1._SP,0._SP),& +& BS_cpl_mat,BS_blk_dim(ik1),Vi(frag_pointer(2)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+1),1) + call cgemv('t',BS_blk_dim(ik1),BS_blk_dim(ik2),(1._SP,0._SP),& +& BS_cpl_mat,BS_blk_dim(ik1),Vi(frag_pointer(1)+BS_K_dim+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(2)+1),1) + call cgemv('n',BS_blk_dim(ik1),BS_blk_dim(ik2),(-1._SP,0._SP),& +& conjg(BS_cpl_mat),BS_blk_dim(ik1),Vi(frag_pointer(2)+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(1)+BS_K_dim+1),1) + call cgemv('t',BS_blk_dim(ik1),BS_blk_dim(ik2),(-1._SP,0._SP),& +& conjg(BS_cpl_mat),BS_blk_dim(ik1),Vi(frag_pointer(1)+1),1,(1._SP,0._SP),& +& Vo(frag_pointer(2)+BS_K_dim+1),1) +#endif + endif ! endif ! @@ -206,7 +281,23 @@ subroutine K_BSmat_by_V(iq,iter,Vi,Vo) enddo endif ! + ! Filling the anti-resonant and anti-coupling parts ! + if (BS_K_coupling.and..not.use_fragments.and.db_load) then + ! + ! If Coupling the half lower part of K must be filled + ! + ! Anti-resonant + ! + forall(i1=BS_K_dim+1:BS_H_dim,i2=BS_K_dim+1:BS_H_dim) & +& BS_mat(i1,i2)=-conjg(BS_mat(i1-BS_K_dim,i2-BS_K_dim)) + ! + ! Anti-coupling + ! + forall(i1=BS_K_dim+1:BS_H_dim,i2=1:BS_K_dim) & +& BS_mat(i1,i2)=-conjg(BS_mat(i1-BS_K_dim,i2+BS_K_dim)) + ! + endif ! ! No need to do any matrix X vector multiplication ! diff --git a/src/bse/K_Haydock.F b/src/bse/K_Haydock.F index 892bc1e332..0427d84522 100644 --- a/src/bse/K_Haydock.F +++ b/src/bse/K_Haydock.F @@ -69,6 +69,7 @@ subroutine K_Haydock(iq,W) ! if (BS_K_Coupling) then ! + allocate(Cf(Max_iterations+1),Vn(2*BS_K_dim),Vnm1(2*BS_K_dim),Vnp1(2*BS_K_dim)) ! else allocate(Vn(BS_K_dim),Vnm1(BS_K_dim),Vnp1(BS_K_dim)) @@ -82,6 +83,7 @@ subroutine K_Haydock(iq,W) ! if (BS_K_Coupling) then ! + io_err=ioBSS_Haydock(ID,1,2*BS_K_dim,reached_treshold,Af(1),Bf(1),Vnm1,Vn,Vnp1=Vnp1,Cf=Cf(1)) ! else io_err=ioBSS_Haydock(ID,1,BS_K_dim,reached_treshold,Af(1),Bf(1),Vnm1,Vn) @@ -93,8 +95,8 @@ subroutine K_Haydock(iq,W) it_on_disk=io_err call io_control(ACTION=OP_RD_CL,SEC=(/1,2/),ID=ID) if (BS_K_Coupling) then - ! - ! + io_err=ioBSS_Haydock(ID,it_on_disk,2*BS_K_dim,reached_treshold,Af(:it_on_disk),& +& Bf(:it_on_disk+1),Vnm1,Vn,Vnp1=Vnp1,Cf=Cf(:it_on_disk+1)) else io_err=ioBSS_Haydock(ID,it_on_disk,BS_K_dim,reached_treshold,Af(:it_on_disk),& & Bf(:it_on_disk+1),Vnm1,Vn) @@ -132,6 +134,25 @@ subroutine K_Haydock(iq,W) ! if (BS_K_Coupling) then ! + ! Before starting iterate one needs + ! + ! | Vn > = |q->0>/(0| F (H |q->0>)) + ! | Vnp1 > = H |Vn> + ! + ItParity = (-1._SP)**it_on_disk + if (it_on_disk==0) then + Cf= (0.,0.) + Vn(:BS_K_dim) = BSS_rhoq0(:BS_K_dim) *BS_eh_f(:BS_K_dim) + Vn(BS_K_dim+1:)=-BSS_rhoq0(BS_K_dim+1:)*BS_eh_f(:BS_K_dim) + call K_BSmat_by_V(iq,iter=it_on_disk,Vi=Vn,Vo=Vnp1) + Haydock_v0_mod=sqrt(fdot_product(Vn,Vnp1)) + Cf(1) = dot_product(BSS_rhoq0,Vn)/Haydock_v0_mod + Vnp1(:) = Vnp1(:)/Haydock_v0_mod + Vn(:) = Vn(:)/Haydock_v0_mod + else + Haydock_v0_mod = Bf(1) + Bf(1) = 0._SP + endif ! else ! @@ -184,14 +205,45 @@ subroutine K_Haydock(iq,W) ! else ! + ! A(n) = (that is ), + ! = 0 by symmetry + ! + ItParity=(-1._SP)**it + Af(it)=0._SP + ! + !|Vn+1> = |Vn+1> - A(n)|Vn> - B(n)|Vn-1> + ! + call V_by_V_plus_V(2*BS_K_dim,-Af(it)*(1._SP,0._SP),Vn,Vnp1) + call V_by_V_plus_V(2*BS_K_dim,-Bf(it)*(1._SP,0._SP),Vnm1,Vnp1) + ! + ! |Vn-1> = |Vn> + ! |Vn> = |Vn+1> + ! + Vnm1=Vn + Vn =Vnp1 + ! + !|Vn+1> = H |Vn> + ! + call K_BSmat_by_V(iq,iter=it-it_on_disk,Vi=Vn,Vo=Vnp1) + ! + ! B(n+1)= ^(1/2) (that is ^(1/2)) + ! = (2*Re())^(1/2) by symmetry, + ! where the dot_product is just on eh pair space + ! + Bf(it+1)=sqrt(fdot_product(Vn,Vnp1)) + ! + !|Vn> =|Vn+1> / B(n+1) + ! + forall(i1=1:2*BS_K_dim) Vn(i1)=Vn(i1)/Bf(it+1) + forall(i1=1:2*BS_K_dim) Vnp1(i1)=Vnp1(i1)/Bf(it+1) + Cf(it+1)= dot_product(BSS_rhoq0,Vn)/Haydock_v0_mod ! endif ! if (it>2) then ! if (BS_K_Coupling) then - ! - ! + call build_L_and_check_convergence(Af(:it),Bf(:it+1),it,Cf(:it)) else call build_L_and_check_convergence(Af(:it),Bf(:it+1),it) end if @@ -219,8 +271,10 @@ subroutine K_Haydock(iq,W) ! call io_control(ACTION=OP_WR_CL,SEC=(/1,2/),ID=ID) if (BS_K_Coupling) then - ! - ! + Bf(1) = Haydock_v0_mod + io_err=ioBSS_Haydock(ID,it,2*BS_K_dim,reached_treshold,Af(:it),Bf(:it+1), & +& Vnm1,Vn,Vnp1=Vnp1,Cf=Cf(:it+1)) + Bf(1) = 0._SP else io_err=ioBSS_Haydock(ID,it,BS_K_dim,reached_treshold,Af(:it),Bf(:it+1),Vnm1,Vn) end if @@ -237,8 +291,7 @@ subroutine K_Haydock(iq,W) ! CLEAN ! if (allocated(BS_mat)) deallocate(BS_mat) - ! - ! + if (BS_K_Coupling) deallocate(Cf) deallocate(Vn,Vnm1,Vnp1) ! ! Either if BS_mat is allocated (no fragments) or is not allocated @@ -276,9 +329,7 @@ subroutine build_L_and_check_convergence(Af,Bf,it,Cf) UseTerminator = .false. ! if (it > MIN_ITER.and.Haydock_terminator) UseTerminator = .true. - ! - ! - ! + if (BS_K_Coupling) allocate(X_i(2,it)) if (Useterminator) then if (.not.BS_K_Coupling) then Av1 = 0._SP @@ -290,8 +341,14 @@ subroutine build_L_and_check_convergence(Af,Bf,it,Cf) Av1 = Av1/it Av2 = Av2/it else - ! - ! + Av1 = 0._SP + Av2 = 0._SP + do i1 = 2,it+1,2 + Av1 = Av1 + Bf(i1) + Av2 = Av2 + Bf(i1+1) + end do + Av1 = Av1/(it/2+mod(it,2)) + Av2 = Av2/(it/2) end if end if ! @@ -299,8 +356,8 @@ subroutine build_L_and_check_convergence(Af,Bf,it,Cf) X_o=(0._SP,0._SP) if (UseTerminator) then if (.not.BS_K_Coupling) X_t = terminator(W%p(iw),Av1,Av2,RES) - ! - ! + if (BS_K_Coupling.and.(mod(it,2)==0)) X_t = terminator(W%p(iw),Av1,Av2,COUPL) + if (BS_K_Coupling.and.(mod(it,2)==1)) X_t = terminator(W%p(iw),Av2,Av1,COUPL) endif X_o(:,it)=1._SP/(W%p(iw)-Af(it)-Bf(it+1)**2*X_t) do i1=it-1,1,-1 @@ -309,8 +366,20 @@ subroutine build_L_and_check_convergence(Af,Bf,it,Cf) do i1=it-2,1,-1 X_o(2,i1)=1._SP/(W%p(iw)-Af(i1)-Bf(i1+1)**2*X_o(2,i1+1)) enddo - ! - ! + if (BS_K_Coupling) then + X_i(:,1) = X_o(:,1) + X_i(:,2) = Bf(2)*X_o(:,2)*X_i(:,1) + YbyX(:) = X_i(:,1)*Cf(1) + X_i(:,2)*Cf(2) + do i1 = 2,it-1 + X_i(1,i1+1) = -Bf(i1+1)*X_o(1,i1+1)*X_i(1,i1) + YbyX(1) = YbyX(1) + Cf(i1+1)*X_i(1,i1+1) + end do + do i1 = 2,it-2 + X_i(2,i1+1) = -Bf(i1+1)*X_o(2,i1+1)*X_i(2,i1) + YbyX(2) = YbyX(2) + Cf(i1+1)*X_i(2,i1+1) + end do + X_o(:,1) = YbyX(:) + endif X_epsilon(1,iw)=W%p(iw) X_epsilon(2,iw)=-X_o(1,1)*Co+1._SP X_epsilon(4,iw)=-X_o(2,1)*Co+1._SP @@ -319,8 +388,7 @@ subroutine build_L_and_check_convergence(Af,Bf,it,Cf) if (Haydock_treshold<0._SP) reached_treshold=reached_treshold+& & abs(X_epsilon(2,iw)-X_epsilon(4,iw))/abs(X_epsilon(2,iw))/real(W%n(1)) enddo - ! - ! + if (BS_K_Coupling) deallocate(X_i) end subroutine ! function fdot_product(v,w) @@ -342,8 +410,10 @@ function terminator(x,c1,c2,which) ! select case(which) case (1) - ! - ! + f = x**2 - c1**2 + c2**2 + g = 2._SP*x*c2**2 + terminator = f + sqrt(f**2 -2._SP*x*g) + terminator = terminator/g case (2) f = (x-c1) terminator = f + sqrt(f**2 -4._SP*c2**2.) diff --git a/src/bse/K_diagonalization.F b/src/bse/K_diagonalization.F index 350ab54c2e..279f21876d 100644 --- a/src/bse/K_diagonalization.F +++ b/src/bse/K_diagonalization.F @@ -293,10 +293,8 @@ subroutine K_diagonalization(iq,W) ! distinction between the resonant and antiresonant e/h Green's ! functions. ! - if (K_is_not_hermitian) X_epsilon(2,:)=& -& X_epsilon(2,:)-BS_R(i1)/(W%p(:)-BS_E(i1)) - if (.not.K_is_not_hermitian) X_epsilon(2,:)=& -& X_epsilon(2,:)-BS_R(i1)*conjg(BS_R(i1))/(W%p(:)-BS_E(i1)) + if (K_is_not_hermitian) X_epsilon(2,:)=X_epsilon(2,:)-BS_R(i1)/(W%p(:)-BS_E(i1)) + if (.not.K_is_not_hermitian) X_epsilon(2,:)=X_epsilon(2,:)-BS_R(i1)*conjg(BS_R(i1))/(W%p(:)-BS_E(i1)) ! call live_timing(steps=1) enddo diff --git a/src/bse/K_output_file.F b/src/bse/K_output_file.F index 239187bcb6..953130d004 100644 --- a/src/bse/K_output_file.F +++ b/src/bse/K_output_file.F @@ -23,7 +23,7 @@ subroutine K_output_file(iq,mode) ! use pars, ONLY:pi,SP use units, ONLY:HARTREE - use BS, ONLY:BSS_description, BSS_uses_GreenF,BSS_long_gauge,BSS_Vnl_included,& + use BS, ONLY:BSS_description, BSS_uses_GreenF,BSS_Vnl_included,& & BSS_n_descs,BS_K_coupling,BS_K_is_ALDA,BSS_n_freqs,BSS_q0 use X_m, ONLY:X_epsilon,eval_alpha,eps_2_alpha,O_eels use com, ONLY:msg,of_open_close @@ -76,7 +76,7 @@ subroutine K_output_file(iq,mode) endif ! call X_write_q_plus_g(1,BSS_q0,ig=1) - call X_write_messages_before_headers(iq,BSS_uses_GreenF,BSS_Vnl_included,BSS_long_gauge,ordering) + call X_write_messages_before_headers(iq,BSS_uses_GreenF,BSS_Vnl_included,ordering) ! ! Description of the run and headers for output files: ! diff --git a/src/bse/K_solvers.F b/src/bse/K_solvers.F index 9ee1430343..6cb7affce1 100644 --- a/src/bse/K_solvers.F +++ b/src/bse/K_solvers.F @@ -23,11 +23,11 @@ subroutine K_solvers(iq,Ken,Xk,q,X_static,W_bss) ! use pars, ONLY:SP,pi,schlen use memory_m, ONLY:mem_est - use X_m, ONLY:X_t,X_epsilon,X_alloc,DIP_iq_dot_r,& + use X_m, ONLY:X_t,X_epsilon,X_alloc,DIP_q_dot_iR,& & alpha_dim,eps_2_alpha,X_duplicate use BS, ONLY:BSS_n_freqs,BS_K_dim,BSS_rhoq0,BS_eh_f,BS_not_const_eh_f,& & BS_bands,BSS_q0,BS_eh_table,BSS_mode,BS_eh_E,& -& BSS_n_descs,BS_K_coupling,BSS_Vnl_included,BSS_long_gauge +& BSS_n_descs,BS_K_coupling,BSS_Vnl_included use com, ONLY:msg,warning use drivers, ONLY:l_bs_fxc use frequency, ONLY:w_samp @@ -83,7 +83,6 @@ subroutine K_solvers(iq,Ken,Xk,q,X_static,W_bss) X_oscillators%q0=BSS_q0 call Dipole_driver(Ken,Xk,X_oscillators,BSS_q0) BSS_Vnl_included=X_oscillators%Vnl_included - BSS_long_gauge=X_oscillators%long_gauge do i1=1,BS_K_dim ik =BS_eh_table(i1,1) iv =BS_eh_table(i1,2) @@ -94,11 +93,11 @@ subroutine K_solvers(iq,Ken,Xk,q,X_static,W_bss) ! ! iq . = - conjg( iq . ) ! - BSS_rhoq0(i1)=-conjg(DIP_iq_dot_r(ic,iv,ik,i_sp)) + BSS_rhoq0(i1)=-conjg(DIP_q_dot_iR(ic,iv,ik,i_sp)) ! ! minus comes from the occupation factor ! - if (BS_K_coupling) BSS_rhoq0(BS_K_dim+i1)=DIP_iq_dot_r(ic,iv,ik,i_sp) + if (BS_K_coupling) BSS_rhoq0(BS_K_dim+i1)=DIP_q_dot_iR(ic,iv,ik,i_sp) ! X_epsilon(3,:)=X_epsilon(3,:)-BSS_rhoq0(i1)*conjg(BSS_rhoq0(i1))*BS_eh_f(i1)/& & (W_bss%p(:)-BS_eh_E(i1)) @@ -118,8 +117,8 @@ subroutine K_solvers(iq,Ken,Xk,q,X_static,W_bss) ! ! Solvers ! - run_Haydock =index(BSS_mode,'h')/=0 run_Diago =index(BSS_mode,'d')/=0 + run_Haydock =index(BSS_mode,'h')/=0 ! ! if (run_Haydock.and.associated(Ken%W)) then @@ -145,6 +144,6 @@ subroutine K_solvers(iq,Ken,Xk,q,X_static,W_bss) ! deallocate(BSS_rhoq0,X_epsilon) call mem_est("BSS_rhoq0") - if (iq==1) call X_alloc('OptOsc') + if (iq==1) call X_alloc('DIP_q_dot_iR') ! end subroutine diff --git a/src/bz_ops/k_ibz2bz.F b/src/bz_ops/k_ibz2bz.F index 75a25c39d4..3267a02f23 100644 --- a/src/bz_ops/k_ibz2bz.F +++ b/src/bz_ops/k_ibz2bz.F @@ -59,6 +59,7 @@ subroutine k_ibz2bz(k,units,FORCE_BZ) ! ...then the pts ! allocate(k%ptbz(k%nbz,3)) + ! do i1=1,k%nbz ! ! k_bz is in iku diff --git a/src/common/.objects b/src/common/.objects index 1731a43a91..dae65e1a2f 100644 --- a/src/common/.objects +++ b/src/common/.objects @@ -1 +1 @@ -objs = eval_G_minus_G.o FREQUENCIES_damping.o FREQUENCIES_setup.o FREQUENCIES_coarse_grid.o crystal_lattice.o Kramers_Kronig.o OCCUPATIONS_Fermi.o OCCUPATIONS_Extend.o QP_state_extract.o pol_fit.o +objs = eval_G_minus_G.o FREQUENCIES_damping.o FREQUENCIES_setup.o FREQUENCIES_coarse_grid.o FREQUENCIES_Green_Function.o FREQUENCIES_reset.o crystal_lattice.o Kramers_Kronig.o Integrate.o QPartilize.o OCCUPATIONS_Fermi.o OCCUPATIONS_Extend.o QP_state_extract.o pol_fit.o diff --git a/src/common/FREQUENCIES_Green_Function.F b/src/common/FREQUENCIES_Green_Function.F new file mode 100644 index 0000000000..b8fccbc9c0 --- /dev/null +++ b/src/common/FREQUENCIES_Green_Function.F @@ -0,0 +1,43 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine FREQUENCIES_Green_Function(iqp,W,E) + ! + use pars, ONLY:SP + use R_lattice, ONLY:nkibz + use electrons, ONLY:n_bands,spin,n_sp_pol + use frequency, ONLY:w_samp + use QP_m, ONLY:QP_table,QP_G_dr,QP_G_er,QP_G_zoom_er + implicit none + ! + integer :: iqp + type(w_samp) :: W + real(SP) :: E(n_bands,nkibz,n_sp_pol) + ! + if (allocated(QP_G_zoom_er)) then + W%er=QP_G_zoom_er(iqp,:) + else + W%er=QP_G_er+E(QP_table(iqp,1),QP_table(iqp,3),spin(QP_table(iqp,:))) + endif + W%dr=QP_G_dr + call FREQUENCIES_setup(W) + ! +end subroutine diff --git a/src/common/FREQUENCIES_reset.F b/src/common/FREQUENCIES_reset.F new file mode 100644 index 0000000000..1f41c911b0 --- /dev/null +++ b/src/common/FREQUENCIES_reset.F @@ -0,0 +1,46 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine FREQUENCIES_reset(Xw) + ! + ! Input + !------- + ! integer :: npts + ! real(SP) :: bg_pt(npts),cg_percentual + ! + ! Deallocate: + ! bg_pt(:), cg_pt(:), cg_index_bg(:), + ! X_poles_tab(:), rg_index_bg(:), bg_npts(:) + ! + use frequency, ONLY:w_samp,bg_npts,cg_npts,cg_pt,rg_index_bg,cg_index_bg + use memory_m , ONLY:mem_est + use X_m, ONLY:X_poles_tab + implicit none + ! + type(w_samp) :: Xw + ! + if(associated(Xw%p)) deallocate(Xw%p) + call mem_est("W-p") + if (allocated(rg_index_bg)) deallocate(rg_index_bg) + deallocate(X_poles_tab,bg_npts,cg_pt,cg_index_bg) + call mem_est("X_poles_tab RGi BGn CGp CGi") + ! +end subroutine diff --git a/src/common/Integrate.F b/src/common/Integrate.F new file mode 100644 index 0000000000..e0901534f4 --- /dev/null +++ b/src/common/Integrate.F @@ -0,0 +1,43 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +function Integrate(F,W,N) + ! + use pars, ONLY:SP + implicit none + integer, intent(in) :: N + real(SP), intent(in) :: W(N),F(N) + real(SP) :: Integrate + ! + ! Work Space + ! + integer :: i_w + real(SP) :: M,Q + ! + Integrate=0._SP + ! + do i_w=1,N-1 + M=(F(i_w)-F(i_w+1))/(W(i_w)-W(i_w+1)) + Q=F(i_w)-M*W(i_w) + Integrate=Integrate+1./2.*M*(W(i_w+1)**2.-W(i_w)**2.)+Q*(W(i_w+1)-W(i_w)) + enddo + ! +end function diff --git a/src/common/QPartilize.F b/src/common/QPartilize.F new file mode 100644 index 0000000000..15c5f888ec --- /dev/null +++ b/src/common/QPartilize.F @@ -0,0 +1,65 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine QPartilize(N,G,W,E,Z,dG_step) + ! + use pars, ONLY:SP + implicit none + integer, intent(in) :: N + real(SP), intent(inout):: dG_step + real(SP), intent(in) :: W(N) + complex(SP), intent(in) :: G(N) + complex(SP), intent(out) :: E,Z + ! + ! Work Space + ! + integer ::i_w,n_poles,dG_steps,ip_larger_R,ip + integer,parameter::n_max_poles=100 + real(SP) ::R(n_max_poles),Eqp(n_max_poles),Width(n_max_poles) + complex(SP) ::beta(n_max_poles),G_m1(N),A(N) + ! + G_m1(:)=1./G(:) + A(:) =W(:)-G_m1(:) + ! + if (dG_step>0.) dG_steps=int(dG_step/real(W(2)-W(1))) + ! + n_poles=0 + do i_w=1,N-1 + if (real(G_m1(i_w))*real(G_m1(i_w+1))<0.) then + n_poles=n_poles+1 + Eqp(n_poles) =real(W(i_w)) + Width(n_poles) =aimag(A(i_w)) + R(n_poles) =abs(aimag(G(i_w))) + beta(n_poles)=(A(min(i_w+dG_steps,N))-A(i_w))/real(W(min(i_w+dG_steps,N))-W(i_w)) + endif + if (n_poles==n_max_poles) exit + enddo + ip_larger_R=1 + if (n_poles>1) then + do ip=2,n_poles + if (R(ip)>R(ip_larger_R)) ip_larger_R=ip + enddo + endif + ! + Z=1./(1.-beta(ip_larger_R)) + E=Eqp(ip_larger_R)+Z*CMPLX(0.,Width(ip_larger_R)) + ! +end subroutine diff --git a/src/interface/barriers.F b/src/interface/barriers.F index 679ce7d19f..6e2780e940 100644 --- a/src/interface/barriers.F +++ b/src/interface/barriers.F @@ -89,7 +89,6 @@ subroutine barriers() ! (those runlevels can be on if a non-gpl input file is being reading) !==================================================================== ! - call switch_off_runlevel('life',on_name=' ') call switch_off_runlevel('acfdt',on_name=' ') call switch_off_runlevel('bs_fxc',on_name=' ') ! @@ -104,7 +103,17 @@ subroutine barriers() endif ! ! + if (l_life) then + call switch_off_runlevel('all',on_name='life em1d '//trim(always_runlevels)) + goto 1 + endif ! + if (gw0_raxis) then + on_runlevels='gw0 em1d el_el el_ph HF_and_locXC' + on_runlevels='gw0 em1d el_el HF_and_locXC' + call switch_off_runlevel('all',on_name=trim(on_runlevels)//' '//trim(always_runlevels)) + goto 1 + endif ! if (gw0_ppa) then on_runlevels='gw0 ppa em1d el_el el_ph HF_and_locXC' @@ -122,7 +131,7 @@ subroutine barriers() endif if (l_bss.or.l_bse) then on_runlevels='tdhf bs_fxc bse bss em1s em1d ppa optics' - on_runlevels='tdhf bse bss em1s ppa optics' + on_runlevels='tdhf bse bss em1s em1d ppa optics' call switch_off_runlevel('all',on_name=trim(on_runlevels)//' '//trim(always_runlevels)) goto 1 endif diff --git a/src/interface/init.F b/src/interface/init.F index a2f80c59af..99f7e9c38e 100644 --- a/src/interface/init.F +++ b/src/interface/init.F @@ -63,11 +63,12 @@ integer function init(en,q,k,X,Xw,instr,lnstr,CLOSE_Gs,FINALIZE) ! ! Work Space ! - integer :: io_err,ioWF_err,io_X_err(4),io_ID,ioBS_err,& + integer :: io_err,ioWF_err,io_X_err(4),io_ID,ioBS_err,io_KB_PP_err,& & ioBS_Fxc_err,ioQINDX_err,io_SC_E_err,io_SC_V_err,io_COLLISIONS_err - integer, external :: ioX,ioOSTNTS,ioGROT,ioQINDX,ioRIM,& + integer, external :: ioX,io_DIPOLES,ioGROT,ioQINDX,ioRIM,& & io_HF_and_locXC,ioQP,ioBS,ioDB1,ioKB_PP,& & ioCOL_CUT + logical :: OSTNTS_Vnl_included ! ! ! @@ -225,20 +226,29 @@ integer function init(en,q,k,X,Xw,instr,lnstr,CLOSE_Gs,FINALIZE) ! if (associated(qp%table)) nullify(qp%table) ! + ! Green Functions + ! + call io_control(ACTION=OP_RD_CL,COM=NONE,SEC=(/1,2/),MODE=DUMP,ID=io_ID) + io_err=ioQP('G',qp,io_ID) + ! + ! W + ! + call io_control(ACTION=OP_RD_CL,COM=NONE,SEC=(/1,2/),MODE=DUMP,ID=io_ID) + io_err=ioQP('W',qp,io_ID) ! ! ostnts ! call io_control(ACTION=OP_RD_CL,COM=NONE,SEC=(/1/),MODE=DUMP,ID=io_ID) - io_err=ioOSTNTS(X(3),en,io_ID) + io_err=io_DIPOLES(X(3),en,io_ID) + ! + OSTNTS_Vnl_included=io_err==0.and.X(3)%Vnl_included ! ! kb_pp ! call io_control(ACTION=OP_RD_CL,COM=NONE,SEC=(/1/),MODE=DUMP,ID=io_ID) - io_err=ioKB_PP(io_ID) - ! - ! If the KB non local factors are found set to T the Vnl_included logical + io_KB_PP_err=ioKB_PP(io_ID) ! - if (io_err==0) forall(i1=1:4) X(i1)%Vnl_included=.true. + if (io_err/=0) OSTNTS_Vnl_included=io_KB_PP_err==0 ! ! I transfer to all X types the X(3) used in the previous io's ! @@ -251,6 +261,13 @@ integer function init(en,q,k,X,Xw,instr,lnstr,CLOSE_Gs,FINALIZE) if (io_X_err(i1)>0) X(i1)%iq(1)=io_X_err(i1)+1 enddo ! + ! The GLOBAL vcalue of %Vnl_included is decided on the basis of the contents + ! of db.OSTENTS OR on the presence of the KB_PP. This means that if the + ! response functions DBs were made in presence of KB_PP and later this + ! DB is deleted the X dbs will be recalculated + ! + forall(i1=1:4) X(i1)%Vnl_included=OSTNTS_Vnl_included + ! ! bs ! call io_control(ACTION=OP_RD_CL,COM=NONE,MODE=DUMP,SEC=(/1/),ID=io_ID) @@ -453,14 +470,16 @@ integer function init(en,q,k,X,Xw,instr,lnstr,CLOSE_Gs,FINALIZE) call initactivate(1,'BoseCut LongPath') call initactivate(1,'BSSmod BEnRange BDmRange BDmERef BEnSteps BLongDir') if (index(BSS_mode,'d')/=0) call initactivate(1,'WRbsWF') - if (index(BSS_mode,'h')/=0) call initactivate(1,'BSHayTrs') - if (index(BSS_mode,'h')/=0) call initactivate(1,'BSHayTer') + if (index(BSS_mode,'h')/=0) call initactivate(1,'BSHayTrs BSHayTer') if (index(BSS_mode,'t')/=0) call initactivate(1,'FxcGRLc FxcSVdig FxcMEStps FxcCausal') if (index(BSS_mode,'i')/=0) call initactivate(1,'SkipFullInv') endif ! ! DFT ! +#if defined _TWO_LEVELS + call initactivate(1,'RT3Levels RTEnLevels RTDephasing RTLifeTimes RTDipoles') +#endif ! ! ! RAS @@ -579,6 +598,7 @@ subroutine logicalson l_alda_fxc=runlevel_is_on('alda_fxc') l_lrc_fxc=runlevel_is_on('lrc_fxc') l_td_hf=runlevel_is_on('tdhf') + l_life=runlevel_is_on('life') ! ! l_cohsex=runlevel_is_on('cohsex') @@ -681,7 +701,8 @@ end subroutine varsetup2 subroutine read_command_line(rstr,init_) ! use stderr, ONLY:string_split - use it_m, ONLY:V_RL,V_kpt,V_sc,V_qp,V_io,V_general,V_resp,V_real_time + use it_m, ONLY:V_RL,V_kpt,V_sc,V_qp,V_io,V_general,V_resp,& +& V_real_time,V_all implicit none integer :: init_ character(*):: rstr @@ -733,6 +754,7 @@ subroutine read_command_line(rstr,init_) ! V_general=6 ! V_resp=7 ! V_real_time=8 + ! V_all=99 ! if ( trim(rstr_piece(i1)) == 'infver' ) then select case (trim(rstr_piece(i1+1))) @@ -752,6 +774,8 @@ subroutine read_command_line(rstr,init_) infile_verbosity=V_resp case ('rt') infile_verbosity=V_real_time + case ('all') + infile_verbosity=V_all end select endif ! diff --git a/src/interface/init_load.F b/src/interface/init_load.F index 6302a27876..88bd26f311 100644 --- a/src/interface/init_load.F +++ b/src/interface/init_load.F @@ -24,10 +24,11 @@ subroutine init_load(defs,en,q,k,X,Xw) use pars, ONLY:schlen use electrons, ONLY:levels,nel,filled_tresh use frequency, ONLY:w_samp - use it_m, ONLY:it,initdefs,E_unit,G_unit,T_unit,Bfield_unit,Efield_unit,& + use it_m, ONLY:it,initdefs,E_unit,G_unit,T_unit,Bfield_unit,& & Time_unit,I_unit,Angle_unit,& & V_RL,V_kpt,V_sc,V_qp,V_io,V_general,V_resp,V_real_time - use X_m, ONLY:X_t,long_path,q_plus_G_direction,Q_Shift_Order + use X_m, ONLY:X_t,q_plus_G_direction,Q_Shift_Order + use fields, ONLY:global_gauge,grid_path use QP_m, ONLY:QP_cg_percent,QP_G_damp,QP_solver,& & QP_n_G_bands,QP_ng_Sx,& & QP_G_er,QP_G_dr,QP_Sc_steps,GWo_iterations,SC_band_mixing,& @@ -46,6 +47,9 @@ subroutine init_load(defs,en,q,k,X,Xw) use ACFDT, ONLY:ACFDT_n_lambda,ACFDT_n_freqs,ACFDT_E_range use com, ONLY:more_io_path,com_path use functions, ONLY:bose_E_cut +#if defined _TWO_LEVELS + use real_time, ONLY:l_3levels,en_levels,lambda_diag,lambda_offdiag,dip +#endif ! implicit none type(initdefs)::defs @@ -72,6 +76,7 @@ subroutine init_load(defs,en,q,k,X,Xw) call it('r',defs,'rim_cut', '[R RIM CUT] Coulomb interaction') call it('r',defs,'tdhf', '[R TDHF] Unscreened e-h interaction: TD Hartree-Fock') call it('r',defs,'cohsex', '[R Xp] COlumb Hole Screened EXchange') + call it('r',defs,'life', '[R GW] GoWo Quasiparticle lifetimes') ! ! ! @@ -88,20 +93,27 @@ subroutine init_load(defs,en,q,k,X,Xw) call it(defs,'NonPDirs','[X/BSS] Non periodic chartesian directions (X,Y,Z,XY...)',non_periodic_directions,verb_level=V_resp) call it(defs,'IkSigLim','[KPT] QP K-points indices range',QP_states_k,verb_level=V_kpt) call it(defs,'IkXLim', '[KPT] X grid last k-point index',nXkibz,verb_level=V_kpt) + call it(defs,'Nelectro','Electrons number',nel,verb_level=V_general) ! ! S_xc ! + call it(defs,'LifeTrCG', '[GW] [o/o] Lifetime transitions reduction',QP_cg_percent) call it(defs,'EXXRLvcs', '[XX] Exchange RL components',QP_ng_Sx,G_unit) call it(defs,'GbndRnge', '[GW] G[W] bands range',QP_n_G_bands) call it(defs,'GDamping', '[GW] G[W] damping',QP_G_damp,E_unit) call it(defs,'GDmRnge', '[GW] G_gw damping range',QP_G_dr,E_unit) call it(defs,'dScStep', '[GW] Energy step to evalute Z factors',QP_dSc_delta,E_unit) - ! + call it(defs,'DysSolver','[GW] Dyson Equation solver (`n`,`s`,`g`)',QP_solver,protect=.FALSE.) + call it(defs,'GEnSteps', '[GW] G_gw energy steps',QP_Sc_steps) + call it(defs,'GEnRnge', '[GW] G_gw energy range (centered in the bare energy)',QP_G_er,E_unit) + call it('f',defs,'NewtDchk', '[F GW] Test dSc/dw convergence',verb_level=V_qp) + call it('f',defs,'ExtendOut', '[F GW] Print all variables in the output file',verb_level=V_qp) ! ! ! BSE ! call it(defs,'BSresKmod', '[BSK] Resonant Kernel mode. (`x`;`c`;`d`)',BS_res_mode) + call it(defs,'BScplKmod', '[BSK] Coupling Kernel mode. (`x`;`c`;`d`;`u`)',BS_cpl_mode) call it(defs,'BSEBands','[BSK] Bands range',BS_bands) call it(defs,'BSENGBlk','[BSK] Screened interaction block size',BS_n_g_W,G_unit) call it(defs,'BSENGexx','[BSK] Exchange components',BS_n_g_exch,G_unit) diff --git a/src/io/.objects b/src/io/.objects index 09a8d4b7ff..c181de7afe 100644 --- a/src/io/.objects +++ b/src/io/.objects @@ -1 +1 @@ -objs = ver_is_gt_or_eq.o IO_and_Messaging_switch.o ioRIM.o ioGROT.o ioOSTNTS.o ioQINDX.o ioQP.o ioX.o ioDB1.o ioWF.o ioKB_PP.o io_HF_and_locXC.o ioBS.o ioBSS_diago.o ioBSS_Haydock.o io_elemental.o io_bulk.o io_header.o ioE_RIM.o ioCOL_CUT.o ioDB1_reload.o variables_BS.o IO_Setup.o Fragments_Restart.o Fragments_Synchronize.o +objs = ver_is_gt_or_eq.o IO_and_Messaging_switch.o ioRIM.o ioGROT.o io_DIPOLES.o ioQINDX.o ioQP.o ioX.o ioDB1.o ioWF.o ioKB_PP.o io_HF_and_locXC.o ioBS.o ioBSS_diago.o ioBSS_Haydock.o io_elemental.o io_bulk.o io_header.o ioE_RIM.o ioCOL_CUT.o ioDB1_reload.o variables_BS.o IO_Setup.o Fragments_Restart.o Fragments_Synchronize.o diff --git a/src/io/ioBSS_Haydock.F b/src/io/ioBSS_Haydock.F index 080cb69ca6..7f2f2f8441 100644 --- a/src/io/ioBSS_Haydock.F +++ b/src/io/ioBSS_Haydock.F @@ -43,7 +43,7 @@ integer function ioBSS_Haydock(ID,it,BS_H_dim,reached_treshold,Af,Bf,Vnm1,Vn,Vnp ! if (any((/io_sec(ID,:)==1/))) then ! - ioBSS_Haydock=io_header(ID,XC_KIND="K_E force") + ioBSS_Haydock=io_header(ID,XC_KIND="K_E force",GAUGE=.TRUE.,IMPOSE_GAUGE=.TRUE.) ! call io_elemental(ID,VAR="PARS",VAR_SZ=6,MENU=0) call io_elemental(ID,VAR=& diff --git a/src/io/ioDB1.F b/src/io/ioDB1.F index ef7feae5b0..666a00925f 100644 --- a/src/io/ioDB1.F +++ b/src/io/ioDB1.F @@ -40,7 +40,7 @@ integer function ioDB1(E,k,ID) & DUMP,io_mode,write_is_on,io_code_version,ver_is_gt_or_eq use QP_m, ONLY:QP_nk,QP_nb,QP_n_G_bands use BS, ONLY:BS_bands - use xc_functionals, ONLY:GS_xc_KIND,GS_xc_FUNCTIONAL,XC_USED_LIBXC + use xc_functionals, ONLY:GS_xc_KIND,GS_xc_FUNCTIONAL use global_XC, ONLY:setup_global_XC,EXT_NONE ! implicit none @@ -152,14 +152,6 @@ integer function ioDB1(E,k,ID) & VAR=' WF G-vectors :') endif ! - ! Check (n)s.db1 was generated with a a2y linked or not to libxc - ! If this is the case and yambo reading it is not, the program should stop - ! and viceversa if a2y was not linked to libxc and yambo it is, the program should stop - ! - l_use_libxc = .false. - if (l_use_libxc.and.(GS_xc_KIND.ne.XC_USED_LIBXC)) & - &call error("ns.db1 not compatible with yambo executable [LIBXC]") - ! ! Atomic data ! if (ver_is_gt_or_eq(ID,(/3,0,4/))) then diff --git a/src/io/ioOSTNTS.F b/src/io/ioOSTNTS.F index f4111f390b..e69de29bb2 100644 --- a/src/io/ioOSTNTS.F +++ b/src/io/ioOSTNTS.F @@ -1,144 +0,0 @@ -! -! Copyright (C) 2000-2010 A. Marini and the YAMBO team -! http://www.yambo-code.org -! -! This file is distributed under the terms of the GNU -! General Public License. You can redistribute it and/or -! modify it under the terms of the GNU General Public -! License as published by the Free Software Foundation; -! either version 2, or (at your option) any later version. -! -! This program is distributed in the hope that it will -! be useful, but WITHOUT ANY WARRANTY; without even the -! implied warranty of MERCHANTABILITY or FITNESS FOR A -! PARTICULAR PURPOSE. See the GNU General Public License -! for more details. -! -! You should have received a copy of the GNU General Public -! License along with this program; if not, write to the Free -! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, -! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. -! -integer function ioOSTNTS(X,Xen,ID) - ! - use pars, ONLY:SP,schlen - use units, ONLY:HARTREE - use X_m, ONLY:X_t,DIP_iR_or_P,X_alloc,Dipole_bands_ordered - use electrons, ONLY:levels,n_sp_pol - use R_lattice, ONLY:nXkibz,q0_def_norm - use matrix_operate,ONLY:mat_c2r,mat_r2c - use IO_m, ONLY:io_connect,io_disconnect,io_sec,& -& io_elemental,io_status,io_bulk,& -& read_is_on,write_is_on,io_header,ver_is_gt_or_eq - use global_XC, ONLY:Dipole_WF_xc_string,loaded_WF_xc_string - ! - implicit none - type(X_t)::X - type(levels)::Xen - integer ::ID - ! - ! Work Space - ! - integer :: i1,ixyz,sec_size,i_spin - integer :: db_nbm,db_nbf,db_nb(2) - character(schlen) :: VAR_name - real(SP),allocatable :: DIP_iR_or_P_disk(:,:,:) - ! - ioOSTNTS=io_connect(desc='ostnts',type=2,ID=ID) - if (ioOSTNTS/=0) goto 1 - ! - sec_size=8 - ! - ! In V 3.0.8 q0_def_norm is stored in case a longitudinal - ! gauge grid has already redefined the q0 norm - ! - if (ver_is_gt_or_eq(ID,(/3,0,9/))) sec_size=10 - if (ver_is_gt_or_eq(ID,(/3,0,15/))) sec_size=11 - ! - if (any((/io_sec(ID,:)==1/))) then - ! - ioOSTNTS=io_header(ID,R_LATT=.true.,WF=.true.,IMPOSE_SN=.true.,T_EL=.true.) - if (ioOSTNTS/=0) goto 1 - ! - call io_elemental(ID,VAR="PARS",VAR_SZ=sec_size,MENU=0) - call io_elemental(ID,DB_I1=db_nb,& -& VAR=" X band range :",I1=X%ib,CHECK=.true.,OP=(/">=","<="/)) - call io_elemental(ID,UNIT=HARTREE,& -& VAR=" X e/h energy range [ev]:",R1=X%ehe,CHECK=.true.,OP=(/">=","<="/)) - call io_elemental(ID,DB_I0=db_nbm,& -& VAR=" Metallic bands :",I0=Xen%nbm,CHECK=.true.,OP=(/"=="/)) - call io_elemental(ID,DB_I0=db_nbf,& -& VAR=" Filled bands :",I0=Xen%nbf,CHECK=.true.,OP=(/"=="/)) - call io_elemental(ID,& -& VAR=" RL vectors in the sum :",I0=X%ngostnts,WARN=.true.,OP=(/"<="/)) - call io_elemental(ID,& -& VAR=" [r,Vnl] included :",L0=X%Vnl_included,CHECK=.true.,OP=(/"=="/)) - ! - if (ver_is_gt_or_eq(ID,(/3,0,15/))) then - call io_elemental(ID,& -& VAR=" Transitions ordered :",L0=Dipole_bands_ordered,CHECK=.true.,OP=(/"=="/)) - endif - ! - if (ver_is_gt_or_eq(ID,(/3,0,9/))) then - call io_elemental(ID,& -& VAR=" Longitudinal Gauge :",L0=X%long_gauge,CHECK=.true.,OP=(/"=="/)) - call io_elemental(ID,& -& VAR=" Field momentum norm :",R0=q0_def_norm,CHECK=.true.,OP=(/"=="/)) - endif - ! - if (ver_is_gt_or_eq(ID,(/3,0,15/))) then - ! - call io_elemental(ID,VAR="",VAR_SZ=0,MENU=0) - ! - ! Wavefunctions xc - ! - call io_elemental(ID,VAR='WAVE_FUNC_XC',CH0="",VAR_SZ=1,MENU=0) - call io_elemental(ID,DB_CH0=Dipole_WF_xc_string,CH0=loaded_WF_xc_string,& -& VAR=' Wavefunctions :',CHECK=.true.,OP=(/"=="/)) - call io_elemental(ID,VAR="",VAR_SZ=0) - else - ! - call io_elemental(ID,VAR="",VAR_SZ=0) - ! - endif - ! - ioOSTNTS=io_status(ID) - if (ioOSTNTS/=0) goto 1 - endif - ! - ! On disk the size is DIP_iR_or_P(3,db_nb(2),db_nbm,nXkibz) - ! - sec_size=3*db_nb(2)*db_nbm - if (any((/io_sec(ID,:)==2/))) then - ! - allocate(DIP_iR_or_P_disk(db_nb(2),db_nbm,2)) - ! - if (read_is_on(ID)) call X_alloc('OptGrad',(/3,X%ib(2),Xen%nbm,nXkibz/)) - ! - do i1=1,nXkibz - ! - do ixyz=1,3 - ! - do i_spin=1,n_sp_pol - ! - if (write_is_on(ID)) call mat_c2r(DIP_iR_or_P(ixyz,:,:,i1,i_spin),DIP_iR_or_P_disk) - ! - write (VAR_name,'(3(a,i4.4))') 'OptGrad_k_',i1,'_xyz_',ixyz,'_spin_',i_spin - call io_bulk(ID,VAR=trim(VAR_name),VAR_SZ=shape(DIP_iR_or_P_disk)) - call io_bulk(ID,R3=DIP_iR_or_P_disk) - ! - if (read_is_on(ID)) call mat_r2c(DIP_iR_or_P_disk,DIP_iR_or_P(ixyz,:,:,i1,i_spin)) - ! - enddo - ! - enddo - ! - enddo - ! - deallocate(DIP_iR_or_P_disk) - ! - endif - ! -1 call io_disconnect(ID=ID) - ! -end function diff --git a/src/io/ioQP.F b/src/io/ioQP.F index 7029b287aa..86aee05041 100644 --- a/src/io/ioQP.F +++ b/src/io/ioQP.F @@ -74,18 +74,26 @@ integer function ioQP(what,qp,ID) ! reload to apply QP corrections. In this case QP_DB_kind must be ! present to be transferred in the energy type by mod_qp_ctl ! +#if defined _NETCDF_IO if (variable_is_found(ID,"QP_DB_kind")==1.or.write_is_on(ID)) then +#endif call io_elemental(ID,VAR="QP_DB_kind",VAR_SZ=1,MENU=0) call io_elemental(ID,I0=QP_DB_kind) call io_elemental(ID,VAR="",VAR_SZ=0,MENU=0) +#if defined _NETCDF_IO else if (variable_is_found(ID,"QP_DB_kind")< 0.and.read_is_on(ID)) then QP_DB_kind=0 ioQP=IO_INCOMPATIBLE_VAR return endif +#endif + ! endif ! l_CUTOFF=.FALSE. +#if defined_TEST_2 + l_CUTOFF=.TRUE. +#endif ! if (what=="W") then ioQP=io_header(ID,QPTS=.true.,R_LATT=.true.,WF=.true.,T_EL=.true.,D_LATT=.true.,XC_KIND="Xd",CUTOFF=l_CUTOFF) @@ -131,6 +139,23 @@ integer function ioQP(what,qp,ID) ! call io_elemental(ID,VAR="",VAR_SZ=0,MENU=0) ! + if (index(what,".G")/=0.or.what=="G") then + ! + ! Tough the GF's parameters are already stored in the QP_descriptions + ! I need then on-fly to, eventually, rebuild the Green's function + ! + var_size=6 + if (ver_is_gt_or_eq(ID,revision=530)) var_size=7 + call io_elemental(ID,VAR="SE_OPERATOR_PARAMETERS",VAR_SZ=var_size,MENU=0) + call io_elemental(ID,I0=qp%GreenF_n_steps) + call io_elemental(ID,R1=QP_G_er) ! <- This is, actually, not needed as the + ! full frequency dependence is stored in sec 3 + call io_elemental(ID,R1=QP_G_dr) + call io_elemental(ID,L0=GF_is_causal) + if (ver_is_gt_or_eq(ID,revision=530)) call io_elemental(ID,R0=QP_G_Zoom_treshold) + call io_elemental(ID,VAR="",VAR_SZ=0,MENU=0) + ! + endif ! ioQP=io_status(ID) if (.not.DB_is_OK(ID)) goto 1 @@ -172,7 +197,6 @@ integer function ioQP(what,qp,ID) endif else if (what=="G") then call io_elemental(ID,CH0=qp%description(i1),VAR='',CHECK=.true.,OP=(/"=="/)) - if (index(qp%description(i1),'Zoom')>0.and.io_status(ID)<0) io_status(ID)=IO_NO_BINDING_ERROR else if (index(qp%description(i1),'QP')>0) then call io_elemental(ID,CH0=qp%description(i1),VAR='') @@ -304,6 +328,51 @@ integer function ioQP(what,qp,ID) ! endif ! + if (index(what,".G")/=0.or.what=="G") then + ! + allocate(qp_DATA_disk(qp%n_states,qp%GreenF_n_steps,6)) + ! + ! Real axis Self Energy & Green Function + !========================================= + ! + if (write_is_on(ID)) then + do i1=1,qp%n_states + qp_DATA_disk(i1,:,1)=real(qp%S_total(i1,:)) + qp_DATA_disk(i1,:,2)=aimag(qp%S_total(i1,:)) + qp_DATA_disk(i1,:,3)=real(qp%GreenF(i1,:)) + qp_DATA_disk(i1,:,4)=aimag(qp%GreenF(i1,:)) + qp_DATA_disk(i1,:,5)=real(qp%GreenF_W(i1,:)) + qp_DATA_disk(i1,:,6)=aimag(qp%GreenF_W(i1,:)) + enddo + endif + ! + call io_bulk(ID,VAR="SE_Operator",VAR_SZ=(/qp%n_states,qp%GreenF_n_steps,2/)) + call io_bulk(ID,R3=qp_DATA_disk(:,:,1:2)) + call io_bulk(ID,VAR="Green_Functions",VAR_SZ=(/qp%n_states,qp%GreenF_n_steps,2/)) + call io_bulk(ID,R3=qp_DATA_disk(:,:,3:4)) + call io_bulk(ID,VAR="Green_Functions_Energies",VAR_SZ=(/qp%n_states,qp%GreenF_n_steps,2/)) + call io_bulk(ID,R3=qp_DATA_disk(:,:,5:6)) + ! + if (read_is_on(ID)) then + allocate(qp%S_total(qp%n_states,qp%GreenF_n_steps)) + call mem_est("qp_S_total",(/size(qp%S_total)/)) + allocate(qp%GreenF(qp%n_states,qp%GreenF_n_steps)) + call mem_est("qp_GreenF",(/size(qp%GreenF)/)) + allocate(qp%GreenF_W(qp%n_states,qp%GreenF_n_steps)) + call mem_est("qp_GreenF_W",(/size(qp%GreenF_W)/)) + do i1=1,qp%n_states + qp%S_total(i1,:) =cmplx(qp_DATA_disk(i1,:,1),qp_DATA_disk(i1,:,2),SP) + qp%GreenF(i1,:) =cmplx(qp_DATA_disk(i1,:,3),qp_DATA_disk(i1,:,4),SP) + qp%GreenF_W(i1,:)=cmplx(qp_DATA_disk(i1,:,5),qp_DATA_disk(i1,:,6),SP) + enddo + endif + ! + deallocate(qp_DATA_disk) + ! + endif + ! + ioQP=io_status(ID) + if (.not.DB_is_OK(ID)) goto 1 ! endif ! diff --git a/src/io/ioX.F b/src/io/ioX.F index 4707e07451..c8c598d623 100644 --- a/src/io/ioX.F +++ b/src/io/ioX.F @@ -44,8 +44,9 @@ integer function ioX(X,Xw,ID) !Work Space ! integer ::sec_size,i1,iq,i_err,io_com_save - logical ::different_db_RL_order,l_CUTOFF + logical ::different_db_RL_order,l_CUTOFF,local_long_gauge character (schlen) ::ch,db_desc + character (2) ::local_XC_KIND real(SP) ::Wd_disk(2) real(SP),allocatable::X_disk(:,:,:),W_disk(:,:),RL_vecs_disk(:,:) ! @@ -72,10 +73,11 @@ integer function ioX(X,Xw,ID) ! l_CUTOFF=.FALSE. ! - if (X%whoami==1) ioX=io_header(ID,QPTS=.true.,R_LATT=.true.,WF=.true.,T_EL=.true.,XC_KIND="Xx",CUTOFF=l_CUTOFF) - if (X%whoami==2) ioX=io_header(ID,QPTS=.true.,R_LATT=.true.,WF=.true.,T_EL=.true.,XC_KIND="Xs",CUTOFF=l_CUTOFF) - if (X%whoami==3) ioX=io_header(ID,QPTS=.true.,R_LATT=.true.,WF=.true.,T_EL=.true.,XC_KIND="Xd",CUTOFF=l_CUTOFF) - if (X%whoami==4) ioX=io_header(ID,QPTS=.true.,R_LATT=.true.,WF=.true.,T_EL=.true.,XC_KIND="Xp",CUTOFF=l_CUTOFF) + if (X%whoami==1) local_XC_KIND="Xx" + if (X%whoami==2) local_XC_KIND="Xs" + if (X%whoami==3) local_XC_KIND="Xd" + if (X%whoami==4) local_XC_KIND="Xp" + ioX=io_header(ID,QPTS=.true.,R_LATT=.true.,WF=.true.,T_EL=.true.,XC_KIND=local_XC_KIND,CUTOFF=l_CUTOFF,GAUGE=.TRUE.) if (ioX/=0) goto 1 ! ! PARS_1 @@ -127,10 +129,11 @@ integer function ioX(X,Xw,ID) & VAR=" Rl vectors in the sum :",I0=X%ngostnts,WARN=.true.,OP=(/"<="/)) call io_elemental(ID,& & VAR=" [r,Vnl] included :",L0=X%Vnl_included,WARN=.true.,OP=(/"=="/)) - if (ver_is_gt_or_eq(ID,(/3,0,9/))) then - call io_elemental(ID,& -& VAR=" Longitudinal Gauge :",L0=X%long_gauge,WARN=.true.,OP=(/"=="/)) + ! + if (.not.ver_is_gt_or_eq(ID,(/3,0,9/)).and.read_is_on(ID)) then + call io_elemental(ID,VAR=" Longitudinal Gauge :",L0=local_long_gauge) endif + ! call io_elemental(ID,& & VAR=" Field direction :",R1=X%q0,CHECK=.true.,OP=(/"==","==","=="/)) ! diff --git a/src/io/io_DIPOLES.F b/src/io/io_DIPOLES.F new file mode 100644 index 0000000000..6eeb8f5bab --- /dev/null +++ b/src/io/io_DIPOLES.F @@ -0,0 +1,131 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +integer function io_DIPOLES(X,Xen,ID) + ! + use pars, ONLY:SP,schlen + use units, ONLY:HARTREE + use X_m, ONLY:X_t,DIP_iR,X_alloc,Dipole_bands_ordered,& +& Dipole_uses_shifted_grids + use electrons, ONLY:levels,n_sp_pol + use R_lattice, ONLY:nXkibz,q0_def_norm + use matrix_operate,ONLY:mat_c2r,mat_r2c + use IO_m, ONLY:io_connect,io_disconnect,io_sec,io_control,& +& io_elemental,io_status,io_bulk,IO_OUTDATED_DB,& +& read_is_on,write_is_on,io_header,io_mode,io_com,io_action + use global_XC, ONLY:Dipole_WF_xc_string,loaded_WF_xc_string + ! + implicit none + type(X_t)::X + type(levels)::Xen + integer ::ID + ! + ! Work Space + ! + integer :: i1,ixyz,sec_size,i_spin,IDP + integer :: db_nbm,db_nbf,db_nb(2) + character(schlen) :: VAR_name + logical :: local_long_gauge + real(SP),allocatable :: DIP_iR_disk(:,:,:) + ! + io_DIPOLES=io_connect(desc='dipoles',type=2,ID=ID) + ! + if (io_DIPOLES/=0) goto 1 + ! + if (any((/io_sec(ID,:)==1/))) then + ! + io_DIPOLES=io_header(ID,R_LATT=.true.,WF=.true.,IMPOSE_SN=.true.,T_EL=.true.) + if (io_DIPOLES/=0) goto 1 + ! + sec_size=16 + ! + call io_elemental(ID,VAR="PARS",VAR_SZ=sec_size,MENU=0) + call io_elemental(ID,DB_I1=db_nb,& +& VAR=" X band range :",I1=X%ib,CHECK=.true.,OP=(/">=","<="/)) + call io_elemental(ID,UNIT=HARTREE,& +& VAR=" X e/h energy range [ev]:",R1=X%ehe,CHECK=.true.,OP=(/">=","<="/)) + call io_elemental(ID,DB_I0=db_nbm,& +& VAR=" Metallic bands :",I0=Xen%nbm,CHECK=.true.,OP=(/"=="/)) + call io_elemental(ID,DB_I0=db_nbf,& +& VAR=" Filled bands :",I0=Xen%nbf,CHECK=.true.,OP=(/"=="/)) + call io_elemental(ID,& +& VAR=" RL vectors in the sum :",I0=X%ngostnts,WARN=.true.,OP=(/"<="/)) + call io_elemental(ID,& +& VAR=" [r,Vnl] included :",L0=X%Vnl_included,CHECK=.true.,OP=(/"=="/)) + call io_elemental(ID,& +& VAR=" Transitions ordered :",L0=Dipole_bands_ordered,CHECK=.true.,OP=(/"=="/)) + call io_elemental(ID,& +& VAR=" Using shifted grids :",L0=Dipole_uses_shifted_grids,CHECK=.true.,OP=(/"=="/)) + call io_elemental(ID,& +& VAR=" Field momentum norm :",R0=q0_def_norm,CHECK=.true.,OP=(/"=="/)) + ! + call io_elemental(ID,VAR="",VAR_SZ=0,MENU=0) + ! + ! Wavefunctions xc + ! + call io_elemental(ID,VAR='WAVE_FUNC_XC',CH0="",VAR_SZ=1,MENU=0) + call io_elemental(ID,DB_CH0=Dipole_WF_xc_string,CH0=loaded_WF_xc_string,& +& VAR=' Wavefunctions :',CHECK=.true.,OP=(/"=="/)) + call io_elemental(ID,VAR="",VAR_SZ=0) + ! + io_DIPOLES=io_status(ID) + if (io_DIPOLES/=0) goto 1 + endif + ! + if (.not.Dipole_bands_ordered) db_nbm=db_nb(2) + ! + ! On disk the size is DIP_iR_or_P(3,db_nb(2),db_nbm,nXkibz) + ! + sec_size=3*db_nb(2)*db_nbm + if (any((/io_sec(ID,:)==2/))) then + ! + allocate(DIP_iR_disk(db_nb(2),db_nbm,2)) + ! + if(read_is_on(ID)) call X_alloc('DIP_iR',(/3,X%ib(2),db_nbm,nXkibz/)) + ! + do i1=1,nXkibz + ! + do ixyz=1,3 + ! + do i_spin=1,n_sp_pol + ! + if (write_is_on(ID)) call mat_c2r(DIP_iR(ixyz,:,:,i1,i_spin),DIP_iR_disk) + ! + write (VAR_name,'(3(a,i4.4))') 'DIP_iR_k_',i1,'_xyz_',ixyz,'_spin_',i_spin + call io_bulk(ID,VAR=trim(VAR_name),VAR_SZ=shape(DIP_iR_disk)) + call io_bulk(ID,R3=DIP_iR_disk) + ! + if (read_is_on(ID)) call mat_r2c(DIP_iR_disk,DIP_iR(ixyz,:,:,i1,i_spin)) + ! + enddo + ! + enddo + ! + enddo + ! + deallocate(DIP_iR_disk) + ! + endif + ! + ! +1 call io_disconnect(ID=ID) + ! +end function diff --git a/src/io/io_HF_and_locXC.F b/src/io/io_HF_and_locXC.F index 811cb2d2ce..a0bd3ca028 100644 --- a/src/io/io_HF_and_locXC.F +++ b/src/io/io_HF_and_locXC.F @@ -51,6 +51,8 @@ integer function io_HF_and_locXC(ID) ! io_HF_and_locXC=io_header(ID,QPTS=.true.,R_LATT=.true.,WF=.true.,IMPOSE_SN=.true.,T_EL=.true.,XC_KIND="G_WF force") ! + if (io_HF_and_locXC/=0) goto 1 + ! call io_elemental(ID,VAR="PARS",VAR_SZ=8,MENU=0) ! call io_elemental(ID,I0=QP_nb,CHECK=.true.,OP=(/"<="/)) diff --git a/src/io/io_header.F b/src/io/io_header.F index 32680ea03a..2a4508a495 100644 --- a/src/io/io_header.F +++ b/src/io/io_header.F @@ -19,10 +19,12 @@ ! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, ! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. ! -integer function io_header(ID,QPTS,R_LATT,WF,IMPOSE_SN,T_EL,KPTS,D_LATT,XC_KIND,CUTOFF) +integer function io_header(ID,QPTS,R_LATT,WF,IMPOSE_SN,T_EL,KPTS,D_LATT,XC_KIND,& +& CUTOFF,GAUGE,IMPOSE_GAUGE) ! use pars, ONLY:SP,schlen,lchlen use units, ONLY:HARTREE,Kelvin + use fields, ONLY:global_gauge use drivers, ONLY:list_dbs use LOGO, ONLY:code_version,code_revision use com, ONLY:msg,warning @@ -43,7 +45,7 @@ integer function io_header(ID,QPTS,R_LATT,WF,IMPOSE_SN,T_EL,KPTS,D_LATT,XC_KIND, ! implicit none integer :: ID - logical, optional :: QPTS,R_LATT,WF,IMPOSE_SN,T_EL,KPTS,D_LATT,CUTOFF + logical, optional :: QPTS,R_LATT,WF,IMPOSE_SN,T_EL,KPTS,D_LATT,CUTOFF,GAUGE,IMPOSE_GAUGE character(*),optional :: XC_KIND ! ! Work Space @@ -54,7 +56,7 @@ integer function io_header(ID,QPTS,R_LATT,WF,IMPOSE_SN,T_EL,KPTS,D_LATT,XC_KIND, real(SP) :: local_zero(3),D_LATT_vec_disk(3),D_LATT_vec(3),& & save_Tel(2) logical :: WARN,xc_kind_force - character(schlen) :: xc_string_kinds(10) + character(schlen) :: xc_string_kinds(10),DB_gauge ! io_header=0 ! @@ -299,6 +301,29 @@ integer function io_header(ID,QPTS,R_LATT,WF,IMPOSE_SN,T_EL,KPTS,D_LATT,XC_KIND, if (io_header/=0) goto 1 endif ! + ! GAUGE + ! + if (present(GAUGE).and.ver_is_gt_or_eq(ID,revision=625) ) then + call io_elemental(ID,VAR='GAUGE',CH0="",VAR_SZ=1,MENU=MENU) + MENU=0 + if (present(IMPOSE_GAUGE)) then + if (IMPOSE_GAUGE) then + call io_elemental(ID,CH0=global_gauge,& +& VAR=' Global Gauge :',CHECK=.true.,OP=(/"=="/)) + else + call io_elemental(ID,CH0=global_gauge,& +& VAR=' Global Gauge :',WARN=.true.,OP=(/"=="/)) + endif + else + call io_elemental(ID,CH0=global_gauge,& +& VAR=' Global Gauge :',WARN=.true.,OP=(/"=="/)) + endif + call io_elemental(ID,VAR="",VAR_SZ=0,MENU=0) + ! + io_header=io_status(ID) + if (io_header/=0) goto 1 + endif + ! 1 io_mode(ID)=MODE ! contains diff --git a/src/io/variables_BS.F b/src/io/variables_BS.F index fdcb995f0d..fa26d196ab 100644 --- a/src/io/variables_BS.F +++ b/src/io/variables_BS.F @@ -31,7 +31,7 @@ integer function variables_BS(ID,local_description,DB,X) & BS_n_g_exch,BS_n_g_W,BS_eh_en,BS_eh_win,BS_W_is_diagonal,& & BS_K_is_ALDA,BS_cpl_K_exchange,& & BS_cpl_K_corr - use IO_m, ONLY: io_elemental,io_status,ver_is_gt_or_eq + use IO_m, ONLY: io_elemental,io_status,ver_is_gt_or_eq,read_is_on implicit none ! integer ::ID @@ -43,7 +43,7 @@ integer function variables_BS(ID,local_description,DB,X) ! integer ::var_size,i_local_desc logical ::res_corr_disk,res_exch_disk,cpl_disk,& -& cpl_corr_disk,cpl_exch_disk,td_hf_disk +& cpl_corr_disk,cpl_exch_disk,td_hf_disk,local_long_gauge ! ! PARS1 ! @@ -152,14 +152,15 @@ integer function variables_BS(ID,local_description,DB,X) & VAR=" |[r,Vnl] included :",L0=X%Vnl_included,WARN=.true.,OP=(/"=="/),& & DESCRIPTION=local_description(i_local_desc)) ! - if (ver_is_gt_or_eq(ID,(/3,0,9/))) then + if (ver_is_gt_or_eq(ID,(/3,0,9/)).and..not.ver_is_gt_or_eq(ID,revision=626)) then ! ! v 3.0.8: Inclusion of %long_gauge ! i_local_desc=i_local_desc+1 call io_elemental(ID,& -& VAR=" |Longitudinal Gauge :",L0=X%long_gauge,WARN=.true.,OP=(/"=="/),& +& VAR=" |Longitudinal Gauge :",L0=local_long_gauge,& & DESCRIPTION=local_description(i_local_desc)) + ! endif ! i_local_desc=i_local_desc+1 diff --git a/src/modules/.objects b/src/modules/.objects index e78f7cd521..963d10c48a 100644 --- a/src/modules/.objects +++ b/src/modules/.objects @@ -1 +1,16 @@ -objs = mod_pars.o mod_units.o mod_stderr.o mod_par_proc.o mod_wrapper.o mod_collision.o mod_drivers.o mod_FFT.o mod_timing.o mod_logo.o mod_com.o mod_memory.o mod_R_lattice.o mod_par_indexes.o mod_wave_func.o mod_matrix_operate.o mod_D_lattice.o mod_frequency.o mod_wave_func.o mod_vec_operate.o mod_electrons.o mod_X.o mod_functions.o mod_zeros.o mod_pseudo.o mod_BS.o mod_QP.o mod_TDDFT.o mod_ACFDT.o mod_IO.o mod_fragments.o mod_xc_functionals.o mod_global_XC.o mod_X_output.o std_presets.o +#if defined _LIBXC +LIBXC_objects=mod_libxc_int.o +#endif +objs = mod_pars.o mod_units.o mod_stderr.o mod_par_proc.o \ + mod_wrapper.o mod_fields.o mod_collision.o mod_drivers.o mod_FFT.o \ + mod_timing.o mod_logo.o \ + mod_com.o mod_memory.o mod_R_lattice.o \ + mod_par_indexes.o mod_wave_func.o \ + mod_matrix_operate.o mod_D_lattice.o mod_frequency.o\ + mod_wave_func.o mod_vec_operate.o mod_electrons.o mod_X.o \ + mod_functions.o mod_zeros.o mod_pseudo.o \ + mod_BS.o mod_QP.o mod_TDDFT.o mod_ACFDT.o \ + mod_IO.o mod_fragments.o mod_xc_functionals.o mod_global_XC.o \ + mod_X_output.o std_presets.o \ + $(LIBXC_objects) + diff --git a/src/modules/mod_BS.F b/src/modules/mod_BS.F index 5d1ccbc68a..584f9bfddb 100644 --- a/src/modules/mod_BS.F +++ b/src/modules/mod_BS.F @@ -91,7 +91,6 @@ module BS logical :: BSS_uses_RIM logical :: BSS_Vnl_included logical :: BSS_uses_GreenF - logical :: BSS_long_gauge logical :: BSS_uses_partial_diago complex(SP),allocatable :: BSS_rhoq0(:) ! diff --git a/src/modules/mod_IO.F b/src/modules/mod_IO.F index fa9728b450..5d74c1dafe 100644 --- a/src/modules/mod_IO.F +++ b/src/modules/mod_IO.F @@ -58,7 +58,7 @@ module IO_m ! I/O errors ! integer, parameter:: IO_NO_DATABASE=-1,IO_INCOMPATIBLE_VAR=-2,IO_GENERIC_ERROR=-3,& -& IO_NO_BINDING_ERROR=-4,IO_NO_ERROR=0 +& IO_NO_BINDING_ERROR=-4,IO_NO_ERROR=0,IO_OUTDATED_DB=-5 ! ! Units ! @@ -127,9 +127,9 @@ subroutine io_bulk(ID,VAR,VAR_SZ,I0,I1,I2,I3,R0,R1,R2,R3,C1,L1,L3,IPOS) ! end subroutine ! - integer function io_header(ID,QPTS,R_LATT,WF,IMPOSE_SN,T_EL,KPTS,D_LATT,XC_KIND,CUTOFF) + integer function io_header(ID,QPTS,R_LATT,WF,IMPOSE_SN,T_EL,KPTS,D_LATT,XC_KIND,CUTOFF,GAUGE,IMPOSE_GAUGE) integer :: ID - logical,optional :: QPTS,R_LATT,WF,IMPOSE_SN,T_EL,KPTS,D_LATT,CUTOFF + logical,optional :: QPTS,R_LATT,WF,IMPOSE_SN,T_EL,KPTS,D_LATT,CUTOFF,GAUGE,IMPOSE_GAUGE character(*),optional :: XC_KIND end function ! diff --git a/src/modules/mod_R_lattice.F b/src/modules/mod_R_lattice.F index 8d6a07bdea..a38b0c0afe 100644 --- a/src/modules/mod_R_lattice.F +++ b/src/modules/mod_R_lattice.F @@ -104,6 +104,7 @@ module R_lattice real(SP) :: box_length(3) real(SP) :: cyl_cut character(schlen) :: cut_geometry + character(schlen) :: cut_description complex(SP),allocatable :: bare_qpg(:,:) logical :: CUTOFF_plus_RIM real(SP) :: cyl_vr_save @@ -192,6 +193,7 @@ subroutine cutoff_presets() box_length=0. cyl_length=0. cut_geometry='none' + cut_description='none' CUTOFF_plus_RIM=.false. end subroutine ! diff --git a/src/modules/mod_X.F b/src/modules/mod_X.F index 03b02033ce..9e6cf21e60 100644 --- a/src/modules/mod_X.F +++ b/src/modules/mod_X.F @@ -40,12 +40,13 @@ module X_m ! logical :: Vnl_commutator_warning ! - ! Dipole M.E.: gauges and kinetic energy + ! Dipole M.E. ! logical :: Dipole_bands_ordered - logical :: Eval_P_and_P2_only - complex(SP), allocatable :: DIP_iq_dot_r(:,:,:,:) - complex(SP), allocatable :: DIP_iR_or_P(:,:,:,:,:) + logical :: Dipole_uses_shifted_grids + complex(SP), allocatable :: DIP_q_dot_iR(:,:,:,:) + complex(SP), allocatable :: DIP_iR(:,:,:,:,:) + complex(SP), allocatable :: DIP_P(:,:,:,:,:) complex(SP), allocatable :: P_square(:,:,:,:) ! ! Absorption & Polarizability @@ -64,7 +65,7 @@ module X_m logical :: half_X_mat_only logical :: self_detect_E_range logical :: use_X_RIM - character(lchlen) :: long_path + ! type X_t integer :: whoami ! 1:Xo 2:em1s 3:em1d 4:pp integer :: ng @@ -78,7 +79,6 @@ module X_m real(SP) :: ppaE complex(SP):: Wd logical :: Vnl_included - logical :: long_gauge character(1) :: ordering character(schlen):: f_xc end type @@ -113,13 +113,17 @@ subroutine X_alloc(what,d) if (allocated(X_mat)) return allocate(X_mat(d(1),d(2),d(3)),stat=err) call mem_est(what,(/product(d)/),errors=(/err/)) - case('OptGrad') - if (allocated(DIP_iR_or_P)) return - allocate(DIP_iR_or_P(d(1),d(2),d(3),d(4),n_sp_pol),stat=err) + case('DIP_iR') + if (allocated(DIP_iR)) return + allocate(DIP_iR(d(1),d(2),d(3),d(4),n_sp_pol),stat=err) + call mem_est(what,(/product(d)*n_sp_pol/),errors=(/err/)) + case('DIP_P') + if (allocated(DIP_P)) return + allocate(DIP_P(d(1),d(2),d(3),d(4),n_sp_pol),stat=err) call mem_est(what,(/product(d)*n_sp_pol/),errors=(/err/)) - case('OptOsc') - if (allocated(DIP_iq_dot_r)) return - allocate(DIP_iq_dot_r(d(1),d(2),d(3),n_sp_pol),stat=err) + case('DIP_q_dot_iR') + if (allocated(DIP_q_dot_iR)) return + allocate(DIP_q_dot_iR(d(1),d(2),d(3),n_sp_pol),stat=err) call mem_est(what,(/product(d)*n_sp_pol/),errors=(/err/)) case('P_square') if (allocated(P_square)) return @@ -139,12 +143,14 @@ subroutine X_alloc(what,d) if (.not.allocated(X_mat)) return deallocate(X_mat) ! - case('OptGrad') - if (.not.allocated(DIP_iR_or_P)) return - deallocate(DIP_iR_or_P) - case('OptOsc') - if (.not.allocated(DIP_iq_dot_r)) return - deallocate(DIP_iq_dot_r) + case('DIP_iR') + if (.not.allocated(DIP_iR)) return + case('DIP_P') + if (.not.allocated(DIP_P)) return + deallocate(DIP_iR) + case('DIP_q_dot_iR') + if (.not.allocated(DIP_q_dot_iR)) return + deallocate(DIP_q_dot_iR) case('P_square') if (.not.allocated(P_square)) return deallocate(P_square) @@ -162,7 +168,6 @@ subroutine X_duplicate(Xi,Xo) Xo%ngostnts=Xi%ngostnts Xo%ordering=Xi%ordering Xo%Vnl_included=Xi%Vnl_included - Xo%long_gauge=Xi%long_gauge Xo%ehe=Xi%ehe Xo%q0=Xi%q0 Xo%cg_percentual=Xi%cg_percentual @@ -179,7 +184,6 @@ subroutine X_reset(X,type) X%iq=(/1,-1/) X%whoami=0 X%Vnl_included=.false. - X%long_gauge=.false. if (present(type)) X%whoami=type X%ehe=(/-1.,-1./)/HARTREE X%q0=(/1.,0.,0./) diff --git a/src/modules/mod_X_output.F b/src/modules/mod_X_output.F index 2beaf6fa26..6ded822fe6 100644 --- a/src/modules/mod_X_output.F +++ b/src/modules/mod_X_output.F @@ -25,6 +25,7 @@ module X_output use stderr, ONLY:intc use com, ONLY:msg use X_m, ONLY:alpha_dim + use fields, ONLY:global_gauge implicit none ! character(15) :: headers(7) @@ -32,11 +33,11 @@ module X_output ! contains ! - subroutine X_write_messages_before_headers(iq,GreenF,Vnl_included,long_gauge,ordering) + subroutine X_write_messages_before_headers(iq,GreenF,Vnl_included,ordering) use drivers, ONLY:l_chi,l_bss use global_XC, ONLY:X_E_xc_string,X_WF_xc_string,K_E_xc_string,K_WF_xc_string integer :: iq - logical :: Vnl_included,long_gauge,GreenF + logical :: Vnl_included,GreenF character(1) :: ordering ! WS character(schlen) :: message @@ -63,9 +64,11 @@ subroutine X_write_messages_before_headers(iq,GreenF,Vnl_included,long_gauge,ord ! Gauges ! if (iq==1) then - message='- The gauge is transverse -' - if (Vnl_included) message='- The gauge is transverse + [r,Vnl] -' - if (long_gauge) message='- The gauge is longitudinal -' + message='- Using the Length Gauge -' + if (trim(global_gauge)=='velocity') message='- Using the Velocity Gauge -' + call msg('o eps eel alpha','# ',trim(message),INDENT=0) + message='- [r,Vnl] is *NOT* included -' + if (Vnl_included) message='- [r,Vnl] *is* included -' call msg('o eps eel alpha','# ',trim(message),INDENT=0) endif ! diff --git a/src/modules/mod_com.F b/src/modules/mod_com.F index 7719175a90..875006872a 100644 --- a/src/modules/mod_com.F +++ b/src/modules/mod_com.F @@ -183,6 +183,8 @@ subroutine of_open_close(of_name,mode) integer :: i2,i3,file_index character(lchlen):: local_file_name ! + of_name=trim(of_name) + ! if (len_trim(of_name)==0) return ! if (present(mode)) then diff --git a/src/modules/mod_drivers.F b/src/modules/mod_drivers.F index 090fb687fd..677b692fc4 100644 --- a/src/modules/mod_drivers.F +++ b/src/modules/mod_drivers.F @@ -79,6 +79,12 @@ module drivers logical :: l_sc_cohsex logical :: l_sc_contains_lda ! + ! OEP approximations + ! + logical :: l_oep_exact + logical :: l_oep_kli + logical :: l_oep_slater + ! ! The XC types (Fxc) ! logical :: l_alda_fxc diff --git a/src/modules/mod_fields.F b/src/modules/mod_fields.F new file mode 100644 index 0000000000..4e826c742f --- /dev/null +++ b/src/modules/mod_fields.F @@ -0,0 +1,35 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +module fields + ! + use pars, ONLY:SP,schlen,lchlen + ! + ! Main Gauge + ! + character(schlen) :: global_gauge + ! + ! + ! Path for the shifted grids + ! + character(lchlen) :: grid_path + ! +end module diff --git a/src/modules/mod_functions.F b/src/modules/mod_functions.F index aa64b2afab..6ebb14edc5 100644 --- a/src/modules/mod_functions.F +++ b/src/modules/mod_functions.F @@ -86,13 +86,64 @@ pure function bose_decay(E) ! end function ! + ! Lifetime functions e2et/h2ht + !------------------------------ + ! + ! Gamma_n = 2 i \sum_m { -/+ i Im[e^-1(e_n -e_m) (2-f+bose_f) <- e2et + ! + i Im[e^-1(e_m -e_n) (f+bose_f) <- h2ht } + ! + ! where - for T-ordered theory, + for causal (finite Tel) + ! + function e2et(is,os,E,F) + ! + use electrons, ONLY:levels + use drivers, ONLY:Finite_Tel + integer ::is(2),os(2),e2et + type(levels) ::E + real(SP) :: F + real(SP) :: dE !ws + e2et=0 + ! + ! "Electron 2 Electron" decay + ! + dE=E%E(is(1),is(2),1)-E%E(os(1),os(2),1) + ! + F= -(2.-E%f(os(1),os(2),1)+bose_f(dE)) + if (Finite_Tel) F=(2.-E%f(os(1),os(2),1)+bose_f(dE)) + ! + if (dE>0..and.abs(F)>epsilon(1.)) e2et=1 + if (e2et==0) F=0. + end function + ! + function h2ht(is,os,E,F) + ! + use electrons, ONLY:levels + use drivers, ONLY:Finite_Tel + integer ::is(2),os(2),h2ht + type(levels) ::E + real(SP) :: F + ! + ! Work Space + ! + real(SP) :: dE + h2ht=0 + ! + !"Hole 2 Hole" decay + ! + dE=E%E(os(1),os(2),1)-E%E(is(1),is(2),1) + ! + F=E%f(os(1),os(2),1)+bose_f(dE) + ! + if (dE>0..and.abs(F)>epsilon(1.)) h2ht=1 + if (h2ht==0) F=0. + end function ! logical function K_scatter(E,Ep,Er,Dr) ! real(SP)::E,Ep,Er(2),Dr(2),D K_scatter=.false. D=Dr(1)+(E-Er(1))/(Er(2)-Er(1))*(Dr(2)-Dr(1)) - D=max(D,0.) + D=max(D,0._SP) K_scatter=abs(E-Ep)<=D end function ! diff --git a/src/modules/mod_global_XC.F b/src/modules/mod_global_XC.F index 7ff9c306fb..284294f2e2 100644 --- a/src/modules/mod_global_XC.F +++ b/src/modules/mod_global_XC.F @@ -79,7 +79,13 @@ module global_XC SC_HF = 302, & ! Hartree-Fock SC_COHSEX = 303, & ! Coulomb-Hole Screened-eXchange SC_EXXC = 304, & ! OEP-EXX + LDA correlation - SC_SRPA = 305 ! OEP-COHSEX (static RPA) + SC_SRPA = 305, & ! OEP-COHSEX (static RPA) + SC_EXX_KLI = 311, & ! OEP-EXX (KLI apprx) + SC_EXXC_KLI = 314, & ! OEP-EXX (KLI apprx) + LDA correlation + SC_SRPA_KLI = 315, & ! OEP-COHSEX (KLI apprx) + SC_EXX_SLT = 321, & ! OEP-EXX (SLT apprx) + SC_EXXC_SLT = 324, & ! OEP-EXX (SLT apprx) + LDA correlation + SC_SRPA_SLT = 325 ! OEP-COHSEX (SLT apprx) ! ! Self Energies integer, public, parameter :: & diff --git a/src/modules/mod_libxc_int.F b/src/modules/mod_libxc_int.F new file mode 100644 index 0000000000..a68c1c0697 --- /dev/null +++ b/src/modules/mod_libxc_int.F @@ -0,0 +1,122 @@ +! +! Copyright (C) 2000-2010 M. Gruning and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +module libxc_int + ! + ! Interface between libxc and yambo + ! + use xc_functionals, ONLY: N_Available_LDA + use xc_f90_lib_m + use xc_f90_types_m + use libxc_funcs_m + implicit none + ! + type xc_fnctl_t + integer :: family + integer :: kind + integer :: id + integer :: spin_channels + integer :: flags + type(xc_f90_pointer_t) :: conf + type(xc_f90_pointer_t) :: info + end type xc_fnctl_t + integer, parameter :: & + XC_FACTOR = 1000, & + XC_NOT_AVAILABLE = 999999, & + NOXC = 0, & + XC_LDA_C_KP = 99 !This is not defined within libxc + ! + ! This matrix of integers maps XC numbering of Yambo for LDA into libxc + ! Note that though with the same name the parameters (defined in libxc_funcs_m) + ! correspond in principle to a different number than in mod_xc_functionals + ! + integer, parameter, dimension(N_Available_LDA) :: LDA_Y2LIBXC =& +& (/XC_LDA_X, & !1 + XC_LDA_C_WIGNER, & !2 + XC_LDA_C_RPA, & !3 + XC_LDA_C_HL, & !4 + XC_LDA_C_GL, & !5 + XC_LDA_C_XALPHA, & !6 + XC_LDA_C_VWN, & !7 + XC_LDA_C_PZ, & !8 + XC_LDA_C_OB_PZ, & !9 + XC_LDA_C_PW, & !10 + XC_LDA_C_OB_PW, & !11 + XC_GGA_C_LYP, & !12 + XC_LDA_C_2D_AMGB, & !13 + XC_LDA_C_KP, & !14 + XC_LDA_XC_TETER93, & !15 + XC_LDA_C_VWN_RPA, & !16 + XC_LDA_C_PZ_MOD, & !17 + XC_LDA_C_PW_MOD, & !18 + XC_LDA_C_vBH/) + ! + contains + ! + integer function XC_y2libxc(Yambo_Func,Yambo_kind) + integer Yambo_Func,Yambo_kind, Tmp_Func, libxc_kind + ! + ! This function map the information from Yambo_Func & Yambo_kind + ! in one single number = [X] * 1000 + [C] following the libxc conventions + ! Status(22/07/10): Only LDA supported + ! + XC_y2libxc = 0 + if (Yambo_kind == NOXC) return + Tmp_Func = LDA_Y2LIBXC(Yambo_Func) + libxc_kind = Yambo_kind -1 ! Note that the kind in Yambo are shifted by +1 wrt libxc + if (Tmp_Func==XC_LDA_XC_TETER93) libxc_kind = XC_CORRELATION ! since this functional contains already X + if ((libxc_kind) == XC_EXCHANGE.or.(libxc_kind) == XC_EXCHANGE_CORRELATION)& + &XC_y2libxc = XC_LDA_X*1000 + if (libxc_kind == XC_CORRELATION.or.libxc_kind == XC_EXCHANGE_CORRELATION)& + &XC_y2libxc =XC_y2libxc + Tmp_Func + end function XC_y2libxc + ! + subroutine xc_setup_fnctl(f,fun,nsp) + ! + use xc_f90_types_m + use xc_f90_lib_m, ONLY:xc_f90_func_init,xc_f90_info_family,& + xc_f90_info_kind,xc_f90_info_flags + implicit none + integer, intent(in) :: fun, nsp + type(xc_fnctl_t), intent(out) :: f(2) + integer :: ixc + type(xc_f90_pointer_t) :: p,info + ! + call xc_f90_func_init(p,info,NOXC,1) + ! + f(1)%id = fun/XC_FACTOR + f(2)%id = fun-f(1)%id*XC_FACTOR + ! + do ixc=1,2 + if (f(ixc)%id==0) cycle + call xc_f90_func_init(p,info,f(ixc)%id,nsp) + f(ixc)%family=xc_f90_info_family(info) + f(ixc)%kind=xc_f90_info_kind(info) + f(ixc)%id=fun + f(ixc)%spin_channels=nsp + f(ixc)%flags=xc_f90_info_flags(info) + f(ixc)%conf=p + f(ixc)%info=info + enddo + ! + end subroutine xc_setup_fnctl + ! +end module libxc_int diff --git a/src/modules/mod_xc_functionals.F b/src/modules/mod_xc_functionals.F index 8db40b73a4..41aa926f4c 100644 --- a/src/modules/mod_xc_functionals.F +++ b/src/modules/mod_xc_functionals.F @@ -1,7 +1,7 @@ ! ! Copyright (C) 2003 M. Marques, A. Castro, A. Rubio, G. Bertsch ! -! Copyright (C) 2000-2010 A. Marini, M. Gruning and the YAMBO team +! Copyright (C) 2000-2010 A. Marini and the YAMBO team ! http://www.yambo-code.org ! ! This file is distributed under the terms of the GNU @@ -33,9 +33,6 @@ module xc_functionals real(SP), public, allocatable:: E_xc(:) real(SP), public, allocatable:: F_xc(:) ! - ! whether interface to dft code used libxc (for compatibility issues) - ! - integer, public, parameter :: XC_USED_LIBXC = 4 integer, public, parameter :: & XC_UNPOLARIZED = 1, & ! Spin unpolarized XC_POLARIZED = 2 ! Spin polarized @@ -46,11 +43,23 @@ module xc_functionals ! ! Kinds integer, public, parameter :: & + XC_NOT_AVAILABLE = -1, & NONE = 0, & XC_EXCHANGE = 1, & XC_CORRELATION = 2, & XC_EXCHANGE_CORRELATION = 3 ! + ! Families + integer, public, parameter :: & + XC_FAMILY_UNKNOWN = -1, & + XC_FAMILY_NONE = 0, & + XC_FAMILY_LDA = 1, & + XC_FAMILY_GGA = 2, & + XC_FAMILY_MGGA = 4, & + XC_FAMILY_LCA = 8, & + XC_FAMILY_OEP = 16, & + XC_FAMILY_HYB_GGA = 32 + ! ! the LDAs integer, public, parameter :: & XC_LDA_X = 1, & ! Exchange @@ -66,7 +75,12 @@ module xc_functionals XC_LDA_C_OB_PW = 11, & ! Ortiz & Ballone (PW) XC_LDA_C_LYP = 12, & ! Lee, Yang, & Parr LDA XC_LDA_C_AMGB = 13, & ! Attacalite et al - XC_LDA_C_KP = 14 ! Kurth & Perdew + XC_LDA_C_KP = 14, & ! Kurth & Perdew + XC_LDA_XC_TETER93 = 15, & ! Teter Parametrization(93) + XC_LDA_C_VWN_RPA = 16, & ! Vosko, Wilk, & Nussair (RPA) + XC_LDA_C_PZ_MOD = 17, & ! Perdew & Zunger (Modified) + XC_LDA_C_PW_MOD = 18, & ! Perdew & Wang (Modified) + XC_LDA_C_vBH = 19 ! von Barth & Hedin ! ! the GGAs integer, public, parameter :: & @@ -79,8 +93,48 @@ module xc_functionals XC_MGGA_X_TPSS = 201, & ! Perdew, Tao, Staroverov & Scuseria exchange XC_MGGA_C_TPSS = 202 ! Perdew, Tao, Staroverov & Scuseria correlation ! + ! strings with XC description + ! + integer, public, parameter :: N_Available_LDA = 19 + character(*), parameter, dimension( N_Available_LDA) :: XC_LDA_DESCRIPTION =& +& (/"Exchange only ",& + "Wigner parametrization ",& + "Random Phase Approximation ",& + "Hedin & Lundqvist ",& + "Gunnarson & Lundqvist ",& + "Slaters Xalpha ",& + "Vosko, Wilk, & Nussair ",& + "Perdew & Zunger ",& + "Ortiz & Ballone(PZ) ",& + "Perdew & Wang ",& + "Ortiz & Ballone(PW) ",& + "Lee, Yang, & Parr LDA ",& + "Attacalite et al ",& + "Kurth & Perdew ",& + "Teter Parametrization(93) ",& + "Vosko, Wilk, & Nussair(RPA) ",& + "Perdew & Zunger (Modified) ",& + "Perdew & Wang (Modified) ",& + "von Barth & Hedin "/) + ! contains ! + logical function xc_OnlyInLibxc(func) + integer func + xc_OnlyInLibxc = .false. + if (func==XC_LDA_XC_TETER93.or. & + func==XC_LDA_C_RPA.or. & + func==XC_LDA_C_GL.or.& + func==XC_LDA_C_VWN.or.& + func==XC_LDA_C_VWN_RPA.or.& + func==XC_LDA_C_OB_PZ.or.& + func==XC_LDA_C_OB_PW.or.& + func==XC_LDA_C_PZ_MOD.or.& + func== XC_LDA_C_vBH.or.& + func==XC_LDA_C_AMGB.or.& + func==XC_LDA_C_PW_MOD) xc_OnlyInLibxc=.true. + end function xc_OnlyInLibxc + ! character(lchlen) function xc_string(kind,functional) ! integer :: kind,functional @@ -88,7 +142,7 @@ character(lchlen) function xc_string(kind,functional) integer :: i1,fun ! select case (kind) - case(-1,NONE) + case(NONE) ch_kind=' ' case(XC_EXCHANGE) ch_kind='(x)' @@ -98,28 +152,15 @@ character(lchlen) function xc_string(kind,functional) ch_kind='(xc)' end select ! - select case (functional) - case(-1,3,5,7,9) + select case(functional) + case (XC_NOT_AVAILABLE) xc_string='Unknown '//trim(ch_kind) - case(NONE) - xc_string='none' - case(XC_LDA_X) - xc_string='Exchange only '//trim(ch_kind) - case(XC_LDA_C_WIGNER) - xc_string='Wigner '//trim(ch_kind) - case(XC_LDA_C_HL) - xc_string='Hedin & Lundqvist '//trim(ch_kind) - case(XC_LDA_C_XALPHA) - xc_string='Slaters Xalpha '//trim(ch_kind) - case(XC_LDA_C_PZ) - xc_string='Perdew & Zunger '//trim(ch_kind) - case(XC_LDA_C_PW) - xc_string='Perdew & Wang '//trim(ch_kind) - case(XC_LDA_C_KP) - xc_string='Kurth & Perdew '//trim(ch_kind) - end select - ! + case (NONE) + xc_string='None '//trim(ch_kind) + case (XC_LDA_X:N_Available_LDA) + xc_string=trim(XC_LDA_DESCRIPTION(functional))//' '//trim(ch_kind) + end select + ! end function -! - ! + ! end module xc_functionals diff --git a/src/modules/std_presets.F b/src/modules/std_presets.F index f57660658a..01aa429efe 100644 --- a/src/modules/std_presets.F +++ b/src/modules/std_presets.F @@ -25,9 +25,10 @@ subroutine std_presets(INSTR,IND,OD,JS,COM_DIR) use pars, ONLY:SP use units, ONLY:HARTREE,FS2AUT,AU2Vmm use LOGO, ONLY:ID_logo,ID_logo_stderr + use fields, ONLY:global_gauge,grid_path use X_m, ONLY:current_iq,self_detect_E_range,half_X_mat_only,eps_2_alpha,& -& alpha_dim,use_X_RIM,long_path,X_RIM_nkpts,Dipole_bands_ordered,& -& eval_alpha,Vnl_commutator_warning,Eval_P_and_P2_only,& +& alpha_dim,use_X_RIM,X_RIM_nkpts,Dipole_bands_ordered,& +& eval_alpha,Vnl_commutator_warning,Dipole_uses_shifted_grids,& & q_plus_G_direction,Q_Shift_Order use QP_m, ONLY:QP_dSc_steps,QP_n_W_freqs,QP_G_Zoom_treshold,& & QP_dSc_test,QP_solver,QP_G_damp,QP_dSc_delta,& @@ -62,7 +63,7 @@ subroutine std_presets(INSTR,IND,OD,JS,COM_DIR) & BS_res_mode,BS_K_dim,BS_cpl_mode,BS_not_const_eh_f,& & BSS_mode,BSS_n_freqs,BSS_n_descs,BSS_er,BSS_dr,& & BSS_q0,Haydock_treshold,BS_K_is_ALDA,BSS_uses_RIM,BSS_damp_reference,& -& BS_columns,BSS_Vnl_included,BSS_long_gauge,BSS_uses_GreenF +& BS_columns,BSS_Vnl_included,BSS_uses_GreenF use TDDFT, ONLY:FXC_description,FXC_type,FXC_n_descs,FXC_n_g_corr,& & FXC_per_memstps,FXC_LRC_alpha,FXC_SVD_digits,& & FXC_is_causal @@ -71,6 +72,9 @@ subroutine std_presets(INSTR,IND,OD,JS,COM_DIR) use memory_m, ONLY:mem_reset use zeros, ONLY:zero_norm,k_iku_zero,k_rlu_zero,G_iku_zero,zero_dfl use xc_functionals, ONLY:GS_xc_FUNCTIONAL,GS_xc_KIND +#if defined _TWO_LEVELS + use real_time, ONLY:l_3levels,en_levels,dip,lambda_diag,lambda_offdiag +#endif ! implicit none ! @@ -274,19 +278,17 @@ subroutine std_presets(INSTR,IND,OD,JS,COM_DIR) self_detect_E_range=.false. half_X_mat_only=.false. Dipole_bands_ordered=.true. + Dipole_uses_shifted_grids=.false. use_X_RIM=.false. eps_2_alpha=1._SP alpha_dim='adim' - long_path='none' + global_gauge='length' + grid_path='none' eval_alpha=.FALSE. Vnl_commutator_warning=.FALSE. q_plus_G_direction=0. Q_Shift_Order=1 ! - !Dipole - ! - Eval_P_and_P2_only=.FALSE. - ! !QPm ! QP_n_states=0 @@ -345,7 +347,6 @@ subroutine std_presets(INSTR,IND,OD,JS,COM_DIR) Haydock_treshold=-0.02 BS_columns=0 BSS_Vnl_included=.FALSE. - BSS_long_gauge=.FALSE. BSS_uses_GreenF=.FALSE. ! ! TDDFT @@ -372,6 +373,16 @@ subroutine std_presets(INSTR,IND,OD,JS,COM_DIR) ! ! ! +#if defined _TWO_LEVELS + ! + l_3levels =.FALSE. + en_levels =0._SP + lambda_offdiag=0._SP + lambda_diag =0._SP + dip =0._SP + ! +#endif + ! ! ! end subroutine diff --git a/src/parser/mod_itm.F b/src/parser/mod_itm.F index 7174827ceb..c49ad99950 100644 --- a/src/parser/mod_itm.F +++ b/src/parser/mod_itm.F @@ -48,6 +48,7 @@ module it_m integer, parameter :: V_general=6 integer, parameter :: V_resp=7 integer, parameter :: V_real_time=8 + integer, parameter :: V_all=99 ! integer, parameter :: maxrnlvls=50 integer, parameter :: maxflags =50 @@ -948,6 +949,7 @@ logical function Iverbose_enough(var_name,nvar,var_val) if (it_verbose_Ilevel(i1)==0) cycle if (trim(it_verbose_Ivars(i1))==trim(var_name))then if(all(var_val==it_verbose_Idefs(i1,:nvar)).and.& +& infile_verbosity/=V_all.and.& & infile_verbosity/=it_verbose_Ilevel(i1)) return endif enddo @@ -966,6 +968,7 @@ logical function Rverbose_enough(var_name,nvar,var_val) if (it_verbose_Rlevel(i1)==0) cycle if (trim(it_verbose_Rvars(i1))==trim(var_name))then if(all(var_val==it_verbose_Rdefs(i1,:nvar)).and.& +& infile_verbosity/=V_all.and.& & infile_verbosity/=it_verbose_Rlevel(i1)) return endif enddo @@ -984,6 +987,7 @@ logical function Cverbose_enough(var_name,nvar,var_val) if (it_verbose_Clevel(i1)==0) cycle if (trim(it_verbose_Cvars(i1))==trim(var_name))then if(var_val(1)==it_verbose_Cdefs(i1).and.& +& infile_verbosity/=V_all.and.& & infile_verbosity/=it_verbose_Clevel(i1)) return endif enddo @@ -999,6 +1003,7 @@ logical function Chverbose_enough(var_name,var_val) if (it_verbose_Chlevel(i1)==0) cycle if (trim(it_verbose_Chvars(i1))==trim(var_name))then if(var_val==it_verbose_Chdefs(i1).and.& +& infile_verbosity/=V_all.and.& & infile_verbosity/=it_verbose_Chlevel(i1)) return endif enddo @@ -1016,6 +1021,7 @@ logical function Fverbose_enough(var_name,ivar) if (it_verbose_Flevel(i1)==0) cycle if (trim(it_verbose_Fvars(i1))==trim(var_name))then if(fstatus(ivar)>0.and.& +& infile_verbosity/=V_all.and.& & infile_verbosity/=it_verbose_Flevel(i1)) return endif enddo diff --git a/src/pol_function/Dipole_driver.F b/src/pol_function/Dipole_driver.F index e04bfd8af6..11e318d035 100644 --- a/src/pol_function/Dipole_driver.F +++ b/src/pol_function/Dipole_driver.F @@ -21,16 +21,17 @@ ! subroutine Dipole_driver(Xen,Xk,X,field_dir) ! - use drivers, ONLY:l_optics,l_real_time - use pars, ONLY:SP - use com, ONLY:warning + use drivers, ONLY:l_optics,l_real_time,l_sc_run + use pars, ONLY:SP,cI + use com, ONLY:warning,error use drivers, ONLY:l_col_cut use electrons, ONLY:levels,n_sp_pol use par_proc_m, ONLY:pp_redux_wait use D_lattice, ONLY:nsym,i_time_rev,dl_sop,sop_inv use R_lattice, ONLY:bz_samp,q0_def_norm,q_norm,bare_qpg - use X_m, ONLY:DIP_iq_dot_r,X_alloc,X_t,DIP_iR_or_P,long_path,Dipole_bands_ordered,& -& Vnl_commutator_warning,Eval_P_and_P2_only,P_square + use fields, ONLY:grid_path,global_gauge + use X_m, ONLY:DIP_q_dot_iR,X_alloc,X_t,DIP_iR,Dipole_bands_ordered,& +& Vnl_commutator_warning,Dipole_uses_shifted_grids,DIP_P use IO_m, ONLY:io_control,OP_RD_CL,OP_WR_CL,VERIFY,REP use wave_func, ONLY:wf_ng,WF_load,WF_free use vec_operate, ONLY:v_norm @@ -44,15 +45,15 @@ subroutine Dipole_driver(Xen,Xk,X,field_dir) ! Work Space ! integer :: ik,i1,ic,iv,is,i_spin - real(SP) :: field_dir_rot(3) + real(SP) :: field_dir_rot(3),E_m_Ep logical :: use_trans_gauge ! !I/O ! integer :: ID,io_err - integer, external :: ioOSTNTS + integer, external :: io_DIPOLES ! - if (allocated(DIP_iq_dot_r)) return + if (allocated(DIP_q_dot_iR)) return ! ! Transitions ordering ! @@ -64,17 +65,15 @@ subroutine Dipole_driver(Xen,Xk,X,field_dir) ! in real-time simulations P and P^2 are ALWAYS calculated => no band ordering ! ! - ! Check first if ostnts DB is already done + Dipole_uses_shifted_grids=trim(grid_path)/='none' ! - X%ngostnts=wf_ng ! + ! Check first if ostnts DB is already done ! - ! Swith off I/O in case gradient is calculated - if (Eval_P_and_P2_only) call IO_and_Messaging_switch("-io_in -io_out") + X%ngostnts=wf_ng ! call io_control(ACTION=OP_RD_CL,COM=REP,SEC=(/1,2/),MODE=VERIFY,ID=ID) - io_err=ioOSTNTS(X,Xen,ID) - ! + io_err=io_DIPOLES(X,Xen,ID) ! if (io_err.ne.0) then ! @@ -84,12 +83,13 @@ subroutine Dipole_driver(Xen,Xk,X,field_dir) ! in the new basis. So all transitions are needed ! if (Dipole_bands_ordered) then - call X_alloc('OptGrad',(/3,X%ib(2),Xen%nbm,Xk%nibz/)) + call X_alloc('DIP_iR',(/3,X%ib(2),Xen%nbm,Xk%nibz/)) else - call X_alloc('OptGrad',(/3,X%ib(2),X%ib(2),Xk%nibz/)) + call X_alloc('DIP_iR',(/3,X%ib(2),X%ib(2),Xk%nibz/)) endif - DIP_iR_or_P=(0.,0.) + DIP_iR=(0.,0.) ! + ! ! call WF_load(0,1,X%ib,(/1,Xk%nibz/),space='G',title='-Oscillators/G space') ! @@ -102,18 +102,16 @@ subroutine Dipole_driver(Xen,Xk,X,field_dir) ! call WF_free() ! - ! Write the DIP_iR_or_P to its database (Only in the case of "R") - ! call io_control(ACTION=OP_WR_CL,COM=REP,SEC=(/1,2/),ID=ID) - io_err=ioOSTNTS(X,Xen,ID) + io_err=io_DIPOLES(X,Xen,ID) ! end if ! ! ! The field direction and the gamma-point norm must be renormalized here in case the - ! oscillator strengths have been calculated using the long. gauge. + ! oscillator strengths have been calculated using the shifted grids. ! In this case q0_def_norm is not the default one but corresponds - ! to the norm of the longitudinal shift. + ! to the norm of the grid shift. ! field_dir=field_dir*q0_def_norm/v_norm(field_dir) q_norm(1)=q0_def_norm @@ -126,10 +124,10 @@ subroutine Dipole_driver(Xen,Xk,X,field_dir) Vnl_commutator_warning=.TRUE. endif ! - ! Calculate the q-dependent oscillators (OptOsc) + ! Calculate the q-dependent oscillators ! - call X_alloc('OptOsc',(/X%ib(2),Xen%nbm,Xk%nbz/)) - DIP_iq_dot_r = (0.0_SP,0.0_SP) + call X_alloc('DIP_q_dot_iR',(/X%ib(2),Xen%nbm,Xk%nbz/)) + DIP_q_dot_iR = (0.0_SP,0.0_SP) ! do i1 = 1,Xk%nbz ! @@ -147,21 +145,19 @@ subroutine Dipole_driver(Xen,Xk,X,field_dir) ! ! DIP_iq_dot_r = i q . < v k | r | c k > ! - DIP_iq_dot_r(ic,iv,i1,i_spin) = & -& dot_product( field_dir_rot, DIP_iR_or_P(:,ic,iv,ik,i_spin) ) - if ( is > nsym/(i_time_rev+1) ) DIP_iq_dot_r(ic,iv,i1,i_spin) = & -& dot_product( DIP_iR_or_P(:,ic,iv,ik,i_spin), field_dir_rot ) + DIP_q_dot_iR(ic,iv,i1,i_spin) = dot_product( field_dir_rot, DIP_iR(:,ic,iv,ik,i_spin) ) + if ( is > nsym/(i_time_rev+1) ) DIP_q_dot_iR(ic,iv,i1,i_spin) = dot_product( DIP_iR(:,ic,iv,ik,i_spin), field_dir_rot ) ! - ! ! enddo - DIP_iq_dot_r(iv,iv,i1,:) = (1.0_SP,0.0_SP) + DIP_q_dot_iR(iv,iv,i1,:) = (1.0_SP,0.0_SP) enddo + ! enddo call pp_redux_wait() ! ! Clean up ! - call X_alloc('OptGrad') ! Used in O_transverse/O_longitudinal/ioOSTNTS + call X_alloc('DIP_iR') ! end subroutine diff --git a/src/pol_function/Dipole_transverse.F b/src/pol_function/Dipole_transverse.F index 33b793ab35..97b34eae84 100644 --- a/src/pol_function/Dipole_transverse.F +++ b/src/pol_function/Dipole_transverse.F @@ -20,6 +20,8 @@ ! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. ! subroutine Dipole_transverse(Xen,Xk,X) + ! + ! This routine returns and . ! use pars, ONLY:SP,pi use timing, ONLY:live_timing @@ -29,10 +31,11 @@ subroutine Dipole_transverse(Xen,Xk,X) use D_lattice, ONLY:alat,n_atomic_species,n_atoms_species use pseudo, ONLY:pp_kbv_dim,pp_n_l_comp,PP_free use R_lattice, ONLY:g_vec,bz_samp - use X_m, ONLY:X_t,DIP_iR_or_P,Dipole_bands_ordered,Eval_P_and_P2_only,P_square + use X_m, ONLY:X_t,DIP_iR,Dipole_bands_ordered,P_square use IO_m, ONLY:io_control,VERIFY,REP,OP_RD,RD_CL,RD use memory_m, ONLY:mem_est use wave_func, ONLY:wf,wf_ng,wf_state, WF_load + use wrapper, ONLY:M_by_V ! implicit none type(bz_samp), intent(in) :: Xk @@ -72,6 +75,8 @@ subroutine Dipole_transverse(Xen,Xk,X) enddo allocate(kbv(wf_ng,pp_kbv_dim,4),stat=alloc_err) call mem_est("KBV",(/size(kbv)/),errors=(/alloc_err/)) + else + X%Vnl_included=.false. endif ! ! Set up band limits @@ -125,10 +130,9 @@ subroutine Dipole_transverse(Xen,Xk,X) Ev_m_Ec=Xen%E(iv,ik,E_i_spin)-Xen%E(ic,ik,E_i_spin) if (associated(Xen%Eo)) Ev_m_Ec=Xen%Eo(iv,ik,E_i_spin)-Xen%Eo(ic,ik,E_i_spin) ! - if (any((/-Ev_m_Ec0.0_SP,& -& -Ev_m_Ec>X%ehe(2).and.X%ehe(2)>0.0_SP/))) cycle + if (any((/-Ev_m_Ec0.0_SP,-Ev_m_Ec>X%ehe(2).and.X%ehe(2)>0.0_SP/))) cycle ! - if (abs(Ev_m_Ec)<=1.E-5_SP.and..not.Eval_P_and_P2_only) cycle + if (abs(Ev_m_Ec)<=1.E-5_SP) cycle ! ivfft=wf_state(iv,ik,i_spin) ! @@ -147,23 +151,19 @@ subroutine Dipole_transverse(Xen,Xk,X) ! ! rho(4) = \sum_G u^*_{iv ik}(G) ( k + G )^2 u_{ic ik}(G) ! -#if defined _DOUBLE - call zgemv('C',wf_ng,wf_dim,(1._SP,0._SP),wf_xyz,wf_ng,wf(:,icfft),1,(0._SP,0._SP),rho,1) -#else - call cgemv('C',wf_ng,wf_dim,(1._SP,0._SP),wf_xyz,wf_ng,wf(:,icfft),1,(0._SP,0._SP),rho,1) -#endif + call M_by_V('C',wf_ng,wf_dim,(1._SP,0._SP),wf_xyz,wf_ng,wf(:,icfft),1,(0._SP,0._SP),rho,1) ! if (ioKB_err==0) call Dipole_kb_sum(ivfft,icfft,rho(:3),kbv) ! - if ( .not.Eval_P_and_P2_only ) rho(:)=rho(:)/Ev_m_Ec - ! ! as [x,p_x]=i we get [x,H] = [x,p^2/2]= i p_x ! ! DIP_iR_or_P(c,v) = i = i /(Ec-Ev) = i /(Ec-Ev) = ! = - / (Ec-Ev) = /(Ev-Ec) + ! + rho(1:3)=rho(1:3)/Ev_m_Ec ! ! - DIP_iR_or_P(:,ic,iv,ik,i_spin)=rho(:3) + DIP_iR(:,ic,iv,ik,i_spin)=rho(:3) ! ! enddo ! conduction band loop @@ -188,7 +188,7 @@ subroutine Dipole_transverse(Xen,Xk,X) endif ! do i_spin=1,n_sp_pol - call pp_redux_wait(DIP_iR_or_P(:,:,:,:,i_spin)) + call pp_redux_wait(DIP_iR(:,:,:,:,i_spin)) enddo ! end subroutine Dipole_transverse diff --git a/src/pol_function/O_driver.F b/src/pol_function/O_driver.F index 5347dc5c75..593d8e33e5 100644 --- a/src/pol_function/O_driver.F +++ b/src/pol_function/O_driver.F @@ -30,11 +30,11 @@ subroutine O_driver(Xen,Xk,q,wv,X) use memory_m, ONLY:mem_est use stderr, ONLY:intc,set_real_printed_length use drivers, ONLY:l_bs_fxc,l_alda_fxc - use frequency, ONLY:w_samp,rg_index_bg,bg_npts,cg_pt,cg_index_bg + use frequency, ONLY:w_samp use electrons, ONLY:levels,BZ_RIM_tot_nkpts use R_lattice, ONLY:bz_samp,q_norm,bare_qpg use com, ONLY:msg,of_open_close - use X_m, ONLY:X_t,X_epsilon,X_mat,X_alloc,X_poles_tab,X_fxc,& + use X_m, ONLY:X_t,X_epsilon,X_mat,X_alloc,X_fxc,& & use_X_RIM,eval_alpha,eps_2_alpha,O_eels use par_proc_m, ONLY:pp_redux_wait use wave_func, ONLY:WF_free @@ -180,7 +180,7 @@ subroutine O_driver(Xen,Xk,q,wv,X) ! ! Unfortunately some of the variables need in this second bunch of messages is setup only in X_os ! - call X_write_messages_before_headers(iq,associated(Xen%GreenF),X%Vnl_included,X%long_gauge,X%ordering) + call X_write_messages_before_headers(iq,associated(Xen%GreenF),X%Vnl_included,X%ordering) ! ! Titles ! @@ -271,17 +271,17 @@ subroutine O_driver(Xen,Xk,q,wv,X) ! call X_alloc('X') ! - deallocate(X_epsilon,wv%p) - call mem_est("W-p") + deallocate(X_epsilon) + ! + call FREQUENCIES_reset(wv) + ! if (allocated(X_fxc)) deallocate(X_fxc) - if (allocated(rg_index_bg)) deallocate(rg_index_bg) - deallocate(X_poles_tab,bg_npts,cg_pt,cg_index_bg) - call mem_est("X_poles_tab RGi BGn CGp CGi") + ! enddo ! ! CLEAN ! - call X_alloc('OptOsc') + call X_alloc('DIP_q_dot_iR') call WF_free() call pp_redux_wait call set_real_printed_length() diff --git a/src/pol_function/X_em1.F b/src/pol_function/X_em1.F index 2e3e48e6a3..cc73d2391c 100644 --- a/src/pol_function/X_em1.F +++ b/src/pol_function/X_em1.F @@ -25,9 +25,9 @@ integer function X_em1(Xen,Xk,q,X,Xw,APPEND_NO_VERIFY) ! use pars, ONLY:SP use drivers, ONLY:l_alda_fxc,l_gw0 - use X_m, ONLY:X_t,X_alloc,X_mat,X_poles_tab,self_detect_E_range + use X_m, ONLY:X_t,X_alloc,X_mat,self_detect_E_range use memory_m, ONLY:mem_est - use frequency, ONLY:w_samp,rg_index_bg,bg_npts,cg_pt,cg_index_bg + use frequency, ONLY:w_samp use R_lattice, ONLY:bz_samp use electrons, ONLY:levels use par_proc_m, ONLY:pp_redux_wait @@ -82,6 +82,19 @@ integer function X_em1(Xen,Xk,q,X,Xw,APPEND_NO_VERIFY) ! call X_pre_setup(Xen,X) ! + ! In the next lines Yambo will VERIFY the em1d database + ! to check if this iq has been alrady done. + ! When self_detect_E_range=.TRUE. however the Xw%er setup is + ! done only in X_os and the VERIFY fails. This is why the procedure + ! must be repeated here (Andrea fixed Bug by C. Attaccalite on October 2008) + ! + if (self_detect_E_range) then + call X_eh_setup(-X%iq(1),X,Xen,Xk,minmax_ehe) + Xw%er=minmax_ehe + endif + ! + ! Build frequency range only if Xw%p was not already allocated + ! like in LifeTimes calculations or when self_detect_E_range=.TRUE. (real axis GW) ! call FREQUENCIES_setup(Xw) ! @@ -127,20 +140,14 @@ integer function X_em1(Xen,Xk,q,X,Xw,APPEND_NO_VERIFY) ! ! CLEAN (1) ! - if (associated(Xw%p)) then - deallocate(Xw%p) - call mem_est("W-p") - endif - if (allocated(rg_index_bg)) deallocate(rg_index_bg) - deallocate(X_poles_tab,bg_npts,cg_pt,cg_index_bg) - call mem_est("X_poles_tab RGi BGn CGp CGi") + call FREQUENCIES_reset(Xw) ! enddo ! ! CLEAN (2) ! call X_alloc('X') - call X_alloc('OptOsc') + call X_alloc('DIP_q_dot_iR') if (.not.APPEND_NO_VERIFY.or.X%iq(2)==q%nibz) call WF_free() call pp_redux_wait ! diff --git a/src/pol_function/X_os.F b/src/pol_function/X_os.F index e84754c2fe..b80e9bdfe3 100644 --- a/src/pol_function/X_os.F +++ b/src/pol_function/X_os.F @@ -38,7 +38,7 @@ subroutine X_os(Xo,iq,fr,Xen,Xk,Xw,X) use R_lattice, ONLY:g_rot,qindx_X,bz_samp,q_norm use wave_func, ONLY:WF_load use memory_m, ONLY:mem_est - use X_m, ONLY:X_t,X_poles,DIP_iq_dot_r,current_iq,X_poles_tab,& + use X_m, ONLY:X_t,X_poles,DIP_q_dot_iR,current_iq,X_poles_tab,& & self_detect_E_range,half_X_mat_only,use_X_RIM implicit none type(levels) :: Xen @@ -228,7 +228,7 @@ subroutine X_os(Xo,iq,fr,Xen,Xk,Xw,X) isc%rhotw(g1)=rhotw_save(g2) enddo if (is>nsym/(i_time_rev+1)) isc%rhotw=conjg(isc%rhotw) - isc%rhotw(1)=-conjg(DIP_iq_dot_r(ic,iv,ikbz,i_spin)) + isc%rhotw(1)=-conjg(DIP_q_dot_iR(ic,iv,ikbz,i_spin)) ! ! ! Filling the upper triangular part of the residual here ! @@ -317,6 +317,7 @@ subroutine X_os(Xo,iq,fr,Xen,Xk,Xw,X) enddo endif ! + ! ! CLEAN !======= ! diff --git a/src/qp/.objects b/src/qp/.objects index 6680fc8e45..f49bd2908a 100644 --- a/src/qp/.objects +++ b/src/qp/.objects @@ -1,6 +1 @@ -GW_REAL_AXIS_objects = - - - - -objs = QP_driver.o QP_newton.o QP_ppa_cohsex.o $(GW_REAL_AXIS_objects) QP_descriptions.o QP_report_and_write.o QP_of.o XCo_driver.o XCo_Hartree_Fock.o XCo_local.o +objs = QP_driver.o QP_newton.o QP_ppa_cohsex.o QP_secant.o QP_real_axis.o QP_W2Sc.o QP_life_transitions.o secant.o bracket.o QP_descriptions.o QP_report_and_write.o QP_Green_Function.o QP_of.o XCo_driver.o XCo_Hartree_Fock.o XCo_local.o diff --git a/src/qp/QP_Green_Function.F b/src/qp/QP_Green_Function.F new file mode 100644 index 0000000000..0590270461 --- /dev/null +++ b/src/qp/QP_Green_Function.F @@ -0,0 +1,149 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine QP_Green_Function(qp,E,GFdb_read_err) + ! + ! This routine calcualtes a modified definition of the GF that + ! yields an always positive spectral function. This is because the + ! Im(G) can be easily post-processed. + ! + ! So regardless of the temperature the G defined here is + ! + ! G = conjg( G_causal ) + ! + use pars, ONLY:SP,pi + use com, ONLY:msg + use units, ONLY:HARTREE + use QP_m, ONLY:QP_t,QP_Sc,QP_Vnl_xc,QP_Vxc,QP_table,& +& QP_n_states,use_GreenF_to_eval_QP,QP_Sc_steps,QP_G_dr,& +& QP_dSc_delta,QP_G_Zoom_treshold,& +& QP_G_amplitude_integral,QP_G_zoom_er,use_GreenF_Zoom + use electrons, ONLY:spin,levels + use frequency, ONLY:w_samp,W_reset + implicit none + ! + type(levels) ::E + type(QP_t) ::qp + integer ::GFdb_read_err + ! + ! Work Space + ! + type(w_samp) ::W + integer ::iqp,ib,ik,i_spin,iw,iw_start,iw_end + real(SP) ::Eo,Amp_Peak,old_min_max_step(2),new_min_max_step(2),W_step + ! + real(SP), external :: Integrate + ! + ! The Green's function + ! + call W_reset(W) + W%n =QP_Sc_steps + W%dr =QP_G_dr + old_min_max_step =(/1000.,0./) + new_min_max_step =(/1000.,0./) + ! + do iqp=1,QP_n_states + ! + ! Frequency range + ! + if (GFdb_read_err<0) then + call FREQUENCIES_Green_Function(iqp,W,E%E) + qp%GreenF_W(iqp,:)=W%p(:) + endif + ! + ! Frequency step + ! + W_step=real(qp%GreenF_W(iqp,2)-qp%GreenF_W(iqp,1)) + if (W_stepold_min_max_step(2)) old_min_max_step(2)=W_step + ! + ! Bare energy + ! + ib =QP_table(iqp,1) + ik =QP_table(iqp,3) + i_spin=spin(QP_table(iqp,:)) + ! + Eo=E%E(ib,ik,i_spin) + if (associated(E%Eo)) Eo=E%Eo(ib,ik,i_spin) + ! + if (GFdb_read_err<0) then + if ( allocated(QP_Vnl_xc) ) then + qp%S_total(iqp,:)=QP_Sc(iqp,:)+QP_Vnl_xc(iqp)-QP_Vxc(iqp) + else + qp%S_total(iqp,:)=QP_Sc(iqp,:) + endif + endif + ! + do iw=1,W%n(1) + ! + qp%GreenF(iqp,iw)=1./(real(qp%GreenF_W(iqp,iw))-Eo-qp%S_total(iqp,iw)) + ! + enddo + ! + ! Spectral Function Integral + ! + QP_G_amplitude_integral(iqp)=Integrate(aimag(qp%GreenF(iqp,:)),real(qp%GreenF_W(iqp,:)),W%n(1))/pi + ! + ! Search for a finer energy range + ! + if (use_GreenF_Zoom.and.GFdb_read_err==0) then + ! + Amp_Peak=maxval(aimag(qp%GreenF(iqp,:))) + iw_start=0 + do iw=1,W%n(1) + if (iw_start==0.and.aimag(qp%GreenF(iqp,iw))>Amp_Peak*QP_G_Zoom_treshold/100.) iw_start=iw + if (iw_start/=0) exit + enddo + iw_end=0 + do iw=W%n(1),1,-1 + if (iw_end==0.and.aimag(qp%GreenF(iqp,iw))>Amp_Peak*QP_G_Zoom_treshold/100.) iw_end=iw + if (iw_end/=0) exit + enddo + QP_G_zoom_er(iqp,:)=(/real(qp%GreenF_W(iqp,iw_start)),real(qp%GreenF_W(iqp,iw_end))/) + qp%GreenF(iqp,:)=(0.,0.) + ! + ! Zoomed Frequency step + ! + W_step=(QP_G_zoom_er(iqp,2)-QP_G_zoom_er(iqp,1))/QP_Sc_steps + if (W_stepnew_min_max_step(2)) new_min_max_step(2)=W_step + ! + call FREQUENCIES_Green_Function(iqp,W,E%E) + qp%GreenF_W(iqp,:)=W%p(:) + ! + endif + ! + if (.not.use_GreenF_to_eval_QP) cycle + ! + ! Extraction of the QP properties directly from the Green's function + !==================================================================== + ! + call QPartilize(W%n(1),qp%GreenF(iqp,:),real(qp%GreenF_W(iqp,:)),qp%E(iqp),qp%Z(iqp),QP_dSc_delta) + ! + enddo + ! + if ( use_GreenF_Zoom .and. GFdb_read_err==0) then + call msg('nr','[Green zoom] Frequency step (current) [ev]:',old_min_max_step*HARTREE) + call msg('r', '[Green zoom] (zoom) [ev]:',new_min_max_step*HARTREE) + endif + ! +end subroutine + diff --git a/src/qp/QP_W2Sc.F b/src/qp/QP_W2Sc.F new file mode 100644 index 0000000000..f5e174e223 --- /dev/null +++ b/src/qp/QP_W2Sc.F @@ -0,0 +1,107 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine QP_W2Sc(iqbz,k,E,Xw,Sc_W) + ! + ! Performs the complex Hilber transofrmation corresponding to + ! + ! \int dw' G(w-w')W(w') + ! + use pars, ONLY:SP + use QP_m, ONLY:QP_t,QP_W,QP_Sc,QP_n_states,QP_G_damp,QP_n_G_bands,& +& QP_solver_state,QP_n_W_freqs,QP_solver,QP_table + use drivers, ONLY:Finite_Tel + use frequency, ONLY:w_samp + use par_proc_m, ONLY:pp_redux_wait,pp_indexes,pp_indexes_reset + use par_indexes_m, ONLY:par_indexes + use functions, ONLY:bose_f,bose_decay + use electrons, ONLY:levels + use R_lattice, ONLY:qindx_S,bz_samp + ! + implicit none + type(bz_samp)::k + type(levels) ::E + type(w_samp) ::Xw,Sc_W(QP_n_states) + integer ::iqbz + ! + ! WorkSpace + ! + integer :: i1,i2,ib,is(2),os(2) + type(pp_indexes):: px + real(SP) :: tsign + complex(SP) :: QP_W_here(QP_n_W_freqs) + complex(SP), allocatable :: dSc(:) + ! + tsign=-1. + if (Finite_Tel) tsign=1. + ! + call pp_indexes_reset(px) + call par_indexes(px,(/QP_n_G_bands(2)/),low_range=(/QP_n_G_bands(1)/)) + ! + do i1=1,QP_n_states + ! + if (allocated(QP_solver_state)) then + if(QP_solver_state(i1)<=0) cycle + endif + ! + allocate(dSc(Sc_W(i1)%n(1))) + ! + dSc=(0._SP,0._SP) + ! + is=(/QP_table(i1,1),QP_table(i1,3)/) ! (nk) QP + os(2)=k%sstar(qindx_S(is(2),iqbz,1),1) ! (nk) intermediate state + ! + do ib=QP_n_G_bands(1),QP_n_G_bands(2) + if (.not.px%element_1D(ib)) cycle + ! + os(1)=ib + ! + ! 1st term: (2-f_os+fbose) + ! + forall(i2=1:Xw%n(1)) QP_W_here(i2)=QP_W(i1,ib,i2)*& +& (2.-E%f(os(1),os(2),1)+bose_f(real(Xw%p(i2))))*& +& bose_decay(real(Xw%p(i2))) + ! + call Kramers_Kronig(QP_W_here,real(Xw%p(:)),QP_n_W_freqs,dSc,& +& real(Sc_W(i1)%p(:))-cmplx(0.,tsign,SP)*aimag(Sc_W(i1)%p(:)),& +& Sc_W(i1)%n(1),& +& E%E(os(1),os(2),1)+tsign*cmplx(0.,QP_G_damp,SP)) + ! + ! 2nd term: (f_os+fbose) + ! + forall(i2=1:Xw%n(1)) QP_W_here(i2)=QP_W(i1,ib,i2)*& +& (E%f(os(1),os(2),1)+bose_f(real(Xw%p(i2))))*& +& bose_decay(real(Xw%p(i2))) + ! + call Kramers_Kronig(-QP_W_here,real(Xw%p(:)),QP_n_W_freqs,dSc,& +& -real(Sc_W(i1)%p(:))+cmplx(0.,1.,SP)*aimag(Sc_W(i1)%p(:)),Sc_W(i1)%n(1),& +& -E%E(os(1),os(2),1)-cmplx(0.,QP_G_damp,SP)) + ! + enddo + ! + call pp_redux_wait(dSc) + forall(i2=1:Sc_W(i1)%n(1)) QP_Sc(i1,i2)=QP_Sc(i1,i2)+dSc(i2) + ! + deallocate(dSc) + ! + enddo + ! +end subroutine diff --git a/src/qp/QP_descriptions.F b/src/qp/QP_descriptions.F index e626084fe2..813f4a72e0 100644 --- a/src/qp/QP_descriptions.F +++ b/src/qp/QP_descriptions.F @@ -49,6 +49,7 @@ subroutine QP_descriptions(qp,X,Xw,Update) integer, external :: QP_state_extract ! if (.not.Update) n_descs_save=qp%n_descs + if ( Update) qp%n_descs =n_descs_save ! if (allocated(RIM_qpg)) then qp%n_descs=qp%n_descs+1 diff --git a/src/qp/QP_driver.F b/src/qp/QP_driver.F index ea064915f7..6431810641 100644 --- a/src/qp/QP_driver.F +++ b/src/qp/QP_driver.F @@ -67,8 +67,6 @@ subroutine QP_driver(X,Xen,Xk,en,k,q,Xw) !============================================================== if (.not.allocated(QP_Vnl_xc).and..not.allocated(QP_Vxc).and.l_el_corr.and..not.l_life) return ! - QP_solver='n' - ! ! Head message !============== call QP_reset(qp) @@ -79,8 +77,24 @@ subroutine QP_driver(X,Xen,Xk,en,k,q,Xw) call section('*','Dyson equation: Newton solver') else if (trim(QP_solver)=='s') then ! + ! When using secant no SC available + ! + GWo_iterations=0 + ! + write (qp%description(1),'(a)') ' GW solver : Secant' + call section('*','Dyson equation: non perturbative secant method') ! + else if (trim(QP_solver)=='g') then + write (qp%description(1),'(a)') ' GW solver : Full Green`s function' + call section('*','Dyson equation: full Green`s function') + call msg('r', '[Green] Sc steps :',QP_Sc_steps) + call msg('r', '[Green] Sc energy range (centered in the bare value) [ev]:',QP_G_er*HARTREE) + call msg('rn','[Green] Sc damping range [ev]:',QP_G_dr*HARTREE) ! + GF_is_causal=Finite_Tel.or.l_ph_corr + ! + else if (.not.l_life) then + return ! endif qp%n_descs=2 @@ -99,6 +113,13 @@ subroutine QP_driver(X,Xen,Xk,en,k,q,Xw) QP_DB_kind=SE_COHSEX ! else + ! + qp%description(2)=' == Real Axis GW ==' + ! + ! Here I am not considering the case where ELPH + GW is used. + ! For this case I need to create a new global KIND. + ! + QP_DB_kind=SE_GoWo ! ! endif @@ -116,6 +137,14 @@ subroutine QP_driver(X,Xen,Xk,en,k,q,Xw) qp%nb=QP_nb qp%n_states=QP_n_states ! + ! In lifetimes calculations the X db may be not + ! present. So I need to define some variables that + ! must be correctly written in the QP_description(s) + ! + if (l_life) then + call X_pre_setup(Xen,X) + if (X%ng_db==0) X%ng_db=X%ng + endif ! ! Local Allocations !=================== @@ -155,9 +184,54 @@ subroutine QP_driver(X,Xen,Xk,en,k,q,Xw) return endif ! + ! Green Functions I/O + !===================== + ! + if (trim(QP_solver)=='g') then + ! + call io_control(ACTION=OP_RD_CL,COM=REP,SEC=(/1,2/),MODE=VERIFY,ID=ID_G) + io_G_err=ioQP('G',qp,ID_G) + ! + if (io_G_err==0.or.io_G_err==IO_NO_BINDING_ERROR) then + ! + if (use_GreenF_to_eval_QP) use_GreenF_to_eval_QP=(.not.io_QP_err==0).and.(io_G_err==0) + if (use_GreenF_Zoom) use_GreenF_Zoom=io_G_err==IO_NO_BINDING_ERROR.and..not.use_GreenF_to_eval_QP + ! + call io_control(ACTION=OP_RD_CL,COM=NONE,SEC=(/1,2,3/),ID=ID_G) + io_G_err=ioQP('G',qp,ID_G) + ! + if (use_GreenF_to_eval_QP) then + call QP_Green_Function(qp,en,0) + call QP_report_and_write(k,qp,en,-1) + call local_free() + return + endif + ! + if (use_GreenF_Zoom) then + allocate(QP_G_zoom_er(QP_n_states,2)) + call mem_est("QP_G_zoom_er",(/size(QP_G_zoom_er)/),(/SP/)) + else + QP_G_Zoom_treshold=0. + endif + call QP_Green_Function(qp,en,0) + if (.not.use_GreenF_Zoom) then + call QP_report_and_write(k,qp,en,0) + call local_free() + return + endif + ! + else + QP_G_Zoom_treshold=0. + endif + ! + endif ! endif ! + ! Updated descriptions (needed in case of Zoomed GFs) + !======================= + ! + call QP_descriptions(qp,X,Xw,.TRUE.) ! call mem_est("QP_Sc",(/size(QP_Sc)/)) ! @@ -179,10 +253,40 @@ subroutine QP_driver(X,Xen,Xk,en,k,q,Xw) endif endif ! + ! SECANT + !-------- + if (trim(QP_solver)=='s') call QP_secant(X,Xen,Xk,en,k,q,qp,Xw) ! ! GREEN`s FUNCTIONS !------------------- ! + if (trim(QP_solver)=='g') then + ! + if (l_el_corr.and.(.not.l_ppa.and..not.l_cohsex)) call QP_real_axis(X,Xen,Xk,en,k,q,qp,Xw,0) + ! + + ! + call QP_Green_Function(qp,en,-1) + ! + call QP_report_and_write(k,qp,en,-1) + call local_free() + return + ! + endif + ! + ! LIFETIMES + ! + if (l_life) then + ! + call QP_real_axis(X,Xen,Xk,en,k,q,qp,Xw,0) + ! + qp%Z=1. + ! + call QP_report_and_write(k,qp,en,-1) + call local_free() + return + ! + endif ! ! Reporting ! diff --git a/src/qp/QP_life_transitions.F b/src/qp/QP_life_transitions.F new file mode 100644 index 0000000000..1fc4a521a5 --- /dev/null +++ b/src/qp/QP_life_transitions.F @@ -0,0 +1,123 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +integer function QP_life_transitions(iqibz_in,E,k,q,X_life_W) + ! + use pars, ONLY:SP + use memory_m, ONLY:mem_est + use QP_m, ONLY:QP_t,QP_n_G_bands,QP_n_states,QP_table,& +& QP_cg_percent + use functions, ONLY:bose_f,h2ht,e2et + use com, ONLY:error + use frequency, ONLY:w_samp,bg_npts,cg_npts,cg_pt,rg_index_bg + use electrons, ONLY:levels + use R_lattice, ONLY:qindx_S,bz_samp + ! + implicit none + type(levels) ::E + type(bz_samp)::q,k + type(w_samp) ::X_life_W + integer ::iqibz_in + ! + ! Work Space + ! + integer :: iqbz,iqibz,ob,ib,ik,ok,ob_max,ob_min,i1,& +& n_of_transitions(q%nibz),i_or,i2 + real(SP) :: F + real(SP),allocatable :: X_freqs(:,:) + real(SP),external :: FREQUENCIES_damping + ! + ! : + ! /:\ iqibz + ! : + ! (ik,ib ) --<--: + ! \_ + ! |\ + ! (ok,ob) + ! + ! 2 conditions: + ! + ! [1] E(ib,ik)-E(ob,ok)>0 & 2-f(ok,ob)+bose(E(ik,ib)-E(ok,ob))>0 + ! [2] E(ib,ik)-E(ob,ok)<0 & f(ok,ob)+bose(E(ok,ob)-E(ik,ib))>0 + ! + QP_life_transitions=0 + ! + do i1=1,2 + ! + n_of_transitions=0 + ob_min=QP_n_G_bands(2) + ob_max=-1 + ! + do iqbz=1,q%nbz + iqibz=q%sstar(iqbz,1) + do i2=1,QP_n_states + ib=QP_table(i2,1) + ik=QP_table(i2,3) + ok=k%sstar(qindx_S(ik,iqbz,1),1) + do ob=QP_n_G_bands(1),QP_n_G_bands(2) + i_or=IOR(e2et((/ib,ik/),(/ob,ok/),E,F),h2ht((/ib,ik/),(/ob,ok/),E,F)) + if (i_or==0) cycle + n_of_transitions(iqibz)=n_of_transitions(iqibz)+1 + if (i1==2) X_freqs(iqibz,n_of_transitions(iqibz))=abs(E%E(ib,ik,1)-E%E(ob,ok,1)) + ob_min=min(ob_min,ob) + ob_max=max(ob_max,ob) + enddo + enddo + enddo + if (i1==2) QP_n_G_bands=(/ob_min,ob_max/) + if (i1==2) exit + ! + if (maxval(n_of_transitions)==0) then + call error('All virtual transitions are null. Change QP states') + else if (any(n_of_transitions==0)) then + call error('One or more Q-point virtual transitions are null. Change QP states') + endif + ! + allocate(X_freqs(q%nibz,maxval(n_of_transitions))) + X_freqs=0. + enddo + QP_n_G_bands=(/ob_min,ob_max/) + ! + QP_life_transitions=maxval(n_of_transitions) + ! + if (iqibz_in<0) return + ! + ! Here I reduce the frequencies + ! + call FREQUENCIES_coarse_grid('Life',X_freqs(iqibz_in,:n_of_transitions(iqibz_in)),& +& n_of_transitions(iqibz_in),QP_cg_percent) + X_life_W%n=cg_npts + ! + X_life_W%er=(/minval(cg_pt),maxval(cg_pt)/) + ! + allocate(X_life_W%p(cg_npts)) + do i1=1,cg_npts + ! + X_life_W%p(i1)=cg_pt(i1)+FREQUENCIES_damping(X_life_W,cg_pt(i1))*cmplx(0.,1.,SP) + ! + enddo + ! + ! Clean + ! + deallocate(bg_npts,cg_pt,rg_index_bg) + call mem_est("BGn CGp RGi") + ! +end function diff --git a/src/qp/QP_newton.F b/src/qp/QP_newton.F index a879643772..aad56fd0b9 100644 --- a/src/qp/QP_newton.F +++ b/src/qp/QP_newton.F @@ -72,6 +72,7 @@ subroutine QP_newton(X,Xen,Xk,en,k,q,qp,Xw) call QP_ppa_cohsex(X,Xk,en,k,q,qp,Xw,iter) else ! + call QP_real_axis(X,Xen,Xk,en,k,q,qp,Xw,iter) ! endif else diff --git a/src/qp/QP_of.F b/src/qp/QP_of.F index 6975e84462..bcd8e4eb24 100644 --- a/src/qp/QP_of.F +++ b/src/qp/QP_of.F @@ -152,6 +152,52 @@ subroutine QP_of(qp,en,QPdb_read_err) enddo call of_open_close(file_name) ! + else if (trim(QP_solver)=='g') then + ! + call QP_states_simmetrize(en,state_is_2do=state_is_2do) + ! + call set_real_printed_length(f_length=12,g_length=12) + ! + do i1=1,QP_n_states + ! + write (G_Sc_name,'(2(a,i3.3))') 'G_Sc_band_',QP_table(i1,1),'_k_',QP_table(i1,3) + call of_open_close(G_Sc_name,'ot') + call msg('o G_Sc','# GW [Green`s function & Self-Energy]') + call msg('o G_Sc','#') + do i2=2,qp%n_descs + call msg('o G_Sc','# ',trim(qp%description(i2)),INDENT=0) + enddo + call msg('o G_Sc','#') + call msg('o G_Sc','# Spectral function at this point integrates to:',QP_G_amplitude_integral(i1),& +& INDENT=0) + call msg('o G_Sc','#') + titles(1)='Energy' + titles(2:3)=(/' Re[G] ','|Im(G)|'/) + titles(4)= 'Re(S_tot)' + titles(5:6)=(/'|Im(S_c)|',' Re(S_c) '/) + if (report_Vnlxc) then + call msg('o G_Sc','#',titles(:6),INDENT=0,USE_TABS=.TRUE.) + else + call msg('o G_Sc','#',titles(:5),INDENT=0,USE_TABS=.TRUE.) + endif + call msg('o G_Sc','#') + do i2=1,QP_Sc_steps + values=0. + values(1) =real(qp%GreenF_W(i1,i2))*HARTREE + values(2:3)=(/real(qp%GreenF(i1,i2)),aimag(qp%GreenF(i1,i2))/)/HARTREE + values(4:5)=(/real(qp%S_total(i1,i2)),aimag(qp%S_total(i1,i2))/)*HARTREE + values(6)=values(4) + if (report_Vnlxc) then + values(6)=values(6)-real(QP_Vnl_xc(i1)-QP_Vxc(i1))*HARTREE + call msg('o G_Sc','',values(:6),INDENT=-2,USE_TABS=.TRUE.) + else + call msg('o G_Sc','',values(:5),INDENT=-2,USE_TABS=.TRUE.) + endif + enddo + call of_open_close(G_Sc_name) + enddo + ! + call set_real_printed_length() ! endif ! diff --git a/src/qp/QP_real_axis.F b/src/qp/QP_real_axis.F new file mode 100644 index 0000000000..fea5e9a570 --- /dev/null +++ b/src/qp/QP_real_axis.F @@ -0,0 +1,439 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine QP_real_axis(X,Xen,Xk,en,k,q,qp,Xw,GW_iter) + ! + ! This routine calculates + ! + ! [1] Real axis GW S.E. + ! [2] Electronic Lifetimes + ! + use pars, ONLY:SP,schlen,pi,IP + use units, ONLY:HARTREE + use drivers, ONLY:l_life + use electrons, ONLY:levels,spin + use frequency, ONLY:w_samp,cg_index_bg,cg_npts,W_reset + use timing, ONLY:live_timing + use com, ONLY:msg + use par_proc_m, ONLY:pp_redux_wait,pp_indexes,myid,pp_indexes_reset + use par_indexes_m, ONLY:par_indexes + use collision, ONLY:ggwinfo,collision_reset + use IO_m, ONLY:io_control,REP,VERIFY,NONE,RD,RD_CL,OP_RD,& +& OP_WR_CL,OP_APP_WR_CL,RD_CL_IF_END,OP_RD_CL + use QP_m, ONLY:QP_t,QP_n_G_bands,QP_nb,QP_dSc_steps,& +& QP_Sc,QP_n_states,QP_G_damp,QP_table,QP_dSc_delta,& +& QP_W,QP_solver,QP_W_er,QP_W_dr,QP_n_W_freqs,& +& QP_G_dr,QP_Sc_steps + use X_m, ONLY:X_alloc,X_mat,X_t + use functions, ONLY:e2et,h2ht,bose_decay + use memory_m, ONLY:mem_est + use R_lattice, ONLY:qindx_S,bz_samp + use D_lattice, ONLY:nsym,i_space_inv,i_time_rev + use wave_func, ONLY:WF_load,WF_free + use stderr, ONLY:intc + ! + implicit none + type(levels) ::en,Xen + type(bz_samp)::Xk,k,q + type(X_t) ::X + type(QP_t) ::qp + type(w_samp) ::Xw + integer ::GW_iter + ! + ! WorkSpace + ! + type(ggwinfo) ::isc,iscp + type(pp_indexes) ::px + type(w_samp) ::Sc_W(qp%n_states),X_life_W(q%nibz) + integer :: i1,i2,i3,i_or,iqbz,iqibz,ib,i_err,iq_to_start,iqs,& +& iv4(4),io_err,XID,WID,QP_n_W_freqs_redux + complex(SP) :: lrhotw(X%ng),X_mat_ws(X%ng,X%ng) +#if defined _DOUBLE + complex(SP) :: zdotc,zdotu +#else + complex(SP) :: cdotc,cdotu +#endif + real(SP) ::life_Fe,life_Fh + integer, allocatable::life_W_table(:,:) + integer, external ::ioX,QP_state_extract,ioQP,X_em1 + ! + integer, external ::QP_life_transitions + ! + character(schlen) ::ch,Sc_name + logical ::X_is_TR_rotated + ! + ! Presets + ! + call collision_reset(isc) + call collision_reset(iscp) + call pp_indexes_reset(px) + if (l_life) then + ! + do i1=1,q%nibz + call W_reset(X_life_W(i1)) + X_life_W(i1)%dr=Xw%dr + enddo + ! + else + do i1=1,qp%n_states + call W_reset(Sc_W(i1)) + enddo + endif + ! + call k_expand(k) + ! + ! ALLOCATION (W/Sc) + ! + if (.not.l_life) then + QP_n_W_freqs=Xw%n(1) + QP_n_W_freqs_redux=Xw%n(1) + allocate(QP_W(QP_n_states,QP_n_G_bands(2),QP_n_W_freqs),stat=i_err) + call mem_est("QP_W",(/size(QP_W)/),errors=(/i_err/)) + endif + ! + !Sc Energy points (1 type each QP state !) + ! + if (trim(QP_solver)=='n') then + do i1=1,qp%n_states + Sc_W(i1)%n=QP_dSc_steps + allocate(Sc_W(i1)%p(Sc_W(i1)%n(1))) + forall (i2=1:QP_dSc_steps) Sc_W(i1)%p(i2)=qp%E_bare(i1)+(i2-1)*QP_dSc_delta + enddo + else if (trim(QP_solver)=='g') then + ! + ! I need to put to 0. the G damping embodied in QP_W2Sc + ! + QP_G_damp=0. + ! + do i1=1,qp%n_states + Sc_W(i1)%n =QP_Sc_steps + Sc_W(i1)%dr=QP_G_dr + call FREQUENCIES_Green_Function(i1,Sc_W(i1),en%E) + QP_Sc_steps=Sc_W(i1)%n(1) + enddo + endif + ! + Sc_name="G"//trim(intc(GW_iter)) + ! + if (GW_iter==0) call section('+',trim(Sc_name)//"W0 on the real axis") + if (GW_iter> 0) call section('=',trim(Sc_name)//"W0 on the real axis") + ! + if (.not.l_life) call msg('r', '[GW] Bands range :',QP_n_G_bands) + call msg('r', '[GW] G damping [ev]:',QP_G_damp*HARTREE) + call msg('r','') + iv4=(/1,1,0,0/) + do while(QP_state_extract(iv4)>0) + write (ch,'(4(a,i3.3))') 'QP @ K ',iv4(1),' - ',iv4(2),' : b ',iv4(3),' - ',iv4(4) + call msg('r',trim(ch)) + enddo + call msg('r','') + ! + ! W DB + ! + call io_control(ACTION=OP_RD,COM=REP,SEC=(/1,2/),MODE=VERIFY,ID=WID) + io_err=ioQP('W',qp,WID) + ! + iq_to_start=1 + if (io_err>0) iq_to_start=io_err + ! + if (io_err==0) then + if (QP_solver=="n".or.QP_solver=="g") then + QP_Sc=cmplx(0.,0.,SP) + call live_timing(trim(Sc_name)//'W0',q%nbz) + do iqbz=1,q%nbz + if (iqbz< q%nbz) call io_control(ACTION=RD,COM=NONE,SEC=(/2+iqbz/),ID=WID) + if (iqbz==q%nbz) call io_control(ACTION=RD_CL,COM=NONE,SEC=(/2+iqbz/),ID=WID) + io_err=ioQP('W',qp,WID) + ! + ! Xw%p is stored in ioX, not in ioW ! + ! + Xw%er=QP_W_er + Xw%dr=QP_W_dr + Xw%n=QP_n_W_freqs + call FREQUENCIES_setup(Xw) + ! + call QP_W2Sc(iqbz,k,en,Xw,Sc_W) + call live_timing(steps=1) + enddo + deallocate(QP_W) + call mem_est("QP_W") + call live_timing() + return + endif + ! + if (QP_solver=="s") return + ! + endif + ! + ! Lifetimes Transitions + ! + if (l_life) then + ! + call section('=','Lifetimes Transitions Selector') + ! + QP_n_W_freqs=QP_life_transitions(-1,en,k,q,X_life_W(1)) + allocate(life_W_table(q%nibz,QP_n_W_freqs)) + call mem_est("life_W_table",(/size(life_W_table)/),(/IP/)) + life_W_table=0 + QP_n_W_freqs_redux=-1 + do iqibz=1,q%nibz + ! + X_life_W(iqibz)%dr=Xw%dr + ! + i1=QP_life_transitions(iqibz,en,k,q,X_life_W(iqibz)) + life_W_table(iqibz,:size(cg_index_bg))=cg_index_bg(:) + ! + QP_n_W_freqs_redux=max(QP_n_W_freqs_redux,cg_npts) + deallocate(cg_index_bg) + call mem_est("CGi") + ! + enddo + do iqibz=1,q%nibz + X%iq=iqibz + i_err=X_em1(Xen,Xk,q,X,X_life_W(iqibz),iqibz>1) + if (i_err==0) exit + enddo + endif + ! + ! WFs + ! + call WF_load(X%ng,maxval(qindx_S(:,:,2)),& +& (/QP_n_G_bands(1),max(QP_n_G_bands(2),QP_nb)/),(/1,k%nibz/),title='-Sc') + ! + ! Test the spatial Inversion + ! + call WF_spatial_invertion(en,Xk) + ! + ! ALLOCATION (isc) + ! + call X_alloc('X',(/X%ng,X%ng,QP_n_W_freqs_redux/)) + ! + allocate(isc%gamp(X%ng,X%ng),isc%rhotw(X%ng),iscp%rhotw(X%ng),stat=i_err) + call mem_est("ISC-GAMP",(/X%ng**2+2*X%ng/),errors=(/i_err/)) + ! + ! MPI index + call par_indexes(px,(/QP_n_states,QP_n_G_bands(2)/),low_range=(/1,QP_n_G_bands(1)/)) + ! + QP_Sc=cmplx(0.,0.,SP) + isc%ngrho=X%ng + isc%iqref=0 + ! + ! Main Loop + !=========== + ! + call pp_redux_wait + call live_timing(trim(Sc_name)//'W0',px%n_of_element_1D(myid+1)*(q%nbz-iq_to_start+1)) + ! + do iqbz=iq_to_start,q%nbz + ! + if (.not.l_life) QP_W=(0.,0.) + ! + isc%qs(2:)=(/q%sstar(iqbz,1),q%sstar(iqbz,2)/) + iqibz=isc%qs(2) + iqs =isc%qs(3) + ! + if (iqibz/=isc%iqref) then + call scatterGamp(isc,'c') + ! + ! X I/O + ! + if (iqbz==iq_to_start) call io_control(ACTION=OP_RD,COM=NONE,& +& SEC=(/1,2,2*iqibz+1/),ID=XID) + if (iqbz/=iq_to_start) call io_control(ACTION=RD_CL_IF_END,COM=NONE,& +& SEC=(/2*iqibz,2*iqibz+1/),ID=XID) + if (q%nbz==1) call io_control(ACTION=OP_RD_CL,COM=NONE,SEC=(/1,2,3/),ID=XID) + ! + if (l_life) then + io_err=ioX(X,X_life_W(iqibz),XID) + call X_delta_part(X,X_life_W(iqibz),isc%gamp) + else + ! + io_err=ioX(X,Xw,XID) + ! + ! Xw%er/dr/n are not known here and are not read from ioX so + ! I need to redefine them in terms of Xw%p (read from ioX) + ! + Xw%er=(/real(Xw%p(1)),real(Xw%p(Xw%n(1)))/) + Xw%dr=(/aimag(Xw%p(1)),aimag(Xw%p(Xw%n(1)))/) + call X_delta_part(X,Xw,isc%gamp) + ! + endif + ! + if (l_life) QP_n_W_freqs=0 + ! + X_is_TR_rotated=.false. + ! + endif + ! + if (iqs>nsym/(i_time_rev+1) .and. i_space_inv == 0 .and..not.X_is_TR_rotated) then + X_is_TR_rotated=.true. + do i3=1,Xw%n(1) + forall(i1=1:X%ng,i2=1:X%ng) X_mat_ws(i2,i1)=X_mat(i1,i2,i3) + X_mat(:,:,i3)=X_mat_ws(:,:) + enddo + endif + ! + do i1=1,QP_n_states + ! + isc%is=(/QP_table(i1,1),QP_table(i1,3),1,spin(QP_table(i1,:))/) + isc%os(2:)=(/k%sstar(qindx_S(isc%is(2),iqbz,1),:),spin(QP_table(i1,:))/) + iscp%is=(/QP_table(i1,2),QP_table(i1,3),1,spin(QP_table(i1,:))/) + ! + isc%qs(1)=qindx_S(QP_table(i1,3),iqbz,2) + iscp%qs=isc%qs + ! + do ib=QP_n_G_bands(1),QP_n_G_bands(2) + ! + isc%os(1)=ib + i_or=IOR(e2et(isc%is(:2),isc%os(:2),en,life_Fe),& +& h2ht(isc%is(:2),isc%os(:2),en,life_Fh)) + if (l_life) QP_n_W_freqs=QP_n_W_freqs+i_or + ! + if (.not.px%element_2D(i1,ib)) cycle + ! + iscp%os=isc%os + call live_timing(steps=1) + ! + if (l_life.and.i_or==0) cycle + ! + call scatterBamp(isc) + iscp%rhotw=isc%rhotw + ! + if (any(isc%is/=iscp%is)) call scatterBamp(iscp) + ! + if (l_life) then + ! + i2=life_W_table(iqibz,QP_n_W_freqs) + do i3=1,X%ng +#if defined _DOUBLE + lrhotw(i3)=zdotu(X%ng,isc%rhotw,1,X_mat(1,i3,i2),1) +#else + lrhotw(i3)=cdotu(X%ng,isc%rhotw,1,X_mat(1,i3,i2),1) +#endif + enddo + ! + ! To compensate the Tel/w divergence of the Bose function at finite + ! Tel I multiply the X_mat function by [w/(Tel*Bose_E_cut)]^2 + ! + ! Note that the same procedure is applied in QP_W2Sc when S_c is + ! calculated + ! +#if defined _DOUBLE + qp%E(i1)=qp%E(i1)-(0.,2.)*pi*(life_Fe+life_Fh)*& +& bose_decay( real(X_life_W(iqibz)%p(i2)) )*& +& zdotc(X%ng,iscp%rhotw,1,lrhotw,1) +#else + qp%E(i1)=qp%E(i1)-(0.,2.)*pi*(life_Fe+life_Fh)*& +& bose_decay( real(X_life_W(iqibz)%p(i2)) )*& +& cdotc(X%ng,iscp%rhotw,1,lrhotw,1) +#endif + ! + else ! .not.l_life + ! + do i2=1,QP_n_W_freqs + ! + do i3=1,X%ng +#if defined _DOUBLE + lrhotw(i3)=zdotu(X%ng,isc%rhotw,1,X_mat(1,i3,i2),1) +#else + lrhotw(i3)=cdotu(X%ng,isc%rhotw,1,X_mat(1,i3,i2),1) +#endif + enddo + ! +#if defined _DOUBLE + QP_W(i1,ib,i2)=-2.*zdotc(X%ng,iscp%rhotw,1,lrhotw,1) +#else + QP_W(i1,ib,i2)=-2.*cdotc(X%ng,iscp%rhotw,1,lrhotw,1) +#endif + enddo + ! + endif + ! + enddo + enddo + ! + call pp_redux_wait + call pp_redux_wait(QP_W) + ! + ! To be written in W + ! + if (.not.l_life) then + ! + QP_W_er=Xw%er + QP_W_dr=Xw%dr + QP_n_W_freqs=Xw%n(1) + ! + ! Only if I am doing secant solver of Dyson I store on disk + ! the W m.e. + ! + if (trim(QP_solver)=='n'.or.trim(QP_solver)=='g') then + call QP_W2Sc(iqbz,k,en,Xw,Sc_W) + else + ! + if (iqbz==1) then + call io_control(ACTION=OP_WR_CL,COM=REP,SEC=(/1,2,3/),ID=WID) + else + call io_control(ACTION=OP_APP_WR_CL,COM=REP,SEC=(/2+iqbz/),ID=WID) + endif + io_err=ioQP('W',qp,WID) + ! + endif + ! + endif + ! + enddo + ! + call live_timing() + call pp_redux_wait() + ! + ! CLEAN + ! + deallocate(isc%gamp,isc%rhotw,iscp%rhotw) + call mem_est("ISC-GAMP") + ! + if (l_life) then + deallocate(life_W_table) + call mem_est("life_W_table") + call pp_redux_wait(qp%E) + else + if (trim(QP_solver)=='n'.or.trim(QP_solver)=='g') then + deallocate(QP_W) + call mem_est("QP_W") + endif + endif + if (l_life) then + do i1=1,q%nibz + call W_reset(X_life_W(i1)) + enddo + else + do i1=1,qp%n_states + call W_reset(Sc_W(i1)) + enddo + endif + ! + call X_alloc('X') + call WF_free() + call collision_reset(isc) + call collision_reset(iscp) + call pp_indexes_reset(px) + call pp_redux_wait() + ! +end subroutine diff --git a/src/qp/QP_report_and_write.F b/src/qp/QP_report_and_write.F index 2ed2cbc7c7..a8abd61e4a 100644 --- a/src/qp/QP_report_and_write.F +++ b/src/qp/QP_report_and_write.F @@ -158,6 +158,29 @@ subroutine QP_report_and_write(k,qp,en,QPdb_read_err) endif ! ! + if (l_life) then + ! + ! Note the definition of Lifetime is + ! + ! tau = 1/(2 pi Sigma_2) => tau[fs] = 1/(2 pi Sigma_2[H]) / HA2FSm1 + ! + ! now HA2FSm1 = HARTREE*241.798474*10^-3 = > tau[fs]= 1/( 2 Sigma_2[H] HARTREE ) * + ! 10^3/( pi 241.798474 ) + ! + ! SHADOWS + ! ======= + ! + n_shadows=4 + shadow(:2)=(/qp%E_bare(iqp)*HARTREE,real(qp%Z(iqp))/) + shadow(3)=aimag(qp%E(iqp))*HARTREE/1.E-3 + shadow(4)=0. + if (shadow(3)>0.) shadow(4)=1./(HA2FSm1*2*pi*aimag(qp%E(iqp))) + shadow_string(2)=' Z=' + shadow_fmt(2) ='(a,f4.2)' + shadow_string(3)=' Gm=' + shadow_string(4)=' Gf=' + endif + ! ! Write to Report File !====================== ! diff --git a/src/qp/QP_secant.F b/src/qp/QP_secant.F new file mode 100644 index 0000000000..caa9f42801 --- /dev/null +++ b/src/qp/QP_secant.F @@ -0,0 +1,227 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine QP_secant(X,Xen,Xk,en,k,q,qp,Xw) + ! + !Procedure for non-perturbative solution of Dyson equation. + !Here it`s solved real part of Dyson Eq. with a fixed + !shift from real axis + ! + ! E1 = Eo + Sx - Vxc + Re[Sc(E1)] (1) + ! E1 is REAL + ! + ! Re(E) = E1 + Im[Sc(E1)] + Re[Sc(E)] + ! Im(E) = Im[Sc(E)] + ! + !After convergence is reached Sc(w) is analtically extendend + !outside real axe using first order Taylor expansion around E1. + !In this way it`s found the solution of Complex Equation + ! + ! E = E1 + 1/(1-dSc(E1)) Im[Sc(E1)] + ! + !The procedure for finding root of (1) is divided into two steps: + ! + ! 1. BRACKETING: starting for an initial guess of the + ! QP-energy [Eqp] (given by E_lda) the routine gives a range + ! of energy where the solution is expected to be + ! + ! 2. SECANT SOLUTION: see "Numerical Recipes (Fortran version)", Pag. 248 + ! + use pars, ONLY:SP,IP + use units, ONLY:HARTREE + use memory_m, ONLY:mem_est + use timing, ONLY:live_timing + use electrons, ONLY:levels + use R_lattice, ONLY:bz_samp + use frequency, ONLY:w_samp,W_reset + use X_m, ONLY:X_t + use par_proc_m, ONLY:pp_redux_wait + use QP_m, ONLY:QP_t,QP_Vxc,QP_Vnl_xc,QP_Sc,QP_n_W_freqs,QP_W,QP_n_states,& +& QP_solver_state,QP_dSc_steps,QP_dSc,& +& QP_W_dr,QP_W_er,QP_dSc_delta + use IO_m, ONLY:io_control,NONE,RD,RD_CL,OP_RD + ! + implicit none + type(levels)::en,Xen + type(bz_samp)::Xk,k,q + type(X_t):: X + type(QP_t):: qp + type(w_samp)::Xw + ! + ! Work Space + ! + integer :: i1,i2,ID,io_err,iqbz + type(w_samp) :: Sc_W(qp%n_states) + real(SP) :: Eqp_raxis(qp%n_states,2),dQP_raxis(qp%n_states,2) + complex(SP) :: Z(QP_dSc_steps-1),Eqp(QP_dSc_steps-1) + real(SP),parameter :: root_acc=0.001/HARTREE + ! + integer, external :: ioQP,secant,bracket + ! + ! + !I need to store on disk the m. elemnts + ! + call QP_real_axis(X,Xen,Xk,en,k,q,qp,Xw,0) + ! + !I prepare the Sc energy type + ! + do i1=1,qp%n_states + call W_reset(Sc_W(i1)) + Sc_W(i1)%n=1 + allocate(Sc_W(i1)%p(1)) + enddo + ! + !SOLVER STATE + !============= + ! + ! 1 : updating Eqp_raxis(?,1) + ! 2 : updating Eqp_raxis(?,2) + ! + ! <0: Eqp_raxis(?,iabs(SOLVER_STATE)) converged + ! + allocate(QP_solver_state(qp%n_states)) + call mem_est("QP_solver_state",(/qp%n_states/),(/IP/)) + ! + !Brackets + ! + forall(i1=1:qp%n_states) Eqp_raxis(i1,1)=qp%E_bare(i1) + forall(i1=1:qp%n_states) Eqp_raxis(i1,2)=Eqp_raxis(i1,1)+1./HARTREE + ! + QP_solver_state=1 + call W2Sc_local_call() + ! + QP_solver_state=2 + call W2Sc_local_call() + ! + call live_timing('[SECANT] Brackets',QP_n_states,SERIAL=.true.) + ! + do while(any(QP_solver_state/=0)) + do i1=1,qp%n_states + i2=QP_solver_state(i1) + QP_solver_state(i1)=bracket(Eqp_raxis(i1,1),dQP_raxis(i1,1),& +& Eqp_raxis(i1,2),dQP_raxis(i1,2)) + if (i2/=0.and.QP_solver_state(i1)==0) call live_timing(steps=1) + enddo + call W2Sc_local_call() + enddo + call live_timing() + call pp_redux_wait() + ! + ! Root + ! + QP_solver_state=1 + call live_timing('[SECANT] Root(s)',QP_n_states,SERIAL=.true.) + do while(any(QP_solver_state>=0)) + do i1=1,qp%n_states + i2=QP_solver_state(i1) + QP_solver_state(i1)=secant(Eqp_raxis(i1,1),dQP_raxis(i1,1),& +& Eqp_raxis(i1,2),dQP_raxis(i1,2),root_acc) + if (i2>0.and.QP_solver_state(i1)<0) call live_timing(steps=1) + enddo + call W2Sc_local_call() + enddo + call live_timing() + call pp_redux_wait() + ! + ! Analytic continuation + ! + call live_timing('[SECANT] Analytic continuation',q%nbz,SERIAL=.true.) + do i1=1,qp%n_states + Eqp_raxis(i1,1)=Eqp_raxis(i1,iabs(QP_solver_state(i1))) + call W_reset(Sc_W(i1)) + Sc_W(i1)%n=QP_dSc_steps + allocate(Sc_W(i1)%p(QP_dSc_steps)) + forall (i2=1:QP_dSc_steps) Sc_W(i1)%p(i2)=Eqp_raxis(i1,1)+(i2-1)*QP_dSc_delta + enddo + deallocate(QP_solver_state) + call mem_est("QP_solver_state") + ! + call W2Sc_local_call() + ! + do i1=1,QP_n_states + do i2=1,QP_dSc_steps-1 + QP_dSc(i1,i2)=(QP_Sc(i1,i2+1)-QP_Sc(i1,i2))/QP_dSc_delta + Z(i2)=1./(1.-QP_dSc(i1,i2)) + Eqp(i2)=Eqp_raxis(i1,1)+cmplx(0.,real(Z(i2))*aimag(QP_Sc(i1,1)),SP) + enddo + qp%E(i1)=Eqp(1) + qp%Z(i1)=Z(1) + enddo + call live_timing() + ! + !CLEAN + ! + do i1=1,QP_n_states + call W_reset(Sc_W(i1)) + enddo + deallocate(QP_W) + call mem_est("QP_W") + call pp_redux_wait() + ! + contains + ! + subroutine W2Sc_local_call() + ! + ! Here I call the QP_W2Sc to routine to calculate the Sc SE + ! corresponding to Eqp_raxis(i1,QP_solver_state(i1)) energies + ! + ! The driver QP_solver_state tells if I am updating the 1st or the 2nd + ! approx to the final QP energy. + ! + if (allocated(QP_solver_state)) then + do i1=1,QP_n_states + if (QP_solver_state(i1)>0) Sc_W(i1)%p(1)=Eqp_raxis(i1,QP_solver_state(i1)) + enddo + endif + QP_Sc=(0._SP,0._SP) + do iqbz=1,q%nbz + ! + if (iqbz==1) call io_control(ACTION=OP_RD,COM=NONE,SEC=(/1,2,3/),ID=ID) + if (iqbz> 1.and.iqbz0) dQP_raxis(i1,QP_solver_state(i1))=& +& ( qp%E_bare(i1)+real(QP_Sc(i1,1))+QP_Vnl_xc(i1)-QP_Vxc(i1) )-& +& Eqp_raxis(i1,QP_solver_state(i1)) + enddo + endif + ! + end subroutine + ! +end subroutine diff --git a/src/qp/XCo_local.F b/src/qp/XCo_local.F index 8624594a96..f5e9320d4f 100644 --- a/src/qp/XCo_local.F +++ b/src/qp/XCo_local.F @@ -63,7 +63,7 @@ subroutine XCo_local(E,Xk) l_call_the_driver=.TRUE. ! if (l_call_the_driver) then - call xc_lda_driver(E,Xk,WF_KIND,WF_xc_functional,1) + call XC_potential_driver(E,Xk,WF_KIND,WF_xc_functional,1) endif ! ! diff --git a/src/qp/bracket.F b/src/qp/bracket.F new file mode 100644 index 0000000000..f3ca1a0356 --- /dev/null +++ b/src/qp/bracket.F @@ -0,0 +1,43 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +integer function bracket(x1,f1,x2,f2) + ! + use pars,ONLY: SP + ! + implicit none + real(SP) :: f1,f2,x1,x2 + ! + ! Work Space + ! + real(SP),parameter :: precision=1.E-5,factor=1.6 + ! + bracket=0 + if (f1*f2G) diff --git a/src/xc_functionals/.objects b/src/xc_functionals/.objects index 51bea3e0f0..38481845df 100644 --- a/src/xc_functionals/.objects +++ b/src/xc_functionals/.objects @@ -1 +1,11 @@ -objs = mod_xc_costants.o xc_rpa_kp.o xcspol.o xcpzca.o xcwign.o xchelu.o xcxalp.o el_density.o xc_lda_driver.o +LIBXC_objects = XC_lda_driver.o +#if defined _BOLTZMANN +BOLTZMANN_objects = el_current.o +#endif +#if defined _LIBXC +LIBXC_objects = XC_libxc_driver.o +#endif +objs = mod_xc_costants.o xc_rpa_kp.o \ + xcspol.o xcpzca.o xcwign.o xchelu.o xcxalp.o \ + spin_density.o el_density.o XC_potential_driver.o\ + $(LIBXC_objects) $(BOLTZMANN_objects) diff --git a/src/xc_functionals/XC_lda_driver.F b/src/xc_functionals/XC_lda_driver.F new file mode 100644 index 0000000000..4807151f8f --- /dev/null +++ b/src/xc_functionals/XC_lda_driver.F @@ -0,0 +1,175 @@ +! +! Copyright (C) 2000-2010 A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine XC_lda_driver(en,Xk,KIND,FUNCTIONAL,ORDER) + ! + ! These sources have been slightly adpated from Abinit 4.6.5. + ! + ! xcspol: XC_LDA_C_PW + ! LDA or LSD, Teter Pade parametrization (4/93, published in S. Goedecker, + ! M. Teter, J. Huetter, Phys.Rev.B54, 1703 (1996)), + ! which reproduces Perdew-Wang (which reproduces Ceperley-Alder!). + ! + ! xcpzca: XC_LDA_C_PZ + ! LDA, Perdew-Zunger-Ceperley-Alder (no spin-polarization) + ! + ! xcwign: XC_LDA_C_WIGNER + ! LDA, Wigner functional (no spin-polarization) + ! + ! ! xchelu: XC_LDA_C_HL + ! LDA, Hedin-Lundqvist functional (no spin-polarization) + ! + ! xcxalp: XC_LDA_C_XALPHA + ! LDA, "X-alpha" functional (no spin-polarization) + ! + ! Note that, in order to get LDA exchange and correlation separated + ! + ! Exx = -efac*rsm1 + ! Vxx = -vfac*rsm1 + ! + ! in xcpzca xcwign xchelu xcxalp. While in xcspol xctetr there is not + ! such a type of distinction. + ! + ! Andrea (9/8/06) + ! + ! In case of spin polarized systems *REMEMBER* that: + ! + ! whatever(1) == whatever UP + ! whatever(2) == whatever DOWN + ! + ! Yambo have been considering the wrong order up to revision 335 (November 2008) + ! + use pars, ONLY:SP,DP,pi + use R_lattice, ONLY:bz_samp + use D_lattice, ONLY:DL_vol + use com, ONLY:msg,warning + use electrons, ONLY:levels,n_spin + use FFT_m, ONLY:fft_size + use xc_functionals,ONLY:XC_EXCHANGE,XC_EXCHANGE_CORRELATION,XC_CORRELATION,& +& E_xc,F_xc,V_xc,XC_LDA_X,& +& XC_LDA_C_PW,xc_string,XC_LDA_C_PZ,& +& XC_LDA_C_WIGNER,XC_LDA_C_HL,XC_LDA_C_XALPHA,& +& XC_LDA_C_KP,two_spin_density,NONE + implicit none + type(levels) ::en + type(bz_samp)::Xk + integer ::KIND,FUNCTIONAL,ORDER + ! + ! Work Space + ! + integer :: i1 + ! + ! XC routines arrays + ! + real(SP) :: rho_sp(fft_size) + real(DP) :: rspts(fft_size),zeta(fft_size),rho(fft_size) + ! + ! note that multidimensional dvxc & vxc is allowed only in xcspol + ! + real(DP) :: dvxc(fft_size,3),exc(fft_size),vxc(fft_size,n_spin) + real(DP) :: dvx(fft_size,3), ex(fft_size), vx(fft_size,n_spin) + ! + ! Init + ! + zeta=0._DP + exc=0._DP + vxc=0._DP + dvxc=0._DP + ex=0._DP + vx=0._DP + dvx=0._DP + ! + ! Electronic density & Rs + ! + call el_density(en,Xk,rho_sp,.false.) + rho=rho_sp*real(fft_size,DP)/real(DL_vol,DP) + forall (i1=1:fft_size) rspts(i1)=(3._DP/4._DP/pi/rho(i1))**(1._DP/3._DP) + ! + ! + ! Hartree + ! + if (KIND==NONE) then + if (order==0) E_xc=0. + if (order==1) V_xc=0. + if (order==2) F_xc=0. + return + endif + ! + ! EXX = X_alpha when alpha=3/4. HF when alpha=1. + ! + if (KIND/=XC_EXCHANGE_CORRELATION.or.FUNCTIONAL==XC_LDA_X) then + call xcxalp(1._DP,ex,fft_size,n_spin,ORDER,rspts,zeta,vx,dvx) + if (KIND==XC_EXCHANGE.or.FUNCTIONAL==XC_LDA_X) then + exc=ex + vxc=vx + dvxc=dvx + endif + endif + ! + ! CORRELATION + ! + if (KIND/=XC_EXCHANGE) then + ! + select case (FUNCTIONAL) + ! + case (XC_LDA_C_PW) + call xcspol(dvxc,exc,fft_size,n_spin,ORDER,rspts,vxc,zeta) + ! + case (XC_LDA_C_PZ) + call xcpzca(dvxc(:,1),exc,fft_size,ORDER,rho,rspts,vxc(:,1)) + ! + case (XC_LDA_C_WIGNER) + call xcwign(dvxc(:,1),exc,fft_size,ORDER,rho,rspts,vxc(:,1)) + ! + case (XC_LDA_C_HL) + call xchelu(dvxc(:,1),exc,fft_size,ORDER,rspts,vxc(:,1)) + ! + case (XC_LDA_C_XALPHA) + call xcxalp(3._DP/2._DP,exc,fft_size,n_spin,ORDER,rspts,zeta,vxc,dvxc) + ! + case (XC_LDA_C_KP) + call xc_rpa_kp(rspts,exc) + ! + end select + ! + ! + endif + ! + ! CORRELATION ONLY + ! + if (KIND==XC_CORRELATION) then + exc =exc-ex + vxc =vxc-vx + dvxc=dvxc-dvx + endif + ! + ! Allocation and Transfer + ! + select case (ORDER) + case(0) + E_xc=exc + case(1) + V_xc=vxc + case(2) + F_xc=dvxc(:,1) + end select + ! +end subroutine diff --git a/src/xc_functionals/XC_libxc_driver.F b/src/xc_functionals/XC_libxc_driver.F new file mode 100644 index 0000000000..cc51dbd1cc --- /dev/null +++ b/src/xc_functionals/XC_libxc_driver.F @@ -0,0 +1,103 @@ +! +! Copyright (C) 2000-2010 M. Gruning and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine XC_libxc_driver(en,Xk,KIND,FUNCTIONAL,ORDER) + ! + ! In case yambo is linked with libxc, this routine use the libxc + ! routine to evaluate the xc functional (28/02/2010: only LDA) + ! + use pars, ONLY:SP,DP,PI + use R_lattice, ONLY:bz_samp + use D_lattice, ONLY:DL_vol + use com, ONLY:msg,warning + use electrons, ONLY:levels,n_spin + use FFT_m, ONLY:fft_size + use libxc_int, ONLY:xc_fnctl_t,XC_y2libxc,xc_setup_fnctl + use xc_functionals,ONLY:E_xc,F_xc,V_xc,two_spin_density,& +& xc_string,XC_LDA_C_KP,XC_FAMILY_LDA + implicit none + type(levels) ::en + type(bz_samp)::Xk + integer ::KIND,FUNCTIONAL,ORDER + ! + ! Work Space + ! + integer :: i1,i2,ixc,LIBXC_FUNC + ! + ! XC routines arrays + ! + real(SP) :: rho_sp(fft_size,n_spin),rspts(fft_size) + real(DP) :: tmp_rho(n_spin,fft_size) + ! + ! note that multidimensional dvxc & vxc is allowed only in xcspol + ! + real(DP) :: tmp_dvxc(fft_size),tmp_vxc(n_spin,fft_size),tmp_exc(fft_size) + ! + ! Functional and Functional infos + ! + type(xc_fnctl_t) :: fnctl(2) + ! + !=============================== + ! Convert and Initialize + !=============================== + ! + tmp_exc=0._DP + tmp_vxc=0._DP + tmp_dvxc=0._DP + ! + LIBXC_FUNC = XC_y2libxc(FUNCTIONAL,KIND) + call xc_setup_fnctl(fnctl,LIBXC_FUNC,n_spin) + ! + !=============================== + ! Electronic/spin density + !=============================== + ! + call el_density(en,Xk,rho_sp,.false.) + ! + forall(i1=1:fft_size,i2=1:n_spin) tmp_rho(i2,i1)=rho_sp(i1,i2)*real(fft_size,DP)/real(DL_vol,DP) + ! + !=========================== + ! Evaluate the xc functional + !=========================== + ! + if (order==0) E_xc=0. + if (order==1) V_xc=0. + if (order==2) F_xc=0. + ! + ! + do ixc = 1,2 + if (fnctl(ixc)%id == 0) cycle + if (fnctl(ixc)%family==XC_FAMILY_LDA) then + select case(ORDER) + case(0) + call xc_f90_lda_exc(fnctl(ixc)%conf, fft_size, tmp_rho(1,1), tmp_exc(1)) + E_xc = E_xc + tmp_exc + case(1) + call xc_f90_lda_vxc(fnctl(ixc)%conf, fft_size, tmp_rho(1,1), tmp_vxc(1,1)) + forall(i1=1:fft_size,i2=1:n_spin) V_xc(i1,i2) = V_xc(i1,i2) + tmp_vxc(i2,i1) + case(2) + call xc_f90_lda_fxc(fnctl(ixc)%conf, fft_size, tmp_rho(1,1), tmp_dvxc(1)) + F_xc = F_xc + tmp_dvxc + end select + end if + enddo + ! +end subroutine diff --git a/src/xc_functionals/XC_potential_driver.F b/src/xc_functionals/XC_potential_driver.F new file mode 100644 index 0000000000..b6ea534128 --- /dev/null +++ b/src/xc_functionals/XC_potential_driver.F @@ -0,0 +1,76 @@ +! +! Copyright (C) 2000-2010 M. Gruning and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine XC_potential_driver(en,Xk,KIND,FUNCTIONAL,ORDER) + ! + ! Wrapper for calculation of the XC potential. If available libxc, + ! will be used instead of internal routines + ! + use pars, ONLY:SP,DP + use R_lattice, ONLY:bz_samp + use electrons, ONLY:levels + use com, ONLY:msg,warning + use xc_functionals,ONLY:XC_EXCHANGE_CORRELATION,XC_LDA_C_PW,xc_string,& + & XC_NOT_AVAILABLE,xc_OnlyInLibxc + implicit none + type(levels) ::en + type(bz_samp)::Xk + integer ::KIND,FUNCTIONAL,ORDER +! +! COM +! +#if defined _LIBXC + if (KIND>0.and.FUNCTIONAL>0) then +#else + if(xc_OnlyInLibxc(FUNCTIONAL)) then + call warning('XC functional not available. Default used instead') + call warning('By linking Yambo to libxc, this functional is available') + call msg('rsn','[xc] Requested functional: ', xc_string(KIND,FUNCTIONAL)) + KIND=XC_EXCHANGE_CORRELATION + FUNCTIONAL=XC_LDA_C_PW + call msg('rsn','[xc] Functional used: ', xc_string(KIND,FUNCTIONAL)) + elseif (KIND>0.and.FUNCTIONAL>0) then +#endif + call msg('rsn','[xc] Functional ',xc_string(KIND,FUNCTIONAL)) + else if (FUNCTIONAL==XC_NOT_AVAILABLE) then + KIND=XC_EXCHANGE_CORRELATION + FUNCTIONAL=XC_LDA_C_PW + call msg('rsn','[xc] Functional unknown. Used ',xc_string(KIND,FUNCTIONAL)) + endif + ! +#if defined _LIBXC + ! + ! Maybe also (spin)density should be calculated in the wrapper, since it does + ! not use anything of libxc. + ! + call msg('rsn','[xc] LIBXC used to calculate xc functional ') + ! + call XC_libxc_driver(& + ! +#else + call XC_lda_driver(& +#endif +& en,Xk,KIND,FUNCTIONAL,ORDER& +& ) + ! +end subroutine XC_potential_driver + + diff --git a/src/xc_functionals/el_current.F b/src/xc_functionals/el_current.F new file mode 100644 index 0000000000..acd307c73e --- /dev/null +++ b/src/xc_functionals/el_current.F @@ -0,0 +1,46 @@ +! +! Copyright (C) 2000-2010 D. Sangalli and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine el_current(en,Xk,ik,J,J_diam,l_directions) + ! + ! Electronic current: j = (q/m)* < [p - (q/c) *A] > + ! (q=-1,m=1) ==> = -1 < [p + (1/137.)*A] > + ! + use pars, ONLY:SP + use electrons, ONLY:levels,n_spin,n_sp_pol,Boltz_f + use R_lattice, ONLY:bz_samp + use D_lattice, ONLY:DL_sop + use FFT_m, ONLY:fft_size,fft_rot_r + use wave_func, ONLY:wf_state,wf,wf_x,wf_y,wf_z + use units, ONLY:SPEED_OF_LIGHT + ! + implicit none + type(bz_samp)::Xk + type(levels) ::en + ! + integer :: ik + logical :: l_directions(3) + real(SP) :: J(3,n_sp_pol,fft_size) + real(SP) :: J_diam(3,n_sp_pol,fft_size) + ! + ! + ! +end subroutine diff --git a/src/xc_functionals/spin_density.F b/src/xc_functionals/spin_density.F new file mode 100644 index 0000000000..620e967dac --- /dev/null +++ b/src/xc_functionals/spin_density.F @@ -0,0 +1,53 @@ +! +! Copyright (C) 2000-2010 D. De Fausti and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine spin_density(en,Xk) + ! + ! Electronic Spin density for electronic Spinors + ! + ! SD = Sum_I (PSI^I)* Sigma_Pauli* PSI^I I=occupied states + ! + ! PSI^I = spinor + ! Sigma_Pauli=the 3 Pauli matrices + ! + use pars, ONLY:SP + use R_lattice, ONLY:bz_samp + use electrons, ONLY:levels,n_spin,n_spinor + use D_lattice, ONLY:nsym + use FFT_m, ONLY:fft_size,fft_rot_r + use wave_func, ONLY:wf_state,wf + use xc_functionals, ONLY:two_spin_density + implicit none + type(levels) ::en + type(bz_samp)::Xk + ! + ! Work Space + ! + integer :: i1,i2,ifft_up,ifft_dn + real(SP):: cv(fft_size,3) + ! + two_spin_density=0. + cv=0. + ! + if (n_spin==1) return + ! + ! +end subroutine diff --git a/src/xc_functionals/xc_lda_driver.F b/src/xc_functionals/xc_lda_driver.F index fe65d7adcb..e69de29bb2 100644 --- a/src/xc_functionals/xc_lda_driver.F +++ b/src/xc_functionals/xc_lda_driver.F @@ -1,185 +0,0 @@ -! -! Copyright (C) 2000-2010 A. Marini and the YAMBO team -! http://www.yambo-code.org -! -! This file is distributed under the terms of the GNU -! General Public License. You can redistribute it and/or -! modify it under the terms of the GNU General Public -! License as published by the Free Software Foundation; -! either version 2, or (at your option) any later version. -! -! This program is distributed in the hope that it will -! be useful, but WITHOUT ANY WARRANTY; without even the -! implied warranty of MERCHANTABILITY or FITNESS FOR A -! PARTICULAR PURPOSE. See the GNU General Public License -! for more details. -! -! You should have received a copy of the GNU General Public -! License along with this program; if not, write to the Free -! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, -! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. -! -subroutine xc_lda_driver(en,Xk,KIND,FUNCTIONAL,ORDER) - ! - ! These sources have been slightly adpated from Abinit 4.6.5. - ! - ! xcspol: XC_LDA_C_PW - ! LDA or LSD, Teter Pade parametrization (4/93, published in S. Goedecker, - ! M. Teter, J. Huetter, Phys.Rev.B54, 1703 (1996)), - ! which reproduces Perdew-Wang (which reproduces Ceperley-Alder!). - ! - ! xcpzca: XC_LDA_C_PZ - ! LDA, Perdew-Zunger-Ceperley-Alder (no spin-polarization) - ! - ! xcwign: XC_LDA_C_WIGNER - ! LDA, Wigner functional (no spin-polarization) - ! - ! ! xchelu: XC_LDA_C_HL - ! LDA, Hedin-Lundqvist functional (no spin-polarization) - ! - ! xcxalp: XC_LDA_C_XALPHA - ! LDA, "X-alpha" functional (no spin-polarization) - ! - ! Note that, in order to get LDA exchange and correlation separated - ! - ! Exx = -efac*rsm1 - ! Vxx = -vfac*rsm1 - ! - ! in xcpzca xcwign xchelu xcxalp. While in xcspol xctetr there is not - ! such a type of distinction. - ! - ! Andrea (9/8/06) - ! - ! In case of spin polarized systems *REMEMBER* that: - ! - ! whatever(1) == whatever UP - ! whatever(2) == whatever DOWN - ! - ! Yambo have been considering the wrong order up to revision 335 (November 2008) - ! - use pars, ONLY:SP,DP,pi - use R_lattice, ONLY:bz_samp - use D_lattice, ONLY:DL_vol - use com, ONLY:msg,warning - use electrons, ONLY:levels,n_spin - use FFT_m, ONLY:fft_size - use xc_functionals,ONLY:XC_EXCHANGE,XC_EXCHANGE_CORRELATION,XC_CORRELATION,& -& E_xc,F_xc,V_xc,XC_LDA_X,& -& XC_LDA_C_PW,xc_string,XC_LDA_C_PZ,& -& XC_LDA_C_WIGNER,XC_LDA_C_HL,XC_LDA_C_XALPHA,& -& XC_LDA_C_KP,two_spin_density,NONE - implicit none - type(levels) ::en - type(bz_samp)::Xk - integer ::KIND,FUNCTIONAL,ORDER - ! - ! Work Space - ! - integer :: i1 - ! - ! XC routines arrays - ! - real(SP) :: rho_sp(fft_size) - real(DP) :: rspts(fft_size),zeta(fft_size),rho(fft_size) - ! - ! note that multidimensional dvxc & vxc is allowed only in xcspol - ! - real(DP) :: dvxc(fft_size,3),exc(fft_size),vxc(fft_size,n_spin) - real(DP) :: dvx(fft_size,3), ex(fft_size), vx(fft_size,n_spin) - ! - ! Init - ! - zeta=0._DP - exc=0._DP - vxc=0._DP - dvxc=0._DP - ex=0._DP - vx=0._DP - dvx=0._DP - ! - ! COM - ! - if (KIND>0.and.FUNCTIONAL>0) then - call msg('rsn','[xc] Functional ',xc_string(KIND,FUNCTIONAL)) - else if (KIND<0.or.FUNCTIONAL<0) then - KIND=XC_EXCHANGE_CORRELATION - FUNCTIONAL=XC_LDA_C_PW - call msg('rsn','[xc] Functional unknown. Used ',xc_string(KIND,FUNCTIONAL)) - endif - ! - ! Electronic density & Rs - ! - call el_density(en,Xk,rho_sp,.false.) - rho=rho_sp*real(fft_size,DP)/real(DL_vol,DP) - forall (i1=1:fft_size) rspts(i1)=(3._DP/4._DP/pi/rho(i1))**(1._DP/3._DP) - ! - ! - ! Hartree - ! - if (KIND==NONE) then - if (order==0) E_xc=0. - if (order==1) V_xc=0. - if (order==2) F_xc=0. - return - endif - ! - ! EXX = X_alpha when alpha=3/4. HF when alpha=1. - ! - if (KIND/=XC_EXCHANGE_CORRELATION.or.FUNCTIONAL==XC_LDA_X) then - call xcxalp(1._DP,ex,fft_size,n_spin,ORDER,rspts,zeta,vx,dvx) - if (KIND==XC_EXCHANGE.or.FUNCTIONAL==XC_LDA_X) then - exc=ex - vxc=vx - dvxc=dvx - endif - endif - ! - ! CORRELATION - ! - if (KIND/=XC_EXCHANGE) then - ! - select case (FUNCTIONAL) - ! - case (XC_LDA_C_PW) - call xcspol(dvxc,exc,fft_size,n_spin,ORDER,rspts,vxc,zeta) - ! - case (XC_LDA_C_PZ) - call xcpzca(dvxc(:,1),exc,fft_size,ORDER,rho,rspts,vxc(:,1)) - ! - case (XC_LDA_C_WIGNER) - call xcwign(dvxc(:,1),exc,fft_size,ORDER,rho,rspts,vxc(:,1)) - ! - case (XC_LDA_C_HL) - call xchelu(dvxc(:,1),exc,fft_size,ORDER,rspts,vxc(:,1)) - ! - case (XC_LDA_C_XALPHA) - call xcxalp(3._DP/2._DP,exc,fft_size,n_spin,ORDER,rspts,zeta,vxc,dvxc) - ! - case (XC_LDA_C_KP) - call xc_rpa_kp(rspts,exc) - ! - end select - ! - ! - endif - ! - ! CORRELATION ONLY - ! - if (KIND==XC_CORRELATION) then - exc =exc-ex - vxc =vxc-vx - dvxc=dvxc-dvx - endif - ! - ! Allocation and Transfer - ! - select case (ORDER) - case(0) - E_xc=exc - case(1) - V_xc=vxc - case(2) - F_xc=dvxc(:,1) - end select - ! -end subroutine diff --git a/ypp/.objects b/ypp/.objects index 6f02eb44cf..4599461984 100644 --- a/ypp/.objects +++ b/ypp/.objects @@ -4,4 +4,4 @@ -objs = mod_YPP.o ypp_i.o ypp_init.o ypp_init_load.o k_grids.o excitons_driver.o electrons_driver.o plot_check_and_launch.o plot_gnuplot.o plot_xcrysden.o excitons_sort_and_report.o $(SPIN_objects) $(RAS_objects) +objs = mod_YPP.o ypp_i.o ypp_init.o ypp_init_load.o k_grids.o excitons_driver.o electrons_driver.o plot_check_and_launch.o plot_gnuplot.o plot_xcrysden.o plot_cube.o excitons_sort_and_report.o $(SPIN_objects) $(ELPH_objects) $(RAS_objects) $(RT_objects) $(magnetic_objects) $(boltz_objects) diff --git a/ypp/electrons_driver.F b/ypp/electrons_driver.F index 7aa1231f85..e137a4382d 100644 --- a/ypp/electrons_driver.F +++ b/ypp/electrons_driver.F @@ -31,7 +31,7 @@ subroutine electrons_driver(Xk,Xen) use electrons, ONLY:levels,n_spin,n_spinor,spin,n_sp_pol use QP_m, ONLY:QP_table,QP_n_states,QP_state use YPP, ONLY:l_density,l_mag,v2plot,output_fname,plot_dim,use_xcrysden,& -& use_gnuplot,nr,l_sp_wf,deg_energy,mag_dir,l_norm_to_one,& +& use_gnuplot,use_cube,nr,l_sp_wf,deg_energy,mag_dir,l_norm_to_one,& & plot_title,l_dos,dos_broadening,DOS_bands,l_angular_momentum,& & l_position,l_current use com, ONLY:msg,of_open_close,warning @@ -94,7 +94,7 @@ subroutine electrons_driver(Xk,Xen) titles(2)='DOS [up]' titles(3)='DOS [dn]' titles(4)='DOS [up+dn]' - call msg('o dos','#',titles(:2),INDENT=0,USE_TABS=.true.) + call msg('o dos','#',titles(:4),INDENT=0,USE_TABS=.true.) endif call msg('o dos',"#") ! @@ -138,11 +138,16 @@ subroutine electrons_driver(Xk,Xen) v2plot=el_den ! ch_ws='density_'//trim(intc(plot_dim)) + if (use_cube) output_fname=trim(ch_ws)//'d.cube' if (use_xcrysden) output_fname=trim(ch_ws)//'d.xsf' if (use_gnuplot) output_fname=trim(ch_ws)//'d' ! - call of_open_close(trim(output_fname),'ot') - call msg('o den',"#") + if (use_cube) then + call of_open_close(trim(output_fname),'o') + else + call of_open_close(trim(output_fname),'ot') + call msg('o den',"#") + endif ! plot_title='density' call plot_check_and_launch(.false.) @@ -191,11 +196,16 @@ subroutine electrons_driver(Xk,Xen) ch_ws='sp_wf_k'//trim(intc(ik))//'_b'//trim(intc(ib))//'_'//trim(intc(plot_dim)) ! ! + if (use_cube) output_fname=trim(ch_ws)//'d.cube' if (use_xcrysden) output_fname=trim(ch_ws)//'d.xsf' if (use_gnuplot) output_fname=trim(ch_ws)//'d' ! - call of_open_close(trim(output_fname),'ot') - call msg('o wf',"#") + if (use_cube) then + call of_open_close(trim(output_fname),'o') + else + call of_open_close(trim(output_fname),'ot') + call msg('o wf',"#") + endif ! call plot_check_and_launch(.false.) ! diff --git a/ypp/excitons_driver.F b/ypp/excitons_driver.F index cd0684bc69..55e3e82c99 100644 --- a/ypp/excitons_driver.F +++ b/ypp/excitons_driver.F @@ -31,7 +31,7 @@ subroutine excitons_driver(k,Xk,en,Xen,q) use YPP, ONLY:l_sort,l_exc_wf,l_spin,ncell,lambda,r_hole,& & nr,v2plot,nr_tot,state_ctl,BS_H_dim,n_lambda,& & excitons_sort_and_report,l_amplitude,l_free_hole,& -& output_fname,use_xcrysden,use_gnuplot,plot_dim,l_norm_to_one +& output_fname,use_xcrysden,use_gnuplot,use_cube,plot_dim,l_norm_to_one use IO_m, ONLY:io_control,OP_RD_CL,DUMP,NONE,REP,VERIFY use BS, ONLY:BS_K_dim,BS_bands,BS_q,& & BS_eh_table,BS_mat,BS_K_coupling,BS_cpl_mode,& @@ -391,14 +391,19 @@ subroutine excitons_driver(k,Xk,en,Xen,q) ! ! PLOT ! + if (use_cube) output_fname='exc_'//trim(intc(plot_dim))//'d_'//trim(intc(lambda))//'.cube' if (use_xcrysden) output_fname='exc_'//trim(intc(plot_dim))//'d_'//trim(intc(lambda))//'.xsf' if (use_gnuplot) output_fname='exc_'//trim(intc(plot_dim))//'d_'//trim(intc(lambda)) ! - call of_open_close(trim(output_fname),'ot') - do i1=1,BSS_n_descs - call msg('o exc',"#",trim(BSS_description(i1)),INDENT=0) - enddo - call msg('o exc',"#") + if (use_cube) then + call of_open_close(trim(output_fname),'o') + else + call of_open_close(trim(output_fname),'ot') + do i1=1,BSS_n_descs + call msg('o exc',"#",trim(BSS_description(i1)),INDENT=0) + enddo + call msg('o exc',"#") + endif ! call plot_check_and_launch(.false.) ! diff --git a/ypp/excitons_sort_and_report.F b/ypp/excitons_sort_and_report.F index a24b33d026..563769650e 100644 --- a/ypp/excitons_sort_and_report.F +++ b/ypp/excitons_sort_and_report.F @@ -26,12 +26,14 @@ subroutine excitons_sort_and_report(BS_H_dim,Xk,E,BS_R,BS_E,BS_E_degs,& use units, ONLY:HARTREE use electrons, ONLY:spin_occ,levels,spin use stderr, ONLY:intc - use BS, ONLY:BSS_description,BSS_n_descs,BS_eh_table + use BS, ONLY:BSS_description,BSS_n_descs,BS_eh_table,BS_bands use com, ONLY:msg,of_open_close use YPP, ONLY:deg_energy,lambda use R_lattice, ONLY:d3k_factor,q0_def_norm,bz_samp use D_lattice, ONLY:n_atomic_species,n_atoms_species,atom_pos,Z_species use vec_operate, ONLY:sort + use QP_CTL_m, ONLY:QP_apply + ! implicit none ! integer ::BS_H_dim @@ -104,9 +106,13 @@ subroutine excitons_sort_and_report(BS_H_dim,Xk,E,BS_R,BS_E,BS_E_degs,& allocate(S_indx(BS_H_dim),amp_trans(BS_H_dim,2)) ! ! Sort the weights - ! + ! call sort(arrin=A_weight,indx=S_indx) ! + ! Apply quasi-particle correction if presents + ! + call QP_apply(BS_bands,E,Xk,3,msg_fmt='s') + ! ! report on file the weights and the amplitude ! of the excitonic state... ! @@ -175,7 +181,7 @@ subroutine excitons_sort_and_report(BS_H_dim,Xk,E,BS_R,BS_E,BS_E_degs,& if (A_weight(neh)<0.05) cycle ! call msg('o weight','',(/real(iv,SP),real(ic,SP),real(ikibz,SP),real(is,SP),& -& A_weight(neh),amp_trans(amp_n_trans,1)/),& +& A_weight(neh),amp_trans(amp_n_trans,1)*HARTREE/),& & INDENT=-2,USE_TABS=.FALSE.) enddo ! diff --git a/ypp/k_grids.F b/ypp/k_grids.F index e533363d4f..0b48cf27df 100644 --- a/ypp/k_grids.F +++ b/ypp/k_grids.F @@ -202,10 +202,10 @@ subroutine k_grids(en,k,Xk,q) ! if (l_long_grid) then ! - ! Generate shifted k-points set for longitudinal gauge calculation + ! Generate shifted k-points set for dipole calculation ! - call section('*',"== Longitudinal gauge grids generator ==") - !=========================================================== + call section('*',"== Dipole grids generator ==") + !================================================ ! call parser('KShift1',qlong(1,:)) call parser('KShift2',qlong(2,:)) diff --git a/ypp/mod_YPP.F b/ypp/mod_YPP.F index 974226f3ad..63c6cb5573 100644 --- a/ypp/mod_YPP.F +++ b/ypp/mod_YPP.F @@ -21,7 +21,7 @@ ! module YPP ! - use pars, ONLY:SP + use pars, ONLY:SP,DP use pars, ONLY:lchlen use pars, ONLY:schlen implicit none @@ -60,6 +60,7 @@ module YPP logical ::l_norm_to_one logical ::use_gnuplot logical ::use_xcrysden + logical ::use_cube character(1) ::p_format character(1) ::mag_dir character(3) ::p_dir diff --git a/ypp/plot_check_and_launch.F b/ypp/plot_check_and_launch.F index 5cc7c9eb0a..17b0870e4f 100644 --- a/ypp/plot_check_and_launch.F +++ b/ypp/plot_check_and_launch.F @@ -22,7 +22,7 @@ subroutine plot_check_and_launch(check_only) ! use pars, ONLY:SP - use YPP, ONLY:p_format,p_dir,plot_dim,use_gnuplot,use_xcrysden + use YPP, ONLY:p_format,p_dir,plot_dim,use_gnuplot,use_xcrysden,use_cube use com, ONLY:warning use D_lattice, ONLY:a use vec_operate, ONLY:v_is_zero @@ -44,10 +44,11 @@ subroutine plot_check_and_launch(check_only) ! use_gnuplot=p_format=='g' use_xcrysden=p_format=='x' + use_cube=p_format=='c' switch_to_xcrysden=.false. switch_to_gnuplot=.false. ! - if (.not.use_gnuplot.and..not.use_xcrysden) use_gnuplot=.true. + if (.not.use_gnuplot.and..not.use_xcrysden.and..not.use_cube) use_gnuplot=.true. ! use_1d(1) =p_dir(:1)=='1' use_1d(2) =p_dir(:1)=='2' @@ -91,6 +92,8 @@ subroutine plot_check_and_launch(check_only) endif ! if (plot_dim==1.and.use_xcrysden) switch_to_gnuplot=.true. + if (plot_dim==1.and.use_cube) switch_to_gnuplot=.true. + if (plot_dim==2.and.use_cube) switch_to_gnuplot=.true. if (plot_dim==3.and.use_gnuplot) switch_to_xcrysden=.true. ! if (switch_to_xcrysden.and.check_only) call warning('Switching to XCrySDen') @@ -102,6 +105,7 @@ subroutine plot_check_and_launch(check_only) if (check_only) return ! if (use_gnuplot) call plot_gnuplot(use_1d,use_2d,.FALSE.) + if (use_cube) call plot_cube() if (use_xcrysden) then if (plot_dim==2) call plot_gnuplot(use_1d,use_2d,.TRUE.) call plot_xcrysden(use_2d) diff --git a/ypp/plot_cube.F b/ypp/plot_cube.F new file mode 100644 index 0000000000..12c1370cd1 --- /dev/null +++ b/ypp/plot_cube.F @@ -0,0 +1,104 @@ +! +! Copyright (C) 2000-2010 D. Varsano, A. Marini and the YAMBO team +! http://www.yambo-code.org +! +! This file is distributed under the terms of the GNU +! General Public License. You can redistribute it and/or +! modify it under the terms of the GNU General Public +! License as published by the Free Software Foundation; +! either version 2, or (at your option) any later version. +! +! This program is distributed in the hope that it will +! be useful, but WITHOUT ANY WARRANTY; without even the +! implied warranty of MERCHANTABILITY or FITNESS FOR A +! PARTICULAR PURPOSE. See the GNU General Public License +! for more details. +! +! You should have received a copy of the GNU General Public +! License along with this program; if not, write to the Free +! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, +! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt. +! +subroutine plot_cube() + ! + use pars, ONLY:SP,schlen + use FFT_m, ONLY:fft_dim + use LOGO, ONLY:code_version + use com, ONLY:msg + use YPP, ONLY:nr,v2plot,ncell,r_hole,plot_dim,l_free_hole,& +& l_norm_to_one,WF_multiplier,l_exc_wf,plot_title,output_string + use D_lattice, ONLY:n_atomic_species,n_atoms_species,a,atom_pos,Z_species,n_atoms + use timing, ONLY:live_timing + ! + implicit none + ! + ! Work Space... + ! + integer :: i1,i2,i3,ir,is,ia + real(SP) :: rv(3),max_ + character(schlen) :: ch + ! + call msg(output_string,'CUBE FILE') + call msg(output_string,'Generated with YPP',code_version) + if (.not.l_free_hole.and.l_exc_wf) then + write(ch,'(i5,3f10.5)') n_atoms*ncell(1)*ncell(2)*ncell(3)+1,0._SP,0._SP,0._SP + else + write(ch,'(i5,3f10.5)') n_atoms*ncell(1)*ncell(2)*ncell(3),0._SP,0._SP,0._SP + endif + call msg(output_string,'',ch) + do i1=1,3 + write(ch,'(i5,3f10.5)') nr(i1),a(i1,:)/fft_dim(i1) + call msg(output_string,'',ch) + enddo + ! + ! write the atoms position + ! + if (.not.l_free_hole.and.l_exc_wf) then + write(ch,'(2i5,3f10.5)') -1,-1,r_hole + call msg(output_string,'',ch,INDENT=0,USE_TABS=.FALSE.) + endif + + do is=1,n_atomic_species + do ia=1,n_atoms_species(is) + do i1=0,ncell(1)-1 + do i2=0,ncell(2)-1 + do i3=0,ncell(3)-1 + rv(1)=atom_pos(1,ia,is)+i1*a(1,1)+i2*a(2,1)+i3*a(3,1) + rv(2)=atom_pos(2,ia,is)+i1*a(1,2)+i2*a(2,2)+i3*a(3,2) + rv(3)=atom_pos(3,ia,is)+i1*a(1,3)+i2*a(2,3)+i3*a(3,3) + write(ch,'(2i5,3f10.5)') Z_species(is), Z_species(is),rv(:) + call msg(output_string,'',ch,INDENT=0,USE_TABS=.FALSE.) + enddo + enddo + enddo + enddo + enddo + ! + write(ch,'(2i5)') 1,1 + call msg(output_string,'',ch,INDENT=0,USE_TABS=.FALSE.) + ! + ir = 0 + max_=maxval(v2plot) + if (.not.l_norm_to_one) max_=1. + ! + if (len_trim(plot_title)==0) then + call live_timing('3D Plot',nr(1),SERIAL=.true.) + else + call live_timing('3D Plot of '//trim(plot_title),nr(1),SERIAL=.true.) + endif + ! + do i1 = 0, nr(1)-1 + do i2 = 0, nr(2)-1 + do i3 = 0, nr(3)-1 + ir = 1 + i1 + i2*nr(1) + i3*nr(1)*nr(2) + call msg(output_string,'',v2plot(ir)/max_*WF_multiplier) + enddo + enddo + ! + call live_timing(steps=1) + ! + enddo + ! + call live_timing() + ! +end subroutine diff --git a/ypp/ypp_init.F b/ypp/ypp_init.F index 5aeee9fcd7..ac0a358c7c 100644 --- a/ypp/ypp_init.F +++ b/ypp/ypp_init.F @@ -29,7 +29,7 @@ subroutine ypp_init(instr,lnstr,FINALIZE) use it_m, ONLY:initdefs,initmode,ofiles_append,& & initinfio,infile,infile_dump,initactivate,& & nrnlvls,rnlvls,runlevel_is_on,& -& infile_verbosity,V_general,V_qp +& infile_verbosity,V_general,V_qp,V_all use drivers, ONLY:infile_editing use com, ONLY:file_exists use stderr, ONLY:string_split @@ -133,6 +133,7 @@ subroutine ypp_init(instr,lnstr,FINALIZE) if (i1<50) then if( trim(rstr_piece(i1)) == 'infver' .and. trim(rstr_piece(i1+1))=='gen' ) infile_verbosity=V_general if( trim(rstr_piece(i1)) == 'infver' .and. trim(rstr_piece(i1+1))=='qp' ) infile_verbosity=V_qp + if( trim(rstr_piece(i1)) == 'infver' .and. trim(rstr_piece(i1+1))=='all' ) infile_verbosity=V_all endif ! do i2=1,nrnlvls @@ -230,6 +231,7 @@ subroutine ypp_init(instr,lnstr,FINALIZE) endif ! if (l_excitons) then + if (l_amplitude) call QP_ctl_switch('G') call initactivate(1,"States") if (l_exc_wf.or.l_amplitude) call initactivate(1,"Degen_Step") if (l_exc_wf.and..not.l_free_hole) call initactivate(1,"Cells Hole Dimension") diff --git a/ypp/ypp_init_load.F b/ypp/ypp_init_load.F index 87fe086d13..92592c651d 100644 --- a/ypp/ypp_init_load.F +++ b/ypp/ypp_init_load.F @@ -44,13 +44,13 @@ subroutine ypp_init_load(defs) call it('r',defs,'plot', '[R] Plot') call it('r',defs,'density', '[R] Density') call it('r',defs,'wavefunction', '[R] Wavefunction') + call QP_ctl_load(defs,3) ! ! ! DOS ! call it(defs,'DOS_broad', 'Broadening of the DOS',dos_broadening,E_unit) call it(defs,'DOS_bands', 'DOS bands',dos_bands) - call QP_ctl_load(defs,3) ! ! BZ grids ! @@ -81,9 +81,9 @@ subroutine ypp_init_load(defs) ! ! p_dir: plot cut in the a1,a2,a3 basis ! - ! p_format: (g)nuplot/(x)crysden + ! p_format: (c)ube/(g)nuplot/(x)crysden ! - call it(defs,'Format', 'Output format [(g)nuplot/(x)crysden]',p_format) + call it(defs,'Format', 'Output format [(c)ube/(g)nuplot/(x)crysden]',p_format) call it(defs,'Direction', '[rlu] [1/2/3] for 1d or [12/13/23] for 2d [123] for 3D',p_dir) ! !