# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 2
#     language: python
#     name: python2
# ---

import quippy

from quippy import descriptors

at = quippy.diamond(3.5, 6)
at.positions

at.Z

# 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"

# !quip descriptor_str="--help"

# +
# Pair distance descriptor between carbon atoms, with a 1.5 AA cutoff
# -

desc = descriptors.Descriptor("distance_2b Z1=6 Z2=6 cutoff=4")

desc.n_dim # number of dimensions (1, since it's a scalar distance)

# Get the connectivity
at.set_cutoff(desc.cutoff()+1.0)
at.calc_connect()

# How many instances of the descriptor can we use in the Atoms object?
desc.count(at)

d = desc.calc(at)
d # the actual pairwise distances (One array is the cutoff, the other the distances)

# %matplotlib inline
import matplotlib.pyplot as plt
plt.hist(d.descriptor)
plt.show()

n_desc, n_cross = desc.descriptor_sizes(at)
print("n_desc=%d n_cross=%d" % (n_desc,n_cross))

res = desc.calc(at, grad=True)
list(res)

res.descriptor[:, :10]

res.grad[:,:10]

res.index[:,:10]

# +
# Many-body descriptor: SOAP
# -

desc = descriptors.Descriptor("soap cutoff=3 l_max=4 n_max=4 atom_sigma=0.5 n_Z=1 Z={6} ")

at.set_cutoff(desc.cutoff())
at.calc_connect()

desc.descriptor_sizes(at) # Note one descriptor for each atom in the structure

desc.n_dim

desc.calc(at)

# Add a hydrogen atom, re-calculate the descriptor
at.add_atoms((0.2, 0.2, 0.2), (1))

at.calc_connect()
desc.descriptor_sizes(at) 
# Note no change in the size, as the SOAP descriptor above is just about carbons

# New H and C SOAP
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}")

desc.descriptor_sizes(at) # Note the size has gone up

desc.calc(at)

# #### We can get the same descriptor using the QUIP program, in which the input is an extended xyz file and the specification as before

at.write("diamondH.xyz")

# !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}"

# +
# ^ Each vector is on its own line, starting with a DESC
# -

# Three-site benzene monomer
desc_monomer = descriptors.Descriptor('general_monomer signature={6 6 6} atom_ordercheck=F')

desc_monomer.n_dim # Three distances -> dimensionality 3

benzat = quippy.AtomsList('benzene_frames.xyz') # some benzenes

# +
from ase.visualize import view
from ase.atoms import Atoms as AseAtoms 

firstBenzene = benzat[0]
view(AseAtoms(firstBenzene), viewer='x3d')
# -

 # Doesn't look much like benzene to me, but what do I know?

benzat.pos.shape

for i in range(len(benzat)):
    benzat[i].set_cutoff(desc_monomer.cutoff()+0.5)
    benzat[i].calc_connect()

res = desc_monomer.calc(benzat, grad=True)[0]

res.grad.shape



StarcoderdataJupyterScriptsDedupFiltered