mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-29 06:35:28 -04:00
xTB spin polarisation Hamiltonian native implementation
This commit is contained in:
Juerg Hutter
parent
3d5583cc72
commit
9838f79063
15 changed files with 803 additions and 13 deletions
80
data/xTB_spin_pol_params
Normal file
80
data/xTB_spin_pol_params
Normal file
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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., &
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 "
|
||||
|
|
|
|||
|
|
@ -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, &
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 ...
|
||||
|
|
|
|||
496
src/xtb_spinpol.F
Normal file
496
src/xtb_spinpol.F
Normal file
|
|
@ -0,0 +1,496 @@
|
|||
!--------------------------------------------------------------------------------------------------!
|
||||
! CP2K: A general program to perform molecular dynamics simulations !
|
||||
! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
|
||||
! !
|
||||
! 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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue