# DQC: a Python program package for Differentiable Quantum Chemistry

M. F. Kasim  Present address: Machine Discovery Ltd, muhammad@machine-
discovery.com Department of Physics, Clarendon Laboratory, University of
Oxford, Parks Road, Oxford OX1 3PU, UK  S. Lehtola  Molecular Sciences
Software Institute, Blacksburg, Virginia 24061, United States  S. M. Vinko
<EMAIL_ADDRESS>Department of Physics, Clarendon Laboratory,
University of Oxford, Parks Road, Oxford OX1 3PU, UK  Central Laser Facility,
STFC Rutherford Appleton Laboratory, Didcot OX11 0QX, UK

###### Abstract

Automatic differentiation represents a paradigm shift in scientific
programming, where evaluating both functions and their derivatives is required
for most applications. By removing the need to explicitly derive expressions
for gradients, development times can be be shortened, and calculations
simplified. For these reasons, automatic differentiation has fueled the rapid
growth of a variety of sophisticated machine learning techniques over the past
decade, but is now also increasingly showing its value to support ab initio
simulations of quantum systems, and enhance computational quantum chemistry.
Here we present an open-source differentiable quantum chemistry simulation
code, DQC, and explore applications facilitated by automatic differentiation:
(1) calculating molecular perturbation properties; (2) reoptimizing a basis
set for hydrocarbons; (3) checking the stability of self-consistent field wave
functions; and (4) predicting molecular properties via alchemical
perturbations.

##  I Introduction

Automatic differentiation is a collection of techniques used to evaluate, up
to machine precision, the derivative of a function specified by a computer
program. It allows software developers to focus solely on designing the best
model for a given problem, without having to worry about implementing any
derivatives of the model with respect to its various mathematical parameters.
It has already had a transformative effect in machine learning, enabling the
development of many new techniques over the past decade, such as batch
normalization ioffe2015batchnorm , attention layers vaswani2017attention , and
unique neural network architectures schutt2018schnet ; ronneberger2015unet .

Automatic differentiation is still relatively new in the context of
computational sciences, but is already showing promise across a diverse set of
applications including tensor networks liao2019differentiable-tensor-network ,
computational fluid dynamics schenck2018spnets-diff-cfd , and molecular
dynamics simulations schoenholz2020jaxmd . Automatic differentiation is also
growing increasingly popular in quantum chemistry, where it has been used to
optimize molecular basis sets tamayo2018automatic-diffiqult , to calculate
higher derivatives of various exchange-correlation (xc) functionals
ekstrom2010arbitrary-xc  of density functional theory (DFT)
hohenberg1964inhomogeneous ; kohn1965self-ks , and to determine arbitrary-
order nuclear coordinate derivatives of electronic energies
abbott2021arbitrary-geom-qc .

Automatic differentiation is also an essential stepping stone to enable direct
integration of quantum chemistry methods with machine learning models and
their training. In this context, a differentiable implementation of DFT was
recently used to learn the xc functional li2021kohn-sham-regularizer  from
accurate reference calculations within the density matrix renormalization
group approach, or from a mixture of computational and experimental data
kasim2021learning-xc , showing a promising new approach to developing
transferable and robust xc functionals via deep learning.

Although differentiation in quantum chemistry can be done via finite-
difference schemes, calculating the gradient of a parameter with respect to
some other input parameter can be time consuming as one has to run the
simulation as many times as the number of input parameters. Moreover, finite-
difference methods are prone to numerical instabilities and are very sensitive
to the chosen step-size. Automatic differentiation addresses these challenges
by calculating the analytical gradient via the chain rule, eliminating both
the need to run the simulation many times, and the need for step-size tuning.

To address the growing need for automatic differentiation in quantum
chemistry, we introduce Differentiable Quantum Chemistry (DQC), a DFT and
Hartree–Fock (HF) hartree1928wave  simulation code. DQC is written in Python
using PyTorch neurips2019-pytorch  and xitorch kasim2020xitorch . While
PyTorch provides gradient calculations for elementary operations such as
matrix multiplication and explicit eigendecomposition, xitorch provides
gradient calculations for functional operations such as root finding,
optimization, and implicit eigendecomposition. The use of PyTorch and xitorch
in DQC enables various applications in quantum chemistry that would otherwise
be far more difficult to pursue. For example, the work by Kasim and Vinko
kasim2021learning-xc  on learning the xc functional directly from a set of
heterogeneous experimental data and calculated density profiles within the
constraints of Kohn–Sham DFT already mentioned above was enabled by the use of
DQC.

We begin this paper by outlining the basic theory behind DQC in section II.
The implementation is described in section III. Applications that exemplify
the present approach are given in section IV. We discuss the challenges in the
use of automatic differentiation techniques in computational science in
section V, and summarize our work in section VI.

##  II Methods

Quantum chemical calculations of the electronic structure typically require
the evaluation of abstract functionals, such as root finding for self-
consistent field (SCF) iterations, and minimization for geometry optimizations
and the direct minimization approach to the SCF problem. We discuss the
handling of these functionals in the context of DQC in this section.

###  II.1 SCF iterations

DQC is based on the use of a Gaussian basis set within the linear combination
of atomic orbitals approach (LCAO). For simplicity, we will only present the
theory for the spin-restricted formalism, as the spin-unrestricted (and
restricted open-shell) formalisms are analogous. The SCF iterations proceed in
the LCAO approach by solving Roothaan’s equation roothaan1951new ;
lehtola2020overview

$\mathbf{F(D)C}=\mathbf{SCE},$ (1)

where $\mathbf{F}(\mathbf{D})$ is the Fock matrix which is a function of the
density matrix $\mathbf{D}$, $\mathbf{C}$ is the orbital matrix, $\mathbf{S}$
is the overlap matrix, and $\mathbf{E}$ is a diagonal matrix containing the
orbital energies. The generalized eigenvalue equation in (LABEL:roothaan) can
be reduced to the normal form by orthogonalizing the overlap matrix
$\mathbf{S}$ lowdin1956quantum , and operating in the linearly independent
basis spanned by the eigenvectors of $\mathbf{S}$ with large enough
eigenvalues; any ill-conditioning in the overlap matrix can be removed by
choosing a well-defined sub basis with the help of a pivoted Cholesky
decomposition lehtola2019curing . The density matrix $\mathbf{D}$, required to
construct the Fock matrix $\mathbf{F}$, can be obtained by

$\mathbf{D}=\mathbf{CNC^{\dagger}}$ (2)

with $\mathbf{N}$ a diagonal matrix containing the occupation numbers of the
orbitals. As discussed in ref. 19, the SCF procedure requires repeatedly
solving (LABEL:roothaan,_dm) until self-consistency is achieved.

The Roothaan iteration in (LABEL:roothaan) can be written mathematically as
the process of finding a vector $\mathbf{y}$ such that

$\mathbf{y}=\mathbf{f}(\mathbf{y},\boldsymbol{\theta}),$ (3)

where $\mathbf{f}$ is a function that takes the vector $\mathbf{y}$ and other
parameters $\boldsymbol{\theta}$, and returns the expected value of vector
$\mathbf{y}$. In DQC, the parameter $\mathbf{y}$ is represented by the Fock
matrix $\mathbf{F}$, so the function $\mathbf{f}$ is the procedure that solves
equation (1), calculating the density matrix from equation (2) and
constructing back the Fock matrix from the density matrix. The parameters
$\boldsymbol{\theta}$ represent the other variables involved in the
calculation, such as the overlap matrix $\mathbf{S}$ and the occupation number
matrix $\mathbf{N}$. The algorithm and gradient calculation for
(LABEL:roothaan,_self-consistent) are already available in PyTorch and
xitorch.

###  II.2 Direct minimization

An alternative approach to solving the self-consistent field equations is to
directly find the orbitals that minimize the total energy $\mathcal{E}$ by the
use of optimization algorithms head1988optimization . The energy $\mathcal{E}$
can be calculated from the density matrix $\mathbf{D}$ which can be obtained
in turn from the orbital coefficients $\mathbf{C}$ using (LABEL:dm). This
relation makes the energy a function of the orbital matrix,
$\mathcal{E}(\mathbf{C})$. However, as the orbitals must remain orthonormal,
$\mathbf{C^{\dagger}SC}=\mathbf{I}$, we introduce a new variable $\mathbf{Q}$,
defined in terms of its relation with $\mathbf{C}$ as

$\displaystyle\mathbf{Q}$ $\displaystyle=\mathbf{\hat{Q}R}$ (4)
$\displaystyle\mathbf{C}$ $\displaystyle=\mathbf{S}^{-1/2}\mathbf{\hat{Q}}.$
(5)

uq-relation-1 is the QR decomposition of the matrix $\mathbf{Q}$ into an
orthogonal matrix $\mathbf{\hat{Q}}$ and an upper triangular matrix
$\mathbf{R}$. uq-relation-2 involves the inverse square root of the overlap
matrix $\mathbf{S}$, which can be computed using the eigendecomposition of
$\mathbf{S}$. The energy can then be parameterized by $\mathbf{Q}$, reducing
the direct minimization problem to an unbounded optimization problem:

$\mathbf{Q^{*}}=\arg\min_{\mathbf{Q}}\mathcal{E}(\mathbf{Q}).$ (6)

The gradient of the energy with respect to $\mathbf{Q}$, i.e.
$\partial\mathcal{E}/\partial\mathbf{Q}$, is required for an efficient
solution. It is automatically computed by PyTorch and xitorch.

Once the optimum $\mathbf{Q}^{*}$ in equation (6) is found, $\mathbf{Q}^{*}$
can still be differentiated with respect to any other variables, such as the
nuclear coordinates or the occupation number matrix $\mathbf{N}$. This is made
possible by the gradient of the optimization functional provided by xitorch.

We note that our approach of using QR decomposition is slightly different from
common direct minimization techniques that use the matrix exponential of a
skew-Hermitian matrix such as the geometric direct minimization algorithm of
ref. 23.

Table 1: Execution speed comparison between DQC (SCF iterations) and PySCF. Cases | DQC | PySCF
---|---|---
H2O (HF/cc-pVDZ) | $96\ \mathrm{ms}$ | $245\ \mathrm{ms}$
H2O (PW92/cc-pVDZ) | $530\ \mathrm{ms}$ | $430\ \mathrm{ms}$
C4H5N (HF/cc-pVTZ) | $108\ \mathrm{s}$ | $17\ \mathrm{s}$
C4H5N (PW92/cc-pVTZ) | $101\ \mathrm{s}$ | $25\ \mathrm{s}$
C4H5N (density fit PW92/cc-pVTZ) | $30\ \mathrm{s}$ | $22\ \mathrm{s}$
C6H8O6 (density fit PW92/cc-pVDZ) | $87\ \mathrm{s}$ | $57\ \mathrm{s}$

##  III Implementation

DQC implements restricted and unrestricted HF and Kohn–Sham DFT calculations
without periodic boundary conditions. The energy can be minimized using either
SCF iterations with Broyden’s good method broyden1965class  for Fock matrix
updates, or with direct minimization using gradient descent with momentum
pearlmutter1991gd-momentum ; more elaborate SCF convergence accelerators will
be implemented at a later stage of the project. All parameters are
differentiable with respect to any other parameter present in the calculation,
including nuclear coordinates, atomic numbers, electron occupation numbers, as
well as the basis set exponents and contraction coefficients. The
differentiability of these parameters allows for the exploration of new
applications with DQC, a few of which are presented in what follows.

Although execution speed is not a top priority for DQC, the program shows good
overall performance. speed compares the running times of DQC with PySCF
sun2020recent-pyscf , an established quantum chemistry code. For small
systems, we find that DQC is as efficient as PySCF. However, for moderate-
sized molecules (e.g. C4H5N), DQC is about 4–6 times slower than PySCF. This
slowdown can be attributed to DQC currently not taking advantage of the
symmetry of the 2-electron integral tensor, which can reduce the tensor size
by roughly a factor of 8 when the basis functions are real-valued. In
contrast, DQC is only about 50% slower than PySCF for systems with density
fitting. The difference is mainly caused by the higher number of SCF
iterations required in DQC over PySCF, due to the use of different SCF
convergence acceleration schemes in the two codes.

Figure 1: Unrestricted PW92/cc-pVDZ perdew1992accurate-pw92 ;
dunning1989gaussian-ccbasis  energy for the hydrogen molecule cation
H${}_{2}^{+}$ as a function of the internuclear distance calculated via the
SCF (using Broyden’s method broyden1965class ) and direct minimization
approaches. Note that the SCF iterations do not converge for internuclear
distances larger than 17 Bohr.

##  IV Applications

###  IV.1 Direct minimization

The automatic availability of gradients makes it easy to implement the direct
minimization method in DQC. Direct minimization is known to be more robust
than SCF iterations in finding a converged solution lehtola2020overview , and
is particularly useful in difficult cases such as for molecules near their
dissociation limits. We illustrate this in variational-vs-scf, where we show
the PW92/cc-pVDZ perdew1992accurate-pw92 ; dunning1989gaussian-ccbasis  total
energy of H${}_{2}^{+}$ as a function of the internuclear distance. The SCF
method does not converge for internuclear distances greater than around 17
Bohr. In contrast, the direct minimization method continues to converge even
for substantially larger distances in excess of 40 Bohr.

###  IV.2 Checking SCF stability

One difficulty in SCF calculations is that the calculation can converge onto a
saddle point, which corresponds to an excited state lehtola2020overview . If
the calculation has converged onto a local minimum, the Hessian of the energy
with respect to orbital rotations must be positive semidefinite. Saddle point
solutions, instead, correspond to one or more negative eigenvalues of the
orbital Hessian. To maintain the orthonormality of the orbitals, the Hessian
is calculated based on the $\mathbf{Q}$ variable. Conveniently, we do not need
to derive the expression for the Hessian matrix
$\partial^{2}\mathcal{E}/\partial\mathbf{Q}^{2}$ as it is automatically
obtained by PyTorch.

Moreover, only the lowest eigenvalue needs to be calculated for the stability
check, so the full Hessian matrix does not need to be constructed. We obtain
the lowest eigenvalue using Davidson’s iterative algorithm
davidson1975iterative , as implemented in xitorch. Note that the gradients of
the lowest energy eigenvalue with respect to all other parameters in the
calculation are automatically available, which may prove useful in further
future applications.

Table 2: Lowest eigenvalue of the HF/pc-1 orbital Hessian matrix. All of the calculations correspond to stationary points of the HF energy. Molecules | Distance | Energy | Lowest eigenvalue
---|---|---|---
| (Bohr) | ($E_{h}$) | ($E_{h}$)
O2 (ground) | 2.0 | $-149.54$ | $-6\times 10^{-14}$
O2 (excited) | 2.0 | $-149.17$ | $-0.73$
BeH (ground) | 2.5 | $-15.14$ | $-4\times 10^{-14}$
BeH (excited) | 2.5 | $-15.01$ | $-0.25$
CH (ground) | 2.0 | $-38.259$ | $-5\times 10^{-8}$
CH (excited) | 2.0 | $-38.256$ | $-0.07$

We show some examples of SCF stability checks for HF/pc-1 jensen2001
calculations on diatomic molecules in scf-stability. As can be seen from the
data, the lowest eigenvalues of the orbital Hessian are very close to zero for
the ground state, while the excited states yield large negative eigenvalues.
This illustrates that it is straightforward to check whether or not a state
corresponds to a true minimum with DQC.

###  IV.3 Molecular properties

A key advantage of writing quantum chemistry software with automatic
differentiation is that calculations of molecular properties can be
implemented efficiently: all we need to know is how a specific property
relates to some derivative expression. For example, vibrational modes and
frequencies of a molecule can be written as

$\mathbf{q}_{\rm vib},\omega_{\rm
vib}=\mathrm{eig}\left(\frac{\partial^{2}\mathcal{E}}{\partial\mathbf{X}^{2}}\right),$
(7)

where $\mathbf{X}$ is an $n\times 3$ matrix containing the nuclear coordinates
of the $n$ atoms, $\mathrm{eig}$ is the eigendecomposition, $\mathbf{q}_{\rm
vib}$ is one of the vibrational modes, and $\omega_{\rm vib}$ is its
frequency. The expression shown in (LABEL:vibr) is all that is needed to
calculate the vibrational characteristics of the molecule; the explicit form
of the derivative expression is not required.

Table 3: Perturbative properties of H2O from HF/cc-pVDZ calculations. The middle column presents the values calculated in DQC while the last column shows the CCCBDB values NIST_CCCBDB . The IR intensities and Raman activities are presented at the frequency of 1800 cm-1. Properties | DQC | CCCBDB
---|---|---
IR intensities (km/mol) | 80.69 | 80.70
Raman activities (A4/amu) | 4.79 | 4.79
Dipole (D) | $-2.044$ | $-2.044$
Quadrupolexx (DA) | $-7.008$ | $-7.008$

Another useful example is the calculation of the electric dipole and
quadrupole moments of a molecule. The electric dipole moment is given by

$\boldsymbol{\mu}=-\frac{\partial\mathcal{E}}{\partial\mathbf{E}}+\sum_{i}Z_{i}\mathbf{x}_{i},$
(8)

while the electric quadrupole tensor is given by

$\mathbf{M}=-2\frac{\partial\mathcal{E}}{\partial\nabla\mathbf{E}}+\sum_{i}Z_{i}\mathbf{x}_{i}\mathbf{x}_{i}^{T}.$
(9)

In both expressions, $\mathcal{E}$ is the total energy of the molecule,
$\mathbf{E}$ is the electric field, $\mathbf{Z}_{i}$ is the atomic number of
$i$-th atom, and $\mathbf{x}_{i}$ are its coordinates.

Having these vibrational and electric multipole properties readily available
makes obtaining the Raman vibrational spectrum straightforward. For example,
the intensity of the infrared vibrational transition for a normal mode
$\mathbf{q}$ at a frequency $\omega$ is given by

$I_{\rm IR}=\sum_{i}\left(\sum_{jk}\frac{\partial\mu_{i}}{\partial
X_{jk}}q_{jk}\right)^{2},$ (10)

where $\mathbf{\mu}$ is the electric dipole moment and $\mathbf{X}$ is the
matrix of atomic coordinates.

Similarly, the Raman activity at the same frequency and normal mode is
proportional to p2007analytic-raman

$\displaystyle I_{\rm Raman}$
$\displaystyle=5\left[\mathrm{tr}\left(\boldsymbol{\alpha}\right)\right]^{2}+7\gamma;$
(11) $\displaystyle\alpha_{ij}$
$\displaystyle=\sum_{kl}\frac{\partial^{2}\mu_{i}}{\partial E_{j}\partial
X_{kl}}q_{kl};$ (12) $\displaystyle\gamma$
$\displaystyle=\sum_{ij}(1-\delta_{ij})\left[\frac{1}{4}(\alpha_{ii}-\alpha_{jj})^{2}+\frac{3}{2}\alpha_{ij}^{2}\right],$
(13)

where $\delta_{ij}$ is the Kronecker delta.

Figure 2: Basis set truncation errors for the cc-pVDZ and cc-pVTZ basis sets
for a range of hydrocarbons, compared with the truncation error for a
reoptimized cc-pVDZ basis set. The reference energies are computed in the cc-
pV5Z basis set. None of the molecules shown in the figure were used in the
basis set optimization.

The HF/cc-pVDZ perturbation properties of H2O, found using the expressions
above, are displayed in Table 3. The values are in excellent agreement with
data from the the computational chemistry comparison and benchmark database
(CCCBDB) NIST_CCCBDB , even though DQC does not have any explicit code to
calculate the gradients required for these properties.

###  IV.4 Basis set optimization

The differentiability of DQC with respect to the basis set parameters enables
the optimization of system-specific basis sets. Here we show this capability
by optimizing a basis set for hydrocarbons within Kohn–Sham DFT using the PW92
functional. We reoptimized the cc-pVDZ basis dunning1989gaussian-ccbasis  for
a training set of molecules consisting of CH (methylidyne), CH3 (methyl
radical), CH4 (methane), C2H2 (acetylene), and C2H4 (ethylene). The accuracy
of the reoptimized basis was then tested on a set of hydrocarbons that were
not included in the training set. The results are shown in basis-opt. The
reoptimization of the cc-pVDZ basis set leads to a marked decrease of the
total energies of all the molecules in the test set, as the cc-pVXZ basis set
series is designed to capture correlation energies instead of polarization
energies dunning1989gaussian-ccbasis , and as hydrocarbons are chemically
similar.

###  IV.5 Alchemical perturbation

One of the differentiable quantities in DQC is the atomic number, allowing us
to perform alchemical perturbation studies to predict the properties of
molecules without actually needing to simulate them balawender2019exploring-
alchemical ; von2020alchemical-dft . As a simple example, we will show here
that some properties of diatomic molecules CO (atomic numbers 6 and 8) and BF
(atomic numbers 5 and 9) can be estimated directly from the properties of N2
(atomic number 7) and its alchemical perturbations. To do this, we parametrize
the atomic number of the atoms in the diatomic molecule by a continuous
variable $\lambda$, so that the atoms have atomic numbers $Z_{l}=7+\lambda$
and $Z_{r}=7-\lambda$. The molecules N2, CO, and BF thus correspond to
$\lambda=0$, $\lambda=1$, and $\lambda=2$, respectively.

The properties we aim to predict are the bond length $s^{*}$ and the energy
$E^{*}$ at the equilibrium position, which can be mathematically expressed as

$\displaystyle s^{*}(\lambda)$ $\displaystyle=\arg\min_{s}E(s,\lambda),$ (14)
$\displaystyle E^{*}(\lambda)$ $\displaystyle=E(s^{*},\lambda).$ (15)

Performing HF/pc-1 calculations, we evaluate the equilibrium distance and the
energy at the equilibrium position in two ways. The first way is to optimize
the geometry for various fixed values of $\lambda$, and calculate the above
properties directly. The calculations were performed separately, but with the
same basis (pc-1 nitrogen basis on all atoms) for different molecules, in
analogy to the procedure used in refs 33 and 34. The second way is to estimate
those properties using a Taylor series expansion to second order in $\lambda$:

$\displaystyle s^{*}(\lambda)$ $\displaystyle\approx
s^{*}(0)+\lambda\frac{\partial
s^{*}}{\partial\lambda}(0)+\frac{1}{2}\lambda^{2}\frac{\partial^{2}s^{*}}{\partial\lambda^{2}}(0)$
(16) $\displaystyle E^{*}(\lambda)$ $\displaystyle\approx
E^{*}(0)+\lambda\frac{\partial
E^{*}}{\partial\lambda}(0)+\frac{1}{2}\lambda^{2}\frac{\partial^{2}E^{*}}{\partial\lambda^{2}}(0).$
(17)

As $s^{*}$ and $E^{*}$ are calculated at equilibrium, calculating the
perturbation terms requires propagating the gradient through the optimization
process. However, automatic differentiation makes the propagation trivial, as
it is automatically handled by the optimization routine in xitorch.

Figure 3: Comparison of the properties of CO and BF obtained via an exact
calculation and an estimation from the second-order gradients via alchemical
perturbations, employing the nitrogen pc-1 basis on all atoms.

The results obtained via these two approaches are compared in alchemy. As we
can see from these results, the properties of CO and BF can be estimated
accurately from alchemical perturbations of N2. The estimated equilibrium
distances for CO and BF differ by $-0.0004$ and 0.024 Bohr, respectively,
while the estimated equilibrium energies deviate by 33 and 587 m$E_{h}$,
respectively. This shows that the equilibrium position of new molecules can be
estimated well with the alchemical gradient calculated by DQC, without
actually having to calculate those molecules.

##  V Discussion

Implementing quantum chemistry with automatic differentiation libraries is a
promising way to accelerate simulation workflows and to enable novel
applications. However, the implementation comes with several challenges. An
overarching challenge is that the automatic differentiation library used here,
PyTorch, is primarily designed for deep learning rather than for scientific
computing. As deep learning only focuses on low-order derivatives, accessing
high-order gradients that are commonly required for scientific applications
can be difficult.

Detecting numerical instabilities in high-order gradient calculations can also
be demanding. Instabilities that produce NaN (not-a-number) in PyTorch are
relatively straightforward to manage with its debugging feature since version
1.8, but other instabilities that do not produce a NaN can be challenging to
detect.

Another challenge is debugging and profiling higher-order gradient
calculations. As the gradient is generated automatically, it is hard to find
the slowest part of the code, or the part that requires the most memory,
because this information is not readily provided by the available profiling
tools.

Besides higher-order gradient calculations, using automatic differentiation
libraries also poses unique challenges in terms of code optimization. For
example, quantum chemistry codes usually save the two-electron integrals on
disk due to their large size and process them only in blocks small enough to
fit easily into memory. This scheme can only be used in DQC if the gradients
with respect to the nuclear positions and the basis are not required. To the
authors’ best knowledge, there is currently no obvious structure in PyTorch to
allow gradient-propagating operations to work with large tensors saved on
disk.

##  VI Conclusions

Implementing quantum chemical calculations using automatic differentiation
liberates us from needing to derive analytical gradient expressions. With
gradients automatically generated by the program, software developers can
focus on designing better, more detailed computational models, and on applying
them to problems at hand. We have shown how automatic differentiation allows
us to easily explore various applications in quantum chemistry, and are
confident that further exploration of this approach will unveil new
applications that have not been considered to date.

## Code availability

The Differentiable Quantum Chemistry (DQC) code can be found at
https://github.com/diffqc/dqc/. The repository that contains the applications
presented in this paper is located at https://github.com/diffqc/dqc-apps/.

###### Acknowledgements.

M.F.K. and S.M.V. acknowledge support from the UK EPSRC grant EP/P015794/1 and
the Royal Society. S.M.V. is a Royal Society University Research Fellow. The
authors declare no conflict of interest.

## References

  * [1] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In International conference on machine learning, pages 448–456. PMLR, 2015.
  * [2] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  arXiv preprint arXiv:1706.03762, 2017.
  * [3] Kristof T Schütt, Huziel E Sauceda, P-J Kindermans, Alexandre Tkatchenko, and K-R Müller.  Schnet–a deep learning architecture for molecules and materials.  J. Chem. Phys., 148(24):241722, 2018.
  * [4] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * [5] Hai-Jun Liao, Jin-Guo Liu, Lei Wang, and Tao Xiang.  Differentiable programming tensor networks.  Physical Review X, 9(3):031041, 2019.
  * [6] Connor Schenck and Dieter Fox.  Spnets: Differentiable fluid dynamics for deep neural networks.  In Conference on Robot Learning, pages 317–335. PMLR, 2018.
  * [7] Samuel Schoenholz and Ekin Dogus Cubuk.  Jax md: a framework for differentiable physics.  Advances in Neural Information Processing Systems, 33, 2020.
  * [8] Teresa Tamayo-Mendoza, Christoph Kreisbeck, Roland Lindh, and Alán Aspuru-Guzik.  Automatic differentiation in quantum chemistry with applications to fully variational hartree–fock.  ACS central science, 4(5):559–566, 2018.
  * [9] Ulf Ekström, Lucas Visscher, Radovan Bast, Andreas J Thorvaldsen, and Kenneth Ruud.  Arbitrary-order density functional response theory from automatic differentiation.  Journal of chemical theory and computation, 6(7):1971–1980, 2010\.
  * [10] Pierre Hohenberg and Walter Kohn.  Inhomogeneous electron gas.  Phys. Rev., 136(3B):B864, 1964.
  * [11] Walter Kohn and Lu Jeu Sham.  Self-consistent equations including exchange and correlation effects.  Phys. Rev., 140(4A):A1133, 1965.
  * [12] Adam S Abbott, Boyi Z Abbott, Justin M Turney, and Henry F Schaefer III.  Arbitrary-order derivatives of quantum chemical methods via automatic differentiation.  The Journal of Physical Chemistry Letters, 12(12):3232–3239, 2021\.
  * [13] Li Li, Stephan Hoyer, Ryan Pederson, Ruoxi Sun, Ekin D Cubuk, Patrick Riley, and Kieron Burke.  Kohn-Sham equations as regularizer: Building prior knowledge into machine-learned physics.  Phys. Rev. Lett., 126(3):036401, 2021.
  * [14] Muhammad F Kasim and Sam M Vinko.  Learning the exchange-correlation functional from nature with fully differentiable density functional theory.  arXiv preprint arXiv:2102.04229, 2021.
  * [15] Douglas Rayne Hartree.  The wave mechanics of an atom with a non-coulomb central field. part ii. some results and discussion.  In Mathematical Proceedings of the Cambridge Philosophical Society, volume 24, pages 111–132. Cambridge University Press, 1928.
  * [16] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32, pages 8026–8037. Curran Associates, Inc., 2019.
  * [17] Muhammad F Kasim and Sam M Vinko.  $\xi$-torch: differentiable scientific computing library.  arXiv preprint arXiv:2010.01921, 2020.
  * [18] Clemens Carel Johannes Roothaan.  New developments in molecular orbital theory.  Rev. Mod. Phys., 23(2):69, 1951.
  * [19] Susi Lehtola, Frank Blockhuys, and Christian Van Alsenoy.  An overview of self-consistent field calculations within finite basis sets.  Molecules, 25(5):1218, 2020.
  * [20] Per-Olov Löwdin.  Quantum theory of cohesive properties of solids.  Adv. Phys., 5(17):1–171, 1956.
  * [21] Susi Lehtola.  Curing basis set overcompleteness with pivoted Cholesky decompositions.  J. Chem. Phys., 151(24):241102, 2019.
  * [22] Martin Head-Gordon and John A Pople.  Optimization of wave function and geometry in the finite basis Hartree–Fock method.  J. Phys. Chem., 92(11):3063–3069, 1988.
  * [23] Troy Van Voorhis and Martin Head-Gordon.  A geometric approach to direct minimization.  Molecular Physics, 100(11):1713–1721, 2002.
  * [24] Charles G Broyden.  A class of methods for solving nonlinear simultaneous equations.  Mathematics of computation, 19(92):577–593, 1965.
  * [25] Barak A Pearlmutter.  Gradient descent: Second-order momentum and saturating error.  Advances in neural information processing systems, pages 887–894, 1991.
  * [26] Qiming Sun, Xing Zhang, Samragni Banerjee, Peng Bao, Marc Barbry, Nick S Blunt, Nikolay A Bogdanov, George H Booth, Jia Chen, Zhi-Hao Cui, et al.  Recent developments in the PySCF program package.  J. Chem. Phys, 153(2):024109, 2020.
  * [27] John P Perdew and Yue Wang.  Accurate and simple analytic representation of the electron-gas correlation energy.  Phys. Rev. B, 45(23):13244, 1992.
  * [28] Thom H Dunning Jr.  Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen.  J. Chem. Phys., 90(2):1007–1023, 1989.
  * [29] E. R. Davidson.  The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real symmetric matrices.  Journal of Computational Physics, 17(1):87–94, 1975.
  * [30] Frank Jensen.  Polarization consistent basis sets: Principles.  J. Chem. Phys., 115:9113–9125, 2001.
  * [31] Editor: R. D. Johnson III.  NIST Computational Chemistry Comparison and Benchmark Database. NIST Standard Reference Database Number 101 (Release 21, August 2020), [Online]. Available: https://dx.doi.org/10.18434/T47C7Z. National Institute of Standards and Technology, Gaithersburg, MD., 2020.
  * [32] Darragh P. O’Neill, Mihály Kállay, and Jürgen Gauss.  Analytic evaluation of raman intensities in coupled-cluster theory.  Molecular Physics, 105(19-22):2447–2453, 2007.
  * [33] Robert Balawender, Michael Lesiuk, Frank De Proft, Christian Van Alsenoy, and Paul Geerlings.  Exploring chemical space with alchemical derivatives: alchemical transformations of h through ar and their ions as a proof of concept.  Physical Chemistry Chemical Physics, 21(43):23865–23879, 2019.
  * [34] Guido Falk von Rudorff and O Anatole von Lilienfeld.  Alchemical perturbation density functional theory.  Physical Review Research, 2(2):023220, 2020.


# The Wind in Our Sails: Developing a Reusable and Maintainable Dutch Maritime
History Knowledge Graph

Stijn Schouten  ,  Victor de Boer  Vrije Universiteit Amsterdam
Department of Computer ScienceDe Boelelaan 1111AmsterdamThe Netherlands
<EMAIL_ADDRESS><EMAIL_ADDRESS>,  Lodewijk Petram  Huygens
InstituutOudezijds Achterburgwal 185Amsterdamthe Netherlands
<EMAIL_ADDRESS>and  Marieke van Erp  KNAW Humanities Cluster
DHLabOudezijds Achterburgwal 185Amsterdamthe Netherlands
<EMAIL_ADDRESS>
(2021)

###### Abstract.

Digital sources are more prevalent than ever but effectively using them can be
challenging. One core challenge is that digitized sources are often
distributed, thus forcing researchers to spend time collecting, interpreting,
and aligning different sources. A knowledge graph can accelerate research by
providing a single connected source of truth that humans and machines can
query. During two design-test cycles, we convert four data sets from the
historical maritime domain into a knowledge graph. The focus during these
cycles is on creating a sustainable and usable approach that can be adopted in
other linked data conversion efforts. Furthermore, our knowledge graph is
available for maritime historians and other interested users to investigate
the daily business of the Dutch East India Company through a unified portal.

Knowledge Graph, Digital Humanities, Maritime history

††journalyear: 2021††copyright: rightsretained††conference: Proceedings of the
11th Knowledge Capture Conference; December 2–3, 2021; Virtual Event,
USA††booktitle: Proceedings of the 11th Knowledge Capture Conference (K-CAP
’21), December 2–3, 2021, Virtual Event, USA††doi:
10.1145/3460210.3493548††isbn: 978-1-4503-8457-5/21/12††ccs: Applied computing
Arts and humanities††ccs: Information systems Graph-based database models

##  1\. Introduction

Recent large-scale digitization efforts have increased the volume of digitally
available sources enormously (Claire Warwick et al., 2012). This is a huge
boon to digital humanities research, a multidisciplinary field, where
researchers utilize computational processing power to develop and present
insights that were not possible before in traditional humanities research
(Berry, 2012). However, using these digitized sources in large-scale research
is not always straightforward. When data sets are available as standalone
files, researchers have to go through the laborious process of collecting the
different files, interpreting and refining the data, defining relations
between different entities, and finally creating the queries and filtering
operations needed to answer their research questions. This can be a complex
process; for one thing, interpretation of variables may differ between
researchers (Haslhofer et al., 2018) and the process is often repeated by
every new researcher who wants to work with those data sets.

In this project, we create a knowledge graph of data related to the Dutch East
India Company (Vereenigde Oostindische Compagnie, abbreviated as VOC). The VOC
is often considered the first multinational company. It was a large Dutch
trading company that was founded in 1602 and remained in business until the
end of the eighteenth century (Gelderblom et al., 2013). The VOC kept records
of its ships, their cargo, and the sailors and soldiers who ventured to Asia.
This information is contained in four large data sets, three of which are
available online from the websites of the Dutch National Archive (VOC
Opvarenden, or VOC Crew members) and the Huygens Institute (Dutch Asiatic
Shipping, Bookkeeper General Batavia and Places). The last data set was
provided to us by the researchers from the Huygens Institute. The datasets are
described in more detail in Section 4. Making these data sets available in a
well-structured knowledge graph allows digital humanities scholars to
accelerate research with these rich and unique data (Hughes, 2012). While many
humanities researchers work with structured data such as databases and
spreadsheets, linked data is still uncharted territory for many. This is
partly due to unfamiliarity with the technology, and its at times steep
learning curve to start using it (Blanke et al., 2012).

The goal of this project is to design and develop a usable and sustainable
knowledge graph, that will be hosted by the Huygens Institute, the research
institute that manages three of the four data sets used in this project, and
that aims to drive data-driven research. We adopt an iterative design process;
in each iteration we capture additional knowledge from domain experts and the
data sets. We are not the first to create linked data for historical resources
(cf. (Blanke et al., 2012; Hyvönen et al., 2020)) or even the first to do this
for the historical maritime domain (cf. (de Boer et al., 2015)) However, our
work is the first to focus on sustainability of the knowledge graph, as well
as in bringing together four Dutch historical maritime data sets.

The challenges in our design process are that the approach should be
sustainable and that the designed artifact should be useful for researchers
and domain experts. Measuring and determining the most sustainable approach is
not straightforward and can be subjective. We therefore define four quality
attributes that can be used to measure, steer, and lead the design to a
sustainable solution. These design requirements are set in cooperation with
the project stakeholders and are as follows:

  * •

_Transparency:_ Clear and well-described design decisions and rules.

  * •

_Accessability:_ Results and intermediate results should be accessible.

  * •

_Re-usability:_ Re-usable methods and semantic models.

  * •

_Maintability:_ Future-proof knowledge graph and infrastructure.

We evaluate our knowledge graph using competency questions following (Azzaoui
et al., 2013).

##  2\. Related work

This project is on the intersection of computer science and humanities,
therefore we discuss relevant related work from both disciplines, starting
with best practices in knowledge modelling from the computer science/semantic
web domain.

The design guidelines of developing an accessible, transparent, reusable, and
maintainable knowledge graph are closely related to the FAIR (Wilkinson et
al., 2016) principles. Waagmeester et al. (Waagmeester et al., 2020) provide
us with tips and pointers on how to develop a knowledge graph according to the
FAIR principles. Dimou (Dimou, 2020) dives more into the technical side of
knowledge graph engineering and inspires further research and comparison of
these approaches. Finally, Jovanovik (Jovanovik, 2016) brings these together
in his linked data application development methodology, which focuses both on
FAIR principles and the technical implementation.

Our work is not the first to design a knowledge graph for the historical
maritime domain. De Boer et al. (de Boer et al., 2015) started in 2014 with
Dutch Ships and Sailors (DSS). The researchers involved in this project
connected four different data sets in a first effort to unify the scattered
Dutch historical maritime data. Additional work in this domain was carried out
by Entjes (Entjes, 2015) in 2015. Both projects score high in terms of
_transparency_ and _accessibility_ , with design and conversion decisions
being documented and reusable mapping scripts made available. In terms of re-
usability, the verdict is less clear. DSS uses elements from multiple general-
purpose ontologies and data models, and combines these with concepts and
classes from a specific DSS model to describe the maritime activities. This
design decision has definite upsides as it is generally recommended to apply
widely used models and ontologies to enable interoperability with other graphs
on the web.111https://www.w3.org/standards/semanticweb/ontology The re-use
potential within the (maritime) history domain is low, however: it is unlikely
that the dss:Ship concept is used between projects outside of the DSS sphere.

As a use case for the OntoME platform Beretta et al. (Beretta et al., 2019)
developed another ontology for Dutch maritime data on ships and sailors in
2019.222http://ontome.net OntoME is aimed at maintaining and distributing
ontologies for digital heritage projects, and the approach of Beretta et al.
therefore scores high in terms of maintainability. Furthermore, their ontology
is based on the CIDOC CRM (Cidoc, 2003), a widely accepted and highly detailed
model for digital cultural heritage. Concepts in this model, such as voyages,
have a better chance of being re-used in other semantic web projects, thus
offering great alignment possibilities. It is accessible and available for re-
use on the OntoME platform.333http://ontome.net/project/28 The transparency
and accessibility design principles of Beretta et al. and, more generally, the
OntoME platform, leave room for improvement. Design decisions are not openly
published, making it difficult for outsiders to understand and follow the
thought process during the design.

Each of the discussed projects has its strengths and weaknesses. We adapt some
of their strong aspects in our project, such as the transparent methodology of
DSS and the model maintainability of OntoME.

##  3\. Research method and design

###  3.1. Design decisions

In this project, the CIDOC CRM is selected as the ontology of choice. In
contrast to the aforementioned DSS project, this will allow for further
interoperability with other Digital Humanities Linked Data solutions. CIDOC
CRM is a de facto domain standard and a safe choice; as it is for example also
already used in the OntoME project.

We use GraphDB for data
conversion.444https://graphdb.ontotext.com/documentation/standard/release-
notes.html GraphDB is transparent due to its interface; users can directly
intuitively influence the mapping. It is reusable due to the ability to import
and export mapping tables between projects. Maintainability also scores high
due to the active development of the application. Other conversion tools were
tested and examined as well. These include but are not limited to the tools of
the Datalagend ecosphere555https://iisg.amsterdam/en/clariah which we found to
be well capable of quickly producing high volumes of RDF data with minimal
user input, but too rigid for modelling and linking entities. The RML-based
editors, we investigated such as RMLEditor666https://rml.io/tools/rmleditor/
places a limit on file size and requests for an extension were not answered.
Other RML editors (such as CARML,777https://github.com/carml/carml, and
Yarrrml (matey)888https://rml.io/yarrrml/matey/) lack a workable visual or GUI
representation. This downside makes the recommendation to use such an approach
difficult for users outside computer and information sciences. Therefore, we
deemed GraphDB the most suitable for our research goals.

###  3.2. Iterative design

We select design science as our central design methodology, as the primary
goal of our project is the development of an artifact: the knowledge graph. In
Design Science research, several research methods and frameworks exist (cf.
(Wieringa, 2014; Hevner et al., 2004)). The main takeaway of these theories is
the design-test cycle. In this iterative cycle, the researcher uses the
designed artifact to solve the given problem while grounding the artifact in a
socio-technical foundation. The artifact is evaluated and the evaluation
results are used in the second iteration of the design-test cycle. The cycle
continues until a certain predefined threshold is met.

We complete two design cycles. This allows for a review and evaluation period
while still being time-efficient. In the first design-test cycle, we implement
the approach of Jovanovik (Jovanovik, 2016) with the goal of re-using some of
the steps during the second design-test cycle. During the first design-test
cycle we also conduct an open and unstructured interview with a domain expert
from the Huygens Institute. The domain expert is responsible for some and
involved in most of the curation efforts of data sets used during in this
project. An interview does not only help with understanding the data sets but
also creates an understanding of the use case, how domain experts plan to use
the artifact, and what they expect from it. Not all questions can be asked and
answered during a single interview; some questions may arise during the cycle.
Due to the pandemic, in person collaborations at the research institute’s
office were not possible; we therefore keep close contact via email. Finally,
the documentation provided by the authors of the previous historical maritime
conversions ((de Boer et al., 2015; Entjes, 2015; Beretta et al., 2019)) is
examined during desk research.

The first review period has two goals. First, errors, shortcomings, and other
issues are investigated. Second, the review period enables the start and
direction of the second design iteration. An unstructured discussion session
with digital humanities experts and the domain expert completes the review
period.

Finally, in the second review period, we test the design again and conduct a
workshop with several stakeholders from the Huygens Institute and the
University of Amsterdam. In this workshop, the developed artifact can be user-
tested, thus validating the research results. Aside from the workshop, a short
data story is written to illustrate the usefulness of the created knowledge
graph (presented in Section 6).

We evaluate the _sustainability_ aspect of our project through the design
guidelines. _Usefulness_ can be measured with the completion of the competency
questions (Azzaoui et al., 2013). Our competency questions come from two
sources: the domain expert of the Huygens Institute and a related
project(Entjes, 2015). In the list of competency questions below, (D)
indicates that the domain expert submitted the question; questions derived
from the related project are denoted with (R):

  * •

(D) Which VOC Chamber was accountable for the largest number of slaves
transported?

  * •

(D) How were the shipping routes in Asia divided between the various VOC
chambers (for example, did ships from a specific chamber predominantly sail on
certain routes)?

  * •

(D) What was the average value of cargo on VOC return voyages per crew
member/ship’s ton, and how did this evolve over time?

  * •

(D) To what extent was the value per ships’ ton on return voyages correlated
with the skipper’s track record (how much experience, i.e., how many previous
voyages / how quickly did a ship get to Asia on an outbound voyage)?

  * •

(R) How did the price of goods change over time?

  * •

(R) Perform a search for specific information in large data sets, such as
names;

  * •

(R) Can external factors, such as war, be linked to the intensity of shipping?

##  4\. First iteration

As discussed in the methodology section, we derive the structure of the
iterations from the work of Jovanonik (Jovanovik, 2016).

Figure 1. Example of mapping branch

###  4.1. Data sets

We include three VOC related data sets in the first design iteration:

  * •

_Dutch-Asiatic Shipping (DAS)_ This data set contains information about all
voyages from the Netherlands to Asia and back made by the VOC and its
predecessors. Some core details are available for every voyage, such as the
VOC Chamber that outfitted and contracted the ship, and the places and dates
of departure and arrival.999http://resources.huygens.knaw.nl/das Retrieved on:
13-05-2021 During these voyages, events occurred, such as a storm or a capture
of an enemy vessel. These events are stored in the voyParticulars column.
Voyages had a master and were made by ships that have a name, type and
tonnage, and were built at shipyards. Ships could make multiple journeys and
sometimes have their name changed. DAS contains the information of 4,700
voyages from the Netherlands to Asia and of 3,400 return voyages.

  * •

_Bookkeeper-General Batavia (BGB)_ In Batavia (now Jakarta in Indonesia), the
bookkeeper and his clerks kept track of all imports and exports over sea. For
a sizable portion of the eighteenth century, these records have survived. They
were entered into a database in the past decade and published
online.101010https://bgb.huygens.knaw.nl Retrieved on: 13-05-2021 Similarly to
the voyages in DAS, voyages in BGB leave and arrive on a specific date to and
from a specific port. There are also overlapping voyages between the different
data sets. These overlapping voyages were identified by previous research. The
voyages transported several goods with a specific value, unit of measurement,
and specification. Voyages partook on ships that have different names over
time. Note that in comparison to DAS, less information about ships is captured
in this data set. BGB consists of 55 volumes of the 18th century containing
18,722 records.

  * •

_Places_ The Places data set is a single table containing information about
geographical places related to DAS and BGB. It contains toponym attestations
from DAS and BGB that are linked to a standardized place name (toponym
followed by a space and the ISO 3166-1 alpha 2 country code). GeoNames and
Wikidata URIs, as well as geographical coordinates are
available.111111https://www.geonames.org/; https://www.wikidata.org/ The
places dataset consists 4,252 instances.

###  4.2. Data and ontology modeling

In the data and ontology modeling phase, we aim to map data concepts to our
base ontology (CIDOC CRM) and add new attributes to that model if one of these
concepts does not fit the model.

We initially identify four main components: voyages, ships, persons, and
companies. Their relations include, but are not limited to, a voyage taking
place on a ship, a sailor participating on a voyage, and the VOC chamber
outfitting the ship. With these high-level concepts defined, we can specify
the lower-level concepts.

Voyages form the starting point of the lower-level, detailed concepts. A
voyage is the journey that a ship undertakes, moving from one place to
another. In this journey, sailors, soldiers, and other personnel are employed
to ensure that the ship’s cargo safely reaches its destination.

We aim to re-use CIDOC CRM classes and properties where possible. If no good
fit is available, we re-use classes and properties of Beretta et al. (Beretta
et al., 2019). Finally, we create new classes and properties accommodate the
envisioned data model. These newly created properties and classes are
integrated into the CIDOC CRM model as a subclass of existing classes, thus
making them interoperable with other CIDOC CRM extensions and models.

For reasons of readability, we do not display the entire ontology model, but
instead show a branch starting from a voyage in Figure 1. The complete mapping
figures can be found on the
GitHub.121212https://github.com/stijnschout3n/thewindinoursails

We initially defined the voyage instance as a subclass of CRM:E7 Activity. We
later changed this definition to a subclass of CRM:E9 Move to accommodate the
linking of goods (physical objects) to the voyage via the property CRM:P25
moved. We define ships as a subclass of CRM:E22 man-made object. The property
CRM:P8 took place on or within is used as property between the voyages and the
ships. (Incidentally, the property description of CRM:P8 took place on or
within, even uses this relation as an example: ‘It describes a period that can
be located with respect to the space defined by an CRM: E19 Physical Object
such as a ship or a building.’) We define the name of the ship through the
property CRM:P1 is identified by, which links the ship instance to a CRM:E41
Appellation instance. Other definitions and classifications, such as ship
types and productions, can be found in Figure 1.

###  4.3. Linked Data Transformation

In the transformation phase, we transform the acquired knowledge into RDF. In
short, a table structure is mapped to the different positions within a RDF
triple. We summarize our most important design decisions:

Time:

Time plays a central role in the maritime model. It relates to many concepts.
Such as the arrival or departure of ships. In the initial phases of the
conversion, the dates were sanitized and formatted according to the ISO
8601-1:2019 standard, which does not allow for uncertainties. The reasoning
behind this decision is as follows: more complex date formats are difficult to
calculate with and can therefore lead to lower usability of the artifact. In
consultation with the domain expert, we decided to add the original dates
(which sometimes state uncertainties) as an extra attribute for researchers to
consult while retaining data usability. Dates were formatted, according to CRM
standards, as URIs and not as XSD:Date literals. This decision helps link
events together that happened during the same time-span.

URIs:

We create URIs for unique (and unambiguous) entities. Identifying entities
across different data sets is rather tricky(Christen, 2012) and it is not
captured in this research project. Fortunately, other researchers have
performed data resolution and linking on these data sets (Hendriks et al.,
2020; Petram et al., ). For example, the BGB data set is linked at some
entries to DAS voyages. In these cases, the DAS voyage entity is reused,
effectively linking the two data sets. For goods shipped by the VOC, we create
unique URIs for every good that was carried on a specific voyage. This is to
allow for linking voyage-specific information (such as quantity or price) to
the good.

###  4.4. Data interaction

After converting and connecting the different data sets, we make our knowledge
graph accessible. We do this such that researchers are able to query the
created knowledge graph and wrangle the retrieved data to achieve their
research goal. Our knowledge graph can be accessed via
https://dutchshipsandsailors.nl/.

During the project, we also experimented with the Metaphacts
platform131313https://metaphacts.com/product which is well integrated with
GraphDB and comes in a pre-configured Docker container. Second, the Metaphacts
platform hosts very usable graph exploration tools. These tools can help in
answering the defined competency questions. Due to licensing issue, we cannot
expose our data through Metaphacts, but invite readers to try out this
platform, as well as exploring the data through the triplestore of their
choice.

Figure 2. Visual integration of the four data sets

###  4.5. Review of the initial design

Modeling data sets into an existing ontology is a delicate process, even with
high-level ontologies, such as CIDOC CRM. There is a continuous push and pull
between modeling from a data perspective and modeling from an ontology
perspective. In a perfect world, the ontology and the data are a match. In
this case, there exists a class and property for every envisioned concept.
Nevertheless, this is not the case when using a general ontology, as it
defeats its purpose. Therefore, the envisioned concepts should be fitted
within the scope of the ontology.

Differences in semantics between classes is not unique to this project. The
domain experts foresaw this and started developing an extension for the CRM
ontology. This extension, focusing on social, economic, and legal life, allows
for modeling concepts outside of the museum artifact
realm.141414http://ontome.net/project/64 Last visited: 16-09-2021

We now take a step back to a more meta-level by checking the research
questions’ progress. Any forthcoming discrepancies can be fixed in the second
design iteration.

Usability: The solution is approaching usable compared to the research goals.
In Subsection 4.4, some of the competency questions could be answered with the
interface. However, not all competency questions can be answered at this
moment. The first design iteration did not include the very extensive VOC
Opvarenden (VOCOP) data set. This data set contains rich records of the
sailors and soldiers who crewed the VOC ships. These records will be included
in the second design iteration, thus fulfilling the usability requirement.

Sustainability: This is not as strictly defined as the usability requirement.
We can only review the approach and the decisions made during the first design
iteration and determine where the sustainability can be improved. The
approach, thus far, is _transparent_ due to the extensive documentation. The
triple store allows for an _accessible_ environment where the knowledge graph
can be accessed interactively. By reviewing different approaches and methods,
the current approach seems reasonably _reusable_ as well. As we use a CRM
extension in the second iteration, we are actively promoting reusability.

##  5\. Second iteration

In the second design iteration, we process the results of the first review
period, and try to resolve any shortcomings regarding the research questions.
This iteration will be noticeably compacter due to our improved understanding
of the source data, making mere improvements rather than additions, and having
fewer data sets to convert.

The design thus far consists of voyages (DAS), cargo moved on voyages (BGB),
and the places where the voyages called (Places). However, several competency
questions are related to the people involved in these activities. In the
second design iteration, we add persons to the design to solve these
competency questions.

The VOC kept extensive records of its contracted crew members (‘Opvarenden’ in
Dutch). For the late seventeenth and eighteenth centuries, many pay ledgers
have been preserved. Pay ledger entries include personal data, ranks, wages,
and references to voyages. The voyages overlap with the voyages of the DAS
data set. The information of 774,200 crew members was preserved in the VOC
Opvarenden data set.

The core of the model does not change in the second iteration of the design.
The changes in the second iteration are more apparent in the details. One of
the more critical changes is related to the structure of the voyages. In
summary of the previous iteration: a voyage is linked to a starting point via
a departing activity and ends with the arrival of the voyage at its
destination. Persons and goods are directly linked to the voyage; we assume
that they were involved during the whole duration of the voyage. However, in
the VOCOP data set, this is not the case. Persons could change ships during a
stopover, effectively participating in two voyages. The first iteration could
not adequately convey this concept.

To correctly convey this concept, we introduce a new class. The leg class is,
similar to the voyage, a phase. It indicates a continuous part of the voyage
from and to a port of call. A voyage can consist of multiple legs. The time
between the end of the first leg and the beginning of the second leg is the
stopover. Finally, the sailors are linked to the departing and arriving
properties to indicate (dis)embarkment activities.

We choose the CIDOC CRM extension for economic, social, and legal life (SDHSS
prefix) for mapping the VOCOP data set. To relate persons to a voyage, the
_participation_ concept is used. Other mappings include ranks and their
descriptions.

The changes made during the second design iteration aim to improve the
_reusability_ of the artifact by utilizing an extended ontology. Creating
mapping figures and tables improves _transparency_ concerns. Converting the
VOCOP data set ensures that the _usability_ aspect as defined by the
competency questions is met. To improve the _maintainability_ aspect of the
artifact, the OntoME platform is proposed. The power of OntoME is in the reuse
of ontologies. The application has similarities with software development
platforms such as Github in that it allows users to collaborate and discuss
aspects concerning the creation and utilization of classes and properties.
Furthermore, other researchers can easily reuse the captured knowledge within
the ontology for their own projects.

Changes in the modeling have a cascading effect on the front-end part of the
artifact. In the second design iteration, these issues are examined and fixed.
Furthermore, the VOCOP data set was also integrated in the knowledge graph. An
example of how the different sources can be explored in a single screen
through the Metaphacts platform is shown in Figure 2.

Community interaction is not necessarily a research goal. However, it can have
clear benefits. As (Waagmeester et al., 2020) points out, sustainable growth
of the knowledge graph can be realized with community involvement. A proposal
that captures this issue is the utilization of the Geovistory
platform.151515https://docs.geovistory.com/ This platform can interact with a
knowledge graph and hosts a text editor that can recognize entities and
propose linkage to existing entities within the knowledge graph. Academics and
citizen-scientists could for example, with approval, enrich the knowledge
graph by browsing through the available data and adding connections within the
graph.

##  6\. Evaluation

With the completion of the second design iteration, most of the issues
discovered during the review period are resolved. These include, but are not
limited to, better utilization of the CIDOC CRM model and its extensions,
ontology management, community involvement, and the conversion of data on crew
member (VOCOP). In the second review, we explore these improvements and
related them to the primary research goal. The review consists of a workshop
with five domain experts from the Huygens Institute and the University of
Amsterdam. The expertise of the domain experts covered maritime history,
digital data management and data science, all with an affinity for (digital)
humanities research.

In the workshop, multiple discussions where brought up to challenge design
decisions made during the research project. These topics include, but are not
limited to, the utilization of CIDOC CRM, the conceptual mapping of a voyage,
and the proposed infrastructure. The feedback was overall positive but also
provided constructive criticism and included recommendations to put more focus
on data provenance in future iterations.

Initially, one of the goals of this project was to create a usable knowledge
graph in the field of Dutch historical maritime data. Discussion on what makes
a knowledge graph useful stems mainly from the perspective of maritime
historians, as they are the primary stakeholders. They provided the competency
questions which are used to test and improve the artifact. The answer to the
solvability of these questions can be found in Table 1.

However, the developed artifact is not only intended to be used by maritime
historians. Other interested parties can be drawn in by their research of an
event or a long-lost forefather who had a remarkable VOC career. Determining
if the knowledge graph is useful for these users is challenging as there is no
clear picture of who these people are and what their research questions
entail. Nevertheless, we found one subject and their research question is
explored in a short data
story.161616https://github.com/stijnschout3n/thewindinoursails/blob/main/Datastory.pdf
In short, the data story briefly dives into the life of a sailor, Hendrik
Prins. Quantitative and qualitative data analytics are performed in order to
portray the personal life and career of Hendrik Prins. The developed knowledge
graph makes the process of combining data sources straightforward and quick.

Question | Answer
---|---
Which VOC Chamber was accountable |  Can be answered by querying voyages that transported slaves
for the highest number of slaves transported?
How were the shipping routes in Asia divided between |  Can be answered. Already a preconfigured component in Metaphacts
the various VOC chambers (for example, did ships
from a specific chamber only sailed on certain routes)?
What was the average value of cargo on |  Can be answered. Goods, participants and prices are available
VOC return voyages per crew member/ship’s ton,
and how did this evolve over time?
To what extent was the value per ships’ ton on return |  Can be answered with data research, the components are available in the knowledge graph. Data paper provides an example.
voyages correlated with the skipper’s track record
(how much experience, i.e., how many previous voyages
/ how quickly did a ship get to Asia on an outbound voyage)?
How did the price of goods changed over time? | Can be answered, goods are linked
to voyages which have dates
Search for specific information in |  Can be answered with search components of Metaphacts
large data sets, such as names.
Can external factors, such as war, |  If the external factors are captured in the same ontology, then yes.
be linked to the intensity of shipping?
Table 1. Competency questions including statement regarding solvability

##  7\. Discussion and Conclusion

In this paper, we presented a reusable and maintainable knowledge graph for
Dutch maritime history. The knowledge graph was created using an iterative
design and development process that involved the conversion and linking of
four different data sets concerning the Dutch East India Company’s business
dealings. The resulting knowledge graph consists of more than 21 Million RDF
triples. The dataset is available through an online triple
store171717https://semanticweb.cs.vu.nl/eots/ and can be queried through a
SPARQL endpoint181818https://semanticweb.cs.vu.nl/eots/sparql. The data itself
is made available under open licences191919The Knowledge Graph and scripts to
produce it is published as four separate datasets through the following DOIs:
Places dataset: https://doi.org/10.5281/zenodo.5509814 ; VOCOPV dataset:
https://doi.org/10.5281/zenodo.5509820 DAS dataset:
https://doi.org/10.5281/zenodo.5507139 ; BGB data
https://doi.org/10.5281/zenodo.5507126. The data model and knowledge graph
were evaluated using competency questions provided specifically for this
project by a domain expert as well as competency questions formulated in
earlier work for this domain. Furthermore, the model and knowledge graph were
reviewed by five different experts with complementary knowledge regarding the
domain and computational analysis in the context of digital humanities
research.

The path for future work, appending more data sets and improving the model, is
paved by the used methodology. By following the framework of Jovanovik
(Jovanovik, 2016), which focused on reusable components, other researchers can
complete further iterations. Our documentation is made available to aid
researchers in understanding the decisions made during the design. Other data
sets could include the Northern
Musterrols202020https://www.noordelijkscheepvaartmuseum.nl/nl/collectie-en-
onderzoek/onderzoek/monsterrollen or the Archangel and Elbing
voyages212121http://resources.huygens.knaw.nl/archangel/app/inleiding in the
Baltic sea and similar data sets from other countries.

It should also be noted that even with many changes in the second design
iteration, the processing time of the second design iteration was much shorter
than the first iteration. Whereas the first iteration took around six weeks to
complete, we finished the second design iteration in two weeks. There is a
steep learning curve in understanding the concepts, tooling, and methods used.
Hopefully, other researchers will initiate new iterations and utilize the
results of this research project to face a less steep learning curve.

There are limitations to our project. In the second design iteration, the
proposed infrastructure is expanded with the use of OntoME and Geovistory.
However, the two extensions are not integrated entirely as they are web-based
and, therefore, difficult to integrate without direct support. In future
iterations, the proposed architecture, including extensions, could be
implemented in a production web server. However, this responsibility is
outside the scope of this research project and lies with the curators and
creators of the data.

During the workshop with domain experts, another limitation was discussed:
provenance. Data-altering improvements can also be captured in the artifact.
If these changes are stored in the knowledge graph, other researchers could
backtrack and understand the changes made in the data, thus improving the
transparency of the artifact. However, provenance information should be stored
in a separate data layer to avoid graph pollution.

The developed artifact, a knowledge graph, aims to support historians and
other users in their research. While the methodology also aims to develop the
artifact sustainably, such that the artifact can be used and expanded over the
years, we can conclude that the artifact can be deemed usable, as illustrated
by Table 1. Close to all competency questions could be answered by using the
knowledge graph. Due to the unavailability of data regarding events that had
an impact on a voyage but did not on the ship, the last competency question
(’Can external factors, such as war, be linked to the intensity of shipping?’)
could not be answered. Such research questions provide a useful insight into
the types of research questions that domain experts are interested in and a
solid direction for future extensions of the model and knowledge graph.

Furthermore, the participants of the workshop further confirmed the usability.
A sidestep has been made by evaluating the artifact with the data story.
During this evaluation, it was apparent that the artifact was useful for a
broader range of people even though only a single evaluation was performed.

The results regarding the sustainability of the methods used and the developed
artifact are multifaceted. There is no single answer that defines a
sustainable artifact and process. In the related project section, we discussed
the sustainability of the approaches taken by earlier project. We used these
insights to shape our project in a more sustainable way.

Developing a sustainable knowledge graph consists of making design decisions.
On a high or detailed level, design decisions influence the sustainability
aspects. Picking between these design options can often lead to dilemmas where
there is no clear answer. The choice between data usability and data
correctness is one of these dilemmas. Do the data need to be accurate, even if
this means the data will be difficult to process? In this project, we chose to
store both original and cleaned dates. We faced another dilemma during the
first review period. In this review, we examined the utilization of CIDOC CRM.
The intention to re-use existing concepts of an ontology can conflict with the
formal definition or use-case of that ontology. This in turn diminishes the
sustainability as different projects might use the concept differently.

The two design cycles make the project stand out as the development of the
artifact benefits immensely from receiving feedback. Design decisions made
during the cycles are dependent on the researchers involved. Personal beliefs,
experience, and other factors might steer the artifact in a different
direction. Therefore, the transparent and accessible nature of the project is
of utmost importance. It means that other researchers can not only continue
with further iterations but also backtrack decisions and take a different
path. Extended documentation, along with code an a Docker instance, can be
found on GitHub.222222https://github.com/stijnschout3n/thewindinoursails

###### Acknowledgements.

The authors want to thank Francesco Beretta, Vincent Alamercerey for help with
CIDOC CRM and participants in the workshop of KNAW HuC, the Huygens Institute
and the University of Amsterdam.

## References

  * (1)
  * Azzaoui et al. (2013) Kamal Azzaoui, Edgar Jacoby, Stefan Senger, Emiliano Cuadrado Rodríguez, Mabel Loza, Barbara Zdrazil, Marta Pinto, Antony J. Williams, Victor de la Torre, Jordi Mestres, Manuel Pastor, Olivier Taboureau, Matthias Rarey, Christine Chichester, Steve Pettifer, Niklas Blomberg, Lee Harland, Bryn Williams-Jones, and Gerhard F. Ecker. 2013.  Scientific competency questions as the basis for semantically enriched open pharmacological space development.  Drug Discovery Today 18, 17-18 (Sept. 2013), 843–852.  https://doi.org/10.1016/j.drudis.2013.05.008
  * Beretta et al. (2019) Francesco Beretta, Vincent Alamercery, Sebastiaan Derks, Lodewijk Petram, and Jonas Schneider. 2019\.  Geohistorical FAIR data: data integration and Interoperability using the OntoME platform.  Time Machine Conference 2019. (Oct. 2019).  https://halshs.archives-ouvertes.fr/halshs-02314003 Poster.
  * Berry (2012) David M. Berry (Ed.). 2012\.  Understanding digital humanities.  Palgrave Macmillan, Houndmills, Basingstoke, Hampshire ; New York.  OCLC: ocn701020028.
  * Blanke et al. (2012) Tobias Blanke, Gabriel Bodard, Michael Bryant, Stuart Dunn, Mark Hedges, Michael Jackson, and David Scott. 2012.  Linked data for humanities research—The SPQR experiment. In 2012 6th IEEE international conference on digital ecosystems and technologies (DEST). IEEE, 1–6.
  * Christen (2012) Peter Christen. 2012\.  Data matching: concepts and techniques for record linkage, entity resolution, and duplicate detection.  Springer, Berlin ; New York.  OCLC: ocn809643173.
  * Cidoc (2003) Crm Cidoc. 2003\.  The CIDOC conceptual reference model.  (2003).
  * Claire Warwick et al. (2012) Claire Warwick, Melissa M. Terras, and Julianne Nyhan (Eds.). 2012.  Digital humanities in practice: Digitisation and digital resources in the humanities.  Facet Publishing in association with UCL Centre for Digital Humanities, London.  OCLC: ocn707962666.
  * de Boer et al. (2015) Victor de Boer, Matthias van Rossum, and Rik Hoekstra. 2015\.  The Dutch Ships and Sailors Project.  DHCommons Journal 1 (July 2015).
  * Dimou (2020) A. Dimou. 2020.  Knowledge Graphs and Big Data Processing. Lecture Notes in Computer Science, Vol. 12072.  Springer International Publishing, Cham.  https://doi.org/10.1007/978-3-030-53199-7
  * Entjes (2015) Jeroen Entjes. 2015\.  Linking Maritime Datasets to Dutch Ships and Sailors Cloud - Case studies on Archangelvaart and Elbing.  Ph.D. Dissertation. VU Amsterdam.  http://www.victordeboer.com/wp-content/uploads/2015/08/jeroen_entjes_final_thesis.pdf
  * Gelderblom et al. (2013) Oscar Gelderblom, Abe de Jong, and Joost Jonker. 2013\.  The Formative Years of the Modern Corporation: The Dutch East India Company VOC, 1602–1623.  The Journal of Economic History 73, 4 (Dec. 2013), 1050–1076.  https://doi.org/10.1017/S0022050713000879
  * Haslhofer et al. (2018) Bernhard Haslhofer, Antoine Isaac, and Rainer Simon. 2018\.  Knowledge Graphs in the Libraries and Digital Humanities Domain.  arXiv:1803.03198 [cs] (2018), 1–8.  https://doi.org/10.1007/978-3-319-63962-8_291-1 arXiv: 1803.03198.
  * Hendriks et al. (2020) B Hendriks, P Groth, and M van Erp. 2020.  Recognising and linking entities in old dutchtext: a case study on voc notary records.  Collect & Connect Leiden (2020).
  * Hevner et al. (2004) Hevner, March, Park, and Ram. 2004\.  Design Science in Information Systems Research.  MIS Quarterly 28, 1 (2004), 75.  https://doi.org/10.2307/25148625
  * Hughes (2012) Lorna M. Hughes (Ed.). 2012\.  Evaluating and measuring the value, use and impact of digital collections.  Facet Publ, London.  OCLC: 774166752.
  * Hyvönen et al. (2020) Eero Hyvönen et al. 2020\.  Linked open data infrastructure for digital humanities in finland.  Proceedings of Digital Humanities in Nordic Countries (DHN 2020) (2020).
  * Jovanovik (2016) Milos Jovanovik. 2016\.  Linked Data Application Development Methodology.  Ph.D. Dissertation.
  * (19) Lodewijk Petram, Marijn Koolen, Rutger van Koert, Melvin Wevers, and Jelle van Lottum.  Data on the Maritime Workforce of the Dutch East India Company in the 18th Century.  to be published (????).
  * Waagmeester et al. (2020) Andra Waagmeester, Gregory Stupp, Sebastian Burgstaller-Muehlbacher, Benjamin M Good, Malachi Griffith, Obi L Griffith, Kristina Hanspers, Henning Hermjakob, Toby S Hudson, Kevin Hybiske, Sarah M Keating, Magnus Manske, Michael Mayers, Daniel Mietchen, Elvira Mitraka, Alexander R Pico, Timothy Putman, Anders Riutta, Nuria Queralt-Rosinach, Lynn M Schriml, Thomas Shafee, Denise Slenter, Ralf Stephan, Katherine Thornton, Ginger Tsueng, Roger Tu, Sabah Ul-Hasan, Egon Willighagen, Chunlei Wu, and Andrew I Su. 2020.  Wikidata as a knowledge graph for the life sciences.  eLife 9 (March 2020), e52614.  https://doi.org/10.7554/eLife.52614
  * Wieringa (2014) Roel J Wieringa. 2014\.  Design science methodology for information systems and software engineering.  Springer.
  * Wilkinson et al. (2016) Mark D. Wilkinson, Michel Dumontier, IJsbrand Jan Aalbersberg, Gabrielle Appleton, Myles Axton, Arie Baak, Niklas Blomberg, Jan-Willem Boiten, Luiz Bonino da Silva Santos, Philip E. Bourne, Jildau Bouwman, Anthony J. Brookes, Tim Clark, Mercè Crosas, Ingrid Dillo, Olivier Dumon, Scott Edmunds, Chris T. Evelo, Richard Finkers, Alejandra Gonzalez-Beltran, Alasdair J.G. Gray, Paul Groth, Carole Goble, Jeffrey S. Grethe, Jaap Heringa, Peter A.C ’t Hoen, Rob Hooft, Tobias Kuhn, Ruben Kok, Joost Kok, Scott J. Lusher, Maryann E. Martone, Albert Mons, Abel L. Packer, Bengt Persson, Philippe Rocca-Serra, Marco Roos, Rene van Schaik, Susanna-Assunta Sansone, Erik Schultes, Thierry Sengstag, Ted Slater, George Strawn, Morris A. Swertz, Mark Thompson, Johan van der Lei, Erik van Mulligen, Jan Velterop, Andra Waagmeester, Peter Wittenburg, Katherine Wolstencroft, Jun Zhao, and Barend Mons. 2016\.  The FAIR Guiding Principles for scientific data management and stewardship.  Scientific Data 3, 1 (Dec. 2016), 160018\.  https://doi.org/10.1038/sdata.2016.18