Electric field perturbation doesn't work on ccsd energy calculation

#1

Hi,

I recently encountered a problem with adding electric field on a ccsd energy calculation.
The problem is that when I used the following input, the perturbation would be already added in SCF calculation, which is not what I want. (I just need the unrelaxed orbitals from SCF calculation and conduct ccsd energy calculation with perturbation.)

molecule mol {
0 1
C .0000030150 .0000036400 .0000139900
H 1.0962824850 .0002186300 .0088333700
H -.3657450550 1.0334718700 .0088329700
H -.3725226950 -.5271066900 .8861702200
H -.3580328250 -.5066056500 -.9039205000

symmetry c1
no_reorient
no_com
}

set {
basis d-aug-cc-pVDZ
guess sad
scf_type direct
freeze_core true
df_basis_scf d-aug-cc-pvdz-jkfit
e_convergence 5.0E-10
d_convergence 5.0E-10
r_convergence 5.0E-9
}

pert = 3.7794522500e-5
set perturb_h true
set perturb_with dipole
set perturb_dipole [$pert, 0, 0]
energy(‘ccsd’)

However, when I change the input to the following, the calculation didn’t take the perturbation at all, the output was the same as the one without any perturbation setting…


scf_e, scf_wfn = energy(‘scf’, return_wfn=True)
pert = 3.7794522500e-5
set perturb_h true
set perturb_with dipole
set perturb_dipole [$pert, 0, 0]
energy(‘ccsd’,ref_wfn=scf_wfn)

So does anyone know how to add the perturbation properly to a ccsd energy calculation? Thanks a lot!

#2

You can add the perturbation to the fock matrix of the reference wavefunction:

scf_e, scf_wfn = energy('scf', return_wfn=True)
pert = 3.7794522500e-5
mints = psi4.core.MintsHelper(scf_wfn)
dipole = mints.ao_dipole()
dipole[0].scale(pert) # dipole_x perturbation
Perturb_matrix = scf_wfn.Fa()
Perturb_matrix.add(dipole[0])
Perturb_matrix.print_out() # checking
scf_wfn.Fa().print_out() # checking
energy('ccsd',ref_wfn=scf_wfn)

This should work as the CCSD module should grab the AO basis fock matrix from the reference wavefunction.

#3

Just tried it. This works! Thanks a lot.

#4

As a note @andysim and I fixed this in the next Psi4 version (I believe?).

Molecule in uniform electric field
#5

Do you mean 1.2 version? I used Psi4 1.2 version.

#6

Ah apologies, I read through this too quickly. The perturbations are only added at the SCF level and propagated onwards, this is not exactly typical where you want to break that cycle. So nothing really to fix here.

BTW 1.2 is not out yet, 1.1 is the current version.

#7

Hi, I am now getting another question about this procedure. It seems that when I change the cc_type to “df”, the syntax doesn’t work anymore. The error message looks like:

==> 3-index integrals <==

PSIO_ERROR: Can’t find TOC Entry (Q|mn) Integrals
PSIO_ERROR: unit = 97, errval = 13
PSIO_ERROR: 13 (no such TOC entry)
Traceback (most recent call last):
File “/opt/software/psi4/psi4-github-install/bin/psi4”, line 262, in
exec(content)
File “”, line 48, in
File “/opt/software/psi4/psi4-github-install/lib//psi4/driver/driver.py”, line 467, in energy
wfn = procedures[‘energy’][lowername](lowername, molecule=molecule, **kwargs)
File “/opt/software/psi4/psi4-github-install/lib//psi4/driver/procrouting/proc.py”, line 754, in select_ccsd
return func(name, **kwargs)
File “/opt/software/psi4/psi4-github-install/lib//psi4/driver/procrouting/proc.py”, line 3891, in run_fnodfcc
fnocc_wfn = core.fnocc(ref_wfn)

RuntimeError:
Fatal Error: PSIO Error
Error occurred in file: /opt/software/psi4/psi4-github/psi4/src/psi4/libpsio/error.cc on line: 129
The most recent 5 function calls were:

psi::PsiException::PsiException(std::string, char const*, int)
psi::psio_error(unsigned long, unsigned long)
psi::PSIO::read(unsigned long, char const*, char*, unsigned long, psi::psio_address, psi::psio_address*)

Could you help me with that? And also what if I want to use “cc_type=cd”?

Best!

#8

The problem is that the CCSD procedure expects DF_INTS_IO to be set to SAVE. Usually this is done by the driver before the SCF procedure, but since you ran the SCF procedure separately, it was never set. You should be able to fix the issue with

set DF_INTS_IO SAVE

This option should also fix any issues with cc_type cd as well.

Problem when passing a reference wavefunction to DF-CCSD in the' fnocc' module
#9

Hi, thank you for your reply! I tried to add the syntax in my set session, but it seems that for cc_type df calculation, the perturbation was not added to the reference wavefunction. The result was the same as the one without perturbation. And for cc_type cd calculation, Psi4 gave some segmentation error when it got to the 3-index integral part… Is there anything I set wrong? Could you help me with this? Thanks a lot!

set {
basis d-aug-cc-pVDZ
guess sad
scf_type df
cc_type df
freeze_core true
df_basis_cc aug-cc-pvqz-ri
df_basis_scf aug-cc-pvqz-jkfit
e_convergence 5.0E-10
d_convergence 5.0E-10
r_convergence 5.0E-9
DF_INTS_IO SAVE
}

scf_e, scf_wfn = energy(‘scf’, return_wfn=True)
pert = 1.8897261250e-5
pert2 = 1.8897261250e-5
mints = psi4.core.MintsHelper(scf_wfn)
dipole = mints.ao_dipole()
dipole[0].scale(pert) # dipole_x perturbation
dipole[1].scale(pert2) # dipole_x perturbation
Perturb_matrix = scf_wfn.Fa()
Perturb_matrix.add(dipole[0])
Perturb_matrix.add(dipole[1])
energy(‘ccsd’,ref_wfn=scf_wfn)
energy = psi4.get_variable(‘CCSD TOTAL ENERGY’)
print_out("\nYY - CC TOTAL ENERGY [a.u.] {:25.15f}\n".format(energy))

set {
basis d-aug-cc-pVDZ
guess sad
scf_type cd
cc_type cd
freeze_core true
e_convergence 5.0E-10
d_convergence 5.0E-10
r_convergence 5.0E-9
CHOLESKY_TOLERANCE 1.0E-5
DF_INTS_IO SAVE
}

scf_e, scf_wfn = energy(‘scf’, return_wfn=True)
pert = -3.7794522500e-5
mints = psi4.core.MintsHelper(scf_wfn)
dipole = mints.ao_dipole()
dipole[0].scale(pert) # dipole_x perturbation
Perturb_matrix = scf_wfn.Fa()
Perturb_matrix.add(dipole[0])
energy(‘ccsd’,ref_wfn=scf_wfn)
energy = psi4.get_variable(‘CCSD TOTAL ENERGY’)
print_out("\nYY - CC TOTAL ENERGY [a.u.] {:25.15f}\n".format(energy))

#10

In the df-ccsd procedure, the fock matrix is not read from the reference wavefunction. Instead, it is built from the core hamiltonian (in every iteration) which is again built from scratch. So, this would require a change in the source code where the core Hamiltonian (H) is read from the reference wavefunction. This way, you can change H after the scf and before the df-ccsd procedure just like you modified the fock matrix before. So, inside T1Fock() function in df_t1_transformation.cc file in fnocc module, maybe you can try replacing

H = mints->so_kinetic();
H->add(mints->so_potential());

by

H = reference_wavefunction->H();

Then add the perturbation to H after the scf from the python side. @deprince can verify if this approach would work.

#11

I tested this, and I think that change would work, but note the following:

  1. the lines to be changed are in psi4/src/fnocc/df_ccsd.cc (lines 709-710).

  2. On the python side, as @ashutosh said, be sure to add the perturbation to scf_wfn.H, not scf_wfn.Fa, because the DF-CCSD code in fnocc wants to build the two-electron part of the Fock matrix itself.

#12

Hi, Thank you for your reply! I looked into the source code. The only thing I found in T1Fock() function related to the Hamiltonian is " double ** hp = H->pointer();". I tried to put “H = reference_wavefunction->H();” in front of it, but got this error message during compiling:

/home/urs/psi4/psi4/src/psi4/fnocc/df_t1_transformation.cc:130:9: error: invalid use of member function ‘std::shared_ptrpsi::Wavefunction psi::Wavefunction::reference_wavefunction() const’ (did you forget the ‘()’ ?)
H = reference_wavefunction->H();
^~~~~~~~~~~~~~~~~~~~~~
/home/urs/psi4/psi4/src/psi4/fnocc/df_t1_transformation.cc:130:31: error: base operand of ‘->’ is not a pointer
H = reference_wavefunction->H();

#13

Don’t change anything in df_t1_transformation.cc. Change lines 709-710 of df_ccsd.cc. Also, reference_wavefunction_ should have an underscore after it.

#14

I mistakenly posted the wrong file that needs to get changed. Please follow @deprince’s suggestions.

#15

Thanks for your reply!

#16

Thanks for your reply! I changed the source code and recompiled, but the calculation gave a wrong energy compared to the one calculated with CONV ccsd. Is there anything wrong with my input?

scf_e, scf_wfn = energy(‘scf’, return_wfn=True)
pert = -1.8897261250e-5
pert2 = -1.8897261250e-5
mints = psi4.core.MintsHelper(scf_wfn)
dipole = mints.ao_dipole()
dipole[0].scale(pert) # dipole_x perturbation
dipole[1].scale(pert2) # dipole_x perturbation
Perturb_hamiltonian = scf_wfn.H()
Perturb_hamiltonian.add(dipole[0])
Perturb_hamiltonian.add(dipole[1])
energy(‘ccsd’,ref_wfn=scf_wfn)
energy = psi4.get_variable(‘CCSD TOTAL ENERGY’)
print_out("\nYY - CC TOTAL ENERGY [a.u.] {:25.15f}\n".format(energy))

And I think that @ashutosh mentioned that the core hamiltonian will be built from scratch in every iteration. Will they be the same for every iteration? If not, does that mean it’s wrong for me to put the same “reference_wavefunction_->H()” with perturbation for each iteration?

By the way, the “cc_type cd” is still not working, which dies after the “DFCoupledCluster::WriteBanner()” function…

Best