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 <div><img alt="xx.png" width="1299px" height="450px" src="cid:682d68ea-a32c-4c16-a8e3-bf0a6813a321" /> </div><div class="gmail_quote"><div dir="auto" class="gmail_attr">在2023年6月11日星期日 UTC+8 00:39:06<Jürg Hutter> 写道： </div><blockquote class="gmail_quote" style="margin: 0 0 0 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">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: <a href data-email-masked rel="nofollow">cp...@googlegroups.com</a> <<a href data-email-masked rel="nofollow">cp...@googlegroups.com</a>> on behalf of jin xi <<a href data-email-masked rel="nofollow">jinxi.i...@gmail.com</a>> 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 -.<a href="tel:(365)%20745-0550" value="+13657450550" target="_blank" rel="nofollow">3657450550</a> 1.0334718700 .0088329700 H -.3725226950 -.5271066900 .8861702200 H -.3580328250 -.<a href="tel:(506)%20605-6500" value="+15066056500" target="_blank" rel="nofollow">5066056500</a> -.<a href="tel:(903)%20920-5000" value="+19039205000" target="_blank" rel="nofollow">9039205000</a> 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 <a href data-email-masked rel="nofollow">cp2k+uns...@googlegroups.com</a><mailto:<a href data-email-masked rel="nofollow">cp2k+uns...@googlegroups.com</a>>. To view this discussion on the web visit <a href="https://groups.google.com/d/msgid/cp2k/780660d3-ae84-4beb-b3c5-6c12afee7d66n%40googlegroups.com" target="_blank" rel="nofollow" data-saferedirecturl="https://www.google.com/url?hl=zh-CN&q=https://groups.google.com/d/msgid/cp2k/780660d3-ae84-4beb-b3c5-6c12afee7d66n%2540googlegroups.com&source=gmail&ust=1686539068148000&usg=AOvVaw0CHZm2wg-1YOQZO6fdfDvW">https://groups.google.com/d/msgid/cp2k/780660d3-ae84-4beb-b3c5-6c12afee7d66n%40googlegroups.com</a><<a href="https://groups.google.com/d/msgid/cp2k/780660d3-ae84-4beb-b3c5-6c12afee7d66n%40googlegroups.com?utm_medium=email&utm_source=footer" target="_blank" rel="nofollow" data-saferedirecturl="https://www.google.com/url?hl=zh-CN&q=https://groups.google.com/d/msgid/cp2k/780660d3-ae84-4beb-b3c5-6c12afee7d66n%2540googlegroups.com?utm_medium%3Demail%26utm_source%3Dfooter&source=gmail&ust=1686539068148000&usg=AOvVaw3gv1G1EVNSimn9xAPM6cQm">https://groups.google.com/d/msgid/cp2k/780660d3-ae84-4beb-b3c5-6c12afee7d66n%40googlegroups.com?utm_medium=email&utm_source=footer</a>>. </blockquote></div> -- 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 <a href="mailto:cp2k+unsubscribe@googlegroups.com">cp2k+unsubscribe@googlegroups.com</a>. To view this discussion on the web visit <a href="https://groups.google.com/d/msgid/cp2k/840acdea-b25a-4606-a7d0-b0efc86d7c22n%40googlegroups.com?utm_medium=email&utm_source=footer">https://groups.google.com/d/msgid/cp2k/840acdea-b25a-4606-a7d0-b0efc86d7c22n%40googlegroups.com</a>.