From 9838f79063eb5dc115d6559fee9eec3e928e6030 Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Mon, 13 Jul 2026 10:17:38 +0200 Subject: [PATCH] xTB spin polarisation Hamiltonian native implementation --- data/xTB_spin_pol_params | 80 +++++++ src/CMakeLists.txt | 1 + src/cp_control_types.F | 16 +- src/cp_control_utils.F | 23 ++ src/input_cp2k_tb.F | 21 ++ src/qs_diis.F | 3 + src/qs_energy_types.F | 3 +- src/qs_environment.F | 10 +- src/qs_initial_guess.F | 7 +- src/qs_ks_utils.F | 2 +- src/xtb_coulomb.F | 7 + src/xtb_ks_matrix.F | 9 +- src/xtb_parameters.F | 107 +++++++++ src/xtb_spinpol.F | 496 +++++++++++++++++++++++++++++++++++++++ src/xtb_types.F | 31 ++- 15 files changed, 803 insertions(+), 13 deletions(-) create mode 100644 data/xTB_spin_pol_params create mode 100644 src/xtb_spinpol.F diff --git a/data/xTB_spin_pol_params b/data/xTB_spin_pol_params new file mode 100644 index 0000000000..50f52e32f1 --- /dev/null +++ b/data/xTB_spin_pol_params @@ -0,0 +1,80 @@ +# Spin paramters for gfn1-xTB (units Eh) +# +#High-throughput screening of spin states for transition metal +#complexes with spin-polarized extended tight-binding methods +#Hagen Neugebauer, Benedikt Baedorf, Sebastian Ehlert, Andreas Hansen, Stefan Grimme +#J Comput Chem. 44:2120-2129 (2023) +# +# element Wss Wsp Wpp Wsd Wpd Wdd +1 0.071550 0.000000 0.000000 0.000000 0.000000 0.000000 +2 0.614675 0.033587 0.125800 0.000000 0.000000 0.000000 +3 0.017775 0.013937 0.018050 0.000000 0.000000 0.000000 +4 0.022850 0.018612 0.017575 0.000000 0.000000 0.000000 +5 0.027325 0.022037 0.019600 0.000000 0.000000 0.000000 +6 0.030200 0.025025 0.022725 0.000000 0.000000 0.000000 +7 0.033000 0.027475 0.025475 0.000000 0.000000 0.000000 +8 0.035100 0.029500 0.027825 0.000000 0.000000 0.000000 +9 0.036900 0.031200 0.029900 0.000000 0.000000 0.000000 +10 0.055008 0.012830 0.022600 0.011925 0.016737 0.080725 +11 0.015100 0.013337 0.023025 0.000000 0.000000 0.000000 +12 0.016500 0.013175 0.017400 0.000000 0.000000 0.000000 +13 0.018250 0.013837 0.014000 0.008175 0.011637 0.012875 +14 0.019525 0.015000 0.014350 0.008437 0.011637 0.014075 +15 0.020550 0.016112 0.014900 0.009300 0.011975 0.014825 +16 0.021325 0.017012 0.015500 0.009987 0.012137 0.014950 +17 0.021825 0.017712 0.016075 0.010987 0.012612 0.015075 +18 0.342475 0.077825 0.120675 0.015500 0.023075 0.051925 +19 0.010650 0.010900 0.016375 0.000000 0.000000 0.000000 +20 0.011800 0.010387 0.013350 0.005562 0.003512 0.010200 +21 0.012725 0.010912 0.013850 0.004787 0.002412 0.012525 +22 0.013525 0.011225 0.014675 0.004350 0.001975 0.013900 +23 0.014075 0.011512 0.015275 0.004037 0.001725 0.014900 +24 0.015175 0.012450 0.021225 0.004150 0.001662 0.013875 +25 0.015000 0.011787 0.016725 0.003550 0.001325 0.016525 +26 0.015400 0.011925 0.017850 0.003300 0.001162 0.017125 +27 0.015825 0.012037 0.018700 0.003137 0.001050 0.017750 +28 0.016150 0.012175 0.019700 0.002987 0.000950 0.018300 +29 0.017150 0.013175 0.030375 0.002775 0.000650 0.017475 +30 0.016850 0.012312 0.021450 0.000000 0.000000 0.000000 +31 0.017225 0.012787 0.013400 0.008525 0.013000 0.015775 +32 0.017550 0.013375 0.013575 0.008112 0.012825 0.017525 +33 0.017750 0.013762 0.013600 0.007987 0.012387 0.017550 +34 0.017975 0.014087 0.013625 0.008162 0.011962 0.017200 +35 0.018100 0.014375 0.013725 0.008275 0.011762 0.016675 +36 0.299025 0.066587 0.101875 0.012575 0.021287 0.048300 +37 0.009550 0.009600 0.016725 0.000000 0.000000 0.000000 +38 0.010650 0.009237 0.012525 0.000000 0.000000 0.000000 +39 0.011425 0.009487 0.012300 0.006725 0.003987 0.009725 +40 0.011950 0.009612 0.013525 0.006150 0.003075 0.010725 +41 0.012575 0.010262 0.019075 0.006062 0.002887 0.010475 +42 0.012925 0.010500 0.022225 0.005562 0.002362 0.010925 +43 0.013150 0.010662 0.024725 0.005112 0.002025 0.011300 +44 0.013375 0.010750 0.027500 0.004750 0.001662 0.011625 +45 0.013525 0.010912 0.032025 0.004400 0.001425 0.011875 +46 0.018975 0.023937 0.180200 0.002087 0.001487 0.011325 +47 0.013925 0.011100 0.039800 0.003887 0.001012 0.012400 +48 0.013850 0.010500 0.019650 0.000000 0.000000 0.000000 +49 0.014125 0.010550 0.011575 0.005062 0.009375 0.010100 +50 0.014300 0.010912 0.011675 0.004600 0.009125 0.011875 +51 0.014525 0.011125 0.011650 0.004375 0.008725 0.012525 +52 0.014525 0.011237 0.011550 0.004137 0.008162 0.012250 +53 0.014575 0.011337 0.011450 0.004450 0.008312 0.012825 +54 0.255850 0.055587 0.085625 0.004662 0.013337 0.037350 +55 0.008200 0.008575 0.015300 0.000000 0.000000 0.000000 +56 0.009275 0.008200 0.011250 0.000000 0.000000 0.000000 +57 0.009925 0.008412 0.011400 0.005925 0.003312 0.009025 +72 0.012175 0.009625 0.012600 0.007637 0.004187 0.010425 +73 0.012325 0.009575 0.013375 0.007137 0.003475 0.010925 +74 0.012500 0.009562 0.014450 0.006725 0.002950 0.011225 +75 0.012600 0.009662 0.014800 0.006275 0.002600 0.011450 +76 0.012600 0.009200 0.020600 0.005950 0.002075 0.011550 +77 0.012725 0.009275 0.021000 0.005700 0.001912 0.011650 +78 0.013075 0.010212 0.033575 0.005550 0.001812 0.011150 +79 0.013150 0.009962 0.053000 0.005287 0.001462 0.011175 +80 0.013025 0.009187 0.029250 0.000000 0.000000 0.000000 +81 0.013275 0.009112 0.010725 0.000000 0.000000 0.000000 +82 0.013475 0.009350 0.010975 0.000000 0.000000 0.000000 +83 0.013625 0.009537 0.010975 0.000000 0.000000 0.000000 +84 0.013725 0.009625 0.010850 0.000000 0.000000 0.000000 +85 0.013775 0.009737 0.010725 0.002612 0.007362 0.011925 +86 0.254400 0.050400 0.080625 0.001087 0.011050 0.035175 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 89417a3cf5..f805f677f7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -884,6 +884,7 @@ list( xray_diffraction.F xtb_qresp.F xtb_coulomb.F + xtb_spinpol.F xtb_eeq.F xtb_ehess.F xtb_ehess_force.F diff --git a/src/cp_control_types.F b/src/cp_control_types.F index 4470a26447..ff28629fb2 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -278,6 +278,7 @@ MODULE cp_control_types INTEGER :: vdw_type = -1 CHARACTER(LEN=default_path_length) :: parameter_file_path = "" CHARACTER(LEN=default_path_length) :: parameter_file_name = "" + CHARACTER(LEN=default_path_length) :: spinpol_param_file_name = "" ! CHARACTER(LEN=default_path_length) :: dispersion_parameter_file = "" REAL(KIND=dp) :: epscn = 0.0_dp @@ -296,6 +297,7 @@ MODULE cp_control_types ! LOGICAL :: xb_interaction = .FALSE. LOGICAL :: do_nonbonded = .FALSE. + LOGICAL :: do_spinpol = .FALSE. LOGICAL :: coulomb_interaction = .FALSE. LOGICAL :: coulomb_lr = .FALSE. LOGICAL :: tb3_interaction = .FALSE. @@ -310,7 +312,11 @@ MODULE cp_control_types DIMENSION(:, :), POINTER :: kab_param => NULL() INTEGER, DIMENSION(:, :), POINTER :: kab_types => NULL() INTEGER :: kab_nval = 0 - REAL, DIMENSION(:), POINTER :: kab_vals => NULL() + REAL(KIND=dp), DIMENSION(:), POINTER :: kab_vals => NULL() + ! + INTEGER, DIMENSION(:), POINTER :: spinpol_type => NULL() + REAL(KIND=dp), DIMENSION(:, :), & + POINTER :: spinpol_vals => NULL() ! TYPE(pair_potential_p_type), POINTER :: nonbonded => NULL() REAL(KIND=dp) :: eps_pair = 0.0_dp @@ -1303,6 +1309,8 @@ CONTAINS NULLIFY (xtb_control%kab_types) NULLIFY (xtb_control%nonbonded) NULLIFY (xtb_control%rcpair) + NULLIFY (xtb_control%spinpol_type) + NULLIFY (xtb_control%spinpol_vals) END SUBROUTINE xtb_control_create @@ -1329,6 +1337,12 @@ CONTAINS IF (ASSOCIATED(xtb_control%nonbonded)) THEN CALL pair_potential_p_release(xtb_control%nonbonded) END IF + IF (ASSOCIATED(xtb_control%spinpol_type)) THEN + DEALLOCATE (xtb_control%spinpol_type) + END IF + IF (ASSOCIATED(xtb_control%spinpol_vals)) THEN + DEALLOCATE (xtb_control%spinpol_vals) + END IF DEALLOCATE (xtb_control) END IF END SUBROUTINE xtb_control_release diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index abe607899e..c017d89aa7 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -1579,6 +1579,9 @@ CONTAINS ELSE qs_control%xtb_control%do_ewald = (qs_control%periodicity /= 0) END IF + ! Spin Polarisation + CALL section_vals_val_get(xtb_section, "SPIN_POLARISATION", & + l_val=qs_control%xtb_control%do_spinpol) ! vdW CALL section_vals_val_get(xtb_section, "VDW_POTENTIAL", explicit=explicit) IF (explicit) THEN @@ -1630,6 +1633,9 @@ CONTAINS CPABORT("GFN type") END SELECT END IF + ! + CALL section_vals_val_get(xtb_parameter, "SPINPOL_PARAM_FILE_NAME", & + c_val=qs_control%xtb_control%spinpol_param_file_name) ! D3 Dispersion CALL section_vals_val_get(xtb_parameter, "DISPERSION_RADIUS", & r_val=qs_control%xtb_control%rcdisp) @@ -1853,6 +1859,23 @@ CONTAINS END DO END IF + ! Spin Polarisation + CALL section_vals_val_get(xtb_parameter, "SPIN_POL_PARAM", n_rep_val=n_rep) + IF (n_rep > 0) THEN + ALLOCATE (qs_control%xtb_control%spinpol_type(n_rep)) + ALLOCATE (qs_control%xtb_control%spinpol_vals(6, n_rep)) + DO j = 1, n_rep + CALL section_vals_val_get(xtb_parameter, "SPIN_POL_PARAM", i_rep_val=j, c_vals=clist) + READ (clist(1), '(I3)') qs_control%xtb_control%spinpol_type(j) + READ (clist(2), '(F20.8)') qs_control%xtb_control%spinpol_vals(1, j) + READ (clist(3), '(F20.8)') qs_control%xtb_control%spinpol_vals(2, j) + READ (clist(4), '(F20.8)') qs_control%xtb_control%spinpol_vals(3, j) + READ (clist(5), '(F20.8)') qs_control%xtb_control%spinpol_vals(4, j) + READ (clist(6), '(F20.8)') qs_control%xtb_control%spinpol_vals(5, j) + READ (clist(7), '(F20.8)') qs_control%xtb_control%spinpol_vals(6, j) + END DO + END IF + IF (qs_control%xtb_control%gfn_type == 0) THEN CALL section_vals_val_get(xtb_parameter, "SRB_PARAMETER", r_vals=scal) qs_control%xtb_control%ksrb = scal(1) diff --git a/src/input_cp2k_tb.F b/src/input_cp2k_tb.F index ba3342df1b..4af089e227 100644 --- a/src/input_cp2k_tb.F +++ b/src/input_cp2k_tb.F @@ -243,6 +243,12 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="SPIN_POLARISATION", & + description="Use the spin polarisation Hamiltonian for gfn1/2", & + usage="SPIN_POLARISATION 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="COULOMB_INTERACTION", & description="Use Coulomb interaction terms (electrostatics + TB3); for debug only", & usage="COULOMB_INTERACTION T", default_l_val=.TRUE., lone_keyword_l_val=.TRUE.) @@ -438,6 +444,14 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="SPINPOL_PARAM_FILE_NAME", & + description="Specify file that contains parameters for "// & + "xTB spin polarisation Hamiltonian", & + usage="SPINPOL_PARAM_FILE_NAME filename", & + n_var=1, type_of_var=char_t, default_c_val="xTB_spin_pol_params") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="DISPERSION_PARAMETER_FILE", & description="Specify file that contains the atomic dispersion "// & "parameters for the D3 method", & @@ -525,6 +539,13 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="SPIN_POL_PARAM", & + description="Specifies the spin polarisation parameters for kind A.", & + usage="SPIN_POL_PARAM atomtype Wss Wsp Wpp Wsd Wpd Wdd", repeats=.TRUE., & + n_var=-1, type_of_var=char_t) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="XB_RADIUS", & description="Specifies the radius [Bohr] of the XB pair interaction in xTB.", & usage="XB_RADIUS 20.0 ", repeats=.FALSE., & diff --git a/src/qs_diis.F b/src/qs_diis.F index 5b89b5f419..cbc8b1596b 100644 --- a/src/qs_diis.F +++ b/src/qs_diis.F @@ -571,6 +571,9 @@ CONTAINS TYPE(mp_para_env_type), POINTER :: para_env CALL timeset(routineN, handle) + IF (ls_scf_env%do_pao) THEN + CPABORT("LS_SCF%LS_DIIS not compatible with PAO") + END IF nspin = ls_scf_env%nspins diis_step = .FALSE. my_nmixing = 2 diff --git a/src/qs_energy_types.F b/src/qs_energy_types.F index bb58784a8d..f527355b98 100644 --- a/src/qs_energy_types.F +++ b/src/qs_energy_types.F @@ -80,7 +80,8 @@ MODULE qs_energy_types surf_dipole = 0.0_dp, & embed_corr = 0.0_dp, & ! correction for embedding potential xtb_xb_inter = 0.0_dp, & ! correction for halogen bonding within GFN1-xTB - xtb_nonbonded = 0.0_dp ! correction for nonbonded interactions within GFN1-xTB + xtb_nonbonded = 0.0_dp, & ! correction for nonbonded interactions within GFN1-xTB + xtb_spinpol = 0.0_dp ! spin-polarised Hamiltonian energy within GFN1/2-xTB REAL(KIND=dp), DIMENSION(:), POINTER :: ddapc_restraint => NULL() END TYPE qs_energy_type diff --git a/src/qs_environment.F b/src/qs_environment.F index 08194e09b6..1743d76cea 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -218,7 +218,9 @@ MODULE qs_environment USE transport, ONLY: transport_env_create USE xtb_parameters, ONLY: init_xtb_basis,& xtb_parameters_init,& - xtb_parameters_set + xtb_parameters_set,& + xtb_spinpol_ext,& + xtb_spinpol_init USE xtb_potentials, ONLY: xtb_pp_radius USE xtb_types, ONLY: allocate_xtb_atom_param,& set_xtb_atom_param @@ -1238,6 +1240,12 @@ CONTAINS CALL xtb_parameters_init(qs_kind%xtb_parameter, gfn_type, element_symbol, & xtb_control%parameter_file_path, xtb_control%parameter_file_name, & para_env) + IF (xtb_control%do_spinpol) THEN + CALL xtb_spinpol_init(qs_kind%xtb_parameter, gfn_type, element_symbol, & + xtb_control%parameter_file_path, xtb_control%spinpol_param_file_name, & + para_env) + CALL xtb_spinpol_ext(qs_kind%xtb_parameter, gfn_type, xtb_control) + END IF ! set dependent parameters CALL xtb_parameters_set(qs_kind%xtb_parameter) ! Generate basis set diff --git a/src/qs_initial_guess.F b/src/qs_initial_guess.F index 0df0b262a5..d5767ed9a6 100644 --- a/src/qs_initial_guess.F +++ b/src/qs_initial_guess.F @@ -1259,7 +1259,6 @@ CONTAINS pdiag(:) = 0.0_dp ALLOCATE (sdiag(nao)) - sdiag(:) = 0.0_dp IF (has_unit_metric) THEN sdiag(:) = 1.0_dp @@ -1341,13 +1340,14 @@ CONTAINS isgfa = first_sgf(atom_a) IF (z == 1 .AND. nsgf == 2) THEN ! Hydrogen 2s basis - pdiag(isgfa) = 1.0_dp + pdiag(isgfa) = 1.0_dp/REAL(nspin, dp) pdiag(isgfa + 1) = 0.0_dp ELSE DO isgf = 1, nsgf na = naox(isgf) la = laox(isgf) occ = REAL(occupation(la + 1), dp)/REAL(2*la + 1, dp) + occ = occ/REAL(nspin, dp) pdiag(isgfa + isgf - 1) = occ END DO END IF @@ -1456,8 +1456,11 @@ CONTAINS END DO DO ispin = 1, nspin IF (nelectron_spin(ispin) /= 0) THEN + rscale = SUM(pdiag)/REAL(nelectron_spin(ispin), dp) matrix_p => pmat(ispin)%matrix + pdiag = rscale*pdiag CALL dbcsr_set_diag(matrix_p, pdiag) + pdiag = pdiag/rscale END IF END DO ELSE diff --git a/src/qs_ks_utils.F b/src/qs_ks_utils.F index aa83ab16a1..00c65cf2cd 100644 --- a/src/qs_ks_utils.F +++ b/src/qs_ks_utils.F @@ -1088,7 +1088,7 @@ CONTAINS "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit END IF ELSE -!ZMP to print some variables at each step + !ZMP to print some variables at each step IF (dft_control%apply_external_density) THEN WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & "DOING ZMP CALCULATION FROM EXTERNAL DENSITY " diff --git a/src/xtb_coulomb.F b/src/xtb_coulomb.F index a95fdcede7..16de64133a 100644 --- a/src/xtb_coulomb.F +++ b/src/xtb_coulomb.F @@ -73,6 +73,7 @@ MODULE xtb_coulomb sap_int_type USE virial_methods, ONLY: virial_pair_force USE virial_types, ONLY: virial_type + USE xtb_spinpol, ONLY: build_xtb_spinpol USE xtb_types, ONLY: get_xtb_atom_param,& xtb_atom_type #include "./base/base_uses.f90" @@ -644,6 +645,12 @@ CONTAINS DEALLOCATE (zeffk, xgamma) END IF + IF (xtb_control%do_spinpol) THEN + CALL qs_rho_get(rho, rho_ao_kp=matrix_p) + CALL build_xtb_spinpol(qs_env, ks_matrix, matrix_p, energy, & + sap_int, calculate_forces, just_energy) + END IF + ! QMMM IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN CALL build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, & diff --git a/src/xtb_ks_matrix.F b/src/xtb_ks_matrix.F index b9725df4fe..133f8c786c 100644 --- a/src/xtb_ks_matrix.F +++ b/src/xtb_ks_matrix.F @@ -431,8 +431,9 @@ CONTAINS energy%qmmm_el = energy%qmmm_el + pc_ener END IF - energy%total = energy%core + energy%hartree + energy%efield + energy%qmmm_el + & - energy%repulsive + energy%dispersion + energy%dftb3 + energy%kTS + energy%total = energy%core + energy%repulsive + & + energy%hartree + energy%xtb_spinpol + energy%efield + & + energy%qmmm_el + energy%dispersion + energy%dftb3 + energy%kTS iounit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DETAILED_ENERGY", & extension=".scfLog") @@ -442,6 +443,10 @@ CONTAINS "Zeroth order Hamiltonian energy: ", energy%core, & "Charge fluctuation energy: ", energy%hartree, & "London dispersion energy: ", energy%dispersion + IF (dft_control%qs_control%xtb_control%do_spinpol) THEN + WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") & + "Spin polarisation correction: ", energy%xtb_spinpol + END IF IF (dft_control%qs_control%xtb_control%xb_interaction) THEN WRITE (UNIT=iounit, FMT="(T9,A,T60,F20.10)") & "Correction for halogen bonding: ", energy%xtb_xb_inter diff --git a/src/xtb_parameters.F b/src/xtb_parameters.F index 67bb511bd2..5d0f343518 100644 --- a/src/xtb_parameters.F +++ b/src/xtb_parameters.F @@ -165,6 +165,7 @@ MODULE xtb_parameters ! *** Public data types *** PUBLIC :: xtb_parameters_init, xtb_parameters_set, init_xtb_basis, xtb_set_kab + PUBLIC :: xtb_spinpol_init, xtb_spinpol_ext PUBLIC :: metal, early3d, pp_gfn0 CONTAINS @@ -468,6 +469,112 @@ CONTAINS END SUBROUTINE xtb1_parameters_init +! ************************************************************************************************** +!> \brief ... +!> \param param ... +!> \param gfn_type ... +!> \param element_symbol ... +!> \param parameter_file_path ... +!> \param spinpol_param_file_name ... +!> \param para_env ... +! ************************************************************************************************** + SUBROUTINE xtb_spinpol_init(param, gfn_type, element_symbol, parameter_file_path, spinpol_param_file_name, & + para_env) + + TYPE(xtb_atom_type), POINTER :: param + INTEGER, INTENT(IN) :: gfn_type + CHARACTER(LEN=2), INTENT(IN) :: element_symbol + CHARACTER(LEN=*), INTENT(IN) :: parameter_file_path, & + spinpol_param_file_name + TYPE(mp_para_env_type), POINTER :: para_env + + CHARACTER(len=default_string_length) :: filename + INTEGER :: zin, znum + LOGICAL :: at_end + TYPE(cp_parser_type) :: parser + + SELECT CASE (gfn_type) + CASE (0) + CPABORT("gfn_type = 0: No spin polarisation possible!") + CASE (1) + ! OK + CASE (2) + CPABORT("gfn_type = 2 not yet supported") + CASE DEFAULT + CPABORT("Wrong gfn_type") + END SELECT + + filename = ADJUSTL(TRIM(parameter_file_path))//ADJUSTL(TRIM(spinpol_param_file_name)) + CALL parser_create(parser, filename, apply_preprocessing=.FALSE., para_env=para_env) + znum = 0 + param%wall = 0.0_dp + CALL get_ptable_info(element_symbol, znum) + DO + at_end = .FALSE. + CALL parser_get_next_line(parser, 1, at_end) + IF (at_end) EXIT + CALL parser_get_object(parser, zin) + IF (zin == znum) THEN + CALL parser_get_object(parser, param%wall(1, 1)) + CALL parser_get_object(parser, param%wall(1, 2)) + CALL parser_get_object(parser, param%wall(2, 2)) + CALL parser_get_object(parser, param%wall(1, 3)) + CALL parser_get_object(parser, param%wall(2, 3)) + CALL parser_get_object(parser, param%wall(3, 3)) + param%wall(2, 1) = param%wall(1, 2) + param%wall(3, 1) = param%wall(1, 3) + param%wall(3, 2) = param%wall(2, 3) + END IF + END DO + CALL parser_release(parser) + + END SUBROUTINE xtb_spinpol_init + +! ************************************************************************************************** +!> \brief ... +!> \param param ... +!> \param gfn_type ... +!> \param xtb_control ... +! ************************************************************************************************** + SUBROUTINE xtb_spinpol_ext(param, gfn_type, xtb_control) + TYPE(xtb_atom_type), POINTER :: param + INTEGER, INTENT(IN) :: gfn_type + TYPE(xtb_control_type), INTENT(IN), POINTER :: xtb_control + + INTEGER :: i + + SELECT CASE (gfn_type) + CASE (0) + CPABORT("gfn_type = 0: No spin polarisation possible!") + CASE (1) + ! OK + CASE (2) + CPABORT("gfn_type = 2 not yet supported") + CASE DEFAULT + CPABORT("Wrong gfn_type") + END SELECT + + IF (param%defined) THEN + IF (ASSOCIATED(xtb_control%spinpol_type)) THEN + DO i = 1, SIZE(xtb_control%spinpol_type) + IF (xtb_control%spinpol_type(i) == param%z) THEN + param%wall(1, 1) = xtb_control%spinpol_vals(1, i) + param%wall(1, 2) = xtb_control%spinpol_vals(2, i) + param%wall(2, 2) = xtb_control%spinpol_vals(3, i) + param%wall(1, 3) = xtb_control%spinpol_vals(4, i) + param%wall(2, 3) = xtb_control%spinpol_vals(5, i) + param%wall(3, 3) = xtb_control%spinpol_vals(6, i) + param%wall(2, 1) = param%wall(1, 2) + param%wall(3, 1) = param%wall(1, 3) + param%wall(3, 2) = param%wall(2, 3) + EXIT + END IF + END DO + END IF + END IF + + END SUBROUTINE xtb_spinpol_ext + ! ************************************************************************************************** !> \brief Read atom parameters for xTB Hamiltonian from input file !> \param param ... diff --git a/src/xtb_spinpol.F b/src/xtb_spinpol.F new file mode 100644 index 0000000000..5a6f431da1 --- /dev/null +++ b/src/xtb_spinpol.F @@ -0,0 +1,496 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2026 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Calculation of Spin Polarisation contributions in xTB +!> \author JGH +! ************************************************************************************************** +MODULE xtb_spinpol + USE atomic_kind_types, ONLY: atomic_kind_type,& + get_atomic_kind,& + get_atomic_kind_set + USE atprop_types, ONLY: atprop_type + USE cell_types, ONLY: cell_type + USE cp_control_types, ONLY: dft_control_type + USE cp_dbcsr_api, ONLY: dbcsr_add,& + dbcsr_get_block_p,& + dbcsr_iterator_blocks_left,& + dbcsr_iterator_next_block,& + dbcsr_iterator_start,& + dbcsr_iterator_stop,& + dbcsr_iterator_type,& + dbcsr_p_type, dbcsr_type + USE kinds, ONLY: dp + USE kpoint_types, ONLY: get_kpoint_info,& + kpoint_type + USE message_passing, ONLY: mp_para_env_type + USE mulliken, ONLY: ao_charges + USE particle_types, ONLY: particle_type + USE qs_energy_types, ONLY: qs_energy_type + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type + USE qs_force_types, ONLY: qs_force_type + USE qs_kind_types, ONLY: get_qs_kind,& + get_qs_kind_set,& + qs_kind_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 sap_kind_types, ONLY: sap_int_type + USE virial_types, ONLY: virial_type + 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_spinpol' + + PUBLIC :: build_xtb_spinpol + +CONTAINS + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param ks_matrix ... +!> \param matrix_p ... +!> \param energy ... +!> \param calculate_forces ... +!> \param just_energy ... +! ************************************************************************************************** + SUBROUTINE build_xtb_spinpol(qs_env, ks_matrix, matrix_p, energy, & + sap_int, calculate_forces, just_energy) + + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_p + TYPE(qs_energy_type), POINTER :: energy + TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int + LOGICAL, INTENT(in) :: calculate_forces, just_energy + + CHARACTER(len=*), PARAMETER :: routineN = 'build_xtb_spinpol' + + INTEGER :: atom_a, handle, iatom, ikind, is, na, ns, nspins, & + natom, natorb, nimg, nkind, nsgf, lmax, la, lb, lma, & + icol, irow, ia, ib, jkind, jatom, i, ic, nb, & + atom_i, atom_j, iac + INTEGER, DIMENSION(25) :: lao + INTEGER, DIMENSION(3) :: cellind + INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index + INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of + LOGICAL :: use_virial, defined, found + REAL(KIND=dp) :: espin, dr, fval, fi, fo + REAL(KIND=dp), DIMENSION(3) :: rij, fij + REAL(KIND=dp), DIMENSION(5) :: pal + REAL(KIND=dp), DIMENSION(3, 3) :: wall + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: docg + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, bocg, wab + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wabk + REAL(KIND=dp), DIMENSION(:, :), POINTER :: aksb, bksb, sblock, pamat, pbmat, dsblock + REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(atprop_type), POINTER :: atprop + TYPE(cell_type), POINTER :: cell + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p_kp, matrix_s, matrix_s_kp + TYPE(dbcsr_type), POINTER :: s_matrix + TYPE(dbcsr_iterator_type) :: iter + TYPE(dft_control_type), POINTER :: dft_control + TYPE(kpoint_type), POINTER :: kpoints + 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(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(neighbor_list_iterator_p_type), & + DIMENSION(:), POINTER :: nl_iterator + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: n_list + TYPE(virial_type), POINTER :: virial + TYPE(xtb_atom_type), POINTER :: xtb_kind + + CALL timeset(routineN, handle) + + energy%xtb_spinpol = 0.0_dp + + CALL get_qs_env(qs_env, dft_control=dft_control) + nspins = dft_control%nspins + nimg = dft_control%nimages + + IF (nspins == 2) THEN + + CALL get_qs_env(qs_env, & + qs_kind_set=qs_kind_set, & + particle_set=particle_set, & + atomic_kind_set=atomic_kind_set, & + cell=cell, & + virial=virial, & + atprop=atprop) + + CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & + kind_of=kind_of, & + atom_of_kind=atom_of_kind) + + use_virial = .FALSE. + IF (calculate_forces) THEN + use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) + END IF + + CALL get_qs_env(qs_env, nkind=nkind, natom=natom) + CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf) + CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env) + + ! expand parameters + ALLOCATE(wabk(nsgf, nsgf, nkind)) + wabk = 0.0_dp + DO ikind = 1, nkind + CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind) + CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, wall=wall) + DO ia = 1, natorb + la = lao(ia) + 1 + DO ib = 1, natorb + lb = lao(ib) + 1 + wabk(ia, ib, ikind) = wall(la, lb) + END DO + END DO + END DO + + ! Calculate charges + ALLOCATE (aocg(nsgf, natom), bocg(nsgf, natom)) + aocg = 0.0_dp + IF (nimg > 1) THEN + matrix_s_kp => matrix_s(:, :) + matrix_p_kp => matrix_p(1:1, :) + CALL ao_charges(matrix_p_kp, matrix_s_kp, aocg, para_env) + matrix_p_kp => matrix_p(2:2, :) + CALL ao_charges(matrix_p_kp, matrix_s_kp, bocg, para_env) + ELSE + s_matrix => matrix_s(1, 1)%matrix + p_matrix => matrix_p(1:1, 1) + CALL ao_charges(p_matrix, s_matrix, aocg, para_env) + p_matrix => matrix_p(2:2, 1) + CALL ao_charges(p_matrix, s_matrix, bocg, para_env) + END IF + + ! calculate energy + DO ikind = 1, nkind + CALL get_atomic_kind(atomic_kind_set(ikind), natom=na) + CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind) + CALL get_xtb_atom_param(xtb_kind, defined=defined, natorb=natorb) + IF (.NOT. defined .OR. natorb < 1) CYCLE + ALLOCATE(docg(natorb), wab(natorb, natorb)) + wab(1:natorb, 1:natorb) = wabk(1:natorb, 1:natorb, ikind) + DO iatom = 1, na + atom_a = atomic_kind_set(ikind)%atom_list(iatom) + docg = 0.0_dp + docg(1:natorb) = aocg(1:natorb, atom_a) - bocg(1:natorb, atom_a) + espin = 0.5_dp * DOT_PRODUCT(docg, MATMUL(wab, docg)) + energy%xtb_spinpol = energy%xtb_spinpol + espin + IF (atprop%energy) THEN + atprop%atecoul(iatom) = atprop%atecoul(iatom) + espin + END IF + END DO + DEALLOCATE(docg, wab) + END DO + + ! Forces and Virial + IF (calculate_forces) THEN + CALL get_qs_env(qs_env=qs_env, force=force) + NULLIFY (cell_to_index) + IF (nimg > 1) THEN + NULLIFY (kpoints) + CALL get_qs_env(qs_env=qs_env, kpoints=kpoints) + CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) + END IF + IF (nimg == 1) THEN + ! no k-points; all matrices have been transformed to periodic bsf + CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix) + DO WHILE (dbcsr_iterator_blocks_left(iter)) + CALL dbcsr_iterator_next_block(iter, irow, icol, sblock) + ikind = kind_of(irow) + atom_i = atom_of_kind(irow) + jkind = kind_of(icol) + atom_j = atom_of_kind(icol) + + CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, & + row=irow, col=icol, block=pamat, found=found) + CPASSERT(found) + CALL dbcsr_get_block_p(matrix=matrix_p(2, 1)%matrix, & + row=irow, col=icol, block=pbmat, found=found) + CPASSERT(found) + + na = SIZE(pamat,1) + nb = SIZE(pamat,2) + + fval = 1.0_dp + + DO i = 1, 3 + CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, & + row=irow, col=icol, block=dsblock, found=found) + CPASSERT(found) + + CALL fupdate(fi, pamat, pbmat, dsblock, na, nb, fval, & + wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), & + aocg(1:na,irow), aocg(1:nb,icol), bocg(1:na,irow), bocg(1:nb,icol)) + + force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi + force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi + END DO + + END DO + CALL dbcsr_iterator_stop(iter) + ! use dsint list + IF (use_virial) THEN + CPASSERT(ASSOCIATED(sap_int)) + DO ikind = 1, nkind + DO jkind = 1, nkind + iac = ikind + nkind*(jkind - 1) + IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE + DO ia = 1, sap_int(iac)%nalist + IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE + iatom = sap_int(iac)%alist(ia)%aatom + DO ic = 1, sap_int(iac)%alist(ia)%nclist + jatom = sap_int(iac)%alist(ia)%clist(ic)%catom + rij = sap_int(iac)%alist(ia)%clist(ic)%rac + dr = SQRT(SUM(rij(:)**2)) + IF (dr > 1.e-6_dp) THEN + dsint => sap_int(iac)%alist(ia)%clist(ic)%acint + icol = MAX(iatom, jatom) + irow = MIN(iatom, jatom) + CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, & + row=irow, col=icol, block=pamat, found=found) + CPASSERT(found) + CALL dbcsr_get_block_p(matrix=matrix_p(2, 1)%matrix, & + row=irow, col=icol, block=pbmat, found=found) + CPASSERT(found) + fval = 1.0_dp + IF (irow == iatom) fval = -1.0_dp + DO i = 1, 3 + CALL fupdate(fi, pamat, pbmat, dsint(:, :, i), na, nb, fval, & + wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), & + aocg(1:na,irow), aocg(1:nb,icol), & + bocg(1:na,irow), bocg(1:nb,icol)) + fij(i) = fi + END DO + fi = 1.0_dp + IF (iatom == jatom) fi = 0.5_dp + CALL virial_pair_force(virial%pv_virial, fi, fij, rij) + END IF + END DO + END DO + END DO + END DO + END IF + ELSE + NULLIFY (n_list) + CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) + CALL neighbor_list_iterator_create(nl_iterator, n_list) + 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=cellind) + + dr = SQRT(SUM(rij**2)) + IF (iatom == jatom .AND. dr < 1.0e-6_dp) CYCLE + + icol = MAX(iatom, jatom) + irow = MIN(iatom, jatom) + + ic = cell_to_index(cellind(1), cellind(2), cellind(3)) + CPASSERT(ic > 0) + + atom_i = atom_of_kind(iatom) + atom_j = atom_of_kind(jatom) + ! + CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, & + row=irow, col=icol, block=pamat, found=found) + CPASSERT(found) + CALL dbcsr_get_block_p(matrix=matrix_p(2, ic)%matrix, & + row=irow, col=icol, block=pbmat, found=found) + CPASSERT(found) + + na = SIZE(pamat,1) + nb = SIZE(pamat,2) + + fval = 1.0_dp + IF (irow == iatom) fval = -1.0_dp + + fij = 0.0_dp + DO i = 1, 3 + CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, & + row=irow, col=icol, block=dsblock, found=found) + CPASSERT(found) + + CALL fupdate(fi, pamat, pbmat, dsblock, na, nb, fval, & + wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), & + aocg(1:na,irow), aocg(1:nb,icol), bocg(1:na,irow), bocg(1:nb,icol)) + + force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi + force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi + fij(i) = fi + END DO + IF (use_virial) THEN + fi = 1.0_dp + IF (iatom == jatom) fi = 0.5_dp + CALL virial_pair_force(virial%pv_virial, fi, fij, rij) + END IF + + END DO + CALL neighbor_list_iterator_release(nl_iterator) + + END IF + END IF + + ! KS matrix + IF (.NOT. just_energy) THEN + IF (nimg > 1) THEN + CALL get_qs_env(qs_env=qs_env, kpoints=kpoints) + CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) + END IF + IF (nimg == 1) THEN + ! no k-points; all matrices have been transformed to periodic bsf + CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix) + DO WHILE (dbcsr_iterator_blocks_left(iter)) + CALL dbcsr_iterator_next_block(iter, irow, icol, sblock) + CALL dbcsr_get_block_p(matrix=ks_matrix(1, 1)%matrix, & + row=irow, col=icol, block=aksb, found=found) + CPASSERT(found) + CALL dbcsr_get_block_p(matrix=ks_matrix(2, 1)%matrix, & + row=irow, col=icol, block=bksb, found=found) + CPASSERT(found) + na = SIZE(aksb,1) + nb = SIZE(aksb,2) + ikind = kind_of(irow) + jkind = kind_of(icol) + fval = 0.5_dp + CALL ksupdate(aksb, bksb, sblock, na, nb, fval, & + wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), & + aocg(1:na,irow), aocg(1:nb,icol), bocg(1:na,irow), bocg(1:nb,icol)) + END DO + CALL dbcsr_iterator_stop(iter) + ELSE + CALL get_qs_env(qs_env=qs_env, sab_orb=n_list) + CALL neighbor_list_iterator_create(nl_iterator, n_list) + 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=cellind) + + icol = MAX(iatom, jatom) + irow = MIN(iatom, jatom) + + ic = cell_to_index(cellind(1), cellind(2), cellind(3)) + CPASSERT(ic > 0) + + ikind = kind_of(iatom) + jkind = kind_of(jatom) + + CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, & + row=irow, col=icol, block=sblock, found=found) + CPASSERT(found) + CALL dbcsr_get_block_p(matrix=ks_matrix(1, ic)%matrix, & + row=irow, col=icol, block=aksb, found=found) + CPASSERT(found) + CALL dbcsr_get_block_p(matrix=ks_matrix(2, ic)%matrix, & + row=irow, col=icol, block=bksb, found=found) + CPASSERT(found) + + na = SIZE(aksb,1) + nb = SIZE(aksb,2) + fval = 0.5_dp + CALL ksupdate(aksb, bksb, sblock, na, nb, fval, & + wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), & + aocg(1:na,irow), aocg(1:nb,icol), bocg(1:na,irow), bocg(1:nb,icol)) + END DO + CALL neighbor_list_iterator_release(nl_iterator) + END IF + + END IF + + DEALLOCATE (wabk) + DEALLOCATE (aocg, bocg) + END IF + + CALL timestop(handle) + + END SUBROUTINE build_xtb_spinpol + + SUBROUTINE ksupdate(aksb, bksb, sb, na, nb, fval, & + wabi, wabj, qai, qaj, qbi, qbj) + + REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aksb, bksb + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sb + INTEGER, INTENT(IN) :: na, nb + REAL(KIND=dp), INTENT(IN) :: fval + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: wabi + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: wabj + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: qai + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: qaj + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: qbi + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: qbj + + INTEGER :: ia, ib + REAL(KIND=dp), DIMENSION(na) :: dqa, wa + REAL(KIND=dp), DIMENSION(nb) :: dqb, wb + REAL(KIND=dp), DIMENSION(na, nb) :: wqab + + dqa = qai - qbi + dqb = qaj - qbj + wa = MATMUL(wabi, dqa) + wb = MATMUL(wabj, dqb) + DO ib=1,nb + DO ia=1,na + wqab(ia,ib) = fval * sb(ia,ib) * (wa(ia) + wb(ib)) + END DO + END DO + + aksb = aksb + wqab + bksb = bksb - wqab + + END SUBROUTINE ksupdate + + SUBROUTINE fupdate(fij, pa, pb, ds, na, nb, fval, & + wabi, wabj, qai, qaj, qbi, qbj) + + REAL(KIND=dp), INTENT(OUT) :: fij + REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: pa, pb + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ds + INTEGER, INTENT(IN) :: na, nb + REAL(KIND=dp), INTENT(IN) :: fval + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: wabi + REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: wabj + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: qai + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: qaj + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: qbi + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: qbj + + INTEGER :: ia, ib + REAL(KIND=dp), DIMENSION(na) :: dqa, wa, dpsa + REAL(KIND=dp), DIMENSION(nb) :: dqb, wb, dpsb + REAL(KIND=dp), DIMENSION(na, nb) :: dpab + + dqa = qai - qbi + dqb = qaj - qbj + wa = MATMUL(wabi, dqa) + wb = MATMUL(wabj, dqb) + dpab = pa - pb + dpsa = 0.0_dp + dpsb = 0.0_dp + DO ib=1,nb + DO ia=1,na + dpsa(ia) = dpsa(ia) + dpab(ia, ib)*ds(ia, ib) + dpsb(ib) = dpsb(ib) + dpab(ia, ib)*ds(ia, ib) + END DO + END DO + + fij = SUM(wa*dpsa) + SUM(wb*dpsb) + + END SUBROUTINE fupdate + +END MODULE xtb_spinpol diff --git a/src/xtb_types.F b/src/xtb_types.F index 9ec859eae6..5075595868 100644 --- a/src/xtb_types.F +++ b/src/xtb_types.F @@ -69,6 +69,7 @@ MODULE xtb_types REAL(KIND=dp), DIMENSION(5) :: kappa = -1.0_dp REAL(KIND=dp), DIMENSION(5) :: hen = -1.0_dp REAL(KIND=dp), DIMENSION(5) :: zeta = -1.0_dp + REAL(KIND=dp), DIMENSION(3, 3) :: wall = -1.0_dp ! spin polarisation ! gfn0 params REAL(KIND=dp) :: en = -1.0_dp REAL(KIND=dp) :: kqat2 = -1.0_dp @@ -127,6 +128,7 @@ CONTAINS xtb_parameter%occupation = 0 xtb_parameter%kpoly = 0.0_dp xtb_parameter%kappa = 0.0_dp + xtb_parameter%wall = 0.0_dp xtb_parameter%hen = 0.0_dp xtb_parameter%zeta = 0.0_dp xtb_parameter%en = 0.0_dp @@ -180,6 +182,7 @@ CONTAINS !> \param lval ... !> \param kpoly ... !> \param kappa ... +!> \param wall ... !> \param hen ... !> \param zeta ... !> \param xi ... @@ -195,7 +198,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, & rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, & - hen, zeta, xi, kappa0, alpg, occupation, electronegativity, chmax, & + wall, hen, zeta, xi, kappa0, alpg, occupation, electronegativity, chmax, & en, kqat2, kcn, kq) TYPE(xtb_atom_type), POINTER :: xtb_parameter @@ -210,7 +213,10 @@ CONTAINS REAL(KIND=dp), INTENT(OUT), OPTIONAL :: rcut, rcov, kx, eta, xgamma, alpha, zneff INTEGER, INTENT(OUT), OPTIONAL :: nshell INTEGER, DIMENSION(5), INTENT(OUT), OPTIONAL :: nval, lval - REAL(KIND=dp), DIMENSION(5), INTENT(OUT), OPTIONAL :: kpoly, kappa, hen, zeta + REAL(KIND=dp), DIMENSION(5), INTENT(OUT), OPTIONAL :: kpoly, kappa + REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), & + OPTIONAL :: wall + REAL(KIND=dp), DIMENSION(5), INTENT(OUT), OPTIONAL :: hen, zeta REAL(KIND=dp), INTENT(OUT), OPTIONAL :: xi, kappa0, alpg INTEGER, DIMENSION(5), INTENT(OUT), OPTIONAL :: occupation REAL(KIND=dp), INTENT(OUT), OPTIONAL :: electronegativity, chmax, en, kqat2 @@ -243,6 +249,7 @@ CONTAINS IF (PRESENT(occupation)) occupation = xtb_parameter%occupation IF (PRESENT(kpoly)) kpoly = xtb_parameter%kpoly IF (PRESENT(kappa)) kappa = xtb_parameter%kappa + IF (PRESENT(wall)) wall(1:3, 1:3) = xtb_parameter%wall(1:3, 1:3) IF (PRESENT(hen)) hen = xtb_parameter%hen IF (PRESENT(zeta)) zeta = xtb_parameter%zeta IF (PRESENT(chmax)) chmax = xtb_parameter%chmax @@ -280,6 +287,7 @@ CONTAINS !> \param lval ... !> \param kpoly ... !> \param kappa ... +!> \param wall ... !> \param hen ... !> \param zeta ... !> \param xi ... @@ -295,7 +303,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE set_xtb_atom_param(xtb_parameter, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, & rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, & - hen, zeta, xi, kappa0, alpg, electronegativity, occupation, chmax, & + wall, hen, zeta, xi, kappa0, alpg, electronegativity, occupation, chmax, & en, kqat2, kcn, kq) TYPE(xtb_atom_type), POINTER :: xtb_parameter @@ -309,7 +317,10 @@ CONTAINS REAL(KIND=dp), INTENT(IN), OPTIONAL :: rcut, rcov, kx, eta, xgamma, alpha, zneff INTEGER, INTENT(IN), OPTIONAL :: nshell INTEGER, DIMENSION(5), INTENT(IN), OPTIONAL :: nval, lval - REAL(KIND=dp), DIMENSION(5), INTENT(IN), OPTIONAL :: kpoly, kappa, hen, zeta + REAL(KIND=dp), DIMENSION(5), INTENT(IN), OPTIONAL :: kpoly, kappa + REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), & + OPTIONAL :: wall + REAL(KIND=dp), DIMENSION(5), INTENT(IN), OPTIONAL :: hen, zeta REAL(KIND=dp), INTENT(IN), OPTIONAL :: xi, kappa0, alpg, electronegativity INTEGER, DIMENSION(5), INTENT(IN), OPTIONAL :: occupation REAL(KIND=dp), INTENT(IN), OPTIONAL :: chmax, en, kqat2 @@ -341,6 +352,7 @@ CONTAINS IF (PRESENT(occupation)) xtb_parameter%occupation = occupation IF (PRESENT(kpoly)) xtb_parameter%kpoly = kpoly IF (PRESENT(kappa)) xtb_parameter%kappa = kappa + IF (PRESENT(wall)) xtb_parameter%wall(1:3, 1:3) = wall(1:3, 1:3) IF (PRESENT(hen)) xtb_parameter%hen = hen IF (PRESENT(zeta)) xtb_parameter%zeta = zeta IF (PRESENT(chmax)) xtb_parameter%chmax = chmax @@ -370,9 +382,10 @@ CONTAINS CHARACTER(LEN=default_string_length) :: aname, bb INTEGER :: i, io_unit, m, natorb, nshell INTEGER, DIMENSION(5) :: lval, nval, occupation - LOGICAL :: defined + LOGICAL :: defined, have_sp REAL(dp) :: zeff REAL(KIND=dp) :: alpha, en, eta, xgamma, zneff + REAL(KIND=dp), DIMENSION(3, 3) :: wall REAL(KIND=dp), DIMENSION(5) :: hen, kappa, kpoly, zeta TYPE(cp_logger_type), POINTER :: logger @@ -394,6 +407,10 @@ CONTAINS CALL get_xtb_atom_param(xtb_parameter, nshell=nshell, lval=lval, nval=nval, occupation=occupation) CALL get_xtb_atom_param(xtb_parameter, kpoly=kpoly, kappa=kappa, hen=hen, zeta=zeta) CALL get_xtb_atom_param(xtb_parameter, electronegativity=en, xgamma=xgamma, eta=eta, alpha=alpha, zneff=zneff) + wall = 0.0_dp + CALL get_xtb_atom_param(xtb_parameter, wall=wall) + have_sp = .FALSE. + IF (SUM(ABS(wall)) /= 0.0_dp) have_sp = .TRUE. bb = " " WRITE (UNIT=io_unit, FMT="(/,A,T67,A14)") " xTB parameters: ", TRIM(aname) @@ -413,6 +430,10 @@ CONTAINS (kappa(i), i=1, nshell) WRITE (UNIT=io_unit, FMT="(T16,A,T71,F10.3)") "3rd Order constant", xgamma WRITE (UNIT=io_unit, FMT="(T16,A,T61,2F10.3)") "Repulsion potential [Z,alpha]", zneff, alpha + IF (have_sp) THEN + WRITE (UNIT=io_unit, FMT="(T16,A,T51,3F10.4)") "Spin Polarisation Wss sp pp", wall(1, 1), wall(1, 2), wall(2, 2) + WRITE (UNIT=io_unit, FMT="(T16,A,T51,3F10.4)") " Wsd pd dd", wall(1, 3), wall(2, 3), wall(3, 3) + END IF ELSE WRITE (UNIT=io_unit, FMT="(T55,A)") "Parameters are not defined" END IF