diff --git a/.github/workflows/build_conda.yml b/.github/workflows/build_conda.yml new file mode 100644 index 0000000..29a726d --- /dev/null +++ b/.github/workflows/build_conda.yml @@ -0,0 +1,81 @@ +name: Conda + +on: + workflow_dispatch: + pull_request: + push: + branches: + - master + - develop + tags: + - "v**" + release: + types: + - published + +jobs: + check-version-strings: + runs-on: ubuntu-latest + container: python:3.11-slim-bullseye + + steps: + - name: Checkout repo + uses: actions/checkout@v3 + - name: Install dependencies + run: | + pip install pyyaml + - name: Test versions + run: | + python versiontest.py + + build-anaconda: + needs: check-version-strings + runs-on: ubuntu-latest + container: continuumio/miniconda3:23.10.0-1 + + steps: + - name: Checkout repo + uses: actions/checkout@v3 + - name: Build + run: | + rm pyproject.toml + conda install conda-build + conda install -c ifilot pyqint pylebedev + conda install -c conda-forge tqdm + conda build . --output-folder conda-bld/ + - name: Archive packages + uses: actions/upload-artifact@v3 + with: + name: anaconda-packages + path: conda-bld/noarch/* + + deploy-anaconda: + name: Publish Anaconda / Windows + if: startsWith(github.ref, 'refs/tags/v') + needs: build-anaconda + runs-on: ubuntu-latest + environment: + name: anaconda + url: https://anaconda.org/ifilot/pyqint + steps: + - name: Checkout repo + uses: actions/checkout@v3 + - name: Set-up miniconda + uses: conda-incubator/setup-miniconda@v2 + with: + activate-environment: test + environment-file: environment.yml + python-version: 3.8 + auto-activate-base: false + - name: Retrieve packages + uses: actions/download-artifact@v3 + with: + name: anaconda-packages + path: packages + - name: publish-to-conda + shell: bash -l {0} + env: + INPUT_ANACONDATOKEN: ${{ secrets.ANACONDA_TOKEN }} + run: | + export ANACONDA_API_TOKEN=$INPUT_ANACONDATOKEN + anaconda upload packages/*.tar.bz2 \ No newline at end of file diff --git a/.github/workflows/build_pypi.yml b/.github/workflows/build_pypi.yml new file mode 100644 index 0000000..0f7728d --- /dev/null +++ b/.github/workflows/build_pypi.yml @@ -0,0 +1,91 @@ +name: PyPI + +on: + workflow_dispatch: + pull_request: + push: + branches: + - master + - develop + tags: + - "v**" + release: + types: + - published + +jobs: + check-version-strings: + name: Check version strings + runs-on: ubuntu-latest + container: python:3.11-slim-bullseye + + steps: + - name: Checkout repo + uses: actions/checkout@v3 + - name: Install dependencies + run: | + pip install pyyaml + - name: Test versions + run: | + python versiontest.py + + build-pip: + name: Build pip + needs: [check-version-strings] + runs-on: ubuntu-latest + container: python:3.11-slim-bullseye + + steps: + - uses: actions/checkout@v4 + - name: Build dependencies + run: | + pip install virtualenv + virtualenv venv + . venv/bin/activate + pip install build + - name: Build WHL + run: | + . venv/bin/activate + python -m build + - uses: actions/upload-artifact@v3 + with: + path: ./dist/*.whl + name: whl + + test-pip: + name: Perform unit testing + needs: [build-pip] + runs-on: ubuntu-latest + container: python:3.11-slim-bullseye + + steps: + - uses: actions/checkout@v4 + - uses: actions/download-artifact@v3 + with: + name: whl + path: dist + - name: Install PyDFT + run: | + pip install virtualenv + virtualenv venv + . venv/bin/activate + pip install build pytest dist/*.whl + - name: Perform unit tests + run: | + . venv/bin/activate + python -m pytest tests/*.py --verbose + + upload_pypi: + needs: [build-pip, test-pip] + runs-on: ubuntu-latest + environment: pypi + permissions: + id-token: write + if: startsWith(github.ref, 'refs/tags/v') + steps: + - uses: actions/download-artifact@v3 + with: + name: whl + path: dist + + - uses: pypa/gh-action-pypi-publish@release/v1 \ No newline at end of file diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 0000000..b7d09fa --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,64 @@ +name: Docs + +on: + push: + branches: [ "docs" ] + pull_request: + branches: [ "docs" ] + +jobs: + build-docs: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - name: Install packages + run: sudo apt update && sudo apt install -y git build-essential curl wget python3 python3-pip texlive texlive-pictures texlive-latex-extra pdf2svg poppler-utils netpbm imagemagick ghostscript + - name: Install Python packages + run: sudo pip install sphinx sphinx-rtd-theme sphinxcontrib-tikz + - name: Build documentation + run: cd docs && make html + - name: Upload math result for job 1 + uses: actions/upload-artifact@v3 + with: + name: html-docs + path: ./docs/_build/html + + deploy: + runs-on: ubuntu-latest + needs: build-docs + permissions: + contents: read + packages: write + + steps: + - name: Checkout repository + uses: actions/checkout@v3 + + - name: Download docs + uses: actions/download-artifact@v3 + with: + name: html-docs + path: html-docs + + - name: Extract metadata (tags, labels) for Docker + id: meta + uses: docker/metadata-action@v4 + with: + images: | + ghcr.io/ifilot/pydft-docs-nginx + + - name: Log in to the Container registry + uses: docker/login-action@v2 + with: + registry: ghcr.io + username: ${{ github.actor }} + password: ${{ secrets.GITHUB_TOKEN }} + + - name: Build and push Docker images + uses: docker/build-push-action@v4 + with: + context: . + file: doc.Dockerfile + push: true + tags: ${{ steps.meta.outputs.tags }} + labels: ${{ steps.meta.outputs.labels }} \ No newline at end of file diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0b7badf --- /dev/null +++ b/.gitignore @@ -0,0 +1,152 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ +*.xyz + +# Spynx documentation +docs/build/* + +# build folder +build/* +dist/* +wheelhouse/* +test.py +nose_debug +anaconda-upload/*.tar.bz2 +!pydft/molecules/*.xyz +.vscode diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml new file mode 100644 index 0000000..a890cd1 --- /dev/null +++ b/.gitlab-ci.yml @@ -0,0 +1,204 @@ +variables: + PIP_CACHE_DIR: "$CI_PROJECT_DIR/.cache/pip" + +################################################################################ +# Check that all three versions match +################################################################################ +check-versions: + image: ubuntu:22.04 + cache: + paths: + - .cache/pip + before_script: + - apt update && apt install -y python3 python3-pip + - python3 --version ; pip3 --version # for debugging purposes + - pip install virtualenv + - virtualenv venv + - source venv/bin/activate + script: + - python3 versiontest.py + +################################################################################ +# PIP +################################################################################ +build-pip: + image: ubuntu:22.04 + cache: + paths: + - .cache/pip + needs: + - job: check-versions + before_script: + - apt update && apt install -y python3 python3-pip + - python3 --version ; pip3 --version # for debugging purposes + - pip3 install virtualenv + - virtualenv venv + - source venv/bin/activate + - pip3 install build + script: + - python3 -m build + artifacts: + paths: + - dist/* + +test-pip: + image: ubuntu:22.04 + needs: + - job: build-pip + artifacts: true + before_script: + - apt update && apt install -y python3 python3-pip + - python3 --version ; pip3 --version # for debugging purposes + - pip3 install virtualenv + - virtualenv venv + - source venv/bin/activate + - pip3 install dist/pydft-*.whl + - pip3 install numpy scipy matplotlib pytest pyqint pylebedev mendeleev tqdm pytest + script: + - python3 -m pytest tests/* --verbose + +code-coverage: + image: ubuntu:22.04 + cache: + paths: + - .cache/pip + needs: + - job: test-pip + before_script: + - apt update && apt install -y python3 python3-pip + - python3 --version ; pip3 --version # for debugging purposes + - pip3 install virtualenv + - virtualenv venv + - source venv/bin/activate + - pip3 install numpy scipy matplotlib pytest pyqint pylebedev mendeleev tqdm coverage pytest + script: + - coverage run -m pytest tests/* + - coverage xml pydft/*.py + - coverage report pydft/*.py + coverage: '/(?i)total.*? (100(?:\.0+)?\%|[1-9]?\d(?:\.\d+)?\%)$/' + artifacts: + reports: + coverage_report: + coverage_format: cobertura + path: coverage.xml + +deploy-pip: + image: ubuntu:22.04 + cache: + paths: + - .cache/pip + before_script: + - apt update && apt install -y python3 python3-pip + - python3 --version ; pip3 --version # for debugging purposes + - pip3 install virtualenv + - virtualenv venv + - source venv/bin/activate + - pip3 install twine + needs: + - job: build-pip + artifacts: true + - job: test-pip + script: + - twine upload -u __token__ -p ${pypi_token} dist/* + only: + - tags + except: + - branches + +deployment-pip-verification: + image: ubuntu:22.04 + cache: + paths: + - .cache/pip + before_script: + - apt update && apt install -y python3 python3-pip + - python3 --version ; pip3 --version # for debugging purposes + - pip3 install virtualenv + - virtualenv venv + - source venv/bin/activate + - pip3 install numpy scipy matplotlib pytest pydft + needs: + - job: deploy-pip + script: + - python3 -m pytest tests/* + only: + - tags + except: + - branches + +################################################################################ +# ANACONDA +################################################################################ +build-conda: + image: continuumio/miniconda3:23.10.0-1 + needs: + - job: check-versions + before_script: + - rm pyproject.toml # ignore pyproject.toml file + - conda install conda-build + - conda install -c ifilot pyqint pylebedev # needed for testing + - conda install -c conda-forge tqdm # needed for testing + script: + - conda build . --output-folder conda-bld/ + artifacts: + paths: + - conda-bld/noarch/* + +deploy-anaconda: + image: continuumio/miniconda3:23.10.0-1 + before_script: + - rm pyproject.toml # ignore pyproject.toml file + - conda install anaconda-client + needs: + - job: build-conda + artifacts: true + script: + - anaconda -t ${anaconda_token} upload -u ${anaconda_user} conda-bld/noarch/pydft-*.tar.bz2 + only: + - tags + except: + - branches + +deployment-anaconda-verification: + image: continuumio/miniconda3:23.10.0-1 + before_script: + - conda install -c ifilot pyqint pylebedev pydft pytest + needs: + - job: deploy-anaconda + script: + - python -m pytest tests/* + only: + - tags + except: + - branches + +################################################################################ +# Documentation +################################################################################ + +# build the Sphinx documentation +build-docs: + image: ghcr.io/ifilot/sphinx-qm:v0.2.0 + rules: + - if: '$CI_COMMIT_BRANCH == "docs"' + script: + - cd docs + - make html + artifacts: + paths: + - docs/_build/html/ + expire_in: 1 hour + +# copy the documentation over to Gitlab pages +pages: # do not change this name + rules: + - if: '$CI_COMMIT_BRANCH == "docs"' + needs: + - job: build-docs + artifacts: true + script: + - cp -rv docs/_build/html public + artifacts: + paths: + - public + expire_in: 1 week \ No newline at end of file diff --git a/COMPILATION.md b/COMPILATION.md new file mode 100644 index 0000000..9312262 --- /dev/null +++ b/COMPILATION.md @@ -0,0 +1,27 @@ +# Compilation details + +## Compilation and testing under Linux Debian / Ubuntu + +Compile locally + +```python +python3 -m build +``` + +and install it locally + +```python +pip3 install dist/pydft--py3-none-any.whl +``` + +and finally test it + +```python +python3 -m pytest tests/* +``` + +## Outputting content of Python Wheel file + +```bash +unzip -l dist/*.whl +``` \ No newline at end of file diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..b789a17 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,92 @@ +# Contributing + +When contributing to this repository, please first discuss the change you wish to make via issue, +email, or any other method with the owners of this repository before making a change. + +Please note we have a code of conduct, please follow it in all your interactions with the project. + +## Pull Request Process + +1. Ensure any install or build dependencies are removed before the end of the layer when doing a + build. +2. Update the README.md with details of changes to the interface, this includes new environment + variables, exposed ports, useful file locations and container parameters. +3. Increase the version numbers in any examples files and the README.md to the new version that this + Pull Request would represent. The versioning scheme we use is [SemVer](http://semver.org/). +4. You may merge the Pull Request in once you have the sign-off of two other developers, or if you + do not have permission to do that, you may request the second reviewer to merge it for you. + +## Code of Conduct + +### Our Pledge + +In the interest of fostering an open and welcoming environment, we as +contributors and maintainers pledge to making participation in our project and +our community a harassment-free experience for everyone, regardless of age, body +size, disability, ethnicity, gender identity and expression, level of experience, +nationality, personal appearance, race, religion, or sexual identity and +orientation. + +### Our Standards + +Examples of behavior that contributes to creating a positive environment +include: + +* Using welcoming and inclusive language +* Being respectful of differing viewpoints and experiences +* Gracefully accepting constructive criticism +* Focusing on what is best for the community +* Showing empathy towards other community members + +Examples of unacceptable behavior by participants include: + +* The use of sexualized language or imagery and unwelcome sexual attention or +advances +* Trolling, insulting/derogatory comments, and personal or political attacks +* Public or private harassment +* Publishing others' private information, such as a physical or electronic + address, without explicit permission +* Other conduct which could reasonably be considered inappropriate in a + professional setting + +### Our Responsibilities + +Project maintainers are responsible for clarifying the standards of acceptable +behavior and are expected to take appropriate and fair corrective action in +response to any instances of unacceptable behavior. + +Project maintainers have the right and responsibility to remove, edit, or +reject comments, commits, code, wiki edits, issues, and other contributions +that are not aligned to this Code of Conduct, or to ban temporarily or +permanently any contributor for other behaviors that they deem inappropriate, +threatening, offensive, or harmful. + +### Scope + +This Code of Conduct applies both within project spaces and in public spaces +when an individual is representing the project or its community. Examples of +representing a project or community include using an official project e-mail +address, posting via an official social media account, or acting as an appointed +representative at an online or offline event. Representation of a project may be +further defined and clarified by project maintainers. + +### Enforcement + +Instances of abusive, harassing, or otherwise unacceptable behavior may be +reported by contacting the project team at [INSERT EMAIL ADDRESS]. All +complaints will be reviewed and investigated and will result in a response that +is deemed necessary and appropriate to the circumstances. The project team is +obligated to maintain confidentiality with regard to the reporter of an incident. +Further details of specific enforcement policies may be posted separately. + +Project maintainers who do not follow or enforce the Code of Conduct in good +faith may face temporary or permanent repercussions as determined by other +members of the project's leadership. + +### Attribution + +This Code of Conduct is adapted from the [Contributor Covenant][homepage], version 1.4, +available at [http://contributor-covenant.org/version/1/4][version] + +[homepage]: http://contributor-covenant.org +[version]: http://contributor-covenant.org/version/1/4/ diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..b4998ac --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1 @@ +include pydft/molecules/*.xyz \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..a204605 --- /dev/null +++ b/README.md @@ -0,0 +1,66 @@ +# PyDFT + +[![build](https://github.com/ifilot/pydft/actions/workflows/build_pypi.yml/badge.svg)](https://github.com/ifilot/pydft/actions/workflows/build_pypi.yml) +[![build](https://github.com/ifilot/pydft/actions/workflows/build_conda.yml/badge.svg)](https://github.com/ifilot/pydft/actions/workflows/build_conda.yml) +[![Anaconda-Server Badge](https://anaconda.org/ifilot/pydft/badges/version.svg)](https://anaconda.org/ifilot/pydft) +[![PyPI](https://img.shields.io/pypi/v/pydft?color=green)](https://pypi.org/project/pytessel/) +[![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) + +Python based Density Functional Theory code for educational purposes. The +documentation of PyDFT can be found [here](https://ifilot.pages.tue.nl/pydft/). + +## Purpose + +PyDFT is a pure-Python package for performing localized-orbital DFT calculations +using Gaussian Type Orbitals. PyDFT currently supports LDA and PBE +exchange-correlation functionals. The purpose of PyDFT is mainly to serve as an +educational tool to explain the inner workings of a DFT calculation. This +program is not intended for professional calculations. It is not particularly +fast nor offers a lot of features that more mature open-source of commercial +packages offer. It does offer a unique insight into a working code and a +considerable effort was made in documenting everything. + +## Installation + +This code depends on a few other packages. To install this code and its +dependencies, run the following one-liner from Anaconda prompt + +```bash +conda install -c ifilot pydft pyqint pylebedev pytessel +``` + +## Usage + +### Check the current version + +```python +import pydft +print(pydft.__version__) +``` + +### Performing a simple calculation + +```python +from pydft import MoleculeBuilder, DFT + +CO = MoleculeBuilder().get_molecule("CO") +dft = DFT(CO, basis='sto3g', verbose=True) +en = dft.scf(1e-4) +print("Total electronic energy: %f Ht" % en) +``` + +## Community guidelines + +* Contributions to PyDFT are always welcome and appreciated. Before doing so, + please first read the [CONTRIBUTING](CONTRIBUTING.md) guide. +* For reporting issues or problems with the software, you are kindly invited to + to open a [new issue with the bug label](https://github.com/ifilot/pydft/issues/new?labels=bug). +* If you seek support in using PyDFT, please + [open an issue with the question](https://github.com/ifilot/pydft/issues/new?labels=question) + label. +* If you wish to contact the developers, please send an e-mail to ivo@ivofilot.nl. + +## License + +Unless otherwise stated, all code in this repository is provided under the GNU +General Public License version 3. \ No newline at end of file diff --git a/build.sh b/build.sh new file mode 100644 index 0000000..fec5047 --- /dev/null +++ b/build.sh @@ -0,0 +1 @@ +$PYTHON setup.py install # Python command to install the script. \ No newline at end of file diff --git a/doc.Dockerfile b/doc.Dockerfile new file mode 100644 index 0000000..000052c --- /dev/null +++ b/doc.Dockerfile @@ -0,0 +1,2 @@ +FROM nginx:latest +COPY html-docs/ /usr/share/nginx/html/ \ No newline at end of file diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d4bb2cb --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/_static/css/custom.css b/docs/_static/css/custom.css new file mode 100644 index 0000000..aa5e1e7 --- /dev/null +++ b/docs/_static/css/custom.css @@ -0,0 +1,11 @@ +.tight-table td { + white-space: normal !important; +} + +.wy-side-nav-search { + background: #def4fa; +} + +.icon.icon-home { + color: #1292b3; +} diff --git a/docs/_static/img/favicon.ico b/docs/_static/img/favicon.ico new file mode 100644 index 0000000..2a6a7d4 Binary files /dev/null and b/docs/_static/img/favicon.ico differ diff --git a/docs/_static/img/pydft_128px.png b/docs/_static/img/pydft_128px.png new file mode 100644 index 0000000..fbe36d5 Binary files /dev/null and b/docs/_static/img/pydft_128px.png differ diff --git a/docs/_static/img/pydft_192px.png b/docs/_static/img/pydft_192px.png new file mode 100644 index 0000000..c9ff974 Binary files /dev/null and b/docs/_static/img/pydft_192px.png differ diff --git a/docs/_static/img/pydft_256px.png b/docs/_static/img/pydft_256px.png new file mode 100644 index 0000000..bd6d30b Binary files /dev/null and b/docs/_static/img/pydft_256px.png differ diff --git a/docs/_static/img/pydft_500px.png b/docs/_static/img/pydft_500px.png new file mode 100644 index 0000000..ce4a1be Binary files /dev/null and b/docs/_static/img/pydft_500px.png differ diff --git a/docs/_static/img/pydft_logo_128px.png b/docs/_static/img/pydft_logo_128px.png new file mode 100644 index 0000000..ca64169 Binary files /dev/null and b/docs/_static/img/pydft_logo_128px.png differ diff --git a/docs/_static/img/pydft_logo_512px.png b/docs/_static/img/pydft_logo_512px.png new file mode 100644 index 0000000..73f393b Binary files /dev/null and b/docs/_static/img/pydft_logo_512px.png differ diff --git a/docs/_static/img/pydft_logo_full_512px.png b/docs/_static/img/pydft_logo_full_512px.png new file mode 100644 index 0000000..9a596d9 Binary files /dev/null and b/docs/_static/img/pydft_logo_full_512px.png differ diff --git a/docs/_static/img/user_interface/01-electron-density.png b/docs/_static/img/user_interface/01-electron-density.png new file mode 100644 index 0000000..19e9007 Binary files /dev/null and b/docs/_static/img/user_interface/01-electron-density.png differ diff --git a/docs/_static/img/user_interface/02-electron-density-gradient.png b/docs/_static/img/user_interface/02-electron-density-gradient.png new file mode 100644 index 0000000..f297339 Binary files /dev/null and b/docs/_static/img/user_interface/02-electron-density-gradient.png differ diff --git a/docs/_static/img/user_interface/03-matrices.png b/docs/_static/img/user_interface/03-matrices.png new file mode 100644 index 0000000..265d9db Binary files /dev/null and b/docs/_static/img/user_interface/03-matrices.png differ diff --git a/docs/_static/img/user_interface/04-kohn-sham-orbitals.png b/docs/_static/img/user_interface/04-kohn-sham-orbitals.png new file mode 100644 index 0000000..a3ae125 Binary files /dev/null and b/docs/_static/img/user_interface/04-kohn-sham-orbitals.png differ diff --git a/docs/_static/img/user_interface/05-orbital-contour-plot.png b/docs/_static/img/user_interface/05-orbital-contour-plot.png new file mode 100644 index 0000000..86d3d23 Binary files /dev/null and b/docs/_static/img/user_interface/05-orbital-contour-plot.png differ diff --git a/docs/_static/img/user_interface/06-becke-grid.png b/docs/_static/img/user_interface/06-becke-grid.png new file mode 100644 index 0000000..8f3df0d Binary files /dev/null and b/docs/_static/img/user_interface/06-becke-grid.png differ diff --git a/docs/_static/img/user_interface/06-grid-points.png b/docs/_static/img/user_interface/06-grid-points.png new file mode 100644 index 0000000..84cd621 Binary files /dev/null and b/docs/_static/img/user_interface/06-grid-points.png differ diff --git a/docs/_static/img/user_interface/07-spherical-harmonics-projection.png b/docs/_static/img/user_interface/07-spherical-harmonics-projection.png new file mode 100644 index 0000000..8520751 Binary files /dev/null and b/docs/_static/img/user_interface/07-spherical-harmonics-projection.png differ diff --git a/docs/_static/img/user_interface/08-correlation-potential.png b/docs/_static/img/user_interface/08-correlation-potential.png new file mode 100644 index 0000000..2db7665 Binary files /dev/null and b/docs/_static/img/user_interface/08-correlation-potential.png differ diff --git a/docs/_static/img/user_interface/08-exchange-potential.png b/docs/_static/img/user_interface/08-exchange-potential.png new file mode 100644 index 0000000..594185e Binary files /dev/null and b/docs/_static/img/user_interface/08-exchange-potential.png differ diff --git a/docs/_templates/breadcrumbs.html b/docs/_templates/breadcrumbs.html new file mode 100644 index 0000000..6184fee --- /dev/null +++ b/docs/_templates/breadcrumbs.html @@ -0,0 +1,4 @@ +
+ Gitlab +
+
diff --git a/docs/api.rst b/docs/api.rst new file mode 100644 index 0000000..eb113f9 --- /dev/null +++ b/docs/api.rst @@ -0,0 +1,21 @@ +.. _api: +.. index:: api + +API +=== + +.. automodule:: pydft + + DFT + --- + + .. autoclass:: DFT + :members: + :special-members: __init__ + + MolecularGrid + ------------- + + .. autoclass:: MolecularGrid + :members: + :special-members: __init__ \ No newline at end of file diff --git a/docs/background.rst b/docs/background.rst new file mode 100644 index 0000000..a881a4f --- /dev/null +++ b/docs/background.rst @@ -0,0 +1,42 @@ +.. _background: +.. index:: Background + +Background +========== + +:program:`PyDFT` is a pure-Python package for performing DFT calculations, +extending upon the functionality of `PyQInt `_ +and leveraging a number of other packages such as +`PyLebedev `_ and +`PyTessel `_ for quadrature on the unit sphere +and generation of isosurfaces, respectively. + +.. tip:: + + More information on the inner workings of :program:`PyDFT` can be obtained + from the textbook "Elements of Electronic Structure Theory" (specifically + chapter 4), which is freely available `via this website `_. + +Molecular decomposition +----------------------- + +Solving the integrals involved in the electronic structure calculation is handled +by means of numerical integration, also termed quadrature. The quadratures are +solved by decomposing the molecule into so-called "fuzzy" cells as documented +in the work of Becke. + +Hartree potential +----------------- + +Electron-electron repulsion is handled by calculating the Hartree potential +by means of solving Poisson's equation. This equation is solved per fuzzy cell, +as detailed in the seminal paper of Becke. + + +Exchange-correlation functions +------------------------------ + +:program:`PyDFT` currently supports two exchange-correlation functions: + +* LDA: Slater exchange + SVWN5 for the correlation +* PBE: The (standard) Perdew-Burke-Ernzerhof exchange-correlation functional \ No newline at end of file diff --git a/docs/community_guidelines.rst b/docs/community_guidelines.rst new file mode 100644 index 0000000..734c38f --- /dev/null +++ b/docs/community_guidelines.rst @@ -0,0 +1,12 @@ +.. _community_guidelines: +.. index:: Community Guidelines + +Community guidelines +******************** + +* Contributions to :program:`PyDFT` are always welcome and appreciated. Before doing + so, please first read the `CONTRIBUTING `_ + guide. +* For reporting issues or problems with the software, you are kindly invited to + `to open a new issue on Gitlab `_. +* If you wish to contact the developers, please send an e-mail to i.a.w.filot@tue.nl. diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..1a7728d --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,75 @@ +# Configuration file for the Sphinx documentation builder. +# +# This file only contains a selection of the most common options. For a full +# list see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +sys.path.insert(0, os.path.abspath('..')) +import sphinx_rtd_theme + +# -- Project information ----------------------------------------------------- + +project = 'PyDFT' +copyright = '2021, Inorganic Materials and Catalysis' +author = 'Ivo Filot' + + +# -- General configuration --------------------------------------------------- + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'sphinx.ext.mathjax', + 'sphinx.ext.autosectionlabel', + 'sphinx_rtd_theme', + 'sphinx.ext.autodoc', + 'sphinx.ext.napoleon' +] + +suppress_warnings = ['autosectionlabel.*'] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'sphinx_rtd_theme' + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +master_doc = 'index' +html_static_path = ['_static'] +# html_theme_options = { +# 'display_version': True, +# 'analytics_id': 'G-H71EPP6GVB' +# } +html_logo = "_static/img/pydft_128px.png" +html_favicon = "_static/img/favicon.ico" +html_css_files = [ + "https://cdnjs.cloudflare.com/ajax/libs/font-awesome/6.1.1/css/all.min.css" +] + +# other options +html_show_sourcelink = False + +def setup(app): + app.add_css_file('css/custom.css') diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 0000000..90a1a01 --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,52 @@ +PyDFT: pure-python density functional theory +============================================ + +.. image:: https://anaconda.org/ifilot/pydft/badges/version.svg + :target: https://anaconda.org/ifilot/pydft +.. image:: https://img.shields.io/pypi/v/pydft?color=green + :target: https://pypi.org/project/pydft/ +.. image:: https://gitlab.tue.nl/ifilot/pydft/badges/master/pipeline.svg + :target: https://gitlab.tue.nl/ifilot/pydft/-/commits/master +.. image:: https://img.shields.io/badge/License-GPLv3-blue.svg + :target: https://www.gnu.org/licenses/gpl-3.0 + +:program:`PyDFT` is a pure-Python package for performing localized-orbital DFT +calculations using Gaussian Type Orbitals. + +.. image:: _static/img/pydft_logo_full_512px.png + +:program:`PyDFT` currently supports LDA and PBE +exchange-correlation functionals. The purpose of :program:`PyDFT` is mainly to +serve as an educational tool to explain the inner workings of a DFT +calculation. This program is not intended for professional calculations. It is +not particularly fast nor offers a lot of features that more mature +open-source of commercial packages offer. It does offer a unique insight into +a working code and a considerable effort was made in documenting everything. + +.. tip:: + + More information on the inner workings of :program:`PyDFT` can be obtained + from the textbook "Elements of Electronic Structure Theory" (specifically + chapter 4), which is freely available `via this website `_. + +:program:`PyDFT` has been developed at the Eindhoven University of Technology, +Netherlands. :program:`PyDFT` and its development are hosted on `github +`_. Bugs and feature +requests are ideally submitted via the `github issue tracker +`_. + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + installation + background + user_interface + api + community_guidelines + +Indices and tables +------------------ + +* :ref:`genindex` +* :ref:`search` diff --git a/docs/installation.rst b/docs/installation.rst new file mode 100644 index 0000000..4095700 --- /dev/null +++ b/docs/installation.rst @@ -0,0 +1,731 @@ +.. _installation: +.. index:: Installation + +Installation +============ + +:program:`PyDFT` is available via both Anaconda as well as PyPi. If you are using +Windows and have little experience with Python, we recommend to use Anaconda. +Anaconda is a complete Python package containing many modules and useful programs. +Anaconda can be obtained `via this link `_. + +.. note:: + + If you run into trouble installing :program:`PyDFT` in Anaconda, an easy + solution can be to create a separate environment for Anaconda. Please consult + the "Troubleshooting" section as seen below. + +Anaconda +-------- + +Open a Anaconda command prompt and run the following command: + +.. code:: bash + + conda install -c ifilot pydft pyqint pylebedev pytessel + conda install -c conda-forge mendeleev + +PyPi +---- + +Open a terminal and run the following command: + +.. code:: bash + + pip install pydft pyqint pylebedev pytessel mendeleev + +Testing +------- + +To test that your installation is working, you can run the following snippet +of code + +.. code:: python + + import pydft + print(pydft.__version__) + +The version number should be returned. + +Simple calculation +------------------ + +To perform a more involved calculation, one can run the following small +example code: + +.. code:: python + + from pydft import MoleculeBuilder,DFT + + co = MoleculeBuilder().get_molecule("CO") + dft = DFT(co, basis='sto3g') + en = dft.scf(1e-4) + print("Total electronic energy: %f Ht" % en) + +which should return the following result: + +.. code:: + + Total electronic energy: -111.147096 Ht + +Troubleshooting +--------------- + +The Anaconda packaging system can sometimes be quite finicky and +packages conflict with each other. A common scenario wherein this happens is +when you have installed a combination of Conda and PyPi packages. +A way to work around this issue is to create a separate environment and only +use that environment for the electronic resources associated with this project. + +.. note:: + + The result of each command is explicitly shown. Please do not let this + overwhelm you. You do not have to explicitly check every line. These + results are merely shown so that you know what to expect. + +To create the new environment (called :code:`pydft-env`), run:: + + conda create -n pydft-env python=3.11 + +which will give the following result:: + + ## Package Plan ## + + environment location: C:\Users\iawfi\anaconda3\envs\pydft-env + + added / updated specs: + - python=3.11 + + + The following packages will be downloaded: + + package | build + ---------------------------|----------------- + bzip2-1.0.8 | hcfcfb64_5 122 KB conda-forge + ca-certificates-2023.11.17 | h56e8100_0 151 KB conda-forge + libexpat-2.5.0 | h63175ca_1 135 KB conda-forge + libffi-3.4.2 | h8ffe710_5 41 KB conda-forge + libsqlite-3.44.2 | hcfcfb64_0 833 KB conda-forge + libzlib-1.2.13 | hcfcfb64_5 54 KB conda-forge + openssl-3.2.0 | hcfcfb64_0 7.8 MB conda-forge + pip-23.3.1 | pyhd8ed1ab_0 1.3 MB conda-forge + python-3.11.6 |h2628c8c_0_cpython 17.3 MB conda-forge + setuptools-68.2.2 | pyhd8ed1ab_0 454 KB conda-forge + tk-8.6.13 | h5226925_1 3.3 MB conda-forge + tzdata-2023c | h71feb2d_0 115 KB conda-forge + ucrt-10.0.22621.0 | h57928b3_0 1.2 MB conda-forge + vc-14.3 | h64f974e_17 17 KB conda-forge + vc14_runtime-14.36.32532 | hdcecf7f_17 722 KB conda-forge + vs2015_runtime-14.36.32532 | h05e6639_17 17 KB conda-forge + wheel-0.41.3 | pyhd8ed1ab_0 57 KB conda-forge + xz-5.2.6 | h8d14728_0 213 KB conda-forge + ------------------------------------------------------------ + Total: 33.8 MB + + The following NEW packages will be INSTALLED: + + bzip2 conda-forge/win-64::bzip2-1.0.8-hcfcfb64_5 + ca-certificates conda-forge/win-64::ca-certificates-2023.11.17-h56e8100_0 + libexpat conda-forge/win-64::libexpat-2.5.0-h63175ca_1 + libffi conda-forge/win-64::libffi-3.4.2-h8ffe710_5 + libsqlite conda-forge/win-64::libsqlite-3.44.2-hcfcfb64_0 + libzlib conda-forge/win-64::libzlib-1.2.13-hcfcfb64_5 + openssl conda-forge/win-64::openssl-3.2.0-hcfcfb64_0 + pip conda-forge/noarch::pip-23.3.1-pyhd8ed1ab_0 + python conda-forge/win-64::python-3.11.6-h2628c8c_0_cpython + setuptools conda-forge/noarch::setuptools-68.2.2-pyhd8ed1ab_0 + tk conda-forge/win-64::tk-8.6.13-h5226925_1 + tzdata conda-forge/noarch::tzdata-2023c-h71feb2d_0 + ucrt conda-forge/win-64::ucrt-10.0.22621.0-h57928b3_0 + vc conda-forge/win-64::vc-14.3-h64f974e_17 + vc14_runtime conda-forge/win-64::vc14_runtime-14.36.32532-hdcecf7f_17 + vs2015_runtime conda-forge/win-64::vs2015_runtime-14.36.32532-h05e6639_17 + wheel conda-forge/noarch::wheel-0.41.3-pyhd8ed1ab_0 + xz conda-forge/win-64::xz-5.2.6-h8d14728_0 + + + Proceed ([y]/n)? y + + + Downloading and Extracting Packages + + Preparing transaction: done + Verifying transaction: done + Executing transaction: done + # + # To activate this environment, use + # + # $ conda activate pydft-env + # + # To deactivate an active environment, use + # + # $ conda deactivate + +Next, we will active the environment:: + + conda activate pydft-env + +You will see that your command line now starts with :code:`(pydft-env)` instead +of :code:`(base)`. + +We can now install the required packages into the environment:: + + conda install -c ifilot pyqint pylebedev pydft pytessel + +You will see a response similar to the one as seen below:: + + ## Package Plan ## + + environment location: C:\Users\iawfi\anaconda3\envs\pydft-env + + added / updated specs: + - pydft + - pylebedev + - pyqint + - pytessel + + + The following packages will be downloaded: + + package | build + ---------------------------|----------------- + colorama-0.4.6 | pyhd8ed1ab_0 25 KB conda-forge + intel-openmp-2023.2.0 | h57928b3_50497 2.4 MB conda-forge + libblas-3.9.0 | 20_win64_mkl 4.8 MB conda-forge + libcblas-3.9.0 | 20_win64_mkl 4.8 MB conda-forge + libhwloc-2.9.3 |default_haede6df_1009 2.5 MB conda-forge + libiconv-1.17 | h8ffe710_0 698 KB conda-forge + liblapack-3.9.0 | 20_win64_mkl 4.8 MB conda-forge + libxml2-2.11.6 | hc3477c8_0 1.6 MB conda-forge + mkl-2023.2.0 | h6a75c08_50497 138.0 MB conda-forge + numpy-1.26.2 | py311h0b4df5a_0 6.8 MB conda-forge + pthreads-win32-2.9.1 | hfa6e2cd_3 141 KB conda-forge + pydft-0.2.4.1 | pyh4f56d60_0 52 KB ifilot + pyqint-0.14.0.1 | py311hcfd9ee6_0 286 KB ifilot + pytessel-1.1.0 | py311h1d48e73_0 55 KB ifilot + python_abi-3.11 | 4_cp311 7 KB conda-forge + scipy-1.11.4 | py311h0b4df5a_0 14.2 MB conda-forge + tbb-2021.10.0 | h91493d7_2 153 KB conda-forge + tqdm-4.66.1 | pyhd8ed1ab_0 87 KB conda-forge + ------------------------------------------------------------ + Total: 181.1 MB + + The following NEW packages will be INSTALLED: + + colorama conda-forge/noarch::colorama-0.4.6-pyhd8ed1ab_0 + intel-openmp conda-forge/win-64::intel-openmp-2023.2.0-h57928b3_50497 + libblas conda-forge/win-64::libblas-3.9.0-20_win64_mkl + libcblas conda-forge/win-64::libcblas-3.9.0-20_win64_mkl + libhwloc conda-forge/win-64::libhwloc-2.9.3-default_haede6df_1009 + libiconv conda-forge/win-64::libiconv-1.17-h8ffe710_0 + liblapack conda-forge/win-64::liblapack-3.9.0-20_win64_mkl + libxml2 conda-forge/win-64::libxml2-2.11.6-hc3477c8_0 + mkl conda-forge/win-64::mkl-2023.2.0-h6a75c08_50497 + numpy conda-forge/win-64::numpy-1.26.2-py311h0b4df5a_0 + pthreads-win32 conda-forge/win-64::pthreads-win32-2.9.1-hfa6e2cd_3 + pydft ifilot/noarch::pydft-0.2.4.1-pyh4f56d60_0 + pylebedev ifilot/noarch::pylebedev-1.0.0-pyh1d129d4_0 + pyqint ifilot/win-64::pyqint-0.14.0.1-py311hcfd9ee6_0 + pytessel ifilot/win-64::pytessel-1.1.0-py311h1d48e73_0 + python_abi conda-forge/win-64::python_abi-3.11-4_cp311 + scipy conda-forge/win-64::scipy-1.11.4-py311h0b4df5a_0 + tbb conda-forge/win-64::tbb-2021.10.0-h91493d7_2 + tqdm conda-forge/noarch::tqdm-4.66.1-pyhd8ed1ab_0 + + + Proceed ([y]/n)? y + + + Downloading and Extracting Packages + + Preparing transaction: done + Verifying transaction: done + Executing transaction: done + +Finally, you can install the IDE Spyder using:: + + conda install spyder matplotlib scipy pandas openpyxl mendeleev + +This might take a while (the environment needs to resolve all dependencies), +but you should see a response similar to the one below:: + + ## Package Plan ## + + environment location: C:\Users\iawfi\anaconda3\envs\pydft-env + + added / updated specs: + - matplotlib + - openpyxl + - pandas + - scipy + - spyder + + + The following packages will be downloaded: + + package | build + ---------------------------|----------------- + alabaster-0.7.13 | pyhd8ed1ab_0 18 KB conda-forge + arrow-1.3.0 | pyhd8ed1ab_0 98 KB conda-forge + astroid-3.0.1 | py311h1ea47a8_0 498 KB conda-forge + asttokens-2.4.1 | pyhd8ed1ab_0 28 KB conda-forge + atomicwrites-1.4.1 | pyhd8ed1ab_0 12 KB conda-forge + attrs-23.1.0 | pyh71513ae_1 54 KB conda-forge + autopep8-2.0.4 | pyhd8ed1ab_0 45 KB conda-forge + babel-2.13.1 | pyhd8ed1ab_0 6.6 MB conda-forge + bcrypt-4.0.1 | py311hc37eb10_1 144 KB conda-forge + beautifulsoup4-4.12.2 | pyha770c72_0 112 KB conda-forge + binaryornot-0.4.4 | py_1 370 KB conda-forge + black-23.10.1 | py311h1ea47a8_0 369 KB conda-forge + bleach-6.1.0 | pyhd8ed1ab_0 128 KB conda-forge + brotli-1.1.0 | hcfcfb64_1 19 KB conda-forge + brotli-bin-1.1.0 | hcfcfb64_1 20 KB conda-forge + brotli-python-1.1.0 | py311h12c1d0e_1 315 KB conda-forge + certifi-2023.11.17 | pyhd8ed1ab_0 155 KB conda-forge + cffi-1.16.0 | py311ha68e1ae_0 290 KB conda-forge + chardet-5.2.0 | py311h1ea47a8_1 278 KB conda-forge + charset-normalizer-3.3.2 | pyhd8ed1ab_0 46 KB conda-forge + click-8.1.7 | win_pyh7428d3b_0 83 KB conda-forge + cloudpickle-3.0.0 | pyhd8ed1ab_0 24 KB conda-forge + comm-0.1.4 | pyhd8ed1ab_0 11 KB conda-forge + contourpy-1.2.0 | py311h005e61a_0 201 KB conda-forge + cookiecutter-2.5.0 | pyhca7485f_0 97 KB conda-forge + cryptography-41.0.5 | py311h28e9c30_0 1.1 MB conda-forge + cycler-0.12.1 | pyhd8ed1ab_0 13 KB conda-forge + debugpy-1.8.0 | py311h12c1d0e_1 3.7 MB conda-forge + decorator-5.1.1 | pyhd8ed1ab_0 12 KB conda-forge + defusedxml-0.7.1 | pyhd8ed1ab_0 23 KB conda-forge + diff-match-patch-20230430 | pyhd8ed1ab_0 40 KB conda-forge + dill-0.3.7 | pyhd8ed1ab_0 86 KB conda-forge + docstring-to-markdown-0.13 | pyhd8ed1ab_0 31 KB conda-forge + docutils-0.20.1 | py311h1ea47a8_2 950 KB conda-forge + entrypoints-0.4 | pyhd8ed1ab_0 9 KB conda-forge + et_xmlfile-1.1.0 | pyhd8ed1ab_0 10 KB conda-forge + exceptiongroup-1.2.0 | pyhd8ed1ab_0 20 KB conda-forge + executing-2.0.1 | pyhd8ed1ab_0 27 KB conda-forge + flake8-6.1.0 | pyhd8ed1ab_0 109 KB conda-forge + fonttools-4.45.1 | py311ha68e1ae_0 2.3 MB conda-forge + freetype-2.12.1 | hdaf720e_2 498 KB conda-forge + gettext-0.21.1 | h5728263_0 5.3 MB conda-forge + glib-2.78.1 | h12be248_1 495 KB conda-forge + glib-tools-2.78.1 | h12be248_1 141 KB conda-forge + gst-plugins-base-1.22.7 | h001b923_0 1.9 MB conda-forge + gstreamer-1.22.7 | hb4038d2_0 1.8 MB conda-forge + icu-72.1 | h63175ca_0 12.6 MB conda-forge + idna-3.5 | pyhd8ed1ab_0 48 KB conda-forge + imagesize-1.4.1 | pyhd8ed1ab_0 10 KB conda-forge + importlib-metadata-6.8.0 | pyha770c72_0 25 KB conda-forge + importlib_metadata-6.8.0 | hd8ed1ab_0 9 KB conda-forge + importlib_resources-6.1.1 | pyhd8ed1ab_0 29 KB conda-forge + inflection-0.5.1 | pyh9f0ad1d_0 9 KB conda-forge + intervaltree-3.1.0 | pyhd8ed1ab_1 27 KB conda-forge + ipykernel-6.26.0 | pyha63f2e9_0 114 KB conda-forge + ipython-8.18.0 | pyh5737063_0 577 KB conda-forge + isort-5.12.0 | pyhd8ed1ab_1 72 KB conda-forge + jaraco.classes-3.3.0 | pyhd8ed1ab_0 11 KB conda-forge + jedi-0.19.1 | pyhd8ed1ab_0 822 KB conda-forge + jellyfish-1.0.3 | py311hc37eb10_0 191 KB conda-forge + jinja2-3.1.2 | pyhd8ed1ab_1 99 KB conda-forge + jsonschema-4.20.0 | pyhd8ed1ab_0 70 KB conda-forge + jsonschema-specifications-2023.11.1| pyhd8ed1ab_0 15 KB conda-forge + jupyter_client-8.6.0 | pyhd8ed1ab_0 103 KB conda-forge + jupyter_core-5.5.0 | py311h1ea47a8_0 108 KB conda-forge + jupyterlab_pygments-0.3.0 | pyhd8ed1ab_0 18 KB conda-forge + keyring-24.3.0 | py311h1ea47a8_0 92 KB conda-forge + kiwisolver-1.4.5 | py311h005e61a_1 55 KB conda-forge + krb5-1.20.1 | heb0366b_0 701 KB conda-forge + lcms2-2.15 | he9d350c_2 486 KB conda-forge + lerc-4.0.0 | h63175ca_0 190 KB conda-forge + libbrotlicommon-1.1.0 | hcfcfb64_1 69 KB conda-forge + libbrotlidec-1.1.0 | hcfcfb64_1 32 KB conda-forge + libbrotlienc-1.1.0 | hcfcfb64_1 241 KB conda-forge + libclang-16.0.6 |default_heb8d277_2 35 KB conda-forge + libclang13-16.0.6 |default_hc80b9e7_2 22.1 MB conda-forge + libdeflate-1.19 | hcfcfb64_0 150 KB conda-forge + libglib-2.78.1 | h16e383f_1 2.5 MB conda-forge + libjpeg-turbo-2.1.5.1 | hcfcfb64_1 672 KB conda-forge + libogg-1.3.4 | h8ffe710_1 34 KB conda-forge + libpng-1.6.39 | h19919ed_0 336 KB conda-forge + libsodium-1.0.18 | h8d14728_1 697 KB conda-forge + libspatialindex-1.9.3 | h39d44d4_4 437 KB conda-forge + libtiff-4.6.0 | h4554b19_1 766 KB conda-forge + libvorbis-1.3.7 | h0e60522_0 267 KB conda-forge + libwebp-1.3.2 | hcfcfb64_1 69 KB conda-forge + libwebp-base-1.3.2 | hcfcfb64_0 263 KB conda-forge + libxcb-1.15 | hcd874cb_0 947 KB conda-forge + m2w64-gcc-libgfortran-5.3.0| 6 342 KB conda-forge + m2w64-gcc-libs-5.3.0 | 7 520 KB conda-forge + m2w64-gcc-libs-core-5.3.0 | 7 214 KB conda-forge + m2w64-gmp-6.1.0 | 2 726 KB conda-forge + m2w64-libwinpthread-git-5.0.0.4634.697f757| 2 31 KB conda-forge + markdown-it-py-3.0.0 | pyhd8ed1ab_0 63 KB conda-forge + markupsafe-2.1.3 | py311ha68e1ae_1 29 KB conda-forge + matplotlib-3.8.2 | py311h1ea47a8_0 9 KB conda-forge + matplotlib-base-3.8.2 | py311h6e989c2_0 7.3 MB conda-forge + matplotlib-inline-0.1.6 | pyhd8ed1ab_0 12 KB conda-forge + mccabe-0.7.0 | pyhd8ed1ab_0 11 KB conda-forge + mdurl-0.1.0 | pyhd8ed1ab_0 13 KB conda-forge + mistune-3.0.2 | pyhd8ed1ab_0 64 KB conda-forge + more-itertools-10.1.0 | pyhd8ed1ab_0 52 KB conda-forge + msys2-conda-epoch-20160418 | 1 3 KB conda-forge + munkres-1.1.4 | pyh9f0ad1d_0 12 KB conda-forge + mypy_extensions-1.0.0 | pyha770c72_0 10 KB conda-forge + nbclient-0.8.0 | pyhd8ed1ab_0 63 KB conda-forge + nbconvert-7.11.0 | pyhd8ed1ab_0 8 KB conda-forge + nbconvert-core-7.11.0 | pyhd8ed1ab_0 183 KB conda-forge + nbconvert-pandoc-7.11.0 | pyhd8ed1ab_0 7 KB conda-forge + nbformat-5.9.2 | pyhd8ed1ab_0 98 KB conda-forge + nest-asyncio-1.5.8 | pyhd8ed1ab_0 11 KB conda-forge + numpydoc-1.5.0 | pyhd8ed1ab_0 46 KB conda-forge + openjpeg-2.5.0 | h3d672ee_3 231 KB conda-forge + openpyxl-3.1.2 | py311ha68e1ae_1 635 KB conda-forge + packaging-23.2 | pyhd8ed1ab_0 48 KB conda-forge + pandas-2.1.3 | py311hf63dbb6_0 13.2 MB conda-forge + pandoc-3.1.3 | h57928b3_0 17.8 MB conda-forge + pandocfilters-1.5.0 | pyhd8ed1ab_0 11 KB conda-forge + paramiko-3.3.1 | pyhd8ed1ab_0 155 KB conda-forge + parso-0.8.3 | pyhd8ed1ab_0 69 KB conda-forge + pathspec-0.11.2 | pyhd8ed1ab_0 38 KB conda-forge + pcre2-10.42 | h17e33f8_0 860 KB conda-forge + pexpect-4.8.0 | pyh1a96a4e_2 48 KB conda-forge + pickleshare-0.7.5 | py_1003 9 KB conda-forge + pillow-10.0.1 | py311hd926f49_1 44.7 MB conda-forge + pkgutil-resolve-name-1.3.10| pyhd8ed1ab_1 11 KB conda-forge + platformdirs-4.0.0 | pyhd8ed1ab_0 19 KB conda-forge + pluggy-1.3.0 | pyhd8ed1ab_0 22 KB conda-forge + ply-3.11 | py_1 44 KB conda-forge + prompt-toolkit-3.0.41 | pyha770c72_0 264 KB conda-forge + prompt_toolkit-3.0.41 | hd8ed1ab_0 7 KB conda-forge + psutil-5.9.5 | py311ha68e1ae_1 504 KB conda-forge + pthread-stubs-0.4 | hcd874cb_1001 6 KB conda-forge + ptyprocess-0.7.0 | pyhd3deb0d_0 16 KB conda-forge + pure_eval-0.2.2 | pyhd8ed1ab_0 14 KB conda-forge + pycodestyle-2.11.1 | pyhd8ed1ab_0 34 KB conda-forge + pycparser-2.21 | pyhd8ed1ab_0 100 KB conda-forge + pydocstyle-6.3.0 | pyhd8ed1ab_0 39 KB conda-forge + pyflakes-3.1.0 | pyhd8ed1ab_0 57 KB conda-forge + pygments-2.17.2 | pyhd8ed1ab_0 840 KB conda-forge + pylint-3.0.2 | pyhd8ed1ab_0 338 KB conda-forge + pylint-venv-3.0.3 | pyhd8ed1ab_0 11 KB conda-forge + pyls-spyder-0.4.0 | pyhd8ed1ab_0 10 KB conda-forge + pynacl-1.5.0 | py311hd53affc_3 1.2 MB conda-forge + pyparsing-3.1.1 | pyhd8ed1ab_0 87 KB conda-forge + pyqt-5.15.9 | py311h125bc19_5 3.7 MB conda-forge + pyqt5-sip-12.12.2 | py311h12c1d0e_5 78 KB conda-forge + pyqtwebengine-5.15.9 | py311h5a77453_5 123 KB conda-forge + pysocks-1.7.1 | pyh0701188_6 19 KB conda-forge + python-dateutil-2.8.2 | pyhd8ed1ab_0 240 KB conda-forge + python-fastjsonschema-2.19.0| pyhd8ed1ab_0 221 KB conda-forge + python-lsp-black-1.3.0 | pyhd8ed1ab_0 12 KB conda-forge + python-lsp-jsonrpc-1.1.2 | pyhd8ed1ab_0 14 KB conda-forge + python-lsp-server-1.9.0 | pyhd8ed1ab_0 7 KB conda-forge + python-lsp-server-base-1.9.0| pyhd8ed1ab_0 60 KB conda-forge + python-slugify-8.0.1 | pyhd8ed1ab_2 15 KB conda-forge + python-tzdata-2023.3 | pyhd8ed1ab_0 140 KB conda-forge + pytoolconfig-1.2.5 | pyhd8ed1ab_0 21 KB conda-forge + pytz-2023.3.post1 | pyhd8ed1ab_0 183 KB conda-forge + pywin32-306 | py311h12c1d0e_2 5.8 MB conda-forge + pywin32-ctypes-0.2.2 | py311h1ea47a8_1 56 KB conda-forge + pyyaml-6.0.1 | py311ha68e1ae_1 171 KB conda-forge + pyzmq-25.1.1 | py311h9250fbb_2 480 KB conda-forge + qdarkstyle-3.2 | pyhd8ed1ab_0 612 KB conda-forge + qstylizer-0.2.2 | pyhd8ed1ab_0 17 KB conda-forge + qt-main-5.15.8 | h2c8576c_12 56.9 MB conda-forge + qt-webengine-5.15.8 | h5b1ea0b_0 62.1 MB conda-forge + qtawesome-1.2.3 | pyhd8ed1ab_0 1.4 MB conda-forge + qtconsole-5.5.1 | pyhd8ed1ab_0 7 KB conda-forge + qtconsole-base-5.5.1 | pyha770c72_0 98 KB conda-forge + qtpy-2.4.1 | pyhd8ed1ab_0 60 KB conda-forge + referencing-0.31.0 | pyhd8ed1ab_0 37 KB conda-forge + requests-2.31.0 | pyhd8ed1ab_0 55 KB conda-forge + rich-13.7.0 | pyhd8ed1ab_0 180 KB conda-forge + rope-1.11.0 | pyhd8ed1ab_1 145 KB conda-forge + rpds-py-0.13.1 | py311hc37eb10_0 178 KB conda-forge + rtree-1.1.0 | py311hcacb13a_0 62 KB conda-forge + sip-6.7.12 | py311h12c1d0e_0 581 KB conda-forge + six-1.16.0 | pyh6c4a22f_0 14 KB conda-forge + snowballstemmer-2.2.0 | pyhd8ed1ab_0 57 KB conda-forge + sortedcontainers-2.4.0 | pyhd8ed1ab_0 26 KB conda-forge + soupsieve-2.5 | pyhd8ed1ab_1 36 KB conda-forge + sphinx-7.2.6 | pyhd8ed1ab_0 1.2 MB conda-forge + sphinxcontrib-applehelp-1.0.7| pyhd8ed1ab_0 29 KB conda-forge + sphinxcontrib-devhelp-1.0.5| pyhd8ed1ab_0 24 KB conda-forge + sphinxcontrib-htmlhelp-2.0.4| pyhd8ed1ab_0 32 KB conda-forge + sphinxcontrib-jsmath-1.0.1 | pyhd8ed1ab_0 10 KB conda-forge + sphinxcontrib-qthelp-1.0.6 | pyhd8ed1ab_0 26 KB conda-forge + sphinxcontrib-serializinghtml-1.1.9| pyhd8ed1ab_0 28 KB conda-forge + spyder-5.5.0 | py311h1ea47a8_3 11.2 MB conda-forge + spyder-kernels-2.5.0 | win_pyh7428d3b_0 80 KB conda-forge + stack_data-0.6.2 | pyhd8ed1ab_0 26 KB conda-forge + text-unidecode-1.3 | pyhd8ed1ab_1 64 KB conda-forge + textdistance-4.5.0 | pyhd8ed1ab_0 28 KB conda-forge + three-merge-0.1.1 | pyh9f0ad1d_0 8 KB conda-forge + tinycss2-1.2.1 | pyhd8ed1ab_0 23 KB conda-forge + toml-0.10.2 | pyhd8ed1ab_0 18 KB conda-forge + tomli-2.0.1 | pyhd8ed1ab_0 16 KB conda-forge + tomlkit-0.12.3 | pyha770c72_0 36 KB conda-forge + tornado-6.3.3 | py311ha68e1ae_1 826 KB conda-forge + traitlets-5.13.0 | pyhd8ed1ab_0 107 KB conda-forge + types-python-dateutil-2.8.19.14| pyhd8ed1ab_0 21 KB conda-forge + typing-extensions-4.8.0 | hd8ed1ab_0 10 KB conda-forge + typing_extensions-4.8.0 | pyha770c72_0 34 KB conda-forge + ujson-5.8.0 | py311h12c1d0e_0 47 KB conda-forge + urllib3-2.1.0 | pyhd8ed1ab_0 83 KB conda-forge + watchdog-3.0.0 | py311h1ea47a8_1 151 KB conda-forge + wcwidth-0.2.12 | pyhd8ed1ab_0 32 KB conda-forge + webencodings-0.5.1 | pyhd8ed1ab_2 15 KB conda-forge + whatthepatch-1.0.5 | pyhd8ed1ab_0 17 KB conda-forge + win_inet_pton-1.1.0 | pyhd8ed1ab_6 8 KB conda-forge + xorg-libxau-1.0.11 | hcd874cb_0 50 KB conda-forge + xorg-libxdmcp-1.1.3 | hcd874cb_0 66 KB conda-forge + yaml-0.2.5 | h8ffe710_2 62 KB conda-forge + yapf-0.40.1 | pyhd8ed1ab_0 172 KB conda-forge + zeromq-4.3.5 | h63175ca_0 4.0 MB conda-forge + zipp-3.17.0 | pyhd8ed1ab_0 19 KB conda-forge + zstd-1.5.5 | h12be248_0 335 KB conda-forge + ------------------------------------------------------------ + Total: 318.3 MB + + The following NEW packages will be INSTALLED: + + alabaster conda-forge/noarch::alabaster-0.7.13-pyhd8ed1ab_0 + arrow conda-forge/noarch::arrow-1.3.0-pyhd8ed1ab_0 + astroid conda-forge/win-64::astroid-3.0.1-py311h1ea47a8_0 + asttokens conda-forge/noarch::asttokens-2.4.1-pyhd8ed1ab_0 + atomicwrites conda-forge/noarch::atomicwrites-1.4.1-pyhd8ed1ab_0 + attrs conda-forge/noarch::attrs-23.1.0-pyh71513ae_1 + autopep8 conda-forge/noarch::autopep8-2.0.4-pyhd8ed1ab_0 + babel conda-forge/noarch::babel-2.13.1-pyhd8ed1ab_0 + bcrypt conda-forge/win-64::bcrypt-4.0.1-py311hc37eb10_1 + beautifulsoup4 conda-forge/noarch::beautifulsoup4-4.12.2-pyha770c72_0 + binaryornot conda-forge/noarch::binaryornot-0.4.4-py_1 + black conda-forge/win-64::black-23.10.1-py311h1ea47a8_0 + bleach conda-forge/noarch::bleach-6.1.0-pyhd8ed1ab_0 + brotli conda-forge/win-64::brotli-1.1.0-hcfcfb64_1 + brotli-bin conda-forge/win-64::brotli-bin-1.1.0-hcfcfb64_1 + brotli-python conda-forge/win-64::brotli-python-1.1.0-py311h12c1d0e_1 + certifi conda-forge/noarch::certifi-2023.11.17-pyhd8ed1ab_0 + cffi conda-forge/win-64::cffi-1.16.0-py311ha68e1ae_0 + chardet conda-forge/win-64::chardet-5.2.0-py311h1ea47a8_1 + charset-normalizer conda-forge/noarch::charset-normalizer-3.3.2-pyhd8ed1ab_0 + click conda-forge/noarch::click-8.1.7-win_pyh7428d3b_0 + cloudpickle conda-forge/noarch::cloudpickle-3.0.0-pyhd8ed1ab_0 + comm conda-forge/noarch::comm-0.1.4-pyhd8ed1ab_0 + contourpy conda-forge/win-64::contourpy-1.2.0-py311h005e61a_0 + cookiecutter conda-forge/noarch::cookiecutter-2.5.0-pyhca7485f_0 + cryptography conda-forge/win-64::cryptography-41.0.5-py311h28e9c30_0 + cycler conda-forge/noarch::cycler-0.12.1-pyhd8ed1ab_0 + debugpy conda-forge/win-64::debugpy-1.8.0-py311h12c1d0e_1 + decorator conda-forge/noarch::decorator-5.1.1-pyhd8ed1ab_0 + defusedxml conda-forge/noarch::defusedxml-0.7.1-pyhd8ed1ab_0 + diff-match-patch conda-forge/noarch::diff-match-patch-20230430-pyhd8ed1ab_0 + dill conda-forge/noarch::dill-0.3.7-pyhd8ed1ab_0 + docstring-to-mark~ conda-forge/noarch::docstring-to-markdown-0.13-pyhd8ed1ab_0 + docutils conda-forge/win-64::docutils-0.20.1-py311h1ea47a8_2 + entrypoints conda-forge/noarch::entrypoints-0.4-pyhd8ed1ab_0 + et_xmlfile conda-forge/noarch::et_xmlfile-1.1.0-pyhd8ed1ab_0 + exceptiongroup conda-forge/noarch::exceptiongroup-1.2.0-pyhd8ed1ab_0 + executing conda-forge/noarch::executing-2.0.1-pyhd8ed1ab_0 + flake8 conda-forge/noarch::flake8-6.1.0-pyhd8ed1ab_0 + fonttools conda-forge/win-64::fonttools-4.45.1-py311ha68e1ae_0 + freetype conda-forge/win-64::freetype-2.12.1-hdaf720e_2 + gettext conda-forge/win-64::gettext-0.21.1-h5728263_0 + glib conda-forge/win-64::glib-2.78.1-h12be248_1 + glib-tools conda-forge/win-64::glib-tools-2.78.1-h12be248_1 + gst-plugins-base conda-forge/win-64::gst-plugins-base-1.22.7-h001b923_0 + gstreamer conda-forge/win-64::gstreamer-1.22.7-hb4038d2_0 + icu conda-forge/win-64::icu-72.1-h63175ca_0 + idna conda-forge/noarch::idna-3.5-pyhd8ed1ab_0 + imagesize conda-forge/noarch::imagesize-1.4.1-pyhd8ed1ab_0 + importlib-metadata conda-forge/noarch::importlib-metadata-6.8.0-pyha770c72_0 + importlib_metadata conda-forge/noarch::importlib_metadata-6.8.0-hd8ed1ab_0 + importlib_resourc~ conda-forge/noarch::importlib_resources-6.1.1-pyhd8ed1ab_0 + inflection conda-forge/noarch::inflection-0.5.1-pyh9f0ad1d_0 + intervaltree conda-forge/noarch::intervaltree-3.1.0-pyhd8ed1ab_1 + ipykernel conda-forge/noarch::ipykernel-6.26.0-pyha63f2e9_0 + ipython conda-forge/noarch::ipython-8.18.0-pyh5737063_0 + isort conda-forge/noarch::isort-5.12.0-pyhd8ed1ab_1 + jaraco.classes conda-forge/noarch::jaraco.classes-3.3.0-pyhd8ed1ab_0 + jedi conda-forge/noarch::jedi-0.19.1-pyhd8ed1ab_0 + jellyfish conda-forge/win-64::jellyfish-1.0.3-py311hc37eb10_0 + jinja2 conda-forge/noarch::jinja2-3.1.2-pyhd8ed1ab_1 + jsonschema conda-forge/noarch::jsonschema-4.20.0-pyhd8ed1ab_0 + jsonschema-specif~ conda-forge/noarch::jsonschema-specifications-2023.11.1-pyhd8ed1ab_0 + jupyter_client conda-forge/noarch::jupyter_client-8.6.0-pyhd8ed1ab_0 + jupyter_core conda-forge/win-64::jupyter_core-5.5.0-py311h1ea47a8_0 + jupyterlab_pygmen~ conda-forge/noarch::jupyterlab_pygments-0.3.0-pyhd8ed1ab_0 + keyring conda-forge/win-64::keyring-24.3.0-py311h1ea47a8_0 + kiwisolver conda-forge/win-64::kiwisolver-1.4.5-py311h005e61a_1 + krb5 conda-forge/win-64::krb5-1.20.1-heb0366b_0 + lcms2 conda-forge/win-64::lcms2-2.15-he9d350c_2 + lerc conda-forge/win-64::lerc-4.0.0-h63175ca_0 + libbrotlicommon conda-forge/win-64::libbrotlicommon-1.1.0-hcfcfb64_1 + libbrotlidec conda-forge/win-64::libbrotlidec-1.1.0-hcfcfb64_1 + libbrotlienc conda-forge/win-64::libbrotlienc-1.1.0-hcfcfb64_1 + libclang conda-forge/win-64::libclang-16.0.6-default_heb8d277_2 + libclang13 conda-forge/win-64::libclang13-16.0.6-default_hc80b9e7_2 + libdeflate conda-forge/win-64::libdeflate-1.19-hcfcfb64_0 + libglib conda-forge/win-64::libglib-2.78.1-h16e383f_1 + libjpeg-turbo conda-forge/win-64::libjpeg-turbo-2.1.5.1-hcfcfb64_1 + libogg conda-forge/win-64::libogg-1.3.4-h8ffe710_1 + libpng conda-forge/win-64::libpng-1.6.39-h19919ed_0 + libsodium conda-forge/win-64::libsodium-1.0.18-h8d14728_1 + libspatialindex conda-forge/win-64::libspatialindex-1.9.3-h39d44d4_4 + libtiff conda-forge/win-64::libtiff-4.6.0-h4554b19_1 + libvorbis conda-forge/win-64::libvorbis-1.3.7-h0e60522_0 + libwebp conda-forge/win-64::libwebp-1.3.2-hcfcfb64_1 + libwebp-base conda-forge/win-64::libwebp-base-1.3.2-hcfcfb64_0 + libxcb conda-forge/win-64::libxcb-1.15-hcd874cb_0 + m2w64-gcc-libgfor~ conda-forge/win-64::m2w64-gcc-libgfortran-5.3.0-6 + m2w64-gcc-libs conda-forge/win-64::m2w64-gcc-libs-5.3.0-7 + m2w64-gcc-libs-co~ conda-forge/win-64::m2w64-gcc-libs-core-5.3.0-7 + m2w64-gmp conda-forge/win-64::m2w64-gmp-6.1.0-2 + m2w64-libwinpthre~ conda-forge/win-64::m2w64-libwinpthread-git-5.0.0.4634.697f757-2 + markdown-it-py conda-forge/noarch::markdown-it-py-3.0.0-pyhd8ed1ab_0 + markupsafe conda-forge/win-64::markupsafe-2.1.3-py311ha68e1ae_1 + matplotlib conda-forge/win-64::matplotlib-3.8.2-py311h1ea47a8_0 + matplotlib-base conda-forge/win-64::matplotlib-base-3.8.2-py311h6e989c2_0 + matplotlib-inline conda-forge/noarch::matplotlib-inline-0.1.6-pyhd8ed1ab_0 + mccabe conda-forge/noarch::mccabe-0.7.0-pyhd8ed1ab_0 + mdurl conda-forge/noarch::mdurl-0.1.0-pyhd8ed1ab_0 + mistune conda-forge/noarch::mistune-3.0.2-pyhd8ed1ab_0 + more-itertools conda-forge/noarch::more-itertools-10.1.0-pyhd8ed1ab_0 + msys2-conda-epoch conda-forge/win-64::msys2-conda-epoch-20160418-1 + munkres conda-forge/noarch::munkres-1.1.4-pyh9f0ad1d_0 + mypy_extensions conda-forge/noarch::mypy_extensions-1.0.0-pyha770c72_0 + nbclient conda-forge/noarch::nbclient-0.8.0-pyhd8ed1ab_0 + nbconvert conda-forge/noarch::nbconvert-7.11.0-pyhd8ed1ab_0 + nbconvert-core conda-forge/noarch::nbconvert-core-7.11.0-pyhd8ed1ab_0 + nbconvert-pandoc conda-forge/noarch::nbconvert-pandoc-7.11.0-pyhd8ed1ab_0 + nbformat conda-forge/noarch::nbformat-5.9.2-pyhd8ed1ab_0 + nest-asyncio conda-forge/noarch::nest-asyncio-1.5.8-pyhd8ed1ab_0 + numpydoc conda-forge/noarch::numpydoc-1.5.0-pyhd8ed1ab_0 + openjpeg conda-forge/win-64::openjpeg-2.5.0-h3d672ee_3 + openpyxl conda-forge/win-64::openpyxl-3.1.2-py311ha68e1ae_1 + packaging conda-forge/noarch::packaging-23.2-pyhd8ed1ab_0 + pandas conda-forge/win-64::pandas-2.1.3-py311hf63dbb6_0 + pandoc conda-forge/win-64::pandoc-3.1.3-h57928b3_0 + pandocfilters conda-forge/noarch::pandocfilters-1.5.0-pyhd8ed1ab_0 + paramiko conda-forge/noarch::paramiko-3.3.1-pyhd8ed1ab_0 + parso conda-forge/noarch::parso-0.8.3-pyhd8ed1ab_0 + pathspec conda-forge/noarch::pathspec-0.11.2-pyhd8ed1ab_0 + pcre2 conda-forge/win-64::pcre2-10.42-h17e33f8_0 + pexpect conda-forge/noarch::pexpect-4.8.0-pyh1a96a4e_2 + pickleshare conda-forge/noarch::pickleshare-0.7.5-py_1003 + pillow conda-forge/win-64::pillow-10.0.1-py311hd926f49_1 + pkgutil-resolve-n~ conda-forge/noarch::pkgutil-resolve-name-1.3.10-pyhd8ed1ab_1 + platformdirs conda-forge/noarch::platformdirs-4.0.0-pyhd8ed1ab_0 + pluggy conda-forge/noarch::pluggy-1.3.0-pyhd8ed1ab_0 + ply conda-forge/noarch::ply-3.11-py_1 + prompt-toolkit conda-forge/noarch::prompt-toolkit-3.0.41-pyha770c72_0 + prompt_toolkit conda-forge/noarch::prompt_toolkit-3.0.41-hd8ed1ab_0 + psutil conda-forge/win-64::psutil-5.9.5-py311ha68e1ae_1 + pthread-stubs conda-forge/win-64::pthread-stubs-0.4-hcd874cb_1001 + ptyprocess conda-forge/noarch::ptyprocess-0.7.0-pyhd3deb0d_0 + pure_eval conda-forge/noarch::pure_eval-0.2.2-pyhd8ed1ab_0 + pycodestyle conda-forge/noarch::pycodestyle-2.11.1-pyhd8ed1ab_0 + pycparser conda-forge/noarch::pycparser-2.21-pyhd8ed1ab_0 + pydocstyle conda-forge/noarch::pydocstyle-6.3.0-pyhd8ed1ab_0 + pyflakes conda-forge/noarch::pyflakes-3.1.0-pyhd8ed1ab_0 + pygments conda-forge/noarch::pygments-2.17.2-pyhd8ed1ab_0 + pylint conda-forge/noarch::pylint-3.0.2-pyhd8ed1ab_0 + pylint-venv conda-forge/noarch::pylint-venv-3.0.3-pyhd8ed1ab_0 + pyls-spyder conda-forge/noarch::pyls-spyder-0.4.0-pyhd8ed1ab_0 + pynacl conda-forge/win-64::pynacl-1.5.0-py311hd53affc_3 + pyparsing conda-forge/noarch::pyparsing-3.1.1-pyhd8ed1ab_0 + pyqt conda-forge/win-64::pyqt-5.15.9-py311h125bc19_5 + pyqt5-sip conda-forge/win-64::pyqt5-sip-12.12.2-py311h12c1d0e_5 + pyqtwebengine conda-forge/win-64::pyqtwebengine-5.15.9-py311h5a77453_5 + pysocks conda-forge/noarch::pysocks-1.7.1-pyh0701188_6 + python-dateutil conda-forge/noarch::python-dateutil-2.8.2-pyhd8ed1ab_0 + python-fastjsonsc~ conda-forge/noarch::python-fastjsonschema-2.19.0-pyhd8ed1ab_0 + python-lsp-black conda-forge/noarch::python-lsp-black-1.3.0-pyhd8ed1ab_0 + python-lsp-jsonrpc conda-forge/noarch::python-lsp-jsonrpc-1.1.2-pyhd8ed1ab_0 + python-lsp-server conda-forge/noarch::python-lsp-server-1.9.0-pyhd8ed1ab_0 + python-lsp-server~ conda-forge/noarch::python-lsp-server-base-1.9.0-pyhd8ed1ab_0 + python-slugify conda-forge/noarch::python-slugify-8.0.1-pyhd8ed1ab_2 + python-tzdata conda-forge/noarch::python-tzdata-2023.3-pyhd8ed1ab_0 + pytoolconfig conda-forge/noarch::pytoolconfig-1.2.5-pyhd8ed1ab_0 + pytz conda-forge/noarch::pytz-2023.3.post1-pyhd8ed1ab_0 + pywin32 conda-forge/win-64::pywin32-306-py311h12c1d0e_2 + pywin32-ctypes conda-forge/win-64::pywin32-ctypes-0.2.2-py311h1ea47a8_1 + pyyaml conda-forge/win-64::pyyaml-6.0.1-py311ha68e1ae_1 + pyzmq conda-forge/win-64::pyzmq-25.1.1-py311h9250fbb_2 + qdarkstyle conda-forge/noarch::qdarkstyle-3.2-pyhd8ed1ab_0 + qstylizer conda-forge/noarch::qstylizer-0.2.2-pyhd8ed1ab_0 + qt-main conda-forge/win-64::qt-main-5.15.8-h2c8576c_12 + qt-webengine conda-forge/win-64::qt-webengine-5.15.8-h5b1ea0b_0 + qtawesome conda-forge/noarch::qtawesome-1.2.3-pyhd8ed1ab_0 + qtconsole conda-forge/noarch::qtconsole-5.5.1-pyhd8ed1ab_0 + qtconsole-base conda-forge/noarch::qtconsole-base-5.5.1-pyha770c72_0 + qtpy conda-forge/noarch::qtpy-2.4.1-pyhd8ed1ab_0 + referencing conda-forge/noarch::referencing-0.31.0-pyhd8ed1ab_0 + requests conda-forge/noarch::requests-2.31.0-pyhd8ed1ab_0 + rich conda-forge/noarch::rich-13.7.0-pyhd8ed1ab_0 + rope conda-forge/noarch::rope-1.11.0-pyhd8ed1ab_1 + rpds-py conda-forge/win-64::rpds-py-0.13.1-py311hc37eb10_0 + rtree conda-forge/win-64::rtree-1.1.0-py311hcacb13a_0 + sip conda-forge/win-64::sip-6.7.12-py311h12c1d0e_0 + six conda-forge/noarch::six-1.16.0-pyh6c4a22f_0 + snowballstemmer conda-forge/noarch::snowballstemmer-2.2.0-pyhd8ed1ab_0 + sortedcontainers conda-forge/noarch::sortedcontainers-2.4.0-pyhd8ed1ab_0 + soupsieve conda-forge/noarch::soupsieve-2.5-pyhd8ed1ab_1 + sphinx conda-forge/noarch::sphinx-7.2.6-pyhd8ed1ab_0 + sphinxcontrib-app~ conda-forge/noarch::sphinxcontrib-applehelp-1.0.7-pyhd8ed1ab_0 + sphinxcontrib-dev~ conda-forge/noarch::sphinxcontrib-devhelp-1.0.5-pyhd8ed1ab_0 + sphinxcontrib-htm~ conda-forge/noarch::sphinxcontrib-htmlhelp-2.0.4-pyhd8ed1ab_0 + sphinxcontrib-jsm~ conda-forge/noarch::sphinxcontrib-jsmath-1.0.1-pyhd8ed1ab_0 + sphinxcontrib-qth~ conda-forge/noarch::sphinxcontrib-qthelp-1.0.6-pyhd8ed1ab_0 + sphinxcontrib-ser~ conda-forge/noarch::sphinxcontrib-serializinghtml-1.1.9-pyhd8ed1ab_0 + spyder conda-forge/win-64::spyder-5.5.0-py311h1ea47a8_3 + spyder-kernels conda-forge/noarch::spyder-kernels-2.5.0-win_pyh7428d3b_0 + stack_data conda-forge/noarch::stack_data-0.6.2-pyhd8ed1ab_0 + text-unidecode conda-forge/noarch::text-unidecode-1.3-pyhd8ed1ab_1 + textdistance conda-forge/noarch::textdistance-4.5.0-pyhd8ed1ab_0 + three-merge conda-forge/noarch::three-merge-0.1.1-pyh9f0ad1d_0 + tinycss2 conda-forge/noarch::tinycss2-1.2.1-pyhd8ed1ab_0 + toml conda-forge/noarch::toml-0.10.2-pyhd8ed1ab_0 + tomli conda-forge/noarch::tomli-2.0.1-pyhd8ed1ab_0 + tomlkit conda-forge/noarch::tomlkit-0.12.3-pyha770c72_0 + tornado conda-forge/win-64::tornado-6.3.3-py311ha68e1ae_1 + traitlets conda-forge/noarch::traitlets-5.13.0-pyhd8ed1ab_0 + types-python-date~ conda-forge/noarch::types-python-dateutil-2.8.19.14-pyhd8ed1ab_0 + typing-extensions conda-forge/noarch::typing-extensions-4.8.0-hd8ed1ab_0 + typing_extensions conda-forge/noarch::typing_extensions-4.8.0-pyha770c72_0 + ujson conda-forge/win-64::ujson-5.8.0-py311h12c1d0e_0 + urllib3 conda-forge/noarch::urllib3-2.1.0-pyhd8ed1ab_0 + watchdog conda-forge/win-64::watchdog-3.0.0-py311h1ea47a8_1 + wcwidth conda-forge/noarch::wcwidth-0.2.12-pyhd8ed1ab_0 + webencodings conda-forge/noarch::webencodings-0.5.1-pyhd8ed1ab_2 + whatthepatch conda-forge/noarch::whatthepatch-1.0.5-pyhd8ed1ab_0 + win_inet_pton conda-forge/noarch::win_inet_pton-1.1.0-pyhd8ed1ab_6 + xorg-libxau conda-forge/win-64::xorg-libxau-1.0.11-hcd874cb_0 + xorg-libxdmcp conda-forge/win-64::xorg-libxdmcp-1.1.3-hcd874cb_0 + yaml conda-forge/win-64::yaml-0.2.5-h8ffe710_2 + yapf conda-forge/noarch::yapf-0.40.1-pyhd8ed1ab_0 + zeromq conda-forge/win-64::zeromq-4.3.5-h63175ca_0 + zipp conda-forge/noarch::zipp-3.17.0-pyhd8ed1ab_0 + zstd conda-forge/win-64::zstd-1.5.5-h12be248_0 + + + Proceed ([y]/n)? y + + + Downloading and Extracting Packages + + Preparing transaction: done + Verifying transaction: done + Executing transaction: done + + +.. note:: + + * When you open the Spyder IDE, make sure you select the right one. If you + have multiple installations of Spyder, the specific environment of Spyder + will be mentioned between parentheses. + * You can only open one Spyder window at the same time. Please make sure you + have closed all other Spyder windows (specifically those corresponding + to a different environment) before opening :code:`Spyder (pydft-env)`. + * After installing :program:`PyDFT` into its own environment, please revisit + the "Simple calculation" section on this page and check that your + installation is working. \ No newline at end of file diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..922152e --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/scripts/00-start.py b/docs/scripts/00-start.py new file mode 100644 index 0000000..2042a98 --- /dev/null +++ b/docs/scripts/00-start.py @@ -0,0 +1,10 @@ +import numpy as np +from pydft import MoleculeBuilder, DFT + +# perform DFT calculation on the CO molecule +co = MoleculeBuilder().from_name("CO") +dft = DFT(co, basis='sto3g', verbose=False) +en = dft.scf(1e-6) + +# print total energy +print("Total electronic energy: %12.6f Ht" % en) \ No newline at end of file diff --git a/docs/scripts/00-verbose.py b/docs/scripts/00-verbose.py new file mode 100644 index 0000000..26d627a --- /dev/null +++ b/docs/scripts/00-verbose.py @@ -0,0 +1,7 @@ +import numpy as np +from pydft import MoleculeBuilder, DFT + +# perform DFT calculation on the CO molecule +co = MoleculeBuilder().from_name("CO") +dft = DFT(co, basis='sto3g', verbose=True) +en = dft.scf(1e-5) \ No newline at end of file diff --git a/docs/scripts/00-xc.py b/docs/scripts/00-xc.py new file mode 100644 index 0000000..4c1dd8b --- /dev/null +++ b/docs/scripts/00-xc.py @@ -0,0 +1,13 @@ +import numpy as np +from pydft import MoleculeBuilder, DFT + +# perform DFT calculation on the CO molecule +co = MoleculeBuilder().from_name("CO") + +# use SVWM XC functional +dft = DFT(co, basis='sto3g', verbose=False, functional='svwn5') +print('SVWN: ', dft.scf(1e-5), 'Ht') + +# use PBE XC functional +dft = DFT(co, basis='sto3g', verbose=False, functional='pbe') +print('PBE: ', dft.scf(1e-5), 'Ht') \ No newline at end of file diff --git a/docs/scripts/01-electron-density.py b/docs/scripts/01-electron-density.py new file mode 100644 index 0000000..3961267 --- /dev/null +++ b/docs/scripts/01-electron-density.py @@ -0,0 +1,38 @@ +import numpy as np +from pydft import MoleculeBuilder, DFT +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable + +# perform DFT calculation on the CO molecule +co = MoleculeBuilder().from_name("CO") +dft = DFT(co, basis='sto3g', verbose=False) +en = dft.scf(1e-6) + +# generate grid of points and calculate the electron density for these points +sz = 4 # size of the domain +npts = 100 # number of sampling points per cartesian direction + +# produce meshgrid for the xz-plane +x = np.linspace(-sz/2,sz/2,npts) +zz, xx = np.meshgrid(x, x, indexing='ij') +gridpoints = np.zeros((zz.shape[0], xx.shape[1], 3)) +gridpoints[:,:,0] = xx +gridpoints[:,:,2] = zz +gridpoints[:,:,1] = np.zeros_like(xx) # set y-values to 0 +gridpoints = gridpoints.reshape((-1,3)) + +# calculate (logarithmic) scalar field and convert if back to an 2D array +density = dft.get_density_at_points(gridpoints) +density = np.log10(density.reshape((npts, npts))) + +# build contour plot +fig, ax = plt.subplots(1,1, dpi=144, figsize=(4,4)) +im = ax.contourf(x, x, density, levels=np.linspace(-3,3,13, endpoint=True)) +ax.contour(x, x, density, colors='black', levels=np.linspace(-3,3,13, endpoint=True)) +ax.set_aspect('equal', 'box') +ax.set_xlabel('x-coordinate [a.u.]') +ax.set_ylabel('z-coordinate [a.u.]') +divider = make_axes_locatable(ax) +cax = divider.append_axes('right', size='5%', pad=0.05) +fig.colorbar(im, cax=cax, orientation='vertical') +ax.set_title('Electron density') \ No newline at end of file diff --git a/docs/scripts/02-electron-density-gradient.py b/docs/scripts/02-electron-density-gradient.py new file mode 100644 index 0000000..48abf06 --- /dev/null +++ b/docs/scripts/02-electron-density-gradient.py @@ -0,0 +1,38 @@ +import numpy as np +from pydft import MoleculeBuilder, DFT +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable + +# perform DFT calculation on the CO molecule +co = MoleculeBuilder().from_name("CO") +dft = DFT(co, basis='sto3g', verbose=False) +en = dft.scf(1e-6) + +# generate grid of points and calculate the electron density for these points +sz = 4 # size of the domain +npts = 100 # number of sampling points per cartesian direction + +# produce meshgrid for the xz-plane +x = np.linspace(-sz/2,sz/2,npts) +zz, xx = np.meshgrid(x, x, indexing='ij') +gridpoints = np.zeros((zz.shape[0], xx.shape[1], 3)) +gridpoints[:,:,0] = xx +gridpoints[:,:,2] = zz +gridpoints[:,:,1] = np.zeros_like(xx) # set y-values to 0 +gridpoints = gridpoints.reshape((-1,3)) + +# calculate (logarithmic) scalar field and convert if back to an 2D array +gradient = np.linalg.norm(dft.get_gradient_at_points(gridpoints), axis=1) +gradient = np.log10(gradient.reshape((npts, npts))) + +# build contour plot +fig, ax = plt.subplots(1,1, dpi=144, figsize=(4,4)) +im = ax.contourf(x, x, gradient, levels=np.linspace(-2,4,7, endpoint=True)) +ax.contour(x, x, gradient, colors='black', levels=np.linspace(-2,4,7, endpoint=True)) +ax.set_aspect('equal', 'box') +ax.set_xlabel('x-coordinate [a.u.]') +ax.set_ylabel('z-coordinate [a.u.]') +divider = make_axes_locatable(ax) +cax = divider.append_axes('right', size='5%', pad=0.05) +fig.colorbar(im, cax=cax, orientation='vertical') +ax.set_title('Electron density gradient magnitude') \ No newline at end of file diff --git a/docs/scripts/03-energy-decomposition.py b/docs/scripts/03-energy-decomposition.py new file mode 100644 index 0000000..ea0f28f --- /dev/null +++ b/docs/scripts/03-energy-decomposition.py @@ -0,0 +1,34 @@ +import numpy as np +from pydft import MoleculeBuilder, DFT + +co = MoleculeBuilder().from_name("CO") +dft = DFT(co, basis='sto3g', verbose=False) +en = dft.scf(1e-4) +print("Total electronic energy: %12.6f Ht" % en) +print() + +# retrieve molecular matrices +res = dft.get_data() +P = res['P'] +T = res['T'] +V = res['V'] +J = res['J'] + +# calculate energy terms +Et = np.einsum('ji,ij', P, T) +Ev = np.einsum('ji,ij', P, V) +Ej = 0.5 * np.einsum('ji,ij', P, J) +Ex = res['Ex'] +Ec = res['Ec'] +Exc = res['Exc'] +Enuc = res['enucrep'] + +print('Kinetic energy: %12.6f Ht' % Et) +print('Nuclear attraction: %12.6f Ht' % Ev) +print('Electron-electron repulsion: %12.6f Ht' % Ej) +print('Exchange energy: %12.6f Ht' % (Ex)) +print('Correlation energy: %12.6f Ht' % (Ec)) +print('Exchange-correlation energy: %12.6f Ht' % (Exc)) +print('Nucleus-nucleus repulsion: %12.6f Ht' % (Enuc)) +print() +print('Sum: %12.6f Ht' % (Et + Ev + Ej + Exc + Enuc)) diff --git a/docs/scripts/03-matrices.py b/docs/scripts/03-matrices.py new file mode 100644 index 0000000..362d5b1 --- /dev/null +++ b/docs/scripts/03-matrices.py @@ -0,0 +1,50 @@ +import numpy as np +from pydft import MoleculeBuilder, DFT +import matplotlib.pyplot as plt + +def main(): + # perform DFT calculation on the CO molecule + co = MoleculeBuilder().from_name("CO") + dft = DFT(co, basis='sto3g', verbose=False) + en = dft.scf(1e-6) + + # build list of basis functions + labels = [] + for a in co.atoms: + for o in ['1s', '2s', '2px', '2py', '2pz']: + labels.append('%s - %s' % (a[0],o)) + + fig, ax = plt.subplots(1, 2, dpi=144, figsize=(8,4)) + plot_matrix(ax[0], dft.get_data()['S'], xlabels=labels, ylabels=labels, title='Overlap matrix') + plot_matrix(ax[1], dft.get_data()['F'], xlabels=labels, ylabels=labels, title='Hamiltonian matrix') + +def plot_matrix(ax, mat, xlabels=None, ylabels=None, title = None, xlabelrot=90): + """ + Produce plot of matrix + """ + ax.imshow(mat, vmin=-1, vmax=1, cmap='PiYG') + for i in range(mat.shape[0]): + for j in range(mat.shape[1]): + ax.text(i, j, '%.2f' % mat[j,i], ha='center', va='center', + fontsize=7) + ax.set_xticks([]) + ax.set_yticks([]) + ax.hlines(np.arange(1, mat.shape[0])-0.5, -0.5, mat.shape[0] - 0.5, + color='black', linestyle='--', linewidth=1) + ax.vlines(np.arange(1, mat.shape[0])-0.5, -0.5, mat.shape[0] - 0.5, + color='black', linestyle='--', linewidth=1) + + ax.set_xticks(np.arange(0, mat.shape[0])) + if xlabels is not None: + ax.set_xticklabels(xlabels, rotation=xlabelrot) + ax.set_yticks(np.arange(0, mat.shape[0])) + + if ylabels is not None: + ax.set_yticklabels(ylabels, rotation=0) + ax.tick_params(axis='both', which='major', labelsize=7) + + if title is not None: + ax.set_title(title) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/docs/scripts/04-kohn-sham-orbitals.py b/docs/scripts/04-kohn-sham-orbitals.py new file mode 100644 index 0000000..daeeac1 --- /dev/null +++ b/docs/scripts/04-kohn-sham-orbitals.py @@ -0,0 +1,50 @@ +import numpy as np +from pydft import MoleculeBuilder, DFT +import matplotlib.pyplot as plt + +def main(): + # perform DFT calculation on the CO molecule + co = MoleculeBuilder().from_name("CO") + dft = DFT(co, basis='sto3g', verbose=False) + en = dft.scf(1e-6) + + # build list of basis functions + labels = [] + for a in co.atoms: + for o in ['1s', '2s', '2px', '2py', '2pz']: + labels.append('%s - %s' % (a[0],o)) + + fig, ax = plt.subplots(1, 1, dpi=144, figsize=(4,4)) + plot_matrix(ax, dft.get_data()['C'], xlabels=[str(i) for i in np.arange(1,11)], + ylabels=labels, title='Coefficient matrix') + +def plot_matrix(ax, mat, xlabels=None, ylabels=None, title = None, xlabelrot=90): + """ + Produce plot of matrix + """ + ax.imshow(mat, vmin=-1, vmax=1, cmap='PiYG') + for i in range(mat.shape[0]): + for j in range(mat.shape[1]): + ax.text(i, j, '%.2f' % mat[j,i], ha='center', va='center', + fontsize=7) + ax.set_xticks([]) + ax.set_yticks([]) + ax.hlines(np.arange(1, mat.shape[0])-0.5, -0.5, mat.shape[0] - 0.5, + color='black', linestyle='--', linewidth=1) + ax.vlines(np.arange(1, mat.shape[0])-0.5, -0.5, mat.shape[0] - 0.5, + color='black', linestyle='--', linewidth=1) + + ax.set_xticks(np.arange(0, mat.shape[0])) + if xlabels is not None: + ax.set_xticklabels(xlabels, rotation=xlabelrot) + ax.set_yticks(np.arange(0, mat.shape[0])) + + if ylabels is not None: + ax.set_yticklabels(ylabels, rotation=0) + ax.tick_params(axis='both', which='major', labelsize=7) + + if title is not None: + ax.set_title(title) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/docs/scripts/05-orbital-contour-plot.py b/docs/scripts/05-orbital-contour-plot.py new file mode 100644 index 0000000..1438498 --- /dev/null +++ b/docs/scripts/05-orbital-contour-plot.py @@ -0,0 +1,54 @@ +import numpy as np +from pydft import MoleculeBuilder, DFT +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable + +def main(): + # perform DFT calculation on the CO molecule + co = MoleculeBuilder().from_name("CO") + dft = DFT(co, basis='sto3g', verbose=False) + en = dft.scf(1e-6) + + # grab molecular orbital energies and coefficients + orbc = dft.get_data()['C'] + orbe = dft.get_data()['orbe'] + + # generate grid of points and calculate the electron density for these points + sz = 4 # size of the domain + npts = 100 # number of sampling points per cartesian direction + + # produce meshgrid for the xz-plane + x = np.linspace(-sz/2,sz/2,npts) + zz, xx = np.meshgrid(x, x, indexing='ij') + gridpoints = np.zeros((zz.shape[0], xx.shape[1], 3)) + gridpoints[:,:,0] = xx + gridpoints[:,:,2] = zz + gridpoints[:,:,1] = np.zeros_like(xx) # set y-values to 0 + gridpoints = gridpoints.reshape((-1,3)) + + # grab a copy of the MolecularGrid object + molgrid = dft.get_molgrid_copy() + + fig, ax = plt.subplots(2, 5, dpi=144, figsize=(16,6)) + for i in range(len(orbe)): + m_ax = ax[i//5, i%5] + + # calculate field + field = molgrid.get_amplitude_at_points(gridpoints, orbc[:,i]).reshape((npts, npts)) + + # plot field + levels = np.linspace(-2,2,17, endpoint=True) + im = m_ax.contourf(x, x, field, levels=levels, cmap='PiYG') + m_ax.contour(x, x, field, colors='black') + m_ax.set_aspect('equal', 'box') + m_ax.set_xlabel('x-coordinate [a.u.]') + m_ax.set_ylabel('z-coordinate [a.u.]') + divider = make_axes_locatable(m_ax) + cax = divider.append_axes('right', size='5%', pad=0.05) + fig.colorbar(im, cax=cax, orientation='vertical') + m_ax.set_title('MO %i: %6.4f Ht' % (i+1, orbe[i])) + + plt.tight_layout() + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/docs/scripts/06-becke-grid.py b/docs/scripts/06-becke-grid.py new file mode 100644 index 0000000..8840bd8 --- /dev/null +++ b/docs/scripts/06-becke-grid.py @@ -0,0 +1,40 @@ +from pydft import MolecularGrid +from pyqint import MoleculeBuilder +import numpy as np +import matplotlib.pyplot as plt + +# construct molecule +mol = MoleculeBuilder().from_name('benzene') +cgfs, atoms = mol.build_basis('sto3g') + +# construct molecular grid +molgrid = MolecularGrid(atoms, cgfs) + +# produce grid of sampling points to calculate the atomic +# weight coefficients for +N = 150 +sz = 8 +x = np.linspace(-sz,sz,N) +xv,yv = np.meshgrid(x,x) +points = np.array([[x,y,0] for x,y in zip(xv.flatten(),yv.flatten())]) + +# calculate the atomic weights +mweights = molgrid.calculate_weights_at_points(points, k=3) + +# plot the atomic weights +plt.figure(dpi=144) +plt.imshow(np.max(mweights,axis=0).reshape((N,N)), + extent=(-sz,sz,-sz,sz), interpolation='bicubic') +plt.xlabel('x [a.u.]') +plt.ylabel('y [a.u.]') +plt.colorbar() +plt.grid(linestyle='--', color='black', alpha=0.5) + +# add the atoms to the plot +r = np.zeros((len(atoms), 3)) +for i,at in enumerate(atoms): + r[i] = at[0] +plt.scatter(r[0:6,0], r[0:6,1], s=50.0, color='grey', edgecolor='black') +plt.scatter(r[6:12,0], r[6:12,1], s=50.0, color='white', edgecolor='black') + +plt.tight_layout() \ No newline at end of file diff --git a/docs/scripts/06-grid-points.py b/docs/scripts/06-grid-points.py new file mode 100644 index 0000000..bf17d94 --- /dev/null +++ b/docs/scripts/06-grid-points.py @@ -0,0 +1,31 @@ +from pydft import MolecularGrid +from pyqint import MoleculeBuilder +import matplotlib.pyplot as plt + +# construct molecule +mol = MoleculeBuilder().from_name('co') +cgfs, atoms = mol.build_basis('sto3g') + +# construct molecular grid +molgrid = MolecularGrid(atoms, cgfs, nshells=16, nangpts=74) +molgrid.initialize() + +# obtain the grid points +gridpoints = molgrid.get_grid_coordinates() + +# plot the atomic weights +colors = '#DD0000', '#222222' +fig = plt.figure(dpi=300, figsize=(8,6)) +ax = fig.add_subplot(projection='3d') +for i in range(0, len(gridpoints)): + ax.scatter(gridpoints[i][:,0], gridpoints[i][:,1], gridpoints[i][:,2], + s = 1.5, alpha=0.5, color=colors[i]) + +# set axes +ax.set_xlim(-5,5) +ax.set_ylim(-5,5) +ax.set_zlim(-5,5) +ax.set_xlabel('x [a.u.]') +ax.set_ylabel('y [a.u.]') +ax.set_zlabel('z [a.u.]') +ax.set_box_aspect(aspect=None, zoom=0.8) \ No newline at end of file diff --git a/docs/scripts/07-hartee-potential.py b/docs/scripts/07-hartee-potential.py new file mode 100644 index 0000000..1062828 --- /dev/null +++ b/docs/scripts/07-hartee-potential.py @@ -0,0 +1,21 @@ +import numpy as np +from pydft import MoleculeBuilder, DFT +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable + +def main(): + # perform DFT calculation on the CO molecule + co = MoleculeBuilder().from_name("CO") + dft = DFT(co, basis='sto3g', verbose=False) + en = dft.scf(1e-6) + + # get the projection onto the spherical harmonics + lmprojection = dft.get_molgrid_copy().get_rho_lm_atoms() + + fig, ax = plt.subplots(1, 2, dpi=144, figsize=(6,4)) + plot_spherical_harmonic_projection(ax[1], fig, lmprojection[0], title='O') + plot_spherical_harmonic_projection(ax[0], fig, lmprojection[1], title='C') + plt.tight_layout() + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/docs/scripts/07-spherical-harmonics-projection.py b/docs/scripts/07-spherical-harmonics-projection.py new file mode 100644 index 0000000..daacd43 --- /dev/null +++ b/docs/scripts/07-spherical-harmonics-projection.py @@ -0,0 +1,38 @@ +import numpy as np +from pydft import MoleculeBuilder, DFT +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable + +def main(): + # perform DFT calculation on the CO molecule + co = MoleculeBuilder().from_name("CO") + dft = DFT(co, basis='sto3g', verbose=False) + en = dft.scf(1e-6) + + # get the projection onto the spherical harmonics + lmprojection = dft.get_molgrid_copy().get_rho_lm_atoms() + + fig, ax = plt.subplots(1, 2, dpi=144, figsize=(6,4)) + plot_spherical_harmonic_projection(ax[1], fig, lmprojection[0], title='O') + plot_spherical_harmonic_projection(ax[0], fig, lmprojection[1], title='C') + plt.tight_layout() + + +def plot_spherical_harmonic_projection(ax, fig, mat, minval=5e-2, title=None): + im = ax.imshow(np.flip(mat, axis=0), origin='lower', cmap='PiYG', vmin=-minval,vmax=minval) + ax.grid(linestyle='--', color='black', alpha=0.8) + ax.set_xticks([-0.5,0.5,3.5,8.5,15.5], labels=['s','p','d','f','g']) + ax.set_yticks(np.arange(0, mat.shape[0], 5)-0.5, + labels=np.arange(0, mat.shape[0], 5)) + ax.set_ylabel('Radial points $r_{i}$') + ax.set_xlim(-.5,15) + ax.set_xlabel('Spherical harmonics (lm)') + divider = make_axes_locatable(ax) + cax = divider.append_axes('right', size='5%', pad=0.05) + fig.colorbar(im, cax=cax, orientation='vertical') + + if title is not None: + ax.set_title(title) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/docs/scripts/08-correlation-potential.py b/docs/scripts/08-correlation-potential.py new file mode 100644 index 0000000..4c8a8a8 --- /dev/null +++ b/docs/scripts/08-correlation-potential.py @@ -0,0 +1,47 @@ +from pydft import DFT +from pyqint import MoleculeBuilder +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable + +# perform DFT calculation on the CO molecule +co = MoleculeBuilder().from_name("CO") +dft = DFT(co, basis='sto3g', verbose=False) +en = dft.scf(1e-6) + +# grab molecular orbital energies and coefficients +orbc = dft.get_data()['C'] +orbe = dft.get_data()['orbe'] + +# generate grid of points and calculate the electron density for these points +sz = 4 # size of the domain +npts = 100 # number of sampling points per cartesian direction + +# produce meshgrid for the xz-plane +x = np.linspace(-sz/2,sz/2,npts) +zz, xx = np.meshgrid(x, x, indexing='ij') +gridpoints = np.zeros((zz.shape[0], xx.shape[1], 3)) +gridpoints[:,:,0] = xx +gridpoints[:,:,2] = zz +gridpoints[:,:,1] = np.zeros_like(xx) # set y-values to 0 +gridpoints = gridpoints.reshape((-1,3)) + +# grab a copy of the MolecularGrid object +molgrid = dft.get_molgrid_copy() + +# construct exchange potential +field = molgrid.get_correlation_potential_at_points(gridpoints, dft.get_data()['P']) +field = field.reshape((npts, npts)) + +# plot field +fig, ax = plt.subplots(1, 1, dpi=144, figsize=(4,4)) +levels = np.linspace(-0.2,0.2,101, endpoint=True) +im = ax.contourf(x, x, field, levels=levels, cmap='PiYG') +ax.contour(x, x, field, levels=levels, colors='black', linewidths=0.5) +ax.set_aspect('equal', 'box') +ax.set_xlabel('x-coordinate [a.u.]') +ax.set_ylabel('z-coordinate [a.u.]') +divider = make_axes_locatable(ax) +cax = divider.append_axes('right', size='5%', pad=0.05) +fig.colorbar(im, cax=cax, orientation='vertical') +ax.set_title('Correlation potential') \ No newline at end of file diff --git a/docs/scripts/08-exchange-potential.py b/docs/scripts/08-exchange-potential.py new file mode 100644 index 0000000..6383284 --- /dev/null +++ b/docs/scripts/08-exchange-potential.py @@ -0,0 +1,47 @@ +from pydft import DFT +from pyqint import MoleculeBuilder +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable + +# perform DFT calculation on the CO molecule +co = MoleculeBuilder().from_name("CO") +dft = DFT(co, basis='sto3g', verbose=False) +en = dft.scf(1e-6) + +# grab molecular orbital energies and coefficients +orbc = dft.get_data()['C'] +orbe = dft.get_data()['orbe'] + +# generate grid of points and calculate the electron density for these points +sz = 4 # size of the domain +npts = 100 # number of sampling points per cartesian direction + +# produce meshgrid for the xz-plane +x = np.linspace(-sz/2,sz/2,npts) +zz, xx = np.meshgrid(x, x, indexing='ij') +gridpoints = np.zeros((zz.shape[0], xx.shape[1], 3)) +gridpoints[:,:,0] = xx +gridpoints[:,:,2] = zz +gridpoints[:,:,1] = np.zeros_like(xx) # set y-values to 0 +gridpoints = gridpoints.reshape((-1,3)) + +# grab a copy of the MolecularGrid object +molgrid = dft.get_molgrid_copy() + +# construct exchange potential +field = molgrid.get_exchange_potential_at_points(gridpoints, dft.get_data()['P']) +field = field.reshape((npts, npts)) + +# plot field +fig, ax = plt.subplots(1, 1, dpi=144, figsize=(4,4)) +levels = np.linspace(-6,6,25, endpoint=True) +im = ax.contourf(x, x, field, levels=levels, cmap='PiYG') +ax.contour(x, x, field, levels=levels, colors='black', linewidths=0.5) +ax.set_aspect('equal', 'box') +ax.set_xlabel('x-coordinate [a.u.]') +ax.set_ylabel('z-coordinate [a.u.]') +divider = make_axes_locatable(ax) +cax = divider.append_axes('right', size='5%', pad=0.05) +fig.colorbar(im, cax=cax, orientation='vertical') +ax.set_title('Exchange potential') \ No newline at end of file diff --git a/docs/user_interface.rst b/docs/user_interface.rst new file mode 100644 index 0000000..6a3907a --- /dev/null +++ b/docs/user_interface.rst @@ -0,0 +1,349 @@ +.. _user_interface: +.. index:: userinterface + +User Interface +============== + +:program:`PyDFT` is an educational pure-Python module to perform simple density +functional theory (DFT) calculations using a Gaussian basis set. Uniquely, +:program:`PyDFT` exposes many routines that are normally hidden from the user, +which are potentially interesting for the user to understand how a DFT calculations +is performed. We will provide a demonstration here. + +.. note:: + + Every code snippet listed below can be run stand-alone. This does imply + that every code snippet re-performed the initial self-consistent field + calculation. If that is undesirable, the reader is invited to copy + the relevant parts of the code (highlighted) and place these in a single + file. + +Building molecules +------------------ + +Molecules can be built either directly using the :code:`Molecule` class or via +a :code:`MoleculeBuilder` convenience routine. + +.. important:: + + The :code:`Molecule` and :code:`MoleculeBuilder` classes are not implemented + in :program:`PyDFT` but are obtained from the :program:`PyQInt` module. + +Molecule class +############## + +To manually build a molecule, one first need to construct a :code:`Molecule` +object after which one or more atoms can be assigned to the molecule. If no +:code:`unit` is specified, it is assumed that all coordinates are given in +atomic units (i.e. Bohr units). + +.. code-block:: python + :linenos: + + from pyqint import Molecule + + mol = Molecule('co') + mol.add_atom('C', 0.0, 0.0, 0.0, unit='angstrom') + mol.add_atom('O', 0.0, 0.0, 1.2, unit='angstrom') + +MoleculeBuilder class +##################### + +Alternatively, a Molecule can be constructed from the :code:`MoleculeBuilder` +class which uses a name. For more information about the :code:`MoleculeBuilder` +class, please consult the +`PyQInt documentation on the MoleculeBuilder class `_. + +.. code-block:: python + :linenos: + + from pyqint import Molecule + + co = MoleculeBuilder().from_name("CO") + +Calculating electronic energy +----------------------------- + +To start, we perform a high-level calculation of the electronic structure +of the carbon-monoxide molecule using the PBE exchange-correlation functional. +To perform this calculation, we first have to construct a +:class:`pydft.DFT` object which requires a molecule as its input. Next, we use +the :meth:`pydft.DFT.scf` routine to start the self-consistent field calculation. + +.. literalinclude:: scripts/00-start.py + :language: python + :linenos: + +Performing this calculation shows that the total electronic energy for this +system corresponds to:: + + Total electronic energy: -111.147096 Ht + +Showing the electronic steps +############################ + +.. literalinclude:: scripts/00-verbose.py + :language: python + :linenos: + +Executing the script above yields the following output:: + + 001 | Energy: -179.237380 | 0.0020 ms + 002 | Energy: -117.709405 | 0.1130 ms + 003 | Energy: -118.058343 | 0.0900 ms + 004 | Energy: -107.368999 | 0.0860 ms + 005 | Energy: -117.406523 | 0.0868 ms + 006 | Energy: -112.897655 | 0.0880 ms + 007 | Energy: -115.401985 | 0.0870 ms + 008 | Energy: -107.642758 | 0.0860 ms + 009 | Energy: -112.323123 | 0.0870 ms + 010 | Energy: -111.477884 | 0.0860 ms + 011 | Energy: -111.214223 | 0.0860 ms + 012 | Energy: -111.139379 | 0.0850 ms + 013 | Energy: -111.147186 | 0.0860 ms + 014 | Energy: -111.147098 | 0.0860 ms + 015 | Energy: -111.147096 | 0.0850 ms + 016 | Energy: -111.147096 | 0.0895 ms + 017 | Energy: -111.147096 | 0.0855 ms + 018 | Energy: -111.147096 | 0.0860 ms + Stopping SCF cycle, convergence reached. + +Different exchange-correlation functionals +########################################## + +To use a different exchange-correlation functional, we can use the + +.. literalinclude:: scripts/00-xc.py + :language: python + :linenos: + +which yields the following total electronic energies for the :code:`SVWN5` and +:code:`PBE` exchange-correlation functions:: + + SVWN: -111.14709591483225 Ht + PBE: -111.65660426438342 Ht + +Electron density and its gradient +--------------------------------- + +The :class:`pydft.DFT` class exposes its internal matrices and a number of +useful functions which we can readily use to interpret its operation. Let us +start by generating a plot of the electron density using the method +:meth:`pydft.DFT.get_density_at_points`. + +.. literalinclude:: scripts/01-electron-density.py + :language: python + :linenos: + :emphasize-lines: 25 + +Running the above script yields the electron density. + +.. figure:: _static/img/user_interface/01-electron-density.png + +Using largely the same code, we can also readily build the electron density +gradient magnitude using :meth:`pydft.DFT.get_gradient_at_points`. + +.. literalinclude:: scripts/02-electron-density-gradient.py + :language: python + :linenos: + :emphasize-lines: 25 + +.. figure:: _static/img/user_interface/02-electron-density-gradient.png + +Self-consistent field matrices +------------------------------ + +To obtain any of the matrices used in the self-consistent field procedure, +we can invoke the :meth:`pydft.DFT.get_data` method. For example, to visualize +the overlap matrix :math:`\mathbf{S}` and the Fock matrix +:math:`\mathbf{F}`, we can use the script as found below. + +.. literalinclude:: scripts/03-matrices.py + :language: python + :linenos: + :emphasize-lines: 18,19 + +.. figure:: _static/img/user_interface/03-matrices.png + +Energy decomposition +#################### + +The molecular matrices can be used to perform a so-called energy decomposition, +i.e., decompose the total electronic energy into the kinetic, nuclear attraction, +electron-electron repulsion and exchange-correlation energy. + +.. literalinclude:: scripts/03-energy-decomposition.py + :language: python + :linenos: + +The above script yields the following output:: + + Total electronic energy: -111.147096 Ht + + Kinetic energy: 110.216045 Ht + Nuclear attraction: -304.930390 Ht + Electron-electron repulsion: 75.597401 Ht + Exchange energy: -12.055579 Ht + Correlation energy: -1.232665 Ht + Exchange-correlation energy: -13.288244 Ht + Nucleus-nucleus repulsion: 21.258092 Ht + + Sum: -111.147096 Ht + +Kohn-Sham orbitals and their visualization +------------------------------------------ + +Within the Kohn-Sham approximation, the set of molecular orbitals that minimizes +the total electronic energy are found via the diagonalization of the Fock +matrix. + +Coefficient matrix +################## + +These Kohn-Sham orbitals are encoded in the coefficient matrix which +is extracted via the script below. + +.. literalinclude:: scripts/04-kohn-sham-orbitals.py + :language: python + :linenos: + +.. figure:: _static/img/user_interface/04-kohn-sham-orbitals.png + +Contour plots +############# + +These orbitals can be readily visualized, as exemplified in the code below. Here, +we make use of the :meth:`pydft.MolecularGrid.get_amplitude_at_points` method. +To use this method, we first have to retrieve the :class:`pydft.MolecularGrid` +object from the :class:`pydft.DFT` class via its :meth:`pydft.get_molgrid_copy` +method. + +.. literalinclude:: scripts/05-orbital-contour-plot.py + :language: python + :linenos: + :emphasize-lines: 37 + +.. figure:: _static/img/user_interface/05-orbital-contour-plot.png + +.. seealso:: + + To get a better impression of the three-dimensional shape of the molecular + orbitals, it is recommended to produce isosurfaces rather than contour + plots. + +Isosurfaces +########### + +Build isosurfaces + +Analyzing the Becke grid +------------------------ + +All numerical integrations are performed by means of Gauss-Chebychev and Lebedev +quadrature using the Becke grid. + +Atomic fuzzy cells +################## + +It is possible to produce a contour plot of +the fuzzy cells (molecular weights) on a plane using the +:meth:`pydft.MolecularGrid.calculate_weights_at_points` method. An example +is provided below. + +.. literalinclude:: scripts/06-becke-grid.py + :language: python + :linenos: + :emphasize-lines: 22 + +In the script above, we color every point by the maximum value for each of the +atomic weights. When this maximum value is one, this implies that the grid point +belongs to a single atom. When multiple atoms 'share' a grid point, the maximum +value among the atomic weights will be lower than one. + +.. image:: _static/img/user_interface/06-becke-grid.png + +.. note:: + + Producing such a contour plot is only meaningful for planar molecules such + as benzene. For more complex molecules such as methane, it is rather + difficult to make sense of the fuzzy cells upon projection on a plane. + +Grid points +########### + +To obtain the set of grid points on which the numerical integration (quadrature), +is performed, we can invoke the :meth:`pydft.MolecularGrid.get_grid_coordinates` +method. + +.. literalinclude:: scripts/06-grid-points.py + :language: python + :linenos: + :emphasize-lines: 14 + +.. image:: _static/img/user_interface/06-grid-points.png + +Hartree potential +----------------- + +The Hartree potential is calculated from the electron density distribution +using the Poisson equation. + +Projection onto spherical harmonics +################################### + +The Poisson equation is solved on the **atomic grid** by first projecting +the electron density onto the spherical harmonics such that the following +equation holds. + +.. math:: + + \rho(r_{k}, \Omega) = \sum_{lm} \rho_{klm} Y_{lm} + +where :math:`\rho(r_{k}, \Omega)` is the electron density at radius +:math:`r_{k}`, :math:`\rho_{klm}` the linear expansion coefficient at position +:math:`r_{k}` for the spherical harmonic :math:`lm` and :math:`Y_{lm}` a +spherical harmonic with quantum numbers :math:`l` and :math:`m`. + +We can readily visualize and interpret this projection using the +:meth:`pydft.MolecularGrid.get_rho_lm_atoms` as shown by the script below. + +.. literalinclude:: scripts/07-spherical-harmonics-projection.py + :language: python + :linenos: + :emphasize-lines: 13 + +.. image:: _static/img/user_interface/07-spherical-harmonics-projection.png + +Exchange-correlation functional +------------------------------- + +The exchange and correlation potentials depend on the electron density and once +the latter has been established, we can readily visualize these potentials. + +Visualizing exchange potential +############################## + +To construct the exchange potential at arbitrary grid points, we can use the +script as shown below which uses the :meth:`pydft.MolecularGrid.get_exchange_potential_at_points` +method. + +.. literalinclude:: scripts/08-exchange-potential.py + :language: python + :linenos: + :emphasize-lines: 13 + +.. image:: _static/img/user_interface/08-exchange-potential.png + +Visualizing correlation potential +################################# + +In a similar fashion as for the exchange potential, we can use the method +:meth:`pydft.MolecularGrid.get_correlation_potential_at_points` to obtain the +correlation potential field. + +.. literalinclude:: scripts/08-correlation-potential.py + :language: python + :linenos: + :emphasize-lines: 13 + +.. image:: _static/img/user_interface/08-correlation-potential.png \ No newline at end of file diff --git a/examples/.gitignore b/examples/.gitignore new file mode 100644 index 0000000..2211df6 --- /dev/null +++ b/examples/.gitignore @@ -0,0 +1 @@ +*.txt diff --git a/examples/benzene_atomic_cells.py b/examples/benzene_atomic_cells.py new file mode 100644 index 0000000..089ca8b --- /dev/null +++ b/examples/benzene_atomic_cells.py @@ -0,0 +1,40 @@ +# -*- coding: utf-8 -*- +import numpy as np +import matplotlib.pyplot as plt +import os, sys + +# add a reference to load the module +ROOT = os.path.dirname(__file__) +sys.path.insert(1, os.path.join(ROOT, '..')) + +from pydft import MoleculeBuilder, MolecularGrid + +def main(): + # construct molecule + mol = MoleculeBuilder().from_name('benzene') + cgfs, atoms = mol.build_basis('sto3g') + + # construct molecular grid + molgrid = MolecularGrid(atoms, cgfs) + + # produce grid of sampling points to calculate the atomic + # weight coefficients for + N = 100 + sz = 8 + x = np.linspace(-sz,sz,N) + xv,yv = np.meshgrid(x,x) + points = np.array([[x,y,0] for x,y in zip(xv.flatten(),yv.flatten())]) + + # calculate the atomic weights + mweights = molgrid.calculate_weights_at_points(points, k=3) + + # plot the atomic weights + plt.imshow(np.max(mweights,axis=0).reshape((N,N)), + extent=(-sz,sz,-sz,sz)) + plt.xlabel('x') + plt.xlabel('y') + plt.colorbar() + plt.tight_layout() + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/examples/bonding_analysis_co.py b/examples/bonding_analysis_co.py new file mode 100644 index 0000000..b797c2a --- /dev/null +++ b/examples/bonding_analysis_co.py @@ -0,0 +1,48 @@ +import matplotlib.pyplot as plt +import numpy as np +from mpl_toolkits.axes_grid1 import make_axes_locatable +import os,sys + +# add a reference to load the module +ROOT = os.path.dirname(__file__) +sys.path.insert(1, os.path.join(ROOT, '..')) + +from pydft import MoleculeBuilder, DFT + +mol_builder = MoleculeBuilder() +mol = mol_builder.from_name('co') + +# construct dft object +dft = DFT(mol, basis='sto3g') +energy = dft.scf() +C = dft.get_data()['C'] + +molgrid = dft.get_molgrid_copy() +she1 = molgrid.get_spherical_harmonic_expansion_of_amplitude(C[:,4], radial_factor=True) +she2 = molgrid.get_spherical_harmonic_expansion_of_amplitude(C[:,5], radial_factor=True) + +fig, ax = plt.subplots(1,4,dpi=144, figsize=(10,4)) + +for i in range(0,4): + she = she1 if i < 2 else she2 + limit = max(abs(np.min(she[i % 2])), abs(np.max(she[i % 2])) ) + im = ax[i].imshow(she[i % 2,:,0:16], origin='lower', vmin=-limit, vmax=limit, cmap='BrBG') + + title = '$\psi_{4}$' if i < 2 else '$\psi_{5}$' + title += '/ C' if i % 2 == 0 else '/ O' + ax[i].set_title(title) + + ax[i].set_ylabel(r'$r_{i}^{2} \rho_{lm} (r_{i})$') + ax[i].set_xlabel('$Y_{lm}$') + + ax[i].set_xticks([-0.5,0.5,3.5,8.5]) + ax[i].set_xticklabels(['s','p','d','f']) + ax[i].set_xticks(np.arange(0,16), minor=True) + ax[i].grid(linestyle='--', color='black', alpha=0.5) + + # add colorbars + divider = make_axes_locatable(ax[i]) + cax = divider.append_axes('right', size='5%', pad=0.05) + fig.colorbar(im, cax=cax, orientation='vertical') + +plt.tight_layout() \ No newline at end of file diff --git a/examples/co.py b/examples/co.py new file mode 100644 index 0000000..f00b1c5 --- /dev/null +++ b/examples/co.py @@ -0,0 +1,18 @@ +# -*- coding: utf-8 -*- +import os,sys + +# add a reference to load the module +ROOT = os.path.dirname(__file__) +sys.path.insert(1, os.path.join(ROOT, '..')) + +from pydft import MoleculeBuilder, DFT + +# +# Example: Calculate total electronic energy for CO using standard +# settings. +# + +CO = MoleculeBuilder().from_name("CO") +dft = DFT(CO, basis='sto3g') +en = dft.scf(1e-4) +print("Total electronic energy: %f Ht" % en) \ No newline at end of file diff --git a/examples/co_screening.py b/examples/co_screening.py new file mode 100644 index 0000000..dd61e68 --- /dev/null +++ b/examples/co_screening.py @@ -0,0 +1,38 @@ +# -*- coding: utf-8 -*- +import os,sys +import numpy as np +import time + +# add a reference to load the module +ROOT = os.path.dirname(__file__) +sys.path.insert(1, os.path.join(ROOT, '..')) + +from pydft import MoleculeBuilder, DFT + +# +# Example: Calculate total electronic energy for CO using standard +# settings. +# + +CO = MoleculeBuilder().from_name("CO") + +rshells = [8,16,32,64,128] +ang = [6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810] + +energies = np.zeros((len(rshells), len(ang))) +timestats = np.zeros((len(rshells), len(ang))) + +for i,nr in enumerate(rshells): + for j,a in enumerate(ang): + start_time = time.time() + dft = DFT(CO, basis='sto3g', nshells=nr, nangpts=a) + en = dft.scf(1e-4) + end_time = time.time() + elapsed_time = end_time - start_time + timestats[i,j] = elapsed_time + energies[i,j] = en + + print(nr, a, en, elapsed_time) + +np.savetxt('energies.txt', energies) +np.savetxt('timestats.txt', timestats) \ No newline at end of file diff --git a/examples/co_xc.py b/examples/co_xc.py new file mode 100644 index 0000000..41f9f9f --- /dev/null +++ b/examples/co_xc.py @@ -0,0 +1,17 @@ +# -*- coding: utf-8 -*- +from pydft import MoleculeBuilder,DFT + +# +# Example: Calculate total electronic energy for CO using different +# XC-functionals. +# + +CO = MoleculeBuilder().from_name("CO") +dft = DFT(CO, basis='sto3g', functional='svwn5') +en = dft.scf(1e-4) +print("Total electronic energy (SVWN5): %f Ht" % en) + +CO = MoleculeBuilder().from_name("CO") +dft = DFT(CO, basis='sto3g', functional='pbe') +en = dft.scf(1e-4) +print("Total electronic energy (PBE): %f Ht" % en) \ No newline at end of file diff --git a/examples/dft_vs_hf.py b/examples/dft_vs_hf.py new file mode 100644 index 0000000..665fc51 --- /dev/null +++ b/examples/dft_vs_hf.py @@ -0,0 +1,16 @@ +# -*- coding: utf-8 -*- +import os,sys + +# add a reference to load the module +ROOT = os.path.dirname(__file__) +sys.path.insert(1, os.path.join(ROOT, '..')) + +from pydft import MoleculeBuilder, DFT + +co = MoleculeBuilder().from_name("CO") +dft = DFT(co, basis='sto3g') +en = dft.scf(1e-4) +print("Total electronic energy: %f Ht" % en) + +res = dft.get_data() +print(res.keys()) \ No newline at end of file diff --git a/examples/h2.py b/examples/h2.py new file mode 100644 index 0000000..b9a6560 --- /dev/null +++ b/examples/h2.py @@ -0,0 +1,18 @@ +# -*- coding: utf-8 -*- +import os,sys + +# add a reference to load the module +ROOT = os.path.dirname(__file__) +sys.path.insert(1, os.path.join(ROOT, '..')) + +from pydft import MoleculeBuilder, DFT + +# +# Example: Calculate total electronic energy for CO using standard +# settings. +# + +CO = MoleculeBuilder().from_name("H2") +dft = DFT(CO, basis='sto3g') +en = dft.scf(1e-4) +print("Total electronic energy: %f Ht" % en) \ No newline at end of file diff --git a/meta.yaml b/meta.yaml new file mode 100644 index 0000000..c38c083 --- /dev/null +++ b/meta.yaml @@ -0,0 +1,42 @@ +package: + name: "pydft" + version: "0.6.1" + +source: + path: . + +build: + noarch: python + +requirements: + build: + - python>=3.9 + + host: + - pip + - python>=3.9 + - setuptools + - numpy + + run: + - python>=3.9 + - scipy + - numpy + +test: + requires: + - numpy + - nose + - scipy + - conda + - packaging + - python>=3.9 + source_files: + - tests/*.py + +about: + home: https://gitlab.tue.nl/ifilot/pydft + license: GPL3 + license_family: GPL + summary: Python package for performing simple DFT calculations + description: See the package README.md for more information. \ No newline at end of file diff --git a/pydft/__init__.py b/pydft/__init__.py new file mode 100644 index 0000000..9f51685 --- /dev/null +++ b/pydft/__init__.py @@ -0,0 +1,7 @@ +from .dft import DFT +from .atomicgrid import AtomicGrid +from .moleculargrid import MolecularGrid +from pyqint import MoleculeBuilder +from pyqint import Molecule + +from ._version import __version__ \ No newline at end of file diff --git a/pydft/_version.py b/pydft/_version.py new file mode 100644 index 0000000..8411e55 --- /dev/null +++ b/pydft/_version.py @@ -0,0 +1 @@ +__version__ = '0.6.1' diff --git a/pydft/angulargrid.py b/pydft/angulargrid.py new file mode 100644 index 0000000..f6301ac --- /dev/null +++ b/pydft/angulargrid.py @@ -0,0 +1,61 @@ +# -*- coding: utf-8 -*- + +import os +import numpy as np +import pylebedev + +class AngularGrid: + def __init__(self): + """ + Construct the class + """ + + # check if PyLebedev is available + self.__leblib = pylebedev.PyLebedev() + + self.__lebedev_orders = [3, 5, 7, 9, 11, 13, 15, 17, 19, + 21, 23, 25, 27, 29, 31, 35, 41, 47, + 53, 59, 65, 71, 77, 83, 89, 95, 101, + 107, 113, 119, 125, 131] + + # lebedev number of points + self.__numpoints = [6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, + 230, 266, 302, 350, 434, 590, 770, 974, 1202, + 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, + 4334, 4802, 5294, 5810] + + self.__coeff = dict() + + def get_coefficients(self, numpoints:int): + """ + Get the coefficients for Lebedev quadrature on the unit sphere + by number of points of that order + """ + if numpoints in self.__numpoints: + # load the data points from file if they do not yet exist + if numpoints not in self.__coeff.keys(): + order = int(self.__lebedev_orders[self.__numpoints.index(numpoints)]) + self.__add_coefficients(order) + + # return datapoints + return self.__coeff[str(numpoints)] + else: + raise Exception('There is no Lebedev order with %s number of points' % numpoints) + + def get_dataset_sizes(self): + """ + Get the number of points on the unit sphere for each of the + Lebedev orders + """ + return [int(key) for key in self.__coeff.keys()] + + def __add_coefficients(self, order): + """ + Add coefficients from file and store these in a dictionary + """ + points, weights = self.__leblib.get_points_and_weights(order) + angles, weights = self.__leblib.get_points_and_weights(order, solid_angles = True) + + data = np.vstack([angles[:,0], angles[:,1], points[:,0], points[:,1], points[:,2], weights]).transpose() + key = str(len(data)) + self.__coeff[key] = data \ No newline at end of file diff --git a/pydft/atomicgrid.py b/pydft/atomicgrid.py new file mode 100644 index 0000000..144d4c0 --- /dev/null +++ b/pydft/atomicgrid.py @@ -0,0 +1,470 @@ +# -*- coding: utf-8 -*- + +import numpy as np +from scipy.interpolate import CubicSpline +from . import bragg_slater +from .angulargrid import AngularGrid +from .spherical_harmonics import spherical_harmonic + +class AtomicGrid: + def __init__(self, at, nshells:int=32, nangpts:int=110, lmax:int=8): + """ + Initialize the class + + at: atom + nshells: number of radial shells + nangpts: number of angular points + """ + self.__atom = at + self.__cp = at[0] # center of the atomic grid (= position of the atom) + self.__aidx = at[1] # element id + self.__nshells = nshells + + # set coefficients and parameters + self.__set_bragg_slater_radius() + self.__build_chebychev_grid() + self.__lebedev_coeff = self.__load_lebedev_coefficients(nangpts) + + # store the positions of the angular points, their weights for the + # Lebedev integration and the value of the spherical harmonics at + # these points + self.__usangles = [np.array([p[0], p[1]]) for p in self.__lebedev_coeff] + self.__angpts = [np.array([p[2], p[3], p[4]]) for p in self.__lebedev_coeff] + self.__wangpts = np.array([p[5] for p in self.__lebedev_coeff]) + self.__lmax = lmax + self.__build_ylm() + self.__build_weight_grid() + self.__mweights = np.ones_like(self.__wgrid) + + def get_charge(self): + """ + Get the atomic charge + """ + return self.__aidx + + def get_radial_grid(self): + """ + Return radial grid + """ + return self.__rr + + def get_full_grid(self): + """ + Get all the grid points as a list of positions (Nx3) matrix + """ + return np.multiply.outer(self.__rr, self.__angpts) + + def set_density(self, density): + """ + Set the density at the grid points + """ + self.__edens = density.reshape((len(self.__rr), len(self.__angpts))) + + def set_gradient(self, gradient): + """ + Set the density at the grid points + """ + self.__gradient = gradient.reshape((len(self.__rr), len(self.__angpts), 3)) + + def get_local_hartree_potential(self): + """ + Returns the local Hartree potential for only this atomic grid and + not considering any interactions of the other atomic grids + + Note: + * This function should only be used after htpot has been generated + using either calculate_coulomb_energy() or variants of this function + * This function should only be used for educational purposes + """ + return self.__htpot.flatten() + + def build_hartree_potential(self): + """ + Construct the Hartree potential at the grid points including cubic + interpolation to obtain the Hartree potential expansion coefficients + at other positions + """ + # for each radial grid point, construct the set of spherical harmonic + # coefficients that describes the density in that radial shell + self.__calculate_density_coefficients() + + # from the radial coefficients (rho_lm), construct the spherical + # harmonic coefficients that describe the Hartree potential (ulm) + # in that radial shell + self.__calculate_hartree_potential_coefficients() + + # build cubic interpolation + self.__ulm_splines = [] + rr = np.flip(self.__rr) + for i in range((self.__lmax+1)**2): + self.__ulm_splines.append(CubicSpline(rr, np.flip(self.__ulm[:,i]))) + + def get_ulm(self): + """ + Get Hartree potential coefficients + """ + return self.__ulm + + def get_rho_lm(self): + """ + Get electron density coefficients + """ + return self.__rho_lm + + def get_ylm(self): + """ + Get values for the spherical harmonics + """ + return self.__ylm + + def calculate_interpolated_ulm(self, rr): + """ + Calculate interpolated Hartree potential expansion coefficients + for all points r in the array rr + """ + ulmintp = np.zeros((len(self.__ulm_splines), len(rr))) + + # loop over lm pairs and interpolate value for Ulm at points rr + for i in range((self.__lmax+1)**2): + ulmintp[i,:] = self.__ulm_splines[i](rr) + + return ulmintp + + def set_molecular_weights(self, mweights): + """ + Set the molecular weights (Becke weights) at the grid points + """ + self.__mweights = mweights.reshape((len(self.__rr), len(self.__angpts))) + + def get_gridpoints(self): + """ + Return list of grid points + """ + return self.__gridpoints + + def get_weights(self): + """ + Get the weights on the atomic grid for the quadrature + """ + return self.__wgrid.flatten() + + def get_weights_angular_points(self): + """ + Get the weights for the angular points only + """ + return self.__wangpts + + def get_becke_weights(self): + """ + Return the fuzzy cell weights (Becke weights) + """ + return self.__mweights + + def get_density(self): + """ + Get the number of electrons for this atomic grid + """ + return self.__edens.flatten() + + def get_gradient(self): + """ + Get the gradient of the electron density for this atomic grid + """ + return self.__gradient.reshape((len(self.__rr) * len(self.__angpts), 3)) + + def get_gradient_squared(self): + """ + Get the modulus squared of the gradient for this atomic grid + """ + grad = self.get_gradient() + #return np.linalg.norm(grad, axis=1) + return np.einsum('ij,ij->i', grad, grad) + + def count_electrons(self): + """ + Sum over the electron density to obtain the number of electrons for + this atomic cell + """ + return np.einsum('ij,ij,ij', self.__edens, self.__wgrid, self.__mweights) + + def get_dfa_exchange(self): + """ + Get density functional approximation of the exchange energy for this atomic cell + """ + alpha = 2.0 / 3.0 + fac = -2.25 * alpha * np.power(0.75 / np.pi, 1.0 / 3.0) + return fac * np.einsum('ij,ij,ij', np.power(self.__edens, 4/3), self.__wgrid, self.__mweights) + + def get_dfa_kinetic(self): + """ + Get density functional approximation of the kinetic energy for this atomic cell + """ + Ckin = 3.0 / 40.0 * (3 / np.pi)**(2/3) * (2 * np.pi)**2 + return Ckin * np.einsum('ij,ij,ij', np.power(self.__edens, 5/3), self.__wgrid, self.__mweights) + + def get_dfa_nuclear_local(self): + """ + Get density functional approximation of the kinetic energy for this atomic cell + """ + npot = -self.__aidx / self.__rr # build nuclear attraction potential + return np.einsum('i,ij,ij,ij', npot, self.__edens, self.__wgrid, self.__mweights) + + def perform_spherical_harmonic_expansion(self, f): + """ + Calculate the expansion coefficients using a basis set composed of + spherical harmonics given a atomic orbital + + The expansion is performed over the function f which needs to be a matrix + of Nk x Nj elements where Nk is the number of radial points and Nj is + the number of angular points. + """ + return np.einsum('ij,kj,j,kj->ki', + self.__ylm, + f, + self.__wangpts, + self.__mweights) * 4.0 * np.pi + + def get_bragg_slater_radius(self) -> float: + """ + Get the Bragg-Slater radius + + Returns + ------- + float + Bragg-Slater radius + """ + return self.__rm + + def __calculate_density_coefficients(self): + """ + Calculate the expansion coefficients using a basis set composed of + spherical harmonics given a atomic orbital + + Note that we explicitly multiply with the Becke weights so that we + only integrate over the electron density assigned for this fuzzy cell. + That means that in any back-transformation, we should explicitly not + multiply with the Becke weights + """ + self.__rho_lm = self.perform_spherical_harmonic_expansion(self.__edens) + + def __calculate_hartree_potential_coefficients(self): + """ + Calculate the Hartree potential using the finite difference approximation + """ + # construct base matrix + Mbase = self.__build_finite_difference_matrix(self.__rr, 0.5 * self.__rm) + self.__ulm = np.zeros_like(self.__rho_lm) + + lmctr = 0 # indexing counter + for l in range(0, self.__lmax+1): + for m in range(-l, l+1): + M = Mbase.copy() # create copy of base matrix + N = M.shape[0] # grid size + 2 + b = np.zeros(N) + + for i in range(0, N): + # set first element for first lm value + if i == 0: + if l == 0: + b[i] = np.sqrt(4. * np.pi) * self.count_electrons() + continue + + # skip last element + if i == N-1: + continue + + M[i,i] -= l * (l+1) / self.__rr[i-1]**2 + b[i] = -4.0 * np.pi * self.__rr[i-1] * self.__rho_lm[i-1, lmctr] + + self.__ulm[:,lmctr] = np.linalg.solve(M,b)[1:-1] + + # increment lm counter + lmctr += 1 + + def calculate_nuclear_attraction(self): + """ + Calculate the nuclear attraction given the electron density + """ + npot = -self.__aidx / self.__rr # build nuclear attraction potential + return np.einsum('ij,i,ij', self.__wgrid, npot, self.__edens) + + def calculate_coulomb_energy(self): + """ + Calculate the electron-electron repulsion for this atomic grid + + This function is only for educational purposes, for any evaluation + of the coulomb energy, extensive interpolation over all the atomic + grids is necessary + """ + # build Hartree potential + pot = 1.0 / self.__rr + self.__htpot = np.einsum('ij,jk,i,ik->ik', self.__ulm, self.__ylm, pot, self.__edens) + + return np.einsum('ij,ij', self.__htpot, self.__wgrid) + + def calculate_coulomb_energy_interpolation(self): + """ + Calculate the electron-electron repulsion for this atomic grid + using cubic spline interpolation + + This function is only for educational purposes, for any evaluation + of the coulomb energy, extensive interpolation over all the atomic + grids is necessary + """ + # build Hartree potential + pot = 1.0 / self.__rr + ulmintp = self.calculate_interpolated_ulm(self.__rr).transpose() + self.__htpot = np.einsum('ij,jk,i,ik->ik', ulmintp, self.__ylm, pot, self.__edens) + + return np.einsum('ij,ij', self.__htpot, self.__wgrid) + + def __build_weight_grid(self): + """ + Calculates the weight factors on the grid and also builds the list + of coordinates of the grid points + """ + radial_weights = np.array([wgc * r**2 for wgc,r in zip(self.__wgcs, self.__rr)]) + self.__wgrid = np.outer(4.0 * np.pi * radial_weights, self.__wangpts) + + self.__gridpoints = [] + for rc,r in enumerate(self.__rr): + [self.__gridpoints.append(r * angpt + self.__cp) for angpt in self.__angpts] + + def __load_lebedev_coefficients(self, order:int): + """ + Load Lebedev coefficients from data file and store these as a class + variable in a dictionary. Each set of angular points is stored as + a Nx4 matrix wherein the first three per row correspond to the + position on the unit sphere and the last row to the weight as used + in the Lebedev integration. + """ + ag = AngularGrid() + return ag.get_coefficients(order) + + def __build_chebychev_grid(self): + """ + Build Chebychev grid of concentric spheres + """ + f = np.pi / float(self.__nshells+1) + z = np.array(range(1, self.__nshells+1)) + x = np.cos(f * z) + rm = self.__rm * 0.5 + + self.__rr = rm * (1.0 + x) / (1.0 - x) + self.__wgcs = f * np.sin(f * z)**2 / np.sqrt(1.0 - x**2) * 2.0 * rm / (1.0 - x)**2 + + def __build_ylm(self): + """ + Calculate the value for spherical harmonics at the angular points + """ + self.__ylm = np.zeros(((self.__lmax+1)**2, len(self.__angpts))) + lmctr = 0 + for l in range(0, self.__lmax+1): + for m in range(-l, l+1): + ylm = [spherical_harmonic(l, m, p[0], p[1]) for p in self.__usangles] + self.__ylm[lmctr,:] = np.array(ylm) + lmctr += 1 + + def __build_finite_difference_matrix(self, r, rm): + """ + Construct the finite difference matrix to solve for Ulm + + r : vector of distances + rm: half of the Bragg-Slater radius + """ + N = len(r) + A = np.zeros((N+2, N+2)) + h = 1.0 / float(N+1) + + for i in range(0, N+2): + c1 = 0.0 + c2 = 0.0 + + if i > 0 and i < N+1: + c1 = self.__dzdrsq(r[i-1], rm) + c2 = self.__d2zdr2(r[i-1], rm) + + if i == 0: + A[0,0] = 1.0 + continue + + if i == 1: + c1 /= 12.0 * h * h + c2 /= 12.0 * h + A[i,0] = 11.0 * c1 -3.0 * c2 + A[i,1] = -20.0 * c1 - 10.0 * c2 + A[i,2] = 6.0 * c1 + 18.0 * c2 + A[i,3] = 4.0 * c1 - 6.0 * c2 + A[i,4] = -1.0 * c1 + 1.0 * c2 + continue + + if i == 2: + c1 /= 12.0 * h * h + c2 /= 60.0 * h + A[i,0] = -1.0 * c1 + 3.0 * c2 + A[i,1] = 16.0 * c1 - 30.0 * c2 + A[i,2] = -30.0 * c1 - 20.0 * c2 + A[i,3] = 16.0 * c1 + 60.0 * c2 + A[i,4] = -1.0 * c1 - 15.0 * c2 + A[i,5] = 0.0 * c1 + 2.0 * c2 + continue + + if i == N-1: + c1 /= 12.0 * h * h + c2 /= 60.0 * h + A[i,N-4] = 0.0 * c1 - 2.0 * c2 + A[i,N-3] = -1.0 * c1 + 15.0 * c2 + A[i,N-2] = 16.0 * c1 - 60.0 * c2 + A[i,N-1] = -30.0 * c1 + 20.0 * c2 + A[i,N] = 16.0 * c1 + 30.0 * c2 + A[i,N+1] = -1.0 * c1 - 3.0 * c2 + continue + + if i == N: + c1 /= 12.0 * h * h + c2 /= 12.0 * h + A[i,N-3] = -1.0 * c1 - 1.0 * c2 + A[i,N-2] = 4.0 * c1 + 6.0 * c2 + A[i,N-1] = 6.0 * c1 - 18.0 * c2 + A[i,N] = -20.0 * c1 + 10.0 * c2 + A[i,N+1] = 11.0 * c1 + 3.0 * c2 + continue + + if i == N+1: + A[i,i] = 1.0 + continue + + c1 /= 180.0 * h * h + c2 /= 60.0 * h + A[i,i-3] = 2.0 * c1 - 1.0 * c2 + A[i,i-2] = -27.0 * c1 + 9.0 * c2 + A[i,i-1] = 270.0 * c1 - 45.0 * c2 + A[i,i] = -490.0 * c1 + A[i,i+1] = 270.0 * c1 + 45.0 * c2 + A[i,i+2] = -27.0 * c1 - 9.0 * c2 + A[i,i+3] = 2.0 * c1 + 1.0 * c2 + + return A + + def __set_bragg_slater_radius(self): + """ + Set the Bragg-Slater radius for the atom + """ + # note that the bragg_slater.BSRADII list starts with hydrogen at + # index 0 whereas aidx variable sets Hydrogen at index 1, Helium at + # index 2, and so on + self.__rm = bragg_slater.BSRADII[self.__aidx-1] + + def __d2zdr2(self, r, rm): + """ + Calculate second derivative of z-grid towards (regular) r-grid + """ + nom = rm * rm * (rm + 3.0 * r) + denom = 2.0 * np.pi * (rm * r)**(3/2) * (rm + r) ** 2 + return nom / denom + + def __dzdrsq(self, r, rm): + """ + Calculate square of first derivative of z-grid towards (regular) r-grid + """ + return rm / (np.pi * np.pi * r * (rm + r) * (rm + r)) \ No newline at end of file diff --git a/pydft/bragg_slater.py b/pydft/bragg_slater.py new file mode 100644 index 0000000..6ba23d7 --- /dev/null +++ b/pydft/bragg_slater.py @@ -0,0 +1,18 @@ +# -*- coding: utf-8 -*- + +# Bragg-Slater radii in angstrom +radii = [0.35, 2.00, 1.45, 1.05, 0.85, 0.70, 0.65, 0.60, 0.50, + 2.25, 1.80, 1.50, 1.25, 1.10, 1.00, 1.00, 1.00, 2.50, + 2.20, 1.80, 1.60, 1.40, 1.35, 1.40, 1.40, 1.40, 1.35, + 1.35, 1.35, 1.35, 1.30, 1.25, 1.15, 1.15, 1.15, 2.75, + 2.35, 2.00, 1.80, 1.55, 1.45, 1.45, 1.35, 1.30, 1.35, + 1.40, 1.60, 1.55, 1.55, 1.45, 1.45, 1.40, 1.40, 3.00, + 2.60, 2.15, 1.95, 1.85, 1.85, 1.85, 1.85, 1.85, 1.85, + 1.80, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.55, + 1.45, 1.35, 1.35, 1.30, 1.35, 1.35, 1.35, 1.50, 1.90, + 1.80, 1.60, 1.90, 1.65, 3.25, 2.80, 2.15, 1.95, 1.80, + 1.80, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75, + 1.75, 1.75, 1.75, 1.75, 1.55, 1.55] + +ANGSTROM2BOHR = 1.88973 +BSRADII = [radius * ANGSTROM2BOHR for radius in radii] \ No newline at end of file diff --git a/pydft/dft.py b/pydft/dft.py new file mode 100644 index 0000000..f526c8a --- /dev/null +++ b/pydft/dft.py @@ -0,0 +1,437 @@ +# -*- coding: utf-8 -*- + +from .moleculargrid import MolecularGrid +from pyqint import PyQInt, Molecule, cgf +import numpy as np +import time +from packaging import version +from copy import deepcopy + +# couple of hardcoded variables for the DIIS algorithm +SUBSPACE_LENGTH = 3 +SUBSPACE_START = 4 + +class DFT(): + def __init__(self, + mol:Molecule, + basis:str|list[cgf] = 'sto3g', + functional:str = 'svwn5', + nshells:int = 32, + nangpts:int = 110, + lmax:int = 8, + verbose:bool = False): + """ + Constructs the DFT class + + Parameters + ---------- + mol : Molecule + molecule + basis : str, optional + basis set, can be either a string or a list of cgf objects, by default 'sto3g' + functional : str, optional + exchange-correlation function, valid options are :code:`swn5` + and :code:`pbe`, by default 'svwn5' + nshells : int, optional + number of radial shells, by default 32 + nangpts : int, optional + number of angular sampling points, by default 110 + lmax : int, optional + maximum value of l in the spherical harmonic expansion, by default 8 + verbose : bool, optional + whether to provide verbose output, by default False + """ + self.__mol = mol + self.__integrator = PyQInt() + self.__basis = basis + self.__verbose = verbose + self.__time_stats = {} + self.__itermax = 100 + self.__nshells = nshells + self.__nangpts = nangpts + self.__lmax = lmax + self.__functional = functional + + def get_data(self) -> dict: + """ + Get relevant memory objects, only valid after successful SCF calculation + + Returns + ------- + dict + Dictionary containing relevant memory objects (see below) + + Notes + ----- + The dictionary contains the following elements: + + * :code:`S`: overlap matrix + * :code:`T`: kinetic energy matrix + * :code:`V`: nuclear attraction matrix + * :code:`C`: coefficient matrix + * :code:`J`: Hartree matrix + * :code:`P`: density matrix + * :code:`XC`: exchange-correlation matrix + * :code:`F`: Fock-matrix + * :code:`Exc`: exchange-correlation energy + * :code:`Ex`: exchange energy + * :code:`Ec`: correlation energy + * :code:`energies`: all energies + * :code:`energy`: last energy + * :code:`orbc`: coeffient matrix (duplicate) + * :code:`orbe`: molecular orbital eigenvalues + * :code:`enucrep`: electrostatic repulsion of the nuclei + + """ + data = { + 'S' : self.__S, # overlap matrix + 'T' : self.__T, # kinetic energy matrix + 'V' : self.__V, # nuclear attraction matrix + 'C' : self.__C, # coefficient matrix + 'J' : self.__J, # Hartree matrix + 'P' : self.__P, # density matrix + 'XC' : self.__XC, # exchange-correlation matrix + 'F' : self.__F, # Fock-matrix + 'Exc': self.__Exc, # exchange-correlation energy + 'Ex': self.__Ex, # exchange energy + 'Ec': self.__Ec, # correlation energy + 'energies': self.__energies, # all energy + 'energy': self.__energies[-1], # last energy + 'orbc': self.__C, # coeffient matrix (duplicate) + 'orbe': self.__e, # molecular orbital eigenvalues + 'enucrep': self.__enuc, # electrostatic repulsion of the nuclei + } + + return data + + def get_molgrid_copy(self) -> MolecularGrid: + """ + Get a copy of the underlying molgrid object + + Returns: + MolecularGrid: copy of the :class:`MolecularGrid` object + """ + return deepcopy(self.__molgrid) + + def get_density_at_points(self, spoints:np.ndarray) -> np.ndarray: + """ + Get the electron density at provided points + + Parameters + ---------- + spoints : ndarray + Set of sampling points (:math:`X \\times 3` array) + + Returns + ------- + ndarray + Electron density scalar field (:math:`N \\times 1` array) + """ + if len(spoints.shape) != 2: + raise Exception('Grid points need to be supplied as a Nx3 array.') + + return self.__molgrid.get_density_at_points(spoints, self.__P) + + def get_gradient_at_points(self, spoints:np.ndarray) -> np.ndarray: + """ + Get the electron density gradient at provided points + + Parameters + ---------- + spoints : ndarray + Set of sampling points (:math:`X \\times 3` array) + + Returns + ------- + ndarray + Electron density gradient vector field (:math:`N \\times 3` array) + """ + if len(spoints.shape) != 2: + raise Exception('Grid points need to be supplied as a Nx3 array.') + + return self.__molgrid.get_gradient_at_points(spoints, self.__P) + + def scf(self, tol:float=1e-5) -> float: + """ + Perform the self-consistent field procedure + + Parameters + ---------- + tol : float, optional + electronic convergence criterion, by default 1e-5 + + Returns + ------- + float + total electronic energy (in Hartrees) + """ + # construct stagnant matrices + self.__setup() + + # create empty P matrix as initial guess + self.__P = np.zeros_like(self.__S) + + # start SCF iterative procedure + nitfin = 0 + for niter in range(0, self.__itermax): + start = time.time() + energy = self.__iterate(niter, + giis=True if nitfin == 0 else False, + mix=0.9) + self.__energies.append(energy) + stop = time.time() + itertime = stop - start + self.__time_stats['iterations'].append(itertime) + + if self.__verbose: + print('%03i | Energy: %12.6f | %0.4f ms' % (niter+1, energy, itertime)) + + if niter > 2: + ediff = np.abs(energy - self.__energies[-2]) + if ediff < tol: + # terminate giis self-convergence and continue with mixing + nitfin += 1 + + if nitfin < 3: + continue + + # terminate self-convergence cycle + if self.__verbose: + print("Stopping SCF cycle, convergence reached.") + + # update density matrix from last found coefficient matrix + self.__P = self.__calculate_P() + break + + return energy + + def get_construction_times(self) -> dict: + """ + Return construct times dictionary from MolecularGrid to get insights + into the time to construct particular properties + + Returns + ------- + dict + dictionary of construction times for the MolecularGrid objects + """ + return self.__molgrid.construct_times + + def __iterate(self, niter, giis=True, mix=0.9): + """ + Perform single-step iteration + """ + # calculate J and XC matrices based on the current electron + # density estimate as captured in the density matrix P + if niter > SUBSPACE_START and giis: + try: + diis_coeff = self.__calculate_diis_coefficients(self.__evs_diis) + self.__F = self.__extrapolate_fock_from_diis_coefficients(self.__fmats_diis, diis_coeff) + self.__Fprime = self.__X.transpose().dot(self.__F).dot(self.__X) + self.__e, self.__Cprime = np.linalg.eigh(self.__Fprime) + self.__C = self.__X.dot(self.__Cprime) + self.__P = self.__calculate_P() + except np.linalg.LinAlgError: # set giis to False if not working + giis = False + + # calculate J and XC matrices based on the current electron + # density estimate as captured in the density matrix P + if np.any(self.__P): + self.__molgrid.build_density(self.__P, normalize=True) + self.__J = self.__calculate_J() + self.__XC, self.__Exc = self.__calculate_XC() + + # calculate Fock matrix + self.__F = self.__H + self.__J + self.__XC + + # perform unitary transformation on Fock matrix + self.__Fprime = self.__X.transpose().dot(self.__F).dot(self.__X) + + # diagonalize Fock matrix + try: + self.__e, self.__Cprime = np.linalg.eigh(self.__Fprime) + except np.linalg.LinAlgError: + print('Error: eigenvalue convergence failed in iteration: %i' % niter) + print('F:', self.__Fprime) + print('H:', self.__H) + print('J:', self.__J) + print('XC:', self.__XC) + print('P:', self.__P) + raise np.linalg.LinAlgError + + # back-transform + self.__C = self.__X.dot(self.__Cprime) + + # calculate total electronic energy + M = self.__T + self.__V + 0.5 * self.__J + P = self.__calculate_P() + energy = np.einsum('ji,ij', P, M) + self.__Exc + self.__enuc + + # for the first few iterations, build a new density + # matrix from the coefficients, else, resort to the DIIS + # algorithm + if niter <= SUBSPACE_START or not giis: + P = self.__calculate_P() + self.__P = (1.0 - mix) * self.__P + mix * P # use linear mixing + + # calculate DIIS coefficients + e = (self.__F.dot(self.__P.dot(self.__S)) - \ + self.__S.dot(self.__P.dot(self.__F))).flatten() # calculate error vector + self.__enorm = np.linalg.norm(e) # store error vector norm + self.__fmats_diis.append(self.__F) # add Fock matrix to list + self.__pmat_diis.append(self.__P) # add density matrix to list + self.__evs_diis.append(e) + + # prune size of the old Fock, density and error vector lists + # only SUBSPACE_LENGTH iterations are used to guess the new + # solution + if len(self.__fmats_diis) > SUBSPACE_LENGTH: + self.__fmats_diis = self.__fmats_diis[-SUBSPACE_LENGTH:] + self.__pmat_diis = self.__pmat_diis[-SUBSPACE_LENGTH:] + self.__evs_diis = self.__evs_diis[-SUBSPACE_LENGTH:] + + return energy + + def __setup(self): + """ + Construct the bare classes and matrices necessary to start the + self-consistent-field procedure + """ + # construct basis functions and nuclei + if issubclass(type(self.__basis), str): # if a basis set name is given + self.__cgfs, self.__nuclei = self.__mol.build_basis(self.__basis) + else: # either assume a list of CGFs objects is given + self.__cgfs = self.__basis + self.__nuclei = self.__mol.get_nuclei() + + # build molecular grid + self.__molgrid = MolecularGrid(self.__nuclei, + self.__cgfs, + nshells=self.__nshells, + nangpts=self.__nangpts, + lmax=self.__lmax, + functional=self.__functional) + self.__molgrid.initialize() # molecular grid uses late initialization + + # build one-electron matrices; because these matrices are Hermetian, + # we only have to evaluate the upper half and then simply copy the + # upper half to the lower half + N = len(self.__cgfs) + self.__S = np.zeros((N,N)) + self.__T = np.zeros_like(self.__S) + self.__V = np.zeros_like(self.__S) + for j in range(0,N): + cgf1 = self.__cgfs[j] + for i in range(j,N): + cgf2 = self.__cgfs[i] + self.__S[j,i] = self.__integrator.overlap(cgf1, cgf2) + self.__T[j,i] = self.__integrator.kinetic(cgf1, cgf2) + + # in principle, the nuclear attraction could also be directly + # obtained from the molecular grid, but analytical evaluation + # is faster + for k,nucleus in enumerate(self.__nuclei): + vjik = self.__integrator.nuclear(cgf1, cgf2, nucleus[0], nucleus[1]) + self.__V[j,i] += vjik + + # copy upper triangle elements to lower triangle + if i != j: + self.__S[i,j] = self.__S[j,i] + self.__T[i,j] = self.__T[j,i] + self.__V[i,j] = self.__V[j,i] + + # build single-electron matrix + self.__H = self.__T + self.__V + + # diagonalize S and use it to construct the unitary transformation + # matrix that orthonormalizes the basis set + s, U = np.linalg.eigh(self.__S) + self.__X = U.dot(np.diag(1.0/np.sqrt(s))) + + # create empty matrices for the coulomb and exchange-correlation + # energies + self.__J = np.zeros_like(self.__S) + self.__XC = np.zeros_like(self.__S) + self.__Exc = 0.0 + + # create zero density matrix + self.__P = np.zeros_like(self.__S) + + # calculate nuclear repulsion + self.__enuc = 0.0 + for i in range(0, len(self.__nuclei)): + for j in range(i+1, len(self.__nuclei)): + r = np.linalg.norm(self.__nuclei[i][0] - self.__nuclei[j][0]) + self.__enuc += self.__nuclei[i][1] * self.__nuclei[j][1] / r + + # build containers to store per-iteration data + self.__energies = [] + self.__time_stats['iterations'] = [] + self.__fmats_diis = [] + self.__pmat_diis = [] + self.__evs_diis = [] + + def __calculate_J(self): + """ + Calculate the coulombic interaction matrix using the + molecular grid + """ + return self.__molgrid.calculate_coulombic_matrix() + + def __calculate_P(self): + """ + Calculate density matrix from current coefficient matrix + """ + N = len(self.__cgfs) + P = np.zeros_like(self.__S) + nelec = np.sum([nucleus[1] for nucleus in self.__nuclei]) + for i in range(0,N): + for j in range(0,N): + for k in range(0,int(nelec/2)): + P[i,j] += 2.0 * self.__C[i,k] * self.__C[j,k] + + return P + + def __calculate_XC(self): + """ + Calculate the exchange-correlation matrix and the + exchange-correlation energy + """ + X, self.__Ex = self.__molgrid.calculate_exchange() + C, self.__Ec = self.__molgrid.calculate_correlation() + + return X+C, self.__Ex + self.__Ec + + def __calculate_diis_coefficients(self, evs_diis): + """ + Calculate the DIIS coefficients + """ + B = np.zeros((len(evs_diis)+1, len(evs_diis)+1)) + B[-1,:] = -1 + B[:,-1] = -1 + B[-1,-1]= 0 + + rhs = np.zeros((len(evs_diis)+1, 1)) + rhs[-1,-1] = -1 + + for i in range(len(evs_diis)): + for j in range(i+1): + B[i,j] = np.dot(evs_diis[i].transpose(), evs_diis[j]) + B[j,i] = B[i,j] + + *diis_coeff, _ = np.linalg.solve(B,rhs) + + return diis_coeff + + def __extrapolate_fock_from_diis_coefficients(self, fmats_diis, diis_coeff): + """ + Extrapolate the Fock matrix from the DIIS coefficients + """ + norbs = fmats_diis[-1].shape[0] + fguess = np.zeros((norbs,norbs)) + + for i in range(len(fmats_diis)): + fguess += fmats_diis[i]*diis_coeff[i] + + return fguess + \ No newline at end of file diff --git a/pydft/moleculargrid.py b/pydft/moleculargrid.py new file mode 100644 index 0000000..ec437ed --- /dev/null +++ b/pydft/moleculargrid.py @@ -0,0 +1,919 @@ +# -*- coding: utf-8 -*- + +import numpy as np +from .atomicgrid import AtomicGrid +from . import bragg_slater +from .spherical_harmonics import spherical_harmonic +from pyqint import PyQInt +import pyqint as pq +import time +from .xcfunctionals import Functionals + +class MolecularGrid: + def __init__(self, + atoms:list, + cgfs:list[pq.cgf], + nshells:int=32, + nangpts:int=110, + lmax:int=8, + functional:str='svwn5'): + """ + Construct MolecularGrid + + Parameters + ---------- + atoms : list + list of atoms and their charges + cgfs : list[pq.cgf] + list of contracted Gaussian functions (basis set) + nshells : int, optional + number of radial shells, by default 32 + nangpts : int, optional + number of angular sampling points per shell, by default 110 + lmax : int, optional + maximum value for l for projection of spherical harmonics, by default 8 + functional : str, optional + exchange-correlation functional, by default 'svwn5' + """ + # keep track of build times + self.construct_times = {} + + # set properties + self.__atoms = atoms + self.__nelec = np.sum([nuc[1] for nuc in self.__atoms]) + self.__lmax = lmax + self.__nshells = nshells + self.__nangpts = nangpts + self.__basis = cgfs + self.__functionals = Functionals(functional) + self.__is_initialized = False + + def initialize(self): + """ + Initialize the molecular grid + + This is a relatively expensive procedure and thus the user can + delay the initialization of the molecular grid until. + """ + if self.__is_initialized: + return + + self.__build_molecular_grid() + self.__build_amplitudes() + self.__is_initialized = True + + def build_density(self, P:np.ndarray, normalize:bool=False): + """ + Build the electron density from the density matrix + + Parameters + ---------- + P : np.ndarray + density matrix + normalize : bool, optional + whether to normalize the density matrix based on the total number + of electrons, by default False + """ + # this function requires the amplitudes and molecular grid + # to be built + self.initialize() + + + # calculate densities at each grid point + self.__densities = np.einsum('ijk,jl,ilk->ik', + self.__amplitudes, + P, + self.__amplitudes) + + # also build the gradient of the density + self.__gradients = 2.0 * np.einsum('ijk,jl,ilkm->ikm', + self.__amplitudes, + P, + self.__ampgrads) + + # perform optional normalization + if normalize: + nelec = np.sum(self.__mgw * self.__densities.flatten()) + cn = nelec / self.__nelec # normalization constant + self.__densities *= cn + self.__gradients *= cn + + # and place the densities back into the atomic grids + for i,atgrid in enumerate(self.__atomgrids): + atgrid.set_density(self.__densities[i,:]) + atgrid.set_gradient(self.__gradients[i,:,:]) + atgrid.build_hartree_potential() + + def get_becke_weights(self) -> np.ndarray: + """ + Get the Becke coefficients + + Returns + ------- + np.array + :math:`(N \\times G \\times 3)` array with :math:`N` number of atoms + and :math:`G` number of grid points + """ + return self.__mweights + + def get_grid_coordinates(self) -> list[np.array]: + """ + Get the grid coordinates + + Returns + ------- + list[np.array] + list over atoms with for each atom a :math:`(G \\times 1)` array of + electron densities where :math:`G` is the number of grid points + per atom + """ + return self.__gridpoints + + def get_hartree_potential(self) -> np.array: + """ + Get the Hartree potential at each grid cell + + Returns + ------- + np.array + :math:`(G \\times 1)` array with :math:`G` number of molecular + grid points + """ + return self.__ugpts + + def get_densities(self) -> list[np.array]: + """ + Get the densities at the grid points + + Returns + ------- + list[np.array] + list over atoms with for each atom a :math:`(G \\times 1)` array of + electron densities where :math:`G` is the number of grid points + per atom + """ + return [atgrid.get_density() for atgrid in self.__atomgrids] + + def get_gradients(self) -> list[np.array]: + """ + Get the gradient magnitude of the electron density field + + Returns + ------- + list[np.array] + list over atoms with for each atom a :math:`(G \\times 1)` array of + electron gradient magnitudes where :math:`G` is the number of grid + points per atom + """ + return [atgrid.get_gradient() for atgrid in self.__atomgrids] + + def get_rho_lm_atoms(self) -> np.array: + """ + Get the electron density projected onto spherical harmonics per atom + + Returns + ------- + np.array + :math:`\\rho_{n,lm}` where :math:`n` loops over the atoms and the + pairs :math:`lm` over the spherical harmonics + """ + rho_lm = np.ndarray((len(self.__atoms), self.__nshells, (self.__lmax+1)**2)) + + for i,at in enumerate(self.__atomgrids): + rho_lm[i,:,:] = (at.get_rho_lm().T * at.get_radial_grid()**2).T + + return rho_lm + + def get_atomic_grid(self, i:int) -> AtomicGrid: + """ + Get an atomic grid of atom :math:`i` + + Parameters + ---------- + i : int + atom index + + Returns + ------- + AtomicGrid + :class:`pydft.AtomicGrid` of atom :math:`i` + """ + return self.__atomgrids[i] + + def count_electrons(self) -> float: + """ + Count number of electrons in the system + + Returns + ------- + float + Total number of electrons + """ + return np.sum([atgrid.count_electrons() for atgrid in self.__atomgrids]) + + def count_electrons_from_rho_lm(self) -> float: + """ + Count number of electrons from the spherical harmonic expansion; this + function mainly acts to verify that the spherical harmonic expansion + works correctly + + Returns + ------- + float + Total number of electrons + + Notes + ----- + The Becke weights are already used when generating the density + coefficients and hence for the back-transformation, we should *not* + multiply with the Becke weights + """ + # get spherical harmonic expansion coefficients per atomic grid + ngpts = self.__nshells * self.__nangpts # gridpoints per atom + rho_lm = np.ndarray((len(self.__atoms), self.__nshells, (self.__lmax+1)**2)) + + for i,at in enumerate(self.__atomgrids): + rho_lm[i,:,:] = at.get_rho_lm() + + # reshape the rho_lm array such that spherical harmonic coefficients + # are set for each grid point and only store the coefficients for + # each atomic grid associated with its own atom + rho_lm_gpts = np.ndarray((len(self.__atoms), (self.__lmax+1)**2, ngpts)) + ylm = np.zeros_like(rho_lm_gpts) + for i,at in enumerate(self.__atomgrids): + for j in range(0, (self.__lmax+1)**2): + rho_lm_gpts[i,j,:] = np.outer(rho_lm[i,:,j],np.ones(self.__nangpts)).flatten() + ylm[i,j,:] = self.__ylmgpts[i,j,i*ngpts:(i+1)*ngpts] + + return np.einsum('ijk,ijk,ik', rho_lm_gpts, + ylm, + np.array([an.get_weights() for an in self.__atomgrids])) + + def calculate_dfa_nuclear_attraction_local(self) -> float: + """ + Get density functional approximation for the nuclear attraction energy + using only the local potential for each atomic grid and ignoring the + influence of atomic attraction due to other atoms + + Returns + ------- + float + Approximate nuclear attraction energy (Hartree) + + Notes + ----- + This function will of course not give an accurate result and is only + meant for educational purposes + """ + return np.sum([atgrid.get_dfa_nuclear_local() for atgrid in self.__atomgrids]) + + def calculate_dfa_nuclear_attraction_full(self) -> np.array: + """ + Get density functional approximation for the nuclear attraction energy + + Returns + ------- + float + Total nuclear attraction energy (Hartree) + """ + return np.einsum('i,i,i', self.__npot, + np.array([an.get_density() for an in self.__atomgrids]).flatten(), + self.__mgw) + + def calculate_dfa_coulomb(self) -> float: + """ + Get density functional approximation for the coulombic repulsion + between electrons + + Returns + ------- + float + Total coulombic repulsion (Hartree) + """ + # for each grid point, collect the spherical harmonic expansion coefficients + # of the hartree potential for all the atoms + ulmgpts = np.ndarray((len(self.__atoms), (self.__lmax+1)**2, len(self.__rgridpoints[0]))) + for i,at in enumerate(self.__atomgrids): + ulmgpts[i,:,:] = at.calculate_interpolated_ulm(self.__rgridpoints[i]) + + # for each grid point determine the Hartree potential from the + # spherical harmonic coefficients with respect to each atomic center + self.__ugpts = np.einsum('ijk,ik,ijk->k', ulmgpts, self.__rigridpoints, self.__ylmgpts) + + return 0.5 * np.einsum('i,i,i', self.__ugpts, + np.array([an.get_density() for an in self.__atomgrids]).flatten(), + self.__mgw) + + def calculate_dfa_coulomb_no_interpolation(self) -> float: + """ + Calculate the Coulombic self-repulsion in the limit of the sum of atomic centers, + without any interpolation of the contribution of the other atoms + + Returns + ------- + float + Approximate coulombic repulsion (Hartree) + """ + return 0.5 * np.sum([at.calculate_coulomb_energy() for at in self.__atomgrids]) + + def calculate_coulombic_matrix(self) -> np.array: + r""" + Build the coulomb matrix :math:`\mathbf{J}` + + Returns + ------- + np.array + Coulomb matrix :math:`\mathbf{J}` + """ + # for each grid point, collect the spherical harmonic expansion coefficients + # of the hartree potential for all the atoms + ulmgpts = np.ndarray((len(self.__atoms), (self.__lmax+1)**2, len(self.__rgridpoints[0]))) + for i,at in enumerate(self.__atomgrids): + ulmgpts[i,:,:] = at.calculate_interpolated_ulm(self.__rgridpoints[i]) + + # for each grid point determine the Hartree potential from the + # spherical harmonic coefficients with respect to each atomic center + self.__ugpts = np.einsum('ijk,ik,ijk->k', ulmgpts, self.__rigridpoints, self.__ylmgpts) + + # construct coulombic repulsion matrix by integrating the interaction + # of the hartree potential with the basis function amplitudes + N = len(self.__basis) + J = np.zeros((N,N)) + for i in range(0, N): + for j in range(i, N): + J[i,j] = np.einsum('i,i,i,i', self.__ugpts, + self.__fullgrid_amplitudes[i,:], + self.__fullgrid_amplitudes[j,:], + self.__mgw) + if i != j: + J[j,i] = J[i,j] + + return J + + def calculate_dfa_kinetic(self) -> float: + """ + Get density functional approximation for the kinetic energy + + Returns + ------- + float + Total kinetic energy (Hartree) + """ + return np.sum([atgrid.get_dfa_kinetic() for atgrid in self.__atomgrids]) + + def calculate_dfa_exchange(self): + """ + Get density functional approximation for the exchange energy + + Returns + ------- + float + Total exchange energy (Hartree) + """ + return np.sum([atgrid.get_dfa_exchange() for atgrid in self.__atomgrids]) + + def calculate_exchange(self) -> (np.ndarray,float): + r""" + Build the exchange matrix and give the exchange energy for the molecule + + Returns + ------- + np.ndarray,float + Exchange matrix :math:`\mathbf{X}` and exchange energy + """ + + dens = np.array([atgrid.get_density() for atgrid in self.__atomgrids]).flatten() + if self.__functionals.is_gga(): + grad = np.array([atgrid.get_gradient_squared() for atgrid in self.__atomgrids]).flatten() + fx, vfx = self.__functionals.calc_x(dens, grad) + else: + fx, vfx = self.__functionals.calc_x(dens) + + # calculate exchange energy + ex = np.einsum('i,i,i', self.__mgw, fx, dens) + + # build matrix and pre-calculate weights + X = np.zeros((len(self.__basis), len(self.__basis))) + + # exchange parameters + for i in range(0, len(self.__basis)): + for j in range(i, len(self.__basis)): + X[i,j] = np.einsum('i,i,i,i', vfx, + self.__fullgrid_amplitudes[i,:], + self.__fullgrid_amplitudes[j,:], + self.__mgw) + + if i != j: + X[j,i] = X[i,j] + + return X, ex + + def calculate_correlation(self) -> (np.ndarray, float): + """ + Build the correlation matrix and give the correlation energy for the molecule + + Returns + ------- + np.ndarray,float + Correlation matrix and exchange energy + """ + dens = np.array([atgrid.get_density() for atgrid in self.__atomgrids]).flatten() + if self.__functionals.is_gga(): + grad = np.array([atgrid.get_gradient_squared() for atgrid in self.__atomgrids]).flatten() + fc, vfc = self.__functionals.calc_c(dens, grad) + else: + fc, vfc = self.__functionals.calc_c(dens) + + # calculate correlation energy + ec = np.einsum('i,i,i', self.__mgw, fc, dens) + + # build matrix and pre-calculate weights + C = np.zeros((len(self.__basis), len(self.__basis))) + + # exchange parameters + for i in range(0, len(self.__basis)): + for j in range(i, len(self.__basis)): + C[i,j] = np.einsum('i,i,i,i', vfc, + self.__fullgrid_amplitudes[i,:], + self.__fullgrid_amplitudes[j,:], + self.__mgw) + + if i != j: + C[j,i] = C[i,j] + + return C, ec + + def get_density_at_points(self, + spoints:np.ndarray, + P:np.ndarray) -> np.array: + """ + Calculate the electron density at the points as given by spoints + using density matrix P + + Parameters + ---------- + spoints : np.ndarray + Array of grid points :math:`(N \\times 3)` with :math:`N` the number + of grid points + P : np.ndarray + Density matrix :math:`(K \\times K)` with :math:`K` the number of + basis functions + + Returns + ------- + np.array + Electron density specified at grid points :math:`(N \\times 1)` + with :math:`N` the number of grid points + """ + # build the amplitudes at the specified points + amps = np.ndarray((len(self.__basis), len(spoints))) + for i,cgf in enumerate(self.__basis): + amps[i,:] = np.array([cgf.get_amp(p) for p in spoints]) + + dens = np.einsum('ik,ij,jk->k', amps, P, amps) + + return dens + + def get_amplitude_at_points(self, + spoints:np.ndarray, + c:np.ndarray) -> np.ndarray: + """ + Calculate the wave function amplitude at the points as given by + spoints using solution vector c + + Parameters + ---------- + spoints : np.ndarray + Array of grid points :math:`(N \\times 3)` with :math:`N` the number + of grid points + c : np.ndarray + Coefficient vector :math:`\\vec{c}` of a single molecular orbital + + Returns + ------- + np.array + Electron density specified at grid points :math:`(N \\times 1)` + with :math:`N` the number of grid points + """ + # build the amplitudes at the specified points + amps = np.ndarray((len(self.__basis), len(spoints))) + for i,cgf in enumerate(self.__basis): + amps[i,:] = np.array([cgf.get_amp(p) for p in spoints]) + + wfamp = np.einsum('ik,i->k', amps, c) + + return wfamp + + def get_gradient_at_points(self, + spoints:np.ndarray, + P:np.ndarray): + """ + Calculate the gradient of the electron density at the points as + given by spoints using density matrix :math:`\\mathbf{P}` + + Parameters + ---------- + spoints : np.ndarray + Array of grid points :math:`(N \\times 3)` with :math:`N` the number + of grid points + P : np.ndarray + Density matrix :math:`(K \\times K)` with :math:`K` the number of + basis functions + + Returns + ------- + np.array + Electron density gradients at grid points :math:`(N \\times 3)` + with :math:`N` the number of grid points + """ + # build the amplitudes at the specified points + amps = np.ndarray((len(self.__basis), len(spoints))) + grads = np.ndarray((len(self.__basis), len(spoints), 3)) + for i,cgf in enumerate(self.__basis): + amps[i,:] = np.array([cgf.get_amp(p) for p in spoints]) + grads[i,:,:] = np.array([cgf.get_grad(p) for p in spoints]) + + # + # note that the 'ij' derivative component is equal + # to the 'ji' derivative component + # + #t1 = np.einsum('ij,ikl,jk->kl', P, grads, amps) + #t2 = np.einsum('ij,jkl,ik->kl', P, grads, amps) + #return t1 + t2 + + return 2.0 * np.einsum('ij,ikl,jk->kl', P, grads, amps) + + def calculate_coulomb_potential_at_points(self, + pts:np.ndarray) -> np.ndarray: + r""" + Collect the spherical harmonic expansion coefficients of the hartree + potential for all the atoms at the grid points + + Parameters + ---------- + pts : np.ndarray + Array of grid points :math:`(N \\times 3)` with :math:`N` the number + of grid points + + Returns + ------- + np.array + Coulomb potential at grid points :math:`(\left[l_{\\textrm{max}}+1\\right]^{2} \\times N)` + with :math:`N` the number of grid points + """ + # for each grid point, + ulmgpts = np.zeros(((self.__lmax+1)**2, len(pts))) + for i,at in enumerate(self.__atomgrids): + rgridpts = np.linalg.norm(pts - self.__atoms[i][0], axis=1) + igpts = 1.0 / rgridpts + val = at.calculate_interpolated_ulm(rgridpts) + for j in range(0, len(val)): + val[j,:] *= igpts + ulmgpts += val + + return ulmgpts + + def get_exchange_potential_at_points(self, + pts:np.ndarray, + P:np.ndarray) -> np.ndarray: + """ + Get the exchange potential at points pts + + Parameters + ---------- + pts : np.ndarray + Array of grid points :math:`(N \\times 3)` with :math:`N` the number + of grid points + P : np.ndarray + Density matrix :math:`(K \\times K)` with :math:`K` the number of + basis functions + + Returns + ------- + np.ndarray + Exchange potential :math:`\\nu_{x}` at grid points :math:`(N \\times 1)` + with :math:`N` the number of grid points + """ + dens = self.get_density_at_points(pts, P) + if self.__functionals.is_gga(): + grad = np.linalg.norm(self.get_gradient_at_points(pts, P), axis=1) + fx, vfx = self.__functionals.calc_x(dens, grad) + else: + fx, vfx = self.__functionals.calc_x(dens) + + return vfx + + def get_correlation_potential_at_points(self, + pts:np.ndarray, + P:np.ndarray) -> np.ndarray: + """ + Get the correlation potential at points pts + + Parameters + ---------- + pts : np.ndarray + Array of grid points :math:`(N \\times 3)` with :math:`N` the number + of grid points + P : np.ndarray + Density matrix :math:`(K \\times K)` with :math:`K` the number of + basis functions + + Returns + ------- + np.ndarray + Correlation potential :math:`\\nu_{c}` at grid points :math:`(N \\times 1)` + with :math:`N` the number of grid points + """ + dens = self.get_density_at_points(pts, P) + if self.__functionals.is_gga(): + grad = np.linalg.norm(self.get_gradient_at_points(pts, P), axis=1) + fc, vfc = self.__functionals.calc_c(dens, grad) + else: + fc, vfc = self.__functionals.calc_c(dens) + + return vfc + + def get_spherical_harmonic_expansion_of_amplitude(self, + c:np.ndarray, + radial_factor=False) -> list[np.ndarray]: + """ + Get the spherical harmonic expansion representation of a wavefunction + amplitude using solution vector c + + Parameters + ---------- + c : np.ndarray + Wave function expansion coefficient + radial_factor : bool, optional + whether to multiply result by :math:`r^{2}`, by default False + + Returns + ------- + list[np.ndarray] + List of spherical harmonic expansions per atom stored as + :math:`(N_{r} \\times N_{lm})) arrays where :math:`N_{r}` is the + number of radial points and :math:`N_{lm}` the set of spherical + harmonics used in the expansion. + """ + she = [] + for at in self.__atomgrids: + atgpts = at.get_full_grid() + nr, na, _dummy = atgpts.shape + amps = self.get_amplitude_at_points(atgpts.reshape(-1,3), c) + she_at = at.perform_spherical_harmonic_expansion(amps.reshape((nr,na))) + + if radial_factor: + r2 = at.get_radial_grid()**2 + she_at = (she_at.T * r2).T + + she.append(she_at) + + return np.array(she) + + def calculate_weights_at_points(self, + points:np.ndarray, + k:int=3) -> np.ndarray: + """ + Custom function that (re-)calculates weights at given points + + Parameters + ---------- + points : np.ndarray + :math:`N \\times 3` array of grid points + k : int, optional + smoothing value, by default 3 + + Returns + ------- + np.ndarray + :math:`N_{a} \\times N` array with :math:`N_{a}` the number of + atoms and :math:`N` the number of grid points + """ + mweights = np.zeros((len(self.__atoms), len(points))) + smats = np.ndarray((len(self.__atoms), len(self.__atoms), len(points))) + + # loop over atoms + for i in range(0, len(self.__atoms)): + R1 = np.array(self.__atoms[i][0]) # position of atom 1 + rm1 = bragg_slater.BSRADII[self.__atoms[i][1]-1] # bragg-slater radius + for j in range(i+1, len(self.__atoms)): + R2 = np.array(self.__atoms[j][0]) # position of atom 2 + rm2 = bragg_slater.BSRADII[self.__atoms[j][1]-1] # bragg-slater radius + Rij = np.linalg.norm(R1 - R2) # distance between the two atoms + chi = rm1 / rm2 # fraction of bragg-slater radii + + # calculate the eliptical coordinate mu at the point r on the + # grid with respect to the two atoms + muv = (np.linalg.norm(points - R1, axis=1) - \ + np.linalg.norm(points - R2, axis=1)) / Rij + + # calculate the fuzzy cell coefficient between the two atoms i and j + smats[i,j,:] = self.__sk(self.__vij(muv, chi), k) + + # automatically calculate the for j,i by taking one minus the result + smats[j,i,:] = np.ones(len(smats[i,j,:])) - smats[i,j,:] + + # construct the cell functions for the points + Pn = np.ones((len(self.__atoms), len(points))) + for i in range(0, len(self.__atoms)): + for j in range(0, len(self.__atoms)): + if j != i: + Pn[i,:] = np.multiply(Pn[i,:], smats[i,j,:]) + + # calculate total of cell function for each atom at each grid point + Pt = np.einsum('ij->j', Pn) + + # and normalize each point so that the sum equals unity + for i in range(0, len(self.__atoms)): + mweights[i,:] = np.divide(Pn[i,:], Pt) # cell function for each atom + + return mweights + + def __build_molecular_grid(self): + """ + Build the molecular grid from the atomic grids + """ + + # build atomic grids + st = time.time() + self.__atomgrids = [] + for atom in self.__atoms: + self.__atomgrids.append(AtomicGrid(atom, + self.__nshells, + self.__nangpts, + self.__lmax) + ) + self.construct_times['atomic_grids'] = time.time() - st + + # assign to each gridpoint in the atomic grid a weight in the molecular + # grid + st = time.time() + self.__gridpoints = [] + self.__mweights = np.zeros((len(self.__atoms), + self.__nshells * int(self.__nangpts))) + + for g,atgrid in enumerate(self.__atomgrids): + grid = np.array(atgrid.get_gridpoints()) + self.__gridpoints.append(grid) + smats = np.ndarray((len(self.__atoms), len(self.__atoms), len(grid))) + for i in range(0, len(self.__atoms)): + R1 = np.array(self.__atoms[i][0]) # position of atom 1 + rm1 = bragg_slater.BSRADII[self.__atoms[i][1]-1] # bragg-slater radius + for j in range(i+1, len(self.__atoms)): + R2 = np.array(self.__atoms[j][0]) # position of atom 2 + rm2 = bragg_slater.BSRADII[self.__atoms[j][1]-1] # bragg-slater radius + Rij = np.linalg.norm(R1 - R2) # distance between the two atoms + chi = rm1 / rm2 # fraction of bragg-slater radii + + # calculate the confocal elliptical coordinate mu at the point r on the + # grid with respect to the two atoms + muv = (np.linalg.norm(grid - R1, axis=1) - \ + np.linalg.norm(grid - R2, axis=1)) / Rij + + # calculate the fuzzy cell coefficient between the two atoms i and j + smats[i,j,:] = self.__sk(self.__vij(muv, chi), 3) + + # automatically calculate the for j,i by taking one minus the result + smats[j,i,:] = np.ones(len(smats[i,j,:])) - smats[i,j,:] + + # for each atom, construct the cell function based on the product + # of the fuzzy cells coefficient for that cell with each other cell + Pn = np.ones((len(self.__atoms), len(grid))) + for i in range(0, len(self.__atoms)): + for j in range(0, len(self.__atoms)): + if j != i: + Pn[i,:] = np.multiply(Pn[i,:], smats[i,j,:]) + + # calculate total of cell function for each atom at each grid point + Pt = np.einsum('ij->j', Pn) + + # and normalize each point so that the sum equals unity + for i in range(0, len(self.__atoms)): + Pn[i,:] = np.divide(Pn[i,:], Pt) # cell function for each atom + + # store the result as a private variable + self.__mweights[g,:] = Pn[g,:] + + # and store the molecular weight functions into the atomic grids + atgrid.set_molecular_weights(self.__mweights[g]) + self.construct_times['fuzzy_cell_decomposition'] = time.time() - st + + # for each grid point in the molecular grid, store the distance to all the + # nuclei in the grid + st = time.time() + self.__rgridpoints = np.zeros((len(self.__atoms), np.prod(self.__mweights.shape))) + self.__xgridpoints = np.zeros_like(self.__rgridpoints) + self.__ygridpoints = np.zeros_like(self.__rgridpoints) + self.__zgridpoints = np.zeros_like(self.__rgridpoints) + gpts = np.array(np.vstack(self.__gridpoints)) # global grid points + for i,at in enumerate(self.__atoms): + ap = at[0] + # build relative coordinates + self.__xgridpoints[i,:] = np.subtract(gpts[:,0], ap[0]) + self.__ygridpoints[i,:] = np.subtract(gpts[:,1], ap[1]) + self.__zgridpoints[i,:] = np.subtract(gpts[:,2], ap[2]) + + # calculate radial distance between global gridpoint and atom i + xy = np.add(np.power(self.__xgridpoints[i,:], 2), + np.power(self.__ygridpoints[i,:], 2)) + xyz = np.add(xy, np.power(self.__zgridpoints[i,:], 2)) + self.__rgridpoints[i,:] = np.sqrt(xyz) + + # build unit sphere angles + self.__theta_gridpoints = np.arctan2(self.__ygridpoints, + self.__xgridpoints) + self.__phi_gridpoints = np.arccos(np.divide(self.__zgridpoints, + self.__rgridpoints)) + + # calculate the nuclear potential at each grid point due to the other nuclei + self.__rigridpoints = np.divide(1.0, self.__rgridpoints) + + self.__npot = np.einsum('i,ij->j', + [-at[1] for at in self.__atoms], + self.__rigridpoints) + self.construct_times['nuclear_distance_and_potential'] = time.time() - st + + # calculate values for the spherical harmonics for each grid point + # with respect to each atom as central point and for each value + # of l,m (thus a rank-3 tensor) + st = time.time() + self.__ylmgpts = np.ndarray((len(self.__atoms), + (self.__lmax+1)**2, + np.prod(self.__mweights.shape))) + for i,at in enumerate(self.__atoms): + lmctr = 0 + for l in range(0, self.__lmax+1): + for m in range(-l, l+1): + self.__ylmgpts[i,lmctr,:] = spherical_harmonic(l, m, \ + self.__theta_gridpoints[i,:], + self.__phi_gridpoints[i,:]) + lmctr += 1 + self.construct_times['spherical_harmonics'] = time.time() - st + + # collect complete weights for the full molecular grid + weights = np.array([an.get_weights() for an in self.__atomgrids]).flatten() + self.__mgw = np.multiply(weights, self.__mweights.flatten()) + + def __build_amplitudes(self): + """ + Precalculate the basis function amplitudes at the selected grid + points + """ + # calculate the amplitudes of the basis functions and store these + # per atomic grid, per basis function and per points in the atomic + # grid (i.e. in a rank-3 tensor) + st = time.time() + integrator = PyQInt() + + self.__amplitudes = np.ndarray((len(self.__atoms), len(self.__basis), len(self.__gridpoints[0]))) + for i,at in enumerate(self.__atoms): + for j,cgf in enumerate(self.__basis): + self.__amplitudes[i,j,:] = integrator.plot_wavefunction(np.array(self.__gridpoints[i]), [1], [cgf]) + + # also build amplitude gradients + self.__ampgrads = np.ndarray((len(self.__atoms), len(self.__basis), len(self.__gridpoints[0]), 3)) + for i,at in enumerate(self.__atoms): + for j,cgf in enumerate(self.__basis): + self.__ampgrads[i,j,:,:] = integrator.plot_gradient(np.array(self.__gridpoints[i]), [1], [cgf]) + + # reassemble the amplitudes per atomic grid into one for the complete + # molecular grid + self.__fullgrid_amplitudes = np.ndarray((len(self.__basis), np.prod(self.__mweights.shape))) + for i,cgf in enumerate(self.__basis): + self.__fullgrid_amplitudes[i,:] = np.hstack([self.__amplitudes[a,i,:] for a in range(len(self.__atoms))]) + self.construct_times['basis_function_amplitudes'] = time.time() - st + + def __step(self, mu): + """ + Becke fuzzy grid cut-off function + """ + if mu <= 0.0: + return 1.0 + else: + return 0.0 + + def __vij(self, mu, chi): + """ + Becke aij parameter for heteroatomic interactions + """ + uij = (chi - 1.0) / (chi + 1.0) + aij = uij / (uij**2 - 1.0) + aij = min(aij, 0.5) + aij = max(aij, -0.5) + return mu + aij * (1.0 - mu**2) + + def __sk(self, mu, n): + """ + Becke fuzzy grid cutoff profile + """ + return 0.5 * (1.0 - self.__fuzzy(mu,n)) + + def __fuzzy(self, mu, n): + """ + Becke fuzzy grid iterative function for adaptive sharpening of + the transition + """ + if n == 0: + return mu + else: + return self.__fuzzy(3.0/2.0 * mu - 0.5 * mu**3, n-1) \ No newline at end of file diff --git a/pydft/spherical_harmonics.py b/pydft/spherical_harmonics.py new file mode 100644 index 0000000..ca2f5dd --- /dev/null +++ b/pydft/spherical_harmonics.py @@ -0,0 +1,37 @@ +# -*- coding: utf-8 -*- + +from scipy.special import sph_harm +import numpy as np + +def spherical_harmonic(l, m, theta, phi): + """ + Calculate value of spherical harmonic function + + l: angular quantum number + m: magnetic quantum number + theta: azimuthal angle in radians + phi: polar angle in radians + """ + if m < 0: + val = np.sqrt(2) * np.imag(sph_harm(np.abs(m), l, theta, phi)) + elif m > 0: + val = np.sqrt(2) * np.real(sph_harm(m, l, theta, phi)) + else: + val = np.real(sph_harm(m, l, theta, phi)) + + return val + +def spherical_harmonic_cart(l, m, p): + """ + Calculate the value of the spherical harmonic depending on the position + p in Cartesian coordinates. This function assumes that the position p + lies on the unit sphere + + l: angular quantum number + m: magnetic quantum number + p: position three-vector on the unit sphere + """ + theta = np.arctan2(p[1], p[0]) + phi = np.arccos(p[2]) + + return spherical_harmonic(l, m, theta, phi) \ No newline at end of file diff --git a/pydft/xcfunctionals.py b/pydft/xcfunctionals.py new file mode 100644 index 0000000..cb65ba0 --- /dev/null +++ b/pydft/xcfunctionals.py @@ -0,0 +1,216 @@ +# -*- coding: utf-8 -*- + +import numpy as np + +class Functionals: + """ + Class holding all exchange-correlation functionals + + Note that all the functions ending in _deriv are the derivative + of the energy with respect to rho. + + To calculate the potential, nu, one has to calculate + deltaf/deltarho = df/drho * rho + f + """ + def __init__(self, functional = 'svwn5'): + functionals = { + 'svwn5': [self.__slater, self.__vwn5, + self.__slater_deriv, self.__vwn5_deriv, False], + 'pbe': [self.__pbe_x, self.__pbe_c, + self.__pbe_x_deriv, self.__pbe_c_deriv, True] + } + + if functional not in functionals.keys(): + raise Exception('Illegal XC-functional requested.') + else: + # store exchange and correlation functional + self.__xf = functionals[functional][0] + self.__cf = functionals[functional][1] + self.__dxf = functionals[functional][2] + self.__dcf = functionals[functional][3] + self.__gga = functionals[functional][4] + + def is_gga(self): + """ + Returns whether the loaded XC functional is of GGA type + """ + return self.__gga + + def calc_x(self, rho, grad=None): + """ + Calculate exchange using specified correlation potential + """ + with np.errstate(divide='ignore', invalid='ignore'): + rho = self.__parse_dens(rho) + if self.__gga: + grad = self.__parse_dens(grad) + ex = self.__xf(rho, grad) + fx = self.__dxf(rho, grad) + else: + ex = self.__xf(rho) + fx = self.__dxf(rho) + + return ex, (fx * rho + ex) + + def calc_c(self, rho, grad=None): + """ + Calculate correlation using specified correlation potential + """ + with np.errstate(divide='ignore', invalid='ignore'): + rho = self.__parse_dens(rho) + if self.__gga: + ex = self.__cf(rho, grad) + fx = self.__dcf(rho, grad) + else: + ex = self.__cf(rho) + fx = self.__dcf(rho) + + ex = np.nan_to_num(ex) + fx = np.nan_to_num(fx) + + return ex, (fx * rho + ex) + + def __parse_dens(self, dens): + return np.maximum(dens, 1e-12) + + def __slater(self, dens): + """ + Slater exchange functional + """ + return -(3/4) * (3 / np.pi)**(1/3) * np.power(dens, 1./3.) + + def __slater_deriv(self, dens): + """ + Slater exchange functional derivative towards rho + """ + return -(3 / np.pi)**(1/3) / 4 / dens**(2/3) + + def __vwn5(self, dens): + # Note that equation E.27 is incorrect in the book of Parr and Yang + # but is correct in the original paper of Vosko, Wilk and Nusair + # https://cdnsciencepub.com/doi/pdf/10.1139/p80-159 + A = 0.0621814 + x0 = -0.409286 + b = 13.0720 + c = 42.7198 + + rs = (3. / 4. / np.pi / dens)**(1./3.) + + x = rs**(1/2) + X = x**2 + b * x + c + X0 = x0**2 + b * x0 + c + Q = (4 * c - b**2)**(1/2) + atan = np.arctan(Q / (2*x+b)) + + return A/2 * (np.log(x**2/X) + 2*b/Q * atan - b*x0/X0 * (np.log((x-x0)**2 / X) + 2 * (b + 2*x0) / Q * atan)) + + def __vwn5_deriv(self, dens): + """ + Derivative of the VWN5 correlation functional towards the density + """ + A = 0.0621814 + x0 = -0.409286 + b = 13.0720 + c = 42.7198 + X0 = x0**2 + b * x0 + c + Q = (4 * c - b**2)**(1/2) + + rs = (3. / 4. / np.pi / dens)**(1./3.) + x = rs**(1/2) + + T1 = (2*c + b*x)/(x*(c + x*(b + x))) + T2 = (4*b)/(Q**2 + (b + 2*x)**2) + T3 = (b*x0*(-((b + 2*x)/(c + x*(b + x))) + 2/(x - x0) - (4*(b + 2*x0))/(Q**2 + (b + 2*x)**2))) / X0 + + dfdx = A/2 * (T1 - T2 - T3) + dxdrho = -((1/dens)**(7/6)/(2*2**(1/3)*3**(5/6)*np.pi**(1/6))) + return dfdx * dxdrho + + def __pbe_x(self, rho, gamma): + """ + PBE exchange functional for spin-unpolarized density + """ + R = 0.804 + mu = 0.2195149727645171 + + c1 = 0.738558766382 + c2 = 0.0192920212964 + c3 = 0.0261211729852 + + return rho**(1/3) * (-c1 - c2 * gamma * mu * R \ + / (c3 * gamma * mu + R * rho**(8/3))) + + def __pbe_x_deriv(self, rho, gamma): + """ + PBE exchange derivative for spin-unpolarized density using reparametrization + """ + R = 0.804 + mu = 0.2195149727645171 + + c1 = 0.246186255461 + c2 = 18.8495559215 + c3 = 65.9734457254 + c4 = 360.809905669 + c5 = 38.2831200025 + + t1 = gamma**2 * mu**2 * (-c1 - c1 * R) + t2 = gamma * mu * R * (-c2 + c3 * R) * rho**(8/3) + t3 = -c4 * R**2 * rho**(16/3) + t4 = (rho**(2/3) * (gamma * mu + c5 * R * rho**(8/3))**2) + + return (t1 + t2 + t3) / t4 + + def __pbe_x_deriv_numerical(self, rho, gamma): + """ + PBE Exchange derivative using numerical approximation + """ + # use finite difference discretization to approximate result + dx = 1e-5 + return (self.pbe_x_deriv_simplified(rho + dx, gamma) - self.pbe_x_deriv_simplified(rho - dx, gamma)) / (2. * dx) + + def __pbe_c(self, rho, gamma): + """ + Reparametrization of the PBE correlation functional for spin-unpolarized density + """ + c1 = -0.0621814 + c2 = 0.008243319792565314 + c3 = 0.3720033616668558 + c4 = 0.13838902240433937 + c5 = 0.049771773124349425 + c6 = 0.011795838478103926 + c7 = 0.031090690869654897 + c8 = 7.341562356353513 + c9 = 0.13621078885675922 + c10 = 0.3720033616668558 + c11 = 0.049771773124349425 + c12 = 2.0000005873362636 + c13 = 0.2651378776729259 + + irho = 1/rho + irho12 = np.sqrt(irho) + irho16 = (irho)**(1/6) + irho13 = irho16**2 + irho23 = irho13**2 + + c3irho16 = c3*irho16 + c6irho23 = c6*irho23 + c13irho13 = c13*irho13 + c4irho13 = c4*irho13 + c5irho12 = c5*irho12 + c9gamma = c9 * gamma + + t1 = (c1 - c2*irho13) + t2 = 1 + 1/(c3irho16 + c4irho13 + c5irho12 + c6irho23) + t3 = (-1. + 1.*(1 + 1/(c3irho16 + c4irho13 + c5irho12 + c6irho23))**(c12 + c13irho13)) + t4 = (-1. + 1.*(1 + 1/(c10*irho16 + c4irho13 + c11*irho12 + c6irho23))**(c12 + c13irho13)) + t5 = 1. + 1 / ((c8*rho**(7/3))/gamma + (c9gamma)/(t4*(c9gamma + t3*rho**(7/3)))) + + return t1*np.log(t2) + c7*np.log(t5) + + def __pbe_c_deriv(self, rho, gamma): + """ + PBE correlation derivative using numerical approximation + """ + # use finite difference discretization to approximate result + dx = 1e-5 + return (self.__pbe_c(rho + dx, gamma) - self.__pbe_c(rho - dx, gamma)) / (2. * dx) \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..91cb94f --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,46 @@ +[build-system] +requires = ["setuptools", "setuptools-scm"] +build-backend = "setuptools.build_meta" + +[project] +name = "pydft" +version = "0.6.1" +authors = [ + { name="Ivo Filot", email="i.a.w.filot@tue.nl" } +] +maintainers = [ + { name="Ivo Filot", email="i.a.w.filot@tue.nl" }, +] +description = "pydft" +readme = "README.md" +requires-python = ">=3.9" +license = {text = "GPL-3.0-or-later"} +classifiers = [ + "Programming Language :: Python :: 3", + "Operating System :: OS Independent", +] +dependencies = [ + "scipy", + "numpy", + "pyqint", + "pylebedev", + "mendeleev", +] + +[tool.setuptools] +include-package-data = true + +[tool.setuptools.packages.find] +where = ["."] +include = ["pydft*"] + +[tool.conda.environment] +name = "demo" +channels = [ + "ifilot", + "conda-forge", +] + +[project.urls] +"Homepage" = "https://gitlab.tue.nl/ifilot/pydft" +"Bug Tracker" = "https://gitlab.tue.nl/ifilot/pydft/-/issues" \ No newline at end of file diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..628db6e --- /dev/null +++ b/setup.py @@ -0,0 +1,50 @@ +from setuptools import Extension, setup +import os +import sys +import re + +PKG = "pydft" +VERSIONFILE = os.path.join(os.path.dirname(__file__), PKG, "_version.py") +verstr = "unknown" +try: + verstrline = open(VERSIONFILE, "rt").read() +except EnvironmentError: + pass # Okay, there is no version file. +else: + VSRE = r"^__version__ = ['\"]([^'\"]*)['\"]" + mo = re.search(VSRE, verstrline, re.M) + if mo: + verstr = mo.group(1) + else: + print(r"Unable to find version in %s" % (VERSIONFILE,)) + raise RuntimeError(r"If %s.py exists, it is required to be well-formed" % (VERSIONFILE,)) + +with open("README.md", "r", encoding="utf-8") as fh: + long_description = fh.read() + +def package_files(directory): + paths = [] + for (path, directories, filenames) in os.walk(directory): + for filename in filenames: + paths.append(os.path.join('..', path, filename)) + return paths + +setup( + name=PKG, + version=verstr, + author="Ivo Filot", + author_email="ivo@ivofilot.nl", + description="Python package for performing simple DFT calculations", + long_description=long_description, + long_description_content_type="text/markdown", + url="https://gitlab.tue.nl/ifilot/pydft", + packages=[PKG], + include_package_data=True, + classifiers=[ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: MIT License", + "Operating System :: POSIX", + ], + python_requires='>=3.5', + zip_safe=False +) diff --git a/tests/test_bragg_slater.py b/tests/test_bragg_slater.py new file mode 100644 index 0000000..d58119b --- /dev/null +++ b/tests/test_bragg_slater.py @@ -0,0 +1,56 @@ +import unittest +import numpy as np +from pydft import MoleculeBuilder, MolecularGrid + +class TestCustomBasisSet(unittest.TestCase): + + def test_bragg_slater_bf3(self): + """ + Test grabbing correct Bragg-Slater radii + """ + # construct molecule + mol = MoleculeBuilder().from_name('bf3') + cgfs, atoms = mol.build_basis('sto3g') + + # construct molecular grid + molgrid = MolecularGrid(atoms, cgfs) + molgrid.initialize() + + ANGSTROM2BOHR = 1.88973 + + # test Bragg-Slater radius for B + np.testing.assert_almost_equal( + molgrid.get_atomic_grid(0).get_bragg_slater_radius(), + 0.85 * ANGSTROM2BOHR) + + # test Bragg-Slater radius for F + np.testing.assert_almost_equal( + molgrid.get_atomic_grid(1).get_bragg_slater_radius(), + 0.50 * ANGSTROM2BOHR) + + def test_bragg_slater_ch4(self): + """ + Test grabbing correct Bragg-Slater radii + """ + # construct molecule + mol = MoleculeBuilder().from_name('ch4') + cgfs, atoms = mol.build_basis('sto3g') + + # construct molecular grid + molgrid = MolecularGrid(atoms, cgfs) + molgrid.initialize() + + ANGSTROM2BOHR = 1.88973 + + # test Bragg-Slater radius for B + np.testing.assert_almost_equal( + molgrid.get_atomic_grid(0).get_bragg_slater_radius(), + 0.70 * ANGSTROM2BOHR) + + # test Bragg-Slater radius for F + np.testing.assert_almost_equal( + molgrid.get_atomic_grid(1).get_bragg_slater_radius(), + 0.35 * ANGSTROM2BOHR) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/test_custom_basisset.py b/tests/test_custom_basisset.py new file mode 100644 index 0000000..3d639e6 --- /dev/null +++ b/tests/test_custom_basisset.py @@ -0,0 +1,33 @@ +import unittest +import numpy as np +from pydft import MoleculeBuilder, DFT +from pyqint import cgf + +class TestCustomBasisSet(unittest.TestCase): + + def test_h2(self): + """ + Test DFT calculation of hydrogen molecule + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('h2') + + cgfs = [] + for n in mol.get_nuclei(): + _cgf = cgf(n[0]) + + _cgf.add_gto(0.154329, 3.425251, 0, 0, 0) + _cgf.add_gto(0.535328, 0.623914, 0, 0, 0) + _cgf.add_gto(0.444635, 0.168855, 0, 0, 0) + + cgfs.append(_cgf) + + # construct dft object + dft = DFT(mol, basis=cgfs) + energy = dft.scf() + + answer = -1.1570136 + np.testing.assert_almost_equal(energy, answer, 4) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_dft.py b/tests/test_dft.py new file mode 100644 index 0000000..8e3f661 --- /dev/null +++ b/tests/test_dft.py @@ -0,0 +1,148 @@ +import unittest +import numpy as np +from pydft import MoleculeBuilder, DFT + +class TestDFT(unittest.TestCase): + + def test_helium(self): + """ + Test DFT calculation of Helium atom + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('He') + + # construct dft object + dft = DFT(mol, basis='sto3g') + energy = dft.scf() + + answer = -2.809567 + np.testing.assert_almost_equal(energy, answer, 4) + + def test_h2o(self): + """ + Test DFT calculation of water molecule + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('H2O') + + # construct dft object + dft = DFT(mol, basis='sto3g') + energy = dft.scf() + + answer = -74.9393412503195 + np.testing.assert_almost_equal(energy, answer, 4) + + def test_co(self): + """ + Test DFT calculation of CO molecule + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('CO') + + # construct dft object + dft = DFT(mol, basis='sto3g') + energy = dft.scf() + + answer = -111.1458218142885 + np.testing.assert_almost_equal(energy, answer, 4) + + def test_bf3(self): + """ + Test DFT calculation of BF3 molecule + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('bf3') + + # construct dft object + dft = DFT(mol, basis='sto3g') + energy = dft.scf() + + answer = -318.2583234959651 + np.testing.assert_almost_equal(energy, answer, 4) + + def test_ch4(self): + """ + Test DFT calculation of methane molecule + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('ch4') + + # construct dft object + dft = DFT(mol, basis='sto3g') + energy = dft.scf() + + answer = -39.801509613192096 + np.testing.assert_almost_equal(energy, answer, 4) + + def test_co2(self): + """ + Test DFT calculation of co2 molecule + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('co2') + + # construct dft object + dft = DFT(mol, basis='sto3g') + energy = dft.scf() + + answer = -185.0106604454506 + np.testing.assert_almost_equal(energy, answer, 4) + + def test_ethylene(self): + """ + Test DFT calculation of an ethylene molecule + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('ethylene') + + # construct dft object + dft = DFT(mol, basis='sto3g') + energy = dft.scf() + + answer = -77.16054913818746 + np.testing.assert_almost_equal(energy, answer, 4) + + def test_h2(self): + """ + Test DFT calculation of hydrogen molecule + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('h2') + + # construct dft object + dft = DFT(mol, basis='sto3g') + energy = dft.scf() + + answer = -1.1570136 + np.testing.assert_almost_equal(energy, answer, 4) + + def test_lih(self): + """ + Test DFT calculation of LiH molecule + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('lih') + + # construct dft object + dft = DFT(mol, basis='sto3g') + energy = dft.scf() + + answer = -7.86536764070593 + np.testing.assert_almost_equal(energy, answer, 4) + + # def test_benzene(self): + # """ + # Test DFT calculation of benzene molecule + # """ + # mol_builder = MoleculeBuilder() + # mol = mol_builder.from_name('benzene') + + # # construct dft object + # dft = DFT(mol, basis='sto3g') + # energy = dft.scf() + + # answer = -228.0695756259556 + # np.testing.assert_almost_equal(energy, answer, 4) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_energy_decomposition.py b/tests/test_energy_decomposition.py new file mode 100644 index 0000000..b914a5b --- /dev/null +++ b/tests/test_energy_decomposition.py @@ -0,0 +1,72 @@ +# -*- coding: utf-8 -*- + +import numpy as np +import unittest +import sys, os + +# add a reference to load the PyPWDFT module +sys.path.append(os.path.join(os.path.dirname(__file__), '..')) + +from pydft import DFT, MoleculeBuilder +import pydft + +class TestEnergyDecomposition(unittest.TestCase): + """ + Test existence of dictionary keys + """ + + def test_energy_decomposition(self): + """ + Test energy decomposition of a molecule + """ + co = MoleculeBuilder().from_name("CO") + dft = DFT(co, basis='sto3g') + en = dft.scf(1e-4) + print("Total electronic energy: %f Ht" % en) + + # retrieve molecular matrices + res = dft.get_data() + P = res['P'] + T = res['T'] + V = res['V'] + J = res['J'] + + # calculate energy terms + Et = np.einsum('ji,ij', P, T) + Ev = np.einsum('ji,ij', P, V) + Ej = 0.5 * np.einsum('ji,ij', P, J) + Ex = res['Ex'] + Ec = res['Ec'] + Exc = res['Exc'] + Enuc = res['enucrep'] + + # print('Kinetic energy: %12.6f' % Et) + # print('Nuclear attraction: %12.6f' % Ev) + # print('Electron-electron repulsion: %12.6f' % Ej) + # print('Exchange energy: %12.6f' % (Ex)) + # print('Correlation energy: %12.6f' % (Ec)) + # print('Exchange-correlation energy: %12.6f' % (Exc)) + # print('Nucleus-nucleus repulsion: %12.6f' % (Enuc)) + + esum = Et + Ev + Ej + Exc + Enuc + np.testing.assert_almost_equal(esum, en, decimal=5) + + def test_key_existence(self): + """ + Test whether all keys are in dictionary + """ + h2 = MoleculeBuilder().from_name("H2") + dft = DFT(h2, basis='sto3g') + en = dft.scf(1e-4) + mats = dft.get_data() + + keys = ['S', 'T', 'V', 'C', 'J', 'P', 'XC', 'F', 'Exc', 'Ex', 'Ec', 'energies', 'energy', 'orbc', 'orbe', 'enucrep'] + + for key in keys: + assert key in mats + + # also test for false positives + assert 'fakekey' not in mats + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/test_moleculargrid.py b/tests/test_moleculargrid.py new file mode 100644 index 0000000..0f5081f --- /dev/null +++ b/tests/test_moleculargrid.py @@ -0,0 +1,53 @@ +# -*- coding: utf-8 -*- + +from pydft import MoleculeBuilder, MolecularGrid +import numpy as np +import unittest + +class TestMolecularGrid(unittest.TestCase): + """ + Perform a quick test to verify whether all DFT routines + are properly working + """ + + def test_weight_functions(self): + """ + Test DFT calculation of Helium atom + """ + # construct molecule + mol = MoleculeBuilder().from_name('bf3') + cgfs, atoms = mol.build_basis('sto3g') + + # construct molecular grid + molgrid = MolecularGrid(atoms, cgfs) + + # produce grid of sampling points to calculate the atomic + # weight coefficients for + N = 10 + sz = 8 + x = np.linspace(-sz,sz,N) + xv,yv = np.meshgrid(x,x) + points = np.array([[x,y,0] for x,y in zip(xv.flatten(),yv.flatten())]) + + # calculate the atomic weights + mweights = molgrid.calculate_weights_at_points(points, k=3) + + # verify results + np.testing.assert_almost_equal(mweights[0,40], + 7.438664028803838e-06, 4) + np.testing.assert_almost_equal(mweights[0,50], + 0.0023096748652996083, 4) + np.testing.assert_almost_equal(mweights[0,60], + 0.03962756655918886, 4) + np.testing.assert_almost_equal(mweights[0,70], + 0.08886587061825403, 4) + np.testing.assert_almost_equal(mweights[0,80], + 0.04752230028517094, 4) + + np.testing.assert_almost_equal(mweights[1,70], + 0.41031416799944265, 4) + np.testing.assert_almost_equal(mweights[1,80], + 0.7991044370004883, 4) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/test_version.py b/tests/test_version.py new file mode 100644 index 0000000..a3a05be --- /dev/null +++ b/tests/test_version.py @@ -0,0 +1,13 @@ +import unittest +import pydft + +class TestVersion(unittest.TestCase): + + def test_version(self): + strver = pydft.__version__.split('.') + self.assertTrue(strver[0].isnumeric()) + self.assertTrue(strver[1].isnumeric()) + self.assertTrue(strver[2].isnumeric()) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/test_xc.py b/tests/test_xc.py new file mode 100644 index 0000000..93e54dc --- /dev/null +++ b/tests/test_xc.py @@ -0,0 +1,33 @@ +import unittest +from pydft import MoleculeBuilder, DFT +import numpy as np + +class TestXC(unittest.TestCase): + + def test_svwn5(self): + """ + Test DFT calculation of Helium atom + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('He') + + # construct dft object + dft = DFT(mol, basis='sto3g', functional='svwn5') + energy = dft.scf() + + answer = -2.809567 + np.testing.assert_almost_equal(energy, answer, 4) + + def test_pbe(self): + """ + Test DFT calculation of Helium atom + """ + mol_builder = MoleculeBuilder() + mol = mol_builder.from_name('He') + + # construct dft object + dft = DFT(mol, basis='sto3g', functional='pbe') + energy = dft.scf() + + answer = -2.8301366521757787 + np.testing.assert_almost_equal(energy, answer, 4) \ No newline at end of file diff --git a/versiontest.py b/versiontest.py new file mode 100644 index 0000000..acc3cde --- /dev/null +++ b/versiontest.py @@ -0,0 +1,70 @@ +# -*- coding: utf-8 -*- + +import os +import re + +ROOT = os.path.dirname(__file__) + +def main(): + version_projecttoml = get_version_projecttoml() + version_versionpy = get_version_versionpy() + version_metayaml = get_version_metayaml() + + try: + for i in range(0,3): + assert version_projecttoml[i] == version_versionpy[i] + assert version_versionpy[i] == version_metayaml[i] + except Exception as e: + print(e) + raise Exception('Invalid version strings encountered') + +def get_version_projecttoml(): + """ + Extract the version string from the pyproject.toml file + """ + pattern = re.compile(r'^version\s*=\s*"(\d+\.\d+.\d+)"\s*$') + + f = open(os.path.join(ROOT, 'pyproject.toml')) + lines = f.readlines() + for line in lines: + match = re.match(pattern, line) + if match: + version = match.groups(1)[0] + return [int(i) for i in version.split('.')] + + return None + +def get_version_versionpy(): + """ + Extract the version string from the _version.py file + """ + pattern = re.compile(r'^__version__\s*=\s*\'(\d+\.\d+.\d+)\'\s*$') + + f = open(os.path.join(ROOT, 'pydft', '_version.py')) + lines = f.readlines() + for line in lines: + match = re.match(pattern, line) + if match: + version = match.groups(1)[0] + return [int(i) for i in version.split('.')] + + return None + +def get_version_metayaml(): + """ + Extract the version string from the pyproject.toml file + """ + pattern = re.compile(r'^\s*version\s*:\s*"(\d+\.\d+.\d+)"\s*$') + + f = open(os.path.join(ROOT, 'meta.yaml')) + lines = f.readlines() + for line in lines: + match = re.match(pattern, line) + if match: + version = match.groups(1)[0] + return [int(i) for i in version.split('.')] + + return None + +if __name__ == '__main__': + main()