[CP2K-user] [CP2K:18909] molecular polarizability calculated via cp2k is different from psi4 and gaussian09
jin xi
jinxi.insilico at gmail.com
Sat Jun 10 03:08:35 UTC 2023
Dear CP2K users,
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.
[image: 1.png]
[image: 2.png]
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.
#####==== cp2k ch4.inp ====#####
&GLOBAL
PROJECT CH4
RUN_TYPE ENERGY
PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep ! Electronic structure method (DFT,...)
&DFT
BASIS_SET_FILE_NAME ./d-aug-cc-pvdz
BASIS_SET_FILE_NAME ./d-aug-cc-pvdz-jkfit
POTENTIAL_FILE_NAME POTENTIAL
RESTART_FILE_NAME CH4-RESTART.wfn
&MGRID
CUTOFF 400
NGRIDS 6
REL_CUTOFF 60
&END MGRID
&POISSON ! Solver requested for non periodic
calculations
PERIODIC NONE
POISSON_SOLVER WAVELET ! Type of solver
&END POISSON
&QS ! Parameters needed to set up the Quickstep
framework
METHOD GAPW ! Method: gaussian and augmented plane
waves
EPS_DEFAULT 1.0E-10
&END QS
&SCF
SCF_GUESS RESTART
EPS_SCF 1.0E-6
MAX_SCF 100
&OUTER_SCF
MAX_SCF 50
EPS_SCF 1.0E-6
&END OUTER_SCF
&OT
PRECONDITIONER FULL_KINETIC
MINIMIZER CG
&END OT
&END SCF
&XC ! Parametes needed to compute the electronic
exchange potential
&XC_FUNCTIONAL
&LYP
SCALE_C 0.81
&END
&BECKE88
SCALE_X 0.72
&END
&VWN
FUNCTIONAL_TYPE VWN3
SCALE_C 0.19
&END
&XALPHA
SCALE_X 0.08
&END
&END XC_FUNCTIONAL
&HF
&SCREENING
EPS_SCHWARZ 1.0E-10
&END
&MEMORY
MAX_MEMORY 100
&END
FRACTION 0.20
&END
&END XC
&PERIODIC_EFIELD
INTENSITY 1.8897261250E-5
POLARISATION -1.0 0.0 0.0 ！change the polarisation for x+ x- y+ y-
z+ z-
&END PERIODIC_EFIELD
&PRINT
&MOMENTS
PERIODIC FALSE
REFERENCE COM
&EACH
MD 1
&END EACH
FILENAME =dipole.txt
&END MOMENTS
&END PRINT
&END DFT
&SUBSYS
&CELL
ABC 10. 10. 10.
PERIODIC NONE ! Non periodic calculations. That's why
the POISSON section is needed
&END CELL
&TOPOLOGY ! Section used to center the atomic
coordinates in the given box. Useful for big molecules
&CENTER_COORDINATES
&END
&END
&COORD
C 0.000003015000 0.000003640000 0.000013990000
H 1.096282485000 0.000218630000 0.008833370000
H -0.365745055000 1.033471870000 0.008832970000
H -0.372522695000 -0.527106690000 0.886170220000
H -0.358032825000 -0.506605650000 -0.903920500000
&END COORD
&KIND H ! potential and basis for
H
BASIS_SET d-aug-cc-pVDZ
#BASIS_SET AUX_FIT d-aug-cc-pVDZ-JKFIT
POTENTIAL ALL
&END KIND
&KIND C ! potential and basis for C
BASIS_SET d-aug-cc-pVDZ
#BASIS_SET AUX_FIT d-aug-cc-pVDZ-JKFIT
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL
[image: 3.png]
#####==== gaussian CH4.gif ====#####
%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 -.9039205000
symmetry 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]))
--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+unsubscribe at googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/780660d3-ae84-4beb-b3c5-6c12afee7d66n%40googlegroups.com.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20230609/ad7aa66e/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 3.png
Type: image/png
Size: 86309 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20230609/ad7aa66e/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 2.png
Type: image/png
Size: 32740 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20230609/ad7aa66e/attachment-0004.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 1.png
Type: image/png
Size: 78987 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20230609/ad7aa66e/attachment-0005.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: input.rar
Type: application/x-rar
Size: 377251 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20230609/ad7aa66e/attachment-0001.bin>
More information about the CP2K-user
mailing list