SebWouters / CheMPS2

CheMPS2: a spin-adapted implementation of DMRG for ab initio quantum chemistry
GNU General Public License v2.0
70 stars 34 forks source link

Questions about output files #22

Closed nesquik91 closed 8 years ago

nesquik91 commented 8 years ago

Dear CheMPS2 member,

Could I make each alpha and beta orbital information about molecule with CheMPS2 DMRG(SCF,CI,CASPT2)?

Also, I wonder how to obtain molecular total density (alpha + beta density) and spin density (alpha-beta) information with CheMPS2.

Thank you for the help.

SebWouters commented 8 years ago

Hi nesquick91,

Do you mean: Can I get DMRG-SCF orbital coefficients out with respect to the atomic orbital basis?

In case you are using psi4 the answer should be very simple: Call the molden writer after calling energy('dmrg'). @loriab or @jturney : Can you verify that this dumps out the latest wfn->Ca() instead of the Hartree-Fock orbitals?

The DMRG-SCF molecular total density with respect to the molecular orbitals can be obtained by calling

  theDMRG->calc2DMandCorrelations();

and then

  theDMRG->get2DM()->get1RDM_HAM(const int orb1, const int orb2);

The DMRG-CASPT2 molecular total density is not yet implemented. The spin density can not be obtained in CheMPS2, because the code works directly in the averaged spin-ensemble. So density matrices for a spin S wavefunction are in fact

  Gamma = \frac{1}{2S+1} \sum_{S^z} Gamma(S^z)

Spin densities are therefore always zero.

Would it help you if the dmrg.cc plugin allows to print the DMRG-SCF total density with respect to the atomic orbital basis? @loriab or @jturney : Everything in the dmrg.cc plugin is with respect to the SO basis currently. How do I obtain the AO->SO transfo?

Best wishes, Sebastian

nesquik91 commented 8 years ago

Dear Sebastian,

Thank you for your kind reply.

I could obtain .molden format file with "set molden_write true"

here is my "set molden_write true" location with example DMRG-CI input file (http://sebwouters.github.io/CheMPS2/interfaces.html#psi4-dmrg-plugin)

molecule H2O { 0 1 O 0.000000000 0.00 0.000000000 H 0.790689766 0.00 0.612217330 H -0.790689766 0.00 0.612217330 units angstrom }

set molden_write true plugin_load("dmrgscf/dmrgscf.so")

set basis cc-pVDZ set reference rhf set e_convergence 1e-13 set d_convergence 1e-13 set ints_tolerance 0.0

set dmrgscf wfn_irrep 0 set dmrgscf wfn_multp 1 set dmrgscf frozen_docc [ 1 , 0 , 0 , 0 ] set dmrgscf active [ 5 , 0 , 4 , 2 ]

set dmrgscf dmrg_states [ 200, 500, 1000, 1000 ] set dmrgscf dmrg_econv [ 1e-8, 1e-8, 1e-8, 1e-8 ] set dmrgscf dmrg_maxsweeps [ 5, 5, 5, 100 ] set dmrgscf dmrg_noiseprefactors [ 0.03, 0.03, 0.03, 0.0 ] set dmrgscf dmrg_print_corr false set dmrgscf mps_chkpt false

set dmrgscf dmrgscf_max_iter 1

scf()

plugin("dmrgscf.so")

By the way is it possible to use molden_write option with DMRG-CI or DMRG-CASPT2 method? Because psi4 site (http://www.psicode.org/psi4manual/master/autodoc_glossary_options_c.html#term-molden-write-scf) has only mentioned only for specific case.

Also, is FCIDUMP format file is necessary for DMRG-* calculation? Is that so, I wonder how to set input file to consider the file within calculation.

Thank you again.

SebWouters commented 8 years ago

Hi nesquick91,

@loriab @jturney : how does one invoke the moldenwriter at an instance of the wfn object?

If you use psi4 do to DMRG calculations (by calling the dmrg or dmrg-scf plugin) you don't need to generate any fcidump files.

@nesquik91 : in which format would you like to obtain the total density matrix (with respect to AO, MO, ...)?

Best wishes, Sebastian

kannon92 commented 8 years ago

Hey Sebastian,

I do not know how the new psi4 changes will affect the MoldenWrite.

We use the class MoldenWriter, located in libmints/writer.cc. If you have a wavefunction, I believe you can just pass it to the MoldenWriter.

Here is some example code where we call this:

    boost::shared_ptr<MoldenWriter> molden(new      MoldenWriter(Process::environment.wavefunction()));
    std::string filename = get_writer_file_prefix() + ".molden";

    if(remove(filename.c_str()) == 0){
        outfile->Printf("\n  Remove previous molden file named %s.", filename.c_str());
    }
    outfile->Printf("\n  Write molden file to %s.", filename.c_str());
    molden->write(filename, Ca, Ca, diag_F, diag_F, occupation, occupation);
nesquik91 commented 8 years ago

Dear Sebastian,

Actually, my interest is obtain each alpha and beta orbital coefficients and electron densities with DMRG.

So, I guess total density matrix of AO can be a starter..

loriab commented 8 years ago

Should just be:

E, wfn = energy(mtd, return_wfn=True)
molden('moldenfile', wfn)
nesquik91 commented 8 years ago

Dear Sebastian, I've re-installed dmrg plugin with your update dmrg.cc file. But while using

set dmrg dmrg_molden set dmrg dmrg_density_ao

only alpha molecular orbital is recorded in the molden format file.

Could you check this please?

here is my input file

memory 2000 Mb molecule Fe { 0 5 symmetry c1 Fe 0.0 0.0 0.0 units au }

sys.path.insert(0, './..') import dmrg

set basis def2-tzvpp set reference uhf set scf_type PK

set dmrg wfn_irrep 0 set dmrg wfn_multp 5 set dmrg frozen_docc [ 9 ] set dmrg active [ 9 ]

set dmrg dmrg_states [ 2500, 3000, 3000 ] set dmrg dmrg_e_convergence [ 1e-10, 1e-10, 1e-10 ] set dmrg dmrg_maxsweeps [ 5, 5, 10 ] set dmrg dmrg_noiseprefactors [ 0.05, 0.05, 0.0 ] set dmrg dmrg_print_corr false set dmrg dmrg_chkpt false

set dmrg d_convergence 1e-10 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 150 set dmrg dmrg_which_root 1 set dmrg dmrg_state_avg true set dmrg dmrg_active_space no # INPUT; NO; LOC set dmrg dmrg_loc_random false

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

energy('dmrg')

SebWouters commented 8 years ago

Dear nesquick91,

In CheMPS2 alpha and beta orbitals are always equal. So what I call in dmrg.cc (line 1004) is

https://github.com/SebWouters/CheMPS2/blob/d6d416c26cdaac3eedb3eadb3cdb2c1cf864a572/integrals/psi4plugins/dmrg.cc#L1004
molden->write( filename, wfn->Ca(), wfn->Ca(), sp_energies, sp_energies, occupation, occupation );

So if psi4 only prints the alpha orbital coefficients, it is because psi4 detects that the first and second matrices ( wfn->Ca() ) are equal and prints only one of them.

@loriab: Can you confirm?

Best wishes, Sebastian

loriab commented 8 years ago

Looks to be so (https://github.com/psi4/psi4/blob/master/src/lib/libmints/writer.cc#L393) if you follow the SameOcc variable.

nesquik91 commented 8 years ago

Dear Sebastian,

Then..... is it possible to calculate open-shell system via chemps2?? I've tried to calculate Fe atom as you can see it in my input file..

SebWouters commented 8 years ago

Dear nequick91,

It surely is possible.

CheMPS2 works in the spin-ensemble. If you study a wavefunction with spin S (e.g. triplet) then CheMPS2 calculates all 2S+1 (e.g. three) wavefunctions in the multiplet together. It does so because it implicitly factorizes MPS tensors into reduced MPS tensors and Clebsch-Gordan coefficients. It only calculates the reduced MPS tensors and never takes the CG's into account. The energies you so obtain are exactly the same as if you would use two alpha electrons more compared to the number of beta electrons (assuming no quintet etc. is found instead). For the density matrices, I return the spin-summed reduced density matrices, i.e.

Gamma = sum_{Sz} Gamma(Sz} / (2S+1).

They can be easily obtained without ever having to generate a specific wavefunction in the multiplet. Does this help you understand what's going on?

Best wishes, Sebastian

SebWouters commented 8 years ago

@nesquik91

E.g. for the 2-RDM:

Gamma_{ijkl} = 1 / ( 2 S + 1 ) * sum_{sigma, tau, sz} < Psi(sz) | a^+_{i,sigma} a^+_{j,tau} a_{l,tau} a_{k,sigma} | Psi(sz) > 

Best wishes, Sebastian

nesquik91 commented 8 years ago

Dear Sebastian,

Thank you for the kind and detailed explanation. That's really helpful to me.

Thank you again for the many help.

SebWouters commented 8 years ago

@nesquik91

I worked around the problem of spin densities and defined the spin-density for spin-ensembles as follows:

Gamma^{spin}ij = 3 / ( (S+1)(2S+1) ) sum{Sz} Sz < S Sz | a^+{i,up} a{j,up} - a^+{i,down} a{j,down} | S Sz >.

It corresponds exactly for spin singlet, doublet, and triplet to the spin-projection-symmetry-based version where there are 0, 1, or 2 up electrons more than down electrons. For higher spin values, there will be a difference between the spin-projection spin densities and the spin-ensemble spin densities, but the normalization is chosen so that trace(Gamma^{spin}) = 2S.

https://github.com/SebWouters/CheMPS2/commit/e8b0c4a4a8b003d0caa753f52028c7a37c219b84

When I see a chance, I will update the dmrg.cc plugin accordingly.

nesquik91 commented 8 years ago

Dear Sebastian,

Sounds great! Thank you for your effort.

I will look foward to your update.