Merge pull request #21 from lcmd-epfl/local-spahm

Local spahm

Author: briling

.gitignore

@@ -1,4 +1,129 @@
-__pycache__
-*.pyc
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
 *.so
-*.swp
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+.python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/

qstack/compound.py

@@ -282,3 +282,13 @@ def fragment_partitioning(fragments, prop_atom_inp, normalize=True):
         return props_frag
     else:
         return props_frag[0]
+
+
+def make_atom(q, basis):
+    mol = gto.Mole()
+    mol.atom = q + " 0.0 0.0 0.0"
+    mol.charge = 0
+    mol.spin = data.elements.ELEMENTS_PROTON[q] % 2
+    mol.basis = basis
+    mol.build()
+    return mol

qstack/math/matrix.py

@@ -15,3 +15,9 @@ def from_tril(mat_tril):
     mat[ind] = mat_tril
     mat = mat + mat.T - np.diag(np.diag(mat))
     return mat
+
+def sqrtm(m, eps=1e-13):
+    e, b = np.linalg.eigh(m)
+    e[abs(e) < eps] = 0.0
+    sm = b @ np.diag(np.sqrt(e)) @ b.T
+    return (sm+sm.T)*0.5

qstack/math/wigner.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+
+import sys
+import sympy as sp
+from sympy import symbols, Symbol, simplify, expand, cancel, expand_trig, Ynm, Ynm_c, Matrix, poly, zeros
+from .xyz_integrals_sym import xyz as xyzint
+
+
+# variables
+x,y,z    = symbols('x y z')
+x1,y1,z1 = symbols('x1 y1 z1')
+# coefficients
+xx,xy,xz = symbols('xx xy xz')
+yx,yy,yz = symbols('yx yy yz')
+zx,zy,zz = symbols('zx zy zz')
+
+def real_Y_correct_phase(l, m, theta, phi):
+  # returns real spherical harmonic in Condon--Shortley phase convention
+  # (sympy's Znm uses some other convention)
+  ym1 = Ynm  (l, -abs(m), theta, phi)
+  ym2 = Ynm_c(l, -abs(m), theta, phi)
+  if m==0:
+    return ym1
+  elif m<0:
+    return sp.I / sp.sqrt(2) * (ym1 - ym2)
+  elif m>0:
+    return 1 / sp.sqrt(2) * (ym1 + ym2)
+
+def get_polynom_Y(l, m):
+  # rewrites a real spherical harmonic as a polynom of x,y,z
+  theta = Symbol("theta", real=True)
+  phi = Symbol("phi", real=True)
+  r = Symbol('r', nonnegative=True)
+  expr = real_Y_correct_phase(l,m, theta, phi) * r**l
+  expr = expand(expr, func=True)
+  expr = expr.rewrite(sp.cos)#.simplify().trigsimp()
+  expr = expand_trig(expr)
+  expr = cancel(expr)
+  expr = expr.subs({r: sp.sqrt(x*x+y*y+z*z), phi: sp.atan2(y,x), theta: sp.atan2(sp.sqrt(x*x+y*y),z)})
+  if m!=0:
+    expr = cancel(expr).simplify()
+  expr = expr.subs({x*x+y*y: 1-z*z,
+                   3*x*x+3*y*y : 3-3*z*z })
+  return expr
+
+def xyzint_wrapper(knm, integrals_xyz_dict):
+  k,n,m = knm
+  if k%2 or n%2 or m%2:
+    return 0
+  else:
+    knm = tuple(sorted([k//2, n//2, m//2], reverse=True))
+    if not knm in integrals_xyz_dict.keys():
+      integrals_xyz_dict[knm] = xyzint(*knm)
+    return integrals_xyz_dict[knm]
+
+def product_Y(Y1,Y2):
+  # computes the product of two spherical harmonics
+  # and returns coefficients and a list of powers
+  prod = Y1 * Y2
+  prod = prod.expand().cancel()
+  prod = poly(prod, gens=[x,y,z])
+  return Matrix(prod.coeffs()), prod.monoms()
+
+
+def print_wigner(D):
+  for l,d in enumerate(D):
+    for m1 in range(-l,l+1):
+      for m2 in range(-l,l+1):
+        print('D[%d][% d,% d] = '%(l,m1,m2), d[m1,m2])
+    print()
+
+def compute_wigner(lmax):
+
+  Y     = [ [0]*(2*l+1) for l in range(lmax+1)]
+  Y_rot = [ [0]*(2*l+1) for l in range(lmax+1)]
+  for l in range(lmax+1):
+    for m in range(-l,l+1):
+      # spherical harmonic
+      Y[l][m] = get_polynom_Y(l, m)
+      # rotated spherical harmonic
+      Y_rot[l][m] = Y[l][m].subs({x: x1, y:y1, z:z1}).subs({x1:xx*x+xy*y+xz*z, y1:yx*x+yy*y+yz*z, z1:zx*x+zy*y+zz*z})
+
+
+  D = [zeros(2*l+1,2*l+1) for l in range(lmax+1)]
+  integrals_xyz_dict = {}
+  for l in range(lmax+1):
+    for m1 in range(-l,l+1):
+      for m2 in range(-l,l+1):
+        coefs, pows = product_Y(Y[l][m2], Y_rot[l][m1])
+        mono_integrals = [xyzint_wrapper(p,integrals_xyz_dict) for p in pows]
+        D[l][m1,m2] = coefs.dot(mono_integrals).factor() .subs({zx**2+zy**2: 1-zz**2, xx**2+xy**2:1-xz**2, yx**2+yy**2:1-yz**2}).simplify()
+  return D
+
+
+if __name__ == "__main__":
+  if len(sys.argv)<2:
+    lmax = 2
+  else:
+    lmax = int(sys.argv[1])
+
+  D = compute_wigner(lmax)
+  print_wigner(D)
+

qstack/math/xyz_integrals_float.py

@@ -0,0 +1,37 @@
+#!/usr/bin/env python3
+
+import sys
+
+def xyz(n,m,k):
+# computes the integral of x^2k y^2n z^2m over a sphere
+  k,n,m = sorted([k,n,m], reverse=True)
+  # k>=n>=m
+  if n==0:
+    xyz = 2.0 * (1.0 - (2.0*k-1.0)/(2.0*k+1.0))
+  else:
+    xyz = (2*k-1) * I23(n,m,k)
+  return xyz
+
+def I23(n,m,k):
+  I23 = 0.0
+  for l in range(0, n+m+2):
+    I23 = I23 + (-1)**l * trinomial( n+m+1, n+m+1-l, l) / (2.0*l+2.0*k-1.0)
+  I23 = I23 / ( (2*n+1) * 2**(2*n+2*m) )
+  for l in range(1, n+2):
+    I23 = I23 * (2*n+3-2*l) / (2*m-1+2*l)
+  return I23
+
+def trinomial(k1,k2,k3):
+# (k1+k2+k3)! / (k1! * k2! * k3!)
+  k1,k2,k3 = sorted([k1,k2,k3], reverse=True)
+  trinom = 1.0
+  for k in range(1,k2+1):
+    trinom = trinom * (k+k1) / k
+  for k in range(1,k3+1):
+    trinom = trinom * (k+k1+k2) / k
+  return trinom
+
+if __name__ == "__main__":
+  k,n,m = map(int, sys.argv[1:4])
+  print("%.15f π"%xyz(k,n,m))
+

qstack/math/xyz_integrals_sym.py

@@ -0,0 +1,38 @@
+#!/usr/bin/env python3
+
+import sys
+import sympy
+
+def xyz(n,m,k):
+# computes the integral of x^2k y^2n z^2m over a sphere
+  k,n,m = sorted([k,n,m], reverse=True)
+  # k>=n>=m
+  if n==0:
+    K = sympy.symbols('K')
+    xyz = (2 * (1 - (2*K-1)/(2*K+1))).subs(K,k)
+  else:
+    xyz = (2*k-1) * I23(n,m,k)
+  return xyz * sympy.pi
+
+def I23(n,m,k):
+  I23 = 0.0
+  K = sympy.symbols('K')
+  for l in range(0, n+m+2):
+    I23 = I23 + (-1)**l * trinomial( n+m+1, n+m+1-l, l) / (2*l+2*K-1)
+  I23 = I23.subs(K,k)
+  I23 = I23 / ( (2*n+1) * 2**(2*n+2*m) )
+  for l in range(1, n+2):
+    I23 = I23 * (2*n+3-2*l) / (2*m-1+2*l)
+  return I23
+
+def trinomial(k1,k2,k3):
+# (k1+k2+k3)! / (k1! * k2! * k3!)
+  k1,k2,k3 = sorted([k1,k2,k3])
+  trinom = sympy.FallingFactorial(k1+k2+k3, k3) / (sympy.factorial(k1)*sympy.factorial(k2))
+  return trinom
+
+if __name__ == "__main__":
+  k,n,m = map(int, sys.argv[1:4])
+  x = xyz(k,n,m)
+  print("%.15f ="%x, x)
+

qstack/spahm/rho/01_bind.py

@@ -0,0 +1,30 @@
+#!/usr/bin/env python3
+
+import os
+import argparse
+import numpy as np
+
+parser = argparse.ArgumentParser(description='This script generates a single file per atomic species binding all the feature-vector representations obtained using the single-structures representations generator script (1_DMbRep.py)')
+parser.add_argument('--X-dir',   required = True,                     type = str,            dest='XDir', help = 'The path to the directory containing all the single-structure X representations.')
+parser.add_argument('--species', default = ["C", "H", "N", "O", "S"], type = str, nargs='+', dest='Species', help = 'The species contained in the DatBase.')
+args = parser.parse_args()
+
+species = list(sorted(args.Species))
+X_files = sorted(os.listdir(args.XDir))
+
+X_atom = {}
+for e in species:
+    X_atom[e] = []
+
+for f in X_files :
+    X = np.load(args.XDir+'/'+f, allow_pickle=True)
+    for q,x in X :
+        X_atom[q].append(x)
+    print(f)
+
+for e in X_atom.keys() :
+    X_array = np.array(X_atom[e])
+    name_out = 'X_' + e + '_' + args.XDir
+    if name_out[-1] == '/': name_out = name_out[:-1]
+    np.save(name_out, X_array)
+    print('# Atoms for ', e, '=', X_array.shape)