deleting the several RC_* keywords in SE and introducing several section for

COULOMB, EXCHANGE and SCREENING. Added different taper functions for the
several terms and reset 2 regtests for numerics.


svn-origin-rev: 8322
This commit is contained in:
Teodoro Laino 2009-03-17 16:28:25 +00:00
parent d135a747ce
commit cf0dda02de
38 changed files with 840 additions and 389 deletions

View file

@ -85,9 +85,10 @@ MODULE cp_control_types
LOGICAL :: scp
INTEGER :: integral_screening, periodic_type
REAL(KIND = dp) :: delta
REAL(KIND = dp) :: rc_interaction
REAL(KIND = dp) :: rc_coulomb
REAL(KIND = dp) :: rc_range
! Parameters controlling the evaluation of the integrals
REAL(KIND = dp) :: cutoff_cou, taper_cou, range_cou
REAL(KIND = dp) :: cutoff_exc, taper_exc, range_exc
REAL(KIND = dp) :: taper_scr , range_scr
END TYPE semi_empirical_control_type
! *****************************************************************************

View file

@ -376,7 +376,7 @@ CPSourceFileRef,&
INTEGER :: itmp, j, jj, k, n_rep, n_var, &
nrep, stat
INTEGER, DIMENSION(:), POINTER :: tmplist
LOGICAL :: failure, was_present
LOGICAL :: explicit, failure, was_present
REAL(dp) :: tmpsqrt, value
TYPE(enumeration_type), POINTER :: enum
TYPE(keyword_type), POINTER :: keyword
@ -626,17 +626,43 @@ CPSourceFileRef,&
l_val=qs_control%se_control%analytical_gradients,error=error)
CALL section_vals_val_get(se_section,"FORCE_KDSO-D_EXCHANGE",&
l_val=qs_control%se_control%force_kdsod_EX,error=error)
CALL section_vals_val_get(se_section,"RC_INTERACTION",&
r_val=qs_control%se_control%rc_interaction,error=error)
CALL section_vals_val_get(se_section,"RC_COULOMB",&
r_val=qs_control%se_control%rc_coulomb,error=error)
CALL section_vals_val_get(se_section,"RC_RANGE",&
r_val=qs_control%se_control%rc_range,error=error)
CALL section_vals_val_get(se_section,"SCP",&
l_val=qs_control%se_control%scp,error=error)
! Integral Screening
CALL section_vals_val_get(se_section,"INTEGRAL_SCREENING",&
i_val=qs_control%se_control%integral_screening,error=error)
! Coulomb
CALL section_vals_val_get(se_section,"COULOMB%CUTOFF",&
r_val=qs_control%se_control%cutoff_cou,error=error)
qs_control%se_control%taper_cou = qs_control%se_control%cutoff_cou
CALL section_vals_val_get(se_section,"COULOMB%RC_TAPER",&
explicit=explicit, error=error)
IF (explicit) THEN
CALL section_vals_val_get(se_section,"COULOMB%RC_TAPER",&
r_val=qs_control%se_control%taper_cou, error=error)
END IF
CALL section_vals_val_get(se_section,"COULOMB%RC_RANGE",&
r_val=qs_control%se_control%range_cou,error=error)
! Exchange
CALL section_vals_val_get(se_section,"EXCHANGE%CUTOFF",&
r_val=qs_control%se_control%cutoff_exc,error=error)
qs_control%se_control%taper_exc = qs_control%se_control%cutoff_exc
CALL section_vals_val_get(se_section,"EXCHANGE%RC_TAPER",&
explicit=explicit, error=error)
IF (explicit) THEN
CALL section_vals_val_get(se_section,"EXCHANGE%RC_TAPER",&
r_val=qs_control%se_control%taper_exc, error=error)
END IF
CALL section_vals_val_get(se_section,"EXCHANGE%RC_RANGE",&
r_val=qs_control%se_control%range_exc,error=error)
! Screening (only if the integral scheme is of dumped type)
IF (qs_control%se_control%integral_screening==do_se_IS_kdso_d) THEN
CALL section_vals_val_get(se_section,"SCREENING%RC_TAPER",&
r_val=qs_control%se_control%taper_scr, error=error)
CALL section_vals_val_get(se_section,"SCREENING%RC_RANGE",&
r_val=qs_control%se_control%range_scr,error=error)
END IF
! SCP
CALL section_vals_val_get(se_section,"SCP",&
l_val=qs_control%se_control%scp,error=error)
! Periodic Type Calculation
CALL section_vals_val_get(se_section,"PERIODIC",&
i_val=qs_control%se_control%periodic_type,error=error)

View file

@ -2284,25 +2284,6 @@ CONTAINS
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="RC_INTERACTION",&
description="Atomic Cutoff Radius for Exchange and Core Hamiltonian terms.",&
usage="RC_INTERACTION {real} ",unit_str="bohr",default_r_val=24._dp, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="RC_COULOMB",&
description="Atomic Cutoff Radius for Coulomb interactions.",&
usage="RC_COULOMB {real} ",unit_str="bohr",default_r_val=24._dp, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="RC_RANGE",&
description="Range of cutoff switch function (tapering): 0.5*(1-TANH((r-r0)/RC_RANGE)), "//&
"where r0=2*RC_COULOMB-20*RC_RANGE. ",&
usage="RC_RANGE {real} ",unit_str="bohr",default_r_val=0.0_dp, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="INTEGRAL_SCREENING",&
description="Specifies the functional form for the ",&
usage="INTEGRAL_SCREENING (KDSO|KDSO-D|SLATER)",&
@ -2342,13 +2323,25 @@ CONTAINS
CALL keyword_create(keyword, name="FORCE_KDSO-D_EXCHANGE",&
description="This keywords forces the usage of the KDSO-D integral screening "//&
"for the Exchange integrals.",default_l_val=.FALSE.,lone_keyword_l_val=.TRUE.,&
"for the Exchange integrals (default is to apply the screening only to the "//&
"Coulomb integrals.",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)
NULLIFY(subsection)
CALL create_coulomb_section(subsection,error)
CALL section_add_subsection(section, subsection, error=error)
CALL section_release(subsection,error=error)
CALL create_exchange_section(subsection,error)
CALL section_add_subsection(section, subsection, error=error)
CALL section_release(subsection,error=error)
CALL create_screening_section(subsection,error)
CALL section_add_subsection(section, subsection, error=error)
CALL section_release(subsection,error=error)
CALL create_neighbor_lists_section(subsection,error)
CALL section_add_subsection(section, subsection, error=error)
CALL section_release(subsection,error=error)
@ -2364,6 +2357,149 @@ CONTAINS
END SUBROUTINE create_se_control_section
! *****************************************************************************
!> \brief Create the COULOMB se section
!> \param section the section to create
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
!> \author Teodoro Laino [tlaino]
!> \date 03.2009
! *****************************************************************************
SUBROUTINE create_coulomb_section(section,error)
TYPE(section_type), POINTER :: section
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'create_coulomb_section', &
routineP = moduleN//':'//routineN
LOGICAL :: failure
TYPE(keyword_type), POINTER :: keyword
failure=.FALSE.
CPPrecondition(.NOT.ASSOCIATED(section),cp_failure_level,routineP,error,failure)
IF (.NOT. failure) THEN
CALL section_create(section,name="COULOMB",&
description="Setup parameters for the evaluation of the COULOMB term in SE "//&
"calculations.", n_keywords=0, n_subsections=1, repeats=.FALSE., required=.FALSE.,&
error=error)
NULLIFY(keyword)
CALL keyword_create(keyword, name="CUTOFF",&
description="Atomic Cutoff Radius Cutoff for the evaluation of the Coulomb integrals. ",&
usage="CUTOFF {real} ",unit_str="bohr",default_r_val=24.0_dp, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="RC_TAPER",&
description="Atomic Cutoff Radius Cutoff for Tapering Coulomb integrals. "//&
"If not specified it assumes the same value specified for the CUTOFF.",&
usage="RC_TAPER {real} ",unit_str="bohr",type_of_var=real_t, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="RC_RANGE",&
description="Range of cutoff switch function (tapering): 0.5*(1-TANH((r-r0)/RC_RANGE)), "//&
"where r0=2*RC_TAPER-20*RC_RANGE.",&
usage="RC_RANGE {real} ",unit_str="bohr",default_r_val=0.0_dp, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
END IF
END SUBROUTINE create_coulomb_section
! *****************************************************************************
!> \brief Create the EXCHANGE se section
!> \param section the section to create
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
!> \author Teodoro Laino [tlaino]
!> \date 03.2009
! *****************************************************************************
SUBROUTINE create_exchange_section(section,error)
TYPE(section_type), POINTER :: section
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'create_exchange_section', &
routineP = moduleN//':'//routineN
LOGICAL :: failure
TYPE(keyword_type), POINTER :: keyword
failure=.FALSE.
CPPrecondition(.NOT.ASSOCIATED(section),cp_failure_level,routineP,error,failure)
IF (.NOT. failure) THEN
CALL section_create(section,name="EXCHANGE",&
description="Setup parameters for the evaluation of the EXCHANGE and "//&
" core Hamiltonian terms in SE calculations.", n_keywords=0, n_subsections=1,&
repeats=.FALSE., required=.FALSE.,error=error)
NULLIFY(keyword)
CALL keyword_create(keyword, name="CUTOFF",&
description="Atomic Cutoff Radius Cutoff for the evaluation of the Exchange integrals. ",&
usage="CUTOFF {real} ",unit_str="bohr",default_r_val=24.0_dp, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="RC_TAPER",&
description="Atomic Cutoff Radius Cutoff for Tapering Exchange integrals. "//&
"If not specified it assumes the same value specified for the CUTOFF.",&
usage="RC_TAPER {real} ",unit_str="bohr",type_of_var=real_t, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="RC_RANGE",&
description="Range of cutoff switch function (tapering): 0.5*(1-TANH((r-r0)/RC_RANGE)), "//&
"where r0=2*RC_TAPER-20*RC_RANGE.",&
usage="RC_RANGE {real} ",unit_str="bohr",default_r_val=0.0_dp, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
END IF
END SUBROUTINE create_exchange_section
! *****************************************************************************
!> \brief Create the SCREENING se section
!> \param section the section to create
!> \param error variable to control error logging, stopping,...
!> see module cp_error_handling
!> \author Teodoro Laino [tlaino]
!> \date 03.2009
! *****************************************************************************
SUBROUTINE create_screening_section(section,error)
TYPE(section_type), POINTER :: section
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'create_screening_section', &
routineP = moduleN//':'//routineN
LOGICAL :: failure
TYPE(keyword_type), POINTER :: keyword
failure=.FALSE.
CPPrecondition(.NOT.ASSOCIATED(section),cp_failure_level,routineP,error,failure)
IF (.NOT. failure) THEN
CALL section_create(section,name="SCREENING",&
description="Setup parameters for the tapering of the Coulomb/Exchange Screening in "//&
"KDSO-D integral scheme,", n_keywords=0, n_subsections=1, repeats=.FALSE., required=.FALSE.,&
error=error)
NULLIFY(keyword)
CALL keyword_create(keyword, name="RC_TAPER",&
description="Atomic Cutoff Radius Cutoff for Tapering the screening term. ",&
usage="RC_TAPER {real} ",unit_str="bohr",default_r_val=12.0_dp, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
CALL keyword_create(keyword, name="RC_RANGE",&
description="Range of cutoff switch function (tapering): 0.5*(1-TANH((r-r0)/RC_RANGE)), "//&
"where r0=2*RC_TAPER-20*RC_RANGE.",&
usage="RC_RANGE {real} ",unit_str="bohr",default_r_val=0.0_dp, error=error)
CALL section_add_keyword(section,keyword,error=error)
CALL keyword_release(keyword,error=error)
END IF
END SUBROUTINE create_screening_section
! *****************************************************************************
!> \brief Create the print se section
!> \param section the section to create

View file

@ -75,15 +75,17 @@ MODULE nddo_methods
semi_empirical_si_type
USE semi_empirical_types, ONLY: get_se_param,&
se_int_control_type,&
se_taper_type,&
semi_empirical_type,&
setup_se_int_control_type
USE semi_empirical_utils, ONLY: get_se_type
USE semi_empirical_utils, ONLY: finalize_se_taper,&
get_se_type,&
initialize_se_taper
USE sparse_matrix_types, ONLY: &
allocate_matrix_set, cp_sm_scale_and_add, cp_sm_sm_trace, &
create_dbm_from_sm, deallocate_dbm_from_sm, deallocate_matrix, &
get_block_node, mpsum_dbm, real_matrix_p_type, real_matrix_type, &
set_matrix, transfer_matrix
USE taper_types, ONLY: taper_type
USE timings, ONLY: timeset,&
timestop
#include "cp_common_uses.h"
@ -346,20 +348,21 @@ CONTAINS
TYPE(qs_force_type), DIMENSION(:), &
POINTER :: force
TYPE(se_int_control_type) :: se_int_control
TYPE(se_taper_type), POINTER :: se_taper
TYPE(semi_empirical_control_type), &
POINTER :: se_control
TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b
TYPE(taper_type), POINTER :: taper
failure = .FALSE.
enuclear = 0.0_dp
NULLIFY(dft_control,cell,force,particle_set,se_control,taper,atomic_kind_set)
NULLIFY(dft_control,cell,force,particle_set,se_control,se_taper,atomic_kind_set)
CALL timeset(routineN,handle)
CPPrecondition(ASSOCIATED(qs_env),cp_failure_level,routineP,error,failure)
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, taper=taper,&
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper,&
energy=energy, error=error)
CALL initialize_se_taper(se_taper,coulomb=.TRUE.,error=error)
! Parameters
se_control => dft_control%qs_control%se_control
anag = se_control%analytical_gradients
@ -421,7 +424,7 @@ CONTAINS
! Core-Core energy term
CALL corecore (se_kind_a,se_kind_b,rij,enuc=enuc,itype=itype,anag=anag,&
se_int_control=se_int_control, taper=taper, error=error)
se_int_control=se_int_control, se_taper=se_taper, error=error)
enuclear = enuclear + enuc
! Residual integral (1/R^3) correction
@ -440,7 +443,7 @@ CONTAINS
atom_b = atom_of_kind(jatom)
CALL dcorecore (se_kind_a,se_kind_b,rij,denuc=force_ab,itype=itype,delta=delta,&
anag=anag,se_int_control=se_int_control,taper=taper,error=error)
anag=anag,se_int_control=se_int_control,se_taper=se_taper,error=error)
! Residual integral (1/R^3) correction
IF (se_int_control%do_ewald_r3) THEN
@ -478,6 +481,7 @@ CONTAINS
energy%core_overlap = enuclear
END IF
CALL finalize_se_taper(se_taper,error=error)
CALL timestop(handle)
END SUBROUTINE se_core_core_interaction
@ -532,23 +536,24 @@ CONTAINS
TYPE(real_matrix_type), POINTER :: ks_scp, pscp
TYPE(scp_environment_type), POINTER :: scp_env
TYPE(se_int_control_type) :: se_int_control
TYPE(se_taper_type), POINTER :: se_taper
TYPE(semi_empirical_control_type), &
POINTER :: se_control
TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b
TYPE(taper_type), POINTER :: taper
failure=.FALSE.
CALL timeset(routineN,handle)
IF ( .NOT. failure ) THEN
NULLIFY(dft_control,cell,force,particle_set,diagmat_ks,diagmat_p,&
se_control,taper,pscp,ks_scp,scp_env)
se_control,se_taper,pscp,ks_scp,scp_env)
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, taper=taper,&
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper,&
para_env=para_env, sab_2c=sab_2c, atomic_kind_set=atomic_kind_set, &
particle_set=particle_set, error=error)
! Parameters
CALL initialize_se_taper(se_taper,coulomb=.TRUE.,error=error)
se_control => dft_control%qs_control%se_control
anag = se_control%analytical_gradients
scp_nddo = se_control%scp
@ -733,28 +738,28 @@ CONTAINS
IF ( nspins == 1 ) THEN
CALL fock2_1el (se_kind_a,se_kind_b,rij, ksa_block_a, ksb_block_a,&
pa_a, pb_a, ecore=ecore2, itype=itype, anag=anag, se_int_control=se_int_control,&
taper=taper, store_int_env=store_int_env, error=error)
se_taper=se_taper, store_int_env=store_int_env, error=error)
ELSE IF ( nspins == 2 ) THEN
CALL fock2_1el (se_kind_a,se_kind_b,rij, ksa_block_a, ksb_block_a,&
pa_block_a, pb_block_a, ecore=ecore2, itype=itype, anag=anag,&
se_int_control=se_int_control, taper=taper, store_int_env=store_int_env,&
se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env,&
error=error)
CALL fock2_1el (se_kind_a,se_kind_b,rij, ksa_block_b, ksb_block_b,&
pa_b, pb_b, ecore=ecore2, itype=itype, anag=anag, se_int_control=se_int_control,&
taper=taper, store_int_env=store_int_env, error=error)
se_taper=se_taper, store_int_env=store_int_env, error=error)
END IF
! Coulomb Terms
IF ( nspins == 1 ) THEN
CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b,&
ksb_block_a, factor=0.5_dp, anag=anag, se_int_control=se_int_control, taper=taper,&
ksb_block_a, factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,&
store_int_env=store_int_env, error=error)
ELSE IF ( nspins == 2 ) THEN
CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b,&
ksb_block_a, factor=1.0_dp, anag=anag, se_int_control=se_int_control, taper=taper,&
ksb_block_a, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,&
store_int_env=store_int_env, error=error)
CALL fock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_b, p_block_tot_b,&
ksb_block_b, factor=1.0_dp, anag=anag, se_int_control=se_int_control, taper=taper,&
ksb_block_b, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper,&
store_int_env=store_int_env, error=error)
END IF
@ -766,14 +771,14 @@ CONTAINS
force_ab = 0.0_dp
IF ( nspins == 1 ) THEN
CALL dfock2_1el (se_kind_a,se_kind_b,rij, pa_a, pb_a, itype=itype, anag=anag,&
se_int_control=se_int_control, taper=taper, force=force_ab,&
se_int_control=se_int_control, se_taper=se_taper, force=force_ab,&
delta=delta, error=error)
ELSE IF ( nspins == 2 ) THEN
CALL dfock2_1el (se_kind_a,se_kind_b,rij, pa_block_a, pb_block_a, itype=itype, anag=anag,&
se_int_control=se_int_control, taper=taper, force=force_ab, delta=delta,&
se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,&
error=error)
CALL dfock2_1el (se_kind_a,se_kind_b,rij, pa_b, pb_b, itype=itype, anag=anag,&
se_int_control=se_int_control, taper=taper, force=force_ab, delta=delta,&
se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,&
error=error)
END IF
@ -791,15 +796,15 @@ CONTAINS
force_ab = 0._dp
IF ( nspins == 1 ) THEN
CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp,&
anag=anag, se_int_control=se_int_control, taper=taper, force=force_ab, delta=delta,&
anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,&
error=error)
ELSE IF ( nspins == 2 ) THEN
CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp,&
anag=anag, se_int_control=se_int_control, taper=taper, force=force_ab, delta=delta,&
anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,&
error=error)
CALL dfock2C(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp,&
anag=anag, se_int_control=se_int_control, taper=taper, force=force_ab, delta=delta,&
anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta,&
error=error)
END IF
IF ( switch ) THEN
@ -853,7 +858,7 @@ CONTAINS
CALL mp_sum(ecore2,para_env%group)
energy%hartree = ecore2 - energy%core
END IF
CALL finalize_se_taper(se_taper,error=error)
CALL timestop(handle)
END SUBROUTINE build_fock_matrix_coulomb
@ -1592,22 +1597,22 @@ CONTAINS
TYPE(qs_force_type), DIMENSION(:), &
POINTER :: force
TYPE(se_int_control_type) :: se_int_control
TYPE(se_taper_type), POINTER :: se_taper
TYPE(semi_empirical_control_type), &
POINTER :: se_control
TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b
TYPE(taper_type), POINTER :: taper
failure=.FALSE.
CALL timeset(routineN,handle)
NULLIFY(dft_control,cell,force,particle_set,se_control,taper)
CALL get_qs_env(qs_env=qs_env,dft_control=dft_control,cell=cell,taper=taper,&
NULLIFY(dft_control,cell,force,particle_set,se_control,se_taper)
CALL get_qs_env(qs_env=qs_env,dft_control=dft_control,cell=cell,se_taper=se_taper,&
para_env=para_env,error=error)
! set values for tapering function
CALL initialize_se_taper(se_taper,exchange=.TRUE.,error=error)
se_control => dft_control%qs_control%se_control
anag = se_control%analytical_gradients
nspins=dft_control%nspins
CPPrecondition(nspins<3,cp_failure_level,routineP,error,failure)
CPPrecondition(ASSOCIATED(matrix_p),cp_failure_level,routineP,error,failure)
CPPrecondition(SIZE(ks_matrix)>0,cp_failure_level,routineP,error,failure)
@ -1717,16 +1722,16 @@ CONTAINS
! Exchange Terms
IF ( nspins == 1 ) THEN
CALL fock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_tot, p_block_a, ks_block_a,&
factor=0.5_dp, anag=anag, se_int_control=se_int_control, taper=taper, store_int_env=store_int_env,&
error=error)
factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
store_int_env=store_int_env, error=error)
ELSE IF ( nspins == 2 ) THEN
CALL fock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_tot, p_block_a, ks_block_a,&
factor=1.0_dp, anag=anag, se_int_control=se_int_control, taper=taper, store_int_env=store_int_env,&
error=error)
factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
store_int_env=store_int_env, error=error)
CALL fock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_tot, p_block_b, ks_block_b,&
factor=1.0_dp, anag=anag, se_int_control=se_int_control, taper=taper, store_int_env=store_int_env,&
error=error)
factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
store_int_env=store_int_env, error=error)
END IF
IF(calculate_forces) THEN
atom_a = atom_of_kind(iatom)
@ -1734,15 +1739,15 @@ CONTAINS
force_ab = 0.0_dp
IF ( nspins == 1 ) THEN
CALL dfock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_tot, p_block_a,&
factor=0.5_dp, anag=anag, se_int_control=se_int_control, taper=taper, force=force_ab,&
factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab,&
delta=delta, error=error)
ELSE IF ( nspins == 2 ) THEN
CALL dfock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_tot, p_block_a,&
factor=1.0_dp, anag=anag, se_int_control=se_int_control, taper=taper, force=force_ab,&
factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab,&
delta=delta, error=error)
CALL dfock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_tot, p_block_b,&
factor=1.0_dp, anag=anag, se_int_control=se_int_control, taper=taper, force=force_ab,&
factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab,&
delta=delta, error=error)
END IF
IF (switch) THEN
@ -1771,6 +1776,7 @@ CONTAINS
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
END IF
CALL finalize_se_taper(se_taper,error=error)
CALL timestop(handle)
END SUBROUTINE build_fock_matrix_exchange
@ -1781,7 +1787,7 @@ CONTAINS
!> \date 04.2008 [tlaino]
! *****************************************************************************
SUBROUTINE fock2_1el (sepi, sepj, rij, ksi_block, ksj_block, pi_block, pj_block,&
ecore, itype, anag, se_int_control, taper, store_int_env, error)
ecore, itype, anag, se_int_control, se_taper, store_int_env, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
REAL(KIND=dp), &
@ -1800,7 +1806,7 @@ CONTAINS
INTEGER, INTENT(IN) :: itype
LOGICAL, INTENT(IN) :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(semi_empirical_si_type), POINTER :: store_int_env
TYPE(cp_error_type), INTENT(inout) :: error
@ -1814,7 +1820,7 @@ CONTAINS
failure = .FALSE.
! Compute integrals
CALL rotnuc (sepi, sepj, rij, e1b=e1b, e2a=e2a, itype=itype, anag=anag,&
se_int_control=se_int_control, taper=taper, store_int_env=store_int_env,&
se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env,&
error=error)
!
! Add the electron-nuclear attraction term for atom sepi
@ -1861,7 +1867,7 @@ CONTAINS
!> \date 04.2008 [tlaino]
! *****************************************************************************
SUBROUTINE dfock2_1el (sepi, sepj, rij, pi_block, pj_block, itype, anag,&
se_int_control, taper, force, delta, error)
se_int_control, se_taper, force, delta, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
REAL(KIND=dp), &
@ -1873,7 +1879,7 @@ CONTAINS
INTEGER, INTENT(IN) :: itype
LOGICAL, INTENT(IN) :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
REAL(KIND=dp), DIMENSION(3), &
INTENT(INOUT) :: force
REAL(KIND=dp), INTENT(IN) :: delta
@ -1890,7 +1896,7 @@ CONTAINS
failure = .FALSE.
! Compute integrals
CALL drotnuc (sepi, sepj, rij, de1b=de1b, de2a=de2a, itype=itype, anag=anag,&
se_int_control=se_int_control, taper=taper, delta=delta, error=error)
se_int_control=se_int_control, se_taper=se_taper, delta=delta, error=error)
!
! Add the electron-nuclear attraction term for atom sepi
!
@ -2009,7 +2015,7 @@ CONTAINS
!> \date 04.2008 [tlaino]
! *****************************************************************************
SUBROUTINE fock2C(sepi, sepj, rij, switch, pi_tot, fi_mat, pj_tot, fj_mat, &
factor, anag, se_int_control, taper, store_int_env, error)
factor, anag, se_int_control, se_taper, store_int_env, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
LOGICAL, INTENT(IN) :: switch
@ -2026,7 +2032,7 @@ CONTAINS
REAL(KIND=dp), INTENT(IN) :: factor
LOGICAL, INTENT(IN) :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(semi_empirical_si_type), POINTER :: store_int_env
TYPE(cp_error_type), INTENT(inout) :: error
@ -2042,12 +2048,12 @@ CONTAINS
failure = .FALSE.
! Evaluate integrals
IF (.NOT.switch) THEN
CALL rotint (sepi,sepj, rij,w,anag=anag,se_int_control=se_int_control,taper=taper,&
store_int_env=store_int_env, error=error)
CALL rotint (sepi,sepj, rij,w,anag=anag,se_int_control=se_int_control,&
se_taper=se_taper,store_int_env=store_int_env, error=error)
ELSE
irij = -rij
CALL rotint (sepj,sepi,irij,w,anag=anag,se_int_control=se_int_control,taper=taper,&
store_int_env=store_int_env, error=error)
CALL rotint (sepj,sepi,irij,w,anag=anag,se_int_control=se_int_control,&
se_taper=se_taper,store_int_env=store_int_env, error=error)
END IF
kr = 0
natorb(1) = sepi%natorb
@ -2099,7 +2105,7 @@ CONTAINS
!> \date 04.2008 [tlaino]
! *****************************************************************************
SUBROUTINE dfock2C(sepi, sepj, rij, switch, pi_tot, pj_tot, factor, anag,&
se_int_control, taper, force, delta, error)
se_int_control, se_taper, force, delta, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
LOGICAL, INTENT(IN) :: switch
@ -2108,7 +2114,7 @@ CONTAINS
REAL(KIND=dp), INTENT(IN) :: factor
LOGICAL, INTENT(IN) :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
REAL(KIND=dp), DIMENSION(3), &
INTENT(INOUT) :: force
REAL(KIND=dp), INTENT(IN) :: delta
@ -2128,11 +2134,11 @@ CONTAINS
! Evaluate integrals' derivatives
IF (.NOT.switch) THEN
CALL drotint (sepi,sepj, rij,dw,delta,anag=anag,se_int_control=se_int_control,&
taper=taper,error=error)
se_taper=se_taper,error=error)
ELSE
irij = -rij
CALL drotint (sepj,sepi,irij,dw,delta,anag=anag,se_int_control=se_int_control,&
taper=taper,error=error)
se_taper=se_taper,error=error)
END IF
kr = 0
@ -2183,7 +2189,7 @@ CONTAINS
!> \date 04.2008 [tlaino]
! *****************************************************************************
SUBROUTINE fock2E(sepi, sepj, rij, switch, isize, pi_tot, pi_mat, fi_mat, factor,&
anag, se_int_control, taper, store_int_env, error)
anag, se_int_control, se_taper, store_int_env, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
LOGICAL, INTENT(IN) :: switch
@ -2199,7 +2205,7 @@ CONTAINS
REAL(KIND=dp), INTENT(IN) :: factor
LOGICAL, INTENT(IN) :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(semi_empirical_si_type), POINTER :: store_int_env
TYPE(cp_error_type), INTENT(inout) :: error
@ -2215,12 +2221,12 @@ CONTAINS
failure = .FALSE.
! Evaluate integrals
IF (.NOT.switch) THEN
CALL rotint (sepi,sepj, rij,w,anag=anag,se_int_control=se_int_control,taper=taper,&
store_int_env=store_int_env, error=error)
CALL rotint (sepi,sepj, rij,w,anag=anag,se_int_control=se_int_control,&
se_taper=se_taper,store_int_env=store_int_env, error=error)
ELSE
irij = -rij
CALL rotint (sepj,sepi,irij,w,anag=anag,se_int_control=se_int_control,taper=taper,&
store_int_env=store_int_env, error=error)
CALL rotint (sepj,sepi,irij,w,anag=anag,se_int_control=se_int_control,&
se_taper=se_taper,store_int_env=store_int_env, error=error)
END IF
kr = 0
natorb(1) = sepi%natorb
@ -2266,7 +2272,7 @@ CONTAINS
!> \date 04.2008 [tlaino]
! *****************************************************************************
SUBROUTINE dfock2E(sepi, sepj, rij, switch, isize, pi_tot, pi_mat, factor, anag,&
se_int_control, taper, force, delta, error)
se_int_control, se_taper, force, delta, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rij
LOGICAL, INTENT(IN) :: switch
@ -2279,7 +2285,7 @@ CONTAINS
REAL(KIND=dp), INTENT(IN) :: factor
LOGICAL, INTENT(IN) :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
REAL(KIND=dp), DIMENSION(3), &
INTENT(INOUT) :: force
REAL(KIND=dp), INTENT(IN) :: delta
@ -2300,11 +2306,11 @@ CONTAINS
! Evaluate integrals' derivatives
IF (.NOT.switch) THEN
CALL drotint (sepi,sepj, rij,dw,delta,anag=anag,se_int_control=se_int_control,&
taper=taper,error=error)
se_taper=se_taper,error=error)
ELSE
irij = -rij
CALL drotint (sepj,sepi,irij,dw,delta,anag=anag,se_int_control=se_int_control,&
taper=taper,error=error)
se_taper=se_taper,error=error)
END IF
kr = 0

View file

@ -43,6 +43,7 @@ MODULE qmmm_se_energy
rotnuc
USE semi_empirical_types, ONLY: get_se_param,&
se_int_control_type,&
se_taper_type,&
semi_empirical_create,&
semi_empirical_release,&
semi_empirical_type,&
@ -55,7 +56,6 @@ MODULE qmmm_se_energy
real_matrix_p_type,&
replicate_matrix_structure,&
set_matrix
USE taper_types, ONLY: taper_type
USE termination, ONLY: stop_program
USE timings, ONLY: timeset,&
timestop
@ -104,15 +104,15 @@ CONTAINS
TYPE(real_matrix_p_type), DIMENSION(:), &
POINTER :: matrix_s
TYPE(se_int_control_type) :: se_int_control
TYPE(se_taper_type), POINTER :: se_taper
TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm
TYPE(taper_type), POINTER :: taper
CALL timeset(routineN,handle)
NULLIFY(logger)
logger => cp_error_get_logger(error)
NULLIFY (matrix_s, atomic_kind_set)
NULLIFY (se_kind_a, se_kind_mm, taper)
NULLIFY (se_kind_a, se_kind_mm, se_taper)
CALL build_qs_neighbor_lists(qs_env,para_env,force_env_section=qs_env%input,error=error)
CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, error=error)
CALL build_overlap_matrix(qs_env,para_env,matrix_s=matrix_s,&
@ -123,7 +123,7 @@ CONTAINS
error=error)
CALL set_qs_env(qs_env=qs_env,matrix_s=matrix_s,error=error)
CALL get_qs_env(qs_env=qs_env,&
taper=taper,&
se_taper=se_taper,&
atomic_kind_set=atomic_kind_set,&
ks_qmmm_env=ks_qmmm_env_loc,&
dft_control=dft_control,error=error)
@ -190,7 +190,7 @@ CONTAINS
qmmm_env%qm_atom_index(iatom),&
enuclear,&
itype,&
taper,&
se_taper,&
se_int_control,&
anag,&
qmmm_env%spherical_cutoff,&
@ -210,7 +210,7 @@ CONTAINS
qmmm_env%qm_atom_index(iatom),&
enuclear,&
itype,&
taper,&
se_taper,&
se_int_control,&
anag,&
qmmm_env%spherical_cutoff,&
@ -245,7 +245,7 @@ CONTAINS
! *****************************************************************************
SUBROUTINE build_se_qmmm_matrix_low(h_block_a, se_kind_a, se_kind_mm, potentials,&
mm_particles, mm_charges, mm_el_pot_radius, mm_atom_index, num_mm_atoms,&
mm_cell, IndQM, enuclear, itype, taper, se_int_control, anag, &
mm_cell, IndQM, enuclear, itype, se_taper, se_int_control, anag, &
qmmm_spherical_cutoff, error)
REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block_a
@ -261,7 +261,7 @@ CONTAINS
INTEGER, INTENT(IN) :: IndQM
REAL(KIND=dp), INTENT(INOUT) :: enuclear
INTEGER, INTENT(IN) :: itype
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
LOGICAL, INTENT(IN) :: anag
REAL(KIND=dp), INTENT(IN) :: qmmm_spherical_cutoff(2)
@ -300,9 +300,9 @@ CONTAINS
END IF
IF (ABS(se_kind_mm%zeff)<=EPSILON(0.0_dp)) CYCLE
CALL rotnuc (se_kind_a, se_kind_mm, rij, itype=itype, e1b=e1b, anag=anag,&
se_int_control=se_int_control, taper=taper, error=error)
se_int_control=se_int_control, se_taper=se_taper, error=error)
CALL corecore(se_kind_a, se_kind_mm, rij, itype=itype, enuc=enuc, anag=anag,&
se_int_control=se_int_control, taper=taper, error=error)
se_int_control=se_int_control, se_taper=se_taper, error=error)
enuclear = enuclear + enuc
! Contribution to the iatom block
! Computation of the QMMM core matrix

View file

@ -36,6 +36,7 @@ MODULE qmmm_se_forces
drotnuc
USE semi_empirical_types, ONLY: get_se_param,&
se_int_control_type,&
se_taper_type,&
semi_empirical_create,&
semi_empirical_release,&
semi_empirical_type,&
@ -44,7 +45,6 @@ MODULE qmmm_se_forces
se_param_set_default
USE sparse_matrix_types, ONLY: get_block_node,&
real_matrix_p_type
USE taper_types, ONLY: taper_type
USE termination, ONLY: stop_program
USE timings, ONLY: timeset,&
timestop
@ -98,17 +98,17 @@ CONTAINS
TYPE(real_matrix_p_type), DIMENSION(:), &
POINTER :: matrix_p
TYPE(se_int_control_type) :: se_int_control
TYPE(se_taper_type), POINTER :: se_taper
TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_mm
TYPE(taper_type), POINTER :: taper
failure = .FALSE.
CALL timeset(routineN,handle)
IF (calc_force) THEN
NULLIFY (rho, atomic_kind_set, taper)
NULLIFY (rho, atomic_kind_set, se_taper)
NULLIFY (se_kind_a, se_kind_mm)
CALL get_qs_env(qs_env=qs_env,&
rho=rho,&
taper=taper,&
se_taper=se_taper,&
atomic_kind_set=atomic_kind_set,&
ks_qmmm_env=ks_qmmm_env_loc,&
dft_control=dft_control,error=error)
@ -176,7 +176,7 @@ CONTAINS
para_env,&
Forces,&
Forces_QM(:,iqm),&
taper,&
se_taper,&
se_int_control,&
anag,&
delta,&
@ -199,7 +199,7 @@ CONTAINS
para_env,&
Forces_added_charges,&
Forces_QM(:,iqm),&
taper,&
se_taper,&
se_int_control,&
anag,&
delta,&
@ -245,7 +245,7 @@ CONTAINS
! *****************************************************************************
SUBROUTINE deriv_se_qmmm_matrix_low(p_block_a, se_kind_a, se_kind_mm,&
potentials, mm_particles, mm_charges, mm_el_pot_radius, mm_atom_index, &
num_mm_atoms, mm_cell, IndQM, itype, para_env, forces, forces_qm, taper,&
num_mm_atoms, mm_cell, IndQM, itype, para_env, forces, forces_qm, se_taper,&
se_int_control, anag, delta, qmmm_spherical_cutoff, error)
REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block_a
@ -263,7 +263,7 @@ CONTAINS
REAL(KIND=dp), DIMENSION(:, :), POINTER :: forces
REAL(KIND=dp), DIMENSION(:), &
INTENT(INOUT) :: forces_qm
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
LOGICAL, INTENT(IN) :: anag
REAL(KIND=dp), INTENT(IN) :: delta, &
@ -304,9 +304,11 @@ CONTAINS
IF (ABS(se_kind_mm%zeff)<=EPSILON(0.0_dp)) CYCLE
! Integrals derivatives involving QM - MM atoms
CALL drotnuc(se_kind_a, se_kind_mm, rij, itype=itype, de1b=de1b,&
se_int_control=se_int_control, anag=anag, delta=delta, taper=taper, error=error)
se_int_control=se_int_control, anag=anag, delta=delta, &
se_taper=se_taper, error=error)
CALL dcorecore(se_kind_a, se_kind_mm, rij, itype=itype, denuc=denuc,&
se_int_control=se_int_control, anag=anag, delta=delta, taper=taper, error=error)
se_int_control=se_int_control, anag=anag, delta=delta, &
se_taper=se_taper, error=error)
! Nucler - Nuclear term
force_ab(1:3)=-denuc(1:3)
! Force contribution from the QMMM Hamiltonian

View file

@ -164,8 +164,8 @@ MODULE qs_environment
USE semi_empirical_mpole_types, ONLY: nddo_mpole_type
USE semi_empirical_store_int_types, ONLY: semi_empirical_si_create,&
semi_empirical_si_type
USE taper_types, ONLY: taper_create,&
taper_type
USE semi_empirical_types, ONLY: se_taper_create,&
se_taper_type
USE termination, ONLY: stop_memory,&
stop_program
USE timings, ONLY: timeset,&
@ -403,13 +403,13 @@ CONTAINS
TYPE(rel_control_type), POINTER :: rel_control
TYPE(rho0_mpole_type), POINTER :: rho0_mpole
TYPE(scf_control_type), POINTER :: scf_control
TYPE(se_taper_type), POINTER :: se_taper
TYPE(section_vals_type), POINTER :: dft_section, et_becke_section, &
et_coupling_section, et_ddapc_section, ewald_section, poisson_section, &
print_section, qs_section, se_section, vdw_section, xc_fun_section
TYPE(semi_empirical_control_type), &
POINTER :: se_control
TYPE(semi_empirical_si_type), POINTER :: se_store_int_env
TYPE(taper_type), POINTER :: taper
CALL timeset(routineN,handle)
NULLIFY(logger)
@ -428,7 +428,7 @@ CONTAINS
END IF
! Initialise the Quickstep environment
NULLIFY (mos, taper, mos_aux_fit)
NULLIFY (mos, se_taper, mos_aux_fit)
NULLIFY (dft_control)
NULLIFY (energy)
NULLIFY (force)
@ -681,22 +681,22 @@ CONTAINS
se_section => section_vals_get_subs_vals(qs_section,"SE",error=error)
se_control => dft_control%qs_control%se_control
! Make interaction radius compatible with unit cell
rc = se_control%rc_interaction
rc = se_control%cutoff_exc
IF(cell%perd(1)==1) rc = MIN(rc,0.25_dp*plane_distance(1,0,0,cell))
IF(cell%perd(2)==1) rc = MIN(rc,0.25_dp*plane_distance(0,1,0,cell))
IF(cell%perd(3)==1) rc = MIN(rc,0.25_dp*plane_distance(0,0,1,cell))
IF (rc /= se_control%rc_interaction) THEN
IF (rc /= se_control%cutoff_exc) THEN
IF (output_unit>0) THEN
WRITE(output_unit,*)
WRITE(output_unit,'(A,T34,A)') " WARNING|",&
" Integral interaction cutoff has been redefined"
WRITE(output_unit,'(A,T71,F10.3)') " WARNING| Old value [a.u.]",&
se_control%rc_interaction
se_control%cutoff_exc
WRITE(output_unit,'(A,T71,F10.3)') " WARNING| New value [a.u.]",rc
WRITE(output_unit,*)
END IF
END IF
se_control%rc_interaction = rc
se_control%cutoff_exc = rc
SELECT CASE ( dft_control%qs_control%method_id)
CASE DEFAULT
CASE (do_method_rm1,do_method_am1,do_method_mndo,do_method_pm3,&
@ -737,8 +737,11 @@ CONTAINS
dft_control%qs_control%method_id,error)
END IF
! Taper function
CALL taper_create(taper, se_control%rc_coulomb, se_control%rc_range, error)
CALL set_qs_env(qs_env, taper=taper, error=error)
CALL se_taper_create(se_taper, se_control%integral_screening, &
se_control%taper_cou, se_control%range_cou, &
se_control%taper_exc, se_control%range_exc, &
se_control%taper_scr, se_control%range_scr, error)
CALL set_qs_env(qs_env, se_taper=se_taper, error=error)
! Store integral environment
CALL semi_empirical_si_create(se_store_int_env, se_section, error)
CALL set_qs_env(qs_env, se_store_int_env=se_store_int_env, error=error)

View file

@ -159,10 +159,10 @@ MODULE qs_environment_types
nddo_mpole_type
USE semi_empirical_store_int_types, ONLY: semi_empirical_si_release,&
semi_empirical_si_type
USE semi_empirical_types, ONLY: se_taper_release,&
se_taper_type
USE sparse_matrix_types, ONLY: deallocate_matrix_set,&
real_matrix_p_type
USE taper_types, ONLY: taper_release,&
taper_type
USE task_list_types, ONLY: deallocate_task_list,&
task_list_type
USE timings, ONLY: timeset,&
@ -372,7 +372,7 @@ MODULE qs_environment_types
TYPE(ewald_environment_type),POINTER :: ewald_env
TYPE(ewald_pw_type),POINTER :: ewald_pw
! Semi-empirical types
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(semi_empirical_si_type),POINTER :: se_store_int_env
TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole
TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env
@ -411,7 +411,8 @@ CONTAINS
VH_3c_list,rho0_s_rs,rho0_s_gs,requires_mo_derivs,mo_derivs,mo_derivs_aux_fit,tddfpt_control,&
neighbor_list_id,kg_sub_pw_env,linres_control,xas_env,virial,cp_ddapc_env,cp_ddapc_ewald,&
outer_scf_history,outer_scf_ihistory,ep_qs_env,x_data,et_coupling,dftb_potential,results,&
scp_env,taper,se_store_int_env,se_nddo_mpole,se_nonbond_env, wfn_fitting_env, dispersion_env, error)
scp_env,se_taper,se_store_int_env,se_nddo_mpole,se_nonbond_env, wfn_fitting_env, &
dispersion_env, error)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(atomic_kind_type), DIMENSION(:), &
OPTIONAL, POINTER :: atomic_kind_set
@ -527,7 +528,7 @@ CONTAINS
TYPE(cp_result_type), OPTIONAL, POINTER :: results
TYPE(scp_environment_type), OPTIONAL, &
POINTER :: scp_env
TYPE(taper_type), OPTIONAL, POINTER :: taper
TYPE(se_taper_type), OPTIONAL, POINTER :: se_taper
TYPE(semi_empirical_si_type), OPTIONAL, &
POINTER :: se_store_int_env
TYPE(nddo_mpole_type), OPTIONAL, POINTER :: se_nddo_mpole
@ -623,7 +624,7 @@ CONTAINS
IF (PRESENT(tddfpt_control)) tddfpt_control => qs_env%dft_control%tddfpt_control
IF (PRESENT(linres_control)) linres_control => qs_env%linres_control
IF (PRESENT(virial)) virial => qs_env%virial
IF (PRESENT(taper)) taper => qs_env%taper
IF (PRESENT(se_taper)) se_taper => qs_env%se_taper
IF (PRESENT(se_store_int_env)) se_store_int_env => qs_env%se_store_int_env
IF (PRESENT(se_nddo_mpole)) se_nddo_mpole => qs_env%se_nddo_mpole
IF (PRESENT(se_nonbond_env)) se_nonbond_env => qs_env%se_nonbond_env
@ -774,7 +775,7 @@ CONTAINS
NULLIFY (qs_env%task_list_aux_fit)
NULLIFY (qs_env%task_list_soft)
NULLIFY (qs_env%scp_env)
NULLIFY (qs_env%taper)
NULLIFY (qs_env%se_taper)
NULLIFY (qs_env%se_store_int_env)
NULLIFY (qs_env%se_nddo_mpole)
NULLIFY (qs_env%se_nonbond_env)
@ -824,7 +825,8 @@ CONTAINS
rhoz_set,rhoz_tot,ecoul_1c,VH_3c_list,requires_mo_derivs,mo_derivs,mo_derivs_aux_fit,&
neighbor_list_id,kg_sub_pw_env,linres_control,xas_env,virial,cp_ddapc_env,cp_ddapc_ewald,&
outer_scf_history,outer_scf_ihistory,ep_qs_env,x_data,et_coupling,dftb_potential,&
scp_env,taper,se_store_int_env,se_nddo_mpole,se_nonbond_env,wfn_fitting_env,dispersion_env,error)
scp_env,se_taper,se_store_int_env,se_nddo_mpole,se_nonbond_env,wfn_fitting_env,&
dispersion_env,error)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(cell_type), OPTIONAL, POINTER :: cell, super_cell, cell_ref
@ -916,7 +918,7 @@ CONTAINS
DIMENSION(:, :), OPTIONAL, POINTER :: dftb_potential
TYPE(scp_environment_type), OPTIONAL, &
POINTER :: scp_env
TYPE(taper_type), OPTIONAL, POINTER :: taper
TYPE(se_taper_type), OPTIONAL, POINTER :: se_taper
TYPE(semi_empirical_si_type), OPTIONAL, &
POINTER :: se_store_int_env
TYPE(nddo_mpole_type), OPTIONAL, POINTER :: se_nddo_mpole
@ -1163,9 +1165,9 @@ CONTAINS
CALL scp_env_release(qs_env%scp_env,error=error)
qs_env%scp_env => scp_env
END IF
IF (PRESENT(taper)) THEN
CALL taper_release(qs_env%taper,error=error)
qs_env%taper => taper
IF (PRESENT(se_taper)) THEN
CALL se_taper_release(qs_env%se_taper,error=error)
qs_env%se_taper => se_taper
END IF
IF (PRESENT(se_store_int_env)) THEN
CALL semi_empirical_si_release(qs_env%se_store_int_env,error=error)
@ -1533,8 +1535,8 @@ CONTAINS
IF (ASSOCIATED(qs_env%scp_env)) THEN
CALL scp_env_release ( qs_env % scp_env, error )
END IF
IF (ASSOCIATED(qs_env%taper)) THEN
CALL taper_release(qs_env%taper, error)
IF (ASSOCIATED(qs_env%se_taper)) THEN
CALL se_taper_release(qs_env%se_taper, error)
END IF
IF (ASSOCIATED(qs_env%se_store_int_env)) THEN
CALL semi_empirical_si_release(qs_env%se_store_int_env, error)

View file

@ -531,7 +531,7 @@ CONTAINS
CALL get_gto_basis_set(gto_basis_set=orb_basis_set,&
kind_radius=kind_radius)
kind_radius = MAX (kind_radius,se_control%rc_interaction)
kind_radius = MAX (kind_radius,se_control%cutoff_exc)
CALL set_gto_basis_set(gto_basis_set=orb_basis_set,&
kind_radius=kind_radius)

View file

@ -774,7 +774,7 @@ CONTAINS
ALLOCATE (c_radius(nkind),STAT=stat)
CPPostcondition(stat==0,cp_failure_level,routineP,error,failure)
semi_present = .TRUE.
c_radius = qs_env%dft_control%qs_control%se_control%rc_coulomb
c_radius = qs_env%dft_control%qs_control%se_control%cutoff_cou
! Build the neighbor lists for the Hartree terms
CALL build_neighbor_lists(sab_2c,atom,cell,r_skin,maxatom_local,&

View file

@ -21,7 +21,8 @@ MODULE semi_empirical_int_ana
do_method_pdg,&
do_method_pm3,&
do_method_pm6,&
do_method_undef
do_method_undef,&
do_se_IS_kdso_d
USE kinds, ONLY: dp
USE physcon, ONLY: angstrom,&
evolt
@ -42,11 +43,11 @@ MODULE semi_empirical_int_ana
rotmat_type,&
se_int_control_type,&
se_int_screen_type,&
se_taper_type,&
semi_empirical_type,&
setup_se_int_control_type
USE taper_types, ONLY: dtaper_eval,&
taper_eval,&
taper_type
taper_eval
#include "cp_common_uses.h"
IMPLICIT NONE
@ -86,7 +87,7 @@ CONTAINS
!> \author Teodoro Laino [tlaino] - Zurich University
! *****************************************************************************
RECURSIVE SUBROUTINE rotnuc_ana (sepi,sepj,rijv,itype,e1b,e2a,de1b,de2a,&
se_int_control,taper,error)
se_int_control,se_taper,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
INTEGER, INTENT(IN) :: itype
@ -95,7 +96,7 @@ CONTAINS
REAL(dp), DIMENSION(3, 45), &
INTENT(OUT), OPTIONAL :: de1b, de2a
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'rotnuc_ana', &
@ -141,7 +142,7 @@ CONTAINS
END IF
END IF
CALL dcore_nucint_ana(sepi,sepj,rij,core=core,dcore=dcore,itype=itype,taper=taper,&
CALL dcore_nucint_ana(sepi,sepj,rij,core=core,dcore=dcore,itype=itype,se_taper=se_taper,&
se_int_control=se_int_control,lgrad=lgrad,error=error)
! Copy parameters over to arrays for do loop.
@ -324,7 +325,7 @@ CONTAINS
! Possibly debug the analytical values versus the numerical ones
IF (debug_this_module) THEN
CALL check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, taper, e1b, e2a, de1b, de2a, error)
CALL check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, e1b, e2a, de1b, de2a, error)
END IF
END IF
END SUBROUTINE rotnuc_ana
@ -354,7 +355,7 @@ CONTAINS
!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the core-core part
!> \author Teodoro Laino [tlaino] - Zurich University
! *****************************************************************************
RECURSIVE SUBROUTINE corecore_ana (sepi,sepj,rijv,itype,enuc,denuc,se_int_control,taper,error)
RECURSIVE SUBROUTINE corecore_ana (sepi,sepj,rijv,itype,enuc,denuc,se_int_control,se_taper,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
INTEGER, INTENT(IN) :: itype
@ -362,7 +363,7 @@ CONTAINS
REAL(dp), DIMENSION(3), INTENT(OUT), &
OPTIONAL :: denuc
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'corecore_ana', &
@ -389,12 +390,12 @@ CONTAINS
rij = SQRT(rij)
CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE.,&
integral_screening=se_int_control%integral_screening)
CALL dssss_nucint_ana(sepi,sepj,rij,ssss=ssss,dssss=dssss,itype=itype,taper=taper,&
CALL dssss_nucint_ana(sepi,sepj,rij,ssss=ssss,dssss=dssss,itype=itype,se_taper=se_taper,&
se_int_control=se_int_control_off,lgrad=l_denuc,error=error)
! In case let's compute the short-range part of the (ss|ss) integral
IF (se_int_control%shortrange) THEN
CALL dssss_nucint_ana(sepi,sepj,rij,ssss=ssss_sr,dssss=dssss_sr,itype=itype,&
taper=taper, se_int_control=se_int_control, lgrad=l_denuc, error=error)
se_taper=se_taper, se_int_control=se_int_control, lgrad=l_denuc, error=error)
ELSE
ssss_sr = ssss
dssss_sr = dssss
@ -639,7 +640,7 @@ CONTAINS
END IF
! Debug statement
IF (debug_this_module) THEN
CALL check_dcorecore_ana(sepi, sepj, rijv, itype, se_int_control, taper, enuc, denuc, error)
CALL check_dcorecore_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, enuc, denuc, error)
END IF
ENDIF
END SUBROUTINE corecore_ana
@ -826,13 +827,13 @@ CONTAINS
!> for computing integrals
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
SUBROUTINE dssss_nucint_ana(sepi, sepj, rij, ssss, dssss, itype, taper, se_int_control,&
SUBROUTINE dssss_nucint_ana(sepi, sepj, rij, ssss, dssss, itype, se_taper, se_int_control,&
lgrad, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), INTENT(IN) :: rij
REAL(dp), INTENT(OUT) :: ssss, dssss
INTEGER, INTENT(IN) :: itype
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
LOGICAL, INTENT(IN) :: lgrad
TYPE(cp_error_type), INTENT(inout) :: error
@ -850,12 +851,19 @@ CONTAINS
ft = 1.0_dp
dft= 0.0_dp
IF (itype/=do_method_pchg) THEN
ft = taper_eval(taper, rij, error)
dft= dtaper_eval(taper, rij, error)
ft = taper_eval(se_taper%taper, rij, error)
dft= dtaper_eval(se_taper%taper, rij, error)
END IF
se_int_screen%ft = ft
se_int_screen%dft = dft
! Evaluate additional taper function for dumped integrals
IF (se_int_control%integral_screening==do_se_IS_kdso_d) THEN
se_int_screen%ft = 1.0_dp
se_int_screen%dft = 0.0_dp
IF (itype/=do_method_pchg) THEN
se_int_screen%ft = taper_eval(se_taper%taper_add, rij, error)
se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij, error)
END IF
END IF
! Value of the integrals for sp shell
CALL nucint_sp_num(sepi, sepj, rij, ssss=ssss, itype=itype, se_int_control=se_int_control,&
se_int_screen=se_int_screen, error=error)
@ -874,7 +882,7 @@ CONTAINS
! Debug Procedure.. Check valifity of analytical gradients of nucint
IF (debug_this_module.AND.lgrad) THEN
CALL check_dssss_nucint_ana(sepi,sepj,rij,dssss,itype,se_int_control,taper=taper,error=error)
CALL check_dssss_nucint_ana(sepi,sepj,rij,dssss,itype,se_int_control,se_taper=se_taper,error=error)
END IF
END SUBROUTINE dssss_nucint_ana
@ -903,13 +911,13 @@ CONTAINS
!> for computing integrals
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
SUBROUTINE dcore_nucint_ana ( sepi, sepj, rij, core, dcore, itype, taper, &
SUBROUTINE dcore_nucint_ana ( sepi, sepj, rij, core, dcore, itype, se_taper, &
se_int_control, lgrad, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), INTENT(IN) :: rij
REAL(dp), DIMENSION(10, 2), INTENT(OUT) :: core, dcore
INTEGER, INTENT(IN) :: itype
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
LOGICAL, INTENT(IN) :: lgrad
TYPE(cp_error_type), INTENT(inout) :: error
@ -927,11 +935,18 @@ CONTAINS
ft = 1.0_dp
dft= 0.0_dp
IF (itype/=do_method_pchg) THEN
ft = taper_eval(taper, rij, error)
dft= dtaper_eval(taper, rij, error)
ft = taper_eval(se_taper%taper, rij, error)
dft= dtaper_eval(se_taper%taper, rij, error)
END IF
! Evaluate additional taper function for dumped integrals
IF (se_int_control%integral_screening==do_se_IS_kdso_d) THEN
se_int_screen%ft = 1.0_dp
se_int_screen%dft = 0.0_dp
IF (itype/=do_method_pchg) THEN
se_int_screen%ft = taper_eval(se_taper%taper_add, rij, error)
se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij, error)
END IF
END IF
se_int_screen%ft = ft
se_int_screen%dft = dft
! Value of the integrals for sp shell
CALL nucint_sp_num(sepi, sepj, rij, core=core, itype=itype,&
@ -975,7 +990,7 @@ CONTAINS
! Debug Procedure.. Check valifity of analytical gradients of nucint
IF (debug_this_module.AND.lgrad) THEN
CALL check_dcore_nucint_ana(sepi,sepj,rij,dcore,itype,se_int_control,taper=taper,error=error)
CALL check_dcore_nucint_ana(sepi,sepj,rij,dcore,itype,se_int_control,se_taper=se_taper,error=error)
END IF
END SUBROUTINE dcore_nucint_ana
@ -1177,7 +1192,7 @@ CONTAINS
!> for computing integrals
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
RECURSIVE SUBROUTINE rotint_ana (sepi,sepj,rijv,w,dw,se_int_control,taper,error)
RECURSIVE SUBROUTINE rotint_ana (sepi,sepj,rijv,w,dw,se_int_control,se_taper,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
REAL(dp), DIMENSION(2025), INTENT(OUT), &
@ -1185,7 +1200,7 @@ CONTAINS
REAL(dp), DIMENSION(3, 2025), &
INTENT(OUT), OPTIONAL :: dw
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'rotint_ana', &
@ -1224,7 +1239,7 @@ CONTAINS
CALL rotmat (sepi, sepj, rijv, rij, ij_matrix, do_derivatives=lgrad, do_invert=invert, error=error)
! Compute integrals in diatomic frame as well their derivatives (if requested)
CALL dterep_ana(sepi,sepj,rij,rep,rep_d,taper,se_int_control,lgrad=lgrad,error=error)
CALL dterep_ana(sepi,sepj,rij,rep,rep_d,se_taper,se_int_control,lgrad=lgrad,error=error)
IF (lgrad) THEN
drij(1) = rijv(1)/rij
@ -1241,7 +1256,7 @@ CONTAINS
ii = sepi%natorb
kk = sepj%natorb
! First step in rotation of integrals
CALL rot_2el_2c_first(sepi, sepj, rijv, se_int_control, taper, invert, ii, kk, rep, logv, ij_matrix,&
CALL rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, rep, logv, ij_matrix,&
v, lgrad, rep_d, v_d, logv_d, drij, error)
! Integrals if requested
@ -1339,7 +1354,7 @@ CONTAINS
IF (debug_this_module) THEN
! Check value of integrals
CALL check_rotint_ana(sepi,sepj,rijv,w,se_int_control=se_int_control,taper=taper,error=error)
CALL check_rotint_ana(sepi,sepj,rijv,w,se_int_control=se_int_control,se_taper=se_taper,error=error)
END IF
END IF
@ -1486,7 +1501,7 @@ CONTAINS
IF (debug_this_module) THEN
! Check derivatives
CALL check_rotint_ana(sepi,sepj,rijv,dw=dw,se_int_control=se_int_control,taper=taper,error=error)
CALL check_rotint_ana(sepi,sepj,rijv,dw=dw,se_int_control=se_int_control,se_taper=se_taper,error=error)
END IF
END IF
CALL rotmat_release(ij_matrix, error)
@ -1509,13 +1524,13 @@ CONTAINS
!> 03.2008 created [tlaino]
!> \author Teodoro Laino [tlaino] - Zurich University
! *****************************************************************************
RECURSIVE SUBROUTINE dterep_ana (sepi, sepj, rij, rep, rep_d, taper, se_int_control,&
lgrad, error)
RECURSIVE SUBROUTINE dterep_ana (sepi, sepj, rij, rep, rep_d, se_taper, &
se_int_control, lgrad, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(KIND=dp), INTENT(IN) :: rij
REAL(KIND=dp), DIMENSION(491), &
INTENT(OUT) :: rep, rep_d
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
LOGICAL, INTENT(IN) :: lgrad
TYPE(cp_error_type), INTENT(inout) :: error
@ -1531,16 +1546,18 @@ CONTAINS
failure = .FALSE.
! Compute the tapering function and its derivatives
ft = taper_eval(taper,rij,error)
ft = taper_eval(se_taper%taper,rij,error)
dft = 0.0_dp
ft1 = ft
IF (lgrad) THEN
ft1 = 1.0_dp
dft = dtaper_eval(taper,rij,error)
dft = dtaper_eval(se_taper%taper,rij,error)
END IF
! Evaluate additional taper function for dumped integrals
IF (se_int_control%integral_screening==do_se_IS_kdso_d) THEN
se_int_screen%ft = taper_eval(se_taper%taper_add, rij, error)
se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij, error)
END IF
se_int_screen%ft = ft
se_int_screen%dft = dft
! Integral Values for sp shells only
CALL terep_sp_num(sepi, sepj, rij, rep, se_int_control=se_int_control, &
@ -1586,7 +1603,7 @@ CONTAINS
! Possibly debug 2el 2cent integrals and derivatives
IF (debug_this_module) THEN
CALL check_dterep_ana(sepi, sepj, rij, rep, rep_d, se_int_control, taper=taper,&
CALL check_dterep_ana(sepi, sepj, rij, rep, rep_d, se_int_control, se_taper=se_taper,&
lgrad=lgrad, error=error)
END IF
END SUBROUTINE dterep_ana

View file

@ -135,7 +135,7 @@ END SUBROUTINE check_rotmat_der
!> 04.2008 created [tlaino]
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
SUBROUTINE rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, taper, invert, ii, kk, v_d, error)
SUBROUTINE rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, v_d, error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_int_utils, ONLY: rotmat
@ -145,8 +145,8 @@ SUBROUTINE rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, taper, inver
rotmat_type,&
rotmat_create,&
rotmat_release,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
USE semi_empirical_int_utils, ONLY: rot_2el_2c_first
#include "cp_common_uses.h"
IMPLICIT NONE
@ -154,7 +154,7 @@ SUBROUTINE rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, taper, inver
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rijv
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
LOGICAL, INTENT(IN) :: invert
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
INTEGER, INTENT(IN) :: ii, kk
REAL(KIND=dp), DIMENSION(3, 45, 45), &
INTENT(IN) :: v_d
@ -198,14 +198,14 @@ SUBROUTINE rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, taper, inver
error=error)
! Compute integrals in diatomic frame
CALL terep_num(sepi,sepj,r,rep,taper=taper,se_int_control=se_int_control,error=error)
CALL terep_num(sepi,sepj,r,rep,se_taper=se_taper,se_int_control=se_int_control,error=error)
IF (i==1) THEN
CALL rot_2el_2c_first(sepi, sepj, r0, se_int_control, taper, invert, ii, kk, rep, logv, matrix,&
CALL rot_2el_2c_first(sepi, sepj, r0, se_int_control, se_taper, invert, ii, kk, rep, logv, matrix,&
v_p, lgrad=.FALSE., error=error)
END IF
IF (i==2) THEN
CALL rot_2el_2c_first(sepi, sepj, r0, se_int_control, taper, invert, ii, kk, rep, logv, matrix,&
CALL rot_2el_2c_first(sepi, sepj, r0, se_int_control, se_taper, invert, ii, kk, rep, logv, matrix,&
v_m, lgrad=.FALSE., error=error)
END IF
CALL rotmat_release(matrix, error)
@ -232,13 +232,13 @@ END SUBROUTINE rot_2el_2c_first_debug
!> 04.2008 created [tlaino]
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
SUBROUTINE check_dssss_nucint_ana (sepi,sepj,r,dssss,itype,se_int_control,taper,error)
SUBROUTINE check_dssss_nucint_ana (sepi,sepj,r,dssss,itype,se_int_control,se_taper,error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_int_num, ONLY: ssss_nucint_num
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
@ -246,7 +246,7 @@ SUBROUTINE check_dssss_nucint_ana (sepi,sepj,r,dssss,itype,se_int_control,taper,
REAL(dp), INTENT(IN) :: dssss
INTEGER, INTENT(IN) :: itype
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: moduleN='semi_empirical_int_debug',&
@ -259,9 +259,9 @@ SUBROUTINE check_dssss_nucint_ana (sepi,sepj,r,dssss,itype,se_int_control,taper,
delta = 1.0E-8_dp
od = 0.5_dp/delta
rn = r + delta
CALL ssss_nucint_num(sepi,sepj,rn,ssssp,itype,taper,se_int_control,error=error)
CALL ssss_nucint_num(sepi,sepj,rn,ssssp,itype,se_taper,se_int_control,error=error)
rn = r - delta
CALL ssss_nucint_num(sepi,sepj,rn,ssssm,itype,taper,se_int_control,error=error)
CALL ssss_nucint_num(sepi,sepj,rn,ssssm,itype,se_taper,se_int_control,error=error)
nssss = od * (ssssp - ssssm)
! check
WRITE(*,*)"DEBUG::"//routineP
@ -279,13 +279,13 @@ END SUBROUTINE check_dssss_nucint_ana
!> 04.2008 created [tlaino]
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
SUBROUTINE check_dcore_nucint_ana (sepi,sepj,r,dcore,itype,se_int_control,taper,error)
SUBROUTINE check_dcore_nucint_ana (sepi,sepj,r,dcore,itype,se_int_control,se_taper,error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_int_num, ONLY: core_nucint_num
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
@ -293,7 +293,7 @@ SUBROUTINE check_dcore_nucint_ana (sepi,sepj,r,dcore,itype,se_int_control,taper,
REAL(dp), DIMENSION(10, 2), INTENT(IN) :: dcore
INTEGER, INTENT(IN) :: itype
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: moduleN='semi_empirical_int_debug',&
@ -307,9 +307,9 @@ SUBROUTINE check_dcore_nucint_ana (sepi,sepj,r,dcore,itype,se_int_control,taper,
delta = 1.0E-8_dp
od = 0.5_dp/delta
rn = r + delta
CALL core_nucint_num(sepi,sepj,rn,corep,itype,taper,se_int_control,error=error)
CALL core_nucint_num(sepi,sepj,rn,corep,itype,se_taper,se_int_control,error=error)
rn = r - delta
CALL core_nucint_num(sepi,sepj,rn,corem,itype,taper,se_int_control,error=error)
CALL core_nucint_num(sepi,sepj,rn,corem,itype,se_taper,se_int_control,error=error)
ncore = od * (corep - corem)
! check
WRITE(*,*)"DEBUG::"//routineP
@ -365,21 +365,21 @@ END FUNCTION check_value
!> 04.2008 created [tlaino]
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
SUBROUTINE check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, taper, e1b, e2a, de1b, de2a, error)
SUBROUTINE check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, e1b, e2a, de1b, de2a, error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_int_num, ONLY: rotnuc_num,&
drotnuc_num
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
INTEGER, INTENT(IN) :: itype
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
REAL(dp), DIMENSION(45), INTENT(IN), &
OPTIONAL :: e1b, e2a
REAL(dp), DIMENSION(3, 45), &
@ -404,7 +404,7 @@ SUBROUTINE check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, taper, e1b
delta = 1.0E-5_dp
! Check value of integrals
WRITE(*,*)"DEBUG::"//routineP
CALL rotnuc_num(sepi,sepj,rijv,e1b2,e2a2,itype,se_int_control,taper=taper,error=error)
CALL rotnuc_num(sepi,sepj,rijv,e1b2,e2a2,itype,se_int_control,se_taper=se_taper,error=error)
IF (l_e1b) THEN
DO j = 1, 45
IF (.NOT.check_value(e1b2(j), e1b(j), delta, 0.1_dp)) THEN
@ -425,7 +425,7 @@ SUBROUTINE check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, taper, e1b
! Check derivatives
IF (lgrad) THEN
CALL drotnuc_num(sepi,sepj,rijv,de1b2,de2a2,itype,delta=delta,&
se_int_control=se_int_control,taper=taper,error=error)
se_int_control=se_int_control,se_taper=se_taper,error=error)
IF (l_de1b) THEN
DO i = 1, 3
DO j = 1, 45
@ -463,14 +463,14 @@ END SUBROUTINE check_drotnuc_ana
!> 04.2007 created [tlaino]
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
SUBROUTINE check_dcorecore_ana(sepi, sepj, rijv, itype,se_int_control, taper, enuc, denuc, error)
SUBROUTINE check_dcorecore_ana(sepi, sepj, rijv, itype,se_int_control, se_taper, enuc, denuc, error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_int_num, ONLY: corecore_num,&
dcorecore_num
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
@ -480,7 +480,7 @@ SUBROUTINE check_dcorecore_ana(sepi, sepj, rijv, itype,se_int_control, taper, en
REAL(dp), DIMENSION(3), INTENT(IN), &
OPTIONAL :: denuc
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', &
@ -495,14 +495,14 @@ SUBROUTINE check_dcorecore_ana(sepi, sepj, rijv, itype,se_int_control, taper, en
delta = 1.0E-7_dp
! check
IF (PRESENT(enuc)) THEN
CALL corecore_num(sepi,sepj,rijv, enuc_num,itype,se_int_control,taper,error)
CALL corecore_num(sepi,sepj,rijv, enuc_num,itype,se_int_control,se_taper,error)
IF (.NOT.check_value(enuc, enuc_num, delta, 0.001_dp)) THEN
WRITE(*,*)"ERROR for CORE-CORE energy value (numerical different from analytical)!!"
STOP
END IF
END IF
IF (PRESENT(denuc)) THEN
CALL dcorecore_num(sepi,sepj,rijv,denuc_num,itype,delta,se_int_control,taper,error)
CALL dcorecore_num(sepi,sepj,rijv,denuc_num,itype,delta,se_int_control,se_taper,error)
DO j = 1, 3
IF (.NOT.check_value(denuc(j), denuc_num(j), delta, 0.001_dp)) THEN
WRITE(*,*)"ERROR for CORE-CORE energy derivative value (numerical different from analytical). DENUC(j), j::", j
@ -520,13 +520,13 @@ END SUBROUTINE check_dcorecore_ana
!> 04.2007 created [tlaino]
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
SUBROUTINE check_dterep_ana (sepi,sepj,r,ri,dri,se_int_control,taper,lgrad,error)
SUBROUTINE check_dterep_ana (sepi,sepj,r,ri,dri,se_int_control,se_taper,lgrad,error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_int_num, ONLY: terep_num
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
@ -534,7 +534,7 @@ SUBROUTINE check_dterep_ana (sepi,sepj,r,ri,dri,se_int_control,taper,lgrad,error
REAL(dp), DIMENSION(491), INTENT(IN) :: ri, dri
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
LOGICAL, INTENT(IN) :: lgrad
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', &
@ -548,12 +548,12 @@ SUBROUTINE check_dterep_ana (sepi,sepj,r,ri,dri,se_int_control,taper,lgrad,error
delta = 1.0E-8_dp
od = 0.5_dp/delta
rn = r
CALL terep_num(sepi,sepj,rn,ri0,taper,se_int_control,error=error)
CALL terep_num(sepi,sepj,rn,ri0,se_taper,se_int_control,error=error)
IF (lgrad) THEN
rn = r + delta
CALL terep_num(sepi,sepj,rn,rip,taper,se_int_control,error=error)
CALL terep_num(sepi,sepj,rn,rip,se_taper,se_int_control,error=error)
rn = r - delta
CALL terep_num(sepi,sepj,rn,rim,taper,se_int_control,error=error)
CALL terep_num(sepi,sepj,rn,rim,se_taper,se_int_control,error=error)
nri = od * (rip - rim)
END IF
! check
@ -584,14 +584,14 @@ END SUBROUTINE check_dterep_ana
!> 04.2008 created [tlaino]
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
SUBROUTINE check_rotint_ana(sepi,sepj,rijv,w,dw,se_int_control,taper,error)
SUBROUTINE check_rotint_ana(sepi,sepj,rijv,w,dw,se_int_control,se_taper,error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_int_num, ONLY: rotint_num,&
drotint_num
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
@ -601,7 +601,7 @@ SUBROUTINE check_rotint_ana(sepi,sepj,rijv,w,dw,se_int_control,taper,error)
REAL(dp), DIMENSION(3, 2025), &
INTENT(IN), OPTIONAL :: dw
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: moduleN = 'semi_empirical_int_debug', &
@ -617,7 +617,7 @@ SUBROUTINE check_rotint_ana(sepi,sepj,rijv,w,dw,se_int_control,taper,error)
WRITE(*,*)"DEBUG::"//routineP
IF (PRESENT(w)) THEN
w2=0.0_dp
CALL rotint_num(sepi,sepj,rijv,w2,se_int_control,taper=taper,error=error)
CALL rotint_num(sepi,sepj,rijv,w2,se_int_control,se_taper=se_taper,error=error)
DO j = 1, 2025
IF (.NOT.check_value(w(j), w2(j), delta, 0.1_dp)) THEN
WRITE(*,*)"ERROR for integral value W(j), j::",j
@ -630,8 +630,8 @@ SUBROUTINE check_rotint_ana(sepi,sepj,rijv,w,dw,se_int_control,taper,error)
! First of all let's decide if the value we get for delta is compatible
! with a reasonable value of the integral.. (compatible if the value of the
! integral is greater than 1.0E-6)
CALL drotint_num(sepi,sepj,rijv,dw2,delta=delta,se_int_control=se_int_control,taper=taper,error=error)
CALL rotint_num(sepi,sepj,rijv,w2,se_int_control=se_int_control,taper=taper,error=error)
CALL drotint_num(sepi,sepj,rijv,dw2,delta=delta,se_int_control=se_int_control,se_taper=se_taper,error=error)
CALL rotint_num(sepi,sepj,rijv,w2,se_int_control=se_int_control,se_taper=se_taper,error=error)
DO i = 1, 3
DO j = 1, 2025
IF ((ABS(w2(j))>delta).AND.(ABS(dw2(i,j))>delta*10)) THEN

View file

@ -34,12 +34,13 @@ END INTERFACE check_rotmat_der
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
INTERFACE check_dssss_nucint_ana
SUBROUTINE check_dssss_nucint_ana (sepi,sepj,r,dssss,itype,se_int_control,taper,error)
SUBROUTINE check_dssss_nucint_ana (sepi,sepj,r,dssss,itype,se_int_control,&
se_taper,error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
@ -47,7 +48,7 @@ INTERFACE check_dssss_nucint_ana
REAL(dp), INTENT(IN) :: dssss
INTEGER, INTENT(IN) :: itype
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
END SUBROUTINE check_dssss_nucint_ana
END INTERFACE check_dssss_nucint_ana
@ -61,12 +62,13 @@ END INTERFACE check_dssss_nucint_ana
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
INTERFACE check_dcore_nucint_ana
SUBROUTINE check_dcore_nucint_ana (sepi,sepj,r,dcore,itype,se_int_control,taper,error)
SUBROUTINE check_dcore_nucint_ana (sepi,sepj,r,dcore,itype,se_int_control,&
se_taper,error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
@ -74,7 +76,7 @@ INTERFACE check_dcore_nucint_ana
REAL(dp), DIMENSION(10, 2), INTENT(IN) :: dcore
INTEGER, INTENT(IN) :: itype
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
END SUBROUTINE check_dcore_nucint_ana
END INTERFACE check_dcore_nucint_ana
@ -88,19 +90,20 @@ END INTERFACE check_dcore_nucint_ana
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
INTERFACE check_drotnuc_ana
SUBROUTINE check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, taper, e1b, e2a, de1b, de2a, error)
SUBROUTINE check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, se_taper,&
e1b, e2a, de1b, de2a, error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
INTEGER, INTENT(IN) :: itype
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
REAL(dp), DIMENSION(45), INTENT(IN), &
OPTIONAL :: e1b, e2a
REAL(dp), DIMENSION(45, 3), &
@ -118,12 +121,13 @@ END INTERFACE check_drotnuc_ana
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
INTERFACE check_dcorecore_ana
SUBROUTINE check_dcorecore_ana(sepi, sepj, rijv, itype,se_int_control, taper, enuc, denuc, error)
SUBROUTINE check_dcorecore_ana(sepi, sepj, rijv, itype,se_int_control,&
se_taper, enuc, denuc, error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
@ -133,7 +137,7 @@ INTERFACE check_dcorecore_ana
REAL(dp), DIMENSION(3), INTENT(IN), &
OPTIONAL :: denuc
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
END SUBROUTINE check_dcorecore_ana
@ -148,19 +152,20 @@ END INTERFACE check_dcorecore_ana
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
INTERFACE rot_2el_2c_first_debug
SUBROUTINE rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, taper, invert, ii, kk, v_d, error)
SUBROUTINE rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, se_taper,&
invert, ii, kk, v_d, error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rijv
LOGICAL, INTENT(IN) :: invert
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
INTEGER, INTENT(IN) :: ii, kk
REAL(KIND=dp), DIMENSION(45, 45, 3), &
INTENT(IN) :: v_d
@ -177,12 +182,12 @@ END INTERFACE rot_2el_2c_first_debug
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
INTERFACE check_dterep_ana
SUBROUTINE check_dterep_ana (sepi,sepj,r,ri,dri,se_int_control,taper,lgrad,error)
SUBROUTINE check_dterep_ana (sepi,sepj,r,ri,dri,se_int_control,se_taper,lgrad,error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
@ -190,7 +195,7 @@ INTERFACE check_dterep_ana
REAL(dp), DIMENSION(491), INTENT(IN) :: ri, dri
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
LOGICAL, INTENT(IN) :: lgrad
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
END SUBROUTINE check_dterep_ana
END INTERFACE check_dterep_ana
@ -204,12 +209,12 @@ END INTERFACE check_dterep_ana
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
INTERFACE check_rotint_ana
SUBROUTINE check_rotint_ana(sepi,sepj,rijv,w,dw,se_int_control,taper,error)
SUBROUTINE check_rotint_ana(sepi,sepj,rijv,w,dw,se_int_control,se_taper,error)
USE f77_blas
USE kinds, ONLY: dp
USE semi_empirical_types, ONLY: semi_empirical_type,&
se_int_control_type
USE taper_types, ONLY: taper_type
se_int_control_type,&
se_taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
TYPE(semi_empirical_type), POINTER :: sepi, sepj
@ -219,7 +224,7 @@ INTERFACE check_rotint_ana
REAL(dp), DIMENSION(2025, 3), &
INTENT(IN), OPTIONAL :: dw
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
END SUBROUTINE check_rotint_ana
END INTERFACE check_rotint_ana

View file

@ -20,7 +20,7 @@ MODULE semi_empirical_int_num
do_method_pdg,&
do_method_pm3,&
do_method_pm6,&
do_method_undef
do_method_undef,do_se_IS_kdso_d
USE kinds, ONLY: dp
USE physcon, ONLY: angstrom,&
evolt
@ -38,9 +38,9 @@ MODULE semi_empirical_int_num
se_int_control_type,&
se_int_screen_type,&
semi_empirical_type,&
setup_se_int_control_type
USE taper_types, ONLY: taper_eval,&
taper_type
setup_se_int_control_type,&
se_taper_type
USE taper_types, ONLY: taper_eval
#include "cp_common_uses.h"
#include "semi_empirical_int_args.h"
IMPLICIT NONE
@ -70,12 +70,12 @@ CONTAINS
!> written by Ernest R. Davidson, Indiana University.
!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : major rewriting
! *****************************************************************************
SUBROUTINE rotint_num (sepi,sepj,rijv,w,se_int_control,taper,error)
SUBROUTINE rotint_num (sepi,sepj,rijv,w,se_int_control,se_taper,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
REAL(dp), DIMENSION(2025), INTENT(OUT) :: w
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'rotint_num', &
@ -106,7 +106,7 @@ CONTAINS
CALL rotmat (sepi, sepj, rijv, rij, ij_matrix, do_derivatives=.FALSE., error=error)
! Compute Integrals in Diatomic Frame
CALL terep_num(sepi,sepj,rij,rep,taper=taper,se_int_control=se_int_control,error=error)
CALL terep_num(sepi,sepj,rij,rep,se_taper=se_taper,se_int_control=se_int_control,error=error)
! Rotate Integrals
ii = sepi%natorb
@ -120,7 +120,7 @@ CONTAINS
END DO
! First step in rotation of integrals
CALL rot_2el_2c_first(sepi, sepj, rijv, se_int_control, taper, .FALSE., ii, kk, rep, &
CALL rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, .FALSE., ii, kk, rep, &
logv, ij_matrix, v, lgrad=.FALSE.,error=error)
! Second step in rotation of integrals
@ -224,11 +224,11 @@ CONTAINS
!> 03.2008 created [tlaino]
!> \author Teodoro Laino [tlaino] - Zurich University
! *****************************************************************************
SUBROUTINE terep_num(sepi, sepj, rij, rep, taper, se_int_control, error)
SUBROUTINE terep_num(sepi, sepj, rij, rep, se_taper, se_int_control, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), INTENT(IN) :: rij
REAL(dp), DIMENSION(491), INTENT(OUT) :: rep
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(cp_error_type), INTENT(inout) :: error
@ -240,8 +240,11 @@ CONTAINS
TYPE(se_int_screen_type) :: se_int_screen
failure = .FALSE.
ft = taper_eval(taper, rij, error)
se_int_screen%ft = ft
ft = taper_eval(se_taper%taper, rij, error)
! In case of dumped integrals compute an additional taper term
IF (se_int_control%integral_screening==do_se_IS_kdso_d) THEN
se_int_screen%ft = taper_eval(se_taper%taper_add, rij, error)
END IF
! Contribution from sp shells
CALL terep_sp_num(sepi, sepj, rij, rep, se_int_control, se_int_screen, ft,&
@ -412,14 +415,14 @@ CONTAINS
!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : major rewriting
!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the core-core part
! *****************************************************************************
SUBROUTINE rotnuc_num(sepi,sepj,rijv,e1b,e2a,itype,se_int_control,taper,error)
SUBROUTINE rotnuc_num(sepi,sepj,rijv,e1b,e2a,itype,se_int_control,se_taper,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
REAL(dp), DIMENSION(45), INTENT(OUT), &
OPTIONAL :: e1b, e2a
INTEGER, INTENT(IN) :: itype
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'rotnuc_num', &
@ -445,7 +448,7 @@ CONTAINS
CALL rotmat (sepi, sepj, rijv, rij, ij_matrix, do_derivatives=.FALSE., error=error)
! Compute Integrals in Diatomic Frame
CALL core_nucint_num(sepi,sepj,rij,core=core,itype=itype,taper=taper,&
CALL core_nucint_num(sepi,sepj,rij,core=core,itype=itype,se_taper=se_taper,&
se_int_control=se_int_control,error=error)
! Copy parameters over to arrays for do loop.
@ -523,13 +526,13 @@ CONTAINS
!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : major rewriting
!> Teodoro Laino [tlaino] - University of Zurich 04.2008 : splitted from rotnuc
! *****************************************************************************
SUBROUTINE corecore_num(sepi,sepj,rijv,enuc,itype,se_int_control,taper,error)
SUBROUTINE corecore_num(sepi,sepj,rijv,enuc,itype,se_int_control,se_taper,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
REAL(dp), INTENT(OUT) :: enuc
INTEGER, INTENT(IN) :: itype
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'corecore_num', &
@ -549,11 +552,11 @@ CONTAINS
rij = SQRT(rij)
CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE.,&
integral_screening=se_int_control%integral_screening)
CALL ssss_nucint_num(sepi,sepj,rij,ssss=ssss,itype=itype,taper=taper,&
CALL ssss_nucint_num(sepi,sepj,rij,ssss=ssss,itype=itype,se_taper=se_taper,&
se_int_control=se_int_control_off,error=error)
! In case let's compute the short-range part of the (ss|ss) integral
IF (se_int_control%shortrange) THEN
CALL ssss_nucint_num(sepi,sepj,rij,ssss=ssss_sr,itype=itype,taper=taper,&
CALL ssss_nucint_num(sepi,sepj,rij,ssss=ssss_sr,itype=itype,se_taper=se_taper,&
se_int_control=se_int_control,error=error)
ELSE
ssss_sr = ssss
@ -704,12 +707,12 @@ CONTAINS
!> 03.2008 created [tlaino]
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
SUBROUTINE ssss_nucint_num(sepi, sepj, rij, ssss, itype, taper, se_int_control, error)
SUBROUTINE ssss_nucint_num(sepi, sepj, rij, ssss, itype, se_taper, se_int_control, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), INTENT(IN) :: rij
REAL(dp), INTENT(OUT) :: ssss
INTEGER, INTENT(IN) :: itype
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(cp_error_type), INTENT(inout) :: error
@ -723,9 +726,18 @@ CONTAINS
failure = .FALSE.
! Computing Tapering function
ft = 1.0_dp
IF (itype /=do_method_pchg) ft = taper_eval (taper, rij, error)
se_int_screen%ft = ft
ft = 1.0_dp
IF (itype /=do_method_pchg) THEN
ft = taper_eval(se_taper%taper, rij, error)
END IF
! In case of dumped integrals compute an additional taper term
IF (se_int_control%integral_screening==do_se_IS_kdso_d) THEN
se_int_screen%ft = 1.0_dp
IF (itype /=do_method_pchg) THEN
se_int_screen%ft = taper_eval(se_taper%taper_add, rij, error)
END IF
END IF
! Contribution from the sp shells
CALL nucint_sp_num(sepi, sepj, rij, ssss=ssss, itype=itype,&
@ -753,12 +765,12 @@ CONTAINS
!> 03.2008 created [tlaino]
!> \author Teodoro Laino - Zurich University
! *****************************************************************************
SUBROUTINE core_nucint_num(sepi, sepj, rij, core, itype, taper, se_int_control, error)
SUBROUTINE core_nucint_num(sepi, sepj, rij, core, itype, se_taper, se_int_control, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), INTENT(IN) :: rij
REAL(dp), DIMENSION(10, 2), INTENT(OUT) :: core
INTEGER, INTENT(IN) :: itype
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(cp_error_type), INTENT(inout) :: error
@ -772,9 +784,18 @@ CONTAINS
failure = .FALSE.
! Computing the Tapering function
ft = 1.0_dp
IF (itype /=do_method_pchg) ft = taper_eval (taper, rij, error)
se_int_screen%ft = ft
ft = 1.0_dp
IF (itype /=do_method_pchg) THEN
ft = taper_eval(se_taper%taper, rij, error)
END IF
! In case of dumped integrals compute an additional taper term
IF (se_int_control%integral_screening==do_se_IS_kdso_d) THEN
se_int_screen%ft = 1.0_dp
IF (itype /=do_method_pchg) THEN
se_int_screen%ft = taper_eval(se_taper%taper_add, rij, error)
END IF
END IF
! Contribution from the sp shells
CALL nucint_sp_num(sepi, sepj, rij, core=core, itype=itype,&
@ -978,14 +999,14 @@ CONTAINS
! *****************************************************************************
!> \brief Numerical Derivatives for rotint
! *****************************************************************************
SUBROUTINE drotint_num(sepi,sepj,r,dw,delta,se_int_control,taper, error)
SUBROUTINE drotint_num(sepi,sepj,r,dw,delta,se_int_control,se_taper, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: r
REAL(dp), DIMENSION(3, 2025), &
INTENT(OUT) :: dw
REAL(dp), INTENT(IN) :: delta
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'drotint_num', &
@ -1001,9 +1022,9 @@ CONTAINS
DO i=1,3
rr=r
rr(i)=rr(i)+delta
CALL rotint_num(sepi,sepj,rr,wp,se_int_control,taper=taper,error=error)
CALL rotint_num(sepi,sepj,rr,wp,se_int_control,se_taper=se_taper,error=error)
rr(i)=rr(i)-2._dp*delta
CALL rotint_num(sepi,sepj,rr,wm,se_int_control,taper=taper,error=error)
CALL rotint_num(sepi,sepj,rr,wm,se_int_control,se_taper=se_taper,error=error)
DO j = 1, nsize
dw(i,j) = od * (wp(j) - wm(j))
END DO
@ -1014,7 +1035,7 @@ CONTAINS
! *****************************************************************************
!> \brief Numerical Derivatives for rotnuc
! *****************************************************************************
SUBROUTINE drotnuc_num(sepi,sepj,r,de1b,de2a,itype,delta,se_int_control,taper,error)
SUBROUTINE drotnuc_num(sepi,sepj,r,de1b,de2a,itype,delta,se_int_control,se_taper,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: r
REAL(dp), DIMENSION(3, 45), &
@ -1022,7 +1043,7 @@ CONTAINS
INTEGER, INTENT(IN) :: itype
REAL(dp), INTENT(IN) :: delta
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'drotnuc_num', &
@ -1041,9 +1062,9 @@ CONTAINS
DO i=1,3
rr=r
rr(i)=rr(i)+delta
CALL rotnuc_num(sepi,sepj,rr,e1p,e2p,itype,se_int_control,taper=taper,error=error)
CALL rotnuc_num(sepi,sepj,rr,e1p,e2p,itype,se_int_control,se_taper=se_taper,error=error)
rr(i)=rr(i)-2._dp*delta
CALL rotnuc_num(sepi,sepj,rr,e1m,e2m,itype,se_int_control,taper=taper,error=error)
CALL rotnuc_num(sepi,sepj,rr,e1m,e2m,itype,se_int_control,se_taper=se_taper,error=error)
IF (l_de1b) THEN
DO j = 1, sepi%atm_int_size
de1b(i,j) = od * (e1p(j) - e1m(j))
@ -1060,14 +1081,14 @@ CONTAINS
! *****************************************************************************
!> \brief Numerical Derivatives for corecore
! *****************************************************************************
SUBROUTINE dcorecore_num(sepi,sepj,r,denuc,itype,delta,se_int_control,taper,error)
SUBROUTINE dcorecore_num(sepi,sepj,r,denuc,itype,delta,se_int_control,se_taper,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: r
REAL(dp), DIMENSION(3), INTENT(OUT) :: denuc
INTEGER, INTENT(IN) :: itype
REAL(dp), INTENT(IN) :: delta
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'dcorecore_num', &
@ -1081,9 +1102,9 @@ CONTAINS
DO i=1,3
rr=r
rr(i)=rr(i)+delta
CALL corecore_num(sepi,sepj,rr,enucp,itype,se_int_control,taper=taper,error=error)
CALL corecore_num(sepi,sepj,rr,enucp,itype,se_int_control,se_taper=se_taper,error=error)
rr(i)=rr(i)-2._dp*delta
CALL corecore_num(sepi,sepj,rr,enucm,itype,se_int_control,taper=taper,error=error)
CALL corecore_num(sepi,sepj,rr,enucm,itype,se_int_control,se_taper=se_taper,error=error)
denuc(i) = od * ( enucp - enucm )
END DO
END SUBROUTINE dcorecore_num

View file

@ -22,8 +22,8 @@ MODULE semi_empirical_int_utils
USE semi_empirical_types, ONLY: rotmat_type,&
se_int_control_type,&
se_int_screen_type,&
se_taper_type,&
semi_empirical_type
USE taper_types, ONLY: taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
@ -105,7 +105,7 @@ CONTAINS
itype, dcharg_int_nri, error)
! Handle the derivative of the Tapering term
IF (se_int_control%integral_screening==do_se_IS_kdso_d) THEN
IF ((se_int_control%integral_screening==do_se_IS_kdso_d).AND.(itype==do_method_pchg)) THEN
dfs = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,&
se_int_control%integral_screening, .FALSE., itype, dcharg_int_nri_fs,&
error)
@ -228,6 +228,7 @@ CONTAINS
dij = 0.0_dp
dkl = 0.0_dp
add = 0.0_dp
fact_screen = 0.0_dp
DO l1 = l1min, l1max
IF (l1 /= 0) THEN
dij = sepi%cs(lij)*fact_ij
@ -243,7 +244,7 @@ CONTAINS
ccc = clm_sp(ij, l1, m) * clm_sp(kl, l2, m)
IF (ABS(ccc) > EPSILON(0.0_dp)) THEN
mm = ABS (m)
chrg = eval(r, l1, l2, clm_sp(ij, l1, m), clm_sp(kl, l2, m), dij, dkl, add, error)
chrg = eval(r, l1, l2, clm_sp(ij, l1, m), clm_sp(kl, l2, m), dij, dkl, add, fact_screen, error)
s1 = s1 + chrg
END IF
END DO
@ -927,7 +928,7 @@ CONTAINS
itype, dcharg_int_ri, error)
! Handle the derivative of the Tapering term
IF (se_int_control%integral_screening==do_se_IS_kdso_d) THEN
IF ((se_int_control%integral_screening==do_se_IS_kdso_d).AND.(itype==do_method_pchg)) THEN
dfs = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,&
se_int_control%integral_screening, .FALSE., itype, dcharg_int_ri_fs,&
error)
@ -1053,6 +1054,7 @@ CONTAINS
dij = 0.0_dp
dkl = 0.0_dp
add = 0.0_dp
fact_screen = 0.0_dp
DO l1 = l1min, l1max
IF (l1 /= 0) THEN
dij = sepi%cs(lij)
@ -1068,7 +1070,7 @@ CONTAINS
ccc = clm_d(ij, l1, m) * clm_d(kl, l2, m)
IF (ABS(ccc) > EPSILON(0.0_dp)) THEN
mm = ABS (m)
chrg = eval(r, l1, l2, mm, dij, dkl, add, error)
chrg = eval(r, l1, l2, mm, dij, dkl, add, fact_screen, error)
s1 = s1 + chrg * ccc
END IF
END DO
@ -1884,13 +1886,13 @@ CONTAINS
!> \author Teodoro Laino [tlaino] - University of Zurich
!> \date 04.2008 [tlaino]
! *****************************************************************************
RECURSIVE SUBROUTINE rot_2el_2c_first (sepi, sepj, rijv, se_int_control, taper,&
RECURSIVE SUBROUTINE rot_2el_2c_first (sepi, sepj, rijv, se_int_control, se_taper,&
invert, ii, kk, rep, logv, ij_matrix, v, lgrad, rep_d, v_d, logv_d, drij,&
error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rijv
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
LOGICAL, INTENT(IN) :: invert
INTEGER, INTENT(IN) :: ii, kk
REAL(KIND=dp), DIMENSION(491), &
@ -2211,7 +2213,7 @@ CONTAINS
END DO
END DO
IF (debug_this_module) THEN
CALL rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, taper, invert, ii, kk, v_d, error)
CALL rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, v_d, error)
END IF
END IF
END SUBROUTINE rot_2el_2c_first

View file

@ -30,8 +30,8 @@ MODULE semi_empirical_integrals
rotnuc_num
USE semi_empirical_store_int_types, ONLY: semi_empirical_si_type
USE semi_empirical_types, ONLY: se_int_control_type,&
se_taper_type,&
semi_empirical_type
USE taper_types, ONLY: taper_type
#include "cp_common_uses.h"
IMPLICIT NONE
@ -49,13 +49,13 @@ CONTAINS
!> \author Teodoro Laino [tlaino] - University of Zurich
!> \date 05.2008
! *****************************************************************************
SUBROUTINE rotint(sepi,sepj,rij,w,anag,se_int_control,taper,store_int_env,error)
SUBROUTINE rotint(sepi,sepj,rij,w,anag,se_int_control,se_taper,store_int_env,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rij
REAL(dp), DIMENSION(2025), INTENT(OUT) :: w
LOGICAL :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(semi_empirical_si_type), POINTER :: store_int_env
TYPE(cp_error_type), INTENT(inout) :: error
@ -86,9 +86,9 @@ CONTAINS
END IF
! Compute Integrals
IF (anag) THEN
CALL rotint_ana (sepi,sepj,rij,w,se_int_control=se_int_control,taper=taper,error=error)
CALL rotint_ana (sepi,sepj,rij,w,se_int_control=se_int_control,se_taper=se_taper,error=error)
ELSE
CALL rotint_num (sepi,sepj,rij,w,se_int_control=se_int_control,taper=taper,error=error)
CALL rotint_num (sepi,sepj,rij,w,se_int_control=se_int_control,se_taper=se_taper,error=error)
END IF
! Store integrals if we did not go overflow
IF (.NOT.buffer_overflow) THEN
@ -162,17 +162,17 @@ CONTAINS
END IF
ELSE
IF (anag) THEN
CALL rotint_ana (sepi,sepj,rij,w,se_int_control=se_int_control,taper=taper,error=error)
CALL rotint_ana (sepi,sepj,rij,w,se_int_control=se_int_control,se_taper=se_taper,error=error)
ELSE
CALL rotint_num (sepi,sepj,rij,w,se_int_control=se_int_control,taper=taper,error=error)
CALL rotint_num (sepi,sepj,rij,w,se_int_control=se_int_control,se_taper=se_taper,error=error)
END IF
END IF
END IF
ELSE
IF (anag) THEN
CALL rotint_ana (sepi,sepj,rij,w,se_int_control=se_int_control,taper=taper,error=error)
CALL rotint_ana (sepi,sepj,rij,w,se_int_control=se_int_control,se_taper=se_taper,error=error)
ELSE
CALL rotint_num (sepi,sepj,rij,w,se_int_control=se_int_control,taper=taper,error=error)
CALL rotint_num (sepi,sepj,rij,w,se_int_control=se_int_control,se_taper=se_taper,error=error)
END IF
END IF
END SUBROUTINE rotint
@ -182,7 +182,7 @@ CONTAINS
!> \author Teodoro Laino [tlaino] - University of Zurich
!> \date 05.2008
! *****************************************************************************
SUBROUTINE rotnuc (sepi,sepj,rij,e1b,e2a,itype,anag,se_int_control,taper,store_int_env,error)
SUBROUTINE rotnuc (sepi,sepj,rij,e1b,e2a,itype,anag,se_int_control,se_taper,store_int_env,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rij
REAL(dp), DIMENSION(45), INTENT(OUT), &
@ -190,7 +190,7 @@ CONTAINS
INTEGER, INTENT(IN) :: itype
LOGICAL, INTENT(IN) :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(semi_empirical_si_type), OPTIONAL, &
POINTER :: store_int_env
TYPE(cp_error_type), INTENT(inout) :: error
@ -230,10 +230,10 @@ CONTAINS
! Compute Integrals
IF (anag) THEN
CALL rotnuc_ana (sepi, sepj, rij, e1b=e1b, e2a=e2a, itype=itype,&
se_int_control=se_int_control, taper=taper, error=error)
se_int_control=se_int_control, se_taper=se_taper, error=error)
ELSE
CALL rotnuc_num (sepi, sepj, rij, e1b=e1b, e2a=e2a, itype=itype,&
se_int_control=se_int_control, taper=taper, error=error)
se_int_control=se_int_control, se_taper=se_taper, error=error)
END IF
! Store integrals if we did not go overflow
IF (.NOT.buffer_overflow) THEN
@ -313,20 +313,20 @@ CONTAINS
ELSE
IF (anag) THEN
CALL rotnuc_ana (sepi, sepj, rij, e1b=e1b, e2a=e2a, itype=itype,&
se_int_control=se_int_control, taper=taper, error=error)
se_int_control=se_int_control, se_taper=se_taper, error=error)
ELSE
CALL rotnuc_num (sepi, sepj, rij, e1b=e1b, e2a=e2a, itype=itype,&
se_int_control=se_int_control, taper=taper, error=error)
se_int_control=se_int_control, se_taper=se_taper, error=error)
END IF
END IF
END IF
ELSE
IF (anag) THEN
CALL rotnuc_ana (sepi, sepj, rij, e1b=e1b, e2a=e2a, itype=itype,&
se_int_control=se_int_control, taper=taper, error=error)
se_int_control=se_int_control, se_taper=se_taper, error=error)
ELSE
CALL rotnuc_num (sepi, sepj, rij, e1b=e1b, e2a=e2a, itype=itype,&
se_int_control=se_int_control, taper=taper, error=error)
se_int_control=se_int_control, se_taper=se_taper, error=error)
END IF
END IF
@ -340,14 +340,14 @@ CONTAINS
!> \author Teodoro Laino [tlaino] - University of Zurich
!> \date 04.2008
! *****************************************************************************
SUBROUTINE corecore (sepi,sepj,rij,enuc,itype,anag,se_int_control,taper,error)
SUBROUTINE corecore (sepi,sepj,rij,enuc,itype,anag,se_int_control,se_taper,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rij
REAL(dp), INTENT(OUT) :: enuc
INTEGER, INTENT(IN) :: itype
LOGICAL, INTENT(IN) :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'corecore', &
@ -355,10 +355,10 @@ CONTAINS
IF (anag) THEN
CALL corecore_ana (sepi,sepj,rij,enuc=enuc,itype=itype,se_int_control=se_int_control,&
taper=taper, error=error)
se_taper=se_taper, error=error)
ELSE
CALL corecore_num (sepi,sepj,rij,enuc=enuc,itype=itype,se_int_control=se_int_control,&
taper=taper, error=error)
se_taper=se_taper, error=error)
END IF
END SUBROUTINE corecore
@ -368,7 +368,7 @@ CONTAINS
!> \author Teodoro Laino [tlaino] - University of Zurich
!> \date 04.2008
! *****************************************************************************
SUBROUTINE drotint (sepi,sepj,rij,dw,delta,anag,se_int_control,taper, error)
SUBROUTINE drotint (sepi,sepj,rij,dw,delta,anag,se_int_control,se_taper, error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rij
REAL(dp), DIMENSION(3, 2025), &
@ -376,16 +376,16 @@ CONTAINS
REAL(dp), INTENT(IN) :: delta
LOGICAL, INTENT(IN) :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'drotint', &
routineP = moduleN//':'//routineN
IF (anag) THEN
CALL rotint_ana(sepi,sepj,rij,dw=dw, se_int_control=se_int_control,taper=taper,error=error)
CALL rotint_ana(sepi,sepj,rij,dw=dw, se_int_control=se_int_control,se_taper=se_taper,error=error)
ELSE
CALL drotint_num(sepi,sepj,rij,dw,delta,se_int_control=se_int_control,taper=taper,error=error)
CALL drotint_num(sepi,sepj,rij,dw,delta,se_int_control=se_int_control,se_taper=se_taper,error=error)
END IF
END SUBROUTINE drotint
@ -395,7 +395,7 @@ CONTAINS
!> \author Teodoro Laino [tlaino] - University of Zurich
!> \date 04.2008
! *****************************************************************************
SUBROUTINE drotnuc (sepi,sepj,rij,de1b,de2a,itype,delta,anag,se_int_control,taper,error)
SUBROUTINE drotnuc (sepi,sepj,rij,de1b,de2a,itype,delta,anag,se_int_control,se_taper,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rij
REAL(dp), DIMENSION(3, 45), &
@ -404,7 +404,7 @@ CONTAINS
REAL(dp), INTENT(IN) :: delta
LOGICAL, INTENT(IN) :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'drotnuc', &
@ -412,10 +412,10 @@ CONTAINS
IF (anag) THEN
CALL rotnuc_ana (sepi, sepj, rij, de1b=de1b, de2a=de2a, itype=itype,&
se_int_control=se_int_control,taper=taper,error=error)
se_int_control=se_int_control,se_taper=se_taper,error=error)
ELSE
CALL drotnuc_num (sepi, sepj, rij, de1b=de1b, de2a=de2a, itype=itype,&
delta=delta,se_int_control=se_int_control,taper=taper,error=error)
delta=delta,se_int_control=se_int_control,se_taper=se_taper,error=error)
END IF
END SUBROUTINE drotnuc
@ -425,7 +425,7 @@ CONTAINS
!> \author Teodoro Laino [tlaino] - University of Zurich
!> \date 04.2008
! *****************************************************************************
SUBROUTINE dcorecore (sepi,sepj,rij,denuc,itype,delta,anag,se_int_control,taper,error)
SUBROUTINE dcorecore (sepi,sepj,rij,denuc,itype,delta,anag,se_int_control,se_taper,error)
TYPE(semi_empirical_type), POINTER :: sepi, sepj
REAL(dp), DIMENSION(3), INTENT(IN) :: rij
REAL(dp), DIMENSION(3), INTENT(OUT) :: denuc
@ -433,7 +433,7 @@ CONTAINS
REAL(dp), INTENT(IN) :: delta
LOGICAL, INTENT(IN) :: anag
TYPE(se_int_control_type), INTENT(IN) :: se_int_control
TYPE(taper_type), POINTER :: taper
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'dcorecore', &
@ -441,10 +441,10 @@ CONTAINS
IF (anag) THEN
CALL corecore_ana (sepi,sepj,rij,denuc=denuc,itype=itype,se_int_control=se_int_control,&
taper=taper, error=error)
se_taper=se_taper, error=error)
ELSE
CALL dcorecore_num (sepi,sepj,rij,denuc=denuc,delta=delta,itype=itype,&
se_int_control=se_int_control,taper=taper, error=error)
se_int_control=se_int_control,se_taper=se_taper, error=error)
END IF
END SUBROUTINE dcorecore

View file

@ -15,14 +15,10 @@ MODULE semi_empirical_types
cp_print_key_should_output,&
cp_print_key_unit_nr
USE f77_blas
USE input_constants, ONLY: do_method_am1,&
do_method_mndo,&
do_method_mndod,&
do_method_pdg,&
do_method_pm3,&
do_method_pm6,&
do_method_pnnl,&
do_method_rm1
USE input_constants, ONLY: &
do_method_am1, do_method_mndo, do_method_mndod, do_method_pdg, &
do_method_pm3, do_method_pm6, do_method_pnnl, do_method_rm1, &
do_se_IS_kdso_d
USE input_section_types, ONLY: section_vals_type
USE kinds, ONLY: default_string_length,&
dp
@ -33,6 +29,9 @@ MODULE semi_empirical_types
semi_empirical_expns3_release
USE semi_empirical_mpole_types, ONLY: semi_empirical_mpole_p_release,&
semi_empirical_mpole_p_type
USE taper_types, ONLY: taper_create,&
taper_release,&
taper_type
USE termination, ONLY: stop_program
#include "cp_common_uses.h"
@ -121,11 +120,25 @@ MODULE semi_empirical_types
INTEGER :: integral_screening
END TYPE se_int_control_type
! *****************************************************************************
!> \brief Store the value of the tapering function and possibly its derivative
!> for screened integrals
! *****************************************************************************
TYPE se_int_screen_type
! Store the value of the tapering function and possibly its derivative
REAL(KIND=dp) :: ft, dft
END TYPE se_int_screen_type
! *****************************************************************************
!> \brief Taper type use in semi-empirical calculations
! *****************************************************************************
TYPE se_taper_type
TYPE(taper_type), POINTER :: taper
TYPE(taper_type), POINTER :: taper_cou
TYPE(taper_type), POINTER :: taper_exc
! This taper is for KDSO-D integrals
TYPE(taper_type), POINTER :: taper_add
END TYPE se_taper_type
PUBLIC :: semi_empirical_type,&
semi_empirical_create,&
semi_empirical_release,&
@ -137,7 +150,10 @@ MODULE semi_empirical_types
write_se_param,&
se_int_control_type,&
setup_se_int_control_type,&
se_int_screen_type
se_int_screen_type,&
se_taper_type,&
se_taper_release,&
se_taper_create
CONTAINS
@ -505,6 +521,66 @@ CONTAINS
se_int_control%integral_screening = integral_screening
END SUBROUTINE setup_se_int_control_type
! *****************************************************************************
!> \brief Creates the taper type used in SE calculations
!> \author Teodoro Laino [tlaino] - 03.2009
! *****************************************************************************
SUBROUTINE se_taper_create(se_taper, integral_screening, taper_cou,&
range_cou, taper_exc, range_exc, taper_scr, range_scr, error)
TYPE(se_taper_type), POINTER :: se_taper
INTEGER, INTENT(IN) :: integral_screening
REAL(KIND=dp), INTENT(IN) :: taper_cou, range_cou, &
taper_exc, range_exc, &
taper_scr, range_scr
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'se_taper_create', &
routineP = moduleN//':'//routineN
INTEGER :: stat
LOGICAL :: failure
failure = .FALSE.
CPPrecondition(.NOT.ASSOCIATED(se_taper),cp_failure_level,routineP,error,failure)
ALLOCATE(se_taper, stat=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
NULLIFY(se_taper%taper)
NULLIFY(se_taper%taper_cou)
NULLIFY(se_taper%taper_exc)
NULLIFY(se_taper%taper_add)
! Create the sub-typo taper
CALL taper_create(se_taper%taper_cou, taper_cou, range_cou, error)
CALL taper_create(se_taper%taper_exc, taper_exc, range_exc, error)
IF (integral_screening==do_se_IS_kdso_d) THEN
CALL taper_create(se_taper%taper_add, taper_scr, range_scr, error)
END IF
END SUBROUTINE se_taper_create
! *****************************************************************************
!> \brief Releases the taper type used in SE calculations
!> \author Teodoro Laino [tlaino] - 03.2009
! *****************************************************************************
SUBROUTINE se_taper_release(se_taper, error)
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'se_taper_release', &
routineP = moduleN//':'//routineN
INTEGER :: stat
LOGICAL :: failure
failure = .FALSE.
IF (ASSOCIATED(se_taper)) THEN
CALL taper_release(se_taper%taper_cou, error=error)
CALL taper_release(se_taper%taper_exc, error=error)
CALL taper_release(se_taper%taper_add, error=error)
DEALLOCATE(se_taper, stat=stat)
CPPrecondition(stat==0,cp_failure_level,routineP,error,failure)
END IF
END SUBROUTINE se_taper_release
! *****************************************************************************
!> \brief Writes the semi-empirical type
!> \par History

View file

@ -29,7 +29,8 @@ MODULE semi_empirical_utils
pm6_default_parameter,&
pnnl_default_parameter,&
rm1_default_parameter
USE semi_empirical_types, ONLY: semi_empirical_type
USE semi_empirical_types, ONLY: se_taper_type,&
semi_empirical_type
USE termination, ONLY: stop_program
#include "cp_common_uses.h"
@ -39,9 +40,56 @@ MODULE semi_empirical_utils
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_utils'
PUBLIC :: init_se_param, se_param_set_default, get_se_type
PUBLIC :: init_se_param, se_param_set_default, get_se_type,&
initialize_se_taper, finalize_se_taper
CONTAINS
! *****************************************************************************
!> \brief Initializes the semi-empirical taper for a chunk calculation
!> \author Teodoro Laino [tlaino] - 03.2009
! *****************************************************************************
SUBROUTINE initialize_se_taper(se_taper, coulomb, exchange, error)
TYPE(se_taper_type), POINTER :: se_taper
LOGICAL, INTENT(IN), OPTIONAL :: coulomb, exchange
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'initialize_se_taper', &
routineP = moduleN//':'//routineN
LOGICAL :: check, failure, l_coulomb, &
l_exchange
failure = .FALSE.
check = PRESENT(coulomb).NEQV.PRESENT(exchange)
CPPrecondition(check,cp_failure_level,routineP,error,failure)
check = .NOT.ASSOCIATED(se_taper%taper)
CPPrecondition(check,cp_failure_level,routineP,error,failure)
l_coulomb = .FALSE.
l_exchange = .FALSE.
IF (PRESENT(coulomb)) l_coulomb = coulomb
IF (PRESENT(exchange)) l_exchange = exchange
IF (l_coulomb) se_taper%taper => se_taper%taper_cou
IF (l_exchange) se_taper%taper => se_taper%taper_exc
END SUBROUTINE initialize_se_taper
! *****************************************************************************
!> \brief Finalizes the semi-empirical taper for a chunk calculation
!> \author Teodoro Laino [tlaino] - 03.2009
! *****************************************************************************
SUBROUTINE finalize_se_taper(se_taper, error)
TYPE(se_taper_type), POINTER :: se_taper
TYPE(cp_error_type), INTENT(inout) :: error
CHARACTER(len=*), PARAMETER :: routineN = 'finalize_se_taper', &
routineP = moduleN//':'//routineN
LOGICAL :: check, failure
failure = .FALSE.
check = ASSOCIATED(se_taper%taper)
CPPrecondition(check,cp_failure_level,routineP,error,failure)
NULLIFY(se_taper%taper)
END SUBROUTINE finalize_se_taper
! *****************************************************************************
!> \brief Initialize semi_empirical type

View file

@ -84,7 +84,7 @@ CONTAINS
! *****************************************************************************
!> \brief Taper functions
! *****************************************************************************
FUNCTION taper_eval (taper,rij, error) RESULT(ft)
FUNCTION taper_eval (taper, rij, error) RESULT(ft)
TYPE(taper_type), POINTER :: taper
REAL(KIND=dp), INTENT(IN) :: rij
TYPE(cp_error_type), INTENT(inout) :: error

View file

@ -221,3 +221,5 @@ stp6.inp
zn_1w_debug.inp
# more resets after change on GUESS
scan_zn_1w_cp2k.inp
# numerics
md.inp

View file

@ -15,9 +15,14 @@
&QS
METHOD PM6
&SE
RC_INTERACTION [angstrom] 50.0
RC_COULOMB [angstrom] 50.0
RC_RANGE [angstrom] 0.0
&COULOMB
CUTOFF [angstrom] 50.0
RC_RANGE [angstrom] 0.0
&END
&EXCHANGE
CUTOFF [angstrom] 50.0
RC_RANGE [angstrom] 0.0
&END
&END
&END QS
&SCF

View file

@ -7,7 +7,13 @@
&QS
METHOD PM3
&SE
RC_RANGE [bohr] 1.0
&COULOMB
RC_RANGE [bohr] 1.0
&END
&EXCHANGE
RC_TAPER [bohr] 24.0
RC_RANGE [bohr] 1.0
&END
&END SE
&END QS
&SCF

View file

@ -6,7 +6,12 @@
MULTIPLICITY 1
&QS
&SE
RC_RANGE 0.1
&COULOMB
RC_RANGE [bohr] 0.1
&END
&EXCHANGE
RC_RANGE [bohr] 0.1
&END
&END
METHOD PM6
EXTRAPOLATION ASPC

View file

@ -15,9 +15,14 @@
&QS
METHOD PM6
&SE
RC_INTERACTION [angstrom] 50.0
RC_COULOMB [angstrom] 50.0
RC_RANGE [angstrom] 0.0
&COULOMB
CUTOFF [angstrom] 50.0
RC_RANGE [angstrom] 0.0
&END
&EXCHANGE
CUTOFF [angstrom] 50.0
RC_RANGE [angstrom] 0.0
&END
&END
&END QS
&SCF

View file

@ -6,7 +6,13 @@
MULTIPLICITY 1
&QS
&SE
RC_RANGE 10.
&COULOMB
RC_RANGE [bohr] 10.0
&END
&EXCHANGE
RC_TAPER [bohr] 24.0
RC_RANGE 10.0
&END
&END
METHOD PM6
EXTRAPOLATION ASPC

View file

@ -7,8 +7,14 @@
METHOD MNDO
&SE
PERIODIC EWALD
RC_COULOMB [angstrom] 5.0
RC_RANGE 0.01
&COULOMB
CUTOFF [angstrom] 5.0
RC_RANGE 0.01
&END
&EXCHANGE
RC_TAPER [angstrom] 5.0
RC_RANGE 0.01
&END
&PRINT
&EWALD_INFO
&END

View file

@ -7,8 +7,14 @@
METHOD MNDO
&SE
PERIODIC EWALD
RC_COULOMB [angstrom] 15.0
RC_RANGE 0.001
&COULOMB
CUTOFF [angstrom] 15.0
RC_RANGE 0.001
&END
&EXCHANGE
RC_TAPER [angstrom] 15.0
RC_RANGE 0.001
&END
&PRINT
&EWALD_INFO
&END

View file

@ -7,8 +7,14 @@
METHOD MNDO
&SE
PERIODIC EWALD
RC_COULOMB [angstrom] 11.1
RC_RANGE 0.01
&COULOMB
CUTOFF [angstrom] 11.1
RC_RANGE 0.01
&END
&EXCHANGE
RC_TAPER [angstrom] 11.1
RC_RANGE 0.01
&END
&PRINT
&EWALD_INFO
&END

View file

@ -7,9 +7,15 @@
METHOD PM6
&SE
ANALYTICAL_GRADIENTS T
RC_INTERACTION [angstrom] 4.9325
RC_COULOMB [angstrom] 12.0
RC_RANGE 1.0
&COULOMB
CUTOFF [angstrom] 12.0
RC_RANGE 1.0
&END
&EXCHANGE
CUTOFF [angstrom] 4.9325
RC_TAPER [angstrom] 12.0
RC_RANGE 1.0
&END
&END
&END QS
&SCF

View file

@ -6,9 +6,19 @@
METHOD MNDO
&SE
INTEGRAL_SCREENING KDSO-D
RC_COULOMB [angstrom] 6.0
RC_INTERACTION [angstrom] 6.0
RC_RANGE [angstrom] 0.1
&COULOMB
CUTOFF [angstrom] 6.0
RC_RANGE [angstrom] 0.1
&END
&EXCHANGE
CUTOFF [angstrom] 6.0
RC_TAPER [angstrom] 6.0
RC_RANGE [angstrom] 0.1
&END
&SCREENING
RC_TAPER [angstrom] 6.0
RC_RANGE [angstrom] 0.1
&END
&PRINT
&EWALD_INFO
&END

View file

@ -7,8 +7,14 @@
METHOD MNDO
&SE
PERIODIC EWALD_R3
RC_COULOMB [angstrom] 15.0
RC_RANGE 0.001
&COULOMB
CUTOFF [angstrom] 15.0
RC_RANGE 0.001
&END
&EXCHANGE
RC_TAPER [angstrom] 15.0
RC_RANGE 0.001
&END
&PRINT
&EWALD_INFO
&END

View file

@ -6,8 +6,14 @@
METHOD MNDO
&SE
PERIODIC EWALD
RC_COULOMB [angstrom] 15.0
RC_RANGE 0.001
&COULOMB
CUTOFF [angstrom] 15.0
RC_RANGE 0.001
&END
&EXCHANGE
RC_TAPER [angstrom] 15.0
RC_RANGE 0.001
&END
&PRINT
&EWALD_INFO
&END

View file

@ -6,9 +6,19 @@
METHOD MNDO
&SE
INTEGRAL_SCREENING KDSO-D
RC_COULOMB [angstrom] 6.0
RC_INTERACTION [angstrom] 50.0
RC_RANGE [angstrom] 0.1
&COULOMB
CUTOFF [angstrom] 6.0
RC_RANGE [angstrom] 0.1
&END
&EXCHANGE
CUTOFF [angstrom] 50.0
RC_TAPER [angstrom] 6.0
RC_RANGE [angstrom] 0.1
&END
&SCREENING
RC_TAPER [angstrom] 6.0
RC_RANGE [angstrom] 0.1
&END
&PRINT
&EWALD_INFO
&END

View file

@ -6,9 +6,14 @@
METHOD MNDO
&SE
INTEGRAL_SCREENING KDSO
RC_COULOMB [angstrom] 50.0
RC_INTERACTION [angstrom] 50.0
RC_RANGE [angstrom] 0.1
&COULOMB
CUTOFF [angstrom] 50.0
RC_RANGE [angstrom] 0.1
&END
&EXCHANGE
CUTOFF [angstrom] 50.0
RC_RANGE [angstrom] 0.1
&END
&PRINT
&EWALD_INFO
&END

View file

@ -7,9 +7,15 @@
EPS_DEFAULT 1.0E-10
&SE
ANALYTICAL_GRADIENTS T
RC_INTERACTION [angstrom] 4.9325
RC_COULOMB [angstrom] 12.0
RC_RANGE 1.0
&COULOMB
CUTOFF [angstrom] 12.0
RC_RANGE 1.0
&END
&EXCHANGE
CUTOFF [angstrom] 4.9325
RC_TAPER [angstrom] 12.0
RC_RANGE 1.0
&END
&END
&END QS
&SCF

View file

@ -6,8 +6,14 @@
METHOD MNDO
&SE
PERIODIC EWALD
RC_COULOMB [angstrom] 15.0
RC_RANGE 0.001
&COULOMB
CUTOFF [angstrom] 15.0
RC_RANGE 0.001
&END
&EXCHANGE
RC_TAPER [angstrom] 15.0
RC_RANGE 0.001
&END
&PRINT
&EWALD_INFO OFF
&END

View file

@ -14,7 +14,9 @@
&QS
METHOD PM6
&SE
RC_COULOMB [angstrom] 4.175
&COULOMB
CUTOFF [angstrom] 4.175
&END
&END
&END QS
&SCF

View file

@ -26,3 +26,5 @@ H2O-MNDO.inp
O2.inp
# bug fix
LM-KSDO-D.inp
# numerics
LM-KSDO-D.inp