C++ VectorType::adjoint方法代码示例

typename GaussianProcess<TScalarType>::VectorType
GaussianProcess<TScalarType>::Predict(const typename GaussianProcess<TScalarType>::VectorType &x){
    Initialize();
    CheckInputDimension(x, "GaussianProcess::Predict: ");
    VectorType Kx;
    ComputeKernelVector(x, Kx);
    return (Kx.adjoint() * m_RegressionVectors).adjoint();
}

int main(int argc, char* argv[]){

    // gather inputs
    if (argc<2){
        cout << "please specify input" << endl;
        exit(0);
    }
    InputManager manager(argv[1]);
    int nstep = atoi( manager["CPMD"]["nstep"].c_str() );
    double dt = atof( manager["CPMD"]["dt"].c_str() );
    double emass = atof( manager["CPMD"]["emass"].c_str() );
    double thres = atof( manager["CPMD"]["conv_thr"].c_str() );
    double L = atof( manager["DFT"]["L"].c_str() );
    double Ecut = atof( manager["DFT"]["Ecut"].c_str() );
    int nbasis = atoi( manager["DFT"]["maxBasis"].c_str() );
    int nx = atoi( manager["DFT"]["ngrid"].c_str() );
    string traj=manager["output"]["trajectory"]; // Output files

    // Put a hydrogen ion at the origin
    ParticlePool pPool(1);
    ParticleSet gPset(pPool.myParticles());
    
    gPset.ptcls[0]->name="H"; 
    gPset.ptcls[0]->m=1836;
    gPset.ptcls[0]->q=1;
    cout << gPset.str() << endl;
    
    // choose a basis set
    BasisSet waveFunctionBasis(nbasis,L);
    BasisSet densityBasis(pow(2,_DFT_DIM)*nbasis,L);
    waveFunctionBasis.initPlaneWaves(Ecut);
    densityBasis.initPlaneWaves(2*Ecut);
    cout << "Number of Wave Function Basis: " << waveFunctionBasis.size() << endl;
    cout << "Number of Density Basis: " << densityBasis.size() << endl;
    
    // construct external potential and Hamiltonian
    ExternalPotential   Vext(&densityBasis,gPset);
    Vext.initGrid(L,nx);
    Vext.updatePlaneWaves(); // obtain expansion coefficients for Vext
    Hamiltonian H(&waveFunctionBasis);
    H.update(Vext); // obtain expansion coefficients for Hamiltonian
    
    // diagonalize the Hamiltonian
    Eigen::SelfAdjointEigenSolver<MatrixType> eigensolver(*H.myHam());
    ComplexType lambda=eigensolver.eigenvalues()[0]; // Lagrange Multiplier
    cout << "The lowest eigenvalue of H is: " << lambda.real() << endl;
    VectorType c = eigensolver.eigenvectors().col(0);
    
    // mess up the electronic ground state a little
    RealType E,norm;
    c = c.Random(c.size())/100.+c;
    c /= c.norm();
    E = ( c.adjoint()*(*H.myHam())*c )(0,0).real();
    cout << "Energy after slight electronic perturbation: " << E << endl;
    cout << "Starting energy minimization " << endl;
    cout << "(Energy, Norm)" << endl; 
    // See that CP can bring back the ground state
    VectorType oldc=c;
    RealType oldE=E;
    for (int istep=0;istep<nstep;istep++){
    
        // update Hamiltonian
        Vext.updatePlaneWaves();
        H.update(Vext);
    
        // SHAKE it for the correct lambda
        VectorType uc; // uncontraint c
        uc = 2*c-oldc-pow(dt,2)/emass*( (*H.myHam())*c -lambda*c );
        ComplexType usigma(-1.0,0.0);
        for (int i=0;i<uc.size();i++){
            usigma += conj( uc[i] )*uc[i];
        }
        ComplexType dot(0.0,0.0);
        for (int i=0;i<uc.size();i++){
            dot += conj( uc[i] )*conj( c[i] )/(ComplexType)emass;
        }
        lambda = usigma/dot;
        
        // move the electrons
        for (int i=0;i<uc.size();i++){
            c[i] = uc[i]-conj( c[i] )*lambda*(ComplexType)(dt/emass);
        }
        
        // report properties
        norm = c.norm();
        E = ( c.adjoint()*(*H.myHam())*c )(0,0).real();
        cout << "( " << E << "," << norm << ")" << endl; 
        
        if (abs(E-oldE)<thres) break;
        
        oldc = c; 
        oldE = E;
    }
    
    // ---------- we are now ready to do MD ----------
    
    // create lowest KS orbital and its associated density
    Function waveFunction(&waveFunctionBasis);
    Density density(&densityBasis);
    
//.........这里部分代码省略.........