vaspwfc.py 文件源码

def gvectors(self, ikpt=1, gamma=False):
        '''
        Generate the G-vectors that satisfies the following relation
            (G + k)**2 / 2 < ENCUT
        '''
        assert 1 <= ikpt  <= self._nkpts,  'Invalid kpoint index!'

        kvec = self._kvecs[ikpt-1]
        # fx, fy, fz = [fftfreq(n) * n for n in self._ngrid]
        # fftfreq in scipy.fftpack is a little different with VASP frequencies
        fx = [ii if ii < self._ngrid[0] / 2 + 1 else ii - self._ngrid[0]
                for ii in range(self._ngrid[0])]
        fy = [jj if jj < self._ngrid[1] / 2 + 1 else jj - self._ngrid[1]
                for jj in range(self._ngrid[1])]
        fz = [kk if kk < self._ngrid[2] / 2 + 1 else kk - self._ngrid[2]
                for kk in range(self._ngrid[2])]
        if gamma:
            # parallel gamma version of VASP WAVECAR exclude some planewave
            # components, -DwNGZHalf
            kgrid = np.array([(fx[ii], fy[jj], fz[kk])
                              for kk in range(self._ngrid[2])
                              for jj in range(self._ngrid[1])
                              for ii in range(self._ngrid[0])
                              if (
                                  (fz[kk] > 0) or
                                  (fz[kk] == 0 and fy[jj] > 0) or
                                  (fz[kk] == 0 and fy[jj] == 0 and fx[ii] >= 0)
                              )], dtype=float)
        else:
            kgrid = np.array([(fx[ii], fy[jj], fz[kk])
                              for kk in range(self._ngrid[2])
                              for jj in range(self._ngrid[1])
                              for ii in range(self._ngrid[0])], dtype=float)

        # Kinetic_Energy = (G + k)**2 / 2
        # HSQDTM    =  hbar**2/(2*ELECTRON MASS)
        KENERGY = HSQDTM * np.linalg.norm(
                    np.dot(kgrid + kvec[np.newaxis,:] , TPI*self._Bcell), axis=1
                )**2
        # find Gvectors where (G + k)**2 / 2 < ENCUT
        Gvec = kgrid[np.where(KENERGY < self._encut)[0]]

        if self._lsoc:
                assert Gvec.shape[0] == self._nplws[ikpt - 1] / 2, \
                       'No. of planewaves not consistent for an SOC WAVECAR! %d %d %d' % \
                       (Gvec.shape[0], self._nplws[ikpt -1], np.prod(self._ngrid))
        else:
            assert Gvec.shape[0] == self._nplws[ikpt - 1], 'No. of planewaves not consistent! %d %d %d' % \
                    (Gvec.shape[0], self._nplws[ikpt -1], np.prod(self._ngrid))

        return np.asarray(Gvec, dtype=int)