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

- clean-ups, comments and merging 

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1829 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
cavazzon 2005-04-22 15:23:19 +00:00
parent e1412a35e0
commit 193516156e
44 changed files with 1077 additions and 1627 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CPV/Makefile
View File

@ -9,7 +9,6 @@ bessel.o \
berryion.o \
bforceion.o \
blacs.o \
bmesh.o \
brillouin.o \
cg.o \
cg_sub.o \

14
CPV/bmesh.f90
View File

@ -27,14 +27,12 @@
!=----------------------------------------------------------------------------=!
SUBROUTINE bands_procgrid_init( )
!
! This subroutine initialize the processor grid, among which the electronic
! bands are distributed. Note that bands are not distributed all across the
! execution, but only in some subroutines, and only if it is convenient to
! do so.
!
!
! This subroutine initialize the processor grid, among which the electronic
! bands are distributed. Note that bands are not distributed all across the
! execution, but only in some subroutines, and only if it is convenient to
! do so.
!
USE mp_global, ONLY: mpime, nproc, group
USE processors_grid_module, ONLY: grid_init

4
CPV/chargedensity.f90
View File

@ -306,7 +306,7 @@
!=----------------------------------------------------------------------=!
! BEGIN manual
SUBROUTINE rhoofr (gv, kp, c0, cdesc, fi, rhoe, desc, box)
SUBROUTINE rhoofr (kp, c0, cdesc, fi, rhoe, desc, box)
! this routine computes:
! rhoe = normalized electron density in real space
@ -336,7 +336,6 @@
USE mp_global, ONLY: mpime
USE turbo, ONLY: tturbo, nturbo, turbo_states
USE cell_module, ONLY: boxdimensions
USE cp_types, ONLY: recvecs
USE wave_types, ONLY: wave_descriptor
USE charge_types, ONLY: charge_descriptor
USE io_global, ONLY: stdout
@ -347,7 +346,6 @@
! ... declare subroutine arguments
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
COMPLEX(dbl) :: c0(:,:,:,:)
TYPE (boxdimensions), INTENT(IN) :: box

41
CPV/chargemix.f90
View File

@ -27,7 +27,7 @@
! SUBROUTINE allocate_charge_mix(ng_l)
! SUBROUTINE deallocate_charge_mix
! SUBROUTINE charge_mix_print_info(unit)
! SUBROUTINE newrho(rhoe,nfi,tcel,drho,gv)
! SUBROUTINE newrho(rhoe,nfi,tcel,drho)
! SUBROUTINE invgen(aa,dimaa)
! ----------------------------------------------
! END manual
@ -141,7 +141,7 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE newrho(rhoe, drho, gv, nfi)
SUBROUTINE newrho(rhoe, drho, nfi)
! (describe briefly what this routine does...)
! ----------------------------------------------
@ -150,7 +150,8 @@
USE fft, ONLY: pfwfft, pinvfft
USE ions_base, ONLY: nsp
USE cell_base, ONLY: tpiba2
USE cp_types, ONLY: recvecs
USE reciprocal_vectors, ONLY: gstart, gzero, g
USE gvecp, ONLY: ngm
USE wave_base, ONLY: scalw
USE mp_global, ONLY: group
USE io_global, ONLY: stdout
@ -159,28 +160,24 @@
IMPLICIT NONE
! ... declare subroutine arguments
TYPE (recvecs), INTENT(IN) :: gv
REAL(dbl), INTENT(INOUT) :: rhoe(:,:,:)
REAL(dbl), INTENT(OUT) :: drho
INTEGER, INTENT(IN) :: nfi
! ... declare other variables
COMPLEX(dbl) :: dr
COMPLEX(dbl) :: rhoout(gv%ng_l)
COMPLEX(dbl) :: rhoout(ngm)
REAL(dbl) :: g02, g12, ar, den, num, rsc
REAL(dbl) :: alpha(daamax)
REAL(dbl), ALLOCATABLE :: aa(:,:)
REAL(dbl), ALLOCATABLE :: rho_old(:,:,:)
INTEGER :: ns, sp, is, ism, i, ig, gstart
INTEGER :: ns, sp, is, ism, i, ig
LOGICAL, SAVE :: tfirst = .TRUE.
INTEGER, SAVE :: dimaa, dimaaold, nrho_t, ierr
! ... end of declarations
! ----------------------------------------------
gstart = 1
IF( gv%gzero ) gstart = 2
IF( nfi /= 0 .AND. tfirst ) THEN
CALL errore(' newrho ', ' not initialized ', nfi )
@ -188,19 +185,19 @@
ELSE IF( nfi == 0 )THEN
IF( tfirst ) THEN
CALL allocate_charge_mix( gv%ng_l )
CALL allocate_charge_mix( ngm )
END IF
! ... define array chmix = A * G^2 / (G^2 + G_0^2) and metric = (G^2 + G_1^2) / G^2
g02 = g0chmix2 / tpiba2
g12 = g1met2 / tpiba2
IF(gv%gzero) THEN
IF(gzero) THEN
chmix(1) = 0.0d0
metric(1) = 0.0d0
END IF
DO ig = gstart, gv%ng_l
chmix(ig) = achmix * gv%hg_l(ig) / (gv%hg_l(ig)+g02)
metric(ig) = (gv%hg_l(ig)+g12) / gv%hg_l(ig)
DO ig = gstart, ngm
chmix(ig) = achmix * g(ig) / (g(ig)+g02)
metric(ig) = (g(ig)+g12) / g(ig)
END DO
tfirst = .FALSE.
@ -235,18 +232,18 @@
WRITE( stdout, fmt='( 3X,"charge mixing of order 1")' )
DO ig = gstart, gv%ng_l
DO ig = gstart, ngm
dr = rhoout(ig) - rho(ig,1)
rr(ig,1) = dr
rhoout(ig) = rho(ig,1) + chmix(ig) * dr
rho(ig,is) = rhoout(ig)
END DO
IF( gv%gzero ) THEN
IF( gzero ) THEN
rhoout(1) = rho(1,1)
rr(1,1) = (0.d0,0.d0)
END IF
IF( daamax /= 1 )THEN
rsc = scalw(gv%gzero, rr(:,1), rr(:,1), metric)
rsc = scalw(gzero, rr(:,1), rr(:,1), metric)
aa_save(1, 1) = rsc
END IF
@ -261,10 +258,10 @@
WRITE( stdout, fmt='( 3X,"charge mixing of order ",I2)' ) dimaa
DO ig = gstart, gv%ng_l
DO ig = gstart, ngm
rr(ig,ism) = rhoout(ig) - rho(ig,ism)
END DO
IF(gv%gzero) THEN
IF(gzero) THEN
rr(1,ism) = (0.d0, 0.d0)
END IF
@ -277,7 +274,7 @@
! ... Compute new matrix A
DO i = 1, dimaa
rsc = scalw(gv%gzero,rr(:,i),rr(:,ism),metric)
rsc = scalw(gzero,rr(:,i),rr(:,ism),metric)
aa(i,ism)= rsc
aa(ism,i)= rsc
END DO
@ -295,14 +292,14 @@
DEALLOCATE( aa, STAT=ierr )
IF( ierr /= 0 ) CALL errore(' newrho ', ' deallocating aa ', ierr)
DO ig = gstart, gv%ng_l
DO ig = gstart, ngm
rhoout(ig) = (0.d0,0.d0)
DO i = 1, dimaa
rhoout(ig) = rhoout(ig) + alpha(i) * ( rho(ig,i) + chmix(ig) * rr(ig,i) )
END DO
rho(ig,is) = rhoout(ig)
END DO
IF(gv%gzero) THEN
IF(gzero) THEN
rhoout(1) = rho(1,1)
END IF

28
CPV/chi2.f90
View File

@ -34,21 +34,21 @@
IF( ALLOCATED( rhochi ) ) DEALLOCATE(rhochi)
end subroutine deallocate_chi2
SUBROUTINE PRINTCHI2(gv,box)
SUBROUTINE PRINTCHI2(box)
USE cell_base, ONLY: tpiba
use cp_types, ONLY: recvecs
USE cell_module, only: boxdimensions
use mp, ONLY: mp_sum
USE mp_global, ONLY: nproc, mpime, group
USE io_files, ONLY: chiunit, chifile
USE reciprocal_vectors, ONLY: gstart, gx
USE gvecp, ONLY: ngm
IMPLICIT NONE
!---------------------------------------------------ARGUMENT
type (recvecs), intent(in) :: gv
type (boxdimensions), intent(in) :: box
REAL(dbl) GX,GY,GZ
REAL(dbl) GXR,GYR,GZR
!---------------------------------------------------COMMON
@ -69,28 +69,28 @@
!=======================================================================
CHI = 0.0d0
DO IG=gv%gstart,gv%ng_l
DO IG = gstart, ngm
ctmp = CONJG(RHOCHI(IG))*VLOCAL(IG)
do i=1,3
GX = gv%gx_l(i,IG)*TPIBA
GXR = gx(i,IG)*TPIBA
do j=1,3
GY = gv%gx_l(j,IG)*TPIBA
GYR = gx(j,IG)*TPIBA
do k=1,3
GZ = gv%gx_l(k,IG)*TPIBA
CHI(i,j,k) = CHI(i,j,k) - GX*GY*GZ*ctmp
GZR = gx(k,IG)*TPIBA
CHI(i,j,k) = CHI(i,j,k) - GXR*GYR*GZR*ctmp
end do
end do
end do
END DO
DO IG=gv%gstart,gv%ng_l
DO IG = gstart, ngm
ctmp = RHOCHI(IG)*CONJG(VLOCAL(IG))
do i=1,3
GX = -gv%gx_l(i,IG)*TPIBA
GXR = -gx(i,IG)*TPIBA
do j=1,3
GY = -gv%gx_l(j,IG)*TPIBA
GYR = -gx(j,IG)*TPIBA
do k=1,3
GZ = -gv%gx_l(k,IG)*TPIBA
CHI(i,j,k) = CHI(i,j,k) - GX*GY*GZ*ctmp
GZR = -gx(k,IG)*TPIBA
CHI(i,j,k) = CHI(i,j,k) - GXR*GYR*GZR*ctmp
end do
end do
end do

43
CPV/cp_fpmd.f90
View File

@ -310,16 +310,21 @@ end module qrl_mod
end subroutine
!-----------------------------------------------------------------------
subroutine gcalb(b1b,b2b,b3b)
subroutine gcalb( alatb, b1b_ , b2b_ , b3b_ )
!-----------------------------------------------------------------------
!
use control_flags, only: iprint
use gvecb
!
implicit none
real(kind=8) b1b(3),b2b(3),b3b(3)
real(kind=8), intent(in) :: alatb, b1b_ (3), b2b_ (3), b3b_ (3)
real(kind=8) :: b1b(3), b2b(3), b3b(3)
!
integer i, i1,i2,i3,ig
b1b = b1b_ * alatb
b2b = b2b_ * alatb
b3b = b3b_ * alatb
!
! calculation of gxb(3,ngbx)
!
@ -334,7 +339,7 @@ end module qrl_mod
enddo
!
return
end
end subroutine
!-------------------------------------------------------------------------
@ -1010,7 +1015,7 @@ END SUBROUTINE
!-------------------------------------------------------------------------
subroutine gcal(b1,b2,b3,gmax)
subroutine gcal( alat, b1_ , b2_ , b3_ , gmax )
!-----------------------------------------------------------------------
! calculates the values of g-vectors to be assigned to the lattice
! points generated in subroutine ggen. these values are derived
@ -1024,18 +1029,23 @@ END SUBROUTINE
!
! the g's are in units of 2pi/a.
!
use constants, only: tpi
use control_flags, only: iprint
use reciprocal_vectors, only: g, gx, mill_l
use gvecp, only: ngm
use gvecw, only: ngw
use gvecw, only: ggp, agg => ecutz, sgg => ecsig, e0gg => ecfix
use cell_base, only: tpiba2
implicit none
!
real(kind=8) b1(3),b2(3),b3(3), gmax
real(kind=8) alat, b1_(3),b2_(3),b3_(3), gmax
real(kind=8), external :: erf
real(kind=8) b1(3),b2(3),b3(3), tpiba2
!
integer i1,i2,i3,ig
b1 = b1_ * alat
b2 = b2_ * alat
b3 = b3_ * alat
!
! calculation of gx(3,ng)
!
@ -1050,6 +1060,8 @@ END SUBROUTINE
g(ig)=gx(1,ig)**2 + gx(2,ig)**2 + gx(3,ig)**2
if(g(ig).gt.gmax) gmax=g(ig)
enddo
tpiba2 = ( tpi / alat ) ** 2
!
do ig=1,ngw
ggp(ig) = g(ig) + &
@ -1078,25 +1090,24 @@ END SUBROUTINE
INTEGER :: i
REAL(dbl) :: alatb, b1b(3),b2b(3),b3b(3)
alatb = alat/nr1*nr1b
tpibab = 2.d0*pi/alatb
alatb = alat / nr1*nr1b
tpibab = 2.d0*pi / alatb
do i=1,3
a1b(i)=a1(i)/nr1*nr1b
a2b(i)=a2(i)/nr2*nr2b
a3b(i)=a3(i)/nr3*nr3b
enddo
omegab=omega/nr1*nr1b/nr2*nr2b/nr3*nr3b
!
call recips(a1b,a2b,a3b,b1b,b2b,b3b)
b1b = b1b * alatb
b2b = b2b * alatb
b3b = b3b * alatb
call gcalb(b1b,b2b,b3b)
call recips( a1b, a2b, a3b, b1b, b2b, b3b )
!
call gcalb( alatb, b1b, b2b, b3b )
!
do i=1,3
ainvb(1,i)=b1b(i)/alatb
ainvb(2,i)=b2b(i)/alatb
ainvb(3,i)=b3b(i)/alatb
ainvb(1,i)=b1b(i)
ainvb(2,i)=b2b(i)
ainvb(3,i)=b3b(i)
end do
RETURN

12
CPV/cpr.f90
View File

@ -92,7 +92,7 @@
use gvecs, only: ngs
use gvecb, only: ngb
use gvecw, only: ngw
use reciprocal_vectors, only: gstart
use reciprocal_vectors, only: gstart, mill_l
use ions_base, only: na, nat, pmass, nax, nsp, rcmax
use ions_base, only: ind_srt, ions_cofmass, ions_kinene, ions_temp, ions_thermal_stress
use ions_base, only: ions_vrescal, fricp, greasp, iforce, ions_shiftvar
@ -137,6 +137,7 @@
electrons_nose_nrg, electrons_nose_shiftvar, electrons_nosevel, electrons_noseupd
use from_scratch_module, only: from_scratch
use from_restart_module, only: from_restart
use phase_factors_module, only: strucf
! wavefunctions
!
@ -189,7 +190,7 @@
real(kind=8) &
& tempp, fccc, savee, saveh, savep, &
& enthal, epot, epre, enow, econs, econt, &
& ettt, ccc, bigr, dt2bye, dt2hbe
& ettt, ccc, bigr, dt2bye
real(kind=8) ekinc0, ekinp, ekinpr, ekincm, ekinc
real(kind=8) temps(nsx)
real(kind=8) ekinh, temphc, temp1, temp2, randy
@ -252,13 +253,12 @@
call init_dimensions( )
dt2bye = dt2/emass
dt2hbe = dt2by2/emass
etot_out = 0.0d0
tfirst = .true.
tlast = .false.
nacc = 5
call init ( ibrav, ndr, nbeg, tfirst, tau0, taus )
call init_geometry ( )
WRITE( stdout,*) ' out from init'
@ -465,7 +465,7 @@
!
! strucf calculates the structure factor sfac
!
call strucf(ei1,ei2,ei3,sfac)
call strucf( sfac, ei1, ei2, ei3, mill_l, ngs )
if (thdyn) call formf(tfirst,eself)
if( tefield ) then
@ -683,7 +683,7 @@
if( thdyn ) then
hold = h
h = hnew
call newinit(ibrav)
call newinit( h )
call newnlinit
else
hold = h

146
CPV/cptypes.f90
View File

@ -40,11 +40,6 @@
! BEGIN manual
! TYPE DEFINITIONS
!! ... Structure Factor
TYPE structure_factor
COMPLEX(dbl), POINTER :: s(:) ! S(G) = sum_I exp(i G dot R_I)
END TYPE
!! ... pseudopotential
TYPE pseudo
!! ... local part
@ -73,53 +68,11 @@
END TYPE pseudo
!! ... arrays of reciprocal space vectors (G and k+G)
TYPE recvecs
LOGICAL gzero
LOGICAL hspace
INTEGER gstart ! 2 if hg_l(1) = 0, 1 otherwise
INTEGER ng_l ! local number of G vectors
INTEGER ng_g ! global number of G vectors
INTEGER ngw_l
INTEGER ngw_g
!! ... quantities related to G vectors
REAL(dbl), POINTER :: hg_l(:) ! length squared of G vectors
REAL(dbl), POINTER :: gx_l(:,:) ! components of G vectors
! first index: x,y,z
! second index: G vector
INTEGER, POINTER :: mill(:,:) ! indices for FFT
! first index: x,y,z
! second index: G vector
INTEGER, POINTER :: ig(:) ! Global indices of the G vectors
! index: G vector
REAL(dbl), POINTER :: khgcutz_l(:,:) ! smooth cutoff factor index: G vector
! first index: G vector
! second index: k point
!! ... quantities related to k+G vectors
REAL(dbl), POINTER :: kg_mask_l(:,:) ! cutoff mask
REAL(dbl), POINTER :: khg_l(:,:) ! length squared of k+G
REAL(dbl), POINTER :: kgx_l(:,:,:) ! components of k+G
! first index: G vector
! second index: x,y,z
! third index: k point
REAL(dbl) :: bi1(3), bi2(3), bi3(3) ! Initial reciprocal lattice base (used to count the g's)
REAL(dbl) :: b1(3), b2(3), b3(3) ! Actual reciprocal lattice base (used to determine g2 )
END TYPE recvecs
! ----------------------------------------------
! END manual
PUBLIC :: pseudo, recvecs, pseudo_ncpp
PUBLIC :: pseudo, pseudo_ncpp
PUBLIC :: allocate_pseudo, deallocate_pseudo
PUBLIC :: allocate_recvecs, deallocate_rvecs
! end of module-scope declarations
! ----------------------------------------------
@ -129,103 +82,6 @@
! subroutines
SUBROUTINE allocate_recvecs(gv,ngl,ngg,ngwl,ngwg, tk,nk)
INTEGER, INTENT(IN) :: ngl,ngg,ngwl,ngwg,nk
LOGICAL, INTENT(IN) :: tk
TYPE (recvecs), INTENT(OUT) :: gv
INTEGER :: ierr
gv%gstart = -100000
gv%gzero = .TRUE.
gv%ng_l = ngl
gv%ng_g = ngg
gv%ngw_l = ngwl
gv%ngw_g = ngwg
!ALLOCATE(gv%hg_l(ngl), STAT=ierr)
!IF( ierr /= 0 ) CALL errore(' allocate_recvecs ', ' allocating %hg_l ', ierr )
NULLIFY( gv%hg_l )
!ALLOCATE(gv%gx_l(ngl,3), STAT=ierr)
!IF( ierr /= 0 ) CALL errore(' allocate_recvecs ', ' allocating %gx_l ', ierr )
NULLIFY( gv%gx_l )
IF (tk) THEN
gv%hspace = .FALSE.
ELSE
gv%hspace = .TRUE.
END IF
IF (tk) THEN
ALLOCATE(gv%kg_mask_l(ngwl,nk), STAT=ierr)
ELSE
ALLOCATE(gv%kg_mask_l(1,1), STAT=ierr)
END IF
IF( ierr /= 0 ) CALL errore(' allocate_recvecs ', ' allocating %kg_mask_l ', ierr )
ALLOCATE(gv%kgx_l(3,ngwl,nk), STAT=ierr)
IF( ierr /= 0 ) CALL errore(' allocate_recvecs ', ' allocating %kgx_l ', ierr )
ALLOCATE(gv%khg_l(ngwl,nk), STAT=ierr)
IF( ierr /= 0 ) CALL errore(' allocate_recvecs ', ' allocating %khg_l ', ierr )
ALLOCATE(gv%khgcutz_l(ngwl,nk), STAT=ierr)
IF( ierr /= 0 ) CALL errore(' allocate_recvecs ', ' allocating %khgcutz_l ', ierr )
NULLIFY( gv%mill )
NULLIFY( gv%ig )
!ALLOCATE(gv%mill(3,ngl), STAT=ierr)
!IF( ierr /= 0 ) CALL errore(' allocate_recvecs ', ' allocating %mill ', ierr )
!ALLOCATE(gv%ig(ngl), STAT=ierr)
!IF( ierr /= 0 ) CALL errore(' allocate_recvecs ', ' allocating %ig ', ierr )
RETURN
END SUBROUTINE allocate_recvecs
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE deallocate_rvecs(gv)
TYPE (recvecs), INTENT(INOUT) :: gv
INTEGER :: ierr
gv%gstart = 0
gv%ng_l = 0
gv%ng_g = 0
gv%ngw_l = 0
gv%ngw_g = 0
IF(ASSOCIATED(gv%hg_l)) THEN
NULLIFY( gv%hg_l )
END IF
IF(ASSOCIATED(gv%gx_l)) THEN
NULLIFY( gv%gx_l )
END IF
IF(ASSOCIATED(gv%kg_mask_l)) THEN
DEALLOCATE(gv%kg_mask_l, STAT=ierr)
IF( ierr /= 0 ) CALL errore(' deallocate_rvecs ', ' deallocating %kg_mask_l ', ierr )
END IF
IF(ASSOCIATED(gv%kgx_l)) THEN
DEALLOCATE(gv%kgx_l, STAT=ierr)
IF( ierr /= 0 ) CALL errore(' deallocate_rvecs ', ' deallocating %kgx_l ', ierr )
END IF
IF(ASSOCIATED(gv%khg_l)) THEN
DEALLOCATE(gv%khg_l, STAT=ierr)
IF( ierr /= 0 ) CALL errore(' deallocate_rvecs ', ' deallocating %khg_l ', ierr )
END IF
IF(ASSOCIATED(gv%khgcutz_l)) THEN
DEALLOCATE(gv%khgcutz_l, STAT=ierr)
IF( ierr /= 0 ) CALL errore(' deallocate_rvecs ', ' deallocating %khgcutz_l ', ierr )
END IF
IF(ASSOCIATED(gv%mill)) THEN
NULLIFY( gv%mill )
END IF
IF(ASSOCIATED(gv%ig)) THEN
NULLIFY( gv%ig )
END IF
RETURN
END SUBROUTINE deallocate_rvecs
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE allocate_pseudo(ps,nsp,ng,ngw,nk,lnlx,ngh,tcc)

52
CPV/diis.f90
View File

@ -27,10 +27,10 @@
! SUBROUTINE diis_print_info(unit)
! SUBROUTINE converg(c,gemax,cnorm,tprint)
! SUBROUTINE converg_kp(c,weight,gemax,cnorm,tprint)
! SUBROUTINE hesele(svar2,gv,vpp,wsg,wnl,eigr,vps,ik)
! SUBROUTINE simupd(gemax,doions,gv,c0,cgrad,svar0,svarm,svar2,etot,f, &
! SUBROUTINE hesele(svar2,vpp,wsg,wnl,eigr,vps,ik)
! SUBROUTINE simupd(gemax,doions,c0,cgrad,svar0,svarm,svar2,etot,f, &
! eigr,vps,wsg,wnl,treset,istate,cnorm,eold,ndiis,nowv)
! SUBROUTINE simupd_kp(gemax,doions,gv,c0,cgrad,svar0,svarm,svar2, &
! SUBROUTINE simupd_kp(gemax,doions,c0,cgrad,svar0,svarm,svar2, &
! etot,f,eigr,vps,wsg,wnl, treset,istate,cnorm, &
! eold,ndiis,nowv)
! SUBROUTINE updis(ndiis,nowv,nsize,max_diis,iact)
@ -306,7 +306,7 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE hesele(svar2, gv, vpp, wsg, wnl, eigr, sfac, vps, ik)
SUBROUTINE hesele(svar2, vpp, wsg, wnl, eigr, sfac, vps, ik)
! (describe briefly what this routine does...)
! ----------------------------------------------
@ -315,16 +315,16 @@
USE cell_base, ONLY: tpiba2
USE gvecw, ONLY: tecfix
USE pseudopotential, ONLY: tl, l2ind, nspnl, lm1x
USE cp_types, ONLY: recvecs
USE ions_base, ONLY: nsp, na
USE mp_global, ONLY: group
USE mp, ONLY: mp_sum, mp_max
USE reciprocal_vectors, ONLY: gstart, gzero
USE reciprocal_space_mesh, ONLY: gkmask_l, gkcutz_l, gk_l
IMPLICIT NONE
! ... declare subroutine arguments
REAL(dbl) :: svar2
TYPE (recvecs) :: gv
REAL(dbl) :: wnl(:,:,:,:), wsg(:,:), vps(:,:)
REAL(dbl) :: vpp(:)
COMPLEX(dbl) :: sfac(:,:)
@ -350,9 +350,9 @@
! ... local potential
vp = 0.0d0
IF( .NOT. PRESENT(ik) ) THEN
IF(gv%gzero) THEN
IF(gzero) THEN
DO i = 1, nsp
vp = vp + REAL( sfac(i,1) ) * vps(1,i)
vp = vp + REAL( sfac(1,i) ) * vps(1,i)
END DO
END IF
CALL mp_sum(vp, group)
@ -382,26 +382,24 @@
IF( PRESENT(ik) ) ikk = ik
IF(tecfix) THEN
vpp(:) = vpp(:) + ftpi * gv%khgcutz_l( 1:SIZE(vpp), ikk )
vpp(:) = vpp(:) + ftpi * gkcutz_l( 1:SIZE(vpp), ikk )
ELSE
vpp(:) = vpp(:) + ftpi * gv%khg_l( 1:SIZE(vpp), ikk )
vpp(:) = vpp(:) + ftpi * gk_l( 1:SIZE(vpp), ikk )
END IF
WHERE ( ABS( vpp ) .LT. hthrs_diis ) vpp = hthrs_diis
vpp = 1.0d0 / vpp
IF ( PRESENT(ik) ) vpp(:) = svar2 * gv%kg_mask_l(:,ik)
IF ( PRESENT(ik) ) vpp(:) = svar2 * gkmask_l(:,ik)
RETURN
END SUBROUTINE hesele
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE fermi_diis( ent, gv, kp, occ, nb, nel, eig, wke, efermi, sume, temp)
SUBROUTINE fermi_diis( ent, kp, occ, nb, nel, eig, wke, efermi, sume, temp)
USE brillouin, ONLY: kpoints
USE electrons_module, ONLY: fermi_energy
USE cp_types, ONLY: recvecs
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
REAL(dbl) :: occ(:,:,:)
REAL(dbl) :: wke(:,:,:)
@ -424,7 +422,7 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE simupd(gemax,doions,gv,c0,cgrad,cdesc,svar0,svarm,svar2,etot,f, &
SUBROUTINE simupd(gemax,doions,c0,cgrad,cdesc,svar0,svarm,svar2,etot,f, &
eigr,sfac,vps,wsg,wnl,treset,istate,cnorm,eold,ndiis,nowv)
! this routine performs a DIIS step
@ -432,11 +430,12 @@
! ----------------------------------------------
! ... declare modules
USE cp_types, ONLY: recvecs
USE wave_types, ONLY: wave_descriptor
USE mp_global, ONLY: group
USE io_global, ONLY: stdout
USE mp, ONLY: mp_sum, mp_max
USE gvecw, ONLY: ngw
USE reciprocal_vectors, ONLY: gstart, gzero
IMPLICIT NONE
@ -449,7 +448,6 @@
! hthrs_diis (REAL(dbl) ) minimum value for a hessel matrix elm.
! etot (REAL(dbl) ) Kohn-Sham energy
TYPE (recvecs) :: gv
COMPLEX(dbl), INTENT(INOUT) :: c0(:,:,:), cgrad(:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
COMPLEX(dbl) :: sfac(:,:)
@ -462,11 +460,11 @@
INTEGER :: ndiis, nowv
! ... declare other variables
COMPLEX(dbl), ALLOCATABLE :: cm(:,:) ! cm(gv%ngw_l,c0(1)%nb_l)
COMPLEX(dbl), ALLOCATABLE :: cm(:,:) ! cm(ngw,c0(1)%nb_l)
REAL(dbl), ALLOCATABLE :: bc(:,:) ! bc(max_diis+1,max_diis+1)
REAL(dbl), ALLOCATABLE :: vc(:) ! vc(max_diis+1)
REAL(dbl), ALLOCATABLE :: fm1(:)
REAL(dbl), ALLOCATABLE :: vpp(:) ! vpp(gv%ngw_l)
REAL(dbl), ALLOCATABLE :: vpp(:) ! vpp(ngw)
REAL(dbl) cnorm
REAL(dbl) var2,rc0rc0,ff,vvpp,fff,rri
REAL(dbl) rrj,rii,rij,r1,r2
@ -553,13 +551,13 @@
CALL update_diis_buffers( c0, cgrad, cdesc, nowv)
! ... calculate the new Hessian
CALL hesele(svar2,gv,vpp,wsg,wnl,eigr,sfac,vps)
CALL hesele(svar2,vpp,wsg,wnl,eigr,sfac,vps)
! ... set up DIIS matrix
ff=1.0d0
bc(1:nsize,1:nsize) = 0.0d0
DO l=gv%gstart,gv%ngw_l
DO l=gstart,ngw
vvpp=vpp(l)*vpp(l)
DO k=1, cdesc%nbl( 1 )
IF(oqnr_diis ) ff = fm1(k)
@ -576,7 +574,7 @@
END DO
END DO
IF(gv%gzero) THEN
IF(gzero) THEN
DO k = 1, cdesc%nbl( 1 )
IF(oqnr_diis ) ff = fm1(k)
fff=ff*ff*vpp(1)*vpp(1)
@ -641,7 +639,7 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE simupd_kp(gemax,doions,gv,kp,c0,cgrad,cdesc,svar0,svarm,svar2, &
SUBROUTINE simupd_kp(gemax,doions,kp,c0,cgrad,cdesc,svar0,svarm,svar2, &
etot,occ,eigr,sfac,vps,wsg,wnl, treset,istate,cnorm, &
eold,ndiis,nowv)
@ -650,7 +648,6 @@
! ----------------------------------------------
! ... declare modules
USE cp_types, ONLY: recvecs
USE wave_types, ONLY: wave_descriptor
USE brillouin, ONLY: kpoints
USE mp_global, ONLY: group
@ -668,7 +665,6 @@
! hthrs_diis (REAL(dbl) ) minimum value for a Hessian matrix element
! etot (REAL(dbl) ) Kohn-Sham energy
TYPE (recvecs) gv
TYPE (kpoints) kp
COMPLEX(dbl), INTENT(INOUT) :: c0(:,:,:), cgrad(:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
@ -775,13 +771,13 @@
END WHERE
! ... calculate the new Hessian
CALL hesele(svar2,gv,vpp,wsg,wnl,eigr,sfac,vps,ik)
CALL hesele(svar2,vpp,wsg,wnl,eigr,sfac,vps,ik)
! ... set up DIIS matrix
ff=1.0d0
bc(1:nsize,1:nsize) = CMPLX(0.0d0,0.0d0)
DO l=1,gv%ngw_l
DO l=1, cdesc%ngwl
vvpp = vpp(l)*vpp(l)
DO k=1, cdesc%nbl( 1 )
IF(oqnr_diis) ff = fm1(k)
@ -820,7 +816,7 @@
c0(:,:,ik) = CMPLX(0.0d0,0.0d0)
DO i=1,nsize-1
DO j=1,cdesc%nbl( 1 )
DO k=1,gv%ngw_l
DO k=1,cdesc%ngwl
cm(k,j) = cm(k,j) + CONJG(vc(i)) * grade(k,j,ik,i)
c0(k,j,ik) = c0(k,j,ik) + CONJG(vc(i)) * parame(k,j,ik,i)
END DO

4
CPV/electrons.f90
View File

@ -148,7 +148,7 @@
IMPLICIT NONE
INTEGER :: i, n1, n2, n3, ierr
INTEGER :: i, ierr
IF( band_first ) THEN
CALL errore(' bmeshset ',' module not initialized ',0)
@ -162,8 +162,10 @@
IF( mpime < MOD( n_emp, nproc ) ) n_emp_l( i ) = n_emp_l( i ) + 1
END DO
IF( ALLOCATED( ib_owner ) ) DEALLOCATE( ib_owner )
ALLOCATE( ib_owner( 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)
IF( ierr/=0 ) CALL errore( ' bmeshset ',' allocating ib_local ', ierr)

95
CPV/emptystates.f90
View File

@ -75,7 +75,7 @@
!=----------------------------------------------------------------------------=!
!
LOGICAL FUNCTION readempty( c_emp, wempt, gv )
LOGICAL FUNCTION readempty( c_emp, wempt )
! ... This subroutine reads empty states from unit emptyunit
@ -85,13 +85,12 @@
USE mp, ONLY: mp_bcast
USE mp_wave, ONLY: splitwf
USE io_files, ONLY: scradir
USE cp_types, ONLY: recvecs
USE reciprocal_vectors, ONLY: ig_l2g
IMPLICIT none
COMPLEX(dbl), INTENT(INOUT) :: c_emp(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: wempt
TYPE (recvecs) :: gv
LOGICAL :: exst
INTEGER :: ierr, ig, i, ik, nl, ispin
@ -153,7 +152,7 @@
READ(emptyunit) ( ctmp(ig), ig = 1, MIN( SIZE(ctmp), ngw_rd ) )
END IF
IF( i <= ne(ispin) ) THEN
CALL splitwf(c_emp(:,i,ik,ispin), ctmp, gv%ngw_l, gv%ig, mpime, nproc, root)
CALL splitwf(c_emp(:,i,ik,ispin), ctmp, wempt%ngwl, ig_l2g, mpime, nproc, root)
END IF
END DO
END DO
@ -175,7 +174,7 @@
!=----------------------------------------------------------------------------=!
!
SUBROUTINE writeempty( c_emp, wempt, gv )
SUBROUTINE writeempty( c_emp, wempt )
! ... This subroutine writes empty states to unit emptyunit
@ -183,13 +182,12 @@
USE mp_global, ONLY: mpime, nproc, group, root
USE mp_wave, ONLY: mergewf
USE io_files, ONLY: scradir
USE cp_types, ONLY: recvecs
USE reciprocal_vectors, ONLY: ig_l2g
COMPLEX(dbl), INTENT(IN) :: c_emp(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: wempt
TYPE (recvecs) :: gv
INTEGER :: ig, i, ik, nl, ne(2), ngwm_g, nk, ispin, nspin
INTEGER :: ig, i, ik, nl, ne(2), ngwm_g, nk, ispin, nspin, ngw
LOGICAL :: exst
COMPLEX(dbl), ALLOCATABLE :: ctmp(:)
!
@ -201,6 +199,7 @@
nk = wempt%nkl
nspin = wempt%nspin
ngwm_g = wempt%ngwt
ngw = wempt%ngwl
ne = 0
ne( 1 : nspin ) = wempt%nbl( 1 : nspin )
@ -222,7 +221,7 @@
DO ik = 1, nk
DO i = 1, ne(ispin)
ctmp = 0.0d0
CALL MERGEWF( c_emp(:,i,ik,ispin), ctmp(:), gv%ngw_l, gv%ig, mpime, nproc, root)
CALL MERGEWF( c_emp(:,i,ik,ispin), ctmp(:), ngw, ig_l2g, mpime, nproc, root)
IF( mpime == 0 ) THEN
WRITE (emptyunit) ( ctmp(ig), ig=1, ngwm_g )
END IF
@ -361,19 +360,19 @@
!
!================================================================
!
SUBROUTINE randomizza(c_occ, wfill, ampre, c_emp, wempt, gv, kp)
SUBROUTINE randomizza(c_occ, wfill, ampre, c_emp, wempt, kp)
USE wave_types, ONLY: wave_descriptor
USE reciprocal_vectors, ONLY: ig_l2g, ngw, ngwt
USE reciprocal_vectors, ONLY: ig_l2g
USE mp_global, ONLY: mpime, nproc, root
USE mp_wave, ONLY: splitwf
USE brillouin, ONLY: kpoints
USE cp_types, ONLY: recvecs
USE control_flags, ONLY: force_pairing
USE reciprocal_space_mesh, ONLY: gkmask_l
USE wave_functions, ONLY: wave_rand_init
! ... Arguments
COMPLEX(dbl), INTENT(INOUT) :: c_occ(:,:,:,:), c_emp(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: wfill, wempt
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
REAL(dbl) :: ampre
@ -382,41 +381,33 @@
! ... Locals
INTEGER :: ig_local
INTEGER :: ib, ik, ispin, i, ig, ntest, j, ispin_wfc
INTEGER :: ngw, ngwt
INTEGER :: ib, ik, ispin, ispin_wfc
LOGICAL :: tortho = .FALSE.
REAL(dbl) :: rranf1, rranf2
COMPLEX(dbl), ALLOCATABLE :: pwt( : )
!
ntest = gv%ngw_g / 4
! ... Subroutine body
! ... initialize the wave functions in such a way that the values
! ... of the components are independent on the number of processors
ngwt = wfill%ngwt
ngw = wfill%ngwl
ALLOCATE( pwt( ngwt ) )
DO ispin = 1, wempt%nspin
ispin_wfc = ispin
IF( force_pairing ) ispin_wfc = 1
DO ik = 1, wempt%nkl
c_emp(:,:,ik,ispin) = 0.0d0
DO ib = 1, wempt%nbl( ispin )
pwt( 1 ) = 0.0d0
DO ig = 2, ntest
rranf1 = 0.5d0 - rranf()
rranf2 = rranf()
pwt( ig ) = ampre * CMPLX(rranf1, rranf2)
END DO
CALL splitwf ( c_emp( :, ib, ik, ispin ), pwt, ngw, ig_l2g, mpime, nproc, 0 )
END DO
CALL wave_rand_init( c_emp( :, :, ik, ispin ) )
IF ( .NOT. kp%gamma_only ) THEN
! .. . set to zero all elements outside the cutoff sphere
DO ib = 1, wempt%nbl( ispin )
c_emp(:,ib,ik,ispin) = c_emp(:,ib,ik,ispin) * gv%kg_mask_l(:,ik)
c_emp(:,ib,ik,ispin) = c_emp(:,ib,ik,ispin) * gkmask_l(:,ik)
END DO
END IF
IF ( gv%gzero ) THEN
IF ( wempt%gzero ) THEN
c_emp(1,:,ik,ispin) = (0.0d0, 0.0d0)
END IF
CALL gram_empty( ispin, tortho, c_occ(:,:,ik,ispin_wfc), wfill, c_emp(:,:,ik,ispin), wempt )
@ -432,7 +423,7 @@
!=======================================================================
!
SUBROUTINE EMPTY_SD( tortho, atoms, gv, c_occ, wfill, c_emp, wempt, kp, vpot, eigr, ps)
SUBROUTINE EMPTY_SD( tortho, atoms, c_occ, wfill, c_emp, wempt, kp, vpot, eigr, ps)
USE wave_types, ONLY: wave_descriptor
USE wave_functions, ONLY: cp_kinetic_energy, crot, dft_kinetic_energy, fixwave
@ -456,8 +447,9 @@
USE atoms_type_module, ONLY: atoms_type
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE cp_types, ONLY: recvecs, pseudo
USE cp_types, ONLY: pseudo
USE control_flags, ONLY: force_pairing
USE reciprocal_space_mesh, ONLY: gkmask_l, gkx_l, gk_l
IMPLICIT NONE
@ -466,7 +458,6 @@
COMPLEX(dbl), INTENT(INOUT) :: c_occ(:,:,:,:), c_emp(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: wfill, wempt
TYPE (atoms_type), INTENT(INOUT) :: atoms ! ions structure
TYPE (recvecs), INTENT(IN) :: gv
REAL (dbl), INTENT(IN) :: vpot(:,:,:,:)
TYPE (kpoints), INTENT(IN) :: kp
TYPE (pseudo), INTENT(IN) :: ps
@ -477,9 +468,9 @@
!
INTEGER :: i, k, j, iter, ik, nk
INTEGER :: ngw, ngw_g, nspin, ispin, ispin_wfc
INTEGER :: nspin, ispin, ispin_wfc
INTEGER :: n_occ( wfill%nspin )
INTEGER :: ig, iprinte, iks, nrl, jl
INTEGER :: ig, iprinte, iks, nrl, jl, ngw
REAL(dbl) :: dek, ekinc, ekinc_old
REAL(dbl) :: ampre
REAL(dbl), ALLOCATABLE :: dt2bye( : )
@ -487,7 +478,7 @@
TYPE (projector) :: fnle( SIZE(c_emp, 3), SIZE(c_emp, 4) )
COMPLEX(dbl), ALLOCATABLE :: eforce(:,:,:,:), cp_emp(:,:,:,:)
REAL(dbl), ALLOCATABLE :: fi(:,:,:)
LOGICAL :: gamma, gzero
LOGICAL :: gamma
LOGICAL, SAVE :: exst
!
! ... SUBROUTINE BODY
@ -496,9 +487,7 @@
nk = wfill%nkl
nspin = wfill%nspin
ngw = wfill%ngwl
ngw_g = wfill%ngwt
n_occ = wfill%nbt
gzero = wfill%gzero
gamma = wfill%gamma
ampre = 0.001d0
@ -513,10 +502,10 @@
IF( ionode ) WRITE( stdout,56)
exst = readempty( c_emp, wempt, gv )
exst = readempty( c_emp, wempt )
! .. WRITE( stdout, * )' DEBUG empty 1 ', exst
IF( .NOT. exst ) THEN
CALL randomizza(c_occ, wfill, ampre, c_emp, wempt, gv, kp)
CALL randomizza(c_occ, wfill, ampre, c_emp, wempt, kp)
END IF
dt2bye = delt * delt / pmss
@ -537,17 +526,17 @@
DO ik = 1, kp%nkpt
CALL nlsm1( ispin, ps%wnl(:,:,:,ik), atoms, eigr, c_emp(:,:,ik,ispin), wempt, &
gv%khg_l(:,ik), gv%kgx_l(:,:,ik), fnle(ik,ispin))
gk_l(:,ik), gkx_l(:,:,ik), fnle(ik,ispin))
CALL dforce_all( ispin, c_emp(:,:,:,ispin), wempt, fi(:,:,ispin), eforce(:,:,:,ispin), &
gv, vpot(:,:,:,ispin), fnle(:,ispin), eigr, ps, ik)
vpot(:,:,:,ispin), fnle(:,ispin), eigr, ps, ik)
! ... Steepest descent
DO i = 1, n_emp
cp_emp(:,i,ik,ispin) = c_emp(:,i,ik,ispin) + dt2bye(:) * eforce(:, i, ik, ispin)
END DO
CALL fixwave( ispin, cp_emp(:,:,ik,ispin), wempt, gv%kg_mask_l(:,ik) )
CALL fixwave( ispin, cp_emp(:,:,ik,ispin), wempt, gkmask_l(:,ik) )
IF (tortho) THEN
@ -563,7 +552,7 @@
END DO
ekinc = ekinc + cp_kinetic_energy( ispin, cp_emp(:,:,:,ispin), c_emp(:,:,:,ispin), wempt, &
kp, gv%kg_mask_l, pmss, delt)
kp, gkmask_l, pmss, delt)
END DO SPIN_LOOP
@ -589,9 +578,9 @@
END DO ITERATIONS
CALL empty_eigs(tortho, atoms, c_emp, wempt, fi, vpot, eforce, fnle, ps, eigr, gv, kp)
CALL empty_eigs(tortho, atoms, c_emp, wempt, fi, vpot, eforce, fnle, ps, eigr, kp)
CALL writeempty( c_emp, wempt, gv )
CALL writeempty( c_emp, wempt )
CALL deallocate_projector(fnle)
DEALLOCATE( eforce )
@ -614,7 +603,7 @@
!
!=----------------------------------------------------------------------------=!
SUBROUTINE empty_eigs(tortho, atoms, c_emp, wempt, fi, vpot, eforce, fnle, ps, eigr, gv, kp)
SUBROUTINE empty_eigs(tortho, atoms, c_emp, wempt, fi, vpot, eforce, fnle, ps, eigr, kp)
USE wave_types, ONLY: wave_descriptor
USE wave_functions, ONLY: crot
@ -631,7 +620,8 @@
USE mp, ONLY: mp_sum
USE mp_global, ONLY: mpime, nproc, group
USE atoms_type_module, ONLY: atoms_type
USE cp_types, ONLY: recvecs, pseudo
USE cp_types, ONLY: pseudo
USE reciprocal_space_mesh, ONLY: gkx_l, gk_l
IMPLICIT NONE
@ -640,7 +630,6 @@
COMPLEX(dbl), INTENT(inout) :: c_emp(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: wempt
TYPE(atoms_type), INTENT(INOUT) :: atoms ! ions structure
TYPE (recvecs), INTENT(in) :: gv
REAL (dbl), INTENT(in) :: vpot(:,:,:,:), fi(:,:,:)
COMPLEX (dbl) :: eforce(:,:,:,:)
TYPE (kpoints), INTENT(in) :: kp
@ -653,9 +642,9 @@
!
INTEGER kk, i, k, j, iopt, iter, nwh, ik, nk, ibl
INTEGER ngw, ngw_g, n_occ, nspin, ispin
INTEGER ngw, n_occ, nspin, ispin
INTEGER ig, iprinte, iks, nrl, jl
LOGICAL gamma, gzero
LOGICAL gamma
REAL(dbl), ALLOCATABLE :: gam(:,:)
COMPLEX(dbl), ALLOCATABLE :: cgam(:,:)
@ -666,8 +655,6 @@
nspin = wempt%nspin
nk = wempt%nkl
ngw = wempt%ngwl
ngw_g = wempt%ngwt
gzero = wempt%gzero
gamma = wempt%gamma
!
@ -681,11 +668,11 @@
DO ik = 1, nk
CALL nlsm1( ispin, ps%wnl(:,:,:,ik), atoms, eigr, c_emp(:,:,ik,ispin), wempt, &
gv%khg_l(:,ik), gv%kgx_l(:,:,ik), fnle(ik,ispin))
gk_l(:,ik), gkx_l(:,:,ik), fnle(ik,ispin))
END DO
! ... Calculate | dH / dpsi(j) >
CALL dforce_all( ispin, c_emp(:,:,:,ispin), wempt, fi(:,:,ispin), eforce(:,:,:,ispin), gv, &
CALL dforce_all( ispin, c_emp(:,:,:,ispin), wempt, fi(:,:,ispin), eforce(:,:,:,ispin), &
vpot(:,:,:,ispin), fnle(:,ispin), eigr, ps)
DO ik = 1, kp%nkpt

39
CPV/exch_corr.f90
View File

@ -27,13 +27,15 @@
CONTAINS
!=----------------------------------------------------------------------------=!
SUBROUTINE v2gc(v2xc, grho, rhoer, vpot, gv)
SUBROUTINE v2gc(v2xc, grho, rhoer, vpot)
USE fft
USE fft_base, ONLY: dfftp
USE cell_base, ONLY: tpiba
USE cp_types
USE mp_global
USE reciprocal_vectors, ONLY: gstart, gx
USE gvecp, ONLY: ngm
!
implicit none
!
@ -41,7 +43,6 @@
REAL(dbl), intent(in) :: v2xc(:,:,:,:,:)
REAL(dbl), intent(in) :: grho(:,:,:,:,:)
REAL(dbl), intent(in) :: rhoer(:,:,:,:)
type (recvecs), intent(in) :: GV
!
integer :: ig, ipol, nxl, nyl, nzl, i, j, k, is, js, nspin
integer :: ldx, ldy, ldz
@ -62,8 +63,8 @@
!fac = REAL(nspin)
fac = 1.0d0
ALLOCATE( vtemp( gv%ng_l ) )
ALLOCATE( vtemp_pol( gv%ng_l ) )
ALLOCATE( vtemp( ngm ) )
ALLOCATE( vtemp_pol( ngm ) )
DO js = 1, nspin
!
@ -84,8 +85,8 @@
!
CALL pfwfft( vtemp_pol, psi )
!
DO ig = gv%gstart, gv%ng_l
vtemp(ig) = vtemp(ig) + vtemp_pol(ig) * CMPLX( 0.d0, tpiba * gv%gx_l( ipol, ig ) )
DO ig = gstart, ngm
vtemp(ig) = vtemp(ig) + vtemp_pol(ig) * CMPLX( 0.d0, tpiba * gx( ipol, ig ) )
END DO
!
END DO
@ -163,19 +164,19 @@
!=----------------------------------------------------------------------------=!
SUBROUTINE stress_xc( dexc, strvxc, sfac, vxc, grho, v2xc, &
gagx_l, gv, tnlcc, rhocp, box)
gagx_l, tnlcc, rhocp, box)
use ions_base, only: nsp
USE cell_module, only: boxdimensions
USE cell_base, ONLY: tpiba
USE cp_types, ONLY: recvecs
USE funct, ONLY: igcx, igcc
use ions_base, only: nsp
USE cell_module, only: boxdimensions
USE cell_base, ONLY: tpiba
USE funct, ONLY: igcx, igcc
USE reciprocal_vectors, ONLY: gstart, g
USE gvecp, ONLY: ngm
IMPLICIT NONE
! -- ARGUMENT
type (recvecs), intent(in) :: gv
type (boxdimensions), intent(in) :: box
LOGICAL :: tnlcc(:)
COMPLEX(dbl) :: vxc(:,:)
@ -206,18 +207,18 @@
IF ( ANY( tnlcc ) ) THEN
DO ig = gv%gstart, gv%ng_l
DO ig = gstart, ngm
tex1 = (0.0_dbl , 0.0_dbl)
DO is=1,nsp
IF ( tnlcc(is) ) THEN
tex1 = tex1 + sfac( is, ig ) * CMPLX(rhocp(ig,is))
tex1 = tex1 + sfac( ig, is ) * CMPLX(rhocp(ig,is))
END IF
END DO
tex2 = 0.0_dbl
DO ispin = 1, nspin
tex2 = tex2 + CONJG( vxc(ig, ispin) )
END DO
tex3 = REAL(tex1 * tex2) / SQRT(gv%hg_l(ig)) / tpiba
tex3 = REAL(tex1 * tex2) / SQRT( g( ig ) ) / tpiba
dexc = dexc + tex3 * gagx_l(:,ig)
END DO
dexc = dexc * 2.0_dbl * omega
@ -241,10 +242,9 @@
!=----------------------------------------------------------------------------=!
SUBROUTINE exch_corr_energy(rhoetr, rhoetg, grho, vpot, sxc, vxc, v2xc, gv)
SUBROUTINE exch_corr_energy(rhoetr, rhoetg, grho, vpot, sxc, vxc, v2xc)
USE kinds, ONLY: dbl
USE cp_types, ONLY: recvecs
USE grid_dimensions, ONLY: nr1l, nr2l, nr3l
USE funct, ONLY: igcx, igcc
@ -255,7 +255,6 @@
REAL(dbl) :: sxc ! E_xc energy
REAL(dbl) :: vxc ! SUM ( v(r) * rho(r) )
REAL (dbl) :: v2xc(:,:,:,:,:)
TYPE (recvecs), INTENT(IN) :: gv
REAL(dbl) :: ddot
INTEGER :: nspin, nnr, ispin, j, k
@ -271,7 +270,7 @@
!
IF( ( igcx > 0 ) .OR. ( igcc > 0 ) ) THEN
! ... vpot additional term for gradient correction
CALL v2gc( v2xc, grho, rhoetr, vpot, gv )
CALL v2gc( v2xc, grho, rhoetr, vpot )
END If
!

10
CPV/fft.f90
View File

@ -123,7 +123,6 @@
WRITE( stdout, fmt ="(/,3X,'PC3FFT TIMINGS')" )
WRITE( stdout,910)
WRITE( stdout,999) ( ( fft_timing(i,j), i = 1, 4), j = 1, 2 )
DEALLOCATE( fft_timing )
910 FORMAT(' FFTXP FFTYP FFTZP TRASP FFTXW FFTYW FFTZW TRASW')
999 FORMAT(8(F9.3))
@ -256,15 +255,6 @@
CALL psi2c( zwrk(1), SIZE( zwrk ), c(1), c(1), ng, 10 )
!DO ig=1,ng
! is = ( np( ig ) - 1 ) / dfftp%nr3x + 1
! k = MOD( ( np( ig ) - 1 ), dfftp%nr3x ) + 1
! WRITE( 201, fmt="( 4I6,4D22.14 )" ) ig, is, k, np( ig ), &
! C( ig ), zwrk( np( ig ) )
! ! C( ig ), zwrk( k + (is-1) * dfftp%nr3x )
!END DO
!CLOSE( 201 )
DEALLOCATE( zwrk )

5
CPV/fftdrv.f90
View File

@ -96,11 +96,6 @@
ldz = dfft%nr3x
nz_l = dfft%npp( mpime + 1 )
IF( .NOT. ALLOCATED( fft_timing ) ) THEN
ALLOCATE( fft_timing( 4, 2 ) )
fft_timing = 0.0d0
END IF
IF( FFT_MODE == FFT_MODE_POTE ) THEN
ns_l = dfft%nsp( mpime + 1 )
ELSE

86
CPV/forces.f90
View File

@ -13,19 +13,19 @@
! SISSA, Trieste, Italy - 1997-99
! Last modified: Wed Oct 13 15:28:58 MDT; 1999
! ----------------------------------------------
! SUBROUTINE dforce_p(ik,ib,c,df,da,f,gv,v,fnl,eigr,ps)
! SUBROUTINE dforce_p(ik,ib,c,df,da,f,v,fnl,eigr,ps)
! SUBROUTINE dforce1(co,ce,dco,dce,fio,fie,hg,v,psi_stored)
! SUBROUTINE dforce2(fi,fip1,df,da,gv,fnl,eigr,wsg,wnl)
! SUBROUTINE dforce2(fi,fip1,df,da,fnl,eigr,wsg,wnl)
!
! SUBROUTINE dforce_d(ik,ib,c,df,f,hg,v,fnl,eigr,ps)
! SUBROUTINE dforce1_d(co,dco,fi,gv,v)
! SUBROUTINE dforce2_d(i,fi,df,gv,fnl,eigr,wsg,wnl)
! SUBROUTINE dforce1_d(co,dco,fi,v)
! SUBROUTINE dforce2_d(i,fi,df,fnl,eigr,wsg,wnl)
!
! SUBROUTINE dforce1_kp(ib,ik,c0,df,fi,gv,v)
! SUBROUTINE dforce1_kp(ib,ik,c0,df,fi,gv,v)
! SUBROUTINE dforce2_kp(ib,ik,fi,df,gv,fnlk,eigr,wsg,wnl)
! SUBROUTINE dforce1_kp(ib,ik,c0,df,fi,v)
! SUBROUTINE dforce1_kp(ib,ik,c0,df,fi,v)
! SUBROUTINE dforce2_kp(ib,ik,fi,df,fnlk,eigr,wsg,wnl)
!
! SUBROUTINE dforce_all(c0,cgrad,gv,vpot,fnl,eigr,ps)
! SUBROUTINE dforce_all(c0,cgrad,vpot,fnl,eigr,ps)
!
! INTERFACE dforce
! ----------------------------------------------
@ -195,16 +195,16 @@
!=----------------------------------------------------------------------------=!
subroutine dforce_p( ik, ib, c, cdesc, f, df, da, gv, v, fnl, eigr, ps )
subroutine dforce_p( ik, ib, c, cdesc, f, df, da, v, fnl, eigr, ps )
USE wave_types, ONLY: wave_descriptor
USE turbo, ONLY: tturbo, nturbo, turbo_states
USE gvecw, ONLY: tecfix
USE reciprocal_space_mesh, ONLY: gkcutz_l, gkx_l, gk_l
implicit none
integer ik,ib
COMPLEX(dbl), INTENT(IN) :: c(:,:,:)
type (wave_descriptor), INTENT(IN) :: cdesc
COMPLEX(dbl) :: df(:), da(:)
type (recvecs) :: gv
REAL (dbl) :: v(:,:,:), f(:,:)
REAL (dbl) :: fnl(:,:,:)
type (pseudo) :: ps
@ -213,9 +213,9 @@
INTEGER :: istate
istate = (ib+1)/2
IF(tecfix) THEN
hg => gv%khgcutz_l(:,ik)
hg => gkcutz_l(:,ik)
ELSE
hg => gv%khg_l(:,ik)
hg => gk_l(:,ik)
END IF
IF( tturbo .AND. ( istate <= nturbo ) ) THEN
CALL dforce1(c(:,ib,ik), c(:,ib+1,ik), df, da, &
@ -225,30 +225,30 @@
f(ib,ik), f(ib+1,ik), hg, v)
END IF
CALL dforce2(f(ib,ik), f(ib+1,ik), df, da, fnl(:,:,ib), &
fnl(:,:,ib+1), gv%khg_l(:,ik), gv%kgx_l(:,:,ik), eigr, &
fnl(:,:,ib+1), gk_l(:,ik), gkx_l(:,:,ik), eigr, &
ps%wsg, ps%wnl(:,:,:,ik))
return
end subroutine
subroutine dforce_p_s(ib, c, cdesc, f, df, da, gv, v, fnl, eigr, wsg, wnl)
subroutine dforce_p_s(ib, c, cdesc, f, df, da, v, fnl, eigr, wsg, wnl)
USE wave_types, ONLY: wave_descriptor
USE turbo, ONLY: tturbo, nturbo, turbo_states
USE gvecw, ONLY: tecfix
USE reciprocal_space_mesh, ONLY: gkcutz_l, gkx_l, gk_l
IMPLICIT NONE
INTEGER, INTENT(IN) :: ib
COMPLEX(dbl), INTENT(IN) :: c(:,:)
type (wave_descriptor), INTENT(IN) :: cdesc
COMPLEX(dbl), INTENT(OUT) :: df(:), da(:)
TYPE (recvecs), INTENT(IN) :: gv
REAL (dbl), INTENT(IN) :: v(:,:,:), fnl(:,:,:), wnl(:,:,:), wsg(:,:), f(:)
COMPLEX(dbl), INTENT(IN) :: eigr(:,:)
REAL(dbl), POINTER :: hg(:)
INTEGER :: istate
istate = (ib+1)/2
IF(tecfix) THEN
hg => gv%khgcutz_l(:,1)
hg => gkcutz_l(:,1)
ELSE
hg => gv%khg_l(:,1)
hg => gk_l(:,1)
END IF
IF(tturbo.AND.(istate.LE.nturbo)) THEN
CALL dforce1(c(:,ib), c(:,ib+1), df, da, &
@ -257,7 +257,7 @@
CALL dforce1(c(:,ib), c(:,ib+1), df, da, f(ib), f(ib+1), hg, v)
END IF
CALL dforce2(f(ib), f(ib+1), df, da, fnl(:,:,ib), &
fnl(:,:,ib+1), gv%khg_l(:,1), gv%kgx_l(:,:,1), eigr, wsg, wnl)
fnl(:,:,ib+1), gk_l(:,1), gkx_l(:,:,1), eigr, wsg, wnl)
return
end subroutine
@ -337,7 +337,6 @@
IMPLICIT NONE
! ... declare subroutine arguments
TYPE (recvecs) gv
REAL(dbl), INTENT(IN) :: wnl(:,:,:)
COMPLEX(dbl), INTENT(IN) :: eigr(:,:)
REAL(dbl), INTENT(IN) :: fi, fnl(:,:),wsg(:,:)
@ -391,52 +390,52 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE dforce_d(ik, ib, c, cdesc, f, df, gv, v, fnl, eigr, ps)
SUBROUTINE dforce_d(ik, ib, c, cdesc, f, df, v, fnl, eigr, ps)
USE wave_types, ONLY: wave_descriptor
USE gvecw, ONLY: tecfix
USE reciprocal_space_mesh, ONLY: gkcutz_l, gkx_l, gk_l
IMPLICIT NONE
integer ik, ib
COMPLEX(dbl), INTENT(IN) :: c(:,:,:)
type (wave_descriptor), INTENT(IN) :: cdesc
COMPLEX(dbl) :: df(:)
type (recvecs) :: gv
REAL (dbl) :: v(:,:,:), f(:,:)
REAL (dbl) :: fnl(:,:,:)
type (pseudo) :: ps
COMPLEX(dbl) :: eigr (:,:)
REAL(dbl), POINTER :: hg(:)
IF(tecfix) THEN
hg => gv%khgcutz_l(:,ik)
hg => gkcutz_l(:,ik)
ELSE
hg => gv%khg_l(:,ik)
hg => gk_l(:,ik)
END IF
CALL dforce1_d(c(:,ib,ik), df, f(ib,ik), hg, v)
CALL dforce2_d(f(ib,ik), df, fnl(:,:,ib), gv%khg_l(:,ik), &
gv%kgx_l(:,:,ik), eigr, ps%wsg, ps%wnl(:,:,:,ik))
CALL dforce2_d(f(ib,ik), df, fnl(:,:,ib), gk_l(:,ik), &
gkx_l(:,:,ik), eigr, ps%wsg, ps%wnl(:,:,:,ik))
return
end subroutine
subroutine dforce_d_s(ib, c, cdesc, f, df, gv, v, fnl, eigr, wsg, wnl)
subroutine dforce_d_s(ib, c, cdesc, f, df, v, fnl, eigr, wsg, wnl)
USE wave_types, ONLY: wave_descriptor
USE gvecw, ONLY: tecfix
USE turbo, ONLY: tturbo, nturbo, turbo_states
USE reciprocal_space_mesh, ONLY: gkcutz_l, gkx_l, gk_l
IMPLICIT NONE
INTEGER, INTENT(IN) :: ib
COMPLEX(dbl), INTENT(IN) :: c(:,:)
type (wave_descriptor), INTENT(IN) :: cdesc
COMPLEX(dbl), INTENT(OUT) :: df(:)
TYPE (recvecs), INTENT(IN) :: gv
REAL (dbl), INTENT(IN) :: v(:,:,:), fnl(:,:,:), wnl(:,:,:), wsg(:,:), f(:)
COMPLEX(dbl) :: eigr (:,:)
REAL(dbl), POINTER :: hg(:)
IF(tecfix) THEN
hg => gv%khgcutz_l(:,1)
hg => gkcutz_l(:,1)
ELSE
hg => gv%khg_l(:,1)
hg => gk_l(:,1)
END IF
CALL dforce1_d(c(:,ib), df, f(ib), hg, v)
CALL dforce2_d(f(ib), df, fnl(:,:,ib), gv%khg_l(:,1), &
gv%kgx_l(:,:,1), eigr, wsg, wnl)
CALL dforce2_d(f(ib), df, fnl(:,:,ib), gk_l(:,1), &
gkx_l(:,:,1), eigr, wsg, wnl)
return
end subroutine
@ -444,20 +443,21 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE dforce_kp( ik, ib, c0, cdesc, f, df, gv, v, fnlk, eigr, ps )
SUBROUTINE dforce_kp( ik, ib, c0, cdesc, f, df, v, fnlk, eigr, ps )
! (describe briefly what this routine does...)
! ----------------------------------------------
USE wave_types, ONLY: wave_descriptor
USE reciprocal_space_mesh, ONLY: gkmask_l
USE gvecw, ONLY: tecfix
USE reciprocal_space_mesh, ONLY: gkcutz_l, gkx_l, gk_l
IMPLICIT NONE
INTEGER, INTENT(IN) :: ib, ik
COMPLEX(dbl), INTENT(INOUT) :: c0(:,:,:)
type (wave_descriptor), INTENT(IN) :: cdesc
TYPE (recvecs) :: gv
COMPLEX(dbl) :: df(:)
REAL (dbl) :: v(:,:,:), f(:,:)
COMPLEX(dbl) :: fnlk(:,:,:)
@ -465,15 +465,14 @@
TYPE (pseudo) :: ps
IF(tecfix) THEN
CALL dforce1_kp( c0(:,ib,ik), df, f(ib,ik), gv%khgcutz_l(:,ik), v)
CALL dforce1_kp( c0(:,ib,ik), df, f(ib,ik), gkcutz_l(:,ik), v)
ELSE
CALL dforce1_kp( c0(:,ib,ik), df, f(ib,ik), gv%khg_l(:,ik), v)
CALL dforce1_kp( c0(:,ib,ik), df, f(ib,ik), gk_l(:,ik), v)
END IF
CALL dforce2_kp( f(ib,ik), df, fnlk(:,:,ib), gv%khg_l(:,ik), &
gv%kgx_l(:,:,ik), eigr, ps%wsg, ps%wnl(:,:,:,ik) )
CALL dforce2_kp( f(ib,ik), df, fnlk(:,:,ib), gk_l(:,ik), &
gkx_l(:,:,ik), eigr, ps%wsg, ps%wnl(:,:,:,ik) )
c0(:,ib,ik) = c0(:,ib,ik) * gv%kg_mask_l(:,ik)
df(:) = df(:) * gv%kg_mask_l(:,ik)
df(:) = df(:) * gkmask_l(:,ik)
RETURN
END SUBROUTINE dforce_kp
@ -594,7 +593,7 @@
END SUBROUTINE
SUBROUTINE dforce_all( ispin, c, cdesc, f, cgrad, gv, vpot, fnl, eigr, ps, ik)
SUBROUTINE dforce_all( ispin, c, cdesc, f, cgrad, vpot, fnl, eigr, ps, ik)
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
IMPLICIT NONE
@ -604,7 +603,6 @@
type (wave_descriptor), INTENT(IN) :: cdesc
TYPE (pseudo), INTENT(IN) :: ps
COMPLEX(dbl) :: eigr(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (projector) :: fnl(:)
REAL(dbl), INTENT(IN) :: vpot(:,:,:), f(:,:)
COMPLEX(dbl), INTENT(OUT) :: cgrad(:,:,:)
@ -627,18 +625,18 @@
ELSE
IF( .NOT. cdesc%gamma ) THEN
DO ib = 1, nb
CALL dforce(ikk, ib, c, cdesc, f, cgrad(:,ib,ikk), gv, vpot(:,:,:), fnl(ikk)%c, eigr, ps)
CALL dforce(ikk, ib, c, cdesc, f, cgrad(:,ib,ikk), vpot(:,:,:), fnl(ikk)%c, eigr, ps)
END DO
ELSE
! ... Process two states at the same time
DO ib = 1, ( nb - MOD(nb,2) ), 2
CALL dforce(ikk, ib, c, cdesc, f, cgrad(:,ib,ikk), cgrad(:,ib+1,1), &
gv, vpot(:,:,:), fnl(ikk)%r, eigr, ps)
vpot(:,:,:), fnl(ikk)%r, eigr, ps)
END DO
! ... Account for an odd number of states
ib = nb
IF( MOD(ib,2) .GT. 0 ) THEN
CALL dforce(ikk, ib, c, cdesc, f, cgrad(:,ib,ikk), gv, vpot(:,:,:), fnl(ikk)%r, eigr, ps)
CALL dforce(ikk, ib, c, cdesc, f, cgrad(:,ib,ikk), vpot(:,:,:), fnl(ikk)%r, eigr, ps)
END IF
END IF
END IF

27
CPV/fromscra.f90
View File

@ -33,7 +33,7 @@
! ----------------------------------------------
! BEGIN manual
SUBROUTINE from_scratch_fpmd( gv, kp, ps, rhoe, desc, cm, c0, cdesc, &
SUBROUTINE from_scratch_fpmd( kp, ps, rhoe, desc, cm, c0, cdesc, &
eigr, ei1, ei2, ei3, sfac, fi, ht, atoms, fnl, vpot, edft )
! (describe briefly what this routine does...)
@ -43,13 +43,13 @@
! ... declare modules
USE kinds
USE cp_types, ONLY: recvecs, pseudo
USE cp_types, ONLY: pseudo
USE wave_types, ONLY: wave_descriptor
USE atoms_type_module, ONLY: atoms_type
USE wave_functions, ONLY: gram, fixwave
USE wave_base, ONLY: wave_steepest
USE charge_density, ONLY: rhoofr
USE phase_factors_module, ONLY: strucf
USE phase_factors_module, ONLY: strucf, phfacs
USE cell_module, only: boxdimensions
USE electrons_module, ONLY: nspin, pmss
USE cp_electronic_mass, ONLY: emass
@ -67,6 +67,9 @@
USE charge_types, ONLY: charge_descriptor
USE time_step, ONLY: delt
USE runcp_module, ONLY: runcp_ncpp
use grid_dimensions, only: nr1, nr2, nr3
USE reciprocal_vectors, ONLY: mill_l
USE gvecp, ONLY: ngm
IMPLICIT NONE
@ -77,7 +80,6 @@
COMPLEX(dbl) :: ei1(:,:)
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
TYPE (recvecs) :: gv
TYPE (kpoints) :: kp
REAL(dbl) :: rhoe(:,:,:,:)
COMPLEX(dbl) :: sfac(:,:)
@ -108,11 +110,12 @@
atoms%for = 0.0d0
CALL strucf(sfac, atoms, eigr, ei1, ei2, ei3, gv)
edft%enl = nlrh_m(cm, cdesc, ttforce, atoms, fi, gv, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL rhoofr(gv, kp, cm, cdesc, fi, rhoe, desc, ht)
CALL phfacs( ei1, ei2, ei3, eigr, mill_l, atoms%taus, nr1, nr2, nr3, atoms%nat )
CALL strucf( sfac, ei1, ei2, ei3, mill_l, ngm )
edft%enl = nlrh_m(cm, cdesc, ttforce, atoms, fi, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL rhoofr(kp, cm, cdesc, fi, rhoe, desc, ht)
CALL vofrhos(ttprint, rhoe, desc, tfor, thdyn, ttforce, atoms, &
gv, kp, fnl, vpot, ps, cm, cdesc, fi, eigr, ei1, ei2, ei3, sfac, timepre, ht, edft)
kp, fnl, vpot, ps, cm, cdesc, fi, eigr, ei1, ei2, ei3, sfac, timepre, ht, edft)
CALL debug_energies( edft ) ! DEBUG
@ -120,7 +123,7 @@
IF( .NOT. force_pairing ) THEN
CALL runcp_ncpp( cm, cm, c0, cdesc, gv, kp, ps, vpot, eigr, &
CALL runcp_ncpp( cm, cm, c0, cdesc, kp, ps, vpot, eigr, &
fi, fnl, vdum, gam, cgam, fromscra = .TRUE. )
ELSE
@ -169,7 +172,8 @@ SUBROUTINE from_scratch_cp( sfac, eigr, ei1, ei2, ei3, bec, becdr, tfirst, eself
USE cpr_subroutines, ONLY: compute_stress, print_atomic_var, print_lambda, elec_fakekine
USE core, ONLY: nlcc_any
USE gvecw, ONLY: ngw
USE reciprocal_vectors, ONLY: gstart
USE reciprocal_vectors, ONLY: gstart, mill_l
USE gvecs, ONLY: ngs
USE wave_base, ONLY: wave_steepest
USE cvan, ONLY: nvb
USE ions_nose, ONLY: xnhp0, xnhpm, vnhp
@ -180,6 +184,7 @@ SUBROUTINE from_scratch_cp( sfac, eigr, ei1, ei2, ei3, bec, becdr, tfirst, eself
USE ensemble_dft, ONLY: tens, compute_entropy
USE runcp_module, ONLY: runcp_uspp
USE electrons_base, ONLY: f, nspin
USE phase_factors_module, ONLY: strucf
COMPLEX(kind=8) :: eigr(:,:), ei1(:,:), ei2(:,:), ei3(:,:)
COMPLEX(kind=8) :: eigrb(:,:)
@ -264,7 +269,7 @@ SUBROUTINE from_scratch_cp( sfac, eigr, ei1, ei2, ei3, bec, becdr, tfirst, eself
lambdam = lambda
!
call phfac( tau0, ei1, ei2, ei3, eigr )
call strucf( ei1, ei2, ei3, sfac )
call strucf( sfac, ei1, ei2, ei3, mill_l, ngs )
call formf( tfirst, eself )
IF( tefield ) THEN

243
CPV/gsmesh.f90
View File

@ -25,7 +25,6 @@
! ----------------------------------------------
! routines in this module:
! SUBROUTINE gmeshset(ngw,ng)
! SUBROUTINE gindexset(gv,nr1,nr2,nr3)
! INTEGER FUNCTION owner_of_gvec(ig)
! SUBROUTINE newg(gv,kp,htm1)
! ----------------------------------------------
@ -34,16 +33,35 @@
! ... declare included modules
USE kinds
USE cp_types
USE cell_base, ONLY: tpiba, tpiba2
IMPLICIT NONE
SAVE
REAL(dbl) :: b1(3) = 0.0d0
REAL(dbl) :: b2(3) = 0.0d0
REAL(dbl) :: b3(3) = 0.0d0
!! ... quantities related to k+G vectors
REAL(dbl), ALLOCATABLE, TARGET :: gkcutz_l(:,:) ! smooth cutoff factor index: G vector
! first index: G vector
! second index: k point
REAL(dbl), ALLOCATABLE :: gkmask_l(:,:) ! cutoff mask
REAL(dbl), ALLOCATABLE, TARGET :: gk_l(:,:) ! length squared of k+G
REAL(dbl), ALLOCATABLE :: gkx_l(:,:,:) ! components of k+G
! first index: G vector
! second index: x,y,z
! third index: k point
PRIVATE
PUBLIC :: recvecs_units, newg, gindexset, gmeshinfo, gindex_closeup
PUBLIC :: recvecs_units, newg, gmeshinfo, gindex_closeup
PUBLIC :: b1, b2, b3
PUBLIC :: gkmask_l, gkcutz_l, gkx_l, gk_l
! ... end of module-scope declarations
! ----------------------------------------------
@ -52,8 +70,6 @@
CONTAINS
!=----------------------------------------------------------------------------=!
! subroutines
SUBROUTINE recvecs_units( alat )
USE constants, ONLY: tpi
@ -65,9 +81,12 @@
RETURN
END SUBROUTINE
! ----------------------------------------------
SUBROUTINE gmeshinfo( )
SUBROUTINE gmeshinfo( )
! (describe briefly what this routine does...)
! ----------------------------------------------
@ -121,136 +140,48 @@
RETURN
END SUBROUTINE
END SUBROUTINE
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE gindexset(gv, b1, b2, b3)
! (describe briefly what this routine does...)
! ----------------------------------------------
! ... declare modules
USE mp_global, ONLY: mpime
USE io_global, ONLY: stdout
USE fft_types, ONLY: fft_dlay_descriptor
USE fft_base, ONLY: dfftp
USE stick_base, ONLY: stown => sticks_owner
USE reciprocal_vectors, only: &
ngw_g => ngwt, &
ngw_l => ngw , &
ngw_lx => ngwx, &
ng_g => ngmt, &
ng_l => ngm , &
ng_lx => ngmx, &
gstart, &
gzero, &
ig_l2g, &
mill_l, g, gx
IMPLICIT NONE
! ... declare subroutine arguments
TYPE (recvecs) gv
REAL(dbl), INTENT(IN) :: b1(3), b2(3), b3(3)
! ... declare other variables
INTEGER ig, ng, i, j, k, is
! ... end of declarations
! ----------------------------------------------
ng = 0
gv%gzero = gzero
gv%gstart = gstart
gv%bi1 = b1
gv%bi2 = b2
gv%bi3 = b3
gv%b1 = b1
gv%b2 = b2
gv%b3 = b3
gv%hg_l => g
gv%ig => ig_l2g( 1:ng_l )
gv%mill => mill_l
gv%gx_l => gx
IF( ng_l /= gv%ng_l ) THEN
WRITE( stdout,*) ' MSG: gv%ng_l = ',gv%ng_l,' ng = ', ng_l
CALL errore(' gindexset ',' inconsistent ng ', 2)
END IF
ng = ng_l
IF( gv%gzero .and. ( .not. ( gv%gstart == 2 ) ) ) THEN
CALL errore(' gindexset ',' gzero and gstart are inconsistent ', 3)
END IF
IF( gv%gstart < 1 .or. gv%gstart > 2 ) THEN
CALL errore(' gindexset ',' gstart out of range ', 4)
END IF
200 FORMAT(3I5)
RETURN
END SUBROUTINE gindexset
! ----------------------------------------------
! ----------------------------------------------
INTEGER FUNCTION owner_of_gvec(mill)
! (describe briefly what this routine does...)
! ----------------------------------------------
! ... declare modules
USE stick_base, ONLY: stown => sticks_owner
INTEGER FUNCTION owner_of_gvec(mill)
IMPLICIT NONE
USE stick_base, ONLY: stown => sticks_owner
! ... declare function arguments
INTEGER, INTENT(IN) :: mill(:)
IMPLICIT NONE
! ... end of declarations
! ----------------------------------------------
INTEGER, INTENT(IN) :: mill(:)
owner_of_gvec = stown( mill(1), mill(2) )
owner_of_gvec = stown( mill(1), mill(2) )
RETURN
END FUNCTION owner_of_gvec
RETURN
END FUNCTION owner_of_gvec
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE newg(gv,kp,htm1)
! this routine computes the squared modulus, Cartesian components and
! smooth cutoff masks of G vectors, from their Miller indices and the
! current cell shape. G vectors are expressed in units of 2*pi/alat.
SUBROUTINE newg( kp, htm1 )
! this routine computes the squared modulus, and Ciartesian components
! of G vectors, from their Miller indices and current cell shape.
! G vectors are expressed in units of 2*pi/alat.
! In generic-k-points calculations, squared modulus and components of
! k+G vectors are computed too, together with cutoff masks
! ----------------------------------------------
! ... declare modules
USE gvecw, ONLY: tecfix, gcfix, gcsig, gcutz, gcutw
USE cell_module, ONLY: alat
USE brillouin, ONLY: kpoints
USE reciprocal_vectors, only: &
ngw_g => ngwt, &
ngw_l => ngw , &
ngw_lx => ngwx, &
ng_g => ngmt, &
ng_l => ngm , &
ng_lx => ngmx, &
gstart, &
gzero, &
ig_l2g, &
mill_l
USE gvecw, ONLY: ngw, tecfix, gcfix, gcsig, gcutz, gcutw
USE gvecp, ONLY: ngm
USE reciprocal_vectors, only: g, gx, mill_l
USE cell_module, ONLY: alat
USE brillouin, ONLY: kpoints
IMPLICIT NONE
! ... declare subroutine argumentsz
TYPE (recvecs), INTENT(INOUT) :: gv
TYPE (kpoints), INTENT(IN) :: kp
REAL(dbl), INTENT(IN) :: htm1(3,3)
REAL(dbl) :: xk(3, SIZE(kp%xk, 2) )
@ -260,9 +191,9 @@
EXTERNAL erf, erfc
! ... declare other variables
INTEGER ig,i,j,k,ik
INTEGER isign
REAL(dbl) :: b1(3), b2(3), b3(3)
INTEGER :: ig, i, j, k, ik, nkp
INTEGER :: isign
REAL(dbl) :: gmax
! ... end of declarations
! ----------------------------------------------
@ -270,67 +201,57 @@
b1 = htm1(:,1)
b2 = htm1(:,2)
b3 = htm1(:,3)
gv%b1 = htm1(:,1)
gv%b2 = htm1(:,2)
gv%b3 = htm1(:,3)
DO ig = 1, gv%ng_l
call gcal( alat, b1, b2, b3, gmax )
i = gv%mill(1,ig)
j = gv%mill(2,ig)
k = gv%mill(3,ig)
! ... compute components of G
gv%gx_l(1,ig) = b1(1)*i + b2(1)*j + b3(1)*k
gv%gx_l(2,ig) = b1(2)*i + b2(2)*j + b3(2)*k
gv%gx_l(3,ig) = b1(3)*i + b2(3)*j + b3(3)*k
gv%gx_l(1,ig) = gv%gx_l(1,ig) * alat
gv%gx_l(2,ig) = gv%gx_l(2,ig) * alat
gv%gx_l(3,ig) = gv%gx_l(3,ig) * alat
! ... compute squared length
gv%hg_l(ig) = gv%gx_l(1,ig)**2 + gv%gx_l(2,ig)**2 + gv%gx_l(3,ig)**2
END DO
nkp = kp%nkpt
IF( .NOT. ALLOCATED( gk_l ) ) ALLOCATE( gk_l( ngw, nkp ) )
IF( .NOT. ALLOCATED( gkx_l ) ) ALLOCATE( gkx_l( 3, ngw, nkp ) )
IF( .NOT. ALLOCATED( gkcutz_l ) ) ALLOCATE( gkcutz_l( ngw, nkp ) )
IF( .NOT. ALLOCATED( gkmask_l ) ) ALLOCATE( gkmask_l ( ngw, nkp ) )
IF( SIZE( gk_l ) /= ( ngw * nkp ) ) THEN
DEALLOCATE( gk_l )
DEALLOCATE( gkx_l )
DEALLOCATE( gkcutz_l )
DEALLOCATE( gkmask_l )
ALLOCATE( gk_l( ngw, nkp ) )
ALLOCATE( gkx_l( 3, ngw, nkp ) )
ALLOCATE( gkcutz_l( ngw, nkp ) )
ALLOCATE( gkmask_l( ngw, nkp ) )
END IF
IF( kp%scheme == 'gamma' ) THEN
DO ik = 1, kp%nkpt
gv%khg_l( 1:SIZE(gv%khg_l,1) ,ik) = gv%hg_l( 1:SIZE(gv%khg_l,1) )
gv%kgx_l( 1:3, 1:SIZE(gv%khg_l,1), ik) = gv%gx_l( 1:3, 1:SIZE(gv%khg_l,1) )
gk_l( 1:ngw ,ik) = g( 1:ngw )
gkx_l( 1:3, 1:ngw, ik) = gx( 1:3, 1:ngw )
END DO
ELSE
DO ik = 1, kp%nkpt
! ... Bring k-points in the unscaled reciprocal space
xk(:,ik) = kp%xk(:,ik)
!xk(1,ik) = htm1(1,1)*kp%xk(1,ik) + htm1(1,2)*kp%xk(2,ik) + htm1(1,3)*kp%xk(3,ik)
!xk(1,ik) = xk(1,ik) * alat
!xk(2,ik) = htm1(2,1)*kp%xk(1,ik) + htm1(2,2)*kp%xk(2,ik) + htm1(2,3)*kp%xk(3,ik)
!xk(2,ik) = xk(2,ik) * alat
!xk(3,ik) = htm1(3,1)*kp%xk(1,ik) + htm1(3,2)*kp%xk(2,ik) + htm1(3,3)*kp%xk(3,ik)
!xk(3,ik) = xk(3,ik) * alat
!WRITE( stdout,*) alat, xk(1,ik), xk(2,ik), xk(3,ik)
! ... compute components of G
DO ig = 1, SIZE( gv%khg_l, 1 )
DO ig = 1, ngw
! ... compute components of G+k
gv%kgx_l(1,ig,ik) = gv%gx_l(1,ig) + xk(1,ik)
gv%kgx_l(2,ig,ik) = gv%gx_l(2,ig) + xk(2,ik)
gv%kgx_l(3,ig,ik) = gv%gx_l(3,ig) + xk(3,ik)
gkx_l(1,ig,ik) = gx(1,ig) + xk(1,ik)
gkx_l(2,ig,ik) = gx(2,ig) + xk(2,ik)
gkx_l(3,ig,ik) = gx(3,ig) + xk(3,ik)
! ... compute squared length
gv%khg_l(ig,ik) = gv%kgx_l(1,ig,ik)**2 + gv%kgx_l(2,ig,ik)**2 + gv%kgx_l(3,ig,ik)**2
! ... compute cutoff mask for k+G
IF( gv%khg_l(ig, ik) .LE. gcutw) THEN
gv%kg_mask_l(ig, ik) = 1.0d0
gk_l(ig,ik) = gkx_l(1,ig,ik)**2 + gkx_l(2,ig,ik)**2 + gkx_l(3,ig,ik)**2
! compute cutoff mask for k+G
IF( gk_l(ig, ik) <= gcutw ) THEN
gkmask_l(ig, ik) = 1.0d0
ELSE
gv%kg_mask_l(ig, ik) = 0.0d0
gkmask_l(ig, ik) = 0.0d0
END IF
END DO
@ -339,18 +260,18 @@
IF(tecfix) THEN
DO ik = 1, kp%nkpt
DO ig = 1, SIZE( gv%khgcutz_l, 1 )
DO ig = 1, ngw
! ... compute smooth cutoff G+k vectors
gv%khgcutz_l(ig,ik) = erf((gv%khg_l(ig,ik) - gcfix)/gcsig)
gv%khgcutz_l(ig,ik) = gv%khg_l(ig,ik) + gcutz * ( 1.0d0 + gv%khgcutz_l(ig,ik))
gkcutz_l(ig,ik) = erf( ( gk_l(ig,ik) - gcfix ) / gcsig )
gkcutz_l(ig,ik) = gk_l(ig,ik) + gcutz * ( 1.0d0 + gkcutz_l(ig,ik) )
END DO
END DO
END IF
RETURN
END SUBROUTINE newg
!=----------------------------------------------------------------------------=!
SUBROUTINE gindex_closeup
IMPLICIT NONE

9
CPV/guess.f90
View File

@ -28,7 +28,7 @@
! routines in this module:
! SUBROUTINE guess_setup(tguess_inp)
! SUBROUTINE guessc0(tk,c0,cm)
! SUBROUTINE guessrho(rho,cm,c0,occ,gv,ht)
! SUBROUTINE guessrho(rho,cm,c0,occ,ht)
! SUBROUTINE ucalc_kp(a,b,ngw,gzero,n,lambda)
! SUBROUTINE ucalc(a,b,ngw,gzero,n,lambda)
! ----------------------------------------------
@ -248,7 +248,7 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE guessrho(rho, desc, cm, c0, cdesc, occ, gv, kp, ht )
SUBROUTINE guessrho(rho, desc, cm, c0, cdesc, occ, kp, ht )
! (describe briefly what this routine does...)
! ----------------------------------------------
@ -265,7 +265,6 @@
TYPE (charge_descriptor), INTENT(IN) :: desc
COMPLEX(dbl), INTENT(IN) :: c0(:,:,:,:), cm(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
TYPE (boxdimensions), INTENT(IN) :: ht
REAL(dbl), INTENT(IN) :: occ(:,:,:)
@ -286,12 +285,12 @@
IF( tfirst ) THEN
ALLOCATE( rho_save( nx, ny, nz, nspin ) )
CALL rhoofr(gv, kp, cm, cdesc, occ, rho_save, desc, ht)
CALL rhoofr(kp, cm, cdesc, occ, rho_save, desc, ht)
tfirst = .FALSE.
END IF
ALLOCATE( rho0( nx, ny, nz, nspin ) )
CALL rhoofr(gv, kp, c0, cdesc, occ, rho0, desc, ht)
CALL rhoofr(kp, c0, cdesc, occ, rho0, desc, ht)
rho = 2.0d0 * rho0 - rho_save

335
CPV/init.f90
View File

@ -31,45 +31,35 @@
! BEGIN manual
SUBROUTINE init0s(gv, kp, ps, atoms_m, atoms_0, atoms_p, wfill, &
wempt, ht_m, ht, fnl, eigr, ei1, ei2, ei3, nspin)
SUBROUTINE init0s(kp, ps, atoms_m, atoms_0, atoms_p, wfill, &
wempt, ht_m, ht, eigr, ei1, ei2, ei3)
! this routine handles data initialization
! ----------------------------------------------
! END manual
! ... declare modules
use mp_global, only: nproc
USE parameters, ONLY: nspinx
USE phase_factors_module, ONLY: strucf
USE cp_types
USE atoms_type_module, ONLY: atoms_type
USE time_step, ONLY: delt
USE cell_module, ONLY: cell_init, get_lattice_vectors
USE cell_base, ONLY: omega, alat
USE electrons_module, ONLY: electron_mass_init, band_init, n_emp, bmeshset
USE electrons_base, ONLY: nupdwn
USE reciprocal_space_mesh, ONLY: newg
USE reciprocal_vectors, ONLY: ngwt, gstart, gzero, ngm, ngmt, ngw
USE pseudopotential, ONLY: formf, nsanl, ngh, pseudopotential_init
USE ions_module, ONLY: atoms_init
USE ions_base, ONLY: nsp, na, nat
USE pseudo_projector, ONLY: allocate_projector, projector
USE diis, ONLY: allocate_diis, delt_diis
USE cell_module, only: boxdimensions
USE mp_global, ONLY: mpime, root
USE brillouin, ONLY: kpoints
USE wave_types, ONLY: wave_descriptor, wave_descriptor_init, &
wave_descriptor_info
USE descriptors_module, ONLY: get_local_dims, get_global_dims
USE control_flags, ONLY: nbeg, tbeg, timing, t_diis, iprsta
USE input_parameters, ONLY: rd_ht
USE turbo, ONLY: tturbo, allocate_turbo
USE ions_base, ONLY: tau_srt, tau_units, ind_srt, if_pos, atm
USE fft_base, ONLY: dfftp
USE grid_dimensions, ONLY: nr1, nr2, nr3
USE problem_size, ONLY: cpsizes
USE reciprocal_space_mesh, ONLY : gindexset
USE cp_types, ONLY: pseudo
use mp_global, only: nproc
USE phase_factors_module, ONLY: strucf, phfacs
USE atoms_type_module, ONLY: atoms_type
USE cell_module, ONLY: cell_init, get_lattice_vectors
USE cell_module, only: boxdimensions
USE electrons_module, ONLY: electron_mass_init, band_init, n_emp
USE electrons_base, ONLY: nspin, nupdwn
USE reciprocal_vectors, ONLY: ngwt, gstart, gzero, ngm, ngmt, ngw, mill_l
USE pseudopotential, ONLY: formf, nsanl, ngh, pseudopotential_init
USE ions_base, ONLY: nsp, na, nat
USE ions_module, ONLY: atoms_init
USE brillouin, ONLY: kpoints
USE wave_types, ONLY: wave_descriptor, wave_descriptor_init, wave_descriptor_info
USE control_flags, ONLY: nbeg, tbeg, timing, iprsta
USE turbo, ONLY: tturbo, allocate_turbo
USE ions_base, ONLY: ind_srt, if_pos, atm
USE fft_base, ONLY: dfftp
USE grid_dimensions, ONLY: nr1, nr2, nr3
USE ions_positions, ONLY: taus
USE reciprocal_space_mesh, ONLY: newg
IMPLICIT NONE
@ -78,112 +68,76 @@
TYPE (atoms_type) :: atoms_0, atoms_p, atoms_m
TYPE (wave_descriptor) :: wfill, wempt
TYPE (pseudo) :: ps
TYPE (projector) :: fnl(:,:)
COMPLEX(dbl) :: eigr(:,:)
COMPLEX(dbl) :: ei1(:,:)
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
TYPE (recvecs) :: gv
TYPE (kpoints) :: kp
TYPE (boxdimensions) :: ht_m, ht
INTEGER :: nspin
! ... declare other variables
REAL(dbl) :: s1, s2, s3, s4, s5
REAL(dbl) :: a1(3), a2(3), a3(3)
real(dbl) :: b1(3), b2(3), b3(3)
INTEGER :: i, ispin, isym
LOGICAL :: tk
COMPLEX(dbl), ALLOCATABLE :: sfac(:,:)
REAL(dbl) :: rat1, rat2, rat3
INTEGER :: neupdwn( nspinx )
INTEGER :: neupdwn( nspin )
! end of declarations
! ----------------------------------------------
s1 = cclock()
call get_lattice_vectors( a1, a2, a3 )
call recips( a1, a2, a3, b1, b2, b3 )
! ... arrange for reciprocal lattice vectors
tk = .NOT. ( kp%scheme == 'gamma' )
CALL allocate_recvecs(gv, ngm, ngmt, ngw, ngwt, tk, kp%nkpt)
! ... Initialize cell variables
CALL gindexset( gv, b1, b2, b3 )
CALL cell_init( ht, a1, a2, a3 )
CALL cell_init( ht_m, a1, a2, a3 )
! ... set the bands mesh
!
CALL bmeshset( )
! ... compute reciprocal lattice vectors
CALL newg(kp, ht%m1)
CALL cpsizes( nproc )
!
CALL cpflush( )
! ... initialize atomic configuration (should be called after metric_init)
CALL atoms_init( atoms_m, atoms_0, atoms_p, taus, ind_srt, if_pos, atm, ht%hmat )
! ... compute structure factors
ALLOCATE( sfac( ngm, nsp ) )
CALL phfacs( ei1, ei2, ei3, eigr, mill_l, atoms_0%taus, nr1, nr2, nr3, atoms_0%nat )
CALL strucf( sfac, ei1, ei2, ei3, mill_l, ngm )
! ... Allocate + Initialize pseudopotentials
CALL pseudopotential_init(ps, na, nsp, gv, kp)
s2 = cclock()
! ... Initialize simulation cell
! ... if tbeg = .TRUE. cell parameters have been specified in the input file
! ... (see card 'TBEG')
! ... But the g-space grid is genereted according to celldm
IF( tbeg ) THEN
CALL cell_init( ht, rd_ht )
CALL cell_init( ht_m, rd_ht )
ELSE
CALL cell_init( ht, a1, a2, a3 )
CALL cell_init( ht_m, a1, a2, a3 )
END IF
! ... initialize atomic configuration (should be called after metric_init)
CALL atoms_init( atoms_m, atoms_0, atoms_p, tau_srt, ind_srt, if_pos, atm, tau_units, alat, ht )
! ... compute reciprocal lattice vectors
CALL newg(gv, kp, ht%m1)
s3 = cclock()
! ... compute structure factors
ALLOCATE( sfac( nsp, ngm ) )
CALL strucf( sfac, atoms_0, eigr, ei1, ei2, ei3, gv )
! ... Allocate + Initialize pseudopotentials
CALL pseudopotential_init(ps, na, nsp, kp)
s4 = cclock()
! ... compute local form factors
CALL formf(ht, gv, kp, ps)
! ... compute local form factors
CALL formf(ht, kp, ps)
s5 = cclock()
IF(ionode) THEN
WRITE( stdout,'(/," ggen (sec) : ",F8.3)') (s2-s1)
WRITE( stdout,'( " newg (sec) : ",F8.3)') (s3-s2)
WRITE( stdout,'( " strucf (sec) : ",F8.3)') (s4-s3)
WRITE( stdout,'( " formf (sec) : ",F8.3)') (s5-s4)
END IF
tk = .NOT. ( kp%scheme == 'gamma' )
isym = 0
IF( tk ) isym = 1
! empty states, always same number of spin up and down states
neupdwn( 1:nspin ) = n_emp
CALL wave_descriptor_init( wfill, gv%ngw_l, gv%ngw_g, nupdwn, nupdwn, &
kp%nkpt, kp%nkpt, nspin, isym, gv%gzero )
CALL wave_descriptor_init( wempt, gv%ngw_l, gv%ngw_g, neupdwn, neupdwn, &
kp%nkpt, kp%nkpt, nspin, isym, gv%gzero )
CALL wave_descriptor_init( wfill, ngw, ngwt, nupdwn, nupdwn, &
kp%nkpt, kp%nkpt, nspin, isym, gzero )
CALL wave_descriptor_init( wempt, ngw, ngwt, neupdwn, neupdwn, &
kp%nkpt, kp%nkpt, nspin, isym, gzero )
IF( iprsta > 2 ) THEN
CALL wave_descriptor_info( wfill, 'wfill', stdout )
CALL wave_descriptor_info( wempt, 'wempt', stdout )
END IF
! ... if tturbo=.TRUE. some data is stored in memory instead of being
! ... recalculated (see card 'TURBO')
! ... if tturbo=.TRUE. some data is stored in memory instead of being
! ... recalculated (see card 'TURBO')
IF( tturbo ) THEN
CALL allocate_turbo( dfftp%nr1x, dfftp%nr2x, dfftp%npl )
ENDIF
@ -198,87 +152,57 @@
! BEGIN manual
SUBROUTINE init1s(gv, kp, ps, atoms_m, atoms_0, atoms_p, cm, c0, wfill, &
ce, wempt, ht_m, ht, fnl, eigr, occ)
SUBROUTINE init1s(kp, ps, cm, c0, wfill, ce, wempt, fnl, eigr, occ)
! this routine handles data initialization
! ----------------------------------------------
! END manual
! ... declare modules
USE phase_factors_module, ONLY: strucf
USE wave_init, ONLY: pw_atomic_init
USE cp_types
USE atoms_type_module, ONLY: atoms_type
USE time_step, ONLY: delt
USE cell_module, ONLY: cell_init, get_lattice_vectors, alat
USE electrons_module, ONLY: electron_mass_init, band_init, nbnd
USE reciprocal_space_mesh, ONLY: newg
USE reciprocal_vectors, ONLY: ngwt, gstart, gzero, ngm, ngmt, ngw
USE pseudopotential, ONLY: formf, nsanl, ngh, pseudopotential_init
USE ions_module, ONLY: atoms_init
USE ions_base, ONLY: nsp, na, nat
USE pseudo_projector, ONLY: allocate_projector, projector
USE diis, ONLY: allocate_diis, delt_diis
USE cell_module, only: boxdimensions
USE mp_global, ONLY: mpime, root
USE brillouin, ONLY: kpoints
USE wave_types, ONLY: wave_descriptor
USE descriptors_module, ONLY: get_local_dims, get_global_dims
USE control_flags, ONLY: nbeg, tbeg, timing, t_diis
USE cell_module, ONLY: alat
USE electrons_module, ONLY: electron_mass_init, band_init, nbnd
USE reciprocal_vectors, ONLY: g
USE gvecw, ONLY: ngw
USE pseudopotential, ONLY: nsanl, ngh
USE pseudo_projector, ONLY: allocate_projector, projector
USE diis, ONLY: allocate_diis
USE brillouin, ONLY: kpoints
USE wave_types, ONLY: wave_descriptor
USE wave_init, ONLY: pw_atomic_init
USE control_flags, ONLY: nbeg, t_diis
IMPLICIT NONE
! ... declare subroutine arguments
TYPE (atoms_type) :: atoms_0, atoms_p, atoms_m
COMPLEX(dbl) :: cm(:,:,:,:), c0(:,:,:,:), ce(:,:,:,:)
TYPE (wave_descriptor) :: wfill, wempt
TYPE (pseudo) :: ps
REAL(dbl) :: occ(:,:,:)
TYPE (projector) :: fnl(:,:)
COMPLEX(dbl) :: eigr(:,:)
TYPE (recvecs) :: gv
TYPE (kpoints) :: kp
TYPE (boxdimensions) :: ht_m, ht
! ... declare other variables
REAL(dbl) :: s1, s2, s3, s4, s5
REAL(dbl) :: a1(3), a2(3), a3(3)
INTEGER :: i
LOGICAL :: tk
REAL(dbl), ALLOCATABLE :: hg_g(:) ! squared length
INTEGER, ALLOCATABLE :: mill(:,:) ! Miller index, axis x, y, z
COMPLEX(dbl), ALLOCATABLE :: sfac(:,:)
! end of declarations
! ----------------------------------------------
s1 = cclock()
! ... initialize bands
CALL band_init( occ )
! ... initialize wave functions
CALL pw_atomic_init(nbeg, cm, c0, wfill, ce, wempt, gv, kp, eigr )
CALL pw_atomic_init(nbeg, cm, c0, wfill, ce, wempt, kp, eigr )
! ... initialize the electronic fictitious mass and time step for
! ... electron dynamics. Note that for 'diis' the electronic time step is
! ... different from that of ions (delt), this is because in the diis
! ... the time step for electrons is simply a convergence parameter.
IF (t_diis) THEN
CALL electron_mass_init(alat, gv%hg_l, gv%ngw_l)
ELSE
CALL electron_mass_init(alat, gv%hg_l, gv%ngw_l)
END IF
! ... initialize the electronic fictitious mass and time step for
! ... electron dynamics.
CALL electron_mass_init(alat, g, ngw)
! ... initialize nonlocal pseudopotentials coefficients
! ... initialize nonlocal pseudopotentials coefficients
CALL allocate_projector(fnl, nsanl, nbnd, ngh, kp%gamma_only)
IF(t_diis) THEN
! ... arrange for DIIS minimization
CALL allocate_diis(gv%ngw_l, nbnd, kp%nkpt)
! ... arrange for DIIS minimization
CALL allocate_diis(ngw, nbnd, kp%nkpt)
END IF
CALL cpflush ! flush output streams
@ -330,7 +254,7 @@
USE control_flags, ONLY: tdipole
USE berry_phase, ONLY: berry_setup
USE real_space_mesh, ONLY: realspace_procgrid_init
USE bands_mesh, ONLY: bands_procgrid_init
USE electrons_module, ONLY: bmeshset
implicit none
!
@ -345,11 +269,10 @@
3X,'------------------------------------' )
END IF
!
! ... Initialize processor grid for parallel linear algebra
! used electronic states lagrange multiplier matrixes
! ... Initialize bands indexes for parallel linear algebra
! ... (distribute bands to processors)
!
CALL bands_procgrid_init( )
CALL bmeshset( )
!
! ... Initialize (global) real and compute global reciprocal dimensions
@ -490,38 +413,48 @@
!-----------------------------------------------------------------------
subroutine init ( ibrav, ndr, nbeg, tfirst, tau0, taus )
subroutine init_geometry ( )
!-----------------------------------------------------------------------
!
use control_flags, only: iprint, thdyn
use io_global, only: stdout
use ions_base, only: na, nsp, natx
use cell_base, only: ainv, a1, a2, a3, r_to_s, s_to_r
use cell_base, only: hold, h
use cp_restart, only: cp_read_cell
USE input_parameters, ONLY: trd_ht
use control_flags, only: iprint, thdyn, ndr, nbeg
use io_global, only: stdout
use ions_base, only: na, nsp, nat, natx, tau_srt
use cell_base, only: a1, a2, a3, r_to_s
use cell_base, only: ibrav, ainv, h, hold, tcell_base_init
USE ions_positions, ONLY: tau0, taus
use cp_restart, only: cp_read_cell
implicit none
! input/output
integer ibrav, ndr, nbeg
logical tfirst
real(kind=8) tau0(3,natx), taus(3,natx)
! local
integer i, j
real(kind=8) gvel(3,3)
real(kind=8) xnhh0(3,3),xnhhm(3,3),vnhh(3,3),velh(3,3)
!
!
! taus = scaled, tau0 = alat units
!
CALL r_to_s( tau0, taus, na, nsp, ainv )
!
if( nbeg >= 0 ) then
!
! local
!
integer :: i, j
real(kind=8) :: gvel(3,3)
real(kind=8) :: xnhh0(3,3), xnhhm(3,3), vnhh(3,3), velh(3,3)
IF( .NOT. tcell_base_init ) &
CALL errore( ' init_geometry ', ' cell_base_init has not been call yet! ', 1 )
!
! Scale positions that have been read from standard input
! according to the cell given in the standard input too
! taus = scaled, tau0 = atomic units
!
tau0 = 0.0d0
taus = 0.0d0
tau0 ( 1:3 , 1:nat ) = tau_srt ( 1:3 , 1:nat )
CALL r_to_s( tau0, taus, na, nsp, ainv )
!
! if trd_ht = .true. the geometry is given in the standard input even if
! we are restarting a previous run
!
if( ( nbeg >= 0 ) .and. ( .not. trd_ht ) ) then
!
! read only h and hold from file ndr
!
CALL cp_read_cell( ndr, ' ', .TRUE., h, hold, velh, gvel, xnhh0, xnhhm, vnhh )
h = TRANSPOSE( h )
@ -535,11 +468,9 @@
WRITE( stdout,*)
else
!
! with variable-cell we use h to describe the cell
! geometry is set to the cell parameters read from stdin ( a1, a2, a3 )
!
do i = 1, 3
h(i,1) = a1(i)
h(i,2) = a2(i)
@ -549,56 +480,52 @@
hold = h
end if
!
! ==============================================================
! ==== generate true g-space ====
! ==============================================================
!
call newinit( ibrav )
call newinit( h )
!
!
344 format(' ibrav = ',i4,' cell parameters ',/)
345 format(3(4x,f10.5))
return
end
end subroutine init_geometry
!-----------------------------------------------------------------------
subroutine newinit(ibrav)
!-----------------------------------------------------------------------
! re-initialization of lattice parameters and g-space vectors.
! Note that direct and reciprocal lattice primitive vectors
! a1,a2,a3, ainv, and corresponding quantities for small boxes
! are recalculated according to the value of cell parameter h
!
use control_flags, only: iprsta
use io_global, only: stdout
use cell_base, only: h, a1, a2, a3, omega, alat, cell_base_reinit
!
subroutine newinit( h )
!
! re-initialization of lattice parameters and g-space vectors.
! Note that direct and reciprocal lattice primitive vectors
! a1,a2,a3, ainv, and corresponding quantities for small boxes
! are recalculated according to the value of cell parameter h
!
use cell_base, only: a1, a2, a3, omega, alat, cell_base_reinit
!
implicit none
integer ibrav
!
! local
integer :: i, j
!
real(kind=8) :: h(3,3)
! local
!
real(kind=8) :: gmax, b1(3), b2(3), b3(3)
!
!
!
! re-initialize the cell base module with the new geometry
!
CALL cell_base_reinit( TRANSPOSE( h ) )
!
!
call recips( a1, a2, a3, b1, b2, b3 )
b1 = b1 * alat
b2 = b2 * alat
b3 = b3 * alat
call gcal( b1, b2, b3, gmax )
!
! ==============================================================
! generation of little box g-vectors
! ==============================================================
!
call gcal( alat, b1, b2, b3, gmax )
!
! generation of little box g-vectors
!
call newgb( a1, a2, a3, omega, alat )
!
!
return
end
end subroutine newinit

13
CPV/input.f90
View File

@ -835,8 +835,6 @@
USE potentials, ONLY: potential_init
USE kohn_sham_states, ONLY: ks_states_init
USE electrons_module, ONLY: electrons_setup
USE ions_positions, ONLY: tau0
USE ions_base, ONLY: tau_srt, tions_base_init
USE electrons_base, ONLY: electrons_base_initval
USE ensemble_dft, ONLY: ensemble_initval
@ -956,6 +954,8 @@
CALL electrons_base_initval( nelec, nelup, neldw, nbnd, nspin, occupations, f_inp )
CALL electrons_setup( empty_states_nbnd, emass, emass_cutoff, nkstot )
CALL ensemble_initval &
( occupations, n_inner, fermi_energy, rotmass, occmass, rotation_damping, &
@ -967,11 +967,7 @@
IF( program_name == 'FPMD' ) THEN
CALL pseudopotential_setup( ntyp, tpstab_inp, pstab_size_inp, ion_radius )
CALL electrons_setup( empty_states_nbnd, emass, emass_cutoff, nkstot )
!
o_diis_inp = .TRUE.
oqnr_diis_inp = .TRUE.
tolene_inp = etot_conv_thr
@ -987,11 +983,6 @@
CALL guess_setup( diis_chguess )
CALL charge_mix_setup(diis_achmix, diis_g0chmix, diis_nchmix, diis_g1chmix)
ELSE
tau0 = 0.0d0
tau0 ( 1:3 , 1:nat ) = tau_srt ( 1:3 , 1:nat )
END IF

49
CPV/ions.f90
View File

@ -469,47 +469,25 @@
return
END SUBROUTINE ions_setup
! BEGIN manual -------------------------------------------------------------
! -------------------------------------------------------------
SUBROUTINE atoms_init(atoms_m, atoms_0, atoms_p, tau_srt, ind_srt, &
if_pos, atml, aunits, alat, ht_0)
SUBROUTINE atoms_init(atoms_m, atoms_0, atoms_p, stau, ind_srt, if_pos, atml, h)
! Allocate and fill the three atoms structure using position an cell
! read from stdin. These initial values could be overwritten in module
! restart.f90 with the values read from restart file.
! This subroutine is called from modules init.f90
! --------------------------------------------------------------------------
! END manual ---------------------------------------------------------------
USE cell_module, ONLY: boxdimensions, s_to_r, r_to_s
USE constants, ONLY: angstrom_au
! Allocate and fill the three atoms structure using scaled position an cell
IMPLICIT NONE
!
TYPE (atoms_type) :: atoms_0, atoms_p, atoms_m
TYPE (boxdimensions) :: ht_0
REAL(dbl), INTENT(IN) :: tau_srt(:,:)
CHARACTER(LEN=*), INTENT(IN) :: aunits
REAL(dbl), INTENT(IN) :: h( 3, 3 )
REAL(dbl), INTENT(IN) :: stau(:,:)
CHARACTER(LEN=3), INTENT(IN) :: atml(:)
REAL(dbl), INTENT(IN) :: alat
INTEGER, INTENT(IN) :: ind_srt( : )
INTEGER, INTENT(IN) :: if_pos( :, : )
REAL(dbl), ALLOCATABLE :: stau(:,:)
LOGICAL , ALLOCATABLE :: ismb(:,:)
REAL(dbl) :: p( 3 )
INTEGER :: ia, isa
LOGICAL, SAVE :: atoms_init_first = .true.
LOGICAL, ALLOCATABLE :: ismb(:,:)
INTEGER :: ia, isa
IF( .NOT. atoms_init_first ) THEN
CALL errore(' atoms_init ',' atoms objects already initialized ', 1)
END IF
ALLOCATE( stau( 3, nat ), ismb( 3, nat ) )
DO ia = 1, nat
CALL r_to_s( tau_srt(:,ia), stau(:,ia), ht_0)
END DO
ALLOCATE( ismb( 3, nat ) )
ismb = .TRUE.
DO isa = 1, nat
@ -517,18 +495,15 @@
IF( if_pos( 1, ia ) == 0 ) ismb( 1, isa ) = .FALSE.
IF( if_pos( 2, ia ) == 0 ) ismb( 2, isa ) = .FALSE.
IF( if_pos( 3, ia ) == 0 ) ismb( 3, isa ) = .FALSE.
! WRITE(6,*) 'DEBUG atoms_init ', ismb(1:3,isa)
END DO
CALL atoms_type_init(atoms_m, stau, ismb, atml, pmass, na, nsp, ht_0%hmat)
CALL atoms_type_init(atoms_0, stau, ismb, atml, pmass, na, nsp, ht_0%hmat)
CALL atoms_type_init(atoms_p, stau, ismb, atml, pmass, na, nsp, ht_0%hmat)
CALL atoms_type_init(atoms_m, stau, ismb, atml, pmass, na, nsp, h)
CALL atoms_type_init(atoms_0, stau, ismb, atml, pmass, na, nsp, h)
CALL atoms_type_init(atoms_p, stau, ismb, atml, pmass, na, nsp, h)
CALL print_scaled_positions( atoms_0, stdout, 'from standard input')
DEALLOCATE( stau, ismb )
atoms_init_first = .FALSE.
DEALLOCATE( ismb )
RETURN
END SUBROUTINE atoms_init

18
CPV/ksstates.f90
View File

@ -216,7 +216,7 @@
! ----------------------------------------------
! BEGIN manual
SUBROUTINE ks_states(cf, wfill, ce, wempt, occ, gv, kp, ps, vpot, eigr, fnl)
SUBROUTINE ks_states(cf, wfill, ce, wempt, occ, kp, ps, vpot, eigr, fnl)
! (describe briefly what this routine does...)
! ----------------------------------------------
@ -242,7 +242,6 @@
TYPE (wave_descriptor), INTENT(IN) :: wfill, wempt
TYPE (pseudo), INTENT(IN) :: ps
COMPLEX(dbl) :: eigr(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
REAL(dbl), INTENT(IN) :: occ(:,:,:)
TYPE (projector) :: fnl(:,:)
@ -287,7 +286,7 @@
ALLOCATE( eforce( ngw, nb_l, nk ))
CALL dforce_all( ispin, cf(:,:,:,ispin_wfc), wfill, occ(:,:,ispin), eforce, &
gv, vpot(:,:,:,ispin), fnl(:,ispin), eigr, ps)
vpot(:,:,:,ispin), fnl(:,ispin), eigr, ps)
CALL kohn_sham( ispin, cf(:,:,:,ispin_wfc), wfill, eforce, kp )
@ -309,7 +308,7 @@
ALLOCATE( eforce( ngw, nb_l, nk ))
CALL dforce_all( ispin, ce(:,:,:,ispin), wempt, fi, eforce, gv, vpot(:,:,:,ispin), &
CALL dforce_all( ispin, ce(:,:,:,ispin), wempt, fi, eforce, vpot(:,:,:,ispin), &
fnl(:,ispin), eigr, ps)
CALL kohn_sham( ispin, ce(:,:,:,ispin), wempt, eforce, kp )
@ -382,7 +381,7 @@
! ----------------------------------------------
! BEGIN manual
SUBROUTINE ks_states_force_pairing(cf, wfill, ce, wempt, occ, gv, kp, ps, vpot, eigr, fnl)
SUBROUTINE ks_states_force_pairing(cf, wfill, ce, wempt, occ, kp, ps, vpot, eigr, fnl)
! (describe briefly what this routine does...)
! ----------------------------------------------
@ -409,7 +408,6 @@
TYPE (wave_descriptor), INTENT(IN) :: wfill, wempt
TYPE (pseudo), INTENT(IN) :: ps
COMPLEX(dbl) :: eigr(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
REAL(dbl), INTENT(IN) :: occ(:,:,:)
TYPE (projector) :: fnl(:,:)
@ -453,9 +451,9 @@
ALLOCATE( eforce( ngw, nb, 1, 2 ) )
CALL dforce_all( 1, cf(:,:,:,1), wfill, occ(:,:,1), eforce(:,:,:,1), &
gv, vpot(:,:,:,1), fnl(:,1), eigr, ps)
vpot(:,:,:,1), fnl(:,1), eigr, ps)
CALL dforce_all( 2, cf(:,:,:,1), wfill, occ(:,:,2), eforce(:,:,:,2), &
gv, vpot(:,:,:,2), fnl(:,2), eigr, ps)
vpot(:,:,:,2), fnl(:,2), eigr, ps)
DO i = 1, nupdwn(2)
eforce(:,i,1,1) = occ(i,1,1) * eforce(:,i,1,1) + occ(i,1,2) * eforce(:,i,1,2)
@ -484,12 +482,12 @@
ALLOCATE( eforce( ngw, nb_l, 1, 1 ))
CALL dforce_all( 1, ce(:,:,:,1), wempt, fi, eforce(:,:,:,1), gv, vpot(:,:,:,1), &
CALL dforce_all( 1, ce(:,:,:,1), wempt, fi, eforce(:,:,:,1), vpot(:,:,:,1), &
fnl(:,1), eigr, ps)
CALL kohn_sham( 1, ce(:,:,:,1), wempt, eforce(:,:,:,1), kp )
CALL dforce_all( 2, ce(:,:,:,2), wempt, fi, eforce(:,:,:,1), gv, vpot(:,:,:,2), &
CALL dforce_all( 2, ce(:,:,:,2), wempt, fi, eforce(:,:,:,1), vpot(:,:,:,2), &
fnl(:,2), eigr, ps)
CALL kohn_sham( 2, ce(:,:,:,2), wempt, eforce(:,:,:,1), kp )

100
CPV/main.f90
View File

@ -69,7 +69,7 @@
! ... declare modules
USE kinds
USE phase_factors_module, ONLY : strucf
USE phase_factors_module, ONLY : strucf, phfacs
USE restart_file, ONLY : writefile, readfile
USE parameters, ONLY: nacx, nspinx
USE runcp_module, ONLY: runcp, runcp_force_pairing
@ -89,7 +89,7 @@
USE atoms_type_module, ONLY: atoms_type, deallocate_atoms_type
USE print_out_module, ONLY: print_legend, printout, print_time, print_sfac, &
printacc
USE cell_module, ONLY: movecell, press, boxdimensions, updatecell, get_celldm
USE cell_module, ONLY: movecell, press, boxdimensions, updatecell
USE empty_states, ONLY: empty
USE polarization, ONLY: ddipole
USE energies, ONLY: dft_energy_type, debug_energies
@ -146,12 +146,13 @@
USE reciprocal_space_mesh, ONLY: gmeshinfo, newg, gindex_closeup
!
USE reciprocal_vectors, ONLY: &
gcutw, & ! Wave function cut-off ( units of (2PI/alat)^2 => tpiba2 )
gcutp, & ! Potentials and Charge density cut-off ( same units )
gcuts, & ! Smooth mesh Potentials and Charge density cut-off ( same units )
gkcut, & ! Wave function augmented cut-off (take into account all G + k_i , same units)
ngw, & !
ngm, & !
mill_l, & ! G-vectors generators
gcutw, & ! Wave function cut-off ( units of (2PI/alat)^2 => tpiba2 )
gcutp, & ! Potentials and Charge density cut-off ( same units )
gcuts, & ! Smooth mesh Potentials and Charge density cut-off ( same units )
gkcut, & ! Wave function augmented cut-off (take into account all G + k_i , same units)
ngw, & !
ngm, & !
ngs
!
USE recvecs_subroutines, ONLY: recvecs_init
@ -181,7 +182,6 @@
INTEGER :: nnrg
INTEGER :: n1, n2, n3
INTEGER :: n1s, n2s, n3s
INTEGER :: ngm_ , ngw_ , ngs_
REAL(dbl) :: ekinc, ekcell, ekinp, erhoold, maxfion
REAL(dbl) :: derho
@ -241,9 +241,6 @@
!
COMPLEX(dbl), ALLOCATABLE :: sfac(:,:)
! reciprocal lattice and reciprocal vectors
TYPE (recvecs) :: gv ! reciprocal lattice
! cell geometry
TYPE (boxdimensions) :: ht_m, ht_0, ht_p ! cell metrics
@ -253,8 +250,6 @@
REAL(dbl), ALLOCATABLE :: vpot(:,:,:,:)
REAL(dbl) :: vnosee, vnosep
REAL(dbl) :: celldm( 6 )
INTEGER :: ibrav
INTEGER :: lds_wfc
REAL(dbl), EXTERNAL :: cclock
@ -277,6 +272,9 @@
CALL init_dimensions( )
CALL cpsizes( nproc )
CALL cpflush( )
timepre = 0.0_dbl
timernl = 0.0_dbl
timerho = 0.0_dbl
@ -310,14 +308,18 @@
ALLOCATE( ei2( -nr2:nr2, nat ) )
ALLOCATE( ei3( -nr3:nr3, nat ) )
! ... get information
!
CALL get_celldm(ibrav, celldm)
! ==================================================================
! initialize system geometry, cell and positions
! ==================================================================
! ... initialization routines
CALL init_geometry( )
! ==================================================================
! initialize variables and types
! ==================================================================
!
CALL init0s(gv, kp, ps, atoms_m, atoms_0, atoms_p, wfill, &
wempt, ht_m, ht_0, fnl, eigr, ei1, ei2, ei3, nspin)
CALL init0s(kp, ps, atoms_m, atoms_0, atoms_p, wfill, &
wempt, ht_m, ht_0, eigr, ei1, ei2, ei3 )
lds_wfc = wfill%lds
IF( force_pairing ) lds_wfc = 1
@ -334,8 +336,7 @@
ce = 0.0d0
fi = 0.0d0
CALL init1s(gv, kp, ps, atoms_m, atoms_0, atoms_p, cm, c0, wfill, &
ce, wempt, ht_m, ht_0, fnl, eigr, fi )
CALL init1s(kp, ps, cm, c0, wfill, ce, wempt, fnl, eigr, fi )
CALL print_legend( )
@ -344,8 +345,7 @@
dfftp%nr1, dfftp%nr2, dfftp%npl, dfftp%nr1x, dfftp%nr2x, &
dfftp%npl, nspin )
ALLOCATE( sfac( atoms_0%nsp, gv%ng_l ) )
ALLOCATE( sfac( ngm , atoms_0%nsp ) )
ALLOCATE( vpot( dfftp%nr1x, dfftp%nr2x, dfftp%npl, nspin), STAT=ierr)
@ -357,7 +357,7 @@
! ... create a new configuration from scratch
!
ttprint = .true.
CALL from_scratch(gv, kp, ps, rhoe, desc, cm, c0, wfill, eigr, ei1, ei2, ei3, &
CALL from_scratch(kp, ps, rhoe, desc, cm, c0, wfill, eigr, ei1, ei2, ei3, &
sfac, fi, ht_0, atoms_0, fnl, vpot, edft )
CALL printout(nfi, atoms_0, ekinc, ekcell, ttprint, &
@ -369,10 +369,10 @@
! ... (Fortran I/O unit number ndr, file fort.<ndr>)
CALL readfile( nfi, tps, c0, cm, wfill, fi, &
atoms_0, atoms_m, avgs, taui, cdmi, ibrav, celldm, ht_m, ht_0, rhoe, &
desc, vpot, gv, kp)
atoms_0, atoms_m, avgs, taui, cdmi, ht_m, ht_0, rhoe, &
desc, vpot, kp)
CALL from_restart( nfi, avgs, gv, kp, ps, rhoe, desc, cm, c0, wfill, eigr, &
CALL from_restart( nfi, avgs, kp, ps, rhoe, desc, cm, c0, wfill, eigr, &
ei1, ei2, ei3, sfac, fi, ht_m, ht_0, atoms_m, atoms_0, fnl, vpot, edft)
velh = ht_m%hvel
@ -426,21 +426,22 @@
IF( thdyn ) THEN
! ... the simulation cell isn't fixed, recompute the reciprocal lattice
CALL newg(gv, kp, ht_0%m1)
CALL newg(kp, ht_0%m1)
END IF
IF(memchk) CALL memstat(1)
IF( tfor .OR. thdyn ) THEN
! ... ionic positions aren't fixed, recompute structure factors
CALL strucf(sfac, atoms_0, eigr, ei1, ei2, ei3, gv)
CALL phfacs( ei1, ei2, ei3, eigr, mill_l, atoms_0%taus, nr1, nr2, nr3, atoms_0%nat )
CALL strucf( sfac, ei1, ei2, ei3, mill_l, ngm )
END IF
IF(memchk) CALL memstat(2)
IF( thdyn ) THEN
! ... recompute local pseudopotential Fourier expansion
CALL formf(ht_0, gv, kp, ps)
CALL formf(ht_0, kp, ps)
END IF
IF(memchk) CALL memstat(3)
@ -450,29 +451,29 @@
IF( ttdiis .AND. t_diis_simple ) THEN
! ... perform DIIS minimization on electronic states
CALL runsdiis(ttprint, rhoe, desc, atoms_0, gv, kp, &
CALL runsdiis(ttprint, rhoe, desc, atoms_0, kp, &
ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, wfill, thdyn, ht_0, fi, ei, &
fnl, vpot, doions, edft )
ELSE IF (ttdiis .AND. t_diis_rot) THEN
! ... perform DIIS minimization with wavefunctions rotation
IF(nspin.GT.1) CALL errore(' cpmain ',' lsd+diis not allowed ',0)
CALL rundiis(ttprint, rhoe, desc, atoms_0, gv, kp, &
CALL rundiis(ttprint, rhoe, desc, atoms_0, kp, &
ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, wfill, thdyn, ht_0, fi, ei, &
fnl, vpot, doions, edft )
ELSE IF ( tconjgrad ) THEN
! ... on entry c0 should contain the wavefunctions to be optimized
CALL runcg(tortho, ttprint, rhoe, desc, atoms_0, gv, kp, &
CALL runcg(tortho, ttprint, rhoe, desc, atoms_0, kp, &
ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, wfill, thdyn, ht_0, fi, ei, &
fnl, vpot, doions, edft, ekin_maxiter, etot_conv_thr )
! ... on exit c0 and cp both contain the updated wave function
! ... cm are overwritten (used as working space)
ELSE IF ( tsteepdesc ) THEN
CALL runsd(tortho, ttprint, ttforce, rhoe, desc, atoms_0, gv, kp, &
CALL runsd(tortho, ttprint, ttforce, rhoe, desc, atoms_0, kp, &
ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, wfill, thdyn, ht_0, fi, ei, &
fnl, vpot, doions, edft, ekin_maxiter, ekin_conv_thr )
ELSE IF ( tconjgrad_ion%active ) THEN
CALL runcg_ion(nfi, tortho, ttprint, rhoe, desc, atoms_p, atoms_0, &
atoms_m, gv, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, wfill, thdyn, ht_0, fi, ei, &
atoms_m, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, wfill, thdyn, ht_0, fi, ei, &
fnl, vpot, doions, edft, tconvthrs%derho, tconvthrs%force, tconjgrad_ion%nstepix, &
tconvthrs%ekin, tconjgrad_ion%nstepex )
! ... when ions are being relaxed by this subroutine they
@ -486,7 +487,7 @@
! ... compute nonlocal pseudopotential
atoms_0%for = 0.0d0
edft%enl = nlrh_m( c0, wfill, ttforce, atoms_0, fi, gv, kp, fnl, ps%wsg, ps%wnl, eigr)
edft%enl = nlrh_m( c0, wfill, ttforce, atoms_0, fi, kp, fnl, ps%wsg, ps%wnl, eigr)
IF(memchk) CALL memstat(5)
@ -494,7 +495,7 @@
timernl = s5 - s4
! ... compute the new charge density "rhoe"
CALL rhoofr(gv, kp, c0, wfill, fi, rhoe, desc, ht_0)
CALL rhoofr(kp, c0, wfill, fi, rhoe, desc, ht_0)
! CALL printrho(nfi, rhoe, atoms_0, ht_0) ! DEBUG
IF(memchk) CALL memstat(6)
@ -505,7 +506,7 @@
! ... vofrhos compute the new DFT potential "vpot", and energies "edft",
! ... ionc forces "fion" and stress "pail".
CALL vofrhos(ttprint, rhoe, desc, tfor, thdyn, ttforce, atoms_0, &
gv, kp, fnl, vpot, ps, c0, wfill, fi, eigr, ei1, ei2, ei3, sfac, timepre, ht_0, edft)
kp, fnl, vpot, ps, c0, wfill, fi, eigr, ei1, ei2, ei3, sfac, timepre, ht_0, edft)
! .. WRITE( stdout,*) 'DEBUG MAIN', atoms_0%for( 1:3, 1:atoms_0%nat )
! CALL debug_energies( edft ) ! DEBUG
@ -528,12 +529,12 @@
if ( force_pairing ) then
! unpaired electron is assumed of spinup and in highest
! index band; and put equal for paired wf spin up and down
CALL runcp_force_pairing(ttprint, tortho, tsde, cm, c0, cp, wfill, gv, &
CALL runcp_force_pairing(ttprint, tortho, tsde, cm, c0, cp, wfill, &
kp, ps, vpot, eigr, fi, ekincs, timerd, &
timeorto, ht_0, ei, fnl, vnosee )
! ekincs(2) = 0
ELSE
CALL runcp(ttprint, tortho, tsde, cm, c0, cp, wfill, gv, &
CALL runcp(ttprint, tortho, tsde, cm, c0, cp, wfill, &
kp, ps, vpot, eigr, fi, ekincs, timerd, &
timeorto, ht_0, ei, fnl, vnosee )
endif
@ -609,7 +610,7 @@
! ... Here find Empty states eigenfunctions and eigenvalues
IF ( ttempst ) THEN
CALL empty(tortho, atoms_0, gv, c0, wfill, ce, wempt, kp, vpot, eigr, ps )
CALL empty(tortho, atoms_0, c0, wfill, ce, wempt, kp, vpot, eigr, ps )
END IF
! ... dipole
@ -625,12 +626,12 @@
! ... Optical properties
IF( ttoptical ) THEN
CALL opticalp(nfi, ht_0, atoms_0, c0, wfill, fi, ce, wempt, vpot, fnl, eigr, ps, gv, kp)
CALL opticalp(nfi, ht_0, atoms_0, c0, wfill, fi, ce, wempt, vpot, fnl, eigr, ps, kp)
END IF
IF( self_interaction /= 0 ) THEN
IF ( nat_localisation > 0 .AND. ttprint ) THEN
CALL localisation( cp( : , nupdwn(1), 1, 1 ), atoms_0, gv, kp, ht_0, desc)
CALL localisation( cp( : , nupdwn(1), 1, 1 ), atoms_0, kp, ht_0, desc)
END IF
END IF
@ -725,8 +726,8 @@
! ... write the restart file
IF( ttsave .OR. ttexit ) THEN
CALL writefile( nfi, tps, c0, cm, wfill, &
fi, atoms_0, atoms_m, avgs, taui, cdmi, ibrav, celldm, &
ht_m, ht_0, rhoe, desc, vpot, gv, kp)
fi, atoms_0, atoms_m, avgs, taui, cdmi, &
ht_m, ht_0, rhoe, desc, vpot, kp)
END IF
@ -756,14 +757,14 @@
IF(tksout) THEN
IF ( force_pairing ) THEN
CALL ks_states_force_pairing(c0, wfill, ce, wempt, fi, gv, kp, ps, vpot, eigr, fnl)
CALL ks_states_force_pairing(c0, wfill, ce, wempt, fi, kp, ps, vpot, eigr, fnl)
ELSE
CALL ks_states(c0, wfill, ce, wempt, fi, gv, kp, ps, vpot, eigr, fnl)
CALL ks_states(c0, wfill, ce, wempt, fi, kp, ps, vpot, eigr, fnl)
END IF
END IF
IF(tprnsfac) THEN
CALL print_sfac(gv, rhoe, desc, sfac)
CALL print_sfac(rhoe, desc, sfac)
END IF
! ... report statistics
@ -812,7 +813,6 @@
IF( ierr /= 0 ) CALL errore(' cpmain ', ' deallocating fnl ', ierr)
CALL deallocate_pseudo(ps)
CALL deallocate_rvecs(gv)
CALL optical_closeup()
CALL gindex_closeup

24
CPV/nl.f90
View File

@ -43,7 +43,7 @@
! ----------------------------------------------
! BEGIN manual
REAL(dbl) FUNCTION nlrh_m(c0, cdesc, tforce, atoms, occ, gv, kp, fnl, wsg, wnl, eigr)
REAL(dbl) FUNCTION nlrh_m(c0, cdesc, tforce, atoms, occ, kp, fnl, wsg, wnl, eigr)
! this routine computes:
! fnl: Kleinman-Bylander pseudopotential terms (see nlsm1,nlsm1_kp)
@ -73,7 +73,6 @@
! ... declare subroutine arguments
COMPLEX(dbl) :: eigr(:,:) ! exp(i G dot r)
TYPE (recvecs), INTENT(IN) :: gv ! G and k vectors
TYPE (kpoints), INTENT(IN) :: kp ! G and k vectors
COMPLEX(dbl), INTENT(INOUT) :: c0(:,:,:,:) ! wave functions
TYPE (wave_descriptor), INTENT(IN) :: cdesc ! wave functions descriptor
@ -90,7 +89,7 @@
! end of declarations
! ----------------------------------------------
CALL nl_projector_m(c0, cdesc, atoms, gv, kp, fnl, wsg, wnl, eigr)
CALL nl_projector_m(c0, cdesc, atoms, kp, fnl, wsg, wnl, eigr)
nlrh_m = nl_energy_m(fnl, atoms, occ, kp, wsg)
! WRITE( stdout,*) 'DEBUG nlrh ', SUM( fnl(1,1)%r ), SUM( wsg ), SUM( wnl )
IF( tforce ) THEN
@ -98,7 +97,7 @@
ispin_wfc = ispin
IF( force_pairing ) ispin_wfc = 1
CALL nl_ionic_forces_v( ispin, c0( :, :, :, ispin_wfc ), cdesc, atoms, occ(:,:,ispin), &
gv, kp, fnl(:,ispin), wsg, wnl, eigr)
kp, fnl(:,ispin), wsg, wnl, eigr)
END DO
END IF
! .. WRITE( stdout,*) 'DEBUG NLRH:',atoms%for(:,1:atoms%nat)
@ -109,7 +108,7 @@
! ----------------------------------------------
! BEGIN manual
SUBROUTINE nl_projector_m(c0, cdesc, atoms, gv, kp, fnl, wsg, wnl, eigr)
SUBROUTINE nl_projector_m(c0, cdesc, atoms, kp, fnl, wsg, wnl, eigr)
! this routine computes:
! fnl: Kleinman-Bylander pseudopotential terms (see nlsm1,nlsm1_kp)
@ -123,13 +122,13 @@
USE pseudo_projector, ONLY: projector
USE atoms_type_module, ONLY: atoms_type
USE control_flags, ONLY: force_pairing
USE reciprocal_space_mesh, ONLY: gkx_l, gk_l
IMPLICIT NONE
! ... declare subroutine arguments
COMPLEX(dbl) :: eigr(:,:) ! exp(i G dot r)
TYPE(atoms_type), INTENT(INOUT) :: atoms ! ions structure
TYPE (recvecs), INTENT(IN) :: gv ! G and k vectors
TYPE (kpoints), INTENT(IN) :: kp ! G and k vectors
COMPLEX(dbl), INTENT(INOUT) :: c0(:,:,:,:) ! wave functions
TYPE (wave_descriptor), INTENT(IN) :: cdesc ! wave functions
@ -151,7 +150,7 @@
DO ik = 1, cdesc%nkl
! WRITE( stdout,*) 'DEBUG nl_projector ', SUM( eigr ), SUM( c0( :, :, ik, ispin_wfc ) )
CALL nlsm1_s( ispin, wnl(:,:,:,ik), atoms, eigr, c0(:, :, ik, ispin_wfc), cdesc, &
gv%khg_l(:,ik), gv%kgx_l(:,:,ik), fnl(ik, ispin))
gk_l(:,ik), gkx_l(:,:,ik), fnl(ik, ispin))
END DO
END DO
RETURN
@ -208,7 +207,7 @@
! ----------------------------------------------
! BEGIN manual
SUBROUTINE nl_ionic_forces_m(c0, cdesc, atoms, occ, gv, kp, fnl, wsg, wnl, eigr)
SUBROUTINE nl_ionic_forces_m(c0, cdesc, atoms, occ, kp, fnl, wsg, wnl, eigr)
! this routine computes:
! fnl/fnlk: Kleinman-Bylander pseudopotential terms (see nlsm1,nlsm1_kp)
@ -240,7 +239,6 @@
! ... declare subroutine arguments
COMPLEX(dbl) :: eigr(:,:) ! exp(i G dot r)
TYPE (recvecs), INTENT(IN) :: gv ! G and k vectors
TYPE (kpoints), INTENT(IN) :: kp ! G and k vectors
COMPLEX(dbl), INTENT(INOUT) :: c0(:,:,:,:) ! wave functions
TYPE (wave_descriptor), INTENT(IN) :: cdesc ! wave functions desc.
@ -256,7 +254,7 @@
ispin_wfc = ispin
IF( force_pairing ) ispin_wfc = 1
!
CALL nl_ionic_forces_v( ispin, c0( :, :, :, ispin_wfc), cdesc, atoms, occ(:,:,ispin), gv, &
CALL nl_ionic_forces_v( ispin, c0( :, :, :, ispin_wfc), cdesc, atoms, occ(:,:,ispin), &
kp, fnl(:,ispin), wsg, wnl, eigr)
!
END DO
@ -267,7 +265,7 @@
! ----------------------------------------------
! BEGIN manual
SUBROUTINE nl_ionic_forces_v( ispin, c0, cdesc, atoms, occ, gv, kp, fnl, wsg, wnl, eigr)
SUBROUTINE nl_ionic_forces_v( ispin, c0, cdesc, atoms, occ, kp, fnl, wsg, wnl, eigr)
! this routine computes:
! nonlocal potential contribution to forces on ions
@ -289,13 +287,13 @@
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector, allocate_projector, deallocate_projector
USE atoms_type_module, ONLY: atoms_type
USE reciprocal_space_mesh, ONLY: gkx_l, gk_l
IMPLICIT NONE
! ... declare subroutine arguments
COMPLEX(dbl) :: eigr(:,:) ! exp(i G dot r)
TYPE(atoms_type), INTENT(INOUT) :: atoms ! ions structure
TYPE (recvecs), INTENT(IN) :: gv ! G and G+k vectors
TYPE (kpoints), INTENT(IN) :: kp ! K points
COMPLEX(dbl), INTENT(INOUT) :: c0(:,:,:) ! wave functions
TYPE (wave_descriptor), INTENT(IN) :: cdesc ! wave functions desc.
@ -320,7 +318,7 @@
CARTE: DO k = 1, 3 ! x,y,z directions
CALL allocate_projector( dfnl, nsanl, cdesc%nbl( ispin ), ngh, cdesc%gamma )
CALL nlsm2_s( ispin, wnl(:,:,:,ik), atoms, eigr, c0(:,:,ik), cdesc, &
gv%khg_l(:,ik), gv%kgx_l(:,:,ik), dfnl, k)
gk_l(:,ik), gkx_l(:,:,ik), dfnl, k)
CHANN: DO igh = 1, ngh
BANDE: DO ib = 1, cdesc%nbl( ispin )
isa=0

42
CPV/nl_base.f90
View File

@ -25,9 +25,9 @@
! routines in this module:
! REAL(dbl) FUNCTION ene_nl(fnl,wsg,occ,ngh,nspnl,na)
! REAL(dbl) FUNCTION ene_nl_kp(fnlk,wk,wsg,occ,nk,ngh,nspnl,na)
! SUBROUTINE dedh_nls(gv,wnl,wnla,eigr,ll,auxc,gwork,na)
! SUBROUTINE dedh_nlp(gv,gagx,wnl,wnla,eigr,ll,kk,auxc,gwork,na)
! SUBROUTINE dedh_nld(gv,gagx,wnl,wnla,eigr,ll,kk,dm,auxc,gwork,na)
! SUBROUTINE dedh_nls(wnl,wnla,eigr,ll,auxc,gwork,na)
! SUBROUTINE dedh_nlp(gagx,wnl,wnla,eigr,ll,kk,auxc,gwork,na)
! SUBROUTINE dedh_nld(gagx,wnl,wnla,eigr,ll,kk,dm,auxc,gwork,na)
! ----------------------------------------------
! END manual
@ -122,7 +122,7 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE dedh_nls( gv, wnl, wnla, eigr, ll, auxc, gwork, na )
SUBROUTINE dedh_nls( wnl, wnla, eigr, ll, auxc, gwork, na )
! this routine computes the nonlocal contribution of s-orbitals to
! cell derivatives of total energy
@ -131,10 +131,10 @@
! orbitals: s(r) = 1
! ----------------------------------------------
USE cp_types, ONLY: recvecs
USE gvecw, ONLY: ngw
USE reciprocal_vectors, ONLY: gstart
! ... declare subroutine arguments
TYPE (recvecs), INTENT(IN) :: gv
REAL(dbl), INTENT(IN) :: wnl(:,:)
REAL(dbl), INTENT(IN) :: wnla(:,:)
COMPLEX(dbl), INTENT(IN) :: eigr(:,:)
@ -151,14 +151,14 @@
! ----------------------------------------------
IF(ll.GT.0) THEN
ALLOCATE(gwtmp(gv%ngw_l))
ALLOCATE(gwtmp(ngw))
gwtmp(1) = 0.0d0
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
gwtmp(ig) = gwork(ig) * (wnl(ig,ll)-wnla(ig,ll))
END DO
DO ia = 1, na
auxc(1,ia) = CMPLX( 0.0d0 )
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
auxc(ig,ia) = gwtmp(ig) * eigr(ig,ia)
END DO
END DO
@ -173,7 +173,7 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE dedh_nlp( gv, gagx, wnl, wnla, eigr, ll, kk, auxc, gwork, na )
SUBROUTINE dedh_nlp( gagx, wnl, wnla, eigr, ll, kk, auxc, gwork, na )
! this routine computes the nonlocal contribution of p-orbitals to
! cell derivatives of total energy
@ -184,10 +184,10 @@
! p_z(r) = z/r
! ----------------------------------------------
USE cp_types, ONLY: recvecs
USE reciprocal_vectors, ONLY: gstart
USE gvecw, ONLY: ngw
! ... declare subroutine arguments
TYPE (recvecs), INTENT(IN) :: gv
REAL(dbl), INTENT(IN) :: gagx(:,:)
REAL(dbl), INTENT(IN) :: wnl(:,:)
REAL(dbl), INTENT(IN) :: wnla(:,:)
@ -205,15 +205,15 @@
! ----------------------------------------------
IF(ll.GT.0) THEN
ALLOCATE(gwtmp(gv%ngw_l))
ALLOCATE(gwtmp(ngw))
gwtmp(1) = 0.0d0
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
gwtmp(ig) = gagx(kk,ig) * gwork(ig)* &
( 3.d0 * wnl(ig,ll) - wnla(ig,ll) )
END DO
DO ia = 1, na
auxc(1,ia) = CMPLX( 0.0d0 )
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
auxc(ig,ia) = uimag * gwtmp(ig) * eigr(ig,ia)
END DO
END DO
@ -227,7 +227,7 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE dedh_nld( gv, gagx, wnl, wnla, eigr, ll, kk, dm, auxc, gwork, na )
SUBROUTINE dedh_nld( gagx, wnl, wnla, eigr, ll, kk, dm, auxc, gwork, na )
! this routine computes the nonlocal contribution of d-orbitals to
! cell derivatives of total energy
@ -236,10 +236,10 @@
! orbitals: d_m(r) = gkl(x,y,z,r,m)
! ----------------------------------------------
USE cp_types, ONLY: recvecs
USE reciprocal_vectors, ONLY: gstart
USE gvecw, ONLY: ngw
! ... declare subroutine arguments
TYPE (recvecs), INTENT(IN) :: gv
REAL(dbl), INTENT(IN) :: gagx(:,:)
REAL(dbl), INTENT(IN) :: wnl(:,:)
REAL(dbl), INTENT(IN) :: wnla(:,:)
@ -257,16 +257,16 @@
! ----------------------------------------------
IF(ll.GT.0) THEN
ALLOCATE(gwtmp(gv%ngw_l))
ALLOCATE(gwtmp(ngw))
gwtmp(1) = 0.0d0
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
gwtmp(ig) = gwork(ig) * gagx(kk,ig) * &
( 5.d0 * wnl(ig,ll) - wnla(ig,ll) )
gwtmp(ig) = gwtmp(ig) - 2.0d0/3.0d0 * dm * wnl(ig,ll)
END DO
DO ia= 1 , na
auxc(1,ia) = CMPLX( 0.0d0 )
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
auxc(ig,ia) = - gwtmp(ig) * eigr(ig,ia)
END DO
END DO

46
CPV/nlcc.f90
View File

@ -27,12 +27,13 @@
!=----------------------------------------------------------------------------=!
subroutine add_core_charge( rhoetg, rhoetr, sfac, ps, gv, nsp)
subroutine add_core_charge( rhoetg, rhoetr, sfac, ps, nsp)
!=----------------------------------------------------------------------------=!
USE fft, ONLY: pinvfft
USE cp_types, ONLY: pseudo, recvecs
USE cp_types, ONLY: pseudo
use electrons_base, only: nspin
use gvecp, only: ngm
implicit none
@ -40,29 +41,26 @@
COMPLEX(dbl) :: rhoetg(:)
REAL(dbl) :: rhoetr(:,:,:)
type (pseudo), intent(in) :: ps
type (recvecs), intent(in) :: gv
COMPLEX(dbl), INTENT(IN) :: sfac(:,:)
COMPLEX(dbl), ALLOCATABLE :: vtemp(:)
REAL(dbl) :: fac
integer :: is, ig, ng
integer :: is, ig
ng = gv%ng_l
ALLOCATE( vtemp( ng ) )
ALLOCATE( vtemp( ngm ) )
vtemp = CMPLX( 0.0d0 )
fac = 1.0d0 / DBLE( nspin )
DO is = 1, nsp
if( ps%tnlcc( is ) ) then
do ig = 1, ng
vtemp(ig) = vtemp(ig) + fac * sfac( is, ig ) * CMPLX(ps%rhoc1(ig,is),0.0d0)
do ig = 1, ngm
vtemp(ig) = vtemp(ig) + fac * sfac( ig, is ) * CMPLX(ps%rhoc1(ig,is),0.0d0)
end do
endif
end do
rhoetg( 1:ng ) = rhoetg( 1:ng ) + vtemp( 1:ng )
rhoetg( 1:ngm ) = rhoetg( 1:ngm ) + vtemp( 1:ngm )
! rhoetr = 1.0d0 * rhoetr + INVFFT( vtemp )
CALL pinvfft( rhoetr, vtemp, 1.0d0 )
@ -76,7 +74,7 @@
!=----------------------------------------------------------------------------=!
SUBROUTINE core_charge_forces( fion, vxc, rhoc1, tnlcc, atoms, ht, ei1, ei2, ei3, gv, kp )
SUBROUTINE core_charge_forces( fion, vxc, rhoc1, tnlcc, atoms, ht, ei1, ei2, ei3, kp )
!=----------------------------------------------------------------------------=!
! This subroutine computes the non local core correction
@ -87,14 +85,14 @@
USE cell_module, ONLY: boxdimensions
USE brillouin, ONLY: kpoints
USE atoms_type_module, ONLY: atoms_type
USE cp_types, ONLY: recvecs
USE grid_dimensions, ONLY: nr1, nr2, nr3
USE reciprocal_vectors, ONLY: mill_l, gstart, gx
USE gvecp, ONLY: ngm
USE ions_base, ONLY: nat
IMPLICIT NONE
TYPE (atoms_type), INTENT(IN) :: atoms ! atomic positions
TYPE (recvecs), INTENT(IN) :: gv ! reciprocal space vectors
TYPE (boxdimensions), INTENT(IN) :: ht ! cell parameters
TYPE (kpoints), INTENT(IN) :: kp ! k points
COMPLEX(dbl) :: ei1( -nr1:nr1, nat) !
@ -106,7 +104,7 @@
COMPLEX(dbl) :: vxc(:,:) ! XC potential
INTEGER :: ig, ig1, ig2, ig3, isa, ia, is, ispin, nspin
COMPLEX(dbl) :: gx, gy, gz, tx, ty, tz, teigr, cxc
COMPLEX(dbl) :: gxc, gyc, gzc, tx, ty, tz, teigr, cxc
COMPLEX(dbl), ALLOCATABLE :: ftmp(:,:)
REAL(dbl) :: cost
@ -116,13 +114,13 @@
ALLOCATE( ftmp( 3, atoms%nat ) )
ftmp = CMPLX( 0.0d0, 0.0d0 )
DO ig = gv%gstart, gv%ng_l
ig1 = gv%mill(1,ig)
ig2 = gv%mill(2,ig)
ig3 = gv%mill(3,ig)
GX = CMPLX(0.D0, gv%gx_l(1,ig))
GY = CMPLX(0.D0, gv%gx_l(2,ig))
GZ = CMPLX(0.D0, gv%gx_l(3,ig))
DO ig = gstart, ngm
ig1 = mill_l(1,ig)
ig2 = mill_l(2,ig)
ig3 = mill_l(3,ig)
GXC = CMPLX(0.D0, gx(1,ig))
GYC = CMPLX(0.D0, gx(2,ig))
GZC = CMPLX(0.D0, gx(3,ig))
isa = 1
DO is = 1, atoms%nsp
IF ( tnlcc(is) ) THEN
@ -130,9 +128,9 @@
DO ispin = 1, nspin
CXC = CXC + rhoc1( ig, is ) * CONJG( vxc( ig, ispin ) )
END DO
TX = CXC * GX / DBLE( nspin )
TY = CXC * GY / DBLE( nspin )
TZ = CXC * GZ / DBLE( nspin )
TX = CXC * GXC / DBLE( nspin )
TY = CXC * GYC / DBLE( nspin )
TZ = CXC * GZC / DBLE( nspin )
DO IA = 1, atoms%na(is)
teigr = ei1( ig1, isa ) * ei2( ig2, isa ) * ei3( ig3, isa )
ftmp( 1, isa ) = ftmp( 1, isa ) + TEIGR * TX

41
CPV/optical.f90
View File

@ -58,7 +58,7 @@
RETURN
END SUBROUTINE
SUBROUTINE opticalp(nfi, box, atoms, c0, wfill, occ, ce, wempt, vpot, fnl, eigr, ps, gv, kp)
SUBROUTINE opticalp(nfi, box, atoms, c0, wfill, occ, ce, wempt, vpot, fnl, eigr, ps, kp)
USE cp_types
USE cell_module, ONLY: boxdimensions
@ -78,6 +78,8 @@
USE atoms_type_module, ONLY: atoms_type
USE io_files, ONLY: dielecunit, dielecfile
USE control_flags, ONLY: force_pairing
USE reciprocal_vectors, ONLY: gx, g
USE reciprocal_space_mesh, ONLY: gkx_l
INTEGER, INTENT(IN) :: nfi
TYPE(boxdimensions), INTENT(IN) :: box
@ -85,7 +87,6 @@
COMPLEX(dbl), INTENT(IN) :: c0(:,:,:,:)
COMPLEX(dbl), INTENT(INOUT) :: ce(:,:,:,:)
TYPE(wave_descriptor), INTENT(IN) :: wempt, wfill
TYPE(recvecs), INTENT(IN) :: gv
REAL(dbl), INTENT(IN) :: occ(:,:,:)
TYPE(projector) :: fnl(:,:)
REAL (dbl), INTENT(in) :: vpot(:,:,:,:)
@ -133,9 +134,9 @@
nb_l = wfill%nbl( ispin )
ALLOCATE( eforce( ngw, nb_l, nk ) )
CALL nlsm1( ispin, ps%wnl(:,:,:,1), atoms, eigr, cf(:,:,1,ispin), wfill, gv%hg_l, &
gv%gx_l, fnl(1,ispin))
CALL dforce_all( ispin, cf(:,:,:,ispin), wfill, occ(:,:,ispin), eforce, gv, vpot(:,:,:,ispin), &
CALL nlsm1( ispin, ps%wnl(:,:,:,1), atoms, eigr, cf(:,:,1,ispin), wfill, g, &
gx, fnl(1,ispin))
CALL dforce_all( ispin, cf(:,:,:,ispin), wfill, occ(:,:,ispin), eforce, vpot(:,:,:,ispin), &
fnl(:,ispin), eigr, ps)
CALL kohn_sham( ispin, cf(:,:,:,ispin), wfill, eforce, kp )
DEALLOCATE( eforce )
@ -150,9 +151,9 @@
ALLOCATE( eforce( ngw, nb_l, nk ) )
CALL allocate_projector(fnle, nsanl, nb_l, ngh, kp%gamma_only)
CALL nlsm1( ispin, ps%wnl(:,:,:,1), atoms, eigr, ce(:,:,1,ispin), wempt, gv%hg_l, &
gv%gx_l, fnle(1))
CALL dforce_all( ispin, ce(:,:,:,ispin), wempt, ff, eforce, gv, vpot(:,:,:,ispin), &
CALL nlsm1( ispin, ps%wnl(:,:,:,1), atoms, eigr, ce(:,:,1,ispin), wempt, g, &
gx, fnle(1))
CALL dforce_all( ispin, ce(:,:,:,ispin), wempt, ff, eforce, vpot(:,:,:,ispin), &
fnle, eigr, ps)
CALL kohn_sham( ispin, ce(:,:,:,ispin), wempt, eforce, kp )
CALL deallocate_projector(fnle)
@ -225,11 +226,11 @@
IF( gamma_symmetry ) THEN
curr = 0.0d0
DO ig = 1, ngw
curr(1) = curr(1) + gv%gx_l(1,ig) * &
curr(1) = curr(1) + gx(1,ig) * &
AIMAG( ce( ig, cie, ik, ispin ) * CONJG( cf( ig, cif, ik, ispin ) )
curr(2) = curr(2) + gv%gx_l(2,ig) * &
curr(2) = curr(2) + gx(2,ig) * &
AIMAG( ce( ig, cie, ik, ispin ) * CONJG( cf( ig, cif, ik, ispin ) )
curr(3) = curr(3) + gv%gx_l(3,ig) * &
curr(3) = curr(3) + gx(3,ig) * &
AIMAG( ce( ig, cie, ik, ispin ) * CONJG( cf( ig, cif, ik, ispin ) )
END DO
! parallel sum of curr
@ -256,11 +257,11 @@
! IF( gamma_symmetry ) THEN
! curr = 0.0d0
! DO ig = 1, ngw
! curr(1) = curr(1) + gv%gx_l(1,ig) * &
! 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) + gv%gx_l(2,ig) * &
! 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) + gv%gx_l(3,ig) * &
! 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
@ -297,11 +298,11 @@
IF( gamma_symmetry ) THEN
curr = 0.0d0
DO ig = 1, ngw
curr(1) = curr(1) + gv%gx_l(1,ig) * &
curr(1) = curr(1) + gx(1,ig) * &
AIMAG( ce( ig, ie, ik, ispin ) * CONJG( cf( ig, if, ik, ispin ) ) )
curr(2) = curr(2) + gv%gx_l(2,ig) * &
curr(2) = curr(2) + gx(2,ig) * &
AIMAG( ce( ig, ie, ik, ispin ) * CONJG( cf( ig, if, ik, ispin ) ) )
curr(3) = curr(3) + gv%gx_l(3,ig) * &
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, group )
@ -309,11 +310,11 @@
ELSE
ccurr = 0.0d0
DO ig = 1, ngw
ccurr(1) = ccurr(1) + gv%kgx_l(1, ig, ik) * &
ccurr(1) = ccurr(1) + gkx_l(1, ig, ik) * &
ce( ig, ie, ik, ispin ) * CONJG( cf( ig, if, ik, ispin ) )
ccurr(2) = ccurr(2) + gv%kgx_l(2, ig, ik) * &
ccurr(2) = ccurr(2) + gkx_l(2, ig, ik) * &
ce( ig, ie, ik, ispin ) * CONJG( cf( ig, if, ik, ispin ) )
ccurr(3) = ccurr(3) + gv%kgx_l(3, ig, ik) * &
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, group )

352
CPV/phasefactor.f90
View File

@ -12,62 +12,62 @@
! BEGIN manual
!=----------------------------------------------------------------------------=!
MODULE phase_factors_module
MODULE phase_factors_module
!=----------------------------------------------------------------------------=!
! (describe briefly what this module does...)
! ----------------------------------------------
! routines in this module:
! SUBROUTINE phfacs(eigr,atoms)
! SUBROUTINE strucf(atoms,eigr,gv)
! SUBROUTINE phfacs
! SUBROUTINE strucf
! ----------------------------------------------
! END manual
! ... declare modules
USE kinds
IMPLICIT NONE
SAVE
PRIVATE
PUBLIC :: strucf
PUBLIC :: strucf, phfacs
! ... --+ end of module-scope declarations +--
!=----------------------------------------------------------------------------=!
CONTAINS
CONTAINS
!=----------------------------------------------------------------------------=!
!=----------------------------------------------------------------------------=!
! BEGIN manual
SUBROUTINE phfacs( ei1, ei2, ei3, eigr, mill, taus, nat )
SUBROUTINE phfacs( ei1, ei2, ei3, eigr, mill, taus, nr1, nr2, nr3, nat )
! this routine computes the phase factors
!
! ei1(ix,ia) = exp (i ix G_1 dot R(ia))
! ei2(iy,ia) = exp (i iy G_2 dot R(ia))
! ei3(iz,ia) = exp (i iz G_3 dot R(ia))
!
! eigr(ig,ia) = exp(i G dot R(ia)) =
! = ei1(ix,ia) * ei2(iy,ia) * ei3(iz,ia)
!
! G_1,G_2,G_3 = reciprocal lattice generators
!
! ia = index of ion
! ig = index of G vector
! ix,iy,iz = Miller indices
! ----------------------------------------------
! this routine computes the phase factors
!
! ei1(ix,ia) = exp (i ix G_1 dot R(ia))
! ei2(iy,ia) = exp (i iy G_2 dot R(ia))
! ei3(iz,ia) = exp (i iz G_3 dot R(ia))
!
! eigr(ig,ia) = exp(i G dot R(ia)) =
! = ei1(ix,ia) * ei2(iy,ia) * ei3(iz,ia)
!
! G_1,G_2,G_3 = reciprocal lattice generators
!
! ia = index of ion
! ig = index of G vector
! ix,iy,iz = Miller indices
! ----------------------------------------------
! END manual
USE kinds, ONLY: dbl
USE constants, ONLY: tpi
USE grid_dimensions, ONLY: nr1, nr2, nr3
! ... declare subroutine arguments
! ... declare subroutine arguments
INTEGER, INTENT(IN) :: nat
INTEGER, INTENT(IN) :: nr1, nr2, nr3
COMPLEX(dbl) :: ei1( -nr1 : nr1, nat )
COMPLEX(dbl) :: ei2( -nr2 : nr2, nat )
COMPLEX(dbl) :: ei3( -nr3 : nr3, nat )
@ -75,48 +75,53 @@
REAL(dbl) :: taus( 3, nat )
INTEGER :: mill( :, : )
! ... declare other variables
! ... declare other variables
COMPLEX(dbl) :: ctep1, ctep2, ctep3, ctem1, ctem2, ctem3
REAL(dbl) :: ar1, ar2, ar3
INTEGER :: i, j, k, isa
INTEGER :: ig, ig1, ig2, ig3, ngw
! ... --+ end of declarations +--
! ... --+ end of declarations +--
if(nr1.lt.3) call errore(' phfacs ',' nr1 too small ',nr1)
if(nr2.lt.3) call errore(' phfacs ',' nr2 too small ',nr2)
if(nr3.lt.3) call errore(' phfacs ',' nr3 too small ',nr3)
DO isa = 1, nat
! ... Miller index = 0: exp(i 0 dot R(ia)) = 1
! ... Miller index = 0: exp(i 0 dot R(ia)) = 1
ei1( 0, isa ) = CMPLX( 1.d0, 0.d0 )
ei2( 0, isa ) = CMPLX( 1.d0, 0.d0 )
ei3( 0, isa ) = CMPLX( 1.d0, 0.d0 )
! ... let R_1,R_2,R_3 be the direct lattice generators,
! ... G_1,G_2,G_3 the reciprocal lattice generators
! ... by definition G_i dot R_j = 2 pi delta_{ij}
! ... ionic coordinates are in units of R_1,R_2,R_3
! ... then G_i dot R(ia) = 2 pi R(ia)_i
! ... let R_1,R_2,R_3 be the direct lattice generators,
! ... G_1,G_2,G_3 the reciprocal lattice generators
! ... by definition G_i dot R_j = 2 pi delta_{ij}
! ... ionic coordinates are in units of R_1,R_2,R_3
! ... then G_i dot R(ia) = 2 pi R(ia)_i
ar1 = tpi * taus(1,isa) ! G_1 dot R(ia)
ar2 = tpi * taus(2,isa) ! G_2 dot R(ia)
ar3 = tpi * taus(3,isa) ! G_3 dot R(ia)
! ... Miller index = 1: exp(i G_i dot R(ia))
! ... Miller index = 1: exp(i G_i dot R(ia))
ctep1 = CMPLX( cos( ar1 ), -sin( ar1 ) )
ctep2 = CMPLX( cos( ar2 ), -sin( ar2 ) )
ctep3 = CMPLX( cos( ar3 ), -sin( ar3 ) )
! ... Miller index = -1: exp(-i G_i dot R(ia))
! ... Miller index = -1: exp(-i G_i dot R(ia))
ctem1 = CONJG(ctep1)
ctem2 = CONJG(ctep2)
ctem3 = CONJG(ctep3)
! ... Miller index > 0: exp(i N G_i dot R(ia)) =
! ... = exp(i G_i dot R(ia)) exp(i (N-1) G_i dot R(ia))
! ... Miller index < 0: exp(-i N G_i dot R(ia)) =
! ... = exp(-i G_i dot R(ia)) exp(-i (N-1) G_i dot R(ia))
! ... Miller index > 0: exp(i N G_i dot R(ia)) =
! ... = exp(i G_i dot R(ia)) exp(i (N-1) G_i dot R(ia))
! ... Miller index < 0: exp(-i N G_i dot R(ia)) =
! ... = exp(-i G_i dot R(ia)) exp(-i (N-1) G_i dot R(ia))
DO i = 1, nr1 - 1
ei1( i, isa ) = ei1( i - 1, isa ) * ctep1
@ -148,13 +153,13 @@
END DO
RETURN
END SUBROUTINE phfacs
END SUBROUTINE phfacs
!=----------------------------------------------------------------------------=!
! BEGIN manual
SUBROUTINE strucf( sfac, atoms, eigr, ei1, ei2, ei3, gv )
SUBROUTINE strucf( sfac, ei1, ei2, ei3, mill, ngm )
! this routine computes the structure factors
!
@ -174,49 +179,49 @@
! ----------------------------------------------
! END manual
USE atoms_type_module, ONLY: atoms_type
USE cp_types, ONLY: recvecs
USE grid_dimensions, ONLY: nr1, nr2, nr3
USE ions_base, ONLY: nat
USE kinds, ONLY: dbl
USE ions_base, ONLY: nat, na, nsp
use grid_dimensions, only: nr1, nr2, nr3
IMPLICIT NONE
! ... declare subroutine arguments
TYPE (atoms_type) :: atoms
COMPLEX(dbl) :: eigr(:,:)
! ... declare subroutine arguments
!
COMPLEX(dbl) :: ei1( -nr1 : nr1, nat )
COMPLEX(dbl) :: ei2( -nr2 : nr2, nat )
COMPLEX(dbl) :: ei3( -nr3 : nr3, nat )
TYPE (recvecs) :: gv
INTEGER :: mill( :, : )
INTEGER :: ngm
COMPLEX(dbl), INTENT(OUT) :: sfac(:,:)
! ... declare other variables
! ... declare other variables
!
INTEGER :: is, ig, ia, ig1, ig2, ig3, isa
! ... --+ end of declarations +--
CALL phfacs( ei1, ei2, ei3, eigr, gv%mill, atoms%taus, atoms%nat )
call start_clock( 'strucf' )
DO ig = 1, gv%ng_l
ig1 = gv%mill( 1, ig )
ig2 = gv%mill( 2, ig )
ig3 = gv%mill( 3, ig )
DO ig = 1, ngm
ig1 = mill( 1, ig )
ig2 = mill( 2, ig )
ig3 = mill( 3, ig )
isa = 1
DO is = 1, atoms%nsp
sfac( is, ig ) = CMPLX (0.0_dbl, 0.0_dbl)
DO ia = 1, atoms%na(is)
sfac( is, ig ) = sfac( is, ig ) + &
DO is = 1, nsp
sfac( ig, is ) = CMPLX (0.0d0, 0.0d0)
DO ia = 1, na(is)
sfac( ig, is ) = sfac( ig, is ) + &
ei1( ig1, isa ) * ei2( ig2, isa ) * ei3( ig3, isa )
isa = isa + 1
END DO
END DO
END DO
call stop_clock( 'strucf' )
RETURN
END SUBROUTINE strucf
!=----------------------------------------------------------------------------=!
END MODULE phase_factors_module
END MODULE phase_factors_module
!=----------------------------------------------------------------------------=!
@ -231,164 +236,39 @@
! The value of n1,n2,n3 for a vector g is supplied by arrays mill_l
! calculated in ggen .
!
use constants, only: pi, fpi
use parameters, only: natx, nsx
use control_flags, only: iprint, iprsta
use control_flags, only: iprsta
use io_global, only: stdout
use ions_base, only: nsp, na, nat
use cell_base, only: ainv
use cell_base, only: ainv, r_to_s
use grid_dimensions, only: nr1, nr2, nr3
use gvecw, only: ngw
use reciprocal_vectors, only: mill_l
use gvecw, only: ngw
use phase_factors_module, only: phfacs
!
implicit none
real(kind=8) tau0(3,natx)
real(kind=8) tau0(3,nat)
!
complex(kind=8) ei1(-nr1:nr1,nat), ei2(-nr2:nr2,nat), &
& ei3(-nr3:nr3,nat), eigr(ngw,nat)
!
integer i,j,k, ia, is, ig, isa
real(kind=8) taup(3), s, ar1,ar2,ar3
complex(kind=8) ctep1,ctep2,ctep3,ctem1,ctem2,ctem3
integer :: i, isa
real(kind=8), allocatable :: taus(:,:)
!
if(nr1.lt.3) call errore(' phfac ',' nr1 too small ',nr1)
if(nr2.lt.3) call errore(' phfac ',' nr1 too small ',nr2)
if(nr3.lt.3) call errore(' phfac ',' nr1 too small ',nr3)
allocate( taus(3,nat) )
!
if(iprsta.gt.3) then
WRITE( stdout,*) ' phfac: tau0 '
WRITE( stdout,*) ( ( tau0(i,isa), i=1, 3 ), isa=1, SUM(na(1:nsp)) )
endif
isa = 0
do is=1,nsp
do ia=1,na(is)
isa = isa + 1
do i=1,3
s=0.d0
do j=1,3
s=s+ainv(i,j)*tau0(j,isa)
end do
taup(i)=s
!
! tau0=x1*a1+x2*a2+x3*a3 => taup(1)=x1=tau0*b1 and so on
!
end do
!
ar1=2.d0*pi*taup(1)
ctep1=cmplx(cos(ar1),-sin(ar1))
ctem1=conjg(ctep1)
ei1( 0,isa)=cmplx(1.d0,0.d0)
ei1( 1,isa)=ctep1
ei1(-1,isa)=ctem1
do i=2,nr1-1
ei1( i,isa)=ei1( i-1,isa)*ctep1
ei1(-i,isa)=ei1(-i+1,isa)*ctem1
end do
!
ar2=2.d0*pi*taup(2)
ctep2=cmplx(cos(ar2),-sin(ar2))
ctem2=conjg(ctep2)
ei2( 0,isa)=cmplx(1.d0,0.d0)
ei2( 1,isa)=ctep2
ei2(-1,isa)=ctem2
do j=2,nr2-1
ei2( j,isa)=ei2( j-1,isa)*ctep2
ei2(-j,isa)=ei2(-j+1,isa)*ctem2
end do
!
ar3=2.d0*pi*taup(3)
ctep3=cmplx(cos(ar3),-sin(ar3))
ctem3=conjg(ctep3)
ei3( 0,isa)=cmplx(1.d0,0.d0)
ei3( 1,isa)=ctep3
ei3(-1,isa)=ctem3
do k=2,nr3-1
ei3( k,isa)=ei3( k-1,isa)*ctep3
ei3(-k,isa)=ei3(-k+1,isa)*ctem3
end do
!
end do
end do
!
if(iprsta.gt.4) then
WRITE( stdout,*)
if(nsp.gt.1) then
isa = 0
do is=1,nsp
WRITE( stdout,'(33x,a,i4)') ' ei1, ei2, ei3 (is)',is
do ig=1,4
WRITE( stdout,'(6f9.4)') &
& ei1(ig,1+isa),ei2(ig,1+isa),ei3(ig,1+isa)
end do
WRITE( stdout,*)
isa = isa + na(is)
end do
else
do ia=1,na(1)
WRITE( stdout,'(33x,a,i4)') ' ei1, ei2, ei3 (ia)',ia
do ig=1,4
WRITE( stdout,'(6f9.4)') &
& ei1(ig,ia),ei2(ig,ia),ei3(ig,ia)
end do
WRITE( stdout,*)
end do
endif
endif
!
! calculation of eigr(g,ia,is)=e^(-ig.r(ia,is))
!
do isa=1,nat
do ig=1,ngw
i = mill_l(1,ig)
j = mill_l(2,ig)
k = mill_l(3,ig)
eigr(ig,isa) = ei1(i,isa) * ei2(j,isa) * ei3(k,isa)
end do
end do
CALL r_to_s( tau0, taus, na, nsp, ainv )
CALL phfacs( ei1, ei2, ei3, eigr, mill_l, taus, nr1, nr2, nr3, nat )
deallocate( taus )
!
return
end
!-------------------------------------------------------------------------
subroutine strucf (ei1,ei2,ei3,sfac)
!-----------------------------------------------------------------------
! computes the structure factor sfac(ngs,nsp) and returns it in "pseu"
!
!
use gvecs, only: ngs
use constants, only: pi, fpi
use grid_dimensions, only: nr1, nr2, nr3
use ions_base, only: nsp, na, nat
use reciprocal_vectors, only: mill_l
!
implicit none
complex(kind=8) ei1(-nr1:nr1,nat), ei2(-nr2:nr2,nat), &
& ei3(-nr3:nr3,nat), sfac(ngs,nsp)
integer is, ig, ia, i, j, k, isa
!
call start_clock( 'strucf' )
isa = 1
do is=1,nsp
do ig=1,ngs
sfac(ig,is)=(0.,0.)
end do
do ia=1,na(is)
do ig=1,ngs
i = mill_l( 1, ig )
j = mill_l( 2, ig )
k = mill_l( 3, ig )
sfac(ig,is)=sfac(ig,is) + ei1(i,isa) *ei2(j,isa) *ei3(k,isa)
end do
isa = isa + 1
end do
end do
!
call stop_clock( 'strucf' )
return
end subroutine
!-----------------------------------------------------------------------
subroutine phbox(taub,eigrb)
!-----------------------------------------------------------------------
@ -397,23 +277,21 @@
! Uses the same logic for fast calculation as in phfac (see below)
!
use io_global, only: stdout
use control_flags, only: iprint, iprsta
use constants, only: pi, fpi
use parameters, only: natx, nsx
use control_flags, only: iprsta
use ions_base, only: nsp, na, nat
use gvecb, only: ngb, mill_b
use cell_base, only: ainv
use cell_base, only: r_to_s
use small_box, only: ainvb
use smallbox_grid_dimensions, only: nr1b, nr2b, nr3b
use phase_factors_module, only: phfacs
!
implicit none
real(kind=8) :: taub(3,natx)
real(kind=8) :: taub(3,nat)
complex(kind=8) :: eigrb(ngb,nat)
! local
integer :: i,j,k, is, ia, ig, isa
real(kind=8) taup(3),s,ar1,ar2,ar3
complex(kind=8), allocatable:: ei1b(:,:), ei2b(:,:), ei3b(:,:)
complex(kind=8) ctep1,ctep2,ctep3,ctem1,ctem2,ctem3
real(kind=8), allocatable :: taus(:,:)
!
if(nr1b.lt.3) call errore(' phbox ',' nr1b too small ',nr1b)
if(nr2b.lt.3) call errore(' phbox ',' nr2b too small ',nr2b)
@ -422,67 +300,14 @@
allocate(ei1b(-nr1b:nr1b,nat))
allocate(ei2b(-nr2b:nr2b,nat))
allocate(ei3b(-nr3b:nr3b,nat))
allocate( taus( 3, nat ) )
!
if(iprsta.gt.3) then
WRITE( stdout,*) ' phbox: taub '
WRITE( stdout,*) ( (taub(i,isa), i=1, 3 ), isa=1, SUM(na(1:nsp)) )
endif
do isa=1,nat
do i=1,3
s=0.d0
do j=1,3
s = s + ainvb(i,j)*taub(j,isa)
end do
taup(i)=s
end do
!
ar1=2.d0*pi*taup(1)
ctep1=cmplx(cos(ar1),-sin(ar1))
ctem1=conjg(ctep1)
ei1b( 0,isa)=cmplx(1.d0,0.d0)
ei1b( 1,isa)=ctep1
ei1b(-1,isa)=ctem1
do i=2,nr1b-1
ei1b( i,isa)=ei1b( i-1,isa)*ctep1
ei1b(-i,isa)=ei1b(-i+1,isa)*ctem1
end do
!
ar2=2.d0*pi*taup(2)
ctep2=cmplx(cos(ar2),-sin(ar2))
ctem2=conjg(ctep2)
ei2b( 0,isa)=cmplx(1.d0,0.d0)
ei2b( 1,isa)=ctep2
ei2b(-1,isa)=ctem2
do j=2,nr2b-1
ei2b( j,isa)=ei2b( j-1,isa)*ctep2
ei2b(-j,isa)=ei2b(-j+1,isa)*ctem2
end do
!
ar3=2.d0*pi*taup(3)
ctep3=cmplx(cos(ar3),-sin(ar3))
ctem3=conjg(ctep3)
ei3b( 0,isa)=cmplx(1.d0,0.d0)
ei3b( 1,isa)=ctep3
ei3b(-1,isa)=ctem3
do k=2,nr3b-1
ei3b( k,isa)=ei3b( k-1,isa)*ctep3
ei3b(-k,isa)=ei3b(-k+1,isa)*ctem3
end do
!
end do
!
! calculation of eigrb(g,ia,is)=e^(-ig.r(ia,is))
!
do isa=1,nat
do ig=1,ngb
i = mill_b(1,ig)
j = mill_b(2,ig)
k = mill_b(3,ig)
eigrb(ig,isa) = ei1b(i,isa) * ei2b(j,isa) * ei3b(k,isa)
end do
end do
CALL r_to_s( taub, taus, na, nsp, ainvb )
CALL phfacs( ei1b, ei2b, ei3b, eigrb, mill_b, taus, nr1b, nr2b, nr3b, nat )
!
if(iprsta.gt.4) then
WRITE( stdout,*)
@ -512,6 +337,7 @@
deallocate(ei3b)
deallocate(ei2b)
deallocate(ei1b)
deallocate( taus )
!
return
end subroutine

189
CPV/potentials.f90
View File

@ -97,26 +97,26 @@
!! ...................................................................... !!
!! ...................................................................... !!
SUBROUTINE vofmean( gv, sfac, rhops, rhoeg )
SUBROUTINE vofmean( sfac, rhops, rhoeg )
USE constants, ONLY: fpi
USE cell_base, ONLY: tpiba2, tpiba
USE mp, ONLY: mp_sum
USE mp_global, ONLY: nproc, mpime, group, root
USE io_global, ONLY: ionode
USE cp_types, ONLY: recvecs
USE reciprocal_vectors, ONLY: gstart, gx, g
USE gvecp, ONLY: ngm
TYPE (recvecs), INTENT(IN) :: gv
REAL(dbl), INTENT(IN) :: RHOPS(:,:)
COMPLEX(dbl), INTENT(IN) :: RHOEG(:)
COMPLEX(dbl), INTENT(IN) :: sfac(:,:)
COMPLEX(dbl) fpi_tpiba2, rp, vcg
COMPLEX(dbl) :: fpi_tpiba2, rp, vcg
REAL(dbl) :: gxt, dr, r
REAL(dbl), ALLOCATABLE :: vrmean(:)
REAL(dbl) dr, gx, r
INTEGER gstart, ig, is, iasse, ipiano1, ipiano2
INTEGER ir
INTEGER :: ig, is, iasse, ipiano1, ipiano2
INTEGER :: ir
IF( (vhasse.NE.'X') .AND. (vhasse.NE.'Y') .AND. (vhasse.NE.'Z') ) THEN
CALL errore( ' vofmean ', ' wrong asse ',0)
@ -128,8 +128,6 @@
CALL errore( ' vofmean ', ' wrong nr ',0)
END IF
gstart = 1
IF(gv%gzero) gstart = 2
fpi_tpiba2 = CMPLX(FPI/TPIBA2,0.0d0)
dr = ( vhrmax - vhrmin ) / vhnr
@ -152,17 +150,17 @@
DO ir = 1, vhnr
vrmean(ir) = 0.0d0
r = vhrmin + (ir-1)*dr
DO ig = gstart, gv%ng_l
DO ig = gstart, ngm
rp = (0.D0,0.D0)
DO is = 1, SIZE( sfac, 1 )
rp = rp + sfac( is, ig ) * rhops( ig, is )
DO is = 1, SIZE( sfac, 2 )
rp = rp + sfac( ig, is ) * rhops( ig, is )
END DO
IF((gv%gx_l(ipiano1,IG).EQ.0.d0).AND. &
(gv%gx_l(ipiano2,IG).EQ.0.d0))THEN
vcg = fpi_tpiba2 * (rhoeg(ig) + rp) / gv%hg_l(ig)
gx = gv%gx_l(iasse, ig) * tpiba
vrmean(ir) = vrmean(ir) + REAL(vcg) * COS(gx*r)
vrmean(ir) = vrmean(ir) - AIMAG(vcg) * SIN(gx*r)
IF((gx(ipiano1,IG).EQ.0.d0).AND. &
(gx(ipiano2,IG).EQ.0.d0))THEN
vcg = fpi_tpiba2 * (rhoeg(ig) + rp) / g(ig)
gxt = gx(iasse, ig) * tpiba
vrmean(ir) = vrmean(ir) + REAL(vcg) * COS(gxt*r)
vrmean(ir) = vrmean(ir) - AIMAG(vcg) * SIN(gxt*r)
END IF
END DO
vrmean(ir) = 2.0d0 * vrmean(ir)
@ -184,14 +182,14 @@
! -------------------------------------------------------------------------
SUBROUTINE kspotential( tprint, tforce, tstress, rhoe, desc, &
atoms, gv, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht, fi, fnl, vpot, edft, timepre )
atoms, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht, fi, fnl, vpot, edft, timepre )
USE charge_density, ONLY: rhoofr
USE nl, ONLY: nlrh_m
USE energies, ONLY: dft_energy_type
USE cell_module, ONLY: boxdimensions
USE brillouin, ONLY: kpoints
USE cp_types, ONLY: pseudo, recvecs
USE cp_types, ONLY: pseudo
USE atoms_type_module, ONLY: atoms_type
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
@ -208,7 +206,6 @@
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
COMPLEX(dbl) :: eigr(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
TYPE (boxdimensions), INTENT(INOUT) :: ht
REAL(dbl) :: fi(:,:,:)
@ -220,11 +217,11 @@
REAL(dbl), INTENT(OUT) :: timepre
TYPE (charge_descriptor), INTENT(IN) :: desc
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fi, gv, kp, fnl, ps%wsg, ps%wnl, eigr)
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fi, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL rhoofr(gv, kp, c0, cdesc, fi, rhoe, desc, ht)
CALL rhoofr(kp, c0, cdesc, fi, rhoe, desc, ht)
CALL vofrhos(tprint, rhoe, desc, tforce, tstress, tforce, atoms, gv, &
CALL vofrhos(tprint, rhoe, desc, tforce, tstress, tforce, atoms, &
kp, fnl, vpot, ps, c0, cdesc, fi, eigr, ei1, ei2, ei3, sfac, timepre, ht, edft)
RETURN
@ -236,7 +233,7 @@
! BEGIN manual
SUBROUTINE vofrhos(tprint, rhoe, desc, tfor, thdyn, tforce, atoms, &
gv, kp, fnl, vpot, ps, c0, cdesc, fi, eigr, ei1, ei2, ei3, sfac, timepre, box, edft)
kp, fnl, vpot, ps, c0, cdesc, fi, eigr, ei1, ei2, ei3, sfac, timepre, box, edft)
! this routine computes:
! ekin = dft_kinetic term of the DFT functional (see dft_kinetic_energy)
@ -274,9 +271,11 @@
USE charge_types, ONLY: charge_descriptor
USE control_flags, ONLY: tscreen, tchi2, iprsta
USE io_global, ONLY: ionode
USE cp_types, ONLY: recvecs, pseudo
USE cp_types, ONLY: pseudo
USE io_global, ONLY: stdout
USE sic_module, ONLY: self_interaction, si_epsilon
USE reciprocal_vectors, ONLY: gx
USE gvecp, ONLY: ngm
IMPLICIT NONE
@ -292,7 +291,6 @@
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
COMPLEX(dbl) :: eigr(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
COMPLEX(dbl), INTENT(IN) :: c0(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
@ -383,7 +381,7 @@
IF(tchi2) THEN
CALL allocate_chi2(gv%ng_l)
CALL allocate_chi2(ngm)
END IF
ALLOCATE( rhoetr( nr1x, nr2x, nr3x, nspin) )
@ -400,8 +398,8 @@
ALLOCATE( v2xc( 1, 1, 1, 1, 1 ) )
END IF
ALLOCATE( rhoeg(gv%ng_l, nspin) )
ALLOCATE( rhoetg(gv%ng_l, nspin) )
ALLOCATE( rhoeg(ngm, nspin) )
ALLOCATE( rhoetg(ngm, nspin) )
edft%self_sxc = 0.d0
edft%sxc = 0.d0
@ -439,7 +437,7 @@
edft%ekin = 0.0_dbl
edft%emkin = 0.0_dbl
edft%ekin = dft_kinetic_energy(c0, cdesc, gv, kp, tecfix, fi, rsum, edft%emkin)
edft%ekin = dft_kinetic_energy(c0, cdesc, kp, tecfix, fi, rsum, edft%emkin)
IF(tprint) THEN
CALL checkrho(rhoe, desc, rsum, omega)
@ -477,7 +475,7 @@
! ... add core correction
! ... rhoetg = rhoeg + cc
! ... rhoetr = rhoe + cc
CALL add_core_charge( rhoetg(:,ispin), rhoetr(:,:,:,ispin), sfac, ps, gv, atoms%nsp )
CALL add_core_charge( rhoetg(:,ispin), rhoetr(:,:,:,ispin), sfac, ps, atoms%nsp )
ELSE
! ... no core correction
@ -493,7 +491,7 @@
END IF
IF(tgc) THEN
CALL gradrho( rhoetg(:,ispin), grho(:,:,:,:,ispin), gv%gx_l )
CALL gradrho( rhoetg(:,ispin), grho(:,:,:,:,ispin), gx )
END IF
END DO
@ -505,7 +503,7 @@
! ... Self-interaction correction --- Excor Part
!In any case calculate the Excor part with rhoetr
CALL exch_corr_energy(rhoetr, rhoetg, grho, vpot, sxcp, vxc, v2xc, gv)
CALL exch_corr_energy(rhoetr, rhoetg, grho, vpot, sxcp, vxc, v2xc)
IF( ttsic ) THEN
IF( ionode ) THEN
@ -537,7 +535,7 @@
self_grho(:,:,:,:,2) = 0.D0
ENDIF
CALL exch_corr_energy(self_rho, rhoetg, self_grho, self_vpot, &
self_sxcp, self_vxc, self_v2xc, gv)
self_sxcp, self_vxc, self_v2xc)
vpot(:,:,:,1) = vpot(:,:,:,1) - self_vpot(:,:,:,1)
vpot(:,:,:,2) = vpot(:,:,:,2) + self_vpot(:,:,:,1)
IF (tgc) THEN
@ -570,7 +568,7 @@
! write(stdout,*)'DA XC SELF_RHO2',self_rho(:,:,:,2)
CALL exch_corr_energy(self_rho, rhoetg, self_grho, self_vpot, &
self_sxcp, self_vxc, self_v2xc, gv)
self_sxcp, self_vxc, self_v2xc)
vpot (:,:,:,2) = self_vpot(:,:,:,2) + self_vpot(:,:,:,1)
vpot (:,:,:,1) = 0.d0
@ -613,7 +611,7 @@
END DO
! ... now rhoetg contains the xc potential
IF (ttforce) THEN
CALL core_charge_forces(fion, rhoetg, ps%rhoc1, ps%tnlcc, atoms, box, ei1, ei2, ei3, gv, kp )
CALL core_charge_forces(fion, rhoetg, ps%rhoc1, ps%tnlcc, atoms, box, ei1, ei2, ei3, kp )
END IF
END IF
@ -629,9 +627,9 @@
! ... local part of the pseudopotential and its energy contribution (eps)
! ... Self-interaction correction --- Hartree part
ALLOCATE( vloc( gv%ng_l ) )
ALLOCATE( vloc( ngm ) )
CALL vofloc(ttscreen, ttforce, edft%ehte, edft%ehti, ehp, &
eps, vloc, rhoeg, fion, atoms, ps%rhops, ps%vps, gv, kp, eigr, &
eps, vloc, rhoeg, fion, atoms, ps%rhops, ps%vps, kp, eigr, &
ei1, ei2, ei3, sfac, box, desc, ps%ap )
! edft%ehte = REAL ( ehtep )
@ -642,7 +640,7 @@
! Delta_Esic to Hartree is ALWAYS = -EH(rhoup-rhodw)
ALLOCATE(self_vloc(gv%ng_l), self_rhoeg(gv%ng_l), STAT = ierr)
ALLOCATE(self_vloc(ngm), self_rhoeg(ngm), STAT = ierr)
IF( ierr /= 0 )&
CALL errore(' vofrhos ', ' allocating self_vloc, self_rhoeg ', ierr)
@ -652,7 +650,7 @@
! working on the total charge density
CALL self_vofloc(ttscreen, self_ehtep, self_vloc, self_rhoeg, &
gv, kp, box, desc)
kp, box, desc)
CALL pinvfft(self_vpot(:,:,:,1), self_vloc(:))
self_vpot(:,:,:,1) = si_epsilon * self_vpot(:,:,:,1)
@ -738,7 +736,7 @@
! ... compute stress tensor
IF( ttstress .AND. kp%gamma_only ) THEN
s8 = cclock()
CALL pstress( strvxc, rhoeg, rhoetg, pail, desr, gv, fnl, &
CALL pstress( strvxc, rhoeg, rhoetg, pail, desr, fnl, &
ps, c0, cdesc, fi, eigr, sfac, grho, v2xc, box, edft)
timepre = cclock() - s8
END IF
@ -753,11 +751,11 @@
IF(tprint.AND.tvhmean) THEN
IF(nspin.GT.2) CALL errore(' vofrho ',' spin + vmean ',nspin)
IF(nspin==1) CALL vofmean( gv, sfac, ps%rhops, rhoeg(:,1) )
IF(nspin==1) CALL vofmean( sfac, ps%rhops, rhoeg(:,1) )
IF(nspin==2) THEN
ALLOCATE(rho12(gv%ng_l))
ALLOCATE(rho12(ngm))
rho12 (:) = rhoeg(:,1)+rhoeg(:,2)
CALL vofmean( gv, sfac, ps%rhops, rho12)
CALL vofmean( sfac, ps%rhops, rho12)
DEALLOCATE(rho12)
END IF
END IF
@ -921,19 +919,19 @@
! BEGIN manual
SUBROUTINE vofloc(tscreen, ttforce, ehte, ehti, eh, eps, vloc, rhoeg, &
fion, atoms, rhops, vps, gv, kp, eigr, ei1, ei2, ei3, sfac, ht, desc, ap)
fion, atoms, rhops, vps, kp, eigr, ei1, ei2, ei3, sfac, ht, desc, ap)
! this routine computes:
! omega = ht%deth
! rho_e(ig) = (sum over ispin) rhoeg(ig,ispin)
! rho_I(ig) = (sum over is) sfac(is,ig) * rhops(ig,is)
! vloc_h(ig) = fpi / ( gv%hg_l(ig) * tpiba2 ) * { rho_e(ig) + rho_I(ig) }
! vloc_h(ig) = fpi / ( g(ig) * tpiba2 ) * { rho_e(ig) + rho_I(ig) }
! vloc_ps(ig) = (sum over is) sfac(is,ig) * vps(ig,is)
!
! Eps = Fact * omega * (sum over ig) CMPLX( rho_e(ig) ) * vloc_ps(ig)
! if Gamma symmetry Fact = 2 else Fact = 1
!
! Eh = Fact * omega * (sum over ig) * fpi / ( gv%hg_l(ig) * tpiba2 ) *
! Eh = Fact * omega * (sum over ig) * fpi / ( g(ig) * tpiba2 ) *
! { rho_e(ig) + rho_I(ig) } * conjugate { rho_e(ig) + rho_I(ig) }
! if Gamma symmetry Fact = 1 else Fact = 1/2
!
@ -941,13 +939,13 @@
! vloc(ig) = vloc_h(ig) + vloc_ps(ig)
!
! Local contribution to the forces on the ions
! eigrx(ig,isa) = ei1( gv%mill(1,ig), isa)
! eigry(ig,isa) = ei2( gv%mill(2,ig), isa)
! eigrz(ig,isa) = ei3( gv%mill(3,ig), isa)
! fpibg = fpi / ( gv%hg_l(ig) * tpiba2 )
! eigrx(ig,isa) = ei1( mill(1,ig), isa)
! eigry(ig,isa) = ei2( mill(2,ig), isa)
! eigrz(ig,isa) = ei3( mill(3,ig), isa)
! fpibg = fpi / ( g(ig) * tpiba2 )
! tx_h(ig,is) = fpibg * rhops(ig, is) * CONJG( rho_e(ig) + rho_I(ig) )
! tx_ps(ig,is) = vps(ig,is) * CONJG( rho_e(ig) )
! gx(ig) = CMPLX(0.D0, gv%gx_l(1,ig)) * tpiba
! gx(ig) = CMPLX(0.D0, gx(1,ig)) * tpiba
! fion(x,isa) = fion(x,isa) +
! Fact * omega * ( sum over ig, ispin) (tx_h(ig,is) + tx_ps(ig,is)) *
! gx(ig) * eigrx(ig,isa) * eigry(ig,isa) * eigrz(ig,isa)
@ -962,18 +960,19 @@
USE brillouin, ONLY: kpoints
USE charge_types, ONLY: charge_descriptor
USE atoms_type_module, ONLY: atoms_type
USE cp_types, ONLY: recvecs
USE pseudo_types, ONLY: pseudo_ncpp
USE io_global, ONLY: stdout
USE grid_dimensions, ONLY: nr1, nr2, nr3
USE reciprocal_vectors, ONLY: mill_l
USE ions_base, ONLY: nat
USE gvecp, ONLY: ngm
USE reciprocal_vectors, ONLY: gstart, gx, g
IMPLICIT NONE
! ... Arguments
TYPE (atoms_type) :: atoms
TYPE (recvecs), INTENT(in) :: gv
TYPE (pseudo_ncpp), INTENT(IN) :: ap(:)
TYPE (kpoints), INTENT(in) :: kp
TYPE (boxdimensions), INTENT(in) :: ht
@ -996,7 +995,7 @@
INTEGER :: is, ia, isa, ig, ig1, ig2, ig3, nspin, ispin
REAL(dbl) :: fpibg, cost, omega
COMPLEX(dbl) :: cxc, rhet, rhog, vp, rp, gx, gy, gz
COMPLEX(dbl) :: cxc, rhet, rhog, vp, rp, gxc, gyc, gzc
COMPLEX(dbl) :: teigr, cnvg, cvn, tx, ty, tz, vscreen
COMPLEX(dbl), ALLOCATABLE :: ftmp(:,:)
COMPLEX(dbl), ALLOCATABLE :: screen_coul(:)
@ -1010,8 +1009,8 @@
END IF
IF( tscreen ) THEN
ALLOCATE( screen_coul( gv%ng_l ) )
CALL cluster_bc( screen_coul, gv%hg_l, ht, desc )
ALLOCATE( screen_coul( ngm ) )
CALL cluster_bc( screen_coul, g, ht, desc )
END IF
!=======================================================================
@ -1024,19 +1023,19 @@
ehte = 0.0d0
ehti = 0.0d0
DO IG = gv%gstart, gv%ng_l
DO IG = gstart, ngm
RP = (0.D0,0.D0)
VP = (0.D0,0.D0)
DO IS = 1, atoms%nsp
VP = VP + sfac( is, IG ) * vps( ig, is )
RP = RP + sfac( is, IG ) * rhops( ig, is )
VP = VP + sfac( ig, is ) * vps( ig, is )
RP = RP + sfac( ig, is ) * rhops( ig, is )
END DO
! WRITE( stdout,*) 'vp ',ig, vp ! DEBUG
! WRITE( stdout,*) 'rp ',ig, rp ! DEBUG
! WRITE( stdout,*) 'rhoeg ',ig, SUM( RHOEG( ig, : ) ) ! DEBUG
! WRITE( stdout,*) 'mill ',ig, gv%mill(1,ig),gv%mill(2,ig),gv%mill(3,ig)
! WRITE( stdout,*) 'mill ',ig, mill(1,ig),mill(2,ig),mill(3,ig)
RHET = SUM( RHOEG( ig, : ) )
RHOG = RHET + RP
@ -1045,9 +1044,9 @@
! WRITE( stdout,*) 'rhog ',ig, rhog ! DEBUG
IF( tscreen ) THEN
FPIBG = fpi / ( gv%hg_l(ig) * tpiba2 ) + screen_coul(ig)
FPIBG = fpi / ( g(ig) * tpiba2 ) + screen_coul(ig)
ELSE
FPIBG = fpi / ( gv%hg_l(ig) * tpiba2 )
FPIBG = fpi / ( g(ig) * tpiba2 )
END IF
! WRITE( stdout,*) 'fpibg ',ig, fpibg ! DEBUG
@ -1059,19 +1058,19 @@
ehti = ehti + fpibg * REAL( rp * CONJG(rp))
IF(TTFORCE) THEN
ig1 = GV%mill(1,IG)
ig2 = GV%mill(2,IG)
ig3 = GV%mill(3,IG)
GX = CMPLX(0.D0,gv%gx_l(1,IG))
GY = CMPLX(0.D0,gv%gx_l(2,IG))
GZ = CMPLX(0.D0,gv%gx_l(3,IG))
ig1 = mill_l(1,IG)
ig2 = mill_l(2,IG)
ig3 = mill_l(3,IG)
GXC = CMPLX(0.D0,gx(1,IG))
GYC = CMPLX(0.D0,gx(2,IG))
GZC = CMPLX(0.D0,gx(3,IG))
isa = 1
DO IS = 1, atoms%nsp
CNVG = RHOPS(IG,is) * FPIBG * CONJG(rhog)
CVN = VPS(ig, is) * CONJG(rhet)
TX = (CNVG+CVN) * GX
TY = (CNVG+CVN) * GY
TZ = (CNVG+CVN) * GZ
TX = (CNVG+CVN) * GXC
TY = (CNVG+CVN) * GYC
TZ = (CNVG+CVN) * GZC
DO IA = 1, atoms%na(is)
TEIGR = ei1(IG1,ISA) * ei2(IG2,ISA) * ei3(IG3,ISA)
ftmp(1,ISA) = ftmp(1,ISA) + TEIGR*TX
@ -1095,7 +1094,7 @@
END IF
! ... G = 0 element
IF ( gv%gstart == 2 ) THEN
IF ( gstart == 2 ) THEN
vp = (0.D0,0.D0)
rp = (0.D0,0.D0)
IF( tscreen ) THEN
@ -1104,8 +1103,8 @@
vscreen = 0.0d0
END IF
DO IS = 1, atoms%nsp
vp = vp + sfac(is,1) * vps(1, is)
rp = rp + sfac(is,1) * rhops(1, is)
vp = vp + sfac( 1, is) * vps(1, is)
rp = rp + sfac( 1, is) * rhops(1, is)
END DO
rhet = SUM( rhoeg(1, :) )
rhog = rhet + rp
@ -1370,7 +1369,7 @@
! ----------------------------------------------
! BEGIN manual
SUBROUTINE self_vofloc(tscreen, ehte, vloc, rhoeg, gv, kp, ht, desc)
SUBROUTINE self_vofloc(tscreen, ehte, vloc, rhoeg, kp, ht, desc)
! adds the hartree part of the self interaction
!
@ -1381,13 +1380,13 @@
USE cell_module, ONLY: boxdimensions
USE brillouin, ONLY: kpoints
USE charge_types, ONLY: charge_descriptor
USE cp_types, ONLY: recvecs
USE cell_base, ONLY: tpiba2
USE gvecp, ONLY: ngm
USE reciprocal_vectors, ONLY: gstart, g
IMPLICIT NONE
! ... Arguments
TYPE (recvecs), INTENT(in) :: gv
TYPE (kpoints), INTENT(in) :: kp
TYPE (boxdimensions), INTENT(in) :: ht
TYPE (charge_descriptor), INTENT(IN) :: desc
@ -1408,8 +1407,8 @@
IF( tscreen ) THEN
ALLOCATE( screen_coul( gv%ng_l ) )
CALL cluster_bc( screen_coul, gv%hg_l, ht, desc )
ALLOCATE( screen_coul( ngm ) )
CALL cluster_bc( screen_coul, g, ht, desc )
END IF
!=======================================================================
@ -1419,13 +1418,13 @@
omega = ht%deth
ehte = 0.D0
DO IG = gv%gstart, gv%ng_l
DO IG = gstart, ngm
rhog = rhoeg(ig)
IF( tscreen ) THEN
FPIBG = fpi / ( gv%hg_l(ig) * tpiba2 ) + screen_coul(ig)
FPIBG = fpi / ( g(ig) * tpiba2 ) + screen_coul(ig)
ELSE
FPIBG = fpi / ( gv%hg_l(ig) * tpiba2 )
FPIBG = fpi / ( g(ig) * tpiba2 )
END IF
vloc(ig) = fpibg * rhog
@ -1436,7 +1435,7 @@
! ... G = 0 element
IF ( gv%gstart == 2 ) THEN
IF ( gstart == 2 ) THEN
IF( tscreen ) THEN
vscreen = screen_coul(1)
ELSE
@ -1460,7 +1459,7 @@
SUBROUTINE localisation( wfc, atoms_m, gv, kp, ht, desc)
SUBROUTINE localisation( wfc, atoms_m, kp, ht, desc)
! adds the hartree part of the self interaction
!
@ -1477,8 +1476,9 @@
USE sic_module, ONLY: rad_localisation, pos_localisation
USE ions_base, ONLY: ind_srt
USE fft_base, ONLY: dfftp
USE cp_types, ONLY: recvecs
USE cell_base, ONLY: tpiba2
USE reciprocal_vectors, ONLY: gstart, g
USE gvecp, ONLY: ngm
IMPLICIT NONE
@ -1486,7 +1486,6 @@
COMPLEX(dbl), INTENT(IN) :: wfc(:)
TYPE (atoms_type), INTENT(in) :: atoms_m
TYPE (recvecs), INTENT(in) :: gv
TYPE (kpoints), INTENT(in) :: kp
TYPE (boxdimensions), INTENT(in) :: ht
TYPE (charge_descriptor), INTENT(IN) :: desc
@ -1511,8 +1510,8 @@
IF( .FALSE. ) THEN
ALLOCATE( screen_coul( gv%ng_l ) )
CALL cluster_bc( screen_coul, gv%hg_l, ht, desc )
ALLOCATE( screen_coul( ngm ) )
CALL cluster_bc( screen_coul, g, ht, desc )
END IF
@ -1529,7 +1528,7 @@
ALLOCATE( density( nr1x, nr2x, nr3x ) )
ALLOCATE( psi( nr1x, nr2x, nr3x ) )
ALLOCATE( cpsi( nr1x, nr2x, nr3x ) )
ALLOCATE( k_density( gv%ng_l ) )
ALLOCATE( k_density( ngm ) )
CALL pw_invfft( cpsi(:,:,:), wfc(:), wfc(:) )
psi = REAL( cpsi, dbl )
@ -1591,13 +1590,13 @@
! ... G /= 0 elements
DO IG = gv%gstart, gv%ng_l
DO IG = gstart, ngm
rhog = k_density(ig)
IF( .FALSE. ) THEN
FPIBG = fpi / ( gv%hg_l(ig) * tpiba2 ) + screen_coul(ig)
FPIBG = fpi / ( g(ig) * tpiba2 ) + screen_coul(ig)
ELSE
FPIBG = fpi / ( gv%hg_l(ig) * tpiba2 )
FPIBG = fpi / ( g(ig) * tpiba2 )
END IF
ehte = ehte + fpibg * REAL(rhog * CONJG(rhog))
@ -1606,7 +1605,7 @@
! ... G = 0 element
IF ( gv%gstart == 2 ) THEN
IF ( gstart == 2 ) THEN
IF( .FALSE. ) THEN
vscreen = screen_coul(1)
ELSE

22
CPV/print_out.f90
View File

@ -583,15 +583,15 @@
!=----------------------------------------------------------------------------=!
SUBROUTINE print_sfac( gv, rhoe, desc, sfac )
SUBROUTINE print_sfac( rhoe, desc, sfac )
USE cp_types, ONLY: recvecs
USE mp_global, ONLY: mpime, nproc, group
USE mp, ONLY: mp_max, mp_get, mp_put
USE fft, ONLY : pfwfft, pinvfft
USE charge_types, ONLY: charge_descriptor
USE reciprocal_vectors, ONLY: ig_l2g, gx, g
USE gvecp, ONLY: ngm
TYPE (recvecs), INTENT(IN) :: gv
TYPE (charge_descriptor), INTENT(IN) :: desc
REAL(dbl), INTENT(IN) :: rhoe(:,:,:,:)
COMPLEX(dbl), INTENT(IN) :: sfac(:,:)
@ -605,10 +605,10 @@
COMPLEX(dbl), ALLOCATABLE :: sfac_rcv(:,:)
nspin = SIZE(rhoe,4)
nsp = SIZE(sfac,1)
ngx_l = gv%ng_l
nsp = SIZE(sfac,2)
ngx_l = ngm
CALL mp_max(ngx_l, group)
ALLOCATE(rhoeg(gv%ng_l,nspin))
ALLOCATE(rhoeg(ngm,nspin))
ALLOCATE(hg_rcv(ngx_l))
ALLOCATE(gx_rcv(3,ngx_l))
ALLOCATE(ig_rcv(ngx_l))
@ -623,15 +623,15 @@
END IF
DO ip = 1, nproc
CALL mp_get(ng, gv%ng_l, mpime, ionode_id, ip-1, ip)
CALL mp_get(hg_rcv(:), gv%hg_l(:), mpime, ionode_id, ip-1, ip)
CALL mp_get(gx_rcv(:,:), gv%gx_l(:,:), mpime, ionode_id, ip-1, ip)
CALL mp_get(ig_rcv(:), gv%ig(:), mpime, ionode_id, ip-1, ip)
CALL mp_get(ng, ngm, mpime, ionode_id, ip-1, ip)
CALL mp_get(hg_rcv(:), g(:), mpime, ionode_id, ip-1, ip)
CALL mp_get(gx_rcv(:,:), gx(:,:), mpime, ionode_id, ip-1, ip)
CALL mp_get(ig_rcv(:), ig_l2g(:), mpime, ionode_id, ip-1, ip)
DO ispin = 1, nspin
CALL mp_get( rhoeg_rcv(:,ispin), rhoeg(:,ispin), mpime, ionode_id, ip-1, ip)
END DO
DO is = 1, nsp
CALL mp_get( sfac_rcv(:,is), sfac(is,:), mpime, ionode_id, ip-1, ip)
CALL mp_get( sfac_rcv(:,is), sfac(:,is), mpime, ionode_id, ip-1, ip)
END DO
IF( ionode ) THEN
DO ig = 1, ng

67
CPV/pseudopot.f90
View File

@ -29,9 +29,9 @@
! INTEGER FUNCTION l2ind( lw, is )
! SUBROUTINE pseudopotential_setup(nsp,psfile)
! SUBROUTINE deallocate_pseudopotential
! SUBROUTINE formf(ht,gv,ps)
! SUBROUTINE nlin(gv,kp,wnl)
! SUBROUTINE nlin_stress(wnla,gv)
! SUBROUTINE formf(ht,ps)
! SUBROUTINE nlin(kp,wnl)
! SUBROUTINE nlin_stress(wnla)
! SUBROUTINE pseudo_wave_info(oc_out,nchan_out,mesh_out,dx_out, &
! rw_out,rps_out)
! REAL(dbl) FUNCTION zvpseudo(is)
@ -41,7 +41,7 @@
! ... declare modules
USE kinds
USE parameters, ONLY: cp_lmax
USE cp_types, ONLY: pseudo, allocate_pseudo, recvecs
USE cp_types, ONLY: pseudo, allocate_pseudo
USE splines, ONLY: spline_data
USE read_pseudo_module_fpmd, ONLY: nspnl, l2ind
@ -250,17 +250,18 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE pseudopotential_init( ps, na, nsp, gv, kp )
SUBROUTINE pseudopotential_init( ps, na, nsp, kp )
! ... declare modules
USE brillouin, ONLY: kpoints
USE pseudo_types, ONLY: pseudo_ncpp, pseudo_upf
USE read_pseudo_module_fpmd, ONLY: ap
USE gvecp, ONLY: ngm
USE gvecw, ONLY: ngw
IMPLICIT NONE
TYPE (pseudo) :: ps
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
INTEGER, INTENT(IN) :: na(:), nsp
@ -273,7 +274,7 @@
tcc = .FALSE.
WHERE ( ap(1:nsp)%tnlcc ) tcc = .TRUE.
CALL allocate_pseudo(ps, nsp, gv%ng_l, gv%ngw_l, kp%nkpt, lnlx, ngh, tcc)
CALL allocate_pseudo(ps, nsp, ngm, ngw, kp%nkpt, lnlx, ngh, tcc)
ps%ap => ap
@ -337,16 +338,16 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE formf(ht, gv, kp, ps)
SUBROUTINE formf(ht, kp, ps)
! this routine computes:
! the form factors of:
! pseudopotential (gv%vps)
! ionic pseudocharge (gv%rhops)
! core corrections to the pseudopotential (gv%cc(:)%rhoc1)
! pseudopotential (vps)
! ionic pseudocharge (rhops)
! core corrections to the pseudopotential (cc(:)%rhoc1)
! the derivatives with respect to cell degrees of freedom of:
! pseudopotential form factors (gv%dvps)
! core correction to the pseudopotential (gv%cc(:)%rhocp)
! pseudopotential form factors (dvps)
! core correction to the pseudopotential (cc(:)%rhocp)
! ----------------------------------------------
! ... declare modules
@ -359,12 +360,12 @@
USE pseudo_base, ONLY: formfn_base
USE pseudo_base, ONLY: corecor_base
USE pseudo_base, ONLY: rhops_base
USE reciprocal_vectors, ONLY: g
IMPLICIT NONE
! ... declare subroutine arguments
TYPE (pseudo) ps
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
TYPE (boxdimensions), INTENT(IN) :: ht
@ -385,12 +386,12 @@
IF(tpstab .AND. (.NOT.tpstab_first)) THEN
! ... check the consistency of the tab with respect the g vectors
tpstab_first = chkpstab(gv%hg_l, xgtabmax)
tpstab_first = chkpstab(g, xgtabmax)
END IF
IF(tpstab .AND. tpstab_first) THEN
! ... build the pseudopotential tabs
CALL build_pstab(gv)
CALL build_pstab()
tpstab_first = .FALSE.
END IF
@ -398,23 +399,23 @@
! ... handle local part
DO is = 1, nsp
CALL rhops_base(ps%ap(is), gv%hg_l, ps%rhops(:,is), omega)
CALL rhops_base(ps%ap(is), g, ps%rhops(:,is), omega)
IF( ps%ap(is)%tnlcc ) THEN
IF(tpstab) THEN
CALL corecortab_base(gv%hg_l, ps%rhoc1(:,is), ps%rhocp(:,is), &
CALL corecortab_base(g, ps%rhoc1(:,is), ps%rhocp(:,is), &
rhoc1_sp(is), rhocp_sp(is), xgtabmax, omega)
ELSE
CALL corecor_base(ps%ap(is), gv%hg_l, ps%rhoc1(:,is), ps%rhocp(:,is), omega)
CALL corecor_base(ps%ap(is), g, ps%rhoc1(:,is), ps%rhocp(:,is), omega)
END IF
END IF
! ... numeric pseudopotential
IF(tpstab) THEN
CALL formftab_base(gv%hg_l, ps%vps(:,is), ps%dvps(:,is), &
CALL formftab_base(g, ps%vps(:,is), ps%dvps(:,is), &
vps_sp(is), dvps_sp(is), xgtabmax, omega )
ELSE
CALL formfn_base(ps%ap(is), gv%hg_l, ps%vps(:,is), ps%dvps(:,is), omega)
CALL formfn_base(ps%ap(is), g, ps%vps(:,is), ps%dvps(:,is), omega)
END IF
! DEBUG
@ -433,7 +434,7 @@
! ... handle nonlocal part
CALL nlin(gv, kp, ps%wnl)
CALL nlin(kp, ps%wnl)
RETURN
@ -443,7 +444,7 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE build_pstab(gv)
SUBROUTINE build_pstab()
USE ions_base, ONLY: nsp
USE constants, ONLY: pi, fpi
@ -457,10 +458,10 @@
USE pseudo_base, ONLY: corecor_base
USE pseudo_types, ONLY: pseudo_ncpp, pseudo_upf
USE read_pseudo_module_fpmd, ONLY: ap
USE reciprocal_vectors, ONLY: g
IMPLICIT NONE
TYPE (recvecs), INTENT(IN) :: gv
REAL(dbl), ALLOCATABLE :: fintl(:,:)
INTEGER :: ig, is, mmax, lloc, nval, l, ll
REAL(dbl) :: xg, xgmax, xgmin, dxg, res
@ -475,7 +476,7 @@
nval = pstab_size
xgmin = 0.0d0
xgmax = tpiba * SQRT( MAXVAL( gv%hg_l ) )
xgmax = tpiba * SQRT( MAXVAL( g ) )
CALL mp_max(xgmax, group)
xgmax = xgmax + (xgmax-xgmin)
dxg = (xgmax - xgmin) / REAL(nval-1)
@ -540,7 +541,7 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE nlin(gv, kp, wnl)
SUBROUTINE nlin(kp, wnl)
! this routine computes the temporary arrays twnl
! to be used by nlrh and dforce
@ -552,11 +553,11 @@
USE pseudotab_base, ONLY: nlintab_base
USE pseudo_base, ONLY: nlin_base
USE read_pseudo_module_fpmd, ONLY: ap
USE reciprocal_space_mesh, ONLY: gk_l
IMPLICIT NONE
! ... declare subroutine arguments
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
REAL(dbl) :: wnl(:,:,:,:)
@ -572,9 +573,9 @@
DO is = 1, nspnl
DO ik = 1, kp%nkpt
IF( tpstab ) THEN
CALL nlintab_base(gv%khg_l(:,ik), wnl(:,:,is,ik), ap(is)%lnl, wnl_sp(:,is), xgtabmax)
CALL nlintab_base(gk_l(:,ik), wnl(:,:,is,ik), ap(is)%lnl, wnl_sp(:,is), xgtabmax)
ELSE
CALL nlin_base(ap(is), gv%khg_l(:,ik), wnl(:,:,is,ik))
CALL nlin_base(ap(is), gk_l(:,ik), wnl(:,:,is,ik))
END IF
END DO
END DO
@ -584,7 +585,7 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE nlin_stress(wnla,gv)
SUBROUTINE nlin_stress(wnla)
! this routine computes the temporary arrays twnl
! to be used by nlrh and dforce.
@ -595,11 +596,11 @@
USE pseudotab_base, ONLY: nlintab_base
USE pseudo_base, ONLY: nlin_stress_base
USE read_pseudo_module_fpmd, ONLY: ap
USE reciprocal_vectors, ONLY: g
IMPLICIT NONE
! ... declare subroutine arguments
TYPE (recvecs) :: gv
REAL(dbl) :: wnla(:,:,:)
INTEGER :: is
@ -608,9 +609,9 @@
wnla = 0.0d0
DO IS = 1, NSPNL
IF (tpstab) THEN
CALL nlintab_base(gv%hg_l, wnla(:,:,is), ap(is)%lnl, wnla_sp(:,is), xgtabmax)
CALL nlintab_base(g, wnla(:,:,is), ap(is)%lnl, wnla_sp(:,is), xgtabmax)
ELSE
CALL nlin_stress_base(ap(is), gv%hg_l, wnla(:,:,is))
CALL nlin_stress_base(ap(is), g, wnla(:,:,is))
END IF
END DO
RETURN

25
CPV/restart.f90
View File

@ -224,11 +224,10 @@
SUBROUTINE writefile_fpmd( nfi, trutime, c0, cm, cdesc, occ, &
atoms_0, atoms_m, acc, taui, cdmi, ibrav, celldm, &
ht_m, ht_0, rho, desc, vpot, gv, kp)
atoms_0, atoms_m, acc, taui, cdmi, &
ht_m, ht_0, rho, desc, vpot, kp)
use electrons_module, only: nspin
use cp_types, only: recvecs
USE cell_module, only: boxdimensions, r_to_s
USE brillouin, only: kpoints
USE wave_types, ONLY: wave_descriptor
@ -243,15 +242,15 @@
USE cell_nose, ONLY: xnhh0, xnhhm, vnhh
USE ions_nose, ONLY: vnhp, xnhp0, xnhpm, nhpcl
USE cp_restart, ONLY: cp_writefile, cp_write_wfc
USE reciprocal_space_mesh, ONLY: b1, b2, b3
IMPLICIT NONE
INTEGER, INTENT(IN) :: nfi, ibrav
INTEGER, INTENT(IN) :: nfi
COMPLEX(dbl), INTENT(IN) :: c0(:,:,:,:), cm(:,:,:,:)
REAL(dbl), INTENT(IN) :: occ(:,:,:)
TYPE (kpoints), INTENT(IN) :: kp
TYPE (boxdimensions), INTENT(IN) :: ht_m, ht_0
TYPE (recvecs), INTENT(IN) :: gv
TYPE (atoms_type), INTENT(IN) :: atoms_0, atoms_m
REAL(dbl), INTENT(IN) :: rho(:,:,:,:)
TYPE (charge_descriptor), INTENT(IN) :: desc
@ -259,7 +258,6 @@
REAL(dbl), INTENT(INOUT) :: vpot(:,:,:,:)
REAL(dbl), INTENT(IN) :: taui(:,:)
REAL(dbl), INTENT(IN) :: celldm(:)
REAL(dbl), INTENT(IN) :: acc(:), cdmi(:)
REAL(dbl), INTENT(IN) :: trutime
@ -284,7 +282,7 @@
CALL cp_writefile( ndw, ' ', .TRUE., nfi, trutime, acc, kp%nkpt, kp%xk, kp%weight, &
ht_0%a, ht_m%a, ht_0%hvel, ht_0%gvel, xnhh0, xnhhm, vnhh, taui, cdmi, &
atoms_0%taus, atoms_0%vels, atoms_m%taus, atoms_m%vels, atoms_0%for, vnhp, &
xnhp0, xnhpm, nhpcl, occ, occ, lambda, lambda, gv%b1, gv%b2, gv%b3, &
xnhp0, xnhpm, nhpcl, occ, occ, lambda, lambda, b1, b2, b3, &
xnhe0, xnhem, vnhe, ekincm, mat_z )
DEALLOCATE( lambda )
@ -308,11 +306,10 @@
!=----------------------------------------------------------------------------=!
SUBROUTINE readfile_fpmd( nfi, trutime, &
c0, cm, cdesc, occ, atoms_0, atoms_m, acc, taui, cdmi, ibrav, celldm, &
ht_m, ht_0, rho, desc, vpot, gv, kp)
c0, cm, cdesc, occ, atoms_0, atoms_m, acc, taui, cdmi, &
ht_m, ht_0, rho, desc, vpot, kp)
use electrons_base, only: nbsp
use cp_types, ONLY: recvecs
USE cell_module, only: boxdimensions, cell_init, r_to_s, s_to_r
USE brillouin, only: kpoints
use parameters, only: npkx, nsx
@ -337,16 +334,15 @@
USE cell_nose, ONLY: xnhh0, xnhhm, vnhh
USE ions_nose, ONLY: vnhp, xnhp0, xnhpm, nhpcl
USE cp_restart, ONLY: cp_readfile, cp_read_wfc
USE reciprocal_space_mesh, ONLY: b1, b2, b3
IMPLICIT NONE
INTEGER, INTENT(INOUT) :: ibrav
INTEGER, INTENT(OUT) :: nfi
COMPLEX(dbl), INTENT(INOUT) :: c0(:,:,:,:), cm(:,:,:,:)
REAL(dbl), INTENT(INOUT) :: occ(:,:,:)
TYPE (kpoints), INTENT(INOUT) :: kp
TYPE (boxdimensions), INTENT(INOUT) :: ht_m, ht_0
TYPE (recvecs), INTENT(INOUT) :: gv
TYPE (atoms_type), INTENT(INOUT) :: atoms_0, atoms_m
REAL(dbl), INTENT(INOUT) :: rho(:,:,:,:)
TYPE (charge_descriptor), INTENT(IN) :: desc
@ -354,7 +350,6 @@
REAL(dbl), INTENT(INOUT) :: vpot(:,:,:,:)
REAL(dbl), INTENT(OUT) :: taui(:,:)
REAL(dbl), INTENT(INOUT) :: celldm(:)
REAL(dbl), INTENT(OUT) :: acc(:), cdmi(:)
REAL(dbl), INTENT(OUT) :: trutime
@ -378,8 +373,8 @@
CALL cp_readfile( ndr, ' ', .TRUE., nfi, trutime, acc, kp%nkpt, kp%xk, kp%weight, &
hp0_ , hm1_ , hvel_ , gvel_ , xnhh0, xnhhm, vnhh, taui, cdmi, &
atoms_0%taus, atoms_0%vels, atoms_m%taus, atoms_m%vels, atoms_0%for, vnhp, &
xnhp0, xnhpm, nhpcl, occ, occ, lambda_ , lambda_ , gv%b1, gv%b2, &
gv%b3, xnhe0, xnhem, vnhe, ekincm, mat_z_ , tens )
xnhp0, xnhpm, nhpcl, occ, occ, lambda_ , lambda_ , b1, b2, &
b3, xnhe0, xnhem, vnhe, ekincm, mat_z_ , tens )
DEALLOCATE( lambda_ )

39
CPV/restart_sub.f90
View File

@ -46,7 +46,8 @@ CONTAINS
ef_tune, wf_options, ef_potential
USE core, ONLY: nlcc_any
USE gvecw, ONLY: ngw
USE reciprocal_vectors, ONLY: gstart
USE gvecs, ONLY: ngs
USE reciprocal_vectors, ONLY: gstart, mill_l
USE wave_base, ONLY: wave_verlet
USE cvan, ONLY: nvb
USE ions_nose, ONLY: xnhp0, xnhpm, vnhp
@ -58,6 +59,7 @@ CONTAINS
USE time_step, ONLY: delt
use electrons_nose, only: xnhe0, xnhem, vnhe
USE cell_nose, ONLY: xnhh0, xnhhm, vnhh, cell_nosezero
use phase_factors_module, only: strucf
COMPLEX(kind=8) :: eigr(:,:), ei1(:,:), ei2(:,:), ei3(:,:)
COMPLEX(kind=8) :: eigrb(:,:)
@ -128,7 +130,7 @@ CONTAINS
end if
call phfac( tau0, ei1, ei2, ei3, eigr )
call strucf( ei1, ei2, ei3, sfac )
call strucf( sfac, ei1, ei2, ei3, mill_l, ngs )
call formf( tfirst, eself )
call calbec ( 1, nsp, eigr, c0, bec )
if (tpre) call caldbec( 1, nsp, eigr, c0 )
@ -282,6 +284,10 @@ CONTAINS
use io_global, only: stdout
use cell_base, only: s_to_r
use cpr_subroutines, only: print_atomic_var
USE reciprocal_vectors, ONLY: gstart, mill_l
USE gvecs, ONLY: ngs
use phase_factors_module, only: strucf
IMPLICIT NONE
@ -317,7 +323,7 @@ CONTAINS
endif
!
call phfac(tau0,ei1,ei2,ei3,eigr)
call strucf(ei1,ei2,ei3,sfac)
call strucf( sfac, ei1, ei2, ei3, mill_l, ngs )
call formf(tfirst,eself)
call calbec (1,nsp,eigr,c0,bec)
if (tpre) call caldbec(1,nsp,eigr,c0)
@ -328,14 +334,14 @@ CONTAINS
!=----------------------------------------------------------------------------=!
SUBROUTINE from_restart_fpmd( nfi, acc, gv, kp, ps, rhoe, desc, cm, c0, cdesc, &
SUBROUTINE from_restart_fpmd( nfi, acc, kp, ps, rhoe, desc, cm, c0, cdesc, &
eigr, ei1, ei2, ei3, sfac, fi, ht_m, ht_0, atoms_m, atoms_0, fnl, vpot, edft)
! this routine recreates the starting configuration from a restart file
! ... declare modules
USE kinds, ONLY: dbl
USE phase_factors_module, ONLY: strucf
USE phase_factors_module, ONLY: strucf, phfacs
USE time_step, ONLY: delt
USE reciprocal_space_mesh, ONLY: newg
USE charge_density, ONLY: rhoofr
@ -346,7 +352,7 @@ CONTAINS
USE ions_module, ONLY: set_reference_positions, print_scaled_positions, &
constraints_setup, set_velocities
USE energies, ONLY: dft_energy_type
USE cp_types, ONLY: recvecs, pseudo
USE cp_types, ONLY: pseudo
USE pseudopotential, ONLY: formf
USE cell_module, only: boxdimensions, gethinv, alat
USE cell_base, ONLY: r_to_s, s_to_r
@ -371,6 +377,9 @@ CONTAINS
USE charge_types, ONLY: charge_descriptor
USE ions_base, ONLY: vel_srt, tau_units
USE runcp_module, ONLY: runcp_ncpp
use grid_dimensions, only: nr1, nr2, nr3
USE reciprocal_vectors, ONLY: mill_l
USE gvecp, ONLY: ngm
IMPLICIT NONE
@ -384,7 +393,6 @@ CONTAINS
COMPLEX(dbl) :: ei1(:,:)
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
TYPE (recvecs) :: gv
TYPE (kpoints) :: kp
COMPLEX(dbl), INTENT(INOUT) :: cm(:,:,:,:), c0(:,:,:,:)
REAL(dbl) :: fi(:,:,:)
@ -415,11 +423,11 @@ CONTAINS
! ... with the g vectors modules
! ... if tbeg is false, ht_0 is read from the restart file and now
! ... we have to compute the inverse and the volume of the cell,
! ... together with the new reciprocal vectors (gv)
! ... together with the new reciprocal vectors
IF ( .NOT. tbeg ) THEN
CALL newg( gv, kp, ht_0%m1 )
CALL newg( kp, ht_0%m1 )
CALL newgb( ht_0%hmat(:,1), ht_0%hmat(:,2), ht_0%hmat(:,3), ht_0%omega, alat )
!
@ -499,8 +507,9 @@ CONTAINS
! ... computes form factors and initializes nl-pseudop. according
! ... to starting cell (from ndr or again fort.10)
CALL strucf(sfac, atoms_0, eigr, ei1, ei2, ei3, gv)
CALL formf(ht_0, gv, kp, ps)
CALL phfacs( ei1, ei2, ei3, eigr, mill_l, atoms_0%taus, nr1, nr2, nr3, atoms_0%nat )
CALL strucf( sfac, ei1, ei2, ei3, mill_l, ngm )
CALL formf(ht_0, kp, ps)
IF( tzeroe .OR. tzerop ) THEN
@ -512,16 +521,16 @@ CONTAINS
atoms_0%for = 0.0d0
edft%enl = nlrh_m(c0, cdesc, ttforce, atoms_0, fi, gv, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL rhoofr(gv, kp, c0, cdesc, fi, rhoe, desc, ht_0)
CALL vofrhos( ( iprsta > 1), rhoe, desc, tfor, thdyn, ttforce, atoms_0, gv, kp, fnl, vpot, ps, &
edft%enl = nlrh_m(c0, cdesc, ttforce, atoms_0, fi, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL rhoofr(kp, c0, cdesc, fi, rhoe, desc, ht_0)
CALL vofrhos( ( iprsta > 1), rhoe, desc, tfor, thdyn, ttforce, atoms_0, kp, fnl, vpot, ps, &
c0, cdesc, fi, eigr, ei1, ei2, ei3, sfac, timepre, ht_0, edft)
IF( tzeroe ) THEN
IF( tcarpar .AND. ( .NOT. force_pairing ) ) THEN
CALL runcp_ncpp( cm, cm, c0, cdesc, gv, kp, ps, vpot, eigr, fi, fnl, vdum, &
CALL runcp_ncpp( cm, cm, c0, cdesc, kp, ps, vpot, eigr, fi, fnl, vdum, &
gam, cgam, restart = .TRUE.)
IF(tortho) THEN

24
CPV/runcg.f90
View File

@ -54,7 +54,7 @@
! -----------------------------------------------------------------------
! BEGIN manual
SUBROUTINE runcg_new(tortho, tprint, rhoe, desc, atoms_0, gv, kp, &
SUBROUTINE runcg_new(tortho, tprint, rhoe, desc, atoms_0, kp, &
ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht0, occ, ei, &
fnl, vpot, doions, edft, maxnstep, cgthr )
@ -87,6 +87,7 @@
USE charge_types, ONLY: charge_descriptor
USE control_flags, ONLY: force_pairing
USE environment, ONLY: start_cclock_val
USE reciprocal_space_mesh, ONLY: gkmask_l
IMPLICIT NONE
@ -103,7 +104,6 @@
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
COMPLEX(dbl) :: sfac(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
TYPE (boxdimensions), INTENT(INOUT) :: ht0
REAL(dbl) :: occ(:,:,:)
@ -183,7 +183,7 @@
s1 = cclock()
CALL kspotential( ttprint, ttforce, ttstress, rhoe, desc, &
atoms_0, gv, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0, occ, fnl, vpot, edft, timepre )
atoms_0, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0, occ, fnl, vpot, edft, timepre )
s2 = cclock()
@ -193,7 +193,7 @@
! ... |d H / dPsi_j > = H |Psi_j> - Sum{i} <Psi_i|H|Psi_j> |Psi_i>
CALL dforce_all( ispin, c0(:,:,:,ispin), cdesc, occ(:,:,ispin), cp(:,:,:,ispin), &
gv, vpot(:,:,:,ispin), fnl(:,ispin), eigr, ps)
vpot(:,:,:,ispin), fnl(:,ispin), eigr, ps)
! ... Project the gradient
IF( gamma_symmetry ) THEN
@ -242,7 +242,7 @@
! perform line minimization in the direction of "hacca"
CALL CGLINMIN(emin, demin, tbad, edft, cp, c0, cdesc, occ, vpot, rhoe, desc, hacca, &
atoms_0, ht0, fnl, ps, eigr, ei1, ei2, ei3, sfac, gv, kp)
atoms_0, ht0, fnl, ps, eigr, ei1, ei2, ei3, sfac, kp)
! CALL print_energies( edft )
s5 = cclock()
@ -258,7 +258,7 @@
cp = c0 + cm
CALL fixwave( cp, cdesc, gv%kg_mask_l )
CALL fixwave( cp, cdesc, gkmask_l )
IF( tortho ) THEN
CALL ortho( c0, cp, cdesc, pmss, emass )
@ -271,7 +271,7 @@
ekinc = 0.0d0
DO ispin = 1, nspin
ekinc = ekinc + cp_kinetic_energy( ispin, cp(:,:,:,ispin), c0(:,:,:,ispin), cdesc, kp, &
gv%kg_mask_l, pmss, delt)
gkmask_l, pmss, delt)
END DO
IF( iter > 1 ) THEN
dek = ekinc - ekinc_old
@ -317,7 +317,7 @@
DO ispin = 1, nspin
CALL dforce_all( ispin, c0(:,:,:,ispin), cdesc, occ(:,:,ispin), hacca(:,:,:,ispin), &
gv, vpot(:,:,:,ispin), fnl(:,ispin), eigr, ps)
vpot(:,:,:,ispin), fnl(:,ispin), eigr, ps)
nb_g( ispin ) = cdesc%nbt( ispin )
@ -358,7 +358,7 @@
! ---------------------------------------------------------------------- !
SUBROUTINE CGLINMIN(emin, ediff, tbad, edft, cp, c, cdesc, occ, vpot, rhoe, desc, hacca, &
atoms, ht, fnl, ps, eigr, ei1, ei2, ei3, sfac, gv, kp)
atoms, ht, fnl, ps, eigr, ei1, ei2, ei3, sfac, kp)
! ... declare modules
@ -374,6 +374,7 @@
USE potentials, ONLY: kspotential
USE atoms_type_module, ONLY: atoms_type
USE charge_types, ONLY: charge_descriptor
USE reciprocal_space_mesh, ONLY: gkmask_l
IMPLICIT NONE
@ -392,7 +393,6 @@
COMPLEX(dbl) :: ei1(:,:)
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
TYPE (boxdimensions), INTENT(INOUT) :: ht
REAL(dbl) :: occ(:,:,:)
@ -595,11 +595,11 @@
cp = c + hstep * hacca
CALL fixwave( cp, cdesc, gv%kg_mask_l )
CALL fixwave( cp, cdesc, gkmask_l )
CALL gram( cp, cdesc )
CALL kspotential( ttprint, ttforce, ttstress, rhoe, desc, &
atoms, gv, kp, ps, eigr, ei1, ei2, ei3, sfac, cp, cdesc, tcel, ht, occ, fnl, vpot, edft, timepre )
atoms, kp, ps, eigr, ei1, ei2, ei3, sfac, cp, cdesc, tcel, ht, occ, fnl, vpot, edft, timepre )
cgenergy = edft%etot

18
CPV/runcg_ion.f90
View File

@ -31,7 +31,7 @@
! BEGIN manual
SUBROUTINE runcg_ion(nfi, tortho, tprint, rhoe, desc, atomsp, atoms0, atomsm, &
gv, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht, occ, ei, &
kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht, occ, ei, &
fnl, vpot, doions, edft, etol, ftol, maxiter, sdthr, maxnstep )
! this routine computes the equilibrium ionic positions via conjugate gradient
@ -49,7 +49,7 @@
USE io_global, ONLY: stdout
USE cell_module, ONLY: boxdimensions, s_to_r, r_to_s
USE brillouin, ONLY: kpoints
USE cp_types, ONLY: recvecs, pseudo
USE cp_types, ONLY: pseudo
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
USE potentials, ONLY: kspotential
@ -79,7 +79,6 @@
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
COMPLEX(dbl) :: sfac(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
TYPE (boxdimensions), INTENT(INOUT) :: ht
REAL(dbl) :: occ(:,:,:)
@ -163,7 +162,7 @@
s1 = cclock()
old_clock_value = s1
CALL runsd(ttortho, ttprint, ttforce, rhoe, desc, atoms0, gv, kp, &
CALL runsd(ttortho, ttprint, ttforce, rhoe, desc, atoms0, kp, &
ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht, occ, ei, &
fnl, vpot, doions, edft, maxnstep, sdthr )
@ -192,7 +191,7 @@
WRITE( stdout,fmt="(/,8X,'cgion: iter',I5,' line minimization along gradient starting')") iter
CALL CGLINMIN(fret, edft, cp, c0, cm, cdesc, occ, ei, vpot, rhoe, desc, xi, atomsp, atoms0, &
ht, fnl, ps, eigr, ei1, ei2, ei3, sfac, gv, kp, maxnstep, sdthr, displ)
ht, fnl, ps, eigr, ei1, ei2, ei3, sfac, kp, maxnstep, sdthr, displ)
IF( tbad ) THEN
! displ = displ * 2.0d0
@ -206,7 +205,7 @@
IF( ionode ) WRITE( stdout, fmt='(8X,"cgion: bad step")') ! perform steepest descent
displ = displ / 2.0d0
CALL runsd(ttortho, ttprint, ttforce, rhoe, desc, atoms0, gv, kp, &
CALL runsd(ttortho, ttprint, ttforce, rhoe, desc, atoms0, kp, &
ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht, occ, ei, &
fnl, vpot, doions, edft, maxnstep, sdthr )
@ -298,12 +297,12 @@
! ---------------------------------------------------------------------- !
SUBROUTINE CGLINMIN(emin, edft, cp, c0, cm, cdesc, occ, ei, vpot, &
rhoe, desc, hacca, atomsp, atoms0, ht, fnl, ps, eigr, ei1, ei2, ei3, sfac, gv, kp, &
rhoe, desc, hacca, atomsp, atoms0, ht, fnl, ps, eigr, ei1, ei2, ei3, sfac, kp, &
maxnstep, sdthr, displ)
! ... declare modules
USE cp_types, ONLY: recvecs, pseudo
USE cp_types, ONLY: pseudo
USE wave_types, ONLY: wave_descriptor
USE brillouin, ONLY: kpoints
USE pseudo_projector, ONLY: projector
@ -336,7 +335,6 @@
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
COMPLEX(dbl) :: sfac(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
TYPE (boxdimensions), INTENT(INOUT) :: ht
REAL(dbl) :: occ(:,:,:)
@ -561,7 +559,7 @@
! ... Calculate Forces (fion) and DFT Total Energy (edft) for the new ionic
! ... positions (atomsp)
CALL runsd(ttortho, ttprint, ttforce, rhoe, desc, atomsp, gv, kp, &
CALL runsd(ttortho, ttprint, ttforce, rhoe, desc, atomsp, kp, &
ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht, occ, ei, &
fnl, vpot, doions, edft, maxnstep, sdthr )

69
CPV/runcp.f90
View File

@ -27,7 +27,7 @@
! ----------------------------------------------
! BEGIN manual
SUBROUTINE runcp( ttprint, tortho, tsde, cm, c0, cp, cdesc, gv, kp, ps, &
SUBROUTINE runcp( ttprint, tortho, tsde, cm, c0, cp, cdesc, kp, ps, &
vpot, eigr, fi, ekinc, timerd, timeorto, ht, ei, fnl, vnosee )
! This subroutine performs a Car-Parrinello or Steepest-Descent step
@ -49,7 +49,7 @@
USE wave_functions, ONLY : rande, cp_kinetic_energy, gram
USE wave_base, ONLY : frice
USE wave_base, ONLY: hpsi
USE cp_types, ONLY: recvecs, pseudo
USE cp_types, ONLY: pseudo
USE cell_module, ONLY: boxdimensions
USE time_step, ONLY: delt
USE forces, ONLY: dforce
@ -60,6 +60,7 @@
USE control_flags, ONLY: tnosee
USE control_flags, ONLY: tdamp
USE wave_constrains, ONLY: update_lambda
USE reciprocal_space_mesh, ONLY: gkmask_l
IMPLICIT NONE
@ -70,7 +71,6 @@
TYPE (wave_descriptor), INTENT(IN) :: cdesc
TYPE (pseudo), INTENT(IN) :: ps
COMPLEX(dbl) :: eigr(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
REAL(dbl), INTENT(IN) :: fi(:,:,:)
TYPE (boxdimensions), INTENT(IN) :: ht
@ -110,7 +110,7 @@
! Compute electronic forces and move electrons
CALL runcp_ncpp( cm, c0, cp, cdesc, gv, kp, ps, vpot, eigr, fi, fnl, vnosee, &
CALL runcp_ncpp( cm, c0, cp, cdesc, kp, ps, vpot, eigr, fi, fnl, vnosee, &
gam, cgam, lambda = ttprint )
! Compute eigenstate
@ -143,7 +143,7 @@
DO is = 1, cdesc%nspin
ekinc(is) = cp_kinetic_energy( is, cp(:,:,:,is), cm(:,:,:,is), cdesc, kp, &
gv%kg_mask_l, pmss, delt)
gkmask_l, pmss, delt)
END DO
DEALLOCATE( cgam, gam, STAT=ierr)
@ -160,7 +160,7 @@
! ----------------------------------------------
! BEGIN manual
SUBROUTINE runcp_ncpp( cm, c0, cp, cdesc, gv, kp, ps, &
SUBROUTINE runcp_ncpp( cm, c0, cp, cdesc, kp, ps, &
vpot, eigr, fi, fnl, vnosee, gam, cgam, lambda, fromscra, diis, restart )
! This subroutine performs a Car-Parrinello or Steepest-Descent step
@ -179,7 +179,7 @@
USE electrons_module, ONLY: pmss
USE cp_electronic_mass, ONLY: emass
USE wave_base, ONLY: frice, wave_steepest, wave_verlet
USE cp_types, ONLY: recvecs, pseudo
USE cp_types, ONLY: pseudo
USE time_step, ONLY: delt
USE forces, ONLY: dforce
USE brillouin, ONLY: kpoints
@ -187,6 +187,7 @@
USE wave_constrains, ONLY: update_lambda
USE control_flags, ONLY: tnosee, tdamp, tsde
USE pseudo_projector, ONLY: projector
USE reciprocal_space_mesh, ONLY: gkmask_l
IMPLICIT NONE
@ -198,7 +199,6 @@
TYPE (wave_descriptor), INTENT(IN) :: cdesc
TYPE (pseudo), INTENT(IN) :: ps
COMPLEX(dbl) :: eigr(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
REAL(dbl), INTENT(IN) :: fi(:,:,:)
TYPE (projector) :: fnl(:,:)
@ -291,12 +291,12 @@
DO i = 1, nb, 2
IF ( cdesc%gamma ) THEN
CALL dforce( ik, i, c0(:,:,:,is), cdesc, fi(:,:,is), c2, c3, gv, vpot(:,:,:,is), &
CALL dforce( ik, i, c0(:,:,:,is), cdesc, fi(:,:,is), c2, c3, vpot(:,:,:,is), &
fnl(ik, is)%r(:,:,:), eigr, ps )
ELSE
CALL dforce( ik, i, c0(:,:,:,is), cdesc, fi(:,:,is), c2, gv, vpot(:,:,:,is), &
CALL dforce( ik, i, c0(:,:,:,is), cdesc, fi(:,:,is), c2, vpot(:,:,:,is), &
fnl(ik, is)%c, eigr, ps )
CALL dforce( ik, i+1, c0(:,:,:,is), cdesc, fi(:,:,is), c3, gv, vpot(:,:,:,is), &
CALL dforce( ik, i+1, c0(:,:,:,is), cdesc, fi(:,:,is), c3, vpot(:,:,:,is), &
fnl(ik, is)%c, eigr, ps )
END IF
@ -327,8 +327,8 @@
CALL wave_verlet( cp(:,i+1,ik,is), c0(:,i+1,ik,is), svar1, svar2, svar3, c3 )
END IF
IF( .NOT. cdesc%gamma ) THEN
cp(:,i,ik,is) = cp(:,i,ik,is) * gv%kg_mask_l(:,ik)
cp(:,i+1,ik,is) = cp(:,i+1,ik,is) * gv%kg_mask_l(:,ik)
cp(:,i,ik,is) = cp(:,i,ik,is) * gkmask_l(:,ik)
cp(:,i+1,ik,is) = cp(:,i+1,ik,is) * gkmask_l(:,ik)
ELSE
IF( cdesc%gzero ) cp(1,i,ik,is) = REAL( cp(1,i,ik,is), dbl )
IF( cdesc%gzero ) cp(1,i+1,ik,is) = REAL( cp(1,i+1,ik,is), dbl )
@ -339,10 +339,10 @@
IF( MOD(nx,2) /= 0) THEN
nb = nx
IF ( cdesc%gamma ) THEN
CALL dforce( ik, nb, c0(:,:,:,is), cdesc, fi(:,:,is), c2, gv, vpot(:,:,:,is), &
CALL dforce( ik, nb, c0(:,:,:,is), cdesc, fi(:,:,is), c2, vpot(:,:,:,is), &
fnl(ik,is)%r(:,:,:), eigr, ps )
ELSE
CALL dforce( ik, nb, c0(:,:,:,is), cdesc, fi(:,:,is), c2, gv, vpot(:,:,:,is), &
CALL dforce( ik, nb, c0(:,:,:,is), cdesc, fi(:,:,is), c2, vpot(:,:,:,is), &
fnl(ik,is)%c, eigr, ps )
END IF
IF( tlam ) THEN
@ -364,7 +364,7 @@
CALL wave_verlet( cp(:,nb,ik,is), c0(:,nb,ik,is), svar1, svar2, svar3, c2 )
END IF
IF( .NOT. cdesc%gamma ) THEN
cp(:,nb,ik,is) = cp(:,nb,ik,is) * gv%kg_mask_l(:,ik)
cp(:,nb,ik,is) = cp(:,nb,ik,is) * gkmask_l(:,ik)
ELSE
IF( cdesc%gzero ) cp(1,nb,ik,is) = REAL( cp(1,nb,ik,is), dbl )
END IF
@ -386,11 +386,10 @@
!cdesc is the desciptor for the wf
!gv g vector
!eigr==e^ig*r f is the occupation number
!fnl if the factor non local
SUBROUTINE runcp_force_pairing(ttprint, tortho, tsde, cm, c0, cp, cdesc, gv, kp, ps, &
SUBROUTINE runcp_force_pairing(ttprint, tortho, tsde, cm, c0, cp, cdesc, kp, ps, &
vpot, eigr, fi, ekinc, timerd, timeorto, ht, ei, fnl, vnosee)
! same as runcp, except that electrons are paired forcedly
@ -409,7 +408,7 @@
USE wave_functions, ONLY : rande, cp_kinetic_energy, gram
USE wave_base, ONLY: frice, wave_steepest, wave_verlet
USE wave_base, ONLY: hpsi
USE cp_types, ONLY: recvecs, pseudo
USE cp_types, ONLY: pseudo
USE cell_module, ONLY: boxdimensions
USE time_step, ONLY: delt
USE forces, ONLY: dforce
@ -422,6 +421,7 @@
USE constants, ONLY: au
USE io_global, ONLY: ionode
USE wave_constrains, ONLY: update_lambda
USE reciprocal_space_mesh, ONLY: gkmask_l
IMPLICIT NONE
@ -432,7 +432,6 @@
TYPE (wave_descriptor), INTENT(IN) :: cdesc
TYPE (pseudo), INTENT(IN) :: ps
COMPLEX(dbl) :: eigr(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
REAL(dbl), INTENT(INOUT) :: fi(:,:,:)
TYPE (boxdimensions), INTENT(IN) :: ht
@ -573,9 +572,9 @@
! dforce calcola la forza c2 e c3 sulle bande i e i+1 (sono reali => ne fa due alla volta)
! per il vpot (da potential ed e' il potetnziale di KS) in spin up e in down
!
CALL dforce( ik, i, c0(:,:,:,1), cdesc, fi(:,:,1), c2, c3, gv, vpot(:,:,:,1), &
CALL dforce( ik, i, c0(:,:,:,1), cdesc, fi(:,:,1), c2, c3, vpot(:,:,:,1), &
fnl(ik, 1)%r(:,:,:), eigr, ps )
CALL dforce( ik, i, c0(:,:,:,1), cdesc, fi(:,:,2), c4, c5, gv, vpot(:,:,:,2), &
CALL dforce( ik, i, c0(:,:,:,1), cdesc, fi(:,:,2), c4, c5, vpot(:,:,:,2), &
fnl(ik, 2)%r(:,:,:), eigr, ps )
!
! accoppia c2 e c3 da vpot (spin 1) stato i e i+1
@ -598,16 +597,16 @@
! raddoppia (questa pesa un 30/40% sul conto)
! qui sono in C => ogni FFT serve per una wf complessa che ha bisogno di due componenti
CALL dforce( ik, i, c0(:,:,:,1), cdesc, fi(:,:,1), c2, gv, vpot(:,:,:,1), &
CALL dforce( ik, i, c0(:,:,:,1), cdesc, fi(:,:,1), c2, vpot(:,:,:,1), &
fnl(ik, 1)%c, eigr, ps )
CALL dforce( ik, i, c0(:,:,:,1), cdesc, fi(:,:,2), c4, gv, vpot(:,:,:,2), &
CALL dforce( ik, i, c0(:,:,:,1), cdesc, fi(:,:,2), c4, vpot(:,:,:,2), &
fnl(ik, 2)%c, eigr, ps )
c2 = occup(i, ik)*c2 + occdown(i, ik)*c4
CALL dforce( ik, i+1, c0(:,:,:,1), cdesc, fi(:,:,1), c3, gv, vpot(:,:,:,1), &
CALL dforce( ik, i+1, c0(:,:,:,1), cdesc, fi(:,:,1), c3, vpot(:,:,:,1), &
fnl(ik, 1)%c, eigr, ps )
CALL dforce( ik, i+1, c0(:,:,:,1), cdesc, fi(:,:,2), c5, gv, vpot(:,:,:,2), &
CALL dforce( ik, i+1, c0(:,:,:,1), cdesc, fi(:,:,2), c5, vpot(:,:,:,2), &
fnl(ik, 2)%c, eigr, ps )
c3 = occup(i+1, ik)*c3 + occdown(i+1, ik)*c5
@ -660,8 +659,8 @@
CALL wave_verlet( cp(:,i+1,ik,1), c0(:,i+1,ik,1), svar1, svar2, svar3, c3 )
END IF
IF( .NOT. cdesc%gamma ) THEN
cp(:,i,ik,1) = cp(:,i,ik,1) * gv%kg_mask_l(:,ik)
cp(:,i+1,ik,1) = cp(:,i+1,ik,1) * gv%kg_mask_l(:,ik)
cp(:,i,ik,1) = cp(:,i,ik,1) * gkmask_l(:,ik)
cp(:,i+1,ik,1) = cp(:,i+1,ik,1) * gkmask_l(:,ik)
ELSE
IF( cdesc%gzero ) cp(1,i,ik,1) = REAL( cp(1,i,ik,1), dbl )
IF( cdesc%gzero ) cp(1,i+1,ik,1) = REAL( cp(1,i+1,ik,1), dbl )
@ -680,15 +679,15 @@
! per questo conservo in ei_t(:,:,2) il suo autovalore
!
IF ( cdesc%gamma ) THEN
CALL dforce( ik, nb, c0(:,:,:,1), cdesc, fi(:,:,1), c2, gv, vpot(:,:,:,1), &
CALL dforce( ik, nb, c0(:,:,:,1), cdesc, fi(:,:,1), c2, vpot(:,:,:,1), &
fnl(ik,1)%r(:,:,:), eigr, ps )
CALL dforce( ik, nb, c0(:,:,:,1), cdesc, fi(:,:,2), c3, gv, vpot(:,:,:,2), &
CALL dforce( ik, nb, c0(:,:,:,1), cdesc, fi(:,:,2), c3, vpot(:,:,:,2), &
fnl(ik,2)%r(:,:,:), eigr, ps )
c2 = occup(nb, ik)*c2 + occdown(nb, ik)*c3
ELSE
CALL dforce( ik, nb, c0(:,:,:,1), cdesc, fi(:,:,1), c2, gv, vpot(:,:,:,1), &
CALL dforce( ik, nb, c0(:,:,:,1), cdesc, fi(:,:,1), c2, vpot(:,:,:,1), &
fnl(ik,1)%c, eigr, ps )
CALL dforce( ik, nb, c0(:,:,:,1), cdesc, fi(:,:,2), c3, gv, vpot(:,:,:,2), &
CALL dforce( ik, nb, c0(:,:,:,1), cdesc, fi(:,:,2), c3, vpot(:,:,:,2), &
fnl(ik,2)%c, eigr, ps )
c2 = occup(nb, ik)*c2 + occdown(nb, ik)*c3
END IF
@ -716,7 +715,7 @@
CALL wave_verlet( cp(:,nb,ik,1), c0(:,nb,ik,1), svar1, svar2, svar3, c2 )
END IF
IF( .NOT. cdesc%gamma ) THEN
cp(:,nb,ik,1) = cp(:,nb,ik,1) * gv%kg_mask_l(:,ik)
cp(:,nb,ik,1) = cp(:,nb,ik,1) * gkmask_l(:,ik)
ELSE
IF( cdesc%gzero ) cp(1,nb,ik,1) = REAL( cp(1,nb,ik,1), dbl )
END IF
@ -804,9 +803,9 @@
! Compute fictitious kinetic energy of the electrons at time t
ekinc(1) = cp_kinetic_energy( 1, cp(:,:,:,1), cm(:,:,:,1), cdesc, kp, &
gv%kg_mask_l, pmss, delt)
gkmask_l, pmss, delt)
ekinc(2) = cp_kinetic_energy( 2, cp(:,:,:,1), cm(:,:,:,1), cdesc, kp, &
gv%kg_mask_l, pmss, delt)
gkmask_l, pmss, delt)
DEALLOCATE( ei_t, svar3, c2, c3, c4, c5, cgam, gam, occup, occdown, STAT=ierr)

72
CPV/rundiis.f90
View File

@ -28,7 +28,7 @@
! ----------------------------------------------
! BEGIN manual
SUBROUTINE rundiis(tprint, rhoe, desc, atoms, gv, kp, &
SUBROUTINE rundiis(tprint, rhoe, desc, atoms, kp, &
ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cgrad, cdesc, tcel, ht0, fi, eig, &
fnl, vpot, doions, edft )
@ -82,7 +82,7 @@
USE electrons_module, ONLY: pmss
USE time_step, ONLY: delt
USE wave_functions, ONLY: gram, proj, crot
USE phase_factors_module, ONLY: strucf
USE phase_factors_module, ONLY: strucf, phfacs
USE charge_mix
USE charge_density, ONLY: rhoofr
USE guess
@ -99,6 +99,9 @@
USE atoms_type_module, ONLY: atoms_type
USE charge_types, ONLY: charge_descriptor
USE control_flags, ONLY: force_pairing
use grid_dimensions, only: nr1, nr2, nr3
USE reciprocal_vectors, ONLY: mill_l
USE gvecp, ONLY: ngm
IMPLICIT NONE
@ -115,7 +118,6 @@
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
TYPE (charge_descriptor) :: desc
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
TYPE (boxdimensions), INTENT(INOUT) :: ht0
REAL(dbl) :: fi(:,:,:)
@ -193,12 +195,14 @@
! ... starting guess on the wavefunctions
! CALL guessrho(rhoe, cm, c0, fi, gv, kp, ht0)
! CALL guessrho(rhoe, cm, c0, fi, kp, ht0)
CALL strucf(sfac, atoms, eigr, ei1, ei2, ei3, gv)
CALL rhoofr(gv, kp, c0, cdesc, fi, rhoe, desc, ht0)
CALL newrho(rhoe(:,:,:,1), drho, gv, 0) ! memorize density
CALL strucf(sfac, atoms, eigr, ei1, ei2, ei3, gv)
CALL phfacs( ei1, ei2, ei3, eigr, mill_l, atoms%taus, nr1, nr2, nr3, atoms%nat )
CALL strucf( sfac, ei1, ei2, ei3, mill_l, ngm )
CALL rhoofr( kp, c0, cdesc, fi, rhoe, desc, ht0)
CALL newrho(rhoe(:,:,:,1), drho, 0) ! memorize density
CALL phfacs( ei1, ei2, ei3, eigr, mill_l, atoms%taus, nr1, nr2, nr3, atoms%nat )
CALL strucf( sfac, ei1, ei2, ei3, mill_l, ngm )
CALL guessc0( .NOT. kp%gamma_only, c0, cm, cdesc)
! ... Initialize the rotation index srot
@ -233,13 +237,13 @@
END IF
! ... self consistent energy
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fi, gv, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL rhoofr(gv, kp, c0, cdesc, fi, rhoe, desc, ht0)
CALL vofrhos(.FALSE., rhoe, desc, tforce, tstress, tforce, atoms, gv, &
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fi, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL rhoofr( kp, c0, cdesc, fi, rhoe, desc, ht0)
CALL vofrhos(.FALSE., rhoe, desc, tforce, tstress, tforce, atoms, &
kp, fnl, vpot, ps, c0, cdesc, fi, eigr, ei1, ei2, ei3, sfac, timepre, ht0, edft)
! ... density upgrade
CALL newrho(rhoe(:,:,:,1), drho, gv, idiis)
CALL newrho(rhoe(:,:,:,1), drho, idiis)
IF (ionode) WRITE( stdout,45) idiis, edft%etot, drho
dene = abs(edft%etot - etot_m)
etot_m = edft%etot
@ -247,8 +251,8 @@
45 FORMAT('etot drho ',i3,1x,2(1x,f18.10))
! ... recalculate potential
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fi, gv, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL vofrhos(.FALSE., rhoe, desc, tforce, tstress, tforce, atoms, gv, kp, fnl, &
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fi, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL vofrhos(.FALSE., rhoe, desc, tforce, tstress, tforce, atoms, kp, fnl, &
vpot, ps, c0, cdesc, fi, eigr, ei1, ei2, ei3, sfac, timepre, ht0, edft)
IF( idiis /= 1 )THEN
@ -265,9 +269,9 @@
! ... calculate lambda_i,j=<c_i| H |c_j>
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, gv, kp, fnl, ps%wsg, ps%wnl, eigr)
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL dforce_all( 1, c0(:,:,:,1), cdesc, fi(:,:,1), cgrad(:,:,:,1), gv, vpot(:,:,:,1), &
CALL dforce_all( 1, c0(:,:,:,1), cdesc, fi(:,:,1), cgrad(:,:,:,1), vpot(:,:,:,1), &
fnl(:,1), eigr, ps)
IF(.NOT.kp%gamma_only) THEN
@ -283,8 +287,8 @@
call adjef_s(eig(1,1,1),fi(1,1,1),efermi,nel, cdesc%nbl( 1 ),temp_elec,sume)
call enthropy_s(fi(1,1,1),temp_elec,cdesc%nbl(1),edft%ent)
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, gv, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL dforce_all( 1, c0(:,:,:,1), cdesc, fi(:,:,1), cgrad(:,:,:,1), gv, &
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL dforce_all( 1, c0(:,:,:,1), cdesc, fi(:,:,1), cgrad(:,:,:,1), &
vpot(:,:,:,1), fnl(:,1), eigr, ps)
DO ik = 1, kp%nkpt
@ -296,9 +300,9 @@
ELSE
! ... DIIS on c0 at FIXED potential
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, gv, kp, fnl, ps%wsg, ps%wnl, eigr)
edft%enl = nlrh_m(c0, cdesc, tforce, atoms, fs, kp, fnl, ps%wsg, ps%wnl, eigr)
CALL dforce_all( 1, c0(:,:,:,1), cdesc, fi(:,:,1), cgrad(:,:,:,1), gv, &
CALL dforce_all( 1, c0(:,:,:,1), cdesc, fi(:,:,1), cgrad(:,:,:,1), &
vpot(:,:,:,1), fnl(:,1), eigr, ps)
IF( kp%gamma_only ) THEN
@ -333,11 +337,11 @@
svar2 = -1.d0
svar3_0 = delt * delt / pmss(1)
IF(.NOT.kp%gamma_only) THEN
CALL simupd_kp(ekinc,doions,gv,kp,c0(:,:,:,1),cgrad(:,:,:,1),cdesc, svar1,svar2, &
CALL simupd_kp(ekinc,doions,kp,c0(:,:,:,1),cgrad(:,:,:,1),cdesc, svar1,svar2, &
svar3_0,edft%etot,fs(:,:,1),eigr,sfac,ps%vps,ps%wsg,ps%wnl,treset_diis, &
istate,cnorm,eold,ndiis,nowv)
ELSE
CALL simupd(ekinc,doions,gv,c0(:,:,:,1),cgrad(:,:,:,1),cdesc, svar1,svar2, &
CALL simupd(ekinc,doions,c0(:,:,:,1),cgrad(:,:,:,1),cdesc, svar1,svar2, &
svar3_0,edft%etot,fs(:,1,1),eigr,sfac,ps%vps,ps%wsg,ps%wnl, &
treset_diis,istate,cnorm,eold,ndiis,nowv)
END IF
@ -395,7 +399,7 @@
! ----------------------------------------------
! BEGIN manual
SUBROUTINE runsdiis(tprint, rhoe, desc, atoms, gv, kp, &
SUBROUTINE runsdiis(tprint, rhoe, desc, atoms, kp, &
ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cgrad, cdesc, tcel, ht0, fi, eig, &
fnl, vpot, doions, edft )
@ -447,7 +451,6 @@
USE electrons_module, ONLY: ei, pmss
USE time_step, ONLY: delt
USE wave_functions, ONLY: gram, proj, update_wave_functions
USE phase_factors_module, ONLY: strucf
USE cp_types
USE diis
USE cell_module, ONLY: boxdimensions
@ -478,7 +481,6 @@
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
COMPLEX(dbl) :: sfac(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
TYPE (boxdimensions), INTENT(INOUT) :: ht0
REAL(dbl) :: fi(:,:,:)
@ -546,7 +548,7 @@
END IF
CALL kspotential( .FALSE., tforce, tstress, rhoe, desc, &
atoms, gv, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0, fi, fnl, vpot, edft, timepre )
atoms, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0, fi, fnl, vpot, edft, timepre )
s0 = cclock()
seconds_per_iter = (s0 - old_clock_value)
@ -562,7 +564,7 @@
isteep = isteep + 1
cm = c0
CALL runcp(.FALSE., tortho, tsde, cm, c0, &
cgrad, cdesc, gv, kp, ps, vpot, eigr, fi, &
cgrad, cdesc, kp, ps, vpot, eigr, fi, &
ekinc, timerd, timeorto, ht0, ei, fnl, 0.0d0 )
CALL update_wave_functions(cm, c0, cgrad, cdesc)
@ -591,7 +593,7 @@
! ... so on).
CALL dforce_all( ispin, c0(:,:,:,ispin), cdesc, fi(:,:,ispin), cgrad(:,:,:,ispin), &
gv, vpot(:,:,:,ispin), fnl(:,ispin), eigr, ps)
vpot(:,:,:,ispin), fnl(:,ispin), eigr, ps)
IF(.NOT.kp%gamma_only) THEN
DO ik = 1, kp%nkpt
@ -603,12 +605,12 @@
s4 = cclock()
IF(.NOT.kp%gamma_only) THEN
CALL simupd_kp(ekinc(ispin), ddoions(ispin), gv, kp, c0(:,:,:,ispin), &
CALL simupd_kp(ekinc(ispin), ddoions(ispin), kp, c0(:,:,:,ispin), &
cgrad(:,:,:,ispin), cdesc, svar1, svar2, svar3_0, edft%etot, fi(:,:,ispin), &
eigr, sfac, ps%vps, ps%wsg, ps%wnl, ttreset_diis(ispin), istate, &
cnorm, eold, ndiis, nowv)
ELSE
CALL simupd(ekinc(ispin), ddoions(ispin), gv, c0(:,:,:,ispin), cgrad(:,:,:,ispin), cdesc, &
CALL simupd(ekinc(ispin), ddoions(ispin), c0(:,:,:,ispin), cgrad(:,:,:,ispin), cdesc, &
svar1, svar2, svar3_0, edft%etot, fi(:,1,ispin), eigr, sfac, &
ps%vps, ps%wsg, ps%wnl, ttreset_diis(ispin), istate, cnorm, &
eold, ndiis, nowv)
@ -645,7 +647,7 @@
END IF
IF ( tprint ) THEN
CALL diis_eigs(.TRUE., atoms, c0, cdesc, fi, vpot, cgrad, fnl, ps, eigr, gv, kp)
CALL diis_eigs(.TRUE., atoms, c0, cdesc, fi, vpot, cgrad, fnl, ps, eigr, kp)
END IF
cgrad = c0
@ -694,7 +696,7 @@
! ----------------------------------------------
SUBROUTINE diis_eigs(tortho, atoms, c, cdesc, fi, vpot, eforce, fnle, ps, eigr, gv, kp)
SUBROUTINE diis_eigs(tortho, atoms, c, cdesc, fi, vpot, eforce, fnle, ps, eigr, kp)
USE kinds
USE cp_types
USE wave_types
@ -713,6 +715,7 @@
USE mp, ONLY: mp_sum
USE mp_global, ONLY: mpime, nproc, group
USE atoms_type_module, ONLY: atoms_type
USE reciprocal_space_mesh, ONLY: gkx_l, gk_l
IMPLICIT NONE
@ -720,7 +723,6 @@
COMPLEX(dbl), INTENT(inout) :: c(:,:,:,:)
COMPLEX(dbl), INTENT(inout) :: eforce(:,:,:,:)
TYPE (wave_descriptor), INTENT(in) :: cdesc
TYPE (recvecs), INTENT(in) :: gv
REAL (dbl), INTENT(in) :: vpot(:,:,:,:), fi(:,:,:)
TYPE (kpoints), INTENT(in) :: kp
TYPE (pseudo), INTENT(INOUT) :: ps
@ -755,11 +757,11 @@
DO ispin = 1, nspin
DO ik = 1, SIZE( c, 3 )
CALL nlsm1( ispin, ps%wnl(:,:,:,ik), atoms, eigr, c(:,:,ik,ispin), cdesc, &
gv%khg_l(:,ik), gv%kgx_l(:,:,ik), fnle(ik,ispin))
gk_l(:,ik), gkx_l(:,:,ik), fnle(ik,ispin))
END DO
! ... Calculate | dH / dpsi(j) >
CALL dforce_all( ispin, c(:,:,:,ispin), cdesc, fi(:,:,ispin), eforce(:,:,:,ispin), gv, &
CALL dforce_all( ispin, c(:,:,:,ispin), cdesc, fi(:,:,ispin), eforce(:,:,:,ispin), &
vpot(:,:,:,ispin), fnle(:,ispin), eigr, ps)
DO ik = 1, kp%nkpt

47
CPV/runsd.f90
View File

@ -29,7 +29,7 @@
! -----------------------------------------------------------------------
! BEGIN manual
SUBROUTINE runsd(tortho, tprint, tforce, rhoe, desc, atoms_0, gv, kp, &
SUBROUTINE runsd(tortho, tprint, tforce, rhoe, desc, atoms_0, kp, &
ps, eigr, ei1, ei2, ei3, sfac, c0, cm, cp, cdesc, tcel, ht0, occ, ei, &
fnl, vpot, doions, edft, maxnstep, sdthr )
@ -37,22 +37,25 @@
! END manual
! ... declare modules
USE energies, ONLY: dft_energy_type, print_energies
USE wave_functions, ONLY: update_wave_functions
USE check_stop, ONLY: check_stop_now
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE cell_module, ONLY: boxdimensions
USE brillouin, ONLY: kpoints
USE cp_types, ONLY: recvecs, pseudo
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
USE potentials, ONLY: kspotential
USE atoms_type_module, ONLY: atoms_type
USE runcp_module, ONLY: runcp
USE phase_factors_module, ONLY : strucf
USE charge_types, ONLY: charge_descriptor
USE control_flags, ONLY: force_pairing
USE energies, ONLY: dft_energy_type, print_energies
USE wave_functions, ONLY: update_wave_functions
USE check_stop, ONLY: check_stop_now
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE cell_module, ONLY: boxdimensions
USE brillouin, ONLY: kpoints
USE cp_types, ONLY: pseudo
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
USE potentials, ONLY: kspotential
USE atoms_type_module, ONLY: atoms_type
USE runcp_module, ONLY: runcp
USE phase_factors_module, ONLY: strucf, phfacs
USE charge_types, ONLY: charge_descriptor
USE control_flags, ONLY: force_pairing
use grid_dimensions, only: nr1, nr2, nr3
USE reciprocal_vectors, ONLY: mill_l
USE gvecp, ONLY: ngm
IMPLICIT NONE
@ -70,7 +73,6 @@
COMPLEX(dbl) :: ei1(:,:)
COMPLEX(dbl) :: ei2(:,:)
COMPLEX(dbl) :: ei3(:,:)
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
TYPE (boxdimensions), INTENT(INOUT) :: ht0
REAL(dbl) :: occ(:,:,:)
@ -123,19 +125,20 @@
old_clock_value = cclock()
CALL strucf(sfac, atoms_0, eigr, ei1, ei2, ei3, gv)
CALL phfacs( ei1, ei2, ei3, eigr, mill_l, atoms_0%taus, nr1, nr2, nr3, atoms_0%nat )
CALL strucf( sfac, ei1, ei2, ei3, mill_l, ngm )
STEEPEST_DESCENT: DO iter = 1, maxnstep
s1 = cclock()
CALL kspotential( ttprint, ttforce, ttstress, rhoe, desc, &
atoms_0, gv, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0, occ, fnl, vpot, edft, timepre )
atoms_0, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0, occ, fnl, vpot, edft, timepre )
s2 = cclock()
CALL runcp(ttprint, ttortho, ttsde, cm, c0, cp, &
cdesc, gv, kp, ps, vpot, eigr, occ, ekincs, timerd, &
cdesc, kp, ps, vpot, eigr, occ, ekincs, timerd, &
timeorto, ht0, ei, fnl, vnosee)
ekinc = SUM( ekincs )
@ -174,7 +177,7 @@
IF( tforce ) THEN
atoms_0%for = 0.0d0
CALL kspotential( ttprint, tforce, ttstress, rhoe, desc, &
atoms_0, gv, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0, occ, fnl, vpot, edft, timepre )
atoms_0, kp, ps, eigr, ei1, ei2, ei3, sfac, c0, cdesc, tcel, ht0, occ, fnl, vpot, edft, timepre )
IF(ionode ) THEN
WRITE( stdout,fmt="(12X,'runsd: fion and edft calculated = ',F14.6)") edft%etot
END IF

9
CPV/smcp.f90
View File

@ -73,7 +73,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
use gvecs, only: ngs
use gvecb, only: ngb
use gvecw, only: ngw
use reciprocal_vectors, only: ng0 => gstart
use reciprocal_vectors, only: ng0 => gstart, mill_l
use ions_base, only: na, nat, pmass, nax, nsp, rcmax
use ions_base, only: ind_srt, ions_vel, ions_cofmass, ions_kinene, ions_temp
use ions_base, only: ions_thermal_stress, ions_vrescal, fricp, greasp, iforce
@ -146,6 +146,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
USE from_restart_module, ONLY: from_restart
USE runcp_module, ONLY: runcp_uspp
use phase_factors_module, only: strucf
USE cp_main_variables, ONLY: ei1, ei2, ei3, eigr, sfac, irb, taub, eigrb, &
rhog, rhor, rhos, rhoc, becdr, bephi, becp, ema0bg, allocate_mainvar
@ -616,7 +617,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
!
! strucf calculates the structure factor sfac
!
call strucf(ei1,ei2,ei3,sfac)
call strucf( sfac, ei1, ei2, ei3, mill_l, ngs )
IF(sm_k == 1) call formf(tfirst,eself)
@ -944,7 +945,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
!
! strucf calculates the structure factor sfac
!
call strucf(ei1,ei2,ei3,sfac)
call strucf( sfac, ei1, ei2, ei3, mill_l, ngs )
if (thdyn) call formf(tfirst,eself)
!
IF(sm_k==1) THEN
@ -1228,7 +1229,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
if( thdyn ) then
hold = h
h = hnew
call newinit(ibrav)
call newinit( h )
call newnlinit
else
hold = h

2
CPV/smd.f90
View File

@ -129,7 +129,7 @@ subroutine sminit (ibrav,celldm, ecut, ecutw,ndr,nbeg, &
! ==== generate true g-space ====
! ==============================================================
!
call newinit( ibrav )
call newinit( h )
!
344 format(' ibrav = ',i4,' cell parameters ',/)
345 format(3(4x,f10.5))

189
CPV/stress.f90
View File

@ -58,25 +58,27 @@
! BEGIN manual
SUBROUTINE pstress( strvxc, rhoeg, vxc, pail, desr, &
gv, fnl, ps, c0, cdesc, occ, eigr, sfac, grho, v2xc, box, edft)
fnl, ps, c0, cdesc, occ, eigr, sfac, grho, v2xc, box, edft)
! this routine computes stress tensor from dft total energy
! ----------------------------------------------
! END manual
! ... declare modules
USE cp_types, ONLY: recvecs, pseudo
USE cell_module, ONLY: boxdimensions
USE energies, ONLY: dft_energy_type
USE ions_base, ONLY: nsp
USE mp_global, ONLY: mpime, nproc, group
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 cp_types, ONLY: pseudo
USE cell_module, ONLY: boxdimensions
USE energies, ONLY: dft_energy_type
USE ions_base, ONLY: nsp
USE mp_global, ONLY: mpime, nproc, group
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
USE control_flags, ONLY: iprsta
USE control_flags, ONLY: iprsta
USE reciprocal_vectors, ONLY: gx
USE gvecp, ONLY: ngm
IMPLICIT NONE
@ -87,7 +89,6 @@
COMPLEX(dbl), INTENT(IN) :: sfac(:,:)
REAL(dbl), INTENT(IN) :: occ(:,:,:)
TYPE (pseudo), INTENT(IN) :: ps
TYPE (recvecs), INTENT(IN) :: gv
COMPLEX(dbl), INTENT(IN) :: c0(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
TYPE (boxdimensions), INTENT(IN) :: box
@ -118,10 +119,10 @@
! ... compute G_alpha * G_beta
ALLOCATE(gagx_l(6,gv%ng_l))
ALLOCATE(gagx_l(6,ngm))
DO k = 1, 6
DO ig = 1, gv%ng_l
gagx_l(k,ig) = gv%gx_l(alpha(k),ig) * gv%gx_l(beta(k),ig) * tpiba2
DO ig = 1, ngm
gagx_l(k,ig) = gx(alpha(k),ig) * gx(beta(k),ig) * tpiba2
END DO
END DO
@ -129,17 +130,17 @@
! ... compute kinetic energy contribution
CALL stress_kin(dekin, c0, cdesc, occ, gagx_l, gv)
CALL stress_kin(dekin, c0, cdesc, occ, gagx_l)
IF(timing) s2 = cclock()
! ... compute hartree energy contribution
CALL stress_har(deht, ehr, sfac, ps, rhoeg, gagx_l, gv, box)
CALL stress_har(deht, ehr, sfac, ps, rhoeg, gagx_l, box)
IF(timing) s3 = cclock()
! ... compute exchange & correlation energy contribution
CALL stress_xc(dexc, strvxc, sfac, vxc, grho, v2xc, gagx_l, gv, &
CALL stress_xc(dexc, strvxc, sfac, vxc, grho, v2xc, gagx_l, &
ps%tnlcc, ps%rhocp, box)
IF(timing) s4 = cclock()
@ -151,12 +152,12 @@
IF(timing) s5 = cclock()
CALL pseudo_stress(deps, edft%epseu, gv, gagx_l, sfac, ps%dvps, rhoeg, box)
CALL pseudo_stress(deps, edft%epseu, gagx_l, sfac, ps%dvps, rhoeg, box)
IF(timing) s6 = cclock()
! ... compute enl (non-local) contribution
CALL stress_nl(denl, gagx_l, gv, c0, cdesc, occ, eigr, ps%wsg,fnl, &
CALL stress_nl(denl, gagx_l, c0, cdesc, occ, eigr, ps%wsg,fnl, &
ps%wnl(:,:,:,1), edft%enl)
IF(timing) s7 = cclock()
@ -231,7 +232,7 @@
! BEGIN manual
SUBROUTINE stress_nl(denl, gagx_l, gv, c0, cdesc, occ, eigr, wsg, fnl, wnl, enl)
SUBROUTINE stress_nl(denl, gagx_l, c0, cdesc, occ, eigr, wsg, fnl, wnl, enl)
! this routine computes nl part of the stress tensor from dft total energy
! ----------------------------------------------
@ -247,13 +248,13 @@
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
USE cell_base, ONLY: tpiba2
USE cp_types, ONLY: recvecs
USE control_flags, ONLY: force_pairing
USE reciprocal_vectors, ONLY: gstart, gzero, g
USE reciprocal_space_mesh, ONLY: gkx_l, gk_l
IMPLICIT NONE
! ... declare subroutine arguments
TYPE (recvecs), INTENT(IN) :: gv
REAL(dbl), INTENT(IN) :: occ(:,:,:)
COMPLEX(dbl), INTENT(IN) :: c0(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
@ -297,17 +298,17 @@
me = mpime + 1
nspin = cdesc%nspin
IF(gv%gzero) THEN
IF(gzero) THEN
denl = - enl * dalbe
ELSE
denl = 0.0_dbl
END IF
ngw = gv%ngw_l
ngw = cdesc%ngwl
! ... initialize array wnla
ALLOCATE(wnla(ngw, lnlx, nsp))
CALL nlin_stress(wnla,gv)
CALL nlin_stress(wnla)
ALLOCATE(gwork(ngw))
ALLOCATE(gwtmp(ngw))
@ -328,37 +329,37 @@
IF (ts) THEN
igh=igh+1
gspha(1) = 0.0d0
CALL spharm(gspha, gv%kgx_l(:,:,1), gv%khg_l(:,1), ngw, 0, 0)
DO ig = gv%gstart, ngw
gspha(ig) = gspha(ig) / (gv%hg_l(ig)*tpiba2)
CALL spharm(gspha, gkx_l(:,:,1), gk_l(:,1), ngw, 0, 0)
DO ig = gstart, ngw
gspha(ig) = gspha(ig) / (g(ig)*tpiba2)
END DO
DO kk = 1, 6
fnls = 0.0d0
iss = 1
gwork(1) = 0.0d0
DO ig = gv%gstart, ngw
DO ig = gstart, ngw
gwork(ig) = gagx_l(kk,ig) * gspha(ig)
END DO
DO is = 1, nspnl
ll = l2ind(1,is)
IF(ll.GT.0) THEN
ALLOCATE(auxc(gv%ngw_l,na(is)))
ALLOCATE(auxc(ngw,na(is)))
gwtmp(1) = 0.0d0
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
gwtmp(ig) = gwork(ig) * (wnl(ig,ll,is)-wnla(ig,ll,is))
END DO
DO ia = 1, na(is)
auxc(1,ia) = CMPLX(0.0d0,0.0d0)
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
auxc(ig,ia) = csign(0) * gwtmp(ig) * eigr(ig,ia+iss-1)
END DO
END DO
CALL DGEMM( 'T', 'N', na(is), nx, 2*gv%ngw_l, 1.0d0, auxc(1,1), &
2*gv%ngw_l, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, 0.0d0, fnls(iss,1), nsanl )
CALL DGEMM( 'T', 'N', na(is), nx, 2*ngw, 1.0d0, auxc(1,1), &
2*ngw, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, 0.0d0, fnls(iss,1), nsanl )
DEALLOCATE(auxc)
@ -386,9 +387,9 @@
igh=igh+1
gspha(1) = 0.0d0
CALL spharm(gspha, gv%kgx_l(:,:,1), gv%khg_l(:,1), ngw, 1, m-2)
DO ig = gv%gstart, ngw
gspha(ig) = gspha(ig) / (gv%hg_l(ig)*tpiba2)
CALL spharm(gspha, gkx_l(:,:,1), gk_l(:,1), ngw, 1, m-2)
DO ig = gstart, ngw
gspha(ig) = gspha(ig) / (g(ig)*tpiba2)
END DO
DO kk=1,6
@ -396,7 +397,7 @@
fnls = 0.0d0
iss = 1
gwork(1) = 0.0d0
DO ig=gv%gstart,gv%ngw_l
DO ig = gstart, ngw
gwork(ig)= gagx_l(kk,ig) * gspha(ig)
END DO
@ -404,22 +405,22 @@
ll = l2ind(2,is)
IF(ll.GT.0) THEN
ALLOCATE(auxc(gv%ngw_l,na(is)))
ALLOCATE(auxc(ngw,na(is)))
gwtmp(1) = 0.0d0
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
gwtmp(ig) = gwork(ig) * ( 3.d0 * wnl(ig,ll,is) - wnla(ig,ll,is) )
END DO
DO ia = 1, na(is)
auxc(1,ia) = CMPLX(0.0d0,0.0d0)
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
auxc(ig,ia) = csign(1) * gwtmp(ig) * eigr(ig,ia+iss-1)
END DO
END DO
CALL DGEMM( 'T', 'N', na(is), nx, 2*gv%ngw_l, 1.0d0, &
auxc(1,1),2*gv%ngw_l, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, &
CALL DGEMM( 'T', 'N', na(is), nx, 2*ngw, 1.0d0, &
auxc(1,1),2*ngw, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, &
0.0d0, fnls(iss,1), nsanl )
DEALLOCATE(auxc)
@ -461,9 +462,9 @@
igh=igh+1
gspha(1) = 0.0d0
CALL spharm(gspha, gv%kgx_l(:,:,1), gv%khg_l(:,1), ngw, 2, m-3)
DO ig = gv%gstart, ngw
gspha(ig) = gspha(ig) / (gv%hg_l(ig)*tpiba2)
CALL spharm(gspha, gkx_l(:,:,1), gk_l(:,1), ngw, 2, m-3)
DO ig = gstart, ngw
gspha(ig) = gspha(ig) / (g(ig)*tpiba2)
END DO
DO kk=1,6
@ -471,7 +472,7 @@
fnls = 0.0d0
iss = 1
gwork(1) = 0.0d0
DO ig=gv%gstart,gv%ngw_l
DO ig=gstart,ngw
gwork(ig)= gagx_l(kk,ig) * gspha(ig)
END DO
@ -479,26 +480,26 @@
ll = l2ind(3,is)
IF(ll.GT.0) THEN
ALLOCATE(auxc(gv%ngw_l,na(is)))
ALLOCATE(auxc(ngw,na(is)))
gwtmp(1) = 0.0d0
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
gwtmp(ig) = gwork(ig) * ( 5.d0 * wnl(ig,ll,is) - wnla(ig,ll,is) )
END DO
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
gwtmp(ig) = gwtmp(ig) - 2.0d0/3.0d0 * dm(kk,m) * wnl(ig,ll,is)
END DO
DO ia= 1 , na(is)
auxc(1,ia) = CMPLX(0.0d0,0.0d0)
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
auxc(ig,ia) = csign(2) * gwtmp(ig) * eigr(ig,ia+iss-1)
END DO
END DO
CALL DGEMM( 'T', 'N', na(is), nx, 2*gv%ngw_l, 1.0d0, &
auxc(1,1), 2*gv%ngw_l, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, &
CALL DGEMM( 'T', 'N', na(is), nx, 2*ngw, 1.0d0, &
auxc(1,1), 2*ngw, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, &
0.0d0, fnls(iss,1), nsanl )
DEALLOCATE(auxc)
@ -545,9 +546,9 @@
igh=igh+1
gspha(1) = 0.0d0
CALL spharm(gspha, gv%kgx_l(:,:,1), gv%khg_l(:,1), ngw, 3, m-4)
DO ig = gv%gstart, ngw
gspha(ig) = gspha(ig) / (gv%hg_l(ig)*tpiba2)
CALL spharm(gspha, gkx_l(:,:,1), gk_l(:,1), ngw, 3, m-4)
DO ig = gstart, ngw
gspha(ig) = gspha(ig) / (g(ig)*tpiba2)
END DO
DO kk=1,6
@ -555,7 +556,7 @@
fnls = 0.0d0
iss = 1
gwork(1) = 0.0d0
DO ig=gv%gstart,gv%ngw_l
DO ig = gstart, ngw
gwork(ig)= gagx_l(kk,ig) * gspha(ig)
END DO
@ -564,36 +565,36 @@
ll = l2ind(4,is)
IF(ll > 0) THEN
ALLOCATE(auxc(gv%ngw_l,na(is)))
ALLOCATE(auxc(ngw,na(is)))
gwtmp(1) = 0.0d0
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
gwtmp(ig) = gwork(ig) * ( 7.d0 * wnl(ig,ll,is) - wnla(ig,ll,is) )
END DO
al = alpha(kk)
be = beta(kk)
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
fg = 0.0d0
gmod = SQRT( gv%hg_l(ig) )
gmod = SQRT( g(ig) )
DO s = 1, 3
fg = fg + 3.0d0/5.0d0 * fm(be,s,s,m) * gv%kgx_l(al,ig,1) / gmod
fg = fg + 3.0d0/5.0d0 * fm(be,s,s,m) * gkx_l(al,ig,1) / gmod
END DO
DO s = 1, 3
fg = fg + 6.0d0/5.0d0 * fm(be,s,al,m) * gv%kgx_l(s,ig,1) / gmod
fg = fg + 6.0d0/5.0d0 * fm(be,s,al,m) * gkx_l(s,ig,1) / gmod
END DO
gwtmp(ig) = gwtmp(ig) - fg * wnl(ig,ll,is)
END DO
DO ia= 1 , na(is)
auxc(1,ia) = CMPLX(0.0d0,0.0d0)
DO ig = gv%gstart, gv%ngw_l
DO ig = gstart, ngw
auxc(ig,ia) = csign(3) * gwtmp(ig) * eigr(ig,ia+iss-1)
END DO
END DO
CALL DGEMM( 'T', 'N', na(is), nx, 2*gv%ngw_l, 1.0d0, &
auxc(1,1), 2*gv%ngw_l, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, &
CALL DGEMM( 'T', 'N', na(is), nx, 2*ngw, 1.0d0, &
auxc(1,1), 2*ngw, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, &
0.0d0, fnls(iss,1), nsanl)
DEALLOCATE(auxc)
@ -646,18 +647,18 @@
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE pseudo_stress(deps, epseu, gv, gagx_l, sfac, dvps, rhoeg, ht)
SUBROUTINE pseudo_stress(deps, epseu, gagx_l, sfac, dvps, rhoeg, ht)
! (describe briefly what this routine does...)
! ----------------------------------------------
! ... declare modules
USE cell_module, only: boxdimensions
USE ions_base, ONLY: nsp
USE cp_types, ONLY: recvecs
USE cell_module, only: boxdimensions
USE ions_base, ONLY: nsp
USE reciprocal_vectors, ONLY: gstart, gzero
USE gvecp, ONLY: ngm
! ... declare subroutine arguments
TYPE (recvecs), INTENT(IN) :: gv
TYPE (boxdimensions), INTENT(IN) :: ht
REAL(dbl), INTENT(OUT) :: deps(:)
REAL(dbl), INTENT(IN) :: gagx_l(:,:)
@ -679,12 +680,12 @@
depst = (0.d0,0.d0)
DO is = 1, nsp
DO ig = gv%gstart, gv%ng_l
DO ig = gstart, ngm
rhets = rhoeg(ig, 1)
IF( nspin > 1) THEN
rhets = rhets + rhoeg(ig, 2)
END IF
rhets = 2.d0 * sfac( is, ig ) * dvps(ig,is) * CONJG(rhets)
rhets = 2.d0 * sfac( ig, is ) * dvps(ig,is) * CONJG(rhets)
depst(1) = depst(1) + rhets * gagx_l(1,ig)
depst(2) = depst(2) + rhets * gagx_l(2,ig)
depst(3) = depst(3) + rhets * gagx_l(3,ig)
@ -694,7 +695,7 @@
END DO
END DO
IF(gv%gzero) THEN
IF(gzero) THEN
deps = 2.0_dbl * omega * REAL(depst) - epseu * dalbe
ELSE
deps = 2.0_dbl * omega * REAL(depst)
@ -709,7 +710,7 @@
! BEGIN manual
SUBROUTINE stress_kin(dekin, c0, cdesc, occ, gagx_l, gv)
SUBROUTINE stress_kin(dekin, c0, cdesc, occ, gagx_l)
! this routine computes the kinetic energy contribution to the stress
! tensor
@ -724,8 +725,8 @@
USE gvecw, ONLY: tecfix, gcsig, gcfix, gcutz
USE wave_types, ONLY: wave_descriptor
USE constants, ONLY: pi
USE cp_types, ONLY: recvecs
USE control_flags, ONLY: force_pairing
USE reciprocal_vectors, ONLY: gstart, g
IMPLICIT NONE
@ -734,7 +735,6 @@
COMPLEX(dbl), INTENT(IN) :: c0(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: cdesc
REAL(dbl), INTENT(IN) :: occ(:,:,:)
TYPE (recvecs), INTENT(IN) :: gv
REAL(dbl) gagx_l(:,:)
! ... declare other variables
@ -749,9 +749,9 @@
dekin = 0.0_dbl
cost1 = 2.0_dbl / SQRT(pi)
ALLOCATE( arg( cdesc%ldg ) )
DO ig = gv%gstart, cdesc%ngwl
DO ig = gstart, cdesc%ngwl
IF(tecfix) THEN
arg(ig) = 1.0_dbl + cost1 * gcutz * exp(-((gv%hg_l(ig)-gcfix)/gcsig)**2)/gcsig
arg(ig) = 1.0_dbl + cost1 * gcutz * exp(-((g(ig)-gcfix)/gcsig)**2)/gcsig
ELSE
arg(ig) = 1.0_dbl
END IF
@ -763,7 +763,7 @@
IF( force_pairing ) ispin_wfc = 1
DO ib = 1, cdesc%nbl( ispin )
sk = 0.0_dbl
DO ig = gv%gstart, cdesc%ngwl
DO ig = gstart, cdesc%ngwl
scg = arg(ig) * CONJG( c0(ig,ib,1,ispin_wfc) ) * c0(ig,ib,1,ispin_wfc)
sk(1) = sk(1) + scg * gagx_l(1,ig)
sk(2) = sk(2) + scg * gagx_l(2,ig)
@ -785,20 +785,21 @@
!== COMPUTES HARTREE ENERGY CONTRIBUTION ==
!=======================================================================
SUBROUTINE STRESS_HAR(DEHT, EHR, sfac, PS, RHOEG, GAgx_L, GV, box )
SUBROUTINE STRESS_HAR(DEHT, EHR, sfac, PS, RHOEG, GAgx_L, box )
use ions_base, only: nsp
USE cell_module, only: boxdimensions
use mp_global, ONLY: mpime, nproc
USE constants, ONLY: fpi
USE cell_base, ONLY: tpiba2
USE cp_types, ONLY: recvecs, pseudo, pseudo_ncpp
use ions_base, only: nsp
USE cell_module, only: boxdimensions
use mp_global, ONLY: mpime, nproc
USE constants, ONLY: fpi
USE cell_base, ONLY: tpiba2
USE cp_types, ONLY: pseudo, pseudo_ncpp
USE reciprocal_vectors, ONLY: gstart, g
USE gvecp, ONLY: ngm
IMPLICIT NONE
!---------------------------------------------------ARGUMENT
type (recvecs) :: gv
type (boxdimensions) :: box
TYPE (pseudo), INTENT(IN) :: ps
REAL(dbl) :: DEHT(:), EHR, GAgx_L(:,:)
@ -829,14 +830,14 @@
DEHC = (0.D0,0.D0)
DEHT = 0.D0
DO IG = gv%gstart, gv%NG_L
DO IG = gstart, ngm
RHOP = (0.D0,0.D0)
RHOPR= (0.D0,0.D0)
DO IS = 1, NSP
RHOP = RHOP + sfac( is, IG ) * ps%RHOPS(IG,is)
RHOPR = RHOPR + sfac( is, IG ) * ps%RHOPS(IG,is) * ps%ap(is)%RAGGIO**2 * 0.5D0
RHOP = RHOP + sfac( IG, is ) * ps%RHOPS(IG,is)
RHOPR = RHOPR + sfac( IG, is ) * ps%RHOPS(IG,is) * ps%ap(is)%RAGGIO**2 * 0.5D0
END DO
HGM1 = 1.D0 / gv%HG_L(IG) / TPIBA2
HGM1 = 1.D0 / g(IG) / TPIBA2
RHET = 0.0_dbl
DO ispin = 1, nspin
RHET = RHET + RHOEG(ig,ispin)

11
CPV/wave.f90
View File

@ -71,7 +71,7 @@
! ----------------------------------------------
REAL(dbl) FUNCTION dft_kinetic_energy(c0, cdesc, gv, kp, tecfix, f, rsum, xmkin)
REAL(dbl) FUNCTION dft_kinetic_energy(c0, cdesc, kp, tecfix, f, rsum, xmkin)
! This function compute the Total Quanto-Mechanical Kinetic Energy of the Kohn-Sham
! wave function
@ -81,15 +81,14 @@
USE brillouin, ONLY: kpoints
USE wave_types, ONLY: wave_descriptor
USE electrons_module, ONLY: pmss
USE cp_types, ONLY: recvecs
USE control_flags, ONLY: force_pairing
USE reciprocal_space_mesh, ONLY: gkcutz_l, gk_l
IMPLICIT NONE
COMPLEX(dbl), INTENT(IN) :: c0(:,:,:,:) ! wave functions coefficients
TYPE (wave_descriptor), INTENT(IN) :: cdesc ! descriptor of c0
REAL(dbl), INTENT(IN) :: f(:,:,:) ! occupation numbers
TYPE (recvecs), INTENT(IN) :: gv ! reciprocal space vectors
TYPE (kpoints), INTENT(IN) :: kp ! k points
LOGICAL, INTENT(IN) :: tecfix ! Constant Cut-off is used
REAL(dbl), INTENT(OUT) :: rsum(:) ! charge density
@ -104,7 +103,7 @@
IF( ( cdesc%nkl > SIZE( c0, 3 ) ) .OR. &
( cdesc%nkl > SIZE( kp%weight ) ) .OR. &
( cdesc%nkl > SIZE( gv%khg_l, 2 ) ) .OR. &
( cdesc%nkl > SIZE( gk_l, 2 ) ) .OR. &
( cdesc%nkl > SIZE( f, 2 ) ) ) &
CALL errore( ' dft_kinetic_energy ', ' wrong arrays sizes ', 1 )
@ -126,9 +125,9 @@
END IF
IF( tecfix ) THEN
gmod2 => gv%khgcutz_l(:,ik)
gmod2 => gkcutz_l(:,ik)
ELSE
gmod2 => gv%khg_l(:,ik)
gmod2 => gk_l(:,ik)
ENDIF
xkin = xkin + fact * &

40
CPV/waveinit.f90
View File

@ -15,10 +15,10 @@
! (describe briefly what this module does...)
! ----------------------------------------------
! routines in this module:
! SUBROUTINE pw_rand_init(nbeg,cm,c0,gv)
! SUBROUTINE calphi(pslm,philm,gv,eigr,nchan)
! SUBROUTINE calpslm(pslm,gv,nchan,ns,np,nd)
! SUBROUTINE pw_atomic_init(nbeg,cm,c0,gv,eigr)
! SUBROUTINE pw_rand_init(nbeg,cm,c0)
! SUBROUTINE calphi(pslm,philm,eigr,nchan)
! SUBROUTINE calpslm(pslm,nchan,ns,np,nd)
! SUBROUTINE pw_atomic_init(nbeg,cm,c0,eigr)
! ----------------------------------------------
! END manual
@ -42,7 +42,7 @@
! subroutines
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE pw_rand_init(nbeg, cm, c0, wfill, ce, wempt, gv, kp)
SUBROUTINE pw_rand_init(nbeg, cm, c0, wfill, ce, wempt, kp)
! this routine sets the initial wavefunctions at random
! ----------------------------------------------
@ -55,9 +55,9 @@
USE mp_global, ONLY: mpime, nproc, root
USE reciprocal_vectors, ONLY: ig_l2g, ngw, ngwt
USE brillouin, ONLY: kpoints
USE cp_types, ONLY: recvecs
USE io_base, ONLY: stdout
USE control_flags, ONLY: force_pairing
USE reciprocal_space_mesh, ONLY: gkmask_l
IMPLICIT NONE
@ -66,7 +66,6 @@
! ... declare subroutine arguments
COMPLEX(dbl), INTENT(OUT) :: cm(:,:,:,:), c0(:,:,:,:), ce(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: wfill, wempt
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
INTEGER, INTENT(IN) :: nbeg
@ -121,32 +120,13 @@
DO ispin = 1, nspin
DO ik = 1, SIZE( cm, 3)
CALL wave_rand_init( cm( :, :, ik, ispin ) )
! ntest = ngwt / 4
! IF( ntest < wfill%nbt( ispin ) ) THEN
! ntest = ngwt
! END IF
!! ... assign random values to wave functions
! DO ib = 1, SIZE( cm, 2 )
! pwt( : ) = 0.0d0
! DO ig = 3, ntest
! rranf1 = 0.5d0 - rranf()
! rranf2 = rranf()
! pwt( ig ) = ampre * DCMPLX(rranf1, rranf2)
! END DO
! CALL splitwf ( cm( :, ib, ik, ispin ), pwt, ngw, ig_l2g, mpime, nproc, 0 )
! END DO
IF ( .NOT. kp%gamma_only ) THEN
! .. . set to zero all elements outside the cutoff sphere
DO ib = 1, SIZE( cm, 2 )
cm(:, ib, ik, ispin) = cm(:, ib, ik, ispin) * gv%kg_mask_l(:,ik)
cm(:, ib, ik, ispin) = cm(:, ib, ik, ispin) * gkmask_l(:,ik)
END DO
END IF
END DO
! DO ik = 1, SIZE( cm, 3 )
! IF ( gv%gzero ) THEN
! cm( 1, :, ik, ispin ) = (0.0d0, 0.0d0)
! END IF
! END DO
END DO
DEALLOCATE( pwt )
@ -185,7 +165,7 @@
SUBROUTINE pw_atomic_init(nbeg, cm, c0, wfill, ce, wempt, gv, kp, eigr)
SUBROUTINE pw_atomic_init(nbeg, cm, c0, wfill, ce, wempt, kp, eigr)
! (describe briefly what this routine does...)
! ----------------------------------------------
@ -197,13 +177,11 @@
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE brillouin, ONLY: kpoints
USE cp_types, ONLY: recvecs
IMPLICIT NONE
! ... declare subroutine arguments
INTEGER, INTENT(IN) :: nbeg
TYPE (recvecs), INTENT(IN) :: gv
TYPE (kpoints), INTENT(IN) :: kp
COMPLEX(dbl), INTENT(OUT) :: cm(:,:,:,:), c0(:,:,:,:), ce(:,:,:,:)
TYPE (wave_descriptor), INTENT(IN) :: wfill, wempt
@ -214,7 +192,7 @@
! ... end of declarations
! ----------------------------------------------
CALL pw_rand_init(nbeg, cm, c0, wfill, ce, wempt, gv, kp)
CALL pw_rand_init(nbeg, cm, c0, wfill, ce, wempt, kp)
IF( nbeg < 0 ) THEN
IF( ionode ) THEN

2
Modules/fft_base.f90
View File

@ -45,7 +45,7 @@
INTEGER, ALLOCATABLE :: stmask(:)
REAL(dbl), ALLOCATABLE :: fft_timing(:,:)
REAL(dbl) :: fft_timing( 4, 2 ) = 0.0d0
!=----------------------------------------------------------------------=!