scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

- fix for NLCC contribution to forces and stress 

- more BGL porting
- clean-ups


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3020 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
cavazzon 2006-04-18 07:33:11 +00:00
parent aa6c592052
commit 25b02439aa
19 changed files with 534 additions and 388 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CPV/cg_sub.f90
View File

@ -14,7 +14,7 @@
lambdap, lambda )
use kinds, only: dp
use control_flags, only: iprint, thdyn, tpre, tbuff, iprsta, &
use control_flags, only: iprint, thdyn, tpre, iprsta, &
tfor, tvlocw, taurdr, tprnfor
use control_flags, only: ndr, ndw, nbeg, nomore, tsde, tortho, tnosee, &
tnosep, trane, tranp, tsdp, tcp, tcap, ampre, amprp, tnoseh

172
CPV/chargedensity.f90
View File

@ -32,7 +32,7 @@
! end of module-scope declarations
! ----------------------------------------------
PUBLIC :: rhoofr, fillgrad
PUBLIC :: rhoofr, fillgrad, checkrho
INTERFACE rhoofr
MODULE PROCEDURE rhoofr_fpmd, rhoofr_cp
@ -359,6 +359,7 @@
END SUBROUTINE rhoofr_fpmd
!=----------------------------------------------------------------------=!
!-----------------------------------------------------------------------
SUBROUTINE rhoofr_cp (nfi,c,irb,eigrb,bec,rhovan,rhor,rhog,rhos,enl,ekin)
!-----------------------------------------------------------------------
@ -374,7 +375,7 @@
! e_v = sum_i,ij rho_i,ij d^ion_is,ji
!
USE kinds, ONLY: DP
USE control_flags, ONLY: iprint, tbuff, iprsta, thdyn, tpre, trhor
USE control_flags, ONLY: iprint, iprsta, thdyn, tpre, trhor
USE ions_base, ONLY: nat
USE gvecp, ONLY: ngm
USE gvecs, ONLY: ngs, nps, nms
@ -511,35 +512,24 @@
END DO
CALL invfft('Wave',psis,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
! wavefunctions in unit 21
!
IF(tbuff) WRITE(21,iostat=ios) psis
iss1=ispin(i)
sa1=f(i)/omega
IF (i.NE.n) THEN
iss2=ispin(i+1)
sa2=f(i+1)/omega
iss1 = ispin(i)
sa1 = f(i) / omega
IF ( i .NE. n ) THEN
iss2 = ispin(i+1)
sa2 = f(i+1) / omega
ELSE
iss2=iss1
sa2=0.0
iss2 = iss1
sa2 = 0.0
END IF
DO ir=1,nnrsx
rhos(ir,iss1)=rhos(ir,iss1) + sa1*( DBLE(psis(ir)))**2
rhos(ir,iss2)=rhos(ir,iss2) + sa2*(AIMAG(psis(ir)))**2
END DO
!
! buffer 21
!
IF(tbuff) THEN
IF(ios.NE.0) CALL errore(' rhoofr',' error in writing unit 21',ios)
ENDIF
DO ir = 1, nnrsx
rhos(ir,iss1) = rhos(ir,iss1) + sa1 * ( DBLE(psis(ir)))**2
rhos(ir,iss2) = rhos(ir,iss2) + sa2 * (AIMAG(psis(ir)))**2
END DO
!
END DO
!
IF(tbuff) REWIND 21
!
! smooth charge in g-space is put into rhog(ig)
!
IF(nspin.EQ.1)THEN
@ -597,39 +587,15 @@
rhor(ir,isdw)=AIMAG(psi(ir))
END DO
ENDIF
IF (dft_is_meta()) CALL kedtauofr_meta(c, psi, psis) ! METAGGA
!
IF(iprsta.GE.3)THEN
DO iss=1,nspin
rsumg(iss)=omega*DBLE(rhog(1,iss))
rsumr(iss)=SUM(rhor(:,iss))*omega/DBLE(nr1*nr2*nr3)
END DO
IF ( gstart /= 2 ) THEN
!
! in the parallel case, only one processor has G=0 !
!
DO iss=1,nspin
rsumg(iss)=0.0
END DO
END IF
CALL mp_sum( rsumg( 1:nspin ), intra_image_comm )
CALL mp_sum( rsumr( 1:nspin ), intra_image_comm )
IF ( nspin == 1 ) THEN
WRITE( stdout, 10) rsumg(1), rsumr(1)
ELSE
WRITE( stdout, 20) rsumg(1), rsumr(1), rsumg(2), rsumr(2)
ENDIF
ENDIF
!
IF ( dft_is_meta() ) CALL kedtauofr_meta( c, psi, psis ) ! METAGGA
!
! add vanderbilt contribution to the charge density
! drhov called before rhov because input rho must be the smooth part
!
IF (tpre) CALL drhov(irb,eigrb,rhovan,rhog,rhor)
IF ( tpre ) CALL drhov( irb, eigrb, rhovan, rhog, rhor )
!
CALL rhov(irb,eigrb,rhovan,rhog,rhor)
CALL rhov( irb, eigrb, rhovan, rhog, rhor )
rhopr = rhor
@ -639,31 +605,18 @@
!
! here to check the integral of the charge density
!
!
IF( iprsta .GE. 2 ) THEN
CALL checkrho(nnrx,nspin,rhor,rmin,rmax,rsum,rnegsum)
rnegsum=rnegsum*omega/DBLE(nr1*nr2*nr3)
rsum=rsum*omega/DBLE(nr1*nr2*nr3)
WRITE( stdout,'(a,4(1x,f12.6))') &
IF( ( iprsta >= 2 ) .OR. ( nfi == 0 ) .OR. &
( MOD(nfi, iprint) == 0 ) .AND. ( .NOT. tcg ) ) THEN
IF( iprsta >= 2 ) THEN
CALL checkrho( nnrx, nspin, rhor, rmin, rmax, rsum, rnegsum )
rnegsum = rnegsum * omega / DBLE(nr1*nr2*nr3)
rsum = rsum * omega / DBLE(nr1*nr2*nr3)
WRITE( stdout,'(a,4(1x,f12.6))') &
& ' rhoofr: rmin rmax rnegsum rsum ',rmin,rmax,rnegsum,rsum
END IF
!
IF( nfi == 0 .OR. MOD(nfi, iprint) == 0 .AND. .NOT. tcg ) THEN
DO iss=1,nspin
rsumg(iss)=omega*DBLE(rhog(1,iss))
rsumr(iss)=SUM(rhor(:,iss),1)*omega/DBLE(nr1*nr2*nr3)
END DO
IF (gstart.NE.2) THEN
! in the parallel case, only one processor has G=0 !
DO iss=1,nspin
rsumg(iss)=0.0
END DO
END IF
CALL mp_sum( rsumg( 1:nspin ), intra_image_comm )
CALL mp_sum( rsumr( 1:nspin ), intra_image_comm )
CALL sum_charge( rsumg, rsumr )
IF ( nspin == 1 ) THEN
WRITE( stdout, 10) rsumg(1), rsumr(1)
@ -688,7 +641,38 @@
!
RETURN
END SUBROUTINE rhoofr_cp
CONTAINS
!
!
SUBROUTINE sum_charge( rsumg, rsumr )
!
REAL(DP), INTENT(OUT) :: rsumg( : )
REAL(DP), INTENT(OUT) :: rsumr( : )
INTEGER :: iss
!
DO iss=1,nspin
rsumg(iss)=omega*DBLE(rhog(1,iss))
rsumr(iss)=SUM(rhor(:,iss),1)*omega/DBLE(nr1*nr2*nr3)
END DO
IF (gstart.NE.2) THEN
! in the parallel case, only one processor has G=0 !
DO iss=1,nspin
rsumg(iss)=0.0
END DO
END IF
CALL mp_sum( rsumg( 1:nspin ), intra_image_comm )
CALL mp_sum( rsumr( 1:nspin ), intra_image_comm )
RETURN
END SUBROUTINE
!-----------------------------------------------------------------------
END SUBROUTINE rhoofr_cp
!-----------------------------------------------------------------------
@ -758,6 +742,44 @@
END SUBROUTINE fillgrad
!
!----------------------------------------------------------------------
SUBROUTINE checkrho(nnr,nspin,rhor,rmin,rmax,rsum,rnegsum)
!----------------------------------------------------------------------
!
! check \int rho(r)dr and the negative part of rho
!
USE kinds, ONLY: DP
USE mp, ONLY: mp_sum
USE mp_global, ONLY: intra_image_comm
IMPLICIT NONE
INTEGER nnr, nspin
REAL(DP) rhor(nnr,nspin), rmin, rmax, rsum, rnegsum
!
REAL(DP) roe
INTEGER ir, iss
!
rsum =0.0
rnegsum=0.0
rmin =100.
rmax =0.0
DO iss = 1, nspin
DO ir = 1, nnr
roe = rhor(ir,iss)
rsum = rsum + roe
IF ( roe < 0.0 ) rnegsum = rnegsum + roe
rmax = MAX( rmax, roe )
rmin = MIN( rmin, roe )
END DO
END DO
CALL mp_sum( rsum, intra_image_comm )
CALL mp_sum( rnegsum, intra_image_comm )
RETURN
END SUBROUTINE checkrho
!=----------------------------------------------------------------------=!
END MODULE charge_density
!=----------------------------------------------------------------------=!

188
CPV/cplib.f90
View File

@ -501,7 +501,7 @@
! sum_i,ij d^q_i,ij (-i)**l beta_i,i(g)
! e^-ig.r_i < beta_i,j | c_n >}
USE kinds, ONLY: dp
USE control_flags, ONLY: iprint, tbuff
USE control_flags, ONLY: iprint
USE gvecs
USE gvecw, ONLY: ngw
USE cvan, ONLY: ish
@ -546,26 +546,14 @@
!
ci=(0.0,1.0)
!
IF (.NOT.tbuff) THEN
!
psi (:) = (0.d0, 0.d0)
DO ig=1,ngw
psi (:) = (0.d0, 0.d0)
DO ig=1,ngw
psi(nms(ig))=CONJG(c(ig)-ci*ca(ig))
psi(nps(ig))=c(ig)+ci*ca(ig)
END DO
!
CALL invfft('Wave',psi,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
END DO
CALL invfft('Wave',psi,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
!
ELSE
!
! read psi from buffer 21
!
READ(21,iostat=ios) psi
IF(ios.NE.0) CALL errore &
& (' dforce',' error in reading unit 21',ios)
!
ENDIF
!
iss1=ispin(i)
!
! the following avoids a potential out-of-bounds error
@ -796,26 +784,27 @@
END DO
END DO
END DO
!
IF (nvb.EQ.0) RETURN
!
IF ( nvb == 0 ) RETURN
ALLOCATE( v( nnr ) )
ALLOCATE( qv( nnrb ) )
ALLOCATE( dqgbt( ngb, 2 ) )
ci=(0.,1.)
!
IF(nspin.EQ.1) THEN
! ------------------------------------------------------------------
! nspin=1 : two fft at a time, one per atom, if possible
! ------------------------------------------------------------------
ci =( 0.0d0, 1.0d0 )
IF( nspin == 1 ) THEN
!
! nspin=1 : two fft at a time, one per atom, if possible
!
DO i=1,3
DO j=1,3
!
v(:) = (0.d0, 0.d0)
!
iss=1
isa=1
DO is=1,nvb
#ifdef __PARA
DO ia=1,na(is)
@ -827,9 +816,9 @@
#endif
dqgbt(:,:) = (0.d0, 0.d0)
IF (ia.EQ.na(is)) nfft=1
!
! nfft=2 if two ffts at the same time are performed
!
!
! nfft=2 if two ffts at the same time are performed
!
DO ifft=1,nfft
DO iv=1,nh(is)
DO jv=iv,nh(is)
@ -848,9 +837,9 @@
END DO
END DO
END DO
!
! add structure factor
!
!
! add structure factor
!
qv(:) = (0.d0, 0.d0)
IF(nfft.EQ.2) THEN
DO ig=1,ngb
@ -869,36 +858,37 @@
ENDIF
!
CALL invfft('Box',qv,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx,isa)
!
! qv = US contribution in real space on box grid
! for atomic species is, real(qv)=atom ia, imag(qv)=atom ia+1
!
! add qv(r) to v(r), in real space on the dense grid
!
CALL box2grid(irb(1,isa),1,qv,v)
!
! qv = US contribution in real space on box grid
! for atomic species is, real(qv)=atom ia, imag(qv)=atom ia+1
!
! add qv(r) to v(r), in real space on the dense grid
!
CALL box2grid( irb(1,isa), 1, qv, v )
IF (nfft.EQ.2) CALL box2grid(irb(1,isa+1),2,qv,v)
15 isa=isa+nfft
15 isa = isa + nfft
!
END DO
END DO
!
DO ir=1,nnr
drhor(ir,iss,i,j)=drhor(ir,iss,i,j)+DBLE(v(ir))
drhor(ir,iss,i,j) = drhor(ir,iss,i,j) + DBLE(v(ir))
END DO
!
CALL fwfft('Dense', v,nr1,nr2,nr3,nr1x,nr2x,nr3x)
CALL fwfft( 'Dense', v, nr1, nr2, nr3, nr1x, nr2x, nr3x )
!
DO ig=1,ng
drhog(ig,iss,i,j)=drhog(ig,iss,i,j)+v(np(ig))
drhog(ig,iss,i,j) = drhog(ig,iss,i,j) + v(np(ig))
END DO
!
ENDDO
ENDDO
!
ELSE
! ------------------------------------------------------------------
! nspin=2: two fft at a time, one for spin up and one for spin down
! ------------------------------------------------------------------
!
! nspin=2: two fft at a time, one for spin up and one for spin down
!
isup=1
isdw=2
DO i=1,3
@ -929,9 +919,9 @@
END DO
END DO
END DO
!
! add structure factor
!
!
! add structure factor
!
qv(:) = (0.d0, 0.d0)
DO ig=1,ngb
qv(npb(ig))= eigrb(ig,isa)*dqgbt(ig,1) &
@ -941,14 +931,16 @@
END DO
!
CALL invfft('Box',qv,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx,isa)
!
! qv is the now the US augmentation charge for atomic species is
! and atom ia: real(qv)=spin up, imag(qv)=spin down
!
! add qv(r) to v(r), in real space on the dense grid
!
!
! qv is the now the US augmentation charge for atomic species is
! and atom ia: real(qv)=spin up, imag(qv)=spin down
!
! add qv(r) to v(r), in real space on the dense grid
!
CALL box2grid2(irb(1,isa),qv,v)
25 isa=isa+1
!
25 isa = isa + 1
!
END DO
END DO
!
@ -2674,7 +2666,7 @@
COMPLEX(DP), ALLOCATABLE :: v(:), vs(:)
REAL(DP), ALLOCATABLE :: gagb(:,:)
!
REAL(DP) :: fion1( 3, nat )
REAL(DP), ALLOCATABLE :: fion1( :, : )
REAL(DP), ALLOCATABLE :: stmp( :, : )
!
COMPLEX(DP), ALLOCATABLE :: self_vloc(:)
@ -2682,6 +2674,7 @@
REAL(DP) :: self_ehtet, fpibg
LOGICAL :: ttsic
REAL(DP) :: detmp( 3, 3 ), desr( 6 ), deps( 6 )
REAL(DP) :: detmp2( 3, 3 )
REAL(DP) :: deht( 6 )
!
INTEGER, DIMENSION(6), PARAMETER :: alpha = (/ 1,2,3,2,3,3 /)
@ -2709,15 +2702,11 @@
!
ttsic = ( ABS( self_interaction ) /= 0 )
!
IF( tpre .AND. ttsic ) &
& CALL errore( ' vofrho ', ' in cplib tpre and ttsic are not possible ', 1 )
IF( ttsic ) ALLOCATE( self_vloc( ng ) )
!
! first routine in which fion is calculated: annihilation
!
fion = 0.d0
fion1 = 0.d0
!
! forces on ions, ionic term in real space
!
@ -2727,7 +2716,7 @@
!
CALL r_to_s( tau0, stmp, na, nsp, ainv )
!
CALL vofesr( 0, esr, desr, fion, stmp, tpre, h )
CALL vofesr( 1, esr, desr, fion, stmp, tpre, h )
!
call mp_sum( fion, intra_image_comm )
!
@ -2755,7 +2744,6 @@
drhotmp( 1:ng, :, : ) = drhotmp( 1:ng, :, : ) + drhog( 1:ng, 2, :, : )
ENDIF
END IF
!
! calculation local potential energy
!
@ -2769,17 +2757,23 @@
epseu = wz * DBLE(SUM(vtemp))
IF (gstart == 2) epseu=epseu-vtemp(1)
CALL mp_sum( epseu, intra_image_comm )
epseu = epseu * omega
!
IF( tpre ) THEN
!
CALL denps_box( drhotmp, sfac, vtemp, dps )
!
CALL pseudo_stress( deps, 0.0d0, gagb, sfac, dvps, rhog, omega )
!
call mp_sum( deps, intra_image_comm )
!
DO k = 1, 6
detmp( alpha(k), beta(k) ) = deps(k)
detmp( beta(k), alpha(k) ) = detmp( alpha(k), beta(k) )
END DO
!
dps = dps + MATMUL( detmp(:,:), TRANSPOSE( ainv(:,:) ) )
!
END IF
@ -2791,9 +2785,7 @@
!
self_ehtet = 0.d0
!
IF( ALLOCATED( self_vloc ) ) THEN
self_vloc = 0.d0
END IF
IF( ttsic ) self_vloc = 0.d0
DO is=1,nsp
DO ig=1,ngs
@ -2842,11 +2834,16 @@
END IF
IF(tpre) DEALLOCATE(drhotmp)
! ===================================================================
! forces on ions, ionic term in reciprocal space
! -------------------------------------------------------------------
IF( tprnfor .OR. tfor .OR. tpre) &
& CALL force_ps(rhotmp,rhog,vtemp,ei1,ei2,ei3,fion1)
!
! forces on ions, ionic term in reciprocal space
!
ALLOCATE( fion1( 3, nat ) )
!
fion1 = 0.d0
!
IF( tprnfor .OR. tfor .OR. tpre) THEN
CALL force_ps(rhotmp,rhog,vtemp,ei1,ei2,ei3,fion1)
END IF
!
! calculation hartree + local pseudo potential
@ -2930,8 +2927,10 @@
fion = fion + fion1
END IF
DEALLOCATE( fion1 )
!
IF( ALLOCATED( self_vloc ) ) DEALLOCATE( self_vloc )
IF( ttsic ) DEALLOCATE( self_vloc )
!
! ===================================================================
! fourier transform of total potential to r-space (dense grid)
@ -3093,42 +3092,3 @@
!
END SUBROUTINE vofrho
!
!----------------------------------------------------------------------
SUBROUTINE checkrho(nnr,nspin,rhor,rmin,rmax,rsum,rnegsum)
!----------------------------------------------------------------------
!
! check \int rho(r)dr and the negative part of rho
!
USE kinds, ONLY: DP
USE mp, ONLY: mp_sum
USE mp_global, ONLY: intra_image_comm
IMPLICIT NONE
INTEGER nnr, nspin
REAL(DP) rhor(nnr,nspin), rmin, rmax, rsum, rnegsum
!
REAL(DP) roe
INTEGER ir, iss
!
rsum =0.0
rnegsum=0.0
rmin =100.
rmax =0.0
DO iss = 1, nspin
DO ir = 1, nnr
roe = rhor(ir,iss)
rsum = rsum + roe
IF ( roe < 0.0 ) rnegsum = rnegsum + roe
rmax = MAX( rmax, roe )
rmin = MIN( rmin, roe )
END DO
END DO
CALL mp_sum( rsum, intra_image_comm )
CALL mp_sum( rnegsum, intra_image_comm )
RETURN
END SUBROUTINE checkrho
!______________________________________________________________________

2
CPV/cpr.f90
View File

@ -13,7 +13,7 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
!
USE kinds, ONLY : DP
USE constants, ONLY : bohr_radius_angs, amu_au
USE control_flags, ONLY : iprint, isave, thdyn, tpre, tbuff, &
USE control_flags, ONLY : iprint, isave, thdyn, tpre, &
iprsta, tfor, tvlocw, &
taurdr, tprnfor, tsdc, lconstrain, lwf, &
lneb, lcoarsegrained, ndr, ndw, nomore, &

56
CPV/cprsub.f90
View File

@ -183,7 +183,7 @@ subroutine nlfh( bec, dbec, lambda )
!
USE kinds, ONLY : DP
use cvan, ONLY : nvb, ish
use uspp, ONLY : nhsa => nkb, qq
use uspp, ONLY : nkb, qq
use uspp_param, ONLY : nh, nhm
use ions_base, ONLY : na
use electrons_base, ONLY : nbspx, nbsp, nudx, nspin, nupdwn, iupdwn
@ -195,11 +195,12 @@ subroutine nlfh( bec, dbec, lambda )
!
implicit none
real(DP), intent(in) :: bec( nhsa, nbsp ), dbec( nhsa, nbsp, 3, 3 )
real(DP), intent(in) :: bec( nkb, nbsp ), dbec( nkb, nbsp, 3, 3 )
real(DP), intent(in) :: lambda( nudx, nudx, nspin )
!
integer :: i, j, ii, jj, inl, iv, jv, ia, is, iss, nss, istart
real(DP) :: fpre(3,3)
INTEGER :: i, j, ii, jj, inl, iv, jv, ia, is, iss, nss, istart
INTEGER :: jnl
REAL(DP) :: fpre(3,3), TT, T1, T2
!
REAL(DP), ALLOCATABLE :: tmpbec(:,:), tmpdh(:,:), temp(:,:)
!
@ -208,6 +209,7 @@ subroutine nlfh( bec, dbec, lambda )
fpre(:,:) = 0.d0
do ii=1,3
do jj=1,3
do is=1,nvb
do ia=1,na(is)
!
@ -241,30 +243,27 @@ subroutine nlfh( bec, dbec, lambda )
if(nh(is).gt.0)then
temp = 0.d0
!
call MXMA &
& (tmpdh,1,nudx,tmpbec,1,nhm,temp,1,nudx,nss,nh(is),nss)
!
! do j=1,nss
! do i=1,nss
! temp(i,j)=temp(i,j)*lambda(i,j,iss)
! end do
! end do
! fpre(ii,jj)=fpre(ii,jj)+2.*SUM(temp(1:nss,1:nss))
!
DO j = 1, nss
DO i = 1, nss
fpre(ii,jj) = fpre(ii,jj) + 2.0d0*temp(i,j)*lambda(i,j,iss)
END DO
END DO
do j=1,nss
do i=1,nss
temp(i,j)=temp(i,j)*lambda(i,j,iss)
end do
end do
fpre(ii,jj)=fpre(ii,jj)+2.0d0*SUM(temp(1:nss,1:nss))
endif
!
end do
!
end do
!
end do
!
end do
!
end do
do i=1,3
@ -274,14 +273,17 @@ subroutine nlfh( bec, dbec, lambda )
enddo
enddo
DEALLOCATE ( tmpbec, tmpdh, temp )
IF( iprsta >= 2 ) THEN
WRITE( stdout,*)
WRITE( stdout,*) "constraints contribution to stress"
WRITE( stdout,5555) ((-fpre(i,j),j=1,3),i=1,3)
fpre = MATMUL( fpre, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0
WRITE( stdout,5555) ((fpre(i,j),j=1,3),i=1,3)
WRITE( stdout,*)
END IF
!
DEALLOCATE ( tmpbec, tmpdh, temp )
5555 FORMAT(1x,f12.5,1x,f12.5,1x,f12.5/ &
& 1x,f12.5,1x,f12.5,1x,f12.5/ &
@ -319,8 +321,7 @@ subroutine nlinit
use core, ONLY : rhocb, nlcc_any, allocate_core
use constants, ONLY : pi, fpi
use ions_base, ONLY : na, nsp
use uspp, ONLY : aainit, beta, qq, dvan, nhtol, nhtolm, indv,&
nhsa => nkb, nhsavb=>nkbus
use uspp, ONLY : aainit, beta, qq, dvan, nhtol, nhtolm, indv
use uspp_param, ONLY : kkbeta, qqq, nqlc, betar, lmaxq, dion,&
nbeta, nbetam, lmaxkb, lll, nhm, nh, tvanp
use atom, ONLY : mesh, r, rab, nlcc, numeric
@ -533,14 +534,18 @@ subroutine dqvan2b(ngy,iv,jv,is,ylm,dylm,dqg)
ivs=indv(iv,is)
jvs=indv(jv,is)
!
if (ivs >= jvs) then
ijvs = ivs*(ivs-1)/2 + jvs
else
ijvs = jvs*(jvs-1)/2 + ivs
end if
!
! ijvs is the packed index for (ivs,jvs)
!
ivl=nhtolm(iv,is)
jvl=nhtolm(jv,is)
!
if (ivl > nlx .OR. jvl > nlx) &
call errore (' qvan2 ', ' wrong dimensions (2)', MAX(ivl,jvl))
!
@ -578,15 +583,16 @@ subroutine dqvan2b(ngy,iv,jv,is,ylm,dylm,dqg)
!
! sig= (-i)^l
!
sig=(0.,-1.)**(l-1)
sig=sig*ap(lp,ivl,jvl)
! WRITE(200,'(3I4,2D14.6)') ijvs,l,is,SUM(qradb(:,ijvs,l,is))
sig = (0.0d0,-1.0d0)**(l-1)
!
sig = sig * ap( lp, ivl, jvl )
!
do ij=1,3
do ii=1,3
do ig=1,ngy
dqg(ig,ii,ij) = dqg(ig,ii,ij) + sig * &
& ( ylm(ig,lp) * dqrad(ig,ijvs,l,is,ii,ij) - &
& dylm(ig,lp,ii,ij)*qradb(ig,ijvs,l,is) ) ! SEGNO
& dylm(ig,lp,ii,ij) * qradb(ig,ijvs,l,is) ) ! SEGNO
end do
end do
end do

105
CPV/exch_corr.f90
View File

@ -181,29 +181,29 @@
IF ( ANY( tnlcc ) ) THEN
DO ig = gstart, ngm
tex1 = (0.0_DP , 0.0_DP)
tex1 = ( 0.0d0 , 0.0d0 )
DO is=1,nsp
IF ( tnlcc(is) ) THEN
tex1 = tex1 + sfac( ig, is ) * CMPLX(rhocp(ig,is), 0.d0)
END IF
END DO
tex2 = 0.0_DP
tex2 = 0.0d0
DO ispin = 1, nspin
tex2 = tex2 + CONJG( vxc(ig, ispin) )
END DO
tex3 = DBLE(tex1 * tex2) / SQRT( g( ig ) ) / tpiba
dcc = dcc + tex3 * gagb(:,ig)
END DO
dcc = dcc * 2.0_DP * omega
dcc = dcc * 2.0d0 * omega
! DEBUG
! DO k=1,6
! detmp(alpha(k),beta(k)) = dcc(k)
! detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k))
! END DO
! detmp = MATMUL( detmp(:,:), box%m1(:,:) )
! WRITE( stdout,*) "derivative of e(xc) - nlcc part"
! WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
!DO k=1,6
! detmp(alpha(k),beta(k)) = dcc(k)
! detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k))
!END DO
!detmp = MATMUL( detmp(:,:), box%m1(:,:) )
!WRITE( stdout,*) "derivative of e(xc) - nlcc part"
!WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
END IF
@ -391,10 +391,13 @@
self_exc = 0.d0
!
if( dft_is_meta() ) then
!
call tpssmeta( nnr, nspin, gradr, rhor, kedtaur, exc )
!
else
!
CALL exch_corr_cp(nnr, nspin, gradr, rhor, exc)
!
!
IF ( ttsic ) THEN
CALL exch_corr_cp(nnr, nspin, self_gradr, self_rho, self_exc)
self_exc = sic_alpha * (exc - self_exc)
@ -415,14 +418,8 @@
!
dxc = 0.0d0
!
if (tpre) then
if ( tpre ) then
!
IF( nlcc_any ) CALL denlcc( nnr, nspin, rhor, sfac, drhocg, dcc )
!
! DEBUG
!
! write (stdout,*) "derivative of e(xc) - nlcc part"
! write (stdout,5555) ((dcc(i,j),j=1,3),i=1,3)
!
do iss = 1, nspin
do j=1,3
@ -432,26 +429,15 @@
end do
end do
end do
drc = 0.0d0
IF( nlcc_any ) THEN
do j=1,3
do i=1,3
do ir=1,nnr
drc(i,j) = drc(i,j) + rhor( ir, iss ) * rhoc( ir ) * ainv(j,i)
end do
end do
end do
END IF
dxc = dxc - drc * 1.0d0 / nspin
end do
!
dxc = dxc * omega / ( nr1*nr2*nr3 )
!
call mp_sum ( dxc, intra_image_comm )
!
do j=1,3
do i=1,3
dxc(i,j) = dxc(i,j) + exc * ainv(j,i)
do j = 1, 3
do i = 1, 3
dxc( i, j ) = dxc( i, j ) + exc * ainv( j, i )
end do
end do
!
@ -460,21 +446,19 @@
! write (stdout,*) "derivative of e(xc)"
! write (stdout,5555) ((dxc(i,j),j=1,3),i=1,3)
!
dxc = dxc + dcc
IF( iprsta >= 2 ) THEN
DO i=1,3
DO j=1,3
detmp(i,j)=exc*ainv(j,i)
END DO
END DO
WRITE( stdout,*) "derivative of e(xc) - diag - kbar"
detmp = -1.0d0 * MATMUL( detmp, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0
WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
END IF
!
end if
!
IF( iprsta >= 2 ) THEN
DO i=1,3
DO j=1,3
detmp(i,j)=exc*ainv(j,i)
END DO
END DO
WRITE( stdout,*) "derivative of e(xc) - diag - kbar"
detmp = -1.0d0 * MATMUL( detmp, TRANSPOSE( h ) ) / omega * au_gpa * 10.0d0
WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
END IF
!
!
if (dft_is_gradient()) then
!
@ -496,11 +480,10 @@
IF( ttsic ) THEN
!
!
IF (dft_is_gradient()) then
call gradh( nspin, self_gradr, self_rhog, self_rho, dexc)
!
gradr(:,:, 1) = (1.d0 - sic_alpha ) * gradr(:,:, 1)
gradr(:,:, 2) = (1.d0 - sic_alpha ) * gradr(:,:, 2) + &
& sic_alpha * ( self_gradr(:,:,1) + self_gradr(:,:,2) )
@ -516,6 +499,36 @@
!
ENDIF
IF( tpre ) THEN
!
dcc = 0.0d0
!
IF( nlcc_any ) CALL denlcc( nnr, nspin, rhor, sfac, drhocg, dcc )
!
! DEBUG
!
! write (stdout,*) "derivative of e(xc) - nlcc part"
! write (stdout,5555) ((dcc(i,j),j=1,3),i=1,3)
!
dxc = dxc + dcc
!
do iss = 1, nspin
drc = 0.0d0
IF( nlcc_any ) THEN
do j=1,3
do i=1,3
do ir=1,nnr
drc(i,j) = drc(i,j) + rhor( ir, iss ) * rhoc( ir ) * ainv(j,i)
end do
end do
end do
call mp_sum ( drc, intra_image_comm )
END IF
dxc = dxc - drc * ( 1.0d0 / nspin ) * omega / ( nr1*nr2*nr3 )
end do
!
END IF
IF( ALLOCATED( gradr ) ) DEALLOCATE( gradr )
5555 format(1x,f12.5,1x,f12.5,1x,f12.5/ &

2
CPV/inner_loop.f90
View File

@ -20,7 +20,7 @@
! declares modules
USE kinds, ONLY: dp
USE control_flags, ONLY: iprint, thdyn, tpre, tbuff, iprsta, &
USE control_flags, ONLY: iprint, thdyn, tpre, iprsta, &
tfor, tvlocw, taurdr, &
tprnfor, ndr, ndw, nbeg, nomore, &
tsde, tortho, tnosee, tnosep, trane, &

6
CPV/move_electrons.f90
View File

@ -170,12 +170,6 @@ SUBROUTINE move_electrons( nfi, tfirst, tlast, b1, b2, b3, fion, &
lambdam = lambda
END IF
!
IF( force_pairing ) THEN
lambda( 1:nupdwn(2), 1:nupdwn(2), 2 ) = lambda(1:nupdwn(2), 1:nupdwn(2), 1 )
lambdam( 1:nupdwn(2), 1:nupdwn(2), 2 ) = lambdam(1:nupdwn(2), 1:nupdwn(2), 1 )
lambdap( 1:nupdwn(2), 1:nupdwn(2), 2 ) = lambdap(1:nupdwn(2), 1:nupdwn(2), 1 )
ENDIF
!
! ... calphi calculates phi
! ... the electron mass rises with g**2
!

17
CPV/nl_base.f90
View File

@ -494,9 +494,9 @@
!
!
call start_clock( 'calbec' )
call nlsm1(n,nspmn,nspmx,eigr,c,bec)
call nlsm1( n, nspmn, nspmx, eigr, c, bec )
!
if (iprsta.gt.2) then
if ( iprsta > 2 ) then
WRITE( stdout,*)
do is=1,nspmx
if(nspmx.gt.1) then
@ -590,6 +590,8 @@ SUBROUTINE caldbec( ngw, nkb, n, nspmn, nspmx, eigr, c, dbec, tred )
ixi = 1
signre = -1.0
signim = 1.0
else
CALL errore(' caldbec ', ' l not implemented ', ABS( l ) )
endif
!
do ia=1,na(is)
@ -636,13 +638,13 @@ subroutine dennl( bec, denl )
USE kinds, ONLY : DP
use cvan, only : ish
use uspp_param, only : nh
use uspp, only : nkb, dvan
use uspp, only : nkb, dvan, deeq
use cdvan, ONLY : drhovan, dbec
use ions_base, only : nsp, na
use cell_base, only : h
use io_global, only : stdout
!
use electrons_base, only : n => nbsp, ispin, f, nspin
use electrons_base, only : n => nbsp, ispin, f, nspin, ispin
use reciprocal_vectors, only : gstart
implicit none
@ -677,17 +679,18 @@ subroutine dennl( bec, denl )
enddo
enddo
end do
dsum=0.d0
!
do iss=1,nspin
dsum=0.d0
do k=1,3
do j=1,3
drhovan(ijv,isa,iss,j,k)=dsums(iss,j,k)
dsum(j,k)=dsum(j,k)+dsums(iss,j,k)
enddo
enddo
if(iv.ne.jv) dsum=2.d0*dsum
denl = denl + dsum * dvan(jv,iv,is)
end do
if(iv.ne.jv) dsum=2.d0*dsum
denl = denl + dsum*dvan(jv,iv,is)
end do
end do
end do

10
CPV/nlcc.f90
View File

@ -464,10 +464,10 @@
end do
end if
!
call invfft('Box',qv,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx,isa)
!
! note that a factor 1/2 is hidden in fac if nspin=2
!
call invfft('Box',qv,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx,ia+isa)
!
! note that a factor 1/2 is hidden in fac if nspin=2
!
do iss=1,nspin
fcc(ix,ia+isa) = fcc(ix,ia+isa) + fac * &
& boxdotgrid(irb(1,ia +isa),1,qv,vxc(1,iss))
@ -495,7 +495,7 @@
!
!-----------------------------------------------------------------------
subroutine set_cc(irb,eigrb,rhoc)
subroutine set_cc( irb, eigrb, rhoc )
!-----------------------------------------------------------------------
!
! Calculate core charge contribution in real space, rhoc(r)

8
CPV/ortho_base.f90
View File

@ -935,8 +935,8 @@ CONTAINS
!
IF ( iprsta > 2 ) THEN
WRITE( stdout,*)
DO is=1,nsp
IF(nsp.GT.1) THEN
DO is = 1, nvb
IF( nvb .GT. 1 ) THEN
WRITE( stdout,'(33x,a,i4)') ' updatc: bec (is)',is
WRITE( stdout,'(8f9.4)') &
& ((bec(ish(is)+(iv-1)*na(is)+1,i+istart-1),iv=1,nh(is)),i=1,nss)
@ -1070,8 +1070,8 @@ CONTAINS
WRITE( stdout,*) 'in calphi sqrt(emtot)=',SQRT(emtot)
WRITE( stdout,*)
DO is=1,nsp
IF(nsp > 1) THEN
DO is = 1, nvb
IF( nvb > 1 ) THEN
WRITE( stdout,'(33x,a,i4)') ' calphi: bec (is)',is
WRITE( stdout,'(8f9.4)') &
& ((bec(ish(is)+(iv-1)*na(is)+1,i),iv=1,nh(is)),i=1,n)

56
CPV/pseudopot.f90
View File

@ -727,8 +727,9 @@ CONTAINS
!
IF( ALLOCATED( betagx ) ) DEALLOCATE( betagx )
IF( ALLOCATED( dbetagx ) ) DEALLOCATE( dbetagx )
!
ALLOCATE( betagx ( mmx, nhm, nsp ) )
IF (tpre) ALLOCATE( dbetagx( mmx, nhm, nsp ) )
IF ( tpre ) ALLOCATE( dbetagx( mmx, nhm, nsp ) )
!
do is = 1, nsp
@ -738,6 +739,7 @@ CONTAINS
allocate( djl ( kkbeta( is ) ) )
allocate( jltmp( kkbeta( is ) ) )
end if
!
allocate( fint ( kkbeta( is ) ) )
allocate( jl ( kkbeta( is ) ) )
!
@ -758,21 +760,29 @@ CONTAINS
!
i0 = 1
if ( r(1,is) < 1.0d-8 ) i0 = 2
! special case q=0
!
if ( xg < 1.0d-8 ) then
!
! special case q=0
!
if (l == 1) then
! Note that dj_1/dx (x=0) = 1/3
jltmp(:) = 1.0d0/3.d0
else
jltmp(:) = 0.0d0
end if
!
else
!
call sph_bes (kkbeta(is)+1-i0, r(i0,is), xg, ltmp-1, jltmp )
!
end if
!
do ir = i0, kkbeta(is)
xrg = r(ir,is) * xg
djl(ir) = jltmp(ir) * xrg - l * jl(ir)
end do
!
if ( i0 == 2 ) djl(1) = djl(2)
!
endif
@ -804,6 +814,7 @@ CONTAINS
!
deallocate(jl)
deallocate(fint)
!
if (tpre) then
deallocate(jltmp)
deallocate(djl)
@ -874,7 +885,6 @@ CONTAINS
ALLOCATE( qrl( nr, nbeta(is)*(nbeta(is)+1)/2, nqlc(is)) )
!
call fill_qrl ( is, qrl )
! qrl = 0.0d0
!
do l = 1, nqlc( is )
!
@ -1152,7 +1162,6 @@ CONTAINS
do ig=1,ngb
do l=1,nqlc(is)
qradb(ig,ijv,l,is)= c*qradx(ig,ijv,l,is)
WRITE(100,'(4I4,2D14.6)') ig,ijv,l,is,qradb(ig,ijv,l,is)
enddo
enddo
enddo
@ -1354,11 +1363,6 @@ CONTAINS
do j=1,3
dbeta( 1, iv, is, i, j ) = -0.5 * beta( 1, iv, is ) * ainv( j, i ) &
& - c * dylm( 1, lp, i, j ) * betagl ! SEGNO
!dbeta(1,iv,is,i,j)=-0.5*beta(1,iv,is)*ainv(j,i) !DEBUG stress ok!
!dbeta(1,iv,is,i,j)=-c*dylm(1,lp,i,j)*betagl !DEBUG stress ok!
!dbeta(1,iv,is,i,j)=0.0d0 !DEBUG stress
enddo
enddo
do ig = gstart, ngw
@ -1375,11 +1379,6 @@ CONTAINS
& - c * dylm( ig, lp, i, j ) * betagl & ! SEGNO
& - c * ylm ( ig, lp ) * dbetagl * gx( i, ig ) / g( ig ) &
& * ( gx( 1, ig ) * ainv( j, 1 ) + gx( 2, ig ) * ainv( j, 2 ) + gx( 3, ig ) * ainv( j, 3 ) )
!dbeta(ig,iv,is,i,j)=- 0.5d0 * beta( ig, iv, is ) * ainv( j, i ) !DEBUG stress ok!
!dbeta(ig,iv,is,i,j)=- c * dylm( ig, lp, i, j ) * betagl !DEBUG stress ok!
!dbeta( ig, iv, is, i, j ) = &
! - c * ylm ( ig, lp ) * dbetagl * gx( i, ig ) / g( ig ) &
! * ( gx( 1, ig ) * ainv( j, 1 ) + gx( 2, ig ) * ainv( j, 2 ) + gx( 3, ig ) * ainv( j, 3 ) )
end do
end do
end do
@ -1451,7 +1450,10 @@ CONTAINS
do iv= 1,nbeta(is)
do jv = iv, nbeta(is)
ijv = jv*(jv-1)/2 + iv
do ig=1,ngb
do l=1,nqlc(is)
qradb(1,ijv,l,is) = c * qradx(1,ijv,l,is)
end do
do ig=2,ngb
gg=gb(ig)*tpibab*tpibab/refg
jj=int(gg)+1
do l=1,nqlc(is)
@ -1495,7 +1497,7 @@ CONTAINS
allocate(dqgbs(ngb,3,3))
dqrad(:,:,:,:,:,:) = 0.d0
!
call dylmr2_(lmaxq*lmaxq, ngb, gxb, gb, ainv, dylmb)
call dylmr2_( lmaxq*lmaxq, ngb, gxb, gb, ainv, dylmb )
!
do is=1,nvb
!
@ -1503,7 +1505,8 @@ CONTAINS
do jv=iv,nbeta(is)
ijv = jv*(jv-1)/2 + iv
do l=1,nqlc(is)
do ig=1,ngb
dqradb(1,ijv,l,is) = dqradx(1,ijv,l,is)
do ig=2,ngb
gg=gb(ig)*tpibab*tpibab/refg
jj=int(gg)+1
if(jj.ge.mmx) then
@ -1516,18 +1519,15 @@ CONTAINS
enddo
do i=1,3
do j=1,3
dqrad(1,ijv,l,is,i,j) = &
-qradb(1,ijv,l,is) * ainv(j,i)
! dqrad(1,ijv,l,is,i,j) = 0.d0 ! DEBUG
dqrad(1,ijv,l,is,i,j) = - qradb(1,ijv,l,is) * ainv(j,i)
do ig=2,ngb
dqrad(ig,ijv,l,is,i,j) = &
& -qradb(ig,ijv,l,is)*ainv(j,i) &
& -c*dqradb(ig,ijv,l,is)* &
& - qradb(ig,ijv,l,is)*ainv(j,i) &
& - c * dqradb(ig,ijv,l,is)* &
& gxb(i,ig)/gb(ig)* &
& (gxb(1,ig)*ainv(j,1)+ &
& gxb(2,ig)*ainv(j,2)+ &
& gxb(3,ig)*ainv(j,3))
! dqrad(ig,ijv,l,is,i,j) = -qradb(ig,ijv,l,is)*ainv(j,i) ! DEBUG
enddo
enddo
enddo
@ -1611,6 +1611,7 @@ CONTAINS
allocate( djl ( kkbeta( is ) ) )
allocate( jltmp( kkbeta( is ) ) )
end if
!
allocate( fint ( kkbeta( is ) ) )
allocate( jl ( kkbeta( is ) ) )
!
@ -1677,6 +1678,7 @@ CONTAINS
!
deallocate(jl)
deallocate(fint)
!
if (tpre) then
deallocate(jltmp)
deallocate(djl)
@ -1779,6 +1781,8 @@ CONTAINS
IF ( kkbeta(is) > dim1 ) &
CALL errore ('fill_qrl', 'bad 1st dimension for array qrl', 1)
!
qrl = 0.0d0
!
do iv = 1, nbeta(is)
!
do jv = iv, nbeta(is)
@ -1793,8 +1797,10 @@ CONTAINS
lmin = ABS (lll(jv,is) - lll(iv,is)) + 1
lmax = lll(jv,is) + lll(iv,is) + 1
! WRITE( stdout, * ) ' is, jv, iv = ', is, jv, iv
! WRITE( stdout, * ) ' lll jv, iv = ', lll(jv,is), lll(iv,is)
! WRITE( stdout, * ) 'QRL is, jv, iv = ', is, jv, iv
! WRITE( stdout, * ) 'QRL lll jv, iv = ', lll(jv,is), lll(iv,is)
! WRITE( stdout, * ) 'QRL lmin, lmax = ', lmin, lmax
! WRITE( stdout, * ) '---------------- '
IF ( lmin < 1 .OR. lmax > dim3) THEN
WRITE( stdout, * ) ' lmin, lmax = ', lmin, lmax

184
CPV/runcp.f90
View File

@ -17,6 +17,7 @@
PUBLIC :: runcp, runcp_force_pairing
PUBLIC :: runcp_uspp, runcp_uspp_force_pairing, runcp_ncpp
PUBLIC :: runcp_uspp_bgl
!=----------------------------------------------------------------------------=!
CONTAINS
@ -593,7 +594,7 @@
fromscra, restart )
!
use wave_base, only: wave_steepest, wave_verlet
use control_flags, only: tbuff, lwf, tsde
use control_flags, only: lwf, tsde
!use uspp, only : nhsa=> nkb, betae => vkb, rhovan => becsum, deeq
use uspp, only : deeq, betae => vkb
use reciprocal_vectors, only : gstart
@ -687,20 +688,187 @@
deallocate(c2)
deallocate(c3)
!
!==== end of loop which updates electronic degrees of freedom
!
! buffer for wavefunctions is unit 21
!
if(tbuff) rewind 21
END SUBROUTINE runcp_uspp
!
!
!=----------------------------------------------------------------------------=!
!
!
SUBROUTINE runcp_uspp_bgl( nfi, fccc, ccc, ema0bg, dt2bye, rhos, bec, c0, cm, &
fromscra, restart )
!
use wave_base, only: my_wave_steepest, my_wave_verlet
use control_flags, only: lwf, tsde
use uspp, only: deeq, betae => vkb
use reciprocal_vectors, only: gstart
use electrons_base, only: n=>nbsp, nspin
use wannier_subroutines, only: ef_potential
use efield_module, only: dforce_efield, tefield, dforce_efield2, tefield2
use gvecw, only: ngw
use smooth_grid_dimensions, only: nr1s, nr2s, nr3s, nr1sx, nr2sx, nr3sx, nnrsx
USE fft_base, ONLY: dffts
USE mp_global, ONLY: me_image, nogrp, me_ogrp
USE parallel_include
use task_groups
!
IMPLICIT NONE
integer, intent(in) :: nfi
real(8) :: fccc, ccc
real(8) :: ema0bg(:), dt2bye
real(8) :: rhos(:,:)
real(8) :: bec(:,:)
complex(8) :: c0(:,:), cm(:,:)
logical, optional, intent(in) :: fromscra
logical, optional, intent(in) :: restart
!
real(8) :: verl1, verl2, verl3
real(8), allocatable:: emadt2(:)
real(8), allocatable:: emaver(:)
complex(8), allocatable:: c2(:), c3(:)
integer :: i, index_in, index
integer :: iflag, ierr
logical :: ttsde
iflag = 0
IF( PRESENT( fromscra ) ) THEN
IF( fromscra ) iflag = 1
END IF
IF( PRESENT( restart ) ) THEN
IF( restart ) iflag = 2
END IF
!
! ... set verlet variables
!
verl1 = 2.0d0 * fccc
verl2 = 1.0d0 - verl1
verl3 = 1.0d0 * fccc
allocate(c2(ngw))
allocate(c3(ngw))
ALLOCATE( emadt2( ngw ) )
ALLOCATE( emaver( ngw ) )
ccc = fccc * dt2bye
emadt2 = dt2bye * ema0bg
emaver = emadt2 * verl3
IF( iflag == 0 ) THEN
ttsde = tsde
ELSE IF( iflag == 1 ) THEN
ttsde = .TRUE.
ELSE IF( iflag == 2 ) THEN
ttsde = .FALSE.
END IF
if( lwf ) then
!
call ef_potential( nfi, rhos, bec, deeq, betae, c0, cm, emadt2, emaver, verl1, verl2, c2, c3 )
!
else
IF (.NOT.(ALLOCATED(tg_c2))) ALLOCATE(tg_c2((NOGRP+1)*ngw))
IF (.NOT.(ALLOCATED(tg_c3))) ALLOCATE(tg_c3((NOGRP+1)*ngw))
!---------------------------------------------------------------
!This loop is parallelized accross the eigenstates as well as
!in the FFT, similar to rhoofr
!---------------------------------------------------------------
!------------------------------------
!The potential in rhos
!is distributed accros all processors
!We need to redistribute it so that
!it is completely contained in the
!processors of an orbital TASK-GROUP
!------------------------------------
!
recv_cnt(1) = dffts%npp(NOLIST(1)+1)*nr1sx*nr2sx
recv_displ(1) = 0
DO i = 2, NOGRP
recv_cnt(i) = dffts%npp(NOLIST(i)+1)*nr1sx*nr2sx
recv_displ(i) = recv_displ(i-1) + recv_cnt(i)
ENDDO
IF (.NOT.ALLOCATED(tg_rhos)) ALLOCATE(tg_rhos( (NOGRP+1)*nr1sx*nr2sx*maxval(dffts%npp),nspin))
tg_c3(:) = 0D0
tg_c3(:) = 0D0
tg_rhos(:,:) = 0D0
DO i = 1, nspin
CALL MPI_Allgatherv(rhos(1,i), dffts%npp(me_image+1)*nr1sx*nr2sx, MPI_DOUBLE_PRECISION, &
tg_rhos(1,i), recv_cnt, recv_displ, MPI_DOUBLE_PRECISION, ME_OGRP, IERR)
ENDDO
do i = 1, n, 2*NOGRP ! 2*NOGRP eigenvalues are treated at each iteration
!----------------------------------------------------------------
!The input coefficients to dforce cover eigenstates i:i+2*NOGRP-1
!Thus, in dforce the dummy arguments for c0(1,i,1,1) and
!c0(1,i+1,1,1) hold coefficients for eigenstates i,i+2*NOGRP-2,2
!and i+1,i+2*NOGRP...for example if NOGRP is 4 then we would have
!1-3-5-7 and 2-4-6-8
!----------------------------------------------------------------
call dforce( bec, betae, i, c0(1,i), c0(1,i+1), tg_c2, tg_c3, tg_rhos)
!-------------------------------------------------------
!C. Bekas: This is not implemented yet! I need to see it
!-------------------------------------------------------
if( tefield ) then
CALL errore( ' runcp_uspp ', ' electric field on BGL not implemented yet ', 1 )
end if
IF( iflag == 2 ) THEN
DO index = 1, 2 * NOGRP, 2
cm(:,i+index-1) = c0(:,i+index-1)
cm(:,i+index) = c0(:,i+index)
ENDDO
END IF
index_in = 1
DO index = 1, 2*NOGRP, 2
IF (tsde) THEN
CALL my_wave_steepest( cm(:, i+index-1 ), c0(:, i+index-1 ), emaver, tg_c2, ngw, index_in )
CALL my_wave_steepest( cm(:, i+index), c0(:, i+index), emaver, tg_c3, ngw, index_in )
ELSE
CALL my_wave_verlet( cm(:, i+index-1 ), c0(:, i+index-1 ), &
verl1, verl2, emaver, tg_c2, ngw, index_in )
CALL my_wave_verlet( cm(:, i+index), c0(:, i+index ), &
verl1, verl2, emaver, tg_c3, ngw, index_in )
ENDIF
if ( gstart == 2 ) then
cm(1, i+index-1)=cmplx(real(cm(1, i+index-1)),0.0)
cm(1,i+index)=cmplx(real(cm(1,i+index)),0.0)
end if
index_in = index_in+1
ENDDO ! End loop accross 2*NOGRP current eigenstates
end do ! End loop accross eigenstates
end if
DEALLOCATE( emadt2 )
DEALLOCATE( emaver )
deallocate(c2)
deallocate(c3)
END SUBROUTINE runcp_uspp_bgl
!
!=----------------------------------------------------------------------------=!
!=----------------------------------------------------------------------------=!
SUBROUTINE runcp_uspp_force_pairing( nfi, fccc, ccc, ema0bg, dt2bye, rhos, bec, c0, cm, &
intermed, fromscra, restart )
!
USE wave_base, ONLY: wave_steepest, wave_verlet
USE control_flags, ONLY: tbuff, lwf, tsde
USE control_flags, ONLY: lwf, tsde
! use uspp, only : nhsa=> nkb, betae => vkb, rhovan => becsum, deeq
USE uspp, ONLY : deeq, betae => vkb
USE reciprocal_vectors, ONLY : gstart

16
CPV/smcp.f90
View File

@ -14,7 +14,7 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
! ... main subroutine for SMD by Yosuke Kanai
!
USE kinds, ONLY : DP
USE control_flags, ONLY : iprint, isave, thdyn, tpre, tbuff, &
USE control_flags, ONLY : iprint, isave, thdyn, tpre, &
iprsta, tfor, tvlocw, &
taurdr, tprnfor, tsdc, ndr, ndw, nbeg, &
nomore, tsde, tortho, tnosee, tnosep, &
@ -228,12 +228,6 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
smpm = smd_p -1
!
!
IF(tbuff) THEN
WRITE(stdout,*) "TBUFF set to .FALSE., Option not implemented with SMD"
tbuff = .FALSE.
ENDIF
!
!
! Opening files
!
!
@ -600,10 +594,6 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
CALL runcp_uspp( nfi, fccc(sm_k), ccc(sm_k), ema0bg, dt2bye, rhos, &
rep_el(sm_k)%bec, rep_el(sm_k)%cm, rep_el(sm_k)%c0, fromscra = .TRUE. )
!
! buffer for wavefunctions is unit 21
!
IF(tbuff) REWIND 21
!
! nlfq needs deeq calculated in newd
!
IF ( tfor .OR. tprnfor ) CALL nlfq(rep_el(sm_k)%cm,eigr, &
@ -924,10 +914,6 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
CALL runcp_uspp( nfi, fccc(sm_k), ccc(sm_k), ema0bg, dt2bye, rhos, &
rep_el(sm_k)%bec, rep_el(sm_k)%c0, rep_el(sm_k)%cm )
!
! buffer for wavefunctions is unit 21
!
IF(tbuff) REWIND 21
!
!----------------------------------------------------------------------
! contribution to fion due to lambda
!----------------------------------------------------------------------

42
CPV/stress.f90
View File

@ -560,9 +560,7 @@
! ----------------------------------------------
SUBROUTINE pseudo_stress( deps, epseu, gagb, sfac, dvps, rhoeg, omega )
!
!
USE ions_base, ONLY: nsp
USE reciprocal_vectors, ONLY: gstart
USE gvecs, ONLY: ngs
@ -578,7 +576,7 @@
INTEGER :: ig,k,is, ispin
COMPLEX(DP) :: rhets, depst(6)
!
!
depst = (0.d0,0.d0)
DO is = 1, nsp
@ -606,7 +604,6 @@
RETURN
END SUBROUTINE pseudo_stress
! ----------------------------------------------
@ -752,23 +749,30 @@
SUBROUTINE stress_debug(dekin, deht, dexc, desr, deps, denl, htm1)
USE io_global, ONLY: stdout
USE mp_global, ONLY: intra_image_comm
USE mp, ONLY: mp_sum
REAL(DP) :: dekin(:), deht(:), dexc(:), desr(:), deps(:), denl(:)
REAL(DP) :: detot( 6 ), htm1(3,3)
REAL(DP) :: detmp(3,3)
INTEGER :: k, i, j
detot = dekin + deht + dexc + desr + deps + denl
WRITE( stdout,106) detot
WRITE( stdout,100) dekin
WRITE( stdout,101) deht
WRITE( stdout,102) dexc
WRITE( stdout,103) desr
WRITE( stdout,104) deps
WRITE( stdout,105) denl
DO k=1,6
detmp(alpha(k),beta(k)) = detot(k)
detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k))
END DO
CALL mp_sum( detmp, intra_image_comm )
detmp = MATMUL( detmp(:,:), htm1(:,:) )
WRITE( stdout,*) "derivative of e(tot)"
WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
DO k=1,6
detmp(alpha(k),beta(k)) = dekin(k)
detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k))
END DO
CALL mp_sum( detmp, intra_image_comm )
detmp = MATMUL( detmp(:,:), htm1(:,:) )
WRITE( stdout,*) "derivative of e(kin)"
WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
@ -777,6 +781,7 @@
detmp(alpha(k),beta(k)) = deht(k) + desr(k)
detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k))
END DO
CALL mp_sum( detmp, intra_image_comm )
detmp = MATMUL( detmp(:,:), htm1(:,:) )
WRITE( stdout,*) "derivative of e(electrostatic)"
WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
@ -785,6 +790,7 @@
detmp(alpha(k),beta(k)) = deht(k)
detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k))
END DO
CALL mp_sum( detmp, intra_image_comm )
detmp = MATMUL( detmp(:,:), htm1(:,:) )
WRITE( stdout,*) "derivative of e(h)"
WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
@ -793,6 +799,7 @@
detmp(alpha(k),beta(k)) = desr(k)
detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k))
END DO
CALL mp_sum( detmp, intra_image_comm )
detmp = MATMUL( detmp(:,:), htm1(:,:) )
WRITE( stdout,*) "derivative of e(sr)"
WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
@ -801,6 +808,7 @@
detmp(alpha(k),beta(k)) = deps(k)
detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k))
END DO
CALL mp_sum( detmp, intra_image_comm )
detmp = MATMUL( detmp(:,:), htm1(:,:) )
WRITE( stdout,*) "derivative of e(ps)"
WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
@ -809,6 +817,7 @@
detmp(alpha(k),beta(k)) = denl(k)
detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k))
END DO
CALL mp_sum( detmp, intra_image_comm )
detmp = MATMUL( detmp(:,:), htm1(:,:) )
WRITE( stdout,*) "derivative of e(nl)"
WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
@ -817,20 +826,15 @@
detmp(alpha(k),beta(k)) = dexc(k)
detmp(beta(k),alpha(k)) = detmp(alpha(k),beta(k))
END DO
CALL mp_sum( detmp, intra_image_comm )
detmp = MATMUL( detmp(:,:), htm1(:,:) )
WRITE( stdout,*) "derivative of e(xc)"
WRITE( stdout,5555) ((detmp(i,j),j=1,3),i=1,3)
5555 format(1x,f12.5,1x,f12.5,1x,f12.5/ &
& 1x,f12.5,1x,f12.5,1x,f12.5/ &
& 1x,f12.5,1x,f12.5,1x,f12.5//)
100 FORMAT(' dekin :',6F12.4)
101 FORMAT(' deht :',6F12.4)
102 FORMAT(' dexc :',6F12.4)
103 FORMAT(' desr :',6F12.4)
104 FORMAT(' deps :',6F12.4)
105 FORMAT(' denl :',6F12.4)
106 FORMAT(' detot :',6F12.4)
RETURN
END SUBROUTINE stress_debug

45
CPV/wf.f90
View File

@ -3572,7 +3572,7 @@ SUBROUTINE dforce_field( bec, deeq, betae, i, c, ca, df, da, v, v1 )
! ... e^-ig.r_i < beta_i,j | c_n > }
!
USE kinds, ONLY : DP
USE control_flags, ONLY : iprint, tbuff
USE control_flags, ONLY : iprint
USE gvecs, ONLY : nms, nps
USE gvecw, ONLY : ngw
USE cvan, ONLY : ish
@ -3612,25 +3612,15 @@ SUBROUTINE dforce_field( bec, deeq, betae, i, c, ca, df, da, v, v1 )
END DO
ENDIF
!
IF (.NOT.tbuff) THEN
!
psi( 1:nnrsx ) = 0.D0
DO ig=1,ngw
psi(nms(ig))=CONJG(c(ig)-CI*ca(ig))
psi(nps(ig))=c(ig)+CI*ca(ig)
END DO
!
CALL invfft('Wave',psi,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
!
ELSE
!
! read psi from buffer 21
!
READ(21,iostat=ios) psi
IF(ios.NE.0) CALL errore( 'dforce', 'error in reading unit 21', ios )
!
ENDIF
!
!
psi( 1:nnrsx ) = 0.D0
DO ig=1,ngw
psi(nms(ig))=CONJG(c(ig)-CI*ca(ig))
psi(nps(ig))=c(ig)+CI*ca(ig)
END DO
!
CALL invfft('Wave',psi,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
!
iss1=ispin(i)
!
! the following avoids a potential out-of-bounds error
@ -3866,7 +3856,7 @@ SUBROUTINE rhoiofr( nfi, c, irb, eigrb, bec, &
!
USE constants, ONLY : bohr_radius_angs
USE kinds, ONLY : DP
USE control_flags, ONLY : iprint, tbuff, iprsta, thdyn, tpre, trhor
USE control_flags, ONLY : iprint, iprsta, thdyn, tpre, trhor
USE ions_base, ONLY : nax, nat, nsp, na
USE cell_base, ONLY : a1, a2, a3
USE recvecs_indexes, ONLY : np, nM
@ -3889,6 +3879,7 @@ SUBROUTINE rhoiofr( nfi, c, irb, eigrb, bec, &
USE uspp_param, ONLY : nh, nhm
USE uspp, ONLY : nkb
USE fft_module, ONLY : fwfft, invfft
USE charge_density, ONLY : checkrho
USE mp, ONLY : mp_sum
!
IMPLICIT NONE
@ -4015,9 +4006,6 @@ SUBROUTINE rhoiofr( nfi, c, irb, eigrb, bec, &
!
CALL invfft('Wave',psis,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
!
! wavefunctions in unit 21
!
IF(tbuff) WRITE(21,iostat=ios) psis
! iss1=ispin(i)
iss1=1
sa1=f(i)/omega
@ -4033,17 +4021,8 @@ SUBROUTINE rhoiofr( nfi, c, irb, eigrb, bec, &
rhos(ir,iss2)=rhos(ir,iss2) + sa2*(AIMAG(psis(ir)))**2
END DO
!
! buffer 21
!
IF(tbuff) THEN
IF(ios.NE.0) CALL errore &
& (' rhoofr',' error in writing unit 21',ios)
ENDIF
!
! end do
!
IF(tbuff) REWIND 21
!
! smooth charge in g-space is put into rhog(ig)
!
IF(nspin.EQ.1)THEN

3
Modules/control_flags.f90
View File

@ -53,7 +53,7 @@ MODULE control_flags
tscreen, gamma_only, force_pairing, tchi2
!
PUBLIC :: fix_dependencies, check_flags
PUBLIC :: tbuff, tvlocw, trhor, thdyn, iprsta
PUBLIC :: tvlocw, trhor, thdyn, iprsta
PUBLIC :: twfcollect, printwfc
PUBLIC :: tuspp
PUBLIC :: program_name
@ -66,7 +66,6 @@ MODULE control_flags
! 'CPVC' cp
! 'SMCP' smcp
!
LOGICAL :: tbuff = .FALSE. ! save wfc on unit 21 (only cp, never used)
LOGICAL :: tvlocw = .FALSE. ! write potential to unit 46 (only cp, seldom used)
LOGICAL :: trhor = .FALSE. ! read rho from unit 47 (only cp, seldom used)
!

5
Modules/read_upf.f90
View File

@ -418,6 +418,11 @@ subroutine read_pseudo_nl (upf, iunps)
call scan_begin (iunps, "QFCOEF", .false.)
read (iunps,*,err=106,end=106) &
( ( upf%qfcoef(i,lp,nb,mb), i=1,upf%nqf ), lp=1,upf%nqlc )
do i = 1, upf%nqf
do lp = 1, upf%nqlc
upf%qfcoef(i,lp,mb,nb) = upf%qfcoef(i,lp,nb,mb)
end do
end do
call scan_end (iunps, "QFCOEF")
end if

3
PW/stres_loc.f90
View File

@ -61,7 +61,8 @@ subroutine stres_loc (sigmaloc)
* vloc (igtongl (ng), nt) * fact
enddo
enddo
! WRITE( stdout,*) ' evloc ', evloc, evloc*omega
!
! WRITE( 6,*) ' evloc ', evloc, evloc*omega ! DEBUG
!
do nt = 1, ntyp
! dvloc contains dV_loc(G)/dG