fix cycle dependency for FCP

This commit is contained in:
Satoshi Hagiwara 2020-12-30 14:28:52 +09:00 committed by MinoruOtani
parent f7d2974c33
commit eb945e2b75
2 changed files with 31 additions and 44 deletions

View File

@ -152,7 +152,7 @@ CONTAINS
grad_in, fcell, iforceh, felec, & grad_in, fcell, iforceh, felec, &
energy_thr, grad_thr, cell_thr, fcp_thr, & energy_thr, grad_thr, cell_thr, fcp_thr, &
energy_error, grad_error, cell_error, fcp_error, & energy_error, grad_error, cell_error, fcp_error, &
lmovecell, lfcp, fcp_cap, step_accepted, stop_bfgs, istep ) lmovecell, lfcp, fcp_cap, fcp_hess, step_accepted, stop_bfgs, istep )
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! !
! ... list of input/output arguments : ! ... list of input/output arguments :
@ -193,6 +193,7 @@ CONTAINS
LOGICAL, INTENT(IN) :: lmovecell LOGICAL, INTENT(IN) :: lmovecell
LOGICAL, INTENT(IN) :: lfcp ! include FCP, or not ? LOGICAL, INTENT(IN) :: lfcp ! include FCP, or not ?
REAL(DP), INTENT(IN) :: fcp_cap ! capacitance for FCP REAL(DP), INTENT(IN) :: fcp_cap ! capacitance for FCP
REAL(DP), INTENT(IN) :: fcp_hess ! hessian for FCP
REAL(DP), INTENT(OUT) :: energy_error, grad_error, cell_error, fcp_error REAL(DP), INTENT(OUT) :: energy_error, grad_error, cell_error, fcp_error
LOGICAL, INTENT(OUT) :: step_accepted, stop_bfgs LOGICAL, INTENT(OUT) :: step_accepted, stop_bfgs
INTEGER, INTENT(OUT) :: istep INTEGER, INTENT(OUT) :: istep
@ -283,7 +284,7 @@ CONTAINS
IF ( lmovecell ) fname="enthalpy" IF ( lmovecell ) fname="enthalpy"
IF ( lfcp ) fname = "grand-" // TRIM(fname) IF ( lfcp ) fname = "grand-" // TRIM(fname)
! !
CALL read_bfgs_file( pos, grad, energy, n, lfcp, fcp_cap ) CALL read_bfgs_file( pos, grad, energy, n, lfcp, fcp_hess )
! !
scf_iter = scf_iter + 1 scf_iter = scf_iter + 1
istep = scf_iter istep = scf_iter
@ -411,7 +412,7 @@ CONTAINS
tr_min_hit = 1 tr_min_hit = 1
END IF END IF
! !
CALL reset_bfgs( n, lfcp, fcp_cap ) CALL reset_bfgs( n, lfcp, fcp_hess )
! !
step(:) = - ( inv_hess(:,:) .times. grad(:) ) step(:) = - ( inv_hess(:,:) .times. grad(:) )
if (lmovecell) FORALL( i=1:3, j=1:3) step( n-NADD+j+3*(i-1) ) = step( n-NADD+j+3*(i-1) )*iforceh(i,j) if (lmovecell) FORALL( i=1:3, j=1:3) step( n-NADD+j+3*(i-1) ) = step( n-NADD+j+3*(i-1) )*iforceh(i,j)
@ -454,7 +455,7 @@ CONTAINS
! !
CALL check_wolfe_conditions( lwolfe, energy, grad ) CALL check_wolfe_conditions( lwolfe, energy, grad )
! !
CALL update_inverse_hessian( pos, grad, n, lfcp, fcp_cap ) CALL update_inverse_hessian( pos, grad, n, lfcp, fcp_hess )
! !
END IF END IF
! compute new search direction and store NR step length ! compute new search direction and store NR step length
@ -477,7 +478,7 @@ CONTAINS
WRITE( UNIT = stdout, & WRITE( UNIT = stdout, &
FMT = '(5X,"uphill step: resetting bfgs history",/)' ) FMT = '(5X,"uphill step: resetting bfgs history",/)' )
! !
CALL reset_bfgs( n, lfcp, fcp_cap ) CALL reset_bfgs( n, lfcp, fcp_hess )
step(:) = - ( inv_hess(:,:) .times. grad(:) ) step(:) = - ( inv_hess(:,:) .times. grad(:) )
if (lmovecell) FORALL( i=1:3, j=1:3) step( n-NADD+j+3*(i-1) ) = step( n-NADD+j+3*(i-1) )*iforceh(i,j) if (lmovecell) FORALL( i=1:3, j=1:3) step( n-NADD+j+3*(i-1) ) = step( n-NADD+j+3*(i-1) )*iforceh(i,j)
! !
@ -496,7 +497,7 @@ CONTAINS
! !
ELSE ELSE
! !
CALL compute_trust_radius( lwolfe, energy, grad, n, lfcp, fcp_cap ) CALL compute_trust_radius( lwolfe, energy, grad, n, lfcp, fcp_hess )
! !
END IF END IF
! !
@ -648,32 +649,22 @@ CONTAINS
END SUBROUTINE bfgs END SUBROUTINE bfgs
! !
!------------------------------------------------------------------------ !------------------------------------------------------------------------
SUBROUTINE reset_bfgs( n, lfcp, fcp_cap ) SUBROUTINE reset_bfgs( n, lfcp, fcp_hess )
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! ... inv_hess in re-initialized to the initial guess ! ... inv_hess in re-initialized to the initial guess
! ... defined as the inverse metric ! ... defined as the inverse metric
! !
INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: n
LOGICAL, INTENT(IN) :: lfcp LOGICAL, INTENT(IN) :: lfcp
REAL(DP), INTENT(IN) :: fcp_cap REAL(DP), INTENT(IN) :: fcp_hess
!
REAL(DP) :: helec
! !
inv_hess = inv_metric inv_hess = inv_metric
! !
IF ( lfcp ) THEN IF ( lfcp ) THEN
! !
CALL fcp_hessian(helec) IF ( fcp_hess > eps4 ) THEN
! !
IF ( fcp_cap > eps4 ) THEN inv_hess(n,n) = fcp_hess
!
helec = MIN( fcp_cap, helec )
!
END IF
!
IF ( helec > eps4 ) THEN
!
inv_hess(n,n) = helec
! !
END IF END IF
! !
@ -684,7 +675,7 @@ CONTAINS
END SUBROUTINE reset_bfgs END SUBROUTINE reset_bfgs
! !
!------------------------------------------------------------------------ !------------------------------------------------------------------------
SUBROUTINE read_bfgs_file( pos, grad, energy, n, lfcp, fcp_cap ) SUBROUTINE read_bfgs_file( pos, grad, energy, n, lfcp, fcp_hess )
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! !
IMPLICIT NONE IMPLICIT NONE
@ -693,12 +684,11 @@ CONTAINS
REAL(DP), INTENT(INOUT) :: grad(:) REAL(DP), INTENT(INOUT) :: grad(:)
INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: n
LOGICAL, INTENT(IN) :: lfcp LOGICAL, INTENT(IN) :: lfcp
REAL(DP), INTENT(IN) :: fcp_cap REAL(DP), INTENT(IN) :: fcp_hess
REAL(DP), INTENT(INOUT) :: energy REAL(DP), INTENT(INOUT) :: energy
! !
INTEGER :: iunbfgs INTEGER :: iunbfgs
LOGICAL :: file_exists LOGICAL :: file_exists
REAL(DP) :: helec
! !
! !
INQUIRE( FILE = bfgs_file, EXIST = file_exists ) INQUIRE( FILE = bfgs_file, EXIST = file_exists )
@ -739,18 +729,9 @@ CONTAINS
! !
IF ( lfcp ) THEN IF ( lfcp ) THEN
! !
! replace diagonal element of FCP IF ( fcp_hess > eps4 ) THEN
CALL fcp_hessian(helec)
! !
IF ( fcp_cap > eps4 ) THEN inv_hess(n,n) = fcp_hess
!
helec = MIN( fcp_cap, helec )
!
END IF
!
IF ( helec > eps4 ) THEN
!
inv_hess(n,n) = helec
! !
END IF END IF
! !
@ -806,8 +787,7 @@ CONTAINS
! !
END SUBROUTINE write_bfgs_file END SUBROUTINE write_bfgs_file
! !
!------------------------------------------------------------------------ SUBROUTINE update_inverse_hessian( pos, grad, n, lfcp, fcp_hess )
SUBROUTINE update_inverse_hessian( pos, grad, n, lfcp, fcp_cap )
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! !
IMPLICIT NONE IMPLICIT NONE
@ -816,7 +796,7 @@ CONTAINS
REAL(DP), INTENT(IN) :: grad(:) REAL(DP), INTENT(IN) :: grad(:)
INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: n
LOGICAL, INTENT(IN) :: lfcp LOGICAL, INTENT(IN) :: lfcp
REAL(DP), INTENT(IN) :: fcp_cap REAL(DP), INTENT(IN) :: fcp_hess
INTEGER :: info INTEGER :: info
! !
REAL(DP), ALLOCATABLE :: y(:), s(:) REAL(DP), ALLOCATABLE :: y(:), s(:)
@ -839,7 +819,7 @@ CONTAINS
& "behaviour in update_inverse_hessian")' ) & "behaviour in update_inverse_hessian")' )
WRITE( stdout, '( 5X," resetting bfgs history",/)' ) WRITE( stdout, '( 5X," resetting bfgs history",/)' )
! !
CALL reset_bfgs( n, lfcp, fcp_cap ) CALL reset_bfgs( n, lfcp, fcp_hess )
! !
RETURN RETURN
! !
@ -935,7 +915,7 @@ CONTAINS
END FUNCTION gradient_wolfe_condition END FUNCTION gradient_wolfe_condition
! !
!------------------------------------------------------------------------ !------------------------------------------------------------------------
SUBROUTINE compute_trust_radius( lwolfe, energy, grad, n, lfcp, fcp_cap ) SUBROUTINE compute_trust_radius( lwolfe, energy, grad, n, lfcp, fcp_hess )
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
! !
IMPLICIT NONE IMPLICIT NONE
@ -945,7 +925,7 @@ CONTAINS
REAL(DP), INTENT(IN) :: grad(:) REAL(DP), INTENT(IN) :: grad(:)
INTEGER, INTENT(IN) :: n INTEGER, INTENT(IN) :: n
LOGICAL, INTENT(IN) :: lfcp LOGICAL, INTENT(IN) :: lfcp
REAL(DP), INTENT(IN) :: fcp_cap REAL(DP), INTENT(IN) :: fcp_hess
! !
REAL(DP) :: a REAL(DP) :: a
LOGICAL :: ltest LOGICAL :: ltest
@ -987,7 +967,7 @@ CONTAINS
WRITE( UNIT = stdout, & WRITE( UNIT = stdout, &
FMT = '(5X,"small trust_radius: resetting bfgs history",/)' ) FMT = '(5X,"small trust_radius: resetting bfgs history",/)' )
! !
CALL reset_bfgs( n, lfcp, fcp_cap ) CALL reset_bfgs( n, lfcp, fcp_hess )
step(:) = - ( inv_hess(:,:) .times. grad(:) ) step(:) = - ( inv_hess(:,:) .times. grad(:) )
! !
nr_step_length = scnorm(step) nr_step_length = scnorm(step)

View File

@ -24,7 +24,7 @@ SUBROUTINE move_ions( idone, ions_status )
!! Coefficients for potential and wavefunctions extrapolation are !! Coefficients for potential and wavefunctions extrapolation are
!! no longer computed here but in update_pot. !! no longer computed here but in update_pot.
! !
USE constants, ONLY : e2, eps6, ry_kbar USE constants, ONLY : e2, eps4, eps6, ry_kbar
USE io_global, ONLY : stdout USE io_global, ONLY : stdout
USE io_files, ONLY : tmp_dir, prefix USE io_files, ONLY : tmp_dir, prefix
USE kinds, ONLY : DP USE kinds, ONLY : DP
@ -64,7 +64,7 @@ SUBROUTINE move_ions( idone, ions_status )
LOGICAL :: step_accepted, exst LOGICAL :: step_accepted, exst
REAL(DP), ALLOCATABLE :: pos(:), grad(:) REAL(DP), ALLOCATABLE :: pos(:), grad(:)
REAL(DP) :: h(3,3), fcell(3,3)=0.d0, epsp1 REAL(DP) :: h(3,3), fcell(3,3)=0.d0, epsp1
REAL(DP) :: relec, felec, capacitance, tot_charge_ REAL(DP) :: relec, felec, helec, capacitance, tot_charge_
LOGICAL :: conv_ions LOGICAL :: conv_ions
CHARACTER(LEN=320) :: filebfgs CHARACTER(LEN=320) :: filebfgs
! !
@ -114,6 +114,13 @@ SUBROUTINE move_ions( idone, ions_status )
felec = (ef - fcp_mu) felec = (ef - fcp_mu)
CALL fcp_capacitance( capacitance, -1.0_DP ) CALL fcp_capacitance( capacitance, -1.0_DP )
tot_charge_ = tot_charge tot_charge_ = tot_charge
!
! Make hessian for FCP
CALL fcp_hessian( helec )
IF ( capacitance > eps4 ) THEN
helec = MIN( capacitance, helec )
END IF
!
END IF END IF
! !
IF ( ANY( if_pos(:,:) == 1 ) .OR. lmovecell .OR. lfcp ) THEN IF ( ANY( if_pos(:,:) == 1 ) .OR. lmovecell .OR. lfcp ) THEN
@ -121,7 +128,7 @@ SUBROUTINE move_ions( idone, ions_status )
filebfgs = TRIM(tmp_dir) // TRIM(prefix) // '.bfgs' filebfgs = TRIM(tmp_dir) // TRIM(prefix) // '.bfgs'
CALL bfgs( filebfgs, pos, h, relec, etot, grad, fcell, iforceh, felec, epse, & CALL bfgs( filebfgs, pos, h, relec, etot, grad, fcell, iforceh, felec, epse, &
epsf, epsp1, fcp_eps, energy_error, gradient_error, cell_error, fcp_error, & epsf, epsp1, fcp_eps, energy_error, gradient_error, cell_error, fcp_error, &
lmovecell, lfcp, capacitance, step_accepted, conv_ions, istep ) lmovecell, lfcp, capacitance, helec, step_accepted, conv_ions, istep )
! !
ELSE ELSE
! !