From abdbd495f708da2f7f3fa73b7a748807f04bc136 Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Wed, 22 Jan 2025 14:18:16 +0100 Subject: [PATCH] Efield implementation for xTB (#3883) --- src/CMakeLists.txt | 1 + src/cp_control_types.F | 1 + src/cp_control_utils.F | 3 + src/eeq_method.F | 56 +-- src/efield_tb_methods.F | 18 +- src/input_cp2k_tb.F | 6 + src/qs_environment_types.F | 23 +- src/qs_scf_post_gpw.F | 2 +- src/qs_scf_post_se.F | 2 +- src/qs_scf_post_tb.F | 58 +-- src/response_solver.F | 6 + src/xtb_eeq.F | 123 ++++++- src/xtb_matrices.F | 23 +- src/xtb_potentials.F | 9 +- src/xtb_qresp.F | 404 +++++++++++++++++++++ tests/QS/regtest-dft-vdw-corr-4/TEST_FILES | 2 +- tests/TEST_DIRS | 2 + tests/xTB/regtest-2/TEST_FILES | 4 +- tests/xTB/regtest-gfn0-d/TEST_FILES | 19 + tests/xTB/regtest-gfn0-d/ch2o_t01.inp | 50 +++ tests/xTB/regtest-gfn0-d/ch2o_t02.inp | 49 +++ tests/xTB/regtest-gfn0-d/ch2o_t03.inp | 49 +++ tests/xTB/regtest-gfn0-d/ch2o_t04.inp | 50 +++ tests/xTB/regtest-gfn0-d/ch2o_t05.inp | 48 +++ tests/xTB/regtest-gfn0-d/ch2o_t06.inp | 49 +++ tests/xTB/regtest-gfn0-d/ch2o_t07.inp | 42 +++ tests/xTB/regtest-gfn0-d/ch2o_t08.inp | 45 +++ tests/xTB/regtest-gfn0-d/ch2o_t09.inp | 46 +++ tests/xTB/regtest-gfn0-d/ch2o_t10.inp | 46 +++ tests/xTB/regtest-gfn1-d/TEST_FILES | 19 + tests/xTB/regtest-gfn1-d/ch2o_t01.inp | 48 +++ tests/xTB/regtest-gfn1-d/ch2o_t02.inp | 48 +++ tests/xTB/regtest-gfn1-d/ch2o_t03.inp | 50 +++ tests/xTB/regtest-gfn1-d/ch2o_t04.inp | 47 +++ tests/xTB/regtest-gfn1-d/ch2o_t05.inp | 47 +++ tests/xTB/regtest-gfn1-d/ch2o_t06.inp | 49 +++ tests/xTB/regtest-gfn1-d/ch2o_t07.inp | 42 +++ tests/xTB/regtest-gfn1-d/ch2o_t08.inp | 45 +++ tests/xTB/regtest-gfn1-d/ch2o_t09.inp | 46 +++ tests/xTB/regtest-gfn1-d/ch2o_t10.inp | 46 +++ 40 files changed, 1642 insertions(+), 81 deletions(-) create mode 100644 src/xtb_qresp.F create mode 100644 tests/xTB/regtest-gfn0-d/TEST_FILES create mode 100644 tests/xTB/regtest-gfn0-d/ch2o_t01.inp create mode 100644 tests/xTB/regtest-gfn0-d/ch2o_t02.inp create mode 100644 tests/xTB/regtest-gfn0-d/ch2o_t03.inp create mode 100644 tests/xTB/regtest-gfn0-d/ch2o_t04.inp create mode 100644 tests/xTB/regtest-gfn0-d/ch2o_t05.inp create mode 100644 tests/xTB/regtest-gfn0-d/ch2o_t06.inp create mode 100644 tests/xTB/regtest-gfn0-d/ch2o_t07.inp create mode 100644 tests/xTB/regtest-gfn0-d/ch2o_t08.inp create mode 100644 tests/xTB/regtest-gfn0-d/ch2o_t09.inp create mode 100644 tests/xTB/regtest-gfn0-d/ch2o_t10.inp create mode 100644 tests/xTB/regtest-gfn1-d/TEST_FILES create mode 100644 tests/xTB/regtest-gfn1-d/ch2o_t01.inp create mode 100644 tests/xTB/regtest-gfn1-d/ch2o_t02.inp create mode 100644 tests/xTB/regtest-gfn1-d/ch2o_t03.inp create mode 100644 tests/xTB/regtest-gfn1-d/ch2o_t04.inp create mode 100644 tests/xTB/regtest-gfn1-d/ch2o_t05.inp create mode 100644 tests/xTB/regtest-gfn1-d/ch2o_t06.inp create mode 100644 tests/xTB/regtest-gfn1-d/ch2o_t07.inp create mode 100644 tests/xTB/regtest-gfn1-d/ch2o_t08.inp create mode 100644 tests/xTB/regtest-gfn1-d/ch2o_t09.inp create mode 100644 tests/xTB/regtest-gfn1-d/ch2o_t10.inp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a448c47176..1c97dbc507 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 diff --git a/src/cp_control_types.F b/src/cp_control_types.F index d1e7e4c0d7..b12bbaf160 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -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 diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index 7cfa1e2ed9..4399eae7a9 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -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 diff --git a/src/eeq_method.F b/src/eeq_method.F index 152c49111c..b03cf7401d 100644 --- a/src/eeq_method.F +++ b/src/eeq_method.F @@ -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 diff --git a/src/efield_tb_methods.F b/src/efield_tb_methods.F index 54cfa72c39..98c336153f 100644 --- a/src/efield_tb_methods.F +++ b/src/efield_tb_methods.F @@ -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 diff --git a/src/input_cp2k_tb.F b/src/input_cp2k_tb.F index 1c17de8ece..bdb16cf7f2 100644 --- a/src/input_cp2k_tb.F +++ b/src/input_cp2k_tb.F @@ -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) diff --git a/src/qs_environment_types.F b/src/qs_environment_types.F index 5448fc3e63..b3de47c912 100644 --- a/src/qs_environment_types.F +++ b/src/qs_environment_types.F @@ -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) diff --git a/src/qs_scf_post_gpw.F b/src/qs_scf_post_gpw.F index 2a7e0c1446..58a007b9d4 100644 --- a/src/qs_scf_post_gpw.F +++ b/src/qs_scf_post_gpw.F @@ -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 diff --git a/src/qs_scf_post_se.F b/src/qs_scf_post_se.F index 126e30706c..f0e26bf5c6 100644 --- a/src/qs_scf_post_se.F +++ b/src/qs_scf_post_se.F @@ -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 diff --git a/src/qs_scf_post_tb.F b/src/qs_scf_post_tb.F index c898e0c71e..bb4ac4ef40 100644 --- a/src/qs_scf_post_tb.F +++ b/src/qs_scf_post_tb.F @@ -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))') & diff --git a/src/response_solver.F b/src/response_solver.F index 30f07e48c1..c60a714762 100644 --- a/src/response_solver.F +++ b/src/response_solver.F @@ -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) diff --git a/src/xtb_eeq.F b/src/xtb_eeq.F index b7fc6a9861..4432cca847 100644 --- a/src/xtb_eeq.F +++ b/src/xtb_eeq.F @@ -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 diff --git a/src/xtb_matrices.F b/src/xtb_matrices.F index 513733507d..454ea2d0d4 100644 --- a/src/xtb_matrices.F +++ b/src/xtb_matrices.F @@ -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 diff --git a/src/xtb_potentials.F b/src/xtb_potentials.F index ff9a0770b2..fb07187854 100644 --- a/src/xtb_potentials.F +++ b/src/xtb_potentials.F @@ -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 diff --git a/src/xtb_qresp.F b/src/xtb_qresp.F new file mode 100644 index 0000000000..b1aab49250 --- /dev/null +++ b/src/xtb_qresp.F @@ -0,0 +1,404 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2025 CP2K developers group ! +! ! +! 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 + diff --git a/tests/QS/regtest-dft-vdw-corr-4/TEST_FILES b/tests/QS/regtest-dft-vdw-corr-4/TEST_FILES index d239df0fc7..d33122c15f 100644 --- a/tests/QS/regtest-dft-vdw-corr-4/TEST_FILES +++ b/tests/QS/regtest-dft-vdw-corr-4/TEST_FILES @@ -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 diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index a15e5e7eb9..f6fc4b9e35 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -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 diff --git a/tests/xTB/regtest-2/TEST_FILES b/tests/xTB/regtest-2/TEST_FILES index d657df7964..cfd612c561 100644 --- a/tests/xTB/regtest-2/TEST_FILES +++ b/tests/xTB/regtest-2/TEST_FILES @@ -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 diff --git a/tests/xTB/regtest-gfn0-d/TEST_FILES b/tests/xTB/regtest-gfn0-d/TEST_FILES new file mode 100644 index 0000000000..0edfb5b30f --- /dev/null +++ b/tests/xTB/regtest-gfn0-d/TEST_FILES @@ -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 diff --git a/tests/xTB/regtest-gfn0-d/ch2o_t01.inp b/tests/xTB/regtest-gfn0-d/ch2o_t01.inp new file mode 100644 index 0000000000..f5f27f2948 --- /dev/null +++ b/tests/xTB/regtest-gfn0-d/ch2o_t01.inp @@ -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 diff --git a/tests/xTB/regtest-gfn0-d/ch2o_t02.inp b/tests/xTB/regtest-gfn0-d/ch2o_t02.inp new file mode 100644 index 0000000000..027a038ace --- /dev/null +++ b/tests/xTB/regtest-gfn0-d/ch2o_t02.inp @@ -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 diff --git a/tests/xTB/regtest-gfn0-d/ch2o_t03.inp b/tests/xTB/regtest-gfn0-d/ch2o_t03.inp new file mode 100644 index 0000000000..9e6a1284e0 --- /dev/null +++ b/tests/xTB/regtest-gfn0-d/ch2o_t03.inp @@ -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 diff --git a/tests/xTB/regtest-gfn0-d/ch2o_t04.inp b/tests/xTB/regtest-gfn0-d/ch2o_t04.inp new file mode 100644 index 0000000000..79957208eb --- /dev/null +++ b/tests/xTB/regtest-gfn0-d/ch2o_t04.inp @@ -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 diff --git a/tests/xTB/regtest-gfn0-d/ch2o_t05.inp b/tests/xTB/regtest-gfn0-d/ch2o_t05.inp new file mode 100644 index 0000000000..29ba4a3f04 --- /dev/null +++ b/tests/xTB/regtest-gfn0-d/ch2o_t05.inp @@ -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 diff --git a/tests/xTB/regtest-gfn0-d/ch2o_t06.inp b/tests/xTB/regtest-gfn0-d/ch2o_t06.inp new file mode 100644 index 0000000000..e1f27930f5 --- /dev/null +++ b/tests/xTB/regtest-gfn0-d/ch2o_t06.inp @@ -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 diff --git a/tests/xTB/regtest-gfn0-d/ch2o_t07.inp b/tests/xTB/regtest-gfn0-d/ch2o_t07.inp new file mode 100644 index 0000000000..8e37095e3b --- /dev/null +++ b/tests/xTB/regtest-gfn0-d/ch2o_t07.inp @@ -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 diff --git a/tests/xTB/regtest-gfn0-d/ch2o_t08.inp b/tests/xTB/regtest-gfn0-d/ch2o_t08.inp new file mode 100644 index 0000000000..d9117a536f --- /dev/null +++ b/tests/xTB/regtest-gfn0-d/ch2o_t08.inp @@ -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 diff --git a/tests/xTB/regtest-gfn0-d/ch2o_t09.inp b/tests/xTB/regtest-gfn0-d/ch2o_t09.inp new file mode 100644 index 0000000000..e0eccb0f32 --- /dev/null +++ b/tests/xTB/regtest-gfn0-d/ch2o_t09.inp @@ -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 diff --git a/tests/xTB/regtest-gfn0-d/ch2o_t10.inp b/tests/xTB/regtest-gfn0-d/ch2o_t10.inp new file mode 100644 index 0000000000..5b5f0bd445 --- /dev/null +++ b/tests/xTB/regtest-gfn0-d/ch2o_t10.inp @@ -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 diff --git a/tests/xTB/regtest-gfn1-d/TEST_FILES b/tests/xTB/regtest-gfn1-d/TEST_FILES new file mode 100644 index 0000000000..d004b112da --- /dev/null +++ b/tests/xTB/regtest-gfn1-d/TEST_FILES @@ -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 diff --git a/tests/xTB/regtest-gfn1-d/ch2o_t01.inp b/tests/xTB/regtest-gfn1-d/ch2o_t01.inp new file mode 100644 index 0000000000..d067e8f636 --- /dev/null +++ b/tests/xTB/regtest-gfn1-d/ch2o_t01.inp @@ -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 diff --git a/tests/xTB/regtest-gfn1-d/ch2o_t02.inp b/tests/xTB/regtest-gfn1-d/ch2o_t02.inp new file mode 100644 index 0000000000..d067e8f636 --- /dev/null +++ b/tests/xTB/regtest-gfn1-d/ch2o_t02.inp @@ -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 diff --git a/tests/xTB/regtest-gfn1-d/ch2o_t03.inp b/tests/xTB/regtest-gfn1-d/ch2o_t03.inp new file mode 100644 index 0000000000..b646fed40a --- /dev/null +++ b/tests/xTB/regtest-gfn1-d/ch2o_t03.inp @@ -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 diff --git a/tests/xTB/regtest-gfn1-d/ch2o_t04.inp b/tests/xTB/regtest-gfn1-d/ch2o_t04.inp new file mode 100644 index 0000000000..8b8dd88997 --- /dev/null +++ b/tests/xTB/regtest-gfn1-d/ch2o_t04.inp @@ -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 diff --git a/tests/xTB/regtest-gfn1-d/ch2o_t05.inp b/tests/xTB/regtest-gfn1-d/ch2o_t05.inp new file mode 100644 index 0000000000..8b8dd88997 --- /dev/null +++ b/tests/xTB/regtest-gfn1-d/ch2o_t05.inp @@ -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 diff --git a/tests/xTB/regtest-gfn1-d/ch2o_t06.inp b/tests/xTB/regtest-gfn1-d/ch2o_t06.inp new file mode 100644 index 0000000000..d06428df8c --- /dev/null +++ b/tests/xTB/regtest-gfn1-d/ch2o_t06.inp @@ -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 diff --git a/tests/xTB/regtest-gfn1-d/ch2o_t07.inp b/tests/xTB/regtest-gfn1-d/ch2o_t07.inp new file mode 100644 index 0000000000..09051ef2eb --- /dev/null +++ b/tests/xTB/regtest-gfn1-d/ch2o_t07.inp @@ -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 diff --git a/tests/xTB/regtest-gfn1-d/ch2o_t08.inp b/tests/xTB/regtest-gfn1-d/ch2o_t08.inp new file mode 100644 index 0000000000..9dc0510de0 --- /dev/null +++ b/tests/xTB/regtest-gfn1-d/ch2o_t08.inp @@ -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 diff --git a/tests/xTB/regtest-gfn1-d/ch2o_t09.inp b/tests/xTB/regtest-gfn1-d/ch2o_t09.inp new file mode 100644 index 0000000000..c5b29eef10 --- /dev/null +++ b/tests/xTB/regtest-gfn1-d/ch2o_t09.inp @@ -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 diff --git a/tests/xTB/regtest-gfn1-d/ch2o_t10.inp b/tests/xTB/regtest-gfn1-d/ch2o_t10.inp new file mode 100644 index 0000000000..49574cd6be --- /dev/null +++ b/tests/xTB/regtest-gfn1-d/ch2o_t10.inp @@ -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