Calculations with cubic-scaling BigDFT code
This material is adapted from work by Luigi Genovese, Laura Ratcliff and others
In this tutorial we will show some basic calculations of some molecular systems performed with the traditional direct-minization algorithm of BigDFT code. We will base our examples on a system of water molecules.
Let’s first start by copying the xyz file associated to the water molecules. We will use for this the XYZReader class:
!curl https://gitlab.com/max-centre/benchmarks/-/raw/master/BigDFT/H2O/CPU/posinp.xyz -o H2O-32.xyz 2> /dev/null
Be sure to set the user!
user = 'eugenovesel'
#if you are in a Google colab session
try:
from google.colab import drive
in_colab = True
except:
in_colab = False
print('We are in a google colab session: ',in_colab)
if in_colab:
import subprocess
print(subprocess.check_output(['pip', 'install', '-U','pybigdft','pyfutile','py3Dmol','remotemanager']))
We are in a google colab session: False
# this sets up some internal variables of the remotemanager plugin
import remotemanager
remotemanager.Logger.level = 'debug'
remotemanager.Logger.path = 'H2O-vega'
After some session setup, let’s start by defining the System:
from BigDFT.IO import XYZReader
from BigDFT.Systems import System
from BigDFT.Fragments import Fragment
sys = System()
with XYZReader("H2O-32.xyz") as ifile:
for iat, atom in enumerate(ifile):
if iat % 3 == 0:
frag = Fragment()
sys['HOH:'+str(int(iat/3)+1)] = frag
frag.append(atom)
sys.cell=ifile.cell
sys.display()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
<BigDFT.Visualization.InlineVisualizer at 0x7f1d8c306730>
Single water molecule calculation
Let us now start by a simple example of one single water molecule. We take as reference a system made by one of the molecues of the total assembly
choice='HOH:16'
H2O = System({choice:sys[choice]})
H2O.display()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
<BigDFT.Visualization.InlineVisualizer at 0x7f1d68213bb0>
We save this system in a xyz file for future reuse.
from BigDFT.IO import write_xyz
with open('H2O-1.xyz','w') as ofile:
write_xyz(system=H2O,ofile=ofile)
Remote calculation of a dataset of different runs of the water molecule
We would like to employ the remote machine to run different calculations of the water molecule. We will perform single point calculations with:
LDA functional, in gas phase PBE functional, in gas phase PBE functional, in implicit solvent PBE0 functional, in gas phase (with CPU only) Let us now construct the remote function we would like to use for the calculations. Such function employs some of the functionalities of the BigDFT code and the role of each line is documented. Of course this serves only as an overview of the BigDFT functionalities. We encourage to have a look at the other tutorials of the code to inspect more the way in which such functionalities can be tuned.
def calculate_single_point(run_name, filename, functional, solvent=False, run_directory='Calculations', gpu=False):
from BigDFT import Calculators as C, Inputfiles as I
from futile.Utils import create_tarball
from os.path import join
# create a reasonable set of input parameters for the molecule
inp = I.Inputfile()
# self-explanatory line
inp.set_xc(functional)
# grid spacing (like a plane-wave cutoff, the lower the more precise)
inp.set_hgrid(0.4)
# Get the pseudpotential files
inp.set_psp_nlcc()
# include soft-sphere implicit model if required
if solvent:
inp.set_implicit_solvent(solvent='water',minres=1.0e-8)
if gpu:
inp.setdefault('psolver',{}).setdefault('setup',{}).update({'accel': 'CUDA'})
# create the calculator
code = C.SystemCalculator(skip=True)
# launch the calculation and retrieve the logfile instance
log = code.run(input=inp, name=run_name, posinp=filename, run_dir=run_directory)
# we create a tarfile that contains the logfile as well as the information about the timing
create_tarball(run_name+'.tar.gz',files=[join(run_directory,'log-'+run_name+'.yaml'),
join(run_directory,'time-'+run_name+'.yaml')])
# return the energy
return log.energy
Such function can be executed locally (we encourage you to do so in a workstation where bigdft-suite is installed). Here, we create a dataset of such remote functions to execute the pool of jobs into Vega, from this notebook.
from remotemanager import Dataset
if in_colab:
from remotemanager import URL
url = URL()
else:
from remotemanager.connection.computers import Vega
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In[10], line 6
4 url = URL()
5 else:
----> 6 from remotemanager.connection.computers import Vega
ImportError: cannot import name 'Vega' from 'remotemanager.connection.computers' (/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/remotemanager/connection/computers/__init__.py)
As per the remotemanager tutorial, we set up the information about the vega setup. It is important also to set the information about the gpu script. Such information can be provided in a configuration file that is associate to the computer.
import yaml
# vega yaml specification:
vega_spec=yaml.load("""
modules:
- SciPy-bundle/2022.05-intel-2022a
- CMake/3.23.1-GCCcore-11.3.0
- Meson/0.62.1-GCCcore-11.3.0
modules_gpu:
- SciPy-bundle/2020.03-foss-2020a-Python-3.8.2
- fosscuda
- MKL
- CMake/3.18.4-GCCcore-10.2.0
- Meson/0.55.3-GCCcore-10.2.0
- GSL/2.6-GCC-10.2.0
- Doxygen/1.8.20-GCCcore-10.2.0
extra: |
export PREFIX=/ceph/hpc/data/d2021-135-users/softwares/BigDFT/binaries/install
source $PREFIX/bin/bigdftvars.sh
export OMPI_MCA_orte_base_help_aggregate=0
export BIGDFT_MPIRUN=mpirun
export FUTILE_PROFILING_DEPTH=0
export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK
extra_gpu: |
export PREFIX=/ceph/hpc/data/d2021-135-users/softwares/BigDFT/binaries-cuda/install
source $PREFIX/bin/bigdftvars.sh
export OMPI_MCA_orte_base_help_aggregate=0
export BIGDFT_MPIRUN=mpirun
export FUTILE_PROFILING_DEPTH=0
export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK
""", Loader=yaml.Loader)
We save the specification file on the disk for reuse in the following tutorial
with open('vega.yaml','w') as ofile:
yaml.dump(vega_spec,ofile)
if not in_colab:
url = Vega(host='vega')
for attribute, value in vega_spec.items():
setattr(url, attribute, value)
!wget https://raw.githubusercontent.com/BigDFT-group/bigdft-school/main/packaging/install.py 2> /dev/null
import install
if in_colab:
install.data('data/H2O-CS.tar.xz')
from os.path import join, relpath
from BigDFT.RemoteRunners import computer_runner, RemoteDataset
local_dir = 'H2O-single'
remote_dir = f'/ceph/hpc/data/d2021-135-users/bigdft-school/{user}/H2O-single' #use your own directory for new calculations
force_reinitialization = not in_colab # in case you rerun this cell multiple times
filename = 'H2O-1.xyz'
all_runs = Dataset(function = calculate_single_point,
url = url,
name = 'H2O-single',
remote_dir = remote_dir,
local_dir = local_dir,
extra_files_send=filename,
block_reinit = force_reinitialization)
#we collect the run names for future reuse
run_names = []
# the we append the runs we would like to perform
for functional, solvent in [('LDA', False), ('PBE', False), ('PBE', True), ('PBE0', False)]:
#identify a name for the run
run_name = '-'.join([filename.split('.')[0],functional] + (['PCM'] if solvent else []))
# we then have the arguments of the function to submit
func_kwargs = dict(run_name=run_name, filename=filename, functional=functional, solvent=solvent)
# we decide to execute such calculations with cpu, 4 MPI and 4 OMP threads
run_args = {}
run_args['jobname'] = run_name
run_args['mpi'] = 4
run_args['omp'] = 4
run_args['time'] = 2 # minutes
all_runs.append_run(args=func_kwargs,
extra_files_recv=[run_name+'.tar.gz'],
**run_args)
run_names.append(run_name)
#for slow connections
# all_runs.url.timeout=120
# all_runs.url.max_timeouts=20
# all_runs.url.raise_errors=False #to override ssh error messages
all_runs.run()
all_runs.run_cmds
The jobs are now submitted on vega. In case you would like to verify, log in on the machine and verify the squeue command.
#you can also run this:
all_runs.url.cmd('squeue --me')
print(all_runs.is_finished)
all_runs.fetch_results(ignore_errors=True)
energies = all_runs.results
energies
Analysis of the runs
After those runs have been executed we can compare different quantities among the various runs, for illustrative purposes. We will highlight:
The Density of States
The Molecule dipole
The profiling of the different timing categories
We start by unpacking the received data
def extract_results(directory):
from futile.Utils import file_list
import tarfile
from os.path import join
from os import system
for archive in file_list(directory,suffix='.tar.gz',exclude='files',include_directory_path=True):
arch = tarfile.open(archive)
arch.extractall(path=directory)
arch.close()
#patch the yaml files from a mistake in the case of low profiling depth
system('sed -i s/^\ *\:\ null/\ \ \ \ null/g '+join(directory,'Calculations','time*'))
system('sed -i s/^\ *\:\ null/\ \ \ \ null/g '+join(directory,'Calculations','data*','time*'))
# the data are unpacked in a "Calculations" subdirectory (see above arguments)
extract_results('H2O-single')
from BigDFT.Logfiles import Logfile as L
logs = {name.lstrip('H2O-1-'): L(join('H2O-single','Calculations','log-'+name+'.yaml'))
for name in run_names}
Density of States
from BigDFT.DoS import DoS
dos = DoS(logfiles_dict=logs)
dos.plot().legend(loc='best')
Where we see the localization effect of the Hybrid functional with respect to the energies of the semilocal ones.
Total Dipole
We do a pretty printing via dataframes, just for fun. This is optional of course.
def dipole_info(log):
data = {coord+ ' (AU)': log.dipole[i] for i, coord in enumerate(['x', 'y', 'z'])}
data['Norm (Debye)'] = log.log['Electric Dipole Moment (Debye)']['norm(P)']
return data
from pandas import DataFrame
data = {name: dipole_info(log)for name, log in logs.items()}
table = DataFrame(data).T
table.round(2)
Here we see the polarizability effect induced by the environment.
Time to solution
We can also inspect the different timings for this runs, to understand how the computational resources have been used in the different approximations
times = {name: join('H2O-single','Calculations','time-'+name+'.yaml')
for name in run_names}
def plot_timings(times):
from futile.Time import TimeData
from matplotlib import pyplot as plt
tt=TimeData(*times.values(),static=True,plottype='Seconds')
#tweak the name of the labels
fig=plt.gcf()
ax=fig.get_axes()
_=ax[0].set_xticklabels(times.keys(),rotation=90)#,size=self.fontsize)
plot_timings(times)
We see that the PCM computation and the PBE0 calcuations take the same amount of time for the SCF cycle. Those data are only partially considering the time for the entire run which also takes into account the initialisation and the finalization
def get_total_time_info(time):
import yaml
with open(time) as ifile:
dt = yaml.load(ifile,Loader=yaml.Loader)
data = {k: v[1] for k, v in dt['SUMMARY'].items()}
data.update({'WFN_OPT %': dt['SUMMARY']['WFN_OPT'][0]})
return data
table=DataFrame({name:get_total_time_info(time) for name, time in times.items()}).T
table.round(3)
GPU acceleration
We have the possiblity of accelerating the computation using GPU. We start by performing a rather insignificant accelration on a single molecule, and then we will move to the larger 32 H2O system. Let’s increment our dataset:
#identify a name for the run
run_name = 'H2O1-PBE0-GPU'
# we then have the arguments of the function to submit
func_kwargs = dict(run_name=run_name, filename=filename, functional='PBE0')
run_args = {}
run_args['gpu'] = True # enable GPU
run_args['jobname'] = run_name
run_args['mpi'] = 1
run_args['omp'] = 4
run_args['time'] = 2
all_runs.append_run(args=func_kwargs,
extra_files_recv=[run_name+'.tar.gz'],
**run_args)
all_runs.run()
all_runs.run_cmds
all_runs.results
Exercise: perform the same analysis as above to verify what has changed.
Run of the full system
We will now crate a dataset of larger runs with the system identified from the full molecule. We will run the calculation with:
PBE, CPU only PBE0, CPU only PBE0, GPU accelerated We will employ one single node of the VEGA architecture, knowing that the two architectures are different
local_dir = 'H2O-many'
remote_dir = f'/ceph/hpc/data/d2021-135-users/bigdft-school/{user}/H2O-many' #use your own directory for new calculations
filename = 'H2O-32.xyz'
large_runs = Dataset(function = calculate_single_point,
url = url,
name = 'H2O-many',
extra_files_send=filename,
remote_dir = remote_dir,
local_dir = local_dir,
block_reinit = force_reinitialization)
#we collect the run names for future reuse
large_run_names = []
# the we append the runs we would like to perform
for functional, gpu in [('PBE', False), ('PBE0', False), ('PBE0', True)]:
#identify a name for the run
run_name = '-'.join([filename.split('.')[0],functional] + (['GPU'] if gpu else []))
# we then have the arguments of the function to submit
func_kwargs = dict(run_name=run_name, filename=filename, functional=functional, gpu=gpu)
# the full system has 128 KS orbitals
run_args = {}
run_args['jobname'] = run_name
run_args['mpi'] = 32 if gpu else 128
run_args['omp'] = 4 if gpu else 2 # the cpu node has 256 cores
run_args['time'] = '00:10:00'
run_args['gpu'] = gpu
large_runs.append_run(args=func_kwargs,
extra_files_recv=[run_name+'.tar.gz'],
**run_args)
large_run_names.append(run_name)
#for slow connections
# large_runs.url.timeout=120
# large_runs.url.max_timeouts=20
# large_runs.url.raise_errors=False #to override ssh error messages
large_runs.run()
large_runs.run_cmds
large_runs.url.cmd('squeue --me')
large_runs.fetch_results()
large_energies = large_runs.results
large_energies
extract_results('H2O-many')
large_logs = {name.lstrip('H2O-32-'): L(join('H2O-many','Calculations','log-'+name+'.yaml'))
for name in reversed(large_run_names)}
large_dos = DoS(logfiles_dict=large_logs)
large_dos.plot().legend(loc='best')
reordered = [large_run_names[i] for i in [1,2,0]]
large_times = {name: join('H2O-many','Calculations','time-'+name+'.yaml')
for name in reordered}
plot_timings(large_times)
We see that with the usage of one single node, we get similar timings between the GPU-accelerated PBE0 and the traditionale PBE calculation on a GPU architecture