1-RDM, 2-RDM and canonical orbitals

[cuanto]

Hi Sebastian,

I am interested in QTAIM and energy partitions over these domains. I need the 1-RDM and 2-RDM as well the canonical orbitals. To do it, I added new options in dmrg plugin, for example in order to print 1-RDM and 2-RDM

/* Print 1-RDM in MO basis */
if (( dmrg_print_opm ) && ( nIterations > 0)){

    std::ofstream opm;
    std::string filename_opm = get_writer_file_prefix( wfn->molecule()->name() ) + ".opm";
    outfile->Printf( "Write 1-RDM file to %s. \n", filename_opm.c_str() );
    outfile->Printf( "Tolerance for printing 1-RDM values %e. \n", dmrg_print_tol );
    opm.open( filename_opm.c_str() , ios::out | ios::trunc );

    int ij;
    double tr = 0.0;
    for (int i = 0; i < nOrbDMRG; i++ ){
        for (int j = 0; j < nOrbDMRG; j++ ){
            ij = i*nOrbDMRG+j;
            if ( i == j ) {
                tr += DMRG1DM[ij];
            }
            if ( abs( DMRG1DM[ij] ) > dmrg_print_tol ) {
                opm << std::setiosflags(ios::fixed) << std::setprecision(12) << i << " " 
                                                                             << j << " " 
                                                                             << DMRG1DM[ij] << endl;
            }
        }
    } 

    opm.close();
    outfile->Printf( "tr(1-RDM) = %20.12lf\n", tr );
}

/* Print 2-RDM in MO basis */
if (( dmrg_print_tpm ) && ( nIterations > 0)){

    std::ofstream tpm;
    std::string filename_tpm = get_writer_file_prefix( wfn->molecule()->name() ) + ".tpm";
    outfile->Printf( "Write 2-RDM file to %s. \n", filename_tpm.c_str() );
    outfile->Printf( "Tolerance for printing 2-RDM values %e. \n", dmrg_print_tol );
    tpm.open( filename_tpm.c_str() , ios::out | ios::trunc );

    int ij;
    int ijkl;
    double tr = 0.0;
    for (int i = 0; i < nOrbDMRG; i++ ){
        for (int j = 0; j < nOrbDMRG; j++ ){
            ij = i*nOrbDMRG*nOrbDMRG*nOrbDMRG+j*nOrbDMRG*nOrbDMRG+i*nOrbDMRG+j; 
            tr += DMRG2DM[ij];
            for (int k = 0; k < nOrbDMRG; k++ ){
                for (int l = 0; l < nOrbDMRG; l++ ){
                    /* ijkl -> ikjl permute to obtain 2-RDM in chemist format */
                    ijkl = i*nOrbDMRG*nOrbDMRG*nOrbDMRG+j*nOrbDMRG*nOrbDMRG+k*nOrbDMRG+l;
                    if ( abs ( DMRG2DM[ijkl] ) > dmrg_print_tol ){
                       tpm << std::setiosflags(ios::fixed) << std::setprecision(12) << i << " " 
                                                                                    << k << " " 
                                                                                    << j << " " 
                                                                                    << l << " " 
                                                                                    << DMRG2DM[ijkl] << endl;
                    }
                }
            }   
        }
    } 

    tpm.close();
    outfile->Printf( "tr(2-RDM) = %20.12lf\n", tr );

}

I run it without problems using the following example

sys.path.insert(0, '/home/cuanto/src/chemps2/integrals/psi4plugins/dmrg') import dmrg

molecule { O H 1 1.00 H 1 1.00 2 103.1 }

set basis cc-pvdz set reference rhf set scf_type PK set e_convergence 1e-12 set d_convergence 1e-12

scf_e, scf_wfn = energy('scf', return_wfn=True) molden(scf_wfn, 'h2o.mol')

set dmrg wfn_irrep 0 set dmrg wfn_multp 1 set dmrg frozen_docc [1, 0, 0, 0] set dmrg active [3, 0, 1, 2]

set dmrg dmrg_states [ 500, 1000, 1000 ] set dmrg dmrg_e_convergence [ 1e-10, 1e-10, 1e-10 ] set dmrg dmrg_maxsweeps [ 15, 25, 110 ] set dmrg dmrg_noiseprefactors [ 0.05, 0.05, 0.0 ] set dmrg dmrg_dvdson_rtol [ 1e-5, 1e-5, 1e-9 ] set dmrg dmrg_print_corr true set dmrg dmrg_chkpt false

set dmrg d_convergence 1e-6 set dmrg dmrg_store_unit true set dmrg dmrg_do_diis true set dmrg dmrg_diis_branch 1e-2 set dmrg dmrg_store_diis true

set dmrg dmrg_max_iter 100 set dmrg dmrg_which_root 1 # Ground state set dmrg dmrg_state_avg false set dmrg dmrg_active_space INPUT # INPUT; NO; LOC set dmrg dmrg_loc_random false

set dmrg dmrg_molden true # Converged DMRG-SCF pseudocanonical orbitals set dmrg dmrg_density_ao false # Converged DMRG-SCF density in the AO basis

set dmrg dmrg_print_opm true set dmrg dmrg_print_tpm true set dmrg dmrg_print_tol 1e-12

dmrg_e, dmrg_wfn = energy('dmrg', return_wfn=True) molden(dmrg_wfn, 'h2o_dmrg.mol')

for example the 1-RDM read

0 0 1.991989837168 0 1 0.010598819956 0 2 -0.003445956475 1 0 0.010598819956 1 1 1.982327263575 1 2 0.024202834780 2 0 -0.003445956475 2 1 0.024202834780 2 2 0.025841279686 3 3 1.999368821701 4 4 1.973666420668 4 5 -0.013880082508 5 4 -0.013880082508 5 5 0.026806377203

but I think the indexes are wrong or in some strange order, for example orbital 2 has a low population, so If i use it to integrate (of course I have to add the core part) I can not reproduce the energy value.

By now I can do it using pyscf and the order is what I expect.

Next question concern the canonical orbitals, if I use these lines

dmrg_e, dmrg_wfn = energy('dmrg', return_wfn=True) molden(dmrg_wfn, 'h2o_dmrg.mol')

return the HF orbitals, If i do DMRG-CI I could use it to integrate in my real space code, but if I do DMRG-CASSCF I need the canonical rotated orbitals. I would like to add a new order to print it like

/* Print molden file with canonical orbitals */ if (( dmrg_molden_ca ) && ( nIterations > 0)){

    boost::shared_ptr<MoldenWriter> molden_ca( new MoldenWriter( wfn ) );
    std::string filename_ca = get_writer_file_prefix( wfn->molecule()->name() ) + ".canonical.molden";
    outfile->Printf( "Write canonical molden file to %s. \n", filename_ca.c_str() );
    molden_ca->write( filename_ca, var1, var2, var3, var4, var5, var6);

}

Could you help me, what I should pass in var1, var2, .... etc ?

Best Regards, Jose Luis

[SebWouters]

Hi Jose,

For the first part, the DMRG 1-RDM and DMRG 2-RDM are with respect to the CASSCF (DMRG-SCF) active space orbitals. If you want the 1-RDM and 2-RDM in the pseudocanonical orbitals, you should add the printer after https://github.com/SebWouters/CheMPS2/blob/9fd492793c91aeb5c3e04854aca28622abdecdfb/integrals/psi4plugins/dmrg.cc#L974 (lines 974 - 992).

1. With which integrals do you contract the RDMs to obtain the energy? Since the plugin computes correct DMRG-SCF and DMRG-CASPT2 energies, I think there's just a convention misunderstanding.

2. With

   set dmrg frozen_docc [1, 0, 0, 0] set dmrg active [3, 0, 1, 2]

active space orbital number 3 (or 2 in c++ count) will be the orbital of irrep A1 with the highest single particle energy still included in the active space. With the specifications above, four A1 orbitals can be occupied (the lowest one ic certainly doubly occupied as it's frozen_docc). Take a look at the lines "NORB = "; "NOCC = "; "NDMRG = "; and "NVIRT = " in the output. Per irrep, the latter three add up to the NORB value. So when you count the active space orbitals from 0 to 5, their irreps are A1, A1, A1, B1, B2, B2. When you add the counter for single-particle energies, they would be 2A1, 3A1, 4A1, 1B1, 1B2, 2B2, however, because the 1A1 orbital is explicitly doubly occupied and not taken into the active space.

The pseudocanonical orbitals are then printed to a molden file in the snippet https://github.com/SebWouters/CheMPS2/blob/9fd492793c91aeb5c3e04854aca28622abdecdfb/integrals/psi4plugins/dmrg.cc#L994 (lines 994 - 1028). The variables are just the alpha and beta coefficients, single-particle energies, and occupations, respectively.

Does this fix your problems?

Best wishes, Sebastian

[cuanto]

Hi Sebastian,

I solved the problem, I think.

By now I use C1 symmetry to get the 1-RDM and 2-RDM in the same order as in the molden file, I need to change the print loops to take this in to account. It is my first time using c++ and psi4, so it is not trivial for me.

Concerning the orbitals, I put the printer calls before the pseudocanonical rotation, including a new printer for molden, as my code internally construct it. For some other tasks I want to do, I need to recover the HF occupancies inside the plugin, I was reading the psi4 doc about the wavefunction object but I can not find in which var occupancies are stored, how can easily access it ?

As you say I think it was a misunderstanding between the names for the orbitals.

With these changes I could recover the energy in my program !!!! :-)

Thanks for your time and patience.

Best regargs, Jose Luis

[SebWouters]

Hi Jose,

You're welcome to ask questions.

That's good news!

HF occupancies for RHF and ROHF can be obtained with

int * docc       = wfn->doccpi();
int * socc       = wfn->soccpi();

These are arrays which contain per irrep the number of doubly occupied (docc) and singly occupied (socc) orbitals. So for RHF, socc[ i ] == 0 for all i and sum_i docc[ i ] = N_electrons / 2.

Best wishes, Sebastian

[cuanto]

Hi Sebastian,

Thanks for all, by the moment all is working in point group C1, I will read carefully the psi4 API about the irreps and how to operate with it, If i find some trouble I will ask. :-)

Best wishes, Jose Luis