C++ Array2::GetDim方法代码示例

void Wf_return::setVals(Array2 <dcomplex> & vals ,
                        Array1 <doublevar> &  p) {
  
  // here we extract amplitude and phase, and their gradients and laplacians,
  // from the complex value and its gradient and derivative
  
  is_complex=1;
  cvals=vals;
  int ntype=vals.GetDim(1);
  int nwf=vals.GetDim(0);
  for (int w=0; w< nwf; w++) {
    amp(w,0)=vals(w,0).real();
    phase(w,0)=p(w);
    
    doublevar sum_ii=0;
    doublevar sum_ri=0;
    if (ntype>=4) {
      for (int i=1; i<4; i++) {
        amp(w,i)=vals(w,i).real();
        phase(w,i)=vals(w,i).imag();
        sum_ii+=phase(w,i)*phase(w,i);
        sum_ri+=amp(w,i)*phase(w,i);
      }
    }
    
    if (ntype>=5) {
      amp(w,4)=vals(w,4).real()+sum_ii;
      phase(w,4)=vals(w,4).imag()-2*sum_ri;
    }
  }
  
}

void Minimization_function::init_overlap_ma_R(Array2<doublevar> &overlap_ma_R){
  assert(overlap_ma_R.GetDim(0)==overlap_ma_R.GetDim(1));
  _overlap_ma_R.Resize(overlap_ma_R.GetDim(0),overlap_ma_R.GetDim(0));
  for(int k=0;k<overlap_ma_R.GetDim(0);k++)
    for(int l=0;l<=k;l++)
      _overlap_ma_R(l,k)=_overlap_ma_R(k,l)=overlap_ma_R(k,l);
}

void recenter(Array2 <doublevar> & pos, Array1 <doublevar> & mass) {
  int natoms=pos.GetDim(0);
  assert(pos.GetDim(0)==mass.GetDim(0));
  assert(pos.GetDim(1)==3);

  Array1 <doublevar> center(3,0.0); //center of mass
  doublevar tot=0.0;
  for(int at=0; at< natoms; at++) {
    for(int d=0; d< 3; d++) {
      center(d)+=mass(at)*pos(at,d);
    }
    tot+=mass(at);
  }

  for(int d=0; d< 3; d++) {
    center(d)/=tot;
  }

  for(int at=0; at< natoms; at++) {
    for(int d=0; d< 3; d++) {
     pos(at,d)-=center(d);
     }
  }

}

inline void output_array(Array2 <doublevar> & arr) {
  for(int i=0; i< arr.GetDim(0); i++) {
    for(int j=0; j < arr.GetDim(1); j++) {
      cout << arr(i,j) << "  ";
    }
    cout << endl;
  }
}

void Cutoff_cusp::calcLap(
  const Array1 <doublevar> & r,
  Array2 <doublevar> & symvals,
  const int startfill
)
{
  assert(r.GetDim(0) >=5);
  assert(symvals.GetDim(0) >= 1+startfill);
  assert(symvals.GetDim(1) >= 5);
  if(r(0) < rcut) {
    doublevar zz=r(0)*rcutinv;
    doublevar zz2=zz*zz;
    doublevar pp=zz-zz2+zz*zz2/3;
    doublevar pade=1./(1+gamma*pp);
    symvals(startfill,0)=cusp*rcut*(pp*pade-pade0);
    doublevar pade2=pade*pade;
    doublevar ppd=1.-2.*zz+zz2;
    doublevar ppdd=-2.+2.*zz;
    doublevar dadr=ppd*pade2/r(0);
    doublevar dadr2=ppdd*pade2*rcutinv
                    -2.*gamma*ppd*ppd*pade2*pade*rcutinv
                    +2.*dadr;

    for(int i=1; i< 4; i++) {
      symvals(startfill, i)=cusp*dadr*r(i+1);
    }
    symvals(startfill, 4)=cusp*dadr2;
  }
  else {
    for(int i=0; i< 5; i++)
      symvals(startfill, i)=0;
  }
}

void MO_matrix_blas::writeorb(ostream & os, Array2 <doublevar> & rotation, Array1 <int>  &moList) {


  int nmo_write=moList.GetDim(0);
  assert(rotation.GetDim(0)==nmo_write);
  assert(rotation.GetDim(1)==nmo_write);
  os.precision(15);
  int counter=0;
  for(int m=0; m < nmo_write; m++) {
    int mo=moList(m);
    for(int ion=0; ion<centers.size(); ion++)
    {
      int f=0;

      for(int n=0; n< centers.nbasis(ion); n++)
      {
        int fnum=centers.basis(ion,n);
        int imax=basis(fnum)->nfunc();

        for(int i=0; i<imax; i++)
        {
          os << mo+1 << "  "   << f+1 << "   " << ion+1 << "   " << counter+1 << endl;
          f++;  //keep a total of functions on center
          counter++;
        } //i
      } //n
    }  //ion
  }
  os << "COEFFICIENTS\n";
  ifstream orbin(orbfile.c_str());
  rotate_orb(orbin, os, rotation, moList, totbasis);
  orbin.close();


}

doublevar Wannier_method::eval_tstep(Array3 <dcomplex> & eikr, Array2 <doublevar> & Rgen,
                                     Array2 <doublevar> & Rgen_save, Array2 <doublevar> & deriv, doublevar tstep,
                                     Array2 <doublevar> & R) {
    int norb=Rgen.GetDim(0);
    for(int ii=0; ii< norb; ii++) {
        for(int jj=ii+1; jj < norb; jj++) {
            Rgen(ii,jj)=Rgen_save(ii,jj)-tstep*deriv(ii,jj);
        }
    }
    return evaluate_local(eikr,Rgen,R);
}

void MO_matrix_standard::writeorb(ostream & os, 
                                  Array2 <doublevar> & rotation, 
                                  Array1 <int>  &moList) {

  int nmo_write=moList.GetDim(0);
  assert(rotation.GetDim(0)==nmo_write);
  assert(rotation.GetDim(1)==nmo_write);
  cout << "nmo_write " << nmo_write << endl;
 
  int counter=0;
  for(int m=0; m < nmo_write; m++) {
    int mo=moList(m);
    for(int ion=0; ion<centers.size(); ion++)
    {
      int f=0;

      for(int n=0; n< centers.nbasis(ion); n++)
      {
        int fnum=centers.basis(ion,n);
        int imax=basis(fnum)->nfunc();

        for(int i=0; i<imax; i++)
        {
          os << mo+1 << "  "   << f+1 << "   " << ion+1 << "   " 
             << counter+1 << endl;
          f++;  //keep a total of functions on center
          counter++;
        } //i
      } //n
    }  //ion
  }
  os << "COEFFICIENTS\n";

  Array2 <doublevar> rotated_mo(nmo_write, totbasis);
  rotated_mo=0;
  for(int m=0; m< nmo_write; m++) {
    for(int m2=0; m2 < nmo_write; m2++) {
      for(int f=0; f< totbasis; f++) {
        rotated_mo(m, f)+=rotation(m,m2)*moCoeff(m2,f);
      }
    }
  }

  int counter2=1;
  for(int m=0; m < nmo_write; m++) {
    for(int f=0; f< totbasis; f++) {
      os << rotated_mo(m, f) << "   ";
      if(counter2 % 5 ==0) os << endl;
    }
  }

}

void MO_matrix_standard::updateLap(
  Sample_point * sample, int e,  int listnum,
  Array2 <doublevar> & newvals
  //!< The return: in form (MO, [val, grad, lap])
){
  Array2 <doublevar> basisvals(totbasis, 5);
  int centermax=centers.size();
  centers.updateDistance(e, sample);
  Basis_function * tempbasis;
  Array1 <doublevar> R(5);
  int currfunc=0;
  newvals=0;
  //cout << "centermax " << centermax << endl;
  for(int ion=0; ion < centermax; ion++) {
    centers.getDistance(e, ion, R);
    //cout << "nbasis " << centers.nbasis(ion) << endl;
    for(int n=0; n< centers.nbasis(ion); n++)
    {

      tempbasis=basis(centers.basis(ion,n));
      tempbasis->calcLap(R, basisvals, currfunc);
      currfunc+=tempbasis->nfunc();
    }
  }

  doublevar c;
  int nmo_list=moLists(listnum).GetDim(0);
  int valscale=newvals.GetDim(1);
  int basisscale=basisvals.GetDim(1);
  int mocoeffscale=moCoeff.GetDim(1);
  for(int m=0; m < nmo_list; m++) {
    int mo=moLists(listnum)(m);
    int effmo=mo*mocoeffscale;
    int effm=m*valscale;
    int efff;
    for(int f=0; f< totbasis; f++) {
      //c=moCoeff(mo,f);
      //cout << " c " << c << endl;
      c=moCoeff.v[effmo+f];
      //cout << "c2 " << c << endl;
      efff=f*basisscale;
      for(int d=0; d< 5; d++) {
        //cout << "basisvals " << basisvals(f,d) << endl;
        //cout << "basisvals2 " << basisvals.v[efff+d] << endl;
        //
        //newvals(m,d)+=c*basisvals(f,d);
        //newvals.v[effm+d]+=c*basisvals(f,d);
        newvals.v[effm+d]+=c*basisvals.v[efff+d];
      }
    }
  }
}

//----------------------------------------------------------------------
//
//
//
doublevar Wannier_method::evaluate_local(const Array3 <dcomplex> & eikr,
        Array2 <doublevar> & Rgen, Array2 <doublevar> & R) {
    int norb=Rgen.GetDim(0);
    Array2 <doublevar> gvec(3,3);
    sys->getPrimRecipLattice(gvec);
    Array1 <doublevar> gnorm(3);
    gnorm=0;
    for(int d=0; d < 3; d++) {
        for(int d1=0; d1 < 3; d1++) {
            gnorm(d)+=gvec(d,d1)*gvec(d,d1);
        }
        gnorm(d)=sqrt(gnorm(d));
    }

    Array2 <dcomplex> tmp(norb,norb),tmp2(norb,norb);

    make_rotation_matrix(Rgen,R);
    doublevar func=0;
    for(int d=0; d< 3; d++) {
        tmp=0.0;
        for(int i=0; i< norb; i++) {
            for(int j=0; j< norb; j++) {
                for(int k=0; k< norb; k++) {
                    tmp(i,k)+=eikr(d,i,j)*R(k,j);
                }
            }
        }
        tmp2=0;
        for(int i=0; i< norb; i++) {
            for(int j=0; j< norb; j++) {
                for(int k=0; k< norb; k++) {
                    tmp2(i,k)+=R(i,j)*tmp(j,k);
                }
            }
        }
        //cout << "========for d=" << d << endl;

        for(int i=0; i< norb; i++)  {
            doublevar f=  -log(norm(tmp2(i,i)))/(gnorm(d)*gnorm(d));
            func+=f;

            // cout << setw(9) << f;
        }
    }
    //cout << endl;
    return func/(3*norb);

}

void Rgaussian_function::calcLap(
  const Array1 <doublevar> & r,
  Array2 <doublevar> & symvals,
  const int startfill){

  assert(symvals.GetDim(0) >= nmax+startfill);
  assert(r.GetDim(0) >= 5);
  doublevar exponent;
  //cout << "v_l for l=" << l << "   atom=" << at << endl;
  int index=startfill;

  for(int i=0; i< nmax; i++) {
    doublevar v_l=0, dv_l=0, d2v_l=0;

    for(int j=0; j< numExpansion(i); j++) {
      exponent=gaussexp(i, j)*r(1);
      exponent=min(exponent,(doublevar) 60.0);
      exponent=exp(-exponent);
      doublevar n=radiusexp(i,j);
      v_l+= gausscoeff(i, j)
            *pow(r(0),n )
            *exponent;
      dv_l+=gausscoeff(i,j)*exponent*pow(r(0), n-1)*( n
                                                      -2*r(1)*gaussexp(i,j));

      //note that the laplacian is untested..
      d2v_l+=gausscoeff(i,j)*pow(r(0), n-2)*( n*(n-1)-2*gaussexp(i,j)*n*r(1)
                                              -2*(n+1)*gaussexp(i,j)*r(1)
                                              +4*gaussexp(i,j)*gaussexp(i,j)
                                              *r(1)*r(1))*exponent;

    }
    dv_l/=r(0);

    symvals(index,0) = v_l;
    symvals(index,1) = dv_l*r(2);
    symvals(index,2) = dv_l*r(3);
    symvals(index,3) = dv_l*r(4);
    symvals(index,4) = d2v_l+2*dv_l;
    index++;
  }
  
}

void create_local_basis(string & centerin, string & basisin,
			Array1 <Center> & centers,
                        Array1 <Contracted_gaussian >  & basis ) {

  Array2 <doublevar> centerpos;
  vector <string> labels;


  read_centerpos(centerin, centerpos, labels);

  read_basis(basisin, basis);


  int ncenters=centerpos.GetDim(0);
  int nbasis=basis.GetDim(0);

  centers.Resize(ncenters);
  for(int cen=0; cen < ncenters; cen++) {
    for(int d=0; d< 3; d++) centers(cen).pos(d)=centerpos(cen,d);
    centers(cen).label=labels[cen];
  }



  vector <vector <int> > cenbasis, basiscen;
  cenbasis.resize(ncenters);
  basiscen.resize(nbasis);
  for(int cen=0; cen < ncenters; cen++) {
    int foundbas=0;
    for(int bas=0; bas < nbasis; bas++) {

      if(centers(cen).label==basis(bas).center_name) {
        foundbas=1;
        cenbasis[cen].push_back(bas);
        basiscen[bas].push_back(cen);
      }
    }
    if(!foundbas) 
      cout << "WARNING!  Couldn't find a basis set for center "
           << centers(cen).label << endl;
  }

  

  for(int bas=0; bas < nbasis; bas++) {
    int ncen=basiscen[bas].size();
    basis(bas).center.Resize(ncen);
    for(int c=0; c< ncen; c++) {
      basis(bas).center(c)=basiscen[bas][c];
      //cout << "basis " << bas << " ->  center " << basis(bas).center(c) << endl;
    }
  }

  for(int cen=0; cen < ncenters; cen++) {
    int nbas=cenbasis[cen].size();
    centers(cen).basis.Resize(nbas);
    for(int b=0; b< nbas; b++) {
      centers(cen).basis(b)=cenbasis[cen][b];
      //cout << "center " << cen << " -> basis " << centers(cen).basis(b) << endl;
    }
  }



  for(int bas=0; bas < nbasis; bas++) {
    int nalpha=basis(bas).alpha.GetDim(0);
    //cout << "basis " << bas << endl;
    for(int a =0; a < nalpha; a++) {
        int totL=basis(bas).lvals(0)+basis(bas).lvals(1)+basis(bas).lvals(2);
        doublevar exponent=basis(bas).alpha(a);
        doublevar fac=sqrt(2.*exponent/pi);
        doublevar feg=4.*exponent;
        doublevar feg2=feg*feg;
        doublevar norm=0;
        switch(totL)
        {
        case 0:
          norm=sqrt(2.*feg*fac);
          break;
        case 1:
          norm=sqrt(2.*feg2*fac/3.);
          break;
        case 2:
          norm=sqrt(2.*feg*feg2*fac/15.);
          break;
        case 3:
          norm=sqrt(2.*feg2*feg2*fac/105.);
          break;
        default:
          norm=0;
          error("Unknown symmetry in Cubic_spline::readbasis! Shouldn't be here!");
        }
        basis(bas).coeff(a)*=norm;
    }
  }


}

void MO_matrix_blas::updateLap(
  Sample_point * sample, 
  int e, 
  int listnum,
  Array2 <doublevar> & newvals) {

#ifdef USE_BLAS
  int centermax=centers.size();

  assert(e < sample->electronSize());
  assert(newvals.GetDim(1) >=5);


  Array1 <doublevar> R(5);
  static Array2 <doublevar> symmvals_temp(maxbasis,5);
  static Array2 <doublevar> newvals_T;
  Array2 <doublevar> & moCoefftmp(moCoeff_list(listnum));
  int totbasis=moCoefftmp.GetDim(0);
  int nmo_list=moCoefftmp.GetDim(1);
  newvals_T.Resize(5, nmo_list);
  newvals_T=0.0;

  int b;
  int totfunc=0;
  
  centers.updateDistance(e, sample);

  static Array1 <MOBLAS_CalcObjLap> calcobjs(totbasis);
  int ncalcobj=0;
  

  for(int ion=0; ion < centermax; ion++)  {
    centers.getDistance(e, ion, R);
    for(int n=0; n< centers.nbasis(ion); n++) {
      b=centers.basis(ion, n);
      if(R(0) < obj_cutoff(b)) {
        basis(b)->calcLap(R, symmvals_temp);
        int imax=nfunctions(b);
        for(int i=0; i< imax; i++) {
          if(R(0) < cutoff(totfunc)) {
            calcobjs(ncalcobj).moplace=moCoefftmp.v+totfunc*nmo_list;
            
            for(int j=0; j< 5; j++) {
              calcobjs(ncalcobj).sval[j]=symmvals_temp(i,j);
              //cblas_daxpy(nmo_list,symmvals_temp(i,j),
              //            moCoefftmp.v+totfunc*nmo_list,1,
              //            newvals_T.v+j*nmo_list,1);
            }           
            ncalcobj++;
          }
          totfunc++;
        }
      }
      else {
        totfunc+=nfunctions(b);
      }
    }
  }

  for(int i=0; i< ncalcobj; i++) { 
    for(int j=0; j< 5; j++) { 
      cblas_daxpy(nmo_list,calcobjs(i).sval[j],calcobjs(i).moplace,1,
          newvals_T.v+j*nmo_list,1);
    }
    //cblas_dger(CblasRowMajor,5,nmo_list,1.,
    //      calcobjs(i).sval,1,
    //      calcobjs(i).moplace,1,
    //      newvals_T.v,nmo_list);
  }
  

  for(int m=0; m < nmo_list; m++) {
    for(int j=0; j< 5; j++) {
      newvals(m,j)=newvals_T(j,m);
    }
  }

#else
  error("BLAS libraries not compiled in, so you cannot use BLAS_MO");
#endif

}

void get_onebody_parms_diffspin(vector<string> & words, vector<string> & atomnames,
				Array2 <doublevar> & unique_parameters,
				int s,
				Array1 <int> & nparms,
				vector <string> & parm_labels,
				Array2 <int> & linear_parms,
                                int & counter,
				Array1 <int> & parm_centers) {
  vector < vector < string> > parmtxt;
  vector <string> parmtmp;
  unsigned int pos=0;
  if(s==0){
    while(readsection(words, pos,parmtmp, "LIKE_COEFFICIENTS"))
      parmtxt.push_back(parmtmp);
    if(parmtxt.size()==0)
      error("Didn't find LIKE COEFFICIENTS in one of threebody spin Jastrow section");
  }
  else{
    while(readsection(words, pos,parmtmp, "UNLIKE_COEFFICIENTS"))
      parmtxt.push_back(parmtmp);
    if(parmtxt.size()==0)
      error("Didn't find UNLIKE COEFFICIENTS in one of threebody spin Jastrow section");
  }
    
  int nsec=parmtxt.size();
  int maxsize=0;
  nparms.Resize(nsec);
  for(int i=0; i< nsec; i++)
    nparms(i)=parmtxt[i].size()-1;

  for(int i=0; i< nsec; i++)
    if(maxsize < nparms(i)) maxsize=nparms(i);

  unique_parameters.Resize(nsec, maxsize);
  linear_parms.Resize(nsec, maxsize);
  //cout <<"spin "<<s<<" counter "<<counter<<endl;
  for(int i=0; i< nsec; i++) {
    for(int j=0; j< nparms(i); j++) {
      
      unique_parameters(i,j)=atof(parmtxt[i][j+1].c_str());
      linear_parms(i,j)=counter++;
    }
  }

  parm_labels.resize(nsec);
  for(int i=0; i< nsec; i++)
    parm_labels[i]=parmtxt[i][0];

  int natoms=atomnames.size();
  parm_centers.Resize(natoms);
  parm_centers=-1;

  for(int i=0; i< nsec; i++) {
    int found=0;
    for(int j=0; j< natoms; j++) {
      if(atomnames[j]==parm_labels[i]) {
        parm_centers(j)=i;
        found=1;
      }
    }
    if(!found)
      error("Couldn't find a matching center for ", parm_labels[i]);
  }

  doublevar scale;
  if(readvalue(words, pos=0, scale, "SCALE")) {
    for(int i=0; i< unique_parameters.GetDim(0); i++) {
      for(int j=0; j< nparms(i); j++) {
        unique_parameters(i,j)*=scale;
      }
    }
  }
}

void Properties_final_average::showSummary(ostream & os, Array2 <Average_generator*>  avg_gen) {
  int nwf=avg.GetDim(1);
  int naux=aux_energy.GetDim(0);
  int autocorr_depth=autocorr.GetDim(1);
  os << "Threw out the first " << threw_out << " blocks in the averaging." 
     << endl;

  if(show_autocorr) {
    os << "Autocorrelation of energy " << endl;
    for(int i=0; i< autocorr_depth; i++) {
      os << "  en_ac_d" << i << "   " << autocorr(0,i) 
         << " +/- " << sqrt(autocorrerr(0,i)) << endl;
    }
    
    os << endl << endl;
  }
   
  using namespace  Properties_types;

  doublevar indep_points=sqrt(avgvar(0,0))/sqrt(err(0,0));
  int field=os.precision()+5;
  for(int w=0; w< nwf; w++) {
    for(int i=0; i < NUM_QUANTITIES; i++) {
      string tmp=names[i];
      append_number(tmp, w);
      os <<  setw(15)<< left  <<tmp << "    "
         << right << setw(field) << avg(i,w)
         << " +/- " << setw(field) << sqrt(err(i,w))  
	 << " (sigma " << setw(field) << sqrt(avgvar(i,w)) << " ) " << endl;
    }

  }

  if(nwf > 1) {
    os << "Energy differences " << endl;
    for(int w =1; w < nwf; w++) 
      os << "  diff_" << w << "-0    " << diff_energy(w)
         << "  +/-   " << sqrt(diff_energyerr(w)) << endl;
  }
  
  if(naux > 0) 
    os << endl << endl << "Auxillary differences " << endl;
  assert(aux_size.GetDim(0)==naux);
  int n_cvg=aux_diff.GetDim(1);
  for(int i=0; i< naux; i++) {
    for(int w=0; w< n_cvg; w++) {
      os << "final_auxdiff" << i << "-" << w << "   "
      << aux_diff(i,w)/aux_size(i) << "  +/-   " 
      << sqrt(aux_differr(i,w))/aux_size(i) 
      << "   (sigma " << sqrt(aux_differr(i,w))*indep_points << " ) "
      << endl;
    }
  }
  os << endl;
  
 
  assert(avg_gen.GetDim(1)==avgavg.GetDim(1) && avgavg.GetDim(1)==avgerr.GetDim(1));
  for(int w=0; w< nwf; w++) { 
    for(int i=0; i< avg_gen.GetDim(1); i++) { 
      avg_gen(w,i)->write_summary(avgavg(w,i),avgerr(w,i), os);
    }
  }
 
  os << "--------Properties differences-------------\n";
  for(int w=1; w< nwf; w++) { 
    for(int i=0; i< avg_gen.GetDim(1); i++) { 
      avg_gen(w,i)->write_summary(avgdiffavg(w,i),avgdifferr(w,i), os);
    }
  }
  
  doublevar totpoints=indep_points*indep_points;
  os << "approximate number of independent points: " 
     << totpoints << endl;
  if(totpoints> totweight) {
    os << "Warning!  The estimated number of independent points is \n"
      " _greater_ than the points in the log file.  Error bars are \n"
      "probably unreliable!\n";
  }

}