scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Again an improvement of the optimization procedure in neb calculations: 

the new version of quick-min estimates the optimal time step on the basis of
the approximate force constant along the displacement vector (this is only done
after some optimization steps).
In all test cases the algorithm is at least two times faster.
Nevertheless the algorithm is not yet optimal.
C.S.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@759 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
sbraccia 2004-04-02 09:45:55 +00:00
parent a3e02b9d06
commit d729a0c79c
5 changed files with 218 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Modules/formats.f90
View File

@ -35,11 +35,11 @@ MODULE formats
scf_fmt_para = "(5X,'cpu = ',I2,'; tcpu = ',F8.2," // &
& "'; self-consistency for image ', I3)", &
run_output = "(/,5X,'iteration:',I4,4X,'E activation ='," // &
& " F6.3,' eV',4X,'error =',F8.4,' eV / bohr'/)", &
& " F6.3,' eV',4X,'error =',F8.4,' eV / A'/)", &
run_output_T_const = "(/,5X,'iteration:',I4,4X,'temperature ='," // &
& " F8.2,' K',4X,'forces =',F8.4,' eV / bohr')", &
& " F8.2,' K',4X,'forces =',F8.4,' eV / A')", &
final_output = "(5X,'image: ',I2,' E tot = ',F16.8," // &
& "' eV error = ',F8.4,' eV / bohr')"
& "' eV error = ',F8.4,' eV / A')"
!
CHARACTER (LEN=*), PARAMETER :: &
stringfmt = "(5X,A,T35,' = ',A)"

151
Modules/minimization_routines.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001-2003 PWSCF group
! Copyright (C) 2003-2004 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
@ -10,10 +10,14 @@
MODULE minimization_routines
!---------------------------------------------------------------------------
!
! ... This module contains all subroutines and functions needed for
! ... the optimization of the reaction path (NEB calculations)
! ... Written by Carlo Sbraccia ( 04-11-2003 )
!
USE kinds, ONLY : DP
USE constants, ONLY : AU, eV_to_kelvin, eps32
USE neb_variables, ONLY : pos, ds, grad, norm_grad
USE basic_algebra_routines, ONLY : norm
USE neb_variables, ONLY : ds, pos, grad, norm_grad
USE basic_algebra_routines
!
IMPLICIT NONE
!
@ -32,10 +36,10 @@ MODULE minimization_routines
!
IF ( norm_grad(index) >= eps32 ) THEN
!
pos(:,index) = pos(:,index) - ds * grad(:,index)
pos(:,index) = pos(:,index) - ds(index) * grad(:,index)
!
END IF
!
!
RETURN
!
END SUBROUTINE steepest_descent
@ -47,15 +51,16 @@ MODULE minimization_routines
SUBROUTINE velocity_Verlet_first_step( index )
!----------------------------------------------------------------------
!
USE neb_variables, ONLY : vel, mass
USE neb_variables, ONLY : vel, mass
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: index
!
!
vel(:,index) = vel(:,index) - ds / 2.D0 * grad(:,index) / mass(:)
pos(:,index) = pos(:,index) + ds * vel(:,index)
!
vel(:,index) = vel(:,index) - &
ds(index) / 2.D0 * grad(:,index) / mass(:)
pos(:,index) = pos(:,index) + ds(index) * vel(:,index)
!
RETURN
!
@ -66,7 +71,7 @@ MODULE minimization_routines
SUBROUTINE velocity_Verlet_second_step( index )
!----------------------------------------------------------------------
!
USE neb_variables, ONLY : vel, mass, damp, ldamped_dyn, lmol_dyn
USE neb_variables, ONLY : vel, mass, damp, ldamped_dyn, lmol_dyn
!
IMPLICIT NONE
!
@ -76,11 +81,12 @@ MODULE minimization_routines
IF ( ldamped_dyn ) THEN
!
vel(:,index) = damp * ( vel(:,index) - &
ds / 2.D0 * grad(:,index) / mass(:) )
ds(index) / 2.D0 * grad(:,index) / mass(:) )
!
ELSE IF ( lmol_dyn ) THEN
!
vel(:,index) = vel(:,index) - ds / 2.D0 * grad(:,index) / mass(:)
vel(:,index) = vel(:,index) - &
ds(index) / 2.D0 * grad(:,index) / mass(:)
!
END IF
!
@ -92,31 +98,61 @@ MODULE minimization_routines
SUBROUTINE quick_min_second_step( index )
!----------------------------------------------------------------------
!
USE neb_variables, ONLY : vel, mass, dim
USE neb_variables, ONLY : pos_old, grad_old, vel, mass, dim
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: index
REAL (KIND=DP), DIMENSION(dim) :: force_versor
REAL (KIND=DP) :: vel_component
INTEGER, INTENT(IN) :: index
REAL (KIND=DP) :: force_versor(dim)
REAL (KIND=DP) :: vel_component
REAL (KIND=DP), ALLOCATABLE :: y(:), s(:)
LOGICAL, SAVE :: new_step = .FALSE.
!
!
ds(index) = optimal_time_step( index )
!
vel(:,index) = vel(:,index) - ds / 2.D0 * grad(:,index) / mass(:)
vel(:,index) = vel(:,index) - &
ds(index) / 2.D0 * grad(:,index) / mass(:)
!
IF ( norm_grad(index) >= eps32 ) THEN
!
force_versor = - grad(:,index) / norm_grad(index)
!
vel_component = DOT_PRODUCT( vel(:,index) , force_versor )
vel_component = ( vel(:,index) .dot. force_versor )
!
IF ( vel_component > 0.D0 ) THEN
IF ( vel_component > 0.D0 .OR. new_step ) THEN
!
vel(:,index) = vel_component * force_versor
!
new_step = .FALSE.
!
ELSE
!
! PRINT '(/5X,"IMAGE = ",I2," resetting velocity"/)', index
!
! ... an approximate newton-raphson step is performed
!
vel(:,index) = 0.D0
!
ALLOCATE( y( dim ), s( dim ) )
!
y = grad(:,index) - grad_old(:,index)
s = pos(:,index) - pos_old(:,index)
!
ds(index) = 1.D0
!
! PRINT '(5X,"projection = ",F7.4)',( s .dot. grad(:,index) ) / &
! ( norm( s ) * norm_grad(index) )
!
grad(:,index) = 2.D0 * s * &
ABS( ( s .dot. grad(:,index) ) / ( y .dot. s ) )
!
! PRINT '(5X,"step length: ",F12.8)', norm( grad(:,index) )
!
new_step = .TRUE.
!
DEALLOCATE( y, s )
!
END IF
!
ELSE
@ -125,6 +161,11 @@ MODULE minimization_routines
!
END IF
!
! ... pos_old and grad_old are updated here
!
pos_old(:,index) = pos(:,index)
grad_old(:,index) = grad(:,index)
!
RETURN
!
END SUBROUTINE quick_min_second_step
@ -135,7 +176,7 @@ MODULE minimization_routines
SUBROUTINE thermalization( N_in , N_fin )
!----------------------------------------------------------------------
!
USE neb_variables, ONLY : vel, mass, temp, temp_req, deg_of_freedom
USE neb_variables, ONLY : vel, mass, temp, temp_req, deg_of_freedom
!
IMPLICIT NONE
!
@ -148,8 +189,8 @@ MODULE minimization_routines
!
DO image = N_in, N_fin
!
local_temp = DOT_PRODUCT( mass(:) * vel(:,image) , &
vel(:,image) ) / REAL( deg_of_freedom )
local_temp = ( ( mass(:) * vel(:,image) ) .dot. vel(:,image) ) / &
REAL( deg_of_freedom )
!
temp = temp + local_temp
!
@ -167,7 +208,7 @@ MODULE minimization_routines
!
DO image = N_in, N_fin
!
temp = temp + DOT_PRODUCT( mass(:) * vel(:,image) , vel(:,image) )
temp = temp + ( ( mass(:) * vel(:,image) ) .dot. vel(:,image) )
!
END DO
!
@ -181,4 +222,68 @@ MODULE minimization_routines
!
END SUBROUTINE thermalization
!
! ... this routine computes the optimal time step using an approximation
! ... of the second derivative along the last displacement vector
!
!------------------------------------------------------------------------
FUNCTION optimal_time_step( index )
!------------------------------------------------------------------------
!
USE constants, ONLY : eps8
USE neb_variables, ONLY : dim, pos_old, grad_old, istep_neb
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: index
REAL(KIND=DP) :: optimal_time_step
REAL(KIND=DP), ALLOCATABLE :: y(:), s(:)
REAL(KIND=DP) :: projection
INTEGER, PARAMETER :: steps_at_fixed_ds = 10
! for the first "steps_at_fixed_ds" steps the update of the time
! step is disabled
REAL(KIND=DP), PARAMETER :: ds_min = 0.4D0, &
ds_max = 4.0D0
! minimum and maximum allowed time-steps
!
!
ALLOCATE( y( dim ), s( dim ) )
!
y = grad(:,index) - grad_old(:,index)
s = pos(:,index) - pos_old(:,index)
!
IF ( istep_neb < steps_at_fixed_ds ) THEN
!
optimal_time_step = ds(index)
!
ELSE
!
IF ( ABS( y .dot. s ) > eps8 ) THEN
!
projection = ( s .dot. grad(:,index) ) / &
( norm( s ) * norm_grad(index) )
!
! PRINT '(5X,"projection = ",F7.4)', projection
!
optimal_time_step = ABS( projection ) * &
SQRT( 2.D0 * ( s .dot. s ) / ABS( y .dot. s ) )
!
ELSE
!
optimal_time_step = ds(index)
!
END IF
!
! PRINT '(5X,"before: ",F12.8)', optimal_time_step
!
optimal_time_step = MAX( ds_min, optimal_time_step )
optimal_time_step = MIN( ds_max, optimal_time_step )
!
! PRINT '(5X,"after: ",F12.8)', optimal_time_step
!
END IF
!
DEALLOCATE( y, s )
!
END FUNCTION optimal_time_step
!
END MODULE minimization_routines

148
Modules/neb_base.f90
View File

@ -34,18 +34,20 @@ MODULE neb_base
!-----------------------------------------------------------------------
!
USE input_parameters, ONLY : pos, nat, restart_mode, calculation, &
minimization_scheme, climbing, nstep
minimization_scheme, climbing, nstep, ds
USE io_files, ONLY : prefix, iunneb, neb_file, &
dat_file, int_file, xyz_file, axsf_file
USE cell_base, ONLY : alat
USE neb_variables, ONLY : pos_ => pos, &
USE neb_variables, ONLY : ds_ => ds, &
pos_ => pos, &
climbing_ => climbing, &
pos_old, grad_old, &
vel, num_of_images, dim, PES, PES_gradient, &
elastic_gradient, tangent, grad, norm_grad, &
error, mass, free_minimization, CI_scheme, &
optimization, k, k_min, k_max, Emax_index, &
VEC_scheme, ds, neb_thr, lquick_min , &
ldamped_dyn, lmol_dyn, nstep_neb, istep_neb
VEC_scheme, neb_thr, lquick_min, lmol_dyn, &
ldamped_dyn, nstep_neb, istep_neb
USE neb_variables, ONLY : neb_dyn_allocation
USE parser, ONLY : int_to_char
USE io_routines, ONLY : read_restart
@ -53,17 +55,17 @@ MODULE neb_base
USE io_global, ONLY : ionode
!
IMPLICIT NONE
CHARACTER(LEN=2) :: prog ! ... specify the calling program
!
! ... input variables
!
CHARACTER(LEN=2) :: prog
! ... specify the calling program
!
! ... local variables
!
INTEGER :: i
REAL (KIND=DP), ALLOCATABLE :: d_R(:)
!
! ... end of local variables
!
!
! ... NEB internal variables are set
!
@ -102,11 +104,13 @@ MODULE neb_base
!
! ... all other arrays are initialized
!
ds_ = ds
PES = 0.D0
PES_gradient = 0.D0
elastic_gradient = 0.D0
tangent = 0.D0
grad = 0.D0
grad_old = 0.D0
norm_grad = 0.D0
error = 0.D0
mass = 1.D0
@ -188,6 +192,10 @@ MODULE neb_base
!
END IF
!
! ... pos_old is initialized
!
pos_old = pos_
!
! ... the actual number of degrees of freedom is computed
!
CALL compute_deg_of_freedom()
@ -253,41 +261,6 @@ MODULE neb_base
!
!
!------------------------------------------------------------------------
SUBROUTINE compute_tangent()
!------------------------------------------------------------------------
!
USE neb_variables, ONLY : num_of_images, tangent
USE basic_algebra_routines, ONLY : norm
!
IMPLICIT NONE
!
! ... local variables
!
INTEGER :: image
!
! ... end of local variables
!
!
tangent = 0
!
DO image = 2, ( num_of_images - 1 )
!
! ... tangent to the path ( normalized )
!
!!! tangent(:,image) = path_tangent( image )
!!! workaround for ifc8 compiler internal error
CALL path_tangent_( image, tangent(:,image) )
!
tangent(:,image) = tangent(:,image) / norm( tangent(:,image) )
!
END DO
!
RETURN
!
END SUBROUTINE compute_tangent
!
!
!------------------------------------------------------------------------
SUBROUTINE elastic_constants()
!------------------------------------------------------------------------
!
@ -302,13 +275,15 @@ MODULE neb_base
!
! ... local variables
!
INTEGER :: i
REAL (KIND=DP) :: F_ortho_max, F_ortho_max_i, &
F_para_max_i, F_para_max
REAL (KIND=DP) :: delta_E
REAL (KIND=DP) :: norm_grad_V, norm_grad_V_min, norm_grad_V_max
!
! ... end of local variables
INTEGER :: i
REAL(KIND=DP) :: F_ortho_max, F_ortho_max_i, &
F_para_max_i, F_para_max, rescale_coeff
REAL(KIND=DP), PARAMETER :: rescale_coeff_min = 0.3D0, &
rescale_coeff_mAX = 3.D0
! minimum allowed rescaling coefficient ( 30% )
! maximum allowed rescaling coefficient ( 300% )
REAL(KIND=DP) :: delta_E
REAL(KIND=DP) :: norm_grad_V, norm_grad_V_min, norm_grad_V_max
!
!
IF ( VEC_scheme == "energy-weighted" ) THEN
@ -376,13 +351,16 @@ MODULE neb_base
!
END DO
!
PRINT '(/5X,"F_ortho_max = ",F10.6)', F_ortho_max
PRINT '(5X,"F_para_max = ",F10.6)', F_para_max
PRINT '(5X,"ALPHA = ",F10.6)', F_ortho_max / F_para_max
rescale_coeff = MAX( ( F_ortho_max / F_para_max ), rescale_coeff_min )
rescale_coeff = MIN( rescale_coeff, rescale_coeff_max )
!
k = k * F_ortho_max / F_para_max
k_max = k_max * F_ortho_max / F_para_max
k_min = k_min * F_ortho_max / F_para_max
PRINT '(/5X,"F_ortho_max = ",F10.6)', F_ortho_max
PRINT '( 5X,"F_para_max = ",F10.6)', F_para_max
PRINT '( 5X,"ALPHA = ",F10.6)', rescale_coeff
!
k = k * rescale_coeff
k_max = k_max * rescale_coeff
k_min = k_min * rescale_coeff
!
RETURN
!
@ -406,8 +384,6 @@ MODULE neb_base
!
INTEGER :: i
!
! ... end of local variables
!
!
CALL elastic_constants()
!
@ -479,8 +455,6 @@ MODULE neb_base
INTEGER :: i
REAL (KIND=DP) :: V_p, V_h, V_n
!
! ... end of local variables
!
!
climbing = .FALSE.
free_minimization = .FALSE.
@ -532,8 +506,6 @@ MODULE neb_base
INTEGER :: N_in, N_fin
INTEGER :: i
!
! ... end of local variables
!
!
err = 0.D0
!
@ -565,6 +537,39 @@ MODULE neb_base
END SUBROUTINE compute_error
!
!
!------------------------------------------------------------------------
SUBROUTINE compute_tangent()
!------------------------------------------------------------------------
!
USE neb_variables, ONLY : num_of_images, tangent
USE basic_algebra_routines, ONLY : norm
!
IMPLICIT NONE
!
! ... local variables
!
INTEGER :: image
!
!
tangent = 0
!
DO image = 2, ( num_of_images - 1 )
!
! ... tangent to the path ( normalized )
!
!!! tangent(:,image) = path_tangent( image )
!!! workaround for ifc8 compiler internal error
CALL path_tangent_( image, tangent(:,image) )
!
tangent(:,image) = tangent(:,image) / norm( tangent(:,image) )
!
END DO
!
RETURN
!
END SUBROUTINE compute_tangent
!
!
!-----------------------------------------------------------------------
!!! FUNCTION path_tangent( index )
!!! workaround for ifc8 compiler internal error
@ -587,8 +592,6 @@ MODULE neb_base
REAL (KIND=DP) :: abs_next, abs_previous
REAL (KIND=DP) :: delta_V_max, delta_V_min
!
! ... end of local variables
!
!
V_previous = PES( index - 1 )
V_actual = PES( index )
@ -637,10 +640,9 @@ MODULE neb_base
!!! workaround for ifc8 compiler internal error
END SUBROUTINE path_tangent_
!
!
!-----------------------------------------------------------------------
!------------------------------------------------------------------------
SUBROUTINE born_oppenheimer_PES( flag, stat )
!-----------------------------------------------------------------------
!------------------------------------------------------------------------
!
USE neb_variables, ONLY : num_of_images, Emax_index, Emin, Emax, &
PES, PES_gradient, suspended_image
@ -657,8 +659,6 @@ MODULE neb_base
INTEGER :: i, image
INTEGER :: N_in, N_fin
!
! ... end of local variables
!
!
IF ( flag ) THEN
!
@ -728,7 +728,7 @@ MODULE neb_base
!
! ... external functions
!
REAL (kind=DP), EXTERNAL :: get_clock
REAL (KIND=DP), EXTERNAL :: get_clock
!
!
conv_neb = .FALSE.
@ -812,6 +812,10 @@ MODULE neb_base
!
END IF
!
! ... istep_neb is updated after a self-consistency step
!
istep_neb = istep_neb + 1
!
IF ( .NOT. stat ) THEN
!
conv_neb = .FALSE.
@ -862,8 +866,6 @@ MODULE neb_base
!
CALL write_dat_files()
!
istep_neb = istep_neb + 1
!
! ... informations are written on the standard output
!
IF ( ionode ) THEN

10
Modules/neb_variables.f90
View File

@ -37,7 +37,6 @@ MODULE neb_variables
! achieved
REAL (KIND=DP) :: &
neb_thr, &! convergence threshold for NEB
ds, &! minimization step
damp, &! damp coefficient
temp, &! actual temperature ( average over images )
temp_req ! required temperature
@ -47,11 +46,14 @@ MODULE neb_variables
ldamped_dyn, &! .TRUE. if minimization_scheme = "damped-dyn"
lmol_dyn ! .TRUE. if minimization_scheme = "mol-dyn"
REAL (KIND=DP), ALLOCATABLE :: &
ds(:), &! minimization step
pos(:,:), &!
pos_old(:,:), &!
vel(:,:), &!
PES(:), &!
PES_gradient(:,:), &!
grad(:,:), &!
grad_old(:,:), &!
norm_grad(:), &!
k(:), &!
elastic_gradient(:), &!
@ -77,7 +79,9 @@ MODULE neb_variables
IMPLICIT NONE
!
!
ALLOCATE( ds( num_of_images ) )
ALLOCATE( pos( dim, num_of_images ) )
ALLOCATE( pos_old( dim, num_of_images ) )
!
IF ( lquick_min .OR. ldamped_dyn .OR. lmol_dyn ) THEN
!
@ -86,6 +90,7 @@ MODULE neb_variables
END IF
!
ALLOCATE( grad( dim, num_of_images ) )
ALLOCATE( grad_old( dim, num_of_images ) )
ALLOCATE( norm_grad( num_of_images ) )
ALLOCATE( PES( num_of_images ) )
ALLOCATE( PES_gradient( dim, num_of_images ) )
@ -107,9 +112,12 @@ MODULE neb_variables
IMPLICIT NONE
!
!
IF ( ALLOCATED( ds ) ) DEALLOCATE( ds )
IF ( ALLOCATED( pos ) ) DEALLOCATE( pos )
IF ( ALLOCATED( pos_old ) ) DEALLOCATE( pos_old )
IF ( ALLOCATED( vel ) ) DEALLOCATE( vel )
IF ( ALLOCATED( grad ) ) DEALLOCATE( grad )
IF ( ALLOCATED( grad_old ) ) DEALLOCATE( grad_old )
IF ( ALLOCATED( norm_grad ) ) DEALLOCATE( norm_grad )
IF ( ALLOCATED( PES ) ) DEALLOCATE( PES )
IF ( ALLOCATED( PES_gradient ) ) DEALLOCATE( PES_gradient )

8
PW/input.f90
View File

@ -28,7 +28,7 @@ SUBROUTINE iosys()
USE bp, ONLY : nppstr_ => nppstr, &
gdir_ => gdir, &
lberry_ => lberry
USE brilz, ONLY : at, alat, omega, &
USE cell_base, ONLY : at, alat, omega, &
celldm_ => celldm, &
ibrav_ => ibrav
USE ions_base, ONLY : ntyp_ => nsp
@ -107,7 +107,6 @@ SUBROUTINE iosys()
optimization_ => optimization, &
damp_ => damp, &
temp_req_ => temp_req, &
ds_ => ds, &
k_max_ => k_max, &
k_min_ => k_min, &
neb_thr_ => neb_thr, &
@ -168,7 +167,7 @@ SUBROUTINE iosys()
tempw, tolp, upscale, potential_extrapolation, &
CI_scheme, VEC_scheme, minimization_scheme, &
num_of_images, optimization, damp, temp_req, &
ds, k_max, k_min, neb_thr, &
k_max, k_min, neb_thr, &
trust_radius_max, trust_radius_min, &
trust_radius_ini, trust_radius_end, &
w_1, w_2, lbfgs_ndim
@ -797,7 +796,6 @@ SUBROUTINE iosys()
optimization_ = optimization
damp_ = damp
temp_req_ = temp_req
ds_ = ds
k_max_ = k_max
k_min_ = k_min
neb_thr_ = neb_thr
@ -1011,7 +1009,7 @@ SUBROUTINE read_cards( psfile, atomic_positions_ )
!-----------------------------------------------------------------------
!
USE wvfct, ONLY : gamma_only
USE brilz, ONLY : at, ibrav, symm_type, celldm
USE cell_base, ONLY : at, ibrav, symm_type, celldm
USE basis, ONLY : nat, ntyp, ityp, tau, atm
USE klist, ONLY : nks
USE ktetra, ONLY : nk1_ => nk1, &