From ab3ad02c4a5053bf959adc569c05b5699668c348 Mon Sep 17 00:00:00 2001 From: giannozz Date: Mon, 11 Dec 2017 10:41:41 +0000 Subject: [PATCH] Added new version of pw2gw with documentation (Olivia Pulci, Davide Grassano, Adriano Mosca Conte); version number in documentation updated. Maybe ready? git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@14038 c92efa57-630b-4861-b058-cf58834340f0 --- CPV/Doc/user_guide.tex | 2 +- Doc/developer_man.tex | 2 +- Doc/user_guide.tex | 2 +- NEB/Doc/user_guide.tex | 2 +- PHonon/Doc/developer_man.tex | 2 +- PHonon/Doc/user_guide.tex | 2 +- PP/Doc/INPUT_pw2gw.html | 265 +++++++++++++++++++++++++++ PP/Doc/Makefile | 1 + PP/Doc/user_guide.tex | 10 +- PP/src/pw2gw.f90 | 346 ++++++++++++++++++++--------------- dev-tools/release.sh | 2 +- 11 files changed, 479 insertions(+), 157 deletions(-) create mode 100644 PP/Doc/INPUT_pw2gw.html diff --git a/CPV/Doc/user_guide.tex b/CPV/Doc/user_guide.tex index 1286fca70d..e4e427adf7 100644 --- a/CPV/Doc/user_guide.tex +++ b/CPV/Doc/user_guide.tex @@ -1,5 +1,5 @@ \documentclass[12pt,a4paper]{article} -\def\version{6.2} +\def\version{6.2.1} \def\qe{{\sc Quantum ESPRESSO}} \usepackage{html} diff --git a/Doc/developer_man.tex b/Doc/developer_man.tex index 24841f1871..58bc40e91c 100644 --- a/Doc/developer_man.tex +++ b/Doc/developer_man.tex @@ -1,5 +1,5 @@ \documentclass[12pt,a4paper]{article} -\def\version{6.2} +\def\version{6.2.1} \def\qe{{\sc Quantum ESPRESSO}} \def\qeforge{\texttt{qe-forge.org}} \textwidth = 17cm diff --git a/Doc/user_guide.tex b/Doc/user_guide.tex index ce853aa09f..9b296770ac 100644 --- a/Doc/user_guide.tex +++ b/Doc/user_guide.tex @@ -1,5 +1,5 @@ \documentclass[12pt,a4paper]{article} -\def\version{6.2} +\def\version{6.2.1} \def\qe{{\sc Quantum ESPRESSO}} \usepackage{html} diff --git a/NEB/Doc/user_guide.tex b/NEB/Doc/user_guide.tex index e1d648a96c..ab79c65de9 100644 --- a/NEB/Doc/user_guide.tex +++ b/NEB/Doc/user_guide.tex @@ -1,5 +1,5 @@ \documentclass[12pt,a4paper]{article} -\def\version{6.2} +\def\version{6.2.1} \def\qe{{\sc Quantum ESPRESSO}} \def\NEB{\texttt{PWneb}} % to be decided \usepackage{html} diff --git a/PHonon/Doc/developer_man.tex b/PHonon/Doc/developer_man.tex index 1041ac6f60..7e79fc2f13 100644 --- a/PHonon/Doc/developer_man.tex +++ b/PHonon/Doc/developer_man.tex @@ -1,5 +1,5 @@ \documentclass[12pt,a4paper]{article} -\def\version{6.2} +\def\version{6.2.1} \def\qe{{\sc Quantum ESPRESSO}} \def\qeforge{\texttt{qe-forge.org}} \textwidth = 17cm diff --git a/PHonon/Doc/user_guide.tex b/PHonon/Doc/user_guide.tex index e0dc5a9130..cd4fbb8b9d 100644 --- a/PHonon/Doc/user_guide.tex +++ b/PHonon/Doc/user_guide.tex @@ -1,5 +1,5 @@ \documentclass[12pt,a4paper]{article} -\def\version{6.2} +\def\version{6.2.1} \def\qe{{\sc Quantum ESPRESSO}} \usepackage{html} diff --git a/PP/Doc/INPUT_pw2gw.html b/PP/Doc/INPUT_pw2gw.html new file mode 100644 index 0000000000..d7d2c268bd --- /dev/null +++ b/PP/Doc/INPUT_pw2gw.html @@ -0,0 +1,265 @@ + + + + + +pw2gw.x: input description + + + + +
+

Input File Description

+

Program: + pw2gw.x / PWscf / Quantum Espresso (version: 6.x) +

+
+
+

TABLE OF CONTENTS

+
+ + +

INTRODUCTION

+

&INPUTPP

+
+prefix | outdir | what | qplda | vxcdiag | vkb | Emin | Emax | DeltaE +
+ +
+
+
+

INTRODUCTION

+
+Purpose of pw2gw.x: Optical properties in single-particle approach (Fermi Golden Rule).
+Interface with GW and excitonic codes.
+
+The code computes and writes ("matrixelements" file) the optical matrix elemenents in the 
+dipole approximation.
+
+The code computes the imaginary part of the dielectric tensor xx, yy and zz ("epsX.dat", "epsY.dat",
+ "epsZ.dat") and the average ("epsTOT.dat")
+
+
+
+ + +

Namelist: &INPUTPP +

+ + + + + + +
prefixCHARACTER
+prefix of files saved by program pw.x
+
+ + + + + + + + + + + +
outdirCHARACTER
Default: +value of the ESPRESSO_TMPDIR environment variable if set; +current directory ('./') otherwise + 
+directory containing the input data, i.e. the same as in pw.x
+
+ + + + + + + + + + + + + +
whatCHARACTER
Default: +'gw' + 
+'gw'    = Calculate dipole optical matrix elements (use for norm-conserving pseudopotentials)
+          and imaginary part of the dielectric function.
+          If qplda = .true., write the interface file "QPLDA" to GW and BSE codes (chisig, dpforexc).
+'gmaps' = write g-maps for each processor in a file "fort.'100 + processor number'"
+
+ + + + + + + + + + + + + + + + +
qpldaLOGICAL
Default: +.false. + 
+  .true.  = Write the interface file "QPLDA" to GW and BSE codes (chisig, dpforexc).
+  .false. = Skip this step.
+
+ + + + + + + + + + + + + + + + +
vxcdiagLOGICAL
Default: +.false. + 
+  .true.  = Calculates the expectation value of the exchange and correlation potential
+            on all the Kohn-Sham states and write it into the "vxcdiag.dat" file.
+  .false. = Skip this step.
+
+ + + + + + + + + + + + + + + + +
vkbLOGICAL
Default: +.false. + 
+  .true.  = Use Kleyman-Bylander projectors to write additional informatio into fort.15 file
+            (Still in development)
+  .false. = Skip this step.
+
+ + + + + + + + + + + + + + + + +
EminREAL
Default: +0.0 + 
+  Starting photon energy for which the dielectric function is calculated (in eV)
+
+ + + + + + + + + + + + + + + + +
EmaxREAL
Default: +30.0 + 
+  Highest photon energy for which the dielectric function is calculated (in eV)
+
+ + + + + + + + + + + + + + + + +
DeltaEREAL
Default: +0.05 + 
+  Energy step with which the dielectric function is calculated (in eV)
+
+ + + + diff --git a/PP/Doc/Makefile b/PP/Doc/Makefile index 252e66ed8c..747d851f6c 100644 --- a/PP/Doc/Makefile +++ b/PP/Doc/Makefile @@ -65,6 +65,7 @@ link_on_main_doc: ln -fs ../PP/Doc/INPUT_pw_export.html . ; \ ln -fs ../PP/Doc/INPUT_pw_export.xml . ; \ ln -fs ../PP/Doc/INPUT_pw_export.txt . ; \ + ln -fs ../PP/Doc/INPUT_pw2gw.html . ; \ ln -fs ../PP/Doc/INPUT_pw2bgw.html . ; \ ln -fs ../PP/Doc/INPUT_pw2bgw.xml . ; \ ln -fs ../PP/Doc/INPUT_pw2bgw.txt . ; \ diff --git a/PP/Doc/user_guide.tex b/PP/Doc/user_guide.tex index d8f128ec52..7133246639 100644 --- a/PP/Doc/user_guide.tex +++ b/PP/Doc/user_guide.tex @@ -1,5 +1,5 @@ \documentclass[12pt,a4paper]{article} -\def\version{6.2} +\def\version{6.2.1} \def\qe{{\sc Quantum ESPRESSO}} \usepackage{html} @@ -74,6 +74,8 @@ \section{People and terms of use} Stefano de Gironcoli, Andrea Dal Corso (SISSA), Paolo Giannozzi (Univ. Udine), and many others. We mention in particular: \begin{itemize} + \item Olivia Pulci, Adriano Mosca Conte, Davide Grassano (RomaII) + for the \texttt{pw2gw} utility; \item Andrea Benassi (SISSA) for the \texttt{epsilon} utility, Tae Yun Kim and Cheol-Hwan Park (Seoul National University) for fixes to it; @@ -313,8 +315,10 @@ \subsection{Interfaces to/from other code} versa. See files \texttt{Doc/INPUT\_pw2bgw.*} and \texttt{Doc/INPUT\_bgw2pw.*} for input data documentation. -Undocumented code \texttt{pw2gw.x} converts data files from \pwx\ to -a format suitable for usage by another GW code. +Code \texttt{pw2gw.x} converts data files from \pwx\ to a format suitable +for usage by another GW code, computes optical properties in single-particle +approach (Fermi Golden Rule). See file \texttt{Doc/INPUT\_pw2gw.html} +for input data documentation. Code \texttt{pw\_export.x}, not compiled by default, is an interface to other codes, documented in \texttt{Doc/INPUT\_pw\_export.*} diff --git a/PP/src/pw2gw.f90 b/PP/src/pw2gw.f90 index 87abad6783..70c7526c42 100644 --- a/PP/src/pw2gw.f90 +++ b/PP/src/pw2gw.f90 @@ -10,6 +10,17 @@ ! Please note just symmorphic symm. op. have to be used ! Use input option of pw.x: force_symmorphic=.TRUE. +! Update 20 November 2017 (Davide Grassano, Adriano Mosca Conte) +! Added spin-orbit case for nspin == 4 +! Added flag to input file: +! Emin: Starting energy of eps spectra +! Emax: Max energy of eps spectra +! DeltaE: Energy steps for eps sectra +! qplda: disable/enable priting of qplda file (default = false) +! vkb: disable/enable priting of fort.15 file (default = false) +! vxcdiag: disable/enable priting of vxcdiag.dat file (default = false) +! Using proper units(tpiba) for k-vectors in k+G sum before calculating matrixelements + !----------------------------------------------------------------------- PROGRAM pw2gw !----------------------------------------------------------------------- @@ -25,6 +36,7 @@ PROGRAM pw2gw USE mp_pools, ONLY : kunit USE environment,ONLY : environment_start, environment_end USE us, ONLY : spline_ps + USE kinds, ONLY : DP ! IMPLICIT NONE ! @@ -33,11 +45,12 @@ PROGRAM pw2gw ! INTEGER :: ios INTEGER :: kunittmp - LOGICAL :: use_gmaps + LOGICAL :: use_gmaps, qplda, vkb, vxcdiag CHARACTER(len=20) :: what - CHARACTER(len=30) :: when + REAL(kind=DP) :: Emin, Emax, DeltaE - NAMELIST / inputpp / prefix, outdir, what, use_gmaps + NAMELIST / inputpp / prefix, outdir, what, use_gmaps, Emin, Emax, DeltaE, & + qplda, vkb, vxcdiag ! ! initialise environment ! @@ -52,7 +65,13 @@ PROGRAM pw2gw CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir ) IF ( trim( outdir ) == ' ' ) outdir = './' what = 'gw' + qplda = .false. + vkb = .false. + vxcdiag = .false. use_gmaps = .false. + Emin = 0.0 + Emax = 30.0 + DeltaE = 0.05 ios = 0 IF ( ionode ) THEN @@ -71,6 +90,13 @@ PROGRAM pw2gw CALL mp_bcast(tmp_dir, ionode_id, world_comm ) CALL mp_bcast( what, ionode_id, world_comm ) CALL mp_bcast( use_gmaps, ionode_id, world_comm ) + CALL mp_bcast( qplda, ionode_id, world_comm ) + CALL mp_bcast( vkb, ionode_id, world_comm ) + CALL mp_bcast( vxcdiag, ionode_id, world_comm ) + CALL mp_bcast( Emin, ionode_id, world_comm ) + CALL mp_bcast( Emax, ionode_id, world_comm ) + CALL mp_bcast( DeltaE, ionode_id, world_comm ) + ! spline_ps = .false. @@ -85,9 +111,10 @@ PROGRAM pw2gw kunittmp = 1 #endif ! + !WRITE(*,*) what, qplda, vxcdiag, Emin, Emax, DeltaE IF( trim( what ) == 'gw' ) THEN - CALL compute_gw ( use_gmaps ) - ELSE + CALL compute_gw ( Emin, Emax, DeltaE, use_gmaps, qplda, vkb, vxcdiag ) + ELSEIF( trim( what ) == 'gmaps' ) THEN CALL write_gmaps ( kunittmp ) ENDIF ! @@ -98,7 +125,7 @@ PROGRAM pw2gw END PROGRAM pw2gw -SUBROUTINE compute_gw( use_gmaps ) +SUBROUTINE compute_gw( omegamin, omegamax, d_omega, use_gmaps, qplda, vkb, vxcdiag_f ) ! This routine creates the QPLDA and the matrixelements ! tform = .false. UNFORMATTED QPLDA @@ -136,7 +163,8 @@ SUBROUTINE compute_gw( use_gmaps ) IMPLICIT NONE - LOGICAL, INTENT(in) :: use_gmaps + LOGICAL, INTENT(in) :: use_gmaps, qplda, vkb, vxcdiag_f + REAL(kind=DP), INTENT( in) :: omegamax, omegamin, d_omega INTEGER :: ii(16), ngw, nkpt, ig, ik, ir, n, i,j,k, io = 98, iband1, iband2 INTEGER :: npw, omax, o, iproc @@ -148,7 +176,7 @@ SUBROUTINE compute_gw( use_gmaps ) REAL(kind=sgl), ALLOCATABLE :: xk_s(:,:), eig_s(:,:), focc_s(:,:) REAL(kind=DP):: g2max, a1(3), a2(3), a3(3),norm, xkgk(3), rrhotwx(3), delta REAL(kind=DP):: alpha, egap, halfalpha, Df, const, dummy - REAL(kind=DP), PARAMETER :: omegamax = 30.0 + !REAL(kind=DP), PARAMETER :: omegamax = 30.0 REAL(kind=DP), ALLOCATABLE:: gsort(:), eig(:,:), focc(:,:), kpg(:,:), omegatt(:), omeg(:) REAL(kind=DP), ALLOCATABLE:: pp1(:,:), pp2(:,:), pp3(:,:) REAL(kind=DP), ALLOCATABLE:: epsx(:,:), epsy(:,:), epsz(:,:) @@ -163,6 +191,7 @@ SUBROUTINE compute_gw( use_gmaps ) INTEGER :: igwx, igwxx, comm, ierr, ig_max, igwx_r INTEGER :: igwx_p(nproc) INTEGER, ALLOCATABLE :: igk_l2g(:) + INTEGER :: npol ! REAL(kind=DP), ALLOCATABLE :: vkb0(:), djl(:), vec_tab(:), vec_tab_d2y(:) INTEGER :: nb, nt, size_tab, size_tab_d2y, ipw, l @@ -173,17 +202,18 @@ SUBROUTINE compute_gw( use_gmaps ) INTEGER :: istatus( MPI_STATUS_SIZE ) #endif ! - IF( nspin > 1 ) CALL errore('pw2gw','Spin polarization not implemented',1) - IF( npool > 1 ) CALL errore('pw2gw','parallel run with pools not allowed yet',1) + IF( nspin == 4 ) WRITE(*,*) "nspin = 4" + IF( nspin == 2 ) CALL errore('pw2gw','Spin polarization not implemented',1) + IF( npool > 1 ) CALL errore('pw2gw','parallel run with pools not allowed yet',1) ! ! IF( mpime == 0 ) THEN IF (t_form) THEN WRITE (6,'(//" writing LDA info on unit 98 FORMATTED")') - OPEN (io, FILE='QPLDA',STATUS='unknown',FORM='FORMATTED') + IF( qplda) OPEN (io, FILE='QPLDA',STATUS='unknown',FORM='FORMATTED') ELSE WRITE (6,'(//" writing LDA info on unit io UNFORMATTED")') - OPEN (io, FILE='QPLDA',STATUS='unknown',FORM='UNFORMATTED') + IF( qplda) OPEN (io, FILE='QPLDA',STATUS='unknown',FORM='UNFORMATTED') ENDIF WRITE (6,'(//" writing matrixelements on unit 98 FORMATTED")') OPEN (90, FILE='matrixelements',STATUS='unknown',FORM='FORMATTED') @@ -195,10 +225,10 @@ SUBROUTINE compute_gw( use_gmaps ) titleo(2)='test version' IF( mpime == 0 ) THEN IF (t_form) THEN - WRITE (io,'(A80/A80)') titleo(1), titleo(2) + IF( qplda) WRITE (io,'(A80/A80)') titleo(1), titleo(2) ELSE - WRITE (io) titleo(1) - WRITE (io) titleo(2) + IF( qplda) WRITE (io) titleo(1) + IF( qplda) WRITE (io) titleo(2) ENDIF ! WRITE(6,*) 'qplda title' @@ -215,10 +245,10 @@ SUBROUTINE compute_gw( use_gmaps ) ii(:) = 0 IF (t_form) THEN ii(1)=0 - IF( mpime == 0 ) WRITE (io,'(16I5)') ii + IF( mpime == 0 .AND. qplda ) WRITE (io,'(16I5)') ii ELSE ii(1)=1 - IF( mpime == 0 ) WRITE (io) ii + IF( mpime == 0 .AND. qplda ) WRITE (io) ii ENDIF ! WRITE(6,'(16I5)') ii @@ -235,12 +265,12 @@ SUBROUTINE compute_gw( use_gmaps ) IF( mpime == 0 ) THEN ! IF (t_form) THEN - WRITE (io,'(3E26.18)') a1, a2, a3 + IF( qplda) WRITE (io,'(3E26.18)') a1, a2, a3 ELSE IF (t_single) THEN - WRITE (io) a1_s, a2_s, a3_s + IF( qplda) WRITE (io) a1_s, a2_s, a3_s ELSE - WRITE (io) a1, a2, a3 + IF( qplda) WRITE (io) a1, a2, a3 ENDIF ENDIF ! @@ -265,17 +295,17 @@ SUBROUTINE compute_gw( use_gmaps ) WRITE(6,*)'nrot=',nsym WRITE(6,'(3E26.18)') (((float(s(i,j,k)),j=1,3),i=1,3),k=1,nsym) IF (t_form) THEN - WRITE (io,'(I2)') nsym - WRITE (io,'(3E26.18)') (((float(s(i,j,k)),j=1,3),i=1,3),k=1,nsym) + IF( qplda) WRITE (io,'(I2)') nsym + IF( qplda) WRITE (io,'(3E26.18)') (((float(s(i,j,k)),j=1,3),i=1,3),k=1,nsym) IF (ii(3) == 1) THEN ! READ (10,1020) ((VOFFSET(I,J),I=1,3),J=1,NOP) ! WRITE (6,'(//" Run program CNVNSY to convert QPLDA file first.")') CALL errore('pw2gw','non-symmorphic translation vectors',ii(3)) ENDIF ELSE - WRITE (io) nsym + IF( qplda) WRITE (io) nsym IF (t_single) THEN - WRITE (io) (((float(s(i,j,k)),j=1,3),i=1,3),k=1,nsym) + IF( qplda) WRITE (io) (((float(s(i,j,k)),j=1,3),i=1,3),k=1,nsym) ELSE WRITE (io) (((dble(s(i,j,k)),j=1,3),i=1,3),k=1,nsym) ENDIF @@ -397,13 +427,13 @@ SUBROUTINE compute_gw( use_gmaps ) IF (t_form) THEN IF( mpime == 0 ) THEN - WRITE (io,'(I12)') igwxx - WRITE (io,'(3I5)') (in1_tmp(ig),in2_tmp(ig),in3_tmp(ig),ig=1,igwxx) + IF( qplda) WRITE (io,'(I12)') igwxx + IF( qplda) WRITE (io,'(3I5)') (in1_tmp(ig),in2_tmp(ig),in3_tmp(ig),ig=1,igwxx) ENDIF ELSE IF( mpime == 0 ) THEN - WRITE (io) igwxx - WRITE (io) (in1_tmp(ig),in2_tmp(ig),in3_tmp(ig),ig=1,igwxx) + IF( qplda) WRITE (io) igwxx + IF( qplda) WRITE (io) (in1_tmp(ig),in2_tmp(ig),in3_tmp(ig),ig=1,igwxx) ENDIF ENDIF ! @@ -430,95 +460,99 @@ SUBROUTINE compute_gw( use_gmaps ) IF( mpime == 0 ) THEN IF (t_form) THEN - WRITE (io,'(I12)') nkpt - WRITE (io,'(3E26.18)') ((xk_s(i,ik),i=1,3),ik=1,nkpt) + IF( qplda) WRITE (io,'(I12)') nkpt + IF( qplda) WRITE (io,'(3E26.18)') ((xk_s(i,ik),i=1,3),ik=1,nkpt) ELSE - WRITE (io) nkpt - WRITE(6,*) 'nkpt',nkpt + IF( qplda) WRITE (io) nkpt + IF( qplda) WRITE(6,*) 'nkpt',nkpt IF(t_single) THEN - WRITE (io) ((xk_s(i,ik),i=1,3),ik=1,nkpt) + IF( qplda) WRITE (io) ((xk_s(i,ik),i=1,3),ik=1,nkpt) ELSE - WRITE (io) ((xk(i,ik),i=1,3),ik=1,nkpt) + IF( qplda) WRITE (io) ((xk(i,ik),i=1,3),ik=1,nkpt) ENDIF ENDIF WRITE(6,'(1x,3f10.6)') ( (xk_s(i,ik),i=1,3), ik=1,nkpt) ENDIF + + IF( vkb) THEN ! -------------------------- ! vkb0 ! -------------------------- - DO ik=1,nkpt - npw = ngk(ik) - WRITE(15,*) "npw", npw - ALLOCATE(vkb0(1:npw)) - - size_tab=size(tab,1) - size_tab_d2y=size(tab_d2y,1) - - ALLOCATE(vec_tab(1:size_tab)) - if(.not.allocated(vec_tab_d2y)) ALLOCATE(vec_tab_d2y(1:size_tab_d2y)) - - DO nt = 1, ntyp - DO nb = 1, upf(nt)%nbeta - vkb0(:) = 0.0_dp - vec_tab(:) = 0.0_dp - vec_tab_d2y(:) = 0.0_dp - vec_tab(:) = tab(:,nb,nt) - IF(spline_ps) vec_tab_d2y(:) = tab_d2y(:,nb,nt) - CALL gen_us_vkb0(ik,npw,vkb0,size_tab,vec_tab,spline_ps,vec_tab_d2y) - WRITE(15,*) "---------------DEBUG-VKB0----------------------" - WRITE(15,*) "ik= ", ik - WRITE(15,*) "nt= ", nt - WRITE(15,*) "nb= ", nb - WRITE(15,*) "l= ", upf(nt)%lll(nb) - WRITE (15,'(8f15.9)') vkb0 - WRITE(15,*) "--------------END-DEBUG------------------------" -! WRITE(io) vkb0 + DO ik=1,nkpt + npw = ngk(ik) + WRITE(15,*) "npw", npw + ALLOCATE(vkb0(1:npw)) + + size_tab=size(tab,1) + size_tab_d2y=size(tab_d2y,1) + + ALLOCATE(vec_tab(1:size_tab)) + if(.not.allocated(vec_tab_d2y)) ALLOCATE(vec_tab_d2y(1:size_tab_d2y)) + + DO nt = 1, ntyp + DO nb = 1, upf(nt)%nbeta + vkb0(:) = 0.0_dp + vec_tab(:) = 0.0_dp + vec_tab_d2y(:) = 0.0_dp + vec_tab(:) = tab(:,nb,nt) + IF(spline_ps) vec_tab_d2y(:) = tab_d2y(:,nb,nt) + CALL gen_us_vkb0(ik,npw,vkb0,size_tab,vec_tab,spline_ps,vec_tab_d2y) + WRITE(15,*) "---------------DEBUG-VKB0----------------------" + WRITE(15,*) "ik= ", ik + WRITE(15,*) "nt= ", nt + WRITE(15,*) "nb= ", nb + WRITE(15,*) "l= ", upf(nt)%lll(nb) + WRITE (15,'(8f15.9)') vkb0 + WRITE(15,*) "--------------END-DEBUG------------------------" + ! WRITE(io) vkb0 + ENDDO ENDDO - ENDDO - DEALLOCATE(vkb0) - DEALLOCATE(vec_tab) - IF(allocated(vec_tab_d2y)) DEALLOCATE(vec_tab_d2y) + DEALLOCATE(vkb0) + DEALLOCATE(vec_tab) + IF(allocated(vec_tab_d2y)) DEALLOCATE(vec_tab_d2y) - ENDDO + ENDDO !--------------------------- ! djl !--------------------------- - DO ik=1,nkpt - npw = ngk(ik) + DO ik=1,nkpt + npw = ngk(ik) - ALLOCATE(djl(1:npw)) - - size_tab=size(tab,1) - size_tab_d2y=size(tab_d2y,1) - - ALLOCATE(vec_tab(1:size_tab)) - IF(.not. allocated(vec_tab_d2y)) ALLOCATE(vec_tab_d2y(1:size_tab_d2y)) - DO nt = 1, ntyp - DO nb = 1, upf(nt)%nbeta - djl(:) = 0.0_dp - vec_tab(:) = 0.0_dp - vec_tab_d2y(:) = 0.0_dp - vec_tab(:) = tab(:,nb,nt) - IF(spline_ps) vec_tab_d2y(:) = tab_d2y(:,nb,nt) - CALL gen_us_djl(ik,npw,djl,size_tab,vec_tab,spline_ps,vec_tab_d2y) -! WRITE(0,*) "---------------DEBUG-----------------------" -! WRITE(0,*) "spline: ", spline_ps -! WRITE(0,*) "ik= ", ik -! WRITE(0,*) "nt= ", nt -! WRITE(0,*) "nb= ", nb -! WRITE(0,*) "l= ", upf(nt)%lll(nb) -! WRITE (0,'(8f15.9)') djl -! WRITE(0,*) "--------------END-DEBUG------------------------" -! WRITE(io) djl + ALLOCATE(djl(1:npw)) + + size_tab=size(tab,1) + size_tab_d2y=size(tab_d2y,1) + + ALLOCATE(vec_tab(1:size_tab)) + IF(.not. allocated(vec_tab_d2y)) ALLOCATE(vec_tab_d2y(1:size_tab_d2y)) + DO nt = 1, ntyp + DO nb = 1, upf(nt)%nbeta + djl(:) = 0.0_dp + vec_tab(:) = 0.0_dp + vec_tab_d2y(:) = 0.0_dp + vec_tab(:) = tab(:,nb,nt) + IF(spline_ps) vec_tab_d2y(:) = tab_d2y(:,nb,nt) + CALL gen_us_djl(ik,npw,djl,size_tab,vec_tab,spline_ps,vec_tab_d2y) + ! WRITE(0,*) "---------------DEBUG-----------------------" + ! WRITE(0,*) "spline: ", spline_ps + ! WRITE(0,*) "ik= ", ik + ! WRITE(0,*) "nt= ", nt + ! WRITE(0,*) "nb= ", nb + ! WRITE(0,*) "l= ", upf(nt)%lll(nb) + ! WRITE (0,'(8f15.9)') djl + ! WRITE(0,*) "--------------END-DEBUG------------------------" + ! WRITE(io) djl + ENDDO ENDDO - ENDDO - DEALLOCATE(djl) - DEALLOCATE(vec_tab) - IF(allocated(vec_tab_d2y)) DEALLOCATE(vec_tab_d2y) + DEALLOCATE(djl) + DEALLOCATE(vec_tab) + IF(allocated(vec_tab_d2y)) DEALLOCATE(vec_tab_d2y) - ENDDO + ENDDO + + ENDIF !vkb !----------------------- !----------------------- @@ -530,10 +564,10 @@ SUBROUTINE compute_gw( use_gmaps ) eig(:,:) = et(:,1:nkpt)*0.5d0 IF (t_form) THEN - IF( mpime == 0 ) WRITE (io,'(I12)') n - IF( mpime == 0 ) WRITE (io,'(3E26.18)') ((eig(i,ik),ik=1,nkpt),i=1,n) + IF( mpime == 0 .AND. qplda) WRITE (io,'(I12)') n + IF( mpime == 0 .AND. qplda) WRITE (io,'(3E26.18)') ((eig(i,ik),ik=1,nkpt),i=1,n) ELSE - IF( mpime == 0 ) WRITE (io) n + IF( mpime == 0 .AND. qplda) WRITE (io) n IF(t_single) THEN DO ik=1,nkpt DO i=1,n @@ -541,10 +575,10 @@ SUBROUTINE compute_gw( use_gmaps ) ENDDO ENDDO WRITE(6,*) 'nbndsi=',n - IF( mpime == 0 ) WRITE (io) ((eig_s(i,ik),ik=1,nkpt),i=1,n) + IF( mpime == 0 .AND. qplda) WRITE (io) ((eig_s(i,ik),ik=1,nkpt),i=1,n) ELSE WRITE(6,*) 'nbndsi=',n - IF( mpime == 0 ) WRITE (io) ((eig(i,ik),ik=1,nkpt),i=1,n) + IF( mpime == 0 .AND. qplda) WRITE (io) ((eig(i,ik),ik=1,nkpt),i=1,n) ENDIF ENDIF @@ -567,12 +601,12 @@ SUBROUTINE compute_gw( use_gmaps ) IF( mpime == 0 ) THEN IF (t_form) THEN - WRITE (io,'(3E26.18)') ((focc(i,ik), ik=1,nkpt), i=1,n) + IF( qplda) WRITE (io,'(3E26.18)') ((focc(i,ik), ik=1,nkpt), i=1,n) ELSE IF(t_single) THEN - WRITE (io) ((focc_s(i,ik),ik=1,nkpt),i=1,n) + IF( qplda) WRITE (io) ((focc_s(i,ik),ik=1,nkpt),i=1,n) ELSE - WRITE (io) ((focc(i,ik),ik=1,nkpt),i=1,n) + IF( qplda) WRITE (io) ((focc(i,ik),ik=1,nkpt),i=1,n) ENDIF ENDIF ENDIF @@ -590,9 +624,9 @@ SUBROUTINE compute_gw( use_gmaps ) ! ! omax = nbnd*6 - omax = 600 + omax = floor((omegamax - omegamin)/d_omega) WRITE(6,*) 'io sono omax = ', omax - alpha = omegamax/omax + alpha = d_omega WRITE(6,*) 'alpha = ', alpha halfalpha= alpha*.5 WRITE(6,*) 'halfalpha = ', halfalpha @@ -606,7 +640,7 @@ SUBROUTINE compute_gw( use_gmaps ) ALLOCATE( pp1(nkpt,omax+1), pp2(nkpt,omax+1), pp3(nkpt,omax+1), omegatt(omax+1) ) DO o = 1, omax + 1 - omeg(o) = (o-1)*alpha + omeg(o) = omegamin + (o-1)*alpha ENDDO epstx(:)=0.0 @@ -621,7 +655,16 @@ SUBROUTINE compute_gw( use_gmaps ) WRITE(6,*) pp1(1,1), epstx(1) - ALLOCATE ( c0(igwx), c0_s(igwx), kpg(3,igwx), c0_m(igwx,n), c0_tmp(igwxx) ) + npol=1 + IF (nspin == 4) npol=2 + write(*,*) "igwx, npw", igwx, npw + ALLOCATE ( c0(igwx*npol), c0_s(igwx*npol), kpg(3,igwx), c0_m(igwx*npol,n), & + c0_tmp(igwxx*npol) ) + c0 = 0.0 + c0_s = 0.0 + kpg = 0.0 + c0_m = 0.0 + c0_tmp = 0.0 !IF (gamma_only) ALLOCATE ( c0_gamma(2*igwx-1), c0_gamma_s(2*igwx-1) ) CALL cryst_to_cart (nks, xk, bg, +1) IF (ionode) WRITE(6,*) 'Costruisco le psi ed il matrixelements' @@ -670,11 +713,12 @@ SUBROUTINE compute_gw( use_gmaps ) c0 (:) = 0.d0 DO ig=1,npw c0(igk_k(ig,ik)) = evc(ig,i) + IF (nspin == 4) c0(igk_k(ig,ik)+igwx) = evc(ig+npwx,i) ENDDO c0_m(:,i)=c0(:) c0_tmp = c0_tmp_dp - IF( mpime == 0 ) WRITE(io) c0_tmp ! c0_s + IF( mpime == 0 .AND. qplda) WRITE(io) c0_tmp ! c0_s DEALLOCATE( c0_tmp_dp ) @@ -686,7 +730,7 @@ SUBROUTINE compute_gw( use_gmaps ) ! k + g thet must be in 2piba units kpg(:,:) = 0.d0 DO ig=1,npw - kpg(:,igk_k(ig,ik))= xk_s(:,ik)+g(:,igk_k(ig,ik)) + kpg(:,igk_k(ig,ik))= xk(:,ik)+g(:,igk_k(ig,ik)) ENDDO DO iband1 = 1,n @@ -704,6 +748,12 @@ SUBROUTINE compute_gw( use_gmaps ) rhotwx(1) = rhotwx(1) + xkgk(1) * ctemp rhotwx(2) = rhotwx(2) + xkgk(2) * ctemp rhotwx(3) = rhotwx(3) + xkgk(3) * ctemp + IF (nspin == 4) THEN + ctemp= conjg(c0_m(ig+igwx,iband1))*c0_m(ig+igwx,iband2) + rhotwx(1) = rhotwx(1) + xkgk(1) * ctemp + rhotwx(2) = rhotwx(2) + xkgk(2) * ctemp + rhotwx(3) = rhotwx(3) + xkgk(3) * ctemp + ENDIF ENDDO CALL mp_sum( rhotwx, world_comm ) @@ -736,41 +786,43 @@ SUBROUTINE compute_gw( use_gmaps ) ! CALCULATE POTENTIAL MATRIX ELEMNTS - OPEN (UNIT=313,FILE="vxcdiag.dat",STATUS="UNKNOWN") - WRITE(313,*) "# K BND " - ALLOCATE ( vxc(dfftp%nnr,nspin) ) - CALL v_xc (rho, rho_core, rhog_core, etxc, vtxc, vxc) - DO ik=1,nkpt - npw = ngk(ik) - CALL davcio( evc, 2*nwordwfc, iunwfc, ik, -1 ) - DO iband1 = 1, nbnd - psic(:) = (0.d0, 0.d0) - DO ig = 1, npw - psic(nl(igk_k(ig,ik))) = evc(ig,iband1) - ENDDO - - CALL invfft ('Dense', psic, dfftp) - vxcdiag = 0.0d0 - !norma = 0.0d0 - DO ir = 1, dfftp%nnr - vxcdiag = vxcdiag + vxc(ir,nspin) * & - ( dble(psic (ir) ) **2 + aimag(psic (ir) ) **2) - ! norma = norma + ( DBLE(psic (ir) ) **2 + AIMAG(psic (ir) ) **2) / (nr1*nr2*nr3) - ENDDO - ! PG: this is the correct integral - 27/8/2010 - vxcdiag = vxcdiag * rytoev / (dfftp%nr1*dfftp%nr2*dfftp%nr3) - CALL mp_sum( vxcdiag, world_comm ) !, intra_pool_comm ) - ! ONLY FOR DEBUG! - !IF (norma /= 1.0) THEN - ! WRITE(*,*) "norma =", norma - ! WRITE(*,*) "nrxx =", nrxx - ! STOP - !ENDIF - WRITE(313,"(i1,2x,i1,4x,f18.14)") ik, iband1, vxcdiag - ENDDO - ENDDO - DEALLOCATE ( vxc ) - CLOSE (313) + IF( vxcdiag_f) THEN + OPEN (UNIT=313,FILE="vxcdiag.dat",STATUS="UNKNOWN") + WRITE(313,*) "# K BND " + ALLOCATE ( vxc(dfftp%nnr,nspin) ) + CALL v_xc (rho, rho_core, rhog_core, etxc, vtxc, vxc) + DO ik=1,nkpt + npw = ngk(ik) + CALL davcio( evc, 2*nwordwfc, iunwfc, ik, -1 ) + DO iband1 = 1, nbnd + psic(:) = (0.d0, 0.d0) + DO ig = 1, npw + psic(nl(igk_k(ig,ik))) = evc(ig,iband1) + ENDDO + + CALL invfft ('Dense', psic, dfftp) + vxcdiag = 0.0d0 + !norma = 0.0d0 + DO ir = 1, dfftp%nnr + vxcdiag = vxcdiag + vxc(ir,nspin) * & + ( dble(psic (ir) ) **2 + aimag(psic (ir) ) **2) + ! norma = norma + ( DBLE(psic (ir) ) **2 + AIMAG(psic (ir) ) **2) / (nr1*nr2*nr3) + ENDDO + ! PG: this is the correct integral - 27/8/2010 + vxcdiag = vxcdiag * rytoev / (dfftp%nr1*dfftp%nr2*dfftp%nr3) + CALL mp_sum( vxcdiag, world_comm ) !, intra_pool_comm ) + ! ONLY FOR DEBUG! + !IF (norma /= 1.0) THEN + ! WRITE(*,*) "norma =", norma + ! WRITE(*,*) "nrxx =", nrxx + ! STOP + !ENDIF + WRITE(313,"(i1,2x,i1,4x,f18.14)") ik, iband1, vxcdiag + ENDDO + ENDDO + DEALLOCATE ( vxc ) + CLOSE (313) + ENDIF ! ! @@ -830,7 +882,7 @@ SUBROUTINE compute_gw( use_gmaps ) DEALLOCATE ( epsx, epsy, epsz ) DEALLOCATE ( epstx, epsty, epstz ) ! - IF( mpime == 0 ) CLOSE(io) + IF( mpime == 0 .AND. qplda) CLOSE(io) IF( mpime == 0 ) CLOSE(90) ! END SUBROUTINE compute_gw diff --git a/dev-tools/release.sh b/dev-tools/release.sh index 3371607f9f..2dd296e9c3 100755 --- a/dev-tools/release.sh +++ b/dev-tools/release.sh @@ -1,7 +1,7 @@ #!/bin/sh -x version=6.2.1 -revision=14035 +revision=14038 user=giannozz tempdir=$HOME/tempdir