@@ -681,7 +681,7 @@ void Forces::cal_force_nl(ModuleBase::matrix& forcenl)
681681 for (int ik = 0 ;ik < GlobalC::kv.nks ;ik++)
682682 {
683683 if (GlobalV::NSPIN==2 ) GlobalV::CURRENT_SPIN = GlobalC::kv.isk [ik];
684- GlobalC::wf. npw = GlobalC::kv.ngk [ik];
684+ const int nbasis = GlobalC::kv.ngk [ik];
685685 // generate vkb
686686 if (GlobalC::ppcell.nkb > 0 )
687687 {
@@ -697,9 +697,11 @@ void Forces::cal_force_nl(ModuleBase::matrix& forcenl)
697697 {
698698 for (int i=0 ;i<nkb;i++)
699699 {
700- for (int ig=0 ; ig<GlobalC::wf.npw ; ig++)
700+ const std::complex <double >* ppsi = &(GlobalC::wf.psi [0 ](ik, ib, 0 ));
701+ const std::complex <double >* pvkb = &(GlobalC::ppcell.vkb (i, 0 ));
702+ for (int ig=0 ; ig<nbasis; ig++)
701703 {
702- becp (i,ib) += GlobalC::wf. evc [ik](ib,ig) * conj ( GlobalC::ppcell. vkb (i,ig) );
704+ becp (i,ib) += ppsi[ig] * conj ( pvkb[ig] );
703705 }
704706 }
705707 }
@@ -715,17 +717,17 @@ void Forces::cal_force_nl(ModuleBase::matrix& forcenl)
715717 {
716718 if (ipol==0 )
717719 {
718- for (int ig=0 ; ig<GlobalC::wf. npw ; ig++)
720+ for (int ig=0 ; ig<nbasis ; ig++)
719721 vkb1 (i, ig) = GlobalC::ppcell.vkb (i, ig) * ModuleBase::NEG_IMAG_UNIT * GlobalC::pw.get_G_cartesian_projection (GlobalC::wf.igk (ik, ig), 0 );
720722 }
721723 if (ipol==1 )
722724 {
723- for (int ig=0 ; ig<GlobalC::wf. npw ; ig++)
725+ for (int ig=0 ; ig<nbasis ; ig++)
724726 vkb1 (i, ig) = GlobalC::ppcell.vkb (i, ig) * ModuleBase::NEG_IMAG_UNIT * GlobalC::pw.get_G_cartesian_projection (GlobalC::wf.igk (ik, ig), 1 );
725727 }
726728 if (ipol==2 )
727729 {
728- for (int ig=0 ; ig<GlobalC::wf. npw ; ig++)
730+ for (int ig=0 ; ig<nbasis ; ig++)
729731 vkb1 (i, ig) = GlobalC::ppcell.vkb (i, ig) * ModuleBase::NEG_IMAG_UNIT * GlobalC::pw.get_G_cartesian_projection (GlobalC::wf.igk (ik, ig), 2 );
730732 }
731733 }
@@ -737,9 +739,11 @@ void Forces::cal_force_nl(ModuleBase::matrix& forcenl)
737739 if (GlobalC::wf.wg (ik, ib) < ModuleBase::threshold_wg) continue ;
738740 for (int i=0 ; i<nkb; i++)
739741 {
740- for (int ig=0 ; ig<GlobalC::wf.npw ; ig++)
742+ const std::complex <double >* ppsi = &(GlobalC::wf.psi [0 ](ik, ib, 0 ));
743+ const std::complex <double >* pvkb1 = &(vkb1 (i, 0 ));
744+ for (int ig=0 ; ig<nbasis; ig++)
741745 {
742- dbecp (i,ib, ipol) += conj ( vkb1 (i,ig) ) * GlobalC::wf. evc [ik](ib,ig) ;
746+ dbecp (i,ib, ipol) += conj ( pvkb1[ig] ) * ppsi[ig] ;
743747 }
744748 }
745749 }
0 commit comments