comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Bug fixes, consistent parameters, disabled CN term (forces not debugged)

Browse source
This commit is contained in:
Juerg Hutter 2018-12-17 14:28:45 +01:00
parent 2c6ccf5aa7
commit 841ab0f2eb
10 changed files with 424 additions and 263 deletions
Download patch file Download diff file Expand all files Collapse all files

10
data/xTB_parameters
View file

@ -160,7 +160,7 @@
4d 14.328052 0.0 -2.305768 0.853415
La 0.495969 1.500000 1.311200 48.337542 5d -44.208463 -6.396752 -5.872226 3.000000
6s -3.988102 0.0 -6.500000 1.492677
6p 40.847293 -1.500000 -0.727921 1.350000
6p 40.847293 -1.500000 -0.727921 1.350000
Ce 0.350000 1.200000 0.839861 30.638143 5d -36.440945 -5.245538 -5.032003 3.000000
6s 6.148475 0.0 -6.275363 1.553483
6p 42.873822 -1.500000 0.291196 1.380859
@ -201,8 +201,8 @@
6s 26.045692 0.0 -6.224540 1.857761
6p 42.541738 -1.500000 -0.301351 1.437991
Lu 0.249975 1.200000 0.936321 76.041625 5d -30.990414 -2.895435 -3.900750 2.899993
6s 27.703793 0.0 -6.220305 1.883118
6p 42.514065 -1.500000 -0.350730 1.442752
6s 27.703793 0.0 -6.220305 1.883118
6p 42.514065 -1.500000 -0.350730 1.442752
Hf 0.269977 0.847011 0.853744 55.222897 5d -21.116286 -1.485678 -4.360558 2.466693
6s 15.014122 0.0 -5.910623 2.039390
6p 22.898249 -1.500000 -2.814338 1.450000
@ -242,5 +242,5 @@
6p -12.804436 -3.402676 -10.031093 2.392024
5d 16.836546 0.0 -0.852571 1.380239
Rn 0.798170 -0.838429 0.998641 86.000000 6s -22.139701 0.0 -18.381647 3.520683
6p -20.539955 -2.380762 -10.236606 2.535389
5d 17.249637 0.0 -0.973687 1.418875
6p -20.539955 -2.380762 -10.236606 2.535389
5d 17.249637 0.0 -0.973687 1.418875

1
src/common/mathlib.F
View file

@ -394,6 +394,7 @@ CONTAINS
! Diagonalize the matrix a
info = 0
IF (divide_and_conquer) THEN
CALL dsyevd("V", "U", n, a, n, eigval, work, lwork, iwork, liwork, info)
ELSE

9
src/mulliken.F
View file

@ -709,7 +709,6 @@ CONTAINS
INTEGER :: blk, handle, i, iblock_col, iblock_row, &
ispin, j, nspin
LOGICAL :: found
REAL(kind=dp) :: mult
REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, s_block
TYPE(dbcsr_iterator_type) :: iter
@ -730,18 +729,12 @@ CONTAINS
IF (.NOT. (ASSOCIATED(s_block) .AND. ASSOCIATED(p_block))) CYCLE
IF (iblock_row == iatom) THEN
IF (iblock_row == iblock_col) THEN
mult = 1.0_dp ! avoid double counting of diagonal blocks
ELSE
mult = 2.0_dp
ENDIF
DO j = 1, SIZE(p_block, 2)
DO i = 1, SIZE(p_block, 1)
charges(i) = charges(i)+mult*p_block(i, j)*s_block(i, j)
charges(i) = charges(i)+p_block(i, j)*s_block(i, j)
END DO
END DO
ELSEIF (iblock_col == iatom) THEN
mult = 2.0_dp
DO i = 1, SIZE(p_block, 1)
DO j = 1, SIZE(p_block, 2)
charges(j) = charges(j)+p_block(i, j)*s_block(i, j)

1
src/qs_diis.F
View file

@ -441,6 +441,7 @@ CONTAINS
! Solve the linear DIIS equation system
ev(1:nb1) = 0.0_dp
CALL diamat_all(b(1:nb1, 1:nb1), ev(1:nb1))
a(1:nb1, 1:nb1) = b(1:nb1, 1:nb1)

4
src/qs_dispersion_pairpot.F
View file

@ -1069,6 +1069,8 @@ CONTAINS
DO iatom = 1, natom
ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
END DO
ELSE
ALLOCATE (dcnum(1))
END IF
CALL d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, ghost, floating, atomnumber, &
@ -1657,6 +1659,8 @@ CONTAINS
DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
END DO
DEALLOCATE (dcnum)
ELSE
DEALLOCATE (dcnum)
END IF
END IF

12
src/qs_environment.F
View file

@ -660,9 +660,8 @@ CONTAINS
NULLIFY (tmp_basis_set)
CALL init_xtb_basis(qs_kind%xtb_parameter, tmp_basis_set, ngauss)
CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, "ORB")
! cutoff radius
CALL get_gto_basis_set(tmp_basis_set, kind_radius=qs_kind%xtb_parameter%rcut)
! potential
zeff_correction = 0.0_dp
CALL init_potential(qs_kind%all_potential, itype="BARE", &
zeff=qs_kind%xtb_parameter%zeff, zeff_correction=zeff_correction)
qs_kind%xtb_parameter%zeff = qs_kind%xtb_parameter%zeff-zeff_correction
@ -813,6 +812,15 @@ CONTAINS
! *** Initialize the atomic interaction radii ***
CALL init_interaction_radii(dft_control%qs_control, atomic_kind_set, qs_kind_set)
!
IF (dft_control%qs_control%method_id == do_method_xtb) THEN
! cutoff radius
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind, basis_set=tmp_basis_set)
CALL get_gto_basis_set(tmp_basis_set, kind_radius=qs_kind%xtb_parameter%rcut)
END DO
END IF
IF (.NOT. be_silent) THEN
CALL write_pgf_orb_radii("orb", atomic_kind_set, qs_kind_set, subsys_section)

95
src/xtb_coulomb.F
View file

@ -127,10 +127,9 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'build_xtb_coulomb', &
routineP = moduleN//':'//routineN
INTEGER :: atom_i, atom_j, blk, ewald_type, handle, &
i, ia, iatom, ic, icol, ikind, img, &
irow, is, j, jatom, jkind, natom, ni, &
nimg, nj, nkind, nmat
INTEGER :: atom_i, atom_j, blk, ewald_type, handle, i, ia, iatom, ic, icol, ikind, img, &
irow, is, j, jatom, jkind, la, lb, natom, ni, nimg, nj, nkind, nmat, za, zb
INTEGER, DIMENSION(25) :: laoa, laob
INTEGER, DIMENSION(3) :: cellind, periodic
INTEGER, DIMENSION(:), POINTER :: atom_of_kind, kind_of
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
@ -162,13 +161,14 @@ CONTAINS
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(virial_type), POINTER :: virial
TYPE(xtb_atom_type), POINTER :: xtb_kind
TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind
CALL timeset(routineN, handle)
NULLIFY (gamma_matrix, matrix_p, matrix_s, virial, atprop, dft_control)
CALL get_qs_env(qs_env, &
qs_kind_set=qs_kind_set, &
gamma_matrix=gamma_matrix, &
particle_set=particle_set, &
cell=cell, &
@ -208,6 +208,10 @@ CONTAINS
nj = SIZE(gblock, 2)
gcint(1:ni) = MATMUL(gblock, charges(jatom, 1:nj))
gchrg(iatom, 1:ni, 1) = gchrg(iatom, 1:ni, 1)+gcint(1:ni)
IF (iatom /= jatom) THEN
gcint(1:nj) = MATMUL(charges(iatom, 1:ni), gblock)
gchrg(jatom, 1:nj, 1) = gchrg(jatom, 1:nj, 1)+gcint(1:nj)
END IF
IF (calculate_forces) THEN
DO i = 2, 4
NULLIFY (gblock)
@ -217,13 +221,15 @@ CONTAINS
nj = SIZE(gblock, 2)
gcint(1:ni) = MATMUL(gblock, charges(jatom, 1:nj))
gchrg(iatom, 1:ni, i) = gchrg(iatom, 1:ni, i)+gcint(1:ni)
IF (iatom /= jatom) THEN
gcint(1:nj) = MATMUL(charges(iatom, 1:ni), gblock)
gchrg(jatom, 1:nj, i) = gchrg(jatom, 1:nj, i)-gcint(1:nj)
END IF
END DO
END IF
ENDDO
CALL dbcsr_iterator_stop(iter)
energy%hartree = energy%hartree+ecsr
IF (calculate_forces .AND. use_virial) THEN
CALL dbcsr_iterator_start(iter, gamma_matrix(1)%matrix)
DO WHILE (dbcsr_iterator_blocks_left(iter))
@ -285,12 +291,14 @@ CONTAINS
rij = particle_set(iatom)%r-particle_set(jatom)%r
rij = pbc(rij, cell)
dr = SQRT(SUM(rij(:)**2))
gmcharge(iatom, 1) = gmcharge(iatom, 1)+mcharge(jatom)/dr
gmcharge(jatom, 1) = gmcharge(jatom, 1)+mcharge(iatom)/dr
DO i = 2, nmat
gmcharge(iatom, i) = gmcharge(iatom, i)+rij(i-1)*mcharge(jatom)/dr**3
gmcharge(jatom, i) = gmcharge(jatom, i)-rij(i-1)*mcharge(iatom)/dr**3
END DO
IF (dr > 1.e-6_dp) THEN
gmcharge(iatom, 1) = gmcharge(iatom, 1)+mcharge(jatom)/dr
gmcharge(jatom, 1) = gmcharge(jatom, 1)+mcharge(iatom)/dr
DO i = 2, nmat
gmcharge(iatom, i) = gmcharge(iatom, i)+rij(i-1)*mcharge(jatom)/dr**3
gmcharge(jatom, i) = gmcharge(jatom, i)-rij(i-1)*mcharge(iatom)/dr**3
END DO
END IF
END DO
END DO
END DO
@ -300,7 +308,6 @@ CONTAINS
! global sum of gamma*p arrays
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
force=force, para_env=para_env)
CALL mp_sum(gmcharge(:, 1), para_env%group)
CALL mp_sum(gchrg(:, :, 1), para_env%group)
@ -322,7 +329,8 @@ CONTAINS
DO iatom = 1, natom
ikind = kind_of(iatom)
CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
CALL get_xtb_atom_param(xtb_kind, nshell=ni)
CALL get_xtb_atom_param(xtb_kind, lmax=ni)
ni = ni+1
ecsr = ecsr+SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, 1))
END DO
energy%hartree = energy%hartree+0.5_dp*ecsr
@ -331,7 +339,8 @@ CONTAINS
CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
DO ikind = 1, SIZE(local_particles%n_el)
CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
CALL get_xtb_atom_param(xtb_kind, nshell=ni)
CALL get_xtb_atom_param(xtb_kind, lmax=ni)
ni = ni+1
DO ia = 1, local_particles%n_el(ikind)
iatom = local_particles%list(ikind)%array(ia)
atprop%atecoul(iatom) = atprop%atecoul(iatom)+ &
@ -347,10 +356,18 @@ CONTAINS
ikind = kind_of(iatom)
atom_i = atom_of_kind(iatom)
CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
CALL get_xtb_atom_param(xtb_kind, nshell=ni)
CALL get_xtb_atom_param(xtb_kind, lmax=ni)
! short range
ni = ni+1
DO i = 1, 3
fij(i) = SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, i+1))+ &
gmcharge(iatom, i+1)*mcharge(iatom)
fij(i) = SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, i+1))
END DO
force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i)-fij(1)
force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i)-fij(2)
force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i)-fij(3)
! long range
DO i = 1, 3
fij(i) = gmcharge(iatom, i+1)*mcharge(iatom)
END DO
force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i)-fij(1)
force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i)-fij(2)
@ -398,12 +415,27 @@ CONTAINS
CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
row=irow, col=icol, block=sblock, found=found)
CPASSERT(found)
! atomic parameters
CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
ni = SIZE(sblock, 1)
nj = SIZE(sblock, 2)
ALLOCATE (gcij(ni, nj))
DO i = 1, ni
DO j = 1, nj
gcij(i, j) = 0.5_dp*(gchrg(iatom, i, 1)+gchrg(jatom, j, 1))
IF (irow == iatom) THEN
la = laoa(i)+1
lb = laob(j)+1
gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1)+gchrg(jatom, lb, 1))
ELSE
la = laoa(j)+1
lb = laob(i)+1
gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1)+gchrg(jatom, lb, 1))
END IF
END DO
END DO
gmij = 0.5_dp*(gmcharge(iatom, 1)+gmcharge(jatom, 1))
@ -412,7 +444,8 @@ CONTAINS
CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
row=irow, col=icol, block=ksblock, found=found)
CPASSERT(found)
ksblock = ksblock-gcij*sblock-gmij*sblock
ksblock = ksblock-gcij*sblock
ksblock = ksblock-gmij*sblock
END DO
IF (calculate_forces) THEN
@ -420,7 +453,10 @@ CONTAINS
atom_i = atom_of_kind(iatom)
jkind = kind_of(jatom)
atom_j = atom_of_kind(jatom)
IF (irow == jatom) gmij = -gmij
IF (irow == jatom) THEN
gmij = -gmij
gcij = -gcij
END IF
NULLIFY (pblock)
CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
row=irow, col=icol, block=pblock, found=found)
@ -430,10 +466,17 @@ CONTAINS
CALL dbcsr_get_block_p(matrix=matrix_s(1+i, ic)%matrix, &
row=irow, col=icol, block=dsblock, found=found)
CPASSERT(found)
fi = -2.0_dp*(SUM(pblock*dsblock*gcij)+gmij*SUM(pblock*dsblock))
fij(i) = 0.0_dp
! short range
fi = -2.0_dp*SUM(pblock*dsblock*gcij)
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
fij(i) = fij(i)+fi
! long range
fi = -2.0_dp*gmij*SUM(pblock*dsblock)
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) = fij(i)+fi
END DO
IF (use_virial) THEN
CALL virial_pair_force(virial%pv_virial, 1._dp, fij, rij)
@ -463,8 +506,8 @@ CONTAINS
CALL get_xtb_atom_param(xtb_kind, xgamma=xgamma(ikind), zeff=zeffk(ikind))
END DO
! Diagonal 3rd order correction (DFTB3)
!deb CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, &
!deb calculate_forces, just_energy)
CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, &
calculate_forces, just_energy)
DEALLOCATE (zeffk, xgamma)
DEALLOCATE (gmcharge, gchrg, atom_of_kind, kind_of)

500
src/xtb_matrices.F
View file

@ -30,30 +30,26 @@ MODULE xtb_matrices
cp_print_key_unit_nr
USE cp_para_types, ONLY: cp_para_env_type
USE dbcsr_api, ONLY: &
convert_offsets_to_sizes, dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, &
dbcsr_finalize, dbcsr_get_block_p, dbcsr_multiply, dbcsr_p_type, dbcsr_type, &
dbcsr_type_antisymmetric, dbcsr_type_symmetric
dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_finalize, &
dbcsr_get_block_p, dbcsr_multiply, dbcsr_p_type, dbcsr_type, dbcsr_type_antisymmetric, &
dbcsr_type_symmetric
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: default_string_length,&
dp
USE kinds, ONLY: dp
USE kpoint_types, ONLY: get_kpoint_info,&
kpoint_type
USE message_passing, ONLY: mp_sum
USE mulliken, ONLY: ao_charges
USE particle_methods, ONLY: get_particle_set
USE particle_types, ONLY: particle_type
USE qs_dispersion_pairpot, ONLY: d3_cnumber,&
dcnum_type
USE qs_dispersion_types, ONLY: qs_atom_dispersion_type,&
qs_dispersion_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_force_types, ONLY: qs_force_type
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type
USE qs_ks_types, ONLY: get_ks_env,&
qs_ks_env_type,&
@ -107,23 +103,25 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'build_xtb_matrices', &
routineP = moduleN//':'//routineN
INTEGER :: after, atom_a, atom_b, handle, i, iatom, ic, icol, ikind, img, ir, irow, iw, j, &
jatom, jkind, la, lb, n1, n2, na, natom, natorb_a, natorb_b, nb, nderivatives, nimg, &
nkind, nmat, nsa, nsb, nshell, za, zb
INTEGER :: after, atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, &
iw, j, jatom, jkind, katom, kkind, la, lb, lmaxa, lmaxb, n1, n2, na, natom, natorb_a, &
natorb_b, nb, nderivatives, nimg, nkind, nmat, nsa, nsb, nshell, za, zb
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, atomnumber, kind_of
INTEGER, DIMENSION(25) :: laoa, laob, lval, naoa, naob
INTEGER, DIMENSION(3) :: cell
INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
LOGICAL :: defined, floating_a, found, ghost_a, &
omit_headers, use_virial
LOGICAL :: defined, diagblock, floating_a, found, &
ghost_a, h2scase, omit_headers, &
use_virial
LOGICAL, ALLOCATABLE, DIMENSION(:) :: floating, ghost
REAL(KIND=dp) :: alphaa, alphab, derepij, dfij, dr, ena, enb, erep, erepij, etaa, etab, f0, &
f1, fen, fij, foab, k2sh, kab, kcnd, kcnp, kcns, kd, ken, kf, kg, kp, ks, ksp, kx2, kxr, &
rcova, rcovab, rcovb, rcut, rcuta, rcutb, rrab, xgammaa, xgammab, zneffa, zneffb
REAL(KIND=dp) :: alphaa, alphab, derepij, dfij, dfp, dhij, dr, drk, drx, ena, enb, erep, &
erepij, etaa, etab, f0, f1, fen, fhua, fhub, fij, foab, hij, k2sh, kab, kcnd, kcnp, kcns, &
kd, ken, kf, kg, kia, kjb, kp, ks, ksp, kx2, kxr, rcova, rcovab, rcovb, rcut, rcuta, &
rcutb, rrab, xgammaa, xgammab, zneffa, zneffb
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dgamma, huckel, kabset
REAL(KIND=dp), DIMENSION(0:3) :: kcnl, kl
REAL(KIND=dp), DIMENSION(3) :: force_ab, force_rr, rij
REAL(KIND=dp), DIMENSION(3) :: fdik, force_ab, force_rr, rij, rik
REAL(KIND=dp), DIMENSION(5) :: dpia, dpib, hena, henb, kappaa, kappab, &
kpolya, kpolyb, pia, pib
REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, gblock, pblock, sblock
@ -141,7 +139,6 @@ CONTAINS
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_orb
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_atom_dispersion_type), POINTER :: disp_a
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
@ -176,10 +173,6 @@ CONTAINS
CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
nderivatives = 0
IF (calculate_forces) nderivatives = 1
CALL setup_matrices2(qs_env, nderivatives, nimg, matrix_s, "OVERLAP", sab_orb)
CALL setup_matrices2(qs_env, 0, nimg, matrix_h, "CORE HAMILTONIAN", sab_orb)
CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
! global parameters (Table 2 from Ref.)
ks = xtb_control%ks
@ -212,10 +205,8 @@ CONTAINS
IF (calculate_forces) THEN
NULLIFY (rho, force, matrix_w)
CALL get_qs_env(qs_env=qs_env, &
rho=rho, &
matrix_w_kp=matrix_w, &
virial=virial, &
force=force)
rho=rho, matrix_w_kp=matrix_w, &
virial=virial, force=force)
CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
IF (SIZE(matrix_p, 1) == 2) THEN
@ -257,6 +248,17 @@ CONTAINS
basis_type_b="ORB", &
sab_nl=sab_orb)
END IF
CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
! initialize H matrix
CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
DO img = 1, nimg
ALLOCATE (matrix_h(1, img)%matrix)
CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
name="HAMILTONIAN MATRIX")
CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
END DO
CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
! Calculate coordination numbers
! needed for effective atomic energy levels (Eq. 12)
@ -269,18 +271,19 @@ CONTAINS
DO iatom = 1, natom
ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
END DO
ELSE
ALLOCATE (dcnum(1))
END IF
nkind = SIZE(atomic_kind_set)
ALLOCATE (ghost(nkind), floating(nkind), atomnumber(nkind))
DO ikind = 1, nkind
CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost_a, floating=floating_a)
ghost(ikind) = ghost_a
floating(ikind) = floating_a
atomnumber(ikind) = za
END DO
CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
CALL d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, ghost, floating, atomnumber, &
calculate_forces, .FALSE.)
@ -297,6 +300,9 @@ CONTAINS
kl(1) = kp
kl(2) = kd
kl(3) = 0.0_dp
!deb
kcnl = 0.0_dp
!deb
ALLOCATE (huckel(5, natom), kind_of(natom))
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
DO iatom = 1, natom
@ -308,7 +314,6 @@ CONTAINS
huckel(i, iatom) = hena(i)*(1._dp+kcnl(lval(i))*cnumbers(iatom))
END DO
END DO
DEALLOCATE (kind_of)
! Calculate KAB parameters
! precalculate the KAB for all kind pairs (Table 2)
@ -341,11 +346,11 @@ CONTAINS
! atomic parameters
CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, xgamma=xgammaa, &
nshell=nsa, alpha=alphaa, zneff=zneffa, kpoly=kpolya, kappa=kappaa, hen=hena, &
electronegativity=ena)
lmax=lmaxa, nshell=nsa, alpha=alphaa, zneff=zneffa, kpoly=kpolya, &
kappa=kappaa, hen=hena, electronegativity=ena)
CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, xgamma=xgammab, &
nshell=nsb, alpha=alphab, zneff=zneffb, kpoly=kpolyb, kappa=kappab, hen=henb, &
electronegativity=enb)
lmax=lmaxb, nshell=nsb, alpha=alphab, zneff=zneffb, kpoly=kpolyb, &
kappa=kappab, hen=henb, electronegativity=enb)
IF (nimg == 1) THEN
ic = 1
@ -385,103 +390,256 @@ CONTAINS
END IF
! Calculate Pi = Pia * Pib (Eq. 11)
! factor 0.01 from original Grimme code, not in paper?
rcovab = rcova+rcovb
rrab = SQRT(dr/rcovab)
DO i = 1, nsa
pia(i) = 1._dp+0.01_dp*kpolya(i)*rrab
dpia(i) = 0.5_dp*0.01_dp*kpolya(i)/rrab
pia(i) = 1._dp+kpolya(i)*rrab
END DO
DO i = 1, nsb
pib(i) = 1._dp+0.01_dp*kpolyb(i)*rrab
dpib(i) = 0.5_dp*0.01_dp*kpolyb(i)/rrab
pib(i) = 1._dp+kpolyb(i)*rrab
END DO
IF (calculate_forces) THEN
IF (dr > 1.e-6_dp) THEN
drx = 0.5_dp/rrab/rcovab
ELSE
drx = 0.0_dp
END IF
dpia(1:nsa) = drx*kpolya(1:nsa)
dpib(1:nsb) = drx*kpolyb(1:nsb)
END IF
IF (iatom == jatom .AND. dr < 0.001_dp) THEN
! diagonal block
! we use <a|H0|a> = Hla
! diagonal block
diagblock = .FALSE.
IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
!
! Eq. 10
!
! get KAB
kab = kabset(ikind, jkind)
! get Fen = (1+ken*deltaEN^2)
fen = 1.0_dp+ken*(ena-enb)**2
!
DO j = 1, natorb_b
lb = laob(j)
nb = naob(j)
DO i = 1, natorb_a
la = laoa(i)
na = naoa(i)
fblock(i, i) = fblock(i, i)+huckel(na, iatom)
IF (diagblock .AND. i == j) THEN
fblock(i, i) = fblock(i, i)+huckel(na, iatom)
ELSE
kia = kl(la)
kjb = kl(lb)
h2scase = .FALSE.
IF (zb == 1 .AND. nb == 2) THEN
kjb = k2sh
h2scase = .TRUE.
END IF
IF (za == 1 .AND. na == 2) THEN
kia = k2sh
h2scase = .TRUE.
END IF
hij = 0.5_dp*(huckel(na, iatom)+huckel(nb, jatom))*pia(na)*pib(nb)
IF (h2scase) THEN
fij = 0.5_dp*(kia+kjb)*hij
ELSE
IF ((la == 0 .AND. lb == 1) .OR. (la == 1 .AND. lb == 0)) THEN
fij = ksp*hij*kab*fen
ELSE
fij = 0.5_dp*(kia+kjb)*hij*kab*fen
END IF
END IF
IF (iatom <= jatom) THEN
fblock(i, j) = fblock(i, j)+fij*sblock(i, j)
ELSE
fblock(j, i) = fblock(j, i)+fij*sblock(j, i)
END IF
END IF
END DO
IF (calculate_forces) THEN
!qxtb
! Derivative wrt coordination number is missing here
!qxtb
END IF
ELSE
! off-diagonal block
! Eq. 10
!
! get KAB
kab = kabset(ikind, jkind)
! get Fen = (1+ken*deltaEN^2)
! original code reads: (1-0.01*ken*den^2)
! fen = 1.0_dp+ken*(ena-enb)**2
fen = 1.0_dp - 0.01_dp*ken*(ena-enb)**2
END DO
IF (calculate_forces) THEN
f0 = 1.0_dp
IF (irow == iatom) f0 = -1.0_dp
! Derivative wrt coordination number
fhua = 0.0_dp
fhub = 0.0_dp
DO j = 1, natorb_b
lb = laob(j)
nb = naob(j)
DO i = 1, natorb_a
la = laoa(i)
na = naoa(i)
fij = 0.5_dp*(kl(la)+kl(lb))*0.5_dp*(huckel(na, iatom)+huckel(nb, jatom))*pia(na)*pib(nb)
IF (iatom <= jatom) THEN
fblock(i, j) = fblock(i, j)+fij*sblock(i, j)*kab*fen
IF (diagblock .AND. i == j) THEN
fhua = fhua+pblock(i, i)*kcnl(la)*0.5_dp
fhub = fhub+pblock(i, i)*kcnl(lb)*0.5_dp
ELSE
fblock(j, i) = fblock(j, i)+fij*sblock(j, i)*kab*fen
kia = kl(la)
kjb = kl(lb)
h2scase = .FALSE.
IF (zb == 1 .AND. nb == 2) THEN
kjb = k2sh
h2scase = .TRUE.
END IF
IF (za == 1 .AND. na == 2) THEN
kia = k2sh
h2scase = .TRUE.
END IF
hij = 0.5_dp*pia(na)*pib(nb)
IF (h2scase) THEN
fij = 0.5_dp*(kia+kjb)*hij
ELSE
IF ((la == 0 .AND. lb == 1) .OR. (la == 1 .AND. lb == 0)) THEN
fij = ksp*hij*kab*fen
ELSE
fij = 0.5_dp*(kia+kjb)*hij*kab*fen
END IF
END IF
IF (iatom <= jatom) THEN
fhua = fhua+2.0_dp*fij*sblock(i, j)*pblock(i, j)*kcnl(la)
fhub = fhub+2.0_dp*fij*sblock(i, j)*pblock(i, j)*kcnl(lb)
ELSE
fhua = fhua+2.0_dp*fij*sblock(j, i)*pblock(j, i)*kcnl(la)
fhub = fhub+2.0_dp*fij*sblock(j, i)*pblock(j, i)*kcnl(lb)
END IF
END IF
END DO
END DO
IF (calculate_forces) THEN
!qxtb
! Derivative wrt coordination number is missing here
!qxtb
IF (diagblock) THEN
fhua = 0.5_dp*fhua
fhub = 0.5_dp*fhub
END IF
fhua = -1.0_dp*fhua
fhub = -1.0_dp*fhub
! iatom
atom_a = atom_of_kind(iatom)
DO i = 1, dcnum(iatom)%neighbors
katom = dcnum(iatom)%nlist(i)
kkind = kind_of(katom)
rik = dcnum(iatom)%rik(:, i)
drk = SQRT(SUM(rik(:)**2))
fdik(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
atom_c = atom_of_kind(katom)
force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a)-fdik(:)
force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c)+fdik(:)
IF (use_virial) THEN
CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
END IF
IF (atprop%stress) THEN
CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fdik, rik)
CALL virial_pair_force(atprop%atstress(:, :, katom), -0.5_dp, fdik, rik)
END IF
END DO
! iatom
atom_b = atom_of_kind(jatom)
DO i = 1, dcnum(jatom)%neighbors
katom = dcnum(jatom)%nlist(i)
kkind = kind_of(katom)
rik = dcnum(jatom)%rik(:, i)
drk = SQRT(SUM(rik(:)**2))
fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
atom_c = atom_of_kind(katom)
force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b)-fdik(:)
force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c)+fdik(:)
IF (use_virial) THEN
CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
END IF
IF (atprop%stress) THEN
CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fdik, rik)
CALL virial_pair_force(atprop%atstress(:, :, katom), -0.5_dp, fdik, rik)
END IF
END DO
IF (diagblock) THEN
force_ab = 0._dp
ELSE
! force from R dendent Huckel element
n1 = SIZE(fblock, 1)
n2 = SIZE(fblock, 2)
f0 = 1.0_dp
IF (irow == iatom) f0 = -1.0_dp
ALLOCATE (dfblock(n1, n2))
dfblock = 0.0_dp
DO j = 1, natorb_b
lb = laob(j)
nb = naob(j)
DO i = 1, natorb_a
la = laoa(i)
na = naoa(i)
dfij = 0.5_dp*(kl(la)+kl(lb))*0.5_dp*(huckel(na, iatom)+huckel(nb, jatom))*dpia(na)*pib(nb)+ &
0.5_dp*(kl(la)+kl(lb))*0.5_dp*(huckel(na, iatom)+huckel(nb, jatom))*pia(na)*dpib(nb)
IF (iatom <= jatom) THEN
dfblock(i, j) = fblock(i, j)+dfij*sblock(i, j)*kab*fen
kia = kl(la)
kjb = kl(lb)
h2scase = .FALSE.
IF (zb == 1 .AND. nb == 2) THEN
kjb = k2sh
h2scase = .TRUE.
END IF
IF (za == 1 .AND. na == 2) THEN
kia = k2sh
h2scase = .TRUE.
END IF
dhij = 0.5_dp*(huckel(na, iatom)+huckel(nb, jatom))*(dpia(na)*pib(nb)+pia(na)*dpib(nb))
IF (h2scase) THEN
dfij = 0.5_dp*(kia+kjb)*dhij
ELSE
dfblock(j, i) = fblock(j, i)+dfij*sblock(j, i)*kab*fen
IF ((la == 0 .AND. lb == 1) .OR. (la == 1 .AND. lb == 0)) THEN
dfij = ksp*dhij*kab*fen
ELSE
dfij = 0.5_dp*(kia+kjb)*dhij*kab*fen
END IF
END IF
IF (iatom <= jatom) THEN
dfblock(i, j) = dfblock(i, j)+dfij*sblock(i, j)
ELSE
dfblock(j, i) = dfblock(j, i)+dfij*sblock(j, i)
END IF
END DO
END DO
dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
DO ir = 1, 3
foab = 2.0_dp*dfp*rij(ir)/dr
! force from overlap matrix contribution to H
DO j = 1, natorb_b
lb = laob(j)
nb = naob(j)
DO i = 1, natorb_a
la = laoa(i)
na = naoa(i)
fij = 0.5_dp*(kl(la)+kl(lb))*0.5_dp*(huckel(na, iatom)+huckel(nb, jatom))*pia(na)*pib(nb)
foab = foab+2.0_dp*f0*dfblock(j, i)*pblock(j, i)*rij(ir)/dr
foab = foab+2.0_dp*f0*kab*fen*fij*dsblocks(i)%block(j, i)*pblock(j, i)
kia = kl(la)
kjb = kl(lb)
h2scase = .FALSE.
IF (zb == 1 .AND. nb == 2) THEN
kjb = k2sh
h2scase = .TRUE.
END IF
IF (za == 1 .AND. na == 2) THEN
kia = k2sh
h2scase = .TRUE.
END IF
hij = 0.5_dp*(huckel(na, iatom)+huckel(nb, jatom))*pia(na)*pib(nb)
IF (h2scase) THEN
fij = 0.5_dp*(kia+kjb)*hij
ELSE
IF ((la == 0 .AND. lb == 1) .OR. (la == 1 .AND. lb == 0)) THEN
fij = ksp*hij*kab*fen
ELSE
fij = 0.5_dp*(kia+kjb)*hij*kab*fen
END IF
END IF
IF (iatom <= jatom) THEN
foab = foab+2.0_dp*fij*dsblocks(ir+1)%block(i, j)*pblock(i, j)
ELSE
foab = foab+2.0_dp*fij*dsblocks(ir+1)%block(j, i)*pblock(j, i)
END IF
END DO
END DO
force_ab(ir) = foab
END DO
IF (use_virial) THEN
CALL virial_pair_force(virial%pv_virial, -f0, force_ab, rij)
IF (atprop%stress) THEN
f1 = 0.5_dp*f0
CALL virial_pair_force(atprop%atstress(:, :, iatom), -f1, force_ab, rij)
CALL virial_pair_force(atprop%atstress(:, :, jatom), -f1, force_ab, rij)
END IF
END IF
DEALLOCATE (dfblock)
END IF
IF (use_virial) THEN
CALL virial_pair_force(virial%pv_virial, -f0, force_ab, rij)
IF (atprop%stress) THEN
f1 = 0.5_dp*f0
CALL virial_pair_force(atprop%atstress(:, :, iatom), -f1, force_ab, rij)
CALL virial_pair_force(atprop%atstress(:, :, jatom), -f1, force_ab, rij)
END IF
END IF
END IF
IF (calculate_forces) THEN
@ -496,12 +654,20 @@ CONTAINS
CALL get_xtb_atom_param(xtb_atom_a, rcut=rcuta)
CALL get_xtb_atom_param(xtb_atom_b, rcut=rcutb)
rcut = rcuta+rcutb
CALL gamma_rab_sr(gblock, dr, nsa, kappaa, etaa, nsb, kappab, etab, kg, rcut)
IF (irow == iatom) THEN
CALL gamma_rab_sr(gblock, dr, lmaxa+1, kappaa, etaa, lmaxb+1, kappab, etab, kg, rcut)
ELSE
CALL gamma_rab_sr(gblock, dr, lmaxb+1, kappab, etab, lmaxa+1, kappaa, etaa, kg, rcut)
END IF
IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
n1 = SIZE(gblock, 1)
n2 = SIZE(gblock, 2)
ALLOCATE (dgamma(n1, n2))
CALL dgamma_rab_sr(dgamma, dr, nsa, kappaa, etaa, nsb, kappab, etab, kg, rcut)
IF (irow == iatom) THEN
CALL dgamma_rab_sr(dgamma, dr, lmaxa+1, kappaa, etaa, lmaxb+1, kappab, etab, kg, rcut)
ELSE
CALL dgamma_rab_sr(dgamma, dr, lmaxb+1, kappab, etab, lmaxa+1, kappaa, etaa, kg, rcut)
END IF
DO i = 1, 3
CPASSERT(ASSOCIATED(dgblocks(i+1)%block))
IF (irow == iatom) THEN
@ -540,9 +706,6 @@ CONTAINS
END IF
END IF
END IF
!qxtb
! calculate halogen correction term here
!qxtb
END IF
END DO
@ -552,17 +715,15 @@ CONTAINS
CALL dbcsr_finalize(gamma_matrix(i)%matrix)
ENDDO
CALL set_ks_env(ks_env, gamma_matrix=gamma_matrix)
DO i = 1, SIZE(matrix_s, 1)
DO img = 1, nimg
CALL dbcsr_finalize(matrix_s(i, img)%matrix)
END DO
ENDDO
DO i = 1, SIZE(matrix_h, 1)
DO img = 1, nimg
CALL dbcsr_finalize(matrix_h(i, img)%matrix)
END DO
ENDDO
!qxtb
! calculate halogen correction term here
!qxtb
! set repulsive energy
CALL mp_sum(erep, para_env%group)
energy%repulsive = erep
@ -574,6 +735,8 @@ CONTAINS
DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
END DO
DEALLOCATE (dcnum)
ELSE
DEALLOCATE (dcnum)
END IF
! deallocate Huckel parameters
@ -619,6 +782,7 @@ CONTAINS
"DFT%PRINT%AO_MATRICES/OVERLAP")
END IF
DEALLOCATE (kind_of)
IF (calculate_forces) THEN
IF (SIZE(matrix_p, 1) == 2) THEN
DO img = 1, nimg
@ -649,7 +813,7 @@ CONTAINS
INTEGER :: atom_a, handle, iatom, ikind, img, &
iounit, is, ispin, natom, natorb, &
nkind, ns, nspins
INTEGER, DIMENSION(25) :: nao
INTEGER, DIMENSION(25) :: lao
INTEGER, DIMENSION(5) :: occ
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: aocharge, mcharge
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: charges
@ -711,14 +875,14 @@ CONTAINS
DO ikind = 1, nkind
CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
CALL get_xtb_atom_param(xtb_kind, natorb=natorb, nao=nao, occupation=occ)
CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
ALLOCATE (aocharge(natorb))
DO iatom = 1, natom
atom_a = atomic_kind_set(ikind)%atom_list(iatom)
CALL ao_charges(matrix_p, matrix_s, aocharge, atom_a, para_env)
charges(atom_a, :) = REAL(occ(:), KIND=dp)
DO is = 1, natorb
ns = nao(is)
ns = lao(is)+1
charges(atom_a, ns) = charges(atom_a, ns)-aocharge(is)
END DO
mcharge(atom_a) = SUM(charges(atom_a, :))
@ -778,10 +942,10 @@ CONTAINS
!>
!> \param gmat ...
!> \param rab ...
!> \param nsa ...
!> \param nla ...
!> \param kappaa ...
!> \param etaa ...
!> \param nsb ...
!> \param nlb ...
!> \param kappab ...
!> \param etab ...
!> \param kg ...
@ -790,13 +954,13 @@ CONTAINS
!> 10.2018 JGH
!> \version 1.1
! **************************************************************************************************
SUBROUTINE gamma_rab_sr(gmat, rab, nsa, kappaa, etaa, nsb, kappab, etab, kg, rcut)
SUBROUTINE gamma_rab_sr(gmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: gmat
REAL(dp), INTENT(IN) :: rab
INTEGER, INTENT(IN) :: nsa
INTEGER, INTENT(IN) :: nla
REAL(dp), DIMENSION(:), INTENT(IN) :: kappaa
REAL(dp), INTENT(IN) :: etaa
INTEGER, INTENT(IN) :: nsb
INTEGER, INTENT(IN) :: nlb
REAL(dp), DIMENSION(:), INTENT(IN) :: kappab
REAL(dp), INTENT(IN) :: etab, kg, rcut
@ -807,10 +971,10 @@ CONTAINS
REAL(KIND=dp) :: fcut, r, rk, x
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eta
ALLOCATE (eta(nsa, nsb))
ALLOCATE (eta(nla, nlb))
DO j = 1, nsb
DO i = 1, nsa
DO j = 1, nlb
DO i = 1, nla
eta(i, j) = 1._dp/(etaa*(1._dp+kappaa(i)))+1._dp/(etab*(1._dp+kappab(j)))
eta(i, j) = 2._dp/eta(i, j)
END DO
@ -823,6 +987,7 @@ CONTAINS
! do nothing
ELSE
rk = rab**kg
eta = eta**(-kg)
IF (rab < rcut-rsmooth) THEN
fcut = 1.0_dp
ELSE
@ -830,7 +995,8 @@ CONTAINS
x = r/rsmooth
fcut = -10._dp*x**5+15._dp*x**4-10._dp*x**3+1._dp
END IF
gmat(:, :) = gmat(:, :)+fcut*((1._dp/(rk+eta(:, :)**(-kg)))**(1._dp/kg)-1.0_dp/rab)
gmat(:, :) = gmat(:, :)+fcut*(1._dp/(rk+eta(:, :)))**(1._dp/kg)
gmat(:, :) = gmat(:, :)-fcut*1.0_dp/rab
END IF
DEALLOCATE (eta)
@ -846,10 +1012,10 @@ CONTAINS
!>
!> \param dgmat ...
!> \param rab ...
!> \param nsa ...
!> \param nla ...
!> \param kappaa ...
!> \param etaa ...
!> \param nsb ...
!> \param nlb ...
!> \param kappab ...
!> \param etab ...
!> \param kg ...
@ -858,13 +1024,13 @@ CONTAINS
!> 10.2018 JGH
!> \version 1.1
! **************************************************************************************************
SUBROUTINE dgamma_rab_sr(dgmat, rab, nsa, kappaa, etaa, nsb, kappab, etab, kg, rcut)
SUBROUTINE dgamma_rab_sr(dgmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: dgmat
REAL(dp), INTENT(IN) :: rab
INTEGER, INTENT(IN) :: nsa
INTEGER, INTENT(IN) :: nla
REAL(dp), DIMENSION(:), INTENT(IN) :: kappaa
REAL(dp), INTENT(IN) :: etaa
INTEGER, INTENT(IN) :: nsb
INTEGER, INTENT(IN) :: nlb
REAL(dp), DIMENSION(:), INTENT(IN) :: kappab
REAL(dp), INTENT(IN) :: etab, kg, rcut
@ -875,10 +1041,10 @@ CONTAINS
REAL(KIND=dp) :: dfcut, fcut, r, rk, x
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eta
ALLOCATE (eta(nsa, nsb))
ALLOCATE (eta(nla, nlb))
DO j = 1, nsb
DO i = 1, nsa
DO j = 1, nlb
DO i = 1, nla
eta(i, j) = 1._dp/(etaa*(1._dp+kappaa(i)))+1._dp/(etab*(1._dp+kappab(j)))
eta(i, j) = 2._dp/eta(i, j)
END DO
@ -888,8 +1054,9 @@ CONTAINS
! on site terms
dgmat(:, :) = 0.0_dp
ELSEIF (rab > rcut) THEN
! do nothing
dgmat(:, :) = 0.0_dp
ELSE
eta = eta**(-kg)
rk = rab**kg
IF (rab < rcut-rsmooth) THEN
fcut = 1.0_dp
@ -899,107 +1066,16 @@ CONTAINS
x = r/rsmooth
fcut = -10._dp*x**5+15._dp*x**4-10._dp*x**3+1._dp
dfcut = -50._dp*x**4+60._dp*x**3-30._dp*x**2
dfcut = dfcut/rsmooth
END IF
dgmat(:, :) = dgmat(:, :)+dfcut*((1._dp/(rk+eta(:, :)**(-kg)))**(1._dp/kg)-1.0_dp/rab)
dgmat(:, :) = dgmat(:, :)-fcut*((1._dp/(rk+eta(:, :)**(-kg)))**(1._dp/kg-1.0_dp)* &
(1._dp/(rk+eta(:, :)**(-kg)))**2*rab/rk- &
1.0_dp/rab**2)
dgmat(:, :) = dfcut*(1._dp/(rk+eta(:, :)))**(1._dp/kg)
dgmat(:, :) = dgmat(:, :)-dfcut/rab+fcut/rab**2
dgmat(:, :) = dgmat(:, :)-fcut/(rk+eta(:, :))*(1._dp/(rk+eta(:, :)))**(1._dp/kg)*rk/rab
END IF
DEALLOCATE (eta)
END SUBROUTINE dgamma_rab_sr
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param nderivative ...
!> \param nimg ...
!> \param matrices ...
!> \param mnames ...
!> \param sab_nl ...
! **************************************************************************************************
SUBROUTINE setup_matrices2(qs_env, nderivative, nimg, matrices, mnames, sab_nl)
TYPE(qs_environment_type), POINTER :: qs_env
INTEGER, INTENT(IN) :: nderivative, nimg
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrices
CHARACTER(LEN=*) :: mnames
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: sab_nl
CHARACTER(len=*), PARAMETER :: routineN = 'setup_matrices2', &
routineP = moduleN//':'//routineN
CHARACTER(1) :: symmetry_type
CHARACTER(LEN=default_string_length) :: matnames
INTEGER :: i, img, natom, neighbor_list_id, nkind, &
nmat, nsgf
INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
NULLIFY (particle_set, atomic_kind_set)
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set, &
dbcsr_dist=dbcsr_dist, &
neighbor_list_id=neighbor_list_id)
nkind = SIZE(atomic_kind_set)
natom = SIZE(particle_set)
CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
ALLOCATE (first_sgf(natom))
ALLOCATE (last_sgf(natom))
CALL get_particle_set(particle_set, qs_kind_set, &
first_sgf=first_sgf, &
last_sgf=last_sgf)
nmat = 0
IF (nderivative == 0) nmat = 1
IF (nderivative == 1) nmat = 4
IF (nderivative == 2) nmat = 10
CPASSERT(nmat > 0)
ALLOCATE (row_blk_sizes(natom))
CALL convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
CALL dbcsr_allocate_matrix_set(matrices, nmat, nimg)
! Up to 2nd derivative take care to get the symmetries correct
DO img = 1, nimg
DO i = 1, nmat
IF (i .GT. 1) THEN
matnames = TRIM(mnames)//" DERIVATIVE MATRIX xTB"
symmetry_type = dbcsr_type_antisymmetric
IF (i .GT. 4) symmetry_type = dbcsr_type_symmetric
ELSE
symmetry_type = dbcsr_type_symmetric
matnames = TRIM(mnames)//" MATRIX xTB"
END IF
ALLOCATE (matrices(i, img)%matrix)
CALL dbcsr_create(matrix=matrices(i, img)%matrix, &
name=TRIM(matnames), &
dist=dbcsr_dist, matrix_type=symmetry_type, &
row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
nze=0, mutable_work=.TRUE.)
CALL cp_dbcsr_alloc_block_from_nbl(matrices(i, img)%matrix, sab_nl)
END DO
END DO
DEALLOCATE (first_sgf)
DEALLOCATE (last_sgf)
DEALLOCATE (row_blk_sizes)
END SUBROUTINE setup_matrices2
! **************************************************************************************************
!> \brief ...
@ -1044,8 +1120,8 @@ CONTAINS
ikind = kind_of(iatom)
qs_kind => qs_kind_set(ikind)
CALL get_qs_kind(qs_kind, xtb_parameter=xtb_parameter)
CALL get_xtb_atom_param(xtb_parameter, nshell=ns)
row_blk_sizes(iatom) = ns
CALL get_xtb_atom_param(xtb_parameter, lmax=ns)
row_blk_sizes(iatom) = ns+1
ENDDO
CALL dbcsr_allocate_matrix_set(gammat, nmat)

30
src/xtb_parameters.F
View file

@ -151,9 +151,6 @@ CONTAINS
found = .TRUE.
CALL parser_get_object(parser, param%eta)
CALL parser_get_object(parser, param%xgamma)
!xtb
param%xgamma = param%xgamma/evolt
!xtb
CALL parser_get_object(parser, param%alpha)
CALL parser_get_object(parser, param%zneff)
DO i = 1, 5
@ -257,16 +254,16 @@ CONTAINS
IF (ia >= 49 .AND. ia <= 57) THEN
nshell = nshell+1
k = nshell
param%nval(i) = ia-48
param%nval(k) = ia-48
SELECT CASE (label (2:2))
CASE ("s", "S")
param%lval(i) = 0
param%lval(k) = 0
CASE ("p", "P")
param%lval(i) = 1
param%lval(k) = 1
CASE ("d", "D")
param%lval(i) = 2
param%lval(k) = 2
CASE ("f", "F")
param%lval(i) = 3
param%lval(k) = 3
CASE DEFAULT
CPABORT("xTB PARAMETER ERROR")
END SELECT
@ -359,6 +356,7 @@ CONTAINS
routineP = moduleN//':'//routineN
INTEGER :: i, is, l, na
REAL(KIND=dp), DIMENSION(5) :: kp
IF (param%defined) THEN
! AO to shell pointer
@ -384,6 +382,22 @@ CONTAINS
END IF
! orbital energies [evolt] -> [a.u.]
param%hen = param%hen/evolt
! some forgotten scaling parameters (not in orig. paper)
param%xgamma = 0.1_dp*param%xgamma
param%kpoly(:) = 0.01*param%kpoly(:)
! we have 1/6 g * q**3 (not 1/3)
param%xgamma = 2.0_dp*param%xgamma
! kappa is read for each shell, but depends on l, reindex here
kp(:) = param%kappa(:)
param%kappa = 0._dp
DO is = 1, param%nshell
l = param%lval(is)
IF (param%kappa(l+1) /= 0._dp) THEN
! test that we don't have ambigous kappa values
CPASSERT(ABS(param%kappa(l+1)-kp(is)) < 1.e-12_dp)
ENDIF
param%kappa(l+1) = kp(is)
END DO
END IF
END SUBROUTINE xtb_parameters_set

25
tests/xTB/regtest-1/ch2o.inp
View file

@ -1,8 +1,16 @@
&FORCE_EVAL
&DFT
&PRINT
&AO_MATRICES
# OVERLAP
# CORE_HAMILTONIAN
# KOHN_SHAM_MATRIX
&END
&END
&QS
METHOD xTB
&xTB
DO_EWALD T
&PARAMETER
DISPERSION_PARAMETER_FILE dftd3.dat
&END PARAMETER
@ -18,7 +26,12 @@
&END QS
&SCF
SCF_GUESS MOPAC
MAX_SCF 20
MAX_SCF 50
EPS_SCF 1.e-8
&MIXING
METHOD DIRECT_P_MIXING
ALPHA 0.4
&END
&END SCF
&POISSON
&EWALD
@ -42,6 +55,14 @@
&END FORCE_EVAL
&GLOBAL
PROJECT CH2O
RUN_TYPE ENERGY
# RUN_TYPE ENERGY_FORCE
RUN_TYPE DEBUG
PRINT_LEVEL LOW
&END GLOBAL
&DEBUG
DEBUG_FORCES T
DEBUG_STRESS_TENSOR F
DX 0.0001
EPS_NO_ERROR_CHECK 0.000001
STOP_ON_MISMATCH F
&END DEBUG