Skip to content

Commit

Permalink
Starting from this version, the lattice zero-point energy will be rem…
Browse files Browse the repository at this point in the history 
…oved from the total energy value.
  • Loading branch information
@wangy2014
wangy2014 committed May 5, 2022
1 parent 83c9d4c commit 18ac767
Show file tree
Hide file tree
Showing 7 changed files with 300 additions and 136 deletions.
100 changes: 84 additions & 16 deletions MST/src/BookKeepingModule.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,10 @@ module BookKeepingModule
writeRow, &
getHeadLineValue, &
getLastColumnValue
!
interface writeHeadLine
module procedure writeHL0, writeHL1
end interface
!
interface writeRow
module procedure writeRow0, writeRow1, writeRow2
Expand All @@ -28,11 +32,12 @@ module BookKeepingModule
integer (kind=IntKind) :: num_rows = 0
integer (kind=IntKind) :: num_columns = 0
integer (kind=IntKind) :: num_columns_old = 0
integer (kind=IntKind), parameter :: MaxLen = 100
!
character (len=12) :: ColumnTitle(MaxColumns)
character (len=12) :: ColumnValue(MaxColumns)
character (len=40) :: HeadLineTitle(MaxHeadLines)
character (len=60) :: HeadLineValue(MaxHeadLines)
character (len=MaxLen) :: HeadLineTitle(MaxHeadLines)
character (len=MaxLen) :: HeadLineValue(MaxHeadLines)
!
integer (kind=IntKind) :: ColumnLength(MaxColumns)
!
Expand Down Expand Up @@ -84,21 +89,31 @@ subroutine initBookKeeping(fpath,short_label,myid)
! ----------------------------------------------------------------
LOOP_1: do
read(funit,'(a)',iostat=status)text
text = adjustl(text)
if (status < 0 .or. rflag > 2) then
exit LOOP_1
else if (text(1:25) == '=========================') then
rflag = rflag + 1
if (rflag == 0) then
rflag = 1
endif
else if (text(1:25) == '-------------------------') then
rflag = rflag + 1
if (rflag == 1) then
rflag = 2
endif
else if (rflag == 0) then
num_head_lines = num_head_lines + 1
id = index(text,':')
if (id < 2 .or. id > len(text)) then
call WarningHandler('initBookKeeping','Invalid headline',text)
exit LOOP_1
if (text(1:1) /= '*') then
id = index(text,':')
if (id < 2 .or. id > len(text)) then
call WarningHandler('initBookKeeping','Invalid headline',text)
exit LOOP_1
endif
HeadLineTitle(num_head_lines) = adjustl(text(:id-1))
HeadLineValue(num_head_lines) = adjustl(text(id+1:))
else
HeadLineTitle(num_head_lines) = trim(text)
HeadLineValue(num_head_lines) = ' '
endif
HeadLineTitle(num_head_lines) = adjustl(text(:id-1))
HeadLineValue(num_head_lines) = adjustl(text(id+1:))
else if (rflag == 1) then
! ----------------------------------------------------------
call setString(text)
Expand Down Expand Up @@ -238,7 +253,7 @@ function getLastColumnValue(title) result(v)
!
character (len=*), intent(in) :: title
!
character (len=60) :: v
character (len=MaxLen) :: v
!
integer (kind=IntKind) :: i
!
Expand Down Expand Up @@ -266,15 +281,17 @@ end function getLastColumnValue
! *******************************************************************
!
! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine writeHeadLine(title,value)
subroutine writeHL0(title,value)
! ===================================================================
implicit none
!
character (len=*), intent(in) :: title
character (len=*), intent(in) :: value
!
character (len=30) :: tt
character (len=48) :: tv
integer (kind=IntKind), parameter :: colon_position = 38
!
character (len=colon_position-1) :: tt
character (len=50) :: tv
!
if (FileExist) then
return
Expand All @@ -293,7 +310,57 @@ subroutine writeHeadLine(title,value)
HeadLineValue(num_head_lines) = tv
call FlushFile(funit)
!
end subroutine writeHeadLine
end subroutine writeHL0
! ===================================================================
!
! *******************************************************************
!
! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine writeHL1(title,num_appears)
! ===================================================================
implicit none
!
character (len=*), intent(in) :: title
integer, intent(in), optional :: num_appears
integer (kind=IntKind) :: np, i, j, slen
!
character (len=MaxLen) :: tt
!
if (FileExist) then
return
else if (present(num_appears)) then
if (num_appears < 1) then
call ErrorHandler('writeHeadLine', &
'Number of appears of the text pattern < 1',num_appears)
endif
np = num_appears
else
np = 1
endif
!
slen = len_trim(title)
tt = ' '
do i = 1, np
j = (i-1)*slen
if (j+slen > MaxLen) then
exit
endif
tt(j+1:j+slen) = trim(title)
enddo
!
write(funit,'(a)')trim(tt)
!
num_head_lines = num_head_lines + 1
!
if (num_head_lines > MaxHeadLines) then
call ErrorHandler('writeHeadLine', &
'Number of headlines exceed upper limit',MaxHeadLines)
endif
HeadLineTitle(num_head_lines) = trim(tt)
HeadLineValue(num_head_lines) = ' '
call FlushFile(funit)
!
end subroutine writeHL1
! ===================================================================
!
! *******************************************************************
Expand All @@ -316,7 +383,8 @@ subroutine insertColumn(title,value)
'Number of columns exceed the upper limit',MaxColumns)
endif
!
nt = len_trim(title)
! nt = len_trim(title)
nt = len(title)
nv = len_trim(value)
n = max(nt,nv)
!
Expand Down
52 changes: 28 additions & 24 deletions MST/src/PotentialGenerationModule.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2958,22 +2958,22 @@ subroutine printMadelungShiftTable(iter,fu)
!
if ((isLSMS().or.isKKR()) .and. max_shells > 0 .and. printNeighbor) then
if (max_shells == 1) then
write(fu,'(87(''=''))')
write(fu,'(a)')' Atom Index Q Qmt dQ Vmad Local Energy dQ(1nn)'
! write(fu,'(87(''-''))')
write(fu,'(107(''=''))')
write(fu,'(a)')' Atom Index Q Qmt dQ Vmad Local Energy dQ(1nn) LocalEnergy-ZptE'
! write(fu,'(107(''-''))')
else
write(fu,'(99(''=''))')
write(fu,'(a)')' Atom Index Q Qmt dQ Vmad Local Energy dQ(1nn) dQ(2nn)'
! write(fu,'(99(''-''))')
write(fu,'(119(''=''))')
write(fu,'(a)')' Atom Index Q Qmt dQ Vmad Local Energy dQ(1nn) dQ(2nn) LocalEnergy-ZptE'
! write(fu,'(119(''-''))')
endif
else if (na == 1) then
write(fu,'(75(''=''))')
write(fu,'(a)')' Atom Index Q Qmt dQ Vmad Local Energy'
! write(fu,'(75(''-''))')
write(fu,'(95(''=''))')
write(fu,'(a)')' Atom Index Q Qmt dQ Vmad Local Energy LocalEnergy-ZptE'
! write(fu,'(95(''-''))')
else
write(fu,'(80(''=''))')
write(fu,'(a)')' Atom Site Species Q Qmt dQ Vmad Local Energy'
! write(fu,'(80(''-''))')
write(fu,'(100(''=''))')
write(fu,'(a)')' Atom Site Species Q Qmt dQ Vmad Local Energy LocalEnergy-ZptE'
! write(fu,'(100(''-''))')
endif
global_table_line => getGlobalTableLine()
Q_Table => getGlobalOnSiteElectronTableOld()
Expand All @@ -2982,48 +2982,52 @@ subroutine printMadelungShiftTable(iter,fu)
if ((isLSMS().or.isKKR()) .and. max_shells > 0 .and. printNeighbor) then
if (max_shells == 1) then
do ig = 1, GlobalNumAtoms
write(fu,'(2x,a3,2x,i6,4(2x,f9.5),6x,f12.5,2x,f10.5)') &
write(fu,'(2x,a3,2x,i6,4(2x,f9.5),6x,f12.5,2x,f10.5,6x,f12.5)') &
getAtomName(ig), ig, &
Q_Table(ig),Qmt_Table(ig), &
Q_Table(ig)-getAtomicNumber(ig), &
MadelungShiftTable(ig),atom_en(1,ig), &
NeighborChargeTable(1,ig)
NeighborChargeTable(1,ig), &
atom_en(1,ig)-atom_en(3,ig)
enddo
write(fu,'(87(''=''))')
write(fu,'(107(''=''))')
else
do ig = 1, GlobalNumAtoms
write(fu,'(2x,a3,2x,i6,4(2x,f9.5),6x,f12.5,2(2x,f10.5))') &
write(fu,'(2x,a3,2x,i6,4(2x,f9.5),6x,f12.5,2(2x,f10.5),6x,f12.5)') &
getAtomName(ig), ig, &
Q_Table(ig),Qmt_Table(ig), &
Q_Table(ig)-getAtomicNumber(ig), &
MadelungShiftTable(ig),atom_en(1,ig), &
NeighborChargeTable(1,ig), &
NeighborChargeTable(2,ig)
NeighborChargeTable(2,ig), &
atom_en(1,ig)-atom_en(3,ig)
enddo
write(fu,'(99(''=''))')
write(fu,'(119(''=''))')
endif
else if (na == 1) then
do ig = 1, GlobalNumAtoms
write(fu,'(2x,a3,2x,i6,4(2x,f9.5),6x,f12.5)') &
write(fu,'(2x,a3,2x,i6,4(2x,f9.5),6x,f12.5,6x,f12.5)') &
getAtomName(ig), ig, &
Q_Table(ig),Qmt_Table(ig), &
Q_Table(ig)-getAtomicNumber(ig), &
MadelungShiftTable(ig),atom_en(1,ig)
MadelungShiftTable(ig),atom_en(1,ig), &
atom_en(1,ig)-atom_en(3,ig)
enddo
write(fu,'(75(''=''))')
write(fu,'(95(''=''))')
else
do ig = 1, GlobalNumAtoms
do ia = 1, getNumAlloyElements(ig)
! write(fu,'(2x,a3,2x,i6,4(2x,f20.16))')getAtomName(ig,ia), ig, ia, &
lig = global_table_line(ig) + ia
write(fu,'(2x,a3,2x,i6,2x,i2,4(2x,f9.5),6x,f12.5)') &
write(fu,'(2x,a3,2x,i6,2x,i2,4(2x,f9.5),6x,f12.5,6x,f12.5)') &
getAtomName(ig,ia), ig, ia, &
Q_Table(lig),Qmt_Table(lig), &
Q_Table(lig)-getAtomicNumber(ig,ia), &
MadelungShiftTable(lig),atom_en(1,ig)
MadelungShiftTable(lig),atom_en(1,ig),&
atom_en(1,ig)-atom_en(3,ig)
enddo
enddo
write(fu,'(80(''=''))')
write(fu,'(100(''=''))')
endif
if (fu == 6) then
write(6,'(/)')
Expand Down
12 changes: 6 additions & 6 deletions MST/src/SystemModule.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -781,12 +781,12 @@ end subroutine copyCharArray2String
endif
!
if (isDataStorageExisting('Energy Pressure')) then
EnPres => getDataStorage('Energy Pressure',2,NumAtoms,RealMark)
EnPres => getDataStorage('Energy Pressure',4,NumAtoms,RealMark)
else
! ----------------------------------------------------------------
call createDataStorage('Energy Pressure',2*NumAtoms,RealType)
call createDataStorage('Energy Pressure',4*NumAtoms,RealType)
! ----------------------------------------------------------------
EnPres => getDataStorage('Energy Pressure',2,NumAtoms,RealMark)
EnPres => getDataStorage('Energy Pressure',4,NumAtoms,RealMark)
endif
!
if (isDataStorageExisting('RMS Info')) then
Expand Down Expand Up @@ -1614,12 +1614,12 @@ subroutine setEnPres(i,pen)
! ===================================================================
implicit none
integer (kind=IntKind), intent(in) :: i
real (kind=RealKind), intent(in) :: pen(2)
real (kind=RealKind), intent(in) :: pen(4)
!
if (i<1 .or. i>NumAtoms) then
call ErrorHandler('setForce','Invalid atom index',i)
endif
EnPres(1:2,i) = pen(1:2)
EnPres(1:4,i) = pen(1:4)
!
end subroutine setEnPres
! ===================================================================
Expand Down Expand Up @@ -2210,7 +2210,7 @@ subroutine resetSystemMovie()
endif
!
Force(1:3,1:NumAtoms) = ZERO
EnPres(1:2,1:NumAtoms) = ZERO
EnPres(1:4,1:NumAtoms) = ZERO
RMS(1:2,1:NumAtoms) = ZERO
if ( nspin > 2 ) then
RMS(3:4,1:NumAtoms) = ZERO
Expand Down
Loading

0 comments on commit 18ac767

Please sign in to comment.