band-unfolding / banduppy

Python version ofthe BandUP code
GNU General Public License v3.0
17 stars 7 forks source link

Different number of plane waves ncnt vs nplane #19

Closed dapshenay closed 1 year ago

dapshenay commented 1 year ago

Dear Developers,

I had a sudden problem on one of calculations using VASP. For one of k-points in list there appeared an error:

RuntimeError: *** error - computed ncnt=18134 != input nplane=18133

Could you, please, help to resolve it?

Best regards, Dmitry

Below is the full output:

RuntimeError Traceback (most recent call last) Cell In [46], line 5 3 if not os.path.exists(os.path.join(wd, bands_pickle_file)): 4 print("reading bandstructure") ----> 5 bands = banduppy.BandStructure(code="vasp", spinor=spinorbit, fPOS = "POSCAR", fWAV = "WAVECAR") 6 pickle.dump(bands,open(bands_pickle_file,"wb")) 7 print ("reading bandstructure - done")

File ~/nfs_sm/prog/opt/miniforge3/envs/phono3py/lib/python3.9/site-packages/irrep/bandstructure.py:136, in BandStructure.init(self, fWAV, fWFK, prefix, fPOS, Ecut, IBstart, IBend, kplist, spinor, code, EF, onlysym, spin_channel, refUC, shiftUC, search_cell) 133 if spin_channel=='down' : spin_channel='dw' 135 if code == "vasp": --> 136 self.__init_vasp( 137 fWAV, fPOS, Ecut, IBstart, IBend, kplist, spinor, EF=EF, onlysym=onlysym, refUC=refUC, shiftUC=shiftUC, search_cell=search_cell 138 ) 139 elif code == "abinit": 140 self.__init_abinit( 141 fWFK, Ecut, IBstart, IBend, kplist, EF=EF, onlysym=onlysym, refUC=refUC, shiftUC=shiftUC, search_cell=search_cell 142 )

File ~/nfs_sm/prog/opt/miniforge3/envs/phono3py/lib/python3.9/site-packages/irrep/bandstructure.py:280, in BandStructure.__init_vasp(self, fWAV, fPOS, Ecut, IBstart, IBend, kplist, spinor, EF, onlysym, refUC, shiftUC, search_cell) 278 kplist = np.array([k for k in kplist if k >= 0 and k < NK]) 279 # print (kplist) --> 280 self.kpoints = [ 281 Kpoint( 282 ik, 283 NBin, 284 IBstart, 285 IBend, 286 Ecut, 287 Ecut0, 288 self.RecLattice, 289 symmetries_SG=self.spacegroup.symmetries, 290 spinor=self.spinor, 291 WCF=WCF, 292 ) 293 for ik in kplist 294 ]

File ~/nfs_sm/prog/opt/miniforge3/envs/phono3py/lib/python3.9/site-packages/irrep/bandstructure.py:281, in (.0) 278 kplist = np.array([k for k in kplist if k >= 0 and k < NK]) 279 # print (kplist) 280 self.kpoints = [ --> 281 Kpoint( 282 ik, 283 NBin, 284 IBstart, 285 IBend, 286 Ecut, 287 Ecut0, 288 self.RecLattice, 289 symmetries_SG=self.spacegroup.symmetries, 290 spinor=self.spinor, 291 WCF=WCF, 292 ) 293 for ik in kplist 294 ]

File ~/nfs_sm/prog/opt/miniforge3/envs/phono3py/lib/python3.9/site-packages/irrep/kpoint.py:197, in Kpoint.init(self, ik, NBin, IBstart, IBend, Ecut, Ecut0, RecLattice, symmetriesSG, spinor, code, kpt, npw, fWFK, WCF, prefix, kptxml, flag, usepaw, eigenval, spin_channel, IBstartE) 194 self.symmetries_SG = symmetries_SG # lazy_property needs it 196 if code.lower() == "vasp": --> 197 self.WF, self.ig = self.__init_vasp( 198 WCF, ik, NBin, IBstart, IBend, Ecut, Ecut0 199 ) 200 elif code.lower() == "abinit": 201 self.WF, self.ig = self.__init_abinit( 202 fWFK, 203 ik, (...) 212 usepaw=usepaw, 213 )

File ~/nfs_sm/prog/opt/miniforge3/envs/phono3py/lib/python3.9/site-packages/irrep/kpoint.py:607, in Kpoint.__init_vasp(self, WCF, ik, NBin, IBstart, IBend, Ecut, Ecut0) 604 except BaseException: 605 self.upper = np.NaN --> 607 ig = calc_gvectors( 608 self.K, self.RecLattice, Ecut0, npw, Ecut, spinor=self.spinor 609 ) 610 selectG = np.hstack((ig[3], ig[3] + int(npw / 2))) if self.spinor else ig[3] 611 WF = np.array( 612 [ 613 WCF.record(3 + ik * (NBin + 1) + ib, npw, np.complex64)[selectG] 614 for ib in range(IBstart, IBend) 615 ] 616 )

File ~/nfs_sm/prog/opt/miniforge3/envs/phono3py/lib/python3.9/site-packages/irrep/gvectors.py:146, in calc_gvectors(K, RecLattice, Ecut, nplane, Ecut1, thresh, spinor, nplanemax) 144 else: 145 if ncnt != nplane: --> 146 raise RuntimeError( 147 "*** error - computed ncnt={0} != input nplane={1}".format( 148 ncnt, nplane 149 ) 150 ) 151 igall = np.array(igall, dtype=int) 152 ng = igall.max(axis=0) - igall.min(axis=0)

RuntimeError: *** error - computed ncnt=18134 != input nplane=18133

stepan-tsirkin commented 1 year ago

Hi, this is an issue with the IrRep package.

More precisely, it is because vasp does not write to the WAVECAR which plane waves it uses. Therefore, IrRep while reading the bandstructure tries to reproduce the behavior of VASP in choosing and ordering the plane waves, which fall inside the cut-off sphere. In very rear cases, due to numerical errors, some plane wave may be considered within the sphwere by one code, but not by the other, which happened in your case with one vector.

as a workaround, we introduced a correction coefficient here: https://github.com/stepan-tsirkin/irrep/blob/522460369668f0b36cd02e40cc6647e6c7ae0cf2/irrep/gvectors.py#L32-L38

this correction needs to be really small (of order 1e-6 or so), either positive or negative. unfortunately, it is hard-coded, so you will need to download irrep, modify this line and install from your local repository, which is very inconvenient, because it is needed only for a particular calculation.

@miraola, can we make this correction available as a parameter of bandStructure for the user to modify (instead of hard-coding)

dapshenay commented 1 year ago

Thank you for answer! correction = -1e-7 worked for me.

I'd like also to ask, if it is possible to use manually set Ecut parameter in this case. Then, in the case of error user could reduce Ecut a little bit to discard some plane waves near the boundary of cut-off sphere. In the present case, as far as I understand, Ecut can be reduced after igall was generated. Is it possible to exchange the sequence?

stepan-tsirkin commented 1 year ago

sure, I opened a PR , you may try and comment on it, if you want.

dapshenay commented 1 year ago

Than you! I updated irrep and now I can use the correction inside the function call. In my case: bands = banduppy.BandStructure(code="vasp", spinor=spinorbit, fPOS = "POSCAR", fWAV = "WAVECAR", _correct_Ecut0=-1e-7)

stepan-tsirkin commented 1 year ago

@dapshenay , thanks for checking, now irrep-1.8.1 (with _correct-Ecut0) is on pip .