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

more options for lri basis set fitting

Browse source 
svn-origin-rev: 14864
This commit is contained in:
Dorothea Golze 2014-12-19 09:30:53 +00:00
parent 2a35ab6171
commit 8c1ea80eca
9 changed files with 1061 additions and 229 deletions
Download patch file Download diff file Expand all files Collapse all files

5
src/input_constants.F
View file

@ -827,6 +827,11 @@ MODULE input_constants
do_opt_coeff=2,&
do_opt_exps=3
! LRI basis optimization parameters
INTEGER, PARAMETER, PUBLIC :: do_lri_opt_all=0,&
do_lri_opt_coeff=1,&
do_lri_opt_exps=2
! callgraph parameters
INTEGER, PARAMETER, PUBLIC :: callgraph_none=0,&
callgraph_master=1,&

140
src/input_cp2k_dft.F
View file

@ -47,30 +47,30 @@ MODULE input_cp2k_dft
do_full_density, do_gapw_gcs, do_gapw_gct, do_gapw_log, do_loc_both, &
do_loc_crazy, do_loc_direct, do_loc_homo, do_loc_jacobi, &
do_loc_l1_norm_sd, do_loc_lumo, do_loc_max, do_loc_min, do_loc_none, &
do_method_am1, do_method_dftb, do_method_gapw, do_method_gapw_xc, &
do_method_gpw, do_method_lrigpw, do_method_mndo, do_method_mndod, &
do_method_ofgpw, do_method_pdg, do_method_pm3, do_method_pm6, &
do_method_pnnl, do_method_rm1, do_method_scptb, do_pade, &
do_ppl_analytic, do_ppl_grid, do_pwgrid_ns_fullspace, &
do_pwgrid_ns_halfspace, do_pwgrid_spherical, do_s2_constraint, &
do_s2_restraint, do_se_is_kdso, do_se_is_kdso_d, do_se_is_slater, &
do_se_lr_ewald, do_se_lr_ewald_gks, do_se_lr_ewald_r3, do_se_lr_none, &
do_spin_density, do_taylor, ehrenfest, gaussian, gaussian_env, &
general_roks, high_spin_roks, history_guess, kerker_mix, &
kg_color_dsatur, kg_color_greedy, kg_tnadd_atomic, kg_tnadd_embed, &
ls_2pnt, ls_3pnt, ls_gold, ls_none, mopac_guess, multisec_mix, &
no_excitations, no_guess, no_mix, numerical, oe_gllb, oe_lb, oe_none, &
oe_saop, oe_sic, op_loc_berry, op_loc_boys, op_loc_pipek, orb_dx2, &
orb_dxy, orb_dy2, orb_dyz, orb_dz2, orb_dzx, orb_px, orb_py, orb_pz, &
orb_s, ot_algo_irac, ot_algo_taylor_or_diag, ot_chol_irac, &
ot_lwdn_irac, ot_mini_broyden, 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, ot_precond_solver_default, &
ot_precond_solver_direct, ot_precond_solver_inv_chol, &
ot_precond_solver_update, outer_scf_becke_constraint, &
outer_scf_ddapc_constraint, outer_scf_none, &
outer_scf_optimizer_bisect, outer_scf_optimizer_diis, &
do_lri_opt_all, do_lri_opt_coeff, do_lri_opt_exps, do_method_am1, &
do_method_dftb, do_method_gapw, do_method_gapw_xc, do_method_gpw, &
do_method_lrigpw, do_method_mndo, do_method_mndod, do_method_ofgpw, &
do_method_pdg, do_method_pm3, do_method_pm6, do_method_pnnl, &
do_method_rm1, do_method_scptb, do_pade, do_ppl_analytic, do_ppl_grid, &
do_pwgrid_ns_fullspace, do_pwgrid_ns_halfspace, do_pwgrid_spherical, &
do_s2_constraint, do_s2_restraint, do_se_is_kdso, do_se_is_kdso_d, &
do_se_is_slater, do_se_lr_ewald, do_se_lr_ewald_gks, &
do_se_lr_ewald_r3, do_se_lr_none, do_spin_density, do_taylor, &
ehrenfest, gaussian, gaussian_env, general_roks, high_spin_roks, &
history_guess, kerker_mix, kg_color_dsatur, kg_color_greedy, &
kg_tnadd_atomic, kg_tnadd_embed, ls_2pnt, ls_3pnt, ls_gold, ls_none, &
mopac_guess, multisec_mix, no_excitations, no_guess, no_mix, &
numerical, oe_gllb, oe_lb, oe_none, oe_saop, oe_sic, op_loc_berry, &
op_loc_boys, op_loc_pipek, orb_dx2, orb_dxy, orb_dy2, orb_dyz, &
orb_dz2, orb_dzx, orb_px, orb_py, orb_pz, orb_s, ot_algo_irac, &
ot_algo_taylor_or_diag, ot_chol_irac, ot_lwdn_irac, ot_mini_broyden, &
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, &
ot_precond_solver_default, ot_precond_solver_direct, &
ot_precond_solver_inv_chol, ot_precond_solver_update, &
outer_scf_becke_constraint, outer_scf_ddapc_constraint, &
outer_scf_none, outer_scf_optimizer_bisect, outer_scf_optimizer_diis, &
outer_scf_optimizer_none, outer_scf_optimizer_sd, &
outer_scf_s2_constraint, outer_scf_scp, plus_u_lowdin, &
plus_u_mulliken, plus_u_mulliken_charges, pulay_mix, pw_interp, &
@ -5415,6 +5415,7 @@ CONTAINS
LOGICAL :: failure
TYPE(keyword_type), POINTER :: keyword
TYPE(section_type), POINTER :: subsection
failure=.FALSE.
@ -5426,7 +5427,7 @@ CONTAINS
n_keywords=1, n_subsections=0, repeats=.FALSE., required=.FALSE.,&
error=error)
NULLIFY(keyword)
NULLIFY(keyword,subsection)
CALL keyword_create(keyword, name="ACCURACY", &
description="Target accuracy for the objective function (RHOEND)",&
@ -5449,9 +5450,98 @@ CONTAINS
CALL section_add_keyword(section, keyword, error=error)
CALL keyword_release(keyword, error=error)
CALL keyword_create(keyword, name="CONDITION_WEIGHT", &
description="This keyword allows to give different weight "//&
"factors to the condition number (LOG(cond) is used).",&
usage="CONDITION_WEIGHT 1.0E-4", default_r_val=1.0E-6_dp,&
error=error)
CALL section_add_keyword(section, keyword, error=error)
CALL keyword_release(keyword, error=error)
CALL keyword_create(keyword, name="USE_CONDITION_NUMBER",&
description="Determines whether condition number should be part "//&
"of optimization or not",&
usage="USE_CONDITION_NUMBER",&
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="GEOMETRIC_SEQUENCE",&
description="Exponents are assumed to be a geometric squence. "//&
"Only the minimal and maximal exponents of one set are optimized and "//&
"the other exponents are obtained by geometric progression.",&
usage="GEOMETRIC_SEQUENCE",&
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="DEGREES_OF_FREEDOM",&
description="Specifies the degrees of freedom in the basis "//&
"optimization.",&
usage="DEGREES_OF_FREEDOM ALL", &
enum_c_vals=s2a( "ALL", "COEFFICIENTS","EXPONENTS"),&
enum_desc=s2a("Set all parameters in the basis to be variable.",&
"Set all coefficients in the basis set to be variable.",&
"Set all exponents in the basis to be variable."),&
enum_i_vals=(/do_lri_opt_all, do_lri_opt_coeff, do_lri_opt_exps/),&
default_i_val=do_lri_opt_exps,error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL create_constrain_exponents_section(subsection,error)
CALL section_add_subsection(section, subsection, error=error)
CALL section_release(subsection,error=error)
END IF
END SUBROUTINE create_optimize_lri_basis_section
! *****************************************************************************
!> \brief input section for constraints for auxiliary basis set optimization
!> \param section the section to create
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
!> \author Dorothea Golze [11.2014]
! *****************************************************************************
SUBROUTINE create_constrain_exponents_section(section,error)
TYPE(section_type), POINTER :: section
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: &
routineN = 'create_constrain_exponents_section', &
routineP = moduleN//':'//routineN
LOGICAL :: failure
TYPE(keyword_type), POINTER :: keyword
failure=.FALSE.
IF (.NOT. failure) THEN
CALL section_create(section,name="CONSTRAIN_EXPONENTS",&
description="specifies constraints for the exponents of the "//&
"lri auxiliary basis sets in the optimization.",&
n_keywords=1, n_subsections=0, repeats=.FALSE., required=.FALSE.,&
error=error)
NULLIFY(keyword)
CALL keyword_create(keyword, name="SCALE",&
description="Defines the upper and lower boundaries as "//&
"(1+scale)*exp and (1-scale)*exp. Fermi-like constraint "//&
"function",&
usage="SCALE 0.3", default_r_val=0.3_dp, required=.FALSE., error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="FERMI_EXP",&
description="Exponent in the fermi-like constraint function. ",&
usage="FERMI_EXP 2.63", default_r_val=2.63391_dp, required=.FALSE., error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
ENDIF
END SUBROUTINE create_constrain_exponents_section
! *****************************************************************************
!> \brief creates the multigrid
!> \param section ...

74
src/lri_environment_methods.F
View file

@ -339,7 +339,7 @@ CONTAINS
ra(:) = pbc(particle_set(iatom)%r, cell)
rb(:) = pbc(particle_set(jatom)%r, cell)
! calculate integrals (a,b,fa) and (a,b,fb)
! calculate integrals (aa,bb)
CALL lri_int_aabb(lriir%soaabb,obasa,obasb,rab,ra,rb,lri_env%debug,&
lriir%dmax_aabb,error)
@ -1723,4 +1723,76 @@ CONTAINS
END SUBROUTINE output_debug_info
! *****************************************************************************
!> \brief SVD of auxiliary overlap matrix
!> \param lri_env ...
!> \param qs_env ...
!> \param error ...
!> \note not called, routine for TESTING, to be deleted soon
! *****************************************************************************
SUBROUTINE svd_of_s(lri_env,qs_env,error)
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cp_error_type), INTENT(INOUT) :: error
CHARACTER(LEN=*), PARAMETER :: routineN = 'svd_of_s', &
routineP = moduleN//':'//routineN
INTEGER :: i, iac, iatom, ikind, ilist, &
info, jatom, jkind, &
jneighbor, lwork, m, n, &
nkind, nlist, nneighbor, stat
INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
LOGICAL :: failure
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sig, work
REAL(KIND=dp), DIMENSION(3) :: rab
REAL(KIND=dp), DIMENSION(:, :), POINTER :: s_temp, u, vt
TYPE(lri_int_type), POINTER :: lrii
TYPE(lri_list_type), POINTER :: lri_ints
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: nl_iterator
TYPE(neighbor_list_set_p_type), &
DIMENSION(:), POINTER :: soo_list
NULLIFY(soo_list, lrii, lri_ints, nl_iterator, s_temp, u , vt)
soo_list => lri_env%soo_list
lri_ints => lri_env%lri_ints
CALL neighbor_list_iterator_create(nl_iterator,soo_list)
DO WHILE (neighbor_list_iterate(nl_iterator)==0)
CALL get_iterator_info(nl_iterator,ikind=ikind,jkind=jkind,&
nlist=nlist,ilist=ilist,nnode=nneighbor,inode=jneighbor,&
iatom=iatom,jatom=jatom,r=rab)
iac = ikind + nkind*(jkind - 1)
lrii => lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
! re-inverse sinv to get S
CALL invmat(lrii%sinv,stat,error)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
m=SIZE(lrii%sinv,1)
n=m
lwork = m*(6+4*m)+m
!mxm matrix
ALLOCATE(s_temp(m,m),u(m,m),vt(m,m),sig(m),iwork(8*m),work(lwork))
s_temp(:,:) = lrii%sinv
! and change this back
CALL invmat(lrii%sinv,stat,error)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
! do SVD
CALL DGESDD('A',m,m,s_temp,m,sig,u,m,vt,m,work,lwork,iwork,info)
WRITE(*,*) "iatom, jatom", iatom, jatom
DO i=1,m
WRITE(*,*) i, sig(i)
ENDDO
DEALLOCATE(s_temp, u, vt, sig, iwork, work)
END DO
CALL neighbor_list_iterator_release(nl_iterator)
END SUBROUTINE
END MODULE lri_environment_methods

18
src/lri_environment_types.F
View file

@ -84,6 +84,8 @@ MODULE lri_environment_types
INTEGER :: nfa
! number of spherical fit basis functions (bi)
INTEGER :: nfb
! condition number of overlap matrix
REAL(KIND=dp) :: cond_num
! integrals (a,b,ai)
REAL(KIND=dp), DIMENSION(:,:,:), POINTER :: abaint
! integrals (a,b,bi)
@ -161,10 +163,6 @@ MODULE lri_environment_types
REAL(KIND=dp), DIMENSION(:,:), POINTER :: orb_ovlp
END TYPE lri_bas_overlap_type
TYPE lri_gcc_p_type
REAL(KIND=dp), DIMENSION(:,:,:), POINTER :: gcc_orig
END TYPE lri_gcc_p_type
! *****************************************************************************
TYPE lri_environment_type
@ -175,8 +173,6 @@ MODULE lri_environment_types
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: orb_basis
! lri (fit) basis set
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: ri_basis
! holds the original contraction coeff of the lri basis
TYPE(lri_gcc_p_type), DIMENSION(:), POINTER :: ri_gcc_orig
! orb_basis neighborlist, LRI integrals
TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: soo_list
! local RI integrals
@ -283,7 +279,6 @@ CONTAINS
NULLIFY(lri_env%orb_basis)
NULLIFY(lri_env%ri_basis)
NULLIFY(lri_env%ri_gcc_orig)
NULLIFY(lri_env%soo_list)
NULLIFY(lri_env%lri_ints)
@ -332,15 +327,6 @@ CONTAINS
DEALLOCATE(lri_env%ri_basis,stat=stat)
CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error)
END IF
IF(ASSOCIATED(lri_env%ri_gcc_orig)) THEN
nkind = SIZE(lri_env%ri_gcc_orig)
DO ikind=1,nkind
DEALLOCATE(lri_env%ri_gcc_orig(ikind)%gcc_orig,stat=stat)
CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error)
END DO
DEALLOCATE(lri_env%ri_gcc_orig,stat=stat)
CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error)
END IF
IF (ASSOCIATED(lri_env%soo_list)) THEN
DO i=1,SIZE(lri_env%soo_list)
CALL deallocate_neighbor_list_set(lri_env%soo_list(i)%neighbor_list_set)

676
src/lri_optimize_ri_basis.F
View file

@ -26,7 +26,11 @@ MODULE lri_optimize_ri_basis
cp_print_key_should_output,&
cp_print_key_unit_nr
USE cp_para_types, ONLY: cp_para_env_type
USE input_section_types, ONLY: section_vals_get_subs_vals,&
USE input_constants, ONLY: do_lri_opt_all,&
do_lri_opt_coeff,&
do_lri_opt_exps
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: default_path_length,&
@ -41,7 +45,11 @@ MODULE lri_optimize_ri_basis
lri_int_type,&
lri_list_type,&
lri_rhoab_type
USE mathconstants, ONLY: pi
USE lri_optimize_ri_basis_types, ONLY: create_lri_opt,&
deallocate_lri_opt,&
get_original_gcc,&
lri_opt_type,&
orthonormalize_gcc
USE memory_utilities, ONLY: reallocate
USE message_passing, ONLY: mp_sum
USE powell, ONLY: opt_state_type,&
@ -88,9 +96,7 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'optimize_lri_basis', &
routineP = moduleN//':'//routineN
INTEGER :: iunit, n10, nkind, stat
LOGICAL :: failure
REAL(KIND=dp), DIMENSION(:), POINTER :: x
INTEGER :: iunit, nkind
TYPE(atomic_kind_type), DIMENSION(:), &
POINTER :: atomic_kind_set
TYPE(cp_dbcsr_p_type), DIMENSION(:), &
@ -99,14 +105,15 @@ CONTAINS
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(lri_density_type), POINTER :: lri_density
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_opt_type), POINTER :: lri_opt
TYPE(opt_state_type) :: opt_state
TYPE(qs_rho_type), POINTER :: rho_struct
TYPE(section_vals_type), POINTER :: dft_section, input, &
lri_optbas_section
NULLIFY(atomic_kind_set, dft_section, lri_density, lri_env, &
lri_optbas_section, rho_struct)
NULLIFY(input, logger, para_env, x)
lri_opt,lri_optbas_section, rho_struct)
NULLIFY(input, logger, para_env)
CALL get_qs_env(qs_env,atomic_kind_set=atomic_kind_set,input=input,&
lri_env=lri_env,lri_density=lri_density,nkind=nkind,&
@ -127,44 +134,35 @@ CONTAINS
ENDIF
! *** initialization
CALL init_optimization(lri_env,lri_optbas_section,&
opt_state,x,nkind,iunit,error)
CALL create_lri_opt(lri_opt,error)
CALL init_optimization(lri_env,lri_opt,lri_optbas_section,&
opt_state,lri_opt%x,lri_opt%zet_init,nkind,iunit,error)
CALL calculate_lri_overlap_aabb(lri_env,qs_env,error)
n10 = MAX(opt_state%maxfun/100,1)
! *** start optimize
! *** ======================= START optimization =====================
opt_state%state = 0
DO
IF ( opt_state%state == 2 ) THEN
CALL calc_lri_integrals_get_objective(lri_env,lri_density,qs_env,&
opt_state,pmatrix,para_env,&
atomic_kind_set,nkind,x,error)
lri_opt,opt_state,pmatrix,para_env,&
atomic_kind_set,nkind,lri_opt%x,error)
ENDIF
IF ( opt_state%state == -1 ) EXIT
! *** ensure that exponents will be positive
x = SQRT(x)
CALL powell_optimize (opt_state%nvar, x, opt_state)
x= x**2._dp
IF ( opt_state%nf == 2 .AND. opt_state%state ==2 .AND. iunit > 0 ) THEN
WRITE(iunit,'(/," POWELL| Initial value of function",T61,F20.10)') opt_state%f
END IF
IF ( MOD(opt_state%nf,n10) == 0 .AND. opt_state%nf > 1 .AND. iunit > 0 ) THEN
WRITE(iunit,'(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
INT(REAL(opt_state%nf,dp)/REAL(opt_state%maxfun,dp)*100._dp), opt_state%fopt
END IF
CALL powell_optimize (opt_state%nvar, lri_opt%x, opt_state)
CALL update_exponents(lri_env,lri_opt,lri_opt%x,lri_opt%zet_init,nkind,error)
CALL print_optimization_update(opt_state,lri_opt,iunit)
ENDDO
! *** ======================= END optimization =======================
! *** get final optimized parameters
opt_state%state = 8
CALL powell_optimize (opt_state%nvar, lri_opt%x, opt_state)
CALL update_exponents(lri_env,lri_opt,lri_opt%x,lri_opt%zet_init,nkind,error)
x = SQRT(x)
CALL powell_optimize (opt_state%nvar, x, opt_state)
x= x**2._dp
CALL write_optimized_lri_basis(lri_env,dft_section,nkind,x,&
CALL write_optimized_lri_basis(lri_env,dft_section,nkind,lri_opt,&
atomic_kind_set,error)
IF ( iunit > 0 ) THEN
@ -176,103 +174,56 @@ CONTAINS
CALL cp_print_key_finished_output(iunit,logger,input,&
"PRINT%PROGRAM_RUN_INFO", error=error)
DEALLOCATE (x,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
CALL deallocate_lri_opt(lri_opt,error)
END SUBROUTINE optimize_lri_basis
! *****************************************************************************
!> \brief calculate objective only...
!> \param qs_env ...
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
!> \note not called; testing only, to be removed later
! *****************************************************************************
SUBROUTINE calc_objective_only(qs_env,error)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cp_error_type), INTENT(INOUT) :: error
CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_objective_only', &
routineP = moduleN//':'//routineN
INTEGER :: iunit, nkind
REAL(KIND=dp) :: objf
TYPE(atomic_kind_type), DIMENSION(:), &
POINTER :: atomic_kind_set
TYPE(cp_dbcsr_p_type), DIMENSION(:), &
POINTER :: pmatrix
TYPE(cp_logger_type), POINTER :: logger
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(lri_density_type), POINTER :: lri_density
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(qs_rho_type), POINTER :: rho_struct
TYPE(section_vals_type), POINTER :: input
NULLIFY(atomic_kind_set, lri_density, lri_env, para_env, rho_struct)
CALL get_qs_env(qs_env,atomic_kind_set=atomic_kind_set,input=input,&
lri_env=lri_env,lri_density=lri_density,nkind=nkind,&
rho=rho_struct,para_env=para_env,error=error)
CALL qs_rho_get(rho_struct, rho_ao=pmatrix, error=error)
objf=0._dp
logger => cp_error_get_logger(error)
iunit=cp_print_key_unit_nr(logger,input,"PRINT%PROGRAM_RUN_INFO",&
extension=".opt",error=error)
CALL calculate_lri_overlap_aabb(lri_env,qs_env,error)
CALL calculate_lri_integrals(lri_env,qs_env,calculate_forces=.FALSE.,error=error)
CALL calculate_avec(lri_env,lri_density,qs_env,pmatrix,error=error)
CALL calculate_objective(lri_env,lri_density,pmatrix,para_env,objf,error)
IF ( iunit > 0 ) THEN
WRITE(iunit,'("OBJF",T71,F20.15)') objf
ENDIF
CALL cp_print_key_finished_output(iunit,logger,input,&
"PRINT%PROGRAM_RUN_INFO", error=error)
END SUBROUTINE calc_objective_only
! *****************************************************************************
!> \brief calculates the lri integrals and coefficients with the new exponents
!> of the lri basis sets and calculates the objective function
!> \param lri_env ...
!> \brief initialize optimization parameter
!> \param lri_env lri environment
!> \param lri_opt optimization environment
!> \param lri_optbas_section ...
!> \param opt_state ...
!> \param x parameters to be optimized, i.e. exponents of the lri basis set
!> \param nkind ...
!> \param iunit ...
!> \param opt_state state of the optimizer
!> \param x parameters to be optimized, i.e. exponents and contraction coeffs
!> of the lri basis set
!> \param zet_init initial values of the exponents
!> \param nkind number of atom kinds
!> \param iunit output unit
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE init_optimization(lri_env,lri_optbas_section,opt_state,x,nkind,&
iunit,error)
SUBROUTINE init_optimization(lri_env,lri_opt,lri_optbas_section,opt_state,&
x,zet_init,nkind,iunit,error)
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_opt_type), POINTER :: lri_opt
TYPE(section_vals_type), POINTER :: lri_optbas_section
TYPE(opt_state_type) :: opt_state
REAL(KIND=dp), DIMENSION(:), POINTER :: x
REAL(KIND=dp), DIMENSION(:), POINTER :: x, zet_init
INTEGER, INTENT(IN) :: nkind, iunit
TYPE(cp_error_type), INTENT(INOUT) :: error
CHARACTER(LEN=*), PARAMETER :: routineN = 'init_optimization', &
routineP = moduleN//':'//routineN
INTEGER :: ikind, iset, n, nset, stat
INTEGER, DIMENSION(:), POINTER :: npgf
INTEGER :: ikind, iset, ishell, n, nset, &
stat
INTEGER, DIMENSION(:), POINTER :: npgf, nshell
LOGICAL :: failure
REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
REAL(KIND=dp), DIMENSION(:, :, :), &
POINTER :: gcc_orig
TYPE(gto_basis_set_type), POINTER :: fbas
NULLIFY(fbas,npgf,zet)
NULLIFY(fbas, gcc_orig, npgf, nshell, zet)
failure = .FALSE.
ALLOCATE(lri_env%ri_gcc_orig(nkind),STAT=stat)
ALLOCATE(lri_opt%ri_gcc_orig(nkind),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
CALL section_vals_val_get(lri_optbas_section,"ACCURACY",&
r_val=opt_state%rhoend, error=error)
CALL section_vals_val_get(lri_optbas_section,"STEP_SIZE",&
r_val=opt_state%rhobeg, error=error)
CALL section_vals_val_get(lri_optbas_section,"MAX_FUN",&
i_val=opt_state%maxfun, error=error)
! *** get parameters
CALL get_optimization_parameter(lri_env,lri_opt,lri_optbas_section,&
opt_state,error)
opt_state%nvar =0
opt_state%nf = 0
@ -280,25 +231,68 @@ CONTAINS
opt_state%unit = iunit
! *** init exponents
n = 0
DO ikind=1,nkind
fbas => lri_env%ri_basis(ikind)%gto_basis_set
CALL get_gto_basis_set(gto_basis_set=fbas,&
npgf=npgf,nset=nset,zet=zet)
DO iset =1,nset
opt_state%nvar = opt_state%nvar + npgf(iset)
CALL reallocate(x,1,opt_state%nvar)
x(n+1:n+npgf(iset)) = zet(:,iset)
n = n + npgf(iset)
IF(lri_opt%opt_exps) THEN
n = 0
DO ikind=1,nkind
fbas => lri_env%ri_basis(ikind)%gto_basis_set
CALL get_gto_basis_set(gto_basis_set=fbas,&
npgf=npgf,nset=nset,zet=zet)
DO iset =1,nset
IF(lri_opt%use_geometric_seq.AND.npgf(iset)>2) THEN
opt_state%nvar = opt_state%nvar + 2
CALL reallocate(x,1,opt_state%nvar)
x(n+1) = MAXVAL(zet(1:npgf(iset),iset))
x(n+2) = MINVAL(zet(1:npgf(iset),iset))
n = n + 2
ELSE
opt_state%nvar = opt_state%nvar + npgf(iset)
CALL reallocate(x,1,opt_state%nvar)
x(n+1:n+npgf(iset)) = zet(1:npgf(iset),iset)
n = n + npgf(iset)
ENDIF
lri_opt%nexp=lri_opt%nexp + npgf(iset)
ENDDO
ENDDO
ENDDO
! *** constraints on exponents
IF(lri_opt%use_constraints) THEN
ALLOCATE(zet_init(SIZE(x)),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
zet_init(:)=x
ELSE
x(:)=SQRT(x)
ENDIF
ENDIF
! *** get the original gcc without normalization factor
DO ikind=1,nkind
fbas => lri_env%ri_basis(ikind)%gto_basis_set
CALL get_original_gcc(lri_env%ri_gcc_orig(ikind)%gcc_orig,fbas,error)
CALL get_original_gcc(lri_opt%ri_gcc_orig(ikind)%gcc_orig,fbas,&
lri_opt,error)
ENDDO
! *** init coefficients
IF(lri_opt%opt_coeffs) THEN
DO ikind=1,nkind
fbas => lri_env%ri_basis(ikind)%gto_basis_set
gcc_orig => lri_opt%ri_gcc_orig(ikind)%gcc_orig
CALL get_gto_basis_set(gto_basis_set=fbas,&
npgf=npgf,nset=nset,nshell=nshell,zet=zet)
! *** Gram Schmidt orthonormalization
CALL orthonormalize_gcc(gcc_orig,fbas,lri_opt,error)
n=opt_state%nvar
DO iset=1,nset
DO ishell=1,nshell(iset)
opt_state%nvar = opt_state%nvar + npgf(iset)
CALL reallocate(x,1,opt_state%nvar)
x(n+1:n+npgf(iset))=gcc_orig(1:npgf(iset),ishell,iset)
lri_opt%ncoeff = lri_opt%ncoeff + npgf(iset)
n = n + npgf(iset)
ENDDO
ENDDO
ENDDO
ENDIF
IF(iunit > 0 ) THEN
WRITE(iunit,'(/," POWELL| Accuracy",T69,ES12.5)') opt_state%rhoend
WRITE(iunit,'(" POWELL| Initial step size",T69,ES12.5)') opt_state%rhobeg
@ -310,28 +304,266 @@ CONTAINS
END SUBROUTINE init_optimization
! *****************************************************************************
!> \brief read input for optimization
!> \param lri_env lri environment
!> \param lri_opt optimization environmnet
!> \param lri_optbas_section ...
!> \param opt_state state of the optimizer
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE get_optimization_parameter(lri_env,lri_opt,lri_optbas_section,&
opt_state,error)
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_opt_type), POINTER :: lri_opt
TYPE(section_vals_type), POINTER :: lri_optbas_section
TYPE(opt_state_type) :: opt_state
TYPE(cp_error_type), INTENT(INOUT) :: error
CHARACTER(LEN=*), PARAMETER :: routineN = 'get_optimization_parameter', &
routineP = moduleN//':'//routineN
INTEGER :: degree_freedom
TYPE(section_vals_type), POINTER :: constrain_exp_section
NULLIFY(constrain_exp_section)
! *** parameter for POWELL optimizer
CALL section_vals_val_get(lri_optbas_section,"ACCURACY",&
r_val=opt_state%rhoend, error=error)
CALL section_vals_val_get(lri_optbas_section,"STEP_SIZE",&
r_val=opt_state%rhobeg, error=error)
CALL section_vals_val_get(lri_optbas_section,"MAX_FUN",&
i_val=opt_state%maxfun, error=error)
! *** parameters which are optimized, i.e. exps or coeff or both
CALL section_vals_val_get(lri_optbas_section,"DEGREES_OF_FREEDOM",&
i_val=degree_freedom, error=error)
SELECT CASE(degree_freedom)
CASE(do_lri_opt_all)
lri_opt%opt_coeffs=.TRUE.
lri_opt%opt_exps=.TRUE.
CASE(do_lri_opt_coeff)
lri_opt%opt_coeffs=.TRUE.
CASE(do_lri_opt_exps)
lri_opt%opt_exps=.TRUE.
CASE DEFAULT
CALL cp_assert(.FALSE.,cp_fatal_level,cp_assertion_failed,&
routineP,"No initialization available?????",&
only_ionode=.TRUE.)
END SELECT
! *** restraint
CALL section_vals_val_get(lri_optbas_section,"USE_CONDITION_NUMBER",&
l_val=lri_opt%use_condition_number,error=error)
CALL section_vals_val_get(lri_optbas_section,"CONDITION_WEIGHT",&
r_val=lri_opt%cond_weight,error=error)
CALL section_vals_val_get(lri_optbas_section,"GEOMETRIC_SEQUENCE",&
l_val=lri_opt%use_geometric_seq,error=error)
! *** get constraint info
constrain_exp_section => section_vals_get_subs_vals(lri_optbas_section,&
"CONSTRAIN_EXPONENTS",error=error)
CALL section_vals_get(constrain_exp_section,explicit=lri_opt%use_constraints,&
error=error)
IF(lri_opt%use_constraints) THEN
CALL section_vals_val_get(constrain_exp_section,"SCALE",&
r_val=lri_opt%scale_exp, error=error)
CALL section_vals_val_get(constrain_exp_section,"FERMI_EXP",&
r_val=lri_opt%fermi_exp, error=error)
ENDIF
END SUBROUTINE get_optimization_parameter
! *****************************************************************************
!> \brief update exponents after optimization step
!> \param lri_env lri environment
!> \param lri_opt optimization environment
!> \param x optimization parameters
!> \param zet_init initial values of the exponents
!> \param nkind number of atomic kinds
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE update_exponents(lri_env,lri_opt,x,zet_init,nkind,error)
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_opt_type), POINTER :: lri_opt
REAL(KIND=dp), DIMENSION(:), POINTER :: x, zet_init
INTEGER, INTENT(IN) :: nkind
TYPE(cp_error_type), INTENT(INOUT) :: error
CHARACTER(LEN=*), PARAMETER :: routineN = 'update_exponents', &
routineP = moduleN//':'//routineN
INTEGER :: ikind, iset, ishell, n, nset, &
nvar_exp, stat
INTEGER, DIMENSION(:), POINTER :: npgf, nshell
LOGICAL :: failure
REAL(KIND=dp) :: zet_max, zet_min
REAL(KIND=dp), DIMENSION(:), POINTER :: zet, zet_trans
REAL(KIND=dp), DIMENSION(:, :, :), &
POINTER :: gcc_orig
TYPE(gto_basis_set_type), POINTER :: fbas
NULLIFY(fbas,gcc_orig,npgf,nshell,zet_trans, zet)
failure = .FALSE.
! nvar_exp: number of exponents that are variables
nvar_exp= SIZE(x) - lri_opt%ncoeff
ALLOCATE(zet_trans(nvar_exp),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
! *** update exponents
IF(lri_opt%opt_exps) THEN
IF (lri_opt%use_constraints) THEN
zet => x(1:nvar_exp)
CALL transfer_exp(lri_opt,zet,zet_init,zet_trans,nvar_exp,error)
ELSE
zet_trans(:) = x(1:nvar_exp)**2.0_dp
ENDIF
n=0
DO ikind =1,nkind
fbas => lri_env%ri_basis(ikind)%gto_basis_set
CALL get_gto_basis_set(gto_basis_set=fbas,npgf=npgf,nset=nset)
DO iset=1,nset
IF(lri_opt%use_geometric_seq.AND.npgf(iset)>2) THEN
zet_max = MAXVAL(zet_trans(n+1:n+2))
zet_min = MINVAL(zet_trans(n+1:n+2))
zet => fbas%zet(1:npgf(iset),iset)
CALL geometric_progression(zet,zet_max,zet_min,npgf(iset))
n = n + 2
ELSE
fbas%zet(1:npgf(iset),iset) = zet_trans(n+1:n+npgf(iset))
n = n + npgf(iset)
ENDIF
ENDDO
ENDDO
ENDIF
! *** update coefficients
IF(lri_opt%opt_coeffs) THEN
n=nvar_exp
DO ikind=1,nkind
fbas => lri_env%ri_basis(ikind)%gto_basis_set
gcc_orig => lri_opt%ri_gcc_orig(ikind)%gcc_orig
CALL get_gto_basis_set(gto_basis_set=fbas,&
nshell=nshell,npgf=npgf,nset=nset)
DO iset=1,nset
DO ishell=1,nshell(iset)
gcc_orig(1:npgf(iset),ishell,iset) = x(n+1:n+npgf(iset))
n = n + npgf(iset)
ENDDO
ENDDO
! *** Gram Schmidt orthonormalization
CALL orthonormalize_gcc(gcc_orig,fbas,lri_opt,error)
ENDDO
ENDIF
DEALLOCATE(zet_trans,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END SUBROUTINE update_exponents
! *****************************************************************************
!> \brief employ Fermi constraint, transfer exponents
!> \param lri_opt optimization environment
!> \param zet untransferred exponents
!> \param zet_init intial value of the eponents
!> \param zet_trans transferred exponents
!> \param nvar number of optimized exponents
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE transfer_exp(lri_opt,zet,zet_init,zet_trans,nvar,error)
TYPE(lri_opt_type), POINTER :: lri_opt
REAL(KIND=dp), DIMENSION(:), POINTER :: zet, zet_init, zet_trans
INTEGER, INTENT(IN) :: nvar
TYPE(cp_error_type), INTENT(INOUT) :: error
CHARACTER(LEN=*), PARAMETER :: routineN = 'transfer_exp', &
routineP = moduleN//':'//routineN
INTEGER :: stat
LOGICAL :: failure
REAL(KIND=dp) :: a
REAL(KIND=dp), DIMENSION(:), POINTER :: zet_max, zet_min
failure = .FALSE.
ALLOCATE(zet_max(nvar),zet_min(nvar),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
zet_min(:) = zet_init(:)*(1.0_dp-lri_opt%scale_exp)
zet_max(:) = zet_init(:)*(1.0_dp+lri_opt%scale_exp)
a=lri_opt%fermi_exp
zet_trans= zet_min + (zet_max-zet_min)/(1+EXP(-a*(zet-zet_init)))
DEALLOCATE(zet_max,zet_min,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END SUBROUTINE transfer_exp
! *****************************************************************************
!> \brief complete geometric sequence
!> \param zet all exponents of the set
!> \param zet_max maximal exponent of the set
!> \param zet_min minimal exponent of the set
!> \param nexp number of exponents of the set
! *****************************************************************************
SUBROUTINE geometric_progression(zet,zet_max,zet_min,nexp)
REAL(KIND=dp), DIMENSION(:), POINTER :: zet
REAL(KIND=dp), INTENT(IN) :: zet_max, zet_min
INTEGER, INTENT(IN) :: nexp
CHARACTER(LEN=*), PARAMETER :: routineN = 'geometric_progression', &
routineP = moduleN//':'//routineN
INTEGER :: i, n
REAL(KIND=dp) :: q
n=nexp-1
q=(zet_min/zet_max)**(1._dp/REAL(n,dp))
DO i=1,nexp
zet(i)=zet_max*q**(i-1)
ENDDO
END SUBROUTINE geometric_progression
! *****************************************************************************
!> \brief calculates the lri integrals and coefficients with the new exponents
!> of the lri basis sets and calculates the objective function
!> \param lri_env ...
!> \param lri_env lri environment
!> \param lri_density ...
!> \param qs_env ...
!> \param opt_state ...
!> \param lri_opt optimization enviroment
!> \param opt_state state of the optimizer
!> \param pmatrix density matrix
!> \param para_env ...
!> \param atomic_kind_set ...
!> \param nkind ...
!> \param x parameters to be optimized, i.e. exponents of the lri basis set
!> \param nkind number of atomic kinds
!> \param x parameters to be optimized, i.e. exponents and contraction coeffs
!> of the lri basis set
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE calc_lri_integrals_get_objective(lri_env,lri_density,qs_env,&
opt_state,pmatrix,para_env,&
lri_opt,opt_state,pmatrix,para_env,&
atomic_kind_set,nkind,x,error)
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_density_type), POINTER :: lri_density
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(lri_opt_type), POINTER :: lri_opt
TYPE(opt_state_type) :: opt_state
TYPE(cp_dbcsr_p_type), DIMENSION(:), &
POINTER :: pmatrix
@ -346,31 +578,27 @@ CONTAINS
routineN = 'calc_lri_integrals_get_objective', &
routineP = moduleN//':'//routineN
INTEGER :: ikind, ipgf, iset, n, nset
INTEGER :: ikind, nset
INTEGER, DIMENSION(:), POINTER :: npgf
TYPE(gto_basis_set_type), POINTER :: fbas
NULLIFY(fbas,npgf)
!*** sort the exponents and build new transformation matrices sphi
n = 0
!*** build new transformation matrices sphi with new exponents
DO ikind =1,nkind
fbas => lri_env%ri_basis(ikind)%gto_basis_set
CALL get_gto_basis_set(gto_basis_set=fbas,npgf=npgf,nset=nset)
DO iset =1,nset
DO ipgf=1,npgf(iset)
fbas%zet(ipgf,iset) = x(n+ipgf)
ENDDO
n = n + npgf(iset)
ENDDO
!build new sphi
fbas%gcc = lri_env%ri_gcc_orig(ikind)%gcc_orig
fbas%gcc = lri_opt%ri_gcc_orig(ikind)%gcc_orig
CALL init_orb_basis_set(fbas,error)
ENDDO
CALL lri_basis_init(lri_env,atomic_kind_set,error)
CALL calculate_lri_integrals(lri_env,qs_env,calculate_forces=.FALSE.,error=error)
CALL calculate_avec(lri_env,lri_density,qs_env,pmatrix,error=error)
CALL calculate_objective(lri_env,lri_density,pmatrix,para_env,&
IF(lri_opt%use_condition_number) THEN
CALL get_condition_number_of_overlap(lri_env,error)
ENDIF
CALL calculate_objective(lri_env,lri_density,lri_opt,pmatrix,para_env,&
opt_state%f,error)
@ -380,19 +608,21 @@ CONTAINS
!> \brief calculates the objective function defined as integral of the square
!> of rhoexact - rhofit, i.e. integral[(rhoexact-rhofit)**2]
!> rhoexact is the exact pair density and rhofit the lri pair density
!> \param lri_env ...
!> \param lri_env lri environment
!> \param lri_density ...
!> \param lri_opt optimization environment
!> \param pmatrix density matrix
!> \param para_env ...
!> \param fobj objective function
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE calculate_objective(lri_env,lri_density,pmatrix,para_env,&
SUBROUTINE calculate_objective(lri_env,lri_density,lri_opt,pmatrix,para_env,&
fobj,error)
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(lri_density_type), POINTER :: lri_density
TYPE(lri_opt_type), POINTER :: lri_opt
TYPE(cp_dbcsr_p_type), DIMENSION(:), &
POINTER :: pmatrix
TYPE(cp_para_env_type), POINTER :: para_env
@ -429,6 +659,7 @@ CONTAINS
nkind = lri_env%lri_ints%nkind
nspin = SIZE(pmatrix)
fobj = 0._dp
lri_opt%rho_diff = 0._dp
DO ispin = 1, nspin
@ -532,7 +763,12 @@ CONTAINS
obj_ab = 2.0_dp*(rhoexact_sq - 2._dp*rhomix + rhofit_sq)
ENDIF
fobj = fobj + obj_ab
IF(lri_opt%use_condition_number) THEN
fobj = fobj + obj_ab + lri_opt%cond_weight*LOG(lrii%cond_num)
lri_opt%rho_diff = lri_opt%rho_diff + obj_ab
ELSE
fobj = fobj + obj_ab
ENDIF
ENDDO
@ -548,69 +784,149 @@ CONTAINS
END SUBROUTINE calculate_objective
! *****************************************************************************
!> \brief primitive Cartesian Gaussian functions are normalized. The normalization
!> factor is included in the Gaussian contraction coefficients.
!> Division by this factor to get the original gcc.
!> \param gcc_orig ...
!> \param gto_basis_set ...
!> \brief get condition number of overlap matrix
!> \param lri_env lri environment
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE get_original_gcc(gcc_orig,gto_basis_set,error)
REAL(KIND=dp), DIMENSION(:, :, :), &
POINTER :: gcc_orig
TYPE(gto_basis_set_type), POINTER :: gto_basis_set
SUBROUTINE get_condition_number_of_overlap(lri_env,error)
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(LEN=*), PARAMETER :: routineN = 'get_original_gcc', &
CHARACTER(LEN=*), PARAMETER :: &
routineN = 'get_condition_number_of_overlap', &
routineP = moduleN//':'//routineN
INTEGER :: ipgf, iset, ishell, l, &
maxpgf, maxshell, nset, stat
INTEGER :: handle, iac, iatom, ikind, ilist, info, jatom, jkind, &
jneighbor, lwork, nfa, nfb, nkind, nlist, nn, nneighbor, stat
LOGICAL :: failure
REAL(KIND=dp) :: expzet, gcca, prefac, zeta
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag, off_diag, tau
REAL(KIND=dp), DIMENSION(:), POINTER :: work
REAL(KIND=dp), DIMENSION(:, :), POINTER :: smat
TYPE(lri_int_type), POINTER :: lrii
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: nl_iterator
TYPE(neighbor_list_set_p_type), &
DIMENSION(:), POINTER :: soo_list
failure = .FALSE.
maxpgf = SIZE(gto_basis_set%gcc,1)
maxshell = SIZE(gto_basis_set%gcc,2)
nset = SIZE(gto_basis_set%gcc,3)
failure = .FALSE.
CALL timeset(routineN,handle)
NULLIFY(lrii, nl_iterator, smat, soo_list)
ALLOCATE(gcc_orig(maxpgf,maxshell,nset),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
DO iset=1,gto_basis_set%nset
DO ishell=1,gto_basis_set%nshell(iset)
l = gto_basis_set%l(ishell,iset)
expzet = 0.25_dp*REAL(2*l + 3,dp)
prefac = 2.0_dp**l*(2.0_dp/pi)**0.75_dp
DO ipgf=1,gto_basis_set%npgf(iset)
gcca = gto_basis_set%gcc(ipgf,ishell,iset)
zeta = gto_basis_set%zet(ipgf,iset)
gcc_orig(ipgf,ishell,iset) = gcca/(prefac*zeta**expzet)
END DO
END DO
END DO
soo_list => lri_env%soo_list
END SUBROUTINE get_original_gcc
nkind = lri_env%lri_ints%nkind
CALL neighbor_list_iterator_create(nl_iterator,soo_list)
DO WHILE (neighbor_list_iterate(nl_iterator)==0)
CALL get_iterator_info(nl_iterator,ikind=ikind,jkind=jkind,iatom=iatom,&
jatom=jatom,nlist=nlist,ilist=ilist,nnode=nneighbor,inode=jneighbor)
iac = ikind + nkind*(jkind - 1)
IF(.NOT.ASSOCIATED(lri_env%lri_ints%lri_atom(iac)%lri_node)) CYCLE
lrii => lri_env%lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
nfa=lrii%nfa
nfb=lrii%nfb
nn = nfa + nfb
! build the overlap matrix
IF(iatom == jatom) THEN
ALLOCATE(smat(nfa,nfa),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ELSE
ALLOCATE(smat(nn,nn),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ENDIF
smat(1:nfa,1:nfa) = lri_env%bas_ovlp(ikind)%ri_ovlp(1:nfa,1:nfa)
IF(iatom /= jatom) THEN
nn = nfa+nfb
smat(1:nfa,nfa+1:nn) = lrii%sab(1:nfa,1:nfb)
smat(nfa+1:nn,1:nfa) = TRANSPOSE(lrii%sab(1:nfa,1:nfb))
smat(nfa+1:nn,nfa+1:nn) = lri_env%bas_ovlp(jkind)%ri_ovlp(1:nfb,1:nfb)
ENDIF
IF(iatom==jatom) nn=nfa
ALLOCATE(diag(nn),off_diag(nn-1),tau(nn-1),work(1),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
diag=0.0_dp
off_diag=0.0_dp
tau=0.0_dp
work=0.0_dp
lwork = -1
! get lwork
CALL DSYTRD('U',nn,smat,nn,diag,off_diag,tau,work,lwork,info)
lwork=INT(work(1))
CALL reallocate(work,1,lwork)
! get the eigenvalues
CALL DSYTRD('U',nn,smat,nn,diag,off_diag,tau,work,lwork,info)
CALL DSTERF(nn,diag,off_diag,info)
lrii%cond_num=MAXVAL(ABS(diag))/MINVAL(ABS(diag))
DEALLOCATE(diag,off_diag,smat,tau,work)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
END DO
CALL neighbor_list_iterator_release(nl_iterator)
CALL timestop(handle)
END SUBROUTINE get_condition_number_of_overlap
! *****************************************************************************
!> \brief print recent information on optimization
!> \param opt_state state of the optimizer
!> \param lri_opt optimization environment
!> \param iunit ...
! *****************************************************************************
SUBROUTINE print_optimization_update(opt_state,lri_opt,iunit)
TYPE(opt_state_type) :: opt_state
TYPE(lri_opt_type), POINTER :: lri_opt
INTEGER, INTENT(IN) :: iunit
CHARACTER(LEN=*), PARAMETER :: routineN = 'print_optimization_update', &
routineP = moduleN//':'//routineN
INTEGER :: n10
n10 = MAX(opt_state%maxfun/100,1)
IF ( opt_state%nf == 2 .AND. opt_state%state ==2 .AND. iunit > 0 ) THEN
WRITE(iunit,'(/," POWELL| Initial value of function",T61,F20.10)') opt_state%f
END IF
IF ( MOD(opt_state%nf,n10) == 0 .AND. opt_state%nf > 1 .AND. iunit > 0 ) THEN
WRITE(iunit,'(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
INT(REAL(opt_state%nf,dp)/REAL(opt_state%maxfun,dp)*100._dp), opt_state%fopt
ENDIF
IF(lri_opt%use_condition_number) THEN
IF ( MOD(opt_state%nf,n10) == 0 .AND. opt_state%nf > 1 .AND. iunit > 0 ) THEN
WRITE(iunit,'(" POWELL| Recent value of function without condition nr.",T61,F20.10)') &
lri_opt%rho_diff
ENDIF
ENDIF
END SUBROUTINE print_optimization_update
! *****************************************************************************
!> \brief write optimized LRI basis set to file
!> \param lri_env ...
!> \param dft_section ...
!> \param nkind ...
!> \param xopt ...
!> \param lri_opt ...
!> \param atomic_kind_set ...
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE write_optimized_lri_basis(lri_env,dft_section,nkind,xopt,&
SUBROUTINE write_optimized_lri_basis(lri_env,dft_section,nkind,lri_opt,&
atomic_kind_set,error)
TYPE(lri_environment_type), POINTER :: lri_env
TYPE(section_vals_type), POINTER :: dft_section
INTEGER, INTENT(IN) :: nkind
REAL(KIND=dp), DIMENSION(:), POINTER :: xopt
TYPE(lri_opt_type), POINTER :: lri_opt
TYPE(atomic_kind_type), DIMENSION(:), &
POINTER :: atomic_kind_set
TYPE(cp_error_type), INTENT(INOUT) :: error
@ -620,7 +936,7 @@ CONTAINS
CHARACTER(LEN=default_path_length) :: filename
INTEGER :: cc_l, ikind, ipgf, iset, &
ishell, n, nset, output_file
ishell, nset, output_file
INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
INTEGER, DIMENSION(:, :), POINTER :: l
REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
@ -632,20 +948,6 @@ CONTAINS
NULLIFY(fbas,gcc_orig,l,lmax,lmin,logger,npgf,nshell,print_key,zet)
!*** sort the exponents
n = 0
DO ikind =1,nkind
fbas => lri_env%ri_basis(ikind)%gto_basis_set
CALL get_gto_basis_set(gto_basis_set=fbas,&
npgf=npgf,nset=nset,zet=zet)
DO iset =1,nset
DO ipgf=1,npgf(iset)
zet(ipgf,iset) = xopt(n+ipgf)
ENDDO
n = n + npgf(iset)
ENDDO
ENDDO
!*** do the printing
print_key => section_vals_get_subs_vals(dft_section,&
"PRINT%OPTIMIZE_LRI_BASIS",&
@ -670,7 +972,7 @@ CONTAINS
DO ikind =1, nkind
fbas => lri_env%ri_basis(ikind)%gto_basis_set
gcc_orig => lri_env%ri_gcc_orig(ikind)%gcc_orig
gcc_orig => lri_opt%ri_gcc_orig(ikind)%gcc_orig
CALL get_gto_basis_set(gto_basis_set=fbas,&
l=l, lmax=lmax, lmin=lmin,&
npgf=npgf,nshell=nshell,&

285
src/lri_optimize_ri_basis_types.F Normal file
View file

@ -0,0 +1,285 @@
!-----------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright (C) 2000 - 2014 CP2K developers group !
!-----------------------------------------------------------------------------!
! *****************************************************************************
!> \brief sets the environment for optimization of exponents and contraction
!> coefficients of the lri auxiliary
!> lri : local resolution of the identity
!> \par History
!> created Dorothea Golze [12.2014]
!> \authors Dorothea Golze
! *****************************************************************************
MODULE lri_optimize_ri_basis_types
USE basis_set_types, ONLY: get_gto_basis_set,&
gto_basis_set_type
USE kinds, ONLY: dp
USE mathconstants, ONLY: pi
#include "./common/cp_common_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_optimize_ri_basis_types'
PUBLIC :: lri_opt_type
PUBLIC :: create_lri_opt, deallocate_lri_opt, get_original_gcc,&
orthonormalize_gcc
! *****************************************************************************
TYPE lri_gcc_p_type
! gcc without normalization factor
REAL(KIND=dp), DIMENSION(:,:,:), POINTER :: gcc_orig
END TYPE lri_gcc_p_type
TYPE lri_subset_type
! amount of l quantum numbers per set
INTEGER :: nl
! number of contraction per l quantum number for a given set
INTEGER, DIMENSION(:), POINTER :: ncont_l
END TYPE lri_subset_type
TYPE lri_opt_type
LOGICAL :: opt_exps
LOGICAL :: opt_coeffs
LOGICAL :: use_condition_number
LOGICAL :: use_geometric_seq
LOGICAL :: use_constraints
INTEGER :: nexp
INTEGER :: ncoeff
REAL(KIND=dp) :: cond_weight
REAL(KIND=dp) :: scale_exp
REAL(KIND=dp) :: fermi_exp
REAL(KIND=dp) :: rho_diff
! array holding the variables that are optimized
REAL(KIND=dp), DIMENSION(:), POINTER :: x
! intial exponents
REAL(KIND=dp), DIMENSION(:), POINTER :: zet_init
! holds the original contraction coeff of the lri basis
TYPE(lri_gcc_p_type), DIMENSION(:), POINTER :: ri_gcc_orig
TYPE(lri_subset_type), DIMENSION(:), POINTER :: subset
END TYPE lri_opt_type
! *****************************************************************************
CONTAINS
! *****************************************************************************
!> \brief creates lri_opt
!> \param lri_opt optimization environment
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE create_lri_opt(lri_opt,error)
TYPE(lri_opt_type), POINTER :: lri_opt
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'create_lri_opt', &
routineP = moduleN//':'//routineN
INTEGER :: stat
LOGICAL :: failure
ALLOCATE(lri_opt,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
NULLIFY(lri_opt%ri_gcc_orig)
NULLIFY(lri_opt%subset)
NULLIFY(lri_opt%x)
NULLIFY(lri_opt%zet_init)
lri_opt%opt_exps = .FALSE.
lri_opt%opt_coeffs = .FALSE.
lri_opt%use_condition_number = .FALSE.
lri_opt%use_geometric_seq = .FALSE.
lri_opt%use_constraints = .FALSE.
lri_opt%nexp=0
lri_opt%ncoeff=0
END SUBROUTINE create_lri_opt
! *****************************************************************************
!> \brief deallocates lri_opt
!> \param lri_opt optimization environment
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE deallocate_lri_opt(lri_opt,error)
TYPE(lri_opt_type), POINTER :: lri_opt
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_lri_opt', &
routineP = moduleN//':'//routineN
INTEGER :: i, stat
LOGICAL :: failure
failure=.FALSE.
IF(ASSOCIATED(lri_opt)) THEN
IF(ASSOCIATED(lri_opt%subset)) THEN
DO i=1,SIZE(lri_opt%subset)
DEALLOCATE(lri_opt%subset(i)%ncont_l,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ENDDO
DEALLOCATE(lri_opt%subset,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ENDIF
IF(ASSOCIATED(lri_opt%x)) THEN
DEALLOCATE(lri_opt%x,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ENDIF
IF(ASSOCIATED(lri_opt%zet_init)) THEN
DEALLOCATE(lri_opt%zet_init,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ENDIF
IF(ASSOCIATED(lri_opt%ri_gcc_orig)) THEN
DO i=1,SIZE(lri_opt%ri_gcc_orig)
DEALLOCATE(lri_opt%ri_gcc_orig(i)%gcc_orig,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ENDDO
DEALLOCATE(lri_opt%ri_gcc_orig,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ENDIF
DEALLOCATE(lri_opt,STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ENDIF
END SUBROUTINE deallocate_lri_opt
! *****************************************************************************
!> \brief primitive Cartesian Gaussian functions are normalized. The normalization
!> factor is included in the Gaussian contraction coefficients.
!> Division by this factor to get the original gcc.
!> \param gcc_orig original contraction coefficient
!> \param gto_basis_set gaussian type basis set
!> \param lri_opt optimization environment
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE get_original_gcc(gcc_orig,gto_basis_set,lri_opt,error)
REAL(KIND=dp), DIMENSION(:, :, :), &
POINTER :: gcc_orig
TYPE(gto_basis_set_type), POINTER :: gto_basis_set
TYPE(lri_opt_type), POINTER :: lri_opt
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(LEN=*), PARAMETER :: routineN = 'get_original_gcc', &
routineP = moduleN//':'//routineN
INTEGER :: il, ipgf, iset, ishell, l, &
maxpgf, maxshell, nl, nset, &
stat
INTEGER, DIMENSION(:), POINTER :: lmax, lmin, ncont_l
LOGICAL :: failure
REAL(KIND=dp) :: expzet, gcca, prefac, zeta
failure = .FALSE.
maxpgf = SIZE(gto_basis_set%gcc,1)
maxshell = SIZE(gto_basis_set%gcc,2)
nset = SIZE(gto_basis_set%gcc,3)
ALLOCATE(gcc_orig(maxpgf,maxshell,nset),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
gcc_orig = 0.0_dp
DO iset=1,gto_basis_set%nset
DO ishell=1,gto_basis_set%nshell(iset)
l = gto_basis_set%l(ishell,iset)
expzet = 0.25_dp*REAL(2*l + 3,dp)
prefac = 2.0_dp**l*(2.0_dp/pi)**0.75_dp
DO ipgf=1,gto_basis_set%npgf(iset)
gcca = gto_basis_set%gcc(ipgf,ishell,iset)
zeta = gto_basis_set%zet(ipgf,iset)
gcc_orig(ipgf,ishell,iset) = gcca/(prefac*zeta**expzet)
END DO
END DO
END DO
IF(lri_opt%opt_coeffs) THEN
! **** get number of contractions per quantum number
CALL get_gto_basis_set(gto_basis_set=gto_basis_set,&
lmax=lmax,lmin=lmin)
ALLOCATE(lri_opt%subset(nset),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
DO iset=1,gto_basis_set%nset
nl=lmax(iset)-lmin(iset)+1
lri_opt%subset(iset)%nl=nl
il=1
ALLOCATE(lri_opt%subset(iset)%ncont_l(nl),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
ncont_l => lri_opt%subset(iset)%ncont_l
ncont_l=1
DO ishell=2,gto_basis_set%nshell(iset)
l = gto_basis_set%l(ishell,iset)
IF(l==gto_basis_set%l(ishell-1,iset)) THEN
ncont_l(il)=ncont_l(il) + 1
ELSE
il = il + 1
ncont_l(il)=1
ENDIF
ENDDO
END DO
ENDIF
END SUBROUTINE get_original_gcc
! *****************************************************************************
!> \brief orthonormalize contraction coefficients using Gram-Schmidt
!> \param gcc contraction coefficient
!> \param gto_basis_set gaussian type basis set
!> \param lri_opt optimization environment
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
! *****************************************************************************
SUBROUTINE orthonormalize_gcc(gcc,gto_basis_set,lri_opt,error)
REAL(KIND=dp), DIMENSION(:, :, :), &
POINTER :: gcc
TYPE(gto_basis_set_type), POINTER :: gto_basis_set
TYPE(lri_opt_type), POINTER :: lri_opt
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(LEN=*), PARAMETER :: routineN = 'orthonormalize_gcc', &
routineP = moduleN//':'//routineN
INTEGER :: il, iset, ishell, ishell1, &
ishell2, istart, nset
INTEGER, DIMENSION(:), POINTER :: nshell
LOGICAL :: failure
REAL(KIND=dp) :: gs_scale
failure = .FALSE.
CALL get_gto_basis_set(gto_basis_set=gto_basis_set,nset=nset,nshell=nshell)
DO iset=1,nset
istart = 1
DO il=1,lri_opt%subset(iset)%nl
DO ishell1=istart,istart+lri_opt%subset(iset)%ncont_l(il)-2
DO ishell2=ishell1+1,istart+lri_opt%subset(iset)%ncont_l(il)-1
gs_scale=DOT_PRODUCT(gcc(:,ishell2,iset),gcc(:,ishell1,iset))/&
DOT_PRODUCT(gcc(:,ishell1,iset),gcc(:,ishell1,iset))
gcc(:,ishell2,iset)=gcc(:,ishell2,iset)-&
gs_scale*gcc(:,ishell1,iset)
ENDDO
ENDDO
istart=istart+lri_opt%subset(iset)%ncont_l(il)
ENDDO
DO ishell=1,gto_basis_set%nshell(iset)
gcc(:,ishell,iset)=gcc(:,ishell,iset)/&
SQRT(DOT_PRODUCT(gcc(:,ishell,iset),gcc(:,ishell,iset)))
END DO
ENDDO
END SUBROUTINE orthonormalize_gcc
END MODULE lri_optimize_ri_basis_types

24
tests/QS/regtest-lrigpw/BASIS_LRI
View file

@ -264,3 +264,27 @@
0.235808199465 1.0 1.0 1.0
1 0 0 1 1
0.093521836600 1.0
O LRI_SZV-GTH_contract
6
2 0 0 7 3
24.031909411024 0.136381346028 -0.487008842022 0.371951484702
9.531086324462 -0.008551530264 -0.336891840071 0.464174122125
3.780041151565 0.120145433842 0.353938939468 -0.369361638096
1.499169204968 0.476023668732 -0.608592982400 -0.503256055757
0.594572443793 0.664249286619 0.192519561986 -0.051733664727
0.235808199466 0.318856064555 0.251367060276 0.358770747752
0.093521836600 0.444294088012 0.231125299456 0.353704648706
2 1 1 5 2
9.531086324462 0.293811752557 -0.487403705598
3.780041151565 0.473570777123 0.303086024199
1.499169204968 0.568069258815 -0.593766876464
0.594572443793 0.469346353287 0.373506422283
0.235808199466 0.382644340034 0.422504838590
2 2 2 1 1
3.780041151565 1.000000000000
2 2 2 1 1
1.499169204968 1.000000000000
2 2 2 1 1
0.594572443793 1.000000000000
2 3 3 1 1
1.499169204968 1.000000000000

67
tests/QS/regtest-lrigpw/O2_opt_lribas_contract.inp Normal file
View file

@ -0,0 +1,67 @@
&GLOBAL
PROJECT O2_opt_lribas_contract
RUN_TYPE ENERGY
PRINT_LEVEL low
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
CHARGE 0
MULTIPLICITY 3
UNRESTRICTED_KOHN_SHAM T
BASIS_SET_FILE_NAME BASIS_LRI
POTENTIAL_FILE_NAME ../../../data/POTENTIAL
&MGRID
NGRIDS 1
CUTOFF 150
&END MGRID
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END
&QS
METHOD GPW
&OPTIMIZE_LRI_BASIS
ACCURACY 1.0E-10
STEP_SIZE 0.005
MAX_FUN 4
DEGREES_OF_FREEDOM ALL
GEOMETRIC_SEQUENCE
USE_CONDITION_NUMBER
CONDITION_WEIGHT 1.0E-4
&CONSTRAIN_EXPONENTS
SCALE 0.3
&END
&END
&END QS
&SCF
SCF_GUESS atomic
MAX_SCF 2
&DIAGONALIZATION
&END
&END SCF
&XC
&XC_FUNCTIONAL PADE
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 12.0 12.0 12.0
PERIODIC NONE
&END CELL
&COORD
O 0.000000 0.000000 0.000000
O 0.000000 0.000000 1.210000
&END COORD
&TOPOLOGY
&CENTER_COORDINATES
&END
&END
&KIND O
BASIS_SET SZV-MOLOPT-GTH
POTENTIAL GTH-PADE-q6
LRI_BASIS_SET LRI_SZV-GTH_contract
&END KIND
&END SUBSYS
&END FORCE_EVAL

1
tests/QS/regtest-lrigpw/TEST_FILES
View file

@ -9,4 +9,5 @@ H2He_tz2p_lri.inp 1 6e-12
H2_tz2p_lri_diag.inp 1 8e-12
H2_tz2p_lri_ot.inp 1 4e-12
O2_opt_lribas.inp 46
O2_opt_lribas_contract.inp 46
O2_debug_ints.inp 0