Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Alignment with VASP #1

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Alignment with VASP #1

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4f021f6
Fix interface with VASP
martin-schlipf Feb 7, 2023
519ce88
Update INCAR reader to simplify code
martin-schlipf Feb 7, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
150 changes: 67 additions & 83 deletions extpot.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
module mextpot
use prec
implicit none

type extpot
logical :: lextpot, dextpot
integer :: nx, ny, nz ! should be equal to the dimension of potential grid
Expand All @@ -22,12 +22,23 @@ module mextpot
! parameters for real space VLM projection
integer, parameter :: KMAX=20, di_smear=5, NGRID=50
real(q), parameter :: Delta = 0.1_q

contains

function EXTPT_LEXTPOT()
logical :: EXTPT_LEXTPOT
EXTPT_LEXTPOT=EXTPT%lextpot
end function EXTPT_LEXTPOT

function EXTPT_DEXTPOT()
logical :: EXTPT_DEXTPOT
EXTPT_DEXTPOT=EXTPT%dextpot
end function EXTPT_DEXTPOT

! read INCAR FILE for external potential information
subroutine EXTPT_READER(NGXF,NGYF,NGZF,IU0,IU5,IU6)
use base
use reader_tags
use vaspxml
implicit none

Expand All @@ -48,41 +59,13 @@ subroutine EXTPT_READER(NGXF,NGYF,NGZF,IU0,IU5,IU6)
EXTPT%nx = -1
EXTPT%ny = -1
EXTPT%nz = -1
CALL RDATAB(lopen,INCAR,IU5,'LEXTPOT','=','#',';','L', & ! calling rdatab subroutine for reading INCAR and searching for LEXTPOT
& IDUM,RDUM,CDUM,EXTPT%lextpot,CHARAC,N,1,IERR)
IF (((IERR/=0).AND.(IERR/=3)).OR. &
& ((IERR==0).AND.(N<1))) THEN
IF (IU0>=0) &
WRITE(IU0,*)'Error reading item ''LEXTPOT'' from file INCAR.'
EXTPT%lextpot = .FALSE.
RETURN
ENDIF
CALL RDATAB(lopen,INCAR,IU5,'DEXTPOT','=','#',';','L', & ! calling rdatab subroutine for reading INCAR and searching for DEXTPOT
& IDUM,RDUM,CDUM,EXTPT%dextpot,CHARAC,N,1,IERR)
IF (((IERR/=0).AND.(IERR/=3)).OR. &
& ((IERR==0).AND.(N<1))) THEN
IF (IU0>=0) &
WRITE(IU0,*)'Error reading item ''DEXTPOT'' from file INCAR.'
EXTPT%dextpot = .FALSE.
RETURN
ENDIF
CALL XML_INCAR('LEXTPOT','L',IDUM,RDUM,CDUM,EXTPT%lextpot,CHARAC,N)
CALL XML_INCAR('DEXTPOT','L',IDUM,RDUM,CDUM,EXTPT%dextpot,CHARAC,N)

CALL PROCESS_INCAR(INCAR_F, 'LEXTPOT', EXTPT%lextpot)
CALL PROCESS_INCAR(INCAR_F, 'DEXTPOT', EXTPT%dextpot)
if( EXTPT%lextpot ) then
CALL RDATAB(lopen,INCAR,IU5,'FCORR','=','#',';','F', &
IDUM,EXTPT%FCORR,CDUM,LDUM,CHARAC,N,1,IERR)
IF (((IERR/=0).AND.(IERR/=3)).OR. &
& ((IERR==0).AND.(N<1))) THEN
IF (IU0>=0) &
WRITE(IU0,*)'Error reading item ''FCORR'' from file INCAR.'
EXTPT%FCORR = 0.67_q
RETURN
ENDIF
CALL XML_INCAR('FCORR','F',IDUM,EXTPT%FCORR,CDUM,LDUM,CHARAC,N)
CALL PROCESS_INCAR(INCAR_F, 'FCORR', EXTPT%FCORR)
endif

if( EXTPT%lextpot ) then
if( EXTPT%lextpot ) then
open(unit=IU_POT,file='EXTPOT',status='old')
read(IU_POT,*) EXTPT%nx, EXTPT%ny, EXTPT%nz
! check grid size consistency
Expand All @@ -95,7 +78,7 @@ subroutine EXTPT_READER(NGXF,NGYF,NGZF,IU0,IU5,IU6)
write(IU0,*) ' NX, NY, NZ:', EXTPT%nx, EXTPT%ny, EXTPT%nx
stop
endif
allocate(EXTPT_GRID(EXTPT%nx, EXTPT%ny, EXTPT%nz))
allocate(EXTPT_GRID(EXTPT%nx, EXTPT%ny, EXTPT%nz))
! follow the tradition of CHGCAR format, ix be the fastest index
read(IU_POT,*) (((EXTPT_GRID(ix,iy,iz), ix=1,EXTPT%nx), &
iy=1,EXTPT%ny), iz=1,EXTPT%nz)
Expand All @@ -110,7 +93,7 @@ end subroutine EXTPT_READER
! In realmode, CVTOT allocated as complex array, but used as real grid
!***********************************************************************

SUBROUTINE EXTERNAL_POT_ADD(GRIDC, LATT_CUR, CVTOT)
SUBROUTINE EXTPT_EXTERNAL_POT_ADD(GRIDC, LATT_CUR, CVTOT)

USE prec
USE base
Expand All @@ -129,15 +112,44 @@ SUBROUTINE EXTERNAL_POT_ADD(GRIDC, LATT_CUR, CVTOT)

call EXTERNAL_POT_ADD_(GRIDC, LATT_CUR, CVTOT(1))

end subroutine EXTERNAL_POT_ADD
end subroutine EXTPT_EXTERNAL_POT_ADD

subroutine EXTPT_EXTERNAL_POT_ADD_PAW(POTAE, POT, NDIM, NCDIJ, LMMAX, NIP)

USE prec
USE base
USE lattice
USE mpimy
USE mgrid
USE poscar
USE constant
implicit none

INTEGER NDIM, NCDIJ, LMMAX, NIP
REAL(q) POT(:,:,:), POTAE(:,:,:)

integer ISP, LM, II

do ISP = 1, NCDIJ
do LM = 1, LMMAX
do II = 1, NDIM
POTAE(II,LM,ISP) = POTAE(II,LM,ISP) + VLM(II,LM,NIP) !+1.0_q
POT(II,LM,ISP) = POT(II,LM,ISP) + VLM(II,LM,NIP) !+1.0_q
enddo
enddo
enddo

end subroutine EXTPT_EXTERNAL_POT_ADD_PAW



!*********************** SUBROUTINE CALC_VLM ***************************
! This subroutine calculates the LM component of the external potential
! This subroutine calculates the LM component of the external potential
! This code was written assuming realmode and NGZhalf = .True.
!***********************************************************************

SUBROUTINE CALC_VLM(LATT_CUR, GRIDC, T_INFO, P, LMDIM, WDES, IU0)
SUBROUTINE EXTPT_CALC_VLM(LATT_CUR, GRIDC, T_INFO, P, LMDIM, WDES, IU0)

use prec
use base
use lattice
Expand Down Expand Up @@ -181,7 +193,7 @@ SUBROUTINE CALC_VLM(LATT_CUR, GRIDC, T_INFO, P, LMDIM, WDES, IU0)

! copy the external potential to CVEXT, change normal layout to VASP layout
CVEXT = 0.0_q
call EXTERNAL_POT_ADD_(GRIDC, LATT_CUR, CVEXT(1))
call EXTERNAL_POT_ADD_(GRIDC, LATT_CUR, CVEXT(1))

! set up the maximum radial grid dimension, copied from SET_DD_PAW
NDIM=0
Expand Down Expand Up @@ -223,7 +235,7 @@ SUBROUTINE CALC_VLM(LATT_CUR, GRIDC, T_INFO, P, LMDIM, WDES, IU0)
, IRDMAX, INDMAX, .True.)

! start projection
VLMFULL = 0.0_q
VLMFULL = 0.0_q
! do ion loops
ion1: do NI=1,T_INFO%NIONS
RI(:) = T_INFO%POSION(1:3,NI) ! ion position in direct coordinate
Expand Down Expand Up @@ -288,7 +300,7 @@ SUBROUTINE CALC_VLM(LATT_CUR, GRIDC, T_INFO, P, LMDIM, WDES, IU0)
LM=L*L+M+1
VLM(:,LM,NIP) = VLMFULL(:,LM,NI)
enddo
enddo
enddo
enddo ion2
CALLMPI( M_sum_d(WDES%COMM_INTER, VLM, NDIM*LMMAX*WDES%NIONS) )

Expand Down Expand Up @@ -316,14 +328,14 @@ SUBROUTINE CALC_VLM(LATT_CUR, GRIDC, T_INFO, P, LMDIM, WDES, IU0)
deallocate(VLMFULL, CVEXT)

return
END SUBROUTINE CALC_VLM
END SUBROUTINE EXTPT_CALC_VLM

!*************************** subroutine SETYLM_FINE_GRID **************************
! this subroutine setup the fine grid for the real space projection integeration
! I made up this layout to parallel the real space projection code
! Use LPAR to choose whether to distritute the local fine grid points or not
! Note:
! In CALC_VLM, we can not parallelize over uniform grid (GRIDC), so set LPAR = .Ture.
! In CALC_VLM, we can not parallelize over uniform grid (GRIDC), so set LPAR = .Ture.
! to parallelize over the local fine grid
! In DE_DVRN, parallelize over GRIDC, so set LPAR = .False., making the local
! fine grid layout serial.
Expand All @@ -344,7 +356,7 @@ subroutine SETYLM_FINE_GRID( GRIDC, NGRID, YLM, LYMAX, LMMAX, &
real(q) :: YLM(IRDMAX,LMMAX), DIST(IRDMAX)
integer :: ix, iy, iz, NODE_ME, NCPU, IND
logical :: LPAR

#ifdef MPI
NCPU = GRIDC%COMM%NCPU
NODE_ME = GRIDC%COMM%NODE_ME
Expand Down Expand Up @@ -475,7 +487,7 @@ end subroutine BSPLINE_INTERPOLATE_COEFF
! Calculate the 1-D B-spline coefficients for a "delta" grid
!***********************************************************************
subroutine CALC_Sdelta_COEFF(N, s, di_smear)
use prec
use prec
use constant
implicit none
integer :: N, di_smear
Expand Down Expand Up @@ -582,7 +594,7 @@ END FUNCTION RK
!************************ FUNCTION FSWITCH *****************************
!
! This function calculates the switch function f(r) for V(r)
! In principle, can be anything as long as:
! In principle, can be anything as long as:
! f(RMAX) = 1
! f'(RMAX) = 0
! f(RMAX+Delta) = 0
Expand All @@ -605,34 +617,6 @@ FUNCTION FSWITCH( R, RMAX, Delta )
return
end function FSWITCH

subroutine ADD_EXTERNAL_POTENTIAL(POTAE, POT, NDIM, NCDIJ, LMMAX, NIP)

USE prec
USE base
USE lattice
USE mpimy
USE mgrid
USE poscar
USE constant
implicit none

INTEGER NDIM, NCDIJ, LMMAX, NIP
REAL(q) POT(:,:,:), POTAE(:,:,:)

integer ISP, LM, II

do ISP = 1, NCDIJ
do LM = 1, LMMAX
do II = 1, NDIM
POTAE(II,LM,ISP) = POTAE(II,LM,ISP) + VLM(II,LM,NIP) !+1.0_q
POT(II,LM,ISP) = POT(II,LM,ISP) + VLM(II,LM,NIP) !+1.0_q
enddo
enddo
enddo

end subroutine ADD_EXTERNAL_POTENTIAL


!************************* ENERGY_DERIVATIVE ***************************
! This subroutine augments the electron density to give dE/dV(r)
! Theoretically it should be the AE density
Expand All @@ -643,7 +627,7 @@ end subroutine ADD_EXTERNAL_POTENTIAL
! defined in GRID_SOFT
! CHTOT: tilde(n) + hat(n) as input, correct dE/dV as output
!***********************************************************************
SUBROUTINE DE_DVRN(WDES, GRID_SOFT, GRIDC, SOFT_TO_C, LATT_CUR, &
SUBROUTINE EXTPT_DE_DVRN(WDES, GRID_SOFT, GRIDC, SOFT_TO_C, LATT_CUR, &
P, T_INFO, LMDIM, CRHODE, CHTOT_, CHDEN, IRDMAX )
use prec
use base
Expand Down Expand Up @@ -817,7 +801,7 @@ SUBROUTINE DE_DVRN(WDES, GRID_SOFT, GRIDC, SOFT_TO_C, LATT_CUR, &
l_loop: DO L =1,P(NT)%LMAX
LMP=LM
lp_loop: DO LP=L,P(NT)%LMAX

if ( .not. LINFO ) then
! setup QIJ
call SETQIJ(GRIDC, L,LP,P(NT),QIJ, NDIM)
Expand Down Expand Up @@ -864,7 +848,7 @@ SUBROUTINE DE_DVRN(WDES, GRID_SOFT, GRIDC, SOFT_TO_C, LATT_CUR, &
!*****************************
! compute and write DERINFO
!*****************************
DO IC=ISTART,IEND-1
DO IC=ISTART,IEND-1
JLL = JL(IC)
JLM = JS(IC)

Expand All @@ -875,14 +859,14 @@ SUBROUTINE DE_DVRN(WDES, GRID_SOFT, GRIDC, SOFT_TO_C, LATT_CUR, &
do INDF = 1, INDFMAX
R = DISTF(INDF)*RCUT ! in Angstrom
if ( R <= RCUT ) then
RVECF(1) = XSF(INDF)*R
RVECF(2) = YSF(INDF)*R
RVECF(1) = XSF(INDF)*R
RVECF(2) = YSF(INDF)*R
RVECF(3) = ZSF(INDF)*R
RVECF = matmul( transpose(LATT_CUR%B(1:3,1:3)), RVECF(1:3) )
DRVEC = RVECF - RVEC ! in direct
if(abs(DRVEC(1)*NXMAX)<di_smear+1 .and. &
abs(DRVEC(2)*NYMAX)<di_smear+1 .and. &
abs(DRVEC(3)*NZMAX)<di_smear+1) then
abs(DRVEC(3)*NZMAX)<di_smear+1) then
RVECF(:) = RVECF(:) + RI(:) ! grid coordinate in direct
val = bspline( RVECF, SGRID, NXMAX, NYMAX, NZMAX )
VLMK(K) = VLMK(K) + dV* FSWITCH(R,RMAX,Delta)* &
Expand Down Expand Up @@ -947,7 +931,7 @@ SUBROUTINE DE_DVRN(WDES, GRID_SOFT, GRIDC, SOFT_TO_C, LATT_CUR, &
close( ofile )

return
END SUBROUTINE DE_DVRN
END SUBROUTINE EXTPT_DE_DVRN

end module mextpot

Expand Down