Commit 2b4e4632 authored by Daniel Wortmann's avatar Daniel Wortmann

Finalized basic Pulay implementation

parent 41570f84
...@@ -331,6 +331,12 @@ CONTAINS ...@@ -331,6 +331,12 @@ CONTAINS
input%imix = 7 input%imix = 7
CASE ("Pulay") CASE ("Pulay")
input%imix = 9 input%imix = 9
CASE ("pPulay")
input%imix = 11
CASE ("rPulay")
input%imix = 13
CASE ("aPulay")
input%imix = 15
CASE DEFAULT CASE DEFAULT
STOP 'Error: unknown mixing scheme selected!' STOP 'Error: unknown mixing scheme selected!'
END SELECT END SELECT
......
...@@ -880,6 +880,9 @@ ...@@ -880,6 +880,9 @@
<xsd:enumeration value="Broyden2"/> <xsd:enumeration value="Broyden2"/>
<xsd:enumeration value="Anderson"/> <xsd:enumeration value="Anderson"/>
<xsd:enumeration value="Pulay"/> <xsd:enumeration value="Pulay"/>
<xsd:enumeration value="rPulay"/>
<xsd:enumeration value="pPulay"/>
<xsd:enumeration value="aPulay"/>
</xsd:restriction> </xsd:restriction>
</xsd:simpleType> </xsd:simpleType>
......
...@@ -411,11 +411,9 @@ CONTAINS ...@@ -411,11 +411,9 @@ CONTAINS
field2 = field field2 = field
! mix input and output densities ! mix input and output densities
CALL timestart("mixing")
CALL mix_charge(field2,DIMENSION,mpi,(iter==input%itmax.OR.judft_was_argument("-mix_io")),& CALL mix_charge(field2,DIMENSION,mpi,(iter==input%itmax.OR.judft_was_argument("-mix_io")),&
stars,atoms,sphhar,vacuum,input,& stars,atoms,sphhar,vacuum,input,&
sym,cell,noco,oneD,archiveType,inDen,outDen,results) sym,cell,noco,oneD,archiveType,inDen,outDen,results)
CALL timestop("mixing")
IF(mpi%irank == 0) THEN IF(mpi%irank == 0) THEN
WRITE (6,FMT=8130) iter WRITE (6,FMT=8130) iter
......
...@@ -29,6 +29,7 @@ contains ...@@ -29,6 +29,7 @@ contains
use m_umix use m_umix
USE m_kerker USE m_kerker
use m_pulay use m_pulay
use m_a_pulay
use m_types_mixvector use m_types_mixvector
USE m_distance USE m_distance
use m_mixing_history use m_mixing_history
...@@ -60,35 +61,7 @@ contains ...@@ -60,35 +61,7 @@ contains
INTEGER :: it,maxiter INTEGER :: it,maxiter
MPI0_a: IF( mpi%irank == 0 ) THEN CALL timestart("Charge Density Mixing")
!determine type of mixing:
!imix=0:straight, imix=o broyden first, imix=5:broyden second
!imix=:generalozed anderson mixing
select case( input%imix )
case( 0 )
WRITE( 6, fmt='(a,2f10.5)' ) 'STRAIGHT MIXING',input%alpha
IF (input%jspins.EQ.1) WRITE (6,FMT='(a,2f10.5)')&
& 'charge density mixing parameter:',input%alpha
IF (input%jspins.EQ.2) WRITE (6,FMT='(a,2f10.5)')&
& 'spin density mixing parameter:',input%alpha*input%spinf
case( 3 )
write( 6, fmt='(a,f10.5)' ) 'BROYDEN FIRST MIXING',input%alpha
case( 5 )
write( 6, fmt='(a,f10.5)' ) 'BROYDEN SECOND MIXING',input%alpha
case( 7 )
write( 6, fmt='(a,f10.5)' ) 'ANDERSON GENERALIZED',input%alpha
case ( 9 )
write( 6, fmt='(a,f10.5)' ) 'PULAY MIXING',input%alpha
case default
call juDFT_error( "mix: input%imix =/= 0,3,5,7,9 ", calledby ="mix" )
end select
if ( input%jspins == 2 .and. input%imix /= 0 ) then
write( 6, '(''WARNING : for QUASI-NEWTON METHODS SPINF=1'')' )
end if
ENDIF MPI0_a
l_densitymatrix=.FALSE. l_densitymatrix=.FALSE.
IF (atoms%n_u>0) THEN IF (atoms%n_u>0) THEN
l_densitymatrix=.NOT.input%ldaulinmix l_densitymatrix=.NOT.input%ldaulinmix
...@@ -99,6 +72,8 @@ contains ...@@ -99,6 +72,8 @@ contains
if (mpi%irank.ne.0) inden%mmpmat=0.0 if (mpi%irank.ne.0) inden%mmpmat=0.0
ENDIF ENDIF
ENDIF ENDIF
CALL timestart("Reading of distances")
CALL mixvector_init(mpi%mpi_comm,l_densitymatrix,oneD,input,vacuum,noco,sym,stars,cell,sphhar,atoms) CALL mixvector_init(mpi%mpi_comm,l_densitymatrix,oneD,input,vacuum,noco,sym,stars,cell,sphhar,atoms)
CALL mixing_history_open(mpi,input%maxiter) CALL mixing_history_open(mpi,input%maxiter)
...@@ -107,24 +82,58 @@ contains ...@@ -107,24 +82,58 @@ contains
CALL mixing_history(input%imix,maxiter,inden,outden,sm,fsm,it) CALL mixing_history(input%imix,maxiter,inden,outden,sm,fsm,it)
CALL distance(mpi%irank,cell%vol,input%jspins,fsm(it),inDen,outDen,results,fsm_Mag) CALL distance(mpi%irank,cell%vol,input%jspins,fsm(it),inDen,outDen,results,fsm_Mag)
CALL timestop("Reading of distances")
! KERKER PRECONDITIONER ! KERKER PRECONDITIONER
IF( input%preconditioning_param /= 0 ) THEN IF( input%preconditioning_param /= 0 ) THEN
CALL timestart("Preconditioner")
CALL kerker(field, DIMENSION, mpi, & CALL kerker(field, DIMENSION, mpi, &
stars, atoms, sphhar, vacuum, input, sym, cell, noco, & stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, inDen, outDen, fsm(it) ) oneD, inDen, outDen, fsm(it) )
!Store modified density in history !Store modified density in history
CALL mixing_history_store(fsm(it)) CALL mixing_history_store(fsm(it))
CALL timestop("Preconditioner")
END IF END IF
CALL timestart("Mixing")
!mixing of the densities !mixing of the densities
if(input%imix==0.or.it==1) CALL stmix(atoms,input,noco,fsm(it),fsm_mag,sm(it)) SELECT CASE(input%imix)
!if(it>1.and.input%imix==9) CALL pulay(input%alpha,fsm,sm) CASE(0)
if(it>1.and.(input%imix==3.or.input%imix==5.or.input%imix==7)) Call broyden(input%alpha,fsm,sm) IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,f10.5)' ) &
!PRINT *,"ATTENTION Broyden replaced by Pulay" 'STRAIGHT MIXING: alpha=',input%alpha," spin-mixing=",MERGE(input%alpha*input%spinf,0.,input%jspins>1)
IF (it>1.and.input%imix==9) THEN CALL stmix(atoms,input,noco,fsm(it),fsm_mag,sm(it))
CALL pulay(input%alpha,fsm,sm) CASE(3,5)
if (it==input%maxiter) call mixing_history_limit(1) CALL judft_error("Broyden 1/2 method not implemented")
endif CASE(7)
IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,i0)' ) &
'GENERALIZED ANDERSON MIXING: alpha=',input%alpha," History-length=",it-1
Call broyden(input%alpha,fsm,sm)
CASE(9)
IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,i0,a,i0)' ) &
'PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
CALL pulay(input%alpha,fsm,sm,0)
CASE(11)
IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,i0,a,i0)' ) &
'PERIODIC PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
CALL pulay(input%alpha,fsm,sm,input%maxiter)
CASE(13)
IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,i0,a,i0)' ) &
'RESTARTED PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
CALL pulay(input%alpha,fsm,sm,0)
IF (it==input%maxiter) CALL mixing_history_limit(0) !Restarting Pulay
CASE(15)
IF (mpi%irank==0) WRITE( 6, fmt='(a,f10.5,a,i0,a,i0)' ) &
'ADAPTED PULAY MIXING: alpha=',input%alpha," History-length=",it-1,"/",input%maxiter
CALL a_pulay(input%alpha,fsm,sm)
IF (input%maxiter<4) CALL judft_error("History length too short for aPulay")
IF (it==4) CALL mixing_history_limit(0)
CASE DEFAULT
CALL judft_error("Unknown Mixing schema")
END SELECT
CALL timestop("Mixing")
CALL timestart("Postprocessing")
!extracte mixed density !extracte mixed density
inDen%pw=0.0;inDen%mt=0.0 inDen%pw=0.0;inDen%mt=0.0
IF (ALLOCATED(inDen%vacz)) inden%vacz=0.0 IF (ALLOCATED(inDen%vacz)) inden%vacz=0.0
...@@ -171,6 +180,9 @@ contains ...@@ -171,6 +180,9 @@ contains
IF (l_writehistory.AND.input%imix.NE.0) CALL mixing_history_close(mpi) IF (l_writehistory.AND.input%imix.NE.0) CALL mixing_history_close(mpi)
CALL timestop("Postprocessing")
CALL timestop("Charge Density Mixing")
END SUBROUTINE mix_charge END SUBROUTINE mix_charge
END MODULE m_mix END MODULE m_mix
...@@ -8,6 +8,7 @@ mix/broyden_history.F90 ...@@ -8,6 +8,7 @@ mix/broyden_history.F90
mix/mixing_history.F90 mix/mixing_history.F90
mix/distance.F90 mix/distance.F90
mix/pulay.F90 mix/pulay.F90
mix/a_pulay.F90
#mix/broyden2.F90 #mix/broyden2.F90
#mix/brysh1.f90 #mix/brysh1.f90
#mix/brysh2.f90 #mix/brysh2.f90
......
...@@ -28,6 +28,11 @@ CONTAINS ...@@ -28,6 +28,11 @@ CONTAINS
TYPE(t_mixvector),allocatable :: u_store(:),v_store(:) TYPE(t_mixvector),allocatable :: u_store(:),v_store(:)
hlen=size(fm) hlen=size(fm)
IF (hlen<2) THEN !Do a simple mixing step
sm(hlen)=sm(hlen)+alpha*fm(hlen)
RETURN
ENDIF
ALLOCATE(u_store(hlen-2),v_store(hlen-2)) ALLOCATE(u_store(hlen-2),v_store(hlen-2))
do it=1,hlen-2 do it=1,hlen-2
call u_store(it)%alloc() call u_store(it)%alloc()
......
...@@ -80,7 +80,7 @@ CONTAINS ...@@ -80,7 +80,7 @@ CONTAINS
if (.not.allocated(sm_store)) THEN if (.not.allocated(sm_store)) THEN
allocate(sm_store(maxiter),fsm_store(maxiter)) allocate(sm_store(maxiter),fsm_store(maxiter))
endif endif
IF (iter_stored+1==maxiter.AND.imix.NE.9) iter_stored=0 !This is a broyden method which has to IF (iter_stored+1==maxiter.AND.imix<8) iter_stored=0 !This is a broyden method which has to
!be reset as soon as maxiter is reached !be reset as soon as maxiter is reached
it=iter_stored+1 it=iter_stored+1
allocate(sm(it),fsm(it)) allocate(sm(it),fsm(it))
......
...@@ -2,30 +2,36 @@ MODULE m_pulay ...@@ -2,30 +2,36 @@ MODULE m_pulay
USE m_juDFT USE m_juDFT
CONTAINS CONTAINS
SUBROUTINE pulay(alpha,fm,sm) SUBROUTINE pulay(alpha,fm,sm,simple_steps)
USE m_types_mat USE m_types_mat
USE m_types_mixvector USE m_types_mixvector
USE m_mixing_history
IMPLICIT NONE IMPLICIT NONE
REAL,INTENT(IN) :: alpha
TYPE(t_mixvector),INTENT(IN) :: fm(:) TYPE(t_mixvector),INTENT(IN) :: fm(:)
TYPE(t_mixvector),INTENT(INOUT) :: sm(:) TYPE(t_mixvector),INTENT(INOUT) :: sm(:)
INTEGER,INTENT(IN) :: simple_steps
REAL,INTENT(IN) :: alpha
! Locals ! Locals
INTEGER :: h_len,n,nn INTEGER :: h_len,n,nn
REAL,ALLOCATABLE :: b(:) REAL,ALLOCATABLE :: b(:)
TYPE(t_mat) :: a TYPE(t_mat) :: a
TYPE(t_mixvector),ALLOCATABLE :: df(:),ds(:) TYPE(t_mixvector),ALLOCATABLE :: df(:),ds(:)
TYPE(t_mixvector) ::mdf
h_len=SIZE(fm)-1 h_len=SIZE(fm)-1
sm(h_len+1)=sm(h_len+1)+alpha*fm(h_len+1) !Simple mixing IF (h_len==0) THEN
sm(h_len+1)=sm(h_len+1)+alpha*fm(h_len+1) !Simple mixing
PRINT *,"len:",h_len RETURN !No history present
DO n=1,h_len ENDIF
PRINT *,"D:",sm(n).dot.sm(n),fm(n).dot.fm(n) IF (simple_steps>0) THEN
ENDDO IF (MOD(h_len,simple_steps).NE.0) THEN
IF (h_len==1) RETURN !Simple mixing step of periodic Pulay
sm(h_len+1)=sm(h_len+1)+alpha*fm(h_len+1)
RETURN
ENDIF
ENDIF
CALL a%alloc(.TRUE.,h_len,h_len) CALL a%alloc(.TRUE.,h_len,h_len)
ALLOCATE(df(h_len),ds(h_len),b(h_len)) ALLOCATE(df(h_len),ds(h_len),b(h_len))
...@@ -33,25 +39,22 @@ CONTAINS ...@@ -33,25 +39,22 @@ CONTAINS
DO n=1,h_len DO n=1,h_len
df(n)=fm(n+1)-fm(n) df(n)=fm(n+1)-fm(n)
ds(n)=sm(n+1)-sm(n) ds(n)=sm(n+1)-sm(n)
b(n)=df(n).dot.fm(h_len+1) mdf=df(n).apply_metric()
b(n)=mdf.dot.fm(h_len+1)
DO nn=1,n DO nn=1,n
a%data_r(n,nn)=df(n).dot.df(nn) a%data_r(n,nn)=mdf.dot.df(nn)
a%data_r(nn,n)=a%data_r(n,nn) a%data_r(nn,n)=a%data_r(n,nn)
ENDDO ENDDO
ENDDO ENDDO
call a%inverse() call a%inverse()
PRINT *,"b0:",b
b=MATMUL(a%data_r,b) b=MATMUL(a%data_r,b)
sm(h_len+1)=sm(h_len+1)+alpha*fm(h_len+1) !Simple mixing
PRINT *,"b1:",b
DO n=1,h_len DO n=1,h_len
sm(h_len+1)=sm(h_len+1)-b(n)*(ds(n)+alpha*df(n)) sm(h_len+1)=sm(h_len+1)-b(n)*(ds(n)+alpha*df(n))
ENDDO ENDDO
IF (simple_steps>0) CALL mixing_history_limit(0)
END SUBROUTINE pulay END SUBROUTINE pulay
END MODULE m_pulay END MODULE m_pulay
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