From c42694680362a5216b256a4932a01641908feba8 Mon Sep 17 00:00:00 2001 From: "Oleg E. Peil" Date: Thu, 19 Nov 2015 12:11:55 +0100 Subject: [PATCH] Added implementation for NORMION = True The mapping for NORMION = True has been implemented. Also, the orthogonalization loop has been fixed. First of all, orthogonalization should be done separately for each block map 'bl_map'. Second, one has to take into account that the orbital dimensions of the block matrix can vary from block to block. To make that the overlap matrix is non-singular one, thus, has to pass to 'orthogonalize_projector_matrix()' only a view of a submatrix of 'pmat' corresponding to the current block. Two tests to check the simplest cases have been added. --- python/vasp/proj_group.py | 21 +++++++++++++----- .../test/_proj_group/projortho_normion.out | 22 +++++++++++++++++++ python/vasp/test/_proj_group/test_one_site.py | 18 +++++++++++++++ python/vasp/test/_proj_group/test_two_site.py | 21 +++++++++++++++++- 4 files changed, 76 insertions(+), 6 deletions(-) create mode 100644 python/vasp/test/_proj_group/projortho_normion.out diff --git a/python/vasp/proj_group.py b/python/vasp/proj_group.py index 51f2b3ba..488b1ece 100644 --- a/python/vasp/proj_group.py +++ b/python/vasp/proj_group.py @@ -138,8 +138,19 @@ class ProjectorGroup: # Determine the dimension of the projector matrix # and map the blocks to the big matrix if self.normion: -# TODO: add the case of 'normion = True' - raise NotImplementedError("'NORMION = True' is not yet implemented") + block_maps = [] + ndim = 0 + for ish in self.ishells: + _shell = self.shells[ish] + nion, ns, nk, nlm, nb_max = _shell.proj_win.shape + ndim = max(ndim, nlm) + for ion in xrange(nion): + i1_bl = 0 + i2_bl = nlm + block = {'bmat_range': (i1_bl, i2_bl)} + block['shell_ion'] = (ish, ion) + bl_map = [block] + block_maps.append(bl_map) else: block_maps = [] @@ -165,17 +176,17 @@ class ProjectorGroup: for ik in xrange(nk): nb = self.ib_win[ik, isp, 1] - self.ib_win[ik, isp, 0] + 1 # Combine all projectors of the group to one block projector - p_mat[:, :] = 0.0j # !!! Clean-up from the last k-point! for bl_map in block_maps: + p_mat[:, :] = 0.0j # !!! Clean-up from the last k-point and block! for ibl, block in enumerate(bl_map): i1, i2 = block['bmat_range'] ish, ion = block['shell_ion'] shell = self.shells[ish] p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb] # Now orthogonalize the obtained block projector - p_orth, overl, eig = self.orthogonalize_projector_matrix(p_mat) + ibl_max = i2 + p_orth, overl, eig = self.orthogonalize_projector_matrix(p_mat[:ibl_max, :nb]) # Distribute projectors back using the same mapping - for bl_map in block_maps: for ibl, block in enumerate(bl_map): i1, i2 = block['bmat_range'] ish, ion = block['shell_ion'] diff --git a/python/vasp/test/_proj_group/projortho_normion.out b/python/vasp/test/_proj_group/projortho_normion.out new file mode 100644 index 00000000..ff02073f --- /dev/null +++ b/python/vasp/test/_proj_group/projortho_normion.out @@ -0,0 +1,22 @@ +density matrix: [[[[ 0.53073815 0.00000001 -0. -0.00000002 0. ] + [ 0.00000001 0.53074869 0. 0.00000002 -0.00000001] + [-0. 0. 0.54606387 -0. -0.0000001 ] + [-0.00000002 0.00000002 -0. 0.5308456 -0. ] + [ 0. -0.00000001 -0.0000001 -0. 0.54513241]] + + [[ 0.53004807 0.00000002 0. 0.00000001 -0. ] + [ 0.00000002 0.53003083 0. -0.00000001 -0. ] + [ 0. 0. 0.54606372 0. 0.00000009] + [ 0.00000001 -0.00000001 0. 0.52993693 0. ] + [-0. -0. 0.00000009 0. 0.54513254]]]] +overlap matrix: [[[[ 1. -0. -0. -0. -0.] + [-0. 1. 0. 0. -0.] + [-0. 0. 1. 0. -0.] + [-0. 0. 0. 1. 0.] + [-0. -0. -0. 0. 1.]] + + [[ 1. 0. -0. 0. 0.] + [ 0. 1. 0. -0. -0.] + [-0. 0. 1. -0. 0.] + [ 0. -0. -0. 1. -0.] + [ 0. -0. 0. -0. 1.]]]] diff --git a/python/vasp/test/_proj_group/test_one_site.py b/python/vasp/test/_proj_group/test_one_site.py index 8e3416b1..7d1ea2ef 100644 --- a/python/vasp/test/_proj_group/test_one_site.py +++ b/python/vasp/test/_proj_group/test_one_site.py @@ -24,6 +24,7 @@ class TestProjectorGroup(mytest.MyTestCase): Scenarios: - **test** that orthogonalization is correct + - **test** that NORMION = True gives the same result """ def setUp(self): conf_file = _rpath + 'example.cfg' @@ -55,4 +56,21 @@ class TestProjectorGroup(mytest.MyTestCase): expected_file = _rpath + 'projortho.out' self.assertFileEqual(testout, expected_file) +# Scenario 2 + def test_ortho_normion(self): + self.proj_gr.normion = True + self.proj_gr.orthogonalize() + + dens_mat, overl = self.proj_sh.density_matrix(self.el_struct) + + testout = _rpath + 'projortho.out.test' + with open(testout, 'wt') as f: + f.write("density matrix: %s\n"%(dens_mat)) + f.write("overlap matrix: %s\n"%(overl)) + + self.assertEqual(overl, np.eye(5)) + + expected_file = _rpath + 'projortho.out' + self.assertFileEqual(testout, expected_file) + diff --git a/python/vasp/test/_proj_group/test_two_site.py b/python/vasp/test/_proj_group/test_two_site.py index 0c4599d2..de7cca84 100644 --- a/python/vasp/test/_proj_group/test_two_site.py +++ b/python/vasp/test/_proj_group/test_two_site.py @@ -25,7 +25,8 @@ class TestProjectorGroupTwoSite(mytest.MyTestCase): ProjectorGroup(sh_pars, proj_raw) Scenarios: - - **test** that orthogonalization is correct + - **test** that orthogonalization with NORMION = False is correct + - **test** that orthogonalization with NORMION = True is correct """ def setUp(self): conf_file = _rpath + 'example_two_site.cfg' @@ -58,4 +59,22 @@ class TestProjectorGroupTwoSite(mytest.MyTestCase): expected_file = _rpath + 'projortho_2site.out' self.assertFileEqual(testout, expected_file) +# Scenario 2 + def test_ortho_normion(self): + self.proj_gr.normion = True + self.proj_gr.orthogonalize() + + dens_mat, overl = self.proj_sh.density_matrix(self.el_struct) + + testout = _rpath + 'projortho_normion.out.test' + with open(testout, 'wt') as f: + f.write("density matrix: %s\n"%(dens_mat)) + f.write("overlap matrix: %s\n"%(overl)) + + self.assertEqual(overl[0, 0, ...], np.eye(5)) + self.assertEqual(overl[0, 1, ...], np.eye(5)) + + expected_file = _rpath + 'projortho_normion.out' + self.assertFileEqual(testout, expected_file) +