Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Re-enable cp2k example #44

Merged
merged 2 commits into from
Mar 26, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/docs.yml
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ jobs:
- roy-gch
- sample-selection
- gaas-map
# - batch-cp2k
- batch-cp2k

steps:
- uses: actions/checkout@v4
Expand Down
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
production/
parameters/

cp2k.out
cp2k.inp
cp2k_shell.ssmp
File renamed without changes.
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
channels:
- conda-forge
dependencies:
- python=3.11
- pip
Expand Down
File renamed without changes.
Original file line number Diff line number Diff line change
Expand Up @@ -8,34 +8,46 @@
using `CP2K <https://www.cp2k.org>`_ using its `reftraj functionality
<https://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/MD/REFTRAJ.html>`_. The inputs are a
set of structures in :download:`example.xyz` using the DFT parameters defined in
:download:`reftraj_template.cp2k` importing basis set and pseudopotentials from the
local CP2K installation. The reference DFT parameters are taken from `Cheng et al. Ab
initio thermodynamics of liquid and solid water 2019
:download:`reftraj_template.cp2k`. The reference DFT parameters are taken from `Cheng et
al. Ab initio thermodynamics of liquid and solid water 2019
<https://www.pnas.org/doi/10.1073/pnas.1815117116>`_. Due to the small size of the test
structure and convergence issues, we have decreased the size of the ``CUTOFF_RADIUS``
from :math:`6.0\,\mathrm{Å}` to :math:`3.0\,\mathrm{Å}`. For actual production
calculations adapt the template!

To run this example, we use a bare executable called with ``cp2k``. If you want to use
another version you can either adjust the the names within this example or link your
binary with a different name to ``cp2k``.
"""

# %%
# We start the example by importing the required packages.


import os
import re
import platform
import subprocess
from os.path import basename, splitext
from typing import List, Union

import ase.io
import ase.visualize.plot
import matplotlib.pyplot as plt
import numpy as np
from ase.calculators.cp2k import CP2K
import requests


# %%
#
# Install CP2K
# ------------
#
# We'll need a working installation of cp2k. The best way to do so depends on your
# platform, here are some possible solutions, but feel free to replace them with another
# installation method.

if platform.system() == "Linux":
# use conda on Linux
subprocess.run(["conda", "install", "cp2k", "-c", "conda-forge", "-y"], check=True)
elif platform.system() == "Darwin":
# use homebrew on macOS
subprocess.run(["brew", "install", "cp2k"], check=True)
else:
print("no known way to install cp2k, skipping installation")


# %%
Expand Down Expand Up @@ -136,17 +148,25 @@ def write_cp2k_in(


# %%
#
# We will now download basis set files from CP2K website. Depending on your CP2K
# installation, this might not be necessary!


def mkdir_force(*args, **kwargs) -> None:
"""Warpper to ``os.mkdir``.
def download_parameter(file):
path = os.path.join("parameters", file)

if not os.path.exists(path):
url = f"https://raw.githubusercontent.com/cp2k/cp2k/support/v2024.1/data/{file}"
response = requests.get(url)
response.raise_for_status()
with open(path, "wb") as f:
f.write(response.content)

The function does not raise an error if the directory already exists.
"""
try:
os.mkdir(*args, **kwargs)
except OSError:
pass

os.makedirs("parameters", exist_ok=True)
for file in ["GTH_BASIS_SETS", "BASIS_ADMM", "POTENTIAL", "dftd3.dat", "t_c_g.dat"]:
download_parameter(file)


# %%
Expand Down Expand Up @@ -200,11 +220,9 @@ def mkdir_force(*args, **kwargs) -> None:
# directory named ``H4O2`` because our dataset consists only of a single structure with
# two water molecules.

mkdir_force(project_directory)

for stoichiometry, frames in frames_dict.items():
current_directory = f"{project_directory}/{stoichiometry}"
mkdir_force(current_directory)
os.makedirs(current_directory, exist_ok=True)

write_cp2k_in(
f"{current_directory}/in.cp2k",
Expand Down Expand Up @@ -281,50 +299,5 @@ def mkdir_force(*args, **kwargs) -> None:

new_frames += frames_dft

new_fname = f"{splitext(basename(write_to_file))[0]}_dft.xyz"
new_fname = f"{os.path.splitext(os.path.basename(write_to_file))[0]}_dft.xyz"
ase.io.write(f"{project_directory}/{new_fname}", new_frames)

# %%
# Perform calculations using ASE calculator
# -----------------------------------------
# Above we performed the calculations using an external bash script. ASE also provides a
# calculator class that we can use the perform the calculations with our input file
# without a detour of writing files to disk.
#
# To use the ASE calculator together with a custom input script this requires some
# adjustments. First the name of the executable that has the exact name ``cp2k_shell``.
# We create a symlink to follow this requirement.

# %%
# Next, we load the input file abd remove ``GLOBAL`` section because from it

inp = open("./production/H4O2/in.cp2k", "r").read()
inp = re.sub(
f"{re.escape('&GLOBAL')}.*?{re.escape('&END GLOBAL')}", "", inp, flags=re.DOTALL
)

# %%
# Afterwards we define the :py:class:`ase.calculators.cp2k.CP2K`` calculator. Note that
# we disable all parameters because we want to use all options from our input file

calc = CP2K(
inp=inp,
max_scf=None,
cutoff=None,
xc=None,
force_eval_method=None,
basis_set=None,
pseudo_potential=None,
basis_set_file=None,
potential_file=None,
stress_tensor=False,
poisson_solver=None,
print_level=None,
)

# %%
# We now load a new structure, add the calculator and perform the computation.

atoms = ase.io.read("example.xyz")
atoms.set_calculator(calc)
# atoms.get_potential_energy()
Original file line number Diff line number Diff line change
Expand Up @@ -52,9 +52,9 @@
&END FORCES
&END PRINT
&DFT
BASIS_SET_FILE_NAME GTH_BASIS_SETS
BASIS_SET_FILE_NAME BASIS_ADMM
POTENTIAL_FILE_NAME POTENTIAL
BASIS_SET_FILE_NAME ../../parameters/GTH_BASIS_SETS
BASIS_SET_FILE_NAME ../../parameters/BASIS_ADMM
POTENTIAL_FILE_NAME ../../parameters/POTENTIAL
&MGRID
CUTOFF 400
&END MGRID
Expand Down Expand Up @@ -97,7 +97,7 @@
&INTERACTION_POTENTIAL
POTENTIAL_TYPE TRUNCATED
CUTOFF_RADIUS 3.0
T_C_G_DATA t_c_g.dat
T_C_G_DATA ../../parameters/t_c_g.dat
&END
&HF_INFO
&END HF_INFO
Expand All @@ -109,7 +109,7 @@
R_CUTOFF 15
LONG_RANGE_CORRECTION TRUE
REFERENCE_FUNCTIONAL revPBE0
PARAMETER_FILE_NAME dftd3.dat
PARAMETER_FILE_NAME ../../parameters/dftd3.dat
&END
&END
&XC_GRID
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,6 @@

for i in $(find ./production/ -mindepth 1 -type d); do
cd $i
cp2k -i in.cp2k
cp2k.ssmp -i in.cp2k
cd -
done