mirror of https://gitlab.com/QEF/q-e.git
old path_....f90 routines has been moved to path_...._pre.f90. Only
used by CP at the moment git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@7033 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
parent
49a80d0d7f
commit
443d256c7c
File diff suppressed because it is too large
Load Diff
|
@ -0,0 +1,52 @@
|
|||
!
|
||||
! Copyright (C) 2003-2005 Quantum ESPRESSO 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,
|
||||
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||
!
|
||||
!----------------------------------------------------------------------------
|
||||
MODULE path_formats_pre
|
||||
!----------------------------------------------------------------------------
|
||||
!
|
||||
! ... this module contains the I/O formats used by all "path"-routines
|
||||
!
|
||||
! ... Written by Carlo Sbraccia ( 2003-2005 )
|
||||
!
|
||||
CHARACTER (LEN=*), PARAMETER :: &
|
||||
lattice_vectors = "(3(2X,F14.10),/,3(2X,F14.10),/,3(2X,F14.10))"
|
||||
!
|
||||
CHARACTER (LEN=*), PARAMETER :: &
|
||||
restart_first = "(3(2X,F18.12),3(2X,F18.12),3(2X,I1))", &
|
||||
restart_others = "(3(2X,F18.12),3(2X,F18.12))"
|
||||
!
|
||||
CHARACTER (LEN=*), PARAMETER :: &
|
||||
quick_min = "(9(2X,F18.12))", &
|
||||
energy = "(2X,F18.10)"
|
||||
!
|
||||
CHARACTER (LEN=*), PARAMETER :: &
|
||||
dat_fmt = "(3(2X,F16.10))", &
|
||||
int_fmt = "(2(2X,F16.10))", &
|
||||
xyz_fmt = "(A2,3(2X,F14.10))", &
|
||||
axsf_fmt = "(A2,6(2X,F14.10))"
|
||||
!
|
||||
CHARACTER (LEN=*), PARAMETER :: &
|
||||
scf_iter_fmt = "(/,5X,30('-'),(1X,'iteration ',I3,1X),30('-'),/)", &
|
||||
scf_fmt = "(5X,'tcpu = ',F8.1," // &
|
||||
& "' self-consistency for image ', I3)", &
|
||||
scf_fmt_para = "(5X,'cpu = ',I2,' tcpu = ',F8.1," // &
|
||||
& "' self-consistency for image ', I3)"
|
||||
!
|
||||
CHARACTER (LEN=*), PARAMETER :: &
|
||||
run_info = "(5X,'image',8X,'energy (eV)',8X,'error (eV/A)',8X,'frozen',/)"
|
||||
!
|
||||
CHARACTER (LEN=*), PARAMETER :: &
|
||||
run_output = "(5X,I5,4X,F15.7,10X,F10.6,12X,L1)"
|
||||
!
|
||||
CHARACTER (LEN=*), PARAMETER :: &
|
||||
summary_fmt = "(5X,A,T35,' = ',1X,A)"
|
||||
!
|
||||
CHARACTER (LEN=*), PARAMETER :: &
|
||||
final_fmt = "(/,5X,75('-'),/)"
|
||||
!
|
||||
END MODULE path_formats_pre
|
File diff suppressed because it is too large
Load Diff
|
@ -0,0 +1,499 @@
|
|||
!
|
||||
! Copyright (C) 2003-2006 Quantum ESPRESSO 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,
|
||||
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||
!
|
||||
!--------------------------------------------------------------------------
|
||||
MODULE path_opt_routines_pre
|
||||
!---------------------------------------------------------------------------
|
||||
!
|
||||
! ... This module contains all subroutines and functions needed for
|
||||
! ... the optimisation of the reaction path (NEB and SMD calculations)
|
||||
!
|
||||
! ... Written by Carlo Sbraccia ( 2003-2006 )
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE constants, ONLY : eps8, eps16
|
||||
USE path_variables_pre, ONLY : ds, pos, grad
|
||||
USE io_global, ONLY : meta_ionode, meta_ionode_id
|
||||
USE mp, ONLY : mp_bcast
|
||||
!
|
||||
USE basic_algebra_routines
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
PRIVATE
|
||||
!
|
||||
PUBLIC :: langevin, steepest_descent, quick_min, broyden, broyden2
|
||||
!
|
||||
CONTAINS
|
||||
!
|
||||
!----------------------------------------------------------------------
|
||||
SUBROUTINE langevin( idx )
|
||||
!----------------------------------------------------------------------
|
||||
!
|
||||
USE path_variables_pre, ONLY : lang
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
INTEGER, INTENT(IN) :: idx
|
||||
!
|
||||
IF ( meta_ionode ) THEN
|
||||
!
|
||||
pos(:,idx) = pos(:,idx) - ds*grad(:,idx) + lang(:,idx)
|
||||
!
|
||||
END IF
|
||||
!
|
||||
CALL mp_bcast( pos, meta_ionode_id )
|
||||
!
|
||||
RETURN
|
||||
!
|
||||
END SUBROUTINE langevin
|
||||
!
|
||||
!----------------------------------------------------------------------
|
||||
SUBROUTINE steepest_descent( idx )
|
||||
!----------------------------------------------------------------------
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
INTEGER, INTENT(IN) :: idx
|
||||
!
|
||||
IF ( meta_ionode ) THEN
|
||||
!
|
||||
pos(:,idx) = pos(:,idx) - ds*ds*grad(:,idx)
|
||||
!
|
||||
END IF
|
||||
!
|
||||
CALL mp_bcast( pos, meta_ionode_id )
|
||||
!
|
||||
RETURN
|
||||
!
|
||||
END SUBROUTINE steepest_descent
|
||||
!
|
||||
!----------------------------------------------------------------------
|
||||
SUBROUTINE quick_min( idx, istep )
|
||||
!----------------------------------------------------------------------
|
||||
!
|
||||
! ... projected Verlet algorithm
|
||||
!
|
||||
USE path_variables_pre, ONLY : dim1, posold
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
INTEGER, INTENT(IN) :: idx, istep
|
||||
!
|
||||
REAL(DP), ALLOCATABLE :: vel(:), force_versor(:), step(:)
|
||||
REAL(DP) :: projection, norm_grad, norm_vel, norm_step
|
||||
!
|
||||
REAL(DP), PARAMETER :: max_step = 0.6_DP ! in bohr
|
||||
!
|
||||
!
|
||||
IF ( meta_ionode ) THEN
|
||||
!
|
||||
ALLOCATE( vel( dim1 ), force_versor( dim1 ), step( dim1 ) )
|
||||
!
|
||||
vel(:) = pos(:,idx) - posold(:,idx)
|
||||
!
|
||||
norm_grad = norm( grad(:,idx) )
|
||||
!
|
||||
norm_vel = norm( vel(:) )
|
||||
!
|
||||
IF ( norm_grad > eps16 .AND. norm_vel > eps16 ) THEN
|
||||
!
|
||||
force_versor(:) = - grad(:,idx) / norm_grad
|
||||
!
|
||||
projection = ( vel(:) .dot. force_versor(:) )
|
||||
!
|
||||
vel(:) = MAX( 0.0_DP, projection ) * force_versor(:)
|
||||
!
|
||||
ELSE
|
||||
!
|
||||
vel(:) = 0.0_DP
|
||||
!
|
||||
END IF
|
||||
!
|
||||
posold(:,idx) = pos(:,idx)
|
||||
!
|
||||
step(:) = vel(:) - ds*ds*grad(:,idx)
|
||||
!
|
||||
norm_step = norm( step(:) )
|
||||
!
|
||||
step(:) = step(:) / norm_step
|
||||
!
|
||||
pos(:,idx) = pos(:,idx) + step(:) * MIN( norm_step, max_step )
|
||||
!
|
||||
DEALLOCATE( vel, force_versor, step )
|
||||
!
|
||||
END IF
|
||||
!
|
||||
CALL mp_bcast( pos, meta_ionode_id )
|
||||
CALL mp_bcast( posold, meta_ionode_id )
|
||||
!
|
||||
RETURN
|
||||
!
|
||||
END SUBROUTINE quick_min
|
||||
!
|
||||
! ... Broyden (rank one) optimisation
|
||||
!
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE broyden()
|
||||
!-----------------------------------------------------------------------
|
||||
!
|
||||
USE control_flags, ONLY : lsmd
|
||||
USE io_files, ONLY : broy_file, iunbroy, iunpath
|
||||
USE path_variables_pre, ONLY : dim1, frozen, tangent, nim => num_of_images
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
REAL(DP), ALLOCATABLE :: t(:), g(:), s(:,:)
|
||||
INTEGER :: k, i, j, I_in, I_fin
|
||||
REAL(DP) :: s_norm, coeff, norm_g
|
||||
REAL(DP) :: J0
|
||||
LOGICAL :: exists
|
||||
!
|
||||
REAL(DP), PARAMETER :: step_max = 0.6_DP
|
||||
INTEGER, PARAMETER :: broyden_ndim = 5
|
||||
!
|
||||
!
|
||||
! ... starting guess for the inverse Jacobian
|
||||
!
|
||||
J0 = ds*ds
|
||||
!
|
||||
ALLOCATE( g( dim1*nim ) )
|
||||
ALLOCATE( s( dim1*nim, broyden_ndim ) )
|
||||
ALLOCATE( t( dim1*nim ) )
|
||||
!
|
||||
g(:) = 0.0_DP
|
||||
t(:) = 0.0_DP
|
||||
!
|
||||
DO i = 1, nim
|
||||
!
|
||||
I_in = ( i - 1 )*dim1 + 1
|
||||
I_fin = i * dim1
|
||||
!
|
||||
IF ( frozen(i) ) CYCLE
|
||||
!
|
||||
IF ( lsmd ) t(I_in:I_fin) = tangent(:,i)
|
||||
!
|
||||
g(I_in:I_fin) = grad(:,i)
|
||||
!
|
||||
END DO
|
||||
!
|
||||
norm_g = MAXVAL( ABS( g ) )
|
||||
!
|
||||
IF ( norm_g == 0.0_DP ) RETURN
|
||||
!
|
||||
IF ( meta_ionode ) THEN
|
||||
!
|
||||
! ... open the file containing the broyden's history
|
||||
!
|
||||
INQUIRE( FILE = broy_file, EXIST = exists )
|
||||
!
|
||||
IF ( exists ) THEN
|
||||
!
|
||||
OPEN( UNIT = iunbroy, FILE = broy_file, STATUS = "OLD" )
|
||||
!
|
||||
READ( UNIT = iunbroy , FMT = * ) i
|
||||
!
|
||||
! ... if the number of images is changed the broyden history is
|
||||
! ... reset and the algorithm starts from scratch
|
||||
!
|
||||
exists = ( i == nim )
|
||||
!
|
||||
END IF
|
||||
!
|
||||
IF ( exists ) THEN
|
||||
!
|
||||
READ( UNIT = iunbroy , FMT = * ) k
|
||||
READ( UNIT = iunbroy , FMT = * ) s(:,:)
|
||||
!
|
||||
k = k + 1
|
||||
!
|
||||
ELSE
|
||||
!
|
||||
s(:,:) = 0.0_DP
|
||||
!
|
||||
k = 1
|
||||
!
|
||||
END IF
|
||||
!
|
||||
CLOSE( UNIT = iunbroy )
|
||||
!
|
||||
! ... Broyden's update
|
||||
!
|
||||
IF ( k > broyden_ndim ) THEN
|
||||
!
|
||||
! ... the Broyden's subspace is swapped and the projection of
|
||||
! ... s along the current tangent is removed (this last thing
|
||||
! ... in the smd case only, otherwise t = 0.0_DP)
|
||||
!
|
||||
k = broyden_ndim
|
||||
!
|
||||
DO j = 1, k - 1
|
||||
!
|
||||
s(:,j) = s(:,j+1) - t(:) * ( s(:,j+1) .dot. t(:) )
|
||||
!
|
||||
END DO
|
||||
!
|
||||
END IF
|
||||
!
|
||||
s(:,k) = - J0 * g(:)
|
||||
!
|
||||
IF ( k > 1 ) THEN
|
||||
!
|
||||
DO j = 1, k - 2
|
||||
!
|
||||
s(:,k) = s(:,k) + ( s(:,j) .dot. s(:,k) ) / &
|
||||
( s(:,j) .dot. s(:,j) ) * s(:,j+1)
|
||||
!
|
||||
END DO
|
||||
!
|
||||
coeff = ( s(:,k-1) .dot. ( s(:,k-1) - s(:,k) ) )
|
||||
!
|
||||
IF ( coeff > eps8 ) THEN
|
||||
!
|
||||
s(:,k) = ( s(:,k-1) .dot. s(:,k-1) ) / coeff * s(:,k)
|
||||
!
|
||||
END IF
|
||||
!
|
||||
END IF
|
||||
!
|
||||
IF ( ( s(:,k) .dot. g(:) ) > 0.0_DP ) THEN
|
||||
!
|
||||
! ... uphill step : history reset
|
||||
!
|
||||
WRITE( UNIT = iunpath, &
|
||||
FMT = '(/,5X,"broyden uphill step : history is reset",/)' )
|
||||
!
|
||||
k = 1
|
||||
!
|
||||
s(:,:) = 0.0_DP
|
||||
s(:,k) = - J0 * g(:)
|
||||
!
|
||||
END IF
|
||||
!
|
||||
s_norm = norm( s(:,k) )
|
||||
!
|
||||
s(:,k) = s(:,k) / s_norm * MIN( s_norm, step_max )
|
||||
!
|
||||
! ... save the file containing the history
|
||||
!
|
||||
OPEN( UNIT = iunbroy, FILE = broy_file )
|
||||
!
|
||||
WRITE( UNIT = iunbroy, FMT = * ) nim
|
||||
WRITE( UNIT = iunbroy, FMT = * ) k
|
||||
WRITE( UNIT = iunbroy, FMT = * ) s
|
||||
!
|
||||
CLOSE( UNIT = iunbroy )
|
||||
!
|
||||
! ... broyden's step
|
||||
!
|
||||
pos(:,1:nim) = pos(:,1:nim) + RESHAPE( s(:,k), (/ dim1, nim /) )
|
||||
!
|
||||
END IF
|
||||
!
|
||||
CALL mp_bcast( pos, meta_ionode_id )
|
||||
!
|
||||
DEALLOCATE( t )
|
||||
DEALLOCATE( g )
|
||||
DEALLOCATE( s )
|
||||
!
|
||||
RETURN
|
||||
!
|
||||
END SUBROUTINE broyden
|
||||
!
|
||||
! ... Broyden (rank one) optimisation - second attempt
|
||||
!
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE broyden2()
|
||||
!-----------------------------------------------------------------------
|
||||
#define DEBUG
|
||||
!
|
||||
USE control_flags, ONLY : lsmd
|
||||
USE io_files, ONLY : broy_file, iunbroy, iunpath
|
||||
USE path_variables_pre, ONLY : dim1, frozen, tangent, nim => num_of_images
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
REAL(DP), PARAMETER :: step_max = 0.6_DP
|
||||
INTEGER, PARAMETER :: broyden_ndim = 5
|
||||
!
|
||||
REAL(DP), ALLOCATABLE :: dx(:,:), df(:,:), x(:), f(:)
|
||||
REAL(DP), ALLOCATABLE :: x_last(:), f_last(:), mask(:)
|
||||
REAL(DP), ALLOCATABLE :: b(:,:), c(:), work(:)
|
||||
INTEGER, ALLOCATABLE :: iwork(:)
|
||||
!
|
||||
REAL(DP) :: x_norm, gamma0, J0, d2, d2_estimate
|
||||
LOGICAL :: exists
|
||||
INTEGER :: i, I_in, I_fin, info, j, niter
|
||||
!
|
||||
! ... starting guess for the inverse Jacobian
|
||||
!
|
||||
J0 = ds*ds
|
||||
!
|
||||
ALLOCATE( dx( dim1*nim, broyden_ndim ), df( dim1*nim, broyden_ndim ) )
|
||||
ALLOCATE( x( dim1*nim ), f( dim1*nim ) )
|
||||
ALLOCATE( x_last( dim1*nim ), f_last( dim1*nim ), mask( dim1*nim ) )
|
||||
!
|
||||
! define mask to skip frozen images
|
||||
!
|
||||
mask(:) = 0.0_DP
|
||||
DO i = 1, nim
|
||||
I_in = ( i - 1 )*dim1 + 1
|
||||
I_fin = i * dim1
|
||||
IF ( frozen(i) ) CYCLE
|
||||
mask(I_in:I_fin) = 1.0_DP
|
||||
END DO
|
||||
!
|
||||
! copy current positions and gradients in local arrays
|
||||
!
|
||||
DO i = 1, nim
|
||||
I_in = ( i - 1 )*dim1 + 1
|
||||
I_fin = i * dim1
|
||||
f(I_in:I_fin) =-grad(:,i)
|
||||
x(I_in:I_fin) = pos(:,i)
|
||||
END DO
|
||||
!
|
||||
! only meta_ionode execute this part
|
||||
!
|
||||
IF ( meta_ionode ) THEN
|
||||
d2 = DOT_PRODUCT( f(:),mask(:)*f(:) )
|
||||
#ifdef DEBUG
|
||||
WRITE (*,*) " CURRENT ACTUAL D2 = ", d2
|
||||
#endif
|
||||
!
|
||||
! ... open the file containing the broyden's history
|
||||
!
|
||||
INQUIRE( FILE = broy_file, EXIST = exists )
|
||||
!
|
||||
IF ( exists ) THEN
|
||||
!
|
||||
OPEN( UNIT = iunbroy, FILE = broy_file, STATUS = "OLD" )
|
||||
!
|
||||
READ( UNIT = iunbroy , FMT = * ) i
|
||||
!
|
||||
! ... if the number of images is changed the broyden history is
|
||||
! ... reset and the algorithm starts from scratch
|
||||
!
|
||||
exists = ( i == nim )
|
||||
!
|
||||
END IF
|
||||
!
|
||||
IF ( exists ) THEN
|
||||
!
|
||||
READ( UNIT = iunbroy , FMT = * ) niter, d2_estimate
|
||||
READ( UNIT = iunbroy , FMT = * ) df(:,:), dx(:,:)
|
||||
READ( UNIT = iunbroy , FMT = * ) f_last(:), x_last(:)
|
||||
niter = min(broyden_ndim, niter + 1)
|
||||
!
|
||||
if (d2 > 2.0_DP * d2_estimate ) then
|
||||
#ifdef DEBUG
|
||||
write (*,*) " bad D2 estimate ... reset history "
|
||||
#endif
|
||||
niter = 1
|
||||
df(:,:) = 0.0_DP
|
||||
dx(:,:) = 0.0_DP
|
||||
end if
|
||||
ELSE
|
||||
!
|
||||
df(:,:) = 0.0_DP
|
||||
dx(:,:) = 0.0_DP
|
||||
niter = 0
|
||||
!
|
||||
END IF
|
||||
CLOSE( UNIT = iunbroy )
|
||||
!
|
||||
! ... Broyden's update
|
||||
!
|
||||
! shift previous history, automatically discarding oldest iterations
|
||||
!
|
||||
DO i = broyden_ndim, 2, -1
|
||||
df(:,i) = df(:,i-1)
|
||||
dx(:,i) = dx(:,i-1)
|
||||
END DO
|
||||
!
|
||||
! and update it with last increment
|
||||
!
|
||||
IF (niter > 0 ) THEN
|
||||
df(:,1) = f(:) - f_last(:)
|
||||
dx(:,1) = x(:) - x_last(:)
|
||||
END IF
|
||||
! save for later use
|
||||
f_last(:) = f(:)
|
||||
x_last(:) = x(:)
|
||||
!
|
||||
x(:) = 0.0_DP
|
||||
IF ( niter > 0 ) THEN
|
||||
!
|
||||
ALLOCATE (b(niter,niter), c(niter), work(niter), iwork(niter))
|
||||
!
|
||||
! create the matrix and the right-hand side of the liner system
|
||||
!
|
||||
b(:,:) = 0.0_DP
|
||||
c(:) = 0.0_DP
|
||||
DO i = 1,niter
|
||||
DO j = 1,niter
|
||||
b(i,j) = DOT_PRODUCT(df(:,i),mask(:)*df(:,j))
|
||||
END DO
|
||||
c(i) = DOT_PRODUCT(f(:),mask(:)*df(:,i))
|
||||
END DO
|
||||
!
|
||||
! solve the linear system
|
||||
!
|
||||
CALL DSYTRF( 'U', niter, b, niter, iwork, work, niter, info )
|
||||
CALL errore( 'broyden', 'factorization', abs(info) )
|
||||
CALL DSYTRI( 'U', niter, b, niter, iwork, work, info )
|
||||
CALL errore( 'broyden', 'DSYTRI', abs(info) )
|
||||
FORALL( i = 1:niter, j = 1:niter, j > i ) b(j,i) = b(i,j)
|
||||
!
|
||||
! set the best correction vector and gradient
|
||||
!
|
||||
DO i = 1, niter
|
||||
gamma0 = DOT_PRODUCT( b(1:niter,i), c(1:niter) )
|
||||
call DAXPY(dim1*nim, -gamma0, dx(:,i),1, x,1)
|
||||
call DAXPY(dim1*nim, -gamma0, df(:,i),1, f,1)
|
||||
END DO
|
||||
!
|
||||
DEALLOCATE (b,c,work,iwork)
|
||||
!
|
||||
END IF
|
||||
d2 = DOT_PRODUCT( f(:), mask(:)*f(:) )
|
||||
x(:) = mask(:) * ( x(:) + J0 * f(:) )
|
||||
x_norm = norm(x)
|
||||
x(:) = x_last(:) + x(:) * min ( 1.0_DP, step_max/x_norm)
|
||||
#ifdef DEBUG
|
||||
WRITE (*,*) " ESTIMATED NEXT D2 = ", d2
|
||||
IF (x_norm > step_max) &
|
||||
WRITE (*,*) " x_norm = ", x_norm, step_max
|
||||
#endif
|
||||
!
|
||||
! ... save the file containing the history
|
||||
!
|
||||
OPEN( UNIT = iunbroy, FILE = broy_file )
|
||||
!
|
||||
WRITE( UNIT = iunbroy, FMT = * ) nim
|
||||
WRITE( UNIT = iunbroy, FMT = * ) niter, d2
|
||||
WRITE( UNIT = iunbroy, FMT = * ) df(:,:), dx(:,:)
|
||||
WRITE( UNIT = iunbroy, FMT = * ) f_last(:), x_last(:)
|
||||
!
|
||||
CLOSE( UNIT = iunbroy )
|
||||
!
|
||||
! ... copy broyden's step on the position array ...
|
||||
!
|
||||
pos(:,1:nim) = RESHAPE( x, (/ dim1, nim /) )
|
||||
!
|
||||
END IF
|
||||
!
|
||||
! ... and distribute it
|
||||
!
|
||||
CALL mp_bcast( pos, meta_ionode_id )
|
||||
!
|
||||
DEALLOCATE( df, dx, f, x, f_last, x_last, mask )
|
||||
!
|
||||
RETURN
|
||||
!
|
||||
END SUBROUTINE broyden2
|
||||
!
|
||||
END MODULE path_opt_routines_pre
|
|
@ -0,0 +1,308 @@
|
|||
!
|
||||
! Copyright (C) 2003-2006 Quantum ESPRESSO 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,
|
||||
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||
!
|
||||
!
|
||||
!---------------------------------------------------------------------------
|
||||
MODULE path_reparametrisation_pre
|
||||
!---------------------------------------------------------------------------
|
||||
!
|
||||
! ... This module contains all subroutines and functions needed for
|
||||
! ... the reparametrisation of the path in the string method
|
||||
!
|
||||
! ... Written by Carlo Sbraccia ( 2003-2006 )
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
USE io_files, ONLY : iunpath
|
||||
USE io_global, ONLY : meta_ionode, meta_ionode_id
|
||||
USE mp, ONLY : mp_bcast
|
||||
!
|
||||
USE basic_algebra_routines
|
||||
!
|
||||
PRIVATE
|
||||
!
|
||||
PUBLIC :: reparametrise, spline_interpolation
|
||||
!
|
||||
INTERFACE spline_interpolation
|
||||
!
|
||||
MODULE PROCEDURE spline_interpolation_1D, spline_interpolation_2D
|
||||
!
|
||||
END INTERFACE
|
||||
!
|
||||
CONTAINS
|
||||
!
|
||||
! ... reparametrisation routines in real space
|
||||
!
|
||||
!------------------------------------------------------------------------
|
||||
SUBROUTINE reparametrise()
|
||||
!------------------------------------------------------------------------
|
||||
!
|
||||
USE path_variables_pre, ONLY : pos
|
||||
USE path_variables_pre, ONLY : nim => num_of_images
|
||||
USE path_variables_pre, ONLY : climbing
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
INTEGER :: i, ni, nf
|
||||
!
|
||||
!
|
||||
IF ( meta_ionode ) THEN
|
||||
!
|
||||
IF ( ANY( climbing(:) ) ) THEN
|
||||
!
|
||||
ni = 1
|
||||
!
|
||||
DO i = 2, nim
|
||||
!
|
||||
IF ( .NOT. climbing(i) ) CYCLE
|
||||
!
|
||||
nf = i
|
||||
!
|
||||
CALL spline_interpolation( pos, ni, nf )
|
||||
!
|
||||
ni = nf
|
||||
!
|
||||
END DO
|
||||
!
|
||||
nf = nim
|
||||
!
|
||||
CALL spline_interpolation( pos, ni, nf )
|
||||
!
|
||||
ELSE
|
||||
!
|
||||
ni = 1
|
||||
nf = nim
|
||||
!
|
||||
CALL spline_interpolation( pos, ni, nf )
|
||||
!
|
||||
END IF
|
||||
!
|
||||
END IF
|
||||
!
|
||||
CALL mp_bcast( pos, meta_ionode_id )
|
||||
!
|
||||
RETURN
|
||||
!
|
||||
END SUBROUTINE reparametrise
|
||||
!
|
||||
!--------------------------------------------------------------------
|
||||
SUBROUTINE spline_interpolation_1D( vec, ni, nf, nim )
|
||||
!--------------------------------------------------------------------
|
||||
!
|
||||
USE splinelib, ONLY : dosplineint
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
REAL(DP), INTENT(INOUT) :: vec(:)
|
||||
INTEGER, INTENT(IN) :: ni, nf
|
||||
INTEGER, INTENT(IN), OPTIONAL :: nim
|
||||
!
|
||||
INTEGER :: i, j
|
||||
INTEGER :: nio, nfo
|
||||
REAL(DP) :: delta, length
|
||||
REAL(DP), ALLOCATABLE :: new_vec(:)
|
||||
REAL(DP), ALLOCATABLE :: old_mesh(:), new_mesh(:)
|
||||
!
|
||||
!
|
||||
IF ( PRESENT( nim ) ) THEN
|
||||
!
|
||||
nio = 1
|
||||
nfo = nim
|
||||
!
|
||||
ELSE
|
||||
!
|
||||
nio = ni
|
||||
nfo = nf
|
||||
!
|
||||
END IF
|
||||
!
|
||||
! ... cubic spline interpolation
|
||||
!
|
||||
ALLOCATE( new_vec( ni:nf ) )
|
||||
!
|
||||
ALLOCATE( old_mesh( nio:nfo ) )
|
||||
ALLOCATE( new_mesh( ni:nf ) )
|
||||
!
|
||||
old_mesh(:) = 0.0_DP
|
||||
new_mesh(:) = 0.0_DP
|
||||
!
|
||||
DO i = nio, nfo - 1
|
||||
!
|
||||
old_mesh(i+1) = old_mesh(i) + ABS( vec(i+1) - vec(i) )
|
||||
!
|
||||
END DO
|
||||
!
|
||||
length = old_mesh(nfo)
|
||||
!
|
||||
delta = length / DBLE( nf - ni )
|
||||
!
|
||||
DO j = 0, nf - ni
|
||||
!
|
||||
new_mesh(j+ni) = DBLE(j) * delta
|
||||
!
|
||||
END DO
|
||||
!
|
||||
old_mesh(:) = old_mesh(:) / length
|
||||
new_mesh(:) = new_mesh(:) / length
|
||||
!
|
||||
CALL dosplineint( old_mesh(:), vec(nio:nfo), new_mesh(:), new_vec(:) )
|
||||
!
|
||||
vec(ni:nf) = new_vec(:)
|
||||
!
|
||||
DEALLOCATE( new_vec, old_mesh, new_mesh )
|
||||
!
|
||||
RETURN
|
||||
!
|
||||
END SUBROUTINE spline_interpolation_1D
|
||||
!
|
||||
!--------------------------------------------------------------------
|
||||
SUBROUTINE spline_interpolation_2D( vec, ni, nf, nim )
|
||||
!--------------------------------------------------------------------
|
||||
!
|
||||
USE splinelib, ONLY : dosplineint
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
REAL(DP), INTENT(INOUT) :: vec(:,:)
|
||||
INTEGER, INTENT(IN) :: ni, nf
|
||||
INTEGER, INTENT(IN), OPTIONAL :: nim
|
||||
!
|
||||
INTEGER :: i, j
|
||||
INTEGER :: nio, nfo
|
||||
INTEGER :: dim1
|
||||
REAL(DP) :: delta, length
|
||||
REAL(DP), ALLOCATABLE :: new_vec(:,:)
|
||||
REAL(DP), ALLOCATABLE :: old_mesh(:), new_mesh(:)
|
||||
!
|
||||
!
|
||||
dim1 = SIZE( vec, 1 )
|
||||
!
|
||||
IF ( PRESENT( nim ) ) THEN
|
||||
!
|
||||
nio = 1
|
||||
nfo = nim
|
||||
!
|
||||
ELSE
|
||||
!
|
||||
nio = ni
|
||||
nfo = nf
|
||||
!
|
||||
END IF
|
||||
!
|
||||
! ... cubic spline interpolation
|
||||
!
|
||||
ALLOCATE( new_vec( dim1, ni:nf ) )
|
||||
!
|
||||
ALLOCATE( old_mesh( nio:nfo ) )
|
||||
ALLOCATE( new_mesh( ni:nf ) )
|
||||
!
|
||||
old_mesh(:) = 0.0_DP
|
||||
new_mesh(:) = 0.0_DP
|
||||
!
|
||||
DO i = nio, nfo - 1
|
||||
!
|
||||
old_mesh(i+1) = old_mesh(i) + norm( vec(:,i+1) - vec(:,i) )
|
||||
!
|
||||
END DO
|
||||
!
|
||||
length = old_mesh(nfo)
|
||||
!
|
||||
delta = length / DBLE( nf - ni )
|
||||
!
|
||||
DO j = 0, nf - ni
|
||||
!
|
||||
new_mesh(j+ni) = DBLE(j) * delta
|
||||
!
|
||||
END DO
|
||||
!
|
||||
old_mesh(:) = old_mesh(:) / length
|
||||
new_mesh(:) = new_mesh(:) / length
|
||||
!
|
||||
CALL dosplineint( old_mesh(:), vec(:,nio:nfo), new_mesh(:), new_vec(:,:) )
|
||||
!
|
||||
vec(:,ni:nf) = new_vec(:,:)
|
||||
!
|
||||
DEALLOCATE( new_vec, old_mesh, new_mesh )
|
||||
!
|
||||
RETURN
|
||||
!
|
||||
END SUBROUTINE spline_interpolation_2D
|
||||
!
|
||||
!--------------------------------------------------------------------
|
||||
SUBROUTINE cubic_interpolation( ni, nf )
|
||||
!--------------------------------------------------------------------
|
||||
!
|
||||
USE path_variables_pre, ONLY : dim1, pos
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
INTEGER, INTENT(IN) :: ni, nf
|
||||
!
|
||||
INTEGER :: i, j
|
||||
REAL(DP) :: r, delta, x
|
||||
REAL(DP), ALLOCATABLE :: a(:,:), b(:,:), c(:,:), d(:,:), t(:,:), s(:)
|
||||
!
|
||||
ALLOCATE( a( dim1, ni:nf-1 ) )
|
||||
ALLOCATE( b( dim1, ni:nf-1 ) )
|
||||
ALLOCATE( c( dim1, ni:nf-1 ) )
|
||||
ALLOCATE( d( dim1, ni:nf-1 ) )
|
||||
ALLOCATE( t( dim1, ni:nf ) )
|
||||
ALLOCATE( s( ni:nf ) )
|
||||
!
|
||||
t(:,ni) = pos(:,ni+1) - pos(:,ni)
|
||||
t(:,nf) = pos(:,nf) - pos(:,nf-1)
|
||||
!
|
||||
DO i = ni+1, nf - 1
|
||||
!
|
||||
t(:,i) = ( pos(:,i+1) - pos(:,i-1) ) / 2.0_DP
|
||||
!
|
||||
END DO
|
||||
!
|
||||
s(ni) = 0.0_DP
|
||||
!
|
||||
DO i = ni, nf - 1
|
||||
!
|
||||
r = norm( pos(:,i+1) - pos(:,i) )
|
||||
!
|
||||
s(i+1) = s(i) + r
|
||||
!
|
||||
! ... cubic interpolation
|
||||
!
|
||||
a(:,i) = 2.0_DP *( pos(:,i) - pos(:,i+1) ) / r**3 + &
|
||||
( t(:,i) + t(:,i+1) ) / r**2
|
||||
!
|
||||
b(:,i) = 3.0_DP *( pos(:,i+1) - pos(:,i) ) / r**2 - &
|
||||
( 2.0_DP*t(:,i) + t(:,i+1) ) / r
|
||||
!
|
||||
c(:,i) = t(:,i)
|
||||
!
|
||||
d(:,i) = pos(:,i)
|
||||
!
|
||||
END DO
|
||||
!
|
||||
i = ni
|
||||
!
|
||||
delta = s(nf) / DBLE( nf - ni )
|
||||
!
|
||||
DO j = ni, nf
|
||||
!
|
||||
r = DBLE( j - ni ) * delta
|
||||
!
|
||||
IF ( r >= s(i+1) .AND. i < nf - 1 ) i = i + 1
|
||||
!
|
||||
x = r - s(i)
|
||||
!
|
||||
pos(:,j) = a(:,i)*x**3 + b(:,i)*x**2 + c(:,i)*x + d(:,i)
|
||||
!
|
||||
END DO
|
||||
!
|
||||
DEALLOCATE( a, b, c, d, t, s )
|
||||
!
|
||||
RETURN
|
||||
!
|
||||
END SUBROUTINE cubic_interpolation
|
||||
!
|
||||
END MODULE path_reparametrisation_pre
|
|
@ -0,0 +1,153 @@
|
|||
!
|
||||
! Copyright (C) 2003-2006 Quantum ESPRESSO 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,
|
||||
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||
!
|
||||
!--------------------------------------------------------------------------
|
||||
MODULE path_variables_pre
|
||||
!---------------------------------------------------------------------------
|
||||
!
|
||||
! ... This module contains all variables needed by path optimisations
|
||||
!
|
||||
! ... Written by Carlo Sbraccia ( 2003-2006 )
|
||||
!
|
||||
USE kinds, ONLY : DP
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
SAVE
|
||||
!
|
||||
! ... "general" variables :
|
||||
!
|
||||
LOGICAL :: &
|
||||
conv_path ! .TRUE. when "path" convergence has been
|
||||
! achieved
|
||||
LOGICAL :: &
|
||||
first_last_opt, &! if .TRUE. the first and the last image
|
||||
! are optimised too.
|
||||
use_masses, &! if .TRUE. mass weighted coordinates are
|
||||
! used
|
||||
fixed_tan, &! if. TRUE. the projection is done using the
|
||||
! tangent of the average path
|
||||
use_freezing, &! if .TRUE. images are optimised according
|
||||
! to their error (see frozen array)
|
||||
tune_load_balance ! if .TRUE. the load balance for image
|
||||
! parallelisation is tuned at
|
||||
! runtime
|
||||
INTEGER :: &
|
||||
dim1, &! dimension of the configuration space
|
||||
num_of_images, &! number of images
|
||||
deg_of_freedom, &! number of degrees of freedom
|
||||
! ( dim1 - #( of fixed coordinates ) )
|
||||
pending_image ! last image for which scf has not been
|
||||
! achieved
|
||||
REAL(DP) :: &
|
||||
ds, &! the optimization step
|
||||
path_thr, &! convergence threshold
|
||||
temp_req, &! required temperature
|
||||
activation_energy, &! forward activatation energy
|
||||
err_max, &! the largest error
|
||||
path_length ! length of the path
|
||||
LOGICAL :: &
|
||||
lsteep_des = .FALSE., &! .TRUE. if opt_scheme = "sd"
|
||||
lquick_min = .FALSE., &! .TRUE. if opt_scheme = "quick-min"
|
||||
lbroyden = .FALSE., &! .TRUE. if opt_scheme = "broyden"
|
||||
lbroyden2 = .FALSE., &! .TRUE. if opt_scheme = "broyden2"
|
||||
llangevin = .FALSE. ! .TRUE. if opt_scheme = "langevin"
|
||||
INTEGER :: &
|
||||
istep_path, &! iteration in the optimization procedure
|
||||
nstep_path ! maximum number of iterations
|
||||
!
|
||||
! ... "general" real space arrays
|
||||
!
|
||||
REAL(DP), ALLOCATABLE :: &
|
||||
pes(:), &! the potential enrgy along the path
|
||||
error(:) ! the error from the true MEP
|
||||
REAL(DP), ALLOCATABLE :: &
|
||||
pos(:,:), &! reaction path
|
||||
grad_pes(:,:), &! gradients acting on the path
|
||||
tangent(:,:) ! tangent to the path
|
||||
LOGICAL, ALLOCATABLE :: &
|
||||
frozen(:) ! .TRUE. if the image or mode has not
|
||||
! to be optimized
|
||||
!
|
||||
! ... "neb specific" variables :
|
||||
!
|
||||
LOGICAL, ALLOCATABLE :: &
|
||||
climbing(:) ! .TRUE. if the image is required to climb
|
||||
CHARACTER(LEN=20) :: &
|
||||
CI_scheme ! Climbing Image scheme
|
||||
INTEGER :: &
|
||||
Emax_index ! index of the image with the highest energy
|
||||
!
|
||||
REAL (DP) :: &
|
||||
k_max, &!
|
||||
k_min, &!
|
||||
Emax, &!
|
||||
Emin !
|
||||
!
|
||||
! ... real space arrays
|
||||
!
|
||||
REAL(DP), ALLOCATABLE :: &
|
||||
elastic_grad(:), &! elastic part of the gradients
|
||||
mass(:), &! atomic masses
|
||||
k(:) ! elastic constants
|
||||
REAL(DP), ALLOCATABLE :: &
|
||||
posold(:,:), &! old positions (for the quick-min)
|
||||
grad(:,:), &!
|
||||
lang(:,:) ! langevin random force
|
||||
!
|
||||
CONTAINS
|
||||
!
|
||||
!----------------------------------------------------------------------
|
||||
SUBROUTINE path_allocation()
|
||||
!----------------------------------------------------------------------
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
ALLOCATE( pos( dim1, num_of_images ) )
|
||||
!
|
||||
ALLOCATE( posold( dim1, num_of_images ) )
|
||||
ALLOCATE( grad( dim1, num_of_images ) )
|
||||
ALLOCATE( grad_pes( dim1, num_of_images ) )
|
||||
ALLOCATE( tangent( dim1, num_of_images ) )
|
||||
!
|
||||
ALLOCATE( pes( num_of_images ) )
|
||||
ALLOCATE( k( num_of_images ) )
|
||||
ALLOCATE( error( num_of_images ) )
|
||||
ALLOCATE( climbing( num_of_images ) )
|
||||
ALLOCATE( frozen( num_of_images ) )
|
||||
!
|
||||
ALLOCATE( mass( dim1 ) )
|
||||
ALLOCATE( elastic_grad( dim1 ) )
|
||||
!
|
||||
ALLOCATE( lang( dim1, num_of_images ) )
|
||||
!
|
||||
END SUBROUTINE path_allocation
|
||||
!
|
||||
!
|
||||
!----------------------------------------------------------------------
|
||||
SUBROUTINE path_deallocation()
|
||||
!----------------------------------------------------------------------
|
||||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
IF ( ALLOCATED( pos ) ) DEALLOCATE( pos )
|
||||
IF ( ALLOCATED( posold ) ) DEALLOCATE( posold )
|
||||
IF ( ALLOCATED( grad ) ) DEALLOCATE( grad )
|
||||
IF ( ALLOCATED( pes ) ) DEALLOCATE( pes )
|
||||
IF ( ALLOCATED( grad_pes ) ) DEALLOCATE( grad_pes )
|
||||
IF ( ALLOCATED( k ) ) DEALLOCATE( k )
|
||||
IF ( ALLOCATED( mass ) ) DEALLOCATE( mass )
|
||||
IF ( ALLOCATED( elastic_grad ) ) DEALLOCATE( elastic_grad )
|
||||
IF ( ALLOCATED( tangent ) ) DEALLOCATE( tangent )
|
||||
IF ( ALLOCATED( error ) ) DEALLOCATE( error )
|
||||
IF ( ALLOCATED( climbing ) ) DEALLOCATE( climbing )
|
||||
IF ( ALLOCATED( frozen ) ) DEALLOCATE( frozen )
|
||||
IF ( ALLOCATED( lang ) ) DEALLOCATE( lang )
|
||||
!
|
||||
END SUBROUTINE path_deallocation
|
||||
!
|
||||
END MODULE path_variables_pre
|
Loading…
Reference in New Issue