[CP2K-user] [CP2K:18914] molecular polarizability calculated via cp2k is different from psi4 and gaussian09

[jin xi]

Hi Jürg,

Thanks for your answer. I modifield the input file as you suggested, and I 
also reduced the EPS_SCF 1.0E-8 of the &SCF section. Happily, the 
calculation results are consistent with the psi4 and gaussian09 results.

Regards
Jin

[image: xx.png]
在2023年6月11日星期日 UTC+8 00:39:06<Jürg Hutter> 写道：

> Hi
>
> your setup is almost correct. In order to do a fully non-periodic 
> calculation you also
> have to specify a non-periodic electrical field. Instead of
> &PERIODIC_EFIELD you need &EFIELD
>
> For GAPW calculations it is also better to use PRECONDITIONER FULL_ALL.
>
> regards
> JH
>
> ________________________________________
> From: cp... at googlegroups.com <cp... at googlegroups.com> on behalf of jin xi 
> <jinxi.i... at gmail.com>
> Sent: Saturday, June 10, 2023 5:08 AM
> To: cp2k
> Subject: [CP2K:18909] molecular polarizability calculated via cp2k is 
> different from psi4 and gaussian09
>
> 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.
> [1.png]
> [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
> [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 <(365)%20745-0550> 1.0334718700 .0088329700
> H -.3725226950 -.5271066900 .8861702200
> H -.3580328250 -.5066056500 <(506)%20605-6500> -.9039205000 
> <(903)%20920-5000>
>
> 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+uns... at googlegroups.com<mailto:cp2k+uns... 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
> <
> https://groups.google.com/d/msgid/cp2k/780660d3-ae84-4beb-b3c5-6c12afee7d66n%40googlegroups.com?utm_medium=email&utm_source=footer
> >.
>

-- 
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/840acdea-b25a-4606-a7d0-b0efc86d7c22n%40googlegroups.com.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20230610/5ac74ae6/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: xx.png
Type: image/png
Size: 51260 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20230610/5ac74ae6/attachment-0001.png>