comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-27 13:45:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Efield implementation for xTB (#3883)

Browse source
This commit is contained in:
Juerg Hutter 2025-01-22 14:18:16 +01:00 • committed by GitHub
parent 9a7de708c1
commit abdbd495f7
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
40 changed files with 1642 additions and 81 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/CMakeLists.txt
View file

@ -848,6 +848,7 @@ list(
xas_tdp_utils.F
xas_tp_scf.F
xray_diffraction.F
xtb_qresp.F
xtb_coulomb.F
xtb_eeq.F
xtb_ehess.F

1
src/cp_control_types.F
View file

@ -169,6 +169,7 @@ MODULE cp_control_types
LOGICAL :: coulomb_lr = .FALSE.
LOGICAL :: tb3_interaction = .FALSE.
LOGICAL :: check_atomic_charges = .FALSE.
LOGICAL :: var_dipole = .FALSE.
!
REAL(KIND=dp) :: xb_radius = 0.0_dp
REAL(KIND=dp) :: coulomb_sr_cut = 0.0_dp

3
src/cp_control_utils.F
View file

@ -1455,6 +1455,8 @@ CONTAINS
qs_control%xtb_control%coulomb_lr = .FALSE.
qs_control%xtb_control%tb3_interaction = .FALSE.
qs_control%xtb_control%check_atomic_charges = .FALSE.
CALL section_vals_val_get(xtb_section, "VARIATIONAL_DIPOLE", &
l_val=qs_control%xtb_control%var_dipole)
CASE (1)
! For debugging purposes
CALL section_vals_val_get(xtb_section, "COULOMB_INTERACTION", &
@ -1466,6 +1468,7 @@ CONTAINS
! Check for bad atomic charges
CALL section_vals_val_get(xtb_section, "CHECK_ATOMIC_CHARGES", &
l_val=qs_control%xtb_control%check_atomic_charges)
qs_control%xtb_control%var_dipole = .FALSE.
CASE (2)
CPABORT("gfn2-xtb tbd")
CASE DEFAULT

56
src/eeq_method.F
View file

@ -118,33 +118,43 @@ CONTAINS
!> \param qs_env ...
!> \param iounit ...
!> \param print_level ...
!> \param ext ...
! **************************************************************************************************
SUBROUTINE eeq_print(qs_env, iounit, print_level)
SUBROUTINE eeq_print(qs_env, iounit, print_level, ext)
TYPE(qs_environment_type), POINTER :: qs_env
INTEGER, INTENT(IN) :: iounit, print_level
LOGICAL, INTENT(IN) :: ext
CHARACTER(LEN=2) :: element_symbol
INTEGER :: enshift_type, iatom, ikind, natom
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charges
REAL(KIND=dp), DIMENSION(:), POINTER :: charges
TYPE(cell_type), POINTER :: cell
TYPE(eeq_solver_type) :: eeq_sparam
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
MARK_USED(print_level)
CALL get_qs_env(qs_env, natom=natom, particle_set=particle_set)
ALLOCATE (charges(natom))
CALL get_qs_env(qs_env, cell=cell)
! enforce en shift method 1 (original/molecular)
! method 2 from paper on PBC seems not to work
enshift_type = 1
!IF (ALL(cell%perd == 0)) enshift_type = 1
CALL eeq_charges(qs_env, charges, eeq_sparam, 2, enshift_type)
CALL get_qs_env(qs_env, natom=natom, particle_set=particle_set, cell=cell)
IF (ext) THEN
NULLIFY (charges)
CALL get_qs_env(qs_env, eeq=charges)
CPASSERT(ASSOCIATED(charges))
enshift_type = 0
ELSE
ALLOCATE (charges(natom))
! enforce en shift method 1 (original/molecular)
! method 2 from paper on PBC seems not to work
enshift_type = 1
!IF (ALL(cell%perd == 0)) enshift_type = 1
CALL eeq_charges(qs_env, charges, eeq_sparam, 2, enshift_type)
END IF
IF (iounit > 0) THEN
IF (enshift_type == 1) THEN
IF (enshift_type == 0) THEN
WRITE (UNIT=iounit, FMT="(/,T2,A)") "EEQ Charges (External)"
ELSE IF (enshift_type == 1) THEN
WRITE (UNIT=iounit, FMT="(/,T2,A)") "EEQ Charges (Parametrization 2019 (Molecules))"
ELSE IF (enshift_type == 2) THEN
WRITE (UNIT=iounit, FMT="(/,T2,A)") "EEQ Charges (Parametrization 2019 (Crystals))"
@ -164,7 +174,7 @@ CONTAINS
END IF
DEALLOCATE (charges)
IF (.NOT. ext) DEALLOCATE (charges)
END SUBROUTINE eeq_print
@ -344,7 +354,7 @@ CONTAINS
LOGICAL :: do_ewald, use_virial
LOGICAL, ALLOCATABLE, DIMENSION(:) :: default_present
REAL(KIND=dp) :: ala, alb, alpha, cn, ctot, dr, dr2, drk, &
elag, esg, gam2, gama, grc, kappa, &
elag, esg, fe, gam2, gama, grc, kappa, &
qlam, qq, qq1, qq2, rcut, scn, sgamma, &
subcells, totalcharge, xi
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: c_radius, cnumbers, gam, qlag
@ -530,6 +540,8 @@ CONTAINS
dr2 = SUM(rij**2)
dr = SQRT(dr2)
IF (dr > rcut .OR. dr < 1.E-6_dp) CYCLE
fe = 1.0_dp
IF (iatom == jatom) fe = 0.5_dp
IF (response_only) THEN
qq = -qlag(iatom)*charges(jatom)
ELSE
@ -550,13 +562,13 @@ CONTAINS
gradient(:, iatom) = gradient(:, iatom) + fdik(:)
gradient(:, jatom) = gradient(:, jatom) - fdik(:)
IF (use_virial) THEN
CALL virial_pair_force(stress, -1._dp, fdik, rij)
CALL virial_pair_force(stress, -fe, fdik, rij)
END IF
fdik(:) = qq2*grc*rij(:)/dr
gradient(:, iatom) = gradient(:, iatom) - fdik(:)
gradient(:, jatom) = gradient(:, jatom) + fdik(:)
IF (use_virial) THEN
CALL virial_pair_force(stress, 1._dp, fdik, rij)
CALL virial_pair_force(stress, fe, fdik, rij)
END IF
END DO
CALL neighbor_list_iterator_release(nl_iterator)
@ -1479,7 +1491,6 @@ CONTAINS
CALL para_env%sum(potential(1:na))
CALL cp_fm_matvec(fm_mat, charges, potential, alpha=1.0_dp, beta=1.0_dp)
eeqn = 0.5_dp*SUM(charges(1:na)*potential(1:na)) + SUM(charges(1:na)*rhs(1:na))
! potential(1:na) = potential(1:na) + rhs(1:na)
potential(1:ns) = potential(1:ns) + rhs(1:ns)
ALLOCATE (mvec(ns))
CALL cp_fm_matvec(mmat, potential, mvec, alpha=-1.0_dp, beta=0.0_dp)
@ -1550,7 +1561,7 @@ CONTAINS
END DO
ef_energy = -0.25_dp*omega/twopi*ef_energy
ELSE
ef_energy = -SUM(fpolvec(:)*qi(:))
ef_energy = SUM(fpolvec(:)*qi(:))
END IF
ELSE IF (dft_control%apply_efield) THEN
@ -1616,7 +1627,7 @@ CONTAINS
DO idir = 1, 3
kvec(:) = twopi*cell%h_inv(idir, :)
kr = SUM(kvec(:)*ria(:))
efr(ia) = efr(ia) - kr*fpolvec(idir)
efr(ia) = efr(ia) + kr*fpolvec(idir)
END DO
END DO
END IF
@ -1756,6 +1767,7 @@ CONTAINS
LOGICAL :: dfield, use_virial
REAL(KIND=dp) :: q
REAL(KIND=dp), DIMENSION(3) :: fa, fieldpol, ria
REAL(KIND=dp), DIMENSION(3, 3) :: pve
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(dft_control_type), POINTER :: dft_control
@ -1786,23 +1798,25 @@ CONTAINS
CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
CALL get_qs_env(qs_env=qs_env, force=force)
pve = 0.0_dp
DO ikind = 1, nkind
force(ikind)%efield = 0.0_dp
DO ia = 1, local_particles%n_el(ikind)
iatom = local_particles%list(ikind)%array(ia)
q = charges(iatom) - qlag(iatom)
fa(1:3) = -q*fieldpol(1:3)
fa(1:3) = q*fieldpol(1:3)
atom_a = atom_of_kind(iatom)
force(ikind)%efield(1:3, atom_a) = fa
IF (use_virial) THEN
ria = particle_set(ia)%r
ria = pbc(ria, cell)
CALL virial_pair_force(virial%pv_virial, -0.5_dp, fa, ria)
CALL virial_pair_force(virial%pv_virial, -0.5_dp, ria, fa)
CALL virial_pair_force(pve, -0.5_dp, fa, ria)
CALL virial_pair_force(pve, -0.5_dp, ria, fa)
END IF
END DO
CALL para_env%sum(force(ikind)%efield)
END DO
virial%pv_virial = virial%pv_virial + pve
END SUBROUTINE eeq_efield_force_periodic

18
src/efield_tb_methods.F
View file

@ -345,7 +345,7 @@ CONTAINS
END DO
END DO
qi = AIMAG(LOG(zi))
energy%efield = -SUM(fpolvec(:)*qi(:))
energy%efield = SUM(fpolvec(:)*qi(:))
IF (.NOT. just_energy) THEN
CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
@ -364,7 +364,7 @@ CONTAINS
CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
IF (para_env%mepos == 0) THEN
DO ia = 1, natom
charge = -mcharge(ia)
charge = mcharge(ia)
iatom = atom_of_kind(ia)
ikind = kind_of(ia)
force(ikind)%efield(:, iatom) = fieldpol(:)*charge
@ -409,9 +409,9 @@ CONTAINS
CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
row=irow, col=icol, block=ks_block, found=found)
CPASSERT(found)
ks_block = ks_block + 0.5_dp*fdir*s_block
ks_block = ks_block - 0.5_dp*fdir*s_block
END DO
IF (calculate_forces) THEN
IF (calculate_forces .AND. irow /= icol) THEN
ikind = kind_of(irow)
jkind = kind_of(icol)
atom_a = atom_of_kind(irow)
@ -425,7 +425,7 @@ CONTAINS
CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, 1)%matrix, &
row=irow, col=icol, BLOCK=ds_block, found=found)
CPASSERT(found)
fij(idir) = fij(idir) + SUM(p_block*ds_block)
fij(idir) = fij(idir) - SUM(p_block*ds_block)
END DO
END DO
force(ikind)%efield(1:3, atom_a) = force(ikind)%efield(1:3, atom_a) + fdir*fij(1:3)
@ -471,10 +471,10 @@ CONTAINS
kvec(:) = twopi*cell%h_inv(idir, :)
dd = SUM(kvec(:)*ria(:))
zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
fdir = fdir - fpolvec(idir)*AIMAG(LOG(zdeta))
dd = SUM(kvec(:)*rib(:))
zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
fdir = fdir + fpolvec(idir)*AIMAG(LOG(zdeta))
fdir = fdir - fpolvec(idir)*AIMAG(LOG(zdeta))
END DO
fi = 1.0_dp
IF (iatom == jatom) fi = 0.5_dp
@ -535,7 +535,7 @@ CONTAINS
CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
row=irow, col=icol, block=ks_block, found=found)
CPASSERT(found)
ks_block = ks_block + 0.5_dp*fdir*s_block
ks_block = ks_block - 0.5_dp*fdir*s_block
END DO
IF (calculate_forces) THEN
atom_a = atom_of_kind(iatom)
@ -549,7 +549,7 @@ CONTAINS
CALL dbcsr_get_block_p(matrix=matrix_s(idir + 1, ic)%matrix, &
row=irow, col=icol, BLOCK=ds_block, found=found)
CPASSERT(found)
fij(idir) = fij(idir) + SUM(p_block*ds_block)
fij(idir) = fij(idir) - SUM(p_block*ds_block)
END DO
END DO
IF (irow == iatom) fij = -fij

6
src/input_cp2k_tb.F
View file

@ -216,6 +216,12 @@ CONTAINS
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="VARIATIONAL_DIPOLE", &
description="gfn0-xTB use dipole definition from energy derivative.", &
usage="VARIATIONAL_DIPOLE T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="EPS_PAIRPOTENTIAL", &
description="Accuracy for the repulsive pair potential.", &
usage="EPS_PAIRPOTENTIAL 1.0E-8", default_r_val=1.0e-10_dp)

23
src/qs_environment_types.F
View file

@ -308,6 +308,8 @@ MODULE qs_environment_types
TYPE(pw_r3d_rs_type), POINTER :: spin_embed_pot => NULL()
! Polarizability tensor
TYPE(polar_env_type), POINTER :: polar_env => NULL()
! EEQ charges
REAL(KIND=dp), DIMENSION(:), POINTER :: eeq => NULL()
! Resp charges
REAL(KIND=dp), DIMENSION(:), POINTER :: rhs => NULL()
REAL(KIND=dp) :: total_zeff_corr = -1.0_dp, surface_dipole_moment = -1.0_dp
@ -475,6 +477,7 @@ CONTAINS
!> \param spin_embed_pot ...
!> \param polar_env ...
!> \param mos_last_converged ... [SGh]
!> \param eeq ...
!> \param rhs ...
!> \date 23.01.2002
!> \author MK
@ -506,7 +509,8 @@ CONTAINS
rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, &
WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, &
s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, &
gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, &
eeq, rhs)
TYPE(qs_environment_type), INTENT(IN) :: qs_env
TYPE(atomic_kind_type), DIMENSION(:), OPTIONAL, &
POINTER :: atomic_kind_set
@ -639,7 +643,7 @@ CONTAINS
TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: embed_pot, spin_embed_pot
TYPE(polar_env_type), OPTIONAL, POINTER :: polar_env
TYPE(mo_set_type), DIMENSION(:), OPTIONAL, POINTER :: mos_last_converged
REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: rhs
REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eeq, rhs
TYPE(rho0_mpole_type), POINTER :: rho0_m
@ -706,6 +710,9 @@ CONTAINS
! Polarisability tensor
IF (PRESENT(polar_env)) polar_env => qs_env%polar_env
! EEQ charges
IF (PRESENT(eeq)) eeq => qs_env%eeq
! Resp charges
IF (PRESENT(rhs)) rhs => qs_env%rhs
@ -1023,6 +1030,7 @@ CONTAINS
!> \param spin_embed_pot ...
!> \param polar_env ...
!> \param mos_last_converged ... [SGh]
!> \param eeq ...
!> \param rhs ...
!> \date 23.01.2002
!> \author MK
@ -1043,7 +1051,7 @@ CONTAINS
do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, &
harris_env, gcp_env, mp2_env, bs_env, kg_env, force, &
kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, &
spin_embed_pot, polar_env, mos_last_converged, rhs)
spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
TYPE(qs_environment_type), INTENT(INOUT) :: qs_env
TYPE(cell_type), OPTIONAL, POINTER :: super_cell
@ -1119,7 +1127,7 @@ CONTAINS
TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: embed_pot, spin_embed_pot
TYPE(polar_env_type), OPTIONAL, POINTER :: polar_env
TYPE(mo_set_type), DIMENSION(:), OPTIONAL, POINTER :: mos_last_converged
REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: rhs
REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: eeq, rhs
TYPE(qs_subsys_type), POINTER :: subsys
@ -1343,6 +1351,9 @@ CONTAINS
IF (PRESENT(WannierCentres)) qs_env%WannierCentres => WannierCentres
IF (PRESENT(kpoints)) CALL set_ks_env(qs_env%ks_env, kpoints=kpoints)
! EEQ charges
IF (PRESENT(eeq)) qs_env%eeq => eeq
! Resp charges
IF (PRESENT(rhs)) qs_env%rhs => rhs
@ -1597,6 +1608,8 @@ CONTAINS
END DO
DEALLOCATE (qs_env%WannierCentres)
END IF
! EEQ charges
IF (ASSOCIATED(qs_env%eeq)) DEALLOCATE (qs_env%eeq)
! Resp charges
IF (ASSOCIATED(qs_env%rhs)) DEALLOCATE (qs_env%rhs)
@ -1814,6 +1827,8 @@ CONTAINS
END DO
DEALLOCATE (qs_env%WannierCentres)
END IF
! EEQ charges
IF (ASSOCIATED(qs_env%eeq)) DEALLOCATE (qs_env%eeq)
! Resp charges
IF (ASSOCIATED(qs_env%rhs)) DEALLOCATE (qs_env%rhs)

2
src/qs_scf_post_gpw.F
View file

@ -2674,7 +2674,7 @@ CONTAINS
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%EEQ_CHARGES", extension=".eeq", log_filename=.FALSE.)
print_level = 1
CALL eeq_print(qs_env, unit_nr, print_level)
CALL eeq_print(qs_env, unit_nr, print_level, ext=.FALSE.)
CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MULLIKEN")
END IF

2
src/qs_scf_post_se.F
View file

@ -566,7 +566,7 @@ CONTAINS
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%EEQ_CHARGES", &
extension=".eeq", log_filename=.FALSE.)
CALL eeq_print(qs_env, unit_nr, print_level=1)
CALL eeq_print(qs_env, unit_nr, print_level=1, ext=.FALSE.)
CALL cp_print_key_finished_output(unit_nr, logger, print_key)
END IF

58
src/qs_scf_post_tb.F
View file

@ -48,9 +48,7 @@ MODULE qs_scf_post_tb
USE cp_result_methods, ONLY: cp_results_erase,&
put_results
USE cp_result_types, ONLY: cp_result_type
USE eeq_input, ONLY: eeq_solver_type
USE eeq_method, ONLY: eeq_charges,&
eeq_print
USE eeq_method, ONLY: eeq_print
USE input_constants, ONLY: ot_precond_full_all
USE input_section_types, ONLY: section_get_ival,&
section_get_ivals,&
@ -138,6 +136,7 @@ MODULE qs_scf_post_tb
USE task_list_methods, ONLY: generate_qs_task_list
USE task_list_types, ONLY: allocate_task_list,&
task_list_type
USE xtb_qresp, ONLY: build_xtb_qresp
USE xtb_types, ONLY: get_xtb_atom_param,&
xtb_atom_type
#include "./base/base_uses.f90"
@ -173,17 +172,17 @@ CONTAINS
CHARACTER(LEN=6) :: ana
CHARACTER(LEN=default_string_length) :: aname
INTEGER :: after, enshift_type, handle, homo, iat, iatom, ikind, img, ispin, iw, nat, natom, &
INTEGER :: after, gfn_type, handle, homo, iat, iatom, ikind, img, ispin, iw, nat, natom, &
nkind, nlumo_stm, nlumos, nspins, print_level, unit_nr
LOGICAL :: do_cube, do_kpoints, explicit, has_homo, &
omit_headers, print_it, rebuild
LOGICAL :: do_cube, do_kpoints, explicit, gfn0, &
has_homo, omit_headers, print_it, &
rebuild, vdip
REAL(KIND=dp) :: zeff
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: echarge, mcharge
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mcharge
REAL(KIND=dp), DIMENSION(2, 2) :: homo_lumo
REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
REAL(KIND=dp), DIMENSION(:), POINTER :: echarge, mo_eigenvalues
REAL(KIND=dp), DIMENSION(:, :), POINTER :: charges
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: unoccupied_evals_stm
TYPE(cp_fm_type), DIMENSION(:), POINTER :: unoccupied_orbs_stm
TYPE(cp_fm_type), POINTER :: mo_coeff
@ -192,7 +191,6 @@ CONTAINS
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_p, matrix_s
TYPE(dbcsr_type), POINTER :: mo_coeff_deriv
TYPE(dft_control_type), POINTER :: dft_control
TYPE(eeq_solver_type) :: eeq_sparam
TYPE(kpoint_type), POINTER :: kpoints
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(mp_para_env_type), POINTER :: para_env
@ -213,10 +211,23 @@ CONTAINS
logger => cp_get_default_logger()
gfn0 = .FALSE.
vdip = .FALSE.
CALL get_qs_env(qs_env, dft_control=dft_control)
SELECT CASE (TRIM(tb_type))
CASE ("DFTB")
CASE ("xTB")
gfn_type = dft_control%qs_control%xtb_control%gfn_type
gfn0 = (gfn_type == 0)
vdip = dft_control%qs_control%xtb_control%var_dipole
CASE DEFAULT
CPABORT("unknown TB type")
END SELECT
CPASSERT(ASSOCIATED(qs_env))
NULLIFY (dft_control, rho, para_env, matrix_s, matrix_p)
NULLIFY (rho, para_env, matrix_s, matrix_p)
CALL get_qs_env(qs_env, scf_env=scf_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
dft_control=dft_control, rho=rho, natom=natom, para_env=para_env, &
rho=rho, natom=natom, para_env=para_env, &
particle_set=particle_set, do_kpoints=do_kpoints, matrix_s_kp=matrix_s)
nspins = dft_control%nspins
CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
@ -345,7 +356,7 @@ CONTAINS
IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
unit_nr = cp_print_key_unit_nr(logger, print_section, "EEQ_CHARGES", &
extension=".eeq", log_filename=.FALSE.)
CALL eeq_print(qs_env, unit_nr, print_level)
CALL eeq_print(qs_env, unit_nr, print_level, ext=gfn0)
CALL cp_print_key_finished_output(unit_nr, logger, print_key)
END IF
@ -382,18 +393,15 @@ CONTAINS
unit_nr = cp_print_key_unit_nr(logger, print_section, "MOMENTS", &
extension=".data", middle_name="tb_dipole", log_filename=.FALSE.)
moments_section => section_vals_get_subs_vals(print_section, "MOMENTS")
IF (tb_type == "xTB" .AND. dft_control%qs_control%xtb_control%gfn_type == 0) THEN
enshift_type = dft_control%qs_control%xtb_control%enshift_type
IF (enshift_type == 0) THEN
CALL get_qs_env(qs_env, cell=cell)
enshift_type = 1
IF (.NOT. ALL(cell%perd == 0)) enshift_type = 2
IF (gfn0) THEN
NULLIFY (echarge)
CALL get_qs_env(qs_env, eeq=echarge)
CPASSERT(ASSOCIATED(echarge))
IF (vdip) THEN
CALL build_xtb_qresp(qs_env, mcharge)
mcharge(1:natom) = echarge(1:natom) - mcharge(1:natom)
END IF
ALLOCATE (echarge(natom))
echarge = 0.0_dp
CALL eeq_charges(qs_env, echarge, eeq_sparam, 1, enshift_type)
CALL tb_dipole(qs_env, moments_section, unit_nr, echarge)
DEALLOCATE (echarge)
CALL tb_dipole(qs_env, moments_section, unit_nr, mcharge)
ELSE
CALL tb_dipole(qs_env, moments_section, unit_nr, mcharge)
END IF
@ -836,7 +844,7 @@ CONTAINS
IF (unit_nr > 0) THEN
WRITE (unit_nr, '(/,T2,A,T31,A50)') &
'TB_DIPOLE| Dipole type', ADJUSTR(TRIM(dipole_type))
WRITE (unit_nr, "(T2,A,T33,3F16.8)") "TB_DIPOLE| Reference Point [Bohr]", rcc
WRITE (unit_nr, "(T2,A,T30,3(1X,F16.8))") "TB_DIPOLE| Ref. Point [Bohr]", rcc
WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &
'TB_DIPOLE| Moment [a.u.]', dipole(1:3)
WRITE (unit_nr, '(T2,A,T30,3(1X,F16.8))') &

6
src/response_solver.F
View file

@ -2498,6 +2498,12 @@ CONTAINS
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Total Force ", fodeb
DEALLOCATE (ftot1, ftot2, ftot3)
END IF
IF (debug_stress .AND. use_virial) THEN
stdeb = fconv*(virial%pv_virial - sttot)
CALL para_env%sum(stdeb)
IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
'STRESS| Stress Response ', one_third_sum_diag(stdeb), det_3x3(stdeb)
END IF
IF (do_ex) THEN
CALL dbcsr_deallocate_matrix_set(mpa)

123
src/xtb_eeq.F
View file

@ -60,7 +60,7 @@ MODULE xtb_eeq
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_eeq'
PUBLIC :: xtb_eeq_calculation, xtb_eeq_forces
PUBLIC :: xtb_eeq_calculation, xtb_eeq_forces, xtb_eeq_lagrange
CONTAINS
@ -220,14 +220,17 @@ CONTAINS
!> \param qs_env ...
!> \param charges ...
!> \param dcharges ...
!> \param qlagrange ...
!> \param cnumbers ...
!> \param dcnum ...
!> \param eeq_sparam ...
! **************************************************************************************************
SUBROUTINE xtb_eeq_forces(qs_env, charges, dcharges, cnumbers, dcnum, eeq_sparam)
SUBROUTINE xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
TYPE(qs_environment_type), POINTER :: qs_env
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges, dcharges, cnumbers
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: charges, dcharges
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: qlagrange
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cnumbers
TYPE(dcnum_type), DIMENSION(:), INTENT(IN) :: dcnum
TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
@ -240,7 +243,7 @@ CONTAINS
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
LOGICAL :: defined, do_ewald, use_virial
REAL(KIND=dp) :: ala, alb, alpha, cn, ctot, dr, dr2, drk, &
elag, esg, gam2, gama, grc, kappa, &
elag, esg, fe, gam2, gama, grc, kappa, &
qlam, qq, qq1, qq2, rcut, scn, sgamma, &
totalcharge, xi
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: gam
@ -393,7 +396,7 @@ CONTAINS
IF (do_ewald) THEN
CALL get_qs_env(qs_env, sab_tbe=sab_tbe)
CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
rcut = 2.0_dp*rcut
rcut = 1.0_dp*rcut
CALL neighbor_list_iterator_create(nl_iterator, sab_tbe)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
@ -404,6 +407,8 @@ CONTAINS
dr2 = SUM(rij**2)
dr = SQRT(dr2)
IF (dr > rcut .OR. dr < 1.E-6_dp) CYCLE
fe = 1.0_dp
IF (iatom == jatom) fe = 0.5_dp
qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
gama = gab(ikind, jkind)
gam2 = gama*gama
@ -415,13 +420,13 @@ CONTAINS
force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + fdik(:)
force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) - fdik(:)
IF (use_virial) THEN
CALL virial_pair_force(virial%pv_virial, 1._dp, fdik, rij)
CALL virial_pair_force(virial%pv_virial, fe, fdik, rij)
END IF
fdik(:) = qq2*grc*rij(:)/dr
force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) - fdik(:)
force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) + fdik(:)
IF (use_virial) THEN
CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rij)
CALL virial_pair_force(virial%pv_virial, -fe, fdik, rij)
END IF
END DO
CALL neighbor_list_iterator_release(nl_iterator)
@ -482,10 +487,114 @@ CONTAINS
END IF
END IF
! return Lagrange multipliers
qlagrange(1:natom) = qlag(1:natom)
DEALLOCATE (gab, chrgx, dchia, qlag)
CALL timestop(handle)
END SUBROUTINE xtb_eeq_forces
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param dcharges ...
!> \param qlagrange ...
!> \param eeq_sparam ...
! **************************************************************************************************
SUBROUTINE xtb_eeq_lagrange(qs_env, dcharges, qlagrange, eeq_sparam)
TYPE(qs_environment_type), POINTER :: qs_env
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: dcharges
REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: qlagrange
TYPE(eeq_solver_type), INTENT(IN) :: eeq_sparam
CHARACTER(len=*), PARAMETER :: routineN = 'xtb_eeq_lagrange'
INTEGER :: handle, ikind, iunit, jkind, natom, nkind
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
LOGICAL :: defined, do_ewald
REAL(KIND=dp) :: ala, alb, elag, gama, qlam, totalcharge
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: gam
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gab
REAL(KIND=dp), DIMENSION(:), POINTER :: qlag
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(dft_control_type), POINTER :: dft_control
TYPE(ewald_environment_type), POINTER :: ewald_env
TYPE(ewald_pw_type), POINTER :: ewald_pw
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
TYPE(xtb_control_type), POINTER :: xtb_control
CALL timeset(routineN, handle)
iunit = cp_logger_get_default_unit_nr()
CALL get_qs_env(qs_env, &
qs_kind_set=qs_kind_set, &
atomic_kind_set=atomic_kind_set, &
particle_set=particle_set, &
cell=cell, &
dft_control=dft_control)
CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
xtb_control => dft_control%qs_control%xtb_control
totalcharge = dft_control%charge
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
atom_of_kind=atom_of_kind, kind_of=kind_of)
! gamma[a,b]
ALLOCATE (gab(nkind, nkind), gam(nkind))
gab = 0.0_dp
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
IF (.NOT. defined) CYCLE
CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
gam(ikind) = gama
DO jkind = 1, nkind
CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
IF (.NOT. defined) CYCLE
CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
!
gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
!
END DO
END DO
ALLOCATE (qlag(natom))
do_ewald = xtb_control%do_ewald
CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
IF (do_ewald) THEN
CALL get_qs_env(qs_env=qs_env, &
ewald_env=ewald_env, ewald_pw=ewald_pw)
CALL eeq_solver(qlag, qlam, elag, &
particle_set, kind_of, cell, -dcharges, gam, gab, &
para_env, blacs_env, dft_control, eeq_sparam, &
ewald=do_ewald, ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
ELSE
CALL eeq_solver(qlag, qlam, elag, &
particle_set, kind_of, cell, -dcharges, gam, gab, &
para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
END IF
! return Lagrange multipliers
qlagrange(1:natom) = qlag(1:natom)
DEALLOCATE (gab, qlag)
CALL timestop(handle)
END SUBROUTINE xtb_eeq_lagrange
END MODULE xtb_eeq

23
src/xtb_matrices.F
View file

@ -57,7 +57,8 @@ MODULE xtb_matrices
USE qs_dispersion_types, ONLY: qs_dispersion_type
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
qs_environment_type,&
set_qs_env
USE qs_force_types, ONLY: qs_force_type
USE qs_integral_utils, ONLY: basis_set_list_setup,&
get_memory_usage
@ -158,14 +159,14 @@ CONTAINS
REAL(KIND=dp) :: dfp, dhij, dr, drk, drx, eeq_energy, ef_energy, enonbonded, enscale, erep, &
esrb, etaa, etab, f0, f1, f2, fhua, fhub, fhud, foab, fqa, fqb, hij, kf, qlambda, rcova, &
rcovab, rcovb, rrab
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charges, cnumbers, dcharges
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charges, cnumbers, dcharges, qlagrange
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, dqhuckel, huckel, owork
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
REAL(KIND=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_ab, rij, rik
REAL(KIND=dp), DIMENSION(5) :: dpia, dpib, hena, henb, kpolya, kpolyb, &
pia, pib
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:), POINTER :: eeq_q, set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
scon_a, scon_b, wblock, zeta, zetb
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
@ -668,7 +669,8 @@ CONTAINS
! EEQ forces (response and direct)
IF (calculate_forces) THEN
CALL para_env%sum(dcharges)
CALL xtb_eeq_forces(qs_env, charges, dcharges, cnumbers, dcnum, eeq_sparam)
ALLOCATE (qlagrange(natom))
CALL xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
END IF
kf = xtb_control%kf
@ -699,6 +701,17 @@ CONTAINS
CALL para_env%sum(erep)
energy%repulsive = erep
! save EEQ charges
NULLIFY (eeq_q)
CALL get_qs_env(qs_env, eeq=eeq_q)
IF (ASSOCIATED(eeq_q)) THEN
CPASSERT(SIZE(eeq_q) == natom)
ELSE
ALLOCATE (eeq_q(natom))
eeq_q(1:natom) = charges(1:natom)
END IF
CALL set_qs_env(qs_env, eeq=eeq_q)
! deallocate coordination numbers
CALL cnumber_release(cnumbers, dcnum, calculate_forces)
@ -713,7 +726,7 @@ CONTAINS
! deallocate charges
DEALLOCATE (charges)
IF (calculate_forces) THEN
DEALLOCATE (dcharges)
DEALLOCATE (dcharges, qlagrange)
END IF
! AO matrix outputs

9
src/xtb_potentials.F
View file

@ -222,7 +222,7 @@ CONTAINS
dr0, drk, enta, entb, esrbij, etasrb, &
f1, fhua, fhub, gscal, ksrb, rra0, &
rrb0, shift
REAL(KIND=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_rr, rij, rik
REAL(KIND=dp), DIMENSION(3) :: fdik, force_rr, rij, rik
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(atprop_type), POINTER :: atprop
TYPE(neighbor_list_iterator_p_type), &
@ -309,11 +309,10 @@ CONTAINS
rik = dcnum(iatom)%rik(:, i)
drk = SQRT(SUM(rik(:)**2))
IF (drk > 1.e-3_dp) THEN
fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fdika(:)
force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fdika(:)
fdik(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
force(ikind)%repulsive(:, atom_a) = force(ikind)%repulsive(:, atom_a) - fdik(:)
force(kkind)%repulsive(:, atom_c) = force(kkind)%repulsive(:, atom_c) + fdik(:)
IF (use_virial) THEN
fdik = fdika + fdikb
CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
END IF
END IF

404
src/xtb_qresp.F Normal file
View file

@ -0,0 +1,404 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Calculation of charge response in xTB (EEQ only)
!> Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov
!> JCTC 13, 1989-2009, (2017)
!> DOI: 10.1021/acs.jctc.7b00118
!> \author JGH
! **************************************************************************************************
MODULE xtb_qresp
USE ai_contraction, ONLY: block_add,&
contraction
USE ai_overlap, ONLY: overlap_ab
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind_set
USE basis_set_types, ONLY: gto_basis_set_p_type,&
gto_basis_set_type
USE cp_control_types, ONLY: dft_control_type,&
xtb_control_type
USE cp_dbcsr_api, ONLY: dbcsr_add,&
dbcsr_get_block_p,&
dbcsr_p_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE eeq_input, ONLY: eeq_solver_type
USE input_constants, ONLY: vdw_pairpot_dftd4
USE kinds, ONLY: dp
USE kpoint_types, ONLY: get_kpoint_info,&
kpoint_type
USE message_passing, ONLY: mp_para_env_type
USE orbital_pointers, ONLY: ncoset
USE particle_types, ONLY: particle_type
USE qs_dispersion_cnum, ONLY: cnumber_init,&
cnumber_release,&
dcnum_type
USE qs_dispersion_types, ONLY: qs_dispersion_type
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_integral_utils, ONLY: basis_set_list_setup,&
get_memory_usage
USE qs_kind_types, ONLY: get_qs_kind,&
qs_kind_type
USE qs_ks_types, ONLY: get_ks_env,&
qs_ks_env_type
USE qs_neighbor_list_types, ONLY: get_iterator_info,&
neighbor_list_iterate,&
neighbor_list_iterator_create,&
neighbor_list_iterator_p_type,&
neighbor_list_iterator_release,&
neighbor_list_set_p_type
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE xtb_eeq, ONLY: xtb_eeq_calculation,&
xtb_eeq_lagrange
USE xtb_hcore, ONLY: gfn0_huckel,&
gfn0_kpair
USE xtb_types, ONLY: get_xtb_atom_param,&
xtb_atom_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_qresp'
PUBLIC :: build_xtb_qresp
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param qresp ...
! **************************************************************************************************
SUBROUTINE build_xtb_qresp(qs_env, qresp)
TYPE(qs_environment_type), POINTER :: qs_env
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: qresp
INTEGER :: gfn_type
TYPE(dft_control_type), POINTER :: dft_control
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
gfn_type = dft_control%qs_control%xtb_control%gfn_type
qresp = 0.0_dp
SELECT CASE (gfn_type)
CASE (0)
CALL build_gfn0_qresp(qs_env, qresp)
CASE (1)
CPABORT("QRESP: gfn_type = 1 not available")
CASE (2)
CPABORT("QRESP: gfn_type = 2 not available")
CASE DEFAULT
CPABORT("QRESP: Unknown gfn_type")
END SELECT
END SUBROUTINE build_xtb_qresp
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param qresp ...
! **************************************************************************************************
SUBROUTINE build_gfn0_qresp(qs_env, qresp)
TYPE(qs_environment_type), POINTER :: qs_env
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: qresp
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_gfn0_qresp'
INTEGER :: atom_a, atom_b, handle, i, iatom, ic, icol, ikind, img, irow, iset, j, jatom, &
jkind, jset, ldsab, lmaxa, lmaxb, maxder, n1, n2, na, natom, natorb_a, natorb_b, nb, &
ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, nsetb, sgfa, sgfb, za, zb
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
INTEGER, DIMENSION(3) :: cell
INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
npgfb, nsgfa, nsgfb
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: defined, diagblock, do_nonbonded, found
REAL(KIND=dp) :: dr, drx, eeq_energy, ef_energy, etaa, &
etab, f2, fqa, fqb, hij, qlambda, &
rcova, rcovab, rcovb, rrab
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charges, cnumbers, dcharges
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dhuckel, dqhuckel, huckel, owork
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
REAL(KIND=dp), DIMENSION(3) :: rij
REAL(KIND=dp), DIMENSION(5) :: hena, henb, kpolya, kpolyb, pia, pib
REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock, rpgfa, rpgfb, scon_a, scon_b, &
zeta, zetb
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_w
TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
TYPE(dft_control_type), POINTER :: dft_control
TYPE(eeq_solver_type) :: eeq_sparam
TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
TYPE(kpoint_type), POINTER :: kpoints
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: nl_iterator
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
TYPE(xtb_control_type), POINTER :: xtb_control
CALL timeset(routineN, handle)
NULLIFY (logger)
logger => cp_get_default_logger()
NULLIFY (matrix_p, atomic_kind_set, qs_kind_set, sab_orb, ks_env)
CALL get_qs_env(qs_env=qs_env, &
ks_env=ks_env, &
energy=energy, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
para_env=para_env, &
dft_control=dft_control, &
sab_orb=sab_orb)
nkind = SIZE(atomic_kind_set)
xtb_control => dft_control%qs_control%xtb_control
do_nonbonded = xtb_control%do_nonbonded
nimg = dft_control%nimages
nderivatives = 0
maxder = ncoset(nderivatives)
NULLIFY (particle_set)
CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
natom = SIZE(particle_set)
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
atom_of_kind=atom_of_kind, kind_of=kind_of)
NULLIFY (rho, matrix_w)
CALL get_qs_env(qs_env=qs_env, rho=rho)
CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
IF (SIZE(matrix_p, 1) == 2) THEN
DO img = 1, nimg
CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
END DO
END IF
NULLIFY (cell_to_index)
IF (nimg > 1) THEN
CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
END IF
! set up basis set lists
ALLOCATE (basis_set_list(nkind))
CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
! Calculate coordination numbers
! needed for effective atomic energy levels
! code taken from D3 dispersion energy
CALL cnumber_init(qs_env, cnumbers, dcnum, 2, .TRUE.)
ALLOCATE (charges(natom), dcharges(natom))
charges = 0.0_dp
CALL xtb_eeq_calculation(qs_env, charges, cnumbers, eeq_sparam, eeq_energy, ef_energy, qlambda)
dcharges = qlambda/REAL(para_env%num_pe, KIND=dp)
CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
IF (dispersion_env%pp_type == vdw_pairpot_dftd4 .AND. dispersion_env%ext_charges) THEN
IF (ASSOCIATED(dispersion_env%dcharges)) THEN
dcharges(1:natom) = dcharges(1:natom) + dispersion_env%dcharges(1:natom)
ELSE
CPWARN("gfn0-xTB variational dipole DFTD4 contribution missing")
END IF
END IF
! Calculate Huckel parameters
CALL gfn0_huckel(qs_env, cnumbers, charges, huckel, dhuckel, dqhuckel, .TRUE.)
! Calculate KAB parameters and electronegativity correction
CALL gfn0_kpair(qs_env, kijab)
! loop over all atom pairs with a non-zero overlap (sab_orb)
CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
iatom=iatom, jatom=jatom, r=rij, cell=cell)
CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
IF (.NOT. defined .OR. natorb_a < 1) CYCLE
CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
IF (.NOT. defined .OR. natorb_b < 1) CYCLE
dr = SQRT(SUM(rij(:)**2))
! atomic parameters
CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
lmax=lmaxa, nshell=nsa, kpoly=kpolya, hen=hena)
CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
lmax=lmaxb, nshell=nsb, kpoly=kpolyb, hen=henb)
IF (nimg == 1) THEN
ic = 1
ELSE
ic = cell_to_index(cell(1), cell(2), cell(3))
CPASSERT(ic > 0)
END IF
icol = MAX(iatom, jatom)
irow = MIN(iatom, jatom)
NULLIFY (pblock)
CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
row=irow, col=icol, block=pblock, found=found)
CPASSERT(ASSOCIATED(pblock))
! overlap
basis_set_a => basis_set_list(ikind)%gto_basis_set
IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
basis_set_b => basis_set_list(jkind)%gto_basis_set
IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
atom_a = atom_of_kind(iatom)
atom_b = atom_of_kind(jatom)
! basis ikind
first_sgfa => basis_set_a%first_sgf
la_max => basis_set_a%lmax
la_min => basis_set_a%lmin
npgfa => basis_set_a%npgf
nseta = basis_set_a%nset
nsgfa => basis_set_a%nsgf_set
rpgfa => basis_set_a%pgf_radius
set_radius_a => basis_set_a%set_radius
scon_a => basis_set_a%scon
zeta => basis_set_a%zet
! basis jkind
first_sgfb => basis_set_b%first_sgf
lb_max => basis_set_b%lmax
lb_min => basis_set_b%lmin
npgfb => basis_set_b%npgf
nsetb = basis_set_b%nset
nsgfb => basis_set_b%nsgf_set
rpgfb => basis_set_b%pgf_radius
set_radius_b => basis_set_b%set_radius
scon_b => basis_set_b%scon
zetb => basis_set_b%zet
ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
ALLOCATE (sint(natorb_a, natorb_b, maxder))
sint = 0.0_dp
DO iset = 1, nseta
ncoa = npgfa(iset)*ncoset(la_max(iset))
n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
sgfa = first_sgfa(1, iset)
DO jset = 1, nsetb
IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
ncob = npgfb(jset)*ncoset(lb_max(jset))
n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
sgfb = first_sgfb(1, jset)
CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
rij, sab=oint(:, :, 1))
! Contraction
CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
END DO
END DO
! Calculate Pi = Pia * Pib (Eq. 11)
rcovab = rcova + rcovb
rrab = SQRT(dr/rcovab)
pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
! diagonal block
diagblock = .FALSE.
IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
!
f2 = 1.0_dp
IF (iatom /= jatom) f2 = 2.0_dp
! Derivative wrt coordination number
fqa = 0.0_dp
fqb = 0.0_dp
IF (diagblock) THEN
DO i = 1, natorb_a
na = naoa(i)
fqa = fqa + pblock(i, i)*dqhuckel(na, iatom)
END DO
dcharges(iatom) = dcharges(iatom) + fqa
ELSE
DO j = 1, natorb_b
nb = naob(j)
DO i = 1, natorb_a
na = naoa(i)
hij = 0.5_dp*pia(na)*pib(nb)
drx = f2*hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)
IF (iatom <= jatom) THEN
fqa = fqa + drx*pblock(i, j)*dqhuckel(na, iatom)
fqb = fqb + drx*pblock(i, j)*dqhuckel(nb, jatom)
ELSE
fqa = fqa + drx*pblock(j, i)*dqhuckel(na, iatom)
fqb = fqb + drx*pblock(j, i)*dqhuckel(nb, jatom)
END IF
END DO
END DO
dcharges(iatom) = dcharges(iatom) + fqa
dcharges(jatom) = dcharges(jatom) + fqb
END IF
DEALLOCATE (oint, owork, sint)
END DO
CALL neighbor_list_iterator_release(nl_iterator)
! EEQ response
CALL para_env%sum(dcharges)
CALL xtb_eeq_lagrange(qs_env, dcharges, qresp, eeq_sparam)
! deallocate coordination numbers
CALL cnumber_release(cnumbers, dcnum, .TRUE.)
! deallocate Huckel parameters
DEALLOCATE (huckel, dhuckel, dqhuckel)
! deallocate KAB parameters
DEALLOCATE (kijab)
! deallocate charges
DEALLOCATE (charges, dcharges)
DEALLOCATE (basis_set_list)
IF (SIZE(matrix_p, 1) == 2) THEN
DO img = 1, nimg
CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
beta_scalar=-1.0_dp)
END DO
END IF
CALL timestop(handle)
END SUBROUTINE build_gfn0_qresp
END MODULE xtb_qresp

2
tests/QS/regtest-dft-vdw-corr-4/TEST_FILES
View file

@ -9,5 +9,5 @@ pbe_dftd4_stress.inp 31 1.0E-07
ta1.inp 0
ta2.inp 0
ta3.inp 0
ta4.inp 31 1.0E-07 2.83597961983E+04
ta4.inp 31 1.0E-07 2.83598724318E+04
#EOF

2
tests/TEST_DIRS
View file

@ -17,6 +17,7 @@ QS/regtest-debug-5 libint
QS/regtest-gpw-4
ATOM/regtest-1
xTB/regtest-gfn0 libdftd4
xTB/regtest-gfn0-d libdftd4
xTB/regtest-gfn0-m1 libdftd4
xTB/regtest-gfn0-m2 libdftd4
xTB/regtest-gfn0-m3 libdftd4
@ -39,6 +40,7 @@ QS/regtest-ot
QS/regtest-ec-force libint
QS/regtest-embed libint
xTB/regtest-5
xTB/regtest-gfn1-d
QS/regtest-gpw-2-3
QS/regtest-p-efield
SIRIUS/regtest-1 sirius

4
tests/xTB/regtest-2/TEST_FILES
View file

@ -8,8 +8,8 @@ H2O-field-gopt.inp 1 3e-11
H2O-field-gopt-lsd.inp 1 4e-12 -5.76970797174834
H2O-field.inp 0 4e-13
H2O-field-lsd.inp 1 4e-14 -5.76959233315201
HF-field.inp 1 1e-12 -5.62813329327978
HF-field-gopt.inp 1 5e-09 -5.66367831097363
HF-field.inp 1 1e-12 -5.70162277773935
HF-field-gopt.inp 1 5e-09 -5.65593182181073
HF-field-debug.inp 0 1e-12
HF-dfilter-debug.inp 0 1e-12
HF-dfield-gopt.inp 1 1e-09 -5.66065888870620

19
tests/xTB/regtest-gfn0-d/TEST_FILES Normal file
View file

@ -0,0 +1,19 @@
# runs are executed in the same order as in this file
# the second field tells which test should be run in order to compare with the last available output
# e.g. 0 means do not compare anything, running is enough
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
ch2o_t01.inp 0
ch2o_t02.inp 0
ch2o_t03.inp 0
#
ch2o_t04.inp 0
ch2o_t05.inp 0
ch2o_t06.inp 0
#
ch2o_t07.inp 0
##ch2o_t08.inp 0
#
##ch2o_t09.inp 0
##ch2o_t10.inp 0
#EOF

50
tests/xTB/regtest-gfn0-d/ch2o_t01.inp Normal file
View file

@ -0,0 +1,50 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE T
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&EFIELD
&END EFIELD
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
GFN_TYPE 0
VARIATIONAL_DIPOLE
VDW_POTENTIAL DFTD3
&END XTB
&END QS
&SCF
SCF_GUESS NONE
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 20.0 20.0 20.0
PERIODIC NONE
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

49
tests/xTB/regtest-gfn0-d/ch2o_t02.inp Normal file
View file

@ -0,0 +1,49 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE T
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&EFIELD
&END EFIELD
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
GFN_TYPE 0
VARIATIONAL_DIPOLE
&END XTB
&END QS
&SCF
SCF_GUESS NONE
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 20.0 20.0 20.0
PERIODIC NONE
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

49
tests/xTB/regtest-gfn0-d/ch2o_t03.inp Normal file
View file

@ -0,0 +1,49 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES T
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&EFIELD
INTENSITY 0.02
POLARISATION -1 0 0
&END EFIELD
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
GFN_TYPE 0
&END XTB
&END QS
&SCF
SCF_GUESS NONE
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

50
tests/xTB/regtest-gfn0-d/ch2o_t04.inp Normal file
View file

@ -0,0 +1,50 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE T
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&PERIODIC_EFIELD
&END PERIODIC_EFIELD
&PRINT
&MOMENTS ON
PERIODIC .TRUE.
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
CHECK_ATOMIC_CHARGES OFF
GFN_TYPE 0
VARIATIONAL_DIPOLE
VDW_POTENTIAL DFTD3
&END XTB
&END QS
&SCF
SCF_GUESS NONE
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 10.0 10.0 10.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

48
tests/xTB/regtest-gfn0-d/ch2o_t05.inp Normal file
View file

@ -0,0 +1,48 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE T
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&PERIODIC_EFIELD
&END PERIODIC_EFIELD
&PRINT
&MOMENTS ON
PERIODIC .TRUE.
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
GFN_TYPE 0
VARIATIONAL_DIPOLE
&END XTB
&END QS
&SCF
SCF_GUESS NONE
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 10.0 10.0 10.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

49
tests/xTB/regtest-gfn0-d/ch2o_t06.inp Normal file
View file

@ -0,0 +1,49 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES T
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&PERIODIC_EFIELD
INTENSITY 0.05
POLARISATION -1 0 0
&END PERIODIC_EFIELD
&PRINT
&MOMENTS ON
PERIODIC T
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
GFN_TYPE 0
&END XTB
&END QS
&SCF
SCF_GUESS NONE
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

42
tests/xTB/regtest-gfn0-d/ch2o_t07.inp Normal file
View file

@ -0,0 +1,42 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR T
EPS_NO_ERROR_CHECK 1.E-4
MAX_RELATIVE_ERROR 0.10
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
STRESS_TENSOR ANALYTICAL
&DFT
&QS
METHOD xTB
&XTB
GFN_TYPE 0
&END XTB
&END QS
&SCF
SCF_GUESS NONE
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

45
tests/xTB/regtest-gfn0-d/ch2o_t08.inp Normal file
View file

@ -0,0 +1,45 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR T
EPS_NO_ERROR_CHECK 1.E-4
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
STRESS_TENSOR ANALYTICAL
&DFT
&PERIODIC_EFIELD
INTENSITY 0.05
POLARISATION -1 0 0
&END PERIODIC_EFIELD
&QS
METHOD xTB
&XTB
GFN_TYPE 0
&END XTB
&END QS
&SCF
SCF_GUESS NONE
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

46
tests/xTB/regtest-gfn0-d/ch2o_t09.inp Normal file
View file

@ -0,0 +1,46 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES T
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
EPS_NO_ERROR_CHECK 1.E-4
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
STRESS_TENSOR ANALYTICAL
&DFT
&PERIODIC_EFIELD
DISPLACEMENT_FIELD T
INTENSITY 0.05
POLARISATION -1 0 0
&END PERIODIC_EFIELD
&QS
METHOD xTB
&XTB
GFN_TYPE 0
&END XTB
&END QS
&SCF
SCF_GUESS NONE
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

46
tests/xTB/regtest-gfn0-d/ch2o_t10.inp Normal file
View file

@ -0,0 +1,46 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR T
EPS_NO_ERROR_CHECK 1.E-4
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
STRESS_TENSOR ANALYTICAL
&DFT
&PERIODIC_EFIELD
DISPLACEMENT_FIELD T
INTENSITY 0.05
POLARISATION -1 0 0
&END PERIODIC_EFIELD
&QS
METHOD xTB
&XTB
GFN_TYPE 0
&END XTB
&END QS
&SCF
SCF_GUESS NONE
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

19
tests/xTB/regtest-gfn1-d/TEST_FILES Normal file
View file

@ -0,0 +1,19 @@
# runs are executed in the same order as in this file
# the second field tells which test should be run in order to compare with the last available output
# e.g. 0 means do not compare anything, running is enough
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
ch2o_t01.inp 0
ch2o_t02.inp 0
ch2o_t03.inp 0
#
ch2o_t04.inp 0
ch2o_t05.inp 0
ch2o_t06.inp 0
#
ch2o_t07.inp 0
ch2o_t08.inp 0
#
#ch2o_t09.inp 0
#ch2o_t10.inp 0
#EOF

48
tests/xTB/regtest-gfn1-d/ch2o_t01.inp Normal file
View file

@ -0,0 +1,48 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE T
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&EFIELD
&END EFIELD
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
GFN_TYPE 1
&END XTB
&END QS
&SCF
SCF_GUESS MOPAC
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 20.0 20.0 20.0
PERIODIC NONE
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

48
tests/xTB/regtest-gfn1-d/ch2o_t02.inp Normal file
View file

@ -0,0 +1,48 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE T
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&EFIELD
&END EFIELD
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
GFN_TYPE 1
&END XTB
&END QS
&SCF
SCF_GUESS MOPAC
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 20.0 20.0 20.0
PERIODIC NONE
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

50
tests/xTB/regtest-gfn1-d/ch2o_t03.inp Normal file
View file

@ -0,0 +1,50 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES T
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
EPS_NO_ERROR_CHECK 1.E-4
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&EFIELD
INTENSITY 0.02
POLARISATION -1 0 0
&END EFIELD
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
GFN_TYPE 1
&END XTB
&END QS
&SCF
SCF_GUESS MOPAC
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 5.0 5.0 5.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

47
tests/xTB/regtest-gfn1-d/ch2o_t04.inp Normal file
View file

@ -0,0 +1,47 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE T
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&PERIODIC_EFIELD
&END PERIODIC_EFIELD
&PRINT
&MOMENTS ON
PERIODIC .TRUE.
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
GFN_TYPE 1
&END XTB
&END QS
&SCF
SCF_GUESS MOPAC
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 10.0 10.0 10.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

47
tests/xTB/regtest-gfn1-d/ch2o_t05.inp Normal file
View file

@ -0,0 +1,47 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE T
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&PERIODIC_EFIELD
&END PERIODIC_EFIELD
&PRINT
&MOMENTS ON
PERIODIC .TRUE.
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
GFN_TYPE 1
&END XTB
&END QS
&SCF
SCF_GUESS MOPAC
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 10.0 10.0 10.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

49
tests/xTB/regtest-gfn1-d/ch2o_t06.inp Normal file
View file

@ -0,0 +1,49 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES T
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
&DFT
&PERIODIC_EFIELD
INTENSITY 0.05
POLARISATION -1 0 0
&END PERIODIC_EFIELD
&PRINT
&MOMENTS ON
PERIODIC T
REFERENCE COM
&END MOMENTS
&END PRINT
&QS
METHOD xTB
&XTB
GFN_TYPE 1
&END XTB
&END QS
&SCF
SCF_GUESS MOPAC
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

42
tests/xTB/regtest-gfn1-d/ch2o_t07.inp Normal file
View file

@ -0,0 +1,42 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR T
EPS_NO_ERROR_CHECK 1.E-4
MAX_RELATIVE_ERROR 0.10
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
STRESS_TENSOR ANALYTICAL
&DFT
&QS
METHOD xTB
&XTB
GFN_TYPE 1
&END XTB
&END QS
&SCF
SCF_GUESS MOPAC
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

45
tests/xTB/regtest-gfn1-d/ch2o_t08.inp Normal file
View file

@ -0,0 +1,45 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR T
EPS_NO_ERROR_CHECK 1.E-4
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
STRESS_TENSOR ANALYTICAL
&DFT
&PERIODIC_EFIELD
INTENSITY 0.05
POLARISATION -1 0 0
&END PERIODIC_EFIELD
&QS
METHOD xTB
&XTB
GFN_TYPE 1
&END XTB
&END QS
&SCF
SCF_GUESS MOPAC
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

46
tests/xTB/regtest-gfn1-d/ch2o_t09.inp Normal file
View file

@ -0,0 +1,46 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES T
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR F
EPS_NO_ERROR_CHECK 1.E-4
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
STRESS_TENSOR ANALYTICAL
&DFT
&PERIODIC_EFIELD
DISPLACEMENT_FIELD T
INTENSITY 0.05
POLARISATION -1 0 0
&END PERIODIC_EFIELD
&QS
METHOD xTB
&XTB
GFN_TYPE 1
&END XTB
&END QS
&SCF
SCF_GUESS MOPAC
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL

46
tests/xTB/regtest-gfn1-d/ch2o_t10.inp Normal file
View file

@ -0,0 +1,46 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT CH2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
DE 0.0002
DEBUG_DIPOLE F
DEBUG_FORCES F
DEBUG_POLARIZABILITY F
DEBUG_STRESS_TENSOR T
EPS_NO_ERROR_CHECK 1.E-4
STOP_ON_MISMATCH T
&END DEBUG
&FORCE_EVAL
STRESS_TENSOR ANALYTICAL
&DFT
&PERIODIC_EFIELD
DISPLACEMENT_FIELD T
INTENSITY 0.05
POLARISATION -1 0 0
&END PERIODIC_EFIELD
&QS
METHOD xTB
&XTB
GFN_TYPE 1
&END XTB
&END QS
&SCF
SCF_GUESS MOPAC
&END SCF
&END DFT
&SUBSYS
&CELL
ABC 8.0 8.0 8.0
&END CELL
&COORD
O 0.051368 0.000000 0.000000
C 1.278612 0.000000 0.000000
H 1.870460 0.939607 0.000000
H 1.870460 -0.939607 0.000000
&END COORD
&END SUBSYS
&END FORCE_EVAL