diff --git a/src/qs_dftb3_methods.F b/src/qs_dftb3_methods.F index 130315e627..e0643cdb11 100644 --- a/src/qs_dftb3_methods.F +++ b/src/qs_dftb3_methods.F @@ -225,10 +225,12 @@ CONTAINS fij(i) = -gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i)) END IF END DO - CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij) + fi = 1.0_dp + IF (iatom == jatom) fi = 0.5_dp + CALL virial_pair_force(virial%pv_virial, fi, fij, rij) IF (atprop%stress) THEN - CALL virial_pair_force(atprop%atstress(:, :, irow), -0.5_dp, fij, rij) - CALL virial_pair_force(atprop%atstress(:, :, icol), -0.5_dp, fij, rij) + CALL virial_pair_force(atprop%atstress(:, :, irow), fi*0.5_dp, fij, rij) + CALL virial_pair_force(atprop%atstress(:, :, icol), fi*0.5_dp, fij, rij) END IF END IF END DO @@ -281,10 +283,12 @@ CONTAINS fij(i) = fi END DO IF (use_virial) THEN - CALL virial_pair_force(virial%pv_virial, 1._dp, fij, rij) + fi = 1.0_dp + IF (iatom == jatom) fi = 0.5_dp + CALL virial_pair_force(virial%pv_virial, fi, fij, rij) IF (atprop%stress) THEN - CALL virial_pair_force(atprop%atstress(:, :, iatom), 0.5_dp, fij, rij) - CALL virial_pair_force(atprop%atstress(:, :, jatom), 0.5_dp, fij, rij) + CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij) + CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij) END IF END IF diff --git a/src/qs_dftb_coulomb.F b/src/qs_dftb_coulomb.F index c0c5d247ab..83bb2df4c5 100644 --- a/src/qs_dftb_coulomb.F +++ b/src/qs_dftb_coulomb.F @@ -231,10 +231,12 @@ CONTAINS END IF END DO IF (use_virial) THEN - CALL virial_pair_force(virial%pv_virial, 1._dp, fij, rij) + fi = 1.0_dp + IF (iatom == jatom) fi = 0.5_dp + CALL virial_pair_force(virial%pv_virial, fi, fij, rij) IF (atprop%stress) THEN - CALL virial_pair_force(atprop%atstress(:, :, iatom), 0.5_dp, fij, rij) - CALL virial_pair_force(atprop%atstress(:, :, jatom), 0.5_dp, fij, rij) + CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij) + CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij) END IF END IF END IF @@ -433,10 +435,12 @@ CONTAINS fij(i) = -2.0_dp*gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i)) END IF END DO - CALL virial_pair_force(virial%pv_virial, -1._dp, fij, rij) + fi = 1.0_dp + IF (iatom == jatom) fi = 0.5_dp + CALL virial_pair_force(virial%pv_virial, fi, fij, rij) IF (atprop%stress) THEN - CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, fij, rij) - CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, fij, rij) + CALL virial_pair_force(atprop%atstress(:, :, iatom), 0.5_dp*fi, fij, rij) + CALL virial_pair_force(atprop%atstress(:, :, jatom), 0.5_dp*fi, fij, rij) END IF END IF END DO @@ -492,10 +496,12 @@ CONTAINS fij(i) = fi END DO IF (use_virial) THEN - CALL virial_pair_force(virial%pv_virial, 1._dp, fij, rij) + fi = 1.0_dp + IF (iatom == jatom) fi = 0.5_dp + CALL virial_pair_force(virial%pv_virial, fi, fij, rij) IF (atprop%stress) THEN - CALL virial_pair_force(atprop%atstress(:, :, iatom), 0.5_dp, fij, rij) - CALL virial_pair_force(atprop%atstress(:, :, jatom), 0.5_dp, fij, rij) + CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij) + CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij) END IF END IF END IF @@ -737,7 +743,7 @@ CONTAINS dsblock = dsblock-dsblockm dsblock = dsblock/(2.0_dp*ddr) - clist%acint(1:n1, 1:n2, i) = dsblock(1:n1, 1:n2) + clist%acint(1:n1, 1:n2, i) = -dsblock(1:n1, 1:n2) ENDDO DEALLOCATE (dsblock, dsblockm) END IF diff --git a/src/qs_dftb_matrices.F b/src/qs_dftb_matrices.F index 9ded2e77a4..0ab94a1283 100644 --- a/src/qs_dftb_matrices.F +++ b/src/qs_dftb_matrices.F @@ -368,6 +368,9 @@ CONTAINS END IF ENDDO IF (use_virial) THEN +!deb iatom=jatom f0=0.5*f0 + IF (iatom == jatom) f0 = 0.5_dp*f0 +!deb CALL virial_pair_force(virial%pv_virial, -f0, force_ab, rij) CALL virial_pair_force(virial%pv_virial, -f0, force_w, rij) IF (atprop%stress) THEN @@ -413,10 +416,12 @@ CONTAINS force(jkind)%repulsive(:, atom_b) = & force(jkind)%repulsive(:, atom_b)+force_rr(:) IF (use_virial) THEN - CALL virial_pair_force(virial%pv_virial, -1._dp, force_rr, rij) + f0 = -1.0_dp + IF (iatom == jatom) f0 = -0.5_dp + CALL virial_pair_force(virial%pv_virial, f0, force_rr, rij) IF (atprop%stress) THEN - CALL virial_pair_force(atprop%atstress(:, :, iatom), -0.5_dp, force_rr, rij) - CALL virial_pair_force(atprop%atstress(:, :, jatom), -0.5_dp, force_rr, rij) + CALL virial_pair_force(atprop%atstress(:, :, iatom), f0*0.5_dp, force_rr, rij) + CALL virial_pair_force(atprop%atstress(:, :, jatom), f0*0.5_dp, force_rr, rij) END IF END IF END IF diff --git a/src/xtb_coulomb.F b/src/xtb_coulomb.F index 42e542b7a9..cb200d8a08 100644 --- a/src/xtb_coulomb.F +++ b/src/xtb_coulomb.F @@ -486,10 +486,10 @@ CONTAINS DO j = 1, nj la = laoa(i)+1 lb = laob(j)+1 - gcij(i, j) = 0.5_dp*(gchrg(irow, la, 1)+gchrg(icol, lb, 1)) + gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1)+gchrg(jatom, lb, 1)) END DO END DO - gmij = 0.5_dp*(gmcharge(irow, 1)+gmcharge(icol, 1)) + gmij = 0.5_dp*(gmcharge(iatom, 1)+gmcharge(jatom, 1)) icol = MAX(iatom, jatom) irow = MIN(iatom, jatom) NULLIFY (pblock) @@ -506,7 +506,7 @@ CONTAINS f1 = -2.0_dp*SUM(TRANSPOSE(pblock)*dsint(:, :, i)*gcij) f2 = -2.0_dp*gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i)) END IF - fij(i) = fij(i)+f1+f2 + fij(i) = f1+f2 END DO DEALLOCATE (gcij) fi = 1.0_dp @@ -576,7 +576,7 @@ CONTAINS IF (calculate_forces) THEN atom_i = atom_of_kind(iatom) atom_j = atom_of_kind(jatom) - IF (irow == jatom) THEN + IF (irow /= iatom) THEN gmij = -gmij gcij = -gcij END IF diff --git a/tests/DFTB/regtest-scc-2/TEST_FILES b/tests/DFTB/regtest-scc-2/TEST_FILES index 028d0abc50..862e8dccb8 100644 --- a/tests/DFTB/regtest-scc-2/TEST_FILES +++ b/tests/DFTB/regtest-scc-2/TEST_FILES @@ -4,7 +4,7 @@ base1.inp 1 1.0E-12 -4.08611632954554 base2.inp 1 1.0E-12 -4.08611630192434 base3.inp 1 1.0E-12 -4.08519715577319 -str1.inp 1 1.0E-12 -32.67545045484233 +str1.inp 1 1.0E-12 -4.08612900850098 str2.inp 1 1.0E-12 -32.67544573650541 str3.inp 1 1.0E-12 -4.08521550626912 #EOF diff --git a/tests/DFTB/regtest-scc-2/str1.inp b/tests/DFTB/regtest-scc-2/str1.inp index c2a0a8082d..e0d6b8c7d2 100644 --- a/tests/DFTB/regtest-scc-2/str1.inp +++ b/tests/DFTB/regtest-scc-2/str1.inp @@ -1,4 +1,4 @@ -@SET NREP 2 +@SET NREP 1 &FORCE_EVAL STRESS_TENSOR ANALYTICAL &DFT diff --git a/tests/DFTB/regtest-scc/TEST_FILES b/tests/DFTB/regtest-scc/TEST_FILES index 931f0be4c0..8ee9fa158e 100644 --- a/tests/DFTB/regtest-scc/TEST_FILES +++ b/tests/DFTB/regtest-scc/TEST_FILES @@ -47,11 +47,11 @@ h2o_disp2.inp 1 4e-09 - h2o_disp3.inp 1 3e-09 -130.76542980403056 # DFTB3 h2o-6.inp 1 4e-09 -131.25337238175840 -h2o-7.inp 1 1e-08 -131.26536731088510 +h2o-7.inp 1 1e-08 -131.26536755467674 # Properties stress h2o-atprop1.inp 1 4e-09 -130.77763673429416 -h2o-atprop2.inp 1 8e-09 -131.68099551186123 -h2o-atprop3.inp 1 8e-09 -131.33465952358338 +h2o-atprop2.inp 1 8e-09 -131.68099550861825 +h2o-atprop3.inp 1 8e-09 -131.33465894123017 # c_kp1.inp 1 1.0E-14 -13.76936383448730 c_kp2.inp 1 7e-14 -13.76786043591555 diff --git a/tests/xTB/regtest-1/TEST_FILES b/tests/xTB/regtest-1/TEST_FILES index 5d7640790c..efacf0f749 100644 --- a/tests/xTB/regtest-1/TEST_FILES +++ b/tests/xTB/regtest-1/TEST_FILES @@ -10,11 +10,13 @@ ch2o_smear.inp 1 1.0E-12 -7.19650182 tmol.inp 1 1.0E-12 -41.90853885813245 h2.inp 1 1.0E-12 -1.03458111733093 h2o-md.inp 1 1.0E-12 -185.15731412756392 -h2o_str.inp 1 1.0E-12 -5.76530887622206 -h2o-atprop.inp 1 4.0E-08 -185.17649920422642 +h2o_str.inp 1 1.0E-12 -46.12964016870387 +h2o-atprop.inp 1 1.0E-10 -185.17668389962344 h2o-atprop0.inp 1 1.0E-12 -187.45499307089042 si_geo.inp 1 1.0E-12 -14.55678348854957 si_kp.inp 1 1.0E-12 -14.75431646498735 h2o_dimer.inp 1 1.0E-12 -11.54506384130837 AdeThyvdW.inp 1 1.0E-12 -58.89501875211175 +ice.inp 1 1.0E-10 -370.501223771625689 +ice2.inp 1 1.0E-10 -46.30633640965694 #EOF diff --git a/tests/xTB/regtest-1/h2o_str.inp b/tests/xTB/regtest-1/h2o_str.inp index 199cff7a1d..2487bb7db9 100644 --- a/tests/xTB/regtest-1/h2o_str.inp +++ b/tests/xTB/regtest-1/h2o_str.inp @@ -1,3 +1,4 @@ +@SET NREP 2 &FORCE_EVAL STRESS_TENSOR ANALYTICAL &DFT @@ -5,38 +6,34 @@ METHOD xTB &xTB DO_EWALD T - COULOMB_INTERACTION T - TB3_INTERACTION T - &PARAMETER - DISPERSION_PARAMETER_FILE dftd3.dat - &END PARAMETER &END + EPS_DEFAULT 1.E-14 &END QS &SCF SCF_GUESS MOPAC &MIXING METHOD DIRECT_P_MIXING - ALPHA 0.10 + ALPHA 0.25 &END - MAX_SCF 2000 + MAX_SCF 200 + EPS_SCF 1.e-10 &END SCF &KPOINTS - SCHEME NONE - #SCHEME GAMMA - #SCHEME MONKHORST-PACK 1 1 1 + SCHEME NONE &END KPOINTS &POISSON &EWALD - ALPHA 2.0 - EWALD_TYPE SPME - GMAX 50 - O_SPLINE 6 + ALPHA 3.0 &END EWALD &END POISSON &END DFT &SUBSYS + &TOPOLOGY + MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP} + &END &CELL ABC 7.0 7.0 7.0 + MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP} &END CELL &COORD O -4.583 5.333 1.560 @@ -48,14 +45,13 @@ &GLOBAL PROJECT h2o RUN_TYPE DEBUG -# RUN_TYPE ENERGY_FORCE PRINT_LEVEL LOW &END GLOBAL &DEBUG DEBUG_FORCES F DEBUG_STRESS_TENSOR T - DX 0.001 - EPS_NO_ERROR_CHECK 0.000001 - STOP_ON_MISMATCH F + DX 0.0002 + EPS_NO_ERROR_CHECK 0.0001 + STOP_ON_MISMATCH T &END DEBUG diff --git a/tests/xTB/regtest-1/hf.inp b/tests/xTB/regtest-1/hf.inp index 0074e7e15d..54453707fb 100644 --- a/tests/xTB/regtest-1/hf.inp +++ b/tests/xTB/regtest-1/hf.inp @@ -62,5 +62,5 @@ DEBUG_STRESS_TENSOR T DX 0.002 EPS_NO_ERROR_CHECK 0.000001 - STOP_ON_MISMATCH F + STOP_ON_MISMATCH T &END DEBUG diff --git a/tests/xTB/regtest-1/ice.inp b/tests/xTB/regtest-1/ice.inp new file mode 100644 index 0000000000..3606e09a8d --- /dev/null +++ b/tests/xTB/regtest-1/ice.inp @@ -0,0 +1,85 @@ +@SET NREP 2 +&GLOBAL + PRINT_LEVEL low + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + DEBUG_FORCES no + DEBUG_STRESS_TENSOR yes + EPS_NO_ERROR_CHECK 1.0E-4 + DX 2.E-4 + STOP_ON_MISMATCH yes +&END DEBUG + +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR analytical + &PRINT + &STRESS_TENSOR ON + &END + &END + &DFT + &KPOINTS + SCHEME NONE + &END + &QS + EPS_DEFAULT 1.0E-12 + EXTRAPOLATION ASPC + EXTRAPOLATION_ORDER 3 + METHOD xTB + &xTB + DO_EWALD T + &END xTB + &END QS + &SCF + EPS_SCF 1.0E-6 + MAX_SCF 10 + SCF_GUESS MOPAC + &OT on + MINIMIZER DIIS + PRECONDITIONER FULL_SINGLE_INVERSE + &END OT + &OUTER_SCF on + EPS_SCF 1.0E-6 + MAX_SCF 10 + &END OUTER_SCF + &END SCF + &END DFT + &SUBSYS + &TOPOLOGY + MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP} + &END + &CELL + ABC 6.358 6.358 6.358 + MULTIPLE_UNIT_CELL ${NREP} ${NREP} ${NREP} + &END CELL + &COORD + O 3.4613786876 5.6912049480 0.4543727289 H2O + H 2.8183735096 5.2012524560 -0.0760381983 H2O + H 2.9114406133 6.4013882133 0.9299789968 H2O + O 0.2726839173 5.8131769039 3.6524426114 H2O + H -0.3305134720 5.2542646859 3.0350424519 H2O + H 0.8608072875 6.3736221304 3.0591507235 H2O + O 0.3733209371 2.7092162634 0.4720541405 H2O + H -0.1034698773 2.0668534417 -0.1013800786 H2O + H -0.2883205804 3.2313214887 0.9861387459 H2O + O 5.1679371793 4.2946104941 2.0104883056 H2O + H 4.5349860728 3.7358082901 2.5469828636 H2O + H 4.6439911293 4.8921754874 1.4197090565 H2O + O 1.9804166604 1.0500358045 2.0298628573 H2O + H 2.5408437547 1.7234516156 2.4956106850 H2O + H 1.3203099947 1.5937364881 1.5079432584 H2O + O 1.8865307172 4.1962271939 5.1674455966 H2O + H 1.1978491058 4.7123566273 4.6845538753 H2O + H 1.4552978348 3.6225372508 5.7958754719 H2O + O 3.5550539674 2.8142815783 3.4327698835 H2O + H 2.9552982643 3.3909531502 4.0032656392 H2O + H 3.9895494155 2.1557772081 4.0448159552 H2O + O 5.1227218630 1.1112395839 5.1795427775 H2O + H 5.7513589539 0.5921262315 4.6175072767 H2O + H 4.6567143202 0.4540581626 5.7714857593 H2O + &END COORD + &END SUBSYS +&END FORCE_EVAL + diff --git a/tests/xTB/regtest-1/ice2.inp b/tests/xTB/regtest-1/ice2.inp new file mode 100644 index 0000000000..f600f7067c --- /dev/null +++ b/tests/xTB/regtest-1/ice2.inp @@ -0,0 +1,71 @@ +&GLOBAL + PRINT_LEVEL low + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + DEBUG_FORCES no + DEBUG_STRESS_TENSOR yes + EPS_NO_ERROR_CHECK 1.0E-4 + DX 1.E-3 + STOP_ON_MISMATCH yes +&END DEBUG + +&FORCE_EVAL + METHOD Quickstep + STRESS_TENSOR analytical + &PRINT + &STRESS_TENSOR ON + &END + &END + &DFT + &KPOINTS + SCHEME GAMMA + &END + &QS + EPS_DEFAULT 1.0E-14 + EXTRAPOLATION ASPC + EXTRAPOLATION_ORDER 3 + METHOD xTB + &xTB + DO_EWALD T + &END xTB + &END QS + &SCF + EPS_SCF 1.0E-9 + MAX_SCF 100 + &END SCF + &END DFT + &SUBSYS + &CELL + ABC 6.358 6.358 6.358 + &END CELL + &COORD + O 3.4613786876 5.6912049480 0.4543727289 H2O + H 2.8183735096 5.2012524560 -0.0760381983 H2O + H 2.9114406133 6.4013882133 0.9299789968 H2O + O 0.2726839173 5.8131769039 3.6524426114 H2O + H -0.3305134720 5.2542646859 3.0350424519 H2O + H 0.8608072875 6.3736221304 3.0591507235 H2O + O 0.3733209371 2.7092162634 0.4720541405 H2O + H -0.1034698773 2.0668534417 -0.1013800786 H2O + H -0.2883205804 3.2313214887 0.9861387459 H2O + O 5.1679371793 4.2946104941 2.0104883056 H2O + H 4.5349860728 3.7358082901 2.5469828636 H2O + H 4.6439911293 4.8921754874 1.4197090565 H2O + O 1.9804166604 1.0500358045 2.0298628573 H2O + H 2.5408437547 1.7234516156 2.4956106850 H2O + H 1.3203099947 1.5937364881 1.5079432584 H2O + O 1.8865307172 4.1962271939 5.1674455966 H2O + H 1.1978491058 4.7123566273 4.6845538753 H2O + H 1.4552978348 3.6225372508 5.7958754719 H2O + O 3.5550539674 2.8142815783 3.4327698835 H2O + H 2.9552982643 3.3909531502 4.0032656392 H2O + H 3.9895494155 2.1557772081 4.0448159552 H2O + O 5.1227218630 1.1112395839 5.1795427775 H2O + H 5.7513589539 0.5921262315 4.6175072767 H2O + H 4.6567143202 0.4540581626 5.7714857593 H2O + &END COORD + &END SUBSYS +&END FORCE_EVAL +