Repository: Atomistica/atomistica (language: python)

Issue:
Missing stress calculation with DFTB
I'm using the ASE interface with the `native.TightBinding` potential, and querying stresses from the `ase.Atoms` object leads to the following `RuntimeError`:
```
RuntimeError: Traceback (most recent call last)
  File "build/potentials_factory_f90.f90", line 223
    *wpot_per_at* argument present but not supported by potential TightBinding.
```
It looks to me that the `dense_notb_energy_and_forces`, which is ultimately called by the calculator, supports the `wpot_per_at` argument (although it does not seem to do anything with it), so perhaps the check on whether `wpot_per_at` is bound is extraneous?

Here's an MWE to reproduce the error: [dftb_stresses.zip](https://github.com/Atomistica/atomistica/files/8341251/dftb_stresses.zip). Thanks for taking the time.




Discussion / hints:
Here's the patch I applied to get rid of the error by generating python bindings with per atom quantities, but it causes a segfault. [patch.txt](https://github.com/Atomistica/atomistica/files/8769021/patch.txt)

Fix patch:
diff --git a/ChangeLog.md b/ChangeLog.md
index fd93e75d..b534cf2a 100644
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -1,6 +1,11 @@
 Change log
 ==========
 
+v1.1.0 (23Jun22)
+----------------
+
+- ENH: Per atom virial for NOTB
+
 v1.0.4 (24Apr22)
 ----------------
 
diff --git a/requirements.txt b/requirements.txt
index b0fe3134..7725f9ff 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,1 +1,2 @@
 ase>=3.15.0
+pytest
diff --git a/src/notb/dense/dense_forces.f90 b/src/notb/dense/dense_forces.f90
index 9f467493..0c9d43ba 100644
--- a/src/notb/dense/dense_forces.f90
+++ b/src/notb/dense/dense_forces.f90
@@ -50,21 +50,23 @@ module dense_force_calc
   !**********************************************************************
   ! Force calculation including k-space summation
   !**********************************************************************
-  subroutine forces(p, nl, tb, db, f, wpot, wpot_per_bond, error)
+  subroutine forces(p, nl, tb, db, f, wpot, wpot_per_at, wpot_per_bond, error)
     implicit none
 
-    type(particles_t), intent(inout)   :: p
-    type(neighbors_t), intent(in)      :: nl
-    type(dense_hamiltonian_t), intent(in)  :: tb
-    type(materials_t), intent(in)      :: db
-    real(DP), intent(inout)            :: f(3, p%maxnatloc)
-    real(DP), intent(inout)            :: wpot(3, 3)
+    type(particles_t),         intent(inout) :: p
+    type(neighbors_t),         intent(in)    :: nl
+    type(dense_hamiltonian_t), intent(in)    :: tb
+    type(materials_t),         intent(in)    :: db
+    real(DP),                  intent(inout) :: f(3, p%maxnatloc)
+    real(DP),                  intent(inout) :: wpot(3, 3)
 #ifdef LAMMPS
-    real(DP), intent(inout), optional  :: wpot_per_bond(6, nl%neighbors_size)
+    real(DP),        optional, intent(inout) :: wpot_per_at(6, p%maxnatloc)
+    real(DP),        optional, intent(inout) :: wpot_per_bond(6, nl%neighbors_size)
 #else
-    real(DP), intent(inout), optional  :: wpot_per_bond(3, 3, nl%neighbors_size)
+    real(DP),        optional, intent(inout) :: wpot_per_at(3, 3, p%maxnatloc)
+    real(DP),        optional, intent(inout) :: wpot_per_bond(3, 3, nl%neighbors_size)
 #endif
-    integer, intent(inout), optional   :: error
+    integer,         optional, intent(inout) :: error
 
     ! ---
 
@@ -108,9 +110,9 @@ subroutine forces(p, nl, tb, db, f, wpot, wpot_per_bond, error)
     list(0:7,  8) = [7,  1,  2,  3,  4,  5,  6,  7             ] ! pd
     list(0:10, 9) = [10, 1,  2,  3,  4,  5,  6,  7,  8,  9,  10] ! spd
 
-    lo     = (/0,1,1,1,2,2,2,2,2/)
+    lo = (/0,1,1,1,2,2,2,2,2/)
 
-    w      = 0.0_DP
+    w = 0.0_DP
 
     error_loc  = ERROR_NONE
 
@@ -120,7 +122,7 @@ subroutine forces(p, nl, tb, db, f, wpot, wpot_per_bond, error)
     !$omp& private(I, Imu, mu, J, Jnu, nu, kk, l, m, ni) &
     !$omp& private(noI, noJ, nr, q, rho_ImuJnu, rIJ) &
     !$omp& private(Sdiff_ij, Sdiff_ji, wij) &
-    !$omp& shared(db, f, list, lo, nl, p, tb, wpot_per_bond) &
+    !$omp& shared(db, f, list, lo, nl, p, tb, wpot_per_at, wpot_per_bond) &
     !$omp& shared(tb_at, tb_rho, tb_e) &
     !$omp& reduction(+:w) reduction(+:error_loc)
 
@@ -262,6 +264,11 @@ subroutine forces(p, nl, tb, db, f, wpot, wpot_per_bond, error)
                    if (present(wpot_per_bond)) then
                       SUM_VIRIAL(wpot_per_bond, ni, wij)
                    endif
+                   if (present(wpot_per_at)) then
+                      wij = wij/2
+                      SUM_VIRIAL(wpot_per_at, I, wij)
+                      SUM_VIRIAL(wpot_per_at, J, wij)
+                   endif
 
                 endif I_lt_J
 
diff --git a/src/notb/dense/dense_hs.f90 b/src/notb/dense/dense_hs.f90
index f80e650e..61075645 100644
--- a/src/notb/dense/dense_hs.f90
+++ b/src/notb/dense/dense_hs.f90
@@ -505,9 +505,9 @@ subroutine calc_matrix_elements(this, db, el_i, el_j, my_r_ij, &
     m  = max(noi, noj)
     nr = list(0, m)                    ! number of non-vanishing matrix elements
     a  = interval(db%HS(eli, elj), my_r_ij, error)  ! determine the interval dr can be found in (within xs)
-    PASS_ERROR(error)
+    PASS_ERROR_WITH_INFO("eli = " // eli // ", elj = " // elj, error)
     b  = interval(db%HS(elj, eli), my_r_ij, error)  ! determine the interval dr can be found in (within xs)
-    PASS_ERROR(error)
+    PASS_ERROR_WITH_INFO("elj = " // elj // ", eli = " // eli, error)
     do q = 1, nr
        bo = list(q, m)   ! bond order
        if (bo <= 0) then
diff --git a/src/notb/dense/dense_notb.f90 b/src/notb/dense/dense_notb.f90
index 0dee036c..19500720 100644
--- a/src/notb/dense/dense_notb.f90
+++ b/src/notb/dense/dense_notb.f90
@@ -23,6 +23,7 @@
 !   shared:directory
 !   dependencies:../materials.f90,dense_hamiltonian_type.f90,dense_hamiltonian.f90,c_dense_hamiltonian.cpp,dense_forces.f90,dense_hs.f90,dense_repulsion.f90,solver/dense_occupation.f90,solver/dense_solver_lapack.f90,solver/dense_solver_cp.f90,solver/dense_solver_dispatch.f90,dense_scc.f90,analysis/dense_bonds.f90
 !   classtype:dense_notb_t classname:TightBinding interface:potentials
+!   features:per_at,per_bond
 ! @endmeta
 
 #include "macros.inc"
@@ -627,7 +628,7 @@ subroutine dense_notb_energy_and_forces(this, p, nl, epot, f, wpot, q, &
     !   e_bs     = Band-structure energy
     !   e_rep    = Repulsive energy
     !   e_atomic = Energy of the individual charge neutral atoms
-    call forces(p, nl, this%tb, this%mat, f, wpot, wpot_per_bond, error=ierror)
+    call forces(p, nl, this%tb, this%mat, f, wpot, wpot_per_at, wpot_per_bond, error=ierror)
     PASS_ERROR_AND_STOP_TIMER("dense_notb_energy_and_forces", ierror)
     epot  = epot + e_bs(this%solver, this%tb) + &
          e_rep(this%tb, this%mat, p, nl) - e_atomic(this%tb, p) 
diff --git a/src/support/nonuniform_spline.f90 b/src/support/nonuniform_spline.f90
index 76d6d63a..e4c878a3 100644
--- a/src/support/nonuniform_spline.f90
+++ b/src/support/nonuniform_spline.f90
@@ -459,6 +459,10 @@ integer function nonuniform_spline_interval(this, x, ierror)
 
     ! ---
 
+    if (this%n < 0) then
+      RAISE_ERROR("Spline not initialized.", ierror)
+    endif
+
     !
     ! if x out of interpolation region
     !


Test patch:
diff --git a/.github/workflows/run_tests.yml b/.github/workflows/run_tests.yml
index b6d3a921..89bdfc23 100644
--- a/.github/workflows/run_tests.yml
+++ b/.github/workflows/run_tests.yml
@@ -59,7 +59,7 @@ jobs:
     - name: Run Python binding tests
       run: |
         cd tests
-        MIO='../mio-1-1' DFTB3='../3ob-3-1' python run_tests.py
+        MIO='../mio-1-1' DFTB3='../3ob-3-1' pytest -s --verbose
 
     - name: Run standalone examples
       run: |
diff --git a/tests/bulk_properties.py b/tests/test_bulk_properties.py
similarity index 90%
rename from tests/bulk_properties.py
rename to tests/test_bulk_properties.py
index 0531d082..3e2628f5 100755
--- a/tests/bulk_properties.py
+++ b/tests/test_bulk_properties.py
@@ -40,7 +40,7 @@
 from ase.units import GPa
 
 from atomistica import *
-from atomistica.tests import test_cubic_elastic_constants
+from atomistica.tests import test_cubic_elastic_constants as cubic_elastic_constants
 
 from ase.lattice.cubic import Diamond, BodyCenteredCubic
 from ase.lattice.cubic import FaceCenteredCubic, SimpleCubic
@@ -183,8 +183,8 @@ class TestElasticConstants(unittest.TestCase):
 
     def test_elastic_constants(self):
         for pot, par, mats in tests:
-             test_cubic_elastic_constants(mats, pot, par, sx, dev_thres,
-                                         test=self)
+             cubic_elastic_constants(mats, pot, par, sx, dev_thres,
+                                     test=self)
 
 ###
 
@@ -192,24 +192,24 @@ def test_elastic_constants(self):
     nok = 0
     nfail = 0
     for pot, par, mats in tests:
-        if len(sys.argv) > 1:
-            found = False
-            if par is not None:
-                for keyword in sys.argv[1:]:
-                    if '__ref__' in par:
-                        if par['__ref__'].lower().find(keyword.lower()) != -1:
-                            found = True
-            try:
-                potname = pot.__name__
-            except:
-                potname = pot.__class__.__name__
-            for keyword in sys.argv[1:]:
-                if potname.lower().find(keyword.lower()) != -1:
-                    found = True
-            if not found:
-                continue
-
-        _nok, _nfail = test_cubic_elastic_constants(mats, pot, par, sx, dev_thres)
+        #if len(sys.argv) > 1:
+        #    found = False
+        #    if par is not None:
+        #        for keyword in sys.argv[1:]:
+        #            if '__ref__' in par:
+        #                if par['__ref__'].lower().find(keyword.lower()) != -1:
+        #                    found = True
+        #    try:
+        #        potname = pot.__name__
+        #    except:
+        #        potname = pot.__class__.__name__
+        #    for keyword in sys.argv[1:]:
+        #        if potname.lower().find(keyword.lower()) != -1:
+        #            found = True
+        #    if not found:
+        #        continue
+
+        _nok, _nfail = cubic_elastic_constants(mats, pot, par, sx, dev_thres)
         nok += _nok
         nfail += _nfail
     print('{0} tests passed, {1} tests failed.'.format(nok, nfail))
diff --git a/tests/coulomb.py b/tests/test_coulomb.py
similarity index 100%
rename from tests/coulomb.py
rename to tests/test_coulomb.py
diff --git a/tests/dftb3.py b/tests/test_dftb3.py
similarity index 95%
rename from tests/dftb3.py
rename to tests/test_dftb3.py
index 83debc18..b06ed715 100644
--- a/tests/dftb3.py
+++ b/tests/test_dftb3.py
@@ -25,8 +25,6 @@
 Test the DFTB3 parametrizations.
 """
 
-from __future__ import print_function
-
 import os
 import sys
 import unittest
@@ -97,7 +95,7 @@
 
 ###
 
-def run_dftb3_test(test=None, tol=0.05):
+def test_dftb3(test=None, tol=0.05):
     mio_database_folder = os.getenv('MIO')
     if mio_database_folder is None:
         raise RuntimeError('Please use environment variable MIO to specify path to mio Slater-Koster tables.')
@@ -281,20 +279,3 @@ def run_dftb3_test(test=None, tol=0.05):
             test.assertTrue(success_DFTB3)
             test.assertTrue(success_DFTB3_XH)
             test.assertTrue(success_DFTB3_3ob_XH)
-
-###
-
-class TestMIO(unittest.TestCase):
-
-   def test_dftb3(self):
-
-       if os.getenv('MIO') is None:
-           print('Skipping DFTB3 test. Specify path to mio Slater-Koster ' \
-                 'tables in MIO environment variable if you want to run it.')
-       else:
-           run_dftb3_test(self)
-
-###
-
-if __name__ == '__main__':
-    run_dftb3_test()
diff --git a/tests/dimers.py b/tests/test_dimers.py
similarity index 100%
rename from tests/dimers.py
rename to tests/test_dimers.py
diff --git a/tests/eam_special_cases.py b/tests/test_eam_special_cases.py
similarity index 93%
rename from tests/eam_special_cases.py
rename to tests/test_eam_special_cases.py
index 6ff8eb6d..eea471fd 100755
--- a/tests/eam_special_cases.py
+++ b/tests/test_eam_special_cases.py
@@ -36,7 +36,7 @@
 import ase.io as io
 
 from atomistica import TabulatedAlloyEAM
-from atomistica.tests import test_forces
+from atomistica.tests import test_forces as forces
 
 ###
 
@@ -59,7 +59,7 @@ def test_dense_forces(self):
             a = orig_a.copy()
             a.set_cell(fac*a.cell, scale_atoms=True)
             a.calc = c
-            ffd, f0, maxdf = test_forces(a, dx=dx)
+            ffd, f0, maxdf = forces(a, dx=dx)
             if maxdf > tol:
                 nfail += 1
                 print("forces .failed.")
diff --git a/tests/forces_and_virial.py b/tests/test_forces_and_virial.py
similarity index 91%
rename from tests/forces_and_virial.py
rename to tests/test_forces_and_virial.py
index 7d181b2a..b86f911a 100755
--- a/tests/forces_and_virial.py
+++ b/tests/test_forces_and_virial.py
@@ -18,10 +18,6 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 # ======================================================================
-#! /usr/bin/env python
-
-
-from __future__ import print_function
 
 import math
 import sys
@@ -40,9 +36,9 @@
 
 import atomistica.native as native
 from atomistica import *
-from atomistica.tests import test_forces, test_potential, test_virial
-
-# import ase_ext_io as io
+from atomistica.tests import test_forces as forces
+from atomistica.tests import test_potential as potential
+from atomistica.tests import test_virial as virial
 
 ###
 
@@ -206,26 +202,26 @@ def assign_charges(a, els):
 
 ###
 
-def run_forces_and_virial_test(test=None):
+def test_forces_and_virial(test=None):
     nok = 0
     nfail = 0
     for pot, par, mats in tests:
-        if len(sys.argv) > 1:
-            found = False
-            if par is not None:
-                for keyword in sys.argv[1:]:
-                    if '__ref__' in par:
-                        if par['__ref__'].lower().find(keyword.lower()) != -1:
-                            found = True
-            try:
-                potname = pot.__name__
-            except:
-                potname = pot.__class__.__name__
-            for keyword in sys.argv[1:]:
-                if potname.lower().find(keyword.lower()) != -1:
-                    found = True
-            if not found:
-                continue
+        #if len(sys.argv) > 1:
+        #    found = False
+        #    if par is not None:
+        #        for keyword in sys.argv[1:]:
+        #            if '__ref__' in par:
+        #                if par['__ref__'].lower().find(keyword.lower()) != -1:
+        #                    found = True
+        #    try:
+        #        potname = pot.__name__
+        #    except:
+        #        potname = pot.__class__.__name__
+        #    for keyword in sys.argv[1:]:
+        #        if potname.lower().find(keyword.lower()) != -1:
+        #            found = True
+        #    if not found:
+        #        continue
 
         try:
             potname = pot.__name__
@@ -279,7 +275,7 @@ def run_forces_and_virial_test(test=None):
                         print('--- using random mask ---')
                     c.set_mask(mask)
 
-                    ffd, f0, maxdf = test_forces(a, dx=dx)
+                    ffd, f0, maxdf = forces(a, dx=dx)
         
                     if test is None:
                         if abs(maxdf) < tol:
@@ -305,7 +301,7 @@ def run_forces_and_virial_test(test=None):
                       test.assertTrue(abs(maxdf) < tol,
                                         msg=errmsg+'; forces')
 
-                    sfd, s0, maxds = test_virial(a, de=dx)
+                    sfd, s0, maxds = virial(a, de=dx)
 
                     if test is None:
                         if abs(maxds) < tol:
@@ -328,7 +324,7 @@ def run_forces_and_virial_test(test=None):
                         test.assertTrue(abs(maxds) < tol,
                                         msg=errmsg+'; virial')
 
-                    pfd, p0, maxdp = test_potential(a, dq=dx)
+                    pfd, p0, maxdp = potential(a, dq=dx)
 
                     if test is None:
                         if abs(maxdp) < tol:
@@ -355,14 +351,3 @@ def run_forces_and_virial_test(test=None):
     if test is None:
         print('{0} tests passed, {1} tests failed.'.format(nok, nfail))
 
-###
-
-class TestForcesAndVirial(unittest.TestCase):
-
-    def test_forces_and_virial(self):
-        run_forces_and_virial_test(self)
-
-###
-
-if __name__ == '__main__':
-    run_forces_and_virial_test()
diff --git a/tests/test_io.py b/tests/test_io.py
index 2480ea82..e37c5fc0 100755
--- a/tests/test_io.py
+++ b/tests/test_io.py
@@ -10,31 +10,23 @@
 from atomistica import Rebo2, Rebo2Scr
 
 
-class IOTest(unittest.TestCase):
-    def test_energy(self):
-        vac = 8
-        dist_min = 1.2
-
-        atoms = Atoms('CC', positions=[[0, 0, 0], [dist_min, 0, 0]])
-        atoms.center(vacuum=vac)
-
-        for calc in [Rebo2(), Rebo2Scr()]:
-            atoms.calc = calc
-            energy = atoms.get_potential_energy()
-            forces_ac = atoms.get_forces()
-            stress = atoms.get_stress()
-
-            fname = 'structure.traj'
-            atoms.write(fname)
-            atoms = io.read(fname)
-
-            self.assertTrue(np.abs(energy - atoms.get_potential_energy())
-                            < 1e-10)
-            self.assertTrue((np.abs(forces_ac - atoms.get_forces())
-                             < 1e-10).all())
-            self.assertTrue((np.abs(stress - atoms.get_stress())
-                             < 1e-10).all())
-
-
-if __name__ == '__main__':
-    unittest.main()
+def test_energy():
+    vac = 8
+    dist_min = 1.2
+
+    atoms = Atoms('CC', positions=[[0, 0, 0], [dist_min, 0, 0]])
+    atoms.center(vacuum=vac)
+
+    for calc in [Rebo2(), Rebo2Scr()]:
+        atoms.calc = calc
+        energy = atoms.get_potential_energy()
+        forces_ac = atoms.get_forces()
+        stress = atoms.get_stress()
+
+        fname = 'structure.traj'
+        atoms.write(fname)
+        atoms = io.read(fname)
+
+        assert np.abs(energy - atoms.get_potential_energy()) < 1e-10
+        assert (np.abs(forces_ac - atoms.get_forces()) < 1e-10).all()
+        assert (np.abs(stress - atoms.get_stress()) < 1e-10).all()
diff --git a/tests/mask.py b/tests/test_mask.py
similarity index 100%
rename from tests/mask.py
rename to tests/test_mask.py
diff --git a/tests/mio.py b/tests/test_mio.py
similarity index 100%
rename from tests/mio.py
rename to tests/test_mio.py
diff --git a/tests/neighbor_list.py b/tests/test_neighbor_list.py
similarity index 100%
rename from tests/neighbor_list.py
rename to tests/test_neighbor_list.py
diff --git a/tests/pbc.py b/tests/test_pbc.py
similarity index 100%
rename from tests/pbc.py
rename to tests/test_pbc.py
diff --git a/tests/rebo2_molecules.py b/tests/test_rebo2_molecules.py
similarity index 94%
rename from tests/rebo2_molecules.py
rename to tests/test_rebo2_molecules.py
index a20c2ac0..c0f3793d 100755
--- a/tests/rebo2_molecules.py
+++ b/tests/test_rebo2_molecules.py
@@ -18,7 +18,6 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 # ======================================================================
-#! /usr/bin/env python
 
 
 """Check energies against the databases published in
@@ -111,7 +110,7 @@ def molecule(mol):
 
 ###
 
-def run_rebo2_molecules_test(test=None):
+def test_rebo2_molecules(test=None):
     for potname, c, reference_database in [ 
         ( 'Rebo2', atomistica.Rebo2(),
           Brenner_et_al_CH ),
@@ -126,9 +125,9 @@ def run_rebo2_molecules_test(test=None):
         nfailed = 0
     
         for mol, edft, de in reference_database:
-            if len(sys.argv) > 1:
-                if mol not in sys.argv[1:]:
-                    continue
+            #if len(sys.argv) > 1:
+            #    if mol not in sys.argv[1:]:
+            #        continue
 
             eref = edft-de
 
@@ -160,15 +159,3 @@ def run_rebo2_molecules_test(test=None):
         if test is None:
             print('{0} molecule tests ok, {1} molecule tests failed.' \
                 .format(nok, nfailed))
-
-###
-
-class TestREBO2Molecules(unittest.TestCase):
-
-    def test_rebo2_molecules(self):
-        run_rebo2_molecules_test(self)
-
-###
-
-if __name__ == '__main__':
-    run_rebo2_molecules_test()
diff --git a/tests/surface_properties.py b/tests/test_surface_properties.py
similarity index 92%
rename from tests/surface_properties.py
rename to tests/test_surface_properties.py
index 2371c0eb..994ae2ab 100755
--- a/tests/surface_properties.py
+++ b/tests/test_surface_properties.py
@@ -41,7 +41,7 @@
     from ase.atoms import string2symbols
 
 from atomistica import *
-from atomistica.tests import test_surface_energies
+from atomistica.tests import test_surface_energies as surface_energies
 
 from ase.lattice.cubic import Diamond, FaceCenteredCubic
 from ase.lattice.compounds import B1, B2, B3, L1_2
@@ -294,31 +294,6 @@ def dia_100_2x1(sym, a0):
 
 ###
 
-class TestSurfaceEnergies(unittest.TestCase):
-
-    def test_surface_energies(self):
-        for pot, par, mats in tests:
-            test_surface_energies(mats, pot, par, nx, vacuum, test=self)
-
-###
-
-if __name__ == '__main__':
+def test_surface_energies():
     for pot, par, mats in tests:
-        if len(sys.argv) > 1:
-            found = False
-            if par is not None:
-                for keyword in sys.argv[1:]:
-                    if '__ref__' in par:
-                        if par['__ref__'].lower().find(keyword.lower()) != -1:
-                            found = True
-            try:
-                potname = pot.__name__
-            except:
-                potname = pot.__class__.__name__
-            for keyword in sys.argv[1:]:
-                if potname.lower().find(keyword.lower()) != -1:
-                    found = True
-            if not found:
-                continue
-
-        test_surface_energies(mats, pot, par, nx, vacuum, dump=True)
+        surface_energies(mats, pot, par, nx, vacuum)
diff --git a/tests/run_tests.py b/tests/test_tb_stresses.py
old mode 100755
new mode 100644
similarity index 52%
rename from tests/run_tests.py
rename to tests/test_tb_stresses.py
index b9229bc1..28f10d69
--- a/tests/run_tests.py
+++ b/tests/test_tb_stresses.py
@@ -18,26 +18,46 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 # ======================================================================
-#! /usr/bin/env python
 
+"""
+Test computation of stresses from tight binding.
+"""
 
+import os
+import sys
 import unittest
 
-from coulomb import *
-from eam_special_cases import *
-from neighbor_list import *
-from pbc import *
-from bulk_properties import *
-from surface_properties import *
-from forces_and_virial import *
-from rebo2_molecules import *
-from mask import *
-from mio import *
-from dftb3 import *
-from dimers import *
-from test_io import IOTest
-
-###
-
-unittest.main()
+import numpy as np
+
+import ase.io as io
+from ase.build import molecule
+from ase.optimize import FIRE
+
+import atomistica.native as native
+from atomistica import Atomistica
+from atomistica.tests import test_virial as virial
+
+def test_stresses(test=None, de=1e-6):
+    database_folder = os.getenv('MIO')
+    if database_folder is None:
+        raise RuntimeError('Please use environment variable MIO to specify path to mio Slater-Koster tables.')
+
+    calc = Atomistica(
+        [ native.TightBinding(
+        database_folder = database_folder,
+        SolverLAPACK = dict(electronic_T=0.001),
+        ) ],
+        avgn = 1000
+        )
+
+    a = io.read('aC_small.cfg')
+    a.calc = calc
+    s = a.get_stress() * a.get_volume()
+    s_at = a.get_stresses()
+
+    np.testing.assert_array_almost_equal(s, s_at.sum(axis=0))
+
+    sfd, s0, maxds = virial(a, de=de)
+    np.testing.assert_array_almost_equal(s0, sfd)
+