Question on calculation of Dipole

[kaypu]

Dear quantunwise staff：
      Here is the molecule system in the attachment，we calculate the dipole in Gaussian and ATK, respectively.
For atk（GGA-PBE，DZP）：dipole     (e*bohr)    = 5.3662629780200824e-14, -7.290162499260885e-13, -2.6478867127832016e-12
For Gaussian（B3LYP，6-31g*）：dipole     (debye)  9.03655361E-01, -8.41922597E-01, -1.92566377E-02

PS: The structure optimization is performed in Gaussian(B3LYP，6-31g*). Then, put the relaxed molecular configuration to atk and Gaussian for electronic structure calculation.

why such a big different?

[mlee]

Can you describe how to calculate the dipole from your script in detail? 

[kaypu]

the script is as follow (https://docs.quantumatk.com/manual/Types/ElectronDensity/ElectronDensity.html):
***************************************
# import an electron density
path = u'/home/Electron Density.nc'
density = nlread(path, ElectronDensity)[0]

# Get the shape of the data.
ni, nj, nk = density.shape()
# Find the volume elements.
dX, dY, dZ = density.volumeElement()
length_unit = dX.unit()
# Calculate the volume of the volume element.
dV = numpy.dot(dX, numpy.cross(dY,dZ)) * length_unit**3
# Get a reference for the data of the density.
grid_data = density[:,:,:]*Units.e

# Calculate the absoute density sum.
density_abs_sum = abs(grid_data).sum()

# Calculate center of mass
Mi = sum([ i*abs(grid_data[i,:,:]).sum() for i in range(ni)])/density_abs_sum
Mj = sum([ j*abs(grid_data[:,j,:]).sum() for j in range(nj)])/density_abs_sum
Mk = sum([ k*abs(grid_data[:,:,k]).sum() for k in range(nk)])/density_abs_sum
center_of_mass = Mi*dX + Mj*dY+ Mk*dZ

# Calculate monopole
density_sum = grid_data.sum()
m0 = density_sum*dV

#calculate dipole
dens_i = sum([ (i-Mi)*grid_data[i,:,:].sum() for i in range(ni)])
dens_j = sum([ (j-Mj)*grid_data[:,j,:].sum() for j in range(nj)])
dens_k = sum([ (k-Mk)*grid_data[:,:,k].sum() for k in range(nk)])
m1 = (dX*dens_i + dY*dens_j + dZ*dens_k)*dV

# Calculate quadropoles (3 x_i *x_j-r*r delta_ij)*dens
dens_ii = sum([ (i-Mi)**2 * grid_data[i,:,:].sum() for i in range(ni)])
dens_jj = sum([ (j-Mj)**2 * grid_data[:,j,:].sum() for j in range(nj)])
dens_kk = sum([ (k-Mk)**2 * grid_data[:,:,k].sum() for k in range(nk)])
dens_ij = sum([ (i-Mi)*(j-Mj)*grid_data[i,j,:].sum()
             for i in range(ni) for j in range(nj)])
dens_ik = sum([ (i-Mi)*(k-Mk)*grid_data[i,:,k].sum()
             for i in range(ni) for k in range(nk)])
dens_jk = sum([ (j-Mj)*(k-Mk)*grid_data[:,j,k].sum()
             for j in range(nj) for k in range(nk)])

m2 = (dX.outer(dX)*dens_ii + dY.outer(dY)*dens_jj + dZ.outer(dZ)*dens_kk\
   + (dX.outer(dY)+dY.outer(dX))*dens_ij + (dX.outer(dZ)+dZ.outer(dX))*dens_ik\
   + (dY.outer(dZ)+dZ.outer(dY))*dens_jk) * dV

r2 = sum(m2)
m2 = 3*m2 - numpy.identity(3)*r2

print("center of mass (bohr)  = ", center_of_mass)

print("monopole   (e)         = ", m0)
print("dipole     (e*bohr)    = ", m1)
print("quadropole (e*bohr**2) = ", m2)
***********************************************************

[Petr Khomyakov]

Could you try using GGA-PBE in Gaussian to see what dipole values you get, and them compare to the QATK results. Otherwise, this is a comparison between 2 calculations done with very different settings?