Dear CP2K users,<br /><br />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.<br /><img alt="1.png" width="534px" height="299px" src="cid:55f3a647-f51d-4392-81b5-17550d919047" /><div><img alt="2.png" width="534px" height="142px" src="cid:20c5d2a6-2cd4-4e40-b2a0-71afaa9b0354" /></div><div>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? <br /><br />Thanks in advanced.<br /></div><div><br /></div><div>#####==== cp2k ch4.inp ====#####</div><div>&GLOBAL<br /> PROJECT CH4<br /> RUN_TYPE ENERGY<br /> PRINT_LEVEL LOW<br />&END GLOBAL<br /><br />&FORCE_EVAL<br /> METHOD Quickstep ! Electronic structure method (DFT,...)<br /> &DFT<br /> BASIS_SET_FILE_NAME ./d-aug-cc-pvdz<br /> BASIS_SET_FILE_NAME ./d-aug-cc-pvdz-jkfit<br /> POTENTIAL_FILE_NAME POTENTIAL<br /> RESTART_FILE_NAME CH4-RESTART.wfn<br /> &MGRID<br /> CUTOFF 400<br /> NGRIDS 6<br /> REL_CUTOFF 60<br /> &END MGRID<br /> &POISSON ! Solver requested for non periodic calculations<br /> PERIODIC NONE<br /> POISSON_SOLVER WAVELET ! Type of solver<br /> &END POISSON<br /> &QS ! Parameters needed to set up the Quickstep framework<br /> METHOD GAPW ! Method: gaussian and augmented plane waves <br /> EPS_DEFAULT 1.0E-10<br /> &END QS<br /> &SCF<br /> SCF_GUESS RESTART<br /> EPS_SCF 1.0E-6<br /> MAX_SCF 100<br /> &OUTER_SCF<br /> MAX_SCF 50<br /> EPS_SCF 1.0E-6<br /> &END OUTER_SCF<br /> &OT<br /> PRECONDITIONER FULL_KINETIC<br /> MINIMIZER CG<br /> &END OT<br /> &END SCF<br /> &XC ! Parametes needed to compute the electronic exchange potential <br /> &XC_FUNCTIONAL<br /> &LYP<br /> SCALE_C 0.81<br /> &END <br /> &BECKE88<br /> SCALE_X 0.72<br /> &END<br /> &VWN<br /> FUNCTIONAL_TYPE VWN3<br /> SCALE_C 0.19<br /> &END<br /> &XALPHA<br /> SCALE_X 0.08<br /> &END<br /> &END XC_FUNCTIONAL<br /> &HF<br /> &SCREENING<br /> EPS_SCHWARZ 1.0E-10<br /> &END<br /> &MEMORY<br /> MAX_MEMORY 100<br /> &END<br /> FRACTION 0.20<br /> &END<br /> &END XC<br /> &PERIODIC_EFIELD<br /> INTENSITY 1.8897261250E-5<br /> POLARISATION -1.0 0.0 0.0 ！change the polarisation for x+ x- y+ y- z+ z-<br /> &END PERIODIC_EFIELD<br /> &PRINT<br /> &MOMENTS<br /> PERIODIC FALSE<br /> REFERENCE COM<br /> &EACH<br /> MD 1<br /> &END EACH<br /> FILENAME =dipole.txt<br /> &END MOMENTS<br /> &END PRINT<br /> &END DFT<br /><br /> &SUBSYS<br /> &CELL<br /> ABC 10. 10. 10.<br /> PERIODIC NONE ! Non periodic calculations. That's why the POISSON section is needed<br /> &END CELL<br /> &TOPOLOGY ! Section used to center the atomic coordinates in the given box. Useful for big molecules<br /> &CENTER_COORDINATES<br /> &END<br /> &END<br /> &COORD<br /> C 0.000003015000 0.000003640000 0.000013990000<br /> H 1.096282485000 0.000218630000 0.008833370000<br /> H -0.365745055000 1.033471870000 0.008832970000<br /> H -0.372522695000 -0.527106690000 0.886170220000<br /> H -0.358032825000 -0.506605650000 -0.903920500000<br /> &END COORD<br /> &KIND H ! potential and basis for H <br /> BASIS_SET d-aug-cc-pVDZ<br /> #BASIS_SET AUX_FIT d-aug-cc-pVDZ-JKFIT<br /> POTENTIAL ALL<br /> &END KIND<br /> &KIND C ! potential and basis for C<br /> BASIS_SET d-aug-cc-pVDZ<br /> #BASIS_SET AUX_FIT d-aug-cc-pVDZ-JKFIT<br /> POTENTIAL ALL<br /> &END KIND<br /> &END SUBSYS<br />&END FORCE_EVAL<br /></div><div><img alt="3.png" width="534px" height="260px" src="cid:1b6ffb49-efee-4990-9db3-66a1e3b18b81" /><br /></div><div>#####==== gaussian CH4.gif ====#####</div><div>%chk=C:\Users\hp\Desktop\spectra\gaussian_methanol\demo\demo.chk<br />#p b3lyp/gen nosymm polar<br /><br />demo<br /><br />0 1<br /> C 0.00000301 0.00000364 0.00001399<br /> H 1.09628248 0.00021863 0.00883337<br /> H -0.36574505 1.03347187 0.00883297<br /> H -0.37252270 -0.52710669 0.88617022<br /> H -0.35803282 -0.50660565 -0.90392050<br /><br />H 0<br />S 4 1.00<br /> 13.0100000 0.0196850<br /> 1.9620000 0.1379770<br /> 0.4446000 0.4781480<br /> 0.1220000 0.5012400<br />S 1 1.00<br /> 0.1220000 1.0000000<br />S 1 1.00<br /> 0.0297400 1.0000000<br />S 1 1.00<br /> 0.00725 1.000000<br />P 1 1.00<br /> 0.7270000 1.0000000<br />P 1 1.00<br /> 0.1410000 1.0000000<br />P 1 1.00<br /> 0.02730 1.000000<br />****<br />C 0<br />S 9 1.00<br /> 6665.0000000 0.0006920<br /> 1000.0000000 0.0053290<br /> 228.0000000 0.0270770<br /> 64.7100000 0.1017180<br /> 21.0600000 0.2747400<br /> 7.4950000 0.4485640<br /> 2.7970000 0.2850740<br /> 0.5215000 0.0152040<br /> 0.1596000 -0.0031910<br />S 9 1.00<br /> 6665.0000000 -0.0001460<br /> 1000.0000000 -0.0011540<br /> 228.0000000 -0.0057250<br /> 64.7100000 -0.0233120<br /> 21.0600000 -0.0639550<br /> 7.4950000 -0.1499810<br /> 2.7970000 -0.1272620<br /> 0.5215000 0.5445290<br /> 0.1596000 0.5804960<br />S 1 1.00<br /> 0.1596000 1.0000000<br />S 1 1.00<br /> 0.0469000 1.0000000<br />S 1 1.00<br /> 0.0138 1.000000<br />P 4 1.00<br /> 9.4390000 0.0381090<br /> 2.0020000 0.2094800<br /> 0.5456000 0.5085570<br /> 0.1517000 0.4688420<br />P 1 1.00<br /> 0.1517000 1.0000000<br />P 1 1.00<br /> 0.0404100 1.0000000<br />P 1 1.00<br /> 0.0108 1.000000<br />D 1 1.00<br /> 0.5500000 1.0000000<br />D 1 1.00<br /> 0.1510000 1.0000000<br />D 1 1.00<br /> 0.0415 1.000000<br />****<br /><br />#####==== psi4 CH4.in ====#####<br />memory 12000 mb<br /><br />molecule mol {<br />0 1<br />C<span style="white-space: pre;"> </span>.0000030150<span style="white-space: pre;"> </span>.0000036400<span style="white-space: pre;"> </span>.0000139900<br />H<span style="white-space: pre;"> </span>1.0962824850<span style="white-space: pre;"> </span>.0002186300<span style="white-space: pre;"> </span>.0088333700<br />H<span style="white-space: pre;"> </span>-.3657450550<span style="white-space: pre;"> </span>1.0334718700<span style="white-space: pre;"> </span>.0088329700<br />H<span style="white-space: pre;"> </span>-.3725226950<span style="white-space: pre;"> </span>-.5271066900<span style="white-space: pre;"> </span>.8861702200<br />H<span style="white-space: pre;"> </span>-.3580328250<span style="white-space: pre;"> </span>-.5066056500<span style="white-space: pre;"> </span>-.9039205000<br /><br />symmetry c1<br />no_reorient<br />no_com<br />}<br /><br />set {<br />basis d-aug-cc-pVDZ<br />guess sad<br />scf_type direct<br />df_basis_scf d-aug-cc-pvdz-jkfit<br />e_convergence 1.0E-10<br />d_convergence 1.0E-10<br />}<br /><br />import math<br />dipole_NE = []<br />quadrupole_NE = []<br />e, wfn = energy('B3LYP',return_wfn=True)<br /><br />dpx = psi4.get_variable('SCF DIPOLE X')<br />dpy = psi4.get_variable('SCF DIPOLE Y')<br />dpz = psi4.get_variable('SCF DIPOLE Z')<br />dptot = math.sqrt(dpx*dpx + dpy*dpy + dpz*dpz)<br />dipole_NE.append(dpx)<br />dipole_NE.append(dpy)<br />dipole_NE.append(dpz)<br />dipole_NE.append(dptot)<br />oeprop(wfn,'QUADRUPOLE',title="B3LYP")<br />qpxx = psi4.get_variable('B3LYP QUADRUPOLE XX')<br />qpxy = psi4.get_variable('B3LYP QUADRUPOLE XY')<br />qpyy = psi4.get_variable('B3LYP QUADRUPOLE YY')<br />qpxz = psi4.get_variable('B3LYP QUADRUPOLE XZ')<br />qpyz = psi4.get_variable('B3LYP QUADRUPOLE YZ')<br />qpzz = psi4.get_variable('B3LYP QUADRUPOLE ZZ')<br />quadrupole_NE.append(qpxx)<br />quadrupole_NE.append(qpxy)<br />quadrupole_NE.append(qpyy)<br />quadrupole_NE.append(qpxz)<br />quadrupole_NE.append(qpyz)<br />quadrupole_NE.append(qpzz)<br /><br />#--------------------------------------------------------------------<br />mints = psi4.core.MintsHelper(wfn.basisset())<br />ao_overlap = mints.ao_overlap()<br />Co = wfn.Ca()<br />import numpy<br />numpy.set_printoptions(threshold=numpy.nan)<br />Da_so = wfn.Da()<br />np.save('molecule0001_ao_overlap_E0.npy',ao_overlap)<br />np.save('molecule0001_so_density_E0.npy',Da_so)<br />np.save('molecule0001_mo_coefficients_E0.npy',Co)<br />#--------------------------------------------------------------------<br /><br />#perturbation<br />pert = 1.8897261250e-5<br />lambdas = [pert, -pert]<br /><br />dipole_EX = []<br />dipole_EY = []<br />dipole_EZ = []<br />set perturb_h true<br />set perturb_with dipole<br />sign = 'P'<br />for l in lambdas:<br /> set perturb_dipole [$l, 0, 0]<br /> set scf guess read<br /> e, wfn = energy('B3LYP',return_wfn=True)<br /> dpx = psi4.get_variable('SCF DIPOLE X')<br /> dpy = psi4.get_variable('SCF DIPOLE Y')<br /> dpz = psi4.get_variable('SCF DIPOLE Z')<br /> dipole_EX.append(dpx)<br /> dipole_EX.append(dpy)<br /> dipole_EX.append(dpz)<br />#--------------------------------------------------------------------<br /> mints = psi4.core.MintsHelper(wfn.basisset())<br /> ao_overlap = mints.ao_overlap()<br /> Co = wfn.Ca()<br /> Da_so = wfn.Da()<br /> Da_mo = Matrix.triplet(wfn.Ca(), Da_so, wfn.Ca(), True, False, False)<br /> np.save('molecule0001_ao_overlap_E' + sign + 'X.npy',ao_overlap)<br /> np.save('molecule0001_so_density_E' + sign + 'X.npy',Da_so)<br /> np.save('molecule0001_mo_coefficients_E' + sign + 'X.npy',Co)<br />#--------------------------------------------------------------------<br />#--------------------------------------------------------------------<br /> set perturb_dipole [0, $l, 0]<br /> set scf guess read<br /> e, wfn = energy('B3LYP',return_wfn=True)<br /> oeprop(wfn,'DIPOLE')<br /> dpx = psi4.get_variable('SCF DIPOLE X')<br /> dpy = psi4.get_variable('SCF DIPOLE Y')<br /> dpz = psi4.get_variable('SCF DIPOLE Z')<br /> dipole_EY.append(dpx)<br /> dipole_EY.append(dpy)<br /> dipole_EY.append(dpz)<br />#--------------------------------------------------------------------<br /> mints = psi4.core.MintsHelper(wfn.basisset())<br /> ao_overlap = mints.ao_overlap()<br /> Co = wfn.Ca()<br /> Da_so = wfn.Da()<br /> np.save('molecule0001_ao_overlap_E' + sign + 'Y.npy',ao_overlap)<br /> np.save('molecule0001_so_density_E' + sign + 'Y.npy',Da_so)<br /> np.save('molecule0001_mo_coefficients_E' + sign + 'Y.npy',Co)<br />#--------------------------------------------------------------------<br />#--------------------------------------------------------------------<br /> set perturb_dipole [0, 0, $l]<br /> set scf guess read<br /> e, wfn = energy('B3LYP',return_wfn=True)<br /> oeprop(wfn,'DIPOLE')<br /> dpx = psi4.get_variable('SCF DIPOLE X')<br /> dpy = psi4.get_variable('SCF DIPOLE Y')<br /> dpz = psi4.get_variable('SCF DIPOLE Z')<br /> dipole_EZ.append(dpx)<br /> dipole_EZ.append(dpy)<br /> dipole_EZ.append(dpz)<br />#--------------------------------------------------------------------<br /> mints = psi4.core.MintsHelper(wfn.basisset())<br /> ao_overlap = mints.ao_overlap()<br /> Co = wfn.Ca()<br /> Da_so = wfn.Da()<br /> np.save('molecule0001_ao_overlap_E' + sign + 'Z.npy',ao_overlap)<br /> np.save('molecule0001_so_density_E' + sign + 'Z.npy',Da_so)<br /> np.save('molecule0001_mo_coefficients_E' + sign + 'Z.npy',Co)<br /> <br /> sign = 'N'<br />#--------------------------------------------------------------------<br />#--------------------------------------------------------------------<br /><br />pl_11 = -(dipole_EX[0] - dipole_EX[3])/(2.0*pert*psi_dipmom_au2debye)<br />pl_12 = -(dipole_EX[1] - dipole_EX[4])/(2.0*pert*psi_dipmom_au2debye)<br />pl_13 = -(dipole_EX[2] - dipole_EX[5])/(2.0*pert*psi_dipmom_au2debye)<br /><br />pl_21 = -(dipole_EY[0] - dipole_EY[3])/(2.0*pert*psi_dipmom_au2debye)<br />pl_22 = -(dipole_EY[1] - dipole_EY[4])/(2.0*pert*psi_dipmom_au2debye)<br />pl_23 = -(dipole_EY[2] - dipole_EY[5])/(2.0*pert*psi_dipmom_au2debye)<br /><br />pl_31 = -(dipole_EZ[0] - dipole_EZ[3])/(2.0*pert*psi_dipmom_au2debye)<br />pl_32 = -(dipole_EZ[1] - dipole_EZ[4])/(2.0*pert*psi_dipmom_au2debye)<br />pl_33 = -(dipole_EZ[2] - dipole_EZ[5])/(2.0*pert*psi_dipmom_au2debye)<br /><br />print_out("\nYY-Dipole Moment [Debye] x+ x- y+ y- z+ z-")<br />print_out("\n")<br />psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EX[0],dipole_EX[1],dipole_EX[2]))<br />psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EX[3],dipole_EX[4],dipole_EX[5]))<br />psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EY[0],dipole_EY[1],dipole_EY[2]))<br />psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EY[3],dipole_EY[4],dipole_EY[5]))<br />psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EZ[0],dipole_EZ[1],dipole_EZ[2]))<br />psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(dipole_EZ[3],dipole_EZ[4],dipole_EZ[5]))<br /><br />print_out("\nYY-Polarizability tensor [a.u.]")<br />print_out("\n")<br />psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(pl_11,pl_12,pl_13))<br />psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(pl_21,pl_22,pl_23))<br />psi4.print_out("{:22.15f}{:25.15f}{:25.15f}\n".format(pl_31,pl_32,pl_33))<br /><br />print_out("\nYY-Dipole Moment(Debye)")<br />print_out("\n")<br />psi4.print_out("{}{:24.15f}{}{:22.15f}{}{:22.15f}\n".format("X",dipole_NE[0],"<span style="white-space: pre;"> </span>Y",dipole_NE[1],"<span style="white-space: pre;"> </span>Z",dipole_NE[2]))<br />psi4.print_out("{}{:22.15f}\n".format("Tol",dipole_NE[3]))<br /><br />print_out("\nYY-Quadrupole Moments (Debye-Ang)")<br />print_out("\n")<br />psi4.print_out("{}{:22.15f}{}{:22.15f}{}{:22.15f}\n".format("XX",quadrupole_NE[0],"<span style="white-space: pre;"> </span>XY",quadrupole_NE[1],"<span style="white-space: pre;"> </span>YY",quadrupole_NE[2]))<br />psi4.print_out("{}{:22.15f}{}{:22.15f}{}{:22.15f}\n".format("XZ",quadrupole_NE[3],"<span style="white-space: pre;"> </span>YZ",quadrupole_NE[4],"<span style="white-space: pre;"> </span>ZZ",quadrupole_NE[5]))<br /></div>
<p></p>
-- <br />
You received this message because you are subscribed to the Google Groups "cp2k" group.<br />
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>.<br />
To view this discussion on the web visit <a href="https://groups.google.com/d/msgid/cp2k/780660d3-ae84-4beb-b3c5-6c12afee7d66n%40googlegroups.com?utm_medium=email&utm_source=footer">https://groups.google.com/d/msgid/cp2k/780660d3-ae84-4beb-b3c5-6c12afee7d66n%40googlegroups.com</a>.<br />