Fixed some bugs in NEB. Added a new field in the neb restart file.

To use old restart files the following gawk-script can be used:

BEGIN{ level = 0 }
{
if ( $1 == "Image:" ) {
   print ;
   if ( $2 > level ) {
      level = $2 ;
      getline ;
      printf "%2s,  F \n", $1 ;
      }
   }
else { print }
}

C.S.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1088 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
sbraccia 2004-07-20 15:25:45 +00:00
parent 3baf4ce3bc
commit 122c392223
5 changed files with 120 additions and 74 deletions

View File

@ -28,7 +28,8 @@ MODULE io_routines
PES_gradient, suspended_image, Emax, &
Emin, Emax_index, istep_neb, nstep_neb, &
lquick_min , ldamped_dyn, lmol_dyn, &
reset_vel, k_min, k_max, pos_old, grad_old
k_min, k_max, pos_old, grad_old, &
reset_vel, vel_zeroed
USE io_global, ONLY : ionode, ionode_id
USE mp, ONLY : mp_bcast
!
@ -36,7 +37,7 @@ MODULE io_routines
!
! ... local variables
!
INTEGER :: i, j, ia
INTEGER :: i, j, ia, ierr
CHARACTER (LEN=256) :: input_line
LOGICAL, EXTERNAL :: matches
!
@ -92,7 +93,7 @@ MODULE io_routines
END IF
!
READ( UNIT = iunrestart, FMT = * )
READ( UNIT = iunrestart, FMT = * ) frozen(1)
READ( UNIT = iunrestart, FMT = * ) frozen(1), vel_zeroed(1)
READ( UNIT = iunrestart, FMT = * ) PES(1)
!
ia = 0
@ -124,7 +125,7 @@ MODULE io_routines
DO i = 2, num_of_images
!
READ( UNIT = iunrestart, FMT = * )
READ( UNIT = iunrestart, FMT = * ) frozen(i)
READ( UNIT = iunrestart, FMT = * ) frozen(i), vel_zeroed(i)
READ( UNIT = iunrestart, FMT = * ) PES(i)
!
DO j = 1, dim, 3
@ -154,11 +155,11 @@ MODULE io_routines
!
END DO
!
IF ( .NOT. reset_vel .AND. &
READ( UNIT = iunrestart, FMT = '(256A)', IOSTAT = ierr ) input_line
!
IF ( .NOT. reset_vel .AND. ( ierr == 0 ) .AND. &
( lquick_min .OR. ldamped_dyn .OR. lmol_dyn ) ) THEN
!
READ( UNIT = iunrestart, FMT = '(256A)' ) input_line
!
IF ( matches( "QUICK-MIN FIELDS", input_line ) ) THEN
!
! ... optional fields
@ -250,7 +251,7 @@ MODULE io_routines
PES_gradient, dim, suspended_image, &
lquick_min , ldamped_dyn, lmol_dyn, &
istep_neb, nstep_neb, k_min, k_max, &
pos_old, grad_old
pos_old, grad_old, vel_zeroed
USE formats, ONLY : energy, restart_first, restart_others, &
quick_min
USE io_global, ONLY : ionode
@ -291,7 +292,8 @@ MODULE io_routines
DO i = 1, num_of_images
!
WRITE( UNIT = iunrestart, FMT = '("Image: ",I4)' ) i
WRITE( UNIT = iunrestart, FMT = '(L1)' ) frozen(i)
WRITE( UNIT = iunrestart, FMT = '(2(L1,X))' ) frozen(i), &
vel_zeroed(i)
WRITE( UNIT = iunrestart, FMT = energy ) PES(i)
!
ia = 0
@ -387,7 +389,8 @@ MODULE io_routines
DO i = 1, num_of_images
!
WRITE( UNIT = iunrestart, FMT = '("Image: ",I4)' ) i
WRITE( UNIT = iunrestart, FMT = '(L1)' ) frozen(i)
WRITE( UNIT = iunrestart, FMT = '(2(L1,X))' ) frozen(i), &
vel_zeroed(i)
WRITE( UNIT = iunrestart, FMT = energy ) PES(i)
!
ia = 0

View File

@ -5,6 +5,8 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
#define USE_SMART_STEP
!#define DEBUG_SMART_STEP
!
!--------------------------------------------------------------------------
MODULE minimization_routines
@ -15,7 +17,7 @@ MODULE minimization_routines
! ... Written by Carlo Sbraccia ( 04-11-2003 )
!
USE kinds, ONLY : DP
USE constants, ONLY : AU, eV_to_kelvin, eps32
USE constants, ONLY : au, eV_to_kelvin, eps32
USE neb_variables, ONLY : ds, pos, grad, norm_grad
USE basic_algebra_routines
!
@ -34,7 +36,7 @@ MODULE minimization_routines
INTEGER, INTENT(IN) :: index
!
!
IF ( norm_grad(index) >= eps32 ) THEN
IF ( norm_grad(index) > eps32 ) THEN
!
pos(:,index) = pos(:,index) - ds(index) * grad(:,index)
!
@ -51,16 +53,29 @@ MODULE minimization_routines
SUBROUTINE velocity_Verlet_first_step( index )
!----------------------------------------------------------------------
!
USE neb_variables, ONLY : vel, mass
USE neb_variables, ONLY : vel, mass, vel_zeroed
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: index
!
!
vel(:,index) = vel(:,index) - &
ds(index) / 2.D0 * grad(:,index) / mass(:)
pos(:,index) = pos(:,index) + ds(index) * vel(:,index)
IF ( vel_zeroed(index) ) THEN
!
vel(:,index) = vel(:,index) - 0.5D0 * grad(:,index) / mass(:)
!
pos(:,index) = pos(:,index) + vel(:,index)
!
vel_zeroed(index) = .FALSE.
!
ELSE
!
vel(:,index) = vel(:,index) - &
0.5D0 * ds(index) * grad(:,index) / mass(:)
!
pos(:,index) = pos(:,index) + ds(index) * vel(:,index)
!
END IF
!
RETURN
!
@ -81,12 +96,12 @@ MODULE minimization_routines
IF ( ldamped_dyn ) THEN
!
vel(:,index) = damp * ( vel(:,index) - &
ds(index) / 2.D0 * grad(:,index) / mass(:) )
0.5D0 * ds(index) * grad(:,index) / mass(:) )
!
ELSE IF ( lmol_dyn ) THEN
!
vel(:,index) = vel(:,index) - &
ds(index) / 2.D0 * grad(:,index) / mass(:)
0.5D0 * ds(index) * grad(:,index) / mass(:)
!
END IF
!
@ -98,8 +113,9 @@ MODULE minimization_routines
SUBROUTINE quick_min_second_step( index )
!----------------------------------------------------------------------
!
USE constants, ONLY : eps8
USE neb_variables, ONLY : pos_old, grad_old, vel, mass, dim
USE constants, ONLY : eps8
USE neb_variables, ONLY : pos_old, grad_old, vel, &
mass, dim, vel_zeroed
!
IMPLICIT NONE
!
@ -108,13 +124,11 @@ MODULE minimization_routines
REAL (KIND=DP) :: vel_component
REAL (KIND=DP), ALLOCATABLE :: y(:), s(:)
!
!
! ds(index) = optimal_time_step( index )
!
vel(:,index) = vel(:,index) - &
ds(index) / 2.D0 * grad(:,index) / mass(:)
0.5D0 * ds(index) * grad(:,index) / mass(:)
!
IF ( norm_grad(index) >= eps32 ) THEN
IF ( norm_grad(index) > eps32 ) THEN
!
force_versor = - grad(:,index) / norm_grad(index)
!
@ -126,10 +140,16 @@ MODULE minimization_routines
!
ELSE
!
! PRINT '(/5X,"IMAGE = ",I2," resetting velocity"/)', index
#if defined (DEBUG_SMART_STEP) && defined (USE_SMART_STEP)
PRINT '(/5X,"IMAGE = ",I2," resetting velocity"/)', index
#endif
!
vel(:,index) = 0.D0
!
!
vel_zeroed(index) = .TRUE.
!
#if defined (USE_SMART_STEP)
!
! ... an approximate newton-raphson step is performed
!
IF ( norm( pos_old(:,index) ) > 0.D0 .AND. &
@ -140,22 +160,26 @@ MODULE minimization_routines
y = grad(:,index) - grad_old(:,index)
s = pos(:,index) - pos_old(:,index)
!
! PRINT '(5X,"projection = ",F7.4)',( s .dot. grad(:,index) ) / &
! ( norm( s ) * norm_grad(index) )
# if defined (DEBUG_SMART_STEP)
PRINT '(5X,"projection = ",F7.4)', &
( s .dot. grad(:,index) ) / ( norm( s ) * norm_grad(index) )
# endif
!
IF ( ABS( y .dot. s ) > eps8 ) THEN
!
ds(index) = 1.D0
!
grad(:,index) = 2.D0 * s * &
ABS( ( s .dot. grad(:,index) ) / ( y .dot. s ) )
!
END IF
!
! PRINT '(5X,"step length: ",F12.8)', norm( grad(:,index) )
!
# if defined (DEBUG_SMART_STEP)
PRINT '(5X,"step length: ",F12.8)', norm( grad(:,index) )
# endif
!
DEALLOCATE( y, s )
!
#endif
!
END IF
!
END IF
@ -207,7 +231,7 @@ MODULE minimization_routines
!
#if defined ( _DEBUG ) || ( _DEBUG_THERMALIZATION )
!
PRINT *, "temperature before thermalization", temp * AU * eV_to_kelvin
PRINT *, "temperature before thermalization", temp * au * eV_to_kelvin
!
temp = 0.D0
!
@ -219,7 +243,7 @@ MODULE minimization_routines
!
temp = temp / REAL( ( N_fin - N_in ) * deg_of_freedom )
!
PRINT *, "temperature after thermalization", temp * AU * eV_to_kelvin
PRINT *, "temperature after thermalization", temp * au * eV_to_kelvin
!
#endif
!

View File

@ -5,6 +5,7 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
#define USE_ELASTIC_CONSTANTS_RESCALING
!#define DEBUG_ELASTIC_CONSTANTS
!
!-----------------------------------------------------------------------
@ -102,7 +103,7 @@ MODULE neb_base
optimization, k, k_min, k_max, Emax_index, &
VEC_scheme, neb_thr, lquick_min, lmol_dyn, &
ldamped_dyn, nstep_neb, istep_neb, &
suspended_image
suspended_image, vel_zeroed
USE neb_variables, ONLY : neb_dyn_allocation
USE parser, ONLY : int_to_char
USE io_routines, ONLY : read_restart
@ -185,6 +186,8 @@ MODULE neb_base
!
vel = 0.D0
!
vel_zeroed = .FALSE.
!
END IF
!
! ... initial path is read ( restart_mode == "restart" )
@ -322,8 +325,6 @@ MODULE neb_base
"restart_mode", TRIM( restart_mode )
WRITE( UNIT = iunneb, FMT = stringfmt ) &
"CI_scheme", TRIM( CI_scheme )
WRITE( UNIT = iunneb, FMT = stringfmt ) &
"VEC_scheme", TRIM( VEC_scheme )
WRITE( UNIT = iunneb, FMT = stringfmt ) &
"minimization_scheme", TRIM( minimization_scheme )
WRITE( UNIT = iunneb, FMT = stringfmt ) &
@ -408,7 +409,7 @@ MODULE neb_base
USE constants, ONLY : pi, eps32
USE neb_variables, ONLY : pos, num_of_images, Emax, Emin, &
k_max, k_min, k, PES, PES_gradient, &
VEC_scheme, elastic_gradient, tangent
elastic_gradient, tangent
USE supercell, ONLY : pbc
USE basic_algebra_routines
!
@ -419,7 +420,8 @@ MODULE neb_base
INTEGER :: i
REAL(KIND=DP) :: F_ortho_max, F_ortho_max_i, &
F_para_max_i, F_para_max, rescale_coeff
REAL(KIND=DP) :: delta_E
REAL(KIND=DP) :: delta_E, delta_norm_grad
REAL(KIND=DP) :: k_sum, k_diff
REAL(KIND=DP) :: norm_grad_V, norm_grad_V_min, norm_grad_V_max
!
! ... local parameters
@ -433,54 +435,58 @@ MODULE neb_base
!
!
rescale_coeff = k_max / k_min
!
k_min = MAX( k_min, k_minimal )
k_max = MAX( k_max, k_minimal * rescale_coeff )
!
IF ( VEC_scheme == "energy-weighted" ) THEN
!
delta_E = Emax - Emin
!
IF ( delta_E <= eps32 ) THEN
!
k = k_min
!
RETURN
!
END IF
delta_E = Emax - Emin
!
k_sum = k_max + k_min
k_diff = k_max - k_min
!
k(:) = k_min
!
IF ( delta_E > eps32 ) THEN
!
DO i = 1, num_of_images
!
k(i) = 0.5D0 * ( ( k_max + k_min ) - ( k_max - k_min ) * &
COS( pi * ( PES(i) - Emin ) / delta_E ) )
k(i) = 0.5D0 * ( k_sum - k_diff * &
COS( pi * ( PES(i) - Emin ) / delta_E ) )
!
END DO
!
ELSE
END IF
!
norm_grad_V_min = + 1.0D32
norm_grad_V_max = - 1.0D32
!
DO i = 1, num_of_images
!
norm_grad_V = norm( PES_gradient(:,i) )
!
norm_grad_V_min = + 1.0D32
norm_grad_V_max = - 1.0D32
IF ( norm_grad_V < norm_grad_V_min ) norm_grad_V_min = norm_grad_V
IF ( norm_grad_V > norm_grad_V_max ) norm_grad_V_max = norm_grad_V
!
END DO
!
delta_norm_grad = norm_grad_V_max - norm_grad_V_min
!
IF ( delta_norm_grad > eps32 ) THEN
!
DO i = 1, num_of_images
!
norm_grad_V = norm( PES_gradient(:,i) )
!
IF ( norm_grad_V < norm_grad_V_min ) norm_grad_V_min = norm_grad_V
IF ( norm_grad_V > norm_grad_V_max ) norm_grad_V_max = norm_grad_V
!
END DO
!
DO i = 1, num_of_images
!
norm_grad_V = norm( PES_gradient(:,i) )
!
k(i) = 0.5D0 * ( ( k_max + k_min ) - ( k_max - k_min ) * &
COS( pi * ( norm_grad_V - norm_grad_V_min ) / &
( norm_grad_V_max - norm_grad_V_min ) ) )
k(i) = k(i) + 0.5D0 * ( k_sum - k_diff * &
COS( pi * ( norm_grad_V - norm_grad_V_min ) / &
delta_norm_grad ) )
!
END DO
END DO
!
END IF
!
k(:) = 0.5D0 * k(:)
!
F_ortho_max = 0.D0
F_para_max = 0.D0
!
@ -503,10 +509,14 @@ MODULE neb_base
rescale_coeff = MAX( ( F_ortho_max / F_para_max ), rescale_coeff_min )
rescale_coeff = MIN( rescale_coeff, rescale_coeff_max )
!
#if defined (USE_ELASTIC_CONSTANTS_RESCALING)
!
k = k * rescale_coeff
k_max = k_max * rescale_coeff
k_min = k_min * rescale_coeff
!
#endif
!
#if defined (DEBUG_ELASTIC_CONSTANTS)
!
PRINT '(/5X,"F_ortho_max = ",F10.6 )', F_ortho_max
@ -1067,10 +1077,15 @@ MODULE neb_base
END IF
!
! ... the programs checks if the maximum number of iterations has
! ... been reached or if the user has required a soft exit
! ... been reached
!
IF ( istep_neb >= nstep_neb ) THEN
!
IF ( ionode ) &
WRITE( UNIT = iunneb, &
FMT = '(/,5X,"NEB: reached the maximum number of ", &
& "steps")' )
!
suspended_image = 0
!
EXIT minimization

View File

@ -27,7 +27,9 @@ MODULE neb_variables
climbing(:), &! .TRUE. if the image is required to climb
free_minimization(:), &! .TRUE. if the image is required to be
! optimized (no springs, no projections)
frozen(:) ! .TRUE. if the image is not optimized
frozen(:), &! .TRUE. if the image is not optimized
vel_zeroed(:) ! .TRUE. if the velocity of this image has
! been reset
CHARACTER (LEN=20) :: &
CI_scheme, &! Climbing Image scheme
VEC_scheme ! Variable Elastic Constant scheme
@ -71,8 +73,8 @@ MODULE neb_variables
Emax, &!
Emin !
!
INTEGER :: istep_neb
INTEGER :: nstep_neb
INTEGER :: istep_neb
INTEGER :: nstep_neb
!
CONTAINS
!
@ -90,7 +92,8 @@ MODULE neb_variables
!
IF ( lquick_min .OR. ldamped_dyn .OR. lmol_dyn ) THEN
!
ALLOCATE( vel( dim, num_of_images ) )
ALLOCATE( vel( dim, num_of_images ) )
ALLOCATE( vel_zeroed( num_of_images ) )
!
END IF
!
@ -137,6 +140,7 @@ MODULE neb_variables
IF ( ALLOCATED( climbing ) ) DEALLOCATE( climbing )
IF ( ALLOCATED( free_minimization ) ) DEALLOCATE( free_minimization )
IF ( ALLOCATED( frozen ) ) DEALLOCATE( frozen )
IF ( ALLOCATED( vel_zeroed ) ) DEALLOCATE( vel_zeroed )
!
END SUBROUTINE neb_deallocation
!

View File

@ -311,7 +311,7 @@ PROGRAM images_interpolator
DO i = 1, new_num_of_images
!
WRITE( UNIT = iunrestart, FMT = '("Image: ",I4)' ) i
WRITE( UNIT = iunrestart, FMT = '("F")' )
WRITE( UNIT = iunrestart, FMT = '("F, F")' )
WRITE( UNIT = iunrestart, FMT = energy ) new_V(i)
!
ia = 0