I am trying to calculate molecular polarizability of a single CH4 molecular via cp2k software using the B3LYP functional and the d-aug-cc-pVDZ basis set. The strategy is as follows. By changing the polarization direction of the electric field, the corresponding molecular dipole moment is obtained, and then the molecular polarizability is calculated by differential method. To verify the correctness of the calculations, I also used PSI4 1.1 and Gaussian09 software to calculate the molecular polarizability at the same functional and basis set, and the results are as follows.
It is not difficult to see that the main diagonal results of the molecular polarizability calculated by cp2k are inconsistent with the calculations of psi4 1.1 and Gaussian09. I tried changing parameters like &MGRID, EPS_default and &XC, but none worked. The input files for each software are as follows. Does anyone know what the reason is?
Thanks in advanced.
%chk=C:\Users\hp\Desktop\spectra\gaussian_methanol\demo\demo.chk
#p b3lyp/gen nosymm polar
demo
0 1
C 0.00000301 0.00000364 0.00001399
H 1.09628248 0.00021863 0.00883337
H -0.36574505 1.03347187 0.00883297
H -0.37252270 -0.52710669 0.88617022
H -0.35803282 -0.50660565 -0.90392050
H 0
S 4 1.00
13.0100000 0.0196850
1.9620000 0.1379770
0.4446000 0.4781480
0.1220000 0.5012400
S 1 1.00
0.1220000 1.0000000
S 1 1.00
0.0297400 1.0000000
S 1 1.00
0.00725 1.000000
P 1 1.00
0.7270000 1.0000000
P 1 1.00
0.1410000 1.0000000
P 1 1.00
0.02730 1.000000
****
C 0
S 9 1.00
6665.0000000 0.0006920
1000.0000000 0.0053290
228.0000000 0.0270770
64.7100000 0.1017180
21.0600000 0.2747400
7.4950000 0.4485640
2.7970000 0.2850740
0.5215000 0.0152040
0.1596000 -0.0031910
S 9 1.00
6665.0000000 -0.0001460
1000.0000000 -0.0011540
228.0000000 -0.0057250
64.7100000 -0.0233120
21.0600000 -0.0639550
7.4950000 -0.1499810
2.7970000 -0.1272620
0.5215000 0.5445290
0.1596000 0.5804960
S 1 1.00
0.1596000 1.0000000
S 1 1.00
0.0469000 1.0000000
S 1 1.00
0.0138 1.000000
P 4 1.00
9.4390000 0.0381090
2.0020000 0.2094800
0.5456000 0.5085570
0.1517000 0.4688420
P 1 1.00
0.1517000 1.0000000
P 1 1.00
0.0404100 1.0000000
P 1 1.00
0.0108 1.000000
D 1 1.00
0.5500000 1.0000000
D 1 1.00
0.1510000 1.0000000
D 1 1.00
0.0415 1.000000
****
#####==== psi4 CH4.in ====#####
memory 12000 mb
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 -.
9039205000symmetry c1
no_reorient
no_com
}
set {
basis d-aug-cc-pVDZ
guess sad
scf_type direct
df_basis_scf d-aug-cc-pvdz-jkfit
e_convergence 1.0E-10
d_convergence 1.0E-10
}
import math
dipole_NE = []
quadrupole_NE = []
e, wfn = energy('B3LYP',return_wfn=True)
dpx = psi4.get_variable('SCF DIPOLE X')
dpy = psi4.get_variable('SCF DIPOLE Y')
dpz = psi4.get_variable('SCF DIPOLE Z')
dptot = math.sqrt(dpx*dpx + dpy*dpy + dpz*dpz)
dipole_NE.append(dpx)
dipole_NE.append(dpy)
dipole_NE.append(dpz)
dipole_NE.append(dptot)
oeprop(wfn,'QUADRUPOLE',title="B3LYP")
qpxx = psi4.get_variable('B3LYP QUADRUPOLE XX')
qpxy = psi4.get_variable('B3LYP QUADRUPOLE XY')
qpyy = psi4.get_variable('B3LYP QUADRUPOLE YY')
qpxz = psi4.get_variable('B3LYP QUADRUPOLE XZ')
qpyz = psi4.get_variable('B3LYP QUADRUPOLE YZ')
qpzz = psi4.get_variable('B3LYP QUADRUPOLE ZZ')
quadrupole_NE.append(qpxx)
quadrupole_NE.append(qpxy)
quadrupole_NE.append(qpyy)
quadrupole_NE.append(qpxz)
quadrupole_NE.append(qpyz)
quadrupole_NE.append(qpzz)
#--------------------------------------------------------------------
mints = psi4.core.MintsHelper(wfn.basisset())
ao_overlap = mints.ao_overlap()
Co = wfn.Ca()
import numpy
numpy.set_printoptions(threshold=numpy.nan)
Da_so = wfn.Da()
np.save('molecule0001_ao_overlap_E0.npy',ao_overlap)
np.save('molecule0001_so_density_E0.npy',Da_so)
np.save('molecule0001_mo_coefficients_E0.npy',Co)
#--------------------------------------------------------------------
#perturbation
pert = 1.8897261250e-5
lambdas = [pert, -pert]
dipole_EX = []
dipole_EY = []
dipole_EZ = []
set perturb_h true
set perturb_with dipole
sign = 'P'
for l in lambdas:
set perturb_dipole [$l, 0, 0]
set scf guess read
e, wfn = energy('B3LYP',return_wfn=True)
dpx = psi4.get_variable('SCF DIPOLE X')
dpy = psi4.get_variable('SCF DIPOLE Y')
dpz = psi4.get_variable('SCF DIPOLE Z')
dipole_EX.append(dpx)
dipole_EX.append(dpy)
dipole_EX.append(dpz)
#--------------------------------------------------------------------
mints = psi4.core.MintsHelper(wfn.basisset())
ao_overlap = mints.ao_overlap()
Co = wfn.Ca()
Da_so = wfn.Da()
Da_mo = Matrix.triplet(wfn.Ca(), Da_so, wfn.Ca(), True, False, False)
np.save('molecule0001_ao_overlap_E' + sign + 'X.npy',ao_overlap)
np.save('molecule0001_so_density_E' + sign + 'X.npy',Da_so)
np.save('molecule0001_mo_coefficients_E' + sign + 'X.npy',Co)
#--------------------------------------------------------------------
#--------------------------------------------------------------------
set perturb_dipole [0, $l, 0]
set scf guess read
e, wfn = energy('B3LYP',return_wfn=True)
oeprop(wfn,'DIPOLE')
dpx = psi4.get_variable('SCF DIPOLE X')
dpy = psi4.get_variable('SCF DIPOLE Y')
dpz = psi4.get_variable('SCF DIPOLE Z')
dipole_EY.append(dpx)
dipole_EY.append(dpy)
dipole_EY.append(dpz)
#--------------------------------------------------------------------
mints = psi4.core.MintsHelper(wfn.basisset())
ao_overlap = mints.ao_overlap()
Co = wfn.Ca()
Da_so = wfn.Da()
np.save('molecule0001_ao_overlap_E' + sign + 'Y.npy',ao_overlap)
np.save('molecule0001_so_density_E' + sign + 'Y.npy',Da_so)
np.save('molecule0001_mo_coefficients_E' + sign + 'Y.npy',Co)
#--------------------------------------------------------------------
#--------------------------------------------------------------------
set perturb_dipole [0, 0, $l]
set scf guess read
e, wfn = energy('B3LYP',return_wfn=True)
oeprop(wfn,'DIPOLE')
dpx = psi4.get_variable('SCF DIPOLE X')
dpy = psi4.get_variable('SCF DIPOLE Y')
dpz = psi4.get_variable('SCF DIPOLE Z')
dipole_EZ.append(dpx)
dipole_EZ.append(dpy)
dipole_EZ.append(dpz)
#--------------------------------------------------------------------
mints = psi4.core.MintsHelper(wfn.basisset())
ao_overlap = mints.ao_overlap()
Co = wfn.Ca()
Da_so = wfn.Da()
np.save('molecule0001_ao_overlap_E' + sign + 'Z.npy',ao_overlap)
np.save('molecule0001_so_density_E' + sign + 'Z.npy',Da_so)
np.save('molecule0001_mo_coefficients_E' + sign + 'Z.npy',Co)
sign = 'N'
#--------------------------------------------------------------------
#--------------------------------------------------------------------
pl_11 = -(dipole_EX[0] - dipole_EX[3])/(2.0*pert*psi_dipmom_au2debye)
pl_12 = -(dipole_EX[1] - dipole_EX[4])/(2.0*pert*psi_dipmom_au2debye)
pl_13 = -(dipole_EX[2] - dipole_EX[5])/(2.0*pert*psi_dipmom_au2debye)
pl_21 = -(dipole_EY[0] - dipole_EY[3])/(2.0*pert*psi_dipmom_au2debye)
pl_22 = -(dipole_EY[1] - dipole_EY[4])/(2.0*pert*psi_dipmom_au2debye)
pl_23 = -(dipole_EY[2] - dipole_EY[5])/(2.0*pert*psi_dipmom_au2debye)
pl_31 = -(dipole_EZ[0] - dipole_EZ[3])/(2.0*pert*psi_dipmom_au2debye)
pl_32 = -(dipole_EZ[1] - dipole_EZ[4])/(2.0*pert*psi_dipmom_au2debye)
pl_33 = -(dipole_EZ[2] - dipole_EZ[5])/(2.0*pert*psi_dipmom_au2debye)
print_out("\nYY-Dipole Moment [Debye] x+ x- y+ y- z+ z-")
print_out("\n")
psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EX[0],dipole_EX[1],dipole_EX[2]))
psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EX[3],dipole_EX[4],dipole_EX[5]))
psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EY[0],dipole_EY[1],dipole_EY[2]))
psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EY[3],dipole_EY[4],dipole_EY[5]))
psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EZ[0],dipole_EZ[1],dipole_EZ[2]))
psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EZ[3],dipole_EZ[4],dipole_EZ[5]))
print_out("\nYY-Polarizability tensor [a.u.]")
print_out("\n")
psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(pl_11,pl_12,pl_13))
psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(pl_21,pl_22,pl_23))
psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(pl_31,pl_32,pl_33))
print_out("\nYY-Dipole Moment(Debye)")
print_out("\n")
psi4.print_out("{}{:24.15f}{}{:22.15f}{}{:22.15f}\n".format("X",dipole_NE[0],"
Y",dipole_NE[1],"
Z",dipole_NE[2]))
psi4.print_out("{}{:22.15f}\n".format("Tol",dipole_NE[3]))
print_out("\nYY-Quadrupole Moments (Debye-Ang)")
print_out("\n")
psi4.print_out("{}{:22.15f}{}{:22.15f}{}{:22.15f}\n".format("XX",quadrupole_NE[0],"
XY",quadrupole_NE[1],"
YY",quadrupole_NE[2]))
psi4.print_out("{}{:22.15f}{}{:22.15f}{}{:22.15f}\n".format("XZ",quadrupole_NE[3],"
YZ",quadrupole_NE[4],"
ZZ",quadrupole_NE[5]))