Skip to content

GitLab

Commit b610e50f authored by Marco Govoni's avatar Marco Govoni
Browse files

Init 836.

parents
Expand all Hide whitespace changes
Inline Side-by-side
Coulomb_kernel/Makefile 0 → 100644
View file @ b610e50f
# Makefile for Coulomb_kernel
include ../../make.inc
# location of needed modules
MODFLAGS= $(MOD_FLAG)../../iotk/src $(MOD_FLAG)../../Modules $(MOD_FLAG)../../LAXlib \
$(MOD_FLAG)../../FFTXlib $(MOD_FLAG)../../PW/src \
$(MOD_FLAG).
IFLAGS=
COULOMB_KERNEL_OBJS = \
store_sqvc.o
PWOBJS = ../../PW/src/libpw.a
QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a
TLDEPS= bindir mods libs pw
all : tldeps lib_coulomb_kernel.a
tldeps:
test -n "$(TLDEPS)" && ( cd ../.. ; $(MAKE) $(MFLAGS) $(TLDEPS) || exit 1) || :
lib_coulomb_kernel.a : $(COULOMB_KERNEL_OBJS)
$(AR) $(ARFLAGS) $@ $?
$(RANLIB) $@
clean :
- /bin/rm -f *.x *.o *.a *~ *_tmp.f90 *.d *.mod *.i *.L
- /bin/rm -f lib_coulomb_kernel.a
title :
@echo " "
@echo "#######################"
@echo "### Coulomb_kernel ####"
@echo "#######################"
@echo " "
include make.depend
Coulomb_kernel/make.depend 0 → 100644
View file @ b610e50f
store_sqvc.o : ../../Modules/cell_base.o
store_sqvc.o : ../../Modules/constants.o
store_sqvc.o : ../../Modules/io_global.o
store_sqvc.o : ../../Modules/kind.o
store_sqvc.o : ../../Modules/recvec.o
store_sqvc.o : ../../PW/src/exx.o
store_sqvc.o : ../../PW/src/pwcom.o
store_vspx.o : ../../Modules/cell_base.o
store_vspx.o : ../../Modules/constants.o
store_vspx.o : ../../Modules/kind.o
store_vspx.o : ../../Modules/recvec.o
store_vspx.o : ../../PW/src/exx.o
store_vspx.o : ../../PW/src/pwcom.o
Coulomb_kernel/store_sqvc.f90 0 → 100644
View file @ b610e50f
!
! Copyright (C) 2015-2016 M. Govoni
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
! This file is part of WEST.
!
! Contributors to this file:
! Marco Govoni
!
!-----------------------------------------------------------------------
SUBROUTINE store_sqvc(sqvc_tmp,numg,singularity_removal_mode,div)
!-----------------------------------------------------------------------
!
! This routine computes results of applying sqVc to a potential
! associated with vector q. Coulomb cutoff technique can be used
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE constants, ONLY : fpi, e2, tpi, pi
USE cell_base, ONLY : tpiba2, tpiba,omega,at,alat,bg
USE gvect, ONLY : g, gstart
USE coulomb_vcut_module, ONLY : vcut_init, vcut_type, vcut_info, &
vcut_get, vcut_spheric_get, vcut_destroy
USE mp, ONLY : mp_sum
USE mp_global, ONLY : intra_bgrp_comm
USE mp_world, ONLY : mpime,world_comm,nproc
USE gvecw, ONLY : ecutwfc
USE control_flags, ONLY : gamma_only
USE random_numbers, ONLY : randy
!
IMPLICIT NONE
!
! I/O
!
INTEGER,INTENT(IN) :: numg
REAL(DP),INTENT(OUT) :: sqvc_tmp(numg)
INTEGER,INTENT(IN) :: singularity_removal_mode
REAL(DP),INTENT(OUT) :: div
!
! Workspace
!
REAL(DP) :: gnorm2,q(3),x,nq(3),ecutvcut,atws(3,3),dq(3),xq(3),alpha
INTEGER :: ig, ipol, iq1, iq2, iq3, i1, i2, i3
LOGICAL :: on_double_grid
REAL(DP),PARAMETER :: eps=1.d-6
REAL(DP) :: grid_factor = 8.d0/7.d0
TYPE(vcut_type) :: vcut
LOGICAL :: try_ort_div=.TRUE.,i_am_ort
REAL(DP) :: prod(3,3), qhelp, edge(3), qbz(3), rand, qmo, vbz, vhelp
REAL(DP) :: intcounter
!
CALL start_clock('storesqvc')
!
! ======
! BODY
! ======
!
SELECT CASE(singularity_removal_mode)
!
CASE(1)
!
! In this case we use the spherical region
!
sqvc_tmp(1)=0._DP
DO ig = gstart,numg
gnorm2 = ( g(1,ig)*g(1,ig) + g(2,ig)*g(2,ig) + g(3,ig)*g(3,ig) ) * tpiba2
sqvc_tmp(ig) = SQRT(e2*fpi/gnorm2)
ENDDO
!
CASE(2)
!
! In this case we use Gygi-Baldereschi method
!
nq(1)=1._DP
nq(2)=1._DP
nq(3)=1._DP
!
sqvc_tmp(1)=0._DP
DO ig = gstart,numg
q(:) = g(:,ig)
gnorm2 = SUM(q(:)**2) * tpiba2
on_double_grid = .TRUE.
DO ipol = 1,3
x = 0.5_DP*( q(1)*at(1,ipol)+q(2)*at(2,ipol)+q(3)*at(3,ipol) )*nq(ipol)
on_double_grid = on_double_grid .AND. (ABS(x-NINT(x))<eps)
ENDDO
IF(on_double_grid) THEN
sqvc_tmp(ig)=0._DP
ELSE
sqvc_tmp(ig)=SQRT(e2*fpi*grid_factor/gnorm2)
ENDIF
ENDDO
!
CASE(3)
!
! In this case we use CUT_WS
!
!
nq(1)=1._DP
nq(2)=1._DP
nq(3)=1._DP
ecutvcut = 0.7_DP
!
! build the superperiodicity direct lattice
!
atws = alat * at
!
DO ipol=1,3
atws(:,ipol) = atws(:,ipol) * nq(ipol)
ENDDO
!
CALL vcut_init( vcut, atws, ecutvcut )
!
DO ig = 1,numg
!
q(:) = g(:,ig) * tpiba
!
sqvc_tmp( ig ) = DSQRT( vcut_get(vcut,q) )
!
ENDDO
!
CALL vcut_destroy(vcut)
!
END SELECT
!
! ======
! HEAD
! ======
!
SELECT CASE(singularity_removal_mode)
!
CASE(1)
!
! In this case we use the spherical region
!
div = ( (6._DP * pi * pi / omega )**(1._DP/3._DP) ) / ( 2._DP * pi * pi ) * fpi * e2
WRITE(stdout,"(5X,'Spherical div = ',es14.6)") div
!
IF( try_ort_div ) THEN
!
! Redefine div according to orthorombic cell
!
prod = 0._DP
DO i1 = 1, 3
DO i2 = 1, 3
DO i3 = 1, 3
prod(i1,i2) = prod(i1,i2) + bg(i3,i1) * bg(i3,i2) * tpiba2
ENDDO
ENDDO
ENDDO
!
i_am_ort = .TRUE.
DO i1 = 1, 3
DO i2 = 1, 3
IF( i1 == i2 ) CYCLE
IF( ABS( prod(i1,i2)) > 1.d-8 ) THEN
WRITE(stdout,"(5X,'Warning: non orthorombic RL')")
i_am_ort = .FALSE.
ENDIF
ENDDO
ENDDO
!
IF ( i_am_ort ) THEN
!
edge(1) = SQRT(prod(1,1)) / 2._DP
edge(2) = SQRT(prod(2,2)) / 2._DP
edge(3) = SQRT(prod(3,3)) / 2._DP
!
qhelp = MIN( edge(1),edge(2),edge(3) )
vbz = tpi**3 / omega
vhelp = fpi / 3._DP * qhelp**3
!
rand=randy(mpime)
div = 0._DP
intcounter = 0
!
DO i1 = 1, 100000
qmo=0._DP
DO i2 = 1, 3
qbz(i2) = randy() * edge(i2)
qmo = qmo + qbz(i2)**2
ENDDO
qmo = SQRT( qmo )
IF( qmo < qhelp ) CYCLE
div = div + 1._DP/qmo/qmo
intcounter = intcounter + 1._DP
ENDDO
!
div = div * ( vbz - vhelp ) / intcounter
div = div + fpi * qhelp
div = div * fpi * e2 / ( tpi * tpi * tpi )
!
div = div / REAL(nproc,KIND=DP)
CALL mp_sum(div,world_comm)
!
WRITE(stdout,"(5X,'Orthorombic div = ',es14.6)") div
!
ENDIF
!
ENDIF
!
CASE(2)
!
! In this case we use Gygi-Baldereschi method
!
alpha = 10._DP / ecutwfc
!
nq(1)=1._DP
nq(2)=1._DP
nq(3)=1._DP
!
dq(1)=1._DP/nq(1)
dq(2)=1._DP/nq(2)
dq(3)=1._DP/nq(3)
!
div = 0._DP
!
DO iq1=1,1,INT(nq(1))
DO iq2=1,1,INT(nq(2))
DO iq3=1,1,INT(nq(3))
xq(:) = bg(:,1)*(iq1-1)*dq(1) + bg(:,2)*(iq2-1)*dq(2) + bg(:,3)*(iq3-1)*dq(3)
!
DO ig = gstart,numg
q(:) = g(:,ig) + xq(:)
gnorm2 = SUM(q(:)**2) * tpiba2
on_double_grid = .TRUE.
DO ipol = 1,3
x = 0.5_DP*( q(1)*at(1,ipol)+q(2)*at(2,ipol)+q(3)*at(3,ipol) )*nq(ipol)
on_double_grid = on_double_grid .AND. (ABS(x-NINT(x))<eps)
ENDDO
IF(.NOT.on_double_grid) THEN
div = div - EXP( -alpha * gnorm2 ) / gnorm2
ENDIF
ENDDO
!
ENDDO
ENDDO
ENDDO
!
CALL mp_sum( div, intra_bgrp_comm )
!
IF( gamma_only ) THEN
div = div * grid_factor * e2 * fpi / omega * 2._DP
ELSE
div = div * grid_factor * e2 * fpi / omega
ENDIF
!
div = div + nq(1)*nq(2)*nq(3) * e2 / SQRT( alpha * pi )
!
CASE(3)
!
! In this case we use CUT_WS
!
div = 0._DP
!
END SELECT
!
CALL stop_clock('storesqvc')
!
END SUBROUTINE
!
!------------------------------------------------------------------
SUBROUTINE store_vcspecial_H(vc_tmp,numg)
!
! vcspecail_H
!
USE kinds, ONLY : DP
USE cell_base, ONLY : tpiba, at,omega
USE pwcom, ONLY : npw,npwx
USE gvect, ONLY : g,gstart
USE constants, ONLY : fpi, e2, pi
USE exx, ONLY : exxdiv
!
IMPLICIT NONE
!
! I/O
!
INTEGER,INTENT(IN) :: numg
REAL(DP),INTENT(OUT) :: vc_tmp(numg)
!
!Local variables
INTEGER :: ig !Counters
REAL(DP) :: q(3), qq, x
REAL(DP) :: grid_factor = 8.d0/7.d0
REAL(DP) :: eps = 1.d-6
LOGICAL :: on_double_grid
REAL(DP) :: kost
!
kost = e2*fpi
!
IF(gstart==2) vc_tmp(1) = -exxdiv
!
!$OMP PARALLEL SHARED(numg,g,tpiba,at,vc_tmp,kost,gstart) PRIVATE(ig,q,qq,on_double_grid,x)
!$OMP DO
DO ig=gstart,numg
!
q(:)= g(:,ig)
!
q = q * tpiba
!
qq = SUM(q(:)**2)
!
on_double_grid = .TRUE.
x= 0.5_DP/tpiba*(q(1)*at(1,1)+q(2)*at(2,1)+q(3)*at(3,1))!*nq1
on_double_grid = on_double_grid .AND. (ABS(x-NINT(x))<eps)
x= 0.5_DP/tpiba*(q(1)*at(1,2)+q(2)*at(2,2)+q(3)*at(3,2))!*nq2
on_double_grid = on_double_grid .AND. (ABS(x-NINT(x))<eps)
x= 0.5_DP/tpiba*(q(1)*at(1,3)+q(2)*at(2,3)+q(3)*at(3,3))!*nq3
on_double_grid = on_double_grid .AND. (ABS(x-NINT(x))<eps)
!
IF( on_double_grid ) THEN
!
vc_tmp(ig) = -kost / qq
!
ELSE
!
vc_tmp(ig)= kost/qq / 7._DP
!
ENDIF
!
ENDDO
!$OMP ENDDO
!$OMP END PARALLEL
!
END SUBROUTINE
!
!
!
SUBROUTINE store_sqvc_sphcut(sqvc_tmp,numg,rcut)
!
! This routine computes results of applying sqVc to a potential
! associated with vector q. Coulomb cutoff technique can be used
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE constants, ONLY : fpi, e2, tpi
USE cell_base, ONLY : tpiba2, tpiba
USE gvect, ONLY : g, gstart
!
IMPLICIT NONE
!
! I/O
!
INTEGER,INTENT(IN) :: numg
REAL(DP),INTENT(OUT) :: sqvc_tmp(numg)
REAL(DP),INTENT(IN) :: rcut
!
! Workspace
!
REAL(DP) :: gnorm2,gnorm
INTEGER :: ig
!
!
CALL start_clock('storesqvc')
!
!
IF(gstart==2) sqvc_tmp(1)=SQRT(tpi*e2*rcut*rcut)
!
!$OMP PARALLEL private(ig,gnorm2,gnorm)
!$OMP DO
DO ig = gstart,numg
gnorm2 = ( g(1,ig)*g(1,ig) + g(2,ig)*g(2,ig) + g(3,ig)*g(3,ig) ) * tpiba2
gnorm = SQRT(gnorm2)
sqvc_tmp(ig) = SQRT(e2*fpi/gnorm2 * (1._DP-COS(gnorm*rcut)))
ENDDO
!$OMP ENDDO
!$OMP END PARALLEL
!
CALL stop_clock('storesqvc')
!
END SUBROUTINE
DFPT_kernel/Makefile 0 → 100644
View file @ b610e50f
# Makefile for DFPT_kernel
include ../../make.inc
# location of needed modules
MODFLAGS= $(MOD_FLAG)../../iotk/src $(MOD_FLAG)../../Modules $(MOD_FLAG)../../LAXlib \
$(MOD_FLAG)../../FFTXlib $(MOD_FLAG)../../PW/src \
$(MOD_FLAG)../Modules \
$(MOD_FLAG)../Tools \
$(MOD_FLAG)../FFT_kernel \
$(MOD_FLAG).
IFLAGS=
DFPT_KERNEL_OBJS = \
apply_sternheimerop_to_m_wfcs.o \
precondition_m_wfcts.o \
linsolve_sternheimer_m_wfcts.o \
dfpt.o
PWOBJS = ../../PW/src/libpw.a
QEMODS = ../../Modules/libqemod.a ../../FFTXlib/libqefft.a ../../LAXlib/libqela.a
TLDEPS= bindir mods libs pw
all : title tldeps lib_dfpt_kernel.a
tldeps:
test -n "$(TLDEPS)" && ( cd ../.. ; $(MAKE) $(MFLAGS) $(TLDEPS) || exit 1) || :
lib_dfpt_kernel.a : $(DFPT_KERNEL_OBJS)
$(AR) $(ARFLAGS) $@ $?
$(RANLIB) $@
clean :
- /bin/rm -f *.x *.o *.a *~ *_tmp.f90 *.d *.mod *.i *.L
- /bin/rm -f lib_dfpt_kernel.a
title :
@echo " "
@echo "####################"
@echo "### DFPT_kernel ####"
@echo "####################"
@echo " "
include make.depend
DFPT_kernel/apply_sternheimerop_to_m_wfcs.f90 0 → 100644
View file @ b610e50f
!
! Copyright (C) 2015-2016 M. Govoni
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
! This file is part of WEST.
!
! Contributors to this file:
! Marco Govoni
!
!-----------------------------------------------------------------------
SUBROUTINE apply_sternheimerop_to_m_wfcs(nbndval, psi, hpsi, e, alpha, m)
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE pwcom, ONLY : npw,npwx
USE wvfct, ONLY : nbnd
USE becmod, ONLY : bec_type, becp, calbec
USE uspp, ONLY : nkb, vkb
USE mp_global, ONLY : intra_pool_comm,intra_bgrp_comm
USE mp, ONLY : mp_sum
USE wavefunctions_module, ONLY : evc
USE noncollin_module, ONLY : npol, noncolin
USE westcom, ONLY : l_kinetic_only
!
IMPLICIT NONE
!
! I/O
!
INTEGER, INTENT(IN) :: nbndval, m
! input: the number of val bands
! input: the number of bands
REAL(DP), INTENT(IN) :: e (m)
COMPLEX(DP), INTENT(IN) :: alpha
! input: the eigenvalue
COMPLEX(DP), INTENT(IN) :: psi (npwx*npol,m)
COMPLEX(DP), INTENT(OUT) :: hpsi (npwx*npol,m)
! input: the vector
! output: the operator applied to the vector
!
! Workspace
!
INTEGER :: ibnd, ig
COMPLEX(DP) :: za
!
CALL start_clock ('stern')
!
! compute the product of the hamiltonian with the h vector
!
hpsi=(0.0_DP,0.0_DP)
!
IF(l_kinetic_only) THEN
CALL k_psi( npwx, npw, m, psi, hpsi )
ELSE
CALL h_psi( npwx, npw, m, psi, hpsi )
ENDIF
!
! then we compute the operator H-epsilon S
!
DO ibnd=1,m
za = CMPLX( -e(ibnd), 0._DP, KIND=DP )
CALL ZAXPY(npw,za,psi(1,ibnd),1,hpsi(1,ibnd),1)
ENDDO
IF(npol==2) THEN
DO ibnd=1,m
za = CMPLX( -e(ibnd), 0._DP, KIND=DP )
CALL ZAXPY(npw,za,psi(npwx+1,ibnd),1,hpsi(npwx+1,ibnd),1)
ENDDO
ENDIF
!
CALL apply_alpha_pv_to_m_wfcs(nbndval,m,psi,hpsi,alpha)
!
CALL stop_clock ('stern')
!
END SUBROUTINE
DFPT_kernel/dfpt.f90 0 → 100644
View file @ b610e50f
!
! Copyright (C) 2015-2016 M. Govoni
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
! This file is part of WEST.
!
! Contributors to this file:
! Marco Govoni
!
!-----------------------------------------------------------------------
SUBROUTINE dfpt (m,dvg,dng,tr2)
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE wvfct, ONLY : nbnd,g2kin,et
USE fft_base, ONLY : dfftp,dffts
USE gvect, ONLY : nl,gstart,ig_l2g
USE wavefunctions_module, ONLY : evc,psic
USE gvecs, ONLY : ngms
USE gvecw, ONLY : gcutw
USE mp, ONLY : mp_sum,mp_barrier,mp_bcast
USE mp_global, ONLY : inter_image_comm,inter_pool_comm,my_image_id
USE fft_at_gamma, ONLY : single_fwfft_gamma,single_invfft_gamma,double_fwfft_gamma,double_invfft_gamma
USE fft_at_k, ONLY : single_fwfft_k,single_invfft_k
USE buffers, ONLY : get_buffer
USE noncollin_module, ONLY : noncolin,npol
USE bar, ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
USE pwcom, ONLY : current_spin,wk,nks,nelup,neldw,isk,g,igk_k,ngm,tpiba2,xk,omega,npw,npwx,lsda,nkstot,&