Commit 87d45ae0 authored by Daniel Wortmann's avatar Daniel Wortmann

Modified experimental parallel pulay method.

parent c645cc1b
MODULE m_tlmplm_cholesky
use m_judft
IMPLICIT NONE
!*********************************************************************
! sets up the local Hamiltonian, i.e. the Hamiltonian in the
......
......@@ -6,6 +6,7 @@
MODULE m_io_matrix
USE m_types
USE m_judft
IMPLICIT NONE
private
INTEGER:: mode=1
......
......@@ -125,8 +125,6 @@ contains
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
......
......@@ -2,72 +2,104 @@ MODULE m_a_pulay
USE m_juDFT
PRIVATE
REAL :: distance(3)
REAL :: alpha(4)=(/0.0,0.0,0.0,0.0/)
PUBLIC a_pulay
CONTAINS
SUBROUTINE a_pulay(alpha_in,fm,sm)
SUBROUTINE a_pulay(alpha,fm,sm)
USE m_pulay
USE m_types_mat
USE m_types_mixvector
USE m_mixing_history
IMPLICIT NONE
REAL,INTENT(IN) :: alpha_in
REAL,INTENT(IN) :: alpha
TYPE(t_mixvector),INTENT(IN) :: fm(:)
TYPE(t_mixvector),INTENT(INOUT) :: sm(:)
TYPE(t_mixvector) :: mdf
TYPE(t_mat) :: a
INTEGER :: n,nn
IF (alpha(1)==0.0) alpha=(/1.,2.,0.5,1./)*alpha_in
REAL :: fac(2)=(/1.0,0.5/)
INTEGER,parameter :: parallel_images=2
INTEGER,PARAMETER :: maxiter=4,mode=2
INTEGER :: hlen,image,i,ii,local_length
INTEGER :: local_hist(maxiter*parallel_images)
SELECT CASE(SIZE(sm))
hlen=SIZE(sm) !Number of results in history
image=MOD(hlen-1,parallel_images)+1
PRINT *,"HI:",hlen,image
IF (hlen<parallel_images) THEN
!create inital parallel images by simple mixing of first result
sm(hlen)=sm(1)+(alpha*fac(image))*fm(1)
RETURN
ENDIF
SELECT CASE(mode)
CASE(1)
mdf=fm(1)%apply_metric()
sm(1)=sm(1)+alpha(1)*fm(1) !Simple mixing
CASE(2)
mdf=fm(2)%apply_metric()
distance(1)=fm(2).dot.mdf
sm(2)=sm(1)+alpha(2)*fm(1) !Simple mixing with double alpha
CASE(3)
mdf=fm(3)%apply_metric()
distance(2)=fm(3).dot.mdf
sm(3)=sm(1)+alpha(3)*fm(1) !Simple mixing with half alpha
CASE(4)
!Find best distance
mdf=fm(4)%apply_metric()
distance(3)=fm(4).dot.mdf
!IF (distance(2)==MINVAL(distance)) alpha(4)=2*alpha(1)
!IF (distance(3)==MINVAL(distance)) alpha(4)=0.5*alpha(1)
!WRITE(*,'(a,3f9.4)'),"A-Pulay:",distance
!WRITE(*,*) alpha(1)
CALL pulay(alpha_in,fm,sm,0)
!!$ CALL a%alloc(.TRUE.,5,5)
!!$ DO n=1,4
!!$ mdf=fm(n)%apply_metric()
!!$ DO nn=1,n
!!$ a%data_r(n,nn)=mdf.dot.fm(nn)
!!$ a%data_r(nn,n)=a%data_r(n,nn)
!!$ ENDDO
!!$ ENDDO
!!$ a%data_r(:,5)=1.0
!!$ a%data_r(5,:)=1.0
!!$ a%data_r(5,5)=.0
!!$ PRINT *,"a:",a%data_r(1,1)
!!$ PRINT *,"a:",a%data_r(2,2)
!!$ PRINT *,"a:",a%data_r(3,3)
!!$ PRINT *,"a:",a%data_r(4,4)
!!$
!!$ CALL a%inverse()
!!$ PRINT *,"p:",a%data_r(:4,5)
!!$ mdf=0.0*mdf
!!$ DO n=1,4
!!$ mdf=mdf+a%data_r(n,5)*(sm(n)+alpha(n)*fm(n))
!!$ ENDDO
!!$ sm(4)=mdf
local_length=0
DO i=image,hlen,parallel_images
local_length=local_length+1
local_hist(local_length)=i
ENDDO
IF (local_length>maxiter) THEN !recreate to enable cross-image mixing
local_length=hlen+1-parallel_images
local_hist(1)=image
local_hist(2:local_length)=(/(i,i=parallel_images+1,hlen)/)
END IF
sm(hlen)=simple_pulay(alpha*fac(image),fm(local_hist(:local_length)),sm(local_hist(:local_length)))
IF (hlen==parallel_images*(maxiter+1)) CALL mixing_history_limit(parallel_images)
case(2)
local_length=hlen-image
local_hist(1:local_length)=(/(i,i=1,local_length)/)
local_length=local_length+1
local_hist(local_length)=hlen
sm(hlen)=simple_pulay(alpha*fac(image),fm(local_hist(:local_length)),sm(local_hist(:local_length)))
IF (hlen==parallel_images*(maxiter)) CALL mixing_history_limit(parallel_images)
END SELECT
PRINT *,"P:",local_hist(:local_length)
END SUBROUTINE a_pulay
FUNCTION simple_pulay(alpha,fm,sm)RESULT(sm_out)
USE m_types_mixvector
IMPLICIT NONE
REAL,INTENT(IN) :: alpha
TYPE(t_mixvector),INTENT(IN) :: fm(:)
TYPE(t_mixvector),INTENT(IN) :: sm(:)
TYPE(t_mixvector) :: sm_out
TYPE(t_mixvector) :: mdf
REAL,ALLOCATABLE :: a(:,:),b(:),work(:)
INTEGER,ALLOCATABLE :: ipiv(:)
INTEGER :: n,nn,info,lwork,hlen
hlen=SIZE(sm)
IF (hlen==1) THEN
sm_out=sm(1)+alpha*fm(1)
RETURN
ENDIF
ALLOCATE(a(hlen+1,hlen+1),b(hlen+1),ipiv(hlen+1),work((hlen+1)**2))
DO n=1,hlen
mdf=fm(n)%apply_metric()
DO nn=1,n
a(n,nn)=mdf.dot.fm(nn)
a(nn,n)=a(n,nn)
ENDDO
ENDDO
a(:,hlen+1)=1.0
a(hlen+1,:)=1.0
a(hlen+1,hlen+1)=0.0
b=0.0;b(hlen+1)=1.
CALL DSYSV( 'Upper', hlen+1, 1, a, SIZE(a,1), ipiv, b, SIZE(b,1), work, SIZE(work), INFO )
sm_out=0.0*mdf
DO n=1,hlen
sm_out=sm_out+b(n)*(sm(n)+alpha*fm(n))
ENDDO
END FUNCTION SIMPLE_PULAY
END MODULE m_a_pulay
......@@ -4,6 +4,7 @@
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_mpi_bc_xcpot
use m_judft
CONTAINS
SUBROUTINE mpi_bc_xcpot(xcpot,mpi)
USE m_types
......
MODULE m_types_mat
USE m_judft
IMPLICIT NONE
PRIVATE
!<This is the basic type to store and manipulate real/complex rank-2 matrices
!!
!! In its simple implementation here, it contains a fields for the matrix-size and
......@@ -20,6 +20,7 @@ MODULE m_types_mat
PROCEDURE :: transpose=>t_mat_transpose !> transpose the matrix
PROCEDURE :: from_packed=>t_mat_from_packed !> initialized from a packed-storage matrix
PROCEDURE :: inverse =>t_mat_inverse !> invert the matrix
PROCEDURE :: linear_problem => t_mat_lproblem !> Solve linear equation
PROCEDURE :: to_packed=>t_mat_to_packed !> convert to packed-storage matrix
PROCEDURE :: clear => t_mat_clear !> set data arrays to zero
PROCEDURE :: copy => t_mat_copy !> copy into another t_mat (overloaded for t_mpimat)
......@@ -30,9 +31,41 @@ MODULE m_types_mat
PROCEDURE :: free => t_mat_free !> dealloc the data (overloaded for t_mpimat)
PROCEDURE :: add_transpose => t_mat_add_transpose!> add the tranpose/Hermitian conjg. without the diagonal (overloaded for t_mpimat)
END type t_mat
PUBLIC t_mat
CONTAINS
SUBROUTINE t_mat_lproblem(mat,vec)
IMPLICIT NONE
CLASS(t_mat),INTENT(IN) :: mat
TYPE(t_mat),INTENT(INOUT) :: vec
INTEGER:: lwork,info
REAL,ALLOCATABLE:: work(:)
INTEGER,allocatable::ipiv(:)
IF ((mat%l_real.NE.vec%l_real).OR.(mat%matsize1.NE.mat%matsize2).OR.(mat%matsize1.NE.vec%matsize1)) &
CALL judft_error("Invalid matices in t_mat_lproblem")
IF (mat%l_real) THEN
IF (ALL(ABS(mat%data_r-TRANSPOSE(mat%data_r))<1E-8)) THEN
!Matrix is symmetric
CALL DPOSV( 'Upper', mat%matsize1, vec%matsize2, mat%data_r, mat%matsize1, vec%data_r, vec%matsize1, INFO )
IF (INFO>0) THEN
!Matrix was not positive definite
lwork=-1;ALLOCATE(work(1))
CALL DSYSV( 'Upper', mat%matsize1, vec%matsize2, mat%data_r, mat%matsize1, IPIV, vec%data_r, vec%matsize1, WORK, LWORK,INFO )
lwork=INT(work(1))
DEALLOCATE(work);ALLOCATE(ipiv(mat%matsize1),work(lwork))
CALL DSYSV( 'Upper', mat%matsize1, vec%matsize2, mat%data_r, mat%matsize1, IPIV, vec%data_r, vec%matsize1, WORK, LWORK,INFO )
IF (info.NE.0) CALL judft_error("Could not solve linear equation, matrix singular")
END IF
ELSE
CALL judft_error("TODO: mode not implemented in t_mat_lproblem")
END IF
ELSE
CALL judft_error("TODO: mode not implemented in t_mat_lproblem")
ENDIF
END SUBROUTINE t_mat_lproblem
SUBROUTINE t_mat_free(mat)
CLASS(t_mat),INTENT(INOUT)::mat
IF (ALLOCATED(mat%data_c)) DEALLOCATE(mat%data_c)
......
......@@ -4,6 +4,7 @@
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_mkgxyz3
USE m_judft
!.....------------------------------------------------------------------
!c by use of cartesian x,y,z components of charge density gradients,
!c make the quantities
......
......@@ -4,6 +4,7 @@
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_vmt_xc
USE m_judft
!.....------------------------------------------------------------------
! Calculate the GGA xc-potential in the MT-spheres
!.....------------------------------------------------------------------
......
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