Commit 244d0af6 authored by Gregor Michalicek's avatar Gregor Michalicek

Introduce Wrapper for Broyden history IO.

...this surely is not yet the final version.
parent 55c9ea91
......@@ -20,6 +20,7 @@ io/cdnpot_io_hdf.F90
io/cdnpot_io_common.F90
io/cdn_io.F90
io/pot_io.F90
io/broyd_io.F90
io/rw_inp.f90
io/rw_noco.f90
io/r_inpXML.F90
......
MODULE m_broyd_io
USE m_types
USE m_cdn_io
IMPLICIT NONE
CONTAINS
SUBROUTINE readLastIterInAndDiffDen(hybrid,vecLen,nextIter,alpha,inDenVec,diffDenVec)
TYPE(t_hybrid), INTENT(IN) :: hybrid
INTEGER, INTENT(IN) :: vecLen
INTEGER, INTENT(OUT) :: nextIter
REAL, INTENT(OUT) :: alpha
REAL, INTENT(OUT) :: inDenVec(vecLen), diffDenVec(vecLen)
INTEGER :: mode
CALL getIOMode(mode)
! At the moment broyden IO is mode independent
IF (hybrid%l_calhf) THEN
OPEN (57,file='hf_broyd',form='unformatted',status='unknown')
ELSE
OPEN (57,file='broyd',form='unformatted',status='unknown')
ENDIF
READ(57) nextIter,alpha,diffDenVec,inDenVec
CLOSE(57)
END SUBROUTINE readLastIterInAndDiffDen
SUBROUTINE writeLastIterInAndDiffDen(hybrid,vecLen,nextIter,alpha,inDenVec,diffDenVec)
TYPE(t_hybrid), INTENT(IN) :: hybrid
INTEGER, INTENT(IN) :: nextIter,vecLen
REAL, INTENT(IN) :: alpha
REAL, INTENT(IN) :: inDenVec(vecLen), diffDenVec(vecLen)
INTEGER :: mode
CALL getIOMode(mode)
! At the moment broyden IO is mode independent
IF (hybrid%l_calhf) THEN
OPEN (57,file='hf_broyd',form='unformatted',status='unknown')
ELSE
OPEN (57,file='broyd',form='unformatted',status='unknown')
ENDIF
WRITE(57) nextIter,alpha,diffDenVec,inDenVec
CLOSE(57)
END SUBROUTINE writeLastIterInAndDiffDen
SUBROUTINE readUVec(input,hybrid,vecLen,relIter,currentIter,uVec)
TYPE(t_input), INTENT(IN) :: input
TYPE(t_hybrid), INTENT(IN) :: hybrid
INTEGER, INTENT(IN) :: vecLen, relIter, currentIter
REAL, INTENT(OUT) :: uVec(vecLen)
INTEGER :: mode, npos
INTEGER*8 :: recLen
CALL getIOMode(mode)
! At the moment broyden IO is mode independent
recLen=(vecLen+1)*8
IF (hybrid%l_calhf) THEN
OPEN (59,file='hf_broyd.'//CHAR(input%imix+48),access='direct',&
recl=recLen,form='unformatted',status='unknown')
ELSE
OPEN (59,file='broyd.'//CHAR(input%imix+48),access='direct',&
recl=recLen,form='unformatted',status='unknown')
ENDIF
npos=currentIter+relIter-1
IF (currentIter.GT.input%maxiter+1) npos = MOD(currentIter-2,input%maxiter)+1
READ(59,rec=2*npos-1) uVec(:vecLen)
CLOSE(59)
END SUBROUTINE readUVec
SUBROUTINE writeUVec(input,hybrid,vecLen,currentIter,uVec)
TYPE(t_input), INTENT(IN) :: input
TYPE(t_hybrid), INTENT(IN) :: hybrid
INTEGER, INTENT(IN) :: vecLen, currentIter
REAL, INTENT(IN) :: uVec(vecLen)
INTEGER :: mode, npos
INTEGER*8 :: recLen
CALL getIOMode(mode)
! At the moment broyden IO is mode independent
recLen=(vecLen+1)*8
IF (hybrid%l_calhf) THEN
OPEN (59,file='hf_broyd.'//CHAR(input%imix+48),access='direct',&
recl=recLen,form='unformatted',status='unknown')
ELSE
OPEN (59,file='broyd.'//CHAR(input%imix+48),access='direct',&
recl=recLen,form='unformatted',status='unknown')
ENDIF
npos=currentIter-1
IF (currentIter.GT.input%maxiter+1) npos = MOD(currentIter-2,input%maxiter)+1
WRITE(59,rec=2*npos-1) uVec(:vecLen)
CLOSE(59)
END SUBROUTINE writeUVec
SUBROUTINE readVVec(input,hybrid,vecLen,relIter,currentIter,dfivi,vVec)
TYPE(t_input), INTENT(IN) :: input
TYPE(t_hybrid), INTENT(IN) :: hybrid
INTEGER, INTENT(IN) :: vecLen, relIter, currentIter
REAL, INTENT(OUT) :: dfivi
REAL, INTENT(OUT) :: vVec(vecLen)
INTEGER :: mode, npos
INTEGER*8 :: recLen
CALL getIOMode(mode)
! At the moment broyden IO is mode independent
recLen=(vecLen+1)*8
IF (hybrid%l_calhf) THEN
OPEN (59,file='hf_broyd.'//CHAR(input%imix+48),access='direct',&
recl=recLen,form='unformatted',status='unknown')
ELSE
OPEN (59,file='broyd.'//CHAR(input%imix+48),access='direct',&
recl=recLen,form='unformatted',status='unknown')
ENDIF
npos=currentIter+relIter-1
IF (currentIter.GT.input%maxiter+1) npos = MOD(currentIter-2,input%maxiter)+1
READ(59,rec=2*npos) vVec(:vecLen), dfivi
CLOSE(59)
END SUBROUTINE readVVec
SUBROUTINE writeVVec(input,hybrid,vecLen,currentIter,dfivi,vVec)
TYPE(t_input), INTENT(IN) :: input
TYPE(t_hybrid), INTENT(IN) :: hybrid
INTEGER, INTENT(IN) :: vecLen, currentIter
REAL, INTENT(IN) :: dfivi
REAL, INTENT(IN) :: vVec(vecLen)
INTEGER :: mode, npos
INTEGER*8 :: recLen
CALL getIOMode(mode)
! At the moment broyden IO is mode independent
recLen=(vecLen+1)*8
IF (hybrid%l_calhf) THEN
OPEN (59,file='hf_broyd.'//CHAR(input%imix+48),access='direct',&
recl=recLen,form='unformatted',status='unknown')
ELSE
OPEN (59,file='broyd.'//CHAR(input%imix+48),access='direct',&
recl=recLen,form='unformatted',status='unknown')
ENDIF
npos=currentIter-1
IF (currentIter.GT.input%maxiter+1) npos = MOD(currentIter-2,input%maxiter)+1
WRITE(59,rec=2*npos) vVec(:vecLen), dfivi
CLOSE(59)
END SUBROUTINE writeVVec
LOGICAL FUNCTION initBroydenHistory(input,hybrid, vecLen)
! Initializes a Broyden history
! returns true if there already exists a Broyden history
TYPE(t_input), INTENT(IN) :: input
TYPE(t_hybrid), INTENT(IN) :: hybrid
INTEGER, INTENT(IN) :: vecLen
INTEGER*8 :: recLen
LOGICAL :: l_exist
INQUIRE (file='broyd.'//CHAR(input%imix+48),exist=l_exist)
recLen=(vecLen+1)*8
IF (hybrid%l_calhf) THEN
OPEN (59,file='hf_broyd.'//CHAR(input%imix+48),access='direct',&
recl=recLen,form='unformatted',status='unknown')
ELSE
OPEN (59,file='broyd.'//CHAR(input%imix+48),access='direct',&
recl=recLen,form='unformatted',status='unknown')
ENDIF
CLOSE(59)
initBroydenHistory = l_exist
END FUNCTION initBroydenHistory
END MODULE m_broyd_io
......@@ -35,6 +35,7 @@ MODULE m_cdn_io
PUBLIC readStepfunction, writeStepfunction
PUBLIC setStartingDensity, readPrevEFermi, deleteDensities
PUBLIC storeStructureIfNew
PUBLIC getIOMode
PUBLIC CDN_INPUT_DEN_const, CDN_OUTPUT_DEN_const
PUBLIC CDN_ARCHIVE_TYPE_CDN1_const, CDN_ARCHIVE_TYPE_NOCO_const
PUBLIC CDN_ARCHIVE_TYPE_CDN_const
......@@ -76,7 +77,7 @@ MODULE m_cdn_io
CHARACTER(LEN=19) :: timeStampString
CHARACTER(LEN=15) :: distanceString
CALL getMode(mode)
CALL getIOMode(mode)
WRITE(*,*) 'Available densities info:'
WRITE(*,*)
......@@ -168,7 +169,7 @@ MODULE m_cdn_io
fermiEnergy = 0.0
l_qfix = .FALSE.
CALL getMode(mode)
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
#ifdef CPP_HDF
......@@ -385,7 +386,7 @@ MODULE m_cdn_io
CHARACTER(LEN=10) :: timeString
CHARACTER(LEN=10) :: zone
CALL getMode(mode)
CALL getIOMode(mode)
CALL DATE_AND_TIME(dateString,timeString,zone)
READ(dateString,'(i8)') date
READ(timeString,'(i6)') time
......@@ -648,7 +649,7 @@ MODULE m_cdn_io
LOGICAL :: l_qfix, l_exist
CHARACTER(LEN=30) :: archiveName
CALL getMode(mode)
CALL getIOMode(mode)
eFermiPrev = 0.0
l_error = .FALSE.
......@@ -702,7 +703,7 @@ MODULE m_cdn_io
INTEGER :: currentStructureIndex,currentStepfunctionIndex
INTEGER :: readDensityIndex, lastDensityIndex
CALL getMode(mode)
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
#ifdef CPP_HDF
......@@ -773,7 +774,7 @@ MODULE m_cdn_io
INTEGER :: currentStructureIndex,currentStepfunctionIndex
INTEGER :: readDensityIndex, lastDensityIndex
CALL getMode(mode)
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
#ifdef CPP_HDF
......@@ -824,7 +825,7 @@ MODULE m_cdn_io
INTEGER(HID_T) :: fileID
#endif
CALL getMode(mode)
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
#ifdef CPP_HDF
......@@ -881,7 +882,7 @@ MODULE m_cdn_io
izmax = 0
igz = 0
CALL getMode(mode)
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
#ifdef CPP_HDF
......@@ -953,7 +954,7 @@ MODULE m_cdn_io
izmax = 0
igz = 0
CALL getMode(mode)
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
......@@ -1058,7 +1059,7 @@ MODULE m_cdn_io
ifftd=size(stars%ufft)
CALL getMode(mode)
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
#ifdef CPP_HDF
......@@ -1117,7 +1118,7 @@ MODULE m_cdn_io
ioStatus = 0
ifftd = 27*stars%mx1*stars%mx2*stars%mx3
CALL getMode(mode)
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
......@@ -1208,7 +1209,7 @@ MODULE m_cdn_io
WRITE(*,'(a,i0,a)') 'Using density ', sdIndex, ' as starting density.'
CALL getMode(mode)
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
#ifdef CPP_HDF
......@@ -1266,7 +1267,7 @@ MODULE m_cdn_io
startNumber = -1
endNumber = -1
IF (TRIM(ADJUSTL(ddString)).EQ.'allbutlast') THEN
CALL getMode(mode)
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
......@@ -1300,7 +1301,7 @@ MODULE m_cdn_io
END IF
END IF
CALL getMode(mode)
CALL getIOMode(mode)
IF(mode.EQ.CDN_HDF5_MODE) THEN
INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
......@@ -1348,7 +1349,7 @@ MODULE m_cdn_io
END SUBROUTINE deleteDensities
SUBROUTINE getMode(mode)
SUBROUTINE getIOMode(mode)
INTEGER, INTENT(OUT) :: mode
mode = CDN_DIRECT_MODE
......@@ -1363,7 +1364,7 @@ MODULE m_cdn_io
! WRITE(*,*) 'Falling back to direct access.'
#endif
END IF
END SUBROUTINE getMode
END SUBROUTINE getIOMode
LOGICAL FUNCTION isDensityFilePresent(archiveType)
......@@ -1378,7 +1379,7 @@ MODULE m_cdn_io
INTEGER(HID_T) :: fileID
#endif
CALL getMode(mode)
CALL getIOMode(mode)
IF (mode.EQ.CDN_HDF5_MODE) THEN
INQUIRE(FILE='cdn.hdf',EXIST=l_exist)
......@@ -1425,7 +1426,7 @@ MODULE m_cdn_io
LOGICAL :: l_exist
INTEGER :: mode
CALL getMode(mode)
CALL getIOMode(mode)
IF (mode.EQ.CDN_HDF5_MODE) THEN
#ifdef CPP_HDF
......
......@@ -23,6 +23,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
USE m_juDFT
USE m_constants
USE m_cdn_io
USE m_broyd_io
USE m_brysh1
USE m_stmix
USE m_broyden
......@@ -51,7 +52,7 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
!Local Scalars
REAL fix,intfac,vacfac
INTEGER i,imap,js,mit,irecl
INTEGER i,imap,js
INTEGER mmap,mmaph,nmaph,nmap,mapmt,mapvac,mapvac2
INTEGER iofl,n_u_keep
LOGICAL l_exist,l_ldaU, l_densityMatrixPresent
......@@ -114,26 +115,18 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
! LDA+U (end)
ALLOCATE (sm(mmap),fsm(mmap))
INQUIRE (file='broyd.'//CHAR(input%imix+48),exist=l_exist)
DO i = 1,6
dist(i) = 0.0
END DO
dist(:) = 0.0
!determine type of mixing:
!imix=0:straight, imix=o broyden first, imix=5:broyden second
!imix=:generalozed anderson mixing
mit = 0
IF (input%imix.EQ.0) THEN
WRITE (16,FMT='(a,2f10.5)') 'STRAIGHT MIXING',input%alpha
ELSE IF (input%imix.EQ.3) THEN
IF ( .NOT.l_exist) mit = 1
WRITE (16,FMT='(a,f10.5)') 'BROYDEN FIRST MIXING',input%alpha
ELSE IF (input%imix.EQ.5) THEN
IF (.NOT.l_exist) mit = 1
WRITE (16,FMT='(a,f10.5)') 'BROYDEN SECOND MIXING',input%alpha
ELSE IF (input%imix.EQ.7) THEN
IF (.NOT.l_exist) mit = 1
WRITE (16,FMT='(a,f10.5)') 'ANDERSON GENERALIZED',input%alpha
ELSE
CALL juDFT_error("mix: input%imix =/= 0,3,5,7 ",calledby ="mix")
......@@ -156,26 +149,12 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
!store the difference fsm - sm in fsm
fsm(:nmap) = fsm(:nmap) - sm(:nmap)
! open files for broyden
irecl=(nmap+1)*8
IF (input%imix.GE.3) THEN
IF (hybrid%l_calhf) THEN
OPEN (57,file='hf_broyd',form='unformatted',status='unknown')
OPEN (59,file='hf_broyd.'//CHAR(input%imix+48),access='direct',&
recl=irecl,form='unformatted',status='unknown')
ELSE
OPEN (57,file='broyd',form='unformatted',status='unknown')
OPEN (59,file='broyd.'//CHAR(input%imix+48),access='direct',&
recl=irecl,form='unformatted',status='unknown')
ENDIF
END IF
!mixing of the densities
IF (input%imix.EQ.0) THEN
CALL stmix(atoms,input,noco, nmap,nmaph,fsm, sm)
ELSE
CALL broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
mmap,nmaph,mapmt,mapvac2,nmap,fsm,mit,sm)
hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fsm,sm)
END IF
!initiatlize mixed density and extract it with brysh2 call
......@@ -289,11 +268,6 @@ SUBROUTINE mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,&
inDen%iter = inDen%iter + 1
IF (input%imix.GT.0) THEN
CLOSE (57)
CLOSE (59)
END IF
7900 FORMAT (/,'----> distance of charge densities for spin ',i2,' it=',i5,':',f13.6,' me/bohr**3')
7901 FORMAT (/,'----> HF distance of charge densities for spin ',i2,' it=',i5,':',f13.6,' me/bohr**3')
8000 FORMAT (/,'----> distance of charge densities for it=',i5,':', f13.6,' me/bohr**3')
......
......@@ -12,12 +12,13 @@ MODULE m_broyden
!################################################################
CONTAINS
SUBROUTINE broyden(cell,stars,atoms,vacuum,sphhar,input,noco,oneD,sym,&
mmap,nmaph,mapmt,mapvac2,nmap,fm,mit,sm,lpot)
hybrid,mmap,nmaph,mapmt,mapvac2,nmap,fm,sm,lpot)
#include"cpp_double.h"
USE m_metric
USE m_types
USE m_broyd_io
IMPLICIT NONE
......@@ -30,11 +31,11 @@ CONTAINS
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_hybrid),INTENT(IN) :: hybrid
! Scalar Arguments
INTEGER, INTENT (IN) :: mmap,nmap
INTEGER, INTENT (IN) :: mapmt,mapvac2
INTEGER, INTENT (INOUT) :: mit
INTEGER, INTENT (IN) :: mapmt,mapvac2
LOGICAL,OPTIONAL,INTENT(IN) :: lpot
! Array Arguments
......@@ -42,9 +43,9 @@ CONTAINS
REAL, INTENT (INOUT) :: sm(nmap)
! Local Scalars
INTEGER :: i,it,k,nit,npos,iread,nmaph
INTEGER :: i,it,k,nit,npos,iread,nmaph, mit
REAL :: bm,dfivi,fmvm,smnorm,vmnorm,alphan
LOGICAL :: l_pot
LOGICAL :: l_pot, l_exist
REAL, PARAMETER :: one=1.0, zero=0.0
! Local Arrays
......@@ -74,12 +75,15 @@ CONTAINS
vm = 0.0
am = 0.0
mit = 0
l_exist = initBroydenHistory(input,hybrid,nmap) ! returns true if there already exists a Broyden history
IF(.NOT.l_exist) mit = 1
IF (mit.NE.1) THEN
! load input charge density (sm1) and difference of
! in and out charge densities (fm1) from previous iteration (m-1)
REWIND 57
READ (57) mit, alphan, (fm1(i),i=1,nmap), (sm1(i),i=1,nmap)
CALL readLastIterInAndDiffDen(hybrid,nmap,mit,alphan,sm1(:nmap),fm1(:nmap))
IF (ABS(input%alpha-alphan) > 0.0001) THEN
WRITE (6,*) 'mixing parameter has been changed; reset'
WRITE (6,*) 'broyden algorithm or set alpha to',alphan
......@@ -93,10 +97,9 @@ CONTAINS
END IF
! save F_m and rho_m for next iteration
REWIND 57
nit = mit +1
IF (nit > input%maxiter+1) nit = 1
WRITE (57) nit,input%alpha,fm,sm
CALL writeLastIterInAndDiffDen(hybrid,nmap,nit,input%alpha,sm,fm)
IF (mit.EQ.1) THEN
! update for rho for mit=1 is straight mixing
......@@ -110,8 +113,9 @@ CONTAINS
END DO
iread = MIN(mit-1,input%maxiter+1)
DO it = 2, iread
READ (59,rec=(it-1)*2-1) (ui(i),i=1,nmap)
READ (59,rec=(it-1)*2) (vi(i),i=1,nmap),dfivi
CALL readUVec(input,hybrid,nmap,it-mit,mit,ui)
CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)
am(it) = CPP_BLAS_sdot(nmap,vi,1,fm1,1)
! calculate um(:) = -am(it)*ui(:) + um
CALL CPP_BLAS_saxpy(nmap,-am(it),ui,1,um,1)
......@@ -133,8 +137,8 @@ CONTAINS
! generate vm = alpha*sm1 - \sum <ui|w|sm1> vi
vm(:) = input%alpha * fm1(:)
DO it = 2,iread
READ (59,rec=2*(it-1)-1) (ui(i),i=1,nmap)
READ (59,rec=2*(it-1)) (vi(i),i=1,nmap), dfivi
CALL readUVec(input,hybrid,nmap,it-mit,mit,ui)
CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)
bm = CPP_BLAS_sdot(nmap,ui,1,fm1,1)
! calculate vm(:) = -bm*vi(:) + vm
CALL CPP_BLAS_saxpy(nmap,-bm,vi,1,vm,1)
......@@ -172,7 +176,7 @@ CONTAINS
mmap,nmaph,mapmt,mapvac2,fm1,vm,l_pot)
DO it = 2,iread
READ (59,rec=2*(it-1)) (vi(i),i=1,nmap), dfivi
CALL readVVec(input,hybrid,nmap,it-mit,mit,dfivi,vi)
! calculate vm(:) = -am(it)*dfivi*vi(:) + vm
CALL CPP_BLAS_saxpy(nmap,-am(it)*dfivi,vi,1,vm,1)
END DO
......@@ -193,8 +197,8 @@ CONTAINS
npos=mit-1
IF (mit.GT.input%maxiter+1) npos = MOD(mit-2,input%maxiter)+1
WRITE (59,rec=2*npos-1) (um(i),i=1,nmap)
WRITE (59,rec=2*npos) (vm(i),i=1,nmap), dfivi
CALL writeUVec(input,hybrid,nmap,mit,um)
CALL writeVVec(input,hybrid,nmap,mit,dfivi,vm)
! update rho(m+1)
! calculate <fm|w|vm>
......@@ -203,7 +207,6 @@ CONTAINS
CALL CPP_BLAS_saxpy(nmap,one-fmvm,um,1,sm,1)
END IF
mit = mit + 1
DEALLOCATE (fm1,sm1,ui,um,vi,vm,am)
END SUBROUTINE broyden
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment