I wish to use the BFD pseudopotentials’ basis along with the ECP for geometry optimization in Psi4. I do not quite understand how to design the input. I created the requisite basis set file with the ECP information as shown here:
Although I am trying to use ECPs for C,H etc.
But that doesn’t seem to work by importing just the basis file.
Kindly let me know if it is possible. I am very new to PSI4 and CC in general. So my question might be very trivial. Unfortunately, I do not know how to proceed.
You shouldn’t need to create your own basis set file, as the psi4 basis library has the ecp info. Was there a particular reason you were constructing your own? What is the input file that you tried? And the output that failed?
The below works fine for me, so without seeing the error you’re getting, I suspect a formatting error for the basis set you composed yourself. Especially make sure the GBS/ECP interface looks like this. Check where the **** lines are and note that in general, comment lines ! stuff can’t appear within a basis set.
molecule mol {
Te
H 1 1.0
H 1 1.0 2 90.0
}
set {
scf_type cd
cc_type cd
#basis bfd-da
basis def2-tzvp
freeze_core true
cholesky_tolerance 1e-9
}
optimize('ccsd(t)')
Thank you for your response. Here is my basis file attached:
!Pseudopotential basis
!References:
!Energy-consistent pseudopotentials for quantum Monte Carlo calculations Burkatzki M, Filippi C, Dolg M, J Chem Phys 2007; 126(23)
!H.ECP.BFD.5s1p.3s1p…
!C.ECP.BFD.11s11p1d.3s3p1d…
!
H 0
S 3 1.00
6.46417546 0.06815644
1.13891461 0.36325459
0.28003249 0.75241006
S 1 1.00
0.05908405 1.00000000
S 1 1.00
0.02974000 1.00000000
P 1 1.00
0.51368060 1.00000000
C 0
S 9 1.00
0.05134400 0.01399100
0.10261900 0.16985200
0.20510000 0.39752900
0.40992400 0.38036900
0.81929700 0.18011300
1.63749400 -0.03351200
3.27279100 -0.12149900
6.54118700 0.01517600
13.07359400 -0.00070500
S 1 1.00
0.12785200 1.00000000
S 1 1.00
0.04690000 1.00000000
P 9 1.00
0.02928100 0.00178700
0.05854700 0.05042600
0.11706300 0.19163400
0.23406400 0.30266700
0.46800300 0.28986800
0.93575700 0.21097900
1.87101600 0.11202400
3.74103500 0.05442500
7.48007600 0.02193100
P 1 1.00
0.14916100 1.00000000
P 1 1.00
0.04041000 1.00000000
D 1 1.00
0.56116000 1.00000000
!ECP
H 0
H-ECP 0 0
s
3
1 25.000000000000 1.000000000000
2 9.368618758833 -8.228005709676
3 10.821821902641 25.000000000000
C 0
C-ECP 1 2
s
3
1 8.359738210000 4.000000000000
2 3.938312580000 -19.175373230000
3 4.483618880000 33.438952850000
p
1
2 5.029916370000 22.551641910000
Yes, I think your basis set file could use some **** separators. Follow the pattern in the def2-tzvp.gbs you originally linked to to reformat your basis file. The parser is forgiving of white space differences within a line, but it does like atom blocks to be separated just so. Don’t forget the “spherical” at the top and all the **** separators. When it reads in correctly, it should look like the below. And you should also check the accounting info printed to make sure it matches expectations.
Name: DEF2-TZVP
Role: ORBITAL
Keyword: BASIS
atoms 1 entry TE line 2156 (ECP: line 3414) file /Users/loriab/linux/psihub/hrw-labfork/objdir-dev27/stage/usr/local/psi4/share/psi4/basis/def2-tzvp.gbs
atoms 2-3 entry H line 15 file /Users/loriab/linux/psihub/hrw-labfork/objdir-dev27/stage/usr/local/psi4/share/psi4/basis/def2-tzvp.gbs
...
==> Primary Basis <==
Basis Set: DEF2-TZVP
Blend: DEF2-TZVP
Number of shells: 24
Number of basis function: 62
Number of Cartesian functions: 71
Spherical Harmonics?: true
Max angular momentum: 3
Core potential: DEF2-TZVP
Number of shells: 4
Number of ECP primitives: 18
Number of ECP core electrons: 28
Max angular momentum: 3
!References:
!Energy-consistent pseudopotentials for quantum Monte Carlo calculations Burkatzki M, Filippi C, Dolg M, J Chem Phys 2007; 126(23)
!H.ECP.BFD.5s1p.3s1p…
!C.ECP.BFD.11s11p1d.3s3p1d…
!
****
H 0
S 3 1.00
6.46417546 0.06815644
1.13891461 0.36325459
0.28003249 0.75241006
S 1 1.00
0.05908405 1.00000000
S 1 1.00
0.02974000 1.00000000
P 1 1.00
0.51368060 1.00000000
****
C 0
S 9 1.00
0.05134400 0.01399100
0.10261900 0.16985200
0.20510000 0.39752900
0.40992400 0.38036900
0.81929700 0.18011300
1.63749400 -0.03351200
3.27279100 -0.12149900
6.54118700 0.01517600
13.07359400 -0.00070500
S 1 1.00
0.12785200 1.00000000
S 1 1.00
0.04690000 1.00000000
P 9 1.00
0.02928100 0.00178700
0.05854700 0.05042600
0.11706300 0.19163400
0.23406400 0.30266700
0.46800300 0.28986800
0.93575700 0.21097900
1.87101600 0.11202400
3.74103500 0.05442500
7.48007600 0.02193100
P 1 1.00
0.14916100 1.00000000
P 1 1.00
0.04041000 1.00000000
D 1 1.00
0.56116000 1.00000000
****
!ECP
H 0
H-ECP 0 0
s
3
1 25.000000000000 1.000000000000
2 9.368618758833 -8.228005709676
3 10.821821902641 25.000000000000
C 0
C-ECP 1 2
s
3
1 8.359738210000 4.000000000000
2 3.938312580000 -19.175373230000
3 4.483618880000 33.438952850000
p
1
2 5.029916370000 22.551641910000
Your basis file worked for me. Is it not finding the basis file (place alongside input file or in library or set PSIPATH)? You’ll have to paste an error, I think.
I do not get any errors. But the energy value is not what I expect. I expect an energy around -26.2 Hartree while I get energy of -65.4 Hartree (for 4 carbon and 6 hydrogen atoms. Here is my input:
memory 38 gb
molecule buta {
units bohr
C -1.1579352682 -0.7336520997 0.0000000000
C 1.1579352682 0.7336520997 0.0000000000
C -3.4664188656 0.2981874225 0.0000000000
C 3.4664188656 -0.2981874225 0.0000000000
H -3.7066332742 2.3264227662 0.0000000000
H 3.7066332742 -2.3264227662 0.0000000000
H -0.9631933359 -2.7721807808 0.0000000000
H 0.9631933359 2.7721807808 0.0000000000
H -5.1469502967 -0.8555507642 0.0000000000
H 5.1469502967 0.8555507642 0.0000000000
}
set {
scf_type cd
cc_type cd
basis bfd-da
freeze_core true
}
I get -26.. And yes, there’s an accounting output for the ECP part (see below SCF calc). Note that you’ll only get the ECP part for the orbital, not the fitting, basis. I also did a ccsd(t) single point and it got -26.685125334.
SCF output:
=> Loading Basis Set <=
Name: MONIKA
Role: ORBITAL
Keyword: BASIS
atoms 1-4 entry C line 23 (ECP: line 65) file /Users/loriab/linux/psihub/hrw-labfork/objdir-dev27/stage/usr/local/psi4/share/psi4/basis/monika.gbs
atoms 5-10 entry H line 11 (ECP: line 58) file /Users/loriab/linux/psihub/hrw-labfork/objdir-dev27/stage/usr/local/psi4/share/psi4/basis/monika.gbs
There are an even number of electrons - assuming singlet.
Specify the multiplicity in the molecule input block.
---------------------------------------------------------
SCF
by Justin Turney, Rob Parrish, Andy Simmonett
and Daniel Smith
RHF Reference
1 Threads, 500 MiB Core
---------------------------------------------------------
==> Geometry <==
Molecular point group: c2h
Full point group: C2h
Geometry (in Bohr), charge = 0, multiplicity = 1:
Center X Y Z Mass
------------ ----------------- ----------------- ----------------- -----------------
C -1.157935268200 -0.733652099700 0.000000000000 12.000000000000
C 1.157935268200 0.733652099700 0.000000000000 12.000000000000
C -3.466418865600 0.298187422500 0.000000000000 12.000000000000
C 3.466418865600 -0.298187422500 0.000000000000 12.000000000000
H -3.706633274200 2.326422766200 0.000000000000 1.007825032070
H 3.706633274200 -2.326422766200 0.000000000000 1.007825032070
H -0.963193335900 -2.772180780800 0.000000000000 1.007825032070
H 0.963193335900 2.772180780800 0.000000000000 1.007825032070
H -5.146950296700 -0.855550764200 0.000000000000 1.007825032070
H 5.146950296700 0.855550764200 0.000000000000 1.007825032070
Running in c2h symmetry.
Rotational constants: A = 1.40755 B = 0.14913 C = 0.13484 [cm^-1]
Rotational constants: A = 42197.41255 B = 4470.69223 C = 4042.41066 [MHz]
Nuclear repulsion = 56.421406026135578
Charge = 0
Multiplicity = 1
Electrons = 22
Nalpha = 11
Nbeta = 11
==> Algorithm <==
SCF Algorithm Type is DF.
DIIS enabled.
MOM disabled.
Fractional occupation disabled.
Guess Type is SAD.
Energy threshold = 1.00e-06
Density threshold = 1.00e-06
Integral threshold = 0.00e+00
==> Primary Basis <==
Basis Set: MONIKA
Blend: MONIKA
Number of shells: 52
Number of basis function: 104
Number of Cartesian functions: 108
Spherical Harmonics?: true
Max angular momentum: 2
Core potential: MONIKA
Number of shells: 14
Number of ECP primitives: 34
Number of ECP core electrons: 8
Max angular momentum: 1
=> Loading Basis Set <=
Name: (MONIKA AUX)
Role: JKFIT
Keyword: DF_BASIS_SCF
atoms 1-4 entry C line 203 file /Users/loriab/linux/psihub/hrw-labfork/objdir-dev27/stage/usr/local/psi4/share/psi4/basis/def2-qzvpp-jkfit.gbs
atoms 5-10 entry H line 23 file /Users/loriab/linux/psihub/hrw-labfork/objdir-dev27/stage/usr/local/psi4/share/psi4/basis/def2-qzvpp-jkfit.gbs
==> Pre-Iterations <==
-------------------------------------------------------
Irrep Nso Nmo Nalpha Nbeta Ndocc Nsocc
-------------------------------------------------------
Ag 39 39 0 0 0 0
Bg 13 13 0 0 0 0
Au 13 13 0 0 0 0
Bu 39 39 0 0 0 0
-------------------------------------------------------
Total 104 104 11 11 11 0
-------------------------------------------------------
==> Integral Setup <==
==> DFJK: Density-Fitted J/K Matrices <==
J tasked: Yes
K tasked: Yes
wK tasked: No
OpenMP threads: 1
Integrals threads: 1
Memory (MB): 375
Algorithm: Core
Integral Cache: NONE
Schwarz Cutoff: 1E-12
Fitting Condition: 1E-12
=> Auxiliary Basis Set <=
Basis Set: (MONIKA AUX)
Blend: DEF2-QZVPP-JKFIT
Number of shells: 136
Number of basis function: 408
Number of Cartesian functions: 476
Spherical Harmonics?: true
Max angular momentum: 4
Minimum eigenvalue in the overlap matrix is 3.3968229837E-05.
Using Symmetric Orthogonalization.
SCF Guess: Superposition of Atomic Densities via on-the-fly atomic UHF.
==> Iterations <==
Total Energy Delta E RMS |[F,P]|
@DF-RHF iter 0: -17.82780749195197 -1.78278e+01 7.96994e-02
Occupation by irrep:
Ag Bg Au Bu
DOCC [ 5, 1, 1, 4 ]
@DF-RHF iter 1: -24.03109731170365 -6.20329e+00 2.47349e-02
@DF-RHF iter 2: -25.25523881675596 -1.22414e+00 2.21317e-02 DIIS
@DF-RHF iter 3: -26.11523154484755 -8.59993e-01 3.61683e-03 DIIS
@DF-RHF iter 4: -26.14679182961096 -3.15603e-02 6.77490e-04 DIIS
@DF-RHF iter 5: -26.14872937169044 -1.93754e-03 2.18405e-04 DIIS
@DF-RHF iter 6: -26.14891784465957 -1.88473e-04 3.66714e-05 DIIS
@DF-RHF iter 7: -26.14892714371793 -9.29906e-06 7.91455e-06 DIIS
@DF-RHF iter 8: -26.14892753753453 -3.93817e-07 1.95824e-06 DIIS
@DF-RHF iter 9: -26.14892756202436 -2.44898e-08 4.08377e-07 DIIS
CCSD(T) sp:
psi4.core.clean()
psi4.set_memory('2 GB')
buta = psi4.geometry("""
units bohr
C -1.1579352682 -0.7336520997 0.0000000000
C 1.1579352682 0.7336520997 0.0000000000
C -3.4664188656 0.2981874225 0.0000000000
C 3.4664188656 -0.2981874225 0.0000000000
H -3.7066332742 2.3264227662 0.0000000000
H 3.7066332742 -2.3264227662 0.0000000000
H -0.9631933359 -2.7721807808 0.0000000000
H 0.9631933359 2.7721807808 0.0000000000
H -5.1469502967 -0.8555507642 0.0000000000
H 5.1469502967 0.8555507642 0.0000000000
""")
# E = psi4.energy('hf/monika')
# print(E)
psi4.set_options({
'scf_type': 'cd',
'cc_type': 'cd',
'basis': 'monika',
'freeze_core': 'true'})
E = psi4.energy('ccsd(t)')
print(E)
I use the same input and I get the following:
==> Algorithm <==
SCF Algorithm Type is CD.
DIIS enabled.
MOM disabled.
Fractional occupation disabled.
Guess Type is SAD.
Energy threshold = 1.00e-08
Density threshold = 1.00e-08
Integral threshold = 0.00e+00
==> Primary Basis <==
Basis Set: BFD-DA
Blend: BFD-DA
Number of shells: 52
Number of basis function: 104
Number of Cartesian functions: 108
Spherical Harmonics?: true
Max angular momentum: 2
==> Pre-Iterations <==
Irrep Nso Nmo Nalpha Nbeta Ndocc Nsocc
A 104 104 0 0 0 0
==> Integral Setup <==
==> CDJK: Cholesky-decomposed J/K Matrices <==
J tasked: Yes
K tasked: Yes
wK tasked: No
OpenMP threads: 1
Integrals threads: 1
Memory (MB): 1430
Algorithm: Core
Integral Cache: SAVE
Schwarz Cutoff: 1E-12
Cholesky tolerance: 1.00E-04
No. Cholesky vectors: 387
Minimum eigenvalue in the overlap matrix is 3.3968229837E-05.
Using Symmetric Orthogonalization.
SCF Guess: Superposition of Atomic Densities via on-the-fly atomic UHF.
==> Iterations <==
Total Energy Delta E RMS |[F,P]|
@RHF iter 0: -67.96525475603923 -6.79653e+01 1.84152e-02 @RHF iter 1: -65.33663418209393 2.62862e+00 3.52024e-03 @RHF iter 2: -65.43133232694314 -9.46981e-02 2.32514e-03 DIIS @RHF iter 3: -65.46103871277732 -2.97064e-02 7.21491e-04 DIIS @RHF iter 4: -65.47038234245221 -9.34363e-03 3.14163e-04 DIIS @RHF iter 5: -65.47238664039669 -2.00430e-03 1.15699e-04 DIIS @RHF iter 6: -65.47278276981467 -3.96129e-04 3.97016e-05 DIIS @RHF iter 7: -65.47282366581705 -4.08960e-05 1.00503e-05 DIIS @RHF iter 8: -65.47282733775010 -3.67193e-06 4.32637e-06 DIIS @RHF iter 9: -65.47282880348095 -1.46573e-06 2.70288e-06 DIIS @RHF iter 10: -65.47282974953623 -9.46055e-07 1.71150e-06 DIIS @RHF iter 11: -65.47283024146479 -4.91929e-07 8.08923e-07 DIIS @RHF iter 12: -65.47283029325962 -5.17948e-08 3.20860e-07 DIIS @RHF iter 13: -65.47283029955256 -6.29294e-09 1.39299e-07 DIIS @RHF iter 14: -65.47283030158739 -2.03484e-09 9.63266e-08 DIIS @RHF iter 15: -65.47283030290248 -1.31509e-09 8.20982e-08 DIIS @RHF iter 16: -65.47283030401977 -1.11729e-09 7.64326e-08 DIIS @RHF iter 17: -65.47283030613318 -2.11341e-09 6.82422e-08 DIIS @RHF iter 18: -65.47283031202079 -5.88761e-09 3.74745e-08 DIIS @RHF iter 19: -65.47283031385726 -1.83647e-09 1.05419e-08 DIIS @RHF iter 20: -65.47283031392502 -6.77574e-11 3.35938e-09 DIIS
Also, when I use the previous input for TiH2 with the def2-tzvp input; I still don’t get the core potential read in; is there some problem with my installation or am I executing it incorrectly?
Wouldn’t hurt to update. psi4 update psi4 -c psi4/label/dev -c psi4 then post your conda list. That should work on Linux. If you’re on Mac I may need to give you an extra line.