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

- empty states calculation implemented also for CP 

- some more merging
- some dependency removed
- clean-ups


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3190 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
cavazzon 2006-06-22 10:05:15 +00:00
parent c853015e9c
commit c7898d6453
30 changed files with 1477 additions and 1433 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CPV/Makefile
View File

@ -35,7 +35,6 @@ environment.o \
exch_corr.o \
fftdrv.o \
fft.o \
fnl.o \
forces.o \
fromscra.o \
greenf.o \
@ -59,6 +58,7 @@ move_electrons.o \
nl_base.o \
nlcc.o \
nl.o \
optical.o \
ortho_base.o \
ortho.o \
para.o \

4
CPV/cg_sub.f90
View File

@ -311,7 +311,7 @@
call prefor(eigr,betae)!ATTENZIONE
do i=1,n,2
call dforce(bec,betae,i,c0(1,i,1),c0(1,i+1,1),c2,c3,rhos)
call dforce(bec,betae,i,c0(1,i,1),c0(1,i+1,1),c2,c3,rhos,ispin,f,n,nspin)
if(tefield .and. (evalue.ne.0.d0)) then
call dforceb(c0, i, betae, ipolp, bec ,ctabin(1,1,ipolp), gqq, gqqm, qmat, deeq, df)
c2(1:ngw)=c2(1:ngw)+evalue*df(1:ngw)
@ -792,7 +792,7 @@
call prefor(eigr,betae)
do i=1,n,2
call dforce(bec,betae,i,c0(1,i,1),c0(1,i+1,1),c2,c3,rhos)
call dforce(bec,betae,i,c0(1,i,1),c0(1,i+1,1),c2,c3,rhos,ispin,f,n,nspin)
if(tefield.and.(evalue .ne. 0.d0)) then
call dforceb &
(c0, i, betae, ipolp, bec ,ctabin(1,1,ipolp), gqq, gqqm, qmat, deeq, df)

40
CPV/cp_restart.f90
View File

@ -852,6 +852,7 @@ MODULE cp_restart
INTEGER, ALLOCATABLE :: ityp_(:)
INTEGER, ALLOCATABLE :: isrt_(:)
REAL(DP), ALLOCATABLE :: tau_(:,:)
REAL(DP), ALLOCATABLE :: occ_(:)
INTEGER, ALLOCATABLE :: if_pos_(:,:)
CHARACTER(LEN=256) :: psfile_(ntypx)
CHARACTER(LEN=80) :: pos_unit
@ -1152,7 +1153,7 @@ MODULE cp_restart
CALL iotk_scan_dat( iunpun, "NUMBER_OF_SPIN_COMPONENTS", nspin_ )
!
IF ( ( nspin_ /= nspin ) .OR. &
( nbnd_ /= nbnd ) .OR. ( NINT( nelec_ ) /= nelt ) ) THEN
( nbnd_ < nbnd ) .OR. ( NINT( nelec_ ) /= nelt ) ) THEN
!
attr = "electron, bands or spin do not match"
ierr = 30
@ -1184,11 +1185,15 @@ MODULE cp_restart
ik_eff = ik + ( ispin - 1 ) * nk
!
IF ( ionode ) THEN
!
ALLOCATE( occ_ ( nbnd_ ) )
!
CALL iotk_scan_dat( iunpun, &
"OCC" // TRIM( cspin ), occ0(:,ispin), &
"OCC" // TRIM( cspin ), occ_ , &
FOUND = found, IERR = ierr )
!
occ0( 1:nbnd, ispin ) = occ_ ( 1:nbnd )
!
IF( ierr /= 0 ) THEN
attr = "error reading OCC"
GOTO 100
@ -1197,7 +1202,7 @@ MODULE cp_restart
occ0(:,ispin) = occ0(:,ispin) * wk_
!
CALL iotk_scan_dat( iunpun, &
"OCCM" // TRIM( cspin ), occm(:,ispin ), &
"OCCM" // TRIM( cspin ), occ_ , &
FOUND = found, IERR = ierr )
!
IF ( .NOT. found ) THEN
@ -1208,7 +1213,7 @@ MODULE cp_restart
!
ELSE
!
occm(:,ispin) = occm(:,ispin) * wk_
occm( 1:nbnd, ispin ) = occ_ ( 1:nbnd ) * wk_
!
END IF
!
@ -1391,12 +1396,11 @@ MODULE cp_restart
STATUS = 'OLD', FORM = 'UNFORMATTED' )
!
READ( 10 ) lambda0
READ( 10 ) lambdam
!
CLOSE( UNIT = 10 )
!
ELSE
!
lambda0 = 0.0d0
!
END IF
!
END IF
@ -1404,9 +1408,31 @@ MODULE cp_restart
DO ispin = 1, nspin
DO ib = 1, SIZE( lambda0, 2 )
CALL mp_bcast( lambda0( :, ib, ispin), ionode_id, intra_image_comm )
END DO
END DO
!
IF ( ionode ) THEN
!
IF ( found ) THEN
!
READ( 10 ) lambdam
!
CLOSE( UNIT = 10 )
!
ELSE
!
lambdam = 0.0d0
!
END IF
!
END IF
!
DO ispin = 1, nspin
DO ib = 1, SIZE( lambda0, 2 )
CALL mp_bcast( lambdam( :, ib, ispin), ionode_id, intra_image_comm )
END DO
END DO
!
s1 = cclock()
!

221
CPV/cplib.f90
View File

@ -279,7 +279,7 @@
!
!-------------------------------------------------------------------------
SUBROUTINE dforce ( bec, betae, i, c, ca, df, da, v )
SUBROUTINE dforce ( bec, betae, i, c, ca, df, da, v, ispin, f, n, nspin )
!-----------------------------------------------------------------------
!computes: the generalized force df=CMPLX(dfr,dfi) acting on the i-th
! electron state at the gamma point of the brillouin zone
@ -297,7 +297,6 @@
USE uspp_param, ONLY: nhm, nh
USE smooth_grid_dimensions, ONLY: nr1s, nr2s, nr3s, &
nr1sx, nr2sx, nr3sx, nnrsx
USE electrons_base, ONLY: n => nbsp, ispin, f, nspin
USE constants, ONLY: pi, fpi
USE ions_base, ONLY: nsp, na, nat
USE gvecw, ONLY: ggp
@ -308,21 +307,28 @@
!
IMPLICIT NONE
!
COMPLEX(8) betae(ngw,nhsa), c(ngw), ca(ngw), df(ngw), da(ngw)
REAL(8) bec(nhsa,n), v(nnrsx,nspin)
INTEGER i
! local variables
INTEGER, INTENT(IN) :: i
INTEGER, INTENT(IN) :: n, nspin
INTEGER, INTENT(IN) :: ispin( n )
REAL(DP), INTENT(IN) :: f( n )
!
COMPLEX(DP) betae(ngw,nhsa), c(ngw), ca(ngw), df(ngw), da(ngw)
REAL(DP) bec(nhsa,n), v(nnrsx,nspin)
!
! local variables
!
INTEGER iv, jv, ia, is, isa, ism, ios, iss1, iss2, ir, ig, inl, jnl
REAL(8) fi, fip, dd
COMPLEX(8) fp,fm,ci
REAL(8) af(nhsa), aa(nhsa) ! automatic arrays
COMPLEX(8) dtemp(ngw) !
COMPLEX(8), ALLOCATABLE :: psi(:)
REAL(DP) fi, fip, dd
COMPLEX(DP) fp,fm,ci
REAL(DP) af(nhsa), aa(nhsa) ! automatic arrays
COMPLEX(DP), ALLOCATABLE :: dtemp(:)
COMPLEX(DP), ALLOCATABLE :: psi(:)
!
!
CALL start_clock( 'dforce' )
!
ALLOCATE( psi( nnrsx ) )
ALLOCATE( dtemp( ngw ) )
!
! important: if n is odd => c(*,n+1)=0.
!
@ -416,18 +422,14 @@
END DO
END DO
!
DO ig=1,ngw
dtemp(ig)=(0.,0.)
END DO
dtemp = 0.0d0
CALL MXMA &
& (betae,1,2*ngw,af,1,nhsa,dtemp,1,2*ngw,2*ngw,nhsa,1)
DO ig=1,ngw
df(ig)=df(ig)+dtemp(ig)
END DO
!
DO ig=1,ngw
dtemp(ig)=(0.,0.)
END DO
dtemp = 0.0d0
CALL MXMA &
& (betae,1,2*ngw,aa,1,nhsa,dtemp,1,2*ngw,2*ngw,nhsa,1)
DO ig=1,ngw
@ -435,6 +437,7 @@
END DO
ENDIF
DEALLOCATE( dtemp )
DEALLOCATE( psi )
!
CALL stop_clock( 'dforce' )
@ -442,6 +445,7 @@
RETURN
END SUBROUTINE dforce
!
!-----------------------------------------------------------------------
SUBROUTINE dotcsc( eigr, cp, ngw, n )
!-----------------------------------------------------------------------
@ -515,6 +519,68 @@
END SUBROUTINE dotcsc
!-----------------------------------------------------------------------
SUBROUTINE dotcsv( csv, eigr, c, v, ngw )
!-----------------------------------------------------------------------
!
USE kinds, ONLY: DP
USE ions_base, ONLY: na, nsp, nat
USE io_global, ONLY: stdout
USE reciprocal_vectors, ONLY: gstart
USE cvan, ONLY: ish, nvb
USE uspp, ONLY: nkb, qq
USE uspp_param, ONLY: nh
USE mp, ONLY: mp_sum
USE mp_global, ONLY: intra_image_comm
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: ngw
COMPLEX(DP) :: eigr(ngw,nat), c(ngw), v(ngw)
REAL(DP), INTENT(OUT) :: csv
! local variables
COMPLEX(DP) temp(ngw) ! automatic array
REAL(DP), ALLOCATABLE :: bec(:), bev(:)
INTEGER ig,is,ia,iv,jv,inl,jnl
!
ALLOCATE(bec(nkb))
ALLOCATE(bev(nkb))
!
! < beta | c > is real. only the i lowest:
!
CALL nlsm1(1,1,nvb,eigr,c,bec)
CALL nlsm1(1,1,nvb,eigr,v,bev)
!
DO ig=1,ngw
temp(ig)=CONJG(c(ig))*v(ig)
END DO
csv = 2.0d0 * DBLE(SUM(temp))
IF (gstart == 2) csv = csv - DBLE(temp(1))
CALL mp_sum( csv, intra_image_comm )
DO is=1,nvb
DO iv=1,nh(is)
DO jv=1,nh(is)
DO ia=1,na(is)
inl=ish(is)+(iv-1)*na(is)+ia
jnl=ish(is)+(jv-1)*na(is)+ia
csv = csv + qq(iv,jv,is)*bec(inl)*bev(jnl)
END DO
END DO
END DO
END DO
!
DEALLOCATE(bec)
DEALLOCATE(bev)
!
RETURN
END SUBROUTINE dotcsv
!-----------------------------------------------------------------------
SUBROUTINE drhov(irb,eigrb,rhovan,rhog,rhor)
!-----------------------------------------------------------------------
@ -1098,6 +1164,127 @@
RETURN
END SUBROUTINE gracsc
!-------------------------------------------------------------------------
SUBROUTINE smooth_csv( c, v, ngwx, csv, n )
!-----------------------------------------------------------------------
USE gvecw, ONLY: ngw
USE kinds, ONLY: DP
USE reciprocal_vectors, ONLY: gstart
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: ngwx, n
REAL(DP) :: c( 2, ngwx )
REAL(DP) :: v( 2, ngwx, n )
REAL(DP) :: csv( n )
INTEGER :: k, ig
REAL(DP), ALLOCATABLE :: temp(:)
!
! calculate csv(k)=<c|v(k)>
!
ALLOCATE( temp( ngw ) )
DO k = 1, n
DO ig = 1, ngw
temp(ig) = v(1,ig,k) * c(1,ig) + v(2,ig,k) * c(2,ig)
END DO
csv(k) = 2.0d0 * SUM(temp)
IF (gstart == 2) csv(k) = csv(k) - temp(1)
END DO
DEALLOCATE( temp )
!
RETURN
END SUBROUTINE smooth_csv
!-------------------------------------------------------------------------
SUBROUTINE grabec( becc, nkbx, betae, c, ngwx )
!-----------------------------------------------------------------------
!
! on output: bec(i) is recalculated
!
USE uspp, ONLY : nkb, nhsavb=>nkbus
USE gvecw, ONLY: ngw
USE kinds, ONLY: DP
USE reciprocal_vectors, ONLY: gstart
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nkbx, ngwx
COMPLEX(DP) :: betae( ngwx, nkb )
REAL(DP) :: becc( nkbx ), c( 2, ngwx )
INTEGER :: ig, inl
REAL(DP), ALLOCATABLE :: temp(:)
!
ALLOCATE( temp( ngw ) )
!
! calculate becc=<c|beta>
!
DO inl=1,nhsavb
DO ig=1,ngw
temp(ig)=c(1,ig)* DBLE(betae(ig,inl))+ &
& c(2,ig)*AIMAG(betae(ig,inl))
END DO
becc(inl)=2.*SUM(temp)
IF (gstart == 2) becc(inl)= becc(inl)-temp(1)
END DO
DEALLOCATE( temp )
RETURN
END SUBROUTINE grabec
!-------------------------------------------------------------------------
SUBROUTINE bec_csv( becc, becv, nkbx, csv, n )
!-----------------------------------------------------------------------
! requires in input the updated becc and becv(k)
! on output: csv is updated
!
USE ions_base, ONLY: na
USE cvan, ONLY :nvb, ish
USE uspp, ONLY : nkb, nhsavb=>nkbus, qq
USE uspp_param, ONLY: nh
USE kinds, ONLY: DP
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nkbx, n
REAL(DP) :: becc( nkbx )
REAL(DP) :: becv( nkbx, n )
REAL(DP) :: csv( n )
INTEGER :: k, is, iv, jv, ia, inl, jnl
REAL(DP) :: rsum
! calculate csv(k) = csv(k) + <c| SUM_nm |beta(n)><beta(m)|v(k)>, k<i
!
DO k=1,n
rsum=0.
DO is=1,nvb
DO iv=1,nh(is)
DO jv=1,nh(is)
IF(ABS(qq(iv,jv,is)).GT.1.e-5) THEN
DO ia=1,na(is)
inl=ish(is)+(iv-1)*na(is)+ia
jnl=ish(is)+(jv-1)*na(is)+ia
rsum = rsum + qq(iv,jv,is)*becc(inl)*becv(jnl,k)
END DO
ENDIF
END DO
END DO
END DO
csv(k)=csv(k)+rsum
END DO
!
RETURN
END SUBROUTINE bec_csv
!-------------------------------------------------------------------------
SUBROUTINE gram( betae, bec, nkbx, cp, ngwx, n )
!-----------------------------------------------------------------------

14
CPV/cpr.f90
View File

@ -138,6 +138,7 @@ SUBROUTINE cprmain( tau_out, fion_out, etot_out )
USE wave_functions, ONLY : elec_fakekine
USE mp, ONLY : mp_bcast
USE mp_global, ONLY : intra_image_comm
USE empty_states, ONLY : empty_cp
!
IMPLICIT NONE
!
@ -454,8 +455,8 @@ SUBROUTINE cprmain( tau_out, fion_out, etot_out )
!
IF ( tortho ) THEN
!
CALL ortho( eigr, cm(:,:,1), phi(:,:,1), &
lambda, bigr, iter, ccc, bephi, becp )
CALL ortho( eigr, cm(:,:,1), phi(:,:,1), ngw, lambda, SIZE(lambda,1), &
bigr, iter, ccc, bephi, becp, nbsp, nspin, nupdwn, iupdwn )
!
ELSE
!
@ -531,7 +532,7 @@ SUBROUTINE cprmain( tau_out, fion_out, etot_out )
!
IF ( .NOT. tcg ) THEN
!
CALL elec_fakekine( ekinc0, ema0bg, emass, c0, cm, ngw, nbsp, delt )
CALL elec_fakekine( ekinc0, ema0bg, emass, c0, cm, ngw, nbsp, 1, delt )
!
ekinc = ekinc0
!
@ -593,6 +594,9 @@ SUBROUTINE cprmain( tau_out, fion_out, etot_out )
!
CALL cp_eigs( nfi, lambdap, lambda )
!
! ... Compute empty states
!
CALL empty_cp ( nfi, c0, rhos )
!
END IF
!
@ -643,7 +647,7 @@ SUBROUTINE cprmain( tau_out, fion_out, etot_out )
econt = econt + electrons_nose_nrg( xnhe0, vnhe, qne, ekincw )
IF ( tnoseh ) &
econt = econt + cell_nose_nrg( qnh, xnhh0, vnhh, temph, iforceh )
!
!
tps = tps + delt * au_ps
!
CALL printout_new( nfi, tfirst, tfile, ttprint, tps, hold, stress, &
@ -802,7 +806,7 @@ SUBROUTINE cprmain( tau_out, fion_out, etot_out )
END DO main_loop
!
!===================== end of main loop of molecular dynamics ===============
!
!
! ... Here copy relevant physical quantities into the output arrays/variables
!
etot_out = etot

6
CPV/dealloc.f90
View File

@ -55,7 +55,6 @@ SUBROUTINE deallocate_modules_var()
USE ions_nose, ONLY : ions_nose_deallocate
USE metagga, ONLY : deallocate_metagga
USE ncpp, ONLY : deallocate_ncpp
USE pseudo_projector, ONLY : fnl, projector, deallocate_projector
USE task_groups, ONLY : deallocate_groups
USE ions_positions, ONLY : deallocate_ions_positions
USE guess, ONLY : guess_closeup
@ -73,11 +72,6 @@ SUBROUTINE deallocate_modules_var()
IF ( ALLOCATED( dqgb ) ) DEALLOCATE( dqgb )
IF ( ALLOCATED( dbeta ) ) DEALLOCATE( dbeta )
!
IF ( ALLOCATED( fnl ) ) THEN
CALL deallocate_projector( fnl )
DEALLOCATE( fnl )
END IF
!
CALL deallocate_mainvar()
CALL deallocate_ions_positions()
CALL deallocate_cvan()

4
CPV/eigs0.f90
View File

@ -6,7 +6,7 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
subroutine eigs0( tprint, nspin, nupdwn, iupdwn, lf, f, nx, lambda, nudx )
subroutine eigs0( ei, tprint, nspin, nupdwn, iupdwn, lf, f, nx, lambda, nudx )
!-----------------------------------------------------------------------
! computes eigenvalues (wr) of the real symmetric matrix lambda
! Note that lambda as calculated is multiplied by occupation numbers
@ -15,7 +15,6 @@
use kinds, only : DP
use io_global, only : stdout
use constants, only : autoev
use electrons_module, only : ei
use parallel_toolkit, only : dspev_drv
USE sic_module, only : self_interaction
@ -24,6 +23,7 @@
logical, intent(in) :: tprint, lf
integer, intent(in) :: nspin, nx, nudx, nupdwn(nspin), iupdwn(nspin)
real(DP), intent(in) :: lambda( nudx, nudx, nspin ), f( nx )
real(DP), intent(out) :: ei( nudx, nspin )
! local variables
real(DP), allocatable :: lambdar(:), wr(:)
real(DP) zr(1)

31
CPV/electrons.f90
View File

@ -13,7 +13,8 @@
USE parameters, ONLY: nspinx
USE parallel_toolkit, ONLY: pdspev_drv, dspev_drv, pzhpev_drv, zhpev_drv
USE electrons_base, ONLY: nbnd, nbndx, nbsp, nbspx, nspin, nel, nelt, &
nupdwn, iupdwn, telectrons_base_initval, f
nupdwn, iupdwn, telectrons_base_initval, f, &
nudx
USE cp_electronic_mass, ONLY: ecutmass => emass_cutoff
USE cp_electronic_mass, ONLY: emass
USE cp_electronic_mass, ONLY: emass_precond
@ -141,18 +142,22 @@
END IF
DO i = 1, nspin
!
IF( i > nspinx ) CALL errore( ' bmeshset ',' spin too large ', i)
!
nb_l( i ) = nupdwn( i ) / nproc_image
IF( me_image < MOD( nupdwn( i ), nproc_image ) ) nb_l( i ) = nb_l( i ) + 1
n_emp_l( i ) = n_emp / nproc_image
IF( me_image < MOD( n_emp, nproc_image ) ) n_emp_l( i ) = n_emp_l( i ) + 1
!
n_emp_l( i ) = nupdwn_emp( i ) / nproc_image
IF( me_image < MOD( nupdwn_emp( i ), nproc_image ) ) n_emp_l( i ) = n_emp_l( i ) + 1
!
END DO
IF( ALLOCATED( ib_owner ) ) DEALLOCATE( ib_owner )
ALLOCATE( ib_owner( nbndx ), STAT=ierr)
ALLOCATE( ib_owner( MAX( n_emp, nbndx ) ), STAT=ierr)
IF( ierr/=0 ) CALL errore( ' bmeshset ',' allocating ib_owner ', ierr)
IF( ALLOCATED( ib_local ) ) DEALLOCATE( ib_local )
ALLOCATE( ib_local( nbndx ), STAT=ierr)
ALLOCATE( ib_local( MAX( n_emp, nbndx ) ), STAT=ierr)
IF( ierr/=0 ) CALL errore( ' bmeshset ',' allocating ib_local ', ierr)
! here define the association between processors and electronic states
@ -160,7 +165,7 @@
ib_local = 0
ib_owner = -1
DO i = 1, nbndx
DO i = 1, MAX( n_emp, nbndx )
ib_local( i ) = ( i - 1 ) / nproc_image ! local index of the i-th band
ib_owner( i ) = MOD( ( i - 1 ), nproc_image ) ! owner of th i-th band
IF( me_image <= ib_owner( i ) ) THEN
@ -190,6 +195,11 @@
CALL errore( ' electrons_setup ', ' electrons_base not initialized ', 1 )
n_emp = n_emp_
!
! assure that the number of empty states is an even number
!
n_emp = n_emp + MOD( n_emp, 2 )
!
nupdwn_emp(1) = n_emp
iupdwn_emp(1) = 1
@ -198,11 +208,8 @@
iupdwn_emp(2) = 1 + n_emp
END IF
IF( n_emp > nbndx ) &
CALL errore( ' electrons_setup ', ' too many empty states ', 1 )
IF( ALLOCATED( ei ) ) DEALLOCATE( ei )
ALLOCATE( ei( nbnd, nspin ), STAT=ierr)
ALLOCATE( ei( nudx, nspin ), STAT=ierr)
IF( ierr/=0 ) CALL errore( ' electrons ',' allocating ei ',ierr)
ei = 0.0_DP
@ -623,9 +630,9 @@
REAL(DP) :: lambda( :, :, : ), lambdap( :, :, : )
if( .not. tens ) then
call eigs0( .false. , nspin, nupdwn, iupdwn, .true. , f, nx, lambda, nudx )
call eigs0( ei, .false. , nspin, nupdwn, iupdwn, .true. , f, nx, lambda, nudx )
else
call eigs0( .false. , nspin, nupdwn, iupdwn, .false. , f, nx, lambdap, nudx )
call eigs0( ei, .false. , nspin, nupdwn, iupdwn, .false. , f, nx, lambdap, nudx )
endif
WRITE( stdout, * )

872
CPV/emptystates.f90
View File

File diff suppressed because it is too large Load Diff

35
CPV/forces.f90
View File

@ -177,28 +177,28 @@
SUBROUTINE dforce( ib, iss, c, f, df, da, v, eigr, bec, nupdwn, iupdwn )
SUBROUTINE dforce( ib, c, f, df, da, v, eigr, bec, n, noffset )
!
USE reciprocal_vectors, ONLY: ggp, g, gx
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: ib, iss ! band and spin index
INTEGER, INTENT(IN) :: ib ! band and spin index
COMPLEX(DP), INTENT(IN) :: c(:,:)
COMPLEX(DP), INTENT(OUT) :: df(:), da(:)
REAL (DP), INTENT(IN) :: v(:), bec(:,:), f(:)
COMPLEX(DP), INTENT(IN) :: eigr(:,:)
INTEGER, INTENT(IN) :: nupdwn(:), iupdwn(:)
INTEGER, INTENT(IN) :: n, noffset ! number of bands, and band index offset
!
COMPLEX(DP), ALLOCATABLE :: dum( : )
!
INTEGER :: ig, in
INTEGER :: in
!
IF( ib > nupdwn( iss ) ) CALL errore( ' dforce ', ' ib out of range ', 1 )
IF( ib > n ) CALL errore( ' dforce ', ' ib out of range ', 1 )
!
in = iupdwn( iss ) + ib - 1
in = noffset + ib - 1
!
IF( ib == nupdwn( iss ) ) THEN
IF( ib == n ) THEN
!
ALLOCATE( dum( SIZE( da ) ) )
!
@ -223,33 +223,32 @@
! ----------------------------------------------
SUBROUTINE dforce_all( iss, c, f, cgrad, vpot, eigr, bec, nupdwn, iupdwn )
SUBROUTINE dforce_all( c, f, cgrad, vpot, eigr, bec, n, noffset )
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: iss
COMPLEX(DP), INTENT(INOUT) :: c(:,:)
REAL(DP), INTENT(IN) :: vpot(:), f(:)
COMPLEX(DP), INTENT(OUT) :: cgrad(:,:)
COMPLEX(DP), INTENT(IN) :: eigr(:,:)
REAL(DP), INTENT(IN) :: bec(:,:)
INTEGER, INTENT(IN) :: nupdwn(:), iupdwn(:)
INTEGER, INTENT(IN) :: n, noffset
INTEGER :: ib
!
IF( nupdwn( iss ) > 0 ) THEN
IF( n > 0 ) THEN
!
! Process two states at the same time
!
DO ib = 1, nupdwn( iss )-1, 2
CALL dforce( ib, iss, c, f, cgrad(:,ib), cgrad(:,ib+1), &
vpot, eigr, bec, nupdwn, iupdwn )
DO ib = 1, n-1, 2
CALL dforce( ib, c, f, cgrad(:,ib), cgrad(:,ib+1), &
vpot, eigr, bec, n, noffset )
END DO
!
IF( MOD( nupdwn( iss ), 2 ) /= 0 ) THEN
ib = nupdwn( iss )
CALL dforce( ib, iss, c, f, cgrad(:,ib), cgrad(:,ib), &
vpot, eigr, bec, nupdwn, iupdwn )
IF( MOD( n, 2 ) /= 0 ) THEN
ib = n
CALL dforce( ib, c, f, cgrad(:,ib), cgrad(:,ib), &
vpot, eigr, bec, n, noffset )
END IF
!
END IF

12
CPV/fromscra.f90
View File

@ -64,9 +64,9 @@ CONTAINS
USE ions_nose, ONLY : xnhp0, xnhpm, vnhp
USE cp_electronic_mass, ONLY : emass
USE efield_module, ONLY : tefield, efield_berry_setup, &
berry_energy, dforce_efield, &
berry_energy, &
tefield2, efield_berry_setup2, &
berry_energy2, dforce_efield2
berry_energy2
USE cg_module, ONLY : tcg
USE ensemble_dft, ONLY : tens, compute_entropy
USE runcp_module, ONLY : runcp_uspp, runcp_ncpp, runcp_uspp_force_pairing, &
@ -78,11 +78,10 @@ CONTAINS
USE wave_base, ONLY : wave_steepest
USE wavefunctions_module, ONLY : c0, cm, phi => cp
USE wave_types, ONLY : wave_descriptor
USE wave_functions, ONLY : fixwave, wave_rand_init, elec_fakekine
USE wave_functions, ONLY : wave_rand_init, elec_fakekine
USE nl, ONLY : nlrh_m
USE charge_density, ONLY : rhoofr
USE potentials, ONLY : vofrhos
USE forces, ONLY : dforce_all
USE grid_dimensions, ONLY : nr1, nr2, nr3
USE print_out_module, ONLY : printout
@ -322,7 +321,8 @@ CONTAINS
if( tortho ) then
CALL ortho( eigr, c0(:,:,1), phi(:,:,1), lambda, bigr, iter, ccc, bephi, becp )
CALL ortho( eigr, c0(:,:,1), phi(:,:,1), ngw, lambda, SIZE(lambda,1), &
bigr, iter, ccc, bephi, becp, nbsp, nspin, nupdwn, iupdwn )
else
CALL gram( vkb, bec, nkb, c0, ngw, nbsp )
endif
@ -369,7 +369,7 @@ CONTAINS
vnhh (:,:) =0.
velh (:,:)=(h(:,:)-hold(:,:))/delt
!
CALL elec_fakekine( ekincm, ema0bg, emass, c0, cm, ngw, nbsp, delt )
CALL elec_fakekine( ekincm, ema0bg, emass, c0, cm, ngw, nbsp, 1, delt )
xnhe0=0.
xnhem=0.

2
CPV/init.f90
View File

@ -179,7 +179,7 @@
!
IF( program_name == 'FPMD' ) THEN
!
CALL cpsizes( nproc_image )
CALL cpsizes( )
!
END IF
!

2
CPV/init_run.f90
View File

@ -63,7 +63,6 @@ SUBROUTINE init_run()
USE gvecp, ONLY : ecutp
USE funct, ONLY : dft_is_meta
USE metagga, ONLY : crosstaus, dkedtaus, gradwfc
USE pseudo_projector, ONLY : fnl, projector
USE pseudopotential, ONLY : pseudopotential_indexes, nsanl
!
USE efcalc, ONLY : clear_nbeg
@ -84,7 +83,6 @@ SUBROUTINE init_run()
USE printout_base, ONLY : printout_base_init
USE print_out_module, ONLY : print_legend
USE wave_types, ONLY : wave_descriptor_info
USE pseudo_projector, ONLY : allocate_projector
!
IMPLICIT NONE
!

4
CPV/inner_loop.f90
View File

@ -184,7 +184,7 @@
DO i= 1, n, 2
CALL dforce( bec, betae, i, c0(1,i), &
c0(1,i+1), h0c0(1,i), h0c0(1,i+1), &
rhos)
rhos, ispin, f, n, nspin )
END DO
! calculates the Hamiltonian matrix in the basis {c0}
@ -322,7 +322,7 @@
h0c0( :, : )= 0.D0
DO i= 1, n, 2
CALL dforce( bec, betae, i, c0(1,i), c0(1,i+1), &
h0c0(1,i), h0c0(1,i+1), rhos )
h0c0(1,i), h0c0(1,i+1), rhos, ispin, f, n, nspin )
END DO
DO is= 1, nspin
nss= nupdwn( is )

146
CPV/ksstates.f90
View File

@ -27,7 +27,7 @@ MODULE kohn_sham_states
INTEGER, ALLOCATABLE :: indx_ksout_emp(:,:) ! (state inds, spin indxs)
INTEGER, ALLOCATABLE :: n_ksout_emp(:) ! (spin indxs)
PUBLIC :: ks_states_init, kohn_sham, ks_states_closeup
PUBLIC :: ks_states_init, ks_states_closeup
PUBLIC :: n_ksout, indx_ksout, ks_states, tksout
PUBLIC :: ks_states_force_pairing
@ -97,59 +97,6 @@ CONTAINS
! ----------------------------------------------
SUBROUTINE kohn_sham(ispin, c, cdesc, eforces, nupdwn )
! ... declare modules
USE kinds
USE wave_functions, ONLY: crot
USE wave_constrains, ONLY: update_lambda
USE wave_types, ONLY: wave_descriptor
USE electrons_module, ONLY: nb_l
IMPLICIT NONE
! ... declare subroutine arguments
COMPLEX(DP), INTENT(INOUT) :: c(:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
INTEGER, INTENT(IN) :: ispin
INTEGER, INTENT(IN) :: nupdwn(:)
COMPLEX(DP) :: eforces(:,:)
! ... declare other variables
INTEGER :: ib, nb_g, nrl
REAL(DP), ALLOCATABLE :: gam(:,:)
REAL(DP), ALLOCATABLE :: eig(:)
LOGICAL :: tortho = .TRUE.
! ... end of declarations
nb_g = nupdwn( ispin )
IF( nb_g < 1 ) THEN
eforces = 0.0d0
ELSE
nrl = nb_l( ispin )
ALLOCATE( eig( nb_g ) )
ALLOCATE( gam( nrl, nb_g ) )
DO ib = 1, nb_g
CALL update_lambda( ib, gam, c(:,:), cdesc, eforces(:,ib) )
END DO
CALL crot( ispin, c(:,:), cdesc, gam, eig )
DEALLOCATE( gam, eig )
END IF
RETURN
! ...
END SUBROUTINE kohn_sham
! ----------------------------------------------
SUBROUTINE ks_states(cf, wfill, occ, vpot, eigr, bec )
@ -159,11 +106,11 @@ CONTAINS
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE wave_types, ONLY: wave_descriptor, wave_descriptor_init
USE wave_functions, ONLY: kohn_sham
USE forces
USE pseudo_projector, ONLY: projector
USE control_flags, ONLY: force_pairing
USE electrons_base, ONLY: nupdwn, iupdwn, nspin
USE electrons_module, ONLY: n_emp, nupdwn_emp, iupdwn_emp
USE electrons_module, ONLY: n_emp, nupdwn_emp, iupdwn_emp, nb_l, n_emp_l
USE empty_states, ONLY: readempty
IMPLICIT NONE
@ -176,7 +123,7 @@ CONTAINS
REAL (DP) :: vpot(:,:)
! ... declare other variables
INTEGER :: i, ib, ig, ngw, nb_g, nb_l, ispin, iks
INTEGER :: i, ib, ig, ngw, nb_g, ispin, iks
INTEGER :: ispin_wfc
LOGICAL :: tortho = .TRUE.
LOGICAL :: exist
@ -185,7 +132,7 @@ CONTAINS
CHARACTER(LEN=10), DIMENSION(2) :: spin_name
COMPLEX(DP), ALLOCATABLE :: eforce(:,:)
COMPLEX(DP), ALLOCATABLE :: ce(:,:,:)
COMPLEX(DP), ALLOCATABLE :: ce(:,:)
REAL(DP), ALLOCATABLE :: fi(:)
CHARACTER (LEN=6), EXTERNAL :: int_to_char
@ -200,19 +147,17 @@ CONTAINS
DO ispin = 1, nspin
nb_l = wfill%nbl( ispin )
ispin_wfc = ispin
IF( force_pairing ) ispin_wfc = 1
IF( nb_l > 0 ) THEN
IF( nupdwn( ispin ) > 0 ) THEN
ALLOCATE( eforce( ngw, nb_l ) )
ALLOCATE( eforce( ngw, nupdwn( ispin ) ) )
CALL dforce_all( ispin, cf(:,:,ispin_wfc), occ(:,ispin), eforce, &
vpot(:,ispin), eigr, bec, nupdwn, iupdwn )
CALL dforce_all( cf(:,:,ispin_wfc), occ(:,ispin), eforce, &
vpot(:,ispin), eigr, bec, nupdwn(ispin), iupdwn(ispin) )
CALL kohn_sham( ispin, cf(:,:,ispin_wfc), wfill, eforce, nupdwn )
CALL kohn_sham( ispin, cf(:,:,ispin_wfc), wfill, eforce, nupdwn, nb_l )
DEALLOCATE( eforce )
@ -228,28 +173,26 @@ CONTAINS
CALL wave_descriptor_init( wempt, ngw, wfill%ngwt, nupdwn_emp, nupdwn_emp, &
1, 1, nspin, 'gamma' , wfill%gzero )
ALLOCATE( ce( ngw, n_emp, nspin ) )
ALLOCATE( ce( ngw, n_emp * nspin ) )
exist = readempty( ce, wempt )
exist = readempty( ce, n_emp * nspin )
IF( .NOT. exist ) &
CALL errore( ' ks_states ', ' empty states file not found', 1 )
DO ispin = 1, nspin
nb_l = wempt%nbl( ispin )
IF( nupdwn_emp( ispin ) > 0 ) THEN
IF( nb_l > 0 ) THEN
ALLOCATE( fi( nupdwn_emp( ispin ) ) )
fi = 2.0d0 / nspin
ALLOCATE( fi( nb_l ) )
fi( 1:nb_l ) = 2.0d0 / nspin
ALLOCATE( eforce( ngw, nupdwn_emp( ispin ) ) )
ALLOCATE( eforce( ngw, nb_l ) )
CALL dforce_all( ce(:,iupdwn_emp(ispin):), fi(:), eforce, &
vpot(:,ispin), eigr, bec, nupdwn_emp(ispin), iupdwn_emp(ispin) )
CALL dforce_all( ispin, ce(:,:,ispin), fi(:), eforce, &
vpot(:,ispin), eigr, bec, nupdwn_emp, iupdwn_emp )
CALL kohn_sham( ispin, ce(:,:,ispin), wempt, eforce, nupdwn_emp )
CALL kohn_sham( ispin, ce(:,iupdwn_emp(ispin):), wempt, eforce, nupdwn_emp, n_emp_l )
DEALLOCATE( eforce )
DEALLOCATE( fi )
@ -276,17 +219,18 @@ CONTAINS
SUBROUTINE print_all_states(cf, wfill, ce, wempt )
USE kinds, ONLY: DP
USE mp_global, ONLY: intra_image_comm
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE wave_types, ONLY: wave_descriptor
USE control_flags, ONLY: force_pairing
USE kinds, ONLY: DP
USE mp_global, ONLY: intra_image_comm
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE wave_types, ONLY: wave_descriptor
USE control_flags, ONLY: force_pairing
USE electrons_module, ONLY: iupdwn_emp
IMPLICIT NONE
! ... declare subroutine arguments
COMPLEX(DP), INTENT(INOUT) :: cf(:,:,:), ce(:,:,:)
COMPLEX(DP), INTENT(INOUT) :: cf(:,:,:), ce(:,:)
TYPE (wave_descriptor), INTENT(IN) :: wfill, wempt
! ... declare other variables
@ -331,7 +275,7 @@ CONTAINS
IF( ( iks > 0 ) .AND. ( iks <= wempt%nbt( ispin ) ) ) THEN
file_name = TRIM( ks_emp_file ) // &
& trim(spin_name(ispin)) // trim( int_to_char( iks ) )
CALL print_ks_states( ce(:,iks,ispin), file_name )
CALL print_ks_states( ce(:,iupdwn_emp(ispin)+iks-1), file_name )
END IF
END DO
END IF
@ -355,10 +299,10 @@ CONTAINS
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE wave_types, ONLY: wave_descriptor, wave_descriptor_init
USE wave_functions, ONLY: kohn_sham
USE forces
USE pseudo_projector, ONLY: projector
USE electrons_base, ONLY: iupdwn, nupdwn, nspin
USE electrons_module, ONLY: n_emp, iupdwn_emp, nupdwn_emp
USE electrons_module, ONLY: n_emp, iupdwn_emp, nupdwn_emp, nb_l, n_emp_l
USE empty_states, ONLY: readempty
IMPLICIT NONE
@ -370,10 +314,10 @@ CONTAINS
REAL(DP), INTENT(IN) :: occ(:,:), bec(:,:)
REAL (DP) :: vpot(:,:)
COMPLEX(DP), ALLOCATABLE :: ce(:,:,:)
COMPLEX(DP), ALLOCATABLE :: ce(:,:)
! ... declare other variables
INTEGER :: i, ib, ig, ngw, nb_g, nb_l, iks, nb, ispin
INTEGER :: i, ib, ig, ngw, nb_g, iks, nb, ispin
LOGICAL :: tortho = .TRUE.
LOGICAL :: exist
CHARACTER(LEN=4) :: nom
@ -405,8 +349,8 @@ CONTAINS
ALLOCATE( eforce( ngw, nb, 2 ) )
CALL dforce_all( 1, cf(:,:,1), occ(:,1), eforce(:,:,1), vpot(:,1), eigr, bec, nupdwn, iupdwn )
CALL dforce_all( 2, cf(:,:,1), occ(:,2), eforce(:,:,2), vpot(:,2), eigr, bec, nupdwn, iupdwn )
CALL dforce_all( cf(:,:,1), occ(:,1), eforce(:,:,1), vpot(:,1), eigr, bec, nupdwn(1), iupdwn(1) )
CALL dforce_all( cf(:,:,1), occ(:,2), eforce(:,:,2), vpot(:,2), eigr, bec, nupdwn(2), iupdwn(2) )
DO i = 1, nupdwn(2)
eforce(:,i,1) = occ(i,1) * eforce(:,i,1) + occ(i,2) * eforce(:,i,2)
@ -415,7 +359,7 @@ CONTAINS
eforce(:,i,1) = occ(i,1) * eforce(:,i,1)
END DO
CALL kohn_sham( 1, cf(:,:,1), wfill, eforce(:,:,1), nupdwn )
CALL kohn_sham( 1, cf(:,:,1), wfill, eforce(:,:,1), nupdwn, nb_l )
DEALLOCATE( eforce )
@ -428,29 +372,27 @@ CONTAINS
CALL wave_descriptor_init( wempt, ngw, wfill%ngwt, nupdwn_emp, nupdwn_emp, &
1, 1, nspin, 'gamma' , wfill%gzero )
nb_l = wempt%nbl( 1 )
ALLOCATE( ce( ngw, n_emp * nspin ) )
ALLOCATE( ce( ngw, n_emp, nspin ) )
exist = readempty( ce, wempt )
exist = readempty( ce, n_emp * nspin )
IF( .NOT. exist ) &
CALL errore( ' ks_states ', ' empty states file not found', 1 )
IF( nb_l > 0 ) THEN
IF( n_emp > 0 ) THEN
ALLOCATE( fi( nb_l ) )
ALLOCATE( fi( n_emp ) )
fi = 2.0d0
ALLOCATE( eforce( ngw, nb_l, 1 ))
ALLOCATE( eforce( ngw, n_emp, 1 ))
CALL dforce_all( 1, ce(:,:,1), fi(:), eforce(:,:,1), vpot(:,1), eigr, bec, nupdwn_emp, iupdwn_emp )
CALL dforce_all( ce(:,:), fi(:), eforce(:,:,1), vpot(:,1), eigr, bec, nupdwn_emp(1), iupdwn_emp(1) )
CALL kohn_sham( 1, ce(:,:,1), wempt, eforce(:,:,1), nupdwn_emp )
CALL kohn_sham( 1, ce(:,:), wempt, eforce(:,:,1), nupdwn_emp, n_emp_l )
CALL dforce_all( 2, ce(:,:,2), fi(:), eforce(:,:,1), vpot(:,2), eigr, bec, nupdwn_emp, iupdwn_emp )
CALL dforce_all( ce(:,iupdwn_emp(2):), fi(:), eforce(:,:,1), vpot(:,2), eigr, bec, nupdwn_emp(2), iupdwn_emp(2) )
CALL kohn_sham( 2, ce(:,:,2), wempt, eforce(:,:,1), nupdwn_emp )
CALL kohn_sham( 2, ce(:,iupdwn_emp(2):), wempt, eforce(:,:,1), nupdwn_emp, n_emp_l )
DEALLOCATE( eforce )
DEALLOCATE( fi )

2
CPV/main.f90
View File

@ -506,7 +506,7 @@
! ... Here find Empty states eigenfunctions and eigenvalues
!
IF ( ttempst ) THEN
CALL empty( tortho, atoms0, c0, wfill, vpot, eigr )
CALL empty( nfi, tortho, c0, wfill, vpot, eigr, bec )
END IF
! ... dipole

106
CPV/nl.f90
View File

@ -8,13 +8,15 @@
#include "f_defs.h"
!
!
!----------------------------------------------
MODULE nl
!----------------------------------------------
USE kinds
USE spherical_harmonics
USE nl_base
IMPLICIT NONE
SAVE
@ -24,6 +26,7 @@
PUBLIC :: nlrh_m
!----------------------------------------------
CONTAINS
!----------------------------------------------
@ -40,7 +43,6 @@
! ... include modules
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
USE atoms_type_module, ONLY: atoms_type
USE control_flags, ONLY: force_pairing
USE pseudopotential, ONLY: nspnl
@ -69,8 +71,8 @@
REAL(DP), ALLOCATABLE :: btmp( :, :, : )
REAL(DP), ALLOCATABLE :: fion( :, : )
! end of declarations
! ----------------------------------------------
! ... end of declarations
!
DO iss = 1, cdesc%nspin
!
@ -110,102 +112,6 @@
RETURN
END FUNCTION nlrh_m
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE nlsm2_s( ispin, wnl, atoms, eigr, c, cdesc, g2, gx, dfnl, kk)
! this routine computes the derivatives of the Kleinman-Bylander
! factors fnl, to be used for Hellmann-Feynman forces evaluation
!
! ... declare modules
USE pseudopotential, ONLY: nspnl
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
USE atoms_type_module, ONLY: atoms_type
USE cell_base, ONLY: tpiba
USE uspp_param, only: nh, lmaxkb
USE uspp, only: nhtol, nhtolm, indv
IMPLICIT NONE
! ... declare subroutine arguments
COMPLEX(DP), INTENT(IN) :: c(:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
TYPE (projector), INTENT(OUT) :: dfnl
TYPE(atoms_type), INTENT(INOUT) :: atoms ! ions structure
REAL(DP), INTENT(IN) :: wnl(:,:,:), g2(:), gx(:,:)
COMPLEX(DP), INTENT(IN) :: eigr(:,:)
INTEGER, INTENT(IN) :: kk
INTEGER, INTENT(IN) :: ispin
! ... declare other variables
REAL(DP), ALLOCATABLE :: gwork(:,:)
INTEGER :: is, ia, igh, isa, ig, iss, ll, l, m, ngw, nb, lda, ldw, ldf
INTEGER :: iy, ih, iv
COMPLEX(DP), ALLOCATABLE :: auxc(:,:), gxtmp(:)
COMPLEX(DP), PARAMETER :: ONE = (1.0d0,0.0d0), ZERO = (0.0d0,0.0d0)
! ... (-i) * i^l
COMPLEX(DP), PARAMETER :: csign(0:3) = (/ (0.0d0,-1.0d0), &
(1.0d0,0.0d0), (0.0d0,1.0d0), (-1.0d0,0.0d0) /)
ngw = cdesc%ngwl
nb = cdesc%nbl( ispin )
IF( cdesc%gamma ) THEN
lda = 2*SIZE(c, 1)
ldw = 2*SIZE(c, 1)
ldf = SIZE(dfnl%r, 1) * SIZE(dfnl%r, 2)
dfnl%r = 0.0d0
ELSE
lda = SIZE(c, 1)
ldw = SIZE(c, 1)
ldf = SIZE(dfnl%c, 1) * SIZE(dfnl%c, 2)
dfnl%c = 0.0d0
END IF
ALLOCATE(gwork(ngw, (lmaxkb+1)**2 ), gxtmp(ngw))
CALL ylmr2( (lmaxkb+1)**2, ngw, gx, g2, gwork )
!
DO iy = 1, (lmaxkb+1)**2
gwork(1:ngw,iy) = tpiba * gx(kk,1:ngw) * gwork(1:ngw,iy)
END DO
iss = 1
SPECS: DO is = 1, nspnl
ALLOCATE(auxc(ngw,atoms%na(is)))
LM: DO ih = 1, nh( is )
iv = indv ( ih, is )
iy = nhtolm( ih, is )
ll = nhtol ( ih, is ) + 1
l = ll - 1
igh = ih
! write( 6, * ) 'DEBUG = ', SUM( wnl( :, iv, is ) ), SUM( gwork( :, iy ) )
gxtmp(1:ngw) = csign(l) * wnl(1:ngw,iv,is) * gwork(1:ngw,iy)
DO ia = 1, atoms%na(is)
auxc(1:ngw,ia) = gxtmp(1:ngw) * eigr(1:ngw,iss + ia - 1)
END DO
IF( cdesc%gamma ) THEN
CALL DGEMM('T', 'N', atoms%na(is), nb, 2*ngw, 1.0d0, auxc(1,1), 2*ngw, &
c(1,1), ldw, 0.0d0, dfnl%r(iss,igh,1), ldf)
ELSE
CALL ZGEMM('C', 'N', atoms%na(is), nb, ngw, one, auxc(1,1), ngw, &
c(1,1), ldw, zero, dfnl%c(iss,igh,1), ldf)
END IF
END DO LM
DEALLOCATE(auxc)
iss = iss + atoms%na(is)
END DO SPECS
IF( cdesc%gamma ) CALL DSCAL(size(dfnl%r),2.0d0,dfnl%r(1,1,1),1)
!write( 6, * ) 'DEBUG ==== ', SUM( dfnl%r )
DEALLOCATE(gwork, gxtmp)
RETURN
END SUBROUTINE nlsm2_s
!----------------------------------------------
END MODULE nl

186
CPV/nl_base.f90
View File

@ -7,94 +7,6 @@
!
#include "f_defs.h"
!
!----------------------------------
MODULE nl_base
!----------------------------------
USE kinds, ONLY : DP
IMPLICIT NONE
SAVE
PRIVATE
PUBLIC :: ene_nl
!----------------------------------
CONTAINS
!----------------------------------
REAL(DP) FUNCTION ene_nl(fnl, wsg, fi, nspnl, na)
! this function computes and returns the nonlocal contribution to total
! energy, for both Gamma point and Kpoints calculations
!
! ene_nl = (sum over ia,ib,igh,is) f(ib) wsg(igh,is) fnl(ia,igh,ib)**2
!
! fi(ib,1) = occupation numbers
! fnl(ia,igh,ib) = Kleinman-Bylander factor (see nlsm1)
! wsg(igh,is) = inverse denominator in KB's formula
!
! ia = index of ion
! ib = index of band
! igh = index of orbital
! is = index of atomic species
USE pseudo_projector, ONLY: projector
USE control_flags, ONLY: gamma_only
TYPE (projector), INTENT(IN) :: fnl
REAL(DP), INTENT(IN) :: wsg(:,:), fi(:)
INTEGER, INTENT(IN) :: nspnl
INTEGER, INTENT(IN) :: na(:)
INTEGER :: igh, isa, is, ia, ib, nb, ngh
REAL(DP) :: enl, fsum
COMPLEX(DP) :: tt
IF( gamma_only ) THEN
nb = MIN( SIZE( fnl%r, 3), SIZE(fi) )
ngh = SIZE( fnl%r, 2)
ELSE
nb = MIN( SIZE( fnl%c, 3), SIZE(fi) )
ngh = SIZE( fnl%c, 2)
END IF
enl=0.d0
DO igh = 1, ngh
DO ib = 1, nb
isa = 0
DO is = 1, nspnl
fsum = 0.0d0
IF( gamma_only ) THEN
DO ia = 1, na(is)
fsum = fsum + fnl%r(isa+ia,igh,ib)**2
END DO
ELSE
DO ia = 1, na(is)
tt = fnl%c(isa+ia,igh,ib)
fsum = fsum + DBLE( CONJG(tt) * tt )
END DO
END IF
enl = enl + fi(ib) * wsg(igh, is) * fsum
isa = isa + na(is)
END DO
END DO
END DO
ene_nl = enl
RETURN
END FUNCTION ene_nl
!----------------------------------
END MODULE nl_base
!----------------------------------
!-----------------------------------------------------------------------
subroutine nlsm1 ( n, nspmn, nspmx, eigr, c, becp )
!-----------------------------------------------------------------------
@ -126,14 +38,15 @@
real(DP), intent(in) :: eigr( 2, ngw, nat ), c( 2, ngw, n )
real(DP), intent(out) :: becp( nkb, n )
!
complex(DP), allocatable :: wrk2(:,:)
!
integer :: isa, ig, is, iv, ia, l, ixr, ixi, inl, i
integer :: isa, ig, is, iv, ia, l, ixr, ixi, inl, i, nhx
real(DP) :: signre, signim, arg
real(DP), allocatable :: becps( :, : )
real(DP), allocatable :: wrk2( :, :, : )
!
call start_clock( 'nlsm1' )
allocate( wrk2( ngw, nax ) )
allocate( wrk2( 2, ngw, nax ) )
isa = 0
do is = 1, nspmn - 1
@ -141,6 +54,13 @@
end do
do is = nspmn, nspmx
!
IF( nproc_image > 1 ) THEN
nhx = nh( is ) * na( is )
IF( MOD( nhx, 2 ) /= 0 ) nhx = nhx + 1
ALLOCATE( becps( nhx, n ) )
becps = 0.0d0
END IF
!
do iv = 1, nh( is )
!
@ -173,36 +93,55 @@
! q = 0 component (with weight 1.0)
!
if (gstart == 2) then
wrk2(1,ia)= &
CMPLX(signre*beta(1,iv,is)*eigr(ixr,1,ia+isa),signim*beta(1,iv,is)*eigr(ixi,1,ia+isa) )
wrk2( 1, 1, ia ) = signre*beta(1,iv,is)*eigr(ixr,1,ia+isa)
wrk2( 2, 1, ia ) = signim*beta(1,iv,is)*eigr(ixi,1,ia+isa)
end if
!
! q > 0 components (with weight 2.0)
!
do ig = gstart, ngw
arg = 2.0*beta(ig,iv,is)
wrk2(ig,ia) = &
CMPLX(signre*arg*eigr(ixr,ig,ia+isa),signim*arg*eigr(ixi,ig,ia+isa) )
arg = 2.0d0 * beta(ig,iv,is)
wrk2( 1, ig, ia ) = signre*arg*eigr(ixr,ig,ia+isa)
wrk2( 2, ig, ia ) = signim*arg*eigr(ixi,ig,ia+isa)
end do
!
end do
inl=ish(is)+(iv-1)*na(is)+1
!
call MXMA( wrk2, 2*ngw, 1, c, 1, 2*ngw, becp(inl,1), 1, nkb, &
na(is), 2*ngw, n )
IF( nproc_image > 1 ) THEN
inl=(iv-1)*na(is)+1
! CALL MXMA( wrk2, 2*ngw, 1, c, 1, 2*ngw, becps( inl, 1 ), 1, nhx, na(is), 2*ngw, n )
CALL DGEMM( 'T', 'N', na(is), n, 2*ngw, 1.0d0, wrk2, 2*ngw, c, 2*ngw, 0.0d0, becps( inl, 1 ), nhx )
! subroutine mxma (a,na,iad,b,nb,ibd,c,nc,icd,nar,nac,nbc)
! call DGEMM(mode1,mode2,nar,nbc,nac,1.d0,a,lda,b,ldb,0.d0,c,icd)
ELSE
inl=ish(is)+(iv-1)*na(is)+1
! call MXMA( wrk2, 2*ngw, 1, c, 1, 2*ngw, becp( inl, 1 ), 1, nkb, na(is), 2*ngw, n )
CALL DGEMM( 'T', 'N', na(is), n, 2*ngw, 1.0d0, wrk2, 2*ngw, c, 2*ngw, 0.0d0, becp( inl, 1 ), nkb )
END IF
end do
IF( nproc_image > 1 ) THEN
inl=ish(is)+1
do i=1,n
CALL mp_sum( becp( inl : (inl + na(is)*nh(is) - 1), i ), intra_image_comm )
!
inl = ish(is) + 1
!
CALL mp_sum( becps, intra_image_comm )
do i = 1, n
do iv = inl , ( inl + na(is) * nh(is) - 1 )
becp( iv, i ) = becps( iv - inl + 1, i )
end do
end do
DEALLOCATE( becps )
END IF
isa = isa + na(is)
isa = isa + na(is)
end do
@ -245,8 +184,8 @@
real(DP), intent(out) :: becdr(nkb,n,3)
logical, intent(in) :: tred
!
real(DP), allocatable :: gk(:)
complex(DP), allocatable :: wrk2(:,:)
real(DP), allocatable :: gk(:)
real(DP), allocatable :: wrk2(:,:,:)
!
integer :: ig, is, iv, ia, k, l, ixr, ixi, inl, isa, i
real(DP) :: signre, signim, arg
@ -254,7 +193,7 @@
call start_clock( 'nlsm2' )
allocate( gk( ngw ) )
allocate( wrk2( ngw, nax ) )
allocate( wrk2( 2, ngw, nax ) )
becdr = 0.d0
!
@ -296,20 +235,21 @@
endif
!
do ia=1,na(is)
! q = 0 component (with weight 1.0)
if (gstart == 2) then
! q = 0 component (with weight 1.0)
wrk2(1,ia) = &
CMPLX (signre*gk(1)*beta(1,iv,is)*eigr(ixr,1,ia+isa),signim*gk(1)*beta(1,iv,is)*eigr(ixi,1,ia+isa) )
! q > 0 components (with weight 2.0)
wrk2(1,1,ia) = signre*gk(1)*beta(1,iv,is)*eigr(ixr,1,ia+isa)
wrk2(2,1,ia) = signim*gk(1)*beta(1,iv,is)*eigr(ixi,1,ia+isa)
end if
! q > 0 components (with weight 2.0)
do ig=gstart,ngw
arg = 2.0*gk(ig)*beta(ig,iv,is)
wrk2(ig,ia) = &
CMPLX (signre*arg*eigr(ixr,ig,ia+isa),signim*arg*eigr(ixi,ig,ia+isa) )
wrk2(1,ig,ia) = signre*arg*eigr(ixr,ig,ia+isa)
wrk2(2,ig,ia) = signim*arg*eigr(ixi,ig,ia+isa)
end do
end do
inl=ish(is)+(iv-1)*na(is)+1
call MXMA(wrk2,2*ngw,1,c,1,2*ngw,becdr(inl,1,k),1, nkb,na(is),2*ngw,n)
! call MXMA(wrk2,2*ngw,1,c,1,2*ngw,becdr(inl,1,k),1, nkb,na(is),2*ngw,n)
CALL DGEMM( 'T', 'N', na(is), n, 2*ngw, 1.0d0, wrk2, 2*ngw, c, 2*ngw, 0.0d0, becdr( inl, 1, k ), nkb )
end do
isa = isa + na(is)
@ -317,9 +257,7 @@
end do
IF( tred .AND. ( nproc_image > 1 ) ) THEN
DO i = 1, n
CALL mp_sum( becdr(:,i,k), intra_image_comm )
END DO
CALL mp_sum( becdr(:,:,k), intra_image_comm )
END IF
end do
@ -559,12 +497,12 @@ SUBROUTINE caldbec( ngw, nkb, n, nspmn, nspmx, eigr, c, dbec, tred )
real(DP), intent(out) :: dbec( nkb, n, 3, 3 )
logical, intent(in) :: tred
!
complex(DP), allocatable :: wrk2(:,:)
real(DP), allocatable :: wrk2(:,:,:)
!
integer :: ig, is, iv, ia, l, ixr, ixi, inl, i, j, ii, isa
real(DP) :: signre, signim, arg
!
allocate( wrk2( ngw, nax ) )
allocate( wrk2( 2, ngw, nax ) )
!
!
do j=1,3
@ -605,19 +543,19 @@ SUBROUTINE caldbec( ngw, nkb, n, nspmn, nspmx, eigr, c, dbec, tred )
do ia=1,na(is)
if (gstart == 2) then
! q = 0 component (with weight 1.0)
wrk2(1,ia)= &
CMPLX(signre*dbeta(1,iv,is,i,j)*eigr(ixr,1,ia+isa),signim*dbeta(1,iv,is,i,j)*eigr(ixi,1,ia+isa) )
wrk2(1,1,ia)= signre*dbeta(1,iv,is,i,j)*eigr(ixr,1,ia+isa)
wrk2(2,1,ia)= signim*dbeta(1,iv,is,i,j)*eigr(ixi,1,ia+isa)
end if
! q > 0 components (with weight 2.0)
do ig = gstart, ngw
arg = 2.0*dbeta(ig,iv,is,i,j)
wrk2(ig,ia) = &
CMPLX(signre*arg*eigr(ixr,ig,ia+isa),signim*arg*eigr(ixi,ig,ia+isa) )
wrk2(1,ig,ia) = signre*arg*eigr(ixr,ig,ia+isa)
wrk2(2,ig,ia) = signim*arg*eigr(ixi,ig,ia+isa)
end do
end do
inl=ish(is)+(iv-1)*na(is)+1
call MXMA(wrk2,2*ngw,1,c,1,2*ngw,dbec(inl,1,i,j),1, &
& nkb,na(is),2*ngw,n)
! call MXMA(wrk2,2*ngw,1,c,1,2*ngw,dbec(inl,1,i,j),1,nkb,na(is),2*ngw,n)
CALL DGEMM( 'T', 'N', na(is), n, 2*ngw, 1.0d0, wrk2, 2*ngw, c, 2*ngw, 0.0d0, dbec( inl, 1, i, j ), nkb )
end do
if( ( nproc_image > 1 ) .AND. tred ) then
inl=ish(is)+1

506
CPV/optical.f90
View File

@ -5,415 +5,272 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
#include "f_defs.h"
MODULE optical_properties
USE kinds
#include "f_defs.h"
MODULE optical_properties
USE kinds, ONLY: DP
IMPLICIT NONE
SAVE
PRIVATE
INTEGER :: nfreq
REAL(DP) :: maxdie ! Hartree units
REAL(DP) :: ddie ! Hartree units
REAL(DP) :: temperature ! Kelvin
REAL(DP), PARAMETER :: small_freq = 1.d-6 ! Hartree units
INTEGER :: nfreq
REAL(DP) :: maxdie = 10.d0 ! Hartree units
REAL(DP) :: ddie = 0.01d0 ! Hartree units
REAL(DP) :: temperature = 300.0d0 ! Kelvin
COMPLEX(DP), ALLOCATABLE :: dielec_total(:)
REAL(DP), ALLOCATABLE :: sigma_total(:)
INTEGER, ALLOCATABLE :: n_total(:)
PUBLIC :: optical_setup, optical_closeup, opticalp, write_dielec
PUBLIC :: opticalp
CONTAINS
SUBROUTINE optical_setup(woptical, noptical, boptical)
USE constants, ONLY: au
REAL(DP), INTENT(IN) :: woptical
REAL(DP), INTENT(IN) :: boptical
INTEGER, INTENT(IN) :: noptical
IF( noptical < 1 ) THEN
CALL errore(' optical_properties: optical_setup ',' noptical out of range ',noptical)
END IF
IF( woptical < small_freq ) THEN
CALL errore(' optical_properties: optical_setup ',' woptical out of range ',INT(woptical))
END IF
nfreq = noptical
maxdie = woptical / au
ddie = maxdie / DBLE(nfreq)
temperature = boptical
ALLOCATE( dielec_total(nfreq), sigma_total(nfreq), n_total(nfreq) )
dielec_total = 0.0d0
sigma_total = 0.0d0
n_total = 0
RETURN
END SUBROUTINE optical_setup
CONTAINS
SUBROUTINE optical_closeup
IF( ALLOCATED( dielec_total ) ) DEALLOCATE( dielec_total )
IF( ALLOCATED( sigma_total ) ) DEALLOCATE( sigma_total )
IF( ALLOCATED( n_total ) ) DEALLOCATE( n_total )
RETURN
END SUBROUTINE optical_closeup
SUBROUTINE opticalp(nfi, box, atoms, c0, wfill, occ, ce, wempt, vpot, eigr, bec )
SUBROUTINE opticalp(nfi, omega, c0, wfill, ce, wempt, vpot, eigr, bec )
USE cell_module, ONLY: boxdimensions
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector, allocate_projector, deallocate_projector
USE pseudopotential, ONLY: nsanl
USE nl, ONLY: nlsm1_s
USE forces, ONLY: dforce_all
USE brillouin, ONLY: kpoints, kp
USE electrons_module, ONLY: ei, ei_emp
USE electrons_base, ONLY: nupdwn, iupdwn
USE kohn_sham_states, ONLY: kohn_sham
USE constants, ONLY: au, pi, k_boltzman_au, au_to_ohmcmm1
USE cell_base, ONLY: tpiba2
USE mp, ONLY: mp_sum
USE mp_global, ONLY: intra_image_comm
USE io_global, ONLY: ionode
USE atoms_type_module, ONLY: atoms_type
USE io_files, ONLY: dielecunit, dielecfile
USE control_flags, ONLY: force_pairing
USE wave_types, ONLY: wave_descriptor
USE pseudopotential, ONLY: nsanl
USE forces, ONLY: dforce_all
USE electrons_module, ONLY: ei, ei_emp, n_emp, nupdwn_emp, iupdwn_emp, nb_l, n_emp_l
USE electrons_base, ONLY: nupdwn, iupdwn, nspin
USE wave_functions, ONLY: kohn_sham
USE constants, ONLY: autoev, pi, k_boltzman_au, au_to_ohmcmm1, eps8
USE cell_base, ONLY: tpiba2
USE mp, ONLY: mp_sum
USE mp_global, ONLY: intra_image_comm
USE io_global, ONLY: ionode, stdout
USE control_flags, ONLY: force_pairing
USE reciprocal_vectors, ONLY: gx, g
USE reciprocal_space_mesh, ONLY: gkx_l
USE pseudopotential, ONLY: nspnl
USE uspp_param, ONLY: nhm
USE uspp, ONLY: nkb
USE pseudopotential, ONLY: nspnl
USE uspp_param, ONLY: nhm
USE uspp, ONLY: nkb
INTEGER, INTENT(IN) :: nfi
TYPE(boxdimensions), INTENT(IN) :: box
TYPE(atoms_type), INTENT(INOUT) :: atoms ! ions structure
COMPLEX(DP), INTENT(IN) :: c0(:,:,:,:)
COMPLEX(DP), INTENT(INOUT) :: ce(:,:,:,:)
INTEGER, INTENT(IN) :: nfi
REAL(DP), INTENT(IN) :: omega
COMPLEX(DP), INTENT(IN) :: c0(:,:,:)
COMPLEX(DP), INTENT(INOUT) :: ce(:,:,:)
TYPE(wave_descriptor), INTENT(IN) :: wempt, wfill
REAL(DP), INTENT(IN) :: occ(:,:,:)
REAL (DP), INTENT(in) :: vpot(:,:)
REAL (DP) :: bec(:,:)
COMPLEX(DP) :: eigr(:,:)
INTEGER :: nspin, ispin, ngw, nb_l, nk, ngw_g
INTEGER :: ie, if, nf, ne, ik, idie, ig, ierr
COMPLEX(DP), ALLOCATABLE :: eforce(:,:,:)
INTEGER :: iss, ngw
INTEGER :: ie, if, nf, ne, idie, ig, ierr, i, nbt
COMPLEX(DP), ALLOCATABLE :: eforce(:,:)
REAL (DP), ALLOCATABLE :: bece(:,:)
REAL(DP) :: curr(3), currt, wef
REAL(DP) :: curr(3), currt, wef, w, wg2, p2
COMPLEX(DP) :: ccurr(3), ccurrt
COMPLEX(DP), ALLOCATABLE :: diet(:), cf(:,:,:,:)
COMPLEX(DP), ALLOCATABLE :: diet(:), cf(:,:,:)
INTEGER :: cif, cie
REAL(DP), ALLOCATABLE :: sigma(:), fi(:), eig(:), ff(:,:)
REAL(DP), ALLOCATABLE :: sigma(:), fi(:), eig(:), ff(:), dipole(:,:,:,:)
REAL(DP), ALLOCATABLE :: epsilon2(:), dos(:)
INTEGER, ALLOCATABLE :: ndiet(:)
REAL(DP) :: FACT, ef, beta
LOGICAL :: gamma_symmetry, gzero
INTEGER :: nupdwn_emp( wfill%nspin )
INTEGER :: iupdwn_emp( wfill%nspin )
REAL(DP) :: FACT, ef, beta, sumrule(3)
! ... SUBROUTINE BODY
CALL errore( ' opticalp ', ' not working ', 1 )
! ... SUBROUTINE BODY
IF( force_pairing ) &
CALL errore( ' opticalp ', ' force_pairing not implemented ', 1 )
ALLOCATE( cf( SIZE(c0,1), SIZE(c0,2), SIZE(c0,3), SIZE(c0,4) ) )
cf = c0
nk = wfill%nkl
nspin = wfill%nspin
nfreq = INT( maxdie / ddie + 1.0d0 )
!
beta = 1.0d0 / ( k_boltzman_au * temperature )
!
fact = 2.0d0 * pi / ( 3.0d0 * omega )
ngw = wfill%ngwl
nbt = wfill%nbl( 1 ) + wempt%nbl( 1 )
IF( nspin > 1 ) nbt = MAX( nbt, wfill%nbl( 2 ) + wempt%nbl( 2 ) )
ALLOCATE( diet(nfreq), ndiet(nfreq), sigma(nfreq) )
ALLOCATE( cf( SIZE(c0,1), SIZE(c0,2), SIZE(c0,3) ) )
ALLOCATE( dipole( 3, nbt, nbt, nspin ) )
ALLOCATE( epsilon2( nfreq ) )
ALLOCATE( dos( nfreq ) )
diet = 0.0d0
sigma = 0.0d0
dielec_total = 0.0d0
sigma_total = 0.0d0
ndiet = 0
dipole = 0.0d0
epsilon2 = 0.0d0
dos = 0.0d0
DO ispin = 1, nspin
cf = c0
ngw = wfill%ngwl
nb_l = wfill%nbl( ispin )
DO iss = 1, nspin
ALLOCATE( eforce( ngw, nb_l, nk ) )
!
CALL nlsm1 ( nb_l, 1, nspnl, eigr, cf(1,1,1,ispin), bec )
ALLOCATE( ff( MAX( nupdwn( iss ), nupdwn_emp( iss ) ) ) )
CALL dforce_all( ispin, cf(:,:,1,ispin), occ(:,1,ispin), eforce(:,:,1), &
vpot(:,ispin), eigr, bec, nupdwn, iupdwn )
ff = 2.0d0 / nspin
CALL kohn_sham( ispin, cf(:,:,1,ispin), wfill, eforce(:,:,1) )
ALLOCATE( eforce( ngw, nupdwn( iss ) ) )
!
CALL nlsm1 ( nupdwn( iss ), 1, nspnl, eigr, cf(1,1,iss), bec )
ngw = wempt%ngwl
nb_l = wempt%nbl( ispin )
nupdwn_emp = nb_l
iupdwn_emp(1) = 1
IF( nspin == 2 ) iupdwn_emp(2) = 1+nb_l
CALL dforce_all( cf(:,:,iss), ff, eforce, vpot(:,iss), eigr, bec, &
nupdwn(iss), iupdwn(iss) )
CALL kohn_sham( iss, cf(:,:,iss), wfill, eforce, nupdwn, nb_l )
ALLOCATE( ff( nb_l, nk ) )
DO ik = 1, nk
ff( 1:nb_l, ik ) = 2.0d0 / nspin
END DO
DEALLOCATE( eforce )
ALLOCATE( bece( nkb, nb_l ) )
!
CALL nlsm1 ( nb_l, 1, nspnl, eigr, ce(1,1,1,ispin), bece )
!
CALL dforce_all( ispin, ce(:,:,1,ispin), ff( :, 1), eforce(:,:,1), &
vpot(:,ispin), eigr, bece, nupdwn_emp, iupdwn_emp )
!
CALL kohn_sham( ispin, ce(:,:,1,ispin), wempt, eforce(:,:,1) )
ALLOCATE( bece( nkb, nupdwn_emp( iss ) ) )
ALLOCATE( eforce( ngw, nupdwn_emp( iss ) ) )
!
CALL nlsm1 ( nupdwn_emp( iss ), 1, nspnl, eigr, ce(1,1,iss), bece )
!
CALL dforce_all( ce(:,:,iss), ff, eforce, vpot(:,iss), eigr, bece, &
nupdwn_emp(iss), iupdwn_emp(iss) )
!
CALL kohn_sham( iss, ce(:,:,iss), wempt, eforce, nupdwn_emp, n_emp_l )
DEALLOCATE( bece )
DEALLOCATE( eforce )
DEALLOCATE( ff )
DEALLOCATE( bece )
DEALLOCATE( ff )
END DO
fact = 2.0d0 * pi / ( 3.0d0 * box%deth )
gamma_symmetry = wfill%gamma
ierr = 0
IF( ionode ) THEN
OPEN(UNIT=dielecunit, FILE=dielecfile, STATUS='unknown', POSITION='append', IOSTAT=ierr)
END IF
CALL mp_sum( ierr, intra_image_comm )
IF( ierr /= 0 ) &
CALL errore(' opticalp ', ' opening file '//TRIM(dielecfile), 1 )
#if defined __CONDUCTIVITY
WRITE( stdout, * ) 'epsilon: '
WRITE( dielecunit, * ) 'STEP: ',nfi
sumrule = 0.0d0
DO ispin = 1, nspin
DO iss = 1, nspin
nf = wfill%nbl( ispin )
ne = wempt%nbl( ispin )
ngw = wfill%ngwl
nk = kp%nkpt
nf = wfill%nbl( iss )
ne = wempt%nbl( iss )
ALLOCATE( fi( nf + ne ), eig( nf + ne ) )
DO ik = 1, nk
ef = ( ei_emp(1,iss) + ei(nf,iss) ) / 2.0
ef = ( ei_emp(1,ik,ispin) + ei(nf,ik,ispin) ) / 2.0
DO if = 1, nf
fi( if ) = 1.0 / nspin / ( exp( beta * (ei(if,ik,ispin) - ef) ) + 1 )
eig( if ) = ei( if, ik, ispin )
IF( ionode ) WRITE( dielecunit, fmt = '(I4,2F12.6)' ) if, fi(if), eig(if)
END DO
DO ie = nf+1, ne+nf
fi( ie ) = 1.0 / nspin / ( exp( beta * (ei_emp(ie-nf,ik,ispin) - ef) ) + 1 )
eig( ie ) = ei_emp( ie-nf, ik, ispin )
IF( ionode ) WRITE( dielecunit, fmt = '(I4,2F12.6)' ) ie, fi(ie), eig(ie)
END DO
DO if = 1, nf
fi( if ) = 2.0 / nspin
eig( if ) = ei( if, iss )
! IF( ionode ) WRITE( stdout, fmt = '(I4,2F12.6)' ) if, fi(if), eig(if)
END DO
DO if = 1, (nf + ne - 1)
DO ie = if + 1, (nf + ne)
DO ie = nf+1, ne+nf
fi( ie ) = 0.0d0
eig( ie ) = ei_emp( ie-nf, iss )
! IF( ionode ) WRITE( stdout, fmt = '(I4,2F12.6)' ) ie, fi(ie), eig(ie)
END DO
DO if = 1, nf
!
DO ie = nf + 1, (nf + ne)
! frequencies in atomic units
!
wef = eig(ie) - eig(if)
!
! discretize the frequency
!
idie = wef / ddie + 1
IF( (wef > small_freq) .AND. (idie <= nfreq) ) THEN
IF( wef > eps8 ) THEN
IF( ie <= nf ) THEN
cie = ie
ELSE
cie = ie-nf
END IF
IF( if <= nf ) THEN
cif = if
ELSE
cif = if-nf
END IF
IF( gamma_symmetry ) THEN
curr = 0.0d0
ccurr = 0.0d0
DO ig = 1, ngw
curr(1) = curr(1) + gx(1,ig) * &
AIMAG( ce( ig, cie, ik, ispin ) * CONJG( cf( ig, cif, ik, ispin ) )
curr(2) = curr(2) + gx(2,ig) * &
AIMAG( ce( ig, cie, ik, ispin ) * CONJG( cf( ig, cif, ik, ispin ) )
curr(3) = curr(3) + gx(3,ig) * &
AIMAG( ce( ig, cie, ik, ispin ) * CONJG( cf( ig, cif, ik, ispin ) )
ccurr(1) = ccurr(1) + gx(1,ig) * &
ce( ig, cie, iss ) * CONJG( cf( ig, cif, iss ) )
ccurr(2) = ccurr(2) + gx(2,ig) * &
ce( ig, cie, iss ) * CONJG( cf( ig, cif, iss ) )
ccurr(3) = ccurr(3) + gx(3,ig) * &
ce( ig, cie, iss ) * CONJG( cf( ig, cif, iss ) )
END DO
! parallel sum of curr
CALL mp_sum( curr, intra_image_comm )
!
CALL mp_sum( ccurr, intra_image_comm )
!
! the factor 4.0d0 accounts for gamma symmetry
currt = 4.0d0 * (fi(if)-fi(ie)) * ( curr(1)**2 + curr(2)**2 + curr(3)**2 )
!
wg2 = 4.0d0
!
curr = AIMAG( ccurr )
!
!dipole( :, if, ie, iss ) = wg2 * tpiba2 * DBLE( ccurr(:) * CONJG( ccurr(:) ) )
!dipole( :, ie, if, iss ) = wg2 * tpiba2 * DBLE( ccurr(:) * CONJG( ccurr(:) ) )
dipole( :, if, ie, iss ) = wg2 * tpiba2 * curr(:)**2
dipole( :, ie, if, iss ) = wg2 * tpiba2 * curr(:)**2
!
p2 = DBLE( dipole( 1, if, ie, iss ) + dipole( 2, if, ie, iss ) + dipole( 3, if, ie, iss ) )
!
!
currt = wg2 * (fi(if)-fi(ie)) * ( curr(1)**2 + curr(2)**2 + curr(3)**2 )
currt = currt * tpiba2 / wef
!
! update dielectric tensor
diet(idie) = diet(idie) + CMPLX(0.0d0, currt) / wef
sigma(idie) = sigma(idie) + currt
ndiet(idie) = ndiet(idie) + 1
END IF
!
IF( idie <= nfreq ) THEN
diet(idie) = diet(idie) + CMPLX(0.0d0, currt) / wef
sigma(idie) = sigma(idie) + currt
ndiet(idie) = ndiet(idie) + 1
epsilon2(idie) = epsilon2(idie) + 4.0d0 * pi * pi * p2 * fi(if) / wef**2 / 3.0d0
END IF
sumrule = sumrule + fi(if) * 2.0d0 * dipole( :, if, ie, iss ) / wef
END IF
END DO
END DO
! DO if = 1, nf
! DO ie = 1, ne
! ! frequencies in atomic units
! wef = ei_emp(ie,ik,ispin) - ei(if,ik,ispin)
! ! discretize the frequency
! idie = wef / ddie + 1
! IF( (wef .GT. small_freq) .AND. (idie .LT. nfreq) ) THEN
! IF( gamma_symmetry ) THEN
! curr = 0.0d0
! DO ig = 1, ngw
! curr(1) = curr(1) + gx(1,ig) * &
! AIMAG( ce(ik,ispin)%w(ig, ie) * CONJG( cf(ik,ispin)%w(ig, if) ) )
! curr(2) = curr(2) + gx(2,ig) * &
! AIMAG( ce(ik,ispin)%w(ig, ie) * CONJG( cf(ik,ispin)%w(ig, if) ) )
! curr(3) = curr(3) + gx(3,ig) * &
! AIMAG( ce(ik,ispin)%w(ig, ie) * CONJG( cf(ik,ispin)%w(ig, if) ) )
! END DO
! ! parallel sum of curr
! CALL mp_sum( curr, group )
! ! the factor 4.0d0 accounts for gamma symmetry
! currt = 4.0d0 * ( curr(1)**2 + curr(2)**2 + curr(3)**2 )
! currt = currt * tpiba2 ! / wef
! ! update dielectric tensor
! diet(idie) = diet(idie) + CMPLX(0.0d0, currt) ! / wef
! sigma(idie) = sigma(idie) + currt
! ndiet(idie) = ndiet(idie) + 1
! END IF
! END IF
! END DO
! END DO
!
END DO
!
END DO
DEALLOCATE( fi, eig )
END DO
#else
WRITE( dielecunit, * ) 'STEP: ',nfi
DO ispin = 1, nspin
nf = wfill%nbl( ispin )
ne = wempt%nbl( ispin )
ngw = wfill%ngwl
nk = kp%nkpt
DO ik = 1, nk
DO ie = 1, ne
DO if = 1, nf
wef = ei_emp(ie,ik,ispin) - ei(if,ik,ispin)
IF( gamma_symmetry ) THEN
curr = 0.0d0
DO ig = 1, ngw
curr(1) = curr(1) + gx(1,ig) * &
AIMAG( ce( ig, ie, ik, ispin ) * CONJG( cf( ig, if, ik, ispin ) ) )
curr(2) = curr(2) + gx(2,ig) * &
AIMAG( ce( ig, ie, ik, ispin ) * CONJG( cf( ig, if, ik, ispin ) ) )
curr(3) = curr(3) + gx(3,ig) * &
AIMAG( ce( ig, ie, ik, ispin ) * CONJG( cf( ig, if, ik, ispin ) ) )
END DO
CALL mp_sum( curr, intra_image_comm )
currt = 2.0d0 * ( curr(1)**2 + curr(2)**2 + curr(3)**2 )
ELSE
ccurr = 0.0d0
DO ig = 1, ngw
ccurr(1) = ccurr(1) + gkx_l(1, ig, ik) * &
ce( ig, ie, ik, ispin ) * CONJG( cf( ig, if, ik, ispin ) )
ccurr(2) = ccurr(2) + gkx_l(2, ig, ik) * &
ce( ig, ie, ik, ispin ) * CONJG( cf( ig, if, ik, ispin ) )
ccurr(3) = ccurr(3) + gkx_l(3, ig, ik) * &
ce( ig, ie, ik, ispin ) * CONJG( cf( ig, if, ik, ispin ) )
END DO
CALL mp_sum( ccurr, intra_image_comm )
ccurrt = ccurr(1)*CONJG(ccurr(1)) + ccurr(2)*CONJG(ccurr(2)) + ccurr(3)*CONJG(ccurr(3))
WRITE( dielecunit ,100 ) ispin, ik, ie, if, wef, ccurrt
100 FORMAT(4I5,1D14.6,3X,2D14.6)
END IF
END DO
END DO
END DO
END DO
diet = fact * diet
sigma = fact * sigma
#endif
! ... write to stdout
ierr = 0
IF( ionode ) THEN
CLOSE(UNIT=dielecunit, IOSTAT=ierr)
WRITE( stdout,*) ' sumrule = ', sumrule(1)
WRITE( stdout,*) ' sumrule = ', sumrule(2)
WRITE( stdout,*) ' sumrule = ', sumrule(3)
WRITE( stdout,*) ' 2p/omega = ', 2.0d0 * pi / omega
WRITE( stdout,*) ' sumrule = ', sumrule(1) * 2.0d0 * pi / omega
WRITE( stdout,*) ' sumrule = ', sumrule(2) * 2.0d0 * pi / omega
WRITE( stdout,*) ' sumrule = ', sumrule(3) * 2.0d0 * pi / omega
DO i = 1, nfreq
IF( ndiet(i) > 0 ) THEN
w = (DBLE(i)-0.5d0) * ddie
WRITE( 100, * ) w * autoev, epsilon2(i) / omega
END IF
END DO
!WRITE( stdout,40)
!WRITE( stdout,50)
!DO I = 1, SIZE(diet)
! w = (DBLE(i)-0.5d0) * ddie
! WRITE( stdout,110) w * autoev, sigma(i) * au_to_ohmcmm1 / ddie, ndiet(i)
! ! WRITE(stdout,100) &
! ! w * autoev, diet(i) / w / w, sigma(i) * au_to_ohmcmm1 / w, ndiet(i)
!END DO
END IF
CALL mp_sum( ierr, intra_image_comm )
IF( ierr /= 0 ) &
CALL errore(' opticalp ', ' opening file '//TRIM(dielecfile), 1 )
! accumulate statistical values
WHERE( ndiet > 0 )
dielec_total = fact * diet
sigma_total = fact * sigma
n_total = ndiet
END WHERE
DEALLOCATE( diet, ndiet, sigma )
DEALLOCATE( cf )
RETURN
END SUBROUTINE opticalp
SUBROUTINE write_dielec (nfi, tm)
USE constants, ONLY: au, au_to_ohmcmm1
USE io_files, ONLY: dielecunit, dielecfile
USE io_global, ONLY: ionode
USE mp_global, ONLY: intra_image_comm
USE mp, ONLY: mp_sum
INTEGER, INTENT(IN) :: nfi
REAL(DP), INTENT(IN) :: tm
REAL(DP) :: w
INTEGER :: i, ierr
#if defined __CONDUCTIVITY
ierr = 0
IF( ionode ) THEN
OPEN(UNIT=dielecunit, FILE=dielecfile, STATUS='unknown', POSITION='append', IOSTAT=ierr)
END IF
CALL mp_sum( ierr, intra_image_comm )
IF( ierr /= 0 ) &
CALL errore(' write_dielec ', ' opening file '//TRIM(dielecfile), 1 )
WRITE( dielecunit, 30 ) nfi, tm
DO I = 1, SIZE(dielec_total)
w = (DBLE(i)-0.5d0) * ddie
! WRITE(dielecunit,100) &
! w * au, dielec_total(i) / w / w, sigma_total(i) * au_to_ohmcmm1 / w, n_total(i)
WRITE(dielecunit,100) &
w * au, dielec_total(i), sigma_total(i) * au_to_ohmcmm1, n_total(i)
END DO
ierr = 0
IF( ionode ) THEN
CLOSE(UNIT=dielecunit, IOSTAT=ierr)
END IF
CALL mp_sum( ierr, intra_image_comm )
IF( ierr /= 0 ) &
CALL errore(' write_dielec ', ' closing file '//TRIM(dielecfile), 1 )
! ... write to stdout
WRITE( stdout,40)
WRITE( stdout,50)
DO I = 1, SIZE(dielec_total)
w = (DBLE(i)-0.5d0) * ddie
! WRITE( stdout,110) w * au, sigma_total(i) * au_to_ohmcmm1 / w, n_total(i)
WRITE( stdout,110) w * au, sigma_total(i) * au_to_ohmcmm1 / ddie, n_total(i)
END DO
#endif
DEALLOCATE( dipole )
DEALLOCATE( dos )
30 FORMAT(2X,'STEP:',I7,1X,F10.2)
40 FORMAT(/,3X,'Frequency dependent electronic conductivity',/)
@ -421,8 +278,11 @@
90 FORMAT(3X,'Dielectric function')
100 FORMAT(3X,F12.6,3D16.8,I5)
110 FORMAT(3X,F12.6,D16.8,I5)
RETURN
END SUBROUTINE write_dielec
END MODULE optical_properties
RETURN
END SUBROUTINE opticalp
END MODULE optical_properties

14
CPV/ortho.f90
View File

@ -214,7 +214,8 @@
SUBROUTINE ortho_cp( eigr, cp, phi, x0, diff, iter, ccc, bephi, becp )
SUBROUTINE ortho_cp( eigr, cp, phi, ngwx, x0, nudx, diff, iter, ccc, &
bephi, becp, nbsp, nspin, nupdwn, iupdwn )
!
! input = cp (non-orthonormal), beta
! input = phi |phi>=s'|c0>
@ -232,7 +233,7 @@
USE cvan, ONLY: ish, nvb
USE uspp, ONLY: nkb, qq
USE uspp_param, ONLY: nh
USE electrons_base, ONLY: nbsp, nbspx, nudx, nspin, nupdwn, iupdwn, f
USE electrons_base, ONLY: f
USE gvecw, ONLY: ngw
USE control_flags, ONLY: iprint, iprsta, ortho_max
USE control_flags, ONLY: force_pairing
@ -240,19 +241,20 @@
!
IMPLICIT NONE
!
COMPLEX(DP) :: cp(ngw,nbsp), phi(ngw,nbsp), eigr(ngw,nat)
INTEGER :: ngwx, nudx, nbsp, nspin
INTEGER :: nupdwn( nspin ), iupdwn( nspin )
COMPLEX(DP) :: cp(ngwx,nbsp), phi(ngwx,nbsp), eigr(ngwx,nat)
REAL(DP) :: x0( nudx, nudx, nspin ), diff, ccc
INTEGER :: iter
REAL(DP) :: bephi(nkb,nbsp), becp(nkb,nbsp)
!
REAL(DP), ALLOCATABLE :: xloc(:,:)
REAL(DP), ALLOCATABLE:: qbephi(:,:), qbecp(:,:)
REAL(DP), ALLOCATABLE :: qbephi(:,:), qbecp(:,:)
INTEGER :: ngwx, nkbx
INTEGER :: nkbx
INTEGER :: istart, nss, ifail, i, j, iss, iv, jv, ia, is, inl, jnl
INTEGER :: nspin_sub
ngwx = ngw
nkbx = nkb
!
! calculation of becp and bephi

6
CPV/ortho_base.f90
View File

@ -920,9 +920,7 @@ CONTAINS
CALL DGEMM( 'N', 'T', nss, nkbus, nss, 1.0d0, x0, nudx, &
bephi( 1, istart ), nkbx, 0.0d0, wtemp, nss )
! CALL MXMA(x0,1,nudx,bephi(1,istart),nkb,1,wtemp,1,nss,nss,nss,nkbus)
!
!
DO i = 1, nss
DO inl = 1, nkbus
bec( inl, i + istart - 1 ) = wtemp( i, inl ) + becp( inl, i + istart - 1 )
@ -953,7 +951,7 @@ CONTAINS
!
IF ( ccc /= 1.0d0 ) THEN
DO j=1,nss
CALL DSCAL(nss,1.0/ccc,x0(1,j),1)
CALL DSCAL(nss,1.0d0/ccc,x0(1,j),1)
END DO
END IF
!

107
CPV/problem_size.f90
View File

@ -14,35 +14,28 @@
CONTAINS
SUBROUTINE cpsizes(nproc)
SUBROUTINE cpsizes()
USE kinds
USE parameters, ONLY: ndmx
use ions_base, only: nsp, nax, nat
use electrons_base, only: nx => nbnd, nspin
use electrons_module, only: n_emp
USE parameters, ONLY: ndmx
use ions_base, only: nat, nsp
use electrons_base, only: nx => nbnd, nspin
use electrons_module, only: n_emp
use reciprocal_vectors, only: ngwx, ngmx, ngmt
use uspp_param, only: nhm, lmaxkb
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE fft_base, ONLY: dfftp, dffts
use gvecb, only: ngb
use uspp_param, only: nhm
use uspp, only: nkb
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE fft_base, ONLY: dfftp, dffts
implicit none
integer NGWXM_EMP,nhm_EMP,NSAX_EMP,NCHAINX, NPROC
integer NNR1X, NNR2X,NNR3X,LM1X
integer NR1X, NR2X, NR3X, nr1_l,nr2_l,nr3_l
integer NSAX,NAMX,NAFX
integer nr1x, nr2x, nr3x, nr1_l, nr2_l, nr3_l
integer nbyte
integer nbyte_alloc
integer itmp
integer nk
integer :: csiz = 1
integer :: lsiz = 4
integer :: isiz = 4
integer :: dsiz = 8
integer :: zsiz = 16
nr1_l = dfftp%nr1x
nr2_l = dfftp%nr2x
nr3_l = dfftp%npl
@ -51,84 +44,44 @@
nr2x = dfftp%nr2x
nr3x = dfftp%nr3x
nk = 1
nhm_EMP = nhm
NGWXM_EMP = ngwx
NNR1X = 2*NR1X-1
NNR2X = 2*NR2X-1
NNR3X = 2*NR3X-1
NSAX = NSP*NAX
NAMX = NSAX
NAFX = NSAX
LM1X = lmaxkb
nbyte = 0
nbyte_alloc = 0
NSAX_EMP = NSAX
NCHAINX = 1
! ... Atoms type
nbyte = nbyte + 3 * ( 4 * isiz + nsp * isiz + 4 * csiz + 2 * ( nsp * isiz ) + nsp * dsiz + &
3 * ( 3 * dsiz * nat ) + 3 * lsiz * nat + isiz * nat + 3 * lsiz + nsp * dsiz + dsiz )
nbyte = nbyte + 8* 3 * 14 * nat
! ... HG_L mill gx_L
nbyte = nbyte + 8 * ngmx
nbyte = nbyte + 3 * 8 * ngmx
nbyte = nbyte + 8 * ngmx * 3
! ... GVEC
nbyte = nbyte + 8 * 10 * ngb
nbyte = nbyte + 8 * 13 * ngmx
! ... FNL
nbyte = nbyte + 8 * NSAX * NX * nhm * nspin
! ... WNL RHOPS VPS GNL RW RPS VR
nbyte = nbyte + 8*ndmx*NSP*3
nbyte = nbyte + 8*ngwx*nhm*NSP
nbyte = nbyte + 8*NSP*ngmx
nbyte = nbyte + 8*NSP*ngmx
nbyte = nbyte + 8*ndmx*NSP
nbyte = nbyte + 8*ndmx*NSP*LM1X
nbyte = nbyte + 8*ndmx*NSP*LM1X
! ... Pseudo
nbyte = nbyte + 8 * 5 * nkb * nx * nspin
! ... C0 CM CP C_EMP
nbyte = nbyte + 16 * ngwx * nx * nspin * nk
nbyte = nbyte + 3 * 16 * ngwx * nx * nspin * nk
nbyte = nbyte + 3 * 16 * ngwx * nx * nspin
nbyte = nbyte + 3 * 16 * ngwx * n_emp * nspin
! ... ei1 ei2 ei3, eigr, sfac
nbyte = nbyte + 3 * 16 * NNR1X * NAX * NSP
nbyte = nbyte + 16 * ngwx * nax * nsp
nbyte = nbyte + 3 * 16 * MAX( nr1x, nr2x, nr3x ) * nat
nbyte = nbyte + 16 * ngwx * nat
nbyte = nbyte + 16 * ngmx * nsp
! ... rhoe and vpot (nr1_l, nr2_l, nr3_l, nspin)
nbyte = nbyte + ( 8 + 16 )*NR1_L*NR2_L*NR3_L*nspin
! ... rhoe and vpot ( nr1_l, nr2_l, nr3_l, nspin )
nbyte = nbyte + ( 8 + 16 ) * NR1_L * NR2_L * NR3_L * nspin
! ... TEMPORARY ALLOCATED MEMORY
itmp = 2 * 8 * ngmt
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
! ... ortho
itmp = 8 * ( 8*NX*NX+2*NX )
itmp = 8 * 8 * NX * NX
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
! ... nlsm1 e nlsm2
itmp = 8 * 2*ngwx*NSAX + 8*NSAX*NX*nhm
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
! ... eigsnew
itmp = 8 * ( 3*NX + NX*NX + N_EMP + NSAX_EMP*N_EMP*nhm_EMP &
& + NX*(NX+1)/2 + N_EMP*N_EMP + N_EMP + NX )
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
! ... phfac
itmp = 8 * 2* ( 2*NR1X + 2*NR2X + 2*NR3X )
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
! ... pvofrho & pstress
itmp = 8 * (2*3*NAX*NSP+3*NAX*NSP+NR1_L*NR2_L*NR3_L*8 + &
& NSAX*NX*nhm*6 + 6*ngmx + ndmx + 6*ngwx + &
& ngwx*nhm*NSP + 2*ngwx*NSAX)
itmp = 8 * ( NR1_L * NR2_L * NR3_L * 8 + &
& nat * NX * nhm * 6 + 6 * ngmx + ndmx + 6*ngwx + &
& ngwx*nhm*nsp + 2*ngwx*nat )
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
! ... formf
itmp = 8 * 3 * ndmx + 8 * 3 * ngmx
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
IF(ionode) THEN
WRITE( stdout,10) nbyte + nbyte_alloc
END IF

9
CPV/restart_sub.f90
View File

@ -289,8 +289,8 @@ MODULE from_restart_module
!
IF ( tortho ) THEN
!
CALL ortho( eigr, cm(:,:,1), phi(:,:,1), lambda, bigr, iter, &
dt2bye, bephi, becp )
CALL ortho( eigr, cm(:,:,1), phi(:,:,1), ngw, lambda, SIZE(lambda,1), &
bigr, iter, dt2bye, bephi, becp, nbsp, nspin, nupdwn, iupdwn )
!
n_spin_start = nspin
IF( force_pairing ) n_spin_start = 1
@ -380,7 +380,7 @@ MODULE from_restart_module
!
lambdam = lambda
!
CALL elec_fakekine( ekincm, ema0bg, emass, c0, cm, ngw, nbsp, delt )
CALL elec_fakekine( ekincm, ema0bg, emass, c0, cm, ngw, nbsp, 1, delt )
!
xnhe0 = 0.D0
xnhem = 0.D0
@ -489,7 +489,6 @@ MODULE from_restart_module
USE phase_factors_module, ONLY : strucf, phfacs
USE time_step, ONLY : delt
USE charge_density, ONLY : rhoofr
USE wave_functions, ONLY : fixwave
USE wave_base, ONLY : wave_verlet, rande_base
USE electrons_module, ONLY : occn_info
USE electrons_base, ONLY : nspin
@ -502,12 +501,10 @@ MODULE from_restart_module
USE cell_base, ONLY : r_to_s, s_to_r
USE nl, ONLY : nlrh_m
USE potentials, ONLY : vofrhos
USE forces, ONLY : dforce_all
USE orthogonalize, ONLY : ortho
USE io_global, ONLY : ionode, ionode_id
USE io_global, ONLY : stdout
USE wave_types, ONLY : wave_descriptor
USE pseudo_projector, ONLY : projector
USE control_flags, ONLY : tcarpar, nbeg, tranp, amprp, tfor, tsdp, &
thdyn, tsdc, tbeg, tsde, tortho, tzeroe, &
tzerop, tzeroc, taurdr, tv0rd, nv0rd, &

79
CPV/runcg.f90
View File

@ -66,7 +66,7 @@
USE electrons_base, ONLY: nupdwn, iupdwn
USE cp_electronic_mass, ONLY: emass
USE cp_main_variables, ONLY: ema0bg
USE wave_functions, ONLY: cp_kinetic_energy, proj, fixwave
USE wave_functions, ONLY: elec_fakekine, proj
USE wave_base, ONLY: dotp, hpsi
USE wave_constrains, ONLY: update_lambda
USE check_stop, ONLY: check_stop_now
@ -116,12 +116,12 @@
REAL(DP), ALLOCATABLE :: dt2bye( : )
REAL(DP) :: gg, ggo, ekinc_old, emin, demin, dek, dt2fact
REAL(DP) :: gg, ggo, ekinc_old, emin, demin, dek, dt2fact, eks
COMPLEX(DP) :: lambda
COMPLEX(DP), ALLOCATABLE :: hacca(:,:,:)
INTEGER :: ib, ibl, ik, ispin, ngw, nfi_l, nspin, isteep, i
INTEGER :: ib, ibl, ik, iss, ngw, nfi_l, nspin, isteep, i
INTEGER :: nk, iter, ierr
LOGICAL :: gamma_symmetry
LOGICAL :: tbad
@ -179,43 +179,43 @@
s2 = cclock()
DO ispin = 1, nspin
DO iss = 1, nspin
! ... Calculate wave functions gradient (temporarely stored in cp)
! ... |d H / dPsi_j > = H |Psi_j> - Sum{i} <Psi_i|H|Psi_j> |Psi_i>
CALL dforce_all( ispin, c0(:,:,ispin), occ(:,ispin), cp(:,:,ispin), &
vpot(:,ispin), eigr, bec, nupdwn, iupdwn )
CALL dforce_all( c0(:,:,iss), occ(:,iss), cp(:,:,iss), &
vpot(:,iss), eigr, bec, nupdwn(iss), iupdwn(iss) )
! ... Project the gradient
CALL proj( ispin, cp(:,:,ispin), cdesc, c0(:,:,ispin), cdesc )
CALL proj( iss, cp(:,:,iss), cdesc, c0(:,:,iss), cdesc )
END DO
s3 = cclock()
! ... Calculate new direction hacca for the line minimization
DO ispin = 1, nspin
DO iss = 1, nspin
DO ik = 1, nk
DO i = 1, nb( ispin )
cp(:, i, ispin) = cp(:, i, ispin) * dt2bye(:) * dt2fact
DO i = 1, nb( iss )
cp(:, i, iss) = cp(:, i, iss) * dt2bye(:) * dt2fact
IF( iter > 1 ) THEN
IF( gamma_symmetry ) THEN
ggo = dotp( gzero, cm(:, i, ispin), cm(:, i, ispin) )
ggo = dotp( gzero, cm(:, i, iss), cm(:, i, iss) )
ELSE
ggo = dotp( cm(:, i, ispin), cm(:, i, ispin) )
ggo = dotp( cm(:, i, iss), cm(:, i, iss) )
END IF
cm(:, i, ispin) = cp(:, i, ispin) - cm(:, i, ispin)
cm(:, i, iss) = cp(:, i, iss) - cm(:, i, iss)
IF( gamma_symmetry ) THEN
gg = dotp( gzero, cm(:, i, ispin), cp(:, i, ispin))
gg = dotp( gzero, cm(:, i, iss), cp(:, i, iss))
ELSE
gg = dotp( cm(:, i, ispin), cp(:, i, ispin))
gg = dotp( cm(:, i, iss), cp(:, i, iss))
END IF
lambda = gg / ggo
hacca(:, i, ispin) = cp(:, i, ispin) + lambda * hacca(:, i, ispin)
hacca(:, i, iss) = cp(:, i, iss) + lambda * hacca(:, i, iss)
ELSE
hacca(:, i, ispin) = cp(:, i, ispin)
hacca(:, i, iss) = cp(:, i, iss)
END IF
END DO
END DO
@ -246,23 +246,24 @@
cp = c0 + cm
DO ispin = 1, nspin
CALL fixwave( ispin, cp(:,:,ispin), cdesc )
DO iss = 1, nspin
IF( gzero ) cp( 1, :, iss ) = DBLE( cp( 1, :, iss ) )
END DO
IF( tortho ) THEN
CALL ortho( c0, cp, cdesc )
ELSE
DO ispin = 1, nspin
CALL gram( vkb, bec, nkb, cp(1,1,ispin), SIZE(cp,1), cdesc%nbt( ispin ) )
DO iss = 1, nspin
CALL gram( vkb, bec, nkb, cp(1,1,iss), SIZE(cp,1), cdesc%nbt( iss ) )
END DO
END IF
END IF
ekinc = 0.0d0
DO ispin = 1, nspin
ekinc = ekinc + cp_kinetic_energy( ispin, cp(:,:,ispin), c0(:,:,ispin), cdesc, ema0bg, emass, delt)
DO iss = 1, nspin
CALL elec_fakekine( eks, ema0bg, emass, cp(:,:,iss), c0(:,:,iss), ngw, nb( iss ), 1, delt )
ekinc = ekinc + eks
END DO
IF( iter > 1 ) THEN
dek = ekinc - ekinc_old
@ -314,28 +315,28 @@
END IF
IF( tprint ) THEN
DO ispin = 1, nspin
DO iss = 1, nspin
CALL dforce_all( ispin, c0(:,:,ispin), occ(:,ispin), hacca(:,:,ispin), &
vpot(:,ispin), eigr, bec, nupdwn, iupdwn )
CALL dforce_all( c0(:,:,iss), occ(:,iss), hacca(:,:,iss), &
vpot(:,iss), eigr, bec, nupdwn(iss), iupdwn(iss) )
nb_g( ispin ) = cdesc%nbt( ispin )
nb_g( iss ) = cdesc%nbt( iss )
IF( gamma_symmetry ) THEN
ALLOCATE(cgam(1,1), gam( nb_l( ispin ), nb_g( ispin ) ), STAT=ierr)
ALLOCATE(cgam(1,1), gam( nb_l( iss ), nb_g( iss ) ), STAT=ierr)
ELSE
ALLOCATE(cgam(nb_l( ispin ),nb_g( ispin )), gam(1,1), STAT=ierr)
ALLOCATE(cgam(nb_l( iss ),nb_g( iss )), gam(1,1), STAT=ierr)
END IF
IF( ierr/=0 ) CALL errore(' runcg ', ' allocating gam ',ierr)
DO ik = 1, nk
DO i = 1, nb( ispin )
DO i = 1, nb( iss )
IF( gamma_symmetry ) THEN
CALL update_lambda( i, gam, c0(:,:,ispin), cdesc, hacca(:,i,ispin) )
CALL update_lambda( i, gam, c0(:,:,iss), cdesc, hacca(:,i,iss) )
ELSE
CALL update_lambda( i, cgam, c0(:,:,ispin), cdesc, hacca(:,i,ispin) )
CALL update_lambda( i, cgam, c0(:,:,iss), cdesc, hacca(:,i,iss) )
END IF
END DO
CALL eigs( nb( ispin ), gam, tortho, occ(:,ispin), ei(:,ispin) )
CALL eigs( nb( iss ), gam, tortho, occ(:,iss), ei(:,iss) )
END DO
DEALLOCATE( cgam, gam, STAT=ierr )
IF( ierr/=0 ) CALL errore(' runcg ', ' deallocating gam ',ierr)
@ -364,7 +365,6 @@
USE wave_types
USE energies, ONLY: dft_energy_type
USE wave_functions, ONLY: fixwave
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE cell_module, ONLY: boxdimensions
@ -408,7 +408,7 @@
REAL(DP) :: x, p, v, w, e, fw, fv, xm, tol1, tol2, a, b, etemp, d
REAL(DP) :: fx, xmin, brent, eold, tol
LOGICAL :: tbrent
INTEGER :: iter, ispin
INTEGER :: iter, iss
!
! ... SUBROUTINE BODY
@ -589,9 +589,12 @@
cp = c + hstep * hacca
DO ispin = 1, cdesc%nspin
CALL fixwave( ispin, cp(:,:,ispin), cdesc )
CALL gram( vkb, bec, nkb, cp(1,1,ispin), SIZE(cp,1), cdesc%nbt( ispin ) )
DO iss = 1, cdesc%nspin
IF( cdesc%gzero ) cp( 1, :, iss ) = DBLE( cp( 1, :, iss ) )
END DO
!
DO iss = 1, cdesc%nspin
CALL gram( vkb, bec, nkb, cp(1,1,iss), SIZE(cp,1), cdesc%nbt( iss ) )
END DO

60
CPV/runcp.f90
View File

@ -44,13 +44,12 @@
USE electrons_module, ONLY: eigs, nb_l
USE cp_electronic_mass, ONLY: emass
USE cp_main_variables, ONLY: ema0bg
USE wave_functions, ONLY : cp_kinetic_energy
USE wave_functions, ONLY : elec_fakekine
USE wave_base, ONLY: hpsi
USE cell_module, ONLY: boxdimensions
USE time_step, ONLY: delt
USE orthogonalize, ONLY: ortho
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
USE wave_constrains, ONLY: update_lambda
USE uspp, ONLY : vkb, nkb
@ -73,7 +72,7 @@
! ... declare other variables
REAL(DP) :: s1, s2, s3, s4
INTEGER :: nx, nb_lx, ierr, is
INTEGER :: nx, nb_lx, ierr, is, nb, ngw
COMPLEX(DP), ALLOCATABLE :: cgam(:,:,:)
REAL(DP), ALLOCATABLE :: gam(:,:,:)
@ -85,8 +84,10 @@
s1 = cclock()
ngw = SIZE( cp, 1 )
nb_lx = MAX( nb_l(1), nb_l(2) )
nb_lx = MAX( nb_lx, 1 )
!
IF( cdesc%gamma ) THEN
ALLOCATE( cgam(1,1,1), gam( nb_lx, SIZE( c0, 2 ), cdesc%nspin ), STAT=ierr)
ELSE
@ -131,7 +132,8 @@
! Compute fictitious kinetic energy of the electrons at time t
DO is = 1, cdesc%nspin
ekinc(is) = cp_kinetic_energy( is, cp(:,:,is), cm(:,:,is), cdesc, ema0bg, emass, delt)
nb = cdesc%nbt( is )
CALL elec_fakekine( ekinc( is ), ema0bg, emass, cp(:,:,is), cm(:,:,is), ngw, nb, 1, 2.0d0 * delt )
END DO
DEALLOCATE( cgam, gam, STAT=ierr)
@ -171,7 +173,6 @@
USE wave_types, ONLY: wave_descriptor
USE wave_constrains, ONLY: update_lambda
USE control_flags, ONLY: tsde
USE pseudo_projector, ONLY: projector
IMPLICIT NONE
@ -247,7 +248,7 @@
DO i = 1, nb, 2
CALL dforce( i, is, c0(:,:,is), fi(:,is), c2, c3, vpot(:,is), eigr, bec, nupdwn, iupdwn )
CALL dforce( i, c0(:,:,is), fi(:,is), c2, c3, vpot(:,is), eigr, bec, nupdwn(is), iupdwn(is) )
IF( tlam ) THEN
CALL update_lambda( i, gam( :, :,is), c0(:,:,is), cdesc, c2 )
@ -278,7 +279,7 @@
nb = nx
CALL dforce( nx, is, c0(:,:,is), fi(:,is), c2, c3, vpot(:,is), eigr, bec, nupdwn, iupdwn )
CALL dforce( nx, c0(:,:,is), fi(:,is), c2, c3, vpot(:,is), eigr, bec, nupdwn(is), iupdwn(is) )
IF( tlam ) THEN
CALL update_lambda( nb, gam( :, :,is), c0(:,:,is), cdesc, c2 )
@ -331,7 +332,7 @@
USE electrons_base, ONLY: iupdwn, nupdwn, nspin
USE cp_electronic_mass, ONLY: emass
USE cp_main_variables, ONLY: ema0bg
USE wave_functions, ONLY : cp_kinetic_energy
USE wave_functions, ONLY : elec_fakekine
USE wave_base, ONLY: wave_steepest, wave_verlet
USE wave_base, ONLY: hpsi
USE cell_module, ONLY: boxdimensions
@ -443,8 +444,8 @@
DO i = 1, nb, 2
!
CALL dforce( i, 2, c0(:,:,1), fi(:,1), c2, c3, vpot(:,1), eigr, bec, nupdwn, iupdwn )
CALL dforce( i, 2, c0(:,:,1), fi(:,1), c4, c5, vpot(:,2), eigr, bec, nupdwn, iupdwn )
CALL dforce( i, c0(:,:,1), fi(:,1), c2, c3, vpot(:,1), eigr, bec, nupdwn(2), iupdwn(2) )
CALL dforce( i, c0(:,:,1), fi(:,1), c4, c5, vpot(:,2), eigr, bec, nupdwn(2), iupdwn(2) )
!
c2 = occup(i )* (c2 + c4)
c3 = occup(i+1)* (c3 + c5)
@ -476,8 +477,8 @@
!
nb = n_unp - 1
!
CALL dforce( nb, 2, c0(:,:,1), fi(:,1), c2, c3, vpot(:,1), eigr, bec, nupdwn, iupdwn )
CALL dforce( nb, 2, c0(:,:,1), fi(:,2), c4, c5, vpot(:,2), eigr, bec, nupdwn, iupdwn )
CALL dforce( nb, c0(:,:,1), fi(:,1), c2, c3, vpot(:,1), eigr, bec, nupdwn(2), iupdwn(2) )
CALL dforce( nb, c0(:,:,1), fi(:,2), c4, c5, vpot(:,2), eigr, bec, nupdwn(2), iupdwn(2) )
c2 = occup(nb)* (c2 + c4)
@ -495,7 +496,7 @@
END IF
!
CALL dforce( n_unp, 1, c0(:,:,1), fi(:,1), c2, c3, vpot(:,1), eigr, bec, nupdwn, iupdwn )
CALL dforce( n_unp, c0(:,:,1), fi(:,1), c2, c3, vpot(:,1), eigr, bec, nupdwn(1), iupdwn(1) )
intermed = -2.d0 * sum( c2 * conjg( c0(:, n_unp, 1 ) ) )
IF ( gstart == 2 ) THEN
@ -563,8 +564,8 @@
! Compute fictitious kinetic energy of the electrons at time t
ekinc(1) = cp_kinetic_energy( 1, cp(:,:,1), cm(:,:,1), cdesc, ema0bg, emass, delt)
ekinc(2) = cp_kinetic_energy( 2, cp(:,:,1), cm(:,:,1), cdesc, ema0bg, emass, delt)
CALL elec_fakekine( ekinc( 1 ), ema0bg, emass, cp(:,:,1), cm(:,:,1), ngw, cdesc%nbt(1), 1, 2.0d0 * delt )
CALL elec_fakekine( ekinc( 2 ), ema0bg, emass, cp(:,:,1), cm(:,:,1), ngw, cdesc%nbt(2), 1, 2.0d0 * delt )
DEALLOCATE( ei_t, svar3, c2, c3, c4, c5, cgam, gam, occup, occdown, STAT=ierr)
@ -585,7 +586,7 @@
!use uspp, only : nhsa=> nkb, betae => vkb, rhovan => becsum, deeq
use uspp, only : deeq, betae => vkb
use reciprocal_vectors, only : gstart
use electrons_base, only : n=>nbsp
use electrons_base, only : n=>nbsp, ispin, f, nspin
use wannier_subroutines, only: ef_potential
use efield_module, only: dforce_efield, tefield, dforce_efield2, tefield2
@ -645,7 +646,7 @@
call ef_potential( nfi, rhos, bec, deeq, betae, c0, cm, emadt2, emaver, verl1, verl2, c2, c3 )
else
do i=1,n,2
call dforce(bec,betae,i,c0(1,i),c0(1,i+1),c2,c3,rhos)
call dforce(bec,betae,i,c0(1,i),c0(1,i+1),c2,c3,rhos,ispin,f,n,nspin)
if( tefield ) then
CALL dforce_efield ( bec, i, c0, c2, c3, rhos)
end if
@ -855,13 +856,14 @@
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: lwf, tsde
! use uspp, only : nhsa=> nkb, betae => vkb, rhovan => becsum, deeq
USE uspp, ONLY : deeq, betae => vkb
USE reciprocal_vectors, ONLY : gstart
USE wannier_subroutines, ONLY: ef_potential
USE efield_module, ONLY: dforce_efield, tefield
USE wave_base, ONLY : wave_steepest, wave_verlet
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
USE wannier_subroutines, ONLY : ef_potential
USE efield_module, ONLY : dforce_efield, tefield
USE electrons_base, ONLY : ispin, nspin, f, n=>nbsp
!
USE gvecw, ONLY: ngw
!
@ -958,8 +960,8 @@
DO i = 1, npair, 2
!
CALL dforce(bec,betae,i,c0(1,i),c0(1,i+1),c2,c3,rhos(1,1))
CALL dforce(bec,betae,i,c0(1,i),c0(1,i+1),c4,c5,rhos(1,2))
CALL dforce(bec,betae,i,c0(1,i),c0(1,i+1),c2,c3,rhos(1,1),ispin,f,n,nspin)
CALL dforce(bec,betae,i,c0(1,i),c0(1,i+1),c4,c5,rhos(1,2),ispin,f,n,nspin)
!
c2 = occ( i )*(c2 + c4)
c3 = occ(i+1)*(c3 + c5)
@ -988,8 +990,8 @@
npair = n_unp - 1
!
CALL dforce(bec,betae,npair,c0(1,npair),c0(1,nbspx),c2,c3,rhos(1,1))
CALL dforce(bec,betae,npair,c0(1,npair),c0(1,nbspx),c4,c5,rhos(1,2))
CALL dforce(bec,betae,npair,c0(1,npair),c0(1,nbspx),c2,c3,rhos(1,1),ispin,f,n,nspin)
CALL dforce(bec,betae,npair,c0(1,npair),c0(1,nbspx),c4,c5,rhos(1,2),ispin,f,n,nspin)
!
c2 = c2 + c4
!
@ -1018,7 +1020,7 @@
! "TRUE" ONLY WHEN THE POT is NORM_CONSERVING
!
CALL dforce( bec, betae, n_unp, c0(1,n_unp), c0(1,n_unp), c2, c3, rhos(1,1) )
CALL dforce( bec, betae, n_unp, c0(1,n_unp), c0(1,n_unp), c2, c3, rhos(1,1),ispin,f,n,nspin )
!
intermed = - 2.d0 * sum(c2 * conjg(c0(:,n_unp)))
IF ( gstart == 2 ) THEN

14
CPV/rundiis.f90
View File

@ -266,7 +266,7 @@ CONTAINS
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, bec, becdr, eigr)
CALL dforce_all( 1, c0(:,:,1), fi(:,1), cgrad(:,:,1), vpot(:,1), eigr, bec, nupdwn, iupdwn )
CALL dforce_all( c0(:,:,1), fi(:,1), cgrad(:,:,1), vpot(:,1), eigr, bec, nupdwn(1), iupdwn(1) )
CALL proj( 1, cgrad(:,:,1), cdesc, c0(:,:,1), cdesc, lambda )
CALL crot( 1, c0(:,:,1), cdesc, lambda, eig(:,1) )
@ -275,7 +275,7 @@ CONTAINS
call entropy_s(fi(1,1),temp_elec,cdesc%nbl(1),edft%ent)
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, bec, becdr, eigr)
CALL dforce_all( 1, c0(:,:,1), fi(:,1), cgrad(:,:,1), vpot(:,1), eigr, bec, nupdwn, iupdwn )
CALL dforce_all( c0(:,:,1), fi(:,1), cgrad(:,:,1), vpot(:,1), eigr, bec, nupdwn(1), iupdwn(1) )
DO ib = 1, cdesc%nbl( 1 )
cgrad(:,ib,1) = cgrad(:,ib,1) + eig(ib,1)*c0(:,ib,1)
@ -286,7 +286,7 @@ CONTAINS
! ... DIIS on c0 at FIXED potential
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, bec, becdr, eigr)
CALL dforce_all( 1, c0(:,:,1), fi(:,1), cgrad(:,:,1), vpot(:,1), eigr, bec, nupdwn, iupdwn )
CALL dforce_all( c0(:,:,1), fi(:,1), cgrad(:,:,1), vpot(:,1), eigr, bec, nupdwn(1), iupdwn(1) )
CALL proj( 1, cgrad(:,:,1), cdesc, c0(:,:,1), cdesc, lambda)
@ -533,8 +533,8 @@ CONTAINS
! ... of 1, ( row 1 to PE 1, row 2 to PE 2, .. row nproc_image+1 to PE 1 and
! ... so on).
CALL dforce_all( ispin, c0(:,:,ispin), fi(:,ispin), cgrad(:,:,ispin), &
vpot(:,ispin), eigr, bec, nupdwn, iupdwn )
CALL dforce_all( c0(:,:,ispin), fi(:,ispin), cgrad(:,:,ispin), &
vpot(:,ispin), eigr, bec, nupdwn(ispin), iupdwn(ispin) )
CALL proj( ispin, cgrad(:,:,ispin), cdesc, c0(:,:,ispin), cdesc, lambda)
@ -657,8 +657,8 @@ CONTAINS
CALL nlsm1( n, 1, nspnl, eigr, c(1,1,ispin), bec )
! ... Calculate | dH / dpsi(j) >
CALL dforce_all( ispin, c(:,:,ispin), fi(:,ispin), eforce(:,:,ispin), &
vpot(:,ispin), eigr, bec, nupdwn, iupdwn )
CALL dforce_all( c(:,:,ispin), fi(:,ispin), eforce(:,:,ispin), &
vpot(:,ispin), eigr, bec, nupdwn(ispin), iupdwn(ispin) )
! ... Calculate Eij = < psi(i) | H | psi(j) > = < psi(i) | dH / dpsi(j) >
DO i = 1, n

20
CPV/smcp.f90
View File

@ -29,6 +29,7 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
esr, print_energies
USE electrons_base, ONLY : nbspx, nbsp, f, nspin, nudx
USE electrons_base, ONLY : nel, iupdwn, nupdwn
USE electrons_module, ONLY : ei
USE gvecp, ONLY : ngm
USE gvecs, ONLY : ngs
USE gvecb, ONLY : ngb
@ -93,7 +94,7 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
USE orthogonalize, ONLY : ortho
USE orthogonalize_base, ONLY : updatc, calphi
use charge_density, only : rhoofr
USE wave_functions, ONLY : wave_rand_init, elec_fakekine2
USE wave_functions, ONLY : wave_rand_init, elec_fakekine
!
#if ! defined __NOSMD
!
@ -611,8 +612,9 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
! ==========================================================
!
IF(tortho) THEN
CALL ortho (eigr,rep_el(sm_k)%c0,rep_el(sm_k)%phi,rep_el(sm_k)%lambda, &
& bigr,iter,ccc(sm_k),bephi,becp)
CALL ortho ( eigr, rep_el(sm_k)%c0, rep_el(sm_k)%phi, ngw, &
rep_el(sm_k)%lambda, SIZE( rep_el(sm_k)%lambda, 1 ), &
bigr, iter, ccc(sm_k), bephi, becp, nbsp, nspin, nupdwn, iupdwn)
ELSE
CALL gram( vkb, rep_el(sm_k)%bec, nkb, rep_el(sm_k)%c0, ngw, nbsp )
!
@ -780,7 +782,7 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
! ======================================================
!
ekincm(sm_k)=0.0
CALL elec_fakekine2( ekincm(sm_k), ema0bg, emass, rep_el(sm_k)%c0, rep_el(sm_k)%cm, ngw, nbsp, delt )
CALL elec_fakekine( ekincm(sm_k), ema0bg, emass, rep_el(sm_k)%c0, rep_el(sm_k)%cm, ngw, nbsp, 1, delt )
xnhe0(sm_k)=0.
xnhem(sm_k)=0.
@ -1138,9 +1140,9 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
!---------------------------------------------------------------------------
!
IF(tortho) THEN
CALL ortho &
& (eigr,rep_el(sm_k)%cm,rep_el(sm_k)%phi,rep_el(sm_k)%lambda, &
& bigr,iter,ccc(sm_k),bephi,becp)
CALL ortho ( eigr, rep_el(sm_k)%cm, rep_el(sm_k)%phi, ngw, &
rep_el(sm_k)%lambda, SIZE( rep_el(sm_k)%lambda, 1 ), &
bigr, iter, ccc(sm_k), bephi, becp, nbsp, nspin, nupdwn, iupdwn )
ELSE
CALL gram( vkb, rep_el(sm_k)%bec, nkb, rep_el(sm_k)%cm, ngw, nbsp )
IF(iprsta.GT.4) CALL dotcsc(eigr,rep_el(sm_k)%cm, ngw, nbsp )
@ -1189,7 +1191,7 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
!
! fake electronic kinetic energy
!
CALL elec_fakekine2( ekinc0(sm_k), ema0bg, emass, rep_el(sm_k)%c0, rep_el(sm_k)%cm, ngw, nbsp, delt )
CALL elec_fakekine( ekinc0(sm_k), ema0bg, emass, rep_el(sm_k)%c0, rep_el(sm_k)%cm, ngw, nbsp, 1, delt )
!
! ... previous ekinc
!
@ -1228,7 +1230,7 @@ SUBROUTINE smdmain( tau, fion_out, etot_out, nat_out )
!
!
IF(MOD(nfi-1,iprint).EQ.0 .OR. (nfi.EQ.(nomore))) THEN
CALL eigs0(.true.,nspin,nupdwn,iupdwn,.true.,f,nbspx,rep_el(sm_k)%lambda,nudx)
CALL eigs0( ei, .true.,nspin,nupdwn,iupdwn,.true.,f,nbspx,rep_el(sm_k)%lambda,nudx)
WRITE( stdout,*)
ENDIF
!

2
CPV/stress.f90
View File

@ -69,7 +69,6 @@
USE mp_global, ONLY: intra_image_comm
USE mp, ONLY: mp_sum
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
USE cell_base, ONLY: tpiba2
USE io_global, ONLY: ionode
USE exchange_correlation, ONLY: stress_xc
@ -276,7 +275,6 @@
USE mp_global, ONLY: intra_image_comm, root_image, me_image
USE io_global, ONLY: stdout
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
USE cell_base, ONLY: tpiba2, omega
USE control_flags, ONLY: force_pairing
USE reciprocal_vectors, ONLY: gstart, gzero, g, gx

4
CPV/wannier.f90
View File

@ -457,7 +457,7 @@ MODULE wannier_subroutines
USE electric_field_module, ONLY : field_tune, e_tuned, par, rel1, rel2
USE wannier_module, ONLY : rhos1, rhos2, wfc
USE smooth_grid_dimensions, ONLY : nnrsx
USE electrons_base, ONLY : nbsp, nspin, nupdwn
USE electrons_base, ONLY : nbsp, nspin, nupdwn, f, ispin
USE cell_base, ONLY : ainv, a1, a2, a3
USE reciprocal_vectors, ONLY : gstart
USE control_flags, ONLY : tsde
@ -555,7 +555,7 @@ MODULE wannier_subroutines
END DO
CALL dforce_field(bec,deeq,betae,i,c0(1,i),c0(1,i+1),c2,c3,rhos1,rhos2)
ELSE
CALL dforce(bec,betae,i,c0(1,i),c0(1,i+1),c2,c3,rhos)
CALL dforce(bec,betae,i,c0(1,i),c0(1,i+1),c2,c3,rhos,ispin,f,nbsp,nspin)
END IF
IF(tsde) THEN
CALL wave_steepest( cm(:, i ), c0(:, i ), emadt2, c2 )

400
CPV/wave.f90
View File

@ -8,6 +8,99 @@
#include "f_defs.h"
!=----------------------------------------------------------------------------=!
MODULE wave_constrains
!=----------------------------------------------------------------------------=!
! ... include modules
USE kinds
IMPLICIT NONE
SAVE
PRIVATE
PUBLIC :: interpolate_lambda, update_lambda
INTERFACE update_lambda
MODULE PROCEDURE update_rlambda, update_clambda
END INTERFACE
!=----------------------------------------------------------------------------=!
CONTAINS
!=----------------------------------------------------------------------------=!
SUBROUTINE interpolate_lambda( lambdap, lambda, lambdam )
IMPLICIT NONE
REAL(DP) :: lambdap(:,:,:), lambda(:,:,:), lambdam(:,:,:)
!
! interpolate new lambda at (t+dt) from lambda(t) and lambda(t-dt):
!
lambdap= 2.d0*lambda - lambdam
lambdam=lambda
lambda =lambdap
RETURN
END SUBROUTINE interpolate_lambda
SUBROUTINE update_rlambda( i, lambda, c0, cdesc, c2 )
USE electrons_module, ONLY: ib_owner, ib_local
USE mp_global, ONLY: me_image, intra_image_comm
USE mp, ONLY: mp_sum
USE wave_base, ONLY: hpsi
USE wave_types, ONLY: wave_descriptor
IMPLICIT NONE
REAL(DP) :: lambda(:,:)
COMPLEX(DP) :: c0(:,:), c2(:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
INTEGER :: i
!
REAL(DP), ALLOCATABLE :: prod(:)
INTEGER :: ibl
!
ALLOCATE( prod( SIZE( c0, 2 ) ) )
prod = hpsi( cdesc%gzero, c0(:,:), c2 )
CALL mp_sum( prod, intra_image_comm )
IF( me_image == ib_owner( i ) ) THEN
ibl = ib_local( i )
lambda( ibl, : ) = prod( : )
END IF
DEALLOCATE( prod )
RETURN
END SUBROUTINE update_rlambda
SUBROUTINE update_clambda( i, lambda, c0, cdesc, c2 )
USE electrons_module, ONLY: ib_owner, ib_local
USE mp_global, ONLY: me_image, intra_image_comm
USE mp, ONLY: mp_sum
USE wave_base, ONLY: hpsi
USE wave_types, ONLY: wave_descriptor
IMPLICIT NONE
COMPLEX(DP) :: lambda(:,:)
COMPLEX(DP) :: c0(:,:), c2(:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
INTEGER :: i
!
COMPLEX(DP), ALLOCATABLE :: prod(:)
INTEGER :: ibl
!
ALLOCATE( prod( SIZE( c0, 2 ) ) )
prod = hpsi( cdesc%gzero, c0(:,:), c2 )
CALL mp_sum( prod, intra_image_comm )
IF( me_image == ib_owner( i ) ) THEN
ibl = ib_local( i )
lambda( ibl, : ) = prod( : )
END IF
DEALLOCATE( prod )
RETURN
END SUBROUTINE update_clambda
!=----------------------------------------------------------------------------=!
END MODULE wave_constrains
!=----------------------------------------------------------------------------=!
!=----------------------------------------------------------------------------=!
MODULE wave_functions
!=----------------------------------------------------------------------------=!
@ -20,7 +113,7 @@
PRIVATE
PUBLIC :: crot, proj, fixwave
PUBLIC :: crot, proj
INTERFACE crot
MODULE PROCEDURE crot_gamma
END INTERFACE
@ -28,50 +121,34 @@
MODULE PROCEDURE proj_gamma, proj2
END INTERFACE
PUBLIC :: cp_kinetic_energy, elec_fakekine2, elec_fakekine
PUBLIC :: update_wave_functions, wave_rand_init
PUBLIC :: elec_fakekine
PUBLIC :: update_wave_functions, wave_rand_init, kohn_sham
!=----------------------------------------------------------------------------=!
CONTAINS
!=----------------------------------------------------------------------------=!
SUBROUTINE fixwave ( ispin, c, cdesc )
USE wave_types, ONLY: wave_descriptor
IMPLICIT NONE
COMPLEX(DP), INTENT(INOUT) :: c(:,:)
INTEGER, INTENT(IN) :: ispin
TYPE (wave_descriptor), INTENT(IN) :: cdesc
!
INTEGER :: i
IF( cdesc%gzero ) THEN
DO i = 1, cdesc%nbl( ispin )
c( 1, i ) = DBLE( c( 1, i ) )
END DO
END IF
RETURN
END SUBROUTINE fixwave
!=----------------------------------------------------------------------------=!
subroutine elec_fakekine( ekincm, ema0bg, emass, c0, cm, ngw, n, delt )
use mp, only: mp_sum
use mp_global, only: intra_image_comm
use reciprocal_vectors, only: gstart
use wave_base, only: wave_speed2
real(8), intent(out) :: ekincm
real(8), intent(in) :: ema0bg(:), delt, emass
complex(8), intent(in) :: c0(:,:,:), cm(:,:,:)
integer, intent(in) :: ngw, n
real(8), allocatable :: emainv(:)
real(8) :: ftmp
integer :: i
subroutine elec_fakekine( ekincm, ema0bg, emass, c0, cm, ngw, n, noff, delt )
!
! This subroutine computes the CP(fake) wave functions kinetic energy
use mp, only : mp_sum
use mp_global, only : intra_image_comm
use reciprocal_vectors, only : gstart
use wave_base, only : wave_speed2
!
integer, intent(in) :: ngw ! number of plane wave coeff.
integer, intent(in) :: n ! number of bands
integer, intent(in) :: noff ! offset for band index
real(DP), intent(out) :: ekincm
real(DP), intent(in) :: ema0bg( ngw ), delt, emass
complex(DP), intent(in) :: c0( ngw, n ), cm( ngw, n )
!
real(DP), allocatable :: emainv(:)
real(DP) :: ftmp
integer :: i
ALLOCATE( emainv( ngw ) )
emainv = 1.0d0 / ema0bg
@ -79,9 +156,8 @@
if( gstart == 2 ) ftmp = 0.5d0
ekincm=0.0d0
do i=1,n
ekincm = ekincm + 2.0d0 * &
wave_speed2( c0(:,i,1), cm(:,i,1), emainv, ftmp )
do i = noff, n + noff - 1
ekincm = ekincm + 2.0d0 * wave_speed2( c0(:,i), cm(:,i), emainv, ftmp )
end do
ekincm = ekincm * emass / ( delt * delt )
@ -92,95 +168,6 @@
end subroutine elec_fakekine
!=----------------------------------------------------------------------------=!
subroutine elec_fakekine2( ekincm, ema0bg, emass, c0, cm, ngw, n, delt )
use mp, only: mp_sum
use mp_global, only: intra_image_comm
use reciprocal_vectors, only: gstart
use wave_base, only: wave_speed2
real(8), intent(out) :: ekincm
real(8), intent(in) :: ema0bg(:), delt, emass
complex(8), intent(in) :: c0(:,:), cm(:,:)
integer, intent(in) :: ngw, n
real(8), allocatable :: emainv(:)
real(8) :: ftmp
integer :: i
ALLOCATE( emainv( ngw ) )
emainv = 1.0d0 / ema0bg
ftmp = 1.0d0
if( gstart == 2 ) ftmp = 0.5d0
ekincm=0.0d0
do i=1,n
ekincm = ekincm + 2.0d0 * &
wave_speed2( c0(:,i), cm(:,i), emainv, ftmp )
end do
ekincm = ekincm * emass / ( delt * delt )
CALL mp_sum( ekincm, intra_image_comm )
DEALLOCATE( emainv )
return
end subroutine elec_fakekine2
!=----------------------------------------------------------------------------=!
REAL(DP) FUNCTION cp_kinetic_energy( ispin, cp, cm, cdesc, ema0bg, emass, delt)
! (describe briefly what this routine does...)
! if ekinc_fp will hold the full electron kinetic energy (paired and unpaired) and
! the function returns the paired electrons' kinetic energy only
! ----------------------------------------------
! ... declare modules
USE mp, ONLY: mp_sum
USE mp_global, ONLY: intra_image_comm
USE wave_types, ONLY: wave_descriptor
USE wave_base, ONLY: wave_speed2
IMPLICIT NONE
! ... declare subroutine arguments
COMPLEX(DP), INTENT(IN) :: cp(:,:), cm(:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
INTEGER, INTENT( IN ) :: ispin
REAL(DP), INTENT(IN) :: delt
REAL(DP), INTENT(IN) :: ema0bg(:), emass
! ... declare other variables
REAL(DP) ekinc, dt2, fact
INTEGER ib, j
REAL(DP), ALLOCATABLE :: emainv(:)
! ... end of declarations
! ----------------------------------------------
dt2 = delt * delt
ekinc = 0.d0
fact = 1.0d0
IF( cdesc%gzero ) fact = 0.5d0
ALLOCATE( emainv( SIZE( cp, 1 ) ) )
emainv = 1.0d0 / ema0bg
DO ib = 1, cdesc%nbl( ispin )
ekinc = ekinc + wave_speed2( cp(:,ib), cm(:,ib), emainv, fact )
END DO
IF( cdesc%gamma ) ekinc = ekinc * 2.0d0
CALL mp_sum( ekinc, intra_image_comm )
cp_kinetic_energy = ekinc * emass / (4.0d0 * dt2)
DEALLOCATE( emainv )
RETURN
END FUNCTION cp_kinetic_energy
!=----------------------------------------------------------------------------=!
SUBROUTINE update_wave_functions(cm, c0, cp, cdesc)
@ -196,16 +183,14 @@
COMPLEX(DP), INTENT(OUT) :: cm(:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
INTEGER :: ispin, ik, nspin
INTEGER :: ispin, nspin
nspin = cdesc%nspin
IF( force_pairing ) nspin = 1
DO ispin = 1, nspin
DO ik = 1, cdesc%nkl
cm(:,:,ispin) = c0(:,:,ispin)
c0(:,:,ispin) = cp(:,:,ispin)
END DO
cm(:,:,ispin) = c0(:,:,ispin)
c0(:,:,ispin) = cp(:,:,ispin)
END DO
RETURN
@ -516,100 +501,59 @@
END SUBROUTINE wave_rand_init
SUBROUTINE kohn_sham(ispin, c, cdesc, eforces, nupdwn, nupdwnl )
! ... declare modules
USE kinds
USE wave_constrains, ONLY: update_lambda
USE wave_types, ONLY: wave_descriptor
IMPLICIT NONE
! ... declare subroutine arguments
COMPLEX(DP), INTENT(INOUT) :: c(:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
INTEGER, INTENT(IN) :: ispin
INTEGER, INTENT(IN) :: nupdwn(:) ! number of upper and down states
INTEGER, INTENT(IN) :: nupdwnl(:) ! local (to the processor) number of up and down states
COMPLEX(DP) :: eforces(:,:)
! ... declare other variables
INTEGER :: ib, nb_g, nrl
REAL(DP), ALLOCATABLE :: gam(:,:)
REAL(DP), ALLOCATABLE :: eig(:)
LOGICAL :: tortho = .TRUE.
! ... end of declarations
nb_g = nupdwn( ispin )
nrl = nupdwnl( ispin )
IF( nb_g < 1 ) THEN
eforces = 0.0d0
ELSE
ALLOCATE( eig( nb_g ) )
ALLOCATE( gam( nrl, nb_g ) )
DO ib = 1, nb_g
CALL update_lambda( ib, gam, c(:,:), cdesc, eforces(:,ib) )
END DO
CALL crot( ispin, c(:,:), cdesc, gam, eig )
DEALLOCATE( gam, eig )
END IF
RETURN
! ...
END SUBROUTINE kohn_sham
!=----------------------------------------------------------------------------=!
END MODULE wave_functions
!=----------------------------------------------------------------------------=!
!=----------------------------------------------------------------------------=!
MODULE wave_constrains
!=----------------------------------------------------------------------------=!
! ... include modules
USE kinds
IMPLICIT NONE
SAVE
PRIVATE
PUBLIC :: interpolate_lambda, update_lambda
INTERFACE update_lambda
MODULE PROCEDURE update_rlambda, update_clambda
END INTERFACE
!=----------------------------------------------------------------------------=!
CONTAINS
!=----------------------------------------------------------------------------=!
SUBROUTINE interpolate_lambda( lambdap, lambda, lambdam )
IMPLICIT NONE
REAL(DP) :: lambdap(:,:,:), lambda(:,:,:), lambdam(:,:,:)
!
! interpolate new lambda at (t+dt) from lambda(t) and lambda(t-dt):
!
lambdap= 2.d0*lambda - lambdam
lambdam=lambda
lambda =lambdap
RETURN
END SUBROUTINE interpolate_lambda
SUBROUTINE update_rlambda( i, lambda, c0, cdesc, c2 )
USE electrons_module, ONLY: ib_owner, ib_local
USE mp_global, ONLY: me_image, intra_image_comm
USE mp, ONLY: mp_sum
USE wave_base, ONLY: hpsi
USE wave_types, ONLY: wave_descriptor
IMPLICIT NONE
REAL(DP) :: lambda(:,:)
COMPLEX(DP) :: c0(:,:), c2(:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
INTEGER :: i
!
REAL(DP), ALLOCATABLE :: prod(:)
INTEGER :: ibl
!
ALLOCATE( prod( SIZE( c0, 2 ) ) )
prod = hpsi( cdesc%gzero, c0(:,:), c2 )
CALL mp_sum( prod, intra_image_comm )
IF( me_image == ib_owner( i ) ) THEN
ibl = ib_local( i )
lambda( ibl, : ) = prod( : )
END IF
DEALLOCATE( prod )
RETURN
END SUBROUTINE update_rlambda
SUBROUTINE update_clambda( i, lambda, c0, cdesc, c2 )
USE electrons_module, ONLY: ib_owner, ib_local
USE mp_global, ONLY: me_image, intra_image_comm
USE mp, ONLY: mp_sum
USE wave_base, ONLY: hpsi
USE wave_types, ONLY: wave_descriptor
IMPLICIT NONE
COMPLEX(DP) :: lambda(:,:)
COMPLEX(DP) :: c0(:,:), c2(:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
INTEGER :: i
!
COMPLEX(DP), ALLOCATABLE :: prod(:)
INTEGER :: ibl
!
ALLOCATE( prod( SIZE( c0, 2 ) ) )
prod = hpsi( cdesc%gzero, c0(:,:), c2 )
CALL mp_sum( prod, intra_image_comm )
IF( me_image == ib_owner( i ) ) THEN
ibl = ib_local( i )
lambda( ibl, : ) = prod( : )
END IF
DEALLOCATE( prod )
RETURN
END SUBROUTINE update_clambda
!=----------------------------------------------------------------------------=!
END MODULE wave_constrains
!=----------------------------------------------------------------------------=!