From bf641204ca1b5288476d7cfdfa482dedb85a439d Mon Sep 17 00:00:00 2001 From: rubengerritsen Date: Fri, 26 Mar 2021 13:50:57 +0100 Subject: [PATCH 1/2] add numerical polarizability --- pyvotca/__init__.py | 1 + pyvotca/examples/run_polarizability.py | 30 +++++ pyvotca/numerical_polarizability.py | 150 +++++++++++++++++++++++++ 3 files changed, 181 insertions(+) create mode 100755 pyvotca/examples/run_polarizability.py create mode 100644 pyvotca/numerical_polarizability.py diff --git a/pyvotca/__init__.py b/pyvotca/__init__.py index 8aad76e..d97f5f1 100644 --- a/pyvotca/__init__.py +++ b/pyvotca/__init__.py @@ -4,6 +4,7 @@ from .__version__ import __version__ from .numerical_gradient import NumericalGradient +from .numerical_polarizability import NumericalPolarizability from .xtp import DFTGWBSE from .molecule import Molecule from .visualization import Visualization diff --git a/pyvotca/examples/run_polarizability.py b/pyvotca/examples/run_polarizability.py new file mode 100755 index 0000000..264b4c7 --- /dev/null +++ b/pyvotca/examples/run_polarizability.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python +"""Example to perform a polarizability calculation.""" +import numpy as np + +from pyvotca import DFTGWBSE, Molecule, NumericalPolarizability + + +def run_polarizability(): + mol = Molecule() + + # define a molecule + mol.add_atom("C", 0, 0, 0) + mol.add_atom("O", 1.3, 0.0, 0.0) + + # initialize a dftgwbse object + votca = DFTGWBSE(mol, threads=6) + + # setup the polarizability calculator + pol = NumericalPolarizability(votca, dE=0.002) + # run all simulations with different directions of the field + #pol.run_permut() + # perform the actual polarizability calculation + pol_tensor = pol.calc_polarizability("dft_tot") + + print(pol_tensor) + + +if __name__ == "__main__": + run_polarizability() + diff --git a/pyvotca/numerical_polarizability.py b/pyvotca/numerical_polarizability.py new file mode 100644 index 0000000..ad899f5 --- /dev/null +++ b/pyvotca/numerical_polarizability.py @@ -0,0 +1,150 @@ +"""Numerical polarizability. + +API +--- +.. autoclass:: NumericalPolarizability + +""" + +import numpy as np +import os +import itertools +import h5py + +from .molecule import Molecule +from .utils import BOHR2ANG +from .xtp import DFTGWBSE + +__all__ = ["NumericalPolarizability"] + + +class NumericalPolarizability: + """Compute the polarizability numerically using XTP.""" + def __init__(self, + xtp: DFTGWBSE, + dE: float = 0.002, + path_to_simulations: str = './polar_experiments/'): + self.xtp = xtp + self.dE = dE + self.path = path_to_simulations + + def mkfldr(self, path: str): + try: + os.mkdir(path) + except OSError: + print("Creation of the directory %s failed" % path) + + def gen_name(self, name, E) -> str: + return f"{name}_E{E[0]:.4f}E{E[1]:.4f}E{E[2]:.4f}" + + def run_permut(self): + """Run a VOTCA simulation for every displacement of the electric field with displacement dE in atomic units.""" + # compute all necessary directions of the field (permutations) + b = [self.dE, -self.dE, 0.0] + a = [b, b, b] + perm = [element for element in itertools.product(*a)] + perm = [s for s in perm if all(s) == False] + # create a folder to contain all the results from the experiments + if (not os.path.exists(self.path)): + self.mkfldr(self.path) + counter = 1 + for E in perm: + self.xtp.options.dftpackage = { + "package": { + "use_external_field": "true", + "externalfield": f"{E[0]} {E[1]} {E[2]}" + } + } + name = self.gen_name(self.xtp.mol.name, E) + xtp_withField = DFTGWBSE(self.xtp.mol, + threads=self.xtp.threads, + options=self.xtp.options, + jobname=name, + jobdir=self.path) + print( + f"Running XTP for field: X{E[0]:7.3f} Y{E[1]:7.3f} Z{E[2]:7.3f} {counter}/{len(perm)}" + ) + xtp_withField.run() + counter += 1 + + def getEnergiesFromOrbFile(self, kind: str, i: int, j: int, k: int): + """ Read the energies from the orb file for a certain kind of particle/excitation. + + i,j,k represents the direction of the electric field 1,0,0 corresponds to dE,0,0 etc. + OUTPUT: numpy array with the energies of the particle/excitation kind. + """ + nameOrbFile = self.path + self.gen_name( + self.xtp.mol.name, np.array([i * self.dE, j * self.dE, k * self.dE + ])) + ".orb" + f = h5py.File(nameOrbFile, 'r') + orb = f['QMdata'] + if (kind == 'BSE_singlet' or kind == 'BSE_triplet' or kind == 'QPdiag'): + return (np.array(orb[kind]['eigenvalues'][()])) + elif (kind == 'QPpert'): + return (np.array(orb['QPpert_energies'][()])) + elif (kind == 'dft_tot'): + return (np.array(orb.attrs['qm_energy'])) + else: + raise Exception( + 'Invalid kind, kind of particle/excitation should be one of:' + + ' BSE_singlet, BSE_triplet, QPdiag, QPpert and dft_tot\n' + + f'the value of kind was: {kind}') + + def calc_polarizability(self, kind: str, energyLevel=None) -> np.ndarray: + """ Computes the polarizability tensor for a particle/excitation type. + + INPUT: kind of particle/excitation (choices: BSE_singlet, BSE_triplet, + QPdiag, QPpert and dft_tot) and optionally the energy level if + not provided all energy levels will be returned + OUTPUT: numpy array of polarizability tensors. + """ + # Get all the energies from the files (encoding mop = -1 0 1) + ooo = self.getEnergiesFromOrbFile(kind, 0, 0, 0) + moo = self.getEnergiesFromOrbFile(kind, -1, 0, 0) + poo = self.getEnergiesFromOrbFile(kind, 1, 0, 0) + omo = self.getEnergiesFromOrbFile(kind, 0, -1, 0) + opo = self.getEnergiesFromOrbFile(kind, 0, 1, 0) + oom = self.getEnergiesFromOrbFile(kind, 0, 0, -1) + oop = self.getEnergiesFromOrbFile(kind, 0, 0, 1) + + mmo = self.getEnergiesFromOrbFile(kind, -1, -1, 0) + pmo = self.getEnergiesFromOrbFile(kind, 1, -1, 0) + mpo = self.getEnergiesFromOrbFile(kind, -1, 1, 0) + ppo = self.getEnergiesFromOrbFile(kind, 1, 1, 0) + + mom = self.getEnergiesFromOrbFile(kind, -1, 0, -1) + pom = self.getEnergiesFromOrbFile(kind, 1, 0, -1) + mop = self.getEnergiesFromOrbFile(kind, -1, 0, 1) + pop = self.getEnergiesFromOrbFile(kind, 1, 0, 1) + + omm = self.getEnergiesFromOrbFile(kind, 0, -1, -1) + opm = self.getEnergiesFromOrbFile(kind, 0, 1, -1) + omp = self.getEnergiesFromOrbFile(kind, 0, -1, 1) + opp = self.getEnergiesFromOrbFile(kind, 0, 1, 1) + + # Compute second derivatives of the energy with respect to the field + # strength + coeff = (-1. / (self.dE**2)) + axx = coeff * (moo - 2 * ooo + poo) + ayy = coeff * (omo - 2 * ooo + opo) + azz = coeff * (oom - 2 * ooo + oop) + axy = 0.25 * coeff * (mmo - pmo - mpo + ppo) + axz = 0.25 * coeff * (mom - pom - mop + pop) + ayz = 0.25 * coeff * (omm - opm - omp + opp) + + # Create polarization tensor + polar_m = np.squeeze(np.array([[axx, axy, axz], [axy, ayy, ayz], + [axz, ayz, azz]]), + axis=-1) + + # if the entries of the matrix are arrays we need to reshape to get an + # array of matrices + if (axx.size > 1): + polar_m = np.moveaxis(polar_m, [2], [0]) + else: # there can't be an energyLevel + energyLevel = None + + if (energyLevel is None): + return polar_m + else: + return polar_m[energyLevel] From e958d3043a625c8f537b7a5c45b8cca7c5a9ae3e Mon Sep 17 00:00:00 2001 From: rubengerritsen Date: Fri, 26 Mar 2021 14:54:02 +0100 Subject: [PATCH 2/2] minor update --- pyvotca/examples/run_polarizability.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyvotca/examples/run_polarizability.py b/pyvotca/examples/run_polarizability.py index 264b4c7..a172a98 100755 --- a/pyvotca/examples/run_polarizability.py +++ b/pyvotca/examples/run_polarizability.py @@ -13,7 +13,7 @@ def run_polarizability(): mol.add_atom("O", 1.3, 0.0, 0.0) # initialize a dftgwbse object - votca = DFTGWBSE(mol, threads=6) + votca = DFTGWBSE(mol) # setup the polarizability calculator pol = NumericalPolarizability(votca, dE=0.002)