comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Dimer Method: first version working with SD. Cleaning and reorganization of the

Browse source 
CG section, gopt_para_type. Getting rid of useless variables. Introduced the
possibility to have different kind of LS (still not fully implemented). Added
1 regtest for Dimer Method with SD. Common message to all geo_opt methods in case
of non-convergence. Prettifying.


svn-origin-rev: 6833
This commit is contained in:
Teodoro Laino 2008-01-20 23:25:00 +00:00
parent b166b67f0b
commit ea42cd70aa
25 changed files with 519 additions and 284 deletions
Download patch file Download diff file Expand all files Collapse all files

23
src/bfgs_optimizer.F
View file

@ -39,6 +39,7 @@ MODULE bfgs_optimizer
USE global_types, ONLY: global_environment_type
USE gopt_f_methods, ONLY: check_converg,&
print_geo_opt_header,&
print_geo_opt_nc,&
write_cycle_infos,&
write_final_info,&
write_geo_traj
@ -221,24 +222,10 @@ CONTAINS
its = 0
CALL write_cycle_infos(output_unit,its,etot,wildcard=" BFGS")
DO
its = its + 1
DO its = 1, maxiter
gopt_env%its = its
CALL cp_iterate(logger%iter_info,iter_nr=its,error=error)
IF (its > maxiter ) THEN
IF (output_unit > 0) THEN
WRITE(UNIT=output_unit,FMT="(/,T2,A)") REPEAT("*",79)
WRITE (UNIT=output_unit,FMT="(T2,A,T20,A,T78,A)")&
"***","MAXIMUM NUMBER OF OPTIMIZATION STEPS REACHED","***"
WRITE (UNIT=output_unit,FMT="(T2,A,T27,A,T78,A)")&
"***","EXITING GEOMETRY OPTIMIZATION","***"
WRITE(UNIT=output_unit,FMT="(T2,A)") REPEAT("*",79)
CALL m_flush(output_unit)
END IF
EXIT
END IF
! check for an external exit command
CALL external_control(should_stop,"GEO",globenv,error)
IF(should_stop) EXIT
@ -331,9 +318,13 @@ CONTAINS
IF (etot < emin) emin = etot
CALL update_trust_rad(rat,rad,step,ediff)
CALL check_converg(ndf,dr,g,output_unit,conv,gopt_param)
IF (conv) EXIT
IF (conv.OR.(its==maxiter)) EXIT
END DO
IF(its == maxiter .AND. (.NOT.conv))THEN
CALL print_geo_opt_nc(gopt_env, output_unit)
END IF
! Write final particle information, if converged
CALL cp_iterate(logger%iter_info,last=.TRUE.,iter_nr=its,error=error)
CALL write_final_info(conv, its, gopt_env, x0, gopt_env%force_env%para_env%mepos,&

50
src/cg_optimizer.F
View file

@ -39,6 +39,7 @@ MODULE cg_optimizer
USE gopt_f_methods, ONLY: check_converg,&
check_rot_conv,&
print_geo_opt_header,&
print_geo_opt_nc,&
write_cycle_infos,&
write_final_info,&
write_geo_traj,&
@ -49,7 +50,6 @@ MODULE cg_optimizer
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type
USE kinds, ONLY: dp
USE machine, ONLY: m_flush
USE particle_list_types, ONLY: particle_list_type
USE particle_types, ONLY: particle_type
USE termination, ONLY: external_control
@ -177,13 +177,12 @@ CONTAINS
routineP = moduleN//':'//routineN
CHARACTER(LEN=5) :: wildcard
INTEGER :: brent_max_iter, handle, its, &
max_steep_steps, maxiter, stat
INTEGER :: handle, its, max_steep_steps, &
maxiter, stat
LOGICAL :: conv, failure, &
Fletcher_Reeves, should_stop
REAL(KIND=dp) :: brack_limit, brent_tol, emin, &
eold, opt_energy, res_lim, &
step
REAL(KIND=dp) :: emin, eold, opt_energy, &
res_lim
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xold
REAL(KIND=dp), DIMENSION(:), POINTER :: g, h, xi
TYPE(cp_logger_type), POINTER :: logger
@ -191,26 +190,17 @@ CONTAINS
TYPE(section_vals_type), POINTER :: root_section
CALL timeset(routineN,handle)
failure = .FALSE.
NULLIFY(logger,g, h, xi)
root_section => force_env%root_section
para_env => force_env%para_env
logger => cp_error_get_logger(error)
IF (gopt_env%dimer_rotation) CALL cp_unimplemented_error(fromWhere=routineP, &
message="DIMER method not yet working with CG.", &
error=error, error_level=cp_failure_level)
conv = .FALSE.
failure = .FALSE.
IF (.NOT.failure) THEN
maxiter = gopt_param%max_iter
max_steep_steps = gopt_param%max_steep_steps
Fletcher_Reeves = gopt_param%Fletcher_Reeves
brent_tol = gopt_param%brent_tol
brent_max_iter = gopt_param%brent_max_iter
brack_limit = gopt_param%brack_limit
res_lim = gopt_param%restart_limit
step = gopt_param%initial_step
ALLOCATE(g (SIZE(x0)), stat=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(h (SIZE(x0)), stat=stat)
@ -224,10 +214,10 @@ CONTAINS
CALL cp_eval_at(gopt_env, x0, opt_energy, xi, gopt_env%force_env%para_env%mepos,&
gopt_env%force_env%para_env, error)
g = -xi
h = g
xi = h
emin = HUGE(0.0_dp)
g = -xi
h = g
xi = h
emin = HUGE(0.0_dp)
! Main Loop
wildcard = " SD"
IF (.NOT.gopt_env%dimer_rotation) THEN
@ -242,8 +232,9 @@ CONTAINS
xold = x0
eold = opt_energy
emin = MIN(emin, opt_energy)
CALL cg_linmin(gopt_env, x0, xi, opt_energy, step, output_unit, gopt_param, globenv, error)
! Line minimization
CALL cg_linmin(gopt_env, x0, xi, opt_energy, output_unit, gopt_param, globenv, error)
! Check for an external exit command
CALL external_control(should_stop,"GEO",globenv,error)
@ -281,18 +272,7 @@ CONTAINS
xi = h
END DO
IF(its == maxiter .AND. (.NOT.conv))THEN
IF(output_unit>0)THEN
WRITE(UNIT=output_unit,FMT="(/,T2,A)")&
"*** MAXIMUM NUMBER OF OPTIMIZATION STEPS REACHED ***"
IF (.NOT.gopt_env%dimer_rotation) THEN
WRITE(UNIT=output_unit,FMT="(T2,A)")&
"*** EXITING GEOMETRY OPTIMIZATION ***"
ELSE
WRITE(UNIT=output_unit,FMT="(T2,A)")&
"*** EXITING ROTATION OPTIMIZATION ***"
END IF
CALL m_flush(output_unit)
END IF
CALL print_geo_opt_nc(gopt_env, output_unit)
END IF
! Write final particle information, if converged
CALL cp_iterate(logger%iter_info,last=.TRUE.,iter_nr=its,error=error)

210
src/cg_utils.F
View file

@ -29,8 +29,9 @@ MODULE cg_utils
USE gopt_param_types, ONLY: gopt_param_type
USE input_constants, ONLY: default_minimization_method_id,&
default_ts_method_id,&
do_rotate_dimer,&
do_translate_dimer
ls_2pnt,&
ls_fit,&
ls_gold
USE kinds, ONLY: dp
USE mathconstants, ONLY: pi
USE memory_utilities, ONLY: reallocate
@ -71,12 +72,12 @@ CONTAINS
!! 10.2005 created [tlaino]
!!
!!*** **********************************************************************
SUBROUTINE cg_linmin(gopt_env, xvec, xi, opt_energy, step, output_unit,&
gopt_param, globenv, error)
SUBROUTINE cg_linmin(gopt_env, xvec, xi, opt_energy, output_unit, gopt_param,&
globenv, error)
TYPE(gopt_f_type), POINTER :: gopt_env
REAL(KIND=dp), DIMENSION(:), POINTER :: xvec, xi
REAL(KIND=dp) :: opt_energy, step
REAL(KIND=dp), INTENT(INOUT) :: opt_energy
INTEGER :: output_unit
TYPE(gopt_param_type), POINTER :: gopt_param
TYPE(global_environment_type), POINTER :: globenv
@ -85,45 +86,56 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'cg_linmin', &
routineP = moduleN//':'//routineN
INTEGER :: brent_max_iter, handle
INTEGER :: handle
LOGICAL :: failure
REAL(KIND=dp) :: brack_limit, brent_tol
failure = .FALSE.
CALL timeset(routineN,handle)
gopt_env%do_line_search = .TRUE.
SELECT CASE(gopt_env%type_id)
CASE(default_minimization_method_id)
IF (gopt_param%use_fit_line_search) THEN
CALL linminE(gopt_env, xvec, xi, opt_energy, gopt_param%brent_tol,&
gopt_param%brent_max_iter, gopt_param%brack_limit, &
step, output_unit, gopt_param, globenv, error)
ELSE
CALL linmin(gopt_env, xvec, xi, opt_energy, gopt_param%brent_tol,&
gopt_param%brent_max_iter, gopt_param%brack_limit, &
step, output_unit, globenv, error)
END IF
CASE(default_ts_method_id)
SELECT CASE(gopt_env%dimer_env%status_dimer)
CASE(do_rotate_dimer)
CALL rotmin(gopt_env, gopt_env%dimer_env, xvec, xi, error)
CASE(do_translate_dimer)
! To Fix
STOP
SELECT CASE(gopt_param%cg_ls%type_id)
CASE(ls_fit)
CALL linmin_fit(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brent_tol,&
gopt_param%cg_ls%brent_max_iter, gopt_param%cg_ls%brack_limit, &
gopt_param%cg_ls%initial_step, output_unit, gopt_param, globenv, error)
CASE(ls_gold)
CALL linmin_gold(gopt_env, xvec, xi, opt_energy, gopt_param%cg_ls%brent_tol,&
gopt_param%cg_ls%brent_max_iter, gopt_param%cg_ls%brack_limit, &
gopt_param%cg_ls%initial_step, output_unit, globenv, error)
CASE DEFAULT
CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure)
CALL cp_unimplemented_error(fromWhere=routineP, &
message="Line Search type not yet implemented in CG.", &
error=error, error_level=cp_failure_level)
END SELECT
CASE(default_ts_method_id)
SELECT CASE(gopt_param%cg_ls%type_id)
CASE(ls_2pnt)
IF (gopt_env%dimer_rotation) THEN
CALL rotmin_2pnt(gopt_env, gopt_env%dimer_env, xvec, xi, opt_energy, gopt_param,&
error)
ELSE
CALL tslmin_2pnt(gopt_env, gopt_env%dimer_env, xvec, xi, opt_energy, gopt_param,&
error)
END IF
CASE DEFAULT
CALL cp_unimplemented_error(fromWhere=routineP, &
message="Line Search type not yet implemented in CG for TS search.", &
error=error, error_level=cp_failure_level)
END SELECT
END SELECT
gopt_env%do_line_search = .FALSE.
CALL timestop(handle)
END SUBROUTINE cg_linmin
!!****f* cg_utils/rotmin [1.0] *
!!****f* cg_utils/tslmin_2pnt [1.0] *
!!
!! NAME
!! rotmin
!! tslmin_2pnt
!!
!! FUNCTION
!! Rotational minimization for the Dimer MEthod
!! Translational minimization for the Dimer Method - 2pnt LS
!!
!! NOTES
!!
@ -139,13 +151,85 @@ CONTAINS
!!
!!
!!*** **********************************************************************
SUBROUTINE rotmin(gopt_env, dimer_env, x0, theta, error)
SUBROUTINE tslmin_2pnt(gopt_env, dimer_env, x0, tls_vec, opt_energy, gopt_param, error)
TYPE(gopt_f_type), POINTER :: gopt_env
TYPE(dimer_env_type), POINTER :: dimer_env
REAL(KIND=dp), DIMENSION(:), POINTER :: x0, tls_vec
REAL(KIND=dp), INTENT(INOUT) :: opt_energy
TYPE(gopt_param_type), POINTER :: gopt_param
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'tslmin_2pnt', &
routineP = moduleN//':'//routineN
INTEGER :: handle, stat
LOGICAL :: failure
REAL(KIND=dp) :: dx, dx_min, norm_tls_vec, &
opt_energy2
REAL(KIND=dp), DIMENSION(:), POINTER :: tls_norm
failure = .FALSE.
CALL timeset(routineN,handle)
norm_tls_vec = SQRT(DOT_PRODUCT(tls_vec,tls_vec))
IF (norm_tls_vec/=0.0_dp) THEN
ALLOCATE(tls_norm(SIZE(tls_vec)),stat=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
tls_norm = tls_vec / norm_tls_vec
dimer_env%tsl%tls_vec => tls_norm
dx = 0.01_dp
x0 = x0 + dx * tls_norm
CALL cp_eval_at(gopt_env, x0, opt_energy2, master=gopt_env%force_env%para_env%mepos,&
para_env=gopt_env%force_env%para_env, error=error)
dx_min = -opt_energy/(opt_energy2-opt_energy)*dx -dx
x0 = x0 + dx_min * tls_norm
! To remove later..
CALL cp_eval_at(gopt_env, x0, opt_energy, master=gopt_env%force_env%para_env%mepos,&
para_env=gopt_env%force_env%para_env, error=error)
DEALLOCATE(tls_norm,stat=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
ELSE
! Do Nothing, since.. if the effective force is 0 means that we are already
! in the saddle point..
END IF
CALL timestop(handle)
END SUBROUTINE tslmin_2pnt
!!****f* cg_utils/rotmin_2pnt [1.0] *
!!
!! NAME
!! rotmin_2pnt
!!
!! FUNCTION
!! Rotational minimization for the Dimer Method - 2 pnt LS
!!
!! NOTES
!!
!!
!! INPUTS
!! - error: variable to control error logging, stopping,...
!! see module cp_error_handling
!!
!! AUTHOR
!! Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
!!
!! MODIFICATION HISTORY
!!
!!
!!*** **********************************************************************
SUBROUTINE rotmin_2pnt(gopt_env, dimer_env, x0, theta, opt_energy, gopt_param, error)
TYPE(gopt_f_type), POINTER :: gopt_env
TYPE(dimer_env_type), POINTER :: dimer_env
REAL(KIND=dp), DIMENSION(:), POINTER :: x0, theta
REAL(KIND=dp), INTENT(INOUT) :: opt_energy
TYPE(gopt_param_type), POINTER :: gopt_param
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'rotmin', &
CHARACTER(len=*), PARAMETER :: routineN = 'rotmin_2pnt', &
routineP = moduleN//':'//routineN
INTEGER :: handle, stat
@ -153,26 +237,27 @@ CONTAINS
REAL(KIND=dp) :: a0, a1, angle, b1, &
curvature0, curvature1, &
curvature2, dCdp, f
REAL(KIND=dp), DIMENSION(:), POINTER :: nvec
REAL(KIND=dp), DIMENSION(:), POINTER :: work
failure = .FALSE.
CALL timeset(routineN,handle)
ALLOCATE(nvec(SIZE(dimer_env%nvec)),stat=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
nvec = dimer_env%nvec
curvature0 = dimer_env%curvature
dCdp = dimer_env%dCdp
curvature0 = dimer_env%rot%curvature
dCdp = dimer_env%rot%dCdp
b1 = 0.5_dp*dCdp
angle = -0.5_dp*ATAN(dCdp/(2.0_dp*ABS(curvature0)))
dimer_env%angle1 = angle
IF (angle>dimer_env%angle_tol) THEN
dimer_env%rot%angle1 = angle
IF (angle>dimer_env%rot%angle_tol) THEN
! Allocate a local work array
ALLOCATE(work(SIZE(dimer_env%nvec)),stat=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
work = dimer_env%nvec
! Rotating the dimer of dtheta degrees
CALL rotate_dimer(dimer_env%nvec,theta,angle,error=error)
! Re-compute energy, gradients and rotation vector for new R1
CALL cp_eval_at(gopt_env, x0, f, master=gopt_env%force_env%para_env%mepos,&
para_env=gopt_env%force_env%para_env, error=error)
curvature1 = dimer_env%curvature
curvature1 = dimer_env%rot%curvature
a1 = (curvature0 - curvature1 + b1 * SIN(2.0_dp*angle))/(1.0_dp - COS(2.0_dp*angle))
a0 = 2.0_dp*(curvature0-a1)
angle = 0.5_dp*ATAN(b1/a1)
@ -181,26 +266,37 @@ CONTAINS
angle = angle + pi/2.0_dp
curvature2 = a0/2.0_dp + a1 * COS(2.0_dp*angle)+b1 * SIN(2.0_dp*angle)
END IF
dimer_env%angle2 = angle
dimer_env%curvature = curvature2
dimer_env%rot%angle2 = angle
dimer_env%rot%curvature = curvature2
! Rotating the dimer the optimized (in plane) vector position
dimer_env%nvec = nvec
dimer_env%nvec = work
CALL rotate_dimer(dimer_env%nvec,theta,angle,error=error)
! Evaluate (by interpolation) the norm of the rotational force in the
! minimum of the rotational search (this is for print-out only)
work = dimer_env%rot%g1
work = SIN(dimer_env%rot%angle1-dimer_env%rot%angle2)/SIN(dimer_env%rot%angle1)*dimer_env%rot%g1+&
SIN(dimer_env%rot%angle2)/SIN(dimer_env%rot%angle1)*dimer_env%rot%g1p+&
(1.0_dp-COS(dimer_env%rot%angle2)-SIN(dimer_env%rot%angle2)*TAN(dimer_env%rot%angle1/2.0_dp))*&
dimer_env%rot%g0
work = -2.0_dp*(work-dimer_env%rot%g0)
work = work - DOT_PRODUCT(work, dimer_env%nvec) * dimer_env%nvec
opt_energy = SQRT(DOT_PRODUCT(work, work))
DEALLOCATE(work,stat=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
DEALLOCATE(nvec,stat=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
dimer_env%angle2 = angle
dimer_env%rot%angle2 = angle
CALL timestop(handle)
END SUBROUTINE rotmin
END SUBROUTINE rotmin_2pnt
!!****f* cg_utils/LinMinE [1.0] *
!!****f* cg_utils/linmin_fit [1.0] *
!!
!! NAME
!! LinMinE
!! linmin_fit
!!
!! FUNCTION
!! Line Minimization routine Enhanced
!! Line Minimization - Fit
!!
!! NOTES
!! Given as input the vector XVEC and XI, finds the scalar
@ -218,7 +314,7 @@ CONTAINS
!! 10.2005 created [tlaino]
!!
!!*** **********************************************************************
SUBROUTINE LinMinE(gopt_env, xvec, xi, opt_energy, brent_tol, brent_max_iter, &
SUBROUTINE linmin_fit(gopt_env, xvec, xi, opt_energy, brent_tol, brent_max_iter, &
brack_limit, step, output_unit, gopt_param, globenv, error)
TYPE(gopt_f_type), POINTER :: gopt_env
REAL(KIND=dp), DIMENSION(:), POINTER :: xvec, xi
@ -230,7 +326,7 @@ CONTAINS
TYPE(global_environment_type), POINTER :: globenv
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'LinMinE', &
CHARACTER(len=*), PARAMETER :: routineN = 'linmin_fit', &
routineP = moduleN//':'//routineN
INTEGER :: handle, loc_iter, odim, stat
@ -298,15 +394,15 @@ CONTAINS
END IF
CALL timestop(handle)
END SUBROUTINE LinMinE
END SUBROUTINE linmin_fit
!!****f* cg_utils/LinMin [1.0] *
!!****f* cg_utils/linmin_gold [1.0] *
!!
!! NAME
!! LinMin
!! linmin_gold
!!
!! FUNCTION
!! Line Minimization routine
!! Line Minimization routine - GOLD
!!
!! NOTES
!! Given as input the vector XVEC and XI, finds the scalar
@ -324,7 +420,7 @@ CONTAINS
!! 10.2005 created [tlaino]
!!
!!*** **********************************************************************
SUBROUTINE LinMin(gopt_env, xvec, xi, opt_energy, brent_tol, brent_max_iter, &
SUBROUTINE linmin_gold(gopt_env, xvec, xi, opt_energy, brent_tol, brent_max_iter, &
brack_limit, step, output_unit, globenv, error)
TYPE(gopt_f_type), POINTER :: gopt_env
REAL(KIND=dp), DIMENSION(:), POINTER :: xvec, xi
@ -335,7 +431,7 @@ CONTAINS
TYPE(global_environment_type), POINTER :: globenv
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'LinMin', &
CHARACTER(len=*), PARAMETER :: routineN = 'linmin_gold', &
routineP = moduleN//':'//routineN
INTEGER :: handle, stat
@ -372,7 +468,7 @@ CONTAINS
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
CALL timestop(handle)
END SUBROUTINE LinMin
END SUBROUTINE linmin_gold
!!****f* cg_utils/cg_mnbrak [1.0] *
!!

8
src/cp_lbfgs_geo.F
View file

@ -40,6 +40,7 @@ MODULE cp_lbfgs_geo
force_env_type
USE global_types, ONLY: global_environment_type
USE gopt_f_methods, ONLY: print_geo_opt_header,&
print_geo_opt_nc,&
write_final_info
USE gopt_f_types, ONLY: gopt_f_type
USE gopt_param_types, ONLY: gopt_param_type
@ -144,8 +145,13 @@ CONTAINS
END IF
CALL write_restart(force_env=force_env,root_section=root_section,&
force_env_section=force_env_section,error=error)
IF (its==gopt_param%max_iter) EXIT
END DO
IF(its == gopt_param%max_iter .AND. (.NOT.converged))THEN
CALL print_geo_opt_nc(gopt_env, output_unit)
END IF
! Write final particle information, if converged
CALL cp_iterate(logger%iter_info,last=.TRUE.,iter_nr=its,error=error)
CALL write_final_info(converged, its, gopt_env, optimizer%x, optimizer%master,&
@ -153,6 +159,8 @@ CONTAINS
force_env_section, error)
CALL cp_opt_gopt_release(optimizer, error=error)
CALL cp_print_key_finished_output(output_unit,logger,geo_section,&
"PRINT%PROGRAM_RUN_INFO", error=error)
END IF
CALL cp_rm_iter_level(logger%iter_info,ihandle,error)
CALL timestop(0.0_dp,handle)

80
src/dimer_methods.F
View file

@ -34,11 +34,8 @@ MODULE dimer_methods
USE geo_opt, ONLY: cp_rot_opt
USE gopt_f_types, ONLY: gopt_f_type
USE input_constants, ONLY: do_first_rotation_step,&
do_init_dimer,&
do_rotate_dimer,&
do_second_rotation_step,&
do_third_rotation_step,&
do_translate_dimer
do_third_rotation_step
USE kinds, ONLY: dp
USE particle_list_types, ONLY: particle_list_type
USE timings, ONLY: timeset,&
@ -102,67 +99,76 @@ CONTAINS
CPPostcondition(ASSOCIATED(dimer_env),cp_failure_level,routineP,error,failure)
! Possibly rotate Dimer or just compute Gradients of point 0 for Translation
SELECT CASE(dimer_env%status_dimer)
CASE(do_rotate_dimer)
SELECT CASE(dimer_env%rotation_step)
IF (gopt_env%dimer_rotation) THEN
SELECT CASE(dimer_env%rot%rotation_step)
CASE(do_first_rotation_step,do_third_rotation_step)
eval_analytical = .TRUE.
IF ((dimer_env%rotation_step==do_third_rotation_step).AND.(dimer_env%interpolate_gradient)) THEN
IF ((dimer_env%rot%rotation_step==do_third_rotation_step).AND.(dimer_env%rot%interpolate_gradient)) THEN
eval_analytical = .FALSE.
END IF
IF (eval_analytical) THEN
! Compute energy, gradients and rotation vector for R1
CALL cp_eval_at_ts_low(gopt_env, x, 1, dimer_env, calc_force, f1, dimer_env%g1, error)
CALL cp_eval_at_ts_low(gopt_env, x, 1, dimer_env, calc_force, f1, dimer_env%rot%g1, error)
ELSE
angle1 = dimer_env%angle1
angle2 = dimer_env%angle2
dimer_env%g1 = SIN(angle1-angle2)/SIN(angle1)*dimer_env%g1+&
SIN(angle2)/SIN(angle1)*dimer_env%g1p+&
(1.0_dp-COS(angle2)-SIN(angle2)*TAN(angle1/2.0_dp))*dimer_env%g0
angle1 = dimer_env%rot%angle1
angle2 = dimer_env%rot%angle2
dimer_env%rot%g1 = SIN(angle1-angle2)/SIN(angle1)*dimer_env%rot%g1+&
SIN(angle2)/SIN(angle1)*dimer_env%rot%g1p+&
(1.0_dp-COS(angle2)-SIN(angle2)*TAN(angle1/2.0_dp))*dimer_env%rot%g0
END IF
! Determine the theta vector (i.e. the search direction for line minimization)
gradient = -2.0_dp*(dimer_env%g1-dimer_env%g0)
gradient = -2.0_dp*(dimer_env%rot%g1-dimer_env%rot%g0)
gradient = gradient - DOT_PRODUCT(gradient, dimer_env%nvec) * dimer_env%nvec
norm = SQRT(DOT_PRODUCT(gradient, gradient))
f = norm
IF (norm/=0.0_dp) gradient = gradient/norm
! Compute curvature and derivative of the curvature w.r.t. the rotational angle
dimer_env%curvature = DOT_PRODUCT(dimer_env%g1-dimer_env%g0,dimer_env%nvec)/dimer_env%dr
dimer_env%dCdp = 2.0_dp*DOT_PRODUCT(dimer_env%g1-dimer_env%g0,gradient)/dimer_env%dr
dimer_env%rot%curvature = DOT_PRODUCT(dimer_env%rot%g1-dimer_env%rot%g0,dimer_env%nvec)/dimer_env%dr
dimer_env%rot%dCdp = 2.0_dp*DOT_PRODUCT(dimer_env%rot%g1-dimer_env%rot%g0,gradient)/dimer_env%dr
dimer_env%rotation_step = do_second_rotation_step
dimer_env%rot%rotation_step = do_second_rotation_step
gradient = -gradient
CASE(do_second_rotation_step)
! Compute energy, gradients and rotation vector for R1
CALL cp_eval_at_ts_low(gopt_env, x, 1, dimer_env, calc_force, f1, dimer_env%g1p, error)
dimer_env%curvature = DOT_PRODUCT(dimer_env%g1p-dimer_env%g0,dimer_env%nvec)/dimer_env%dr
dimer_env%rotation_step = do_third_rotation_step
CALL cp_eval_at_ts_low(gopt_env, x, 1, dimer_env, calc_force, f1, dimer_env%rot%g1p, error)
dimer_env%rot%curvature = DOT_PRODUCT(dimer_env%rot%g1p-dimer_env%rot%g0,dimer_env%nvec)/dimer_env%dr
dimer_env%rot%rotation_step = do_third_rotation_step
! Determine the theta vector (i.e. the search direction for line minimization)
! This is never used for getting a new theta but may is consistent in order to
! give back the right value of f
gradient = -2.0_dp*(dimer_env%rot%g1p-dimer_env%rot%g0)
gradient = gradient - DOT_PRODUCT(gradient, dimer_env%nvec) * dimer_env%nvec
norm = SQRT(DOT_PRODUCT(gradient, gradient))
f = norm
IF (norm/=0.0_dp) gradient = gradient/norm
END SELECT
CASE(do_init_dimer, do_translate_dimer)
ELSE
! Compute energy, gradients and rotation vector for R0
CALL cp_eval_at_ts_low(gopt_env, x, 0, dimer_env, calc_force, f, dimer_env%g0, error)
CALL cp_eval_at_ts_low(gopt_env, x, 0, dimer_env, calc_force, f, dimer_env%rot%g0, error)
IF (dimer_env%status_dimer==do_init_dimer) THEN
dimer_env%status_dimer = do_rotate_dimer
! The dimer is rotated only when we are out of the translation line search
IF (.NOT.gopt_env%do_line_search) THEN
CALL cp_rot_opt(gopt_env%gopt_dimer_env, x, gopt_env%gopt_dimer_param ,&
gopt_env%gopt_dimer_env%geo_section, error)
dimer_env%status_dimer = do_translate_dimer
dimer_env%rot%rotation_step = do_first_rotation_step
END IF
! Correcting gradients for Translation
IF (calc_force) THEN
IF (dimer_env%curvature>0) THEN
gradient = - DOT_PRODUCT(dimer_env%g0,dimer_env%nvec)*dimer_env%nvec
ELSE
gradient = dimer_env%g0 - 2.0_dp*DOT_PRODUCT(dimer_env%g0,dimer_env%nvec)*dimer_env%nvec
END IF
IF (dimer_env%rot%curvature>0) THEN
gradient = - DOT_PRODUCT(dimer_env%rot%g0,dimer_env%nvec)*dimer_env%nvec
ELSE
gradient = dimer_env%rot%g0 - 2.0_dp*DOT_PRODUCT(dimer_env%rot%g0,dimer_env%nvec)*dimer_env%nvec
END IF
! To Fix the pseudo-energy
f = HUGE(0.0_dp)
END SELECT
IF (.NOT.gopt_env%do_line_search) THEN
f = SQRT(DOT_PRODUCT(gradient,gradient))
ELSE
f = -DOT_PRODUCT(gradient,dimer_env%tsl%tls_vec)
END IF
END IF
CALL timestop(0.0_dp,handle)
END SUBROUTINE cp_eval_at_ts
@ -220,7 +226,7 @@ CONTAINS
END DO
!Compute energy and forces
CALL force_env_calc_energy_force(gopt_env%force_env,calc_force=PRESENT(gradient),error=error)
CALL force_env_calc_energy_force(gopt_env%force_env,calc_force=calc_force,error=error)
! Possibly take the potential energy
IF (PRESENT(f)) THEN

70
src/dimer_types.F
View file

@ -25,8 +25,7 @@
MODULE dimer_types
USE f77_blas
USE global_types, ONLY: global_environment_type
USE input_constants, ONLY: do_first_rotation_step,&
do_init_dimer
USE input_constants, ONLY: do_first_rotation_step
USE input_section_types, ONLY: section_vals_type,&
section_vals_val_get
USE kinds, ONLY: dp
@ -42,7 +41,12 @@ MODULE dimer_types
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dimer_types'
INTEGER, PRIVATE, SAVE :: last_dimer_id=0
PUBLIC :: dimer_env_type, dimer_env_create, dimer_env_retain, dimer_env_release
PUBLIC :: dimer_env_type,&
dimer_rotational_type,&
dimer_translational_type ,&
dimer_env_create,&
dimer_env_retain,&
dimer_env_release
!!****s* dimer_types/dimer_env_type [1.0] *
!!
@ -65,14 +69,30 @@ MODULE dimer_types
!! none
!!
!! SOURCE
!***************************************************************************
!***************************************************************************
! Type containing all informations abour the rotation of the Dimer
TYPE dimer_rotational_type
! Rotational parameters
INTEGER :: rotation_step
LOGICAL :: interpolate_gradient
REAL(KIND=dp) :: angle_tol, angle1, angle2, dCdp, curvature
REAL(KIND=dp), POINTER, DIMENSION(:) :: g0, g1, g1p
END TYPE dimer_rotational_type
! Type containing all informations abour the translation of the Dimer
TYPE dimer_translational_type
! Translational parameters
REAL(KIND=dp), POINTER, DIMENSION(:) :: tls_vec
END TYPE dimer_translational_type
! Dimer Environment Type
TYPE dimer_env_type
INTEGER :: ref_count, id_nr
INTEGER :: status_dimer, rotation_step
LOGICAL :: interpolate_gradient
REAL(KIND=dp) :: dr, angle_tol
REAL(KIND=dp), POINTER, DIMENSION(:) :: nvec, g0, g1, g1p
REAL(KIND=dp) :: angle1, angle2, dCdp, curvature
REAL(KIND=dp) :: dr
REAL(KIND=dp), POINTER, DIMENSION(:) :: nvec
TYPE(dimer_rotational_type) :: rot
TYPE(dimer_translational_type) :: tsl
END TYPE dimer_env_type
CONTAINS
@ -121,27 +141,27 @@ CONTAINS
last_dimer_id=last_dimer_id+1
dimer_env%id_nr=last_dimer_id
! Setup NVEC
NULLIFY(dimer_env%nvec)
NULLIFY(dimer_env%nvec, dimer_env%rot%g0, dimer_env%rot%g1, dimer_env%rot%g1p,&
dimer_env%tsl%tls_vec)
! Allocate the working arrays
ALLOCATE(dimer_env%nvec(natom*3),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(dimer_env%g0(natom*3),stat=stat)
ALLOCATE(dimer_env%rot%g0(natom*3),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(dimer_env%g1(natom*3),stat=stat)
ALLOCATE(dimer_env%rot%g1(natom*3),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ALLOCATE(dimer_env%g1p(natom*3),stat=stat)
ALLOCATE(dimer_env%rot%g1p(natom*3),stat=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
CALL random_numbers(dimer_env%nvec, globenv%gaussian_rng_stream, error)
norm = SQRT(SUM(dimer_env%nvec**2))
CPPostcondition(norm>0.0_dp,cp_failure_level,routineP,error,failure)
dimer_env%nvec = dimer_env%nvec / norm
dimer_env%status_dimer = do_init_dimer
dimer_env%rotation_step= do_first_rotation_step
dimer_env%nvec = dimer_env%nvec / norm
dimer_env%rot%rotation_step= do_first_rotation_step
CALL section_vals_val_get(dimer_section,"DR",r_val=dimer_env%dr,error=error)
CALL section_vals_val_get(dimer_section,"INTERPOLATE_GRADIENT",&
l_val=dimer_env%interpolate_gradient,error=error)
l_val=dimer_env%rot%interpolate_gradient,error=error)
CALL section_vals_val_get(dimer_section,"ANGLE_TOLERANCE",&
r_val=dimer_env%angle_tol,error=error)
r_val=dimer_env%rot%angle_tol,error=error)
END IF
END SUBROUTINE dimer_env_create
!***************************************************************************
@ -228,18 +248,20 @@ CONTAINS
DEALLOCATE(dimer_env%nvec, stat=stat)
CPPostcondition(stat==0,cp_warning_level,routineP,error,failure)
END IF
IF (ASSOCIATED(dimer_env%g0)) THEN
DEALLOCATE(dimer_env%g0, stat=stat)
IF (ASSOCIATED(dimer_env%rot%g0)) THEN
DEALLOCATE(dimer_env%rot%g0, stat=stat)
CPPostcondition(stat==0,cp_warning_level,routineP,error,failure)
END IF
IF (ASSOCIATED(dimer_env%g1)) THEN
DEALLOCATE(dimer_env%g1, stat=stat)
IF (ASSOCIATED(dimer_env%rot%g1)) THEN
DEALLOCATE(dimer_env%rot%g1, stat=stat)
CPPostcondition(stat==0,cp_warning_level,routineP,error,failure)
END IF
IF (ASSOCIATED(dimer_env%g1p)) THEN
DEALLOCATE(dimer_env%g1p, stat=stat)
IF (ASSOCIATED(dimer_env%rot%g1p)) THEN
DEALLOCATE(dimer_env%rot%g1p, stat=stat)
CPPostcondition(stat==0,cp_warning_level,routineP,error,failure)
END IF
! No need to deallocate tls_vec (just a pointer to aother local array)
NULLIFY(dimer_env%tsl%tls_vec)
DEALLOCATE(dimer_env, stat=stat)
CPPostcondition(stat==0,cp_warning_level,routineP,error,failure)
END IF

6
src/gopt_f77_methods.F
View file

@ -117,8 +117,10 @@
ALLOCATE(gradient_ts(particles%n_els*3),stat=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
! Real calculation of energy and forces for transition state optimization
CALL cp_eval_at_ts (gopt_env, x, f_ts, gradient_ts, calc_force=PRESENT(gradient), error=error)
! Real calculation of energy and forces for transition state optimization:
! When doing dimer methods forces have to be always computed since the function
! to minimize is not the energy but the effective force
CALL cp_eval_at_ts (gopt_env, x, f_ts, gradient_ts, calc_force=.TRUE., error=error)
! Possibly take the potential energy
IF (PRESENT(f)) f = f_ts

46
src/gopt_f_methods.F
View file

@ -66,7 +66,8 @@ MODULE gopt_f_methods
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gopt_f_methods'
PUBLIC :: gopt_f_create_x0, write_cycle_infos, check_converg, write_final_info,&
write_geo_traj, print_geo_opt_header, check_rot_conv, write_rot_cycle_infos
write_geo_traj, print_geo_opt_header, check_rot_conv, write_rot_cycle_infos,&
print_geo_opt_nc
CONTAINS
@ -233,7 +234,7 @@ CONTAINS
WRITE(UNIT=output_unit,FMT="(T2,A,T47,A)")&
" Optimization Method = ",wildcard
WRITE(UNIT=output_unit,FMT="(T2,A,F20.10)")&
" Local Curvature = ",dimer_env%curvature
" Local Curvature = ",dimer_env%rot%curvature
WRITE(UNIT=output_unit,FMT="(T2,A,F20.10)")&
" Total Rotational Force = ",etot
IF (PRESENT(ediff)) THEN
@ -413,16 +414,16 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'check_rot_conv', &
routineP = moduleN//':'//routineN
conv = (dimer_env%angle1< dimer_env%angle_tol)
conv = (dimer_env%rot%angle1< dimer_env%rot%angle_tol)
IF(output_unit>0)THEN
WRITE(UNIT=output_unit,FMT="(/,T2,A)")&
" Convergence check :"
WRITE(UNIT=output_unit,FMT="(T2,A,F16.10)")&
" Predicted angle step size = ",dimer_env%angle1
" Predicted angle step size = ",dimer_env%rot%angle1
WRITE(UNIT=output_unit,FMT="(T2,A,F16.10)")&
" Effective angle step size = ",dimer_env%angle2
" Effective angle step size = ",dimer_env%rot%angle2
WRITE(UNIT=output_unit,FMT="(T2,A,F16.10)")&
" Conv. limit for angle step size =",dimer_env%angle_tol
" Conv. limit for angle step size =",dimer_env%rot%angle_tol
IF(conv)THEN
WRITE(UNIT=output_unit,FMT="(T2,2A)")&
" Convergence in angle step size =",&
@ -637,4 +638,37 @@ CONTAINS
END IF
END SUBROUTINE print_geo_opt_header
!!****** gopt_f_methods/print_geo_opt_nc [1.0] *
!!
!! NAME
!! print_geo_opt_nc
!!
!! FUNCTION
!!
!!
!! AUTHOR
!! Teodoro Laino [tlaino] - University of Zurich - 01.2008
!!
!! MODIFICATION HISTORY
!!
!! SOURCE
!***************************************************************************
SUBROUTINE print_geo_opt_nc(gopt_env,output_unit)
TYPE(gopt_f_type), POINTER :: gopt_env
INTEGER, INTENT(IN) :: output_unit
IF(output_unit>0)THEN
WRITE(UNIT=output_unit,FMT="(/,T2,A)")&
"*** MAXIMUM NUMBER OF OPTIMIZATION STEPS REACHED ***"
IF (.NOT.gopt_env%dimer_rotation) THEN
WRITE(UNIT=output_unit,FMT="(T2,A)")&
"*** EXITING GEOMETRY OPTIMIZATION ***"
ELSE
WRITE(UNIT=output_unit,FMT="(T2,A)")&
"*** EXITING ROTATION OPTIMIZATION ***"
END IF
CALL m_flush(output_unit)
END IF
END SUBROUTINE print_geo_opt_nc
END MODULE gopt_f_methods

4
src/gopt_f_types.F
View file

@ -82,7 +82,7 @@ MODULE gopt_f_types
TYPE gopt_f_type
INTEGER :: ref_count, id_nr, its
INTEGER :: type_id, ts_method_id
LOGICAL :: dimer_rotation
LOGICAL :: dimer_rotation, do_line_search
CHARACTER(LEN=default_string_length) :: label
TYPE(force_env_type), POINTER :: force_env
TYPE(global_environment_type), POINTER :: globenv
@ -148,6 +148,7 @@ CONTAINS
gopt_env%its=0
gopt_env%label="GEO_OPT"
gopt_env%dimer_rotation=.FALSE.
gopt_env%do_line_search=.FALSE.
CALL force_env_retain(force_env, error=error)
gopt_env%force_env => force_env
gopt_env%globenv => globenv
@ -176,7 +177,6 @@ CONTAINS
gopt_env%gopt_dimer_env%dimer_env => gopt_env%dimer_env
gopt_env%gopt_dimer_env%label = "ROT_OPT"
gopt_env%gopt_dimer_env%dimer_rotation = .TRUE.
END SELECT
END IF
END IF

38
src/gopt_param_types.F
View file

@ -45,7 +45,7 @@ MODULE gopt_param_types
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gopt_param_types'
INTEGER, PRIVATE, SAVE :: last_gopt_f_id=0, last_gopt_param_id=0
PUBLIC :: gopt_param_type
PUBLIC :: gopt_param_type, cg_ls_param_type
PUBLIC :: gopt_param_create, gopt_param_read, gopt_param_retain, &
gopt_param_release
@ -76,17 +76,23 @@ MODULE gopt_param_types
!!
!! SOURCE
!***************************************************************************
TYPE cg_ls_param_type
INTEGER :: brent_max_iter, type_id
REAL(KIND = dp) :: brent_tol, brack_limit, initial_step
END TYPE cg_ls_param_type
TYPE gopt_param_type
LOGICAL :: Fletcher_Reeves
LOGICAL :: use_fit_line_search
INTEGER :: id_nr, ref_count, method_id, type_id
INTEGER :: ts_method_id
INTEGER :: max_f_per_iter, max_iter, max_h_rank
INTEGER :: brent_max_iter, max_steep_steps
REAL(KIND = dp) :: wanted_proj_gradient, wanted_rel_f_error
REAL(KIND = dp) :: max_dr, max_force, rms_dr, rms_force
REAL(KIND = dp) :: brent_tol, brack_limit, restart_limit, initial_step
REAL(KIND = dp) :: dimer_angle_tol
LOGICAL :: Fletcher_Reeves
INTEGER :: id_nr, ref_count
INTEGER :: method_id, type_id
INTEGER :: ts_method_id
INTEGER :: max_f_per_iter, max_iter, max_h_rank
INTEGER :: max_steep_steps
REAL(KIND = dp) :: restart_limit
REAL(KIND = dp) :: wanted_proj_gradient, wanted_rel_f_error
REAL(KIND = dp) :: max_dr, max_force, rms_dr, rms_force
REAL(KIND = dp) :: dimer_angle_tol
TYPE(cg_ls_param_type) :: cg_ls
END TYPE gopt_param_type
!!***
!****************************************************************************
@ -204,13 +210,13 @@ CONTAINS
! Do nothing
CASE(default_cg_method_id)
CALL section_vals_val_get(gopt_section,"CG%FLETCHER_REEVES",l_val=gopt_param%Fletcher_Reeves,error=error)
CALL section_vals_val_get(gopt_section,"CG%USE_FIT_LINE_SEARCH",l_val=gopt_param%use_fit_line_search,error=error)
CALL section_vals_val_get(gopt_section,"CG%MAX_STEEP_STEPS",i_val=gopt_param%max_steep_steps,error=error)
CALL section_vals_val_get(gopt_section,"CG%BRENT_TOL",r_val=gopt_param%brent_tol,error=error)
CALL section_vals_val_get(gopt_section,"CG%INITIAL_STEP",r_val=gopt_param%initial_step,error=error)
CALL section_vals_val_get(gopt_section,"CG%BRENT_MAX_ITER",i_val=gopt_param%brent_max_iter,error=error)
CALL section_vals_val_get(gopt_section,"CG%BRACK_LIMIT",r_val=gopt_param%brack_limit,error=error)
CALL section_vals_val_get(gopt_section,"CG%RESTART_LIMIT",r_val=gopt_param%restart_limit,error=error)
CALL section_vals_val_get(gopt_section,"CG%LINE_SEARCH%TYPE",i_val=gopt_param%cg_ls%type_id,error=error)
CALL section_vals_val_get(gopt_section,"CG%LINE_SEARCH%BRENT_TOL",r_val=gopt_param%cg_ls%brent_tol,error=error)
CALL section_vals_val_get(gopt_section,"CG%LINE_SEARCH%INITIAL_STEP",r_val=gopt_param%cg_ls%initial_step,error=error)
CALL section_vals_val_get(gopt_section,"CG%LINE_SEARCH%BRENT_MAX_ITER",i_val=gopt_param%cg_ls%brent_max_iter,error=error)
CALL section_vals_val_get(gopt_section,"CG%LINE_SEARCH%BRACK_LIMIT",r_val=gopt_param%cg_ls%brack_limit,error=error)
END SELECT
SELECT CASE(gopt_param%type_id)

8
src/input_constants.F
View file

@ -408,7 +408,7 @@ MODULE input_constants
INTEGER, PARAMETER, PUBLIC :: ot_mini_sd=1,ot_mini_cg=2,ot_mini_diis=3
INTEGER, PARAMETER, PUBLIC :: ot_algo_taylor_or_diag=1,ot_algo_irac=2
INTEGER, PARAMETER, PUBLIC :: ot_chol_irac=1,ot_poly_irac=2,ot_lwdn_irac=3
INTEGER, PARAMETER, PUBLIC :: ot_ls_none=1,ot_ls_2pnt=2,ot_ls_3pnt=3,ot_ls_gold=4
INTEGER, PARAMETER, PUBLIC :: ls_none=1,ls_2pnt=2,ls_3pnt=3,ls_gold=4,ls_fit=5
INTEGER, PARAMETER, PUBLIC :: ot_precond_none=0, &
ot_precond_full_single=1, &
ot_precond_full_kinetic=2, &
@ -430,10 +430,6 @@ MODULE input_constants
INTEGER, PUBLIC, PARAMETER :: none_ts_method_id=0,&
default_dimer_method_id=1
INTEGER, PUBLIC, PARAMETER :: do_init_dimer=0,&
do_rotate_dimer=1,&
do_translate_dimer=2
INTEGER, PUBLIC, PARAMETER :: do_first_rotation_step=0,&
do_second_rotation_step=1,&
do_third_rotation_step=2
@ -549,8 +545,6 @@ MODULE input_constants
lr_precond_s_inverse=3,&
lr_precond_kinetic=4, &
lr_precond_full_all=6
INTEGER, PARAMETER, PUBLIC :: lr_ls_none=0, lr_ls_2pnt=1,&
lr_ls_3pnt=2, lr_ls_gold=3
INTEGER, PARAMETER, PUBLIC :: xc_pbe_orig=11,xc_rev_pbe=12

8
src/input_cp2k_dft.F
View file

@ -4552,11 +4552,11 @@ CONTAINS
" Whereas the 2PNT minimizer is almost always OK, 3PNT might be needed for systems"//&
" in which successive OT CG steps do not decrease the total energy.",&
usage="LINESEARCH GOLD",&
default_i_val=ot_ls_2pnt,&
default_i_val=ls_2pnt,&
enum_c_vals=s2a( "NONE", "2PNT", "3PNT","GOLD"),&
enum_desc=s2a("take fixed lenght steps","extrapolate based on 2 points", &
"... or on 3 points","perform 1D golden section search of the minimum (very expensive)"),&
enum_i_vals=(/ot_ls_none,ot_ls_2pnt,ot_ls_3pnt,ot_ls_gold/),&
enum_i_vals=(/ls_none,ls_2pnt,ls_3pnt,ls_gold/),&
error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
@ -5026,12 +5026,12 @@ CONTAINS
description="1D line search algorithm,"//&
" in increasing order of robustness and cost.",&
usage="LINESEARCH GOLD",&
default_i_val=ot_ls_none,&
default_i_val=ls_none,&
enum_c_vals=s2a( "NONE", "2PNT", "3PNT","GOLD"),&
enum_desc=s2a("take fixed lenght steps","extrapolate based on 2 points", &
"... or on 3 points","perform 1D golden section search of the minimum "//&
"(WARNING this option doesnt work yet)"),&
enum_i_vals=(/lr_ls_none,lr_ls_2pnt,lr_ls_3pnt,lr_ls_gold/),&
enum_i_vals=(/ls_none,ls_2pnt,ls_3pnt,ls_gold/),&
error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)

102
src/input_cp2k_motion.F
View file

@ -347,7 +347,7 @@ CONTAINS
CALL section_release(print_key,error=error)
CALL cp_print_key_section_create(print_key,"TRANSLATION_VECTOR",&
description="Dumps the traslation vector applied along an MD (if any). Useful"//&
description="Dumps the translation vector applied along an MD (if any). Useful"//&
" for postprocessing of QMMM trajectories in which the QM fragment is continuously"//&
" centered in the QM box",&
print_level=high_print_level, common_iter_levels=1,&
@ -942,8 +942,8 @@ CONTAINS
CALL keyword_create(keyword, name="INTERPOLATE_GRADIENT",&
description="This keyword controls the interpolation of the gradient whenever possible"//&
" during the optimization of the Dimer. The use of this keywords saves 1 evaluation "//&
" of energy/forces.",&
usage="INTERPOLATE_GRADIENT {logical}",default_l_val=.TRUE.,error=error)
" of energy/forces.", usage="INTERPOLATE_GRADIENT {logical}",default_l_val=.TRUE.,&
lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(subsection,keyword,error=error)
CALL keyword_release(keyword,error=error)
@ -1061,47 +1061,18 @@ CONTAINS
LOGICAL :: failure
TYPE(keyword_type), POINTER :: keyword
TYPE(section_type), POINTER :: subsection
failure=.FALSE.
IF (.NOT.failure) THEN
! create the CG subsection
NULLIFY(section,keyword)
NULLIFY(section,subsection,keyword)
CALL section_create(section,name="CG",&
description="Provides parameters to tune the conjugate gradient optimization",&
n_keywords=0, n_subsections=1, repeats=.FALSE., required=.FALSE.,&
error=error)
CALL keyword_create(keyword, name="BRENT_TOL",&
description="Tolerance requested during Brent line search in Conjugate Gradients Optimization.",&
usage="BRENT_TOL {real}",&
default_r_val=0.01_dp,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="INITIAL_STEP",&
description="Initial step size used, e.g. for bracketing or minimizers. "//&
"Might need to be reduced for systems with close contacts",&
usage="INITIAL_STEP {real}",&
default_r_val=0.2_dp,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="BRACK_LIMIT",&
description="Limit in 1D bracketing during line search in Conjugate Gradients Optimization.",&
usage="BRACK_LIMIT {real}",&
default_r_val=100.0_dp,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="BRENT_MAX_ITER",&
description="Maximum number of iterations in brent algorithm "// &
"(used for the line search in Conjugated Gradients Optimization)",&
usage="BRENT_MAX_ITER {integer}",&
default_i_val=100,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="MAX_STEEP_STEPS",&
description="Maximum number of steepest descent steps before starting the"//&
" conjugate gradients optimization.",&
@ -1120,14 +1091,6 @@ CONTAINS
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="USE_FIT_LINE_SEARCH",&
description="Use instead of bracketing the minima and than applying the Brent"//&
" algorithm, a least mean square to interpolate the minima in the line search.",&
usage="USE_FIT_LINE_SEARCH",&
default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="FLETCHER_REEVES",&
description="Uses FLETCHER-REEVES instead of POLAK-RIBIERE when using Conjugate Gradients",&
usage="FLETCHER-REEVES",&
@ -1135,6 +1098,61 @@ CONTAINS
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
! Line Search section
CALL section_create(subsection,name="LINE_SEARCH",&
description="Provides parameters to tune the line search during the conjugate gradient optimization",&
n_keywords=0, n_subsections=1, repeats=.FALSE., required=.FALSE.,&
error=error)
CALL keyword_create(keyword, name="TYPE",&
description="1D line search algorithm to be used with the CG optimizer,"//&
" in increasing order of robustness and cost. ",&
usage="TYPE GOLD",&
default_i_val=ls_gold,&
enum_c_vals=s2a( "NONE","2PNT","3PNT","GOLD","FIT"),&
enum_desc=s2a("take fixed lenght steps",&
"extrapolate based on 2 points", &
"extrapolate based on on 3 points",&
"perform 1D golden section search of the minimum (very expensive)",&
"perform 1D fit of a parabola on several evaluation of energy "//&
"(very expensive and more robust vs numerical noise)"),&
enum_i_vals=(/ls_none,ls_2pnt,ls_3pnt,ls_gold,ls_fit/),&
error=error)
CALL section_add_keyword(subsection,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="INITIAL_STEP",&
description="Initial step size used, e.g. for bracketing or minimizers. "//&
"Might need to be reduced for systems with close contacts",&
usage="INITIAL_STEP {real}",&
default_r_val=0.2_dp,error=error)
CALL section_add_keyword(subsection,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="BRACK_LIMIT",&
description="Limit in 1D bracketing during line search in Conjugate Gradients Optimization.",&
usage="BRACK_LIMIT {real}",&
default_r_val=100.0_dp,error=error)
CALL section_add_keyword(subsection,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="BRENT_TOL",&
description="Tolerance requested during Brent line search in Conjugate Gradients Optimization.",&
usage="BRENT_TOL {real}",&
default_r_val=0.01_dp,error=error)
CALL section_add_keyword(subsection,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="BRENT_MAX_ITER",&
description="Maximum number of iterations in brent algorithm "// &
"(used for the line search in Conjugated Gradients Optimization)",&
usage="BRENT_MAX_ITER {integer}",&
default_i_val=100,error=error)
CALL section_add_keyword(subsection,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL section_add_subsection(section,subsection,error=error)
CALL section_release(subsection,error=error)
END IF
END SUBROUTINE create_cg_section

12
src/qs_linres_methods.F
View file

@ -46,8 +46,8 @@ MODULE qs_linres_methods
USE cp_sm_fm_interactions, ONLY: cp_sm_fm_multiply
USE global_types, ONLY: global_environment_type
USE input_constants, ONLY: &
do_loc_none, lr_cg, lr_diis, lr_ls_2pnt, lr_ls_3pnt, lr_ls_gold, &
lr_ls_none, lr_precond_none, lr_sd, op_loc_berry, state_loc_all
do_loc_none, lr_cg, lr_diis, lr_precond_none, lr_sd, ls_2pnt, ls_3pnt, &
ls_gold, ls_none, op_loc_berry, state_loc_all
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
@ -736,13 +736,13 @@ CONTAINS
!
SELECT CASE(linres_control%ls_method)
CASE(lr_ls_gold)
CASE(ls_gold)
CALL stop_program("do_lr_line_search","GOLD NYI")
CASE(lr_ls_3pnt)
CASE(ls_3pnt)
CALL do_lr_ls_3pnt(p_env,vectors,linres_control%ds_min,nspins,error=error)
CASE(lr_ls_2pnt)
CASE(ls_2pnt)
CALL do_lr_ls_2pnt(p_env,vectors,linres_control%ds_min,nspins,error=error)
CASE(lr_ls_none)
CASE(ls_none)
ds=linres_control%ds_min
DO ispin = 1,nspins
CALL cp_fm_scale_and_add(1.0_dp,vectors(ispin)%matrix,ds,&

14
src/qs_linres_module.F
View file

@ -32,9 +32,9 @@ MODULE qs_linres_module
cp_print_key_unit_nr
USE global_types, ONLY: global_environment_type
USE input_constants, ONLY: &
lr_cg, lr_chemshift, lr_diis, lr_gtensor, lr_ls_2pnt, lr_ls_3pnt, &
lr_ls_gold, lr_ls_none, lr_none, lr_precond_full_all, lr_precond_none, &
lr_precond_s_inverse, lr_precond_single, lr_sd
lr_cg, lr_chemshift, lr_diis, lr_gtensor, lr_none, &
lr_precond_full_all, lr_precond_none, lr_precond_s_inverse, &
lr_precond_single, lr_sd, ls_2pnt, ls_3pnt, ls_gold, ls_none
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
@ -363,16 +363,16 @@ CONTAINS
"LINRES| Optimization algorithm"," Steepnd Descent"
IF(linres_control%opt_method==lr_diis) WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")&
"LINRES| Optimization algorithm"," DIIS"
IF(linres_control%ls_method==lr_ls_none) &
IF(linres_control%ls_method==ls_none) &
WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")&
"LINRES| Line Search"," NONE"
IF(linres_control%ls_method==lr_ls_2pnt) &
IF(linres_control%ls_method==ls_2pnt) &
WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")&
"LINRES| Line Search"," 2pnt"
IF(linres_control%ls_method==lr_ls_3pnt) &
IF(linres_control%ls_method==ls_3pnt) &
WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")&
"LINRES| Line Search"," 3pnt"
IF(linres_control%ls_method==lr_ls_gold) &
IF(linres_control%ls_method==ls_gold) &
WRITE (UNIT=output_unit,FMT="(T2,A,T60,A)")&
"LINRES| Line Search"," gold"
IF(linres_control%preconditioner_type==lr_precond_none) &

12
src/qs_ot_scf.F
View file

@ -35,8 +35,8 @@ MODULE qs_ot_scf
USE cp_sm_fm_interactions, ONLY: cp_sm_fm_multiply
USE f77_blas
USE input_constants, ONLY: &
ot_algo_irac, ot_algo_taylor_or_diag, ot_chol_irac, ot_ls_2pnt, &
ot_ls_3pnt, ot_ls_gold, ot_ls_none, ot_lwdn_irac, ot_mini_cg, &
ls_2pnt, ls_3pnt, ls_gold, ls_none, ot_algo_irac, &
ot_algo_taylor_or_diag, ot_chol_irac, ot_lwdn_irac, ot_mini_cg, &
ot_mini_diis, ot_mini_sd, ot_poly_irac, ot_precond_full_all, &
ot_precond_full_kinetic, ot_precond_full_single, &
ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, &
@ -181,13 +181,13 @@ CONTAINS
CALL section_vals_val_get(ot_section,"SAFER_DIIS",l_val=qs_ot_env(1)%settings%safer_diis,error=error)
CALL section_vals_val_get(ot_section,"LINESEARCH",i_val=ls_method,error=error)
SELECT CASE(ls_method)
CASE (ot_ls_none)
CASE (ls_none)
qs_ot_env(1)%settings%line_search_method="NONE"
CASE (ot_ls_2pnt)
CASE (ls_2pnt)
qs_ot_env(1)%settings%line_search_method="2PNT"
CASE (ot_ls_3pnt)
CASE (ls_3pnt)
qs_ot_env(1)%settings%line_search_method="3PNT"
CASE (ot_ls_gold)
CASE (ls_gold)
qs_ot_env(1)%settings%line_search_method="GOLD"
CALL section_vals_val_get(ot_section,"GOLD_TARGET",r_val=qs_ot_env(1)%settings%gold_target,error=error)
CASE DEFAULT

6
src/qs_p_env_methods.F
View file

@ -60,8 +60,8 @@ MODULE qs_p_env_methods
USE hartree_local_methods, ONLY: init_coulomb_local
USE hartree_local_types, ONLY: hartree_local_create
USE input_constants, ONLY: lr_cg,&
lr_ls_none,&
lr_precond_none
lr_precond_none,&
ls_none
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type
@ -374,7 +374,7 @@ CONTAINS
! error=error)
END IF
END IF
IF(linres_control%ls_method /= lr_ls_none) THEN
IF(linres_control%ls_method /= ls_none) THEN
! Initialize the required variables for the line search
END IF

3
tests/Fist/aala.inp
View file

@ -34,8 +34,5 @@
&GEO_OPT
OPTIMIZER CG
MAX_ITER 500
&BFGS
RESTART
&END
&END
&END

2
tests/Fist/regtest-3/TEST_FILES
View file

@ -17,3 +17,5 @@ cu_eam_2.inp 2
# Fumi-Tosi Potential
NaCl.inp 2
NaCl-H2O.inp 2
#Dimer Method
water_2_TS.inp 25

77
tests/Fist/regtest-3/water_2_TS.inp Normal file
View file

@ -0,0 +1,77 @@
&FORCE_EVAL
METHOD FIST
&MM
&FORCEFIELD
parm_file_name ../sample_pot/water.pot
parmtype CHM
&CHARGE
ATOM OT
CHARGE -0.8476
&END CHARGE
&CHARGE
ATOM HT
CHARGE 0.4238
&END CHARGE
&END FORCEFIELD
&POISSON
&EWALD
EWALD_TYPE spme
ALPHA .44
GMAX 24
O_SPLINE 6
&END EWALD
&END POISSON
&END MM
&SUBSYS
&CELL
ABC 15.0 15.0 15.0
UNIT ANGSTROM
&END CELL
&TOPOLOGY
COORD_FILE_NAME ../sample_xyz/water_2_TS.xyz
COORDINATE xyz
&DUMP_PSF
&END
&GENERATE
CREATE_MOLECULES
&END
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT water_2_TS
RUN_TYPE GEO_OPT
#TRACE
&END GLOBAL
&MOTION
&GEO_OPT
TYPE TRANSITION_STATE
MAX_ITER 20
OPTIMIZER CG
&CG
MAX_STEEP_STEPS 20
&LINE_SEARCH
TYPE 2PNT
&END
&END
&TRANSITION_STATE
METHOD DIMER
&DIMER
DR 0.0001
ANGLE_TOLERANCE [deg] 5
INTERPOLATE_GRADIENT
&ROT_OPT
OPTIMIZER CG
MAX_ITER 1000
&CG
MAX_STEEP_STEPS 1000
&LINE_SEARCH
TYPE 2PNT
&END
&END
&END
&END
&END
&END GEO_OPT
&END MOTION

3
tests/Fist/regtest-5/mfi.inp
View file

@ -115,9 +115,6 @@
RUN_TYPE MD
&END GLOBAL
&MOTION
&GEO_OPT
MINIMIZER CG
&END GEO_OPT
&MD
STEPS 20
TEMPERATURE 300

6
tests/Fist/regtest-6/cubane_0.inp
View file

@ -43,17 +43,11 @@
PRINT_LEVEL DEBUG
PROJECT CUBANE
RUN_TYPE GEO_OPT
#RUN_TYPE MD
&END GLOBAL
&MOTION
&GEO_OPT
MAX_ITER 50
&END
&MD
ENSEMBLE NPT_F
STEPS 5
TIMESTEP 0.5
&END
&PRINT
&TRAJECTORY
FILENAME __STD_OUT__

8
tests/Fist/sample_xyz/water_2_TS.xyz Normal file
View file

@ -0,0 +1,8 @@
6
saddle poit as generated by Gamess 6-21G HF
OT 1.2832581715 0.0333112343 0.1267660428
HT 0.7166821936 0.5090935353 -0.4965092418
HT 0.6521873036 -0.4652492545 0.6643731224
OT -1.5084053639 -0.0184355291 -0.0721384111
HT -2.1191020938 0.5718208217 0.3900714787
HT -2.0277721948 -0.6293358052 -0.6125629915

4
tests/QS/regtest-gpw-4/H2-geo-4.inp
View file

@ -46,7 +46,9 @@
MAX_ITER 2
&CG
MAX_STEEP_STEPS 1
USE_FIT_LINE_SEARCH
&LINE_SEARCH
TYPE FIT
&END
&END CG
&END GEO_OPT
&END MOTION

3
tests/TEST_TYPES
View file

@ -1,4 +1,4 @@
24
25
Total energy:!3
POTENTIAL ENERGY!4
Total energy \[eV\]:!4
@ -23,6 +23,7 @@ VIB| !3
PW exchange energy!5
Total Spread!5
DFT+U energy:!3
Local Curvature =!4
#
# these are the tests the can be selected for regtesting.
# do regtest will grep for test_grep (first column) and look if the numeric value