Skip to content

Commit

Permalink
Bug fixes, consistent parameters, disabled CN term (forces not debugged)
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Mar 19, 2019
1 parent 2c6ccf5 commit 841ab0f
Show file tree
Hide file tree
Showing 10 changed files with 423 additions and 262 deletions.
10 changes: 5 additions & 5 deletions data/xTB_parameters
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 change: 1 addition & 0 deletions src/common/mathlib.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -394,6 +394,7 @@ SUBROUTINE diamat_all(a, eigval, dac)

! 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
Expand Down
9 changes: 1 addition & 8 deletions src/mulliken.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -709,7 +709,6 @@ SUBROUTINE ao_charges_1(p_matrix, s_matrix, charges, iatom, para_env)
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

Expand All @@ -730,18 +729,12 @@ SUBROUTINE ao_charges_1(p_matrix, s_matrix, charges, iatom, para_env)
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)
Expand Down
1 change: 1 addition & 0 deletions src/qs_diis.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -441,6 +441,7 @@ SUBROUTINE qs_diis_b_step(diis_buffer, mo_array, kc, sc, delta, error_max, &

! 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)
Expand Down
4 changes: 4 additions & 0 deletions src/qs_dispersion_pairpot.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1069,6 +1069,8 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
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, &
Expand Down Expand Up @@ -1657,6 +1659,8 @@ SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculat
DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
END DO
DEALLOCATE (dcnum)
ELSE
DEALLOCATE (dcnum)
END IF
END IF

Expand Down
12 changes: 10 additions & 2 deletions src/qs_environment.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -660,9 +660,8 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
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
Expand Down Expand Up @@ -813,6 +812,15 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell

! *** 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)
Expand Down
95 changes: 69 additions & 26 deletions src/xtb_coulomb.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,10 +127,9 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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
Expand Down Expand Up @@ -162,13 +161,14 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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, &
Expand Down Expand Up @@ -208,6 +208,10 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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)
Expand All @@ -217,13 +221,15 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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))
Expand Down Expand Up @@ -285,12 +291,14 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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
Expand All @@ -300,7 +308,6 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
! 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)
Expand All @@ -322,7 +329,8 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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
Expand All @@ -331,7 +339,8 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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)+ &
Expand All @@ -347,10 +356,18 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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))
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) = SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, i+1))+ &
gmcharge(iatom, i+1)*mcharge(iatom)
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)
Expand Down Expand Up @@ -398,12 +415,27 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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))
Expand All @@ -412,15 +444,19 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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
ikind = kind_of(iatom)
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)
Expand All @@ -430,10 +466,17 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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) = 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) = fi
fij(i) = fij(i)+fi
END DO
IF (use_virial) THEN
CALL virial_pair_force(virial%pv_virial, 1._dp, fij, rij)
Expand Down Expand Up @@ -463,8 +506,8 @@ SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, &
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)
Expand Down

0 comments on commit 841ab0f

Please sign in to comment.