Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • 26-refactoring-the-infrastructure
  • MPGCQM_workshop
  • app
  • backup-05-08
  • develop
  • kappa_sigma_notebook
  • master
  • staging
8 results

Target

Select target project
  • nomad-lab/analytics
1 result
Select Git revision
  • 26-refactoring-the-infrastructure
  • MPGCQM_workshop
  • app
  • backup-05-08
  • develop
  • kappa_sigma_notebook
  • master
  • staging
8 results
Show changes
View open merge request
Hide whitespace changes
Inline Side-by-side
Showing
with 2 additions and 12212 deletions
3rdparty/gap/doc/gap_fit.rst deleted 100644 → 0
View file @ 23bcf97f
.. HQ XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
.. HQ X
.. HQ X quippy: Python interface to QUIP atomistic simulation library
.. HQ X
.. HQ X Copyright James Kermode 2010
.. HQ X
.. HQ X These portions of the source code are released under the GNU General
.. HQ X Public License, version 2, http://www.gnu.org/copyleft/gpl.html
.. HQ X
.. HQ X If you would like to license the source code under different terms,
.. HQ X please contact James Kermode, james.kermode@gmail.com
.. HQ X
.. HQ X When using this software, please cite the following reference:
.. HQ X
.. HQ X http://www.jrkermode.co.uk/quippy
.. HQ X
.. HQ XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
The gap_fit program
=====================================================
Main options
------------
.. autofunction:: quippy.gap_fit_parse_command_line
GAP options
-----------
.. autofunction:: quippy.gap_fit_parse_gap_str
`sparse_method` options are:
- RANDOM: default, chooses n_sparse random datapoints
- PIVOT: based on the full covariance matrix finds the n_sparse "pivoting" points
- CLUSTER: based on the full covariance matrix performs a k-medoid clustering into n_sparse clusters, returning the medoids
- UNIFORM: makes a histogram of the data based on n_sparse and returns a data point from each bin
- KMEANS: k-means clustering based on the data points
- COVARIANCE: greedy data point selection based on the sparse covariance matrix, to minimise the GP variance of all datapoints
- UNIQ: selects unique datapoints from the dataset
- FUZZY: fuzzy k-means clustering
- FILE: reads sparse points from a file
- INDEX_FILE: reads indices of sparse points from a file
- CUR_COVARIANCE: CUR, based on the full covariance matrix
- CUR_POINTS: CUR, based on the datapoints
3rdparty/gap/doc/gap_fitting_tutorial.ipynb deleted 100644 → 0
View file @ 23bcf97f
Source diff could not be displayed: it is too large. Options to address this: view the blob.
3rdparty/gap/doc/index.rst deleted 100644 → 0
View file @ 23bcf97f
.. HQ XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
.. HQ X
.. HQ X GAP: Gaussian Approximation Potential
.. HQ X
.. HQ X Copyright Albert Bartok-Partay, Gabor Csanyi 2010-2019
.. HQ X
.. HQ X gc121@cam.ac.uk
.. HQ X
.. HQ X
.. HQ XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
.. gap documentation master file
GAP and SOAP documentation
=============================
.. module:: gap
These are the documentation pages for Gaussian Approximation Potential
(GAP) code. GAP is a plugin for QUIP, which itself can be used as a
plugin to LAMMPS or called from ASE. For installation, see the QUIP
documentation pages. The information below is on the usage of GAP.
The purpose of the GAP code is to fit interatomic potentials and then
use them for molecular simulation.
Contents
========
.. toctree::
:maxdepth: 1
gap_fit.rst
tutorials.rst
Indices and tables
==================
* :ref:`genindex`
* :ref:`modindex`
* :ref:`search`
* `Module Source Code <_modules/index.html>`_
3rdparty/gap/doc/modcontents.py deleted 100644 → 0
View file @ 23bcf97f
import sphinx
if sphinx.__version__ < '1.0.1':
raise RuntimeError("Sphinx 1.0.1 or newer is required")
import inspect
import pydoc
def process_docstring(app, what, name, obj, options, lines):
if what == 'module':
print 'Adding contents listing for module %s' % name
title = 'Module contents for :mod:`%s`:' % name
lines.append(title)
classes = module_classes(obj)
if classes:
lines.append("""
.. rubric:: Classes
.. autosummary::
%s
""" % '\n'.join([' %s' % cls for cls in classes]))
funcs = module_functions(obj)
if funcs:
lines.append("""
.. rubric:: Functions
.. autosummary::
%s
""" % '\n'.join([' %s' % func for func in funcs]))
# FIXME add documentation for module attributes
attributes = module_attributes(obj)
if attributes:
lines.append("""
.. rubric:: Attributes
%s
""" % attributes_table(obj, attributes))
def attributes_table(mod, attributes):
MAX_LEN=100
lines = ['']
values = [ str(getattr(mod, attr)) for attr in attributes ]
for i, value in enumerate(values):
if len(value) > MAX_LEN:
values[i] = '---'
attributes = [':attr:`%s`' % attr for attr in attributes ]
col1 = max(len(attr) for attr in attributes)
col2 = max(len(value) for value in values)
fmt = '%%-%ds %%-%ds' % (col1, col2)
head = fmt % ('='*col1, '='*col2)
lines.append(head)
lines.append(fmt % ('Name', 'Value'))
lines.append(head)
for attr, value in zip(attributes, values):
lines.append(fmt % (attr, value))
lines.append(head)
lines.append('')
return '\n'.join(lines)
def module_functions(mod):
if hasattr(mod, '__alldoc__'):
allsymbols = mod.__alldoc__
elif hasattr(mod, '__all__'):
allsymbols = mod.__all__
else:
allsymbols = [name for (name, obj) in inspect.getmembers(mod)]
return [name for name in allsymbols if
not name.startswith('_') and
inspect.isfunction(getattr(mod, name)) and
pydoc.getdoc(getattr(mod, name))]
def module_classes(mod):
if hasattr(mod, '__alldoc__'):
allsymbols = mod.__alldoc__
elif hasattr(mod, '__all__'):
allsymbols = mod.__all__
else:
allsymbols = [name for (name, obj) in inspect.getmembers(mod)]
return [name for name in allsymbols if
not name.startswith('_') and
inspect.isclass(getattr(mod, name)) and
pydoc.getdoc(getattr(mod, name))]
def module_attributes(mod):
if hasattr(mod, '__alldoc__'):
allsymbols = mod.__alldoc__
elif hasattr(mod, '__all__'):
allsymbols = mod.__all__
else:
allsymbols = [name for (name, obj) in inspect.getmembers(mod)]
return [name for name in allsymbols if
not name.startswith('_') and
not callable(getattr(mod, name)) and
pydoc.getdoc(getattr(mod, name))]
def setup(app):
print('modcontents extension loading')
app.connect('autodoc-process-docstring', process_docstring)
3rdparty/gap/doc/quippy-descriptor-tutorial.ipynb deleted 100644 → 0
View file @ 23bcf97f
%% Cell type:markdown id: tags:
# Computing descriptors from python
This tutorial shows how to use quippy to calculate various descriptors for a given atomic structure, contained in an Atoms object.
%% Cell type:code id: tags:
``` python
import quippy
```
%% Cell type:code id: tags:
``` python
from quippy import descriptors
```
%% Cell type:markdown id: tags:
Create a configuration:
%% Cell type:code id: tags:
``` python
at = quippy.diamond(3.5, 6)
at.positions
```
%% Output
array([[ 0. , 0. , 0. ],
[ 0.875, 0.875, 0.875],
[ 1.75 , 1.75 , 0. ],
[ 2.625, 2.625, 0.875],
[ 1.75 , 0. , 1.75 ],
[ 2.625, 0.875, 2.625],
[ 0. , 1.75 , 1.75 ],
[ 0.875, 2.625, 2.625]])
%% Cell type:code id: tags:
``` python
at.Z
```
%% Output
FortranArray([6, 6, 6, 6, 6, 6, 6, 6], dtype=int32)
%% Cell type:markdown id: tags:
All descriptors are instantiated by a call to `Descriptor()`, which takes the descriptor initialisation string as its only argument. For a list of available descriptors and their parameters, see the following list, auto-generated using the following command: ```quip descriptor_str="--help"```
bispectrum_so4
bispectrum_so3
behler
distance_2b
coordination
angle_3b
co_angle_3b
co_distance_2b
cosnx
trihis
water_monomer
water_dimer
A2_dimer
AB_dimer
bond_real_space
atom_real_space
power_so3
power_so4
soap
AN_monomer
general_monomer
general_dimer
general_trimer
rdf
as_distance_2b
molecule_lo_d
alex
com_dimer
distance_Nb
%% Cell type:markdown id: tags:
## A simple descriptor: pairwise distances
%% Cell type:markdown id: tags:
Here we use a simple pair distance between carbon atoms, with a cutoff of 1.5 Angstrom. There are several descriptors that can do this, one is ``distance_2b``, which takes a cutoff argument, and two Z values to specify the atom types. Alternatively, the ``distance_Nb`` descriptor could also do this, with ``order=2``, and it takes a string of Zs to specify the atom types. This is more general, ``order=3`` is a triangle-like three-body descriptor of the three sides of a triangle of 3 atoms.
%% Cell type:code id: tags:
``` python
desc = descriptors.Descriptor("distance_2b Z1=6 Z2=6 cutoff=4")
```
%% Cell type:markdown id: tags:
The descriptor dimension is the length of the descriptor vector. For the scalar distance this is 1.
%% Cell type:code id: tags:
``` python
desc.n_dim # number of dimensions
```
%% Output
1
%% Cell type:markdown id: tags:
This descriptor is very simple: it is scalar (dimension 1), and hence only has a single permutation.
%% Cell type:code id: tags:
``` python
desc.n_perm # number of permutations
```
%% Output
1
%% Cell type:code id: tags:
``` python
desc.permutations() # array of permutation arrays
```
%% Output
array([[1]], dtype=int32)
%% Cell type:markdown id: tags:
Many descriptors rely on the neighbour connectivity, so we need to call `calc_connect`, after setting the Atoms cutoff with a skin of 1 A:
%% Cell type:code id: tags:
``` python
at.set_cutoff(desc.cutoff()+1.0)
at.calc_connect()
```
%% Cell type:markdown id: tags:
We can now calculate how many instances of this descriptor are found in an `Atoms` (or `ASEAtoms`) object:
%% Cell type:code id: tags:
``` python
desc.count(at)
```
%% Output
368
%% Cell type:markdown id: tags:
This also works transparently for iterables (such as lists), returning a list of the counts.
%% Cell type:markdown id: tags:
We can also calculate the actual descriptor values -- in this case, the list of pairwise distances:
%% Cell type:code id: tags:
``` python
d = desc.calc(at)
d
```
%% Output
{'cutoff': array([ 0.30420743, 0.30420743, 1. , 1. , 1. ,
0.30420743, 0.30420743, 0.30420743, 1. , 1. ,
1. , 0.30420743, 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 0.30420743, 1. , 1. ,
1. , 1. , 1. , 0.30420743, 1. ,
0.30420743, 0.30420743, 1. , 1. , 1. ,
1. , 1. , 0.30420743, 1. , 0.30420743,
1. , 0.30420743, 0.30420743, 1. , 0.30420743,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 0.30420743, 1. , 0.30420743,
1. , 0.30420743, 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 0.30420743, 0.30420743, 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 0.30420743, 1. , 1. , 0.30420743,
1. , 1. , 0.30420743, 0.30420743, 1. ,
1. , 1. , 0.30420743, 1. , 1. ,
1. , 0.30420743, 0.30420743, 1. , 1. ,
1. , 1. , 1. , 0.30420743, 0.30420743,
1. , 1. , 1. , 1. , 0.30420743,
0.30420743, 1. , 1. , 1. , 1. ,
1. , 1. , 0.30420743, 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 0.30420743, 0.30420743, 1. ,
1. , 1. , 1. , 0.30420743, 1. ,
0.30420743, 1. , 1. , 0.30420743, 0.30420743,
1. , 0.30420743, 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
0.30420743, 0.30420743, 1. , 1. , 1. ,
0.30420743, 1. , 1. , 0.30420743, 1. ,
1. , 0.30420743, 1. , 1. , 1. ,
1. , 1. , 0.30420743, 1. , 0.30420743,
1. , 1. , 1. , 0.30420743, 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 0.30420743, 1. , 0.30420743,
1. , 0.30420743, 1. , 0.30420743, 1. ,
1. , 1. , 1. , 1. , 1. ,
0.30420743, 1. , 0.30420743, 1. , 1. ,
0.30420743, 1. , 1. , 1. , 1. ,
1. , 0.30420743, 1. , 0.30420743, 0.30420743,
1. , 1. , 0.30420743, 0.30420743, 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 0.30420743, 1. , 1. , 1. ,
0.30420743, 0.30420743, 1. , 0.30420743, 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 0.30420743, 1. ,
0.30420743, 1. , 0.30420743, 1. , 1. ,
1. , 1. , 1. , 0.30420743, 1. ,
1. , 0.30420743, 1. , 0.30420743, 1. ,
1. , 1. , 0.30420743, 1. , 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 0.30420743, 1. , 0.30420743,
1. , 0.30420743, 1. , 1. , 1. ,
0.30420743, 0.30420743, 1. , 1. , 1. ,
1. , 1. , 0.30420743, 1. , 0.30420743,
1. , 1. , 0.30420743, 1. , 1. ,
1. , 1. , 1. , 0.30420743, 1. ,
0.30420743, 1. , 0.30420743, 1. , 1. ,
1. , 1. , 1. , 1. , 0.30420743,
0.30420743, 1. , 1. , 0.30420743, 1. ,
1. , 1. , 1. , 1. , 1. ,
1. , 1. , 0.30420743, 0.30420743, 1. ,
0.30420743, 1. , 1. , 1. , 1. ,
1. , 1. , 1. , 1. , 0.30420743,
0.30420743, 1. , 1. , 1. , 0.30420743,
1. , 1. , 1. , 1. , 1. ,
0.30420743, 1. , 1. , 0.30420743, 1. ,
0.30420743, 1. , 1. , 1. , 1. ,
1. , 1. , 0.30420743, 1. , 0.30420743,
1. , 0.30420743, 1. , 1. , 1. ,
1. , 1. , 1. ]),
'descriptor': array([[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 3.5 ],
[ 1.51554446],
[ 2.47487373],
[ 1.51554446],
[ 2.47487373],
[ 1.51554446],
[ 2.47487373],
[ 1.51554446],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 3.5 ],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 3.5 ],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 3.5 ],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.47487373],
[ 1.51554446],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 2.47487373],
[ 1.51554446],
[ 3.81403658],
[ 2.47487373],
[ 1.51554446],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 2.90204669],
[ 2.90204669],
[ 3.81403658],
[ 3.81403658],
[ 1.51554446],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 1.51554446],
[ 2.47487373],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 1.51554446],
[ 3.5 ],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 2.90204669],
[ 3.81403658],
[ 3.5 ],
[ 2.47487373],
[ 3.81403658],
[ 3.5 ],
[ 2.47487373],
[ 3.81403658],
[ 3.5 ],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 3.5 ],
[ 3.81403658],
[ 3.5 ],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 1.51554446],
[ 2.47487373],
[ 1.51554446],
[ 3.5 ],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 1.51554446],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 2.90204669],
[ 2.47487373],
[ 1.51554446],
[ 3.5 ],
[ 2.90204669],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.5 ],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 3.5 ],
[ 3.5 ],
[ 2.47487373],
[ 3.81403658],
[ 3.5 ],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 1.51554446],
[ 3.5 ],
[ 3.81403658],
[ 3.5 ],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 2.90204669],
[ 3.5 ],
[ 1.51554446],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 2.90204669],
[ 3.81403658],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ]])}
%% Cell type:markdown id: tags:
Notice that the first array is called ``cutoff``, and it supplies the value of a cutoff function, implicit in all descriptors, which takes the descriptor value to zero as the atoms approach the cutoff, i.e. in this case as the distance between the two atoms approaches the cutoff. It is more complicated for three-body and higher-body descriptors, but the end result is always a descriptor which changes smoothly with atomic positions.
%% Cell type:markdown id: tags:
Here is a histogram of the resulting descriptor array, i.e. of the interatomic distances
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
plt.hist(d.descriptor)
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
Calculate size of descriptor data:
%% Cell type:code id: tags:
``` python
n_desc, n_cross = desc.descriptor_sizes(at)
print("n_desc=%d n_cross=%d" % (n_desc,n_cross))
```
%% Output
n_desc=368 n_cross=736
%% Cell type:markdown id: tags:
`n_desc` is number of descriptors, `n_cross` is number of gradients
%% Cell type:markdown id: tags:
Gradients are returned in a `DescriptorCalcResult` object ( if the ``grad=True`` option is set ) which has the following elements:
- `descriptor` = the descriptor values (equivalent to `desc.calc(at)`)
- `grad` = the gradients
- `index` = the indices (`i_desc`, `i_atom`, `i_coord`)
%% Cell type:code id: tags:
``` python
res = desc.calc(at, grad=True)
list(res)
```
%% Output
['descriptor', 'cutoff_grad', 'grad', 'index', 'cutoff']
%% Cell type:code id: tags:
``` python
res.descriptor[:,:10]
```
%% Output
array([[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 3.5 ],
[ 1.51554446],
[ 2.47487373],
[ 1.51554446],
[ 2.47487373],
[ 1.51554446],
[ 2.47487373],
[ 1.51554446],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 3.5 ],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 3.5 ],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 3.5 ],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.47487373],
[ 1.51554446],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 2.47487373],
[ 1.51554446],
[ 3.81403658],
[ 2.47487373],
[ 1.51554446],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 2.90204669],
[ 2.90204669],
[ 3.81403658],
[ 3.81403658],
[ 1.51554446],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 1.51554446],
[ 2.47487373],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 1.51554446],
[ 3.5 ],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 2.90204669],
[ 3.81403658],
[ 3.5 ],
[ 2.47487373],
[ 3.81403658],
[ 3.5 ],
[ 2.47487373],
[ 3.81403658],
[ 3.5 ],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 3.5 ],
[ 3.81403658],
[ 3.5 ],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 1.51554446],
[ 2.47487373],
[ 1.51554446],
[ 3.5 ],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 1.51554446],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 2.90204669],
[ 2.47487373],
[ 1.51554446],
[ 3.5 ],
[ 2.90204669],
[ 2.47487373],
[ 3.81403658],
[ 2.47487373],
[ 3.5 ],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 3.5 ],
[ 3.5 ],
[ 2.47487373],
[ 3.81403658],
[ 3.5 ],
[ 2.47487373],
[ 2.90204669],
[ 3.5 ],
[ 2.47487373],
[ 2.90204669],
[ 2.47487373],
[ 1.51554446],
[ 3.5 ],
[ 3.81403658],
[ 3.5 ],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 3.81403658],
[ 2.90204669],
[ 3.5 ],
[ 2.90204669],
[ 3.5 ],
[ 1.51554446],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ],
[ 3.81403658],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 2.90204669],
[ 1.51554446],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 2.90204669],
[ 3.81403658],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 3.81403658],
[ 1.51554446],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.47487373],
[ 2.90204669],
[ 1.51554446],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 2.90204669],
[ 3.81403658],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ],
[ 3.5 ]])
%% Cell type:code id: tags:
``` python
res.grad[:,:10]
```
%% Output
array([[[-0.22941573],
[ 0.6882472 ],
[ 0.6882472 ]],
[[ 0.22941573],
[-0.6882472 ],
[-0.6882472 ]],
[[ 0.6882472 ],
[-0.22941573],
[ 0.6882472 ]],
...,
[[ 0. ],
[ 1. ],
[ 0. ]],
[[-0. ],
[-0. ],
[-1. ]],
[[ 0. ],
[ 0. ],
[ 1. ]]])
%% Cell type:code id: tags:
``` python
res.index[:,:10]
```
%% Output
array([[ 1, 1],
[ 1, 2],
[ 2, 1],
...,
[367, 8],
[368, 8],
[368, 8]], dtype=int32)
%% Cell type:markdown id: tags:
## A many-body descriptor: SOAP
%% Cell type:code id: tags:
``` python
desc = descriptors.Descriptor("soap cutoff=3 l_max=4 n_max=4 atom_sigma=0.5 n_Z=1 Z={6} ")
```
%% Cell type:code id: tags:
``` python
at.set_cutoff(desc.cutoff())
at.calc_connect()
```
%% Cell type:markdown id: tags:
There are now only 8 descriptors, because SOAP produces one for each atom in the structure
%% Cell type:code id: tags:
``` python
desc.descriptor_sizes(at)
```
%% Output
(8, 258)
%% Cell type:markdown id: tags:
But each descriptor now is a long vector, because it encodes the entire environment of the atom up to the cutoff. The length of the vector depends on ```l_max``` and ```n_max``` and also on the number of atom types.
%% Cell type:code id: tags:
``` python
desc.n_dim
```
%% Output
51
%% Cell type:code id: tags:
``` python
desc.calc(at)
```
%% Output
{'cutoff': array([ 1., 1., 1., 1., 1., 1., 1., 1.]),
'descriptor': array([[ 4.73572570e-01, 5.11159766e-06, 1.01308559e-06,
1.60812107e-03, 5.47553526e-04, 6.27690415e-01,
4.95953867e-05, 9.82948612e-06, 1.54957207e-02,
5.47660232e-03, 4.15981924e-01, 2.40600155e-04,
4.76854071e-05, 7.46577370e-02, 2.73883482e-02,
3.18341025e-01, -5.92864850e-07, -1.17486739e-07,
-1.41512045e-03, 1.81695390e-03, 2.98357233e-01,
-4.06747903e-06, -8.06043591e-07, -9.64209426e-03,
1.28503115e-02, 1.06996282e-01, 3.43814942e-08,
6.81242239e-09, 6.22655008e-04, 3.01463911e-03,
-2.81631715e-03, 5.36228874e-09, 1.06263310e-09,
1.15339605e-05, -1.66628637e-05, -2.63952342e-03,
3.67891552e-08, 7.29042791e-09, 7.85880344e-05,
-1.17847232e-04, -1.33866723e-03, -4.39778747e-10,
-8.71386960e-11, -7.17707326e-06, -3.90981365e-05,
8.37426278e-06, 2.81263730e-12, 5.57301933e-13,
4.13635000e-08, 2.53540179e-07, 0.00000000e+00],
[ 2.49140590e-01, 3.06897443e-04, 2.24691086e-04,
2.63645993e-03, 8.44257436e-04, 5.55127085e-01,
2.97737414e-03, 2.17965866e-03, 2.54342326e-02,
8.40520121e-03, 6.18458197e-01, 1.44425393e-02,
1.05720969e-02, 1.22683598e-01, 4.18410093e-02,
2.47745044e-01, -3.90932807e-05, -3.07806268e-05,
-1.98138033e-03, 2.35398492e-03, 3.90335370e-01,
-2.68180541e-04, -2.11137369e-04, -1.35142491e-02,
1.65884021e-02, 1.23178658e-01, 2.48989464e-06,
2.10833239e-06, 7.59273692e-04, 3.41878216e-03,
-2.19112484e-03, 3.46295943e-07, 2.69352586e-07,
1.61901122e-05, -2.16062298e-05, -3.45223263e-03,
2.37559580e-06, 1.84760358e-06, 1.10427087e-04,
-1.52258473e-04, -1.54068059e-03, -3.11918586e-08,
-2.60914232e-08, -8.76903344e-06, -4.43838310e-05,
9.63517840e-06, 1.95376148e-10, 1.61445692e-10,
5.06386375e-08, 2.88103389e-07, 0.00000000e+00],
[ 2.23838398e-01, 3.78089923e-05, 4.06539273e-05,
1.92860775e-03, 7.66122563e-04, 5.09433866e-01,
3.53632752e-04, 4.29949407e-04, 1.85669919e-02,
7.76143527e-03, 5.79710332e-01, 1.65484402e-03,
2.28688255e-03, 8.93771594e-02, 3.93431430e-02,
2.92339646e-01, -1.55904376e-04, 4.02508059e-04,
-1.89111550e-03, 3.67450327e-03, 4.70463509e-01,
-1.01397193e-03, 3.22651119e-03, -1.28156879e-02,
2.67826166e-02, 1.90902163e-01, 4.60348034e-04,
3.74809897e-03, 1.39702795e-03, 1.25434798e-02,
-2.58972445e-03, 1.32417523e-06, -3.56120051e-06,
1.55903764e-05, -3.32865245e-05, -4.16765521e-03,
8.60690247e-06, -2.85979893e-05, 1.05611567e-04,
-2.42520968e-04, -2.39161731e-03, -5.59006314e-06,
-4.74862442e-05, -1.67590764e-05, -1.59583514e-04,
1.49810595e-05, 3.39496223e-08, 3.00874317e-07,
1.00726608e-07, 1.01530792e-06, 0.00000000e+00],
[ 2.21912827e-01, 3.22794357e-04, 2.74608154e-04,
2.71347565e-03, 9.21261052e-04, 5.22685782e-01,
3.10761607e-03, 2.61564838e-03, 2.61192465e-02,
9.15687879e-03, 6.15557988e-01, 1.49623108e-02,
1.24681640e-02, 1.25730679e-01, 4.55125688e-02,
2.60247510e-01, -3.16196242e-04, -5.90933824e-04,
-2.70412292e-03, 2.39722324e-03, 4.33440876e-01,
-2.09668723e-03, -3.79982381e-03, -1.80490450e-02,
1.69318863e-02, 1.52602189e-01, 6.07534416e-04,
2.09763755e-03, 4.23793685e-03, 3.79617175e-03,
-2.29869885e-03, 2.76019211e-06, 5.15137054e-06,
2.24848968e-05, -2.19757551e-05, -3.82847099e-03,
1.83031786e-05, 3.31247321e-05, 1.49959011e-04,
-1.55178180e-04, -1.90621225e-03, -7.49516446e-06,
-2.58557011e-05, -5.12063357e-05, -4.87638925e-05,
1.19056128e-05, 4.62340053e-08, 1.59350048e-07,
3.09463332e-07, 3.13393204e-07, 0.00000000e+00],
[ 2.23838398e-01, 3.78089923e-05, 4.06539273e-05,
1.92860775e-03, 7.66122563e-04, 5.09433866e-01,
3.53632752e-04, 4.29949407e-04, 1.85669919e-02,
7.76143527e-03, 5.79710332e-01, 1.65484402e-03,
2.28688255e-03, 8.93771594e-02, 3.93431430e-02,
2.92339646e-01, -1.55904376e-04, 4.02508059e-04,
-1.89111550e-03, 3.67450327e-03, 4.70463509e-01,
-1.01397193e-03, 3.22651119e-03, -1.28156879e-02,
2.67826166e-02, 1.90902163e-01, 4.60348034e-04,
3.74809897e-03, 1.39702795e-03, 1.25434798e-02,
-2.58972445e-03, 1.32417523e-06, -3.56120051e-06,
1.55903764e-05, -3.32865245e-05, -4.16765521e-03,
8.60690247e-06, -2.85979893e-05, 1.05611567e-04,
-2.42520968e-04, -2.39161731e-03, -5.59006314e-06,
-4.74862442e-05, -1.67590764e-05, -1.59583514e-04,
1.49810595e-05, 3.39496223e-08, 3.00874317e-07,
1.00726608e-07, 1.01530792e-06, 0.00000000e+00],
[ 2.21912827e-01, 3.22794357e-04, 2.74608154e-04,
2.71347565e-03, 9.21261052e-04, 5.22685782e-01,
3.10761607e-03, 2.61564838e-03, 2.61192465e-02,
9.15687879e-03, 6.15557988e-01, 1.49623108e-02,
1.24681640e-02, 1.25730679e-01, 4.55125688e-02,
2.60247510e-01, -3.16196242e-04, -5.90933824e-04,
-2.70412292e-03, 2.39722324e-03, 4.33440876e-01,
-2.09668723e-03, -3.79982381e-03, -1.80490450e-02,
1.69318863e-02, 1.52602189e-01, 6.07534416e-04,
2.09763755e-03, 4.23793685e-03, 3.79617175e-03,
-2.29869885e-03, 2.76019211e-06, 5.15137054e-06,
2.24848968e-05, -2.19757551e-05, -3.82847099e-03,
1.83031786e-05, 3.31247321e-05, 1.49959011e-04,
-1.55178180e-04, -1.90621225e-03, -7.49516446e-06,
-2.58557011e-05, -5.12063357e-05, -4.87638925e-05,
1.19056128e-05, 4.62340053e-08, 1.59350048e-07,
3.09463332e-07, 3.13393204e-07, 0.00000000e+00],
[ 2.23838398e-01, 3.78089923e-05, 4.06539273e-05,
1.92860775e-03, 7.66122563e-04, 5.09433866e-01,
3.53632752e-04, 4.29949407e-04, 1.85669919e-02,
7.76143527e-03, 5.79710332e-01, 1.65484402e-03,
2.28688255e-03, 8.93771594e-02, 3.93431430e-02,
2.92339646e-01, -1.55904376e-04, 4.02508059e-04,
-1.89111550e-03, 3.67450327e-03, 4.70463509e-01,
-1.01397193e-03, 3.22651119e-03, -1.28156879e-02,
2.67826166e-02, 1.90902163e-01, 4.60348034e-04,
3.74809897e-03, 1.39702795e-03, 1.25434798e-02,
-2.58972445e-03, 1.32417523e-06, -3.56120051e-06,
1.55903764e-05, -3.32865245e-05, -4.16765521e-03,
8.60690247e-06, -2.85979893e-05, 1.05611567e-04,
-2.42520968e-04, -2.39161731e-03, -5.59006314e-06,
-4.74862442e-05, -1.67590764e-05, -1.59583514e-04,
1.49810595e-05, 3.39496223e-08, 3.00874317e-07,
1.00726608e-07, 1.01530792e-06, 0.00000000e+00],
[ 2.21912827e-01, 3.22794357e-04, 2.74608154e-04,
2.71347565e-03, 9.21261052e-04, 5.22685782e-01,
3.10761607e-03, 2.61564838e-03, 2.61192465e-02,
9.15687879e-03, 6.15557988e-01, 1.49623108e-02,
1.24681640e-02, 1.25730679e-01, 4.55125688e-02,
2.60247510e-01, -3.16196242e-04, -5.90933824e-04,
-2.70412292e-03, 2.39722324e-03, 4.33440876e-01,
-2.09668723e-03, -3.79982381e-03, -1.80490450e-02,
1.69318863e-02, 1.52602189e-01, 6.07534416e-04,
2.09763755e-03, 4.23793685e-03, 3.79617175e-03,
-2.29869885e-03, 2.76019211e-06, 5.15137054e-06,
2.24848968e-05, -2.19757551e-05, -3.82847099e-03,
1.83031786e-05, 3.31247321e-05, 1.49959011e-04,
-1.55178180e-04, -1.90621225e-03, -7.49516446e-06,
-2.58557011e-05, -5.12063357e-05, -4.87638925e-05,
1.19056128e-05, 4.62340053e-08, 1.59350048e-07,
3.09463332e-07, 3.13393204e-07, 0.00000000e+00]])}
%% Cell type:markdown id: tags:
Note that the ```cutoff``` array is now all 1, because SOAP takes the cutoff into account when it computes the descriptor vector
%% Cell type:markdown id: tags:
We now add a hydrogen atom to the structure
%% Cell type:code id: tags:
``` python
at.add_atoms((0.2,0.2,0.2),(1))
```
%% Cell type:code id: tags:
``` python
at.calc_connect()
desc.descriptor_sizes(at)
```
%% Output
(9, 285)
%% Cell type:markdown id: tags:
The descriptor sizes did not change! This is because the descriptor was set up above to only look at carbon atoms (Z=6). We need a new descriptor that takes account of H atoms.
%% Cell type:code id: tags:
``` python
desc = descriptors.Descriptor("soap cutoff=3 l_max=4 n_max=4 atom_sigma=0.5 n_Z=2 Z={1 6} n_species=2 species_Z={1 6}")
```
%% Cell type:markdown id: tags:
The specification of which atoms are used as SOAP centers is separate to the specification of which atoms are taken into account in the environment. The ``n_Z=2 Z={1 6}`` options specify that both H and C are to be taken as SOAP centers. The ```n_species=2 species_Z={1 6}``` options specify the atom types to be taken into account in the environment.
Now the dimensionality of the SOAP descriptor goes up:
%% Cell type:code id: tags:
``` python
desc.n_dim
```
%% Output
181
%% Cell type:markdown id: tags:
And the size also increases to 9 (and the cross terms go up too), reflecting the fact that there are 9 atoms altogether.
%% Cell type:code id: tags:
``` python
desc.descriptor_sizes(at)
```
%% Output
(9, 285)
%% Cell type:code id: tags:
``` python
desc.calc(at)
```
%% Output
{'cutoff': array([ 1., 1., 1., 1., 1., 1., 1., 1., 1.]),
'descriptor': array([[ 8.77557915e-02, 6.67877090e-06, 1.32368918e-06, ...,
5.39505036e-08, 3.31313515e-07, 0.00000000e+00],
[ 1.77188829e-03, 3.74461483e-04, 2.74157244e-04, ...,
5.81842224e-08, 3.57313054e-07, 0.00000000e+00],
[ 1.24011279e-04, 4.57421978e-05, 4.91840662e-05, ...,
5.81830007e-08, 3.57305552e-07, 0.00000000e+00],
...,
[ 1.24011279e-04, 4.57421978e-05, 4.91840662e-05, ...,
5.81830007e-08, 3.57305552e-07, 0.00000000e+00],
[ 3.76907877e-04, 4.04885591e-04, 3.44444944e-04, ...,
5.79365059e-08, 3.55791811e-07, 0.00000000e+00],
[ 1.35807048e-01, 0.00000000e+00, 0.00000000e+00, ...,
3.44738745e-07, 1.39187207e-07, 0.00000000e+00]])}
%% Cell type:markdown id: tags:
Note that the same descriptor can be obtained not via python, but from the unix command line, using the ```quip``` program. The input needs to be an extended XYZ file, and the descriptor specification is identical to the above.
So let's write out an xyz file:
%% Cell type:code id: tags:
``` python
at.write("diamondH.xyz")
```
%% Cell type:markdown id: tags:
The following command will print all the descriptor vectors, each vector on a line, which starts with the string ```DESC```
```quip atoms_filename=diamondH.xyz descriptor_str="soap cutoff=3 l_max=4 n_max=4 atom_sigma=0.5 n_Z=2 Z={1 6} n_species=2 species_Z={1 6}" ```
%% Cell type:markdown id: tags:
## General_monomer example: benzene
%% Cell type:markdown id: tags:
Now for a more interesting descriptor: a three-site benzene monomer
%% Cell type:code id: tags:
``` python
desc_monomer = descriptors.Descriptor('general_monomer signature={6 6 6} atom_ordercheck=F')
```
%% Cell type:markdown id: tags:
Three intramolecular distances = dimensionality 3:
%% Cell type:code id: tags:
``` python
desc_monomer.n_dim
```
%% Output
3
%% Cell type:markdown id: tags:
And a lot more permutations now:
%% Cell type:code id: tags:
``` python
desc_monomer.n_perm
```
%% Output
6
%% Cell type:code id: tags:
``` python
desc_monomer.permutations()
```
%% Output
array([[1, 2, 3],
[2, 1, 3],
[1, 3, 2],
[3, 1, 2],
[2, 3, 1],
[3, 2, 1]], dtype=int32)
%% Cell type:markdown id: tags:
Load some test configurations
%% Cell type:code id: tags:
``` python
benzat = quippy.AtomsList('benzene_frames.xyz')
```
%% Cell type:code id: tags:
``` python
benzat.pos.shape
```
%% Output
(40, 3, 648)
%% Cell type:code id: tags:
``` python
for i in range(len(benzat)):
benzat[i].set_cutoff(desc_monomer.cutoff()+0.5)
benzat[i].calc_connect()
```
%% Cell type:markdown id: tags:
Calling `calc` or any of the other methods on an `AtomsList` returns a list of results. Here we look at the results for the first configuration:
%% Cell type:code id: tags:
``` python
res = desc_monomer.calc(benzat, grad=True)[0]
```
%% Cell type:code id: tags:
``` python
res.grad.shape
```
%% Output
(648, 3, 3)
3rdparty/gap/doc/tutorials.rst deleted 100644 → 0
View file @ 23bcf97f
Tutorials
=========
The following tutorials are currently available:
.. toctree::
:maxdepth: 1
gap_fitting_tutorial.ipynb
quippy-descriptor-tutorial.ipynb
3rdparty/gap/error.inc deleted 100644 → 0
View file @ 23bcf97f
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X Portions of this code were written by
! H0 X Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X Copyright 2006-2010.
! H0 X
! H0 X These portions of the source code are released under the GNU General
! H0 X Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X If you would like to license the source code under different terms,
! H0 X please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X Portions of this code were written by Noam Bernstein as part of
! H0 X his employment for the U.S. Government, and are not subject
! H0 X to copyright in the USA.
! H0 X
! H0 X
! H0 X When using this software, please cite the following reference:
! H0 X
! H0 X http://www.libatoms.org
! H0 X
! H0 X Additional contributions by
! H0 X Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!X
!X Error handling, see error.f95 for the functions called in these macros.
!X
!X Error passing works as follows:
!X - *error* needs to be intent(out) and optional
!X - all functions that receive *error* as an argument must call INIT_ERROR(error)
!X - RAISE_ERROR is used whenever an error occurs. If *error* is not present,
!X the program execution will be terminated immediately. If *error* is
!X present it will be set to some value not equal ERROR_NONE and the execution
!X of the subroutine will be stopped.
!X - PASS_ERROR is used after a function or subroutine that returns error, i.e.
!X call sub(..., error=error)
!X PASS_ERROR(error)
!X If no error occurs (i.e. error==ERROR_NONE), execution will proceed as
!X usual. If an error occured, the current function will be terminated after
!X the location of the error is passed to the error module.
!X If the calling routine handles the error itself, rather than passing
!X it up with PASS_ERROR(), CLEAR_ERROR() should be used to clear the error
!X info stack
!X - PASS_ERROR_WITH_INFO is like PASS_ERROR, just an additional string can be
!X provided describing the error, or parameters.
!X - HANDLE_ERROR will print the error history and stop execution of the program
!X after an error occured.
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#define INIT_ERROR(error) if (present(error)) then ; error = ERROR_NONE ; endif
#define ASSERT(condition, message, error) if (.not. (condition)) then ; RAISE_ERROR(message, error) ; endif
#define RAISE_ERROR(message, error) if (.true.) then ; call push_error_with_info(message, __FILE__, __LINE__) ; if (present(error)) then ; error = ERROR_UNSPECIFIED ; return ; else ; call error_abort(error) ; endif ; endif
#define RAISE_ERROR_WITH_KIND(kind, message, error) if (.true.) then ; call push_error_with_info(message, __FILE__, __LINE__, kind) ; if (present(error)) then ; error = kind ; return ; else ; call error_abort(error) ; endif ; endif
#define PASS_ERROR(error) if (present(error)) then ; if (error /= ERROR_NONE) then ; call push_error(__FILE__, __LINE__) ; return ; endif ; endif
#define PASS_ERROR_WITH_INFO(message, error) if (present(error)) then ; if (error /= ERROR_NONE) then ; call push_error_with_info(message, __FILE__, __LINE__) ; return ; endif ; endif
#define HANDLE_ERROR(error) if (error /= ERROR_NONE) then ; call push_error(__FILE__, __LINE__) ; call error_abort(error) ; endif
#define CLEAR_ERROR(error) call error_clear_stack()
#define PRINT_LINE_NUMBER if(.true.) then; print "('LINE ' i0)",__LINE__; endif
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!X
!X MPI errors
!X
!X MPI error string are obtained using mpi_error_string and then pushed
!X onto the error stack.
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#define PASS_MPI_ERROR(mperror, error) if (mperror /= MPI_SUCCESS) then ; call push_MPI_error(mperror, __FILE__, __LINE__) ; if (present(error)) then ; error = ERROR_MPI ; return ; else ; call error_abort(error) ; endif ; endif
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!X
!X MPI BCAST errors
!X
!X Extension of error handling macros to MPI cases where processes
!X perform different tasks. If an error occurs on one process it will
!X be broadcast to all others before the error is propagated
!X upwards. Replace RAISE_ERROR with BCAST_RAISE_ERROR and PASS_ERROR
!X with BCAST_PASS_ERROR. Additionally, BCAST_CHECK_ERROR must be
!X called on the processes in which no error has occured. See
!X CInOutput read() for an example usage of these macros.
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#define BCAST_ASSERT(condition, message, error, mpi) if (.not. (condition)) then ; BCAST_RAISE_ERROR(message, error, mpi) ; endif
#define BCAST_RAISE_ERROR(message, error, mpi) if (present(error)) call bcast(mpi, error); RAISE_ERROR(message, error)
#define BCAST_RAISE_ERROR_WITH_KIND(kind, message, error, mpi) if (present(error)) then; error = kind; call bcast(mpi, error); end if; RAISE_ERROR_WITH_KIND(kind, message, error)
#define BCAST_PASS_ERROR(error, mpi) if (present(error)) then; if (error /= ERROR_NONE) call bcast(mpi, error); endif; PASS_ERROR(error)
#define BCAST_CHECK_ERROR(error, mpi) if (present(error)) then; call bcast(mpi, error); if (error /= ERROR_NONE) then; RAISE_ERROR_WITH_KIND(error, "An error occured on another MPI process", error); endif; endif
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!X
!X Delayed errors - for OpenMP loops
!X
!X A subroutine currently in an OpenMP section cannot be quit using
!X the *return* statement. Hence, the error flag is set using
!X RAISE_DELAYED_ERROR and TRACE_DELAYED_ERROR. After the OpenMP section
!X has finished, INVOKE_DELAYED_ERROR will raise the error and exit
!X the current subroutine if an error occured in the OpenMP section.
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#define RAISE_DELAYED_ERROR(message, error_loc) if (error_loc == ERROR_NONE) then ; call push_error_with_info(message, __FILE__, __LINE__) ; error_loc = ERROR_UNSPECIFIED ; endif
#define TRACE_DELAYED_ERROR(error_loc) if (error_loc /= ERROR_NONE) then ; call push_error(__FILE__, __LINE__) ; endif
#define TRACE_DELAYED_ERROR_WITH_INFO(message, error_loc) if (error_loc /= ERROR_NONE) then ; call push_error_with_info(message, __FILE__, __LINE__) ; endif
#define INVOKE_DELAYED_ERROR(error_loc, error) if (error_loc /= ERROR_NONE) then ; call push_error(__FILE__, __LINE__) ; if (present(error)) then ; error = error_loc ; else ; call error_abort(error) ; endif ; endif
3rdparty/gap/gap_fit.f95 deleted 100644 → 0
View file @ 23bcf97f
! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! HND X
! HND X libAtoms+QUIP: atomistic simulation library
! HND X
! HND X Portions of this code were written by
! HND X Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! HND X Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! HND X
! HND X Copyright 2006-2010.
! HND X
! HND X Not for distribution
! HND X
! HND X Portions of this code were written by Noam Bernstein as part of
! HND X his employment for the U.S. Government, and are not subject
! HND X to copyright in the USA.
! HND X
! HND X When using this software, please cite the following reference:
! HND X
! HND X http://www.libatoms.org
! HND X
! HND X Additional contributions by
! HND X Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! HND X
! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
program gap_fit_program
use libatoms_module
use descriptors_module
use gp_predict_module
use gp_fit_module
use clustering_module
use gap_fit_module
implicit none
type(gap_fit) :: main_gap_fit
call system_initialise(verbosity=PRINT_NORMAL, enable_timing=.false.)
call gap_fit_parse_command_line(main_gap_fit)
call gap_fit_parse_gap_str(main_gap_fit)
if(main_gap_fit%do_core) call read(main_gap_fit%quip_string, trim(main_gap_fit%core_param_file), keep_lf=.true.)
call add_template_string(main_gap_fit) ! if descriptor requires a template xyz file and this is provided, write to a string and add to descriptor_str
call read_descriptors(main_gap_fit) ! initialises descriptors from the descriptor_str and sets max_cutoff according to that.
call read_fit_xyz(main_gap_fit) ! reads in xyz into an array of atoms objects. sets cutoff and does calc_connect on each frame
call print('XYZ file read')
call get_species_xyz(main_gap_fit) ! counts the number of species present in the xyz file.
call add_multispecies_gaps(main_gap_fit)
call parse_config_type_sigma(main_gap_fit)
call parse_config_type_n_sparseX(main_gap_fit)
if(any(main_gap_fit%add_species)) then ! descriptor_str might have changed. reinitialises descriptors from the descriptor_str and sets max_cutoff according to that.
call read_descriptors(main_gap_fit)
endif
call print('Multispecies support added where requested')
call fit_n_from_xyz(main_gap_fit) ! counts number of energies, forces, virials. computes number of descriptors and gradients.
call set_baselines(main_gap_fit) ! sets e0 etc.
call fit_data_from_xyz(main_gap_fit) ! converts atomic neighbourhoods (bond neighbourhoods etc.) do descriptors, and feeds those to the GP
call print('Cartesian coordinates transformed to descriptors')
if(main_gap_fit%sparsify_only_no_fit) then
call system_finalise()
stop
end if
call enable_timing()
call system_timer('GP sparsify')
call gp_covariance_sparse(main_gap_fit%my_gp)
call initialise(main_gap_fit%gp_sp,main_gap_fit%my_gp)
call gap_fit_print_xml(main_gap_fit,main_gap_fit%gp_file,main_gap_fit%sparseX_separate_file)
call system_timer('GP sparsify')
call system_finalise()
end program gap_fit_program
3rdparty/gap/gap_fit_module.f95 deleted 100644 → 0
View file @ 23bcf97f
! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! HND X
! HND X libAtoms+QUIP: atomistic simulation library
! HND X
! HND X Portions of this code were written by
! HND X Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! HND X Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! HND X
! HND X Copyright 2006-2010.
! HND X
! HND X Not for distribution
! HND X
! HND X Portions of this code were written by Noam Bernstein as part of
! HND X his employment for the U.S. Government, and are not subject
! HND X to copyright in the USA.
! HND X
! HND X When using this software, please cite the following reference:
! HND X
! HND X http://www.libatoms.org
! HND X
! HND X Additional contributions by
! HND X Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! HND X
! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#include "error.inc"
module gap_fit_module
use error_module
use libatoms_module
use descriptors_module
use gp_predict_module
use gp_fit_module
use fox_wxml
use potential_module
implicit none
integer, parameter :: SPARSE_LENGTH = 10000
integer, parameter :: THETA_LENGTH = 10000
integer, parameter :: E0_ISOLATED = 1
integer, parameter :: E0_AVERAGE = 2
#ifdef GAP_VERSION
integer, parameter, private :: gap_version = GAP_VERSION
#else
integer, parameter, private :: gap_version = huge(1)
#endif
type gap_fit
!% everything from the command line
type(Atoms), dimension(:), allocatable :: at
character(len=STRING_LENGTH) :: at_file='', core_ip_args = '', e0_str, local_property0_str, &
energy_parameter_name, local_property_parameter_name, force_parameter_name, virial_parameter_name, &
hessian_parameter_name, config_type_parameter_name, sigma_parameter_name, &
config_type_sigma_string, core_param_file, gp_file, template_file
character(len=10240) :: command_line = ''
real(dp), dimension(total_elements) :: e0, local_property0
real(dp) :: max_cutoff
real(dp), dimension(4) :: default_sigma
real(dp) :: sparse_jitter, e0_offset, hessian_delta
integer :: e0_method = E0_ISOLATED
logical :: do_core = .false., do_copy_at_file, has_config_type_sigma, sigma_per_atom = .true.
logical :: sparsify_only_no_fit = .false.
integer :: n_frame, n_coordinate, n_ener, n_force, n_virial, n_hessian, n_local_property, min_save, n_species
type(extendable_str) :: quip_string
type(Potential) :: core_pot
type(gpFull) :: my_gp
type(gpSparse) :: gp_sp
type(descriptor), dimension(:), allocatable :: my_descriptor
character(len=STRING_LENGTH), dimension(99) :: gap_str
real(dp), dimension(:), allocatable :: delta, f0, theta_uniform, zeta, unique_hash_tolerance, unique_descriptor_tolerance !, theta
real(dp), dimension(:,:), allocatable :: sigma
integer, dimension(:), allocatable :: n_sparseX, sparse_method, target_type, n_cross, n_descriptors, species_Z, covariance_type
integer, dimension(:,:), allocatable :: config_type_n_sparseX
character(len=STRING_LENGTH), dimension(:), allocatable :: theta_file, sparse_file, theta_fac_string, config_type, config_type_n_sparseX_string, print_sparse_index
logical, dimension(:), allocatable :: mark_sparse_atoms, add_species, has_theta_fac, has_theta_uniform, has_theta_file, has_zeta
logical :: sparseX_separate_file
logical :: sparse_use_actual_gpcov
logical :: has_template_file, has_e0, has_local_property0, has_e0_offset
endtype gap_fit
private
public :: fit_n_from_xyz
public :: fit_data_from_xyz
public :: e0_from_xyz
public :: w_Z_from_xyz
public :: gap_fit
public :: gap_fit_print_xml
public :: file_print_xml
! public :: print_sparse
public :: set_baselines
public :: parse_config_type_sigma
public :: parse_config_type_n_sparseX
public :: read_fit_xyz
public :: read_descriptors
public :: get_species_xyz
public :: add_multispecies_gaps
public :: add_template_string
public :: gap_fit_parse_command_line
public :: gap_fit_parse_gap_str
contains
subroutine gap_fit_parse_command_line(this)
!% This subroutine parses the main command line options.
type(gap_fit), intent(inout), target :: this
type(Dictionary) :: params
character(len=STRING_LENGTH), pointer :: at_file, e0_str, local_property0_str, &
core_param_file, core_ip_args, &
energy_parameter_name, local_property_parameter_name, force_parameter_name, &
virial_parameter_name, hessian_parameter_name, &
config_type_parameter_name, sigma_parameter_name, config_type_sigma_string, &
gp_file, template_file
character(len=STRING_LENGTH) :: gap_str, verbosity, sparse_method_str, covariance_type_str, e0_method
logical, pointer :: sigma_per_atom, do_copy_at_file, sparseX_separate_file, sparse_use_actual_gpcov
logical :: do_ip_timing, has_sparse_file, has_theta_uniform, has_at_file, has_gap, has_default_sigma
logical, pointer :: sparsify_only_no_fit
real(dp), pointer :: e0_offset, sparse_jitter, hessian_delta
real(dp), dimension(:), pointer :: default_sigma
integer :: rnd_seed
at_file => this%at_file
e0_str => this%e0_str
local_property0_str => this%local_property0_str
e0_offset => this%e0_offset
default_sigma => this%default_sigma
sparse_jitter => this%sparse_jitter
hessian_delta => this%hessian_delta
core_param_file => this%core_param_file
core_ip_args => this%core_ip_args
energy_parameter_name => this%energy_parameter_name
local_property_parameter_name => this%local_property_parameter_name
force_parameter_name => this%force_parameter_name
virial_parameter_name => this%virial_parameter_name
hessian_parameter_name => this%hessian_parameter_name
config_type_parameter_name => this%config_type_parameter_name
sigma_parameter_name => this%sigma_parameter_name
config_type_sigma_string => this%config_type_sigma_string
sigma_per_atom => this%sigma_per_atom
do_copy_at_file => this%do_copy_at_file
sparseX_separate_file => this%sparseX_separate_file
sparse_use_actual_gpcov => this%sparse_use_actual_gpcov
gp_file => this%gp_file
template_file => this%template_file
sparsify_only_no_fit => this%sparsify_only_no_fit
call initialise(params)
call param_register(params, 'at_file', '//MANDATORY//', at_file, has_value_target = has_at_file, help_string="XYZ file with fitting configurations")
call param_register(params, 'gap', '//MANDATORY//', gap_str, has_value_target = has_gap, help_string="Initialisation string for GAPs")
call param_register(params, 'e0', '0.0', e0_str, has_value_target = this%has_e0, &
help_string="Atomic energy value to be subtracted from energies before fitting (and added back on after prediction). &
& Specifiy a single number (used for all species) or by species: {Ti:-150.0:O:-320}. energy = core + GAP + e0")
call param_register(params, 'local_property0', '0.0', local_property0_str, has_value_target = this%has_local_property0, &
help_string="Local property value to be subtracted from the local property before fitting (and added back on after prediction). &
& Specifiy a single number (used for all species) or by species: {H:20.0:Cl:35.0}.")
call param_register(params, 'e0_offset', '0.0', e0_offset, has_value_target = this%has_e0_offset, &
help_string="Offset of baseline. If zero, the offset is the average atomic energy of the input data or the e0 specified manually.")
call param_register(params, 'e0_method','isolated',e0_method, &
help_string="Method to determine e0, if not explicitly specified. Possible options: isolated (default, each atom &
present in the XYZ needs to have an isolated representative, with a valid energy), average (e0 is the average of &
all total energies across the XYZ)")
call param_register(params, 'default_sigma', '//MANDATORY//', default_sigma, has_value_target = has_default_sigma, &
help_string="error in [energies forces virials hessians]")
call param_register(params, 'sparse_jitter', "1.0e-10", sparse_jitter, &
help_string="intrisic error of atomic/bond energy, used to regularise the sparse covariance matrix")
call param_register(params, 'hessian_delta', "1.0e-2", hessian_delta, &
help_string="Delta to use in numerical differentiation when obtaining second derivative for the Hessian covariance")
call param_register(params, 'core_param_file', 'quip_params.xml', core_param_file, &
help_string="QUIP XML file for a potential to subtract from data (and added back after prediction)")
call param_register(params, 'core_ip_args', '', core_ip_args, has_value_target = this%do_core, &
help_string=" QUIP init string for a potential to subtract from data (and added back after prediction)")
call param_register(params, 'energy_parameter_name', 'energy', energy_parameter_name, &
help_string="Name of energy property in the at_file that describes the data")
call param_register(params, 'local_property_parameter_name', 'local_property', local_property_parameter_name, &
help_string="Name of local_property in the at_file that describes the data")
call param_register(params, 'force_parameter_name', 'force', force_parameter_name, &
help_string="Name of force property in the at_file that describes the data")
call param_register(params, 'virial_parameter_name', 'virial', virial_parameter_name, &
help_string="Name of virial property in the at_file that describes the data")
call param_register(params, 'hessian_parameter_name', 'hessian', hessian_parameter_name, &
help_string="Name of hessian property in the at_file that describes the data")
call param_register(params, 'config_type_parameter_name', 'config_type', config_type_parameter_name, &
help_string="Identifier of property determining the type of input data in the at_file")
call param_register(params, 'sigma_parameter_name', 'sigma', sigma_parameter_name, &
help_string="sigma parameters (error hyper) for a given configuration in the database. &
Overrides the command line sigmas. In the XYZ, it must be prepended by energy_, force_, virial_ or hessian_")
call param_register(params, 'config_type_sigma', '', config_type_sigma_string, has_value_target = this%has_config_type_sigma, &
help_string="What sigma values to choose for each type of data. Format: {type:energy:force:virial:hessian}")
call param_register(params, 'sigma_per_atom', 'T', sigma_per_atom, &
help_string="Interpretation of the energy and virial sigmas specified in >>default_sigma<< and >>config_type_sigma<<. &
If >>T<<, they are interpreted as per-atom errors, and the variance will be scaled according to the number of atoms in the configuration. &
If >>F<< they are treated as absolute errors and no scaling is performed. &
NOTE: sigmas specified on a per-configuration basis (see >>sigma_parameter_name<<) are always absolute.")
call param_register(params, 'do_copy_at_file', 'T', do_copy_at_file, &
help_string="Do copy the at_file into the GAP XML file (should be set to False for NetCDF input).")
call param_register(params, 'sparse_separate_file', 'T', sparseX_separate_file, &
help_string="Save sparse coordinates data in separate file")
call param_register(params, 'sparse_use_actual_gpcov', 'F', sparse_use_actual_gpcov, &
help_string="Use actual GP covariance for sparsification methods")
call param_register(params, 'gp_file', 'gp_new.xml', gp_file, help_string="output XML file")
call param_register(params, 'verbosity', 'NORMAL', verbosity, &
help_string="Verbosity control. Options: NORMAL, VERBOSE, NERD, ANAL.")
call param_register(params, "rnd_seed", "-1", rnd_seed, &
help_string="Random seed.")
call param_register(params, 'do_ip_timing', 'F', do_ip_timing, &
help_string="To enable or not timing of the interatomic potential.")
call param_register(params, 'template_file', 'template.xyz', template_file, has_value_target=this%has_template_file, &
help_string="Template XYZ file for initialising object")
call param_register(params, 'sparsify_only_no_fit', 'F', sparsify_only_no_fit, &
help_string="If true, sparsification is done, but no fitting. print the sparse index by adding print_sparse_index=file.dat to the descriptor string.")
if (.not. param_read_args(params, command_line=this%command_line)) then
call print("gap_fit")
call system_abort('Exit: Mandatory argument(s) missing...')
endif
call print_title("Input parameters")
call param_print(params)
call print_title("")
call finalise(params)
if( len_trim(this%gp_file) > 216 ) then ! The filename's length is limited to 255 char.s in some filesystem.
! Without this check, the fit would run but produce a core file and only a temporary xml file.
! The limit is set to 216 as the sparse file can be 39 characters longer.
call system_abort("gp_file's name "//this%gp_file//" is too long. Please start the fit again with a shorter name.")
endif
if(do_ip_timing) call enable_timing()
select case(verbosity)
case ("NORMAL")
call verbosity_push(PRINT_NORMAL)
case ("VERBOSE")
call verbosity_push(PRINT_VERBOSE)
case ("NERD")
call verbosity_push(PRINT_NERD)
case ("ANAL")
call verbosity_push(PRINT_ANAL)
case default
call system_abort("confused by verbosity " // trim(verbosity))
end select
select case(lower_case(e0_method))
case ("isolated")
this%e0_method = E0_ISOLATED
case ("average")
this%e0_method = E0_AVERAGE
case default
call system_abort("confused by e0_method " // trim(e0_method))
end select
if (rnd_seed >= 0) call system_set_random_seeds(rnd_seed)
call print_title('Gaussian Approximation Potentials - Database fitting')
call print('')
call print('Initial parsing of command line arguments finished.')
call split_string(gap_str,':;','{}',this%gap_str(:),this%n_coordinate,matching=.true.)
call print('Found '//this%n_coordinate//' GAPs.')
endsubroutine gap_fit_parse_command_line
subroutine set_baselines(this)
type(gap_fit), intent(inout) :: this
integer :: i
this%e0 = 0.0_dp
if( count( (/this%has_e0, this%has_e0_offset/) ) > 1 ) then
call print_warning('Both e0 and e0_offset has been specified. That means your atomic energy is e0 + e0_offset')
endif
if( this%has_e0 ) then
call parse_atomtype_value_str(this%e0_str,this%e0)
else
call e0_from_xyz(this) ! calculates the average atomic energy so it can be subtracted later.
endif
if( this%has_e0_offset ) this%e0 = this%e0 + this%e0_offset
if( .not. this%has_e0 ) then
do i = 1, size(this%e0)
if( all(i/=this%species_Z) ) this%e0(i) = 0.0_dp
enddo
call print('E0/atom = '//this%e0)
endif
if( this%has_local_property0 ) then
call parse_atomtype_value_str(this%local_property0_str,this%local_property0)
this%e0 = 0.0_dp
else
this%local_property0 = 0.0_dp
endif
endsubroutine set_baselines
subroutine parse_atomtype_value_str(this,values,error)
character(len=STRING_LENGTH), intent(in) :: this
real(dp), dimension(total_elements), intent(out) :: values
integer, intent(out), optional :: error
integer :: n_string_array, i, z
character(len=STRING_LENGTH), dimension(2*total_elements) :: string_array
INIT_ERROR(error)
call split_string(this,':','{}',string_array(:),n_string_array,matching=.true.)
if(n_string_array == 1) then
values = string_to_real(trim(string_array(1)))
elseif(mod(n_string_array,2) == 0) then
values = 0.0_dp
do i = 1, n_string_array / 2
z = atomic_number(trim( string_array((i-1)*2+1) ))
if( z==0 ) then
RAISE_ERROR("parse_atomtype_value_str: invalid atomic symbol "//trim(string_array((i-1)*2+1)),error)
endif
values(z) = string_to_real(trim( string_array(2*i) ))
enddo
else
RAISE_ERROR("parse_atomtype_value_str: number of fields is an odd number. It must be a list of pairs of values, such as {Ti:-150.4:O:-345.1}",error)
endif
endsubroutine parse_atomtype_value_str
subroutine gap_fit_parse_gap_str(this)
!% This subroutine parses the options given in the gap string, for each GAP.
type(gap_fit), intent(inout), target :: this
type(Dictionary) :: params
integer :: i_coordinate
real(dp) :: delta, f0, theta_uniform, zeta, unique_hash_tolerance, unique_descriptor_tolerance
integer :: n_sparseX, sparse_method, covariance_type
character(len=STRING_LENGTH) :: config_type_n_sparseX_string, theta_fac_string, theta_file, sparse_file, print_sparse_index, &
covariance_type_str, sparse_method_str
logical :: mark_sparse_atoms, add_species, has_sparse_file
allocate(this%delta(this%n_coordinate))
allocate(this%f0(this%n_coordinate))
allocate(this%n_sparseX(this%n_coordinate))
allocate(this%config_type_n_sparseX_string(this%n_coordinate))
allocate(this%theta_fac_string(this%n_coordinate))
allocate(this%theta_uniform(this%n_coordinate))
allocate(this%theta_file(this%n_coordinate))
allocate(this%has_theta_fac(this%n_coordinate))
allocate(this%has_theta_uniform(this%n_coordinate))
allocate(this%has_theta_file(this%n_coordinate))
allocate(this%sparse_file(this%n_coordinate))
allocate(this%mark_sparse_atoms(this%n_coordinate))
allocate(this%sparse_method(this%n_coordinate))
allocate(this%add_species(this%n_coordinate))
allocate(this%covariance_type(this%n_coordinate))
allocate(this%zeta(this%n_coordinate))
allocate(this%has_zeta(this%n_coordinate))
allocate(this%print_sparse_index(this%n_coordinate))
allocate(this%unique_hash_tolerance(this%n_coordinate))
allocate(this%unique_descriptor_tolerance(this%n_coordinate))
do i_coordinate = 1, this%n_coordinate
call initialise(params)
call param_register(params, 'delta', "//MANDATORY//", delta, &
help_string="Set the standard deviation of the Gaussian process. Typically this would be &
set to the standard deviation (i.e. root mean square) of the function &
that is approximated with the Gaussian process.")
call param_register(params, 'f0', '0.0', f0, &
help_string="Set the mean of the Gaussian process. Defaults to 0.")
call param_register(params, 'n_sparse', "0", n_sparseX, &
help_string="Number of sparse points to use in the sparsification of the Gaussian process")
call param_register(params, 'config_type_n_sparse', '', config_type_n_sparseX_string, &
help_string="Number of sparse points in each config type. Format: {type1:50:type2:100}")
call param_register(params, 'sparse_method', 'RANDOM', sparse_method_str, &
help_string="Sparsification method. RANDOM(default), PIVOT, CLUSTER, UNIFORM, KMEANS, COVARIANCE, NONE, FUZZY, FILE, &
INDEX_FILE, CUR_COVARIANCE, CUR_POINTS")
call param_register(params, 'theta_fac', '1.0', theta_fac_string, has_value_target = this%has_theta_fac(i_coordinate), &
help_string="Set the width of Gaussians for the ARD_SE and PP kernel by multiplying the range of each descriptor by theta_fac. &
Can be a single number or different for each dimension. For multiple theta_fac separate each value by whitespaces.")
call param_register(params, 'theta_uniform', '0.0', theta_uniform, has_value_target = this%has_theta_uniform(i_coordinate), &
help_string="Set the width of Gaussians for the ARD_SE and PP kernel, same in each dimension.")
call param_register(params, 'theta_file', '', theta_file, has_value_target = this%has_theta_file(i_coordinate), &
help_string="Set the width of Gaussians for the ARD_SE kernel from a file. &
There should be as many real numbers as the number of dimensions, in a single line")
call param_register(params, 'sparse_file', '', sparse_file, has_value_target = has_sparse_file, &
help_string="Sparse points from a file. If sparse_method=FILE, descriptor values (real) listed in a text file, one &
& >>element<< per line. If sparse_method=INDEX_FILE, 1-based index of sparse points, one per line.")
call param_register(params, 'mark_sparse_atoms', 'F', mark_sparse_atoms, &
help_string="Reprints the original xyz file after sparsification process. &
sparse propery added, true for atoms associated with a sparse point.")
call param_register(params, 'add_species', 'F', add_species, &
help_string="Create species-specific descriptor, using the descriptor string as a template.")
call param_register(params, 'covariance_type', "//MANDATORY//", covariance_type_str, &
help_string="Type of covariance function to use. Available: ARD_SE, DOT_PRODUCT, BOND_REAL_SPACE, PP (piecewise polynomial)")
!call param_register(params, 'theta', '1.0', main_gap_fit%theta(i_coordinate), &
!help_string="Width of Gaussians for use with bond real space covariance.")
call param_register(params, 'zeta', '1.0', zeta, has_value_target = this%has_zeta(i_coordinate), &
help_string="Exponent of soap type dot product covariance kernel")
call param_register(params, 'print_sparse_index', '', print_sparse_index, &
help_string="If given, after determinining the sparse points, their 1-based indices are appended to this file")
call param_register(params, 'unique_hash_tolerance', '1.0e-10', unique_hash_tolerance, &
help_string="Hash tolerance when filtering out duplicate data points")
call param_register(params, 'unique_descriptor_tolerance', '1.0e-10', unique_descriptor_tolerance, &
help_string="Descriptor tolerance when filtering out duplicate data points")
if (.not. param_read_line(params, this%gap_str(i_coordinate), ignore_unknown=.true., task='main program gap_str('//i_coordinate//')')) then
call system_abort("main program failed to parse gap string ("//i_coordinate//")='"//trim(this%gap_str(i_coordinate))//"'")
endif
call finalise(params)
this%delta(i_coordinate) = delta
this%f0(i_coordinate) = f0
this%n_sparseX(i_coordinate) = n_sparseX
this%config_type_n_sparseX_string(i_coordinate) = config_type_n_sparseX_string
this%theta_fac_string(i_coordinate) = theta_fac_string
this%theta_uniform(i_coordinate) = theta_uniform
this%theta_file(i_coordinate) = theta_file
this%sparse_file(i_coordinate) = sparse_file
this%mark_sparse_atoms(i_coordinate) = mark_sparse_atoms
this%add_species(i_coordinate) = add_species
this%zeta(i_coordinate) = zeta
this%print_sparse_index(i_coordinate) = print_sparse_index
this%unique_hash_tolerance(i_coordinate) = unique_hash_tolerance
this%unique_descriptor_tolerance(i_coordinate) = unique_descriptor_tolerance
select case(lower_case(trim(sparse_method_str)))
case('random')
this%sparse_method(i_coordinate) = GP_SPARSE_RANDOM
case('pivot')
this%sparse_method(i_coordinate) = GP_SPARSE_PIVOT
case('cluster')
this%sparse_method(i_coordinate) = GP_SPARSE_CLUSTER
case('uniform')
this%sparse_method(i_coordinate) = GP_SPARSE_UNIFORM
case('kmeans')
this%sparse_method(i_coordinate) = GP_SPARSE_KMEANS
case('covariance')
this%sparse_method(i_coordinate) = GP_SPARSE_COVARIANCE
case('uniq')
call system_abort("sparse method UNIQ is no longer in use. Use NONE instead." )
case('fuzzy')
this%sparse_method(i_coordinate) = GP_SPARSE_FUZZY
case('file')
this%sparse_method(i_coordinate) = GP_SPARSE_FILE
case('index_file')
this%sparse_method(i_coordinate) = GP_SPARSE_INDEX_FILE
case('cur_covariance')
this%sparse_method(i_coordinate) = GP_SPARSE_CUR_COVARIANCE
case('cur_points')
this%sparse_method(i_coordinate) = GP_SPARSE_CUR_POINTS
case('none')
this%sparse_method(i_coordinate) = GP_SPARSE_NONE
case default
call system_abort("unknown sparse method "//trim(sparse_method_str))
endselect
if( has_sparse_file ) then
if( this%sparse_method(i_coordinate) /= GP_SPARSE_FILE .and. &
this%sparse_method(i_coordinate) /= GP_SPARSE_INDEX_FILE ) then
call system_abort('"sparse_file" specified in command line, but sparse method not "file" or "index_file"')
endif
endif
select case(lower_case(trim(covariance_type_str)))
case('none')
call system_abort("covariance type cannot be"//trim(covariance_type_str))
this%covariance_type(i_coordinate) = COVARIANCE_NONE
case('ard_se')
this%covariance_type(i_coordinate) = COVARIANCE_ARD_SE
case('dot_product')
this%covariance_type(i_coordinate) = COVARIANCE_DOT_PRODUCT
case('bond_real_space')
this%covariance_type(i_coordinate) = COVARIANCE_BOND_REAL_SPACE
case('pp')
this%covariance_type(i_coordinate) = COVARIANCE_PP
case default
call system_abort("unknown covariance type"//trim(covariance_type_str)//". Available: ARD_SE, DOT_PRODUCT, BOND_REAL_SPACE, PP (piecewise polynomial)")
endselect
enddo
call print('Descriptors have been parsed')
endsubroutine gap_fit_parse_gap_str
subroutine read_fit_xyz(this)
type(gap_fit), intent(inout) :: this
type(cinoutput) :: xyzfile
integer :: n_con
logical :: file_exists
if( allocated(this%at) ) then
do n_con = 1, this%n_frame
call finalise(this%at(n_con))
enddo
deallocate(this%at)
this%n_frame = 0
endif
inquire(file=this%at_file, exist=file_exists)
if( .not. file_exists ) then
call system_abort("read_fit_xyz: at_file "//this%at_file//" could not be found")
endif
call initialise(xyzfile,this%at_file)
this%n_frame = xyzfile%n_frame
allocate(this%at(this%n_frame))
do n_con = 1, this%n_frame
call read(xyzfile,this%at(n_con),frame=n_con-1)
call set_cutoff(this%at(n_con), this%max_cutoff)
call calc_connect(this%at(n_con))
enddo
call finalise(xyzfile)
endsubroutine read_fit_xyz
subroutine read_descriptors(this)
type(gap_fit), intent(inout) :: this
integer :: i
this%max_cutoff = 0.0_dp
if(allocated(this%my_descriptor)) then
do i = 1, size(this%my_descriptor)
call finalise(this%my_descriptor(i))
enddo
deallocate(this%my_descriptor)
endif
allocate(this%my_descriptor(this%n_coordinate))
do i = 1, this%n_coordinate
call initialise(this%my_descriptor(i),this%gap_str(i))
if( this%max_cutoff < cutoff(this%my_descriptor(i)) ) this%max_cutoff = cutoff(this%my_descriptor(i))
enddo
endsubroutine read_descriptors
subroutine fit_n_from_xyz(this)
type(gap_fit), intent(inout) :: this
type(Atoms) :: at
integer :: n_con
logical :: has_ener, has_force, has_virial, has_hessian, has_local_property
real(dp) :: ener, virial(3,3)
real(dp), pointer, dimension(:,:) :: f, hessian_eigenvector_j
real(dp), pointer, dimension(:) :: local_property
integer :: i, j, k
integer :: n_descriptors, n_cross, n_hessian
allocate(this%n_cross(this%n_coordinate))
allocate(this%n_descriptors(this%n_coordinate))
this%n_cross = 0
this%n_descriptors = 0
this%n_ener = 0
this%n_force = 0
this%n_virial = 0
this%n_hessian = 0
do n_con = 1, this%n_frame
has_ener = get_value(this%at(n_con)%params,this%energy_parameter_name,ener)
has_force = assign_pointer(this%at(n_con),this%force_parameter_name, f)
has_virial = get_value(this%at(n_con)%params,this%virial_parameter_name,virial)
has_hessian = get_value(this%at(n_con)%params,"n_"//this%hessian_parameter_name,n_hessian)
has_local_property = assign_pointer(this%at(n_con),this%local_property_parameter_name, local_property)
if( has_ener ) then
this%n_ener = this%n_ener + 1
endif
if( has_force ) then
this%n_force = this%n_force + this%at(n_con)%N*3
endif
if( has_virial ) then
this%n_virial = this%n_virial + 6
endif
if( has_hessian ) then
this%n_hessian = this%n_hessian + n_hessian
at = this%at(n_con)
endif
if( has_local_property ) then
this%n_local_property = this%n_local_property + this%at(n_con)%N
endif
if( has_local_property .and. ( has_ener .or. has_force .or. has_virial .or. has_hessian ) ) then
call system_abort("fit_n_from_xyz: local_property and (energy or force or virial or hessian) present in configuration, currently not allowed.")
endif
do i = 1, this%n_coordinate
call descriptor_sizes(this%my_descriptor(i),this%at(n_con),n_descriptors,n_cross)
if( has_force ) then
this%n_cross(i) = this%n_cross(i) + n_cross*3
endif
if( has_virial ) then
this%n_cross(i) = this%n_cross(i) + n_cross*6
endif
this%n_descriptors(i) = this%n_descriptors(i) + n_descriptors
if( has_hessian ) then
do j = 1, n_hessian
if( .not. assign_pointer(this%at(n_con),trim(this%hessian_parameter_name)//j, hessian_eigenvector_j) ) &
call system_abort("fit_n_from_xyz: could not find the "//j//"th of "//n_hessian//" hessian eigenvector")
hessian_eigenvector_j = hessian_eigenvector_j / sqrt( sum(hessian_eigenvector_j**2) )
do k = -1, 1, 2
at%pos = this%at(n_con)%pos + k * this%hessian_delta * hessian_eigenvector_j
call set_cutoff(at,this%max_cutoff)
call calc_connect(at)
call descriptor_sizes(this%my_descriptor(i),at,n_descriptors,n_cross)
this%n_descriptors(i) = this%n_descriptors(i) + n_descriptors
this%n_cross(i) = this%n_cross(i) + n_descriptors
enddo
enddo
endif
enddo
call finalise(at)
enddo
call print_title("Report on number of descriptors found")
do i = 1, this%n_coordinate
call print("---------------------------------------------------------------------")
call print("Descriptor: "//this%gap_str(i))
call print("Number of descriptors: "//this%n_descriptors(i))
call print("Number of partial derivatives of descriptors: "//this%n_cross(i))
enddo
call print_title("")
end subroutine fit_n_from_xyz
subroutine fit_data_from_xyz(this,error)
type(gap_fit), intent(inout) :: this
integer, optional, intent(out) :: error
type(inoutput) :: theta_inout
type(descriptor_data) :: my_descriptor_data
type(Atoms) :: at
integer :: d
integer :: n_con
logical :: has_ener, has_force, has_virial, has_hessian, has_local_property, &
has_config_type, has_energy_sigma, has_force_sigma, has_virial_sigma, has_hessian_sigma, &
has_force_atom_sigma
real(dp) :: ener, ener_core, my_cutoff, energy_sigma, force_sigma, virial_sigma, hessian_sigma, grad_covariance_cutoff, &
use_force_sigma
real(dp), dimension(3) :: pos
real(dp), dimension(3,3) :: virial, virial_core
real(dp), dimension(:), allocatable :: theta, theta_fac, hessian, hessian_core, grad_data
real(dp), dimension(:), pointer :: force_atom_sigma
real(dp), dimension(:,:), pointer :: f, hessian_eigenvector_i, f_hessian
real(dp), dimension(:), pointer :: local_property
real(dp), dimension(:,:), allocatable :: f_core
integer, dimension(:,:), allocatable :: force_loc, permutations
integer :: ie, i, j, n, k, l, i_coordinate, n_hessian, n_energy_sigma, n_force_sigma, n_force_atom_sigma, n_hessian_sigma, &
n_virial_sigma, n_descriptors
integer, dimension(:), allocatable :: xloc, hessian_loc, local_property_loc
integer, dimension(3,3) :: virial_loc
integer :: i_config_type, n_config_type, n_theta_fac
character(len=STRING_LENGTH) :: config_type
character(len=THETA_LENGTH) :: theta_string
character(len=STRING_LENGTH), dimension(:), allocatable :: theta_string_array
INIT_ERROR(error)
my_cutoff = 0.0_dp
call gp_setParameters(this%my_gp,this%n_coordinate,this%n_ener+this%n_local_property,this%n_force+this%n_virial+this%n_hessian,this%sparse_jitter)
do i_coordinate = 1, this%n_coordinate
d = descriptor_dimensions(this%my_descriptor(i_coordinate))
call gp_setParameters(this%my_gp,i_coordinate, d, this%n_descriptors(i_coordinate), this%n_cross(i_coordinate), this%delta(i_coordinate), this%f0(i_coordinate), &
covariance_type=this%covariance_type(i_coordinate) )
call gp_addDescriptor(this%my_gp,i_coordinate,trim(this%gap_str(i_coordinate)))
allocate(permutations(d,descriptor_n_permutations(this%my_descriptor(i_coordinate))))
call descriptor_permutations(this%my_descriptor(i_coordinate),permutations)
call gp_setPermutations(this%my_gp,i_coordinate,permutations)
deallocate(permutations)
my_cutoff = max(my_cutoff,cutoff(this%my_descriptor(i_coordinate)))
enddo
call print_title("Report on number of target properties found in training XYZ:")
call print("Number of target energies (property name: "//trim(this%energy_parameter_name)//") found: "//this%n_ener)
call print("Number of target local_properties (property name: "//trim(this%local_property_parameter_name)//") found: "//this%n_local_property)
call print("Number of target forces (property name: "//trim(this%force_parameter_name)//") found: "//this%n_force)
call print("Number of target virials (property name: "//trim(this%virial_parameter_name)//") found: "//this%n_virial)
call print("Number of target Hessian eigenvalues (property name: "//trim(this%hessian_parameter_name)//") found: "//this%n_hessian)
call print_title("End of report")
if( this%do_core ) call Initialise(this%core_pot, args_str=this%core_ip_args, param_str=string(this%quip_string))
n_energy_sigma = 0
n_force_sigma = 0
n_force_atom_sigma = 0
n_hessian_sigma = 0
n_virial_sigma = 0
do n_con = 1, this%n_frame
has_ener = get_value(this%at(n_con)%params,this%energy_parameter_name,ener)
has_force = assign_pointer(this%at(n_con),this%force_parameter_name, f)
has_virial = get_value(this%at(n_con)%params,this%virial_parameter_name,virial)
has_hessian = get_value(this%at(n_con)%params,"n_"//this%hessian_parameter_name,n_hessian)
has_config_type = get_value(this%at(n_con)%params,this%config_type_parameter_name,config_type)
has_local_property = assign_pointer(this%at(n_con),this%local_property_parameter_name,local_property)
has_energy_sigma = get_value(this%at(n_con)%params,'energy_'//trim(this%sigma_parameter_name),energy_sigma)
has_force_sigma = get_value(this%at(n_con)%params,'force_'//trim(this%sigma_parameter_name),force_sigma)
has_virial_sigma = get_value(this%at(n_con)%params,'virial_'//trim(this%sigma_parameter_name),virial_sigma)
has_hessian_sigma = get_value(this%at(n_con)%params,'hessian_'//trim(this%sigma_parameter_name),hessian_sigma)
has_force_atom_sigma = assign_pointer(this%at(n_con),'force_atom_'//trim(this%sigma_parameter_name),force_atom_sigma)
if( has_hessian ) then
allocate(hessian(n_hessian))
do i = 1, n_hessian
if( .not. get_value(this%at(n_con)%params,trim(this%hessian_parameter_name)//i,hessian(i)) ) &
call system_abort("fit_data_from_xyz: did not find "//i//"th of "//n_hessian//" hessian element" )
enddo
endif
if( has_config_type ) then
config_type = trim(config_type)
else
config_type = "default"
endif
if( .not. allocated(this%config_type) ) call system_abort('config_type not allocated')
n_config_type = 0
do i_config_type = 1, size(this%config_type)
if( trim(this%config_type(i_config_type)) == trim(config_type) ) n_config_type = i_config_type
enddo
if( n_config_type == 0 ) then ! get the number of the "default" type as default
do i_config_type = 1, size(this%config_type)
if( trim(this%config_type(i_config_type)) == "default" ) n_config_type = i_config_type
enddo
endif
if( this%do_core ) then
allocate( f_core(3,this%at(n_con)%N) )
ener_core = 0.0_dp
f_core = 0.0_dp
virial_core = 0.0_dp
if( this%at(n_con)%cutoff < max(cutoff(this%core_pot),my_cutoff) ) then
call set_cutoff(this%at(n_con), max(cutoff(this%core_pot),my_cutoff))
call calc_connect(this%at(n_con))
endif
if(has_virial .and. has_force) then
call calc(this%core_pot,this%at(n_con),energy=ener_core,force=f_core,virial=virial_core)
elseif(has_force) then
call calc(this%core_pot,this%at(n_con),energy=ener_core,force=f_core)
elseif(has_virial) then
call calc(this%core_pot,this%at(n_con),energy=ener_core,virial=virial_core)
else
call calc(this%core_pot,this%at(n_con),energy=ener_core)
end if
if(has_hessian) then
allocate( hessian_core(n_hessian), f_hessian(3,this%at(n_con)%N) )
hessian_core = 0.0_dp
at = this%at(n_con)
call set_cutoff(at, cutoff(this%core_pot))
do i = 1, n_hessian
if( .not. assign_pointer(this%at(n_con),trim(this%hessian_parameter_name)//i, hessian_eigenvector_i) ) &
call system_abort("fit_data_from_xyz: could not find "//i//"th of "//n_hessian//" hessian eigenvector.")
hessian_eigenvector_i = hessian_eigenvector_i / sqrt( sum(hessian_eigenvector_i**2) )
do j = -1, 1, 2
at%pos = this%at(n_con)%pos + j * this%hessian_delta * hessian_eigenvector_i
call calc_connect(at)
call calc(this%core_pot,at,force = f_hessian)
hessian_core(i) = hessian_core(i) + j * sum(f_hessian*hessian_eigenvector_i) / 2.0_dp / this%hessian_delta
enddo
enddo
call finalise(at)
hessian = hessian - hessian_core
deallocate(hessian_core, f_hessian)
endif
if(has_ener) ener = ener - ener_core
if(has_force) f = f - f_core
if(has_virial) virial = virial - virial_core
deallocate(f_core)
endif
if(has_ener) then
do i = 1, this%at(n_con)%N
ener = ener - this%e0(this%at(n_con)%Z(i))
enddo
endif
if(has_local_property) then
do i = 1, this%at(n_con)%N
local_property(i) = local_property(i) - this%local_property0(this%at(n_con)%Z(i))
enddo
endif
if( has_ener .and. has_local_property ) then
RAISE_ERROR("fit_data_from_xyz: energy and local_property both present in configuration, currently not allowed.",error)
endif
if( this%at(n_con)%cutoff < my_cutoff ) then
call set_cutoff(this%at(n_con),my_cutoff)
call calc_connect(this%at(n_con))
endif
if( .not. has_energy_sigma ) then
if( this%sigma_per_atom ) then
energy_sigma = this%sigma(1,n_config_type)*sqrt(1.0_dp * this%at(n_con)%N)
else
energy_sigma = this%sigma(1,n_config_type)
endif
else
n_energy_sigma = n_energy_sigma + 1
endif
if( .not. has_force_sigma ) then
force_sigma = this%sigma(2,n_config_type)
else
n_force_sigma = n_force_sigma + 1
endif
if( .not. has_virial_sigma ) then
if( this%sigma_per_atom ) then
virial_sigma = this%sigma(3,n_config_type)*sqrt(1.0_dp * this%at(n_con)%N)
else
virial_sigma = this%sigma(3,n_config_type)
endif
else
n_virial_sigma = n_virial_sigma + 1
endif
if( .not. has_hessian_sigma ) then
hessian_sigma = this%sigma(4,n_config_type)
else
n_hessian_sigma = n_hessian_sigma + 1
endif
if( has_ener ) then
if( energy_sigma .feq. 0.0_dp ) then
RAISE_ERROR("fit_data_from_xyz: too small energy_sigma ("//energy_sigma//"), should be greater than zero",error)
endif
ie = gp_addFunctionValue(this%my_gp,ener, energy_sigma)
elseif( has_local_property ) then
if( energy_sigma .feq. 0.0_dp ) then
RAISE_ERROR("fit_data_from_xyz: too small energy_sigma ("//energy_sigma//"), should be greater than zero",error)
endif
allocate(local_property_loc(this%at(n_con)%N))
do i = 1, this%at(n_con)%N
local_property_loc(i) = gp_addFunctionValue(this%my_gp,local_property(i),energy_sigma/this%at(n_con)%N)
enddo
endif
if(has_force) then
allocate(force_loc(3,this%at(n_con)%N))
do i = 1, this%at(n_con)%N
if (has_force_atom_sigma) then
use_force_sigma = force_atom_sigma(i)
n_force_atom_sigma = n_force_atom_sigma + 1
else
use_force_sigma = force_sigma
endif
if( use_force_sigma .feq. 0.0_dp ) then
RAISE_ERROR("fit_data_from_xyz: too small force_sigma ("//use_force_sigma//"), should be greater than zero",error)
endif
do k = 1, 3
force_loc(k,i) = gp_addFunctionDerivative(this%my_gp,-f(k,i),use_force_sigma)
enddo
enddo
endif
if(has_virial) then
! check if virial is symmetric
if( sum((virial - transpose(virial))**2) .fne. 0.0_dp ) &
call print_warning('virial not symmetric, now symmetrised')
! Now symmetrise matrix
virial = ( virial + transpose(virial) ) / 2.0_dp
if( virial_sigma .feq. 0.0_dp ) then
RAISE_ERROR("fit_data_from_xyz: too small virial_sigma ("//virial_sigma//"), should be greater than zero",error)
endif
do k = 1, 3
do l = k, 3
virial_loc(l,k) = gp_addFunctionDerivative(this%my_gp,-virial(l,k),virial_sigma)
enddo
enddo
endif
if(has_hessian) then
if( hessian_sigma .feq. 0.0_dp ) then
RAISE_ERROR("fit_data_from_xyz: too small hessian_sigma ("//hessian_sigma//"), should be greater than zero",error)
endif
allocate(hessian_loc(n_hessian))
do i = 1, n_hessian
hessian_loc(i) = gp_addFunctionDerivative(this%my_gp,hessian(i),hessian_sigma)
enddo
endif
n_descriptors = 0
do i_coordinate = 1, this%n_coordinate
call calc(this%my_descriptor(i_coordinate),this%at(n_con),my_descriptor_data, &
do_descriptor=.true.,do_grad_descriptor=has_force .or. has_virial)
allocate(xloc(size(my_descriptor_data%x)))
n_descriptors = n_descriptors + size(my_descriptor_data%x)
if( has_ener ) then
do i = 1, size(my_descriptor_data%x)
if( .not. my_descriptor_data%x(i)%has_data) cycle
xloc(i) = gp_addCoordinates(this%my_gp,my_descriptor_data%x(i)%data(:),i_coordinate, &
cutoff_in=my_descriptor_data%x(i)%covariance_cutoff, current_y=ie,config_type=n_config_type)
enddo
elseif( has_local_property ) then
do i = 1, size(my_descriptor_data%x)
if( .not. my_descriptor_data%x(i)%has_data) cycle
xloc(i) = gp_addCoordinates(this%my_gp,my_descriptor_data%x(i)%data(:),i_coordinate, &
cutoff_in=my_descriptor_data%x(i)%covariance_cutoff, current_y=local_property_loc(my_descriptor_data%x(i)%ci(1)),config_type=n_config_type)
enddo
else
do i = 1, size(my_descriptor_data%x)
if( .not. my_descriptor_data%x(i)%has_data) cycle
xloc(i) = gp_addCoordinates(this%my_gp,my_descriptor_data%x(i)%data(:),i_coordinate, &
cutoff_in=my_descriptor_data%x(i)%covariance_cutoff, config_type=n_config_type)
enddo
endif
if(has_force) then
do i = 1, size(my_descriptor_data%x)
do n = lbound(my_descriptor_data%x(i)%ii,1), ubound(my_descriptor_data%x(i)%ii,1)
if( .not. my_descriptor_data%x(i)%has_grad_data(n)) cycle
j = my_descriptor_data%x(i)%ii(n)
do k = 1, 3
call gp_addCoordinateDerivatives(this%my_gp,my_descriptor_data%x(i)%grad_data(:,k,n),i_coordinate, &
force_loc(k,j), xloc(i), dcutoff_in=my_descriptor_data%x(i)%grad_covariance_cutoff(k,n) )
enddo
enddo
enddo
endif
if(has_virial) then
do k = 1, 3
do l = k, 3
do i = 1, size(my_descriptor_data%x)
do n = lbound(my_descriptor_data%x(i)%ii,1), ubound(my_descriptor_data%x(i)%ii,1)
if( .not. my_descriptor_data%x(i)%has_grad_data(n)) cycle
j = my_descriptor_data%x(i)%ii(n)
pos = my_descriptor_data%x(i)%pos(:,n)
call gp_addCoordinateDerivatives(this%my_gp,my_descriptor_data%x(i)%grad_data(:,k,n)*pos(l), i_coordinate, &
virial_loc(l,k), xloc(i), dcutoff_in=my_descriptor_data%x(i)%grad_covariance_cutoff(k,n)*pos(l))
enddo
enddo
enddo
enddo
endif
if(allocated(xloc)) deallocate(xloc)
enddo
if( has_local_property ) then
if( n_descriptors /= this%at(n_con)%N ) then
RAISE_ERROR("fit_data_from_xyz: local_propertyes found in configuration, but number of descriptors do not match &
& the number of atoms. Check your descriptors.",error)
endif
endif
if(allocated(force_loc)) deallocate(force_loc)
if(allocated(local_property_loc)) deallocate(local_property_loc)
if( has_hessian ) then
at = this%at(n_con)
call set_cutoff( at, my_cutoff )
do i_coordinate = 1, this%n_coordinate
allocate( grad_data(descriptor_dimensions(this%my_descriptor(i_coordinate))) )
do i = 1, n_hessian
if( .not. assign_pointer(this%at(n_con),trim(this%hessian_parameter_name)//i, hessian_eigenvector_i) ) &
call system_abort("fit_data_from_xyz: could not find "//i//"th of "//n_hessian//" hessian eigenvector.")
do j = -1, 1, 2
at%pos = this%at(n_con)%pos + j * this%hessian_delta * hessian_eigenvector_i
call calc_connect(at)
call calc(this%my_descriptor(i_coordinate),at,my_descriptor_data, &
do_descriptor=.true.,do_grad_descriptor=.true.)
!hessian_core(i) = hessian_core(i) + j * sum(f_hessian*hessian_eigenvector_i) / 2.0_dp / this%hessian_delta
allocate(xloc(size(my_descriptor_data%x)))
do k = 1, size(my_descriptor_data%x)
if( .not. my_descriptor_data%x(k)%has_data) cycle
xloc(k) = gp_addCoordinates(this%my_gp,my_descriptor_data%x(k)%data(:),i_coordinate, &
cutoff_in=my_descriptor_data%x(k)%covariance_cutoff,config_type=EXCLUDE_CONFIG_TYPE)
!cutoff_in=my_descriptor_data%x(k)%covariance_cutoff,config_type=n_config_type)
grad_data = 0.0_dp
grad_covariance_cutoff = 0.0_dp
do n = lbound(my_descriptor_data%x(k)%ii,1), ubound(my_descriptor_data%x(k)%ii,1)
if( .not. my_descriptor_data%x(k)%has_grad_data(n)) cycle
l = my_descriptor_data%x(k)%ii(n)
grad_data = grad_data + j * matmul(my_descriptor_data%x(k)%grad_data(:,:,n), hessian_eigenvector_i(:,l)) / 2.0_dp / this%hessian_delta
grad_covariance_cutoff = grad_covariance_cutoff + &
dot_product(my_descriptor_data%x(k)%grad_covariance_cutoff(:,n), hessian_eigenvector_i(:,l)) / 2.0_dp / this%hessian_delta
enddo
call gp_addCoordinateDerivatives(this%my_gp, grad_data, i_coordinate, &
hessian_loc(i), xloc(k), dcutoff_in=grad_covariance_cutoff)
enddo !k
deallocate(xloc)
enddo !j = -1, 1, 2
enddo ! i = 1, n_hessian
if(allocated(grad_data)) deallocate(grad_data)
enddo ! i_coordinate = 1, n_coordinate
endif !has_hessian
if(allocated(hessian_loc)) deallocate(hessian_loc)
if(allocated(hessian)) deallocate(hessian)
call finalise(my_descriptor_data)
enddo !n_frame
call print_title("Report on per-configuration/per-atom sigma (error parameter) settings")
call print("Number of per-configuration setting of energy_"//trim(this%sigma_parameter_name)//" found: "//n_energy_sigma)
call print("Number of per-configuration setting of force_"//trim(this%sigma_parameter_name)//" found: "//n_force_sigma)
call print("Number of per-configuration setting of virial_"//trim(this%sigma_parameter_name)//" found: "//n_virial_sigma)
call print("Number of per-configuration setting of hessian_"//trim(this%sigma_parameter_name)//" found: "//n_hessian_sigma)
call print("Number of per-atom setting of force_atom_"//trim(this%sigma_parameter_name)//" found: "//n_force_atom_sigma)
call print_title("End of report")
do i_coordinate = 1, this%n_coordinate
if( count( (/this%has_theta_file(i_coordinate), this%has_theta_uniform(i_coordinate), &
this%has_theta_fac(i_coordinate), this%has_zeta(i_coordinate) /) ) /= 1 ) then
call system_abort("fit_data_from_xyz: only one of theta_file, theta_uniform, theta_fac or zeta may be &
specified for each GAP.")
endif
if( this%covariance_type(i_coordinate) == COVARIANCE_DOT_PRODUCT ) then
if( .not. this%has_zeta(i_coordinate) ) call system_abort("fit_data_from_xyz: covariance type is DOT_PRODUCT but no zeta was specified.")
elseif( this%covariance_type(i_coordinate) == COVARIANCE_ARD_SE .or. this%covariance_type(i_coordinate) == COVARIANCE_PP ) then
if( count( (/this%has_theta_file(i_coordinate), this%has_theta_uniform(i_coordinate), this%has_theta_fac(i_coordinate) /) ) /= 1 ) then
call system_abort("fit_data_from_xyz: covariance type is ARD_SE or PP, so one of theta_file, theta_uniform of theta_fac must be specified")
endif
endif
if( this%has_theta_file(i_coordinate) ) then
allocate(theta_string_array(this%my_gp%coordinate(i_coordinate)%d))
allocate(theta(this%my_gp%coordinate(i_coordinate)%d))
call initialise(theta_inout,trim(this%theta_file(i_coordinate)))
read(theta_inout%unit,'(a)') theta_string
call split_string(theta_string,' :;','{}',theta_string_array,d,matching=.true.)
if(this%my_gp%coordinate(i_coordinate)%d /= d) call system_abort('File '//trim(this%theta_file(i_coordinate))//' does not contain the right number of hyperparameters')
do i = 1, d
theta(i) = string_to_real(trim(theta_string_array(i)))
enddo
call gp_setTheta(this%my_gp,i_coordinate,theta=theta)
deallocate(theta_string_array)
deallocate(theta)
call finalise(theta_inout)
elseif(this%has_theta_uniform(i_coordinate)) then
allocate(theta(this%my_gp%coordinate(i_coordinate)%d))
theta = this%theta_uniform(i_coordinate)
call gp_setTheta(this%my_gp,i_coordinate,theta=theta)
deallocate(theta)
elseif(this%has_theta_fac(i_coordinate)) then
allocate(theta_string_array(this%my_gp%coordinate(i_coordinate)%d))
allocate(theta_fac(this%my_gp%coordinate(i_coordinate)%d))
call split_string(trim(this%theta_fac_string(i_coordinate))," :;",'{}',theta_string_array,n_theta_fac,matching=.true.)
if(n_theta_fac == 1) then
theta_fac = string_to_real(theta_string_array(1))
elseif(n_theta_fac == this%my_gp%coordinate(i_coordinate)%d) then
do i = 1, this%my_gp%coordinate(i_coordinate)%d
theta_fac(i) = string_to_real(theta_string_array(i))
enddo
else
call system_abort("theta_fac can only contain one value or as many as dimensions the descriptor is")
endif
call gp_setThetaFactor(this%my_gp,i_coordinate,theta_fac,useSparseX=.false.)
deallocate(theta_fac)
deallocate(theta_string_array)
elseif( this%has_zeta(i_coordinate) ) then
call gp_setTheta(this%my_gp,i_coordinate,zeta=this%zeta(i_coordinate))
endif
enddo
if( this%do_core ) call Finalise(this%core_pot)
call gp_sparsify(this%my_gp,n_sparseX=this%config_type_n_sparseX,default_all=(this%n_sparseX/=0), &
sparseMethod=this%sparse_method, sparse_file=this%sparse_file, &
use_actual_gpcov=this%sparse_use_actual_gpcov, print_sparse_index = this%print_sparse_index, &
unique_hash_tolerance=this%unique_hash_tolerance, unique_descriptor_tolerance=this%unique_descriptor_tolerance)
end subroutine fit_data_from_xyz
subroutine e0_from_xyz(this)
type(gap_fit), intent(inout) :: this
integer :: n_con, n_ener, i, my_n_neighbours
logical :: has_ener
real(dp) :: ener, ener_core
logical, dimension(total_elements) :: found_Z, found_isolated
if( this%do_core ) call Initialise(this%core_pot, this%core_ip_args, param_str=string(this%quip_string))
n_ener = 0
this%e0 = 0.0_dp
found_isolated = .false.
found_Z = .false.
do n_con = 1, this%n_frame
has_ener = get_value(this%at(n_con)%params,trim(this%energy_parameter_name),ener)
found_Z(this%at(n_con)%Z) = .true.
if( has_ener ) then
ener_core = 0.0_dp
if( this%do_core ) then
if( this%at(n_con)%cutoff < cutoff(this%core_pot) ) then
call set_cutoff(this%at(n_con), cutoff(this%core_pot))
call calc_connect(this%at(n_con))
endif
call calc(this%core_pot,this%at(n_con),energy=ener_core)
endif
select case(this%e0_method)
case(E0_ISOLATED)
if( this%at(n_con)%N == 1 ) then
if( this%at(n_con)%cutoff < this%max_cutoff ) then
call set_cutoff(this%at(n_con), this%max_cutoff)
endif
call calc_connect(this%at(n_con))
if( n_neighbours(this%at(n_con),1,max_dist = this%max_cutoff) == 0 ) then
if( found_isolated(this%at(n_con)%Z(1)) ) then
call system_abort("Found more than one isolated atom configuration, which may be ambiguous.")
endif
this%e0(this%at(n_con)%Z(1)) = ener - ener_core
found_isolated(this%at(n_con)%Z(1)) = .true.
endif
endif
case(E0_AVERAGE)
this%e0 = this%e0 + (ener-ener_core) / this%at(n_con)%N
case default
call system_abort("Unknown e0_method")
endselect
n_ener = n_ener + 1
endif
enddo
select case(this%e0_method)
case(E0_ISOLATED)
if( .not. all(found_isolated .eqv. found_Z) ) then
do i = 1, size(found_Z)
if( found_Z(i) .and. .not. found_isolated(i) ) then
call print("Atom species "//i//" present in teaching XYZ, but not found corresponding isolated &
representative")
endif
enddo
call system_abort("Determination of e0 was requested to be based on isolated atom energies, but not all &
& atom types present in the XYZ had an isolated representative.")
endif
case(E0_AVERAGE)
if( n_ener > 0 ) then
this%e0 = this%e0 / n_ener
else
this%e0 = 0.0_dp
endif
case default
call system_abort("Unknown e0_method")
endselect
if( this%do_core ) call Finalise(this%core_pot)
endsubroutine e0_from_xyz
subroutine w_Z_from_xyz(this)
type(gap_fit), intent(inout) :: this
type(cinoutput) :: xyzfile
type(atoms) :: at
call initialise(xyzfile,this%at_file)
call read(xyzfile,at,frame=0)
!call get_weights(at,this%w_Z)
call finalise(at)
call finalise(xyzfile)
end subroutine w_Z_from_xyz
subroutine gap_fit_print_xml(this,filename,sparseX_separate_file)
use iso_c_binding, only : C_NULL_CHAR
type(gap_fit), intent(in) :: this
character(len=*), intent(in) :: filename
logical, intent(in), optional :: sparseX_separate_file
type(xmlf_t) :: xf
!type(extendable_str) :: gap_string
!type(inoutput) :: gp_inout
character(len=STRING_LENGTH) :: gp_tmp_file, gp_label
integer :: i
integer, dimension(8) :: values
logical :: my_sparseX_separate_file
call date_and_time(values=values)
! Get totally unique label for GAP. This will be used at various places.
write(gp_label,'("GAP_"7(i0,"_")i0)') values
! Unique temporary file
gp_tmp_file = 'tmp_'//trim(gp_label)//'.xml'
! Print GAP part of the potential into the temporary file.
call xml_OpenFile(gp_tmp_file,xf,addDecl=.false.)
call xml_NewElement(xf,"GAP_params")
call xml_AddAttribute(xf,"label",trim(gp_label))
call xml_AddAttribute(xf,"gap_version",""//gap_version)
call xml_NewElement(xf,"GAP_data")
call xml_AddAttribute(xf,"do_core",""//this%do_core)
do i = 1, size(this%e0)
call xml_NewElement(xf,"e0")
call xml_AddAttribute(xf,"Z",""//i)
call xml_AddAttribute(xf,"value",""// (this%e0(i)+this%local_property0(i) ))
call xml_EndElement(xf,"e0")
enddo
call xml_EndElement(xf,"GAP_data")
my_sparseX_separate_file = optional_default(.false., sparseX_separate_file)
! Print GP bit of the potential
if (my_sparseX_separate_file) then
call gp_printXML(this%gp_sp,xf,label=gp_label,sparseX_base_filename=trim(filename)//".sparseX")
else
call gp_printXML(this%gp_sp,xf,label=gp_label)
endif
! Print the command line used for the fitting
if(len(trim(this%command_line))> 0 ) then
call xml_NewElement(xf,"command_line")
call xml_AddCharacters(xf,trim(this%command_line),parsed=.false.)
call xml_EndElement(xf,"command_line")
endif
if(this%do_copy_at_file) then
! Print the fitting configurations used for this potential.
if(len(trim(this%at_file)) > 0 ) call file_print_xml(this%at_file,xf)
endif
call xml_EndElement(xf,"GAP_params")
call xml_Close(xf)
!! Now read back into an extendable string what we have just printed out.
!call read(gap_string, trim(gp_tmp_file), keep_lf=.true.)
!! Initialise the final file
!call initialise(gp_inout,trim(filename),action=OUTPUT)
! Open a unique root element for the xml
!call print('<'//trim(gp_label)//'>',file=gp_inout)
!!call system_command('echo "<'//trim(gp_label)//'>" >>'//trim(filename))
call fwrite_line_to_file(trim(filename)//C_NULL_CHAR,'<'//trim(gp_label)//'>'//C_NULL_CHAR,'w'//C_NULL_CHAR)
if(this%do_core) then
! Create the sum potential xml entry (by hand)
!call print('<Potential label="'//trim(gp_label)//'" init_args="Sum init_args_pot1={'//trim(this%ip_args)//'} init_args_pot2={IP GAP label='//trim(gp_label)//'}"/>',file=gp_inout)
!call system_command('echo "<Potential label=\"'//trim(gp_label)//'\" init_args=\"Sum init_args_pot1={'//trim(this%ip_args)//'} init_args_pot2={IP GAP label='//trim(gp_label)//'}\"/>" >>'//trim(filename))
call fwrite_line_to_file(trim(filename)//C_NULL_CHAR, &
'<Potential label="'//trim(gp_label)//'" init_args="Sum init_args_pot1={'//trim(this%core_ip_args)//'} init_args_pot2={IP GAP label='//trim(gp_label)//'}"/>'//C_NULL_CHAR, &
'a'//C_NULL_CHAR)
! Now add the core potential that was used.
!call print(string(this%quip_string),file=gp_inout)
!call system_command('echo "'//string(this%quip_string)//' >>'//trim(filename))
call fappend_file_to_file(trim(filename)//C_NULL_CHAR,trim(this%core_param_file)//C_NULL_CHAR)
else
call fwrite_line_to_file(trim(filename)//C_NULL_CHAR, &
'<Potential label="'//trim(gp_label)//'" init_args="IP GAP label='//trim(gp_label)//'"/>'//C_NULL_CHAR,'a'//C_NULL_CHAR)
endif
! Add the GAP potential
!call print(string(gap_string),file=gp_inout)
!call system_command('cat '//trim(gp_tmp_file)//' >>'//trim(filename))
call fappend_file_to_file(trim(filename)//C_NULL_CHAR,trim(gp_tmp_file)//C_NULL_CHAR)
! Close the root element
!call print('</'//trim(gp_label)//'>',file=gp_inout)
!call system_command('echo "</'//trim(gp_label)//'>" >>'//trim(filename))
call fwrite_line_to_file(trim(filename)//C_NULL_CHAR,'</'//trim(gp_label)//'>'//C_NULL_CHAR,'a'//C_NULL_CHAR)
!call finalise(gp_inout)
!call finalise(gap_string)
! Delete the temporary file
!call system_command('rm -f '//trim(gp_tmp_file))
call frm_file(trim(gp_tmp_file)//C_NULL_CHAR)
endsubroutine gap_fit_print_xml
subroutine file_print_xml(this,xf)
character(len=*), intent(in) :: this
type(xmlf_t), intent(inout) :: xf
type(inoutput) :: atfile
character(len=10240) :: line
integer :: iostat
call initialise(atfile,trim(this))
call xml_NewElement(xf,"XYZ_data")
call xml_AddNewLine(xf)
do
read(atfile%unit,'(a)',iostat=iostat) line
if(iostat < 0) then
exit
elseif(iostat > 0) then
call system_abort('file_print_xml: unkown error ('//iostat//') while reading '//trim(this))
endif
call xml_AddCharacters(xf,trim(line),parsed=.false.)
call xml_AddNewLine(xf)
enddo
call xml_EndElement(xf,"XYZ_data")
call finalise(atfile)
endsubroutine file_print_xml
! subroutine print_sparse(this)
! type(gap_fit), intent(in) :: this
! type(cinoutput) :: xyzfile, xyzfile_out
! type(atoms) :: at, at_out
!
! integer :: li, ui, n_con
! logical, dimension(:), allocatable :: x
! logical, dimension(:), pointer :: sparse
!
! if(this%do_mark_sparse_atoms) then
!
! allocate(x(this%n_descriptors))
! x = .false.
! x(this%r) = .true.
!
! call initialise(xyzfile,this%at_file)
! call initialise(xyzfile_out,this%mark_sparse_atoms,action=OUTPUT)
!
! li = 0
! ui = 0
! do n_con = 1, xyzfile%n_frame
! call read(xyzfile,at,frame=n_con-1)
! at_out = at
!
! call add_property(at_out,'sparse',.false.,ptr=sparse)
!
! li = ui + 1
! ui = ui + at%N
! if(any( x(li:ui) )) sparse(find_indices(x(li:ui))) = .true.
!
! call write(at_out,xyzfile_out,properties="species:pos:sparse")
! enddo
! call finalise(xyzfile)
! call finalise(xyzfile_out)
! deallocate(x)
!
! endif
!
! endsubroutine print_sparse
subroutine parse_config_type_sigma(this)
type(gap_fit), intent(inout) :: this
character(len=STRING_LENGTH), dimension(99) :: config_type_sigma_fields
integer :: config_type_sigma_num_fields, i_default, i, n_config_type
if( this%has_config_type_sigma ) then
call split_string(this%config_type_sigma_string,' :;','{}',config_type_sigma_fields,config_type_sigma_num_fields,matching=.true.)
n_config_type = config_type_sigma_num_fields / 5
! find "default" if present
i_default = 0
do i = 1, config_type_sigma_num_fields, 5
if( trim(config_type_sigma_fields(i)) == "default" ) i_default = i
enddo
if( i_default == 0 ) then
! no default present in the string, we add it, and it'll be the last one
n_config_type = n_config_type + 1
i_default = n_config_type
config_type_sigma_fields(config_type_sigma_num_fields+1) = "default"
config_type_sigma_fields(config_type_sigma_num_fields+2) = ""//this%default_sigma(1)
config_type_sigma_fields(config_type_sigma_num_fields+3) = ""//this%default_sigma(2)
config_type_sigma_fields(config_type_sigma_num_fields+4) = ""//this%default_sigma(3)
config_type_sigma_fields(config_type_sigma_num_fields+5) = ""//this%default_sigma(4)
config_type_sigma_num_fields = config_type_sigma_num_fields + 5
endif
allocate(this%config_type(n_config_type))
allocate(this%sigma(4,n_config_type))
do i = 1, n_config_type
this%config_type(i) = trim(config_type_sigma_fields(5*(i-1)+1))
this%sigma(1,i) = string_to_real(config_type_sigma_fields(5*(i-1)+2))
this%sigma(2,i) = string_to_real(config_type_sigma_fields(5*(i-1)+3))
this%sigma(3,i) = string_to_real(config_type_sigma_fields(5*(i-1)+4))
this%sigma(4,i) = string_to_real(config_type_sigma_fields(5*(i-1)+5))
enddo
call print('Sparse points and target errors per pre-defined types of configurations')
do i = 1, n_config_type
call print(""//trim(this%config_type(i))//" "//this%sigma(:,i))
enddo
else
allocate(this%config_type(1))
allocate(this%sigma(4,1))
this%config_type(1)= "default"
this%sigma(:,1) = this%default_sigma
endif
endsubroutine parse_config_type_sigma
subroutine parse_config_type_n_sparseX(this)
type(gap_fit), intent(inout) :: this
integer :: i, j, i_default, i_coordinate, i_config_type, config_type_n_sparseX_num_fields, n_config_type, new_config_types
character(len=STRING_LENGTH), dimension(99) :: config_type_n_sparseX_fields
logical :: config_type_present
if( .not. allocated(this%config_type) ) call system_abort('config_type not allocated, call parse_config_type_sigma first')
do i = 1, size(this%config_type)
if( trim(this%config_type(i)) == "default" ) i_default = i
enddo
! Check first if we have more new config types than we had from config_type_sigma
do i_coordinate = 1, this%n_coordinate
if( this%n_sparseX(i_coordinate) == 0 .and. len_trim(this%config_type_n_sparseX_string(i_coordinate)) > 0) then
call split_string(this%config_type_n_sparseX_string(i_coordinate),' :;','{}',config_type_n_sparseX_fields,config_type_n_sparseX_num_fields,matching=.true.)
if( mod(config_type_n_sparseX_num_fields,2) /= 0 ) then
call system_abort("parse_config_type_n_sparseX: config_type_n_sparseX could not be parsed correctly, key/value pairs must always be present")
endif
n_config_type = size(this%config_type)
new_config_types = 0 ! Assume there are no new config_types
do j = 1, config_type_n_sparseX_num_fields, 2 ! loop over config_types in the descriptor string
config_type_present = .false.
do i = 1, n_config_type ! loop over config_types previously set
if( trim(this%config_type(i)) == trim(config_type_n_sparseX_fields(j)) ) config_type_present = .true. ! Found config_type among old ones
enddo
if(.not.config_type_present) new_config_types = new_config_types + 1 ! Increment as it's a genuine new config_type
enddo
if( new_config_types > 0 ) then
call reallocate(this%config_type, n_config_type + new_config_types, copy=.true.)
call reallocate(this%sigma,4,n_config_type + new_config_types, copy=.true.)
i_config_type = n_config_type
do j = 1, config_type_n_sparseX_num_fields, 2 ! loop over config_types in the descriptor string
config_type_present = .false.
do i = 1, n_config_type ! loop over config_types previously set
if( trim(this%config_type(i)) == trim(config_type_n_sparseX_fields(j)) ) config_type_present = .true. ! Found config_type among old ones
enddo
if(.not.config_type_present) then ! it's a genuine new config_type
i_config_type = i_config_type + 1
this%config_type(i_config_type) = trim(config_type_n_sparseX_fields(j))
this%sigma(:,i_config_type) = this%sigma(:,i_default)
endif
enddo
endif
elseif(this%n_sparseX(i_coordinate) > 0 .and. len_trim(this%config_type_n_sparseX_string(i_coordinate)) > 0 .and. len_trim(this%sparse_file(i_coordinate)) ==0 ) then
call system_abort('Confused: cannot specify both n_sparse and config_type_n_sparse')
elseif(this%n_sparseX(i_coordinate) == 0 .and. len_trim(this%config_type_n_sparseX_string(i_coordinate)) == 0 .and. len_trim(this%sparse_file(i_coordinate)) == 0) then
call system_abort('Confused: either n_sparse or config_type_n_sparse has to be specified')
endif
enddo
n_config_type = size(this%config_type)
allocate(this%config_type_n_sparseX(n_config_type,this%n_coordinate))
this%config_type_n_sparseX = 0
do i_coordinate = 1, this%n_coordinate
if( this%n_sparseX(i_coordinate) == 0 .and. len_trim(this%config_type_n_sparseX_string(i_coordinate)) > 0) then
call split_string(this%config_type_n_sparseX_string(i_coordinate),' :;','{}',config_type_n_sparseX_fields,config_type_n_sparseX_num_fields,matching=.true.)
do j = 1, config_type_n_sparseX_num_fields, 2 ! loop over config_types in the descriptor string
do i = 1, n_config_type ! loop over config_types previously set
if( trim(this%config_type(i)) == trim(config_type_n_sparseX_fields(j)) ) &
this%config_type_n_sparseX(i,i_coordinate) = string_to_int( config_type_n_sparseX_fields(j+1) )
enddo
enddo
!this%n_sparseX(i_coordinate) = sum( this%config_type_n_sparseX(:,i_coordinate) )
elseif( this%n_sparseX(i_coordinate) > 0 .and. len_trim(this%config_type_n_sparseX_string(i_coordinate)) == 0) then
this%config_type_n_sparseX(i_default,i_coordinate) = this%n_sparseX(i_coordinate)
endif
enddo
endsubroutine parse_config_type_n_sparseX
subroutine get_species_xyz(this)
type(gap_fit), intent(inout) :: this
integer :: n_con, i
integer, dimension(total_elements) :: species_present
this%n_species = 0
species_present = 0
do n_con = 1, this%n_frame
do i = 1, this%at(n_con)%N
if( all(this%at(n_con)%Z(i) /= species_present) ) then
this%n_species = this%n_species + 1
species_present(this%n_species) = this%at(n_con)%Z(i)
endif
enddo
enddo
allocate(this%species_Z(this%n_species))
this%species_Z = species_present(1:this%n_species)
endsubroutine get_species_xyz
subroutine add_multispecies_gaps(this)
type(gap_fit), intent(inout) :: this
integer :: i_coordinate, i, j, n_gap_str, i_add_species
character(STRING_LENGTH), dimension(:), allocatable :: gap_str_i, new_gap_str
! temporary arrays
real(dp), dimension(:), allocatable :: delta, f0, theta_uniform, zeta, unique_hash_tolerance, unique_descriptor_tolerance
integer, dimension(:), allocatable :: n_sparseX, sparse_method, covariance_type
character(len=STRING_LENGTH), dimension(:), allocatable :: theta_file, sparse_file, theta_fac_string, config_type_n_sparseX_string, print_sparse_index
logical, dimension(:), allocatable :: mark_sparse_atoms, has_theta_fac, has_theta_uniform, has_theta_file, has_zeta
n_gap_str = 0
do i_coordinate = 1, this%n_coordinate
if( this%add_species(i_coordinate) ) then
call print('Old GAP: {'//trim(this%gap_str(i_coordinate))//'}')
call descriptor_str_add_species(this%gap_str(i_coordinate),this%species_Z,gap_str_i)
call reallocate(new_gap_str, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(delta, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(f0, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(n_sparseX, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(config_type_n_sparseX_string, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(theta_fac_string, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(theta_uniform, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(theta_file, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(has_theta_fac, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(has_theta_uniform, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(has_theta_file, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(sparse_file, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(mark_sparse_atoms, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(sparse_method, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(covariance_type, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(zeta, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(has_zeta, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(print_sparse_index, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(unique_hash_tolerance, n_gap_str+size(gap_str_i),copy=.true.)
call reallocate(unique_descriptor_tolerance, n_gap_str+size(gap_str_i),copy=.true.)
do i = 1, size(gap_str_i)
i_add_species = index(gap_str_i(i),'add_species')
if(i_add_species /= 0) then
do j = i_add_species, len_trim(gap_str_i(i))
if( gap_str_i(i)(j:j) == " " ) exit
gap_str_i(i)(j:j) = " "
!gap_str_i(i)(i_add_species:i_add_species+len('add_species')-1) = ' '
enddo
endif
new_gap_str(i+n_gap_str) = trim(gap_str_i(i))
call print('New GAP: {'//trim(gap_str_i(i))//'}')
delta(i+n_gap_str) = this%delta(i_coordinate)
f0(i+n_gap_str) = this%f0(i_coordinate)
n_sparseX(i+n_gap_str) = this%n_sparseX(i_coordinate)
config_type_n_sparseX_string(i+n_gap_str) = this%config_type_n_sparseX_string(i_coordinate)
theta_fac_string(i+n_gap_str) = this%theta_fac_string(i_coordinate)
theta_uniform(i+n_gap_str) = this%theta_uniform(i_coordinate)
theta_file(i+n_gap_str) = this%theta_file(i_coordinate)
has_theta_fac(i+n_gap_str) = this%has_theta_fac(i_coordinate)
has_theta_uniform(i+n_gap_str) = this%has_theta_uniform(i_coordinate)
has_theta_file(i+n_gap_str) = this%has_theta_file(i_coordinate)
sparse_file(i+n_gap_str) = this%sparse_file(i_coordinate)
mark_sparse_atoms(i+n_gap_str) = this%mark_sparse_atoms(i_coordinate)
sparse_method(i+n_gap_str) = this%sparse_method(i_coordinate)
covariance_type(i+n_gap_str) = this%covariance_type(i_coordinate)
zeta(i+n_gap_str) = this%zeta(i_coordinate)
has_zeta(i+n_gap_str) = this%has_zeta(i_coordinate)
print_sparse_index(i+n_gap_str) = this%print_sparse_index(i_coordinate)
unique_hash_tolerance(i+n_gap_str) = this%unique_hash_tolerance(i_coordinate)
unique_descriptor_tolerance(i+n_gap_str) = this%unique_descriptor_tolerance(i_coordinate)
enddo
n_gap_str = n_gap_str + size(gap_str_i)
deallocate(gap_str_i)
else
n_gap_str = n_gap_str + 1
call reallocate(new_gap_str, n_gap_str,copy=.true.)
call reallocate(delta, n_gap_str,copy=.true.)
call reallocate(f0, n_gap_str,copy=.true.)
call reallocate(n_sparseX, n_gap_str,copy=.true.)
call reallocate(config_type_n_sparseX_string, n_gap_str,copy=.true.)
call reallocate(theta_fac_string, n_gap_str,copy=.true.)
call reallocate(theta_uniform, n_gap_str,copy=.true.)
call reallocate(theta_file, n_gap_str,copy=.true.)
call reallocate(has_theta_fac, n_gap_str,copy=.true.)
call reallocate(has_theta_uniform, n_gap_str,copy=.true.)
call reallocate(has_theta_file, n_gap_str,copy=.true.)
call reallocate(sparse_file, n_gap_str,copy=.true.)
call reallocate(mark_sparse_atoms, n_gap_str,copy=.true.)
call reallocate(sparse_method, n_gap_str,copy=.true.)
call reallocate(covariance_type, n_gap_str,copy=.true.)
call reallocate(zeta, n_gap_str,copy=.true.)
call reallocate(has_zeta, n_gap_str,copy=.true.)
call reallocate(print_sparse_index, n_gap_str,copy=.true.)
call reallocate(unique_hash_tolerance, n_gap_str,copy=.true.)
call reallocate(unique_descriptor_tolerance, n_gap_str,copy=.true.)
new_gap_str(n_gap_str) = trim(this%gap_str(i_coordinate))
delta(n_gap_str) = this%delta(i_coordinate)
f0(n_gap_str) = this%f0(i_coordinate)
n_sparseX(n_gap_str) = this%n_sparseX(i_coordinate)
config_type_n_sparseX_string(n_gap_str) = this%config_type_n_sparseX_string(i_coordinate)
theta_fac_string(n_gap_str) = this%theta_fac_string(i_coordinate)
theta_uniform(n_gap_str) = this%theta_uniform(i_coordinate)
theta_file(n_gap_str) = this%theta_file(i_coordinate)
has_theta_fac(n_gap_str) = this%has_theta_fac(i_coordinate)
has_theta_uniform(n_gap_str) = this%has_theta_uniform(i_coordinate)
has_theta_file(n_gap_str) = this%has_theta_file(i_coordinate)
sparse_file(n_gap_str) = this%sparse_file(i_coordinate)
mark_sparse_atoms(n_gap_str) = this%mark_sparse_atoms(i_coordinate)
sparse_method(n_gap_str) = this%sparse_method(i_coordinate)
covariance_type(n_gap_str) = this%covariance_type(i_coordinate)
zeta(n_gap_str) = this%zeta(i_coordinate)
has_zeta(n_gap_str) = this%has_zeta(i_coordinate)
print_sparse_index(n_gap_str) = this%print_sparse_index(i_coordinate)
unique_hash_tolerance(n_gap_str) = this%unique_hash_tolerance(i_coordinate)
unique_descriptor_tolerance(n_gap_str) = this%unique_descriptor_tolerance(i_coordinate)
call print('Unchanged GAP: {'//trim(this%gap_str(i_coordinate))//'}')
endif
enddo
call reallocate(this%delta, n_gap_str)
call reallocate(this%f0, n_gap_str)
call reallocate(this%n_sparseX, n_gap_str)
call reallocate(this%config_type_n_sparseX_string, n_gap_str)
call reallocate(this%theta_fac_string, n_gap_str)
call reallocate(this%theta_uniform, n_gap_str)
call reallocate(this%theta_file, n_gap_str)
call reallocate(this%has_theta_fac, n_gap_str)
call reallocate(this%has_theta_uniform, n_gap_str)
call reallocate(this%has_theta_file, n_gap_str)
call reallocate(this%sparse_file, n_gap_str)
call reallocate(this%mark_sparse_atoms, n_gap_str)
call reallocate(this%sparse_method, n_gap_str)
call reallocate(this%covariance_type, n_gap_str)
call reallocate(this%zeta, n_gap_str)
call reallocate(this%has_zeta, n_gap_str)
call reallocate(this%print_sparse_index, n_gap_str)
call reallocate(this%unique_hash_tolerance, n_gap_str)
call reallocate(this%unique_descriptor_tolerance, n_gap_str)
this%gap_str(1:n_gap_str) = new_gap_str
this%delta = delta
this%f0 = f0
this%n_sparseX = n_sparseX
this%config_type_n_sparseX_string = config_type_n_sparseX_string
this%theta_fac_string = theta_fac_string
this%theta_uniform = theta_uniform
this%theta_file = theta_file
this%has_theta_fac = has_theta_fac
this%has_theta_uniform = has_theta_uniform
this%has_theta_file = has_theta_file
this%sparse_file = sparse_file
this%mark_sparse_atoms = mark_sparse_atoms
this%sparse_method = sparse_method
this%covariance_type = covariance_type
this%zeta = zeta
this%has_zeta = has_zeta
this%print_sparse_index = print_sparse_index
this%unique_hash_tolerance = unique_hash_tolerance
this%unique_descriptor_tolerance = unique_descriptor_tolerance
this%n_coordinate = n_gap_str
if(allocated(delta)) deallocate(delta)
if(allocated(f0)) deallocate(f0)
if(allocated(n_sparseX)) deallocate(n_sparseX)
if(allocated(config_type_n_sparseX_string)) deallocate(config_type_n_sparseX_string)
if(allocated(theta_fac_string)) deallocate(theta_fac_string)
if(allocated(theta_uniform)) deallocate(theta_uniform)
if(allocated(theta_file)) deallocate(theta_file)
if(allocated(has_theta_fac)) deallocate(has_theta_fac)
if(allocated(has_theta_uniform)) deallocate(has_theta_uniform)
if(allocated(has_theta_file)) deallocate(has_theta_file)
if(allocated(sparse_file)) deallocate(sparse_file)
if(allocated(mark_sparse_atoms)) deallocate(mark_sparse_atoms)
if(allocated(sparse_method)) deallocate(sparse_method)
if(allocated(covariance_type)) deallocate(covariance_type)
if(allocated(zeta)) deallocate(zeta)
if(allocated(has_zeta)) deallocate(has_zeta)
if(allocated(print_sparse_index)) deallocate(print_sparse_index)
if(allocated(unique_hash_tolerance)) deallocate(unique_hash_tolerance)
if(allocated(unique_descriptor_tolerance)) deallocate(unique_descriptor_tolerance)
endsubroutine add_multispecies_gaps
subroutine add_template_string(this)
type(gap_fit), intent(inout) :: this
character(len=STRING_LENGTH) :: template_string=' '
character(len=STRING_LENGTH),dimension(:), allocatable :: lines_array
type(inoutput) :: tempfile
integer :: i,n_lines,total_length=0
if( this%has_template_file ) then
call print("adding template string, reading from file "//trim(this%template_file))
call initialise(tempfile,trim(this%template_file))
call read_file(tempfile,lines_array,n_lines)
do i=1,n_lines-1
template_string=trim(template_string)//"{"//trim(lines_array(i))//"};"
total_length = total_length + len_trim(lines_array(i))
end do
template_string=trim(template_string)//"{"//trim(lines_array(n_lines))//"}"
total_length = total_length + len_trim(lines_array(n_lines))
if (total_length .ge. STRING_LENGTH) call system_abort("Template atoms object exceeds maximum string size")
do i=1,len_trim(template_string)
if(template_string(i:i)==' ') then
template_string(i:i)='%'
end if
end do
!call print(template_string)
do i=1,this%n_coordinate
this%gap_str(i) = trim(this%gap_str(i))//" atoms_template_string={"//trim(template_string)//"}"
end do
endif
end subroutine add_template_string
end module gap_fit_module
3rdparty/gap/gapversion deleted 100755 → 0
View file @ 23bcf97f
#!/bin/bash
#
# Determine the gap version from the git reported last changed date
# of file (or the file GAP_VERSION if it exists). The gapversion
# is the UNIX timestap of the most recently changed file.
GAP_FILES="descriptors.f95 descriptors_wrapper.f95 \
make_permutations_v2.f95 gp_predict.f95 \
clustering.f95 gp_fit.f95 \
gap_fit_module.f95 gap_fit.f95"
GAP_ROOT=$(dirname "$0")
# By default only show the date
VERBOSE=false
# Options for interactive use
while getopts 'v' opt; do
case $opt in
v)
VERBOSE=true
;;
\?)
echo "Usage: $0"
echo "'-v' to output name of newest file to stderr"
exit 1
;;
esac
done
function git_date
{
if [ -e "$1" ]
then
DATE=$(git log -n1 --format="%at" "$1")
[[ $? -eq 0 && ! -z $DATE ]] && echo "$DATE" || echo 0
else
echo 0
fi
}
if [ -s "${GAP_ROOT}/GAP_VERSION" ]; then
echo -ne "$(cat "${GAP_ROOT}/GAP_VERSION")"
exit 0
elif [ -d "${GAP_ROOT}/.git" ] || [ -s "${GAP_ROOT}/.git" ]; then
# Sort everything as "DATE filename" list and take the last one
GAP_VERSION=$(
for I in $GAP_FILES
do
if [ -d "${GAP_ROOT}/$(dirname "$I")" ]
then
cd "${GAP_ROOT}/$(dirname "$I")"
echo $(git_date $(basename "$I")) "$I"
cd - >/dev/null
fi
done | sort -n | tail -1 )
# filename to stderr if requested;
# split string with parameter expansion
if [ $VERBOSE == "true" ]; then
echo "${GAP_VERSION##* }" >&2
fi
echo "${GAP_VERSION%% *}"
else
echo "0"
exit 0
fi
3rdparty/gap/gp_fit.f95 deleted 100644 → 0
View file @ 23bcf97f
! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! HND X
! HND X libAtoms+QUIP: atomistic simulation library
! HND X
! HND X Portions of this code were written by
! HND X Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! HND X Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! HND X
! HND X Copyright 2006-2010.
! HND X
! HND X Not for distribution
! HND X
! HND X Portions of this code were written by Noam Bernstein as part of
! HND X his employment for the U.S. Government, and are not subject
! HND X to copyright in the USA.
! HND X
! HND X When using this software, please cite the following reference:
! HND X
! HND X http://www.libatoms.org
! HND X
! HND X Additional contributions by
! HND X Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! HND X
! HND XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!X
!X Gaussian Process module
!X
!% Module for general GP function interpolations.
!% A gp object contains the training set (fitting points and function values),
!% important temporary matrices, vectors and parameters.
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#include "error.inc"
module gp_fit_module
use iso_c_binding, only : C_NULL_CHAR
! use libatoms_module
use error_module
use system_module
use extendable_str_module
use linearalgebra_module
use dictionary_module, only : STRING_LENGTH
use gp_predict_module
use clustering_module
implicit none
private
integer, parameter, public :: EXCLUDE_CONFIG_TYPE = -10
interface gp_sparsify
module procedure gpFull_sparsify_array_config_type
endinterface gp_sparsify
public :: gp_sparsify
contains
subroutine gpCoordinates_sparsify_config_type(this, n_sparseX, default_all, sparseMethod, sparse_file, use_actual_gpcov, print_sparse_index, &
unique_hash_tolerance, unique_descriptor_tolerance, error)
type(gpCoordinates), intent(inout) :: this
integer, dimension(:), intent(in) :: n_sparseX
logical, intent(in) :: default_all
integer, optional, intent(in) :: sparseMethod
character(len=STRING_LENGTH), intent(in), optional :: sparse_file, print_sparse_index
logical, intent(in), optional :: use_actual_gpcov
real(dp), intent(in), optional :: unique_descriptor_tolerance, unique_hash_tolerance
integer, optional, intent(out) :: error
integer :: my_sparseMethod, i, j, li, ui, i_config_type, n_config_type, d, n_x, n_sparse_file
integer, dimension(:), allocatable :: config_type_index, sparseX_index, my_n_sparseX
real(dp), dimension(:,:), allocatable :: x, sparseX_array
real(dp), dimension(:), allocatable :: x_hash
real(dp), pointer, dimension(:,:) :: dm
integer, dimension(:), allocatable :: x_index
logical, dimension(:), allocatable :: x_unique
character(len=STRING_LENGTH) :: my_sparse_file
type(Inoutput) :: inout_sparse_index
logical :: exist_sparse_file
real(dp) :: my_unique_hash_tolerance, my_unique_descriptor_tolerance
INIT_ERROR(error)
my_sparseMethod = optional_default(GP_SPARSE_RANDOM,sparseMethod)
my_sparse_file = optional_default("",sparse_file)
my_unique_hash_tolerance = optional_default(1.0e-10_dp,unique_hash_tolerance)
my_unique_descriptor_tolerance = optional_default(1.0e-10_dp,unique_descriptor_tolerance)
if( .not. this%initialised ) then
RAISE_ERROR('gpCoordinates_sparsify: : object not initialised',error)
endif
d = size(this%x,1)
!n_x = size(this%x,2)
n_x = count(EXCLUDE_CONFIG_TYPE /= this%config_type)
allocate(my_n_sparseX(size(n_sparseX)))
my_n_sparseX = 0
! Remove duplicates
allocate(x_hash(n_x))
allocate(x_index(n_x))
! Compute 1-norm hash on all descriptors that we want to include, and the mapping to the full vector
j = 0
do i = 1, size(this%x,2)
if( this%config_type(i) /= EXCLUDE_CONFIG_TYPE ) then
j = j + 1
x_hash(j) = sum(abs(this%x(:,i)))
x_index(j) = i
endif
enddo
! Sort hashes
call heap_sort(x_hash,i_data=x_index)
! Compare neighbouring hashes. If they're within tolerance, compare the corresponding descriptors using the eucledian norm.
! Update the config type if they're equivalent.
do j = 2, n_x
if( abs( x_hash(j-1) - x_hash(j) ) < my_unique_hash_tolerance ) then
!if ( sum( ( this%x(:,x_index(j))-this%x(:,x_index(j-1)) )**2 ) < my_unique_descriptor_tolerance ) then
if ( maxval( abs( this%x(:,x_index(j))-this%x(:,x_index(j-1)) ) ) < my_unique_descriptor_tolerance ) then
this%config_type(x_index(j-1)) = EXCLUDE_CONFIG_TYPE
endif
endif
enddo
deallocate(x_hash)
deallocate(x_index)
! need to update n_x once duplicates are removed
n_x = count(EXCLUDE_CONFIG_TYPE /= this%config_type)
if(my_sparseMethod == GP_SPARSE_UNIQ) then
RAISE_ERROR('gpCoordinates_sparsify: UNIQ is no longer in use, please use NONE instead.',error)
! allocate(x(d,n_x))
! allocate(x_index(n_x))
! allocate(x_unique(n_x))
! !x = this%x
! !x_index = (/(i,i=1,n_x)/)
! j = 0
! do i = 1, size(this%x,2)
! if( this%config_type(i) /= EXCLUDE_CONFIG_TYPE ) then
! j = j + 1
! x(:,j) = this%x(:,i)
! x_index(j) = i
! endif
! enddo
! call heap_sort(x,i_data=x_index)
! call uniq(x,unique=x_unique)
! this%n_sparseX = count(x_unique)
! call print('UNIQ type sparsification specified. The number of sparse points was changed to '//this%n_sparseX//' from '//n_sparseX//'.')
elseif(my_sparseMethod == GP_SPARSE_NONE) then
allocate(x_index(n_x))
j = 0
do i = 1, size(this%x,2)
if( this%config_type(i) /= EXCLUDE_CONFIG_TYPE ) then
j = j + 1
x_index(j) = i
endif
enddo
this%n_sparseX = n_x
call print('NONE type sparsification specified. The number of sparse points was changed to '//this%n_sparseX//' from '//n_sparseX//'.')
elseif(my_sparseMethod == GP_SPARSE_FILE .or. my_sparseMethod == GP_SPARSE_INDEX_FILE) then
inquire(file=trim(my_sparse_file),exist=exist_sparse_file)
if(.not. exist_sparse_file) then
RAISE_ERROR('gpCoordinates_sparsify: '//trim(my_sparse_file)//' does not exist',error)
endif
call fwc_l(trim(my_sparse_file)//C_NULL_CHAR,n_sparse_file)
if( my_sparseMethod == GP_SPARSE_FILE ) then
if(mod(n_sparse_file,d+1) /= 0) then
RAISE_ERROR('gpCoordinates_sparsify: file '//trim(my_sparse_file)//' contains '//n_sparse_file//" lines, not conforming with descriptor size "//d,error)
endif
this%n_sparseX = n_sparse_file / ( d + 1 )
else
this%n_sparseX = n_sparse_file
endif
else
do i_config_type = 1, size(n_sparseX)
if(default_all) then
if( n_x < sum(n_sparseX) ) then
call print_warning('gpCoordinates_sparsify: number of data points ('//n_x//') less than the number of sparse points ('//sum(n_sparseX)//'), &
number of sparse points changed to '//n_x)
call print_warning('gpCoordinates_sparsify: affected descriptor : '//this%descriptor_str)
my_n_sparseX(1) = n_x
else
my_n_sparseX(1) = sum(n_sparseX)
endif
else
if( n_sparseX(i_config_type) == 0 ) cycle
n_config_type = count(i_config_type == this%config_type)
if( n_config_type < n_sparseX(i_config_type) ) then
call print_warning('gpCoordinates_sparsify: number of data points ('//n_config_type//') less than the number of sparse points ('//n_sparseX(i_config_type)//'), &
number of sparse points changed to '//n_config_type)
call print_warning('gpCoordinates_sparsify: affected descriptor : '//this%descriptor_str)
my_n_sparseX(i_config_type) = n_config_type
else
my_n_sparseX(i_config_type) = n_sparseX(i_config_type)
endif
endif
if(default_all) exit
enddo
this%n_sparseX = sum(my_n_sparseX)
endif
call reallocate(this%sparseX, this%d,this%n_sparseX, zero = .true.)
call reallocate(this%sparseX_index, this%n_sparseX, zero = .true.)
call reallocate(this%map_sparseX_globalSparseX, this%n_sparseX, zero = .true.)
call reallocate(this%alpha, this%n_sparseX, zero = .true.)
call reallocate(this%sparseCutoff, this%n_sparseX, zero = .true.)
this%sparseCutoff = 1.0_dp
if( my_sparseMethod /= GP_SPARSE_FILE .and. my_sparseMethod /= GP_SPARSE_INDEX_FILE) then
ui = 0
do i_config_type = 1, size(my_n_sparseX)
if( my_sparseMethod == GP_SPARSE_NONE) exit
if(default_all) then
allocate(config_type_index(n_x), sparseX_index(this%n_sparseX))
!config_type_index = (/(i,i=1,n_x)/)
j = 0
do i = 1, size(this%x,2)
if( this%config_type(i) /= EXCLUDE_CONFIG_TYPE ) then
j = j + 1
config_type_index(j) = i
endif
enddo
li = 1
ui = this%n_sparseX
n_config_type = n_x
else
if( my_n_sparseX(i_config_type) == 0 ) cycle
n_config_type = count(i_config_type == this%config_type)
allocate(config_type_index(n_config_type),sparseX_index(my_n_sparseX(i_config_type)))
config_type_index = find(i_config_type == this%config_type)
li = ui + 1
ui = ui + my_n_sparseX(i_config_type)
endif
select case(my_sparseMethod)
case(GP_SPARSE_RANDOM)
call fill_random_integer(sparseX_index, n_config_type)
case(GP_SPARSE_PIVOT)
if(this%covariance_type == COVARIANCE_DOT_PRODUCT) then
call pivot(this%x(:,config_type_index), sparseX_index)
else
call pivot(this%x(:,config_type_index), sparseX_index, theta = this%theta)
endif
case(GP_SPARSE_CLUSTER)
if(use_actual_gpcov) then
call print('Started kernel distance matrix calculation')
dm => kernel_distance_matrix(this, config_type_index = config_type_index)
call print('Finished kernel distance matrix calculation')
endif
call print('Started kmedoids clustering')
if(use_actual_gpcov) then
call bisect_kmedoids(dm, my_n_sparseX(i_config_type), med = sparseX_index)
!call bisect_kmedoids(dm, my_n_sparseX(i_config_type), c = c, med = sparseX_index)
else
if(this%covariance_type == COVARIANCE_DOT_PRODUCT) then
call bisect_kmedoids(this%x(:,config_type_index), my_n_sparseX(i_config_type), med = sparseX_index, is_distance_matrix = .false.)
else
call bisect_kmedoids(this%x(:,config_type_index), my_n_sparseX(i_config_type), med = sparseX_index, theta = this%theta, is_distance_matrix = .false.)
endif
endif
call print('Finished kmedoids clustering')
if(use_actual_gpcov) deallocate(dm)
case(GP_SPARSE_UNIFORM)
call select_uniform(this%x(:,config_type_index), sparseX_index)
case(GP_SPARSE_KMEANS)
call print('Started kmeans clustering')
if(this%covariance_type == COVARIANCE_DOT_PRODUCT) then
call cluster_kmeans(this%x(:,config_type_index), sparseX_index)
else
call cluster_kmeans(this%x(:,config_type_index), sparseX_index, theta = this%theta)
endif
call print('Finished kmeans clustering')
case(GP_SPARSE_COVARIANCE)
call sparse_covariance(this,sparseX_index,config_type_index,use_actual_gpcov)
case(GP_SPARSE_FUZZY)
call print('Started fuzzy cmeans clustering')
if(this%covariance_type == COVARIANCE_DOT_PRODUCT) then
call cluster_fuzzy_cmeans(this%x(:,config_type_index), sparseX_index, fuzziness=2.0_dp)
else
call cluster_fuzzy_cmeans(this%x(:,config_type_index), sparseX_index, theta=this%theta,fuzziness=2.0_dp)
endif
call print('Finished fuzzy cmeans clustering')
case(GP_SPARSE_CUR_COVARIANCE)
call print("Started covariance matrix calculation")
dm => kernel_distance_matrix(this, config_type_index=config_type_index, covariance_only = .true.)
call print("Finished covariance matrix calculation")
call print("Started CUR decomposition")
call cur_decomposition(dm, sparseX_index)
call print("Finished CUR decomposition")
deallocate(dm)
case(GP_SPARSE_CUR_POINTS)
call print("Started CUR decomposition")
call cur_decomposition(this%x(:,config_type_index), sparseX_index)
call print("Finished CUR decomposition")
case default
RAISE_ERROR('gpCoordinates_sparsify: '//my_sparseMethod//' method is unknown', error)
endselect
this%sparseX_index(li:ui) = config_type_index(sparseX_index)
deallocate(config_type_index,sparseX_index)
if(default_all) exit
enddo
elseif(my_sparseMethod == GP_SPARSE_INDEX_FILE) then
call print('Started reading sparse indices from file '//trim(my_sparse_file))
call fread_array_i(size(this%sparseX_index),this%sparseX_index(1),trim(my_sparse_file)//C_NULL_CHAR)
call print('Finished reading sparse indices from file, '//size(this%sparseX_index)//' of them.')
endif
if(allocated(this%covarianceDiag_sparseX_sparseX)) deallocate(this%covarianceDiag_sparseX_sparseX)
allocate(this%covarianceDiag_sparseX_sparseX(this%n_sparseX))
if(my_sparseMethod == GP_SPARSE_FILE) then
call print('Started reading sparse descriptors from file '//trim(my_sparse_file))
allocate(sparseX_array(d+1,this%n_sparseX))
call fread_array_d(size(sparseX_array),sparseX_array(1,1),trim(my_sparse_file)//C_NULL_CHAR)
this%sparseCutoff = sparseX_array(1,:)
this%sparseX = sparseX_array(2:,:)
this%covarianceDiag_sparseX_sparseX = 1.0_dp ! only used for COVARIANCE_BOND_REAL_SPACE
deallocate(sparseX_array)
call print('Finished reading sparse descriptors from file, '//size(this%sparseCutoff)//' of them.')
else
if(my_sparseMethod == GP_SPARSE_NONE) this%sparseX_index = x_index
call sort_array(this%sparseX_index)
if(this%covariance_type == COVARIANCE_BOND_REAL_SPACE) then
if(allocated(this%sparseX)) deallocate(this%sparseX)
allocate(this%sparseX(maxval(this%x_size(this%sparseX_index)),this%n_sparseX))
if(allocated(this%sparseX_size)) deallocate(this%sparseX_size)
allocate(this%sparseX_size(this%n_sparseX))
this%sparseX(:,:) = this%x(1:maxval(this%x_size(this%sparseX_index)),this%sparseX_index)
this%sparseX_size = this%x_size(this%sparseX_index)
else
this%sparseX(:,:) = this%x(:,this%sparseX_index)
endif
this%covarianceDiag_sparseX_sparseX = this%covarianceDiag_x_x(this%sparseX_index)
this%sparseCutoff = this%cutoff(this%sparseX_index)
if(present(print_sparse_index)) then
if(len_trim(print_sparse_index) > 0) then
call initialise(inout_sparse_index, trim(print_sparse_index), action=OUTPUT, append=.true.)
call print(""//this%sparseX_index,file=inout_sparse_index)
call finalise(inout_sparse_index)
endif
endif
endif
if(allocated(x)) deallocate(x)
if(allocated(x_index)) deallocate(x_index)
if(allocated(x_unique)) deallocate(x_unique)
if(allocated(my_n_sparseX)) deallocate(my_n_sparseX)
this%sparsified = .true.
endsubroutine gpCoordinates_sparsify_config_type
subroutine gpFull_sparsify_array_config_type(this, n_sparseX, default_all, sparseMethod, sparse_file, use_actual_gpcov, print_sparse_index, &
unique_hash_tolerance, unique_descriptor_tolerance, error)
type(gpFull), intent(inout) :: this
integer, dimension(:,:), intent(in) :: n_sparseX
logical, dimension(:), intent(in) :: default_all
integer, dimension(:), optional, intent(in) :: sparseMethod
character(len=STRING_LENGTH), dimension(:), intent(in), optional :: sparse_file, print_sparse_index
logical, intent(in), optional :: use_actual_gpcov
real(dp), dimension(:), optional, intent(in) :: unique_hash_tolerance, unique_descriptor_tolerance
integer, optional, intent(out) :: error
integer :: i
integer, dimension(:), allocatable :: my_sparseMethod
character(len=STRING_LENGTH), dimension(:), allocatable :: my_sparse_file
INIT_ERROR(error)
if( .not. this%initialised ) then
RAISE_ERROR('gpFull_sparsify_array: object not initialised',error)
endif
allocate(my_sparseMethod(this%n_coordinate))
allocate(my_sparse_file(this%n_coordinate))
my_sparseMethod = optional_default((/ (GP_SPARSE_RANDOM, i=1,this%n_coordinate) /),sparseMethod)
my_sparse_file = optional_default((/ ("", i=1,this%n_coordinate) /),sparse_file)
do i = 1, this%n_coordinate
call gpCoordinates_sparsify_config_type(this%coordinate(i),n_sparseX(:,i), default_all(i), &
sparseMethod=my_sparseMethod(i), sparse_file=my_sparse_file(i), use_actual_gpcov=use_actual_gpcov, &
print_sparse_index = print_sparse_index(i), &
unique_hash_tolerance=unique_hash_tolerance(i), unique_descriptor_tolerance=unique_descriptor_tolerance(i), &
error = error)
enddo
if(allocated(my_sparseMethod)) deallocate(my_sparseMethod)
if(allocated(my_sparse_file)) deallocate(my_sparse_file)
endsubroutine gpFull_sparsify_array_config_type
function kernel_distance_matrix(this, config_type_index, covariance_only) result(k_nn)
type(gpCoordinates), intent(in) :: this
integer, dimension(:), intent(in), optional :: config_type_index
logical, intent(in), optional :: covariance_only
real(dp), pointer, dimension(:,:) :: k_nn ! actually the kernel distance matrix
!real(dp), dimension(:,:), allocatable :: k_nn
real(dp), dimension(:), allocatable :: k_self
logical :: do_kernel_distance
integer :: i, j, n, ii, jj
integer :: stat
call system_timer('kernel_distance_matrix')
do_kernel_distance = .not. optional_default(.false., covariance_only)
if(present(config_type_index)) then
n = size(config_type_index)
else
n = size(this%x,2)
endif
allocate(k_self(n))
allocate(k_nn(n,n), stat=stat)
if(stat /= 0) call system_abort('kernel_distance_matrix: could not allocate matrix.')
!$omp parallel do default(none) shared(this,n,config_type_index,k_self) private(i,ii)
do i = 1, n
if(present(config_type_index)) then
ii = config_type_index(i)
else
ii = i
endif
k_self(i) = gpCoordinates_Covariance(this, i_x = ii, j_x = ii, normalise = .false.)
enddo
do j = 1, n
if(present(config_type_index)) then
jj = config_type_index(j)
else
jj = j
endif
!k_nn(j,j) = 1.0_dp ! normalised kernel self-covariance
k_nn(j,j) = 0.0_dp ! distance to itself = 0
!$omp parallel do default(none) shared(n,this,k_nn,jj,j,k_self,config_type_index,do_kernel_distance) private(i,ii)
do i = j+1, n
if(present(config_type_index)) then
ii = config_type_index(i)
else
ii = i
endif
! kernel covariance
k_nn(j,i) = gpCoordinates_Covariance(this, i_x = ii, j_x = jj, normalise = .false.)
! then normalise
k_nn(j,i) = k_nn(j,i) / sqrt(k_self(i)*k_self(j))
if (do_kernel_distance) then
! now convert to distance
k_nn(j,i) = sqrt(2.0_dp * (1.0_dp - k_nn(j,i)))
endif
! finally, symmetrise
k_nn(i,j) = k_nn(j,i)
enddo ! i
enddo ! j
!dm = sqrt(2.0_dp * (1.0_dp - k_nn))
!do i = 1, n
! do j = i+1, n
! dm(i,j) = sqrt(2.0_dp*(1.0_dp - kij))
! dm(j,i) = dm(i,j)
! end do
!end do
!deallocate(k_nn, k_self)
deallocate(k_self)
call system_timer('kernel_distance_matrix')
end function kernel_distance_matrix
subroutine sparse_covariance(this, index_out, config_type_index, use_actual_gpcov)
type(gpCoordinates), intent(in) :: this
integer, dimension(:), intent(out) :: index_out
integer, dimension(:), intent(in), optional :: config_type_index
logical, intent(in), optional :: use_actual_gpcov
real(dp), dimension(:), allocatable :: score, k_self !, xI_xJ
real(dp), dimension(:,:), allocatable :: k_mn, k_mm_k_m
real(dp), dimension(1,1) :: k_mm
integer :: m, n, i, ii, j, jj, i_p
integer, dimension(1) :: j_loc
logical, dimension(:), allocatable :: not_yet_added
logical :: do_use_actual_gpcov
type(LA_Matrix) :: LA_k_mm
call system_timer('sparse_covariance')
if(present(config_type_index)) then
n = size(config_type_index)
else
n = size(this%x,2)
endif
m = size(index_out)
do_use_actual_gpcov = optional_default(.false., use_actual_gpcov)
if(do_use_actual_gpcov) then
call print("sparse_covariance using actual gpCoordinates_Covariance")
else
call print("sparse_covariance using manual 'covariance'")
endif
allocate(k_mn(m,n), score(n), k_mm_k_m(m,n), k_self(n), not_yet_added(n))
k_mn = 0.0_dp
not_yet_added = .true.
!allocate(xI_xJ(this%d))
j = 1
index_out(j) = 1 !ceiling(ran_uniform() * n)
not_yet_added(index_out(j)) = .false.
k_mm = 1.0_dp+1.0e-6_dp
call initialise(LA_k_mm,k_mm)
!$omp parallel do default(none) shared(this,n,config_type_index,k_self,do_use_actual_gpcov) private(i,ii,i_p)
do i = 1, n
if(present(config_type_index)) then
ii = config_type_index(i)
else
ii = i
endif
if(do_use_actual_gpcov) then
k_self(i) = gpCoordinates_Covariance(this, i_x = ii, j_x = ii, normalise = .false.)
else
if(this%covariance_type == COVARIANCE_BOND_REAL_SPACE) then
elseif(this%covariance_type == COVARIANCE_DOT_PRODUCT) then
k_self(i) = dot_product( this%x(:,ii), this%x(:,ii) )**this%zeta
elseif( this%covariance_type == COVARIANCE_ARD_SE ) then
k_self(i) = 0.0_dp
do i_p = 1, this%n_permutations
!xI_xJ = (this%x(this%permutations(:,i_p),i) - this%x(:,j)) / 4.0_dp
k_self(i) = k_self(i) + exp( -0.5_dp * sum((this%x(this%permutations(:,i_p),ii) - this%x(:,ii))**2) / 16.0_dp )
enddo
elseif( this%covariance_type == COVARIANCE_PP ) then
k_self(i) = 0.0_dp
do i_p = 1, this%n_permutations
!xI_xJ = (this%x(this%permutations(:,i_p),i) - this%x(:,j)) / 4.0_dp
k_self(i) = k_self(i) + covariancePP( sqrt( sum((this%x(this%permutations(:,i_p),ii) - this%x(:,ii))**2) ) / 4.0_dp, PP_Q, this%d)
enddo
endif
endif
enddo
do j = 1, m-1
if(present(config_type_index)) then
jj = config_type_index(index_out(j))
else
jj = index_out(j)
endif
!$omp parallel do default(none) shared(n,this,k_mn,jj,j,LA_k_mm,k_mm_k_m,score,k_self,config_type_index, index_out,do_use_actual_gpcov) private(i,i_p,ii)
do i = 1, n
if(present(config_type_index)) then
ii = config_type_index(i)
else
ii = i
endif
if(do_use_actual_gpcov) then
k_mn(j,i) = gpCoordinates_Covariance(this, i_x = ii, j_x = jj, normalise = .false.)
else
if(this%covariance_type == COVARIANCE_BOND_REAL_SPACE) then
elseif(this%covariance_type == COVARIANCE_DOT_PRODUCT) then
k_mn(j,i) = dot_product( this%x(:,ii), this%x(:,jj) )**this%zeta
elseif( this%covariance_type == COVARIANCE_ARD_SE ) then
k_mn(j,i) = 0.0_dp
do i_p = 1, this%n_permutations
!xI_xJ = (this%x(this%permutations(:,i_p),i) - this%x(:,j)) / 4.0_dp
k_mn(j,i) = k_mn(j,i) + exp( -0.5_dp * sum((this%x(this%permutations(:,i_p),ii) - this%x(:,jj))**2) / 16.0_dp )
enddo
elseif( this%covariance_type == COVARIANCE_PP ) then
k_mn(j,i) = 0.0_dp
do i_p = 1, this%n_permutations
!xI_xJ = (this%x(this%permutations(:,i_p),i) - this%x(:,j)) / 4.0_dp
k_mn(j,i) = k_mn(j,i) + covariancePP( sqrt( sum((this%x(this%permutations(:,i_p),ii) - this%x(:,jj))**2) ) / 4.0_dp, PP_Q, this%d)
enddo
endif
endif
k_mn(j,i) = k_mn(j,i) / sqrt(k_self(i)*k_self(index_out(j)))
call Matrix_Solve(LA_k_mm,k_mn(1:j,i),k_mm_k_m(1:j,i))
score(i) = sum( k_mn(1:j,i) * k_mm_k_m(1:j,i) )
enddo
j_loc = minloc(score, mask=not_yet_added)
jj = j_loc(1)
index_out(j+1) = jj
not_yet_added(jj) = .false.
if(j == 1) then
call print('Initial score: '//score)
endif
call print('Min score: '//minval(score))
!k_mm(1:j_i,j_i+1) = k_mn(1:j_i,j)
!k_mm(j_i+1,1:j_i) = k_mn(1:j_i,j)
!k_mm(j_i+1,j_i+1) = 1.0_dp
call LA_Matrix_Expand_Symmetrically(LA_k_mm,(/k_mn(1:j,jj),1.0_dp+1.0e-6_dp/))
!call initialise(LA_k_mm,k_mm(1:j_i+1,1:j_i+1))
enddo
call print('Final score: '//score)
call print('Min score: '//minval(score))
deallocate(k_mn, score, k_mm_k_m, k_self, not_yet_added)
!if(allocated(xI_xJ)) deallocate(xI_xJ)
call finalise(LA_k_mm)
call system_timer('sparse_covariance')
endsubroutine sparse_covariance
end module gp_fit_module
3rdparty/gap/gp_predict.f95 deleted 100644 → 0
View file @ 23bcf97f
Source diff could not be displayed: it is too large. Options to address this: view the blob.
3rdparty/gap/make_permutations_v2.f95 deleted 100644 → 0
View file @ 23bcf97f
!!$
!!$------------Permutation Generator---Alan Nichol---------------------------------
!!$ Generate permutations of the interatomic distance vector of
!!$ a number of atoms. Information about the symmetries present in the cluster
!!$ is specified as an array called 'equivalents' - this is generated automatically
!!$ when a permutation_data_type is initialised with one or more 'signatures' which
!!$ are integer arrays of the atomic numbers.
!!$
!!$ For systems where not all atoms of the same Z are equivalent, the symmetries can
!!$ be specified beforehand by passing an equivalents array to permutation_data_initialise
!!$ An example for toluene would be the following:
!!$
!!$ Toluene : C_6 H_5 - CH_3
!!$
!!$ Atoms in order
!!$ C1 C2 H C3 H C4 H C5 H C6 H C7 H H H
!!$
!!$ where C1 is the tertiary carbon and C2-C6 go in order about the benzene ring
!!$ C7 is the methyl carbon
!!$
!!$
!!$ 15 Atoms in total, and two permutational symmetries. So equivalents array is 2X15
!!$ Symmetry of the 3 methyl hydrogens is specified by the first row:
!!$ ( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 11, 111)
!!$
!!$ And the symmetry of the benzene ring is specified by the second row:
!!$ ( 0, 1, 2, 3, 4, 0, 0, 33, 44, 11, 22, 0, 0 ,0 ,0)
!!$
!!$--------------------------------------------------------------------------------
!!$--------------------------------------------------------------------------------
#include "error.inc"
module permutation_maker_module
use error_module
use system_module, only : dp, print, optional_default, system_timer, operator(//)
implicit none
type permutation_data_type
integer :: perm_number, n_perms
integer, dimension(:), allocatable :: signature_one, signature_two, signature_three, counter, rank, dist_vec !dist_vec is internal name for descriptor
integer, dimension(:,:), allocatable :: dist_vec_permutations
integer, dimension(:,:,:), allocatable :: perm_array
logical :: internal_swaps_only, initialised
endtype permutation_data_type
!% Overloaded assigment operators for permutation data objects.
private :: permutation_data_assignment
interface assignment(=)
module procedure permutation_data_assignment
end interface assignment(=)
contains
subroutine permutation_data_assignment(to,from)
!this is the slow way to copy this object, because it goes through the motions of initialisation again
! call permutation_data_copy for a much faster alternative
type(permutation_data_type), intent(inout) :: to
type(permutation_data_type), intent(in) :: from
! We do not fail if *from* is unitialised, since overloaded operator
! routines are outside scope of error handling mechanism.
if(.not. from%initialised) then
call permutation_data_finalise(to)
return
end if
call permutation_data_initialise(to,signature_one=from%signature_one,signature_two=from%signature_two,signature_three=from%signature_three,internal_swaps_only=from%internal_swaps_only)
end subroutine permutation_data_assignment
subroutine permutation_data_copy(to,from)
type(permutation_data_type), intent(inout) :: to
type(permutation_data_type), intent(in) :: from
if(.not. from%initialised) then
call permutation_data_finalise(to)
return
end if
allocate(to%counter(size(from%counter)))
allocate(to%rank(size(from%rank)))
allocate(to%dist_vec(size(from%dist_vec)))
allocate(to%dist_vec_permutations(size(from%dist_vec_permutations,1),size(from%dist_vec_permutations,2)))
allocate(to%perm_array(size(from%perm_array,1),size(from%perm_array,2),size(from%perm_array,3)))
if (allocated(from%signature_one)) then
allocate(to%signature_one(size(from%signature_one)))
to%signature_one = from%signature_one
if (allocated(from%signature_two)) then
allocate(to%signature_two(size(from%signature_two)))
to%signature_two = from%signature_two
if (allocated(from%signature_three)) then
allocate(to%signature_three(size(from%signature_three)))
to%signature_three = from%signature_three
end if
end if
end if
to%counter = from%counter
to%rank = from%rank
to%dist_vec = from%dist_vec
to%dist_vec_permutations = from%dist_vec_permutations
to%perm_array = from%perm_array
to%n_perms = from%n_perms
to%internal_swaps_only = from%internal_swaps_only
to%initialised = .true.
to%perm_number = 1
end subroutine permutation_data_copy
subroutine equivalents_row_atoms(equivalents_row,signature,offset,N)
implicit none
integer, dimension(:), allocatable, intent(inout) :: equivalents_row
integer, dimension(:), intent(in) :: signature
integer, dimension(:), allocatable :: scratch_row, equivalents_temp
integer, intent(in) :: offset, N
integer :: z_index, i, j, repeats
allocate(scratch_row(N))
allocate(equivalents_temp(1))
do z_index=1,maxval(signature)
repeats=0
scratch_row=0
do i=1,size(signature)
if (signature(i) == z_index) then
do j=repeats,0,-1
scratch_row(i+offset)=scratch_row(i+offset)+10**(j)
end do
repeats = repeats+1
end if
end do
if (repeats .le. 1) cycle
if (.not. allocated(equivalents_row)) then
allocate(equivalents_row(N))
equivalents_row=scratch_row
else
deallocate(equivalents_temp)
allocate(equivalents_temp(size(equivalents_row)))
equivalents_temp=equivalents_row
deallocate(equivalents_row)
allocate(equivalents_row(size(equivalents_temp)+N))
equivalents_row=(/equivalents_temp,scratch_row/)
end if
end do
end subroutine equivalents_row_atoms
subroutine equivalents_row_monomers(equivalents_row,N,signature,pos_a,pos_b,pos_c)
integer,dimension(:), allocatable, intent(inout) :: equivalents_row
integer, intent(in) :: N, pos_a, pos_b
integer, intent(in), optional :: pos_c
integer :: i
integer, dimension(:), intent(in) :: signature
integer, dimension(:), allocatable :: scratch_row, equivalents_temp
allocate(scratch_row(N))
allocate(equivalents_temp(1))
scratch_row=0
do i=1,size(signature)
scratch_row(i+pos_a) = i
scratch_row(i+pos_b) = i*11
if (present(pos_c)) then
scratch_row(i+pos_c) = i*111
end if
end do
if (.not. allocated(equivalents_row)) then
allocate(equivalents_row(N))
equivalents_row=scratch_row
else
deallocate(equivalents_temp)
allocate(equivalents_temp(size(equivalents_row)))
equivalents_temp=equivalents_row
deallocate(equivalents_row)
allocate(equivalents_row(size(equivalents_temp)+N))
equivalents_row=(/equivalents_temp,scratch_row/)
end if
end subroutine equivalents_row_monomers
subroutine permutation_data_initialise(this,equivalents_input,signature_one,signature_two,signature_three,internal_swaps_only,error)
! better make sure that all these optional things are fed in as key=value, because relying on ordering is asking for trouble
implicit none
type(permutation_data_type) :: this
integer, dimension(:), allocatable :: counter, rank, dist_vec, equivalents_row, scratch_row, &
equivalents_temp, atoms, group
integer, dimension(:), allocatable :: signature
integer, dimension(:), optional :: signature_one, signature_two, signature_three
! integer, dimension(:) :: signature_one
integer, dimension(:,:), optional :: equivalents_input
integer, dimension(:,:), allocatable :: group_array, dist_vec_permutations, equivalents
integer, dimension(:,:,:), allocatable :: perm_array
integer, optional, intent(out) :: error
integer :: repeats, num_groups, N, dist_vec_n_perms, i,j,z_index,max_rank,num_distances, &
num_perms, offset_one, offset_two, offset_three, n_atoms_one, n_atoms_two, n_atoms_three
logical, optional :: internal_swaps_only
logical :: two_monomers_given, three_monomers_given, my_internal_swaps_only
real(dp) :: cutoff
INIT_ERROR(error)
call permutation_data_finalise(this)
two_monomers_given=.false.
three_monomers_given=.false.
! automatic generation of equivalents array based on atomic numbers
if (present(equivalents_input)) then
if (present(signature_one)) then
RAISE_ERROR('permutation_data_initialise: mixing of automatically generated permutations and pre-specified symmetries not well defined', error)
end if
N= size(equivalents_input,2)
num_groups = size(equivalents_input,1)
allocate(equivalents(num_groups,N))
equivalents = equivalents_input
else
if (.not. present(signature_one)) then
RAISE_ERROR('permutation_data_initialise doesnt know which permutations to do, provide a signature or equivalents array', error)
end if
n_atoms_one=size(signature_one)
if(present(signature_two)) then
if ( size(signature_two) .ge. 1) then
two_monomers_given= .true.
n_atoms_two=size(signature_two)
if(present(signature_three)) then
if ( size(signature_three) .ge. 1) then
three_monomers_given= .true.
n_atoms_three=size(signature_three)
end if
end if
end if
end if
my_internal_swaps_only = optional_default(.true., internal_swaps_only)
if (three_monomers_given) then
N = size(signature_one)+size(signature_two)+size(signature_three)
else if (two_monomers_given) then
N = size(signature_one)+size(signature_two)
else
N = size(signature_one)
end if
allocate(scratch_row(N))
allocate(equivalents_temp(1))
if (two_monomers_given .and. .not. three_monomers_given) then
if (n_atoms_one .eq. n_atoms_two) then
if (count(signature_one .ne. signature_two) .eq. 0) then
offset_one=0
offset_two=size(signature_one)
call equivalents_row_monomers(equivalents_row,N,signature_one,offset_one,offset_two)
end if
end if
else if (three_monomers_given) then
if (n_atoms_one .eq. n_atoms_two) then
if (count(signature_one .ne. signature_two) .eq. 0) then
! one and two are equivalent
if (n_atoms_one .eq. n_atoms_three) then
if (count(signature_one .ne. signature_three) .eq. 0) then
! all three are equivalent
offset_one=0
offset_two=size(signature_one)
offset_three =size(signature_one)+size(signature_two)
call equivalents_row_monomers(equivalents_row,N,signature_one,offset_one,offset_two,offset_three)
end if
else
! only one and two are equivalent
offset_one=0
offset_two=size(signature_one)
call equivalents_row_monomers(equivalents_row,N,signature_one,offset_one,offset_two)
end if
end if
else if (n_atoms_one .eq. n_atoms_three) then
if (count(signature_one .ne. signature_three) .eq. 0) then
! only one and three are equivalent
offset_one=0
offset_two=size(signature_one)+size(signature_two)
call equivalents_row_monomers(equivalents_row,N,signature_one,offset_one,offset_two)
end if
else if (n_atoms_two .eq. n_atoms_three) then
if (count(signature_two .ne. signature_three) .eq. 0) then
! only two and three are equivalent
offset_one=size(signature_one)
offset_two=size(signature_one)+size(signature_two)
call equivalents_row_monomers(equivalents_row,N,signature_two,offset_one,offset_two)
end if
end if
end if
! if more than one signature is given, and swapping atoms between monomers is allowed
! then the signatures get concatenated to 'signature'. Otherwise 'signature' is just set
! to refer to signature_one
if(two_monomers_given) then
if (my_internal_swaps_only) then
allocate(signature(size(signature_one)))
signature = signature_one
else if (three_monomers_given) then
allocate(signature(size(signature_one)+size(signature_two)+size(signature_three)))
signature = (/signature_one,signature_two,signature_three/)
else
allocate(signature(size(signature_one)+size(signature_two)))
signature = (/signature_one,signature_two/)
end if
else
allocate(signature(size(signature_one)))
signature = signature_one
end if
offset_one=0
call equivalents_row_atoms(equivalents_row,signature,offset_one,N)
if (two_monomers_given .and. my_internal_swaps_only) then
offset_one=size(signature)
call equivalents_row_atoms(equivalents_row,signature_two,offset_one,N)
if (three_monomers_given .and. my_internal_swaps_only) then
offset_one =size(signature)+size(signature_two)
call equivalents_row_atoms(equivalents_row,signature_three,offset_one,N)
end if
end if
if (.not. allocated(equivalents_row)) then
num_groups=0 ! only identity permutation
else
num_groups=size(equivalents_row)/N
allocate(equivalents(num_groups,N))
! make the equivalents array
equivalents =transpose(reshape(equivalents_row,(/ size(equivalents, 2), size(equivalents, 1) /)))
endif
end if
!!$write(*,*) 'equivalents array'
!!$ do i=1,size(equivalents,1)
!!$ write(*,*) equivalents(i,:)
!!$ end do
!--------- Further Array allocations and Initialisation --------------------------!
allocate(atoms(N))
allocate(group(N))
do i=1,size(atoms)
atoms(i)=i
end do
num_distances = N*(N-1)/2
allocate(dist_vec(num_distances))
allocate(counter(num_groups))
allocate(rank(num_groups))
!make rank vector
do i=1,num_groups
group(:) = equivalents(i,:)
num_perms = num_group_perms(group)
rank(i) = num_perms
end do
max_rank = maxval(rank)
!get total number of permutations
dist_vec_n_perms = 1
do i=1,size(rank)
dist_vec_n_perms = dist_vec_n_perms*rank(i)
end do
!initialise counter
counter=1
allocate(dist_vec_permutations(num_distances,dist_vec_n_perms))
allocate(group_array(max_rank,N))
allocate(perm_array(num_groups,N,max_rank))
!initialise arrays permutations to zero
perm_array = 0
dist_vec_permutations=0
!-------------------------------------------------------------------------!
!make 2D array of permutations of each group and add to 3D array perm_array
!-------------------------------------------------------------------------!
do i = 1,num_groups
group(:) = equivalents(i,:)
group_array = permute_atoms(atoms,group,N,max_rank)!this padded with zeroes in case group is of less than max_rank
do j=1,size(group_array, 1)
perm_array(i,:,j) = group_array(j,:)
end do
end do
!-------------------------------------------------------------------------!
!Now assign relevant stuff to the permutation_data_type
!-------------------------------------------------------------------------!
if (present(signature_one)) then
allocate(this%signature_one(size(signature_one)))
this%signature_one=signature_one
end if
if (two_monomers_given) then
allocate(this%signature_two(size(signature_two)))
this%signature_two=signature_two
end if
if (three_monomers_given) then
allocate(this%signature_three(size(signature_three)))
this%signature_three=signature_three
end if
!- If there's only one permutation then write it out explicitly here
if (num_groups == 0) then
do i=1,num_distances
dist_vec_permutations(i,1)=i
end do
end if
allocate(this%counter(size(counter)))
allocate(this%rank(size(rank)))
allocate(this%perm_array(size(perm_array,1),size(perm_array,2),size(perm_array,3)))
allocate(this%dist_vec(num_distances))
allocate(this%dist_vec_permutations(size(dist_vec_permutations,1),size(dist_vec_permutations,2)))
this%counter=counter
this%rank=rank
this%perm_array=perm_array
this%dist_vec_permutations=dist_vec_permutations
this%perm_number=1
this%n_perms=dist_vec_n_perms
this%initialised=.true.
end subroutine permutation_data_initialise
subroutine permutation_data_finalise(this)
implicit none
type(permutation_data_type) :: this
if (.not. this%initialised) return
if(allocated(this%signature_one)) deallocate(this%signature_one)
if(allocated(this%signature_two)) deallocate(this%signature_two)
if(allocated(this%counter)) deallocate(this%counter)
if(allocated(this%rank)) deallocate(this%rank)
if(allocated(this%perm_array)) deallocate(this%perm_array)
if(allocated(this%dist_vec)) deallocate(this%dist_vec)
if(allocated(this%dist_vec_permutations)) deallocate(this%dist_vec_permutations)
this%initialised = .false.
end subroutine permutation_data_finalise
subroutine add_combined_permutation (counter, perm_array, dist_vec_permutations,perm_number)
implicit none
! this gets called by the subroutine next, it should receive a vector 'counter' from which it
! figures out which permutations to combine. It then asks combine_perms to do so and
! gets the dist_vec vector from do_swaps
integer :: i, num_distances, N, perm_number
integer, dimension(:), intent(inout) :: counter
integer, dimension(:), allocatable :: combo, next_perm, dist_vec
integer, dimension(:,:,:), intent(in) :: perm_array
integer, dimension(:,:) :: dist_vec_permutations
N=size(perm_array,2)
num_distances = N*(N-1)/2
allocate(dist_vec(num_distances))
allocate(combo(N))
allocate(next_perm(N))
combo = perm_array(1,:,counter(1))
do i=1, size(counter)-1
next_perm = perm_array(i+1,:,counter(i+1))
combo = combine_perms(combo,next_perm)
end do
call do_swaps(combo, dist_vec)
dist_vec_permutations(:,perm_number)=dist_vec
deallocate(dist_vec)
deallocate(combo)
deallocate(next_perm)
end subroutine add_combined_permutation
recursive subroutine next(this, m)
implicit none
type(permutation_data_type), intent(inout) :: this
integer :: m, num_groups
num_groups = size(this%counter)
if (m .gt. num_groups) then
call add_combined_permutation(this%counter, this%perm_array, this%dist_vec_permutations, this%perm_number)
this%perm_number=this%perm_number+1
else
do while (this%counter(m) .lt. this%rank(m))
call next(this, m+1)
this%counter(m+1:) = 1
this%counter(m) = this%counter(m) + 1
end do
this%counter(m+1:) = 1
call next(this,m+1)
end if
end subroutine next
function num_group_perms(group)
implicit none
integer :: num_group_perms, n_members
integer, dimension(:) :: group
integer, dimension(8) :: factorial
factorial = (/ 1,2,6,24,120,720,5040,40320 /)
n_members = ceiling(log10(real(maxval(group))))
num_group_perms = factorial(n_members)
return
end function num_group_perms
! This modified from Rosetta Code
recursive subroutine update_matrix(std_perms,n_members,position,i_perm,perm_vec)
implicit none
integer :: n_members, value, position, i_perm
integer, dimension(:,:) :: std_perms
integer, dimension(:) :: perm_vec
if (position > n_members) then
std_perms(i_perm,:) = perm_vec
i_perm=i_perm+1
else
do value = 1, n_members
if (.not. any (perm_vec(:position - 1) == value)) then
perm_vec(position)= value
call update_matrix(std_perms,n_members,position+1,i_perm,perm_vec)
end if
end do
end if
end subroutine update_matrix
function permute_atoms(atoms,group,N,max_rank)
implicit none
integer :: i, j, k, i_perm, n_members, num_perms, atoms_per_monomer
integer, intent(IN) :: N, max_rank
integer, dimension(N) :: atoms, group
integer, dimension(:,:), allocatable :: permute_atoms, std_perms
integer, dimension(:), allocatable :: group_vec, perm_vec, indices, offsets
integer, dimension(1) :: p, q, temp
n_members = ceiling(log10(real(maxval(group))))
num_perms = num_group_perms(group)
allocate(group_vec(n_members))
allocate(perm_vec(n_members))
!allocate(indices(n_members))
allocate(std_perms(num_perms,n_members))
allocate(permute_atoms(max_rank,N))
permute_atoms = 0
if (num_perms .eq. 2) then
!just a pair of equivalent atoms or monomers
permute_atoms(1,:) = atoms
permute_atoms(2,:) = atoms
do i=1,count(group .lt. 10 .and. group .gt. 0)
p = minloc(group, mask=group .ge. i)
q = minloc(group, mask=group .gt. 10*i)
!write(*,*) "equivalent pair"
!write(*,'(2I3)') p,q
temp = permute_atoms(2,p(1))
permute_atoms(2,p(1)) = permute_atoms(2,q(1))
permute_atoms(2,q(1)) = temp(1)
end do
else
!Permutations of groups of >2 atoms, no support for >2 monomers yet
i_perm=1
call update_matrix(std_perms,n_members,1,i_perm,perm_vec)
if (.not. any(group .eq. 2)) then ! permutations of identical atoms
allocate(indices(n_members))
! get indices of equivalent atoms
do i=1,n_members
temp =minloc(group, mask = group .ge. 10**(i-1))
indices(i) = temp(1)
end do
do i=1,size(std_perms,1)
perm_vec = std_perms(i,:)
group_vec = indices(perm_vec)
do j=1,n_members
permute_atoms(i,indices(j)) = group_vec(j)
end do
do j=1,N
if (permute_atoms(i,j) ==0) permute_atoms(i,j) = j
end do
end do
else ! permutations of identical monomers
allocate(offsets(n_members))
atoms_per_monomer = maxval(group, mask = group .lt. 10)
! allocate(indices(n_members*atoms_per_monomer))
do i=1,n_members
! find the 1, 11, 111, etc. with which the monomer starts
temp =minloc(group, mask = group .ge. 10**(i-1))
offsets(i) = temp(1)-1
!!$ do j=1,atoms_per_monomer
!!$ indices(j+(i-1)*n_members) = j + offsets(i)
!!$ end do
end do
do i=1,size(std_perms,1)
perm_vec = std_perms(i,:)
group_vec = offsets(perm_vec)
do j=1,n_members
do k=1,atoms_per_monomer
permute_atoms(i,k+offsets(j)) = k+group_vec(j)
end do
end do
do j=1,N
if (permute_atoms(i,j) ==0) permute_atoms(i,j) = j
end do
end do
end if
end if
return
end function permute_atoms
function combine_perms(vec1,vec2)
implicit none
integer, dimension(:), intent(in) :: vec1, vec2
integer, dimension(:), allocatable :: combine_perms
integer :: j
allocate(combine_perms(size(vec1)))
if (size(vec1) /= size(vec2)) then
write(*,*) "combine_perms received vectors of mismatched lengths"
call exit(1)
end if
do j=1,size(vec1)
combine_perms(j) = vec1(vec2(j))
end do
return
end function combine_perms
subroutine do_swaps(atom_vec, dist_vec)
implicit none
integer :: N, start, finish, length, temp, j, i
integer, dimension(:), intent(in) :: atom_vec
integer, dimension(:) :: dist_vec
integer, dimension(1) :: i_vec
integer, dimension(:), allocatable :: temp_vec, scratch_vec!, do_swaps
integer, dimension(:,:), allocatable :: dist_mat, dist_mat_upper
!initialise vector and matrix
N = size(atom_vec)
allocate(scratch_vec(N))
allocate(temp_vec(N))
do i=1,N
temp_vec(i)=i
end do
do i=1,size(dist_vec)
dist_vec(i)=i
end do
allocate(dist_mat(N,N))
allocate(dist_mat_upper(N,N))
dist_mat=0
dist_mat_upper = 0
start = 1
do i=1,N
finish=start + N-i
dist_mat_upper(i,i+1:N) = dist_vec(start:finish-1)
start = finish
end do
dist_mat = dist_mat_upper + transpose(dist_mat_upper)
do while (any(temp_vec .ne. atom_vec))
i_vec = minloc(temp_vec, temp_vec .ne. atom_vec)
i=i_vec(1)
! keep track of swaps
temp = temp_vec(i)
temp_vec(i) = temp_vec(atom_vec(i))
temp_vec(atom_vec(i)) = temp
! now swap in array - rows then columns
scratch_vec = dist_mat(i,:)
dist_mat(i,:) = dist_mat(atom_vec(i),:)
dist_mat(atom_vec(i),:) = scratch_vec
scratch_vec = dist_mat(:,i)
dist_mat(:,i) = dist_mat(:,atom_vec(i))
dist_mat(:,atom_vec(i)) = scratch_vec
end do
!convert back into vector
start = 1
finish=N-1
do i=1,N-1
dist_vec(start:finish) = dist_mat(i,i+1:N)
start = finish+1
finish = finish+N-i-1
end do
deallocate(temp_vec)
deallocate(scratch_vec)
deallocate(dist_mat)
deallocate(dist_mat_upper)
return
end subroutine do_swaps
end module permutation_maker_module
3rdparty/gap/teach_sparse deleted 100755 → 0
View file @ 23bcf97f
#!/bin/bash
echo
echo The teach-sparse program is now called gap_fit
echo
nomad-FAIR @ d9b72a06
Subproject commit d9b72a06689bc850f46b30e96ad3224f2d761e3a
quip @ b7d59496
Subproject commit b7d59496c970679531a4a22c9c6dc6968aeb69ed
3rdparty/quip/.gitignore deleted 100644 → 0
View file @ 23bcf97f
build
doc/_build
*~
*.pyc
.ipynb_checkpoints
.idea/
GitStatistics
src/GAP
src/ThirdParty
src/libAtoms/*.mod
src/Potentials/*.mod
src/Potentials/*.mod.save
quippy/quippy.f90doc
quippy/quippy.spec
3rdparty/quip/.travis.yml deleted 100644 → 0
View file @ 23bcf97f
---
language: python
python:
- "3.6"
env:
matrix:
- QUIP_ARCH=linux_x86_64_gfortran QUIP_INC=Makefile.nogap.inc
- QUIP_ARCH=linux_x86_64_gfortran QUIP_INC=Makefile.inc
- QUIP_ARCH=linux_x86_64_gfortran_openmp QUIP_INC=Makefile.nogap.inc
- QUIP_ARCH=linux_x86_64_gfortran_openmp QUIP_INC=Makefile.inc
stages:
- test
- name: documentation
# freezing the python2.7 version fot this at first
if: branch = archived_py2
- release
- docker
git:
depth: 1
addons:
apt:
update: true
before_install:
- sudo apt-get install -y
gfortran
libblas-dev
liblapack-dev
openmpi-bin
libopenmpi-dev
netcdf-bin
libnetcdf-dev
libhdf5-serial-dev
python-numpy
# Clone the private gap repository using deploy key
- echo $github_deploy_key | base64 -d > ~/.ssh/github_deploy_key
- chmod 0600 ~/.ssh/github_deploy_key
- eval $(ssh-agent -s)
- ssh-add ~/.ssh/github_deploy_key
- git clone --depth 1 git@github.com:libAtoms/GAP.git src/GAP
# Make the build directory manually for rules
- mkdir -p build/${QUIP_ARCH}
# Copy the rules file (rather than make config)
- cp travis/${QUIP_INC} build/${QUIP_ARCH}/Makefile.inc
install:
- pip install -r quippy/requirements.txt
script: # Compile QUIP, libquip and quippy, if a non-skipped build.
# stdout is redirected as it generates too much output.
# Quippy should have built successfully -> start tests
- make
- make libquip
- make quippy
# Sometimes file limit is 64000 and read_loop reads and writes this many
# files causing the build to time out
- ulimit -n 256
- make test
jobs:
include:
- stage: documentation
# freezing the python2.7 version fot this at first
python:
- "2.7"
env: # serial version of QUIP
- QUIP_ARCH=linux_x86_64_gfortran QUIP_INC=Makefile.inc
install: # packages for building docs
- sudo apt-get install -y
libgsl0-dev
libxpm-dev
pandoc
- pip install
f90wrap==0.2.1
ase==3.17.0
docutils==0.14
sphinx
sphinx-rtd-theme
nbsphinx
numpydoc
# needed to nbconvert ipynb files and to process the rst files
- pip install 'nbconvert[execute]' 'ipython<6'
# quippy is working, install it
- make &> /dev/null
- make install-quippy &> /dev/null
# Install atomeye from src
- export QUIP_ROOT=`pwd`
- git clone https://github.com/jameskermode/AtomEye.git src/AtomEye
- cd src/AtomEye/Python
- python setup.py install
- cd ../../..
script:
- cd doc
- make html
deploy:
provider: pages
local_dir: doc/_build/html
skip-cleanup: true
github-token: $GITHUB_TOKEN # Set in the settings page of your repository, as a secure variable
keep-history: true
on:
branch: archived_py2
- stage: release
name: "Trigger for a new release of GAP"
before_install: skip
install: skip
script: skip
deploy:
provider: script
script: bash travis/trigger_gap_release.sh
on:
branch: public
- stage: docker
name: "Trigger for building docker images"
before_install: skip
install: skip
script: skip
deploy:
provider: script
script: bash travis/trigger_quip_docker.sh
on:
branch: public
notifications:
email:
- quip-developers@eng.cam.ac.uk
3rdparty/quip/Makefile deleted 100644 → 0
View file @ 23bcf97f
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
# H0 X
# H0 X libAtoms+QUIP: atomistic simulation library
# H0 X
# H0 X Portions of this code were written by
# H0 X Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
# H0 X Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
# H0 X
# H0 X Copyright 2006-2015.
# H0 X
# H0 X These portions of the source code are released under the GNU General
# H0 X Public License, version 2, http://www.gnu.org/copyleft/gpl.html
# H0 X
# H0 X If you would like to license the source code under different terms,
# H0 X please contact Gabor Csanyi, gabor@csanyi.net
# H0 X
# H0 X Portions of this code were written by Noam Bernstein as part of
# H0 X his employment for the U.S. Government, and are not subject
# H0 X to copyright in the USA.
# H0 X
# H0 X
# H0 X When using this software, please cite the following reference:
# H0 X
# H0 X http://www.libatoms.org
# H0 X
# H0 X Additional contributions by
# H0 X Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
# H0 X
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
.PHONY: config doc clean deepclean distclean install test quippy doc install-structures install-dtds install-Tools install-build.QUIP_ARCH install-quippy libquip
ifndef QUIP_ARCH
$(error "You need to define the architecture using the QUIP_ARCH variable. Check out the arch/ subdirectory.")
else
include arch/Makefile.${QUIP_ARCH}
endif
# include other makefiles and export env variables
export QUIP_ARCH
ifeq (${QUIP_ROOT},)
QUIP_ROOT=${PWD}
endif
export QUIP_ROOT
export SCRIPT_PATH=${QUIP_ROOT}/bin
export BUILDDIR=${QUIP_ROOT}/build/${QUIP_ARCH}${QUIP_ARCH_SUFFIX}
ifneq ($(findstring mpi, ${QUIP_ARCH}),)
QUIP_MPI_SUFFIX=_mpi
endif
export QUIP_MPI_SUFFIX
-include ${BUILDDIR}/Makefile.inc
# create modules list
MODULES = libAtoms
# add any third party packages first, but after libAtoms in case they want to
# use it
ifeq (${HAVE_THIRDPARTY},1)
THIRDPARTY = ThirdParty
MODULES += ThirdParty
THIRDPARTY_LIBS := libthirdparty.a
ifeq (${HAVE_FX},1)
THIRDPARTY_LIBS += libfx.a
endif
ifeq (${HAVE_SCME},1)
THIRDPARTY_LIBS += libscme.a
endif
ifeq (${HAVE_MTP},1)
THIRDPARTY_LIBS += libmtp.a
endif
ifeq (${HAVE_MBD},1)
THIRDPARTY_LIBS += libmbd.a
endif
ifeq (${HAVE_TTM_NF},1)
THIRDPARTY_LIBS += libttm_nf.a
endif
ifeq (${HAVE_CH4},1)
THIRDPARTY_LIBS += libch4.a
endif
endif
# add GAP modules if we have them - they need to come before other modules, except for libAtoms
ifeq (${HAVE_GAP},1)
MODULES += GAP
GAP += GAP
GAP_PROGRAMS = gap_programs
else
GAP =
GAP_PROGRAMS =
endif
# now add the rest of the modules
MODULES += Potentials Utils Programs FilePot_drivers Structure_processors
# diagnostic
$(info Using QUIP_ARCH=${QUIP_ARCH}, MODULES=${MODULES}, QUIP_ROOT=${QUIP_ROOT})
# the first target
all: ${MODULES} ${GAP_PROGRAMS}
FOX = fox
export FOX_LIBDIR=${QUIP_ROOT}/src/${FOX}/objs.${QUIP_ARCH}/lib
export FOX_INCDIR=${QUIP_ROOT}/src/${FOX}/objs.${QUIP_ARCH}/finclude
# now we can include the config makefile, it needs to come after the default target
include Makefile.config
include Makefile.rules
${BUILDDIR}/Makefile.inc:
@if [ "$(MAKECMDGOALS)" != config ]; then\
echo ;\
echo "${BUILDDIR}/Makefile.inc not found. Perhaps you forgot to run \`make config'?" ;\
echo ;\
exit 1 ;\
fi
# Automatically pull the submodules if the user didn't.
# Remove the Makefile.QUIP if it has failed previously.
src/${FOX}/configure:
@echo "Attempting to automatically clone fox submodule"
rm -f src/${FOX}/Makefile.QUIP
git submodule update --init src/${FOX} || \
{ echo "fox clone failed. Download it manually" ; exit 1 ; }
${FOX}: src/${FOX}/configure src/${FOX}/objs.${QUIP_ARCH}/lib/libFoX_common.a
src/${FOX}/objs.${QUIP_ARCH}/lib/libFoX_common.a:
cp Makefile.fox src/${FOX}/Makefile.QUIP
make -C src/${FOX} -I${PWD} -I${PWD}/arch -I${BUILDDIR} -f Makefile.QUIP
FOX_STATIC_LIBFILES = $(patsubst -l%,${FOX_LIBDIR}/lib%.a,${FOX_LIBS})
FOX_STATIC_LIBFILE_OBJS = $(shell for i in ${FOX_STATIC_LIBFILES}; do ar -t $$i; done | grep \.o)
# general rule to make a module
${MODULES}: ${BUILDDIR}/Makefile.inc ${BUILDDIR} ${FOX}
@echo "********************************************"
@echo ""
@echo " Making $@ "
@echo ""
@echo "********************************************"
rm -f ${BUILDDIR}/Makefile
cp ${PWD}/src/$@/Makefile ${BUILDDIR}/Makefile
${MAKE} -C ${BUILDDIR} QUIP_ROOT=${QUIP_ROOT} VPATH=${PWD}/src/$@ -I${PWD} -I${PWD}/arch
rm ${BUILDDIR}/Makefile
# general rule to make a program in the Programs
# src directory, makes sure everything else is
# built first
Programs/% src/Programs/% : ${MODULES}
@echo "********************************************"
@echo ""
@echo " Making Programs "
@echo ""
@echo "********************************************"
rm -f ${BUILDDIR}/Makefile
cp ${PWD}/src/Programs/Makefile ${BUILDDIR}/Makefile
${MAKE} -C ${BUILDDIR} QUIP_ROOT=${QUIP_ROOT} VPATH=${PWD}/src/Programs -I${PWD} -I${PWD}/arch $(lastword $(subst /, ,$@))
rm ${BUILDDIR}/Makefile
# dependencies between modules
ifeq (${HAVE_GAP},1)
GAP: libAtoms ${FOX}
endif
gap_programs: ${MODULES}
@echo "********************************************"
@echo ""
@echo " Making GAP programs "
@echo ""
@echo "********************************************"
rm -f ${BUILDDIR}/Makefile
cp ${PWD}/src/GAP/Makefile ${BUILDDIR}/Makefile
${MAKE} -C ${BUILDDIR} QUIP_ROOT=${QUIP_ROOT} VPATH=${PWD}/src/GAP -I${PWD} -I${PWD}/arch Programs
rm ${BUILDDIR}/Makefile
ThirdParty: libAtoms
Potentials: libAtoms ${GAP}
Utils: libAtoms ${GAP} Potentials
FilePot_drivers: libAtoms Potentials Utils
Programs: libAtoms ${GAP} Potentials Utils
Tests: libAtoms ${GAP} Potentials Utils
libatoms: libAtoms
libquip: libquip.a
libquip.a: ${MODULES}
LIBQUIP_OBJS="$(shell for i in ${BUILDDIR}/libquiputils.a ${BUILDDIR}/libquip_core.a $(addprefix ${BUILDDIR}/,${GAP_LIBFILE}) ${BUILDDIR}/libatoms.a $(addprefix ${BUILDDIR}/,${THIRDPARTY_LIBS}) ${FOX_STATIC_LIBFILES}; do ar -t $$i; done | grep \.o)" && \
cd ${BUILDDIR} && for i in ${FOX_STATIC_LIBFILES}; do ar -x $$i; done && ar -rcs $@ $$LIBQUIP_OBJS
${BUILDDIR}:
@if [ ! -d build/${QUIP_ARCH}${QUIP_ARCH_SUFFIX} ] ; then mkdir -p build/${QUIP_ARCH}${QUIP_ARCH_SUFFIX} ; fi
quippy: libquip.a ${GAP_PROGRAMS}
@echo "********************************************"
@echo ""
@echo " Making quippy "
@echo ""
@echo "********************************************"
rm -f ${BUILDDIR}/Makefile
cp ${PWD}/quippy/Makefile ${BUILDDIR}/Makefile
${MAKE} -C ${BUILDDIR} QUIP_ROOT=${QUIP_ROOT} -I${PWD} -I${PWD}/arch build
rm ${BUILDDIR}/Makefile
install-quippy: quippy
@echo "********************************************"
@echo ""
@echo " Installing quippy "
@echo ""
@echo "********************************************"
rm -f ${BUILDDIR}/Makefile
cp ${PWD}/quippy/Makefile ${BUILDDIR}/Makefile
${MAKE} -C ${BUILDDIR} QUIP_ROOT=${QUIP_ROOT} -I${PWD} -I${PWD}/arch install
rm ${BUILDDIR}/Makefile
clean-quippy:
rm -f ${BUILDDIR}/Makefile
cp ${PWD}/quippy/Makefile ${BUILDDIR}/Makefile
${MAKE} -C ${BUILDDIR} QUIP_ROOT=${QUIP_ROOT} -I${PWD} -I${PWD}/arch clean
rm ${BUILDDIR}/Makefile
clean: ${BUILDDIR}
-${MAKE} clean-quippy
for mods in ${MODULES} ; do \
echo "clean in $$mods"; \
rm -f ${BUILDDIR}/Makefile ; \
cp ${PWD}/src/$$mods/Makefile ${BUILDDIR}/Makefile ; \
${MAKE} -C ${BUILDDIR} USE_MAKEDEP=0 QUIP_ROOT=${QUIP_ROOT} VPATH=${PWD}/src/$$mods -I${PWD} -I${PWD}/arch clean ; \
done
rm -f ${BUILDDIR}/libquip.a
rm -rf src/${FOX}/objs.${QUIP_ARCH}
deepclean: clean
distclean: clean
rm -rf build
install-structures:
ifeq (${QUIP_STRUCTS_DIR},)
@echo
@echo "QUIP_STRUCTS_DIR must be defined to install structures"
else
${MAKE} -C share/Structures QUIP_STRUCTS_DIR=$(QUIP_STRUCTS_DIR) install
endif
install: ${MODULES} install-structures
ifeq (${QUIP_INSTALLDIR},)
@echo
@echo "'make install' needs QUIP_INSTALLDIR to be defined to install "
@echo "programs"
else
if [ ! -d ${QUIP_INSTALLDIR} ]; then \
echo "make install: QUIP_INSTALLDIR '${QUIP_INSTALLDIR}' doesn't exist or isn't a directory"; \
exit 1; \
else \
for mods in ${MODULES} ; do \
rm -f ${BUILDDIR}/Makefile ;\
cp ${PWD}/src/$$mods/Makefile ${BUILDDIR}/Makefile ;\
${MAKE} -C ${BUILDDIR} QUIP_ROOT=${QUIP_ROOT} -I${PWD} -I${PWD}/arch install ;\
rm ${BUILDDIR}/Makefile ;\
done ;\
fi
endif
test: quippy
#- cd tests
#- make
#- cd ..
${MAKE} -C tests -I${PWD} -I${PWD}/arch -I${BUILDDIR}
GIT_SUBDIRS=src/GAP src/ThirdParty
git_pull_all:
git pull --recurse-submodules
@for d in ${GIT_SUBDIRS}; do if [ -d $$d ]; then pushd $$d; git pull; popd; fi; done
distribution:
./bin/gitversion > GIT_VERSION
git archive HEAD > ../QUIP.distribution.`date +%Y-%m-%d`.tar
tar rvf ../QUIP.distribution.`date +%Y-%m-%d`.tar GIT_VERSION
bzip2 ../QUIP.distribution.`date +%Y-%m-%d`.tar
rm GIT_VERSION
@for d in ${GIT_SUBDIRS}; do if [ -d $$d ]; then pushd $$d; git archive HEAD | bzip2 > ../../$$d.distribution.`date +%Y-%m-%d`.tar.bz2; popd; fi; done
3rdparty/quip/Makefile.config deleted 100644 → 0
View file @ 23bcf97f
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
# H0 X
# H0 X libAtoms+QUIP: atomistic simulation library
# H0 X
# H0 X Portions of this code were written by
# H0 X Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
# H0 X Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
# H0 X
# H0 X Copyright 2006-2010.
# H0 X
# H0 X These portions of the source code are released under the GNU General
# H0 X Public License, version 2, http://www.gnu.org/copyleft/gpl.html
# H0 X
# H0 X If you would like to license the source code under different terms,
# H0 X please contact Gabor Csanyi, gabor@csanyi.net
# H0 X
# H0 X Portions of this code were written by Noam Bernstein as part of
# H0 X his employment for the U.S. Government, and are not subject
# H0 X to copyright in the USA.
# H0 X
# H0 X
# H0 X When using this software, please cite the following reference:
# H0 X
# H0 X http://www.libatoms.org
# H0 X
# H0 X Additional contributions by
# H0 X Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
# H0 X
# H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
SHELL := /bin/bash
config: ${BUILDDIR} arch
@echo
@echo "In the following, hit enter to accept the defaults."
@echo
ifneq (${shell [ -f ${BUILDDIR}/Makefile.inc ] && echo 1 || echo 0},1)
@echo "# Place to override setting elsewhere, in particular things set in Makefile.${QUIP_ARCH}" > ${BUILDDIR}/Makefile.inc
@echo "# look in ${QUIP_ROOT}/arch/Makefile.${QUIP_ARCH} for defaults set by arch" >> ${BUILDDIR}/Makefile.inc
@echo "# " >> ${BUILDDIR}/Makefile.inc
@echo "# F77=${F77}" >> ${BUILDDIR}/Makefile.inc
@echo "# F90=${F90}" >> ${BUILDDIR}/Makefile.inc
@echo "# F95=${F95}" >> ${BUILDDIR}/Makefile.inc
@echo "# CC=${CC}" >> ${BUILDDIR}/Makefile.inc
@echo "# CPLUSPLUS=${CPLUSPLUS}" >> ${BUILDDIR}/Makefile.inc
@echo "# FPP=${FPP}" >> ${BUILDDIR}/Makefile.inc
@echo "# LINKER=${LINKER}" >> ${BUILDDIR}/Makefile.inc
@echo "# LIBTOOL=${LIBTOOL}" >> ${BUILDDIR}/Makefile.inc
@echo "# OPTIM=" >> ${BUILDDIR}/Makefile.inc
@echo "# COPTIM=" >> ${BUILDDIR}/Makefile.inc
ifdef SAMPLE_DEBUG
@echo "# DEBUG=${SAMPLE_DEBUG}" >> ${BUILDDIR}/Makefile.inc
endif
@echo "# DEBUG=${DEBUG}" >> ${BUILDDIR}/Makefile.inc
@echo "# CDEBUG=" >> ${BUILDDIR}/Makefile.inc
endif
# BUILD OPTIONS
# #############
@echo
@echo "BUILD OPTIONS"
@echo "============="
@echo
ifndef MATH_LINKOPTS
@echo
@echo "Please enter the linking options for LAPACK and BLAS libraries:"
@echo " Default: ${DEFAULT_MATH_LINKOPTS}"
@read MATH_LINKOPTS && if [[ $$MATH_LINKOPTS == "" ]] ; then \
echo "MATH_LINKOPTS=${DEFAULT_MATH_LINKOPTS}" >> ${BUILDDIR}/Makefile.inc ; \
else echo "MATH_LINKOPTS=$$MATH_LINKOPTS" >> ${BUILDDIR}/Makefile.inc ; fi
endif
ifeq ($(origin EXTRA_LINKOPTS), undefined)
@echo
@echo "Please enter any other extra linking options:"
@echo " Default: none"
@read EXTRA_LINKOPTS; echo "EXTRA_LINKOPTS=$$EXTRA_LINKOPTS" >> ${BUILDDIR}/Makefile.inc
endif
ifndef USE_MAKEDEP
@echo
@echo "Would you like to use makedep? [y/n]"
@echo " Default: y"
@read USE_MAKEDEP_Q && if [[ $$USE_MAKEDEP_Q == 'n' ]]; then \
echo "USE_MAKEDEP=0" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "USE_MAKEDEP=1" >> ${BUILDDIR}/Makefile.inc; \
echo "MAKEDEP=makedep.py" >> ${BUILDDIR}/Makefile.inc ; fi
endif
# OPTIONAL INTERNAL FEATURES
# ##########################
@echo
@echo "OPTIONAL INTERNAL FEATURES"
@echo "=========================="
@echo
ifndef HAVE_CP2K
@echo
@echo "Would you like to compile with CP2K support? [y/n]"
@echo " Default: n"
@read HAVE_CP2K_Q && if [[ $$HAVE_CP2K_Q == 'y' ]]; then \
echo "HAVE_CP2K=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_CP2K=0" >> ${BUILDDIR}/Makefile.inc ; fi
endif
ifndef HAVE_VASP
@echo
@echo "Would you like to compile with VASP support? [y/n]"
@echo " Default: n"
@read HAVE_VASP_Q && if [[ $$HAVE_VASP_Q == 'y' ]]; then \
echo "HAVE_VASP=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_VASP=0" >> ${BUILDDIR}/Makefile.inc ; fi
endif
ifndef HAVE_TB
@echo
@echo "Would you like to compile with Tight Binding (TB) support? [y/n]"
@echo " Default: n"
@read HAVE_TB_Q && if [[ $$HAVE_TB_Q == 'y' ]]; then \
echo "HAVE_TB=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_TB=0" >> ${BUILDDIR}/Makefile.inc ; fi
endif
ifndef HAVE_PRECON
@echo
@echo "Would you like to compile with preconditioned minimisation support? (not working with Intel compiler) [y/n]"
@echo " Default: y"
@read HAVE_PRECON_Q && if [[ $$HAVE_PRECON_Q == 'n' ]]; then \
echo "HAVE_PRECON=0" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_PRECON=1" >> ${BUILDDIR}/Makefile.inc ; fi
endif
ifndef HAVE_LOTF
@echo
@echo "Would you like to compile with Learn-on-the-fly (LOTF) support? [y/n]"
@echo " Default: n"
@read HAVE_LOTF_Q && if [[ $$HAVE_LOTF_Q == 'y' ]]; then \
echo "HAVE_LOTF=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_LOTF=0" >> ${BUILDDIR}/Makefile.inc ; fi
endif
ifndef HAVE_ONIOM
@echo
@echo "Would you like to compile with ONIOM support? [y/n]"
@echo " Default: n"
@read HAVE_ONIOM_Q && if [[ $$HAVE_ONIOM_Q == 'y' ]]; then \
echo "HAVE_ONIOM=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_ONIOM=0" >> ${BUILDDIR}/Makefile.inc ; fi
endif
ifndef HAVE_LOCAL_E_MIX
@echo
@echo "Would you like to compile with LOCAL_E_MIX support? [y/n]"
@echo " Default: n"
@read HAVE_LOCAL_E_MIX_Q && if [[ $$HAVE_LOCAL_E_MIX_Q == 'y' ]]; then \
echo "HAVE_LOCAL_E_MIX=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_LOCAL_E_MIX=0" >> ${BUILDDIR}/Makefile.inc ; fi
endif
ifndef HAVE_QC
@echo
@echo "Would you like to compile with QuasiContinuum wrapper support? [y/n]"
@echo " Default: n"
@read HAVE_QC_Q && if [[ $$HAVE_QC_Q == 'y' ]]; then \
echo "HAVE_QC=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_QC=0" >> ${BUILDDIR}/Makefile.inc ; fi
endif
# OPTIONAL FEATURES (EXTRA CODE REQUIRED)
# #######################################
@echo
@echo "OPTIONAL FEATURES (EXTRA CODE REQUIRED)"
@echo "======================================="
@echo
ifndef HAVE_GAP
@echo
@if [[ ! -d src/GAP ]] ; then \
echo "GAP sources must be acquired separately and placed in the src/GAP subdirectory."; \
fi
@echo "Would you like to compile with GAP support ? [y/n]"
@echo " Default: n"
@read HAVE_GAP_Q && if [[ $$HAVE_GAP_Q == 'y' ]]; then \
echo "HAVE_GAP=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_GAP=0" >> ${BUILDDIR}/Makefile.inc ; fi
endif
ifndef HAVE_QR
@echo
@echo "Would you like to turn on QR decomposition for some GAP operations? [y/n]"
@echo " Default: y"
@read HAVE_QR_Q && if [[ $$HAVE_QR_Q == 'n' ]]; then \
echo "HAVE_QR=0" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_QR=1" >> ${BUILDDIR}/Makefile.inc ; fi
endif
ifndef HAVE_THIRDPARTY
@echo "Checking for ThirdParty directory..."
@if [[ -d src/ThirdParty ]] ; then \
echo " Found." ; \
echo "HAVE_THIRDPARTY=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo " Not found. " ; \
echo "HAVE_THIRDPARTY=0" >> ${BUILDDIR}/Makefile.inc ; \
fi
endif
ifndef HAVE_FX
@echo
@echo "Would you like to compile with TTM3-F support? [y/n]"
@echo " Default: n"
@read HAVE_FX_Q && \
if [[ $$HAVE_FX_Q == 'y' ]]; then \
if [[ ! -d src/ThirdParty/nttm3f ]] ; then \
echo "TTM3-F sources must be present in the src/ThirdParty directory"; \
else \
echo "HAVE_FX=1" >> ${BUILDDIR}/Makefile.inc ; \
fi \
else \
echo "HAVE_FX=0" >> ${BUILDDIR}/Makefile.inc ; \
fi
endif
ifndef HAVE_SCME
@echo "Would you like to compile with SCME support ? [y/n]"
@echo " Default: n"
@read HAVE_SCME_Q && \
if [[ $$HAVE_SCME_Q == 'y' ]]; then \
if [[ ! -d src/ThirdParty/scme ]] ; then \
echo "SCME sources must be present in the src/ThirdParty directory"; \
else \
echo "HAVE_SCME=1" >> ${BUILDDIR}/Makefile.inc ; \
fi \
else \
echo "HAVE_SCME=0" >> ${BUILDDIR}/Makefile.inc ; \
fi
endif
ifndef HAVE_MTP
@echo "Would you like to compile with MTP support ? [y/n]"
@echo " Default: n"
@read HAVE_MTP && \
if [[ $$HAVE_MTP == 'y' ]]; then \
if [[ ! -d src/ThirdParty/mtp ]] ; then \
echo "MTP sources must be present in the src/ThirdParty directory"; \
else \
echo "HAVE_MTP=1" >> ${BUILDDIR}/Makefile.inc ; \
fi \
else \
echo "HAVE_MTP=0" >> ${BUILDDIR}/Makefile.inc ; \
fi
endif
ifndef HAVE_MBD
@echo "Would you like to compile with many-body dispersion support ? [y/n]"
@echo " Default: n"
@read HAVE_MBD && \
if [[ $$HAVE_MBD == 'y' ]]; then \
if [[ ! -d src/ThirdParty/MBD ]] ; then \
echo "MBD sources must be present in the src/ThirdParty directory"; \
else \
echo "HAVE_MBD=1" >> ${BUILDDIR}/Makefile.inc ; \
fi \
else \
echo "HAVE_MBD=0" >> ${BUILDDIR}/Makefile.inc ; \
fi
endif
ifndef HAVE_TTM_NF
@echo "Would you like to compile with TTM_NF support ? [y/n]"
@echo " Default: n"
@read HAVE_TTM_NF_Q && \
if [[ $$HAVE_TTM_NF_Q == 'y' ]]; then \
if [[ ! -d src/ThirdParty/ttm_nf ]] ; then \
echo "TTM_NF sources must be present in the src/ThirdParty directory"; \
else \
echo "HAVE_TTM_NF=1" >> ${BUILDDIR}/Makefile.inc ; \
fi \
else \
echo "HAVE_TTM_NF=0" >> ${BUILDDIR}/Makefile.inc ; \
fi
endif
ifndef HAVE_CH4
@echo "Would you like to compile with CH4 support ? [y/n]"
@echo " Default: n"
@read HAVE_CH4_Q && \
if [[ $$HAVE_CH4_Q == 'y' ]]; then \
if [[ ! -d src/ThirdParty/scme ]] ; then \
echo "CH4 sources must be present in the src/ThirdParty directory"; \
else \
echo "HAVE_CH4=1" >> ${BUILDDIR}/Makefile.inc ; \
fi \
else \
echo "HAVE_CH4=0" >> ${BUILDDIR}/Makefile.inc ; \
fi
endif
# OPTIONAL EXTERNAL LIBRARIES
#############################
@echo
@echo "OPTIONAL EXTERNAL LIBRARIES"
@echo "==========================="
@echo
ifndef HAVE_NETCDF4
@echo "Would you like to compile with NetCDF4 support? [y/n]"
@echo " Default: n"
@read -r HAVE_NETCDF4_Q && if [[ $$HAVE_NETCDF4_Q == 'y' ]]; then \
echo "Enter link flags for NetCDF4"; \
if NETCDF4_LIBS=$$(pkg-config --libs netcdf 2>/dev/null); then \
echo " Default from pkg-config: $$NETCDF4_LIBS"; \
else \
if NETCDF4_LIBS=$$(nc-config --libs 2>/dev/null); then \
echo " Default from nc-config: $$NETCDF4_LIBS"; \
else \
echo " No default found!"; \
NETCDF4_LIBS=''; \
fi \
fi ;\
read -r NETCDF4_LIBS_Q && if [[ -z $$NETCDF4_LIBS_Q ]]; then \
echo "NETCDF4_LIBS=$$NETCDF4_LIBS" >> ${BUILDDIR}/Makefile.inc; \
else \
echo "NETCDF4_LIBS=$$NETCDF4_LIBS_Q" >> ${BUILDDIR}/Makefile.inc; \
fi ;\
echo "Enter additonal flags for NetCDF (e.g. include directories: -I/usr/lib)"; \
if NETCDF4_FLAGS=$$(pkg-config --cflags netcdf 2>/dev/null); then \
echo " Default from pkg-config: $$NETCDF5_FLAGS"; \
else \
if NETCDF4_FLAGS=$$(nc-config --cflags 2>/dev/null); then \
echo " Default from nc-config: $$NETCDF4_FLAGS"; \
else \
echo " No default found!"; \
NETCDF4_FLAGS=''; \
fi \
fi ;\
read -r NETCDF4_FLAGS_Q && if [[ -z $$NETCDF4_FLAGS_Q ]]; then \
echo "NETCDF4_FLAGS=$$NETCDF4_FLAGS" >> ${BUILDDIR}/Makefile.inc; \
else \
echo "NETCDF4_FLAGS=$$NETCDF4_FLAGS_Q" >> ${BUILDDIR}/Makefile.inc; \
fi ;\
echo "HAVE_NETCDF4=1" >> ${BUILDDIR}/Makefile.inc; \
else \
echo "HAVE_NETCDF4=0" >> ${BUILDDIR}/Makefile.inc; \
fi
endif
ifndef HAVE_MDCORE
@echo
@echo "Please enter directories where MDCORE libraries are kept:"
@echo " Default: no MDCORE present"
@read MDCORE_LIBDIR && if [[ $$MDCORE_LIBDIR ]] ; then \
echo "MDCORE_LIBDIR=$$MDCORE_LIBDIR" >> ${BUILDDIR}/Makefile.inc ; \
echo ; \
echo "Please enter directory where MDCORE include files are kept:" ; \
read MDCORE_INCDIR && echo "MDCORE_INCDIR=$$MDCORE_INCDIR" >> ${BUILDDIR}/Makefile.inc ; \
echo "HAVE_MDCORE=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_MDCORE=0" >> ${BUILDDIR}/Makefile.inc;\
fi
endif
ifndef HAVE_ASAP
@echo
@echo "Please enter directory where the ASAP library is kept:"
@echo " Default: no ASAP library present"
@read ASAP_LIBDIR && if [[ $$ASAP_LIBDIR ]] ; then \
echo "ASAP_LIBDIR=$$ASAP_LIBDIR" >> ${BUILDDIR}/Makefile.inc ; \
echo "HAVE_ASAP=1" >> ${BUILDDIR}/Makefile.inc ; \
else echo "HAVE_ASAP=0" >> ${BUILDDIR}/Makefile.inc ; fi
endif
#ifndef HAVE_KIM
# @echo
# @echo "Please enter directory where the OpenKIM library is kept:"
# @echo " Default: no OpenKIM support"
# @read KIM_DIR && if [[ $$KIM_DIR ]]; then \
# echo "KIM_DIR=$$KIM_DIR" >> ${BUILDDIR}/Makefile.inc ; \
# echo "Turn off KIM test.kim autogeneration capability (useful for OpenKIM pipeline only)? [y/n]"; \
# echo " Default: n"; \
# read KIM_NO_AUTOGEN_Q && if [[ $$KIM_NO_AUTOGEN_Q == 'y' ]]; then \
# echo "KIM_NO_AUTOGENERATE_TEST_KIM=1" >> ${BUILDDIR}/Makefile.inc ; \
# else \
# echo "KIM_NO_AUTOGENERATE_TEST_KIM=0" >> ${BUILDDIR}/Makefile.inc ; fi; \
# echo "HAVE_KIM=1" >> ${BUILDDIR}/Makefile.inc ; \
# else \
# echo "HAVE_KIM=0" >> ${BUILDDIR}/Makefile.inc ; fi;
#endif
ifndef HAVE_CGAL
@echo
@echo "Please enter directories where CGAL libraries are kept:"
@echo " Default: no CGAL present"
@read CGAL_LIBDIR && if [[ $$CGAL_LIBDIR ]] ; then \
echo "CGAL_LIBDIR=$$CGAL_LIBDIR" >> ${BUILDDIR}/Makefile.inc ; \
echo ; \
echo "Please enter directory where CGAL include files are kept:" ; \
read CGAL_INCDIR && echo "CGAL_INCDIR=$$CGAL_INCDIR" >> ${BUILDDIR}/Makefile.inc ; \
echo "HAVE_CGAL=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_CGAL=0" >> ${BUILDDIR}/Makefile.inc;\
fi
endif
ifndef HAVE_METIS
@echo
@echo "Please enter directories where METIS libraries are kept:"
@echo " Default: no METIS present"
@read METIS_LIBDIR && if [[ $$METIS_LIBDIR ]] ; then \
echo "METIS_LIBDIR=$$METIS_LIBDIR" >> ${BUILDDIR}/Makefile.inc ; \
echo ; \
echo "HAVE_METIS=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_METIS=0" >> ${BUILDDIR}/Makefile.inc;\
fi
endif
ifndef HAVE_LMTO_TBE
@echo
@echo "Please enter directories where LMTO_TBE libraries are kept:"
@echo " Default: no LMTO_TBE present"
@read LMTO_TBE_LIBDIR && if [[ $$LMTO_TBE_LIBDIR ]] ; then \
echo "LMTO_TBE_LIBDIR=$$LMTO_TBE_LIBDIR" >> ${BUILDDIR}/Makefile.inc ; \
echo ; \
echo "HAVE_LMTO_TBE=1" >> ${BUILDDIR}/Makefile.inc ; \
else \
echo "HAVE_LMTO_TBE=0" >> ${BUILDDIR}/Makefile.inc;\
fi
endif
ifndef SIZEOF_FORTRAN_T
@/bin/bash ${QUIP_ROOT}/bin/find_sizeof_fortran_t >> ${BUILDDIR}/Makefile.inc
endif