Transition density matrix at CASPT2 level of theory

[Jyotirmoyray]

Sir, In a previous issue I have asked about the transition density matrix at CASSCF level of theory. Can you tell me further whether we can print transition density matrix at CASPT2 level of theory? One more question I have. Is it possible to print oscillator strength for a transition from a lower electronic state to a higher electronic state without performing NACME calculation? If this needs few line code modifications, I will be able to do.

Thank you in advance.

[Jyotirmoyray]

Hello sir, I am still waiting for your reply. More precisely I want to know, can we print transition density matrix at CASPT2 level of theory before symmetrization of the transition density matrix?

[shiozaki]

The answer to the second question is no. Virtually nothing can be gained even if you do that.

The first part is https://github.com/qsimulate-open/bagel/blob/master/src/smith/caspt2grad.cc#L132

Just insert

dtotao->print("AO");

for AO density matrix and

(*d0ms + *(task_->d11()) + *(task_->d1())).print("MO");

for MO density matrix. I will close this issue.

[felixplasser]

Hi, could I just briefly follow up on this. The reason why we are interested in the full non-symmetric 1TDM is because that is what we need in the wavefunction analysis with TheoDORE. Physically speaking, the 1TDM is non-symmetric because the excitation hole and excited electron are in different orbitals. By symmetrising we mix electron and hole. So, for example, when you have a charge transfer state, you can no longer distinguish in the analysis where the excitation comes from and where it goes to if you only have the symmetrised 1TDM.

Just to give a bit of context here. But it's probably still not worth the hassle ...

-Felix

[shiozaki]

Just to clarify - are you saying that these quantities are symmetrized when you print? I'd need to read the code again if so, but it's hard to imagine why these things are symmetrized (dipole integrals are symmetric so only the symmetric part contributes to the transition dipole; but these are calculated for gradient evaluation, and it's hard to imagine I had symmerized for no reason...

[Jyotirmoyray]

Actually I found that 1 particle TDM is symmetric in MO basis but in AO basis 1TDM is not symmetric. (d0ms + (task->d11()) + *(task->d1())).print("MO"); In the above line d0ms is not symmetrized but d1 and d11 are symmetrized in code. In another section of the code all three quantities are symmetrized. https://github.com/qsimulate-open/bagel/blob/a9c1ea15c3e0459c590ede35103dc74fc210fa65/src/smith/caspt2grad.cc#L309

[shiozaki]

You are looking at a place that is different from what I wrote above. The following code is self contained, and I don't think there is any symmetrization. https://github.com/qsimulate-open/bagel/blob/master/src/smith/caspt2grad.cc#L132

\  for (int istate = 1; istate != nstate; ++istate) {
    for (int jstate = 0; jstate != istate; ++jstate) {
      task_->compute_gradient(istate, jstate, make_shared<NacmType>("interstate"), /*nocider=*/true);

      {
        auto d0ms = make_shared<Matrix>(nmobasis, nmobasis);
        if (nact)
          d0ms->add_block(1.0, nclosed, nclosed, nact, nact, task_->d10ms());
        {
          const string dmlabel = "CASPT2 unrelaxed dipole moment: " + to_string(istate) + " - " + to_string(jstate);
          auto dtotao = make_shared<Matrix>(*(task_->smith()->algo()->coeff()) * (*d0ms + *(task_->d11()) + *(task_->d1())) ^ *(task_->smith()->algo()->coeff()));
//          dtotao->print("AO");
//          (*d0ms + *(task_->d11()) + *(task_->d1())).print("MO");
          Dipole dipole(geom_, dtotao, dmlabel);
          auto moment = dipole.compute();
          transition_dipole.push_back(moment);
        }
      }
    }
  }

[shiozaki]

This only loops over half of the i-j combination, but you can adjust to loop to go over all possible pairs.

[Jyotirmoyray]

https://github.com/qsimulate-open/bagel/blob/master/src/smith/caspt2grad.cc#L132 If am not wrong then in the above line "compute_gradient" function calculates 1st order 1tdm (d1) and 2nd order 1tdm(d11). And the function "compute_gradient" is defined on the following line. https://github.com/qsimulate-open/bagel/blob/a9c1ea15c3e0459c590ede35103dc74fc210fa65/src/smith/caspt2grad.cc#L191 Then in the following line d1 and d11 are symmetrized. https://github.com/qsimulate-open/bagel/blob/a9c1ea15c3e0459c590ede35103dc74fc210fa65/src/smith/caspt2grad.cc#L165

So if I am not wrong then https://github.com/qsimulate-open/bagel/blob/master/src/smith/caspt2grad.cc#L132 here symmetrized d1 and d11 are used. Yes d0ms is not symmetrized here.

[Jyotirmoyray]

sorry, https://github.com/qsimulate-open/bagel/blob/a9c1ea15c3e0459c590ede35103dc74fc210fa65/src/smith/caspt2grad.cc#L165 here "compute_gradient" is defined and https://github.com/qsimulate-open/bagel/blob/a9c1ea15c3e0459c590ede35103dc74fc210fa65/src/smith/caspt2grad.cc#L191 here d1 and d11 are symmetrized.

[shiozaki]

No, your comment is incorrect. The line before symmetrization is copy-constructing the quantity and does not affect that code I am talking about.

[Jyotirmoyray]

Sir, if I comment https://github.com/qsimulate-open/bagel/blob/a9c1ea15c3e0459c590ede35103dc74fc210fa65/src/smith/caspt2grad.cc#L191 and https://github.com/qsimulate-open/bagel/blob/a9c1ea15c3e0459c590ede35103dc74fc210fa65/src/smith/caspt2grad.cc#L192 lines and I recompile the code I get different 1 particle TDM in MO basis. I am sending you the input file and the output files. bagel_input.zip bagel_input.zip contains input file and two output files. bagel_tdm_symm.out is obtained without commenting those lines and bagel_tdm_unsymm.out is obtained after commenting those lines. In bagel_tdm_symm.out file the tdm is symmetric and in bagel_tdm_unsymm.out the tdm is not symmetric. Please have a look.

[shiozaki]

Sorry, but I cannot do user support (I am no longer in academia...).

[Jyotirmoyray]

No, your comment is incorrect. The line before symmetrization is copy-constructing the quantity and does not affect that code I am talking about.

I wanted to say that the symmetrization is affecting the code you are talking about.

[felixplasser]

I guess the way forward is to try connecting this to TheoDORE. This might already do what you want, @Jyotirmoyray .

To check things manually, set up a small molecule (H2 or formaldehyde) with a minimal basis set. Then print the 1TDM. In the case of, e.g., the HOMO-LUMO transition only one 1TDM element should strongly deviate from zero - the one whose indices correspond to HOMO and LUMO. So you'd have D_HL = 1 and D_LH = 0. If this is the case, then it works properly in the way we need it.

[shiozaki]

@jwpk1201 Hi Jae Woo - sorry to bother you on this but I started to be not sure if the TRDM code is correct. I am sure you've done careful coding back in the days so could you enlighten me and others on this?

The symmetrization of the d1tmp is not necessary but that of d11tmp comes from the Hylleraas functional <I(1)|E{rs}|J(0)> + <I(0)|E{rs}|J(1)>. For a single state gradient, for which this code was originally written, symmetrization would account for this, but for transition density matrices, I don't think this would work as one wishes because you get <I(1)|E{rs}|J(0)> + <J(0)|E{rs}|I(1)>.

So the question is if smith computes <I(1)|E{rs}|J(0)> + <J(1)|E{rs}|I(0)> directly? Do you remember something?

[jwpk1201]

Hi Toru, I remember that SMITH computes <I(1)|E{rs}|J(0)> + <J(1)|E{rs}|I(0)>. If you want to compute non-symmetric 1TRDM, you will have to carefully modify src/smith/caspt2/MSCASPT2.cc line 159 where you symmetrize 1TRDM to get <I(1)|E{rs}|J(0)> + <J(1)|E{rs}|I(0)>.

[shiozaki]

Makes sense, good to know! (I should have read the code when we wrote a paper back then, but oh well)