From 841ab0f2eb6c628e0f7a24548cdbd60b895dcf67 Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Mon, 17 Dec 2018 14:28:45 +0100 Subject: [PATCH] Bug fixes, consistent parameters, disabled CN term (forces not debugged) --- data/xTB_parameters | 10 +- src/common/mathlib.F | 1 + src/mulliken.F | 9 +- src/qs_diis.F | 1 + src/qs_dispersion_pairpot.F | 4 + src/qs_environment.F | 12 +- src/xtb_coulomb.F | 95 +++++-- src/xtb_matrices.F | 500 ++++++++++++++++++++--------------- src/xtb_parameters.F | 30 ++- tests/xTB/regtest-1/ch2o.inp | 25 +- 10 files changed, 424 insertions(+), 263 deletions(-) diff --git a/data/xTB_parameters b/data/xTB_parameters index ae16d85de4..66687910eb 100644 --- a/data/xTB_parameters +++ b/data/xTB_parameters @@ -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 diff --git a/src/common/mathlib.F b/src/common/mathlib.F index 8e6748a273..9063cb6675 100644 --- a/src/common/mathlib.F +++ b/src/common/mathlib.F @@ -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 diff --git a/src/mulliken.F b/src/mulliken.F index 38ae4004ec..1e01288d1a 100644 --- a/src/mulliken.F +++ b/src/mulliken.F @@ -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) diff --git a/src/qs_diis.F b/src/qs_diis.F index 8967fae334..6832eaa514 100644 --- a/src/qs_diis.F +++ b/src/qs_diis.F @@ -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) diff --git a/src/qs_dispersion_pairpot.F b/src/qs_dispersion_pairpot.F index 6bd8c8e466..36ab1d6b7d 100644 --- a/src/qs_dispersion_pairpot.F +++ b/src/qs_dispersion_pairpot.F @@ -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 diff --git a/src/qs_environment.F b/src/qs_environment.F index c220156a3e..815cb73db3 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -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) diff --git a/src/xtb_coulomb.F b/src/xtb_coulomb.F index 607f8cbfd0..6fa9f099ff 100644 --- a/src/xtb_coulomb.F +++ b/src/xtb_coulomb.F @@ -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) diff --git a/src/xtb_matrices.F b/src/xtb_matrices.F index 158c9f0b07..f7b24aa80d 100644 --- a/src/xtb_matrices.F +++ b/src/xtb_matrices.F @@ -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 = 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) diff --git a/src/xtb_parameters.F b/src/xtb_parameters.F index 0556172165..2b8aaad897 100644 --- a/src/xtb_parameters.F +++ b/src/xtb_parameters.F @@ -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 diff --git a/tests/xTB/regtest-1/ch2o.inp b/tests/xTB/regtest-1/ch2o.inp index b3dccf2c1b..52d138de00 100644 --- a/tests/xTB/regtest-1/ch2o.inp +++ b/tests/xTB/regtest-1/ch2o.inp @@ -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