"""
.. module:: Dip_Dip
:synopsis: Generate the Dipole_dipole interaction matrix.
.. moduleauthor:: D. Wang <dwang5@zoho.com>
"""
from math import *
import numpy as np
from netCDF4 import Dataset
import time
from ewald_dd import *
from utility import *
from supercell import Supercell
[docs]class Dip_dip(Supercell):
"""
The parameters are similar to that of 'charge-charge'.
"""
def __init__(self, n1, n2, nz, lattice):
Supercell.__init__(self, n1, n2, nz, lattice)
self.dpij = np.zeros((self.nsites, 6))
self.dipole_matrix_calculated = False
def write_dipole_matrix(self, fn):
if self.dipole_matrix_calculated == False:
print("Calculate the matrix first ...")
self.generate_dipole_matrix()
self.dipole_matrix_calculated = True
dipm = Dataset(fn, "w", format="NETCDF4")
direction = dipm.createDimension("direction", 6)
ia = dipm.createDimension("ia", None)
directions = dipm.createVariable("direction", np.int32, ("direction"))
ias = dipm.createVariable("ia", np.int32, ("ia"))
# The actual 2-d varable.
matrix = dipm.createVariable('matrix', np.float64, ('ia', 'direction'))
dipm.description = 'Dipole matrix: interaction matrix'
dipm.history = 'Created at ' + time.ctime(time.time())
matrix[:, :] = self.dpij
dipm.close()
def generate_dipole_matrix(self):
pi = 4.0 * atan(1.0)
pi2 = pi * 2.0
tol = 1.0e-12
eta = sqrt(-log(tol))
gcut = 2.0 * eta ** 2
gcut2 = gcut ** 2
eta4 = 1.0 / (4 * eta ** 2)
am = np.zeros((3))
for i in range(3):
for k in range(3):
am[i] += self.a[k, i] ** 2
am[i] = sqrt(am[i])
mg1 = int(gcut * am[0] / pi2) + 1
mg2 = int(gcut * am[1] / pi2) + 1
mg3 = int(gcut * am[2] / pi2) + 1
print('Gcut: ', gcut, ' mg1, mg2, mg3: ', mg1, mg2, mg3)
# Set parameters to be used in the PYBIND11 C++ computation.
dd_set_parameters(self.b, mg1, mg2, mg3, gcut2, eta4)
pos0 = np.zeros(3)
pos0 = self.ixa[0] * self.lattice[0, :] \
+ self.iya[0] * self.lattice[1, :] \
+ self.iza[0] * self.lattice[2, :]
print(pos0)
for ia in range(self.nsites):
print('site: ', ia)
pos = np.zeros(3)
pos = self.ixa[ia] * self.lattice[0, :] \
+ self.iya[ia] * self.lattice[1, :] \
+ self.iza[ia] * self.lattice[2, :]
rx = pos[0]
ry = pos[1]
rz = pos[2]
origin = (rx == 0) and (ry == 0) and (rz == 0)
c = 2.0 * pi / self.celvol
residue = 2.0 * eta ** 3 / (3.0 * sqrt(pi))
dum = np.zeros(6)
# print('Summing over k space')
dd_sum_over_k(dum, rx - pos0[0], ry - pos0[1], rz - pos0[2])
# The 2.0 comes from how the sum on k was done.
self.dpij[ia, :] = dum[:] * c * 2.0
if (origin):
self.dpij[ia, 0:3] -= residue
self.dipole_matrix_calculated = True