@@ -164,13 +164,6 @@ void Stochastic_Iter::itermu(int &iter)
164164 double ne3;
165165 double mu3;
166166
167- // test the domin of mu
168- /* for(mu = -5; mu<5;mu+=0.2)
169- {
170- ne3 = calne();
171- cout<<"mu: "<<mu<<" ; ne: "<<ne3<<endl;
172- }
173- exit(0);*/
174167 while (ne1 > targetne)
175168 {
176169 ne2 = ne1;
@@ -208,8 +201,7 @@ void Stochastic_Iter::itermu(int &iter)
208201 mu2 = mu3;
209202 }
210203 Dne = abs (targetne - ne3);
211- // cout<<setw(20)<<"targetne"<<setw(20)<<"ne"<<setw(20)<<"mu"<<setw(20)<<"Dne"<<endl;
212- // cout<<setw(20)<<targetne<<setw(20)<<ne3<<setw(20)<<mu3<<setw(20)<<Dne<<endl;
204+
213205 count++;
214206 if (count > 60 )
215207 {
@@ -388,11 +380,7 @@ void Stochastic_Iter::sum_stoband(Stochastic_WF& stowf)
388380 pchi = stowf.chi0 [ik].c ;
389381
390382 stoche.calfinalvec (stohchi.hchi_reciprocal , pchi, out, nchip[ik]);
391- // hout = new complex<double> [npwall];
392- // stohchi.hchi_reciprocal(out,hout,nchip);
393- // stok_eband = Diago_CG::ddot_real(npwall, out, hout,false) * DeltaE
394- // + Diago_CG::ddot_real(npwall, out, out,false) * Ebar;
395- // sto_eband += stok_eband * kv.wk[ik];
383+
396384 std::complex <double > *tmpout = out;
397385 for (int ichi = 0 ; ichi < nchip[ik] ; ++ichi)
398386 {
@@ -408,33 +396,8 @@ void Stochastic_Iter::sum_stoband(Stochastic_WF& stowf)
408396 }
409397 tmpout+=npwx;
410398 }
411- // delete [] hout;
412399 }
413- // for(int ip = 0 ; ip < NPOOL ; ++ip)
414- // {
415- // MPI_Barrier(MPI_COMM_WORLD);
416- // if(ip == MY_POOL)
417- // {
418- // if(RANK_IN_POOL==0)
419- // {
420- // cout.clear();
421- // for(int ik = 0; ik < kv.nks; ++ik)
422- // {
423- // complex<double> *tmpout = STO_WF.shchi[ik].c;
424- // for(int ichi = 0; ichi < nchip[ik] ; ++ichi)
425- // {
426- // for(int ig = 0; ig < kv.ngk[ik]; ++ig)
427- // {
428- // if(ig%100==0) cout<<tmpout[ig]<<" ";
429- // }
430- // cout<<endl;
431- // tmpout+=npwx;
432- // }
433- // }
434- // if(MY_RANK!=0) cout.setstate(ios::failbit);
435- // }
436- // }
437- // }
400+
438401 GlobalC::CHR.rho_mpi ();
439402 for (int ir = 0 ; ir < nrxx ; ++ir)
440403 {
@@ -571,237 +534,3 @@ double Stochastic_Iter:: nfdlnfd(double e)
571534 return f * log (f) + (1 -f) * log (1 -f);
572535 }
573536}
574-
575-
576-
577- // void Stochastic_Iter:: test(const int & ik)
578- // {
579- // =====================test============================
580- /*
581- complex<double> *in = new complex<double> [pw.nrxx];
582-
583- complex<double> *chig1 = new complex<double> [wf.npw];
584- complex<double> *chig2 = new complex<double> [wf.npw];
585- ZEROS(in,pw.nrxx);
586- ZEROS(in2,pw.nrxx);*/
587-
588- // ---------------------------------------------------
589- /* //test hermit property of hchi matrix
590- Emin = -1;
591- Emax = 1;
592- Stochastic_hchi:: Emin = this -> Emin;
593- Stochastic_hchi:: Emax = this -> Emax;
594- complex<double> *out = new complex<double> [pw.nrxx];
595- complex<double> *in2 = new complex<double> [pw.nrxx];
596- cout<<"------------------------------------"<<endl;
597- complex<double> cij,cji;
598- double dc;
599- for(int i = 0 ; i < 300 ; ++i)
600- {
601- if( i % 10 == 0)
602- cout<<"We are testing "<<i+1 <<" rows..."<<endl;
603- for(int j = i+1; j < 300 ; ++j)
604- {
605- in2[j] = 1;
606- stohchi.hchi_real(in2,out);
607- cij = out[i];
608- in2[j] = 0;
609- in2[i] = 1;
610- stohchi.hchi_real(in2,out);
611- cji = out[j];
612- in2[i] = 0;
613- dc = abs(conj(cij)-cji);
614- if(dc > 1e-6)
615- {
616- cout<<"(i,j) = ("<<i+1<<" , "<<j+1<<") ; cij = "<<cij<<" ; cji = "<<cji<<endl;
617- }
618-
619- }
620- }
621- cout<<"------------------------------------"<<endl;
622- delete[] out;
623- delete[] in2;*/
624- // ---------------------------------------------------
625-
626- // -------------------------------------------------------
627- // compare hchi_reciprocal and h_psi
628- /* Emin = -1;
629- Emax = 1;
630- Stochastic_hchi:: Emin = this -> Emin;
631- Stochastic_hchi:: Emax = this -> Emax;
632-
633- int m = 5;
634- complex<double> *chig1 = new complex<double> [wf.npwx*m];
635- complex<double> *chig2 = new complex<double> [wf.npwx*m];
636- complex<double> *hchig = new complex<double> [wf.npwx*m];
637- //complex<double> *kswf = &wf.evc[ik](0,0);
638- complex<double> *kswf = STO_WF.chi0[ik].c;
639-
640- stohchi.hchi_reciprocal(kswf,chig1,m);
641- hm.hpw.h_psi( kswf , chig2, m);
642- cout<<"==================ik="<<ik<<"====================================="<<endl;
643- if(MY_RANK==0)
644- for(int ib = 0 ; ib < m ; ++ib)
645- {
646- cout<<wf.npw<<" "<<wf.npwx<<endl;
647- for(int i = 0; i<wf.npw;++i)
648- {
649- if(i % 100 == 0 )
650- cout<<kswf[i+ib*wf.npwx]<<" "<<chig1[i+ib*wf.npwx]<<" "<<chig2[i+ib*wf.npwx]<<endl;
651- }
652- cout<<"-------------------------------------------------------------"<<endl;
653- }
654- cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
655- MPI_Barrier(MPI_COMM_WORLD);
656-
657- if(MY_RANK==1)
658- {
659- cout.clear();
660-
661- for(int ib = 0 ; ib < m ; ++ib)
662- {
663- cout<<wf.npw<<" "<<wf.npwx<<endl;
664- for(int i = 0; i<wf.npw;++i)
665- {
666- if(i % 100 == 0 )
667- cout<<kswf[i+ib*wf.npwx]<<" "<<chig1[i+ib*wf.npwx]<<" "<<chig2[i+ib*wf.npwx]<<endl;
668- }
669- cout<<"-------------------------------------------------------------"<<endl;
670- }
671- cout.setstate(ios::failbit);
672- }
673- MPI_Barrier(MPI_COMM_WORLD);
674- if(MY_RANK==3)
675- {
676- cout.clear();
677-
678- for(int ib = 0 ; ib < m ; ++ib)
679- {
680- cout<<wf.npw<<" "<<wf.npwx<<endl;
681- for(int i = 0; i<wf.npw;++i)
682- {
683- if(i % 100 == 0 )
684- cout<<kswf[i+ib*wf.npwx]<<" "<<chig1[i+ib*wf.npwx]<<" "<<chig2[i+ib*wf.npwx]<<endl;
685- }
686- cout<<"-------------------------------------------------------------"<<endl;
687- }
688- cout.setstate(ios::failbit);
689- }
690- MPI_Barrier(MPI_COMM_WORLD);
691- cout<<"================================================================"<<endl;
692- cout<<endl;
693- delete[] chig1;
694- delete[] chig2;
695- if(ik == kv.nks-1) exit(0);*/
696- // ---------------------------------------------------
697-
698- // ---------------------------------------------------
699- // compare eigen energy
700- /*
701- int * GRA_index = stohchi.GRA_index;
702- complex<double> *chigout = new complex<double> [wf.npw];
703- complex<double> *wave = new complex<double> [pw.nrxx];
704- complex<double> *waveout = new complex<double> [pw.nrxx];
705- Emax = 1000;
706- Emin = 0;
707- Stochastic_hchi:: Emin = this->Emin;
708- Stochastic_hchi:: Emax = this->Emax;
709- double Ebar = (Emax + Emin)/2;
710- double DeltaE = (Emax - Emin)/2;
711- fftw_plan pp=fftw_plan_dft_3d(pw.nx,pw.ny,pw.nz,(fftw_complex *)wave,(fftw_complex *)wave, FFTW_BACKWARD, FFTW_ESTIMATE);
712- for(int ib = 0 ; ib < NBANDS ; ++ib)
713- {
714- complex<double> *kswf = &wf.evc[0](ib,0);
715- hm.hpw.h_psi( kswf , chigout);
716- double energy = 0;
717- double norm1 =0;
718- ZEROS(wave,pw.nrxx);
719- for(int ig = 0 ; ig < wf.npw ; ++ig)
720- {
721- energy += real(conj(kswf[ig]) * chigout[ig]);
722- norm1 += norm (kswf[ig]);
723- wave[GRA_index[ig]] = kswf[ig];
724- }
725- fftw_execute(pp);
726- stohchi.hchi_real(wave,waveout);
727- double energy2 = 0;
728- double norm2 =0;
729- for(int ir = 0 ; ir < pw.nrxx ; ++ir)
730- {
731- energy2 += real(conj(wave[ir]) * waveout[ir]) * DeltaE + Ebar * norm(wave[ir]);
732- norm2 += norm(wave[ir]);
733- }
734-
735-
736- cout<<"BAND "<<ib+1<<" "<<energy/norm1<<" "<<energy2/norm2<<endl;
737- }
738- delete []chigout;
739- delete []wave;
740- delete []waveout;*/
741- // ---------------------------------------------------
742-
743- // ---------------------------------------------------
744- // test ne
745- /*
746- int * GRA_index = stohchi.GRA_index;
747- complex<double> *wave = new complex<double> [pw.nrxx];
748- complex<double> *waveout = new complex<double> [pw.nrxx];
749- Emax = 750;
750- Emin = -100;
751- Stochastic_hchi:: Emin = this->Emin;
752- Stochastic_hchi:: Emax = this->Emax;
753- mu = en.ef;
754- stoche.calcoef(this->nfd);
755- fftw_plan pp=fftw_plan_dft_3d(pw.nx,pw.ny,pw.nz,(fftw_complex *)wave,(fftw_complex *)wave, FFTW_BACKWARD, FFTW_ESTIMATE);
756- for(int ib = 0 ; ib < NBANDS ; ++ib)
757- {
758- ZEROS(wave,pw.nrxx);
759- complex<double> *kswf = &wf.evc[0](ib,0);
760- for(int ig = 0 ; ig < wf.npw ; ++ig)
761- {
762- wave[GRA_index[ig]] = kswf[ig];
763- }
764- fftw_execute(pp);
765- double ne =0;
766- double norm1 = 0;
767- stoche.calresult(stohchi.hchi_real, pw.nrxx, wave, waveout);
768- for(int ir = 0 ; ir < pw.nrxx ; ++ir)
769- {
770- ne += real(conj(wave[ir]) * waveout[ir]);
771- norm1 += norm(wave[ir]);
772- }
773- cout<<"Ne of Band "<<ib+1<<" is "<<ne/norm1 * kv.wk[0]<<endl;
774- }
775- delete []wave;
776- delete []waveout;*/
777-
778- // -------------------------------------------------------------
779- // test orthogonal
780- /* int npw=wf.npw;
781- char transC='C';
782- char transN='N';
783- int nchip = STO_WF.nchip;
784- complex<double> *wave = new complex<double> [nchip*npw];
785- for(int i = 0; i < nchip*npw; ++i)
786- {
787- wave[i]=STO_WF.chi0[0].c[i];
788- }
789- complex<double> *sum = new complex<double> [nchip * nchip];
790- zgemm_(&transC, &transN, &nchip, &nchip, &npw, &ONE, STO_WF.chi0[0].c, &npw, wave, &npw, &ZERO, sum, &nchip);
791- Parallel_Reduce::reduce_complex_double_pool(sum, nchip * nchip);
792- if(MY_RANK!=0) cout.clear();
793- double abs2 = 0;
794- for(int i=0;i<nchip * nchip;++i)
795- {
796- if(i%nchip==int(i%nchip)) continue;
797- abs2+=norm(sum[i]);
798- }
799- cout<<abs2/nchip<<endl;
800- delete [] sum;
801- delete [] wave;
802- if(MY_RANK!=0) cout.setstate(ios::failbit);*/
803- // -------------------------------------------------------------
804-
805-
806- // =====================================================
807- // }
0 commit comments